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