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