Skip to content

Commit 4bbdc69

Browse files
lukeshingleslukeshingles
committed
Update makenucleardecayfile.py
1 parent ec1ba60 commit 4bbdc69

File tree

1 file changed

+167
-130
lines changed

1 file changed

+167
-130
lines changed

‎artisatomic/makenucleardecayfile.py‎

Lines changed: 167 additions & 130 deletions
Original file line numberDiff line numberDiff line change
@@ -5,164 +5,201 @@
55

66
import artistools as at
77
import numpy as np
8-
import pandas as pd
98
import polars as pl
109
import requests
1110

11+
colreplacements = {
12+
"Rad Int.": "intensity",
13+
"Rad Ene.": "radiationenergy_kev",
14+
"Rad subtype": "radsubtype",
15+
"Par. Elevel": "parent_elevel",
16+
"Dec Mode": "decaymode",
17+
"T1/2 (num)": "halflife_s",
18+
}
19+
20+
21+
def process_table(z: int, a: int, outpath: Path, textdata, dfbetaminus: pl.DataFrame, dfalpha: pl.DataFrame):
22+
textdata = textdata.replace("**********", "0.")
23+
# print(textdata)
24+
# match
25+
26+
startindex = textdata.find("<pre>") + len("<pre>")
27+
endindex = textdata.rfind("</pre>")
28+
strtable = textdata[startindex:endindex].strip()
29+
strheader = strtable.strip().split("\n")[0].strip()
30+
# print(strheader)
31+
assert (
32+
strheader
33+
== "A Element Z N Par. Elevel Unc. JPi Dec Mode T1/2 (txt) T1/2 (num) "
34+
" Daughter Radiation Rad subtype Rad Ene. Unc EP Ene. Unc Rad Int. Unc "
35+
" Dose Unc"
36+
)
37+
38+
dfnuclide = pl.read_csv(io.StringIO(strtable), separator="\t", schema_overrides={"Par. Elevel": pl.Utf8})
39+
40+
newcols: list[str] = []
41+
for colname in dfnuclide.columns:
42+
colname = colname.strip()
43+
if colname.startswith("Unc"):
44+
colname = newcols[-1] + " (Unc)"
45+
if colname in colreplacements:
46+
colname = colreplacements[colname]
47+
newcols.append(colname)
48+
dfnuclide.columns = newcols
49+
dfnuclide = dfnuclide.with_columns(pl.col(pl.Utf8).str.strip_chars()).with_columns(
50+
pl.col("parent_elevel").str.replace("0.0", "0"),
51+
pl.col("radiationenergy_kev").cast(pl.Float64),
52+
pl.col("intensity").cast(pl.Float64),
53+
pl.col("halflife_s").cast(pl.Float64),
54+
)
55+
56+
found_groundlevel = False
57+
for (parelevel,), dfdecay in dfnuclide.group_by("parent_elevel"):
58+
assert isinstance(parelevel, str)
59+
try:
60+
is_groundlevel = float(parelevel) == 0.0
61+
except ValueError:
62+
is_groundlevel = False
63+
print(f" parent_Elevel: {parelevel} is_groundlevel: {is_groundlevel}")
64+
if not is_groundlevel:
65+
continue
66+
67+
found_groundlevel = True
68+
69+
dfgammadecays = dfdecay.filter(
70+
(pl.col("Radiation") == "G") & pl.col("radsubtype").is_in(["", "Annihil."]) & (pl.col("intensity") > 0.0)
71+
)
72+
73+
maybedfbetaminusrow = dfbetaminus.filter(pl.col("Z") == z).filter(pl.col("A") == a)
74+
maybedfalpharow = dfalpha.filter(pl.col("Z") == z).filter(pl.col("A") == a).limit(1)
75+
nndc_halflife = None
76+
if not dfgammadecays.is_empty():
77+
nndc_halflife = dfgammadecays["halflife_s"].item(0)
78+
print(f" NNDC half-life: {nndc_halflife:7.1e} s")
79+
80+
if maybedfbetaminusrow.height > 0:
81+
halflife = maybedfbetaminusrow["tau[s]"].item() * math.log(2)
82+
strwarn = (
83+
" WARNING!!!!!!"
84+
if (nndc_halflife is not None and not np.isclose(nndc_halflife, halflife, rtol=0.1))
85+
else ""
86+
)
87+
print(f" betaminusdecays.txt half-life: {halflife:7.1e} s {strwarn}")
88+
89+
if maybedfalpharow.height > 0:
90+
halflife = maybedfalpharow["halflife_s"].item()
91+
strwarn = (
92+
" WARNING!!!!!!"
93+
if (nndc_halflife is not None and not np.isclose(nndc_halflife, halflife, rtol=0.1))
94+
else ""
95+
)
96+
print(f" alphadecays.txt half-life: {halflife:7.1e} s {strwarn}")
97+
98+
e_gamma = (dfgammadecays["radiationenergy_kev"] * dfgammadecays["intensity"] / 100.0).sum()
99+
print(f" NNDC Egamma: {e_gamma:7.1f} keV")
100+
101+
if maybedfbetaminusrow.height > 0:
102+
file_e_gamma = maybedfbetaminusrow["E_gamma[MeV]"].item() * 1000
103+
strwarn = "" if np.isclose(e_gamma, file_e_gamma, rtol=0.1) else " WARNING!!!!!!"
104+
print(f" betaminusdecays.txt Egamma: {file_e_gamma:7.1f} keV {strwarn}")
105+
106+
elif maybedfalpharow.height > 0:
107+
file_e_gamma = maybedfalpharow["E_gamma[MeV]"].item() * 1000
108+
strwarn = "" if np.isclose(e_gamma, file_e_gamma, rtol=0.1) else " WARNING!!!!!!"
109+
print(f" alphadecays.txt Egamma: {file_e_gamma:7.1f} keV {strwarn}")
110+
111+
dfout = pl.DataFrame(
112+
{
113+
"energy_mev": dfgammadecays["radiationenergy_kev"] / 1000.0,
114+
"intensity": dfgammadecays["intensity"] / 100.0,
115+
}
116+
).sort("energy_mev")
117+
if len(dfout) > 0:
118+
with outpath.open("w", encoding="utf-8") as fout:
119+
fout.write(f"{len(dfout)}\n")
120+
for energy_mev, intensity in dfout[["energy_mev", "intensity"]].iter_rows():
121+
fout.write(f"{energy_mev:5.3f} {intensity:6.4f}\n")
122+
123+
print(f"Saved {outpath.name}")
124+
else:
125+
print("empty DataFrame")
126+
if not found_groundlevel:
127+
print(" ERROR! did not find ground level")
128+
12129

