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