MDC09a: /home/dayabay/installs/trunk_2011_04_11_opt/NuWa-trunk/dybgaudi/Production/MDC09a/python/MDC09a/chkIBD15/__init_

00001 #!/usr/bin/env python
00002 #
00003 # to check the quality of the IBD15 MC and to produce some basic distributions.
00004 # --- Wei, Apr 2009
00005 
00006 __all__ = ['plotIbdBasics']
00007 from GaudiPython.GaudiAlgs import GaudiAlgo
00008 from GaudiPython import SUCCESS, FAILURE
00009 from GaudiPython import gbl
00010 
00011 from DybPython.Util import irange
00012 from GaudiKernel import SystemOfUnits as units
00013 import ROOT
00014 
00015 TH1F = gbl.TH1F
00016 TH1D = gbl.TH1D
00017 TH2F = gbl.TH2F
00018 TH2D = gbl.TH2D
00019 
00020 class plotIbdBasics(GaudiAlgo):
00021 
00022     def __init__(self,name):
00023         GaudiAlgo.__init__(self,name)
00024         self.hist = {}
00025         # Counters
00026         self.SpillIn_onGd = 0
00027         self.SpillOut_onGd = 0
00028         self.SpillIn_onH = 0
00029         self.SpillOut_onH = 0
00030 
00031         return
00032     
00033     def initialize(self):
00034         print "Initializing the IBD basic ploter", self.name()
00035         status = GaudiAlgo.initialize(self)
00036         if status.isFailure(): return status
00037         self.target_de_name = '/dd/Structure/AD/db-ade1/db-sst1/db-oil1'
00038         self.gds_de_name = '/dd/Structure/AD/db-gds1'
00039         self.lso_de_name = '/dd/Structure/AD/db-lso1'
00040 
00041         # What services do you need?
00042         self.coorSvc = self.svc('ICoordSysSvc', 'CoordSysSvc')
00043         if not self.coorSvc:
00044             print 'Failed to get CoordSysSvc'
00045             return FAILURE
00046 
00047         self.histSvc = self.svc('ITHistSvc', 'THistSvc')
00048 
00049         self.hist["genRZ"] = TH2D("genRZ", "Generation Vertex R-Z", \
00050                                       100, 0.0, 6.190144, 100, -2.48, 2.48)
00051         status = self.histSvc.regHist('/file1/basics/genRZ', \
00052                                           self.hist["genRZ"])
00053         if status.isFailure(): return status
00054         
00055         self.hist["genXY"] = TH2D("genXY", "Generation Vertex X-Y", \
00056                                       100, -2.48, 2.48, 100, -2.48, 2.48)
00057         status = self.histSvc.regHist('/file1/basics/genXY', \
00058                                           self.hist["genXY"])
00059         if status.isFailure(): return status
00060         
00061         self.hist["HitTime"] = TH1F("HitTime", "Hit Time [ms]",
00062                                     500, 0.0, 2000)
00063         status = self.histSvc.regHist('/file1/basics/HitTime', 
00064                                       self.hist["HitTime"])
00065         if status.isFailure(): return status
00066 
00067         self.hist["pe_E"] = TH2F("pe_E", "pe vs visible E (e+ within GdLS)",
00068                                  100, 0, 10, 700, 0, 1400)
00069         status = self.histSvc.regHist('/file1/basics/pe_E', 
00070                                       self.hist["pe_E"])
00071         if status.isFailure(): return status
00072 
00073         self.hist["pe_E_inLS"] = TH2F("pe_E_inLS", \
00074                                         "pe vs visible E (e+ within LS)", \
00075                                         100, 0, 10, 700, 0, 1400)
00076         status = self.histSvc.regHist('/file1/basics/pe_E_inLS', 
00077                                       self.hist["pe_E_inLS"])
00078         if status.isFailure(): return status
00079 
00080         self.hist["pe_E_LS"] = TH2F("pe_E_LS", \
00081                                         "pe vs visible E (e+ in LS)", \
00082                                         100, 0, 10, 700, 0, 1400)
00083         status = self.histSvc.regHist('/file1/basics/pe_E_LS', 
00084                                       self.hist["pe_E_LS"])
00085         if status.isFailure(): return status
00086 
00087         self.hist["pe_E_all"] = TH2F("pe_E_all", "pe vs visible E (all e+)",
00088                                      100, 0, 10, 700, 0, 1400)
00089         status = self.histSvc.regHist('/file1/basics/pe_E_all', 
00090                                       self.hist["pe_E_all"])
00091         if status.isFailure(): return status
00092         
00093         self.hist["nHits_LS"] = TH1F("nHits_LS", \
00094                                          "nCap Hits in LS", \
00095                                          100, 0, 1400)
00096         status = self.histSvc.regHist('/file1/basics/nHits_LS', 
00097                                       self.hist["nHits_LS"])
00098         if status.isFailure(): return status
00099         
00100         self.hist["nHits_MO"] = TH1F("nHits_MO", \
00101                                          "nCap Hits in MO", \
00102                                          100, 0, 1400)
00103         status = self.histSvc.regHist('/file1/basics/nHits_MO', 
00104                                       self.hist["nHits_MO"])
00105         if status.isFailure(): return status
00106 
00107         self.hist["nHits"] = TH1F("nHits", "Neutron Capture Hits",
00108                                   100, 0, 1400)
00109         status = self.histSvc.regHist('/file1/basics/nHits', 
00110                                       self.hist["nHits"])
00111         if status.isFailure(): return status
00112 
00113         self.hist["nHits_all"] = TH1F("nHits_all", "Neutron Capture Hits",
00114                                   100, 0, 1400)
00115         status = self.histSvc.regHist('/file1/basics/nHits_all', 
00116                                       self.hist["nHits_all"])
00117         if status.isFailure(): return status
00118 
00119         self.hist["nCapPos"] = TH1F("nCapPos", "Neutron Capture Position",
00120                                     310, 0, 6.190144)
00121         status = self.histSvc.regHist('/file1/basics/nCapPos',
00122                                       self.hist["nCapPos"])
00123         if status.isFailure(): return status
00124 
00125         self.hist["eDepInGdLS"] = TH1F("eDepInGdLS", "Deposited Energy [MeV]",
00126                                        20, 0, 20)
00127         status = self.histSvc.regHist('/file1/basics/eDepInGdLS', 
00128                                       self.hist["eDepInGdLS"])
00129         if status.isFailure(): return status
00130 
00131         self.hist["eDepInLS"] = TH1F("eDepInLS", "Deposited Energy [MeV]",
00132                                        20, 0, 20)
00133         status = self.histSvc.regHist('/file1/basics/eDepInLS', 
00134                                       self.hist["eDepInLS"])
00135         if status.isFailure(): return status
00136 
00137         self.hist["drift_GdLS"] = TH1F("drift_GdLS",
00138                                       "Neutron Drift Distance in GdLS [cm]",
00139                                        200, 0, 50)
00140         status = self.histSvc.regHist('/file1/basics/drift_GdLS', 
00141                                       self.hist["drift_GdLS"])
00142         if status.isFailure(): return status
00143 
00144         self.hist["drift_LS"] = TH1F("drift_LS",
00145                                       "Neutron Drift Distance in LS [cm]",
00146                                        200, 0, 50)
00147         status = self.histSvc.regHist('/file1/basics/drift_LS', 
00148                                       self.hist["drift_LS"])
00149         if status.isFailure(): return status
00150         
00151         self.hist["time_GdLS"] = TH1F("time_GdLS",
00152                                       "Neutron Capture time in GdLS [ms]",
00153                                        400, 0, 400)
00154         status = self.histSvc.regHist('/file1/basics/time_GdLS', 
00155                                       self.hist["time_GdLS"])
00156         if status.isFailure(): return status
00157 
00158         self.hist["time_LS"] = TH1F("time_LS",
00159                                       "Neutron Capture Time in LS [ms]",
00160                                        400, 0, 2000)
00161         status = self.histSvc.regHist('/file1/basics/time_LS', 
00162                                       self.hist["time_LS"])
00163         if status.isFailure(): return status
00164 
00165         return SUCCESS
00166 
00167     def execute(self):
00168         print "Executing plotIbdBasics", self.name()
00169         evt = self.evtSvc()
00170         simhdr = evt['/Event/Sim/SimHeader']
00171 #        print "SimHeader: ", simhdr
00172         if simhdr == None:
00173             print "No SimHeader in this ReadOut"
00174             return SUCCESS
00175 
00176         det = self.detSvc(self.target_de_name)
00177         det_gds = self.detSvc(self.gds_de_name)
00178         det_lso = self.detSvc(self.lso_de_name)
00179 
00180         # Unobservables
00181         statshdr = simhdr.unobservableStatistics()
00182         stats = statshdr.stats()
00183 
00184         PID_trk1 = stats["pdgId_Trk1"].sum()
00185         PID_trk2 = stats["pdgId_Trk2"].sum()
00186 
00187         if PID_trk1 != -11 or PID_trk2 != 2112:
00188             print "PID of track 1 is", PID_trk1
00189             print "PID of track 2 is", PID_trk2
00190             print "Not an IBD event."
00191             return SUCCESS
00192 
00193         tGen = stats["t_Trk2"].sum()
00194         xGen = stats["x_Trk2"].sum()
00195         yGen = stats["y_Trk2"].sum()
00196         zGen = stats["z_Trk2"].sum()
00197         
00198         tCap = stats["tEnd_Trk2"].sum()
00199         xCap = stats["xEnd_Trk2"].sum()
00200         yCap = stats["yEnd_Trk2"].sum()
00201         zCap = stats["zEnd_Trk2"].sum()
00202         
00203         # Get underlying DE object
00204         de = self.getDet(self.target_de_name)
00205         de_lso = self.getDet(self.lso_de_name)
00206         de_gds = self.getDet(self.gds_de_name)
00207         if not de:
00208             print 'Failed to get DE',self.target_de_name
00209             return FAILURE
00210         
00211         if not de_lso:
00212             print 'Failed to get DE',self.lso_de_name
00213             return FAILURE
00214         
00215         if not de_gds:
00216             print 'Failed to get DE',self.gds_de_name
00217             return FAILURE
00218         
00219         # Get the AD coordinates of the vertexes
00220         import PyCintex
00221         Gaudi = PyCintex.makeNamespace('Gaudi')
00222         genGlbPoint = Gaudi.XYZPoint(xGen, yGen, zGen)
00223         capGlbPoint = Gaudi.XYZPoint(xCap, yCap, zCap)
00224 #        point = de.geometry().toGlobal(point)
00225         genLclPoint = de.geometry().toLocal(genGlbPoint)
00226         capLclPoint = de.geometry().toLocal(capGlbPoint)
00227         genLclPointLso = de_lso.geometry().toLocal(genGlbPoint)
00228         capLclPointLso = de_lso.geometry().toLocal(capGlbPoint)
00229         genLclPointGds = de_gds.geometry().toLocal(genGlbPoint)
00230         capLclPointGds = de_gds.geometry().toLocal(capGlbPoint)
00231 #        print 'In global coordinate [',gpoint.x(),gpoint.y(),gpoint.z(),']'
00232 
00233         ndrift = ROOT.TVector3(xCap-xGen, yCap-yGen, zCap-zGen)
00234         
00235         capTime = tCap - tGen
00236         capDis = ndrift.Mag()
00237 
00238         R2 = genLclPoint.x()/units.meter * genLclPoint.x()/units.meter + \
00239             genLclPoint.y()/units.meter * genLclPoint.y()/units.meter
00240         
00241         self.hist["genRZ"].Fill(R2, genLclPoint.z()/units.meter)
00242 
00243         self.hist["genXY"].Fill(genLclPoint.x()/units.meter,genLclPoint.y()/units.meter)
00244 
00245         # Find the interesting volumes
00246         genDE = self.coorSvc.coordSysDE(genGlbPoint)
00247         capDE = self.coorSvc.coordSysDE(capGlbPoint)
00248         if not genDE:
00249             print 'Failed to find coordinate system DE for generation'\
00250                 '[',genGlbPoint.x(),genGlbPoint.y(),genGlbPoint.z(),']'
00251             print 'Local: [',genLclPoint.x(),genLclPoint.y(),genLclPoint.z(),']'
00252             return FAILURE
00253         else:
00254             gendmvol = genDE.geometry().belongsToPath(genGlbPoint,-1)
00255 
00256         if not capDE:
00257             print 'Failed to find coordinate system DE for capture'\
00258                 '[',capGlbPoint.x(),capGlbPoint.y(),capGlbPoint.z(),']'
00259             return FAILURE
00260         else:
00261             capdmvol = capDE.geometry().belongsToPath(capGlbPoint,-1)
00262 
00263         import re
00264         genDM = re.split('/', gendmvol).pop()
00265         capDM = re.split('/', capdmvol).pop()
00266         print "Generated in ", genDM
00267         print "Captured in ", capDM
00268 
00269         positronHits = 0
00270         neutronHits = 0
00271         positronTimeCut = 500.
00272         simhits = simhdr.hits()
00273         for detector in simhits.hitDetectors():
00274             hitCollection = simhits.hitsByDetector(detector)
00275             if hitCollection == None:
00276                 print "No hits in ", detector
00277             hits = hitCollection.collection()
00278             for hit in hits:
00279 #                print " PMT", hit.sensDetId(), "hit @ time: ", hit.hitTime()
00280                 self.hist["HitTime"].Fill(hit.hitTime()/units.microsecond)
00281                 if hit.hitTime()/units.nanosecond<positronTimeCut and \
00282                         hit.hitTime()/units.nanosecond < \
00283                         capTime/units.nanosecond:
00284                     positronHits += 1
00285                 else:
00286                     neutronHits += 1
00287 
00288         # Visible energy
00289         vis_trk1 = 1.022 + stats['ke_Trk1'].sum()/units.MeV
00290         self.hist["pe_E_all"].Fill(vis_trk1, positronHits)
00291         if genDM ==  'db-gds1':
00292             self.hist["pe_E"].Fill(vis_trk1, positronHits)
00293         if genDM ==  'db-lso1':
00294             self.hist["pe_E_LS"].Fill(vis_trk1, positronHits)
00295         if re.search('db-lso1',gendmvol):
00296             self.hist["pe_E_inLS"].Fill(vis_trk1, positronHits)
00297 
00298         capTarget = stats["capTarget"].sum()
00299         print "The capture target is ", capTarget
00300 
00301         if capTarget == 1 and \
00302                 capLclPoint.z()/units.m < 1. and \
00303                 capLclPoint.z()/units.m > -1.:
00304             self.hist["nCapPos"].Fill(R2)
00305 
00306         self.hist["nHits_all"].Fill(neutronHits)
00307 
00308         if capDM == 'db-lso1':
00309             self.hist["nHits_LS"].Fill(neutronHits)
00310 
00311         if capDM == 'db-oil1':
00312             self.hist["nHits_MO"].Fill(neutronHits)
00313 
00314         if genDM == 'db-gds1' and capDM == 'db-gds1':
00315             self.hist["nHits"].Fill(neutronHits)
00316             self.hist["time_GdLS"].Fill(capTime/units.microsecond)
00317             self.hist["drift_GdLS"].Fill(capDis/units.cm)
00318 
00319         if genDM == 'db-lso1' and capDM == 'db-lso1':
00320             self.hist["time_LS"].Fill(capTime/units.microsecond)
00321             self.hist["drift_LS"].Fill(capDis/units.cm)
00322 
00323         if genDM == 'db-lso1' and capDM == 'db-oil1':
00324             self.SpillOut_onH += 1
00325 
00326         if genDM == 'db-oil1' and capDM == 'db-lso1':
00327             self.SpillIn_onH += 1
00328 
00329         if genDM == 'db-gds1' and capDM == 'db-lso1':
00330             self.SpillOut_onGd += 1
00331 
00332         if genDM == 'db-lso1' and capDM == 'db-gds1':
00333             self.SpillIn_onGd += 1
00334 
00335         eDepInGdLS = stats["EDepInGdLS"].sum()
00336         self.hist["eDepInGdLS"].Fill(eDepInGdLS/units.MeV)
00337         eDepInLS = stats["EDepInLS"].sum()
00338         self.hist["eDepInLS"].Fill(eDepInLS/units.MeV)
00339         
00340         return SUCCESS
00341     
00342     def finalize(self):
00343 
00344         print "Spill-In of on Gd n: ", self.SpillIn_onGd
00345         print "Spill-Out of n from Gd-LS: ", self.SpillOut_onGd
00346         print "Spill-In of n from OIL to LSO: ", self.SpillIn_onH
00347         print "Spill-Out of n from LSO to OIL: ", self.SpillOut_onH
00348 
00349         print "Finalizing ", self.name()
00350         status = GaudiAlgo.finalize(self)
00351         return status
00352 
00353 def configure():
00354     from DetHelpers.DetHelpersConf import CoordSysSvc
00355     from GaudiSvc.GaudiSvcConf import THistSvc
00356 #    from GaudiSvc.GaudiSvcConf import DetectorDataSvc
00357 
00358     histsvc = THistSvc()
00359     histsvc.Output = ["file1 DATAFILE='IbdBasicPlots.root' OPT='RECREATE' TYP='ROOT' "]
00360 
00361     coorSvc = CoordSysSvc()
00362     coorSvc.OutputLevel = 1
00363 
00364     from Gaudi.Configuration import ApplicationMgr
00365     theApp = ApplicationMgr()
00366     theApp.ExtSvc.append(coorSvc)
00367 
00368     return
00369 
00370 def run(app):
00371     '''Configure and add this algorithm to job'''
00372 #    from Gaudi.Configuration import ApplicationMgr
00373 #    app = ApplicationMgr()
00374     app.ExtSvc += ["THistSvc"]
00375     plotBasics = plotIbdBasics("myBasics")
00376     app.addAlgorithm(plotBasics)
00377     pass
In This Package:

init.py

In This Package:

__init__.py

init.py
