BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
D0Topippim2pi0.cxx
Go to the documentation of this file.
1#include "QCMCFilterAlg/D0Topippim2pi0.h"
2#include <stdlib.h>
3#include <iostream>
4#include <string>
5#include <complex>
6#include <vector>
7#include <math.h>
8#include "TMath.h"
9using namespace std; //::endl;
10
11D0Topippim2pi0::~D0Topippim2pi0() {}
12
13void D0Topippim2pi0::init(){
14
15 //std::cout << "D0Topippim2pi0 ==> Initialization !" << std::endl;
16
17 double mag[35], pha[35];
18 mag[0]= 100.0; pha[0]= 0.0;
19 mag[1]= 7.95507; pha[1]= -0.0687407;
20 mag[2]= 37.5559; pha[2]= -1.74946;
21 mag[3]= 61.2172; pha[3]= 2.98079;
22 mag[4]= 187.79; pha[4]= 2.64471;
23 mag[5]= 385.474; pha[5]= -0.137107;
24 mag[6]= 0.330788; pha[6]= 0.268133;
25 mag[7]= 0.584175; pha[7]= -2.89693;
26 mag[8]= 127.158; pha[8]= -2.47773;
27 mag[9]= 339.914; pha[9]= 2.22856;
28 mag[10]=0.320888; pha[10]=-2.6194;
29 mag[11]=0.366283; pha[11]=-0.26867;
30 mag[12]=14.1344; pha[12]=-0.41164;
31 mag[13]=86.0865; pha[13]=-2.49649;
32 mag[14]=6.1541; pha[14]=-1.18299;
33 mag[15]=56.6067; pha[15]=0.142977;
34 mag[16]=92.3073; pha[16]=-2.15881;
35 mag[17]=80.9453; pha[17]=0.825815;
36 mag[18]=16.9555; pha[18]=-2.98994;
37 mag[19]=9.72524; pha[19]=-1.39929;
38 mag[20]=5.71448; pha[20]=0.271902;
39 mag[21]=21.4195; pha[21]=-1.23701;
40 mag[22]=56.8867; pha[22]=-0.385837;
41 mag[23]=231.626; pha[23]=2.14842;
42 mag[24]=2938.45; pha[24]=-0.693491;
43 mag[25]=7252.7; pha[25]=2.23659;
44 mag[26]=5165.87; pha[26]=0.913557;
45 mag[27]=11508.6; pha[27]=-1.07187;
46 mag[28]=2461.86; pha[28]=1.8709;
47 mag[29]=8757.75; pha[29]=2.40756;
48 mag[30]=19.7413; pha[30]=-1.0753;
49 mag[31]=66.3826; pha[31]=2.34666;
50 mag[32]=11.2904; pha[32]=-0.822345;
51 mag[33]=2.04576; pha[33]=-0.281429;
52 mag[34]=0.57927; pha[34]=2.7182;
53
54 fitpara.clear();
55 for(int i=0; i<35; i++){
56 complex<double> ctemp(mag[i]*cos(pha[i]),mag[i]*sin(pha[i]));
57 fitpara.push_back(ctemp);
58 }
59
60 g_uv.clear();
61 for(int i=0; i<4; i++){
62 for(int j=0; j<4; j++){
63 if(i!=j){
64 g_uv.push_back(0.0);
65 }else if(i<3){
66 g_uv.push_back(-1.0);
67 }else if(i==3){
68 g_uv.push_back(1.0);
69 }
70 }
71 }
72
73 epsilon_uvmn.clear();
74 for(int i=0; i<4; i++){
75 for(int j=0; j<4; j++){
76 for(int k=0; k<4; k++){
77 for(int l=0; l<4; l++){
78 if(i==j || i==k || i==l || j==k || j==l || k==l){
79 epsilon_uvmn.push_back(0.0);
80 }else{
81 if(i==0 && j==1 && k==2 && l==3) epsilon_uvmn.push_back(1.0);
82 if(i==0 && j==1 && k==3 && l==2) epsilon_uvmn.push_back(-1.0);
83 if(i==0 && j==2 && k==1 && l==3) epsilon_uvmn.push_back(-1.0);
84 if(i==0 && j==2 && k==3 && l==1) epsilon_uvmn.push_back(1.0);
85 if(i==0 && j==3 && k==1 && l==2) epsilon_uvmn.push_back(1.0);
86 if(i==0 && j==3 && k==2 && l==1) epsilon_uvmn.push_back(-1.0);
87
88 if(i==1 && j==0 && k==2 && l==3) epsilon_uvmn.push_back(-1.0);
89 if(i==1 && j==0 && k==3 && l==2) epsilon_uvmn.push_back(1.0);
90 if(i==1 && j==2 && k==0 && l==3) epsilon_uvmn.push_back(1.0);
91 if(i==1 && j==2 && k==3 && l==0) epsilon_uvmn.push_back(-1.0);
92 if(i==1 && j==3 && k==0 && l==2) epsilon_uvmn.push_back(-1.0);
93 if(i==1 && j==3 && k==2 && l==0) epsilon_uvmn.push_back(1.0);
94
95 if(i==2 && j==0 && k==1 && l==3) epsilon_uvmn.push_back(1.0);
96 if(i==2 && j==0 && k==3 && l==1) epsilon_uvmn.push_back(-1.0);
97 if(i==2 && j==1 && k==0 && l==3) epsilon_uvmn.push_back(-1.0);
98 if(i==2 && j==1 && k==3 && l==0) epsilon_uvmn.push_back(1.0);
99 if(i==2 && j==3 && k==0 && l==1) epsilon_uvmn.push_back(1.0);
100 if(i==2 && j==3 && k==1 && l==0) epsilon_uvmn.push_back(-1.0);
101
102 if(i==3 && j==0 && k==1 && l==2) epsilon_uvmn.push_back(-1.0);
103 if(i==3 && j==0 && k==2 && l==1) epsilon_uvmn.push_back(1.0);
104 if(i==3 && j==1 && k==0 && l==2) epsilon_uvmn.push_back(1.0);
105 if(i==3 && j==1 && k==2 && l==0) epsilon_uvmn.push_back(-1.0);
106 if(i==3 && j==2 && k==0 && l==1) epsilon_uvmn.push_back(-1.0);
107 if(i==3 && j==2 && k==1 && l==0) epsilon_uvmn.push_back(1.0);
108
109 }
110 }
111 }
112 }
113 }
114
115 _nd = 4;
116 math_pi = 3.1415926f;
117 mass_Pion = 0.13957f;
118
119 rRes = 3.0*0.197321;
120 rD = 5.0*0.197321;
121 m_Pi = 0.139570;
122 m_Pi0 = 0.134977;
123 m2_Pi = m_Pi*m_Pi;
124 m2_Pi0 = m_Pi0*m_Pi0;
125
126 m0_rho7700 = 0.77526;
127 w0_rho7700 = 0.1478;
128
129 m0_rho770p = 0.77511;
130 w0_rho770p = 0.1491;
131
132 m0_rho1450 = 1.465;
133 w0_rho1450 = 0.400;
134
135 m0_f21270 = 1.2755;
136 w0_f21270 = 0.1867;
137
138 m0_a11260 = 1.1337;
139 g1_a11260 = 0.00335;
140 g2_a11260 = 0.0;
141
142 m0_pi1300 = 1.498;
143 w0_pi1300 = 0.590;
144
145 m0_a11420 = 1.411;
146 w0_a11420 = 0.161;
147
148 m0_a11640 = 1.655;
149 w0_a11640 = 0.254;
150
151 m0_a21320 = 1.3186;
152 w0_a21320 = 0.105;
153
154 m0_pi21670 = 1.6706;
155 w0_pi21670 = 0.258;
156
157 m0_pi11400 = 1.354;
158 w0_pi11400 = 0.330;
159
160 m0_h11170 = 1.166;
161 w0_h11170 = 0.375;
162
163 m0_omega = 0.78265;
164 w0_omega = 0.00849;
165
166 m0_phi = 1.019461;
167 w0_phi = 0.004249;
168
169 s0_prod = -5.0;
170
171}
172
173vector<double> D0Topippim2pi0::sum_tensor(vector<double> pa, vector<double> pb)
174{
175 if(pa.size()!=pb.size()){
176 cout<<"error sum tensor"<<endl;
177 exit(0);
178 }
179 vector<double> temp; temp.clear();
180 for(int i=0; i<pa.size(); i++){
181 double sum = pa[i] + pb[i];
182 temp.push_back(sum);
183 }
184 return temp;
185}
186
187double D0Topippim2pi0::contract_11_0(vector<double> pa, vector<double> pb){
188 if(pa.size()!=pb.size() || pa.size()!=4) {
189 cout<<"error contract 11->0"<<endl;
190 exit(0);
191 }
192 double temp = pa[3]*pb[3] - pa[0]*pb[0] - pa[1]*pb[1] - pa[2]*pb[2];
193 return temp;
194
195}
196
197vector<double> D0Topippim2pi0::contract_21_1(vector<double> pa, vector<double> pb){
198 if(pa.size()!=16 || pb.size()!=4) {
199 cout<<"error contract 21->1"<<endl;
200 exit(0);
201 }
202 vector<double> temp; temp.clear();
203 for(int i=0; i<4; i++){
204 double sum = 0;
205 for(int j=0; j<4; j++){
206 int idx = i*4+j;
207 sum += pa[idx]*pb[j]*g_uv[4*j+j];
208 }
209 temp.push_back(sum);
210 }
211 return temp;
212
213}
214
215double D0Topippim2pi0::contract_22_0(vector<double> pa, vector<double> pb){
216 if(pa.size()!=pb.size() || pa.size()!=16) {
217 cout<<"error contract 22->0"<<endl;
218 exit(0);
219 }
220 double temp = 0;
221 for(int i=0; i<4; i++){
222 for(int j=0; j<4; j++){
223 int idx = i*4+j;
224 temp += pa[idx]*pb[idx]*g_uv[4*i+i]*g_uv[4*j+j];
225 }
226 }
227 return temp;
228
229}
230
231vector<double> D0Topippim2pi0::contract_31_2(vector<double> pa, vector<double> pb){
232 if(pa.size()!=64 || pb.size()!=4) {
233 cout<<"error contract 31->2"<<endl;
234 exit(0);
235 }
236 vector<double> temp; temp.clear();
237 for(int i=0; i<16; i++){
238 double sum = 0;
239 for(int j=0; j<4; j++){
240 int idx = i*4+j;
241 sum += pa[idx]*pb[j]*g_uv[4*j+j];
242 }
243 temp.push_back(sum);
244 }
245 return temp;
246
247}
248
249vector<double> D0Topippim2pi0::contract_41_3(vector<double> pa, vector<double> pb){
250 if(pa.size()!=256|| pb.size()!=4) {
251 cout<<"error contract 41->3"<<endl;
252 exit(0);
253 }
254 vector<double> temp; temp.clear();
255 for(int i=0; i<64; i++){
256 double sum = 0;
257 for(int j=0; j<4; j++){
258 int idx = i*4+j;
259 sum += pa[idx]*pb[j]*g_uv[4*j+j];
260 }
261 temp.push_back(sum);
262 }
263 return temp;
264
265}
266
267vector<double> D0Topippim2pi0::contract_42_2(vector<double> pa, vector<double> pb){
268 if(pa.size()!=256|| pb.size()!=16) {
269 cout<<"error contract 42->2"<<endl;
270 exit(0);
271 }
272 vector<double> temp; temp.clear();
273 for(int i=0; i<16; i++){
274 double sum = 0;
275 for(int j=0; j<4; j++){
276 for(int k=0; k<4; k++){
277 int idxa = i*16+j*4+k;
278 int idxb = j*4+k;
279 sum += pa[idxa] * pb[idxb] * g_uv[4*j+j] * g_uv[4*k+k];
280 }
281 }
282 temp.push_back(sum);
283 }
284
285 return temp;
286
287}
288vector<double> D0Topippim2pi0::contract_22_2(vector<double> pa, vector<double> pb){
289 if(pa.size()!=16|| pb.size()!=16) {
290 cout<<"error contract 42->2"<<endl;
291 exit(0);
292 }
293 vector<double> temp; temp.clear();
294 for(int i=0; i<4; i++){
295 for(int j=0; j<4; j++){
296 double sum = 0;
297 for(int k=0; k<4; k++){
298 int idxa = i*4+k;
299 int idxb = j*4+k;
300 sum += pa[idxa] * pb[idxb] * g_uv[4*k+k];
301 }
302 temp.push_back(sum);
303 }
304 }
305
306 return temp;
307
308}
309
310/*
311 vector<double> D0Topippim2pi0::contract_11_2(vector<double> pa, vector<double> pb){
312 if(pa.size()!=pb.size() || pa.size()!=4) {
313 cout<<"error contract 11->2"<<endl;
314 exit(0);
315 }
316 vector<double> temp; temp.clear();
317 for(int i=0; i<4; i++){
318 for(int j=0; j<4; j++){
319 double prod = pa[i]*pb[j];
320 temp.push_back(prod);
321 }
322 }
323 return temp;
324 }
325
326 vector<double> D0Topippim2pi0::contract_22_4(vector<double> pa, vector<double> pb){
327 if(pa.size()!=pb.size() || pa.size()!=16) {
328 cout<<"error contract 22->4"<<endl;
329 exit(0);
330 }
331 vector<double> temp; temp.clear();
332 for(int i=0; i<16; i++){
333 for(int j=0; j<16; j++){
334 double prod = pa[i]*pb[j];
335 temp.push_back(prod);
336 }
337 }
338 return temp;
339 }
340 */
341//OrbitalTensors
342vector<double> D0Topippim2pi0::OrbitalTensors(vector<double> pa, vector<double> pb, vector<double> pc, double r, int rank)
343{
344 if(pa.size()!=4 || pb.size()!=4 || pc.size()!=4) {
345 cout<<"Error: pa, pb, pc"<<endl;
346 exit(0);
347 }
348 if(rank<0){
349 cout<<"Error: L<0 !!!"<<endl;
350 exit(0);
351 }
352
353 // relative momentum
354 vector<double> mr; mr.clear();
355
356 for(int i=0; i<4; i++){
357 double temp = pb[i] - pc[i];
358 mr.push_back(temp);
359 }
360
361 // "Masses" of particles
362 double msa = contract_11_0(pa, pa);
363 double msb = contract_11_0(pb, pb);
364 double msc = contract_11_0(pc, pc);
365
366 // Q^2_abc
367 double top = msa + msb - msc;
368 double Q2abc = top*top/(4.0*msa) - msb;
369
370 // Barrier factors
371 double Q_0 = 0.197321f/r;
372 double Q_02 = Q_0*Q_0;
373 double Q_04 = Q_02*Q_02;
374 // double Q_06 = Q_04*Q_02;
375 // double Q_08 = Q_04*Q_04;
376
377 double Q4abc = Q2abc*Q2abc;
378 // double Q6abc = Q4abc*Q2abc;
379 // double Q8abc = Q4abc*Q4abc;
380
381 double mB1 = sqrt(2.0f/(Q2abc + Q_02));
382 double mB2 = sqrt(13.0f/(Q4abc + (3.0f*Q_02)*Q2abc + 9.0f*Q_04));
383 // mB3 = &sqrt(277.0f/(Q6abc + (6.0f*Q_02)*Q4abc + (45.0f*Q_04)*Q2abc + 225.0f*Q_06));
384 // mB4 = &sqrt(12746.0f/(Q8abc + (10.0f*Q_02)*Q6abc + (135.0f*Q_04)*Q4abc + (1575.0f*Q_06)*Q2abc + 11025.0f*Q_08));
385
386 // Projection Operator 2-rank
387 vector<double> proj_uv; proj_uv.clear();
388 for(int i=0; i<4; i++){
389 for(int j=0; j<4; j++){
390 int idx = i*4 + j;
391 double temp = -g_uv[idx] + pa[i]*pa[j]/msa;
392 proj_uv.push_back(temp);
393 }
394 }
395
396 // Orbital Tensors
397 if(rank==0)
398 {
399 vector<double> t; t.clear();
400 t.push_back(1.0);
401 return t;
402
403 }else if(rank<3)
404 {
405 vector<double> t_u; t_u.clear();
406 vector<double> Bt_u; Bt_u.clear();
407 for(int i=0; i<4; i++){
408 double temp = 0;
409 for(int j=0; j<4; j++){
410 int idx = i*4 + j;
411 temp += -proj_uv[idx]*mr[j]*g_uv[j*4+j];
412 }
413 t_u.push_back(temp);
414 Bt_u.push_back(temp*mB1);
415 }
416 if(rank==1) return Bt_u;
417
418 double t_u2 = contract_11_0(t_u,t_u);
419
420 vector<double> Bt_uv; Bt_uv.clear();
421 for(int i=0; i<4; i++){
422 for(int j=0; j<4; j++){
423 int idx = 4*i + j;
424 double temp = t_u[i]*t_u[j] + (1.0/3.0)*proj_uv[idx]*t_u2;
425 Bt_uv.push_back(temp*mB2);
426 }
427 }
428 if(rank==2) return Bt_uv;
429
430 }else
431 {
432 cout<<"rank>2: please add it by yourself!!!"<<endl;
433 exit(0);
434 }
435}
436
437// projection Tensor
438vector<double> D0Topippim2pi0::ProjectionTensors(vector<double> pa, int rank)
439{
440 if(pa.size()!=4) {
441 cout<<"Error: pa"<<endl;
442 exit(0);
443 }
444 if(rank<0){
445 cout<<"Error: L<0 !!!"<<endl;
446 exit(0);
447 }
448
449 double msa = contract_11_0(pa, pa);
450
451 // Projection Operator 2-rank
452 vector<double> proj_uv; proj_uv.clear();
453 for(int i=0; i<4; i++){
454 for(int j=0; j<4; j++){
455 int idx = i*4 + j;
456 double temp = -g_uv[idx] + pa[i]*pa[j]/msa;
457 proj_uv.push_back(temp);
458 }
459 }
460
461 // Orbital Tensors
462 if(rank==0)
463 {
464 vector<double> t; t.clear();
465 t.push_back(1.0);
466 return t;
467
468 }else if(rank==1)
469 {
470 return proj_uv;
471 }else if(rank==2)
472 {
473 vector<double> proj_uvmn; proj_uvmn.clear();
474 for(int i=0; i<4; i++){
475 for(int j=0; j<4; j++){
476 for(int k=0; k<4; k++){
477 for(int l=0; l<4; l++){
478
479 int idx1_1 = 4*i + k;
480 int idx1_2 = 4*i + l;
481 int idx1_3 = 4*i + j;
482
483 int idx2_1 = 4*j + l;
484 int idx2_2 = 4*j + k;
485 int idx2_3 = 4*k + l;
486
487 double temp = (1.0/2.0)*(proj_uv[idx1_1]*proj_uv[idx2_1] + proj_uv[idx1_2]*proj_uv[idx2_2]) - (1.0/3.0)*proj_uv[idx1_3]*proj_uv[idx2_3];
488 proj_uvmn.push_back(temp);
489 }
490 }
491 }
492 }
493 return proj_uvmn;
494
495 }else
496 {
497 cout<<"rank>2: please add it by yourself!!!"<<endl;
498 exit(0);
499 }
500}
501double D0Topippim2pi0::fundecaymomentum(double mr2, double m1_2, double m2_2){
502 double mr = sqrt(mr2);
503 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;
504 double ret = sqrt(poly)/(2*mr);
505 if(poly<0)
506 //if(sqrt(m1_2) +sqrt(m2_2) > mr)
507 ret = 0.0f;
508 return ret;
509}
510
511complex<double> D0Topippim2pi0::breitwigner(double mx2, double mr, double wr)
512{
513 double output_x = 0;
514 double output_y = 0;
515
516 double mr2 = mr*mr;
517 double diff = mr2-mx2;
518 double denom = diff*diff + wr*wr*mr2;
519 if (wr<0) {
520 output_x = 0;
521 output_y = 0;
522 }
523 else if (wr<10) {
524 output_x = diff/denom;
525 output_y = wr*mr/denom;
526 }
527 else { /* phase space */
528 output_x = 1;
529 output_y = 0;
530 }
531 complex<double> output(output_x,output_y);
532 return output;
533}
534
535/* GS propagator */
536double D0Topippim2pi0::h(double m, double q){
537 double h = 2.0/math_pi*q/m*log((m+2.0*q)/(2.0*mass_Pion));
538 return h;
539}
540
541double D0Topippim2pi0::dh(double m0, double q0){
542 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);
543 return dh;
544}
545
546double D0Topippim2pi0::f(double m0, double sx, double q0, double q){
547 double m = sqrt(sx);
548 double f = m0*m0/(q0*q0*q0)*(q*q*(h(m,q)-h(m0,q0))+(m0*m0-sx)*q0*q0*dh(m0,q0));
549 return f;
550}
551
552double D0Topippim2pi0::d(double m0, double q0){
553 double d = 3.0/math_pi*mass_Pion*mass_Pion/(q0*q0)*log((m0+2.0*q0)/(2.0*mass_Pion)) + m0/(2.0*math_pi*q0) - (mass_Pion*mass_Pion*m0)/(math_pi*q0*q0*q0);
554 return d;
555}
556
557double D0Topippim2pi0::fundecaymomentum2(double mr2, double m1_2, double m2_2){
558 double mr = sqrt(mr2);
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
566double D0Topippim2pi0::wid(double mass, double sa, double sb, double sc, double r, int l){
567 double widm = 1.0;
568 double sa0 = mass*mass;
569 double m = sqrt(sa);
570 double q = fundecaymomentum2(sa,sb,sc);
571 double q0 = fundecaymomentum2(sa0,sb,sc);
572 double z = q*r*r;
573 double z0 = q0*r*r;
574 double F = 0.0;
575 if(l == 0) F = 1.0;
576 if(l == 1) F = sqrt((1.0+z0)/(1.0+z));
577 if(l == 2) F = sqrt((9.0+3.0*z0+z0*z0)/(9.0+3.0*z+z*z));
578 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));
579 if(l == 4) F = sqrt((11025.0+1575.0*z0+135.0*z0*z0+10.0*z0*z0*z0+z0*z0*z0*z0)/(11025.0+1575.0*z+135.0*z*z+10.0*z*z*z+z*z*z*z));
580 double t = sqrt(q/q0);
581 //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);
582 uint i=0;
583 for(i=0; i<(2*l+1); i++) {
584 widm *= t;
585 }
586 widm *= (mass/m*F*F);
587 return widm;
588}
589
590/* for rho0, use GS, rho0->pi+ pi-, only angular momentum 1*/
591complex<double> D0Topippim2pi0::GS(double mx2, double mr, double wr, double m1_2, double m2_2, double r, int l){
592
593 double mr2 = mr*mr;
594 double q = fundecaymomentum(mx2, m1_2, m2_2);
595 double q0 = fundecaymomentum(mr2, m1_2, m2_2);
596 double numer = 1.0+d(mr,q0)*wr/mr;
597 double denom_real = mr2-mx2+wr*f(mr,mx2,q0,q);
598 double denom_imag = mr*wr*wid(mr,mx2,m1_2,m2_2,r,l);//real-i*imag;
599
600 double denom = denom_real*denom_real+denom_imag*denom_imag;
601 double output_x = denom_real*numer/denom;
602 double output_y = denom_imag*numer/denom;
603
604 complex<double> output(output_x,output_y);
605 return output;
606}
607
608complex<double> D0Topippim2pi0::RBW(double mx2, double mr, double wr, double m1_2, double m2_2, double r, int l){
609 double mx = sqrt(mx2);
610 double mr2 = mr*mr;
611 double denom_real = mr2-mx2;
612 double denom_imag = 0;
613 if(m1_2>0 && m2_2>0){
614 denom_imag = mr*wr*wid(mr,mx2,m1_2,m2_2,r,l);//real-i*imag;
615 }else{
616 denom_imag = mr*wr;
617 }
618
619 double denom = denom_real*denom_real+denom_imag*denom_imag;
620 double output_x = denom_real/denom;
621 double output_y = denom_imag/denom;
622
623 complex<double> output(output_x,output_y);
624 return output;
625}
626
627/* build a1260 propagator */
628double D0Topippim2pi0::widT1260(int i, double g1, double g2){
629
630 double wid1[300] = { 0.00100302, 0.0069383, 0.0223132, 0.0504984, 0.093998, 0.154569, 0.233464, 0.331844, 0.450141, 0.589068,
631 0.748192, 0.928578, 1.13001, 1.35227, 1.59548, 1.86005, 2.14633, 2.45252, 2.78199, 3.13055,
632 3.50351, 3.89773, 4.31274, 4.75409, 5.21133, 5.69991, 6.20735, 6.74638, 7.30128, 7.8858,
633 8.50289, 9.14654, 9.82395, 10.5209, 11.2643, 12.0436, 12.8585, 13.692, 14.598, 15.5291,
634 16.5158, 17.5337, 18.6289, 19.7599, 20.9847, 22.2557, 23.5959, 25.0095, 26.5123, 28.0789,
635 29.7542, 31.5143, 33.3769, 35.3462, 37.3911, 39.5988, 41.874, 44.2815, 46.7975, 49.401,
636 52.0553, 54.7753, 57.5932, 60.4542, 63.3049, 66.0665, 68.8987, 71.6282, 74.2613, 76.8713,
637 79.3528, 81.722, 84.1212, 86.227, 88.4243, 90.3478, 92.2478, 94.1483, 95.8541, 97.5086,
638 99.0092, 100.48, 101.861, 103.153, 104.338, 105.576, 106.696, 107.647, 108.761, 109.725,
639 110.625, 111.529, 112.426, 113.01, 113.877, 114.647, 115.086, 115.856, 116.533, 117.076,
640 117.646, 118.25, 118.653, 119.023, 119.554, 119.958, 120.384, 121.036, 121.402, 121.686,
641 122.44, 122.592, 122.979, 123.39, 123.819, 123.957, 124.459, 124.681, 125.071, 125.405,
642 125.769, 125.978, 126.542, 126.817, 127.017, 127.292, 127.765, 127.989, 128.542, 128.66,
643 128.923, 129.094, 129.441, 129.716, 130.23, 130.506, 130.658, 131.12, 131.308, 131.579,
644 131.994, 132.28, 132.594, 132.79, 133.107, 133.589, 133.935, 134.242, 134.484, 134.765,
645 135.208, 135.58, 135.922, 136.236, 136.545, 136.949, 137.216, 137.503, 137.994, 138.35,
646 138.62, 138.912, 139.413, 139.831, 140.137, 140.478, 141, 141.3, 141.807, 142.291,
647 142.864, 143.315, 143.678, 144.215, 144.587, 145.122, 145.8, 145.885, 146.583, 147.226,
648 147.661, 148.187, 148.698, 149.227, 149.832, 150.548, 151.122, 151.674, 152.074, 152.666,
649 153.295, 153.899, 154.661, 155.364, 155.908, 156.495, 157.36, 157.719, 158.533, 159.287,
650 159.79, 160.654, 161.257, 161.93, 162.437, 163.468, 163.957, 164.631, 165.414, 166.203,
651 166.738, 167.61, 168.453, 169.101, 170.111, 170.333, 171.123, 171.958, 173.018, 173.663,
652 174.213, 175.241, 175.579, 176.435, 177.291, 178.071, 178.969, 179.635, 180.118, 181.078,
653 182.007, 182.73, 183.282, 184.161, 184.981, 185.695, 186.506, 187.16, 187.996, 188.439,
654 189.416, 190.104, 190.759, 191.786, 192.331, 193.318, 193.836, 194.981, 195.634, 196.231,
655 196.832, 197.835, 198.608, 199.273, 199.854, 200.695, 201.719, 202.105, 202.958, 203.707,
656 204.306, 205.319, 205.977, 206.875, 207.687, 208.352, 209.04, 209.352, 210.313, 211.322,
657 212.02, 212.458, 213.246, 214.331, 214.923, 215.466, 216.536, 217.346, 217.867, 218.463,
658 219.201, 219.88, 220.829, 221.461, 222.399, 223.068, 223.712, 224.174, 224.837, 225.838,
659 227.019, 227.171, 227.797, 228.663, 229.429, 230.323, 230.845, 231.574, 232.417, 232.677 };
660 double wid2[300] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
661 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
662 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
663 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
664 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
665 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
666 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
667 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
668 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
669 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
670 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
671 1.87136e-06, 1.50063e-05, 5.10425e-05, 0.000122121, 0.000240853, 0.000420318, 0.000675161, 0.0010173, 0.00146434, 0.00203321,
672 0.00273489, 0.0035927, 0.00462579, 0.00584255, 0.00727372, 0.00895462, 0.0108831, 0.013085, 0.0156197, 0.0184865,
673 0.0217078, 0.0253423, 0.0294103, 0.0339191, 0.0389837, 0.0446351, 0.0508312, 0.0577268, 0.0653189, 0.0737049,
674 0.0829819, 0.0930611, 0.104328, 0.116663, 0.130105, 0.144922, 0.16122, 0.179091, 0.198759, 0.220133,
675 0.243916, 0.269803, 0.298861, 0.330061, 0.365741, 0.40437, 0.447191, 0.49501, 0.548576, 0.606445,
676 0.674414, 0.748353, 0.831686, 0.929938, 1.03771, 1.16187, 1.30387, 1.47341, 1.65629, 1.88318,
677 2.14353, 2.44169, 2.79831, 3.2009, 3.65522, 4.16317, 4.69597, 5.2585, 5.85965, 6.44984,
678 7.04202, 7.60113, 8.14571, 8.73195, 9.24537, 9.75717, 10.2093, 10.6731, 11.1487, 11.5819,
679 12.0158, 12.4253, 12.8113, 13.2073, 13.5995, 13.9317, 14.312, 14.6595, 14.9511, 15.2668,
680 15.6092, 15.9349, 16.1873, 16.5049, 16.819, 17.0743, 17.3621, 17.6094, 17.8418, 18.0681,
681 18.3141, 18.5914, 18.8187, 19.0562, 19.2282, 19.4918, 19.7326, 19.9112, 20.134, 20.3386,
682 20.511, 20.6865, 20.8958, 21.0518, 21.2967, 21.44, 21.6361, 21.8012, 21.9523, 22.1736,
683 22.2615, 22.4207, 22.6056, 22.7198, 22.9299, 23.0605, 23.2959, 23.3808, 23.4961, 23.6793,
684 23.7843, 23.9697, 24.0689, 24.1919, 24.405, 24.3898, 24.6018, 24.7294, 24.789, 24.9978,
685 25.0626, 25.1728, 25.2809, 25.3579, 25.5444, 25.5995, 25.7644, 25.8397, 25.9229, 26.095,
686 26.1495, 26.2899, 26.3871, 26.54, 26.6603, 26.7008, 26.7836, 26.907, 26.9653, 26.9969,
687 27.1226, 27.226, 27.3543, 27.4686, 27.4887, 27.6163, 27.6986, 27.7506, 27.7884, 27.8662,
688 27.9886, 28.0573, 28.1238, 28.2612, 28.3209, 28.3457, 28.4392, 28.5086, 28.6399, 28.7603,
689 28.788, 28.8502, 28.9038, 28.9667, 28.975, 29.0032, 29.2681, 29.2392, 29.2572, 29.3364 };
690
691 return wid1[i]*g1+wid2[i]*g2;
692}
693
694double D0Topippim2pi0::anywid1260(double sc, double g1, double g2){
695
696 double smin = (0.13957*3)*(0.13957*3);
697 double dh = 0.01;
698 int od = (sc - 0.18)/dh;
699 double sc_m = 0.18 + od*dh;
700 double widuse = 0;
701 if(sc>=0.18 && sc<=3.17){
702 widuse = ((sc-sc_m)/dh)*(widT1260(od+1,g1,g2)-widT1260(od,g1,g2))+widT1260(od,g1,g2);
703 }else if(sc<0.18 && sc>smin){
704 widuse = ((sc - smin)/(0.18-smin))*widT1260(0,g1,g2);
705 }else if(sc>3.17){
706 widuse = widT1260(299,g1,g2);
707 }else{
708 widuse = 0;
709 }
710 return widuse;
711
712}
713
714complex<double> D0Topippim2pi0::RBWa1260(double mx2, double mr, double g1, double g2){
715
716 double mx = sqrt(mx2);
717 double mr2 = mr*mr;
718 double wid0 = anywid1260(mx2,g1,g2);
719
720 double denom_real = mr2-mx2;
721 double denom_imag = mr*wid0;//real-i*imag;
722
723 double denom = denom_real*denom_real+denom_imag*denom_imag;
724 double output_x = denom_real/denom;
725 double output_y = denom_imag/denom;
726
727 complex<double> output(output_x,output_y);
728 return output;
729
730}
731
732/* build pi1300 propagator */
733double D0Topippim2pi0::widT1300(int i){
734 double wid1[300] = { 0.0702928, 0.399073, 0.991742, 1.82025, 2.85953, 4.08606, 5.48082, 7.02683, 8.70496, 10.5007,
735 12.4053, 14.4026, 16.4831, 18.6423, 20.8642, 23.1544, 25.4896, 27.8703, 30.3015, 32.7861,
736 35.2622, 37.8173, 40.3819, 42.974, 45.5732, 48.2303, 50.8659, 53.5741, 56.28, 59.0242,
737 61.738, 64.5642, 67.377, 70.1605, 73.0155, 75.8849, 78.7611, 81.7366, 84.7156, 87.7527,
738 90.7217, 93.8402, 96.8516, 100.036, 103.168, 106.483, 109.772, 113.098, 116.491, 120.013,
739 123.618, 127.069, 130.983, 134.868, 138.605, 142.625, 147.007, 151.154, 155.625, 160.1,
740 164.776, 169.651, 174.646, 179.669, 185.084, 190.409, 196.147, 201.788, 207.901, 214.041,
741 220.327, 226.505, 233.334, 239.816, 246.878, 253.563, 260.393, 267.453, 274.5, 282.15,
742 289.014, 296.45, 303.808, 311.427, 318.649, 326.965, 334.298, 341.576, 349.715, 356.89,
743 365.029, 372.677, 379.882, 387.677, 395.178, 402.445, 410.353, 418.649, 424.994, 432.156,
744 440.002, 448.394, 454.382, 460.97, 468.446, 475.847, 481.956, 489.729, 496.094, 501.22,
745 509.278, 514.618, 521.06, 528.247, 534.246, 540.312, 547.316, 552.549, 559.193, 566.059,
746 572.882, 578.147, 585.118, 589.989, 596.717, 601.222, 607.749, 613.96, 621.107, 625.218,
747 630.396, 635.57, 641.175, 646.024, 651.984, 657.156, 661.385, 666.804, 672.088, 675.939,
748 681.207, 685.072, 690.63, 694.767, 699.469, 704.1, 709.445, 713.704, 716.909, 720.681,
749 726.12, 730.403, 733.553, 739.123, 742.156, 746.6, 750.027, 753.462, 757.426, 761.595,
750 764.336, 768.251, 772.371, 775.963, 778.886, 781.905, 784.798, 788.825, 792.372, 796.27,
751 800.361, 803.544, 806.544, 808.819, 812.146, 814.989, 819.234, 820.073, 824.067, 828.047,
752 830.277, 833.013, 835.374, 838.463, 840.82, 844.655, 846.391, 849.408, 851.659, 853.977,
753 856.409, 860.029, 862.128, 866.104, 866.864, 869.24, 872.133, 872.591, 876.528, 879.029,
754 880.786, 883.8, 886.065, 887.511, 890.301, 892.086, 894.429, 895.666, 897.961, 900.712,
755 901.559, 904.787, 906.882, 908.034, 911.366, 911.249, 914.274, 916.238, 918.105, 920.585,
756 920.473, 924.468, 923.888, 926.046, 928.648, 930.3, 931.861, 934.253, 934.081, 936.95,
757 938.319, 940.464, 940.539, 943.393, 944.729, 946.944, 947.712, 948.948, 951.026, 952.121,
758 954.114, 955.146, 956.206, 959.056, 960.316, 962.919, 961.946, 964.324, 966.134, 967.689,
759 968.612, 970.357, 972.302, 973.514, 976.512, 975.815, 979.043, 979.486, 981.285, 983.173,
760 983.96, 985.947, 987.447, 988.455, 991.739, 992.1, 993.045, 995.918, 997.377, 999.136,
761 1001.51, 1001.12, 1002.46, 1004.57, 1005.76, 1007.12, 1009.23, 1011.7, 1012.48, 1014.84,
762 1014.21, 1017.28, 1017.22, 1018.95, 1021.8, 1021.94, 1023.22, 1025.13, 1026.01, 1027.8,
763 1030.04, 1030.12, 1031.54, 1033.2, 1034.62, 1035.83, 1037.33, 1037.92, 1038.9, 1041.69 };
764 return wid1[i];
765}
766
767double D0Topippim2pi0::anywid1300(double sc){
768
769 double smin = (0.13957*3)*(0.13957*3);
770 double dh = 0.01;
771 int od = (sc - 0.18)/dh;
772 double sc_m = 0.18 + od*dh;
773 double widuse = 0;
774 if(sc>=0.18 && sc<=3.17){
775 widuse = ((sc-sc_m)/dh)*(widT1300(od+1)-widT1300(od))+widT1300(od);
776 }else if(sc<0.18 && sc>smin){
777 widuse = ((sc - smin)/(0.18-smin))*widT1300(0);
778 }else if(sc>3.17){
779 widuse = widT1300(299);
780 }else{
781 widuse = 0;
782 }
783 return widuse;
784}
785
786complex<double> D0Topippim2pi0::RBWpi1300(double mx2, double mr, double wr){
787
788 double mx = sqrt(mx2);
789 double mr2 = mr*mr;
790 double g1 = wr/anywid1300(mr2);
791 double wid0 = anywid1300(mx2)*g1;
792
793 double denom_real = mr2-mx2;
794 double denom_imag = mr*wid0;//real-i*imag;
795
796 double denom = denom_real*denom_real+denom_imag*denom_imag;
797 double output_x = denom_real/denom;
798 double output_y = denom_imag/denom;
799
800 complex<double> output(output_x,output_y);
801 return output;
802
803}
804
805/* build a1640 propagator */
806double D0Topippim2pi0::widT1640(int i){
807 double wid1[300] = { 1.38316e-05, 0.000403892, 0.00181814, 0.0048161, 0.00982907, 0.0172548, 0.0273979, 0.040567, 0.0569061, 0.0768551,
808 0.100513, 0.128031, 0.159729, 0.195626, 0.236099, 0.280881, 0.330745, 0.386095, 0.446448, 0.511879,
809 0.583827, 0.66167, 0.745453, 0.835386, 0.934317, 1.0386, 1.1513, 1.26975, 1.39901, 1.53362,
810 1.68291, 1.84163, 2.0066, 2.18366, 2.37394, 2.57742, 2.7905, 3.02463, 3.27434, 3.53467,
811 3.80737, 4.10838, 4.41975, 4.76341, 5.12572, 5.51301, 5.91839, 6.36597, 6.8457, 7.33806,
812 7.87328, 8.45901, 9.08869, 9.74744, 10.464, 11.2096, 12.0103, 12.8556, 13.7563, 14.7352,
813 15.7336, 16.7432, 17.8117, 18.9327, 20.0186, 21.1632, 22.3549, 23.5172, 24.6518, 25.7808,
814 26.9103, 28.016, 29.1542, 30.0458, 31.0808, 32.1018, 33.0395, 33.9151, 34.8873, 35.7289,
815 36.5603, 37.2489, 38.023, 38.7983, 39.55, 40.2977, 40.8819, 41.4564, 42.1864, 42.7368,
816 43.3923, 43.8651, 44.4667, 44.8108, 45.3935, 45.9551, 46.2652, 46.8683, 47.1943, 47.6864,
817 48.1666, 48.5599, 48.8894, 49.1867, 49.6234, 49.9326, 50.4594, 50.6707, 51.005, 51.2612,
818 51.7638, 51.8946, 52.3176, 52.5107, 52.7378, 52.9418, 53.4019, 53.3571, 53.7937, 54.137,
819 54.2265, 54.3471, 54.6637, 54.897, 55.2174, 55.1577, 55.7098, 55.8616, 55.8862, 56.2106,
820 56.3357, 56.5165, 56.6819, 56.7906, 56.9814, 57.0507, 57.3059, 57.4898, 57.5848, 57.5792,
821 57.7696, 58.0302, 58.1915, 58.3319, 58.3892, 58.4671, 58.6736, 58.7872, 58.7949, 58.8366,
822 59.0247, 59.0881, 59.2675, 59.479, 59.6261, 59.6111, 59.6055, 59.7286, 59.8806, 60.0424,
823 60.1126, 60.0742, 60.2066, 60.2253, 60.565, 60.6557, 60.7359, 60.6405, 60.6429, 60.8521,
824 60.8098, 61.0699, 61.1678, 61.0329, 61.0522, 61.1792, 61.3671, 61.4394, 61.5152, 61.6122,
825 61.584, 61.711, 61.707, 61.7254, 61.816, 61.9248, 61.9748, 61.9498, 62.0014, 62.0634,
826 62.2929, 62.2349, 62.2101, 62.4434, 62.4281, 62.4166, 62.4905, 62.6055, 62.5097, 62.5994,
827 62.6637, 62.6794, 62.7068, 62.7908, 62.8135, 63.0085, 62.8848, 62.8159, 63.047, 62.8632,
828 63.1119, 63.0864, 63.1423, 63.2334, 63.0695, 63.2902, 63.3719, 63.1882, 63.2649, 63.3338,
829 63.4709, 63.4662, 63.3746, 63.623, 63.6402, 63.5632, 63.6611, 63.6012, 63.5904, 63.7467,
830 63.5535, 63.7792, 63.5213, 63.829, 63.8696, 63.8047, 63.9557, 63.9433, 63.9363, 63.9436,
831 63.9804, 64.0707, 64.0105, 63.96, 64.0437, 64.0235, 64.1795, 64.1377, 64.073, 64.2282,
832 64.2933, 64.4369, 64.3887, 64.2474, 64.2373, 64.3553, 64.425, 64.4401, 64.3197, 64.4212,
833 64.5787, 64.4919, 64.6878, 64.4998, 64.5788, 64.6628, 64.6658, 64.5072, 64.7227, 64.7327,
834 64.4472, 64.6792, 64.7801, 64.5715, 64.7263, 64.8505, 64.7488, 64.6448, 64.8962, 64.8815,
835 64.821, 64.902, 64.8944, 64.8959, 64.8957, 64.7882, 65.0725, 64.8787, 64.797, 65.1112,
836 65.1212, 65.157, 64.9412, 65.2601, 65.0662, 65.0093, 65.0899, 65.1035, 65.0865, 65.3276 };
837 return wid1[i];
838}
839
840double D0Topippim2pi0::anywid1640(double sc){
841
842 double smin = (0.13957*3)*(0.13957*3);
843 double dh = 0.01;
844 int od = (sc - 0.18)/dh;
845 double sc_m = 0.18 + od*dh;
846 double widuse = 0;
847 if(sc>=0.18 && sc<=3.17){
848 widuse = ((sc-sc_m)/dh)*(widT1640(od+1)-widT1640(od))+widT1640(od);
849 }else if(sc<0.18 && sc>smin){
850 widuse = ((sc - smin)/(0.18-smin))*widT1640(0);
851 }else if(sc>3.17){
852 widuse = widT1640(299);
853 }else{
854 widuse = 0;
855 }
856 return widuse;
857}
858
859complex<double> D0Topippim2pi0::RBWa1640(double mx2, double mr, double wr){
860
861 double mx = sqrt(mx2);
862 double mr2 = mr*mr;
863 double g1 = wr/anywid1640(mr2);
864 double wid0 = anywid1640(mx2)*g1;
865
866 double denom_real = mr2-mx2;
867 double denom_imag = mr*wid0;//real-i*imag;
868
869 double denom = denom_real*denom_real+denom_imag*denom_imag;
870 double output_x = denom_real/denom;
871 double output_y = denom_imag/denom;
872
873 complex<double> output(output_x,output_y);
874 return output;
875
876}
877
878// PiPi-S wave K-Matrix
879double D0Topippim2pi0::rho22(double sc){
880 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,
881 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,
882 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,
883 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,
884 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,
885 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,
886 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,
887 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,
888 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,
889 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,
890 0.000103559, 0.000108812, 0.000114267, 0.000119938, 0.000125827, 0.00013194, 0.000138278, 0.000144857, 0.000151681, 0.000158752,
891 0.000166074, 0.000173663, 0.000181521, 0.000189652, 0.000198059, 0.000206761, 0.000215761, 0.000225063, 0.00023467, 0.000244599,
892 0.000254855, 0.00026544, 0.000276357, 0.000287629, 0.00029926, 0.000311253, 0.000323609, 0.000336351, 0.000349483, 0.000363009,
893 0.000376926, 0.000391264, 0.000406029, 0.000421225, 0.000436848, 0.000452921, 0.000469458, 0.000486461, 0.00050393, 0.00052187,
894 0.000540322, 0.000559278, 0.000578746, 0.00059872, 0.000619236, 0.0006403, 0.000661911, 0.000684074, 0.000706799, 0.000730127,
895 0.00075405, 0.000778569, 0.000803686, 0.000829443, 0.000855839, 0.000882879, 0.000910561, 0.000938898, 0.000967939, 0.000997674,
896 0.00102811, 0.00105923, 0.0010911, 0.0011237, 0.00115706, 0.00119117, 0.00122601, 0.00126168, 0.00129815, 0.00133543,
897 0.00137351, 0.00141242, 0.00145219, 0.00149283, 0.00153434, 0.0015767, 0.00161995, 0.00166415, 0.00170928, 0.00175534,
898 0.00180232, 0.00185028, 0.00189924, 0.00194919, 0.00200014, 0.00205207, 0.00210503, 0.0021591, 0.00221421, 0.0022704,
899 0.00232766, 0.00238602, 0.00244554, 0.00250619, 0.00256799, 0.0026309, 0.002695, 0.00276033, 0.00282689, 0.00289467,
900 0.00296367, 0.00303389, 0.00310543, 0.0031783, 0.00325244, 0.0033279, 0.0034046, 0.00348275, 0.00356229, 0.00364322,
901 0.00372555, 0.00380924, 0.00389438, 0.00398104, 0.00406914, 0.00415877, 0.00424985, 0.00434235, 0.00443651, 0.00453224,
902 0.00462954, 0.00472848, 0.00482894, 0.00493102, 0.00503483, 0.00514029, 0.00524749, 0.0053563, 0.00546675, 0.00557905,
903 0.0056931, 0.00580901, 0.0059267, 0.00604613, 0.00616735, 0.00629049, 0.00641557, 0.00654254, 0.00667142, 0.00680216,
904 0.00693472, 0.00706946, 0.00720621, 0.00734497, 0.0074858, 0.00762855, 0.00777338, 0.00792036, 0.00806957, 0.00822087,
905 0.00837426, 0.00852982, 0.0086875, 0.00884756, 0.00900991, 0.00917447, 0.00934137, 0.00951052, 0.00968194, 0.0098558,
906 0.010032, 0.0102108, 0.0103919, 0.0105754, 0.0107612, 0.0109496, 0.0111406, 0.0113343, 0.0115305, 0.0117293,
907 0.0119303, 0.0121343, 0.0123409, 0.0125502, 0.0127623, 0.0129771, 0.0131944, 0.0134145, 0.0136376, 0.0138636,
908 0.0140924, 0.0143241, 0.0145587, 0.0147959, 0.0150363, 0.0152797, 0.0155262, 0.0157758, 0.0160283, 0.0162838,
909 0.0165421, 0.016804, 0.0170691, 0.0173374, 0.0176087, 0.0178835, 0.0181612, 0.0184423, 0.0187269, 0.0190149,
910 0.0193063, 0.0196009, 0.0198991, 0.0202003, 0.0205052, 0.0208137, 0.0211259, 0.0214418, 0.0217611, 0.0220841,
911 0.0224105, 0.0227406, 0.0230746, 0.0234125, 0.0237542, 0.0240996, 0.0244486, 0.0248012, 0.025158, 0.0255188,
912 0.0258837, 0.0262527, 0.0266256, 0.0270025, 0.0273833, 0.027768, 0.0281572, 0.0285505, 0.0289483, 0.0293503,
913 0.0297564, 0.0301665, 0.0305808, 0.0309997, 0.0314231, 0.0318511, 0.0322835, 0.0327205, 0.0331616, 0.0336073,
914 0.0340576, 0.0345128, 0.0349727, 0.0354373, 0.0359066, 0.0363807, 0.0368589, 0.0373419, 0.0378302, 0.0383234,
915 0.0388218, 0.0393252, 0.0398336, 0.040347, 0.0408652, 0.041388, 0.0419165, 0.0424502, 0.0429893, 0.0435338,
916 0.0440833, 0.044638, 0.0451976, 0.0457627, 0.0463338, 0.0469103, 0.047492, 0.0480797, 0.0486729, 0.0492716,
917 0.0498757, 0.0504852, 0.0511009, 0.0517229, 0.0523503, 0.0529838, 0.0536231, 0.0542678, 0.054918, 0.0555743,
918 0.0562372, 0.0569065, 0.0575818, 0.0582634, 0.0589511, 0.0596454, 0.0603451, 0.061051, 0.0617635, 0.0624826,
919 0.0632084, 0.0639409, 0.06468, 0.0654254, 0.0661772, 0.0669346, 0.0676994, 0.0684714, 0.0692503, 0.0700354,
920 0.0708285, 0.0716277, 0.0724347, 0.0732479, 0.0740671, 0.0748947, 0.0757299, 0.0765715, 0.0774207, 0.0782771,
921 0.0791407, 0.0800119, 0.0808897, 0.0817743, 0.0826672, 0.0835684, 0.0844769, 0.0853938, 0.0863179, 0.0872493,
922 0.0881882, 0.0891349, 0.090089, 0.0910523, 0.0920236, 0.093002, 0.0939894, 0.094985, 0.0959887, 0.0970003,
923 0.0980191, 0.0990454, 0.100081, 0.101126, 0.10218, 0.103242, 0.104312, 0.105392, 0.10648, 0.107576,
924 0.10868, 0.109793, 0.110916, 0.112048, 0.113188, 0.114339, 0.115498, 0.116666, 0.117843, 0.119028,
925 0.120223, 0.121427, 0.122641, 0.123865, 0.125098, 0.126342, 0.127595, 0.128857, 0.130128, 0.131409,
926 0.132701, 0.134002, 0.135314, 0.136635, 0.137966, 0.139308, 0.14066, 0.142022, 0.143394, 0.144774,
927 0.146166, 0.14757, 0.148985, 0.15041, 0.151845, 0.153291, 0.154749, 0.156215, 0.157694, 0.159182,
928 0.160682, 0.162194, 0.163718, 0.165251, 0.166797, 0.168354, 0.169921, 0.1715, 0.17309, 0.17469,
929 0.176304, 0.177929, 0.179566, 0.181216, 0.182878, 0.184553, 0.186238, 0.187934, 0.189642, 0.191362,
930 0.193096, 0.194842, 0.196602, 0.198374, 0.200158, 0.201954, 0.203764, 0.205586, 0.207421, 0.209266,
931 0.211124, 0.212997, 0.214882, 0.216783, 0.218697, 0.220624, 0.222565, 0.224518, 0.226486, 0.228466,
932 0.230458, 0.232463, 0.234484, 0.23652, 0.238569, 0.240633, 0.242711, 0.244803, 0.246909, 0.249031,
933 0.251165, 0.253313, 0.255475, 0.257649, 0.259841, 0.262051, 0.264274, 0.266514, 0.268768, 0.271036,
934 0.273319, 0.275618, 0.277932, 0.280259, 0.282602, 0.28496, 0.287338, 0.28973, 0.292138, 0.294563,
935 0.297003, 0.299458, 0.30193, 0.304417, 0.306919, 0.309437, 0.311972, 0.314526, 0.317095, 0.319684,
936 0.322289, 0.324911, 0.327551, 0.330205, 0.332876, 0.335567, 0.338271, 0.340993, 0.343736, 0.346496,
937 0.349272, 0.352065, 0.354878, 0.35771, 0.360561, 0.363426, 0.366311, 0.369212, 0.372128, 0.375067,
938 0.378027, 0.381006, 0.384001, 0.387014, 0.39005, 0.393106, 0.396181, 0.399271, 0.402384, 0.405513,
939 0.408661, 0.41183, 0.41502, 0.418233, 0.421462, 0.424709, 0.42798, 0.43127, 0.434583, 0.437914,
940 0.441267, 0.444637, 0.448022, 0.451434, 0.454868, 0.458328, 0.461805, 0.465302, 0.468821, 0.472364,
941 0.475928, 0.47951, 0.483119, 0.486748, 0.490397, 0.494066, 0.497758, 0.501477, 0.505217, 0.508977,
942 0.512762, 0.516567, 0.520394, 0.524247, 0.528125, 0.532027, 0.535947, 0.53989, 0.543852, 0.547844,
943 0.551863, 0.555904, 0.559966, 0.56406, 0.568177, 0.572312, 0.576471, 0.580662, 0.584875, 0.58911,
944 0.593373, 0.597653, 0.601965, 0.606301, 0.610663, 0.615051, 0.619465, 0.623907, 0.62837, 0.632863,
945 0.637383, 0.641924, 0.646494, 0.651091, 0.655708, 0.660356, 0.665027, 0.669732, 0.674464, 0.679227,
946 0.684016, 0.688827, 0.693664, 0.698532, 0.703428, 0.708353, 0.713307, 0.718283, 0.72329, 0.728322,
947 0.733387, 0.738479, 0.743605, 0.748763, 0.753949, 0.759163, 0.764407, 0.769674, 0.774973, 0.780311,
948 0.78567, 0.791057, 0.796476, 0.801922, 0.8074, 0.812919, 0.818466, 0.824044 };
949
950 double m2 = 0.13957*0.13957;
951 double smin = (0.13957*4)*(0.13957*4);
952 double dh = 0.001;
953 int od = (sc - 0.312)/dh;
954 double sc_m = 0.312 + od*dh;
955 double rhouse = 0;
956 if(sc>=0.312 && sc<1){
957 rhouse = ((sc-sc_m)/dh)*(rho[od+1]-rho[od])+rho[od];
958 }else if(sc<0.312 && sc>=smin){
959 rhouse = ((sc - smin)/(0.312-smin))*rho[0];
960 }else if(sc>=1){
961 // rhouse = (1-16*m2/sc)*(1-16*m2/sc)*(1-16*m2/sc)*(1-16*m2/sc)*(1-16*m2/sc);
962 rhouse = sqrt(1-16*m2/sc);
963 }else{
964 rhouse = 0;
965 }
966 return rhouse;
967}
968
969complex<double> D0Topippim2pi0::rhoMTX(int i, int j, double s){
970
971 double rhoijx;
972 double rhoijy;
973 double mpi = 0.13957;
974 if(i==j && i==0 ){
975 double m2 = 0.13957*0.13957;
976 if((1-(4*m2)/s)>0){
977 rhoijx = sqrt(1.0f - (4*m2)/s);
978 rhoijy = 0;
979 }else{
980 rhoijy = sqrt((4*m2)/s - 1.0f);
981 rhoijx = 0;
982 }
983 }
984 if(i==j && i==1 ){
985 double m2 = 0.493677*0.493677;
986 if((1-(4*m2)/s)>0){
987 rhoijx = sqrt(1.0f - (4*m2)/s);
988 rhoijy = 0;
989 }else{
990 rhoijy = sqrt((4*m2)/s - 1.0f);
991 rhoijx = 0;
992 }
993 }
994 if(i==j && i==2){
995 rhoijx = rho22(s);
996 rhoijy = 0;
997 }
998 if(i==j && i==3){
999 double m2 = 0.547862*0.547862;
1000 if((1-(4*m2)/s)>0){
1001 rhoijx = sqrt(1.0f - (4*m2)/s);
1002 rhoijy = 0;
1003 }else{
1004 rhoijy = sqrt((4*m2)/s - 1.0f);
1005 rhoijx = 0;
1006 }
1007 }
1008 if(i==j && i==4){
1009 double m_1 = 0.547862;
1010 double m_2 = 0.95778;
1011 double mp2 = (m_1+m_2)*(m_1+m_2);
1012 double mm2 = (m_1-m_2)*(m_1-m_2);
1013 if((1-mp2/s)>0){
1014 rhoijx = sqrt(1.0f - mp2/s);
1015 rhoijy = 0;
1016 }else{
1017 rhoijy = sqrt(mp2/s - 1.0f);
1018 rhoijx = 0;
1019 }
1020 }
1021
1022 if(i!=j){
1023 rhoijx = 0;
1024 rhoijy = 0;
1025 }
1026 complex<double> rhoij(rhoijx,rhoijy);
1027 return rhoij;
1028
1029}
1030/* Kij = (sum_a gigj/(ma^2-s) + fij(1-s0)/(s-s0))((s-0.5sampi2)(1-sa0)/(s-sa0)) */
1031complex<double> D0Topippim2pi0::KMTX(int i, int j, double s){
1032
1033 double Kijx;
1034 double Kijy;
1035 double mpi = 0.13957;
1036 double m[5] = { 0.65100, 1.20360, 1.55817, 1.21000, 1.82206};
1037
1038 double g1[5] = { 0.22889,-0.55377, 0.00000,-0.39899,-0.34639};
1039 double g2[5] = { 0.94128, 0.55095, 0.00000, 0.39065, 0.31503};
1040 double g3[5] = { 0.36856, 0.23888, 0.55639, 0.18340, 0.18681};
1041 double g4[5] = { 0.33650, 0.40907, 0.85679, 0.19906,-0.00984};
1042 double g5[5] = { 0.18171,-0.17558,-0.79658,-0.00355, 0.22358};
1043
1044 double f1[5] = { 0.23399, 0.15044,-0.20545, 0.32825, 0.35412};
1045
1046 double eps = 1e-11;
1047
1048 double down[5] = { 0,0,0,0,0};
1049 double upreal[5] = { 0,0,0,0,0};
1050 double upimag[5] = { 0,0,0,0,0};
1051
1052 for(int k=0; k<5; k++){
1053
1054 /* down[k] = (m[k]*m[k]-s)*(m[k]*m[k]-s)+eps*eps;
1055 upreal[k] = (m[k]*m[k]-s)/down[k];
1056 upimag[k] = -1.0f * eps/down[k]; */
1057
1058 double dm2 = m[k]*m[k]-s;
1059 if(fabs(dm2)<eps && dm2<=0) dm2 = -eps;
1060 if(fabs(dm2)<eps && dm2>0) dm2 = eps;
1061 upreal[k] = 1.0f/dm2;
1062 upimag[k] = 0;
1063 }
1064
1065 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]*g5[j]*upreal[4];
1066 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]*g5[j]*upimag[4];
1067
1068 double tmp2 = 0;
1069 if(i==0){
1070 tmp2 = f1[j]*(1+3.92637)/(s+3.92637);
1071 }
1072 if(j==0){
1073 tmp2 = f1[i]*(1+3.92637)/(s+3.92637);
1074 }
1075 double tmp3 = (s-0.5*mpi*mpi)*(1+0.15)/(s+0.15);
1076
1077 Kijx = (tmp1x+tmp2)*tmp3;
1078 Kijy = (tmp1y)*tmp3;
1079
1080 complex<double> Kij(Kijx,Kijy);
1081 return Kij;
1082}
1083
1084complex<double> D0Topippim2pi0::IMTX(int i, int j){
1085
1086 double Iijx;
1087 double Iijy;
1088 if(i==j){
1089 Iijx = 1;
1090 Iijy = 0;
1091 }else{
1092 Iijx = 0;
1093 Iijy = 0;
1094 }
1095 complex<double> Iij(Iijx,Iijy);
1096 return Iij;
1097
1098}
1099
1100/* F = I - i*K*rho */
1101complex<double> D0Topippim2pi0::FMTX(double Kijx, double Kijy, double rhojjx, double rhojjy, int i, int j){
1102
1103 double Fijx;
1104 double Fijy;
1105
1106 double tmpx = rhojjx*Kijx - rhojjy*Kijy;
1107 double tmpy = rhojjx*Kijy + rhojjy*Kijx;
1108
1109 Fijx = IMTX(i,j).real() + tmpy;
1110 Fijy = -tmpx;
1111
1112 complex<double> Fij(Fijx,Fijy);
1113 return Fij;
1114}
1115
1116/* inverse for Matrix (I - i*rho*K) using LUP */
1117double D0Topippim2pi0::FINVMTX(double s, double *FINVx, double *FINVy){
1118
1119 int P[5] = { 0,1,2,3,4};
1120
1121 double Fx[5][5];
1122 double Fy[5][5];
1123
1124 double Ux[5][5];
1125 double Uy[5][5];
1126 double Lx[5][5];
1127 double Ly[5][5];
1128
1129 double UIx[5][5];
1130 double UIy[5][5];
1131 double LIx[5][5];
1132 double LIy[5][5];
1133
1134 for(int k=0; k<5; k++){
1135 double rhokkx = rhoMTX(k,k,s).real();
1136 double rhokky = rhoMTX(k,k,s).imag();
1137 Ux[k][k] = rhokkx;
1138 Uy[k][k] = rhokky;
1139 for(int l=k; l<5; l++){
1140 double Kklx = KMTX(k,l,s).real();
1141 double Kkly = KMTX(k,l,s).imag();
1142 Lx[k][l] = Kklx;
1143 Ly[k][l] = Kkly;
1144 Lx[l][k] = Lx[k][l];
1145 Ly[l][k] = Ly[k][l];
1146 }
1147 }
1148
1149 for(int k=0; k<5; k++){
1150 for(int l=0; l<5; l++){
1151 double Fklx = FMTX(Lx[k][l],Ly[k][l],Ux[l][l],Uy[l][l],k,l).real();
1152 double Fkly = FMTX(Lx[k][l],Ly[k][l],Ux[l][l],Uy[l][l],k,l).imag();
1153 Fx[k][l] = Fklx;
1154 Fy[k][l] = Fkly;
1155 }
1156 }
1157
1158 for(int k=0; k<5; k++){
1159 double tmprM = (Fx[k][k]*Fx[k][k]+Fy[k][k]*Fy[k][k]);
1160 int tmpID = 0;
1161 for(int l=k; l<5; l++){
1162 double tmprF = (Fx[l][k]*Fx[l][k]+Fy[l][k]*Fy[l][k]);
1163 if(tmprM<=tmprF){
1164 tmprM = tmprF;
1165 tmpID = l;
1166 }
1167 }
1168
1169 int tmpP = P[k];
1170 P[k] = P[tmpID];
1171 P[tmpID] = tmpP;
1172
1173 for(int l=0; l<5; l++){
1174
1175 double tmpFx = Fx[k][l];
1176 double tmpFy = Fy[k][l];
1177
1178 Fx[k][l] = Fx[tmpID][l];
1179 Fy[k][l] = Fy[tmpID][l];
1180
1181 Fx[tmpID][l] = tmpFx;
1182 Fy[tmpID][l] = tmpFy;
1183
1184 }
1185
1186 for(int l=k+1; l<5; l++){
1187 double rFkk = Fx[k][k]*Fx[k][k] + Fy[k][k]*Fy[k][k];
1188 double Fxlk = Fx[l][k];
1189 double Fylk = Fy[l][k];
1190 double Fxkk = Fx[k][k];
1191 double Fykk = Fy[k][k];
1192 Fx[l][k] = (Fxlk*Fxkk + Fylk*Fykk)/rFkk;
1193 Fy[l][k] = (Fylk*Fxkk - Fxlk*Fykk)/rFkk;
1194 for(int m=k+1; m<5; m++){
1195 Fx[l][m] = Fx[l][m] - (Fx[l][k]*Fx[k][m] - Fy[l][k]*Fy[k][m]);
1196 Fy[l][m] = Fy[l][m] - (Fx[l][k]*Fy[k][m] + Fy[l][k]*Fx[k][m]);
1197 }
1198 }
1199 }
1200
1201 for(int k=0; k<5; k++){
1202 for(int l=0; l<5 ;l++){
1203 if(k==l){
1204 Lx[k][k] = 1;
1205 Ly[k][k] = 0;
1206 Ux[k][k] = Fx[k][k];
1207 Uy[k][k] = Fy[k][k];
1208 }
1209 if(k>l){
1210 Lx[k][l] = Fx[k][l];
1211 Ly[k][l] = Fy[k][l];
1212 Ux[k][l] = 0;
1213 Uy[k][l] = 0;
1214 }
1215 if(k<l){
1216 Ux[k][l] = Fx[k][l];
1217 Uy[k][l] = Fy[k][l];
1218 Lx[k][l] = 0;
1219 Ly[k][l] = 0;
1220 }
1221 }
1222 }
1223
1224 // calculate Inverse for L and U
1225 for(int k=0; k<5; k++){
1226
1227 LIx[k][k] = 1;
1228 LIy[k][k] = 0;
1229
1230 double rUkk = Ux[k][k]*Ux[k][k] + Uy[k][k]*Uy[k][k];
1231 UIx[k][k] = Ux[k][k]/rUkk;
1232 UIy[k][k] = -1.0f * Uy[k][k]/rUkk ;
1233
1234 for(int l=(k+1); l<5; l++){
1235 LIx[k][l] = 0;
1236 LIy[k][l] = 0;
1237 UIx[l][k] = 0;
1238 UIy[l][k] = 0;
1239 }
1240
1241 for(int l=(k-1); l>=0; l--){ // U-1
1242 double sx = 0;
1243 double c_sx = 0;
1244 double sy = 0;
1245 double c_sy = 0;
1246 for(int m=l+1; m<=k; m++){
1247 sx = sx - c_sx;
1248 double sx_tmp = sx + Ux[l][m]*UIx[m][k] - Uy[l][m]*UIy[m][k];
1249 c_sx = (sx_tmp - sx) - (Ux[l][m]*UIx[m][k] - Uy[l][m]*UIy[m][k]);
1250 sx = sx_tmp;
1251
1252 sy = sy - c_sy;
1253 double sy_tmp = sy + Ux[l][m]*UIy[m][k] + Uy[l][m]*UIx[m][k];
1254 c_sy = (sy_tmp - sy) - (Ux[l][m]*UIy[m][k] + Uy[l][m]*UIx[m][k]);
1255 sy = sy_tmp;
1256 }
1257 UIx[l][k] = -1.0f * (UIx[l][l]*sx - UIy[l][l]*sy);
1258 UIy[l][k] = -1.0f * (UIy[l][l]*sx + UIx[l][l]*sy);
1259 }
1260
1261 for(int l=k+1; l<5; l++){ // L-1
1262 double sx = 0;
1263 double c_sx = 0;
1264 double sy = 0;
1265 double c_sy = 0;
1266 for(int m=k; m<l; m++){
1267 sx = sx - c_sx;
1268 double sx_tmp = sx + Lx[l][m]*LIx[m][k] - Ly[l][m]*LIy[m][k];
1269 c_sx = (sx_tmp - sx) - (Lx[l][m]*LIx[m][k] - Ly[l][m]*LIy[m][k]);
1270 sx = sx_tmp;
1271
1272 sy = sy - c_sy;
1273 double sy_tmp = sy + Lx[l][m]*LIy[m][k] + Ly[l][m]*LIx[m][k];
1274 c_sy = (sy_tmp - sy) - (Lx[l][m]*LIy[m][k] + Ly[l][m]*LIx[m][k]);
1275 sy = sy_tmp;
1276 }
1277 LIx[l][k] = -1.0f * sx;
1278 LIy[l][k] = -1.0f * sy;
1279 }
1280 }
1281
1282 for(int m=0; m<5; m++){
1283 double resX = 0;
1284 double c_resX = 0;
1285 double resY = 0;
1286 double c_resY = 0;
1287 for(int k=0; k<5; k++){
1288 for(int l=0; l<5; l++){
1289 double Plm = 0;
1290 if(P[l] == m) Plm = 1;
1291
1292 resX = resX - c_resX;
1293 double resX_tmp = resX + (UIx[0][k]*LIx[k][l] - UIy[0][k]*LIy[k][l])*Plm;
1294 c_resX = (resX_tmp - resX) - ((UIx[0][k]*LIx[k][l] - UIy[0][k]*LIy[k][l])*Plm);
1295 resX = resX_tmp;
1296
1297 resY = resY - c_resY;
1298 double resY_tmp = resY + (UIx[0][k]*LIy[k][l] + UIy[0][k]*LIx[k][l])*Plm;
1299 c_resY = (resY_tmp - resY) - ((UIx[0][k]*LIy[k][l] + UIy[0][k]*LIx[k][l])*Plm);
1300 resY = resY_tmp;
1301 }
1302 }
1303 FINVx[m] = resX;
1304 FINVy[m] = resY;
1305 }
1306
1307 return 1.0f;
1308}
1309
1310complex<double> D0Topippim2pi0::PVTR(int ID, double s){
1311
1312 double VPix;
1313 double VPiy;
1314 double m[5] = { 0.65100, 1.20360, 1.55817, 1.21000, 1.82206};
1315
1316 double eps = 1e-11;
1317
1318 /* double down = (m[ID]*m[ID]-s)*(m[ID]*m[ID]-s)+eps*eps;
1319 double upreal = (m[ID]*m[ID]-s)/down;
1320 double upimag = -1.0f * eps/down; */
1321
1322 double dm2 = m[ID]*m[ID]-s;
1323
1324 if(fabs(dm2)<eps && dm2<=0) dm2 = -eps;
1325 if(fabs(dm2)<eps && dm2>0) dm2 = eps;
1326
1327 VPix = 1.0f/dm2;
1328 VPiy = 0;
1329
1330 complex<double> VPi(VPix,VPiy);
1331 return VPi;
1332}
1333
1334complex<double> D0Topippim2pi0::Fvector(double sa, double s0, int l){
1335
1336 double outputx = 0;
1337 double outputy = 0;
1338
1339 double FINVx[5] = {0,0,0,0,0};
1340 double FINVy[5] = {0,0,0,0,0};
1341
1342 double tmpFLAG = FINVMTX(sa,FINVx,FINVy);
1343
1344 if(l<5){
1345 double g[5][5] = {{ 0.22889,-0.55377, 0.00000,-0.39899,-0.34639},
1346 { 0.94128, 0.55095, 0.00000, 0.39065, 0.31503},
1347 { 0.36856, 0.23888, 0.55639, 0.18340, 0.18681},
1348 { 0.33650, 0.40907, 0.85679, 0.19906,-0.00984},
1349 { 0.18171,-0.17558,-0.79658,-0.00355, 0.22358}};
1350 double resx = 0;
1351 double c_resx = 0;
1352 double resy = 0;
1353 double c_resy = 0;
1354 double Plx = PVTR(l,sa).real();
1355 double Ply = PVTR(l,sa).imag();
1356 for(int j=0; j<5; j++){
1357 resx = resx - c_resx;
1358 double resx_tmp = resx + (FINVx[j]*g[l][j]*Plx - FINVy[j]*g[l][j]*Ply);
1359 c_resx = (resx_tmp - resx) - (FINVx[j]*g[l][j]*Plx - FINVy[j]*g[l][j]*Ply);
1360 resx = resx_tmp;
1361
1362 resy = resy - c_resy;
1363 double resy_tmp = resy + (FINVx[j]*g[l][j]*Ply + FINVy[j]*g[l][j]*Plx);
1364 c_resy = (resy_tmp - resy) - (FINVx[j]*g[l][j]*Ply + FINVy[j]*g[l][j]*Plx);
1365 resy = resy_tmp;
1366 }
1367 outputx = resx;
1368 outputy = resy;
1369 }else{
1370 int idx = l-5;
1371 double eps = 1e-11;
1372 double ds = sa-s0;
1373 if(fabs(ds)<eps && ds<=0) ds = -eps;
1374 if(fabs(ds)<eps && ds>0) ds = eps;
1375 double tmp = (1-s0)/ds;
1376 outputx = FINVx[idx]*tmp;
1377 outputy = FINVy[idx]*tmp;
1378 }
1379
1380 complex<double> output(outputx,outputy);
1381 return output;
1382}
1383
1384complex<double> D0Topippim2pi0::Amp(vector<double> Pip, vector<double> Pim, vector<double> Pi01, vector<double> Pi02)
1385{
1386
1387 vector<double> PipPim; PipPim.clear();
1388 vector<double> PipPi01; PipPi01.clear();
1389 vector<double> PipPi02; PipPi02.clear();
1390 vector<double> PimPi01; PimPi01.clear();
1391 vector<double> PimPi02; PimPi02.clear();
1392 vector<double> Pi01Pi02; Pi01Pi02.clear();
1393
1394 PipPim = sum_tensor(Pip, Pim);
1395 PipPi01 = sum_tensor(Pip, Pi01);
1396 PipPi02 = sum_tensor(Pip, Pi02);
1397 PimPi01 = sum_tensor(Pim, Pi01);
1398 PimPi02 = sum_tensor(Pim, Pi02);
1399 Pi01Pi02 = sum_tensor(Pi01, Pi02);
1400
1401 vector<double> PipPimPi01; PipPimPi01.clear();
1402 vector<double> PipPimPi02; PipPimPi02.clear();
1403 vector<double> PipPi01Pi02; PipPi01Pi02.clear();
1404 vector<double> PimPi01Pi02; PimPi01Pi02.clear();
1405
1406 PipPimPi01 = sum_tensor(PipPim, Pi01);
1407 PipPimPi02 = sum_tensor(PipPim, Pi02);
1408 PipPi01Pi02 = sum_tensor(PipPi01, Pi02);
1409 PimPi01Pi02 = sum_tensor(PimPi01, Pi02);
1410
1411 vector<double> D0; D0.clear();
1412 D0 = sum_tensor(PipPimPi01, Pi02);
1413
1414 double M2_PipPim = contract_11_0(PipPim, PipPim);
1415 double M2_PipPi01 = contract_11_0(PipPi01, PipPi01);
1416 double M2_PipPi02 = contract_11_0(PipPi02, PipPi02);
1417 double M2_PimPi01 = contract_11_0(PimPi01, PimPi01);
1418 double M2_PimPi02 = contract_11_0(PimPi02, PimPi02);
1419 double M2_Pi01Pi02 = contract_11_0(Pi01Pi02, Pi01Pi02);
1420
1421 double M2_PipPimPi01 = contract_11_0(PipPimPi01, PipPimPi01);
1422 double M2_PipPimPi02 = contract_11_0(PipPimPi02, PipPimPi02);
1423 double M2_PipPi01Pi02 = contract_11_0(PipPi01Pi02, PipPi01Pi02);
1424 double M2_PimPi01Pi02 = contract_11_0(PimPi01Pi02, PimPi01Pi02);
1425 double M2_D0 = contract_11_0(D0, D0);
1426
1427 complex<double> GS_rho770_pm = GS(M2_PipPim, m0_rho7700, w0_rho7700, m2_Pi, m2_Pi, rRes, 1);
1428 complex<double> GS_rho770_p1 = GS(M2_PipPi01, m0_rho770p, w0_rho770p, m2_Pi, m2_Pi0, rRes, 1);
1429 complex<double> GS_rho770_p2 = GS(M2_PipPi02, m0_rho770p, w0_rho770p, m2_Pi, m2_Pi0, rRes, 1);
1430 complex<double> GS_rho770_m1 = GS(M2_PimPi01, m0_rho770p, w0_rho770p, m2_Pi, m2_Pi0, rRes, 1);
1431 complex<double> GS_rho770_m2 = GS(M2_PimPi02, m0_rho770p, w0_rho770p, m2_Pi, m2_Pi0, rRes, 1);
1432
1433 complex<double> RBW_f21270_pm = RBW(M2_PipPim, m0_f21270, w0_f21270, m2_Pi, m2_Pi, rRes, 2);
1434 complex<double> RBW_f21270_00 = RBW(M2_Pi01Pi02, m0_f21270, w0_f21270, m2_Pi0, m2_Pi0, rRes, 2);
1435
1436 complex<double> PiPiS_pm_0 = Fvector(M2_PipPim, s0_prod, 0);
1437 complex<double> PiPiS_00_0 = Fvector(M2_Pi01Pi02, s0_prod, 0);
1438
1439 complex<double> PiPiS_pm_1 = Fvector(M2_PipPim, s0_prod, 1);
1440 complex<double> PiPiS_00_1 = Fvector(M2_Pi01Pi02, s0_prod, 1);
1441
1442 complex<double> PiPiS_pm_5 = Fvector(M2_PipPim, s0_prod, 5);
1443 complex<double> PiPiS_00_5 = Fvector(M2_Pi01Pi02, s0_prod, 5);
1444
1445 complex<double> PiPiS_pm_6 = Fvector(M2_PipPim, s0_prod, 6);
1446 complex<double> PiPiS_00_6 = Fvector(M2_Pi01Pi02, s0_prod, 6);
1447
1448 complex<double> RBW_a11260_p = RBWa1260(M2_PipPi01Pi02, m0_a11260, g1_a11260, g2_a11260);
1449 complex<double> RBW_a11260_m = RBWa1260(M2_PimPi01Pi02, m0_a11260, g1_a11260, g2_a11260);
1450 complex<double> RBW_a11260_01 = RBWa1260(M2_PipPimPi01, m0_a11260, g1_a11260, g2_a11260);
1451 complex<double> RBW_a11260_02 = RBWa1260(M2_PipPimPi02, m0_a11260, g1_a11260, g2_a11260);
1452
1453 complex<double> RBW_a11420_p = RBW(M2_PipPi01Pi02, m0_a11420, w0_a11420,-1,-1,-1,-1);
1454 complex<double> RBW_a11420_m = RBW(M2_PimPi01Pi02, m0_a11420, w0_a11420,-1,-1,-1,-1);
1455 complex<double> RBW_a11420_01 = RBW(M2_PipPimPi01, m0_a11420, w0_a11420,-1,-1,-1,-1);
1456 complex<double> RBW_a11420_02 = RBW(M2_PipPimPi02, m0_a11420, w0_a11420,-1,-1,-1,-1);
1457
1458 complex<double> RBW_omega_01 = RBW(M2_PipPimPi01, m0_omega, w0_omega,-1,-1,-1,-1);
1459 complex<double> RBW_omega_02 = RBW(M2_PipPimPi02, m0_omega, w0_omega,-1,-1,-1,-1);
1460
1461 complex<double> RBW_phi_01 = RBW(M2_PipPimPi01, m0_phi, w0_phi,-1,-1,-1,-1);
1462 complex<double> RBW_phi_02 = RBW(M2_PipPimPi02, m0_phi, w0_phi,-1,-1,-1,-1);
1463
1464 complex<double> RBW_a21320_p = RBW(M2_PipPi01Pi02, m0_a21320, w0_a21320,-1,-1,-1,-1);
1465 complex<double> RBW_a21320_m = RBW(M2_PimPi01Pi02, m0_a21320, w0_a21320,-1,-1,-1,-1);
1466
1467 complex<double> RBW_pi1300_p = RBWpi1300(M2_PipPi01Pi02, m0_pi1300, w0_pi1300);
1468 complex<double> RBW_pi1300_m = RBWpi1300(M2_PimPi01Pi02, m0_pi1300, w0_pi1300);
1469 complex<double> RBW_pi1300_01 = RBWpi1300(M2_PipPimPi01, m0_pi1300, w0_pi1300);
1470 complex<double> RBW_pi1300_02 = RBWpi1300(M2_PipPimPi02, m0_pi1300, w0_pi1300);
1471
1472 complex<double> RBW_h11170_01 = RBW(M2_PipPimPi01, m0_h11170, w0_h11170,-1,-1,-1,-1);
1473 complex<double> RBW_h11170_02 = RBW(M2_PipPimPi02, m0_h11170, w0_h11170,-1,-1,-1,-1);
1474
1475 complex<double> RBW_pi21670_01 = RBW(M2_PipPimPi01, m0_pi21670, w0_pi21670,-1,-1,-1,-1);
1476 complex<double> RBW_pi21670_02 = RBW(M2_PipPimPi02, m0_pi21670, w0_pi21670,-1,-1,-1,-1);
1477
1478 // D->XX Projection
1479 vector<double> Proj1_3p; Proj1_3p.clear();
1480 vector<double> Proj1_3m; Proj1_3m.clear();
1481 vector<double> Proj1_3z1; Proj1_3z1.clear();
1482 vector<double> Proj1_3z2; Proj1_3z2.clear();
1483
1484 Proj1_3p = ProjectionTensors(PipPi01Pi02,1);
1485 Proj1_3m = ProjectionTensors(PimPi01Pi02,1);
1486 Proj1_3z1 = ProjectionTensors(PipPimPi01,1);
1487 Proj1_3z2 = ProjectionTensors(PipPimPi02,1);
1488
1489 vector<double> Proj2_3p; Proj2_3p.clear();
1490 vector<double> Proj2_3m; Proj2_3m.clear();
1491 vector<double> Proj2_3z1; Proj2_3z1.clear();
1492 vector<double> Proj2_3z2; Proj2_3z2.clear();
1493
1494 Proj2_3p = ProjectionTensors(PipPi01Pi02,2);
1495 Proj2_3m = ProjectionTensors(PimPi01Pi02,2);
1496 Proj2_3z1 = ProjectionTensors(PipPimPi01,2);
1497 Proj2_3z2 = ProjectionTensors(PipPimPi02,2);
1498
1499 // X->PP Orbital
1500 vector<double> T1_PipPim; T1_PipPim.clear();
1501 vector<double> T1_PipPi01; T1_PipPi01.clear();
1502 vector<double> T1_PipPi02; T1_PipPi02.clear();
1503 vector<double> T1_PimPi01; T1_PimPi01.clear();
1504 vector<double> T1_PimPi02; T1_PimPi02.clear();
1505 vector<double> T1_Pi01Pi02; T1_Pi01Pi02.clear();
1506
1507 T1_PipPim = OrbitalTensors(PipPim, Pip, Pim, rRes, 1);
1508 T1_PipPi01 = OrbitalTensors(PipPi01, Pip, Pi01, rRes, 1);
1509 T1_PipPi02 = OrbitalTensors(PipPi02, Pip, Pi02, rRes, 1);
1510 T1_PimPi01 = OrbitalTensors(PimPi01, Pim, Pi01, rRes, 1);
1511 T1_PimPi02 = OrbitalTensors(PimPi02, Pim, Pi02, rRes, 1);
1512 T1_Pi01Pi02 = OrbitalTensors(Pi01Pi02, Pi01, Pi02, rRes, 1);
1513
1514 vector<double> T2_PipPim; T2_PipPim.clear();
1515 vector<double> T2_Pi01Pi02; T2_Pi01Pi02.clear();
1516
1517 T2_PipPim = OrbitalTensors(PipPim, Pip, Pim, rRes, 2);
1518 T2_Pi01Pi02 = OrbitalTensors(Pi01Pi02, Pi01, Pi02, rRes, 2);
1519
1520 // X->YP Orbital
1521 vector<double> T1_PipPimPi01; T1_PipPimPi01.clear();
1522 vector<double> T1_PipPimPi02; T1_PipPimPi02.clear();
1523 vector<double> T1_PipPi01Pi02; T1_PipPi01Pi02.clear();
1524 vector<double> T1_PipPi02Pi01; T1_PipPi02Pi01.clear();
1525 vector<double> T1_PimPi01Pi02; T1_PimPi01Pi02.clear();
1526 vector<double> T1_PimPi02Pi01; T1_PimPi02Pi01.clear();
1527 vector<double> T1_PipPi01Pim; T1_PipPi01Pim.clear();
1528 vector<double> T1_PipPi02Pim; T1_PipPi02Pim.clear();
1529 vector<double> T1_PimPi01Pip; T1_PimPi01Pip.clear();
1530 vector<double> T1_PimPi02Pip; T1_PimPi02Pip.clear();
1531 vector<double> T1_Pi01Pi02Pip; T1_Pi01Pi02Pip.clear();
1532 vector<double> T1_Pi01Pi02Pim; T1_Pi01Pi02Pim.clear();
1533
1534 T1_PipPimPi01 = OrbitalTensors(PipPimPi01, PipPim, Pi01, rRes, 1);
1535 T1_PipPimPi02 = OrbitalTensors(PipPimPi02, PipPim, Pi02, rRes, 1);
1536 T1_PipPi01Pi02 = OrbitalTensors(PipPi01Pi02, PipPi01, Pi02, rRes, 1);
1537 T1_PipPi02Pi01 = OrbitalTensors(PipPi01Pi02, PipPi02, Pi01, rRes, 1);
1538 T1_PimPi01Pi02 = OrbitalTensors(PimPi01Pi02, PimPi01, Pi02, rRes, 1);
1539 T1_PimPi02Pi01 = OrbitalTensors(PimPi01Pi02, PimPi02, Pi01, rRes, 1);
1540 T1_PipPi01Pim = OrbitalTensors(PipPimPi01, PipPi01, Pim, rRes, 1);
1541 T1_PipPi02Pim = OrbitalTensors(PipPimPi02, PipPi02, Pim, rRes, 1);
1542 T1_PimPi01Pip = OrbitalTensors(PipPimPi01, PimPi01, Pip, rRes, 1);
1543 T1_PimPi02Pip = OrbitalTensors(PipPimPi02, PimPi02, Pip, rRes, 1);
1544 T1_Pi01Pi02Pip = OrbitalTensors(PipPi01Pi02, Pi01Pi02, Pip, rRes, 1);
1545 T1_Pi01Pi02Pim = OrbitalTensors(PimPi01Pi02, Pi01Pi02, Pim, rRes, 1);
1546
1547 vector<double> T2_PipPimPi01; T2_PipPimPi01.clear();
1548 vector<double> T2_PipPimPi02; T2_PipPimPi02.clear();
1549 vector<double> T2_PipPi01Pi02; T2_PipPi01Pi02.clear();
1550 vector<double> T2_PipPi02Pi01; T2_PipPi02Pi01.clear();
1551 vector<double> T2_PimPi01Pi02; T2_PimPi01Pi02.clear();
1552 vector<double> T2_PimPi02Pi01; T2_PimPi02Pi01.clear();
1553 vector<double> T2_PipPi01Pim; T2_PipPi01Pim.clear();
1554 vector<double> T2_PipPi02Pim; T2_PipPi02Pim.clear();
1555 vector<double> T2_PimPi01Pip; T2_PimPi01Pip.clear();
1556 vector<double> T2_PimPi02Pip; T2_PimPi02Pip.clear();
1557 vector<double> T2_Pi01Pi02Pip; T2_Pi01Pi02Pip.clear();
1558 vector<double> T2_Pi01Pi02Pim; T2_Pi01Pi02Pim.clear();
1559
1560 T2_PipPimPi01 = OrbitalTensors(PipPimPi01, PipPim, Pi01, rRes, 2);
1561 T2_PipPimPi02 = OrbitalTensors(PipPimPi02, PipPim, Pi02, rRes, 2);
1562 T2_PipPi01Pi02 = OrbitalTensors(PipPi01Pi02, PipPi01, Pi02, rRes, 2);
1563 T2_PipPi02Pi01 = OrbitalTensors(PipPi01Pi02, PipPi02, Pi01, rRes, 2);
1564 T2_PimPi01Pi02 = OrbitalTensors(PimPi01Pi02, PimPi01, Pi02, rRes, 2);
1565 T2_PimPi02Pi01 = OrbitalTensors(PimPi01Pi02, PimPi02, Pi01, rRes, 2);
1566 T2_PipPi01Pim = OrbitalTensors(PipPimPi01, PipPi01, Pim, rRes, 2);
1567 T2_PipPi02Pim = OrbitalTensors(PipPimPi02, PipPi02, Pim, rRes, 2);
1568 T2_PimPi01Pip = OrbitalTensors(PipPimPi01, PimPi01, Pip, rRes, 2);
1569 T2_PimPi02Pip = OrbitalTensors(PipPimPi02, PimPi02, Pip, rRes, 2);
1570 T2_Pi01Pi02Pip = OrbitalTensors(PipPi01Pi02, Pi01Pi02, Pip, rRes, 2);
1571 T2_Pi01Pi02Pim = OrbitalTensors(PimPi01Pi02, Pi01Pi02, Pim, rRes, 2);
1572
1573 // D->XX Orbital
1574 vector<double> T1_2pm12; T1_2pm12.clear();
1575 vector<double> T1_2p1m2; T1_2p1m2.clear();
1576 vector<double> T1_2p2m1; T1_2p2m1.clear();
1577
1578 T1_2pm12 = OrbitalTensors(D0, PipPim, Pi01Pi02, rD, 1);
1579 T1_2p1m2 = OrbitalTensors(D0, PipPi01, PimPi02, rD, 1);
1580 T1_2p2m1 = OrbitalTensors(D0, PipPi02, PimPi01, rD, 1);
1581
1582 vector<double> T2_2pm12; T2_2pm12.clear();
1583 vector<double> T2_2p1m2; T2_2p1m2.clear();
1584 vector<double> T2_2p2m1; T2_2p2m1.clear();
1585
1586 T2_2pm12 = OrbitalTensors(D0, PipPim, Pi01Pi02, rD, 2);
1587 T2_2p1m2 = OrbitalTensors(D0, PipPi01, PimPi02, rD, 2);
1588 T2_2p2m1 = OrbitalTensors(D0, PipPi02, PimPi01, rD, 2);
1589
1590 // D->XP Orbital
1591 vector<double> T1_3pm; T1_3pm.clear();
1592 vector<double> T1_3mp; T1_3mp.clear();
1593 vector<double> T1_3z12; T1_3z12.clear();
1594 vector<double> T1_3z21; T1_3z21.clear();
1595
1596 T1_3pm = OrbitalTensors(D0, PipPi01Pi02, Pim, rD, 1);
1597 T1_3mp = OrbitalTensors(D0, PimPi01Pi02, Pip, rD, 1);
1598 T1_3z12 = OrbitalTensors(D0, PipPimPi01, Pi02, rD, 1);
1599 T1_3z21 = OrbitalTensors(D0, PipPimPi02, Pi01, rD, 1);
1600
1601 vector<double> T2_3pm; T2_3pm.clear();
1602 vector<double> T2_3mp; T2_3mp.clear();
1603 vector<double> T2_3z12; T2_3z12.clear();
1604 vector<double> T2_3z21; T2_3z21.clear();
1605
1606 T2_3pm = OrbitalTensors(D0, PipPi01Pi02, Pim, rD, 2);
1607 T2_3mp = OrbitalTensors(D0, PimPi01Pi02, Pip, rD, 2);
1608 T2_3z12 = OrbitalTensors(D0, PipPimPi01, Pi02, rD, 2);
1609 T2_3z21 = OrbitalTensors(D0, PipPimPi02, Pi01, rD, 2);
1610
1611 complex<double> amplitude(0,0);
1612
1613 // D0 -> a1(1260)+ {rho(770)+ pi0[S]} pi-
1614 double SF_Ap_S_Vp1P = contract_11_0(contract_21_1(Proj1_3p, T1_PipPi01), T1_3pm);
1615 double SF_Ap_S_Vp2P = contract_11_0(contract_21_1(Proj1_3p, T1_PipPi02), T1_3pm);
1616
1617 amplitude += fitpara[0]*(SF_Ap_S_Vp1P*RBW_a11260_p*GS_rho770_p1 + SF_Ap_S_Vp2P*RBW_a11260_p*GS_rho770_p2);
1618
1619 // D0 -> a1(1260)+ {rho(770)+ pi0[D]} pi-
1620 double SF_Ap_D_Vp1P = contract_11_0(contract_21_1(T2_PipPi01Pi02, T1_PipPi01), T1_3pm);
1621 double SF_Ap_D_Vp2P = contract_11_0(contract_21_1(T2_PipPi02Pi01, T1_PipPi02), T1_3pm);
1622
1623// cout<<SF_Ap_D_VP_1<<","<<SF_Ap_D_VP_2<<","<<SF_Ap_D_VP_3<<","<<SF_Ap_D_VP_4<<","<<endl;
1624// cout<<"-------"<<endl;
1625 amplitude += fitpara[1]*(SF_Ap_D_Vp1P*RBW_a11260_p*GS_rho770_p1 + SF_Ap_D_Vp2P*RBW_a11260_p*GS_rho770_p2);
1626
1627 // D0 -> a1(1260)+ {f2(1270) pi+ [P]} pi-
1628 double SF_Ap_P_TP = contract_11_0(contract_21_1(contract_42_2(Proj2_3p, T2_Pi01Pi02), T1_Pi01Pi02Pip), T1_3pm);
1629
1630 amplitude += fitpara[2]*(SF_Ap_P_TP*RBW_a11260_p*RBW_f21270_00);
1631
1632 // D0 -> a1(1260)+ {f0 pi+ [P]} pi-
1633 double SF_Ap_P_SP = contract_11_0(T1_3pm, T1_Pi01Pi02Pip);
1634
1635 amplitude += fitpara[3]*(SF_Ap_P_SP*RBW_a11260_p*PiPiS_00_0);
1636 amplitude += fitpara[4]*(SF_Ap_P_SP*RBW_a11260_p*PiPiS_00_1);
1637 amplitude += fitpara[5]*(SF_Ap_P_SP*RBW_a11260_p*PiPiS_00_5);
1638
1639 // D0 -> a1(1260)- { rho(770)- pi0 [S]} pi+
1640 double SF_Am_S_Vm1P = contract_11_0(contract_21_1(Proj1_3m, T1_PimPi01), T1_3mp);
1641 double SF_Am_S_Vm2P = contract_11_0(contract_21_1(Proj1_3m, T1_PimPi02), T1_3mp);
1642
1643 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);
1644
1645 // D0 -> a1(1260)- {rho(770)- pi0[D]} pi+
1646 double SF_Am_D_Vm1P = contract_11_0(contract_21_1(T2_PimPi01Pi02, T1_PimPi01), T1_3mp);
1647 double SF_Am_D_Vm2P = contract_11_0(contract_21_1(T2_PimPi02Pi01, T1_PimPi02), T1_3mp);
1648
1649 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);
1650
1651 // D0 -> a1(1260)- {f2(1270) pi- [P]} pi+
1652 double SF_Am_P_TP = contract_11_0(contract_21_1(contract_42_2(Proj2_3m, T2_Pi01Pi02), T1_Pi01Pi02Pim), T1_3mp);
1653
1654 amplitude += fitpara[6]*fitpara[2]*(SF_Am_P_TP*RBW_a11260_m*RBW_f21270_00);
1655
1656 // D0 -> a1(1260)- {f0 pi- [P]} pi+
1657 double SF_Am_P_SP = contract_11_0(T1_3mp, T1_Pi01Pi02Pim);
1658
1659 amplitude += fitpara[6]*fitpara[3]*(SF_Am_P_SP*RBW_a11260_m*PiPiS_00_0);
1660 amplitude += fitpara[6]*fitpara[4]*(SF_Am_P_SP*RBW_a11260_m*PiPiS_00_1);
1661 amplitude += fitpara[6]*fitpara[5]*(SF_Am_P_SP*RBW_a11260_m*PiPiS_00_5);
1662
1663 // D -> a1(1260)0 pi0
1664 double SF_A01_S_Vp1P = contract_11_0(contract_21_1(Proj1_3z1, T1_PipPi01), T1_3z12);
1665 double SF_A02_S_Vp2P = contract_11_0(contract_21_1(Proj1_3z2, T1_PipPi02), T1_3z21);
1666 double SF_A01_S_Vm1P = contract_11_0(contract_21_1(Proj1_3z1, T1_PimPi01), T1_3z12);
1667 double SF_A02_S_Vm2P = contract_11_0(contract_21_1(Proj1_3z2, T1_PimPi02), T1_3z21);
1668 double SF_A01_S_VzP = contract_11_0(contract_21_1(Proj1_3z1, T1_PipPim), T1_3z12);
1669 double SF_A02_S_VzP = contract_11_0(contract_21_1(Proj1_3z2, T1_PipPim), T1_3z21);
1670
1671 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 + SF_A01_S_Vm1P*RBW_a11260_01*GS_rho770_m1 + SF_A02_S_Vm2P*RBW_a11260_02*GS_rho770_m2);
1672
1673 double SF_A01_D_Vp1P = contract_11_0(contract_21_1(T2_PipPi01Pim, T1_PipPi01), T1_3z12);
1674 double SF_A02_D_Vp2P = contract_11_0(contract_21_1(T2_PipPi02Pim, T1_PipPi02), T1_3z21);
1675 double SF_A01_D_Vm1P = contract_11_0(contract_21_1(T2_PimPi01Pip, T1_PimPi01), T1_3z12);
1676 double SF_A02_D_Vm2P = contract_11_0(contract_21_1(T2_PimPi02Pip, T1_PimPi02), T1_3z21);
1677
1678 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 + SF_A01_D_Vm1P*RBW_a11260_01*GS_rho770_m1 + SF_A02_D_Vm2P*RBW_a11260_02*GS_rho770_m2);
1679
1680 double SF_A01_P_TP = contract_11_0(contract_21_1(contract_42_2(Proj2_3z1,T2_PipPim), T1_PipPimPi01), T1_3z12);
1681 double SF_A02_P_TP = contract_11_0(contract_21_1(contract_42_2(Proj2_3z2,T2_PipPim), T1_PipPimPi02), T1_3z21);
1682
1683 amplitude += fitpara[7]*fitpara[2]*(-1.0)*(SF_A01_P_TP*RBW_a11260_01*RBW_f21270_pm + SF_A02_P_TP*RBW_a11260_02*RBW_f21270_pm);
1684
1685 double SF_A01_P_SP = contract_11_0(T1_3z12, T1_PipPimPi01);
1686 double SF_A02_P_SP = contract_11_0(T1_3z21, T1_PipPimPi02);
1687
1688 amplitude += fitpara[7]*fitpara[3]*(-1.0)*(SF_A01_P_SP*RBW_a11260_01*PiPiS_pm_0 + SF_A02_P_SP*RBW_a11260_02*PiPiS_pm_0);
1689 amplitude += fitpara[7]*fitpara[4]*(-1.0)*(SF_A01_P_SP*RBW_a11260_01*PiPiS_pm_1 + SF_A02_P_SP*RBW_a11260_02*PiPiS_pm_1);
1690 amplitude += fitpara[7]*fitpara[5]*(-1.0)*(SF_A01_P_SP*RBW_a11260_01*PiPiS_pm_5 + SF_A02_P_SP*RBW_a11260_02*PiPiS_pm_5);
1691
1692 // D0 -> a1(1420)+ {rho(770)+ pi0[S]} pi-
1693// amplitude += fitpara[8]*(SF_Ap_S_Vp1P*RBW_a11420_p*GS_rho770_p1 + SF_Ap_S_Vp2P*RBW_a11420_p*GS_rho770_p2);
1694// amplitude += fitpara[9]*(SF_A01_S_Vp1P*RBW_a11420_01*GS_rho770_p1 + SF_A02_S_Vp2P*RBW_a11420_02*GS_rho770_p2 + SF_A01_S_Vm1P*RBW_a11420_01*GS_rho770_m1 + SF_A02_S_Vm2P*RBW_a11420_02*GS_rho770_m2);
1695
1696 // D0 -> a1(1420)+ {f0 pi0+[S]} pi-
1697 amplitude += fitpara[8]*(SF_Ap_P_SP*RBW_a11420_p*PiPiS_00_5);
1698 amplitude += fitpara[9]*(SF_Ap_P_SP*RBW_a11420_p*PiPiS_00_6);
1699
1700 // D0 -> a2(1320)+ {rho+ pi0} pi-
1701 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)), PipPi01Pi02), contract_42_2(Proj2_3p, T2_3pm)), T2_PipPi01Pi02);
1702 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)), PipPi01Pi02), contract_42_2(Proj2_3p, T2_3pm)), T2_PipPi02Pi01);
1703
1704 amplitude += fitpara[10]*(SF_Tp_D_Vp1P*GS_rho770_p1*RBW_a21320_p + SF_Tp_D_Vp2P*GS_rho770_p2*RBW_a21320_p);
1705
1706 // D0 -> a2(1320)- {rho- pi0} pi+
1707 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)), PimPi01Pi02), contract_42_2(Proj2_3m, T2_3mp)), T2_PimPi01Pi02);
1708 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)), PimPi01Pi02), contract_42_2(Proj2_3m, T2_3mp)), T2_PimPi02Pi01);
1709 amplitude += fitpara[11]*(SF_Tm_D_Vm1P*GS_rho770_m1*RBW_a21320_m + SF_Tm_D_Vm2P*GS_rho770_m2*RBW_a21320_m);
1710
1711 // D0 -> h1(1170)0 {rho pi0} pi0
1712 amplitude += fitpara[12]*(SF_A01_S_Vp1P*RBW_h11170_01*GS_rho770_p1 + SF_A02_S_Vp2P*RBW_h11170_02*GS_rho770_p2 - 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*GS_rho770_pm - SF_A02_S_VzP*RBW_h11170_02*GS_rho770_pm);
1713
1714 // D0 -> pi(1300)- {rho(770)- pi0} pi+
1715 double SF_Pm_P_Vm1P = contract_11_0(T1_PimPi01,T1_PimPi01Pi02);
1716 double SF_Pm_P_Vm2P = contract_11_0(T1_PimPi02,T1_PimPi02Pi01);
1717
1718 amplitude += fitpara[13]*(SF_Pm_P_Vm1P*GS_rho770_m1*RBW_pi1300_m + SF_Pm_P_Vm2P*GS_rho770_m2*RBW_pi1300_m);
1719
1720 // D0 -> pi(1300)- {f2 pi-} pi+
1721// double SF_Pm_D_TP = contract_22_0(T2_Pi01Pi02,T2_Pi01Pi02Pim);
1722// amplitude += fitpara[14]*fitpara[13]*(SF_Pm_D_TP*RBW_f21270_00*RBW_pi1300_m);
1723
1724 // D0 -> pi(1300)- {f0 pi-} pi+
1725 amplitude += fitpara[14]*fitpara[13]*(RBW_pi1300_m*PiPiS_00_0);
1726// amplitude += fitpara[15]*fitpara[13]*(RBW_pi1300_m*PiPiS_00_5);
1727 amplitude += fitpara[15]*fitpara[13]*(RBW_pi1300_m*PiPiS_00_6);
1728
1729 // D0 -> pi(1300)+ {rho(770)+ pi0} pi-
1730 double SF_Pp_P_Vp1P = contract_11_0(T1_PipPi01,T1_PipPi01Pi02);
1731 double SF_Pp_P_Vp2P = contract_11_0(T1_PipPi02,T1_PipPi02Pi01);
1732
1733 amplitude += fitpara[16]*(SF_Pp_P_Vp1P*GS_rho770_p1*RBW_pi1300_p + SF_Pp_P_Vp2P*GS_rho770_p2*RBW_pi1300_p);
1734
1735 // D0 -> pi(1300)+ {f2 pi-} pi+
1736// double SF_Pp_D_TP = contract_22_0(T2_Pi01Pi02,T2_Pi01Pi02Pip);
1737// amplitude += fitpara[14]*fitpara[17]*(SF_Pp_D_TP*RBW_f21270_00*RBW_pi1300_p);
1738
1739 // D0 -> pi(1300)+ {f0 pi+} pi-
1740 amplitude += fitpara[14]*fitpara[16]*(RBW_pi1300_p*PiPiS_00_0);
1741// amplitude += fitpara[15]*fitpara[17]*(RBW_pi1300_p*PiPiS_00_5);
1742 amplitude += fitpara[15]*fitpara[16]*(RBW_pi1300_p*PiPiS_00_6);
1743
1744 // D0 -> pi(1300)0 {rho pi} pi0
1745 double SF_P01_P_Vp1P = contract_11_0(T1_PipPi01,T1_PipPi01Pim);
1746 double SF_P02_P_Vp2P = contract_11_0(T1_PipPi02,T1_PipPi02Pim);
1747 double SF_P01_P_Vm1P = contract_11_0(T1_PimPi01,T1_PimPi01Pip);
1748 double SF_P02_P_Vm2P = contract_11_0(T1_PimPi02,T1_PimPi02Pip);
1749
1750 amplitude += fitpara[17]*(SF_P01_P_Vp1P*RBW_pi1300_01*GS_rho770_p1 + SF_P02_P_Vp2P*RBW_pi1300_02*GS_rho770_p2 + SF_P01_P_Vm1P*RBW_pi1300_01*GS_rho770_m1 + SF_P02_P_Vm2P*RBW_pi1300_02*GS_rho770_m2);
1751
1752
1753// double SF_P01_D_TP = contract_22_0(T2_PipPim,T2_PipPimPi01);
1754// double SF_P02_D_TP = contract_22_0(T2_PipPim,T2_PipPimPi02);
1755// 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);
1756
1757 amplitude += fitpara[14]*fitpara[17]*(-1.0)*(RBW_pi1300_01*PiPiS_pm_0 + RBW_pi1300_02*PiPiS_pm_0);
1758// amplitude += fitpara[15]*fitpara[18]*(-1.0)*(RBW_pi1300_01*PiPiS_pm_5 + RBW_pi1300_02*PiPiS_pm_5);
1759 amplitude += fitpara[15]*fitpara[17]*(-1.0)*(RBW_pi1300_01*PiPiS_pm_6 + RBW_pi1300_02*PiPiS_pm_6);
1760
1761 // D0 -> rho+ rho- [S]
1762 double SF_Vp1Vm2_S = contract_11_0(T1_PipPi01, T1_PimPi02);
1763 double SF_Vp2Vm1_S = contract_11_0(T1_PipPi02, T1_PimPi01);
1764
1765 amplitude += fitpara[18]*(SF_Vp1Vm2_S*GS_rho770_p1*GS_rho770_m2 + SF_Vp2Vm1_S*GS_rho770_p2*GS_rho770_m1);
1766
1767 // D0 -> rho+ rho- [P]
1768 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);
1769 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);
1770
1771 amplitude += fitpara[19]*(SF_Vp1Vm2_P*GS_rho770_p1*GS_rho770_m2 + SF_Vp2Vm1_P*GS_rho770_p2*GS_rho770_m1);
1772
1773 // D0 -> rho+ rho- [D]
1774 double SF_Vp1Vm2_D = contract_11_0(contract_21_1(T2_2p1m2,T1_PipPi01), T1_PimPi02);
1775 double SF_Vp2Vm1_D = contract_11_0(contract_21_1(T2_2p2m1,T1_PipPi02), T1_PimPi01);
1776 amplitude += fitpara[20]*(SF_Vp1Vm2_D*GS_rho770_p1*GS_rho770_m2 + SF_Vp2Vm1_D*GS_rho770_p2*GS_rho770_m1);
1777
1778 // D0 -> rho0 (Pi0 Pi0)S
1779 double SF_VpmS12_P = contract_11_0(T1_PipPim,T1_2pm12);
1780
1781 amplitude += fitpara[21]*(SF_VpmS12_P*GS_rho770_pm*PiPiS_00_0);
1782 amplitude += fitpara[22]*(SF_VpmS12_P*GS_rho770_pm*PiPiS_00_5);
1783 amplitude += fitpara[23]*(SF_VpmS12_P*GS_rho770_pm*PiPiS_00_6);
1784
1785 // D0 -> f0f0
1786 //00
1787 amplitude += fitpara[24]*(PiPiS_pm_0*PiPiS_00_0 + PiPiS_00_0*PiPiS_pm_0);
1788 amplitude += fitpara[25]*(PiPiS_pm_0*PiPiS_00_1 + PiPiS_00_0*PiPiS_pm_1);
1789 amplitude += fitpara[26]*(PiPiS_pm_1*PiPiS_00_1 + PiPiS_00_1*PiPiS_pm_1);
1790 amplitude += fitpara[27]*(PiPiS_pm_1*PiPiS_00_5 + PiPiS_00_1*PiPiS_pm_5);
1791 amplitude += fitpara[28]*(PiPiS_pm_5*PiPiS_00_5 + PiPiS_00_5*PiPiS_pm_5);
1792 amplitude += fitpara[29]*(PiPiS_pm_5*PiPiS_00_6 + PiPiS_00_5*PiPiS_pm_6);
1793
1794 // D0 -> f2(1270) f0
1795 double SF_TpmS00_D = contract_22_0(T2_PipPim, T2_2pm12);
1796 double SF_T00Spm_D = contract_22_0(T2_Pi01Pi02, T2_2pm12);
1797
1798 amplitude += fitpara[30]*(SF_TpmS00_D*RBW_f21270_pm*PiPiS_00_5 + SF_T00Spm_D*RBW_f21270_00*PiPiS_pm_5);
1799 amplitude += fitpara[31]*(SF_TpmS00_D*RBW_f21270_pm*PiPiS_00_6 + SF_T00Spm_D*RBW_f21270_00*PiPiS_pm_6);
1800
1801 // D0 -> pi2(1670)0[f2(1270) pi0] pi0
1802 double SF_PT01_S_TP = contract_22_0(contract_42_2(Proj2_3z1, T2_PipPim), T2_3z12);
1803 double SF_PT02_S_TP = contract_22_0(contract_42_2(Proj2_3z2, T2_PipPim), T2_3z21);
1804
1805 amplitude += fitpara[32]*(-1.0)*(SF_PT01_S_TP*RBW_f21270_pm*RBW_pi21670_01 + SF_PT02_S_TP*RBW_f21270_pm*RBW_pi21670_02);
1806
1807 // D -> 1--(rho pi) pi0
1808 double SF_V1_Vz = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, PipPimPi01), T1_PipPimPi01), T1_PipPim), contract_21_1(Proj1_3z1, T1_3z12));
1809 double SF_V1_Vp1 = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, PipPimPi01), T1_PipPi01Pim), T1_PipPi01), contract_21_1(Proj1_3z1, T1_3z12));
1810 double SF_V1_Vm1 = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, PipPimPi01), T1_PimPi01Pip), T1_PimPi01), contract_21_1(Proj1_3z1, T1_3z12));
1811
1812 double SF_V2_Vz = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, PipPimPi02), T1_PipPimPi02), T1_PipPim), contract_21_1(Proj1_3z2, T1_3z21));
1813 double SF_V2_Vp2 = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, PipPimPi02), T1_PipPi02Pim), T1_PipPi02), contract_21_1(Proj1_3z2, T1_3z21));
1814 double SF_V1_Vm2 = contract_11_0(contract_21_1(contract_31_2(contract_41_3(epsilon_uvmn, PipPimPi02), T1_PimPi02Pip), T1_PimPi02), contract_21_1(Proj1_3z2, T1_3z21));
1815
1816
1817 // D0 -> omega(rho pi) pi0
1818 amplitude += (-1.0) * fitpara[33]*(SF_V1_Vp1*RBW_omega_01*GS_rho770_p1 - SF_V1_Vz*RBW_omega_01*GS_rho770_pm - SF_V1_Vm1*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*RBW_omega_02*GS_rho770_m2);
1819
1820
1821 // D0 -> phi(rho pi) pi0
1822 amplitude += (-1.0) * fitpara[34]*(SF_V1_Vp1*RBW_phi_01*GS_rho770_p1 - SF_V1_Vz*RBW_phi_01*GS_rho770_pm - SF_V1_Vm1*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*GS_rho770_m2);
1823
1824 return amplitude;
1825
1826}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
P
double P(RecMdcKalTrack *trk)
Definition: CalibEventSelect.cxx:87
mass
double mass
Definition: CosmicGenerator.cxx:138
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
f1
TFile * f1
Definition: DataBase/tau_mode.c:4
g1
TF1 * g1
Definition: DataBase/tau_mode.c:66
ID
int ID[no]
Definition: DocaEangle/fitDocaeangle.h:2
eps
EvtTensor3C eps(const EvtVector3R &v)
Definition: EvtTensor3C.cc:307
output
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per saves r n generator level $ !Flag for chat level in output
Definition: FoamA.h:89
mpi
const double mpi
Definition: Gam4pikp.cxx:47
s
XmlRpcServer s
Definition: HelloServer.cpp:11
q
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition: KKsem.h:33
t
TTree * t
Definition: binning.cxx:23
D0Topippim2pi0::~D0Topippim2pi0
virtual ~D0Topippim2pi0()
Definition: D0Topippim2pi0.cxx:11
D0Topippim2pi0::Amp
complex< double > Amp(vector< double > Pip, vector< double > Pim, vector< double > Pi01, vector< double > Pi02)
Definition: D0Topippim2pi0.cxx:1384
m2
double double * m2
Definition: qcdloop1.h:75