In This Package:

chain.py

Go to the documentation of this file.

#!/usr/bin/env python

'''
Assemble radioactive decay chains.
'''

from g4data import G4RadData

class G4DataChain(object):
    'A radioactive decay chain built from G4RadioActivity data files'

    # cache data.Isotope instantiations.  Look up via (z,a) tuple,
    # return list of that Isotope and different excitation energies
    _isotopes = {}

    # Precision by which to judge equality of nuclear excited state
    # energies
    energyPrecision = 0.01

    def __init__(self,z,a,correlation_time = None):
        ''' Build a decay chain starting from the given (ground state)
        nuclide with Z and A.  If a correlation time is given it will
        be used to separate the chain into sub-chains that start at
        nuclides with half lives greater than the correlation time.'''

        self.corrTime = correlation_time
        
        self.subChain = []
        self.groundStates = []
        self.head = self.chain(z,a)
        return

    def __str__(self):
        return str(self.head)

    def getIsotope(self,z,a,nucexenergy):
        'Retrieve isotope, return (iso,BOOL), BOOL is True if iso is freshly constructed'

        from data import Isotope

        try:
            isos = G4DataChain._isotopes[(z,a)]
        except KeyError:
            iso = Isotope(z,a,level=nucexenergy)
            G4DataChain._isotopes[(z,a)] = [iso]
            #print 'Making first seen: %s'%iso
            return iso,True

        for iso in isos:
            import math
            if math.fabs(nucexenergy - iso.level) <= G4DataChain.energyPrecision*nucexenergy:
                #print 'Found in cache: %s'%iso
                return iso,False # cache hit
            continue

        # found cached isotopes but wrong energy level
        iso = Isotope(z,a,level=nucexenergy)
        G4DataChain._isotopes[(z,a)].append(iso)
        assert len(G4DataChain._isotopes[(z,a)]) > 1
        #print 'In cache, new energy: %s'%iso
        return iso,True
        

    def chain(self,z,a,nucexenergy = 0.0):
        'Return Isotope coresponding to given inputs'
        import math, decay
        from data import units, Transition

        iso,new = self.getIsotope(z,a,nucexenergy)
        
        if new and nucexenergy == 0.0: self.groundStates.append(iso)
        
        if not new:                  # been there, done that
            #print "Already chained: %s"%iso
            return iso
        #print 'Chaining: %s'%iso

        # If nucleous is excited, deal with it promptly
        if nucexenergy > G4DataChain.energyPrecision*units.keV:
            daughter = self.chain(z,a)
            tran = Transition(iso,daughter,decay.GammaDecay(nucexenergy))
            return iso

        try:
            rd = G4RadData()
            nuc = rd.nuclideData(z,a)
        except ValueError,err:
            print 'No data for z=%d, a=%d, assuming stable (%s)'%(z,a,err)
            return iso

        state = nuc.excitedState(nucexenergy)
        if state is None: 
            print 'No excited state for z=%d, a=%d, E = %f, assuming stable'%(z,a,nucexenergy)
            return iso
        iso.halflife = state.halflife
        if new and iso.halflife > self.corrTime:
            self.subChain.append(iso)

        for dk in state.decays:
            daughter = self.chain(dk.daughterZ,dk.daughterA,dk.daughterExcitation)
            tran = Transition(iso,daughter,dk.radiation,fraction=dk.fraction)

        return iso

if '__main__' == __name__:
    import sys, data
    d = G4DataChain(int(sys.argv[1]),int(sys.argv[2]))
    print 'Built chain'
    pp = data.PrettyPrint(d.head)
    print pp 

---