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