GenData.py

Go to the documentation of this file.

#!/usr/bin/env python
'''
Usage:
'''

from UserTagging.UserTaggingAlg import UserTaggingAlg
from UserTagging.Models import Tag, Data, Int, Float, IntArray, FloatArray
from GaudiPython import SUCCESS, FAILURE
from GaudiPython import gbl, loaddict
from DybPython.Util import irange
import GaudiKernel.SystemOfUnits as units
import PyCintex
 
loaddict("libCLHEPRflx")
loaddict("libHepMCRflx")

class GenData(UserTaggingAlg):
    '''Add Generator Info to UserData'''
        
    # ----------------------------------------    
    def initTagList(self):
        # for details see doc-6020
        #
        # generator info, size of nGen
        # positon from first vtx
        # nVtx is the number of 'real primary' particles of each generator
        # IBD event has nVtx = 1 : nu_e_bar
        # Radioactive event nVtx > 1 for correlated decay events
        # Muon event nVtx = 1 : one out-going muon
        #
        # the first incoming particles of each vtx info, size of nVtxTotal, or 0 (muon)
        # for IBD, pid = [-12] for nu_e
        #          KE = KE for nu_e
        # for radioactive event, pid = [parent nuclei]
        # for radioactive event, t[j] - t[i] tells the correlation time
        # nOut[i] is the number of out-going particles at each primary vertex
        # nOut[i] == 2 for IBD 
        # nIn[i] should be 1, but 0 for muons
        #
        # outgoing particles info, size of nOutTotal
        # for IBD, out_pid = [-11, 2112]
        # for radioactive event, out_pid = [daughter nucliei+decay particles]
        self.addTag('Dummy', ''
            ).addData('GenData', '/Event/UserData/General/GenData'
            ).addInt('nGen', 'nVtxTotal', 'nOutTotal'
            ).addIntArray('genType', 'nVtx'
            ).addFloatArray('genX', 'genY', 'genZ'
            ).addIntArray('nIn', 'nOut', 'in_pid', 'out_pid'
            ).addFloatArray('in_x', 'in_y', 'in_z', 'in_t', 'in_KE', 'in_px', 'in_py', 'in_pz'
            ).addFloatArray('out_x', 'out_y', 'out_z', 'out_t', 'out_KE', 'out_px', 'out_py', 'out_pz')
               
        self.target_de_name = '/dd/Structure/AD/db-ade1/db-sst1/db-oil1'
        
        self.genTypes = {
            "IBD_gds":11, "IBD_lso":12, "IBD_oil":13, "IBD_acrylic":14,           
            "U238_gds":21, "U238_lso":22, "U238_PMT":23, "U238_sst":24,
            "Th232_gds":31, "Th232_lso":32, "Th232_PMT":33, "Th232_sst":34,
            "K40_gds":41, "K40_lso":42, "K40_PMT":43, "K40_sst":44,
            "Co60_sst":54,
        }
        
    # ----------------------------------------
    def check(self, evt):

        readoutHdr = evt['/Event/Readout/ReadoutHeader']
        if not readoutHdr:
            self.warning('cannot find readoutHdr')
            return
        
        genHdrs = readoutHdr.findHeaders(gbl.DayaBay.GenHeader.classID())
        self.SaveGenData(genHdrs)
        
        self.tagIt('Dummy')
        
    # ----------------------------------------
    def SaveGenData(self, genHdrs):

        myData = self.getTag('Dummy').getData('GenData')

        de = self.getDet(self.target_de_name)
        if not de:
            self.info('Failed to get DE' + self.target_de_name)
            return FAILURE
        Gaudi = PyCintex.makeNamespace('Gaudi')
                                
        # number of generators associated with this readout
        # nGen should be one except for pile-up events    
        nGen = len(genHdrs) 
        # total number of 'real primary' particles (vtx)
        nVtxTotal = 0
        # total number of outgoing particles
        nOutTotal = 0
 
        for genHdr in genHdrs:
            genName = genHdr.generatorName()
            self.debug("genName:" + genName)
            myData.append('genType', self.genTypes.get(genName, 0))
            
            nVtx = 0
            genEvt = genHdr.event()
            for vtx in irange(genEvt.vertices_begin(), genEvt.vertices_end()):
                # do stuff with the vertex
                nVtx += 1
                nVtxTotal += 1
                position = vtx.position()
                genGlbPoint = Gaudi.XYZPoint(
                    position.x(), 
                    position.y(),
                    position.z()
                )
                genLclPoint = de.geometry().toLocal(genGlbPoint)
                myData.append('in_x', genLclPoint.x()/units.mm)
                myData.append('in_y', genLclPoint.y()/units.mm)
                myData.append('in_z', genLclPoint.z()/units.mm)
                myData.append('in_t', vtx.position().t()/units.microsecond)
                
                if nVtx == 1:
                    myData.append('genX', genLclPoint.x()/units.mm)
                    myData.append('genY', genLclPoint.y()/units.mm)
                    myData.append('genZ', genLclPoint.z()/units.mm)
                                
                nIn = 0    
                for particle in irange(vtx.particles_in_const_begin(),
                                      vtx.particles_in_const_end()):
                    # do stuff with the incoming particles
                    nIn += 1
                    if nIn == 1:                  
                        myData.append('in_pid', particle.pdg_id())
                        momentum = particle.momentum()
                        myData.append('in_KE', (momentum.e() - momentum.m())/units.MeV)
                        myData.append('in_px', momentum.px()/units.MeV)
                        myData.append('in_py', momentum.py()/units.MeV)
                        myData.append('in_pz', momentum.pz()/units.MeV)
                myData.append('nIn', nIn)
                
                nOut = 0
                for particle in irange(vtx.particles_out_const_begin(),
                                       vtx.particles_out_const_end()):
                    # do stuff with the outgoing particles
                    nOut += 1
                    nOutTotal += 1
                    myData.append('out_pid', particle.pdg_id())
                    momentum = particle.momentum()
                    myData.append('out_KE', (momentum.e() - momentum.m())/units.MeV)
                    myData.append('out_px', momentum.px()/units.MeV)
                    myData.append('out_py', momentum.py()/units.MeV)
                    myData.append('out_pz', momentum.pz()/units.MeV)
                                 
                myData.append('nOut', nOut)
                    
            myData.append('nVtx', nVtx)

        # set more user data
        myData.set("nGen", nGen
            ).set("nVtxTotal", nVtxTotal
            ).set("nOutTotal", nOutTotal)
        
        # print 'genX: ', myData.get('genX')
            
#####  Job Configuration for nuwa.py ########################################

def configure(argv=[]):
    pass

def run(app):
    myAlg = GenData("UserData::GenData")
    app.addAlgorithm(myAlg)
    pass