BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToKSKpPipPim.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: EvtDsToKSKpPipPim.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 "EvtDsToKSKpPipPim.hh"
18#include "../EvtGenBase/EvtConst.hh"
19#include "../EvtGenBase/EvtDecayTable.hh"
20#include "../EvtGenBase/EvtFlatte.hh"
21#include "../EvtGenBase/EvtGenKine.hh"
22#include "../EvtGenBase/EvtPDL.hh"
23#include "../EvtGenBase/EvtParticle.hh"
24#include "../EvtGenBase/EvtPatches.hh"
25#include "../EvtGenBase/EvtReport.hh"
26#include "../EvtGenBase/EvtResonance.hh"
27#include "../EvtGenBase/EvtResonance2.hh"
28#include <fstream>
29#include <stdlib.h>
30#include <string>
31using std::endl;
32
33EvtDsToKSKpPipPim::~EvtDsToKSKpPipPim() {}
34
35void EvtDsToKSKpPipPim::getName( std::string& model_name ) { model_name = "DsToKSKpPipPim"; }
36
37EvtDecayBase* EvtDsToKSKpPipPim::clone() { return new EvtDsToKSKpPipPim; }
38
39void EvtDsToKSKpPipPim::init() {
40 // check that there are 0 arguments
41 checkNArg( 0 );
42 checkNDaug( 4 );
43
44 checkSpinParent( EvtSpinType::SCALAR );
45 checkSpinDaughter( 0, EvtSpinType::SCALAR );
46 checkSpinDaughter( 1, EvtSpinType::SCALAR );
47 checkSpinDaughter( 2, EvtSpinType::SCALAR );
48 checkSpinDaughter( 3, EvtSpinType::SCALAR );
49
50 int mode = 0;
51 phi[1] = -3.1173e-03;
52 phi[2] = 5.6119e+00;
53 phi[3] = -7.0591e-01;
54 phi[4] = -3.1220e+00;
55 phi[5] = 3.5914e+00;
56 phi[6] = 3.8510e+00;
57 phi[7] = -1.0830e+00;
58
59 rho[1] = 1.8516e+00;
60 rho[2] = 2.8585e+00;
61 rho[3] = 1.1065e+00;
62 rho[4] = 5.4084e-01;
63 rho[5] = 1.3816e+00;
64 rho[6] = 1.9177e+00;
65 rho[7] = 2.2119e+00;
66 modetype[0] = 7;
67 modetype[1] = 4;
68 modetype[2] = 3;
69 modetype[3] = 15;
70 modetype[4] = 8;
71 modetype[5] = 11;
72 modetype[6] = 12;
73 modetype[7] = 17;
74 phi[0] = 0;
75 rho[0] = 1;
76
77 // cout << "DsToKSKpPipPim :" << endl;
78 // for (int i=0; i<8; i++) {
79 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
80 // }
81
82 mD = 1.86486;
83 metap = 0.95778;
84 mk0 = 0.497614;
85 mass_Kaon = 0.49368;
86 mass_Eta = 0.547862;
87 mass_Pion = 0.13957;
88 math_pi = 3.1415926;
89 mass_Pion2 = 0.0194797849;
90 mass_2Pion = 0.27914;
91 math_2pi = 6.2831852;
92
93 rD2 = 25.0;
94 rRes2 = 9.0;
95 GK1270 = 8.7000e-02;
96 GK1400 = 1.7400e-01;
97 GKst0 = 4.7300e-02;
98 GKstp = 5.0300e-02;
99 GP1680 = 1.5000e-01;
100 GRho = 1.4780e-01;
101 Ga0_980 = 5.0000e-02;
102 Geta1475 = 9.0000e-02;
103 mK1270 = 1.2720e+00;
104 mK1400 = 1.4030e+00;
105 mKst0 = 8.9555e-01;
106 mKstp = 8.9176e-01;
107 mP1680 = 1.6800e+00;
108 mRho = 7.7526e-01;
109 ma0_980 = 9.9000e-01;
110 meta1475 = 1.4750e+00;
111
112 GS1 = 0.636619783;
113 GS2 = 0.01860182466;
114 GS3 = 0.1591549458; // 1/(2*math_2pi)
115 GS4 = 0.00620060822; // mass_Pion2/math_pi
116
117 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
118 int EE[4][4][4][4] = {
119 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
120 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, -1, 0 } },
121 { { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 1, 0, 0 } },
122 { { 0, 0, 0, 0 }, { 0, 0, 1, 0 }, { 0, -1, 0, 0 }, { 0, 0, 0, 0 } } },
123 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 1, 0 } },
124 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
125 { { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 } },
126 { { 0, 0, -1, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 } } },
127 { { { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, -1, 0, 0 } },
128 { { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 } },
129 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
130 { { 0, 1, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } },
131 { { { 0, 0, 0, 0 }, { 0, 0, -1, 0 }, { 0, 1, 0, 0 }, { 0, 0, 0, 0 } },
132 { { 0, 0, 1, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 } },
133 { { 0, -1, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
134 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } } };
135 for ( int i = 0; i < 4; i++ )
136 {
137 for ( int j = 0; j < 4; j++ )
138 {
139 G[i][j] = GG[i][j];
140 for ( int k = 0; k < 4; k++ )
141 {
142 for ( int l = 0; l < 4; l++ ) { E[i][j][k][l] = EE[i][j][k][l]; }
143 }
144 }
145 }
146}
147
148void EvtDsToKSKpPipPim::initProbMax() { setProbMax( 824000 ); }
149
150void EvtDsToKSKpPipPim::decay( EvtParticle* p ) {
151 /*
152 double maxprob=0;
153 for(int ir=0;ir<=60000000;ir++){
154 p->initializePhaseSpace(getNDaug(),getDaugs());
155 EvtVector4R _ks = p->getDaug(0)->getP4();
156 EvtVector4R _kp = p->getDaug(1)->getP4();
157 EvtVector4R _pip = p->getDaug(2)->getP4();
158 EvtVector4R _pim = p->getDaug(3)->getP4();
159
160 double _Pip[4],_Pim[4],_Ks[4],_Kp[4];
161 _Ks[0] = _ks.get(0); _Kp[0] = _kp.get(0); _Pip[0] = _pip.get(0); _Pim[0] =
162 _pim.get(0); _Ks[1] = _ks.get(1); _Kp[1] = _kp.get(1); _Pip[1] = _pip.get(1); _Pim[1] =
163 _pim.get(1); _Ks[2] = _ks.get(2); _Kp[2] = _kp.get(2); _Pip[2] = _pip.get(2); _Pim[2] =
164 _pim.get(2); _Ks[3] = _ks.get(3); _Kp[3] = _kp.get(3); _Pip[3] = _pip.get(3); _Pim[3] =
165 _pim.get(3);
166
167 double _prob;
168 int g0[8]={1,1,1,1,1,1,1,1};
169 int nstates=10;
170
171 calEvaMy(_Ks, _Kp, _Pip, _Pim, mass, width, rho, phi, g0, modetype, nstates,
172 _prob); if(_prob>maxprob) { maxprob=_prob; printf("ir = %d,maxprob = %.10f,prob =
173 %.10f\n\n",ir,maxprob,_prob);
174 }
175 }
176 printf("maxprob = %.10f\n", maxprob);
177 */
178 p->initializePhaseSpace( getNDaug(), getDaugs() );
179 EvtVector4R ks = p->getDaug( 0 )->getP4();
180 EvtVector4R kp = p->getDaug( 1 )->getP4();
181 EvtVector4R pip = p->getDaug( 2 )->getP4();
182 EvtVector4R pim = p->getDaug( 3 )->getP4();
183
184 double Pip[4], Pim[4], Kp[4], Ks[4];
185 Ks[0] = ks.get( 0 );
186 Kp[0] = kp.get( 0 );
187 Pip[0] = pip.get( 0 );
188 Pim[0] = pim.get( 0 );
189 Ks[1] = ks.get( 1 );
190 Kp[1] = kp.get( 1 );
191 Pip[1] = pip.get( 1 );
192 Pim[1] = pim.get( 1 );
193 Ks[2] = ks.get( 2 );
194 Kp[2] = kp.get( 2 );
195 Pip[2] = pip.get( 2 );
196 Pim[2] = pim.get( 2 );
197 Ks[3] = ks.get( 3 );
198 Kp[3] = kp.get( 3 );
199 Pip[3] = pip.get( 3 );
200 Pim[3] = pim.get( 3 );
201
202 // Ks[0] = 0.6096125554; Ks[1] = 0.2036145752; Ks[2] = -0.2256821642; Ks[3] =
203 // 0.1884652565; Kp[0] = 0.7106732604; Kp[1] = -0.2599714599; Kp[2] = 0.3643877399;
204 // Kp[3] = 0.2469330230; Pip[0] = 0.4177584141; Pip[1] = -0.0562004660; Pip[2] =
205 // 0.2821155804; Pip[3] =-0.0020089709; Pim[0] = 0.3945120988; Pim[1] = -0.1225632042;
206 // Pim[2] = -0.3645928242; Pim[3] = 0.0521817527;
207
208 // Ks[0] = 0.6439438863; Ks[1] = 0.3368405824; Ks[2] = 0.2308441601; Ks[3] = 0.0171298406;
209 // Kp[0] = 0.7146391007; Kp[1] = 0.4061832461; Kp[2] = -0.1841669581; Kp[3] = 0.2609401579;
210 // Pip[0] = 0.2817222681; Pip[1] = -0.1005161765;Pip[2] = -0.1172908044;Pip[3] =
211 // -0.1898077097; Pim[0] = 0.3638598699;Pim[1] = -0.3106229764;Pim[2] = 0.0905020525;Pim[3] =
212 // 0.0907574503;
213
214 double value;
215 int g0[8] = { 1, 1, 1, 1, 1, 1, 1, 1 };
216 int nstates = 8;
217
218 calEvaMy( Ks, Kp, Pip, Pim, mass, width, rho, phi, g0, modetype, nstates, value );
219 setProb( value );
220
221 return;
222}
223
224void EvtDsToKSKpPipPim::Com_Multi( double a1[2], double a2[2], double res[2] ) {
225 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
226 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
227}
228void EvtDsToKSKpPipPim::Com_Divide( double a1[2], double a2[2], double res[2] ) {
229 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / ( a2[0] * a2[0] + a2[1] * a2[1] );
230 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / ( a2[0] * a2[0] + a2[1] * a2[1] );
231}
232double EvtDsToKSKpPipPim::SCADot( double a1[4], double a2[4] ) {
233 double _cal = 0.;
234 _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
235 return _cal;
236}
237double EvtDsToKSKpPipPim::barrier( int l, double sa, double sb, double sc, double r2 ) {
238 double F;
239 double tmp = sa + sb - sc;
240 double q = 0.25 * tmp * tmp / sa - sb;
241 if ( q < 0 ) q = 1e-16;
242 double z = q * r2;
243 if ( l == 1 ) { F = sqrt( 2.0 * z / ( 1.0 + z ) ); }
244 else if ( l == 2 )
245 {
246 double z2 = z * z;
247 F = sqrt( 13.0 * z2 / ( 9.0 + 3.0 * z + z2 ) );
248 }
249 else { F = 1.0; }
250 return F;
251}
252void EvtDsToKSKpPipPim::calt1( double daug1[4], double daug0[4], double t1[4] ) {
253 double p, pq;
254 double pa[4], qa[4];
255 for ( int i = 0; i < 4; i++ )
256 {
257 pa[i] = daug1[i] + daug0[i];
258 qa[i] = daug1[i] - daug0[i];
259 }
260 p = SCADot( pa, pa );
261 pq = SCADot( pa, qa );
262 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - pq / p * pa[i]; }
263}
264void EvtDsToKSKpPipPim::calt2( double daug1[4], double daug0[4], double t2[4][4] ) {
265 double p, r;
266 double pa[4], t1[4];
267 calt1( daug1, daug0, t1 );
268 r = SCADot( t1, t1 ) / 3.0;
269 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug0[i]; }
270 p = SCADot( pa, pa );
271 for ( int i = 0; i < 4; i++ )
272 {
273 for ( int j = 0; j < 4; j++ )
274 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
275 }
276}
277
278double EvtDsToKSKpPipPim::wid( double mass, double sa, double sb, double sc, double r2,
279 int l ) {
280 double widm = 0.;
281 double sa0 = mass * mass;
282 double m = sqrt( sa );
283 double tmp = sb - sc;
284 double tmp1 = sa + tmp;
285 double q = 0.25 * tmp1 * tmp1 / sa - sb;
286 if ( q < 0 ) q = 1e-16;
287 double tmp2 = sa0 + tmp;
288 double q0 = 0.25 * tmp2 * tmp2 / sa0 - sb;
289 if ( q0 < 0 ) q0 = 1e-16;
290 double z = q * r2;
291 double z0 = q0 * r2;
292 double t = q / q0;
293 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
294 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
295 else if ( l == 2 )
296 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
297 return widm;
298}
299double EvtDsToKSKpPipPim::widl1( double mass, double sa, double sb, double sc, double r2 ) {
300 double widm = 0.;
301 double sa0 = mass * mass;
302 double m = sqrt( sa );
303 double tmp = sb - sc;
304 double tmp1 = sa + tmp;
305 double q = 0.25 * tmp1 * tmp1 / sa - sb;
306 if ( q < 0 ) q = 1e-16;
307 double tmp2 = sa0 + tmp;
308 double q0 = 0.25 * tmp2 * tmp2 / sa0 - sb;
309 if ( q0 < 0 ) q0 = 1e-16;
310 double z = q * r2;
311 double z0 = q0 * r2;
312 double F = ( 1 + z0 ) / ( 1 + z );
313 double t = q / q0;
314 widm = t * sqrt( t ) * mass / m * F;
315 return widm;
316}
317void EvtDsToKSKpPipPim::Flatte_rhoab( double sa, double sb, double sc, double rho[2] ) {
318 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
319 if ( q > 0 )
320 {
321 rho[0] = 2 * sqrt( q / sa );
322 rho[1] = 0;
323 }
324 else if ( q < 0 )
325 {
326 rho[0] = 0;
327 rho[1] = 2 * sqrt( -q / sa );
328 }
329}
330void EvtDsToKSKpPipPim::propagatora0980( double mass, double sx, double* sb, double* sc,
331 double prop[2] ) {
332 double unit[2] = { 1.0 };
333 double ci[2] = { 0, 1 };
334 double rho1[2];
335 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
336 double rho2[2];
337 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
338
339 double gKK_a980 = 0.892 * 0.341, gPiEta_a980 = 0.341;
340 double tmp1[2] = { gKK_a980, 0 };
341 double tmp11[2];
342 double tmp2[2] = { gPiEta_a980, 0 };
343 double tmp22[2];
344 Com_Multi( tmp1, rho1, tmp11 );
345 Com_Multi( tmp2, rho2, tmp22 );
346 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
347 double tmp31[2];
348 Com_Multi( tmp3, ci, tmp31 );
349 double tmp4[2] = { mass * mass - sx - tmp31[0], -1.0 * tmp31[1] };
350 Com_Divide( unit, tmp4, prop );
351}
352void EvtDsToKSKpPipPim::propagatorRBW( double mass, double width, double sa, double sb,
353 double sc, double r2, int l, double prop[2] ) {
354 double a[2], b[2];
355 a[0] = 1;
356 a[1] = 0;
357 b[0] = mass * mass - sa;
358 b[1] = -mass * width * wid( mass, sa, sb, sc, r2, l );
359 Com_Divide( a, b, prop );
360}
361
362void EvtDsToKSKpPipPim::KPiSLASS( double sa, double sb, double sc, double prop[2] ) {
363 const double m1430 = 1.441;
364 const double sa0 = 2.076481; // m1410*m1410;
365 const double w1430 = 0.193;
366 const double Lass1 = 0.25 / sa0;
367 double tmp = sb - sc;
368 double tmp1 = sa0 + tmp;
369 double q0 = Lass1 * tmp1 * tmp1 - sb;
370 if ( q0 < 0 ) q0 = 1e-16;
371 double tmp2 = sa + tmp;
372 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
373 double q = sqrt( qs );
374 double width = w1430 * q * m1430 / sqrt( sa * q0 );
375 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
376 if ( temp_R < 0 ) temp_R += math_pi;
377 double deltaR = -109.7 + temp_R; // fiR=-109.7
378 double temp_F =
379 atan( 0.226 * q / ( 2.0 - 3.8194 * qs ) ); // 2.0*0.113=0.226; 0.113*33.8=1.926; a=0.113;
380 // r=-33.8
381 if ( temp_F < 0 ) temp_F += math_pi;
382 double deltaF = 0.1 + temp_F; // fiF=0.1
383 double deltaS = deltaR + 2.0 * deltaF; ///
384 double t1 = 0.96 * sin( deltaF ); // F= 0.96
385 double t2 = sin( deltaR );
386 double CF[2], CS[2];
387 CF[0] = cos( deltaF );
388 CF[1] = sin( deltaF );
389 CS[0] = cos( deltaS );
390 CS[1] = sin( deltaS );
391 prop[0] = t1 * CF[0] + t2 * CS[0];
392 prop[1] = t1 * CF[1] + t2 * CS[1];
393}
394
395void EvtDsToKSKpPipPim::propagatorGS( double mass2, double mass, double width, double sa,
396 double sb, double sc, double r2, double prop[2] ) {
397 double a[2], b[2];
398 double tmp = sb - sc;
399 double tmp1 = sa + tmp;
400 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
401 if ( q2 < 0 ) q2 = 1e-16;
402
403 double tmp2 = mass2 + tmp;
404 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
405 if ( q02 < 0 ) q02 = 1e-16;
406
407 double q = sqrt( q2 );
408 double q0 = sqrt( q02 );
409 double m = sqrt( sa );
410 double q03 = q0 * q02;
411 double tmp3 = log( mass + 2 * q0 ) + 1.2760418309;
412
413 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2760418309 );
414 double h0 = GS1 * q0 / mass * tmp3;
415 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
416 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
417 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
418
419 a[0] = 1.0 + d * width / mass;
420 a[1] = 0.0;
421 b[0] = mass2 - sa + width * f;
422 b[1] = -mass * width * widl1( mass, sa, sb, sc, r2 );
423 Com_Divide( a, b, prop );
424}
425
426void EvtDsToKSKpPipPim::calEvaMy( double* Ks, double* Kp, double* Pip, double* Pim,
427 double* mass1, double* width1, double* amp, double* phase,
428 int* g0, int* modetype, int nstates, double& Result ) {
429 double pD[4], pKpPim[4], pKsPim[4], pKsKpPim[4], pKsKp[4], pKsPipPim[4], pKpPipPim[4],
430 pPipPim[4];
431 for ( int i = 0; i != 4; i++ )
432 {
433
434 pD[i] = Ks[i] + Pip[i] + Kp[i] + Pim[i];
435
436 pKsPim[i] = Ks[i] + Pim[i];
437 pKpPim[i] = Kp[i] + Pim[i];
438 pPipPim[i] = Pip[i] + Pim[i];
439 pKsKp[i] = Kp[i] + Ks[i];
440 pKsKpPim[i] = pKsPim[i] + Kp[i];
441 pKpPipPim[i] = Kp[i] + Pip[i] + Pim[i];
442 pKsPipPim[i] = Ks[i] + Pip[i] + Pim[i];
443 }
444
445 double sPim, sPip, sKs, sKp, sKpPim, sD, sKsPim, sKsKpPim, sKsKp, sKpPipPim, sKsPipPim,
446 sPipPim;
447
448 sD = SCADot( pD, pD );
449 sKs = SCADot( Ks, Ks );
450 sKp = SCADot( Kp, Kp );
451 sPip = SCADot( Pip, Pip );
452 sPim = SCADot( Pim, Pim );
453 sPipPim = SCADot( pPipPim, pPipPim );
454 sKsPim = SCADot( pKsPim, pKsPim );
455 sKpPim = SCADot( pKpPim, pKpPim );
456 sKsKpPim = SCADot( pKsKpPim, pKsKpPim );
457 sKpPipPim = SCADot( pKpPipPim, pKpPipPim );
458 sKsPipPim = SCADot( pKsPipPim, pKsPipPim );
459 sKsKp = SCADot( pKsKp, pKsKp );
460
461 double t1D[4], t1KpPim[4], t1KsPim[4], t1KsKpPim_kspim[4], t1KsKpPip_ksk[4],
462 t1KsKp_PipPim[4], t1PipPim[4], t1KsKpPim_ksk[4], t1KsKpPim_kpim[4], t1KsKp[4];
463 double t2KsKpPim_kspim[4][4], t2KsKpPim_kpim[4][4], t2KpPipPim_kpim[4][4],
464 t2KsPipPim_kspim[4][4], t2KpPipPim_pippim[4][4], t2KsPipPim_pippim[4][4];
465
466 calt1( Ks, Pim, t1KsPim );
467 calt1( Kp, Pim, t1KpPim );
468 calt1( Ks, Kp, t1KsKp );
469 calt1( Pip, Pim, t1PipPim );
470 calt1( pKsPim, Kp, t1KsKpPim_kspim );
471 calt1( pKpPim, Ks, t1KsKpPim_kpim );
472 calt1( pKsKp, Pip, t1KsKpPip_ksk );
473 calt1( pKsKp, Pim, t1KsKpPim_ksk );
474 calt1( pKsKp, pPipPim, t1KsKp_PipPim );
475 calt1( pKsKpPim, Pip, t1D );
476
477 calt2( pKsPim, Kp, t2KsKpPim_kspim );
478 calt2( pKpPim, Ks, t2KsKpPim_kpim );
479 calt2( pKpPim, Pip, t2KpPipPim_kpim );
480 calt2( pKsPim, Pip, t2KsPipPim_kspim );
481 calt2( pPipPim, Ks, t2KsPipPim_pippim );
482 calt2( pPipPim, Kp, t2KpPipPim_pippim );
483
484 double cof[2], amp_tmp[2], amp_PDF[2], PDF[2];
485 amp_PDF[0] = 0;
486 amp_PDF[1] = 0;
487 PDF[0] = 0;
488 PDF[1] = 0;
489 double temp_PDF, tmp1;
490 double pro[2], pro1[2], pro_kspim[2], pro_kpim[2], pro_ksk[2], pro_pippim[2],
491 pro_kskpim_kspim[2], pro_kskpim_kpim[2];
492 double ispro_kspim = 0, ispro_kpim = 0, ispro_pippim = 0;
493 double barr_kskpim_kspim = -1.0, barr_kskpim_kpim = -1.0, barr_kspim = -1.0,
494 barr_kpim = -1.0, barr_pippim = -1.0, barr_kskp_pippim = -1, barr_ds_kskpim = -1,
495 barr_kspippim_kspim2 = -1, barr_kpippim_kpim2 = -1, barr_kspippim_pippim2 = -1,
496 barr_kpippim_pippim2 = -1, barr_ds_kspippim = -1, barr_ds_kpippim = -1, barr_ksk = -1,
497 barr_kspippim_kspim = -1, barr_kpippim_kpim = -1;
498
499 double sKs2[2] = { sKs, mass_Pion * mass_Pion };
500 double sKp2[2] = { sKp, mass_Eta * mass_Eta };
501 double mass2sq, tt1, tt2, t2D[4][4];
502
503 for ( int i = 0; i < nstates; i++ )
504 {
505 temp_PDF = 0;
506 amp_tmp[0] = 0;
507 amp_tmp[1] = 0;
508 tmp1 = 0;
509 cof[0] = amp[i] * cos( phase[i] );
510 cof[1] = amp[i] * sin( phase[i] );
511 mass2sq = mRho * mRho;
512
513 // printf("phase= %i, %.5f, %.5f\n", i, amp[i],phase[i]);
514
515 // Ds -> eta(1475)pi+ eta(1475)->K*Kp,K*+>Ks pi+
516 if ( modetype[i] == 3 )
517 {
518 propagatorRBW( meta1475, Geta1475, sKsKpPim, sKsPim, sKp, rRes2, 1, pro_kskpim_kspim );
519 propagatorRBW( mKstp, GKstp, sKsPim, sKs, sPim, rRes2, 1, pro_kspim );
520 ispro_kspim = 1;
521
522 Com_Multi( pro_kskpim_kspim, pro_kspim, pro );
523
524 barr_kspim = barrier( 1, sKsPim, sKs, sPim, rRes2 );
525 barr_kskpim_kspim = barrier( 1, sKsKpPim, sKsPim, sKp, rRes2 );
526 for ( int a = 0; a < 4; a++ ) { temp_PDF += Kp[a] * t1KsPim[a] * G[a][a]; }
527 tmp1 = barr_kspim * barr_kskpim_kspim * temp_PDF;
528 }
529
530 // Ds -> eta(1475)pi+ eta(1475)->K*Ks,K*0>Kp pi+
531 if ( modetype[i] == 4 )
532 {
533 propagatorRBW( meta1475, Geta1475, sKsKpPim, sKpPim, sKs, rRes2, 1, pro_kskpim_kpim );
534 propagatorRBW( mKst0, GKst0, sKpPim, sKp, sPim, rRes2, 1, pro_kpim );
535 ispro_kpim = 1;
536 Com_Multi( pro_kskpim_kpim, pro_kpim, pro );
537
538 barr_kpim = barrier( 1, sKpPim, sKp, sPim, rRes2 );
539 barr_kskpim_kpim = barrier( 1, sKsKpPim, sKpPim, sKs, rRes2 );
540 for ( int a = 0; a < 4; a++ ) { temp_PDF += Ks[a] * t1KpPim[a] * G[a][a]; }
541 tmp1 = barr_kpim * barr_kskpim_kpim * temp_PDF;
542 }
543
544 // 0- eta1475 -> a0ipi+
545 if ( modetype[i] == 15 )
546 {
547 propagatorRBW( meta1475, Geta1475, sKsKpPim, sKsKp, sPim, rRes2, 0, pro1 );
548 propagatora0980( ma0_980, sKsKp, sKs2, sKp2, pro_ksk );
549 Com_Multi( pro1, pro_ksk, pro );
550 tmp1 = 1.0;
551 }
552
553 // 1+ Ds -> K1270 K+ , K1270 -> Ks rho rho -> pi+ pi-
554 if ( modetype[i] == 7 )
555 {
556 propagatorRBW( mK1270, GK1270, sKsPipPim, sPipPim, sKs, rRes2, 0, pro1 );
557 if ( ispro_pippim == 0 )
558 {
559 propagatorGS( mass2sq, mRho, GRho, sPipPim, sPip, sPim, rRes2, pro_pippim );
560 ispro_pippim = 1;
561 }
562 Com_Multi( pro1, pro_pippim, pro );
563
564 if ( barr_pippim < 0.0 ) barr_pippim = barrier( 1, sPipPim, sPip, sPim, rRes2 );
565 if ( barr_ds_kspippim < 0 ) barr_ds_kspippim = barrier( 1, sD, sKsPipPim, sKp, rD2 );
566 for ( int a = 0; a < 4; a++ )
567 {
568 for ( int j = 0; j < 4; j++ )
569 {
570 temp_PDF += t1D[a] * ( pKsPipPim[a] * pKsPipPim[j] / sKsPipPim - G[a][j] ) *
571 t1PipPim[j] * G[a][a] * G[j][j];
572 }
573 }
574 tmp1 = barr_pippim * barr_ds_kspippim * temp_PDF;
575 }
576
577 // 1+ Ds -> K1270 Kp , K1270 -> K+ rho rho -> pi+ pi-
578 if ( modetype[i] == 8 )
579 {
580 propagatorRBW( mK1270, GK1270, sKpPipPim, sPipPim, sKp, rRes2, 0, pro1 );
581 if ( ispro_pippim == 0 )
582 {
583 propagatorGS( mass2sq, mRho, GRho, sPipPim, sPip, sPim, rRes2, pro_pippim );
584 ispro_pippim = 1;
585 }
586 Com_Multi( pro1, pro_pippim, pro );
587
588 if ( barr_pippim < 0.0 ) barr_pippim = barrier( 1, sPipPim, sPip, sPim, rRes2 );
589 if ( barr_ds_kpippim < 0 ) barr_ds_kpippim = barrier( 1, sD, sKpPipPim, sKs, rD2 );
590 for ( int a = 0; a < 4; a++ )
591 {
592 for ( int j = 0; j < 4; j++ )
593 {
594 temp_PDF += t1D[a] * ( pKpPipPim[a] * pKpPipPim[j] / sKpPipPim - G[a][j] ) *
595 t1PipPim[j] * G[a][a] * G[j][j];
596 }
597 }
598 tmp1 = barr_pippim * barr_ds_kpippim * temp_PDF;
599 }
600
601 // 1+ Ds -> K1400 K+ , K1400 -> pi+ K*- ,K*- ->Ks pi-
602 if ( modetype[i] == 11 )
603 {
604 propagatorRBW( mK1400, GK1400, sKsPipPim, sKsPim, sPip, rRes2, 0, pro1 );
605 if ( ispro_kspim == 0 )
606 {
607 propagatorRBW( mKstp, GKstp, sKsPim, sKs, sPim, rRes2, 1, pro_kspim );
608 ispro_kspim = 1;
609 }
610 Com_Multi( pro1, pro_kspim, pro );
611
612 if ( barr_kspim < 0.0 ) barr_kspim = barrier( 1, sKsPim, sKs, sPim, rRes2 );
613 if ( barr_ds_kspippim < 0.0 ) barr_ds_kspippim = barrier( 1, sD, sKsPipPim, sKp, rD2 );
614 for ( int a = 0; a < 4; a++ )
615 {
616 for ( int j = 0; j < 4; j++ )
617 {
618 temp_PDF += t1D[a] * ( pKsPipPim[a] * pKsPipPim[j] / sKsPipPim - G[a][j] ) *
619 t1KsPim[j] * G[a][a] * G[j][j];
620 }
621 }
622 tmp1 = barr_kspim * barr_ds_kspippim * temp_PDF;
623 }
624
625 // 1+ Ds -> K1400 Ks , K1400 -> pi+ K*0 ,K*0 ->K+ pi-
626 if ( modetype[i] == 12 )
627 {
628 propagatorRBW( mK1400, GK1400, sKpPipPim, sKpPim, sPip, rRes2, 0, pro1 );
629 if ( ispro_kpim == 0 )
630 {
631 propagatorRBW( mKst0, GKst0, sKpPim, sKp, sPim, rRes2, 1, pro_kpim );
632 ispro_kpim = 1;
633 }
634 Com_Multi( pro1, pro_kpim, pro );
635
636 if ( barr_kpim < 0.0 ) barr_kpim = barrier( 1, sKpPim, sKp, sPim, rRes2 );
637 if ( barr_ds_kpippim < 0.0 ) barr_ds_kpippim = barrier( 1, sD, sKpPipPim, sKs, rD2 );
638 for ( int a = 0; a < 4; a++ )
639 {
640 for ( int j = 0; j < 4; j++ )
641 {
642 temp_PDF += t1D[a] * ( pKpPipPim[a] * pKpPipPim[j] / sKpPipPim - G[a][j] ) *
643 t1KpPim[j] * G[a][a] * G[j][j];
644 }
645 }
646 tmp1 = barr_kpim * barr_ds_kpippim * temp_PDF;
647
648 // printf("sD:%.15f,%.15f,%.15f,%.15f\n ",sD,sKpPipPim,sKs,rD2);
649 // printf("aaaaa= %.15f,%.15f,%.15f,%.15f\n", barr_kpim,barr_ds_kpippim,temp_PDF,tmp1);
650 }
651
652 // 1+ Ds -> phi1680 phi1680->kspim
653 if ( modetype[i] == 17 )
654 {
655
656 propagatorRBW( mP1680, GP1680, sKsKpPim, sKsPim, sKp, rRes2, 1, pro1 );
657 propagatorRBW( mKstp, GKstp, sKsPim, sKs, sPim, rRes2, 1, pro_kpim );
658 Com_Multi( pro1, pro_kpim, pro );
659
660 barr_kspim = barrier( 1, sKsPim, sKs, sPim, rRes2 );
661 barr_kskpim_kspim = barrier( 1, sKsKpPim, sKsPim, sKp, rRes2 );
662 barr_ds_kskpim = barrier( 1, sD, sKsKpPim, sPip, rD2 );
663 for ( int a = 0; a < 4; a++ )
664 {
665 for ( int j = 0; j < 4; j++ )
666 {
667 for ( int k = 0; k < 4; k++ )
668 {
669 for ( int l = 0; l < 4; l++ )
670 {
671 temp_PDF += E[a][j][k][l] * pKsKpPim[a] * ( pKsPim[j] - Kp[j] ) * Pip[k] *
672 ( Ks[l] - Pim[l] ) * G[a][a] * G[j][j] * G[k][k] * G[l][l];
673 }
674 }
675 }
676 }
677 tmp1 = barr_kspim * barr_kskpim_kspim * barr_ds_kskpim * temp_PDF;
678 }
679
680 if ( modetype[i] == 3 || modetype[i] == 4 || modetype[i] == 5 || modetype[i] == 7 ||
681 modetype[i] == 8 || modetype[i] == 11 || modetype[i] == 12 || modetype[i] == 17 ||
682 modetype[i] == 15 )
683 {
684 amp_tmp[0] = tmp1 * pro[0];
685 amp_tmp[1] = tmp1 * pro[1];
686 }
687 /*
688 if(modetype[i]==18){
689 amp_tmp[0] = tmp1*pro[0];
690 amp_tmp[1] = tmp1*pro[1];
691
692 }*/
693
694 Com_Multi( amp_tmp, cof, amp_PDF );
695
696 // printf("pro[0],pro[1]:: %.15f, %.15f\n", pro[0],pro[1]);
697 // printf("tmp1:: %.15f\n", tmp1);
698 // printf("amp_tmp[0],amp_tmp[1]:: %.15f, %.15f\n", amp_tmp[0],amp_tmp[1]);
699
700 PDF[0] += amp_PDF[0];
701 PDF[1] += amp_PDF[1];
702
703 // printf("cof= %.15f,%.15f\n",cof[0],cof[1]);
704 // printf("PDF= %.15f,%.15f\n",PDF[0],PDF[1]);
705 // printf("----\n");
706 }
707 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
708 if ( value <= 0 ) value = 1e-20;
709 Result = value;
710 // printf("value = %.15f\n",value);
711}
mass
double mass
Definition CosmicGenerator.cxx:128
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
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
CS
TCrossPart * CS
Definition Mcgpj.cxx:51
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
EvtDsToKSKpPipPim::getName
void getName(std::string &name)
Definition EvtDsToKSKpPipPim.cc:35
EvtDsToKSKpPipPim::decay
void decay(EvtParticle *p)
Definition EvtDsToKSKpPipPim.cc:150
EvtDsToKSKpPipPim::clone
EvtDecayBase * clone()
Definition EvtDsToKSKpPipPim.cc:37
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