Belle II Software  release-06-02-00
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("negative" + Particle + str(Belle2.makeROOTCompatible(inputVariable)), "",
254  nBins,
255  limXmin,
256  limXmax)
257  positiveHistogram = ROOT.TH1F("positive" + Particle + str(Belle2.makeROOTCompatible(inputVariable)), "",
258  nBins,
259  limXmin,
260  limXmax)
261  condition = "MCParticles.m_status%2==1 && abs(MCParticles.m_pdg) == " + particleConditions[Particle][1] + " && "
262  condition = condition + cutOnUpsilonFourS
263  condition = condition + " abs(MCParticles.m_pdg[MCParticles.m_mother - 1])==511 && "
264 
265  factorMultiplication = str()
266 
267  if belleOrBelle2 != "Belle" or (belleOrBelle2 == "Belle" and cutOnUpsilonFourS != ""):
268  factorMultiplication = "*10 "
269 
270  tree.Draw(variablesPlotParamsDict[inputVariable][0] +
271  factorMultiplication +
272  ">> negative" +
273  Particle +
274  str(Belle2.makeROOTCompatible(inputVariable)), condition +
275  particleConditions[Particle][2])
276 
277  tree.Draw(variablesPlotParamsDict[inputVariable][0] +
278  factorMultiplication +
279  ">> positive" +
280  Particle +
281  str(Belle2.makeROOTCompatible(inputVariable)), condition +
282  particleConditions[Particle][3])
283 
284  negativeScalingFactor = negativeHistogram.Integral()
285  positiveScalingFactor = positiveHistogram.Integral()
286 
287  if negativeScalingFactor == 0:
288  negativeScalingFactor = 1
289 
290  if positiveScalingFactor == 0:
291  positiveScalingFactor = 1
292 
293  negativeHistogram.Scale(1 / negativeScalingFactor)
294  positiveHistogram.Scale(1 / positiveScalingFactor)
295 
296  negativeArray = np.zeros((negativeHistogram.GetNbinsX(), 2))
297  positiveArray = np.zeros((positiveHistogram.GetNbinsX(), 2))
298 
299  asymmetryArray = np.zeros((positiveHistogram.GetNbinsX(), 2))
300  asymmetryArrayUncertainties = np.zeros((positiveHistogram.GetNbinsX(), 2))
301 
302  for i in range(0, negativeHistogram.GetNbinsX()):
303  negativeArray[i] = np.array([negativeHistogram.GetBinCenter(i + 1), negativeHistogram.GetBinContent(i + 1)])
304  positiveArray[i] = np.array([positiveHistogram.GetBinCenter(i + 1), positiveHistogram.GetBinContent(i + 1)])
305 
306  numerator = float(positiveArray[i][1] - negativeArray[i][1])
307  denominator = float(positiveArray[i][1] + negativeArray[i][1])
308 
309  if denominator == 0:
310  asymmetryArray[i] = np.array([negativeHistogram.GetBinCenter(i + 1), 0])
311  asymmetryArrayUncertainties[i] = np.array([negativeHistogram.GetBinCenter(i + 1), 0])
312  else:
313  asymmetryArray[i] = np.array([negativeHistogram.GetBinCenter(i + 1), float(numerator / denominator)])
314 
315  uncertainty = 2 * math.sqrt((negativeArray[i][1] * negativeHistogram.GetBinError(i + 1))**2 +
316  (positiveArray[i][1] * positiveHistogram.GetBinError(i + 1))**2) / \
317  (negativeArray[i][1] + positiveArray[i][1])**2
318 
319  asymmetryArrayUncertainties[i] = np.array([negativeHistogram.GetBinCenter(i + 1), uncertainty])
320 
321  fig1 = plt.figure(1, figsize=(11, 11))
322 
323  ax1 = plt.axes([0.19, 0.37, 0.76, 0.60])
324 
325  ax1.hist(
326  negativeArray[:, 0], weights=negativeArray[:, 1], bins=negativeHistogram.GetNbinsX(),
327  histtype='step', edgecolor='b', linewidth=4.5, linestyle='dashed',
328 
329  label=r'$' + particleConditions[Particle][0] + '^{-} $')
330 
331  ax1.hist(positiveArray[:, 0], weights=positiveArray[:, 1], bins=positiveHistogram.GetNbinsX(),
332  histtype='step', edgecolor='r', linewidth=4, alpha=0.9,
333  label=r'$' + particleConditions[Particle][0] + '^{+} $') # hatch='.',
334 
335  p1, = ax1.plot([], label=r'$' + particleConditions[Particle][0] + '^{-} $', linewidth=5.5, linestyle='dashed', c='b')
336  p2, = ax1.plot([], label=r'$' + particleConditions[Particle][0] +
337  '^{+} $', linewidth=5, linestyle='solid', alpha=0.9, c='r')
338 
339  binWidth = negativeHistogram.GetBinWidth(2)
340 
341  binWidth = '{:8.3f}'.format(binWidth)
342 
343  legendLocation = 1
344  ax1.set_ylabel(r'${\rm Fraction\hspace{0.25em} of\hspace{0.25em} Events}\, /\, (\, ' + binWidth + r'\, ' +
345  variablesPlotParamsDict[inputVariable][5] + r')$', fontsize=35)
346  if inputVariable == 'p' or inputVariable == 'pt':
347  legendLocation = 7
348 
349  ax1.yaxis.set_major_formatter(FormatStrFormatter(r'$%.2f$'))
350 
351  ax1.set_xlim(limXmin, limXmax)
352 
353  locs, labels = plt.xticks()
354 
355  empty_string_labels = [''] * len(labels)
356  plt.locator_params(axis='x', nbins=len(labels))
357  ax1.set_xticklabels(empty_string_labels)
358  ax1.tick_params(axis='y', labelsize=37)
359 
360  ax1.legend([p1, p2], [r'$' +
361  particleConditions[Particle][0] +
362  '^{-} $', r'$' +
363  particleConditions[Particle][0] +
364  '^{+} $'], prop={'size': 50}, loc=legendLocation, numpoints=1, handlelength=1)
365  ax1.grid(True)
366 
367  ax2 = plt.axes([0.19, 0.15, 0.76, 0.2])
368 
369  ax2.errorbar(asymmetryArray[:, 0], asymmetryArray[:, 1], xerr=float(negativeHistogram.GetBinWidth(2) / 2),
370  yerr=asymmetryArrayUncertainties[:, 1], elinewidth=2, mew=2, ecolor='k',
371  fmt='o', mfc='k', mec='k', markersize=6, label=r'${\rm Data}$')
372 
373  ax2.set_ylabel(r'$\frac{f^{+}\; -\;\, f^{-}}{f^{+}\; +\;\, f^{-}}$', fontsize=50, labelpad=20)
374  ax2.yaxis.labelpad = 8
375 
376  plt.yticks([-0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75],
377  [r'', r'$-0.5$', r'', r'$0.0$', r'', r'$0.5$', r''], rotation=0, size=25)
378  ax2.set_ylim(-0.75, 0.75)
379  ax2.xaxis.grid(True)
380  plt.axhline(y=0.5, linewidth=2, color='tab:gray', linestyle='-')
381  plt.axhline(y=0.25, linewidth=2, color='tab:gray', linestyle='-')
382  plt.axhline(y=0, linewidth=2, color='tab:gray', linestyle='-')
383  plt.axhline(y=-0.25, linewidth=2, color='tab:gray', linestyle='-')
384  plt.axhline(y=-0.5, linewidth=2, color='tab:gray', linestyle='-')
385 
386  xLabel = variablesPlotParamsDict[inputVariable][4]
387 
388  ax2.xaxis.labelpad = 15
389  plt.locator_params(axis='x', nbins=len(labels))
390  ax2.tick_params(axis='x', labelsize=40)
391 
392  ax2.set_xlim(limXmin, limXmax)
393 
394  if inputVariable == 'mcGenY' and belleOrBelle2 == "Belle2":
395  plt.xticks([-0.06, -0.04, -0.02, 0, 0.02, 0.04, 0.06],
396  ['', r'$-0.04$', '', r'$0.00$', '', r'$0.04$', ''], rotation=0, size=40)
397 
398  ax2.set_xlabel(xLabel, fontsize=50)
399  plt.savefig(workingDirectory + '/GenLevelVariablesPlots' + '/' + belleOrBelle2 + withOrWithoutCut +
400  "_" + Particle + "_" + str(Belle2.makeROOTCompatible(inputVariable)) + '.pdf')
401  fig1.clear()
402 
403  negativeHistogram.Delete()
404  positiveHistogram.Delete()
405 
406 
407 DecZ = "MCParticles.m_decayVertex_z"
408 
409 if belleOrBelle2 == "Belle2":
410 
411  variablesPlotB0tag = {'DecZ': [DecZ, dBw, -0.6, 1.5, r"{\rm " + unitImp + r"}\, "],
412  }
413 
414 else:
415  unitImp = "cm"
416  variablesPlotB0tag = {'DecZ': [DecZ, dBw, -1.5, 1.5, r"{\rm " + unitImp + r"}\, "],
417  }
418 
419 particleCondsB0tag = {
420  'B0tag': [
421  'B0',
422  "511 ",
423  " MCParticles.m_pdg < 0 ",
424  " MCParticles.m_pdg > 0 ",
425  " !(MCParticles.m_pdg[MCParticles.m_firstDaughter-1] == 310 && MCParticles.m_pdg[MCParticles.m_firstDaughter] == 443) && ",
426  r'$(z_{\rm tag}^{\rm dec})^{\rm gen}\ [{\rm ' + unitImp + '}]$'],
427  'B0sig': [
428  'B0',
429  "511 ",
430  " MCParticles.m_pdg < 0 ",
431  " MCParticles.m_pdg > 0 ",
432  " (MCParticles.m_pdg[MCParticles.m_firstDaughter-1] == 310 && MCParticles.m_pdg[MCParticles.m_firstDaughter] == 443) && ",
433  r'$(z_{\rm sig}^{\rm dec})^{\rm gen}\ [{\rm ' + unitImp + '}]$']}
434 
435 
436 def makeZtagDecayPlot(cutOnUpsilonFourS="", asymplot=True):
437  """
438  Plots the distribution of the decay z-vertex of the tag-side B0 meson.
439  """
440 
441  withOrWithoutCut = ""
442 
443  if cutOnUpsilonFourS != "":
444  withOrWithoutCut = "WithUpsilonFourSCut"
445  unitImp = "mm"
446  variablesPlotB0tag['DecZ'] = [DecZ, dBw, -0.6, 1.5, r"{\rm " + unitImp + r"}\, "]
447  particleCondsB0tag['B0tag'][5] = r'$(z_{\rm tag}^{\rm dec})^{\rm gen}\ [{\rm ' + unitImp + r'}]$'
448 
449  for Particle in particleCondsB0tag:
450  inputVariable = "DecZ"
451 
452  print("Plotting " + Particle + " z-decay Vertex")
453 
454  nBins = variablesPlotB0tag[inputVariable][1]
455  limXmin = variablesPlotB0tag[inputVariable][2]
456  limXmax = variablesPlotB0tag[inputVariable][3]
457 
458  negativeHistogram = ROOT.TH1F("negative" + Particle + str(Belle2.makeROOTCompatible(inputVariable)), "",
459  nBins,
460  limXmin,
461  limXmax)
462  positiveHistogram = ROOT.TH1F("positive" + Particle + str(Belle2.makeROOTCompatible(inputVariable)), "",
463  nBins,
464  limXmin,
465  limXmax)
466  condition = "abs(MCParticles.m_pdg) == " + particleCondsB0tag[Particle][1] + " && "
467  condition = condition + cutOnUpsilonFourS
468  condition = condition + particleCondsB0tag[Particle][4]
469 
470  factorMultiplication = str()
471 
472  if belleOrBelle2 != "Belle" or (belleOrBelle2 == "Belle" and cutOnUpsilonFourS != ""):
473  factorMultiplication = "*10 "
474 
475  tree.Draw(variablesPlotB0tag[inputVariable][0] + factorMultiplication + ">> negative" + Particle +
476  str(Belle2.makeROOTCompatible(inputVariable)), condition + particleCondsB0tag[Particle][2])
477 
478  tree.Draw(variablesPlotB0tag[inputVariable][0] + factorMultiplication + ">> positive" + Particle +
479  str(Belle2.makeROOTCompatible(inputVariable)), condition + particleCondsB0tag[Particle][3])
480 
481  negativeScalingFactor = negativeHistogram.Integral()
482  positiveScalingFactor = positiveHistogram.Integral()
483 
484  if negativeScalingFactor == 0:
485  negativeScalingFactor = 1
486 
487  if positiveScalingFactor == 0:
488  positiveScalingFactor = 1
489 
490  negativeHistogram.Scale(1 / negativeScalingFactor)
491  positiveHistogram.Scale(1 / positiveScalingFactor)
492 
493  negativeArray = np.zeros((negativeHistogram.GetNbinsX(), 2))
494  positiveArray = np.zeros((positiveHistogram.GetNbinsX(), 2))
495 
496  asymmetryArray = np.zeros((positiveHistogram.GetNbinsX(), 2))
497  asymmetryArrayUncertainties = np.zeros((positiveHistogram.GetNbinsX(), 2))
498 
499  for i in range(0, negativeHistogram.GetNbinsX()):
500  negativeArray[i] = np.array([negativeHistogram.GetBinCenter(i + 1), negativeHistogram.GetBinContent(i + 1)])
501  positiveArray[i] = np.array([positiveHistogram.GetBinCenter(i + 1), positiveHistogram.GetBinContent(i + 1)])
502 
503  numerator = float(positiveArray[i][1] - negativeArray[i][1])
504  denominator = float(positiveArray[i][1] + negativeArray[i][1])
505 
506  if denominator == 0:
507  asymmetryArray[i] = np.array([negativeHistogram.GetBinCenter(i + 1), 0])
508  asymmetryArrayUncertainties[i] = np.array([negativeHistogram.GetBinCenter(i + 1), 0])
509  else:
510  asymmetryArray[i] = np.array([negativeHistogram.GetBinCenter(i + 1), float(numerator / denominator)])
511 
512  uncertainty = 2 * math.sqrt((negativeArray[i][1] * negativeHistogram.GetBinError(i + 1))**2 + (
513  positiveArray[i][1] * positiveHistogram.GetBinError(i + 1))**2) / (negativeArray[i][1] + positiveArray[i][1])**2
514 
515  asymmetryArrayUncertainties[i] = np.array([negativeHistogram.GetBinCenter(i + 1), uncertainty])
516 
517  fig1 = plt.figure(1, figsize=(11, 11))
518 
519  ax1 = plt.axes([0.19, 0.37, 0.76, 0.60])
520 
521  ax1.hist(
522  negativeArray[:, 0], weights=negativeArray[:, 1], bins=negativeHistogram.GetNbinsX(),
523  histtype='step', edgecolor='b', linewidth=4.5, linestyle='dashed',
524 
525  label=r'$' + particleCondsB0tag[Particle][0] + '^{-} $')
526 
527  ax1.hist(positiveArray[:, 0], weights=positiveArray[:, 1], bins=positiveHistogram.GetNbinsX(),
528  histtype='step', edgecolor='r', linewidth=4, alpha=0.9,
529  label=r'$' + particleCondsB0tag[Particle][0] + '^{+} $') # hatch='.',
530 
531  p1, = ax1.plot([], label=r'$\bar{B}^0$', linewidth=4.5, linestyle='dashed', c='b')
532  p2, = ax1.plot([], label=r'$B^0$', linewidth=4, linestyle='solid', alpha=0.9, c='r')
533 
534  binWidth = negativeHistogram.GetBinWidth(2)
535 
536  binWidth = '{:8.3f}'.format(binWidth)
537 
538  legendLocation = 1
539  ax1.set_ylabel(r'${\rm Fraction\hspace{0.25em} of\hspace{0.25em} Events}\, /\, (\, ' + binWidth + r'\, ' +
540  variablesPlotB0tag[inputVariable][4] + r')$', fontsize=35)
541  if inputVariable == 'p' or inputVariable == 'pt':
542  legendLocation = 7
543 
544  ax1.yaxis.set_major_formatter(FormatStrFormatter(r'$%.2f$'))
545 
546  ax1.set_xlim(limXmin, limXmax)
547 
548  locs, labels = plt.xticks()
549  ax1.tick_params(axis='y', labelsize=37)
550 
551  ax1.legend([p1, p2], [r'$\bar{B}^0$', r'$B^0$'], prop={'size': 50}, loc=legendLocation, numpoints=1, handlelength=1)
552  ax1.grid(True)
553  xLabel = particleCondsB0tag[Particle][5]
554 
555  if not asymplot:
556  ax1.xaxis.labelpad = 15
557  plt.locator_params(axis='x', nbins=len(labels))
558  ax1.tick_params(axis='x', labelsize=40)
559  ax1.set_xlim(limXmin, limXmax)
560  ax1.set_xlabel(xLabel, fontsize=50)
561 
562  else:
563 
564  empty_string_labels = [''] * len(labels)
565  plt.locator_params(axis='x', nbins=len(labels))
566  ax1.set_xticklabels(empty_string_labels)
567 
568  ax2 = plt.axes([0.19, 0.15, 0.76, 0.2])
569 
570  ax2.errorbar(asymmetryArray[:, 0], asymmetryArray[:, 1], xerr=float(negativeHistogram.GetBinWidth(2) / 2),
571  yerr=asymmetryArrayUncertainties[:, 1], elinewidth=2, mew=2, ecolor='k',
572  fmt='o', mfc='k', mec='k', markersize=6, label=r'${\rm Data}$')
573 
574  ax2.set_ylabel(
575  r'$\frac{f^{B^0}\; -\;\, f^{\overline{B}^0}}{f^{B^0}\; +\;\, f^{\overline{B}^0}}$',
576  fontsize=46.5,
577  labelpad=20)
578  ax2.yaxis.labelpad = 0
579 
580  plt.yticks([-0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75],
581  [r'', r'$-0.5$', r'', r'$0.0$', r'', r'$0.5$', r''], rotation=0, size=25)
582  ax2.set_ylim(-0.75, 0.75)
583  ax2.xaxis.grid(True)
584  plt.axhline(y=0.5, linewidth=2, color='tab:gray', linestyle='-')
585  plt.axhline(y=0.25, linewidth=2, color='tab:gray', linestyle='-')
586  plt.axhline(y=0, linewidth=2, color='tab:gray', linestyle='-')
587  plt.axhline(y=-0.25, linewidth=2, color='tab:gray', linestyle='-')
588  plt.axhline(y=-0.5, linewidth=2, color='tab:gray', linestyle='-')
589 
590  ax2.xaxis.labelpad = 15
591  plt.locator_params(axis='x', nbins=len(labels))
592  ax2.tick_params(axis='x', labelsize=40)
593  ax2.set_xlim(limXmin, limXmax)
594  ax2.set_xlabel(xLabel, fontsize=50)
595 
596  plt.savefig(workingDirectory + '/GenLevelVariablesPlots' + '/' +
597  belleOrBelle2 + withOrWithoutCut + "_" + Particle + "_" +
598  str(Belle2.makeROOTCompatible(inputVariable)) + '.pdf')
599  fig1.clear()
600 
601  negativeHistogram.Delete()
602  positiveHistogram.Delete()
603 
604 
605 makePlotsForEachParticleKind("")
606 
607 if belleOrBelle2 == "Belle":
608  makeZtagDecayPlot("")
609  makeZtagDecayPlot(" abs(MCParticles.m_decayVertex_z[MCParticles.m_mother-1]) < 0.04 && ")
610  makePlotsForEachParticleKind(" abs(MCParticles.m_decayVertex_z[MCParticles.m_mother[MCParticles.m_mother - 1]-1]) < 0.04 && ")
611 else:
612  makeZtagDecayPlot("")
Belle2::makeROOTCompatible
std::string makeROOTCompatible(std::string str)
Remove special characters that ROOT dislikes in branch names, e.g.
Definition: MakeROOTCompatible.cc:64