Belle II Software  release-05-02-19
pythonFlavorTaggerEfficiency.py
1 #!/usr/bin/env python3
2 # -*- coding: utf-8 -*-
3 
4 # ********************************** Flavor Tagging *******************************
5 # * This scripts calculates the effective efficiency of the flavor tagger based *
6 # * on python histograms. It serves as a crosscheck for the script *
7 # * flavorTaggerEfficiency.py, which uses root histograms, and produces nicer plots. *
8 # * This script also produces plots for different amounts of true categories *
9 # * for the combiner outputs. It produces also plots for individual categories *
10 # * checking when they are true categories and when they are not true. *
11 # * True here means that the target (or targets) of the category are found in a *
12 # * certain event. For more information check Sec. 4.5.3 in BELLE2-PTHESIS-2018-003 *
13 # *
14 # * Author: Fernando Abudinen *
15 # *
16 # ************************************************************************************
17 
18 import ROOT
19 from ROOT import Belle2
20 from basf2 import B2INFO
21 import flavorTagger as ft
22 from defaultEvaluationParameters import categories, Quiet, rbins
23 import basf2_mva
24 from array import array
25 
26 import numpy as np
27 import matplotlib as mpl
28 mpl.use('Agg')
29 mpl.rcParams.update({'font.size': 22})
30 mpl.rcParams['text.usetex'] = True
31 mpl.rcParams['text.latex.preamble'] = [r"\usepackage{amsmath}"]
32 import matplotlib.pyplot as plt
33 from matplotlib.ticker import FormatStrFormatter
34 import math
35 import sys
36 import glob
37 
38 from ROOT import PyConfig
39 PyConfig.IgnoreCommandLineOptions = True
40 PyConfig.StartGuiThread = False
41 from array import array
42 import IPython
43 from mpl_toolkits.mplot3d import Axes3D
44 from matplotlib.collections import PolyCollection
45 from matplotlib.colors import colorConverter
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)