ATLAS Offline Software
MDTPostProcessing.py
Go to the documentation of this file.
1 #
2 # Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
3 #
4 
5 from ROOT import TH1F
6 from .MdtMonUtils import putBoxMdtGlobal, getTubeLength, MDT2DHWName, MDTTubeEff, MDTFitTDC, putBoxMdtRegional
7 from math import sqrt
8 import numpy as np
9 
10 
11 
12 def make_hits_per_evt(inputs):
13  inputs=list(inputs)
14  i01 = inputs[0][1]
15  EvtOccBCap = i01[1].Clone()
16  EvtOccECap = i01[2].Clone()
17  EvtOccBCap.Reset()
18  EvtOccECap.Reset()
19  EvtOccBCap.SetName("HitsPerEvtInBarrelPerChamber_ADCCut")
20  EvtOccBCap.SetTitle("Avg # hits/evt Barrel, ADCCut")
21  EvtOccECap.SetName("HitsPerEvtInEndCapPerChamber_ADCCut")
22  EvtOccECap.SetTitle("Avg # hits/evt Endcap, ADCCut")
23 
24  VolumeMapBCap = i01[1].Clone()
25  VolumeMapECap = i01[2].Clone()
26  VolumeMapBCap.Reset()
27  VolumeMapECap.Reset()
28  VolumeMapBCap.SetName("VolumeMapBarrel")
29  VolumeMapBCap.SetTitle("Volume Map (#tubes*tubeVol in m^3) Barrel")
30  VolumeMapECap.SetName("VolumeMapEncap")
31  VolumeMapECap.SetTitle("Volume Map (#tubes*tubeVol in m^3) Endcap")
32 
33 
34  OccBCap_norm = i01[4].Clone()
35  OccECap_norm = i01[5].Clone()
36  OccBCap_norm.SetName("HitsPerEventInBarrelPerChamber_onSegm_ADCCut")
37  OccBCap_norm.SetTitle("Avg # hits-on-segm/evt/m^3 Barrel")
38  OccECap_norm.SetName("HitsPerEventInEndCapPerChamber_onSegm_ADCCut")
39  OccECap_norm.SetTitle("Avg # hits-on-segm/evt/m^3 Endcap")
40 
41  hflag=i01[6]
42  geo=hflag.GetMean()
43 
44  size = len(inputs)
45  chamberHits_vec = [_[1][0].GetEntries() for _ in inputs]
46 
47  sorted_chamberHits_vec = np.sort(chamberHits_vec)
48  medianChamberHits = 0
49  den = int(size/2)
50  if (size % 2 == 0 and size > 2):
51  medianChamberHits = (sorted_chamberHits_vec[den] + sorted_chamberHits_vec[den + 1])/size
52  elif(size > 0):
53  medianChamberHits = sorted_chamberHits_vec[den]
54 
55  h_trigger = inputs[0][1][3]
56  nTriggers = int(h_trigger.GetEntries())
57  for i in range(size):
58  hvOff = False
59  hi = inputs[i][1][0]
60  name = hi.GetName()
61  xAxis = name[0]+name[4]+name[3]
62  yAxis = name[1]+name[5]+name[6]
63  if((name[0:3] == "BIR" or name[0:3] == "BIM") and (name[5:7] == "11" or name[5:7] == "15")):
64  yAxis += name[2]
65  # BML[45][AC]13-->BML[56][AC]13
66  if(name[0:3] == "BML" and int(name[3]) > 3 and name[5:7] == "13"):
67  xAxis = name[0]+name[4]+str(int(name[3])+1)
68 
69  if(name[0:3] == "BME"):
70  xAxis = name[0] + name[4] + "4"
71  yAxis = "M13"
72  if(name[0:3] == "BMF"):
73  xAxis = name[0] + name[4] + (str(int(name[3])*2-1))
74 
75  nhits = hi.GetEntries()
76 
77  if (name.startswith("BE") or name.startswith("BIS8")):
78  if (nhits < 0.01*medianChamberHits+0.1):
79  hvOff = True
80  elif(nhits < 0.07 * medianChamberHits + 0.1):
81  hvOff = True
82 
83  tubeRadiusScale = 1
84  tubeLength = getTubeLength(name)
85  if(name[0:3] == "BME" or name[0:3] == "BMG"):
86  tubeRadiusScale = 0.25
87  numTubesInChamber=hi.GetNbinsX()
88  #0.0006881 m2 = pi*tube_r^2
89  chamb_vol = numTubesInChamber*tubeLength*0.0006881*tubeRadiusScale
90 
91 
92  if(len(yAxis) == 4):
93  yAxis2 = yAxis[0]+' '+yAxis[1:3]+' '+yAxis[3]
94  #yAxis2=yAxis[0]+','+yAxis[1:3]+','+yAxis[3]
95  else:
96  yAxis2 = yAxis[0]+' '+yAxis[1:3]
97  #yAxis2=yAxis[0]+','+yAxis[1:3]
98  if(xAxis.startswith("B")):
99  VolumeMapBCap.Fill(xAxis, yAxis2, chamb_vol)
100  if (not hvOff):
101  EvtOccBCap.Fill(xAxis, yAxis2, nhits/(float(nTriggers)))
102  else:
103  VolumeMapECap.Fill(xAxis, yAxis2, chamb_vol)
104  if (not hvOff):
105  EvtOccECap.Fill(xAxis, yAxis2, nhits/(float(nTriggers)))
106 
107 
108  EvtTubeNormOccBCap = EvtOccBCap.Clone()
109  EvtTubeNormOccECap = EvtOccECap.Clone()
110  EvtTubeNormOccBCap.SetName("NormedHitsInBarrelPerChamber_ADCCut")
111  EvtTubeNormOccBCap.SetTitle("Avg # hits/evt/m^3 Barrel, ADCCut")
112  EvtTubeNormOccECap.SetName("NormedHitsInEndCapPerChamber_ADCCut")
113  EvtTubeNormOccECap.SetTitle("Avg # hits/evt/m^3 Endcap, ADCCut")
114 
115  EvtTubeNormOccBCap.Divide(VolumeMapBCap)
116  EvtTubeNormOccECap.Divide(VolumeMapECap)
117 
118  OccBCap_norm.Scale(1./nTriggers)
119  OccBCap_norm.Divide(VolumeMapBCap)
120  OccECap_norm.Scale(1./nTriggers)
121  OccECap_norm.Divide(VolumeMapECap)
122 
123  putBoxMdtGlobal(EvtOccBCap, "B", geo)
124  putBoxMdtGlobal(EvtOccECap, "E", geo)
125  putBoxMdtGlobal(VolumeMapBCap, "B", geo)
126  putBoxMdtGlobal(VolumeMapECap, "E", geo)
127  putBoxMdtGlobal(EvtTubeNormOccBCap, "B", geo)
128  putBoxMdtGlobal(EvtTubeNormOccECap, "E", geo)
129  putBoxMdtGlobal(OccBCap_norm, "B", geo)
130  putBoxMdtGlobal(OccECap_norm, "E", geo)
131 
132 
133  return [EvtOccBCap, EvtOccECap, VolumeMapBCap, VolumeMapECap, EvtTubeNormOccBCap, EvtTubeNormOccECap, OccBCap_norm, OccECap_norm]
134 
135 
136 def make_eff_histo(inputs, ec):
137  ecap = ["BA", "BC", "EA", "EC"]
138  ecap_str= ecap[ec]
139  ecap_fullStr_lower = "mdt"+ecap_str
140  heff = TH1F(ecap_fullStr_lower+"_TUBE_eff",ecap_fullStr_lower+"_TUBE_eff",100,0,1)
141 
142  dencut = 10
143 
144  for itr in inputs:
145  if itr is None:
146  continue
147  hi_num = itr[1][0]
148  hi_den = itr[1][1]
149  nbin=hi_den.GetNbinsX()
150  for ibin in range(nbin):
151  if( hi_den.At(ibin) > dencut ):
152  heff.Fill( (hi_num.At(ibin))/(hi_den.At(ibin)) )
153 
154  return [heff]
155 
157  if inputs[0] is None:
158  return []
159 
160  EffBCap = inputs[0][1][2].Clone()
161  EffECap = inputs[0][1][3].Clone()
162  EffBCap.Reset()
163  EffECap.Reset()
164  EffBCap_N = EffBCap.Clone()
165  EffECap_N = EffECap.Clone()
166 
167  EffBCap.SetName("effsInBarrelPerChamber_ADCCut")
168  EffBCap.SetTitle("effsInBarrelPerChamber_ADCCut")
169  EffECap.SetName("effsInEndCapPerChamber_ADCCut")
170  EffECap.SetTitle("effsInEndCapPerChamber_ADCCut")
171 
172  hflag=inputs[0][1][4]
173  geo=hflag.GetMean()
174 
175  for itr in inputs:
176  if itr is None:
177  continue
178  hi_num = itr[1][0]
179  name_num = hi_num.GetName()
180  hi_den = itr[1][1]
181  name=name_num[0:7]
182  countsML1, countsML2, entriesML1, entriesML2 = MDTTubeEff(name,hi_num,hi_den)
183  ch_name = name[0:7]
184  stateta_IMO_c, statphi_IMO_c, stateta_c, statphi_c, statphi_c2 = MDT2DHWName(ch_name)
185  cut=10
186 
187  if( entriesML1 + entriesML2 > cut ):
188  if(name[0:1]=="B"):
189  EffBCap.Fill(stateta_IMO_c,statphi_IMO_c,countsML1+countsML2)
190  EffBCap_N.Fill(stateta_IMO_c,statphi_IMO_c,entriesML1+entriesML2)
191  else:
192  EffECap.Fill(stateta_IMO_c,statphi_IMO_c,countsML1+countsML2)
193  EffECap_N.Fill(stateta_IMO_c,statphi_IMO_c,entriesML1+entriesML2)
194 
195  EffECap.Divide(EffECap_N)
196  EffBCap.Divide(EffBCap_N)
197 
198  putBoxMdtGlobal(EffBCap, "B", geo)
199  putBoxMdtGlobal(EffECap, "E", geo)
200 
201  return [EffBCap, EffECap]
202 
203 
204 def make_eff_histo_perML(inputs, ec):
205  if inputs[0] is None:
206  return []
207 
208  ecap = ["BA", "BC", "EA", "EC"]
209  ecap_str= ecap[ec]
210 
211 
212  heff_outer = inputs[0][1][2].Clone()
213  heff_middle = inputs[0][1][3].Clone()
214  heff_inner = inputs[0][1][4].Clone()
215  heff_extra = inputs[0][1][5].Clone()
216  hflag=inputs[0][1][6]
217  geo=hflag.GetMean()
218 
219  heff_outer.Reset()
220  heff_outer.SetName("effsIn"+ecap_str+"BAOuterPerMultiLayer_ADCCut")
221  heff_outer.SetTitle("effsIn"+ecap_str+"OuterPerMultiLayer, ADCCut")
222  heff_outer_N = heff_outer.Clone()
223  name = heff_outer.GetName()
224 
225  heff_middle.Reset()
226  heff_middle.SetName("effsIn"+ecap_str+"MiddlePerMultiLayer_ADCCut")
227  heff_middle.SetTitle("effsIn"+ecap_str+"MiddlePerMultiLayer, ADCCut")
228  heff_middle_N = heff_middle.Clone()
229 
230  heff_inner.Reset()
231  heff_inner.SetName("effsIn"+ecap_str+"InnerPerMultiLayer_ADCCut")
232  heff_inner.SetTitle("effsIn"+ecap_str+"InnerPerMultiLayer, ADCCut")
233  heff_inner_N = heff_inner.Clone()
234 
235  heff_extra.Reset()
236  heff_extra.SetName("effsIn"+ecap_str+"ExtraPerMultiLayer_ADCCut")
237  heff_extra.SetTitle("effsIn"+ecap_str+"ExtraPerMultiLayer, ADCCut")
238  heff_extra_N = heff_extra.Clone()
239 
240  ecap_fullStr_lower = "mdt"+ecap_str
241  heffML = TH1F(ecap_fullStr_lower+"_ML_eff",ecap_fullStr_lower+"_ML_eff",50,0,1)
242 
243  h2=inputs[0][1][2].Clone()
244  h2.SetTitle((inputs[0][1][2].GetTitle()).replace("_forpp", ""))
245  h3=inputs[0][1][3].Clone()
246  h3.SetTitle((inputs[0][1][3].GetTitle()).replace("_forpp", ""))
247  if(not ecap_str[0] == "E"):
248  putBoxMdtRegional(h2, ecap_str[0]+"O"+ecap_str[1], geo)
249  putBoxMdtRegional(h3, ecap_str[0]+"M"+ecap_str[1], geo)
250 
251  h4=inputs[0][1][4].Clone()
252  h4.SetTitle((inputs[0][1][4].GetTitle()).replace("_forpp", ""))
253  h5=inputs[0][1][5].Clone()
254  h5.SetTitle((inputs[0][1][5].GetTitle()).replace("_forpp", ""))
255  putBoxMdtRegional(h4, ecap_str[0]+"I"+ecap_str[1], geo)
256  if(not ecap_str[0] == "B"):
257  putBoxMdtRegional(h5, ecap_str[0]+"E"+ecap_str[1], geo)
258 
259  for itr in inputs:
260  if itr is None:
261  continue
262  hi_num = itr[1][0]
263  name_num = hi_num.GetName()
264  hi_den = itr[1][1]
265  name=name_num[0:7]
266  countsML1, countsML2, entriesML1, entriesML2 = MDTTubeEff(name,hi_num,hi_den)
267  ch_name = name[0:7]
268  stateta_IMO_c, statphi_IMO_c, stateta_c, statphi_c, statphi_c2 = MDT2DHWName(ch_name)
269  cut=10
270 
271  if( entriesML2 > cut ):
272  heffML.Fill(countsML2/entriesML2)
273  if( entriesML1 > cut ):
274  heffML.Fill(countsML1/entriesML1)
275 
276  if(ch_name[1:2]=="O"):
277  if( entriesML1 > cut ):
278  heff_outer.Fill(stateta_c, statphi_c, countsML1)
279  heff_outer_N.Fill(stateta_c, statphi_c, entriesML1)
280  if( entriesML2 > cut ):
281  heff_outer.Fill(stateta_c, statphi_c2, countsML2)
282  heff_outer_N.Fill(stateta_c, statphi_c2, entriesML2)
283  if(ch_name[1:2]=="M"):
284  if( entriesML1 > cut ):
285  heff_middle.Fill(stateta_c, statphi_c, countsML1)
286  heff_middle_N.Fill(stateta_c, statphi_c, entriesML1)
287  if( entriesML2 > cut ):
288  heff_middle.Fill(stateta_c, statphi_c2, countsML2)
289  heff_middle_N.Fill(stateta_c, statphi_c2, entriesML2)
290  if(ch_name[1:2]=="I"):
291  if( entriesML1 > cut ):
292  heff_inner.Fill(stateta_c, statphi_c, countsML1)
293  heff_inner_N.Fill(stateta_c, statphi_c, entriesML1)
294  if( entriesML2 > cut ):
295  heff_inner.Fill(stateta_c, statphi_c2, countsML2)
296  heff_inner_N.Fill(stateta_c, statphi_c2, entriesML2)
297  if(ch_name[1:2]=="E"):
298  if( entriesML1 > cut ):
299  heff_extra.Fill(stateta_c, statphi_c, countsML1)
300  heff_extra_N.Fill(stateta_c, statphi_c, entriesML1)
301  if( entriesML2 > cut ):
302  heff_extra.Fill(stateta_c, statphi_c2, countsML2)
303  heff_extra_N.Fill(stateta_c, statphi_c2, entriesML2)
304 
305 
306  heff_outer.Divide(heff_outer_N)
307  heff_middle.Divide(heff_middle_N)
308  heff_inner.Divide(heff_inner_N)
309  heff_extra.Divide(heff_extra_N)
310 
311 
312  if(not ecap_str[0] == "E"):
313  putBoxMdtRegional(heff_middle,ecap_str[0]+"M"+ecap_str[1], geo)
314  putBoxMdtRegional(heff_outer,ecap_str[0]+"O"+ecap_str[1], geo)
315 
316  putBoxMdtRegional(heff_inner,ecap_str[0]+"I"+ecap_str[1], geo)
317  if(not ecap_str[0] == "B"):
318  putBoxMdtRegional(heff_extra,ecap_str[0]+"E"+ecap_str[1], geo)
319 
320 
321  return [heff_outer, heff_middle, heff_inner, heff_extra, heffML, h2, h3, h4, h5]
322 
323 def drift_time_monitoring(inputs, ec):
324 
325  ecap = ["BA", "BC", "EA", "EC"]
326  ecap_str= ecap[ec]
327 
328  size = len(inputs)
329 
330  sumt0_name = "MDT_t0_"+ecap_str
331  sumt0 = TH1F(sumt0_name,sumt0_name,1,0,1)
332  sumt0.SetFillColor(42)
333  sumt0.SetMarkerStyle(20)
334  sumt0.SetMarkerColor(42)
335  sumt0.SetAxisRange(0,300,"y")
336  sumt0.Reset()
337 
338  sumtmax_name = "MDT_tmax_"+ecap_str
339  sumtmax = TH1F(sumtmax_name,sumtmax_name,1,0,1)
340  sumtmax.SetFillColor(42)
341  sumtmax.SetMarkerStyle(20)
342  sumtmax.SetMarkerColor(42)
343  sumtmax.SetAxisRange(0,1500,"y")
344  sumtmax.Reset()
345 
346  sumtdrift_name = "MDT_tdrift_"+ecap_str
347  sumtdrift = TH1F(sumtdrift_name,sumtdrift_name,1,0,1)
348  sumtdrift.SetFillColor(42)
349  sumtdrift.SetMarkerStyle(20)
350  sumtdrift.SetMarkerColor(42)
351  sumtdrift.SetAxisRange(0,1200,"y")
352  sumtdrift.Reset()
353 
354 
355  for i in range(size):
356  currentbin=i+1
357  h = inputs[i][1][0]
358  t0, t0err, tmax, tmaxerr = MDTFitTDC(h)
359 
360  layer=""
361  if(currentbin<17):
362  layer = "In"
363  elif (currentbin<33):
364  layer = "Mid"
365  elif (currentbin<49):
366  layer = "Out"
367  else:
368  layer = "Ext"
369 
370  stPhi = currentbin%16
371  if(layer=="Ext" and (ecap_str=="BA" or ecap_str=="BC")):
372  stPhi = 2*currentbin%16
373 
374  if(stPhi==0):
375  stPhi=16
376 
377  phi = str(stPhi)
378  sector = ecap_str+"_"+layer+"_"+phi
379  if(h.GetEntries()<100):
380  sector+="_OFF"
381  elif(t0err>t0):
382  sector+="_ERR"
383 
384  sumt0.Fill(sector,t0)
385  sumt0.SetBinError(currentbin,t0err)
386  sumtmax.Fill(sector,tmax)
387  sumtmax.SetBinError(currentbin,tmaxerr)
388  if(tmax-t0 < 0 or tmax-t0 > 2000) :
389  sumtdrift.Fill(sector,0)
390  sumtdrift.SetBinError(currentbin,2000)
391  else:
392  sumtdrift.Fill(sector,tmax-t0)
393  sumtdrift.SetBinError(currentbin,sqrt(tmaxerr*tmaxerr + t0err*t0err))
394 
395 
396  sumt0.LabelsDeflate("x")
397  sumtmax.LabelsDeflate("x")
398  sumtdrift.LabelsDeflate("x")
399  if(sumt0.GetNbinsX()>1) :
400  sumt0.LabelsOption("v","x")
401  if(sumtmax.GetNbinsX()>1) :
402  sumtmax.LabelsOption("v","x")
403  if(sumtdrift.GetNbinsX()>1) :
404  sumtdrift.LabelsOption("v","x")
405 
406  return [sumt0, sumtmax, sumtdrift]
407 
408 def MdtGlobalBox(inputs):
409  EvtOccBCap = inputs[0][1][0].Clone()
410  EvtOccBCap.SetTitle((inputs[0][1][0].GetTitle()).replace("_forpp", ""))
411  EvtOccECap = inputs[0][1][1].Clone()
412  EvtOccECap.SetTitle((inputs[0][1][1].GetTitle()).replace("_forpp", ""))
413  hflag=inputs[0][1][2]
414  geo=hflag.GetMean()
415  putBoxMdtGlobal(EvtOccBCap, "B", geo)
416  putBoxMdtGlobal(EvtOccECap, "E", geo)
417  return [EvtOccBCap, EvtOccECap]
418 
419 
420 def test(inputs):
421  """ HitsPerEventInXXPerChamber_[onSegm]_ADCCut """
422  print("hello!")
423  print(inputs)
424  print(len(inputs))
425  return []
replace
std::string replace(std::string s, const std::string &s2, const std::string &s3)
Definition: hcg.cxx:307
MDTPostProcessing.test
def test(inputs)
Definition: MDTPostProcessing.py:420
MDTPostProcessing.make_eff_histo_perML
def make_eff_histo_perML(inputs, ec)
Definition: MDTPostProcessing.py:204
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
MdtMonUtils.MDTTubeEff
def MDTTubeEff(name, hi_num, hi_den)
Definition: MdtMonUtils.py:509
MDTPostProcessing.make_eff_histo
def make_eff_histo(inputs, ec)
Definition: MDTPostProcessing.py:136
MDTPostProcessing.drift_time_monitoring
def drift_time_monitoring(inputs, ec)
Definition: MDTPostProcessing.py:323
plotBeamSpotVxVal.range
range
Definition: plotBeamSpotVxVal.py:195
histSizes.list
def list(name, path='/')
Definition: histSizes.py:38
MDTPostProcessing.MdtGlobalBox
def MdtGlobalBox(inputs)
Definition: MDTPostProcessing.py:408
MDTPostProcessing.make_eff_histo_perChamber
def make_eff_histo_perChamber(inputs)
Definition: MDTPostProcessing.py:156
MdtMonUtils.putBoxMdtRegional
def putBoxMdtRegional(h, xAxis, geo)
Definition: MdtMonUtils.py:19
MdtMonUtils.getTubeLength
def getTubeLength(name)
Definition: MdtMonUtils.py:456
MdtMonUtils.MDTFitTDC
def MDTFitTDC(h)
Definition: MdtMonUtils.py:546
generate::GetEntries
double GetEntries(TH1D *h, int ilow, int ihi)
Definition: rmsFrac.cxx:20
Muon::print
std::string print(const MuPatSegment &)
Definition: MuonTrackSteering.cxx:28
if
if(febId1==febId2)
Definition: LArRodBlockPhysicsV0.cxx:569
MdtMonUtils.MDT2DHWName
def MDT2DHWName(name)
Definition: MdtMonUtils.py:472
str
Definition: BTagTrackIpAccessor.cxx:11
MDTPostProcessing.make_hits_per_evt
def make_hits_per_evt(inputs)
Definition: MDTPostProcessing.py:12
python.TrigEgammaMonitorHelper.TH1F
def TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:24
MdtMonUtils.putBoxMdtGlobal
def putBoxMdtGlobal(h, ecap, geo)
Definition: MdtMonUtils.py:169
readCCLHist.float
float
Definition: readCCLHist.py:83