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