BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToKKpi.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: EvtDsToKKpi.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 "EvtDsToKKpi.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
33EvtDsToKKpi::~EvtDsToKKpi() {}
34
35void EvtDsToKKpi::getName( std::string& model_name ) { model_name = "DsToKKpi"; }
36
37EvtDecayBase* EvtDsToKKpi::clone() { return new EvtDsToKKpi; }
38
39void EvtDsToKKpi::init() {
40 // check that there are 0 arguments
41 checkNArg( 0 );
42 checkNDaug( 3 );
43 checkSpinParent( EvtSpinType::SCALAR );
44 checkSpinDaughter( 0, EvtSpinType::SCALAR );
45 checkSpinDaughter( 1, EvtSpinType::SCALAR );
46 checkSpinDaughter( 2, EvtSpinType::SCALAR );
47
48 phi[1] = 6.1944e+00;
49 phi[2] = 4.7398e+00;
50 phi[3] = 2.9047e+00;
51 phi[4] = 1.0068e+00;
52 phi[5] = 5.8035e-01;
53
54 rho[1] = 1.0963e+00;
55 rho[2] = 2.7818e+00;
56 rho[3] = 1.2570e+00;
57 rho[4] = 7.7351e-01;
58 rho[5] = 5.6277e-01;
59
60 mass[0] = 8.9581e-01;
61 mass[1] = 1.019461e+00;
62 mass[2] = 0.919;
63 mass[3] = 1.4712e+00;
64 mass[4] = 1.7220e+00;
65 mass[5] = 1.3500e+00;
66
67 width[0] = 4.7400e-02;
68 width[1] = 4.2660e-03;
69 width[2] = 0.272;
70 width[3] = 2.7000e-01;
71 width[4] = 1.3500e-01;
72 width[5] = 2.6500e-01;
73
74 modetype[0] = 0;
75 modetype[1] = 1;
76 modetype[2] = 13;
77 modetype[3] = 3;
78 modetype[4] = 4;
79 modetype[5] = 5;
80
81 phi[0] = 0;
82 rho[0] = 1;
83
84 // cout << "DsToKKpi :"<< endl;
85 // for (int i=0; i<6; i++) {
86 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
87 // }
88
89 mD = 1.86486;
90 mDs = 1.9683;
91 rRes = 1.5;
92 rD = 3.0;
93 mkstr = 0.89594;
94 mk0 = 0.497614;
95 mass_Kaon = 0.49368;
96 mass_Pion = 0.13957;
97 mass_Pi0 = 0.1349766;
98 mass_EtaP = 0.95778;
99 mass_Eta = 0.547862;
100 math_pi = 3.1415926;
101 afRatio = 2.04835; // BABAR // afRatio= 1.23202; //BES2
102 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
103 for ( int i = 0; i < 4; i++ )
104 {
105 for ( int j = 0; j < 4; j++ ) { G[i][j] = GG[i][j]; }
106 }
107}
108
109void EvtDsToKKpi::initProbMax() { setProbMax( 35638.0 ); }
110
111void EvtDsToKKpi::decay( EvtParticle* p ) {
112 /*
113 double maxprob = 0.0;
114 cout<<" calculate maxprob for DsToKKpi: "<<endl;
115 for(int ir=0;ir<=60000000;ir++){
116 p->initializePhaseSpace(getNDaug(),getDaugs());
117 EvtVector4R D1 = p->getDaug(0)->getP4();
118 EvtVector4R D2 = p->getDaug(1)->getP4();
119 EvtVector4R D3 = p->getDaug(2)->getP4();
120 double P1[4], P2[4], P3[4];
121 P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
122 P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
123 P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
124 double value;
125 int g0[6]={1,1,1,1,1,1};
126 int nstates=6;
127 calEvaMy( P1, P2, P3, mass, width, rho, phi, g0, modetype, nstates, value);
128 if(value>maxprob) {
129 maxprob=value;
130 cout << "Max PDF = " << ir << " " << maxprob << endl;
131 }
132 }
133 printf("MAXprob = %.10f\n",maxprob);
134 */
135 p->initializePhaseSpace( getNDaug(), getDaugs() );
136 EvtVector4R D1 = p->getDaug( 0 )->getP4();
137 EvtVector4R D2 = p->getDaug( 1 )->getP4();
138 EvtVector4R D3 = p->getDaug( 2 )->getP4();
139 double P1[4], P2[4], P3[4];
140 P1[1] = D1.get( 1 );
141 P1[2] = D1.get( 2 );
142 P1[3] = D1.get( 3 );
143 P2[1] = D2.get( 1 );
144 P2[2] = D2.get( 2 );
145 P2[3] = D2.get( 3 );
146 P3[1] = D3.get( 1 );
147 P3[2] = D3.get( 2 );
148 P3[3] = D3.get( 3 );
149
150 double value;
151 int g0[6] = { 1, 1, 1, 1, 1, 1 };
152 int nstates = 6;
153 calEvaMy( P1, P2, P3, mass, width, rho, phi, g0, modetype, nstates, value );
154 setProb( value );
155 return;
156}
157
158void EvtDsToKKpi::MIP_LineShape( double sa, double pro[2] ) {
159 double m0 = mass[2], cKK = width[2];
160 pro[0] = sqrt( 1 / ( ( ( m0 * m0 ) - sa ) * ( ( m0 * m0 ) - sa ) +
161 ( cKK * m0 * sqrt( 1 - 4 * mass_Kaon * mass_Kaon / sa ) ) *
162 ( cKK * m0 * sqrt( 1 - 4 * mass_Kaon * mass_Kaon / sa ) ) ) );
163 pro[1] = 0;
164}
165void EvtDsToKKpi::calEvaMy( double* pKm, double* pKp, double* pPi, double* mass1,
166 double* width1, double* amp, double* phase, int* g0, int* modetype,
167 int nstates, double& Result ) {
168 double pF0_980[4], pPhi1020[4], pKstr[4], pD[4], p23[4];
169 pKp[0] = sqrt( mass_Kaon * mass_Kaon + pKp[1] * pKp[1] + pKp[2] * pKp[2] + pKp[3] * pKp[3] );
170 pKm[0] = sqrt( mass_Kaon * mass_Kaon + pKm[1] * pKm[1] + pKm[2] * pKm[2] + pKm[3] * pKm[3] );
171 pPi[0] = sqrt( mass_Pion * mass_Pion + pPi[1] * pPi[1] + pPi[2] * pPi[2] + pPi[3] * pPi[3] );
172 for ( int u = 0; u < 4; u++ )
173 {
174 pF0_980[u] = pKp[u] + pKm[u];
175 pPhi1020[u] = pKp[u] + pKm[u];
176 pKstr[u] = pKm[u] + pPi[u];
177 p23[u] = pKp[u] + pPi[u];
178 pD[u] = pKp[u] + pKm[u] + pPi[u];
179 }
180 double cof[2], amp_tmp1[2], amp_tmp2[2], amp_tmp[2], amp_PDF[2], PDF[2];
181 double rrho = 3.0;
182 double ra1 = 3.0;
183 //----------------------------------------------------
184 amp_PDF[0] = 0;
185 amp_PDF[1] = 0;
186 PDF[0] = 0;
187 PDF[1] = 0;
188 double temp_PDF, tmp1, tmp2, temp[2];
189 double pro[2], pro0[2], pro1[2];
190 double t1D[4], t2D[4][4], B[3], t1Tm[4];
191 double t1kstr1[4], t1phi1020[4], t1D1[4], t1D2[4];
192
193 double sD, sF0_980, sF0_1710, sF0_1370, sKp, sKm, sPi, sKstr, sPhi1020, sKstr1430;
194 sF0_980 = SCADot( pF0_980, pF0_980 );
195 sF0_1710 = sF0_980;
196 sF0_1370 = sF0_980;
197 sKstr = SCADot( pKstr, pKstr );
198 sD = SCADot( pD, pD );
199 sKstr1430 = sKstr;
200 sPhi1020 = SCADot( pPhi1020, pPhi1020 );
201 sKp = SCADot( pKp, pKp );
202 sKm = SCADot( pKm, pKm );
203 sPi = SCADot( pPi, pPi );
204
205 calt1( pKm, pPi, t1kstr1 );
206 calt1( pKm, pKp, t1phi1020 );
207 calt1( pKstr, pKp, t1D1 );
208 calt1( pPhi1020, pPi, t1D2 );
209
210 double mDs = sqrt( sD );
211 for ( int i = 0; i < nstates; i++ )
212 {
213 int d = 0;
214 amp_tmp[0] = 0;
215 amp_tmp[1] = 0;
216 amp_tmp1[0] = 0;
217 amp_tmp1[1] = 0;
218 amp_tmp2[0] = 0;
219 amp_tmp2[1] = 0;
220 tmp1 = 0;
221 tmp2 = 0;
222 temp_PDF = 0;
223 cof[0] = amp[i] * cos( phase[i] );
224 cof[1] = amp[i] * sin( phase[i] );
225 if ( modetype[i] == 13 )
226 {
227 // a0(980) and f0(980) mixture
228 double amp_a0[2] = { 0 };
229 double sa0_980 = sF0_980;
230 double sKp2[2] = { sKp, mass_Pion * mass_Pion };
231 double sKma2[2] = { sKm, mass_Eta * mass_Eta };
232 // propagatora0980(mass1[i],sa0_980,sKp2,sKma2,pro);
233 MIP_LineShape( sa0_980, pro );
234 B[0] = barrier( 0, sa0_980, sKp, sKm, rRes );
235 temp_PDF = 1;
236 tmp1 = temp_PDF * B[0];
237 amp_a0[0] = tmp1 * pro[0];
238 amp_a0[1] = tmp1 * pro[1];
239 amp_tmp1[0] = amp_a0[0];
240 amp_tmp1[1] = amp_a0[1];
241 }
242 else if ( modetype[i] == 0 )
243 {
244 // K*(892) K+
245 if ( g0[i] == 0 )
246 {
247 pro[0] = 1;
248 pro[1] = 0;
249 }
250 bool neo = true;
251 if ( neo )
252 {
253 double sBC = SCADot( p23, p23 );
254 // if(g0[i] == 1) propagatorRBW(mass1[i],width1[i],sKstr,sKm,sPi,rRes,1,pro);
255 if ( g0[i] == 1 )
256 propagatorRBWNeoKstr892( mass1[i], width1[i], sKstr, sPi, sKm, rRes, 1, pro );
257 B[0] = barrierNeo( 1, sKstr, sPi, sKm, rRes, mass1[i] );
258 B[1] = barrierNeoDs( 1, sD, sKstr, sKp, rD, mDs, mass1[i] );
259 temp_PDF = ( sBC - sPhi1020 + ( ( sD - sKp ) * ( sKm - sPi ) / ( sKstr ) ) );
260 }
261 else
262 {
263 if ( g0[i] == 1 ) propagatorRBW( mass1[i], width1[i], sKstr, sKm, sPi, rRes, 1, pro );
264 // Com_Multi(pro0,pro1,pro);
265 B[0] = barrier( 1, sKstr, sKm, sPi, rRes );
266 B[1] = barrier( 1, sD, sKstr, sKp, rD );
267 for ( int m = 0; m < 4; m++ )
268 {
269 for ( int j = 0; j < 4; j++ ) { temp_PDF += t1D1[m] * t1kstr1[j] * G[m][j]; }
270 }
271 }
272 tmp1 = temp_PDF * B[0] * B[1];
273 amp_tmp1[0] = tmp1 * pro[0];
274 amp_tmp1[1] = tmp1 * pro[1];
275 }
276 else if ( modetype[i] == 1 )
277 {
278 // Phi1020 Pi v
279 if ( g0[i] == 0 )
280 {
281 pro[0] = 1;
282 pro[1] = 0;
283 }
284 /*if(g0[i] == 1) propagatorRBW(mass1[i],width1[i],sPhi1020,sKm,sKp,rRes,1,pro);
285 B[0]=barrier(1,sPhi1020,sKp,sKm,rRes);
286 B[1]=barrier(1,sD,sPhi1020,sPi,rD);
287 for(int i=0; i<4; i++){
288 for(int j=0; j<4; j++){
289 temp_PDF += t1D2[i]*t1phi1020[j]*G[i][j];
290 }
291 }
292 tmp1 = temp_PDF*B[0]*B[1];*/
293 bool neo = true;
294 if ( neo )
295 {
296 if ( g0[i] == 1 )
297 propagatorRBWNeo( mass1[i], width1[i], sPhi1020, sKm, sKp, rRes, 1, pro );
298 B[0] = barrierNeo( 1, sPhi1020, sKp, sKm, rRes, mass1[i] );
299 B[1] = barrierNeoDs( 1, sD, sPhi1020, sPi, rD, mDs, mass1[i] );
300 double sBC = SCADot( p23, p23 );
301 temp_PDF = ( sKstr - sBC + ( ( sD - sPi ) * ( sKp - sKm ) / ( sKstr ) ) );
302 }
303 else
304 {
305 if ( g0[i] == 1 )
306 propagatorRBW( mass1[i], width1[i], sPhi1020, sKm, sKp, rRes, 1, pro );
307 B[0] = barrier( 1, sPhi1020, sKp, sKm, rRes );
308 B[1] = barrier( 1, sD, sPhi1020, sPi, rD );
309 for ( int m = 0; m < 4; m++ )
310 {
311 for ( int j = 0; j < 4; j++ ) { temp_PDF += t1D2[m] * t1phi1020[j] * G[m][j]; }
312 }
313 }
314 // Com_Multi(pro0,pro1,pro);
315
316 tmp1 = temp_PDF * B[0] * B[1];
317 amp_tmp1[0] = tmp1 * pro[0];
318 amp_tmp1[1] = tmp1 * pro[1];
319 }
320 else if ( modetype[i] == 3 )
321 {
322 // Kstr1430 K S
323 /**
324 * Normal Model
325 * */
326 // if(g0[i] == 1) propagatorRBW(mass1[i],width1[i],sKstr1430,sKm,sPi,rRes,0,pro);
327 // propagatorRBWNeo(mass1[i],width1[i],sKstr1430,sKm,sPi,rRes,0,pro);
328 double sKm2[2] = { sKm, mass_EtaP * mass_EtaP };
329 double sPi2[2] = { sPi, mass_Kaon * mass_Kaon };
330 propagatorKstr1430( mass1[i], sKstr1430, sKm2, sPi2, pro );
331 // Com_Multi(pro0,pro1,pro);
332 B[0] = barrier( 0, sPhi1020, sKp, sKm, rRes );
333 tmp1 = 1 * B[0];
334 amp_tmp1[0] = tmp1 * pro[0];
335 amp_tmp1[1] = tmp1 * pro[1];
336
337 /**
338 * KPi-S wave LASS
339 * */
340 /*amp_tmp1[0]=Amp_KPiS1[0];
341 amp_tmp1[1]=Amp_KPiS1[1];*/
342 }
343 else if ( modetype[i] == 4 )
344 {
345 // f0(1710) Pi+ S
346 if ( g0[i] == 1 )
347 propagatorRBWNeo( mass1[i], width1[i], sF0_1710, sKp, sKm, rRes, 0, pro );
348 // if(g0[i] == 1) propagatorRBW(mass1[i],width1[i],sF0_1710,sKp,sKm,rRes,0,pro);
349 if ( g0[i] == 0 )
350 {
351 pro[0] = 1;
352 pro[1] = 0;
353 }
354 // Com_Multi(pro0,pro1,pro);
355 B[0] = barrier( 0, sF0_1710, sKp, sKm, rRes );
356 temp_PDF = 1;
357 tmp1 = temp_PDF * B[0];
358 amp_tmp1[0] = tmp1 * pro[0];
359 amp_tmp1[1] = tmp1 * pro[1];
360 /*amp_tmp1[0]=Amp_KPiS1[0];
361 amp_tmp1[1]=Amp_KPiS1[1];*/
362 const bool debug89 = false;
363 }
364 else if ( modetype[i] == 5 )
365 {
366 // f0(1370) Pi+ S
367 if ( g0[i] == 1 )
368 propagatorRBWNeo( mass1[i], width1[i], sF0_1370, sKp, sKm, rRes, 0, pro );
369 // if(g0[i] == 1) propagatorRBW(mass1[i],width1[i],sF0_1370,sKp,sKm,rRes,0,pro);
370 if ( g0[i] == 0 )
371 {
372 pro[0] = 1;
373 pro[1] = 0;
374 }
375 // Com_Multi(pro0,pro1,pro);
376 B[0] = barrier( 0, sF0_1370, sKp, sKm, rRes );
377 tmp1 = 1 * B[0];
378 amp_tmp1[0] = tmp1 * pro[0];
379 amp_tmp1[1] = tmp1 * pro[1];
380 }
381 amp_tmp[0] = amp_tmp1[0];
382 amp_tmp[1] = amp_tmp1[1];
383 Com_Multi( amp_tmp, cof, amp_PDF );
384 // printf("Line: %u @ %u modetype[%d]=%d amp_tmp[0]=%f amp_tmp[1]=%f cof[0]=%f cof[1]=%f
385 // \n",__LINE__,__FILE__,i,modetype[i],amp_tmp[0],amp_tmp[1],cof[0],cof[1]);
386 PDF[0] += amp_PDF[0];
387 PDF[1] += amp_PDF[1];
388 double valTmp = amp_PDF[0] * amp_PDF[0] + amp_PDF[1] * amp_PDF[1];
389 }
390
391 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
392 if ( value <= 0 ) { value = 1e-20; }
393 // printf("Line: %u @ %u value=%f\n",__LINE__,__FILE__,value);
394 Result = value;
395}
396
397void EvtDsToKKpi::Com_Multi( double a1[2], double a2[2], double res[2] ) {
398 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
399 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
400}
401void EvtDsToKKpi::Com_Divide( double a1[2], double a2[2], double res[2] ) {
402 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / ( a2[0] * a2[0] + a2[1] * a2[1] );
403 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / ( a2[0] * a2[0] + a2[1] * a2[1] );
404}
405
406double EvtDsToKKpi::SCADot( double a1[4], double a2[4] ) {
407 double _cal = 0.;
408 _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
409 return _cal;
410}
411double EvtDsToKKpi::barrier( int l, double sa, double sb, double sc, double r ) {
412 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
413 if ( q < 0 ) q = 1e-16;
414 double z;
415 z = q * r * r;
416 double F = 0.0;
417 if ( l == 0 ) F = 1;
418 if ( l == 1 ) F = sqrt( 2 * z / ( 1 + z ) );
419 if ( l == 2 ) F = sqrt( 13 * z * z / ( 9 + 3 * z + z * z ) );
420 return F;
421}
422double EvtDsToKKpi::barrierNeo( int l, double sa, double sb, double sc, double r, double mR ) {
423 double pAB = ( ( sa - sb - sc ) * ( sa - sb - sc ) / 4.0 - ( sb * sc ) ) / sa;
424 double pR =
425 ( ( mR * mR - sb - sc ) * ( mR * mR - sb - sc ) / 4.0 - ( sb * sc ) ) / ( mR * mR );
426 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
427 double zAB = pAB * r * r;
428 double zR = pR * r * r;
429 double F = 0.0;
430 if ( l == 0 ) F = 1;
431 if ( l == 1 ) F = sqrt( ( 1 + zR ) / ( 1 + zAB ) );
432 if ( l == 2 ) F = sqrt( ( 9 + 3 * zAB + zAB * zAB ) / ( 9 + 3 * zAB + zAB * zAB ) );
433 return F;
434}
435double EvtDsToKKpi::barrierNeoDs( int l, double sa, double sb, double sc, double r, double mR,
436 double mb ) {
437 double pAB = ( ( sa - sb - sc ) * ( sa - sb - sc ) / 4.0 - ( sb * sc ) ) / sa;
438 double pR =
439 ( ( sa - mb * mb - sc ) * ( sa - mb * mb - sc ) / 4.0 - ( mb * mb * sc ) ) / ( mR * mR );
440 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
441 double zAB = pAB * r * r;
442 double zR = pR * r * r;
443 double F = 0.0;
444 if ( l == 0 ) F = 1;
445 if ( l == 1 ) F = sqrt( ( 1 + zR ) / ( 1 + zAB ) );
446 if ( l == 2 ) F = sqrt( ( 9 + 3 * zAB + zAB * zAB ) / ( 9 + 3 * zAB + zAB * zAB ) );
447 return F;
448}
449//------------------spin-------------------------------------------
450void EvtDsToKKpi::calt1( double daug1[4], double daug2[4], double t1[4] ) {
451 double p, pq;
452 double pa[4], qa[4];
453 for ( int i = 0; i < 4; i++ )
454 {
455 pa[i] = daug1[i] + daug2[i];
456 qa[i] = daug1[i] - daug2[i];
457 }
458 p = SCADot( pa, pa );
459 pq = SCADot( pa, qa );
460 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - pq / p * pa[i]; }
461}
462void EvtDsToKKpi::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
463 double p, r;
464 double pa[4], t1[4];
465 calt1( daug1, daug2, t1 );
466 r = SCADot( t1, t1 );
467 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
468 p = SCADot( pa, pa );
469 for ( int i = 0; i < 4; i++ )
470 {
471 for ( int j = 0; j < 4; j++ )
472 { t2[i][j] = t1[i] * t1[j] - 1.0 / 3 * r * ( G[i][j] - pa[i] * pa[j] / p ); }
473 }
474}
475//-------------------prop--------------------------------------------
476void EvtDsToKKpi::propagator( double mass, double width, double sx, double prop[2] ) {
477 double a[2], b[2];
478 a[0] = 1;
479 a[1] = 0;
480 b[0] = mass * mass - sx;
481 b[1] = -mass * width;
482 Com_Divide( a, b, prop );
483}
484double EvtDsToKKpi::wid( double mass, double sa, double sb, double sc, double r, int l ) {
485 double widm = 0.;
486 double q = 0.;
487 double q0 = 0.;
488 double sa0 = mass * mass;
489 double m = sqrt( sa );
490 q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
491 if ( q < 0 ) q = 1e-16;
492 q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
493 if ( q0 < 0 ) q0 = 1e-16;
494 // double F = barrier(l,sa,sb,sc,r);
495 double z = q * r * r;
496 double z0 = q0 * r * r;
497 double F = 0;
498 if ( l == 0 ) F = 1;
499 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
500 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
501 double t = q / q0;
502 widm = pow( t, l + 0.5 ) * mass / m * F * F;
503 return widm;
504}
505
506void EvtDsToKKpi::Flatte_rhoab( double sa, double sb, double sc, double rho[2] ) {
507 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
508 double b[2];
509 if ( q > 0 )
510 {
511 rho[0] = 2 * sqrt( q / sa );
512 rho[1] = 0;
513 }
514 else if ( q < 0 )
515 {
516 rho[0] = 0;
517 rho[1] = 2 * sqrt( -q / sa );
518 }
519}
520
521void EvtDsToKKpi::propagatorFlatte( double mass, double width, double sx, double* sb,
522 double* sc, double prop[2] ) {
523 double unit[2] = { 1.0 };
524 double ci[2] = { 0, 1 };
525 double rho1[2];
526 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
527 double rho2[2];
528 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
529 double g1_f980 = 0.165, g2_f980 = 0.69465;
530 double tmp1[2] = { g1_f980, 0 };
531 double tmp11[2];
532 double tmp2[2] = { g2_f980, 0 };
533 double tmp22[2];
534 Com_Multi( tmp1, rho1, tmp11 );
535 Com_Multi( tmp2, rho2, tmp22 );
536 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
537 double tmp31[2];
538 Com_Multi( tmp3, ci, tmp31 );
539 double tmp4[2] = { mass * mass - sx - tmp31[0], -1.0 * tmp3[1] };
540 Com_Divide( unit, tmp4, prop );
541}
542void EvtDsToKKpi::propagator980( double mass, double sx, double* sb, double* sc,
543 double prop[2] ) {
544 double unit[2] = { 1.0 };
545 double ci[2] = { 0, 1 };
546 double rho1[2];
547 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
548 double rho2[2];
549 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
550 double gK_f980 = 0.69466, gPi_f980 = 0.165;
551 double tmp1[2] = { gK_f980, 0 };
552 double tmp11[2];
553 double tmp2[2] = { gPi_f980, 0 };
554 double tmp22[2];
555 Com_Multi( tmp1, rho1, tmp11 );
556 Com_Multi( tmp2, rho2, tmp22 );
557 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
558 double tmp31[2];
559 Com_Multi( tmp3, ci, tmp31 );
560 double tmp4[2] = { mass * mass - sx - tmp31[0], -1.0 * tmp31[1] };
561 Com_Divide( unit, tmp4, prop );
562}
563
564void EvtDsToKKpi::propagatora0980( double mass, double sx, double* sb, double* sc,
565 double prop[2] ) {
566 double unit[2] = { 1.0 };
567 double ci[2] = { 0, 1 };
568 double rho1[2];
569 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
570 double rho2[2];
571 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
572 double gKK_a980 = 0.892 * 0.341, gPiEta_a980 = 0.341;
573 double tmp1[2] = { gKK_a980, 0 };
574 double tmp11[2];
575 double tmp2[2] = { gPiEta_a980, 0 };
576 double tmp22[2];
577 Com_Multi( tmp1, rho1, tmp11 );
578 Com_Multi( tmp2, rho2, tmp22 );
579 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
580 double tmp31[2];
581 Com_Multi( tmp3, ci, tmp31 );
582 double tmp4[2] = { mass * mass - sx - tmp31[0], -1.0 * tmp31[1] };
583 Com_Divide( unit, tmp4, prop );
584}
585
586void EvtDsToKKpi::propagatorKstr1430( double mass, double sx, double* sb, double* sc,
587 double prop[2] ) {
588 double unit[2] = { 1.0 };
589 double ci[2] = { 0, 1 };
590 double rho1[2];
591 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
592 double rho2[2];
593 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
594 double gKPi_Kstr1430 = 0.2990, gEtaPK_Kstr1430 = 0.0529;
595 double tmp1[2] = { gKPi_Kstr1430, 0 };
596 double tmp11[2];
597 double tmp2[2] = { gEtaPK_Kstr1430, 0 };
598 double tmp22[2];
599 Com_Multi( tmp1, rho1, tmp11 );
600 Com_Multi( tmp2, rho2, tmp22 );
601 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
602 double tmp31[2];
603 Com_Multi( tmp3, ci, tmp31 );
604 double tmp4[2] = { mass * mass - sx - tmp31[0], -1.0 * tmp31[1] };
605 Com_Divide( unit, tmp4, prop );
606}
607
608void EvtDsToKKpi::propagatorRBW( double mass, double width, double sa, double sb, double sc,
609 double r, int l, double prop[2] ) {
610 double a[2], b[2];
611 a[0] = 1;
612 a[1] = 0;
613 b[0] = mass * mass - sa;
614 b[1] = -mass * width * wid( mass, sa, sb, sc, r, l );
615 Com_Divide( a, b, prop );
616}
617
618void EvtDsToKKpi::propagatorRBWNeo( double mass, double width, double sa, double sb, double sc,
619 double r, int l, double prop[2] ) {
620 double a[2], b[2];
621 a[0] = 1;
622 a[1] = 0;
623 b[0] = mass * mass - sa;
624 double tmp1 = ( sa - sb - sc );
625 double tmp2 = sb * sc;
626 double pAB = sqrt( ( tmp1 * tmp1 / 4.0 - tmp2 ) / sa );
627 double pR =
628 sqrt( ( ( mass * mass - sb - sc ) * ( mass * mass - sb - sc ) / 4.0 - ( sb * sc ) ) /
629 ( mass * mass ) );
630 double fR = sqrt( 1.0 + 1.5 * 1.5 * pR * pR ) / sqrt( 1.0 + 1.5 * 1.5 * pAB * pAB );
631 double power;
632 if ( !l )
633 {
634 power = 1;
635 fR = 1;
636 }
637 else if ( l == 1 ) { power = 3; }
638 double gammaAB = width * pow( pAB / pR, power ) * ( mass / sqrt( sa ) ) * fR * fR;
639 b[1] = -mass * gammaAB;
640 Com_Divide( a, b, prop );
641}
642
643void EvtDsToKKpi::propagatorRBWNeoKstr892( double mass, double width, double sa, double sb,
644 double sc, double r, int l, double prop[2] ) {
645 double a[2], b[2];
646 a[0] = 1;
647 a[1] = 0;
648 b[0] = mass * mass - sa;
649 double tmp1 = ( sa - sb - sc );
650 double tmp2 = sb * sc;
651 double pAB = sqrt( ( tmp1 * tmp1 / 4.0 - tmp2 ) / sa );
652 double pR =
653 sqrt( ( ( mass * mass - sb - sc ) * ( mass * mass - sb - sc ) / 4.0 - ( sb * sc ) ) /
654 ( mass * mass ) );
655 double fR = sqrt( 1.0 + 1.5 * 1.5 * pR * pR ) / sqrt( 1.0 + 1.5 * 1.5 * pAB * pAB );
656 double power;
657 if ( !l ) { power = 1; }
658 else if ( l == 1 ) { power = 3; }
659 double gammaAB = width * pow( pAB / pR, power ) * ( mass / sqrt( sa ) ) * fR * fR;
660 b[1] = -mass * gammaAB;
661 Com_Divide( a, b, prop );
662}
663
664//------------GS---used by rho----------------------------
665double EvtDsToKKpi::h( double m, double q ) {
666 double h = 2 / math_pi * q / m * log( ( m + 2 * q ) / ( 2 * mass_Pion ) );
667 return h;
668}
669double EvtDsToKKpi::dh( double mass, double q0 ) {
670 double dh = h( mass, q0 ) * ( 1.0 / ( 8 * q0 * q0 ) - 1.0 / ( 2 * mass * mass ) ) +
671 1.0 / ( 2 * math_pi * mass * mass );
672 return dh;
673}
674double EvtDsToKKpi::f( double mass, double sx, double q0, double q ) {
675 double m = sqrt( sx );
676 double f = mass * mass / ( pow( q0, 3 ) ) *
677 ( q * q * ( h( m, q ) - h( mass, q0 ) ) +
678 ( mass * mass - sx ) * q0 * q0 * dh( mass, q0 ) );
679 return f;
680}
681double EvtDsToKKpi::d( double mass, double q0 ) {
682 double d = 3.0 / math_pi * mass_Pion * mass_Pion / ( q0 * q0 ) *
683 log( ( mass + 2 * q0 ) / ( 2 * mass_Pion ) ) +
684 mass / ( 2 * math_pi * q0 ) -
685 ( mass_Pion * mass_Pion * mass ) / ( math_pi * pow( q0, 3 ) );
686 return d;
687}
688
689// rho
690void EvtDsToKKpi::propagatorGS( double mass, double width, double sa, double sb, double sc,
691 double r, int l, double prop[2] ) {
692 double a[2], b[2];
693 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
694 double sa0 = mass * mass;
695 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
696 if ( q < 0 ) q = 1e-16;
697 if ( q0 < 0 ) q0 = 1e-16;
698 q = sqrt( q );
699 q0 = sqrt( q0 );
700 a[0] = 1 + d( mass, q0 ) * width / mass;
701 a[1] = 0;
702 b[0] = mass * mass - sa + width * f( mass, sa, q0, q );
703 b[1] = -mass * width * wid( mass, sa, sb, sc, r, l );
704 Com_Divide( a, b, prop );
705}
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
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
EvtDsToKKpi::initProbMax
void initProbMax()
Definition EvtDsToKKpi.cc:109
EvtDsToKKpi::decay
void decay(EvtParticle *p)
Definition EvtDsToKKpi.cc:111
EvtDsToKKpi::getName
void getName(std::string &name)
Definition EvtDsToKKpi.cc:35
EvtDsToKKpi::clone
EvtDecayBase * clone()
Definition EvtDsToKKpi.cc:37
EvtDsToKKpi::~EvtDsToKKpi
virtual ~EvtDsToKKpi()
Definition EvtDsToKKpi.cc:33
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