BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDToKSKSK.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// Environment:
3// This software is part of models developed at BES collaboration
4// based on the EvtGen framework. If you use all or part
5// of it, please give an appropriate acknowledgement.
6//
7// Copyright Information: See EvtGen/BesCopyright
8// Copyright (A) 2006 Ping Rong-Gang @IHEP
9//
10// Module: EvtDToKSKSK.cc
11//
12// Description: D+ -> KS0 KS0 K+
13//
14// Modification history:
15//
16// Liaoyuan Dong Wed Nov 9 00:25:20 2022 Module created
17//------------------------------------------------------------------------
18#include "EvtDToKSKSK.hh"
19#include "../EvtGenBase/EvtComplex.hh"
20#include "../EvtGenBase/EvtDecayTable.hh"
21#include "../EvtGenBase/EvtGenKine.hh"
22#include "../EvtGenBase/EvtPDL.hh"
23#include "../EvtGenBase/EvtParticle.hh"
24#include "../EvtGenBase/EvtPatches.hh"
25#include "../EvtGenBase/EvtReport.hh"
26#include <cmath>
27#include <iostream>
28#include <stdlib.h>
29using namespace std;
30
31EvtDToKSKSK::~EvtDToKSKSK() {}
32
33void EvtDToKSKSK::getName( std::string& model_name ) { model_name = "DToKSKSK"; }
34
35EvtDecayBase* EvtDToKSKSK::clone() { return new EvtDToKSKSK; }
36
37void EvtDToKSKSK::init() {
38 checkNArg( 0 );
39 checkNDaug( 3 );
40 checkSpinParent( EvtSpinType::SCALAR );
41
42 mass[0] = 0.980;
43 mass[1] = 0.980;
44 width[0] = 0.0972;
45 width[1] = 0.0972;
46
47 rho[0] = 1.0;
48 rho[1] = 1.1141e+01;
49 phi[0] = 0.00;
50 phi[1] = -5.6220e+00;
51
52 // spin[0] = 0;
53 // spin[1] = 1;
54 // spin[2] = 1;
55 modetype[0] = 12;
56 modetype[1] = 122;
57
58 // std::cout << "Initializing EvtDToKSKSK" << std::endl;
59 // for (int i=0; i<2; i++) {
60 // cout << i << " rho,phi = " << rho[i] << ", "<< phi[i] << endl;
61 // }
62
63 math_K = 0.493677;
64 GS1 = 0.636619783;
65 GS2 = 0.23459086;
66 GS3 = 0.1591549458;
67 GS4 = 0.078196955;
68
69 mD0 = 1.86966;
70 mK0 = 0.497611;
71 mKa = 0.493677;
72 mPi = 0.13957;
73 mK02 = 0.237616707;
74 mPi2 = 0.01947978;
75 mass_EtaP = 0.95778;
76 mass_Kaon = 0.49368;
77
78 math_pi = 3.1415926;
79 mass_Pion = 0.13957;
80 mass_Pion2 = 0.0194797849;
81 mass_2Pion = 0.27914;
82 math_2pi = 6.2831852;
83 rD2 = 25.0;
84 rRes2 = 9.0;
85 g2 = 0.23;
86
87 rho_omega = 0.00294;
88 phi_omega = -0.02;
89}
90
91void EvtDToKSKSK::initProbMax() { setProbMax( 878.9 ); }
92
93void EvtDToKSKSK::decay( EvtParticle* p ) {
94
95 /* double maxprob = 0.0;
96 for(int ir=0;ir<=60000000;ir++){
97 p->initializePhaseSpace(getNDaug(),getDaugs());
98 EvtVector4R D1 = p->getDaug(0)->getP4();
99 EvtVector4R D2 = p->getDaug(1)->getP4();
100 EvtVector4R D3 = p->getDaug(2)->getP4();
101
102 double P1[4], P2[4], P3[4];
103 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
104 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
105 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
106
107 double value;
108 int g0[3]={1,1,1};
109 int spin[3]={0,1,0};
110 int nstates=2;
111 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates,
112 value); if (value<0) continue; if(value>maxprob) { maxprob=value; cout << "ir= " << ir <<
113 endl; cout << "double P1[4] = {" << P1[0] <<","<< P1[1] <<","<< P1[2] <<","<< P1[3]
114 <<"};"<< endl; cout << "double P2[4] = {" << P2[0] <<","<< P2[1] <<","<< P2[2] <<","<<
115 P2[3] <<"};"<< endl; cout << "double P3[4] = {" << P3[0] <<","<< P3[1] <<","<< P3[2]
116 <<","<< P3[3] <<"};"<< endl; cout << "MAX====> " << maxprob << endl;
117 }
118 }
119 printf("MAXprob = %.10f\n",maxprob);
120
121 */
122 p->initializePhaseSpace( getNDaug(), getDaugs() );
123 EvtVector4R D1 = p->getDaug( 0 )->getP4();
124 EvtVector4R D2 = p->getDaug( 1 )->getP4();
125 EvtVector4R D3 = p->getDaug( 2 )->getP4();
126
127 double P1[4], P2[4], P3[4];
128 P1[0] = D1.get( 0 );
129 P1[1] = D1.get( 1 );
130 P1[2] = D1.get( 2 );
131 P1[3] = D1.get( 3 );
132 P2[0] = D2.get( 0 );
133 P2[1] = D2.get( 1 );
134 P2[2] = D2.get( 2 );
135 P2[3] = D2.get( 3 );
136 P3[0] = D3.get( 0 );
137 P3[1] = D3.get( 1 );
138 P3[2] = D3.get( 2 );
139 P3[3] = D3.get( 3 );
140
141 double value;
142 int g0[2] = { 1, 1 };
143 int spin[2] = { 0, 0 };
144 int nstates = 2;
145 calEva( P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value );
146 setProb( value );
147
148 return;
149}
150
151void EvtDToKSKSK::Com_Multi( double a1[2], double a2[2], double res[2] ) {
152 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
153 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
154}
155void EvtDToKSKSK::Com_Divide( double a1[2], double a2[2], double res[2] ) {
156 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
157 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
158 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
159}
160//------------base---------------------------------
161double EvtDToKSKSK::SCADot( double a1[4], double a2[4] ) {
162 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
163 return _cal;
164}
165double EvtDToKSKSK::barrier( int l, double sa, double sb, double sc, double r, double mass ) {
166 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
167 // if(q < 0) q = 1e-16;
168 double z;
169 z = q * r * r;
170 double sa0;
171 sa0 = mass * mass;
172 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
173 // if(q0 < 0) q0 = 1e-16;
174 double z0 = q0 * r * r;
175 double F = 0.0;
176 if ( l == 0 ) F = 1;
177 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
178 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
179 return F;
180}
181
182void EvtDToKSKSK::calt1( double daug1[4], double daug2[4], double t1[4] ) {
183 double p, pq, tmp;
184 double pa[4], qa[4];
185 for ( int i = 0; i < 4; i++ )
186 {
187 pa[i] = daug1[i] + daug2[i];
188 qa[i] = daug1[i] - daug2[i];
189 }
190 p = SCADot( pa, pa );
191 pq = SCADot( pa, qa );
192 tmp = pq / p;
193 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
194}
195
196void EvtDToKSKSK::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
197 double p, r;
198 double pa[4], t1[4];
199 calt1( daug1, daug2, t1 );
200 r = SCADot( t1, t1 ) / 3.0;
201 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
202 p = SCADot( pa, pa );
203 for ( int i = 0; i < 4; i++ )
204 {
205 for ( int j = 0; j < 4; j++ )
206 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
207 }
208}
209
210void EvtDToKSKSK::propagatorCBW( double mass, double width, double sx, double prop[2] ) {
211 double a[2], b[2];
212 a[0] = 1;
213 a[1] = 0;
214 b[0] = mass * mass - sx;
215 b[1] = -mass * width;
216 Com_Divide( a, b, prop );
217}
218
219double EvtDToKSKSK::wid( double mass2, double mass, double sa, double sb, double sc, double r2,
220 int l ) {
221 double widm = 0.;
222 double m = sqrt( sa );
223 double tmp = sb - sc;
224 double tmp1 = sa + tmp;
225 double q = fabs( 0.25 * tmp1 * tmp1 / sa - sb );
226 double tmp2 = mass2 + tmp;
227 double q0 = fabs( 0.25 * tmp2 * tmp2 / mass2 - sb );
228 double z = q * r2;
229 double z0 = q0 * r2;
230 double t = q / q0;
231
232 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
233 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
234 else if ( l == 2 )
235 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
236 return widm;
237}
238double EvtDToKSKSK::widl1( double mass2, double mass, double sa, double sb, double sc,
239 double r2 ) {
240 double widm = 0.;
241 double m = sqrt( sa );
242 double tmp = sb - sc;
243 double tmp1 = sa + tmp;
244 double q = fabs( 0.25 * tmp1 * tmp1 / sa - sb );
245 double tmp2 = mass2 + tmp;
246 double q0 = fabs( 0.25 * tmp2 * tmp2 / mass2 - sb );
247 double z = q * r2;
248 double z0 = q0 * r2;
249 double F = ( 1 + z0 ) / ( 1 + z );
250 double t = q / q0;
251 widm = t * sqrt( t ) * mass / m * F;
252 return widm;
253}
254void EvtDToKSKSK::propagatorRBW( double mass, double width, double sa, double sb, double sc,
255 double r2, int l, double prop[2] ) {
256 double a[2], b[2];
257 double mass2 = mass * mass;
258
259 a[0] = 1;
260 a[1] = 0;
261 b[0] = mass2 - sa;
262 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
263 Com_Divide( a, b, prop );
264}
265
266void EvtDToKSKSK::propagatorFlatte( double mass, double width, double sa, double sb, double sc,
267 double prop[2] ) {
268 double rhoab, rhoKK;
269 rhoab = 1.0 / sa *
270 sqrt( ( sa - ( 0.547862 + 0.139570 ) * ( 0.547862 + 0.139570 ) ) *
271 ( sa - ( 0.547862 - 0.139570 ) * ( 0.547862 - 0.139570 ) ) );
272 rhoKK = 1.0 / sa *
273 sqrt( ( sa - ( 0.497611 + 0.493677 ) * ( 0.497611 + 0.493677 ) ) *
274 ( sa - ( 0.497611 - 0.493677 ) * ( 0.497611 - 0.493677 ) ) );
275 double a[2], b[2];
276 a[0] = 1;
277 a[1] = 0;
278 b[0] = mass * mass - sa;
279 b[1] = -( 0.324 * 0.324 * rhoab + 1.03 * 0.324 * 0.324 * rhoKK );
280 Com_Divide( a, b, prop );
281}
282
283void EvtDToKSKSK::propagatorGS( double mass, double width, double sa, double sb, double sc,
284 double r2, double prop[2] ) {
285 double a[2], b[2];
286 double mass2 = mass * mass;
287 double tmp = sb - sc;
288 double tmp1 = sa + tmp;
289 double q2 = fabs( 0.25 * tmp1 * tmp1 / sa - sb );
290 double tmp2 = mass2 + tmp;
291 double q02 = fabs( 0.25 * tmp2 * tmp2 / mass2 - sb );
292 double q = sqrt( q2 );
293 double q0 = sqrt( q02 );
294 double m = sqrt( sa );
295 double q03 = q0 * q02;
296 double tmp3 = log( mass + 2 * q0 ) + 0.0087501713;
297 double h = GS1 * q / m * ( log( m + 2 * q ) + 0.0087501713 );
298 double h0 = GS1 * q0 / mass * tmp3;
299 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
300 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
301 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
302 a[0] = 1.0 + d * width / mass;
303 a[1] = 0.0;
304 b[0] = mass2 - sa + width * f;
305 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
306 Com_Divide( a, b, prop );
307}
308
309void EvtDToKSKSK::PiPiSWASS( double sa, double sb, double sc, double prop[2] ) {
310 double tmp = sb - sc;
311 double tmp2 = sa + tmp;
312 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
313 double q = sqrt( qs );
314 double a0 = -0.11 / mass_Pion;
315 prop[0] = 1 / ( 1 + a0 * a0 * q * q );
316 prop[1] = a0 * q / ( 1 + a0 * a0 * q * q );
317}
318
319void EvtDToKSKSK::KPiSLASS( double sa, double sb, double sc, double prop[2] ) {
320 const double m1430 = 1.441;
321 const double sa0 = 2.076481;
322 const double w1430 = 0.193;
323 const double Lass1 = 0.25 / sa0;
324 double tmp = sb - sc;
325 double tmp1 = sa0 + tmp;
326 double q0 = fabs( Lass1 * tmp1 * tmp1 - sb );
327 double tmp2 = sa + tmp;
328 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
329 double q = sqrt( qs );
330 double width = w1430 * q * m1430 / sqrt( sa * q0 );
331 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
332 if ( temp_R < 0 ) temp_R += math_K;
333 double deltaR = -109.7 * math_K / 180.0 + temp_R;
334 double temp_F = atan( 0.226 * q / ( 2.0 - 3.8194 * qs ) );
335 if ( temp_F < 0 ) temp_F += math_K;
336 double deltaF = 0.1 * math_K / 180.0 + temp_F;
337 double deltaS = deltaR + 2.0 * deltaF;
338 double t1 = 0.96 * sin( deltaF );
339 double t2 = sin( deltaR );
340 double CF[2], CS[2];
341 CF[0] = cos( deltaF );
342 CF[1] = sin( deltaF );
343 CS[0] = cos( deltaS );
344 CS[1] = sin( deltaS );
345 prop[0] = t1 * CF[0] + t2 * CS[0];
346 prop[1] = t1 * CF[1] + t2 * CS[1];
347}
348
349void EvtDToKSKSK::Flatte_rhoab( double sa, double sb, double sc, double rho[2] ) {
350 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
351 if ( q > 0 )
352 {
353 rho[0] = 2 * sqrt( q / sa );
354 rho[1] = 0;
355 }
356 else if ( q < 0 )
357 {
358 rho[0] = 0;
359 rho[1] = 2 * sqrt( -q / sa );
360 }
361}
362
363void EvtDToKSKSK::propagatorKstr1430( double mass, double sx, double* sb, double* sc,
364 double prop[2] ) {
365 double unit[2] = { 1.0 };
366 double ci[2] = { 0, 1 };
367 double rho1[2];
368 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
369 double rho2[2];
370 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
371 double gKPi_Kstr1430 = 0.2990, gEtaPK_Kstr1430 = 0.0529;
372 double tmp1[2] = { gKPi_Kstr1430, 0 };
373 double tmp11[2];
374 double tmp2[2] = { gEtaPK_Kstr1430, 0 };
375 double tmp22[2];
376 Com_Multi( tmp1, rho1, tmp11 );
377 Com_Multi( tmp2, rho2, tmp22 );
378 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
379 double tmp31[2];
380 Com_Multi( tmp3, ci, tmp31 );
381 double tmp4[2] = { mass * mass - sx - tmp31[0], -1.0 * tmp31[1] };
382 Com_Divide( unit, tmp4, prop );
383}
384void EvtDToKSKSK::rhoab( double sa, double sb, double sc, double res[2] ) {
385 double tmp = sa + sb - sc;
386 double q = 0.25 * tmp * tmp / sa - sb;
387 if ( q >= 0 )
388 {
389 res[0] = 2.0 * sqrt( q / sa );
390 res[1] = 0.0;
391 }
392 else
393 {
394 res[0] = 0.0;
395 res[1] = 2.0 * sqrt( -q / sa );
396 }
397}
398void EvtDToKSKSK::rho4Pi( double sa, double res[2] ) {
399 double temp = 1.0 - 0.3116765584 / sa;
400 if ( temp >= 0 )
401 {
402 res[0] = sqrt( temp ) / ( 1.0 + exp( 9.8 - 3.5 * sa ) );
403 res[1] = 0.0;
404 }
405 else
406 {
407 res[0] = 0.0;
408 res[1] = sqrt( -temp ) / ( 1.0 + exp( 9.8 - 3.5 * sa ) );
409 }
410}
411
412void EvtDToKSKSK::propagatorsigma500( double sa, double sb, double sc, double prop[2] ) {
413 double f = 0.5843 + 1.6663 * sa;
414 const double M = 0.9264;
415 const double mass2 = 0.85821696;
416 const double mpi2d2 = 0.00973989245;
417 double g1 = f * ( sa - mpi2d2 ) / ( mass2 - mpi2d2 ) * exp( ( mass2 - sa ) / 1.082 );
418 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
419 rhoab( sa, sb, sc, rho1s );
420 rhoab( mass2, sb, sc, rho1M );
421 rho4Pi( sa, rho2s );
422 rho4Pi( mass2, rho2M );
423 Com_Divide( rho1s, rho1M, rho1 );
424 Com_Divide( rho2s, rho2M, rho2 );
425 double a[2], b[2];
426 a[0] = 1.0;
427 a[1] = 0.0;
428 b[0] = mass2 - sa + M * ( g1 * rho1[1] + 0.0024 * rho2[1] );
429 b[1] = -M * ( g1 * rho1[0] + 0.0024 * rho2[0] );
430 Com_Divide( a, b, prop );
431}
432
433void EvtDToKSKSK::getprop( double sa, double sb, double sc, double mass, double width,
434 double prop[2] ) {
435 double prop1[2], prop2[2];
436 propagatorGS( mass, width, sa, sb, sc, 9.0, prop1 );
437 propagatorRBW( 0.783, 0.008, sa, sb, sc, 3.0, 1, prop2 );
438 double coef_omega[2];
439 coef_omega[0] = rho_omega * cos( phi_omega ), coef_omega[1] = rho_omega * sin( phi_omega );
440 double one[2];
441 one[0] = 1;
442 one[1] = 0;
443 double temp[2];
444 Com_Multi( coef_omega, prop2, temp );
445 temp[0] = one[0] + 0.783 * 0.783 * temp[0];
446 temp[1] = one[1] + 0.783 * 0.783 * temp[1];
447 Com_Multi( prop1, temp, prop );
448}
449
450double EvtDToKSKSK::DDalitz( double P1[4], double P2[4], double P3[4], int Ang, double mass ) {
451 double pR[4], pD[4];
452 double temp_PDF, v_re;
453 temp_PDF = 0.0;
454 v_re = 0.0;
455 double B[2], s1, s2, s3, sR, sD;
456 for ( int i = 0; i < 4; i++ )
457 {
458 pR[i] = P1[i] + P2[i];
459 pD[i] = pR[i] + P3[i];
460 }
461 s1 = SCADot( P1, P1 );
462 s2 = SCADot( P2, P2 );
463 s3 = SCADot( P3, P3 );
464 sR = SCADot( pR, pR );
465 sD = SCADot( pD, pD );
466 int G[4][4];
467 for ( int i = 0; i != 4; i++ )
468 {
469 for ( int j = 0; j != 4; j++ )
470 {
471 if ( i == j )
472 {
473 if ( i == 0 ) G[i][j] = 1;
474 else G[i][j] = -1;
475 }
476 else G[i][j] = 0;
477 }
478 }
479 if ( Ang == 0 )
480 {
481 B[0] = 1;
482 B[1] = 1;
483 temp_PDF = 1;
484 }
485 if ( Ang == 1 )
486 {
487 B[0] = barrier( 1, sR, s1, s2, 3.0, mass );
488 B[1] = barrier( 1, sD, sR, s3, 5.0, mD0 );
489 double t1[4], T1[4];
490 calt1( P1, P2, t1 );
491 calt1( pR, P3, T1 );
492 temp_PDF = 0;
493 for ( int i = 0; i < 4; i++ ) { temp_PDF += t1[i] * T1[i] * G[i][i]; }
494 }
495 if ( Ang == 2 )
496 {
497 B[0] = barrier( 2, sR, s1, s2, 3.0, mass );
498 B[1] = barrier( 2, sD, sR, s3, 5.0, mD0 );
499 double t2[4][4], T2[4][4];
500 calt2( P1, P2, t2 );
501 calt2( pR, P3, T2 );
502 temp_PDF = 0;
503 for ( int i = 0; i < 4; i++ )
504 {
505 for ( int j = 0; j < 4; j++ ) { temp_PDF += t2[i][j] * T2[j][i] * G[i][i] * G[j][j]; }
506 }
507 }
508 v_re = temp_PDF * B[0] * B[1];
509 return v_re;
510}
511
512void EvtDToKSKSK::calEva( double* Ks01, double* Ks02, double* Kp, double* mass1,
513 double* width1, double* amp, double* phase, int* g0, int* spin,
514 int* modetype, int nstates, double& Result ) {
515 double numEvents;
516 double P23[4], P13[4], P12[4];
517 // double Ks01[4], Ks02[4], Kp[4], P23[4], P13[4], P12[4];
518 double cof[2], amp_PDF[2], PDF[2];
519 double s12, s23, s13;
520 for ( int i = 0; i < 4; i++ )
521 {
522 P12[i] = Ks01[i] + Ks02[i];
523 P13[i] = Ks01[i] + Kp[i];
524 P23[i] = Ks02[i] + Kp[i];
525 }
526 s12 = SCADot( P12, P12 );
527 s13 = SCADot( P13, P13 );
528 s23 = SCADot( P23, P23 );
529 double s1, s2, s3;
530 s1 = SCADot( Ks01, Ks01 );
531 s2 = SCADot( Ks02, Ks02 );
532 s3 = SCADot( Kp, Kp );
533 double pro[2], temp_PDF, amp_tmp[2];
534 double pro1[2], temp_PDF1, pro2[2], temp_PDF2;
535 double pro3[2], temp_PDF3, pro4[2], temp_PDF4;
536 double pro5[2], temp_PDF5, pro6[2], temp_PDF6;
537 double pro7[2], temp_PDF7, pro8[2], temp_PDF8;
538 double pro9[2], temp_PDF9, pro10[2], temp_PDF10;
539 double Amp_KPiS[2];
540 amp_PDF[0] = 0;
541 amp_PDF[1] = 0;
542 PDF[0] = 0;
543 PDF[1] = 0;
544 amp_tmp[0] = 0;
545 amp_tmp[1] = 0;
546 for ( int i = 0; i < 3; i++ )
547 {
548 amp_tmp[0] = 0;
549 amp_tmp[1] = 0;
550 cof[0] = amp[i] * cos( phase[i] );
551 cof[1] = amp[i] * sin( phase[i] );
552 temp_PDF = 0;
553 if ( modetype[i] == 12 )
554 {
555 temp_PDF1 = DDalitz( Ks01, Kp, Ks02, spin[i], mass1[i] );
556 if ( g0[i] == 1 ) propagatorFlatte( mass1[i], width1[i], s13, s1, s3, pro1 );
557
558 if ( g0[i] == 2 )
559 {
560 double skm2[2] = { s1, mass_EtaP * mass_EtaP };
561 double spi2[2] = { s2, mass_Kaon * mass_Kaon };
562 propagatorKstr1430( mass1[i], s12, skm2, spi2, pro1 );
563 }
564 if ( g0[i] == 4 ) KPiSLASS( s12, s1, s2, pro1 );
565 if ( g0[i] == 0 )
566 {
567 pro1[0] = 1;
568 pro1[1] = 0;
569 }
570
571 temp_PDF2 = DDalitz( Ks02, Kp, Ks01, spin[i], mass1[i] );
572 if ( g0[i] == 1 ) propagatorFlatte( mass1[i], width1[i], s23, s2, s3, pro2 );
573
574 if ( g0[i] == 2 )
575 {
576 double skm2[2] = { s1, mass_EtaP * mass_EtaP };
577 double spi2[2] = { s3, mass_Kaon * mass_Kaon };
578 propagatorKstr1430( mass1[i], s13, skm2, spi2, pro2 );
579 }
580 if ( g0[i] == 4 ) KPiSLASS( s13, s1, s3, pro2 );
581 if ( g0[i] == 0 )
582 {
583 pro2[0] = 1;
584 pro2[1] = 0;
585 }
586 amp_tmp[0] = temp_PDF1 * pro1[0] + temp_PDF2 * pro2[0];
587 amp_tmp[1] = temp_PDF1 * pro1[1] + temp_PDF2 * pro2[1];
588 }
589
590 if ( modetype[i] == 123 )
591 {
592 temp_PDF3 = DDalitz( Ks01, Kp, Ks02, spin[i], mass1[i] );
593 if ( g0[i] == 1 ) propagatorGS( mass1[i], width1[i], s13, s1, s3, 9.0, pro3 );
594 if ( g0[i] == 2 )
595 {
596 double skm2[2] = { s1, mass_EtaP * mass_EtaP };
597 double spi2[2] = { s2, mass_Kaon * mass_Kaon };
598 propagatorKstr1430( mass1[i], s12, skm2, spi2, pro3 );
599 }
600 if ( g0[i] == 4 ) KPiSLASS( s12, s1, s2, pro3 );
601 if ( g0[i] == 0 )
602 {
603 pro1[0] = 1;
604 pro1[1] = 0;
605 }
606
607 temp_PDF4 = DDalitz( Ks02, Kp, Ks01, spin[i], mass1[i] );
608 if ( g0[i] == 1 ) propagatorGS( mass1[i], width1[i], s23, s2, s3, 9.0, pro4 );
609
610 if ( g0[i] == 2 )
611 {
612 double skm2[2] = { s1, mass_EtaP * mass_EtaP };
613 double spi2[2] = { s3, mass_Kaon * mass_Kaon };
614 propagatorKstr1430( mass1[i], s13, skm2, spi2, pro2 );
615 }
616 if ( g0[i] == 4 ) KPiSLASS( s13, s1, s3, pro2 );
617 if ( g0[i] == 0 )
618 {
619 pro2[0] = 1;
620 pro2[1] = 0;
621 }
622 amp_tmp[0] = temp_PDF3 * pro3[0] + temp_PDF4 * pro4[0];
623 amp_tmp[1] = temp_PDF3 * pro3[1] + temp_PDF4 * pro4[1];
624 }
625
626 if ( modetype[i] == 122 )
627 {
628 amp_tmp[0] = 2;
629 amp_tmp[1] = 0;
630 }
631 // cout<<"pdf: "<<i<<", "<<amp_tmp[0]<<", "<<amp_tmp[1]<<endl;
632 Com_Multi( amp_tmp, cof, amp_PDF );
633 PDF[0] += amp_PDF[0];
634 PDF[1] += amp_PDF[1];
635
636 // cout<<"PDF: "<<i<<", "<<amp_PDF[0]<<", "<<amp_PDF[1]<<endl;
637 }
638 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
639 // cout<<"value: "<<value<<endl;
640 if ( value <= 0 ) value = 1e-20;
641
642 Result = value;
643}
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
a0
character *LEPTONflag integer iresonances real zeta5 real a0
Definition Ekhara/src/ekhara/common.h:4
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
one
const DifNumber one
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
CS
TCrossPart * CS
Definition Mcgpj.cxx:51
EvtDToKSKSK::initProbMax
void initProbMax()
Definition EvtDToKSKSK.cc:91
EvtDToKSKSK::getName
void getName(std::string &name)
Definition EvtDToKSKSK.cc:33
EvtDToKSKSK::~EvtDToKSKSK
virtual ~EvtDToKSKSK()
Definition EvtDToKSKSK.cc:31
EvtDToKSKSK::clone
EvtDecayBase * clone()
Definition EvtDToKSKSK.cc:35
EvtDToKSKSK::decay
void decay(EvtParticle *p)
Definition EvtDToKSKSK.cc:93
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
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