ATLAS Offline Software
Loading...
Searching...
No Matches
LArMonTransforms.py
Go to the documentation of this file.
2# Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3#
4#
5# Author: Margherita Spalla (margherita.spalla@cern.ch)
6# Function for post-processing of histograms from LArMonitoring
7#
8# The functions in this file are being implemented through the yaml config file inside the DataQuality package:
9# DataQuality/DataQualityUtils/data/postprocessing/LArMonPostProc.yaml
10#
11
12
13from ROOT import TMath, TH1F, TH2F, TProfile
14import cppyy
15
16def setMaxMin(inputs,maxVal=0,minVal=0,useMax=True,useMin=True):
17 """ This function sets the histogram max/min """
18 assert len(inputs) == 1 , len(inputs) # for testing
19 cl = inputs[0][1][0].Clone()
20
21 if useMax:
22 cl.SetMaximum(maxVal)
23 pass
24
25 if useMin:
26 cl.SetMinimum(minVal)
27 pass
28
29 return [cl]
30
31cppyy.cppdef("""
32void takemax(int sc, const TH2I* inhist, TH2I* outhist) {
33 for (int i = 0; i < outhist->GetSize(); ++i) {
34 if (inhist->At(i) != 0 && sc > outhist->At(i)) {
35 outhist->SetBinContent(i, sc);
36 }
37 }
38}
39""")
40
41cppyy.cppdef("""
42void clearhist(TH2I* hist) {
43 for (int i = 0; i < hist->GetSize(); ++i) {
44 hist->SetBinContent(i,0);
45 hist->SetBinError(i,0);
46 }
47}
48""")
49# @profile
50def fillWithMaxCoverage(inputs,isFtSlotPlot=True):
51 """ For each bin, fill the output with the max filled error code. All histograms should have the same bin content"""
52 inputs = list(inputs) # evaluate this just once
53 assert len(inputs) >= 1 , len(inputs) # for testing
54
55 cl = inputs[0][1][0].Clone()
56 partition = inputs[0][0]['part']
57 if 'side' in inputs[0][0]:
58 partition=partition+inputs[0][0]['side']
59 pass
60 if 'sampling' in inputs[0][0]:
61 sampling = int(inputs[0][0]['sampling'])
62 else:
63 sampling=-1 #a default value, != 0,1,2,3
64 pass
65
66 doFixEmptyBins = (not isFtSlotPlot and (('EMB' in partition and sampling==1) or ('EMEC' in partition and sampling==2) or ('HEC' in partition)))
67
68 #if we only have one and doesn't need bin fix, return it
69 if len(inputs)==1 and not doFixEmptyBins:
70 cl.SetMaximum(4.4)
71 cl.SetMinimum(1)
72 return [cl]
73
74 if len(inputs)>1 :
75 #check for incompatible histograms
76 Nx=cl.GetNbinsX()
77 Ny=cl.GetNbinsY()
78 binsok=True
79 for j in range(1,len(inputs)):
80 if inputs[j][1][0].GetNbinsX()!=Nx:
81 binsok=False
82 break
83 if inputs[j][1][0].GetNbinsY()!=Ny:
84 binsok=False
85 break
86 pass
87 assert binsok , "Histograms must all have the same number of bins"
88
89 #fill with max
90 cl.Clear()
91 cppyy.gbl.clearhist(cl)
92 for i in inputs:
93 statusCode=int(i[0]['sc'])
94 h = i[1][0]
95 cppyy.gbl.takemax(statusCode, h, cl)
96
97
98 cl.SetMaximum(4.4)
99 cl.SetMinimum(1)
100
101 if(not doFixEmptyBins):
102 return [cl]
103
104 #only for some eta-phi plots: fix empty bins
105 if 'EMBA' in partition:
106 ieta_min = 1
107 ieta_max = 447
108 iphi_min = 1
109 iphi_max = 256
110 iphi_step = 4
111 Diphi_toGet = 1
112 Diphi_toSet = [0,2,3]
113 elif 'EMBC' in partition:
114 ieta_min = 4
115 ieta_max = 450
116 iphi_min = 1
117 iphi_max = 256
118 iphi_step = 4
119 Diphi_toGet = 2
120 Diphi_toSet = [0,1,3]
121 elif 'EMECA' in partition:
122 ieta_min = 44
123 ieta_max = 51
124 iphi_min = 1
125 iphi_max = 256
126 iphi_step = 4
127 Diphi_toGet = 1
128 Diphi_toSet = [0,2,3]
129 elif 'EMECC' in partition:
130 ieta_min = 1
131 ieta_max = 7
132 iphi_min = 1
133 iphi_max = 256
134 iphi_step = 4
135 Diphi_toGet = 2
136 Diphi_toSet = [0,1,3]
137 else: #HEC
138 hecetabinmin=[11,11,9,9]
139 hecetabinmax=[14,13,12,12]
140 if partition=='HECA':
141 ieta_min = hecetabinmin[sampling]
142 ieta_max = hecetabinmax[sampling]
143 else:
144 ieta_min = 1
145 ieta_max = hecetabinmax[sampling]-hecetabinmin[sampling]
146 pass
147 iphi_min = 1
148 iphi_max = 64
149 iphi_step = 2
150 Diphi_toGet = 1
151 Diphi_toSet = [0]
152 pass
153
154 iphi=iphi_min
155 while iphi<=iphi_max:
156 ieta=ieta_min
157 while ieta<=ieta_max:
158 content = cl.GetBinContent(ieta,iphi+Diphi_toGet)
159 error = cl.GetBinError(ieta,iphi+Diphi_toGet)
160 for Di in Diphi_toSet:
161 cl.SetBinContent(ieta,iphi+Di,content)
162 cl.SetBinError(ieta,iphi+Di,error)
163 pass
164 ieta=ieta+1
165 pass
166 iphi=iphi+iphi_step
167 pass
168
169 return [cl]
170
171
172
173
174
175def normToEntries(inputs,titleToReplace="",replaceWith=""):
176 """ This function creates TH2F as ratio of input TH2I and the number of entries of histogram 2 (which is supposed to represent the number of events) """
177 assert len(inputs) == 1 , len(inputs)
178 assert len(inputs[0][1]) == 2 , len(inputs[0][1])
179
180 #cl = inputs[0][1][0].Clone()
181 inh=inputs[0][1][0]
182 cl=TH2F(inh.GetName(),inh.GetTitle(),inh.GetXaxis().GetNbins(),inh.GetXaxis().GetBinLowEdge(1),inh.GetXaxis().GetBinUpEdge(inh.GetXaxis().GetNbins()), inh.GetYaxis().GetNbins(),inh.GetYaxis().GetBinLowEdge(1),inh.GetYaxis().GetBinUpEdge(inh.GetYaxis().GetNbins()))
183 clx=cl.GetXaxis()
184 clx.SetTitle(inh.GetXaxis().GetTitle())
185 cl.GetYaxis().SetTitle(inh.GetYaxis().GetTitle())
186
187 Nen = inputs[0][1][1].GetEntries()
188 if Nen!=0:
189 for xbins in range(1,cl.GetNbinsX()+1):
190 for ybins in range(1,cl.GetNbinsY()+1):
191 ibin=inh.GetBin(xbins,ybins)
192 cl.SetBinContent(ibin, 100.*inh.GetBinContent(ibin) / Nen)
193 pass
194
195 if titleToReplace=="":
196 return [cl]
197
198 tit = cl.GetTitle()
199 tit=tit.replace(titleToReplace,replaceWith)
200 cl.SetTitle(tit)
201
202 return [cl]
203
204
205
206def normToEntriesAndSetMin(inputs,minVal=0,maxVal=0,useMax=False,clone=True,titleToReplace="",replaceWith=""):
207 """ This function normalises histogram 1 to the number of entries of histogram 2 (which is supposed to represent the number of events) and sets the histogram max/min """
208
209 assert len(inputs) == 1 , len(inputs)
210 assert len(inputs[0][1]) == 2 , len(inputs[0][1])
211
212 if clone:
213 cl = inputs[0][1][0].Clone()
214 else:
215 inh=inputs[0][1][0]
216 cl=TH2F(inh.GetName()[3:],inh.GetTitle(),inh.GetXaxis().GetNbins(),inh.GetXaxis().GetBinLowEdge(1),inh.GetXaxis().GetBinUpEdge(inh.GetXaxis().GetNbins()), inh.GetYaxis().GetNbins(),inh.GetYaxis().GetBinLowEdge(1),inh.GetYaxis().GetBinUpEdge(inh.GetYaxis().GetNbins()))
217 cl.GetXaxis().SetTitle(inh.GetXaxis().GetTitle())
218 cl.GetYaxis().SetTitle(inh.GetYaxis().GetTitle())
219 for xbins in range(1,cl.GetNbinsX()+1):
220 for ybins in range(1,cl.GetNbinsY()+1):
221 ibin=inh.GetBin(xbins,ybins)
222 cl.SetBinContent(ibin, inh.GetBinContent(ibin))
223
224 Nen = inputs[0][1][1].GetEntries()
225 if Nen!=0:
226 cl.Scale(100./Nen)
227 pass
228
229 cl.SetMinimum(minVal)
230 if useMax:
231 cl.SetMaximum(maxVal)
232 pass
233
234 if titleToReplace=="":
235 return [cl]
236
237 tit = cl.GetTitle()
238 tit=tit.replace(titleToReplace,replaceWith)
239 cl.SetTitle(tit)
240
241 return [cl]
242
243def normToBinAndSetMinMax(inputs,bin_norm=0,minVal=0,maxVal=0,useMax=False, titleToReplace="",replaceWith="",newYaxis=""):
244 """ This function normalises histogram 1 to the content of bin bin_norm (which is supposed to represent the number of events) and sets the histogram max/min """
245
246 assert len(inputs) == 1 , len(inputs)
247 assert len(inputs[0][1]) == 1 , len(inputs[0][1])
248
249 inh=inputs[0][1][0]
250 cl=TH1F(inh.GetName()+"Yield",inh.GetTitle(),inh.GetNbinsX(),inh.GetBinLowEdge(1),inh.GetBinLowEdge(inh.GetNbinsX()))
251 clx=cl.GetXaxis()
252 clx.SetTitle(inh.GetXaxis().GetTitle())
253 for i in range(1,inh.GetNbinsX()+1):
254 clx.SetBinLabel(i,inh.GetXaxis().GetBinLabel(i))
255 cl.GetYaxis().SetTitle(inh.GetYaxis().GetTitle())
256
257 Nen = inputs[0][1][0].GetBinContent(bin_norm)
258 if Nen!=0:
259 for ibin in range(0,inh.GetNbinsX()+1):
260 cl.SetBinContent(ibin,100.*inh.GetBinContent(ibin)/Nen)
261 pass
262
263 cl.SetMinimum(minVal)
264 if useMax:
265 cl.SetMaximum(maxVal)
266 pass
267
268 if newYaxis!="":
269 cl.GetYaxis().SetTitle(newYaxis)
270
271 if titleToReplace=="":
272 return [cl]
273
274 tit = cl.GetTitle()
275 tit=tit.replace(titleToReplace,replaceWith)
276 cl.SetTitle(tit)
277
278 return [cl]
279
280
281
282def divideHist(inputs,titleToReplace="",replaceWith=""):
283 """ This function create a new TH1F from ratio of two ROOT histograms """
284 assert len(inputs) == 1
285 assert len(inputs[0][1]) == 2
286
287 inh=inputs[0][1][0]
288 cl=TH1F(inh.GetName(),inh.GetTitle(),inh.GetNbinsX(),inh.GetBinLowEdge(1),inh.GetBinLowEdge(inh.GetNbinsX()))
289 clx=cl.GetXaxis()
290 clx.SetTitle(inh.GetXaxis().GetTitle())
291 cl.GetYaxis().SetTitle(inh.GetYaxis().GetTitle())
292 inn=inputs[0][1][1]
293 xbins = inh.GetNbinsX()
294 if inn.GetNbinsX() < xbins: xbins=inn.GetNbinsX()
295 for ibin in range(1,xbins+1):
296 if inn.GetBinContent(ibin) > 0:
297 cl.SetBinContent(ibin,1.*inh.GetBinContent(ibin)/inn.GetBinContent(ibin))
298 else:
299 cl.SetBinContent(ibin,1.*inh.GetBinContent(ibin))
300
301 if titleToReplace=="":
302 return [cl]
303
304 tit = cl.GetTitle()
305 tit=tit.replace(titleToReplace,replaceWith)
306 cl.SetTitle(tit)
307
308 return [cl]
309
310
311def digitSummary(inputs,TreshOut=5,TreshSat=5,TreshNull=5):
312
313 cl = TH2F("summary","LArDigit Summary;Status;Partition",4,0.,4.,8,0.,8.)
314 cl.SetDirectory(0)
315 cl.GetYaxis().SetBinLabel(1,"EMBC")
316 cl.GetYaxis().SetBinLabel(2,"EMBA")
317 cl.GetYaxis().SetBinLabel(3,"EMECC")
318 cl.GetYaxis().SetBinLabel(4,"EMECA")
319 cl.GetYaxis().SetBinLabel(5,"HECC")
320 cl.GetYaxis().SetBinLabel(6,"HECA")
321 cl.GetYaxis().SetBinLabel(7,"FCalC")
322 cl.GetYaxis().SetBinLabel(8,"FCalA")
323 cl.GetXaxis().SetBinLabel(1,"OutOfRange")
324 cl.GetXaxis().SetBinLabel(2,"Saturation")
325 cl.GetXaxis().SetBinLabel(3,"Null Digits")
326 cl.GetXaxis().SetBinLabel(4,"Mean Time")
327 cl.GetXaxis().SetLabelSize(0.055)
328 cl.GetYaxis().SetLabelSize(0.055)
329
330 for i in inputs:
331 assert len(i[1])==1, len(i[1])
332
333 if i[1][0].GetEntries()==0:
334 continue
335
336 #y bin
337 pr=i[0]['part']
338 if "EMB" in pr:
339 iy=1
340 elif "EMEC" in pr:
341 iy=3
342 elif "HEC" in pr:
343 iy=5
344 elif "FCal" in pr:
345 iy=7
346 else:
347 iy=-10
348 pass
349
350 if i[0]['side'] == "A":
351 iy=iy+1
352 pass
353
354 #x bin
355 hn=i[0]['histname']
356 if hn=="EnVsTime":
357 ix=4
358 content = 0
359 for xbin in range(1,i[1][0].GetNbinsX()+1):
360 sumx=0
361 for ybin in range(1,i[1][0].GetNbinsY()+1):
362 sumx=sumx+i[1][0].GetBinContent(xbin,ybin)
363 pass
364 content=content+(sumx*i[1][0].GetXaxis().GetBinLowEdge(xbin))
365 pass
366 content=content/float(i[1][0].GetEntries())
367 else:
368 if "OutOfRange" in hn:
369 ix=1
370 treshold=TreshOut
371 elif "Saturation" in hn:
372 ix=2
373 treshold=TreshSat
374 elif "NullDigit" in hn:
375 ix=3
376 treshold=TreshNull
377 else:
378 ix=-1
379 treshold=1e4
380 pass
381
382 content=0
383 for xbin in range(i[1][0].GetNbinsX()):
384 for ybin in range(i[1][0].GetNbinsY()):
385 if i[1][0].GetBinContent(xbin,ybin)>treshold :
386 content=content+1
387 pass
388 pass
389 pass
390 pass
391
392 if ix>0 and iy>0:
393 cl.SetBinContent(ix,iy,content)
394 pass
395
396 pass #end of input loop
397
398 return [cl]
399
400
402 #check that the inputs are ok
403 assert len(inputs)==2 , len(inputs)
404
405 assert ((inputs[0][0]['name']=='average' and inputs[1][0]['name']=='partialSum') or (inputs[1][0]['name']=='average' and inputs[0][0]['name']=='partialSum')) , "Wrong inputs"
406
407 if inputs[1][0]['name']=='partialSum':
408 i_parSum=1
409 i_av=0
410 else:
411 i_parSum=0
412 i_av=1
413
414 inh=inputs[i_parSum][1][0]
415 cl=TH2F(inh.GetName()[3:],inh.GetTitle()[2:12],inh.GetXaxis().GetNbins(),inh.GetXaxis().GetBinLowEdge(1),inh.GetXaxis().GetBinUpEdge(inh.GetXaxis().GetNbins()), inh.GetYaxis().GetNbins(),inh.GetYaxis().GetBinLowEdge(1),inh.GetYaxis().GetBinUpEdge(inh.GetYaxis().GetNbins()))
416 cl.GetXaxis().SetTitle(inh.GetXaxis().GetTitle())
417 cl.GetYaxis().SetTitle(inh.GetYaxis().GetTitle())
418
419 for i in range(1,cl.GetNbinsX()+1):
420 mean1 = inputs[i_av][1][0].GetBinContent(i)
421 sumVar1 = inputs[i_parSum][1][0].GetBinContent(i,i)
422 N = inputs[i_av][1][0].GetBinEntries(i)
423 if N==0:
424 continue
425 for j in range(i+1,cl.GetNbinsY()+1):
426 mean2 = inputs[i_av][1][0].GetBinContent(j)
427 sumVar2 = inputs[i_parSum][1][0].GetBinContent(j,j)
428 if (sumVar1-N*mean1*mean1)*(sumVar2-N*mean2*mean2)==0:
429 continue
430 cor=(inputs[i_parSum][1][0].GetBinContent(i,j)-N*mean1*mean2)/TMath.Sqrt((sumVar1-N*mean1*mean1)*(sumVar2-N*mean2*mean2))
431 cl.SetBinContent(i,j,cor)
432 cl.SetBinContent(j,i,cor)
433
434
435 return [cl]
436
437
438def Mean(inputs):
439
440 cl = TProfile("summary","Summary;Partition;Mean",8,0.,8.,0.,10., "s")
441 cl.SetDirectory(0)
442 cl.GetXaxis().SetBinLabel(1,"EMBA")
443 cl.GetXaxis().SetBinLabel(2,"EMBC")
444 cl.GetXaxis().SetBinLabel(3,"EMECA")
445 cl.GetXaxis().SetBinLabel(4,"EMECC")
446 cl.GetXaxis().SetBinLabel(5,"HECA")
447 cl.GetXaxis().SetBinLabel(6,"HECC")
448 cl.GetXaxis().SetBinLabel(7,"FCalA")
449 cl.GetXaxis().SetBinLabel(8,"FCalC")
450 cl.GetXaxis().SetLabelSize(0.055)
451 cl.GetYaxis().SetLabelSize(0.055)
452
453 for i in inputs:
454 assert len(i[1])==1, len(i[1])
455
456 #x bin
457 pr=i[0]['part']
458 if "EMB" in pr:
459 ix=1
460 elif "EMEC" in pr:
461 ix=3
462 elif "HEC" in pr:
463 ix=5
464 elif "FCal" in pr:
465 ix=7
466 else:
467 ix=-10
468 pass
469 if i[0]['side'] == "C":
470 ix=ix+1
471 pass
472
473 #y bin
474 if i[1][0].GetEntries()==0:
475 continue
476 else:
477 if ix>=0:
478 #filling the profile histogram for each x_bin separately
479 for xbin in range(0,i[1][0].GetNbinsX()):
480 #cl.Fill (ix-1, i[1][0].GetBinContent(xbin)) #for all bins
481 if i[1][0].GetBinContent(xbin) > 0.:
482 cl.Fill (ix-1, i[1][0].GetBinContent(xbin)) #only for non-zero bins
483 pass #end of input loop
484
485 return [cl]
if(febId1==febId2)
TGraphErrors * GetEntries(TH2F *histo)
normToBinAndSetMinMax(inputs, bin_norm=0, minVal=0, maxVal=0, useMax=False, titleToReplace="", replaceWith="", newYaxis="")
normToEntries(inputs, titleToReplace="", replaceWith="")
normToEntriesAndSetMin(inputs, minVal=0, maxVal=0, useMax=False, clone=True, titleToReplace="", replaceWith="")
fillWithMaxCoverage(inputs, isFtSlotPlot=True)
divideHist(inputs, titleToReplace="", replaceWith="")
setMaxMin(inputs, maxVal=0, minVal=0, useMax=True, useMin=True)
digitSummary(inputs, TreshOut=5, TreshSat=5, TreshNull=5)