ATLAS Offline Software
LArMonTransforms.py
Go to the documentation of this file.
1 #
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 
13 from ROOT import TMath, TH1F, TH2F, TProfile
14 import cppyy
15 
16 def 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 
31 cppyy.cppdef("""
32 void 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 
41 cppyy.cppdef("""
42 void 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
50 def 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 
175 def 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 
206 def 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 
243 def 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 
282 def 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 
311 def 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 
438 def 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]
LArMonTransforms.setMaxMin
def setMaxMin(inputs, maxVal=0, minVal=0, useMax=True, useMin=True)
Definition: LArMonTransforms.py:16
LArMonTransforms.digitSummary
def digitSummary(inputs, TreshOut=5, TreshSat=5, TreshNull=5)
Definition: LArMonTransforms.py:311
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
GetEntries
TGraphErrors * GetEntries(TH2F *histo)
Definition: TRTCalib_makeplots.cxx:4019
python.TrigEgammaMonitorHelper.TH2F
def TH2F(name, title, nxbins, bins_par2, bins_par3, bins_par4, bins_par5=None, bins_par6=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:45
LArMonTransforms.fillWithMaxCoverage
def fillWithMaxCoverage(inputs, isFtSlotPlot=True)
Definition: LArMonTransforms.py:50
python.TrigEgammaMonitorHelper.TProfile
def TProfile(*args, **kwargs)
Definition: TrigEgammaMonitorHelper.py:81
plotBeamSpotVxVal.range
range
Definition: plotBeamSpotVxVal.py:195
LArMonTransforms.normToEntries
def normToEntries(inputs, titleToReplace="", replaceWith="")
Definition: LArMonTransforms.py:175
histSizes.list
def list(name, path='/')
Definition: histSizes.py:38
LArMonTransforms.normToEntriesAndSetMin
def normToEntriesAndSetMin(inputs, minVal=0, maxVal=0, useMax=False, clone=True, titleToReplace="", replaceWith="")
Definition: LArMonTransforms.py:206
LArMonTransforms.computeCorrelations
def computeCorrelations(inputs)
Definition: LArMonTransforms.py:401
LArMonTransforms.normToBinAndSetMinMax
def normToBinAndSetMinMax(inputs, bin_norm=0, minVal=0, maxVal=0, useMax=False, titleToReplace="", replaceWith="", newYaxis="")
Definition: LArMonTransforms.py:243
LArMonTransforms.divideHist
def divideHist(inputs, titleToReplace="", replaceWith="")
Definition: LArMonTransforms.py:282
if
if(febId1==febId2)
Definition: LArRodBlockPhysicsV0.cxx:567
python.TrigEgammaMonitorHelper.TH1F
def TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:24
readCCLHist.float
float
Definition: readCCLHist.py:83
LArMonTransforms.Mean
def Mean(inputs)
Definition: LArMonTransforms.py:438