In This Package:

Histogram.py

Go to the documentation of this file.

#!/usr/bin/env python

'''
A job option module and algorithm to histogram kinematics from GenDecay
'''

from DybPython.DybPythonAlg import DybPythonAlg
from GaudiPython import SUCCESS, FAILURE, loaddict
from ROOT import TH1F
from DybPython.Util import irange
import GaudiKernel.SystemOfUnits as units
import math

# Load HepMC classes
loaddict('libHepMCRflx')

class MyPythonAlg(DybPythonAlg):
    def __init__(self,name,default_path=""):
        DybPythonAlg.__init__(self,name)
        if default_path[-1] == '/': default_path = default_path[:-1]
        self.default_path = default_path
        self.data = {}
        return

    def initialize(self):
        status = DybPythonAlg.initialize(self)
        if status.isFailure(): return status
        self.statsSvc = self.svc('IStatisticsSvc','StatisticsSvc')
        if self.statsSvc == None:
            self.error("Failed to initialize StatisticsSvc.")
            return FAILURE        
        return SUCCESS

    def finalize(self):
        return DybPythonAlg.finalize(self)

    def _key2path(self,key):
        if key[0] != '/': 
            key = self.default_path + '/' + key
        return key

    def __setitem__(self,key,value):
        #self.info('MyPthonAlg.__setitem__(%s,%s)'%(key,value))
        key = self._key2path(key)
        if self.data.has_key(key):
            raise KeyError, 'key "%s" already set'%key
        self.statsSvc.put(key,value)
        self.data[key] = value
        return

    def __getitem__(self,key):
        return self.data[self._key2path(key)]

    pass

def pdgid2name(idn):
    _pdgid2name = {
        11: "e-",
        -11: "e+",
        22: "gamma", 
        1000020040: "alpha",
        1000802060: "Hg-206",
        1000812060: "Tl-206",
        1000822060: "Pb-206",
        1000812080: "Tl-208",
        1000822080: "Pb-208",
        1000812100: "Tl-210",
        1000822100: "Pb-210",
        1000832100: "Bi-210",
        1000842100: "Po-210",
        1000822120: "Pb-212",
        1000832120: "Bi-212",
        1000842120: "Po-212",
        1000822140: "Pb-214",
        1000832140: "Bi-214",
        1000842140: "Po-214",
        1000842160: "Po-216",
        1000842180: "Po-218",
        1000862200: "Rn-220",
        1000862220: "Rn-222",
        1000882240: "Ra-224",
        1000882260: "Ra-226",
        1000882280: "Ra-228",
        1000892280: "Ac-228",
        1000902280: "Th-228",
        1000902300: "Th-230",
        1000902320: "Th-232",
        1000902340: "Th-234",
        1000912340: "Pa-234",
        1000922340: "U-234",
        1000922380: "U-238",
        }
    try:
        return _pdgid2name[idn]
    except KeyError:
        return str(idn)
    
    

