Belle II Software development
ecl_energy.py
1
8
9"""
10ECL Energy Validation
11"""
12
13from prompt import ValidationSettings
14import sys
15
16
20
21
22settings = ValidationSettings(name='ECL_Energy',
23 description=__doc__,
24 download_files=['stdout'],
25 expert_config=None)
26
27
28def run_validation(job_path, input_data_path, expert_config):
29 # job_path will be replaced with path/to/calibration_results
30 # input_data_path will be replaced with path/to/data_path used for calibration, e.g. /group/belle2/dataprod/Data/PromptSkim/
31
32 # Verify that output from airflow is OK
33 from ROOT import TH1F, TCanvas, TFile, TGraph, TLine, kRed, gStyle
34 import numpy as np
35 import os
36 from array import array
37
38 if not os.path.exists('plots'):
39 os.makedirs('plots')
40
41 # Add underflow, overflow and fit result
42 gStyle.SetOptStat(1111111)
43 gStyle.SetOptFit(1111)
44
45 # ------------------------------------------------------------------------
46 # Validation histograms
47 ggRatio = TH1F("ggRatio", "Ratio of output to input calib, gamma gammaRatio", 200, 0.9, 1.1)
48
49 eeRatio = TH1F("eeRatio", "Ratio of output to input calib, Bhabha5x5Ratio", 400, 0.9, 1.1)
50 eeggRatio = TH1F("eeggRatio", "Ratio of ee/gg calibrations, thetaID [14,57]Ratio", 200, 0.95, 1.05)
51
52 mumuggRatio = TH1F("mumuggRatio", "Ratio of muon pair/gg calibrations, thetaID [14,57]Ratio", 200, 0.93, 1.03)
53
54 checkMerge = TH1F("checkMerge", "Ratio between final and gg calib valuesRatio", 200, 0., 2.)
55 finalRatio = TH1F("finalRatio", "Ratio of new/existing energy calibrationsratio", 200, 0.95, 1.05)
56
57 myC = TCanvas("myC")
58
59 # ------------------------------------------------------------------------
60 # Gamma Gamma
61 print("\n---------------------------------------- \nGamma gamma comparison: \n\n")
62 gg = TFile(f"{job_path}/ecl_gamma_gamma/0/algorithm_output/eclGammaGammaE_algorithm.root")
63
64 # Entries / luminosity
65 EnVsCrysID = gg.Get("EnVsCrysID")
66 ggEntries = EnVsCrysID.GetEntries()
67 crossSection = 3990000. # fb
68 lumi = ggEntries / crossSection
69 nomLum = 2.36 # fb-1
70 estUnc = 0.5 * np.sqrt(nomLum / lumi)
71
72 print(f"{ggEntries:f} entries in gg EnVsCrysID lumi = {lumi:.2f} fb-1 est uncerta = {estUnc:.2f} \n")
73
74 # Summarize fit status
75 hStatusgg = gg.Get("hStatus")
76 success = 100. * (hStatusgg.GetBinContent(22) + hStatusgg.GetBinContent(14)) / 8736.
77 print(f"\nSummary of Gamma Gamma fit status. {success:.1f} good fits:\n")
78 print(f"16 good fit: {hStatusgg.GetBinContent(22):.4f} \n")
79 print(f" 8 iterations: {hStatusgg.GetBinContent(14):.4f} \n")
80 print(f" 4 at limit: {hStatusgg.GetBinContent(10):.4f} \n")
81 print(f" 3 poor fit: {hStatusgg.GetBinContent(9):.4f} \n")
82 print(f"-1 low stats: {hStatusgg.GetBinContent(5):.4f} \n")
83
84 # Compare output to input (gamma gamma) calibrations
85 StatusVsCrysIDgg = gg.Get("StatusVsCrysID")
86 AverageInitCalibgg = gg.Get("AverageInitCalib")
87 CalibVsCrysIDgg = gg.Get("CalibVsCrysID")
88
89 print("\nCompare gamma gamma output calibration to input for large changes:\n")
90 bigChange = 0
91 minBarrelRatio = 99.
92 minCellID = -1
93 maxBarrelRatio = 0.
94 maxCellID = -1
95 for cellID in range(1, 8736 + 1):
96 status = StatusVsCrysIDgg.GetBinContent(cellID)
97 if status > 0:
98 inputCalib = AverageInitCalibgg.GetBinContent(cellID)
99 outputCalib = CalibVsCrysIDgg.GetBinContent(cellID)
100 ratio = outputCalib / inputCalib
101 ggRatio.Fill(ratio)
102 if(ratio < 0.95 or ratio > 1.05):
103 bigChange += 1
104 print(f"{bigChange:.2f} cellID {cellID:.4f} {ratio:.3f} {status}\n")
105
106 # Look for large deviations in the barrel excluding first and last
107 if(cellID >= 1297 and cellID <= 7632):
108 if(ratio < minBarrelRatio):
109 minBarrelRatio = ratio
110 minCellID = cellID
111 if(ratio > maxBarrelRatio):
112 maxBarrelRatio = ratio
113 maxCellID = cellID
114
115 ggRatio.Fit("gaus")
116 myC.Print("plots/ggRatio.pdf")
117
118 print("\nLargest changes in thetaID [14,57]\n")
119 print(f"ratio = {minBarrelRatio:.3f} in cellID {minCellID:f} \n")
120 print(f"ratio = {maxBarrelRatio:.3f} in cellID {maxCellID:f} \n")
121
122 histMean = ggRatio.GetMean()
123 ggfit = ggRatio.GetFunction("gaus")
124 fitMean = ggfit.GetParameter(1)
125 fitSigma = ggfit.GetParameter(2)
126 print(
127 "\n\nRatio of New/old Gamma gamma calib consts: mean = %.4f mean of Gaussian fit = %.4f, sigma = %.4f\n" %
128 (histMean, fitMean, fitSigma))
129
130 # ------------------------------------------------------------------------
131 # Bhabha
132 print("\n---------------------------------------- \nBhabha comparison: \n\n")
133 ee = TFile(f"{job_path}/ecl_ee5x5/0/algorithm_output/eclee5x5Algorithm.root")
134
135 # Compare output to input calibrations
136 AverageInitCalibee = ee.Get("AverageInitCalib")
137 CalibVsCrysIDee = ee.Get("CalibVsCrysID")
138
139 print("\nCompare Bhabha output calibration to input:\n")
140 bigChange = 0
141 eeCalibrated = 0
142 for cellID in range(1, 8736 + 1):
143 inputCalib = AverageInitCalibee.GetBinContent(cellID)
144 outputCalib = CalibVsCrysIDee.GetBinContent(cellID)
145 if(outputCalib > 0):
146 eeCalibrated += 1
147 ratio = outputCalib / inputCalib
148 eeRatio.Fill(ratio)
149 if(ratio < 0.95 or ratio > 1.05):
150 bigChange += 1
151 print(f" {int(bigChange):2} cellID {cellID:.4f} {ratio:5.3f}\n")
152
153 print(f"Total calibrated crystals = {eeCalibrated:.4f} = {100.0 * eeCalibrated / 8736.0:.1f} \n\n")
154 eeRatio.Fit("gaus")
155 myC.Print("plots/eeRatio.pdf")
156
157 histMean = eeRatio.GetMean()
158 eefit = eeRatio.GetFunction("gaus")
159 fitMean = eefit.GetParameter(1)
160 fitSigma = eefit.GetParameter(2)
161 print(f"\nBhabha mean ratio = {histMean:.4f} mean of Gaussian fit = {fitMean:.4f}, sigma = {fitSigma:.4f}\n")
162
163 # ------------------------------------------------------------------------
164 # Bhabha/gg comparison
165 print("\n---------------------------------------- \nBhabha to gamma comparison: \n\n")
166 xcellID = array('d', [0] * 8736)
167 yratio = array('d', [0] * 8736)
168 nPo = 0
169 for cellID in range(1, 8736 + 1):
170 ggstatus = StatusVsCrysIDgg.GetBinContent(cellID)
171 ggCalib = CalibVsCrysIDgg.GetBinContent(cellID)
172 eeCalib = CalibVsCrysIDee.GetBinContent(cellID)
173 if(ggstatus > 0. and eeCalib > 0.):
174 if(cellID >= 1297 and cellID <= 7632):
175 eeggRatio.Fill(eeCalib / ggCalib)
176 xcellID[nPo] = cellID
177 yratio[nPo] = eeCalib / ggCalib
178 nPo += 1
179
180 # Fit ratio
181 eeggRatio.Fit("gaus")
182 myC.Print("plots/eeggRatio.pdf")
183
184 histMean = eeggRatio.GetMean()
185 eeggfit = eeggRatio.GetFunction("gaus")
186 fitMean = eeggfit.GetParameter(1)
187 fitSigma = eeggfit.GetParameter(2)
188 print(f"\nBhabha to gg ratio mean = {histMean:.4f} mean of Gaussian fit = {fitMean:.4f}, sigma = {fitSigma:.4f}\n")
189
190 # Graph
191 eeggRatiovsCellID = TGraph(nPo, xcellID, yratio)
192 eeggRatiovsCellID.SetName("eeggRatiovsCellID")
193 eeggRatiovsCellID.SetMarkerStyle(20)
194 eeggRatiovsCellID.SetMarkerSize(0.3)
195 eeggRatiovsCellID.SetTitle("Bhabha/gg calibration ratiocellIDratio")
196 eeggRatiovsCellID.Draw("AP")
197 eeggRatiovsCellID.GetXaxis().SetRangeUser(1, 8737)
198 eeggRatiovsCellID.GetYaxis().SetRangeUser(0.95, 1.05)
199 one = TLine(1, 1, 8737, 1)
200 one.SetLineColor(kRed)
201 one.SetLineStyle(9)
202 one.Draw()
203 myC.Print("plots/eeggRatiovsCellID.pdf")
204
205 # ------------------------------------------------------------------------
206 # muon pair
207 print("\n---------------------------------------- \nMuon pair comparison: \n\n")
208 mumu = TFile(f"{job_path}/ecl_mu_mu/0/algorithm_output/eclMuMuEAlgorithm.root")
209
210 # Summarize fit status
211 hStatusmumu = mumu.Get("hStatus")
212 success = 100. * (hStatusmumu.GetBinContent(22) + hStatusmumu.GetBinContent(14)) / 8736.
213 print(r"\nSummary of muon pair fit status. \%.1f good fits:\n", success)
214 print(f"16 good fit: {hStatusmumu.GetBinContent(22):.4f} \n")
215 print(f" 8 iterations: {hStatusmumu.GetBinContent(14):.4f} \n")
216 print(f" 4 at limit: {hStatusmumu.GetBinContent(10):.4f} \n")
217 print(f" 3 poor fit: {hStatusmumu.GetBinContent(9):.4f} \n")
218 print(f" 2 no peak: {hStatusmumu.GetBinContent(8):.4f} \n")
219 print(f"-1 low stats: {hStatusmumu.GetBinContent(5):.4f} \n")
220
221 # ------------------------------------------------------------------------
222 # Muon pair/gg comparison
223 print("\n---------------------------------------- \nMuon pair to gamma comparison: \n\n")
224 StatusVsCrysIDmumu = mumu.Get("StatusVsCrysID")
225 # AverageInitCalibmumu = mumu.Get("AverageInitCalib")
226 CalibVsCrysIDmumu = mumu.Get("CalibVsCrysID")
227
228 nPo = 0
229 for cellID in range(1, 8736 + 1):
230 ggstatus = StatusVsCrysIDgg.GetBinContent(cellID)
231 ggCalib = CalibVsCrysIDgg.GetBinContent(cellID)
232 mumustatus = StatusVsCrysIDmumu.GetBinContent(cellID)
233 mumuCalib = CalibVsCrysIDmumu.GetBinContent(cellID)
234 if(ggstatus > 7.5 and mumustatus > 7.5):
235 if(cellID >= 1297 and cellID <= 7632):
236 mumuggRatio.Fill(mumuCalib / ggCalib)
237 xcellID[nPo] = cellID
238 yratio[nPo] = mumuCalib / ggCalib
239 nPo += 1
240
241 mumuggRatio.Fit("gaus")
242 myC.Print("plots/mumuggRatio.pdf")
243
244 histMean = mumuggRatio.GetMean()
245 mumuggfit = mumuggRatio.GetFunction("gaus")
246 fitMean = mumuggfit.GetParameter(1)
247 fitSigma = mumuggfit.GetParameter(2)
248 print(
249 f"\nMuon pair to gg ratio mean = {histMean:.4f} mean of Gaussian fit = {fitMean:.4f}, sigma = {fitSigma:.4f}\n")
250
251 # Graph of ratio versus cellID
252 mumuggRatiovsCellID = TGraph(nPo, xcellID, yratio)
253 mumuggRatiovsCellID.SetName("mumuggRatiovsCellID")
254 mumuggRatiovsCellID.SetMarkerStyle(20)
255 mumuggRatiovsCellID.SetMarkerSize(0.3)
256 mumuggRatiovsCellID.SetTitle("Muon pair/gg calibration ratiocellIDratio")
257 mumuggRatiovsCellID.Draw("AP")
258 mumuggRatiovsCellID.GetXaxis().SetRangeUser(1, 8737)
259 mumuggRatiovsCellID.GetYaxis().SetRangeUser(0.93, 1.03)
260 one.Draw()
261 myC.Print("plots/mumuggRatiovsCellID.pdf")
262
263 # ------------------------------------------------------------------------
264 # Check merge stage
265 print("\n---------------------------------------- \nFinal merge comparison: \n\n")
266 merge = TFile(f"{job_path}/ecl_merge/0/algorithm_output/ECLCrystalEnergy.root")
267 newPayload = merge.Get("newPayload")
268 existingPayload = merge.Get("existingPayload")
269
270 # Final calib should be equal to gamma gamma if available unchanged otherwise.
271 nPo = 0
272 for cellID in range(1, 8736 + 1):
273 ggstatus = StatusVsCrysIDgg.GetBinContent(cellID)
274 newOldRatio = newPayload.GetBinContent(cellID) / existingPayload.GetBinContent(cellID)
275 newggRatio = newPayload.GetBinContent(cellID) / CalibVsCrysIDgg.GetBinContent(cellID)
276 if(ggstatus > 0.):
277 checkMerge.Fill(newggRatio)
278 finalRatio.Fill(newOldRatio)
279 xcellID[nPo] = cellID
280 yratio[nPo] = newOldRatio
281 nPo += 1
282 else:
283 checkMerge.Fill(newOldRatio)
284
285 checkMerge.Draw()
286 myC.Print("plots/checkMerge.pdf")
287
288 finalRatio.Fit("gaus")
289 myC.Print("plots/finalRatio.pdf")
290
291 histMean = finalRatio.GetMean()
292 finalRatiofit = finalRatio.GetFunction("gaus")
293 fitMean = finalRatiofit.GetParameter(1)
294 fitSigma = finalRatiofit.GetParameter(2)
295 print(f"\nNew to exist ratio mean = {histMean:.4f} mean of Gaussian fit = {fitMean:.4f}, sigma = {fitSigma:.4f}\n")
296
297 # Graph of / old calibration constants
298 newOldRatiovsCellID = TGraph(nPo, xcellID, yratio)
299 newOldRatiovsCellID.SetName("newOldRatiovsCellID")
300 newOldRatiovsCellID.SetMarkerStyle(20)
301 newOldRatiovsCellID.SetMarkerSize(0.3)
302 newOldRatiovsCellID.SetTitle("New/existing calibration ratiocellIDratio")
303 newOldRatiovsCellID.Draw("AP")
304 newOldRatiovsCellID.GetXaxis().SetRangeUser(1, 8737)
305 newOldRatiovsCellID.GetYaxis().SetRangeUser(0.95, 1.05)
306 one.Draw()
307 myC.Print("plots/newOldRatiovsCellID.pdf")
308
309 # ------------------------------------------------------------------------
310 # Write out histograms
311 fout = TFile("validateAirflow.root", "recreate")
312 ggRatio.Write()
313 eeRatio.Write()
314 eeggRatio.Write()
315 mumuggRatio.Write()
316 checkMerge.Write()
317 finalRatio.Write()
318
319 eeggRatiovsCellID.Write()
320 mumuggRatiovsCellID.Write()
321 newOldRatiovsCellID.Write()
322 one.Draw()
323
324 fout.Close()
325
326
327if __name__ == "__main__":
328 run_validation(*sys.argv[1:])
run_validation
Definition: run_validation.py:1