BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDToKSKmpippip.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// Environment:
3// This software is part of models developed at BES collaboration
4// based on the EvtGen framework. If you use all or part
5// of it, please give an appropriate acknowledgement.
6//
7// Copyright Information: See EvtGen/BesCopyright
8// Copyright (A) 2006 Ping Rong-Gang @IHEP
9//
10// Module: EvtDToKSKmpippip.cc
11// the necessary file: EvtDToKSKmpippip.hh
12//
13// Description: D+ -> KS K- pi+ pi+
14//
15// Modification history:
16//
17// Zhenxuan Li Jun. 15, 2023 Module created
18// Jiahao Wang Nov. 03, 2023 modif
19//------------------------------------------------------------------------
20#include "EvtDToKSKmpippip.hh"
21#include "../EvtGenBase/EvtComplex.hh"
22#include "../EvtGenBase/EvtDecayTable.hh"
23#include "../EvtGenBase/EvtGenKine.hh"
24#include "../EvtGenBase/EvtPDL.hh"
25#include "../EvtGenBase/EvtParticle.hh"
26#include "../EvtGenBase/EvtPatches.hh"
27#include "../EvtGenBase/EvtReport.hh"
28#include "../EvtGenBase/EvtVector4R.hh"
29#include <cmath>
30#include <iostream>
31#include <stdlib.h>
32using namespace std;
33
34EvtDToKSKmpippip::~EvtDToKSKmpippip() {}
35
36void EvtDToKSKmpippip::getName( std::string& model_name ) { model_name = "DToKSKmpippip"; }
37
38EvtDecayBase* EvtDToKSKmpippip::clone() { return new EvtDToKSKmpippip; }
39
40void EvtDToKSKmpippip::init() {
41 checkNArg( 0 );
42 checkNDaug( 4 );
43 checkSpinParent( EvtSpinType::SCALAR );
44 // cout << "Initializing EvtD0ToKSKppimpi0" << endl;
45 mass_Pion = 0.13957;
46 mass_Pion_N = 0.134977;
47 mass_Eta = 0.547862;
48 mD = 1.86966;
49 mPi = 0.13957;
50 mKa = 0.493677;
51 math_pi = 3.1415926;
52 rD2 = 25.0;
53 rRes2 = 9.0;
54 GS1 = 0.636619783;
55 GS2 = 0.01860182466;
56 GS3 = 0.1591549458;
57 GS4 = 0.00620060822;
58
59 mK1400 = 1.403;
60 GK1400 = 0.174;
61 m1510 = 1.518;
62 G1510 = 0.073;
63 meta1405 = 1.4139;
64 Geta1405 = 0.048;
65 meta1475 = 1.475;
66 Geta1475 = 0.09;
67 mf1285 = 1.2819;
68 Gf1285 = 0.0227;
69 mKstp = 0.89176;
70 GKstp = 0.0514;
71 mKst0 = 0.89555;
72 GKst0 = 0.0473;
73 ma0_980 = 0.990;
74 Ga0_980 = 0.05;
75
76 rho[0] = 1.0; // 2S
77 phi[0] = 0.0;
78
79 phi[1] = 4.1350e-02;
80 phi[2] = -8.3438e-01;
81 phi[3] = -7.9307e-01;
82 phi[4] = -2.2340e+00;
83 phi[5] = -3.5357e-01;
84 phi[6] = -3.5357e-01;
85 phi[7] = 9.1698e+00;
86 phi[8] = 2.1065e+00;
87 rho[1] = 5.5001e-01; // 1S
88 rho[2] = 8.9374e-01; // 2P
89 rho[3] = 1.8706e+00; // 16S
90 rho[4] = 1.2629e+01; // 201
91 rho[5] = 3.7866e+00; // 201
92 rho[6] = 3.7866e+00; // 101
93 rho[7] = 8.6436e-01; // 101
94 rho[8] = 2.1213e+00; // 101
95
96 modetype[0] = 2;
97 modetype[1] = 2;
98 modetype[2] = 2;
99 modetype[3] = 15;
100 modetype[4] = 1500;
101 modetype[5] = 1200;
102 modetype[6] = 1100;
103 modetype[7] = 9;
104 modetype[8] = 20;
105
106 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
107 int EE[4][4][4][4] = {
108 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
109 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, -1, 0 } },
110 { { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 1, 0, 0 } },
111 { { 0, 0, 0, 0 }, { 0, 0, 1, 0 }, { 0, -1, 0, 0 }, { 0, 0, 0, 0 } } },
112 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 1, 0 } },
113 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
114 { { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 } },
115 { { 0, 0, -1, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 } } },
116 { { { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, -1, 0, 0 } },
117 { { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 } },
118 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
119 { { 0, 1, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } },
120 { { { 0, 0, 0, 0 }, { 0, 0, -1, 0 }, { 0, 1, 0, 0 }, { 0, 0, 0, 0 } },
121 { { 0, 0, 1, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 } },
122 { { 0, -1, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
123 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } } };
124 for ( int i = 0; i < 4; i++ )
125 {
126 for ( int j = 0; j < 4; j++ )
127 {
128 G[i][j] = GG[i][j];
129 for ( int k = 0; k < 4; k++ )
130 {
131 for ( int l = 0; l < 4; l++ ) { E[i][j][k][l] = EE[i][j][k][l]; }
132 }
133 }
134 }
135}
136
137void EvtDToKSKmpippip::initProbMax() { setProbMax( 97200 ); }
138
139void EvtDToKSKmpippip::decay( EvtParticle* p ) {
140 //-----------for max value------------------
141 ///
142 ////
143 /*
144 double maxprob = 0.0;
145 for(int ir=0;ir<=60000000;ir++){
146 p->initializePhaseSpace(getNDaug(),getDaugs());
147 EvtVector4R ks = p->getDaug(0)->getP4();
148 EvtVector4R km = p->getDaug(1)->getP4();
149 EvtVector4R pip1 = p->getDaug(2)->getP4();
150 EvtVector4R pip2 = p->getDaug(3)->getP4();
151 mother_c=EvtPDL::getStdHep(p->getId());
152 int cc;
153 if(mother_c==411){
154 cc=1;
155 } else if(mother_c==-411){
156 cc=-1;
157 }
158 double value;
159 double KS[4],Km[4],Pip1[4],Pip2[4];
160
161 KS[0] = ks.get(0); Km[0] = km.get(0); Pip1[0] = pip1.get(0); Pip2[0]
162 = pip2.get(0); KS[1] = cc*ks.get(1); Km[1] = cc*km.get(1); Pip1[1] = cc*pip1.get(1);
163 Pip2[1] = cc*pip2.get(1); KS[2] = cc*ks.get(2); Km[2] = cc*km.get(2); Pip1[2] =
164 cc*pip1.get(2); Pip2[2] = cc*pip2.get(2); KS[3] = cc*ks.get(3); Km[3] = cc*km.get(3);
165 Pip1[3] = cc*pip1.get(3); Pip2[3] = cc*pip2.get(3); double mass1[9] ={ mKstp, mKstp,
166 mKstp, mf1285, m1510, m1510, m1510, meta1405, meta1475}; double mass2[9] ={
167 mKst0, mKst0, mKst0, ma0_980, ma0_980, mKst0, mKstp, ma0_980, mKstp };
168 double width1[9]={ GKstp, GKstp, GKstp, Gf1285, G1510, G1510, G1510,
169 Geta1405, Geta1475}; double width2[9]={ GKst0, GKst0, GKst0, Ga0_980, Ga0_980,
170 GKst0, GKstp, Ga0_980, GKstp }; int g0[9] ={ 1, 1, 1, 1, 1,
171 1, 1, 1, 1 }; int g1[9] ={ 1, 1, 1, 1,
172 1, 1, 1, 1, 1 }; int g2[9] ={ 0, 1, 2,
173 0, 1, 0, 0, 0, 0 };
174
175 value=calEva(KS,Km,Pip1,Pip2,mass1,mass2,width1,width2,rho,phi,g0,g1,g2,modetype,value);
176 if(value>maxprob) {
177 maxprob=value;
178 std::cout << "Max PDF = " << ir << " prob= " << value << std::endl;
179 }
180 }
181 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
182 return;
183 //-----------------------------------------------
184 */
185 p->initializePhaseSpace( getNDaug(), getDaugs() );
186 EvtVector4R ks = p->getDaug( 0 )->getP4();
187 EvtVector4R km = p->getDaug( 1 )->getP4();
188 EvtVector4R pip1 = p->getDaug( 2 )->getP4();
189 EvtVector4R pip2 = p->getDaug( 3 )->getP4();
190
191 mother_c = EvtPDL::getStdHep( p->getId() );
192 int cc;
193 if ( mother_c == 411 ) { cc = 1; }
194 else if ( mother_c == -411 ) { cc = -1; }
195
196 double KS[4], Km[4], Pip1[4], Pip2[4];
197
198 KS[0] = ks.get( 0 );
199 Km[0] = km.get( 0 );
200 Pip1[0] = pip1.get( 0 );
201 Pip2[0] = pip2.get( 0 );
202 KS[1] = cc * ks.get( 1 );
203 Km[1] = cc * km.get( 1 );
204 Pip1[1] = cc * pip1.get( 1 );
205 Pip2[1] = cc * pip2.get( 1 );
206 KS[2] = cc * ks.get( 2 );
207 Km[2] = cc * km.get( 2 );
208 Pip1[2] = cc * pip1.get( 2 );
209 Pip2[2] = cc * pip2.get( 2 );
210 KS[3] = cc * ks.get( 3 );
211 Km[3] = cc * km.get( 3 );
212 Pip1[3] = cc * pip1.get( 3 );
213 Pip2[3] = cc * pip2.get( 3 );
214
215 // KS[0] = 0.5844057460 ; Km[0] = 0.6644209662 ; Pip1[0] = 0.3998331690 ; Pip2[0] =
216 // 0.2369772689 ; KS[1] =-0.0451579928 ; Km[1] =-0.1621486835 ; Pip1[1] =
217 // 0.2182356619 ; Pip2[1] = 0.0516153672 ; KS[2] =-0.2455917598 ; Km[2] = 0.3728754865 ;
218 // Pip1[2]
219 //=-0.1884855750 ; Pip2[2] =-0.1706218636 ; KS[3] =-0.1776479997 ; Km[3]
220 //=-0.1800275423 ; Pip1[3] = 0.2392343060 ; Pip2[3] = 0.0700168961 ;
221
222 double value;
223 double mass1[9] = { mKstp, mKstp, mKstp, mf1285, m1510, m1510, m1510, meta1405, meta1475 };
224 double mass2[9] = { mKst0, mKst0, mKst0, ma0_980, ma0_980, mKst0, mKstp, ma0_980, mKstp };
225 double width1[9] = { GKstp, GKstp, GKstp, Gf1285, G1510, G1510, G1510, Geta1405, Geta1475 };
226 double width2[9] = { GKst0, GKst0, GKst0, Ga0_980, Ga0_980, GKst0, GKstp, Ga0_980, GKstp };
227 int g0[9] = { 1, 1, 1, 1, 1, 1, 1, 1, 1 };
228 int g1[9] = { 1, 1, 1, 1, 1, 1, 1, 1, 1 };
229 int g2[9] = { 0, 1, 2, 0, 1, 0, 0, 0, 0 };
230
231 value = calEva( KS, Km, Pip1, Pip2, mass1, mass2, width1, width2, rho, phi, g0, g1, g2,
232 modetype, value );
233 setProb( value );
234
235 return;
236}
237
238void EvtDToKSKmpippip::Com_Multi( double a1[2], double a2[2], double res[2] ) {
239 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
240 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
241}
242void EvtDToKSKmpippip::Com_Divide( double a1[2], double a2[2], double res[2] ) {
243 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
244 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
245 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
246}
247//------------base---------------------------------
248double EvtDToKSKmpippip::SCADot( double a1[4], double a2[4] ) {
249 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
250 return _cal;
251}
252
253double EvtDToKSKmpippip::Barrier( double mass2, int l, double sa, double sb, double sc,
254 double r2 ) {
255 double F;
256 double tmp = sa + sb - sc;
257 double q = fabs( 0.25 * tmp * tmp / sa - sb );
258 // if (q < 0) q = 1e-16;
259 double tmp2 = mass2 + sb - sc;
260 double q0 = fabs( 0.25 * tmp2 * tmp2 / mass2 - sb );
261 // if (q0 < 0) q0 = 1e-16;
262 double z = q * r2;
263 double z0 = q0 * r2;
264 if ( l == 1 ) { F = sqrt( ( 1.0 + z0 ) / ( 1.0 + z ) ); }
265 else if ( l == 2 )
266 {
267 double z2 = z * z;
268 double z02 = z0 * z0;
269 F = sqrt( ( 9.0 + 3.0 * z0 + z02 ) / ( 9.0 + 3.0 * z + z2 ) );
270 }
271 else { F = 1.0; }
272 return F;
273}
274
275double EvtDToKSKmpippip::barrier( int l, double sa, double sb, double sc, double r2 ) {
276 double F;
277 double tmp = sa + sb - sc;
278 double q = 0.25 * tmp * tmp / sa - sb;
279 if ( q < 0 ) q = 1e-16;
280 double z = q * r2;
281 if ( l == 1 ) { F = sqrt( 2.0 * z / ( 1.0 + z ) ); }
282 else if ( l == 2 )
283 {
284 double z2 = z * z;
285 F = sqrt( 13.0 * z2 / ( 9.0 + 3.0 * z + z2 ) );
286 }
287 else { F = 1.0; }
288 return F;
289}
290//------------------spin-------------------------------------------
291void EvtDToKSKmpippip::calt1( double daug1[4], double daug2[4], double t1[4] ) {
292 double p, pq, tmp;
293 double pa[4], qa[4];
294 for ( int i = 0; i < 4; i++ )
295 {
296 pa[i] = daug1[i] + daug2[i];
297 qa[i] = daug1[i] - daug2[i];
298 }
299 p = SCADot( pa, pa );
300 pq = SCADot( pa, qa );
301 tmp = pq / p;
302 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
303}
304void EvtDToKSKmpippip::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
305 double p, r;
306 double pa[4], t1[4];
307 calt1( daug1, daug2, t1 );
308 r = SCADot( t1, t1 ) / 3.0;
309 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
310 p = SCADot( pa, pa );
311 for ( int i = 0; i < 4; i++ )
312 {
313 for ( int j = 0; j < 4; j++ )
314 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
315 }
316}
317//-------------------prop--------------------------------------------
318void EvtDToKSKmpippip::Flatte_rhoab( double sa, double sb, double sc, double rho[2] ) {
319 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
320 if ( q > 0 )
321 {
322 rho[0] = 2 * sqrt( q / sa );
323 rho[1] = 0;
324 }
325 else if ( q < 0 )
326 {
327 rho[0] = 0;
328 rho[1] = 2 * sqrt( -q / sa );
329 }
330}
331void EvtDToKSKmpippip::propagator980( double mass2, double mass, double sx, double* sb,
332 double* sc, double prop[2] ) {
333 double unit[2] = { 1.0 };
334 double ci[2] = { 0, 1 };
335 double rho1[2];
336 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
337 double rho2[2];
338 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
339
340 double gK_f980 = 0.69465, gPi_f980 = 0.165;
341 double tmp1[2] = { gK_f980, 0 };
342 double tmp11[2];
343 double tmp2[2] = { gPi_f980, 0 };
344 double tmp22[2];
345 Com_Multi( tmp1, rho1, tmp11 );
346 Com_Multi( tmp2, rho2, tmp22 );
347 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
348 double tmp31[2];
349 Com_Multi( tmp3, ci, tmp31 );
350 double tmp4[2] = { mass2 - sx - tmp31[0], -1.0 * tmp31[1] };
351 Com_Divide( unit, tmp4, prop );
352}
353void EvtDToKSKmpippip::propagatora0980( double mass2, double mass, double sx, double* sb,
354 double* sc, double prop[2] ) {
355 double unit[2] = { 1.0 };
356 double ci[2] = { 0, 1 };
357 double rho1[2];
358 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
359 double rho2[2];
360 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
361
362 double gKK_a980 = 0.892 * 0.341, gPiEta_a980 = 0.341;
363 double tmp1[2] = { gKK_a980, 0 };
364 double tmp11[2];
365 double tmp2[2] = { gPiEta_a980, 0 };
366 double tmp22[2];
367 Com_Multi( tmp1, rho1, tmp11 );
368 Com_Multi( tmp2, rho2, tmp22 );
369 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
370 double tmp31[2];
371 Com_Multi( tmp3, ci, tmp31 );
372 double tmp4[2] = { mass2 - sx - tmp31[0], -1.0 * tmp31[1] };
373 Com_Divide( unit, tmp4, prop );
374}
375//--------------------------------------------------------------------------------------------------------
376void EvtDToKSKmpippip::propagator( double mass2, double mass, double width, double sx,
377 double prop[2] ) {
378 double a[2], b[2];
379 a[0] = 1;
380 a[1] = 0;
381 b[0] = mass2 - sx;
382 b[1] = -mass * width;
383 Com_Divide( a, b, prop );
384}
385double EvtDToKSKmpippip::wid( double mass2, double mass, double sa, double sb, double sc,
386 double r2, int l ) {
387 double widm = 0.;
388 double m = sqrt( sa );
389 double tmp = sb - sc;
390 double tmp1 = sa + tmp;
391 double q = 0.25 * tmp1 * tmp1 / sa - sb;
392 if ( q < 0 ) q = 1e-16;
393 double tmp2 = mass2 + tmp;
394 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
395 if ( q0 < 0 ) q0 = 1e-16;
396 double z = q * r2;
397 double z0 = q0 * r2;
398 double t = q / q0;
399 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
400 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
401 else if ( l == 2 )
402 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
403 return widm;
404}
405double EvtDToKSKmpippip::widl1( double mass2, double mass, double sa, double sb, double sc,
406 double r2 ) {
407 double widm = 0.;
408 double m = sqrt( sa );
409 double tmp = sb - sc;
410 double tmp1 = sa + tmp;
411 double q = 0.25 * tmp1 * tmp1 / sa - sb;
412 if ( q < 0 ) q = 1e-16;
413 double tmp2 = mass2 + tmp;
414 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
415 if ( q0 < 0 ) q0 = 1e-16;
416 double z = q * r2;
417 double z0 = q0 * r2;
418 double F = ( 1 + z0 ) / ( 1 + z );
419 double t = q / q0;
420 widm = t * sqrt( t ) * mass / m * F;
421 return widm;
422}
423void EvtDToKSKmpippip::propagatorRBW( double mass2, double mass, double width, double sa,
424 double sb, double sc, double r2, int l,
425 double prop[2] ) {
426 double a[2], b[2];
427 a[0] = 1;
428 a[1] = 0;
429 b[0] = mass2 - sa;
430 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
431 Com_Divide( a, b, prop );
432}
433void EvtDToKSKmpippip::propagatorRBWl1( double mass2, double mass, double width, double sa,
434 double sb, double sc, double r2, double prop[2] ) {
435 double a[2], b[2];
436 a[0] = 1;
437 a[1] = 0;
438 b[0] = mass2 - sa;
439 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
440 Com_Divide( a, b, prop );
441}
442//------------GS---used by rho----------------------------
443void EvtDToKSKmpippip::propagatorGS( double mass2, double mass, double width, double sa,
444 double sb, double sc, double r2, double prop[2] ) {
445 double a[2], b[2];
446 double tmp = sb - sc;
447 double tmp1 = sa + tmp;
448 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
449 if ( q2 < 0 ) q2 = 1e-16;
450
451 double tmp2 = mass2 + tmp;
452 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
453 if ( q02 < 0 ) q02 = 1e-16;
454
455 double q = sqrt( q2 );
456 double q0 = sqrt( q02 );
457 double m = sqrt( sa );
458 double q03 = q0 * q02;
459 double tmp3 = log( mass + 2 * q0 ) + 1.2760418309; // log(mass_2Pion) = 1.2760418309;
460
461 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2760418309 );
462 double h0 = GS1 * q0 / mass * tmp3;
463 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
464 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
465 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
466
467 a[0] = 1.0 + d * width / mass;
468 a[1] = 0.0;
469 b[0] = mass2 - sa + width * f;
470 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
471 Com_Divide( a, b, prop );
472}
473void EvtDToKSKmpippip::rhoab( double sa, double sb, double sc, double res[2] ) {
474 double tmp = sa + sb - sc;
475 double q = 0.25 * tmp * tmp / sa - sb;
476 if ( q >= 0 )
477 {
478 res[0] = 2.0 * sqrt( q / sa );
479 res[1] = 0.0;
480 }
481 else
482 {
483 res[0] = 0.0;
484 res[1] = 2.0 * sqrt( -q / sa );
485 }
486}
487void EvtDToKSKmpippip::rho4Pi( double sa, double res[2] ) {
488 double temp = 1.0 - 0.3116765584 / sa; // 0.3116765584=0.13957*0.13957*16
489 if ( temp >= 0 )
490 {
491 res[0] = sqrt( temp ) / ( 1.0 + exp( 9.8 - 3.5 * sa ) );
492 res[1] = 0.0;
493 }
494 else
495 {
496 res[0] = 0.0;
497 res[1] = sqrt( -temp ) / ( 1.0 + exp( 9.8 - 3.5 * sa ) );
498 }
499}
500void EvtDToKSKmpippip::propagatorsigma500( double sa, double sb, double sc,
501 double prop[2] ) { // for pipi_s wave
502 double f = 0.5843 + 1.6663 * sa;
503 const double M = 0.9264;
504 const double mass2 = 0.85821696; // M*M
505 const double mpi2d2 = 0.00973989245;
506 double g1 = f * ( sa - mpi2d2 ) / ( mass2 - mpi2d2 ) * exp( ( mass2 - sa ) / 1.082 );
507 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
508 rhoab( sa, sb, sc, rho1s );
509 rhoab( mass2, sb, sc, rho1M );
510 rho4Pi( sa, rho2s );
511 rho4Pi( mass2, rho2M );
512 Com_Divide( rho1s, rho1M, rho1 );
513 Com_Divide( rho2s, rho2M, rho2 );
514 double a[2], b[2];
515 a[0] = 1.0;
516 a[1] = 0.0;
517 b[0] = mass2 - sa + M * ( g1 * rho1[1] + 0.0024 * rho2[1] );
518 b[1] = -M * ( g1 * rho1[0] + 0.0024 * rho2[0] );
519 Com_Divide( a, b, prop );
520}
521void EvtDToKSKmpippip::KPiSLASS( double sa, double sb, double sc, double prop[2] ) {
522 const double m1430 = 1.441;
523 const double sa0 = 2.076481; // m1430*m1430;
524 const double w1430 = 0.193;
525 const double Lass1 = 0.25 / sa0;
526 double tmp = sb - sc;
527 double tmp1 = sa0 + tmp;
528 double q0 = fabs( Lass1 * tmp1 * tmp1 - sb );
529 // if(q0<0) q0 = 1e-16;
530 double tmp2 = sa + tmp;
531 double qs = fabs( 0.25 * tmp2 * tmp2 / sa - sb );
532 double q = sqrt( qs );
533 double width = w1430 * q * m1430 / sqrt( sa * q0 );
534 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
535 if ( temp_R < 0 ) temp_R += math_pi;
536 double deltaR = -109.7 * math_pi / 180.0 + temp_R;
537 double temp_F =
538 atan( 0.226 * q / ( 2.0 - 3.8194 * qs ) ); // 2.0*0.113 = 0.226; -33.8*0.113 = -3.8194
539 if ( temp_F < 0 ) temp_F += math_pi;
540 double deltaF = 0.1 * math_pi / 180.0 + temp_F;
541 double deltaS = deltaR + 2.0 * deltaF;
542 double t1 = 0.96 * sin( deltaF );
543 double t2 = sin( deltaR );
544 double CF[2], CS[2];
545 CF[0] = cos( deltaF );
546 CF[1] = sin( deltaF );
547 CS[0] = cos( deltaS );
548 CS[1] = sin( deltaS );
549 prop[0] = t1 * CF[0] + t2 * CS[0];
550 prop[1] = t1 * CF[1] + t2 * CS[1];
551}
552
553//---------------------------------------------------------------------------------------------------------------------
554double EvtDToKSKmpippip::calEva( double* Km, double* Kp, double* Pip, double* Pi0,
555 double* mass1, double* mass2, double* width1, double* width2,
556 double* amp, double* phase, int* g0, int* g1, int* g2,
557 int* modetype, double Result )
558// r0,r1 just send
559// propagatorGS_Omg(mass1sq,mass1[i],width1[i],srho1,spion1,spi0,rRes2,amp_Omg,phs_Omg,proRho1);
560
561{
562 double cof[2], amp_tmp1[2], amp_tmp2[2], amp_tmp[2], amp_PDF[2], PDF[2];
563 double pKstr0[4], pKstrp[4], pKstrm[4], prhop[4], pphi[4], pK10[4], pK11[4], pK12[4],
564 pK13[4], pK14[4]; // pD[4];
565 double pD[4], pKsPip1[4], pKmPip2[4], pKsPip2[4], pKmPip1[4], pKsKmPip1[4], pKsKmPip2[4],
566 pKsKm[4], pKmPipPip[4], pKsPipPip[4];
567
568 for ( int i = 0; i != 4; i++ )
569 {
570 pD[i] = Km[i] + Pip[i] + Kp[i] + Pi0[i];
571 pKsPip1[i] = Km[i] + Pip[i];
572 pKsPip2[i] = Km[i] + Pi0[i];
573 pKmPip1[i] = Kp[i] + Pip[i];
574 pKmPip2[i] = Kp[i] + Pi0[i];
575 pKsKm[i] = Km[i] + Kp[i];
576 pKsKmPip1[i] = pKsPip1[i] + Kp[i];
577 pKsKmPip2[i] = pKsPip2[i] + Kp[i];
578 pKmPipPip[i] = Kp[i] + Pip[i] + Pi0[i];
579 pKsPipPip[i] = Km[i] + Pip[i] + Pi0[i];
580
581 pphi[i] = Km[i] + Kp[i];
582 prhop[i] = Pip[i] + Pi0[i];
583 pKstr0[i] = Km[i] + Pip[i];
584 pKstrp[i] = Kp[i] + Pi0[i];
585 pKstrm[i] = Km[i] + Pi0[i];
586 pK10[i] = pKstrm[i] + Kp[i];
587 pK11[i] = pKstrm[i] + Pip[i];
588 pK12[i] = pKstrp[i] + Km[i];
589 pK13[i] = pKstr0[i] + Pi0[i]; // can delete, no use, equal to pK11
590 pK14[i] = pphi[i] + Pi0[i];
591 }
592 double spi0, spionp, skaon1, skaon2, sphi, srhop, skstr0, skstrp, skstrm, sk10, sk11, sk12,
593 sk13, sk14; //,sD;
594 double sPip2, sPip1, sKs, sKm, sKsPip1, sKmPip2, sD, sKsPip2, sKmPip1, sKsKmPip1, sKsKmPip2,
595 sKsKm, sKmPipPip, sKsPipPip;
596
597 sKs = SCADot( Km, Km );
598 sKm = SCADot( Kp, Kp );
599 sPip1 = SCADot( Pip, Pip );
600 sPip2 = SCADot( Pi0, Pi0 );
601 sKsPip1 = SCADot( pKsPip1, pKsPip1 );
602 sKsPip2 = SCADot( pKsPip2, pKsPip2 );
603 sKmPip1 = SCADot( pKmPip1, pKmPip1 );
604 sKmPip2 = SCADot( pKmPip2, pKmPip2 );
605 sKsKmPip1 = SCADot( pKsKmPip1, pKsKmPip1 );
606 sKsKmPip2 = SCADot( pKsKmPip2, pKsKmPip2 );
607 sKmPipPip = SCADot( pKmPipPip, pKmPipPip );
608 sKsPipPip = SCADot( pKsPipPip, pKsPipPip );
609 sKsKm = SCADot( pKsKm, pKsKm );
610
611 skaon1 = SCADot( Km, Km );
612 skaon2 = SCADot( Kp, Kp );
613 spionp = SCADot( Pip, Pip );
614 spi0 = SCADot( Pi0, Pi0 );
615
616 sphi = SCADot( pphi, pphi );
617 srhop = SCADot( prhop, prhop );
618 skstr0 = SCADot( pKstr0, pKstr0 );
619 skstrp = SCADot( pKstrp, pKstrp );
620 skstrm = SCADot( pKstrm, pKstrm );
621 sk10 = SCADot( pK10, pK10 );
622 sk11 = SCADot( pK11, pK11 );
623 sk12 = SCADot( pK12, pK12 );
624 // sk13 = SCADot(pK13,pK13);
625 sk14 = SCADot( pK14, pK14 );
626 sD = SCADot( pD, pD );
627 //-------------------------------------------------------------------------------------------------
628 double t1phi[4], t1A1[4], t1A2[4], t1rhop[4], t1kstr1[4], t1kstr2[4], t1kstr3[4], t1kstr4[4],
629 t2k11[4][4], t2k12[4][4], t2k13[4][4], t2k14[4][4], t2k21[4][4];
630 double t1KsPip1[4], t1KmPip2[4], t1KsPip2[4], t1KmPip1[4], t2KsKmPip1_kspip[4][4],
631 t2KsKmPip2_kspip[4][4], t2KsKmPip1_kmpip[4][4], t2KsKmPip2_kmpip[4][4],
632 t1KsKmPip1_kspip[4], t1KsKmPip2_kspip[4], t1KsKmPip1_kmpip[4], t1KsKmPip2_kmpip[4],
633 t1KsKmPip1_ksk[4], t1KsKmPip2_ksk[4], t2KmPipPip_kmpip[4][4], t2KmPipPip_kmpip2[4][4],
634 t2KsPipPip_kspip[4][4], t2KsPipPip_kspip2[4][4], t1KmPipPip_kmpip[4],
635 t1KmPipPip_kmpip2[4], t1KsPipPip_kspip[4], t1KsPipPip_kspip2[4];
636 calt1( Km, Pip, t1KsPip1 );
637 calt1( Km, Pi0, t1KsPip2 );
638 calt1( Kp, Pip, t1KmPip1 );
639 calt1( Kp, Pi0, t1KmPip2 );
640 calt1( pKsPip1, Kp, t1KsKmPip1_kspip );
641 calt1( pKsPip2, Kp, t1KsKmPip2_kspip );
642 calt1( pKmPip1, Km, t1KsKmPip1_kmpip );
643 calt1( pKmPip2, Km, t1KsKmPip2_kmpip );
644
645 calt1( pKmPip1, Pi0, t1KmPipPip_kmpip );
646 calt1( pKmPip2, Pip, t1KmPipPip_kmpip2 );
647 calt1( pKsPip1, Pi0, t1KsPipPip_kspip );
648 calt1( pKsPip2, Pip, t1KsPipPip_kspip2 );
649 calt1( pKsKm, Pip, t1KsKmPip1_ksk );
650 calt1( pKsKm, Pi0, t1KsKmPip2_ksk );
651
652 calt2( pKsPip1, Kp, t2KsKmPip1_kspip );
653 calt2( pKsPip2, Kp, t2KsKmPip2_kspip );
654 calt2( pKmPip1, Km, t2KsKmPip1_kmpip );
655 calt2( pKmPip2, Km, t2KsKmPip2_kmpip );
656 calt2( pKmPip2, Pip, t2KmPipPip_kmpip );
657 calt2( pKmPip1, Pi0, t2KmPipPip_kmpip2 );
658 calt2( pKsPip2, Pip, t2KsPipPip_kspip );
659 calt2( pKsPip1, Pi0, t2KsPipPip_kspip2 );
660
661 /* calt1(Km,Kp,t1phi);
662 calt1(Pip,Pi0,t1rhop);
663 calt1(Km,Pip,t1kstr1);
664 calt1(Kp,Pi0,t1kstr2);
665 calt1(Km,Pi0,t1kstr3);
666 calt1(Kp,Pi0,t1kstr4);
667 calt1(pKstr0,Pi0,t1A1);
668 calt1(pKstrm,Pip,t1A2);
669
670 calt2(pKstr0,Pi0,t2k11);
671 calt2(pKstrm,Pip,t2k12);
672 calt2(pKstrm,Kp,t2k13);
673 calt2(pKstrp,Km,t2k14);
674 calt2(prhop,Km,t2k21);
675 */
676 double B_phi = -1.0, B_rho = -1.0, B_rho_2 = -1.0, B_phirho1 = -1.0, B_phirho2 = -1.0,
677 B_phi_2 = -1.0;
678 double B_Kstp = -1.0, B_Kst0 = -1.0, B_KsKs1 = -1.0, B_KsKs2 = -1.0, B_Kstp_2 = -1.0,
679 B_Kst0_2 = -1.0, B_DA2 = -1.0;
680 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,
681 B_A2_1 = -1.0;
682 double B_D_Arho = -1.0, B_Arho = -1.0, B_K11 = -1.0, B_D11 = -1.0, B_Arho1 = -1.0,
683 B_K14_phi = -1.0, B_A10 = -1.0;
684 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,
685 B_K14_km = -1.0, B_K14_kp = -1.0;
686 double B_K_Kst0pi0 = -1.0, B_D_KK = -1.0, B_K_Kstmpip = -1.0, B_D_K11_2 = -1.0,
687 B_K11_21 = -1.0, B_K11_22 = -1.0, B_A12 = -1.0;
688
689 PDF[0] = 0;
690 PDF[1] = 0;
691 double mass1sq, mass2sq, tt1, tt2, tmp2;
692 double temp_PDF, tmp1, temp[2], tmp3, temp_PDF1;
693 double pro[2], pro0[2], pro1[2], pro2[2], pro3[2], pro4[2];
694 double pro1b[2], pro2b[2], pro_ksk[2];
695 double pro_kspi1[2], pro_kspi2[2], pro_kpi1[2], pro_kpi2[2], pro_kpi2sw[2], pro_kpi1sw[2],
696 pro_kspi2sw[2], pro_kspi1sw[2];
697 double pro_kspi1_0[2], pro_kspi2_0[2], pro_kpi1_0[2], pro_kpi2_0[2];
698 double ispro_kspi1 = 0, ispro_kspi2 = 0, ispro_kpi1 = 0, ispro_kpi2 = 0, ispro_kpi2sw = 0,
699 ispro_kpi1sw = 0, ispro_kspi2sw = 0, ispro_kspi1sw = 0, ispro_ksk = 0;
700 double barr_kskpi1_ksk = -1.0, barr_kskpi2_ksk = -1.0, barr_kskpi1_kspi = -1.0,
701 barr_kskpi2_kspi = -1.0, barr_kskpi1_kpi = -1.0, barr_kskpi2_kpi = -1.0,
702 barr_kspi1 = -1.0, barr_kspi2 = -1.0, barr_kpi1 = -1.0, barr_kpi2 = -1.0,
703 barr_kspi1_kpi2 = -1.0, barr_kspi2_kpi1 = -1.0, barr_ds_kskpi1 = -1.0,
704 barr_ds_kskpi2 = -1.0;
705 double t1A[4], t1V[4], t1D[4], t2D[4][4], B[3];
706 double sKs2[2] = { sKs, mass_Pion * mass_Pion };
707 double sKm2[2] = { sKm, mass_Eta * mass_Eta };
708 double pS[4], pV[4], pA[4];
709 double flag[3];
710
711 double sa[3], sb[3], sc[3];
712 int isKstp = 0.0, isKst0 = 0.0, isKstm = 0.0, isRho = 0.0, isf0 = 0.0, isKmPip_S = 0.0,
713 isKpPi0_S = 0.0, isKmPi0_S = 0.0, isPiPi_S = 0.0, isA0980 = 0.0;
714 double proRho[2], proKstp[2], proKstm[2], proKst0[2], proPiPi_S[2], proKmPip_S[2],
715 proKpPi0_S[2], proKmPi0_S[2], proA0980[2], prof0[2];
716 for ( int i = 0; i < 9; i++ )
717 {
718 flag[0] = g0[i];
719 flag[1] = g1[i];
720 flag[2] = g2[i];
721
722 temp_PDF = 0;
723 temp_PDF1 = 0;
724 cof[0] = amp[i] * cos( phase[i] );
725 cof[1] = amp[i] * sin( phase[i] );
726 mass1sq = mass1[i] * mass1[i]; // avoid parameters input reversed
727 mass2sq = mass2[i] * mass2[i];
728
729 if ( modetype[i] == 2 )
730 {
731 propagatorRBW( mass1sq, mass1[i], width1[i], sKsPip1, sKs, sPip1, rRes2, 1, pro_kspi1 );
732 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip2, sKm, sPip2, rRes2, 1, pro_kpi2 );
733 Com_Multi( pro_kspi1, pro_kpi2, pro1b );
734
735 if ( barr_kspi1 < 0.0 ) barr_kspi1 = barrier( 1, sKsPip1, sKs, sPip1, rRes2 );
736 if ( barr_kpi2 < 0.0 ) barr_kpi2 = barrier( 1, sKmPip2, sKm, sPip2, rRes2 );
737 if ( barr_kspi1_kpi2 < 0.0 ) barr_kspi1_kpi2 = barrier( 1, sD, sKsPip1, sKmPip2, rD2 );
738 if ( g2[i] == 0 )
739 {
740 for ( int a = 0; a < 4; a++ ) { temp_PDF += G[a][a] * t1KsPip1[a] * t1KmPip2[a]; }
741 tmp1 = barr_kspi1 * barr_kpi2 * temp_PDF;
742 }
743 if ( g2[i] == 1 )
744 {
745 calt1( pKsPip1, pKmPip2, t1D );
746 for ( int a = 0; a < 4; a++ )
747 {
748 for ( int j = 0; j < 4; j++ )
749 {
750 tt1 = G[a][a] * G[j][j] * pD[a] * t1D[j];
751 for ( int k = 0; k < 4; k++ )
752 {
753 tt2 = t1KsPip1[k] * G[k][k];
754 for ( int l = 0; l < 4; l++ )
755 { temp_PDF += E[a][j][k][l] * t1KmPip2[l] * G[l][l] * tt1 * tt2; }
756 }
757 }
758 }
759 tmp1 = barr_kspi1 * barr_kpi2 * barr_kspi1_kpi2 * temp_PDF;
760 }
761 if ( g2[i] == 2 )
762 {
763 calt2( pKsPip1, pKmPip2, t2D );
764 for ( int a = 0; a < 4; a++ )
765 {
766 for ( int j = 0; j < 4; j++ )
767 { temp_PDF += t2D[a][j] * t1KsPip1[a] * t1KmPip2[j] * G[a][a] * G[j][j]; }
768 }
769 B[2] = barrier( 2, sD, sKsPip1, sKmPip2, rD2 );
770 tmp1 = barr_kspi1 * barr_kpi2 * B[2] * temp_PDF;
771 }
772 temp_PDF = 0;
773 propagatorRBW( mass1sq, mass1[i], width1[i], sKsPip2, sKs, sPip2, rRes2, 1, pro_kspi2 );
774 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip1, sKm, sPip1, rRes2, 1, pro_kpi1 );
775 Com_Multi( pro_kspi2, pro_kpi1, pro2b );
776
777 if ( barr_kspi2 < 0.0 ) barr_kspi2 = barrier( 1, sKsPip2, sKs, sPip2, rRes2 );
778 if ( barr_kpi1 < 0.0 ) barr_kpi1 = barrier( 1, sKmPip1, sKm, sPip1, rRes2 );
779 if ( barr_kspi2_kpi1 < 0.0 ) barr_kspi2_kpi1 = barrier( 1, sD, sKsPip2, sKmPip1, rD2 );
780 if ( g2[i] == 0 )
781 {
782 for ( int a = 0; a < 4; a++ ) { temp_PDF += G[a][a] * t1KsPip2[a] * t1KmPip1[a]; }
783 tmp2 = barr_kspi2 * barr_kpi1 * temp_PDF;
784 }
785 if ( g2[i] == 1 )
786 {
787 calt1( pKsPip2, pKmPip1, t1D );
788 for ( int a = 0; a < 4; a++ )
789 {
790 for ( int j = 0; j < 4; j++ )
791 {
792 tt1 = G[a][a] * G[j][j] * pD[a] * t1D[j];
793 for ( int k = 0; k < 4; k++ )
794 {
795 tt2 = t1KsPip2[k] * G[k][k];
796 for ( int l = 0; l < 4; l++ )
797 { temp_PDF += E[a][j][k][l] * t1KmPip1[l] * G[l][l] * tt1 * tt2; }
798 }
799 }
800 }
801 tmp2 = barr_kspi2 * barr_kpi1 * barr_kspi2_kpi1 * temp_PDF;
802 }
803 if ( g2[i] == 2 )
804 {
805 calt2( pKsPip2, pKmPip1, t2D );
806 for ( int a = 0; a < 4; a++ )
807 {
808 for ( int j = 0; j < 4; j++ )
809 { temp_PDF += t2D[a][j] * t1KsPip2[a] * t1KmPip1[j] * G[a][a] * G[j][j]; }
810 }
811 B[2] = barrier( 2, sD, sKsPip2, sKmPip1, rD2 );
812 tmp2 = barr_kspi2 * barr_kpi1 * B[2] * temp_PDF;
813 }
814 }
815
816 // text
817 // K*1410 K*
818 if ( modetype[i] == 3 )
819 {
820 propagatorRBW( mass1sq, mass1[i], width1[i], sKsPip1, sKs, sPip1, rRes2, 1, pro_kspi1 );
821 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip2, sKm, sPip2, rRes2, 1, pro_kpi2 );
822 Com_Multi( pro_kspi1, pro_kpi2, pro1b );
823
824 if ( barr_kspi1 < 0.0 ) barr_kspi1 = barrier( 1, sKsPip1, sKs, sPip1, rRes2 );
825 if ( barr_kpi2 < 0.0 ) barr_kpi2 = barrier( 1, sKmPip2, sKm, sPip2, rRes2 );
826 if ( barr_kspi1_kpi2 < 0.0 ) barr_kspi1_kpi2 = barrier( 1, sD, sKsPip1, sKmPip2, rD2 );
827 if ( g2[i] == 0 )
828 {
829 for ( int a = 0; a < 4; a++ ) { temp_PDF += G[a][a] * t1KsPip1[a] * t1KmPip2[a]; }
830 tmp1 = barr_kspi1 * barr_kpi2 * temp_PDF;
831 }
832 if ( g2[i] == 1 )
833 {
834 calt1( pKsPip1, pKmPip2, t1D );
835 for ( int a = 0; a < 4; a++ )
836 {
837 for ( int j = 0; j < 4; j++ )
838 {
839 tt1 = G[a][a] * G[j][j] * pD[a] * t1D[j];
840 for ( int k = 0; k < 4; k++ )
841 {
842 tt2 = t1KsPip1[k] * G[k][k];
843 for ( int l = 0; l < 4; l++ )
844 { temp_PDF += E[a][j][k][l] * t1KmPip2[l] * G[l][l] * tt1 * tt2; }
845 }
846 }
847 }
848 tmp1 = barr_kspi1 * barr_kpi2 * barr_kspi1_kpi2 * temp_PDF;
849 }
850 if ( g2[i] == 2 )
851 {
852 calt2( pKsPip1, pKmPip2, t2D );
853 for ( int a = 0; a < 4; a++ )
854 {
855 for ( int j = 0; j < 4; j++ )
856 { temp_PDF += t2D[a][j] * t1KsPip1[a] * t1KmPip2[j] * G[a][a] * G[j][j]; }
857 }
858 B[2] = barrier( 2, sD, sKsPip1, sKmPip2, rD2 );
859 tmp1 = barr_kspi1 * barr_kpi2 * B[2] * temp_PDF;
860 }
861
862 temp_PDF = 0;
863 propagatorRBW( mass1sq, mass1[i], width1[i], sKsPip2, sKs, sPip2, rRes2, 1, pro_kspi2 );
864 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip1, sKm, sPip1, rRes2, 1, pro_kpi1 );
865 Com_Multi( pro_kspi2, pro_kpi1, pro2b );
866
867 if ( barr_kspi2 < 0.0 ) barr_kspi2 = barrier( 1, sKsPip2, sKs, sPip2, rRes2 );
868 if ( barr_kpi1 < 0.0 ) barr_kpi1 = barrier( 1, sKmPip1, sKm, sPip1, rRes2 );
869 if ( barr_kspi2_kpi1 < 0.0 ) barr_kspi2_kpi1 = barrier( 1, sD, sKsPip2, sKmPip1, rD2 );
870 if ( g2[i] == 0 )
871 {
872 for ( int a = 0; a < 4; a++ ) { temp_PDF += G[a][a] * t1KsPip2[a] * t1KmPip1[a]; }
873 tmp2 = barr_kspi2 * barr_kpi1 * temp_PDF;
874 }
875 if ( g2[i] == 1 )
876 {
877 calt1( pKsPip2, pKmPip1, t1D );
878 for ( int a = 0; a < 4; a++ )
879 {
880 for ( int j = 0; j < 4; j++ )
881 {
882 tt1 = G[a][a] * G[j][j] * pD[a] * t1D[j];
883 for ( int k = 0; k < 4; k++ )
884 {
885 tt2 = t1KsPip2[k] * G[k][k];
886 for ( int l = 0; l < 4; l++ )
887 { temp_PDF += E[a][j][k][l] * t1KmPip1[l] * G[l][l] * tt1 * tt2; }
888 }
889 }
890 }
891 tmp2 = barr_kspi2 * barr_kpi1 * barr_kspi2_kpi1 * temp_PDF;
892 }
893 if ( g2[i] == 2 )
894 {
895 calt2( pKsPip2, pKmPip1, t2D );
896 for ( int a = 0; a < 4; a++ )
897 {
898 for ( int j = 0; j < 4; j++ )
899 { temp_PDF += t2D[a][j] * t1KsPip2[a] * t1KmPip1[j] * G[a][a] * G[j][j]; }
900 }
901 B[2] = barrier( 2, sD, sKsPip2, sKmPip1, rD2 );
902 tmp2 = barr_kspi2 * barr_kpi1 * B[2] * temp_PDF;
903 }
904 }
905
906 /// Kstp Kmi_S_wave
907 if ( modetype[i] == 21 )
908 {
909 if ( ispro_kspi1 == 0 )
910 {
911 propagatorRBW( mass1sq, mass1[i], width1[i], sKsPip1, sKs, sPip1, rRes2, 1,
912 pro_kspi1 );
913 ispro_kspi1 = 1;
914 }
915 if ( ispro_kpi2sw == 0 )
916 {
917 KPiSLASS( sKmPip2, sKm, sPip2, pro_kpi2sw );
918 ispro_kpi2sw = 1;
919 }
920 Com_Multi( pro_kspi1, pro_kpi2sw, pro1b );
921 calt1( pKsPip1, pKmPip2, t1D );
922
923 if ( barr_kspi1 < 0.0 ) barr_kspi1 = barrier( 1, sKsPip1, sKs, sPip1, rRes2 );
924 if ( barr_kspi1_kpi2 < 0.0 ) barr_kspi1_kpi2 = barrier( 1, sD, sKsPip1, sKmPip2, rD2 );
925 for ( int a = 0; a < 4; a++ ) { temp_PDF += t1D[a] * t1KsPip1[a] * G[a][a]; }
926 tmp1 = barr_kspi1 * barr_kspi1_kpi2 * temp_PDF;
927
928 temp_PDF = 0;
929 if ( ispro_kspi2 == 0 )
930 {
931 propagatorRBW( mass1sq, mass1[i], width1[i], sKsPip2, sKs, sPip2, rRes2, 1,
932 pro_kspi2 );
933 ispro_kspi2 = 1;
934 }
935 if ( ispro_kpi1sw == 0 )
936 {
937 KPiSLASS( sKmPip1, sKm, sPip1, pro_kpi1sw );
938 ispro_kpi1sw = 1;
939 }
940 Com_Multi( pro_kspi2, pro_kpi1sw, pro2b );
941
942 calt1( pKsPip2, pKmPip1, t1D );
943 if ( barr_kspi2 < 0.0 ) barr_kspi2 = barrier( 1, sKsPip2, sKs, sPip2, rRes2 );
944 if ( barr_kspi2_kpi1 < 0.0 ) barr_kspi2_kpi1 = barrier( 1, sD, sKsPip2, sKmPip1, rD2 );
945 for ( int a = 0; a < 4; a++ ) { temp_PDF += t1D[a] * t1KsPip2[a] * G[a][a]; }
946 tmp2 = barr_kspi2 * barr_kspi2_kpi1 * temp_PDF;
947 }
948
949 /// Kst0 Kspi_S_wave
950 if ( modetype[i] == 22 )
951 {
952 if ( ispro_kpi2 == 0 )
953 {
954 propagatorRBW( mass1sq, mass1[i], width1[i], sKmPip2, sKm, sPip2, rRes2, 1, pro_kpi2 );
955 ispro_kpi2 = 1;
956 }
957 if ( ispro_kspi1sw == 0 )
958 {
959 KPiSLASS( sKsPip1, sKs, sPip1, pro_kspi1sw );
960 ispro_kspi1sw = 1;
961 }
962 Com_Multi( pro_kpi2, pro_kspi1sw, pro1b );
963
964 calt1( pKsPip1, pKmPip2, t1D );
965 if ( barr_kpi2 < 0.0 ) barr_kpi2 = barrier( 1, sKmPip2, sKm, sPip2, rRes2 );
966 if ( barr_kspi1_kpi2 < 0.0 ) barr_kspi1_kpi2 = barrier( 1, sD, sKsPip1, sKmPip2, rD2 );
967 for ( int a = 0; a < 4; a++ ) { temp_PDF += t1D[a] * t1KmPip2[a] * G[a][a]; }
968 tmp1 = barr_kpi2 * barr_kspi1_kpi2 * temp_PDF;
969
970 temp_PDF = 0;
971 if ( ispro_kpi1 == 0 )
972 {
973 propagatorRBW( mass1sq, mass1[i], width1[i], sKmPip1, sKm, sPip1, rRes2, 1, pro_kpi1 );
974 ispro_kpi1 = 1;
975 }
976 if ( ispro_kspi2sw == 0 )
977 {
978 KPiSLASS( sKsPip2, sKs, sPip2, pro_kspi2sw );
979 ispro_kspi2sw = 1;
980 }
981 Com_Multi( pro_kpi1, pro_kspi2sw, pro2b );
982
983 calt1( pKsPip2, pKmPip1, t1D );
984 if ( barr_kpi1 < 0.0 ) barr_kpi1 = barrier( 1, sKmPip1, sKm, sPip1, rRes2 );
985 if ( barr_kspi2_kpi1 < 0.0 ) barr_kspi2_kpi1 = barrier( 1, sD, sKsPip2, sKmPip1, rD2 );
986 for ( int a = 0; a < 4; a++ ) { temp_PDF += t1D[a] * t1KmPip1[a] * G[a][a]; }
987 tmp2 = barr_kspi2_kpi1 * barr_kpi1 * temp_PDF;
988 }
989
990 // Kspi_swave kpi_swave
991 if ( modetype[i] == 1 )
992 {
993 if ( ispro_kspi1sw == 0 )
994 {
995 KPiSLASS( sKsPip1, sKs, sPip1, pro_kspi1sw );
996 ispro_kspi1sw = 1;
997 }
998 if ( ispro_kpi2sw == 0 )
999 {
1000 KPiSLASS( sKmPip2, sKm, sPip2, pro_kpi2sw );
1001 ispro_kpi2sw = 1;
1002 }
1003 Com_Multi( pro_kspi1sw, pro_kpi2sw, pro1b );
1004 tmp1 = 1;
1005
1006 temp_PDF = 0;
1007 if ( ispro_kspi2sw == 0 )
1008 {
1009 KPiSLASS( sKsPip2, sKs, sPip2, pro_kspi2sw );
1010 ispro_kspi2sw = 1;
1011 }
1012 if ( ispro_kpi1sw == 0 )
1013 {
1014 KPiSLASS( sKmPip1, sKm, sPip1, pro_kpi1sw );
1015 ispro_kpi1sw = 1;
1016 }
1017 Com_Multi( pro_kspi2sw, pro_kpi1sw, pro2b );
1018 tmp2 = 1;
1019 }
1020
1021 // Ds -> eta(1475)pi+ eta(1475)->K*Ks,K*0>Km pi+
1022 if ( modetype[i] == 4 )
1023 {
1024 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKmPip1, sKs, rRes2, 1, pro1 );
1025 if ( ispro_kpi1 == 0 )
1026 {
1027 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip1, sKm, sPip1, rRes2, 1, pro_kpi1 );
1028 ispro_kpi1 = 1;
1029 }
1030 Com_Multi( pro1, pro_kpi1, pro1b );
1031
1032 if ( barr_kpi1 < 0.0 ) barr_kpi1 = barrier( 1, sKmPip1, sKm, sPip1, rRes2 );
1033 barr_kskpi1_kpi = barrier( 1, sKsKmPip1, sKmPip1, sKs, rRes2 );
1034 for ( int a = 0; a < 4; a++ ) { temp_PDF += Km[a] * t1KmPip1[a] * G[a][a]; }
1035 tmp1 = barr_kpi1 * barr_kskpi1_kpi * temp_PDF;
1036
1037 temp_PDF = 0;
1038 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKmPip2, sKs, rRes2, 1, pro1 );
1039 if ( ispro_kpi2 == 0 )
1040 {
1041 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip2, sKm, sPip2, rRes2, 1, pro_kpi2 );
1042 ispro_kpi2 = 1;
1043 }
1044 Com_Multi( pro1, pro_kpi2, pro2b );
1045
1046 if ( barr_kpi2 < 0.0 ) barr_kpi2 = barrier( 1, sKmPip2, sKm, sPip2, rRes2 );
1047 barr_kskpi2_kpi = barrier( 1, sKsKmPip2, sKmPip2, sKs, rRes2 );
1048 for ( int a = 0; a < 4; a++ ) { temp_PDF += Km[a] * t1KmPip2[a] * G[a][a]; }
1049 tmp2 = barr_kpi2 * barr_kskpi2_kpi * temp_PDF;
1050 }
1051
1052 // Ds -> eta(1475)pi+ eta(1475)->K*Km,K*+>Ks pi+
1053 if ( modetype[i] == 14 )
1054 {
1055 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKsPip1, sKm, rRes2, 1, pro1 );
1056 if ( ispro_kspi1 == 0 )
1057 {
1058 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip1, sKs, sPip1, rRes2, 1,
1059 pro_kspi1 );
1060 ispro_kspi1 = 1;
1061 }
1062 Com_Multi( pro1, pro_kspi1, pro1b );
1063
1064 if ( barr_kspi1 < 0.0 ) barr_kspi1 = barrier( 1, sKsPip1, sKs, sPip1, rRes2 );
1065 if ( barr_kskpi1_kspi < 0.0 )
1066 barr_kskpi1_kspi = barrier( 1, sKsKmPip1, sKsPip1, sKm, rRes2 );
1067 for ( int a = 0; a < 4; a++ ) { temp_PDF += Kp[a] * t1KsPip1[a] * G[a][a]; }
1068 tmp1 = barr_kspi1 * barr_kskpi1_kspi * temp_PDF;
1069
1070 temp_PDF = 0;
1071 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKsPip2, sKm, rRes2, 1, pro1 );
1072 if ( ispro_kspi2 == 0 )
1073 {
1074 propagatorRBW( mass1sq, mass2[i], width2[i], sKsPip2, sKs, sPip2, rRes2, 1,
1075 pro_kspi2 );
1076 ispro_kspi2 = 1;
1077 }
1078 Com_Multi( pro1, pro_kspi2, pro2b );
1079
1080 if ( barr_kspi2 < 0.0 ) barr_kspi2 = barrier( 1, sKsPip2, sKs, sPip2, rRes2 );
1081 if ( barr_kskpi2_kspi < 0.0 )
1082 barr_kskpi2_kspi = barrier( 1, sKsKmPip2, sKsPip2, sKm, rRes2 );
1083 for ( int a = 0; a < 4; a++ ) { temp_PDF += Kp[a] * t1KsPip2[a] * G[a][a]; }
1084 tmp2 = barr_kspi2 * barr_kskpi2_kspi * temp_PDF;
1085 }
1086 if ( modetype[i] == 8 )
1087 {
1088 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKsKm, sPip1, rRes2, 0, pro1 );
1089 if ( ispro_ksk == 0 )
1090 {
1091 propagatora0980( mass2sq, mass2[i], sKsKm, sKs2, sKm2, pro_ksk );
1092 ispro_ksk = 1;
1093 }
1094 Com_Multi( pro1, pro_ksk, pro1b );
1095 tmp1 = 1.0;
1096
1097 temp_PDF = 0;
1098 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKsKm, sPip2, rRes2, 0, pro1 );
1099 if ( ispro_ksk == 0 )
1100 {
1101 propagatora0980( mass2sq, mass2[i], sKsKm, sKs2, sKm2, pro_ksk );
1102 ispro_ksk = 1;
1103 }
1104 Com_Multi( pro1, pro_ksk, pro2b );
1105 tmp2 = 1.0;
1106 }
1107
1108 // f1285->a0 pip
1109 if ( modetype[i] == 17 )
1110 {
1111 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKsKm, sPip1, rRes2, 1, pro1 );
1112 if ( ispro_ksk == 0 )
1113 {
1114 propagatora0980( mass2sq, mass2[i], sKsKm, sKs2, sKm2, pro_ksk );
1115 ispro_ksk = 1;
1116 }
1117 Com_Multi( pro1, pro_ksk, pro1b );
1118 calt1( pKsKmPip1, Pi0, t1D );
1119
1120 if ( barr_kskpi1_ksk < 0.0 )
1121 barr_kskpi1_ksk = barrier( 1, sKsKmPip1, sKsKm, sPip1, rRes2 );
1122 if ( barr_ds_kskpi1 < 0.0 ) barr_ds_kskpi1 = barrier( 1, sD, sKsKmPip1, sPip2, rD2 );
1123 for ( int a = 0; a < 4; a++ ) { temp_PDF += G[a][a] * t1D[a] * t1KsKmPip1_ksk[a]; }
1124 tmp1 = barr_kskpi1_ksk * barr_ds_kskpi1 * temp_PDF;
1125
1126 temp_PDF = 0;
1127 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKsKm, sPip2, rRes2, 1, pro1 );
1128 if ( ispro_ksk == 0 )
1129 {
1130 propagatora0980( mass2sq, mass2[i], sKsKm, sKs2, sKm2, pro_ksk );
1131 ispro_ksk = 1;
1132 }
1133 Com_Multi( pro1, pro_ksk, pro2b );
1134 calt1( pKsKmPip2, Pip, t1D );
1135
1136 if ( barr_kskpi2_ksk < 0.0 )
1137 barr_kskpi2_ksk = barrier( 1, sKsKmPip2, sKsKm, sPip2, rRes2 );
1138 if ( barr_ds_kskpi2 < 0.0 ) barr_ds_kskpi2 = barrier( 1, sD, sKsKmPip2, sPip1, rD2 );
1139 for ( int a = 0; a < 4; a++ ) { temp_PDF += G[a][a] * t1D[a] * t1KsKmPip2_ksk[a]; }
1140 tmp2 = barr_kskpi2_ksk * barr_ds_kskpi2 * temp_PDF;
1141 }
1142
1143 // 1+ Ds -> 1550 pi+ , 1550 -> K-K*+ ,K*+ ->Ks pi+
1144 if ( modetype[i] == 5 )
1145 {
1146 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKsPip1, sKm, rRes2, g2[i],
1147 pro1 );
1148 // pro1[0]=1;pro1[1]=0;
1149 if ( ispro_kspi1 == 0 )
1150 {
1151 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip1, sKs, sPip1, rRes2, 1,
1152 pro_kspi1 );
1153 ispro_kspi1 = 1;
1154 }
1155 Com_Multi( pro1, pro_kspi1, pro1b );
1156 calt1( pKsKmPip1, Pi0, t1D );
1157
1158 if ( barr_kspi1 < 0.0 ) barr_kspi1 = barrier( 1, sKsPip1, sKs, sPip1, rRes2 );
1159 if ( barr_ds_kskpi1 < 0.0 ) barr_ds_kskpi1 = barrier( 1, sD, sKsKmPip1, sPip2, rD2 );
1160 barr_kskpi1_kspi = barrier( 2, sKsKmPip1, sKsPip1, sKm, rRes2 );
1161 if ( g2[i] == 0 )
1162 {
1163 for ( int a = 0; a < 4; a++ )
1164 {
1165 for ( int j = 0; j < 4; j++ )
1166 {
1167 temp_PDF += t1D[a] * ( pKsKmPip1[a] * pKsKmPip1[j] / sKsKmPip1 - G[a][j] ) *
1168 t1KsPip1[j] * G[a][a] * G[j][j];
1169 }
1170 }
1171 tmp1 = barr_kspi1 * barr_ds_kskpi1 * temp_PDF;
1172 }
1173 if ( g2[i] == 2 )
1174 {
1175 for ( int a = 0; a < 4; a++ )
1176 {
1177 for ( int j = 0; j < 4; j++ )
1178 { temp_PDF += t1D[a] * t2KsKmPip1_kspip[a][j] * t1KsPip1[j] * G[a][a] * G[j][j]; }
1179 }
1180 tmp1 = barr_kspi1 * barr_ds_kskpi1 * barr_kskpi1_kspi * temp_PDF;
1181 }
1182
1183 temp_PDF = 0;
1184 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKsPip2, sKm, rRes2, g2[i],
1185 pro2 );
1186 // pro2[0]=1;pro2[1]=0;
1187 if ( ispro_kspi2 == 0 )
1188 {
1189 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip2, sKs, sPip2, rRes2, 1,
1190 pro_kspi2 );
1191 ispro_kspi2 = 1;
1192 }
1193 Com_Multi( pro2, pro_kspi2, pro2b );
1194
1195 calt1( pKsKmPip2, Pip, t1D );
1196 if ( barr_kspi2 < 0.0 ) barr_kspi2 = barrier( 1, sKsPip2, sKs, sPip2, rRes2 );
1197 if ( barr_ds_kskpi2 < 0.0 ) barr_ds_kskpi2 = barrier( 1, sD, sKsKmPip2, sPip1, rD2 );
1198 barr_kskpi2_kspi = barrier( 2, sKsKmPip2, sKsPip2, sKm, rRes2 );
1199 if ( g2[i] == 0 )
1200 {
1201 for ( int a = 0; a < 4; a++ )
1202 {
1203 for ( int j = 0; j < 4; j++ )
1204 {
1205 temp_PDF += t1D[a] * ( pKsKmPip2[a] * pKsKmPip2[j] / sKsKmPip2 - G[a][j] ) *
1206 t1KsPip2[j] * G[a][a] * G[j][j];
1207 }
1208 }
1209 tmp2 = barr_kspi2 * barr_ds_kskpi2 * temp_PDF;
1210 }
1211 if ( g2[i] == 2 )
1212 {
1213 for ( int a = 0; a < 4; a++ )
1214 {
1215 for ( int j = 0; j < 4; j++ )
1216 { temp_PDF += t1D[a] * t2KsKmPip2_kspip[a][j] * t1KsPip2[j] * G[a][a] * G[j][j]; }
1217 }
1218 tmp2 = barr_kspi2 * barr_ds_kskpi2 * barr_kskpi2_kspi * temp_PDF;
1219 }
1220 }
1221
1222 // Ds -> 1550 pi+ , 1550 -> KsK*0bar ,K*0bar ->Km pi+
1223 if ( modetype[i] == 26 )
1224 {
1225 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKmPip1, sKs, rRes2, g2[i],
1226 pro1 );
1227 // pro1[0]=1;pro1[1]=0;
1228 if ( ispro_kpi1 == 0 )
1229 {
1230 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip1, sKm, sPip1, rRes2, 1, pro_kpi1 );
1231 ispro_kpi1 = 1;
1232 }
1233 Com_Multi( pro1, pro_kpi1, pro1b );
1234 calt1( pKsKmPip1, Pi0, t1D );
1235
1236 if ( barr_kpi1 < 0.0 ) barr_kpi1 = barrier( 1, sKmPip1, sKm, sPip1, rRes2 );
1237 if ( barr_ds_kskpi1 < 0.0 ) barr_ds_kskpi1 = barrier( 1, sD, sKsKmPip1, sPip2, rD2 );
1238 barr_kskpi1_kpi = barrier( 2, sKsKmPip1, sKmPip1, sKs, rRes2 );
1239 if ( g2[i] == 0 )
1240 {
1241 for ( int a = 0; a < 4; a++ )
1242 {
1243 for ( int j = 0; j < 4; j++ )
1244 {
1245 temp_PDF += t1D[a] * ( pKsKmPip1[a] * pKsKmPip1[j] / sKsKmPip1 - G[a][j] ) *
1246 t1KmPip1[j] * G[a][a] * G[j][j];
1247 }
1248 }
1249 tmp1 = barr_kpi1 * barr_ds_kskpi1 * temp_PDF;
1250 }
1251 if ( g2[i] == 2 )
1252 {
1253 for ( int a = 0; a < 4; a++ )
1254 {
1255 for ( int j = 0; j < 4; j++ )
1256 { temp_PDF += t1D[a] * t2KsKmPip1_kmpip[a][j] * t1KmPip1[j] * G[a][a] * G[j][j]; }
1257 }
1258 tmp1 = barr_kpi1 * barr_ds_kskpi1 * barr_kskpi1_kpi * temp_PDF;
1259 }
1260
1261 temp_PDF = 0;
1262 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKmPip2, sKs, rRes2, g2[i],
1263 pro2 );
1264 // pro2[0]=1;pro2[1]=0;
1265 if ( ispro_kpi2 == 0 )
1266 {
1267 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip2, sKm, sPip2, rRes2, 1, pro_kpi2 );
1268 ispro_kpi2 = 1;
1269 }
1270 Com_Multi( pro2, pro_kpi2, pro2b );
1271
1272 calt1( pKsKmPip2, Pip, t1D );
1273 if ( barr_kpi2 < 0.0 ) barr_kpi2 = barrier( 1, sKmPip2, sKm, sPip2, rRes2 );
1274 if ( barr_ds_kskpi2 < 0.0 ) barr_ds_kskpi2 = barrier( 1, sD, sKsKmPip2, sPip1, rD2 );
1275 barr_kskpi2_kpi = barrier( 2, sKsKmPip2, sKmPip2, sKs, rRes2 );
1276 if ( g2[i] == 0 )
1277 {
1278 for ( int a = 0; a < 4; a++ )
1279 {
1280 for ( int j = 0; j < 4; j++ )
1281 {
1282 temp_PDF += t1D[a] * ( pKsKmPip2[a] * pKsKmPip2[j] / sKsKmPip2 - G[a][j] ) *
1283 t1KmPip2[j] * G[a][a] * G[j][j];
1284 }
1285 }
1286 tmp2 = barr_kpi2 * barr_ds_kskpi2 * temp_PDF;
1287 }
1288 if ( g2[i] == 2 )
1289 {
1290 for ( int a = 0; a < 4; a++ )
1291 {
1292 for ( int j = 0; j < 4; j++ )
1293 { temp_PDF += t1D[a] * t2KsKmPip2_kmpip[a][j] * t1KmPip2[j] * G[a][a] * G[j][j]; }
1294 }
1295 tmp2 = barr_kpi2 * barr_ds_kskpi2 * barr_kskpi2_kpi * temp_PDF;
1296 }
1297 }
1298
1299 // text
1300 // Ds -> eta(1405)pi+ eta(1405)->K*Ks,K*0>Km pi+
1301 if ( modetype[i] == 7 )
1302 {
1303 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKmPip1, sKs, rRes2, 1, pro1 );
1304 if ( ispro_kpi1 == 0 )
1305 {
1306 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip1, sKm, sPip1, rRes2, 1, pro_kpi1 );
1307 ispro_kpi1 = 1;
1308 }
1309 Com_Multi( pro1, pro_kpi1, pro1b );
1310
1311 if ( barr_kpi1 < 0.0 ) barr_kpi1 = barrier( 1, sKmPip1, sKm, sPip1, rRes2 );
1312 barr_kskpi1_kpi = barrier( 1, sKsKmPip1, sKmPip1, sKs, rRes2 );
1313 for ( int a = 0; a < 4; a++ ) { temp_PDF += Km[a] * t1KmPip1[a] * G[a][a]; }
1314 tmp1 = barr_kpi1 * barr_kskpi1_kpi * temp_PDF;
1315
1316 temp_PDF = 0;
1317 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKmPip2, sKs, rRes2, 1, pro1 );
1318 if ( ispro_kpi2 == 0 )
1319 {
1320 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip2, sKm, sPip2, rRes2, 1, pro_kpi2 );
1321 ispro_kpi2 = 1;
1322 }
1323 Com_Multi( pro1, pro_kpi2, pro2b );
1324
1325 if ( barr_kpi2 < 0.0 ) barr_kpi2 = barrier( 1, sKmPip2, sKm, sPip2, rRes2 );
1326 barr_kskpi2_kpi = barrier( 1, sKsKmPip2, sKmPip2, sKs, rRes2 );
1327 for ( int a = 0; a < 4; a++ ) { temp_PDF += Km[a] * t1KmPip2[a] * G[a][a]; }
1328 tmp2 = barr_kpi2 * barr_kskpi2_kpi * temp_PDF;
1329 }
1330
1331 // text
1332 // Ds -> eta(1405)pi+ eta(1405)->K*Km,K*+>Ks pi+
1333 if ( modetype[i] == 10 )
1334 {
1335 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKsPip1, sKm, rRes2, 1, pro1 );
1336 if ( ispro_kspi1 == 0 )
1337 {
1338 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip1, sKs, sPip1, rRes2, 1,
1339 pro_kspi1 );
1340 ispro_kspi1 = 1;
1341 }
1342 Com_Multi( pro1, pro_kspi1, pro1b );
1343
1344 if ( barr_kspi1 < 0.0 ) barr_kspi1 = barrier( 1, sKsPip1, sKs, sPip1, rRes2 );
1345 if ( barr_kskpi1_kspi < 0.0 )
1346 barr_kskpi1_kspi = barrier( 1, sKsKmPip1, sKsPip1, sKm, rRes2 );
1347 for ( int a = 0; a < 4; a++ ) { temp_PDF += Kp[a] * t1KsPip1[a] * G[a][a]; }
1348 tmp1 = barr_kspi1 * barr_kskpi1_kspi * temp_PDF;
1349
1350 temp_PDF = 0;
1351 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKsPip2, sKm, rRes2, 1, pro1 );
1352 if ( ispro_kspi2 == 0 )
1353 {
1354 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip2, sKs, sPip2, rRes2, 1,
1355 pro_kspi2 );
1356 ispro_kspi2 = 1;
1357 }
1358 Com_Multi( pro1, pro_kspi2, pro2b );
1359
1360 if ( barr_kspi2 < 0.0 ) barr_kspi2 = barrier( 1, sKsPip2, sKs, sPip2, rRes2 );
1361 if ( barr_kskpi2_kspi < 0.0 )
1362 barr_kskpi2_kspi = barrier( 1, sKsKmPip2, sKsPip2, sKm, rRes2 );
1363 for ( int a = 0; a < 4; a++ ) { temp_PDF += Kp[a] * t1KsPip2[a] * G[a][a]; }
1364 tmp2 = barr_kspi2 * barr_kskpi2_kspi * temp_PDF;
1365 }
1366
1367 // text
1368 // eta1405->a0 pip
1369 if ( modetype[i] == 9 )
1370 {
1371 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKsKm, sPip1, rRes2, 0, pro1 );
1372 if ( ispro_ksk == 0 )
1373 {
1374 propagatora0980( mass2sq, mass2[i], sKsKm, sKs2, sKm2, pro_ksk );
1375 ispro_ksk = 1;
1376 }
1377 Com_Multi( pro1, pro_ksk, pro1b );
1378 tmp1 = 1.0;
1379
1380 temp_PDF = 0;
1381 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKsKm, sPip2, rRes2, 0, pro1 );
1382 if ( ispro_ksk == 0 )
1383 {
1384 propagatora0980( mass2sq, mass2[i], sKsKm, sKs2, sKm2, pro_ksk );
1385 ispro_ksk = 1;
1386 }
1387 Com_Multi( pro1, pro_ksk, pro2b );
1388 tmp2 = 1.0;
1389 }
1390
1391 // text
1392 // Ds -> eta(1405)pi+ eta(1405)->kpisw-ks
1393 if ( modetype[i] == 133 )
1394 {
1395 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKmPip1, sKs, rRes2, 0, pro1 );
1396 if ( ispro_kpi1sw == 0 )
1397 {
1398 KPiSLASS( sKmPip1, sKm, sPip1, pro_kpi1sw );
1399 ispro_kpi1sw = 1;
1400 }
1401 Com_Multi( pro1, pro_kpi1sw, pro1b );
1402 tmp1 = 1.0;
1403
1404 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKmPip2, sKs, rRes2, 0, pro1 );
1405 if ( ispro_kpi2sw == 0 )
1406 {
1407 KPiSLASS( sKmPip2, sKm, sPip2, pro_kpi2sw );
1408 ispro_kpi2sw = 1;
1409 }
1410 Com_Multi( pro1, pro_kpi2sw, pro2b );
1411 tmp2 = 1.0;
1412 }
1413
1414 // Ds -> eta(1405)pi+ eta(1405)->kpisw-ks
1415 if ( modetype[i] == 13 )
1416 {
1417 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKmPip1, sKs, rRes2, 0, pro1 );
1418 if ( ispro_kpi1sw == 0 )
1419 {
1420 KPiSLASS( sKmPip1, sKm, sPip1, pro_kpi1sw );
1421 ispro_kpi1sw = 1;
1422 }
1423 Com_Multi( pro1, pro_kpi1sw, pro1b );
1424 tmp1 = 1.0;
1425
1426 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKmPip2, sKs, rRes2, 0, pro1 );
1427 if ( ispro_kpi2sw == 0 )
1428 {
1429 KPiSLASS( sKmPip2, sKm, sPip2, pro_kpi2sw );
1430 ispro_kpi2sw = 1;
1431 }
1432 Com_Multi( pro1, pro_kpi2sw, pro2b );
1433 tmp2 = 1.0;
1434 }
1435
1436 // text
1437 // Ds -> eta(1405)pi+ eta(1405)->kspisw-k
1438 if ( modetype[i] == 20 )
1439 {
1440 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKsPip1, sKm, rRes2, 0, pro1 );
1441 if ( ispro_kspi1sw == 0 )
1442 {
1443 KPiSLASS( sKsPip1, sKs, sPip1, pro_kspi1sw );
1444 ispro_kspi1sw = 1;
1445 }
1446 Com_Multi( pro1, pro_kspi1sw, pro1b );
1447 tmp1 = 1.0;
1448
1449 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKsPip2, sKm, rRes2, 0, pro1 );
1450 if ( ispro_kspi2sw == 0 )
1451 {
1452 KPiSLASS( sKsPip2, sKs, sPip2, pro_kspi2sw );
1453 ispro_kspi2sw = 1;
1454 }
1455 Com_Multi( pro1, pro_kspi2sw, pro2b );
1456 tmp2 = 1.0;
1457 }
1458
1459 // text
1460 // Ds -> eta(1475)pi+ eta(1475)->kpisw-kspisw
1461 if ( modetype[i] == 19 )
1462 {
1463 if ( ispro_kspi1sw == 0 )
1464 {
1465 KPiSLASS( sKsPip1, sKs, sPip1, pro_kspi1sw );
1466 ispro_kspi1sw = 1;
1467 }
1468 if ( ispro_kpi1sw == 0 )
1469 {
1470 KPiSLASS( sKmPip1, sKm, sPip1, pro_kpi1sw );
1471 ispro_kpi1sw = 1;
1472 }
1473 Com_Multi( pro_kpi1sw, pro_kspi1sw, pro1b );
1474 tmp1 = 1.0;
1475
1476 if ( ispro_kspi2sw == 0 )
1477 {
1478 KPiSLASS( sKsPip2, sKs, sPip2, pro_kspi2sw );
1479 ispro_kspi2sw = 1;
1480 }
1481 if ( ispro_kpi2sw == 0 )
1482 {
1483 KPiSLASS( sKmPip2, sKm, sPip2, pro_kpi2sw );
1484 ispro_kpi2sw = 1;
1485 }
1486 Com_Multi( pro_kpi2sw, pro_kspi2sw, pro2b );
1487 tmp2 = 1.0;
1488 }
1489
1490 // text
1491 // Ds -> f1420 pi+ , f1420 -> K-K*+ ,K*+ ->Ks pi+
1492 if ( modetype[i] == 11 )
1493 {
1494 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKsPip1, sKm, rRes2, g2[i],
1495 pro1 );
1496 if ( ispro_kspi1 == 0 )
1497 {
1498 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip1, sKs, sPip1, rRes2, 1,
1499 pro_kspi1 );
1500 ispro_kspi1 = 1;
1501 }
1502 Com_Multi( pro1, pro_kspi1, pro1b );
1503 calt1( pKsKmPip1, Pi0, t1D );
1504
1505 if ( barr_kspi1 < 0.0 ) barr_kspi1 = barrier( 1, sKsPip1, sKs, sPip1, rRes2 );
1506 if ( barr_ds_kskpi1 < 0.0 ) barr_ds_kskpi1 = barrier( 1, sD, sKsKmPip1, sPip2, rD2 );
1507 barr_kskpi1_kspi = barrier( 2, sKsKmPip1, sKsPip1, sKm, rRes2 );
1508 if ( g2[i] == 0 )
1509 {
1510 for ( int a = 0; a < 4; a++ )
1511 {
1512 for ( int j = 0; j < 4; j++ )
1513 {
1514 temp_PDF += t1D[a] * ( pKsKmPip1[a] * pKsKmPip1[j] / sKsKmPip1 - G[a][j] ) *
1515 t1KsPip1[j] * G[a][a] * G[j][j];
1516 }
1517 }
1518 tmp1 = barr_kspi1 * barr_ds_kskpi1 * temp_PDF;
1519 }
1520 if ( g2[i] == 2 )
1521 {
1522 for ( int a = 0; a < 4; a++ )
1523 {
1524 for ( int j = 0; j < 4; j++ )
1525 { temp_PDF += t1D[a] * t2KsKmPip1_kspip[a][j] * t1KsPip1[j] * G[a][a] * G[j][j]; }
1526 }
1527 tmp1 = barr_kspi1 * barr_ds_kskpi1 * barr_kskpi1_kspi * temp_PDF;
1528 }
1529
1530 temp_PDF = 0;
1531 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKsPip2, sKm, rRes2, g2[i],
1532 pro1 );
1533 if ( ispro_kspi2 == 0 )
1534 {
1535 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip2, sKs, sPip2, rRes2, 1,
1536 pro_kspi2 );
1537 ispro_kspi2 = 1;
1538 }
1539 Com_Multi( pro1, pro_kspi2, pro2b );
1540
1541 calt1( pKsKmPip2, Pip, t1D );
1542 if ( barr_kspi2 < 0.0 ) barr_kspi2 = barrier( 1, sKsPip2, sKs, sPip2, rRes2 );
1543 if ( barr_ds_kskpi2 < 0.0 ) barr_ds_kskpi2 = barrier( 1, sD, sKsKmPip2, sPip1, rD2 );
1544 barr_kskpi2_kspi = barrier( 2, sKsKmPip2, sKsPip2, sKm, rRes2 );
1545 if ( g2[i] == 0 )
1546 {
1547 for ( int a = 0; a < 4; a++ )
1548 {
1549 for ( int j = 0; j < 4; j++ )
1550 {
1551 temp_PDF += t1D[a] * ( pKsKmPip2[a] * pKsKmPip2[j] / sKsKmPip2 - G[a][j] ) *
1552 t1KsPip2[j] * G[a][a] * G[j][j];
1553 }
1554 }
1555 tmp2 = barr_kspi2 * barr_ds_kskpi2 * temp_PDF;
1556 }
1557 if ( g2[i] == 2 )
1558 {
1559 for ( int a = 0; a < 4; a++ )
1560 {
1561 for ( int j = 0; j < 4; j++ )
1562 { temp_PDF += t1D[a] * t2KsKmPip2_kspip[a][j] * t1KsPip2[j] * G[a][a] * G[j][j]; }
1563 }
1564 tmp2 = barr_kspi2 * barr_ds_kskpi2 * barr_kskpi2_kspi * temp_PDF;
1565 }
1566 }
1567 if ( modetype[i] == 1200 )
1568 {
1569 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKmPip1, sKs, rRes2, g2[i],
1570 pro1 );
1571
1572 pro1[0] = 1;
1573 pro1[1] = 0;
1574 if ( ispro_kpi1 == 0 )
1575 {
1576 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip1, sKm, sPip1, rRes2, 1, pro_kpi1 );
1577 ispro_kpi1 = 1;
1578 }
1579 Com_Multi( pro1, pro_kpi1, pro1b );
1580 calt1( pKsKmPip1, Pi0, t1D );
1581
1582 if ( barr_kpi1 < 0.0 ) barr_kpi1 = barrier( 1, sKmPip1, sKm, sPip1, rRes2 );
1583 if ( barr_ds_kskpi1 < 0.0 ) barr_ds_kskpi1 = barrier( 1, sD, sKsKmPip1, sPip2, rD2 );
1584 barr_kskpi1_kpi = barrier( 2, sKsKmPip1, sKmPip1, sKs, rRes2 );
1585 if ( g2[i] == 0 )
1586 {
1587 for ( int a = 0; a < 4; a++ )
1588 {
1589 for ( int j = 0; j < 4; j++ )
1590 {
1591 temp_PDF += t1D[a] * ( pKsKmPip1[a] * pKsKmPip1[j] / sKsKmPip1 - G[a][j] ) *
1592 t1KmPip1[j] * G[a][a] * G[j][j];
1593 }
1594 }
1595 tmp1 = barr_kpi1 * barr_ds_kskpi1 * temp_PDF;
1596 }
1597 if ( g2[i] == 2 )
1598 {
1599 for ( int a = 0; i < 4; i++ )
1600 {
1601 for ( int j = 0; j < 4; j++ )
1602 { temp_PDF += t1D[a] * t2KsKmPip1_kmpip[a][j] * t1KmPip1[j] * G[a][a] * G[j][j]; }
1603 }
1604 tmp1 = barr_kpi1 * barr_ds_kskpi1 * barr_kskpi1_kpi * temp_PDF;
1605 }
1606
1607 temp_PDF = 0;
1608 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKmPip2, sKs, rRes2, g2[i],
1609 pro2 );
1610 pro2[0] = 1;
1611 pro2[1] = 0;
1612 if ( ispro_kpi2 == 0 )
1613 {
1614 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip2, sKm, sPip2, rRes2, 1, pro_kpi2 );
1615 ispro_kpi2 = 1;
1616 }
1617 Com_Multi( pro2, pro_kpi2, pro2b );
1618
1619 calt1( pKsKmPip2, Pip, t1D );
1620 if ( barr_kpi2 < 0.0 ) barr_kpi2 = barrier( 1, sKmPip2, sKm, sPip2, rRes2 );
1621 if ( barr_ds_kskpi2 < 0.0 ) barr_ds_kskpi2 = barrier( 1, sD, sKsKmPip2, sPip1, rD2 );
1622 barr_kskpi2_kpi = barrier( 2, sKsKmPip2, sKmPip2, sKs, rRes2 );
1623 if ( g2[i] == 0 )
1624 {
1625 for ( int a = 0; a < 4; a++ )
1626 {
1627 for ( int j = 0; j < 4; j++ )
1628 {
1629 temp_PDF += t1D[a] * ( pKsKmPip2[a] * pKsKmPip2[j] / sKsKmPip2 - G[a][j] ) *
1630 t1KmPip2[j] * G[a][a] * G[j][j];
1631 }
1632 }
1633 tmp2 = barr_kpi2 * barr_ds_kskpi2 * temp_PDF;
1634 }
1635 if ( g2[i] == 2 )
1636 {
1637 for ( int a = 0; a < 4; a++ )
1638 {
1639 for ( int j = 0; j < 4; j++ )
1640 { temp_PDF += t1D[a] * t2KsKmPip2_kmpip[a][j] * t1KmPip2[j] * G[a][a] * G[j][j]; }
1641 }
1642 tmp2 = barr_kpi2 * barr_ds_kskpi2 * barr_kskpi2_kpi * temp_PDF;
1643 }
1644 }
1645 // text
1646 // Ds -> f1420(nonResonance) pi+ , f1420 -> K-K*+ ,K*+ ->Ks pi+
1647 if ( modetype[i] == 1100 )
1648 {
1649 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKsPip1, sKm, rRes2, g2[i],
1650 pro1 );
1651 pro1[0] = 1;
1652 pro1[1] = 0;
1653 if ( ispro_kspi1 == 0 )
1654 {
1655 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip1, sKs, sPip1, rRes2, 1,
1656 pro_kspi1 );
1657 ispro_kspi1 = 1;
1658 }
1659 Com_Multi( pro1, pro_kspi1, pro1b );
1660 calt1( pKsKmPip1, Pi0, t1D );
1661
1662 if ( barr_kspi1 < 0.0 ) barr_kspi1 = barrier( 1, sKsPip1, sKs, sPip1, rRes2 );
1663 if ( barr_ds_kskpi1 < 0.0 ) barr_ds_kskpi1 = barrier( 1, sD, sKsKmPip1, sPip2, rD2 );
1664 barr_kskpi1_kspi = barrier( 2, sKsKmPip1, sKsPip1, sKm, rRes2 );
1665 if ( g2[i] == 0 )
1666 {
1667 for ( int a = 0; a < 4; a++ )
1668 {
1669 for ( int j = 0; j < 4; j++ )
1670 {
1671 temp_PDF += t1D[a] * ( pKsKmPip1[a] * pKsKmPip1[j] / sKsKmPip1 - G[a][j] ) *
1672 t1KsPip1[j] * G[a][a] * G[j][j];
1673 }
1674 }
1675 tmp1 = barr_kspi1 * barr_ds_kskpi1 * temp_PDF;
1676 }
1677 if ( g2[i] == 2 )
1678 {
1679 for ( int a = 0; a < 4; a++ )
1680 {
1681 for ( int j = 0; j < 4; j++ )
1682 { temp_PDF += t1D[a] * t2KsKmPip1_kspip[a][j] * t1KsPip1[j] * G[a][a] * G[j][j]; }
1683 }
1684 tmp1 = barr_kspi1 * barr_ds_kskpi1 * barr_kskpi1_kspi * temp_PDF;
1685 }
1686
1687 temp_PDF = 0;
1688 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKsPip2, sKm, rRes2, g2[i],
1689 pro1 );
1690 pro1[0] = 1;
1691 pro1[1] = 0;
1692 if ( ispro_kspi2 == 0 )
1693 {
1694 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip2, sKs, sPip2, rRes2, 1,
1695 pro_kspi2 );
1696 ispro_kspi2 = 1;
1697 }
1698 Com_Multi( pro1, pro_kspi2, pro2b );
1699
1700 calt1( pKsKmPip2, Pip, t1D );
1701 if ( barr_kspi2 < 0.0 ) barr_kspi2 = barrier( 1, sKsPip2, sKs, sPip2, rRes2 );
1702 if ( barr_ds_kskpi2 < 0.0 ) barr_ds_kskpi2 = barrier( 1, sD, sKsKmPip2, sPip1, rD2 );
1703 barr_kskpi2_kspi = barrier( 2, sKsKmPip2, sKsPip2, sKm, rRes2 );
1704 if ( g2[i] == 0 )
1705 {
1706 for ( int a = 0; a < 4; a++ )
1707 {
1708 for ( int j = 0; j < 4; j++ )
1709 {
1710 temp_PDF += t1D[a] * ( pKsKmPip2[a] * pKsKmPip2[j] / sKsKmPip2 - G[a][j] ) *
1711 t1KsPip2[j] * G[a][a] * G[j][j];
1712 }
1713 }
1714 tmp2 = barr_kspi2 * barr_ds_kskpi2 * temp_PDF;
1715 }
1716 if ( g2[i] == 2 )
1717 {
1718 for ( int a = 0; a < 4; a++ )
1719 {
1720 for ( int j = 0; j < 4; j++ )
1721 { temp_PDF += t1D[a] * t2KsKmPip2_kspip[a][j] * t1KsPip2[j] * G[a][a] * G[j][j]; }
1722 }
1723 tmp2 = barr_kspi2 * barr_ds_kskpi2 * barr_kskpi2_kspi * temp_PDF;
1724 }
1725 }
1726
1727 // text
1728 // Ds -> 1420 pi+ , 1420 -> KsK*0bar ,K*0bar ->Km pi+
1729 if ( modetype[i] == 12 )
1730 {
1731 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKmPip1, sKs, rRes2, g2[i],
1732 pro1 );
1733 if ( ispro_kpi1 == 0 )
1734 {
1735 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip1, sKm, sPip1, rRes2, 1, pro_kpi1 );
1736 ispro_kpi1 = 1;
1737 }
1738 Com_Multi( pro1, pro_kpi1, pro1b );
1739 calt1( pKsKmPip1, Pi0, t1D );
1740
1741 if ( barr_kpi1 < 0.0 ) barr_kpi1 = barrier( 1, sKmPip1, sKm, sPip1, rRes2 );
1742 if ( barr_ds_kskpi1 < 0.0 ) barr_ds_kskpi1 = barrier( 1, sD, sKsKmPip1, sPip2, rD2 );
1743 barr_kskpi1_kpi = barrier( 2, sKsKmPip1, sKmPip1, sKs, rRes2 );
1744 if ( g2[i] == 0 )
1745 {
1746 for ( int a = 0; a < 4; a++ )
1747 {
1748 for ( int j = 0; j < 4; j++ )
1749 {
1750 temp_PDF += t1D[a] * ( pKsKmPip1[a] * pKsKmPip1[j] / sKsKmPip1 - G[a][j] ) *
1751 t1KmPip1[j] * G[a][a] * G[j][j];
1752 }
1753 }
1754 tmp1 = barr_kpi1 * barr_ds_kskpi1 * temp_PDF;
1755 }
1756 if ( g2[i] == 2 )
1757 {
1758 for ( int a = 0; i < 4; i++ )
1759 {
1760 for ( int j = 0; j < 4; j++ )
1761 { temp_PDF += t1D[a] * t2KsKmPip1_kmpip[a][j] * t1KmPip1[j] * G[a][a] * G[j][j]; }
1762 }
1763 tmp1 = barr_kpi1 * barr_ds_kskpi1 * barr_kskpi1_kpi * temp_PDF;
1764 }
1765 temp_PDF = 0;
1766 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKmPip2, sKs, rRes2, g2[i],
1767 pro2 );
1768 if ( ispro_kpi2 == 0 )
1769 {
1770 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip2, sKm, sPip2, rRes2, 1, pro_kpi2 );
1771 ispro_kpi2 = 1;
1772 }
1773 Com_Multi( pro2, pro_kpi2, pro2b );
1774
1775 calt1( pKsKmPip2, Pip, t1D );
1776 if ( barr_kpi2 < 0.0 ) barr_kpi2 = barrier( 1, sKmPip2, sKm, sPip2, rRes2 );
1777 if ( barr_ds_kskpi2 < 0.0 ) barr_ds_kskpi2 = barrier( 1, sD, sKsKmPip2, sPip1, rD2 );
1778 barr_kskpi2_kpi = barrier( 2, sKsKmPip2, sKmPip2, sKs, rRes2 );
1779 if ( g2[i] == 0 )
1780 {
1781 for ( int a = 0; a < 4; a++ )
1782 {
1783 for ( int j = 0; j < 4; j++ )
1784 {
1785 temp_PDF += t1D[a] * ( pKsKmPip2[a] * pKsKmPip2[j] / sKsKmPip2 - G[a][j] ) *
1786 t1KmPip2[j] * G[a][a] * G[j][j];
1787 }
1788 }
1789 tmp2 = barr_kpi2 * barr_ds_kskpi2 * temp_PDF;
1790 }
1791 if ( g2[i] == 2 )
1792 {
1793 for ( int a = 0; a < 4; a++ )
1794 {
1795 for ( int j = 0; j < 4; j++ )
1796 { temp_PDF += t1D[a] * t2KsKmPip2_kmpip[a][j] * t1KmPip2[j] * G[a][a] * G[j][j]; }
1797 }
1798 tmp2 = barr_kpi2 * barr_ds_kskpi2 * barr_kskpi2_kpi * temp_PDF;
1799 }
1800 }
1801 if ( modetype[i] == 1500 )
1802 {
1803 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKsKm, sPip1, rRes2, 1, pro1 );
1804 pro1[0] = 1;
1805 pro1[1] = 0;
1806 if ( ispro_ksk == 0 )
1807 {
1808 propagatora0980( mass2sq, mass2[i], sKsKm, sKs2, sKm2, pro_ksk );
1809 ispro_ksk = 1;
1810 }
1811 Com_Multi( pro1, pro_ksk, pro1b );
1812 calt1( pKsKmPip1, Pi0, t1D );
1813
1814 if ( barr_kskpi1_ksk < 0.0 )
1815 barr_kskpi1_ksk = barrier( 1, sKsKmPip1, sKsKm, sPip1, rRes2 );
1816 if ( barr_ds_kskpi1 < 0.0 ) barr_ds_kskpi1 = barrier( 1, sD, sKsKmPip1, sPip2, rD2 );
1817 for ( int a = 0; a < 4; a++ ) { temp_PDF += G[a][a] * t1D[a] * t1KsKmPip1_ksk[a]; }
1818 tmp1 = barr_kskpi1_ksk * barr_ds_kskpi1 * temp_PDF;
1819
1820 temp_PDF = 0;
1821 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKsKm, sPip2, rRes2, 1, pro1 );
1822 pro1[0] = 1;
1823 pro1[1] = 0;
1824 if ( ispro_ksk == 0 )
1825 {
1826 propagatora0980( mass2sq, mass2[i], sKsKm, sKs2, sKm2, pro_ksk );
1827 ispro_ksk = 1;
1828 }
1829 Com_Multi( pro1, pro_ksk, pro2b );
1830 calt1( pKsKmPip2, Pip, t1D );
1831
1832 if ( barr_kskpi2_ksk < 0.0 )
1833 barr_kskpi2_ksk = barrier( 1, sKsKmPip2, sKsKm, sPip2, rRes2 );
1834 if ( barr_ds_kskpi2 < 0.0 ) barr_ds_kskpi2 = barrier( 1, sD, sKsKmPip2, sPip1, rD2 );
1835 for ( int a = 0; a < 4; a++ ) { temp_PDF += G[a][a] * t1D[a] * t1KsKmPip2_ksk[a]; }
1836 tmp2 = barr_kskpi2_ksk * barr_ds_kskpi2 * temp_PDF;
1837 }
1838
1839 // text
1840 /// f1420->a0 pip
1841 if ( modetype[i] == 15 )
1842 {
1843 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKsKm, sPip1, rRes2, 1, pro1 );
1844 if ( ispro_ksk == 0 )
1845 {
1846 propagatora0980( mass2sq, mass2[i], sKsKm, sKs2, sKm2, pro_ksk );
1847 ispro_ksk = 1;
1848 }
1849 Com_Multi( pro1, pro_ksk, pro1b );
1850 calt1( pKsKmPip1, Pi0, t1D );
1851
1852 if ( barr_kskpi1_ksk < 0.0 )
1853 barr_kskpi1_ksk = barrier( 1, sKsKmPip1, sKsKm, sPip1, rRes2 );
1854 if ( barr_ds_kskpi1 < 0.0 ) barr_ds_kskpi1 = barrier( 1, sD, sKsKmPip1, sPip2, rD2 );
1855 for ( int a = 0; a < 4; a++ ) { temp_PDF += G[a][a] * t1D[a] * t1KsKmPip1_ksk[a]; }
1856 tmp1 = barr_kskpi1_ksk * barr_ds_kskpi1 * temp_PDF;
1857
1858 temp_PDF = 0;
1859 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKsKm, sPip2, rRes2, 1, pro1 );
1860 if ( ispro_ksk == 0 )
1861 {
1862 propagatora0980( mass2sq, mass2[i], sKsKm, sKs2, sKm2, pro_ksk );
1863 ispro_ksk = 1;
1864 }
1865 Com_Multi( pro1, pro_ksk, pro2b );
1866 calt1( pKsKmPip2, Pip, t1D );
1867
1868 if ( barr_kskpi2_ksk < 0.0 )
1869 barr_kskpi2_ksk = barrier( 1, sKsKmPip2, sKsKm, sPip2, rRes2 );
1870 if ( barr_ds_kskpi2 < 0.0 ) barr_ds_kskpi2 = barrier( 1, sD, sKsKmPip2, sPip1, rD2 );
1871 for ( int a = 0; a < 4; a++ ) { temp_PDF += G[a][a] * t1D[a] * t1KsKmPip2_ksk[a]; }
1872 tmp2 = barr_kskpi2_ksk * barr_ds_kskpi2 * temp_PDF;
1873 }
1874
1875 // text
1876 /// f1420->kspi_sw
1877 if ( modetype[i] == 24 )
1878 {
1879 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKsPip1, sKm, rRes2, 1, pro1 );
1880 if ( ispro_kspi1sw == 0 )
1881 {
1882 KPiSLASS( sKsPip1, sKs, sPip1, pro_kspi1sw );
1883 ispro_kspi1sw = 1;
1884 }
1885 Com_Multi( pro1, pro_kspi1sw, pro1b );
1886 calt1( pKsKmPip1, Pi0, t1D );
1887
1888 if ( barr_kskpi1_kspi < 0.0 )
1889 barr_kskpi1_kspi = barrier( 1, sKsKmPip1, sKsPip1, sKm, rRes2 );
1890 if ( barr_ds_kskpi1 < 0.0 ) barr_ds_kskpi1 = barrier( 1, sD, sKsKmPip1, sPip2, rD2 );
1891 for ( int a = 0; a < 4; a++ ) { temp_PDF += G[a][a] * t1D[a] * t1KsKmPip1_kspip[a]; }
1892 tmp1 = barr_kskpi1_kspi * barr_ds_kskpi1 * temp_PDF;
1893
1894 temp_PDF = 0;
1895 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKsPip2, sKm, rRes2, 1, pro1 );
1896 if ( ispro_kspi2sw == 0 )
1897 {
1898 KPiSLASS( sKsPip2, sKs, sPip2, pro_kspi2sw );
1899 ispro_kspi2sw = 1;
1900 }
1901 Com_Multi( pro1, pro_kspi2sw, pro2b );
1902 calt1( pKsKmPip2, Pip, t1D );
1903
1904 if ( barr_kskpi2_kspi < 0.0 )
1905 barr_kskpi2_kspi = barrier( 1, sKsKmPip2, sKsPip2, sKm, rRes2 );
1906 if ( barr_ds_kskpi2 < 0.0 ) barr_ds_kskpi2 = barrier( 1, sD, sKsKmPip2, sPip1, rD2 );
1907 for ( int a = 0; a < 4; a++ ) { temp_PDF += G[a][a] * t1D[a] * t1KsKmPip2_kspip[a]; }
1908 tmp2 = barr_kskpi2_kspi * barr_ds_kskpi2 * temp_PDF;
1909 }
1910
1911 // text
1912 /// f1420->kpi_sw
1913 if ( modetype[i] == 27 )
1914 {
1915 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKmPip1, sKs, rRes2, 1, pro1 );
1916 if ( ispro_kspi1sw == 0 )
1917 {
1918 KPiSLASS( sKsPip1, sKm, sPip1, pro_kpi1sw );
1919 ispro_kpi1sw = 1;
1920 }
1921 Com_Multi( pro1, pro_kpi1sw, pro1b );
1922 calt1( pKsKmPip1, Pi0, t1D );
1923
1924 barr_kskpi1_kpi = barrier( 1, sKsKmPip1, sKmPip1, sKs, rRes2 );
1925 if ( barr_ds_kskpi1 < 0.0 ) barr_ds_kskpi1 = barrier( 1, sD, sKsKmPip1, sPip2, rD2 );
1926 for ( int a = 0; a < 4; a++ ) { temp_PDF += G[a][a] * t1D[a] * t1KsKmPip1_kmpip[a]; }
1927 tmp1 = barr_kskpi1_kpi * barr_ds_kskpi1 * temp_PDF;
1928
1929 temp_PDF = 0;
1930 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKmPip2, sKs, rRes2, 1, pro1 );
1931 if ( ispro_kpi2sw == 0 )
1932 {
1933 KPiSLASS( sKsPip2, sKs, sPip2, pro_kpi2sw );
1934 ispro_kpi2sw = 1;
1935 }
1936 Com_Multi( pro1, pro_kpi2sw, pro2b );
1937 calt1( pKsKmPip2, Pip, t1D );
1938
1939 barr_kskpi2_kpi = barrier( 1, sKsKmPip2, sKmPip2, sKs, rRes2 );
1940 if ( barr_ds_kskpi2 < 0.0 ) barr_ds_kskpi2 = barrier( 1, sD, sKsKmPip2, sPip1, rD2 );
1941 for ( int a = 0; a < 4; a++ ) { temp_PDF += G[a][a] * t1D[a] * t1KsKmPip2_kmpip[a]; }
1942 tmp2 = barr_kskpi2_kpi * barr_ds_kskpi2 * temp_PDF;
1943 }
1944
1945 // Ds -> eta(1475)pi+ eta(1475)->kappa Ks ,K*->Km pi+
1946 if ( modetype[i] == 244 )
1947 {
1948 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKsPip1, sKm, rRes2, 0, pro1 );
1949 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip1, sKs, sPip1, rRes2, 0, pro3 );
1950 Com_Multi( pro1, pro3, pro1b );
1951 tmp1 = 1.0;
1952
1953 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKsPip2, sKm, rRes2, 0, pro1 );
1954 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip2, sKs, sPip2, rRes2, 0, pro3 );
1955 Com_Multi( pro1, pro3, pro2b );
1956 tmp2 = 1.0;
1957 }
1958
1959 // Ds -> kst(1410)K- , phi(1680)->K-K*+ ,K*+ ->Ks pi+
1960 if ( modetype[i] == 280 )
1961 {
1962 propagatorRBW( mass1sq, mass1[i], width1[i], sKsPipPip, sKsPip1, sPip2, rRes2, 1, pro1 );
1963 if ( ispro_kspi1 == 0 )
1964 {
1965 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip1, sKs, sPip1, rRes2, 1,
1966 pro_kspi1 );
1967 ispro_kspi1 = 1;
1968 }
1969 Com_Multi( pro1, pro_kspi1, pro1b );
1970 // printf(" %.15f %.15f %.15f\n",B[0],B[1],B[2]);
1971
1972 B[0] = barrier( 1, sKsPip1, sKs, sPip1, rRes2 );
1973 B[1] = barrier( 1, sKsKmPip1, sKsPip1, sKm, rRes2 );
1974 B[2] = barrier( 1, sD, sKsKmPip1, sPip2, rD2 );
1975 for ( int a = 0; a < 4; a++ )
1976 {
1977 for ( int j = 0; j < 4; j++ )
1978 {
1979 for ( int k = 0; k < 4; k++ )
1980 {
1981 for ( int l = 0; l < 4; l++ )
1982 {
1983 temp_PDF += E[a][j][k][l] * pKsKmPip1[a] * ( pKsPip1[j] - Kp[j] ) * Pi0[k] *
1984 ( Km[l] - Pip[l] ) * G[a][a] * G[j][j] * G[k][k] * G[l][l];
1985 }
1986 }
1987 }
1988 }
1989 tmp1 = B[0] * B[1] * B[2] * temp_PDF;
1990
1991 temp_PDF = 0;
1992 propagatorRBW( mass1sq, mass1[i], width1[i], sKsPipPip, sKsPip2, sPip1, rRes2, 1, pro2 );
1993 if ( ispro_kspi2 == 0 )
1994 {
1995 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip2, sKs, sPip2, rRes2, 1,
1996 pro_kspi2 );
1997 ispro_kspi2 = 1;
1998 }
1999 Com_Multi( pro2, pro_kspi2, pro2b );
2000
2001 B[0] = barrier( 1, sKsPip2, sKs, sPip2, rRes2 );
2002 B[1] = barrier( 1, sKsKmPip2, sKsPip2, sKm, rRes2 );
2003 B[2] = barrier( 1, sD, sKsKmPip2, sPip1, rD2 );
2004
2005 for ( int a = 0; a < 4; a++ )
2006 {
2007 for ( int j = 0; j < 4; j++ )
2008 {
2009 for ( int k = 0; k < 4; k++ )
2010 {
2011 for ( int l = 0; l < 4; l++ )
2012 {
2013 temp_PDF += E[a][j][k][l] * pKsKmPip2[a] * ( pKsPip2[j] - Kp[j] ) * Pip[k] *
2014 ( Km[l] - Pi0[l] ) * G[a][a] * G[j][j] * G[k][k] * G[l][l];
2015 }
2016 }
2017 }
2018 }
2019 tmp2 = B[0] * B[1] * B[2] * temp_PDF;
2020 }
2021
2022 // text
2023 // Ds -> kst(1410)KS0 , phi(1680)->K-K*+ ,K*+ ->Km pi+
2024 if ( modetype[i] == 290 )
2025 {
2026 propagatorRBW( mass1sq, mass1[i], width1[i], sKmPipPip, sKmPip1, sPip2, rRes2, 1, pro1 );
2027 if ( ispro_kpi1 == 0 )
2028 {
2029 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip1, sKm, sPip1, rRes2, 1, pro_kpi1 );
2030 ispro_kpi1 = 1;
2031 }
2032 Com_Multi( pro1, pro_kpi1, pro1b );
2033
2034 B[0] = barrier( 1, sKmPip1, sKm, sPip1, rRes2 );
2035 B[1] = barrier( 1, sKsKmPip1, sKmPip1, sKs, rRes2 );
2036 B[2] = barrier( 1, sD, sKsKmPip1, sPip2, rD2 );
2037 for ( int a = 0; a < 4; a++ )
2038 {
2039 for ( int j = 0; j < 4; j++ )
2040 {
2041 for ( int k = 0; k < 4; k++ )
2042 {
2043 for ( int l = 0; l < 4; l++ )
2044 {
2045 temp_PDF += E[a][j][k][l] * pKsKmPip1[a] * ( pKmPip1[j] - Km[j] ) * Pi0[k] *
2046 ( Kp[l] - Pip[l] ) * G[a][a] * G[j][j] * G[k][k] * G[l][l];
2047 }
2048 }
2049 }
2050 }
2051 tmp1 = B[0] * B[1] * B[2] * temp_PDF;
2052
2053 temp_PDF = 0;
2054 propagatorRBW( mass1sq, mass1[i], width1[i], sKmPipPip, sKmPip2, sPip1, rRes2, 1, pro2 );
2055 if ( ispro_kpi2 == 0 )
2056 {
2057 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip2, sKm, sPip2, rRes2, 1, pro_kpi2 );
2058 ispro_kpi2 = 1;
2059 }
2060 Com_Multi( pro2, pro_kpi2, pro2b );
2061 B[0] = barrier( 1, sKmPip2, sKm, sPip2, rRes2 );
2062 B[1] = barrier( 1, sKsKmPip2, sKmPip2, sKs, rRes2 );
2063 B[2] = barrier( 1, sD, sKsKmPip2, sPip1, rD2 );
2064
2065 for ( int a = 0; a < 4; a++ )
2066 {
2067 for ( int j = 0; j < 4; j++ )
2068 {
2069 for ( int k = 0; k < 4; k++ )
2070 {
2071 for ( int l = 0; l < 4; l++ )
2072 {
2073 temp_PDF += E[a][j][k][l] * pKsKmPip2[a] * ( pKmPip2[j] - Km[j] ) * Pip[k] *
2074 ( Kp[l] - Pi0[l] ) * G[a][a] * G[j][j] * G[k][k] * G[l][l];
2075 }
2076 }
2077 }
2078 }
2079
2080 tmp2 = B[0] * B[1] * B[2] * temp_PDF * ( -1 );
2081 }
2082 // text
2083 // Ds -> phi(1680)pi+ , phi(1680)->K-K*+ ,K*+ ->Ks pi+
2084 if ( modetype[i] == 28 )
2085 {
2086 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKsPip1, sKm, rRes2, 1, pro1 );
2087 if ( ispro_kspi1 == 0 )
2088 {
2089 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip1, sKs, sPip1, rRes2, 1,
2090 pro_kspi1 );
2091 ispro_kspi1 = 1;
2092 }
2093 Com_Multi( pro1, pro_kspi1, pro1b );
2094 // printf(" %.15f %.15f %.15f\n",B[0],B[1],B[2]);
2095
2096 B[0] = barrier( 1, sKsPip1, sKs, sPip1, rRes2 );
2097 B[1] = barrier( 1, sKsKmPip1, sKsPip1, sKm, rRes2 );
2098 B[2] = barrier( 1, sD, sKsKmPip1, sPip2, rD2 );
2099 for ( int a = 0; a < 4; a++ )
2100 {
2101 for ( int j = 0; j < 4; j++ )
2102 {
2103 for ( int k = 0; k < 4; k++ )
2104 {
2105 for ( int l = 0; l < 4; l++ )
2106 {
2107 temp_PDF += E[a][j][k][l] * pKsKmPip1[a] * ( pKsPip1[j] - Kp[j] ) * Pi0[k] *
2108 ( Km[l] - Pip[l] ) * G[a][a] * G[j][j] * G[k][k] * G[l][l];
2109 }
2110 }
2111 }
2112 }
2113 tmp1 = B[0] * B[1] * B[2] * temp_PDF;
2114 temp_PDF = 0;
2115 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKsPip2, sKm, rRes2, 1, pro2 );
2116 if ( ispro_kspi2 == 0 )
2117 {
2118 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip2, sKs, sPip2, rRes2, 1,
2119 pro_kspi2 );
2120 ispro_kspi2 = 1;
2121 }
2122 Com_Multi( pro2, pro_kspi2, pro2b );
2123
2124 B[0] = barrier( 1, sKsPip2, sKs, sPip2, rRes2 );
2125 B[1] = barrier( 1, sKsKmPip2, sKsPip2, sKm, rRes2 );
2126 B[2] = barrier( 1, sD, sKsKmPip2, sPip1, rD2 );
2127
2128 for ( int a = 0; a < 4; a++ )
2129 {
2130 for ( int j = 0; j < 4; j++ )
2131 {
2132 for ( int k = 0; k < 4; k++ )
2133 {
2134 for ( int l = 0; l < 4; l++ )
2135 {
2136 temp_PDF += E[a][j][k][l] * pKsKmPip2[a] * ( pKsPip2[j] - Kp[j] ) * Pip[k] *
2137 ( Km[l] - Pi0[l] ) * G[a][a] * G[j][j] * G[k][k] * G[l][l];
2138 }
2139 }
2140 }
2141 }
2142 tmp2 = B[0] * B[1] * B[2] * temp_PDF;
2143 }
2144
2145 // text
2146 // Ds -> phi(1680)pi+ , phi(1680)->K-K*+ ,K*+ ->Km pi+
2147 if ( modetype[i] == 29 )
2148 {
2149 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKmPip1, sKm, rRes2, 1, pro1 );
2150 if ( ispro_kpi1 == 0 )
2151 {
2152 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip1, sKm, sPip1, rRes2, 1, pro_kpi1 );
2153 ispro_kpi1 = 1;
2154 }
2155 Com_Multi( pro1, pro_kpi1, pro1b );
2156
2157 B[0] = barrier( 1, sKmPip1, sKm, sPip1, rRes2 );
2158 B[1] = barrier( 1, sKsKmPip1, sKmPip1, sKs, rRes2 );
2159 B[2] = barrier( 1, sD, sKsKmPip1, sPip2, rD2 );
2160 for ( int a = 0; a < 4; a++ )
2161 {
2162 for ( int j = 0; j < 4; j++ )
2163 {
2164 for ( int k = 0; k < 4; k++ )
2165 {
2166 for ( int l = 0; l < 4; l++ )
2167 {
2168 temp_PDF += E[a][j][k][l] * pKsKmPip1[a] * ( pKmPip1[j] - Km[j] ) * Pi0[k] *
2169 ( Kp[l] - Pip[l] ) * G[a][a] * G[j][j] * G[k][k] * G[l][l];
2170 }
2171 }
2172 }
2173 }
2174 tmp1 = B[0] * B[1] * B[2] * temp_PDF;
2175
2176 temp_PDF = 0;
2177 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKmPip2, sKs, rRes2, 1, pro2 );
2178 if ( ispro_kpi2 == 0 )
2179 {
2180 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip2, sKm, sPip2, rRes2, 1, pro_kpi2 );
2181 ispro_kpi2 = 1;
2182 }
2183 Com_Multi( pro2, pro_kpi2, pro2b );
2184
2185 B[0] = barrier( 1, sKmPip2, sKm, sPip2, rRes2 );
2186 B[1] = barrier( 1, sKsKmPip2, sKmPip2, sKs, rRes2 );
2187 B[2] = barrier( 1, sD, sKsKmPip2, sPip1, rD2 );
2188
2189 for ( int a = 0; a < 4; a++ )
2190 {
2191 for ( int j = 0; j < 4; j++ )
2192 {
2193 for ( int k = 0; k < 4; k++ )
2194 {
2195 for ( int l = 0; l < 4; l++ )
2196 {
2197 temp_PDF += E[a][j][k][l] * pKsKmPip2[a] * ( pKmPip2[j] - Km[j] ) * Pip[k] *
2198 ( Kp[l] - Pi0[l] ) * G[a][a] * G[j][j] * G[k][k] * G[l][l];
2199 }
2200 }
2201 }
2202 }
2203
2204 tmp2 = B[0] * B[1] * B[2] * temp_PDF * ( -1 );
2205 }
2206
2207 //---------------------- for K1270*0 pi+ ,(K*+ -> KS pi+) ---------------------------
2208 if ( modetype[i] == 16 )
2209 {
2210 propagatorRBW( mass1sq, mass1[i], width1[i], sKsPipPip, sKsPip1, sPip2, rRes2, g2[i],
2211 pro1 );
2212 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip1, sKs, sPip1, rRes2, 1, pro_kpi1 );
2213 Com_Multi( pro1, pro_kpi1, pro1b );
2214
2215 B[0] = barrier( 1, sKsPip1, sKs, sPip1, rRes2 );
2216 B[1] = barrier( 1, sD, sKsPipPip, sKm, rD2 );
2217 B[2] = barrier( 2, sKsPipPip, sKsPip1, sPip2, rRes2 );
2218 if ( g2[i] == 0 )
2219 {
2220 for ( int a = 0; a < 4; a++ )
2221 {
2222 for ( int j = 0; j < 4; j++ )
2223 {
2224 temp_PDF += t1D[a] * ( pKsPipPip[a] * pKsPipPip[j] / sKsPipPip - G[a][j] ) *
2225 t1KsPip1[j] * G[a][a] * G[j][j];
2226 }
2227 }
2228 tmp1 = B[0] * B[1] * temp_PDF;
2229 }
2230 if ( g2[i] == 2 )
2231 {
2232 for ( int a = 0; a < 4; a++ )
2233 {
2234 for ( int j = 0; j < 4; j++ )
2235 { temp_PDF += t1D[a] * t2KsPipPip_kspip[a][j] * t1KsPip1[j] * G[a][a] * G[j][j]; }
2236 }
2237 tmp1 = B[0] * B[1] * B[2] * temp_PDF;
2238 }
2239
2240 temp_PDF = 0;
2241 propagatorRBW( mass1sq, mass1[i], width1[i], sKsPipPip, sKsPip2, sPip1, rRes2, g2[i],
2242 pro1 );
2243 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip2, sKs, sPip2, rRes2, 1, pro_kpi2 );
2244 Com_Multi( pro1, pro_kpi2, pro2b );
2245
2246 B[0] = barrier( 1, sKsPip2, sKs, sPip2, rRes2 );
2247 B[1] = barrier( 1, sD, sKsPipPip, sKm, rD2 );
2248 B[2] = barrier( 2, sKsPipPip, sKsPip2, sPip1, rRes2 );
2249 if ( g2[i] == 0 )
2250 {
2251 for ( int a = 0; a < 4; a++ )
2252 {
2253 for ( int j = 0; j < 4; j++ )
2254 {
2255 temp_PDF += t1D[a] * ( pKsPipPip[a] * pKsPipPip[j] / sKsPipPip - G[a][j] ) *
2256 t1KsPip2[j] * G[a][a] * G[j][j];
2257 }
2258 }
2259 tmp2 = B[0] * B[1] * temp_PDF;
2260 }
2261 if ( g2[i] == 2 )
2262 {
2263 for ( int a = 0; a < 4; a++ )
2264 {
2265 for ( int j = 0; j < 4; j++ )
2266 { temp_PDF += t1D[a] * t2KsPipPip_kspip2[a][j] * t1KsPip2[j] * G[a][a] * G[j][j]; }
2267 }
2268 tmp2 = B[0] * B[1] * B[2] * temp_PDF;
2269 }
2270 }
2271 // K*1270tokpiswave:36 is kmpiswave
2272 if ( modetype[i] == 36 )
2273 {
2274 propagatorRBW( mass1sq, mass1[i], width1[i], sKmPipPip, sKmPip1, sPip2, rRes2, 1, pro1 );
2275 if ( ispro_kspi1sw == 0 )
2276 {
2277 KPiSLASS( sKmPip1, sKm, sPip1, pro_kspi1sw );
2278 ispro_kspi1sw = 1;
2279 }
2280 Com_Multi( pro1, pro_kspi1sw, pro1b );
2281 calt1( pKmPipPip, Km, t1D );
2282
2283 if ( barr_kskpi1_kspi < 0.0 )
2284 barr_kskpi1_kspi = barrier( 1, sKmPipPip, sKmPip1, sPip2, rRes2 );
2285 if ( barr_ds_kskpi1 < 0.0 ) barr_ds_kskpi1 = barrier( 1, sD, sKmPipPip, sKs, rD2 );
2286 for ( int a = 0; a < 4; a++ )
2287 {
2288 temp_PDF += G[a][a] * t1D[a] * t1KmPipPip_kmpip[a];
2289 // temp_PDF += G[a][a]*t1D[a]*t1KsKmPip1_kspip[a];
2290 }
2291 tmp1 = barr_kskpi1_kspi * barr_ds_kskpi1 * temp_PDF;
2292
2293 temp_PDF = 0;
2294 propagatorRBW( mass1sq, mass1[i], width1[i], sKmPipPip, sKmPip2, sPip1, rRes2, 1, pro1 );
2295 if ( ispro_kspi2sw == 0 )
2296 {
2297 KPiSLASS( sKmPip2, sKm, sPip2, pro_kspi2sw );
2298 ispro_kspi2sw = 1;
2299 }
2300 Com_Multi( pro1, pro_kspi2sw, pro2b );
2301 calt1( pKmPipPip, Km, t1D );
2302
2303 if ( barr_kskpi2_kspi < 0.0 )
2304 barr_kskpi2_kspi = barrier( 1, sKmPipPip, sKmPip2, sPip1, rRes2 );
2305 if ( barr_ds_kskpi2 < 0.0 ) barr_ds_kskpi2 = barrier( 1, sD, sKmPipPip, sKs, rD2 );
2306 for ( int a = 0; a < 4; a++ )
2307 {
2308 temp_PDF += G[a][a] * t1D[a] * t1KmPipPip_kmpip2[a];
2309 // temp_PDF += G[a][a]*t1D[a]*t1KsKmPip1_kspip[a];
2310 }
2311 tmp2 = barr_kskpi2_kspi * barr_ds_kskpi2 * temp_PDF;
2312 }
2313
2314 // K*1270tokpiswave:36 is kmpiswave
2315 if ( modetype[i] == 46 )
2316 {
2317 propagatorRBW( mass1sq, mass1[i], width1[i], sKsPipPip, sKsPip1, sPip2, rRes2, 1, pro1 );
2318 if ( ispro_kspi1sw == 0 )
2319 {
2320 KPiSLASS( sKsPip1, sKs, sPip1, pro_kspi1sw );
2321 ispro_kspi1sw = 1;
2322 }
2323 Com_Multi( pro1, pro_kspi1sw, pro1b );
2324 calt1( pKsPipPip, Kp, t1D );
2325
2326 if ( barr_kskpi1_kspi < 0.0 )
2327 barr_kskpi1_kspi = barrier( 1, sKsPipPip, sKsPip1, sPip2, rRes2 );
2328 if ( barr_ds_kskpi1 < 0.0 ) barr_ds_kskpi1 = barrier( 1, sD, sKsPipPip, sKm, rD2 );
2329 for ( int a = 0; a < 4; a++ )
2330 {
2331 temp_PDF += G[a][a] * t1D[a] * t1KsPipPip_kspip[a];
2332 // temp_PDF += G[a][a]*t1D[a]*t1KsKmPip1_kspip[a];
2333 }
2334 tmp1 = barr_kskpi1_kspi * barr_ds_kskpi1 * temp_PDF;
2335 temp_PDF = 0;
2336 propagatorRBW( mass1sq, mass1[i], width1[i], sKsPipPip, sKsPip2, sPip1, rRes2, 1, pro1 );
2337 if ( ispro_kspi2sw == 0 )
2338 {
2339 KPiSLASS( sKsPip2, sKs, sPip2, pro_kspi2sw );
2340 ispro_kspi2sw = 1;
2341 }
2342 Com_Multi( pro1, pro_kspi2sw, pro2b );
2343 calt1( pKsPipPip, Kp, t1D );
2344
2345 if ( barr_kskpi2_kspi < 0.0 )
2346 barr_kskpi2_kspi = barrier( 1, sKsPipPip, sKsPip2, sPip1, rRes2 );
2347 if ( barr_ds_kskpi2 < 0.0 ) barr_ds_kskpi2 = barrier( 1, sD, sKsPipPip, sKm, rD2 );
2348 for ( int a = 0; a < 4; a++ )
2349 {
2350 temp_PDF += G[a][a] * t1D[a] * t1KsPipPip_kspip2[a];
2351 // temp_PDF += G[a][a]*t1D[a]*t1KsKmPip1_kspip[a];
2352 }
2353 tmp2 = barr_kskpi2_kspi * barr_ds_kskpi2 * temp_PDF;
2354 }
2355
2356 // text
2357 /// f1420->kpi_sw
2358 if ( modetype[i] == 46 )
2359 {
2360 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKmPip1, sKs, rRes2, 1, pro1 );
2361 if ( ispro_kspi1sw == 0 )
2362 {
2363 KPiSLASS( sKsPip1, sKm, sPip1, pro_kpi1sw );
2364 ispro_kpi1sw = 1;
2365 }
2366 Com_Multi( pro1, pro_kpi1sw, pro1b );
2367 calt1( pKsKmPip1, Pi0, t1D );
2368
2369 barr_kskpi1_kpi = barrier( 1, sKsKmPip1, sKmPip1, sKs, rRes2 );
2370 if ( barr_ds_kskpi1 < 0.0 ) barr_ds_kskpi1 = barrier( 1, sD, sKsKmPip1, sPip2, rD2 );
2371 for ( int a = 0; a < 4; a++ ) { temp_PDF += G[a][a] * t1D[a] * t1KsKmPip1_kmpip[a]; }
2372 tmp1 = barr_kskpi1_kpi * barr_ds_kskpi1 * temp_PDF;
2373
2374 temp_PDF = 0;
2375 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKmPip2, sKs, rRes2, 1, pro1 );
2376 if ( ispro_kpi2sw == 0 )
2377 {
2378 KPiSLASS( sKsPip2, sKs, sPip2, pro_kpi2sw );
2379 ispro_kpi2sw = 1;
2380 }
2381 Com_Multi( pro1, pro_kpi2sw, pro2b );
2382 calt1( pKsKmPip2, Pip, t1D );
2383
2384 barr_kskpi2_kpi = barrier( 1, sKsKmPip2, sKmPip2, sKs, rRes2 );
2385 if ( barr_ds_kskpi2 < 0.0 ) barr_ds_kskpi2 = barrier( 1, sD, sKsKmPip2, sPip1, rD2 );
2386 for ( int a = 0; a < 4; a++ ) { temp_PDF += G[a][a] * t1D[a] * t1KsKmPip2_kmpip[a]; }
2387 tmp2 = barr_kskpi2_kpi * barr_ds_kskpi2 * temp_PDF;
2388 }
2389
2390 // 1+ Ds -> K1270 Ks , K1270 -> pi+ K*0 ,K*0 ->K+ pi-
2391 if ( modetype[i] == 6 )
2392 {
2393 propagatorRBW( mass1sq, mass1[i], width1[i], sKmPipPip, sKmPip1, sPip2, rRes2, g2[i],
2394 pro1 );
2395 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip1, sKm, sPip1, rRes2, 1, pro_kpi1 );
2396 Com_Multi( pro1, pro_kpi1, pro1b );
2397
2398 B[0] = barrier( 1, sKmPip1, sKm, sPip1, rRes2 );
2399 B[1] = barrier( 1, sD, sKmPipPip, sKs, rD2 );
2400 B[2] = barrier( 2, sKmPipPip, sKmPip1, sPip2, rRes2 );
2401 if ( g2[i] == 0 )
2402 {
2403 for ( int a = 0; a < 4; a++ )
2404 {
2405 for ( int j = 0; j < 4; j++ )
2406 {
2407 temp_PDF += t1D[a] * ( pKmPipPip[a] * pKmPipPip[j] / sKmPipPip - G[a][j] ) *
2408 t1KmPip1[j] * G[a][a] * G[j][j];
2409 }
2410 }
2411 tmp1 = B[0] * B[1] * temp_PDF;
2412 }
2413 if ( g2[i] == 2 )
2414 {
2415 for ( int a = 0; a < 4; a++ )
2416 {
2417 for ( int j = 0; j < 4; j++ )
2418 { temp_PDF += t1D[a] * t2KmPipPip_kmpip[a][j] * t1KmPip1[j] * G[a][a] * G[j][j]; }
2419 }
2420 tmp1 = B[0] * B[1] * B[2] * temp_PDF;
2421 }
2422
2423 temp_PDF = 0;
2424 propagatorRBW( mass1sq, mass1[i], width1[i], sKmPipPip, sKmPip2, sPip1, rRes2, g2[i],
2425 pro1 );
2426 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip2, sKm, sPip2, rRes2, 1, pro_kpi2 );
2427 Com_Multi( pro1, pro_kpi2, pro2b );
2428
2429 B[0] = barrier( 1, sKmPip2, sKm, sPip2, rRes2 );
2430 B[1] = barrier( 1, sD, sKmPipPip, sKs, rD2 );
2431 B[2] = barrier( 2, sKmPipPip, sKmPip2, sPip1, rRes2 );
2432 if ( g2[i] == 0 )
2433 {
2434 for ( int a = 0; a < 4; a++ )
2435 {
2436 for ( int j = 0; j < 4; j++ )
2437 {
2438 temp_PDF += t1D[a] * ( pKmPipPip[a] * pKmPipPip[j] / sKmPipPip - G[a][j] ) *
2439 t1KmPip2[j] * G[a][a] * G[j][j];
2440 }
2441 }
2442 tmp2 = B[0] * B[1] * temp_PDF;
2443 }
2444 if ( g2[i] == 2 )
2445 {
2446 for ( int a = 0; a < 4; a++ )
2447 {
2448 for ( int j = 0; j < 4; j++ )
2449 { temp_PDF += t1D[a] * t2KmPipPip_kmpip2[a][j] * t1KmPip2[j] * G[a][a] * G[j][j]; }
2450 }
2451 tmp2 = B[0] * B[1] * B[2] * temp_PDF;
2452 }
2453 }
2454
2455 /// Kstp Kmpi 1430
2456 if ( modetype[i] == 30 )
2457 {
2458 if ( ispro_kspi1 == 0 )
2459 {
2460 propagatorRBW( mass1sq, mass1[i], width1[i], sKsPip1, sKs, sPip1, rRes2, 1,
2461 pro_kspi1 );
2462 ispro_kspi1 = 1;
2463 }
2464 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip2, sKm, sPip2, rRes2, 0, pro_kpi2_0 );
2465
2466 Com_Multi( pro_kspi1, pro_kpi2_0, pro1b );
2467 calt1( pKsPip1, pKmPip2, t1D );
2468
2469 if ( barr_kspi1 < 0.0 ) barr_kspi1 = barrier( 1, sKsPip1, sKs, sPip1, rRes2 );
2470 if ( barr_kspi1_kpi2 < 0.0 ) barr_kspi1_kpi2 = barrier( 1, sD, sKsPip1, sKmPip2, rD2 );
2471 for ( int a = 0; a < 4; a++ ) { temp_PDF += t1D[a] * t1KsPip1[a] * G[a][a]; }
2472 tmp1 = barr_kspi1 * barr_kspi1_kpi2 * temp_PDF;
2473
2474 temp_PDF = 0;
2475 if ( ispro_kspi2 == 0 )
2476 {
2477 propagatorRBW( mass1sq, mass1[i], width1[i], sKsPip2, sKs, sPip2, rRes2, 1,
2478 pro_kspi2 );
2479 ispro_kspi2 = 1;
2480 }
2481
2482 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip1, sKm, sPip1, rRes2, 0, pro_kpi1_0 );
2483 Com_Multi( pro_kspi2, pro_kpi1_0, pro2b );
2484
2485 calt1( pKsPip2, pKmPip1, t1D );
2486 if ( barr_kspi2 < 0.0 ) barr_kspi2 = barrier( 1, sKsPip2, sKs, sPip2, rRes2 );
2487 if ( barr_kspi2_kpi1 < 0.0 ) barr_kspi2_kpi1 = barrier( 1, sD, sKsPip2, sKmPip1, rD2 );
2488 for ( int a = 0; a < 4; a++ ) { temp_PDF += t1D[a] * t1KsPip2[a] * G[a][a]; }
2489 tmp2 = barr_kspi2 * barr_kspi2_kpi1 * temp_PDF;
2490 }
2491
2492 /// Kst0 Kspi 1430
2493 if ( modetype[i] == 31 )
2494 {
2495 if ( ispro_kpi2 == 0 )
2496 {
2497 propagatorRBW( mass1sq, mass1[i], width1[i], sKmPip2, sKm, sPip2, rRes2, 1, pro_kpi2 );
2498 ispro_kpi2 = 1;
2499 }
2500 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip1, sKs, sPip1, rRes2, 0,
2501 pro_kspi1_0 );
2502 Com_Multi( pro_kpi2, pro_kspi1_0, pro1b );
2503
2504 calt1( pKsPip1, pKmPip2, t1D );
2505 if ( barr_kpi2 < 0.0 ) barr_kpi2 = barrier( 1, sKmPip2, sKm, sPip2, rRes2 );
2506 if ( barr_kspi1_kpi2 < 0.0 ) barr_kspi1_kpi2 = barrier( 1, sD, sKsPip1, sKmPip2, rD2 );
2507 for ( int a = 0; a < 4; a++ ) { temp_PDF += t1D[a] * t1KmPip2[a] * G[a][a]; }
2508 tmp1 = barr_kpi2 * barr_kspi1_kpi2 * temp_PDF;
2509
2510 temp_PDF = 0;
2511 if ( ispro_kpi1 == 0 )
2512 {
2513 propagatorRBW( mass1sq, mass1[i], width1[i], sKmPip1, sKm, sPip1, rRes2, 1, pro_kpi1 );
2514 ispro_kpi1 = 1;
2515 }
2516 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip2, sKs, sPip2, rRes2, 0,
2517 pro_kspi2_0 );
2518 Com_Multi( pro_kpi1, pro_kspi2_0, pro2b );
2519
2520 calt1( pKsPip2, pKmPip1, t1D );
2521 if ( barr_kpi1 < 0.0 ) barr_kpi1 = barrier( 1, sKmPip1, sKm, sPip1, rRes2 );
2522 if ( barr_kspi2_kpi1 < 0.0 ) barr_kspi2_kpi1 = barrier( 1, sD, sKsPip2, sKmPip1, rD2 );
2523 for ( int a = 0; a < 4; a++ ) { temp_PDF += t1D[a] * t1KmPip1[a] * G[a][a]; }
2524 tmp2 = barr_kspi2_kpi1 * barr_kpi1 * temp_PDF;
2525 }
2526
2527 // Kspi 1430 kpi 1430
2528 if ( modetype[i] == 32 )
2529 {
2530 propagatorRBW( mass1sq, mass1[i], width1[i], sKmPip1, sKm, sPip1, rRes2, 0, pro_kpi1_0 );
2531 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip2, sKs, sPip2, rRes2, 0,
2532 pro_kspi2_0 );
2533 Com_Multi( pro_kpi1_0, pro_kpi2_0, pro1b );
2534 tmp1 = 1;
2535 temp_PDF = 0;
2536 propagatorRBW( mass1sq, mass1[i], width1[i], sKmPip2, sKm, sPip2, rRes2, 0, pro_kpi2_0 );
2537 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip1, sKs, sPip1, rRes2, 0,
2538 pro_kspi1_0 );
2539 Com_Multi( pro_kspi2_0, pro_kpi1_0, pro2b );
2540 tmp2 = 1;
2541 }
2542
2543 if ( modetype[i] == 1000 )
2544 {
2545 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKsPip1, sKm, rRes2, 1, pro1 );
2546 pro1[0] = 1;
2547 pro1[1] = 0;
2548 if ( ispro_kspi1 == 0 )
2549 {
2550 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip1, sKs, sPip1, rRes2, 1,
2551 pro_kspi1 );
2552 ispro_kspi1 = 1;
2553 }
2554 Com_Multi( pro1, pro_kspi1, pro1b );
2555 // printf(" %.15f %.15f %.15f\n",B[0],B[1],B[2]);
2556
2557 B[0] = barrier( 1, sKsPip1, sKs, sPip1, rRes2 );
2558 B[1] = barrier( 1, sKsKmPip1, sKsPip1, sKm, rRes2 );
2559 B[2] = barrier( 1, sD, sKsKmPip1, sPip2, rD2 );
2560 for ( int a = 0; a < 4; a++ )
2561 {
2562 for ( int j = 0; j < 4; j++ )
2563 {
2564 for ( int k = 0; k < 4; k++ )
2565 {
2566 for ( int l = 0; l < 4; l++ )
2567 {
2568 temp_PDF += E[a][j][k][l] * pKsKmPip1[a] * ( pKsPip1[j] - Kp[j] ) * Pi0[k] *
2569 ( Km[l] - Pip[l] ) * G[a][a] * G[j][j] * G[k][k] * G[l][l];
2570 }
2571 }
2572 }
2573 }
2574 tmp1 = B[0] * B[1] * B[2] * temp_PDF;
2575
2576 temp_PDF = 0;
2577 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKsPip2, sKm, rRes2, 1, pro2 );
2578
2579 pro2[0] = 1;
2580 pro2[1] = 0;
2581 if ( ispro_kspi2 == 0 )
2582 {
2583 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip2, sKs, sPip2, rRes2, 1,
2584 pro_kspi2 );
2585 ispro_kspi2 = 1;
2586 }
2587 Com_Multi( pro2, pro_kspi2, pro2b );
2588
2589 B[0] = barrier( 1, sKsPip2, sKs, sPip2, rRes2 );
2590 B[1] = barrier( 1, sKsKmPip2, sKsPip2, sKm, rRes2 );
2591 B[2] = barrier( 1, sD, sKsKmPip2, sPip1, rD2 );
2592
2593 for ( int a = 0; a < 4; a++ )
2594 {
2595 for ( int j = 0; j < 4; j++ )
2596 {
2597 for ( int k = 0; k < 4; k++ )
2598 {
2599 for ( int l = 0; l < 4; l++ )
2600 {
2601 temp_PDF += E[a][j][k][l] * pKsKmPip2[a] * ( pKsPip2[j] - Kp[j] ) * Pip[k] *
2602 ( Km[l] - Pi0[l] ) * G[a][a] * G[j][j] * G[k][k] * G[l][l];
2603 }
2604 }
2605 }
2606 }
2607 tmp2 = B[0] * B[1] * B[2] * temp_PDF;
2608 }
2609
2610 // text
2611 // Ds -> phi(1680)pi+ , phi(1680)->K-K*+ ,K*+ ->Km pi+
2612 if ( modetype[i] == 2000 )
2613 {
2614 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKmPip1, sKm, rRes2, 1, pro1 );
2615 if ( ispro_kpi1 == 0 )
2616 {
2617
2618 pro1[0] = 1;
2619 pro1[1] = 0;
2620 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip1, sKm, sPip1, rRes2, 1, pro_kpi1 );
2621 ispro_kpi1 = 1;
2622 }
2623 Com_Multi( pro1, pro_kpi1, pro1b );
2624
2625 B[0] = barrier( 1, sKmPip1, sKm, sPip1, rRes2 );
2626 B[1] = barrier( 1, sKsKmPip1, sKmPip1, sKs, rRes2 );
2627 B[2] = barrier( 1, sD, sKsKmPip1, sPip2, rD2 );
2628 for ( int a = 0; a < 4; a++ )
2629 {
2630 for ( int j = 0; j < 4; j++ )
2631 {
2632 for ( int k = 0; k < 4; k++ )
2633 {
2634 for ( int l = 0; l < 4; l++ )
2635 {
2636 temp_PDF += E[a][j][k][l] * pKsKmPip1[a] * ( pKmPip1[j] - Km[j] ) * Pi0[k] *
2637 ( Kp[l] - Pip[l] ) * G[a][a] * G[j][j] * G[k][k] * G[l][l];
2638 }
2639 }
2640 }
2641 }
2642 tmp1 = B[0] * B[1] * B[2] * temp_PDF;
2643
2644 temp_PDF = 0;
2645 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKmPip2, sKs, rRes2, 1, pro2 );
2646 pro2[0] = 1;
2647 pro2[1] = 0;
2648 if ( ispro_kpi2 == 0 )
2649 {
2650 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip2, sKm, sPip2, rRes2, 1, pro_kpi2 );
2651 ispro_kpi2 = 1;
2652 }
2653 Com_Multi( pro2, pro_kpi2, pro2b );
2654
2655 B[0] = barrier( 1, sKmPip2, sKm, sPip2, rRes2 );
2656 B[1] = barrier( 1, sKsKmPip2, sKmPip2, sKs, rRes2 );
2657 B[2] = barrier( 1, sD, sKsKmPip2, sPip1, rD2 );
2658
2659 for ( int a = 0; a < 4; a++ )
2660 {
2661 for ( int j = 0; j < 4; j++ )
2662 {
2663 for ( int k = 0; k < 4; k++ )
2664 {
2665 for ( int l = 0; l < 4; l++ )
2666 {
2667 temp_PDF += E[a][j][k][l] * pKsKmPip2[a] * ( pKmPip2[j] - Km[j] ) * Pip[k] *
2668 ( Kp[l] - Pi0[l] ) * G[a][a] * G[j][j] * G[k][k] * G[l][l];
2669 }
2670 }
2671 }
2672 }
2673
2674 tmp2 = B[0] * B[1] * B[2] * temp_PDF * ( -1 );
2675 }
2676
2677 if ( modetype[i] == 1001 )
2678 {
2679 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKmPip1, sKs, rRes2, 1, pro1 );
2680 pro1[0] = 1;
2681 pro1[1] = 0;
2682 if ( ispro_kpi1 == 0 )
2683 {
2684 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip1, sKm, sPip1, rRes2, 1, pro_kpi1 );
2685 ispro_kpi1 = 1;
2686 }
2687 Com_Multi( pro1, pro_kpi1, pro1b );
2688
2689 if ( barr_kpi1 < 0.0 ) barr_kpi1 = barrier( 1, sKmPip1, sKm, sPip1, rRes2 );
2690 barr_kskpi1_kpi = barrier( 1, sKsKmPip1, sKmPip1, sKs, rRes2 );
2691 for ( int a = 0; a < 4; a++ ) { temp_PDF += Km[a] * t1KmPip1[a] * G[a][a]; }
2692 tmp1 = barr_kpi1 * barr_kskpi1_kpi * temp_PDF;
2693
2694 temp_PDF = 0;
2695 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKmPip2, sKs, rRes2, 1, pro1 );
2696 pro1[0] = 1;
2697 pro1[1] = 0;
2698 if ( ispro_kpi2 == 0 )
2699 {
2700 propagatorRBW( mass2sq, mass2[i], width2[i], sKmPip2, sKm, sPip2, rRes2, 1, pro_kpi2 );
2701 ispro_kpi2 = 1;
2702 }
2703 Com_Multi( pro1, pro_kpi2, pro2b );
2704
2705 if ( barr_kpi2 < 0.0 ) barr_kpi2 = barrier( 1, sKmPip2, sKm, sPip2, rRes2 );
2706 barr_kskpi2_kpi = barrier( 1, sKsKmPip2, sKmPip2, sKs, rRes2 );
2707 for ( int a = 0; a < 4; a++ ) { temp_PDF += Km[a] * t1KmPip2[a] * G[a][a]; }
2708 tmp2 = barr_kpi2 * barr_kskpi2_kpi * temp_PDF;
2709 }
2710
2711 // Ds -> eta(1475)pi+ eta(1475)->K*Km,K*+>Ks pi+
2712 if ( modetype[i] == 2001 )
2713 {
2714 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip1, sKsPip1, sKm, rRes2, 1, pro1 );
2715 pro1[0] = 1;
2716 pro1[1] = 0;
2717
2718 if ( ispro_kspi1 == 0 )
2719 {
2720 propagatorRBW( mass2sq, mass2[i], width2[i], sKsPip1, sKs, sPip1, rRes2, 1,
2721 pro_kspi1 );
2722 ispro_kspi1 = 1;
2723 }
2724 Com_Multi( pro1, pro_kspi1, pro1b );
2725
2726 if ( barr_kspi1 < 0.0 ) barr_kspi1 = barrier( 1, sKsPip1, sKs, sPip1, rRes2 );
2727 if ( barr_kskpi1_kspi < 0.0 )
2728 barr_kskpi1_kspi = barrier( 1, sKsKmPip1, sKsPip1, sKm, rRes2 );
2729 for ( int a = 0; a < 4; a++ ) { temp_PDF += Kp[a] * t1KsPip1[a] * G[a][a]; }
2730 tmp1 = barr_kspi1 * barr_kskpi1_kspi * temp_PDF;
2731
2732 temp_PDF = 0;
2733 propagatorRBW( mass1sq, mass1[i], width1[i], sKsKmPip2, sKsPip2, sKm, rRes2, 1, pro1 );
2734 pro1[0] = 1;
2735 pro1[1] = 0;
2736 if ( ispro_kspi2 == 0 )
2737 {
2738 propagatorRBW( mass1sq, mass2[i], width2[i], sKsPip2, sKs, sPip2, rRes2, 1,
2739 pro_kspi2 );
2740 ispro_kspi2 = 1;
2741 }
2742 Com_Multi( pro1, pro_kspi2, pro2b );
2743
2744 if ( barr_kspi2 < 0.0 ) barr_kspi2 = barrier( 1, sKsPip2, sKs, sPip2, rRes2 );
2745 if ( barr_kskpi2_kspi < 0.0 )
2746 barr_kskpi2_kspi = barrier( 1, sKsKmPip2, sKsPip2, sKm, rRes2 );
2747 for ( int a = 0; a < 4; a++ ) { temp_PDF += Kp[a] * t1KsPip2[a] * G[a][a]; }
2748 tmp2 = barr_kspi2 * barr_kskpi2_kspi * temp_PDF;
2749 }
2750 if ( modetype[i] == 1100 || modetype[i] == 1200 || modetype[i] == 1500 ||
2751 modetype[i] == 1001 || modetype[i] == 2001 || modetype[i] == 1000 ||
2752 modetype[i] == 2000 || modetype[i] == 290 || modetype[i] == 280 ||
2753 modetype[i] == 244 || modetype[i] == 888 || modetype[i] == 133 ||
2754 modetype[i] == 200 || modetype[i] == 99 || modetype[i] == 36 || modetype[i] == 46 ||
2755 modetype[i] == 16 || modetype[i] == 2 || modetype[i] == 21 || modetype[i] == 22 ||
2756 modetype[i] == 8 || modetype[i] == 4 || modetype[i] == 14 || modetype[i] == 5 ||
2757 modetype[i] == 17 || modetype[i] == 3 || modetype[i] == 7 || modetype[i] == 10 ||
2758 modetype[i] == 9 || modetype[i] == 13 || modetype[i] == 20 || modetype[i] == 11 ||
2759 modetype[i] == 12 || modetype[i] == 15 || modetype[i] == 19 || modetype[i] == 26 ||
2760 modetype[i] == 24 || modetype[i] == 28 || modetype[i] == 29 || modetype[i] == 27 ||
2761 modetype[i] == 6 || modetype[i] == 1 || modetype[i] == 30 || modetype[i] == 31 ||
2762 modetype[i] == 32 )
2763 {
2764 amp_tmp[0] = tmp1 * pro1b[0] + tmp2 * pro2b[0];
2765 amp_tmp[1] = tmp1 * pro1b[1] + tmp2 * pro2b[1];
2766 }
2767 if ( modetype[i] == 1111 )
2768 {
2769 amp_tmp[0] = 1;
2770 amp_tmp[1] = 0;
2771 }
2772
2773 Com_Multi( amp_tmp, cof, amp_PDF );
2774 PDF[0] += amp_PDF[0];
2775 PDF[1] += amp_PDF[1];
2776 }
2777 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
2778 if ( value <= 0 ) value = 0.00000000001;
2779 Result = value;
2780 return Result;
2781}
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
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:209
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
EvtDToKSKmpippip::decay
void decay(EvtParticle *p)
Definition EvtDToKSKmpippip.cc:139
EvtDToKSKmpippip::clone
EvtDecayBase * clone()
Definition EvtDToKSKmpippip.cc:38
EvtDToKSKmpippip::~EvtDToKSKmpippip
virtual ~EvtDToKSKmpippip()
Definition EvtDToKSKmpippip.cc:34
EvtDToKSKmpippip::getName
void getName(std::string &name)
Definition EvtDToKSKmpippip.cc:36
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
EvtPDL::getStdHep
static int getStdHep(EvtId id)
Definition EvtPDL.hh:61
EvtParticle::getId
EvtId getId() const
Definition EvtParticle.cc:100
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