11 #include <TLorentzVector.h>
12 #include <generators/treps/UtrepsB.h>
17 UtrepsB::UtrepsB(
void) :
TrepsB()
23 B2DEBUG(20,
" UtrepsB initg Pmodel " << TrepsB::pmodel);
25 if (TrepsB::pmodel == 251) {
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);
38 if (TrepsB::pmodel == 252) {
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);
48 if (TrepsB::pmodel == 253) {
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);
62 double UtrepsB::d2func(
double y)
const
64 return 1. / (9. + 3.*y * y + y * y * y * y);
67 double UtrepsB::qfunc(
double y,
double mm)
const
69 return sqrt(0.25 * y - mm * mm);
73 double UtrepsB::tpform(
double _q2,
double _w)
const
75 double dis = 1. / pow((1.0 + _q2 / (_w * _w)), 2);
76 if (TrepsB::fmodel == 1) {
78 dis = 1. / pow((1.0 + _q2 / (0.77 * 0.77)), 2);
80 if (TrepsB::fmodel == 2) {
82 dis = 1. / pow((1.0 + _q2 / (3.097 * 3.097)), 2);
90 double UtrepsB::tpangd(
double _z,
double _w)
101 if (TrepsB::pmodel == 201) {
103 dcs = 3.*zz * zz - 1.;
107 if (TrepsB::pmodel == 202) {
109 dcs = pow(1. - zz * zz, 2);
113 if (TrepsB::pmodel == 251) {
116 double gtot = 0.1859;
119 double ggg = 2.62e-6;
120 double phas = 22.8 / 180.*3.14159;
123 double v00 = 0., v20 = 0., v22 = 0. ;
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);
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) *
145 rbr = (cos(phas) + imagu * sin(phas)) * rbr;
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);
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);
158 if (TrepsB::pmodel == 252) {
161 double gtot = 0.0814;
164 double ggg = 0.081e-6;
165 double gtotfa = 0.13;
167 double s00 = 0., s20 = 0., s22 = 0. ;
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);
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) *
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) ;
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);
201 if (TrepsB::pmodel == 253) {
203 double s00 = 0., s20 = 0., s22 = 0. ;
206 s00 = b00.get_val(ww);
207 s20 = b20.get_val(ww);
208 s22 = b22.get_val(ww);
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);
220 double gggbr ; gggbr = 0.00152 * 0.000005;
222 double sigr = 8.*3.14159 * mr / ww * gggbr * gtot / (pow(mr * mr - ww * ww, 2) + mr * mr * gtot * gtot)
226 dcs2 = s00 + sigr + s20 * (5. / 4.) * pow(3.*zz * zz - 1, 2) + s22 * (15. / 8.) * pow(1. - zz * zz, 2);
236 int UtrepsB::tpuser(TLorentzVector _pe, TLorentzVector _pp,
243 B2DEBUG(20,
" _pe(" << _pe.Px() <<
"," << _pe.Py() <<
"," << _pe.Pz() <<
") and "
244 <<
" _pp(" << _pp.Px() <<
"," << _pp.Py() <<
"," << _pp.Pz() <<
") are given but not used");
251 if (_npart == 3 && TrepsB::pmodel >= 301 && TrepsB::pmodel <= 399) {
253 int index1 = 0, index2 = 1, index3 = 2;
254 double z, m12, zp, phip, zs, phis, phi0;
256 tpkin3(part, index1, index2, index3, z, m12, zp,
257 phip, zs, phis, phi0);
259 double u = sqrt(1. - z * z);
260 double up = sqrt(1. - zp * zp);
261 double us = sqrt(1. - zs * zs);
263 double wei = 0., weimax = 0.;
267 if (TrepsB::pmodel == 301) {
272 if (TrepsB::pmodel == 302) {
274 wei = z * z *
u *
u * up * up * pow(sin(phip), 2);
277 if (TrepsB::pmodel == 303) {
279 wei =
u *
u * (
u *
u * zp * zp + z * z * up * up - 2.*
u * z * up * zp * cos(phip));
282 if (TrepsB::pmodel == 304) {
284 wei = pow(3.*zs * zs - 1., 2);
287 if (TrepsB::pmodel == 305) {
289 wei =
u *
u * up * up * zp * zp * pow(sin(phip), 2);
292 if (TrepsB::pmodel == 306) {
294 wei = up * up * (z * z * up * up +
u *
u * zp * zp - 2.*z *
u * zp * up * cos(phip));
297 if (TrepsB::pmodel == 307) {
299 wei = pow(3.*zp * zp - 1., 2);
302 if (TrepsB::pmodel == 308) {
304 wei = z * zp - 0.5 *
u * up * cos(phip);
307 if (TrepsB::pmodel == 309) {
309 wei = pow(3.*zs * zs - 1., 2);
313 if (TrepsB::pmodel == 331) {
315 wei =
u *
u * (1. - us * us * pow(cos(phis), 2));
318 if (TrepsB::pmodel == 332) {
320 wei =
u *
u * (1. + z * z * zs * zs -
u *
u * us * us * pow(cos(phis), 2));
323 if (TrepsB::pmodel == 333) {
328 if (TrepsB::pmodel == 334) {
333 if (TrepsB::pmodel == 335) {
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;
340 if (TrepsB::pmodel == 336) {
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);
347 if (weimax * gRandom->Uniform() > wei) iret = 0;
349 B2DEBUG(20,
" $$ 3B $$ " << wei <<
" " << iret);
353 if (_npart == 4 && TrepsB::pmodel >= 401 && TrepsB::pmodel <= 499) {
355 int index1 = 0, index2 = 1, index3 = 2, index4 = 3;
356 double z, m12, zp, phip, zs, phis, m34, zpp, phipp, zss, phiss, phi0;
359 index1, index2, index3, index4, z, m12, zp,
360 phip, zs, phis, m34, zpp, phipp, zss, phiss, phi0) ;