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