ATLAS Offline Software
Loading...
Searching...
No Matches
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
7import ROOT
8import os
9
10import PlotCalibrationGains
11
12import PlotCalibrationHV
13
14from optparse import OptionParser
15
16import numpy as np
17import matplotlib.pyplot as plt
18
19plt.xkcd()
20
21if __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
31 geometry_convertor = PlotCalibrationGains.L1CaloGeometryConvertor()
32 receiver_gains = PlotCalibrationGains.GainReader()
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
GetEntries
TGraphErrors * GetEntries(TH2F *histo)
Definition TRTCalib_makeplots.cxx:4025
Generated on for ATLAS Offline Software by doxygen 1.14.0