Belle II Software  light-2212-foldex
flavorTaggerEfficiency.py
1 #!/usr/bin/env python3
2 # -*- coding: utf-8 -*-
3 
4 
11 
12 # ************** Flavor Tagging **************
13 # * *
14 # * This script calculates the parameters char- *
15 # * acterizing the performance of the flavor *
16 # * tagger. It produces plots for the qr *
17 # * distribution of the combiners and the dis- *
18 # * tributions of the combiner input variables. *
19 # * *
20 # ***********************************************
21 
22 import ROOT
23 from basf2 import B2INFO
24 import flavorTagger as ft
25 from defaultEvaluationParameters import categories, Quiet, r_subsample, r_size
26 from array import array
27 import pickle
28 import math
29 import glob
30 import sys
31 
32 ROOT.gROOT.SetBatch(True)
33 
34 if len(sys.argv) != 3:
35  sys.exit("Must provide 2 arguments: [input_sim_file] or ['input_sim_file*'] with wildcards and [treeName]"
36  )
37 workingFile = sys.argv[1]
38 workingFiles = glob.glob(str(workingFile))
39 treeName = str(sys.argv[2])
40 
41 if len(workingFiles) < 1:
42  sys.exit("No file name or file names " + str(workingFile) + " found.")
43 
44 
45 workingDirectory = '.'
46 
47 #
48 # *******************************************
49 # DETERMINATION OF TOTAL EFFECTIVE EFFICIENCY
50 # *******************************************
51 #
52 
53 # working directory
54 # needs the B0bar_final.root-file
55 # treeName = 'B0tree'
56 
57 
58 methods = []
59 
60 tree = ROOT.TChain(treeName)
61 
62 mcstatus = array('d', [-511.5, 0.0, 511.5])
63 ROOT.TH1.SetDefaultSumw2()
64 
65 for iFile in workingFiles:
66  tree.AddFile(iFile)
67 
68 totalBranches = []
69 for branch in tree.GetListOfBranches():
70  totalBranches.append(branch.GetName())
71 
72 if 'FBDT_qrCombined' in totalBranches:
73  methods.append("FBDT")
74 
75 if 'FANN_qrCombined' in totalBranches:
76  methods.append("FANN")
77 
78 if 'DNN_qrCombined' in totalBranches:
79  methods.append("DNN")
80 
81 usedCategories = []
82 for cat in categories:
83  catBranch = 'qp' + cat
84  if catBranch in totalBranches:
85  usedCategories.append(cat)
86 
87 YmaxForQrPlot = 0
88 
89 total_notTagged = 0
90 
91 for method in methods:
92  # Get error with GetBinError(), set error with SetBinError()
93  # histogram contains the average r in each of 7 bins -> calculation see below
94  histo_avr_r = ROOT.TH1F('Average_r', 'Average r in each of the bins (B0 and B0bar)', int(r_size - 2),
95  r_subsample)
96  histo_avr_rB0 = ROOT.TH1F('Average_rB0', 'Average r in each of the bins (B0)', int(r_size - 2),
97  r_subsample)
98  histo_avr_rB0bar = ROOT.TH1F('Average_rB0bar', 'Average r in each of the bins (B0bar)', int(r_size - 2),
99  r_subsample)
100 
101  # histogram containing the wrong tags calculated by counting Wrong tags
102  histo_mc_NwB0 = ROOT.TH1F('mc_NwB0', 'Average r in each of the bins (B0)', int(r_size - 2),
103  r_subsample)
104  histo_mc_NwB0bar = ROOT.TH1F('mc_NwB0bar', 'Average r in each of the bins (B0bar)', int(r_size - 2),
105  r_subsample)
106 
107  # histogram contains the mean squared of r in each of 7 bins -> calculation see below
108  histo_ms_r = ROOT.TH1F('MS_r', 'Mean squared average of r in each of the bins (B0 and B0bar)', int(r_size - 2),
109  r_subsample)
110  histo_ms_rB0 = ROOT.TH1F('MS_rB0', 'Mean squared average of r in each of the bins (B0)', int(r_size - 2),
111  r_subsample)
112  histo_ms_rB0bar = ROOT.TH1F('MS_rB0bar', 'Mean squared average of r in each of the bins (B0bar)', int(r_size - 2),
113  r_subsample)
114 
115  # histogram with number of entries in for each bin
116  histo_entries_per_bin = ROOT.TH1F(
117  'entries_per_bin',
118  'Events binned in r_subsample according to their r-value for B0 and B0bar prob',
119  int(r_size - 2),
120  r_subsample)
121  histo_entries_per_binB0 = ROOT.TH1F('entries_per_binB0', 'Events binned in r_subsample according '
122  'to their r-value for B0 prob', int(r_size - 2), r_subsample)
123  histo_entries_per_binB0bar = ROOT.TH1F('entries_per_binB0bar',
124  'Events binned in r_subsample according to their r-value '
125  'for B0bar prob', int(r_size - 2), r_subsample)
126  # histogram network output (not qr and not r) for true B0 (signal) - not necessary
127  histo_Cnet_output_B0 = ROOT.TH1F('Comb_Net_Output_B0', 'Combiner network output [not equal to r] '
128  'for true B0 (binning 100)', 100, 0.0, 1.0)
129  # histogram network output (not qr and not r) for true B0bar (background) - not necessary
130  histo_Cnet_output_B0bar = ROOT.TH1F('Comb_Net_Output_B0bar', 'Combiner network output [not equal to r] '
131  'for true B0bar (binning 100)', 100, 0.0, 1.0)
132  # histogram containing the belle paper plot (qr-tagger output for true B0)
133  histo_belleplotB0 = ROOT.TH1F('BellePlot_B0',
134  'BellePlot for true B0 (binning 50)', 50,
135  -1.0, 1.0)
136  # histogram containing the belle paper plot (qr-tagger output for true B0bar)
137  histo_belleplotB0bar = ROOT.TH1F('BellePlot_B0Bar',
138  'BellePlot for true B0Bar (binning 50)',
139  50, -1.0, 1.0)
140 
141  histo_notTaggedEvents = ROOT.TH1F('notTaggedEvents',
142  'Histogram for not tagged events',
143  1, -3.0, -1.0)
144 
145  # calibration plot for B0. If we get a linaer line our MC is fine, than the assumption r ~ 1- 2w is reasonable
146  # expectation is, that for B0 calibration plot: qr=0 half B0 and half B0bar, qr = 1 only B0 and qr = -1
147  # no B0. Inverse for B0bar calibration plot
148  histo_calib_B0 = ROOT.TH1F('Calibration_B0', 'CalibrationPlot for true B0', 100, -1.0, 1.0)
149  # calibration plot for B0bar calibration plot
150  histo_calib_B0bar = ROOT.TH1F('Calibration_B0Bar',
151  'CalibrationPlot for true B0Bar', 100, -1.0,
152  1.0)
153  # belle plot with true B0 and B0bars
154  hallo12 = ROOT.TH1F('BellePlot_NoCut', 'BellePlot_NoCut (binning 100)',
155  100, -1.0, 1.0)
156 
157 
158  diag = ROOT.TF1('diag', 'pol1', -1, 1)
159 
160 
161  # histograms for the efficiency calculation in wrong way
162  histo_m0 = ROOT.TH1F('BellePlot_m0',
163  'BellePlot_m for true B0 (binning 50)', 50, -1.0, 1.0)
164  histo_m1 = ROOT.TH1F('BellePlot_m1',
165  'BellePlot_m for true B0 (binning 50)', 50, -1.0, 1.0)
166  histo_m2 = ROOT.TH1F('BellePlot_m2',
167  'BellePlot_m for true B0Bar (binning 50)', 50, -1.0,
168  1.0)
169 
170  # filling the histograms
171 
172  tree.Draw(method + '_qrCombined>>BellePlot_B0', 'qrMC == 1')
173  tree.Draw(method + '_qrCombined>>BellePlot_B0Bar', 'qrMC == -1')
174  tree.Draw(method + '_qrCombined>>BellePlot_NoCut', 'abs(qrMC) == 1')
175 
176  tree.Draw(method + '_qrCombined>>Calibration_B0', 'qrMC == 1')
177  tree.Draw(method + '_qrCombined>>Calibration_B0Bar', 'qrMC == -1')
178 
179  tree.Draw(method + '_qrCombined>>notTaggedEvents',
180  'abs(qrMC) == 0 && isSignal == 1 && ' +
181  method + '_qrCombined < -1')
182 
183  # filling histograms wrong efficiency calculation
184  tree.Draw(method + '_qrCombined>>BellePlot_m0',
185  'qrMC == 1 && ' + method + '_qrCombined>0')
186  tree.Draw(method + '_qrCombined>>BellePlot_m1',
187  'qrMC == 1 && ' + method + '_qrCombined<0')
188  tree.Draw(method + '_qrCombined>>BellePlot_m2',
189  'qrMC == -1 && ' + method + '_qrCombined>0 ')
190 
191  # filling with abs(qr) in one of 6 bins with its weight
192  # separate calculation for B0 and B0bar
193 
194  tree.Project('Average_r', 'abs(' + method + '_qrCombined)', 'abs(' + method + '_qrCombined)*(abs(qrMC) == 1)')
195  tree.Project('Average_rB0', 'abs(' + method + '_qrCombined)', 'abs(' + method + '_qrCombined)*(qrMC==1)')
196  tree.Project('Average_rB0bar', 'abs(' + method + '_qrCombined)', 'abs(' + method + '_qrCombined)*(qrMC==-1)')
197 
198  tree.Project('MS_r', 'abs(' + method + '_qrCombined)', '(' + method +
199  '_qrCombined*' + method + '_qrCombined)*(abs(qrMC) == 1)')
200  tree.Project('MS_rB0', 'abs(' + method + '_qrCombined)',
201  '(' + method + '_qrCombined*' + method + '_qrCombined)*(qrMC==1)')
202  tree.Project('MS_rB0bar', 'abs(' + method + '_qrCombined)',
203  '(' + method + '_qrCombined*' + method + '_qrCombined)*(qrMC==-1)')
204 
205  # filling with abs(qr) in one of 6 bins
206  tree.Project('entries_per_bin', 'abs(' + method + '_qrCombined)', 'abs(qrMC) == 1')
207  tree.Project('entries_per_binB0', 'abs(' + method + '_qrCombined)', 'qrMC == 1')
208  tree.Project('entries_per_binB0bar', 'abs(' + method + '_qrCombined)', 'qrMC == -1')
209 
210  # filling histograms with number of wrong tags per r-bin
211  tree.Project('mc_NwB0', 'abs(' + method + '_qrCombined)', ' ' + method + '_qrCombined*qrMC < 0 && qrMC == 1')
212  tree.Project('mc_NwB0bar', 'abs(' + method + '_qrCombined)', ' ' + method + '_qrCombined*qrMC < 0 && qrMC == -1')
213 
214  # producing the average r histograms
215  histo_avr_r.Divide(histo_entries_per_bin)
216  histo_avr_rB0.Divide(histo_entries_per_binB0)
217  histo_avr_rB0bar.Divide(histo_entries_per_binB0bar)
218 
219  histo_ms_r.Divide(histo_entries_per_bin)
220  histo_ms_rB0.Divide(histo_entries_per_binB0)
221  histo_ms_rB0bar.Divide(histo_entries_per_binB0bar)
222 
223  # producing the calibration plots
224  # Errors ok
225  histo_calib_B0.Divide(hallo12)
226  histo_calib_B0bar.Divide(hallo12)
227 
228  # Fit for calibration plot
229  print(' ')
230  print('****************************** CALIBRATION CHECK FOR COMBINER USING ' +
231  method + ' ***************************************')
232  print(' ')
233  print('Fit ploynomial of first order to the calibration plot. Expected value ~0.5')
234  print(' ')
235  histo_calib_B0.Fit(diag, 'TEST')
236  print(' ')
237  print('****************************** MEASURED EFFECTIVE EFFICIENCY FOR COMBINER USING ' +
238  method + ' ***********************************')
239  print('* ' +
240  ' *')
241  # get total number of entries
242  total_tagged = histo_entries_per_bin.GetEntries()
243  total_tagged_B0 = histo_entries_per_binB0.GetEntries()
244  total_tagged_B0bar = histo_entries_per_binB0bar.GetEntries()
245  total_notTagged = histo_notTaggedEvents.GetEntries()
246  total_entries = (total_tagged + total_notTagged)
247  # To a good approximation we assume that half of the not tagged B mesons were B0 (B0bar)
248  total_entriesB0 = (total_tagged_B0 + total_notTagged / 2)
249  total_entriesB0bar = (total_tagged_B0bar + total_notTagged / 2)
250 
251  tagging_eff = total_tagged / (total_tagged + total_notTagged)
252  DeltaTagging_eff = math.sqrt(total_tagged * total_notTagged**2 + total_notTagged * total_tagged**2) / (total_entries**2)
253  tot_eff_effB0 = 0
254  tot_eff_effB0bar = 0
255  average_eff_eff = 0
256  uncertainty_eff_effB0 = 0
257  uncertainty_eff_effB0bar = 0
258  uncertainty_eff_effAverage = 0
259  diff_eff_Uncertainty = 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_fractionTotalUncertainty = array('f', [0] * r_size)
264  eventsInBin_B0 = array('f', [0] * r_size)
265  eventsInBin_B0bar = array('f', [0] * r_size)
266  eventsInBin_Total = array('f', [0] * r_size)
267  event_fractionDiff = array('f', [0] * r_size)
268  event_fractionDiffUncertainty = 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  rvalueStdB0 = array('f', [0] * r_size)
273  rvalueStdB0bar = array('f', [0] * r_size)
274  rvalueStdB0Average = array('f', [0] * r_size)
275  wvalue = array('f', [0] * r_size)
276  wvalueUncertainty = array('f', [0] * r_size)
277  wvalueB0 = array('f', [0] * r_size)
278  wvalueB0bar = array('f', [0] * r_size)
279  wvalueB0Uncertainty = array('f', [0] * r_size)
280  wvalueB0barUncertainty = array('f', [0] * r_size)
281  wvalueDiff = array('f', [0] * r_size)
282  wvalueDiffUncertainty = array('f', [0] * r_size)
283  entries = array('f', [0] * r_size)
284  entriesB0 = array('f', [0] * r_size)
285  entriesB0bar = array('f', [0] * r_size)
286  iEffEfficiency = array('f', [0] * r_size)
287  iEffEfficiencyUncertainty = array('f', [0] * r_size)
288  iEffEfficiencyB0Uncertainty = array('f', [0] * r_size)
289  iEffEfficiencyB0barUncertainty = array('f', [0] * r_size)
290  iEffEfficiencyB0UncertaintyFromOutput = array('f', [0] * r_size)
291  iEffEfficiencyB0barUncertaintyFromOutput = array('f', [0] * r_size)
292 
293  iDeltaEffEfficiency = array('f', [0] * r_size)
294  iDeltaEffEfficiencyUncertainty = array('f', [0] * r_size)
295  muParam = array('f', [0] * r_size)
296  muParamUncertainty = array('f', [0] * r_size)
297  # intervallEff = array('f', [0] * r_size)
298 
299  print('* --> DETERMINATION BASED ON MONTE CARLO ' +
300  ' *')
301  print('* ' +
302  ' *')
303  print('* Note: mu = Delta_Effcy/(2*Efficiency). Needed for CP analysis ' +
304  'together with w and Delta_w *')
305  print('* ' +
306  ' *')
307  print('* ------------------------------------------------------------------' +
308  '-------------------------------------------------- *')
309  print('* r-interval <r> Efficiency Delta_Effcy ' +
310  ' mu w Delta_w *')
311  print('* ------------------------------------------------------------------' +
312  '-------------------------------------------------- *')
313  performance = []
314  for i in range(1, r_size):
315  # get the average r-value
316  entries[i] = histo_entries_per_bin.GetBinContent(i)
317  entriesB0[i] = histo_entries_per_binB0.GetBinContent(i)
318  entriesB0bar[i] = histo_entries_per_binB0bar.GetBinContent(i)
319  # fraction of events/all events
320 
321  event_fractionB0[i] = entriesB0[i] / total_entriesB0
322  event_fractionB0bar[i] = entriesB0bar[i] / total_entriesB0bar
323 
324  # event_fractionTotal[i] = (entriesB0[i] + entriesB0bar[i]) / total_entries
325  # event_fractionDiff[i] = (entriesB0[i] - entriesB0bar[i]) / total_entries
326 
327  event_fractionTotal[i] = (event_fractionB0[i] + event_fractionB0bar[i]) / 2
328  event_fractionDiff[i] = event_fractionB0[i] - event_fractionB0bar[i]
329 
330  event_fractionDiffUncertainty[i] = math.sqrt(entriesB0[i] *
331  (total_entriesB0 -
332  entriesB0[i]) /
333  total_entriesB0**3 +
334  entriesB0bar[i] *
335  (total_entriesB0bar -
336  entriesB0bar[i]) /
337  total_entriesB0bar**3)
338 
339  event_fractionTotalUncertainty[i] = event_fractionDiffUncertainty[i] / 2
340 
341  rvalueB0[i] = histo_avr_rB0.GetBinContent(i)
342  rvalueB0bar[i] = histo_avr_rB0bar.GetBinContent(i)
343  rvalueB0Average[i] = histo_avr_r.GetBinContent(i) # (rvalueB0[i] + rvalueB0bar[i]) / 2
344  rvalueStdB0[i] = math.sqrt(histo_ms_rB0.GetBinContent(
345  i) - (histo_avr_rB0.GetBinContent(i))**2) / math.sqrt(entriesB0[i] - 1)
346  rvalueStdB0bar[i] = math.sqrt(histo_ms_rB0bar.GetBinContent(
347  i) - (histo_avr_rB0bar.GetBinContent(i))**2) / math.sqrt(entriesB0bar[i] - 1)
348  rvalueStdB0Average[i] = math.sqrt(rvalueStdB0[i]**2 + rvalueStdB0bar[i]**2) / 2
349  # math.sqrt(histo_ms_r.GetBinContent(i) - (histo_avr_r.GetBinContent(i))**2)
350  # calculate the wrong tag fractin (only true if MC data good)
351 
352  wvalueB0[i] = histo_mc_NwB0.GetBinContent(i) / entriesB0[i] # (1 - rvalueB0[i]) / 2
353  wvalueB0bar[i] = histo_mc_NwB0bar.GetBinContent(i) / entriesB0bar[i] # (1 - rvalueB0bar[i]) / 2
354  wvalueDiff[i] = wvalueB0[i] - wvalueB0bar[i]
355  wvalueB0Uncertainty[i] = math.sqrt(histo_mc_NwB0.GetBinContent(
356  i) * (entriesB0[i] - histo_mc_NwB0.GetBinContent(i)) / (entriesB0[i]**3))
357  wvalueB0barUncertainty[i] = math.sqrt(histo_mc_NwB0bar.GetBinContent(
358  i) * (entriesB0bar[i] - histo_mc_NwB0bar.GetBinContent(i)) / (entriesB0bar[i]**3))
359  # math.sqrt((rvalueStdB0[i] / 2)**2 + (rvalueStdB0bar[i] / 2)**2)
360  wvalueDiffUncertainty[i] = math.sqrt(wvalueB0Uncertainty[i]**2 + wvalueB0barUncertainty[i]**2)
361  wvalue[i] = (wvalueB0[i] + wvalueB0bar[i]) / 2
362  wvalueUncertainty[i] = wvalueDiffUncertainty[i] / 2
363 
364  # Avr Efficiency
365 
366  # iEffEfficiency[i] = (event_fractionB0[i] * rvalueB0[i] * rvalueB0[i] +
367  # event_fractionB0bar[i] * rvalueB0bar[i] * rvalueB0bar[i]) / 2
368 
369  iEffEfficiency[i] = event_fractionTotal[i] * (1 - 2 * wvalue[i])**2
370 
371  iEffEfficiencyUncertainty[i] = (1 - 2 * wvalue[i]) * \
372  math.sqrt((2 * event_fractionTotal[i] * 2 * wvalueUncertainty[i])**2 +
373  (1 - 2 * wvalue[i])**2 * event_fractionTotalUncertainty[i]**2)
374 
375  # iEffEfficiencyUncertainty[i] = rvalueB0Average[i] * \
376  # math.sqrt((2 * total_entries * entries[i] * rvalueStdB0Average[i])**2 +
377  # rvalueB0Average[i]**2 * entries[i] *
378  # (total_entries * (total_entries - entries[i]) +
379  # entries[i] * total_notTagged)) / (total_entries**2)
380 
381  # iEffEfficiencyB0UncertaintyFromOutput[i] = rvalueB0[i] * \
382  # math.sqrt((2 * total_entriesB0 * entriesB0[i] * rvalueStdB0[i])**2 +
383  # rvalueB0[i]**2 * entriesB0[i] *
384  # total_entriesB0 * (total_entriesB0 - entriesB0[i])) / (total_entriesB0**2)
385  # iEffEfficiencyB0barUncertaintyFromOutput[i] = rvalueB0bar[i] * \
386  # math.sqrt((2 * total_entriesB0bar * entriesB0bar[i] * rvalueStdB0bar[i])**2 +
387  # rvalueB0bar[i]**2 * entriesB0bar[i] *
388  # total_entriesB0bar * (total_entriesB0bar - entriesB0bar[i])) / (total_entriesB0bar**2)
389 
390  # iEffEfficiencyUncertainty[i] =
391  # math.sqrt(iEffEfficiencyB0UncertaintyFromOutput[i]**2 +
392  # iEffEfficiencyB0barUncertaintyFromOutput[i]**2) / 2
393 
394  # iEffEfficiency[i] = (event_fractionB0[i] * (1-2*wvalueB0[i])**2 + event_fractionB0bar[i] * (1-2*wvalueB0bar[i])**2)/2
395 
396  average_eff_eff += iEffEfficiency[i]
397  # average_eff_eff += (event_fractionB0[i] * (1-2*wvalueB0[i])**2 + event_fractionB0bar[i] * (1-2*wvalueB0bar[i])**2)/2
398 
399  # Delta Eff
400 
401  iDeltaEffEfficiency[i] = event_fractionB0[i] * (1 - 2 * wvalueB0[i])**2 - \
402  event_fractionB0bar[i] * (1 - 2 * wvalueB0bar[i])**2
403 
404  iEffEfficiencyB0Uncertainty[i] = (1 - 2 * wvalueB0[i]) * \
405  math.sqrt((2 * total_entriesB0 * entriesB0[i] * 2 * wvalueB0Uncertainty[i])**2 +
406  (1 - 2 * wvalueB0[i])**2 * entriesB0[i] *
407  total_entriesB0 * (total_entriesB0 - entriesB0[i])) / (total_entriesB0**2)
408  iEffEfficiencyB0barUncertainty[i] = (1 - 2 * wvalueB0bar[i]) * \
409  math.sqrt((2 * total_entriesB0bar * entriesB0bar[i] * 2 * wvalueB0barUncertainty[i])**2 +
410  (1 - 2 * wvalueB0bar[i])**2 * entriesB0bar[i] *
411  total_entriesB0bar * (total_entriesB0bar - entriesB0bar[i])) / (total_entriesB0bar**2)
412 
413  iDeltaEffEfficiencyUncertainty[i] = math.sqrt(iEffEfficiencyB0Uncertainty[i]**2 + iEffEfficiencyB0barUncertainty[i]**2)
414  # iEffEfficiencyUncertainty[i] = iDeltaEffEfficiencyUncertainty[i]/2
415 
416  diff_eff_Uncertainty = diff_eff_Uncertainty + iDeltaEffEfficiencyUncertainty[i]**2
417 
418  # finally calculating the total effective efficiency
419  tot_eff_effB0 = tot_eff_effB0 + event_fractionB0[i] * (1 - 2 * wvalueB0[i])**2
420  tot_eff_effB0bar = tot_eff_effB0bar + event_fractionB0bar[i] * (1 - 2 * wvalueB0bar[i])**2
421  uncertainty_eff_effAverage = uncertainty_eff_effAverage + iEffEfficiencyUncertainty[i]**2
422  uncertainty_eff_effB0 = uncertainty_eff_effB0 + iEffEfficiencyB0Uncertainty[i]**2
423  uncertainty_eff_effB0bar = uncertainty_eff_effB0bar + iEffEfficiencyB0barUncertainty[i]**2
424  muParam[i] = event_fractionDiff[i] / (2 * event_fractionTotal[i])
425  muParamUncertainty[i] = event_fractionDiffUncertainty[i] / (2 * event_fractionTotal[i]) * math.sqrt(muParam[i]**2 + 1)
426 
427  # intervallEff[i] = event_fractionTotal[i] * rvalueB0Average[i] * rvalueB0Average[i]
428  print('* ' + '{:.3f}'.format(r_subsample[i - 1]) + ' - ' + '{:.3f}'.format(r_subsample[i]) + ' ' +
429  '{:.3f}'.format(rvalueB0Average[i]) + ' +- ' + '{:.4f}'.format(rvalueStdB0Average[i]) + ' ' +
430  '{:.4f}'.format(event_fractionTotal[i]) + ' ' +
431  '{: .4f}'.format(event_fractionDiff[i]) + ' +- ' + '{:.4f}'.format(event_fractionDiffUncertainty[i]) + ' ' +
432  '{: .4f}'.format(muParam[i]) + ' +- ' + '{:.4f}'.format(muParamUncertainty[i]) + ' ' +
433  '{:.4f}'.format(wvalue[i]) + ' +- ' + '{:.4f}'.format(wvalueUncertainty[i]) + ' ' +
434  '{: .4f}'.format(wvalueDiff[i]) + ' +- ' + '{:.4f}'.format(wvalueDiffUncertainty[i]) + ' *')
435 
436  # average_eff_eff = (tot_eff_effB0 + tot_eff_effB0bar) / 2
437  uncertainty_eff_effAverage = math.sqrt(uncertainty_eff_effAverage)
438  uncertainty_eff_effB0 = math.sqrt(uncertainty_eff_effB0)
439  uncertainty_eff_effB0bar = math.sqrt(uncertainty_eff_effB0bar)
440  diff_eff = tot_eff_effB0 - tot_eff_effB0bar
441  diff_eff_Uncertainty = math.sqrt(diff_eff_Uncertainty)
442  print('* --------------------------------------------------------------------------------------------------' +
443  '------------------ *')
444  print('* *')
445  print('* __________________________________________________________________________________________ *')
446  print('* | | *')
447  print('* | TOTAL NUMBER OF TAGGED EVENTS = ' +
448  '{:<24}'.format("%.0f" % total_tagged) + '{:>36}'.format('| *'))
449  print('* | | *')
450  print(
451  '* | TOTAL AVERAGE EFFICIENCY (q=+-1)= ' +
452  '{:.2f}'.format(
453  tagging_eff *
454  100) +
455  " +- " +
456  '{:.2f}'.format(
457  DeltaTagging_eff *
458  100) +
459  ' % | *')
460  print('* | | *')
461  print(
462  '* | TOTAL AVERAGE EFFECTIVE EFFICIENCY (q=+-1)= ' +
463  '{:.6f}'.format(
464  average_eff_eff *
465  100) +
466  " +- " +
467  '{:.6f}'.format(
468  uncertainty_eff_effAverage *
469  100) +
470  ' % | *')
471  print('* | | *')
472  print(
473  '* | TOTAL AVERAGE EFFECTIVE EFFICIENCY ASYMMETRY (q=+-1)= ' +
474  '{:^9.6f}'.format(
475  diff_eff *
476  100) +
477  " +- " +
478  '{:.6f}'.format(
479  diff_eff_Uncertainty *
480  100) +
481  ' % | *')
482  print('* | | *')
483  print('* | B0-TAGGER TOTAL EFFECTIVE EFFICIENCIES: ' +
484  '{:.2f}'.format(tot_eff_effB0 * 100) + " +-" + '{: 4.2f}'.format(uncertainty_eff_effB0 * 100) +
485  ' % (q=+1) ' +
486  '{:.2f}'.format(tot_eff_effB0bar * 100) + " +-" + '{: 4.2f}'.format(uncertainty_eff_effB0bar * 100) +
487  ' % (q=-1) ' + ' | *')
488  print('* | | *')
489  print('* | FLAVOR PERCENTAGE (MC): ' +
490  '{:.2f}'.format(total_tagged_B0 / total_tagged * 100) + ' % (q=+1) ' +
491  '{:.2f}'.format(total_tagged_B0bar / total_tagged * 100) + ' % (q=-1) Diff=' +
492  '{:^5.2f}'.format((total_tagged_B0 - total_tagged_B0bar) / total_tagged * 100) + ' % | *')
493  print('* |__________________________________________________________________________________________| *')
494  print('* *')
495  print('****************************************************************************************************')
496  print('* *')
497 
498  # not that imortant
499  print('* --------------------------------- *')
500  print('* Efficiency Determination - easiest way *')
501  print('* --------------------------------- *')
502  total_tagged_B0 = histo_belleplotB0.GetEntries()
503  total_tagged_B0Bar = histo_belleplotB0bar.GetEntries()
504  total_tagged_wrong = histo_m1.GetEntries()
505  total_tagged_B0Bar_wrong = histo_m2.GetEntries()
506  total_tagged = total_tagged_B0 + total_tagged_B0Bar
507  total_tagged_wrong = total_tagged_wrong + total_tagged_B0Bar_wrong
508 
509  wrong_tag_fraction_B0 = total_tagged_wrong / total_tagged_B0
510  wrong_tag_fraction_B0Bar = total_tagged_B0Bar_wrong / total_tagged_B0Bar
511  wrong_tag_fraction = total_tagged_wrong / total_tagged
512  right_tag_fraction_B0 = 1 - 2 * wrong_tag_fraction_B0
513  right_tag_fraction_B0Bar = 1 - 2 * wrong_tag_fraction_B0Bar
514  right_tag_fraction = 1 - 2 * wrong_tag_fraction
515  wrong_eff_B0 = right_tag_fraction_B0 * right_tag_fraction_B0
516  wrong_eff_B0Bar = right_tag_fraction_B0Bar * right_tag_fraction_B0Bar
517  wrong_eff = right_tag_fraction * right_tag_fraction
518 
519  print('* wrong_tag_fraction for all: ' +
520  '{:.3f}'.format(wrong_tag_fraction * 100) +
521  ' % *')
522  print('* right_tag_fraction for all: ' +
523  '{:.3f}'.format(right_tag_fraction * 100) +
524  ' % *')
525  print('* wrong calculated eff all: ' + '{:.3f}'.format(wrong_eff * 100) +
526  ' % *')
527  print('* *')
528  print('****************************************************************************************************')
529  print('')
530  print('Table For B2TIP')
531  print('')
532  # write out the histograms
533  # histo_avr_r.Write('', ROOT.TObject.kOverwrite)
534  # histo_entries_per_bin.Write('', ROOT.TObject.kOverwrite)
535 
536  # histo_Cnet_output_B0.Write('', ROOT.TObject.kOverwrite)
537  # histo_Cnet_output_B0bar.Write('', ROOT.TObject.kOverwrite)
538  # histo_belleplotB0.Write('', ROOT.TObject.kOverwrite)
539  # histo_belleplotB0bar.Write('', ROOT.TObject.kOverwrite)
540  # histo_calib_B0.Write('', ROOT.TObject.kOverwrite)
541  # histo_calib_B0bar.Write('', ROOT.TObject.kOverwrite)
542 
543  maxB0 = histo_belleplotB0.GetBinContent(histo_belleplotB0.GetMaximumBin())
544  maxB0bar = histo_belleplotB0bar.GetBinContent(histo_belleplotB0bar.GetMaximumBin())
545 
546  Ymax = max(maxB0, maxB0bar)
547  Ymax = Ymax + Ymax / 12
548 
549  if YmaxForQrPlot < Ymax:
550  YmaxForQrPlot = Ymax
551 
552  # produce a pdf
553  ROOT.gStyle.SetOptStat(0)
554  Canvas1 = ROOT.TCanvas('Bla' + method, 'Final Output', 1200, 800)
555  Canvas1.cd() # activate
556  histo_belleplotB0.SetFillColorAlpha(ROOT.kBlue, 0.2)
557  histo_belleplotB0.SetFillStyle(1001)
558  histo_belleplotB0.GetYaxis().SetLabelSize(0.03)
559  histo_belleplotB0.GetYaxis().SetLimits(0, YmaxForQrPlot)
560  histo_belleplotB0.GetYaxis().SetTitleOffset(1.2)
561  histo_belleplotB0.SetLineColor(ROOT.kBlue)
562  histo_belleplotB0bar.SetFillColorAlpha(ROOT.kRed, 1.0)
563  histo_belleplotB0bar.SetFillStyle(3005)
564  histo_belleplotB0bar.SetLineColor(ROOT.kRed)
565  # SetLabelSize etc SetTitle
566 
567  histo_belleplotB0.SetTitle('Final Flavor Tagger Output; #it{qr}-output ; Events'
568  )
569  histo_belleplotB0.SetMinimum(0)
570  histo_belleplotB0.SetMaximum(YmaxForQrPlot)
571  histo_belleplotB0.Draw('hist')
572  histo_belleplotB0bar.Draw('hist same')
573 
574  leg = ROOT.TLegend(0.75, 0.8, 0.9, 0.9)
575  leg.AddEntry(histo_belleplotB0, 'true B0')
576  leg.AddEntry(histo_belleplotB0bar, 'true B0bar')
577  leg.Draw()
578 
579  Canvas1.Update()
580  # IPython.embed()
581  with Quiet(ROOT.kError):
582  Canvas1.SaveAs(workingDirectory + '/' + 'PIC_Belleplot_both' + method + '.pdf')
583 
584  # produce the nice calibration plot
585  Canvas2 = ROOT.TCanvas('Bla2' + method, 'Calibration plot for true B0', 1200, 800)
586  Canvas2.cd() # activate
587  histo_calib_B0.SetFillColorAlpha(ROOT.kBlue, 0.2)
588  histo_calib_B0.SetFillStyle(1001)
589  histo_calib_B0.GetYaxis().SetTitleOffset(1.2)
590  histo_calib_B0.SetLineColor(ROOT.kBlue)
591 
592  histo_calib_B0.SetTitle('Calibration For True B0; #it{qr}-output ; Calibration '
593  )
594  histo_calib_B0.Draw('hist')
595  diag.Draw('SAME')
596  Canvas2.Update()
597  with Quiet(ROOT.kError):
598  Canvas2.SaveAs(workingDirectory + '/' + 'PIC_Calibration_B0' + method + '.pdf')
599 
600  histo_avr_r.Delete()
601  histo_avr_rB0.Delete()
602  histo_avr_rB0bar.Delete()
603  histo_ms_r.Delete()
604  histo_ms_rB0.Delete()
605  histo_ms_rB0bar.Delete()
606  histo_mc_NwB0.Delete()
607  histo_mc_NwB0bar.Delete()
608  histo_notTaggedEvents.Delete()
609  histo_entries_per_bin.Delete()
610  histo_entries_per_binB0.Delete()
611  histo_entries_per_binB0bar.Delete()
612  histo_Cnet_output_B0.Delete()
613  histo_Cnet_output_B0bar.Delete()
614  histo_belleplotB0.Delete()
615  histo_belleplotB0bar.Delete()
616  histo_calib_B0.Delete()
617  histo_calib_B0bar.Delete()
618  hallo12.Delete()
619  histo_m0.Delete()
620  histo_m1.Delete()
621  histo_m2.Delete()
622  Canvas1.Clear()
623  Canvas2.Clear()
624 
625  print(r'\begin{tabularx}{1\textwidth}{@{}r r r r r r r@{}}')
626  print(r'\hline')
627  print(r'$r$- Interval $\enskip$ & $\varepsilon_i\ $ & $\Delta\varepsilon_i\ $ & $w_i \pm \delta w_i\enskip\; $ ' +
628  r' & $\Delta w_i \pm \delta\Delta w_i $& $\varepsilon_{\text{eff}, i} \pm \delta\varepsilon_{\text{eff}, i}\enskip\, $ ' +
629  r' & $\Delta \varepsilon_{\text{eff}, i} \pm \delta\Delta \varepsilon_{\text{eff}, i}\enskip\, $\\ \hline\hline')
630  for i in range(1, r_size):
631  print('$ ' + '{:.3f}'.format(r_subsample[i - 1]) + ' - ' + '{:.3f}'.format(r_subsample[i]) + '$ & $'
632  '{: 6.1f}'.format(event_fractionTotal[i] * 100) + r'$ & $' +
633  '{: 7.2f}'.format(event_fractionDiff[i] * 100) + r'\;$ & $' +
634  '{: 7.2f}'.format(wvalue[i] * 100) + r" \pm " + '{:2.2f}'.format(wvalueUncertainty[i] * 100) + r'\enskip $ & $' +
635  '{: 7.2f}'.format(wvalueDiff[i] * 100) + r" \pm " +
636  '{:2.2f}'.format(wvalueDiffUncertainty[i] * 100) + r'\enskip $ & $' +
637  '{: 8.4f}'.format(iEffEfficiency[i] * 100) + # + '$ & $' +
638  r" \pm " + '{:2.4f}'.format(iEffEfficiencyUncertainty[i] * 100) + r'\, $ & $' +
639  '{: 6.4f}'.format(iDeltaEffEfficiency[i] * 100) + # +
640  r" \pm " + '{:2.4f}'.format(iDeltaEffEfficiencyUncertainty[i] * 100) +
641  r'\enskip\enskip $ \\ ')
642  print(r'\hline\hline')
643  print(r'\multicolumn{1}{r}{Total} & & \multicolumn{5}{r}{ $\varepsilon_\text{eff} = ' +
644  r'\sum_i \varepsilon_i \cdot \langle 1-2w_i\rangle^2 = ' +
645  '{: 6.2f}'.format(average_eff_eff * 100) + r" \pm " + '{: 6.2f}'.format(uncertainty_eff_effAverage * 100) + r'\enskip\, ')
646  print(r'\Delta \varepsilon_\text{eff} = ' +
647  '{: 6.2f}'.format(diff_eff * 100) + r" \pm " + '{: 6.2f}'.format(diff_eff_Uncertainty * 100) + r'\quad\ $ }' +
648  r' \\ ')
649  print(r'\hline')
650  print(r'\end{tabular}')
651  print('')
652 
653  print('')
654  print('Mu-Values for Table')
655  print('')
656 
657  for i in range(1, r_size):
658  print('$ ' + '{:.3f}'.format(r_subsample[i - 1]) + ' - ' + '{:.3f}'.format(r_subsample[i]) + '$ & $'
659  '{: 7.2f}'.format(muParam[i] * 100) + r" \pm " + '{:2.2f}'.format(muParamUncertainty[i] * 100) + r' $ & ')
660 # ************************************************
661 # DETERMINATION OF INDIVIDUAL EFFECTIVE EFFICIENCY
662 # ************************************************
663 
664 # keep in mind:
665 # the individual efficiency is determined on basis of the combiner training.
666 # Whereas the efficiency is determined on basis of the final expert output.
667 
668 print(ft.getEventLevelParticleLists(usedCategories))
669 
670 
671 print('******************************************* MEASURED EFFECTIVE EFFICIENCY FOR INDIVIDUAL CATEGORIES ' +
672  '**********************************************')
673 print('* ' +
674  ' *')
675 # input: Classifier input from event-level. Output of event-level is recalculated for input on combiner-level.
676 # but is re-evaluated under combiner target. Signal is B0, background is B0Bar.
677 categoriesPerformance = []
678 NbinsCategories = 100
679 for category in usedCategories:
680  # histogram of input variable (only signal) - not yet a probability! It's a classifier plot!
681  hist_both = ROOT.TH1F('Both_' + category, 'Input Both (B0) ' +
682  category + ' (binning)', NbinsCategories, -1.0, 1.0)
683  # histogram of input variable (only signal) - not yet a probability! It's a classifier plot!
684  hist_signal = ROOT.TH1F('Signal_' + category, 'Input Signal (B0) ' +
685  category + ' (binning)', NbinsCategories, -1.0, 1.0)
686  # histogram of input variable (only background) - not yet a probability! It's a classifier plot!
687  hist_background = ROOT.TH1F('Background_' + category, 'Input Background (B0bar) ' +
688  category + ' (binning)', NbinsCategories, -1.0, 1.0)
689 
690  # per definition that input is not comparable to the network output, this has to be transformed.
691  # probability output from 0 to 1 (corresponds to net output probability) -> calculation below
692  hist_probB0 = ROOT.TH1F('Probability' + category,
693  'Transformed to probability (B0) (' + category + ')',
694  NbinsCategories, 0.0, 1.0)
695  hist_probB0bar = ROOT.TH1F('ProbabilityB0bar_' + category,
696  'Transformed to probability (B0bar) (' + category + ')',
697  NbinsCategories, 0.0, 1.0)
698  # qp output from -1 to 1 -> transformation below
699  hist_qrB0 = ROOT.TH1F('QR' + category, 'Transformed to qp (B0)(' +
700  category + ')', NbinsCategories, -1.0, 1.0)
701  hist_qrB0bar = ROOT.TH1F('QRB0bar_' + category, 'Transformed to qp (B0bar) (' +
702  category + ')', NbinsCategories, -1.0, 1.0)
703  # histogram for abs(qp), i.e. this histogram contains the r-values -> transformation below
704  # also used to get the number of entries, sorted into 6 bins
705  histo_entries_per_bin = ROOT.TH1F('entries_per_bin_' + category, 'Abs(qp)(B0) (' + category + ')', int(r_size - 2), r_subsample)
706  histo_entries_per_binB0 = ROOT.TH1F('entries_per_bin' + category, 'Abs(qp)(B0) (' +
707  category + ')', int(r_size - 2), r_subsample)
708  histo_entries_per_binB0bar = ROOT.TH1F('entries_per_binB0bar_' + category,
709  'Abs(qp) (B0bar) (' + category + ')', int(r_size - 2), r_subsample)
710 
711  # histogram contains at the end the average r values -> calculation below
712  # sorted into r bins
713  hist_avr_rB0 = ROOT.TH1F('Average_r' + category, 'Average r for B0' +
714  category, int(r_size - 2), r_subsample)
715  hist_avr_rB0bar = ROOT.TH1F('Average_rB0bar_' + category, 'Average r for B0bar' +
716  category, int(r_size - 2), r_subsample)
717 
718  hist_ms_rB0 = ROOT.TH1F('AverageSqrdR' + category, 'Average r sqrd for B0' +
719  category, int(r_size - 2), r_subsample)
720  hist_ms_rB0bar = ROOT.TH1F('AverageSqrdRB0bar_' + category, 'Average r sqrd for B0bar' +
721  category, int(r_size - 2), r_subsample)
722 
723  # ****** TEST OF CALIBRATION ******
724  # for calibration plot we want to have
725  hist_all = ROOT.TH1F('All_' + category, 'Input Signal (B0) and Background (B0Bar)' +
726  category + ' (binning 50)', 50, 0.0, 1.0)
727  hist_calib_B0 = ROOT.TH1F('Calib_' + category, 'Calibration Plot for true B0' +
728  category + ' (binning 50)', 50, 0.0, 1.0)
729 
730  # fill both
731  if category != "KaonNotWeighted" and category != "LambdaNotWeighted":
732  # if category == 'Lambda' or combinerVariable.find('weighted') == -1:
733  tree.Draw('qp' + category + '>>Both_' + category, 'abs(qrMC) == 1.0')
734  # fill signal
735  tree.Draw('qp' + category + '>>Signal_' + category, 'qrMC == 1.0')
736  # fill background
737  tree.Draw('qp' + category + '>>Background_' + category, 'qrMC == -1.0')
738  tree.Draw('qp' + category + '>>All_' + category, 'qrMC!=0')
739  tree.Draw('qp' + category + '>>Calib_' + category, 'qrMC == 1.0')
740 
741  # elif combinerVariable.find('weighted') != -1:
742  # tree.Draw('extraInfo__boWeightedQpOf' + category + '__bc' + '>>Both_' + category, 'abs(qrMC) == 1.0')
743  # tree.Draw('extraInfo__boWeightedQpOf' + category + '__bc' + '>>Signal_' + category, 'qrMC == 1.0')
744  # tree.Draw('extraInfo__boWeightedQpOf' + category + '__bc' + '>>Background_' + category, 'qrMC == -1.0')
745  # tree.Draw('extraInfo__boWeightedQpOf' + category + '__bc' + '>>All_' + category, 'qrMC!=0')
746  # tree.Draw('extraInfo__boWeightedQpOf' + category + '__bc' + '>>Calib_' + category, 'qrMC == 1.0')
747  # category = category + 'W'
748 
749  elif category == "KaonNotWeighted":
750  tree.Draw('extraInfo__boQpOfKaon__bc' + '>>Both_' + category, 'abs(qrMC) == 1.0')
751  tree.Draw('extraInfo__boQpOfKaon__bc' + '>>Signal_' + category, 'qrMC == 1.0')
752  tree.Draw('extraInfo__boQpOfKaon__bc' + '>>Background_' + category, 'qrMC == -1.0')
753  tree.Draw('extraInfo__boQpOfKaon__bc' + '>>All_' + category, 'qrMC!=0')
754  tree.Draw('extraInfo__boQpOfKaon__bc' + '>>Calib_' + category, 'qrMC == 1.0')
755 
756  elif category == "LambdaNotWeighted":
757  tree.Draw('extraInfo__boQpOfLambda__bc' + '>>Both_' + category, 'abs(qrMC) == 1.0')
758  tree.Draw('extraInfo__boQpOfLambda__bc' + '>>Signal_' + category, 'qrMC == 1.0')
759  tree.Draw('extraInfo__boQpOfLambda__bc' + '>>Background_' + category, 'qrMC == -1.0')
760  tree.Draw('extraInfo__boQpOfLambda__bc' + '>>All_' + category, 'qrMC!=0')
761  tree.Draw('extraInfo__boQpOfLambda__bc' + '>>Calib_' + category, 'qrMC == 1.0')
762 
763  hist_calib_B0.Divide(hist_all)
764  # ****** produce the input plots from combiner level ******
765 
766  maxSignal = hist_signal.GetBinContent(hist_signal.GetMaximumBin())
767  maxBackground = hist_background.GetBinContent(hist_background.GetMaximumBin())
768 
769  Ymax = max(maxSignal, maxBackground)
770  Ymax = Ymax + Ymax / 12
771 
772  ROOT.gStyle.SetOptStat(0)
773  with Quiet(ROOT.kError):
774  Canvas = ROOT.TCanvas('Bla', 'TITEL BLA', 1200, 800)
775  Canvas.cd() # activate
776  Canvas.SetLogy()
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.05)
783  hist_signal.GetYaxis().SetTitleSize(0.05)
784  hist_signal.GetXaxis().SetTitleOffset(0.95)
785  hist_signal.GetYaxis().SetTitleOffset(1.1)
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  hist_signal.SetTitle(category + ' category; #it{qp}-Output ; 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  if category == 'MaximumPstar':
803  legend = ROOT.TLegend(0.4, 0.75, 0.6, 0.9)
804  else:
805  legend = ROOT.TLegend(0.6, 0.75, 0.8, 0.9)
806  legend.AddEntry(hist_signal, 'true B0')
807  legend.AddEntry(hist_background, 'true B0bar')
808  legend.SetTextSize(0.05)
809  legend.Draw()
810 
811  Canvas.Update()
812  with Quiet(ROOT.kError):
813  Canvas.SaveAs(workingDirectory + '/' + 'PIC_' + category + '_Input_Combiner.pdf')
814 
815  # ***** TEST OF CALIBRATION ******
816 
817  # initialize some arrays
818  binCounter = int(NbinsCategories + 1)
819  dilutionB02 = array('d', [0] * binCounter)
820  dilutionB0bar2 = array('d', [0] * binCounter)
821  purityB0 = array('d', [0] * binCounter)
822  purityB0bar = array('d', [0] * binCounter)
823  signal = array('d', [0] * binCounter)
824  back = array('d', [0] * binCounter)
825  weight = array('d', [0] * binCounter)
826 
827  for i in range(1, binCounter):
828  # doing the transformation to probabilities
829  signal[i] = hist_signal.GetBinContent(i)
830  back[i] = hist_background.GetBinContent(i)
831  weight[i] = signal[i] + back[i]
832 
833  # avoid dividing by zero
834  if signal[i] + back[i] == 0:
835  purityB0[i] = 0
836  dilutionB02[i] = 0
837  purityB0bar[i] = 0
838  dilutionB0bar2[i] = 0
839  else:
840 
841  purityB0[i] = signal[i] / (signal[i] + back[i])
842  dilutionB02[i] = -1 + 2 * signal[i] / (signal[i] + back[i])
843 
844  purityB0bar[i] = back[i] / (signal[i] + back[i])
845  dilutionB0bar2[i] = -1 + 2 * back[i] / (signal[i] + back[i])
846 
847  # filling histogram with probability from 0 to 1
848  hist_probB0.Fill(purityB0[i], signal[i])
849  hist_probB0bar.Fill(purityB0bar[i], back[i])
850 
851  # filling histogram with qr from -1 to 1
852  hist_qrB0.Fill(dilutionB02[i], signal[i])
853  hist_qrB0bar.Fill(dilutionB0bar2[i], back[i])
854 
855  # filling histogram with abs(qr), i.e. this histogram contains the r-values (not qp)
856  histo_entries_per_binB0.Fill(abs(dilutionB02[i]), signal[i])
857  histo_entries_per_binB0bar.Fill(abs(dilutionB0bar2[i]), back[i])
858  # filling histogram with abs(qr) special weighted - needed for average r calculation
859  hist_avr_rB0.Fill(abs(dilutionB02[i]), abs(dilutionB02[i]) * signal[i])
860  hist_avr_rB0bar.Fill(abs(dilutionB0bar2[i]), abs(dilutionB0bar2[i]) * back[i])
861  # filling histogram with abs(qr)**2 special weighted - needed for std dev of r calculation
862  hist_ms_rB0.Fill(abs(dilutionB02[i]), abs(dilutionB02[i] * dilutionB02[i]) * signal[i])
863  hist_ms_rB0bar.Fill(abs(dilutionB0bar2[i]), abs(dilutionB0bar2[i] * dilutionB02[i]) * back[i])
864 
865  # hist_avr_rB0bar contains now the average r-value
866  hist_avr_rB0.Divide(histo_entries_per_binB0)
867  hist_avr_rB0bar.Divide(histo_entries_per_binB0bar)
868 
869  hist_ms_rB0.Divide(histo_entries_per_binB0)
870  hist_ms_rB0bar.Divide(histo_entries_per_binB0bar)
871  # now calculating the efficiency
872 
873  # calculating number of events
874  total_entriesB0 = total_notTagged / 2
875  total_entriesB0bar = total_notTagged / 2
876  for i in range(1, r_size):
877  total_entriesB0 = total_entriesB0 + histo_entries_per_binB0.GetBinContent(i)
878  total_entriesB0bar = total_entriesB0bar + histo_entries_per_binB0bar.GetBinContent(i)
879  # initializing some arrays
880  tot_eff_effB0 = 0
881  tot_eff_effB0bar = 0
882  uncertainty_eff_effB0 = 0
883  uncertainty_eff_effB0bar = 0
884  uncertainty_eff_effAverage = 0
885  diff_eff_Uncertainty = 0
886  event_fractionB0 = array('f', [0] * r_size)
887  event_fractionB0bar = array('f', [0] * r_size)
888  rvalueB0 = array('f', [0] * r_size)
889  rvalueB0bar = array('f', [0] * r_size)
890  rvalueStdB0 = array('f', [0] * r_size)
891  rvalueStdB0bar = array('f', [0] * r_size)
892  # wvalue = array('f', [0] * r_size)
893  entriesBoth = array('f', [0] * r_size)
894  entriesB0 = array('f', [0] * r_size)
895  entriesB0bar = array('f', [0] * r_size)
896  iEffEfficiencyB0Uncertainty = array('f', [0] * r_size)
897  iEffEfficiencyB0barUncertainty = array('f', [0] * r_size)
898  iDeltaEffEfficiencyUncertainty = array('f', [0] * r_size)
899 
900  for i in range(1, r_size):
901 
902  entriesBoth[i] = entriesB0bar[i] + entriesB0[i]
903  entriesB0[i] = histo_entries_per_binB0.GetBinContent(i)
904  entriesB0bar[i] = histo_entries_per_binB0bar.GetBinContent(i)
905  event_fractionB0[i] = entriesB0[i] / total_entriesB0
906  event_fractionB0bar[i] = entriesB0bar[i] / total_entriesB0bar
907  # print '* Bin ' + str(i) + ' r-value: ' + str(rvalueB0[i]), 'entriesB0: ' +
908  # str(event_fractionB0[i] * 100) + ' % (' + str(entriesB0[i]) + '/' +
909  # str(total_entriesB0) + ')'
910 
911  rvalueB0[i] = hist_avr_rB0.GetBinContent(i)
912  rvalueB0bar[i] = hist_avr_rB0bar.GetBinContent(i)
913 
914  rvalueStdB0[i] = 0
915  rvalueStdB0bar[i] = 0
916 
917  if entriesB0[i] > 1:
918  rvalueStdB0[i] = math.sqrt(abs(hist_ms_rB0.GetBinContent(
919  i) - (hist_avr_rB0.GetBinContent(i))**2)) / math.sqrt(entriesB0[i] - 1)
920 
921  if entriesB0bar[i] > 1:
922  rvalueStdB0bar[i] = math.sqrt(abs(hist_ms_rB0bar.GetBinContent(
923  i) - (hist_avr_rB0bar.GetBinContent(i))**2)) / math.sqrt(entriesB0bar[i] - 1)
924  # wvalue[i] = (1 - rvalueB0[i]) / 2
925 
926  tot_eff_effB0 = tot_eff_effB0 + event_fractionB0[i] * rvalueB0[i] \
927  * rvalueB0[i]
928  tot_eff_effB0bar = tot_eff_effB0bar + event_fractionB0bar[i] * rvalueB0bar[i] \
929  * rvalueB0bar[i]
930 
931  iEffEfficiencyB0Uncertainty[i] = rvalueB0[i] * \
932  math.sqrt((2 * total_entriesB0 * entriesB0[i] * rvalueStdB0[i])**2 +
933  rvalueB0[i]**2 * entriesB0[i] *
934  (total_entriesB0 * (total_entriesB0 - entriesB0[i]) +
935  entriesB0[i] * total_notTagged)) / (total_entriesB0**2)
936  iEffEfficiencyB0barUncertainty[i] = rvalueB0bar[i] * \
937  math.sqrt((2 * total_entriesB0bar * entriesB0bar[i] * rvalueStdB0bar[i])**2 +
938  rvalueB0bar[i]**2 * entriesB0bar[i] *
939  (total_entriesB0bar * (total_entriesB0bar - entriesB0bar[i]) +
940  entriesB0bar[i] * total_notTagged)) / (total_entriesB0bar**2)
941 
942  iDeltaEffEfficiencyUncertainty[i] = math.sqrt(iEffEfficiencyB0Uncertainty[i]**2 + iEffEfficiencyB0barUncertainty[i]**2)
943 
944  diff_eff_Uncertainty = diff_eff_Uncertainty + iDeltaEffEfficiencyUncertainty[i]**2
945 
946  uncertainty_eff_effB0 = uncertainty_eff_effB0 + iEffEfficiencyB0Uncertainty[i]**2
947  uncertainty_eff_effB0bar = uncertainty_eff_effB0bar + iEffEfficiencyB0barUncertainty[i]**2
948 
949  effDiff = tot_eff_effB0 - tot_eff_effB0bar
950  effAverage = (tot_eff_effB0 + tot_eff_effB0bar) / 2
951 
952  uncertainty_eff_effB0 = math.sqrt(uncertainty_eff_effB0)
953  uncertainty_eff_effB0bar = math.sqrt(uncertainty_eff_effB0bar)
954  diff_eff_Uncertainty = math.sqrt(diff_eff_Uncertainty)
955  uncertainty_eff_effAverage = diff_eff_Uncertainty / 2
956  print(
957  '{:<25}'.format("* " + category) + ' B0-Eff=' +
958  '{: 8.2f}'.format(tot_eff_effB0 * 100) + " +-" + '{: 4.2f}'.format(uncertainty_eff_effB0 * 100) +
959  ' %' +
960  ' B0bar-Eff=' +
961  '{: 8.2f}'.format(tot_eff_effB0bar * 100) + " +-" + '{: 4.2f}'.format(uncertainty_eff_effB0bar * 100) +
962  ' %' +
963  ' EffAverage=' +
964  '{: 8.2f}'.format(effAverage * 100) + " +- " + '{:4.2f}'.format(uncertainty_eff_effAverage * 100) + ' %' +
965  ' EffDiff=' +
966  '{: 8.2f}'.format(effDiff * 100) + " +- " + '{:4.2f}'.format(diff_eff_Uncertainty * 100) + ' % *')
967 
968  # hist_signal.Write('', ROOT.TObject.kOverwrite)
969  # hist_background.Write('', ROOT.TObject.kOverwrite)
970  # hist_probB0.Write('', ROOT.TObject.kOverwrite)
971  # hist_probB0bar.Write('', ROOT.TObject.kOverwrite)
972  # hist_qpB0.Write('', ROOT.TObject.kOverwrite)
973  # hist_qpB0bar.Write('', ROOT.TObject.kOverwrite)
974  # hist_absqpB0.Write('', ROOT.TObject.kOverwrite)
975  # hist_absqpB0bar.Write('', ROOT.TObject.kOverwrite)
976  # hist_avr_rB0.Write('', ROOT.TObject.kOverwrite)
977  # hist_avr_rB0bar.Write('', ROOT.TObject.kOverwrite)
978  # hist_all.Write('', ROOT.TObject.kOverwrite)
979  # hist_calib_B0.Write('', ROOT.TObject.kOverwrite)
980  categoriesPerformance.append((category, effAverage, uncertainty_eff_effAverage, effDiff, diff_eff_Uncertainty))
981  with Quiet(ROOT.kError):
982  Canvas.Clear()
983 # if average_eff != 0:
984  # print '* -------------------------------------------------------------------------'
985  # print '* ' + '{: > 8.2f}'.format(average_eff * 100) + ' %' \
986  # + '{:>85}'.format('TOTAL' + ' *')
987 
988 print('* ' +
989  ' *')
990 print('**************************************************************************************************************************' +
991  '************************')
992 print('')
993 print('Table For B2TIP')
994 print('')
995 print(r'\begin{tabular}{ l r r }')
996 print(r'\hline')
997 print(r'Categories & $\varepsilon_\text{eff} \pm \delta\varepsilon_\text{eff} $& ' +
998  r'$\Delta\varepsilon_\text{eff} \pm \delta\Delta\varepsilon_\text{eff}$\\ \hline\hline')
999 for (category, effAverage, uncertainty_eff_effAverage, effDiff, diff_eff_Uncertainty) in categoriesPerformance:
1000  print(
1001  '{:<23}'.format(category) +
1002  ' & $' +
1003  '{: 6.2f}'.format(effAverage * 100) + r" \pm " + '{:4.2f}'.format(uncertainty_eff_effAverage * 100) +
1004  ' $ & $' +
1005  '{: 6.2f}'.format(effDiff * 100) + r" \pm " + '{:4.2f}'.format(diff_eff_Uncertainty * 100) +
1006  r'\ \enskip $ \\')
1007 print(r'\hline')
1008 print(r'\end{tabular}')
1009 B2INFO('qp Output Histograms in pdf format saved at: ' + workingDirectory)
1010 with open("categoriesPerformance.pkl", "wb") as f:
1011  pickle.dump(categoriesPerformance, f)