Belle II Software  light-2212-foldex
genLevelAsymmsImpactParams.py
1 #!/usr/bin/env python3
2 # -*- coding: utf-8 -*-
3 
4 
11 
12 
22 
23 
24 import basf2 as b2
25 import os
26 import glob
27 import sys
28 import math
29 from matplotlib.ticker import FormatStrFormatter
30 import matplotlib.pyplot as plt
31 import ROOT
32 from ROOT import Belle2
33 
34 import numpy as np
35 import matplotlib as mpl
36 mpl.use('Agg')
37 mpl.rcParams.update({'font.size': 22})
38 mpl.rcParams['text.usetex'] = True
39 mpl.rcParams['text.latex.preamble'] = [r"\usepackage{amsmath}"]
40 
41 if len(sys.argv) != 3:
42  sys.exit("Must provide 2 arguments: [Belle or Belle2] [samplesWildCards]"
43  )
44 
45 belleOrBelle2 = sys.argv[1]
46 samplesWildCards = sys.argv[2]
47 
48 
49 workingDirectory = '.'
50 
51 if not b2.find_file('GenLevelVariablesPlots', silent=True):
52  os.mkdir(workingDirectory + '/GenLevelVariablesPlots')
53 
54 sourceFiles = glob.glob(str(samplesWildCards))
55 
56 tree = ROOT.TChain("tree")
57 for Inputfile in sourceFiles:
58  tree.AddFile(Inputfile)
59 
60 
61 ROOT.TH1.SetDefaultSumw2()
62 
63 binsPt = 50
64 
65 HistoBelle2_PtMC = ROOT.TH1F('PtMCBelle2', 'MC Transverse Momentum', binsPt, 0, 5)
66 HistoBelle2_PtFit = ROOT.TH1F('PtFitBelle2', 'TrackFitResult Transverse Momentum', binsPt, 0, 5)
67 
68 HistoBelle2_pseudoPD0MC = ROOT.TH1F('pseudoPD0MCBelle2', 'MC pseudo Momentum for D0', binsPt, 0, 5)
69 HistoBelle2_pseudoPD0Fit = ROOT.TH1F('pseudoPD0FitBelle2', 'Fit pseudo Momentum for D0', binsPt, 0, 5)
70 
71 HistoBelle2_pseudoPZ0MC = ROOT.TH1F('pseudoPZ0MCBelle2', 'MC pseudo Momentum for Z0', binsPt, 0, 5)
72 HistoBelle2_pseudoPZ0Fit = ROOT.TH1F('pseudoPZ0FitBelle2', 'Fit pseudo Momentum for Z0', binsPt, 0, 5)
73 
74 HistoBelle2_d0Error = ROOT.TH1F('d0ErrorBelle2', 'TrackFitResult Error on d0', binsPt, 0, 0.03)
75 HistoBelle2_z0Error = ROOT.TH1F('z0ErrorBelle2', 'TrackFitResult Error on z0', binsPt, 0, 0.03)
76 
77 D0 = 'abs(TrackFitResults.m_tau[0][0])'
78 Z0 = 'abs(TrackFitResults.m_tau[0][3])'
79 X0 = 'sqrt(TrackFitResults.m_tau[0][0]**2 + TrackFitResults.m_tau[0][3]**2)'
80 
81 mcX = "MCParticles.m_productionVertex_x"
82 mcY = "MCParticles.m_productionVertex_y"
83 mcZ = "MCParticles.m_productionVertex_z"
84 
85 mcD0 = "sqrt(MCParticles.m_productionVertex_x*MCParticles.m_productionVertex_x + " + \
86  " MCParticles.m_productionVertex_y*MCParticles.m_productionVertex_y)"
87 mcZ0 = "abs(MCParticles.m_productionVertex_z)"
88 mcX0 = "sqrt(MCParticles.m_productionVertex_x*MCParticles.m_productionVertex_x + " + \
89  "MCParticles.m_productionVertex_y*MCParticles.m_productionVertex_y + " + \
90  "MCParticles.m_productionVertex_z*MCParticles.m_productionVertex_z)"
91 
92 ptMC = 'sqrt(MCParticles.m_momentum_x**2 + MCParticles.m_momentum_y**2)'
93 ptFit = 'abs(0.003*1.5/TrackFitResults.m_tau[0][2])'
94 
95 pMC = 'sqrt(MCParticles.m_momentum_x**2 + MCParticles.m_momentum_y**2 + + MCParticles.m_momentum_z**2)'
96 pFit = '(abs(0.003*1.5/TrackFitResults.m_tau[0][2])*sqrt(1+TrackFitResults.m_tau[0][4]**2))'
97 chiProb = 'TrackFitResults.m_pValue'
98 
99 pseudoPD0MC = 'sqrt(' + ptMC + '**3/(' + pMC + '*(1+(MCParticles.m_mass[0]/' + pMC + ')**2)))'
100 pseudoPD0Fit = 'sqrt(' + ptFit + '**3/(' + pFit + '*(1+(MCParticles.m_mass[0]/' + pFit + ')**2)))'
101 
102 pseudoPZ0MC = 'sqrt(' + ptMC + '**5/(' + pMC + '**3*(1+(MCParticles.m_mass[0]/' + pMC + ')**2)))'
103 pseudoPZ0Fit = 'sqrt(' + ptFit + '**5/(' + pFit + '**3*(1+(MCParticles.m_mass[0]/' + pFit + ')**2)))'
104 
105 
106 dBw = 50
107 pBins = 50
108 fBins = 100
109 
110 particleConditions = {
111  'Electron': ['e', "11 ", " MCParticles.m_pdg > 0 ", " MCParticles.m_pdg < 0 "],
112  'Muon': [r'\mu', "13 ", " MCParticles.m_pdg > 0 ", " MCParticles.m_pdg < 0 "],
113  'Kaon': ['K', "321 ", " MCParticles.m_pdg < 0 ", " MCParticles.m_pdg > 0 "],
114  'Pion': [r'\pi', "211 ", " MCParticles.m_pdg < 0 ", " MCParticles.m_pdg > 0 "],
115  'Proton': ['p', "2212 ", " MCParticles.m_pdg < 0 ", " MCParticles.m_pdg > 0 "]}
116 
117 
118 Particles = ["Electron", "Muon", "Kaon", "Pion", "Proton"]
119 
120 unitImp = "mm"
121 
122 
123 def makePlotsForEachParticleKind(cutOnUpsilonFourS=""):
124  """
125  Makes plots showing the distribution of the impact parameters and vertices
126  for positively and negatively charged particles.
127  """
128 
129  unitImp = "mm"
130 
131  variablesPlotParamsDict = {'ImpactXY': [D0, dBw, 0, 0.5, r'$d_0\ [{\rm ' + unitImp + '}]$', r"{\rm " + unitImp + r"}\, "],
132  'z0': [Z0, dBw, 0, 1.0, r'$z_0\ [{\rm ' + unitImp + '}]$', r"{\rm " + unitImp + r"}\, "],
133  'distance': [X0, dBw, 0, 1.5, r'$\xi_0\ [{\rm ' + unitImp + '}]$', r"{\rm " + unitImp + r"}\, "],
134  'mcGenX': [mcX, dBw, -0.05, 0.15, r'$x_{\rm MC}\ [{\rm ' + unitImp + '}]$',
135  r"{\rm " + unitImp + r"}\, "],
136  'mcGenY': [mcY, dBw, -0.075, 0.075, r'$y_{\rm MC}\ [{\rm ' + unitImp + '}]$',
137  r"{\rm " + unitImp + r"}\, "],
138  'mcGenZ': [mcZ, dBw, -0.6, 1, r'$z_{\rm MC}\ [{\rm ' + unitImp + '}]$',
139  r"{\rm " + unitImp + r"}\, "],
140  'mcGenImpactXY': [mcD0, dBw, 0, 0.1,
141  r'$\sqrt{x_{\rm MC}^2 + y_{\rm MC}^2}\ [{\rm ' + unitImp + '}]$',
142  r"{\rm " + unitImp + r"}\, "],
143  'mcGenz0': [mcZ0, dBw, 0, 1, r'$|z_{\rm MC}|\ [{\rm ' + unitImp + '}]$',
144  r"{\rm " + unitImp + r"}\, "],
145  'mcGendistance': [mcX0, dBw, 0, 1,
146  r'$\sqrt{x_{\rm MC}^2 + y_{\rm MC}^2 + z_{\rm MC}^2}\ [{\rm ' + unitImp + '}]$',
147  r"{\rm " + unitImp + r"}\, "],
148  }
149 
150  if belleOrBelle2 == "Belle":
151  unitImp = "cm"
152 
153  variablesPlotParamsDict = {
154  'ImpactXY': [
155  D0,
156  dBw,
157  0,
158  0.1,
159  r'$d_0\ [{\rm ' + unitImp + '}]$',
160  r"{\rm " + unitImp + r"}\, "],
161  'z0': [
162  Z0,
163  dBw,
164  0,
165  1.0,
166  r'$z_0\ [{\rm ' + unitImp + '}]$',
167  r"{\rm " + unitImp + r"}\, "],
168  'distance': [
169  X0,
170  dBw,
171  0,
172  1.5,
173  r'$\xi_0\ [{\rm ' + unitImp + '}]$',
174  r"{\rm " + unitImp + r"}\, "],
175  'mcGenX': [
176  mcX,
177  dBw,
178  0,
179  0.1,
180  r'$x_{\rm MC}\ [{\rm ' + unitImp + '}]$',
181  r"{\rm " + unitImp + r"}\, "],
182  'mcGenY': [
183  mcY,
184  dBw,
185  -0.005,
186  0.02,
187  r'$y_{\rm MC}\ [{\rm ' + unitImp + '}]$',
188  r"{\rm " + unitImp + r"}\, "],
189  'mcGenZ': [
190  mcZ,
191  dBw,
192  -0.6,
193  1,
194  r'$z_{\rm MC}\ [{\rm ' + unitImp + '}]$',
195  r"{\rm " + unitImp + r"}\, "],
196  'mcGenImpactXY': [
197  mcD0,
198  dBw,
199  0,
200  0.1,
201  r'$\sqrt{x_{\rm MC}^2 + y_{\rm MC}^2}\ [{\rm ' + unitImp + '}]$',
202  r"{\rm " + unitImp + r"}\, "],
203  'mcGenz0': [
204  mcZ0,
205  dBw,
206  0,
207  1,
208  r'$|z_{\rm MC}|\ [{\rm ' + unitImp + '}]$',
209  r"{\rm " + unitImp + r"}\, "],
210  'mcGendistance': [
211  mcX0,
212  dBw,
213  0,
214  1,
215  r'$\sqrt{x_{\rm MC}^2 + y_{\rm MC}^2 + z_{\rm MC}^2}\ [{\rm ' + unitImp + '}]$',
216  r"{\rm " + unitImp + r"}\, "],
217  }
218 
219  withOrWithoutCut = ""
220  if cutOnUpsilonFourS != "":
221  withOrWithoutCut = "WithUpsilonFourSCut"
222  unitImp = "mm"
223  variablesPlotParamsDict = {'ImpactXY': [D0, dBw, 0, 1.0, r'$d_0\ [{\rm ' + unitImp + '}]$', r"{\rm " + unitImp + r"}\, "],
224  'z0': [Z0, dBw, 0, 1.0, r'$z_0\ [{\rm ' + unitImp + '}]$', r"{\rm " + unitImp + r"}\, "],
225  'distance': [X0, dBw, 0, 1.5, r'$\xi_0\ [{\rm ' + unitImp + '}]$', r"{\rm " + unitImp + r"}\, "],
226  'mcGenX': [mcX, dBw, 0, 1, r'$x_{\rm MC}\ [{\rm ' + unitImp + '}]$',
227  r"{\rm " + unitImp + r"}\, "],
228  'mcGenY': [mcY, dBw, -1, 1, r'$y_{\rm MC}\ [{\rm ' + unitImp + '}]$',
229  r"{\rm " + unitImp + r"}\, "],
230  'mcGenZ': [mcZ, dBw, -0.5, 1, r'$z_{\rm MC}\ [{\rm ' + unitImp + '}]$',
231  r"{\rm " + unitImp + r"}\, "],
232  'mcGenImpactXY': [mcD0, dBw, 0, 1.0,
233  r'$\sqrt{x_{\rm MC}^2 + y_{\rm MC}^2}\ [{\rm ' + unitImp + '}]$',
234  r"{\rm " + unitImp + r"}\, "],
235  'mcGenz0': [mcZ0, dBw, 0, 0.4,
236  r'$|z_{\rm MC}|\ [{\rm ' + unitImp + '}]$', r"{\rm " + unitImp + r"}\, "],
237  'mcGendistance': [mcX0, dBw, 0, 1.5,
238  r'$\sqrt{x_{\rm MC}^2 + y_{\rm MC}^2 + z_{\rm MC}^2}\ [{\rm ' +
239  unitImp + '}]$',
240  r"{\rm " + unitImp + r"}\, "],
241  }
242 
243  for Particle in Particles:
244  print(Particle)
245  for inputVariable in variablesPlotParamsDict:
246 
247  print(inputVariable)
248 
249  nBins = variablesPlotParamsDict[inputVariable][1]
250  limXmin = variablesPlotParamsDict[inputVariable][2]
251  limXmax = variablesPlotParamsDict[inputVariable][3]
252 
253  negativeHistogram = ROOT.TH1F(
254  "negative" + Particle + str(Belle2.MakeROOTCompatible.makeROOTCompatible(inputVariable)),
255  "", nBins, limXmin, limXmax)
256  positiveHistogram = ROOT.TH1F(
257  "positive" + Particle + str(Belle2.MakeROOTCompatible.makeROOTCompatible(inputVariable)),
258  "", nBins, limXmin, limXmax)
259  condition = "MCParticles.m_status%2==1 && abs(MCParticles.m_pdg) == " + particleConditions[Particle][1] + " && "
260  condition = condition + cutOnUpsilonFourS
261  condition = condition + " abs(MCParticles.m_pdg[MCParticles.m_mother - 1])==511 && "
262 
263  factorMultiplication = str()
264 
265  if belleOrBelle2 != "Belle" or (belleOrBelle2 == "Belle" and cutOnUpsilonFourS != ""):
266  factorMultiplication = "*10 "
267 
268  tree.Draw(variablesPlotParamsDict[inputVariable][0] +
269  factorMultiplication +
270  ">> negative" +
271  Particle +
272  str(Belle2.MakeROOTCompatible.makeROOTCompatible(inputVariable)), condition +
273  particleConditions[Particle][2])
274 
275  tree.Draw(variablesPlotParamsDict[inputVariable][0] +
276  factorMultiplication +
277  ">> positive" +
278  Particle +
279  str(Belle2.MakeROOTCompatible.makeROOTCompatible(inputVariable)), condition +
280  particleConditions[Particle][3])
281 
282  negativeScalingFactor = negativeHistogram.Integral()
283  positiveScalingFactor = positiveHistogram.Integral()
284 
285  if negativeScalingFactor == 0:
286  negativeScalingFactor = 1
287 
288  if positiveScalingFactor == 0:
289  positiveScalingFactor = 1
290 
291  negativeHistogram.Scale(1 / negativeScalingFactor)
292  positiveHistogram.Scale(1 / positiveScalingFactor)
293 
294  negativeArray = np.zeros((negativeHistogram.GetNbinsX(), 2))
295  positiveArray = np.zeros((positiveHistogram.GetNbinsX(), 2))
296 
297  asymmetryArray = np.zeros((positiveHistogram.GetNbinsX(), 2))
298  asymmetryArrayUncertainties = np.zeros((positiveHistogram.GetNbinsX(), 2))
299 
300  for i in range(0, negativeHistogram.GetNbinsX()):
301  negativeArray[i] = np.array([negativeHistogram.GetBinCenter(i + 1), negativeHistogram.GetBinContent(i + 1)])
302  positiveArray[i] = np.array([positiveHistogram.GetBinCenter(i + 1), positiveHistogram.GetBinContent(i + 1)])
303 
304  numerator = float(positiveArray[i][1] - negativeArray[i][1])
305  denominator = float(positiveArray[i][1] + negativeArray[i][1])
306 
307  if denominator == 0:
308  asymmetryArray[i] = np.array([negativeHistogram.GetBinCenter(i + 1), 0])
309  asymmetryArrayUncertainties[i] = np.array([negativeHistogram.GetBinCenter(i + 1), 0])
310  else:
311  asymmetryArray[i] = np.array([negativeHistogram.GetBinCenter(i + 1), float(numerator / denominator)])
312 
313  uncertainty = 2 * math.sqrt((negativeArray[i][1] * negativeHistogram.GetBinError(i + 1))**2 +
314  (positiveArray[i][1] * positiveHistogram.GetBinError(i + 1))**2) / \
315  (negativeArray[i][1] + positiveArray[i][1])**2
316 
317  asymmetryArrayUncertainties[i] = np.array([negativeHistogram.GetBinCenter(i + 1), uncertainty])
318 
319  fig1 = plt.figure(1, figsize=(11, 11))
320 
321  ax1 = plt.axes([0.19, 0.37, 0.76, 0.60])
322 
323  ax1.hist(
324  negativeArray[:, 0], weights=negativeArray[:, 1], bins=negativeHistogram.GetNbinsX(),
325  histtype='step', edgecolor='b', linewidth=4.5, linestyle='dashed',
326 
327  label=r'$' + particleConditions[Particle][0] + '^{-} $')
328 
329  ax1.hist(positiveArray[:, 0], weights=positiveArray[:, 1], bins=positiveHistogram.GetNbinsX(),
330  histtype='step', edgecolor='r', linewidth=4, alpha=0.9,
331  label=r'$' + particleConditions[Particle][0] + '^{+} $') # hatch='.',
332 
333  p1, = ax1.plot([], label=r'$' + particleConditions[Particle][0] + '^{-} $', linewidth=5.5, linestyle='dashed', c='b')
334  p2, = ax1.plot([], label=r'$' + particleConditions[Particle][0] +
335  '^{+} $', linewidth=5, linestyle='solid', alpha=0.9, c='r')
336 
337  binWidth = negativeHistogram.GetBinWidth(2)
338 
339  binWidth = '{:8.3f}'.format(binWidth)
340 
341  legendLocation = 1
342  ax1.set_ylabel(r'${\rm Fraction\hspace{0.25em} of\hspace{0.25em} Events}\, /\, (\, ' + binWidth + r'\, ' +
343  variablesPlotParamsDict[inputVariable][5] + r')$', fontsize=35)
344  if inputVariable == 'p' or inputVariable == 'pt':
345  legendLocation = 7
346 
347  ax1.yaxis.set_major_formatter(FormatStrFormatter(r'$%.2f$'))
348 
349  ax1.set_xlim(limXmin, limXmax)
350 
351  locs, labels = plt.xticks()
352 
353  empty_string_labels = [''] * len(labels)
354  plt.locator_params(axis='x', nbins=len(labels))
355  ax1.set_xticklabels(empty_string_labels)
356  ax1.tick_params(axis='y', labelsize=37)
357 
358  ax1.legend([p1, p2], [r'$' +
359  particleConditions[Particle][0] +
360  '^{-} $', r'$' +
361  particleConditions[Particle][0] +
362  '^{+} $'], prop={'size': 50}, loc=legendLocation, numpoints=1, handlelength=1)
363  ax1.grid(True)
364 
365  ax2 = plt.axes([0.19, 0.15, 0.76, 0.2])
366 
367  ax2.errorbar(asymmetryArray[:, 0], asymmetryArray[:, 1], xerr=float(negativeHistogram.GetBinWidth(2) / 2),
368  yerr=asymmetryArrayUncertainties[:, 1], elinewidth=2, mew=2, ecolor='k',
369  fmt='o', mfc='k', mec='k', markersize=6, label=r'${\rm Data}$')
370 
371  ax2.set_ylabel(r'$\frac{f^{+}\; -\;\, f^{-}}{f^{+}\; +\;\, f^{-}}$', fontsize=50, labelpad=20)
372  ax2.yaxis.labelpad = 8
373 
374  plt.yticks([-0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75],
375  [r'', r'$-0.5$', r'', r'$0.0$', r'', r'$0.5$', r''], rotation=0, size=25)
376  ax2.set_ylim(-0.75, 0.75)
377  ax2.xaxis.grid(True)
378  plt.axhline(y=0.5, linewidth=2, color='tab:gray', linestyle='-')
379  plt.axhline(y=0.25, linewidth=2, color='tab:gray', linestyle='-')
380  plt.axhline(y=0, linewidth=2, color='tab:gray', linestyle='-')
381  plt.axhline(y=-0.25, linewidth=2, color='tab:gray', linestyle='-')
382  plt.axhline(y=-0.5, linewidth=2, color='tab:gray', linestyle='-')
383 
384  xLabel = variablesPlotParamsDict[inputVariable][4]
385 
386  ax2.xaxis.labelpad = 15
387  plt.locator_params(axis='x', nbins=len(labels))
388  ax2.tick_params(axis='x', labelsize=40)
389 
390  ax2.set_xlim(limXmin, limXmax)
391 
392  if inputVariable == 'mcGenY' and belleOrBelle2 == "Belle2":
393  plt.xticks([-0.06, -0.04, -0.02, 0, 0.02, 0.04, 0.06],
394  ['', r'$-0.04$', '', r'$0.00$', '', r'$0.04$', ''], rotation=0, size=40)
395 
396  ax2.set_xlabel(xLabel, fontsize=50)
397  plt.savefig(workingDirectory + '/GenLevelVariablesPlots' + '/' + belleOrBelle2 + withOrWithoutCut +
398  "_" + Particle + "_" + str(Belle2.MakeROOTCompatible.makeROOTCompatible(inputVariable)) + '.pdf')
399  fig1.clear()
400 
401  negativeHistogram.Delete()
402  positiveHistogram.Delete()
403 
404 
405 DecZ = "MCParticles.m_decayVertex_z"
406 
407 if belleOrBelle2 == "Belle2":
408 
409  variablesPlotB0tag = {'DecZ': [DecZ, dBw, -0.6, 1.5, r"{\rm " + unitImp + r"}\, "],
410  }
411 
412 else:
413  unitImp = "cm"
414  variablesPlotB0tag = {'DecZ': [DecZ, dBw, -1.5, 1.5, r"{\rm " + unitImp + r"}\, "],
415  }
416 
417 particleCondsB0tag = {
418  'B0tag': [
419  'B0',
420  "511 ",
421  " MCParticles.m_pdg < 0 ",
422  " MCParticles.m_pdg > 0 ",
423  " !(MCParticles.m_pdg[MCParticles.m_firstDaughter-1] == 310 && MCParticles.m_pdg[MCParticles.m_firstDaughter] == 443) && ",
424  r'$(z_{\rm tag}^{\rm dec})^{\rm gen}\ [{\rm ' + unitImp + '}]$'],
425  'B0sig': [
426  'B0',
427  "511 ",
428  " MCParticles.m_pdg < 0 ",
429  " MCParticles.m_pdg > 0 ",
430  " (MCParticles.m_pdg[MCParticles.m_firstDaughter-1] == 310 && MCParticles.m_pdg[MCParticles.m_firstDaughter] == 443) && ",
431  r'$(z_{\rm sig}^{\rm dec})^{\rm gen}\ [{\rm ' + unitImp + '}]$']}
432 
433 
434 def makeZtagDecayPlot(cutOnUpsilonFourS="", asymplot=True):
435  """
436  Plots the distribution of the decay z-vertex of the tag-side B0 meson.
437  """
438 
439  withOrWithoutCut = ""
440 
441  if cutOnUpsilonFourS != "":
442  withOrWithoutCut = "WithUpsilonFourSCut"
443  unitImp = "mm"
444  variablesPlotB0tag['DecZ'] = [DecZ, dBw, -0.6, 1.5, r"{\rm " + unitImp + r"}\, "]
445  particleCondsB0tag['B0tag'][5] = r'$(z_{\rm tag}^{\rm dec})^{\rm gen}\ [{\rm ' + unitImp + r'}]$'
446 
447  for Particle in particleCondsB0tag:
448  inputVariable = "DecZ"
449 
450  print("Plotting " + Particle + " z-decay Vertex")
451 
452  nBins = variablesPlotB0tag[inputVariable][1]
453  limXmin = variablesPlotB0tag[inputVariable][2]
454  limXmax = variablesPlotB0tag[inputVariable][3]
455 
456  negativeHistogram = ROOT.TH1F("negative" + Particle + str(Belle2.MakeROOTCompatible.makeROOTCompatible(inputVariable)), "",
457  nBins,
458  limXmin,
459  limXmax)
460  positiveHistogram = ROOT.TH1F("positive" + Particle + str(Belle2.MakeROOTCompatible.makeROOTCompatible(inputVariable)), "",
461  nBins,
462  limXmin,
463  limXmax)
464  condition = "abs(MCParticles.m_pdg) == " + particleCondsB0tag[Particle][1] + " && "
465  condition = condition + cutOnUpsilonFourS
466  condition = condition + particleCondsB0tag[Particle][4]
467 
468  factorMultiplication = str()
469 
470  if belleOrBelle2 != "Belle" or (belleOrBelle2 == "Belle" and cutOnUpsilonFourS != ""):
471  factorMultiplication = "*10 "
472 
473  tree.Draw(variablesPlotB0tag[inputVariable][0] + factorMultiplication + ">> negative" + Particle +
474  str(Belle2.MakeROOTCompatible.makeROOTCompatible(inputVariable)), condition + particleCondsB0tag[Particle][2])
475 
476  tree.Draw(variablesPlotB0tag[inputVariable][0] + factorMultiplication + ">> positive" + Particle +
477  str(Belle2.MakeROOTCompatible.makeROOTCompatible(inputVariable)), condition + particleCondsB0tag[Particle][3])
478 
479  negativeScalingFactor = negativeHistogram.Integral()
480  positiveScalingFactor = positiveHistogram.Integral()
481 
482  if negativeScalingFactor == 0:
483  negativeScalingFactor = 1
484 
485  if positiveScalingFactor == 0:
486  positiveScalingFactor = 1
487 
488  negativeHistogram.Scale(1 / negativeScalingFactor)
489  positiveHistogram.Scale(1 / positiveScalingFactor)
490 
491  negativeArray = np.zeros((negativeHistogram.GetNbinsX(), 2))
492  positiveArray = np.zeros((positiveHistogram.GetNbinsX(), 2))
493 
494  asymmetryArray = np.zeros((positiveHistogram.GetNbinsX(), 2))
495  asymmetryArrayUncertainties = np.zeros((positiveHistogram.GetNbinsX(), 2))
496 
497  for i in range(0, negativeHistogram.GetNbinsX()):
498  negativeArray[i] = np.array([negativeHistogram.GetBinCenter(i + 1), negativeHistogram.GetBinContent(i + 1)])
499  positiveArray[i] = np.array([positiveHistogram.GetBinCenter(i + 1), positiveHistogram.GetBinContent(i + 1)])
500 
501  numerator = float(positiveArray[i][1] - negativeArray[i][1])
502  denominator = float(positiveArray[i][1] + negativeArray[i][1])
503 
504  if denominator == 0:
505  asymmetryArray[i] = np.array([negativeHistogram.GetBinCenter(i + 1), 0])
506  asymmetryArrayUncertainties[i] = np.array([negativeHistogram.GetBinCenter(i + 1), 0])
507  else:
508  asymmetryArray[i] = np.array([negativeHistogram.GetBinCenter(i + 1), float(numerator / denominator)])
509 
510  uncertainty = 2 * math.sqrt((negativeArray[i][1] * negativeHistogram.GetBinError(i + 1))**2 + (
511  positiveArray[i][1] * positiveHistogram.GetBinError(i + 1))**2) / (negativeArray[i][1] + positiveArray[i][1])**2
512 
513  asymmetryArrayUncertainties[i] = np.array([negativeHistogram.GetBinCenter(i + 1), uncertainty])
514 
515  fig1 = plt.figure(1, figsize=(11, 11))
516 
517  ax1 = plt.axes([0.19, 0.37, 0.76, 0.60])
518 
519  ax1.hist(
520  negativeArray[:, 0], weights=negativeArray[:, 1], bins=negativeHistogram.GetNbinsX(),
521  histtype='step', edgecolor='b', linewidth=4.5, linestyle='dashed',
522 
523  label=r'$' + particleCondsB0tag[Particle][0] + '^{-} $')
524 
525  ax1.hist(positiveArray[:, 0], weights=positiveArray[:, 1], bins=positiveHistogram.GetNbinsX(),
526  histtype='step', edgecolor='r', linewidth=4, alpha=0.9,
527  label=r'$' + particleCondsB0tag[Particle][0] + '^{+} $') # hatch='.',
528 
529  p1, = ax1.plot([], label=r'$\bar{B}^0$', linewidth=4.5, linestyle='dashed', c='b')
530  p2, = ax1.plot([], label=r'$B^0$', linewidth=4, linestyle='solid', alpha=0.9, c='r')
531 
532  binWidth = negativeHistogram.GetBinWidth(2)
533 
534  binWidth = '{:8.3f}'.format(binWidth)
535 
536  legendLocation = 1
537  ax1.set_ylabel(r'${\rm Fraction\hspace{0.25em} of\hspace{0.25em} Events}\, /\, (\, ' + binWidth + r'\, ' +
538  variablesPlotB0tag[inputVariable][4] + r')$', fontsize=35)
539  if inputVariable == 'p' or inputVariable == 'pt':
540  legendLocation = 7
541 
542  ax1.yaxis.set_major_formatter(FormatStrFormatter(r'$%.2f$'))
543 
544  ax1.set_xlim(limXmin, limXmax)
545 
546  locs, labels = plt.xticks()
547  ax1.tick_params(axis='y', labelsize=37)
548 
549  ax1.legend([p1, p2], [r'$\bar{B}^0$', r'$B^0$'], prop={'size': 50}, loc=legendLocation, numpoints=1, handlelength=1)
550  ax1.grid(True)
551  xLabel = particleCondsB0tag[Particle][5]
552 
553  if not asymplot:
554  ax1.xaxis.labelpad = 15
555  plt.locator_params(axis='x', nbins=len(labels))
556  ax1.tick_params(axis='x', labelsize=40)
557  ax1.set_xlim(limXmin, limXmax)
558  ax1.set_xlabel(xLabel, fontsize=50)
559 
560  else:
561 
562  empty_string_labels = [''] * len(labels)
563  plt.locator_params(axis='x', nbins=len(labels))
564  ax1.set_xticklabels(empty_string_labels)
565 
566  ax2 = plt.axes([0.19, 0.15, 0.76, 0.2])
567 
568  ax2.errorbar(asymmetryArray[:, 0], asymmetryArray[:, 1], xerr=float(negativeHistogram.GetBinWidth(2) / 2),
569  yerr=asymmetryArrayUncertainties[:, 1], elinewidth=2, mew=2, ecolor='k',
570  fmt='o', mfc='k', mec='k', markersize=6, label=r'${\rm Data}$')
571 
572  ax2.set_ylabel(
573  r'$\frac{f^{B^0}\; -\;\, f^{\overline{B}^0}}{f^{B^0}\; +\;\, f^{\overline{B}^0}}$',
574  fontsize=46.5,
575  labelpad=20)
576  ax2.yaxis.labelpad = 0
577 
578  plt.yticks([-0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75],
579  [r'', r'$-0.5$', r'', r'$0.0$', r'', r'$0.5$', r''], rotation=0, size=25)
580  ax2.set_ylim(-0.75, 0.75)
581  ax2.xaxis.grid(True)
582  plt.axhline(y=0.5, linewidth=2, color='tab:gray', linestyle='-')
583  plt.axhline(y=0.25, linewidth=2, color='tab:gray', linestyle='-')
584  plt.axhline(y=0, linewidth=2, color='tab:gray', linestyle='-')
585  plt.axhline(y=-0.25, linewidth=2, color='tab:gray', linestyle='-')
586  plt.axhline(y=-0.5, linewidth=2, color='tab:gray', linestyle='-')
587 
588  ax2.xaxis.labelpad = 15
589  plt.locator_params(axis='x', nbins=len(labels))
590  ax2.tick_params(axis='x', labelsize=40)
591  ax2.set_xlim(limXmin, limXmax)
592  ax2.set_xlabel(xLabel, fontsize=50)
593 
594  plt.savefig(workingDirectory + '/GenLevelVariablesPlots' + '/' +
595  belleOrBelle2 + withOrWithoutCut + "_" + Particle + "_" +
596  str(Belle2.MakeROOTCompatible.makeROOTCompatible(inputVariable)) + '.pdf')
597  fig1.clear()
598 
599  negativeHistogram.Delete()
600  positiveHistogram.Delete()
601 
602 
603 makePlotsForEachParticleKind("")
604 
605 if belleOrBelle2 == "Belle":
606  makeZtagDecayPlot("")
607  makeZtagDecayPlot(" abs(MCParticles.m_decayVertex_z[MCParticles.m_mother-1]) < 0.04 && ")
608  makePlotsForEachParticleKind(" abs(MCParticles.m_decayVertex_z[MCParticles.m_mother[MCParticles.m_mother - 1]-1]) < 0.04 && ")
609 else:
610  makeZtagDecayPlot("")
Belle2::MakeROOTCompatible::makeROOTCompatible
static std::string makeROOTCompatible(std::string str)
Remove special characters that ROOT dislikes in branch names, e.g.
Definition: MakeROOTCompatible.cc:60