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