00001
00002
00003
00004
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 TH3F = gbl.TH3F
00020 TH3D = gbl.TH3D
00021
00022 class plotIbdBasics(GaudiAlgo):
00023
00024 def __init__(self,name):
00025 GaudiAlgo.__init__(self,name)
00026 self.hist = {}
00027
00028 self.outputstat = open('IbdStat.txt', 'w')
00029
00030
00031 self.nInGdLS = 0
00032 self.nInIAV = 0
00033 self.nInLSO = 0
00034 self.nInOAV = 0
00035 self.nInMO = 0
00036
00037 self.nCap = 0
00038 self.nGen = 0
00039
00040 self.SpillIn_onGd = 0
00041 self.IavSpillIn_onGd = 0
00042 self.SpillOut = 0
00043 self.SpillIn_onH = 0
00044 self.OavSpillIn_onH = 0
00045 self.SpillOut_onH = 0
00046
00047 self.onGdCap = 0
00048 self.onGdGdsCap = 0
00049
00050 self.onGdCapPassCut = 0
00051 self.onGdGdsCapPassCut = 0
00052 self.nPassGdCut = 0
00053
00054 self.onGdCapCutUncertain = 0
00055 self.nGdCutUncertain = 0
00056 self.onGdGdsCapCutUncertain = 0
00057
00058 self.nonGdPassGdCut = 0
00059
00060 self.onHCap = 0
00061 self.onHCapWithinLS = 0
00062 self.onHCapWithinLSPassCut = 0
00063 self.onHCapWithinLSCutUncertainL = 0
00064 self.onHCapWithinLSCutUncertainU = 0
00065
00066 self.nonHCapWithinLSPassHCut = 0
00067 self.nonHCapWithinLSHCutUncertainL = 0
00068 self.nonHCapWithinLSHCutUncertainU = 0
00069
00070 self.onHCapBeyondLS = 0
00071 self.onHCapBeyondLSPassCut = 0
00072 self.onHCapBeyondLSCutUncertainL = 0
00073 self.onHCapBeyondLSCutUncertainU = 0
00074
00075 self.nonHCapBeyondLSPassHCut = 0
00076 self.nonHCapBeyondLSHCutUncertainL = 0
00077 self.nonHCapBeyondLSHCutUncertainU = 0
00078
00079 self.onHCapOavMoPassCut = 0
00080 self.onHCapOavMoCutUncertainL = 0
00081 self.onHCapOavMoCutUncertainU = 0
00082 self.onHCapLS = 0
00083 self.onHCapMO = 0
00084 self.onHCapOAV = 0
00085
00086 self.onHCapPassCut = 0
00087 self.onHCapCutUncertainL = 0
00088 self.onHCapCutUncertainU = 0
00089
00090 self.nPassHCut = 0
00091 self.nHCutUncertainL = 0
00092 self.nHCutUncertainU = 0
00093
00094 self.nGdPassHCut = 0
00095 self.nGdHCutUncertainL = 0
00096 self.nGdHCutUncertainU = 0
00097
00098 return
00099
00100 def initialize(self):
00101 print "Initializing the IBD basic ploter", self.name()
00102 status = GaudiAlgo.initialize(self)
00103 if status.isFailure(): return status
00104 self.target_de_name = '/dd/Structure/AD/db-ade1/db-sst1/db-oil1'
00105 self.gds_de_name = '/dd/Structure/AD/db-gds1'
00106 self.lso_de_name = '/dd/Structure/AD/db-lso1'
00107
00108
00109 self.coorSvc = self.svc('ICoordSysSvc', 'CoordSysSvc')
00110 if not self.coorSvc:
00111 print 'Failed to get CoordSysSvc'
00112 return FAILURE
00113
00114
00115
00116 self.statsSvc = self.svc('IStatisticsSvc','StatisticsSvc')
00117 if self.statsSvc == None:
00118 self.error("Failed to initialize IStat service.")
00119 return FAILURE
00120
00121 self.hist["genRZ"] = TH2F("genRZ", "Generation Vertex R-Z",
00122 100, 0.0, 6.190144, 100, -2.48, 2.48)
00123 status = self.statsSvc.put('/file1/basics/genRZ',
00124 self.hist["genRZ"])
00125 if status.isFailure(): return status
00126
00127 self.hist["genXY"] = TH2F("genXY", "Generation Vertex X-Y",
00128 100, -2.48, 2.48, 100, -2.48, 2.48)
00129 status = self.statsSvc.put('/file1/basics/genXY',
00130 self.hist["genXY"])
00131 if status.isFailure(): return status
00132
00133 self.hist["HitTime"] = TH1F("HitTime", "Hit Time [#mus]",
00134 500, 0.0, 2000)
00135 status = self.statsSvc.put('/file1/basics/HitTime',
00136 self.hist["HitTime"])
00137 if status.isFailure(): return status
00138
00139 self.hist["GdCapHitTime"] = TH1F("GdCapHitTime",
00140 "Gd Capture Hit Time [#mus]",
00141 500, 0.0, 2000)
00142 status = self.statsSvc.put('/file1/basics/GdCapHitTime',
00143 self.hist["GdCapHitTime"])
00144 if status.isFailure(): return status
00145
00146 self.hist["HCapHitTime"] = TH1F("HCapHitTime",
00147 "H Capture Hit Time [#mus]",
00148 500, 0.0, 2000)
00149 status = self.statsSvc.put('/file1/basics/HCapHitTime',
00150 self.hist["HCapHitTime"])
00151 if status.isFailure(): return status
00152
00153 self.hist["pe_E"] = TH2F("pe_E", "PE vs visible E (e+ within GdLS)",
00154 100, 0, 10, 700, 0, 1400)
00155 status = self.statsSvc.put('/file1/basics/pe_E',
00156 self.hist["pe_E"])
00157 if status.isFailure(): return status
00158
00159 self.hist["pe_E_inLS"] = TH2F("pe_E_inLS",
00160 "PE vs visible E (e+ within LS)",
00161 100, 0, 10, 700, 0, 1400)
00162 status = self.statsSvc.put('/file1/basics/pe_E_inLS',
00163 self.hist["pe_E_inLS"])
00164 if status.isFailure(): return status
00165
00166 self.hist["pe_E_LS"] = TH2F("pe_E_LS",
00167 "PE vs visible E (e+ in LS)",
00168 100, 0, 10, 700, 0, 1400)
00169 status = self.statsSvc.put('/file1/basics/pe_E_LS',
00170 self.hist["pe_E_LS"])
00171 if status.isFailure(): return status
00172
00173 self.hist["pe_E_all"] = TH2F("pe_E_all", "PE vs visible E (all e+)",
00174 100, 0, 10, 700, 0, 1400)
00175 status = self.statsSvc.put('/file1/basics/pe_E_all',
00176 self.hist["pe_E_all"])
00177 if status.isFailure(): return status
00178
00179 self.hist["nHits_LS"] = TH1F("nHits_LS",
00180 "nCap Hits in LS",
00181 100, 0, 1400)
00182 status = self.statsSvc.put('/file1/basics/nHits_LS',
00183 self.hist["nHits_LS"])
00184 if status.isFailure(): return status
00185
00186 self.hist["nHits_MO"] = TH1F("nHits_MO",
00187 "nCap Hits in MO",
00188 100, 0, 1400)
00189 status = self.statsSvc.put('/file1/basics/nHits_MO',
00190 self.hist["nHits_MO"])
00191 if status.isFailure(): return status
00192
00193 self.hist["nHits_onGd"] = TH1F("nHits_onGd",
00194 "on-Gd nCap Hits",
00195 100, 0, 1400)
00196 status = self.statsSvc.put('/file1/basics/nHits_onGd',
00197 self.hist["nHits_onGd"])
00198 if status.isFailure(): return status
00199
00200 self.hist["nHits_onH"] = TH1F("nHits_onH",
00201 "on-H nCap Hits",
00202 100, 0, 1400)
00203 status = self.statsSvc.put('/file1/basics/nHits_onH',
00204 self.hist["nHits_onH"])
00205 if status.isFailure(): return status
00206
00207 self.hist["nHits"] = TH1F("nHits", "Neutron Capture Hits in GdLS",
00208 100, 0, 1400)
00209 status = self.statsSvc.put('/file1/basics/nHits',
00210 self.hist["nHits"])
00211 if status.isFailure(): return status
00212
00213 self.hist["nHits_all"] = TH1F("nHits_all", "Neutron Capture Hits",
00214 100, 0, 1400)
00215 status = self.statsSvc.put('/file1/basics/nHits_all',
00216 self.hist["nHits_all"])
00217 if status.isFailure(): return status
00218
00219 self.hist["nGdCapGenPos"] = TH1F("nGdCapGenPos",
00220 "Neutron Gd Capture Generation Position",
00221 310, 0, 6.190144)
00222 status = self.statsSvc.put('/file1/basics/nGdCapGenPos',
00223 self.hist["nGdCapGenPos"])
00224 if status.isFailure(): return status
00225
00226 self.hist["nGdGdsCapGenPos"] = TH1F("nGdGdsCapGenPos",
00227 "Neutron Gd Capture in Gds Generation Position",
00228 310, 0, 6.190144)
00229 status = self.statsSvc.put('/file1/basics/nGdGdsCapGenPos',
00230 self.hist["nGdGdsCapGenPos"])
00231 if status.isFailure(): return status
00232
00233 self.hist["nGdCapPos"] = TH3F("nGdCapPos",
00234 "Neutron Gd Capture Position",
00235 100, -2.5, 2.5,
00236 100, -2.5, 2.5,
00237 100, -2.5, 2.5)
00238 status = self.statsSvc.put('/file1/basics/nGdCapPos',
00239 self.hist["nGdCapPos"])
00240 if status.isFailure(): return status
00241
00242 self.hist["nHCapGenPos"] = TH1F("nHCapGenPos",
00243 "Neutron H Capture Generation Position",
00244 310, 0, 6.190144)
00245 status = self.statsSvc.put('/file1/basics/nHCapGenPos',
00246 self.hist["nHCapGenPos"])
00247 if status.isFailure(): return status
00248
00249
00250 self.hist["nHOilCapGenPos"] = TH1F("nHOilCapGenPos",
00251 "Neutron H Capture Generation Position",
00252 310, 0, 6.190144)
00253 status = self.statsSvc.put('/file1/basics/nHOilCapGenPos',
00254 self.hist["nHOilCapGenPos"])
00255 if status.isFailure(): return status
00256
00257 self.hist["nHLsoCapGenPos"] = TH1F("nHLsoCapGenPos",
00258 "Neutron H Capture Generation Position",
00259 310, 0, 6.190144)
00260 status = self.statsSvc.put('/file1/basics/nHLsoCapGenPos',
00261 self.hist["nHLsoCapGenPos"])
00262 if status.isFailure(): return status
00263
00264
00265 self.hist["eDepInGdLS"] = TH1F("eDepInGdLS", "Deposited Energy [MeV]",
00266 20, 0, 20)
00267 status = self.statsSvc.put('/file1/basics/eDepInGdLS',
00268 self.hist["eDepInGdLS"])
00269 if status.isFailure(): return status
00270
00271 self.hist["eDepInLS"] = TH1F("eDepInLS", "Deposited Energy [MeV]",
00272 20, 0, 20)
00273 status = self.statsSvc.put('/file1/basics/eDepInLS',
00274 self.hist["eDepInLS"])
00275 if status.isFailure(): return status
00276
00277 self.hist["drift_GdLS"] = TH1F("drift_GdLS",
00278 "Neutron Drift Distance in GdLS [cm]",
00279 200, 0, 50)
00280 status = self.statsSvc.put('/file1/basics/drift_GdLS',
00281 self.hist["drift_GdLS"])
00282 if status.isFailure(): return status
00283
00284 self.hist["drift_LS"] = TH1F("drift_LS",
00285 "Neutron Drift Distance in LS [cm]",
00286 200, 0, 50)
00287 status = self.statsSvc.put('/file1/basics/drift_LS',
00288 self.hist["drift_LS"])
00289 if status.isFailure(): return status
00290
00291 self.hist["time_GdLS"] = TH1F("time_GdLS",
00292 "Neutron Capture time in GdLS [#mus]",
00293 400, 0, 400)
00294 status = self.statsSvc.put('/file1/basics/time_GdLS',
00295 self.hist["time_GdLS"])
00296 if status.isFailure(): return status
00297
00298 self.hist["time_LS"] = TH1F("time_LS",
00299 "Neutron Capture Time in LS [#mus]",
00300 400, 0, 2000)
00301 status = self.statsSvc.put('/file1/basics/time_LS',
00302 self.hist["time_LS"])
00303 if status.isFailure(): return status
00304
00305 return SUCCESS
00306
00307 def execute(self):
00308 print "Executing plotIbdBasics", self.name()
00309 evt = self.evtSvc()
00310 simhdr = evt['/Event/Sim/SimHeader']
00311
00312 if simhdr == None:
00313 print "No SimHeader in this ReadOut"
00314 return SUCCESS
00315
00316 det = self.detSvc(self.target_de_name)
00317 det_gds = self.detSvc(self.gds_de_name)
00318 det_lso = self.detSvc(self.lso_de_name)
00319
00320
00321 statshdr = simhdr.unobservableStatistics()
00322 stats = statshdr.stats()
00323
00324 PID_trk1 = stats["pdgId_Trk1"].sum()
00325 PID_trk2 = stats["pdgId_Trk2"].sum()
00326
00327 if PID_trk1 != -11 or PID_trk2 != 2112:
00328 print "PID of track 1 is", PID_trk1
00329 print "PID of track 2 is", PID_trk2
00330 print "Not an IBD event."
00331 return SUCCESS
00332
00333 tGen = stats["t_Trk2"].sum()
00334 xGen = stats["x_Trk2"].sum()
00335 yGen = stats["y_Trk2"].sum()
00336 zGen = stats["z_Trk2"].sum()
00337
00338 tCap = stats["tEnd_Trk2"].sum()
00339 xCap = stats["xEnd_Trk2"].sum()
00340 yCap = stats["yEnd_Trk2"].sum()
00341 zCap = stats["zEnd_Trk2"].sum()
00342
00343
00344 de = self.getDet(self.target_de_name)
00345 de_lso = self.getDet(self.lso_de_name)
00346 de_gds = self.getDet(self.gds_de_name)
00347 if not de:
00348 print 'Failed to get DE',self.target_de_name
00349 return FAILURE
00350
00351 if not de_lso:
00352 print 'Failed to get DE',self.lso_de_name
00353 return FAILURE
00354
00355 if not de_gds:
00356 print 'Failed to get DE',self.gds_de_name
00357 return FAILURE
00358
00359
00360 import PyCintex
00361 Gaudi = PyCintex.makeNamespace('Gaudi')
00362 genGlbPoint = Gaudi.XYZPoint(xGen, yGen, zGen)
00363 capGlbPoint = Gaudi.XYZPoint(xCap, yCap, zCap)
00364
00365 genLclPoint = de.geometry().toLocal(genGlbPoint)
00366 capLclPoint = de.geometry().toLocal(capGlbPoint)
00367 genLclPointLso = de_lso.geometry().toLocal(genGlbPoint)
00368 capLclPointLso = de_lso.geometry().toLocal(capGlbPoint)
00369 genLclPointGds = de_gds.geometry().toLocal(genGlbPoint)
00370 capLclPointGds = de_gds.geometry().toLocal(capGlbPoint)
00371
00372
00373 ndrift = ROOT.TVector3(xCap-xGen, yCap-yGen, zCap-zGen)
00374
00375 capTime = tCap - tGen
00376 capDis = ndrift.Mag()
00377
00378 R2 = genLclPoint.x()/units.meter * genLclPoint.x()/units.meter + \
00379 genLclPoint.y()/units.meter * genLclPoint.y()/units.meter
00380
00381 self.hist["genRZ"].Fill(R2, genLclPoint.z()/units.meter)
00382
00383 self.hist["genXY"].Fill(genLclPoint.x()/units.meter,genLclPoint.y()/units.meter)
00384
00385
00386 genDE = self.coorSvc.coordSysDE(genGlbPoint)
00387 capDE = self.coorSvc.coordSysDE(capGlbPoint)
00388 if not genDE:
00389 print 'Failed to find coordinate system DE for generation'\
00390 '[',genGlbPoint.x(),genGlbPoint.y(),genGlbPoint.z(),']'
00391 print 'Local: [',genLclPoint.x(),genLclPoint.y(),genLclPoint.z(),']'
00392 return FAILURE
00393 else:
00394 self.nGen += 1
00395 gendmvol = genDE.geometry().belongsToPath(genGlbPoint,-1)
00396
00397 if not capDE:
00398 print 'Failed to find coordinate system DE for capture'\
00399 '[',capGlbPoint.x(),capGlbPoint.y(),capGlbPoint.z(),']'
00400 return FAILURE
00401 else:
00402 self.nCap += 1
00403 capdmvol = capDE.geometry().belongsToPath(capGlbPoint,-1)
00404
00405 import re
00406 genDM = re.split('/', gendmvol).pop()
00407 capDM = re.split('/', capdmvol).pop()
00408 print "Generated in ", genDM
00409 print "Captured in ", capDM
00410
00411 positronHits = 0
00412 neutronHits = 0
00413 positronTimeCut = 500.
00414
00415 positronCut = 126.5
00416 positronCutL = 124.0
00417 positronCutR = 129.0
00418
00419 onHLowerCut = 194.
00420 onHLowerCutL = 190.1
00421 onHLowerCutR = 197.9
00422
00423 onHUpperCut = 453.
00424 onHUpperCutL = 443.9
00425 onHUpperCutR = 462.1
00426
00427 onGdCut = 811.9
00428 onGdCutL = 803.8
00429 onGdCutR = 820.0
00430
00431
00432 vis_trk1 = 1.022 + stats['ke_Trk1'].sum()/units.MeV
00433 self.hist["pe_E_all"].Fill(vis_trk1, positronHits)
00434
00435
00436 capTarget = stats["capTarget"].sum()
00437 print "The capture target is ", capTarget
00438
00439
00440 eDepInGdLS = stats["EDepInGdLS"].sum()
00441 self.hist["eDepInGdLS"].Fill(eDepInGdLS/units.MeV)
00442 eDepInLS = stats["EDepInLS"].sum()
00443 self.hist["eDepInLS"].Fill(eDepInLS/units.MeV)
00444
00445 simhits = simhdr.hits()
00446 for detector in simhits.hitDetectors():
00447 hitCollection = simhits.hitsByDetector(detector)
00448 if hitCollection == None:
00449 print "No hits in ", detector
00450 hits = hitCollection.collection()
00451 for hit in hits:
00452
00453 self.hist["HitTime"].Fill(hit.hitTime()/units.microsecond)
00454 if capTarget == 1:
00455 self.hist["HCapHitTime"].Fill(hit.hitTime()/units.microsecond)
00456 if capTarget == 64:
00457 self.hist["GdCapHitTime"].Fill(hit.hitTime()/units.microsecond)
00458 if hit.hitTime()/units.nanosecond<positronTimeCut and \
00459 hit.hitTime()/units.nanosecond<capTime/units.nanosecond:
00460 positronHits += 1
00461 else:
00462 neutronHits += 1
00463
00464 self.hist["nHits_all"].Fill(neutronHits)
00465
00466 if genDM == 'db-gds1':
00467 self.nInGdLS += 1
00468 self.hist["pe_E"].Fill(vis_trk1, positronHits)
00469 if capDM == 'db-gds1':
00470 self.hist["nHits"].Fill(neutronHits)
00471 self.hist["time_GdLS"].Fill(capTime/units.microsecond)
00472 self.hist["drift_GdLS"].Fill(capDis/units.cm)
00473 else:
00474 self.SpillOut += 1
00475
00476 if genDM == 'db-iav1':
00477 self.nInIAV += 1
00478 if capTarget == 64:
00479 self.IavSpillIn_onGd += 1
00480
00481 if genDM == 'db-lso1':
00482 self.nInLSO += 1
00483 self.hist["pe_E_LS"].Fill(vis_trk1, positronHits)
00484 if capDM == 'db-lso1':
00485 self.hist["time_LS"].Fill(capTime/units.microsecond)
00486 self.hist["drift_LS"].Fill(capDis/units.cm)
00487 if capDM == 'db-oil1' or capDM == 'db-oav1':
00488 self.SpillOut_onH += 1
00489 if capTarget == 64:
00490 self.SpillIn_onGd += 1
00491
00492 if genDM == 'db-oav1':
00493 self.nInOAV += 1
00494 if re.search('db-lso1', capdmvol) and capTarget == 1:
00495 self.OavSpillIn_onH += 1
00496 if neutronHits > onHLowerCut and neutronHits < onHUpperCut:
00497 self.onHCapOavMoPassCut += 1
00498 if neutronHits > onHLowerCutL and neutronHits < onHLowerCutR:
00499 self.onHCapOavMoCutUncertainL += 1
00500 if neutronHits > onHUpperCutL and neutronHits < onHUpperCutR:
00501 self.onHCapOavMoCutUncertainU += 1
00502
00503 if genDM == 'db-oil1':
00504 self.nInMO += 1
00505 if re.search('db-lso1', capdmvol) and capTarget == 1:
00506 self.SpillIn_onH += 1
00507 if neutronHits > onHLowerCut and neutronHits < onHUpperCut:
00508 self.onHCapOavMoPassCut += 1
00509 if neutronHits > onHLowerCutL and neutronHits < onHLowerCutR:
00510 self.onHCapOavMoCutUncertainL += 1
00511 if neutronHits > onHUpperCutL and neutronHits < onHUpperCutR:
00512 self.onHCapOavMoCutUncertainU += 1
00513
00514 if re.search('db-lso1',gendmvol):
00515 self.hist["pe_E_inLS"].Fill(vis_trk1, positronHits)
00516
00517 if re.search('db-lso1',capdmvol):
00518 if capTarget == 1:
00519 self.onHCapWithinLS += 1
00520 if neutronHits > onHLowerCut and neutronHits < onHUpperCut:
00521 self.onHCapWithinLSPassCut += 1
00522 if neutronHits > onHLowerCutL and neutronHits < onHLowerCutR:
00523 self.onHCapWithinLSCutUncertainL += 1
00524 if neutronHits > onHUpperCutL and neutronHits < onHUpperCutR:
00525 self.onHCapWithinLSCutUncertainU += 1
00526 else:
00527 if neutronHits > onHLowerCut and neutronHits < onHUpperCut:
00528 self.nonHCapWithinLSPassHCut += 1
00529 if neutronHits > onHLowerCutL and neutronHits < onHLowerCutR:
00530 self.nonHCapWithinLSHCutUncertainL += 1
00531 if neutronHits > onHUpperCutL and neutronHits < onHUpperCutR:
00532 self.nonHCapWithinLSHCutUncertainU += 1
00533
00534 else:
00535 if capTarget == 1:
00536 self.onHCapBeyondLS += 1
00537 if neutronHits > onHLowerCut and neutronHits < onHUpperCut:
00538 self.onHCapBeyondLSPassCut += 1
00539 if neutronHits > onHLowerCutL and neutronHits < onHLowerCutR:
00540 self.onHCapBeyondLSCutUncertainL += 1
00541 if neutronHits > onHUpperCutL and neutronHits < onHUpperCutR:
00542 self.onHCapBeyondLSCutUncertainU += 1
00543 else:
00544 if neutronHits > onHLowerCut and neutronHits < onHUpperCut:
00545 self.nonHCapBeyondLSPassHCut += 1
00546 if neutronHits > onHLowerCutL and neutronHits < onHLowerCutR:
00547 self.nonHCapBeyondLSHCutUncertainL += 1
00548 if neutronHits > onHUpperCutL and neutronHits < onHUpperCutR:
00549 self.nonHCapBeyondLSHCutUncertainU += 1
00550
00551
00552 if neutronHits > onGdCut:
00553 self.nPassGdCut += 1
00554 if capTarget != 64:
00555 print "non-Gd capture pass Gd cut: ", capTarget
00556 self.nonGdPassGdCut += 1
00557 value = ('Non-Gd capTarget: ', capTarget)
00558 theline = str(value)
00559 self.outputstat.write(theline)
00560 self.outputstat.write('\n')
00561
00562 if neutronHits > onGdCutL and neutronHits < onGdCutR:
00563 self.nGdCutUncertain += 1
00564
00565
00566 if neutronHits > onHLowerCut and neutronHits < onHUpperCut:
00567 self.nPassHCut += 1
00568 if capTarget == 64:
00569 self.nGdPassHCut += 1
00570 if neutronHits > onHLowerCutL and neutronHits < onHLowerCutR:
00571 self.nHCutUncertainL += 1
00572 if capTarget == 64:
00573 self.nGdHCutUncertainL += 1
00574 if neutronHits > onHUpperCutL and neutronHits < onHUpperCutR:
00575 self.nHCutUncertainU += 1
00576 if capTarget == 64:
00577 self.nGdHCutUncertainU += 1
00578
00579 if capTarget == 64:
00580 self.onGdCap += 1
00581 self.hist["nHits_onGd"].Fill(neutronHits)
00582 self.hist["nGdCapPos"].Fill(capLclPoint.x()/units.meter,capLclPoint.y()/units.meter,capLclPoint.z()/units.meter)
00583 if capDM == 'db-gds1':
00584 self.onGdGdsCap += 1
00585 if neutronHits > onGdCut:
00586 self.onGdGdsCapPassCut += 1
00587 if neutronHits > onGdCutL and neutronHits < onGdCutR:
00588 self.onGdGdsCapCutUncertain += 1
00589
00590 if genLclPoint.z()/units.m < 1. and genLclPoint.z()/units.m > -1.:
00591 self.hist["nGdCapGenPos"].Fill(R2)
00592 if capDM == 'db-gds1':
00593 self.hist["nGdGdsCapGenPos"].Fill(R2)
00594
00595 if neutronHits > onGdCut:
00596 self.onGdCapPassCut += 1
00597 if neutronHits > onGdCutL and neutronHits < onGdCutR:
00598 self.onGdCapCutUncertain += 1
00599
00600 if capTarget == 1:
00601 self.onHCap += 1
00602 self.hist["nHits_onH"].Fill(neutronHits)
00603
00604 if genLclPoint.z()/units.m < 1. and genLclPoint.z()/units.m > -1.:
00605 self.hist["nHCapGenPos"].Fill(R2)
00606 if capDM == 'db-lso1':
00607 self.hist["nHLsoCapGenPos"].Fill(R2)
00608 if capDM == 'db-oil1' or capDM == 'db-oav1':
00609 self.hist["nHOilCapGenPos"].Fill(R2)
00610
00611 if neutronHits > onHLowerCut and neutronHits < onHUpperCut:
00612 self.onHCapPassCut += 1
00613 if neutronHits > onHLowerCutL and neutronHits < onHLowerCutR:
00614 self.onHCapCutUncertainL += 1
00615 if neutronHits > onHUpperCutL and neutronHits < onHUpperCutR:
00616 self.onHCapCutUncertainU += 1
00617
00618 if capDM == 'db-lso1':
00619 self.hist["nHits_LS"].Fill(neutronHits)
00620 if capTarget == 1:
00621 self.onHCapLS += 1
00622
00623 if capDM == 'db-oav1':
00624 if capTarget == 1:
00625 self.onHCapOAV += 1
00626
00627 if capDM == 'db-oil1':
00628 self.hist["nHits_MO"].Fill(neutronHits)
00629 if capTarget == 1:
00630 self.onHCapMO += 1
00631
00632 return SUCCESS
00633
00634 def finalize(self):
00635 self.outputstat.close()
00636
00637 print "Total n generation: ", self.nGen
00638 print "Total n gen in GdLS: ", self.nInGdLS
00639 print "Total n gen in IAV: ", self.nInIAV
00640 print "Total n gen in LSO: ", self.nInLSO
00641 print "Total n gen in OAV: ", self.nInOAV
00642 print "Total n gen in MO: ", self.nInMO
00643 print "Total n capture: ", self.nCap
00644 print ""
00645
00646 print "Total n-capture on Gd: ", self.onGdCap
00647 print "Total n-capture on Gd in Gds: ", self.onGdGdsCap
00648 print "Spill-In on Gd n from LS: ", self.SpillIn_onGd
00649 print "Spill-In on Gd from IAV: ", self.IavSpillIn_onGd
00650 print "Spill-Out of n from Gd-LS: ", self.SpillOut
00651 print ""
00652
00653 print "On Gd capture pass on-Gd cut: ", self.onGdCapPassCut
00654 print "On Gd capture on-Gd cut uncertainty: ", self.onGdCapCutUncertain
00655 print "On Gd capture in GDS pass on-Gd cut: ", self.onGdGdsCapPassCut
00656 print "On Gd capture in GDS on-Gd cut uncertainty: ", self.onGdGdsCapCutUncertain
00657 print "Total non-Gd capture pass Gd cut: ", self.nonGdPassGdCut
00658 print "Total pass Gd cut: ", self.nPassGdCut
00659 print "Overall Gd cut uncertainty: ", self.nGdCutUncertain
00660 print ""
00661
00662 print "Total n-capture on H: ", self.onHCap
00663 print "Total n-capture on H in LS: ", self.onHCapLS
00664 print "Total n-capture on H within LS: ", self.onHCapWithinLS
00665 print ""
00666
00667 print "Spill-In of n from OIL to LSO: ", self.SpillIn_onH
00668 print "Spill-In of n from OAV to LSO: ", self.OavSpillIn_onH
00669 print "Spill-Out of n from LSO to OAV or OIL: ", self.SpillOut_onH
00670 print ""
00671
00672 print "n capture pass on-H cut: ", self.nPassHCut
00673 print "n capture pass lower on-H cut uncertainty: ", self.nHCutUncertainL
00674 print "n capture pass upper on-H cut uncertainty: ", self.nHCutUncertainU
00675 print ""
00676
00677 print "On H capture pass on-H cut: ", self.onHCapPassCut
00678 print "On H capture lower on-H cut uncertainty: ", self.onHCapCutUncertainL
00679 print "On H capture upper on-H cut uncertainty: ", self.onHCapCutUncertainU
00680 print ""
00681
00682 print "on H nCap within LS pass H cut: ", self.onHCapWithinLSPassCut
00683 print "on H nCap within LS pass lower H cut uncertainty: ", self.onHCapWithinLSCutUncertainL
00684 print "on H nCap within LS pass upper H cut uncertainty: ", self.onHCapWithinLSCutUncertainU
00685 print ""
00686
00687 print "non H nCap within LS pass H cut: ", self.nonHCapWithinLSPassHCut
00688 print "non H nCap within LS pass lower H cut uncertainty: ", self.nonHCapWithinLSHCutUncertainL
00689 print "non H nCap within LS pass upper H cut uncertainty: ", self.nonHCapWithinLSHCutUncertainU
00690 print ""
00691
00692 print "on H nCap beyond LS pass on-H cut: ", self.onHCapBeyondLSPassCut
00693 print "on H nCap beyond LS pass lower on-H cut uncertainty: ", self.onHCapBeyondLSCutUncertainL
00694 print "on H nCap bdyond LS pass upper on-H cut uncertainty: ", self.onHCapBeyondLSCutUncertainU
00695 print ""
00696
00697 print "non H nCap beyond LS pass H cut: ", self.nonHCapBeyondLSPassHCut
00698 print "non H nCap beyond LS pass lower H cut uncertainty: ", self.nonHCapBeyondLSHCutUncertainL
00699 print "non H nCap beyond LS pass upper H cut uncertainty: ", self.nonHCapBeyondLSHCutUncertainU
00700 print ""
00701
00702 print "On Gd capture pass on-H cut: ", self.nGdPassHCut
00703 print "On Gd capture pass lower on-H cut uncertainty: ", self.nGdHCutUncertainL
00704 print "On Gd capture pass upper on-H cut uncertainty: ", self.nGdHCutUncertainU
00705 print ""
00706
00707 print "Total n-capture on H in MO: ", self.onHCapMO
00708 print "Total n-capture on H in OAV: ", self.onHCapOAV
00709 print "on H nCap in OAV&MO pass on-H cut: ", self.onHCapOavMoPassCut
00710 print "on H nCap in OAV&MO pass lower on-H cut uncertainty: ", self.onHCapOavMoCutUncertainL
00711 print "on H nCap in OAV&MO pass upper on-H cut uncertainty: ", self.onHCapOavMoCutUncertainU
00712 print ""
00713
00714 print "Finalizing ", self.name()
00715 status = GaudiAlgo.finalize(self)
00716 return status
00717
00718 def configure():
00719 from DetHelpers.DetHelpersConf import CoordSysSvc
00720
00721
00722 from StatisticsSvc.StatisticsSvcConf import StatisticsSvc
00723 statsSvc = StatisticsSvc()
00724 statsSvc.Output ={"file1":"IbdBasicPlots.root"}
00725
00726
00727
00728
00729 coorSvc = CoordSysSvc()
00730 coorSvc.OutputLevel = 1
00731
00732 from Gaudi.Configuration import ApplicationMgr
00733 theApp = ApplicationMgr()
00734 theApp.ExtSvc.append(coorSvc)
00735
00736 return
00737
00738 def run(app):
00739 '''Configure and add this algorithm to job'''
00740
00741
00742
00743 app.ExtSvc += ["StatisticsSvc"]
00744 plotBasics = plotIbdBasics("myBasics")
00745 app.addAlgorithm(plotBasics)
00746 pass
