ATLAS Offline Software
Loading...
Searching...
No Matches
LArG4ShowerLibFunctions.py
Go to the documentation of this file.
1# Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
2
3
4import logging
5
6__author__ = 'Radist Morse radist.morse@gmail.com'
7
8
9class FourVector() :
10 def __init__(self) :
11 self.x = 0.0
12 self.y = 0.0
13 self.z = 0.0
14 self.e = 0.0
15
16class StoredTestShower() :
17 def __init__(self):
18 self.vertex = FourVector()
19 self.momentum = FourVector()
20 self.rsize = 0.0
21 self.zsize = 0.0
22 self.shower = [] #hits, list of four vectors (which are actually 5-vector in this case)
23
24class StoredEnergyShower() :
25 def __init__(self):
26 self.shower = [] #hits, list of four vectors (which are actually 5-vector in this case)
27 self.egen = 0.0
28 self.rsize = 0.0
29 self.zsize = 0.0
30 def __hash__(self): #needed for making a set
31 return hash(self.egen)
32 def __eq__(self,other): #same here
33 return self.egen == other.egen
34
35class TestShowerLib() :
36 def __init__(self) :
37 self.library = [] # list of StoredTestShower objs
38 self.detector= ""
39 self.particle= ""
40 self.release= ""
41 self.geometry= ""
42 self.geant= ""
43 self.phys= ""
44 self.comment= ""
45 def fromLibs(self,libs) :
46 for lib in libs :
47 if not isinstance(lib,self.__class__):
48 print ("ERROR: Different types of libs")
49 return False
50 self.detector = libs[0].detector
51 self.particle = libs[0].particle
52 self.release = libs[0].release
53 self.geometry = libs[0].geometry
54 self.geant = libs[0].geant
55 self.phys = libs[0].phys
56 for lib in libs :
57 if ( self.detector != lib.detector or
58 self.particle != lib.particle or
59 self.release != lib.release or
60 self.geometry != lib.geometry or
61 self.geant != lib.geant or
62 self.phys != lib.phys ) :
63 print ("ERROR: DIFFERENT LIBS!!!")
64 return False
65 from datetime import datetime
66 self.comment = "ADDED "+str(datetime.now())
67 for lib in libs :
68 self.library += lib.library
69 return True
70 def readFromFile(self,filename) :
71 from ROOT import TFile
72 tfile = TFile(filename)
73 try :
74 ver = int(tfile.Get("version").GetVal())
75 except Exception:
76 print ("Not a TestShowerLib: Broken file")
77 tfile.Close()
78 return False
79 if (ver != 10) : #<<<<<<<<<<<<<<<<<<<<<<-------------- lib ver
80 print ("Not a TestShowerLib")
81 tfile.Close()
82 return False
83 meta = tfile.Get("meta")
84 libr = tfile.Get("library")
85 for event in meta :
86 self.detector=str(event.detector)
87 self.particle=str(event.particle)
88 self.release=str(event.release)
89 self.geometry=str(event.geometry)
90 self.geant=str(event.geantVersion)
91 self.phys=str(event.physicsList)
92 self.comment=str(event.comment)
93
94 state = 0
95
96 #state == 0 : first line of shower header
97 #state == 1 : second line of sower header
98 #state == 2 : hit
99
100 for event in libr : #this is quite unclear, but easy to implement
101 if (state == 0) : #first line
102 curShower = StoredTestShower()
103 curShower.vertex.x = event.x
104 curShower.vertex.y = event.y
105 curShower.vertex.z = event.z
106 curShower.zsize = event.time
107 hitsInCurShower = event.e
108 state = 1
109 elif (state == 1) : #shower header
110 curShower.momentum.x = event.x
111 curShower.momentum.y = event.y
112 curShower.momentum.z = event.z
113 curShower.momentum.e = event.e
114 curShower.rsize = event.time
115 if (hitsInCurShower > 0) :
116 state = 2
117 else :
118 self.library.append(curShower)
119 state = 0
120 elif (state == 2) :
121 hit = FourVector()
122 hit.e = event.e
123 hit.x = event.x
124 hit.y = event.y
125 hit.z = event.z
126 hit.time = event.time
127 curShower.shower.append(hit)
128 hitsInCurShower -= 1
129 if (hitsInCurShower == 0) : #last hit
130 self.library.append(curShower)
131 state = 0
132 tfile.Close()
133 if (state != 0) :
134 print ("FILE CORRUPTED!!")
135 return False
136 return True
137 def writeToFile(self,filename) :
138 from ROOT import TFile,TTree,TParameter
139 from ROOT import gROOT, addressof
140 gROOT.ProcessLine(
141 "struct MyMetaStruct {\
142 Char_t detector[40];\
143 Char_t release[40];\
144 Char_t geometry[40];\
145 Char_t geant[40];\
146 Char_t phys[40];\
147 Char_t comment[400];\
148 Int_t particle;\
149 };" )
150 from ROOT import MyMetaStruct
151 gROOT.ProcessLine(
152 "struct MyStruct {\
153 Float_t x;\
154 Float_t y;\
155 Float_t z;\
156 Float_t e;\
157 Float_t time;\
158 };" )
159 from ROOT import MyStruct
160
161 tfile = TFile(filename,"RECREATE")
162
163 ver = TParameter(int)("version",10) #<<<<<<<<<<<<<<<<<<<<<<-------------- lib ver
164 ver.Write("version")
165
166 meta = TTree()
167 libr = TTree()
168
169 mmstruct = MyMetaStruct()
170
171 mmstruct.detector = "%s" % (str(self.detector))
172 mmstruct.particle = int(self.particle)
173 mmstruct.release = "%s" % (str(self.release))
174 mmstruct.geometry = "%s" % (str(self.geometry))
175 mmstruct.geant = "%s" % (str(self.geant))
176 mmstruct.phys = "%s" % (str(self.phys))
177 mmstruct.comment = "%s" % (str(self.comment))
178
179 meta.Branch("detector",addressof(mmstruct,"detector"),"detector/C")
180 meta.Branch("particle",addressof(mmstruct,"particle"),"particle/I")
181 meta.Branch("release",addressof(mmstruct,"release"),"release/C")
182 meta.Branch("geometry",addressof(mmstruct,"geometry"),"geometry/C")
183 meta.Branch("geantVersion",addressof(mmstruct,"geant"),"geantVersion/C")
184 meta.Branch("physicsList",addressof(mmstruct,"phys"),"physicsList/C")
185 meta.Branch("comment",addressof(mmstruct,"comment"),"physicsList/C")
186
187 meta.Fill()
188
189 mstruct = MyStruct()
190
191 libr.Branch("x",addressof(mstruct,"x"),"x/F")
192 libr.Branch("y",addressof(mstruct,"y"),"y/F")
193 libr.Branch("z",addressof(mstruct,"z"),"z/F")
194 libr.Branch("e",addressof(mstruct,"e"),"e/F")
195 libr.Branch("time",addressof(mstruct,"time"),"time/F")
196
197 for storedShower in self.library :
198 mstruct.x = storedShower.vertex.x
199 mstruct.y = storedShower.vertex.y
200 mstruct.z = storedShower.vertex.z
201 mstruct.time = storedShower.zsize
202 mstruct.e = len(storedShower.shower)
203 libr.Fill()
204 mstruct.x = storedShower.momentum.x
205 mstruct.y = storedShower.momentum.y
206 mstruct.z = storedShower.momentum.z
207 mstruct.e = storedShower.momentum.e
208 mstruct.time = storedShower.rsize
209 libr.Fill()
210 for hit in storedShower.shower:
211 mstruct.e = hit.e
212 mstruct.x = hit.x
213 mstruct.y = hit.y
214 mstruct.z = hit.z
215 mstruct.time = hit.time
216 libr.Fill()
217 meta.Write("meta")
218 libr.Write("library")
219 tfile.Close()
220 def printInfo(self) :
221 pass
222 def drawHits(self):
223 from ROOT import TH3F
224 from math import sqrt,copysign,log10
225 hits = TH3F("HITS","Hits Distrib",50,1,1000,101,-300,300,100,0,500)
226 containmentZ = TH3F("CONTZ","ContZ Distrib",50,1,1000,101,-300,300,100,0,500)
227 containmentR = TH3F("CONTR","ContR Distrib",50,1,1000,101,-300,300,100,0,500)
228 for storedShower in self.library :
229 containmentR.Fill(log10(storedShower.momentum.e)*333,storedShower.rsize,storedShower.zsize/2,10)
230 containmentR.Fill(log10(storedShower.momentum.e)*333,-storedShower.rsize,storedShower.zsize/2,10)
231 containmentZ.Fill(log10(storedShower.momentum.e)*333,0,storedShower.zsize,10)
232 for hit in storedShower.shower :
233 hits.Fill(log10(storedShower.momentum.e)*333,copysign(sqrt(hit.x*hit.x + hit.y*hit.y),hit.x),hit.z)
234 return hits,containmentZ,containmentR
235
236class EtaEnergyShowerLib() :
237 def __init__(self) :
238 self.library = {} # key (float) - eta, value (list) - list of StoredEnergyShower objs
239 self.detector= ""
240 self.particle= ""
241 self.release= ""
242 self.geometry= ""
243 self.geant= ""
244 self.phys= ""
245 self.comment= ""
246 def scaleEnergy(self,scalefactor) :
247 for etabin in self.library.values():
248 for storedShower in etabin :
249 for hit in storedShower.shower :
250 hit.e *= scalefactor
251 self.comment += " SCALED: "+str(scalefactor)
252 def truncate(self,truncate) :
253 showers = []
254 for eta,etabin in self.library.items():
255 for storedShower in etabin :
256 showers += [(eta,storedShower)]
257 if len(showers) <= truncate :
258 print ("WARNING: Size of the library is already less:",truncate,"<",len(showers))
259 return
260 from random import randint
261 while (len(showers) > truncate) :
262 rand = randint(0,len(showers)-1)
263 self.library[showers[rand][0]].remove(showers[rand][1])
264 del showers[rand]
265 return
266 def fromLibs(self,libs) :
267 for lib in libs :
268 if not isinstance(lib,self.__class__):
269 print ("ERROR: Different types of libs")
270 return False
271 self.detector = libs[0].detector
272 self.particle = libs[0].particle
273 self.release = libs[0].release
274 self.geometry = libs[0].geometry
275 self.geant = libs[0].geant
276 self.phys = libs[0].phys
277 self.mineta = libs[0].mineta
278 self.maxeta = libs[0].maxeta
279 self.comment = libs[0].comment
280 etas = set(libs[0].library.keys())
281 for lib in libs :
282 if ( self.detector != lib.detector or
283 self.particle != lib.particle or
284 self.release != lib.release or
285 self.geometry != lib.geometry or
286 self.geant != lib.geant or
287 self.phys != lib.phys or
288 self.mineta != lib.mineta or
289 self.maxeta != lib.maxeta or
290 etas != set(lib.library.keys()) ) :
291 print ("ERROR: DIFFERENT LIBS!!!")
292 return False
293 for lib in libs :
294 for k,v in lib.library.items():
295 self.library.setdefault(k,set()).update(v)
296 for k,v in self.library.items():
297 self.library[k] = list(v)
298 return True
299 def moveEta(self,oldEta,newEta) :
300 if not (oldEta in self.library.keys()) :
301 return False
302 self.library[newEta] = self.library.pop(oldEta)
303 return True
304 def removeEta(self,eta) :
305 if not (eta in self.library.keys()) :
306 return False
307 self.library.pop(eta)
308 return True
309 def readFromFile(self,filename) :
310 from ROOT import TFile
311 #from sets import Set
312 tfile = TFile(filename)
313 try:
314 ver = int(tfile.Get("version").GetVal())
315 except Exception:
316 print ("Not an EtaEnergyLib: Broken file")
317 tfile.Close()
318 return False
319
320 if (ver != 1) : #<<<<<<<<<<<<<<<<<<<<<<-------------- lib ver
321 print ("Not an EtaEnergyLib")
322 tfile.Close()
323 return False
324 meta = tfile.Get("meta")
325 libr = tfile.Get("library")
326
327 for event in meta :
328 self.detector=str(event.detector)
329 self.particle=str(event.particle)
330 self.release=str(event.release)
331 self.geometry=str(event.geometry)
332 self.geant=str(event.geantVersion)
333 self.phys=str(event.physicsList)
334 self.comment=str(event.comment)
335
336 state = 0
337 lastShower = False
338
339 for event in libr : #this is quite unclear, but easy to implement
340 if (state == 0) : #eta bin header
341 showersInCurEta = event.x
342 curEta = round(event.y,4)
343 self.mineta = event.z
344 self.maxeta = event.e
345 self.library[curEta] = []
346 if (showersInCurEta > 0) :
347 state = 1 #go to shower header
348 elif (state == 1) : #shower header
349 hitsInCurShower = event.x
350 rSize = event.y
351 zSize = event.z
352 genEnergy = event.e
353 showersInCurEta -= 1
354 if (showersInCurEta == 0) : #last shower
355 lastShower = True
356 curShower = StoredEnergyShower()
357 curShower.egen = genEnergy
358 curShower.rsize = rSize
359 curShower.zsize = zSize
360 #curShower["hits"] = []
361 if (hitsInCurShower > 0) :
362 state = 2 #go to hits
363 else : #empty shower
364 self.library[curEta].append(curShower)
365 if (lastShower) : #special case of last shower in bin being the empty one
366 lastShower = False
367 state = 0 #next bin
368 elif (state == 2) :
369 hit = FourVector()
370 hit.e = event.e
371 hit.x = event.x
372 hit.y = event.y
373 hit.z = event.z
374 hit.time = event.time
375 curShower.shower.append(hit)
376 hitsInCurShower -= 1
377 if (hitsInCurShower == 0) : #last hit
378 self.library[curEta].append(curShower)
379 if (lastShower) : # end of eta bin
380 lastShower = False
381 state = 0
382 else : #not yet
383 state = 1
384 tfile.Close()
385 if (state != 0) :
386 print ("FILE CORRUPTED!!")
387 return False
388 return True
389 def writeToFile(self,filename) :
390 from ROOT import TFile,TTree,TParameter
391 from ROOT import gROOT, addressof
392 gROOT.ProcessLine(
393 "struct MyMetaStruct {\
394 Char_t detector[40];\
395 Char_t release[40];\
396 Char_t geometry[40];\
397 Char_t geant[40];\
398 Char_t phys[40];\
399 Char_t comment[400];\
400 Int_t particle;\
401 };" )
402 from ROOT import MyMetaStruct
403 gROOT.ProcessLine(
404 "struct MyStruct {\
405 Float_t x;\
406 Float_t y;\
407 Float_t z;\
408 Float_t e;\
409 Float_t time;\
410 };" )
411 from ROOT import MyStruct
412
413 tfile = TFile(filename,"RECREATE")
414
415 ver = TParameter(int)("version",1) #<<<<<<<<<<<<<<<<<<<<<<-------------- lib ver
416 ver.Write("version")
417
418 meta = TTree()
419 libr = TTree()
420
421 mmstruct = MyMetaStruct()
422
423 mmstruct.detector = "%s" % (str(self.detector))
424 mmstruct.particle = int(self.particle)
425 mmstruct.release = "%s" % (str(self.release))
426 mmstruct.geometry = "%s" % (str(self.geometry))
427 mmstruct.geant = "%s" % (str(self.geant))
428 mmstruct.phys = "%s" % (str(self.phys))
429 mmstruct.comment = "%s" % (str(self.comment))
430
431 meta.Branch("detector",addressof(mmstruct,"detector"),"detector/C")
432 meta.Branch("particle",addressof(mmstruct,"particle"),"particle/I")
433 meta.Branch("release",addressof(mmstruct,"release"),"release/C")
434 meta.Branch("geometry",addressof(mmstruct,"geometry"),"geometry/C")
435 meta.Branch("geantVersion",addressof(mmstruct,"geant"),"geantVersion/C")
436 meta.Branch("physicsList",addressof(mmstruct,"phys"),"physicsList/C")
437 meta.Branch("comment",addressof(mmstruct,"comment"),"physicsList/C")
438
439 meta.Fill()
440
441 mstruct = MyStruct()
442
443 libr.Branch("x",addressof(mstruct,"x"),"x/F")
444 libr.Branch("y",addressof(mstruct,"y"),"y/F")
445 libr.Branch("z",addressof(mstruct,"z"),"z/F")
446 libr.Branch("e",addressof(mstruct,"e"),"e/F")
447 libr.Branch("time",addressof(mstruct,"time"),"time/F")
448
449 etas = self.library.keys()
450
451 for eta in sorted(etas) :
452 mstruct.x = len(self.library[eta])
453 mstruct.y = eta
454 mstruct.z = self.mineta
455 mstruct.e = self.maxeta
456 mstruct.time = 0
457 libr.Fill()
458 self.library[eta].sort(key=lambda x: x.egen)
459 for storedShower in self.library[eta] :
460 mstruct.x = len(storedShower.shower)
461 mstruct.y = storedShower.rsize
462 mstruct.z = storedShower.zsize
463 mstruct.e = storedShower.egen
464 mstruct.time = 0
465 libr.Fill()
466 for hit in storedShower.shower:
467 mstruct.e = hit.e
468 mstruct.x = hit.x
469 mstruct.y = hit.y
470 mstruct.z = hit.z
471 mstruct.time = hit.time
472 libr.Fill()
473 meta.Write("meta")
474 libr.Write("library")
475 tfile.Close()
476 def printInfo(self) :
477 print ("VERSION: EtaEnergyLib","PARTICLE:",self.particle,"DETECTOR:",self.detector)
478 print (self.release, self.geometry, self.geant, self.phys)
479 print (self.comment)
480 ebins = [1,2,5,10,20,50,100,200,500,1000]
481 etas = sorted(self.library.keys())
482 print ("Number of etabins:",str(len(etas)))
483 print ("MinEta:",self.mineta,"MaxEta:",self.maxeta)
484 fstot = 0
485 for eta in etas :
486 fstot +=len(self.library[eta])
487 print ("Number of showers:",str(fstot))
488 print ("-"*(12+len(ebins)*8)) #horizontal line
489 infostr = "|etas|ebins|"
490 for ebin in ebins : #header for energy bins
491 infostr += ("<%d" %ebin).rjust(7) #str(ebin).rjust(7)
492 infostr += "|"
493 print (infostr)
494 print ("-"*(12+len(ebins)*8)) #horizontal line
495 for etalow,etahigh in zip(etas,(etas[1:] + [self.maxeta])) : #looping over eta bins
496 prevebin = 0
497 erec = {}
498 egen = {}
499 hits = {}
500 count = {}
501 for ebin in ebins : # for all energy bins
502 count[ebin] = 0
503 erec[ebin] = 0.
504 egen[ebin] = 0.
505 hits[ebin] = 0.
506 for shower in self.library[etalow] :
507 if (shower.egen <= ebin) and (shower.egen > prevebin) :
508 count[ebin] += 1
509 egenshow = shower.egen
510 erecshow = 0
511 for hit in shower.shower :
512 erecshow += hit.e
513 erec[ebin] += erecshow
514 egen[ebin] += egenshow
515 hits[ebin] += len(shower.shower)
516 if (count[ebin] > 0) :
517 hits[ebin] /= count[ebin]
518 prevebin = ebin
519 infostr = "|#" # |
520 infostr+= str(round(etalow,5)).rjust(9) # | eta header
521 infostr+= "|" # |\
522 infostr2 = "|Hits"
523 infostr2+= str(round(etahigh,3)).rjust(6) # | eta header
524 infostr2+= "|" # |
525 infostr3 = "|ErecEgen"
526 infostr3+= " ".rjust(2) # | eta header
527 infostr3+= "|" # |
528 for ebin in ebins :
529 infostr+= str(count[ebin]).rjust(7) #print the number of showers
530 if (egen[ebin] > 0) :
531 infostr2+= ("%.2f" %(hits[ebin])).rjust(7)
532 infostr3+= ("%.5f" %(erec[ebin]/egen[ebin])).rjust(7)
533 else :
534 infostr2+= ("%.2f" %(hits[ebin])).rjust(7)
535 infostr3+= "0.0".rjust(7) #else print "xxx"
536 infostr+="|"
537 infostr2+="|"
538 infostr3+="|"
539 print (infostr)
540 print (infostr2)
541 print (infostr3)
542 print ("-"*(12+len(ebins)*8)) #horizontal line
543 def drawHits(self):
544 from ROOT import TH3F
545 from math import sqrt,copysign,log10
546 hits = TH3F("HITS","Hits Distrib",50,1,1000,101,-300,300,100,0,500)
547 containmentZ = TH3F("CONTZ","ContZ Distrib",50,1,1000,101,-300,300,100,0,500)
548 containmentR = TH3F("CONTR","ContR Distrib",50,1,1000,101,-300,300,100,0,500)
549 for etabin in self.library.values():
550 for storedShower in etabin :
551 containmentR.Fill(log10(storedShower.egen)*333,storedShower.rsize,storedShower.zsize,10)
552 containmentR.Fill(log10(storedShower.egen)*333,-storedShower.rsize,storedShower.zsize,10)
553 containmentZ.Fill(log10(storedShower.egen)*333,0,storedShower.zsize,10)
554 for hit in storedShower.shower :
555 hits.Fill(log10(storedShower.egen)*333,copysign(sqrt(hit.x*hit.x + hit.y*hit.y),hit.x),hit.z)
556 return hits,containmentZ,containmentR
557
558class FCALDistShowerLib() :
559 def __init__(self) :
560 self.library = {} # key (float) - dist, value (list) - list of StoredEnergyShower objs
561 self.detector= ""
562 self.particle= ""
563 self.release= ""
564 self.geometry= ""
565 self.geant= ""
566 self.phys= ""
567 self.comment= ""
568 self.xrod_cent = 0.0
569 self.yrod_cent = 0.0
570 self.step = 0.0
571 def scaleEnergy(self,scalefactor) :
572 for distbin in self.library.values():
573 for storedShower in distbin :
574 for hit in storedShower.shower :
575 hit.e *= scalefactor
576 self.comment += " SCALED: "+str(scalefactor)
577 def truncate(self,truncate) :
578 showers = []
579 for dist,distbin in self.library.items():
580 for storedShower in distbin :
581 showers += [(dist,storedShower)]
582 if len(showers) <= truncate :
583 print ("WARNING: Size of the library is already less:",truncate,"<",len(showers))
584 return
585 from random import randint
586 while (len(showers) > truncate) :
587 rand = randint(0,len(showers)-1)
588 self.library[showers[rand][0]].remove(showers[rand][1])
589 del showers[rand]
590 return
591 def moveDist(self,oldDist,newDist) :
592 if not (oldDist in self.library.keys()) :
593 return False
594 self.library[newDist] = self.library.pop(oldDist)
595 return True
596 def removeDist(self,dist) :
597 if not (dist in self.library.keys()) :
598 return False
599 self.library.pop(dist)
600 return True
601 def fromLibs(self,libs) :
602 for lib in libs :
603 if not isinstance(lib,self.__class__):
604 print ("ERROR: Different types of libs")
605 return False
606 self.detector = libs[0].detector
607 self.particle = libs[0].particle
608 self.release = libs[0].release
609 self.geometry = libs[0].geometry
610 self.geant = libs[0].geant
611 self.phys = libs[0].phys
612 self.comment = libs[0].comment
613 self.xrod_cent = libs[0].xrod_cent
614 self.yrod_cent = libs[0].yrod_cent
615 self.step = libs[0].step
616 dists = set(libs[0].library.keys())
617 for lib in libs :
618 if ( self.detector != lib.detector or
619 self.particle != lib.particle or
620 self.release != lib.release or
621 self.geometry != lib.geometry or
622 self.geant != lib.geant or
623 self.phys != lib.phys or
624 self.xrod_cent != lib.xrod_cent or
625 self.yrod_cent != lib.yrod_cent or
626 self.step != lib.step or
627 dists != set(lib.library.keys()) ) :
628 print ("ERROR: DIFFERENT LIBS!!!")
629 return False
630 for lib in libs :
631 for k,v in lib.library.items():
632 self.library.setdefault(k,set()).update(v)
633 for k,v in self.library.items():
634 self.library[k] = list(v)
635 return True
636 def readFromFile(self,filename) :
637 from ROOT import TFile
638 #from sets import Set
639 tfile = TFile(filename)
640 try:
641 ver = int(tfile.Get("version").GetVal())
642 except Exception:
643 print ("Not an FCALDistEnergyLib: Broken file")
644 tfile.Close()
645 return False
646 if (ver != 4) : #<<<<<<<<<<<<<<<<<<<<<<-------------- lib ver
647 print ("Not an FCALDistEnergyLib")
648 tfile.Close()
649 return False
650 meta = tfile.Get("meta")
651 libr = tfile.Get("library")
652
653 for event in meta :
654 self.detector=str(event.detector)
655 self.particle=str(event.particle)
656 self.release=str(event.release)
657 self.geometry=str(event.geometry)
658 self.geant=str(event.geantVersion)
659 self.phys=str(event.physicsList)
660 self.comment=str(event.comment)
661
662 state = -1
663 lastShower = False
664
665 for event in libr : #this is quite unclear, but easy to implement
666 if (state == -1) : #library header (calculator parameters)
667 self.xrod_cent = event.x
668 self.yrod_cent = event.y
669 self.step = event.z
670 state = 0
671 elif (state == 0) : #eta bin header
672 showersInCurDist = event.x
673 curDist = round(event.y,4)
674 self.library[curDist] = []
675 if (showersInCurDist > 0) :
676 state = 1 #go to shower header
677 elif (state == 1) : #shower header
678 hitsInCurShower = event.x
679 rSize = event.y
680 zSize = event.z
681 genEnergy = event.e
682 showersInCurDist -= 1
683 if (showersInCurDist == 0) : #last shower
684 lastShower = True
685 curShower = StoredEnergyShower()
686 curShower.egen = genEnergy
687 curShower.rsize = rSize
688 curShower.zsize = zSize
689 #curShower["hits"] = []
690 if (hitsInCurShower > 0) :
691 state = 2 #go to hits
692 else : #empty shower
693 self.library[curDist].append(curShower)
694 if (lastShower) : #special case of last shower in bin being the empty one
695 lastShower = False
696 state = 0 #next bin
697 elif (state == 2) :
698 hit = FourVector()
699 hit.e = event.e
700 hit.x = event.x
701 hit.y = event.y
702 hit.z = event.z
703 hit.time = event.time
704 curShower.shower.append(hit)
705 hitsInCurShower -= 1
706 if (hitsInCurShower == 0) : #last hit
707 self.library[curDist].append(curShower)
708 if (lastShower) : # end of eta bin
709 lastShower = False
710 state = 0
711 else : #not yet
712 state = 1
713 tfile.Close()
714 if (state != 0) :
715 print ("FILE CORRUPTED!!")
716 return False
717 return True
718 def writeToFile(self,filename) :
719 from ROOT import TFile,TTree,TParameter
720 from ROOT import gROOT, addressof
721 gROOT.ProcessLine(
722 "struct MyMetaStruct {\
723 Char_t detector[40];\
724 Char_t release[40];\
725 Char_t geometry[40];\
726 Char_t geant[40];\
727 Char_t phys[40];\
728 Char_t comment[400];\
729 Int_t particle;\
730 };" )
731 from ROOT import MyMetaStruct
732 gROOT.ProcessLine(
733 "struct MyStruct {\
734 Float_t x;\
735 Float_t y;\
736 Float_t z;\
737 Float_t e;\
738 Float_t time;\
739 };" )
740 from ROOT import MyStruct
741
742 tfile = TFile(filename,"RECREATE")
743
744 ver = TParameter(int)("version",4) #<<<<<<<<<<<<<<<<<<<<<<-------------- lib ver
745 ver.Write("version")
746
747 meta = TTree()
748 libr = TTree()
749
750 mmstruct = MyMetaStruct()
751
752 mmstruct.detector = "%s" % (str(self.detector))
753 mmstruct.particle = int(self.particle)
754 mmstruct.release = "%s" % (str(self.release))
755 mmstruct.geometry = "%s" % (str(self.geometry))
756 mmstruct.geant = "%s" % (str(self.geant))
757 mmstruct.phys = "%s" % (str(self.phys))
758 mmstruct.comment = "%s" % (str(self.comment))
759
760 meta.Branch("detector",addressof(mmstruct,"detector"),"detector/C")
761 meta.Branch("particle",addressof(mmstruct,"particle"),"particle/I")
762 meta.Branch("release",addressof(mmstruct,"release"),"release/C")
763 meta.Branch("geometry",addressof(mmstruct,"geometry"),"geometry/C")
764 meta.Branch("geantVersion",addressof(mmstruct,"geant"),"geantVersion/C")
765 meta.Branch("physicsList",addressof(mmstruct,"phys"),"physicsList/C")
766 meta.Branch("comment",addressof(mmstruct,"comment"),"physicsList/C")
767
768 meta.Fill()
769
770 mstruct = MyStruct()
771
772 libr.Branch("x",addressof(mstruct,"x"),"x/F")
773 libr.Branch("y",addressof(mstruct,"y"),"y/F")
774 libr.Branch("z",addressof(mstruct,"z"),"z/F")
775 libr.Branch("e",addressof(mstruct,"e"),"e/F")
776 libr.Branch("time",addressof(mstruct,"time"),"time/F")
777
778 mstruct.x = self.xrod_cent
779 mstruct.y = self.yrod_cent
780 mstruct.z = self.step
781 mstruct.e = 0
782 mstruct.time = 0
783 libr.Fill()
784
785 dists = sorted(self.library.keys())
786
787 for dist in dists :
788 mstruct.x = len(self.library[dist])
789 mstruct.y = dist
790 mstruct.z = 0
791 mstruct.e = 0
792 mstruct.time = 0
793 libr.Fill()
794 self.library[dist].sort(key=lambda x: x.egen)
795 for storedShower in self.library[dist] :
796 mstruct.x = len(storedShower.shower)
797 mstruct.y = storedShower.rsize
798 mstruct.z = storedShower.zsize
799 mstruct.e = storedShower.egen
800 mstruct.time = 0
801 libr.Fill()
802 for hit in storedShower.shower:
803 mstruct.e = hit.e
804 mstruct.x = hit.x
805 mstruct.y = hit.y
806 mstruct.z = hit.z
807 mstruct.time = hit.time
808 libr.Fill()
809 meta.Write("meta")
810 libr.Write("library")
811 tfile.Close()
812 def printInfo(self) :
813 print ("VERSION: FCALDistEnergyLib","PARTICLE:",self.particle,"DETECTOR:",self.detector)
814 print (self.release, self.geometry, self.geant, self.phys)
815 print ("xrodcent:",self.xrod_cent,"yrodcent:",self.yrod_cent,"step:",self.step)
816 print (self.comment)
817 ebins = [1,2,5,10,20,50,100,200,500,1000]
818 dists = sorted(self.library.keys())
819 print ("Number of etabins:",str(len(dists)))
820 fstot = 0
821 for dist in dists :
822 fstot +=len(self.library[dist])
823 print ("Number of showers:",str(fstot))
824 print ("-"*(13+len(ebins)*8)) #horizontal line
825 infostr = "|dists|ebins|"
826 for ebin in ebins : #header for energy bins
827 infostr += ("<%d" %ebin).rjust(7) #str(ebin).rjust(7)
828 infostr += "|"
829 print (infostr)
830 print ("-"*(13+len(ebins)*8)) #horizontal line
831 for distlow,disthigh in zip(dists,(dists[1:] + [4.5])) : #looping over eta bins
832 prevebin = 0
833 erec = {}
834 egen = {}
835 hits = {}
836 count = {}
837 for ebin in ebins : # for all energy bins
838 count[ebin] = 0
839 erec[ebin] = 0.
840 egen[ebin] = 0.
841 hits[ebin] = 0.
842 for shower in self.library[distlow] :
843 if (shower.egen <= ebin) and (shower.egen > prevebin) :
844 count[ebin] += 1
845 egenshow = shower.egen
846 erecshow = 0
847 for hit in shower.shower :
848 erecshow += hit.e
849 erec[ebin] += erecshow
850 egen[ebin] += egenshow
851 hits[ebin] += len(shower.shower)
852 if (count[ebin] > 0) :
853 hits[ebin] /= count[ebin]
854 prevebin = ebin
855 infostr = "|#" # |
856 infostr+= str(round(distlow,5)).rjust(10) # | eta header
857 infostr+= "|" # |\
858 infostr2 = "|Hits"
859 infostr2+= str(round(disthigh,3)).rjust(7) # | eta header
860 infostr2+= "|" # |
861 infostr3 = "|ErecEgen"
862 infostr3+= " ".rjust(3) # | eta header
863 infostr3+= "|" # |
864 for ebin in ebins :
865 infostr+= str(count[ebin]).rjust(7) #print the number of showers
866 if (egen[ebin] > 0) :
867 infostr2+= ("%.2f" %(hits[ebin])).rjust(7)
868 infostr3+= ("%.5f" %(erec[ebin]/egen[ebin])).rjust(7)
869 else :
870 infostr2+= ("%.2f" %(hits[ebin])).rjust(7)
871 infostr3+= "0.0".rjust(7) #else print "xxx"
872 infostr+="|"
873 infostr2+="|"
874 infostr3+="|"
875 print (infostr)
876 print (infostr2)
877 print (infostr3)
878 print ("-"*(12+len(ebins)*8)) #horizontal line
879 def drawHits(self):
880 from ROOT import TH3F
881 from math import sqrt,copysign,log10
882 hits = TH3F("HITS","Hits Distrib",50,1,1000,101,-300,300,100,0,500)
883 containmentZ = TH3F("CONTZ","ContZ Distrib",50,1,1000,101,-300,300,100,0,500)
884 containmentR = TH3F("CONTR","ContR Distrib",50,1,1000,101,-300,300,100,0,500)
885 for distbin in self.library.values():
886 for storedShower in distbin :
887 containmentR.Fill(log10(storedShower.egen)*333,storedShower.rsize,storedShower.zsize,10)
888 containmentR.Fill(log10(storedShower.egen)*333,-storedShower.rsize,storedShower.zsize,10)
889 containmentZ.Fill(log10(storedShower.egen)*333,0,storedShower.zsize,10)
890 for hit in storedShower.shower :
891 hits.Fill(log10(storedShower.egen)*333,copysign(sqrt(hit.x*hit.x + hit.y*hit.y),hit.x),hit.z)
892 return hits,containmentZ,containmentR
893
894class FCALDistEtaShowerLib() :
895
896 def __init__(self) :
897 self.library = {} # key (float) - eta, value (dict), key - dist, value - (list) list of StoredEnergyShower objs
898 self.detector= ""
899 self.particle= ""
900 self.release= ""
901 self.geometry= ""
902 self.geant= ""
903 self.phys= ""
904 self.comment= ""
905 self.xrod_cent = 0.0
906 self.yrod_cent = 0.0
907 self.step = 0.0
908 def scaleEnergy(self,scalefactor) :
909 for etabin in self.library.values():
910 for distbin in etabin.values():
911 for storedShower in distbin :
912 for hit in storedShower.shower :
913 hit.e *= scalefactor
914 self.comment += " SCALED: "+str(scalefactor)
915 def truncate(self,truncate,nShowersMin=0) :
916 log = logging.getLogger("FCALDistEtaShowerLib::truncate()")
917 showers = []
918 for eta,etabin in self.library.items():
919 for dist,distbin in etabin.items():
920 log.info("Number of showers in %s %s is %d",str(eta),str(dist),len(distbin))
921 for storedShower in distbin :
922 showers += [(eta, dist, storedShower)]
923 log.info("total number of showers: %d", len(showers))
924 if nShowersMin:
925 log.info("will not remove from eta-dist bins with less then %d showers", nShowersMin)
926 if len(showers) <= truncate :
927 log.warning("Size of the library is already less: %d < %d",truncate,len(showers))
928 return
929 from random import randint
930 while (len(showers) > truncate) :
931 rand = randint(0,len(showers)-1)
932 if len(self.library[showers[rand][0]][showers[rand][1]]) < nShowersMin:
933 continue
934 self.library[showers[rand][0]][showers[rand][1]].remove(showers[rand][2])
935 del showers[rand]
936 return
937 def moveDist(self,oldDist,newDist) :
938 rez = False
939 for eta,etabin in self.library.items():
940 if (oldDist in etabin.keys()) :
941 etabin[newDist] = etabin.pop(oldDist)
942 rez=True
943 return rez
944 def moveEta(self,oldEta,newEta) :
945 if not (oldEta in self.library.keys()) :
946 return False
947 self.library[newEta] = self.library.pop(oldEta)
948 return True
949 def removeDist(self,dist) :
950 rez = False
951 for eta,etabin in self.library.items():
952 if (dist in etabin.keys()) :
953 self.library.pop(dist)
954 rez=True
955 return rez
956 def removeEta(self,eta) :
957 if not (eta in self.library.keys()) :
958 return False
959 self.library.pop(eta)
960 return True
961 def fromLibs(self,libs) :
962 for lib in libs :
963 if not isinstance(lib,self.__class__):
964 print ("ERROR: Different types of libs")
965 return False
966 self.detector = libs[0].detector
967 self.particle = libs[0].particle
968 self.release = libs[0].release
969 self.geometry = libs[0].geometry
970 self.geant = libs[0].geant
971 self.phys = libs[0].phys
972 self.comment = libs[0].comment
973 self.xrod_cent = libs[0].xrod_cent
974 self.yrod_cent = libs[0].yrod_cent
975 self.step = libs[0].step
976 etas = set(libs[0].library.keys())
977 for lib in libs :
978 if ( self.detector != lib.detector or
979 self.particle != lib.particle or
980 self.release != lib.release or
981 self.geometry != lib.geometry or
982 self.geant != lib.geant or
983 self.phys != lib.phys or
984 self.xrod_cent != lib.xrod_cent or
985 self.yrod_cent != lib.yrod_cent or
986 self.step != lib.step or
987 etas != set(lib.library.keys()) ) :
988 print ("ERROR: DIFFERENT LIBS!!!")
989 return False
990 for eta in libs[0].library.keys() :
991 if (set(libs[0].library[eta].keys()) != set(lib.library[eta].keys())) :
992 print ("ERROR: DIFFERENT LIBS!!!")
993 return False
994 for lib in libs :
995 for k,v in lib.library.items():
996 for ki,vi in v.items():
997 self.library.setdefault(k,dict()).setdefault(ki,set()).update(vi)
998 for k,v in self.library.items():
999 for ki,vi in v.items():
1000 self.library[k][ki] = list(vi)
1001 return True
1002 def readFromFile(self,filename) :
1003 log = logging.getLogger("FCALDistEtaShowerLib::readFromFile()")
1004 from ROOT import TFile
1005 #from sets import Set
1006 tfile = TFile(filename)
1007 try:
1008 ver = int(tfile.Get("version").GetVal())
1009 except Exception:
1010 print ("Not an FCALDistEtaEnergyLib: Broken file")
1011 tfile.Close()
1012 return False
1013
1014 if (ver != 5) : #<<<<<<<<<<<<<<<<<<<<<<-------------- lib ver
1015 print ("Not an FCALDistEtaEnergyLib")
1016 tfile.Close()
1017 return False
1018 meta = tfile.Get("meta")
1019 libr = tfile.Get("library")
1020
1021 for event in meta :
1022 self.detector=str(event.detector)
1023 self.particle=str(event.particle)
1024 self.release=str(event.release)
1025 self.geometry=str(event.geometry)
1026 self.geant=str(event.geantVersion)
1027 self.phys=str(event.physicsList)
1028 self.comment=str(event.comment)
1029
1030 state = -1
1031 lastShower = False
1032 lastEta = False
1033
1034 log.debug("dector: %s", str(event.detector))
1035 log.debug("particle: %s", str(event.particle))
1036 log.debug("release: %s", str(event.release))
1037 log.debug("geometry: %s", str(event.geometry))
1038 log.debug("geant ver: %s", str(event.geantVersion))
1039 log.debug("physList: %s", str(event.physicsList))
1040 log.debug("comment: %s", str(event.comment))
1041
1042 for event in libr : #this is quite unclear, but easy to implement, we change the "state" depending on what we are reading
1043 log.debug("-------")
1044 log.debug("x=%f, y=%f, z=%f, e=%f",event.x,event.y,event.z,event.e)
1045 log.debug("beginnnig ev loop. lastShower: %s",str(lastShower))
1046 log.debug("beginnnig ev loop. state: %s",str(state))
1047
1048 if (state == -1) : #library header (calculator parameters)
1049 log.debug("in state=-1")
1050
1051 self.xrod_cent = event.x
1052 self.yrod_cent = event.y
1053 self.step = event.z
1054 state = 0
1055 elif (state == 0) : #eta bin header
1056 log.debug("in state=0")
1057 log.debug("x=distsInCurEta, y=curEta")
1058
1059 distsInCurEta = event.x
1060 curEta = round(event.y,4)
1061 self.library[curEta] = {}
1062 if (distsInCurEta > 0) :
1063 state = 1 #go to dist header
1064 elif (state == 1) :
1065 log.debug("in state=1")
1066 log.debug("x=showersInCurDist, y=curDist")
1067
1068 showersInCurDist = event.x
1069 curDist = round(event.y,4)
1070 self.library[curEta][curDist] = []
1071 distsInCurEta -= 1
1072 if (distsInCurEta == 0) :
1073 lastEta = True
1074 if (showersInCurDist > 0) :
1075 state = 2 #go to shower header
1076 else : #empty dist bin
1077 if (lastEta) : #special case of last eta bin being the empty one
1078 lastEta = False
1079 state = 0
1080 elif (state == 2) :
1081 # writing shower info
1082 log.debug("in state=2")
1083 log.debug("x=hitsInCurShower, y=curShower.rSize, z=curShower.zSize, e=curShower.genEnergy")
1084
1085 hitsInCurShower = event.x
1086 rSize = event.y
1087 zSize = event.z
1088 genEnergy = event.e
1089 showersInCurDist -= 1
1090 if (showersInCurDist == 0) : #last shower
1091 lastShower = True
1092 curShower = StoredEnergyShower()
1093 curShower.egen = genEnergy
1094 curShower.rsize = rSize
1095 curShower.zsize = zSize
1096 #curShower["hits"] = []
1097 if (hitsInCurShower > 0) :
1098 state = 3 #go to hits
1099 else : #empty shower
1100 log.debug("Appending shower to lib pos %s %s",curEta,curDist)
1101
1102 self.library[curEta][curDist].append(curShower)
1103 if (lastShower) : #special case of last shower in bin being the empty one
1104 lastShower = False
1105 if (lastEta) : #double special case: last shower in last eta bin is empty
1106 lastEta = False
1107 state = 0 #next eta bin
1108 else :
1109 state = 1 #next dist bin
1110 elif (state == 3) :
1111
1112 log.debug("in state=3")
1113 log.debug("x=hit.x, y=hit.y, z=hit.z, e=hit.e")
1114
1115 hit = FourVector()
1116 hit.e = event.e
1117 hit.x = event.x
1118 hit.y = event.y
1119 hit.z = event.z
1120 hit.time = event.time
1121 curShower.shower.append(hit)
1122 hitsInCurShower -= 1
1123 if (hitsInCurShower == 0) : #last hit
1124 log.debug("Appending shower+hit to lib pos %s %s",curEta,curDist)
1125
1126 self.library[curEta][curDist].append(curShower)
1127 if (lastShower) : # end of dist bin
1128 lastShower = False
1129 if (lastEta) : #end of eta bin as well
1130 lastEta = False
1131 state = 0 #next eta bin
1132 else :
1133 state = 1 #next dist bin
1134 else : #not yet
1135 state = 2
1136
1137 log.debug("ending ev loop. lastShower: %s", lastShower)
1138 log.debug("ending ev loop. state %s", state)
1139 if log.root.level == logging.DEBUG:
1140 input("Continue? Press Enter.")
1141
1142 tfile.Close()
1143 if (state != 0) : #the last entry should be the last hit of the last shower in the last bin. if not - file is corrupted
1144 print ("FILE CORRUPTED!!")
1145 return False
1146 return True
1147 def writeToFile(self,filename) :
1148 from ROOT import TFile,TTree,TParameter
1149 from ROOT import gROOT, addressof
1150 gROOT.ProcessLine(
1151 "struct MyMetaStruct {\
1152 Char_t detector[40];\
1153 Char_t release[40];\
1154 Char_t geometry[40];\
1155 Char_t geant[40];\
1156 Char_t phys[40];\
1157 Char_t comment[400];\
1158 Int_t particle;\
1159 };" )
1160 from ROOT import MyMetaStruct
1161 gROOT.ProcessLine(
1162 "struct MyStruct {\
1163 Float_t x;\
1164 Float_t y;\
1165 Float_t z;\
1166 Float_t e;\
1167 Float_t time;\
1168 };" )
1169 from ROOT import MyStruct
1170
1171 tfile = TFile(filename,"RECREATE")
1172
1173 ver = TParameter(int)("version",5) #<<<<<<<<<<<<<<<<<<<<<<-------------- lib ver
1174 ver.Write("version")
1175
1176 meta = TTree()
1177 libr = TTree()
1178
1179 mmstruct = MyMetaStruct()
1180
1181 mmstruct.detector = "%s" % (str(self.detector))
1182 mmstruct.particle = int(self.particle)
1183 mmstruct.release = "%s" % (str(self.release))
1184 mmstruct.geometry = "%s" % (str(self.geometry))
1185 mmstruct.geant = "%s" % (str(self.geant))
1186 mmstruct.phys = "%s" % (str(self.phys))
1187 mmstruct.comment = "%s" % (str(self.comment))
1188
1189 meta.Branch("detector",addressof(mmstruct,"detector"),"detector/C")
1190 meta.Branch("particle",addressof(mmstruct,"particle"),"particle/I")
1191 meta.Branch("release",addressof(mmstruct,"release"),"release/C")
1192 meta.Branch("geometry",addressof(mmstruct,"geometry"),"geometry/C")
1193 meta.Branch("geantVersion",addressof(mmstruct,"geant"),"geantVersion/C")
1194 meta.Branch("physicsList",addressof(mmstruct,"phys"),"physicsList/C")
1195 meta.Branch("comment",addressof(mmstruct,"comment"),"physicsList/C")
1196
1197 meta.Fill()
1198
1199 mstruct = MyStruct()
1200
1201 libr.Branch("x",addressof(mstruct,"x"),"x/F")
1202 libr.Branch("y",addressof(mstruct,"y"),"y/F")
1203 libr.Branch("z",addressof(mstruct,"z"),"z/F")
1204 libr.Branch("e",addressof(mstruct,"e"),"e/F")
1205 libr.Branch("time",addressof(mstruct,"time"),"time/F")
1206
1207 etas = sorted(self.library.keys())
1208
1209 mstruct.x = self.xrod_cent
1210 mstruct.y = self.yrod_cent
1211 mstruct.z = self.step
1212 mstruct.e = 0
1213 mstruct.time = 0
1214 libr.Fill()
1215
1216 for eta in etas :
1217 dists = sorted(self.library[eta].keys())
1218 mstruct.x = len(self.library[eta])
1219 mstruct.y = eta
1220 mstruct.z = 0
1221 mstruct.e = 0
1222 mstruct.time = 0
1223 libr.Fill()
1224 for dist in dists :
1225 mstruct.x = len(self.library[eta][dist])
1226 mstruct.y = dist
1227 mstruct.z = 0
1228 mstruct.e = 0
1229 mstruct.time = 0
1230 libr.Fill()
1231 self.library[eta][dist].sort(key=lambda x: x.egen)
1232 for storedShower in self.library[eta][dist] :
1233 mstruct.x = len(storedShower.shower)
1234 mstruct.y = storedShower.rsize
1235 mstruct.z = storedShower.zsize
1236 mstruct.e = storedShower.egen
1237 mstruct.time = 0
1238 libr.Fill()
1239 for hit in storedShower.shower:
1240 mstruct.e = hit.e
1241 mstruct.x = hit.x
1242 mstruct.y = hit.y
1243 mstruct.z = hit.z
1244 mstruct.time = hit.time
1245 libr.Fill()
1246 meta.Write("meta")
1247 libr.Write("library")
1248 tfile.Close()
1249 def printInfo(self) :
1250 print ("VERSION: FCALDistEtaEnergyLib","PARTICLE:",self.particle,"DETECTOR:",self.detector)
1251 print (self.release, self.geometry, self.geant, self.phys)
1252 print ("xrodcent:",self.xrod_cent,"yrodcent:",self.yrod_cent,"step:",self.step)
1253 print (self.comment)
1254 ebins = [1,2,3,4,5,10,20,50,100,200,500,1000]
1255 etas = sorted(self.library.keys())
1256 print ("Number of etabins:",str(len(etas)))
1257 fstot = 0
1258 for etabin in self.library.values():
1259 for distbin in etabin.values():
1260 fstot +=len(distbin)
1261 print ("Number of showers:",str(fstot))
1262 print ("-"*(13+len(ebins)*8)) #horizontal line
1263 infostr = "|dists|ebins|"
1264 for ebin in ebins : #header for energy bins
1265 infostr += ("<%d" %ebin).rjust(7) #str(ebin).rjust(7)
1266 infostr += "|"
1267 print (infostr)
1268 print ("-"*(13+len(ebins)*8)) #horizontal line
1269 for eta in etas :
1270 dists = sorted(self.library[eta].keys())
1271 for distlow,disthigh in zip(dists,(dists[1:] + [4.5])) : #looping over eta bins
1272 prevebin = 0
1273 erec = {}
1274 egen = {}
1275 hits = {}
1276 count = {}
1277 for ebin in ebins : # for all energy bins
1278 count[ebin] = 0
1279 erec[ebin] = 0.
1280 egen[ebin] = 0.
1281 hits[ebin] = 0.
1282 for shower in self.library[eta][distlow] :
1283 if (shower.egen <= ebin) and (shower.egen > prevebin) :
1284 count[ebin] += 1
1285 egenshow = shower.egen
1286 erecshow = 0
1287 for hit in shower.shower :
1288 erecshow += hit.e
1289 erec[ebin] += erecshow
1290 egen[ebin] += egenshow
1291 hits[ebin] += len(shower.shower)
1292 if (count[ebin] > 0) :
1293 hits[ebin] /= count[ebin]
1294 prevebin = ebin
1295 infostr = "|#" # |
1296 infostr+= str(eta).rjust(10) # | eta header
1297 infostr+= "|" # |\
1298 infostr2 = "|Hits"
1299 infostr2+= str(round(distlow,5)).rjust(7) # | eta header
1300 infostr2+= "|" # |
1301 infostr3 = "|E/E"
1302 infostr3+= str(round(disthigh,3)).rjust(8) # | eta header
1303 infostr3+= "|" # |
1304 for ebin in ebins :
1305 infostr+= str(count[ebin]).rjust(7) #print the number of showers
1306 if (egen[ebin] > 0) :
1307 infostr2+= ("%.2f" %(hits[ebin])).rjust(7)
1308 infostr3+= ("%.5f" %(erec[ebin]/egen[ebin])).rjust(7)
1309 else :
1310 infostr2+= ("%.2f" %(hits[ebin])).rjust(7)
1311 infostr3+= "0.0".rjust(7) #else print "xxx"
1312 infostr+="|"
1313 infostr2+="|"
1314 infostr3+="|"
1315 print (infostr)
1316 print (infostr2)
1317 print (infostr3)
1318 print ("-"*(12+len(ebins)*8)) #horizontal line
1319 def drawHits(self):
1320 from ROOT import TH3F
1321 from math import sqrt,copysign,log10
1322 hits = TH3F("HITS","Hits Distrib",50,1,1000,101,-300,300,100,0,500)
1323 containmentZ = TH3F("CONTZ","ContZ Distrib",50,1,1000,101,-300,300,100,0,500)
1324 containmentR = TH3F("CONTR","ContR Distrib",50,1,1000,101,-300,300,100,0,500)
1325
1326 etas = sorted(self.library.keys())
1327 for eta in etas :
1328 dists = sorted(self.library[eta].keys())
1329 for distlow,disthigh in zip(dists,(dists[1:] + [4.5])) : #looping over eta bins
1330 for storedShower in self.library[eta][distlow] :
1331 containmentR.Fill(log10(storedShower.egen)*333,storedShower.rsize,storedShower.zsize,10)
1332 containmentR.Fill(log10(storedShower.egen)*333,-storedShower.rsize,storedShower.zsize,10)
1333 containmentZ.Fill(log10(storedShower.egen)*333,0,storedShower.zsize,10)
1334 for hit in storedShower.shower :
1335 hits.Fill(log10(storedShower.egen)*333,copysign(sqrt(hit.x*hit.x + hit.y*hit.y),hit.x),hit.z)
1336
1337 return hits,containmentZ,containmentR
sort
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
Definition DVL_algorithms.h:554
set
STL class.
Generated on for ATLAS Offline Software by doxygen 1.14.0