BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDToKSpipipi.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// Environment:
3// This software is part of models developed at BES collaboration
4// based on the EvtGen framework. If you use all or part
5// of it, please give an appropriate acknowledgement.
6//
7// Copyright Information: See EvtGen/BesCopyright
8// Copyright (A) 2006 Ping Rong-Gang @IHEP
9//
10// Module: EvtDToKSpipipi.cc
11// the necessary file: EvtDToKSpipipi.hh
12//
13// Description: D+ -> KS0 pi+ pi+ pi-,
14// PHYSICAL REVIEW D 100, 072008 (2019)
15//
16// Modification history:
17//
18// Liaoyuan Dong Jan.15, 2020 Module created
19//
20//------------------------------------------------------------------------
21#include "EvtDToKSpipipi.hh"
22#include "../EvtGenBase/EvtComplex.hh"
23#include "../EvtGenBase/EvtDecayTable.hh"
24#include "../EvtGenBase/EvtGenKine.hh"
25#include "../EvtGenBase/EvtPDL.hh"
26#include "../EvtGenBase/EvtParticle.hh"
27#include "../EvtGenBase/EvtPatches.hh"
28#include "../EvtGenBase/EvtReport.hh"
29#include <stdlib.h>
30
31EvtDToKSpipipi::~EvtDToKSpipipi() {}
32
33void EvtDToKSpipipi::getName( std::string& model_name ) { model_name = "DToKSpipipi"; }
34
35EvtDecayBase* EvtDToKSpipipi::clone() { return new EvtDToKSpipipi; }
36
37void EvtDToKSpipipi::init() {
38 checkNArg( 0 );
39 checkNDaug( 4 );
40 checkSpinParent( EvtSpinType::SCALAR );
41 /*
42 checkSpinDaughter(0,EvtSpinType::SCALAR);
43 checkSpinDaughter(1,EvtSpinType::SCALAR);
44 checkSpinDaughter(3,EvtSpinType::SCALAR);
45 checkSpinDaughter(4,EvtSpinType::SCALAR);
46 */
47 // std::cout << "Initializing EvtDToKSpipipi" << std::endl;
48
49 mK1270 = 1.272;
50 mK1400 = 1.403;
51 GK1270 = 0.09;
52 GK1400 = 0.174;
53 mKstr = 0.89166;
54 mrho = 0.77549;
55 GKstr = 0.0462;
56 Grho = 0.1491;
57 msigma = 0.472;
58 Gsigma = 0.542;
59 phi_omega = -0.02;
60 mK1650 = 1.65;
61 GK1650 = 0.158;
62 rho[0] = 1.0;
63 phi[0] = 0.0;
64
65 ma1 = 1.22;
66 Ga1 = 0.4282;
67 mK1460 = 1.4152;
68 GK1460 = 0.2485;
69 rho_omega = 0.00294;
70
71 phi[1] = -1.55;
72 rho[1] = 0.5843;
73
74 phi[2] = -1.8223;
75 rho[2] = 2.0974;
76
77 phi[3] = -2.6751;
78 rho[3] = 0.46642;
79
80 phi[4] = -2.2429;
81 rho[4] = 0.33334;
82
83 phi[5] = -0.55888;
84 rho[5] = 0.15549;
85
86 phi[6] = -1.8778;
87 rho[6] = 0.94452;
88
89 phi[7] = 2.7724;
90 rho[7] = 0.99411;
91
92 phi[8] = -2.6461;
93 rho[8] = 0.21231;
94
95 phi[9] = -0.95137;
96 rho[9] = 0.29895;
97
98 phi[10] = -3.0843;
99 rho[10] = 3.6361;
100
101 phi[11] = 2.0954;
102 rho[11] = 0.96472;
103
104 phi[12] = -2.4965;
105 rho[12] = 0.48470;
106
107 // for (int i=0; i<13; i++) {
108 // cout << i << " rho,phi = " << rho[i] << ", "<< phi[i] << endl;
109 // }
110
111 mD = 1.86486;
112 rD = 5;
113 metap = 0.95778;
114 mkstr = 0.89594;
115 mk0 = 0.497614;
116 mass_Kaon = 0.49368;
117 mass_Pion = 0.13957;
118 math_pi = 3.1415926;
119
120 pi = 3.1415926;
121 mpi = 0.13957;
122 g1 = 0.5468;
123 g2 = 0.23;
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 EvtDToKSpipipi::initProbMax() { setProbMax( 3700.0 ); }
157
158void EvtDToKSpipipi::decay( EvtParticle* p ) {
159 /*
160 double maxprob = 0.0;
161 for(int ir=0;ir<=60000000;ir++){
162 p->initializePhaseSpace(getNDaug(),getDaugs());
163 EvtVector4R Ks0 = p->getDaug(0)->getP4();
164 EvtVector4R pi1 = p->getDaug(1)->getP4();
165 EvtVector4R pi2 = p->getDaug(2)->getP4();
166 EvtVector4R pi3 = p->getDaug(3)->getP4();
167 double Ks[4],Pip1[4],Pip2[4],Pim[4];
168 Ks[0] = Ks0.get(0); Pip1[0] = pi1.get(0); Pip2[0] = pi2.get(0); Pim[0] = pi3.get(0);
169 Ks[1] = Ks0.get(1); Pip1[1] = pi1.get(1); Pip2[1] = pi2.get(1); Pim[1] = pi3.get(1);
170 Ks[2] = Ks0.get(2); Pip1[2] = pi1.get(2); Pip2[2] = pi2.get(2); Pim[2] = pi3.get(2);
171 Ks[3] = Ks0.get(3); Pip1[3] = pi1.get(3); Pip2[3] = pi2.get(3); Pim[3] = pi3.get(3);
172 double Prob = calPDF(Ks, Pip1, Pip2, Pim);
173 if(Prob>maxprob) {
174 maxprob=Prob;
175 std::cout << "Max PDF = " << ir << " prob= " << Prob << std::endl;
176 }
177 }
178 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
179 */
180 p->initializePhaseSpace( getNDaug(), getDaugs() );
181 EvtVector4R Ks0 = p->getDaug( 0 )->getP4();
182 EvtVector4R pi1 = p->getDaug( 1 )->getP4();
183 EvtVector4R pi2 = p->getDaug( 2 )->getP4();
184 EvtVector4R pi3 = p->getDaug( 3 )->getP4();
185
186 double Ks[4], Pip1[4], Pip2[4], Pim[4];
187 Ks[0] = Ks0.get( 0 );
188 Pip1[0] = pi1.get( 0 );
189 Pip2[0] = pi2.get( 0 );
190 Pim[0] = pi3.get( 0 );
191 Ks[1] = Ks0.get( 1 );
192 Pip1[1] = pi1.get( 1 );
193 Pip2[1] = pi2.get( 1 );
194 Pim[1] = pi3.get( 1 );
195 Ks[2] = Ks0.get( 2 );
196 Pip1[2] = pi1.get( 2 );
197 Pip2[2] = pi2.get( 2 );
198 Pim[2] = pi3.get( 2 );
199 Ks[3] = Ks0.get( 3 );
200 Pip1[3] = pi1.get( 3 );
201 Pip2[3] = pi2.get( 3 );
202 Pim[3] = pi3.get( 3 );
203 double prob = calPDF( Ks, Pip1, Pip2, Pim );
204 setProb( prob );
205 return;
206}
207
208double EvtDToKSpipipi::calPDF( double Ks[], double Pip1[], double Pip2[], double Pim[] ) {
209 EvtComplex PDF[100];
210 double P14[4], P24[4], P34[4], P124[4], P134[4], P234[4];
211 for ( int i = 0; i < 4; i++ )
212 {
213 P14[i] = Ks[i] + Pim[i];
214 P24[i] = Pip1[i] + Pim[i];
215 P34[i] = Pip2[i] + Pim[i];
216 P124[i] = P14[i] + Pip1[i];
217 P134[i] = P14[i] + Pip2[i];
218 P234[i] = P24[i] + Pip2[i];
219 }
220 //----------D->a1Ks--------------
221 //----------a1->rhoPi------------
222 PDF[0] = D2AP_A2VP( Ks, Pip2, Pip1, Pim, 0 ) * getprop( P24, Pip2, ma1, Ga1, 1, 0 ) *
223 getprop( Pip1, Pim, mrho, Grho, 2, 1 ) +
224 D2AP_A2VP( Ks, Pip1, Pip2, Pim, 0 ) * getprop( P34, Pip1, ma1, Ga1, 1, 0 ) *
225 getprop( Pip2, Pim, mrho, Grho, 2, 1 );
226 PDF[1] = D2AP_A2VP( Ks, Pip2, Pip1, Pim, 2 ) * getprop( P24, Pip2, ma1, Ga1, 1, 2 ) *
227 getprop( Pip1, Pim, mrho, Grho, 2, 1 ) +
228 D2AP_A2VP( Ks, Pip1, Pip2, Pim, 2 ) * getprop( P34, Pip1, ma1, Ga1, 1, 2 ) *
229 getprop( Pip2, Pim, mrho, Grho, 2, 1 );
230 //----------a1->sigma pi---------
231 PDF[2] = D2AP_A2SP( Ks, Pip2, Pip1, Pim ) * getprop( P24, Pip2, ma1, Ga1, 1, 1 ) *
232 getprop( Pip1, Pim, msigma, Gsigma, 4, 0 ) +
233 D2AP_A2SP( Ks, Pip1, Pip2, Pim ) * getprop( P34, Pip1, ma1, Ga1, 1, 1 ) *
234 getprop( Pip2, Pim, msigma, Gsigma, 4, 0 );
235 //----------D->a1K finish-----
236 //---------D->K1(1400) pi-----
237 // K1400[S]->K* pi
238 PDF[3] = D2AP_A2VP( Pip2, Pip1, Ks, Pim, 0 ) * getprop( P14, Pip1, mK1400, GK1400, 1, 0 ) *
239 getprop( Ks, Pim, mKstr, GKstr, 1, 1 ) +
240 D2AP_A2VP( Pip1, Pip2, Ks, Pim, 0 ) * getprop( P14, Pip2, mK1400, GK1400, 1, 0 ) *
241 getprop( Ks, Pim, mKstr, GKstr, 1, 1 );
242 // K1400[D]->K* pi
243 PDF[4] = D2AP_A2VP( Pip2, Pip1, Ks, Pim, 2 ) * getprop( P14, Pip1, mK1400, GK1400, 1, 2 ) *
244 getprop( Ks, Pim, mKstr, GKstr, 1, 1 ) +
245 D2AP_A2VP( Pip1, Pip2, Ks, Pim, 2 ) * getprop( P14, Pip2, mK1400, GK1400, 1, 2 ) *
246 getprop( Ks, Pim, mKstr, GKstr, 1, 1 );
247 //-----------------------------
248 //-------K1270[S]->Ksrho-------
249 PDF[5] = D2AP_A2VP( Pip2, Ks, Pip1, Pim, 0 ) * getprop( P24, Ks, mK1270, GK1270, 0, 0 ) *
250 getprop( Pip1, Pim, mrho, Grho, 2, 1 ) +
251 D2AP_A2VP( Pip1, Ks, Pip2, Pim, 0 ) * getprop( P34, Ks, mK1270, GK1270, 0, 0 ) *
252 getprop( Pip2, Pim, mrho, Grho, 2, 1 );
253 //-------D->rhoKAD------------
254 PDF[6] = D2AP_A2VP( Pip2, Ks, Pip1, Pim, 0 ) * getprop( Pip1, Pim, mrho, Grho, 2, 1 ) +
255 D2AP_A2VP( Pip1, Ks, Pip2, Pim, 0 ) * getprop( Pip2, Pim, mrho, Grho, 2, 1 );
256 PDF[7] = D2AP_A2VP( Pip2, Ks, Pip1, Pim, 2 ) * getprop( Pip1, Pim, mrho, Grho, 2, 1 ) +
257 D2AP_A2VP( Pip1, Ks, Pip2, Pim, 2 ) * getprop( Pip2, Pim, mrho, Grho, 2, 1 );
258 //-------D->K1460, K1460->Ks rho---------
259 PDF[8] = D2PP_P2VP( Pip2, Ks, Pip1, Pim ) * getprop( P24, Ks, mK1460, GK1460, 1, 1 ) *
260 getprop( Pip1, Pim, mrho, Grho, 2, 1 ) +
261 D2PP_P2VP( Pip1, Ks, Pip2, Pim ) * getprop( P34, Ks, mK1460, GK1460, 1, 1 ) *
262 getprop( Pip2, Pim, mrho, Grho, 2, 1 );
263 //--------K*PiA (K1650)---------------------
264 PDF[9] = D2AP_A2VP( Pip2, Pip1, Ks, Pim, 0 ) * getprop( P14, Pip1, mK1650, GK1650, 1, 0 ) *
265 getprop( Ks, Pim, mKstr, GKstr, 1, 1 ) +
266 D2AP_A2VP( Pip1, Pip2, Ks, Pim, 0 ) * getprop( P14, Pip2, mK1650, GK1650, 1, 0 ) *
267 getprop( Ks, Pim, mKstr, GKstr, 1, 1 );
268 //-------KsPiPiSPiA-----------------
269 PDF[10] = D2AP_A2SP( Pip2, Ks, Pip1, Pim ) + D2AP_A2SP( Pip1, Ks, Pip2, Pim );
270 //-------KPiS wave------------------
271 EvtComplex corr( 2, 0 );
272 PDF[11] = corr * PHSP( Ks, Pim );
273 // cout<<PHSP(Ks,Pim)<<endl;
274 //-------D->K1460pi, K1460->K*-pi+--------
275 PDF[12] = D2PP_P2VP( Pip2, Pip1, Ks, Pim ) * getprop( P14, Pip1, mK1460, GK1460, 1, 1 ) *
276 getprop( Ks, Pim, mKstr, GKstr, 1, 1 ) +
277 D2PP_P2VP( Pip1, Pip2, Ks, Pim ) * getprop( P14, Pip2, mK1460, GK1460, 1, 1 ) *
278 getprop( Ks, Pim, mKstr, GKstr, 1, 1 );
279 //-------------------------------------------
280 EvtComplex cof;
281 EvtComplex pdf, module;
282 pdf = EvtComplex( 0, 0 );
283 for ( int i = 0; i < 13; i++ )
284 {
285 cof = EvtComplex( rho[i] * cos( phi[i] ), rho[i] * sin( phi[i] ) );
286 pdf = pdf + cof * PDF[i];
287 }
288 module = conj( pdf ) * pdf;
289 double value;
290 value = real( module );
291 return ( value <= 0 ) ? 1e-20 : value;
292}
293
294EvtComplex EvtDToKSpipipi::KPiSFormfactor( double sa, double sb, double sc, double r ) {
295 double m1430 = 1.463;
296 double sa0 = m1430 * m1430;
297 double w1430 = 0.233;
298 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
299 if ( q0 < 0 ) q0 = 1e-16;
300 double qs = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
301 double q = sqrt( qs );
302 double width = w1430 * q * m1430 / sqrt( sa * q0 );
303 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
304 if ( temp_R < 0 ) temp_R += math_pi;
305 double deltaR = -5.31 + temp_R;
306 double temp_F = atan( 2 * 1.07 * q / ( 2 + ( -1.8 ) * 1.07 * qs ) );
307 if ( temp_F < 0 ) temp_F += math_pi;
308 double deltaF = 2.33 + temp_F;
309 EvtComplex cR( cos( deltaR ), sin( deltaR ) );
310 EvtComplex cF( cos( deltaF ), sin( deltaF ) );
311 EvtComplex amp = 0.8 * sin( deltaF ) * cF + sin( deltaR ) * cR * cF * cF;
312 return amp;
313}
314EvtComplex EvtDToKSpipipi::D2AP_A2VP( double P1[], double P2[], double P3[], double P4[],
315 int L ) {
316 double temp_PDF = 0;
317 EvtComplex amp_PDF( 0, 0 );
318 EvtComplex pro[2];
319 double t1V[4], t1D[4], t2A[4][4];
320 double sa[3], sb[3], sc[3], B[3];
321 double pV[4], pA[4], pD[4];
322 for ( int i = 0; i != 4; i++ )
323 {
324 pV[i] = P3[i] + P4[i];
325 pA[i] = pV[i] + P2[i];
326 pD[i] = pA[i] + P1[i];
327 }
328 sa[0] = dot( pV, pV );
329 sb[0] = dot( P3, P3 );
330 sc[0] = dot( P4, P4 );
331 sa[1] = dot( pA, pA );
332 sb[1] = sa[0];
333 sc[1] = dot( P2, P2 );
334 sa[2] = dot( pD, pD );
335 sb[2] = sa[1];
336 sc[2] = dot( P1, P1 );
337 B[0] = barrier( 1, sa[0], sb[0], sc[0], 3 );
338 B[2] = barrier( 1, sa[2], sb[2], sc[2], rD );
339 calt1( P3, P4, t1V );
340 calt1( pA, P1, t1D );
341 if ( L == 0 )
342 {
343 for ( int i = 0; i != 4; i++ )
344 {
345 for ( int j = 0; j != 4; j++ )
346 {
347 temp_PDF +=
348 t1D[i] * ( G[i][j] - pA[i] * pA[j] / sa[1] ) * t1V[j] * ( G[i][i] ) * ( G[j][j] );
349 }
350 }
351 B[1] = 1;
352 }
353 if ( L == 2 )
354 {
355 calt2( pV, P2, t2A );
356 for ( int i = 0; i != 4; i++ )
357 {
358 for ( int j = 0; j != 4; j++ )
359 { temp_PDF += t1D[i] * t2A[i][j] * t1V[j] * ( G[i][i] ) * ( G[j][j] ); }
360 }
361 B[1] = barrier( 2, sa[1], sb[1], sc[1], 3 );
362 }
363 amp_PDF = temp_PDF * B[0] * B[1] * B[2];
364 return amp_PDF;
365}
366EvtComplex EvtDToKSpipipi::D2AP_A2SP( double P1[], double P2[], double P3[], double P4[] ) {
367 double temp_PDF = 0;
368 EvtComplex amp_PDF( 0, 0 );
369 EvtComplex pro;
370 double sa[3], sb[3], sc[3], B[3];
371 double t1D[4], t1A[4];
372 double pS[4], pA[4], pD[4];
373 for ( int i = 0; i != 4; i++ )
374 {
375 pS[i] = P3[i] + P4[i];
376 pA[i] = pS[i] + P2[i];
377 pD[i] = pA[i] + P1[i];
378 }
379 sa[0] = dot( pS, pS );
380 sb[0] = dot( P3, P3 );
381 sc[0] = dot( P4, P4 );
382 sa[1] = dot( pA, pA );
383 sb[1] = sa[0];
384 sc[1] = dot( P2, P2 );
385 sa[2] = dot( pD, pD );
386 sb[2] = sa[1];
387 sc[2] = dot( P1, P1 );
388 B[1] = barrier( 1, sa[1], sb[1], sc[1], 3 );
389 B[2] = barrier( 1, sa[2], sb[2], sc[2], rD );
390 calt1( pA, P1, t1D );
391 calt1( pS, P2, t1A );
392 for ( int i = 0; i != 4; i++ ) { temp_PDF += t1D[i] * t1A[i] * ( G[i][i] ); }
393 amp_PDF = temp_PDF * B[1] * B[2];
394 return amp_PDF;
395}
396EvtComplex EvtDToKSpipipi::D2PP_P2VP( double P1[], double P2[], double P3[], double P4[] ) {
397 double temp_PDF = 0;
398 EvtComplex amp_PDF( 0, 0 );
399 EvtComplex pro;
400 double sa[3], sb[3], sc[3], B[3];
401 double t1V[4];
402 double pV[4], pP[4], pD[4];
403 for ( int i = 0; i != 4; i++ )
404 {
405 pV[i] = P3[i] + P4[i];
406 pP[i] = pV[i] + P2[i];
407 pD[i] = pP[i] + P1[i];
408 }
409 sa[0] = dot( pV, pV );
410 sb[0] = dot( P3, P3 );
411 sc[0] = dot( P4, P4 );
412 sa[1] = dot( pP, pP );
413 sb[1] = sa[0];
414 sc[1] = dot( P2, P2 );
415 sa[2] = dot( pD, pD );
416 sb[2] = sa[1];
417 sc[2] = dot( P1, P1 );
418 B[0] = barrier( 1, sa[0], sb[0], sc[0], 3 );
419 B[1] = barrier( 1, sa[1], sb[1], sc[1], 3 );
420 calt1( P3, P4, t1V );
421 for ( int i = 0; i != 4; i++ ) { temp_PDF += P2[i] * t1V[i] * ( G[i][i] ); }
422 amp_PDF = temp_PDF * B[0] * B[1];
423 return amp_PDF;
424}
425EvtComplex EvtDToKSpipipi::D2VP_V2VP( double P1[], double P2[], double P3[], double P4[] ) {
426 double temp_PDF = 0;
427 EvtComplex amp_PDF( 0, 0 );
428 EvtComplex pro;
429 double sa[3], sb[3], sc[3], B[3];
430 double pV1[4], pV2[4], qV1[4], qV2[4], pD[4];
431 for ( int i = 0; i != 4; i++ )
432 {
433 pV2[i] = P3[i] + P4[i];
434 qV2[i] = P3[i] - P4[i];
435 pV1[i] = pV2[i] + P2[i];
436 qV1[i] = pV2[i] - P2[i];
437 pD[i] = pV1[i] + P1[i];
438 }
439 for ( int i = 0; i != 4; i++ )
440 {
441 for ( int j = 0; j != 4; j++ )
442 {
443 for ( int k = 0; k != 4; k++ )
444 {
445 for ( int l = 0; l != 4; l++ )
446 {
447 temp_PDF += E[i][j][k][l] * pV1[i] * qV1[j] * P1[k] * qV2[l] * ( G[i][i] ) *
448 ( G[j][j] ) * ( G[k][k] ) * ( G[l][l] );
449 }
450 }
451 }
452 }
453 sa[0] = dot( pV2, pV2 );
454 sb[0] = dot( P3, P3 );
455 sc[0] = dot( P4, P4 );
456 sa[1] = dot( pV1, pV1 );
457 sb[1] = sa[0];
458 sc[1] = dot( P2, P2 );
459 sa[2] = dot( pD, pD );
460 sb[2] = sa[1];
461 sc[2] = dot( P1, P1 );
462 B[0] = barrier( 1, sa[0], sb[0], sc[0], 3.0 );
463 B[1] = barrier( 1, sa[1], sb[1], sc[1], 3.0 );
464 B[2] = barrier( 1, sa[2], sb[2], sc[2], rD );
465 amp_PDF = temp_PDF * B[0] * B[1] * B[2];
466 return amp_PDF;
467}
468EvtComplex EvtDToKSpipipi::PHSP( double Km[], double Pip[] ) {
469 EvtComplex amp_PDF( 0, 0 );
470 double sa, sb, sc;
471 double KPi[4];
472 for ( int i = 0; i != 4; i++ ) { KPi[i] = Km[i] + Pip[i]; }
473 sa = dot( KPi, KPi );
474 sb = dot( Km, Km );
475 sc = dot( Pip, Pip );
476 amp_PDF = KPiSFormfactor( sa, sb, sc, 3.0 );
477 return amp_PDF;
478}
479EvtComplex EvtDToKSpipipi::getprop( double daug1[], double daug2[], double mass, double width,
480 int flag, int L ) {
481 // flag = 0, RBW with constant width
482 // flag = 1, RBW with width depends on momentum
483 // flag = 2, GS convolute with RBW_omega
484 // flag = 3, flatte for f_0(980)/a_0(980)
485 // flag = 4, Bugg formula for sigma;
486 // effect radii == 3.0 GeV^-1
487 EvtComplex prop1, prop2, prop;
488 EvtComplex one( 1, 0 );
489 double pR[4];
490 for ( int i = 0; i < 4; i++ ) { pR[i] = daug1[i] + daug2[i]; }
491 double sa, sb, sc;
492 sa = dot( pR, pR );
493 sb = dot( daug1, daug1 );
494 sc = dot( daug2, daug2 );
495 if ( flag == 0 ) return propogator( mass, width, sa );
496 if ( flag == 1 ) return propagatorRBW( mass, width, sa, sb, sc, 3.0, L );
497 if ( flag == 2 )
498 {
499 prop1 = propagatorGS( mass, width, sa, sb, sc, 3.0, L );
500 prop2 = propagatorRBW( 0.783, 0.008, sa, sb, sc, 3.0, L );
501 EvtComplex coef_omega( rho_omega * cos( phi_omega ), rho_omega * sin( phi_omega ) );
502 prop = prop1 * ( one + 0.783 * 0.783 * coef_omega * prop2 );
503 return prop;
504 }
505 if ( flag == 3 )
506 {
507 // Not need for D+ -> Ks 3pi
508 }
509 if ( flag == 4 )
510 {
511 EvtComplex ci( 0, 1 );
512 double f = 0.5843 + 1.6663 * sa;
513 double M = 0.9264;
514 double mpi2 = mass_Pion * mass_Pion;
515 double mass2 = M * M;
516 double g1 = f * ( sa - mpi2 / 2 ) / ( mass2 - mpi2 / 2 ) * exp( ( mass2 - sa ) / 1.082 );
517 EvtComplex rho1s = rhoab( sa, sb, sc );
518 EvtComplex rho1M = rhoab( mass2, sb, sc );
519 EvtComplex rho2s = rho4Pi( sa );
520 EvtComplex rho2M = rho4Pi( mass2 );
521 prop = 1.0 / ( M * M - sa - ci * M * ( g1 * rho1s / rho1M + 0.0024 * rho2s / rho2M ) );
522 return prop;
523 }
524 cout << "Please set a right flag! " << flag << endl;
525 return one;
526}
527EvtComplex EvtDToKSpipipi::rhoab( double sa, double sb, double sc ) {
528 EvtComplex one( 1, 0 );
529 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
530 EvtComplex rho;
531 EvtComplex ci( 0, 1 );
532 if ( q > 0 ) rho = one * sqrt( q / sa );
533 if ( q < 0 ) rho = ci * sqrt( -q / sa );
534 rho = 2.0 * rho;
535 return rho;
536}
537EvtComplex EvtDToKSpipipi::rho4Pi( double sa ) {
538 double mpi = 0.13957;
539 EvtComplex one( 1, 0 );
540 EvtComplex rho( 0, 0 );
541 EvtComplex ci( 0, 1 );
542 double temp = 1 - 16 * mpi * mpi / sa;
543 if ( temp > 0 ) rho = one * sqrt( temp ) / ( 1 + exp( 9.8 - 3.5 * sa ) );
544 if ( temp < 0 ) rho = ci * sqrt( -temp ) / ( 1 + exp( 9.8 - 3.5 * sa ) );
545 return rho;
546}
547
548EvtComplex EvtDToKSpipipi::propogator( double mass, double width, double sx ) const {
549 EvtComplex ci( 0, 1 );
550 EvtComplex prop = 1.0 / ( mass * mass - sx - ci * mass * width );
551 return prop;
552}
553double EvtDToKSpipipi::wid( double mass, double sa, double sb, double sc, double r,
554 int l ) const {
555 double widm( 0. ), q( 0. ), q0( 0. );
556 double sa0 = mass * mass;
557 double m = sqrt( sa );
558 q = Qabcs( sa, sb, sc );
559 q0 = Qabcs( sa0, sb, sc );
560 double z, z0;
561 z = q * r * r;
562 z0 = q0 * r * r;
563 double t = q / q0;
564 double F( 0. );
565 if ( l == 0 ) F = 1;
566 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
567 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
568 widm = pow( t, l + 0.5 ) * mass / m * F * F;
569 return widm;
570}
571double EvtDToKSpipipi::h( double m, double q ) const {
572 double h( 0. );
573 h = 2 / pi * q / m * log( ( m + 2 * q ) / ( 2 * mpi ) );
574 return h;
575}
576double EvtDToKSpipipi::dh( double mass, double q0 ) const {
577 double dh = h( mass, q0 ) * ( 1.0 / ( 8 * q0 * q0 ) - 1.0 / ( 2 * mass * mass ) ) +
578 1.0 / ( 2 * pi * mass * mass );
579 return dh;
580}
581double EvtDToKSpipipi::f( double mass, double sx, double q0, double q ) const {
582 double m = sqrt( sx );
583 double f = mass * mass / ( pow( q0, 3 ) ) *
584 ( q * q * ( h( m, q ) - h( mass, q0 ) ) +
585 ( mass * mass - sx ) * q0 * q0 * dh( mass, q0 ) );
586 return f;
587}
588double EvtDToKSpipipi::d( double mass, double q0 ) const {
589 double d = 3.0 / pi * mpi * mpi / ( q0 * q0 ) * log( ( mass + 2 * q0 ) / ( 2 * mpi ) ) +
590 mass / ( 2 * pi * q0 ) - ( mpi * mpi * mass ) / ( pi * pow( q0, 3 ) );
591 return d;
592}
593EvtComplex EvtDToKSpipipi::propagatorRBW( double mass, double width, double sa, double sb,
594 double sc, double r, int l ) const {
595 EvtComplex ci( 0, 1 );
596 EvtComplex prop =
597 1.0 / ( mass * mass - sa - ci * mass * width * wid( mass, sa, sb, sc, r, l ) );
598 return prop;
599}
600EvtComplex EvtDToKSpipipi::propagatorGS( double mass, double width, double sa, double sb,
601 double sc, double r, int l ) const {
602 EvtComplex ci( 0, 1 );
603 double q = Qabcs( sa, sb, sc );
604 double sa0 = mass * mass;
605 double q0 = Qabcs( sa0, sb, sc );
606 q = sqrt( q );
607 q0 = sqrt( q0 );
608 EvtComplex prop = ( 1 + d( mass, q0 ) * width / mass ) /
609 ( mass * mass - sa + width * f( mass, sa, q0, q ) -
610 ci * mass * width * wid( mass, sa, sb, sc, r, l ) );
611 return prop;
612}
613double EvtDToKSpipipi::Flatte_rhoab( double sa, double sb, double sc ) const {
614 double q = Qabcs( sa, sb, sc );
615 double rho = sqrt( q / sa );
616 return rho;
617}
618EvtComplex EvtDToKSpipipi::propagatorFlatte( double mass, double width, double sx, double* sb,
619 double* sc ) const {
620 EvtComplex ci( 0, 1 );
621 double rho1 = Flatte_rhoab( sx, sb[0], sc[0] );
622 double rho2 = Flatte_rhoab( sx, sb[1], sc[1] );
623 EvtComplex prop = 1.0 / ( mass * mass - sx - ci * ( g1 * g1 * rho1 + g2 * g2 * rho2 ) );
624 return prop;
625}
626
627double EvtDToKSpipipi::dot( double* a1, double* a2 ) const {
628 double dot = 0;
629 for ( int i = 0; i != 4; i++ ) { dot += a1[i] * a2[i] * G[i][i]; }
630 return dot;
631}
632double EvtDToKSpipipi::Qabcs( double sa, double sb, double sc ) const {
633 double Qabcs = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
634 if ( Qabcs < 0 ) Qabcs = 1e-16;
635 return Qabcs;
636}
637double EvtDToKSpipipi::barrier( double l, double sa, double sb, double sc, double r ) const {
638 double q = Qabcs( sa, sb, sc );
639 q = sqrt( q );
640 double z = q * r;
641 z = z * z;
642 double F = 1;
643 if ( l > 2 ) F = 0;
644 if ( l == 0 ) F = 1;
645 if ( l == 1 ) { F = sqrt( ( 2 * z ) / ( 1 + z ) ); }
646 if ( l == 2 ) { F = sqrt( ( 13 * z * z ) / ( 9 + 3 * z + z * z ) ); }
647 return F;
648}
649void EvtDToKSpipipi::calt1( double daug1[], double daug2[], double t1[] ) const {
650 double p, pq;
651 double pa[4], qa[4];
652 for ( int i = 0; i != 4; i++ )
653 {
654 pa[i] = daug1[i] + daug2[i];
655 qa[i] = daug1[i] - daug2[i];
656 }
657 p = dot( pa, pa );
658 pq = dot( pa, qa );
659 for ( int i = 0; i != 4; i++ ) { t1[i] = qa[i] - pq / p * pa[i]; }
660}
661void EvtDToKSpipipi::calt2( double daug1[], double daug2[], double t2[][4] ) const {
662 double p, r;
663 double pa[4], t1[4];
664 calt1( daug1, daug2, t1 );
665 r = dot( t1, t1 );
666 for ( int i = 0; i != 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
667 p = dot( pa, pa );
668 for ( int i = 0; i != 4; i++ )
669 {
670 for ( int j = 0; j != 4; j++ )
671 { t2[i][j] = t1[i] * t1[j] - 1.0 / 3 * r * ( G[i][j] - pa[i] * pa[j] / p ); }
672 }
673}
mass
double mass
Definition CosmicGenerator.cxx:128
pi2
character *LEPTONflag integer iresonances real pi2
Definition Ekhara/src/ekhara/common.h:3
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:209
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
one
const DifNumber one
q
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition KKsem.h:33
EvtDToKSpipipi::getName
void getName(std::string &name)
Definition EvtDToKSpipipi.cc:33
EvtDToKSpipipi::decay
void decay(EvtParticle *p)
Definition EvtDToKSpipipi.cc:158
EvtDToKSpipipi::clone
EvtDecayBase * clone()
Definition EvtDToKSpipipi.cc:35
EvtDToKSpipipi::~EvtDToKSpipipi
virtual ~EvtDToKSpipipi()
Definition EvtDToKSpipipi.cc:31
EvtDecayBase::checkSpinParent
void checkSpinParent(EvtSpinType::spintype sp)
Definition EvtDecayBase.cc:479
EvtDecayBase::setProbMax
void setProbMax(double prbmx)
Definition EvtDecayBase.cc:261
EvtDecayBase::checkNDaug
void checkNDaug(int d1, int d2=-1)
Definition EvtDecayBase.cc:465
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition EvtDecayBase.hh:68
EvtDecayBase::checkNArg
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
Definition EvtDecayBase.cc:449
EvtDecayProb::setProb
void setProb(double prob)
Definition EvtDecayProb.hh:40
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition EvtParticle.cc:108
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition EvtParticle.cc:80
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition EvtParticle.cc:1042
EvtVector4R::get
double get(int i) const
Definition EvtVector4R.hh:171
rb::mpi2
const double mpi2
Definition TConstant.h:28
t
int t()
Definition t.c:1