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