BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToKKpipipi.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: EvtDsToKKpipipi.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/EvtDsToKKpipipi.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
34EvtDsToKKpipipi::~EvtDsToKKpipipi() {}
35
36void EvtDsToKKpipipi::getName(std::string& model_name){
37 model_name="DsToKKpipipi";
38}
39
43
44void EvtDsToKKpipipi::init(){
45 // check that there are 0 arguments
46 checkNArg(0);
47 checkNDaug(5);
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 checkSpinDaughter(4,EvtSpinType::SCALAR);
55
56 int mode = 0;
57//--------------------amp---------------
58 phi[1] = 1.4732;
59 phi[2] = -4.2901;
60 rho[1] = 0.23396;
61 rho[2] = 7.4744;
62//--------------------mass---------------------
63 mass1[0] = 1.019461;
64 mass1[1] = 1.019461;
65 mass1[2] = 1.019461;
66 mass2[0] = 0.77526;
67 mass2[1] = 0.77526;
68 mass2[2] = 0.77526;
69 mass3[0] = 1.23;
70 mass3[1] = 1.23;
71 mass3[2] = 1.23;
72
73 width1[0] = 0.004249;
74 width1[1] = 0.004249;
75 width1[2] = 0.004249;
76 width2[0] = 0.1478;
77 width2[1] = 0.1478;
78 width2[2] = 0.1478;
79 width3[0] = 0.42;
80 width3[1] = 0.42;
81 width3[2] = 0.42;
82
83 modetype[0]= 1;
84 modetype[1]= 3;
85 modetype[2]= 13;
86
87 phi[0] = 0;
88 rho[0] = 1;
89
90 //cout << "DsToKKpipipi :" << endl;
91 //for (int i=0; i<3; i++) {
92 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
93 //}
94
95 mass_Pion = 0.13957;
96 mass_Pion_N = 0.134977;
97 mass_Eta = 0.547862;
98 math_pi = 3.1415926;
99 rD2 = 25.0; // 5*5
100 rRes2 = 9.0; // 3*3
101 GS1 = 0.636619783;
102 GS2 = 0.01860182466;
103 GS3 = 0.1591549458; // 1/(2*math_2pi)
104 GS4 = 0.00620060822; // mass_Pion2/math_pi
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 for (int j=0; j<4; j++) {
126 G[i][j] = GG[i][j];
127 for (int k=0; k<4; k++) {
128 for (int l=0; l<4; l++) {
129 E[i][j][k][l] = EE[i][j][k][l];
130 }
131 }
132 }
133 }
134
135}
136
137void EvtDsToKKpipipi::initProbMax() {
138 setProbMax(585000.0);
139}
140
141void EvtDsToKKpipipi::decay( EvtParticle *p){
142/*
143 double maxprob = 0.0;
144 for(int ir=0;ir<=60000000;ir++){
145 p->initializePhaseSpace(getNDaug(),getDaugs());
146 EvtVector4R _km = p->getDaug(0)->getP4();
147 EvtVector4R _kp = p->getDaug(1)->getP4();
148 EvtVector4R _pip1 = p->getDaug(2)->getP4();
149 EvtVector4R _pip2 = p->getDaug(3)->getP4();
150 EvtVector4R _pim = p->getDaug(4)->getP4();
151
152 double _Km[4], _Kp[4], _Pip1[4], _Pip2[4], _Pim[4];
153 _Km[0] = _km.get(0); _Kp[0] = _kp.get(0); _Pip1[0] = _pip1.get(0); _Pip2[0] = _pip2.get(0); _Pim[0] = _pim.get(0);
154 _Km[1] = _km.get(1); _Kp[1] = _kp.get(1); _Pip1[1] = _pip1.get(1); _Pip2[1] = _pip2.get(1); _Pim[1] = _pim.get(1);
155 _Km[2] = _km.get(2); _Kp[2] = _kp.get(2); _Pip1[2] = _pip1.get(2); _Pip2[2] = _pip2.get(2); _Pim[2] = _pim.get(2);
156 _Km[3] = _km.get(3); _Kp[3] = _kp.get(3); _Pip1[3] = _pip1.get(3); _Pip2[3] = _pip2.get(3); _Pim[3] = _pim.get(3);
157
158 double _prob;
159 int nstates=3;
160 int g0[3]={1,1,1};
161 int g1[3]={1,1,1};
162 int g2[3]={0,0,0};
163 int g3[3]={0,1,0};
164 int g4[3]={0,0,0};
165 calEvaMy(_Km,_Kp,_Pip1,_Pip2,_Pim,mass1,mass2,mass3,width1,width2,width3,rho,phi,modetype,nstates,_prob);
166 if(_prob>maxprob) {
167 maxprob=_prob;
168 cout << "Max PDF = " << ir << " " << _prob << endl;
169 }
170 }
171 cout << "Max!!!!!!!!!!! " << maxprob<< endl;
172*/
173 p->initializePhaseSpace(getNDaug(),getDaugs());
174 EvtVector4R km = p->getDaug(0)->getP4();
175 EvtVector4R kp = p->getDaug(1)->getP4();
176 EvtVector4R pip1 = p->getDaug(2)->getP4();
177 EvtVector4R pip2 = p->getDaug(3)->getP4();
178 EvtVector4R pim = p->getDaug(4)->getP4();
179
180 double Km[4],Kp[4], Pip1[4],Pip2[4],Pim[4];
181 Km[0] = km.get(0); Kp[0] = kp.get(0); Pip1[0] = pip1.get(0); Pip2[0] = pip2.get(0); Pim[0] = pim.get(0);
182 Km[1] = km.get(1); Kp[1] = kp.get(1); Pip1[1] = pip1.get(1); Pip2[1] = pip2.get(1); Pim[1] = pim.get(1);
183 Km[2] = km.get(2); Kp[2] = kp.get(2); Pip1[2] = pip1.get(2); Pip2[2] = pip2.get(2); Pim[2] = pim.get(2);
184 Km[3] = km.get(3); Kp[3] = kp.get(3); Pip1[3] = pip1.get(3); Pip2[3] = pip2.get(3); Pim[3] = pim.get(3);
185
186 double prob;
187 int nstates=3;
188 calEvaMy(Km,Kp,Pip1,Pip2,Pim,mass1,mass2,mass3,width1,width2,width3,rho,phi,modetype,nstates,prob);
189
190 setProb(prob);
191
192 return;
193}
194
195void EvtDsToKKpipipi::Com_Multi(double a1[2], double a2[2], double res[2]){
196 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
197 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
198}
199void EvtDsToKKpipipi::Com_Divide(double a1[2], double a2[2], double res[2]){
200 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
201 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
202 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
203}
204double EvtDsToKKpipipi::SCADot(double a1[4], double a2[4]){
205 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
206 return _cal;
207}
208double EvtDsToKKpipipi::barrier(int l, double sa, double sb, double sc, double r, double mass)
209{
210 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
211 if(q < 0) q = 1e-16;
212 double z;
213 z = q*r*r;
214 double sa0;
215 sa0 = mass*mass;
216 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
217 if(q0 < 0) q0 = 1e-16;
218 double z0 = q0*r*r;
219 double F = 0.0;
220 if(l == 0) F = 1;
221 if(l == 1) F = sqrt((1+z0)/(1+z));
222 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
223 return F;
224}
225//double EvtDsToKKpipipi::barrier(int l, double sa, double sb, double sc, double r2){
226// double F;
227// double tmp = sa+sb-sc;
228// double q = 0.25*tmp*tmp/sa-sb;
229// if (q < 0) q = 1e-16;
230// double z = q*r2;
231// if (l == 1) {
232// F = sqrt(2.0*z/(1.0+z));
233// }
234// else if (l == 2) {
235// double z2 = z*z;
236// F = sqrt(13.0*z2/(9.0+3.0*z+z2));
237// } else {
238// F = 1.0;
239// }
240// return F;
241//}
242void EvtDsToKKpipipi::calt1(double daug1[4], double daug2[4], double t1[4]){
243 double p, pq, tmp;
244 double pa[4], qa[4];
245 for(int i=0; i<4; i++) {
246 pa[i] = daug1[i] + daug2[i];
247 qa[i] = daug1[i] - daug2[i];
248 }
249 p = SCADot(pa,pa);
250 pq = SCADot(pa,qa);
251 tmp = pq/p;
252 for(int i=0; i<4; i++) {
253 t1[i] = qa[i] - tmp*pa[i];
254 }
255}
256void EvtDsToKKpipipi::calt2(double daug1[4], double daug2[4], double t2[4][4]){
257 double p, r;
258 double pa[4], t1[4];
259 calt1(daug1,daug2,t1);
260 r = SCADot(t1,t1)/3.0;
261 for(int i=0; i<4; i++) {
262 pa[i] = daug1[i] + daug2[i];
263 }
264 p = SCADot(pa,pa);
265 for(int i=0; i<4; i++) {
266 for(int j=0; j<4; j++) {
267 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
268 }
269 }
270}
271void EvtDsToKKpipipi::Flatte_rhoab(double sa, double sb, double sc, double rho[2]){
272 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
273 if(q>0) {
274 rho[0]=2* sqrt(q/sa);
275 rho[1]=0;
276 }
277 else if(q<0){
278 rho[0]=0;
279 rho[1]=2*sqrt(-q/sa);
280 }
281}
282void EvtDsToKKpipipi::propagatora0980(double mass2, double mass, double sx, double *sb, double *sc, double prop[2]){
283 double unit[2]={1.0};
284 double ci[2]={0,1};
285 double rho1[2];
286 Flatte_rhoab(sx,sb[0],sc[0],rho1);
287 double rho2[2];
288 Flatte_rhoab(sx,sb[1],sc[1],rho2);
289
290 double gKK_a980=0.892*0.341, gPiEta_a980=0.341;
291 double tmp1[2]={gKK_a980,0};
292 double tmp11[2];
293 double tmp2[2]={gPiEta_a980,0};
294 double tmp22[2];
295 Com_Multi(tmp1,rho1,tmp11);
296 Com_Multi(tmp2,rho2,tmp22);
297 double tmp3[2]={tmp11[0]+tmp22[0],tmp11[1]+tmp22[1]};
298 double tmp31[2];
299 Com_Multi(tmp3, ci,tmp31);
300 double tmp4[2]={mass2-sx-tmp31[0], -1.0*tmp31[1]};
301 Com_Divide( unit,tmp4, prop);
302}
303double EvtDsToKKpipipi::wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l){
304 double widm = 0.;
305 double m = sqrt(sa);
306 double tmp = sb-sc;
307 double tmp1 = sa+tmp;
308 double q = 0.25*tmp1*tmp1/sa-sb;
309 if(q<0) q = 1e-16;
310 double tmp2 = mass2+tmp;
311 double q0 = 0.25*tmp2*tmp2/mass2-sb;
312 if(q0<0) q0 = 1e-16;
313 double z = q*r2;
314 double z0 = q0*r2;
315 double t = q/q0;
316 if(l == 0){widm = sqrt(t)*mass/m;}
317 else if(l == 1){widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
318 else if(l == 2){widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
319 return widm;
320}
321double EvtDsToKKpipipi::widl1(double mass2, double mass, double sa, double sb, double sc, double r2){
322 double widm = 0.;
323 double m = sqrt(sa);
324 double tmp = sb-sc;
325 double tmp1 = sa+tmp;
326 double q = 0.25*tmp1*tmp1/sa-sb;
327 if(q<0) q = 1e-16;
328 double tmp2 = mass2+tmp;
329 double q0 = 0.25*tmp2*tmp2/mass2-sb;
330 if(q0<0) q0 = 1e-16;
331 double z = q*r2;
332 double z0 = q0*r2;
333 double F = (1+z0)/(1+z);
334 double t = q/q0;
335 widm = t*sqrt(t)*mass/m*F;
336 return widm;
337}
338void EvtDsToKKpipipi::propagatorRBW(double mass2, double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2]){
339 double a[2], b[2];
340 a[0] = 1;
341 a[1] = 0;
342 b[0] = mass2-sa;
343 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
344 Com_Divide(a,b,prop);
345}
346void EvtDsToKKpipipi::propagatorNBW(double mass2, double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2]){
347 double a[2], b[2];
348 a[0] = 1;
349 a[1] = 0;
350 b[0] = mass2-sa;
351 b[1] = -mass*width;
352 Com_Divide(a,b,prop);
353}
354void EvtDsToKKpipipi::propagatorRBWl1(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2]){
355 double a[2], b[2];
356 a[0] = 1;
357 a[1] = 0;
358 b[0] = mass2-sa;
359 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
360 Com_Divide(a,b,prop);
361}
362void EvtDsToKKpipipi::propagatorGS(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2]){
363 double a[2], b[2];
364 double tmp = sb-sc;
365 double tmp1 = sa+tmp;
366 double q2 = 0.25*tmp1*tmp1/sa-sb;
367 if(q2<0) q2 = 1e-16;
368
369 double tmp2 = mass2+tmp;
370 double q02 = 0.25*tmp2*tmp2/mass2-sb;
371 if(q02<0) q02 = 1e-16;
372
373 double q = sqrt(q2);
374 double q0 = sqrt(q02);
375 double m = sqrt(sa);
376 double q03 = q0*q02;
377 double tmp3 = log(mass+2*q0)+1.2760418309; // log(mass_2Pion) = 1.2760418309;
378
379 double h = GS1*q/m*(log(m+2*q)+1.2760418309);
380 double h0 = GS1*q0/mass*tmp3;
381 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
382 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
383 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
384
385 a[0] = 1.0+d*width/mass;
386 a[1] = 0.0;
387 b[0] = mass2-sa+width*f;
388 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
389 Com_Divide(a,b,prop);
390}
391void EvtDsToKKpipipi::KPiSLASS(double sa, double sb, double sc, double prop[2]){
392 const double m1430 = 1.441;
393 const double sa0 = 2.076481; // m1430*m1430;
394 const double w1430 = 0.193;
395 const double Lass1 = 0.25/sa0;
396 double tmp = sb-sc;
397 double tmp1 = sa0+tmp;
398 double q0 = Lass1*tmp1*tmp1-sb;
399 if(q0<0) q0 = 1e-16;
400 double tmp2 = sa+tmp;
401 double qs = 0.25*tmp2*tmp2/sa-sb;
402 double q = sqrt(qs);
403 double width = w1430*q*m1430/sqrt(sa*q0);
404 double temp_R = atan(m1430*width/(sa0-sa));
405 if(temp_R<0) temp_R += math_pi;
406 double deltaR = -109.7*math_pi/180.0 + temp_R;
407 double temp_F = atan(0.226*q/(2.0-3.8194*qs)); // 2.0*1.07 = 2.14; 1.8*1.07 = 1.926
408 if(temp_F<0) temp_F += math_pi;
409 double deltaF = 0.1*math_pi/180.0 + temp_F;
410 double deltaS = deltaR + 2.0*deltaF;
411 double t1 = 0.96*sin(deltaF);
412 double t2 = sin(deltaR);
413 double CF[2], CS[2];
414 CF[0] = cos(deltaF);
415 CF[1] = sin(deltaF);
416 CS[0] = cos(deltaS);
417 CS[1] = sin(deltaS);
418 prop[0] = t1*CF[0] + t2*CS[0];
419 prop[1] = t1*CF[1] + t2*CS[1];
420}
421
422void EvtDsToKKpipipi::calEvaMy(double* Km, double* Kp, double* Pip1, double* Pip2, double* Pim, double *mass1, double *mass2, double *mass3, double *width1, double *width2, double *width3, double *amp, double *phase, int* modetype, int nstates, double & Result){
423 double cof[2], amp_tmp[2], amp_PDF[2], PDF[2];
424 //Km[0] = sqrt(0.243719942 + Km[1]*Km[1] + Km[2]*Km[2] + Km[3]*Km[3]);
425 //Kp[0] = sqrt(0.243719942 + Kp[1]*Kp[1] + Kp[2]*Kp[2] + Kp[3]*Kp[3]);
426 //Pip1[0] = sqrt(0.019479785 + Pip1[1]*Pip1[1] + Pip1[2]*Pip1[2] + Pip1[3]*Pip1[3]);
427 //Pip2[0] = sqrt(0.019479785 + Pip2[1]*Pip2[1] + Pip2[2]*Pip2[2] + Pip2[3]*Pip2[3]);
428 //Pim[0] = sqrt(0.019479785 + Pim[1]*Pim[1] + Pim[2]*Pim[2] + Pim[3]*Pim[3]);
429 //Km[0] = 0.603525; Kp[0] = 0.643143; Pip1[0] = 0.225747; Pip2[0] = 0.332362; Pim[0] = 0.215164;
430 //Km[1] = 0.155547; Kp[1] = 0.0828093; Pip1[1] = 0.0581283; Pip2[1] = -0.25147; Pim[1] = 0.09105;
431 //Km[2] = -0.275397; Kp[2] = 0.261094; Pip1[2] = 0.161565; Pip2[2] = 0.150701; Pim[2] = -0.0387834;
432 //Km[3] = -0.143121; Kp[3] = -0.308036; Pip1[3] = 0.0447196; Pip2[3] = -0.0709744; Pim[3] = 0.130466;
433//Km[0] = 0.559664; Kp[0] = 0.621872; Pip1[0] = 0.206021; Pip2[0] = 0.387562; Pim[0] = 0.267452;
434//Km[1] = 0.263265; Kp[1] = -0.20385; Pip1[1] = -0.0550443; Pip2[1] = 0.182318; Pim[1] = 0.0824861;
435//Km[2] = 0.009757; Kp[2] = -0.316095; Pip1[2] = -0.061181; Pip2[2] = -0.309785; Pim[2] = 0.211064;
436//Km[3] = 0.0101476; Kp[3] = 0.0392071; Pip1[3] = -0.127247; Pip2[3] = -0.0389575; Pim[3] = 0.0264334;
437//---------------------------------------------------------------
438 double prho1[4],prho2[4],pphi[4],pa11[4],pa12[4],pD1[4],pD2[4],pD3[4],pD4[4],psigma1[4],psigma2[4],pa13[4],pa14[4];
439 for(int i=0;i!=4;i++){
440 pphi[i] =Km[i] +Kp[i];
441 prho1[i] =Pim[i] +Pip1[i];
442 prho2[i] =Pim[i] +Pip2[i];
443 pa11[i] =prho1[i] +Pip2[i];
444 pa12[i] =prho2[i] +Pip1[i];
445 psigma1[i] =Pim[i] +Pip1[i];
446 psigma2[i] =Pim[i] +Pip2[i];
447 pa13[i] =psigma1[i] +Pip2[i];
448 pa14[i] =psigma2[i] +Pip1[i];
449 pD1[i] =pa11[i] +pphi[i];
450 pD2[i] =pa12[i] +pphi[i];
451 pD3[i] =pa13[i] +pphi[i];
452 pD4[i] =pa14[i] +pphi[i];
453 }
454
455 double skaon1,skaon2,spion1,spion2,spim,sphi,sa11,sa12,srho1,srho2,sD1,sD2,sD3,sD4,ssigma1,ssigma2,sa13,sa14;
456 skaon1 = SCADot(Km, Km);
457 skaon2 = SCADot(Kp, Kp);
458 sphi = SCADot(pphi, pphi);
459
460 spion1 = SCADot(Pip1, Pip1);
461 spion2 = SCADot(Pip2, Pip2);
462 spim = SCADot(Pim, Pim);
463
464 srho1 = SCADot(prho1, prho1);
465 srho2 = SCADot(prho2, prho2);
466 sa11 = SCADot(pa11, pa11);
467 sa12 = SCADot(pa12, pa12);
468
469 sD1 = SCADot(pD1, pD1);
470 sD2 = SCADot(pD2, pD2);
471 sD3 = SCADot(pD3, pD3);
472 sD4 = SCADot(pD4, pD4);
473
474 ssigma1 = SCADot(psigma1,psigma1);
475 ssigma2 = SCADot(psigma2,psigma2);
476 sa13 = SCADot(pa13, pa13);
477 sa14 = SCADot(pa14, pa14);
478//-----------------------------------------------------------------------------------------------
479 double t1phi[4], t1rho1[4], t2a11[4][4], t1sigma1[4], t1a13[4], t1D1[4], t1D3[4], t2D1[4][4], t2D3[4][4],
480 t1rho2[4], t2a12[4][4], t1sigma2[4], t1a14[4], t1D2[4], t1D4[4], t2D2[4][4], t2D4[4][4];
481
482 calt1(Km, Kp, t1phi);
483
484 calt1(Pip1, Pim, t1rho1);
485 calt1(Pip2, Pim, t1rho2);
486 calt1(Pip1, Pim, t1sigma1);
487 calt1(Pip2, Pim, t1sigma2);
488 calt1(psigma1,Pip2, t1a13);
489 calt1(psigma2,Pip1, t1a14);
490
491 calt1(pa11,pphi,t1D1);
492 calt1(pa12,pphi,t1D2);
493 calt1(pa13,pphi,t1D3);
494 calt1(pa14,pphi,t1D4);
495
496 calt2(prho1,Pip2, t2a11);
497 calt2(prho2,Pip1, t2a12);
498
499 calt2(pa11,pphi,t2D1);
500 calt2(pa12,pphi,t2D2);
501 calt2(pa13,pphi,t2D3);
502 calt2(pa14,pphi,t2D4);
503
504 amp_PDF[0] = 0;
505 amp_PDF[1] = 0;
506 PDF[0] = 0;
507 PDF[1] = 0;
508 double tt1,tt2,tt4;
509 double temp_PDF, tmp1, tmp2, amp_tmp1[2], amp_tmp2[2];
510 double pro[2], pro0[2], pro1[2],pro2[2],pro3[2];
511 double mass1sq, mass2sq, mass3sq;
512 for(int i=0; i<nstates; i++){
513 tmp1 = 0;tmp2 = 0;temp_PDF = 0;
514 cof[0] = amp[i]*cos(phase[i]);
515 cof[1] = amp[i]*sin(phase[i]);
516 mass1sq = mass1[i]*mass1[i];
517 mass2sq = mass2[i]*mass2[i];
518 mass3sq = mass3[i]*mass3[i];
519
520 amp_tmp[0] = 0;
521 amp_tmp[1] = 0;
522 amp_tmp1[0] = 0;
523 amp_tmp1[1] = 0;
524 amp_tmp2[0] = 0;
525 amp_tmp2[1] = 0;
526 pro[0] = 0;
527 pro[1] = 0;
528 pro1[0] = 0;
529 pro1[1] = 0;
530 pro2[0] = 0;
531 pro2[1] = 0;
532 pro3[0] = 0;
533 pro3[1] = 0;
534 if(modetype[i]==1){
535 double B_phi=-1.0, B_rho1=-1.0;
536 propagatorRBWl1(mass1sq, mass1[i], width1[i], sphi, skaon1, skaon2, rRes2, pro0);
537 propagatorGS(mass2sq, mass2[i], width2[i], srho1, spion1, spim, rRes2, pro1);
538 propagatorRBW(mass3sq, mass3[i], width3[i], sa11, srho1, spion2, rRes2, 0 , pro2);
539 Com_Multi(pro0,pro1,pro3);
540 Com_Multi(pro2,pro3,pro);
541 if (B_phi<0.0) B_phi = barrier(1, sphi, skaon1, skaon2, rRes2, 1.019461);
542 if (B_rho1<0.0) B_rho1 = barrier(1, srho1, spion1, spim, rRes2, 0.77526);
543 for(int a=0; a<4; a++){
544 for(int j=0; j<4; j++){
545 temp_PDF += (G[a][j]-pa11[a]*pa11[j]/sa11)*t1phi[a]*t1rho1[j]*G[a][a]*G[j][j];
546 }
547 }
548 tmp1 = B_phi*B_rho1*temp_PDF;
549 amp_tmp1[0] = tmp1*pro[0];
550 amp_tmp1[1] = tmp1*pro[1];
551
552 temp_PDF=0;
553 double B_rho2=-1.0;
554 propagatorRBWl1(mass1sq, mass1[i], width1[i], sphi, skaon1, skaon2, rRes2, pro0);
555 propagatorGS(mass2sq, mass2[i], width2[i], srho2, spion2, spim, rRes2, pro1);
556 propagatorRBW(mass3sq, mass3[i], width3[i], sa12, srho2, spion1, rRes2, 0, pro2);
557 Com_Multi(pro0,pro1,pro3);
558 Com_Multi(pro2,pro3,pro);
559 if (B_phi<0.0) B_phi = barrier(1, sphi, skaon1, skaon2, rRes2, 1.019461);
560 if (B_rho2<0.0) B_rho2 = barrier(1, srho2, spion2, spim, rRes2, 0.77526);
561 for(int a=0; a<4; a++){
562 for(int j=0; j<4; j++){
563 temp_PDF += (G[a][j]-pa12[a]*pa12[j]/sa12)*t1phi[a]*t1rho2[j]*G[a][a]*G[j][j];
564 }
565 }
566 tmp2 = B_phi*B_rho2*temp_PDF;
567 amp_tmp2[0] = tmp2*pro[0];
568 amp_tmp2[1] = tmp2*pro[1];
569 }
570 if(modetype[i]==3){
571 double B_phi=-1.0, B_rho1=-1.0, B_phia11=-1.0;
572 propagatorRBWl1(mass1sq, mass1[i], width1[i], sphi, skaon1, skaon2, rRes2, pro0);
573 propagatorGS(mass2sq, mass2[i], width2[i], srho1, spion1, spim, rRes2, pro1);
574 propagatorRBW(mass3sq, mass3[i], width3[i], sa11, srho1, spion2, rRes2, 0 , pro2);
575//printf("mass1[0]:%.6f\n",mass1[0]);
576//printf("mass1[1]:%.6f\n",mass1[1]);
577//printf("width1[0]:%.6f\n",width1[0]);
578//printf("width1[1]:%.6f\n",width1[1]);
579//printf("pro0[0]:%.6f\n",pro0[0]);
580//printf("pro0[1]:%.6f\n",pro0[1]);
581//printf("pro1[0]:%.6f\n",pro1[0]);
582//printf("pro1[1]:%.6f\n",pro1[1]);
583//printf("pro2[0]:%.6f\n",pro2[0]);
584//printf("pro2[1]:%.6f\n",pro2[1]);
585 Com_Multi(pro0,pro1,pro3);
586 Com_Multi(pro2,pro3,pro);
587 if (B_phi<0.0) B_phi = barrier(1, sphi, skaon1, skaon2, rRes2, 1.019461);
588 if (B_rho1<0.0) B_rho1 = barrier(1, srho1, spion1, spim, rRes2, 0.77526);
589 if (B_phia11<0.0) B_phia11 = barrier(1, sD1, sphi, sa11, rD2, 1.9683);
590 for(int w=0; w<4; w++){
591 tt1=pD1[w]*G[w][w];
592 for(int j=0; j<4; j++){
593 tt2=t1D1[j]*G[j][j];
594 for(int m=0; m<4; m++){
595 for(int k=0; k<4; k++){
596 tt4=t1phi[k]*G[k][k];
597 for(int l=0; l<4; l++){
598 temp_PDF += E[w][m][k][l]*G[m][m]*tt1*tt2*(G[m][j]-pa11[m]*pa11[j]/sa11)*tt4*t1rho1[l]*G[l][l];
599 }
600 }
601 }
602 }
603 }
604 tmp1 = B_phi*B_rho1*B_phia11*temp_PDF;//a1(l=0),B_a11=1
605 amp_tmp1[0] = tmp1*pro[0];
606 amp_tmp1[1] = tmp1*pro[1];
607
608 temp_PDF=0;
609 double B_rho2=-1.0, B_phia12=-1.0;
610 propagatorRBWl1(mass1sq, mass1[i], width1[i], sphi, skaon1, skaon2, rRes2, pro0);
611 propagatorGS(mass2sq, mass2[i], width2[i], srho2, spion2, spim, rRes2, pro1);
612 propagatorRBW(mass3sq, mass3[i], width3[i], sa12, srho2, spion1, rRes2, 0, pro2);
613 Com_Multi(pro0,pro1,pro3);
614 Com_Multi(pro2,pro3,pro);
615 if (B_phi<0.0) B_phi = barrier(1, sphi, skaon1, skaon2, rRes2, 1.019461);
616 if (B_rho2<0.0) B_rho2 = barrier(1, srho2, spion2, spim, rRes2, 0.77526);
617 if (B_phia12<0.0) B_phia12 = barrier(1, sD2, sphi, sa12, rD2, 1.9683);
618 for(int w=0; w<4; w++){
619 tt1=pD2[w]*G[w][w];
620 for(int j=0; j<4; j++){
621 tt2=t1D2[j]*G[j][j];
622 for(int m=0; m<4; m++){
623 for(int k=0; k<4; k++){
624 tt4=t1phi[k]*G[k][k];
625 for(int l=0; l<4; l++){
626 temp_PDF += E[w][m][k][l]*G[m][m]*tt1*tt2*(G[m][j]-pa12[m]*pa12[j]/sa12)*tt4*t1rho2[l]*G[l][l];
627 }
628 }
629 }
630 }
631 }
632 tmp2 = B_phi*B_rho2*B_phia12*temp_PDF;//a1(l=0),B_a12=1
633 amp_tmp2[0] = tmp2*pro[0];
634 amp_tmp2[1] = tmp2*pro[1];
635 }
636 if(modetype[i]==13){
637 amp_tmp1[0] = 1.;
638 amp_tmp1[1] = 0.;
639 amp_tmp2[0] = 1.;
640 amp_tmp2[1] = 0.;
641 }
642 amp_tmp[0] = amp_tmp1[0]+amp_tmp2[0];
643 amp_tmp[1] = amp_tmp1[1]+amp_tmp2[1];
644 Com_Multi(amp_tmp,cof,amp_PDF);
645 PDF[0] += amp_PDF[0];
646 PDF[1] += amp_PDF[1];
647 }
648 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
649 //printf("pdf:%.6f\n",value);
650 if(value <=0) value = 1e-20;
651 Result = value;
652// cout<<"!!!!!!!!!!!!"<<setprecision(20)<<value<<endl;
653}
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")
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
EvtDsToKKpipipi::~EvtDsToKKpipipi
virtual ~EvtDsToKKpipipi()
Definition EvtDsToKKpipipi.cc:34
EvtDsToKKpipipi::getName
void getName(std::string &name)
Definition EvtDsToKKpipipi.cc:36
EvtDsToKKpipipi::decay
void decay(EvtParticle *p)
Definition EvtDsToKKpipipi.cc:141
EvtDsToKKpipipi::clone
EvtDecayBase * clone()
Definition EvtDsToKKpipipi.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