MDC09a: /home/dayabay/installs/trunk_2011_04_11_opt/NuWa-trunk/dybgaudi/Production/MDC09a/python/MDC09a/chkIBD/AdPerformance.py Source File

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 PyCintex
00014 import re
00015 import ROOT
00016 
00017 TH1F = gbl.TH1F
00018 TH1D = gbl.TH1D
00019 TH2F = gbl.TH2F
00020 TH2D = gbl.TH2D
00021 TH3F = gbl.TH3F
00022 TH3D = gbl.TH3D
00023 TParameter = gbl.TParameter
00024 
00025 class plotIbdBasics(GaudiAlgo):
00026 
00027     def __init__(self,name):
00028         GaudiAlgo.__init__(self,name)
00029 
00030         return
00031     
00032     def initialize(self):
00033         print "Initializing the IBD basic ploter", self.name()
00034         status = GaudiAlgo.initialize(self)
00035         if status.isFailure(): return status
00036 
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.statsSvc = self.svc('IStatisticsSvc','StatisticsSvc')
00048         if self.statsSvc == None:
00049             self.error("Failed to initialize IStat service.")
00050             return FAILURE        
00051 
00052         # output file
00053         self.outputstat = open('IbdStat.txt', 'w')
00054         
00055         # Counters
00056         self.counter = {}
00057         self.counter["nInGdLS"] = 0.
00058         
00059         self.counter["nInIAV"] = 0.
00060         self.counter["nInLSO"] = 0.
00061         self.counter["nInOAV"] = 0.
00062         self.counter["nInMO"] = 0.
00063         
00064         
00065         self.counter["nCap"] = 0.
00066         self.counter["nGen"] = 0.
00067         
00068         self.counter["SpillIn_onGd"] = 0.
00069         self.counter["LSSpillIn_onGd"] = 0.
00070         self.counter["IavSpillIn_onGd"] = 0.
00071         self.counter["SpillOut"] = 0.
00072         self.counter["OilSpillIn_onH"] = 0.
00073         self.counter["OavSpillIn_onH"] = 0.
00074         self.counter["SpillOut_onH"] = 0.
00075         
00076         self.counter["onGdCap"] = 0.           # total on-Gd cap
00077         self.counter["onGdGdsCap"] = 0.        # total on-Gd cap in GDS
00078 
00079         # Add ratios of capture target in GdLS, Liang Zhan 
00080         ###########################################################
00081         self.counter["onHCapGds"] = 0.        # total on-H cap in GDS
00082         self.counter["onCCapGds"] = 0.        # total on-C cap in GDS
00083         self.counter["onOtherCapGds"] = 0.        # total other cap in GDS
00084         ###########################################################
00085 
00086         # Add coincidence time cut, Liang Zhan
00087         ################################################
00088         self.counter["nGdsCapTime"] = 0.
00089         self.counter["nGdsCapTimePassCut"] = 0.
00090         self.counter["nGdsCapTimeL"] = 0.
00091         self.counter["nGdsCapTimeUncertainL"] = 0.
00092         self.counter["nGdsCapTimeH"] = 0.
00093         self.counter["nGdsCapTimeUncertainH"] = 0.
00094 
00095         self.counter["nLSCapTime"] = 0.
00096         self.counter["nLSCapTimePassCut"] = 0.
00097         self.counter["nLSCapTimeL"] = 0.
00098         self.counter["nLSCapTimeUncertainL"] = 0.
00099         self.counter["nLSCapTimeH"] = 0.
00100         self.counter["nLSCapTimeUncertainH"] = 0.
00101         ################################################
00102         
00103         self.counter["onGdCapPassCut"] = 0.    # total on-Gd cap pass Gd cut
00104         self.counter["onGdGdsCapPassCut"] = 0. # total on-Gd cap in GDS pass Gd cut
00105         self.counter["onGdGdsCapPassEdepCut"] = 0. # total on-Gd cap in GDS pass Gd cut (6 MeV true deposit energy)
00106         self.counter["nPassGdCut"] = 0.        # total pass Gd cut
00107 
00108         self.counter["GdCapDetEff"] = 0.
00109         self.counter["GdCapDetEffUncertain"] = 0.
00110 
00111         self.counter["onGdCapCutUncertain"] = 0.     
00112         self.counter["nGdCutUncertain"] = 0.
00113         self.counter["onGdGdsCapCutUncertain"] = 0.
00114 
00115         self.counter["positronPassCut"] = 0.
00116         self.counter["positronCutUncertain"] = 0.
00117         
00118         self.counter["nonGdPassGdCut"] = 0.
00119 
00120         self.counter["onHCap"] = 0.
00121         self.counter["onHCapWithinLS"] = 0.
00122         self.counter["onHCapWithinLSPassCut"] = 0.
00123         self.counter["onHCapWithinLSCutUncertainL"] = 0.
00124         self.counter["onHCapWithinLSCutUncertainU"] = 0.
00125 
00126         self.counter["HCapDetEff"] = 0.
00127         self.counter["HCapDetEffUncertainL"] = 0.
00128         self.counter["HCapDetEffUncertainU"] = 0.
00129 
00130         self.counter["HCapCanoDetEff"] = 0.
00131         self.counter["HCapCanoDetEffUncertainL"] = 0.
00132         self.counter["HCapCanoDetEffUncertainU"] = 0.
00133 
00134         self.counter["HCapDetContamination"] = 0.
00135         self.counter["HCapDetEffContaminationUncertainL"] = 0.
00136         self.counter["HCapDetEffContaminationUncertainU"] = 0.
00137 
00138         self.counter["nonHCapWithinLSPassHCut"] = 0.
00139         self.counter["nonHCapWithinLSHCutUncertainL"] = 0.
00140         self.counter["nonHCapWithinLSHCutUncertainU"] = 0.
00141 
00142         self.counter["onHCapBeyondLS"] = 0.
00143         self.counter["onHCapBeyondLSPassCut"] = 0.
00144         self.counter["onHCapBeyondLSCutUncertainL"] = 0.
00145         self.counter["onHCapBeyondLSCutUncertainU"] = 0.
00146 
00147         self.counter["nonHCapBeyondLSPassHCut"] = 0.
00148         self.counter["nonHCapBeyondLSHCutUncertainL"] = 0.
00149         self.counter["nonHCapBeyondLSHCutUncertainU"] = 0.
00150 
00151         self.counter["onHCapOavMoPassCut"] = 0.
00152         self.counter["onHCapOavMoCutUncertainL"] = 0.
00153         self.counter["onHCapOavMoCutUncertainU"] = 0.
00154         self.counter["onHCapLS"] = 0.
00155         self.counter["onHCapMO"] = 0.
00156         self.counter["onHCapOAV"] = 0.
00157         
00158         self.counter["onHCapPassCut"] = 0.
00159         self.counter["onHCapCutUncertainL"] = 0.
00160         self.counter["onHCapCutUncertainU"] = 0.
00161         
00162         self.counter["nPassHCut"] = 0.
00163         self.counter["nHCutUncertainL"] = 0.
00164         self.counter["nHCutUncertainU"] = 0.
00165 
00166         self.counter["nGdPassHCut"] = 0.
00167         self.counter["nGdHCutUncertainL"] = 0.
00168         self.counter["nGdHCutUncertainU"] = 0.
00169 
00170         # Histograms
00171         self.hist = {}
00172         
00173         self.hist["genRZ"] = TH2F("genRZ", "Generation Vertex R-Z",
00174                                   100, 0.0, 6.190144, 100, -2.48, 2.48)
00175         status = self.statsSvc.put('/file1/basics/genRZ', 
00176                                       self.hist["genRZ"])
00177         if status.isFailure(): return status
00178         
00179         self.hist["genXY"] = TH2F("genXY", "Generation Vertex X-Y", 
00180                                   100, -2.48, 2.48, 100, -2.48, 2.48)
00181         status = self.statsSvc.put('/file1/basics/genXY', 
00182                                       self.hist["genXY"])
00183         if status.isFailure(): return status
00184         
00185         self.hist["HitTime"] = TH1F("HitTime", "Hit Time [#mus]",
00186                                     10000, 0.0, 10.0)
00187         status = self.statsSvc.put('/file1/basics/HitTime', 
00188                                       self.hist["HitTime"])
00189         if status.isFailure(): return status
00190 
00191         self.hist["LongHitTime"] = TH1F("LongHitTime", "Hit Time [#mus]",
00192                                         2000, 0.0, 2000.)
00193         status = self.statsSvc.put('/file1/basics/LongHitTime', 
00194                                       self.hist["LongHitTime"])
00195         if status.isFailure(): return status
00196 
00197         self.hist["GdCapHitTime"] = TH1F("GdCapHitTime",
00198                                          "Gd n-Capture Hit Time [#mus]",
00199                                          2000, 0.0, 2000)
00200         status = self.statsSvc.put('/file1/basics/GdCapHitTime', 
00201                                       self.hist["GdCapHitTime"])
00202         if status.isFailure(): return status
00203         
00204         self.hist["HCapHitTime"] = TH1F("HCapHitTime",
00205                                         "H n-Capture Hit Time [#mus]",
00206                                         2000, 0.0, 2000)
00207         status = self.statsSvc.put('/file1/basics/HCapHitTime', 
00208                                       self.hist["HCapHitTime"])
00209         if status.isFailure(): return status
00210       
00211         
00212         self.hist["hCapTarget"] = TH1F("hCapTarget",
00213                                         "n-Capture Target-Z",
00214                                         100, 0, 100)
00215         status = self.statsSvc.put('/file1/basics/hCapTarget', 
00216                                       self.hist["hCapTarget"])
00217         if status.isFailure(): return status
00218         
00219         # Add time response, Liang Zhan
00220         #################################################################
00221         self.hist["event_span_time"] = TH1F("event_span_time",
00222                                             "Hit span time for one event [ns]",
00223                                             500, 0, 500)
00224         status = self.statsSvc.put('/file1/basics/event_span_time', 
00225                                       self.hist["event_span_time"])
00226         if status.isFailure(): return status
00227         self.hist["pmt_span_time"] = TH1F("pmt_span_time",
00228                                             "Hit span time for one pmt [ns]",
00229                                             500, 0, 500)
00230         status = self.statsSvc.put('/file1/basics/pmt_span_time', 
00231                                       self.hist["pmt_span_time"])
00232         if status.isFailure(): return status
00233         self.hist["pmt_hit_time"] = TH1F("pmt_hit_time",
00234                                             "pmt first hit time [ns]",
00235                                             500, 0, 500)
00236         status = self.statsSvc.put('/file1/basics/pmt_hit_time', 
00237                                       self.hist["pmt_hit_time"])
00238         if status.isFailure(): return status
00239         #################################################################
00240         
00241         self.hist["pe_E"] = TH2F("pe_E", "PE vs visible E (e+ within GdLS)",
00242                                  100, 0, 10, 700, 0, 1400)
00243         status = self.statsSvc.put('/file1/basics/pe_E', 
00244                                       self.hist["pe_E"])
00245         if status.isFailure(): return status
00246         
00247         self.hist["pe_E_inLS"] = TH2F("pe_E_inLS",
00248                                       "PE vs visible E (e+ within LS)", 
00249                                       100, 0, 10, 700, 0, 1400)
00250         status = self.statsSvc.put('/file1/basics/pe_E_inLS', 
00251                                       self.hist["pe_E_inLS"])
00252         if status.isFailure(): return status
00253 
00254         self.hist["pe_E_LS"] = TH2F("pe_E_LS",
00255                                     "PE vs visible E (e+ in LS)",
00256                                     100, 0, 10, 700, 0, 1400)
00257         status = self.statsSvc.put('/file1/basics/pe_E_LS', 
00258                                       self.hist["pe_E_LS"])
00259         if status.isFailure(): return status
00260 
00261         self.hist["pe_E_all"] = TH2F("pe_E_all", "PE vs visible E (all e+)",
00262                                      100, 0, 10, 700, 0, 1400)
00263         status = self.statsSvc.put('/file1/basics/pe_E_all', 
00264                                       self.hist["pe_E_all"])
00265         if status.isFailure(): return status
00266         
00267         # Add response on position, Liang Zhan
00268         #############################################################
00269         self.hist["pe_yield_R2"] = TH2F("pe_yield_R2", "PE/Evis vs. R2",
00270                                         620, 0, 6.190144, 300, 0, 300) 
00271         status = self.statsSvc.put('/file1/basics/pe_yield_R2', 
00272                                       self.hist["pe_yield_R2"])
00273         if status.isFailure(): return status
00274         self.hist["pe_yield_Z"] = TH2F("pe_yield_Z", "PE/Evis vs. Z",
00275                                         200, -2.5, 2.5, 300, 0, 300) 
00276         status = self.statsSvc.put('/file1/basics/pe_yield_Z', 
00277                                       self.hist["pe_yield_Z"])
00278         if status.isFailure(): return status
00279         self.hist["pe_yield_RZ"] = TH3F("pe_yield_RZ", "PE/Evis vs. RZ",
00280                                        4, 0, 1.55, 8, -1.55, 1.55, 300, 0, 300) 
00281         status = self.statsSvc.put('/file1/basics/pe_yield_RZ', 
00282                                       self.hist["pe_yield_RZ"])
00283         if status.isFailure(): return status
00284         #############################################################
00285         
00286         self.hist["nHits_LS"] = TH1F("nHits_LS", 
00287                                      "nCap Hits in LS",
00288                                      700, 0, 1400)
00289         status = self.statsSvc.put('/file1/basics/nHits_LS', 
00290                                       self.hist["nHits_LS"])
00291         if status.isFailure(): return status
00292         
00293         self.hist["nHits_MO"] = TH1F("nHits_MO", 
00294                                      "nCap Hits in MO",
00295                                      700, 0, 1400)
00296         status = self.statsSvc.put('/file1/basics/nHits_MO', 
00297                                       self.hist["nHits_MO"])
00298         if status.isFailure(): return status
00299         
00300         self.hist["nHits_onGd"] = TH1F("nHits_onGd",
00301                                        "on-Gd nCap Hits",
00302                                        700, 0, 1400)
00303         status = self.statsSvc.put('/file1/basics/nHits_onGd', 
00304                                       self.hist["nHits_onGd"])
00305         if status.isFailure(): return status
00306         
00307         self.hist["nHits_onH"] = TH1F("nHits_onH", 
00308                                       "on-H nCap Hits", 
00309                                       700, 0, 1400)
00310         status = self.statsSvc.put('/file1/basics/nHits_onH', 
00311                                       self.hist["nHits_onH"])
00312         if status.isFailure(): return status
00313         
00314         self.hist["nHits"] = TH1F("nHits", "Neutron Capture Hits in GdLS",
00315                                   700, 0, 1400)
00316         status = self.statsSvc.put('/file1/basics/nHits', 
00317                                       self.hist["nHits"])
00318         if status.isFailure(): return status
00319         
00320         self.hist["nHits_all"] = TH1F("nHits_all", "Neutron Capture Hits",
00321                                   700, 0, 1400)
00322         status = self.statsSvc.put('/file1/basics/nHits_all', 
00323                                       self.hist["nHits_all"])
00324         if status.isFailure(): return status
00325 
00326         self.hist["nGdCapGenPos"] = TH1F("nGdCapGenPos",
00327                                          "Neutron Gd Capture Generation Position",
00328                                          620, 0, 6.190144)
00329         status = self.statsSvc.put('/file1/basics/nGdCapGenPos',
00330                                       self.hist["nGdCapGenPos"])
00331         if status.isFailure(): return status
00332         
00333         self.hist["nGdGdsCapGenPos"] = TH1F("nGdGdsCapGenPos",
00334                                             "Neutron Gd Capture in Gds Generation Position",
00335                                             620, 0, 6.190144)
00336         status = self.statsSvc.put('/file1/basics/nGdGdsCapGenPos',
00337                                       self.hist["nGdGdsCapGenPos"])
00338         if status.isFailure(): return status
00339         
00340         self.hist["nGdCapPos"] = TH3F("nGdCapPos",
00341                                       "Neutron Gd Capture Position",
00342                                       100, -2.5, 2.5,
00343                                       100, -2.5, 2.5,
00344                                       100, -2.5, 2.5)
00345         status = self.statsSvc.put('/file1/basics/nGdCapPos',
00346                                       self.hist["nGdCapPos"])
00347         if status.isFailure(): return status
00348 
00349         self.hist["nGdCapOilPos"] = TH3F("nGdCapOilPos",
00350                                          "Neutron Gd Capture beyond GdLS Position",
00351                                          100, -2.5, 2.5,
00352                                          100, -2.5, 2.5,
00353                                          100, -2.5, 2.5)
00354         status = self.statsSvc.put('/file1/basics/nGdCapOilPos',
00355                                       self.hist["nGdCapOilPos"])
00356         if status.isFailure(): return status
00357 
00358         self.hist["nGdCapOilPos_RZ"] = TH2F("nGdCapOilPos_RZ",
00359                                             "Gd Cap beyond GdLS R-Z",
00360                                             100, 0.0, 6.25,
00361                                             100, -2.5, 2.5)
00362         status = self.statsSvc.put('/file1/basics/nGdCapOilPos_RZ',
00363                                       self.hist["nGdCapOilPos_RZ"])
00364         if status.isFailure(): return status
00365 
00366         # Add Gd and H capture position, Liang Zhan.
00367         ##########################################################
00368         self.hist["nGdGdsCapPos_R"] = TH1F("nGdGdsCapPos_R",
00369                                            "Gd Cap R^{2} position in GdLS",
00370                                            620, 0, 6.190144)
00371         status = self.statsSvc.put('/file1/basics/nGdGdsCapPos_R',
00372                                       self.hist["nGdGdsCapPos_R"])
00373         if status.isFailure(): return status
00374 
00375         self.hist["nHCapPos_R"] = TH1F("nHCapPos_R",
00376                                            "H Cap R^{2} position",
00377                                            620, 0, 6.190144)
00378         status = self.statsSvc.put('/file1/basics/nHCapPos_R',
00379                                       self.hist["nHCapPos_R"])
00380         if status.isFailure(): return status
00381         ##########################################################
00382 
00383         self.hist["nHCapGenPos"] = TH1F("nHCapGenPos",
00384                                         "Neutron H Capture Generation Position",
00385                                         620, 0, 6.190144)
00386         status = self.statsSvc.put('/file1/basics/nHCapGenPos',
00387                                       self.hist["nHCapGenPos"])
00388         if status.isFailure(): return status
00389 
00390 
00391         self.hist["nHOilCapGenPos"] = TH1F("nHOilCapGenPos",
00392                                            "Neutron H Capture Generation Position",
00393                                            620, 0, 6.190144)
00394         status = self.statsSvc.put('/file1/basics/nHOilCapGenPos',
00395                                       self.hist["nHOilCapGenPos"])
00396         if status.isFailure(): return status
00397 
00398         self.hist["nHLsoCapGenPos"] = TH1F("nHLsoCapGenPos",
00399                                            "Neutron H Capture Generation Position",
00400                                            620, 0, 6.190144)
00401         status = self.statsSvc.put('/file1/basics/nHLsoCapGenPos',
00402                                       self.hist["nHLsoCapGenPos"])
00403         if status.isFailure(): return status
00404 
00405 
00406         self.hist["eDepInGdLS"] = TH1F("eDepInGdLS", "Deposited Energy [MeV]",
00407                                        200, 0, 20)
00408         status = self.statsSvc.put('/file1/basics/eDepInGdLS', 
00409                                       self.hist["eDepInGdLS"])
00410         if status.isFailure(): return status
00411 
00412         self.hist["eDepInLS"] = TH1F("eDepInLS", "Deposited Energy [MeV]",
00413                                        200, 0, 20)
00414         status = self.statsSvc.put('/file1/basics/eDepInLS', 
00415                                       self.hist["eDepInLS"])
00416         if status.isFailure(): return status
00417 
00418         self.hist["drift_GdLS"] = TH1F("drift_GdLS",
00419                                       "Neutron Drift Distance in GdLS [cm]",
00420                                        200, 0, 50)
00421         status = self.statsSvc.put('/file1/basics/drift_GdLS', 
00422                                       self.hist["drift_GdLS"])
00423         if status.isFailure(): return status
00424 
00425         self.hist["drift_LS"] = TH1F("drift_LS",
00426                                       "Neutron Drift Distance in LS [cm]",
00427                                        200, 0, 50)
00428         status = self.statsSvc.put('/file1/basics/drift_LS', 
00429                                       self.hist["drift_LS"])
00430         if status.isFailure(): return status
00431         
00432         self.hist["time_GdLS"] = TH1F("time_GdLS",
00433                                       "Neutron Capture time in GdLS [#mus]",
00434                                        400, 0, 400)
00435         status = self.statsSvc.put('/file1/basics/time_GdLS', 
00436                                       self.hist["time_GdLS"])
00437         if status.isFailure(): return status
00438 
00439         self.hist["time_LS"] = TH1F("time_LS",
00440                                       "Neutron Capture Time in LS [#mus]",
00441                                        400, 0, 2000)
00442         status = self.statsSvc.put('/file1/basics/time_LS', 
00443                                       self.hist["time_LS"])
00444         if status.isFailure(): return status
00445 
00446         #sub-volume eff
00447         self.hist["countR"] = TH1F("countR",
00448                                    "counts in each sub-R",
00449                                    16,0, 2.4);
00450         status = self.statsSvc.put('/file1/basics/countR', 
00451                                    self.hist["countR"])
00452         if status.isFailure(): return status
00453         
00454         self.hist["countCutR"] = TH1F("countCutR",
00455                                       "counts above cut in each sub-R",
00456                                       16, 0, 2.4)
00457         status = self.statsSvc.put('/file1/basics/countCutR', 
00458                                    self.hist["countCutR"])
00459         if status.isFailure(): return status
00460         
00461         self.hist["countZ"] = TH1F("countZ",
00462                                    "counts in each sub-Z",
00463                                    16, -1.55, 1.55)
00464         status = self.statsSvc.put('/file1/basics/countZ', 
00465                                    self.hist["countZ"])
00466         if status.isFailure(): return status
00467         
00468         self.hist["countCutZ"] = TH1F("countCutZ",
00469                                       "counts above cut in each sub-Z",
00470                                       16, -1.55, 1.55)
00471         status = self.statsSvc.put('/file1/basics/countCutZ', 
00472                                    self.hist["countCutZ"])
00473         if status.isFailure(): return status
00474         
00475         self.hist["count"] = TH2F("count",
00476                                   "counts in each sub-volume",
00477                                   10, 0, 2.4, 10, -1.55, 1.55 )
00478         status = self.statsSvc.put('/file1/basics/count', 
00479                                    self.hist["count"])
00480         if status.isFailure(): return status
00481         
00482         self.hist["countCut"] = TH2F("countCut",
00483                                      "counts above cut in each sub-volume",
00484                                      10, 0, 2.4, 10, -1.55, 1.55 )
00485         status = self.statsSvc.put('/file1/basics/countCut', 
00486                                    self.hist["countCut"])
00487         if status.isFailure(): return status
00488         
00489         self.hist["countPECut"] = TH2F("countPECut",
00490                                        "counts above PE cut in each sub-volume",
00491                                        10, 0, 2.4, 10, -1.55, 1.55 )
00492         status = self.statsSvc.put('/file1/basics/countPECut', 
00493                                    self.hist["countPECut"])
00494         if status.isFailure(): return status
00495         
00496         self.hist["eff"] = TH2F("eff", "efficiency in each sub-volume",
00497                                 10, 0, 2.4, 10, -1.55, 1.55 )
00498         status = self.statsSvc.put('/file1/basics/eff', 
00499                                    self.hist["eff"])
00500         if status.isFailure(): return status        
00501 
00502         self.hist["effR"] = TH1F("effR", "efficiency in each sub-volume",
00503                                  16, 0, 2.4)
00504         status = self.statsSvc.put('/file1/basics/effR', 
00505                                    self.hist["effR"])
00506         if status.isFailure(): return status
00507         
00508         self.hist["effZ"] = TH1F("effZ", "efficiency in each sub-volume",
00509                                  16, -1.55, 1.55)
00510         
00511         status = self.statsSvc.put('/file1/basics/effZ', 
00512                                    self.hist["effZ"])
00513         if status.isFailure(): return status        
00514 
00515         #sub-volume eff H
00516         self.hist["HcountR"] = TH1F("HcountR",
00517                                     "counts in each sub-R",
00518                                     26, 0, 3.92);
00519         status = self.statsSvc.put('/file1/basics/HcountR', 
00520                                    self.hist["HcountR"])
00521         if status.isFailure(): return status
00522         
00523         self.hist["HcountCutR"] = TH1F("HcountCutR",
00524                                        "counts above cut in each sub-R",
00525                                        26, 0, 3.92)
00526         status = self.statsSvc.put('/file1/basics/HcountCutR', 
00527                                    self.hist["HcountCutR"])
00528         if status.isFailure(): return status
00529         
00530         self.hist["HcountZ"] = TH1F("HcountZ",
00531                                     "counts in each sub-Z",
00532                                     26, -1.98, 1.98)
00533         status = self.statsSvc.put('/file1/basics/HcountZ', 
00534                                    self.hist["HcountZ"])
00535         if status.isFailure(): return status
00536         
00537         self.hist["HcountCutZ"] = TH1F("HcountCutZ",
00538                                       "counts above cut in each sub-Z",
00539                                       26, -1.98, 1.98)
00540         status = self.statsSvc.put('/file1/basics/HcountCutZ', 
00541                                    self.hist["HcountCutZ"])
00542         if status.isFailure(): return status
00543         
00544         self.hist["Hcount"] = TH2F("Hcount",
00545                                    "counts in each sub-volume",
00546                                    12, 0, 3.92, 12, -1.98, 1.98)
00547         status = self.statsSvc.put('/file1/basics/Hcount', 
00548                                    self.hist["Hcount"])
00549         if status.isFailure(): return status
00550         
00551         self.hist["HcountCut"] = TH2F("HcountCut",
00552                                       "counts above cut in each sub-volume",
00553                                       12, 0, 3.92, 12, -1.98, 1.98)
00554         status = self.statsSvc.put('/file1/basics/HcountCut', 
00555                                    self.hist["HcountCut"])
00556         if status.isFailure(): return status
00557         
00558         self.hist["HcountPECut"] = TH2F("HcountPECut",
00559                                         "counts above PE cut in each sub-volume",
00560                                         12, 0, 3.92, 12, -1.98, 1.98)
00561         status = self.statsSvc.put('/file1/basics/HcountPECut', 
00562                                    self.hist["HcountPECut"])
00563         if status.isFailure(): return status
00564         
00565         return SUCCESS
00566 
00567     def execute(self):
00568         print "Executing plotIbdBasics", self.name()
00569         evt = self.evtSvc()
00570         simhdr = evt['/Event/Sim/SimHeader']
00571 #        print "SimHeader: ", simhdr
00572         if simhdr == None:
00573             print "No SimHeader in this ReadOut. Skip."
00574             return SUCCESS
00575 
00576 #        det = self.detSvc(self.target_de_name)
00577 #        det_gds = self.detSvc(self.gds_de_name)
00578 #        det_lso = self.detSvc(self.lso_de_name)
00579 
00580         # Unobservables
00581         statshdr = simhdr.unobservableStatistics()
00582         stats = statshdr.stats()
00583 
00584         PID_trk1 = stats["pdgId_Trk1"].sum()
00585         PID_trk2 = stats["pdgId_Trk2"].sum()
00586 
00587         if PID_trk1 != -11 or PID_trk2 != 2112:
00588             print "PID of track 1 is", PID_trk1
00589             print "PID of track 2 is", PID_trk2
00590             print "Not an IBD event."
00591             return SUCCESS
00592 
00593         tGen = stats["t_Trk2"].sum()
00594         xGen = stats["x_Trk2"].sum()
00595         yGen = stats["y_Trk2"].sum()
00596         zGen = stats["z_Trk2"].sum()
00597         
00598         tCap = stats["tEnd_Trk2"].sum()
00599         xCap = stats["xEnd_Trk2"].sum()
00600         yCap = stats["yEnd_Trk2"].sum()
00601         zCap = stats["zEnd_Trk2"].sum()
00602         
00603         # Get underlying DE object
00604         de = self.getDet(self.target_de_name)
00605         if not de:
00606             print 'Failed to get DE',self.target_de_name
00607             return FAILURE
00608         
00609 #        de_lso = self.getDet(self.lso_de_name)
00610 #        de_gds = self.getDet(self.gds_de_name)
00611 #        if not de_lso:
00612 #            print 'Failed to get DE',self.lso_de_name
00613 #            return FAILURE        
00614 #        if not de_gds:
00615 #            print 'Failed to get DE',self.gds_de_name
00616 #            return FAILURE
00617         
00618         # Get the AD coordinates of the vertexes
00619         Gaudi = PyCintex.makeNamespace('Gaudi')
00620         genGlbPoint = Gaudi.XYZPoint(xGen, yGen, zGen)
00621         capGlbPoint = Gaudi.XYZPoint(xCap, yCap, zCap)
00622 #        point = de.geometry().toGlobal(point)
00623         genLclPoint = de.geometry().toLocal(genGlbPoint)
00624         capLclPoint = de.geometry().toLocal(capGlbPoint)
00625 #        genLclPointLso = de_lso.geometry().toLocal(genGlbPoint)
00626 #        capLclPointLso = de_lso.geometry().toLocal(capGlbPoint)
00627 #        genLclPointGds = de_gds.geometry().toLocal(genGlbPoint)
00628 #        capLclPointGds = de_gds.geometry().toLocal(capGlbPoint)
00629 #        print 'In global coordinate [',gpoint.x(),gpoint.y(),gpoint.z(),']'
00630 
00631         ndrift = ROOT.TVector3(xCap-xGen, yCap-yGen, zCap-zGen)
00632         
00633         capTime = tCap - tGen
00634         capDis = ndrift.Mag()
00635 
00636         R2 = genLclPoint.x()/units.meter * genLclPoint.x()/units.meter + \
00637              genLclPoint.y()/units.meter * genLclPoint.y()/units.meter
00638         
00639         R2Cap = capLclPoint.x()/units.meter * capLclPoint.x()/units.meter + \
00640                 capLclPoint.y()/units.meter * capLclPoint.y()/units.meter
00641 
00642         self.hist["genRZ"].Fill(R2, genLclPoint.z()/units.meter)
00643 
00644         self.hist["genXY"].Fill(genLclPoint.x()/units.meter,genLclPoint.y()/units.meter)
00645 
00646         # Find the interesting volumes
00647         genDE = self.coorSvc.coordSysDE(genGlbPoint)
00648         capDE = self.coorSvc.coordSysDE(capGlbPoint)
00649         if not genDE:
00650             print 'Failed to find coordinate system DE for generation'\
00651                 '[',genGlbPoint.x(),genGlbPoint.y(),genGlbPoint.z(),']'
00652             print 'Local: [',genLclPoint.x(),genLclPoint.y(),genLclPoint.z(),']'
00653             return FAILURE
00654         else:
00655             self.counter["nGen"] += 1
00656             gendmvol = genDE.geometry().belongsToPath(genGlbPoint,-1)
00657 
00658         if not capDE:
00659             print 'Failed to find coordinate system DE for capture'\
00660                 '[',capGlbPoint.x(),capGlbPoint.y(),capGlbPoint.z(),']'
00661             return FAILURE
00662         else:
00663             self.counter["nCap"] += 1
00664             capdmvol = capDE.geometry().belongsToPath(capGlbPoint,-1)
00665 
00666         genDM = re.split('/', gendmvol).pop()
00667         capDM = re.split('/', capdmvol).pop()
00668         print "Generated in ", genDM
00669         print "Captured in ", capDM
00670         
00671         positronHits = 0.
00672         neutronHits = 0.
00673         positronTimeCut = 500.
00674 
00675         # Add coincidence time cut, Liang Zhan
00676         ################################################
00677         timeCutL = 1 #mus
00678         timeCutH = 200 #mus
00679         timeCutUncertain = 0.01 #mus, 10ns
00680         ################################################
00681  
00682         positronCut = 126.2  # 1 MeV
00683         positronCutL = 123.7 # 1 MeV
00684         positronCutR = 128.7 # 1 MeV
00685 
00686 #        onHLowerCut = 194.   # 1.5 MeV
00687 #        onHLowerCutL = 190.1 # 194 (1 - 2%)
00688 #        onHLowerCutR = 197.9 # 194 (1 + 2%)
00689                 
00690 #        onHUpperCut = 453.  # 3.5 MeV
00691 #        onHUpperCutL = 443.9 # 453 (1 - 2%)
00692 #        onHUpperCutR = 462.1 # 453 (1 + 2%)
00693 
00694         onHLowerCut = 175.   # ? MeV
00695         onHLowerCutL = 171.5 # x (1 - 2%)
00696         onHLowerCutR = 178.5 # x (1 + 2%)
00697 
00698         onHUpperCut = 481.4  # 3.5 MeV
00699         onHUpperCutL = 471.8 # x (1 - 2%)
00700         onHUpperCutR = 491.0 # x (1 + 2%)
00701 
00702         onGdCut = 811.  # 6MeV
00703         onGdCutL = 802.9 # 811 (1 - 1%)
00704         onGdCutR = 819.1 # 811 (1 + 1%)
00705 
00706         # Capture target
00707         capTarget = stats["capTarget"].sum()
00708         print "The capture target is ", capTarget
00709 
00710         # Deposit energy
00711         eDepInGdLS = stats["EDepInGdLS"].sum()
00712         self.hist["eDepInGdLS"].Fill(eDepInGdLS/units.MeV)
00713         eDepInLS = stats["EDepInLS"].sum()
00714         self.hist["eDepInLS"].Fill(eDepInLS/units.MeV)
00715 
00716         # For PE hit time, Liang Zhan
00717         firstHit = {} #dictionary for PMT first hit, key is PMT id
00718         lastHit = {} #dictionary for PMT last hit, key is PMT id
00719         hitMap = {} #pmt hit map
00720         firstT = 501 # time of first hit in one event 
00721         lastT = -1 # time of last hit in one event
00722 
00723         simhits = simhdr.hits()
00724         for detector in simhits.hitDetectors():
00725             hitCollection = simhits.hitsByDetector(detector)
00726             if hitCollection == None:
00727                 print "No hits in ", detector
00728             hits = hitCollection.collection()
00729             for hit in hits:
00730 #                print " PMT", hit.sensDetId(), "hit @ time: ", hit.hitTime()
00731                 pmtId = hit.sensDetId() 
00732                 self.hist["HitTime"].Fill(hit.hitTime()/units.microsecond)
00733                 self.hist["LongHitTime"].Fill(hit.hitTime()/units.microsecond)
00734                 if hit.hitTime()/units.nanosecond<positronTimeCut and \
00735                         hit.hitTime()/units.nanosecond<capTime/units.nanosecond:
00736                     positronHits += 1
00737                     hitTime = hit.hitTime()/units.nanosecond
00738                     if hitMap.has_key(pmtId):
00739                         hitMap[pmtId] += 1
00740                     else:
00741                         hitMap[pmtId] = 1.
00742                     if firstT > hitTime:
00743                         firstT = hitTime
00744                     if lastT < hitTime:
00745                         lastT = hitTime
00746                     if not firstHit.has_key(pmtId):
00747                         firstHit[pmtId] = hitTime                     
00748                         lastHit[pmtId] = hitTime                     
00749                     else:
00750                         if firstHit[pmtId] > hitTime:
00751                             firstHit[pmtId] = hitTime   
00752                         if lastHit[pmtId] < hitTime:
00753                             lastHit[pmtId] = hitTime   
00754                 else:
00755                     neutronHits += 1.
00756                     if capTarget == 1:
00757                         self.hist["HCapHitTime"].Fill(hit.hitTime()/units.microsecond)
00758                     if capTarget == 64:
00759                         self.hist["GdCapHitTime"].Fill(hit.hitTime()/units.microsecond)
00760 
00761         # Fill time
00762         if positronHits > 1.5:  # 2 hits at least
00763             self.hist["event_span_time"].Fill(lastT-firstT)
00764         for id in firstHit.keys():
00765             self.hist["pmt_hit_time"].Fill(firstHit[id]-firstT)
00766             if hitMap[id] > 1.5:  # 2 hits at least
00767                 self.hist["pmt_span_time"].Fill(lastHit[id]-firstHit[id])
00768 
00769         # Visible energy
00770         vis_trk1 = 1.022 + stats['ke_Trk1'].sum()/units.MeV
00771         self.hist["pe_E_all"].Fill(vis_trk1, positronHits)
00772 
00773         if genLclPoint.z()/units.m < 1. and genLclPoint.z()/units.m > -1.:
00774             self.hist["pe_yield_R2"].Fill(R2, positronHits/vis_trk1)
00775         if R2 < 1.:
00776             self.hist["pe_yield_Z"].Fill(genLclPoint.z()/units.m, positronHits/vis_trk1)
00777 
00778         self.hist["nHits_all"].Fill(neutronHits)
00779 
00780         if genDM == 'db-gds1':
00781             self.counter["nInGdLS"] += 1
00782             self.hist["pe_E"].Fill(vis_trk1, positronHits)
00783             self.hist["pe_yield_RZ"].Fill(ROOT.sqrt(R2), genLclPoint.z()/units.m, positronHits/vis_trk1)
00784             if capDM == 'db-gds1':
00785                 self.hist["nHits"].Fill(neutronHits)
00786                 self.hist["drift_GdLS"].Fill(capDis/units.cm)
00787                 self.hist["time_GdLS"].Fill(capTime/units.microsecond)
00788                 self.counter["nGdsCapTime"] += 1
00789                 if(capTime/units.microsecond > timeCutL and capTime/units.microsecond < timeCutH):
00790                     self.counter["nGdsCapTimePassCut"] += 1
00791                 if(capTime/units.microsecond < timeCutL):
00792                     self.counter["nGdsCapTimeL"] += 1
00793                 if(capTime/units.microsecond > timeCutH):
00794                     self.counter["nGdsCapTimeH"] += 1
00795                 if(capTime/units.microsecond > timeCutL-timeCutUncertain/2 and capTime/units.microsecond < timeCutL+timeCutUncertain/2):
00796                     self.counter["nGdsCapTimeUncertainL"] += 1
00797                 if(capTime/units.microsecond > timeCutH-timeCutUncertain/2 and capTime/units.microsecond < timeCutH+timeCutUncertain/2):
00798                     self.counter["nGdsCapTimeUncertainH"] += 1
00799             else:
00800                 self.counter["SpillOut"] += 1
00801         
00802         if genDM != 'db-gds1':
00803             if capTarget == 64:
00804                 self.counter["SpillIn_onGd"] += 1
00805                 
00806         if genDM == 'db-iav1':
00807             self.counter["nInIAV"] += 1
00808             if capTarget == 64:
00809                 self.counter["IavSpillIn_onGd"] += 1
00810 
00811         if genDM == 'db-lso1':
00812             self.counter["nInLSO"] += 1
00813             self.hist["pe_E_LS"].Fill(vis_trk1, positronHits)
00814             if capDM == 'db-lso1':
00815                 self.hist["drift_LS"].Fill(capDis/units.cm)
00816                 self.hist["time_LS"].Fill(capTime/units.microsecond)
00817                 self.counter["nLSCapTime"] += 1
00818                 if(capTime/units.microsecond > timeCutL and capTime/units.microsecond < timeCutH):
00819                     self.counter["nLSCapTimePassCut"] += 1
00820                 if(capTime/units.microsecond < timeCutL):
00821                     self.counter["nLSCapTimeL"] += 1
00822                 if(capTime/units.microsecond > timeCutH):
00823                     self.counter["nLSCapTimeH"] += 1
00824                 if(capTime/units.microsecond > timeCutL-timeCutUncertain/2 and capTime/units.microsecond < timeCutL+timeCutUncertain/2):
00825                     self.counter["nLSCapTimeUncertainL"] += 1
00826                 if(capTime/units.microsecond > timeCutH-timeCutUncertain/2 and capTime/units.microsecond < timeCutH+timeCutUncertain/2):
00827                     self.counter["nLSCapTimeUncertainH"] += 1
00828             if capDM == 'db-oil1' or capDM == 'db-oav1':
00829                 self.counter["SpillOut_onH"] += 1
00830             if capTarget == 64:
00831                 self.counter["LSSpillIn_onGd"] += 1
00832 
00833         if genDM == 'db-oav1':
00834             self.counter["nInOAV"] += 1
00835             if re.search('db-lso1', capdmvol) and capTarget == 1:
00836                 self.counter["OavSpillIn_onH"] += 1
00837                 
00838                 if neutronHits > onHLowerCut and neutronHits < onHUpperCut:
00839                     self.counter["onHCapOavMoPassCut"] += 1
00840                 if neutronHits > onHLowerCutL and neutronHits < onHLowerCutR:
00841                     self.counter["onHCapOavMoCutUncertainL"] += 1
00842                 if neutronHits > onHUpperCutL and neutronHits < onHUpperCutR:
00843                     self.counter["onHCapOavMoCutUncertainU"] += 1
00844                     
00845         if genDM == 'db-oil1':
00846             self.counter["nInMO"] += 1
00847             if re.search('db-lso1', capdmvol) and capTarget == 1:
00848                 self.counter["OilSpillIn_onH"] += 1
00849                 if neutronHits > onHLowerCut and neutronHits < onHUpperCut:
00850                     self.counter["onHCapOavMoPassCut"] += 1
00851                 if neutronHits > onHLowerCutL and neutronHits < onHLowerCutR:
00852                     self.counter["onHCapOavMoCutUncertainL"] += 1
00853                 if neutronHits > onHUpperCutL and neutronHits < onHUpperCutR:
00854                     self.counter["onHCapOavMoCutUncertainU"] += 1
00855 
00856         if re.search('db-lso1',gendmvol):
00857             self.hist["pe_E_inLS"].Fill(vis_trk1, positronHits)
00858 
00859         if re.search('db-lso1',capdmvol):
00860             if capTarget == 1:
00861                 self.counter["onHCapWithinLS"] += 1
00862                 if neutronHits > onHLowerCut and neutronHits < onHUpperCut:
00863                     self.counter["onHCapWithinLSPassCut"] += 1
00864                 if neutronHits > onHLowerCutL and neutronHits < onHLowerCutR:
00865                     self.counter["onHCapWithinLSCutUncertainL"] += 1
00866                 if neutronHits > onHUpperCutL and neutronHits < onHUpperCutR:
00867                     self.counter["onHCapWithinLSCutUncertainU"] += 1
00868             else:
00869                 if neutronHits > onHLowerCut and neutronHits < onHUpperCut:
00870                     self.counter["nonHCapWithinLSPassHCut"] += 1
00871                 if neutronHits > onHLowerCutL and neutronHits < onHLowerCutR:
00872                     self.counter["nonHCapWithinLSHCutUncertainL"] += 1
00873                 if neutronHits > onHUpperCutL and neutronHits < onHUpperCutR:
00874                     self.counter["nonHCapWithinLSHCutUncertainU"] += 1
00875                 
00876         else:
00877             if capTarget == 1:
00878                 self.counter["onHCapBeyondLS"] += 1
00879                 if neutronHits > onHLowerCut and neutronHits < onHUpperCut:
00880                     self.counter["onHCapBeyondLSPassCut"] += 1
00881                 if neutronHits > onHLowerCutL and neutronHits < onHLowerCutR:
00882                     self.counter["onHCapBeyondLSCutUncertainL"] += 1
00883                 if neutronHits > onHUpperCutL and neutronHits < onHUpperCutR:
00884                     self.counter["onHCapBeyondLSCutUncertainU"] += 1
00885             else:
00886                 if neutronHits > onHLowerCut and neutronHits < onHUpperCut:
00887                     self.counter["nonHCapBeyondLSPassHCut"] += 1
00888                 if neutronHits > onHLowerCutL and neutronHits < onHLowerCutR:
00889                     self.counter["nonHCapBeyondLSHCutUncertainL"] += 1
00890                 if neutronHits > onHUpperCutL and neutronHits < onHUpperCutR:
00891                     self.counter["nonHCapBeyondLSHCutUncertainU"] += 1
00892             
00893         # Passing on-Gd neutron cuts
00894         if neutronHits > onGdCut:
00895             self.counter["nPassGdCut"] += 1
00896             if capTarget != 64:
00897                 print "non-Gd capture pass Gd cut: ", capTarget
00898                 self.counter["nonGdPassGdCut"] += 1
00899                 value = ('Non-Gd capTarget: ', capTarget)
00900                 theline = str(value)
00901                 self.outputstat.write(theline)
00902                 self.outputstat.write('\n')
00903                 
00904         if neutronHits > onGdCutL and neutronHits < onGdCutR:
00905             self.counter["nGdCutUncertain"] += 1
00906 
00907         # Passing on-H neutron cuts
00908         if neutronHits > onHLowerCut and neutronHits < onHUpperCut:
00909             self.counter["nPassHCut"] += 1
00910             if capTarget == 64:
00911                 self.counter["nGdPassHCut"] += 1
00912         if neutronHits > onHLowerCutL and neutronHits < onHLowerCutR:
00913             self.counter["nHCutUncertainL"] += 1
00914             if capTarget == 64:
00915                 self.counter["nGdHCutUncertainL"] += 1                
00916         if neutronHits > onHUpperCutL and neutronHits < onHUpperCutR:
00917             self.counter["nHCutUncertainU"] += 1
00918             if capTarget == 64:
00919                 self.counter["nGdHCutUncertainU"] += 1
00920         
00921         if capTarget == 64:
00922             self.counter["onGdCap"] += 1
00923             self.hist["nHits_onGd"].Fill(neutronHits)
00924             self.hist["nGdCapPos"].Fill(capLclPoint.x()/units.meter,capLclPoint.y()/units.meter,capLclPoint.z()/units.meter)
00925             if capDM == 'db-gds1':
00926                 self.counter["onGdGdsCap"] += 1
00927                 self.hist["count"].Fill(R2Cap, capLclPoint.z()/units.meter)
00928                 if (eDepInLS/units.MeV + eDepInGdLS/units.MeV) > 6.:
00929                     self.hist["countCut"].Fill(R2Cap, capLclPoint.z()/units.meter)
00930                     self.counter["onGdGdsCapPassEdepCut"] += 1
00931                 
00932                 if capLclPoint.z()/units.m < 1. and capLclPoint.z()/units.m > -1.:
00933                     self.hist["nGdGdsCapPos_R"].Fill(R2Cap)
00934                     self.hist["countR"].Fill(R2Cap)
00935                     if (eDepInLS/units.MeV + eDepInGdLS/units.MeV) > 6.:
00936                         self.hist["countCutR"].Fill(R2Cap)
00937                     
00938                 if R2Cap < 0.64:
00939                     self.hist["countZ"].Fill(capLclPoint.z()/units.m)
00940                     if (eDepInLS/units.MeV + eDepInGdLS/units.MeV) > 6.:
00941                         self.hist["countCutZ"].Fill(capLclPoint.z()/units.m)
00942                         
00943                 if neutronHits > onGdCut:
00944                     self.counter["onGdGdsCapPassCut"] += 1
00945                     self.hist["countPECut"].Fill(R2Cap, capLclPoint.z()/units.meter)
00946                 if neutronHits > onGdCutL and neutronHits < onGdCutR:
00947                     self.counter["onGdGdsCapCutUncertain"] += 1
00948                 if positronHits > positronCut:
00949                     self.counter["positronPassCut"] += 1
00950                 if positronHits > positronCutL and positronHits < positronCutR:
00951                     self.counter["positronCutUncertain"] += 1
00952             else:
00953                 self.hist["nGdCapOilPos"].Fill(capLclPoint.x()/units.meter,capLclPoint.y()/units.meter,capLclPoint.z()/units.meter)
00954                 self.hist["nGdCapOilPos_RZ"].Fill(R2Cap, capLclPoint.z()/units.meter)
00955 
00956             if genLclPoint.z()/units.m < 1. and genLclPoint.z()/units.m > -1.:
00957                 self.hist["nGdCapGenPos"].Fill(R2)
00958                 if capDM == 'db-gds1':
00959                     self.hist["nGdGdsCapGenPos"].Fill(R2)
00960                     
00961             if neutronHits > onGdCut:
00962                 self.counter["onGdCapPassCut"] += 1
00963             if neutronHits > onGdCutL and neutronHits < onGdCutR:
00964                 self.counter["onGdCapCutUncertain"] += 1
00965             
00966         if capTarget == 1:
00967             self.counter["onHCap"] += 1
00968             self.hist["nHits_onH"].Fill(neutronHits)
00969             
00970             if capLclPoint.z()/units.m < 1. and capLclPoint.z()/units.m > -1.:
00971                 self.hist["nHCapPos_R"].Fill(R2Cap)
00972 
00973             if genLclPoint.z()/units.m < 1. and genLclPoint.z()/units.m > -1.:
00974                 self.hist["nHCapGenPos"].Fill(R2)
00975                 if capDM == 'db-lso1':
00976                     self.hist["nHLsoCapGenPos"].Fill(R2)
00977                 if capDM == 'db-oil1' or capDM == 'db-oav1':
00978                     self.hist["nHOilCapGenPos"].Fill(R2)
00979                     
00980             if neutronHits > onHLowerCut and neutronHits < onHUpperCut:
00981                 self.counter["onHCapPassCut"] += 1
00982             if neutronHits > onHLowerCutL and neutronHits < onHLowerCutR:
00983                 self.counter["onHCapCutUncertainL"] += 1
00984             if neutronHits > onHUpperCutL and neutronHits < onHUpperCutR:
00985                 self.counter["onHCapCutUncertainU"] += 1
00986 
00987             if re.search('db-lso1', capdmvol):
00988                 self.hist["Hcount"].Fill(R2Cap, capLclPoint.z()/units.meter)
00989                 if neutronHits > onHLowerCut and neutronHits < onHUpperCut:
00990                     self.hist["HcountCut"].Fill(R2Cap, capLclPoint.z()/units.meter)
00991                     
00992                 if capLclPoint.z()/units.m < 1. and capLclPoint.z()/units.m > -1.:
00993                     self.hist["HcountR"].Fill(R2Cap)
00994                     if neutronHits > onHLowerCut and neutronHits < onHUpperCut:
00995                         self.hist["HcountCutR"].Fill(R2Cap)
00996                     
00997                 if R2Cap < 0.64:
00998                     self.hist["HcountZ"].Fill(capLclPoint.z()/units.m)
00999                     if neutronHits > onHLowerCut and neutronHits < onHUpperCut:
01000                         self.hist["HcountCutZ"].Fill(capLclPoint.z()/units.m)
01001                         
01002         if capDM == 'db-gds1':
01003             self.hist["hCapTarget"].Fill(capTarget)
01004             if capTarget == 1:
01005                 self.counter["onHCapGds"] += 1
01006             elif capTarget == 6:
01007                 self.counter["onCCapGds"] += 1
01008             elif capTarget == 64:
01009                 pass
01010             else:
01011                 self.counter["onOtherCapGds"] += 1
01012 
01013         if capDM == 'db-lso1':
01014             self.hist["nHits_LS"].Fill(neutronHits)
01015             if capTarget == 1:
01016                 self.counter["onHCapLS"] += 1
01017                 
01018         if capDM == 'db-oav1':
01019             if capTarget == 1:
01020                 self.counter["onHCapOAV"] += 1
01021                 
01022         if capDM == 'db-oil1':
01023             self.hist["nHits_MO"].Fill(neutronHits)
01024             if capTarget == 1:
01025                 self.counter["onHCapMO"] += 1
01026 
01027         return SUCCESS
01028     
01029     def finalize(self):
01030         
01031         self.counter["GdCapDetEff"] = self.counter["onGdGdsCapPassCut"]/self.counter["onGdGdsCap"]
01032         self.counter["GdCapDetEffUncertain"] = self.counter["onGdGdsCapCutUncertain"]/self.counter["onGdGdsCap"]/2.
01033         
01034         self.counter["HCapDetEff"] = self.counter["nPassHCut"]/self.counter["onHCapWithinLS"]
01035         self.counter["HCapDetEffUncertainL"] = self.counter["nHCutUncertainL"]/self.counter["onHCapWithinLS"]/2.
01036         self.counter["HCapDetEffUncertainU"] = self.counter["nHCutUncertainU"]/self.counter["onHCapWithinLS"]/2.
01037         
01038         self.counter["HCapCanoDetEff"] = self.counter["onHCapWithinLSPassCut"]/self.counter["onHCapWithinLS"]
01039         self.counter["HCapCanoDetEffUncertainL"] = self.counter["onHCapWithinLSCutUncertainL"]/self.counter["onHCapWithinLS"]/2.
01040         self.counter["HCapCanoDetEffUncertainU"] = self.counter["onHCapWithinLSCutUncertainU"]/self.counter["onHCapWithinLS"]/2.
01041         
01042         self.counter["HCapDetEffContaminationUncertainL"] = (self.counter["nonHCapWithinLSHCutUncertainL"]+self.counter["nonHCapBeyondLSHCutUncertainL"])/self.counter["onHCapWithinLS"]/2.
01043         self.counter["HCapDetEffContaminationUncertainU"] = (self.counter["nonHCapWithinLSHCutUncertainU"]+self.counter["nonHCapBeyondLSHCutUncertainU"])/self.counter["onHCapWithinLS"]/2.
01044         
01045         self.counter["HCapDetContamination"] = (self.counter["nonHCapWithinLSPassHCut"]+self.counter["nonHCapBeyondLSPassHCut"])/self.counter["nPassHCut"]
01046         
01047         self.outputstat.close()
01048 
01049         # register all the counters as TParameters
01050         for name, counter in self.counter.iteritems():
01051             temp = TParameter('double')(name, counter)
01052             status = self.statsSvc.put("/file1/counters/"+name, temp)
01053             if status.isFailure(): return status
01054             print "Final value of "+name+": ", temp.GetVal()
01055 
01056         print ""
01057         print "Total n generation: ", self.counter["nGen"]
01058         print "Total n gen in GdLS: ", self.counter["nInGdLS"]
01059         print "Total n gen in IAV: ", self.counter["nInIAV"]
01060         print "Total n gen in LSO: ", self.counter["nInLSO"]
01061         print "Total n gen in OAV: ", self.counter["nInOAV"]
01062         print "Total n gen in MO: ", self.counter["nInMO"]
01063         print ""
01064 
01065         print "Total n capture: ", self.counter["nCap"]
01066         print "Total n-capture on Gd: ", self.counter["onGdCap"]
01067         print "Total n-capture on Gd in Gds: ", self.counter["onGdGdsCap"]
01068         print "Total n-H capture in GdLS: ", self.counter["onHCapGds"]
01069         print "Total n-C capture in GdLS: ", self.counter["onCCapGds"]
01070         print "Total other capture in GdLS: ", self.counter["onOtherCapGds"]
01071         print ""
01072         
01073         print "Spill-In on Gd n from LS: ", self.counter["LSSpillIn_onGd"]
01074         print "Spill-In on Gd from IAV: ", self.counter["IavSpillIn_onGd"]
01075         print "Spill-Out of n from Gd-LS: ", self.counter["SpillOut"]
01076         print ""
01077         
01078         print "On Gd capture pass on-Gd cut: ", self.counter["onGdCapPassCut"]
01079         print "On Gd capture on-Gd cut uncertainty: ", self.counter["onGdCapCutUncertain"]
01080         print "On Gd capture in GDS pass on-Gd cut: ", self.counter["onGdGdsCapPassCut"]
01081         print "On Gd capture in GDS pass on-Gd cut(true deposit energy): ", self.counter["onGdGdsCapPassEdepCut"]
01082         print "On Gd capture in GDS on-Gd cut uncertainty: ", self.counter["onGdGdsCapCutUncertain"]
01083         print "Total non-Gd capture pass Gd cut: ", self.counter["nonGdPassGdCut"]
01084         print "Total pass Gd cut: ", self.counter["nPassGdCut"]
01085         print "Overall Gd cut uncertainty: ", self.counter["nGdCutUncertain"]
01086         print ""
01087 
01088         print "positron pass cut for Gd capture in GDS: ", self.counter["positronPassCut"]
01089         print "positron pass cut uncertainty for Gd capture in GDS: ", self.counter["positronCutUncertain"]
01090         print ""
01091 
01092         print "on Gd Detection Efficiency: ", self.counter["GdCapDetEff"]
01093         print "on Gd Detection Efficiency Uncertainty: ", self.counter["GdCapDetEffUncertain"]
01094         print ""
01095 
01096         print "Number of capture time in GdLS: ", self.counter["nGdsCapTime"]
01097         print "Number of capture time pass cut (1-200mus) in GdLS: ", self.counter["nGdsCapTimePassCut"]
01098         print "Number of capture time less than 1mus in GdLS: ", self.counter["nGdsCapTimeL"]
01099         print "Uncertain number of capture time less than 1mus in GdLS: ", self.counter["nGdsCapTimeUncertainL"]
01100         print "Number of capture time larger than 200mus in GdLS: ", self.counter["nGdsCapTimeH"]
01101         print "Uncertain number of capture time larger than 200mus in GdLS: ", self.counter["nGdsCapTimeUncertainH"]
01102         print ""
01103 
01104         print "Number of capture time in LS: ", self.counter["nLSCapTime"]
01105         print "Number of capture time pass cut (1-200mus) in LS: ", self.counter["nLSCapTimePassCut"]
01106         print "Number of capture time less than 1mus in LS: ", self.counter["nLSCapTimeL"]
01107         print "Uncertain number of capture time less than 1mus in LS: ", self.counter["nLSCapTimeUncertainL"]
01108         print "Number of capture time larger than 200mus in LS: ", self.counter["nLSCapTimeH"]
01109         print "Uncertain number of capture time larger than 200mus in LS: ", self.counter["nLSCapTimeUncertainH"]
01110         print ""
01111         
01112         print "Total n-capture on H: ", self.counter["onHCap"]
01113         print "Total n-capture on H in LS: ", self.counter["onHCapLS"]
01114         print "Total n-capture on H within LS: ", self.counter["onHCapWithinLS"]
01115         print ""
01116         
01117         print "Spill-In of n from OIL to LSO: ", self.counter["OilSpillIn_onH"]
01118         print "Spill-In of n from OAV to LSO: ", self.counter["OavSpillIn_onH"]
01119         print "Spill-Out of n from LSO to OAV or OIL: ", self.counter["SpillOut_onH"]
01120         print ""
01121         
01122         print "n capture pass on-H cut: ", self.counter["nPassHCut"]
01123         print "n capture pass lower on-H cut uncertainty: ", self.counter["nHCutUncertainL"]
01124         print "n capture pass upper on-H cut uncertainty: ", self.counter["nHCutUncertainU"]
01125         print ""
01126 
01127         print "On H capture pass on-H cut: ", self.counter["onHCapPassCut"]
01128         print "On H capture lower on-H cut uncertainty: ", self.counter["onHCapCutUncertainL"]
01129         print "On H capture upper on-H cut uncertainty: ", self.counter["onHCapCutUncertainU"]
01130         print ""
01131         
01132         print "on H nCap within LS pass H cut: ", self.counter["onHCapWithinLSPassCut"]
01133         print "on H nCap within LS pass lower H cut uncertainty: ", self.counter["onHCapWithinLSCutUncertainL"]
01134         print "on H nCap within LS pass upper H cut uncertainty: ", self.counter["onHCapWithinLSCutUncertainU"]
01135         print ""
01136 
01137         print "non H nCap within LS pass H cut: ", self.counter["nonHCapWithinLSPassHCut"]
01138         print "non H nCap within LS pass lower H cut uncertainty: ", self.counter["nonHCapWithinLSHCutUncertainL"]
01139         print "non H nCap within LS pass upper H cut uncertainty: ", self.counter["nonHCapWithinLSHCutUncertainU"]
01140         print ""
01141 
01142         print "on H nCap beyond LS pass on-H cut: ", self.counter["onHCapBeyondLSPassCut"]
01143         print "on H nCap beyond LS pass lower on-H cut uncertainty: ", self.counter["onHCapBeyondLSCutUncertainL"]
01144         print "on H nCap bdyond LS pass upper on-H cut uncertainty: ", self.counter["onHCapBeyondLSCutUncertainU"]
01145         print ""
01146 
01147         print "non H nCap beyond LS pass H cut: ", self.counter["nonHCapBeyondLSPassHCut"]
01148         print "non H nCap beyond LS pass lower H cut uncertainty: ", self.counter["nonHCapBeyondLSHCutUncertainL"]
01149         print "non H nCap beyond LS pass upper H cut uncertainty: ", self.counter["nonHCapBeyondLSHCutUncertainU"]
01150         print ""
01151 
01152         print "On Gd capture pass on-H cut: ", self.counter["nGdPassHCut"]
01153         print "On Gd capture pass lower on-H cut uncertainty: ", self.counter["nGdHCutUncertainL"]
01154         print "On Gd capture pass upper on-H cut uncertainty: ", self.counter["nGdHCutUncertainU"]
01155         print ""
01156 
01157         print "Total n-capture on H in MO: ", self.counter["onHCapMO"]
01158         print "Total n-capture on H in OAV: ", self.counter["onHCapOAV"]
01159         print "on H nCap in OAV&MO pass on-H cut: ", self.counter["onHCapOavMoPassCut"]
01160         print "on H nCap in OAV&MO pass lower on-H cut uncertainty: ", self.counter["onHCapOavMoCutUncertainL"]
01161         print "on H nCap in OAV&MO pass upper on-H cut uncertainty: ", self.counter["onHCapOavMoCutUncertainU"]
01162         print ""
01163         
01164         print "on H capture detection efficiency: ", self.counter["HCapDetEff"]
01165         print "on H capture detection uncertainty L: ", self.counter["HCapDetEffUncertainL"]
01166         print "on H capture detection uncertainy U: ", self.counter["HCapDetEffUncertainU"]
01167         print ""
01168         
01169         print "on H capture detection efficiency(canonical): ", self.counter["HCapCanoDetEff"]
01170         print "on H capture detection uncertainty(L:canonical):", self.counter["HCapCanoDetEffUncertainL"]
01171         print "on H capture detection uncertainty(U:canonical):", self.counter["HCapCanoDetEffUncertainU"]
01172         
01173         print "non H capture contamination uncertainty(L): ", self.counter["HCapDetEffContaminationUncertainL"]
01174         print "non H capture contamination uncertainty(U): ", self.counter["HCapDetEffContaminationUncertainU"]
01175         
01176         print "non H capture contamination:", self.counter["HCapDetContamination"]
01177 
01178         print "Finalizing ", self.name()
01179         status = GaudiAlgo.finalize(self)
01180         return status
01181 
01182 def configure(argv = []):
01183     import sys, getopt
01184     from time import localtime, gmtime, mktime, strftime, strptime, timezone
01185     opts,args = getopt.getopt(argv,"o:w:")
01186     outputrootfile = 'AdBasicPlots.root'
01187     wallTime = 0
01188     for opt,arg in opts:
01189         if opt == "-o":
01190             outputrootfile = arg
01191         if opt == "-w":
01192             if -1 != arg.find('T'):
01193                 wallTime = int(mktime(strptime(arg,
01194                                                DATETIME_FORMAT)) - timezone)
01195             else:
01196                 wallTime = int(arg)
01197 
01198     print "Your output ROOT files is: ", outputrootfile
01199     from DetHelpers.DetHelpersConf import CoordSysSvc
01200 #    from GaudiSvc.GaudiSvcConf import DetectorDataSvc
01201 
01202     from StatisticsSvc.StatisticsSvcConf import StatisticsSvc
01203     statsSvc = StatisticsSvc()
01204     statsSvc.Output ={"file1":outputrootfile}
01205 
01206     coorSvc = CoordSysSvc()
01207     coorSvc.OutputLevel = 1
01208 
01209     from Gaudi.Configuration import ApplicationMgr
01210     theApp = ApplicationMgr()
01211     theApp.ExtSvc.append(coorSvc)
01212 
01213     return
01214 
01215 def run(app):
01216     '''Configure and add this algorithm to job'''
01217 #    from Gaudi.Configuration import ApplicationMgr
01218 #    app = ApplicationMgr()
01219     app.ExtSvc += ["StatisticsSvc"]
01220     plotBasics = plotIbdBasics("myBasics")
01221     app.addAlgorithm(plotBasics)
01222     pass
In This Package:

AdPerformance.py
