BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToKSKmPipPip.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: EvtDsToKSKmPipPip.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 "EvtDsToKSKmPipPip.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
33EvtDsToKSKmPipPip::~EvtDsToKSKmPipPip() {}
34
35void EvtDsToKSKmPipPip::getName( std::string& model_name ) { model_name = "DsToKSKmPipPip"; }
36
37EvtDecayBase* EvtDsToKSKmPipPip::clone() { return new EvtDsToKSKmPipPip; }
38
39void EvtDsToKSKmPipPip::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] = -1.5871e+00;
52 phi[2] = -6.4308e+00;
53 phi[3] = -4.4006e+00;
54 phi[4] = 4.7540e+00;
55 phi[5] = 4.2785e+00;
56 phi[6] = 6.3523e+00;
57 phi[7] = 6.3523e+00;
58 phi[8] = 5.4902e+00;
59 phi[9] = 5.2024e+00;
60 phi[10] = -3.9776e+00;
61
62 rho[1] = 3.6108e-01;
63 rho[2] = 4.5990e-01;
64 rho[3] = 1.7722e+00;
65 rho[4] = 2.1434e+00;
66 rho[5] = 1.2504e+00;
67 rho[6] = 1.0069e+00;
68 rho[7] = 1.0069e+00;
69 rho[8] = 1.7351e+00;
70 rho[9] = 4.1005e+00;
71 rho[10] = 3.7821e+00;
72
73 modetype[0] = 1;
74 modetype[1] = 2;
75 modetype[2] = 3;
76 modetype[3] = 21;
77 modetype[4] = 22;
78 modetype[5] = 8;
79 modetype[6] = 4;
80 modetype[7] = 14;
81 modetype[8] = 17;
82 modetype[9] = 5;
83 modetype[10] = 20;
84
85 phi[0] = 0;
86 rho[0] = 1;
87
88 // cout << "DsToKSKmPipPip :" << endl;
89 // for (int i=0; i<11; i++) {
90 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
91 // }
92
93 G1510 = 1.0000e+00;
94 GKst0 = 4.7300e-02;
95 GKstp = 5.0300e-02;
96 Ga0_980 = 5.0000e-02;
97 Geta1475 = 9.0000e-02;
98 Gf1285 = 2.2700e-02;
99 m1510 = 1.6200e+00;
100 mKst0 = 8.9555e-01;
101 mKstp = 8.9176e-01;
102 ma0_980 = 9.9000e-01;
103 meta1475 = 1.4750e+00;
104 mf1285 = 1.2819e+00;
105
106 mD = 1.86486;
107 metap = 0.95778;
108 mk0 = 0.497614;
109 mass_Kaon = 0.49368;
110 mass_Eta = 0.547862;
111 mass_Pion = 0.13957;
112 math_pi = 3.1415926;
113 mass_Pion2 = 0.0194797849;
114 mass_2Pion = 0.27914;
115 math_2pi = 6.2831852;
116 rD2 = 25.0; // 5*5
117 rRes2 = 9.0; // 3*3
118 g1 = 0.5468;
119 g2 = 0.23;
120 GS1 = 0.636619783;
121 GS2 = 0.01860182466;
122 GS3 = 0.1591549458; // 1/(2*math_2pi)
123 GS4 = 0.00620060822; // mass_Pion2/math_pi
124
125 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
126 int EE[4][4][4][4] = {
127 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
128 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, -1, 0 } },
129 { { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 1, 0, 0 } },
130 { { 0, 0, 0, 0 }, { 0, 0, 1, 0 }, { 0, -1, 0, 0 }, { 0, 0, 0, 0 } } },
131 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 1, 0 } },
132 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
133 { { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 } },
134 { { 0, 0, -1, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 } } },
135 { { { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, -1, 0, 0 } },
136 { { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 } },
137 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
138 { { 0, 1, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } },
139 { { { 0, 0, 0, 0 }, { 0, 0, -1, 0 }, { 0, 1, 0, 0 }, { 0, 0, 0, 0 } },
140 { { 0, 0, 1, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 } },
141 { { 0, -1, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
142 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } } };
143 for ( int i = 0; i < 4; i++ )
144 {
145 for ( int j = 0; j < 4; j++ )
146 {
147 G[i][j] = GG[i][j];
148 for ( int k = 0; k < 4; k++ )
149 {
150 for ( int l = 0; l < 4; l++ ) { E[i][j][k][l] = EE[i][j][k][l]; }
151 }
152 }
153 }
154}
155
156void EvtDsToKSKmPipPip::initProbMax() {
157 // setProbMax(75000);
158 setProbMax( 95050 );
159 // setProbMax(1.0);
160}
161
162void EvtDsToKSKmPipPip::decay( EvtParticle* p ) {
163 /*
164 double maxprob=0;
165 for(int ir=0;ir<=60000000;ir++){
166 p->initializePhaseSpace(getNDaug(),getDaugs());
167 EvtVector4R _ks = p->getDaug(0)->getP4();
168 EvtVector4R _km = p->getDaug(1)->getP4();
169 EvtVector4R _pip = p->getDaug(2)->getP4();
170 EvtVector4R _pip2 = p->getDaug(3)->getP4();
171
172 double _Pip[4],_Pip2[4],_Ks[4],_Km[4];
173 _Ks[0] = _ks.get(0); _Km[0] = _km.get(0); _Pip[0] = _pip.get(0); _Pip2[0] =
174 _pip2.get(0); _Ks[1] = _ks.get(1); _Km[1] = _km.get(1); _Pip[1] = _pip.get(1); _Pip2[1] =
175 _pip2.get(1); _Ks[2] = _ks.get(2); _Km[2] = _km.get(2); _Pip[2] = _pip.get(2); _Pip2[2] =
176 _pip2.get(2); _Ks[3] = _ks.get(3); _Km[3] = _km.get(3); _Pip[3] = _pip.get(3); _Pip2[3] =
177 _pip2.get(3);
178
179 double _prob;
180 int g0[11]={1,1,1,1,1,1,1,1,1,1,1};
181 int nstates=11;
182
183 calEvaMy(_Ks, _Km, _Pip, _Pip2, mass, width, rho, phi, g0, modetype, nstates, _prob);
184 if(_prob>maxprob) {
185 maxprob=_prob;
186 printf("ir = %d,maxprob = %.10f,prob = %.10f\n\n",ir,maxprob,_prob);
187 }
188 }
189 printf("maxprob = %.10f\n", maxprob);*/
190
191 p->initializePhaseSpace( getNDaug(), getDaugs() );
192 EvtVector4R ks = p->getDaug( 0 )->getP4();
193 EvtVector4R km = p->getDaug( 1 )->getP4();
194 EvtVector4R pip = p->getDaug( 2 )->getP4();
195 EvtVector4R pip2 = p->getDaug( 3 )->getP4();
196
197 double Pip[4], Pip2[4], Km[4], Ks[4];
198 Ks[0] = ks.get( 0 );
199 Km[0] = km.get( 0 );
200 Pip[0] = pip.get( 0 );
201 Pip2[0] = pip2.get( 0 );
202 Ks[1] = ks.get( 1 );
203 Km[1] = km.get( 1 );
204 Pip[1] = pip.get( 1 );
205 Pip2[1] = pip2.get( 1 );
206 Ks[2] = ks.get( 2 );
207 Km[2] = km.get( 2 );
208 Pip[2] = pip.get( 2 );
209 Pip2[2] = pip2.get( 2 );
210 Ks[3] = ks.get( 3 );
211 Km[3] = km.get( 3 );
212 Pip[3] = pip.get( 3 );
213 Pip2[3] = pip2.get( 3 );
214
215 // Ks[0] = 0.6096125554; Ks[1] = 0.2036145752; Ks[2] = -0.2256821642; Ks[3] =
216 // 0.1884652565 ; Km[0] = 0.7106732604; Km[1] = -0.2599714599; Km[2] =
217 // 0.3643877399; Km[3] = 0.2469330230 ; Pip[0] = 0.4177584141; Pip[1] =
218 // -0.0562004660; Pip[2] = 0.2821155804; Pip[3] =-0.0020089709 ; Pip2[0] =
219 // 0.3945120988; Pip2[1] = -0.1225632042; Pip2[2] = -0.3645928242; Pip2[3] = 0.0521817527
220 // ;
221
222 // Ks[0] = 0.6439438863; Ks[1] = 0.3368405824; Ks[2] = 0.2308441601; Ks[3] = 0.0171298406;
223 // Km[0] = 0.7146391007; Km[1] = 0.4061832461; Km[2] = -0.1841669581; Km[3] = 0.2609401579;
224 // Pip[0] = 0.2817222681; Pip[1] = -0.1005161765;Pip[2] = -0.1172908044;Pip[3] =
225 // -0.1898077097; Pip2[0] = 0.3638598699;Pip2[1] = -0.3106229764;Pip2[2] =
226 // 0.0905020525;Pip2[3] = 0.0907574503;
227
228 double value;
229 int g0[11] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
230 int nstates = 11;
231
232 calEvaMy( Ks, Km, Pip, Pip2, mass, width, rho, phi, g0, modetype, nstates, value );
233 setProb( value );
234
235 return;
236}
237
238void EvtDsToKSKmPipPip::Com_Multi( double a1[2], double a2[2], double res[2] ) {
239 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
240 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
241}
242void EvtDsToKSKmPipPip::Com_Divide( double a1[2], double a2[2], double res[2] ) {
243 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / ( a2[0] * a2[0] + a2[1] * a2[1] );
244 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / ( a2[0] * a2[0] + a2[1] * a2[1] );
245}
246double EvtDsToKSKmPipPip::SCADot( double a1[4], double a2[4] ) {
247 double _cal = 0.;
248 _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
249 return _cal;
250}
251double EvtDsToKSKmPipPip::barrier( int l, double sa, double sb, double sc, double r2 ) {
252 double F;
253 double tmp = sa + sb - sc;
254 double q = 0.25 * tmp * tmp / sa - sb;
255 if ( q < 0 ) q = 1e-16;
256 double z = q * r2;
257 if ( l == 1 ) { F = sqrt( 2.0 * z / ( 1.0 + z ) ); }
258 else if ( l == 2 )
259 {
260 double z2 = z * z;
261 F = sqrt( 13.0 * z2 / ( 9.0 + 3.0 * z + z2 ) );
262 }
263 else { F = 1.0; }
264 return F;
265}
266void EvtDsToKSKmPipPip::calt1( double daug1[4], double daug0[4], double t1[4] ) {
267 double p, pq;
268 double pa[4], qa[4];
269 for ( int i = 0; i < 4; i++ )
270 {
271 pa[i] = daug1[i] + daug0[i];
272 qa[i] = daug1[i] - daug0[i];
273 }
274 p = SCADot( pa, pa );
275 pq = SCADot( pa, qa );
276 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - pq / p * pa[i]; }
277}
278void EvtDsToKSKmPipPip::calt2( double daug1[4], double daug0[4], double t2[4][4] ) {
279 double p, r;
280 double pa[4], t1[4];
281 calt1( daug1, daug0, t1 );
282 r = SCADot( t1, t1 ) / 3.0;
283 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug0[i]; }
284 p = SCADot( pa, pa );
285 for ( int i = 0; i < 4; i++ )
286 {
287 for ( int j = 0; j < 4; j++ )
288 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
289 }
290}
291
292double EvtDsToKSKmPipPip::wid( double mass, double sa, double sb, double sc, double r2,
293 int l ) {
294 double widm = 0.;
295 double sa0 = mass * mass;
296 double m = sqrt( sa );
297 double tmp = sb - sc;
298 double tmp1 = sa + tmp;
299 double q = 0.25 * tmp1 * tmp1 / sa - sb;
300 if ( q < 0 ) q = 1e-16;
301 double tmp2 = sa0 + tmp;
302 double q0 = 0.25 * tmp2 * tmp2 / sa0 - sb;
303 if ( q0 < 0 ) q0 = 1e-16;
304 double z = q * r2;
305 double z0 = q0 * r2;
306 double t = q / q0;
307 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
308 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
309 else if ( l == 2 )
310 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
311 return widm;
312}
313void EvtDsToKSKmPipPip::Flatte_rhoab( double sa, double sb, double sc, double rho[2] ) {
314 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
315 if ( q > 0 )
316 {
317 rho[0] = 2 * sqrt( q / sa );
318 rho[1] = 0;
319 }
320 else if ( q < 0 )
321 {
322 rho[0] = 0;
323 rho[1] = 2 * sqrt( -q / sa );
324 }
325}
326void EvtDsToKSKmPipPip::propagatora0980( double mass, double sx, double* sb, double* sc,
327 double prop[2] ) {
328 double unit[2] = { 1.0 };
329 double ci[2] = { 0, 1 };
330 double rho1[2];
331 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
332 double rho2[2];
333 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
334
335 double gKK_a980 = 0.892 * 0.341, gPiEta_a980 = 0.341;
336 double tmp1[2] = { gKK_a980, 0 };
337 double tmp11[2];
338 double tmp2[2] = { gPiEta_a980, 0 };
339 double tmp22[2];
340 Com_Multi( tmp1, rho1, tmp11 );
341 Com_Multi( tmp2, rho2, tmp22 );
342 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
343 double tmp31[2];
344 Com_Multi( tmp3, ci, tmp31 );
345 double tmp4[2] = { mass * mass - sx - tmp31[0], -1.0 * tmp31[1] };
346 Com_Divide( unit, tmp4, prop );
347}
348void EvtDsToKSKmPipPip::propagatorRBW( double mass, double width, double sa, double sb,
349 double sc, double r2, int l, double prop[2] ) {
350 double a[2], b[2];
351 a[0] = 1;
352 a[1] = 0;
353 b[0] = mass * mass - sa;
354 b[1] = -mass * width * wid( mass, sa, sb, sc, r2, l );
355 Com_Divide( a, b, prop );
356}
357void EvtDsToKSKmPipPip::KPiSLASS( double sa, double sb, double sc, double prop[2] ) {
358 const double m1430 = 1.441;
359 const double sa0 = 2.076481; // m1410*m1410;
360 const double w1430 = 0.193;
361 const double Lass1 = 0.25 / sa0;
362 double tmp = sb - sc;
363 double tmp1 = sa0 + tmp;
364 double q0 = Lass1 * tmp1 * tmp1 - sb;
365 if ( q0 < 0 ) q0 = 1e-16;
366 double tmp2 = sa + tmp;
367 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
368 double q = sqrt( qs );
369 double width = w1430 * q * m1430 / sqrt( sa * q0 );
370 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
371 if ( temp_R < 0 ) temp_R += math_pi;
372 double deltaR = -1.915 + temp_R; // fiR=-109.7
373 double temp_F =
374 atan( 0.226 * q / ( 2.0 - 3.8194 * qs ) ); // 2.0*0.113=0.226; 0.113*33.8=1.926; a=0.113;
375 // r=-33.8
376 if ( temp_F < 0 ) temp_F += math_pi;
377 double deltaF = 0.002 + temp_F; // fiF=0.1
378 double deltaS = deltaR + 2.0 * deltaF; ///
379 double t1 = 0.96 * sin( deltaF ); // F= 0.96
380 double t2 = sin( deltaR );
381 double CF[2], CS[2];
382 CF[0] = cos( deltaF );
383 CF[1] = sin( deltaF );
384 CS[0] = cos( deltaS );
385 CS[1] = sin( deltaS );
386 prop[0] = t1 * CF[0] + t2 * CS[0];
387 prop[1] = t1 * CF[1] + t2 * CS[1];
388}
389
390void EvtDsToKSKmPipPip::calEvaMy( double* Ks, double* Km, double* Pip, double* Pip2,
391 double* mass1, double* width1, double* amp, double* phase,
392 int* g0, int* modetype, int nstates, double& Result ) {
393 double pD[4], pKsPip1[4], pKmPip2[4], pKsPip2[4], pKmPip1[4], pKsKmPip1[4], pKsKmPip2[4],
394 pKsKm[4], pKmPipPip[4];
395 for ( int i = 0; i != 4; i++ )
396 {
397 pD[i] = Ks[i] + Pip[i] + Km[i] + Pip2[i];
398 pKsPip1[i] = Ks[i] + Pip[i];
399 pKsPip2[i] = Ks[i] + Pip2[i];
400 pKmPip1[i] = Km[i] + Pip[i];
401 pKmPip2[i] = Km[i] + Pip2[i];
402 pKsKm[i] = Km[i] + Ks[i];
403 pKsKmPip1[i] = pKsPip1[i] + Km[i];
404 pKsKmPip2[i] = pKsPip2[i] + Km[i];
405 pKmPipPip[i] = Km[i] + Pip[i] + Pip2[i];
406 }
407 double sPip2, sPip1, sKs, sKm, sKsPip1, sKmPip2, sD, sKsPip2, sKmPip1, sKsKmPip1, sKsKmPip2,
408 sKsKm, sKmPipPip;
409 sD = SCADot( pD, pD );
410 sKs = SCADot( Ks, Ks );
411 sKm = SCADot( Km, Km );
412 sPip1 = SCADot( Pip, Pip );
413 sPip2 = SCADot( Pip2, Pip2 );
414 sKsPip1 = SCADot( pKsPip1, pKsPip1 );
415 sKsPip2 = SCADot( pKsPip2, pKsPip2 );
416 sKmPip1 = SCADot( pKmPip1, pKmPip1 );
417 sKmPip2 = SCADot( pKmPip2, pKmPip2 );
418 sKsKmPip1 = SCADot( pKsKmPip1, pKsKmPip1 );
419 sKsKmPip2 = SCADot( pKsKmPip2, pKsKmPip2 );
420 sKmPipPip = SCADot( pKmPipPip, pKmPipPip );
421 sKsKm = SCADot( pKsKm, pKsKm );
422
423 double t1KsPip1[4], t1KmPip2[4], t1KsPip2[4], t1KmPip1[4], t2KsKmPip1_kspip[4][4],
424 t2KsKmPip2_kspip[4][4], t2KsKmPip1_kmpip[4][4], t2KsKmPip2_kmpip[4][4],
425 t1KsKmPip1_kspip[4], t1KsKmPip2_kspip[4], t1KsKmPip1_kmpip[4], t1KsKmPip2_kmpip[4],
426 t1KsKmPip1_ksk[4], t1KsKmPip2_ksk[4];
427
428 calt1( Ks, Pip, t1KsPip1 );
429 calt1( Ks, Pip2, t1KsPip2 );
430 calt1( Km, Pip, t1KmPip1 );
431 calt1( Km, Pip2, t1KmPip2 );
432 calt1( pKsPip1, Km, t1KsKmPip1_kspip );
433 calt1( pKsPip2, Km, t1KsKmPip2_kspip );
434 calt1( pKmPip1, Ks, t1KsKmPip1_kmpip );
435 calt1( pKmPip2, Ks, t1KsKmPip2_kmpip );
436 calt1( pKsKm, Pip, t1KsKmPip1_ksk );
437 calt1( pKsKm, Pip2, t1KsKmPip2_ksk );
438
439 calt2( pKsPip1, Km, t2KsKmPip1_kspip );
440 calt2( pKsPip2, Km, t2KsKmPip2_kspip );
441 calt2( pKmPip1, Ks, t2KsKmPip1_kmpip );
442 calt2( pKmPip2, Ks, t2KsKmPip2_kmpip );
443
444 double cof[2], amp_tmp[2], amp_PDF[2], PDF[2];
445 amp_PDF[0] = 0;
446 amp_PDF[1] = 0;
447 PDF[0] = 0;
448 PDF[1] = 0;
449 double temp_PDF, tmp1, tmp2, tmp1_a, tmp2_a, temp_PDF_a, tt1, tt2, tt3, tt4;
450 double pro1b[2], pro2b[2], pro1b_a[2], pro2b_a[2], pro_ksk[2], pro3[2], pro4[2], pro1[2],
451 pro2[2];
452 double pro_kspi1[2], pro_kspi2[2], pro_kpi1[2], pro_kpi2[2], pro_kpi2sw[2], pro_kpi1sw[2],
453 pro_kspi2sw[2], pro_kspi1sw[2];
454 double ispro_kspi1 = 0, ispro_kspi2 = 0, ispro_kpi1 = 0, ispro_kpi2 = 0, ispro_kpi2sw = 0,
455 ispro_kpi1sw = 0, ispro_kspi2sw = 0, ispro_kspi1sw = 0, ispro_ksk = 0;
456 double barr_kskpi1_ksk = -1.0, barr_kskpi2_ksk = -1.0, barr_kskpi1_kspi = -1.0,
457 barr_kskpi2_kspi = -1.0, barr_kskpi1_kpi = -1.0, barr_kskpi2_kpi = -1.0,
458 barr_kspi1 = -1.0, barr_kspi2 = -1.0, barr_kpi1 = -1.0, barr_kpi2 = -1.0,
459 barr_kspi1_kpi2 = -1.0, barr_kspi2_kpi1 = -1.0, barr_ds_kskpi1 = -1.0,
460 barr_ds_kskpi2 = -1.0;
461
462 double t1D[4], t2D[4][4], B[3];
463 double sKs2[2] = { sKs, mass_Pion * mass_Pion };
464 double sKm2[2] = { sKm, mass_Eta * mass_Eta };
465
466 for ( int i = 0; i < nstates; i++ )
467 {
468 amp_tmp[0] = 0;
469 amp_tmp[1] = 0;
470 tmp1 = 0;
471 tmp2 = 0;
472 temp_PDF = 0;
473 tmp1_a = 0;
474 tmp2_a = 0;
475 temp_PDF_a = 0;
476
477 cof[0] = amp[i] * cos( phase[i] );
478 cof[1] = amp[i] * sin( phase[i] );
479
480 // printf("phase= %i, %.5f, %.5f\n", i, amp[i],phase[i]);
481
482 // Kst+ Kst0
483 if ( modetype[i] == 1 )
484 {
485 if ( ispro_kspi1 == 0 )
486 {
487 propagatorRBW( mKstp, GKstp, sKsPip1, sKs, sPip1, rRes2, 1, pro_kspi1 );
488 ispro_kspi1 = 1;
489 }
490 if ( ispro_kpi2 == 0 )
491 {
492 propagatorRBW( mKst0, GKst0, sKmPip2, sKm, sPip2, rRes2, 1, pro_kpi2 );
493 ispro_kpi2 = 1;
494 }
495 Com_Multi( pro_kspi1, pro_kpi2, pro1b );
496
497 if ( barr_kspi1 < 0.0 ) barr_kspi1 = barrier( 1, sKsPip1, sKs, sPip1, rRes2 );
498 if ( barr_kpi2 < 0.0 ) barr_kpi2 = barrier( 1, sKmPip2, sKm, sPip2, rRes2 );
499 for ( int i = 0; i < 4; i++ ) { temp_PDF += G[i][i] * t1KsPip1[i] * t1KmPip2[i]; }
500 tmp1 = barr_kspi1 * barr_kpi2 * temp_PDF;
501 // printf("%.15f, %.15f , %.15f, %.15f \n",
502 // pro_kspi1[0],pro_kpi2[0],pro_kspi1[1],pro_kpi2[1]);
503
504 temp_PDF = 0;
505 if ( ispro_kspi2 == 0 )
506 {
507 propagatorRBW( mKstp, GKstp, sKsPip2, sKs, sPip2, rRes2, 1, pro_kspi2 );
508 ispro_kspi2 = 1;
509 }
510 if ( ispro_kpi1 == 0 )
511 {
512 propagatorRBW( mKst0, GKst0, sKmPip1, sKm, sPip1, rRes2, 1, pro_kpi1 );
513 ispro_kpi1 = 1;
514 }
515 Com_Multi( pro_kspi2, pro_kpi1, pro2b );
516
517 if ( barr_kspi2 < 0.0 ) barr_kspi2 = barrier( 1, sKsPip2, sKs, sPip2, rRes2 );
518 if ( barr_kpi1 < 0.0 ) barr_kpi1 = barrier( 1, sKmPip1, sKm, sPip1, rRes2 );
519 for ( int i = 0; i < 4; i++ ) { temp_PDF += G[i][i] * t1KsPip2[i] * t1KmPip1[i]; }
520 tmp2 = barr_kspi2 * barr_kpi1 * temp_PDF;
521 // printf("%.15f, %.15f , %.15f, %.15f \n",
522 // pro_kspi2[0],pro_kpi1[0],pro_kspi2[1],pro_kpi1[1]);
523 }
524
525 // kstp_kst0_p
526 if ( modetype[i] == 2 )
527 {
528 if ( ispro_kspi1 == 0 )
529 {
530 propagatorRBW( mKstp, GKstp, sKsPip1, sKs, sPip1, rRes2, 1, pro_kspi1 );
531 ispro_kspi1 = 1;
532 }
533 if ( ispro_kpi2 == 0 )
534 {
535 propagatorRBW( mKst0, GKst0, sKmPip2, sKm, sPip2, rRes2, 1, pro_kpi2 );
536 ispro_kpi2 = 1;
537 }
538 Com_Multi( pro_kspi1, pro_kpi2, pro1b );
539
540 if ( barr_kspi1 < 0.0 ) barr_kspi1 = barrier( 1, sKsPip1, sKs, sPip1, rRes2 );
541 if ( barr_kpi2 < 0.0 ) barr_kpi2 = barrier( 1, sKmPip2, sKm, sPip2, rRes2 );
542 if ( barr_kspi1_kpi2 < 0.0 ) barr_kspi1_kpi2 = barrier( 1, sD, sKsPip1, sKmPip2, rD2 );
543 calt1( pKsPip1, pKmPip2, t1D );
544 for ( int i = 0; i < 4; i++ )
545 {
546 for ( int j = 0; j < 4; j++ )
547 {
548 tt1 = G[i][i] * G[j][j] * pD[i] * t1D[j];
549 for ( int k = 0; k < 4; k++ )
550 {
551 tt2 = t1KsPip1[k] * G[k][k];
552 for ( int l = 0; l < 4; l++ )
553 { temp_PDF += E[i][j][k][l] * t1KmPip2[l] * G[l][l] * tt1 * tt2; }
554 }
555 }
556 }
557 tmp1 = barr_kspi1 * barr_kpi2 * barr_kspi1_kpi2 * temp_PDF;
558
559 temp_PDF = 0;
560 if ( ispro_kspi2 == 0 )
561 {
562 propagatorRBW( mKstp, GKstp, sKsPip2, sKs, sPip2, rRes2, 1, pro_kspi2 );
563 ispro_kspi2 = 1;
564 }
565 if ( ispro_kpi1 == 0 )
566 {
567 propagatorRBW( mKst0, GKst0, sKmPip1, sKm, sPip1, rRes2, 1, pro_kpi1 );
568 ispro_kpi1 = 1;
569 }
570 Com_Multi( pro_kspi2, pro_kpi1, pro2b );
571
572 if ( barr_kspi2 < 0.0 ) barr_kspi2 = barrier( 1, sKsPip2, sKs, sPip2, rRes2 );
573 if ( barr_kpi1 < 0.0 ) barr_kpi1 = barrier( 1, sKmPip1, sKm, sPip1, rRes2 );
574 if ( barr_kspi2_kpi1 < 0.0 ) barr_kspi2_kpi1 = barrier( 1, sD, sKsPip2, sKmPip1, rD2 );
575 calt1( pKsPip2, pKmPip1, t1D );
576 for ( int i = 0; i < 4; i++ )
577 {
578 for ( int j = 0; j < 4; j++ )
579 {
580 tt1 = G[i][i] * G[j][j] * pD[i] * t1D[j];
581 for ( int k = 0; k < 4; k++ )
582 {
583 tt2 = t1KsPip2[k] * G[k][k];
584 for ( int l = 0; l < 4; l++ )
585 { temp_PDF += E[i][j][k][l] * t1KmPip1[l] * G[l][l] * tt1 * tt2; }
586 }
587 }
588 }
589 tmp2 = barr_kspi2 * barr_kpi1 * barr_kspi2_kpi1 * temp_PDF;
590 }
591
592 // kstp_kst0_d
593 if ( modetype[i] == 3 )
594 {
595 if ( ispro_kspi1 == 0 )
596 {
597 propagatorRBW( mKstp, GKstp, sKsPip1, sKs, sPip1, rRes2, 1, pro_kspi1 );
598 ispro_kspi1 = 1;
599 }
600 if ( ispro_kpi2 == 0 )
601 {
602 propagatorRBW( mKst0, GKst0, sKmPip2, sKm, sPip2, rRes2, 1, pro_kpi2 );
603 ispro_kpi2 = 1;
604 }
605 Com_Multi( pro_kspi1, pro_kpi2, pro1b );
606
607 if ( barr_kspi1 < 0.0 ) barr_kspi1 = barrier( 1, sKsPip1, sKs, sPip1, rRes2 );
608 if ( barr_kpi2 < 0.0 ) barr_kpi2 = barrier( 1, sKmPip2, sKm, sPip2, rRes2 );
609 calt2( pKsPip1, pKmPip2, t2D );
610 for ( int i = 0; i < 4; i++ )
611 {
612 for ( int j = 0; j < 4; j++ )
613 { temp_PDF += t2D[i][j] * t1KsPip1[i] * t1KmPip2[j] * G[i][i] * G[j][j]; }
614 }
615 B[2] = barrier( 2, sD, sKsPip1, sKmPip2, rD2 );
616 tmp1 = barr_kspi1 * barr_kpi2 * B[2] * temp_PDF;
617
618 temp_PDF = 0;
619 if ( ispro_kspi2 == 0 )
620 {
621 propagatorRBW( mKstp, GKstp, sKsPip2, sKs, sPip2, rRes2, 1, pro_kspi2 );
622 ispro_kspi2 = 1;
623 }
624 if ( ispro_kpi1 == 0 )
625 {
626 propagatorRBW( mKst0, GKst0, sKmPip1, sKm, sPip1, rRes2, 1, pro_kpi1 );
627 ispro_kpi1 = 1;
628 }
629 Com_Multi( pro_kspi2, pro_kpi1, pro2b );
630
631 if ( barr_kspi2 < 0.0 ) barr_kspi2 = barrier( 1, sKsPip2, sKs, sPip2, rRes2 );
632 if ( barr_kpi1 < 0.0 ) barr_kpi1 = barrier( 1, sKmPip1, sKm, sPip1, rRes2 );
633 calt2( pKsPip2, pKmPip1, t2D );
634 for ( int i = 0; i < 4; i++ )
635 {
636 for ( int j = 0; j < 4; j++ )
637 { temp_PDF += t2D[i][j] * t1KsPip2[i] * t1KmPip1[j] * G[i][i] * G[j][j]; }
638 }
639 B[2] = barrier( 2, sD, sKsPip2, sKmPip1, rD2 );
640 tmp2 = barr_kspi2 * barr_kpi1 * B[2] * temp_PDF;
641 }
642
643 /// eta1475->a0 pip
644 if ( modetype[i] == 8 )
645 {
646 propagatorRBW( meta1475, Geta1475, sKsKmPip1, sKsKm, sPip1, rRes2, 0, pro1 );
647 if ( ispro_ksk == 0 )
648 {
649 propagatora0980( ma0_980, sKsKm, sKs2, sKm2, pro_ksk );
650 ispro_ksk = 1;
651 }
652 // printf("%.15f, %.15f , %.15f, %.15f \n", pro1[0],pro_ksk[0],pro1[1],pro_ksk[1]);
653 Com_Multi( pro1, pro_ksk, pro1b );
654 tmp1 = 1.0;
655
656 temp_PDF = 0;
657 propagatorRBW( meta1475, Geta1475, sKsKmPip2, sKsKm, sPip2, rRes2, 0, pro1 );
658 if ( ispro_ksk == 0 )
659 {
660 propagatora0980( ma0_980, sKsKm, sKs2, sKm2, pro_ksk );
661 ispro_ksk = 1;
662 }
663 // printf("%.15f, %.15f , %.15f, %.15f \n", pro1[0],pro_ksk[0],pro1[1],pro_ksk[1]);
664 Com_Multi( pro1, pro_ksk, pro2b );
665 tmp2 = 1.0;
666 }
667
668 // Ds -> eta(1475)pi+ eta(1475)->K*Ks,K*0>Km pi+
669 if ( modetype[i] == 4 )
670 {
671 propagatorRBW( meta1475, Geta1475, sKsKmPip1, sKmPip1, sKs, rRes2, 1, pro1 );
672 if ( ispro_kpi1 == 0 )
673 {
674 propagatorRBW( mKst0, GKst0, sKmPip1, sKm, sPip1, rRes2, 1, pro_kpi1 );
675 ispro_kpi1 = 1;
676 }
677 Com_Multi( pro1, pro_kpi1, pro1b );
678
679 if ( barr_kpi1 < 0.0 ) barr_kpi1 = barrier( 1, sKmPip1, sKm, sPip1, rRes2 );
680 if ( barr_kskpi1_kpi < 0.0 )
681 barr_kskpi1_kpi = barrier( 1, sKsKmPip1, sKmPip1, sKs, rRes2 );
682 for ( int a = 0; a < 4; a++ ) { temp_PDF += Ks[a] * t1KmPip1[a] * G[a][a]; }
683 tmp1 = barr_kpi1 * barr_kskpi1_kpi * temp_PDF;
684
685 temp_PDF = 0;
686 propagatorRBW( meta1475, Geta1475, sKsKmPip2, sKmPip2, sKs, rRes2, 1, pro1 );
687 if ( ispro_kpi2 == 0 )
688 {
689 propagatorRBW( mKst0, GKst0, sKmPip2, sKm, sPip2, rRes2, 1, pro_kpi2 );
690 ispro_kpi2 = 1;
691 }
692 Com_Multi( pro1, pro_kpi2, pro2b );
693
694 if ( barr_kpi2 < 0.0 ) barr_kpi2 = barrier( 1, sKmPip2, sKm, sPip2, rRes2 );
695 if ( barr_kskpi2_kpi < 0.0 )
696 barr_kskpi2_kpi = barrier( 1, sKsKmPip2, sKmPip2, sKs, rRes2 );
697 for ( int a = 0; a < 4; a++ ) { temp_PDF += Ks[a] * t1KmPip2[a] * G[a][a]; }
698 tmp2 = barr_kpi2 * barr_kskpi2_kpi * temp_PDF;
699 }
700
701 // Ds -> eta(1475)pi+ eta(1475)->K*Km,K*+>Ks pi
702 if ( modetype[i] == 14 )
703 {
704 propagatorRBW( meta1475, Geta1475, sKsKmPip1, sKsPip1, sKm, rRes2, 1, pro1 );
705 if ( ispro_kspi1 == 0 )
706 {
707 propagatorRBW( mKstp, GKstp, sKsPip1, sKs, sPip1, rRes2, 1, pro_kspi1 );
708 ispro_kspi1 = 1;
709 }
710 Com_Multi( pro1, pro_kspi1, pro1b );
711
712 if ( barr_kspi1 < 0.0 ) barr_kspi1 = barrier( 1, sKsPip1, sKs, sPip1, rRes2 );
713 if ( barr_kskpi1_kspi < 0.0 )
714 barr_kskpi1_kspi = barrier( 1, sKsKmPip1, sKsPip1, sKm, rRes2 );
715 for ( int a = 0; a < 4; a++ ) { temp_PDF += Km[a] * t1KsPip1[a] * G[a][a]; }
716 tmp1 = barr_kspi1 * barr_kskpi1_kspi * temp_PDF;
717
718 temp_PDF = 0;
719 propagatorRBW( meta1475, Geta1475, sKsKmPip2, sKsPip2, sKm, rRes2, 1, pro1 );
720 if ( ispro_kspi2 == 0 )
721 {
722 propagatorRBW( mKstp, GKstp, sKsPip2, sKs, sPip2, rRes2, 1, pro_kspi2 );
723 ispro_kspi2 = 1;
724 }
725 Com_Multi( pro1, pro_kspi2, pro2b );
726
727 if ( barr_kspi2 < 0.0 ) barr_kspi2 = barrier( 1, sKsPip2, sKs, sPip2, rRes2 );
728 if ( barr_kskpi2_kspi < 0.0 )
729 barr_kskpi2_kspi = barrier( 1, sKsKmPip2, sKsPip2, sKm, rRes2 );
730 for ( int a = 0; a < 4; a++ ) { temp_PDF += Km[a] * t1KsPip2[a] * G[a][a]; }
731 tmp2 = barr_kspi2 * barr_kskpi2_kspi * temp_PDF;
732 }
733
734 // f1285->a0 pip
735 if ( modetype[i] == 17 )
736 {
737 propagatorRBW( mf1285, Gf1285, sKsKmPip1, sKsKm, sPip1, rRes2, 1, pro1 );
738 if ( ispro_ksk == 0 )
739 {
740 propagatora0980( ma0_980, sKsKm, sKs2, sKm2, pro_ksk );
741 ispro_ksk = 1;
742 }
743 Com_Multi( pro1, pro_ksk, pro1b );
744 calt1( pKsKmPip1, Pip2, t1D );
745
746 if ( barr_kskpi1_ksk < 0.0 )
747 barr_kskpi1_ksk = barrier( 1, sKsKmPip1, sKsKm, sPip1, rRes2 );
748 if ( barr_ds_kskpi1 < 0.0 ) barr_ds_kskpi1 = barrier( 1, sD, sKsKmPip1, sPip2, rD2 );
749 for ( int i = 0; i < 4; i++ ) { temp_PDF += G[i][i] * t1D[i] * t1KsKmPip1_ksk[i]; }
750 tmp1 = barr_kskpi1_ksk * barr_ds_kskpi1 * temp_PDF;
751
752 temp_PDF = 0;
753 propagatorRBW( mf1285, Gf1285, sKsKmPip2, sKsKm, sPip2, rRes2, 1, pro1 );
754 if ( ispro_ksk == 0 )
755 {
756 propagatora0980( ma0_980, sKsKm, sKs2, sKm2, pro_ksk );
757 ispro_ksk = 1;
758 }
759 Com_Multi( pro1, pro_ksk, pro2b );
760 calt1( pKsKmPip2, Pip, t1D );
761
762 if ( barr_kskpi2_ksk < 0.0 )
763 barr_kskpi2_ksk = barrier( 1, sKsKmPip2, sKsKm, sPip2, rRes2 );
764 if ( barr_ds_kskpi2 < 0.0 ) barr_ds_kskpi2 = barrier( 1, sD, sKsKmPip2, sPip1, rD2 );
765 for ( int i = 0; i < 4; i++ ) { temp_PDF += G[i][i] * t1D[i] * t1KsKmPip2_ksk[i]; }
766 tmp2 = barr_kskpi2_ksk * barr_ds_kskpi2 * temp_PDF;
767 }
768
769 // Ds -> 1550 pi+ , 1550 -> K-K*+ ,K*+ ->Ks pi+
770 if ( modetype[i] == 5 )
771 {
772 // propagatorRBW(m1510,G1510,sKsKmPip1,sKsPip1,sKm,rRes2,0,pro1);
773 pro1[0] = 1;
774 pro1[1] = 0;
775 if ( ispro_kspi1 == 0 )
776 {
777 propagatorRBW( mKstp, GKstp, sKsPip1, sKs, sPip1, rRes2, 1, pro_kspi1 );
778 ispro_kspi1 = 1;
779 }
780 // printf("%.15f, %.15f , %.15f, %.15f \n",
781 // pro1[0],pro_kspi1[0],pro1[1],pro_kspi1[1]);
782 Com_Multi( pro1, pro_kspi1, pro1b );
783 calt1( pKsKmPip1, Pip2, t1D );
784
785 if ( barr_kspi1 < 0.0 ) barr_kspi1 = barrier( 1, sKsPip1, sKs, sPip1, rRes2 );
786 if ( barr_ds_kskpi1 < 0.0 ) barr_ds_kskpi1 = barrier( 1, sD, sKsKmPip1, sPip2, rD2 );
787 barr_kskpi1_kspi = barrier( 2, sKsKmPip1, sKsPip1, sKm, rRes2 );
788 for ( int i = 0; i < 4; i++ )
789 {
790 for ( int j = 0; j < 4; j++ )
791 {
792 temp_PDF += t1D[i] * ( pKsKmPip1[i] * pKsKmPip1[j] / sKsKmPip1 - G[i][j] ) *
793 t1KsPip1[j] * G[i][i] * G[j][j];
794 }
795 }
796 tmp1 = barr_kspi1 * barr_ds_kskpi1 * temp_PDF;
797
798 temp_PDF = 0;
799 // propagatorRBW(m1510,G1510,sKsKmPip2,sKsPip2,sKm,rRes2,0,pro1);
800 pro1[0] = 1;
801 pro1[1] = 0;
802 if ( ispro_kspi2 == 0 )
803 {
804 propagatorRBW( mKstp, GKstp, sKsPip2, sKs, sPip2, rRes2, 1, pro_kspi2 );
805 ispro_kspi2 = 1;
806 }
807 Com_Multi( pro1, pro_kspi2, pro2b );
808 // printf("%.15f, %.15f , %.15f, %.15f \n",
809 // pro1[0],pro_kspi2[0],pro1[1],pro_kspi2[1]);
810
811 calt1( pKsKmPip2, Pip, t1D );
812 if ( barr_kspi2 < 0.0 ) barr_kspi2 = barrier( 1, sKsPip2, sKs, sPip2, rRes2 );
813 if ( barr_ds_kskpi2 < 0.0 ) barr_ds_kskpi2 = barrier( 1, sD, sKsKmPip2, sPip1, rD2 );
814 barr_kskpi2_kspi = barrier( 2, sKsKmPip2, sKsPip2, sKm, rRes2 );
815 for ( int i = 0; i < 4; i++ )
816 {
817 for ( int j = 0; j < 4; j++ )
818 {
819 temp_PDF += t1D[i] * ( pKsKmPip2[i] * pKsKmPip2[j] / sKsKmPip2 - G[i][j] ) *
820 t1KsPip2[j] * G[i][i] * G[j][j];
821 }
822 }
823 tmp2 = barr_kspi2 * barr_ds_kskpi2 * temp_PDF;
824 }
825
826 /// K*+ (Kpi)s_wave
827 if ( modetype[i] == 21 )
828 {
829 if ( ispro_kspi1 == 0 )
830 {
831 propagatorRBW( mKstp, GKstp, sKsPip1, sKs, sPip1, rRes2, 1, pro_kspi1 );
832 ispro_kspi1 = 1;
833 }
834 KPiSLASS( sKmPip2, sKm, sPip2, pro_kpi2sw );
835 Com_Multi( pro_kspi1, pro_kpi2sw, pro1b );
836 calt1( pKsPip1, pKmPip2, t1D );
837
838 if ( barr_kspi1 < 0.0 ) barr_kspi1 = barrier( 1, sKsPip1, sKs, sPip1, rRes2 );
839 if ( barr_kspi1_kpi2 < 0.0 ) barr_kspi1_kpi2 = barrier( 1, sD, sKsPip1, sKmPip2, rD2 );
840 for ( int a = 0; a < 4; a++ ) { temp_PDF += t1D[a] * t1KsPip1[a] * G[a][a]; }
841 tmp1 = barr_kspi1 * barr_kspi1_kpi2 * temp_PDF;
842
843 temp_PDF = 0;
844 if ( ispro_kspi2 == 0 )
845 {
846 propagatorRBW( mKstp, GKstp, sKsPip2, sKs, sPip2, rRes2, 1, pro_kspi2 );
847 ispro_kspi2 = 1;
848 }
849 KPiSLASS( sKmPip1, sKm, sPip1, pro_kpi1sw );
850 Com_Multi( pro_kspi2, pro_kpi1sw, pro2b );
851
852 calt1( pKsPip2, pKmPip1, t1D );
853 if ( barr_kspi2 < 0.0 ) barr_kspi2 = barrier( 1, sKsPip2, sKs, sPip2, rRes2 );
854 if ( barr_kspi2_kpi1 < 0.0 ) barr_kspi2_kpi1 = barrier( 1, sD, sKsPip2, sKmPip1, rD2 );
855 for ( int a = 0; a < 4; a++ ) { temp_PDF += t1D[a] * t1KsPip2[a] * G[a][a]; }
856 tmp2 = barr_kspi2 * barr_kspi2_kpi1 * temp_PDF;
857 }
858
859 /// Kst0 Kspi_S_wave
860 if ( modetype[i] == 22 )
861 {
862 if ( ispro_kpi2 == 0 )
863 {
864 propagatorRBW( mKst0, GKst0, sKmPip2, sKm, sPip2, rRes2, 1, pro_kpi2 );
865 ispro_kpi2 = 1;
866 }
867 KPiSLASS( sKsPip1, sKs, sPip1, pro_kspi1sw );
868 Com_Multi( pro_kpi2, pro_kspi1sw, pro1b );
869
870 calt1( pKsPip1, pKmPip2, t1D );
871 if ( barr_kpi2 < 0.0 ) barr_kpi2 = barrier( 1, sKmPip2, sKm, sPip2, rRes2 );
872 if ( barr_kspi1_kpi2 < 0.0 ) barr_kspi1_kpi2 = barrier( 1, sD, sKsPip1, sKmPip2, rD2 );
873 for ( int a = 0; a < 4; a++ ) { temp_PDF += t1D[a] * t1KmPip2[a] * G[a][a]; }
874 tmp1 = barr_kpi2 * barr_kspi1_kpi2 * temp_PDF;
875
876 temp_PDF = 0;
877 if ( ispro_kpi1 == 0 )
878 {
879 propagatorRBW( mKst0, GKst0, sKmPip1, sKm, sPip1, rRes2, 1, pro_kpi1 );
880 ispro_kpi1 = 1;
881 }
882 KPiSLASS( sKsPip2, sKs, sPip2, pro_kspi2sw );
883 Com_Multi( pro_kpi1, pro_kspi2sw, pro2b );
884
885 calt1( pKsPip2, pKmPip1, t1D );
886 if ( barr_kpi1 < 0.0 ) barr_kpi1 = barrier( 1, sKmPip1, sKm, sPip1, rRes2 );
887 if ( barr_kspi2_kpi1 < 0.0 ) barr_kspi2_kpi1 = barrier( 1, sD, sKsPip2, sKmPip1, rD2 );
888 for ( int a = 0; a < 4; a++ ) { temp_PDF += t1D[a] * t1KmPip1[a] * G[a][a]; }
889 tmp2 = barr_kspi2_kpi1 * barr_kpi1 * temp_PDF;
890 }
891
892 if ( modetype[i] == 20 )
893 {
894 propagatorRBW( meta1475, Geta1475, sKsKmPip1, sKsPip1, sKm, rRes2, 0, pro1 );
895 KPiSLASS( sKsPip1, sKs, sPip1, pro_kspi1sw );
896 Com_Multi( pro1, pro_kspi1sw, pro1b );
897 tmp1 = 1.0;
898
899 propagatorRBW( meta1475, Geta1475, sKsKmPip2, sKsPip2, sKm, rRes2, 0, pro1 );
900 KPiSLASS( sKsPip2, sKs, sPip2, pro_kspi2sw );
901 Com_Multi( pro1, pro_kspi2sw, pro2b );
902 tmp2 = 1.0;
903 }
904
905 if ( modetype[i] == 1 || modetype[i] == 2 || modetype[i] == 3 || modetype[i] == 8 ||
906 modetype[i] == 17 || modetype[i] == 5 || modetype[i] == 21 || modetype[i] == 22 ||
907 modetype[i] == 4 || modetype[i] == 14 || modetype[i] == 20 )
908 {
909 amp_tmp[0] = tmp1 * pro1b[0] + tmp2 * pro2b[0];
910 amp_tmp[1] = tmp1 * pro1b[1] + tmp2 * pro2b[1];
911 }
912
913 /* if(modetype[i]==20){
914 amp_tmp[0] = tmp1*pro1b[0] + tmp2*pro2b[0];
915 amp_tmp[1] = tmp1*pro1b[1] + tmp2*pro2b[1];
916 } */
917 // printf("pro1b[0],pro2b[0]:: %.15f, %.15f\n", pro1b[0],pro2b[0]);
918 // printf("tmp1,tmp2:: %.15f, %.15f\n", tmp1,tmp2);
919 // printf("amp_tmp[0],amp_tmp[1]:: %.15f, %.15f\n", amp_tmp[0],amp_tmp[1]);
920
921 Com_Multi( amp_tmp, cof, amp_PDF );
922 PDF[0] += amp_PDF[0];
923 PDF[1] += amp_PDF[1];
924
925 // printf("cof= %.15f,%.15f \n",cof[0],cof[1]);
926 // printf("PDF= %.15f,%.15f \n",PDF[0],PDF[1]);
927 // printf("---------------------------------------------------\n");
928 }
929 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
930 if ( value <= 0 ) { value = 1e-20; }
931 Result = value;
932 // printf("value = %.15f\n",value);
933 // printf("---------------------------------------------------\n");
934}
mass
double mass
Definition CosmicGenerator.cxx:128
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
EvtDsToKSKmPipPip::getName
void getName(std::string &name)
Definition EvtDsToKSKmPipPip.cc:35
EvtDsToKSKmPipPip::decay
void decay(EvtParticle *p)
Definition EvtDsToKSKmPipPip.cc:162
EvtDsToKSKmPipPip::clone
EvtDecayBase * clone()
Definition EvtDsToKSKmPipPip.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