7 from prompt
import ValidationSettings
16 settings = ValidationSettings(name=
'ECL Crate',
18 download_files=[
'stdout'],
22 def run_validation(job_path, input_data_path, expert_config, **kwargs):
27 from ROOT
import TH1F, TCanvas, TFile, TGraph, TLine, kRed, gStyle
30 from array
import array
32 if not os.path.exists(
'plots'):
36 gStyle.SetOptStat(1111111)
37 gStyle.SetOptFit(1111)
41 ggRatio = TH1F(
"ggRatio",
"Ratio of output to input calib, gamma gammaRatio", 200, 0.9, 1.1)
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)
46 mumuggRatio = TH1F(
"mumuggRatio",
"Ratio of muon pair/gg calibrations, thetaID [14,57]Ratio", 200, 0.93, 1.03)
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)
55 print(
"\n---------------------------------------- \nGamma gamma comparison: \n\n")
56 gg = TFile(f
"{job_path}/ecl_gamma_gamma/0/algorithm_output/eclGammaGammaE_algorithm.root")
59 EnVsCrysID = gg.Get(
"EnVsCrysID")
60 ggEntries = EnVsCrysID.GetEntries()
61 crossSection = 3990000.
62 lumi = ggEntries/crossSection
65 estUnc = 0.5*np.sqrt(nomLum/lumi)
67 print(
"%f entries in gg EnVsCrysID lumi = %.2f fb-1 est uncerta = %.2f \n" % (ggEntries, lumi, estUnc))
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)))
80 StatusVsCrysIDgg = gg.Get(
"StatusVsCrysID")
81 AverageInitCalibgg = gg.Get(
"AverageInitCalib")
82 CalibVsCrysIDgg = gg.Get(
"CalibVsCrysID")
84 print(
"\nCompare gamma gamma output calibration to input for large changes:\n")
90 for cellID
in range(1, 8736+1):
91 status = StatusVsCrysIDgg.GetBinContent(cellID)
93 inputCalib = AverageInitCalibgg.GetBinContent(cellID)
94 outputCalib = CalibVsCrysIDgg.GetBinContent(cellID)
95 ratio = outputCalib/inputCalib
97 if(ratio < 0.95
or ratio > 1.05):
99 print(
"%.2f cellID %.4f %.3f %s\n" % (bigChange, cellID, ratio, status))
102 if(cellID >= 1297
and cellID <= 7632):
103 if(ratio < minBarrelRatio):
104 minBarrelRatio = ratio
106 if(ratio > maxBarrelRatio):
107 maxBarrelRatio = ratio
111 myC.Print(
"plots/ggRatio.pdf")
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))
117 histMean = ggRatio.GetMean()
118 ggfit = ggRatio.GetFunction(
"gaus")
119 fitMean = ggfit.GetParameter(1)
120 fitSigma = ggfit.GetParameter(2)
122 "\n\nRatio of New/old Gamma gamma calib consts: mean = %.4f mean of Gaussian fit = %.4f, sigma = %.4f\n" %
123 (histMean, fitMean, fitSigma))
127 print(
"\n---------------------------------------- \nBhabha comparison: \n\n")
128 ee = TFile(f
"{job_path}/ecl_ee5x5/0/algorithm_output/eclee5x5Algorithm.root")
131 AverageInitCalibee = ee.Get(
"AverageInitCalib")
132 CalibVsCrysIDee = ee.Get(
"CalibVsCrysID")
134 print(
"\nCompare Bhabha output calibration to input:\n")
137 for cellID
in range(1, 8736+1):
138 inputCalib = AverageInitCalibee.GetBinContent(cellID)
139 outputCalib = CalibVsCrysIDee.GetBinContent(cellID)
142 ratio = outputCalib/inputCalib
144 if(ratio < 0.95
or ratio > 1.05):
146 print(
" %2d cellID %.4f %5.3f\n" % (bigChange, cellID, ratio))
148 print(
"Total calibrated crystals = %.4f = %.1f \n\n" % (eeCalibrated, 100.*eeCalibrated/8736.))
150 myC.Print(
"plots/eeRatio.pdf")
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))
160 print(
"\n---------------------------------------- \nBhabha to gamma comparison: \n\n")
161 xcellID = array(
'd', [0]*8736)
162 yratio = array(
'd', [0]*8736)
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
176 eeggRatio.Fit(
"gaus")
177 myC.Print(
"plots/eeggRatio.pdf")
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))
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)
198 myC.Print(
"plots/eeggRatiovsCellID.pdf")
202 print(
"\n---------------------------------------- \nMuon pair comparison: \n\n")
203 mumu = TFile(f
"{job_path}/ecl_mu_mu/0/algorithm_output/eclMuMuEAlgorithm.root")
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)))
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")
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
236 mumuggRatio.Fit(
"gaus")
237 myC.Print(
"plots/mumuggRatio.pdf")
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))
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)
255 myC.Print(
"plots/mumuggRatiovsCellID.pdf")
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")
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)
271 checkMerge.Fill(newggRatio)
272 finalRatio.Fill(newOldRatio)
273 xcellID[nPo] = cellID
274 yratio[nPo] = newOldRatio
277 checkMerge.Fill(newOldRatio)
280 myC.Print(
"plots/checkMerge.pdf")
282 finalRatio.Fit(
"gaus")
283 myC.Print(
"plots/finalRatio.pdf")
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))
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)
301 myC.Print(
"plots/newOldRatiovsCellID.pdf")
305 fout = TFile(
"validateAirflow.root",
"recreate")
313 eeggRatiovsCellID.Write()
314 mumuggRatiovsCellID.Write()
315 newOldRatiovsCellID.Write()
320 if __name__ ==
"__main__":
321 run_validation(*sys.argv[1:])