BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToKKpipipi.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: EvtDsToKKpipipi.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 "EvtDsToKKpipipi.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 <iomanip>
30#include <stdlib.h>
31#include <string>
32using std::endl;
33
34EvtDsToKKpipipi::~EvtDsToKKpipipi() {}
35
36void EvtDsToKKpipipi::getName( std::string& model_name ) { model_name = "DsToKKpipipi"; }
37
38EvtDecayBase* EvtDsToKKpipipi::clone() { return new EvtDsToKKpipipi; }
39
40void EvtDsToKKpipipi::init() {
41 // check that there are 0 arguments
42 checkNArg( 0 );
43 checkNDaug( 5 );
44
45 checkSpinParent( EvtSpinType::SCALAR );
46 checkSpinDaughter( 0, EvtSpinType::SCALAR );
47 checkSpinDaughter( 1, EvtSpinType::SCALAR );
48 checkSpinDaughter( 2, EvtSpinType::SCALAR );
49 checkSpinDaughter( 3, EvtSpinType::SCALAR );
50 checkSpinDaughter( 4, EvtSpinType::SCALAR );
51
52 int mode = 0;
53 //--------------------amp---------------
54 phi[1] = 1.4732;
55 phi[2] = -4.2901;
56 rho[1] = 0.23396;
57 rho[2] = 7.4744;
58 //--------------------mass---------------------
59 mass1[0] = 1.019461;
60 mass1[1] = 1.019461;
61 mass1[2] = 1.019461;
62 mass2[0] = 0.77526;
63 mass2[1] = 0.77526;
64 mass2[2] = 0.77526;
65 mass3[0] = 1.23;
66 mass3[1] = 1.23;
67 mass3[2] = 1.23;
68
69 width1[0] = 0.004249;
70 width1[1] = 0.004249;
71 width1[2] = 0.004249;
72 width2[0] = 0.1478;
73 width2[1] = 0.1478;
74 width2[2] = 0.1478;
75 width3[0] = 0.42;
76 width3[1] = 0.42;
77 width3[2] = 0.42;
78
79 modetype[0] = 1;
80 modetype[1] = 3;
81 modetype[2] = 13;
82
83 phi[0] = 0;
84 rho[0] = 1;
85
86 // cout << "DsToKKpipipi :" << endl;
87 // for (int i=0; i<3; i++) {
88 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
89 // }
90
91 mass_Pion = 0.13957;
92 mass_Pion_N = 0.134977;
93 mass_Eta = 0.547862;
94 math_pi = 3.1415926;
95 rD2 = 25.0; // 5*5
96 rRes2 = 9.0; // 3*3
97 GS1 = 0.636619783;
98 GS2 = 0.01860182466;
99 GS3 = 0.1591549458; // 1/(2*math_2pi)
100 GS4 = 0.00620060822; // mass_Pion2/math_pi
101
102 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
103 int EE[4][4][4][4] = {
104 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
105 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, -1, 0 } },
106 { { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 1, 0, 0 } },
107 { { 0, 0, 0, 0 }, { 0, 0, 1, 0 }, { 0, -1, 0, 0 }, { 0, 0, 0, 0 } } },
108 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 1, 0 } },
109 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
110 { { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 } },
111 { { 0, 0, -1, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 } } },
112 { { { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, -1, 0, 0 } },
113 { { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 } },
114 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
115 { { 0, 1, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } },
116 { { { 0, 0, 0, 0 }, { 0, 0, -1, 0 }, { 0, 1, 0, 0 }, { 0, 0, 0, 0 } },
117 { { 0, 0, 1, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 } },
118 { { 0, -1, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
119 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } } };
120 for ( int i = 0; i < 4; i++ )
121 {
122 for ( int j = 0; j < 4; j++ )
123 {
124 G[i][j] = GG[i][j];
125 for ( int k = 0; k < 4; k++ )
126 {
127 for ( int l = 0; l < 4; l++ ) { E[i][j][k][l] = EE[i][j][k][l]; }
128 }
129 }
130 }
131}
132
133void EvtDsToKKpipipi::initProbMax() { setProbMax( 585000.0 ); }
134
135void EvtDsToKKpipipi::decay( EvtParticle* p ) {
136 /*
137 double maxprob = 0.0;
138 for(int ir=0;ir<=60000000;ir++){
139 p->initializePhaseSpace(getNDaug(),getDaugs());
140 EvtVector4R _km = p->getDaug(0)->getP4();
141 EvtVector4R _kp = p->getDaug(1)->getP4();
142 EvtVector4R _pip1 = p->getDaug(2)->getP4();
143 EvtVector4R _pip2 = p->getDaug(3)->getP4();
144 EvtVector4R _pim = p->getDaug(4)->getP4();
145
146 double _Km[4], _Kp[4], _Pip1[4], _Pip2[4], _Pim[4];
147 _Km[0] = _km.get(0); _Kp[0] = _kp.get(0); _Pip1[0] = _pip1.get(0); _Pip2[0] =
148 _pip2.get(0); _Pim[0] = _pim.get(0); _Km[1] = _km.get(1); _Kp[1] = _kp.get(1); _Pip1[1] =
149 _pip1.get(1); _Pip2[1] = _pip2.get(1); _Pim[1] = _pim.get(1); _Km[2] = _km.get(2); _Kp[2]
150 = _kp.get(2); _Pip1[2] = _pip1.get(2); _Pip2[2] = _pip2.get(2); _Pim[2] = _pim.get(2);
151 _Km[3] = _km.get(3); _Kp[3] = _kp.get(3); _Pip1[3] = _pip1.get(3); _Pip2[3] =
152 _pip2.get(3); _Pim[3] = _pim.get(3);
153
154 double _prob;
155 int nstates=3;
156 int g0[3]={1,1,1};
157 int g1[3]={1,1,1};
158 int g2[3]={0,0,0};
159 int g3[3]={0,1,0};
160 int g4[3]={0,0,0};
161 calEvaMy(_Km,_Kp,_Pip1,_Pip2,_Pim,mass1,mass2,mass3,width1,width2,width3,rho,phi,modetype,nstates,_prob);
162 if(_prob>maxprob) {
163 maxprob=_prob;
164 cout << "Max PDF = " << ir << " " << _prob << endl;
165 }
166 }
167 cout << "Max!!!!!!!!!!! " << maxprob<< endl;
168 */
169 p->initializePhaseSpace( getNDaug(), getDaugs() );
170 EvtVector4R km = p->getDaug( 0 )->getP4();
171 EvtVector4R kp = p->getDaug( 1 )->getP4();
172 EvtVector4R pip1 = p->getDaug( 2 )->getP4();
173 EvtVector4R pip2 = p->getDaug( 3 )->getP4();
174 EvtVector4R pim = p->getDaug( 4 )->getP4();
175
176 double Km[4], Kp[4], Pip1[4], Pip2[4], Pim[4];
177 Km[0] = km.get( 0 );
178 Kp[0] = kp.get( 0 );
179 Pip1[0] = pip1.get( 0 );
180 Pip2[0] = pip2.get( 0 );
181 Pim[0] = pim.get( 0 );
182 Km[1] = km.get( 1 );
183 Kp[1] = kp.get( 1 );
184 Pip1[1] = pip1.get( 1 );
185 Pip2[1] = pip2.get( 1 );
186 Pim[1] = pim.get( 1 );
187 Km[2] = km.get( 2 );
188 Kp[2] = kp.get( 2 );
189 Pip1[2] = pip1.get( 2 );
190 Pip2[2] = pip2.get( 2 );
191 Pim[2] = pim.get( 2 );
192 Km[3] = km.get( 3 );
193 Kp[3] = kp.get( 3 );
194 Pip1[3] = pip1.get( 3 );
195 Pip2[3] = pip2.get( 3 );
196 Pim[3] = pim.get( 3 );
197
198 double prob;
199 int nstates = 3;
200 calEvaMy( Km, Kp, Pip1, Pip2, Pim, mass1, mass2, mass3, width1, width2, width3, rho, phi,
201 modetype, nstates, prob );
202
203 setProb( prob );
204
205 return;
206}
207
208void EvtDsToKKpipipi::Com_Multi( double a1[2], double a2[2], double res[2] ) {
209 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
210 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
211}
212void EvtDsToKKpipipi::Com_Divide( double a1[2], double a2[2], double res[2] ) {
213 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
214 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
215 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
216}
217double EvtDsToKKpipipi::SCADot( double a1[4], double a2[4] ) {
218 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
219 return _cal;
220}
221double EvtDsToKKpipipi::barrier( int l, double sa, double sb, double sc, double r,
222 double mass ) {
223 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
224 if ( q < 0 ) q = 1e-16;
225 double z;
226 z = q * r * r;
227 double sa0;
228 sa0 = mass * mass;
229 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
230 if ( q0 < 0 ) q0 = 1e-16;
231 double z0 = q0 * r * r;
232 double F = 0.0;
233 if ( l == 0 ) F = 1;
234 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
235 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
236 return F;
237}
238// double EvtDsToKKpipipi::barrier(int l, double sa, double sb, double sc, double r2){
239// double F;
240// double tmp = sa+sb-sc;
241// double q = 0.25*tmp*tmp/sa-sb;
242// if (q < 0) q = 1e-16;
243// double z = q*r2;
244// if (l == 1) {
245// F = sqrt(2.0*z/(1.0+z));
246// }
247// else if (l == 2) {
248// double z2 = z*z;
249// F = sqrt(13.0*z2/(9.0+3.0*z+z2));
250// } else {
251// F = 1.0;
252// }
253// return F;
254// }
255void EvtDsToKKpipipi::calt1( double daug1[4], double daug2[4], double t1[4] ) {
256 double p, pq, tmp;
257 double pa[4], qa[4];
258 for ( int i = 0; i < 4; i++ )
259 {
260 pa[i] = daug1[i] + daug2[i];
261 qa[i] = daug1[i] - daug2[i];
262 }
263 p = SCADot( pa, pa );
264 pq = SCADot( pa, qa );
265 tmp = pq / p;
266 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
267}
268void EvtDsToKKpipipi::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
269 double p, r;
270 double pa[4], t1[4];
271 calt1( daug1, daug2, t1 );
272 r = SCADot( t1, t1 ) / 3.0;
273 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
274 p = SCADot( pa, pa );
275 for ( int i = 0; i < 4; i++ )
276 {
277 for ( int j = 0; j < 4; j++ )
278 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
279 }
280}
281void EvtDsToKKpipipi::Flatte_rhoab( double sa, double sb, double sc, double rho[2] ) {
282 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
283 if ( q > 0 )
284 {
285 rho[0] = 2 * sqrt( q / sa );
286 rho[1] = 0;
287 }
288 else if ( q < 0 )
289 {
290 rho[0] = 0;
291 rho[1] = 2 * sqrt( -q / sa );
292 }
293}
294void EvtDsToKKpipipi::propagatora0980( double mass2, double mass, double sx, double* sb,
295 double* sc, double prop[2] ) {
296 double unit[2] = { 1.0 };
297 double ci[2] = { 0, 1 };
298 double rho1[2];
299 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
300 double rho2[2];
301 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
302
303 double gKK_a980 = 0.892 * 0.341, gPiEta_a980 = 0.341;
304 double tmp1[2] = { gKK_a980, 0 };
305 double tmp11[2];
306 double tmp2[2] = { gPiEta_a980, 0 };
307 double tmp22[2];
308 Com_Multi( tmp1, rho1, tmp11 );
309 Com_Multi( tmp2, rho2, tmp22 );
310 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
311 double tmp31[2];
312 Com_Multi( tmp3, ci, tmp31 );
313 double tmp4[2] = { mass2 - sx - tmp31[0], -1.0 * tmp31[1] };
314 Com_Divide( unit, tmp4, prop );
315}
316double EvtDsToKKpipipi::wid( double mass2, double mass, double sa, double sb, double sc,
317 double r2, int l ) {
318 double widm = 0.;
319 double m = sqrt( sa );
320 double tmp = sb - sc;
321 double tmp1 = sa + tmp;
322 double q = 0.25 * tmp1 * tmp1 / sa - sb;
323 if ( q < 0 ) q = 1e-16;
324 double tmp2 = mass2 + tmp;
325 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
326 if ( q0 < 0 ) q0 = 1e-16;
327 double z = q * r2;
328 double z0 = q0 * r2;
329 double t = q / q0;
330 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
331 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
332 else if ( l == 2 )
333 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
334 return widm;
335}
336double EvtDsToKKpipipi::widl1( double mass2, double mass, double sa, double sb, double sc,
337 double r2 ) {
338 double widm = 0.;
339 double m = sqrt( sa );
340 double tmp = sb - sc;
341 double tmp1 = sa + tmp;
342 double q = 0.25 * tmp1 * tmp1 / sa - sb;
343 if ( q < 0 ) q = 1e-16;
344 double tmp2 = mass2 + tmp;
345 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
346 if ( q0 < 0 ) q0 = 1e-16;
347 double z = q * r2;
348 double z0 = q0 * r2;
349 double F = ( 1 + z0 ) / ( 1 + z );
350 double t = q / q0;
351 widm = t * sqrt( t ) * mass / m * F;
352 return widm;
353}
354void EvtDsToKKpipipi::propagatorRBW( double mass2, double mass, double width, double sa,
355 double sb, double sc, double r2, int l, double prop[2] ) {
356 double a[2], b[2];
357 a[0] = 1;
358 a[1] = 0;
359 b[0] = mass2 - sa;
360 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
361 Com_Divide( a, b, prop );
362}
363void EvtDsToKKpipipi::propagatorNBW( double mass2, double mass, double width, double sa,
364 double sb, double sc, double r2, int l, double prop[2] ) {
365 double a[2], b[2];
366 a[0] = 1;
367 a[1] = 0;
368 b[0] = mass2 - sa;
369 b[1] = -mass * width;
370 Com_Divide( a, b, prop );
371}
372void EvtDsToKKpipipi::propagatorRBWl1( double mass2, double mass, double width, double sa,
373 double sb, double sc, double r2, double prop[2] ) {
374 double a[2], b[2];
375 a[0] = 1;
376 a[1] = 0;
377 b[0] = mass2 - sa;
378 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
379 Com_Divide( a, b, prop );
380}
381void EvtDsToKKpipipi::propagatorGS( double mass2, double mass, double width, double sa,
382 double sb, double sc, double r2, double prop[2] ) {
383 double a[2], b[2];
384 double tmp = sb - sc;
385 double tmp1 = sa + tmp;
386 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
387 if ( q2 < 0 ) q2 = 1e-16;
388
389 double tmp2 = mass2 + tmp;
390 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
391 if ( q02 < 0 ) q02 = 1e-16;
392
393 double q = sqrt( q2 );
394 double q0 = sqrt( q02 );
395 double m = sqrt( sa );
396 double q03 = q0 * q02;
397 double tmp3 = log( mass + 2 * q0 ) + 1.2760418309; // log(mass_2Pion) = 1.2760418309;
398
399 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2760418309 );
400 double h0 = GS1 * q0 / mass * tmp3;
401 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
402 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
403 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
404
405 a[0] = 1.0 + d * width / mass;
406 a[1] = 0.0;
407 b[0] = mass2 - sa + width * f;
408 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
409 Com_Divide( a, b, prop );
410}
411void EvtDsToKKpipipi::KPiSLASS( double sa, double sb, double sc, double prop[2] ) {
412 const double m1430 = 1.441;
413 const double sa0 = 2.076481; // m1430*m1430;
414 const double w1430 = 0.193;
415 const double Lass1 = 0.25 / sa0;
416 double tmp = sb - sc;
417 double tmp1 = sa0 + tmp;
418 double q0 = Lass1 * tmp1 * tmp1 - sb;
419 if ( q0 < 0 ) q0 = 1e-16;
420 double tmp2 = sa + tmp;
421 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
422 double q = sqrt( qs );
423 double width = w1430 * q * m1430 / sqrt( sa * q0 );
424 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
425 if ( temp_R < 0 ) temp_R += math_pi;
426 double deltaR = -109.7 * math_pi / 180.0 + temp_R;
427 double temp_F =
428 atan( 0.226 * q / ( 2.0 - 3.8194 * qs ) ); // 2.0*1.07 = 2.14; 1.8*1.07 = 1.926
429 if ( temp_F < 0 ) temp_F += math_pi;
430 double deltaF = 0.1 * math_pi / 180.0 + temp_F;
431 double deltaS = deltaR + 2.0 * deltaF;
432 double t1 = 0.96 * sin( deltaF );
433 double t2 = sin( deltaR );
434 double CF[2], CS[2];
435 CF[0] = cos( deltaF );
436 CF[1] = sin( deltaF );
437 CS[0] = cos( deltaS );
438 CS[1] = sin( deltaS );
439 prop[0] = t1 * CF[0] + t2 * CS[0];
440 prop[1] = t1 * CF[1] + t2 * CS[1];
441}
442
443void EvtDsToKKpipipi::calEvaMy( double* Km, double* Kp, double* Pip1, double* Pip2,
444 double* Pim, double* mass1, double* mass2, double* mass3,
445 double* width1, double* width2, double* width3, double* amp,
446 double* phase, int* modetype, int nstates, double& Result ) {
447 double cof[2], amp_tmp[2], amp_PDF[2], PDF[2];
448 // Km[0] = sqrt(0.243719942 + Km[1]*Km[1] + Km[2]*Km[2] + Km[3]*Km[3]);
449 // Kp[0] = sqrt(0.243719942 + Kp[1]*Kp[1] + Kp[2]*Kp[2] + Kp[3]*Kp[3]);
450 // Pip1[0] = sqrt(0.019479785 + Pip1[1]*Pip1[1] + Pip1[2]*Pip1[2] + Pip1[3]*Pip1[3]);
451 // Pip2[0] = sqrt(0.019479785 + Pip2[1]*Pip2[1] + Pip2[2]*Pip2[2] + Pip2[3]*Pip2[3]);
452 // Pim[0] = sqrt(0.019479785 + Pim[1]*Pim[1] + Pim[2]*Pim[2] + Pim[3]*Pim[3]);
453 // Km[0] = 0.603525; Kp[0] = 0.643143; Pip1[0] = 0.225747; Pip2[0] = 0.332362; Pim[0] =
454 // 0.215164; Km[1] = 0.155547; Kp[1] = 0.0828093; Pip1[1] = 0.0581283; Pip2[1] = -0.25147;
455 // Pim[1] = 0.09105; Km[2] = -0.275397; Kp[2] = 0.261094; Pip1[2] = 0.161565; Pip2[2] =
456 // 0.150701; Pim[2] = -0.0387834; Km[3] = -0.143121; Kp[3] = -0.308036; Pip1[3] = 0.0447196;
457 // Pip2[3] = -0.0709744; Pim[3] = 0.130466; Km[0] = 0.559664; Kp[0] = 0.621872; Pip1[0] =
458 // 0.206021; Pip2[0] = 0.387562; Pim[0] = 0.267452; Km[1] = 0.263265; Kp[1] = -0.20385;
459 // Pip1[1] = -0.0550443; Pip2[1] = 0.182318; Pim[1] = 0.0824861; Km[2] = 0.009757; Kp[2] =
460 // -0.316095; Pip1[2] = -0.061181; Pip2[2] = -0.309785; Pim[2] = 0.211064; Km[3] =
461 // 0.0101476; Kp[3] = 0.0392071; Pip1[3] = -0.127247; Pip2[3] = -0.0389575; Pim[3] =
462 // 0.0264334;
463 //---------------------------------------------------------------
464 double prho1[4], prho2[4], pphi[4], pa11[4], pa12[4], pD1[4], pD2[4], pD3[4], pD4[4],
465 psigma1[4], psigma2[4], pa13[4], pa14[4];
466 for ( int i = 0; i != 4; i++ )
467 {
468 pphi[i] = Km[i] + Kp[i];
469 prho1[i] = Pim[i] + Pip1[i];
470 prho2[i] = Pim[i] + Pip2[i];
471 pa11[i] = prho1[i] + Pip2[i];
472 pa12[i] = prho2[i] + Pip1[i];
473 psigma1[i] = Pim[i] + Pip1[i];
474 psigma2[i] = Pim[i] + Pip2[i];
475 pa13[i] = psigma1[i] + Pip2[i];
476 pa14[i] = psigma2[i] + Pip1[i];
477 pD1[i] = pa11[i] + pphi[i];
478 pD2[i] = pa12[i] + pphi[i];
479 pD3[i] = pa13[i] + pphi[i];
480 pD4[i] = pa14[i] + pphi[i];
481 }
482
483 double skaon1, skaon2, spion1, spion2, spim, sphi, sa11, sa12, srho1, srho2, sD1, sD2, sD3,
484 sD4, ssigma1, ssigma2, sa13, sa14;
485 skaon1 = SCADot( Km, Km );
486 skaon2 = SCADot( Kp, Kp );
487 sphi = SCADot( pphi, pphi );
488
489 spion1 = SCADot( Pip1, Pip1 );
490 spion2 = SCADot( Pip2, Pip2 );
491 spim = SCADot( Pim, Pim );
492
493 srho1 = SCADot( prho1, prho1 );
494 srho2 = SCADot( prho2, prho2 );
495 sa11 = SCADot( pa11, pa11 );
496 sa12 = SCADot( pa12, pa12 );
497
498 sD1 = SCADot( pD1, pD1 );
499 sD2 = SCADot( pD2, pD2 );
500 sD3 = SCADot( pD3, pD3 );
501 sD4 = SCADot( pD4, pD4 );
502
503 ssigma1 = SCADot( psigma1, psigma1 );
504 ssigma2 = SCADot( psigma2, psigma2 );
505 sa13 = SCADot( pa13, pa13 );
506 sa14 = SCADot( pa14, pa14 );
507 //-----------------------------------------------------------------------------------------------
508 double t1phi[4], t1rho1[4], t2a11[4][4], t1sigma1[4], t1a13[4], t1D1[4], t1D3[4], t2D1[4][4],
509 t2D3[4][4], t1rho2[4], t2a12[4][4], t1sigma2[4], t1a14[4], t1D2[4], t1D4[4], t2D2[4][4],
510 t2D4[4][4];
511
512 calt1( Km, Kp, t1phi );
513
514 calt1( Pip1, Pim, t1rho1 );
515 calt1( Pip2, Pim, t1rho2 );
516 calt1( Pip1, Pim, t1sigma1 );
517 calt1( Pip2, Pim, t1sigma2 );
518 calt1( psigma1, Pip2, t1a13 );
519 calt1( psigma2, Pip1, t1a14 );
520
521 calt1( pa11, pphi, t1D1 );
522 calt1( pa12, pphi, t1D2 );
523 calt1( pa13, pphi, t1D3 );
524 calt1( pa14, pphi, t1D4 );
525
526 calt2( prho1, Pip2, t2a11 );
527 calt2( prho2, Pip1, t2a12 );
528
529 calt2( pa11, pphi, t2D1 );
530 calt2( pa12, pphi, t2D2 );
531 calt2( pa13, pphi, t2D3 );
532 calt2( pa14, pphi, t2D4 );
533
534 amp_PDF[0] = 0;
535 amp_PDF[1] = 0;
536 PDF[0] = 0;
537 PDF[1] = 0;
538 double tt1, tt2, tt4;
539 double temp_PDF, tmp1, tmp2, amp_tmp1[2], amp_tmp2[2];
540 double pro[2], pro0[2], pro1[2], pro2[2], pro3[2];
541 double mass1sq, mass2sq, mass3sq;
542 for ( int i = 0; i < nstates; i++ )
543 {
544 tmp1 = 0;
545 tmp2 = 0;
546 temp_PDF = 0;
547 cof[0] = amp[i] * cos( phase[i] );
548 cof[1] = amp[i] * sin( phase[i] );
549 mass1sq = mass1[i] * mass1[i];
550 mass2sq = mass2[i] * mass2[i];
551 mass3sq = mass3[i] * mass3[i];
552
553 amp_tmp[0] = 0;
554 amp_tmp[1] = 0;
555 amp_tmp1[0] = 0;
556 amp_tmp1[1] = 0;
557 amp_tmp2[0] = 0;
558 amp_tmp2[1] = 0;
559 pro[0] = 0;
560 pro[1] = 0;
561 pro1[0] = 0;
562 pro1[1] = 0;
563 pro2[0] = 0;
564 pro2[1] = 0;
565 pro3[0] = 0;
566 pro3[1] = 0;
567 if ( modetype[i] == 1 )
568 {
569 double B_phi = -1.0, B_rho1 = -1.0;
570 propagatorRBWl1( mass1sq, mass1[i], width1[i], sphi, skaon1, skaon2, rRes2, pro0 );
571 propagatorGS( mass2sq, mass2[i], width2[i], srho1, spion1, spim, rRes2, pro1 );
572 propagatorRBW( mass3sq, mass3[i], width3[i], sa11, srho1, spion2, rRes2, 0, pro2 );
573 Com_Multi( pro0, pro1, pro3 );
574 Com_Multi( pro2, pro3, pro );
575 if ( B_phi < 0.0 ) B_phi = barrier( 1, sphi, skaon1, skaon2, rRes2, 1.019461 );
576 if ( B_rho1 < 0.0 ) B_rho1 = barrier( 1, srho1, spion1, spim, rRes2, 0.77526 );
577 for ( int a = 0; a < 4; a++ )
578 {
579 for ( int j = 0; j < 4; j++ )
580 {
581 temp_PDF += ( G[a][j] - pa11[a] * pa11[j] / sa11 ) * t1phi[a] * t1rho1[j] * G[a][a] *
582 G[j][j];
583 }
584 }
585 tmp1 = B_phi * B_rho1 * temp_PDF;
586 amp_tmp1[0] = tmp1 * pro[0];
587 amp_tmp1[1] = tmp1 * pro[1];
588
589 temp_PDF = 0;
590 double B_rho2 = -1.0;
591 propagatorRBWl1( mass1sq, mass1[i], width1[i], sphi, skaon1, skaon2, rRes2, pro0 );
592 propagatorGS( mass2sq, mass2[i], width2[i], srho2, spion2, spim, rRes2, pro1 );
593 propagatorRBW( mass3sq, mass3[i], width3[i], sa12, srho2, spion1, rRes2, 0, pro2 );
594 Com_Multi( pro0, pro1, pro3 );
595 Com_Multi( pro2, pro3, pro );
596 if ( B_phi < 0.0 ) B_phi = barrier( 1, sphi, skaon1, skaon2, rRes2, 1.019461 );
597 if ( B_rho2 < 0.0 ) B_rho2 = barrier( 1, srho2, spion2, spim, rRes2, 0.77526 );
598 for ( int a = 0; a < 4; a++ )
599 {
600 for ( int j = 0; j < 4; j++ )
601 {
602 temp_PDF += ( G[a][j] - pa12[a] * pa12[j] / sa12 ) * t1phi[a] * t1rho2[j] * G[a][a] *
603 G[j][j];
604 }
605 }
606 tmp2 = B_phi * B_rho2 * temp_PDF;
607 amp_tmp2[0] = tmp2 * pro[0];
608 amp_tmp2[1] = tmp2 * pro[1];
609 }
610 if ( modetype[i] == 3 )
611 {
612 double B_phi = -1.0, B_rho1 = -1.0, B_phia11 = -1.0;
613 propagatorRBWl1( mass1sq, mass1[i], width1[i], sphi, skaon1, skaon2, rRes2, pro0 );
614 propagatorGS( mass2sq, mass2[i], width2[i], srho1, spion1, spim, rRes2, pro1 );
615 propagatorRBW( mass3sq, mass3[i], width3[i], sa11, srho1, spion2, rRes2, 0, pro2 );
616 // printf("mass1[0]:%.6f\n",mass1[0]);
617 // printf("mass1[1]:%.6f\n",mass1[1]);
618 // printf("width1[0]:%.6f\n",width1[0]);
619 // printf("width1[1]:%.6f\n",width1[1]);
620 // printf("pro0[0]:%.6f\n",pro0[0]);
621 // printf("pro0[1]:%.6f\n",pro0[1]);
622 // printf("pro1[0]:%.6f\n",pro1[0]);
623 // printf("pro1[1]:%.6f\n",pro1[1]);
624 // printf("pro2[0]:%.6f\n",pro2[0]);
625 // printf("pro2[1]:%.6f\n",pro2[1]);
626 Com_Multi( pro0, pro1, pro3 );
627 Com_Multi( pro2, pro3, pro );
628 if ( B_phi < 0.0 ) B_phi = barrier( 1, sphi, skaon1, skaon2, rRes2, 1.019461 );
629 if ( B_rho1 < 0.0 ) B_rho1 = barrier( 1, srho1, spion1, spim, rRes2, 0.77526 );
630 if ( B_phia11 < 0.0 ) B_phia11 = barrier( 1, sD1, sphi, sa11, rD2, 1.9683 );
631 for ( int w = 0; w < 4; w++ )
632 {
633 tt1 = pD1[w] * G[w][w];
634 for ( int j = 0; j < 4; j++ )
635 {
636 tt2 = t1D1[j] * G[j][j];
637 for ( int m = 0; m < 4; m++ )
638 {
639 for ( int k = 0; k < 4; k++ )
640 {
641 tt4 = t1phi[k] * G[k][k];
642 for ( int l = 0; l < 4; l++ )
643 {
644 temp_PDF += E[w][m][k][l] * G[m][m] * tt1 * tt2 *
645 ( G[m][j] - pa11[m] * pa11[j] / sa11 ) * tt4 * t1rho1[l] * G[l][l];
646 }
647 }
648 }
649 }
650 }
651 tmp1 = B_phi * B_rho1 * B_phia11 * temp_PDF; // a1(l=0),B_a11=1
652 amp_tmp1[0] = tmp1 * pro[0];
653 amp_tmp1[1] = tmp1 * pro[1];
654
655 temp_PDF = 0;
656 double B_rho2 = -1.0, B_phia12 = -1.0;
657 propagatorRBWl1( mass1sq, mass1[i], width1[i], sphi, skaon1, skaon2, rRes2, pro0 );
658 propagatorGS( mass2sq, mass2[i], width2[i], srho2, spion2, spim, rRes2, pro1 );
659 propagatorRBW( mass3sq, mass3[i], width3[i], sa12, srho2, spion1, rRes2, 0, pro2 );
660 Com_Multi( pro0, pro1, pro3 );
661 Com_Multi( pro2, pro3, pro );
662 if ( B_phi < 0.0 ) B_phi = barrier( 1, sphi, skaon1, skaon2, rRes2, 1.019461 );
663 if ( B_rho2 < 0.0 ) B_rho2 = barrier( 1, srho2, spion2, spim, rRes2, 0.77526 );
664 if ( B_phia12 < 0.0 ) B_phia12 = barrier( 1, sD2, sphi, sa12, rD2, 1.9683 );
665 for ( int w = 0; w < 4; w++ )
666 {
667 tt1 = pD2[w] * G[w][w];
668 for ( int j = 0; j < 4; j++ )
669 {
670 tt2 = t1D2[j] * G[j][j];
671 for ( int m = 0; m < 4; m++ )
672 {
673 for ( int k = 0; k < 4; k++ )
674 {
675 tt4 = t1phi[k] * G[k][k];
676 for ( int l = 0; l < 4; l++ )
677 {
678 temp_PDF += E[w][m][k][l] * G[m][m] * tt1 * tt2 *
679 ( G[m][j] - pa12[m] * pa12[j] / sa12 ) * tt4 * t1rho2[l] * G[l][l];
680 }
681 }
682 }
683 }
684 }
685 tmp2 = B_phi * B_rho2 * B_phia12 * temp_PDF; // a1(l=0),B_a12=1
686 amp_tmp2[0] = tmp2 * pro[0];
687 amp_tmp2[1] = tmp2 * pro[1];
688 }
689 if ( modetype[i] == 13 )
690 {
691 amp_tmp1[0] = 1.;
692 amp_tmp1[1] = 0.;
693 amp_tmp2[0] = 1.;
694 amp_tmp2[1] = 0.;
695 }
696 amp_tmp[0] = amp_tmp1[0] + amp_tmp2[0];
697 amp_tmp[1] = amp_tmp1[1] + amp_tmp2[1];
698 Com_Multi( amp_tmp, cof, amp_PDF );
699 PDF[0] += amp_PDF[0];
700 PDF[1] += amp_PDF[1];
701 }
702 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
703 // printf("pdf:%.6f\n",value);
704 if ( value <= 0 ) value = 1e-20;
705 Result = value;
706 // cout<<"!!!!!!!!!!!!"<<setprecision(20)<<value<<endl;
707}
mass
double mass
Definition CosmicGenerator.cxx:128
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
w
double w
Definition EvtPhokharaDef.hh:28
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
EvtDsToKKpipipi::~EvtDsToKKpipipi
virtual ~EvtDsToKKpipipi()
Definition EvtDsToKKpipipi.cc:34
EvtDsToKKpipipi::getName
void getName(std::string &name)
Definition EvtDsToKKpipipi.cc:36
EvtDsToKKpipipi::decay
void decay(EvtParticle *p)
Definition EvtDsToKKpipipi.cc:135
EvtDsToKKpipipi::clone
EvtDecayBase * clone()
Definition EvtDsToKKpipipi.cc:38
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