BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtD0Topipipi0.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: EvtD0Topipipi0.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 "EvtD0Topipipi0.hh"
21#include "../EvtGenBase/EvtGenKine.hh"
22#include "../EvtGenBase/EvtHelSys.hh"
23#include "../EvtGenBase/EvtPDL.hh"
24#include "../EvtGenBase/EvtParticle.hh"
25#include "../EvtGenBase/EvtPatches.hh"
26#include "../EvtGenBase/EvtReport.hh"
27#include "../EvtGenBase/EvtTensor4C.hh"
28#include "../EvtGenBase/EvtVector3R.hh"
29#include "../EvtGenBase/EvtVector4C.hh"
30#include "../EvtGenBase/EvtVector4R.hh"
31#include <complex>
32#include <fstream>
33#include <iostream>
34#include <map>
35#include <math.h>
36#include <stdlib.h>
37#include <string>
38#include <vector>
39using namespace std;
40#define PI acos( -1 )
41
42EvtD0Topipipi0::~EvtD0Topipipi0() {}
43
44void EvtD0Topipipi0::getName( std::string& model_name ) { model_name = "D0Topipipi0"; }
45
46EvtDecayBase* EvtD0Topipipi0::clone() { return new EvtD0Topipipi0; }
47
48void EvtD0Topipipi0::init() {
49 checkNArg( 2 );
50 checkNDaug( 3 );
51 charm = getArg( 0 );
52 tagmode = getArg( 1 );
53 // std::cout << "Initializing EvtD0Topipipi0: charm= "<<charm<<" tagmode=
54 // "<<tagmode<<std::endl;
55
56 g_uv.clear();
57 for ( int i = 0; i < 4; i++ )
58 {
59 for ( int j = 0; j < 4; j++ )
60 {
61 if ( i != j ) { g_uv.push_back( 0.0 ); }
62 else if ( i < 3 ) { g_uv.push_back( -1.0 ); }
63 else if ( i == 3 ) { g_uv.push_back( 1.0 ); }
64 }
65 }
66
67 epsilon_uvmn.clear();
68 for ( int i = 0; i < 4; i++ )
69 {
70 for ( int j = 0; j < 4; j++ )
71 {
72 for ( int k = 0; k < 4; k++ )
73 {
74 for ( int l = 0; l < 4; l++ )
75 {
76 if ( i == j || i == k || i == l || j == k || j == l || k == l )
77 { epsilon_uvmn.push_back( 0.0 ); }
78 else
79 {
80 if ( i == 0 && j == 1 && k == 2 && l == 3 ) epsilon_uvmn.push_back( 1.0 );
81 if ( i == 0 && j == 1 && k == 3 && l == 2 ) epsilon_uvmn.push_back( -1.0 );
82 if ( i == 0 && j == 2 && k == 1 && l == 3 ) epsilon_uvmn.push_back( -1.0 );
83 if ( i == 0 && j == 2 && k == 3 && l == 1 ) epsilon_uvmn.push_back( 1.0 );
84 if ( i == 0 && j == 3 && k == 1 && l == 2 ) epsilon_uvmn.push_back( 1.0 );
85 if ( i == 0 && j == 3 && k == 2 && l == 1 ) epsilon_uvmn.push_back( -1.0 );
86
87 if ( i == 1 && j == 0 && k == 2 && l == 3 ) epsilon_uvmn.push_back( -1.0 );
88 if ( i == 1 && j == 0 && k == 3 && l == 2 ) epsilon_uvmn.push_back( 1.0 );
89 if ( i == 1 && j == 2 && k == 0 && l == 3 ) epsilon_uvmn.push_back( 1.0 );
90 if ( i == 1 && j == 2 && k == 3 && l == 0 ) epsilon_uvmn.push_back( -1.0 );
91 if ( i == 1 && j == 3 && k == 0 && l == 2 ) epsilon_uvmn.push_back( -1.0 );
92 if ( i == 1 && j == 3 && k == 2 && l == 0 ) epsilon_uvmn.push_back( 1.0 );
93
94 if ( i == 2 && j == 0 && k == 1 && l == 3 ) epsilon_uvmn.push_back( 1.0 );
95 if ( i == 2 && j == 0 && k == 3 && l == 1 ) epsilon_uvmn.push_back( -1.0 );
96 if ( i == 2 && j == 1 && k == 0 && l == 3 ) epsilon_uvmn.push_back( -1.0 );
97 if ( i == 2 && j == 1 && k == 3 && l == 0 ) epsilon_uvmn.push_back( 1.0 );
98 if ( i == 2 && j == 3 && k == 0 && l == 1 ) epsilon_uvmn.push_back( 1.0 );
99 if ( i == 2 && j == 3 && k == 1 && l == 0 ) epsilon_uvmn.push_back( -1.0 );
100
101 if ( i == 3 && j == 0 && k == 1 && l == 2 ) epsilon_uvmn.push_back( -1.0 );
102 if ( i == 3 && j == 0 && k == 2 && l == 1 ) epsilon_uvmn.push_back( 1.0 );
103 if ( i == 3 && j == 1 && k == 0 && l == 2 ) epsilon_uvmn.push_back( 1.0 );
104 if ( i == 3 && j == 1 && k == 2 && l == 0 ) epsilon_uvmn.push_back( -1.0 );
105 if ( i == 3 && j == 2 && k == 0 && l == 1 ) epsilon_uvmn.push_back( -1.0 );
106 if ( i == 3 && j == 2 && k == 1 && l == 0 ) epsilon_uvmn.push_back( 1.0 );
107 }
108 }
109 }
110 }
111 }
112
113 _nd = 3;
114 math_pi = 3.1415926;
115 mass_Pion = 0.13957;
116
117 rRes = 3.0 * 0.197321;
118 rD = 5.0 * 0.197321;
119 m_Pi = mass_Pion;
120 m2_Pi = m_Pi * m_Pi;
121 m_Pi0 = 0.134977;
122 m2_Pi0 = m_Pi0 * m_Pi0;
123
124 m0_rho7700 = 0.77526;
125 w0_rho7700 = 0.1478;
126
127 m0_rho770p = 0.77511;
128 w0_rho770p = 0.1491;
129
130 m0_f21270 = 1.2755;
131 w0_f21270 = 0.1867;
132
133 s0_prod = -5.0;
134
135 fitpara.clear();
136 fitpara.push_back( complex<double>( 100, 0 ) );
137 fitpara.push_back( complex<double>( 69.8939 * cos( 3.14983 ), 69.8939 * sin( 3.14983 ) ) );
138 fitpara.push_back( complex<double>( 58.521 * cos( -2.90685 ), 58.521 * sin( -2.90685 ) ) );
139 fitpara.push_back(
140 complex<double>( 483.035 * cos( -0.679009 ), 483.035 * sin( -0.679009 ) ) );
141 fitpara.push_back(
142 complex<double>( 441.921 * cos( -0.879847 ), 441.921 * sin( -0.879847 ) ) );
143 fitpara.push_back(
144 complex<double>( 1356.95 * cos( -0.206653 ), 1356.95 * sin( -0.206653 ) ) );
145 fitpara.push_back( complex<double>( 559.218 * cos( 0.501728 ), 559.218 * sin( 0.501728 ) ) );
146 fitpara.push_back( complex<double>( 3165.25 * cos( 3.39939 ), 3165.25 * sin( 3.39939 ) ) );
147 fitpara.push_back( complex<double>( 1422.93 * cos( 3.05347 ), 1422.93 * sin( 3.05347 ) ) );
148 fitpara.push_back( complex<double>( 2399.8 * cos( 2.24983 ), 2399.8 * sin( 2.24983 ) ) );
149 fitpara.push_back( complex<double>( 4601.19 * cos( 2.74388 ), 4601.19 * sin( 2.74388 ) ) );
150 fitpara.push_back( complex<double>( 1684.1 * cos( 1.99894 ), 1684.1 * sin( 1.99894 ) ) );
151 fitpara.push_back(
152 complex<double>( 678.674 * cos( -2.510691 ), 678.674 * sin( -2.510691 ) ) );
153 fitpara.push_back( complex<double>( 2.19068 * cos( 0.991805 ), 2.19068 * sin( 0.991805 ) ) );
154
155 return;
156}
157
158void EvtD0Topipipi0::initProbMax() { setProbMax( 60000000 ); }
159
160void EvtD0Topipipi0::decay( EvtParticle* p ) {
161 /*
162 double maxprob = 0.0;
163 for(int ir=0;ir<=60000000;ir++){
164 p->initializePhaseSpace(getNDaug(),getDaugs());
165 for(int i=0; i<_nd; i++){
166 _p4Lab[i]=p->getDaug(i)->getP4Lab();
167 _p4CM[i]=p->getDaug(i)->getP4();
168 }
169 double Prob = AmplitudeSquare(charm, tagmode);
170 if(Prob>maxprob) {
171 maxprob=Prob;
172 std::cout << "Max PDF = " << ir << " prob= " << Prob << std::endl;
173 }
174 }
175 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
176 */
177 p->initializePhaseSpace( getNDaug(), getDaugs() );
178 for ( int i = 0; i < _nd; i++ )
179 {
180 _p4Lab[i] = p->getDaug( i )->getP4Lab();
181 _p4CM[i] = p->getDaug( i )->getP4();
182 }
183 double prob = AmplitudeSquare( charm, tagmode );
184 setProb( prob );
185 return;
186}
187
188void EvtD0Topipipi0::setInput( double* pip, double* pim, double* pi0 ) {
189 p4_Pip.clear();
190 p4_Pim.clear();
191 p4_Pi0.clear();
192 p4_Pip.push_back( pip[0] );
193 p4_Pim.push_back( pim[0] );
194 p4_Pi0.push_back( pi0[0] );
195 p4_Pip.push_back( pip[1] );
196 p4_Pim.push_back( pim[1] );
197 p4_Pi0.push_back( pi0[1] );
198 p4_Pip.push_back( pip[2] );
199 p4_Pim.push_back( pim[2] );
200 p4_Pi0.push_back( pi0[2] );
201 p4_Pip.push_back( pip[3] );
202 p4_Pim.push_back( pim[3] );
203 p4_Pi0.push_back( pi0[3] );
204}
205
206vector<double> EvtD0Topipipi0::sum_tensor( vector<double> pa, vector<double> pb ) {
207 if ( pa.size() != pb.size() )
208 {
209 cout << "error sum tensor" << endl;
210 exit( 0 );
211 }
212 vector<double> temp;
213 temp.clear();
214 for ( int i = 0; i < pa.size(); i++ )
215 {
216 double sum = pa[i] + pb[i];
217 temp.push_back( sum );
218 }
219 return temp;
220}
221
222double EvtD0Topipipi0::contract_11_0( vector<double> pa, vector<double> pb ) {
223 if ( pa.size() != pb.size() || pa.size() != 4 )
224 {
225 cout << "error contract 11->0" << endl;
226 exit( 0 );
227 }
228 double temp = pa[3] * pb[3] - pa[0] * pb[0] - pa[1] * pb[1] - pa[2] * pb[2];
229 return temp;
230}
231
232vector<double> EvtD0Topipipi0::contract_21_1( vector<double> pa, vector<double> pb ) {
233 if ( pa.size() != 16 || pb.size() != 4 )
234 {
235 cout << "error contract 21->1" << endl;
236 exit( 0 );
237 }
238 vector<double> temp;
239 temp.clear();
240 for ( int i = 0; i < 4; i++ )
241 {
242 double sum = 0;
243 for ( int j = 0; j < 4; j++ )
244 {
245 int idx = i * 4 + j;
246 sum += pa[idx] * pb[j] * g_uv[4 * j + j];
247 }
248 temp.push_back( sum );
249 }
250 return temp;
251}
252
253double EvtD0Topipipi0::contract_22_0( vector<double> pa, vector<double> pb ) {
254 if ( pa.size() != pb.size() || pa.size() != 16 )
255 {
256 cout << "error contract 22->0" << endl;
257 exit( 0 );
258 }
259 double temp = 0;
260 for ( int i = 0; i < 4; i++ )
261 {
262 for ( int j = 0; j < 4; j++ )
263 {
264 int idx = i * 4 + j;
265 temp += pa[idx] * pb[idx] * g_uv[4 * i + i] * g_uv[4 * j + j];
266 }
267 }
268 return temp;
269}
270
271vector<double> EvtD0Topipipi0::contract_31_2( vector<double> pa, vector<double> pb ) {
272 if ( pa.size() != 64 || pb.size() != 4 )
273 {
274 cout << "error contract 31->2" << endl;
275 exit( 0 );
276 }
277 vector<double> temp;
278 temp.clear();
279 for ( int i = 0; i < 16; i++ )
280 {
281 double sum = 0;
282 for ( int j = 0; j < 4; j++ )
283 {
284 int idx = i * 4 + j;
285 sum += pa[idx] * pb[j] * g_uv[4 * j + j];
286 }
287 temp.push_back( sum );
288 }
289 return temp;
290}
291
292vector<double> EvtD0Topipipi0::contract_41_3( vector<double> pa, vector<double> pb ) {
293 if ( pa.size() != 256 || pb.size() != 4 )
294 {
295 cout << "error contract 41->3" << endl;
296 exit( 0 );
297 }
298 vector<double> temp;
299 temp.clear();
300 for ( int i = 0; i < 64; i++ )
301 {
302 double sum = 0;
303 for ( int j = 0; j < 4; j++ )
304 {
305 int idx = i * 4 + j;
306 sum += pa[idx] * pb[j] * g_uv[4 * j + j];
307 }
308 temp.push_back( sum );
309 }
310 return temp;
311}
312
313vector<double> EvtD0Topipipi0::contract_42_2( vector<double> pa, vector<double> pb ) {
314 if ( pa.size() != 256 || pb.size() != 16 )
315 {
316 cout << "error contract 42->2" << endl;
317 exit( 0 );
318 }
319 vector<double> temp;
320 temp.clear();
321 for ( int i = 0; i < 16; i++ )
322 {
323 double sum = 0;
324 for ( int j = 0; j < 4; j++ )
325 {
326 for ( int k = 0; k < 4; k++ )
327 {
328 int idxa = i * 16 + j * 4 + k;
329 int idxb = j * 4 + k;
330 sum += pa[idxa] * pb[idxb] * g_uv[4 * j + j] * g_uv[4 * k + k];
331 }
332 }
333 temp.push_back( sum );
334 }
335
336 return temp;
337}
338
339vector<double> EvtD0Topipipi0::contract_22_2( vector<double> pa, vector<double> pb ) {
340 if ( pa.size() != 16 || pb.size() != 16 )
341 {
342 cout << "error contract 42->2" << endl;
343 exit( 0 );
344 }
345 vector<double> temp;
346 temp.clear();
347 for ( int i = 0; i < 4; i++ )
348 {
349 for ( int j = 0; j < 4; j++ )
350 {
351 double sum = 0;
352 for ( int k = 0; k < 4; k++ )
353 {
354 int idxa = i * 4 + k;
355 int idxb = j * 4 + k;
356 sum += pa[idxa] * pb[idxb] * g_uv[4 * k + k];
357 }
358 temp.push_back( sum );
359 }
360 }
361
362 return temp;
363}
364
365/*
366 vector<double> EvtD0Topipipi0::contract_11_2(vector<double> pa, vector<double> pb){
367 if(pa.size()!=pb.size() || pa.size()!=4) {
368 cout<<"error contract 11->2"<<endl;
369 exit(0);
370 }
371 vector<double> temp; temp.clear();
372 for(int i=0; i<4; i++){
373 for(int j=0; j<4; j++){
374 double prod = pa[i]*pb[j];
375 temp.push_back(prod);
376 }
377 }
378 return temp;
379 }
380
381 vector<double> EvtD0Topipipi0::contract_22_4(vector<double> pa, vector<double> pb){
382 if(pa.size()!=pb.size() || pa.size()!=16) {
383 cout<<"error contract 22->4"<<endl;
384 exit(0);
385 }
386 vector<double> temp; temp.clear();
387 for(int i=0; i<16; i++){
388 for(int j=0; j<16; j++){
389 double prod = pa[i]*pb[j];
390 temp.push_back(prod);
391 }
392 }
393 return temp;
394 }
395 */
396
397// OrbitalTensors
398vector<double> EvtD0Topipipi0::OrbitalTensors( vector<double> pa, vector<double> pb,
399 vector<double> pc, double r, int rank ) {
400
401 if ( pa.size() != 4 || pb.size() != 4 || pc.size() != 4 )
402 {
403 cout << "Error: pa, pb, pc" << endl;
404 exit( 0 );
405 }
406 if ( rank < 0 )
407 {
408 cout << "Error: L<0 !!!" << endl;
409 exit( 0 );
410 }
411
412 // relative momentum
413 vector<double> mr;
414 mr.clear();
415
416 for ( int i = 0; i < 4; i++ )
417 {
418 double temp = pb[i] - pc[i];
419 mr.push_back( temp );
420 }
421
422 // "Masses" of particles
423 double msa = contract_11_0( pa, pa );
424 double msb = contract_11_0( pb, pb );
425 double msc = contract_11_0( pc, pc );
426
427 // Q^2_abc
428 double top = msa + msb - msc;
429 double Q2abc = top * top / ( 4.0 * msa ) - msb;
430
431 // Barrier factors
432 double Q_0 = 0.197321f / r;
433 double Q_02 = Q_0 * Q_0;
434 double Q_04 = Q_02 * Q_02;
435 // double Q_06 = Q_04*Q_02;
436 // double Q_08 = Q_04*Q_04;
437
438 double Q4abc = Q2abc * Q2abc;
439 // double Q6abc = Q4abc*Q2abc;
440 // double Q8abc = Q4abc*Q4abc;
441
442 double mB1 = sqrt( 2.0f / ( Q2abc + Q_02 ) );
443 double mB2 = sqrt( 13.0f / ( Q4abc + ( 3.0f * Q_02 ) * Q2abc + 9.0f * Q_04 ) );
444 // mB3 = &sqrt(277.0f/(Q6abc + (6.0f*Q_02)*Q4abc + (45.0f*Q_04)*Q2abc + 225.0f*Q_06));
445 // mB4 = &sqrt(12746.0f/(Q8abc + (10.0f*Q_02)*Q6abc + (135.0f*Q_04)*Q4abc +
446 //(1575.0f*Q_06)*Q2abc + 11025.0f*Q_08));
447
448 // Projection Operator 2-rank
449 vector<double> proj_uv;
450 proj_uv.clear();
451 for ( int i = 0; i < 4; i++ )
452 {
453 for ( int j = 0; j < 4; j++ )
454 {
455 int idx = i * 4 + j;
456 double temp = -g_uv[idx] + pa[i] * pa[j] / msa;
457 proj_uv.push_back( temp );
458 }
459 }
460
461 // Orbital Tensors
462 if ( rank == 0 )
463 {
464 vector<double> t;
465 t.clear();
466 t.push_back( 1.0 );
467 return t;
468 }
469 else if ( rank < 3 )
470 {
471 vector<double> t_u;
472 t_u.clear();
473 vector<double> Bt_u;
474 Bt_u.clear();
475 for ( int i = 0; i < 4; i++ )
476 {
477 double temp = 0;
478 for ( int j = 0; j < 4; j++ )
479 {
480 int idx = i * 4 + j;
481 temp += -proj_uv[idx] * mr[j] * g_uv[j * 4 + j];
482 }
483 t_u.push_back( temp );
484 Bt_u.push_back( temp * mB1 );
485 }
486 if ( rank == 1 ) return Bt_u;
487
488 double t_u2 = contract_11_0( t_u, t_u );
489
490 vector<double> Bt_uv;
491 Bt_uv.clear();
492 for ( int i = 0; i < 4; i++ )
493 {
494 for ( int j = 0; j < 4; j++ )
495 {
496 int idx = 4 * i + j;
497 double temp = t_u[i] * t_u[j] + ( 1.0 / 3.0 ) * proj_uv[idx] * t_u2;
498 Bt_uv.push_back( temp * mB2 );
499 }
500 }
501 if ( rank == 2 ) return Bt_uv;
502 }
503 else
504 {
505 cout << "rank>2: please add it by yourself!!!" << endl;
506 exit( 0 );
507 }
508
509 std::cerr << __FILE__ << ":" << __LINE__ << ": should not reach here" << std::endl;
510 exit( 1 );
511}
512
513// projection Tensor
514vector<double> EvtD0Topipipi0::ProjectionTensors( vector<double> pa, int rank ) {
515
516 if ( pa.size() != 4 )
517 {
518 cout << "Error: pa" << endl;
519 exit( 0 );
520 }
521 if ( rank < 0 )
522 {
523 cout << "Error: L<0 !!!" << endl;
524 exit( 0 );
525 }
526
527 double msa = contract_11_0( pa, pa );
528
529 // Projection Operator 2-rank
530 vector<double> proj_uv;
531 proj_uv.clear();
532 for ( int i = 0; i < 4; i++ )
533 {
534 for ( int j = 0; j < 4; j++ )
535 {
536 int idx = i * 4 + j;
537 double temp = -g_uv[idx] + pa[i] * pa[j] / msa;
538 proj_uv.push_back( temp );
539 }
540 }
541
542 // Orbital Tensors
543 if ( rank == 0 )
544 {
545 vector<double> t;
546 t.clear();
547 t.push_back( 1.0 );
548 return t;
549 }
550 else if ( rank == 1 ) { return proj_uv; }
551 else if ( rank == 2 )
552 {
553 vector<double> proj_uvmn;
554 proj_uvmn.clear();
555 for ( int i = 0; i < 4; i++ )
556 {
557 for ( int j = 0; j < 4; j++ )
558 {
559 for ( int k = 0; k < 4; k++ )
560 {
561 for ( int l = 0; l < 4; l++ )
562 {
563
564 int idx1_1 = 4 * i + k;
565 int idx1_2 = 4 * i + l;
566 int idx1_3 = 4 * i + j;
567
568 int idx2_1 = 4 * j + l;
569 int idx2_2 = 4 * j + k;
570 int idx2_3 = 4 * k + l;
571
572 double temp = ( 1.0 / 2.0 ) * ( proj_uv[idx1_1] * proj_uv[idx2_1] +
573 proj_uv[idx1_2] * proj_uv[idx2_2] ) -
574 ( 1.0 / 3.0 ) * proj_uv[idx1_3] * proj_uv[idx2_3];
575 proj_uvmn.push_back( temp );
576 }
577 }
578 }
579 }
580 return proj_uvmn;
581 }
582 else
583 {
584 cout << "rank>2: please add it by yourself!!!" << endl;
585 exit( 0 );
586 }
587}
588double EvtD0Topipipi0::fundecaymomentum( double mr2, double m1_2, double m2_2 ) {
589 double mr = sqrt( mr2 );
590 double poly = mr2 * mr2 + m1_2 * m1_2 + m2_2 * m2_2 - 2 * m1_2 * mr2 - 2 * m2_2 * mr2 -
591 2 * m1_2 * m2_2;
592 double ret = sqrt( poly ) / ( 2 * mr );
593 if ( poly < 0 )
594 // if(sqrt(m1_2) +sqrt(m2_2) > mr)
595 ret = 0.0f;
596 return ret;
597}
598
599complex<double> EvtD0Topipipi0::breitwigner( double mx2, double mr, double wr ) {
600 double output_x = 0;
601 double output_y = 0;
602
603 double mr2 = mr * mr;
604 double diff = mr2 - mx2;
605 double denom = diff * diff + wr * wr * mr2;
606 if ( wr < 0 )
607 {
608 output_x = 0;
609 output_y = 0;
610 }
611 else if ( wr < 10 )
612 {
613 output_x = diff / denom;
614 output_y = wr * mr / denom;
615 }
616 else
617 { /* phase space */
618 output_x = 1;
619 output_y = 0;
620 }
621 complex<double> output( output_x, output_y );
622 return output;
623}
624
625/* GS propagator */
626double EvtD0Topipipi0::h( double m, double q ) {
627 double h = 2.0 / math_pi * q / m * log( ( m + 2.0 * q ) / ( 2.0 * mass_Pion ) );
628 return h;
629}
630
631double EvtD0Topipipi0::dh( double m0, double q0 ) {
632 double dh = h( m0, q0 ) * ( 1.0 / ( 8.0 * q0 * q0 ) - 1.0 / ( 2.0 * m0 * m0 ) ) +
633 1.0 / ( 2.0 * math_pi * m0 * m0 );
634 return dh;
635}
636
637double EvtD0Topipipi0::f( double m0, double sx, double q0, double q ) {
638 double m = sqrt( sx );
639 double f =
640 m0 * m0 / ( q0 * q0 * q0 ) *
641 ( q * q * ( h( m, q ) - h( m0, q0 ) ) + ( m0 * m0 - sx ) * q0 * q0 * dh( m0, q0 ) );
642 return f;
643}
644
645double EvtD0Topipipi0::d( double m0, double q0 ) {
646 double d = 3.0 / math_pi * mass_Pion * mass_Pion / ( q0 * q0 ) *
647 log( ( m0 + 2.0 * q0 ) / ( 2.0 * mass_Pion ) ) +
648 m0 / ( 2.0 * math_pi * q0 ) -
649 ( mass_Pion * mass_Pion * m0 ) / ( math_pi * q0 * q0 * q0 );
650 return d;
651}
652
653double EvtD0Topipipi0::fundecaymomentum2( double mr2, double m1_2, double m2_2 ) {
654 double mr = sqrt( mr2 );
655 double poly = mr2 * mr2 + m1_2 * m1_2 + m2_2 * m2_2 - 2 * m1_2 * mr2 - 2 * m2_2 * mr2 -
656 2 * m1_2 * m2_2;
657 double ret = poly / ( 4.0f * mr2 );
658 if ( poly < 0 ) ret = 0.0f;
659 return ret;
660}
661
662double EvtD0Topipipi0::wid( double mass, double sa, double sb, double sc, double r, int l ) {
663 double widm = 1.0;
664 double sa0 = mass * mass;
665 double m = sqrt( sa );
666 double q = fundecaymomentum2( sa, sb, sc );
667 double q0 = fundecaymomentum2( sa0, sb, sc );
668 r = r / 0.197321;
669 double z = q * r * r;
670 double z0 = q0 * r * r;
671 double F = 0.0;
672 if ( l == 0 ) F = 1.0;
673 if ( l == 1 ) F = sqrt( ( 1.0 + z0 ) / ( 1.0 + z ) );
674 if ( l == 2 ) F = sqrt( ( 9.0 + 3.0 * z0 + z0 * z0 ) / ( 9.0 + 3.0 * z + z * z ) );
675 if ( l == 3 )
676 F = sqrt( ( 225.0 + 45.0 * z0 + 6.0 * z0 * z0 + z0 * z0 * z0 ) /
677 ( 225.0 + 45.0 * z + 6.0 * z * z + z * z * z ) );
678 if ( l == 4 )
679 F = sqrt(
680 ( 11025.0 + 1575.0 * z0 + 135.0 * z0 * z0 + 10.0 * z0 * z0 * z0 + z0 * z0 * z0 * z0 ) /
681 ( 11025.0 + 1575.0 * z + 135.0 * z * z + 10.0 * z * z * z + z * z * z * z ) );
682 double t = sqrt( q / q0 );
683 // printf("sa0 = %f, sb = %f, sc = %f, q = %f, q0 = %0.15f, qq0 = %f, t = %f
684 // \n",sa0,sb,sc,q,q0,q/q0,t);
685 uint i = 0;
686 for ( i = 0; i < ( 2 * l + 1 ); i++ ) { widm *= t; }
687 widm *= ( mass / m * F * F );
688 return widm;
689}
690
691/* for rho0, use GS, rho0->pi+ pi-, only angular momentum 1*/
692complex<double> EvtD0Topipipi0::GS( double mx2, double mr, double wr, double m1_2, double m2_2,
693 double r, int l ) {
694
695 double mr2 = mr * mr;
696 double q = fundecaymomentum( mx2, m1_2, m2_2 );
697 double q0 = fundecaymomentum( mr2, m1_2, m2_2 );
698 double numer = 1.0 + d( mr, q0 ) * wr / mr;
699 double denom_real = mr2 - mx2 + wr * f( mr, mx2, q0, q );
700 double denom_imag = mr * wr * wid( mr, mx2, m1_2, m2_2, r, l ); // real-i*imag;
701
702 double denom = denom_real * denom_real + denom_imag * denom_imag;
703 double output_x = denom_real * numer / denom;
704 double output_y = denom_imag * numer / denom;
705
706 complex<double> output( output_x, output_y );
707 return output;
708}
709
710complex<double> EvtD0Topipipi0::irho( double mr2, double m1_2, double m2_2 ) {
711 double mr = sqrt( mr2 );
712 double poly = mr2 * mr2 + m1_2 * m1_2 + m2_2 * m2_2 - 2 * m1_2 * mr2 - 2 * m2_2 * mr2 -
713 2 * m1_2 * m2_2;
714 double ret_real = 0;
715 double ret_imag = 0;
716 if ( poly >= 0 )
717 {
718 ret_real = 2.0f * sqrt( poly ) / ( 2.0f * mr2 );
719 ret_imag = 0;
720 }
721 else
722 {
723 ret_real = 0;
724 ret_imag = 2.0f * sqrt( -1.0f * poly ) / ( 2.0f * mr2 );
725 }
726 complex<double> ret( ret_real, ret_imag );
727 return ret;
728}
729
730complex<double> EvtD0Topipipi0::Flatte( double mx2, double mr, double g1, double g2,
731 double m1a, double m1b, double m2a, double m2b ) {
732
733 double mr2 = mr * mr;
734 double diff = mr2 - mx2;
735 double g22 = g2 * g1;
736 complex<double> ps1 = irho( mx2, m1a * m1a, m1b * m1b );
737 complex<double> ps2 = irho( mx2, m2a * m2a, m2b * m2b );
738 complex<double> iws = g1 * ps1 + g22 * ps2; /*mass dependent width */
739
740 double denom_real = diff + iws.imag();
741 double denom_imag = iws.real();
742 double denom = denom_real * denom_real + denom_imag * denom_imag;
743
744 double output_x = denom_real / denom;
745 double output_y = denom_imag / denom;
746
747 complex<double> output( output_x, output_y );
748 return output;
749}
750
751complex<double> EvtD0Topipipi0::RBW( double mx2, double mr, double wr, double m1_2,
752 double m2_2, double r, int l ) {
753 double mx = sqrt( mx2 );
754 double mr2 = mr * mr;
755 double denom_real = mr2 - mx2;
756 double denom_imag = 0;
757 if ( m1_2 > 0 && m2_2 > 0 )
758 {
759 denom_imag = mr * wr * wid( mr, mx2, m1_2, m2_2, r, l ); // real-i*imag;
760 }
761 else { denom_imag = mr * wr; }
762
763 double denom = denom_real * denom_real + denom_imag * denom_imag;
764 double output_x = denom_real / denom;
765 double output_y = denom_imag / denom;
766
767 complex<double> output( output_x, output_y );
768 return output;
769}
770
771// PiPi-S wave K-Matrix
772double EvtD0Topipipi0::rho22( double sc ) {
773 double rho[689] = {
774 3.70024e-18, 8.52763e-15, 1.87159e-13, 1.3311e-12, 5.61842e-12, 1.75224e-11,
775 4.48597e-11, 9.99162e-11, 2.00641e-10, 3.71995e-10, 6.47093e-10, 1.06886e-09,
776 1.69124e-09, 2.58031e-09, 3.8168e-09, 5.49601e-09, 7.72996e-09, 1.06509e-08,
777 1.44078e-08, 1.91741e-08, 2.51445e-08, 3.25345e-08, 4.15946e-08, 5.25949e-08,
778 6.58316e-08, 8.16443e-08, 1.00389e-07, 1.22455e-07, 1.48291e-07, 1.78348e-07,
779 2.1313e-07, 2.53192e-07, 2.99086e-07, 3.51462e-07, 4.10993e-07, 4.78349e-07,
780 5.54327e-07, 6.3972e-07, 7.35316e-07, 8.42099e-07, 9.61004e-07, 1.09295e-06,
781 1.2391e-06, 1.40051e-06, 1.57824e-06, 1.77367e-06, 1.98805e-06, 2.22257e-06,
782 2.47877e-06, 2.7581e-06, 3.06186e-06, 3.39182e-06, 3.74971e-06, 4.137e-06,
783 4.5555e-06, 5.00725e-06, 5.4939e-06, 6.01725e-06, 6.57992e-06, 7.18371e-06,
784 7.83044e-06, 8.52301e-06, 9.26342e-06, 1.00535e-05, 1.08967e-05, 1.17953e-05,
785 1.27514e-05, 1.37679e-05, 1.48482e-05, 1.59943e-05, 1.72088e-05, 1.84961e-05,
786 1.98586e-05, 2.12987e-05, 2.28207e-05, 2.44279e-05, 2.61228e-05, 2.79084e-05,
787 2.97906e-05, 3.17718e-05, 3.38544e-05, 3.60443e-05, 3.8345e-05, 4.07591e-05,
788 4.32903e-05, 4.59459e-05, 4.87285e-05, 5.16403e-05, 5.46887e-05, 5.7878e-05,
789 6.12111e-05, 6.46908e-05, 6.83274e-05, 7.21231e-05, 7.60817e-05, 8.0208e-05,
790 8.45102e-05, 8.89919e-05, 9.36544e-05, 9.85082e-05, 0.000103559, 0.000108812,
791 0.000114267, 0.000119938, 0.000125827, 0.00013194, 0.000138278, 0.000144857,
792 0.000151681, 0.000158752, 0.000166074, 0.000173663, 0.000181521, 0.000189652,
793 0.000198059, 0.000206761, 0.000215761, 0.000225063, 0.00023467, 0.000244599,
794 0.000254855, 0.00026544, 0.000276357, 0.000287629, 0.00029926, 0.000311253,
795 0.000323609, 0.000336351, 0.000349483, 0.000363009, 0.000376926, 0.000391264,
796 0.000406029, 0.000421225, 0.000436848, 0.000452921, 0.000469458, 0.000486461,
797 0.00050393, 0.00052187, 0.000540322, 0.000559278, 0.000578746, 0.00059872,
798 0.000619236, 0.0006403, 0.000661911, 0.000684074, 0.000706799, 0.000730127,
799 0.00075405, 0.000778569, 0.000803686, 0.000829443, 0.000855839, 0.000882879,
800 0.000910561, 0.000938898, 0.000967939, 0.000997674, 0.00102811, 0.00105923,
801 0.0010911, 0.0011237, 0.00115706, 0.00119117, 0.00122601, 0.00126168,
802 0.00129815, 0.00133543, 0.00137351, 0.00141242, 0.00145219, 0.00149283,
803 0.00153434, 0.0015767, 0.00161995, 0.00166415, 0.00170928, 0.00175534,
804 0.00180232, 0.00185028, 0.00189924, 0.00194919, 0.00200014, 0.00205207,
805 0.00210503, 0.0021591, 0.00221421, 0.0022704, 0.00232766, 0.00238602,
806 0.00244554, 0.00250619, 0.00256799, 0.0026309, 0.002695, 0.00276033,
807 0.00282689, 0.00289467, 0.00296367, 0.00303389, 0.00310543, 0.0031783,
808 0.00325244, 0.0033279, 0.0034046, 0.00348275, 0.00356229, 0.00364322,
809 0.00372555, 0.00380924, 0.00389438, 0.00398104, 0.00406914, 0.00415877,
810 0.00424985, 0.00434235, 0.00443651, 0.00453224, 0.00462954, 0.00472848,
811 0.00482894, 0.00493102, 0.00503483, 0.00514029, 0.00524749, 0.0053563,
812 0.00546675, 0.00557905, 0.0056931, 0.00580901, 0.0059267, 0.00604613,
813 0.00616735, 0.00629049, 0.00641557, 0.00654254, 0.00667142, 0.00680216,
814 0.00693472, 0.00706946, 0.00720621, 0.00734497, 0.0074858, 0.00762855,
815 0.00777338, 0.00792036, 0.00806957, 0.00822087, 0.00837426, 0.00852982,
816 0.0086875, 0.00884756, 0.00900991, 0.00917447, 0.00934137, 0.00951052,
817 0.00968194, 0.0098558, 0.010032, 0.0102108, 0.0103919, 0.0105754,
818 0.0107612, 0.0109496, 0.0111406, 0.0113343, 0.0115305, 0.0117293,
819 0.0119303, 0.0121343, 0.0123409, 0.0125502, 0.0127623, 0.0129771,
820 0.0131944, 0.0134145, 0.0136376, 0.0138636, 0.0140924, 0.0143241,
821 0.0145587, 0.0147959, 0.0150363, 0.0152797, 0.0155262, 0.0157758,
822 0.0160283, 0.0162838, 0.0165421, 0.016804, 0.0170691, 0.0173374,
823 0.0176087, 0.0178835, 0.0181612, 0.0184423, 0.0187269, 0.0190149,
824 0.0193063, 0.0196009, 0.0198991, 0.0202003, 0.0205052, 0.0208137,
825 0.0211259, 0.0214418, 0.0217611, 0.0220841, 0.0224105, 0.0227406,
826 0.0230746, 0.0234125, 0.0237542, 0.0240996, 0.0244486, 0.0248012,
827 0.025158, 0.0255188, 0.0258837, 0.0262527, 0.0266256, 0.0270025,
828 0.0273833, 0.027768, 0.0281572, 0.0285505, 0.0289483, 0.0293503,
829 0.0297564, 0.0301665, 0.0305808, 0.0309997, 0.0314231, 0.0318511,
830 0.0322835, 0.0327205, 0.0331616, 0.0336073, 0.0340576, 0.0345128,
831 0.0349727, 0.0354373, 0.0359066, 0.0363807, 0.0368589, 0.0373419,
832 0.0378302, 0.0383234, 0.0388218, 0.0393252, 0.0398336, 0.040347,
833 0.0408652, 0.041388, 0.0419165, 0.0424502, 0.0429893, 0.0435338,
834 0.0440833, 0.044638, 0.0451976, 0.0457627, 0.0463338, 0.0469103,
835 0.047492, 0.0480797, 0.0486729, 0.0492716, 0.0498757, 0.0504852,
836 0.0511009, 0.0517229, 0.0523503, 0.0529838, 0.0536231, 0.0542678,
837 0.054918, 0.0555743, 0.0562372, 0.0569065, 0.0575818, 0.0582634,
838 0.0589511, 0.0596454, 0.0603451, 0.061051, 0.0617635, 0.0624826,
839 0.0632084, 0.0639409, 0.06468, 0.0654254, 0.0661772, 0.0669346,
840 0.0676994, 0.0684714, 0.0692503, 0.0700354, 0.0708285, 0.0716277,
841 0.0724347, 0.0732479, 0.0740671, 0.0748947, 0.0757299, 0.0765715,
842 0.0774207, 0.0782771, 0.0791407, 0.0800119, 0.0808897, 0.0817743,
843 0.0826672, 0.0835684, 0.0844769, 0.0853938, 0.0863179, 0.0872493,
844 0.0881882, 0.0891349, 0.090089, 0.0910523, 0.0920236, 0.093002,
845 0.0939894, 0.094985, 0.0959887, 0.0970003, 0.0980191, 0.0990454,
846 0.100081, 0.101126, 0.10218, 0.103242, 0.104312, 0.105392,
847 0.10648, 0.107576, 0.10868, 0.109793, 0.110916, 0.112048,
848 0.113188, 0.114339, 0.115498, 0.116666, 0.117843, 0.119028,
849 0.120223, 0.121427, 0.122641, 0.123865, 0.125098, 0.126342,
850 0.127595, 0.128857, 0.130128, 0.131409, 0.132701, 0.134002,
851 0.135314, 0.136635, 0.137966, 0.139308, 0.14066, 0.142022,
852 0.143394, 0.144774, 0.146166, 0.14757, 0.148985, 0.15041,
853 0.151845, 0.153291, 0.154749, 0.156215, 0.157694, 0.159182,
854 0.160682, 0.162194, 0.163718, 0.165251, 0.166797, 0.168354,
855 0.169921, 0.1715, 0.17309, 0.17469, 0.176304, 0.177929,
856 0.179566, 0.181216, 0.182878, 0.184553, 0.186238, 0.187934,
857 0.189642, 0.191362, 0.193096, 0.194842, 0.196602, 0.198374,
858 0.200158, 0.201954, 0.203764, 0.205586, 0.207421, 0.209266,
859 0.211124, 0.212997, 0.214882, 0.216783, 0.218697, 0.220624,
860 0.222565, 0.224518, 0.226486, 0.228466, 0.230458, 0.232463,
861 0.234484, 0.23652, 0.238569, 0.240633, 0.242711, 0.244803,
862 0.246909, 0.249031, 0.251165, 0.253313, 0.255475, 0.257649,
863 0.259841, 0.262051, 0.264274, 0.266514, 0.268768, 0.271036,
864 0.273319, 0.275618, 0.277932, 0.280259, 0.282602, 0.28496,
865 0.287338, 0.28973, 0.292138, 0.294563, 0.297003, 0.299458,
866 0.30193, 0.304417, 0.306919, 0.309437, 0.311972, 0.314526,
867 0.317095, 0.319684, 0.322289, 0.324911, 0.327551, 0.330205,
868 0.332876, 0.335567, 0.338271, 0.340993, 0.343736, 0.346496,
869 0.349272, 0.352065, 0.354878, 0.35771, 0.360561, 0.363426,
870 0.366311, 0.369212, 0.372128, 0.375067, 0.378027, 0.381006,
871 0.384001, 0.387014, 0.39005, 0.393106, 0.396181, 0.399271,
872 0.402384, 0.405513, 0.408661, 0.41183, 0.41502, 0.418233,
873 0.421462, 0.424709, 0.42798, 0.43127, 0.434583, 0.437914,
874 0.441267, 0.444637, 0.448022, 0.451434, 0.454868, 0.458328,
875 0.461805, 0.465302, 0.468821, 0.472364, 0.475928, 0.47951,
876 0.483119, 0.486748, 0.490397, 0.494066, 0.497758, 0.501477,
877 0.505217, 0.508977, 0.512762, 0.516567, 0.520394, 0.524247,
878 0.528125, 0.532027, 0.535947, 0.53989, 0.543852, 0.547844,
879 0.551863, 0.555904, 0.559966, 0.56406, 0.568177, 0.572312,
880 0.576471, 0.580662, 0.584875, 0.58911, 0.593373, 0.597653,
881 0.601965, 0.606301, 0.610663, 0.615051, 0.619465, 0.623907,
882 0.62837, 0.632863, 0.637383, 0.641924, 0.646494, 0.651091,
883 0.655708, 0.660356, 0.665027, 0.669732, 0.674464, 0.679227,
884 0.684016, 0.688827, 0.693664, 0.698532, 0.703428, 0.708353,
885 0.713307, 0.718283, 0.72329, 0.728322, 0.733387, 0.738479,
886 0.743605, 0.748763, 0.753949, 0.759163, 0.764407, 0.769674,
887 0.774973, 0.780311, 0.78567, 0.791057, 0.796476, 0.801922,
888 0.8074, 0.812919, 0.818466, 0.824044 };
889
890 double m2 = 0.13957 * 0.13957;
891 double smin = ( 0.13957 * 4 ) * ( 0.13957 * 4 );
892 double dh = 0.001;
893 int od = ( sc - 0.312 ) / dh;
894 double sc_m = 0.312 + od * dh;
895 double rhouse = 0;
896 if ( sc >= 0.312 && sc < 1 )
897 { rhouse = ( ( sc - sc_m ) / dh ) * ( rho[od + 1] - rho[od] ) + rho[od]; }
898 else if ( sc < 0.312 && sc >= smin )
899 { rhouse = ( ( sc - smin ) / ( 0.312 - smin ) ) * rho[0]; }
900 else if ( sc >= 1 )
901 {
902 // rhouse = (1-16*m2/sc)*(1-16*m2/sc)*(1-16*m2/sc)*(1-16*m2/sc)*(1-16*m2/sc);
903 rhouse = sqrt( 1 - 16 * m2 / sc );
904 }
905 else { rhouse = 0; }
906 return rhouse;
907}
908
909complex<double> EvtD0Topipipi0::rhoMTX( int i, int j, double s ) {
910
911 double rhoijx;
912 double rhoijy;
913 double mpi = 0.13957;
914 if ( i == j && i == 0 )
915 {
916 double m2 = 0.13957 * 0.13957;
917 if ( ( 1 - ( 4 * m2 ) / s ) > 0 )
918 {
919 rhoijx = sqrt( 1.0f - ( 4 * m2 ) / s );
920 rhoijy = 0;
921 }
922 else
923 {
924 rhoijy = sqrt( ( 4 * m2 ) / s - 1.0f );
925 rhoijx = 0;
926 }
927 }
928 if ( i == j && i == 1 )
929 {
930 double m2 = 0.493677 * 0.493677;
931 if ( ( 1 - ( 4 * m2 ) / s ) > 0 )
932 {
933 rhoijx = sqrt( 1.0f - ( 4 * m2 ) / s );
934 rhoijy = 0;
935 }
936 else
937 {
938 rhoijy = sqrt( ( 4 * m2 ) / s - 1.0f );
939 rhoijx = 0;
940 }
941 }
942 if ( i == j && i == 2 )
943 {
944 rhoijx = rho22( s );
945 rhoijy = 0;
946 }
947 if ( i == j && i == 3 )
948 {
949 double m2 = 0.547862 * 0.547862;
950 if ( ( 1 - ( 4 * m2 ) / s ) > 0 )
951 {
952 rhoijx = sqrt( 1.0f - ( 4 * m2 ) / s );
953 rhoijy = 0;
954 }
955 else
956 {
957 rhoijy = sqrt( ( 4 * m2 ) / s - 1.0f );
958 rhoijx = 0;
959 }
960 }
961 if ( i == j && i == 4 )
962 {
963 double m_1 = 0.547862;
964 double m_2 = 0.95778;
965 double mp2 = ( m_1 + m_2 ) * ( m_1 + m_2 );
966 double mm2 = ( m_1 - m_2 ) * ( m_1 - m_2 );
967 if ( ( 1 - mp2 / s ) > 0 )
968 {
969 rhoijx = sqrt( 1.0f - mp2 / s );
970 rhoijy = 0;
971 }
972 else
973 {
974 rhoijy = sqrt( mp2 / s - 1.0f );
975 rhoijx = 0;
976 }
977 }
978
979 if ( i != j )
980 {
981 rhoijx = 0;
982 rhoijy = 0;
983 }
984 complex<double> rhoij( rhoijx, rhoijy );
985 return rhoij;
986}
987
988/* Kij = (sum_a gigj/(ma^2-s) + fij(1-s0)/(s-s0))((s-0.5sampi2)(1-sa0)/(s-sa0)) */
989complex<double> EvtD0Topipipi0::KMTX( int i, int j, double s ) {
990
991 double Kijx;
992 double Kijy;
993 double mpi = 0.13957;
994 double m[5] = { 0.65100, 1.20360, 1.55817, 1.21000, 1.82206 };
995
996 double g1[5] = { 0.22889, -0.55377, 0.00000, -0.39899, -0.34639 };
997 double g2[5] = { 0.94128, 0.55095, 0.00000, 0.39065, 0.31503 };
998 double g3[5] = { 0.36856, 0.23888, 0.55639, 0.18340, 0.18681 };
999 double g4[5] = { 0.33650, 0.40907, 0.85679, 0.19906, -0.00984 };
1000 double g5[5] = { 0.18171, -0.17558, -0.79658, -0.00355, 0.22358 };
1001
1002 double f1[5] = { 0.23399, 0.15044, -0.20545, 0.32825, 0.35412 };
1003
1004 double eps = 1e-11;
1005
1006 double down[5] = { 0, 0, 0, 0, 0 };
1007 double upreal[5] = { 0, 0, 0, 0, 0 };
1008 double upimag[5] = { 0, 0, 0, 0, 0 };
1009
1010 for ( int k = 0; k < 5; k++ )
1011 {
1012
1013 /* down[k] = (m[k]*m[k]-s)*(m[k]*m[k]-s)+eps*eps;
1014 upreal[k] = (m[k]*m[k]-s)/down[k];
1015 upimag[k] = -1.0f * eps/down[k]; */
1016
1017 double dm2 = m[k] * m[k] - s;
1018 if ( fabs( dm2 ) < eps && dm2 <= 0 ) dm2 = -eps;
1019 if ( fabs( dm2 ) < eps && dm2 > 0 ) dm2 = eps;
1020 upreal[k] = 1.0f / dm2;
1021 upimag[k] = 0;
1022 }
1023
1024 double tmp1x = g1[i] * g1[j] * upreal[0] + g2[i] * g2[j] * upreal[1] +
1025 g3[i] * g3[j] * upreal[2] + g4[i] * g4[j] * upreal[3] +
1026 g5[i] * g5[j] * upreal[4];
1027 double tmp1y = g1[i] * g1[j] * upimag[0] + g2[i] * g2[j] * upimag[1] +
1028 g3[i] * g3[j] * upimag[2] + g4[i] * g4[j] * upimag[3] +
1029 g5[i] * g5[j] * upimag[4];
1030
1031 double tmp2 = 0;
1032 if ( i == 0 ) { tmp2 = f1[j] * ( 1 + 3.92637 ) / ( s + 3.92637 ); }
1033 if ( j == 0 ) { tmp2 = f1[i] * ( 1 + 3.92637 ) / ( s + 3.92637 ); }
1034 double tmp3 = ( s - 0.5 * mpi * mpi ) * ( 1 + 0.15 ) / ( s + 0.15 );
1035
1036 Kijx = ( tmp1x + tmp2 ) * tmp3;
1037 Kijy = (tmp1y)*tmp3;
1038
1039 complex<double> Kij( Kijx, Kijy );
1040 return Kij;
1041}
1042
1043complex<double> EvtD0Topipipi0::IMTX( int i, int j ) {
1044
1045 double Iijx;
1046 double Iijy;
1047 if ( i == j )
1048 {
1049 Iijx = 1;
1050 Iijy = 0;
1051 }
1052 else
1053 {
1054 Iijx = 0;
1055 Iijy = 0;
1056 }
1057 complex<double> Iij( Iijx, Iijy );
1058 return Iij;
1059}
1060
1061/* F = I - i*K*rho */
1062complex<double> EvtD0Topipipi0::FMTX( double Kijx, double Kijy, double rhojjx, double rhojjy,
1063 int i, int j ) {
1064
1065 double Fijx;
1066 double Fijy;
1067
1068 double tmpx = rhojjx * Kijx - rhojjy * Kijy;
1069 double tmpy = rhojjx * Kijy + rhojjy * Kijx;
1070
1071 Fijx = IMTX( i, j ).real() + tmpy;
1072 Fijy = -tmpx;
1073
1074 complex<double> Fij( Fijx, Fijy );
1075 return Fij;
1076}
1077
1078/* inverse for Matrix (I - i*rho*K) using LUP */
1079double EvtD0Topipipi0::FINVMTX( double s, double* FINVx, double* FINVy ) {
1080
1081 int P[5] = { 0, 1, 2, 3, 4 };
1082
1083 double Fx[5][5];
1084 double Fy[5][5];
1085
1086 double Ux[5][5];
1087 double Uy[5][5];
1088 double Lx[5][5];
1089 double Ly[5][5];
1090
1091 double UIx[5][5];
1092 double UIy[5][5];
1093 double LIx[5][5];
1094 double LIy[5][5];
1095
1096 for ( int k = 0; k < 5; k++ )
1097 {
1098 double rhokkx = rhoMTX( k, k, s ).real();
1099 double rhokky = rhoMTX( k, k, s ).imag();
1100 Ux[k][k] = rhokkx;
1101 Uy[k][k] = rhokky;
1102 for ( int l = k; l < 5; l++ )
1103 {
1104 double Kklx = KMTX( k, l, s ).real();
1105 double Kkly = KMTX( k, l, s ).imag();
1106 Lx[k][l] = Kklx;
1107 Ly[k][l] = Kkly;
1108 Lx[l][k] = Lx[k][l];
1109 Ly[l][k] = Ly[k][l];
1110 }
1111 }
1112
1113 for ( int k = 0; k < 5; k++ )
1114 {
1115 for ( int l = 0; l < 5; l++ )
1116 {
1117 double Fklx = FMTX( Lx[k][l], Ly[k][l], Ux[l][l], Uy[l][l], k, l ).real();
1118 double Fkly = FMTX( Lx[k][l], Ly[k][l], Ux[l][l], Uy[l][l], k, l ).imag();
1119 Fx[k][l] = Fklx;
1120 Fy[k][l] = Fkly;
1121 }
1122 }
1123
1124 for ( int k = 0; k < 5; k++ )
1125 {
1126 double tmprM = ( Fx[k][k] * Fx[k][k] + Fy[k][k] * Fy[k][k] );
1127 int tmpID = 0;
1128 for ( int l = k; l < 5; l++ )
1129 {
1130 double tmprF = ( Fx[l][k] * Fx[l][k] + Fy[l][k] * Fy[l][k] );
1131 if ( tmprM <= tmprF )
1132 {
1133 tmprM = tmprF;
1134 tmpID = l;
1135 }
1136 }
1137
1138 int tmpP = P[k];
1139 P[k] = P[tmpID];
1140 P[tmpID] = tmpP;
1141
1142 for ( int l = 0; l < 5; l++ )
1143 {
1144
1145 double tmpFx = Fx[k][l];
1146 double tmpFy = Fy[k][l];
1147
1148 Fx[k][l] = Fx[tmpID][l];
1149 Fy[k][l] = Fy[tmpID][l];
1150
1151 Fx[tmpID][l] = tmpFx;
1152 Fy[tmpID][l] = tmpFy;
1153 }
1154
1155 for ( int l = k + 1; l < 5; l++ )
1156 {
1157 double rFkk = Fx[k][k] * Fx[k][k] + Fy[k][k] * Fy[k][k];
1158 double Fxlk = Fx[l][k];
1159 double Fylk = Fy[l][k];
1160 double Fxkk = Fx[k][k];
1161 double Fykk = Fy[k][k];
1162 Fx[l][k] = ( Fxlk * Fxkk + Fylk * Fykk ) / rFkk;
1163 Fy[l][k] = ( Fylk * Fxkk - Fxlk * Fykk ) / rFkk;
1164 for ( int m = k + 1; m < 5; m++ )
1165 {
1166 Fx[l][m] = Fx[l][m] - ( Fx[l][k] * Fx[k][m] - Fy[l][k] * Fy[k][m] );
1167 Fy[l][m] = Fy[l][m] - ( Fx[l][k] * Fy[k][m] + Fy[l][k] * Fx[k][m] );
1168 }
1169 }
1170 }
1171
1172 for ( int k = 0; k < 5; k++ )
1173 {
1174 for ( int l = 0; l < 5; l++ )
1175 {
1176 if ( k == l )
1177 {
1178 Lx[k][k] = 1;
1179 Ly[k][k] = 0;
1180 Ux[k][k] = Fx[k][k];
1181 Uy[k][k] = Fy[k][k];
1182 }
1183 if ( k > l )
1184 {
1185 Lx[k][l] = Fx[k][l];
1186 Ly[k][l] = Fy[k][l];
1187 Ux[k][l] = 0;
1188 Uy[k][l] = 0;
1189 }
1190 if ( k < l )
1191 {
1192 Ux[k][l] = Fx[k][l];
1193 Uy[k][l] = Fy[k][l];
1194 Lx[k][l] = 0;
1195 Ly[k][l] = 0;
1196 }
1197 }
1198 }
1199
1200 // calculate Inverse for L and U
1201 for ( int k = 0; k < 5; k++ )
1202 {
1203
1204 LIx[k][k] = 1;
1205 LIy[k][k] = 0;
1206
1207 double rUkk = Ux[k][k] * Ux[k][k] + Uy[k][k] * Uy[k][k];
1208 UIx[k][k] = Ux[k][k] / rUkk;
1209 UIy[k][k] = -1.0f * Uy[k][k] / rUkk;
1210
1211 for ( int l = ( k + 1 ); l < 5; l++ )
1212 {
1213 LIx[k][l] = 0;
1214 LIy[k][l] = 0;
1215 UIx[l][k] = 0;
1216 UIy[l][k] = 0;
1217 }
1218
1219 for ( int l = ( k - 1 ); l >= 0; l-- )
1220 { // U-1
1221 double sx = 0;
1222 double c_sx = 0;
1223 double sy = 0;
1224 double c_sy = 0;
1225 for ( int m = l + 1; m <= k; m++ )
1226 {
1227 sx = sx - c_sx;
1228 double sx_tmp = sx + Ux[l][m] * UIx[m][k] - Uy[l][m] * UIy[m][k];
1229 c_sx = ( sx_tmp - sx ) - ( Ux[l][m] * UIx[m][k] - Uy[l][m] * UIy[m][k] );
1230 sx = sx_tmp;
1231
1232 sy = sy - c_sy;
1233 double sy_tmp = sy + Ux[l][m] * UIy[m][k] + Uy[l][m] * UIx[m][k];
1234 c_sy = ( sy_tmp - sy ) - ( Ux[l][m] * UIy[m][k] + Uy[l][m] * UIx[m][k] );
1235 sy = sy_tmp;
1236 }
1237 UIx[l][k] = -1.0f * ( UIx[l][l] * sx - UIy[l][l] * sy );
1238 UIy[l][k] = -1.0f * ( UIy[l][l] * sx + UIx[l][l] * sy );
1239 }
1240
1241 for ( int l = k + 1; l < 5; l++ )
1242 { // L-1
1243 double sx = 0;
1244 double c_sx = 0;
1245 double sy = 0;
1246 double c_sy = 0;
1247 for ( int m = k; m < l; m++ )
1248 {
1249 sx = sx - c_sx;
1250 double sx_tmp = sx + Lx[l][m] * LIx[m][k] - Ly[l][m] * LIy[m][k];
1251 c_sx = ( sx_tmp - sx ) - ( Lx[l][m] * LIx[m][k] - Ly[l][m] * LIy[m][k] );
1252 sx = sx_tmp;
1253
1254 sy = sy - c_sy;
1255 double sy_tmp = sy + Lx[l][m] * LIy[m][k] + Ly[l][m] * LIx[m][k];
1256 c_sy = ( sy_tmp - sy ) - ( Lx[l][m] * LIy[m][k] + Ly[l][m] * LIx[m][k] );
1257 sy = sy_tmp;
1258 }
1259 LIx[l][k] = -1.0f * sx;
1260 LIy[l][k] = -1.0f * sy;
1261 }
1262 }
1263
1264 for ( int m = 0; m < 5; m++ )
1265 {
1266 double resX = 0;
1267 double c_resX = 0;
1268 double resY = 0;
1269 double c_resY = 0;
1270 for ( int k = 0; k < 5; k++ )
1271 {
1272 for ( int l = 0; l < 5; l++ )
1273 {
1274 double Plm = 0;
1275 if ( P[l] == m ) Plm = 1;
1276
1277 resX = resX - c_resX;
1278 double resX_tmp = resX + ( UIx[0][k] * LIx[k][l] - UIy[0][k] * LIy[k][l] ) * Plm;
1279 c_resX =
1280 ( resX_tmp - resX ) - ( ( UIx[0][k] * LIx[k][l] - UIy[0][k] * LIy[k][l] ) * Plm );
1281 resX = resX_tmp;
1282
1283 resY = resY - c_resY;
1284 double resY_tmp = resY + ( UIx[0][k] * LIy[k][l] + UIy[0][k] * LIx[k][l] ) * Plm;
1285 c_resY =
1286 ( resY_tmp - resY ) - ( ( UIx[0][k] * LIy[k][l] + UIy[0][k] * LIx[k][l] ) * Plm );
1287 resY = resY_tmp;
1288 }
1289 }
1290 FINVx[m] = resX;
1291 FINVy[m] = resY;
1292 }
1293
1294 return 1.0f;
1295}
1296
1297complex<double> EvtD0Topipipi0::PVTR( int ID, double s ) {
1298
1299 double VPix;
1300 double VPiy;
1301 double m[5] = { 0.65100, 1.20360, 1.55817, 1.21000, 1.82206 };
1302
1303 double eps = 1e-11;
1304
1305 /* double down = (m[ID]*m[ID]-s)*(m[ID]*m[ID]-s)+eps*eps;
1306 double upreal = (m[ID]*m[ID]-s)/down;
1307 double upimag = -1.0f * eps/down; */
1308
1309 double dm2 = m[ID] * m[ID] - s;
1310
1311 if ( fabs( dm2 ) < eps && dm2 <= 0 ) dm2 = -eps;
1312 if ( fabs( dm2 ) < eps && dm2 > 0 ) dm2 = eps;
1313
1314 VPix = 1.0f / dm2;
1315 VPiy = 0;
1316
1317 complex<double> VPi( VPix, VPiy );
1318 return VPi;
1319}
1320
1321complex<double> EvtD0Topipipi0::Fvector( double sa, double s0, int l ) {
1322
1323 double outputx = 0;
1324 double outputy = 0;
1325
1326 double FINVx[5] = { 0, 0, 0, 0, 0 };
1327 double FINVy[5] = { 0, 0, 0, 0, 0 };
1328
1329 double tmpFLAG = FINVMTX( sa, FINVx, FINVy );
1330
1331 if ( l < 5 )
1332 {
1333 double g[5][5] = { { 0.22889, -0.55377, 0.00000, -0.39899, -0.34639 },
1334 { 0.94128, 0.55095, 0.00000, 0.39065, 0.31503 },
1335 { 0.36856, 0.23888, 0.55639, 0.18340, 0.18681 },
1336 { 0.33650, 0.40907, 0.85679, 0.19906, -0.00984 },
1337 { 0.18171, -0.17558, -0.79658, -0.00355, 0.22358 } };
1338 double resx = 0;
1339 double c_resx = 0;
1340 double resy = 0;
1341 double c_resy = 0;
1342 double Plx = PVTR( l, sa ).real();
1343 double Ply = PVTR( l, sa ).imag();
1344 for ( int j = 0; j < 5; j++ )
1345 {
1346 resx = resx - c_resx;
1347 double resx_tmp = resx + ( FINVx[j] * g[l][j] * Plx - FINVy[j] * g[l][j] * Ply );
1348 c_resx = ( resx_tmp - resx ) - ( FINVx[j] * g[l][j] * Plx - FINVy[j] * g[l][j] * Ply );
1349 resx = resx_tmp;
1350
1351 resy = resy - c_resy;
1352 double resy_tmp = resy + ( FINVx[j] * g[l][j] * Ply + FINVy[j] * g[l][j] * Plx );
1353 c_resy = ( resy_tmp - resy ) - ( FINVx[j] * g[l][j] * Ply + FINVy[j] * g[l][j] * Plx );
1354 resy = resy_tmp;
1355 }
1356 outputx = resx;
1357 outputy = resy;
1358 }
1359 else
1360 {
1361 int idx = l - 5;
1362 double eps = 1e-11;
1363 double ds = sa - s0;
1364 if ( fabs( ds ) < eps && ds <= 0 ) ds = -eps;
1365 if ( fabs( ds ) < eps && ds > 0 ) ds = eps;
1366 double tmp = ( 1 - s0 ) / ds;
1367 outputx = FINVx[idx] * tmp;
1368 outputy = FINVy[idx] * tmp;
1369 }
1370
1371 complex<double> output( outputx, outputy );
1372 return output;
1373}
1374
1375complex<double> EvtD0Topipipi0::CalD0Amp() { return Amp( p4_Pip, p4_Pim, p4_Pi0 ); }
1376complex<double> EvtD0Topipipi0::CalDbAmp() {
1377
1378 vector<double> cpPip;
1379 cpPip.clear();
1380 vector<double> cpPim;
1381 cpPim.clear();
1382 vector<double> cpPi0;
1383 cpPi0.clear();
1384
1385 cpPip.push_back( -p4_Pim[0] );
1386 cpPim.push_back( -p4_Pip[0] );
1387 cpPi0.push_back( -p4_Pi0[0] );
1388 cpPip.push_back( -p4_Pim[1] );
1389 cpPim.push_back( -p4_Pip[1] );
1390 cpPi0.push_back( -p4_Pi0[1] );
1391 cpPip.push_back( -p4_Pim[2] );
1392 cpPim.push_back( -p4_Pip[2] );
1393 cpPi0.push_back( -p4_Pi0[2] );
1394 cpPip.push_back( p4_Pim[3] );
1395 cpPim.push_back( p4_Pip[3] );
1396 cpPi0.push_back( p4_Pi0[3] );
1397
1398 return ( -1.0 ) * Amp( cpPip, cpPim, cpPi0 );
1399}
1400
1401complex<double> EvtD0Topipipi0::Amp( vector<double> Pip, vector<double> Pim,
1402 vector<double> Pi0 ) {
1403
1404 if ( fitpara.size() != 14 )
1405 {
1406 cout << "Error!!! number of para: " << fitpara.size() << endl;
1407 exit( 0 );
1408 }
1409
1410 vector<double> PipPim;
1411 PipPim.clear();
1412 vector<double> PipPi0;
1413 PipPi0.clear();
1414 vector<double> PimPi0;
1415 PimPi0.clear();
1416
1417 PipPim = sum_tensor( Pip, Pim );
1418 PipPi0 = sum_tensor( Pip, Pi0 );
1419 PimPi0 = sum_tensor( Pim, Pi0 );
1420
1421 vector<double> D0;
1422 D0.clear();
1423 D0 = sum_tensor( PipPim, Pi0 );
1424
1425 double M2_PipPim = contract_11_0( PipPim, PipPim );
1426 double M2_PipPi0 = contract_11_0( PipPi0, PipPi0 );
1427 double M2_PimPi0 = contract_11_0( PimPi0, PimPi0 );
1428
1429 double M2_D0 = contract_11_0( D0, D0 );
1430
1431 complex<double> GS_rho770_z = GS( M2_PipPim, m0_rho7700, w0_rho7700, m2_Pi, m2_Pi, rRes, 1 );
1432 complex<double> GS_rho770_p =
1433 GS( M2_PipPi0, m0_rho770p, w0_rho770p, m2_Pi, m2_Pi0, rRes, 1 );
1434 complex<double> GS_rho770_m =
1435 GS( M2_PimPi0, m0_rho770p, w0_rho770p, m2_Pi, m2_Pi0, rRes, 1 );
1436
1437 complex<double> PiPiS_pm_0 = Fvector( M2_PipPim, s0_prod, 0 );
1438 complex<double> PiPiS_pm_1 = Fvector( M2_PipPim, s0_prod, 1 );
1439 complex<double> PiPiS_pm_2 = Fvector( M2_PipPim, s0_prod, 2 );
1440 complex<double> PiPiS_pm_3 = Fvector( M2_PipPim, s0_prod, 3 );
1441 complex<double> PiPiS_pm_4 = Fvector( M2_PipPim, s0_prod, 4 );
1442 complex<double> PiPiS_pm_5 = Fvector( M2_PipPim, s0_prod, 5 );
1443 complex<double> PiPiS_pm_6 = Fvector( M2_PipPim, s0_prod, 6 );
1444 complex<double> PiPiS_pm_7 = Fvector( M2_PipPim, s0_prod, 7 );
1445 complex<double> PiPiS_pm_8 = Fvector( M2_PipPim, s0_prod, 8 );
1446 complex<double> PiPiS_pm_9 = Fvector( M2_PipPim, s0_prod, 9 );
1447
1448 complex<double> RBW_f21270 = RBW( M2_PipPim, m0_f21270, w0_f21270, m2_Pi, m2_Pi, rRes, 2 );
1449
1450 // X->PP Orbital
1451 vector<double> T1_PipPim;
1452 T1_PipPim.clear();
1453 vector<double> T1_PipPi0;
1454 T1_PipPi0.clear();
1455 vector<double> T1_PimPi0;
1456 T1_PimPi0.clear();
1457
1458 T1_PipPim = OrbitalTensors( PipPim, Pip, Pim, rRes, 1 );
1459 T1_PipPi0 = OrbitalTensors( PipPi0, Pip, Pi0, rRes, 1 );
1460 T1_PimPi0 = OrbitalTensors( PimPi0, Pim, Pi0, rRes, 1 );
1461
1462 vector<double> T2_PipPim;
1463 T2_PipPim.clear();
1464
1465 T2_PipPim = OrbitalTensors( PipPim, Pip, Pim, rRes, 2 );
1466
1467 // D->XP Orbital
1468 vector<double> T1_PipPimPi0;
1469 T1_PipPimPi0.clear();
1470 vector<double> T1_PipPi0Pim;
1471 T1_PipPi0Pim.clear();
1472 vector<double> T1_PimPi0Pip;
1473 T1_PimPi0Pip.clear();
1474
1475 T1_PipPimPi0 = OrbitalTensors( D0, PipPim, Pi0, rD, 1 );
1476 T1_PipPi0Pim = OrbitalTensors( D0, PipPi0, Pim, rD, 1 );
1477 T1_PimPi0Pip = OrbitalTensors( D0, PimPi0, Pip, rD, 1 );
1478
1479 vector<double> T2_PipPimPi0;
1480 T2_PipPimPi0.clear();
1481
1482 T2_PipPimPi0 = OrbitalTensors( D0, PipPim, Pi0, rD, 2 );
1483
1484 complex<double> amplitude( 0, 0 );
1485
1486 // D-> VP
1487 double SF_VpPm = contract_11_0( T1_PipPi0Pim, T1_PipPi0 );
1488 amplitude += fitpara[0] * ( SF_VpPm * GS_rho770_p );
1489
1490 double SF_VmPp = contract_11_0( T1_PimPi0Pip, T1_PimPi0 );
1491 amplitude += fitpara[1] * ( SF_VmPp * GS_rho770_m );
1492
1493 double SF_VzPz = contract_11_0( T1_PipPimPi0, T1_PipPim );
1494 amplitude += fitpara[2] * ( SF_VzPz * GS_rho770_z );
1495
1496 // D-> SP
1497 amplitude += fitpara[3] * ( PiPiS_pm_0 );
1498 amplitude += fitpara[4] * ( PiPiS_pm_1 );
1499 amplitude += fitpara[5] * ( PiPiS_pm_2 );
1500 amplitude += fitpara[6] * ( PiPiS_pm_3 );
1501 amplitude += fitpara[7] * ( PiPiS_pm_4 );
1502 amplitude += fitpara[8] * ( PiPiS_pm_5 );
1503 amplitude += fitpara[9] * ( PiPiS_pm_6 );
1504 amplitude += fitpara[10] * ( PiPiS_pm_7 );
1505 amplitude += fitpara[11] * ( PiPiS_pm_8 );
1506 amplitude += fitpara[12] * ( PiPiS_pm_9 );
1507
1508 // D0 -> TP
1509 double SF_TzPz = contract_22_0( T2_PipPimPi0, T2_PipPim );
1510 amplitude += fitpara[13] * ( SF_TzPz * RBW_f21270 );
1511
1512 return amplitude;
1513}
1514
1515int EvtD0Topipipi0::CalAmp() {
1516 m_AmpD0 = CalD0Amp();
1517 m_AmpDb = CalDbAmp();
1518 return 0;
1519}
1520
1521double EvtD0Topipipi0::mag2( complex<double> x ) {
1522 double temp = x.real() * x.real() + x.imag() * x.imag();
1523 return temp;
1524}
1525
1526double EvtD0Topipipi0::mag2_itf( complex<double> x, complex<double> y ) {
1527 double temp = 2 * ( x.real() * y.real() + x.imag() * y.imag() );
1528 return temp;
1529}
1530
1531double EvtD0Topipipi0::arg( complex<double> x ) {
1532 double temp = atan( x.imag() / x.real() );
1533 if ( x.real() < 0 ) temp = temp + PI;
1534 return temp;
1535}
1536double EvtD0Topipipi0::Get_strongPhase() {
1537 double temp = arg( m_AmpD0 ) - arg( m_AmpDb );
1538 while ( temp < -PI ) { temp += 2.0 * PI; }
1539 while ( temp > PI ) { temp -= 2.0 * PI; }
1540 return temp;
1541}
1542
1543double EvtD0Topipipi0::AmplitudeSquare( int charm, int tagmode ) {
1544
1545 EvtVector4R dp1 = GetDaugMomLab( 0 ), dp2 = GetDaugMomLab( 1 ),
1546 dp3 = GetDaugMomLab( 2 ); // D0 -> pi+ pi- pi0
1547 EvtVector4R mp = dp1 + dp2 + dp3;
1548
1549 double emp = mp.get( 0 );
1550 EvtVector3R boostp3( -mp.get( 1 ) / emp, -mp.get( 2 ) / emp, -mp.get( 3 ) / emp );
1551
1552 EvtVector4R dp1bst = boostTo( dp1, boostp3 );
1553 EvtVector4R dp2bst = boostTo( dp2, boostp3 );
1554 EvtVector4R dp3bst = boostTo( dp3, boostp3 );
1555
1556 double p4pip[4], p4pim[4], p4pi0[4];
1557 for ( int i = 0; i < 3; i++ )
1558 {
1559 p4pip[i] = dp1bst.get( i + 1 );
1560 p4pim[i] = dp2bst.get( i + 1 );
1561 p4pi0[i] = dp3bst.get( i + 1 );
1562 }
1563
1564 p4pip[3] = dp1bst.get( 0 );
1565 p4pim[3] = dp2bst.get( 0 );
1566 p4pi0[3] = dp3bst.get( 0 );
1567
1568 setInput( p4pip, p4pim, p4pi0 );
1569 CalAmp();
1570
1571 complex<double> ampD0, ampDb;
1572 if ( charm > 0 )
1573 {
1574 ampD0 = Get_AmpD0();
1575 ampDb = Get_AmpDb();
1576 }
1577 else
1578 {
1579 ampD0 = Get_AmpDb();
1580 ampDb = Get_AmpD0();
1581 }
1582
1583 double ampsq = 1e-20;
1584 double r_tag = 0, R_tag = 0, delta_tag = 0;
1585
1586 if ( tagmode == 1 || tagmode == 2 || tagmode == 3 )
1587 {
1588 if ( tagmode == 1 )
1589 { // Kpi
1590 r_tag = 0.0586;
1591 R_tag = 1;
1592 delta_tag = 192.1 / 180.0 * 3.1415926;
1593 }
1594 else if ( tagmode == 2 )
1595 { // Kpipi0
1596 r_tag = 0.0441;
1597 R_tag = 0.79;
1598 delta_tag = 196.0 / 180.0 * 3.1415926;
1599 }
1600 else if ( tagmode == 3 )
1601 { // K3pi
1602 r_tag = 0.0550;
1603 R_tag = 0.44;
1604 delta_tag = 161.0 / 180.0 * 3.1415926;
1605 }
1606 complex<double> qcf( r_tag * R_tag * cos( -delta_tag ),
1607 r_tag * R_tag * sin( -delta_tag ) );
1608 complex<double> ampD0_part1 = ampD0 - qcf * ampDb;
1609 ampsq = mag2( ampD0_part1 ) + r_tag * r_tag * ( 1 - R_tag * R_tag ) * ( mag2( ampDb ) );
1610 }
1611 else { ampsq = mag2( ampD0 ); }
1612
1613 return ( ampsq <= 0 ) ? 1e-20 : ampsq;
1614}
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
PI
#define PI
Definition EvtD0ToKKpi0.cc:40
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
EvtD0Topipipi0::getName
void getName(std::string &name)
Definition EvtD0Topipipi0.cc:44
EvtD0Topipipi0::clone
EvtDecayBase * clone()
Definition EvtD0Topipipi0.cc:46
EvtD0Topipipi0::~EvtD0Topipipi0
virtual ~EvtD0Topipipi0()
Definition EvtD0Topipipi0.cc:42
EvtD0Topipipi0::decay
void decay(EvtParticle *p)
Definition EvtD0Topipipi0.cc:160
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