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