BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtD0ToKpipipi.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: EvtD0ToKpipipi.cc
11// the necessary file: EvtD0ToKpipipi.hh
12//
13// Description: D0 -> K- pi+ pi+ pi-,
14// see PHYSICAL REVIEW D 95, 072010 (2017)
15//
16// Modification history:
17//
18// Liaoyuan Dong Jan.15, 2020 Module created
19//
20//------------------------------------------------------------------------
21#include "EvtD0ToKpipipi.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
31EvtD0ToKpipipi::~EvtD0ToKpipipi() {}
32
33void EvtD0ToKpipipi::getName( std::string& model_name ) { model_name = "D0ToKpipipi"; }
34
35EvtDecayBase* EvtD0ToKpipipi::clone() { return new EvtD0ToKpipipi; }
36
37void EvtD0ToKpipipi::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 EvtD0ToKpipipi" << std::endl;
48
49 width[0] = 0.09;
50 width[1] = 0.044183653178315;
51 width[2] = 0.541879469380012;
52 width[3] = 0.148423336450619;
53 mass[0] = 1.272;
54 mass[1] = 0.894781734682169;
55 mass[2] = 1.3622013558915;
56 mass[3] = 0.779143408171384;
57
58 phi[0] = 2.34794687054858;
59 rho[0] = 0.0759345115620669;
60 phi[1] = -2.24641399153466;
61 rho[1] = 0.0383327604903577;
62 phi[2] = 2.48955684856045;
63 rho[2] = 0.0931445480476023;
64
65 phi[3] = 0;
66 rho[3] = 1;
67
68 phi[4] = -2.10558220063012;
69 rho[4] = 0.347041869435286;
70 phi[5] = 1.47445088061872;
71 phi[6] = 3.00243265559304;
72 rho[5] = 0.00965088341753795;
73 rho[6] = 0.120536507325731;
74 phi[7] = -2.45477499325158;
75 rho[7] = 0.101419048440676;
76 phi[8] = -1.35809992343491;
77 rho[8] = 4.28149643321317;
78 phi[9] = -2.45149221243198;
79 rho[9] = 0.339492272598394;
80 phi[10] = -0.17419389225461;
81 rho[10] = -0.143619437541254;
82
83 phi[11] = -2.08744386934208;
84 rho[11] = 0.296286583716349;
85 phi[12] = 0.;
86 rho[12] = 0.;
87 phi[13] = -0.432190571560873;
88 rho[13] = 0.657344690733276;
89 phi[14] = -1.39790294886865;
90 rho[14] = 1.71208007006123;
91 phi[15] = 1.58945300476228;
92 rho[15] = 3.58248347683687;
93 phi[16] = 2.58249107256307;
94 rho[16] = -1.10728829503506;
95 phi[17] = -0.163623135170955;
96 rho[17] = 1.70863070178363;
97 phi[18] = -0.134699023080211;
98 rho[18] = 0.567531283682344;
99 phi[19] = -2.12670610368279;
100 rho[19] = 0.276571752504914;
101 phi[20] = -1.3352622107357;
102 rho[20] = 0.416634203151278;
103
104 phi[21] = -2.91571684221842;
105 rho[21] = 0.423062298489176;
106 phi[22] = 2.4544220004327;
107 rho[22] = 1.4017194038459;
108 phi[23] = -2.23388390670423;
109 rho[23] = 4.11110400629068;
110
111 // for (int i=0; i<24; i++) {
112 // cout << i << " rho,phi = " << rho[i] << ", "<< phi[i] << endl;
113 // }
114
115 mD = 1.86486;
116 rRes = 3.0;
117 rD = 5.0;
118 metap = 0.95778;
119 mkstr = 0.89594;
120 mk0 = 0.497614;
121 mass_Kaon = 0.49368;
122 mass_Pion = 0.13957;
123 mass_Pi0 = 0.1349766;
124 math_pi = 3.1415926;
125
126 pi = 3.1415926;
127 mpi = 0.13957;
128 g1 = 0.5468;
129 g2 = 0.23;
130
131 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
132 int EE[4][4][4][4] = {
133 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
134 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, -1, 0 } },
135 { { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 1, 0, 0 } },
136 { { 0, 0, 0, 0 }, { 0, 0, 1, 0 }, { 0, -1, 0, 0 }, { 0, 0, 0, 0 } } },
137 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 1, 0 } },
138 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
139 { { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 } },
140 { { 0, 0, -1, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 } } },
141 { { { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, -1, 0, 0 } },
142 { { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 } },
143 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
144 { { 0, 1, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } },
145 { { { 0, 0, 0, 0 }, { 0, 0, -1, 0 }, { 0, 1, 0, 0 }, { 0, 0, 0, 0 } },
146 { { 0, 0, 1, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 } },
147 { { 0, -1, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
148 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } } };
149 for ( int i = 0; i < 4; i++ )
150 {
151 for ( int j = 0; j < 4; j++ )
152 {
153 G[i][j] = GG[i][j];
154 for ( int k = 0; k < 4; k++ )
155 {
156 for ( int l = 0; l < 4; l++ ) { E[i][j][k][l] = EE[i][j][k][l]; }
157 }
158 }
159 }
160}
161
162void EvtD0ToKpipipi::initProbMax() { setProbMax( 720.0 ); }
163
164void EvtD0ToKpipipi::decay( EvtParticle* p ) {
165 /*
166 double maxprob = 0.0;
167 for(int ir=0;ir<=60000000;ir++){
168 p->initializePhaseSpace(getNDaug(),getDaugs());
169 EvtVector4R Km0 = p->getDaug(0)->getP4();
170 EvtVector4R pi1 = p->getDaug(1)->getP4();
171 EvtVector4R pi2 = p->getDaug(2)->getP4();
172 EvtVector4R pi3 = p->getDaug(3)->getP4();
173 double Km[4],Pip1[4],Pip2[4],Pim[4];
174 Km[0] = Km0.get(0); Pip1[0] = pi1.get(0); Pip2[0] = pi2.get(0); Pim[0] = pi3.get(0);
175 Km[1] = Km0.get(1); Pip1[1] = pi1.get(1); Pip2[1] = pi2.get(1); Pim[1] = pi3.get(1);
176 Km[2] = Km0.get(2); Pip1[2] = pi1.get(2); Pip2[2] = pi2.get(2); Pim[2] = pi3.get(2);
177 Km[3] = Km0.get(3); Pip1[3] = pi1.get(3); Pip2[3] = pi2.get(3); Pim[3] = pi3.get(3);
178 double Prob = calPDF(Km, Pip1, Pip2, Pim);
179 if(Prob>maxprob) {
180 maxprob=Prob;
181 std::cout << "Max PDF = " << ir << " prob= " << Prob << std::endl;
182 }
183 }
184 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
185 */
186 p->initializePhaseSpace( getNDaug(), getDaugs() );
187 EvtVector4R Km0 = p->getDaug( 0 )->getP4();
188 EvtVector4R pi1 = p->getDaug( 1 )->getP4();
189 EvtVector4R pi2 = p->getDaug( 2 )->getP4();
190 EvtVector4R pi3 = p->getDaug( 3 )->getP4();
191
192 double Km[4], Pip1[4], Pip2[4], Pim[4];
193 Km[0] = Km0.get( 0 );
194 Pip1[0] = pi1.get( 0 );
195 Pip2[0] = pi2.get( 0 );
196 Pim[0] = pi3.get( 0 );
197 Km[1] = Km0.get( 1 );
198 Pip1[1] = pi1.get( 1 );
199 Pip2[1] = pi2.get( 1 );
200 Pim[1] = pi3.get( 1 );
201 Km[2] = Km0.get( 2 );
202 Pip1[2] = pi1.get( 2 );
203 Pip2[2] = pi2.get( 2 );
204 Pim[2] = pi3.get( 2 );
205 Km[3] = Km0.get( 3 );
206 Pip1[3] = pi1.get( 3 );
207 Pip2[3] = pi2.get( 3 );
208 Pim[3] = pi3.get( 3 );
209 double prob = calPDF( Km, Pip1, Pip2, Pim );
210 setProb( prob );
211 return;
212}
213
214double EvtD0ToKpipipi::calPDF( double Km[], double Pip1[], double Pip2[], double Pim[] ) {
215 Km[0] = sqrt( mass_Kaon * mass_Kaon + Km[1] * Km[1] + Km[2] * Km[2] + Km[3] * Km[3] );
216 Pip1[0] = sqrt( mass_Pion * mass_Pion + Pip1[1] * Pip1[1] + Pip1[2] * Pip1[2] +
217 Pip1[3] * Pip1[3] );
218 Pip2[0] = sqrt( mass_Pion * mass_Pion + Pip2[1] * Pip2[1] + Pip2[2] * Pip2[2] +
219 Pip2[3] * Pip2[3] );
220 Pim[0] = sqrt( mass_Pion * mass_Pion + Pim[1] * Pim[1] + Pim[2] * Pim[2] + Pim[3] * Pim[3] );
221
222 EvtComplex PDF[24];
223 int g[3];
224 g[0] = 1;
225 g[1] = 1;
226
227 //-----------D->K*rho--------
228 g[2] = 0;
229 PDF[0] = D2VV( Km, Pip1, Pip2, Pim, g ) + D2VV( Km, Pip2, Pip1, Pim, g );
230 g[2] = 1;
231 PDF[1] = D2VV( Km, Pip1, Pip2, Pim, g ) + D2VV( Km, Pip2, Pip1, Pim, g );
232 g[2] = 2;
233 PDF[2] = D2VV( Km, Pip1, Pip2, Pim, g ) + D2VV( Km, Pip2, Pip1, Pim, g );
234 //-----------D->K*rho finish--
235 //----------D->a1K------------
236 g[2] = 0;
237 PDF[3] = D2AP_A2VP( Km, Pip2, Pip1, Pim, g, 2 ) + D2AP_A2VP( Km, Pip1, Pip2, Pim, g, 2 );
238 g[2] = 2;
239 PDF[4] = D2AP_A2VP( Km, Pip2, Pip1, Pim, g, 2 ) + D2AP_A2VP( Km, Pip1, Pip2, Pim, g, 2 );
240 //----------D->a1K finish-----
241 //--D->K1Pi--K1->K*Pi---------
242 g[2] = 0;
243 PDF[5] = D2AP_A2VP( Pip2, Pim, Km, Pip1, g, 1 ) + D2AP_A2VP( Pip1, Pim, Km, Pip2, g, 1 );
244 g[2] = 2;
245 PDF[6] = D2AP_A2VP( Pip2, Pim, Km, Pip1, g, 1 ) + D2AP_A2VP( Pip1, Pim, Km, Pip2, g, 1 );
246 //--D->K1Pi--K1->K*Pi-finish--
247 //--D->K1Pi--K1->rhoK---------
248 g[2] = 0;
249 PDF[7] = D2AP_A2VP( Pip2, Km, Pip1, Pim, g, 3 ) + D2AP_A2VP( Pip1, Km, Pip2, Pim, g, 3 );
250 //--D->K1Pi--K1->rhoK-finish--
251 //--D->K1Pi--K1->K*0(1430)Pi--
252 PDF[8] = D2AP_A2SP( Pip2, Pim, Km, Pip1, 1 ) + D2AP_A2SP( Pip1, Pim, Km, Pip2, 1 );
253 //--------res finish----------
254 //----------------non-res------------------
255 //--------KPiSrho------------------
256 PDF[9] = D2VS( Pip1, Pim, Km, Pip2, 1, 2 ) + D2VS( Pip2, Pim, Km, Pip1, 1, 2 );
257 //--------K*PiPiS-----------------
258 PDF[10] = D2VS( Km, Pip1, Pip2, Pim, 1, 1 ) + D2VS( Km, Pip2, Pip1, Pim, 1, 1 );
259 //--------K*PiP-------------------
260 PDF[11] = D2PP_P2VP( Pip2, Pim, Km, Pip1, 1 ) + D2PP_P2VP( Pip1, Pim, Km, Pip2, 1 );
261 //--------rhoKA--------------------
262 g[0] = 1;
263 g[1] = 0;
264 g[2] = 0;
265 PDF[12] = 0;
266 g[2] = 2;
267 PDF[13] = D2AP_A2VP( Pip2, Km, Pip1, Pim, g, 3 ) + D2AP_A2VP( Pip1, Km, Pip2, Pim, g, 3 );
268 //-------PHSP-----------------------
269 PDF[14] = PHSP( Km, Pip1 ) + PHSP( Km, Pip2 );
270 //------KPiVPiPiV-----------------------
271 g[0] = 0;
272 g[1] = 0;
273 g[2] = 0;
274 PDF[15] = D2VV( Km, Pip1, Pip2, Pim, g ) + D2VV( Km, Pip2, Pip1, Pim, g );
275 //------KPiVPiPiS----------------------
276 PDF[16] = D2VS( Km, Pip1, Pip2, Pim, 0, 1 ) + D2VS( Km, Pip2, Pip1, Pim, 0, 1 );
277 //------KPiSPiPiV----------------------
278 PDF[17] = D2VS( Pip1, Pim, Km, Pip2, 0, 2 ) + D2VS( Pip2, Pim, Km, Pip1, 0, 2 );
279 //-------------------------------------
280 PDF[18] = D2PP_P2VP( Pip2, Km, Pip1, Pim, 2 ) + D2PP_P2VP( Pip1, Km, Pip2, Pim, 2 );
281 //-------------------------------------
282 PDF[19] = D2VP_V2VP( Pip2, Pim, Km, Pip1, 1 ) + D2VP_V2VP( Pip1, Pim, Km, Pip2, 1 );
283 //-------------------------------------
284 PDF[20] = D2VP_V2VP( Pip2, Km, Pip1, Pim, 2 ) + D2VP_V2VP( Pip1, Km, Pip2, Pim, 2 );
285 //-------------------------------------
286 PDF[21] = D2TS( Km, Pip1, Pip2, Pim, 1 ) + D2TS( Km, Pip2, Pip1, Pim, 1 );
287 PDF[22] = D2TS( Pip1, Pim, Km, Pip2, 2 ) + D2TS( Pip2, Pim, Km, Pip1, 2 );
288 PDF[23] = D2AP_A2SP( Km, Pip2, Pip1, Pim, 2 ) + D2AP_A2SP( Km, Pip1, Pip2, Pim, 2 );
289 //------------------------------------------
290 EvtComplex cof;
291 EvtComplex pdf, module;
292 pdf = EvtComplex( 0, 0 );
293 for ( int i = 0; i != 24; i++ )
294 {
295 cof = EvtComplex( rho[i] * cos( phi[i] ), rho[i] * sin( phi[i] ) );
296 // cout<<i << " " << (cof*PDF[i])<<" : "<<cof<<" "<<PDF[i]<<endl;
297 pdf = pdf + cof * PDF[i];
298 }
299 module = conj( pdf ) * pdf;
300 double value;
301 value = real( module );
302 return ( value <= 0 ) ? 1e-20 : value;
303}
304
305EvtComplex EvtD0ToKpipipi::KPiSFormfactor( double sa, double sb, double sc, double r ) {
306 double m1430 = 1.463;
307 double sa0 = m1430 * m1430;
308 double w1430 = 0.233;
309 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
310 if ( q0 < 0 ) q0 = 1e-16;
311 double qs = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
312 double q = sqrt( qs );
313 double width = w1430 * q * m1430 / sqrt( sa * q0 );
314 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
315 if ( temp_R < 0 ) temp_R += math_pi;
316 double deltaR = -5.31 + temp_R;
317 double temp_F = atan( 2 * 1.07 * q / ( 2 + ( -1.8 ) * 1.07 * qs ) );
318 if ( temp_F < 0 ) temp_F += math_pi;
319 double deltaF = 2.33 + temp_F;
320 EvtComplex cR( cos( deltaR ), sin( deltaR ) );
321 EvtComplex cF( cos( deltaF ), sin( deltaF ) );
322 EvtComplex amp = 0.8 * sin( deltaF ) * cF + sin( deltaR ) * cR * cF * cF;
323 return amp;
324}
325
326EvtComplex EvtD0ToKpipipi::D2VV( double P1[], double P2[], double P3[], double P4[],
327 int g[] ) {
328 double t1V1[4], t1V2[4], t1D[4], t2D[4][4];
329 double temp_PDF = 0;
330 EvtComplex amp_PDF( 0, 0 );
331 EvtComplex pro[2];
332 //---------use g[0],g[1] to define res or non-res, g[2] represent S, P or D wave
333 double sa[3], sb[3], sc[3], B[3];
334 double pV1[4], pV2[4], pD[4];
335 for ( int i = 0; i != 4; i++ )
336 {
337 pV1[i] = P1[i] + P2[i];
338 pV2[i] = P3[i] + P4[i];
339 pD[i] = pV1[i] + pV2[i];
340 }
341 sa[0] = dot( pV1, pV1 );
342 sb[0] = dot( P1, P1 );
343 sc[0] = dot( P2, P2 );
344 sa[1] = dot( pV2, pV2 );
345 sb[1] = dot( P3, P3 );
346 sc[1] = dot( P4, P4 );
347 sa[2] = dot( pD, pD );
348 sb[2] = sa[0];
349 sc[2] = sa[1];
350 if ( g[1] == 1 )
351 { pro[1] = propagatorGS( mass[3], width[3], sa[1], sb[1], sc[1], rRes, 1 ); }
352 if ( g[0] == 1 )
353 { pro[0] = propagatorRBW( mass[1], width[1], sa[0], sb[0], sc[0], rRes, 1 ); }
354 if ( g[0] == 0 ) pro[0] = 1;
355 if ( g[1] == 0 ) pro[1] = 1;
356 B[0] = barrier( 1, sa[0], sb[0], sc[0], rRes );
357 B[1] = barrier( 1, sa[1], sb[1], sc[1], rRes );
358 calt1( P1, P2, t1V1 );
359 calt1( P3, P4, t1V2 );
360 if ( g[2] == 0 )
361 {
362 for ( int i = 0; i != 4; i++ ) { temp_PDF += ( G[i][i] ) * t1V1[i] * t1V2[i]; }
363 B[2] = 1;
364 }
365 if ( g[2] == 1 )
366 {
367 calt1( pV1, pV2, t1D );
368 for ( int i = 0; i != 4; i++ )
369 {
370 for ( int j = 0; j != 4; j++ )
371 {
372 for ( int k = 0; k != 4; k++ )
373 {
374 for ( int l = 0; l != 4; l++ )
375 {
376 temp_PDF += E[i][j][k][l] * pD[i] * t1D[j] * t1V1[k] * t1V2[l] * ( G[i][i] ) *
377 ( G[j][j] ) * ( G[l][l] ) * ( G[k][k] );
378 }
379 }
380 }
381 }
382 B[2] = barrier( 1, sa[2], sb[2], sc[2], rD );
383 }
384 if ( g[2] == 2 )
385 {
386 calt2( pV1, pV2, t2D );
387 for ( int i = 0; i != 4; i++ )
388 {
389 for ( int j = 0; j != 4; j++ )
390 { temp_PDF += t2D[i][j] * t1V1[i] * t1V2[j] * ( G[i][i] ) * ( G[j][j] ); }
391 }
392 B[2] = barrier( 2, sa[2], sb[2], sc[2], rD );
393 }
394 amp_PDF = temp_PDF * B[0] * B[1] * B[2] * pro[0] * pro[1];
395 return amp_PDF;
396}
397
398EvtComplex EvtD0ToKpipipi::D2AP_A2VP( double P1[], double P2[], double P3[], double P4[],
399 int g[], int flag ) {
400 // flag = 1, V = K*, A = K1(1270); flag = 2, V = rho, A = a1(1260);
401 // flag = 3, V = rho, A = K1(1270);
402 double temp_PDF = 0;
403 EvtComplex amp_PDF( 0, 0 );
404 EvtComplex pro[2];
405 double t1V[4], t1D[4], t2A[4][4];
406 double sa[3], sb[3], sc[3], B[3];
407 double pV[4], pA[4], pD[4];
408 for ( int i = 0; i != 4; i++ )
409 {
410 pV[i] = P3[i] + P4[i];
411 pA[i] = pV[i] + P2[i];
412 pD[i] = pA[i] + P1[i];
413 }
414 sa[0] = dot( pV, pV );
415 sb[0] = dot( P3, P3 );
416 sc[0] = dot( P4, P4 );
417 sa[1] = dot( pA, pA );
418 sb[1] = sa[0];
419 sc[1] = dot( P2, P2 );
420 sa[2] = dot( pD, pD );
421 sb[2] = sa[1];
422 sc[2] = dot( P1, P1 );
423 if ( g[0] == 1 )
424 {
425 if ( flag == 1 ) pro[0] = propagatorRBW( mass[1], width[1], sa[0], sb[0], sc[0], rRes, 1 );
426 if ( flag == 2 || flag == 3 )
427 pro[0] = propagatorGS( mass[3], width[3], sa[0], sb[0], sc[0], rRes, 1 );
428 }
429 if ( g[1] == 1 )
430 {
431 if ( flag == 1 ) pro[1] = propogator( mass[0], width[0], sa[1] );
432 if ( flag == 2 )
433 pro[1] = propagatorRBW( mass[2], width[2], sa[1], sb[1], sc[1], rRes, g[2] );
434 if ( flag == 3 ) pro[1] = propogator( mass[0], width[0], sa[1] );
435 }
436 if ( g[0] == 0 ) pro[0] = 1;
437 if ( g[1] == 0 ) pro[1] = 1;
438 B[0] = barrier( 1, sa[0], sb[0], sc[0], rRes );
439 B[2] = barrier( 1, sa[2], sb[2], sc[2], rD );
440 calt1( P3, P4, t1V );
441 calt1( pA, P1, t1D );
442 if ( g[2] == 0 )
443 {
444 for ( int i = 0; i != 4; i++ )
445 {
446 for ( int j = 0; j != 4; j++ )
447 {
448 temp_PDF +=
449 t1D[i] * ( G[i][j] - pA[i] * pA[j] / sa[1] ) * t1V[j] * ( G[i][i] ) * ( G[j][j] );
450 }
451 }
452 B[1] = 1;
453 }
454 if ( g[2] == 2 )
455 {
456 calt2( pV, P2, t2A );
457 for ( int i = 0; i != 4; i++ )
458 {
459 for ( int j = 0; j != 4; j++ )
460 { temp_PDF += t1D[i] * t2A[i][j] * t1V[j] * ( G[i][i] ) * ( G[j][j] ); }
461 }
462 B[1] = barrier( 2, sa[1], sb[1], sc[1], rRes );
463 }
464 amp_PDF = temp_PDF * B[0] * B[1] * B[2] * pro[0] * pro[1];
465 return amp_PDF;
466}
467EvtComplex EvtD0ToKpipipi::D2AP_A2SP( double P1[], double P2[], double P3[], double P4[],
468 int flag ) {
469 // flag = 1, S = K*; flag = 2, S = rho
470 double temp_PDF = 0;
471 EvtComplex amp_PDF( 0, 0 );
472 EvtComplex pro;
473 double sa[3], sb[3], sc[3], B[3];
474 double t1D[4], t1A[4];
475 double pS[4], pA[4], pD[4];
476 for ( int i = 0; i != 4; i++ )
477 {
478 pS[i] = P3[i] + P4[i];
479 pA[i] = pS[i] + P2[i];
480 pD[i] = pA[i] + P1[i];
481 }
482 sa[0] = dot( pS, pS );
483 sb[0] = dot( P3, P3 );
484 sc[0] = dot( P4, P4 );
485 sa[1] = dot( pA, pA );
486 sb[1] = sa[0];
487 sc[1] = dot( P2, P2 );
488 sa[2] = dot( pD, pD );
489 sb[2] = sa[1];
490 sc[2] = dot( P1, P1 );
491 B[1] = barrier( 1, sa[1], sb[1], sc[1], rRes );
492 B[2] = barrier( 1, sa[2], sb[2], sc[2], rD );
493 calt1( pA, P1, t1D );
494 calt1( pS, P2, t1A );
495 for ( int i = 0; i != 4; i++ ) { temp_PDF += t1D[i] * t1A[i] * ( G[i][i] ); }
496 amp_PDF = temp_PDF * B[1] * B[2];
497 if ( flag == 1 ) amp_PDF *= KPiSFormfactor( sa[0], sb[0], sc[0], rRes );
498 return amp_PDF;
499}
500EvtComplex EvtD0ToKpipipi::D2PP_P2VP( double P1[], double P2[], double P3[], double P4[],
501 int flag ) {
502 // flag = 1, V = K*; flag = 2, V = rho
503 double temp_PDF = 0;
504 EvtComplex amp_PDF( 0, 0 );
505 EvtComplex pro;
506 double sa[3], sb[3], sc[3], B[3];
507 double t1V[4];
508 double pV[4], pP[4], pD[4];
509 for ( int i = 0; i != 4; i++ )
510 {
511 pV[i] = P3[i] + P4[i];
512 pP[i] = pV[i] + P2[i];
513 pD[i] = pP[i] + P1[i];
514 }
515 sa[0] = dot( pV, pV );
516 sb[0] = dot( P3, P3 );
517 sc[0] = dot( P4, P4 );
518 sa[1] = dot( pP, pP );
519 sb[1] = sa[0];
520 sc[1] = dot( P2, P2 );
521 sa[2] = dot( pD, pD );
522 sb[2] = sa[1];
523 sc[2] = dot( P1, P1 );
524 B[0] = barrier( 1, sa[0], sb[0], sc[0], rRes );
525 B[1] = barrier( 1, sa[1], sb[1], sc[1], rRes );
526 if ( flag == 1 ) pro = propagatorRBW( mass[1], width[1], sa[0], sb[0], sc[0], rRes, 1 );
527 if ( flag == 2 ) pro = propagatorGS( mass[3], width[3], sa[0], sb[0], sc[0], rRes, 1 );
528 calt1( P3, P4, t1V );
529 for ( int i = 0; i != 4; i++ ) { temp_PDF += P2[i] * t1V[i] * ( G[i][i] ); }
530 amp_PDF = temp_PDF * B[0] * B[1] * pro;
531 return amp_PDF;
532}
533EvtComplex EvtD0ToKpipipi::D2VP_V2VP( double P1[], double P2[], double P3[], double P4[],
534 int flag ) {
535 // flag = 1, (K*Pi)V; flag = 2, (rhoK)V
536 double temp_PDF = 0;
537 EvtComplex amp_PDF( 0, 0 );
538 EvtComplex pro;
539 double sa[3], sb[3], sc[3], B[3];
540 double pV1[4], pV2[4], qV1[4], qV2[4], pD[4];
541 for ( int i = 0; i != 4; i++ )
542 {
543 pV2[i] = P3[i] + P4[i];
544 qV2[i] = P3[i] - P4[i];
545 pV1[i] = pV2[i] + P2[i];
546 qV1[i] = pV2[i] - P2[i];
547 pD[i] = pV1[i] + P1[i];
548 }
549 for ( int i = 0; i != 4; i++ )
550 {
551 for ( int j = 0; j != 4; j++ )
552 {
553 for ( int k = 0; k != 4; k++ )
554 {
555 for ( int l = 0; l != 4; l++ )
556 {
557 temp_PDF += E[i][j][k][l] * pV1[i] * qV1[j] * P1[k] * qV2[l] * ( G[i][i] ) *
558 ( G[j][j] ) * ( G[k][k] ) * ( G[l][l] );
559 }
560 }
561 }
562 }
563 sa[0] = dot( pV2, pV2 );
564 sb[0] = dot( P3, P3 );
565 sc[0] = dot( P4, P4 );
566 sa[1] = dot( pV1, pV1 );
567 sb[1] = sa[0];
568 sc[1] = dot( P2, P2 );
569 sa[2] = dot( pD, pD );
570 sb[2] = sa[1];
571 sc[2] = dot( P1, P1 );
572 if ( flag == 1 ) pro = propagatorRBW( mass[1], width[1], sa[0], sb[0], sc[0], rRes, 1 );
573 if ( flag == 2 ) pro = propagatorGS( mass[3], width[3], sa[0], sb[0], sc[0], rRes, 1 );
574 B[0] = barrier( 1, sa[0], sb[0], sc[0], rRes );
575 B[1] = barrier( 1, sa[1], sb[1], sc[1], rRes );
576 B[2] = barrier( 1, sa[2], sb[2], sc[2], rD );
577 amp_PDF = temp_PDF * B[0] * B[1] * B[2] * pro;
578 return amp_PDF;
579}
580EvtComplex EvtD0ToKpipipi::D2VS( double P1[], double P2[], double P3[], double P4[], int g,
581 int flag ) {
582 // flag = 1, V = K*; flag = 2, V = rho
583 double temp_PDF = 0;
584 EvtComplex amp_PDF( 0, 0 );
585 EvtComplex pro;
586 double sa[3], sb[3], sc[3], B[3];
587 double t1D[4], t1V[4];
588 double pS[4], pV[4], pD[4];
589 for ( int i = 0; i != 4; i++ )
590 {
591 pS[i] = P3[i] + P4[i];
592 pV[i] = P1[i] + P2[i];
593 pD[i] = pS[i] + pV[i];
594 }
595 sa[0] = dot( pS, pS );
596 sb[0] = dot( P3, P3 );
597 sc[0] = dot( P4, P4 );
598 sa[1] = dot( pV, pV );
599 sb[1] = dot( P1, P1 );
600 sc[1] = dot( P2, P2 );
601 sa[2] = dot( pD, pD );
602 sb[2] = sa[0];
603 sc[2] = sa[1];
604 if ( g == 1 )
605 {
606 if ( flag == 2 ) pro = propagatorGS( mass[3], width[3], sa[1], sb[1], sc[1], rRes, 1 );
607 if ( flag == 1 ) pro = propagatorRBW( mass[1], width[1], sa[1], sb[1], sc[1], rRes, 1 );
608 }
609 if ( g == 0 ) pro = 1;
610 B[1] = barrier( 1, sa[1], sb[1], sc[1], rRes );
611 B[2] = barrier( 1, sa[2], sb[2], sc[2], rD );
612 calt1( P1, P2, t1V );
613 calt1( pS, pV, t1D );
614 for ( int i = 0; i != 4; i++ ) { temp_PDF += G[i][i] * t1D[i] * t1V[i]; }
615 amp_PDF = temp_PDF * B[1] * B[2] * pro;
616 if ( flag == 2 ) amp_PDF *= KPiSFormfactor( sa[0], sb[0], sc[0], rRes );
617 return amp_PDF;
618}
619EvtComplex EvtD0ToKpipipi::D2TS( double P1[], double P2[], double P3[], double P4[],
620 int flag ) {
621 // flag = 1, T = K*; flag = 2, T = rho
622 double temp_PDF = 0;
623 EvtComplex amp_PDF( 0, 0 );
624 double sa[3], sb[3], sc[3], B[3];
625 double t2D[4][4], t2T[4][4];
626 double pS[4], pT[4], pD[4];
627 for ( int i = 0; i != 4; i++ )
628 {
629 pS[i] = P3[i] + P4[i];
630 pT[i] = P1[i] + P2[i];
631 pD[i] = pT[i] + pS[i];
632 }
633 sa[0] = dot( pT, pT );
634 sb[0] = dot( P1, P1 );
635 sc[0] = dot( P2, P2 );
636 sa[1] = dot( pS, pS );
637 sb[1] = dot( P3, P3 );
638 sc[1] = dot( P4, P4 );
639 sa[2] = dot( pD, pD );
640 sb[2] = sa[0];
641 sc[2] = sa[1];
642 B[0] = barrier( 2, sa[0], sb[0], sc[0], rRes );
643 B[2] = barrier( 2, sa[2], sb[2], sc[2], rD );
644 calt2( P1, P2, t2T );
645 calt2( pT, pS, t2D );
646 for ( int i = 0; i != 4; i++ )
647 {
648 for ( int j = 0; j != 4; j++ )
649 { temp_PDF += t2D[i][j] * t2T[j][i] * ( G[i][i] ) * ( G[j][j] ); }
650 }
651 amp_PDF = temp_PDF * B[0] * B[2];
652 if ( flag == 2 ) amp_PDF *= KPiSFormfactor( sa[1], sb[1], sc[1], rRes );
653 return amp_PDF;
654}
655EvtComplex EvtD0ToKpipipi::PHSP( double Km[], double Pip[] ) {
656 EvtComplex amp_PDF( 0, 0 );
657 double sa, sb, sc;
658 double KPi[4];
659 for ( int i = 0; i != 4; i++ ) { KPi[i] = Km[i] + Pip[i]; }
660 sa = dot( KPi, KPi );
661 sb = dot( Km, Km );
662 sc = dot( Pip, Pip );
663 amp_PDF = KPiSFormfactor( sa, sb, sc, rRes );
664 return amp_PDF;
665}
666double EvtD0ToKpipipi::calDalEva( double P1[], double P2[], double P3[] ) {
667 // Ks Pi0 Pi0 exchange Pi0 in this case
668 // CLEOc model
669 EvtComplex PDF[7], cof, pdf, module;
670 // Ks Pi0 Pi0
671 double P[3][4];
672 for ( int i = 0; i < 4; i++ )
673 {
674 P[0][i] = P1[i];
675 P[1][i] = P2[i];
676 P[2][i] = P3[i];
677 }
678 PDF[0] = Spin_factor( P[1], P[2], P[0], 0 ) + Spin_factor( P[2], P[1], P[0], 0 );
679 PDF[1] = Spin_factor( P[1], P[2], P[0], 10 ) + Spin_factor( P[2], P[1], P[0], 10 );
680 PDF[2] = Spin_factor( P[1], P[2], P[0], 100 ) + Spin_factor( P[2], P[1], P[0], 100 );
681 //-----------------
682 PDF[3] = Spin_factor( P[0], P[1], P[2], 1 ) + Spin_factor( P[0], P[2], P[1], 1 );
683 PDF[4] = Spin_factor( P[0], P[1], P[2], 11 ) + Spin_factor( P[0], P[2], P[1], 11 );
684 //-----------------
685 PDF[5] = Spin_factor( P[1], P[2], P[0], 2 ) + Spin_factor( P[2], P[1], P[0], 2 );
686 PDF[6] = Spin_factor( P[0], P[1], P[2], 12 ) + Spin_factor( P[0], P[2], P[1], 12 );
687
688 rho[0] = 4.2551e-01;
689 phi[0] = 3.7370e+00;
690 rho[1] = -1.0086e+00;
691 phi[1] = 6.9153e+00;
692 rho[2] = 6.5290e-01;
693 phi[2] = 2.8041e+00;
694 rho[3] = 1.0;
695 phi[3] = 0.0;
696 rho[4] = -2.5079e-01;
697 phi[4] = 6.3740e-01;
698 rho[5] = -4.1640e-01;
699 phi[5] = 7.7850e+00;
700 rho[6] = 2.1281e-01;
701 phi[6] = 5.3474e+00;
702
703 pdf = EvtComplex( 0, 0 );
704 for ( int i = 0; i < 7; i++ )
705 {
706 cof = EvtComplex( rho[i] * cos( phi[i] ), rho[i] * sin( phi[i] ) );
707 pdf += cof * PDF[i];
708 }
709 module = conj( pdf ) * pdf;
710 double value;
711 value = real( module );
712 return ( value <= 0 ) ? 1e-20 : value;
713}
714EvtComplex EvtD0ToKpipipi::Spin_factor( double P1[], double P2[], double P3[], int spin ) {
715 // D-> R P3, R->P1 P2
716 double R[4], s[3], sp2, sD, B[2];
717 for ( int i = 0; i < 4; i++ ) { R[i] = P1[i] + P2[i]; }
718 s[0] = dot( R, R );
719 s[1] = dot( P1, P1 );
720 s[2] = dot( P2, P2 );
721 sp2 = dot( P3, P3 );
722 sD = mD * mD;
723 EvtComplex amp( 0, 0 ), prop;
724 prop = getProp( s, spin );
725 double tmp( 0. );
726 if ( spin % 10 == 0 ) { amp = prop; }
727 else if ( spin % 10 == 1 )
728 {
729 tmp = 0;
730 double T1[4], t1[4];
731 calt1( R, P3, T1 );
732 calt1( P1, P2, t1 );
733 B[0] = barrier( 1, s[0], s[1], s[2], 3.0 );
734 B[1] = barrier( 1, sD, s[0], sp2, 3.0 );
735 for ( int i = 0; i < 4; i++ ) { tmp += T1[i] * t1[i] * G[i][i]; }
736 amp = tmp * prop * B[0] * B[1];
737 }
738 else if ( spin % 10 == 2 )
739 {
740 tmp = 0;
741 double T2[4][4], t2[4][4];
742 calt2( R, P3, T2 );
743 calt2( P1, P2, t2 );
744 B[0] = barrier( 2, s[0], s[1], s[2], 3.0 );
745 B[1] = barrier( 2, sD, s[0], sp2, 3.0 );
746 for ( int i = 0; i < 4; i++ )
747 {
748 for ( int j = 0; j < 4; j++ ) { tmp += T2[i][j] * t2[j][i] * G[j][j] * G[i][i]; }
749 }
750 amp = tmp * prop * B[0] * B[1];
751 }
752 else { cout << "Only S, P, D wave allowed" << endl; }
753 return amp;
754}
755EvtComplex EvtD0ToKpipipi::getProp( double s[], int flag ) {
756 // Ks Pi0 Pi0
757 // 0 980
758 // 10 f0_1370
759 // 20 f0_1500
760 // 2 f2_1270
761 // 1 K*892
762 // 12 K*2_1430
763 // 11 K*1680
764 // 100 f0(500)
765 EvtComplex prop( 1, 0 );
766 double snpi, scpi, snk, sck;
767 snpi = mass_Pi0 * mass_Pi0;
768 scpi = mass_Pion * mass_Pion;
769 sck = mass_Kaon * mass_Kaon;
770 snk = mk0 * mk0;
771 if ( flag == 0 )
772 {
773 double mass, gpipi, gkk;
774 mass = 0.965;
775 gpipi = 0.406;
776 gkk = 2 * gpipi;
777 EvtComplex ci( 0, 1 );
778 EvtComplex rhokk, rhopipi;
779 rhokk = 0.5 * ( rhofactor( s[0], sck ) + rhofactor( s[0], snk ) );
780 rhopipi = 1.0 / 3.0 * ( 2.0 * rhofactor( s[0], scpi ) + rhofactor( s[0], snpi ) );
781 prop = 1.0 / ( mass * mass - s[0] - ci * ( gpipi * gpipi * rhopipi + gkk * gkk * rhokk ) );
782 }
783 if ( flag == 10 ) { prop = propagatorRBW( 1.35, 0.265, s[0], s[1], s[2], 3.0, 0 ); }
784 if ( flag == 20 ) { prop = propagatorRBW( 1.505, 0.109, s[0], s[1], s[2], 3.0, 0 ); }
785 if ( flag == 1 ) { prop = propagatorRBW( 0.896, 0.0503, s[0], s[1], s[2], 3.0, 1 ); }
786 if ( flag == 11 ) { prop = propagatorRBW( 1.717, 0.322, s[0], s[1], s[2], 3.0, 1 ); }
787 if ( flag == 2 ) { prop = propagatorRBW( 1.2751, 0.185, s[0], s[1], s[2], 3.0, 2 ); }
788 if ( flag == 12 ) { prop = propagatorRBW( 1.4324, 0.109, s[0], s[1], s[2], 3.0, 2 ); }
789 if ( flag == 100 )
790 {
791 EvtComplex pole( 0.470, -0.220 );
792 prop = 1.0 / ( s[0] - pole * pole );
793 }
794 return prop;
795}
796EvtComplex EvtD0ToKpipipi::rhofactor( double sx, double sdau ) {
797 EvtComplex one( 1, 0 );
798 EvtComplex ci( 0, 1 );
799 EvtComplex res;
800 double tmp;
801 tmp = 1 - 4 * sdau / sx;
802 if ( tmp > 0 ) res = one * sqrt( tmp );
803 if ( tmp < 0 ) res = ci * sqrt( fabs( tmp ) );
804 return res;
805}
806EvtComplex EvtD0ToKpipipi::propogator( double mass, double width, double sx ) const {
807 EvtComplex ci( 0, 1 );
808 EvtComplex prop = 1.0 / ( mass * mass - sx - ci * mass * width );
809 return prop;
810}
811double EvtD0ToKpipipi::wid( double mass, double sa, double sb, double sc, double r,
812 int l ) const {
813 double widm( 0. ), q( 0. ), q0( 0. );
814 double sa0 = mass * mass;
815 double m = sqrt( sa );
816 q = Qabcs( sa, sb, sc );
817 q0 = Qabcs( sa0, sb, sc );
818 double z, z0;
819 z = q * r * r;
820 z0 = q0 * r * r;
821 double t = q / q0;
822 double F( 0. );
823 if ( l == 0 ) F = 1;
824 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
825 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
826 widm = pow( t, l + 0.5 ) * mass / m * F * F;
827 return widm;
828}
829double EvtD0ToKpipipi::h( double m, double q ) const {
830 double h( 0. );
831 h = 2 / pi * q / m * log( ( m + 2 * q ) / ( 2 * mpi ) );
832 return h;
833}
834double EvtD0ToKpipipi::dh( double mass, double q0 ) const {
835 double dh = h( mass, q0 ) * ( 1.0 / ( 8 * q0 * q0 ) - 1.0 / ( 2 * mass * mass ) ) +
836 1.0 / ( 2 * pi * mass * mass );
837 return dh;
838}
839double EvtD0ToKpipipi::f( double mass, double sx, double q0, double q ) const {
840 double m = sqrt( sx );
841 double f = mass * mass / ( pow( q0, 3 ) ) *
842 ( q * q * ( h( m, q ) - h( mass, q0 ) ) +
843 ( mass * mass - sx ) * q0 * q0 * dh( mass, q0 ) );
844 return f;
845}
846double EvtD0ToKpipipi::d( double mass, double q0 ) const {
847 double d = 3.0 / pi * mpi * mpi / ( q0 * q0 ) * log( ( mass + 2 * q0 ) / ( 2 * mpi ) ) +
848 mass / ( 2 * pi * q0 ) - ( mpi * mpi * mass ) / ( pi * pow( q0, 3 ) );
849 return d;
850}
851EvtComplex EvtD0ToKpipipi::propagatorRBW( double mass, double width, double sa, double sb,
852 double sc, double r, int l ) const {
853 EvtComplex ci( 0, 1 );
854 EvtComplex prop =
855 1.0 / ( mass * mass - sa - ci * mass * width * wid( mass, sa, sb, sc, r, l ) );
856 return prop;
857}
858EvtComplex EvtD0ToKpipipi::propagatorGS( double mass, double width, double sa, double sb,
859 double sc, double r, int l ) const {
860 EvtComplex ci( 0, 1 );
861 double q = Qabcs( sa, sb, sc );
862 double sa0 = mass * mass;
863 double q0 = Qabcs( sa0, sb, sc );
864 q = sqrt( q );
865 q0 = sqrt( q0 );
866 EvtComplex prop = ( 1 + d( mass, q0 ) * width / mass ) /
867 ( mass * mass - sa + width * f( mass, sa, q0, q ) -
868 ci * mass * width * wid( mass, sa, sb, sc, r, l ) );
869 return prop;
870}
871double EvtD0ToKpipipi::Flatte_rhoab( double sa, double sb, double sc ) const {
872 double q = Qabcs( sa, sb, sc );
873 double rho = sqrt( q / sa );
874 return rho;
875}
876EvtComplex EvtD0ToKpipipi::propagatorFlatte( double mass, double width, double sx, double* sb,
877 double* sc ) const {
878 EvtComplex ci( 0, 1 );
879 double rho1 = Flatte_rhoab( sx, sb[0], sc[0] );
880 double rho2 = Flatte_rhoab( sx, sb[1], sc[1] );
881 EvtComplex prop = 1.0 / ( mass * mass - sx - ci * ( g1 * g1 * rho1 + g2 * g2 * rho2 ) );
882 return prop;
883}
884double EvtD0ToKpipipi::dot( double* a1, double* a2 ) const {
885 double dot = 0;
886 for ( int i = 0; i != 4; i++ ) { dot += a1[i] * a2[i] * G[i][i]; }
887 return dot;
888}
889double EvtD0ToKpipipi::Qabcs( double sa, double sb, double sc ) const {
890 double Qabcs = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
891 if ( Qabcs < 0 ) Qabcs = 1e-16;
892 return Qabcs;
893}
894double EvtD0ToKpipipi::barrier( double l, double sa, double sb, double sc, double r ) const {
895 double q = Qabcs( sa, sb, sc );
896 q = sqrt( q );
897 double z = q * r;
898 z = z * z;
899 double F = 1;
900 if ( l > 2 ) F = 0;
901 if ( l == 0 ) F = 1;
902 if ( l == 1 ) { F = sqrt( ( 2 * z ) / ( 1 + z ) ); }
903 if ( l == 2 ) { F = sqrt( ( 13 * z * z ) / ( 9 + 3 * z + z * z ) ); }
904 return F;
905}
906void EvtD0ToKpipipi::calt1( double daug1[], double daug2[], double t1[] ) const {
907 double p, pq;
908 double pa[4], qa[4];
909 for ( int i = 0; i != 4; i++ )
910 {
911 pa[i] = daug1[i] + daug2[i];
912 qa[i] = daug1[i] - daug2[i];
913 }
914 p = dot( pa, pa );
915 pq = dot( pa, qa );
916 for ( int i = 0; i != 4; i++ ) { t1[i] = qa[i] - pq / p * pa[i]; }
917}
918void EvtD0ToKpipipi::calt2( double daug1[], double daug2[], double t2[][4] ) const {
919 double p, r;
920 double pa[4], t1[4];
921 calt1( daug1, daug2, t1 );
922 r = dot( t1, t1 );
923 for ( int i = 0; i != 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
924 p = dot( pa, pa );
925 for ( int i = 0; i != 4; i++ )
926 {
927 for ( int j = 0; j != 4; j++ )
928 { t2[i][j] = t1[i] * t1[j] - 1.0 / 3 * r * ( G[i][j] - pa[i] * pa[j] / p ); }
929 }
930}
P
double P(RecMdcKalTrack *trk)
Definition CalibEventSelect.cxx:84
mass
double mass
Definition CosmicGenerator.cxx:128
pi2
character *LEPTONflag integer iresonances real pi2
Definition Ekhara/src/ekhara/common.h:3
mass_Pion
const double mass_Pion
Definition EvtD0ToKSpipipi0.cc:1561
s
XmlRpcServer s
Definition HelloServer.cpp:11
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
EvtD0ToKpipipi::clone
EvtDecayBase * clone()
Definition EvtD0ToKpipipi.cc:35
EvtD0ToKpipipi::getName
void getName(std::string &name)
Definition EvtD0ToKpipipi.cc:33
EvtD0ToKpipipi::decay
void decay(EvtParticle *p)
Definition EvtD0ToKpipipi.cc:164
EvtD0ToKpipipi::~EvtD0ToKpipipi
virtual ~EvtD0ToKpipipi()
Definition EvtD0ToKpipipi.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::R
complex_t R(double Q2, double M2, double G, double Mp2, double Mm2)
Definition TUtil.h:22
t
int t()
Definition t.c:1