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