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