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