BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToKSKpPipPim.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: EvtDsToKSKpPipPim.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/EvtDsToKSKpPipPim.hh"
27#include <string>
28#include "EvtGenBase/EvtConst.hh"
29#include "EvtGenBase/EvtFlatte.hh"
30#include "EvtGenBase/EvtDecayTable.hh"
31using std::endl;
32
33EvtDsToKSKpPipPim::~EvtDsToKSKpPipPim() {}
34
35void EvtDsToKSKpPipPim::getName(std::string& model_name){
36 model_name="DsToKSKpPipPim";
37}
38
42
43void EvtDsToKSKpPipPim::init(){
44 // check that there are 0 arguments
45 checkNArg(0);
46 checkNDaug(4);
47
48 checkSpinParent(EvtSpinType::SCALAR);
49 checkSpinDaughter(0,EvtSpinType::SCALAR);
50 checkSpinDaughter(1,EvtSpinType::SCALAR);
51 checkSpinDaughter(2,EvtSpinType::SCALAR);
52 checkSpinDaughter(3,EvtSpinType::SCALAR);
53
54 int mode = 0;
55 phi[1] = -3.1173e-03;
56 phi[2] = 5.6119e+00;
57 phi[3] = -7.0591e-01;
58 phi[4] = -3.1220e+00;
59 phi[5] = 3.5914e+00;
60 phi[6] = 3.8510e+00;
61 phi[7] = -1.0830e+00;
62
63 rho[1] = 1.8516e+00;
64 rho[2] = 2.8585e+00;
65 rho[3] = 1.1065e+00;
66 rho[4] = 5.4084e-01;
67 rho[5] = 1.3816e+00;
68 rho[6] = 1.9177e+00;
69 rho[7] = 2.2119e+00;
70 modetype[0]= 7;
71 modetype[1]= 4;
72 modetype[2]= 3;
73 modetype[3]= 15;
74 modetype[4]= 8;
75 modetype[5]= 11;
76 modetype[6]= 12;
77 modetype[7]= 17;
78 phi[0] = 0;
79 rho[0] = 1;
80
81 //cout << "DsToKSKpPipPim :" << endl;
82 //for (int i=0; i<8; i++) {
83 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
84 //}
85
86 mD = 1.86486;
87 metap = 0.95778;
88 mk0 = 0.497614;
89 mass_Kaon = 0.49368;
90 mass_Eta = 0.547862;
91 mass_Pion = 0.13957;
92 math_pi = 3.1415926;
93 mass_Pion2 = 0.0194797849;
94 mass_2Pion = 0.27914;
95 math_2pi = 6.2831852;
96
97 rD2 = 25.0;
98 rRes2 = 9.0;
99 GK1270 = 8.7000e-02;
100 GK1400 = 1.7400e-01;
101 GKst0 = 4.7300e-02;
102 GKstp = 5.0300e-02;
103 GP1680 = 1.5000e-01;
104 GRho = 1.4780e-01;
105 Ga0_980 = 5.0000e-02;
106 Geta1475= 9.0000e-02;
107 mK1270 = 1.2720e+00;
108 mK1400 = 1.4030e+00;
109 mKst0 = 8.9555e-01;
110 mKstp = 8.9176e-01;
111 mP1680 = 1.6800e+00;
112 mRho = 7.7526e-01;
113 ma0_980 = 9.9000e-01;
114 meta1475= 1.4750e+00;
115
116 GS1 = 0.636619783;
117 GS2 = 0.01860182466;
118 GS3 = 0.1591549458; // 1/(2*math_2pi)
119 GS4 = 0.00620060822; // mass_Pion2/math_pi
120
121 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
122 int EE[4][4][4][4] =
123 { { {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
124 {{0,0,0,0}, {0,0,0,0}, {0,0,0,1}, {0,0,-1,0}},
125 {{0,0,0,0}, {0,0,0,-1}, {0,0,0,0}, {0,1,0,0} },
126 {{0,0,0,0}, {0,0,1,0}, {0,-1,0,0}, {0,0,0,0} } },
127 { {{0,0,0,0}, {0,0,0,0}, {0,0,0,-1}, {0,0,1,0} },
128 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
129 {{0,0,0,1}, {0,0,0,0}, {0,0,0,0}, {-1,0,0,0}},
130 {{0,0,-1,0}, {0,0,0,0}, {1,0,0,0}, {0,0,0,0} } },
131 { {{0,0,0,0}, {0,0,0,1}, {0,0,0,0}, {0,-1,0,0}},
132 {{0,0,0,-1}, {0,0,0,0}, {0,0,0,0}, {1,0,0,0} },
133 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
134 {{0,1,0,0}, {-1,0,0,0}, {0,0,0,0}, {0,0,0,0} } },
135 { {{0,0,0,0}, {0,0,-1,0}, {0,1,0,0}, {0,0,0,0} },
136 {{0,0,1,0}, {0,0,0,0}, {-1,0,0,0}, {0,0,0,0} },
137 {{0,-1,0,0}, {1,0,0,0}, {0,0,0,0}, {0,0,0,0} },
138 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} } } };
139 for (int i=0; i<4; i++) {
140 for (int j=0; j<4; j++) {
141 G[i][j] = GG[i][j];
142 for (int k=0; k<4; k++) {
143 for (int l=0; l<4; l++) {
144 E[i][j][k][l] = EE[i][j][k][l];
145 }
146 }
147 }
148 }
149
150}
151
155
156void EvtDsToKSKpPipPim::decay( EvtParticle *p){
157/*
158 double maxprob=0;
159 for(int ir=0;ir<=60000000;ir++){
160 p->initializePhaseSpace(getNDaug(),getDaugs());
161 EvtVector4R _ks = p->getDaug(0)->getP4();
162 EvtVector4R _kp = p->getDaug(1)->getP4();
163 EvtVector4R _pip = p->getDaug(2)->getP4();
164 EvtVector4R _pim = p->getDaug(3)->getP4();
165
166 double _Pip[4],_Pim[4],_Ks[4],_Kp[4];
167 _Ks[0] = _ks.get(0); _Kp[0] = _kp.get(0); _Pip[0] = _pip.get(0); _Pim[0] = _pim.get(0);
168 _Ks[1] = _ks.get(1); _Kp[1] = _kp.get(1); _Pip[1] = _pip.get(1); _Pim[1] = _pim.get(1);
169 _Ks[2] = _ks.get(2); _Kp[2] = _kp.get(2); _Pip[2] = _pip.get(2); _Pim[2] = _pim.get(2);
170 _Ks[3] = _ks.get(3); _Kp[3] = _kp.get(3); _Pip[3] = _pip.get(3); _Pim[3] = _pim.get(3);
171
172 double _prob;
173 int g0[8]={1,1,1,1,1,1,1,1};
174 int nstates=10;
175
176 calEvaMy(_Ks, _Kp, _Pip, _Pim, mass, width, rho, phi, g0, modetype, nstates, _prob);
177 if(_prob>maxprob) {
178 maxprob=_prob;
179 printf("ir = %d,maxprob = %.10f,prob = %.10f\n\n",ir,maxprob,_prob);
180 }
181 }
182 printf("maxprob = %.10f\n", maxprob);
183*/
184 p->initializePhaseSpace(getNDaug(),getDaugs());
185 EvtVector4R ks = p->getDaug(0)->getP4();
186 EvtVector4R kp = p->getDaug(1)->getP4();
187 EvtVector4R pip = p->getDaug(2)->getP4();
188 EvtVector4R pim = p->getDaug(3)->getP4();
189
190 double Pip[4],Pim[4],Kp[4],Ks[4];
191 Ks[0] = ks.get(0); Kp[0] = kp.get(0); Pip[0] = pip.get(0); Pim[0] = pim.get(0);
192 Ks[1] = ks.get(1); Kp[1] = kp.get(1); Pip[1] = pip.get(1); Pim[1] = pim.get(1);
193 Ks[2] = ks.get(2); Kp[2] = kp.get(2); Pip[2] = pip.get(2); Pim[2] = pim.get(2);
194 Ks[3] = ks.get(3); Kp[3] = kp.get(3); Pip[3] = pip.get(3); Pim[3] = pim.get(3);
195
196 //Ks[0] = 0.6096125554; Ks[1] = 0.2036145752; Ks[2] = -0.2256821642; Ks[3] = 0.1884652565;
197 //Kp[0] = 0.7106732604; Kp[1] = -0.2599714599; Kp[2] = 0.3643877399; Kp[3] = 0.2469330230;
198 //Pip[0] = 0.4177584141; Pip[1] = -0.0562004660; Pip[2] = 0.2821155804; Pip[3] =-0.0020089709;
199 //Pim[0] = 0.3945120988; Pim[1] = -0.1225632042; Pim[2] = -0.3645928242; Pim[3] = 0.0521817527;
200
201 //Ks[0] = 0.6439438863; Ks[1] = 0.3368405824; Ks[2] = 0.2308441601; Ks[3] = 0.0171298406;
202 //Kp[0] = 0.7146391007; Kp[1] = 0.4061832461; Kp[2] = -0.1841669581; Kp[3] = 0.2609401579;
203 //Pip[0] = 0.2817222681; Pip[1] = -0.1005161765;Pip[2] = -0.1172908044;Pip[3] = -0.1898077097;
204 //Pim[0] = 0.3638598699;Pim[1] = -0.3106229764;Pim[2] = 0.0905020525;Pim[3] = 0.0907574503;
205
206 double value;
207 int g0[8]={1,1,1,1,1,1,1,1};
208 int nstates=8;
209
210 calEvaMy(Ks,Kp,Pip,Pim,mass,width,rho,phi,g0,modetype,nstates,value);
211 setProb(value);
212
213 return;
214}
215
216void EvtDsToKSKpPipPim::Com_Multi(double a1[2], double a2[2], double res[2]){
217 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
218 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
219}
220void EvtDsToKSKpPipPim::Com_Divide(double a1[2], double a2[2], double res[2]){
221 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/(a2[0]*a2[0]+a2[1]*a2[1]);
222 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/(a2[0]*a2[0]+a2[1]*a2[1]);
223}
224double EvtDsToKSKpPipPim::SCADot(double a1[4], double a2[4])
225{
226 double _cal = 0.;
227 _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
228 return _cal;
229}
230double EvtDsToKSKpPipPim::barrier(int l, double sa, double sb, double sc, double r2)
231{
232 double F;
233 double tmp = sa+sb-sc;
234 double q = 0.25*tmp*tmp/sa-sb;
235 if (q < 0) q = 1e-16;
236 double z = q*r2;
237 if (l == 1) {
238 F = sqrt(2.0*z/(1.0+z));
239 }
240 else if (l == 2) {
241 double z2 = z*z;
242 F = sqrt(13.0*z2/(9.0+3.0*z+z2));
243 } else {
244 F = 1.0;
245 }
246 return F;
247}
248void EvtDsToKSKpPipPim::calt1(double daug1[4], double daug0[4], double t1[4]){
249 double p, pq;
250 double pa[4], qa[4];
251 for(int i=0; i<4; i++){
252 pa[i] = daug1[i] + daug0[i];
253 qa[i] = daug1[i] - daug0[i];
254 }
255 p = SCADot(pa,pa);
256 pq = SCADot(pa,qa);
257 for(int i=0; i<4; i++){
258 t1[i] = qa[i] - pq/p*pa[i];
259 }
260}
261void EvtDsToKSKpPipPim::calt2(double daug1[4], double daug0[4], double t2[4][4]){
262 double p, r;
263 double pa[4], t1[4];
264 calt1(daug1,daug0,t1);
265 r = SCADot(t1,t1)/3.0;
266 for(int i=0; i<4; i++) {
267 pa[i] = daug1[i] + daug0[i];
268 }
269 p = SCADot(pa,pa);
270 for(int i=0; i<4; i++) {
271 for(int j=0; j<4; j++) {
272 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
273 }
274 }
275}
276
277double EvtDsToKSKpPipPim::wid(double mass, double sa, double sb, double sc, double r2, int l){
278 double widm = 0.;
279 double sa0 = mass*mass;
280 double m = sqrt(sa);
281 double tmp = sb-sc;
282 double tmp1 = sa+tmp;
283 double q = 0.25*tmp1*tmp1/sa-sb;
284 if(q<0) q = 1e-16;
285 double tmp2 = sa0+tmp;
286 double q0 = 0.25*tmp2*tmp2/sa0-sb;
287 if(q0<0) q0 = 1e-16;
288 double z = q*r2;
289 double z0 = q0*r2;
290 double t = q/q0;
291 if(l == 0) {widm = sqrt(t)*mass/m;}
292 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
293 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
294 return widm;
295}
296double EvtDsToKSKpPipPim::widl1(double mass, double sa, double sb, double sc, double r2)
297{
298 double widm = 0.;
299 double sa0 = mass*mass;
300 double m = sqrt(sa);
301 double tmp = sb-sc;
302 double tmp1 = sa+tmp;
303 double q = 0.25*tmp1*tmp1/sa-sb;
304 if(q<0) q = 1e-16;
305 double tmp2 = sa0+tmp;
306 double q0 = 0.25*tmp2*tmp2/sa0-sb;
307 if(q0<0) q0 = 1e-16;
308 double z = q*r2;
309 double z0 = q0*r2;
310 double F = (1+z0)/(1+z);
311 double t = q/q0;
312 widm = t*sqrt(t)*mass/m*F;
313 return widm;
314}
315void EvtDsToKSKpPipPim::Flatte_rhoab(double sa, double sb, double sc, double rho[2])
316{
317 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
318 if(q>0) {
319 rho[0]=2* sqrt(q/sa);
320 rho[1]=0;
321 }
322 else if(q<0){
323 rho[0]=0;
324 rho[1]=2*sqrt(-q/sa);
325 }
326}
327void EvtDsToKSKpPipPim::propagatora0980(double mass, double sx, double *sb, double *sc, double prop[2])
328{
329 double unit[2]={1.0};
330 double ci[2]={0,1};
331 double rho1[2];
332 Flatte_rhoab(sx,sb[0],sc[0],rho1);
333 double rho2[2];
334 Flatte_rhoab(sx,sb[1],sc[1],rho2);
335
336 double gKK_a980=0.892*0.341, gPiEta_a980=0.341;
337 double tmp1[2]={gKK_a980,0};
338 double tmp11[2];
339 double tmp2[2]={gPiEta_a980,0};
340 double tmp22[2];
341 Com_Multi(tmp1,rho1,tmp11);
342 Com_Multi(tmp2,rho2,tmp22);
343 double tmp3[2]={tmp11[0]+tmp22[0],tmp11[1]+tmp22[1]};
344 double tmp31[2];
345 Com_Multi(tmp3, ci,tmp31);
346 double tmp4[2]={mass*mass-sx-tmp31[0], -1.0*tmp31[1]};
347 Com_Divide( unit,tmp4, prop);
348}
349void EvtDsToKSKpPipPim::propagatorRBW(double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2]){
350 double a[2], b[2];
351 a[0] = 1;
352 a[1] = 0;
353 b[0] = mass*mass-sa;
354 b[1] = -mass*width*wid(mass,sa,sb,sc,r2,l);
355 Com_Divide(a,b,prop);
356}
357
358
359void EvtDsToKSKpPipPim::KPiSLASS(double sa, double sb, double sc, double prop[2]) {
360 const double m1430 = 1.441;
361 const double sa0 = 2.076481; // m1410*m1410;
362 const double w1430 = 0.193;
363 const double Lass1 = 0.25/sa0;
364 double tmp = sb-sc;
365 double tmp1 = sa0+tmp;
366 double q0 = Lass1*tmp1*tmp1-sb;
367 if(q0<0) q0 = 1e-16;
368 double tmp2 = sa+tmp;
369 double qs = 0.25*tmp2*tmp2/sa-sb;
370 double q = sqrt(qs);
371 double width = w1430*q*m1430/sqrt(sa*q0);
372 double temp_R = atan(m1430*width/(sa0-sa));
373 if(temp_R<0) temp_R += math_pi;
374 double deltaR = -109.7 + temp_R; //fiR=-109.7
375 double temp_F = atan(0.226*q/(2.0-3.8194*qs)); // 2.0*0.113=0.226; 0.113*33.8=1.926; a=0.113; r=-33.8
376 if(temp_F<0) temp_F += math_pi;
377 double deltaF = 0.1 + temp_F; //fiF=0.1
378 double deltaS = deltaR + 2.0*deltaF; ///
379 double t1 = 0.96*sin(deltaF); //F= 0.96
380 double t2 = sin(deltaR);
381 double CF[2], CS[2];
382 CF[0] = cos(deltaF);
383 CF[1] = sin(deltaF);
384 CS[0] = cos(deltaS);
385 CS[1] = sin(deltaS);
386 prop[0] = t1*CF[0] + t2*CS[0];
387 prop[1] = t1*CF[1] + t2*CS[1];
388}
389
390void EvtDsToKSKpPipPim:: propagatorGS(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
391{
392 double a[2], b[2];
393 double tmp = sb-sc;
394 double tmp1 = sa+tmp;
395 double q2 = 0.25*tmp1*tmp1/sa-sb;
396 if(q2<0) q2 = 1e-16;
397
398 double tmp2 = mass2+tmp;
399 double q02 = 0.25*tmp2*tmp2/mass2-sb;
400 if(q02<0) q02 = 1e-16;
401
402 double q = sqrt(q2);
403 double q0 = sqrt(q02);
404 double m = sqrt(sa);
405 double q03 = q0*q02;
406 double tmp3 = log(mass+2*q0)+1.2760418309;
407
408 double h = GS1*q/m*(log(m+2*q)+1.2760418309);
409 double h0 = GS1*q0/mass*tmp3;
410 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
411 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
412 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
413
414 a[0] = 1.0+d*width/mass;
415 a[1] = 0.0;
416 b[0] = mass2-sa+width*f;
417 b[1] = -mass*width*widl1(mass,sa,sb,sc,r2);
418 Com_Divide(a,b,prop);
419}
420
421void EvtDsToKSKpPipPim::calEvaMy(double* Ks,double* Kp, double* Pip, double* Pim, double *mass1, double *width1, double *amp, double *phase,int* g0, int* modetype, int nstates, double & Result)
422{
423 double pD[4],pKpPim[4],pKsPim[4],pKsKpPim[4],pKsKp[4],pKsPipPim[4],pKpPipPim[4],pPipPim[4];
424 for(int i=0;i!=4;i++)
425 {
426
427 pD[i] = Ks[i]+Pip[i]+Kp[i]+Pim[i];
428
429 pKsPim[i] =Ks[i]+Pim[i];
430 pKpPim[i] =Kp[i]+Pim[i];
431 pPipPim[i]=Pip[i]+Pim[i];
432 pKsKp[i] = Kp[i]+Ks[i];
433 pKsKpPim[i] = pKsPim[i] + Kp[i];
434 pKpPipPim[i] = Kp[i]+Pip[i]+Pim[i];
435 pKsPipPim[i] = Ks[i]+Pip[i]+Pim[i];
436
437 }
438
439 double sPim,sPip,sKs,sKp,sKpPim,sD,sKsPim,sKsKpPim,sKsKp,sKpPipPim,sKsPipPim,sPipPim;
440
441 sD = SCADot(pD,pD);
442 sKs = SCADot(Ks,Ks);
443 sKp = SCADot(Kp,Kp);
444 sPip= SCADot(Pip,Pip);
445 sPim= SCADot(Pim,Pim);
446 sPipPim= SCADot(pPipPim,pPipPim);
447 sKsPim = SCADot(pKsPim,pKsPim);
448 sKpPim = SCADot(pKpPim,pKpPim);
449 sKsKpPim = SCADot(pKsKpPim,pKsKpPim);
450 sKpPipPim = SCADot(pKpPipPim,pKpPipPim);
451 sKsPipPim = SCADot(pKsPipPim,pKsPipPim);
452 sKsKp =SCADot(pKsKp,pKsKp);
453
454 double t1D[4],t1KpPim[4],t1KsPim[4],t1KsKpPim_kspim[4],t1KsKpPip_ksk[4],t1KsKp_PipPim[4],t1PipPim[4],t1KsKpPim_ksk[4],t1KsKpPim_kpim[4],t1KsKp[4];
455 double t2KsKpPim_kspim[4][4],t2KsKpPim_kpim[4][4],t2KpPipPim_kpim[4][4],t2KsPipPim_kspim[4][4],t2KpPipPim_pippim[4][4],t2KsPipPim_pippim[4][4];
456
457
458 calt1(Ks,Pim, t1KsPim);
459 calt1(Kp,Pim, t1KpPim);
460 calt1(Ks,Kp, t1KsKp);
461 calt1(Pip,Pim,t1PipPim);
462 calt1(pKsPim,Kp,t1KsKpPim_kspim);
463 calt1(pKpPim,Ks,t1KsKpPim_kpim);
464 calt1(pKsKp,Pip,t1KsKpPip_ksk);
465 calt1(pKsKp,Pim,t1KsKpPim_ksk);
466 calt1(pKsKp,pPipPim,t1KsKp_PipPim);
467 calt1(pKsKpPim,Pip,t1D);
468
469 calt2(pKsPim,Kp,t2KsKpPim_kspim);
470 calt2(pKpPim,Ks,t2KsKpPim_kpim);
471 calt2(pKpPim,Pip,t2KpPipPim_kpim);
472 calt2(pKsPim,Pip,t2KsPipPim_kspim);
473 calt2(pPipPim,Ks,t2KsPipPim_pippim);
474 calt2(pPipPim,Kp,t2KpPipPim_pippim);
475
476 double cof[2],amp_tmp[2],amp_PDF[2], PDF[2];
477 amp_PDF[0] = 0;
478 amp_PDF[1] = 0;
479 PDF[0] = 0;
480 PDF[1] = 0;
481 double temp_PDF,tmp1;
482 double pro[2],pro1[2],pro_kspim[2],pro_kpim[2],pro_ksk[2],pro_pippim[2],pro_kskpim_kspim[2],pro_kskpim_kpim[2];
483 double ispro_kspim=0,ispro_kpim=0,ispro_pippim=0;
484 double barr_kskpim_kspim=-1.0,barr_kskpim_kpim=-1.0,barr_kspim=-1.0,barr_kpim=-1.0,barr_pippim=-1.0,barr_kskp_pippim=-1,barr_ds_kskpim=-1,barr_kspippim_kspim2=-1,barr_kpippim_kpim2=-1,barr_kspippim_pippim2=-1,barr_kpippim_pippim2=-1,barr_ds_kspippim=-1,barr_ds_kpippim=-1,barr_ksk=-1,barr_kspippim_kspim=-1,barr_kpippim_kpim=-1;
485
486 double sKs2[2]={sKs,mass_Pion*mass_Pion};
487 double sKp2[2]={sKp,mass_Eta*mass_Eta};
488 double mass2sq,tt1,tt2,t2D[4][4];
489
490 for(int i=0; i<nstates; i++)
491 {
492 temp_PDF=0; amp_tmp[0]=0; amp_tmp[1]=0; tmp1=0;
493 cof[0] = amp[i]*cos(phase[i]);
494 cof[1] = amp[i]*sin(phase[i]);
495 mass2sq = mRho*mRho;
496
497 //printf("phase= %i, %.5f, %.5f\n", i, amp[i],phase[i]);
498
499
500 //Ds -> eta(1475)pi+ eta(1475)->K*Kp,K*+>Ks pi+
501 if(modetype[i]==3)
502 {
503 propagatorRBW(meta1475,Geta1475,sKsKpPim,sKsPim,sKp,rRes2,1,pro_kskpim_kspim);
504 propagatorRBW(mKstp,GKstp,sKsPim,sKs,sPim,rRes2,1,pro_kspim);
505 ispro_kspim=1;
506
507 Com_Multi(pro_kskpim_kspim,pro_kspim,pro);
508
509 barr_kspim = barrier(1,sKsPim,sKs,sPim,rRes2);
510 barr_kskpim_kspim = barrier(1,sKsKpPim,sKsPim,sKp,rRes2);
511 for(int a=0; a<4; a++)
512 {
513 temp_PDF += Kp[a]*t1KsPim[a]*G[a][a];
514 }
515 tmp1 = barr_kspim*barr_kskpim_kspim*temp_PDF;
516 }
517
518 //Ds -> eta(1475)pi+ eta(1475)->K*Ks,K*0>Kp pi+
519 if(modetype[i]==4)
520 {
521 propagatorRBW(meta1475,Geta1475,sKsKpPim,sKpPim,sKs,rRes2,1,pro_kskpim_kpim);
522 propagatorRBW(mKst0,GKst0,sKpPim,sKp,sPim,rRes2,1,pro_kpim);
523 ispro_kpim=1;
524 Com_Multi(pro_kskpim_kpim,pro_kpim, pro);
525
526 barr_kpim = barrier(1,sKpPim,sKp,sPim,rRes2);
527 barr_kskpim_kpim = barrier(1,sKsKpPim,sKpPim,sKs,rRes2);
528 for(int a=0; a<4; a++)
529 {
530 temp_PDF += Ks[a]*t1KpPim[a]*G[a][a];
531 }
532 tmp1 = barr_kpim*barr_kskpim_kpim*temp_PDF;
533 }
534
535 //0- eta1475 -> a0ipi+
536 if(modetype[i]==15)
537 {
538 propagatorRBW(meta1475,Geta1475,sKsKpPim,sKsKp,sPim,rRes2,0,pro1);
539 propagatora0980(ma0_980,sKsKp,sKs2,sKp2,pro_ksk);
540 Com_Multi(pro1,pro_ksk,pro);
541 tmp1=1.0;
542 }
543
544 //1+ Ds -> K1270 K+ , K1270 -> Ks rho rho -> pi+ pi-
545 if(modetype[i]==7){
546 propagatorRBW(mK1270,GK1270,sKsPipPim,sPipPim,sKs,rRes2,0,pro1);
547 if(ispro_pippim==0){
548 propagatorGS(mass2sq,mRho,GRho,sPipPim,sPip,sPim,rRes2,pro_pippim);
549 ispro_pippim=1;
550 }
551 Com_Multi(pro1,pro_pippim,pro);
552
553 if(barr_pippim<0.0) barr_pippim = barrier(1,sPipPim,sPip,sPim,rRes2);
554 if(barr_ds_kspippim<0) barr_ds_kspippim = barrier(1,sD,sKsPipPim,sKp,rD2);
555 for(int a=0; a<4; a++)
556 {
557 for(int j=0; j<4; j++)
558 {
559 temp_PDF += t1D[a]*(pKsPipPim[a]*pKsPipPim[j]/sKsPipPim-G[a][j])*t1PipPim[j]*G[a][a]*G[j][j];
560 }
561 }
562 tmp1 = barr_pippim*barr_ds_kspippim*temp_PDF;
563 }
564
565 //1+ Ds -> K1270 Kp , K1270 -> K+ rho rho -> pi+ pi-
566 if(modetype[i]==8){
567 propagatorRBW(mK1270,GK1270,sKpPipPim,sPipPim,sKp,rRes2,0,pro1);
568 if(ispro_pippim==0){
569 propagatorGS(mass2sq,mRho,GRho,sPipPim,sPip,sPim,rRes2,pro_pippim);
570 ispro_pippim=1;
571 }
572 Com_Multi(pro1,pro_pippim,pro);
573
574 if(barr_pippim<0.0) barr_pippim = barrier(1,sPipPim,sPip,sPim,rRes2);
575 if(barr_ds_kpippim<0) barr_ds_kpippim = barrier(1,sD,sKpPipPim,sKs,rD2);
576 for(int a=0; a<4; a++)
577 {
578 for(int j=0; j<4; j++)
579 {
580 temp_PDF += t1D[a]*(pKpPipPim[a]*pKpPipPim[j]/sKpPipPim-G[a][j])*t1PipPim[j]*G[a][a]*G[j][j];
581 }
582 }
583 tmp1 = barr_pippim*barr_ds_kpippim*temp_PDF;
584 }
585
586 //1+ Ds -> K1400 K+ , K1400 -> pi+ K*- ,K*- ->Ks pi-
587 if(modetype[i]==11){
588 propagatorRBW(mK1400,GK1400,sKsPipPim,sKsPim,sPip,rRes2,0,pro1);
589 if(ispro_kspim==0){
590 propagatorRBW(mKstp,GKstp,sKsPim,sKs,sPim,rRes2,1,pro_kspim);
591 ispro_kspim=1;
592 }
593 Com_Multi(pro1,pro_kspim,pro);
594
595 if(barr_kspim<0.0) barr_kspim = barrier(1,sKsPim,sKs,sPim,rRes2);
596 if(barr_ds_kspippim<0.0) barr_ds_kspippim = barrier(1,sD,sKsPipPim,sKp,rD2);
597 for(int a=0; a<4; a++)
598 {
599 for(int j=0; j<4; j++)
600 {
601 temp_PDF += t1D[a]*(pKsPipPim[a]*pKsPipPim[j]/sKsPipPim-G[a][j])*t1KsPim[j]*G[a][a]*G[j][j];
602 }
603 }
604 tmp1 = barr_kspim*barr_ds_kspippim*temp_PDF;
605 }
606
607 //1+ Ds -> K1400 Ks , K1400 -> pi+ K*0 ,K*0 ->K+ pi-
608 if(modetype[i]==12){
609 propagatorRBW(mK1400,GK1400,sKpPipPim,sKpPim,sPip,rRes2,0,pro1);
610 if(ispro_kpim==0){
611 propagatorRBW(mKst0,GKst0,sKpPim,sKp,sPim,rRes2,1,pro_kpim);
612 ispro_kpim=1;
613 }
614 Com_Multi(pro1,pro_kpim,pro);
615
616 if(barr_kpim<0.0) barr_kpim = barrier(1,sKpPim,sKp,sPim,rRes2);
617 if(barr_ds_kpippim<0.0) barr_ds_kpippim = barrier(1,sD,sKpPipPim,sKs,rD2);
618 for(int a=0; a<4; a++)
619 {
620 for(int j=0; j<4; j++)
621 {
622 temp_PDF += t1D[a]*(pKpPipPim[a]*pKpPipPim[j]/sKpPipPim-G[a][j])*t1KpPim[j]*G[a][a]*G[j][j];
623 }
624 }
625 tmp1 = barr_kpim*barr_ds_kpippim*temp_PDF;
626
627 //printf("sD:%.15f,%.15f,%.15f,%.15f\n ",sD,sKpPipPim,sKs,rD2);
628 //printf("aaaaa= %.15f,%.15f,%.15f,%.15f\n", barr_kpim,barr_ds_kpippim,temp_PDF,tmp1);
629 }
630
631 //1+ Ds -> phi1680 phi1680->kspim
632 if(modetype[i]==17){
633
634 propagatorRBW(mP1680,GP1680,sKsKpPim,sKsPim,sKp,rRes2,1,pro1);
635 propagatorRBW(mKstp,GKstp,sKsPim,sKs,sPim,rRes2,1,pro_kpim);
636 Com_Multi(pro1,pro_kpim,pro);
637
638 barr_kspim = barrier(1,sKsPim,sKs,sPim,rRes2);
639 barr_kskpim_kspim= barrier(1,sKsKpPim,sKsPim,sKp,rRes2);
640 barr_ds_kskpim = barrier(1,sD,sKsKpPim,sPip,rD2);
641 for(int a=0; a<4; a++)
642 {
643 for(int j=0; j<4; j++)
644 {
645 for(int k=0; k<4; k++)
646 {
647 for(int l=0; l<4; l++)
648 {
649 temp_PDF += E[a][j][k][l]*pKsKpPim[a]*(pKsPim[j]-Kp[j])*Pip[k]*(Ks[l]-Pim[l])*G[a][a]*G[j][j]*G[k][k]*G[l][l];
650 }
651 }
652 }
653 }
654 tmp1 = barr_kspim*barr_kskpim_kspim*barr_ds_kskpim*temp_PDF;
655 }
656
657
658 if(modetype[i]==3||modetype[i]==4||modetype[i]==5||modetype[i]==7||modetype[i]==8||modetype[i]==11||modetype[i]==12||modetype[i]==17||modetype[i]==15){
659 amp_tmp[0] = tmp1*pro[0];
660 amp_tmp[1] = tmp1*pro[1];
661 }
662/*
663 if(modetype[i]==18){
664 amp_tmp[0] = tmp1*pro[0];
665 amp_tmp[1] = tmp1*pro[1];
666
667 }*/
668
669 Com_Multi(amp_tmp,cof,amp_PDF);
670
671 //printf("pro[0],pro[1]:: %.15f, %.15f\n", pro[0],pro[1]);
672 //printf("tmp1:: %.15f\n", tmp1);
673 //printf("amp_tmp[0],amp_tmp[1]:: %.15f, %.15f\n", amp_tmp[0],amp_tmp[1]);
674
675 PDF[0] += amp_PDF[0];
676 PDF[1] += amp_PDF[1];
677
678 //printf("cof= %.15f,%.15f\n",cof[0],cof[1]);
679 //printf("PDF= %.15f,%.15f\n",PDF[0],PDF[1]);
680 //printf("----\n");
681}
682 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
683 if(value <=0) value = 1e-20;
684 Result = value;
685 //printf("value = %.15f\n",value);
686
687}
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")
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
EvtDsToKSKpPipPim::getName
void getName(std::string &name)
Definition EvtDsToKSKpPipPim.cc:35
EvtDsToKSKpPipPim::decay
void decay(EvtParticle *p)
Definition EvtDsToKSKpPipPim.cc:156
EvtDsToKSKpPipPim::clone
EvtDecayBase * clone()
Definition EvtDsToKSKpPipPim.cc:39
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