Belle II Software development
EvtD0Topippim2pi0.cc
1// Model: EvtD0Topippim2pi0
2// This file is an amplitude model for D0 -> pi+ pi- pi0 pi0.
3// The model is from the BESIII Collaboration in Chin. Phys. C 48, 083001 (2024). DOI:  https://doi.org/10.1088/1674-1137/ad3d4d
4//
5// Permission to include these files in basf2 was generously granted by the BESIII Collaboration.
6//
7// Please cite the original reference for any public/published results where this model was used.
8
9#include <EvtGenBase/EvtPatches.hh>
10#include <EvtGenBase/EvtParticle.hh>
11#include <EvtGenBase/EvtGenKine.hh>
12#include <EvtGenBase/EvtPDL.hh>
13#include <EvtGenBase/EvtVector4R.hh>
14#include <EvtGenBase/EvtVector3R.hh>
15#include <EvtGenBase/EvtReport.hh>
16#include <stdlib.h>
17#include <iostream>
18#include <string>
19#include <complex>
20#include <vector>
21#include <math.h>
22#include <TMath.h>
23
24#include <generators/evtgen/EvtGenModelRegister.h>
25#include <generators/evtgen/models/besiii/EvtD0Topippim2pi0.h>
26
27namespace Belle2 {
32
34 B2_EVTGEN_REGISTER_MODEL(EvtD0Topippim2pi0);
35
36 EvtD0Topippim2pi0::~EvtD0Topippim2pi0() {}
37
38 std::string EvtD0Topippim2pi0::getName()
39 {
40 return "D0Topippim2pi0";
41 }
42
43 EvtDecayBase* EvtD0Topippim2pi0::clone()
44 {
45 return new EvtD0Topippim2pi0;
46 }
47
48 void EvtD0Topippim2pi0::init()
49 {
50 checkNArg(2);
51 checkNDaug(4);
52 charm = getArg(0);
53 tagmode = getArg(1);
54
55 double mag[36], pha[36];
56 mag[0] = 100; pha[0] = 0;
57 mag[1] = 6.03691; pha[1] = 0.0111121;
58 mag[2] = 31.0556; pha[2] = -1.58228;
59 mag[3] = 82.5439; pha[3] = -2.17864;
60 mag[4] = 246.483; pha[4] = 0.337741;
61 mag[5] = 659.439; pha[5] = 1.76185;
62 mag[6] = 0.353305; pha[6] = 0.22473;
63 mag[7] = 0.508529; pha[7] = -2.99122;
64 mag[8] = 19.0362; pha[8] = 2.70136;
65 mag[9] = 20.1349; pha[9] = -2.06863;
66 mag[10] = 10.5129; pha[10] = -1.26235;
67 mag[11] = 0.23207; pha[11] = -2.91839;
68 mag[12] = 0.295583; pha[12] = -0.47002;
69 mag[13] = 9.71285; pha[13] = -0.58869;
70 mag[14] = 65.1187; pha[14] = -2.63145;
71 mag[15] = 5.04613; pha[15] = -0.642076;
72 mag[16] = 43.8422; pha[16] = 0.339301;
73 mag[17] = 76.3005; pha[17] = -2.32496;
74 mag[18] = 61.1311; pha[18] = 0.609366;
75 mag[19] = 12.2241; pha[19] = -1.12858;
76 mag[20] = 13.6557; pha[20] = 3.02972;
77 mag[21] = 7.09973; pha[21] = -1.69019;
78 mag[22] = 4.5858; pha[22] = 0.0637721;
79 mag[23] = 8.0728; pha[23] = -1.01323;
80 mag[24] = 8.40155; pha[24] = -1.68028;
81 mag[25] = 40.748; pha[25] = -0.503741;
82 mag[26] = 120.464; pha[26] = 1.72667;
83 mag[27] = 2224.09; pha[27] = -1.04373;
84 mag[28] = 7286.76; pha[28] = 1.72657;
85 mag[29] = 8815.51; pha[29] = -1.10724;
86 mag[30] = 2433.26; pha[30] = 1.79639;
87 mag[31] = 5417.04; pha[31] = 2.67794;
88 mag[32] = 18.3282; pha[32] = -1.3945;
89 mag[33] = 55.535; pha[33] = 2.29339;
90 mag[34] = 1.57835; pha[34] = -0.497796;
91 mag[35] = 0.439629; pha[35] = 2.50596;
92
93 for (int i = 0; i < 36; i++) {
94 std::complex<double> ctemp(mag[i]*cos(pha[i]), mag[i]*sin(pha[i]));
95 fitpara.push_back(ctemp);
96 }
97
98 for (int i = 0; i < 4; i++) {
99 for (int j = 0; j < 4; j++) {
100 if (i != j) {
101 g_uv.push_back(0.0);
102 } else if (i < 3) {
103 g_uv.push_back(-1.0);
104 } else if (i == 3) {
105 g_uv.push_back(1.0);
106 }
107 }
108 }
109
110 for (int i = 0; i < 4; i++) {
111 for (int j = 0; j < 4; j++) {
112 for (int k = 0; k < 4; k++) {
113 for (int l = 0; l < 4; l++) {
114 if (i == j || i == k || i == l || j == k || j == l || k == l) {
115 epsilon_uvmn.push_back(0.0);
116 } else {
117 if (i == 0 && j == 1 && k == 2 && l == 3) epsilon_uvmn.push_back(1.0);
118 if (i == 0 && j == 1 && k == 3 && l == 2) epsilon_uvmn.push_back(-1.0);
119 if (i == 0 && j == 2 && k == 1 && l == 3) epsilon_uvmn.push_back(-1.0);
120 if (i == 0 && j == 2 && k == 3 && l == 1) epsilon_uvmn.push_back(1.0);
121 if (i == 0 && j == 3 && k == 1 && l == 2) epsilon_uvmn.push_back(1.0);
122 if (i == 0 && j == 3 && k == 2 && l == 1) epsilon_uvmn.push_back(-1.0);
123
124 if (i == 1 && j == 0 && k == 2 && l == 3) epsilon_uvmn.push_back(-1.0);
125 if (i == 1 && j == 0 && k == 3 && l == 2) epsilon_uvmn.push_back(1.0);
126 if (i == 1 && j == 2 && k == 0 && l == 3) epsilon_uvmn.push_back(1.0);
127 if (i == 1 && j == 2 && k == 3 && l == 0) epsilon_uvmn.push_back(-1.0);
128 if (i == 1 && j == 3 && k == 0 && l == 2) epsilon_uvmn.push_back(-1.0);
129 if (i == 1 && j == 3 && k == 2 && l == 0) epsilon_uvmn.push_back(1.0);
130
131 if (i == 2 && j == 0 && k == 1 && l == 3) epsilon_uvmn.push_back(1.0);
132 if (i == 2 && j == 0 && k == 3 && l == 1) epsilon_uvmn.push_back(-1.0);
133 if (i == 2 && j == 1 && k == 0 && l == 3) epsilon_uvmn.push_back(-1.0);
134 if (i == 2 && j == 1 && k == 3 && l == 0) epsilon_uvmn.push_back(1.0);
135 if (i == 2 && j == 3 && k == 0 && l == 1) epsilon_uvmn.push_back(1.0);
136 if (i == 2 && j == 3 && k == 1 && l == 0) epsilon_uvmn.push_back(-1.0);
137
138 if (i == 3 && j == 0 && k == 1 && l == 2) epsilon_uvmn.push_back(-1.0);
139 if (i == 3 && j == 0 && k == 2 && l == 1) epsilon_uvmn.push_back(1.0);
140 if (i == 3 && j == 1 && k == 0 && l == 2) epsilon_uvmn.push_back(1.0);
141 if (i == 3 && j == 1 && k == 2 && l == 0) epsilon_uvmn.push_back(-1.0);
142 if (i == 3 && j == 2 && k == 0 && l == 1) epsilon_uvmn.push_back(-1.0);
143 if (i == 3 && j == 2 && k == 1 && l == 0) epsilon_uvmn.push_back(1.0);
144
145 }
146 }
147 }
148 }
149 }
150
151 _nd = 4;
152 math_pi = 3.1415926f;
153 mass_Pion = 0.13957f;
154
155 rRes = 3.0 * 0.197321;
156 rD = 5.0 * 0.197321;
157 m_Pi = 0.139570;
158 m_Pi0 = 0.134977;
159 m2_Pi = m_Pi * m_Pi;
160 m2_Pi0 = m_Pi0 * m_Pi0;
161 m_Ka = 0.493677;
162 m2_Ka = m_Ka * m_Ka;
163
164 m0_f0980 = 0.965;
165 g1_f0980 = 0.165;
166 g2_f0980 = 4.210;
167
168 m0_rho7700 = 0.77526;
169 w0_rho7700 = 0.1478;
170
171 m0_rho770p = 0.77511;
172 w0_rho770p = 0.1491;
173
174 m0_rho1450 = 1.465;
175 w0_rho1450 = 0.400;
176
177 m0_f21270 = 1.2755;
178 w0_f21270 = 0.1867;
179
180 m0_a11260 = 1.1927;
181 g1_a11260 = 0.003777;
182 g2_a11260 = 0.0;
183
184 m0_pi1300 = 1.534;
185 w0_pi1300 = 0.610;
186
187 m0_a11420 = 1.411;
188 w0_a11420 = 0.161;
189
190 m0_a11640 = 1.655;
191 w0_a11640 = 0.254;
192
193 m0_a21320 = 1.3186;
194 w0_a21320 = 0.105;
195
196 m0_pi21670 = 1.6706;
197 w0_pi21670 = 0.258;
198
199 m0_h11170 = 1.166;
200 w0_h11170 = 0.375;
201
202 m0_omega = 0.78265;
203 w0_omega = 0.00849;
204
205 m0_phi = 1.019461;
206 w0_phi = 0.004249;
207
208 s0_prod = -5.0;
209
210 }
211
212 void EvtD0Topippim2pi0::initProbMax()
213 {
214 setProbMax(1098009205);
215 }
216
217 void EvtD0Topippim2pi0::decay(EvtParticle* p)
218 {
219 /*
220 double maxprob = 0.0;
221 for(int ir=0;ir<=60000000;ir++){
222 p->initializePhaseSpace(getNDaug(),getDaugs());
223 for(int i=0; i<_nd; i++){
224 _p4Lab[i]=p->getDaug(i)->getP4Lab();
225 _p4CM[i]=p->getDaug(i)->getP4();
226 }
227 double Prob = AmplitudeSquare(charm, tagmode);
228 if(Prob>maxprob) {
229 maxprob=Prob;
230 }
231 }
232 */
233 p->initializePhaseSpace(getNDaug(), getDaugs());
234 for (int i = 0; i < _nd; i++) {
235 _p4Lab[i] = p->getDaug(i)->getP4Lab();
236 _p4CM[i] = p->getDaug(i)->getP4();
237 }
238 double prob = AmplitudeSquare(charm, tagmode);
239 setProb(prob);
240 return;
241 }
242
243 void EvtD0Topippim2pi0::setInput(double* pip, double* pim, double* pi01, double* pi02)
244 {
245 m_Pip.push_back(pip[0]); m_Pim.push_back(pim[0]); m_Pi01.push_back(pi01[0]); m_Pi02.push_back(pi02[0]);
246 m_Pip.push_back(pip[1]); m_Pim.push_back(pim[1]); m_Pi01.push_back(pi01[1]); m_Pi02.push_back(pi02[1]);
247 m_Pip.push_back(pip[2]); m_Pim.push_back(pim[2]); m_Pi01.push_back(pi01[2]); m_Pi02.push_back(pi02[2]);
248 m_Pip.push_back(pip[3]); m_Pim.push_back(pim[3]); m_Pi01.push_back(pi01[3]); m_Pi02.push_back(pi02[3]);
249 }
250
251 std::vector<double> EvtD0Topippim2pi0::sum_tensor(std::vector<double> pa, std::vector<double> pb)
252 {
253
254 std::vector<double> temp;
255 for (size_t i = 0; i < pa.size(); i++) {
256 double sum = pa[i] + pb[i];
257 temp.push_back(sum);
258 }
259 return temp;
260 }
261
262 double EvtD0Topippim2pi0::contract_11_0(std::vector<double> pa, std::vector<double> pb)
263 {
264 double temp = pa[3] * pb[3] - pa[0] * pb[0] - pa[1] * pb[1] - pa[2] * pb[2];
265 return temp;
266
267 }
268
269 std::vector<double> EvtD0Topippim2pi0::contract_21_1(std::vector<double> pa, std::vector<double> pb)
270 {
271 std::vector<double> temp;
272 for (int i = 0; i < 4; i++) {
273 double sum = 0;
274 for (int j = 0; j < 4; j++) {
275 int idx = i * 4 + j;
276 sum += pa[idx] * pb[j] * g_uv[4 * j + j];
277 }
278 temp.push_back(sum);
279 }
280 return temp;
281
282 }
283
284 double EvtD0Topippim2pi0::contract_22_0(std::vector<double> pa, std::vector<double> pb)
285 {
286 double temp = 0;
287 for (int i = 0; i < 4; i++) {
288 for (int j = 0; j < 4; j++) {
289 int idx = i * 4 + j;
290 temp += pa[idx] * pb[idx] * g_uv[4 * i + i] * g_uv[4 * j + j];
291 }
292 }
293 return temp;
294
295 }
296
297 std::vector<double> EvtD0Topippim2pi0::contract_31_2(std::vector<double> pa, std::vector<double> pb)
298 {
299 std::vector<double> temp;
300 for (int i = 0; i < 16; i++) {
301 double sum = 0;
302 for (int j = 0; j < 4; j++) {
303 int idx = i * 4 + j;
304 sum += pa[idx] * pb[j] * g_uv[4 * j + j];
305 }
306 temp.push_back(sum);
307 }
308 return temp;
309
310 }
311
312 std::vector<double> EvtD0Topippim2pi0::contract_41_3(std::vector<double> pa, std::vector<double> pb)
313 {
314 std::vector<double> temp;
315 for (int i = 0; i < 64; i++) {
316 double sum = 0;
317 for (int j = 0; j < 4; j++) {
318 int idx = i * 4 + j;
319 sum += pa[idx] * pb[j] * g_uv[4 * j + j];
320 }
321 temp.push_back(sum);
322 }
323 return temp;
324
325 }
326
327 std::vector<double> EvtD0Topippim2pi0::contract_42_2(std::vector<double> pa, std::vector<double> pb)
328 {
329 std::vector<double> temp;
330 for (int i = 0; i < 16; i++) {
331 double sum = 0;
332 for (int j = 0; j < 4; j++) {
333 for (int k = 0; k < 4; k++) {
334 int idxa = i * 16 + j * 4 + k;
335 int idxb = j * 4 + k;
336 sum += pa[idxa] * pb[idxb] * g_uv[4 * j + j] * g_uv[4 * k + k];
337 }
338 }
339 temp.push_back(sum);
340 }
341
342 return temp;
343
344 }
345 std::vector<double> EvtD0Topippim2pi0::contract_22_2(std::vector<double> pa, std::vector<double> pb)
346 {
347 std::vector<double> temp;
348 for (int i = 0; i < 4; i++) {
349 for (int j = 0; j < 4; j++) {
350 double sum = 0;
351 for (int k = 0; k < 4; k++) {
352 int idxa = i * 4 + k;
353 int idxb = j * 4 + k;
354 sum += pa[idxa] * pb[idxb] * g_uv[4 * k + k];
355 }
356 temp.push_back(sum);
357 }
358 }
359
360 return temp;
361
362 }
363
364//OrbitalTensors
365 std::vector<double> EvtD0Topippim2pi0::OrbitalTensors(std::vector<double> pa, std::vector<double> pb, std::vector<double> pc,
366 double r, int rank)
367 {
368 // relative momentum
369 std::vector<double> mr;
370
371 for (int i = 0; i < 4; i++) {
372 double temp = pb[i] - pc[i];
373 mr.push_back(temp);
374 }
375
376 // "Masses" of particles
377 double msa = contract_11_0(pa, pa);
378 double msb = contract_11_0(pb, pb);
379 double msc = contract_11_0(pc, pc);
380
381 // Q^2_abc
382 double top = msa + msb - msc;
383 double Q2abc = top * top / (4.0 * msa) - msb;
384
385 // Barrier factors
386 double Q_0 = 0.197321f / r;
387 double Q_02 = Q_0 * Q_0;
388 double Q_04 = Q_02 * Q_02;
389 // double Q_06 = Q_04*Q_02;
390 // double Q_08 = Q_04*Q_04;
391
392 double Q4abc = Q2abc * Q2abc;
393 // double Q6abc = Q4abc*Q2abc;
394 // double Q8abc = Q4abc*Q4abc;
395
396 double mB1 = sqrt(2.0f / (Q2abc + Q_02));
397 double mB2 = sqrt(13.0f / (Q4abc + (3.0f * Q_02) * Q2abc + 9.0f * Q_04));
398 // mB3 = &sqrt(277.0f/(Q6abc + (6.0f*Q_02)*Q4abc + (45.0f*Q_04)*Q2abc + 225.0f*Q_06));
399 // mB4 = &sqrt(12746.0f/(Q8abc + (10.0f*Q_02)*Q6abc + (135.0f*Q_04)*Q4abc + (1575.0f*Q_06)*Q2abc + 11025.0f*Q_08));
400
401 // Projection Operator 2-rank
402 std::vector<double> proj_uv;
403 for (int i = 0; i < 4; i++) {
404 for (int j = 0; j < 4; j++) {
405 int idx = i * 4 + j;
406 double temp = -g_uv[idx] + pa[i] * pa[j] / msa;
407 proj_uv.push_back(temp);
408 }
409 }
410
411 // Orbital Tensors
412 if (rank == 0) {
413 std::vector<double> t;
414 t.push_back(1.0);
415 return t;
416
417 } else if (rank < 3) {
418 std::vector<double> t_u;
419 std::vector<double> Bt_u;
420 for (int i = 0; i < 4; i++) {
421 double temp = 0;
422 for (int j = 0; j < 4; j++) {
423 int idx = i * 4 + j;
424 temp += -proj_uv[idx] * mr[j] * g_uv[j * 4 + j];
425 }
426 t_u.push_back(temp);
427 Bt_u.push_back(temp * mB1);
428 }
429 if (rank == 1) return Bt_u;
430
431 double t_u2 = contract_11_0(t_u, t_u);
432
433 std::vector<double> Bt_uv;
434 for (int i = 0; i < 4; i++) {
435 for (int j = 0; j < 4; j++) {
436 int idx = 4 * i + j;
437 double temp = t_u[i] * t_u[j] + (1.0 / 3.0) * proj_uv[idx] * t_u2;
438 Bt_uv.push_back(temp * mB2);
439 }
440 }
441 if (rank == 2) return Bt_uv;
442
443 }
444 return std::vector<double>();
445 }
446
447// projection Tensor
448 std::vector<double> EvtD0Topippim2pi0::ProjectionTensors(std::vector<double> pa, int rank)
449 {
450 double msa = contract_11_0(pa, pa);
451
452 // Projection Operator 2-rank
453 std::vector<double> proj_uv;
454 for (int i = 0; i < 4; i++) {
455 for (int j = 0; j < 4; j++) {
456 int idx = i * 4 + j;
457 double temp = -g_uv[idx] + pa[i] * pa[j] / msa;
458 proj_uv.push_back(temp);
459 }
460 }
461
462 // Orbital Tensors
463 if (rank == 0) {
464 std::vector<double> t;
465 t.push_back(1.0);
466 return t;
467
468 } else if (rank == 1) {
469 return proj_uv;
470 } else if (rank == 2) {
471 std::vector<double> proj_uvmn;
472 for (int i = 0; i < 4; i++) {
473 for (int j = 0; j < 4; j++) {
474 for (int k = 0; k < 4; k++) {
475 for (int l = 0; l < 4; l++) {
476
477 int idx1_1 = 4 * i + k;
478 int idx1_2 = 4 * i + l;
479 int idx1_3 = 4 * i + j;
480
481 int idx2_1 = 4 * j + l;
482 int idx2_2 = 4 * j + k;
483 int idx2_3 = 4 * k + l;
484
485 double temp = (1.0 / 2.0) * (proj_uv[idx1_1] * proj_uv[idx2_1] + proj_uv[idx1_2] * proj_uv[idx2_2]) -
486 (1.0 / 3.0) * proj_uv[idx1_3] * proj_uv[idx2_3];
487 proj_uvmn.push_back(temp);
488 }
489 }
490 }
491 }
492 return proj_uvmn;
493
494 }
495 return std::vector<double>();
496 }
497 double EvtD0Topippim2pi0::fundecaymomentum(double mr2, double m1_2, double m2_2)
498 {
499 double mr = sqrt(mr2);
500 double poly = mr2 * mr2 + m1_2 * m1_2 + m2_2 * m2_2 - 2 * m1_2 * mr2 - 2 * m2_2 * mr2 - 2 * m1_2 * m2_2;
501 double ret = sqrt(poly) / (2 * mr);
502 if (poly < 0)
503 //if(sqrt(m1_2) +sqrt(m2_2) > mr)
504 ret = 0.0f;
505 return ret;
506 }
507
508 std::complex<double> EvtD0Topippim2pi0::breitwigner(double mx2, double mr, double wr)
509 {
510 double output_x = 0;
511 double output_y = 0;
512
513 double mr2 = mr * mr;
514 double diff = mr2 - mx2;
515 double denom = diff * diff + wr * wr * mr2;
516 if (wr < 0) {
517 output_x = 0;
518 output_y = 0;
519 } else if (wr < 10) {
520 output_x = diff / denom;
521 output_y = wr * mr / denom;
522 } else { /* phase space */
523 output_x = 1;
524 output_y = 0;
525 }
526 std::complex<double> output(output_x, output_y);
527 return output;
528 }
529
530 /* GS propagator */
531 double EvtD0Topippim2pi0::h(double m, double q)
532 {
533 double h = 2.0 / math_pi * q / m * log((m + 2.0 * q) / (2.0 * mass_Pion));
534 return h;
535 }
536
537 double EvtD0Topippim2pi0::dh(double m0, double q0)
538 {
539 double dh = h(m0, q0) * (1.0 / (8.0 * q0 * q0) - 1.0 / (2.0 * m0 * m0)) + 1.0 / (2.0 * math_pi * m0 * m0);
540 return dh;
541 }
542
543 double EvtD0Topippim2pi0::f(double m0, double sx, double q0, double q)
544 {
545 double m = sqrt(sx);
546 double f = m0 * m0 / (q0 * q0 * q0) * (q * q * (h(m, q) - h(m0, q0)) + (m0 * m0 - sx) * q0 * q0 * dh(m0, q0));
547 return f;
548 }
549
550 double EvtD0Topippim2pi0::d(double m0, double q0)
551 {
552 double d = 3.0 / math_pi * mass_Pion * mass_Pion / (q0 * q0) * log((m0 + 2.0 * q0) / (2.0 * mass_Pion)) + m0 /
553 (2.0 * math_pi * q0) - (mass_Pion * mass_Pion * m0) / (math_pi * q0 * q0 * q0);
554 return d;
555 }
556
557 double EvtD0Topippim2pi0::fundecaymomentum2(double mr2, double m1_2, double m2_2)
558 {
559 double poly = mr2 * mr2 + m1_2 * m1_2 + m2_2 * m2_2 - 2 * m1_2 * mr2 - 2 * m2_2 * mr2 - 2 * m1_2 * m2_2;
560 double ret = poly / (4.0f * mr2);
561 if (poly < 0)
562 ret = 0.0f;
563 return ret;
564 }
565
566 double EvtD0Topippim2pi0::wid(double mass, double sa, double sb, double sc, double r, int l)
567 {
568 double widm = 1.0;
569 double sa0 = mass * mass;
570 double m = sqrt(sa);
571 double q = fundecaymomentum2(sa, sb, sc);
572 double q0 = fundecaymomentum2(sa0, sb, sc);
573 r = r / 0.197321;
574 double z = q * r * r;
575 double z0 = q0 * r * r;
576 double F = 0.0;
577 if (l == 0) F = 1.0;
578 if (l == 1) F = sqrt((1.0 + z0) / (1.0 + z));
579 if (l == 2) F = sqrt((9.0 + 3.0 * z0 + z0 * z0) / (9.0 + 3.0 * z + z * z));
580 if (l == 3) F = sqrt((225.0 + 45.0 * z0 + 6.0 * z0 * z0 + z0 * z0 * z0) / (225.0 + 45.0 * z + 6.0 * z * z + z * z * z));
581 if (l == 4) F = sqrt((11025.0 + 1575.0 * z0 + 135.0 * z0 * z0 + 10.0 * z0 * z0 * z0 + z0 * z0 * z0 * z0) /
582 (11025.0 + 1575.0 * z + 135.0 * z * z + 10.0 * z * z * z + z * z * z * z));
583 double t = sqrt(q / q0);
584 //printf("sa0 = %f, sb = %f, sc = %f, q = %f, q0 = %0.15f, qq0 = %f, t = %f \n",sa0,sb,sc,q,q0,q/q0,t);
585 uint i = 0;
586 for (i = 0; i < static_cast<uint>(2 * l + 1); i++) {
587 widm *= t;
588 }
589 widm *= (mass / m * F * F);
590 return widm;
591 }
592
593 /* for rho0, use GS, rho0->pi+ pi-, only angular momentum 1*/
594 std::complex<double> EvtD0Topippim2pi0::GS(double mx2, double mr, double wr, double m1_2, double m2_2, double r, int l)
595 {
596
597 double mr2 = mr * mr;
598 double q = fundecaymomentum(mx2, m1_2, m2_2);
599 double q0 = fundecaymomentum(mr2, m1_2, m2_2);
600 double numer = 1.0 + d(mr, q0) * wr / mr;
601 double denom_real = mr2 - mx2 + wr * f(mr, mx2, q0, q);
602 double denom_imag = mr * wr * wid(mr, mx2, m1_2, m2_2, r, l); //real-i*imag;
603
604 double denom = denom_real * denom_real + denom_imag * denom_imag;
605 double output_x = denom_real * numer / denom;
606 double output_y = denom_imag * numer / denom;
607
608 std::complex<double> output(output_x, output_y);
609 return output;
610 }
611
612 std::complex<double> EvtD0Topippim2pi0::irho(double mr2, double m1_2, double m2_2)
613 {
614 double poly = mr2 * mr2 + m1_2 * m1_2 + m2_2 * m2_2 - 2 * m1_2 * mr2 - 2 * m2_2 * mr2 - 2 * m1_2 * m2_2;
615 double ret_real = 0;
616 double ret_imag = 0;
617 if (poly >= 0) {
618 ret_real = 2.0f * sqrt(poly) / (2.0f * mr2);
619 ret_imag = 0;
620 } else {
621 ret_real = 0;
622 ret_imag = 2.0f * sqrt(-1.0f * poly) / (2.0f * mr2);
623 }
624 std::complex<double> ret(ret_real, ret_imag);
625 return ret;
626 }
627
628 std::complex<double> EvtD0Topippim2pi0::Flatte(double mx2, double mr, double g1, double g2, double m1a, double m1b, double m2a,
629 double m2b)
630 {
631
632 double mr2 = mr * mr;
633 double diff = mr2 - mx2;
634 double g22 = g2 * g1;
635 std::complex<double> ps1 = irho(mx2, m1a * m1a, m1b * m1b);
636 std::complex<double> ps2 = irho(mx2, m2a * m2a, m2b * m2b);
637 std::complex<double> iws = g1 * ps1 + g22 * ps2; /*mass dependent width */
638
639 double denom_real = diff + iws.imag();
640 double denom_imag = iws.real();
641 double denom = denom_real * denom_real + denom_imag * denom_imag;
642
643 double output_x = denom_real / denom;
644 double output_y = denom_imag / denom;
645
646 std::complex<double> output(output_x, output_y);
647 return output;
648
649 }
650
651 std::complex<double> EvtD0Topippim2pi0::RBW(double mx2, double mr, double wr, double m1_2, double m2_2, double r, int l)
652 {
653 double mr2 = mr * mr;
654 double denom_real = mr2 - mx2;
655 double denom_imag = 0;
656 if (m1_2 > 0 && m2_2 > 0) {
657 denom_imag = mr * wr * wid(mr, mx2, m1_2, m2_2, r, l); //real-i*imag;
658 } else {
659 denom_imag = mr * wr;
660 }
661
662 double denom = denom_real * denom_real + denom_imag * denom_imag;
663 double output_x = denom_real / denom;
664 double output_y = denom_imag / denom;
665
666 std::complex<double> output(output_x, output_y);
667 return output;
668 }
669
670 /* build a1260 propagator */
671 double EvtD0Topippim2pi0::widT1260(int i, double g1, double g2)
672 {
673 double wid1[300] = { 0.000765966, 0.00526534, 0.0168824, 0.0381541, 0.0709798, 0.116716, 0.176348, 0.250811, 0.340516, 0.44606,
674 0.567244, 0.704921, 0.859082, 1.02965, 1.21689, 1.42109, 1.64283, 1.88082, 2.13774, 2.41039,
675 2.70334, 3.0139, 3.34224, 3.69257, 4.05705, 4.44768, 4.85524, 5.28944, 5.73832, 6.21287,
676 6.71567, 7.2421, 7.79776, 8.37191, 8.98624, 9.63238, 10.3101, 11.0061, 11.7641, 12.5458,
677 13.3764, 14.2363, 15.1628, 16.1227, 17.1638, 18.2475, 19.3925, 20.603, 21.8923, 23.2399,
678 24.6837, 26.2041, 27.8182, 29.5295, 31.3134, 33.246, 35.249, 37.3791, 39.6217, 41.9618,
679 44.3717, 46.8669, 49.4843, 52.1752, 54.9009, 57.5853, 60.3722, 63.1059, 65.7868, 68.4831,
680 71.0904, 73.6104, 76.1854, 78.5089, 80.9245, 83.0936, 85.2363, 87.3891, 89.3678, 91.2757,
681 93.049, 94.7878, 96.4351, 98.0146, 99.4775, 100.97, 102.344, 103.571, 104.91, 106.134,
682 107.27, 108.404, 109.517, 110.378, 111.45, 112.438, 113.096, 114.086, 114.935, 115.693,
683 116.439, 117.249, 117.822, 118.369, 119.062, 119.631, 120.217, 120.995, 121.537, 121.963,
684 122.856, 123.129, 123.642, 124.165, 124.729, 124.969, 125.572, 125.878, 126.373, 126.837,
685 127.287, 127.575, 128.226, 128.539, 128.808, 129.144, 129.702, 129.988, 130.595, 130.759,
686 131.101, 131.251, 131.637, 131.951, 132.552, 132.804, 132.991, 133.455, 133.674, 133.95,
687 134.351, 134.674, 134.97, 135.174, 135.429, 135.94, 136.292, 136.595, 136.786, 137.014,
688 137.454, 137.782, 138.096, 138.396, 138.629, 139.002, 139.213, 139.419, 139.819, 140.165,
689 140.352, 140.579, 140.992, 141.404, 141.577, 141.786, 142.224, 142.476, 142.807, 143.216,
690 143.697, 144.059, 144.326, 144.735, 144.951, 145.342, 145.887, 145.88, 146.382, 146.848,
691 147.207, 147.553, 147.962, 148.267, 148.676, 149.175, 149.612, 149.996, 150.341, 150.7,
692 151.071, 151.478, 152.081, 152.566, 152.899, 153.305, 153.923, 154.102, 154.68, 155.234,
693 155.518, 156.153, 156.547, 157.044, 157.414, 158.057, 158.427, 158.869, 159.4, 159.95,
694 160.321, 160.917, 161.523, 161.92, 162.626, 162.708, 163.281, 163.842, 164.594, 164.975,
695 165.366, 166.104, 166.274, 166.837, 167.451, 168.012, 168.635, 169.078, 169.339, 170.036,
696 170.683, 171.155, 171.515, 172.134, 172.733, 173.207, 173.756, 174.176, 174.769, 175.036,
697 175.702, 176.163, 176.59, 177.36, 177.654, 178.347, 178.708, 179.455, 180.003, 180.334,
698 180.716, 181.446, 181.945, 182.401, 182.778, 183.317, 184.106, 184.271, 184.896, 185.414,
699 185.818, 186.491, 186.966, 187.533, 188.148, 188.61, 189.022, 189.321, 189.907, 190.607,
700 191.104, 191.341, 191.888, 192.648, 193.027, 193.363, 194.166, 194.692, 194.995, 195.407,
701 195.898, 196.352, 197.005, 197.378, 198.08, 198.479, 198.95, 199.224, 199.647, 200.37,
702 201.205, 201.251, 201.629, 202.244, 202.762, 203.404, 203.65, 204.172, 204.751, 204.817
703 };
704
705 double wid2[300] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
706 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
707 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
708 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
709 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
710 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
711 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
712 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
713 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
714 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
715 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
716 1.87136e-06, 1.50063e-05, 5.10425e-05, 0.000122121, 0.000240853, 0.000420318, 0.000675161, 0.0010173, 0.00146434, 0.00203321,
717 0.00273489, 0.0035927, 0.00462579, 0.00584255, 0.00727372, 0.00895462, 0.0108831, 0.013085, 0.0156197, 0.0184865,
718 0.0217078, 0.0253423, 0.0294103, 0.0339191, 0.0389837, 0.0446351, 0.0508312, 0.0577268, 0.0653189, 0.0737049,
719 0.0829819, 0.0930611, 0.104328, 0.116663, 0.130105, 0.144922, 0.16122, 0.179091, 0.198759, 0.220133,
720 0.243916, 0.269803, 0.298861, 0.330061, 0.365741, 0.40437, 0.447191, 0.49501, 0.548576, 0.606445,
721 0.674414, 0.748353, 0.831686, 0.929938, 1.03771, 1.16187, 1.30387, 1.47341, 1.65629, 1.88318,
722 2.14353, 2.44169, 2.79831, 3.2009, 3.65522, 4.16317, 4.69597, 5.2585, 5.85965, 6.44984,
723 7.04202, 7.60113, 8.14571, 8.73195, 9.24537, 9.75717, 10.2093, 10.6731, 11.1487, 11.5819,
724 12.0158, 12.4253, 12.8113, 13.2073, 13.5995, 13.9317, 14.312, 14.6595, 14.9511, 15.2668,
725 15.6092, 15.9349, 16.1873, 16.5049, 16.819, 17.0743, 17.3621, 17.6094, 17.8418, 18.0681,
726 18.3141, 18.5914, 18.8187, 19.0562, 19.2282, 19.4918, 19.7326, 19.9112, 20.134, 20.3386,
727 20.511, 20.6865, 20.8958, 21.0518, 21.2967, 21.44, 21.6361, 21.8012, 21.9523, 22.1736,
728 22.2615, 22.4207, 22.6056, 22.7198, 22.9299, 23.0605, 23.2959, 23.3808, 23.4961, 23.6793,
729 23.7843, 23.9697, 24.0689, 24.1919, 24.405, 24.3898, 24.6018, 24.7294, 24.789, 24.9978,
730 25.0626, 25.1728, 25.2809, 25.3579, 25.5444, 25.5995, 25.7644, 25.8397, 25.9229, 26.095,
731 26.1495, 26.2899, 26.3871, 26.54, 26.6603, 26.7008, 26.7836, 26.907, 26.9653, 26.9969,
732 27.1226, 27.226, 27.3543, 27.4686, 27.4887, 27.6163, 27.6986, 27.7506, 27.7884, 27.8662,
733 27.9886, 28.0573, 28.1238, 28.2612, 28.3209, 28.3457, 28.4392, 28.5086, 28.6399, 28.7603,
734 28.788, 28.8502, 28.9038, 28.9667, 28.975, 29.0032, 29.2681, 29.2392, 29.2572, 29.3364
735 };
736
737 return wid1[i] * g1 + wid2[i] * g2;
738 }
739
740 double EvtD0Topippim2pi0::anywid1260(double sc, double g1, double g2)
741 {
742
743 double smin = (0.13957 * 3) * (0.13957 * 3);
744 double dh = 0.01;
745 int od = (sc - 0.18) / dh;
746 double sc_m = 0.18 + od * dh;
747 double widuse = 0;
748 if (sc >= 0.18 && sc <= 3.17) {
749 widuse = ((sc - sc_m) / dh) * (widT1260(od + 1, g1, g2) - widT1260(od, g1, g2)) + widT1260(od, g1, g2);
750 } else if (sc < 0.18 && sc > smin) {
751 widuse = ((sc - smin) / (0.18 - smin)) * widT1260(0, g1, g2);
752 } else if (sc > 3.17) {
753 widuse = widT1260(299, g1, g2);
754 } else {
755 widuse = 0;
756 }
757 return widuse;
758
759 }
760
761 std::complex<double> EvtD0Topippim2pi0::RBWa1260(double mx2, double mr, double g1, double g2)
762 {
763
764 double mx = sqrt(mx2);
765 double mr2 = mr * mr;
766 double wid0 = anywid1260(mx2, g1, g2);
767
768 double denom_real = mr2 - mx2;
769 double denom_imag = mx * wid0; //real-i*imag;
770
771 double denom = denom_real * denom_real + denom_imag * denom_imag;
772 double output_x = denom_real / denom;
773 double output_y = denom_imag / denom;
774
775 std::complex<double> output(output_x, output_y);
776 return output;
777
778 }
779
780 /* build pi1300 propagator */
781 double EvtD0Topippim2pi0::widT1300(int i)
782 {
783 double wid1[300] = { 0.032058, 0.181916, 0.451826, 0.828744, 1.30099, 1.85751, 2.48953, 3.18916, 3.9474, 4.75752,
784 5.61567, 6.5142, 7.4487, 8.41743, 9.41268, 10.4375, 11.4806, 12.5426, 13.6259, 14.7324,
785 15.8316, 16.9667, 18.1034, 19.2511, 20.4, 21.5745, 22.7364, 23.9313, 25.1231, 26.331,
786 27.5222, 28.7648, 29.9995, 31.2184, 32.4692, 33.7247, 34.9821, 36.2849, 37.5859, 38.9141,
787 40.208, 41.5699, 42.8807, 44.2708, 45.6339, 47.0794, 48.5124, 49.9607, 51.4382, 52.9726,
788 54.5431, 56.0448, 57.755, 59.4508, 61.0778, 62.8293, 64.7519, 66.559, 68.5168, 70.4768,
789 72.5278, 74.6728, 76.8742, 79.0858, 81.4848, 83.8468, 86.3982, 88.9172, 91.645, 94.4245,
790 97.2676, 100.073, 103.181, 106.146, 109.37, 112.451, 115.574, 118.859, 122.118, 125.663,
791 128.847, 132.328, 135.778, 139.325, 142.708, 146.613, 150.058, 153.466, 157.281, 160.692,
792 164.484, 168.1, 171.52, 175.195, 178.72, 182.098, 185.828, 189.775, 192.768, 196.145,
793 199.837, 203.808, 206.594, 209.711, 213.22, 216.724, 219.613, 223.289, 226.325, 228.679,
794 232.516, 234.97, 237.984, 241.407, 244.237, 247.025, 250.195, 252.545, 255.669, 258.907,
795 262.137, 264.45, 267.707, 270.006, 273.093, 275.135, 278.181, 281.148, 284.513, 286.352,
796 288.876, 291.133, 293.855, 296.082, 298.943, 301.294, 303.401, 305.955, 308.498, 310.318,
797 312.861, 314.59, 317.381, 319.36, 321.528, 323.834, 326.465, 328.516, 330.043, 331.913,
798 334.528, 336.647, 338.083, 340.873, 342.284, 344.527, 346.21, 347.898, 349.829, 351.885,
799 353.297, 355.228, 357.187, 359.115, 360.499, 362.05, 363.388, 365.451, 367.202, 369.099,
800 371.197, 372.755, 374.302, 375.418, 377.114, 378.595, 380.72, 381.193, 383.21, 385.243,
801 386.41, 387.791, 388.97, 390.617, 391.731, 393.692, 394.627, 396.209, 397.339, 398.582,
802 399.836, 401.752, 402.704, 404.86, 405.287, 406.561, 408.001, 408.361, 410.374, 411.672,
803 412.688, 414.259, 415.406, 416.22, 417.659, 418.652, 419.82, 420.586, 421.835, 423.323,
804 423.714, 425.449, 426.524, 427.196, 428.974, 429.086, 430.614, 431.642, 432.744, 434.121,
805 434.083, 436.251, 436.024, 437.267, 438.584, 439.499, 440.402, 441.609, 441.785, 443.319,
806 444.053, 445.303, 445.359, 446.906, 447.629, 448.931, 449.379, 450.143, 451.303, 452.034,
807 453.026, 453.73, 454.424, 455.852, 456.574, 458.04, 457.593, 459.097, 459.977, 461.041,
808 461.527, 462.486, 463.506, 464.275, 465.849, 465.581, 467.364, 467.648, 468.615, 469.74,
809 470.1, 471.406, 472.099, 472.723, 474.476, 474.81, 475.291, 476.741, 477.751, 478.678,
810 479.833, 479.724, 480.47, 481.717, 482.226, 483.068, 484.311, 485.806, 485.987, 487.431,
811 487.108, 488.798, 488.758, 489.889, 491.406, 491.455, 492.192, 493.215, 493.758, 494.763,
812 496.16, 496.003, 496.981, 497.871, 498.674, 499.319, 500.16, 500.455, 501.069, 502.515
813 };
814
815 return wid1[i];
816 }
817
818 double EvtD0Topippim2pi0::anywid1300(double sc)
819 {
820
821 double smin = (0.13957 * 3) * (0.13957 * 3);
822 double dh = 0.01;
823 int od = (sc - 0.18) / dh;
824 double sc_m = 0.18 + od * dh;
825 double widuse = 0;
826 if (sc >= 0.18 && sc <= 3.17) {
827 widuse = ((sc - sc_m) / dh) * (widT1300(od + 1) - widT1300(od)) + widT1300(od);
828 } else if (sc < 0.18 && sc > smin) {
829 widuse = ((sc - smin) / (0.18 - smin)) * widT1300(0);
830 } else if (sc > 3.17) {
831 widuse = widT1300(299);
832 } else {
833 widuse = 0;
834 }
835 return widuse;
836 }
837
838 std::complex<double> EvtD0Topippim2pi0::RBWpi1300(double mx2, double mr, double wr)
839 {
840
841 double mx = sqrt(mx2);
842 double mr2 = mr * mr;
843 double g1 = wr / anywid1300(mr2);
844 double wid0 = anywid1300(mx2) * g1;
845
846 double denom_real = mr2 - mx2;
847 double denom_imag = mx * wid0; //real-i*imag;
848
849 double denom = denom_real * denom_real + denom_imag * denom_imag;
850 double output_x = denom_real / denom;
851 double output_y = denom_imag / denom;
852
853 std::complex<double> output(output_x, output_y);
854 return output;
855
856 }
857
858 /* build a1640 propagator */
859 double EvtD0Topippim2pi0::widT1640(int i)
860 {
861 double wid1[300] = { 0.000266799, 0.0018334, 0.00594952, 0.0136464, 0.0257838, 0.0430639, 0.0660891, 0.0954588, 0.131591, 0.175007,
862 0.225896, 0.284916, 0.352348, 0.42844, 0.513613, 0.608422, 0.713236, 0.827844, 0.953841, 1.09028,
863 1.23896, 1.3999, 1.57261, 1.75997, 1.95828, 2.17416, 2.40274, 2.65003, 2.91017, 3.18896,
864 3.4876, 3.80491, 4.14473, 4.50064, 4.88498, 5.29342, 5.72815, 6.18169, 6.67661, 7.19577,
865 7.75069, 8.33372, 8.96242, 9.62548, 10.3429, 11.1008, 11.9068, 12.7649, 13.684, 14.6561,
866 15.7, 16.8052, 17.9946, 19.2573, 20.5822, 22.02, 23.5335, 25.1379, 26.8401, 28.6285,
867 30.4874, 32.4225, 34.4622, 36.5651, 38.7177, 40.8728, 43.1028, 45.3167, 47.5153, 49.7404,
868 51.913, 54.0368, 56.2129, 58.2269, 60.3061, 62.2158, 64.109, 66.007, 67.7917, 69.5455,
869 71.1901, 72.8026, 74.3817, 75.8687, 77.2941, 78.7569, 80.1295, 81.3587, 82.6853, 83.9199,
870 85.0982, 86.2791, 87.4225, 88.3923, 89.5057, 90.5393, 91.3378, 92.4117, 93.3536, 94.2291,
871 95.0644, 95.9397, 96.6709, 97.3455, 98.1876, 98.8578, 99.5411, 100.383, 101.071, 101.612,
872 102.565, 102.933, 103.563, 104.141, 104.805, 105.172, 105.82, 106.22, 106.798, 107.369,
873 107.889, 108.279, 108.856, 109.296, 109.633, 110.055, 110.639, 110.974, 111.553, 111.775,
874 112.189, 112.364, 112.693, 113.063, 113.737, 113.872, 114.177, 114.604, 114.814, 115.086,
875 115.475, 115.797, 116.066, 116.276, 116.453, 116.959, 117.212, 117.536, 117.714, 117.784,
876 118.206, 118.435, 118.77, 118.937, 119.078, 119.443, 119.54, 119.649, 119.903, 120.219,
877 120.308, 120.402, 120.661, 121.042, 121.051, 121.058, 121.445, 121.57, 121.593, 121.896,
878 122.146, 122.339, 122.519, 122.694, 122.665, 122.943, 123.242, 123.137, 123.276, 123.516,
879 123.635, 123.746, 124.053, 123.942, 124.071, 124.211, 124.4, 124.546, 124.776, 124.811,
880 124.872, 124.814, 125.187, 125.28, 125.261, 125.33, 125.61, 125.509, 125.683, 125.777,
881 125.698, 126.017, 126.035, 126.269, 126.338, 126.377, 126.655, 126.603, 126.602, 126.721,
882 126.917, 126.91, 127.064, 127.141, 127.193, 127.105, 127.27, 127.377, 127.566, 127.461,
883 127.571, 127.747, 127.719, 127.763, 127.834, 128.103, 128.173, 128.217, 128.021, 128.291,
884 128.348, 128.376, 128.528, 128.489, 128.74, 128.767, 128.858, 128.845, 128.925, 128.965,
885 128.947, 129.009, 129.118, 129.165, 129.219, 129.277, 129.356, 129.26, 129.595, 129.474,
886 129.422, 129.679, 129.526, 129.699, 129.702, 129.64, 129.879, 129.773, 129.798, 130.091,
887 130.023, 130.004, 130.118, 129.995, 130.182, 130.233, 130.262, 130.584, 130.462, 130.547,
888 130.551, 130.513, 130.501, 130.627, 130.632, 130.679, 130.908, 130.924, 130.894, 130.986,
889 131.024, 131.075, 130.995, 130.902, 131.228, 131.115, 131.17, 131.257, 131.16, 131.366,
890 131.439, 131.397, 131.386, 131.472, 131.48, 131.564, 131.505, 131.614, 131.573, 131.455
891 };
892 return wid1[i];
893 }
894
895 double EvtD0Topippim2pi0::anywid1640(double sc)
896 {
897
898 double smin = (0.13957 * 3) * (0.13957 * 3);
899 double dh = 0.01;
900 int od = (sc - 0.18) / dh;
901 double sc_m = 0.18 + od * dh;
902 double widuse = 0;
903 if (sc >= 0.18 && sc <= 3.17) {
904 widuse = ((sc - sc_m) / dh) * (widT1640(od + 1) - widT1640(od)) + widT1640(od);
905 } else if (sc < 0.18 && sc > smin) {
906 widuse = ((sc - smin) / (0.18 - smin)) * widT1640(0);
907 } else if (sc > 3.17) {
908 widuse = widT1640(299);
909 } else {
910 widuse = 0;
911 }
912 return widuse;
913 }
914
915 std::complex<double> EvtD0Topippim2pi0::RBWa1640(double mx2, double mr, double wr)
916 {
917
918 double mx = sqrt(mx2);
919 double mr2 = mr * mr;
920 double g1 = wr / anywid1640(mr2);
921 double wid0 = anywid1640(mx2) * g1;
922
923 double denom_real = mr2 - mx2;
924 double denom_imag = mx * wid0; //real-i*imag;
925
926 double denom = denom_real * denom_real + denom_imag * denom_imag;
927 double output_x = denom_real / denom;
928 double output_y = denom_imag / denom;
929
930 std::complex<double> output(output_x, output_y);
931 return output;
932
933 }
934 /* build h11170 propagator */
935 double EvtD0Topippim2pi0::widT1170(int i)
936 {
937 double wid1[300] = { 0.000166175, 0.00161921, 0.00561815, 0.0131971, 0.0252114, 0.0422987, 0.0651321, 0.0943468, 0.130226, 0.173219,
938 0.22388, 0.28243, 0.349257, 0.424971, 0.50948, 0.603832, 0.707638, 0.821261, 0.946341, 1.08278,
939 1.22851, 1.38899, 1.56053, 1.74557, 1.94153, 2.15585, 2.38028, 2.62579, 2.88376, 3.16032,
940 3.452, 3.77036, 4.10553, 4.45572, 4.8352, 5.23722, 5.66246, 6.11599, 6.60404, 7.12086,
941 7.658, 8.23991, 8.85135, 9.50382, 10.2039, 10.9498, 11.7427, 12.5806, 13.4853, 14.4446,
942 15.4604, 16.5243, 17.7062, 18.9394, 20.219, 21.6307, 23.1328, 24.6832, 26.3629, 28.1076,
943 29.9211, 31.8286, 33.8171, 35.8791, 38.0244, 40.1354, 42.376, 44.5605, 46.8186, 48.9738,
944 51.1693, 53.2751, 55.4668, 57.5333, 59.6529, 61.5273, 63.4519, 65.3347, 67.1338, 68.9518,
945 70.5915, 72.1639, 73.7419, 75.2814, 76.7236, 78.3086, 79.6492, 80.8699, 82.228, 83.3878,
946 84.6403, 85.8205, 86.8573, 87.9557, 88.9784, 89.9703, 90.8957, 92.0648, 92.8484, 93.755,
947 94.603, 95.5404, 96.2543, 96.8307, 97.8023, 98.5217, 99.1155, 100.047, 100.693, 101.165,
948 102.153, 102.449, 103.126, 103.777, 104.336, 104.807, 105.412, 105.714, 106.345, 107.009,
949 107.679, 107.941, 108.403, 108.962, 109.354, 109.726, 110.228, 110.788, 111.304, 111.496,
950 111.885, 112.06, 112.414, 112.681, 113.347, 113.465, 113.829, 114.221, 114.53, 114.718,
951 115.217, 115.391, 115.661, 115.903, 116.177, 116.576, 116.951, 117.321, 117.423, 117.454,
952 117.903, 118.201, 118.416, 118.627, 118.788, 119.237, 119.307, 119.397, 119.601, 119.892,
953 119.853, 120.044, 120.214, 120.734, 120.762, 120.81, 121.148, 121.317, 121.209, 121.544,
954 121.786, 121.945, 122.121, 122.401, 122.272, 122.593, 122.968, 122.78, 122.868, 123.262,
955 123.316, 123.424, 123.822, 123.683, 123.577, 123.93, 124.017, 124.27, 124.515, 124.476,
956 124.516, 124.54, 124.82, 124.951, 124.829, 124.911, 125.281, 125.126, 125.511, 125.315,
957 125.354, 125.749, 125.629, 125.968, 126.13, 126.15, 126.372, 126.217, 126.39, 126.451,
958 126.566, 126.618, 126.639, 126.94, 126.874, 126.738, 126.929, 126.999, 127.267, 127.057,
959 127.188, 127.427, 127.424, 127.414, 127.407, 127.749, 127.938, 127.776, 127.687, 128.017,
960 128.04, 128.019, 128.096, 128.156, 128.281, 128.354, 128.496, 128.675, 128.553, 128.666,
961 128.536, 128.63, 128.866, 128.867, 128.908, 128.975, 129.05, 128.894, 129.29, 129.062,
962 129.156, 129.361, 129.105, 129.346, 129.276, 129.368, 129.612, 129.517, 129.494, 129.864,
963 129.677, 129.68, 129.696, 129.629, 129.826, 129.898, 129.984, 130.271, 130.179, 130.301,
964 130.193, 130.208, 130.212, 130.354, 130.244, 130.318, 130.549, 130.554, 130.477, 130.529,
965 130.765, 130.753, 130.709, 130.598, 130.809, 130.67, 130.792, 130.87, 130.857, 131.148,
966 131.145, 130.977, 131.056, 131.149, 131.09, 131.106, 131.156, 131.28, 131.184, 131.122
967 };
968 return wid1[i];
969 }
970
971 double EvtD0Topippim2pi0::anywid1170(double sc)
972 {
973
974 double smin = (0.13957 * 3) * (0.13957 * 3);
975 double dh = 0.01;
976 int od = (sc - 0.18) / dh;
977 double sc_m = 0.18 + od * dh;
978 double widuse = 0;
979 if (sc >= 0.18 && sc <= 3.17) {
980 widuse = ((sc - sc_m) / dh) * (widT1170(od + 1) - widT1170(od)) + widT1170(od);
981 } else if (sc < 0.18 && sc > smin) {
982 widuse = ((sc - smin) / (0.18 - smin)) * widT1170(0);
983 } else if (sc > 3.17) {
984 widuse = widT1170(299);
985 } else {
986 widuse = 0;
987 }
988 return widuse;
989 }
990
991 std::complex<double> EvtD0Topippim2pi0::RBWh11170(double mx2, double mr, double wr)
992 {
993
994 double mx = sqrt(mx2);
995 double mr2 = mr * mr;
996 double g1 = wr / anywid1170(mr2);
997 double wid0 = anywid1170(mx2) * g1;
998
999 double denom_real = mr2 - mx2;
1000 double denom_imag = mx * wid0; //real-i*imag;
1001
1002 double denom = denom_real * denom_real + denom_imag * denom_imag;
1003 double output_x = denom_real / denom;
1004 double output_y = denom_imag / denom;
1005
1006 std::complex<double> output(output_x, output_y);
1007 return output;
1008 }
1009
1010// PiPi-S wave K-Matrix
1011 double EvtD0Topippim2pi0::rho22(double sc)
1012 {
1013 double rho[689] = { 3.70024e-18, 8.52763e-15, 1.87159e-13, 1.3311e-12, 5.61842e-12, 1.75224e-11, 4.48597e-11, 9.99162e-11, 2.00641e-10, 3.71995e-10,
1014 6.47093e-10, 1.06886e-09, 1.69124e-09, 2.58031e-09, 3.8168e-09, 5.49601e-09, 7.72996e-09, 1.06509e-08, 1.44078e-08, 1.91741e-08,
1015 2.51445e-08, 3.25345e-08, 4.15946e-08, 5.25949e-08, 6.58316e-08, 8.16443e-08, 1.00389e-07, 1.22455e-07, 1.48291e-07, 1.78348e-07,
1016 2.1313e-07, 2.53192e-07, 2.99086e-07, 3.51462e-07, 4.10993e-07, 4.78349e-07, 5.54327e-07, 6.3972e-07, 7.35316e-07, 8.42099e-07,
1017 9.61004e-07, 1.09295e-06, 1.2391e-06, 1.40051e-06, 1.57824e-06, 1.77367e-06, 1.98805e-06, 2.22257e-06, 2.47877e-06, 2.7581e-06,
1018 3.06186e-06, 3.39182e-06, 3.74971e-06, 4.137e-06, 4.5555e-06, 5.00725e-06, 5.4939e-06, 6.01725e-06, 6.57992e-06, 7.18371e-06,
1019 7.83044e-06, 8.52301e-06, 9.26342e-06, 1.00535e-05, 1.08967e-05, 1.17953e-05, 1.27514e-05, 1.37679e-05, 1.48482e-05, 1.59943e-05,
1020 1.72088e-05, 1.84961e-05, 1.98586e-05, 2.12987e-05, 2.28207e-05, 2.44279e-05, 2.61228e-05, 2.79084e-05, 2.97906e-05, 3.17718e-05,
1021 3.38544e-05, 3.60443e-05, 3.8345e-05, 4.07591e-05, 4.32903e-05, 4.59459e-05, 4.87285e-05, 5.16403e-05, 5.46887e-05, 5.7878e-05,
1022 6.12111e-05, 6.46908e-05, 6.83274e-05, 7.21231e-05, 7.60817e-05, 8.0208e-05, 8.45102e-05, 8.89919e-05, 9.36544e-05, 9.85082e-05,
1023 0.000103559, 0.000108812, 0.000114267, 0.000119938, 0.000125827, 0.00013194, 0.000138278, 0.000144857, 0.000151681, 0.000158752,
1024 0.000166074, 0.000173663, 0.000181521, 0.000189652, 0.000198059, 0.000206761, 0.000215761, 0.000225063, 0.00023467, 0.000244599,
1025 0.000254855, 0.00026544, 0.000276357, 0.000287629, 0.00029926, 0.000311253, 0.000323609, 0.000336351, 0.000349483, 0.000363009,
1026 0.000376926, 0.000391264, 0.000406029, 0.000421225, 0.000436848, 0.000452921, 0.000469458, 0.000486461, 0.00050393, 0.00052187,
1027 0.000540322, 0.000559278, 0.000578746, 0.00059872, 0.000619236, 0.0006403, 0.000661911, 0.000684074, 0.000706799, 0.000730127,
1028 0.00075405, 0.000778569, 0.000803686, 0.000829443, 0.000855839, 0.000882879, 0.000910561, 0.000938898, 0.000967939, 0.000997674,
1029 0.00102811, 0.00105923, 0.0010911, 0.0011237, 0.00115706, 0.00119117, 0.00122601, 0.00126168, 0.00129815, 0.00133543,
1030 0.00137351, 0.00141242, 0.00145219, 0.00149283, 0.00153434, 0.0015767, 0.00161995, 0.00166415, 0.00170928, 0.00175534,
1031 0.00180232, 0.00185028, 0.00189924, 0.00194919, 0.00200014, 0.00205207, 0.00210503, 0.0021591, 0.00221421, 0.0022704,
1032 0.00232766, 0.00238602, 0.00244554, 0.00250619, 0.00256799, 0.0026309, 0.002695, 0.00276033, 0.00282689, 0.00289467,
1033 0.00296367, 0.00303389, 0.00310543, 0.0031783, 0.00325244, 0.0033279, 0.0034046, 0.00348275, 0.00356229, 0.00364322,
1034 0.00372555, 0.00380924, 0.00389438, 0.00398104, 0.00406914, 0.00415877, 0.00424985, 0.00434235, 0.00443651, 0.00453224,
1035 0.00462954, 0.00472848, 0.00482894, 0.00493102, 0.00503483, 0.00514029, 0.00524749, 0.0053563, 0.00546675, 0.00557905,
1036 0.0056931, 0.00580901, 0.0059267, 0.00604613, 0.00616735, 0.00629049, 0.00641557, 0.00654254, 0.00667142, 0.00680216,
1037 0.00693472, 0.00706946, 0.00720621, 0.00734497, 0.0074858, 0.00762855, 0.00777338, 0.00792036, 0.00806957, 0.00822087,
1038 0.00837426, 0.00852982, 0.0086875, 0.00884756, 0.00900991, 0.00917447, 0.00934137, 0.00951052, 0.00968194, 0.0098558,
1039 0.010032, 0.0102108, 0.0103919, 0.0105754, 0.0107612, 0.0109496, 0.0111406, 0.0113343, 0.0115305, 0.0117293,
1040 0.0119303, 0.0121343, 0.0123409, 0.0125502, 0.0127623, 0.0129771, 0.0131944, 0.0134145, 0.0136376, 0.0138636,
1041 0.0140924, 0.0143241, 0.0145587, 0.0147959, 0.0150363, 0.0152797, 0.0155262, 0.0157758, 0.0160283, 0.0162838,
1042 0.0165421, 0.016804, 0.0170691, 0.0173374, 0.0176087, 0.0178835, 0.0181612, 0.0184423, 0.0187269, 0.0190149,
1043 0.0193063, 0.0196009, 0.0198991, 0.0202003, 0.0205052, 0.0208137, 0.0211259, 0.0214418, 0.0217611, 0.0220841,
1044 0.0224105, 0.0227406, 0.0230746, 0.0234125, 0.0237542, 0.0240996, 0.0244486, 0.0248012, 0.025158, 0.0255188,
1045 0.0258837, 0.0262527, 0.0266256, 0.0270025, 0.0273833, 0.027768, 0.0281572, 0.0285505, 0.0289483, 0.0293503,
1046 0.0297564, 0.0301665, 0.0305808, 0.0309997, 0.0314231, 0.0318511, 0.0322835, 0.0327205, 0.0331616, 0.0336073,
1047 0.0340576, 0.0345128, 0.0349727, 0.0354373, 0.0359066, 0.0363807, 0.0368589, 0.0373419, 0.0378302, 0.0383234,
1048 0.0388218, 0.0393252, 0.0398336, 0.040347, 0.0408652, 0.041388, 0.0419165, 0.0424502, 0.0429893, 0.0435338,
1049 0.0440833, 0.044638, 0.0451976, 0.0457627, 0.0463338, 0.0469103, 0.047492, 0.0480797, 0.0486729, 0.0492716,
1050 0.0498757, 0.0504852, 0.0511009, 0.0517229, 0.0523503, 0.0529838, 0.0536231, 0.0542678, 0.054918, 0.0555743,
1051 0.0562372, 0.0569065, 0.0575818, 0.0582634, 0.0589511, 0.0596454, 0.0603451, 0.061051, 0.0617635, 0.0624826,
1052 0.0632084, 0.0639409, 0.06468, 0.0654254, 0.0661772, 0.0669346, 0.0676994, 0.0684714, 0.0692503, 0.0700354,
1053 0.0708285, 0.0716277, 0.0724347, 0.0732479, 0.0740671, 0.0748947, 0.0757299, 0.0765715, 0.0774207, 0.0782771,
1054 0.0791407, 0.0800119, 0.0808897, 0.0817743, 0.0826672, 0.0835684, 0.0844769, 0.0853938, 0.0863179, 0.0872493,
1055 0.0881882, 0.0891349, 0.090089, 0.0910523, 0.0920236, 0.093002, 0.0939894, 0.094985, 0.0959887, 0.0970003,
1056 0.0980191, 0.0990454, 0.100081, 0.101126, 0.10218, 0.103242, 0.104312, 0.105392, 0.10648, 0.107576,
1057 0.10868, 0.109793, 0.110916, 0.112048, 0.113188, 0.114339, 0.115498, 0.116666, 0.117843, 0.119028,
1058 0.120223, 0.121427, 0.122641, 0.123865, 0.125098, 0.126342, 0.127595, 0.128857, 0.130128, 0.131409,
1059 0.132701, 0.134002, 0.135314, 0.136635, 0.137966, 0.139308, 0.14066, 0.142022, 0.143394, 0.144774,
1060 0.146166, 0.14757, 0.148985, 0.15041, 0.151845, 0.153291, 0.154749, 0.156215, 0.157694, 0.159182,
1061 0.160682, 0.162194, 0.163718, 0.165251, 0.166797, 0.168354, 0.169921, 0.1715, 0.17309, 0.17469,
1062 0.176304, 0.177929, 0.179566, 0.181216, 0.182878, 0.184553, 0.186238, 0.187934, 0.189642, 0.191362,
1063 0.193096, 0.194842, 0.196602, 0.198374, 0.200158, 0.201954, 0.203764, 0.205586, 0.207421, 0.209266,
1064 0.211124, 0.212997, 0.214882, 0.216783, 0.218697, 0.220624, 0.222565, 0.224518, 0.226486, 0.228466,
1065 0.230458, 0.232463, 0.234484, 0.23652, 0.238569, 0.240633, 0.242711, 0.244803, 0.246909, 0.249031,
1066 0.251165, 0.253313, 0.255475, 0.257649, 0.259841, 0.262051, 0.264274, 0.266514, 0.268768, 0.271036,
1067 0.273319, 0.275618, 0.277932, 0.280259, 0.282602, 0.28496, 0.287338, 0.28973, 0.292138, 0.294563,
1068 0.297003, 0.299458, 0.30193, 0.304417, 0.306919, 0.309437, 0.311972, 0.314526, 0.317095, 0.319684,
1069 0.322289, 0.324911, 0.327551, 0.330205, 0.332876, 0.335567, 0.338271, 0.340993, 0.343736, 0.346496,
1070 0.349272, 0.352065, 0.354878, 0.35771, 0.360561, 0.363426, 0.366311, 0.369212, 0.372128, 0.375067,
1071 0.378027, 0.381006, 0.384001, 0.387014, 0.39005, 0.393106, 0.396181, 0.399271, 0.402384, 0.405513,
1072 0.408661, 0.41183, 0.41502, 0.418233, 0.421462, 0.424709, 0.42798, 0.43127, 0.434583, 0.437914,
1073 0.441267, 0.444637, 0.448022, 0.451434, 0.454868, 0.458328, 0.461805, 0.465302, 0.468821, 0.472364,
1074 0.475928, 0.47951, 0.483119, 0.486748, 0.490397, 0.494066, 0.497758, 0.501477, 0.505217, 0.508977,
1075 0.512762, 0.516567, 0.520394, 0.524247, 0.528125, 0.532027, 0.535947, 0.53989, 0.543852, 0.547844,
1076 0.551863, 0.555904, 0.559966, 0.56406, 0.568177, 0.572312, 0.576471, 0.580662, 0.584875, 0.58911,
1077 0.593373, 0.597653, 0.601965, 0.606301, 0.610663, 0.615051, 0.619465, 0.623907, 0.62837, 0.632863,
1078 0.637383, 0.641924, 0.646494, 0.651091, 0.655708, 0.660356, 0.665027, 0.669732, 0.674464, 0.679227,
1079 0.684016, 0.688827, 0.693664, 0.698532, 0.703428, 0.708353, 0.713307, 0.718283, 0.72329, 0.728322,
1080 0.733387, 0.738479, 0.743605, 0.748763, 0.753949, 0.759163, 0.764407, 0.769674, 0.774973, 0.780311,
1081 0.78567, 0.791057, 0.796476, 0.801922, 0.8074, 0.812919, 0.818466, 0.824044
1082 };
1083
1084 double m2 = 0.13957 * 0.13957;
1085 double smin = (0.13957 * 4) * (0.13957 * 4);
1086 double dh = 0.001;
1087 int od = (sc - 0.312) / dh;
1088 double sc_m = 0.312 + od * dh;
1089 double rhouse = 0;
1090 if (sc >= 0.312 && sc < 1) {
1091 rhouse = ((sc - sc_m) / dh) * (rho[od + 1] - rho[od]) + rho[od];
1092 } else if (sc < 0.312 && sc >= smin) {
1093 rhouse = ((sc - smin) / (0.312 - smin)) * rho[0];
1094 } else if (sc >= 1) {
1095 // rhouse = (1-16*m2/sc)*(1-16*m2/sc)*(1-16*m2/sc)*(1-16*m2/sc)*(1-16*m2/sc);
1096 rhouse = sqrt(1 - 16 * m2 / sc);
1097 } else {
1098 rhouse = 0;
1099 }
1100 return rhouse;
1101 }
1102
1103 std::complex<double> EvtD0Topippim2pi0::rhoMTX(int i, int j, double s)
1104 {
1105
1106 double rhoijx;
1107 double rhoijy = 0.0;
1108 if (i == j && i == 0) {
1109 double m2 = 0.13957 * 0.13957;
1110 if ((1 - (4 * m2) / s) > 0) {
1111 rhoijx = sqrt(1.0f - (4 * m2) / s);
1112 rhoijy = 0;
1113 } else {
1114 rhoijy = sqrt((4 * m2) / s - 1.0f);
1115 rhoijx = 0;
1116 }
1117 }
1118 if (i == j && i == 1) {
1119 double m2 = 0.493677 * 0.493677;
1120 if ((1 - (4 * m2) / s) > 0) {
1121 rhoijx = sqrt(1.0f - (4 * m2) / s);
1122 rhoijy = 0;
1123 } else {
1124 rhoijy = sqrt((4 * m2) / s - 1.0f);
1125 rhoijx = 0;
1126 }
1127 }
1128 if (i == j && i == 2) {
1129 rhoijx = rho22(s);
1130 rhoijy = 0;
1131 }
1132 if (i == j && i == 3) {
1133 double m2 = 0.547862 * 0.547862;
1134 if ((1 - (4 * m2) / s) > 0) {
1135 rhoijx = sqrt(1.0f - (4 * m2) / s);
1136 rhoijy = 0;
1137 } else {
1138 rhoijy = sqrt((4 * m2) / s - 1.0f);
1139 rhoijx = 0;
1140 }
1141 }
1142 if (i == j && i == 4) {
1143 double m_1 = 0.547862;
1144 double m_2 = 0.95778;
1145 double mp2 = (m_1 + m_2) * (m_1 + m_2);
1146 if ((1 - mp2 / s) > 0) {
1147 rhoijx = sqrt(1.0f - mp2 / s);
1148 rhoijy = 0;
1149 } else {
1150 rhoijy = sqrt(mp2 / s - 1.0f);
1151 rhoijx = 0;
1152 }
1153 }
1154
1155 if (i != j) {
1156 rhoijx = 0;
1157 rhoijy = 0;
1158 }
1159 std::complex<double> rhoij(rhoijx, rhoijy);
1160 return rhoij;
1161
1162 }
1163 /* Kij = (sum_a gigj/(ma^2-s) + fij(1-s0)/(s-s0))((s-0.5sampi2)(1-sa0)/(s-sa0)) */
1164 std::complex<double> EvtD0Topippim2pi0::KMTX(int i, int j, double s)
1165 {
1166
1167 double Kijx;
1168 double Kijy;
1169 double mpi = 0.13957;
1170 double m[5] = { 0.65100, 1.20360, 1.55817, 1.21000, 1.82206};
1171
1172 double g1[5] = { 0.22889, -0.55377, 0.00000, -0.39899, -0.34639};
1173 double g2[5] = { 0.94128, 0.55095, 0.00000, 0.39065, 0.31503};
1174 double g3[5] = { 0.36856, 0.23888, 0.55639, 0.18340, 0.18681};
1175 double g4[5] = { 0.33650, 0.40907, 0.85679, 0.19906, -0.00984};
1176 double g5[5] = { 0.18171, -0.17558, -0.79658, -0.00355, 0.22358};
1177
1178 double f1[5] = { 0.23399, 0.15044, -0.20545, 0.32825, 0.35412};
1179
1180 double eps = 1e-11;
1181
1182 double upreal[5] = { 0, 0, 0, 0, 0};
1183 double upimag[5] = { 0, 0, 0, 0, 0};
1184
1185 for (int k = 0; k < 5; k++) {
1186
1187 /* down[k] = (m[k]*m[k]-s)*(m[k]*m[k]-s)+eps*eps;
1188 upreal[k] = (m[k]*m[k]-s)/down[k];
1189 upimag[k] = -1.0f * eps/down[k]; */
1190
1191 double dm2 = m[k] * m[k] - s;
1192 if (fabs(dm2) < eps && dm2 <= 0) dm2 = -eps;
1193 if (fabs(dm2) < eps && dm2 > 0) dm2 = eps;
1194 upreal[k] = 1.0f / dm2;
1195 upimag[k] = 0;
1196 }
1197
1198 double tmp1x = g1[i] * g1[j] * upreal[0] + g2[i] * g2[j] * upreal[1] + g3[i] * g3[j] * upreal[2] + g4[i] * g4[j] * upreal[3] + g5[i]
1199 * g5[j] * upreal[4];
1200 double tmp1y = g1[i] * g1[j] * upimag[0] + g2[i] * g2[j] * upimag[1] + g3[i] * g3[j] * upimag[2] + g4[i] * g4[j] * upimag[3] + g5[i]
1201 * g5[j] * upimag[4];
1202
1203 double tmp2 = 0;
1204 if (i == 0) {
1205 tmp2 = f1[j] * (1 + 3.92637) / (s + 3.92637);
1206 }
1207 if (j == 0) {
1208 tmp2 = f1[i] * (1 + 3.92637) / (s + 3.92637);
1209 }
1210 double tmp3 = (s - 0.5 * mpi * mpi) * (1 + 0.15) / (s + 0.15);
1211
1212 Kijx = (tmp1x + tmp2) * tmp3;
1213 Kijy = (tmp1y) * tmp3;
1214
1215 std::complex<double> Kij(Kijx, Kijy);
1216 return Kij;
1217 }
1218
1219 std::complex<double> EvtD0Topippim2pi0::IMTX(int i, int j)
1220 {
1221
1222 double Iijx;
1223 double Iijy;
1224 if (i == j) {
1225 Iijx = 1;
1226 Iijy = 0;
1227 } else {
1228 Iijx = 0;
1229 Iijy = 0;
1230 }
1231 std::complex<double> Iij(Iijx, Iijy);
1232 return Iij;
1233
1234 }
1235
1236 /* F = I - i*K*rho */
1237 std::complex<double> EvtD0Topippim2pi0::FMTX(double Kijx, double Kijy, double rhojjx, double rhojjy, int i, int j)
1238 {
1239
1240 double Fijx;
1241 double Fijy;
1242
1243 double tmpx = rhojjx * Kijx - rhojjy * Kijy;
1244 double tmpy = rhojjx * Kijy + rhojjy * Kijx;
1245
1246 Fijx = IMTX(i, j).real() + tmpy;
1247 Fijy = -tmpx;
1248
1249 std::complex<double> Fij(Fijx, Fijy);
1250 return Fij;
1251 }
1252
1253 /* inverse for Matrix (I - i*rho*K) using LUP */
1254 double EvtD0Topippim2pi0::FINVMTX(double s, double* FINVx, double* FINVy)
1255 {
1256
1257 int P[5] = { 0, 1, 2, 3, 4};
1258
1259 double Fx[5][5];
1260 double Fy[5][5];
1261
1262 double Ux[5][5];
1263 double Uy[5][5];
1264 double Lx[5][5];
1265 double Ly[5][5];
1266
1267 double UIx[5][5];
1268 double UIy[5][5];
1269 double LIx[5][5];
1270 double LIy[5][5];
1271
1272 for (int k = 0; k < 5; k++) {
1273 double rhokkx = rhoMTX(k, k, s).real();
1274 double rhokky = rhoMTX(k, k, s).imag();
1275 Ux[k][k] = rhokkx;
1276 Uy[k][k] = rhokky;
1277 for (int l = k; l < 5; l++) {
1278 double Kklx = KMTX(k, l, s).real();
1279 double Kkly = KMTX(k, l, s).imag();
1280 Lx[k][l] = Kklx;
1281 Ly[k][l] = Kkly;
1282 Lx[l][k] = Lx[k][l];
1283 Ly[l][k] = Ly[k][l];
1284 }
1285 }
1286
1287 for (int k = 0; k < 5; k++) {
1288 for (int l = 0; l < 5; l++) {
1289 double Fklx = FMTX(Lx[k][l], Ly[k][l], Ux[l][l], Uy[l][l], k, l).real();
1290 double Fkly = FMTX(Lx[k][l], Ly[k][l], Ux[l][l], Uy[l][l], k, l).imag();
1291 Fx[k][l] = Fklx;
1292 Fy[k][l] = Fkly;
1293 }
1294 }
1295
1296 for (int k = 0; k < 5; k++) {
1297 double tmprM = (Fx[k][k] * Fx[k][k] + Fy[k][k] * Fy[k][k]);
1298 int tmpID = 0;
1299 for (int l = k; l < 5; l++) {
1300 double tmprF = (Fx[l][k] * Fx[l][k] + Fy[l][k] * Fy[l][k]);
1301 if (tmprM <= tmprF) {
1302 tmprM = tmprF;
1303 tmpID = l;
1304 }
1305 }
1306
1307 int tmpP = P[k];
1308 P[k] = P[tmpID];
1309 P[tmpID] = tmpP;
1310
1311 for (int l = 0; l < 5; l++) {
1312
1313 double tmpFx = Fx[k][l];
1314 double tmpFy = Fy[k][l];
1315
1316 Fx[k][l] = Fx[tmpID][l];
1317 Fy[k][l] = Fy[tmpID][l];
1318
1319 Fx[tmpID][l] = tmpFx;
1320 Fy[tmpID][l] = tmpFy;
1321
1322 }
1323
1324 for (int l = k + 1; l < 5; l++) {
1325 double rFkk = Fx[k][k] * Fx[k][k] + Fy[k][k] * Fy[k][k];
1326 double Fxlk = Fx[l][k];
1327 double Fylk = Fy[l][k];
1328 double Fxkk = Fx[k][k];
1329 double Fykk = Fy[k][k];
1330 Fx[l][k] = (Fxlk * Fxkk + Fylk * Fykk) / rFkk;
1331 Fy[l][k] = (Fylk * Fxkk - Fxlk * Fykk) / rFkk;
1332 for (int m = k + 1; m < 5; m++) {
1333 Fx[l][m] = Fx[l][m] - (Fx[l][k] * Fx[k][m] - Fy[l][k] * Fy[k][m]);
1334 Fy[l][m] = Fy[l][m] - (Fx[l][k] * Fy[k][m] + Fy[l][k] * Fx[k][m]);
1335 }
1336 }
1337 }
1338
1339 for (int k = 0; k < 5; k++) {
1340 for (int l = 0; l < 5 ; l++) {
1341 if (k == l) {
1342 Lx[k][k] = 1;
1343 Ly[k][k] = 0;
1344 Ux[k][k] = Fx[k][k];
1345 Uy[k][k] = Fy[k][k];
1346 }
1347 if (k > l) {
1348 Lx[k][l] = Fx[k][l];
1349 Ly[k][l] = Fy[k][l];
1350 Ux[k][l] = 0;
1351 Uy[k][l] = 0;
1352 }
1353 if (k < l) {
1354 Ux[k][l] = Fx[k][l];
1355 Uy[k][l] = Fy[k][l];
1356 Lx[k][l] = 0;
1357 Ly[k][l] = 0;
1358 }
1359 }
1360 }
1361
1362 // calculate Inverse for L and U
1363 for (int k = 0; k < 5; k++) {
1364
1365 LIx[k][k] = 1;
1366 LIy[k][k] = 0;
1367
1368 double rUkk = Ux[k][k] * Ux[k][k] + Uy[k][k] * Uy[k][k];
1369 UIx[k][k] = Ux[k][k] / rUkk;
1370 UIy[k][k] = -1.0f * Uy[k][k] / rUkk ;
1371
1372 for (int l = (k + 1); l < 5; l++) {
1373 LIx[k][l] = 0;
1374 LIy[k][l] = 0;
1375 UIx[l][k] = 0;
1376 UIy[l][k] = 0;
1377 }
1378
1379 for (int l = (k - 1); l >= 0; l--) { // U-1
1380 double sx = 0;
1381 double c_sx = 0;
1382 double sy = 0;
1383 double c_sy = 0;
1384 for (int m = l + 1; m <= k; m++) {
1385 sx = sx - c_sx;
1386 double sx_tmp = sx + Ux[l][m] * UIx[m][k] - Uy[l][m] * UIy[m][k];
1387 c_sx = (sx_tmp - sx) - (Ux[l][m] * UIx[m][k] - Uy[l][m] * UIy[m][k]);
1388 sx = sx_tmp;
1389
1390 sy = sy - c_sy;
1391 double sy_tmp = sy + Ux[l][m] * UIy[m][k] + Uy[l][m] * UIx[m][k];
1392 c_sy = (sy_tmp - sy) - (Ux[l][m] * UIy[m][k] + Uy[l][m] * UIx[m][k]);
1393 sy = sy_tmp;
1394 }
1395 UIx[l][k] = -1.0f * (UIx[l][l] * sx - UIy[l][l] * sy);
1396 UIy[l][k] = -1.0f * (UIy[l][l] * sx + UIx[l][l] * sy);
1397 }
1398
1399 for (int l = k + 1; l < 5; l++) { // L-1
1400 double sx = 0;
1401 double c_sx = 0;
1402 double sy = 0;
1403 double c_sy = 0;
1404 for (int m = k; m < l; m++) {
1405 sx = sx - c_sx;
1406 double sx_tmp = sx + Lx[l][m] * LIx[m][k] - Ly[l][m] * LIy[m][k];
1407 c_sx = (sx_tmp - sx) - (Lx[l][m] * LIx[m][k] - Ly[l][m] * LIy[m][k]);
1408 sx = sx_tmp;
1409
1410 sy = sy - c_sy;
1411 double sy_tmp = sy + Lx[l][m] * LIy[m][k] + Ly[l][m] * LIx[m][k];
1412 c_sy = (sy_tmp - sy) - (Lx[l][m] * LIy[m][k] + Ly[l][m] * LIx[m][k]);
1413 sy = sy_tmp;
1414 }
1415 LIx[l][k] = -1.0f * sx;
1416 LIy[l][k] = -1.0f * sy;
1417 }
1418 }
1419
1420 for (int m = 0; m < 5; m++) {
1421 double resX = 0;
1422 double c_resX = 0;
1423 double resY = 0;
1424 double c_resY = 0;
1425 for (int k = 0; k < 5; k++) {
1426 for (int l = 0; l < 5; l++) {
1427 double Plm = 0;
1428 if (P[l] == m) Plm = 1;
1429
1430 resX = resX - c_resX;
1431 double resX_tmp = resX + (UIx[0][k] * LIx[k][l] - UIy[0][k] * LIy[k][l]) * Plm;
1432 c_resX = (resX_tmp - resX) - ((UIx[0][k] * LIx[k][l] - UIy[0][k] * LIy[k][l]) * Plm);
1433 resX = resX_tmp;
1434
1435 resY = resY - c_resY;
1436 double resY_tmp = resY + (UIx[0][k] * LIy[k][l] + UIy[0][k] * LIx[k][l]) * Plm;
1437 c_resY = (resY_tmp - resY) - ((UIx[0][k] * LIy[k][l] + UIy[0][k] * LIx[k][l]) * Plm);
1438 resY = resY_tmp;
1439 }
1440 }
1441 FINVx[m] = resX;
1442 FINVy[m] = resY;
1443 }
1444
1445 return 1.0f;
1446 }
1447
1448 std::complex<double> EvtD0Topippim2pi0::PVTR(int ID, double s)
1449 {
1450
1451 double VPix;
1452 double VPiy;
1453 double m[5] = { 0.65100, 1.20360, 1.55817, 1.21000, 1.82206};
1454
1455 double eps = 1e-11;
1456
1457 /* double down = (m[ID]*m[ID]-s)*(m[ID]*m[ID]-s)+eps*eps;
1458 double upreal = (m[ID]*m[ID]-s)/down;
1459 double upimag = -1.0f * eps/down; */
1460
1461 double dm2 = m[ID] * m[ID] - s;
1462
1463 if (fabs(dm2) < eps && dm2 <= 0) dm2 = -eps;
1464 if (fabs(dm2) < eps && dm2 > 0) dm2 = eps;
1465
1466 VPix = 1.0f / dm2;
1467 VPiy = 0;
1468
1469 std::complex<double> VPi(VPix, VPiy);
1470 return VPi;
1471 }
1472
1473 std::complex<double> EvtD0Topippim2pi0::Fvector(double sa, double s0, int l)
1474 {
1475
1476 double outputx = 0;
1477 double outputy = 0;
1478
1479 double FINVx[5] = {0, 0, 0, 0, 0};
1480 double FINVy[5] = {0, 0, 0, 0, 0};
1481
1482 if (l < 5) {
1483 double g[5][5] = {{ 0.22889, -0.55377, 0.00000, -0.39899, -0.34639},
1484 { 0.94128, 0.55095, 0.00000, 0.39065, 0.31503},
1485 { 0.36856, 0.23888, 0.55639, 0.18340, 0.18681},
1486 { 0.33650, 0.40907, 0.85679, 0.19906, -0.00984},
1487 { 0.18171, -0.17558, -0.79658, -0.00355, 0.22358}
1488 };
1489 double resx = 0;
1490 double c_resx = 0;
1491 double resy = 0;
1492 double c_resy = 0;
1493 double Plx = PVTR(l, sa).real();
1494 double Ply = PVTR(l, sa).imag();
1495 for (int j = 0; j < 5; j++) {
1496 resx = resx - c_resx;
1497 double resx_tmp = resx + (FINVx[j] * g[l][j] * Plx - FINVy[j] * g[l][j] * Ply);
1498 c_resx = (resx_tmp - resx) - (FINVx[j] * g[l][j] * Plx - FINVy[j] * g[l][j] * Ply);
1499 resx = resx_tmp;
1500
1501 resy = resy - c_resy;
1502 double resy_tmp = resy + (FINVx[j] * g[l][j] * Ply + FINVy[j] * g[l][j] * Plx);
1503 c_resy = (resy_tmp - resy) - (FINVx[j] * g[l][j] * Ply + FINVy[j] * g[l][j] * Plx);
1504 resy = resy_tmp;
1505 }
1506 outputx = resx;
1507 outputy = resy;
1508 } else {
1509 int idx = l - 5;
1510 double eps = 1e-11;
1511 double ds = sa - s0;
1512 if (fabs(ds) < eps && ds <= 0) ds = -eps;
1513 if (fabs(ds) < eps && ds > 0) ds = eps;
1514 double tmp = (1 - s0) / ds;
1515 outputx = FINVx[idx] * tmp;
1516 outputy = FINVy[idx] * tmp;
1517 }
1518
1519 std::complex<double> output(outputx, outputy);
1520 return output;
1521 }
1522
1523 std::complex<double> EvtD0Topippim2pi0::CalD0Amp()
1524 {
1525 return Amp(m_Pip, m_Pim, m_Pi01, m_Pi02);
1526 }
1527 std::complex<double> EvtD0Topippim2pi0::CalDbAmp()
1528 {
1529
1530 std::vector<double> cpPip;
1531 std::vector<double> cpPim;
1532 std::vector<double> cpPi01;
1533 std::vector<double> cpPi02;
1534
1535 cpPip.push_back(-m_Pim[0]); cpPim.push_back(-m_Pip[0]); cpPi01.push_back(-m_Pi01[0]); cpPi02.push_back(-m_Pi02[0]);
1536 cpPip.push_back(-m_Pim[1]); cpPim.push_back(-m_Pip[1]); cpPi01.push_back(-m_Pi01[1]); cpPi02.push_back(-m_Pi02[1]);
1537 cpPip.push_back(-m_Pim[2]); cpPim.push_back(-m_Pip[2]); cpPi01.push_back(-m_Pi01[2]); cpPi02.push_back(-m_Pi02[2]);
1538 cpPip.push_back(m_Pim[3]); cpPim.push_back(m_Pip[3]); cpPi01.push_back(m_Pi01[3]); cpPi02.push_back(m_Pi02[3]);
1539
1540 return Amp(cpPip, cpPim, cpPi01, cpPi02);
1541 }
1542
1543 std::complex<double> EvtD0Topippim2pi0::Amp(std::vector<double> Pip, std::vector<double> Pim, std::vector<double> Pi01,
1544 std::vector<double> Pi02)
1545 {
1546
1547 auto PipPim = sum_tensor(Pip, Pim);
1548 auto PipPi01 = sum_tensor(Pip, Pi01);
1549 auto PipPi02 = sum_tensor(Pip, Pi02);
1550 auto PimPi01 = sum_tensor(Pim, Pi01);
1551 auto PimPi02 = sum_tensor(Pim, Pi02);
1552 auto Pi01Pi02 = sum_tensor(Pi01, Pi02);
1553
1554 auto PipPimPi01 = sum_tensor(PipPim, Pi01);
1555 auto PipPimPi02 = sum_tensor(PipPim, Pi02);
1556 auto PipPi01Pi02 = sum_tensor(PipPi01, Pi02);
1557 auto PimPi01Pi02 = sum_tensor(PimPi01, Pi02);
1558
1559 auto D0 = sum_tensor(PipPimPi01, Pi02);
1560
1561 double M2_PipPim = contract_11_0(PipPim, PipPim);
1562 double M2_PipPi01 = contract_11_0(PipPi01, PipPi01);
1563 double M2_PipPi02 = contract_11_0(PipPi02, PipPi02);
1564 double M2_PimPi01 = contract_11_0(PimPi01, PimPi01);
1565 double M2_PimPi02 = contract_11_0(PimPi02, PimPi02);
1566 double M2_Pi01Pi02 = contract_11_0(Pi01Pi02, Pi01Pi02);
1567
1568 double M2_PipPimPi01 = contract_11_0(PipPimPi01, PipPimPi01);
1569 double M2_PipPimPi02 = contract_11_0(PipPimPi02, PipPimPi02);
1570 double M2_PipPi01Pi02 = contract_11_0(PipPi01Pi02, PipPi01Pi02);
1571 double M2_PimPi01Pi02 = contract_11_0(PimPi01Pi02, PimPi01Pi02);
1572
1573 std::complex<double> GS_rho770_pm = GS(M2_PipPim, m0_rho7700, w0_rho7700, m2_Pi, m2_Pi, rRes, 1);
1574 std::complex<double> GS_rho770_p1 = GS(M2_PipPi01, m0_rho770p, w0_rho770p, m2_Pi, m2_Pi0, rRes, 1);
1575 std::complex<double> GS_rho770_p2 = GS(M2_PipPi02, m0_rho770p, w0_rho770p, m2_Pi, m2_Pi0, rRes, 1);
1576 std::complex<double> GS_rho770_m1 = GS(M2_PimPi01, m0_rho770p, w0_rho770p, m2_Pi, m2_Pi0, rRes, 1);
1577 std::complex<double> GS_rho770_m2 = GS(M2_PimPi02, m0_rho770p, w0_rho770p, m2_Pi, m2_Pi0, rRes, 1);
1578 std::complex<double> GS_rho1450_m1 = GS(M2_PimPi01, m0_rho1450, w0_rho1450, m2_Pi, m2_Pi0, rRes, 1);
1579 std::complex<double> GS_rho1450_m2 = GS(M2_PimPi02, m0_rho1450, w0_rho1450, m2_Pi, m2_Pi0, rRes, 1);
1580
1581 std::complex<double> FT_f0980_00 = Flatte(M2_Pi01Pi02, m0_f0980, g1_f0980, g2_f0980, m_Pi, m_Pi, m_Ka, m_Ka);
1582 std::complex<double> RBW_f21270_pm = RBW(M2_PipPim, m0_f21270, w0_f21270, m2_Pi, m2_Pi, rRes, 2);
1583 std::complex<double> RBW_f21270_00 = RBW(M2_Pi01Pi02, m0_f21270, w0_f21270, m2_Pi0, m2_Pi0, rRes, 2);
1584
1585 std::complex<double> PiPiS_pm_0 = Fvector(M2_PipPim, s0_prod, 0);
1586 std::complex<double> PiPiS_00_0 = Fvector(M2_Pi01Pi02, s0_prod, 0);
1587
1588 std::complex<double> PiPiS_pm_1 = Fvector(M2_PipPim, s0_prod, 1);
1589 std::complex<double> PiPiS_00_1 = Fvector(M2_Pi01Pi02, s0_prod, 1);
1590
1591 std::complex<double> PiPiS_pm_5 = Fvector(M2_PipPim, s0_prod, 5);
1592 std::complex<double> PiPiS_00_5 = Fvector(M2_Pi01Pi02, s0_prod, 5);
1593
1594 std::complex<double> PiPiS_pm_6 = Fvector(M2_PipPim, s0_prod, 6);
1595 std::complex<double> PiPiS_00_6 = Fvector(M2_Pi01Pi02, s0_prod, 6);
1596
1597 std::complex<double> RBW_a11260_p = RBWa1260(M2_PipPi01Pi02, m0_a11260, g1_a11260, g2_a11260);
1598 std::complex<double> RBW_a11260_m = RBWa1260(M2_PimPi01Pi02, m0_a11260, g1_a11260, g2_a11260);
1599 std::complex<double> RBW_a11260_01 = RBWa1260(M2_PipPimPi01, m0_a11260, g1_a11260, g2_a11260);
1600 std::complex<double> RBW_a11260_02 = RBWa1260(M2_PipPimPi02, m0_a11260, g1_a11260, g2_a11260);
1601
1602 std::complex<double> RBW_a11420_p = RBW(M2_PipPi01Pi02, m0_a11420, w0_a11420, -1, -1, -1, -1);
1603
1604 std::complex<double> RBW_a11640_p = RBWa1640(M2_PipPi01Pi02, m0_a11640, w0_a11640);
1605 std::complex<double> RBW_a11640_m = RBWa1640(M2_PimPi01Pi02, m0_a11640, w0_a11640);
1606
1607 std::complex<double> RBW_omega_01 = RBW(M2_PipPimPi01, m0_omega, w0_omega, -1, -1, -1, -1);
1608 std::complex<double> RBW_omega_02 = RBW(M2_PipPimPi02, m0_omega, w0_omega, -1, -1, -1, -1);
1609
1610 std::complex<double> RBW_phi_01 = RBW(M2_PipPimPi01, m0_phi, w0_phi, -1, -1, -1, -1);
1611 std::complex<double> RBW_phi_02 = RBW(M2_PipPimPi02, m0_phi, w0_phi, -1, -1, -1, -1);
1612
1613 std::complex<double> RBW_a21320_p = RBW(M2_PipPi01Pi02, m0_a21320, w0_a21320, -1, -1, -1, -1);
1614 std::complex<double> RBW_a21320_m = RBW(M2_PimPi01Pi02, m0_a21320, w0_a21320, -1, -1, -1, -1);
1615
1616 std::complex<double> RBW_pi1300_p = RBWpi1300(M2_PipPi01Pi02, m0_pi1300, w0_pi1300);
1617 std::complex<double> RBW_pi1300_m = RBWpi1300(M2_PimPi01Pi02, m0_pi1300, w0_pi1300);
1618 std::complex<double> RBW_pi1300_01 = RBWpi1300(M2_PipPimPi01, m0_pi1300, w0_pi1300);
1619 std::complex<double> RBW_pi1300_02 = RBWpi1300(M2_PipPimPi02, m0_pi1300, w0_pi1300);
1620
1621 std::complex<double> RBW_h11170_01 = RBWh11170(M2_PipPimPi01, m0_h11170, w0_h11170);
1622 std::complex<double> RBW_h11170_02 = RBWh11170(M2_PipPimPi02, m0_h11170, w0_h11170);
1623
1624 std::complex<double> RBW_pi21670_01 = RBW(M2_PipPimPi01, m0_pi21670, w0_pi21670, -1, -1, -1, -1);
1625 std::complex<double> RBW_pi21670_02 = RBW(M2_PipPimPi02, m0_pi21670, w0_pi21670, -1, -1, -1, -1);
1626
1627 // D->XX Projection
1628 auto Proj1_3p = ProjectionTensors(PipPi01Pi02, 1);
1629 auto Proj1_3m = ProjectionTensors(PimPi01Pi02, 1);
1630 auto Proj1_3z1 = ProjectionTensors(PipPimPi01, 1);
1631 auto Proj1_3z2 = ProjectionTensors(PipPimPi02, 1);
1632
1633 auto Proj2_3p = ProjectionTensors(PipPi01Pi02, 2);
1634 auto Proj2_3m = ProjectionTensors(PimPi01Pi02, 2);
1635 auto Proj2_3z1 = ProjectionTensors(PipPimPi01, 2);
1636 auto Proj2_3z2 = ProjectionTensors(PipPimPi02, 2);
1637
1638 // X->PP Orbital
1639 auto T1_PipPim = OrbitalTensors(PipPim, Pip, Pim, rRes, 1);
1640 auto T1_PipPi01 = OrbitalTensors(PipPi01, Pip, Pi01, rRes, 1);
1641 auto T1_PipPi02 = OrbitalTensors(PipPi02, Pip, Pi02, rRes, 1);
1642 auto T1_PimPi01 = OrbitalTensors(PimPi01, Pim, Pi01, rRes, 1);
1643 auto T1_PimPi02 = OrbitalTensors(PimPi02, Pim, Pi02, rRes, 1);
1644 auto T1_Pi01Pi02 = OrbitalTensors(Pi01Pi02, Pi01, Pi02, rRes, 1);
1645
1646 std::vector<double> T2_PipPim;
1647 std::vector<double> T2_Pi01Pi02;
1648
1649 T2_PipPim = OrbitalTensors(PipPim, Pip, Pim, rRes, 2);
1650 T2_Pi01Pi02 = OrbitalTensors(Pi01Pi02, Pi01, Pi02, rRes, 2);
1651
1652 // X->YP Orbital
1653 auto T1_PipPimPi01 = OrbitalTensors(PipPimPi01, PipPim, Pi01, rRes, 1);
1654 auto T1_PipPimPi02 = OrbitalTensors(PipPimPi02, PipPim, Pi02, rRes, 1);
1655 auto T1_PipPi01Pi02 = OrbitalTensors(PipPi01Pi02, PipPi01, Pi02, rRes, 1);
1656 auto T1_PipPi02Pi01 = OrbitalTensors(PipPi01Pi02, PipPi02, Pi01, rRes, 1);
1657 auto T1_PimPi01Pi02 = OrbitalTensors(PimPi01Pi02, PimPi01, Pi02, rRes, 1);
1658 auto T1_PimPi02Pi01 = OrbitalTensors(PimPi01Pi02, PimPi02, Pi01, rRes, 1);
1659 auto T1_PipPi01Pim = OrbitalTensors(PipPimPi01, PipPi01, Pim, rRes, 1);
1660 auto T1_PipPi02Pim = OrbitalTensors(PipPimPi02, PipPi02, Pim, rRes, 1);
1661 auto T1_PimPi01Pip = OrbitalTensors(PipPimPi01, PimPi01, Pip, rRes, 1);
1662 auto T1_PimPi02Pip = OrbitalTensors(PipPimPi02, PimPi02, Pip, rRes, 1);
1663 auto T1_Pi01Pi02Pip = OrbitalTensors(PipPi01Pi02, Pi01Pi02, Pip, rRes, 1);
1664 auto T1_Pi01Pi02Pim = OrbitalTensors(PimPi01Pi02, Pi01Pi02, Pim, rRes, 1);
1665
1666 auto T2_PipPimPi01 = OrbitalTensors(PipPimPi01, PipPim, Pi01, rRes, 2);
1667 auto T2_PipPimPi02 = OrbitalTensors(PipPimPi02, PipPim, Pi02, rRes, 2);
1668 auto T2_PipPi01Pi02 = OrbitalTensors(PipPi01Pi02, PipPi01, Pi02, rRes, 2);
1669 auto T2_PipPi02Pi01 = OrbitalTensors(PipPi01Pi02, PipPi02, Pi01, rRes, 2);
1670 auto T2_PimPi01Pi02 = OrbitalTensors(PimPi01Pi02, PimPi01, Pi02, rRes, 2);
1671 auto T2_PimPi02Pi01 = OrbitalTensors(PimPi01Pi02, PimPi02, Pi01, rRes, 2);
1672 auto T2_PipPi01Pim = OrbitalTensors(PipPimPi01, PipPi01, Pim, rRes, 2);
1673 auto T2_PipPi02Pim = OrbitalTensors(PipPimPi02, PipPi02, Pim, rRes, 2);
1674 auto T2_PimPi01Pip = OrbitalTensors(PipPimPi01, PimPi01, Pip, rRes, 2);
1675 auto T2_PimPi02Pip = OrbitalTensors(PipPimPi02, PimPi02, Pip, rRes, 2);
1676 auto T2_Pi01Pi02Pip = OrbitalTensors(PipPi01Pi02, Pi01Pi02, Pip, rRes, 2);
1677 auto T2_Pi01Pi02Pim = OrbitalTensors(PimPi01Pi02, Pi01Pi02, Pim, rRes, 2);
1678
1679 // D->XX Orbital
1680 auto T1_2pm12 = OrbitalTensors(D0, PipPim, Pi01Pi02, rD, 1);
1681 auto T1_2p1m2 = OrbitalTensors(D0, PipPi01, PimPi02, rD, 1);
1682 auto T1_2p2m1 = OrbitalTensors(D0, PipPi02, PimPi01, rD, 1);
1683
1684 auto T2_2pm12 = OrbitalTensors(D0, PipPim, Pi01Pi02, rD, 2);
1685 auto T2_2p1m2 = OrbitalTensors(D0, PipPi01, PimPi02, rD, 2);
1686 auto T2_2p2m1 = OrbitalTensors(D0, PipPi02, PimPi01, rD, 2);
1687
1688 // D->XP Orbital
1689 auto T1_3pm = OrbitalTensors(D0, PipPi01Pi02, Pim, rD, 1);
1690 auto T1_3mp = OrbitalTensors(D0, PimPi01Pi02, Pip, rD, 1);
1691 auto T1_3z12 = OrbitalTensors(D0, PipPimPi01, Pi02, rD, 1);
1692 auto T1_3z21 = OrbitalTensors(D0, PipPimPi02, Pi01, rD, 1);
1693
1694 auto T2_3pm = OrbitalTensors(D0, PipPi01Pi02, Pim, rD, 2);
1695 auto T2_3mp = OrbitalTensors(D0, PimPi01Pi02, Pip, rD, 2);
1696 auto T2_3z12 = OrbitalTensors(D0, PipPimPi01, Pi02, rD, 2);
1697 auto T2_3z21 = OrbitalTensors(D0, PipPimPi02, Pi01, rD, 2);
1698
1699 std::complex<double> amplitude(0, 0);
1700
1701 // D0 -> a1(1260)+ {rho(770)+ pi0[S]} pi-
1702 double SF_Ap_S_Vp1P = contract_11_0(contract_21_1(Proj1_3p, T1_PipPi01), T1_3pm);
1703 double SF_Ap_S_Vp2P = contract_11_0(contract_21_1(Proj1_3p, T1_PipPi02), T1_3pm);
1704
1705 amplitude += fitpara[0] * (SF_Ap_S_Vp1P * RBW_a11260_p * GS_rho770_p1 + SF_Ap_S_Vp2P * RBW_a11260_p * GS_rho770_p2);
1706
1707 // D0 -> a1(1260)+ {rho(770)+ pi0[D]} pi-
1708 double SF_Ap_D_Vp1P = contract_11_0(contract_21_1(T2_PipPi01Pi02, T1_PipPi01), T1_3pm);
1709 double SF_Ap_D_Vp2P = contract_11_0(contract_21_1(T2_PipPi02Pi01, T1_PipPi02), T1_3pm);
1710
1711 amplitude += fitpara[1] * (SF_Ap_D_Vp1P * RBW_a11260_p * GS_rho770_p1 + SF_Ap_D_Vp2P * RBW_a11260_p * GS_rho770_p2);
1712
1713 // D0 -> a1(1260)+ {f2(1270) pi+ [P]} pi-
1714 double SF_Ap_P_TP = contract_11_0(contract_21_1(contract_42_2(Proj2_3p, T2_Pi01Pi02), T1_Pi01Pi02Pip), T1_3pm);
1715
1716 amplitude += fitpara[2] * (SF_Ap_P_TP * RBW_a11260_p * RBW_f21270_00);
1717
1718 // D0 -> a1(1260)+ {f0 pi+ [P]} pi-
1719 double SF_Ap_P_SP = contract_11_0(T1_3pm, T1_Pi01Pi02Pip);
1720
1721 amplitude += fitpara[3] * (SF_Ap_P_SP * RBW_a11260_p * PiPiS_00_0);
1722 amplitude += fitpara[4] * (SF_Ap_P_SP * RBW_a11260_p * PiPiS_00_5);
1723 amplitude += fitpara[5] * (SF_Ap_P_SP * RBW_a11260_p * PiPiS_00_6);
1724
1725 // D0 -> a1(1260)- { rho(770)- pi0 [S]} pi+
1726 double SF_Am_S_Vm1P = contract_11_0(contract_21_1(Proj1_3m, T1_PimPi01), T1_3mp);
1727 double SF_Am_S_Vm2P = contract_11_0(contract_21_1(Proj1_3m, T1_PimPi02), T1_3mp);
1728
1729 amplitude += fitpara[6] * fitpara[0] * (SF_Am_S_Vm1P * RBW_a11260_m * GS_rho770_m1 + SF_Am_S_Vm2P * RBW_a11260_m * GS_rho770_m2);
1730
1731 // D0 -> a1(1260)- {rho(770)- pi0[D]} pi+
1732 double SF_Am_D_Vm1P = contract_11_0(contract_21_1(T2_PimPi01Pi02, T1_PimPi01), T1_3mp);
1733 double SF_Am_D_Vm2P = contract_11_0(contract_21_1(T2_PimPi02Pi01, T1_PimPi02), T1_3mp);
1734
1735 amplitude += fitpara[6] * fitpara[1] * (SF_Am_D_Vm1P * RBW_a11260_m * GS_rho770_m1 + SF_Am_D_Vm2P * RBW_a11260_m * GS_rho770_m2);
1736
1737 // D0 -> a1(1260)- {f2(1270) pi- [P]} pi+
1738 double SF_Am_P_TP = contract_11_0(contract_21_1(contract_42_2(Proj2_3m, T2_Pi01Pi02), T1_Pi01Pi02Pim), T1_3mp);
1739
1740 amplitude += fitpara[6] * fitpara[2] * (SF_Am_P_TP * RBW_a11260_m * RBW_f21270_00);
1741
1742 // D0 -> a1(1260)- {f0 pi- [P]} pi+
1743 double SF_Am_P_SP = contract_11_0(T1_3mp, T1_Pi01Pi02Pim);
1744
1745 amplitude += fitpara[6] * fitpara[3] * (SF_Am_P_SP * RBW_a11260_m * PiPiS_00_0);
1746 amplitude += fitpara[6] * fitpara[4] * (SF_Am_P_SP * RBW_a11260_m * PiPiS_00_5);
1747 amplitude += fitpara[6] * fitpara[5] * (SF_Am_P_SP * RBW_a11260_m * PiPiS_00_6);
1748
1749 // D -> a1(1260)0 pi0
1750 double SF_A01_S_Vp1P = contract_11_0(contract_21_1(Proj1_3z1, T1_PipPi01), T1_3z12);
1751 double SF_A02_S_Vp2P = contract_11_0(contract_21_1(Proj1_3z2, T1_PipPi02), T1_3z21);
1752 double SF_A01_S_Vm1P = contract_11_0(contract_21_1(Proj1_3z1, T1_PimPi01), T1_3z12);
1753 double SF_A02_S_Vm2P = contract_11_0(contract_21_1(Proj1_3z2, T1_PimPi02), T1_3z21);
1754 double SF_A01_S_VzP = contract_11_0(contract_21_1(Proj1_3z1, T1_PipPim), T1_3z12);
1755 double SF_A02_S_VzP = contract_11_0(contract_21_1(Proj1_3z2, T1_PipPim), T1_3z21);
1756
1757 amplitude += fitpara[7] * fitpara[0] * (SF_A01_S_Vp1P * RBW_a11260_01 * GS_rho770_p1 + SF_A02_S_Vp2P * RBW_a11260_02 * GS_rho770_p2
1758 + SF_A01_S_Vm1P * RBW_a11260_01 * GS_rho770_m1 + SF_A02_S_Vm2P * RBW_a11260_02 * GS_rho770_m2);
1759
1760 double SF_A01_D_Vp1P = contract_11_0(contract_21_1(T2_PipPi01Pim, T1_PipPi01), T1_3z12);
1761 double SF_A02_D_Vp2P = contract_11_0(contract_21_1(T2_PipPi02Pim, T1_PipPi02), T1_3z21);
1762 double SF_A01_D_Vm1P = contract_11_0(contract_21_1(T2_PimPi01Pip, T1_PimPi01), T1_3z12);
1763 double SF_A02_D_Vm2P = contract_11_0(contract_21_1(T2_PimPi02Pip, T1_PimPi02), T1_3z21);
1764
1765 amplitude += fitpara[7] * fitpara[1] * (SF_A01_D_Vp1P * RBW_a11260_01 * GS_rho770_p1 + SF_A02_D_Vp2P * RBW_a11260_02 * GS_rho770_p2
1766 + SF_A01_D_Vm1P * RBW_a11260_01 * GS_rho770_m1 + SF_A02_D_Vm2P * RBW_a11260_02 * GS_rho770_m2);
1767
1768 double SF_A01_P_TP = contract_11_0(contract_21_1(contract_42_2(Proj2_3z1, T2_PipPim), T1_PipPimPi01), T1_3z12);
1769 double SF_A02_P_TP = contract_11_0(contract_21_1(contract_42_2(Proj2_3z2, T2_PipPim), T1_PipPimPi02), T1_3z21);
1770
1771 amplitude += fitpara[7] * fitpara[2] * (-1.0) * (SF_A01_P_TP * RBW_a11260_01 * RBW_f21270_pm + SF_A02_P_TP * RBW_a11260_02 *
1772 RBW_f21270_pm);
1773
1774 double SF_A01_P_SP = contract_11_0(T1_3z12, T1_PipPimPi01);
1775 double SF_A02_P_SP = contract_11_0(T1_3z21, T1_PipPimPi02);
1776
1777 amplitude += fitpara[7] * fitpara[3] * (-1.0) * (SF_A01_P_SP * RBW_a11260_01 * PiPiS_pm_0 + SF_A02_P_SP * RBW_a11260_02 *
1778 PiPiS_pm_0);
1779 amplitude += fitpara[7] * fitpara[4] * (-1.0) * (SF_A01_P_SP * RBW_a11260_01 * PiPiS_pm_5 + SF_A02_P_SP * RBW_a11260_02 *
1780 PiPiS_pm_5);
1781 amplitude += fitpara[7] * fitpara[5] * (-1.0) * (SF_A01_P_SP * RBW_a11260_01 * PiPiS_pm_6 + SF_A02_P_SP * RBW_a11260_02 *
1782 PiPiS_pm_6);
1783
1784 // D0 -> a1(1420)+ {f0 pi0+[S]} pi-
1785 amplitude += fitpara[8] * (SF_Ap_P_SP * RBW_a11420_p * FT_f0980_00);
1786
1787 // D0 -> a1(1640)+ {rho pi[S]} pi-
1788 amplitude += fitpara[9] * (SF_Ap_S_Vp1P * RBW_a11640_p * GS_rho770_p1 + SF_Ap_S_Vp2P * RBW_a11640_p * GS_rho770_p2);
1789
1790 // D0 -> a1(1640)- {rho pi[S]} pi+
1791 amplitude += fitpara[10] * (SF_Am_S_Vm1P * RBW_a11640_m * GS_rho770_m1 + SF_Am_S_Vm2P * RBW_a11640_m * GS_rho770_m2);
1792
1793 // D0 -> a2(1320)+ {rho+ pi0} pi-
1794 double SF_Tp_D_Vp1P = contract_22_0(contract_22_2(contract_31_2(contract_41_3(epsilon_uvmn, contract_21_1(Proj1_3p, T1_PipPi01)),
1795 PipPi01Pi02), contract_42_2(Proj2_3p, T2_3pm)), T2_PipPi01Pi02);
1796 double SF_Tp_D_Vp2P = contract_22_0(contract_22_2(contract_31_2(contract_41_3(epsilon_uvmn, contract_21_1(Proj1_3p, T1_PipPi02)),
1797 PipPi01Pi02), contract_42_2(Proj2_3p, T2_3pm)), T2_PipPi02Pi01);
1798
1799 amplitude += fitpara[11] * (SF_Tp_D_Vp1P * GS_rho770_p1 * RBW_a21320_p + SF_Tp_D_Vp2P * GS_rho770_p2 * RBW_a21320_p);
1800
1801 // D0 -> a2(1320)- {rho- pi0} pi+
1802 double SF_Tm_D_Vm1P = contract_22_0(contract_22_2(contract_31_2(contract_41_3(epsilon_uvmn, contract_21_1(Proj1_3m, T1_PimPi01)),
1803 PimPi01Pi02), contract_42_2(Proj2_3m, T2_3mp)), T2_PimPi01Pi02);
1804 double SF_Tm_D_Vm2P = contract_22_0(contract_22_2(contract_31_2(contract_41_3(epsilon_uvmn, contract_21_1(Proj1_3m, T1_PimPi02)),
1805 PimPi01Pi02), contract_42_2(Proj2_3m, T2_3mp)), T2_PimPi02Pi01);
1806 amplitude += fitpara[12] * (SF_Tm_D_Vm1P * GS_rho770_m1 * RBW_a21320_m + SF_Tm_D_Vm2P * GS_rho770_m2 * RBW_a21320_m);
1807
1808 // D0 -> h1(1170)0 {rho pi0} pi0
1809 amplitude += fitpara[13] * (SF_A01_S_Vp1P * RBW_h11170_01 * GS_rho770_p1 + SF_A02_S_Vp2P * RBW_h11170_02 * GS_rho770_p2 -
1810 SF_A01_S_Vm1P * RBW_h11170_01 * GS_rho770_m1 - SF_A02_S_Vm2P * RBW_h11170_02 * GS_rho770_m2 - SF_A01_S_VzP * RBW_h11170_01 *
1811 GS_rho770_pm - SF_A02_S_VzP * RBW_h11170_02 * GS_rho770_pm);
1812
1813 // D0 -> pi(1300)- {rho(770)- pi0} pi+
1814 double SF_Pm_P_Vm1P = contract_11_0(T1_PimPi01, T1_PimPi01Pi02);
1815 double SF_Pm_P_Vm2P = contract_11_0(T1_PimPi02, T1_PimPi02Pi01);
1816
1817 amplitude += fitpara[14] * (SF_Pm_P_Vm1P * GS_rho770_m1 * RBW_pi1300_m + SF_Pm_P_Vm2P * GS_rho770_m2 * RBW_pi1300_m);
1818
1819 // D0 -> pi(1300)- {f2 pi-} pi+
1820// double SF_Pm_D_TP = contract_22_0(T2_Pi01Pi02,T2_Pi01Pi02Pim);
1821// amplitude += fitpara[14]*fitpara[13]*(SF_Pm_D_TP*RBW_f21270_00*RBW_pi1300_m);
1822
1823 // D0 -> pi(1300)- {f0 pi-} pi+
1824 amplitude += fitpara[15] * fitpara[14] * (RBW_pi1300_m * PiPiS_00_0);
1825// amplitude += fitpara[15]*fitpara[13]*(RBW_pi1300_m*PiPiS_00_5);
1826 amplitude += fitpara[16] * fitpara[14] * (RBW_pi1300_m * PiPiS_00_6);
1827
1828 // D0 -> pi(1300)+ {rho(770)+ pi0} pi-
1829 double SF_Pp_P_Vp1P = contract_11_0(T1_PipPi01, T1_PipPi01Pi02);
1830 double SF_Pp_P_Vp2P = contract_11_0(T1_PipPi02, T1_PipPi02Pi01);
1831
1832 amplitude += fitpara[17] * (SF_Pp_P_Vp1P * GS_rho770_p1 * RBW_pi1300_p + SF_Pp_P_Vp2P * GS_rho770_p2 * RBW_pi1300_p);
1833
1834 // D0 -> pi(1300)+ {f2 pi-} pi+
1835// double SF_Pp_D_TP = contract_22_0(T2_Pi01Pi02,T2_Pi01Pi02Pip);
1836// amplitude += fitpara[14]*fitpara[17]*(SF_Pp_D_TP*RBW_f21270_00*RBW_pi1300_p);
1837
1838 // D0 -> pi(1300)+ {f0 pi+} pi-
1839 amplitude += fitpara[15] * fitpara[17] * (RBW_pi1300_p * PiPiS_00_0);
1840// amplitude += fitpara[15]*fitpara[17]*(RBW_pi1300_p*PiPiS_00_5);
1841 amplitude += fitpara[16] * fitpara[17] * (RBW_pi1300_p * PiPiS_00_6);
1842
1843 // D0 -> pi(1300)0 {rho pi} pi0
1844 double SF_P01_P_Vp1P = contract_11_0(T1_PipPi01, T1_PipPi01Pim);
1845 double SF_P02_P_Vp2P = contract_11_0(T1_PipPi02, T1_PipPi02Pim);
1846 double SF_P01_P_Vm1P = contract_11_0(T1_PimPi01, T1_PimPi01Pip);
1847 double SF_P02_P_Vm2P = contract_11_0(T1_PimPi02, T1_PimPi02Pip);
1848
1849 amplitude += fitpara[18] * (SF_P01_P_Vp1P * RBW_pi1300_01 * GS_rho770_p1 + SF_P02_P_Vp2P * RBW_pi1300_02 * GS_rho770_p2 +
1850 SF_P01_P_Vm1P * RBW_pi1300_01 * GS_rho770_m1 + SF_P02_P_Vm2P * RBW_pi1300_02 * GS_rho770_m2);
1851
1852
1853// double SF_P01_D_TP = contract_22_0(T2_PipPim,T2_PipPimPi01);
1854// double SF_P02_D_TP = contract_22_0(T2_PipPim,T2_PipPimPi02);
1855// amplitude += fitpara[14]*fitpara[18]*(-1.0)*(SF_P01_D_TP*RBW_f21270_pm*RBW_pi1300_01 + SF_P02_D_TP*RBW_f21270_pm*RBW_pi1300_02);
1856
1857 amplitude += fitpara[15] * fitpara[18] * (-1.0) * (RBW_pi1300_01 * PiPiS_pm_0 + RBW_pi1300_02 * PiPiS_pm_0);
1858// amplitude += fitpara[15]*fitpara[18]*(-1.0)*(RBW_pi1300_01*PiPiS_pm_5 + RBW_pi1300_02*PiPiS_pm_5);
1859 amplitude += fitpara[16] * fitpara[18] * (-1.0) * (RBW_pi1300_01 * PiPiS_pm_6 + RBW_pi1300_02 * PiPiS_pm_6);
1860
1861 // D0 -> pi2(1670)0[f2(1270) pi0] pi0
1862 double SF_PT01_S_TP = contract_22_0(contract_42_2(Proj2_3z1, T2_PipPim), T2_3z12);
1863 double SF_PT02_S_TP = contract_22_0(contract_42_2(Proj2_3z2, T2_PipPim), T2_3z21);
1864
1865 amplitude += fitpara[19] * (-1.0) * (SF_PT01_S_TP * RBW_f21270_pm * RBW_pi21670_01 + SF_PT02_S_TP * RBW_f21270_pm * RBW_pi21670_02);
1866
1867 // D0 -> rho+ rho- [S]
1868 double SF_Vp1Vm2_S = contract_11_0(T1_PipPi01, T1_PimPi02);
1869 double SF_Vp2Vm1_S = contract_11_0(T1_PipPi02, T1_PimPi01);
1870
1871 amplitude += fitpara[20] * (SF_Vp1Vm2_S * GS_rho770_p1 * GS_rho770_m2 + SF_Vp2Vm1_S * GS_rho770_p2 * GS_rho770_m1);
1872
1873 // D0 -> rho+ rho- [P]
1874 double SF_Vp1Vm2_P = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, T1_PipPi01), T1_PimPi02), T1_2p1m2), D0);
1875 double SF_Vp2Vm1_P = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, T1_PipPi02), T1_PimPi01), T1_2p2m1), D0);
1876
1877 amplitude += fitpara[21] * (SF_Vp1Vm2_P * GS_rho770_p1 * GS_rho770_m2 + SF_Vp2Vm1_P * GS_rho770_p2 * GS_rho770_m1);
1878
1879 // D0 -> rho+ rho- [D]
1880 double SF_Vp1Vm2_D = contract_11_0(contract_21_1(T2_2p1m2, T1_PipPi01), T1_PimPi02);
1881 double SF_Vp2Vm1_D = contract_11_0(contract_21_1(T2_2p2m1, T1_PipPi02), T1_PimPi01);
1882 amplitude += fitpara[22] * (SF_Vp1Vm2_D * GS_rho770_p1 * GS_rho770_m2 + SF_Vp2Vm1_D * GS_rho770_p2 * GS_rho770_m1);
1883
1884 amplitude += fitpara[23] * (SF_Vp1Vm2_D * GS_rho770_p1 * GS_rho1450_m2 + SF_Vp2Vm1_D * GS_rho770_p2 * GS_rho1450_m1);
1885
1886 // D0 -> rho0 (Pi0 Pi0)S
1887 double SF_VpmS12_P = contract_11_0(T1_PipPim, T1_2pm12);
1888
1889 amplitude += fitpara[24] * (SF_VpmS12_P * GS_rho770_pm * PiPiS_00_0);
1890 amplitude += fitpara[25] * (SF_VpmS12_P * GS_rho770_pm * PiPiS_00_5);
1891 amplitude += fitpara[26] * (SF_VpmS12_P * GS_rho770_pm * PiPiS_00_6);
1892
1893 // D0 -> f0f0
1894 //00
1895 amplitude += fitpara[27] * (PiPiS_pm_0 * PiPiS_00_0 + PiPiS_00_0 * PiPiS_pm_0);
1896 amplitude += fitpara[28] * (PiPiS_pm_0 * PiPiS_00_1 + PiPiS_00_0 * PiPiS_pm_1);
1897 amplitude += fitpara[29] * (PiPiS_pm_1 * PiPiS_00_5 + PiPiS_00_1 * PiPiS_pm_5);
1898 amplitude += fitpara[30] * (PiPiS_pm_5 * PiPiS_00_5 + PiPiS_00_5 * PiPiS_pm_5);
1899 amplitude += fitpara[31] * (PiPiS_pm_5 * PiPiS_00_6 + PiPiS_00_5 * PiPiS_pm_6);
1900
1901 // D0 -> f2(1270) f0
1902 double SF_TpmS00_D = contract_22_0(T2_PipPim, T2_2pm12);
1903 double SF_T00Spm_D = contract_22_0(T2_Pi01Pi02, T2_2pm12);
1904
1905 amplitude += fitpara[32] * (SF_TpmS00_D * RBW_f21270_pm * PiPiS_00_5 + SF_T00Spm_D * RBW_f21270_00 * PiPiS_pm_5);
1906 amplitude += fitpara[33] * (SF_TpmS00_D * RBW_f21270_pm * PiPiS_00_6 + SF_T00Spm_D * RBW_f21270_00 * PiPiS_pm_6);
1907
1908
1909 // D -> 1--(rho pi) pi0
1910 double SF_V1_Vz = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, PipPimPi01), T1_PipPimPi01), T1_PipPim),
1911 contract_21_1(Proj1_3z1, T1_3z12));
1912 double SF_V1_Vp1 = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, PipPimPi01), T1_PipPi01Pim), T1_PipPi01),
1913 contract_21_1(Proj1_3z1, T1_3z12));
1914 double SF_V1_Vm1 = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, PipPimPi01), T1_PimPi01Pip), T1_PimPi01),
1915 contract_21_1(Proj1_3z1, T1_3z12));
1916
1917 double SF_V2_Vz = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, PipPimPi02), T1_PipPimPi02), T1_PipPim),
1918 contract_21_1(Proj1_3z2, T1_3z21));
1919 double SF_V2_Vp2 = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, PipPimPi02), T1_PipPi02Pim), T1_PipPi02),
1920 contract_21_1(Proj1_3z2, T1_3z21));
1921 double SF_V1_Vm2 = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, PipPimPi02), T1_PimPi02Pip), T1_PimPi02),
1922 contract_21_1(Proj1_3z2, T1_3z21));
1923
1924
1925 // D0 -> omega(rho pi) pi0
1926 amplitude += (-1.0) * fitpara[34] * (SF_V1_Vp1 * RBW_omega_01 * GS_rho770_p1 - SF_V1_Vz * RBW_omega_01 * GS_rho770_pm - SF_V1_Vm1 *
1927 RBW_omega_01 * GS_rho770_m1 + SF_V2_Vp2 * RBW_omega_02 * GS_rho770_p2 - SF_V2_Vz * RBW_omega_02 * GS_rho770_pm - SF_V1_Vm2 *
1928 RBW_omega_02 * GS_rho770_m2);
1929
1930
1931 // D0 -> phi(rho pi) pi0
1932 amplitude += (-1.0) * fitpara[35] * (SF_V1_Vp1 * RBW_phi_01 * GS_rho770_p1 - SF_V1_Vz * RBW_phi_01 * GS_rho770_pm - SF_V1_Vm1 *
1933 RBW_phi_01 * GS_rho770_m1 + SF_V2_Vp2 * RBW_phi_02 * GS_rho770_p2 - SF_V2_Vz * RBW_phi_02 * GS_rho770_pm - SF_V1_Vm2 * RBW_phi_02 *
1934 GS_rho770_m2);
1935
1936 return amplitude;
1937
1938 }
1939
1940 int EvtD0Topippim2pi0::CalAmp()
1941 {
1942
1943 m_AmpD0 = CalD0Amp();
1944 m_AmpDb = CalDbAmp();
1945
1946 return 0;
1947 }
1948
1949 double EvtD0Topippim2pi0::mag2(std::complex<double> x)
1950 {
1951 double temp = x.real() * x.real() + x.imag() * x.imag();
1952 return temp;
1953 }
1954
1955 double EvtD0Topippim2pi0::arg(std::complex<double> x)
1956 {
1957 double temp = atan(x.imag() / x.real());
1958 if (x.real() < 0) temp = temp + TMath::Pi();
1959 return temp;
1960 }
1961 double EvtD0Topippim2pi0::Get_strongPhase()
1962 {
1963 double temp = arg(m_AmpD0) - arg(m_AmpDb);
1964 while (temp < -TMath::Pi()) {
1965 temp += 2.0 * TMath::Pi();
1966 }
1967 while (temp > TMath::Pi()) {
1968 temp -= 2.0 * TMath::Pi();
1969 }
1970 return temp;
1971 }
1972
1973 double EvtD0Topippim2pi0::AmplitudeSquare(int Charm, int Tagmode)
1974 {
1975
1976 EvtVector4R dp1 = GetDaugMomLab(0), dp2 = GetDaugMomLab(1), dp3 = GetDaugMomLab(2), dp4 = GetDaugMomLab(3); // pi+ pi- pi0 pi0
1977 EvtVector4R mp = dp1 + dp2 + dp3 + dp4;
1978
1979 double emp = mp.get(0);
1980 EvtVector3R boostp3(-mp.get(1) / emp, -mp.get(2) / emp, -mp.get(3) / emp);
1981
1982 EvtVector4R dp1bst = boostTo(dp1, boostp3);
1983 EvtVector4R dp2bst = boostTo(dp2, boostp3);
1984 EvtVector4R dp3bst = boostTo(dp3, boostp3);
1985 EvtVector4R dp4bst = boostTo(dp4, boostp3);
1986
1987 double p4pip[4], p4pim[4], p4pi01[4], p4pi02[4];
1988 for (int i = 0; i < 3; i++) {
1989 p4pip[i] = dp1bst.get(i + 1);
1990 p4pim[i] = dp2bst.get(i + 1);
1991 p4pi01[i] = dp3bst.get(i + 1);
1992 p4pi02[i] = dp4bst.get(i + 1);
1993 }
1994 p4pip[3] = dp1bst.get(0);
1995 p4pim[3] = dp2bst.get(0);
1996 p4pi01[3] = dp3bst.get(0);
1997 p4pi02[3] = dp4bst.get(0);
1998
1999 setInput(p4pip, p4pim, p4pi01, p4pi02);
2000 CalAmp();
2001 std::complex<double> ampD0, ampDb;
2002 if (Charm > 0) {
2003 ampD0 = Get_AmpD0();
2004 ampDb = Get_AmpDb();
2005 } else {
2006 ampD0 = Get_AmpDb();
2007 ampDb = Get_AmpD0();
2008 }
2009
2010 double ampsq = 1e-20;
2011 double r_tag = 0, R_tag = 0, delta_tag = 0;
2012
2013 if (Tagmode == 1 || Tagmode == 2 || Tagmode == 3) {
2014 if (Tagmode == 1) { // Kpi
2015 r_tag = 0.0586;
2016 R_tag = 1;
2017 delta_tag = 192.1 / 180.0 * 3.1415926;
2018 } else if (Tagmode == 2) { // Kpipi0
2019 r_tag = 0.0441;
2020 R_tag = 0.79;
2021 delta_tag = 196.0 / 180.0 * 3.1415926;
2022 } else if (Tagmode == 3) { // K3pi
2023 r_tag = 0.0550;
2024 R_tag = 0.44;
2025 delta_tag = 161.0 / 180.0 * 3.1415926;
2026 }
2027 std::complex<double> qcf(r_tag * R_tag * cos(-delta_tag), r_tag * R_tag * sin(-delta_tag));
2028 std::complex<double> ampD0_part1 = ampD0 - qcf * ampDb;
2029 ampsq = mag2(ampD0_part1) + r_tag * r_tag * (1 - R_tag * R_tag) * (mag2(ampDb));
2030 } else {
2031 ampsq = mag2(ampD0);
2032 }
2033 return (ampsq <= 0) ? 1e-20 : ampsq;
2034 }
2035
2037}
Belle2::atan
double atan(double a)
atan for double
Definition beamHelpers.h:34
B2_EVTGEN_REGISTER_MODEL
#define B2_EVTGEN_REGISTER_MODEL(classname)
Class to register B2_EVTGEN_REGISTER_MODEL.
Definition EvtGenModelRegister.h:75
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
Belle2
Abstract base class for different kinds of events.
Definition MillepedeAlgorithm.h:17