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