Belle II Software  release-05-02-19
UtrepsB.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2019 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Kiyoshi Hayasaka, Yo Sato *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <TLorentzVector.h>
12 #include <generators/treps/UtrepsB.h>
13 
14 using namespace Belle2;
15 
16 
17 UtrepsB::UtrepsB(void) : TrepsB()
18 {
19 }
20 
21 void UtrepsB::initg(void)
22 {
23  B2DEBUG(20, " UtrepsB initg Pmodel " << TrepsB::pmodel);
24 
25  if (TrepsB::pmodel == 251) {
26  // pi+pi-
27 
28  //parametrization of pi+pi- partial waves
29 
30  static double data001[21] = { 20., 15., 9., 8., 8., 7., 5.3, 5.2, 5.2, 5.2, 5.2, 4., 3., 2., 2., 2., 2., 2., 1., 0., 0. };
31  static double data201[21] = { 1., 5., 5, 5., 4., 4., 4., 4., 3., 3., 3., 2., 1., 1., 1., 1., 0., 0., 0., 1., 1. };
32  static double data221[21] = { 1., 5., 5., 6., 7., 7.18, 6.68, 6.27, 5.91, 5.61, 5., 5., 4., 3., 2., 2., 2., 1., 0., 0., 0. };
33  b00 = Interps(0.3, 2.3, 20, data001);
34  b20 = Interps(0.3, 2.3, 20, data201);
35  b22 = Interps(0.3, 2.3, 20, data221);
36  }
37 
38  if (TrepsB::pmodel == 252) {
39  // K+K-
40  static double data002[15] = { 2., 2., 2., 2., 2., 2., 3., 4., 4.5, 5., 5., 5., 4.5, 4.5, 4.3 };
41  static double data202[15] = { 1., 2., 3., 4., 5., 5., 5., 5., 5., 5., 5., 4., 3.0, 2.2, 1. };
42  static double data222[15] = { 0., 0., 0., 0., 0., 0., -1., -3., -3., -3., -3., -4., -3.5, -3., -2.5 };
43  b00 = Interps(1.0, 2.4, 14, data002);
44  b20 = Interps(1.0, 2.4, 14, data202);
45  b22 = Interps(1.0, 2.4, 14, data222);
46  }
47 
48  if (TrepsB::pmodel == 253) {
49  // ppbar
50  static double data003[9] = { 2.3, 2.3, 7.26, 7.94, 8.53, 7.38, 3.25, 1.98, 2.30};
51  static double data203[9] = { 9.63, 9.63, 10.73, 8.02, 6.18, 3.37, 0.63, 0.10, 0.66};
52  static double data223[9] = { 1.48, 1.48, -4.62, -6.12, -6.78, -5.35, -1.82, -1.02, -1.63};
53  b00 = Interps(2.05, 2.85, 8, data003);
54  b20 = Interps(2.05, 2.85, 8, data203);
55  b22 = Interps(2.05, 2.85, 8, data223);
56  }
57 
58 
59  return;
60 }
61 
62 double UtrepsB::d2func(double y) const
63 {
64  return 1. / (9. + 3.*y * y + y * y * y * y);
65 }
66 
67 double UtrepsB::qfunc(double y, double mm) const
68 {
69  return sqrt(0.25 * y - mm * mm);
70 }
71 
72 //form factor effect
73 double UtrepsB::tpform(double _q2, double _w) const
74 {
75  double dis = 1. / pow((1.0 + _q2 / (_w * _w)), 2);
76  if (TrepsB::fmodel == 1) {
77  // rho form factor
78  dis = 1. / pow((1.0 + _q2 / (0.77 * 0.77)), 2);
79  }
80  if (TrepsB::fmodel == 2) {
81  // J/psi form factor
82  dis = 1. / pow((1.0 + _q2 / (3.097 * 3.097)), 2);
83  }
84 
85  return dis;
86 }
87 
88 
89 //angular distribution for final 2-body case
90 double UtrepsB::tpangd(double _z, double _w)
91 {
92  // no normalization is necessary
93 
94  double ww = _w;
95  double zz = _z;
96 
97  double dcs;
98 
99  dcs = 1.;
100 
101  if (TrepsB::pmodel == 201) {
102  // 2+(0) --> PP
103  dcs = 3.*zz * zz - 1.;
104  return dcs;
105  }
106 
107  if (TrepsB::pmodel == 202) {
108  // 2+(2) --> PP
109  dcs = pow(1. - zz * zz, 2);
110  return dcs;
111  }
112 
113  if (TrepsB::pmodel == 251) {
114  // pi+pi-
115  double mr = 1.2755;
116  double gtot = 0.1859;
117  double br = 0.561;
118  double sr = 3.62;
119  double ggg = 2.62e-6;
120  double phas = 22.8 / 180.*3.14159;
121  double gtot0 = 0.03;
122 
123  double v00 = 0., v20 = 0., v22 = 0. ;
124 
125  if (ww < 0.295) {
126  v00 = v20 = v22 = 0.0 ;
127  } else if (ww < 2.3) {
128  v00 = b00.get_val(ww);
129  v20 = b20.get_val(ww);
130  v22 = b22.get_val(ww);
131  }
132 
133  double dcs2 = 0.;
134 
135  if (ww < 2.3) {
136  double mr0 = 0.98;
137  double speak = 18.3;
138  double mm = 0.1396;
139  double gf = pow(qfunc(ww * ww, mm) / qfunc(mr * mr, mm), 5) *
140  d2func(qfunc(ww * ww, mm) * sr) / d2func(qfunc(mr * mr, mm) * sr);
141  std::complex<double> imagu(0.0, 1.0);
142  std::complex<double> rbr = sqrt(40.*3.14159 * mr / ww * ggg * gtot * br) * gf /
143  ((-ww * ww + mr * mr) - imagu * mr * gtot * gf) *
144  sqrt(389000.);
145  rbr = (cos(phas) + imagu * sin(phas)) * rbr;
146 
147  dcs2 = pow(v00 + v20 * sqrt(5.) * 0.5 * (3.*zz * zz - 1.), 2) +
148  (v22 * v22 + 2.*v22 * rbr.real() + pow(std::abs(rbr), 2)) *
149  15. / 8.*pow(1. - zz * zz, 2)
150  + speak * 0.25 * gtot0 * gtot0 / (pow(ww - mr0, 2) + 0.25 * gtot0 * gtot0);
151  } else {
152  if (abs(zz) < 0.8) dcs2 = 263. / pow(ww, 6) / pow(1. - zz * zz, 2);
153  else dcs2 = 263. / pow(ww, 6) / pow(1. - 0.8 * 0.8, 2);
154  }
155  return dcs2 * 0.5;
156  }
157 
158  if (TrepsB::pmodel == 252) {
159  //K+K-
160  double mr = 1.522;
161  double gtot = 0.0814;
162  double br = 0.45;
163  double sr = 3.62;
164  double ggg = 0.081e-6;
165  double gtotfa = 0.13;
166 
167  double s00 = 0., s20 = 0., s22 = 0. ;
168 
169  if (ww < 1.1) {
170  s00 = s20 = s22 = 0.0 ;
171  } else if (ww < 2.4) {
172  s00 = b00.get_val(ww);
173  s20 = b20.get_val(ww);
174  s22 = b22.get_val(ww);
175  }
176 
177  double dcs2 = 0.;
178 
179  if (ww < 2.4) {
180  double mrfa = 1.29;
181  double speak = 40.;
182  double mm = 0.4937;
183  double gf = pow(qfunc(ww * ww, mm) / qfunc(mr * mr, mm), 5) *
184  d2func(qfunc(ww * ww, mm) * sr) / d2func(qfunc(mr * mr, mm) * sr);
185  std::complex<double> imagu(0.0, 1.0);
186  std::complex<double> rbr = sqrt(40.*3.14159 * mr / ww * ggg * gtot * br) * gf /
187  ((-ww * ww + mr * mr) - imagu * mr * gtot * gf) *
188  sqrt(389000.);
189 
190  dcs2 = s00 + s20 * pow(sqrt(5.) * 0.5 * (3.*zz * zz - 1.), 2) +
191  (s22 + pow(std::abs(rbr), 2)
192  + speak * 0.25 * gtotfa * gtotfa / (pow(ww - mrfa, 2) + 0.25 * gtotfa * gtotfa)) *
193  15. / 8.*pow(1. - zz * zz, 2) ;
194  } else {
195  if (abs(zz) < 0.8) dcs2 = 241. / pow(ww, 6) / pow(1. - zz * zz, 2);
196  else dcs2 = 241. / pow(ww, 6) / pow(1. - 0.8 * 0.8, 2);
197  }
198  return dcs2 * 0.5;
199  }
200 
201  if (TrepsB::pmodel == 253) {
202  //ppbar
203  double s00 = 0., s20 = 0., s22 = 0. ;
204 
205  if (ww < 2.85) {
206  s00 = b00.get_val(ww);
207  s20 = b20.get_val(ww);
208  s22 = b22.get_val(ww);
209  } else {
210  double ww0 = 2.85;
211  s00 = b00.get_val(ww0) * pow(ww0 / ww, 12);
212  s20 = b20.get_val(ww0) * pow(ww0 / ww, 12);
213  s22 = b22.get_val(ww0) * pow(ww0 / ww, 12);
214  }
215  double dcs2 = 0.;
216 
217  // eta_c resonance
218  double mr = 2.984;
219  double gtot = 0.032;
220  double gggbr ; gggbr = 0.00152 * 0.000005;
221 
222  double sigr = 8.*3.14159 * mr / ww * gggbr * gtot / (pow(mr * mr - ww * ww, 2) + mr * mr * gtot * gtot)
223  * 389000.;
224 
225  if (ww > 2.0)
226  dcs2 = s00 + sigr + s20 * (5. / 4.) * pow(3.*zz * zz - 1, 2) + s22 * (15. / 8.) * pow(1. - zz * zz, 2);
227 
228  return dcs2 * 0.5;
229  }
230 
231 
232  return dcs;
233 
234 }
235 
236 int UtrepsB::tpuser(TLorentzVector _pe, TLorentzVector _pp,
237  Part_gen* part, int _npart)
238 {
239  // user decision routine for extra generation conditions.
240  // Return positive integer for the generation, otherwise, this event will
241  // be canceled.
242  // CAUTION!: The 4-momenta of particles are in e+e- cms
243  B2DEBUG(20, " _pe(" << _pe.Px() << "," << _pe.Py() << "," << _pe.Pz() << ") and "
244  << " _pp(" << _pp.Px() << "," << _pp.Py() << "," << _pp.Pz() << ") are given but not used");
245 
246 
247  int iret = 1;
248 
249  // 3-body physics models
250 
251  if (_npart == 3 && TrepsB::pmodel >= 301 && TrepsB::pmodel <= 399) {
252 
253  int index1 = 0, index2 = 1, index3 = 2;
254  double z, m12, zp, phip, zs, phis, phi0;
255 
256  tpkin3(part, index1, index2, index3, z, m12, zp,
257  phip, zs, phis, phi0);
258 
259  double u = sqrt(1. - z * z);
260  double up = sqrt(1. - zp * zp);
261  double us = sqrt(1. - zs * zs);
262 
263  double wei = 0., weimax = 0.;
264 
265 
266 
267  if (TrepsB::pmodel == 301) {
268  // 0- -> VP -> 3P
269  wei = zs * zs;
270  weimax = 1.;
271  }
272  if (TrepsB::pmodel == 302) {
273  // 2+(0) -> VP -> 3P
274  wei = z * z * u * u * up * up * pow(sin(phip), 2);
275  weimax = 1.;
276  }
277  if (TrepsB::pmodel == 303) {
278  // 2+(2) -> VP -> 3P
279  wei = u * u * (u * u * zp * zp + z * z * up * up - 2.*u * z * up * zp * cos(phip));
280  weimax = 2.5;
281  }
282  if (TrepsB::pmodel == 304) {
283  // 0- -> TP -> 3P
284  wei = pow(3.*zs * zs - 1., 2);
285  weimax = 4.;
286  }
287  if (TrepsB::pmodel == 305) {
288  // 2+(0) -> TP -> 3P
289  wei = u * u * up * up * zp * zp * pow(sin(phip), 2);
290  weimax = 1.;
291  }
292  if (TrepsB::pmodel == 306) {
293  // 2+(2) -> TP -> 3P
294  wei = up * up * (z * z * up * up + u * u * zp * zp - 2.*z * u * zp * up * cos(phip));
295  weimax = 2.5;
296  }
297  if (TrepsB::pmodel == 307) {
298  // 2- -> SP -> 3P
299  wei = pow(3.*zp * zp - 1., 2);
300  weimax = 2.;
301  }
302  if (TrepsB::pmodel == 308) {
303  // 2- -> VP -> 3P
304  wei = z * zp - 0.5 * u * up * cos(phip);
305  weimax = 1.5;
306  }
307  if (TrepsB::pmodel == 309) {
308  // 2- -> TP -> 3P
309  wei = pow(3.*zs * zs - 1., 2);
310  weimax = 4.;
311  }
312 
313  if (TrepsB::pmodel == 331) {
314  // 2+(0) -> VP -> PgP
315  wei = u * u * (1. - us * us * pow(cos(phis), 2));
316  weimax = 1.;
317  }
318  if (TrepsB::pmodel == 332) {
319  // 2+(2) -> VP -> PgP
320  wei = u * u * (1. + z * z * zs * zs - u * u * us * us * pow(cos(phis), 2));
321  weimax = 2.;
322  }
323  if (TrepsB::pmodel == 333) {
324  // 0+/- -> Vg -> PPg
325  wei = us * us;
326  weimax = 1.;
327  }
328  if (TrepsB::pmodel == 334) {
329  // 0+ -> Vg (E1) --> llg
330  wei = 1. + zs * zs;
331  weimax = 2.;
332  }
333  if (TrepsB::pmodel == 335) {
334  // 2+(0) -> Vg (E1) --> llg
335  wei = (9.*z * z * z * z - 12.*z * z + 5.) / 24.*(1. + zs * zs)
336  + 0.75 * z * z * (1. - z * z) * (1. - zs * zs) + 0.5 * z * u * zs * us * (-3.*z * z + 2.) * cos(phis)
337  + 0.125 * (1. - z * z) * (3.*z * z - 1.) * (1. - zs * zs) * cos(2.*phis);
338  weimax = 13. / 12. + 0.75 + 1.0 + 0.25;
339  }
340  if (TrepsB::pmodel == 336) {
341  // 2+(2) -> Vg (E1) --> llg
342  wei = 0.125 * pow(1. + z * z, 2) * (1. + zs * zs)
343  + 0.25 * (1. - z * z * z * z) * (1. - zs * zs) - 0.5 * z * u * zs * us * (1. + z * z) * cos(phis)
344  + 0.125 * (1. - z * z * z * z) * (1. - zs * zs) * cos(2.*phis);
345  weimax = 13. / 8.;
346  }
347  if (weimax * gRandom->Uniform() > wei) iret = 0;
348 
349  B2DEBUG(20, " $$ 3B $$ " << wei << " " << iret);
350  }
351 
352 
353  if (_npart == 4 && TrepsB::pmodel >= 401 && TrepsB::pmodel <= 499) {
354 
355  int index1 = 0, index2 = 1, index3 = 2, index4 = 3;
356  double z, m12, zp, phip, zs, phis, m34, zpp, phipp, zss, phiss, phi0;
357 
358  tpkin4(part,
359  index1, index2, index3, index4, z, m12, zp,
360  phip, zs, phis, m34, zpp, phipp, zss, phiss, phi0) ;
361 
362  }
363 
364  return iret ;
365 }
Belle2::Interps
Definition: Sutool.h:58
Belle2::Part_gen
Definition: Particle_array.h:50
Belle2::TrepsB
Definition: Treps3B.h:41
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
prepareAsicCrosstalkSimDB.u
u
merged u1 and u2
Definition: prepareAsicCrosstalkSimDB.py:46
Belle2::UtrepsB::initg
void initg()
initialization of Pparametrization of pi+pi- partial waves
Definition: UtrepsB.cc:21