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