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