Belle II Software  release-08-01-10
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>Yo Sato; yosato@post.kek.jp</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 
21 import ROOT
22 from array import array
23 
24 ROOT.gROOT.SetBatch(True)
25 
26 workingFiles = ["../CPVToolsOutput.root"]
27 treeName = "B0tree"
28 
29 
30 workingDirectory = '.'
31 
32 #
33 # *******************************************
34 # DETERMINATION OF TOTAL EFFECTIVE EFFICIENCY
35 # *******************************************
36 #
37 
38 r_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])
47 r_size = len(r_subsample)
48 average_eff = 0
49 
50 # All possible Categories
51 categories = [
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 
67 methods = []
68 
69 
70 class 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.levellevel = level
77 
78  def __enter__(self):
79  """Enter the context"""
80 
81  self.oldleveloldlevel = ROOT.gErrorIgnoreLevel
82  ROOT.gErrorIgnoreLevel = self.levellevel
83 
84  def __exit__(self, type, value, traceback):
85  """Exit the context"""
86  ROOT.gErrorIgnoreLevel = self.oldleveloldlevel
87 
88 
89 tree = ROOT.TChain(treeName)
90 
91 mcstatus = array('d', [-511.5, 0.0, 511.5])
92 ROOT.TH1.SetDefaultSumw2()
93 
94 for iFile in workingFiles:
95  tree.AddFile(iFile)
96 
97 totalBranches = []
98 for branch in tree.GetListOfBranches():
99  totalBranches.append(branch.GetName())
100 
101 if 'FBDT_qrCombined' in totalBranches:
102  methods.append("FBDT")
103 
104 usedCategories = []
105 for cat in categories:
106  catBranch = 'qp' + cat
107  if catBranch in totalBranches:
108  usedCategories.append(cat)
109 
110 categoriesNtupleList = ''
111 for category in usedCategories:
112  categoriesNtupleList = categoriesNtupleList + "Eff_%s:" % category
113 
114 
115 # Output Validation file
116 outputFile = ROOT.TFile("test6_CPVFlavorTaggerEfficiency.root", "RECREATE")
117 
118 # Values to be watched
119 outputNtuple = 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 
124 outputNtuple.SetAlias('Description', "These are the effective efficiencies of the flavor tagger combiners as well as of " +
125  "the individual tagging efficiencies.")
126 outputNtuple.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.")
130 outputNtuple.SetAlias('Contact', "yosato@post.kek.jp")
131 
132 efficienciesForNtuple = []
133 
134 YmaxForQrPlot = 0
135 
136 for 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  '{:<24}'.format("%.0f" % total_entries) + '{:>36}'.format('| *'))
308  print('* | | *')
309  print('* | TOTAL AVERAGE EFFECTIVE EFFICIENCY (q=+-1)= ' + '{:.2f}'.format(average_eff * 100) +
310  ' % | *')
311  print('* | | *')
312  print('* | B0-TAGGER TOTAL EFFECTIVE EFFICIENCIES: ' +
313  '{:.2f}'.format(tot_eff_effB0 * 100) + ' % (q=+1) ' +
314  '{:.2f}'.format(tot_eff_effB0bar * 100) + ' % (q=-1) EffDiff=' +
315  '{:^5.2f}'.format(diff_eff * 100) + ' % | *')
316  print('* | | *')
317  print('* | FLAVOR PERCENTAGE (MC): ' +
318  '{:.2f}'.format(total_entries_B0 / total_entries * 100) + ' % (q=+1) ' +
319  '{:.2f}'.format(total_entries_B0bar / total_entries * 100) + ' % (q=-1) Diff=' +
320  '{:^5.2f}'.format((total_entries_B0 - total_entries_B0bar) / total_entries * 100) + ' % | *')
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 = ' + '{:<1}'.format("%.0f" % total_entries) + ')'
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  '{:.2f}'.format(
375  tot_eff_effB0 *
376  100) +
377  '% #frac{n_{B^{0}}}{n} = ' +
378  '{:.2f}'.format(
379  total_entries_B0 /
380  total_entries *
381  100) +
382  '%')
383  leg.AddEntry(
384  histo_belleplotB0bar,
385  'true #bar{B}^{0} ' +
386  ' #varepsilon_{eff}(#bar{B}^{0}) = ' +
387  '{:.2f}'.format(
388  tot_eff_effB0bar *
389  100) +
390  '% #frac{n_{#bar{B}^{0}}}{n} = ' +
391  '{:.2f}'.format(
392  total_entries_B0bar /
393  total_entries *
394  100) +
395  '%')
396  leg.AddEntry("", "Avrg. #bf{ #varepsilon_{eff} = " + '{:.2f}'.format(average_eff * 100) +
397  '%} #Delta#varepsilon_{eff} = ' + '{:^5.2f}'.format(diff_eff * 100) + '%')
398  leg.SetTextSize(0.045)
399  leg.Draw()
400 
401  Canvas1.Update()
402 
403  with Quiet(ROOT.kError):
404  Canvas1.SaveAs(workingDirectory + '/' + 'test6_CPVFTqr' + method + '_both.pdf')
405 
406  # Validation Plot 1
407  histo_belleplotBoth.GetXaxis().SetLabelSize(0.04)
408  histo_belleplotBoth.GetYaxis().SetLabelSize(0.04)
409  histo_belleplotBoth.GetYaxis().SetTitleOffset(0.7)
410  histo_belleplotBoth.GetXaxis().SetTitleOffset(0.7)
411  histo_belleplotBoth.GetXaxis().SetTitleSize(0.06)
412  histo_belleplotBoth.GetYaxis().SetTitleSize(0.06)
413 
414  histo_belleplotBoth.GetListOfFunctions().Add(ROOT.TNamed('MetaOptions', 'nostats'))
415  histo_belleplotBoth.GetListOfFunctions().Add(ROOT.TNamed('Description', 'Output of the flavor tagger combiner ' + method))
416  histo_belleplotBoth.GetListOfFunctions().Add(
417  ROOT.TNamed(
418  'Check',
419  'Shape should not change drastically. E.g. Warning if the peak at 0 increases or if the peaks at +-1 decrease.'))
420  histo_belleplotBoth.GetListOfFunctions().Add(ROOT.TNamed('Contact', 'yosato@post.kek.jp'))
421 
422  histo_belleplotBoth.SetTitle(
423  'Flavor tagger output for combiner ' +
424  method +
425  '; #it{qr}_{' +
426  method +
427  '} ; Events (Total = ' +
428  '{:<1}'.format(
429  "%.0f" %
430  total_entries) +
431  ')')
432  histo_belleplotBoth.SetMinimum(0)
433  histo_belleplotBoth.SetMaximum(YmaxForQrPlot)
434  histo_belleplotBoth.SetStats(False)
435  histo_belleplotBoth.Write()
436 
437  # Validation Plot 2
438  histo_belleplotB0.GetYaxis().SetTitleOffset(0.7)
439  histo_belleplotB0.GetXaxis().SetTitleOffset(0.7)
440  histo_belleplotB0.SetLineColor(ROOT.kBlue + 2)
441  histo_belleplotB0.SetTitle(
442  'Flavor tagger output for combiner ' +
443  method +
444  ' for true B^{0}s; #it{qr}_{' +
445  method +
446  '} ; Events (Total = ' +
447  '{:<1}'.format(
448  "%.0f" %
449  histo_belleplotB0.GetEntries()) +
450  ')')
451  histo_belleplotB0.SetStats(False)
452 
453  histo_belleplotB0.GetListOfFunctions().Add(ROOT.TNamed('MetaOptions', 'nostats'))
454  histo_belleplotB0.GetListOfFunctions().Add(
455  ROOT.TNamed('Description', 'Output of the flavor tagger combiner ' + method + ' for true B0s'))
456  histo_belleplotB0.GetListOfFunctions().Add(
457  ROOT.TNamed(
458  'Check',
459  'Shape should not change drastically. E.g. Warning if the peak at 0 increases or if the peak at +1 decreases.'))
460  histo_belleplotB0.GetListOfFunctions().Add(ROOT.TNamed('Contact', 'yosato@post.kek.jp'))
461  histo_belleplotB0.Write()
462 
463  # Validation Plot 3
464  histo_belleplotB0bar.GetXaxis().SetLabelSize(0.04)
465  histo_belleplotB0bar.GetYaxis().SetLabelSize(0.04)
466  histo_belleplotB0bar.GetYaxis().SetTitleOffset(0.7)
467  histo_belleplotB0bar.GetXaxis().SetTitleOffset(0.7)
468  histo_belleplotB0bar.GetXaxis().SetTitleSize(0.06)
469  histo_belleplotB0bar.GetYaxis().SetTitleSize(0.06)
470  histo_belleplotB0bar.SetLineColor(ROOT.kBlue + 2)
471  histo_belleplotB0bar.SetTitle(
472  'Flavor tagger output for combiner ' +
473  method +
474  ' for true #bar{B}^{0}s; #it{qr}_{' +
475  method +
476  '} ; Events (Total = ' +
477  '{:<1}'.format(
478  "%.0f" %
479  histo_belleplotB0bar.GetEntries()) +
480  ')')
481  histo_belleplotB0bar.SetMinimum(0)
482  histo_belleplotB0bar.SetMaximum(YmaxForQrPlot)
483  histo_belleplotB0bar.SetStats(False)
484 
485  histo_belleplotB0bar.GetListOfFunctions().Add(ROOT.TNamed('MetaOptions', 'nostats'))
486  histo_belleplotB0bar.GetListOfFunctions().Add(
487  ROOT.TNamed(
488  'Description',
489  'Output of the flavor tagger combiner ' +
490  method +
491  ' for true B0bars'))
492  histo_belleplotB0bar.GetListOfFunctions().Add(ROOT.TNamed(
493  'Check', 'Shape should not change drastically. E.g. Warning if the peak at 0 increases or if the peak at -1 decreases.'))
494  histo_belleplotB0bar.GetListOfFunctions().Add(ROOT.TNamed('Contact', 'yosato@post.kek.jp'))
495  histo_belleplotB0bar.Write()
496 
497  # IPython.embed()
498 
499  # produce the nice calibration plot
500  with Quiet(ROOT.kError):
501  Canvas2 = ROOT.TCanvas('Bla2', 'Calibration plot for true B0', 1200, 800)
502  Canvas2.cd() # activate
503  Canvas2.SetLeftMargin(0.13)
504  Canvas2.SetRightMargin(0.04)
505  Canvas2.SetTopMargin(0.03)
506  Canvas2.SetBottomMargin(0.14)
507  histo_calib_B0.GetXaxis().SetLabelSize(0.04)
508  histo_calib_B0.GetYaxis().SetLabelSize(0.04)
509  histo_calib_B0.GetYaxis().SetTitleOffset(0.9)
510  histo_calib_B0.GetXaxis().SetTitleSize(0.06)
511  histo_calib_B0.GetYaxis().SetTitleSize(0.06)
512  histo_calib_B0.SetFillColorAlpha(ROOT.kBlue, 0.2)
513  histo_calib_B0.SetFillStyle(1001)
514  histo_calib_B0.GetYaxis().SetTitleOffset(0.9)
515  histo_calib_B0.SetLineColor(ROOT.kBlue)
516 
517  histo_calib_B0.SetTitle('; #it{qr}_{' + method + '} ; Calibration '
518  )
519  histo_calib_B0.Draw('hist')
520  diag.Draw('SAME')
521 
522  leg2 = ROOT.TLegend(0.2, 0.75, 0.63, 0.93)
523  leg2.SetHeader(" y = #it{m}#it{x} + #it{c}", "")
524  leg2.GetListOfPrimitives().First().SetTextAlign(22)
525  leg2.AddEntry(
526  diag,
527  "#it{m} = " +
528  '{:.2f}'.format(
529  diag.GetParameter("p1")) +
530  " #it{c} = " +
531  '{:.2f}'.format(
532  diag.GetParameter("p0")))
533  leg2.SetTextSize(0.05)
534  leg2.Draw()
535 
536  Canvas2.Update()
537  with Quiet(ROOT.kError):
538  Canvas2.SaveAs(workingDirectory + '/' + 'test6_CPVFTCalibration_' + method + '_B0.pdf')
539 
540  # Validation Plot 4
541 
542  histo_calib_B0.GetYaxis().SetTitleOffset(0.7)
543  histo_calib_B0.GetXaxis().SetTitleOffset(0.7)
544  histo_calib_B0.SetLineColor(ROOT.kBlue + 2)
545  histo_calib_B0.SetTitle('Calibration plot for the flavor tagger combiner ' +
546  method + ' ; #it{qr}_{' + method + '} ; Calibration')
547  histo_calib_B0.SetMinimum(-0.2)
548  histo_calib_B0.SetMaximum(+1.2)
549  histo_calib_B0.SetStats(False)
550 
551  histo_calib_B0.GetListOfFunctions().Add(ROOT.TNamed('MetaOptions', 'nostats'))
552  histo_calib_B0.GetListOfFunctions().Add(
553  ROOT.TNamed(
554  'Description',
555  'Calibration plot for the flavor tagger combiner ' +
556  method +
557  ' for true B0s'))
558  histo_calib_B0.GetListOfFunctions().Add(
559  ROOT.TNamed('Check', 'Shape should not change drastically. E.g. warning if the shape stops being linear.'))
560  histo_calib_B0.GetListOfFunctions().Add(ROOT.TNamed('Contact', 'yosato@post.kek.jp'))
561  histo_calib_B0.Write()
562 
563  histo_belleplotBoth.Delete()
564  histo_avr_r.Delete()
565  histo_avr_rB0.Delete()
566  histo_avr_rB0bar.Delete()
567  histo_entries_per_bin.Delete()
568  histo_entries_per_binB0.Delete()
569  histo_entries_per_binB0bar.Delete()
570  histo_Cnet_output_B0.Delete()
571  histo_Cnet_output_B0bar.Delete()
572  histo_belleplotB0.Delete()
573  histo_belleplotB0bar.Delete()
574  histo_calib_B0.Delete()
575  histo_calib_B0bar.Delete()
576  hallo12.Delete()
577  histo_m0.Delete()
578  histo_m1.Delete()
579  histo_m2.Delete()
580  Canvas1.Clear()
581  Canvas2.Clear()
582 
583 
584 # **********************************************
585 # DETERMINATION OF INDIVIDUAL EFFECTIVE EFFICIENCY
586 # **********************************************
587 
588 print('************************* MEASURED EFFECTIVE EFFICIENCY FOR INDIVIDUAL CATEGORIES *********************************')
589 print('* *')
590 # input: Classifier input from event-level. Output of event-level is recalculated for input on combiner-level.
591 # but is re-evaluated under combiner target. Signal is B0, background is B0Bar.
592 
593 for category in usedCategories:
594  # histogram of input variable (only signal) - not yet a probability! It's a classifier plot!
595  hist_signal = ROOT.TH1F('Signal_' + category, 'Input Signal (B0)' +
596  category + ' (binning 50)', 50, -1.0, 1.0)
597  # histogram of input variable (only background) - not yet a probability! It's a classifier plot!
598  hist_background = ROOT.TH1F('Background_' + category, 'Input Background (B0bar)' +
599  category + ' (binning 50)', 50, -1.0, 1.0)
600  hist_both = ROOT.TH1F('qp_' + category, 'Input Background (B0bar)' +
601  category + ' (binning 50)', 100, -1, 1)
602 
603  # per definition that input is not comparable to the network output, this has to be transformed.
604  # probability output from 0 to 1 (corresponds to net output probability) -> calculation below
605  hist_probB0 = ROOT.TH1F('ProbabilityB0_' + category,
606  'Transformed to probability (B0) (' + category + ')',
607  50, 0.0, 1.0)
608  hist_probB0bar = ROOT.TH1F('ProbabilityB0bar_' + category,
609  'Transformed to probability (B0bar) (' + category + ')',
610  50, 0.0, 1.0)
611  # qp output from -1 to 1 -> transformation below
612  hist_qpB0 = ROOT.TH1F('QRB0_' + category, 'Transformed to qp (B0)(' +
613  category + ')', 50, -1.0, 1.0)
614  hist_qpB0bar = ROOT.TH1F('QRB0bar_' + category, 'Transformed to qp (B0bar) (' +
615  category + ')', 50, -1.0, 1.0)
616  # histogram for abs(qp), i.e. this histogram contains the r-values -> transformation below
617  # also used to get the number of entries, sorted into 7 bins
618  hist_absqpB0 = ROOT.TH1F('AbsQRB0_' + category, 'Abs(qp)(B0) (' + category + ')', 7, r_subsample)
619  hist_absqpB0bar = ROOT.TH1F('AbsQRB0bar_' + category, 'Abs(qp) (B0bar) (' + category + ')', 7, r_subsample)
620  # histogram contains at the end the average r values -> calculation below
621  # sorted into 7 bins
622  hist_aver_rB0 = ROOT.TH1F('AverageRB0_' + category, 'A good one (B0)' +
623  category, 7, r_subsample)
624  hist_aver_rB0bar = ROOT.TH1F('AverageRB0bar_' + category, 'A good one (B0bar)' +
625  category, 7, r_subsample)
626  # ****** TEST OF CALIBRATION ******
627  # for calibration plot we want to have
628  hist_all = ROOT.TH1F('All_' + category, 'Input Signal (B0) and Background (B0Bar)' +
629  category + ' (binning 50)', 50, 0.0, 1.0)
630  tree.Draw('qp' + category + '>>All_' + category, 'qrMC!=0')
631  hist_calib_B0 = ROOT.TH1F('Calib_B0_' + category, 'Calibration Plot for true B0' +
632  category + ' (binning 50)', 50, 0.0, 1.0)
633  tree.Draw('qp' + category + '>>Calib_B0_' + category, 'qrMC == 1.0')
634  hist_calib_B0.Divide(hist_all)
635 
636  # fill signal
637  tree.Draw('qp' + category + '>>Signal_' + category, 'qrMC == 1.0')
638  # fill background
639  tree.Draw('qp' + category + '>>Background_' + category, 'qrMC == -1.0'
640  )
641  # fill both
642  tree.Draw('qp' + category + '>>qp_' + category, 'abs(qrMC) == 1.0'
643  )
644 
645  # ***** TEST OF CALIBRATION ******
646 
647  # initialize some arrays
648  purityB0 = array('d', [0] * 51)
649  dilutionB02 = array('d', [0] * 51)
650  purityB0bar = array('d', [0] * 51)
651  dilutionB0bar2 = array('d', [0] * 51)
652  signal = array('d', [0] * 51)
653  back = array('d', [0] * 51)
654  weight = array('d', [0] * 51)
655 
656  for i in range(1, 51):
657  # doing the transformation to probabilities
658  signal[i] = hist_signal.GetBinContent(i)
659  back[i] = hist_background.GetBinContent(i)
660 
661  weight[i] = signal[i] + back[i]
662 
663  # avoid dividing by zero
664  if signal[i] + back[i] == 0:
665  purityB0[i] = 0
666  dilutionB02[i] = 0
667  purityB0bar[i] = 0
668  dilutionB0bar2[i] = 0
669  else:
670  purityB0[i] = signal[i] / (signal[i] + back[i])
671  dilutionB02[i] = -1 + 2 * signal[i] / (signal[i] + back[i])
672 
673  purityB0bar[i] = back[i] / (signal[i] + back[i])
674  dilutionB0bar2[i] = -1 + 2 * back[i] / (signal[i] + back[i])
675 
676  # filling histogram with probability from 0 to 1
677  hist_probB0.Fill(purityB0[i], signal[i])
678  hist_probB0bar.Fill(purityB0bar[i], back[i])
679 
680  # filling histogram with qp from -1 to 1
681  hist_qpB0.Fill(dilutionB02[i], signal[i])
682  hist_qpB0bar.Fill(dilutionB0bar2[i], back[i])
683 
684  # filling histogram with abs(qp), i.e. this histogram contains the r-values (not qp)
685  hist_absqpB0.Fill(abs(dilutionB02[i]), signal[i])
686  hist_absqpB0bar.Fill(abs(dilutionB0bar2[i]), back[i])
687  # filling histogram with abs(qp) special weighted - needed for average r calculation
688  hist_aver_rB0.Fill(abs(dilutionB02[i]), abs(dilutionB02[i]) * signal[i])
689  hist_aver_rB0bar.Fill(abs(dilutionB0bar2[i]), abs(dilutionB0bar2[i]) * back[i])
690 
691  # hist_aver_rB0bar contains now the average r-value
692  hist_aver_rB0.Divide(hist_absqpB0)
693  hist_aver_rB0bar.Divide(hist_absqpB0bar)
694  # now calculating the efficiency
695 
696  # calculating number of events
697  tot_entriesB0 = 0
698  tot_entriesB0bar = 0
699  for i in range(1, r_size):
700  tot_entriesB0 = tot_entriesB0 + hist_absqpB0.GetBinContent(i)
701  tot_entriesB0bar = tot_entriesB0bar + hist_absqpB0bar.GetBinContent(i)
702  # initializing some arrays
703  tot_eff_effB0 = 0
704  tot_eff_effB0bar = 0
705  event_fractionB0 = array('f', [0] * r_size)
706  event_fractionB0bar = array('f', [0] * r_size)
707  rvalueB0 = array('f', [0] * r_size)
708  rvalueB0bar = array('f', [0] * r_size)
709  # wvalue = array('f', [0] * r_size)
710  entriesB0 = array('f', [0] * r_size)
711  entriesB0bar = array('f', [0] * r_size)
712 
713  for i in range(1, r_size):
714  rvalueB0[i] = hist_aver_rB0.GetBinContent(i)
715  rvalueB0bar[i] = hist_aver_rB0bar.GetBinContent(i)
716  # wvalue[i] = (1 - rvalueB0[i]) / 2
717  entriesB0[i] = hist_absqpB0.GetBinContent(i)
718  entriesB0bar[i] = hist_absqpB0bar.GetBinContent(i)
719  event_fractionB0[i] = entriesB0[i] / tot_entriesB0
720  event_fractionB0bar[i] = entriesB0bar[i] / tot_entriesB0bar
721  # print '* Bin ' + str(i) + ' r-value: ' + str(rvalueB0[i]), 'entriesB0: ' +
722  # str(event_fractionB0[i] * 100) + ' % (' + str(entriesB0[i]) + '/' +
723  # str(tot_entriesB0) + ')'
724  tot_eff_effB0 = tot_eff_effB0 + event_fractionB0[i] * rvalueB0[i] \
725  * rvalueB0[i]
726  tot_eff_effB0bar = tot_eff_effB0bar + event_fractionB0bar[i] * rvalueB0bar[i] \
727  * rvalueB0bar[i]
728  effDiff = tot_eff_effB0 - tot_eff_effB0bar
729  effAverage = (tot_eff_effB0 + tot_eff_effB0bar) / 2
730 
731  print(
732  '{:<25}'.format("* " +
733  category) +
734  ' B0-Eff=' +
735  '{: 8.2f}'.format(
736  tot_eff_effB0 *
737  100) +
738  ' %' +
739  ' B0bar-Eff=' +
740  '{: 8.2f}'.format(
741  tot_eff_effB0bar *
742  100) +
743  ' %' +
744  ' EffAverage=' +
745  '{: 8.2f}'.format(effAverage * 100) + ' %' +
746  ' EffDiff=' +
747  '{: 8.2f}'.format(effDiff * 100) + ' % *')
748 
749  # ****** produce the input plots from combiner level ******
750 
751  efficienciesForNtuple.append(float(effAverage * 100))
752 
753  maxSignal = hist_signal.GetBinContent(hist_signal.GetMaximumBin())
754  maxBackground = hist_background.GetBinContent(hist_background.GetMaximumBin())
755 
756  Ymax = max(maxSignal, maxBackground)
757  Ymax = Ymax + Ymax / 12
758 
759  ROOT.gStyle.SetOptStat(0)
760  with Quiet(ROOT.kError):
761  Canvas = ROOT.TCanvas('Bla', 'TITEL BLA', 1200, 800)
762  Canvas.cd() # activate
763  Canvas.SetLogy()
764  Canvas.SetLeftMargin(0.13)
765  Canvas.SetRightMargin(0.04)
766  Canvas.SetTopMargin(0.03)
767  Canvas.SetBottomMargin(0.14)
768  hist_signal.SetFillColorAlpha(ROOT.kBlue, 0.2)
769  hist_signal.SetFillStyle(1001)
770  hist_signal.SetTitleSize(0.1)
771  hist_signal.GetXaxis().SetLabelSize(0.04)
772  hist_signal.GetYaxis().SetLabelSize(0.04)
773  hist_signal.GetXaxis().SetTitleSize(0.06)
774  hist_signal.GetYaxis().SetTitleSize(0.06)
775  hist_signal.GetXaxis().SetLabelSize(0.04)
776  hist_signal.GetYaxis().SetLabelSize(0.04)
777  hist_signal.GetXaxis().SetTitleSize(0.05)
778  hist_signal.GetYaxis().SetTitleSize(0.05)
779  hist_signal.GetXaxis().SetTitleOffset(0.95)
780  hist_signal.GetYaxis().SetTitleOffset(1.1)
781  hist_signal.GetXaxis().SetTitleOffset(1.15)
782  hist_signal.GetYaxis().SetLimits(0, Ymax)
783  hist_signal.SetLineColor(ROOT.kBlue)
784  hist_background.SetFillColorAlpha(ROOT.kRed, 1.0)
785  hist_background.SetFillStyle(3005)
786  hist_background.GetYaxis().SetLimits(0, Ymax)
787  hist_background.SetLineColor(ROOT.kRed)
788 
789  catName = category
790  if category == 'MaximumPstar':
791  catName = 'MaximumP*'
792 
793  hist_signal.SetTitle('; (#it{qp})^{' + catName + '} ; Events')
794  # hist_signal.SetMinimum(0)
795  hist_signal.SetMaximum(Ymax)
796  # hist_background.SetMinimum(0)
797  hist_background.SetMaximum(Ymax)
798 
799  hist_signal.Draw('hist')
800  hist_background.Draw('hist same')
801 
802  l0 = ROOT.TLegend(0.13, 0.65, 0.33, 0.97)
803  l0.SetFillColorAlpha(ROOT.kWhite, 0)
804  l0.AddEntry(hist_signal, ' #varepsilon_{eff}(B^{0}) = ' + '{:.2f}'.format(tot_eff_effB0 * 100) + "%")
805  l0.AddEntry(hist_background, ' #varepsilon_{eff}(#bar{B}^{0}) = ' + '{:.2f}'.format(tot_eff_effB0bar * 100) + "%")
806  l0.AddEntry("", "#bf{#varepsilon_{eff} = " + '{:.2f}'.format(effAverage * 100) + '%}')
807  l0.AddEntry("", '#Delta#varepsilon_{eff} = ' + '{:^5.2f}'.format(effDiff * 100) + '%')
808  l0.SetBorderSize(0)
809  l0.SetTextSize(0.045)
810  l0.Draw()
811 
812  l1 = ROOT.TLegend(0.85, 0.7, 0.96, 0.97)
813  l1.SetFillColorAlpha(ROOT.kWhite, 0.35)
814  l1.AddEntry(hist_signal, 'B^{0}_{MC}')
815  l1.AddEntry(hist_background, '#bar{B}^{0}_{MC}')
816  l1.SetTextSize(0.045)
817  l1.Draw()
818 
819  Canvas.Update()
820  with Quiet(ROOT.kError):
821  Canvas.SaveAs(workingDirectory + '/' + 'test6_CPVFTqp_' + category + '_both.pdf')
822 
823  # Validation Plot 4
824  hist_both.GetXaxis().SetLabelSize(0.04)
825  hist_both.GetYaxis().SetLabelSize(0.04)
826  hist_both.GetYaxis().SetTitleOffset(0.7)
827  hist_both.GetXaxis().SetTitleOffset(0.7)
828  hist_both.GetXaxis().SetTitleSize(0.06)
829  hist_both.GetYaxis().SetTitleSize(0.06)
830 
831  hist_both.GetListOfFunctions().Add(ROOT.TNamed('MetaOptions', 'nostats, logy'))
832  hist_both.GetListOfFunctions().Add(ROOT.TNamed('Description', 'Output of the flavor tagger category ' + catName))
833  hist_both.GetListOfFunctions().Add(
834  ROOT.TNamed('Check', 'Shape should not change drastically. E.g. Warning if there is only a peak at 0.'))
835  hist_both.GetListOfFunctions().Add(ROOT.TNamed('Contact', 'yosato@post.kek.jp'))
836 
837  hist_both.SetTitle(
838  'Flavor tagger output of the category ' +
839  catName +
840  '; #it{qp}_{' +
841  catName +
842  '} ; Events (Total = ' +
843  '{:<1}'.format(
844  "%.0f" %
845  hist_both.GetEntries()) +
846  ')')
847  # hist_both.SetStats(False)
848  hist_both.Write()
849 
850  Canvas.Clear()
851 
852 outputNtuple.Fill(array('f', efficienciesForNtuple))
853 outputNtuple.Write()
854 outputFile.Close()
855 
856 print('* *')
857 print('*******************************************************************************************************************')
oldlevel
the previously set level to be ignored