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