class DecayEnergy(MyPythonAlg):
    '''
    Histogram radioactive decays by parent
    '''

    class ParentHist:
        def __init__(self,particle,alg):
            self.particle = particle
            self.bin=len(alg.parentHist)
            self.alg=alg

            self.name = pdgid2name(particle.pdg_id())

            #print 'DecayEnergy: %d: %s'%(particle.pdg_id(), self.name)

            # Energy spectra
            self.Ealpha = TH1F(self.mkname('Ealpha'),"Alpha kinetic energy (MeV) [%s]"%self.mkname('Ealpha'),2000,0,10)
            alg[self.mkname('Ealpha')] = self.Ealpha
            self.Ebeta = TH1F(self.mkname('Ebeta'),"Beta kinetic energy (MeV) [%s]"%self.mkname('Ebeta'),2000,0,10)
            alg[self.mkname('Ebeta')] = self.Ebeta
            self.Egamma = TH1F(self.mkname('Egamma'),"Gamma kinetic energy (MeV) [%s]"%self.mkname('Egamma'),2000,0,10)
            alg[self.mkname('Egamma')] = self.Egamma
            self.Ecap = TH1F(self.mkname('Ecap'),"Electron capture energy (MeV) [%s]"%self.mkname('Ecap'),2000,0,10)
            alg[self.mkname('Ecap')] = self.Ecap

            # Timing
            self.Talpha = TH1F(self.mkname('Talpha'),"Parent-Alpha time (ns) [%s]"%self.mkname('Talpha'),1000,0,1e9)
            alg[self.mkname('Talpha')] = self.Talpha
            self.Tbeta = TH1F(self.mkname('Tbeta'),"Parent-Beta time (ns) [%s]"%self.mkname('Tbeta'),1000,0,1e9)
            alg[self.mkname('Tbeta')] = self.Tbeta
            self.Tgamma = TH1F(self.mkname('Tgamma'),"Parent-Gamma time (ns) [%s]"%self.mkname('Tgamma'),1000,0,1e9)
            alg[self.mkname('Tgamma')] = self.Tgamma
            self.Tecap = TH1F(self.mkname('Tecap'),"Parent-EC time (ns) [%s]"%self.mkname('Tecap'),1000,0,1e9)
            alg[self.mkname('Tecap')] = self.Tecap

            self.nDecay = TH1F(self.mkname('nDecay'),"Number of decay particles [%s]"%self.mkname('nDecay'),10,0,10)
            alg[self.mkname('nDecay')] = self.nDecay

            return

        def mkname(self,label):
            return self.name+'_'+label

        pass

    def __init__(self,name,dk_type):
        MyPythonAlg.__init__(self,name,'/file1/%s/%s'%(name,dk_type))
        self.dk_type = dk_type
        self.parentHist = {}    # map pdgid to ParentHist
        self.lastRefTime = 0.0
        return

    def initialize(self):
        status = MyPythonAlg.initialize(self)
        if status.isFailure(): return status
        self.info("initializing")
        self['ngamma'] = TH1F("ngamma","Number of Gammas",100,0,100)
        self['necap'] = TH1F("necap","Number of Electron Captures",100,0,100)
        self['alphabeta'] = TH1F("alphabeta","Alpha(<0) or Beta(>0)",2,-1,1)
        self['Ealpha'] = TH1F("Ealpha","Alpha kinetic energy (MeV)",2000,0,10)
        self['Ebeta'] = TH1F("Ebeta","Beta kinetic energy (MeV)",2000,0,10)
        self['Egamma'] = TH1F("Egamma","Gamma kinetic energy (MeV)",2000,0,10)
        self['Ecap'] = TH1F("Ecap","Energy from electron capture (MeV)",2000,0,10)
        self['deltaT'] = TH1F("deltaT","Time between parents (sec)",1000,0,100)
        '''
        histograms to add:  

        * parent nuclide specie
        * per parent alpha, beta & gamma spectra
        * time between parents
        * per parent time from parent to alpha, beta, gamma

        code to add:

        * GenEvent should hold daughter nuclides, take care with
          corelated decays, save only last daughter.

        '''
        return status

    def execute(self):
        evt = self.evtSvc()
        path = "/Event/Gen/GenHeader"
        hdr = evt[path]

        refTimestamp = hdr.timeStamp()
        refSeconds = refTimestamp.GetSeconds()
        #print 'refSeconds=',refSeconds,'self.lastRefTme',self.lastRefTime
        deltaRefSeconds = refSeconds - self.lastRefTime
        self.lastRefTime = refSeconds

        event = hdr.event()
        if event is None:
            error('Failed to get %s'%path)
            return FAILURE


        sigvtx = event.signal_process_vertex()
        if sigvtx is None: 
            error('No signal vertex')
            return FAILURE

        if sigvtx.particles_in_size() == 0:
            error('No input particles in signal process vertex')
            return FAILURE

        parent = sigvtx.particles_in_const_begin().__deref__()
        if parent is None:
            fatal('Unable to deref iterator')
            return FAILURE

        parent_pos = sigvtx.position()
        thisRelSeconds = parent_pos.t()/units.second

        # to fix: get times consistent between More and Gaudi side of things.

        dT = deltaRefSeconds + thisRelSeconds
        #print 'dT=',dT,'(sec) deltaRefSeconds=',deltaRefSeconds,'thisRelSeconds=',thisRelSeconds,'refSeconds=',refSeconds,'self.lastRefTime=',self.lastRefTime
        if dT <= 0.0: logT = 0.0
        else: logT = math.log10(dT)
        self['deltaT'].Fill(logT)

        try:
            parentHist = self.parentHist[parent.pdg_id()]
        except KeyError:
            parentHist = DecayEnergy.ParentHist(parent,self)
            self.parentHist[parent.pdg_id()] = parentHist
            pass            

        ndecay = 0
        ngamma = 0
        necap = 0
        #print "%d vertices:"%event.vertices_size()
        for vertex in irange(event.vertices_begin(),
                             event.vertices_end()):
            pos = vertex.position()
            #print "vertex position = (%f, %f, %f, %e s) %d particles:"%\
            #    (pos.x(), pos.y(), pos.z(), pos.t()/units.second, vertex.particles_out_size())

            deltaT = pos.t()/units.ns

            for particle in irange(vertex.particles_out_const_begin(),
                                   vertex.particles_out_const_end()):
                id = particle.pdg_id()
                mom = particle.momentum()
                mass = particle.generated_mass()
                #print "\t%s 4 momentum = (%f, %f, %f, %f)"%(pdgid2name(id), mom.px(), mom.py(), mom.pz(), mom.e())
                
                ndecay += 1

                ke = (mom.e()-mass)/units.MeV
                # Alphas
                if id == 1000020040: 
                    self['Ealpha'].Fill(ke)
                    parentHist.Ealpha.Fill(ke)
                    parentHist.Talpha.Fill(deltaT)
                    self['alphabeta'].Fill(-0.5)
                    #if ke > 7.0:
                    #    print "BIG ALPHA"
                    pass
                # Betas
                if id == 11 or id == -11: 
                    self['Ebeta'].Fill(ke)
                    parentHist.Ebeta.Fill(ke)
                    parentHist.Tbeta.Fill(deltaT)
                    self['alphabeta'].Fill(+0.5)
                    pass
                # Gammas
                if id == 22: 
                    self['Egamma'].Fill(ke)
                    parentHist.Egamma.Fill(ke)
                    parentHist.Tgamma.Fill(deltaT)
                    ngamma += 1
                    pass
                # EleCapture
                if id == 0:
                    self['Ecap'].Fill(ke)
                    parentHist.Ecap.Fill(ke)
                    parentHist.Tecap.Fill(deltaT)
                    necap += 1
                    pass
                continue
            continue
        #print ngamma,'gammas'
        self['ngamma'].Fill(ngamma)
        self['necap'].Fill(necap)
        parentHist.nDecay.Fill(ndecay)

        return SUCCESS

    def finalize(self):
        return MyPythonAlg.finalize(self)

    pass

def configure(argv=[]):
    try:
        histfile = argv[0]
    except IndexError:
        histfile="gendecay.root"

    from StatisticsSvc.StatisticsSvcConf import StatisticsSvc
    statsSvc = StatisticsSvc()
    statsSvc.Output ={"file1":histfile}
    return

def run(app):
    app.ExtSvc += ["StaticCableSvc", "StatisticsSvc"]

    dke = DecayEnergy("DecayEnergy","Chain")
    app.addAlgorithm(dke)
    return

---