Belle II Software  light-2403-persian
pythonFlavorTaggerEfficiency.py
1 #!/usr/bin/env python3
2 
3 
10 
11 # ********************************** Flavor Tagging *******************************
12 # * This scripts calculates the effective efficiency of the flavor tagger based *
13 # * on python histograms. It serves as a crosscheck for the script *
14 # * flavorTaggerEfficiency.py, which uses root histograms, and produces nicer plots. *
15 # * This script also produces plots for different amounts of true categories *
16 # * for the combiner outputs. It produces also plots for individual categories *
17 # * checking when they are true categories and when they are not true. *
18 # * True here means that the target (or targets) of the category are found in a *
19 # * certain event. For more information check Sec. 4.5.3 in BELLE2-PTHESIS-2018-003 *
20 # *
21 # ************************************************************************************
22 
23 from matplotlib.colors import colorConverter
24 from matplotlib.collections import PolyCollection
25 from ROOT import PyConfig
26 import glob
27 import sys
28 import math
29 import matplotlib.pyplot as plt
30 import ROOT
31 from basf2 import B2INFO
32 from defaultEvaluationParameters import categories, Quiet, rbins
33 from array import array
34 
35 import numpy as np
36 import matplotlib as mpl
37 mpl.use('Agg')
38 mpl.rcParams.update({'font.size': 22})
39 mpl.rcParams['text.usetex'] = True
40 mpl.rcParams['text.latex.preamble'] = [r"\usepackage{amsmath}"]
41 
42 PyConfig.IgnoreCommandLineOptions = True
43 PyConfig.StartGuiThread = False
44 
45 ROOT.RooMsgService.instance().setGlobalKillBelow(ROOT.RooFit.WARNING)
46 
47 
48 def efficiencyCalculator(data, total_notTagged, SimpleOutput=False, v='notVerbose'):
49 
50  # Calculate Efficiency
51 
52  dataB0 = data[np.where(data[:, 0] > 0)]
53  dataB0bar = data[np.where(data[:, 0] < 0)]
54 
55  totaldataB0 = np.absolute(dataB0[:, 1])
56  totaldataB0bar = np.absolute(dataB0bar[:, 1])
57 
58  rvalueB0 = (np.histogram(totaldataB0, rbins, weights=totaldataB0)[0] / np.histogram(totaldataB0, rbins)[0])
59  rvalueB0MeanSquared = (
60  np.histogram(
61  totaldataB0,
62  rbins,
63  weights=totaldataB0 *
64  totaldataB0)[0] /
65  np.histogram(
66  totaldataB0,
67  rbins)[0])
68 
69  entriesB0 = np.histogram(totaldataB0, rbins)[0]
70 
71  rvalueB0bar = (np.histogram(totaldataB0bar, rbins, weights=totaldataB0bar)[0] / np.histogram(totaldataB0bar, rbins)[0])
72  rvalueB0barMeanSquared = (
73  np.histogram(
74  totaldataB0bar,
75  rbins,
76  weights=totaldataB0bar *
77  totaldataB0bar)[0] /
78  np.histogram(
79  totaldataB0bar,
80  rbins)[0])
81  entriesB0bar = np.histogram(totaldataB0bar, rbins)[0]
82 
83  for row in dataB0:
84  if row[0] == 1 and abs(row[1]) > 1:
85  print(row)
86 
87  total_entries_B0 = totaldataB0.shape[0] + total_notTagged / 2
88  tot_eff_eff_B0 = 0
89  event_fraction_B0 = entriesB0.astype(float) / total_entries_B0
90 
91  total_entries_B0bar = totaldataB0bar.shape[0] + total_notTagged / 2
92  tot_eff_eff_B0bar = 0
93  event_fraction_B0bar = entriesB0bar.astype(float) / total_entries_B0bar
94 
95  tot_eff_eff_Avg = 0
96 
97  arrayShape = len(rbins) - 1
98 
99  event_fractionTotal = np.zeros(arrayShape)
100  event_fractionB0 = np.zeros(arrayShape)
101  event_fractionB0bar = np.zeros(arrayShape)
102  event_fractionDiff = np.zeros(arrayShape)
103  event_fractionDiffUncertainty = np.zeros(arrayShape)
104  event_fractionTotalUncertainty = np.zeros(arrayShape)
105 
106  wvalue = np.zeros(arrayShape)
107  wvalueB0 = np.zeros(arrayShape)
108  wvalueB0bar = np.zeros(arrayShape)
109  wvalueDiff = np.zeros(arrayShape)
110  wvalueDiffUncertainty = np.zeros(arrayShape)
111  wvalueUncertainty = np.zeros(arrayShape)
112 
113  rvalueStdB0 = np.zeros(arrayShape)
114  rvalueStdB0bar = np.zeros(arrayShape)
115 
116  muParam = np.zeros(arrayShape)
117  muParamUncertainty = np.zeros(arrayShape)
118 
119  # Wrong tag fraction from MC counts
120 
121  NwrongB0 = np.histogram(
122  np.absolute(
123  data[np.where((np.sign(data[:, 0]) != np.sign(data[:, 1])) & (data[:, 0] == 1))][:, 1]), rbins)[0].astype(float)
124 
125  wvalueB0 = (NwrongB0 / entriesB0)
126 
127  wvalueB0Uncertainty = np.sqrt(NwrongB0 * (entriesB0 - NwrongB0) / (entriesB0**3))
128 
129  NwrongB0bar = np.histogram(
130  np.absolute(
131  data[np.where((np.sign(data[:, 0]) != np.sign(data[:, 1])) & (data[:, 0] == -1))][:, 1]), rbins)[0].astype(float)
132 
133  wvalueB0bar = (NwrongB0bar / entriesB0bar)
134 
135  wvalueB0barUncertainty = np.sqrt(NwrongB0bar * (entriesB0bar - NwrongB0bar) / (entriesB0bar**3))
136 
137  wvalueDiff = wvalueB0 - wvalueB0bar
138  wvalueDiffUncertainty = np.sqrt(wvalueB0Uncertainty**2 + wvalueB0barUncertainty**2)
139 
140  for i in range(0, len(rbins) - 1):
141 
142  event_fractionB0[i] = entriesB0[i] / total_entries_B0
143  event_fractionB0bar[i] = entriesB0bar[i] / total_entries_B0bar
144 
145  event_fractionTotal[i] = (event_fractionB0[i] + event_fractionB0bar[i]) / 2
146  event_fractionDiff[i] = event_fractionB0[i] - event_fractionB0bar[i]
147 
148  event_fractionDiffUncertainty[i] = math.sqrt(entriesB0[i] *
149  (total_entries_B0 -
150  entriesB0[i]) /
151  total_entries_B0**3 +
152  entriesB0bar[i] *
153  (total_entries_B0bar -
154  entriesB0bar[i]) /
155  total_entries_B0bar**3)
156 
157  event_fractionTotalUncertainty[i] = event_fractionDiffUncertainty[i] / 2
158 
159  muParam[i] = event_fractionDiff[i] / (2 * event_fractionTotal[i])
160  muParamUncertainty[i] = event_fractionDiffUncertainty[i] / (2 * event_fractionTotal[i]) * math.sqrt(muParam[i]**2 + 1)
161 
162  rvalueStdB0[i] = math.sqrt(rvalueB0MeanSquared[i] - (rvalueB0[i])**2) / math.sqrt(entriesB0[i] - 1)
163  rvalueStdB0bar[i] = math.sqrt(rvalueB0barMeanSquared[i] - (rvalueB0bar[i])**2) / math.sqrt(entriesB0bar[i] - 1)
164 
165  # wvalueB0[i] = (1 - rvalueB0[i]) / 2
166  # wvalueB0bar[i] = (1 - rvalueB0bar[i]) / 2
167  # wvalueDiff[i] = wvalueB0[i] - wvalueB0bar[i]
168  # wvalueDiffUncertainty[i] = math.sqrt((rvalueStdB0[i] / 2)**2 + (rvalueStdB0bar[i] / 2)**2)
169  wvalue[i] = (wvalueB0[i] + wvalueB0bar[i]) / 2
170  wvalueUncertainty[i] = wvalueDiffUncertainty[i] / 2
171 
172  tot_eff_eff_B0 += event_fraction_B0[i] * (1 - 2 * wvalueB0[i])**2
173  tot_eff_eff_B0bar += event_fraction_B0bar[i] * (1 - 2 * wvalueB0bar[i])**2
174 
175  eff = event_fractionTotal[i] * (1 - 2 * wvalue[i])**2
176  # (event_fraction_B0[i] * rvalueB0[i] * rvalueB0[i] + event_fraction_B0bar[i] * rvalueB0bar[i] * rvalueB0bar[i]) / 2
177  tot_eff_eff_Avg += eff
178 
179  eff2 = (event_fraction_B0[i] * (1 - 2 * wvalueB0[i])**2 + event_fraction_B0bar[i] * (1 - 2 * wvalueB0bar[i])**2) / 2
180 
181  rval = (rvalueB0[i] + rvalueB0bar[i]) / 2
182  if v == 'verbose':
183  print(
184  "r = ",
185  f"{rval: 6.3f}",
186  "Eff1 = ",
187  f"{eff: 6.3f}",
188  "Eff2 = ",
189  f"{eff2: 6.3f}",
190  "w = ",
191  f"{float(wvalue[i] * 100): 8.3f}",
192  "Delta w = ",
193  f"{float(wvalueDiff[i] * 100): 6.3f}" +
194  " +-" +
195  f"{float(wvalueDiffUncertainty[i] * 100): 1.3f}",
196  "epsilon = ",
197  f"{float(event_fractionTotal[i] * 100): 1.3f}",
198  "Delta epsilon = ",
199  f"{float(event_fractionDiff[i] * 100): 1.3f}" +
200  " +-" +
201  f"{float(event_fractionDiffUncertainty[i] * 100): 1.3f}",
202  "mu = ",
203  f"{float(muParam[i] * 100): 1.3f}" +
204  " +-" +
205  f"{float(muParamUncertainty[i] * 100): 1.3f}")
206 
207  tot_eff_eff = (tot_eff_eff_B0 + tot_eff_eff_B0bar) / 2
208  tot_eff_diff = tot_eff_eff_B0 - tot_eff_eff_B0bar
209  efficiency = 100 * tot_eff_eff_Avg
210  efficiencyDiff = 100 * tot_eff_diff
211  efficiencyB0 = 100 * tot_eff_eff_B0
212  efficiencyB0bar = 100 * tot_eff_eff_B0bar
213 
214  if v == 'verbose':
215  print("Total Efficiency from Avr. rB0 and rB0bar ", f"{float(100 * tot_eff_eff): 6.3f}")
216  print("Total Efficiency from Avr. wB0 and wB0bar ", f"{float(100 * tot_eff_eff_Avg): 6.3f}")
217 
218  if SimpleOutput:
219 
220  return efficiency, efficiencyDiff, efficiencyB0, efficiencyB0bar
221 
222  else:
223 
224  return efficiency, efficiencyDiff, efficiencyB0, efficiencyB0bar, wvalue, wvalueB0, \
225  wvalueB0bar, wvalueDiff, wvalueDiffUncertainty, event_fractionB0, event_fractionB0bar, \
226  event_fractionDiffUncertainty, event_fractionTotalUncertainty, muParam, muParamUncertainty
227 
228 
229 def categoriesEfficiencyCalculator(data, v='notVerbose'):
230 
231  # Calculate Efficiency
232 
233  r_subsample = array('d', [
234  0.0,
235  0.1,
236  0.25,
237  0.5,
238  0.625,
239  0.75,
240  0.875,
241  1.0])
242 
243  hist_aver_rB0 = ROOT.TH1F('AverageR' + category, 'A good one (B0)' +
244  category, 6, r_subsample)
245  hist_aver_rB0bar = ROOT.TH1F('AverageRB0bar_' + category, 'A good one (B0bar)' +
246  category, 6, r_subsample)
247 
248  hist_absqrB0 = ROOT.TH1F('AbsQR' + category, 'Abs(qr)(B0) (' + category + ')', 6, r_subsample)
249  hist_absqrB0bar = ROOT.TH1F('AbsQRB0bar_' + category, 'Abs(qr) (B0bar) (' + category + ')', 6, r_subsample)
250 
251  dataB0 = data[np.where(data[:, 0] > 0)]
252  dataB0bar = data[np.where(data[:, 0] < 0)] * (-1)
253 
254  # -----bins
255  bins = list(range(-25, 27, 1))
256  for i in range(0, len(bins)):
257  bins[i] = float(bins[i]) / 25
258  # ------
259 
260  histoB0 = np.histogram(dataB0, bins=bins)[0]
261  histoB0bar = np.histogram(dataB0bar, bins=bins)[0]
262  histoB0bar = histoB0bar[0:len(histoB0bar) - 1]
263  histoB0bar = histoB0bar[::-1]
264 
265  dilutionB0 = np.zeros(50)
266  dilutionB0bar = np.zeros(50)
267 
268  for i in range(0, 50):
269  if (histoB0[i] + histoB0bar[i]) != 0:
270  dilutionB0[i] = -1 + 2 * histoB0[i] / (histoB0[i] + histoB0bar[i])
271  dilutionB0bar[i] = -1 + 2 * histoB0bar[i] / (histoB0[i] + histoB0bar[i])
272  hist_absqrB0.Fill(abs(dilutionB0[i]), histoB0[i])
273  hist_absqrB0bar.Fill(abs(dilutionB0bar[i]), histoB0bar[i])
274  hist_aver_rB0.Fill(abs(dilutionB0[i]), abs(dilutionB0[i]) * histoB0[i])
275  hist_aver_rB0bar.Fill(abs(dilutionB0bar[i]), abs(dilutionB0bar[i]) * histoB0bar[i])
276 
277  hist_aver_rB0.Divide(hist_absqrB0)
278  hist_aver_rB0bar.Divide(hist_absqrB0bar)
279 
280  tot_entriesB0 = dataB0[:, 1].shape[0]
281  tot_entriesB0bar = dataB0bar[:, 1].shape[0]
282 
283  tot_eff_effB0 = 0
284  tot_eff_effB0bar = 0
285  rvalueB0 = np.zeros(rbins.shape[0])
286  rvalueB0bar = np.zeros(rbins.shape[0])
287  entriesB0 = np.zeros(rbins.shape[0])
288  entriesB0bar = np.zeros(rbins.shape[0])
289  event_fractionB0 = np.zeros(rbins.shape[0])
290  event_fractionB0bar = np.zeros(rbins.shape[0])
291 
292  for i in range(1, rbins.shape[0]):
293  rvalueB0[i] = hist_aver_rB0.GetBinContent(i)
294  rvalueB0bar[i] = hist_aver_rB0bar.GetBinContent(i)
295  entriesB0[i] = hist_absqrB0.GetBinContent(i)
296  entriesB0bar[i] = hist_absqrB0bar.GetBinContent(i)
297  event_fractionB0[i] = entriesB0[i] / tot_entriesB0
298  event_fractionB0bar[i] = entriesB0bar[i] / tot_entriesB0bar
299  tot_eff_effB0 = tot_eff_effB0 + event_fractionB0[i] * rvalueB0[i] \
300  * rvalueB0[i]
301  tot_eff_effB0bar = tot_eff_effB0bar + event_fractionB0bar[i] * rvalueB0bar[i] \
302  * rvalueB0bar[i]
303  effDiff = tot_eff_effB0 - tot_eff_effB0bar
304  effAverage = (tot_eff_effB0 + tot_eff_effB0bar) / 2
305 
306  efficiency = 100 * effAverage
307  efficiencyDiff = 100 * effDiff
308  efficiencyB0 = 100 * tot_eff_effB0
309  efficiencyB0bar = 100 * tot_eff_effB0bar
310 
311  return efficiency, efficiencyDiff, efficiencyB0, efficiencyB0bar
312 
313 
314 ROOT.gROOT.SetBatch(True)
315 
316 if len(sys.argv) != 3:
317  sys.exit("Must provide 2 arguments: [input_sim_file] or ['input_sim_file*'] with wildcards and [treeName]"
318  )
319 workingFile = sys.argv[1]
320 workingFiles = glob.glob(str(workingFile))
321 treeName = str(sys.argv[2])
322 
323 if len(workingFiles) < 1:
324  sys.exit("No file name or file names " + str(workingFile) + " found.")
325 
326 workingDirectory = '.'
327 
328 methods = []
329 
330 
331 tree = ROOT.TChain(treeName)
332 
333 mcstatus = array('d', [-511.5, 0.0, 511.5])
334 ROOT.TH1.SetDefaultSumw2()
335 
336 for iFile in workingFiles:
337  tree.AddFile(iFile)
338 
339 totalBranches = []
340 for branch in tree.GetListOfBranches():
341  totalBranches.append(branch.GetName())
342 
343 existsDeltaT = False
344 
345 if 'FBDT_qrCombined' in totalBranches:
346  methods.append("FBDT")
347 
348 if 'FANN_qrCombined' in totalBranches:
349  methods.append("FANN")
350 
351 if 'DNN_qrCombined' in totalBranches:
352  methods.append("DNN")
353 
354 if 'DeltaT' in totalBranches:
355  existsDeltaT = True
356 
357 usedCategories = []
358 for cat in categories:
359  catBranch = 'qp' + cat
360  if catBranch in totalBranches:
361  usedCategories.append(cat)
362 
363 YmaxForQrPlot = 0
364 YmaxForDeltaTPlot = 0
365 
366 total_notTagged = 0
367 
368 for method in methods:
369 
370  if method == "FBDT":
371  label = method
372  elif method == "DNN":
373  label = method
374  elif method == "FANN":
375  label = "MLP"
376 
377  with Quiet(ROOT.kError):
378  qpMaximumPstar = ROOT.RooRealVar('qpMaximumPstar', 'qpMaximumPstar', 0, -2.1, 1.1)
379  qrCombined = ROOT.RooRealVar('' + method + '_qrCombined', '' + method + '_qrCombined', 0, -1.1, 1.1)
380  qrMC = ROOT.RooRealVar('qrMC', 'qrMC', 0, -10.0, 10.0)
381  isSignal = ROOT.RooRealVar('isSignal', 'isSignal', 0, -10.0, 10.0)
382 
383  DeltaT = ROOT.RooRealVar("DeltaT", "DeltaT", 0., -100000., 100000.)
384 
385  argSet = ROOT.RooArgSet(qrCombined, qrMC, isSignal, qpMaximumPstar)
386 
387  if existsDeltaT:
388 
389  argSet.add(DeltaT)
390 
391  for category in usedCategories:
392 
393  exec(
394  f"hasTrueTarget + {category} = " +
395  "ROOT.RooRealVar('hasTrueTarget' + category, 'hasTrueTarget' + category, 0, -2.1, 1.1)"
396  )
397  # exec("%s = %s" % ("isTrueCategory" + category,
398  # "ROOT.RooRealVar('isRightCategory' + category, 'isRightCategory' +
399  # category, 0, -1.0, 1.1)"))
400  exec(f"argSet.add(hasTrueTarget{category})")
401 
402  rooDataSet = ROOT.RooDataSet("data", "data", tree, argSet, "")
403 
404  dataAll = []
405  numberOfTrueCatagories = []
406  deltaT = []
407  # data = []
408 
409  for i in range(rooDataSet.numEntries()):
410  row = rooDataSet.get(i)
411  tqrCombined = row.getRealValue('' + method + '_qrCombined', 0, ROOT.kTRUE)
412  tqrMC = row.getRealValue('qrMC', 0, ROOT.kTRUE)
413  tisSignal = row.getRealValue('isSignal', 0, ROOT.kTRUE)
414 
415  if tisSignal == 1 and tqrCombined < -1 and abs(tqrMC) == 0:
416  total_notTagged += 1
417 
418  # tqpMaximumPstar = row.getRealValue('qpMaximumPstar', 0, ROOT.kTRUE)
419 
420  if tqrCombined >= 1:
421  tqrCombined = 0.999
422 
423  if tqrCombined <= -1:
424  tqrCombined = -0.999
425 
426  tNumberOfTrueCategories = 0
427 
428  for category in usedCategories:
429 
430  if category == "MaximumPstar":
431  continue
432  # tisTrueCategory = row.getRealValue('isRightCategory' + category, 0, ROOT.kTRUE)
433  thasTrueTarget = row.getRealValue('hasTrueTarget' + category, 0, ROOT.kTRUE)
434 
435  if thasTrueTarget > 0:
436  tNumberOfTrueCategories += 1
437  if abs(tqrMC) == 1: # !=0:
438  # data.append([tqrMC, tqrCombined])
439  numberOfTrueCatagories.append(tNumberOfTrueCategories)
440  # if tNumberOfTrueCategories == 0:
441  # dataNoTrueCategory.append([tqrMC, tqrCombined])
442  if existsDeltaT:
443  tDT = row.getRealValue("DeltaT", 0, ROOT.kTRUE)
444  dataAll.append([tqrMC, tqrCombined, tNumberOfTrueCategories, tDT])
445  else:
446  dataAll.append([tqrMC, tqrCombined, tNumberOfTrueCategories])
447 
448  # data = np.array(data)
449  dataAll = np.array(dataAll)
450  dataNoTrueCategory = dataAll[(np.where(dataAll[:, 2] == 0))][:, 0:4]
451  data = dataAll[:, 0:4]
452  numberOfTrueCatagories = np.array(numberOfTrueCatagories)
453  # Calculate Efficiency
454 
455  dataB0 = data[np.where(data[:, 0] > 0)]
456  dataB0bar = data[np.where(data[:, 0] < 0)]
457 
458  efficiency, efficiencyDiff, efficiencyB0, efficiencyB0bar, wvalue, wvalueB0, wvalueB0bar, wvalueDiff, wvalueDiffUncertainty, \
459  event_fractionB0, event_fractionB0bar, event_fractionDiffUncertainty, event_fractionTotalUncertainty, muParam, \
460  muParamUncertainty = efficiencyCalculator(data, total_notTagged, False, 'verbose')
461 
462  efficiency = f'{efficiency:.2f}'
463  efficiencyDiff = f'{efficiencyDiff:.2f}'
464  efficiencyB0 = f'{efficiencyB0:.2f}'
465  efficiencyB0bar = f'{efficiencyB0bar:.2f}'
466 
467  N = data.shape[0]
468  NB0 = dataB0.shape[0] / N
469  NB0 = 100 * NB0
470  NB0 = f'{NB0:.2f}'
471  NB0bar = dataB0bar.shape[0] / N
472  NB0bar = 100 * NB0bar
473  NB0bar = f'{NB0bar:.2f}'
474  NnoTrueCat = dataNoTrueCategory.shape[0]
475  FracNoFTInfo = NnoTrueCat / N
476 
477  print('The total efficiency for B0: ', efficiencyB0, '%')
478  print('The total efficiency for B0bar: ', efficiencyB0bar, '%')
479  print('The total efficiency for method ' + method + ' is: ', efficiency, '%')
480  print('The total efficiency Diff for method ' + method + ' is: ', efficiencyDiff, '%')
481 
482  print('N : ', N)
483  print('Percentage B0/ N: ', NB0, '%')
484  print('Percentage B0bar / N : ', NB0bar, '%')
485 
486  print(" ")
487  print("N without trueCats = ", NnoTrueCat)
488  print("Fraction Of Events without trueCats Info = ", f'{float(FracNoFTInfo * 100): 4.2f}', "% ")
489 
490  # ----------------Delta T Plots -------------------------------------------------------------------------------
491 
492  if existsDeltaT:
493 
494  # -----DeltaTbins
495  DeltaTbins = list(range(-1000, 1000, 1))
496  for i in range(0, len(DeltaTbins)):
497  DeltaTbins[i] = float(DeltaTbins[i]) / 10
498  # ------
499 
500  ax = plt.axes([0.21, 0.15, 0.75, 0.8])
501  y, x, _ = ax.hist(data[:, 3], bins=DeltaTbins, facecolor='r', histtype='stepfilled', edgecolor='r') # histtype='step'
502 
503  YDTmax = y.max()
504  YDTmax = YDTmax - YDTmax / 5
505 
506  if YmaxForDeltaTPlot < YDTmax:
507  YmaxForDeltaTPlot = YDTmax
508 
509  figTmc = plt.figure(1, figsize=(11, 8))
510  axTmc = plt.axes([0.21, 0.15, 0.76, 0.8])
511 
512  axTmc.hist(data[np.where(data[:, 0] == -1)][:, 3], bins=DeltaTbins, histtype='step',
513  edgecolor='b', linewidth=4.5, linestyle='dashed', label=r'$\bar{B}^0_{\rm tag}$') # hatch='.', 'dotted'
514 
515  axTmc.hist(data[np.where(data[:, 0] == 1)][:, 3], bins=DeltaTbins, histtype='step',
516  edgecolor='r', linewidth=4, alpha=0.9, label=r'$B^0_{\rm tag}$') # histtype='step', hatch='\\'
517 
518  p1, = axTmc.plot([], label=r'$B^0_{\rm tag}$', linewidth=4, linestyle='solid', alpha=0.9, c='r')
519  p2, = axTmc.plot([], label=r'$\bar{B}^0_{\rm tag}$', linewidth=4.5, linestyle='dashed', c='b')
520 
521  axTmc.set_ylabel(r'${\rm Number\ \ of\ \ Events /\ (\, 0.1\ ps\, )}$', fontsize=35)
522  axTmc.set_xlabel(r'$\Delta t\ /\ {\rm ps}$', fontsize=35)
523 
524  plt.xticks([-10, -7.5, -5, -2.5, 0, 2.5, 5, 7.5, 10], [r'$-10$', r'', r'$5$',
525  r'', r'$0$', r'', r'$5$', r'', r'$10$'], rotation=0, size=40)
526  axTmc.tick_params(axis='y', labelsize=35)
527  axTmc.legend([p1, p2], [r'$B^0_{\rm tag}$', r'$\bar{B}^0_{\rm tag}$'], prop={'size': 35})
528  axTmc.grid(True)
529  axTmc.set_ylim(0, YmaxForDeltaTPlot)
530  axTmc.set_xlim(-10, 10)
531  plt.savefig(workingDirectory + '/' + method + 'DeltaTMC.pdf')
532  figTmc.clear()
533 
534  figTft = plt.figure(2, figsize=(11, 8))
535  axTft = plt.axes([0.21, 0.15, 0.76, 0.8])
536 
537  axTft.hist(data[np.where(data[:, 1] < 0)][:, 3], bins=DeltaTbins, histtype='step',
538  edgecolor='b', linewidth=4.5, linestyle='dashed', label=r'$\bar{B}^0_{\rm tag}$') # hatch='.', 'dotted'
539 
540  axTft.hist(data[np.where(data[:, 1] > 0)][:, 3], bins=DeltaTbins, histtype='step',
541  edgecolor='r', linewidth=4, alpha=0.9, label=r'$B^0_{\rm tag}$') # histtype='step', hatch='\\'
542 
543  p1, = axTft.plot([], label=r'$B^0_{\rm tag}$', linewidth=4, linestyle='solid', alpha=0.9, c='r')
544  p2, = axTft.plot([], label=r'$\bar{B}^0_{\rm tag}$', linewidth=4.5, linestyle='dashed', c='b')
545 
546  axTft.set_ylabel(r'${\rm Number\ \ of\ \ Events /\ (\, 0.1\ ps\, )}$', fontsize=35)
547  axTft.set_xlabel(r'$\Delta t\ /\ {\rm ps}$', fontsize=35)
548 
549  plt.xticks([-10, -7.5, -5, -2.5, 0, 2.5, 5, 7.5, 10], [r'$-10$', r'', r'$5$',
550  r'', r'$0$', r'', r'$5$', r'', r'$10$'], rotation=0, size=40)
551  axTft.tick_params(axis='y', labelsize=35)
552  axTft.legend([p1, p2], [r'$B^0_{\rm tag}$', r'$\bar{B}^0_{\rm tag}$'], prop={'size': 35})
553  axTft.grid(True)
554  axTft.set_ylim(0, YmaxForDeltaTPlot)
555  axTft.set_xlim(-10, 10)
556  plt.savefig(workingDirectory + '/' + method + 'DeltaTFT.pdf')
557  figTft.clear()
558 
559  # ----------------QR and Calibration Plots ---------------------------------------------------------------------
560 
561  # -----bins
562  bins = list(range(-50, 51, 1))
563  for i in range(0, len(bins)):
564  bins[i] = float(bins[i]) / 50
565  # ------
566 
567  # -----bins2
568  bins2 = list(range(-25, 26, 1))
569  for i in range(0, len(bins2)):
570  bins2[i] = float(bins2[i]) / 25
571  # ------
572 
573  fig1 = plt.figure(2, figsize=(11, 8))
574  ax1 = plt.axes([0.21, 0.15, 0.75, 0.8])
575  y, x, _ = ax1.hist(data[:, 1], bins=bins, facecolor='r', histtype='stepfilled', edgecolor='r') # histtype='step'
576 
577  Ymax = y.max()
578  Ymax = Ymax + Ymax / 4
579 
580  if YmaxForQrPlot < Ymax:
581  YmaxForQrPlot = Ymax
582 
583  ax1.set_ylabel(r'${\rm Number\ \ of\ \ Events /\ (\, 0.02\, )}$', fontsize=35)
584  ax1.set_xlabel(r'$(q\cdot r)_{\rm ' + label + ' }$', fontsize=42)
585  plt.xticks([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1],
586  [r'$-1$', r'', r'$-0.5$', r'', r'$0$', r'', r'$0.5$', r'', r'$1$'], rotation=0, size=20)
587  ax1.tick_params(axis='y', labelsize=20)
588  ax1.grid(True)
589  ax1.set_ylim(0, YmaxForQrPlot)
590  plt.savefig(workingDirectory + '/' + method + 'QR_Output.pdf')
591  fig1.clear()
592 
593  # --------------------------- QR Plot ----------------- ----------------------
594 
595  fig2 = plt.figure(3, figsize=(11, 8))
596  ax2 = plt.axes([0.21, 0.15, 0.76, 0.8])
597 
598  ax2.hist(data[np.where(abs(data[:, 0]) == 1)][:, 1], bins=bins, histtype='step', edgecolor='k',
599  linewidth=4, linestyle='dashdot', alpha=0.9, label=r'${\rm Both}$')
600 
601  ax2.hist(data[np.where(data[:, 0] == 1)][:, 1], bins=bins, histtype='step', edgecolor='r',
602  linewidth=4, alpha=0.9, label=r'$B^0$') # histtype='step', hatch='\\'
603 
604  ax2.hist(data[np.where(data[:, 0] == -1)][:, 1], bins=bins, histtype='step',
605  edgecolor='b', linewidth=4.5, linestyle='dashed', label=r'$\bar{B}^0$') # hatch='.', 'dotted'
606 
607  p3, = ax2.plot([], label=r'${\rm Both}$', linewidth=4, linestyle='dashdot', alpha=0.9, c='k')
608 
609  p2, = ax2.plot([], label=r'$\bar{B}^0$', linewidth=4.5, linestyle='dashed', c='b')
610 
611  p1, = ax2.plot([], label=r'$B^0$', linewidth=4, linestyle='solid', alpha=0.9, c='r')
612 
613  ax2.set_ylabel(r'${\rm Number\ \ of\ \ Events\ /\ (\, 0.02\, )}$', fontsize=35)
614  ax2.set_xlabel(r'$(q\cdot r)_{\rm ' + label + '}$', fontsize=42)
615  # plt.xticks([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1], [r'$-1$', r'',
616  # r'$-0.5$', r'', r'$0$', r'', r'$0.5$', r'', r'$1$'], rotation=0,
617  # size=40)
618  ax2.tick_params(axis='y', labelsize=35)
619  ax2.legend([p3, p2, p1], [r'${\rm Both}$', r'$\bar{B}^0$', r'$B^0$'], prop={'size': 35}, loc=9)
620  ax2.grid(True)
621  plt.xticks([-1, -0.875, -0.75, -0.625, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 0.625, 0.75, 0.875, 1],
622  [r'${-1.0}$', r'', r'', r'', r'${-0.5}$', r'', r'', r'$0$', r'',
623  r'', r'$0.5$', r'', r'', r'', r'$1.0$'], rotation=0, size=35)
624 
625  ax2.set_ylim(0, YmaxForQrPlot)
626  yLocs, yLabels = plt.yticks()
627  plt.savefig(workingDirectory + '/' + method + 'QR_B0bar.pdf')
628  fig2.clear()
629 
630  # --- No true category plot
631 
632  fig2b = plt.figure(4, figsize=(10, 8))
633  ax2b = plt.axes([0.21, 0.15, 0.76, 0.8])
634  p2b = ax2b.hist(dataNoTrueCategory[np.where(dataNoTrueCategory[:, 0] == -1)][:, 1], bins=bins, histtype='step',
635  edgecolor='b', linewidth=3.5, linestyle='dotted', label=r'$\bar{B}^0$') # hatch='.',
636  p1b = ax2b.hist(
637  dataNoTrueCategory[
638  np.where(
639  dataNoTrueCategory[
640  :,
641  0] == 1)][
642  :,
643  1],
644  bins=bins,
645  histtype='step',
646  edgecolor='r',
647  linewidth=2,
648  alpha=0.9,
649  label=r'$B^0$') # histtype='step', hatch='\\'
650  ax2b.set_ylabel(r'${\rm Number\ \ of\ \ Events\ /\ (\, 0.02\, )}$', fontsize=35)
651  ax2b.set_xlabel(r'$(q\cdot r)_{\rm ' + label + '}$', fontsize=42)
652  plt.xticks([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1],
653  [r'$-1$', r'', r'$-0.5$', r'', r'$0$', r'', r'$0.5$', r'', r'$1$'], rotation=0, size=20)
654  ax2b.tick_params(axis='y', labelsize=20)
655  ax2b.legend([r'$\bar{B}^0$', r'$B^0$'], prop={'size': 20})
656  ax2b.grid(True)
657  # ax2b.set_ylim(0, Ymax)
658  plt.savefig(workingDirectory + '/' + method + 'QR_B0barNoTrueCategory.pdf')
659  fig2b.clear()
660 
661  # rbins = np.array([0.0, 0.1, 0.25, 0.5, 0.625, 0.75, 0.875, 1.0])
662 
663  # -----bins
664  # rbins=range(0,100,1)
665  # for i in xrange(len(rbins)):
666  # bins[i]=float(rbins[i])/100
667  # ------
668 
669  # ------------------------ Calibration plots ----------------------##
670 
671  # data_entries = np.histogram(np.absolute(data[:, 1]), rbins)[0]
672 
673  data_entriesB0 = np.histogram(np.absolute(data[np.where(data[:, 0] == 1)][:, 1]), rbins)[0]
674 
675  rvalueB0 = (
676  np.histogram(
677  np.absolute(data[np.where(data[:, 0] == 1)][:, 1]),
678  rbins,
679  weights=np.absolute(data[np.where(data[:, 0] == 1)][:, 1]))[0] / data_entriesB0)
680  delta_rvalueB0 = (
681  np.histogram(
682  np.absolute(data[np.where(data[:, 0] == 1)][:, 1]), rbins, weights=np.absolute(data[np.where(data[:, 0] == 1)][:, 1]) *
683  np.absolute(data[np.where(data[:, 0] == 1)][:, 1]))[0] / data_entriesB0)
684  delta_rvalueB0 = np.sqrt((delta_rvalueB0 - rvalueB0 * rvalueB0) / (data_entriesB0))
685 
686  wrongTaggedB0 = np.histogram(
687  np.absolute(
688  data[np.where((np.sign(data[:, 0]) != np.sign(data[:, 1])) & (data[:, 0] == 1))][:, 1]), rbins)[0].astype(float)
689 
690  wrong_fractionB0 = (wrongTaggedB0 / data_entriesB0)
691 
692  delta_wrong_fractionB0 = np.sqrt(wrongTaggedB0 * (data_entriesB0 - wrongTaggedB0) / (data_entriesB0**3))
693 
694  rmvaluesB0 = 1 - 2 * wrong_fractionB0
695  delta_rmvaluesB0 = 2 * delta_wrong_fractionB0
696 
697  data_entriesB0bar = np.histogram(np.absolute(data[np.where(data[:, 0] == -1)][:, 1]), rbins)[0]
698 
699  rvalueB0bar = (
700  np.histogram(
701  np.absolute(data[np.where(data[:, 0] == -1)][:, 1]),
702  rbins,
703  weights=np.absolute(data[np.where(data[:, 0] == -1)][:, 1]))[0] / data_entriesB0bar)
704  delta_rvalueB0bar = (
705  np.histogram(
706  np.absolute(data[np.where(data[:, 0] == -1)][:, 1]), rbins,
707  weights=np.absolute(data[np.where(data[:, 0] == -1)][:, 1]) *
708  np.absolute(data[np.where(data[:, 0] == -1)][:, 1]))[0] / data_entriesB0bar)
709  delta_rvalueB0bar = np.sqrt((delta_rvalueB0bar - rvalueB0bar * rvalueB0bar) / (data_entriesB0bar))
710 
711  wrongTaggedB0bar = np.histogram(
712  np.absolute(
713  data[np.where((np.sign(data[:, 0]) != np.sign(data[:, 1])) & (data[:, 0] == -1))][:, 1]), rbins)[0].astype(float)
714 
715  wrong_fractionB0bar = (wrongTaggedB0bar / data_entriesB0bar)
716 
717  delta_wrong_fractionB0bar = np.sqrt(wrongTaggedB0bar * (data_entriesB0bar - wrongTaggedB0bar) / (data_entriesB0bar**3))
718 
719  rmvaluesB0bar = 1 - 2 * wrong_fractionB0bar
720  delta_rmvaluesB0bar = 2 * delta_wrong_fractionB0bar
721 
722  rvalue = (rvalueB0 + rvalueB0bar) / 2
723  rmvalues = (rmvaluesB0 + rmvaluesB0bar) / 2
724 
725  delta_rvalue = np.sqrt(delta_rvalueB0**2 + delta_rvalueB0bar**2) / 2
726  delta_rmvalues = np.sqrt(delta_rmvaluesB0**2 + delta_rmvaluesB0bar**2) / 2
727 
728  fig3a = plt.figure(5, figsize=(11, 8))
729  ax3a = plt.axes([0.21, 0.15, 0.74, 0.8]) # plt.axes([0.14, 0.1, 0.8, 0.8])
730 
731  line = range(0, 2, 1)
732  ax3a.plot(line, line, linewidth=4, color='r')
733 
734  (p1ca, capsB0, _) = ax3a.errorbar(rvalueB0, np.absolute(rmvaluesB0), xerr=None, yerr=None,
735  # xerr=delta_rvalueB0, yerr=delta_rmvaluesB0,
736  elinewidth=4, capsize=10, ecolor='r', fmt='o', mfc='r',
737  mec='r', markersize=14) # histtype='step'
738 
739  (p2ca, capsB0bar, _) = ax3a.errorbar(rvalueB0bar, np.absolute(rmvaluesB0bar), xerr=None, yerr=None,
740  # xerr=delta_rvalueB0bar, yerr=delta_rmvaluesB0bar,
741  elinewidth=4, capsize=10, ecolor='b', fmt='o', mfc='b',
742  mec='b', markersize=14) # histtype='step'
743 
744  for cap in capsB0:
745  cap.set_color('r')
746  cap.set_alpha(0.9)
747  cap.set_markeredgewidth(4)
748 
749  for cap in capsB0bar:
750  cap.set_color('b')
751  cap.set_markeredgewidth(4)
752 
753  (p3ca, caps, _) = ax3a.errorbar(rvalue, rmvalues, xerr=None, yerr=None,
754  # xerr=delta_rvalue, yerr=delta_rmvalues,
755  elinewidth=4, capsize=10, ecolor='k', fmt='o', mfc='k', mec='k', markersize=14)
756 
757  for cap in caps:
758  cap.set_color('k')
759  cap.set_markeredgewidth(4)
760 
761  ax3a.legend([p3ca, p2ca, p1ca], [r'${\rm Average}$', r'$\bar{B}^0$', r'$B^0$'], prop={'size': 35}, loc=2)
762 
763  ax3a.set_ylabel(r'$r_{\rm \small MC} = \vert 1-2\cdot w_{\small \rm MC} \vert $', fontsize=42)
764  ax3a.set_xlabel(r'$\langle r_{\rm ' + label + '}$' + r'$\rangle$', fontsize=42)
765  # plt.xticks([0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1], [r'$0$', r'', r'$0.2$',
766  # r'', r'$0.4$', r'', r'$0.6$', r'', r'$0.8$', r'', r'$1.0$'],
767  # rotation=0, size=35)
768  plt.xticks([0, 0.1, 0.25, 0.5, 0.625, 0.75, 0.875, 1],
769  [r'$0$', r'$0.1$', r'$0.25$', r'$0.5$', r'', r'$0.75$', r'', r'$1.0$'], rotation=0, size=35)
770  plt.yticks([0, 0.1, 0.25, 0.5, 0.625, 0.75, 0.875, 1],
771  [r'$0$', r'$0.1$', r'$0.25$', r'$0.5$', r'', r'$0.75$', r'', r'$1.0$'], rotation=0, size=35)
772  ax3a.set_xlim(0, 1)
773  ax3a.set_ylim(0, 1)
774  # ax3.tick_params(axis='y', labelsize=35)
775  ax3a.xaxis.grid(True, linewidth=2)
776  plt.savefig(workingDirectory + '/' + method + 'CalibrationB0B0bar.pdf')
777  fig3a.clear()
778 
779  fig3 = plt.figure(6, figsize=(11, 8))
780  ax3 = plt.axes([0.21, 0.15, 0.74, 0.8]) # plt.axes([0.14, 0.1, 0.8, 0.8])
781 
782  line = range(0, 2, 1)
783  ax3.plot(line, line, linewidth=4, color='r')
784 
785  (_, caps, _) = ax3.errorbar(rvalue, rmvalues, xerr=None, yerr=None,
786  # xerr=delta_rvalue, yerr=delta_rmvalues,
787  elinewidth=4, capsize=10, ecolor='b', fmt='o', markersize=14) # histtype='step'
788  for cap in caps:
789  # cap.set_color('red')
790  cap.set_markeredgewidth(4)
791 
792  ax3.set_ylabel(r'$r_{\rm \small MC} = 1-2\cdot w_{\small \rm MC}$', fontsize=42)
793  ax3.set_xlabel(r'$\langle r_{\rm ' + label + '}$' + r'$\rangle$', fontsize=42)
794  # plt.xticks([0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1], [r'$0$', r'', r'$0.2$',
795  # r'', r'$0.4$', r'', r'$0.6$', r'', r'$0.8$', r'', r'$1.0$'],
796  # rotation=0, size=35)
797  plt.xticks([0, 0.1, 0.25, 0.5, 0.625, 0.75, 0.875, 1],
798  [r'$0$', r'$0.1$', r'$0.25$', r'$0.5$', r'', r'$0.75$', r'', r'$1.0$'], rotation=0, size=35)
799  plt.yticks([0, 0.1, 0.25, 0.5, 0.625, 0.75, 0.875, 1],
800  [r'$0$', r'$0.1$', r'$0.25$', r'$0.5$', r'', r'$0.75$', r'', r'$1.0$'], rotation=0, size=35)
801  ax3.set_xlim(0, 1)
802  ax3.set_ylim(0, 1)
803  # ax3.tick_params(axis='y', labelsize=35)
804  ax3.xaxis.grid(True, linewidth=2)
805  plt.savefig(workingDirectory + '/' + method + 'Calibration.pdf')
806  fig3.clear()
807 
808  # -- R output plot
809 
810  fig4 = plt.figure(7, figsize=(11, 8))
811  ax4 = plt.axes([0.21, 0.15, 0.76, 0.8])
812  ax4.hist(np.absolute(data[:, 1]), bins=bins, facecolor='r', histtype='stepfilled', edgecolor='r') # histtype='step'
813  ax4.set_ylabel(r'${\rm Number\ \ of\ \ Events\ /\ (\, 0.02\, )}$', fontsize=35)
814  ax4.set_xlabel(r'$r_{\rm ' + label + ' }$', fontsize=42)
815  ax4.tick_params(axis='y', labelsize=30)
816  # ax4.axvline(0.1, linestyle='--', color='k', linewidth=2)
817  # ax4.axvline(0.25, linestyle='--', color='k', linewidth=2)
818  # ax4.axvline(0.5, linestyle='--', color='k', linewidth=2)
819  # ax4.axvline(0.625, linestyle='--', color='k', linewidth=2)
820  # ax4.axvline(0.75, linestyle='--', color='k', linewidth=2)
821  # ax4.axvline(0.875, linestyle='--', color='k', linewidth=2)
822  ax4.set_xlim(0, 1)
823  # ax4.set_ylim(0, 1.4)
824  plt.xticks([0, 0.1, 0.25, 0.5, 0.625, 0.75, 0.875, 1],
825  [r'$0$', r'$0.1$', r'$0.25$', r'$0.5$', r'$0.625$', r'$0.75$', r'$0.875$', r'$1$'], rotation=0, size=30)
826  ax4.yaxis.grid(True)
827  ax4.xaxis.grid(True, linewidth=2)
828  plt.savefig(workingDirectory + '/' + method + 'R_Output.pdf')
829  fig4.clear()
830 
831  # ----
832 
833  # QR Plot folded
834 
835  fig2a = plt.figure(8, figsize=(11, 8))
836  ax2a = plt.axes([0.21, 0.15, 0.76, 0.8])
837 
838  ax2a.hist(data[np.where((abs(data[:, 0]) == 1) & (data[:, 1] >= 0))][:, 1], bins=bins, histtype='step', edgecolor='k',
839  linewidth=4, linestyle='dashdot', alpha=0.9, label=r'$q \geq 0 $')
840 
841  ax2a.hist(data[np.where((abs(data[:, 0]) == 1) & (data[:, 1] < 0))][:, 1] * (-1), bins=bins, histtype='step', edgecolor='g',
842  linewidth=4, linestyle='dashed', alpha=0.9, label=r'$q < 0 $')
843 
844  p3a, = ax2a.plot([], label=r'$q \geq 0 $', linewidth=4, linestyle='dashdot', alpha=0.9, c='k')
845 
846  p2a, = ax2a.plot([], label=r'$q < 0 $', linewidth=4.5, linestyle='dashed', c='g')
847 
848  ax2a.set_ylabel(r'${\rm Number\ \ of\ \ Events\ /\ (\, 0.02\, )}$', fontsize=35)
849  ax2a.set_xlabel(r'$(q\cdot q\cdot r)_{\rm ' + label + '}$', fontsize=42)
850  # plt.xticks([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1], [r'$-1$', r'',
851  # r'$-0.5$', r'', r'$0$', r'', r'$0.5$', r'', r'$1$'], rotation=0,
852  # size=40)
853  ax2a.tick_params(axis='y', labelsize=35)
854  ax2a.legend([p3a, p2a], [r'$q \geq 0 $', r'$q < 0 $'], prop={'size': 35}, loc=2)
855  ax2a.grid(True)
856  plt.xticks([0, 0.1, 0.25, 0.5, 0.625, 0.75, 0.875, 1],
857  [r'$0$', r'$0.1$', r'$0.25$', r'$0.5$', r'', r'$0.75$', r'', r'$1.0$'], rotation=0, size=35)
858 
859  ax2a.set_xlim(0, 1)
860  ax2a.set_ylim(0, YmaxForQrPlot)
861  plt.savefig(workingDirectory + '/' + method + 'Folded_QR_Both.pdf')
862  fig2a.clear()
863 
864  # QR Plot folded for each MC flavor
865 
866  # fig2c = plt.figure(2, figsize=(11, 8))
867  # ax2c = plt.axes([0.21, 0.15, 0.76, 0.8])
868  #
869  # nB0bar, nbins, patches = ax2c.hist(data[np.where(data[:, 0] == -1)][:, 1]*(-1), bins=bins, histtype='step',
870  # edgecolor='b', linewidth=4.5, linestyle='dashed', label=r'$\bar{B}^0$') # hatch='.', 'dotted'
871  #
872  # nB0, nbins, patches = ax2c.hist(data[np.where(data[:, 0] == 1)][:, 1], bins=bins, histtype='step', edgecolor='r',
873  # linewidth=4, alpha=0.9, label=r'$B^0$') # histtype='step', hatch='\\'
874  #
875  # p2c, = ax2c.plot([], label=r'$\bar{B}^0$', linewidth=4.5, linestyle='dashed', c='b')
876  #
877  # p1c, = ax2c.plot([], label=r'$B^0$', linewidth=4, linestyle='solid', alpha=0.9, c='r')
878  #
879  # ax2c.set_ylabel(r'${\rm Number\ \ of\ \ Events\ /\ (\, 0.02\, )}$', fontsize=35)
880  # ax2c.set_xlabel(r'$q_{\rm MC}\cdot(q\cdot r)_{\rm ' + label + '}$', fontsize=42)
881  # #plt.xticks([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1],
882  # [r'$-1$', r'', r'$-0.5$', r'', r'$0$', r'', r'$0.5$', r'', r'$1$'], rotation=0, size=40)
883  # ax2c.tick_params(axis='y', labelsize=35)
884  # ax2c.legend([ p2c, p1c], [r'$\bar{B}^0\, (q_{\rm MC} = -1)$', r'$B^0\, (q_{\rm MC} = +1)$'], prop={'size': 35}, loc=2)
885  # ax2c.grid(True)
886  # plt.xticks([-1, -0.875, -0.75, -0.625, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 0.625, 0.75, 0.875, 1],
887  # [r'${-1.0}$', r'', r'', r'', r'${-0.5}$', r'', r'', r'$0$', r'',
888  # r'', r'$0.5$', r'', r'', r'', r'$1.0$'], rotation=0, size=35)
889  #
890  # ax2c.set_ylim(0, YmaxForQrPlot)
891  # plt.savefig(workingDirectory + '/' + method + '_Folded_QR_B0bar.pdf')
892  # fig2c.clear()
893 
894  # QR Plot folded for each MC flavor with Asymmetry -----------------------------
895 
896  fig2e = plt.figure(9, figsize=(11, 11))
897  ax2e = plt.axes([0.21, 0.37, 0.76, 0.60])
898 
899  nB0bar, nbins, patches = ax2e.hist(data[np.where(data[:, 0] == -1)][:, 1] * (-1), bins=bins, histtype='step',
900  edgecolor='b', linewidth=4.5, linestyle='dashed', label=r'$\bar{B}^0$')
901 
902  nB0, nbins, patches = ax2e.hist(data[np.where(data[:, 0] == 1)][:, 1], bins=bins, histtype='step', edgecolor='r',
903  linewidth=4, alpha=0.7, label=r'$B^0$') # histtype='step', hatch='\\'
904 
905  p2c, = ax2e.plot([], label=r'$\bar{B}^0$', linewidth=4, linestyle='dashed', c='b')
906 
907  p1c, = ax2e.plot([], label=r'$B^0$', linewidth=4.5, linestyle='solid', alpha=0.9, c='r')
908 
909  ax2e.set_ylabel(r'${\rm Events\ /\ (\, 0.02\, )}$', fontsize=35)
910 
911  # plt.xticks([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1], [r'$-1$', r'',
912  # r'$-0.5$', r'', r'$0$', r'', r'$0.5$', r'', r'$1$'], rotation=0,
913  # size=40)
914  ax2e.tick_params(axis='y', labelsize=35)
915  ax2e.legend([p2c, p1c], [r'$\bar{B}^0\, (q_{\rm MC} = -1)$', r'$B^0\, (q_{\rm MC} = +1)$'], prop={'size': 35}, loc=2)
916  ax2e.grid(True)
917  plt.xticks([-1, -0.875, -0.75, -0.625, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 0.625, 0.75, 0.875, 1],
918  [r'', r'', r'', r'', r'', r'', r'', r'', r'',
919  r'', r'', r'', r'', r'', r''], rotation=0, size=35)
920 
921  plt.yticks(yLocs, yLabels)
922  ax2e.set_ylim(0, YmaxForQrPlot)
923 
924  # Asymm --
925 
926  np.seterr(divide='ignore', invalid='ignore')
927  asymmetryArray = (nB0 - nB0bar) / (nB0 + nB0bar)
928  asymmetryArrayUncertainties = np.zeros(asymmetryArray.shape)
929 
930  for i in range(0, nB0.shape[0]):
931  denominator = nB0[i] + nB0bar[i]
932  if denominator != 0:
933  asymmetryArrayUncertainties[i] = 2 * math.sqrt(
934  (nB0[i] * math.sqrt(nB0[i])) ** 2 +
935  (nB0bar[i] * math.sqrt(nB0bar[i]))**2) / (denominator)**2
936  # print(asymmetryArray)
937 
938  nBinCenters = 0.5 * (nbins[1:] + nbins[:-1])
939 
940  ax2eA = plt.axes([0.21, 0.15, 0.76, 0.2])
941 
942  # ax2eA.plot(nBinCenters, asymmetryArray, color='k',
943  # linewidth=4, linestyle='dashdot', alpha=0.9, drawstyle='steps')
944 
945  ax2eA.errorbar(nBinCenters, asymmetryArray, xerr=float(0.01),
946  yerr=asymmetryArrayUncertainties, elinewidth=2, mew=2, ecolor='k',
947  fmt='o', mfc='k', mec='k', markersize=3, label=r'${\rm Data}$')
948 
949  ax2eA.set_ylabel(r'$\frac{N^{B^0}\; -\;\, N^{\overline{B}^0}}{N^{B^0}\; +\;\, N^{\overline{B}^0}}$', fontsize=44, labelpad=20)
950  ax2eA.yaxis.labelpad = 0
951  plt.yticks([-0.4, -0.2, 0, 0.2, 0.4],
952  [r'', r'$-0.2$', r'$0.0$', r'$0.2$', r''], rotation=0, size=28)
953 
954  ax2eA.set_ylim(-0.4, 0.4)
955  ax2eA.set_xlim(ax2e.get_xlim())
956  # ax2.tick_params(axis='y', labelsize=30)
957  ax2eA.xaxis.grid(True) # linestyle='--'
958  ax2eA.yaxis.grid(True)
959  # plt.axhline(y= 4, xmin=-1.005, xmax=1.005, linewidth=1, color = 'k', linestyle = '-')
960  # plt.axhline(y= 0.75, linewidth=2, color = 'tab:gray', linestyle = '-')
961  # plt.axhline(y= 2, xmin=limXmin, xmax=limXmax, linewidth=1, color = 'k', linestyle = '-')
962  # plt.axhline(y= 0.5, linewidth=2, color = 'tab:gray', linestyle = '-')
963  # plt.axhline(y= 0.25, linewidth=2, color = 'tab:gray', linestyle = '-')
964  # plt.axhline(y= 0, linewidth=2, color = 'tab:gray', linestyle = '-')
965  # plt.axhline(y= -0.25, linewidth=2, color = 'tab:gray', linestyle = '-')
966  # plt.axhline(y= -0.5, linewidth=2, color = 'tab:gray', linestyle = '-')
967 
968  plt.xticks([-1, -0.875, -0.75, -0.625, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 0.625, 0.75, 0.875, 1],
969  [r'${-1.0}$', r'', r'', r'', r'${-0.5}$', r'', r'', r'$0$', r'',
970  r'', r'$0.5$', r'', r'', r'', r'$1.0$'], rotation=0, size=35)
971 
972  ax2eA.set_xlabel(r'$q_{\rm MC}\cdot(q\cdot r)_{\rm ' + label + '}$', fontsize=42)
973  plt.savefig(workingDirectory + '/' + method + '_Folded_QR_B0barWithAsymm.pdf')
974  fig2e.clear()
975 
976  # R Plot for B0 and B0bar ----------------------------------------------------------
977 
978  fig2d = plt.figure(11, figsize=(11, 8))
979  ax2d = plt.axes([0.21, 0.15, 0.76, 0.8])
980 
981  # ax2d.hist(np.absolute(data[np.where(abs(data[:, 0]) == 1)][:, 1]), bins=bins, histtype='step', edgecolor='k',
982  # linewidth=4, linestyle='dashdot', alpha=0.9, label=r'${\rm Both}$')
983 
984  ax2d.hist(np.absolute(data[np.where(data[:, 0] == -1)][:, 1]), bins=bins, histtype='step',
985  edgecolor='b', linewidth=4.5, linestyle='dashed', label=r'$\bar{B}^0$') # hatch='.', 'dotted'
986 
987  ax2d.hist(np.absolute(data[np.where(data[:, 0] == 1)][:, 1]), bins=bins, histtype='step', edgecolor='r',
988  linewidth=4, alpha=0.9, label=r'$B^0$') # histtype='step', hatch='\\'
989 
990  # p3d, = ax2.plot([], label=r'${\rm Both}$', linewidth=4, linestyle='dashdot', alpha=0.9, c='k')
991 
992  p2d, = ax2d.plot([], label=r'$\bar{B}^0$', linewidth=4.5, linestyle='dashed', c='b')
993 
994  p1d, = ax2d.plot([], label=r'$B^0$', linewidth=4, linestyle='solid', alpha=0.9, c='r')
995 
996  ax2d.set_ylabel(r'${\rm Number\ \ of\ \ Events\ /\ (\, 0.02\, )}$', fontsize=35)
997  ax2d.set_xlabel(r'$r_{\rm ' + label + '}$', fontsize=42)
998  # plt.xticks([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1], [r'$-1$', r'',
999  # r'$-0.5$', r'', r'$0$', r'', r'$0.5$', r'', r'$1$'], rotation=0,
1000  # size=40)
1001  ax2d.tick_params(axis='y', labelsize=35)
1002  ax2d.legend([p2d, p1d], [r'$\bar{B}^0$', r'$B^0$'], prop={'size': 35}, loc=2)
1003  ax2d.grid(True)
1004  plt.xticks([0, 0.1, 0.25, 0.5, 0.625, 0.75, 0.875, 1],
1005  [r'$0$', r'$0.1$', r'$0.25$', r'$0.5$', r'', r'$0.75$', r'', r'$1.0$'], rotation=0, size=35)
1006 
1007  ax2d.set_xlim(0, 1)
1008  ax2d.set_yscale('log', nonposy='clip')
1009  # ax2d.set_ylim(0, YmaxForQrPlot)
1010  plt.savefig(workingDirectory + '/' + method + '_Folded_R_B0bar.pdf')
1011  fig2d.clear()
1012 
1013  # --- Wrong tag fraction plot ----------------------------- ##
1014 
1015  fig5 = plt.figure(12, figsize=(11, 15))
1016  ax5 = plt.axes([0.21, 0.53, 0.76, 0.45])
1017 
1018  rBinCenters = 0.5 * (rbins[1:] + rbins[:-1])
1019  rBinWidths = 0.5 * (rbins[1:] - rbins[:-1])
1020 
1021  eb = ax5.errorbar(rBinCenters, 50 * (wvalueB0bar + wvalueB0), xerr=rBinWidths, yerr=50 * wvalueDiffUncertainty,
1022  label=r'${\rm Average}$', elinewidth=4.5, ecolor='k', fmt='o', mfc='k', mec='k',
1023  mew=0, markersize=0)
1024  eb[-1][0].set_linestyle('dashdot')
1025  eb[-1][0].set_alpha(0.9)
1026 
1027  ebB0bar = ax5.errorbar(rBinCenters, 100 * wvalueB0, xerr=rBinWidths, yerr=50 * wvalueDiffUncertainty,
1028  label=r'$B^0$', elinewidth=4, ecolor='r', fmt='o', mfc='k', mec='k',
1029  mew=0, markersize=0)
1030  ebB0bar[-1][0].set_alpha(0.9)
1031 
1032  ebB0bar = ax5.errorbar(rBinCenters, 100 * wvalueB0bar, xerr=rBinWidths, yerr=50 * wvalueDiffUncertainty,
1033  label=r'$\bar{B}^0$', elinewidth=4.5, ecolor='b', fmt='o', mfc='k', mec='k',
1034  mew=0, markersize=0)
1035  ebB0bar[-1][0].set_linestyle('--')
1036 
1037  ax5.set_ylabel(r'$w\ [\%]$', fontsize=45)
1038  ax5.yaxis.labelpad = 36
1039  ax5.set_ylim(0, 55)
1040  ax5.set_xlim(0, 1)
1041  plt.yticks([0, 10, 20, 30, 40, 50], [r'$0$',
1042  r'$10$', r'$20$', r'$30$', r'$40$', r'$50$'], rotation=0, size=44)
1043  plt.xticks([0, 0.1, 0.25, 0.5, 0.625, 0.75, 0.875, 1], [r'', r'',
1044  r'', r'', r'', r'', r'', r''], rotation=0, size=44)
1045 
1046  p3d, = ax5.plot([], label=r'${\rm Average}$', linewidth=4, linestyle='dashdot', alpha=0.9, c='k')
1047 
1048  p2d, = ax5.plot([], label=r'$\bar{B}^0$', linewidth=4.5, linestyle='dashed', c='b')
1049 
1050  p1d, = ax5.plot([], label=r'$B^0$', linewidth=4, linestyle='solid', alpha=0.9, c='r')
1051 
1052  ax5.legend([p3d, p2d, p1d], [r'${\rm Average}$', r'$\bar{B}^0$', r'$B^0$'], prop={'size': 38}, loc=1)
1053  ax5.grid(True)
1054 
1055  # Asymm plot ----
1056 
1057  ax5A = plt.axes([0.21, 0.08, 0.76, 0.43])
1058 
1059  ax5A.errorbar(rBinCenters, 100 * wvalueDiff, xerr=rBinWidths,
1060  yerr=100 * wvalueDiffUncertainty, elinewidth=3, mew=2, ecolor='k',
1061  fmt='o', mfc='k', mec='k', markersize=3, label=r'${\rm Data}$')
1062 
1063  ax5A.set_ylim(-14, 14)
1064  ax5A.set_xlim(0, 1)
1065  plt.yticks([-14, -12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12, 14],
1066  [r'', r'$-12$', r'', r'$-8$', r'', r'$-4$', r'', r'$0$',
1067  r'', r'$4$', r'', r'$8$', r'', r'$12$', r''], rotation=0, size=44)
1068 
1069  plt.xticks([0, 0.1, 0.25, 0.5, 0.625, 0.75, 0.875, 1],
1070  [r'$0$', r'$0.1$', r'$0.25$', r'$0.5$', r'', r'$0.75$', r'', r'$1.0$'], rotation=0, size=42)
1071 
1072  ax5A.xaxis.grid(True) # linestyle='--'
1073  ax5A.yaxis.grid(True)
1074  ax5A.set_ylabel(r'$\Delta w\ [\%]$', fontsize=45)
1075  ax5A.set_xlabel(r'$r_{\rm ' + label + '}$', fontsize=45)
1076 
1077  plt.savefig(workingDirectory + '/' + method + '_WrongTagFraction.pdf')
1078  fig5.clear()
1079 
1080  # --- Efficiency and Mu plot ----------------------------- ##
1081 
1082  fig6 = plt.figure(13, figsize=(11, 15))
1083  ax6 = plt.axes([0.21, 0.53, 0.76, 0.45])
1084 
1085  rBinCenters = 0.5 * (rbins[1:] + rbins[:-1])
1086  rBinWidths = 0.5 * (rbins[1:] - rbins[:-1])
1087 
1088  eb = ax6.errorbar(rBinCenters, 50 * (event_fractionB0bar + event_fractionB0), xerr=rBinWidths, yerr=50 * muParamUncertainty,
1089  label=r'${\rm Average}$', elinewidth=4.5, ecolor='k', fmt='o', mfc='k', mec='k',
1090  mew=0, markersize=0)
1091  eb[-1][0].set_linestyle('dashdot')
1092  eb[-1][0].set_alpha(0.9)
1093 
1094  ebB0bar = ax6.errorbar(rBinCenters, 100 * event_fractionB0bar, xerr=rBinWidths, yerr=50 * muParamUncertainty,
1095  label=r'$B^0$', elinewidth=4, ecolor='r', fmt='o', mfc='k', mec='k',
1096  mew=0, markersize=0)
1097  ebB0bar[-1][0].set_alpha(0.9)
1098 
1099  ebB0bar = ax6.errorbar(rBinCenters, 100 * event_fractionB0, xerr=rBinWidths, yerr=50 * muParamUncertainty,
1100  label=r'$\bar{B}^0$', elinewidth=4.5, ecolor='b', fmt='o', mfc='k', mec='k',
1101  mew=0, markersize=0)
1102  ebB0bar[-1][0].set_linestyle('--')
1103 
1104  ax6.set_ylabel(r'$\varepsilon\ [\%]$', fontsize=45)
1105  ax6.yaxis.labelpad = 36
1106  ax6.set_ylim(0, 30)
1107  ax6.set_xlim(0, 1)
1108  plt.yticks([0, 5, 10, 15, 20, 25, 30], [r'$0$', r'',
1109  r'$10$', r'', r'$20$', r'', r'$30$'], rotation=0, size=45)
1110  plt.xticks([0, 0.1, 0.25, 0.5, 0.625, 0.75, 0.875, 1], [r'', r'',
1111  r'', r'', r'', r'', r'', r''], rotation=0, size=44)
1112 
1113  p3d, = ax6.plot([], label=r'${\rm Average}$', linewidth=4, linestyle='dashdot', alpha=0.9, c='k')
1114 
1115  p2d, = ax6.plot([], label=r'$\bar{B}^0$', linewidth=4.5, linestyle='dashed', c='b')
1116 
1117  p1d, = ax6.plot([], label=r'$B^0$', linewidth=4, linestyle='solid', alpha=0.9, c='r')
1118 
1119  ax6.legend([p3d, p2d, p1d], [r'${\rm Average}$', r'$\bar{B}^0$', r'$B^0$'], prop={'size': 38}, loc=1)
1120  ax6.grid(True)
1121 
1122  # Asymm plot ----
1123 
1124  ax6A = plt.axes([0.21, 0.08, 0.76, 0.43])
1125 
1126  ax6A.errorbar(rBinCenters, 100 * muParam, xerr=rBinWidths,
1127  yerr=100 * muParamUncertainty, elinewidth=3, mew=2, ecolor='k',
1128  fmt='o', mfc='k', mec='k', markersize=3, label=r'${\rm Data}$')
1129 
1130  ax6A.set_ylim(-4, 4)
1131  ax6A.set_xlim(0, 1)
1132  plt.yticks([-3.5, -3, -2.5, -2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5],
1133  [r'', r'$-3$', r'', r'$-2$', r'', r'$-1$', r'', r'$0$', r'', r'$1$', r'', r'$2$', r'', r'$3$', r''],
1134  rotation=0, size=45)
1135 
1136  plt.xticks([0, 0.1, 0.25, 0.5, 0.625, 0.75, 0.875, 1],
1137  [r'$0$', r'$0.1$', r'$0.25$', r'$0.5$', r'', r'$0.75$', r'', r'$1.0$'], rotation=0, size=42)
1138 
1139  ax6A.xaxis.grid(True) # linestyle='--'
1140  ax6A.yaxis.grid(True)
1141  ax6A.set_ylabel(r'$\mu\ [\%]$', fontsize=45)
1142  ax6A.set_xlabel(r'$r_{\rm ' + label + '}$', fontsize=45)
1143 
1144  ax6A.yaxis.labelpad = 25
1145 
1146  plt.savefig(workingDirectory + '/' + method + '_EfficiencyAndMu.pdf')
1147  fig6.clear()
1148 
1149  # ---- TrueCategories ------------------------------------------------
1150 
1151  fig7 = plt.figure(14, figsize=(11, 8))
1152  ax7 = plt.axes([0.21, 0.15, 0.76, 0.8])
1153  weights = np.ones_like(numberOfTrueCatagories) * 100 / float(len(numberOfTrueCatagories))
1154  hi5, w, _ = ax4.hist(numberOfTrueCatagories, weights=weights, bins=list(
1155  range(0, 14, 1)), facecolor='r', histtype='stepfilled', edgecolor='k')
1156  # print(hi5)
1157  ax7.bar(np.array(range(0, 13, 1)) + 0.1, hi5, 0.8, color='r', ecolor='k')
1158  ax7.set_ylabel(r'${\rm Percentage\ of\ Events\ /\ (\, 0.02\, )}$', fontsize=25)
1159  ax7.set_xlabel(r'${\rm True\ Categories}$', fontsize=25)
1160  plt.xticks([0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5],
1161  [r'$0$', r'$1$', r'$2$', r'$3$', r'$4$', r'$5$', r'$6$', r'$7$', r'$8$'],
1162  rotation=0, horizontalalignment='center', size=25)
1163  ax7.tick_params(axis='y', labelsize=25)
1164  ax7.set_xlim(0, 8)
1165  ax7.set_ylim(0, 40)
1166  # ax7.grid(True)
1167  plt.savefig(workingDirectory + '/' + method + 'TrueCategories.pdf')
1168  fig7.clear()
1169 
1170  fig = plt.figure(figsize=(8, 8))
1171  ax = fig.add_subplot(111, projection='3d')
1172 
1173  # -----bins
1174  bins2 = list(range(-51, 51, 1))
1175  for i in range(0, len(bins2)):
1176  bins2[i] = float(bins2[i]) / 50
1177  # ------
1178  bins2 = np.array(bins2)
1179  # bins2 = np.arange(-1, 0.5, 0.02)
1180 
1181  def cc(arg):
1182  return colorConverter.to_rgba(arg, alpha=0.5)
1183 
1184  vertsa = []
1185  vertsb = []
1186  colorsa = []
1187  colorsb = []
1188  zs = list(range(0, 7, 1))
1189 
1190  Zmax = 0
1191 
1192  iEfficiencies = []
1193 
1194  for z in range(0, 7, 1):
1195  dataTrueCats = dataAll[(np.where(dataAll[:, 2] == z))][:, 0:2]
1196  weightsa = np.ones_like(dataTrueCats[np.where(dataTrueCats[:, 0] == -1)][:, 1]) * \
1197  (hi5[z] * dataAll.shape[0] / 100) / float(len(dataTrueCats[:, 1]))
1198  weightsb = np.ones_like(dataTrueCats[np.where(dataTrueCats[:, 0] == 1)][:, 1]) * \
1199  (hi5[z] * dataAll.shape[0] / 100) / float(len(dataTrueCats[:, 1]))
1200  ya, wa, _ = ax2.hist(dataTrueCats[np.where(dataTrueCats[:, 0] == -1)][:, 1], weights=weightsa, bins=bins, histtype='step',
1201  edgecolor='b', linewidth=3.5, linestyle='dotted', label=r'$\bar{B}^0$')
1202  yb, wb, _ = ax2.hist(dataTrueCats[np.where(dataTrueCats[:, 0] == 1)][:, 1], weights=weightsb, bins=bins, histtype='step',
1203  edgecolor='r', linewidth=2, label=r'$B^0$')
1204  ya = np.insert(ya, ya.shape[0], 0)
1205  ya = np.insert(ya, 0, 0)
1206  yb = np.insert(yb, yb.shape[0], 0)
1207  yb = np.insert(yb, 0, 0)
1208  if ya.max() > Zmax:
1209  Zmax = ya.max()
1210  if yb.max() > Zmax:
1211  Zmax = yb.max()
1212  colorsa.append(cc('b'))
1213  colorsb.append(cc('r'))
1214  vertsa.append(list(zip(bins2, ya)))
1215  vertsb.append(list(zip(bins2, yb)))
1216 
1217  iEfficiency, iEfficiencyDiff, iEfficiencyB0, iEfficiencyB0bar = efficiencyCalculator(dataTrueCats, total_notTagged, True)
1218 
1219  iEfficiencies.append([iEfficiencyB0, iEfficiencyB0bar])
1220 
1221  polya = PolyCollection(vertsa, facecolors=colorsa)
1222  polyb = PolyCollection(vertsb, facecolors=colorsb)
1223  polya.set_alpha(0.5)
1224  polyb.set_alpha(0.5)
1225  ax.add_collection3d(polya, zs=zs, zdir='y')
1226  ax.add_collection3d(polyb, zs=zs, zdir='y')
1227 
1228  # for z in zs:
1229  # ya, wa , _ = ax2b.hist(dataNoTrueCategory[np.where(dataNoTrueCategory[:, 0] == -1)][:, 1], bins=bins, histtype='step',
1230  # edgecolor='b', linewidth=2, label=r'$\bar{B}^0$')
1231  # ax.bar(bins2, bins2, zs=z, zdir = 'y', color='b', ecolor = 'b', alpha=0.75)
1232  # print(bins2)
1233 
1234  ax.set_xlim3d(-1, 1)
1235  ax.set_ylim3d(-0.5, 7)
1236  ax.set_yticks([0, 1, 2, 3, 4, 5, 6])
1237  ax.xaxis.labelpad = 18
1238  ax.yaxis.labelpad = 16
1239  ax.zaxis.labelpad = 22
1240  # ax.tick_params(axis = 'y', right = 'on')
1241  ax.set_yticklabels(sorted({r'$0$', r'$1$', r'$2$', r'$3$', r'$4$', r'$5$', r'$6$'}),
1242  verticalalignment='baseline',
1243  horizontalalignment='left')
1244 
1245  ax.set_xlabel(r'$(q\cdot r)_{\rm ' + label + ' }$')
1246  ax.set_ylabel(r'${\rm True\ Categories}$')
1247  ax.set_zlabel(r'${\rm Number\ \ of\ \ Events\ /\ (\, 0.02\, )}$', rotation=-90)
1248  Zmax = Zmax + Zmax / 5
1249  ax.set_zlim3d(0, Zmax)
1250  p1 = plt.Rectangle((0, 0), 1, 1, fc="b", alpha=0.5)
1251  p2 = plt.Rectangle((0, 0), 1, 1, fc="r", alpha=0.5)
1252  ax.legend([p1, p2], [r'$\bar{B}^0$', r'$B^0$'], prop={'size': 20}, loc=(0.523, 0.687)) # (0.639, 0.71
1253  plt.savefig(workingDirectory + '/3DTrueCategorieQR' + method + '.pdf')
1254  fig.clear()
1255 
1256  iEfficiencies = np.array(iEfficiencies)
1257 
1258  fig = plt.figure(15, figsize=(12, 8))
1259  ax = plt.axes([0.18, 0.15, 0.8, 0.8])
1260  # print(hi5)
1261  width = 0.4
1262  p1 = ax.bar(np.array(range(0, 7, 1)) + width / 4, iEfficiencies[:, 1], width, color='b', ecolor='k', label=r'$\bar{B}^0$')
1263  p2 = ax.bar(np.array(range(0, 7, 1)) + width + width / 4, iEfficiencies[:, 0], width, color='r', ecolor='k', label=r'$B^0$')
1264  ax.set_ylabel(r'$ \varepsilon_{\rm eff}/\ \% $', fontsize=35)
1265  ax.set_xlabel(r'${\rm True\ Categories}$', fontsize=35)
1266  plt.xticks([0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5],
1267  [r'$0$', r'$1$', r'$2$', r'$3$', r'$4$', r'$5$', r'$6$', r'$7$', r'$8$'],
1268  rotation=0, horizontalalignment='center', size=35)
1269  ax.tick_params(axis='y', labelsize=35)
1270  ax.set_xlim(0, 8)
1271  ax.set_ylim(0, 100)
1272  ax.legend([r'$\bar{B}^0$', r'$B^0$'], prop={'size': 35})
1273 
1274  plt.savefig(workingDirectory + '/EfficienciesTrueCats' + method + '.pdf')
1275  fig.clear()
1276 
1277 
1278 NoTargetEfficiencies = []
1279 categoryLabel = []
1280 
1281 categoryLabelsDict = {'Electron': r'${\rm Electron}$',
1282  'IntermediateElectron': r'${\rm Int. Electron}$',
1283  'Muon': r'${\rm Muon}$',
1284  'IntermediateMuon': r'${\rm Int. Muon}$',
1285  'KinLepton': r'${\rm KinLepton}$',
1286  'IntermediateKinLepton': r'${\rm Int. KinLepton}$',
1287  'Kaon': r'${\rm Kaon}$',
1288  'SlowPion': r'${\rm Slow Pion}$',
1289  'FastHadron': r'${\rm Fast Hadron}$',
1290  'Lambda': r'${\rm Lambda}$',
1291  'FSC': r'${\rm FSC}$',
1292  'MaximumPstar': r'${\rm MaximumPstar}$',
1293  'KaonPion': r'${\rm Kaon-Pion}$'}
1294 
1295 # eventLevelParticleLists = [ ('Lambda0:inRoe', 'Lambda',
1296 # 'weightedQpOf(Lambda0:inRoe, isRightCategory(Lambda),
1297 # isRightCategory(Lambda))') ]
1298 
1299 # print(eventLevelParticleLists)
1300 
1301 # eventLevelParticleLists = [
1302 # ('e+:inRoe', 'Electron',
1303 # 'QpOf(e+:inRoe, isRightCategory(Electron), isRightCategory(Electron))'),
1304 # ('e+:inRoe', 'IntermediateElectron',
1305 # 'QpOf(e+:inRoe, isRightCategory(IntermediateElectron), isRightCategory(IntermediateElectron))'),
1306 # ('mu+:inRoe', 'Muon',
1307 # 'QpOf(mu+:inRoe, isRightCategory(Muon), isRightCategory(Muon))'),
1308 # ('mu+:inRoe', 'IntermediateMuon',
1309 # 'QpOf(mu+:inRoe, isRightCategory(IntermediateMuon), isRightCategory(IntermediateMuon))'),
1310 # ('mu+:inRoe', 'KinLepton',
1311 # 'QpOf(mu+:inRoe, isRightCategory(KinLepton), isRightCategory(KinLepton))'),
1312 # ('mu+:inRoe', 'IntermediateKinLepton',
1313 # 'QpOf(mu+:inRoe, isRightCategory(IntermediateKinLepton), isRightCategory(IntermediateKinLepton))'),
1314 # ('K+:inRoe', 'Kaon',
1315 # 'weightedQpOf(K+:inRoe, isRightCategory(Kaon), isRightCategory(Kaon))'),
1316 # ('K+:inRoe', 'KaonPion',
1317 # 'QpOf(K+:inRoe, isRightCategory(KaonPion), isRightCategory(Kaon))'),
1318 # ('pi+:inRoe', 'SlowPion', 'QpOf(pi+:inRoe, isRightCategory(SlowPion), isRightCategory(SlowPion))'),
1319 # ('pi+:inRoe', 'FSC', 'QpOf(pi+:inRoe, isRightCategory(FSC), isRightCategory(SlowPion))'),
1320 # ('pi+:inRoe', 'MaximumPstar',
1321 # 'QpOf(pi+:inRoe, isRightCategory(MaximumPstar), isRightCategory(MaximumPstar))'),
1322 # ('pi+:inRoe', 'FastHadron',
1323 # 'QpOf(pi+:inRoe, isRightCategory(FastHadron), isRightCategory(FastHadron))'),
1324 # ('Lambda0:inRoe', 'Lambda',
1325 # 'weightedQpOf(Lambda0:inRoe, isRightCategory(Lambda), isRightCategory(Lambda))')]
1326 
1327 
1328 qrMC = ROOT.RooRealVar('qrMC', 'qrMC', 0, -1.0, 1.0)
1329 argSet = ROOT.RooArgSet(qrMC)
1330 
1331 for iCategory in usedCategories:
1332  #
1333  with Quiet(ROOT.kError):
1334  exec(
1335  f"hasTrueTarget + {iCategory} = " +
1336  "ROOT.RooRealVar('hasTrueTarget' + iCategory, 'hasTrueTarget' + iCategory, 0, -2.1, 1.1)"
1337  )
1338  exec(
1339  f"qp + {iCategory} = ROOT.RooRealVar('qp' + iCategory, 'hasTrueTarget' + iCategory, 0, -2.1, 1.1)"
1340  )
1341  # exec("%s = %s" % ("isTrueCategory" + iCategory,
1342  # "ROOT.RooRealVar('isRightCategory' + iCategory, 'isRightCategory'
1343  # + iCategory, 0, -2.0, 1.1)"))
1344  exec(f"argSet.add(hasTrueTarget{iCategory})")
1345  exec(f"argSet.add(qp{iCategory})")
1346 rooDataSet = ROOT.RooDataSet("data", "data", tree, argSet, "")
1347 
1348 for category in usedCategories:
1349 
1350  print(category)
1351  # qrCombined=ROOT.RooRealVar('qp' + category, 'qp' + category, 0, -1.0, 1.1)
1352 
1353  data = []
1354 
1355  dataTruth = []
1356 
1357  dataNoTarget = []
1358 
1359  for i in range(rooDataSet.numEntries()):
1360  row = rooDataSet.get(i)
1361  tqp = row.getRealValue('qp' + category, 0, ROOT.kTRUE)
1362  thasTrueTarget = row.getRealValue('hasTrueTarget' + category, 0, ROOT.kTRUE)
1363  tqrMC = row.getRealValue('qrMC', 0, ROOT.kTRUE)
1364 
1365  tNumberOfTrueCategories = 0
1366  noMCAssociated = False
1367 
1368  for iCategory in usedCategories:
1369 
1370  if iCategory == "MaximumPstar":
1371  continue
1372  ihasTrueTarget = row.getRealValue('hasTrueTarget' + iCategory, 0, ROOT.kTRUE)
1373  # print(ihasTrueTarget)
1374 
1375  if ihasTrueTarget > 0:
1376  tNumberOfTrueCategories += 1
1377  if ihasTrueTarget < 0:
1378  noMCAssociated = True
1379  break
1380 
1381  if tqp > 1:
1382  tqp = 1.0
1383 
1384  # if thasTrueTarget < 1:
1385  # print(thasTrueTarget)
1386 
1387  if abs(tqrMC) == 1: # !=0:
1388  data.append([tqrMC, tqp])
1389  if thasTrueTarget > 0:
1390  dataTruth.append([tqrMC, tqp])
1391  if tNumberOfTrueCategories < 1 and not noMCAssociated:
1392  # thasTrueTarget < 1:#if thasTrueTarget == 0 and not noMCAssociated:
1393  if category == "Lambda" and tqp == 0:
1394  continue
1395  # print(thasTrueTarget)
1396  dataNoTarget.append([tqrMC, tqp])
1397 
1398  data = np.array(data)
1399  dataTruth = np.array(dataTruth)
1400  dataNoTarget = np.array(dataNoTarget)
1401 
1402  # print(dataNoTarget)
1403 
1404  noTargetEfficiency, noTargetEfficiencyDiff, noTargetEfficiencyB0, noTargetEfficiencyB0bar = categoriesEfficiencyCalculator(data)
1405  print("Efficiencies for B0, B0bar = ", f'{noTargetEfficiencyB0: 6.2f}',
1406  " %, ", f'{noTargetEfficiencyB0bar: 6.2f}', " %")
1407 
1408  noTargetEfficiency, noTargetEfficiencyDiff, noTargetEfficiencyB0, noTargetEfficiencyB0bar = categoriesEfficiencyCalculator(
1409  dataNoTarget)
1410  print("Efficiencies for B0, B0bar If Not Target = ", f'{noTargetEfficiencyB0: 6.2f}', " %, ",
1411  f'{noTargetEfficiencyB0bar: 6.2f}', " %")
1412 
1413  NoTargetEfficiencies.append([noTargetEfficiencyB0, noTargetEfficiencyB0bar])
1414  categoryLabel.append(categoryLabelsDict[category])
1415 
1416  trueTargetEfficiency, trueTargetEfficiencyDiff, trueTargetEfficiencyB0, \
1417  trueTargetEfficiencyB0bar = categoriesEfficiencyCalculator(dataTruth)
1418  print("Efficiencies for B0, B0bar If True Target= ", f'{trueTargetEfficiencyB0: 6.2f}', " %, ",
1419  f'{trueTargetEfficiencyB0bar: 6.2f}', " %")
1420 
1421  # catLabel = r"${\rm " + category + "}$"
1422  # categoryLabel.append(catLabel)
1423 
1424  # -----bins
1425  bins = list(range(-50, 51, 1))
1426  for i in range(0, len(bins)):
1427  bins[i] = float(bins[i]) / 50
1428  # ------
1429 
1430  title = ''
1431  location = 1
1432  if category != 'Lambda' and category != 'MaximumPstar' and category != 'Kaon':
1433  title = r'$q_{\rm cand}\cdot y_{\rm ' + category + '}$'
1434  elif category == 'MaximumPstar':
1435  title = r'$q_{\rm cand}\cdot y_{{\rm Maximum}\, p^*}$'
1436  location = 9
1437  elif category == 'Kaon':
1438  location = 8
1439  title = r'$(q_{\rm cand}\cdot y_{\rm Kaon})_{\rm eff}$'
1440  elif category == 'Lambda':
1441  title = r'$(q_{\rm \Lambda}\cdot y_{\rm Lambda})_{\rm eff}$'
1442 
1443  if category == 'IntermediateKinLepton':
1444  location = 8
1445  title = r'$q_{\rm cand}\cdot y_{\rm Int.\, Kin.\, Lepton}$'
1446  elif category == 'KinLepton':
1447  title = r'$q_{\rm cand}\cdot y_{\rm Kin.\, Lepton}$'
1448  elif category == 'IntermediateMuon':
1449  title = r'$q_{\rm cand}\cdot y_{\rm Int.\, Muon}$'
1450  elif category == 'IntermediateElectron':
1451  title = r'$q_{\rm cand}\cdot y_{\rm Int.\, Electron}$'
1452  elif category == 'KaonPion':
1453  title = r'$q_{\rm cand}\cdot y_{\rm Kaon-Pion}$'
1454  elif category == 'FastHadron':
1455  location = 8
1456  title = r'$q_{\rm cand}\cdot y_{\rm Fast\, Hadron}$'
1457  elif category == 'SlowPion':
1458  title = r'$q_{\rm cand}\cdot y_{\rm Slow\, Pion}$'
1459 
1460  fig2 = plt.figure(2, figsize=(11, 8))
1461  ax2 = plt.axes([0.21, 0.16, 0.75, 0.81]) # plt.axes([0.14, 0.1, 0.8, 0.8])
1462  ax2.hist(data[np.where(data[:, 0] == -1)][:, 1], bins=bins, histtype='step',
1463  edgecolor='b', linewidth=4.5, linestyle='dashed', label=r'$\bar{B}^0$') # hatch='.',
1464  ax2.hist(
1465  data[np.where(data[:, 0] == 1)][:, 1],
1466  bins=bins,
1467  histtype='step',
1468  edgecolor='r',
1469  linewidth=4,
1470  alpha=0.9,
1471  label=r'$B^0$') # histtype='step', hatch='\\'
1472 
1473  p2, = ax2.plot([], label=r'$\bar{B}^0$', linewidth=4.5, linestyle='dashed', c='b')
1474  p1, = ax2.plot([], label=r'$B^0$', linewidth=4, linestyle='solid', alpha=0.9, c='r')
1475 
1476  ax2.set_ylabel(r'${\rm Number\ \ of\ \ Events\ /\ (\, 0.02\, )}$', fontsize=37)
1477  ax2.set_xlabel(title, fontsize=48)
1478  ax2.set_yscale('log', nonposy='clip')
1479  plt.xticks([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1],
1480  [r'$-1$', r'', r'$-0.5$', r'', r'$0$', r'', r'$0.5$', r'', r'$1$'], rotation=0, size=37)
1481  ax2.tick_params(axis='y', labelsize=37)
1482  ax2.legend([p2, p1], [r'$\bar{B}^0$', r'$B^0$'], prop={'size': 35}, loc=location)
1483  ax2.grid(True)
1484  plt.savefig(workingDirectory + '/' + 'pyPIC_' + category + '_Input_Combiner.pdf')
1485  fig2.clear()
1486 
1487  fig2b = plt.figure(2, figsize=(11, 8))
1488  ax2b = plt.axes([0.21, 0.16, 0.75, 0.81]) # plt.axes([0.14, 0.1, 0.8, 0.8])
1489  ax2b.hist(dataTruth[np.where(dataTruth[:, 0] == -1)][:, 1], bins=bins, histtype='step',
1490  edgecolor='b', linewidth=4.5, linestyle='dashed', label=r'$\bar{B}^0$') # hatch='.',
1491  ax2b.hist(
1492  dataTruth[np.where(dataTruth[:, 0] == 1)][:, 1],
1493  bins=bins,
1494  histtype='step',
1495  edgecolor='r',
1496  linewidth=4,
1497  alpha=0.9,
1498  label=r'$B^0$') # histtype='step', hatch='\\'
1499 
1500  p2, = ax2b.plot([], label=r'$\bar{B}^0$', linewidth=4.5, linestyle='dashed', c='b')
1501  p1, = ax2b.plot([], label=r'$B^0$', linewidth=4, linestyle='solid', alpha=0.9, c='r')
1502 
1503  ax2b.set_ylabel(r'${\rm Number\ \ of\ \ Events\ /\ (\, 0.02\, )}$', fontsize=37)
1504  ax2b.set_xlabel(title, fontsize=48)
1505  ax2b.set_yscale('log', nonposy='clip')
1506  plt.xticks([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1],
1507  [r'$-1$', r'', r'$-0.5$', r'', r'$0$', r'', r'$0.5$', r'', r'$1$'], rotation=0, size=37)
1508  ax2b.tick_params(axis='y', labelsize=37)
1509  ax2b.legend([p2, p1], [r'$\bar{B}^0$', r'$B^0$'], prop={'size': 35}, loc=location)
1510  ax2b.grid(True)
1511  plt.savefig(workingDirectory + '/' + 'pyPIC_' + category + '_Input_CombinerTruth.pdf')
1512  fig2b.clear()
1513 
1514  fig2c = plt.figure(2, figsize=(11, 8))
1515  ax2c = plt.axes([0.21, 0.16, 0.75, 0.81]) # plt.axes([0.14, 0.1, 0.8, 0.8])
1516  ax2c.hist(dataNoTarget[np.where(dataNoTarget[:, 0] == -1)][:, 1], bins=bins, histtype='step',
1517  edgecolor='b', linewidth=4.5, linestyle='dashed', label=r'$\bar{B}^0$') # hatch='.',
1518  ax2c.hist(
1519  dataNoTarget[np.where(dataNoTarget[:, 0] == 1)][:, 1],
1520  bins=bins,
1521  histtype='step',
1522  edgecolor='r',
1523  linewidth=4,
1524  alpha=0.9,
1525  label=r'$B^0$') # histtype='step', hatch='\\'
1526 
1527  p2, = ax2c.plot([], label=r'$\bar{B}^0$', linewidth=4.5, linestyle='dashed', c='b')
1528  p1, = ax2c.plot([], label=r'$B^0$', linewidth=4, linestyle='solid', alpha=0.9, c='r')
1529 
1530  ax2c.set_ylabel(r'${\rm Number\ \ of\ \ Events\ /\ (\, 0.02\, )}$', fontsize=37)
1531  ax2c.set_xlabel(title, fontsize=48)
1532  ax2c.set_yscale('log', nonposy='clip')
1533  plt.xticks([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1],
1534  [r'$-1$', r'', r'$-0.5$', r'', r'$0$', r'', r'$0.5$', r'', r'$1$'], rotation=0, size=37)
1535  ax2c.tick_params(axis='y', labelsize=37)
1536  ax2c.legend([p2, p1], [r'$\bar{B}^0$', r'$B^0$'], prop={'size': 35}, loc=location)
1537  ax2c.grid(True)
1538  plt.savefig(workingDirectory + '/' + 'pyPIC_' + category + '_Input_CombinerNoTarget.pdf')
1539  fig2c.clear()
1540 
1541  percentageCategory = dataTruth.shape[0] / data.shape[0] * 100
1542  percentageCategory = f'{percentageCategory:.2f}'
1543 
1544  print("Percentage of Category " + category + " is = " + percentageCategory + " %")
1545 
1546  # print("Data size = " + str(data.shape[0]))
1547 
1548 NoTargetEfficiencies = np.array(NoTargetEfficiencies)
1549 fig = plt.figure(5, figsize=(29, 18))
1550 ax = plt.axes([0.15, 0.3, 0.81, 0.67])
1551 # print(hi5)
1552 width = 0.4
1553 p1 = ax.bar(np.array(range(0, 13, 1)) + width / 4, NoTargetEfficiencies[:, 1], width, color='b', ecolor='k', label=r'$\bar{B}^0$')
1554 p2 = ax.bar(np.array(range(0, 13, 1)) + width + width / 4, NoTargetEfficiencies[:, 0], width, color='r', ecolor='k', label=r'$B^0$')
1555 ax.set_ylabel(r'$ \varepsilon_{\rm eff}/\ \% $', fontsize=30)
1556 plt.xticks([0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5],
1557  categoryLabel, rotation=50, horizontalalignment='right', size=30)
1558 plt.yticks([5, 10, 15, 20, 25, 50, 100], [r"$5$", r"", r"$15$", r"", r"$25$", r"$50$", r"$100$"], size=30)
1559 ax.set_xlim(0, 13)
1560 ax.set_ylim(0, 100)
1561 ax.legend([r'$\bar{B}^0$', r'$B^0$'], prop={'size': 25})
1562 ax.grid(True)
1563 plt.savefig(workingDirectory + '/EfficienciesNoTarget.pdf')
1564 fig.clear()
1565 
1566 print('* *')
1567 print('*******************************************************************************************************************')
1568 B2INFO('qr Output Histograms in pdf format saved at: ' + workingDirectory)