Skip to content

Commit

Permalink
Decay builder using new DB interface
Browse files Browse the repository at this point in the history
  • Loading branch information
Michael Mendenhall committed Apr 22, 2016
1 parent 4d6042d commit 507cc28
Show file tree
Hide file tree
Showing 3 changed files with 208 additions and 163 deletions.
288 changes: 138 additions & 150 deletions Physics/DecaySpecBuilder.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,10 +2,27 @@
"""@package DecaySpecBuilder
Builds decay generator specification files from ENSDF inputs"""

from ENSDF_Reader import *
from ENSDF_ParsedDB import *
from SMFile import *
from copy import *

class ElementNames:
"""Element name/symbol lookup table"""
def __init__(self, fname="ElementsList.txt"):
self.dat = [(int(x[0]),x[1],x[2]) for x in [l.split() for l in open(fname).readlines()] if len(x)==3]
self.byNum = dict([(e[0],e) for e in self.dat])
self.bySymb = dict([(e[1].lower(),e) for e in self.dat])
self.byName = dict([(e[2].lower(),e) for e in self.dat])

def elNum(self, s):
"""Element Z by name or symbol, case-insensitive"""
s = s.lower()
n = self.bySymb.get(s,None)
n = n if n else self.byName.get(s,None)
return n[0] if n else None

TheElementNames = ElementNames()

class CETable:
"""Conversion electron data table gleaned from a gamma line entry"""

Expand Down Expand Up @@ -38,199 +55,170 @@ def getCE(self):
c = ce[i]
break
if c:
rs = "%g"%c.tofloat()
if c.dv:
rs += "~%g"%c.dv
v,dv = parse_stderr(*c.split())
rs = "%g"%v
if dv:
rs += "~%g"%dv
ces[k] = rs
return ces

class DecaySpecBuilder:
"""Builds decay generator specification files from ENSDF inputs"""

def __init__(self):
def __init__(self,curs):
self.curs = curs
self.sheets = [] # loaded ENSDF datasheets
self.levelcounter = { } # level number assignment count
self.levidx = { } # index of registered levels by level ID
self.parents = [] # collection of all parent specifiers in loaded cards
self.joins = [] # joining points between multiple cards
self.min_tx_prob = 0 # ignore transitions less likely than this

def add_sheet(self, s, linklevel = None):
"""Add parsed ENSDF datasheet; assigns unique global level numbers"""
def add_sheet(self, cid):
"""Add parsed ENSDF datasheet"""
card = ParsedCard(self.curs,cid)
print("\nAdding card:")
card.display()

# possible link into existing level scheme
if linklevel:
assert linklevel.levelID in self.levidx
E0 = linklevel.E.tofloat()
E1 = self.levidx[linklevel.levelID].E.tofloat()
if E0 > E1:
for sh in self.sheets:
sh.shift_E0(E0-E1)
else:
s.shift_E0(E1-E0)
# auto-join multiple branches from same parent
card.shiftE = 0
for p in card.parents:
p.E = p.QP
isJoined = False
for p1 in self.parents:
if p.mass == p1.mass and p.elem == p1.elem and p.J == p1.J:
card.shiftE = p1.card.shiftE + p1.E - p.E
self.joins.append((p,p1))
isJoined = True
print("Joining parents %s"%p)
break
if isJoined: continue
self.parents.append(p)

self.sheets.append(card)

self.sheets.append(s)

for li in s.levidx:
l = s.levels[li]
levid = (l.nucA, l.nucZ)
l.nIn = 0
l.nOut = 0
if hasattr(l,"levelID"):
assert l.levelID[:2] == levid
print("Linking level",l.levelID)
continue
n = self.levelcounter.setdefault(levid, -1) + 1
l.levelID = (l.nucA, l.nucZ, n)
self.levidx[l.levelID] = l
self.levelcounter[levid] += 1
print(l.levelID)

def assign_level_names(self):
"""Assign output names to active levels"""
self.levels = []
for s in self.sheets:
self.levels += s.levels.values()
self.levels += self.parents
self.levels.sort(key=(lambda l: l.E))
isotcounts = {}
for l in self.levels:
a = l.mass
z = TheElementNames.elNum(l.elem)
i = isotcounts.setdefault((a,z), 0)
l.outName = "%i.%i.%i"%(a,z,i)
isotcounts[(a,z)] += 1
for j in self.joins: j[0].outName = j[1].outName

