ATLAS Offline Software
checkHVCorrections.py
Go to the documentation of this file.
1 #!/bin/env python
2 #
3 # Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
4 #
5 
6 
7 import ROOT
8 import os
9 
10 import PlotCalibrationGains
11 
12 import PlotCalibrationHV
13 
14 from optparse import OptionParser
15 
16 import numpy as np
17 import matplotlib.pyplot as plt
18 
19 plt.xkcd()
20 
21 if __name__ == "__main__":
22 
23  print ("Starting checkHVCorrections script")
24 
25  parser = OptionParser()
26 
27  parser.add_option("-f","--InputFile",action="store",type="string",dest="input_file_name",help="Name of input file")
28  (options, args) = parser.parse_args()
29 
30 
33  geometry_convertor.LoadReceiverPPMMap()
34 
35  ROOT.gStyle.SetPalette(1)
36  ROOT.gStyle.SetOptStat(111111)
37  ROOT.gStyle.SetCanvasColor(10)
38 
39  c1 = ROOT.TCanvas('c1','Example',200,10,700,500)
40  c2 = ROOT.TCanvas('c2','Example Partitions',200,10,700,500)
41  c2.Divide(3,2)
42 
43 # now define histograms
44 
45  h_gains_em = PlotCalibrationGains.L1CaloMap("Eta-phi map of EM gains","#eta bin","#phi bin")
46  h_gains_had = PlotCalibrationGains.L1CaloMap("Eta-phi map of HAD gains","#eta bin","#phi bin")
47 
48  h_gains_em_noHV = PlotCalibrationGains.L1CaloMap("Eta-phi map of EM gains without HV corrections","#eta bin","#phi bin")
49  h_gains_had_noHV = PlotCalibrationGains.L1CaloMap("Eta-phi map of HAD gains without HV corrections","#eta bin","#phi bin")
50 
51  h_gains_em_ratio = PlotCalibrationGains.L1CaloMap("ratio of EM gains with and without HV corrections","#eta bin","#phi bin")
52  h_gains_had_ratio = PlotCalibrationGains.L1CaloMap("ratio of HAD gains with and without HV corrections","#eta bin","#phi bin")
53 
54  h1D_gains_em_ratio = PlotCalibrationGains.EmPartitionPlots("ratio of EM gains with and without HV corrections",40,0.,2.5,"gain/gain(no HV)","N")
55  h1D_gains_had_ratio = PlotCalibrationGains.HadPartitionPlots("ratio of HAD gains with and without HV corrections",40,0.,2.5,"gain/gain(no HV)","N")
56 
57  h_corr_em = PlotCalibrationGains.L1CaloMap("Eta-phi map of mean HV corrections in EM layer","#eta bin","#phi bin")
58  h_corr_had = PlotCalibrationGains.L1CaloMap("Eta-phi map of mean HV corrections in HAD layer","#eta bin","#phi bin")
59 
60  h_corrGainPredictor_em = PlotCalibrationGains.L1CaloMap("Eta-phi map of HV corrections from Gain Predictor in EM layer","#eta bin","#phi bin")
61  h_corrGainPredictor_had = PlotCalibrationGains.L1CaloMap("Eta-phi map of HV corrections from Gain Predictor in HAD layer","#eta bin","#phi bin")
62 
63 
64  h1D_corr_em = PlotCalibrationGains.EmPartitionPlots("mean HV corrections in EM layer",40,0.,2.5,"mean HV correction","N")
65  h1D_corr_had = PlotCalibrationGains.HadPartitionPlots("mean HV corrections in HAD layer",40,0.,2.5,"mean HV correction","N")
66 
67  h_corrHVdiff_em = PlotCalibrationGains.L1CaloMap("(ratio-meanHV)/meanHV in EM layer","#eta bin","#phi bin")
68  h_corrHVdiff_had = PlotCalibrationGains.L1CaloMap("(ratio-meanHV)/meanHV in HAD layer","#eta bin","#phi bin")
69 
70  h1D_corrHVdiff_em = PlotCalibrationGains.EmPartitionPlots("(ratio-meanHV)/meanHV in EM layer",40,-0.1,0.1,"(ratio-meanHV)/meanHV","N")
71  h1D_corrHVdiff_had = PlotCalibrationGains.HadPartitionPlots("(ratio-meanHV)/meanHV in HAD layer",40,-0.1,0.1,"(ratio-meanHV)/meanHV","N")
72 
73  h_corrHVdiffGainPredictor_em = PlotCalibrationGains.L1CaloMap("(ratio-GainPredictor)/GainPredictor in EM layer","#eta bin","#phi bin")
74  h_corrHVdiffGainPredictor_had = PlotCalibrationGains.L1CaloMap("(ratio-GainPredictor)/GainPredictor in HAD layer","#eta bin","#phi bin")
75 
76 
77 # input files are define here
78 
79  fileDefault = '/afs/cern.ch/work/j/juraj/public/testarea/21.0.18/CalibrationProcessing/326189_HV/energyscanresults.sqlite'
80  fileNoHV = '/afs/cern.ch/work/j/juraj/public/testarea/21.0.18/CalibrationProcessing/326189/energyscanresults.sqlite'
81 
82  hvCorrFile = '/afs/cern.ch/work/j/juraj/public/testarea/20.7.8.3/HVDumps/hvcorrections_9jun17.sqlite'
83 
84  hvGainPredictor = '/afs/cern.ch/work/j/juraj/public/testarea/20.7.8.3/HVDumps/hvUpdate_9jun17_0p.txt'
85 
86  strategyString = 'GainOneOvEmecFcalLowEta' # defines receiver to TT mapping for EMB/EMEC overlap and hadronic FCAL
87 # strategyString = 'GainOneOvEmbFcalHighEta'
88 
89 
90  receiver_gains.LoadGainsSqlite(fileDefault)
91  receiver_gains.LoadReferenceSqlite(fileNoHV)
92 
93  hv_status = PlotCalibrationHV.L1CaloHVReader(hvCorrFile)
94 
95  strange_channel_file = open('checkHVCorrections.txt','w')
96 
97  gpReceiver,gpCool,gpEta,gpPhi,gpCorrection = np.genfromtxt(hvGainPredictor,unpack=True,dtype='str')
98 
99  gainPredictorCorrections = {}
100 
101  for iii in range(len(gpReceiver)):
102  gainPredictorCorrections[gpReceiver[iii]] = float(gpCorrection[iii])
103 
104  for i_eta in range(-49,45):
105  for i_phi in range(0,64):
106 
107  coolEm = geometry_convertor.getCoolEm(i_eta,i_phi)
108  coolHad = geometry_convertor.getCoolHad(i_eta,i_phi)
109 
110 
111  if not coolEm == '': # there is a channel for this eta-phi
112 
113  gain = receiver_gains.getGain(coolEm)
114  reference_gain = receiver_gains.getReferenceGain(coolEm)
115  receiverEm = geometry_convertor.getReceiverfromPPM(coolEm,strategyString)
116 
117  try:
118  meanEmHVCorrection = (hv_status.GetMeanCorections())[receiverEm]
119  except KeyError:
120  meanEmHVCorrection = 1.
121 
122  try:
123  gainPredictorEmHVCorrection = gainPredictorCorrections[receiverEm]
124  except KeyError:
125  gainPredictorEmHVCorrection = 1.
126 
127  if (not gain == '') and (not reference_gain == ''): # both gains should be available
128 
129  h_gains_em.Fill(i_eta,i_phi,gain)
130  h_gains_em_noHV.Fill(i_eta,i_phi,reference_gain)
131  h_corr_em.Fill(i_eta,i_phi,meanEmHVCorrection)
132  h_corrGainPredictor_em.Fill(i_eta,i_phi,gainPredictorEmHVCorrection)
133 
134  h1D_corr_em.Fill(i_eta,meanEmHVCorrection)
135  if reference_gain > 0.1 :
136  if gain/reference_gain > 0.01: # to make histo look nicer
137  h_gains_em_ratio.Fill(i_eta,i_phi,gain/reference_gain)
138  h1D_gains_em_ratio.Fill(i_eta,gain/reference_gain)
139 
140  if gainPredictorEmHVCorrection > 0.00001: # just to avoid division by zero, usually >=1.
141  relDifference = (gain/reference_gain - gainPredictorEmHVCorrection)/gainPredictorEmHVCorrection
142 
143  if relDifference > 0.001: # to make histo look nicer
144  h_corrHVdiffGainPredictor_em.Fill(i_eta,i_phi,relDifference)
145  if relDifference > 0.05: # print out interesting channels
146  strange_channel_file.write('%s %s %i %i EM corr(RampMaker)= %.4f corr(GainPredictor)= %.4f diff= %.4f \n' % (coolEm , receiverEm, \
147  i_eta , i_phi, gain/reference_gain,gainPredictorEmHVCorrection,relDifference) )
148 
149  if meanEmHVCorrection > 0.00001: # just to avoid division by zero, usually >=1.
150  relDifference = (gain/reference_gain - meanEmHVCorrection)/meanEmHVCorrection
151  if relDifference > 0.01: # to make histo look nicer
152  h_corrHVdiff_em.Fill(i_eta,i_phi,relDifference)
153 
154  h1D_corrHVdiff_em.Fill(i_eta,relDifference)
155 
156  if not coolHad == '': # there is a channel for this eta-phi
157 
158  gain = receiver_gains.getGain(coolHad)
159  reference_gain = receiver_gains.getReferenceGain(coolHad)
160  receiverHad = geometry_convertor.getReceiverfromPPM(coolHad,strategyString)
161 
162  try:
163  meanHadHVCorrection = (hv_status.GetMeanCorections())[receiverHad]
164  except KeyError:
165  meanHadHVCorrection = 1.
166 
167  try:
168  gainPredictorHadHVCorrection = gainPredictorCorrections[receiverHad]
169  except KeyError:
170  gainPredictorHadHVCorrection = 1.
171 
172 
173  if (not gain == '') and (not reference_gain == ''): # both gains should be available
174 
175  h_gains_had.Fill(i_eta,i_phi,gain)
176  h_gains_had_noHV.Fill(i_eta,i_phi,reference_gain)
177  h_corr_had.Fill(i_eta,i_phi,meanHadHVCorrection)
178  h_corrGainPredictor_had.Fill(i_eta,i_phi,gainPredictorHadHVCorrection)
179 
180  h1D_corr_had.Fill(i_eta,meanHadHVCorrection)
181 
182  if reference_gain > 0.1 :
183  if gain/reference_gain > 0.01: # to make histo look nicer
184  h_gains_had_ratio.Fill(i_eta,i_phi,gain/reference_gain)
185  h1D_gains_had_ratio.Fill(i_eta,gain/reference_gain)
186 
187  if gainPredictorHadHVCorrection > 0.00001: # to avoid accidental zeros
188  relDifference = (gain/reference_gain - gainPredictorHadHVCorrection)/gainPredictorHadHVCorrection
189  if relDifference > 0.001: # to make histo look nicer
190  h_corrHVdiffGainPredictor_had.Fill(i_eta,i_phi,relDifference)
191  if relDifference > 0.05: # print interesting channels
192  strange_channel_file.write('%s %s %i %i HAD corr(RampMaker)= %.4f corr(GainPredictor)= %.4f diff= %.4f \n' % (coolHad,receiverHad, \
193  i_eta , i_phi, gain/reference_gain,gainPredictorHadHVCorrection,relDifference) )
194 
195 
196  if meanHadHVCorrection > 0.00001:
197  relDifference = (gain/reference_gain - meanHadHVCorrection)/meanHadHVCorrection
198  if relDifference > 0.01: # to make histo look nicer
199  h_corrHVdiff_had.Fill(i_eta,i_phi,relDifference)
200  h1D_corrHVdiff_had.Fill(i_eta,relDifference)
201 
202 
203  c1.cd()
204  ROOT.gPad.SetLogy(0)
205 
206  h_gains_em.SetMinimum(0.6)
207  h_gains_em.SetMaximum(1.4)
208  h_gains_em.Draw()
209  c1.Print("checkHVCorrections.ps(")
210 
211  h_gains_had.SetMinimum(0.6)
212  h_gains_had.SetMaximum(1.4)
213  h_gains_had.Draw()
214  c1.Print("checkHVCorrections.ps")
215 
216  h_gains_em_noHV.SetMinimum(0.6)
217  h_gains_em_noHV.SetMaximum(1.4)
218  h_gains_em_noHV.Draw()
219  c1.Print("checkHVCorrections.ps")
220 
221  h_gains_had_noHV.SetMinimum(0.6)
222  h_gains_had_noHV.SetMaximum(1.4)
223  h_gains_had_noHV.Draw()
224  c1.Print("checkHVCorrections.ps")
225 
226  h_gains_em_ratio.SetMinimum(1.05)
227  h_gains_em_ratio.SetMaximum(2.2)
228  h_gains_em_ratio.Draw()
229  c1.Print("checkHVCorrections.ps")
230 
231  h_gains_had_ratio.SetMinimum(1.05)
232  h_gains_had_ratio.SetMaximum(2.2)
233  h_gains_had_ratio.Draw()
234  c1.Print("checkHVCorrections.ps")
235 
236  h_corr_em.SetMinimum(1.05)
237  h_corr_em.SetMaximum(2.2)
238  h_corr_em.Draw()
239  c1.Print("checkHVCorrections.ps")
240 
241  h_corr_had.SetMinimum(1.05)
242  h_corr_had.SetMaximum(2.2)
243  h_corr_had.Draw()
244  c1.Print("checkHVCorrections.ps")
245 
246  h_corrGainPredictor_em.SetMinimum(1.0)
247  h_corrGainPredictor_em.SetMaximum(2.1)
248  h_corrGainPredictor_em.Draw()
249  c1.Print("checkHVCorrections.ps")
250 
251  h_corrGainPredictor_had.SetMinimum(1.0)
252  h_corrGainPredictor_had.SetMaximum(2.1)
253  h_corrGainPredictor_had.Draw()
254  c1.Print("checkHVCorrections.ps")
255 
256 
257  h_corrHVdiff_em.SetMinimum(-.1)
258  h_corrHVdiff_em.SetMaximum( .1)
259  h_corrHVdiff_em.Draw()
260  c1.Print("checkHVCorrections.ps")
261 
262  h_corrHVdiff_had.SetMinimum(-.1)
263  h_corrHVdiff_had.SetMaximum( .1)
264  h_corrHVdiff_had.Draw()
265  c1.Print("checkHVCorrections.ps")
266 
267  h_corrHVdiffGainPredictor_em.SetMinimum(-.5)
268  h_corrHVdiffGainPredictor_em.SetMaximum( .5)
269  h_corrHVdiffGainPredictor_em.Draw()
270  c1.Print("checkHVCorrections.ps")
271 
272  h_corrHVdiffGainPredictor_had.SetMinimum(-.5)
273  h_corrHVdiffGainPredictor_had.SetMaximum( .5)
274  h_corrHVdiffGainPredictor_had.Draw()
275  c1.Print("checkHVCorrections.ps")
276 
277  h_corrHVdiffGainPredictor_em.SetMinimum(-.1)
278  h_corrHVdiffGainPredictor_em.SetMaximum( .1)
279  h_corrHVdiffGainPredictor_em.Draw()
280  c1.Print("checkHVCorrections.ps")
281 
282  h_corrHVdiffGainPredictor_had.SetMinimum(-.1)
283  h_corrHVdiffGainPredictor_had.SetMaximum( .1)
284  h_corrHVdiffGainPredictor_had.Draw()
285  c1.Print("checkHVCorrections.ps")
286 
287  h_corrHVdiffGainPredictor_em.SetMinimum(-.05)
288  h_corrHVdiffGainPredictor_em.SetMaximum( .05)
289  h_corrHVdiffGainPredictor_em.Draw()
290  c1.Print("checkHVCorrections.ps")
291 
292  h_corrHVdiffGainPredictor_had.SetMinimum(-.05)
293  h_corrHVdiffGainPredictor_had.SetMaximum( .05)
294  h_corrHVdiffGainPredictor_had.Draw()
295  c1.Print("checkHVCorrections.ps")
296 
297  c2.Clear()
298  c2.Divide(3,2)
299  for i_p in range(0,h1D_gains_em_ratio.nPartitions):
300  c2.cd(i_p+1)
301  if h1D_gains_em_ratio.his_partitions[i_p].GetEntries() > 0:
302  ROOT.gPad.SetLogy()
303  else:
304  ROOT.gPad.SetLogy(0)
305  h1D_gains_em_ratio.his_partitions[i_p].Draw()
306 
307  c2.Print("checkHVCorrections.ps")
308 
309  #c2.cd()
310  c2.Clear()
311  c2.Divide(3,2)
312  for i_p in range(0,h1D_gains_had_ratio.nPartitions):
313  c2.cd(i_p+1)
314  if h1D_gains_had_ratio.his_partitions[i_p].GetEntries() > 0:
315  ROOT.gPad.SetLogy()
316  else:
317  ROOT.gPad.SetLogy(0)
318  h1D_gains_had_ratio.his_partitions[i_p].Draw()
319 
320  c2.Print("checkHVCorrections.ps")
321 
322  c2.Clear()
323  c2.Divide(3,2)
324  for i_p in range(0,h1D_corr_em.nPartitions):
325  c2.cd(i_p+1)
326  if h1D_corr_em.his_partitions[i_p].GetEntries() > 0:
327  ROOT.gPad.SetLogy()
328  else:
329  ROOT.gPad.SetLogy(0)
330  h1D_corr_em.his_partitions[i_p].Draw()
331 
332  c2.Print("checkHVCorrections.ps")
333 
334  #c2.cd()
335  c2.Clear()
336  c2.Divide(3,2)
337  for i_p in range(0,h1D_corr_had.nPartitions):
338  c2.cd(i_p+1)
339  if h1D_corr_had.his_partitions[i_p].GetEntries() > 0:
340  ROOT.gPad.SetLogy()
341  else:
342  ROOT.gPad.SetLogy(0)
343  h1D_corr_had.his_partitions[i_p].Draw()
344 
345  c2.Print("checkHVCorrections.ps")
346  c2.Clear()
347  c2.Divide(3,2)
348  for i_p in range(0,h1D_corrHVdiff_em.nPartitions):
349  c2.cd(i_p+1)
350  if h1D_corrHVdiff_em.his_partitions[i_p].GetEntries() > 0:
351  ROOT.gPad.SetLogy()
352  else:
353  ROOT.gPad.SetLogy(0)
354  h1D_corrHVdiff_em.his_partitions[i_p].Draw()
355 
356  c2.Print("checkHVCorrections.ps")
357 
358  #c2.cd()
359  c2.Clear()
360  c2.Divide(3,2)
361  for i_p in range(0,h1D_corrHVdiff_had.nPartitions):
362  c2.cd(i_p+1)
363  if h1D_corrHVdiff_had.his_partitions[i_p].GetEntries() > 0:
364  ROOT.gPad.SetLogy()
365  else:
366  ROOT.gPad.SetLogy(0)
367  h1D_corrHVdiff_had.his_partitions[i_p].Draw()
368 
369  c2.Print("checkHVCorrections.ps)")
370 
371  os.system("ps2pdf checkHVCorrections.ps")
372 
373  strange_channel_file.close()
374 
PlotCalibrationHV.L1CaloHVReader
Definition: PlotCalibrationHV.py:120
GetEntries
TGraphErrors * GetEntries(TH2F *histo)
Definition: TRTCalib_makeplots.cxx:4019
plotBeamSpotVxVal.range
range
Definition: plotBeamSpotVxVal.py:194
PlotCalibrationGains.L1CaloGeometryConvertor
Definition: PlotCalibrationGains.py:83
PlotCalibrationGains.GainReader
Definition: PlotCalibrationGains.py:296
PlotCalibrationGains.EmPartitionPlots
Definition: PlotCalibrationGains.py:517
Trk::open
@ open
Definition: BinningType.h:40
PlotCalibrationGains.HadPartitionPlots
Definition: PlotCalibrationGains.py:562
PlotCalibrationGains.L1CaloMap
Definition: PlotCalibrationGains.py:37
python.LArMinBiasAlgConfig.float
float
Definition: LArMinBiasAlgConfig.py:65