BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToKKpipi0.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: EvtDsToKKpipi0.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 "EvtDsToKKpipi0.hh"
18#include "../EvtGenBase/EvtConst.hh"
19#include "../EvtGenBase/EvtDecayTable.hh"
20#include "../EvtGenBase/EvtFlatte.hh"
21#include "../EvtGenBase/EvtGenKine.hh"
22#include "../EvtGenBase/EvtPDL.hh"
23#include "../EvtGenBase/EvtParticle.hh"
24#include "../EvtGenBase/EvtPatches.hh"
25#include "../EvtGenBase/EvtReport.hh"
26#include "../EvtGenBase/EvtResonance.hh"
27#include "../EvtGenBase/EvtResonance2.hh"
28#include <fstream>
29#include <iomanip>
30#include <stdlib.h>
31#include <string>
32using std::endl;
33
34EvtDsToKKpipi0::~EvtDsToKKpipi0() {}
35
36void EvtDsToKKpipi0::getName( std::string& model_name ) { model_name = "DsToKKpipi0"; }
37
38EvtDecayBase* EvtDsToKKpipi0::clone() { return new EvtDsToKKpipi0; }
39
40void EvtDsToKKpipi0::init() {
41 // check that there are 0 arguments
42 checkNArg( 0 );
43 checkNDaug( 4 );
44
45 checkSpinParent( EvtSpinType::SCALAR );
46 checkSpinDaughter( 0, EvtSpinType::SCALAR );
47 checkSpinDaughter( 1, EvtSpinType::SCALAR );
48 checkSpinDaughter( 2, EvtSpinType::SCALAR );
49 checkSpinDaughter( 3, EvtSpinType::SCALAR );
50
51 int mode = 0;
52 //--------------------amp---------------
53 phi[1] = -1.45750737501177241029;
54 phi[2] = 1.45694587813248066510;
55 phi[3] = -2.14539596963664980223;
56 phi[4] = -0.51808334601975225553;
57 phi[5] = -1.56645961468228378521;
58 phi[6] = 1.87207969874086810336;
59 phi[7] = -0.92156758917410819265;
60 phi[8] = -0.92156758917410819265;
61 phi[9] = 0.61333399627627738226;
62 phi[10] = 2.95348100053888362737;
63 phi[11] = 2.12806871473955627749;
64 phi[12] = 2.12806871473955627749;
65
66 phi[13] = 2.14745266587422545257;
67 phi[14] = -0.25044816999225272269;
68 phi[15] = -0.25044816999225272269;
69 phi[16] = 1.52246676387907431405;
70 phi[17] = 1.52246676387907431405;
71 //-------------------phase--------------
72 rho[1] = 0.14853586409145513869;
73 rho[2] = 0.06437581911225187525;
74 rho[3] = 0.63500586659756663721;
75 rho[4] = 0.11892962207932455954;
76 rho[5] = 0.06386739319734147102;
77 rho[6] = 1.72035932015034198628;
78 rho[7] = 0.84266300186759934832;
79 rho[8] = 0.84266300186759934832;
80 rho[9] = 5.52084783967448444741;
81 rho[10] = 1.18304173083694408319;
82 rho[11] = 0.37041606275538363491;
83 rho[12] = 0.37041606275538363491;
84
85 rho[13] = 1.99906320501213841112;
86 rho[14] = 0.28474658039299072243;
87 rho[15] = 0.28474658039299072243;
88 rho[16] = 0.10589240624900675414;
89 rho[17] = 0.10589240624900675414;
90 //--------------------mass---------------------
91 mass1[0] = 1.019461;
92 mass1[1] = 1.019461;
93 mass1[2] = 1.019461;
94 mass2[0] = 0.77511;
95 mass2[1] = 0.77511;
96 mass2[2] = 0.77511;
97 mass2[6] = 0.77511;
98 mass2[13] = 0.77511;
99 //---------------------
100 mass1[3] = 0.89555;
101 mass1[4] = 0.89555;
102 mass1[5] = 0.89555;
103 mass1[7] = 0.89555;
104 mass1[14] = 0.89555;
105 mass1[16] = 0.89555;
106 mass2[3] = 0.89166;
107 mass2[4] = 0.89166;
108 mass2[5] = 0.89166;
109 mass1[12] = 0.89166;
110 //----------------------
111 mass1[6] = 1.272;
112 mass2[14] = 1.272;
113 mass2[15] = 1.272;
114 mass2[16] = 1.272;
115 mass2[17] = 1.272;
116 //-----------------------
117 mass2[7] = 1.403;
118 mass2[8] = 1.403;
119 mass1[8] = 0.89166;
120 mass1[11] = 0.89166;
121 mass1[15] = 0.89166;
122 mass1[17] = 0.89166;
123 //----------------------
124 mass2[9] = 1.475;
125 mass1[9] = 0.99;
126 mass1[10] = 0.99;
127 mass1[13] = 0.99;
128 //---------------------
129 mass2[10] = 1.4263;
130 mass2[11] = 1.4263;
131 mass2[12] = 1.4263;
132 //---------------------width---------------------
133 width1[0] = 0.004249;
134 width1[1] = 0.004249;
135 width1[2] = 0.004249;
136 width2[0] = 0.1491;
137 width2[1] = 0.1491;
138 width2[2] = 0.1491;
139 width2[6] = 0.1491;
140 width2[13] = 0.1491;
141 //----------------------
142 width1[3] = 0.0473;
143 width1[4] = 0.0473;
144 width1[5] = 0.0473;
145 width1[7] = 0.0473;
146 width1[14] = 0.0473;
147 width1[16] = 0.0473;
148 width2[3] = 0.0508;
149 width2[4] = 0.0508;
150 width2[5] = 0.0508;
151 width1[12] = 0.0508;
152 //----------------------
153 width1[6] = 0.087;
154 width2[14] = 0.087;
155 width2[15] = 0.087;
156 width2[16] = 0.087;
157 width2[17] = 0.087;
158 //-----------------------
159 width2[7] = 0.174;
160 width2[8] = 0.174;
161 width1[8] = 0.0508;
162 width1[11] = 0.0508;
163 width1[15] = 0.0508;
164 width1[17] = 0.0508;
165 //------------------------
166 width2[9] = 0.09;
167 width1[9] = 0.08;
168 width1[10] = 0.08;
169 width1[13] = 0.08;
170 //---------------------
171 width2[10] = 0.0545;
172 width2[11] = 0.0545;
173 width2[12] = 0.0545;
174
175 modetype[0] = 1;
176 modetype[1] = 1;
177 modetype[2] = 1;
178 modetype[3] = 2;
179 modetype[4] = 2;
180 modetype[5] = 2;
181 modetype[6] = 4;
182 modetype[7] = 3;
183 modetype[8] = 333;
184 modetype[9] = 30;
185 modetype[10] = 31;
186 modetype[11] = 73;
187 modetype[12] = 731;
188 modetype[13] = 50;
189 modetype[14] = 3;
190 modetype[15] = 333;
191 modetype[16] = 3;
192 modetype[17] = 333;
193
194 phi[0] = 0;
195 rho[0] = 1;
196
197 // cout << "DsToKKpipi0 :" << endl;
198 // for (int i=0; i<18; i++) {
199 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << " m1= " << mass1[i] << " m2= "
200 // << mass2[i] << " w1= "
201 // << width1[i] << " w2= " << width2[i] << endl;
202 // }
203
204 mass_Pion = 0.13957;
205 mass_Pion_N = 0.134977;
206 mass_Eta = 0.547862;
207 math_pi = 3.1415926;
208 rD2 = 25.0; // 5*5
209 rRes2 = 9.0; // 3*3
210 GS1 = 0.636619783;
211 GS2 = 0.01860182466;
212 GS3 = 0.1591549458; // 1/(2*math_2pi)
213 GS4 = 0.00620060822; // mass_Pion2/math_pi
214
215 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
216 int EE[4][4][4][4] = {
217 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
218 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, -1, 0 } },
219 { { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 1, 0, 0 } },
220 { { 0, 0, 0, 0 }, { 0, 0, 1, 0 }, { 0, -1, 0, 0 }, { 0, 0, 0, 0 } } },
221 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 1, 0 } },
222 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
223 { { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 } },
224 { { 0, 0, -1, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 } } },
225 { { { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, -1, 0, 0 } },
226 { { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 } },
227 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
228 { { 0, 1, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } },
229 { { { 0, 0, 0, 0 }, { 0, 0, -1, 0 }, { 0, 1, 0, 0 }, { 0, 0, 0, 0 } },
230 { { 0, 0, 1, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 } },
231 { { 0, -1, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
232 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } } };
233 for ( int i = 0; i < 4; i++ )
234 {
235 for ( int j = 0; j < 4; j++ )
236 {
237 G[i][j] = GG[i][j];
238 for ( int k = 0; k < 4; k++ )
239 {
240 for ( int l = 0; l < 4; l++ ) { E[i][j][k][l] = EE[i][j][k][l]; }
241 }
242 }
243 }
244}
245
246void EvtDsToKKpipi0::initProbMax() {
247 // setProbMax(43874577.5);
248 setProbMax( 44000000.0 );
249}
250
251void EvtDsToKKpipi0::decay( EvtParticle* p ) {
252 /*
253 double maxprob = 0.0;
254 for(int ir=0;ir<=60000000;ir++){
255 p->initializePhaseSpace(getNDaug(),getDaugs());
256 EvtVector4R _km = p->getDaug(0)->getP4();
257 EvtVector4R _kp = p->getDaug(1)->getP4();
258 EvtVector4R _pip = p->getDaug(2)->getP4();
259 EvtVector4R _pi0 = p->getDaug(3)->getP4();
260
261 double _Km[4], _Kp[4], _Pip[4], _Pi0[4];
262 _Km[0] = _km.get(0); _Kp[0] = _kp.get(0); _Pip[0] = _pip.get(0); _Pi0[0] = _pi0.get(0);
263 _Km[1] = _km.get(1); _Kp[1] = _kp.get(1); _Pip[1] = _pip.get(1); _Pi0[1] = _pi0.get(1);
264 _Km[2] = _km.get(2); _Kp[2] = _kp.get(2); _Pip[2] = _pip.get(2); _Pi0[2] = _pi0.get(2);
265 _Km[3] = _km.get(3); _Kp[3] = _kp.get(3); _Pip[3] = _pip.get(3); _Pi0[3] = _pi0.get(3);
266
267 double _prob;
268 int g0[18]={1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1};
269 int g1[18]={1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1};
270 int g2[18]={0,1,2,0,1,2,0,0,0,0,1,0,0,1,0,0,2,2};
271 int nstates=18;
272 calEvaMy(_Km,_Kp,_Pip,_Pi0,mass1,mass2,width1,width2,rho,phi,g0,g1,g2,modetype,nstates,_prob);
273 if(_prob>maxprob) {
274 maxprob=_prob;
275 cout << "Max PDF = " << ir << " " << _prob << endl;
276 }
277 }
278 cout << "Max!!!!!!!!!!! " << maxprob<< endl;
279 // setProbMax(maxprob);
280 */
281 p->initializePhaseSpace( getNDaug(), getDaugs() );
282 EvtVector4R km = p->getDaug( 0 )->getP4();
283 EvtVector4R kp = p->getDaug( 1 )->getP4();
284 EvtVector4R pip = p->getDaug( 2 )->getP4();
285 EvtVector4R pi0 = p->getDaug( 3 )->getP4();
286
287 double Km[4], Kp[4], Pip[4], Pi0[4];
288 Km[0] = km.get( 0 );
289 Kp[0] = kp.get( 0 );
290 Pip[0] = pip.get( 0 );
291 Pi0[0] = pi0.get( 0 );
292 Km[1] = km.get( 1 );
293 Kp[1] = kp.get( 1 );
294 Pip[1] = pip.get( 1 );
295 Pi0[1] = pi0.get( 1 );
296 Km[2] = km.get( 2 );
297 Kp[2] = kp.get( 2 );
298 Pip[2] = pip.get( 2 );
299 Pi0[2] = pi0.get( 2 );
300 Km[3] = km.get( 3 );
301 Kp[3] = kp.get( 3 );
302 Pip[3] = pip.get( 3 );
303 Pi0[3] = pi0.get( 3 );
304
305 double prob;
306 int g0[18] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
307 int g1[18] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
308 int g2[18] = { 0, 1, 2, 0, 1, 2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 2, 2 };
309 int nstates = 18;
310 calEvaMy( Km, Kp, Pip, Pi0, mass1, mass2, width1, width2, rho, phi, g0, g1, g2, modetype,
311 nstates, prob );
312
313 setProb( prob );
314 return;
315}
316
317void EvtDsToKKpipi0::Com_Multi( double a1[2], double a2[2], double res[2] ) {
318 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
319 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
320}
321void EvtDsToKKpipi0::Com_Divide( double a1[2], double a2[2], double res[2] ) {
322 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
323 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
324 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
325}
326double EvtDsToKKpipi0::SCADot( double a1[4], double a2[4] ) {
327 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
328 return _cal;
329}
330double EvtDsToKKpipi0::barrier( int l, double sa, double sb, double sc, double r2 ) {
331 double F;
332 double tmp = sa + sb - sc;
333 double q = 0.25 * tmp * tmp / sa - sb;
334 if ( q < 0 ) q = fabs( q );
335 // double z = q*r2;
336 if ( l == 1 ) { F = sqrt( 2.0 / ( 1.0 / r2 + q ) ); }
337 else if ( l == 2 )
338 {
339 // double z2 = z*z;
340 // F = sqrt(13.0*z2/(9.0+3.0*z+z2));
341 F = sqrt( 13.0 / ( 9.0 * ( 1 / ( r2 * r2 ) ) + 3.0 * q * ( 1 / r2 ) + q * q ) );
342 }
343 else { F = 1.0; }
344 return F;
345}
346void EvtDsToKKpipi0::calt1( double daug1[4], double daug2[4], double t1[4] ) {
347 double p, pq, tmp;
348 double pa[4], qa[4];
349 for ( int i = 0; i < 4; i++ )
350 {
351 pa[i] = daug1[i] + daug2[i];
352 qa[i] = daug1[i] - daug2[i];
353 }
354 p = SCADot( pa, pa );
355 pq = SCADot( pa, qa );
356 tmp = pq / p;
357 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
358}
359void EvtDsToKKpipi0::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
360 double p, r;
361 double pa[4], t1[4];
362 calt1( daug1, daug2, t1 );
363 r = SCADot( t1, t1 ) / 3.0;
364 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
365 p = SCADot( pa, pa );
366 for ( int i = 0; i < 4; i++ )
367 {
368 for ( int j = 0; j < 4; j++ )
369 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
370 }
371}
372void EvtDsToKKpipi0::Flatte_rhoab( double sa, double sb, double sc, double rho[2] ) {
373 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
374 if ( q > 0 )
375 {
376 rho[0] = 2 * sqrt( q / sa );
377 rho[1] = 0;
378 }
379 else if ( q < 0 )
380 {
381 rho[0] = 0;
382 rho[1] = 2 * sqrt( -q / sa );
383 }
384}
385void EvtDsToKKpipi0::propagatora0980( double mass2, double mass, double sx, double* sb,
386 double* sc, double prop[2] ) {
387 double unit[2] = { 1.0 };
388 double ci[2] = { 0, 1 };
389 double rho1[2];
390 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
391 double rho2[2];
392 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
393
394 double gKK_a980 = 0.892 * 0.341, gPiEta_a980 = 0.341;
395 double tmp1[2] = { gKK_a980, 0 };
396 double tmp11[2];
397 double tmp2[2] = { gPiEta_a980, 0 };
398 double tmp22[2];
399 Com_Multi( tmp1, rho1, tmp11 );
400 Com_Multi( tmp2, rho2, tmp22 );
401 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
402 double tmp31[2];
403 Com_Multi( tmp3, ci, tmp31 );
404 double tmp4[2] = { mass2 - sx - tmp31[0], -1.0 * tmp31[1] };
405 Com_Divide( unit, tmp4, prop );
406}
407double EvtDsToKKpipi0::wid( double mass2, double mass, double sa, double sb, double sc,
408 double r2, int l ) {
409 double widm = 0.;
410 double m = sqrt( sa );
411 double tmp = sb - sc;
412 double tmp1 = sa + tmp;
413 double q = 0.25 * tmp1 * tmp1 / sa - sb;
414 if ( q < 0 ) q = fabs( q );
415 double tmp2 = mass2 + tmp;
416 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
417 if ( q0 < 0 ) q0 = fabs( q0 );
418 double z = q * r2;
419 double z0 = q0 * r2;
420 double t = q / q0;
421 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
422 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
423 else if ( l == 2 )
424 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
425 return widm;
426}
427double EvtDsToKKpipi0::widl1( double mass2, double mass, double sa, double sb, double sc,
428 double r2 ) {
429 double widm = 0.;
430 double m = sqrt( sa );
431 double tmp = sb - sc;
432 double tmp1 = sa + tmp;
433 double q = 0.25 * tmp1 * tmp1 / sa - sb;
434 if ( q < 0 ) q = fabs( q );
435 double tmp2 = mass2 + tmp;
436 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
437 if ( q0 < 0 ) q0 = fabs( q0 );
438 double z = q * r2;
439 double z0 = q0 * r2;
440 double F = ( 1 + z0 ) / ( 1 + z );
441 double t = q / q0;
442 widm = t * sqrt( t ) * mass / m * F;
443 return widm;
444}
445void EvtDsToKKpipi0::propagatorRBW( double mass2, double mass, double width, double sa,
446 double sb, double sc, double r2, int l, double prop[2] ) {
447 double a[2], b[2];
448 a[0] = 1;
449 a[1] = 0;
450 b[0] = mass2 - sa;
451 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
452 Com_Divide( a, b, prop );
453}
454void EvtDsToKKpipi0::propagatorNBW( double mass2, double mass, double width, double sa,
455 double sb, double sc, double r2, int l, double prop[2] ) {
456 double a[2], b[2];
457 a[0] = 1;
458 a[1] = 0;
459 b[0] = mass2 - sa;
460 b[1] = -mass * width;
461 Com_Divide( a, b, prop );
462}
463void EvtDsToKKpipi0::propagatorRBWl1( double mass2, double mass, double width, double sa,
464 double sb, double sc, double r2, double prop[2] ) {
465 double a[2], b[2];
466 a[0] = 1;
467 a[1] = 0;
468 b[0] = mass2 - sa;
469 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
470 Com_Divide( a, b, prop );
471}
472void EvtDsToKKpipi0::propagatorGS( double mass2, double mass, double width, double sa,
473 double sb, double sc, double r2, double prop[2] ) {
474 double a[2], b[2];
475 double tmp = sb - sc;
476 double tmp1 = sa + tmp;
477 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
478 if ( q2 < 0 ) q2 = fabs( q2 );
479
480 double tmp2 = mass2 + tmp;
481 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
482 if ( q02 < 0 ) q02 = fabs( q02 );
483
484 double q = sqrt( q2 );
485 double q0 = sqrt( q02 );
486 double m = sqrt( sa );
487 double q03 = q0 * q02;
488 double tmp3 = log( mass + 2 * q0 ) + 1.2760418309; // log(mass_2Pion) = 1.2760418309;
489
490 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2760418309 );
491 double h0 = GS1 * q0 / mass * tmp3;
492 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
493 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
494 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
495
496 a[0] = 1.0 + d * width / mass;
497 a[1] = 0.0;
498 b[0] = mass2 - sa + width * f;
499 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
500 Com_Divide( a, b, prop );
501}
502void EvtDsToKKpipi0::KPiSLASS( double sa, double sb, double sc, double prop[2] ) {
503 const double m1430 = 1.441;
504 const double sa0 = 2.076481; // m1430*m1430;
505 const double w1430 = 0.193;
506 const double Lass1 = 0.25 / sa0;
507 double tmp = sb - sc;
508 double tmp1 = sa0 + tmp;
509 double q0 = Lass1 * tmp1 * tmp1 - sb;
510 if ( q0 < 0 ) q0 = 1e-16;
511 double tmp2 = sa + tmp;
512 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
513 double q = sqrt( qs );
514 double width = w1430 * q * m1430 / sqrt( sa * q0 );
515 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
516 if ( temp_R < 0 ) temp_R += math_pi;
517 double deltaR = -109.7 * math_pi / 180.0 + temp_R;
518 double temp_F =
519 atan( 0.226 * q / ( 2.0 - 3.8194 * qs ) ); // 2.0*1.07 = 2.14; 1.8*1.07 = 1.926
520 if ( temp_F < 0 ) temp_F += math_pi;
521 double deltaF = 0.1 * math_pi / 180.0 + temp_F;
522 double deltaS = deltaR + 2.0 * deltaF;
523 double t1 = 0.96 * sin( deltaF );
524 double t2 = sin( deltaR );
525 double CF[2], CS[2];
526 CF[0] = cos( deltaF );
527 CF[1] = sin( deltaF );
528 CS[0] = cos( deltaS );
529 CS[1] = sin( deltaS );
530 prop[0] = t1 * CF[0] + t2 * CS[0];
531 prop[1] = t1 * CF[1] + t2 * CS[1];
532}
533
534void EvtDsToKKpipi0::calEvaMy( double* Km, double* Kp, double* Pip, double* Pi0, double* mass1,
535 double* mass2, double* width1, double* width2, double* amp,
536 double* phase, int* g0, int* g1, int* g2, int* modetype,
537 int nstates, double& Result ) {
538 double cof[2], amp_tmp[2], amp_PDF[2], PDF[2];
539 Km[0] = sqrt( 0.243719942 + Km[1] * Km[1] + Km[2] * Km[2] + Km[3] * Km[3] );
540 Kp[0] = sqrt( 0.243719942 + Kp[1] * Kp[1] + Kp[2] * Kp[2] + Kp[3] * Kp[3] );
541 Pip[0] = sqrt( 0.019479785 + Pip[1] * Pip[1] + Pip[2] * Pip[2] + Pip[3] * Pip[3] );
542 Pi0[0] = sqrt( 0.0182187905 + Pi0[1] * Pi0[1] + Pi0[2] * Pi0[2] + Pi0[3] * Pi0[3] );
543 //---------------------------------------------------------------
544 // Km[0] =0.5096197263 ; Km[1] = 0.0200218487 ; Km[2] = 0.1191241019 ; Km[3] =
545 // 0.037428 ; Kp[0] =0.5781433946 ; Kp[1] = 0.0870387478 ; Kp[2] =-0.2596464631 ;
546 // Kp[3] = 0.124650 ; Pip[0] =0.4738800750 ; Pip[1] =-0.3659960295 ; Pip[2]
547 // =-0.0297276846 ; Pip[3] =-0.265039 ; Pi0[0] =0.4584942092 ; Pi0[1] = 0.3530867237
548 // ; Pi0[2] =-0.2591360024 ; Pi0[3] =-0.013287 ;
549
550 double pKstr0[4], pKstrp[4], pKstrm[4], prhop[4], pphi[4], pK10[4], pK11[4], pK12[4],
551 pK13[4], pK14[4], pD[4];
552 for ( int i = 0; i != 4; i++ )
553 {
554 pphi[i] = Km[i] + Kp[i];
555 prhop[i] = Pip[i] + Pi0[i];
556 pKstr0[i] = Km[i] + Pip[i];
557 pKstrp[i] = Kp[i] + Pi0[i];
558 pKstrm[i] = Km[i] + Pi0[i];
559 pK10[i] = pKstrm[i] + Kp[i];
560 pK11[i] = pKstrm[i] + Pip[i];
561 pK12[i] = pKstrp[i] + Km[i];
562 pK13[i] = pKstr0[i] + Pi0[i]; // can delete, no use, equal to pK11
563 pK14[i] = pphi[i] + Pi0[i];
564 pD[i] = pphi[i] + prhop[i];
565 }
566 double spi0, spionp, skaon1, skaon2, sphi, srhop, skstr0, skstrp, skstrm, sk10, sk11, sk12,
567 sk13, sk14, sD;
568 skaon1 = SCADot( Km, Km );
569 skaon2 = SCADot( Kp, Kp );
570 spionp = SCADot( Pip, Pip );
571 spi0 = SCADot( Pi0, Pi0 );
572 sphi = SCADot( pphi, pphi );
573 srhop = SCADot( prhop, prhop );
574 skstr0 = SCADot( pKstr0, pKstr0 );
575 skstrp = SCADot( pKstrp, pKstrp );
576 skstrm = SCADot( pKstrm, pKstrm );
577 sk10 = SCADot( pK10, pK10 );
578 sk11 = SCADot( pK11, pK11 );
579 sk12 = SCADot( pK12, pK12 );
580 sk13 = SCADot( pK13, pK13 );
581 sk14 = SCADot( pK14, pK14 );
582 sD = SCADot( pD, pD );
583 //-----------------------------------------------------------------------------------------------
584 double t1phi[4], t1A1[4], t1A2[4], t1rhop[4], t1kstr1[4], t1kstr2[4], t1kstr3[4], t1kstr4[4],
585 t2k11[4][4], t2k12[4][4], t2k13[4][4], t2k14[4][4], t2k21[4][4];
586 calt1( Km, Kp, t1phi );
587 calt1( Pip, Pi0, t1rhop );
588 calt1( Km, Pip, t1kstr1 );
589 calt1( Kp, Pi0, t1kstr2 );
590 calt1( Km, Pi0, t1kstr3 );
591 calt1( Kp, Pi0, t1kstr4 );
592
593 calt1( pKstr0, Pi0, t1A1 );
594 calt1( pKstrm, Pip, t1A2 );
595
596 calt2( pKstr0, Pi0, t2k11 );
597 calt2( pKstrm, Pip, t2k12 );
598
599 calt2( pKstrm, Kp, t2k13 );
600 calt2( pKstrp, Km, t2k14 );
601 calt2( prhop, Km, t2k21 );
602
603 double B_phi = -1.0, B_rho = -1.0, B_rho_2 = -1.0, B_phirho1 = -1.0, B_phirho2 = -1.0,
604 B_phi_2 = -1.0;
605 double B_Kstp = -1.0, B_Kst0 = -1.0, B_KsKs1 = -1.0, B_KsKs2 = -1.0, B_Kstp_2 = -1.0,
606 B_Kst0_2 = -1.0, B_DA2 = -1.0;
607 double B_Kstm = -1.0, B_DA1 = -1.0, B_A1 = -1.0, B_A2 = -1.0, B_A1_1 = -1.0, B_K13_1 = -1.0,
608 B_A2_1 = -1.0;
609 double B_D_Arho = -1.0, B_Arho = -1.0, B_K11 = -1.0, B_D11 = -1.0, B_Arho1 = -1.0,
610 B_K14_phi = -1.0, B_A10 = -1.0;
611 double B_K11rho = -1.0, B_D_K11rho = -1.0, B_K14_1 = -1.0, B_K14_2 = -1.0, B_D_K14 = -1.0,
612 B_K14_km = -1.0, B_K14_kp = -1.0;
613 double B_K_Kst0pi0 = -1.0, B_D_KK = -1.0, B_K_Kstmpip = -1.0, B_D_K11_2 = -1.0,
614 B_K11_21 = -1.0, B_K11_22 = -1.0, B_A12 = -1.0;
615
616 PDF[0] = 0;
617 PDF[1] = 0;
618 double mass1sq, mass2sq, tt1, tt2, tt3, tt4;
619 double temp_PDF, tmp1, temp[2], tmp3, temp_PDF1;
620 double pro[2], pro0[2], pro1[2], pro2[2], pro3[2], pro4[2];
621 double t1A[4], t1D[4], t2D[4][4], B[3];
622 int isKstp = 0.0, isKst0 = 0.0, isKstm = 0.0, isRho = 0.0, isf0 = 0.0, isKmPip_S = 0.0,
623 isKpPi0_S = 0.0, isKmPi0_S = 0.0, isPiPi_S = 0.0, isA0980 = 0.0;
624 double proRho[2], proKstp[2], proKstm[2], proKst0[2], proPiPi_S[2], proKmPip_S[2],
625 proKpPi0_S[2], proKmPi0_S[2], proA0980[2], prof0[2];
626 for ( int i = 0; i < nstates; i++ )
627 {
628 temp_PDF = 0;
629 temp_PDF1 = 0;
630 cof[0] = amp[i] * cos( phase[i] );
631 cof[1] = amp[i] * sin( phase[i] );
632 mass1sq = mass1[i] * mass1[i];
633 mass2sq = mass2[i] * mass2[i];
634 if ( modetype[i] == 1 )
635 {
636 // phi-rhop[S/P/D]
637 if ( B_phi < 0.0 ) B_phi = barrier( 1, sphi, skaon1, skaon2, rRes2 );
638 if ( B_rho < 0.0 ) B_rho = barrier( 1, srhop, spionp, spi0, rRes2 );
639 if ( g2[i] == 0 )
640 {
641 for ( int i = 0; i < 4; i++ ) { temp_PDF += G[i][i] * t1phi[i] * t1rhop[i]; }
642 tmp1 = B_phi * B_rho * temp_PDF;
643 }
644 else if ( g2[i] == 1 )
645 {
646 calt1( pphi, prhop, t1D );
647 for ( int i = 0; i < 4; i++ )
648 {
649 tt1 = pD[i] * G[i][i];
650 for ( int j = 0; j < 4; j++ )
651 {
652 tt2 = t1D[j] * G[j][j];
653 for ( int k = 0; k < 4; k++ )
654 {
655 tt3 = t1phi[k] * G[k][k];
656 for ( int l = 0; l < 4; l++ )
657 { temp_PDF += E[i][j][k][l] * tt1 * tt2 * tt3 * t1rhop[l] * G[l][l]; }
658 }
659 }
660 }
661 if ( B_phirho1 < 0.0 )
662 {
663 B_phirho1 = barrier( 1, sD, sphi, srhop, rD2 );
664 tmp1 = B_phi * B_rho * B_phirho1 * temp_PDF;
665 }
666 }
667 else if ( g2[i] == 2 )
668 {
669 calt2( pphi, prhop, t2D );
670 for ( int i = 0; i < 4; i++ )
671 {
672 tt1 = t1phi[i] * G[i][i];
673 for ( int j = 0; j < 4; j++ ) { temp_PDF += t2D[i][j] * t1rhop[j] * G[j][j] * tt1; }
674 }
675 if ( B_phirho2 < 0.0 )
676 {
677 B_phirho2 = barrier( 2, sD, sphi, srhop, rD2 );
678 tmp1 = B_phi * B_rho * B_phirho2 * temp_PDF;
679 }
680 }
681 if ( isRho == 0 )
682 {
683 propagatorGS( mass2sq, mass2[i], width2[i], srhop, spionp, spi0, rRes2,
684 proRho ); // rho as mass2
685 isRho = 1.0;
686 }
687 if ( g0[i] == 1 )
688 { propagatorRBWl1( mass1sq, mass1[i], width1[i], sphi, skaon1, skaon2, rRes2, pro0 ); }
689 else if ( g0[i] == 0 )
690 {
691 pro0[0] = 1;
692 pro0[1] = 0;
693 }
694 if ( g1[i] == 1 )
695 {
696 pro1[0] = proRho[0];
697 pro1[1] = proRho[1];
698 }
699 else if ( g1[i] == 0 )
700 {
701 pro1[0] = 1;
702 pro1[1] = 0;
703 }
704 Com_Multi( pro0, pro1, pro );
705 amp_tmp[0] = tmp1 * pro[0];
706 amp_tmp[1] = tmp1 * pro[1];
707 // cout<<setprecision(20)<<"amp0:"<<amp_tmp[0]<<"amp1:"<<amp_tmp[1]<<endl;
708 }
709 else if ( modetype[i] == 2 )
710 {
711 if ( B_Kst0 < 0.0 ) B_Kst0 = barrier( 1, skstr0, skaon1, spionp, rRes2 );
712 if ( B_Kstp < 0.0 ) B_Kstp = barrier( 1, skstrp, skaon2, spi0, rRes2 );
713 if ( g2[i] == 0 )
714 {
715 for ( int i = 0; i < 4; i++ ) { temp_PDF += G[i][i] * t1kstr1[i] * t1kstr2[i]; }
716 tmp1 = B_Kst0 * B_Kstp * temp_PDF;
717 }
718 else if ( g2[i] == 1 )
719 {
720 calt1( pKstr0, pKstrp, t1D );
721 for ( int i = 0; i < 4; i++ )
722 {
723 tt1 = pD[i] * G[i][i];
724 for ( int j = 0; j < 4; j++ )
725 {
726 tt2 = t1D[j] * G[j][j];
727 for ( int k = 0; k < 4; k++ )
728 {
729 tt3 = t1kstr1[k] * G[k][k];
730 for ( int l = 0; l < 4; l++ )
731 { temp_PDF += E[i][j][k][l] * tt1 * tt2 * tt3 * t1kstr2[l] * G[l][l]; }
732 }
733 }
734 }
735 if ( B_KsKs1 < 0.0 )
736 {
737 B_KsKs1 = barrier( 1, sD, skstr0, skstrp, rD2 );
738 tmp1 = B_Kst0 * B_Kstp * B_KsKs1 * temp_PDF;
739 }
740 }
741 else if ( g2[i] == 2 )
742 {
743 calt2( pKstr0, pKstrp, t2D );
744 for ( int i = 0; i < 4; i++ )
745 {
746 tt1 = t1kstr1[i] * G[i][i];
747 for ( int j = 0; j < 4; j++ ) { temp_PDF += t2D[i][j] * t1kstr2[j] * G[j][j] * tt1; }
748 }
749 if ( B_KsKs2 < 0.0 )
750 {
751 B_KsKs2 = barrier( 2, sD, skstr0, skstrp, rD2 );
752 tmp1 = B_Kst0 * B_Kstp * B_KsKs2 * temp_PDF;
753 }
754 }
755 if ( isKst0 == 0 )
756 {
757 propagatorRBWl1( mass1sq, mass1[i], width1[i], skstr0, skaon1, spionp, rRes2,
758 proKst0 ); // Kst0 as mass1
759 isKst0 = 1;
760 }
761 if ( isKstp == 0 )
762 {
763 propagatorRBWl1( mass2sq, mass2[i], width2[i], skstrp, skaon2, spi0, rRes2,
764 proKstp ); // Kstp as mass2
765 isKstp = 1;
766 }
767 if ( g0[i] == 1 )
768 {
769 pro0[0] = proKst0[0];
770 pro0[1] = proKst0[1];
771 }
772 else if ( g0[i] == 0 )
773 {
774 pro0[0] = 1;
775 pro0[1] = 0;
776 }
777 if ( g1[i] == 1 )
778 {
779 pro1[0] = proKstp[0];
780 pro1[1] = proKstp[1];
781 }
782 else if ( g1[i] == 0 )
783 {
784 pro1[0] = 1;
785 pro1[1] = 0;
786 }
787 Com_Multi( pro0, pro1, pro );
788 amp_tmp[0] = tmp1 * pro[0];
789 amp_tmp[1] = tmp1 * pro[1];
790 // cout<<setprecision(20)<<"amp0:"<<amp_tmp[0]<<"amp1:"<<amp_tmp[1]<<endl;
791 }
792 else if ( modetype[i] == 3 )
793 {
794 // K1400 K1270.((K-Pi+)_VPi0)A[S/D]K+ && ((K-Pi0)_VPi+)A[S/D]K+ //D->AP(1) A->VP(0,2)
795 // V->PP(1)
796 if ( B_DA1 < 0.0 ) B_DA1 = barrier( 1, sD, sk11, skaon2, rD2 );
797 if ( B_Kst0 < 0.0 ) B_Kst0 = barrier( 1, skstr0, skaon1, spionp, rRes2 );
798 calt1( pK11, Kp, t1D );
799 if ( g2[i] == 0 )
800 {
801 for ( int i = 0; i < 4; i++ )
802 {
803 tt1 = t1D[i] * G[i][i];
804 for ( int j = 0; j < 4; j++ )
805 { temp_PDF += tt1 * ( G[i][j] - pK11[i] * pK11[j] / sk11 ) * t1kstr1[j] * G[j][j]; }
806 }
807 tmp1 = B_Kst0 * B_DA1 * temp_PDF;
808 }
809 else if ( g2[i] == 2 )
810 {
811 for ( int i = 0; i < 4; i++ )
812 {
813 tt1 = t1D[i] * G[i][i];
814 for ( int j = 0; j < 4; j++ )
815 { temp_PDF += tt1 * t2k11[i][j] * t1kstr1[j] * G[j][j]; }
816 }
817 if ( B_A2 < 0.0 ) B_A2 = barrier( 2, sk11, skstr0, spi0, rRes2 );
818 tmp1 = B_Kst0 * B_A2 * B_DA1 * temp_PDF;
819 }
820 if ( isKst0 == 0 )
821 {
822 propagatorRBWl1( mass1sq, mass1[i], width1[i], skstr0, skaon1, spionp, rRes2,
823 proKst0 ); // Kst0 also as mass1
824 isKst0 = 1;
825 }
826 if ( g0[i] == 1 )
827 {
828 pro0[0] = proKst0[0];
829 pro0[1] = proKst0[1];
830 }
831 else if ( g0[i] == 0 )
832 {
833 pro0[0] = 1;
834 pro0[1] = 0;
835 }
836 if ( g1[i] == 1 )
837 {
838 propagatorRBW( mass2sq, mass2[i], width2[i], sk11, skstr0, spi0, rRes2, g2[i],
839 pro1 ); // K1270(K*) as mass2
840 }
841 else if ( g1[i] == 0 )
842 {
843 pro1[0] = 1;
844 pro1[1] = 0;
845 }
846 Com_Multi( pro0, pro1, pro );
847 amp_tmp[0] = tmp1 * pro[0];
848 amp_tmp[1] = tmp1 * pro[1];
849 // cout<<setprecision(20)<<"amp0:"<<amp_tmp[0]<<"amp1:"<<amp_tmp[1]<<endl;
850 }
851 else if ( modetype[i] == 333 )
852 {
853 // K1400 K1270..((K-Pi+)_VPi0)A[S/D]K+ && ((K-Pi0)_VPi+)A[S/D]K+ //D->AP(1) A->VP(0,2)
854 // V->PP(1)
855 if ( B_DA1 < 0.0 ) B_DA1 = barrier( 1, sD, sk11, skaon2, rD2 );
856 if ( B_Kstm < 0.0 ) B_Kstm = barrier( 1, skstrm, skaon1, spi0, rRes2 );
857 calt1( pK11, Kp, t1D );
858 if ( g2[i] == 0 )
859 {
860 for ( int i = 0; i < 4; i++ )
861 {
862 tt1 = t1D[i] * G[i][i];
863 for ( int j = 0; j < 4; j++ )
864 { temp_PDF1 += tt1 * ( G[i][j] - pK11[i] * pK11[j] / sk11 ) * t1kstr3[j] * G[j][j]; }
865 }
866 tmp3 = B_Kstm * B_DA1 * temp_PDF1;
867 }
868 else if ( g2[i] == 2 )
869 {
870 for ( int i = 0; i < 4; i++ )
871 {
872 tt1 = t1D[i] * G[i][i];
873 for ( int j = 0; j < 4; j++ )
874 { temp_PDF1 += tt1 * t2k12[i][j] * t1kstr3[j] * G[j][j]; }
875 }
876 if ( B_A1 < 0.0 ) B_A1 = barrier( 2, sk11, skstrm, spionp, rRes2 );
877 tmp3 = B_Kstm * B_A1 * B_DA1 * temp_PDF1;
878 }
879 if ( g0[i] == 1 )
880 { propagatorRBWl1( mass1sq, mass1[i], width1[i], skstrm, skaon1, spi0, rRes2, pro2 ); }
881 else if ( g0[i] == 0 )
882 {
883 pro2[0] = 1;
884 pro2[1] = 0;
885 }
886 if ( g1[i] == 1 )
887 {
888 propagatorRBW( mass2sq, mass2[i], width2[i], sk11, skstrm, spionp, rRes2, g2[i],
889 pro3 );
890 }
891 else if ( g1[i] == 0 )
892 {
893 pro3[0] = 1;
894 pro3[1] = 0;
895 }
896 Com_Multi( pro2, pro3, pro4 );
897 amp_tmp[0] = -1.414 * tmp3 * pro4[0];
898 amp_tmp[1] = -1.414 * tmp3 * pro4[1];
899 // cout<<setprecision(20)<<"amp0:"<<amp_tmp[0]<<"amp1:"<<amp_tmp[1]<<endl;
900 }
901 else if ( modetype[i] == 4 )
902 {
903 // K1270..(K-Rho+)_A[S/D]K+ //D->AP(1) A->VP(0,2) V->PP(1)
904 if ( B_rho < 0.0 ) B_rho = barrier( 1, srhop, spionp, spi0, rRes2 );
905 if ( B_D_Arho < 0.0 ) B_D_Arho = barrier( 1, sD, sk11, skaon2, rD2 );
906 calt1( pK11, Kp, t1D );
907 if ( g2[i] == 0 )
908 {
909 for ( int i = 0; i < 4; i++ )
910 {
911 tt1 = t1D[i] * G[i][i];
912 for ( int j = 0; j < 4; j++ )
913 { temp_PDF += tt1 * ( G[i][j] - pK11[i] * pK11[j] / sk11 ) * t1rhop[j] * G[j][j]; }
914 }
915 tmp1 = B_rho * B_D_Arho * temp_PDF;
916 }
917 else if ( g2[i] == 2 )
918 {
919 for ( int i = 0; i < 4; i++ )
920 {
921 tt1 = t1D[i] * G[i][i];
922 for ( int j = 0; j < 4; j++ )
923 { temp_PDF += tt1 * t2k21[i][j] * t1rhop[j] * G[j][j]; }
924 }
925 if ( B_Arho < 0.0 )
926 {
927 B_Arho = barrier( 2, sk11, srhop, skaon1, rRes2 );
928 tmp1 = B_rho * B_Arho * B_D_Arho * temp_PDF;
929 }
930 }
931 if ( isRho == 0 )
932 {
933 propagatorGS( mass2sq, mass2[i], width2[i], srhop, spionp, spi0, rRes2,
934 proRho ); // rho also as mass2
935 isRho = 1.0;
936 }
937 if ( g0[i] == 1 )
938 propagatorRBW( mass1sq, mass1[i], width1[i], sk11, srhop, skaon1, rRes2, g2[i],
939 pro0 ); // K1270(rho) as mass1
940 else if ( g0[i] == 0 )
941 {
942 pro0[0] = 1;
943 pro0[1] = 0;
944 }
945 if ( g1[i] == 1 )
946 {
947 pro1[0] = proRho[0];
948 pro1[1] = proRho[1];
949 }
950 else if ( g1[i] == 0 )
951 {
952 pro1[0] = 1;
953 pro1[1] = 0;
954 }
955 Com_Multi( pro0, pro1, pro );
956 amp_tmp[0] = tmp1 * pro[0];
957 amp_tmp[1] = tmp1 * pro[1];
958 // cout<<setprecision(20)<<"amp0:"<<amp_tmp[0]<<"amp1:"<<amp_tmp[1]<<endl;
959 }
960 else if ( modetype[i] == 30 )
961 {
962 // eta1475_(a0980_pi0)Pi+ //D->PP(0) P->SP(0) S->PP(0)
963 double sKmm[2] = { skaon1, mass_Pion_N * mass_Pion_N };
964 double sKpp[2] = { skaon2, mass_Eta * mass_Eta };
965 if ( isA0980 == 0 )
966 {
967 propagatora0980( mass1sq, mass1[i], sphi, sKmm, sKpp, proA0980 );
968 isA0980 = 1.0;
969 }
970 if ( g0[i] == 1 )
971 {
972 pro0[0] = proA0980[0];
973 pro0[1] = proA0980[1];
974 }
975 else if ( g0[i] == 0 )
976 {
977 pro0[0] = 1;
978 pro0[1] = 0;
979 }
980 if ( g1[i] == 1 )
981 propagatorRBW( mass2sq, mass2[i], width2[i], sk14, sphi, spi0, rRes2, 0, pro1 );
982 else if ( g1[i] == 0 )
983 {
984 pro1[0] = 1;
985 pro1[1] = 0;
986 }
987 Com_Multi( pro0, pro1, pro );
988 amp_tmp[0] = pro[0];
989 amp_tmp[1] = pro[1];
990 // cout<<setprecision(20)<<"amp0:"<<amp_tmp[0]<<"amp1:"<<amp_tmp[1]<<endl;
991 }
992 else if ( modetype[i] == 31 )
993 {
994 // f1(1420)_(a0980_pi0)Pi+ //D->AP(1) A->SP(1) S->PP(0)
995 if ( B_K14_phi < 0.0 ) B_K14_phi = barrier( 1, sk14, sphi, spi0, rRes2 );
996 if ( B_D_K14 < 0.0 ) B_D_K14 = barrier( 1, sD, sk14, spionp, rD2 );
997 double sKmm[2] = { skaon1, mass_Pion_N * mass_Pion_N };
998 double sKpp[2] = { skaon2, mass_Eta * mass_Eta };
999 calt1( pK14, Pip, t1D );
1000 calt1( pphi, Pi0, t1A );
1001 for ( int w = 0; w < 4; w++ ) { temp_PDF += t1D[w] * t1A[w] * G[w][w]; }
1002 tmp1 = temp_PDF * B_K14_phi * B_D_K14;
1003 if ( isA0980 == 0 )
1004 {
1005 propagatora0980( mass1sq, mass1[i], sphi, sKmm, sKpp, proA0980 );
1006 isA0980 = 1.0;
1007 }
1008 if ( g0[i] == 1 )
1009 {
1010 pro0[0] = proA0980[0];
1011 pro0[1] = proA0980[1];
1012 }
1013 else if ( g0[i] == 0 )
1014 {
1015 pro0[0] = 1;
1016 pro0[1] = 0;
1017 }
1018 if ( g1[i] == 1 )
1019 propagatorRBWl1( mass2sq, mass2[i], width2[i], sk14, sphi, spi0, rRes2, pro1 );
1020 else if ( g1[i] == 0 )
1021 {
1022 pro1[0] = 1;
1023 pro1[1] = 0;
1024 }
1025 Com_Multi( pro0, pro1, pro );
1026 amp_tmp[0] = tmp1 * pro[0];
1027 amp_tmp[1] = tmp1 * pro[1];
1028 // cout<<setprecision(20)<<"amp0:"<<amp_tmp[0]<<"amp1:"<<amp_tmp[1]<<endl;
1029 }
1030 else if ( modetype[i] == 72 )
1031 {
1032 // eta1450_[(K-pi0)s K+ K-] Pi+ //D->PP(0) P->SP(0) S->PP(0)
1033 if ( isKmPi0_S == 0 )
1034 {
1035 KPiSLASS( skstrm, skaon1, spi0, proKmPi0_S );
1036 isKmPi0_S = 1;
1037 }
1038 if ( g0[i] == 1 )
1039 {
1040 pro0[0] = proKmPi0_S[0];
1041 pro0[1] = proKmPi0_S[1];
1042 }
1043 else if ( g0[i] == 0 )
1044 {
1045 pro0[0] = 1;
1046 pro0[1] = 0;
1047 }
1048 if ( g1[i] == 1 )
1049 { propagatorRBW( mass2sq, mass2[i], width2[i], sk14, skstrm, skaon2, rRes2, 0, pro1 ); }
1050 else if ( g1[i] == 0 )
1051 {
1052 pro1[0] = 1;
1053 pro1[1] = 0;
1054 }
1055 Com_Multi( pro0, pro1, pro );
1056 amp_tmp[0] = pro[0];
1057 amp_tmp[1] = pro[1];
1058 // cout<<setprecision(20)<<"amp0:"<<amp_tmp[0]<<"amp1:"<<amp_tmp[1]<<endl;
1059 }
1060 else if ( modetype[i] == 721 )
1061 {
1062 // eta1450_[(K+pi0)s K-] Pi+ //D->PP(0) P->SP(0) S->PP(0)
1063 if ( isKpPi0_S == 0 )
1064 {
1065 KPiSLASS( skstrp, skaon2, spi0, proKpPi0_S );
1066 isKpPi0_S = 1;
1067 }
1068 if ( g0[i] == 1 )
1069 {
1070 pro2[0] = proKpPi0_S[0];
1071 pro2[1] = proKpPi0_S[1];
1072 }
1073 else if ( g0[i] == 0 )
1074 {
1075 pro2[0] = 1;
1076 pro2[1] = 0;
1077 }
1078 if ( g1[i] == 1 )
1079 { propagatorRBW( mass2sq, mass2[i], width2[i], sk14, skstrp, skaon1, rRes2, 0, pro3 ); }
1080 else if ( g1[i] == 0 )
1081 {
1082 pro3[0] = 1;
1083 pro3[1] = 0;
1084 }
1085 Com_Multi( pro2, pro3, pro4 );
1086 amp_tmp[0] = pro4[0];
1087 amp_tmp[1] = pro4[1];
1088 // cout<<setprecision(20)<<"amp0:"<<amp_tmp[0]<<"amp1:"<<amp_tmp[1]<<endl;
1089 }
1090 else if ( modetype[i] == 73 )
1091 {
1092 // f1510..((K-Pi0)_VK+)A[S/D]Pi+ //D->AP(1) A->VP(0,2) V->PP(1)
1093 if ( B_DA2 < 0.0 ) B_DA2 = barrier( 1, sD, sk10, spionp, rD2 );
1094 if ( B_Kstm < 0.0 ) B_Kstm = barrier( 1, skstrm, skaon1, spi0, rRes2 );
1095 calt1( pK10, Pip, t1D );
1096 if ( g2[i] == 0 )
1097 {
1098 for ( int i = 0; i < 4; i++ )
1099 {
1100 tt1 = t1D[i] * G[i][i];
1101 for ( int j = 0; j < 4; j++ )
1102 { temp_PDF += tt1 * ( G[i][j] - pK10[i] * pK10[j] / sk10 ) * t1kstr3[j] * G[j][j]; }
1103 }
1104 tmp1 = B_Kstm * B_DA2 * temp_PDF;
1105 }
1106 else if ( g2[i] == 2 )
1107 {
1108 for ( int i = 0; i < 4; i++ )
1109 {
1110 tt1 = t1D[i] * G[i][i];
1111 for ( int j = 0; j < 4; j++ )
1112 { temp_PDF += tt1 * t2k13[i][j] * t1kstr3[j] * G[j][j]; }
1113 }
1114 if ( B_A10 < 0.0 ) B_A10 = barrier( 2, sk10, skstrm, skaon2, rRes2 );
1115 tmp1 = B_Kstm * B_A10 * B_DA2 * temp_PDF;
1116 }
1117 if ( isKstm == 0 )
1118 {
1119 propagatorRBWl1( mass1sq, mass1[i], width1[i], skstrm, skaon1, spi0, rRes2, proKstm );
1120 isKstm = 1;
1121 }
1122 if ( g0[i] == 1 )
1123 {
1124 pro0[0] = proKstm[0];
1125 pro0[1] = proKstm[1];
1126 }
1127 else if ( g0[i] == 0 )
1128 {
1129 pro0[0] = 1;
1130 pro0[1] = 0;
1131 }
1132 if ( g1[i] == 1 )
1133 {
1134 propagatorRBW( mass2sq, mass2[i], width2[i], sk10, skstrm, skaon2, rRes2, g2[i],
1135 pro1 );
1136 }
1137 else if ( g1[i] == 0 )
1138 {
1139 pro1[0] = 1;
1140 pro1[1] = 0;
1141 }
1142 Com_Multi( pro0, pro1, pro );
1143 amp_tmp[0] = tmp1 * pro[0];
1144 amp_tmp[1] = tmp1 * pro[1];
1145 // cout<<setprecision(20)<<"amp0:"<<amp_tmp[0]<<"amp1:"<<amp_tmp[1]<<endl;
1146 }
1147 else if ( modetype[i] == 731 )
1148 {
1149 // f1510..( ((K+Pi0)_VK-)A[S/D]Pi+ //D->AP(1) A->VP(0,2) V->PP(1)
1150 if ( B_DA2 < 0.0 ) B_DA2 = barrier( 1, sD, sk10, spionp, rD2 );
1151 if ( B_Kstp < 0.0 ) B_Kstp = barrier( 1, skstrp, skaon2, spi0, rRes2 );
1152 calt1( pK10, Pip, t1D );
1153 if ( g2[i] == 0 )
1154 {
1155 for ( int i = 0; i < 4; i++ )
1156 {
1157 tt1 = t1D[i] * G[i][i];
1158 for ( int j = 0; j < 4; j++ )
1159 { temp_PDF1 += tt1 * ( G[i][j] - pK10[i] * pK10[j] / sk10 ) * t1kstr4[j] * G[j][j]; }
1160 }
1161 tmp3 = B_Kstp * B_DA2 * temp_PDF1;
1162 }
1163 else if ( g2[i] == 2 )
1164 {
1165 for ( int i = 0; i < 4; i++ )
1166 {
1167 tt1 = t1D[i] * G[i][i];
1168 for ( int j = 0; j < 4; j++ )
1169 { temp_PDF1 += tt1 * t2k14[i][j] * t1kstr4[j] * G[j][j]; }
1170 }
1171 if ( B_A12 < 0.0 ) B_A12 = barrier( 2, sk12, skstrp, skaon1, rRes2 );
1172 tmp3 = B_Kstp * B_A12 * B_DA2 * temp_PDF1;
1173 }
1174 if ( isKstp == 0 )
1175 {
1176 propagatorRBWl1( mass1sq, mass1[i], width1[i], skstrp, skaon2, spi0, rRes2, proKstp );
1177 isKstp = 1;
1178 }
1179 if ( g0[i] == 1 )
1180 {
1181 pro2[0] = proKstp[0];
1182 pro2[1] = proKstp[1];
1183 }
1184 else if ( g0[i] == 0 )
1185 {
1186 pro2[0] = 1;
1187 pro2[1] = 0;
1188 }
1189 if ( g1[i] == 1 )
1190 {
1191 propagatorRBW( mass2sq, mass2[i], width2[i], sk12, skstrp, skaon1, rRes2, g2[i],
1192 pro3 );
1193 }
1194 else if ( g1[i] == 0 )
1195 {
1196 pro3[0] = 1;
1197 pro3[1] = 0;
1198 }
1199 Com_Multi( pro2, pro3, pro4 );
1200 amp_tmp[0] = -tmp3 * pro4[0];
1201 amp_tmp[1] = -tmp3 * pro4[1];
1202 // cout<<setprecision(20)<<"amp0:"<<amp_tmp[0]<<"amp1:"<<amp_tmp[1]<<endl;
1203 }
1204 else if ( modetype[i] == 50 )
1205 {
1206 // a0980 rho+ //D->SV(1) S->PP(0) V->PP(1)
1207 if ( B_rho < 0.0 ) B_rho = barrier( 1, srhop, spionp, spi0, rRes2 );
1208 if ( B_phirho1 < 0.0 ) B_phirho1 = barrier( 1, sD, sphi, srhop, rD2 );
1209 calt1( pphi, prhop, t1D );
1210 for ( int w = 0; w < 4; w++ ) { temp_PDF += t1D[w] * t1rhop[w] * G[w][w]; }
1211 tmp1 = temp_PDF * B_rho * B_phirho1;
1212 double sKmm[2] = { skaon1, mass_Pion_N * mass_Pion_N };
1213 double sKpp[2] = { skaon2, mass_Eta * mass_Eta };
1214 if ( isA0980 == 0 )
1215 {
1216 propagatora0980( mass1sq, mass1[i], sphi, sKmm, sKpp, proA0980 );
1217 isA0980 = 1.0;
1218 }
1219 if ( isRho == 0 )
1220 {
1221 propagatorGS( mass2sq, mass2[i], width2[i], srhop, spionp, spi0, rRes2,
1222 proRho ); // rho as mass2
1223 isRho = 1.0;
1224 }
1225 if ( g0[i] == 1 )
1226 {
1227 pro0[0] = proA0980[0];
1228 pro0[1] = proA0980[1];
1229 }
1230 else if ( g0[i] == 0 )
1231 {
1232 pro0[0] = 1;
1233 pro0[1] = 0;
1234 }
1235 if ( g1[i] == 1 )
1236 {
1237 pro1[0] = proRho[0];
1238 pro1[1] = proRho[1];
1239 }
1240 else if ( g1[i] == 0 )
1241 {
1242 pro1[0] = 1;
1243 pro1[1] = 0;
1244 }
1245 Com_Multi( pro0, pro1, pro );
1246 amp_tmp[0] = tmp1 * pro[0];
1247 amp_tmp[1] = tmp1 * pro[1];
1248 // cout<<setprecision(20)<<"amp0:"<<amp_tmp[0]<<"amp1:"<<amp_tmp[1]<<endl;
1249 }
1250 Com_Multi( amp_tmp, cof, amp_PDF );
1251 PDF[0] += amp_PDF[0];
1252 PDF[1] += amp_PDF[1];
1253 }
1254 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
1255 if ( value <= 0 ) value = 1e-20;
1256 Result = value;
1257 // cout<<"!!!!!!!!!!!!"<<setprecision(20)<<value<<endl;
1258}
mass
double mass
Definition CosmicGenerator.cxx:128
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
g1
TF1 * g1
Definition DataBase/tau_mode.c:58
w
double w
Definition EvtPhokharaDef.hh:28
unit
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per saves r n generator level $ !Flag for chat level in !Latex Output unit
Definition FoamA.h:90
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
q
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition KKsem.h:33
CS
TCrossPart * CS
Definition Mcgpj.cxx:51
EvtDecayBase::checkSpinDaughter
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
Definition EvtDecayBase.cc:492
EvtDecayBase::checkSpinParent
void checkSpinParent(EvtSpinType::spintype sp)
Definition EvtDecayBase.cc:479
EvtDecayBase::setProbMax
void setProbMax(double prbmx)
Definition EvtDecayBase.cc:261
EvtDecayBase::checkNDaug
void checkNDaug(int d1, int d2=-1)
Definition EvtDecayBase.cc:465
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition EvtDecayBase.hh:68
EvtDecayBase::checkNArg
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
Definition EvtDecayBase.cc:449
EvtDecayProb::setProb
void setProb(double prob)
Definition EvtDecayProb.hh:40
EvtDsToKKpipi0::getName
void getName(std::string &name)
Definition EvtDsToKKpipi0.cc:36
EvtDsToKKpipi0::~EvtDsToKKpipi0
virtual ~EvtDsToKKpipi0()
Definition EvtDsToKKpipi0.cc:34
EvtDsToKKpipi0::decay
void decay(EvtParticle *p)
Definition EvtDsToKKpipi0.cc:251
EvtDsToKKpipi0::clone
EvtDecayBase * clone()
Definition EvtDsToKKpipi0.cc:38
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition EvtParticle.cc:108
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition EvtParticle.cc:80
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition EvtParticle.cc:1042
EvtVector4R::get
double get(int i) const
Definition EvtVector4R.hh:171
t
int t()
Definition t.c:1