def headercomments(self):
"""Generate header comments block"""
h = "#########################################\n"
h += "# Decay scheme generated from ENSDF data:\n#\n"
for s in self.sheets:
rID = s.IDrecord
h += "# '%s' %s:"%(rID.DSID, rID.DATE)
rID = s.idrec
h += "# '%s' %s:"%(rID.DSID, rID.EDATE)
for rH in s.history:
if "CIT" in rH.xvals:
h += " %s"%rH.xvals["CIT"][0][1]
if "CIT" in rH.xdata:
h += " %s"%rH.xdata["CIT"][0][1]
h += "\n"
return h

def levellist(self, smf):
"""Output energy levels to SMFile"""
for l in self.levidx.values():
if not l.nOut and not l.nIn:
continue
self.assign_level_names()

for l in self.levels:

ml = KVMap()
ml.insert("nm", "%i.%i.%i"%l.levelID)
ml.insert("E", "%.2f"%l.E.tofloat())
hl = -1 if l.T.vs == "STABLE" else l.T.tofloat()
if l.T.unit in ["MeV","MEV"]: hl = 0
if hl is None: hl = 0
ml.insert("nm", l.outName)
ml.insert("E", "%.2f"%(l.E + l.card.shiftE))

try: hl = float(l.T)
except: hl = 0
if hl == float("inf"): hl = -1
ml.insert("hl","%g"%hl)

if l.J:
ml.insert("jpi",l.J)
if l.J: ml.insert("jpi",l.J)

smf.insert("level", ml)

def transitionList(self, smf):
"""Output transitions to SMFile"""

def maketx(tx):
try:
m = KVMap()
m.insert("from", tx.fromlvl.outName)
m.insert("to", tx.tolvl.outName)
return m
except: return None

for s in self.sheets:

assert len(s.parents) == 1 # almost always true...
P = tuple(s.parents.values())[0] # parent for file

for l in s.levels.values():
for a in s.alphas:
if a.Ialpha <= self.min_tx_prob: continue
ma = maketx(a)
if not ma: continue
ma.insert("I", "%g"%a.Ibeta)
smf.insert("alpha", ma)

for b in s.betas:
if b.Ibeta <= self.min_tx_prob: continue
mb = maketx(b)
if not mb: continue
mb.insert("I", "%g"%b.Ibeta)
if b.UN and b.UN[1:] == "U":
mb.insert("forbidden", b.UN[:1])
smf.insert("beta", mb)

for g in s.gammas:
if g.Igamma <= self.min_tx_prob: continue

# normalization record for daughter, if available
nucnrm = P.norms.get(l.NUCID, None)
nnBR = nucnrm.BR.tofloat() if nucnrm else None
if nnBR is None: nnBR = 1.

for g in l.gammas:
Ig = g.RI.tofloat()
if Ig and nucnrm: Ig *= nucnrm.NR.tofloat() * nnBR
if not Ig or Ig <= self.min_tx_prob: continue
mg = maketx(g)
if not mg: continue
mg.insert("Igamma", "%g"%g.Igamma)
mg.insert("E", "%g"%g.E)

mg = KVMap()
mg.insert("from", "%i.%i.%i"%l.levelID)
mg.insert("to", "%i.%i.%i"%g.goesto.levelID)
mg.insert("Igamma", "%g"%Ig)
mg.insert("E", "%g"%g.E.tofloat())

ces = CETable(g.xvals).getCE()
for c in ces: mg.insert("CE_"+c, ces[c])
ces = CETable(g.xdata).getCE()
for c in ces: mg.insert("CE_"+c, ces[c])

smf.insert("gamma", mg)
l.nOut += 1
g.goesto.nIn += 1

for f in l.feeders:

mf = KVMap()
mf.insert("to", "%i.%i.%i"%l.levelID)

# identify parent level
pnt = None
if f.RTYPE == "E": pnt = s.findParent(f.nucA, f.nucZ+1)
elif f.RTYPE == "B": pnt = s.findParent(f.nucA, f.nucZ-1)
elif f.RTYPE == "A": pnt = s.findParent(f.nucA+4, f.nucZ+2)
assert pnt
mf.insert("from", "%i.%i.%i"%pnt.asLevel.levelID)


