Belle II Software  release-06-02-00
PhokharaEvtgenPlot.py
1 #!/usr/bin/env python3
2 # -*- coding: utf-8 -*-
3 
4 """
5 <header>
6  <input>PhokharaEvtgenAnalysis.root</input>
7  <contact>Kirill Chilikin (chilikin@lebedev.ru)</contact>
8  <description>Analysis of e+ e- -> J/psi eta_c events.</description>
9 </header>
10 """
11 
12 from basf2 import *
13 import ROOT
14 import math
15 
16 nbins_ratio = 91
17 emin_ratio = 6.05
18 emax_ratio = 10.6
19 e_threshold = 6.077216
20 vacpol_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 
114 def 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 
131 def 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 
142 def 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 
160 def 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 
179 input_file = ROOT.TFile('PhokharaEvtgenAnalysis.root')
180 tree = input_file.Get('tree')
181 output_file = ROOT.TFile('PhokharaEvtgen.root', 'recreate')
182 
183 h_mgamma_exp = ROOT.TH1F('h_mgamma_exp', 'Virtual #gamma mass distribution (theory)', nbins_ratio, emin_ratio, emax_ratio)
184 
185 h_ratio = ROOT.TH1F('h_ratio', 'Virtual #gamma mass distribution / expectation', nbins_ratio, emin_ratio, emax_ratio)
186 h_ratio.GetXaxis().SetTitle('M_{J/#psi#eta_{c}}, GeV/c^{2}')
187 h_ratio.GetYaxis().SetTitle('N / N_{expected}')
188 
189 h_helicity_gamma = ROOT.TH1F('h_helicity_gamma', 'Virtual #gamma helicity angle', 20, -1, 1)
190 h_helicity_gamma.GetXaxis().SetTitle('cos#theta_{#gamma}')
191 h_helicity_gamma.GetYaxis().SetTitle('Events / 0.1')
192 h_helicity_gamma.SetMinimum(0)
193 
194 h_helicity_jpsi = ROOT.TH1F('h_helicity_jpsi', 'J/#psi helicity angle', 20, -1, 1)
195 h_helicity_jpsi.GetXaxis().SetTitle('cos#theta_{J/#psi}')
196 h_helicity_jpsi.GetYaxis().SetTitle('Events / 0.1')
197 h_helicity_jpsi.SetMinimum(0)
198 
199 n = tree.GetEntries()
200 tree.GetEntry(0)
201 ecms = tree.ecms
202 for 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 
221 for 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 
228 h_mgamma_exp.Scale(n / h_mgamma_exp.Integral())
229 
230 for 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 
237 contact = 'Kirill Chilikin (chilikin@lebedev.ru)'
238 functions = h_ratio.GetListOfFunctions()
239 functions.Add(ROOT.TNamed('Description', 'Number of events / theoretical expectation'))
240 functions.Add(ROOT.TNamed('Check', 'Should be consistent with 1'))
241 functions.Add(ROOT.TNamed('Contact', contact))
242 functions.Add(ROOT.TNamed('MetaOptions', 'shifter'))
243 functions = h_helicity_gamma.GetListOfFunctions()
244 functions.Add(ROOT.TNamed('Description', 'Virtual photon helicity angle'))
245 functions.Add(ROOT.TNamed('Check', 'Should be distributed as (1 + cos^2 theta)'))
246 functions.Add(ROOT.TNamed('Contact', contact))
247 functions.Add(ROOT.TNamed('MetaOptions', 'shifter'))
248 functions = h_helicity_jpsi.GetListOfFunctions()
249 functions.Add(ROOT.TNamed('Description', 'J/psi helicity angle'))
250 functions.Add(ROOT.TNamed('Check', 'Should be distributed as (1 + cos^2 theta)'))
251 functions.Add(ROOT.TNamed('Contact', contact))
252 functions.Add(ROOT.TNamed('MetaOptions', 'shifter'))
253 
254 output_file.cd()
255 h_ratio.Write()
256 h_helicity_gamma.Write()
257 h_helicity_jpsi.Write()
258 output_file.Close()