Belle II Software development
test6_CPVFlavorTaggerEfficiency.py
1#!/usr/bin/env python3
2
3
10
11"""
12<header>
13 <input>CPVToolsOutput.root</input>
14 <output>test6_CPVFlavorTaggerEfficiency.root</output>
15 <contact>Frank Meier; frank.meier@duke.edu</contact>
16 <description>This file calculates the effective efficiency of the category based flavor tagger considering the two
17 standard combiners and the individual categories. Validation plots are also produced. </description>
18</header>
19"""
20
21import ROOT
22from array import array
23
24ROOT.gROOT.SetBatch(True)
25
26workingFiles = ["../CPVToolsOutput.root"]
27treeName = "B0tree"
28
29
30workingDirectory = '.'
31
32#
33# *******************************************
34# DETERMINATION OF TOTAL EFFECTIVE EFFICIENCY
35# *******************************************
36#
37
38r_subsample = array('d', [
39 0.0,
40 0.1,
41 0.25,
42 0.5,
43 0.625,
44 0.75,
45 0.875,
46 1.0])
47r_size = len(r_subsample)
48average_eff = 0
49
50# All possible Categories
51categories = [
52 'Electron',
53 'IntermediateElectron',
54 'Muon',
55 'IntermediateMuon',
56 'KinLepton',
57 'IntermediateKinLepton',
58 'Kaon',
59 'KaonPion',
60 'SlowPion',
61 'FSC',
62 'MaximumPstar',
63 'FastHadron',
64 'Lambda']
65
66
67methods = []
68
69
70class Quiet:
71 """Context handler class to quiet errors in a 'with' statement"""
72
73 def __init__(self, level=ROOT.kInfo + 1):
74 """Class constructor"""
75
76 self.level = level
77
78 def __enter__(self):
79 """Enter the context"""
80
81 self.oldlevel = ROOT.gErrorIgnoreLevel
82 ROOT.gErrorIgnoreLevel = self.level
83
84 def __exit__(self, type, value, traceback):
85 """Exit the context"""
86 ROOT.gErrorIgnoreLevel = self.oldlevel
87
88
89tree = ROOT.TChain(treeName)
90
91mcstatus = array('d', [-511.5, 0.0, 511.5])
92ROOT.TH1.SetDefaultSumw2()
93
94for iFile in workingFiles:
95 tree.AddFile(iFile)
96
97totalBranches = []
98for branch in tree.GetListOfBranches():
99 totalBranches.append(branch.GetName())
100
101if 'FBDT_qrCombined' in totalBranches:
102 methods.append("FBDT")
103
104usedCategories = []
105for cat in categories:
106 catBranch = 'qp' + cat
107 if catBranch in totalBranches:
108 usedCategories.append(cat)
109
110categoriesNtupleList = ''
111for category in usedCategories:
112 categoriesNtupleList = categoriesNtupleList + f"Eff_{category}:"
113
114
115# Output Validation file
116outputFile = ROOT.TFile("test6_CPVFlavorTaggerEfficiency.root", "RECREATE")
117
118# Values to be watched
119outputNtuple = ROOT.TNtuple(
120 "FT_Efficiencies",
121 "Effective efficiencies of the flavor tagger combiners as well as of the individual tagging categories.",
122 "Eff_FBDT:DeltaEff_FBDT:" + categoriesNtupleList)
123
124outputNtuple.SetAlias('Description', "These are the effective efficiencies of the flavor tagger combiners as well as of " +
125 "the individual tagging efficiencies.")
126outputNtuple.SetAlias(
127 'Check',
128 "These values should not change drastically. Since the nightly reconstruction validation runs" +
129 "on the same input file (which changes only from release to release), the values between builds should be the same.")
130outputNtuple.SetAlias('Contact', "frank.meier@duke.edu")
131
132efficienciesForNtuple = []
133
134YmaxForQrPlot = 0
135
136for method in methods:
137 # histogram contains the average r in each of 7 bins -> calculation see below
138 histo_avr_r = ROOT.TH1F('Average_r', 'Average r in each of 7 bins (B0 and B0bar)', 7,
139 r_subsample)
140 histo_avr_rB0 = ROOT.TH1F('Average_rB0', 'Average r in each of 7 bins (B0)', 7,
141 r_subsample)
142 histo_avr_rB0bar = ROOT.TH1F('Average_rB0bar', 'Average r in each of 7 bins (B0bar)', 7,
143 r_subsample)
144 # histogram with number of entries in for each bin
145 histo_entries_per_bin = ROOT.TH1F(
146 'entries_per_bin',
147 'Events binned in r_subsample according to their r-value for B0 and B0bar prob',
148 7,
149 r_subsample)
150 histo_entries_per_binB0 = ROOT.TH1F('entries_per_binB0', 'Events binned in r_subsample according '
151 'to their r-value for B0 prob', 7, r_subsample)
152 histo_entries_per_binB0bar = ROOT.TH1F('entries_per_binB0bar',
153 'Events binned in r_subsample according to their r-value '
154 'for B0bar prob', 7, r_subsample)
155 # histogram network output (not qr and not r) for true B0 (signal) - not necessary
156 histo_Cnet_output_B0 = ROOT.TH1F('Comb_Net_Output_B0', 'Combiner network output [not equal to r] '
157 'for true B0 (binning 100)', 100, 0.0, 1.0)
158 # histogram network output (not qr and not r) for true B0bar (background) - not necessary
159 histo_Cnet_output_B0bar = ROOT.TH1F('Comb_Net_Output_B0bar', 'Combiner network output [not equal to r] '
160 'for true B0bar (binning 100)', 100, 0.0, 1.0)
161 # histogram containing the belle paper plot (qr-tagger output for true B0)
162 histo_belleplotB0 = ROOT.TH1F('qr_' + method + '_B0',
163 'BellePlot for true B0 (binning 50)', 50,
164 -1.0, 1.0)
165 # histogram containing the belle paper plot (qr-tagger output for true B0bar)
166 histo_belleplotB0bar = ROOT.TH1F('qr_' + method + '_B0Bar',
167 'BellePlot for true B0Bar (binning 50)',
168 50, -1.0, 1.0)
169
170 # histogram containing the qr plot (qr-tagger output)
171 histo_belleplotBoth = ROOT.TH1F('qr_' + method + '_B0Both',
172 'qr-tagger output (binning 50)',
173 50, -1.0, 1.0)
174 # calibration plot for B0. If we get a linear line our MC is fine, than the assumption r ~ 1- 2w is reasonable
175 # expectation is, that for B0 calibration plot: qr=0 half B0 and half B0bar, qr = 1 only B0 and qr = -1
176 # no B0. Inverse for B0bar calibration plot
177 histo_calib_B0 = ROOT.TH1F('Calibration_' + method + '_B0', 'CalibrationPlot for true B0', 100, -1.0, 1.0)
178 # calibration plot for B0bar calibration plot
179 histo_calib_B0bar = ROOT.TH1F('Calibration_' + method + '_B0Bar',
180 'CalibrationPlot for true B0Bar', 100, -1.0,
181 1.0)
182 # belle plot with true B0 and B0bars
183 hallo12 = ROOT.TH1F('BellePlot_NoCut', 'BellePlot_NoCut (binning 100)',
184 100, -1.0, 1.0)
185
186
187 diag = ROOT.TF1('diag', 'pol1', -1, 1)
188
189
190 # histograms for the efficiency calculation in wrong way
191 histo_m0 = ROOT.TH1F('BellePlot_B0_m0',
192 'BellePlot_m for true B0 (binning 50)', 50, -1.0, 1.0)
193 histo_m1 = ROOT.TH1F('BellePlot_B0_m1',
194 'BellePlot_m for true B0 (binning 50)', 50, -1.0, 1.0)
195 histo_m2 = ROOT.TH1F('BellePlot_B0_m2',
196 'BellePlot_m for true B0Bar (binning 50)', 50, -1.0,
197 1.0)
198
199 # filling the histograms
200
201 tree.Draw(method + '_qrCombined>>qr_' + method + '_B0', 'qrMC == 1')
202 tree.Draw(method + '_qrCombined>>qr_' + method + '_B0Bar', 'qrMC == -1')
203 tree.Draw(method + '_qrCombined>>BellePlot_NoCut', 'abs(qrMC) == 1')
204 tree.Draw(method + '_qrCombined>>qr_' + method + '_B0Both', 'abs(qrMC) == 1')
205
206 tree.Draw(method + '_qrCombined>>Calibration_' + method + '_B0', 'qrMC == 1')
207 tree.Draw(method + '_qrCombined>>Calibration_' + method + '_B0Bar', 'qrMC == -1')
208
209 # filling histograms wrong efficiency calculation
210 tree.Draw(method + '_qrCombined>>BellePlot_B0_m0',
211 'qrMC == 1 && ' + method + '_qrCombined>0')
212 tree.Draw(method + '_qrCombined>>BellePlot_B0_m1',
213 'qrMC == 1 && ' + method + '_qrCombined<0')
214 tree.Draw(method + '_qrCombined>>BellePlot_B0_m2',
215 'qrMC == -1 && ' + method + '_qrCombined>0 ')
216
217 # filling with abs(qr) in one of 7 bins with its weight
218 # separate calculation for B0 and B0bar
219
220 tree.Project('Average_r', 'abs(' + method + '_qrCombined)',
221 'abs(' + method + '_qrCombined)')
222 tree.Project('Average_rB0', 'abs(' + method + '_qrCombined)', 'abs(' + method + '_qrCombined)*(qrMC==1)')
223 tree.Project('Average_rB0bar', 'abs(' + method + '_qrCombined)', 'abs(' + method + '_qrCombined)*(qrMC==-1)')
224
225 # filling with abs(qr) in one of 7 bins
226 tree.Project('entries_per_bin', 'abs(' + method + '_qrCombined)', 'abs(qrMC) == 1')
227 tree.Project('entries_per_binB0', 'abs(' + method + '_qrCombined)', 'qrMC == 1')
228 tree.Project('entries_per_binB0bar', 'abs(' + method + '_qrCombined)', 'qrMC == -1')
229
230 # producing the average r histograms
231 histo_avr_r.Divide(histo_entries_per_bin)
232 histo_avr_rB0.Divide(histo_entries_per_binB0)
233 histo_avr_rB0bar.Divide(histo_entries_per_binB0bar)
234
235 # producing the calibration plots
236 # Errors ok
237 histo_calib_B0.Divide(hallo12)
238 histo_calib_B0bar.Divide(hallo12)
239
240 # Fit for calibration plot
241 print(' ')
242 print('****************** CALIBRATION CHECK FOR COMBINER USING ' + method + ' ***************************************')
243 print(' ')
244 print('Fit polynomial of first order to the calibration plot. Expected value ~0.5')
245 print(' ')
246 histo_calib_B0.Fit(diag, 'TEST')
247 print(' ')
248 print('****************** MEASURED EFFECTIVE EFFICIENCY FOR COMBINER USING ' + method + ' ***************************')
249 print('* *')
250 # get total number of entries
251 total_entries = histo_entries_per_bin.GetEntries()
252 total_entries_B0 = histo_entries_per_binB0.GetEntries()
253 total_entries_B0bar = histo_entries_per_binB0bar.GetEntries()
254 tot_eff_effB0 = 0
255 tot_eff_effB0bar = 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_fractionDiff = array('f', [0] * r_size)
260 rvalueB0 = array('f', [0] * r_size)
261 rvalueB0bar = array('f', [0] * r_size)
262 rvalueB0Average = array('f', [0] * r_size)
263 wvalue = array('f', [0] * r_size)
264 wvalueB0 = array('f', [0] * r_size)
265 wvalueB0bar = array('f', [0] * r_size)
266 wvalueDiff = array('f', [0] * r_size)
267 entries = array('f', [0] * r_size)
268 entriesB0 = array('f', [0] * r_size)
269 entriesB0bar = array('f', [0] * r_size)
270 iEffEfficiency = array('f', [0] * r_size)
271 iDeltaEffEfficiency = array('f', [0] * r_size)
272 # intervallEff = array('f', [0] * r_size)
273 performance = []
274 for i in range(1, r_size):
275 # get the average r-value
276 rvalueB0[i] = histo_avr_rB0.GetBinContent(i)
277 rvalueB0bar[i] = histo_avr_rB0bar.GetBinContent(i)
278 rvalueB0Average[i] = (rvalueB0[i] + rvalueB0bar[i]) / 2
279 # calculate the wrong tag fraction (only true if MC data good)
280 wvalue[i] = (1 - rvalueB0Average[i]) / 2
281 wvalueB0[i] = (1 - rvalueB0[i]) / 2
282 wvalueB0bar[i] = (1 - rvalueB0bar[i]) / 2
283 wvalueDiff[i] = wvalueB0[i] - wvalueB0bar[i]
284 entries[i] = histo_entries_per_bin.GetBinContent(i)
285 entriesB0[i] = histo_entries_per_binB0.GetBinContent(i)
286 entriesB0bar[i] = histo_entries_per_binB0bar.GetBinContent(i)
287 # fraction of events/all events
288 event_fractionTotal[i] = (entriesB0[i] + entriesB0bar[i]) / total_entries
289 event_fractionDiff[i] = (entriesB0[i] - entriesB0bar[i]) / total_entries
290 event_fractionB0[i] = entriesB0[i] / total_entries_B0
291 event_fractionB0bar[i] = entriesB0bar[i] / total_entries_B0bar
292 iEffEfficiency[i] = (event_fractionB0[i] * rvalueB0[i] * rvalueB0[i] +
293 event_fractionB0bar[i] * rvalueB0bar[i] * rvalueB0bar[i]) / 2
294 iDeltaEffEfficiency[i] = event_fractionB0[i] * rvalueB0[i] * \
295 rvalueB0[i] - event_fractionB0bar[i] * rvalueB0bar[i] * rvalueB0bar[i]
296 # finally calculating the total effective efficiency
297 tot_eff_effB0 = tot_eff_effB0 + event_fractionB0[i] * rvalueB0[i] * rvalueB0[i]
298 tot_eff_effB0bar = tot_eff_effB0bar + event_fractionB0bar[i] * rvalueB0bar[i] * rvalueB0bar[i]
299
300 average_eff = (tot_eff_effB0 + tot_eff_effB0bar) / 2
301 diff_eff = tot_eff_effB0 - tot_eff_effB0bar
302 print('* ------------------------------------------------------------------------------------------------ *')
303 print('* *')
304 print('* __________________________________________________________________________________________ *')
305 print('* | | *')
306 print('* | TOTAL NUMBER OF TAGGED EVENTS = ' +
307 f"{f'{total_entries:.0f}':<24}" + f"{'| *':>36}")
308 print('* | | *')
309 print('* | TOTAL AVERAGE EFFECTIVE EFFICIENCY (q=+-1)= ' + f'{average_eff * 100:.2f}' +
310 ' % | *')
311 print('* | | *')
312 print('* | B0-TAGGER TOTAL EFFECTIVE EFFICIENCIES: ' +
313 f'{tot_eff_effB0 * 100:.2f}' + ' % (q=+1) ' +
314 f'{tot_eff_effB0bar * 100:.2f}' + ' % (q=-1) EffDiff=' +
315 f'{diff_eff * 100:^5.2f}' + ' % | *')
316 print('* | | *')
317 print('* | FLAVOR PERCENTAGE (MC): ' +
318 f'{total_entries_B0 / total_entries * 100:.2f}' + ' % (q=+1) ' +
319 f'{total_entries_B0bar / total_entries * 100:.2f}' + ' % (q=-1) Diff=' +
320 f'{(total_entries_B0 - total_entries_B0bar) / total_entries * 100:^5.2f}' + ' % | *')
321 print('* |__________________________________________________________________________________________| *')
322 print('* *')
323 print('****************************************************************************************************')
324 print('* *')
325
326 efficienciesForNtuple.append(float(average_eff * 100))
327 efficienciesForNtuple.append(float(diff_eff * 100))
328
329 maxB0 = histo_belleplotB0.GetBinContent(histo_belleplotB0.GetMaximumBin())
330 maxB0bar = histo_belleplotB0bar.GetBinContent(histo_belleplotB0bar.GetMaximumBin())
331 maxB0Both = histo_belleplotBoth.GetBinContent(histo_belleplotBoth.GetMaximumBin())
332
333 Ymax = max(maxB0, maxB0bar, maxB0Both)
334 Ymax = Ymax + Ymax / 12
335
336 if YmaxForQrPlot < Ymax:
337 YmaxForQrPlot = Ymax
338
339 # produce a pdf
340 ROOT.gStyle.SetOptStat(0)
341 with Quiet(ROOT.kError):
342 Canvas1 = ROOT.TCanvas('Bla', 'Final Output', 1200, 800)
343 Canvas1.cd() # activate
344 Canvas1.SetLeftMargin(0.13)
345 Canvas1.SetRightMargin(0.04)
346 Canvas1.SetTopMargin(0.03)
347 Canvas1.SetBottomMargin(0.14)
348 histo_belleplotB0.SetFillColorAlpha(ROOT.kBlue, 0.2)
349 histo_belleplotB0.SetFillStyle(1001)
350 histo_belleplotB0.GetXaxis().SetLabelSize(0.04)
351 histo_belleplotB0.GetYaxis().SetLabelSize(0.04)
352 histo_belleplotB0.GetYaxis().SetTitleOffset(0.9)
353 histo_belleplotB0.GetXaxis().SetTitleSize(0.06)
354 histo_belleplotB0.GetYaxis().SetTitleSize(0.06)
355 histo_belleplotB0.GetYaxis().SetLimits(0, YmaxForQrPlot)
356 histo_belleplotB0.SetLineColor(ROOT.kBlue)
357 histo_belleplotB0bar.SetFillColorAlpha(ROOT.kRed, 1.0)
358 histo_belleplotB0bar.SetFillStyle(3005)
359 histo_belleplotB0bar.SetLineColor(ROOT.kRed)
360 # SetLabelSize etc SetTitle
361
362 histo_belleplotB0.SetTitle('; #it{qr}_{' + method + '} ; Events (Total = ' + f"{f'{total_entries:.0f}':<1}" + ')'
363 )
364 histo_belleplotB0.SetMinimum(0)
365 histo_belleplotB0.SetMaximum(YmaxForQrPlot)
366 histo_belleplotB0.Draw('hist')
367 histo_belleplotB0bar.Draw('hist same')
368
369 leg = ROOT.TLegend(0.2, 0.7, 0.9, 0.95)
370 leg.AddEntry(
371 histo_belleplotB0,
372 'true B^{0} ' +
373 ' #varepsilon_{eff}(B^{0}) = ' +
374 f'{tot_eff_effB0 * 100:.2f}' +
375 '% #frac{n_{B^{0}}}{n} = ' +
376 f'{total_entries_B0 / total_entries * 100:.2f}' +
377 '%')
378 leg.AddEntry(
379 histo_belleplotB0bar,
380 'true #bar{B}^{0} ' +
381 ' #varepsilon_{eff}(#bar{B}^{0}) = ' +
382 f'{tot_eff_effB0bar * 100:.2f}' +
383 '% #frac{n_{#bar{B}^{0}}}{n} = ' +
384 f'{total_entries_B0bar / total_entries * 100:.2f}' +
385 '%')
386 leg.AddEntry("", "Avrg. #bf{ #varepsilon_{eff} = " + f'{average_eff * 100:.2f}' +
387 '%} #Delta#varepsilon_{eff} = ' + f'{diff_eff * 100:^5.2f}' + '%')
388 leg.SetTextSize(0.045)
389 leg.Draw()
390
391 Canvas1.Update()
392
393 with Quiet(ROOT.kError):
394 Canvas1.SaveAs(workingDirectory + '/' + 'test6_CPVFTqr' + method + '_both.pdf')
395
396 # Validation Plot 1
397 histo_belleplotBoth.GetXaxis().SetLabelSize(0.04)
398 histo_belleplotBoth.GetYaxis().SetLabelSize(0.04)
399 histo_belleplotBoth.GetYaxis().SetTitleOffset(0.7)
400 histo_belleplotBoth.GetXaxis().SetTitleOffset(0.7)
401 histo_belleplotBoth.GetXaxis().SetTitleSize(0.06)
402 histo_belleplotBoth.GetYaxis().SetTitleSize(0.06)
403
404 histo_belleplotBoth.GetListOfFunctions().Add(ROOT.TNamed('MetaOptions', 'nostats'))
405 histo_belleplotBoth.GetListOfFunctions().Add(ROOT.TNamed('Description', 'Output of the flavor tagger combiner ' + method))
406 histo_belleplotBoth.GetListOfFunctions().Add(
407 ROOT.TNamed(
408 'Check',
409 'Shape should not change drastically. E.g. Warning if the peak at 0 increases or if the peaks at +-1 decrease.'))
410 histo_belleplotBoth.GetListOfFunctions().Add(ROOT.TNamed('Contact', 'frank.meier@duke.edu'))
411
412 histo_belleplotBoth.SetTitle(
413 'Flavor tagger output for combiner ' +
414 method +
415 '; #it{qr}_{' +
416 method +
417 '} ; Events (Total = ' +
418 f"{f'{total_entries:.0f}':<1}" +
419 ')')
420 histo_belleplotBoth.SetMinimum(0)
421 histo_belleplotBoth.SetMaximum(YmaxForQrPlot)
422 histo_belleplotBoth.SetStats(False)
423 histo_belleplotBoth.Write()
424
425 # Validation Plot 2
426 histo_belleplotB0.GetYaxis().SetTitleOffset(0.7)
427 histo_belleplotB0.GetXaxis().SetTitleOffset(0.7)
428 histo_belleplotB0.SetLineColor(ROOT.kBlue + 2)
429 histo_belleplotB0.SetTitle(
430 'Flavor tagger output for combiner ' +
431 method +
432 ' for true B^{0}s; #it{qr}_{' +
433 method +
434 '} ; Events (Total = ' +
435 f"{f'{histo_belleplotB0.GetEntries():.0f}':<1}" +
436 ')')
437 histo_belleplotB0.SetStats(False)
438
439 histo_belleplotB0.GetListOfFunctions().Add(ROOT.TNamed('MetaOptions', 'nostats'))
440 histo_belleplotB0.GetListOfFunctions().Add(
441 ROOT.TNamed('Description', 'Output of the flavor tagger combiner ' + method + ' for true B0s'))
442 histo_belleplotB0.GetListOfFunctions().Add(
443 ROOT.TNamed(
444 'Check',
445 'Shape should not change drastically. E.g. Warning if the peak at 0 increases or if the peak at +1 decreases.'))
446 histo_belleplotB0.GetListOfFunctions().Add(ROOT.TNamed('Contact', 'frank.meier@duke.edu'))
447 histo_belleplotB0.Write()
448
449 # Validation Plot 3
450 histo_belleplotB0bar.GetXaxis().SetLabelSize(0.04)
451 histo_belleplotB0bar.GetYaxis().SetLabelSize(0.04)
452 histo_belleplotB0bar.GetYaxis().SetTitleOffset(0.7)
453 histo_belleplotB0bar.GetXaxis().SetTitleOffset(0.7)
454 histo_belleplotB0bar.GetXaxis().SetTitleSize(0.06)
455 histo_belleplotB0bar.GetYaxis().SetTitleSize(0.06)
456 histo_belleplotB0bar.SetLineColor(ROOT.kBlue + 2)
457 histo_belleplotB0bar.SetTitle(
458 'Flavor tagger output for combiner ' +
459 method +
460 ' for true #bar{B}^{0}s; #it{qr}_{' +
461 method +
462 '} ; Events (Total = ' +
463 f"{f'{histo_belleplotB0bar.GetEntries():.0f}':<1}" +
464 ')')
465 histo_belleplotB0bar.SetMinimum(0)
466 histo_belleplotB0bar.SetMaximum(YmaxForQrPlot)
467 histo_belleplotB0bar.SetStats(False)
468
469 histo_belleplotB0bar.GetListOfFunctions().Add(ROOT.TNamed('MetaOptions', 'nostats'))
470 histo_belleplotB0bar.GetListOfFunctions().Add(
471 ROOT.TNamed(
472 'Description',
473 'Output of the flavor tagger combiner ' +
474 method +
475 ' for true B0bars'))
476 histo_belleplotB0bar.GetListOfFunctions().Add(ROOT.TNamed(
477 'Check', 'Shape should not change drastically. E.g. Warning if the peak at 0 increases or if the peak at -1 decreases.'))
478 histo_belleplotB0bar.GetListOfFunctions().Add(ROOT.TNamed('Contact', 'frank.meier@duke.edu'))
479 histo_belleplotB0bar.Write()
480
481 # IPython.embed()
482
483 # produce the nice calibration plot
484 with Quiet(ROOT.kError):
485 Canvas2 = ROOT.TCanvas('Bla2', 'Calibration plot for true B0', 1200, 800)
486 Canvas2.cd() # activate
487 Canvas2.SetLeftMargin(0.13)
488 Canvas2.SetRightMargin(0.04)
489 Canvas2.SetTopMargin(0.03)
490 Canvas2.SetBottomMargin(0.14)
491 histo_calib_B0.GetXaxis().SetLabelSize(0.04)
492 histo_calib_B0.GetYaxis().SetLabelSize(0.04)
493 histo_calib_B0.GetYaxis().SetTitleOffset(0.9)
494 histo_calib_B0.GetXaxis().SetTitleSize(0.06)
495 histo_calib_B0.GetYaxis().SetTitleSize(0.06)
496 histo_calib_B0.SetFillColorAlpha(ROOT.kBlue, 0.2)
497 histo_calib_B0.SetFillStyle(1001)
498 histo_calib_B0.GetYaxis().SetTitleOffset(0.9)
499 histo_calib_B0.SetLineColor(ROOT.kBlue)
500
501 histo_calib_B0.SetTitle('; #it{qr}_{' + method + '} ; Calibration '
502 )
503 histo_calib_B0.Draw('hist')
504 diag.Draw('SAME')
505
506 leg2 = ROOT.TLegend(0.2, 0.75, 0.63, 0.93)
507 leg2.SetHeader(" y = #it{m}#it{x} + #it{c}", "")
508 leg2.GetListOfPrimitives().First().SetTextAlign(22)
509 leg2.AddEntry(
510 diag,
511 "#it{m} = " +
512 f"{diag.GetParameter('p1'):.2f}" +
513 " #it{c} = " +
514 f"{diag.GetParameter('p0'):.2f}")
515 leg2.SetTextSize(0.05)
516 leg2.Draw()
517
518 Canvas2.Update()
519 with Quiet(ROOT.kError):
520 Canvas2.SaveAs(workingDirectory + '/' + 'test6_CPVFTCalibration_' + method + '_B0.pdf')
521
522 # Validation Plot 4
523
524 histo_calib_B0.GetYaxis().SetTitleOffset(0.7)
525 histo_calib_B0.GetXaxis().SetTitleOffset(0.7)
526 histo_calib_B0.SetLineColor(ROOT.kBlue + 2)
527 histo_calib_B0.SetTitle('Calibration plot for the flavor tagger combiner ' +
528 method + ' ; #it{qr}_{' + method + '} ; Calibration')
529 histo_calib_B0.SetMinimum(-0.2)
530 histo_calib_B0.SetMaximum(+1.2)
531 histo_calib_B0.SetStats(False)
532
533 histo_calib_B0.GetListOfFunctions().Add(ROOT.TNamed('MetaOptions', 'nostats'))
534 histo_calib_B0.GetListOfFunctions().Add(
535 ROOT.TNamed(
536 'Description',
537 'Calibration plot for the flavor tagger combiner ' +
538 method +
539 ' for true B0s'))
540 histo_calib_B0.GetListOfFunctions().Add(
541 ROOT.TNamed('Check', 'Shape should not change drastically. E.g. warning if the shape stops being linear.'))
542 histo_calib_B0.GetListOfFunctions().Add(ROOT.TNamed('Contact', 'frank.meier@duke.edu'))
543 histo_calib_B0.Write()
544
545 histo_belleplotBoth.Delete()
546 histo_avr_r.Delete()
547 histo_avr_rB0.Delete()
548 histo_avr_rB0bar.Delete()
549 histo_entries_per_bin.Delete()
550 histo_entries_per_binB0.Delete()
551 histo_entries_per_binB0bar.Delete()
552 histo_Cnet_output_B0.Delete()
553 histo_Cnet_output_B0bar.Delete()
554 histo_belleplotB0.Delete()
555 histo_belleplotB0bar.Delete()
556 histo_calib_B0.Delete()
557 histo_calib_B0bar.Delete()
558 hallo12.Delete()
559 histo_m0.Delete()
560 histo_m1.Delete()
561 histo_m2.Delete()
562 Canvas1.Clear()
563 Canvas2.Clear()
564
565
566# **********************************************
567# DETERMINATION OF INDIVIDUAL EFFECTIVE EFFICIENCY
568# **********************************************
569
570print('************************* MEASURED EFFECTIVE EFFICIENCY FOR INDIVIDUAL CATEGORIES *********************************')
571print('* *')
572# input: Classifier input from event-level. Output of event-level is recalculated for input on combiner-level.
573# but is re-evaluated under combiner target. Signal is B0, background is B0Bar.
574
575for category in usedCategories:
576 # histogram of input variable (only signal) - not yet a probability! It's a classifier plot!
577 hist_signal = ROOT.TH1F('Signal_' + category, 'Input Signal (B0)' +
578 category + ' (binning 50)', 50, -1.0, 1.0)
579 # histogram of input variable (only background) - not yet a probability! It's a classifier plot!
580 hist_background = ROOT.TH1F('Background_' + category, 'Input Background (B0bar)' +
581 category + ' (binning 50)', 50, -1.0, 1.0)
582 hist_both = ROOT.TH1F('qp_' + category, 'Input Background (B0bar)' +
583 category + ' (binning 50)', 100, -1, 1)
584
585 # per definition that input is not comparable to the network output, this has to be transformed.
586 # probability output from 0 to 1 (corresponds to net output probability) -> calculation below
587 hist_probB0 = ROOT.TH1F('ProbabilityB0_' + category,
588 'Transformed to probability (B0) (' + category + ')',
589 50, 0.0, 1.0)
590 hist_probB0bar = ROOT.TH1F('ProbabilityB0bar_' + category,
591 'Transformed to probability (B0bar) (' + category + ')',
592 50, 0.0, 1.0)
593 # qp output from -1 to 1 -> transformation below
594 hist_qpB0 = ROOT.TH1F('QRB0_' + category, 'Transformed to qp (B0)(' +
595 category + ')', 50, -1.0, 1.0)
596 hist_qpB0bar = ROOT.TH1F('QRB0bar_' + category, 'Transformed to qp (B0bar) (' +
597 category + ')', 50, -1.0, 1.0)
598 # histogram for abs(qp), i.e. this histogram contains the r-values -> transformation below
599 # also used to get the number of entries, sorted into 7 bins
600 hist_absqpB0 = ROOT.TH1F('AbsQRB0_' + category, 'Abs(qp)(B0) (' + category + ')', 7, r_subsample)
601 hist_absqpB0bar = ROOT.TH1F('AbsQRB0bar_' + category, 'Abs(qp) (B0bar) (' + category + ')', 7, r_subsample)
602 # histogram contains at the end the average r values -> calculation below
603 # sorted into 7 bins
604 hist_aver_rB0 = ROOT.TH1F('AverageRB0_' + category, 'A good one (B0)' +
605 category, 7, r_subsample)
606 hist_aver_rB0bar = ROOT.TH1F('AverageRB0bar_' + category, 'A good one (B0bar)' +
607 category, 7, r_subsample)
608 # ****** TEST OF CALIBRATION ******
609 # for calibration plot we want to have
610 hist_all = ROOT.TH1F('All_' + category, 'Input Signal (B0) and Background (B0Bar)' +
611 category + ' (binning 50)', 50, 0.0, 1.0)
612 tree.Draw('qp' + category + '>>All_' + category, 'qrMC!=0')
613 hist_calib_B0 = ROOT.TH1F('Calib_B0_' + category, 'Calibration Plot for true B0' +
614 category + ' (binning 50)', 50, 0.0, 1.0)
615 tree.Draw('qp' + category + '>>Calib_B0_' + category, 'qrMC == 1.0')
616 hist_calib_B0.Divide(hist_all)
617
618 # fill signal
619 tree.Draw('qp' + category + '>>Signal_' + category, 'qrMC == 1.0')
620 # fill background
621 tree.Draw('qp' + category + '>>Background_' + category, 'qrMC == -1.0'
622 )
623 # fill both
624 tree.Draw('qp' + category + '>>qp_' + category, 'abs(qrMC) == 1.0'
625 )
626
627 # ***** TEST OF CALIBRATION ******
628
629 # initialize some arrays
630 purityB0 = array('d', [0] * 51)
631 dilutionB02 = array('d', [0] * 51)
632 purityB0bar = array('d', [0] * 51)
633 dilutionB0bar2 = array('d', [0] * 51)
634 signal = array('d', [0] * 51)
635 back = array('d', [0] * 51)
636 weight = array('d', [0] * 51)
637
638 for i in range(1, 51):
639 # doing the transformation to probabilities
640 signal[i] = hist_signal.GetBinContent(i)
641 back[i] = hist_background.GetBinContent(i)
642
643 weight[i] = signal[i] + back[i]
644
645 # avoid dividing by zero
646 if signal[i] + back[i] == 0:
647 purityB0[i] = 0
648 dilutionB02[i] = 0
649 purityB0bar[i] = 0
650 dilutionB0bar2[i] = 0
651 else:
652 purityB0[i] = signal[i] / (signal[i] + back[i])
653 dilutionB02[i] = -1 + 2 * signal[i] / (signal[i] + back[i])
654
655 purityB0bar[i] = back[i] / (signal[i] + back[i])
656 dilutionB0bar2[i] = -1 + 2 * back[i] / (signal[i] + back[i])
657
658 # filling histogram with probability from 0 to 1
659 hist_probB0.Fill(purityB0[i], signal[i])
660 hist_probB0bar.Fill(purityB0bar[i], back[i])
661
662 # filling histogram with qp from -1 to 1
663 hist_qpB0.Fill(dilutionB02[i], signal[i])
664 hist_qpB0bar.Fill(dilutionB0bar2[i], back[i])
665
666 # filling histogram with abs(qp), i.e. this histogram contains the r-values (not qp)
667 hist_absqpB0.Fill(abs(dilutionB02[i]), signal[i])
668 hist_absqpB0bar.Fill(abs(dilutionB0bar2[i]), back[i])
669 # filling histogram with abs(qp) special weighted - needed for average r calculation
670 hist_aver_rB0.Fill(abs(dilutionB02[i]), abs(dilutionB02[i]) * signal[i])
671 hist_aver_rB0bar.Fill(abs(dilutionB0bar2[i]), abs(dilutionB0bar2[i]) * back[i])
672
673 # hist_aver_rB0bar contains now the average r-value
674 hist_aver_rB0.Divide(hist_absqpB0)
675 hist_aver_rB0bar.Divide(hist_absqpB0bar)
676 # now calculating the efficiency
677
678 # calculating number of events
679 tot_entriesB0 = 0
680 tot_entriesB0bar = 0
681 for i in range(1, r_size):
682 tot_entriesB0 = tot_entriesB0 + hist_absqpB0.GetBinContent(i)
683 tot_entriesB0bar = tot_entriesB0bar + hist_absqpB0bar.GetBinContent(i)
684 # initializing some arrays
685 tot_eff_effB0 = 0
686 tot_eff_effB0bar = 0
687 event_fractionB0 = array('f', [0] * r_size)
688 event_fractionB0bar = array('f', [0] * r_size)
689 rvalueB0 = array('f', [0] * r_size)
690 rvalueB0bar = array('f', [0] * r_size)
691 # wvalue = array('f', [0] * r_size)
692 entriesB0 = array('f', [0] * r_size)
693 entriesB0bar = array('f', [0] * r_size)
694
695 for i in range(1, r_size):
696 rvalueB0[i] = hist_aver_rB0.GetBinContent(i)
697 rvalueB0bar[i] = hist_aver_rB0bar.GetBinContent(i)
698 # wvalue[i] = (1 - rvalueB0[i]) / 2
699 entriesB0[i] = hist_absqpB0.GetBinContent(i)
700 entriesB0bar[i] = hist_absqpB0bar.GetBinContent(i)
701 event_fractionB0[i] = entriesB0[i] / tot_entriesB0
702 event_fractionB0bar[i] = entriesB0bar[i] / tot_entriesB0bar
703 # print '* Bin ' + str(i) + ' r-value: ' + str(rvalueB0[i]), 'entriesB0: ' +
704 # str(event_fractionB0[i] * 100) + ' % (' + str(entriesB0[i]) + '/' +
705 # str(tot_entriesB0) + ')'
706 tot_eff_effB0 = tot_eff_effB0 + event_fractionB0[i] * rvalueB0[i] \
707 * rvalueB0[i]
708 tot_eff_effB0bar = tot_eff_effB0bar + event_fractionB0bar[i] * rvalueB0bar[i] \
709 * rvalueB0bar[i]
710 effDiff = tot_eff_effB0 - tot_eff_effB0bar
711 effAverage = (tot_eff_effB0 + tot_eff_effB0bar) / 2
712
713 print(
714 f"{'* ' + category:<25}" +
715 ' B0-Eff=' +
716 f'{tot_eff_effB0 * 100: 8.2f}' +
717 ' %' +
718 ' B0bar-Eff=' +
719 f'{tot_eff_effB0bar * 100: 8.2f}' +
720 ' %' +
721 ' EffAverage=' +
722 f'{effAverage * 100: 8.2f}' + ' %' +
723 ' EffDiff=' +
724 f'{effDiff * 100: 8.2f}' + ' % *')
725
726 # ****** produce the input plots from combiner level ******
727
728 efficienciesForNtuple.append(float(effAverage * 100))
729
730 maxSignal = hist_signal.GetBinContent(hist_signal.GetMaximumBin())
731 maxBackground = hist_background.GetBinContent(hist_background.GetMaximumBin())
732
733 Ymax = max(maxSignal, maxBackground)
734 Ymax = Ymax + Ymax / 12
735
736 ROOT.gStyle.SetOptStat(0)
737 with Quiet(ROOT.kError):
738 Canvas = ROOT.TCanvas('Bla', 'TITEL BLA', 1200, 800)
739 Canvas.cd() # activate
740 Canvas.SetLogy()
741 Canvas.SetLeftMargin(0.13)
742 Canvas.SetRightMargin(0.04)
743 Canvas.SetTopMargin(0.03)
744 Canvas.SetBottomMargin(0.14)
745 hist_signal.SetFillColorAlpha(ROOT.kBlue, 0.2)
746 hist_signal.SetFillStyle(1001)
747 hist_signal.SetTitleSize(0.1)
748 hist_signal.GetXaxis().SetLabelSize(0.04)
749 hist_signal.GetYaxis().SetLabelSize(0.04)
750 hist_signal.GetXaxis().SetTitleSize(0.06)
751 hist_signal.GetYaxis().SetTitleSize(0.06)
752 hist_signal.GetXaxis().SetLabelSize(0.04)
753 hist_signal.GetYaxis().SetLabelSize(0.04)
754 hist_signal.GetXaxis().SetTitleSize(0.05)
755 hist_signal.GetYaxis().SetTitleSize(0.05)
756 hist_signal.GetXaxis().SetTitleOffset(0.95)
757 hist_signal.GetYaxis().SetTitleOffset(1.1)
758 hist_signal.GetXaxis().SetTitleOffset(1.15)
759 hist_signal.GetYaxis().SetLimits(0, Ymax)
760 hist_signal.SetLineColor(ROOT.kBlue)
761 hist_background.SetFillColorAlpha(ROOT.kRed, 1.0)
762 hist_background.SetFillStyle(3005)
763 hist_background.GetYaxis().SetLimits(0, Ymax)
764 hist_background.SetLineColor(ROOT.kRed)
765
766 catName = category
767 if category == 'MaximumPstar':
768 catName = 'MaximumP*'
769
770 hist_signal.SetTitle('; (#it{qp})^{' + catName + '} ; Events')
771 # hist_signal.SetMinimum(0)
772 hist_signal.SetMaximum(Ymax)
773 # hist_background.SetMinimum(0)
774 hist_background.SetMaximum(Ymax)
775
776 hist_signal.Draw('hist')
777 hist_background.Draw('hist same')
778
779 l0 = ROOT.TLegend(0.13, 0.65, 0.33, 0.97)
780 l0.SetFillColorAlpha(ROOT.kWhite, 0)
781 l0.AddEntry(hist_signal, ' #varepsilon_{eff}(B^{0}) = ' + f'{tot_eff_effB0 * 100:.2f}' + "%")
782 l0.AddEntry(hist_background, ' #varepsilon_{eff}(#bar{B}^{0}) = ' + f'{tot_eff_effB0bar * 100:.2f}' + "%")
783 l0.AddEntry("", "#bf{#varepsilon_{eff} = " + f'{effAverage * 100:.2f}' + '%}')
784 l0.AddEntry("", '#Delta#varepsilon_{eff} = ' + f'{effDiff * 100:^5.2f}' + '%')
785 l0.SetBorderSize(0)
786 l0.SetTextSize(0.045)
787 l0.Draw()
788
789 l1 = ROOT.TLegend(0.85, 0.7, 0.96, 0.97)
790 l1.SetFillColorAlpha(ROOT.kWhite, 0.35)
791 l1.AddEntry(hist_signal, 'B^{0}_{MC}')
792 l1.AddEntry(hist_background, '#bar{B}^{0}_{MC}')
793 l1.SetTextSize(0.045)
794 l1.Draw()
795
796 Canvas.Update()
797 with Quiet(ROOT.kError):
798 Canvas.SaveAs(workingDirectory + '/' + 'test6_CPVFTqp_' + category + '_both.pdf')
799
800 # Validation Plot 4
801 hist_both.GetXaxis().SetLabelSize(0.04)
802 hist_both.GetYaxis().SetLabelSize(0.04)
803 hist_both.GetYaxis().SetTitleOffset(0.7)
804 hist_both.GetXaxis().SetTitleOffset(0.7)
805 hist_both.GetXaxis().SetTitleSize(0.06)
806 hist_both.GetYaxis().SetTitleSize(0.06)
807
808 hist_both.GetListOfFunctions().Add(ROOT.TNamed('MetaOptions', 'nostats, logy'))
809 hist_both.GetListOfFunctions().Add(ROOT.TNamed('Description', 'Output of the flavor tagger category ' + catName))
810 hist_both.GetListOfFunctions().Add(
811 ROOT.TNamed('Check', 'Shape should not change drastically. E.g. Warning if there is only a peak at 0.'))
812 hist_both.GetListOfFunctions().Add(ROOT.TNamed('Contact', 'frank.meier@duke.edu'))
813
814 hist_both.SetTitle(
815 'Flavor tagger output of the category ' +
816 catName +
817 '; #it{qp}_{' +
818 catName +
819 '} ; Events (Total = ' +
820 f"{f'{hist_both.GetEntries():.0f}':<1}" +
821 ')')
822 # hist_both.SetStats(False)
823 hist_both.Write()
824
825 Canvas.Clear()
826
827outputNtuple.Fill(array('f', efficienciesForNtuple))
828outputNtuple.Write()
829outputFile.Close()
830
831print('* *')
832print('*******************************************************************************************************************')
test6_CPVFlavorTaggerEfficiency.Quiet
Definition: test6_CPVFlavorTaggerEfficiency.py:70
test6_CPVFlavorTaggerEfficiency.Quiet.oldlevel
oldlevel
the previously set level to be ignored
Definition: test6_CPVFlavorTaggerEfficiency.py:81
test6_CPVFlavorTaggerEfficiency.Quiet.__enter__
def __enter__(self)
Definition: test6_CPVFlavorTaggerEfficiency.py:78
test6_CPVFlavorTaggerEfficiency.Quiet.__exit__
def __exit__(self, type, value, traceback)
Definition: test6_CPVFlavorTaggerEfficiency.py:84
test6_CPVFlavorTaggerEfficiency.Quiet.__init__
def __init__(self, level=ROOT.kInfo+1)
Definition: test6_CPVFlavorTaggerEfficiency.py:73
test6_CPVFlavorTaggerEfficiency.Quiet.level
level
the level to quiet
Definition: test6_CPVFlavorTaggerEfficiency.py:76