# beta, electron capture branches normalization
NB = 1
if nucnrm:
nnNB = nucnrm.NB.tofloat()
if nnNB is None: nnNB = 1.
NB = nnNB * nnBR

if f.RTYPE == "E":
doesEC = False
if f.IE.tofloat(): # electron capture
tprob = NB*f.IE.tofloat()
if tprob > self.min_tx_prob:
mfe = deepcopy(mf)
mfe.insert("I","%g"%tprob)
smf.insert("ecapt", mfe)
doesEC = True
if f.IB.tofloat(): # beta+
tprob = NB*f.IB.tofloat()
if tprob > self.min_tx_prob:
mfp = deepcopy(mf)
mfp.insert("I","%g"%tprob)
smf.insert("beta", mfp)
doesEC = True
if not doesEC: continue
elif f.RTYPE == "B": # beta- decay
tprob = NB*f.IB.tofloat()
if tprob <= self.min_tx_prob: continue
mf.insert("I","%g"%tprob)
smf.insert("beta", mf)
elif f.RTYPE == "A": # alpha decay
tprob = f.IA.tofloat()*100
if tprob <= self.min_tx_prob: continue
mf.insert("I","%g"%tprob)
mf.insert("E","%g"%f.E.tofloat())
smf.insert("alpha", mf)

l.nIn += 1
pnt.asLevel.nOut += 1

smf.insert("gamma", mg)

for e in s.ecapts:
if self.min_tx_prob < e.Iec:
me = maketx(e)
if not me: continue
me.insert("I","%g"%e.Iec)
smf.insert("ecapt", me)
if self.min_tx_prob < e.Ibeta:
me = maketx(e)
if not me: continue
me.insert("I","%g"%e.Ibeta)
smf.insert("beta", me)

if __name__=="__main__":
basedir = "/home/mpmendenhall/Documents/PROSPECT/RefPapers/ENSDF/"

DSB = DecaySpecBuilder()
datdir = "/home/mpmendenhall/Documents/ENSDF/"
conn = sqlite3.connect(datdir+"ENSDF.db")
conn.row_factory = sqlite3.Row # fast name-based access to columns
curs = conn.cursor()

if True:
# single-sheet case
DSB.add_sheet( ENSDF_Reader(basedir + "ENSDF_22Na-22Ne.txt") )
DSB = DecaySpecBuilder(curs)

if False:
DSB.add_sheet( ENSDF_Reader(basedir + "ENSDF_214Bi-214Po.txt") )

s2 = ENSDF_Reader(basedir + "ENSDF_214Po-210Pb.txt")
pl = s2.findParent(214,84).asLevel
pl.levelID = (214,84,0)
DSB.add_sheet(s2, pl)

if False:
DSB.add_sheet( ENSDF_Reader(basedir + "ENSDF_152Eu-152Gd.txt") )
s2 = ENSDF_Reader(basedir + "ENSDF_152Eu-152Sm.txt")
pl = s2.findParent(152,63).asLevel
pl.levelID = (152,63,0)
DSB.add_sheet(s2, pl)
if False: # example with two joined branches
curs.execute("SELECT rowid FROM ENSDF_cards WHERE DSID LIKE '40K %'")
for cid in curs.fetchall(): DSB.add_sheet(cid[0])
elif True:
curs.execute("SELECT rowid FROM ENSDF_cards WHERE DSID LIKE '207BI EC%'")
for cid in curs.fetchall(): DSB.add_sheet(cid[0])

print(DSB.headercomments())
smf = SMFile()
DSB.levellist(smf)

#DSB.min_tx_prob = 1
DSB.transitionList(smf)
DSB.levellist(smf)


print(smf.toString())
1 change: 1 addition & 0 deletions Physics/ENSDF_CardDB.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,6 +87,7 @@ def to_text(self):

if True:
upload_dataset(curs,datdir+"ensdf_150208/")
#upload_cards(curs, open(datdir+"ensdf_150208/ensdf_151208.upd").read())
conn.commit()
else:
curs.execute('select rowid from ENSDF_cards WHERE mass = 12 AND elem = "C"')
Expand Down
Loading

0 comments on commit 507cc28

Please sign in to comment.