Belle II Software development
PhokharaEvtgenPlot.py
1#!/usr/bin/env python3
2
3
10
11"""
12<header>
13 <input>PhokharaEvtgenAnalysis.root</input>
14 <contact>Kirill Chilikin (K.A.Chilikin@inp.nsk.su)</contact>
15 <description>Analysis of e+ e- -> J/psi eta_c events.</description>
16</header>
17"""
18
19import ROOT
20import math
21
22nbins_ratio = 91
23emin_ratio = 6.05
24emax_ratio = 10.6
25e_threshold = 6.077216
26vacpol_coef = [
27 1.0628848921570468,
28 1.0630222501435538,
29 1.0631583712375692,
30 1.0632932870752236,
31 1.0634270280964420,
32 1.0635596236029659,
33 1.0636911038563479,
34 1.0638214919564521,
35 1.0639508161489530,
36 1.0640911763911554,
37 1.0642426021238256,
38 1.0643930465044915,
39 1.0645425321383304,
40 1.0646910829012752,
41 1.0648387157899668,
42 1.0649854532304037,
43 1.0651313149251775,
44 1.0652763219764185,
45 1.0654204887222767,
46 1.0655638350494281,
47 1.0657063782276197,
48 1.0658481370304009,
49 1.0659891235597700,
50 1.0661293555636853,
51 1.0662688482593514,
52 1.0664076184528672,
53 1.0665456763532895,
54 1.0666830378959860,
55 1.0668197165557802,
56 1.0669728086606438,
57 1.0671423358423944,
58 1.0673112348835889,
59 1.0674795135759751,
60 1.0676471855363563,
61 1.0678142640605996,
62 1.0679807559683534,
63 1.0681466739443513,
64 1.0683120303864380,
65 1.0684768312444710,
66 1.0686410883688278,
67 1.0688048112946322,
68 1.0689680113655782,
69 1.0691306935088203,
70 1.0692928706527736,
71 1.0694545452046977,
72 1.0696157337719090,
73 1.0697764383327544,
74 1.0699366709637401,
75 1.0700964353013993,
76 1.0702369187285552,
77 1.0703581196112877,
78 1.0704788617065120,
79 1.0705991500549692,
80 1.0707189957188554,
81 1.0708384034084868,
82 1.0709573776751666,
83 1.0710759291100751,
84 1.0711940619641145,
85 1.0713117824086880,
86 1.0714290944104019,
87 1.0715460079976491,
88 1.0716625248052389,
89 1.0717786546050720,
90 1.0718944008474305,
91 1.0720097668616437,
92 1.0721247620626524,
93 1.0722393895512141,
94 1.0723536543855732,
95 1.0724675594488726,
96 1.0725811137265711,
97 1.0726943178288948,
98 1.0728071805444170,
99 1.0729197043593588,
100 1.0730318937365730,
101 1.0731437509784216,
102 1.0732552845114489,
103 1.0733664964658656,
104 1.0734773888878144,
105 1.0735879699573529,
106 1.0736982394896657,
107 1.0738082055112450,
108 1.0739178697580838,
109 1.0740272359656267,
110 1.0741363077987014,
111 1.0742450867795019,
112 1.0743535805843314,
113 1.0744617885288132,
114 1.0745697181614047,
115 1.0746773707456787,
116 1.0749101058390413,
117 1.0752680569983371]
118
119
120def ratio_measured_ratio(ecms, ecut):
121 alpha = 7.2973525664e-3
122 me = 0.510998928e-3
123 riemann_zeta_3 = 1.2020569032
124 pi = math.pi
125 l_e = math.log(ecms / 2 / ecut)
126 L = 2.0 * math.log(ecms / me)
127 r1 = -2.0 * l_e * (L - 1) + 1.5 * L + pi * pi / 3 - 2
128 r2 = 0.5 * pow(-2.0 * l_e * (L - 1), 2) + \
129 (1.5 * L + pi * pi / 3 - 2) * (-2.0 * l_e * (L - 1)) + L * L * (-l_e / 3 + 11. / 8 - pi * pi / 3) + \
130 L * (2.0 * l_e * l_e / 3 + 10. * l_e / 9 - 203. / 48 + 11. * pi * pi / 12 + 3.0 * riemann_zeta_3) - \
131 (4. * l_e * l_e * l_e / 9 + 10. * l_e * l_e / 9 + 2. * l_e / 9 * (28. / 3 - pi * pi)) - \
132 (pow(L - 2. * l_e, 3) / 18 - 5. / 18 * pow(L - 2. * l_e, 2) + (28. / 3 - pi * pi) * (L - 2. * l_e) / 9)
133 r = 1. + alpha / pi * r1 + alpha * alpha / pi / pi * r2
134 return r
135
136
137def cross_section_ee_mumu(ecms):
138 alpha = 7.2973525664e-3
139 mmu = 0.1056583715
140 conv = 0.389379338 * 1e12 # Gev^2 * fb
141 s = ecms * ecms
142 if (ecms < 2.0 * mmu):
143 return 0
144 return 4.0 * math.pi * alpha * alpha / (3.0 * s) * conv * \
145 math.sqrt(1.0 - 4.0 * mmu * mmu / s) * (1.0 + 2.0 * mmu * mmu / s)
146
147
148def effective_luminosity(ecms, e):
149 alpha = 7.2973525664e-3
150 me = 0.510998928e-3
151 pi = math.pi
152 L = 2.0 * math.log(ecms / me)
153 beta = 2.0 * alpha / pi * (L - 1)
154 x = 1.0 - pow(e / ecms, 2)
155 lum = beta * pow(x, beta - 1) * \
156 (1. + alpha / pi * (pi * pi / 3 - 1. / 2) + 3. * beta / 4 -
157 beta * beta / 24 * (L / 3 + 2. * pi * pi - 37. / 4)) - \
158 beta * (1. - x / 2) + beta * beta / 8 * \
159 (4. * (2. - x) * math.log(1. / x) + (1. + 3. * pow(1. - x, 2)) / x *
160 math.log(1. / (1. - x)) - 6 + x)
161 # Jacobian dx -> dE
162 lum = lum * 2.0 * e / ecms / ecms
163 return lum
164
165
166def effective_luminosity_integral(ecms, emin, emax):
167 if (emax < ecms):
168 npoints = 1000
169 intg = 0
170 for i in range(0, npoints):
171 e = emin + (emax - emin) / npoints * (float(i) + 0.5)
172 if (e > e_threshold):
173 intg = intg + effective_luminosity(ecms, e)
174 intg = intg / npoints
175 return intg
176 # The value of effective_luminosity -> infinity as ecut -> 0,
177 # but its integral converges => the approximate formula is used.
178 elif (emin < ecms):
179 ecut = ecms - emin
180 return ratio_measured_ratio(ecms, ecut) / (emax - emin)
181 else:
182 return 0
183
184
185input_file = ROOT.TFile('PhokharaEvtgenAnalysis.root')
186tree = input_file.Get('tree')
187output_file = ROOT.TFile('PhokharaEvtgen.root', 'recreate')
188
189h_mgamma_exp = ROOT.TH1F('h_mgamma_exp', 'Virtual #gamma mass distribution (theory)', nbins_ratio, emin_ratio, emax_ratio)
190
191h_ratio = ROOT.TH1F('h_ratio', 'Virtual #gamma mass distribution / expectation', nbins_ratio, emin_ratio, emax_ratio)
192h_ratio.GetXaxis().SetTitle('M_{J/#psi#eta_{c}}, GeV/c^{2}')
193h_ratio.GetYaxis().SetTitle('N / N_{expected}')
194
195h_helicity_gamma = ROOT.TH1F('h_helicity_gamma', 'Virtual #gamma helicity angle', 20, -1, 1)
196h_helicity_gamma.GetXaxis().SetTitle('cos#theta_{#gamma}')
197h_helicity_gamma.GetYaxis().SetTitle('Events / 0.1')
198h_helicity_gamma.SetMinimum(0)
199
200h_helicity_jpsi = ROOT.TH1F('h_helicity_jpsi', 'J/#psi helicity angle', 20, -1, 1)
201h_helicity_jpsi.GetXaxis().SetTitle('cos#theta_{J/#psi}')
202h_helicity_jpsi.GetYaxis().SetTitle('Events / 0.1')
203h_helicity_jpsi.SetMinimum(0)
204
205n = tree.GetEntries()
206tree.GetEntry(0)
207ecms = tree.ecms
208for i in range(0, n):
209 tree.GetEntry(i)
210 m = math.sqrt(tree.gamma_e * tree.gamma_e - tree.gamma_px * tree.gamma_px -
211 tree.gamma_py * tree.gamma_py - tree.gamma_pz * tree.gamma_pz)
212 h_ratio.Fill(m)
213 p_gamma = ROOT.TLorentzVector(tree.gamma_px, tree.gamma_py,
214 tree.gamma_pz, tree.gamma_e)
215 p_jpsi = ROOT.TLorentzVector(tree.jpsi_px, tree.jpsi_py,
216 tree.jpsi_pz, tree.jpsi_e)
217 p_lepton = ROOT.TLorentzVector(tree.lepton_px, tree.lepton_py,
218 tree.lepton_pz, tree.lepton_e)
219 v_boost_gamma = -p_gamma.BoostVector()
220 v_boost_jpsi = -p_jpsi.BoostVector()
221 p_jpsi.Boost(v_boost_gamma)
222 h_helicity_gamma.Fill(-math.cos(p_gamma.Vect().Angle(p_jpsi.Vect())))
223 p_lepton.Boost(v_boost_jpsi)
224 p_gamma.Boost(v_boost_jpsi)
225 h_helicity_jpsi.Fill(-math.cos(p_gamma.Vect().Angle(p_lepton.Vect())))
226
227for i in range(nbins_ratio, 0, -1):
228 emin = emin_ratio + (emax_ratio - emin_ratio) / nbins_ratio * (i - 1)
229 emax = emin_ratio + (emax_ratio - emin_ratio) / nbins_ratio * i
230 eff_lumi = effective_luminosity_integral(ecms, emin, emax)
231 xs = cross_section_ee_mumu(h_ratio.GetBinLowEdge(i)) * vacpol_coef[i - 1]
232 h_mgamma_exp.SetBinContent(i, xs * eff_lumi)
233
234h_mgamma_exp.Scale(n / h_mgamma_exp.Integral())
235
236for i in range(nbins_ratio, 0, -1):
237 val = h_ratio.GetBinContent(i)
238 err = h_ratio.GetBinError(i)
239 exp = h_mgamma_exp.GetBinContent(i)
240 h_ratio.SetBinContent(i, val / exp)
241 h_ratio.SetBinError(i, err / exp)
242
243contact = 'Kirill Chilikin (K.A.Chilikin@inp.nsk.su)'
244functions = h_ratio.GetListOfFunctions()
245functions.Add(ROOT.TNamed('Description', 'Number of events / theoretical expectation'))
246functions.Add(ROOT.TNamed('Check', 'Should be consistent with 1'))
247functions.Add(ROOT.TNamed('Contact', contact))
248functions.Add(ROOT.TNamed('MetaOptions', 'shifter'))
249functions = h_helicity_gamma.GetListOfFunctions()
250functions.Add(ROOT.TNamed('Description', 'Virtual photon helicity angle'))
251functions.Add(ROOT.TNamed('Check', 'Should be distributed as (1 + cos^2 theta)'))
252functions.Add(ROOT.TNamed('Contact', contact))
253functions.Add(ROOT.TNamed('MetaOptions', 'shifter'))
254functions = h_helicity_jpsi.GetListOfFunctions()
255functions.Add(ROOT.TNamed('Description', 'J/psi helicity angle'))
256functions.Add(ROOT.TNamed('Check', 'Should be distributed as (1 + cos^2 theta)'))
257functions.Add(ROOT.TNamed('Contact', contact))
258functions.Add(ROOT.TNamed('MetaOptions', 'shifter'))
259
260output_file.cd()
261h_ratio.Write()
262h_helicity_gamma.Write()
263h_helicity_jpsi.Write()
264output_file.Close()