g4data.py

Go to the documentation of this file.

#!/usr/bin/env python

from data import units
import elements

class G4RadData(object):
    
    # Cache all loaded states, map from (z,a) tuple
    states = {}

    def nuclideData(self,z,a):
        '''Return given nuclide data'''
        try:
            return G4RadData.states[(z,a)]
        except KeyError:
            pass
        try:
            nuc = NuclideData(z,a)
        except IOError,err:
            err = 'Illegal Z=%d, A=%d\n%s'%(z,a,err)
            raise ValueError,err
        G4RadData.states[(z,a)] = nuc
        return nuc

class NuclideDecay:
    'Holds decay type, fraction qvalue and daughter excitation energy'
    def __init__(self,type,fraction,qvalue,daughterExcitation=0,daughterZ=0,daughterA=0,rad = None):
        self.type = type
        self.fraction = fraction
        self.qvalue = qvalue
        self.daughterExcitation = daughterExcitation
        self.radiation = rad
        self.daughterZ = daughterZ
        self.daughterA = daughterA
        return

    def __str__(self):
        return '\t\t%s %6.3f%% Q=%8.4f [%8.4f] (%s)'%(self.type,self.fraction*100,self.qvalue,self.daughterExcitation,self.radiation)
    pass

class NuclideState:
    'Holds nuclide halflife, excitation energy and collection of NuclideDecay'
    def __init__(self,name,halflife=None):
        self.name = name
        self.halflife = halflife
        self.excitation = 0.0
        self.decays = []
        return

    def __str__(self):
        ret = ['\tstate %s (%s) [%.2f] %d decays:'%(self.name,self.halflife,self.excitation,len(self.decays))]
        for dk in self.decays:
            ret.append(str(dk))
            continue
        return '\n'.join(ret)
    pass


class NuclideData(object):
    'Hold all decay data about a nuclide from G4RadioActivity data files'
    
    def __init__(self,z,a):
        '''Create data for given Z and A'''
        self.Z = z
        self.A = a
        self.states = []
        self.name = '%d%s'%(self.A,elements.symbol(self.Z))
        self.load()
        return

    def __str__(self):
        ret = ['z%d.a%d %s %d excited state(s):'%(self.Z,self.A,self.name,len(self.states))]
        for s in self.states:
            ret.append(str(s))
            continue
        return '\n'.join(ret)

    def excitedState(self,nucexenergy = 0.0):
        '''Return state with given excitation energy'''
        for s in self.states:
            if nucexenergy == s.excitation: return s
            continue
        return None

    def load(self):
        import os
        dir = os.getenv('G4RADIOACTIVEDATA')
        if dir is None:
            raise EnvironmentError,'No G4RADIOACTIVEDATA variable, can not find data for Z=%d, A=%d'%(self.Z,self.A)
        file = dir + '/z%d.a%d'%(self.Z,self.A)
        file = open(file)
        types_seen = {}
        for line in file.readlines():
            line = line.strip()
            #print line

            # skip comments
            if line[0] == '#': continue

            # parent
            if line[0] == 'P':
                self.states.append(NuclideState(self.name))
                types_seen = {}
                p,ex,hl = line.split()
                self.states[-1].excitation = float(ex)
                self.states[-1].halflife = float(hl)*units.second
                continue

            # warning
            if line[0] == 'W':
                #print line
                continue
                
            # now parsing data
            chunks = line.split()

            # skip first total BR for type
            if not types_seen.has_key(chunks[0]):
                types_seen[chunks[0]] = float(chunks[2])
                continue

            dk = NuclideDecay(chunks[0],float(chunks[2])/100.0,
                              float(chunks[3])*units.keV,float(chunks[1])*units.keV)
            import decay
            if dk.type == 'Alpha':
                dk.daughterZ = self.Z - 2
                dk.daughterA = self.A - 4
                dk.radiation = decay.AlphaDecay(dk.qvalue,self.Z)
            elif dk.type == 'BetaMinus':
                dk.daughterZ = self.Z + 1
                dk.daughterA = self.A
                dk.radiation = decay.BetaDecay(dk.qvalue,self.Z)
            elif dk.type == 'BetaPlus':
                dk.daughterZ = self.Z - 1
                dk.daughterA = self.A
                dk.radiation = decay.BetaDecay(dk.qvalue,-1*self.Z)
            else:
                print 'Warning - not yet supported decay type: %s'%dk.type
                continue

            self.states[-1].decays.append(dk)
            continue
        return

    pass

if '__main__' == __name__:
    import sys
    nd = NuclideData(int(sys.argv[1]),int(sys.argv[2]))
    print nd