In This Package:

ACUNeutronCaptureVertex.py

Go to the documentation of this file.

#!/usr/bin/env python
# 
# Usage :
#       nuwa.py -n -1 -m "CaptureVertex -a CalibrationAxis[A,B,C] -s Scale[>0]" SimAndRecData.root
#
# Comments :
# Scale = 43200 / The Number of the Neutrons Generated at the Specified ACU 
# 43200 is the number of neutrons generated at the ACU per day,
# so multiplying a number by "Scale" gives the number per day. 

# Load DybPython
import math
from DybPython.DybPythonAlg import DybPythonAlg
from GaudiPython import SUCCESS, FAILURE
from GaudiPython import gbl, loaddict
from DybPython.Util import irange

TH1F = gbl.TH1F
TH2F = gbl.TH2F

loaddict("libCLHEPRflx")
loaddict("libHepMCRflx")
Detector = gbl.DayaBay.Detector
AdPmtSensor = gbl.DayaBay.AdPmtSensor
ServiceMode = gbl.ServiceMode
ReconStatus = gbl.ReconStatus


class CaptureVertexAlg(DybPythonAlg):

    #Upload the list of possible capture targets
    Isotopes = [[1,1,'1H'],[12,6,'12C'],[13,6,'13C'],[14,7,'14N'],
                [16,8,'16O'],[28,14,'28Si'],[31,15,'31P'],[32,16,'32S'],
                [35,17,'35Cl'],[37,17,'37Cl'],[50,24,'50Cr'],[52,24,'52Cr'],
                [53,24,'53Cr'],[54,24,'54Cr'],[55,25,'55Mn'],[54,26,'54Fe'],
                [56,26,'56Fe'],[57,26,'57Fe'],[58,26,'58Fe'],[58,28,'58Ni'],
                [60,28,'60Ni'],[62,28,'62Ni'],[107,47,'107Ag'],[155,64,'155Gd'],[157,64,'157Gd']]

    def captureIsotope(instance, a,z,reconE,isotopes = Isotopes):
        for i in range(0,len(isotopes)):
            if ((isotopes[i][0] == a) & (isotopes[i][1] == z)):
                instance.stats["/file0/hcap/isotopes"].Fill(isotopes[i][2],scale)
                if reconE >= 6:
                    instance.stats["/file0/hcap/isotopes6MeV"].Fill(isotopes[i][2],scale)

    def makeIsotopeHistogram(instance, hist,isotopes = Isotopes):
        hist.GetXaxis().SetTitle("Isotopes")
        hist.GetYaxis().SetTitle("Neutron Capture Events/Day")
        # Set Bin Labels for Isotopes
        for i in range(0,len(isotopes)):
            hist.GetXaxis().SetBinLabel(2 + i,isotopes[i][2])

    def __init__(self,name):
        DybPythonAlg.__init__(self,name)
        return

    def initialize(self):
        status = DybPythonAlg.initialize(self)
        if status.isFailure(): return status
        self.info("initializing")

        self.cableSvc = self.svc('ICableSvc','StaticCableSvc')
        if self.cableSvc == None:
            self.error("Failed to get StaticCableSvc")
            return FAILURE

        ###Generate Histograms ##
        ## Hadron Capture events

        hist = TH1F("Events","Total Number of Neutron Capture Events" +
                " for Each Neutron Generated in ACU  "+acuAxis,500, 0.0,5.0)
        hist.GetXaxis().SetTitle("Neutron Captures")
        hist.GetYaxis().SetTitle("Events")
        self.stats["/file0/hcap/nHcap"] = hist
        
        hist = TH1F("Isotopes","The Capture Target for Neutrons"+
                " at ACU  " + acuAxis + " that Generate Signals In AD"
                ,26,0.0,26.0)
        self.makeIsotopeHistogram(hist)
        self.stats["/file0/hcap/isotopes"] = hist

        hist = TH1F("Isotopes","The Capture Target for Neutrons"+
                " at ACU  " + acuAxis + " with Reconstructed E > 6MeV"
                ,26,0.0,26.0)
        self.makeIsotopeHistogram(hist)
        self.stats["/file0/hcap/isotopes6MeV"] = hist

        hist = TH2F("EGamma","Z of the Capture Target vs Total Gamma Energy" +
                " for the Neutrons at ACU  " + acuAxis +
                " that Generate Signals in AD" ,500,0.0,12.0,81,-0.5,80.5 )
        hist.GetXaxis().SetTitle("E/Neutron (MeV)")
        hist.GetYaxis().SetTitle("Z")
        hist.SetMarkerStyle(2)
        self.stats["/file0/hcap/ZvsTotalGamma"] = hist

        hist = TH2F("EGamma","Z of the Capture Target vs Total Gamma Energy" + 
                " for the Neutrons at ACU  " + acuAxis + 
                "  with Reconstructed E > 6MeV" ,500,0.0,12.0,81,-0.5,80.5 )
        hist.GetXaxis().SetTitle("E/Neutron (MeV)")
        hist.GetYaxis().SetTitle("Z")
        hist.SetMarkerStyle(2)
        self.stats["/file0/hcap/ZvsTotalGamma6MeV"] = hist

        hist = TH2F("ZvsA", "Z vs A of target atoms for the Neutrons at ACU "+
                        acuAxis,500, -0.5, 210.5, 500,-0.5,80.5)
        hist.GetXaxis().SetTitle("A")
        hist.GetYaxis().SetTitle("Z")
        hist.SetMarkerStyle(2)
        self.stats["/file0/hcap/ZvsA"]= hist
        
        ##Gamma-ray energy
        hist = TH1F("EGamma","Individual Gamma-Ray Energy for" +
                " the Neutrons at ACU  " + acuAxis + 
                " that Generate Signals in AD",500,0.0,12.0)
        hist.GetXaxis().SetTitle("E (MeV)")
        hist.GetYaxis().SetTitle("Events/Day")
        self.stats["/file0/gamma/individual"] = hist

        hist = TH1F("EGamma","Individual Gamma Ray Energy for the Neutrons"
                + " at ACU  " + acuAxis + 
                " with Reconstructed E > 6MeV",500,0.0,12.0)
        hist.GetXaxis().SetTitle("E (MeV)")
        hist.GetYaxis().SetTitle("Events/Day")
        self.stats["/file0/gamma/individual6MeV"] = hist
        
        hist = TH1F("EGamma", "Total Gamma Ray Energy Emitted from the Capture "
                + "a Neutron Generated in ACU  " + acuAxis
                ,500,0.0,12.0)
        hist.GetXaxis().SetTitle("E/Neutron (MeV)")
        hist.GetYaxis().SetTitle("Events/Day")
        self.stats["/file0/gamma/totalCaptureGamma"] = hist
        
        hist = TH1F("EGamma", "Total Gamma Ray Energy Emitted from the Capture "
                        + "a Neutron Generated in ACU  " + acuAxis +
                        "with ReconEnergy > 6MeV"
                        ,500,0.0,12.0)
        hist.GetXaxis().SetTitle("E/Neutron (MeV)")
        hist.GetYaxis().SetTitle("Events/Day")
        self.stats["/file0/gamma/totalCaptureGamma6MeV"] = hist
        
        hist = TH1F("EGamma","Gamma Ray Energy from Inelastic Collisions"
                        + "of Neutrons Generated In ACU  " + 
                        acuAxis ,500,0.0,12.0)
        hist.GetXaxis().SetTitle("E (MeV)")
        hist.GetYaxis().SetTitle("Events/Day")
        self.stats["/file0/gamma/inelastic"] = hist


        
        return SUCCESS

    def execute(self):
        self.info("executing")
        evt = self.evtSvc()

        recHdr = evt["/Event/Rec/RecHeader"]
        if recHdr == None:
            self.error("Failed to get RecHeader")
            return FAILURE

        simHdr = evt["/Event/Sim/SimHeader"]
        if simHdr == None:
            self.error("Failed to get SimHeader")
            return FAILURE
        
        recResults = recHdr.recResults()
        recTrigger = recResults["AdSimple"]
        if recTrigger.positionStatus() == ReconStatus.kGood:
            if recTrigger.energyStatus() == ReconStatus.kGood:
                reconEnergy = recTrigger.energy()
                phis = simHdr.particleHistory()
                vtx = phis.vertexVector()
                nvtx = vtx.size()
                nHcap = 0       #Total number of hadron capture events
                for i in range(0,nvtx):
                    vsec = vtx[i].secondaries()
                    vproc = vtx[i].process().name()
                    PDG = vtx[i].track().track().particle()
                    totalGammaE = 0 #Total gamma energy for each hadron capture
                    
                    # Hadron capture events analysis
                    if PDG == 2112 and ((vproc == "HadronCapture") or (vproc == "nCapture")):
                        eachGammaE = 0
                        Z = 0
                        A = 0   
                        for j in range(0,vsec.size()):
                            if vsec[j].track().vertices().size() > 0:
                                vertex1 = vsec[j].track().vertices()[0]
                                pdg = vsec[j].track().particle()
                                eachGammaE = vertex1.kineticEnergy()
                                if pdg > 1000000000: #code is 10LZZZAAAI
                                    Z = (int(pdg/10000)%1000)
                                    # The atomic mass of the daughter atom has been increased by 1
                                    # so 1 is subtracted to get A of the parent atom
                                    A = (int(pdg/10)%1000)-1
                                    nHcap += 1
                                    self.captureIsotope(A,Z,reconEnergy)
                                    self.stats["/file0/hcap/ZvsA"].Fill(A,Z)
                                if pdg == 22:
                                    totalGammaE += eachGammaE
                        if ((Z > 0) & (A > 0)):
                            self.stats["/file0/hcap/ZvsTotalGamma"].Fill(totalGammaE,Z)
                            self.stats["/file0/gamma/totalCaptureGamma"].Fill(totalGammaE,scale)        
                            if reconEnergy > 6:
                                self.stats["/file0/hcap/ZvsTotalGamma6MeV"].Fill(totalGammaE,Z)
                                self.stats["/file0/gamma/totalCaptureGamma6MeV"].Fill(totalGammaE,scale)        
                    
                    # Loop over all the vertices as well to get the energy of all gammas                        
                    for j in range(0,vsec.size()):
                        eachGammaE = 0
                        if vsec[j].track().vertices().size() > 0:
                            vertex1 = vsec[j].track().vertices()[0]
                            pdg = vsec[j].track().particle()
                            eachGammaE = vertex1.kineticEnergy()
                            if pdg == 22:
                                self.stats["/file0/gamma/individual"].Fill(eachGammaE,scale)
                                if reconEnergy >= 6:
                                    self.stats["/file0/gamma/individual6MeV"].Fill(eachGammaE,scale)
                                if vertex1.process().name() == "NeutronInelastic":
                                    self.stats["/file0/gamma/inelastic"].Fill(eachGammaE,scale) 
                self.stats["/file0/hcap/nHcap"].Fill(nHcap)
        
        return SUCCESS
        
    def finalize(self):
        self.info("finalizing")
        status = DybPythonAlg.finalize(self)
        return status


#####  Job Configuration for nuwa.py ########################################
acuAxis = 'A'
scale = 1.0

def configure(argv = []):
    """Configure this module with ACU unit""" 
    import sys, getopt
    opts,args = getopt.getopt(argv,"a:s:;")
    for opt,arg in opts:
        if opt == "-a":
            global acuAxis
            acuAxis = arg
            print acuAxis
        if opt == "-s":
            global scale
            scale = float(arg)
            print scale


    from StatisticsSvc.StatisticsSvcConf import StatisticsSvc
    statsSvc = StatisticsSvc()
    statsSvc.Output ={"file0":"CaptureVertex.root"}
    import DataSvc
    DataSvc.Configure()
    
    return

def run(app):
    '''
    Configure and add an algorithm to job
    '''
    app.ExtSvc += ["StaticCableSvc", "StatisticsSvc"]
    captureVertexAlg = CaptureVertexAlg("MyCaptureVertex")
    app.addAlgorithm(captureVertexAlg)
    pass

---