13130
def main():
14131
PYDIR = Path(__file__).parent.resolve()
15-
atomicdata = pd.read_csv(PYDIR / "atomic_properties.txt", sep=r"\s+", comment="#")
16-
elsymbols = ["n", *list(atomicdata["symbol"].values)]
132+
atomicdata = pl.read_csv(PYDIR / "atomic_properties.txt", separator=" ", comment_prefix="#", null_values="NA")
133+
elsymbols = at.get_elsymbolslist()
17134

18135
outfolder = Path(__file__).parent.parent.absolute() / "artis_files" / "data"
19136
outfolder.mkdir(parents=True, exist_ok=True)
20-
pd.options.display.max_rows = 999
21-
pd.options.display.max_columns = 999
22-
23-
dfbetaminus = pl.read_csv(
24-
at.get_config()["path_datadir"] / "betaminusdecays.txt",
25-
separator=" ",
26-
comment_prefix="#",
27-
has_header=False,
28-
new_columns=["A", "Z", "Q[MeV]", "E_gamma[MeV]", "E_elec[MeV]", "E_neutrino[MeV]", "tau[s]"],
29-
).filter(pl.col("Q[MeV]") > 0.0)
30-
31-
dfalpha = pl.read_csv(
32-
at.get_config()["path_datadir"] / "alphadecays.txt",
33-
separator=" ",
34-
comment_prefix="#",
35-
has_header=False,
36-
new_columns=[
37-
"A",
38-
"Z",
39-
"branch_alpha",
40-
"branch_beta",
41-
"halflife_s",
42-
"Q_total_alphadec[MeV]",
43-
"Q_total_betadec[MeV]",
44-
"E_alpha[MeV]",
45-
"E_gamma[MeV]",
46-
"E_beta[MeV]",
47-
],
137+
138+
dfbetaminus = (
139+
pl.read_csv(
140+
at.get_config()["path_datadir"] / "betaminusdecays.txt",
141+
separator=" ",
142+
comment_prefix="#",
143+
has_header=False,
144+
new_columns=["A", "Z", "Q[MeV]", "E_gamma[MeV]", "E_elec[MeV]", "E_neutrino[MeV]", "tau[s]"],
145+
)
146+
.filter(pl.col("Q[MeV]") > 0.0)
147+
.filter(pl.col("tau[s]") > 0.0)
48148
)
49149

50-
nuclist = sorted(list(dfbetaminus.select(["Z", "A"]).iter_rows()) + list(dfalpha.select(["Z", "A"]).iter_rows()))
150+
assert dfbetaminus.height == dfbetaminus.unique(("Z", "A")).height
151+
152+
dfalpha = (
153+
pl.read_csv(
154+
at.get_config()["path_datadir"] / "alphadecays.txt",
155+
separator=" ",
156+
comment_prefix="#",
157+
has_header=False,
158+
new_columns=[
159+
"A",
160+
"Z",
161+
"branch_alpha",
162+
"branch_beta",
163+
"halflife_s",
164+
"Q_total_alphadec[MeV]",
165+
"Q_total_betadec[MeV]",
166+
"E_alpha[MeV]",
167+
"E_gamma[MeV]",
168+
"E_beta[MeV]",
169+
],
170+
)
171+
.filter(pl.col("halflife_s") > 0.0)
172+
.unique(("Z", "A"), keep="last", maintain_order=True)
173+
)
174+
assert dfalpha.height == dfalpha.unique(("Z", "A")).height
51175

52-
colreplacements = {
53-
"Rad Int.": "intensity",
54-
"Rad Ene.": "radiationenergy_kev",
55-
"Rad subtype": "radsubtype",
56-
"Par. Elevel": "parent_elevel",
57-
}
176+
nuclist = sorted(list(dfbetaminus.select(["Z", "A"]).iter_rows()) + list(dfalpha.select(["Z", "A"]).iter_rows()))
58177

59178
for z, a in nuclist:
60179
strnuclide = elsymbols[z].lower() + str(a)
61-
print(f"\n(Z={z}) {strnuclide}")
62180
filename = f"{strnuclide}_lines.txt"
63181
outpath = outfolder / filename
64182
# if outpath.is_file():
65-
# print(f" {filename} already exists. skipping...")
183+
# # print(f" {filename} already exists. skipping...")
66184
# continue
185+
print(f"\n(Z={z}) {strnuclide}")
67186

68187
url = f"https://www.nndc.bnl.gov/nudat3/decaysearchdirect.jsp?nuc={strnuclide}&unc=standard&out=file"
69188

70189
with requests.Session() as s:
71190
textdata = s.get(url).text
72-
textdata = textdata.replace("**********", "0.")
73-
# print(textdata)
74-
# match
75191
if "<pre>" not in textdata:
76192
print(f" no table data returned from {url}")
77-
continue
78-
79-
startindex = textdata.find("<pre>") + len("<pre>")
80-
endindex = textdata.rfind("</pre>")
81-
strtable = textdata[startindex:endindex].strip()
82-
strheader = strtable.strip().split("\n")[0].strip()
83-
# print(strheader)
84-
assert (
85-
strheader
86-
== "A Element Z N Par. Elevel Unc. JPi Dec Mode T1/2 (txt) T1/2 (num) "
87-
" Daughter Radiation Rad subtype Rad Ene. Unc EP Ene. Unc Rad Int. Unc "
88-
" Dose Unc"
89-
)
90-
# dfnuclide = pd.read_fwf(io.StringIO(strtable))
91-
dfnuclide = pd.read_csv(io.StringIO(strtable), delimiter="\t", dtype={"Par. Elevel": str})
92-
93-
newcols = []
94-
for colname in dfnuclide.columns:
95-
colname = colname.strip()
96-
if colname.startswith("Unc"):
97-
colname = newcols[-1] + " (Unc)"
98-
if colname in colreplacements:
99-
colname = colreplacements[colname]
100-
newcols.append(colname)
101-
dfnuclide.columns = newcols
102-
dfnuclide["Dec Mode"] = dfnuclide["Dec Mode"].str.strip()
103-
dfnuclide["Radiation"] = dfnuclide["Radiation"].str.strip()
104-
dfnuclide["radsubtype"] = dfnuclide["radsubtype"].str.strip()
105-
dfnuclide["parent_elevel"] = dfnuclide["parent_elevel"].str.strip()
106-
107-
found_groundlevel = False
108-
for parelevel, dfdecay in dfnuclide.groupby("parent_elevel"):
109-
try:
110-
is_groundlevel = float(parelevel) == 0.0
111-
except ValueError:
112-
is_groundlevel = False
113-
print(f" parent_Elevel: {parelevel} is_groundlevel: {is_groundlevel}")
114-
if not is_groundlevel:
115-
continue
116-
found_groundlevel = True
117-
dfgammadecays = dfdecay.query(
118-
"Radiation == 'G' and (radsubtype == '' or radsubtype == 'Annihil.') and intensity >= 0.15",
119-
inplace=False,
193+
else:
194+
process_table(
195+
z=z,
196+
a=a,
197+
outpath=outpath,
198+
textdata=textdata,
199+
dfbetaminus=dfbetaminus,
200+
dfalpha=dfalpha,
120201
)
121202

122-
maybedfbetaminusrow = dfbetaminus.filter(pl.col("Z") == z).filter(pl.col("A") == a)
123-
maybedfalpharow = dfalpha.filter(pl.col("Z") == z).filter(pl.col("A") == a)
124-
if not dfgammadecays.empty:
125-
print(f" NNDC half-life: {dfgammadecays.iloc[0]['T1/2 (num)']:7.1e} s")
126-
127-
if maybedfbetaminusrow.height > 0:
128-
halflife = maybedfbetaminusrow["tau[s]"].item() * math.log(2)
129-
print(f" betaminusdecays.txt half-life: {halflife:7.1e} s")
130-
if maybedfalpharow.height > 0:
131-
print(f" alphadecays.txt half-life: {maybedfalpharow['halflife_s'].item():7.1e} s")
132-
133-
e_gamma = (dfgammadecays["radiationenergy_kev"] * dfgammadecays["intensity"] / 100.0).sum()
134-
print(f" NNDC Egamma: {e_gamma:7.1f} keV")
135-
136-
if maybedfbetaminusrow.height > 0:
137-
file_e_gamma = maybedfbetaminusrow["E_gamma[MeV]"].item() * 1000
138-
print(f" betaminusdecays.txt Egamma: {file_e_gamma:7.1f} keV")
139-
if not np.isclose(e_gamma, file_e_gamma, rtol=0.1):
140-
print("WARNING!!!!!!")
141-
142-
elif maybedfalpharow.height > 0:
143-
file_e_gamma = maybedfalpharow["E_gamma[MeV]"].item() * 1000
144-
print(f" alphadecays.txt Egamma: {file_e_gamma:7.1f} keV")
145-
if not np.isclose(e_gamma, file_e_gamma, rtol=0.1):
146-
print("WARNING!!!!!!")
147-
148-
dfout = pl.DataFrame(
149-
{
150-
"energy_mev": dfgammadecays.radiationenergy_kev.to_numpy() / 1000.0,
151-
"intensity": dfgammadecays.intensity.to_numpy() / 100.0,
152-
}
153-
).sort("energy_mev")
154-
if len(dfout) > 0:
155-
with outpath.open("w", encoding="utf-8") as fout:
156-
fout.write(f"{len(dfout)}\n")
157-
for energy_mev, intensity in dfout[["energy_mev", "intensity"]].iter_rows():
158-
fout.write(f"{energy_mev:5.3f} {intensity:6.4f}\n")
159-
160-
print(f"Saved {filename}")
161-
else:
162-
print("empty DataFrame")
163-
if not found_groundlevel:
164-
print(" ERROR! did not find ground level")
165-
166203

167204
if __name__ == "__main__":
168205
main()

0 commit comments

Comments
 (0)