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