Belle II Software  release-08-01-10
ecl_energy.py
1 
8 
9 """
10 ECL Energy Validation
11 """
12 
13 from prompt import ValidationSettings
14 import sys
15 
16 
20 
21 
22 settings = ValidationSettings(name='ECL_Energy',
23  description=__doc__,
24  download_files=['stdout'],
25  expert_config=None)
26 
27 
28 def 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} entries in gg EnVsCrysID lumi = {:.2f} fb-1 est uncerta = {:.2f} \n".format(ggEntries, lumi, estUnc))
73 
74  # Summarize fit status
75  hStatusgg = gg.Get("hStatus")
76  success = 100. * (hStatusgg.GetBinContent(22) + hStatusgg.GetBinContent(14)) / 8736.
77  print("\nSummary of Gamma Gamma fit status. %.1f good fits:\n" % (success))
78  print("16 good fit: %.4f \n" % (hStatusgg.GetBinContent(22)))
79  print(" 8 iterations: %.4f \n" % (hStatusgg.GetBinContent(14)))
80  print(" 4 at limit: %.4f \n" % (hStatusgg.GetBinContent(10)))
81  print(" 3 poor fit: %.4f \n" % (hStatusgg.GetBinContent(9)))
82  print("-1 low stats: %.4f \n" % (hStatusgg.GetBinContent(5)))
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("{:.2f} cellID {:.4f} {:.3f} {}\n".format(bigChange, cellID, ratio, status))
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("ratio = {:.3f} in cellID {:f} \n".format(minBarrelRatio, minCellID))
120  print("ratio = {:.3f} in cellID {:f} \n".format(maxBarrelRatio, maxCellID))
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(" %2d cellID %.4f %5.3f\n" % (bigChange, cellID, ratio))
152 
153  print("Total calibrated crystals = {:.4f} = {:.1f} \n\n".format(eeCalibrated, 100. * eeCalibrated / 8736.))
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("\nBhabha mean ratio = {:.4f} mean of Gaussian fit = {:.4f}, sigma = {:.4f}\n".format(histMean, fitMean, fitSigma))
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("\nBhabha to gg ratio mean = {:.4f} mean of Gaussian fit = {:.4f}, sigma = {:.4f}\n".format(histMean, fitMean, fitSigma))
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("16 good fit: %.4f \n" % (hStatusmumu.GetBinContent(22)))
215  print(" 8 iterations: %.4f \n" % (hStatusmumu.GetBinContent(14)))
216  print(" 4 at limit: %.4f \n" % (hStatusmumu.GetBinContent(10)))
217  print(" 3 poor fit: %.4f \n" % (hStatusmumu.GetBinContent(9)))
218  print(" 2 no peak: %.4f \n" % (hStatusmumu.GetBinContent(8)))
219  print("-1 low stats: %.4f \n" % (hStatusmumu.GetBinContent(5)))
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  "\nMuon pair to gg ratio mean = {:.4f} mean of Gaussian fit = {:.4f}, sigma = {:.4f}\n".format(
250  histMean,
251  fitMean,
252  fitSigma))
253 
254  # Graph of ratio versus cellID
255  mumuggRatiovsCellID = TGraph(nPo, xcellID, yratio)
256  mumuggRatiovsCellID.SetName("mumuggRatiovsCellID")
257  mumuggRatiovsCellID.SetMarkerStyle(20)
258  mumuggRatiovsCellID.SetMarkerSize(0.3)
259  mumuggRatiovsCellID.SetTitle("Muon pair/gg calibration ratiocellIDratio")
260  mumuggRatiovsCellID.Draw("AP")
261  mumuggRatiovsCellID.GetXaxis().SetRangeUser(1, 8737)
262  mumuggRatiovsCellID.GetYaxis().SetRangeUser(0.93, 1.03)
263  one.Draw()
264  myC.Print("plots/mumuggRatiovsCellID.pdf")
265 
266  # ------------------------------------------------------------------------
267  # Check merge stage
268  print("\n---------------------------------------- \nFinal merge comparison: \n\n")
269  merge = TFile(f"{job_path}/ecl_merge/0/algorithm_output/ECLCrystalEnergy.root")
270  newPayload = merge.Get("newPayload")
271  existingPayload = merge.Get("existingPayload")
272 
273  # Final calib should be equal to gamma gamma if available unchanged otherwise.
274  nPo = 0
275  for cellID in range(1, 8736 + 1):
276  ggstatus = StatusVsCrysIDgg.GetBinContent(cellID)
277  newOldRatio = newPayload.GetBinContent(cellID) / existingPayload.GetBinContent(cellID)
278  newggRatio = newPayload.GetBinContent(cellID) / CalibVsCrysIDgg.GetBinContent(cellID)
279  if(ggstatus > 0.):
280  checkMerge.Fill(newggRatio)
281  finalRatio.Fill(newOldRatio)
282  xcellID[nPo] = cellID
283  yratio[nPo] = newOldRatio
284  nPo += 1
285  else:
286  checkMerge.Fill(newOldRatio)
287 
288  checkMerge.Draw()
289  myC.Print("plots/checkMerge.pdf")
290 
291  finalRatio.Fit("gaus")
292  myC.Print("plots/finalRatio.pdf")
293 
294  histMean = finalRatio.GetMean()
295  finalRatiofit = finalRatio.GetFunction("gaus")
296  fitMean = finalRatiofit.GetParameter(1)
297  fitSigma = finalRatiofit.GetParameter(2)
298  print("\nNew to exist ratio mean = {:.4f} mean of Gaussian fit = {:.4f}, sigma = {:.4f}\n".format(histMean, fitMean, fitSigma))
299 
300  # Graph of / old calibration constants
301  newOldRatiovsCellID = TGraph(nPo, xcellID, yratio)
302  newOldRatiovsCellID.SetName("newOldRatiovsCellID")
303  newOldRatiovsCellID.SetMarkerStyle(20)
304  newOldRatiovsCellID.SetMarkerSize(0.3)
305  newOldRatiovsCellID.SetTitle("New/existing calibration ratiocellIDratio")
306  newOldRatiovsCellID.Draw("AP")
307  newOldRatiovsCellID.GetXaxis().SetRangeUser(1, 8737)
308  newOldRatiovsCellID.GetYaxis().SetRangeUser(0.95, 1.05)
309  one.Draw()
310  myC.Print("plots/newOldRatiovsCellID.pdf")
311 
312  # ------------------------------------------------------------------------
313  # Write out histograms
314  fout = TFile("validateAirflow.root", "recreate")
315  ggRatio.Write()
316  eeRatio.Write()
317  eeggRatio.Write()
318  mumuggRatio.Write()
319  checkMerge.Write()
320  finalRatio.Write()
321 
322  eeggRatiovsCellID.Write()
323  mumuggRatiovsCellID.Write()
324  newOldRatiovsCellID.Write()
325  one.Draw()
326 
327  fout.Close()
328 
329 
330 if __name__ == "__main__":
331  run_validation(*sys.argv[1:])
run_validation
Definition: run_validation.py:1