Belle II Software development
genLevelAsymmsImpactParams.py
1#!/usr/bin/env python3
2
3
10
11
21
22
23import basf2 as b2
24import os
25import glob
26import sys
27import math
28from matplotlib.ticker import FormatStrFormatter
29import matplotlib.pyplot as plt
30import ROOT
31from ROOT import Belle2
32
33import numpy as np
34import matplotlib as mpl
35mpl.use('Agg')
36mpl.rcParams.update({'font.size': 22})
37mpl.rcParams['text.usetex'] = True
38mpl.rcParams['text.latex.preamble'] = [r"\usepackage{amsmath}"]
39
40if len(sys.argv) != 3:
41 sys.exit("Must provide 2 arguments: [Belle or Belle2] [samplesWildCards]"
42 )
43
44belleOrBelle2 = sys.argv[1]
45samplesWildCards = sys.argv[2]
46
47
48workingDirectory = '.'
49
50if not b2.find_file('GenLevelVariablesPlots', silent=True):
51 os.mkdir(workingDirectory + '/GenLevelVariablesPlots')
52
53sourceFiles = glob.glob(str(samplesWildCards))
54
55tree = ROOT.TChain("tree")
56for Inputfile in sourceFiles:
57 tree.AddFile(Inputfile)
58
59
60ROOT.TH1.SetDefaultSumw2()
61
62binsPt = 50
63
64HistoBelle2_PtMC = ROOT.TH1F('PtMCBelle2', 'MC Transverse Momentum', binsPt, 0, 5)
65HistoBelle2_PtFit = ROOT.TH1F('PtFitBelle2', 'TrackFitResult Transverse Momentum', binsPt, 0, 5)
66
67HistoBelle2_pseudoPD0MC = ROOT.TH1F('pseudoPD0MCBelle2', 'MC pseudo Momentum for D0', binsPt, 0, 5)
68HistoBelle2_pseudoPD0Fit = ROOT.TH1F('pseudoPD0FitBelle2', 'Fit pseudo Momentum for D0', binsPt, 0, 5)
69
70HistoBelle2_pseudoPZ0MC = ROOT.TH1F('pseudoPZ0MCBelle2', 'MC pseudo Momentum for Z0', binsPt, 0, 5)
71HistoBelle2_pseudoPZ0Fit = ROOT.TH1F('pseudoPZ0FitBelle2', 'Fit pseudo Momentum for Z0', binsPt, 0, 5)
72
73HistoBelle2_d0Error = ROOT.TH1F('d0ErrorBelle2', 'TrackFitResult Error on d0', binsPt, 0, 0.03)
74HistoBelle2_z0Error = ROOT.TH1F('z0ErrorBelle2', 'TrackFitResult Error on z0', binsPt, 0, 0.03)
75
76D0 = 'abs(TrackFitResults.m_tau[0][0])'
77Z0 = 'abs(TrackFitResults.m_tau[0][3])'
78X0 = 'sqrt(TrackFitResults.m_tau[0][0]**2 + TrackFitResults.m_tau[0][3]**2)'
79
80mcX = "MCParticles.m_productionVertex_x"
81mcY = "MCParticles.m_productionVertex_y"
82mcZ = "MCParticles.m_productionVertex_z"
83
84mcD0 = "sqrt(MCParticles.m_productionVertex_x*MCParticles.m_productionVertex_x + " + \
85 " MCParticles.m_productionVertex_y*MCParticles.m_productionVertex_y)"
86mcZ0 = "abs(MCParticles.m_productionVertex_z)"
87mcX0 = "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
91ptMC = 'sqrt(MCParticles.m_momentum_x**2 + MCParticles.m_momentum_y**2)'
92ptFit = 'abs(0.003*1.5/TrackFitResults.m_tau[0][2])'
93
94pMC = 'sqrt(MCParticles.m_momentum_x**2 + MCParticles.m_momentum_y**2 + + MCParticles.m_momentum_z**2)'
95pFit = '(abs(0.003*1.5/TrackFitResults.m_tau[0][2])*sqrt(1+TrackFitResults.m_tau[0][4]**2))'
96chiProb = 'TrackFitResults.m_pValue'
97
98pseudoPD0MC = 'sqrt(' + ptMC + '**3/(' + pMC + '*(1+(MCParticles.m_mass[0]/' + pMC + ')**2)))'
99pseudoPD0Fit = 'sqrt(' + ptFit + '**3/(' + pFit + '*(1+(MCParticles.m_mass[0]/' + pFit + ')**2)))'
100
101pseudoPZ0MC = 'sqrt(' + ptMC + '**5/(' + pMC + '**3*(1+(MCParticles.m_mass[0]/' + pMC + ')**2)))'
102pseudoPZ0Fit = 'sqrt(' + ptFit + '**5/(' + pFit + '**3*(1+(MCParticles.m_mass[0]/' + pFit + ')**2)))'
103
104
105dBw = 50
106pBins = 50
107fBins = 100
108
109particleConditions = {
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
117Particles = ["Electron", "Muon", "Kaon", "Pion", "Proton"]
118
119unitImp = "mm"
120
121
122def 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 = f'{binWidth:8.3f}'
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
404DecZ = "MCParticles.m_decayVertex_z"
405
406if belleOrBelle2 == "Belle2":
407
408 variablesPlotB0tag = {'DecZ': [DecZ, dBw, -0.6, 1.5, r"{\rm " + unitImp + r"}\, "],
409 }
410
411else:
412 unitImp = "cm"
413 variablesPlotB0tag = {'DecZ': [DecZ, dBw, -1.5, 1.5, r"{\rm " + unitImp + r"}\, "],
414 }
415
416particleCondsB0tag = {
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
433def 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 = f'{binWidth:8.3f}'
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
602makePlotsForEachParticleKind("")
603
604if 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 && ")
608else:
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