16 from basf2
import B2INFO
17 import flavorTagger
as ft
18 from defaultEvaluationParameters
import categories, Quiet, r_subsample, r_size
19 from array
import array
25 ROOT.gROOT.SetBatch(
True)
27 if len(sys.argv) != 3:
28 sys.exit(
"Must provide 2 arguments: [input_sim_file] or ['input_sim_file*'] with wildcards and [treeName]"
30 workingFile = sys.argv[1]
31 workingFiles = glob.glob(str(workingFile))
32 treeName = str(sys.argv[2])
34 if len(workingFiles) < 1:
35 sys.exit(
"No file name or file names " + str(workingFile) +
" found.")
38 workingDirectory =
'.'
53 tree = ROOT.TChain(treeName)
55 mcstatus = array(
'd', [-511.5, 0.0, 511.5])
56 ROOT.TH1.SetDefaultSumw2()
58 for iFile
in workingFiles:
62 for branch
in tree.GetListOfBranches():
63 totalBranches.append(branch.GetName())
65 if 'FBDT_qrCombined' in totalBranches:
66 methods.append(
"FBDT")
68 if 'FANN_qrCombined' in totalBranches:
69 methods.append(
"FANN")
71 if 'DNN_qrCombined' in totalBranches:
75 for cat
in categories:
76 catBranch =
'qp' + cat
77 if catBranch
in totalBranches:
78 usedCategories.append(cat)
80 if len(usedCategories) > 1:
81 ft.WhichCategories(usedCategories)
87 for method
in methods:
90 histo_avr_r = ROOT.TH1F(
'Average_r',
'Average r in each of the bins (B0 and B0bar)', int(r_size - 2),
92 histo_avr_rB0 = ROOT.TH1F(
'Average_rB0',
'Average r in each of the bins (B0)', int(r_size - 2),
94 histo_avr_rB0bar = ROOT.TH1F(
'Average_rB0bar',
'Average r in each of the bins (B0bar)', int(r_size - 2),
98 histo_mc_NwB0 = ROOT.TH1F(
'mc_NwB0',
'Average r in each of the bins (B0)', int(r_size - 2),
100 histo_mc_NwB0bar = ROOT.TH1F(
'mc_NwB0bar',
'Average r in each of the bins (B0bar)', int(r_size - 2),
104 histo_ms_r = ROOT.TH1F(
'MS_r',
'Mean squared average of r in each of the bins (B0 and B0bar)', int(r_size - 2),
106 histo_ms_rB0 = ROOT.TH1F(
'MS_rB0',
'Mean squared average of r in each of the bins (B0)', int(r_size - 2),
108 histo_ms_rB0bar = ROOT.TH1F(
'MS_rB0bar',
'Mean squared average of r in each of the bins (B0bar)', int(r_size - 2),
112 histo_entries_per_bin = ROOT.TH1F(
114 'Events binned in r_subsample according to their r-value for B0 and B0bar prob',
117 histo_entries_per_binB0 = ROOT.TH1F(
'entries_per_binB0',
'Events binned in r_subsample according '
118 'to their r-value for B0 prob', int(r_size - 2), r_subsample)
119 histo_entries_per_binB0bar = ROOT.TH1F(
'entries_per_binB0bar',
120 'Events binned in r_subsample according to their r-value '
121 'for B0bar prob', int(r_size - 2), r_subsample)
123 histo_Cnet_output_B0 = ROOT.TH1F(
'Comb_Net_Output_B0',
'Combiner network output [not equal to r] '
124 'for true B0 (binning 100)', 100, 0.0, 1.0)
126 histo_Cnet_output_B0bar = ROOT.TH1F(
'Comb_Net_Output_B0bar',
'Combiner network output [not equal to r] '
127 'for true B0bar (binning 100)', 100, 0.0, 1.0)
129 histo_belleplotB0 = ROOT.TH1F(
'BellePlot_B0',
130 'BellePlot for true B0 (binning 50)', 50,
133 histo_belleplotB0bar = ROOT.TH1F(
'BellePlot_B0Bar',
134 'BellePlot for true B0Bar (binning 50)',
137 histo_notTaggedEvents = ROOT.TH1F(
'notTaggedEvents',
138 'Histogram for not tagged events',
144 histo_calib_B0 = ROOT.TH1F(
'Calibration_B0',
'CalibrationPlot for true B0', 100, -1.0, 1.0)
146 histo_calib_B0bar = ROOT.TH1F(
'Calibration_B0Bar',
147 'CalibrationPlot for true B0Bar', 100, -1.0,
150 hallo12 = ROOT.TH1F(
'BellePlot_NoCut',
'BellePlot_NoCut (binning 100)',
154 diag = ROOT.TF1(
'diag',
'pol1', -1, 1)
158 histo_m0 = ROOT.TH1F(
'BellePlot_m0',
159 'BellePlot_m for true B0 (binning 50)', 50, -1.0, 1.0)
160 histo_m1 = ROOT.TH1F(
'BellePlot_m1',
161 'BellePlot_m for true B0 (binning 50)', 50, -1.0, 1.0)
162 histo_m2 = ROOT.TH1F(
'BellePlot_m2',
163 'BellePlot_m for true B0Bar (binning 50)', 50, -1.0,
168 tree.Draw(method +
'_qrCombined>>BellePlot_B0',
'qrMC == 1')
169 tree.Draw(method +
'_qrCombined>>BellePlot_B0Bar',
'qrMC == -1')
170 tree.Draw(method +
'_qrCombined>>BellePlot_NoCut',
'abs(qrMC) == 1')
172 tree.Draw(method +
'_qrCombined>>Calibration_B0',
'qrMC == 1')
173 tree.Draw(method +
'_qrCombined>>Calibration_B0Bar',
'qrMC == -1')
175 tree.Draw(method +
'_qrCombined>>notTaggedEvents',
176 'abs(qrMC) == 0 && isSignal == 1 && ' +
177 method +
'_qrCombined < -1')
180 tree.Draw(method +
'_qrCombined>>BellePlot_m0',
181 'qrMC == 1 && ' + method +
'_qrCombined>0')
182 tree.Draw(method +
'_qrCombined>>BellePlot_m1',
183 'qrMC == 1 && ' + method +
'_qrCombined<0')
184 tree.Draw(method +
'_qrCombined>>BellePlot_m2',
185 'qrMC == -1 && ' + method +
'_qrCombined>0 ')
190 tree.Project(
'Average_r',
'abs(' + method +
'_qrCombined)',
'abs(' + method +
'_qrCombined)*(abs(qrMC) == 1)')
191 tree.Project(
'Average_rB0',
'abs(' + method +
'_qrCombined)',
'abs(' + method +
'_qrCombined)*(qrMC==1)')
192 tree.Project(
'Average_rB0bar',
'abs(' + method +
'_qrCombined)',
'abs(' + method +
'_qrCombined)*(qrMC==-1)')
194 tree.Project(
'MS_r',
'abs(' + method +
'_qrCombined)',
'(' + method +
195 '_qrCombined*' + method +
'_qrCombined)*(abs(qrMC) == 1)')
196 tree.Project(
'MS_rB0',
'abs(' + method +
'_qrCombined)',
197 '(' + method +
'_qrCombined*' + method +
'_qrCombined)*(qrMC==1)')
198 tree.Project(
'MS_rB0bar',
'abs(' + method +
'_qrCombined)',
199 '(' + method +
'_qrCombined*' + method +
'_qrCombined)*(qrMC==-1)')
202 tree.Project(
'entries_per_bin',
'abs(' + method +
'_qrCombined)',
'abs(qrMC) == 1')
203 tree.Project(
'entries_per_binB0',
'abs(' + method +
'_qrCombined)',
'qrMC == 1')
204 tree.Project(
'entries_per_binB0bar',
'abs(' + method +
'_qrCombined)',
'qrMC == -1')
207 tree.Project(
'mc_NwB0',
'abs(' + method +
'_qrCombined)',
' ' + method +
'_qrCombined*qrMC < 0 && qrMC == 1')
208 tree.Project(
'mc_NwB0bar',
'abs(' + method +
'_qrCombined)',
' ' + method +
'_qrCombined*qrMC < 0 && qrMC == -1')
211 histo_avr_r.Divide(histo_entries_per_bin)
212 histo_avr_rB0.Divide(histo_entries_per_binB0)
213 histo_avr_rB0bar.Divide(histo_entries_per_binB0bar)
215 histo_ms_r.Divide(histo_entries_per_bin)
216 histo_ms_rB0.Divide(histo_entries_per_binB0)
217 histo_ms_rB0bar.Divide(histo_entries_per_binB0bar)
221 histo_calib_B0.Divide(hallo12)
222 histo_calib_B0bar.Divide(hallo12)
226 print(
'****************************** CALIBRATION CHECK FOR COMBINER USING ' +
227 method +
' ***************************************')
229 print(
'Fit ploynomial of first order to the calibration plot. Expected value ~0.5')
231 histo_calib_B0.Fit(diag,
'TEST')
233 print(
'****************************** MEASURED EFFECTIVE EFFICIENCY FOR COMBINER USING ' +
234 method +
' ***********************************')
238 total_tagged = histo_entries_per_bin.GetEntries()
239 total_tagged_B0 = histo_entries_per_binB0.GetEntries()
240 total_tagged_B0bar = histo_entries_per_binB0bar.GetEntries()
241 total_notTagged = histo_notTaggedEvents.GetEntries()
242 total_entries = (total_tagged + total_notTagged)
244 total_entriesB0 = (total_tagged_B0 + total_notTagged / 2)
245 total_entriesB0bar = (total_tagged_B0bar + total_notTagged / 2)
247 tagging_eff = total_tagged / (total_tagged + total_notTagged)
248 DeltaTagging_eff = math.sqrt(total_tagged * total_notTagged**2 + total_notTagged * total_tagged**2) / (total_entries**2)
252 uncertainty_eff_effB0 = 0
253 uncertainty_eff_effB0bar = 0
254 uncertainty_eff_effAverage = 0
255 diff_eff_Uncertainty = 0
256 event_fractionB0 = array(
'f', [0] * r_size)
257 event_fractionB0bar = array(
'f', [0] * r_size)
258 event_fractionTotal = array(
'f', [0] * r_size)
259 event_fractionTotalUncertainty = array(
'f', [0] * r_size)
260 eventsInBin_B0 = array(
'f', [0] * r_size)
261 eventsInBin_B0bar = array(
'f', [0] * r_size)
262 eventsInBin_Total = array(
'f', [0] * r_size)
263 event_fractionDiff = array(
'f', [0] * r_size)
264 event_fractionDiffUncertainty = array(
'f', [0] * r_size)
265 rvalueB0 = array(
'f', [0] * r_size)
266 rvalueB0bar = array(
'f', [0] * r_size)
267 rvalueB0Average = array(
'f', [0] * r_size)
268 rvalueStdB0 = array(
'f', [0] * r_size)
269 rvalueStdB0bar = array(
'f', [0] * r_size)
270 rvalueStdB0Average = array(
'f', [0] * r_size)
271 wvalue = array(
'f', [0] * r_size)
272 wvalueUncertainty = array(
'f', [0] * r_size)
273 wvalueB0 = array(
'f', [0] * r_size)
274 wvalueB0bar = array(
'f', [0] * r_size)
275 wvalueB0Uncertainty = array(
'f', [0] * r_size)
276 wvalueB0barUncertainty = array(
'f', [0] * r_size)
277 wvalueDiff = array(
'f', [0] * r_size)
278 wvalueDiffUncertainty = array(
'f', [0] * r_size)
279 entries = array(
'f', [0] * r_size)
280 entriesB0 = array(
'f', [0] * r_size)
281 entriesB0bar = array(
'f', [0] * r_size)
282 iEffEfficiency = array(
'f', [0] * r_size)
283 iEffEfficiencyUncertainty = array(
'f', [0] * r_size)
284 iEffEfficiencyB0Uncertainty = array(
'f', [0] * r_size)
285 iEffEfficiencyB0barUncertainty = array(
'f', [0] * r_size)
286 iEffEfficiencyB0UncertaintyFromOutput = array(
'f', [0] * r_size)
287 iEffEfficiencyB0barUncertaintyFromOutput = array(
'f', [0] * r_size)
289 iDeltaEffEfficiency = array(
'f', [0] * r_size)
290 iDeltaEffEfficiencyUncertainty = array(
'f', [0] * r_size)
291 muParam = array(
'f', [0] * r_size)
292 muParamUncertainty = array(
'f', [0] * r_size)
295 print(
'* --> DETERMINATION BASED ON MONTE CARLO ' +
299 print(
'* Note: mu = Delta_Effcy/(2*Efficiency). Needed for CP analysis ' +
300 'together with w and Delta_w *')
303 print(
'* ------------------------------------------------------------------' +
304 '-------------------------------------------------- *')
305 print(
'* r-interval <r> Efficiency Delta_Effcy ' +
307 print(
'* ------------------------------------------------------------------' +
308 '-------------------------------------------------- *')
310 for i
in range(1, r_size):
312 entries[i] = histo_entries_per_bin.GetBinContent(i)
313 entriesB0[i] = histo_entries_per_binB0.GetBinContent(i)
314 entriesB0bar[i] = histo_entries_per_binB0bar.GetBinContent(i)
317 event_fractionB0[i] = entriesB0[i] / total_entriesB0
318 event_fractionB0bar[i] = entriesB0bar[i] / total_entriesB0bar
323 event_fractionTotal[i] = (event_fractionB0[i] + event_fractionB0bar[i]) / 2
324 event_fractionDiff[i] = event_fractionB0[i] - event_fractionB0bar[i]
326 event_fractionDiffUncertainty[i] = math.sqrt(entriesB0[i] *
331 (total_entriesB0bar -
333 total_entriesB0bar**3)
335 event_fractionTotalUncertainty[i] = event_fractionDiffUncertainty[i] / 2
337 rvalueB0[i] = histo_avr_rB0.GetBinContent(i)
338 rvalueB0bar[i] = histo_avr_rB0bar.GetBinContent(i)
339 rvalueB0Average[i] = histo_avr_r.GetBinContent(i)
340 rvalueStdB0[i] = math.sqrt(histo_ms_rB0.GetBinContent(
341 i) - (histo_avr_rB0.GetBinContent(i))**2) / math.sqrt(entriesB0[i] - 1)
342 rvalueStdB0bar[i] = math.sqrt(histo_ms_rB0bar.GetBinContent(
343 i) - (histo_avr_rB0bar.GetBinContent(i))**2) / math.sqrt(entriesB0bar[i] - 1)
344 rvalueStdB0Average[i] = math.sqrt(rvalueStdB0[i]**2 + rvalueStdB0bar[i]**2) / 2
348 wvalueB0[i] = histo_mc_NwB0.GetBinContent(i) / entriesB0[i]
349 wvalueB0bar[i] = histo_mc_NwB0bar.GetBinContent(i) / entriesB0bar[i]
350 wvalueDiff[i] = wvalueB0[i] - wvalueB0bar[i]
351 wvalueB0Uncertainty[i] = math.sqrt(histo_mc_NwB0.GetBinContent(
352 i) * (entriesB0[i] - histo_mc_NwB0.GetBinContent(i)) / (entriesB0[i]**3))
353 wvalueB0barUncertainty[i] = math.sqrt(histo_mc_NwB0bar.GetBinContent(
354 i) * (entriesB0bar[i] - histo_mc_NwB0bar.GetBinContent(i)) / (entriesB0bar[i]**3))
356 wvalueDiffUncertainty[i] = math.sqrt(wvalueB0Uncertainty[i]**2 + wvalueB0barUncertainty[i]**2)
357 wvalue[i] = (wvalueB0[i] + wvalueB0bar[i]) / 2
358 wvalueUncertainty[i] = wvalueDiffUncertainty[i] / 2
365 iEffEfficiency[i] = event_fractionTotal[i] * (1 - 2 * wvalue[i])**2
367 iEffEfficiencyUncertainty[i] = (1 - 2 * wvalue[i]) * \
368 math.sqrt((2 * event_fractionTotal[i] * 2 * wvalueUncertainty[i])**2 +
369 (1 - 2 * wvalue[i])**2 * event_fractionTotalUncertainty[i]**2)
392 average_eff_eff += iEffEfficiency[i]
397 iDeltaEffEfficiency[i] = event_fractionB0[i] * (1 - 2 * wvalueB0[i])**2 - \
398 event_fractionB0bar[i] * (1 - 2 * wvalueB0bar[i])**2
400 iEffEfficiencyB0Uncertainty[i] = (1 - 2 * wvalueB0[i]) * \
401 math.sqrt((2 * total_entriesB0 * entriesB0[i] * 2 * wvalueB0Uncertainty[i])**2 +
402 (1 - 2 * wvalueB0[i])**2 * entriesB0[i] *
403 total_entriesB0 * (total_entriesB0 - entriesB0[i])) / (total_entriesB0**2)
404 iEffEfficiencyB0barUncertainty[i] = (1 - 2 * wvalueB0bar[i]) * \
405 math.sqrt((2 * total_entriesB0bar * entriesB0bar[i] * 2 * wvalueB0barUncertainty[i])**2 +
406 (1 - 2 * wvalueB0bar[i])**2 * entriesB0bar[i] *
407 total_entriesB0bar * (total_entriesB0bar - entriesB0bar[i])) / (total_entriesB0bar**2)
409 iDeltaEffEfficiencyUncertainty[i] = math.sqrt(iEffEfficiencyB0Uncertainty[i]**2 + iEffEfficiencyB0barUncertainty[i]**2)
412 diff_eff_Uncertainty = diff_eff_Uncertainty + iDeltaEffEfficiencyUncertainty[i]**2
415 tot_eff_effB0 = tot_eff_effB0 + event_fractionB0[i] * (1 - 2 * wvalueB0[i])**2
416 tot_eff_effB0bar = tot_eff_effB0bar + event_fractionB0bar[i] * (1 - 2 * wvalueB0bar[i])**2
417 uncertainty_eff_effAverage = uncertainty_eff_effAverage + iEffEfficiencyUncertainty[i]**2
418 uncertainty_eff_effB0 = uncertainty_eff_effB0 + iEffEfficiencyB0Uncertainty[i]**2
419 uncertainty_eff_effB0bar = uncertainty_eff_effB0bar + iEffEfficiencyB0barUncertainty[i]**2
420 muParam[i] = event_fractionDiff[i] / (2 * event_fractionTotal[i])
421 muParamUncertainty[i] = event_fractionDiffUncertainty[i] / (2 * event_fractionTotal[i]) * math.sqrt(muParam[i]**2 + 1)
424 print(
'* ' +
'{:.3f}'.format(r_subsample[i - 1]) +
' - ' +
'{:.3f}'.format(r_subsample[i]) +
' ' +
425 '{:.3f}'.format(rvalueB0Average[i]) +
' +- ' +
'{:.4f}'.format(rvalueStdB0Average[i]) +
' ' +
426 '{:.4f}'.format(event_fractionTotal[i]) +
' ' +
427 '{: .4f}'.format(event_fractionDiff[i]) +
' +- ' +
'{:.4f}'.format(event_fractionDiffUncertainty[i]) +
' ' +
428 '{: .4f}'.format(muParam[i]) +
' +- ' +
'{:.4f}'.format(muParamUncertainty[i]) +
' ' +
429 '{:.4f}'.format(wvalue[i]) +
' +- ' +
'{:.4f}'.format(wvalueUncertainty[i]) +
' ' +
430 '{: .4f}'.format(wvalueDiff[i]) +
' +- ' +
'{:.4f}'.format(wvalueDiffUncertainty[i]) +
' *')
433 uncertainty_eff_effAverage = math.sqrt(uncertainty_eff_effAverage)
434 uncertainty_eff_effB0 = math.sqrt(uncertainty_eff_effB0)
435 uncertainty_eff_effB0bar = math.sqrt(uncertainty_eff_effB0bar)
436 diff_eff = tot_eff_effB0 - tot_eff_effB0bar
437 diff_eff_Uncertainty = math.sqrt(diff_eff_Uncertainty)
438 print(
'* --------------------------------------------------------------------------------------------------' +
439 '------------------ *')
441 print(
'* __________________________________________________________________________________________ *')
443 print(
'* | TOTAL NUMBER OF TAGGED EVENTS = ' +
444 '{:<24}'.format(
"%.0f" % total_tagged) +
'{:>36}'.format(
'| *'))
447 '* | TOTAL AVERAGE EFFICIENCY (q=+-1)= ' +
458 '* | TOTAL AVERAGE EFFECTIVE EFFICIENCY (q=+-1)= ' +
464 uncertainty_eff_effAverage *
469 '* | TOTAL AVERAGE EFFECTIVE EFFICIENCY ASYMMETRY (q=+-1)= ' +
475 diff_eff_Uncertainty *
479 print(
'* | B0-TAGGER TOTAL EFFECTIVE EFFICIENCIES: ' +
480 '{:.2f}'.format(tot_eff_effB0 * 100) +
" +-" +
'{: 4.2f}'.format(uncertainty_eff_effB0 * 100) +
482 '{:.2f}'.format(tot_eff_effB0bar * 100) +
" +-" +
'{: 4.2f}'.format(uncertainty_eff_effB0bar * 100) +
483 ' % (q=-1) ' +
' | *')
485 print(
'* | FLAVOR PERCENTAGE (MC): ' +
486 '{:.2f}'.format(total_tagged_B0 / total_tagged * 100) +
' % (q=+1) ' +
487 '{:.2f}'.format(total_tagged_B0bar / total_tagged * 100) +
' % (q=-1) Diff=' +
488 '{:^5.2f}'.format((total_tagged_B0 - total_tagged_B0bar) / total_tagged * 100) +
' % | *')
489 print(
'* |__________________________________________________________________________________________| *')
491 print(
'****************************************************************************************************')
495 print(
'* --------------------------------- *')
496 print(
'* Efficiency Determination - easiest way *')
497 print(
'* --------------------------------- *')
498 total_tagged_B0 = histo_belleplotB0.GetEntries()
499 total_tagged_B0Bar = histo_belleplotB0bar.GetEntries()
500 total_tagged_wrong = histo_m1.GetEntries()
501 total_tagged_B0Bar_wrong = histo_m2.GetEntries()
502 total_tagged = total_tagged_B0 + total_tagged_B0Bar
503 total_tagged_wrong = total_tagged_wrong + total_tagged_B0Bar_wrong
505 wrong_tag_fraction_B0 = total_tagged_wrong / total_tagged_B0
506 wrong_tag_fraction_B0Bar = total_tagged_B0Bar_wrong / total_tagged_B0Bar
507 wrong_tag_fraction = total_tagged_wrong / total_tagged
508 right_tag_fraction_B0 = 1 - 2 * wrong_tag_fraction_B0
509 right_tag_fraction_B0Bar = 1 - 2 * wrong_tag_fraction_B0Bar
510 right_tag_fraction = 1 - 2 * wrong_tag_fraction
511 wrong_eff_B0 = right_tag_fraction_B0 * right_tag_fraction_B0
512 wrong_eff_B0Bar = right_tag_fraction_B0Bar * right_tag_fraction_B0Bar
513 wrong_eff = right_tag_fraction * right_tag_fraction
515 print(
'* wrong_tag_fraction for all: ' +
516 '{:.3f}'.format(wrong_tag_fraction * 100) +
518 print(
'* right_tag_fraction for all: ' +
519 '{:.3f}'.format(right_tag_fraction * 100) +
521 print(
'* wrong calculated eff all: ' +
'{:.3f}'.format(wrong_eff * 100) +
524 print(
'****************************************************************************************************')
526 print(
'Table For B2TIP')
539 maxB0 = histo_belleplotB0.GetBinContent(histo_belleplotB0.GetMaximumBin())
540 maxB0bar = histo_belleplotB0bar.GetBinContent(histo_belleplotB0bar.GetMaximumBin())
542 Ymax = max(maxB0, maxB0bar)
543 Ymax = Ymax + Ymax / 12
545 if YmaxForQrPlot < Ymax:
549 ROOT.gStyle.SetOptStat(0)
550 Canvas1 = ROOT.TCanvas(
'Bla' + method,
'Final Output', 1200, 800)
552 histo_belleplotB0.SetFillColorAlpha(ROOT.kBlue, 0.2)
553 histo_belleplotB0.SetFillStyle(1001)
554 histo_belleplotB0.GetYaxis().SetLabelSize(0.03)
555 histo_belleplotB0.GetYaxis().SetLimits(0, YmaxForQrPlot)
556 histo_belleplotB0.GetYaxis().SetTitleOffset(1.2)
557 histo_belleplotB0.SetLineColor(ROOT.kBlue)
558 histo_belleplotB0bar.SetFillColorAlpha(ROOT.kRed, 1.0)
559 histo_belleplotB0bar.SetFillStyle(3005)
560 histo_belleplotB0bar.SetLineColor(ROOT.kRed)
563 histo_belleplotB0.SetTitle(
'Final Flavor Tagger Output; #it{qr}-output ; Events'
565 histo_belleplotB0.SetMinimum(0)
566 histo_belleplotB0.SetMaximum(YmaxForQrPlot)
567 histo_belleplotB0.Draw(
'hist')
568 histo_belleplotB0bar.Draw(
'hist same')
570 leg = ROOT.TLegend(0.75, 0.8, 0.9, 0.9)
571 leg.AddEntry(histo_belleplotB0,
'true B0')
572 leg.AddEntry(histo_belleplotB0bar,
'true B0bar')
577 with Quiet(ROOT.kError):
578 Canvas1.SaveAs(workingDirectory +
'/' +
'PIC_Belleplot_both' + method +
'.pdf')
581 Canvas2 = ROOT.TCanvas(
'Bla2' + method,
'Calibration plot for true B0', 1200, 800)
583 histo_calib_B0.SetFillColorAlpha(ROOT.kBlue, 0.2)
584 histo_calib_B0.SetFillStyle(1001)
585 histo_calib_B0.GetYaxis().SetTitleOffset(1.2)
586 histo_calib_B0.SetLineColor(ROOT.kBlue)
588 histo_calib_B0.SetTitle(
'Calibration For True B0; #it{qr}-output ; Calibration '
590 histo_calib_B0.Draw(
'hist')
593 with Quiet(ROOT.kError):
594 Canvas2.SaveAs(workingDirectory +
'/' +
'PIC_Calibration_B0' + method +
'.pdf')
597 histo_avr_rB0.Delete()
598 histo_avr_rB0bar.Delete()
600 histo_ms_rB0.Delete()
601 histo_ms_rB0bar.Delete()
602 histo_mc_NwB0.Delete()
603 histo_mc_NwB0bar.Delete()
604 histo_notTaggedEvents.Delete()
605 histo_entries_per_bin.Delete()
606 histo_entries_per_binB0.Delete()
607 histo_entries_per_binB0bar.Delete()
608 histo_Cnet_output_B0.Delete()
609 histo_Cnet_output_B0bar.Delete()
610 histo_belleplotB0.Delete()
611 histo_belleplotB0bar.Delete()
612 histo_calib_B0.Delete()
613 histo_calib_B0bar.Delete()
621 print(
r'\begin{tabularx}{1\textwidth}{@{}r r r r r r r@{}}')
623 print(
r'$r$- Interval $\enskip$ & $\varepsilon_i\ $ & $\Delta\varepsilon_i\ $ & $w_i \pm \delta w_i\enskip\; $ ' +
624 r' & $\Delta w_i \pm \delta\Delta w_i $& $\varepsilon_{\text{eff}, i} \pm \delta\varepsilon_{\text{eff}, i}\enskip\, $ ' +
625 r' & $\Delta \varepsilon_{\text{eff}, i} \pm \delta\Delta \varepsilon_{\text{eff}, i}\enskip\, $\\ \hline\hline')
626 for i
in range(1, r_size):
627 print(
'$ ' +
'{:.3f}'.format(r_subsample[i - 1]) +
' - ' +
'{:.3f}'.format(r_subsample[i]) +
'$ & $'
628 '{: 6.1f}'.format(event_fractionTotal[i] * 100) +
r'$ & $' +
629 '{: 7.2f}'.format(event_fractionDiff[i] * 100) +
r'\;$ & $' +
630 '{: 7.2f}'.format(wvalue[i] * 100) +
r" \pm " +
'{:2.2f}'.format(wvalueUncertainty[i] * 100) +
r'\enskip $ & $' +
631 '{: 7.2f}'.format(wvalueDiff[i] * 100) +
r" \pm " +
632 '{:2.2f}'.format(wvalueDiffUncertainty[i] * 100) +
r'\enskip $ & $' +
633 '{: 8.4f}'.format(iEffEfficiency[i] * 100) +
634 r" \pm " +
'{:2.4f}'.format(iEffEfficiencyUncertainty[i] * 100) +
r'\, $ & $' +
635 '{: 6.4f}'.format(iDeltaEffEfficiency[i] * 100) +
636 r" \pm " +
'{:2.4f}'.format(iDeltaEffEfficiencyUncertainty[i] * 100) +
637 r'\enskip\enskip $ \\ ')
638 print(
r'\hline\hline')
639 print(
r'\multicolumn{1}{r}{Total} & & \multicolumn{5}{r}{ $\varepsilon_\text{eff} = ' +
640 r'\sum_i \varepsilon_i \cdot \langle 1-2w_i\rangle^2 = ' +
641 '{: 6.2f}'.format(average_eff_eff * 100) +
r" \pm " +
'{: 6.2f}'.format(uncertainty_eff_effAverage * 100) +
r'\enskip\, ')
642 print(
r'\Delta \varepsilon_\text{eff} = ' +
643 '{: 6.2f}'.format(diff_eff * 100) +
r" \pm " +
'{: 6.2f}'.format(diff_eff_Uncertainty * 100) +
r'\quad\ $ }' +
646 print(
r'\end{tabular}')
650 print(
'Mu-Values for Table')
653 for i
in range(1, r_size):
654 print(
'$ ' +
'{:.3f}'.format(r_subsample[i - 1]) +
' - ' +
'{:.3f}'.format(r_subsample[i]) +
'$ & $'
655 '{: 7.2f}'.format(muParam[i] * 100) +
r" \pm " +
'{:2.2f}'.format(muParamUncertainty[i] * 100) +
r' $ & ')
696 print(ft.eventLevelParticleLists)
699 print(
'******************************************* MEASURED EFFECTIVE EFFICIENCY FOR INDIVIDUAL CATEGORIES ' +
700 '**********************************************')
705 categoriesPerformance = []
706 NbinsCategories = 100
707 for (particleList, category, combinerVariable)
in ft.eventLevelParticleLists:
709 hist_both = ROOT.TH1F(
'Both_' + category,
'Input Both (B0) ' +
710 category +
' (binning)', NbinsCategories, -1.0, 1.0)
712 hist_signal = ROOT.TH1F(
'Signal_' + category,
'Input Signal (B0) ' +
713 category +
' (binning)', NbinsCategories, -1.0, 1.0)
715 hist_background = ROOT.TH1F(
'Background_' + category,
'Input Background (B0bar) ' +
716 category +
' (binning)', NbinsCategories, -1.0, 1.0)
720 hist_probB0 = ROOT.TH1F(
'Probability' + category,
721 'Transformed to probability (B0) (' + category +
')',
722 NbinsCategories, 0.0, 1.0)
723 hist_probB0bar = ROOT.TH1F(
'ProbabilityB0bar_' + category,
724 'Transformed to probability (B0bar) (' + category +
')',
725 NbinsCategories, 0.0, 1.0)
727 hist_qrB0 = ROOT.TH1F(
'QR' + category,
'Transformed to qp (B0)(' +
728 category +
')', NbinsCategories, -1.0, 1.0)
729 hist_qrB0bar = ROOT.TH1F(
'QRB0bar_' + category,
'Transformed to qp (B0bar) (' +
730 category +
')', NbinsCategories, -1.0, 1.0)
733 histo_entries_per_bin = ROOT.TH1F(
'entries_per_bin_' + category,
'Abs(qp)(B0) (' + category +
')', int(r_size - 2), r_subsample)
734 histo_entries_per_binB0 = ROOT.TH1F(
'entries_per_bin' + category,
'Abs(qp)(B0) (' +
735 category +
')', int(r_size - 2), r_subsample)
736 histo_entries_per_binB0bar = ROOT.TH1F(
'entries_per_binB0bar_' + category,
737 'Abs(qp) (B0bar) (' + category +
')', int(r_size - 2), r_subsample)
741 hist_avr_rB0 = ROOT.TH1F(
'Average_r' + category,
'Average r for B0' +
742 category, int(r_size - 2), r_subsample)
743 hist_avr_rB0bar = ROOT.TH1F(
'Average_rB0bar_' + category,
'Average r for B0bar' +
744 category, int(r_size - 2), r_subsample)
746 hist_ms_rB0 = ROOT.TH1F(
'AverageSqrdR' + category,
'Average r sqrd for B0' +
747 category, int(r_size - 2), r_subsample)
748 hist_ms_rB0bar = ROOT.TH1F(
'AverageSqrdRB0bar_' + category,
'Average r sqrd for B0bar' +
749 category, int(r_size - 2), r_subsample)
753 hist_all = ROOT.TH1F(
'All_' + category,
'Input Signal (B0) and Background (B0Bar)' +
754 category +
' (binning 50)', 50, 0.0, 1.0)
755 hist_calib_B0 = ROOT.TH1F(
'Calib_' + category,
'Calibration Plot for true B0' +
756 category +
' (binning 50)', 50, 0.0, 1.0)
759 if category !=
"KaonNotWeighted" and category !=
"LambdaNotWeighted":
761 tree.Draw(
'qp' + category +
'>>Both_' + category,
'abs(qrMC) == 1.0')
763 tree.Draw(
'qp' + category +
'>>Signal_' + category,
'qrMC == 1.0')
765 tree.Draw(
'qp' + category +
'>>Background_' + category,
'qrMC == -1.0')
766 tree.Draw(
'qp' + category +
'>>All_' + category,
'qrMC!=0')
767 tree.Draw(
'qp' + category +
'>>Calib_' + category,
'qrMC == 1.0')
777 elif category ==
"KaonNotWeighted":
778 tree.Draw(
'extraInfo__boQpOfKaon__bc' +
'>>Both_' + category,
'abs(qrMC) == 1.0')
779 tree.Draw(
'extraInfo__boQpOfKaon__bc' +
'>>Signal_' + category,
'qrMC == 1.0')
780 tree.Draw(
'extraInfo__boQpOfKaon__bc' +
'>>Background_' + category,
'qrMC == -1.0')
781 tree.Draw(
'extraInfo__boQpOfKaon__bc' +
'>>All_' + category,
'qrMC!=0')
782 tree.Draw(
'extraInfo__boQpOfKaon__bc' +
'>>Calib_' + category,
'qrMC == 1.0')
784 elif category ==
"LambdaNotWeighted":
785 tree.Draw(
'extraInfo__boQpOfLambda__bc' +
'>>Both_' + category,
'abs(qrMC) == 1.0')
786 tree.Draw(
'extraInfo__boQpOfLambda__bc' +
'>>Signal_' + category,
'qrMC == 1.0')
787 tree.Draw(
'extraInfo__boQpOfLambda__bc' +
'>>Background_' + category,
'qrMC == -1.0')
788 tree.Draw(
'extraInfo__boQpOfLambda__bc' +
'>>All_' + category,
'qrMC!=0')
789 tree.Draw(
'extraInfo__boQpOfLambda__bc' +
'>>Calib_' + category,
'qrMC == 1.0')
791 hist_calib_B0.Divide(hist_all)
794 maxSignal = hist_signal.GetBinContent(hist_signal.GetMaximumBin())
795 maxBackground = hist_background.GetBinContent(hist_background.GetMaximumBin())
797 Ymax = max(maxSignal, maxBackground)
798 Ymax = Ymax + Ymax / 12
800 ROOT.gStyle.SetOptStat(0)
801 with Quiet(ROOT.kError):
802 Canvas = ROOT.TCanvas(
'Bla',
'TITEL BLA', 1200, 800)
805 hist_signal.SetFillColorAlpha(ROOT.kBlue, 0.2)
806 hist_signal.SetFillStyle(1001)
807 hist_signal.SetTitleSize(0.1)
808 hist_signal.GetXaxis().SetLabelSize(0.04)
809 hist_signal.GetYaxis().SetLabelSize(0.04)
810 hist_signal.GetXaxis().SetTitleSize(0.05)
811 hist_signal.GetYaxis().SetTitleSize(0.05)
812 hist_signal.GetXaxis().SetTitleOffset(0.95)
813 hist_signal.GetYaxis().SetTitleOffset(1.1)
814 hist_signal.GetYaxis().SetLimits(0, Ymax)
815 hist_signal.SetLineColor(ROOT.kBlue)
816 hist_background.SetFillColorAlpha(ROOT.kRed, 1.0)
817 hist_background.SetFillStyle(3005)
818 hist_background.GetYaxis().SetLimits(0, Ymax)
819 hist_background.SetLineColor(ROOT.kRed)
821 hist_signal.SetTitle(category +
' category; #it{qp}-Output ; Events')
823 hist_signal.SetMaximum(Ymax)
825 hist_background.SetMaximum(Ymax)
827 hist_signal.Draw(
'hist')
828 hist_background.Draw(
'hist same')
830 if category ==
'MaximumPstar':
831 l = ROOT.TLegend(0.4, 0.75, 0.6, 0.9)
833 l = ROOT.TLegend(0.6, 0.75, 0.8, 0.9)
834 l.AddEntry(hist_signal,
'true B0')
835 l.AddEntry(hist_background,
'true B0bar')
840 with Quiet(ROOT.kError):
841 Canvas.SaveAs(workingDirectory +
'/' +
'PIC_' + category +
'_Input_Combiner.pdf')
846 binCounter = int(NbinsCategories + 1)
847 dilutionB02 = array(
'd', [0] * binCounter)
848 dilutionB0bar2 = array(
'd', [0] * binCounter)
849 purityB0 = array(
'd', [0] * binCounter)
850 purityB0bar = array(
'd', [0] * binCounter)
851 signal = array(
'd', [0] * binCounter)
852 back = array(
'd', [0] * binCounter)
853 weight = array(
'd', [0] * binCounter)
855 for i
in range(1, binCounter):
857 signal[i] = hist_signal.GetBinContent(i)
858 back[i] = hist_background.GetBinContent(i)
859 weight[i] = signal[i] + back[i]
862 if signal[i] + back[i] == 0:
866 dilutionB0bar2[i] = 0
869 purityB0[i] = signal[i] / (signal[i] + back[i])
870 dilutionB02[i] = -1 + 2 * signal[i] / (signal[i] + back[i])
872 purityB0bar[i] = back[i] / (signal[i] + back[i])
873 dilutionB0bar2[i] = -1 + 2 * back[i] / (signal[i] + back[i])
876 hist_probB0.Fill(purityB0[i], signal[i])
877 hist_probB0bar.Fill(purityB0bar[i], back[i])
880 hist_qrB0.Fill(dilutionB02[i], signal[i])
881 hist_qrB0bar.Fill(dilutionB0bar2[i], back[i])
884 histo_entries_per_binB0.Fill(abs(dilutionB02[i]), signal[i])
885 histo_entries_per_binB0bar.Fill(abs(dilutionB0bar2[i]), back[i])
887 hist_avr_rB0.Fill(abs(dilutionB02[i]), abs(dilutionB02[i]) * signal[i])
888 hist_avr_rB0bar.Fill(abs(dilutionB0bar2[i]), abs(dilutionB0bar2[i]) * back[i])
890 hist_ms_rB0.Fill(abs(dilutionB02[i]), abs(dilutionB02[i] * dilutionB02[i]) * signal[i])
891 hist_ms_rB0bar.Fill(abs(dilutionB0bar2[i]), abs(dilutionB0bar2[i] * dilutionB02[i]) * back[i])
894 hist_avr_rB0.Divide(histo_entries_per_binB0)
895 hist_avr_rB0bar.Divide(histo_entries_per_binB0bar)
897 hist_ms_rB0.Divide(histo_entries_per_binB0)
898 hist_ms_rB0bar.Divide(histo_entries_per_binB0bar)
902 total_entriesB0 = total_notTagged / 2
903 total_entriesB0bar = total_notTagged / 2
904 for i
in range(1, r_size):
905 total_entriesB0 = total_entriesB0 + histo_entries_per_binB0.GetBinContent(i)
906 total_entriesB0bar = total_entriesB0bar + histo_entries_per_binB0bar.GetBinContent(i)
910 uncertainty_eff_effB0 = 0
911 uncertainty_eff_effB0bar = 0
912 uncertainty_eff_effAverage = 0
913 diff_eff_Uncertainty = 0
914 event_fractionB0 = array(
'f', [0] * r_size)
915 event_fractionB0bar = array(
'f', [0] * r_size)
916 rvalueB0 = array(
'f', [0] * r_size)
917 rvalueB0bar = array(
'f', [0] * r_size)
918 rvalueStdB0 = array(
'f', [0] * r_size)
919 rvalueStdB0bar = array(
'f', [0] * r_size)
921 entriesBoth = array(
'f', [0] * r_size)
922 entriesB0 = array(
'f', [0] * r_size)
923 entriesB0bar = array(
'f', [0] * r_size)
924 iEffEfficiencyB0Uncertainty = array(
'f', [0] * r_size)
925 iEffEfficiencyB0barUncertainty = array(
'f', [0] * r_size)
926 iDeltaEffEfficiencyUncertainty = array(
'f', [0] * r_size)
928 for i
in range(1, r_size):
930 entriesBoth[i] = entriesB0bar[i] + entriesB0[i]
931 entriesB0[i] = histo_entries_per_binB0.GetBinContent(i)
932 entriesB0bar[i] = histo_entries_per_binB0bar.GetBinContent(i)
933 event_fractionB0[i] = entriesB0[i] / total_entriesB0
934 event_fractionB0bar[i] = entriesB0bar[i] / total_entriesB0bar
939 rvalueB0[i] = hist_avr_rB0.GetBinContent(i)
940 rvalueB0bar[i] = hist_avr_rB0bar.GetBinContent(i)
943 rvalueStdB0bar[i] = 0
946 rvalueStdB0[i] = math.sqrt(abs(hist_ms_rB0.GetBinContent(
947 i) - (hist_avr_rB0.GetBinContent(i))**2)) / math.sqrt(entriesB0[i] - 1)
949 if entriesB0bar[i] > 1:
950 rvalueStdB0bar[i] = math.sqrt(abs(hist_ms_rB0bar.GetBinContent(
951 i) - (hist_avr_rB0bar.GetBinContent(i))**2)) / math.sqrt(entriesB0bar[i] - 1)
954 tot_eff_effB0 = tot_eff_effB0 + event_fractionB0[i] * rvalueB0[i] \
956 tot_eff_effB0bar = tot_eff_effB0bar + event_fractionB0bar[i] * rvalueB0bar[i] \
959 iEffEfficiencyB0Uncertainty[i] = rvalueB0[i] * \
960 math.sqrt((2 * total_entriesB0 * entriesB0[i] * rvalueStdB0[i])**2 +
961 rvalueB0[i]**2 * entriesB0[i] *
962 (total_entriesB0 * (total_entriesB0 - entriesB0[i]) +
963 entriesB0[i] * total_notTagged)) / (total_entriesB0**2)
964 iEffEfficiencyB0barUncertainty[i] = rvalueB0bar[i] * \
965 math.sqrt((2 * total_entriesB0bar * entriesB0bar[i] * rvalueStdB0bar[i])**2 +
966 rvalueB0bar[i]**2 * entriesB0bar[i] *
967 (total_entriesB0bar * (total_entriesB0bar - entriesB0bar[i]) +
968 entriesB0bar[i] * total_notTagged)) / (total_entriesB0bar**2)
970 iDeltaEffEfficiencyUncertainty[i] = math.sqrt(iEffEfficiencyB0Uncertainty[i]**2 + iEffEfficiencyB0barUncertainty[i]**2)
972 diff_eff_Uncertainty = diff_eff_Uncertainty + iDeltaEffEfficiencyUncertainty[i]**2
974 uncertainty_eff_effB0 = uncertainty_eff_effB0 + iEffEfficiencyB0Uncertainty[i]**2
975 uncertainty_eff_effB0bar = uncertainty_eff_effB0bar + iEffEfficiencyB0barUncertainty[i]**2
977 effDiff = tot_eff_effB0 - tot_eff_effB0bar
978 effAverage = (tot_eff_effB0 + tot_eff_effB0bar) / 2
980 uncertainty_eff_effB0 = math.sqrt(uncertainty_eff_effB0)
981 uncertainty_eff_effB0bar = math.sqrt(uncertainty_eff_effB0bar)
982 diff_eff_Uncertainty = math.sqrt(diff_eff_Uncertainty)
983 uncertainty_eff_effAverage = diff_eff_Uncertainty / 2
985 '{:<25}'.format(
"* " + category) +
' B0-Eff=' +
986 '{: 8.2f}'.format(tot_eff_effB0 * 100) +
" +-" +
'{: 4.2f}'.format(uncertainty_eff_effB0 * 100) +
989 '{: 8.2f}'.format(tot_eff_effB0bar * 100) +
" +-" +
'{: 4.2f}'.format(uncertainty_eff_effB0bar * 100) +
992 '{: 8.2f}'.format(effAverage * 100) +
" +- " +
'{:4.2f}'.format(uncertainty_eff_effAverage * 100) +
' %' +
994 '{: 8.2f}'.format(effDiff * 100) +
" +- " +
'{:4.2f}'.format(diff_eff_Uncertainty * 100) +
' % *')
1008 categoriesPerformance.append((category, effAverage, uncertainty_eff_effAverage, effDiff, diff_eff_Uncertainty))
1009 with Quiet(ROOT.kError):
1018 print(
'**************************************************************************************************************************' +
1019 '************************')
1021 print(
'Table For B2TIP')
1023 print(
r'\begin{tabular}{ l r r }')
1025 print(
r'Categories & $\varepsilon_\text{eff} \pm \delta\varepsilon_\text{eff} $& ' +
1026 r'$\Delta\varepsilon_\text{eff} \pm \delta\Delta\varepsilon_\text{eff}$\\ \hline\hline')
1027 for (category, effAverage, uncertainty_eff_effAverage, effDiff, diff_eff_Uncertainty)
in categoriesPerformance:
1029 '{:<23}'.format(category) +
1031 '{: 6.2f}'.format(effAverage * 100) +
r" \pm " +
'{:4.2f}'.format(uncertainty_eff_effAverage * 100) +
1033 '{: 6.2f}'.format(effDiff * 100) +
r" \pm " +
'{:4.2f}'.format(diff_eff_Uncertainty * 100) +
1036 print(
r'\end{tabular}')
1037 B2INFO(
'qp Output Histograms in pdf format saved at: ' + workingDirectory)
1038 with open(
"categoriesPerformance.pkl",
"wb")
as f:
1039 pickle.dump(categoriesPerformance, f)