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