BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToKpPipPimPi0.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: EvtDsToKpPipPimPi0.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 "EvtDsToKpPipPimPi0.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 "TComplex.h"
29#include "TMatrix.h"
30#include "TMatrixD.h"
31#include <fstream>
32#include <stdlib.h>
33#include <string>
34
35using std::endl;
36
37EvtDsToKpPipPimPi0::~EvtDsToKpPipPimPi0() {}
38
39void EvtDsToKpPipPimPi0::getName( std::string& model_name ) { model_name = "DsToKpPipPimPi0"; }
40
41EvtDecayBase* EvtDsToKpPipPimPi0::clone() { return new EvtDsToKpPipPimPi0; }
42
43void EvtDsToKpPipPimPi0::init() {
44 // check that there are 0 arguments
45 checkNArg( 0 );
46 checkNDaug( 4 );
47
48 checkSpinParent( EvtSpinType::SCALAR );
49 checkSpinDaughter( 0, EvtSpinType::SCALAR );
50 checkSpinDaughter( 1, EvtSpinType::SCALAR );
51 checkSpinDaughter( 2, EvtSpinType::SCALAR );
52 checkSpinDaughter( 3, EvtSpinType::SCALAR );
53
54 int mode = 0;
55 phi[1] = -4.1936e+00;
56 phi[2] = 2.4056e+00;
57 phi[3] = 1.8102e+00;
58 phi[4] = -1.6106e+00;
59 phi[5] = -1.6106e+00;
60 phi[6] = 5.1214e+00;
61 phi[7] = 5.1214e+00;
62 phi[8] = 6.0675e-01;
63 phi[9] = -5.2450e+00;
64 phi[10] = -2.8735e+00;
65
66 rho[1] = 1.0981e+00;
67 rho[2] = 4.6322e-01;
68 rho[3] = 1.3138e+00;
69 rho[4] = 1.2709e+00;
70 rho[5] = 1.2709e+00;
71 rho[6] = 1.2953e+00;
72 rho[7] = 1.2953e+00;
73 rho[8] = 2.3583e+00;
74 rho[9] = 7.8907e+00;
75 rho[10] = 6.5982e-01;
76
77 phi[0] = 0;
78 rho[0] = 1;
79 mDsM = 1.9683;
80
81 //------------------------new----------------------
82 mKst0 = 0.89555;
83 mKstp = 0.89166;
84 mrho = 0.7751;
85 mrho0 = 0.7753;
86 mK1270 = 1.289;
87 mK1400 = 1.403;
88 mK1650 = 1.672;
89 mOmega = 0.78266;
90 mA1 = 1.230;
91
92 GKst0 = 0.0473;
93 GKstp = 0.0508;
94 Grho = 0.1491;
95 Grho0 = 0.1478;
96 GK1270 = 0.116;
97 GK1400 = 0.174;
98 GK1650 = 0.158;
99 GOmega = 0.00868;
100 GA1 = 0.42;
101
102 // cout << "DsToKpPipPimPi0 :" << endl;
103 // for (int i=0; i<10; i++) {
104 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
105 // }
106
107 mass_Pion = 0.13957;
108 mass_Pion_N = 0.134977;
109 mass_Eta = 0.547862;
110 mass_Kaon = 0.493677;
111 math_pi = 3.1415926;
112
113 rD2 = 5.0; // 5*5
114 rRes1 = 3.0; // 3*3
115 rRes2 = 9.0; // 3*3
116
117 GS1 = 0.636619783;
118 GS2 = 0.01860182466;
119 GS3 = 0.1591549458; // 1/(2*math_2pi)
120 GS4 = 0.00620060822; // mass_Pion2/math_pi
121
122 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
123 int EE[4][4][4][4] = {
124 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
125 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, -1, 0 } },
126 { { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 1, 0, 0 } },
127 { { 0, 0, 0, 0 }, { 0, 0, 1, 0 }, { 0, -1, 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, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
130 { { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 } },
131 { { 0, 0, -1, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 } } },
132 { { { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, -1, 0, 0 } },
133 { { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 } },
134 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
135 { { 0, 1, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } },
136 { { { 0, 0, 0, 0 }, { 0, 0, -1, 0 }, { 0, 1, 0, 0 }, { 0, 0, 0, 0 } },
137 { { 0, 0, 1, 0 }, { 0, 0, 0, 0 }, { -1, 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, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } } };
140 for ( int i = 0; i < 4; i++ )
141 {
142 for ( int j = 0; j < 4; j++ )
143 {
144 G[i][j] = GG[i][j];
145 for ( int k = 0; k < 4; k++ )
146 {
147 for ( int l = 0; l < 4; l++ ) { E[i][j][k][l] = EE[i][j][k][l]; }
148 }
149 }
150 }
151}
152
153void EvtDsToKpPipPimPi0::initProbMax() { setProbMax( 130000 ); }
154
155void EvtDsToKpPipPimPi0::decay( EvtParticle* p ) {
156 //--------------max
157
158 /* double maxprob=0,maxq1270;
159 for(int ir=0;ir<=60000000;ir++){
160 p->initializePhaseSpace(getNDaug(),getDaugs());
161 EvtVector4R _kp = p->getDaug(0)->getP4();
162 EvtVector4R _pip = p->getDaug(1)->getP4();
163 EvtVector4R _pim = p->getDaug(2)->getP4();
164 EvtVector4R _pi0 = p->getDaug(3)->getP4();
165
166 double _Pip[4],_Pim[4],_Pi0[4],_Kp[4];
167 _Kp[0] = _kp.get(0); _Pip[0] = _pip.get(0); _Pim[0] = _pim.get(0); _Pi0[0] =
168 _pi0.get(0); _Kp[1] = _kp.get(1); _Pip[1] = _pip.get(1); _Pim[1] = _pim.get(1); _Pi0[1] =
169 _pi0.get(1); _Kp[2] = _kp.get(2); _Pip[2] = _pip.get(2); _Pim[2] = _pim.get(2); _Pi0[2] =
170 _pi0.get(2); _Kp[3] = _kp.get(3); _Pip[3] = _pip.get(3); _Pim[3] = _pim.get(3); _Pi0[3] =
171 _pi0.get(3);
172
173 double _prob;
174 double value,q1270;
175 int nstates=11;
176 int modetype[11]= {1,1,2,10,12,11,16,17,14,2,15};
177 int g0[11] = {1,1,1,1,1,1,1,1,1,0,0};
178 int g1[11] = {1,1,1,1,1,1,1,1,1,1,0};
179 int g2[11] = {0,1,1,0,0,0,0,0,1,0,0};
180 double mass1[11] = {mKst0,mKst0,mKstp,mK1270,mKstp,mKst0,mA1,mA1,mOmega,mKst0,mKst0};
181 double mass2[11] = {mrho,mrho,mrho0,mrho,mK1400,mK1400,mrho,mrho,mrho,mrho,mrho};
182 double width1[11] =
183 {GKst0,GKst0,GKstp,GK1270,GKstp,GKst0,GA1,GA1,GOmega,GKst0,GKst0}; double width2[11] =
184 {Grho,Grho,Grho0,Grho,GK1400,GK1400,Grho,Grho,Grho,mrho,mrho};
185
186 //
187 calEvaMy(_Kp,_Pip,_Pim,_Pi0,mass1,mass2,width1,width2,rho,phi,g0,g1,g2,modetype,nstates,_prob);
188 calEvaMy(_Kp,_Pip,_Pim,_Pi0,mass1,mass2,width1,width2,rho,phi,g0,g1,g2,modetype,nstates,_prob,q1270);
189
190 if(_prob>maxprob) {
191 maxprob=_prob;
192 maxq1270=q1270;
193 printf("ir = %d,maxprob = %.10f,prob = %.10f,q1270 =
194 %.10f\n\n",ir,maxprob,_prob,maxq1270);
195 // cout<<"mm: "<<_Kp[0]<<" "<<_Pip[0]<<" "<<_Pim[0]<<" "<<_Pi0[0]<<endl;
196 // cout<<"mm: "<<_Kp[1]<<" "<<_Pip[1]<<" "<<_Pim[1]<<" "<<_Pi0[1]<<endl;
197 // cout<<"mm: "<<_Kp[2]<<" "<<_Pip[2]<<" "<<_Pim[2]<<" "<<_Pi0[2]<<endl;
198 // cout<<"mm: "<<_Kp[3]<<" "<<_Pip[3]<<" "<<_Pim[3]<<" "<<_Pi0[3]<<endl;
199 }
200 }
201 printf("maxprob = %.10f\n", maxprob);
202 */
203
204 p->initializePhaseSpace( getNDaug(), getDaugs() );
205 EvtVector4R kp = p->getDaug( 0 )->getP4();
206 EvtVector4R pip = p->getDaug( 1 )->getP4();
207 EvtVector4R pim = p->getDaug( 2 )->getP4();
208 EvtVector4R pi0 = p->getDaug( 3 )->getP4();
209
210 double Pip[4], Pim[4], Kp[4], Pi0[4];
211
212 //----------------------------------------------------------------------------------------
213 Kp[0] = kp.get( 0 );
214 Pip[0] = pip.get( 0 );
215 Pim[0] = pim.get( 0 );
216 Pi0[0] = pi0.get( 0 );
217 Kp[1] = kp.get( 1 );
218 Pip[1] = pip.get( 1 );
219 Pim[1] = pim.get( 1 );
220 Pi0[1] = pi0.get( 1 );
221 Kp[2] = kp.get( 2 );
222 Pip[2] = pip.get( 2 );
223 Pim[2] = pim.get( 2 );
224 Pi0[2] = pi0.get( 2 );
225 Kp[3] = kp.get( 3 );
226 Pip[3] = pip.get( 3 );
227 Pim[3] = pim.get( 3 );
228 Pi0[3] = pi0.get( 3 );
229 //---------------------------------------------------------------------------------------
230 double value, q1270;
231 int nstates = 11;
232 int modetype[11] = { 1, 1, 2, 10, 12, 11, 16, 17, 14, 2, 15 };
233 int g0[11] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0 };
234 int g1[11] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0 };
235 int g2[11] = { 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0 };
236 double mass1[11] = { mKst0, mKst0, mKstp, mK1270, mKstp, mKst0,
237 mA1, mA1, mOmega, mKst0, mKst0 };
238 double mass2[11] = { mrho, mrho, mrho0, mrho, mK1400, mK1400, mrho, mrho, mrho, mrho, mrho };
239 double width1[11] = { GKst0, GKst0, GKstp, GK1270, GKstp, GKst0,
240 GA1, GA1, GOmega, GKst0, GKst0 };
241 double width2[11] = { Grho, Grho, Grho0, Grho, GK1400, GK1400,
242 Grho, Grho, Grho, mrho, mrho };
243
244 calEvaMy( Kp, Pip, Pim, Pi0, mass1, mass2, width1, width2, rho, phi, g0, g1, g2, modetype,
245 nstates, value, q1270 );
246 setProb( value );
247 return;
248}
249
250double EvtDsToKpPipPimPi0::CalRho4pi( double s ) {
251 if ( s >= 1. ) { return sqrt( ( s - 16. * mass_Pion * mass_Pion ) / s ); }
252 else
253 {
254 double gam = 0.0005;
255 double s2 = s * s;
256 double s3 = s2 * s;
257 double s4 = s2 * s2;
258 double s5 = s4 * s;
259 double s6 = s5 * s;
260
261 gam -= 0.0193 * s;
262 gam += 0.1385 * s2;
263 gam -= 0.2084 * s3;
264 gam -= 0.2974 * s4;
265 gam += 0.1366 * s5;
266 gam += 1.0789 * s6;
267
268 return gam;
269 }
270}
271TComplex EvtDsToKpPipPimPi0::ResonanceSkm( double& m2 ) {
272 double g11 = 0.22889, g12 = -0.55377, g13 = 0.00, g14 = -0.39899, g15 = -0.34639;
273 double g21 = 0.94128, g22 = 0.55095, g23 = 0.00, g24 = 0.39065, g25 = 0.31503;
274 double g31 = 0.36856, g32 = 0.23888, g33 = 0.55639, g34 = 0.18340, g35 = 0.18681;
275 double g41 = 0.33650, g42 = 0.40907, g43 = 0.85679, g44 = 0.19906, g45 = -0.00984;
276 double g51 = 0.18171, g52 = -0.17558, g53 = -0.79658, g54 = -0.00355, g55 = 0.22358;
277 double fs11 = 0.23399, fs12 = 0.15044, fs13 = -0.20545, fs14 = 0.32825, fs15 = 0.35412;
278 double m_1 = 0.65100, m_2 = 1.20360, m_3 = 1.55817, m_4 = 1.21000, m_5 = 1.82206;
279 double ss0 = -3.92637;
280 double sp0 = -3.0; // v0
281 double sA0 = -0.15;
282 double sA = 1.0;
283 const double mass_Eta = 0.547862;
284 const double mass_Etap = 0.95778;
285
286 double km11 = ( g11 * g11 / ( m_1 * m_1 - m2 ) + g21 * g21 / ( m_2 * m_2 - m2 ) +
287 g31 * g31 / ( m_3 * m_3 - m2 ) + g41 * g41 / ( m_4 * m_4 - m2 ) +
288 g51 * g51 / ( m_5 * m_5 - m2 ) + fs11 * ( 1 - ss0 ) / ( m2 - ss0 ) ) *
289 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
290 double km12 = ( g11 * g12 / ( m_1 * m_1 - m2 ) + g21 * g22 / ( m_2 * m_2 - m2 ) +
291 g31 * g32 / ( m_3 * m_3 - m2 ) + g41 * g42 / ( m_4 * m_4 - m2 ) +
292 g51 * g52 / ( m_5 * m_5 - m2 ) + fs12 * ( 1 - ss0 ) / ( m2 - ss0 ) ) *
293 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
294 double km13 = ( g11 * g13 / ( m_1 * m_1 - m2 ) + g21 * g23 / ( m_2 * m_2 - m2 ) +
295 g31 * g33 / ( m_3 * m_3 - m2 ) + g41 * g43 / ( m_4 * m_4 - m2 ) +
296 g51 * g53 / ( m_5 * m_5 - m2 ) + fs13 * ( 1 - ss0 ) / ( m2 - ss0 ) ) *
297 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
298 double km14 = ( g11 * g14 / ( m_1 * m_1 - m2 ) + g21 * g24 / ( m_2 * m_2 - m2 ) +
299 g31 * g34 / ( m_3 * m_3 - m2 ) + g41 * g44 / ( m_4 * m_4 - m2 ) +
300 g51 * g54 / ( m_5 * m_5 - m2 ) + fs14 * ( 1 - ss0 ) / ( m2 - ss0 ) ) *
301 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
302 double km15 = ( g11 * g15 / ( m_1 * m_1 - m2 ) + g21 * g25 / ( m_2 * m_2 - m2 ) +
303 g31 * g35 / ( m_3 * m_3 - m2 ) + g41 * g45 / ( m_4 * m_4 - m2 ) +
304 g51 * g55 / ( m_5 * m_5 - m2 ) + fs15 * ( 1 - ss0 ) / ( m2 - ss0 ) ) *
305 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
306 double km21 = km12, km31 = km13, km41 = km14, km51 = km15;
307 double km23 = ( g12 * g13 / ( m_1 * m_1 - m2 ) + g22 * g23 / ( m_2 * m_2 - m2 ) +
308 g32 * g33 / ( m_3 * m_3 - m2 ) + g42 * g43 / ( m_4 * m_4 - m2 ) +
309 g52 * g53 / ( m_5 * m_5 - m2 ) ) *
310 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
311 double km24 = ( g12 * g14 / ( m_1 * m_1 - m2 ) + g22 * g24 / ( m_2 * m_2 - m2 ) +
312 g32 * g34 / ( m_3 * m_3 - m2 ) + g42 * g44 / ( m_4 * m_4 - m2 ) +
313 g52 * g54 / ( m_5 * m_5 - m2 ) ) *
314 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
315 double km25 = ( g12 * g15 / ( m_1 * m_1 - m2 ) + g22 * g25 / ( m_2 * m_2 - m2 ) +
316 g32 * g35 / ( m_3 * m_3 - m2 ) + g42 * g45 / ( m_4 * m_4 - m2 ) +
317 g52 * g55 / ( m_5 * m_5 - m2 ) ) *
318 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
319 double km32 = km23, km42 = km24, km52 = km25;
320 double km34 = ( g13 * g14 / ( m_1 * m_1 - m2 ) + g23 * g24 / ( m_2 * m_2 - m2 ) +
321 g33 * g34 / ( m_3 * m_3 - m2 ) + g43 * g44 / ( m_4 * m_4 - m2 ) +
322 g53 * g54 / ( m_5 * m_5 - m2 ) ) *
323 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
324 double km35 = ( g13 * g15 / ( m_1 * m_1 - m2 ) + g23 * g25 / ( m_2 * m_2 - m2 ) +
325 g33 * g35 / ( m_3 * m_3 - m2 ) + g43 * g45 / ( m_4 * m_4 - m2 ) +
326 g53 * g55 / ( m_5 * m_5 - m2 ) ) *
327 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
328 double km43 = km34, km53 = km35;
329 double km45 = ( g14 * g15 / ( m_1 * m_1 - m2 ) + g24 * g25 / ( m_2 * m_2 - m2 ) +
330 g34 * g35 / ( m_3 * m_3 - m2 ) + g44 * g45 / ( m_4 * m_4 - m2 ) +
331 g54 * g55 / ( m_5 * m_5 - m2 ) ) *
332 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
333 double km54 = km45;
334 double km22 = ( g12 * g12 / ( m_1 * m_1 - m2 ) + g22 * g22 / ( m_2 * m_2 - m2 ) +
335 g32 * g32 / ( m_3 * m_3 - m2 ) + g42 * g42 / ( m_4 * m_4 - m2 ) +
336 g52 * g52 / ( m_5 * m_5 - m2 ) ) *
337 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
338 double km33 = ( g13 * g13 / ( m_1 * m_1 - m2 ) + g23 * g23 / ( m_2 * m_2 - m2 ) +
339 g33 * g33 / ( m_3 * m_3 - m2 ) + g43 * g43 / ( m_4 * m_4 - m2 ) +
340 g53 * g53 / ( m_5 * m_5 - m2 ) ) *
341 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
342 double km44 = ( g14 * g14 / ( m_1 * m_1 - m2 ) + g24 * g24 / ( m_2 * m_2 - m2 ) +
343 g34 * g34 / ( m_3 * m_3 - m2 ) + g44 * g44 / ( m_4 * m_4 - m2 ) +
344 g54 * g54 / ( m_5 * m_5 - m2 ) ) *
345 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
346 double km55 = ( g15 * g15 / ( m_1 * m_1 - m2 ) + g25 * g25 / ( m_2 * m_2 - m2 ) +
347 g35 * g35 / ( m_3 * m_3 - m2 ) + g45 * g45 / ( m_4 * m_4 - m2 ) +
348 g55 * g55 / ( m_5 * m_5 - m2 ) ) *
349 ( 1 - sA0 ) / ( m2 - sA0 ) * ( m2 - sA * mass_Pion * mass_Pion * 0.5 );
350 TComplex fp11( 16.55, -3.87 );
351 TComplex beta1( 8.89, 3.36 );
352 TComplex fp12( -13.99, 5.12 );
353 TComplex beta2( 16.20, 5.77 );
354 TComplex fp13( -77.34, -80.04 );
355 TComplex beta3( -21.60, 27.41 );
356 TComplex fp14( -28.52, -3.20 );
357 TComplex beta4( -40.15, 35.36 );
358 TComplex fp15( 30.47, 5.70 );
359 TComplex beta5( 30.0, -43.09 );
360 TComplex P1 = fp11 * ( 1 - sp0 ) / ( m2 - sp0 ) + beta1 * g11 / ( m_1 * m_1 - m2 ) +
361 beta2 * g21 / ( m_2 * m_2 - m2 ) + beta3 * g31 / ( m_3 * m_3 - m2 ) +
362 beta4 * g41 / ( m_4 * m_4 - m2 ) + beta5 * g51 / ( m_5 * m_5 - m2 );
363 TComplex P2 = fp12 * ( 1 - sp0 ) / ( m2 - sp0 ) + beta1 * g12 / ( m_1 * m_1 - m2 ) +
364 beta2 * g22 / ( m_2 * m_2 - m2 ) + beta3 * g32 / ( m_3 * m_3 - m2 ) +
365 beta4 * g42 / ( m_4 * m_4 - m2 ) + beta5 * g52 / ( m_5 * m_5 - m2 );
366 TComplex P3 = fp13 * ( 1 - sp0 ) / ( m2 - sp0 ) + beta1 * g13 / ( m_1 * m_1 - m2 ) +
367 beta2 * g23 / ( m_2 * m_2 - m2 ) + beta3 * g33 / ( m_3 * m_3 - m2 ) +
368 beta4 * g43 / ( m_4 * m_4 - m2 ) + beta5 * g53 / ( m_5 * m_5 - m2 );
369 TComplex P4 = fp14 * ( 1 - sp0 ) / ( m2 - sp0 ) + beta1 * g14 / ( m_1 * m_1 - m2 ) +
370 beta2 * g24 / ( m_2 * m_2 - m2 ) + beta3 * g34 / ( m_3 * m_3 - m2 ) +
371 beta4 * g44 / ( m_4 * m_4 - m2 ) + beta5 * g54 / ( m_5 * m_5 - m2 );
372 TComplex P5 = fp15 * ( 1 - sp0 ) / ( m2 - sp0 ) + beta1 * g15 / ( m_1 * m_1 - m2 ) +
373 beta2 * g25 / ( m_2 * m_2 - m2 ) + beta3 * g35 / ( m_3 * m_3 - m2 ) +
374 beta4 * g45 / ( m_4 * m_4 - m2 ) + beta5 * g55 / ( m_5 * m_5 - m2 );
375
376 TMatrixD MI( 5, 5 ), MA( 5, 5 ), MA_invt( 5, 5 ), MB( 5, 5 ), KM( 5, 5 ), RhoA( 5, 5 ),
377 RhoB( 5, 5 ), MRe( 5, 5 ), MIm( 5, 5 );
378 MI.UnitMatrix();
379
380 RhoA.UnitMatrix();
381 RhoB.UnitMatrix();
382 TMatrixDRow( KM, 0 )( 0 ) = km11;
383 TMatrixDRow( KM, 1 )( 0 ) = km21;
384 TMatrixDRow( KM, 2 )( 0 ) = km31;
385 TMatrixDRow( KM, 3 )( 0 ) = km41;
386 TMatrixDRow( KM, 4 )( 0 ) = km51;
387 TMatrixDRow( KM, 0 )( 1 ) = km12;
388 TMatrixDRow( KM, 1 )( 1 ) = km22;
389 TMatrixDRow( KM, 2 )( 1 ) = km32;
390 TMatrixDRow( KM, 3 )( 1 ) = km42;
391 TMatrixDRow( KM, 4 )( 1 ) = km52;
392 TMatrixDRow( KM, 0 )( 2 ) = km13;
393 TMatrixDRow( KM, 1 )( 2 ) = km23;
394 TMatrixDRow( KM, 2 )( 2 ) = km33;
395 TMatrixDRow( KM, 3 )( 2 ) = km43;
396 TMatrixDRow( KM, 4 )( 2 ) = km53;
397 TMatrixDRow( KM, 0 )( 3 ) = km14;
398 TMatrixDRow( KM, 1 )( 3 ) = km24;
399 TMatrixDRow( KM, 2 )( 3 ) = km34;
400 TMatrixDRow( KM, 3 )( 3 ) = km44;
401 TMatrixDRow( KM, 4 )( 3 ) = km54;
402 TMatrixDRow( KM, 0 )( 4 ) = km15;
403 TMatrixDRow( KM, 1 )( 4 ) = km25;
404 TMatrixDRow( KM, 2 )( 4 ) = km35;
405 TMatrixDRow( KM, 3 )( 4 ) = km45;
406 TMatrixDRow( KM, 4 )( 4 ) = km55;
407
408 TMatrixDRow( RhoA, 0 )( 0 ) = sqrt( 1.0 - 4.0 * mass_Pion * mass_Pion / m2 );
409 TMatrixDRow( RhoB, 0 )( 0 ) = 0.0;
410 if ( ( 1.0 - 4.0 * mass_Kaon * mass_Kaon / m2 ) > 0 )
411 {
412 TMatrixDRow( RhoA, 1 )( 1 ) = sqrt( 1.0 - 4.0 * mass_Kaon * mass_Kaon / m2 );
413 TMatrixDRow( RhoB, 1 )( 1 ) = 0.0;
414 }
415 else
416 {
417 TMatrixDRow( RhoA, 1 )( 1 ) = 0.0;
418 TMatrixDRow( RhoB, 1 )( 1 ) = sqrt( 4.0 * mass_Kaon * mass_Kaon / m2 - 1.0 );
419 }
420 TMatrixDRow( RhoA, 2 )( 2 ) = CalRho4pi( m2 );
421 TMatrixDRow( RhoB, 2 )( 2 ) = 0.0;
422 if ( ( 1.0 - 4.0 * mass_Eta * mass_Eta / m2 ) > 0 )
423 {
424 TMatrixDRow( RhoA, 3 )( 3 ) = sqrt( 1.0 - 4.0 * mass_Eta * mass_Eta / m2 );
425 TMatrixDRow( RhoB, 3 )( 3 ) = 0.0;
426 }
427 else
428 {
429 TMatrixDRow( RhoA, 3 )( 3 ) = 0.0;
430 TMatrixDRow( RhoB, 3 )( 3 ) = sqrt( 4.0 * mass_Eta * mass_Eta / m2 - 1.0 );
431 }
432 TMatrixDRow( RhoA, 4 )( 4 ) = 0.0;
433 TMatrixDRow( RhoB, 4 )( 4 ) =
434 sqrt( ( mass_Eta + mass_Etap ) * ( mass_Eta + mass_Etap ) / m2 - 1.0 );
435
436 MA = MI + KM * RhoB;
437 MB = -1.0 * KM * RhoA;
438 MA_invt = MA;
439 MA_invt.Invert();
440 MRe = MA + MB * MA_invt * MB;
441 MRe.Invert();
442 MIm = MA_invt * MB * MRe;
443 TComplex ciR( 1.0, 0.0 );
444 TComplex ciM( 0.0, 1.0 );
445 TComplex amp;
446 amp = ( ciR * TMatrixDRow( MRe, 0 )( 0 ) - ciM * TMatrixDRow( MIm, 0 )( 0 ) ) * P1 +
447 ( ciR * TMatrixDRow( MRe, 0 )( 1 ) - ciM * TMatrixDRow( MIm, 0 )( 1 ) ) * P2 +
448 ( ciR * TMatrixDRow( MRe, 0 )( 2 ) - ciM * TMatrixDRow( MIm, 0 )( 2 ) ) * P3 +
449 ( ciR * TMatrixDRow( MRe, 0 )( 3 ) - ciM * TMatrixDRow( MIm, 0 )( 3 ) ) * P4 +
450 ( ciR * TMatrixDRow( MRe, 0 )( 4 ) - ciM * TMatrixDRow( MIm, 0 )( 4 ) ) * P5;
451 return amp;
452}
453void EvtDsToKpPipPimPi0::KPiSLASS( double sa, double sb, double sc, double prop[2] ) {
454 const double m1430 = 1.463;
455 const double sa0 = 2.140369; // m1430*m1430;
456 const double w1430 = 0.233;
457 const double Lass1 = 0.25 / sa0;
458 double tmp = sb - sc;
459 double tmp1 = sa0 + tmp;
460 double q0 = Lass1 * tmp1 * tmp1 - sb;
461 if ( q0 < 0 ) q0 = 1e-16;
462 double tmp2 = sa + tmp;
463 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
464 double q = sqrt( qs );
465 double width = w1430 * q * m1430 / sqrt( sa * q0 );
466 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
467 if ( temp_R < 0 ) temp_R += math_pi;
468 double deltaR = -5.31 + temp_R;
469 double temp_F = atan( 2.14 * q / ( 2.0 - 1.926 * qs ) ); // 2.0*1.07 = 2.14; 1.8*1.07 = 1.926
470 if ( temp_F < 0 ) temp_F += math_pi;
471 double deltaF = 2.33 + temp_F;
472 double deltaS = deltaR + 2.0 * deltaF;
473 double t1 = 0.8 * sin( deltaF );
474 double t2 = sin( deltaR );
475 double CF[2], CS[2];
476 CF[0] = cos( deltaF );
477 CF[1] = sin( deltaF );
478 CS[0] = cos( deltaS );
479 CS[1] = sin( deltaS );
480 prop[0] = t1 * CF[0] + t2 * CS[0];
481 prop[1] = t1 * CF[1] + t2 * CS[1];
482}
483
484void EvtDsToKpPipPimPi0::Com_Multi( double a1[2], double a2[2], double res[2] ) {
485 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
486 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
487}
488void EvtDsToKpPipPimPi0::Com_Divide( double a1[2], double a2[2], double res[2] ) {
489 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / ( a2[0] * a2[0] + a2[1] * a2[1] );
490 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / ( a2[0] * a2[0] + a2[1] * a2[1] );
491}
492double EvtDsToKpPipPimPi0::SCADot( double a1[4], double a2[4] ) {
493 double _cal = 0.;
494 _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
495 return _cal;
496}
497double EvtDsToKpPipPimPi0::barrier( int l, double sa, double sb, double sc, double r,
498 double mass ) {
499 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
500 // if(q < 0) q = 1e-16;
501 if ( q < 0 ) q = -q;
502
503 double z;
504 z = q * r * r;
505 double sa0;
506 sa0 = mass * mass;
507 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
508 // if(q0 < 0) q0 = 1e-16;
509 if ( q0 < 0 ) q0 = -q0;
510 double z0 = q0 * r * r;
511 double F = 0.0;
512 if ( l == 0 ) F = 1;
513 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
514 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
515 // if(l == 1) F = sqrt((2*r)/(1+z));
516 // if(l == 2) F = sqrt((13*r*r)/(9+3*z+z*z));
517 return F;
518}
519void EvtDsToKpPipPimPi0::calt1( double daug1[4], double daug0[4], double t1[4] ) {
520 double p, pq;
521 double pa[4], qa[4];
522 for ( int i = 0; i < 4; i++ )
523 {
524 pa[i] = daug1[i] + daug0[i];
525 qa[i] = daug1[i] - daug0[i];
526 }
527 p = SCADot( pa, pa );
528 pq = SCADot( pa, qa );
529 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - pq / p * pa[i]; }
530}
531void EvtDsToKpPipPimPi0::calt2( double daug1[4], double daug0[4], double t2[4][4] ) {
532 double p, r;
533 double pa[4], t1[4];
534 calt1( daug1, daug0, t1 );
535 r = SCADot( t1, t1 ) / 3.0;
536 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug0[i]; }
537 p = SCADot( pa, pa );
538 for ( int i = 0; i < 4; i++ )
539 {
540 for ( int j = 0; j < 4; j++ )
541 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
542 }
543}
544
545double EvtDsToKpPipPimPi0::wid( double mass2, double mass, double sa, double sb, double sc,
546 double r, int l ) {
547 double widm = 0.;
548 double m = sqrt( sa );
549 double tmp = sb - sc;
550 double tmp1 = sa + tmp;
551 double q = 0.25 * tmp1 * tmp1 / sa - sb;
552 // if(q<0) q = 1e-16;
553 if ( q < 0 ) q = -q;
554 double tmp2 = mass2 + tmp;
555 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
556 // if(q0<0) q0 = 1e-16;
557 if ( q0 < 0 ) q0 = -q0;
558 double z = q * r * r;
559 double z0 = q0 * r * r;
560 double t = q / q0;
561 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
562 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
563 else if ( l == 2 )
564 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
565 return widm;
566}
567double EvtDsToKpPipPimPi0::widl1( double mass2, double mass, double sa, double sb, double sc,
568 double r ) {
569 double widm = 0.;
570 double m = sqrt( sa );
571 double tmp = sb - sc;
572 double tmp1 = sa + tmp;
573 double q = 0.25 * tmp1 * tmp1 / sa - sb;
574 // if(q<0) q = 1e-16;
575 if ( q < 0 ) q = -q;
576 double tmp2 = mass2 + tmp;
577 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
578 // if(q0<0) q0 = 1e-16;
579 if ( q0 < 0 ) q0 = -q0;
580 double z = q * r * r;
581 double z0 = q0 * r * r;
582 double F = ( 1 + z0 ) / ( 1 + z );
583 double t = q / q0;
584 widm = t * sqrt( t ) * mass / m * F;
585 return widm;
586}
587void EvtDsToKpPipPimPi0::propagatorNBW( double mass2, double mass, double width, double sa,
588 double sb, double sc, double r, int l,
589 double prop[2] ) {
590 double a[2], b[2];
591 a[0] = 1;
592 a[1] = 0;
593 b[0] = mass2 - sa;
594 b[1] = -mass * width;
595 Com_Divide( a, b, prop );
596}
597
598void EvtDsToKpPipPimPi0::propagatorRBW( double mass2, double mass, double width, double sa,
599 double sb, double sc, double r, int l,
600 double prop[2] ) {
601 double a[2], b[2];
602 a[0] = 1;
603 a[1] = 0;
604 b[0] = mass2 - sa;
605 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r, l );
606 Com_Divide( a, b, prop );
607}
608void EvtDsToKpPipPimPi0::propagatorRBWl1( double mass2, double mass, double width, double sa,
609 double sb, double sc, double r, double prop[2] )
610// firstly code for propagator_omg, now could use for propagatorRBW that l=1, and will be more
611// quick.
612{
613 double a[2], b[2];
614 a[0] = 1;
615 a[1] = 0;
616 b[0] = mass2 - sa;
617 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r );
618 Com_Divide( a, b, prop );
619}
620//------------GS---used by rho----------------------------
621void EvtDsToKpPipPimPi0::propagatorGS( double mass2, double mass, double width, double sa,
622 double sb, double sc, double r, double prop[2] ) {
623 double a[2], b[2];
624 double tmp = sb - sc;
625 double tmp1 = sa + tmp;
626 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
627 // if(q2<0) q2 = 1e-16;
628 if ( q2 < 0 ) q2 = -q2;
629
630 double tmp2 = mass2 + tmp;
631 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
632 // if(q02<0) q02 = 1e-16;
633 if ( q02 < 0 ) q02 = -q02;
634
635 double q = sqrt( q2 );
636 double q0 = sqrt( q02 );
637 double m = sqrt( sa );
638 double q03 = q0 * q02;
639 double tmp3 = log( mass + 2 * q0 ) + 1.2760418309; // log(mass_2Pion) = 1.2760418309;
640
641 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2760418309 );
642 double h0 = GS1 * q0 / mass * tmp3;
643 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
644 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
645 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
646
647 a[0] = 1.0 + d * width / mass;
648 a[1] = 0.0;
649 b[0] = mass2 - sa + width * f;
650 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r );
651 Com_Divide( a, b, prop );
652}
653void EvtDsToKpPipPimPi0::calEvaMy( double* Kp, double* Pip, double* Pim, double* Pi0,
654 double* mass1, double* mass2, double* width1,
655 double* width2, double* amp, double* phase, int* g0,
656 int* g1, int* g2, int* modetype, int nstates,
657 double& Result, double& q1270 )
658
659{
660
661 TComplex tmpswave;
662 double spipi( 0. );
663
664 spipi =
665 ( Pip[0] + Pim[0] ) * ( Pip[0] + Pim[0] ) - ( Pip[1] + Pim[1] ) * ( Pip[1] + Pim[1] ) -
666 ( Pip[2] + Pim[2] ) * ( Pip[2] + Pim[2] ) - ( Pip[3] + Pim[3] ) * ( Pip[3] + Pim[3] );
667 tmpswave = ResonanceSkm( spipi );
668 double realpipis = tmpswave.Re();
669 double imagpipis = tmpswave.Im();
670
671 double pKstr0[4], pKstrp[4], prhop[4], prhom[4], prho0[4], pK11[4], pK12[4], pK13[4], pA1[4],
672 pD[4], pomega[4];
673 for ( int i = 0; i != 4; i++ )
674 {
675 prhop[i] = Pip[i] + Pi0[i];
676 prhom[i] = Pim[i] + Pi0[i];
677 prho0[i] = Pim[i] + Pip[i];
678 pKstr0[i] = Kp[i] + Pim[i];
679 pKstrp[i] = Kp[i] + Pi0[i];
680 pK11[i] = pKstr0[i] + Pi0[i];
681 pK12[i] = Kp[i] + Pim[i];
682 pK13[i] = pKstr0[i] + Pip[i];
683 pA1[i] = prhop[i] + Pim[i];
684 pomega[i] = Pim[i] + Pip[i] + Pi0[i];
685 pD[i] = pKstr0[i] + prhop[i];
686 }
687 double spi0, spionp, spionm, skaon, srhop, srhom, srho0, sKst0, sKstp, sk11, sk12, sk13, sa1,
688 sD, somega;
689
690 skaon = SCADot( Kp, Kp );
691 spionp = SCADot( Pip, Pip );
692 spionm = SCADot( Pim, Pim );
693 spi0 = SCADot( Pi0, Pi0 );
694
695 srhop = SCADot( prhop, prhop );
696 srhom = SCADot( prhom, prhom );
697 srho0 = SCADot( prho0, prho0 );
698 somega = SCADot( pomega, pomega );
699
700 sKst0 = SCADot( pKstr0, pKstr0 );
701 sKstp = SCADot( pKstrp, pKstrp );
702 sk11 = SCADot( pK11, pK11 );
703 sk12 = SCADot( pK12, pK12 );
704 sk13 = SCADot( pK13, pK13 );
705 sa1 = SCADot( pA1, pA1 );
706 sD = SCADot( pD, pD );
707 //-------------------------------------------------------------------------------------------------
708 double t2A1[4][4], t2A2[4][4], t2A3[4][4], t1rhom[4], t1rhop[4], t1rho0[4], t1Kst0[4],
709 t1Kstp[4], t2k11[4][4], t2k12[4][4];
710 double t1K1_KPi[4], t2k21[4][4], t2k22[4][4], t1A3[4];
711 calt1( Pip, Pi0, t1rhop );
712 calt1( Pim, Pi0, t1rhom );
713 calt1( Pim, Pip, t1rho0 );
714
715 calt1( Kp, Pim, t1Kst0 );
716 calt1( Kp, Pi0, t1Kstp );
717 calt1( Kp, Pim, t1K1_KPi );
718
719 calt1( prho0, Pi0, t1A3 );
720
721 calt2( prhop, Pim, t2A1 );
722 calt2( prhom, Pip, t2A2 );
723 calt2( prho0, Pi0, t2A3 );
724
725 calt2( pKstr0, Pi0, t2k11 );
726 calt2( pKstrp, Pim, t2k12 );
727
728 calt2( prhom, Kp, t2k21 );
729 calt2( prho0, Kp, t2k22 );
730
731 double B_rhop = -1.0, B_rhom = -1.0, B_Kst0rhop1 = -1.0, B_Kst0rhop2 = -1.0;
732 double B_Kst0 = -1.0, B_K1_KPi = -1.0, B_K1rhop1 = -1.0, B_K1rhop2 = -1.0;
733 double B_Kstp = -1.0, B_K1rho = -1.0, B_D_K1rho = -1.0, B_RhopKs1 = -1.0;
734 double B_K1_Kst0pi0 = -1.0, B_D_K1Pi = -1.0, B_K1_Kstppim = -1.0, B_K1Kstp2 = -1.0;
735 double B_D_A1 = -1.0, B_D_A2 = -1.0, B_rho0 = -1.0, B_Rho0Ks1 = -1.0;
736 double B_K1PiPi = -1.0, B_K14_1 = -1.0, B_K14_2 = -1.0;
737 double B_Kstprho01 = -1.0, B_Kstprho02 = -1.0;
738
739 double cof[2], amp_tmp[2], amp_PDF[2], PDF[2];
740 PDF[0] = 0;
741 PDF[1] = 0;
742 double mass1sq, mass2sq, tt1, tt2, tt21, tt22, tt23, tt3, qqq = 10, qqq0 = 10, nq = 0;
743 double temp_PDF, tmp1, tmp3, temp_PDF1;
744 double pro[2], pro0[2], pro1[2], pro2[2], pro3[2], pro4[2];
745 double t1D[4], t2D[4][4], t1Ds_omega[4], B1_Ds_omega, B1_1V23;
746 int isKstp = 0.0, isKst0 = 0.0, isRhop = 0.0, isRhom = 0.0, isRho0 = 0.0, isKpPim_S = 0.0,
747 isPipPi0_S = 0.0, isf0 = 0.0, isPipPim_S = 0.0;
748 double proRhop[2], proRhom[2], proRho0[2], proKstp[2], proKst0[2], proKPi_S[2], proPiPi_S[2],
749 prof0[2];
750 double temp_PDF11, temp_PDF12, temp_PDF13;
751 double pro0_11[2], pro0_12[2], pro0_13[2];
752 double pro_11[2], pro_12[2], pro_13[2];
753
754 for ( int i = 0; i < nstates; i++ )
755 {
756 temp_PDF = 0;
757 temp_PDF1 = 0;
758 temp_PDF11 = 0;
759 temp_PDF12 = 0;
760 temp_PDF13 = 0;
761 cof[0] = amp[i] * cos( phase[i] );
762 cof[1] = amp[i] * sin( phase[i] );
763 mass1sq = mass1[i] * mass1[i]; // avoid parameters input reversed
764 mass2sq = mass2[i] * mass2[i];
765
766 if ( modetype[i] == 1 )
767 {
768 if ( B_Kst0 < 0.0 ) B_Kst0 = barrier( 1, sKst0, skaon, spionm, rRes1, mass1[i] );
769 if ( B_rhop < 0.0 ) B_rhop = barrier( 1, srhop, spionp, spi0, rRes1, mass2[i] );
770 if ( g2[i] == 0 )
771 {
772 for ( int w = 0; w < 4; w++ ) { temp_PDF += G[w][w] * t1Kst0[w] * t1rhop[w]; }
773 tmp1 = B_Kst0 * B_rhop * temp_PDF;
774 }
775 else if ( g2[i] == 1 )
776 {
777 calt1( pKstr0, prhop, t1D );
778 for ( int w = 0; w < 4; w++ )
779 {
780 tt1 = pD[w] * G[w][w];
781 for ( int j = 0; j < 4; j++ )
782 {
783 tt2 = t1D[j] * G[j][j];
784 for ( int k = 0; k < 4; k++ )
785 {
786 tt3 = t1Kst0[k] * G[k][k];
787 for ( int l = 0; l < 4; l++ )
788 { temp_PDF += E[w][j][k][l] * tt1 * tt2 * tt3 * t1rhop[l] * G[l][l]; }
789 }
790 }
791 }
792 if ( B_Kst0rhop1 < 0.0 )
793 {
794 B_Kst0rhop1 = barrier( 1, sD, sKst0, srhop, rD2, mDsM );
795 tmp1 = B_Kst0 * B_rhop * B_Kst0rhop1 * temp_PDF;
796 }
797 }
798 else if ( g2[i] == 2 )
799 {
800 calt2( pKstr0, prhop, t2D );
801 for ( int w = 0; w < 4; w++ )
802 {
803 tt1 = t1Kst0[w] * G[w][w];
804 for ( int j = 0; j < 4; j++ ) { temp_PDF += t2D[w][j] * t1rhop[j] * G[j][j] * tt1; }
805 }
806 if ( B_Kst0rhop2 < 0.0 )
807 {
808 B_Kst0rhop2 = barrier( 2, sD, sKst0, srhop, rD2, mDsM );
809 tmp1 = B_Kst0 * B_rhop * B_Kst0rhop2 * temp_PDF;
810 }
811 }
812 if ( isRhop == 0 )
813 {
814 propagatorGS( mass2sq, mass2[i], width2[i], srhop, spionp, spi0, rRes1,
815 proRhop ); // rho as mass2
816 isRhop = 1.0;
817 }
818 if ( g0[i] == 1 )
819 { propagatorRBWl1( mass1sq, mass1[i], width1[i], sKst0, skaon, spionm, rRes1, pro0 ); }
820 else if ( g0[i] == 0 )
821 {
822 pro0[0] = 1;
823 pro0[1] = 0;
824 }
825 if ( g1[i] == 1 )
826 {
827 pro1[0] = proRhop[0];
828 pro1[1] = proRhop[1];
829 }
830 else if ( g1[i] == 0 )
831 {
832 pro1[0] = 1;
833 pro1[1] = 0;
834 }
835 Com_Multi( pro0, pro1, pro );
836 amp_tmp[0] = tmp1 * pro[0];
837 amp_tmp[1] = tmp1 * pro[1];
838 }
839 if ( modetype[i] == 2 )
840 {
841 if ( B_Kstp < 0.0 ) B_Kstp = barrier( 1, sKstp, skaon, spi0, rRes1, mass1[i] );
842 if ( B_rho0 < 0.0 ) B_rho0 = barrier( 1, srho0, spionp, spionm, rRes1, mass2[i] );
843 if ( g2[i] == 0 )
844 {
845 for ( int w = 0; w < 4; w++ ) { temp_PDF += G[w][w] * t1Kstp[w] * t1rho0[w]; }
846 tmp1 = B_Kstp * B_rho0 * temp_PDF;
847 }
848 else if ( g2[i] == 1 )
849 {
850 calt1( pKstrp, prho0, t1D );
851 for ( int w = 0; w < 4; w++ )
852 {
853 tt1 = pD[w] * G[w][w];
854 for ( int j = 0; j < 4; j++ )
855 {
856 tt2 = t1D[j] * G[j][j];
857 for ( int k = 0; k < 4; k++ )
858 {
859 tt3 = t1Kstp[k] * G[k][k];
860 for ( int l = 0; l < 4; l++ )
861 { temp_PDF += E[w][j][k][l] * tt1 * tt2 * tt3 * t1rho0[l] * G[l][l]; }
862 }
863 }
864 }
865 if ( B_Kstprho01 < 0.0 )
866 {
867 B_Kstprho01 = barrier( 1, sD, sKstp, srho0, rD2, mDsM );
868 tmp1 = B_Kstp * B_rho0 * B_Kstprho01 * temp_PDF;
869 }
870 }
871 else if ( g2[i] == 2 )
872 {
873 calt2( pKstrp, prho0, t2D );
874 for ( int w = 0; w < 4; w++ )
875 {
876 tt1 = t1Kstp[w] * G[w][w];
877 for ( int j = 0; j < 4; j++ ) { temp_PDF += t2D[w][j] * t1rho0[j] * G[j][j] * tt1; }
878 }
879 if ( B_Kstprho02 < 0.0 )
880 {
881 B_Kstprho02 = barrier( 2, sD, sKstp, srho0, rD2, mDsM );
882 tmp1 = B_Kstp * B_rho0 * B_Kstprho02 * temp_PDF;
883 }
884 }
885 if ( isRho0 == 0 )
886 {
887 propagatorGS( mass2sq, mass2[i], width2[i], srho0, spionp, spionm, rRes1,
888 proRho0 ); // rho as mass2
889 isRho0 = 1.0;
890 }
891 if ( g0[i] == 1 )
892 { propagatorRBWl1( mass1sq, mass1[i], width1[i], sKstp, skaon, spi0, rRes1, pro0 ); }
893 else if ( g0[i] == 0 )
894 {
895 pro0[0] = 1;
896 pro0[1] = 0;
897 }
898 if ( g1[i] == 1 )
899 {
900 pro1[0] = proRho0[0];
901 pro1[1] = proRho0[1];
902 }
903 else if ( g1[i] == 0 )
904 {
905 pro1[0] = 1;
906 pro1[1] = 0;
907 }
908 Com_Multi( pro0, pro1, pro );
909 amp_tmp[0] = tmp1 * pro[0];
910 amp_tmp[1] = tmp1 * pro[1];
911 }
912 else if ( modetype[i] == 6 )
913 {
914 if ( B_D_A1 < 0.0 ) B_D_A1 = barrier( 1, sD, sa1, skaon, rD2, mDsM );
915 if ( B_rhop < 0.0 ) B_rhop = barrier( 1, srhop, spionp, spi0, rRes1, mass2[i] );
916 calt1( pA1, Kp, t1D );
917 if ( g2[i] == 0 )
918 {
919 for ( int w = 0; w < 4; w++ )
920 {
921 tt1 = t1D[w] * G[w][w];
922 for ( int j = 0; j < 4; j++ )
923 { temp_PDF += tt1 * ( G[w][j] - pA1[w] * pA1[j] / sa1 ) * t1rhop[j] * G[j][j]; }
924 }
925 tmp1 = B_rhop * B_D_A1 * temp_PDF;
926 }
927 else if ( g2[i] == 2 )
928 {
929 for ( int w = 0; w < 4; w++ )
930 {
931 tt1 = t1D[w] * G[w][w];
932 for ( int j = 0; j < 4; j++ ) { temp_PDF += tt1 * t2A1[w][j] * t1rhop[j] * G[j][j]; }
933 }
934 if ( B_D_A2 < 0.0 ) B_D_A2 = barrier( 2, sa1, srhop, spionm, rRes1, mass1[i] );
935 tmp1 = B_rhop * B_D_A2 * B_D_A1 * temp_PDF;
936 }
937 if ( isRhop == 0 )
938 {
939 propagatorGS( mass2sq, mass2[i], width2[i], srhop, spionp, spi0, rRes1, proRhop );
940 isRhop = 1;
941 }
942 if ( g0[i] == 1 )
943 propagatorRBW( mass1sq, mass1[i], width1[i], sa1, srhop, spionm, rRes1, g2[i], pro0 );
944 else if ( g0[i] == 0 )
945 {
946 pro0[0] = 1;
947 pro0[1] = 0;
948 }
949 if ( g1[i] == 1 )
950 {
951 pro1[0] = proRhop[0];
952 pro1[1] = proRhop[1];
953 }
954 else if ( g1[i] == 0 )
955 {
956 pro1[0] = 1;
957 pro1[1] = 0;
958 }
959 Com_Multi( pro0, pro1, pro );
960 amp_tmp[0] = tmp1 * pro[0];
961 amp_tmp[1] = tmp1 * pro[1];
962 }
963 else if ( modetype[i] == 7 )
964 {
965 if ( B_D_A1 < 0.0 ) B_D_A1 = barrier( 1, sD, sa1, skaon, rD2, mDsM );
966 if ( B_rhom < 0.0 ) B_rhom = barrier( 1, srhom, spionm, spi0, rRes1, mass2[i] );
967 calt1( pA1, Kp, t1D );
968 if ( g2[i] == 0 )
969 {
970 for ( int w = 0; w < 4; w++ )
971 {
972 tt1 = t1D[w] * G[w][w];
973 for ( int j = 0; j < 4; j++ )
974 { temp_PDF += tt1 * ( G[w][j] - pA1[w] * pA1[j] / sa1 ) * t1rhom[j] * G[j][j]; }
975 }
976 tmp1 = B_rhom * B_D_A1 * temp_PDF;
977 }
978 else if ( g2[i] == 2 )
979 {
980 for ( int w = 0; w < 4; w++ )
981 {
982 tt1 = t1D[w] * G[w][w];
983 for ( int j = 0; j < 4; j++ ) { temp_PDF += tt1 * t2A2[w][j] * t1rhom[j] * G[j][j]; }
984 }
985 if ( B_D_A2 < 0.0 ) B_D_A2 = barrier( 2, sa1, srhom, spionp, rRes1, mass1[i] );
986 tmp1 = B_rhom * B_D_A2 * B_D_A1 * temp_PDF;
987 }
988 if ( isRhom == 0 )
989 {
990 propagatorGS( mass2sq, mass2[i], width2[i], srhom, spionm, spi0, rRes1, proRhom );
991 isRhom = 1;
992 }
993 if ( g0[i] == 1 )
994 propagatorRBW( mass1sq, mass1[i], width1[i], sa1, srhom, spionp, rRes1, g2[i], pro0 );
995 else if ( g0[i] == 0 )
996 {
997 pro0[0] = 1;
998 pro0[1] = 0;
999 }
1000 if ( g1[i] == 1 )
1001 {
1002 pro1[0] = proRhom[0];
1003 pro1[1] = proRhom[1];
1004 }
1005 else if ( g1[i] == 0 )
1006 {
1007 pro1[0] = 1;
1008 pro1[1] = 0;
1009 }
1010 Com_Multi( pro0, pro1, pro );
1011 amp_tmp[0] = -tmp1 * pro[0];
1012 amp_tmp[1] = -tmp1 * pro[1];
1013 }
1014
1015 else if ( modetype[i] == 10 )
1016 {
1017 if ( B_rhom < 0.0 ) B_rhom = barrier( 1, srhom, spionm, spi0, rRes1, mass2[i] );
1018 // printf("B_rhom= %.15f\n",B_rhom);
1019 if ( B_D_K1rho < 0.0 ) B_D_K1rho = barrier( 1, sD, sk11, spionp, rD2, mDsM );
1020 // printf("B_D_K1rho= %.15f\n",B_D_K1rho);
1021 calt1( pK11, Pip, t1D );
1022 if ( g2[i] == 0 )
1023 {
1024 for ( int w = 0; w < 4; w++ )
1025 {
1026 tt1 = t1D[w] * G[w][w];
1027 for ( int j = 0; j < 4; j++ )
1028 { temp_PDF += tt1 * ( G[w][j] - pK11[w] * pK11[j] / sk11 ) * t1rhom[j] * G[j][j]; }
1029 }
1030 tmp1 = B_rhom * B_D_K1rho * temp_PDF;
1031 }
1032 else if ( g2[i] == 2 )
1033 {
1034 for ( int w = 0; w < 4; w++ )
1035 {
1036 tt1 = t1D[w] * G[w][w];
1037 for ( int j = 0; j < 4; j++ )
1038 { temp_PDF += tt1 * t2k21[w][j] * t1rhom[j] * G[j][j]; }
1039 }
1040 if ( B_K1rho < 0.0 )
1041 {
1042 B_K1rho = barrier( 2, sk11, srhom, skaon, rRes1, mass1[i] );
1043 tmp1 = B_rhom * B_K1rho * B_D_K1rho * temp_PDF;
1044 }
1045 }
1046 if ( isRhom == 0 )
1047 {
1048 propagatorGS( mass2sq, mass2[i], width2[i], srhom, spionm, spi0, rRes1,
1049 proRhom ); // rho also as mass2, but for this process the input parameter
1050 // should change position (differ to bofore)).
1051 isRhom = 1.0;
1052 }
1053 qqq = 0.25 * ( sk11 + srhom - skaon ) * ( sk11 + srhom - skaon ) / sk11 - srhom;
1054 qqq0 = 0.25 * ( mK1270 + srhom - skaon ) * ( mK1270 + srhom - skaon ) / mK1270 - srhom;
1055 if ( qqq < 2e-4 ) continue;
1056 if ( g0[i] == 1 )
1057 propagatorRBW( mass1sq, mass1[i], width1[i], sk11, srhom, skaon, rRes1, g2[i],
1058 pro0 ); // K1270(rho) as mass1 which is different from K1270(K*)(as
1059 // mass2)
1060 else if ( g0[i] == 0 )
1061 {
1062 pro0[0] = 1;
1063 pro0[1] = 0;
1064 }
1065 if ( g1[i] == 1 )
1066 {
1067 pro1[0] = proRhom[0];
1068 pro1[1] = proRhom[1];
1069 }
1070 else if ( g1[i] == 0 )
1071 {
1072 pro1[0] = 1;
1073 pro1[1] = 0;
1074 }
1075 Com_Multi( pro0, pro1, pro );
1076 amp_tmp[0] = tmp1 * pro[0];
1077 amp_tmp[1] = tmp1 * pro[1];
1078 }
1079 else if ( modetype[i] == 11 )
1080 {
1081 if ( B_D_K1Pi < 0.0 ) B_D_K1Pi = barrier( 1, sD, sk11, spionp, rD2, mDsM );
1082 if ( B_Kst0 < 0.0 ) B_Kst0 = barrier( 1, sKst0, skaon, spionm, rRes1, mass1[i] );
1083 calt1( pK11, Pip, t1D );
1084 if ( g2[i] == 0 )
1085 {
1086 for ( int w = 0; w < 4; w++ )
1087 {
1088 tt1 = t1D[w] * G[w][w];
1089 for ( int j = 0; j < 4; j++ )
1090 { temp_PDF += tt1 * ( G[w][j] - pK11[w] * pK11[j] / sk11 ) * t1Kst0[j] * G[j][j]; }
1091 }
1092 tmp1 = B_Kst0 * B_D_K1Pi * temp_PDF;
1093 }
1094 else if ( g2[i] == 2 )
1095 {
1096 for ( int w = 0; w < 4; w++ )
1097 {
1098 tt1 = t1D[w] * G[w][w];
1099 for ( int j = 0; j < 4; j++ )
1100 { temp_PDF += tt1 * t2k11[w][j] * t1Kst0[j] * G[j][j]; }
1101 }
1102 if ( B_K1rhop2 < 0.0 ) B_K1rhop2 = barrier( 2, sk11, sKst0, spi0, rRes1, mass2[i] );
1103 tmp1 = B_Kst0 * B_K1rhop2 * B_D_K1Pi * temp_PDF;
1104 }
1105 if ( isKst0 == 0 )
1106 {
1107 propagatorRBWl1( mass1sq, mass1[i], width1[i], sKst0, skaon, spionm, rRes1,
1108 proKst0 ); // Kst0 also as mass1
1109 isKst0 = 1;
1110 }
1111 if ( g0[i] == 1 )
1112 {
1113 pro0[0] = proKst0[0];
1114 pro0[1] = proKst0[1];
1115 }
1116 else if ( g0[i] == 0 )
1117 {
1118 pro0[0] = 1;
1119 pro0[1] = 0;
1120 }
1121 if ( g1[i] == 1 )
1122 {
1123 propagatorRBW( mass2sq, mass2[i], width2[i], sk11, sKst0, spi0, rRes1, g2[i],
1124 pro1 ); // K1270(K*) as mass2
1125 }
1126 else if ( g1[i] == 0 )
1127 {
1128 pro1[0] = 1;
1129 pro1[1] = 0;
1130 }
1131 Com_Multi( pro0, pro1, pro );
1132 amp_tmp[0] = tmp1 * pro[0];
1133 amp_tmp[1] = tmp1 * pro[1];
1134 }
1135 else if ( modetype[i] == 12 )
1136 {
1137 if ( B_D_K1Pi < 0.0 ) B_D_K1Pi = barrier( 1, sD, sk11, spionp, rD2, mDsM );
1138 if ( B_Kstp < 0.0 ) B_Kstp = barrier( 1, sKstp, skaon, spi0, rRes1, mass1[i] );
1139 calt1( pK11, Pip, t1D );
1140 if ( g2[i] == 0 )
1141 {
1142 for ( int w = 0; w < 4; w++ )
1143 {
1144 tt1 = t1D[w] * G[w][w];
1145 for ( int j = 0; j < 4; j++ )
1146 { temp_PDF1 += tt1 * ( G[w][j] - pK11[w] * pK11[j] / sk11 ) * t1Kstp[j] * G[j][j]; }
1147 }
1148 tmp1 = B_Kstp * B_D_K1Pi * temp_PDF1;
1149 }
1150 else if ( g2[i] == 2 )
1151 {
1152 for ( int w = 0; w < 4; w++ )
1153 {
1154 tt1 = t1D[w] * G[w][w];
1155 for ( int j = 0; j < 4; j++ )
1156 { temp_PDF1 += tt1 * t2k12[w][j] * t1Kstp[j] * G[j][j]; }
1157 }
1158 if ( B_K1Kstp2 < 0.0 ) B_K1Kstp2 = barrier( 2, sk11, sKstp, spionm, rRes1, mass2[i] );
1159 tmp1 = B_Kstp * B_K1Kstp2 * B_D_K1Pi * temp_PDF1;
1160 }
1161 if ( isKstp == 0 )
1162 {
1163 propagatorRBWl1( mass1sq, mass1[i], width1[i], sKstp, skaon, spi0, rRes1,
1164 proKstp ); // Kst0 also as mass1
1165 isKstp = 1;
1166 }
1167 if ( g0[i] == 1 )
1168 {
1169 pro0[0] = proKstp[0];
1170 pro0[1] = proKstp[1];
1171 }
1172 else if ( g0[i] == 0 )
1173 {
1174 pro0[0] = 1;
1175 pro0[1] = 0;
1176 }
1177 if ( g1[i] == 1 )
1178 {
1179 propagatorRBW( mass2sq, mass2[i], width2[i], sk11, sKstp, spionm, rRes1, g2[i],
1180 pro1 ); // K1270(K*) as mass2
1181 }
1182 else if ( g1[i] == 0 )
1183 {
1184 pro1[0] = 1;
1185 pro1[1] = 0;
1186 }
1187 Com_Multi( pro0, pro1, pro );
1188 amp_tmp[0] = -tmp1 * pro[0];
1189 amp_tmp[1] = -tmp1 * pro[1];
1190 }
1191 else if ( modetype[i] == 14 )
1192 {
1193 double B1_omega1 = barrier( 1, somega, srhop, spionm, rRes1, mass1[i] );
1194 double B1_omega2 = barrier( 1, somega, srho0, spi0, rRes1, mass1[i] );
1195 double B1_omega3 = barrier( 1, somega, srhom, spionp, rRes1, mass1[i] );
1196 calt1( pomega, Kp, t1Ds_omega );
1197 double B1_1V11 = barrier( 1, srhop, spionp, spi0, rRes1, mrho );
1198 double B1_1V12 = barrier( 1, srho0, spionp, spionm, rRes1, mrho0 );
1199 double B1_1V13 = barrier( 1, srhom, spi0, spionm, rRes1, mrho );
1200 double B1_Ds_omega = barrier( 1, sD, somega, skaon, rD2, mDsM );
1201 for ( int w = 0; w < 4; w++ )
1202 {
1203 tt1 = pomega[w] * G[w][w];
1204 for ( int j = 0; j < 4; j++ )
1205 {
1206 tt21 = ( prhop[j] - Pim[j] ) * G[j][j];
1207 tt22 = ( prho0[j] - Pi0[j] ) * G[j][j];
1208 tt23 = ( prhom[j] - Pip[j] ) * G[j][j];
1209 for ( int k = 0; k < 4; k++ )
1210 {
1211 tt3 = Kp[k] * G[k][k];
1212 for ( int l = 0; l < 4; l++ )
1213 {
1214 temp_PDF11 += E[w][j][k][l] * tt1 * tt21 * tt3 * ( Pip[l] - Pi0[l] ) * G[l][l];
1215 temp_PDF12 += E[w][j][k][l] * tt1 * tt22 * tt3 * ( Pip[l] - Pim[l] ) * G[l][l];
1216 temp_PDF13 += E[w][j][k][l] * tt1 * tt23 * tt3 * ( Pim[l] - Pi0[l] ) * G[l][l];
1217 }
1218 }
1219 }
1220 }
1221 temp_PDF11 = temp_PDF11 * B1_Ds_omega * B1_omega1 * B1_1V11;
1222 temp_PDF12 = temp_PDF12 * B1_Ds_omega * B1_omega2 * B1_1V12;
1223 temp_PDF13 = temp_PDF13 * B1_Ds_omega * B1_omega3 * B1_1V13;
1224 double pro1V11[2], pro1V12[2], pro1V13[2];
1225 propagatorGS( mass2sq, mass2[i], width2[i], srhop, spionp, spi0, rRes1, pro1V11 );
1226 propagatorGS( mass2sq, mass2[i], width2[i], srho0, spionp, spionm, rRes1, pro1V12 );
1227 propagatorGS( mass2sq, mass2[i], width2[i], srhom, spi0, spionm, rRes1, pro1V13 );
1228 if ( g0[i] == 1 )
1229 {
1230 propagatorRBW( mass1sq, mass1[i], width1[i], somega, srhop, spionm, rRes1, 1,
1231 pro0_11 );
1232 propagatorRBW( mass1sq, mass1[i], width1[i], somega, srho0, spi0, rRes1, 1, pro0_12 );
1233 propagatorRBW( mass1sq, mass1[i], width1[i], somega, srhom, spionp, rRes1, 1,
1234 pro0_13 );
1235 }
1236 else if ( g0[i] == 0 )
1237 {
1238 pro0_11[0] = 1;
1239 pro0_11[1] = 0;
1240 pro0_12[0] = 1;
1241 pro0_12[1] = 0;
1242 pro0_13[0] = 1;
1243 pro0_13[1] = 0;
1244 }
1245 Com_Multi( pro0_11, pro1V11, pro_11 );
1246 Com_Multi( pro0_12, pro1V12, pro_12 );
1247 Com_Multi( pro0_13, pro1V13, pro_13 );
1248 amp_tmp[0] =
1249 ( temp_PDF11 * pro_11[0] + temp_PDF12 * pro_12[0] + temp_PDF13 * pro_13[0] );
1250 amp_tmp[1] =
1251 ( temp_PDF11 * pro_11[1] + temp_PDF12 * pro_12[1] + temp_PDF13 * pro_13[1] );
1252 }
1253 if ( modetype[i] == 15 )
1254 {
1255 tmp1 = 1;
1256 if ( isPipPi0_S == 0 )
1257 {
1258 proPiPi_S[0] = realpipis;
1259 proPiPi_S[1] = imagpipis;
1260 isPipPi0_S = 1;
1261 }
1262 if ( isKpPim_S == 0 )
1263 {
1264 KPiSLASS( sKstp, skaon, spi0, proKPi_S );
1265 isKpPim_S = 1;
1266 }
1267 Com_Multi( proKPi_S, proPiPi_S, amp_tmp );
1268 }
1269 else if ( modetype[i] == 16 )
1270 {
1271 if ( B_D_A1 < 0.0 ) B_D_A1 = barrier( 1, sD, sa1, skaon, rD2, mDsM );
1272 if ( B_rhop < 0.0 ) B_rhop = barrier( 1, srhop, spionp, spi0, rRes1, mass2[i] );
1273 calt1( pA1, Kp, t1D );
1274 if ( g2[i] == 0 )
1275 {
1276 for ( int w = 0; w < 4; w++ )
1277 {
1278 tt1 = t1D[w] * G[w][w];
1279 for ( int j = 0; j < 4; j++ )
1280 { temp_PDF += tt1 * ( G[w][j] - pA1[w] * pA1[j] / sa1 ) * t1rhop[j] * G[j][j]; }
1281 }
1282 tmp1 = B_rhop * B_D_A1 * temp_PDF;
1283 }
1284 else if ( g2[i] == 2 )
1285 {
1286 for ( int w = 0; w < 4; w++ )
1287 {
1288 tt1 = t1D[w] * G[w][w];
1289 for ( int j = 0; j < 4; j++ ) { temp_PDF += tt1 * t2A1[w][j] * t1rhop[j] * G[j][j]; }
1290 }
1291 if ( B_D_A2 < 0.0 ) B_D_A2 = barrier( 2, sa1, srhop, spionm, rRes1, mass1[i] );
1292 tmp1 = B_rhop * B_D_A2 * B_D_A1 * temp_PDF;
1293 }
1294 if ( isRhop == 0 )
1295 {
1296 propagatorGS( mass2sq, mass2[i], width2[i], srhop, spionp, spi0, rRes1, proRhop );
1297 isRhop = 1;
1298 }
1299 if ( g0[i] == 1 )
1300 propagatorRBW( mass1sq, mass1[i], width1[i], sa1, srhop, spionm, rRes1, g2[i], pro0 );
1301 else if ( g0[i] == 0 )
1302 {
1303 pro0[0] = 1;
1304 pro0[1] = 0;
1305 }
1306 if ( g1[i] == 1 )
1307 {
1308 pro1[0] = proRhop[0];
1309 pro1[1] = proRhop[1];
1310 }
1311 else if ( g1[i] == 0 )
1312 {
1313 pro1[0] = 1;
1314 pro1[1] = 0;
1315 }
1316 Com_Multi( pro0, pro1, pro );
1317 amp_tmp[0] = tmp1 * pro[0];
1318 amp_tmp[1] = tmp1 * pro[1];
1319 }
1320 else if ( modetype[i] == 17 )
1321 {
1322 if ( B_D_A1 < 0.0 ) B_D_A1 = barrier( 1, sD, sa1, skaon, rD2, mDsM );
1323 if ( B_rhom < 0.0 ) B_rhom = barrier( 1, srhom, spionm, spi0, rRes1, mass2[i] );
1324 calt1( pA1, Kp, t1D );
1325 if ( g2[i] == 0 )
1326 {
1327 for ( int w = 0; w < 4; w++ )
1328 {
1329 tt1 = t1D[w] * G[w][w];
1330 for ( int j = 0; j < 4; j++ )
1331 { temp_PDF += tt1 * ( G[w][j] - pA1[w] * pA1[j] / sa1 ) * t1rhom[j] * G[j][j]; }
1332 }
1333 tmp1 = B_rhom * B_D_A1 * temp_PDF;
1334 }
1335 else if ( g2[i] == 2 )
1336 {
1337 for ( int w = 0; w < 4; w++ )
1338 {
1339 tt1 = t1D[w] * G[w][w];
1340 for ( int j = 0; j < 4; j++ ) { temp_PDF += tt1 * t2A2[w][j] * t1rhom[j] * G[j][j]; }
1341 }
1342 if ( B_D_A2 < 0.0 ) B_D_A2 = barrier( 2, sa1, srhom, spionp, rRes1, mass1[i] );
1343 tmp1 = B_rhom * B_D_A2 * B_D_A1 * temp_PDF;
1344 }
1345 if ( isRhom == 0 )
1346 {
1347 propagatorGS( mass2sq, mass2[i], width2[i], srhom, spionm, spi0, rRes1, proRhom );
1348 isRhom = 1;
1349 }
1350 if ( g0[i] == 1 )
1351 propagatorRBW( mass1sq, mass1[i], width1[i], sa1, srhom, spionp, rRes1, g2[i], pro0 );
1352 else if ( g0[i] == 0 )
1353 {
1354 pro0[0] = 1;
1355 pro0[1] = 0;
1356 }
1357 if ( g1[i] == 1 )
1358 {
1359 pro1[0] = proRhom[0];
1360 pro1[1] = proRhom[1];
1361 }
1362 else if ( g1[i] == 0 )
1363 {
1364 pro1[0] = 1;
1365 pro1[1] = 0;
1366 }
1367 Com_Multi( pro0, pro1, pro );
1368 amp_tmp[0] = -tmp1 * pro[0];
1369 amp_tmp[1] = -tmp1 * pro[1];
1370 }
1371 Com_Multi( amp_tmp, cof, amp_PDF );
1372 PDF[0] += amp_PDF[0];
1373 PDF[1] += amp_PDF[1];
1374 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
1375 if ( value <= 0 ) { value = 1e-20; }
1376 Result = value;
1377 q1270 = qqq;
1378 }
1379}
mass
double mass
Definition CosmicGenerator.cxx:128
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
TComplex
complex< double > TComplex
Definition Dalitz.h:14
g1
TF1 * g1
Definition DataBase/tau_mode.c:58
mass_Pion
const double mass_Pion
Definition EvtD0ToKSpipipi0.cc:1561
w
double w
Definition EvtPhokharaDef.hh:28
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
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
EvtDsToKpPipPimPi0::clone
EvtDecayBase * clone()
Definition EvtDsToKpPipPimPi0.cc:41
EvtDsToKpPipPimPi0::decay
void decay(EvtParticle *p)
Definition EvtDsToKpPipPimPi0.cc:155
EvtDsToKpPipPimPi0::getName
void getName(std::string &name)
Definition EvtDsToKpPipPimPi0.cc:39
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
m2
double double * m2
Definition qcdloop1.h:83
t
int t()
Definition t.c:1