In This Package:

ACUNeutronCapturePosition.py

Go to the documentation of this file.

#!/usr/bin/env python
#
# Usage :
#       nuwa.py -n -1 CapturePosition SimAndRecData.root
#
# 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

XYZPoint = gbl.Gaudi.XYZPoint

class CapturePositionAlg(DybPythonAlg):
    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

        # Make position histograms
        hist = TH2F("CapturePosition","Captured Particle z vs x",
                500,-2500.0,2500.0, 500, 1500.0, 4500.0)
        hist.GetXaxis().SetTitle("x [mm]")
        hist.GetYaxis().SetTitle("z [mm]")
        self.stats["/file0/position/CapturePositionAll"] = hist


        hist = TH2F("CapturePosition","Captured Particle z vs x",
                500,-2500.0,2500.0, 500, 1500.0 , 4500.0)
        hist.GetXaxis().SetTitle("x [mm]")
        hist.GetYaxis().SetTitle("z [mm]")
        self.stats["/file0/position/CapturePosition"] = hist

        hist = TH2F("CapturePosition","Capture Position for particles with E > 6MeV "
                ,500,-2500.0,2500.0, 500, 1500.0 , 4500.0)
        hist.GetXaxis().SetTitle("x [mm]")
        hist.GetYaxis().SetTitle("z [mm]")
        self.stats["/file0/position/CapturePosition6MeV"] = hist
    
        hist = TH2F("Reconstructed Position","Energy range 0~2MeV"
                ,500, -2500.0, 2500.0, 500, -500.0, 3500.0)
        hist.GetXaxis().SetTitle("x [mm]")
        hist.GetYaxis().SetTitle("z [mm]")
        self.stats["/file0/position/ReconPosition1"] = hist
        
        hist = TH2F("Reconstructed Position","Energy range 2~4MeV"
                ,500, -2500.0, 2500.0, 500, -500.0, 3500.0)
        hist.GetXaxis().SetTitle("x [mm]")
        hist.GetYaxis().SetTitle("z [mm]")
        self.stats["/file0/position/ReconPosition2"] = hist

        hist = TH2F("Reconstructed Position","Energy range 4~6MeV"
                ,500, -2500.0, 2500.0, 500, -500.0, 3500.0)
        hist.GetXaxis().SetTitle("x [mm]")
        hist.GetYaxis().SetTitle("z [mm]")
        self.stats["/file0/position/ReconPosition3"] = hist

        hist = TH2F("Reconstructed Position","Energy range 6~8MeV"
                ,500, -2500.0, 2500.0, 500, -500.0, 3500.0)
        hist.GetXaxis().SetTitle("x [mm]")
        hist.GetYaxis().SetTitle("z [mm]")
        self.stats["/file0/position/ReconPosition4"] = hist
        
        hist = TH2F("Reconstructed Position","Energy range > 8 MeV"
                ,500, -2500.0, 2500.0, 500, -500.0, 3500.0)
        hist.GetXaxis().SetTitle("x [mm]")
        hist.GetYaxis().SetTitle("z [mm]")
        self.stats["/file0/position/ReconPosition5"] = hist

        hist = TH1F("reconEnergy", "Reconstructed Energy",500,0.0,10.0)
        hist.GetXaxis().SetTitle("E (MeV)")
        hist.GetYaxis().SetTitle("Events")
        self.stats["/file0/energy/ReconEnergy"] = hist  
        
        return SUCCESS

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

        #The Capture Position of the Particles that Generate Signals in AD
        simHdr = evt["/Event/Sim/SimHeader"]
        if simHdr == None:
            self.error("Failed to get SimHeader")
            return FAILURE
        statsHdr = simHdr.unobservableStatistics()
        simStats = statsHdr.stats()

        x1 = simStats["xEnd_Trk1"].sum()
        y1 = simStats["yEnd_Trk1"].sum()
        z1 = simStats["zEnd_Trk1"].sum()

        globalPos1 = XYZPoint(x1, y1, z1)
        ad1 = self.detSvc("/dd/Structure/AD/db-oil1")
        localPos1 = ad1.geometry().toLocal(globalPos1)
        cx = localPos1.x()
        cz = localPos1.z()
        self.stats["/file0/position/CapturePositionAll"].Fill(cx,cz)

        recHdr = evt["/Event/Rec/RecHeader"]
        if recHdr == None:
            self.error("Failed to get RecHeader")
            return FAILURE
        recResults = recHdr.recResults()
        recTrigger = recResults["AdSimple"]
        if recTrigger.positionStatus() == ReconStatus.kGood:
            reconPos = recTrigger.position()
            rx = reconPos.x()
            rz = reconPos.z()
            self.stats["/file0/position/CapturePosition"].Fill(cx,cz)
            
        if recTrigger.energyStatus() == ReconStatus.kGood:
            reconEnergy = recTrigger.energy()
            self.stats["/file0/energy/ReconEnergy"].Fill(reconEnergy)
            if ((reconEnergy >= 0) & (reconEnergy < 2)):
                self.stats["/file0/position/ReconPosition1"].Fill(rx,rz)

            if ((reconEnergy >= 2) & (reconEnergy < 4)) :
                self.stats["/file0/position/ReconPosition2"].Fill(rx,rz)

            if ((reconEnergy >= 4) & (reconEnergy < 6)) :
                self.stats["/file0/position/ReconPosition3"].Fill(rx,rz)

            if ((reconEnergy >= 6) & (reconEnergy < 8)) :
                self.stats["/file0/position/ReconPosition4"].Fill(rx,rz)
                self.stats["/file0/position/CapturePosition6MeV"].Fill(cx,cz)
            if reconEnergy >= 8 :
                self.stats["/file0/position/ReconPosition5"].Fill(rx,rz)
                self.stats["/file0/position/CapturePosition6MeV"].Fill(cx,cz)

        return SUCCESS
        
    def finalize(self):
        self.info("finalizing")
        status = DybPythonAlg.finalize(self)
        return status


#####  Job Configuration for nuwa.py ########################################

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

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

---