BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToKpipi.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// Environment:
3// This software is part of the EvtGen package developed jointly
4// for the BaBar and CLEO collaborations. If you use all or part
5// of it, please give an appropriate acknowledgement.
6//
7// Copyright Information: See EvtGen/COPYRIGHT
8// Copyright (C) 1998 Caltech, UCSB
9//
10// Module: EvtDsToKpipi.cc
11//
12// Description: Routine to handle three-body decays of Ds+/Ds-
13//
14// Modification history:
15// Liaoyuan Dong Sun Jun 02 01:23:25 2024 Module created
16//------------------------------------------------------------------------
17#include "EvtDsToKpipi.hh"
18#include "../EvtGenBase/EvtConst.hh"
19#include "../EvtGenBase/EvtDecayTable.hh"
20#include "../EvtGenBase/EvtGenKine.hh"
21#include "../EvtGenBase/EvtPDL.hh"
22#include "../EvtGenBase/EvtParticle.hh"
23#include "../EvtGenBase/EvtPatches.hh"
24#include "../EvtGenBase/EvtReport.hh"
25#include "../EvtGenBase/EvtResonance.hh"
26#include "../EvtGenBase/EvtResonance2.hh"
27#include <fstream>
28#include <stdlib.h>
29#include <string>
30using std::endl;
31
32EvtDsToKpipi::~EvtDsToKpipi() {}
33
34void EvtDsToKpipi::getName( std::string& model_name ) { model_name = "DsToKpipi"; }
35
36EvtDecayBase* EvtDsToKpipi::clone() { return new EvtDsToKpipi; }
37
38void EvtDsToKpipi::init() {
39 // check that there are 0 arguments
40 checkNArg( 0 );
41 checkNDaug( 3 );
42 checkSpinParent( EvtSpinType::SCALAR );
43 checkSpinDaughter( 0, EvtSpinType::SCALAR );
44 checkSpinDaughter( 1, EvtSpinType::SCALAR );
45 checkSpinDaughter( 2, EvtSpinType::SCALAR );
46
47 phi[0] = 0;
48 rho[0] = 1;
49 phi[1] = 0;
50 rho[1] = 0;
51 phi[2] = 0;
52 rho[2] = 0;
53 phi[3] = 0;
54 rho[3] = 0;
55 phi[4] = 0;
56 rho[4] = 0;
57 phi[5] = 0;
58 rho[5] = 0;
59 phi[6] = 0;
60 rho[6] = 0;
61 phi[7] = 0;
62 rho[7] = 0;
63
64 phi[1] = -3.47995752;
65 phi[2] = -1.249864467;
66 phi[3] = -0.3067504308;
67 phi[4] = 0.9229242379;
68 phi[5] = -3.278567926;
69 phi[6] = -0.6346408622;
70 phi[7] = 1.762377475;
71
72 rho[1] = 2.463898984;
73 rho[2] = 0.7361782665;
74 rho[3] = 1.90216812;
75 rho[4] = 2.140687169;
76 rho[5] = 0.914684056;
77 rho[6] = 1.116206881;
78 rho[7] = 1.483440555;
79
80 modetype[0] = 23;
81 modetype[1] = 23;
82 modetype[2] = 23;
83 modetype[3] = 23;
84 modetype[4] = 23;
85 modetype[5] = 13;
86 modetype[6] = 13;
87 modetype[7] = 13;
88
89 // cout << "DsToKpipi :" << endl;
90 // for (int i=0; i<8; i++) {
91 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
92 // }
93
94 width[0] = 0.1478;
95 width[1] = 0.400;
96 width[2] = 0.100;
97 width[3] = 0.265;
98 width[4] = 0.270;
99 width[5] = 0.0473;
100 width[6] = 0.232;
101 width[7] = 0.270;
102
103 mass[0] = 0.77526;
104 mass[1] = 1.465;
105 mass[2] = 0.965;
106 mass[3] = 1.35;
107 mass[4] = 1.425;
108 mass[5] = 0.89555;
109 mass[6] = 1.414;
110 mass[7] = 1.432787726;
111
112 mDsM = 1.9683;
113 mD = 1.86486;
114 mDs = 1.9683;
115 // rRes = 9.0;
116 rD = 25.0;
117 metap = 0.95778;
118 mkstr = 0.89594;
119 mk0 = 0.497614;
120 mass_Kaon = 0.49368;
121 mass_Pion = 0.13957;
122 mass_Pi0 = 0.1349766;
123 math_pi = 3.1415926;
124 mKa2 = 0.24371994;
125 mPi2 = 0.01947978;
126 mPi = 0.13957;
127 mKa = 0.49368;
128 // mrho = 0.77549;
129 // Grho = 0.1491;
130 ma0 = 0.99;
131 Ga0 = 0.0756;
132 meta = 0.547862;
133
134 GS1 = 0.636619783;
135 GS2 = 0.01860182466;
136 GS3 = 0.1591549458; // 1/(2*math_2pi)
137 GS4 = 0.00620060822; // mass_Pion2/math_pi
138
139 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
140 for ( int i = 0; i < 4; i++ )
141 {
142 for ( int j = 0; j < 4; j++ ) { G[i][j] = GG[i][j]; }
143 }
144}
145
146void EvtDsToKpipi::initProbMax() { setProbMax( 825.0 ); }
147
148void EvtDsToKpipi::decay( EvtParticle* p ) {
149 /*
150 double maxprob = 0.0;
151 for(int ir=0;ir<=60000000;ir++){
152 p->initializePhaseSpace(getNDaug(),getDaugs());
153 EvtVector4R D1 = p->getDaug(0)->getP4();
154 EvtVector4R D2 = p->getDaug(1)->getP4();
155 EvtVector4R D3 = p->getDaug(2)->getP4();
156
157 double P1[4], P2[4], P3[4];
158 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
159 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
160 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
161
162 double value;
163 int g0[8]={1,1,4,2,500,1,1,2};
164 int spin[8]={1,1,0,0,0,1,1,0};
165 int nstates=8;
166 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
167 if (value<0) continue;
168 if(value>maxprob) {
169 maxprob=value;
170 cout << "ir= " << ir << endl;
171 cout << "double P1[4] = {" << P1[0] <<","<< P1[1] <<","<< P1[2] <<","<< P1[3]
172 <<"};"<< endl; cout << "double P2[4] = {" << P2[0] <<","<< P2[1] <<","<< P2[2] <<","<<
173 P2[3] <<"};"<< endl; cout << "double P3[4] = {" << P3[0] <<","<< P3[1] <<","<< P3[2]
174 <<","<< P3[3] <<"};"<< endl; cout << "MAX====> " << maxprob << endl;
175 }
176 }
177 printf("MAXprob = %.10f\n",maxprob);
178 */
179
180 p->initializePhaseSpace( getNDaug(), getDaugs() );
181 EvtVector4R D1 = p->getDaug( 0 )->getP4();
182 EvtVector4R D2 = p->getDaug( 1 )->getP4();
183 EvtVector4R D3 = p->getDaug( 2 )->getP4();
184
185 double P1[4], P2[4], P3[4];
186 P1[0] = D1.get( 0 );
187 P1[1] = D1.get( 1 );
188 P1[2] = D1.get( 2 );
189 P1[3] = D1.get( 3 );
190 P2[0] = D2.get( 0 );
191 P2[1] = D2.get( 1 );
192 P2[2] = D2.get( 2 );
193 P2[3] = D2.get( 3 );
194 P3[0] = D3.get( 0 );
195 P3[1] = D3.get( 1 );
196 P3[2] = D3.get( 2 );
197 P3[3] = D3.get( 3 );
198
199 // P1[0] =0.863405 ;P1[1] =-0.239724 ;P1[2] =-0.665631 ; P1[3] =-0.0349068 ;
200 // P2[0] =0.209097 ;P2[1] =0.0992859 ;P2[2] =-0.119905 ; P2[3] =-0.00262863 ;
201 // P3[0] =0.941396 ;P3[1] =-0.270638 ;P3[2] =0.890782 ; P3[3] = -0.00293441 ;
202
203 double value;
204 int g0[8] = { 1, 1, 4, 2, 500, 1, 1, 2 };
205 int spin[8] = { 1, 1, 0, 0, 0, 1, 1, 0 };
206 int nstates = 8;
207 calEva( P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value );
208
209 setProb( value );
210
211 return;
212}
213
214void EvtDsToKpipi::Com_Multi( double a1[2], double a2[2], double res[2] ) {
215 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
216 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
217}
218void EvtDsToKpipi::Com_Divide( double a1[2], double a2[2], double res[2] ) {
219 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
220 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
221 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
222}
223
224double EvtDsToKpipi::SCADot( double a1[4], double a2[4] ) {
225 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
226 return _cal;
227}
228double EvtDsToKpipi::barrier( int l, double sa, double sb, double sc, double r, double mass ) {
229 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
230 if ( q < 0 ) q = 1e-16;
231 double z;
232 z = q * r * r;
233 double sa0;
234 sa0 = mass * mass;
235 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
236 if ( q0 < 0 ) q0 = 1e-16;
237 double z0 = q0 * r * r;
238 double F = 0.0;
239 if ( l == 0 ) F = 1;
240 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
241 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
242 return F;
243}
244void EvtDsToKpipi::calt1( double daug1[4], double daug2[4], double t1[4] ) {
245 double p, pq, tmp;
246 double pa[4], qa[4];
247 for ( int i = 0; i < 4; i++ )
248 {
249 pa[i] = daug1[i] + daug2[i];
250 qa[i] = daug1[i] - daug2[i];
251 }
252 p = SCADot( pa, pa );
253 pq = SCADot( pa, qa );
254 tmp = pq / p;
255 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
256}
257void EvtDsToKpipi::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
258 double p, r;
259 double pa[4], t1[4];
260 calt1( daug1, daug2, t1 );
261 r = SCADot( t1, t1 ) / 3.0;
262 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
263 p = SCADot( pa, pa );
264 for ( int i = 0; i < 4; i++ )
265 {
266 for ( int j = 0; j < 4; j++ )
267 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
268 }
269}
270
271//-------------------prop--------------------------------------------
272
273void EvtDsToKpipi::propagatorCBW( double mass, double width, double sx, double prop[2] ) {
274 double a[2], b[2];
275 a[0] = 1;
276 a[1] = 0;
277 b[0] = mass * mass - sx;
278 b[1] = -mass * width;
279 Com_Divide( a, b, prop );
280}
281double EvtDsToKpipi::wid( double mass2, double mass, double sa, double sb, double sc,
282 double r2, int l ) {
283 double widm = 0.;
284 double m = sqrt( sa );
285 double tmp = sb - sc;
286 double tmp1 = sa + tmp;
287 double q = 0.25 * tmp1 * tmp1 / sa - sb;
288 if ( q < 0 ) q = 1e-16;
289 double tmp2 = mass2 + tmp;
290 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
291 if ( q0 < 0 ) q0 = 1e-16;
292 double z = q * r2;
293 double z0 = q0 * r2;
294 double t = q / q0;
295 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
296 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
297 else if ( l == 2 )
298 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
299 return widm;
300}
301
302double EvtDsToKpipi::widl1( double mass2, double mass, double sa, double sb, double sc,
303 double r2 ) {
304 double widm = 0.;
305 double m = sqrt( sa );
306 double tmp = sb - sc;
307 double tmp1 = sa + tmp;
308 double q = 0.25 * tmp1 * tmp1 / sa - sb;
309 if ( q < 0 ) q = 1e-16;
310 double tmp2 = mass2 + tmp;
311 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
312 if ( q0 < 0 ) q0 = 1e-16;
313 double z = q * r2;
314 double z0 = q0 * r2;
315 double F = ( 1 + z0 ) / ( 1 + z );
316 double t = q / q0;
317 widm = t * sqrt( t ) * mass / m * F;
318 return widm;
319}
320void EvtDsToKpipi::propagatorRBW( double mass, double width, double sa, double sb, double sc,
321 double r2, int l, double prop[2] ) {
322 double a[2], b[2];
323 double mass2 = mass * mass;
324
325 a[0] = 1;
326 a[1] = 0;
327 b[0] = mass2 - sa;
328 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
329 Com_Divide( a, b, prop );
330}
331void EvtDsToKpipi::propagatorFlatte( double mass, double width, double sa, double prop[2] ) {
332
333 double q2_Pi, q2_Ka;
334 double rhoPi[2], rhoKa[2];
335
336 q2_Pi = 0.25 * sa - mPi * mPi;
337 q2_Ka = 0.25 * sa - mKa * mKa;
338
339 if ( q2_Pi > 0 )
340 {
341 rhoPi[0] = 2.0 * sqrt( q2_Pi / sa );
342 rhoPi[1] = 0.0;
343 }
344 if ( q2_Pi <= 0 )
345 {
346 rhoPi[0] = 0.0;
347 rhoPi[1] = 2.0 * sqrt( -q2_Pi / sa );
348 }
349
350 if ( q2_Ka > 0 )
351 {
352 rhoKa[0] = 2.0 * sqrt( q2_Ka / sa );
353 rhoKa[1] = 0.0;
354 }
355 if ( q2_Ka <= 0 )
356 {
357 rhoKa[0] = 0.0;
358 rhoKa[1] = 2.0 * sqrt( -q2_Ka / sa );
359 }
360
361 double a[2], b[2];
362 a[0] = 1;
363 a[1] = 0;
364 b[0] = mass * mass - sa + 0.165 * rhoPi[1] + 0.69465 * rhoKa[1];
365 b[1] = -( 0.165 * rhoPi[0] + 0.69465 * rhoKa[0] );
366 Com_Divide( a, b, prop );
367}
368void EvtDsToKpipi::propagatorGS( double mass, double width, double sa, double sb, double sc,
369 double r2, double prop[2] ) {
370 double a[2], b[2];
371 double mass2 = mass * mass;
372 double tmp = sb - sc;
373 double tmp1 = sa + tmp;
374 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
375 if ( q2 < 0 ) q2 = 1e-16;
376
377 double tmp2 = mass2 + tmp;
378 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
379 if ( q02 < 0 ) q02 = 1e-16;
380
381 double q = sqrt( q2 );
382 double q0 = sqrt( q02 );
383 double m = sqrt( sa );
384 double q03 = q0 * q02;
385 double tmp3 = log( mass + 2 * q0 ) + 1.2760418309; // log(mass_2Pion) = 1.2760418309;
386
387 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2760418309 );
388 double h0 = GS1 * q0 / mass * tmp3;
389 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
390 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
391 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
392
393 a[0] = 1.0 + d * width / mass;
394 a[1] = 0.0;
395 b[0] = mass2 - sa + width * f;
396 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
397 Com_Divide( a, b, prop );
398}
399
400void EvtDsToKpipi::Flatte_rhoab( double sa, double sb, double sc, double rho[2] ) {
401 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
402 if ( q > 0 )
403 {
404 rho[0] = 2 * sqrt( q / sa );
405 rho[1] = 0;
406 }
407 else if ( q < 0 )
408 {
409 rho[0] = 0;
410 rho[1] = 2 * sqrt( -q / sa );
411 }
412}
413
414void EvtDsToKpipi::propagatorKstr1430( double mass, double sx, double* sb, double* sc,
415 double prop[2] ) // K*1430 Flatte
416{
417 double unit[2] = { 1.0 };
418 double ci[2] = { 0, 1 };
419 double rho1[2];
420 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
421 double rho2[2];
422 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
423 double gKPi_Kstr1430 = 0.2990, gEtaPK_Kstr1430 = 0.0529;
424 double tmp1[2] = { gKPi_Kstr1430, 0 };
425 double tmp11[2];
426 double tmp2[2] = { gEtaPK_Kstr1430, 0 };
427 double tmp22[2];
428 Com_Multi( tmp1, rho1, tmp11 );
429 Com_Multi( tmp2, rho2, tmp22 );
430 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
431 double tmp31[2];
432 Com_Multi( tmp3, ci, tmp31 );
433 double tmp4[2] = { mass * mass - sx - tmp31[0], -1.0 * tmp31[1] };
434 Com_Divide( unit, tmp4, prop );
435}
436
437double EvtDsToKpipi::DDalitz( double P1[4], double P2[4], double P3[4], int Ang,
438 double mass ) {
439 double pR[4], pD[4];
440 double temp_PDF, v_re;
441 temp_PDF = 0.0;
442 v_re = 0.0;
443 double B[2], s1, s2, s3, sR, sD;
444 for ( int i = 0; i < 4; i++ )
445 {
446 pR[i] = P1[i] + P2[i];
447 pD[i] = pR[i] + P3[i];
448 }
449 s1 = SCADot( P1, P1 );
450 s2 = SCADot( P2, P2 );
451 s3 = SCADot( P3, P3 );
452 sR = SCADot( pR, pR );
453 sD = SCADot( pD, pD );
454 int G[4][4];
455 for ( int i = 0; i != 4; i++ )
456 {
457 for ( int j = 0; j != 4; j++ )
458 {
459 if ( i == j )
460 {
461 if ( i == 0 ) G[i][j] = 1;
462 else G[i][j] = -1;
463 }
464 else G[i][j] = 0;
465 }
466 }
467 if ( Ang == 0 )
468 {
469 B[0] = 1;
470 B[1] = 1;
471 temp_PDF = 1;
472 }
473 if ( Ang == 1 )
474 {
475 B[0] = barrier( 1, sR, s1, s2, 3.0, mass );
476 B[1] = barrier( 1, sD, sR, s3, 5.0, mDsM );
477 // B[0] = Barrier(1,sR,s1,s2,9.0);
478 // B[1] = Barrier(1,sD,sR,s3,25.0);
479 double t1[4], T1[4];
480 calt1( P1, P2, t1 );
481 calt1( pR, P3, T1 );
482 temp_PDF = 0;
483 for ( int i = 0; i < 4; i++ ) { temp_PDF += t1[i] * T1[i] * G[i][i]; }
484 }
485 if ( Ang == 2 )
486 {
487 B[0] = barrier( 2, sR, s1, s2, 3.0, mass );
488 B[1] = barrier( 2, sD, sR, s3, 5.0, mDsM );
489 double t2[4][4], T2[4][4];
490 calt2( P1, P2, t2 );
491 calt2( pR, P3, T2 );
492 temp_PDF = 0;
493 for ( int i = 0; i < 4; i++ )
494 {
495 for ( int j = 0; j < 4; j++ ) { temp_PDF += t2[i][j] * T2[j][i] * G[i][i] * G[j][j]; }
496 }
497 }
498 v_re = temp_PDF * B[0] * B[1];
499 return v_re;
500}
501
502void EvtDsToKpipi::rhoab( double sa, double sb, double sc, double res[2] ) {
503 double tmp = sa + sb - sc;
504 double q = 0.25 * tmp * tmp / sa - sb;
505 if ( q >= 0 )
506 {
507 res[0] = 2.0 * sqrt( q / sa );
508 res[1] = 0.0;
509 }
510 else
511 {
512 res[0] = 0.0;
513 res[1] = 2.0 * sqrt( -q / sa );
514 }
515}
516void EvtDsToKpipi::rho4Pi( double sa, double res[2] ) {
517 double temp = 1.0 - 0.3116765584 / sa; // 0.3116765584=0.13957*0.13957*16
518 if ( temp >= 0 )
519 {
520 res[0] = sqrt( temp ) / ( 1.0 + exp( 9.8 - 3.5 * sa ) );
521 res[1] = 0.0;
522 }
523 else
524 {
525 res[0] = 0.0;
526 res[1] = sqrt( -temp ) / ( 1.0 + exp( 9.8 - 3.5 * sa ) ); //?????????????????????
527 }
528}
529
530void EvtDsToKpipi::propagatorsigma500( double sa, double sb, double sc,
531 double prop[2] ) { // for pipi_s wave
532 double f = 0.5843 + 1.6663 * sa;
533 const double M = 0.9264;
534 const double mass2 = 0.85821696; // M*M
535 const double mpi2d2 = 0.00973989245;
536 double g1 = f * ( sa - mpi2d2 ) / ( mass2 - mpi2d2 ) * exp( ( mass2 - sa ) / 1.082 );
537 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
538 rhoab( sa, sb, sc, rho1s );
539 rhoab( mass2, sb, sc, rho1M );
540 rho4Pi( sa, rho2s );
541 rho4Pi( mass2, rho2M );
542 Com_Divide( rho1s, rho1M, rho1 );
543 Com_Divide( rho2s, rho2M, rho2 );
544 double a[2], b[2];
545 a[0] = 1.0;
546 a[1] = 0.0;
547 b[0] = mass2 - sa + M * ( g1 * rho1[1] + 0.0024 * rho2[1] );
548 b[1] = -M * ( g1 * rho1[0] + 0.0024 * rho2[0] );
549 Com_Divide( a, b, prop );
550}
551
552void EvtDsToKpipi::calEva( double* K, double* Pi1, double* Pi2, double* mass1, double* width1,
553 double* amp, double* phase, int* g0, int* spin, int* modetype,
554 int nstates, double& Result ) {
555 // double Ks0[4], Pic[4], Pi0[4];
556 double rRes = 3.0;
557 double rRess = 9.0;
558 double P23[4], P13[4];
559 double cof[2], amp_PDF[2], PDF[2];
560 // double scpi, sks02, sks01;
561 // double s12, s13, s23;
562 double s13, s23;
563 for ( int i = 0; i < 4; i++ )
564 {
565 // P12[i] = Ks01[i] + Ks02[i];
566 P13[i] = K[i] + Pi2[i];
567 P23[i] = Pi1[i] + Pi2[i];
568 }
569 s13 = SCADot( P13, P13 );
570 s23 = SCADot( P23, P23 );
571 double s1, s2, s3;
572 s1 = SCADot( K, K );
573 s2 = SCADot( Pi1, Pi1 );
574 s3 = SCADot( Pi2, Pi2 );
575 double pro[2], temp_PDF, amp_tmp[2], temp_PDF1, temp_PDF2, pro1[2], pro2[2];
576 // double mass1sq;
577 amp_PDF[0] = 0;
578 amp_PDF[1] = 0;
579 PDF[0] = 0;
580 PDF[1] = 0;
581 amp_tmp[0] = 0;
582 amp_tmp[1] = 0;
583 for ( int i = 0; i < nstates; i++ )
584 {
585 amp_tmp[0] = 0;
586 amp_tmp[1] = 0;
587 // mass1sq = mass1[i]*mass1[i];
588 cof[0] = amp[i] * cos( phase[i] );
589 cof[1] = amp[i] * sin( phase[i] );
590 temp_PDF = 0;
591
592 if ( modetype[i] == 13 )
593 {
594 temp_PDF = DDalitz( K, Pi2, Pi1, spin[i], mass1[i] );
595 if ( g0[i] == 1 )
596 propagatorRBW( mass1[i], width1[i], s13, mKa2, mPi2, rRess, spin[i], pro );
597 if ( g0[i] == 2 )
598 { // K*1430 Flatte
599 // double skm2[2]={s1, mass_EtaP *mass_EtaP};
600 double skm2[2] = { s1, 0.95778 * 0.95778 };
601 double spi2[2] = { s3, mass_Kaon * mass_Kaon };
602 propagatorKstr1430( mass1[i], s13, skm2, spi2, pro );
603 }
604 if ( g0[i] == 0 )
605 {
606 pro[0] = 1;
607 pro[1] = 0;
608 }
609 amp_tmp[0] = temp_PDF * pro[0];
610 amp_tmp[1] = temp_PDF * pro[1];
611 }
612
613 if ( modetype[i] == 23 )
614 {
615 temp_PDF = DDalitz( Pi1, Pi2, K, spin[i], mass1[i] );
616 if ( g0[i] == 4 ) propagatorFlatte( mass1[i], width1[i], s23, pro ); // Only for f0(980)
617 if ( g0[i] == 3 ) propagatorCBW( mass1[i], width1[i], s23, pro );
618 if ( g0[i] == 2 )
619 propagatorRBW( mass1[i], width1[i], s23, mPi2, mPi2, rRess, spin[i], pro );
620 if ( g0[i] == 1 ) propagatorGS( mass1[i], width1[i], s23, mPi2, mPi2, 9.0, pro );
621 if ( g0[i] == 500 ) propagatorsigma500( s23, s2, s3, pro );
622 if ( g0[i] == 0 )
623 {
624 pro[0] = 1;
625 pro[1] = 0;
626 }
627 amp_tmp[0] = temp_PDF * pro[0];
628 amp_tmp[1] = temp_PDF * pro[1];
629 }
630 Com_Multi( amp_tmp, cof, amp_PDF );
631 PDF[0] += amp_PDF[0];
632 PDF[1] += amp_PDF[1];
633 // cout<<"PDF: "<<PDF[0]<<" "<<PDF[1]<<endl;
634 // cout<<"amp_PDF: "<<amp_PDF[0]<<" "<<amp_PDF[1]<<endl;
635 }
636
637 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
638 // cout<<"value = "<<value<<endl;
639 if ( value <= 0 ) value = 1e-20;
640 Result = value;
641 // cout<<"value = "<<value<<" Result = "<<Result<<endl;
642}
mass
double mass
Definition CosmicGenerator.cxx:128
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
g1
TF1 * g1
Definition DataBase/tau_mode.c:58
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:209
unit
*******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 !Latex Output unit
Definition FoamA.h:90
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
EvtDecayBase::checkSpinDaughter
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
Definition EvtDecayBase.cc:492
EvtDecayBase::checkSpinParent
void checkSpinParent(EvtSpinType::spintype sp)
Definition EvtDecayBase.cc:479
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
EvtDsToKpipi::decay
void decay(EvtParticle *p)
Definition EvtDsToKpipi.cc:148
EvtDsToKpipi::getName
void getName(std::string &name)
Definition EvtDsToKpipi.cc:34
EvtDsToKpipi::clone
EvtDecayBase * clone()
Definition EvtDsToKpipi.cc:36
EvtDsToKpipi::~EvtDsToKpipi
virtual ~EvtDsToKpipi()
Definition EvtDsToKpipi.cc:32
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
t
int t()
Definition t.c:1