BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToKKpi.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: EvtDsToKKpi.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/EvtDsToKKpi.hh"
27#include <string>
28#include "EvtGenBase/EvtConst.hh"
29#include "EvtGenBase/EvtFlatte.hh"
30#include "EvtGenBase/EvtDecayTable.hh"
31using std::endl;
32
33EvtDsToKKpi::~EvtDsToKKpi() {}
34
35void EvtDsToKKpi::getName(std::string& model_name){
36 model_name="DsToKKpi";
37}
38
42
43void EvtDsToKKpi::init(){
44 // check that there are 0 arguments
45 checkNArg(0);
46 checkNDaug(3);
47 checkSpinParent(EvtSpinType::SCALAR);
48 checkSpinDaughter(0,EvtSpinType::SCALAR);
49 checkSpinDaughter(1,EvtSpinType::SCALAR);
50 checkSpinDaughter(2,EvtSpinType::SCALAR);
51
52 phi[1] = 6.1944e+00;
53 phi[2] = 4.7398e+00;
54 phi[3] = 2.9047e+00;
55 phi[4] = 1.0068e+00;
56 phi[5] = 5.8035e-01;
57
58 rho[1] = 1.0963e+00;
59 rho[2] = 2.7818e+00;
60 rho[3] = 1.2570e+00;
61 rho[4] = 7.7351e-01;
62 rho[5] = 5.6277e-01;
63
64 mass[0] = 8.9581e-01;
65 mass[1] = 1.019461e+00;
66 mass[2] = 0.919;
67 mass[3] = 1.4712e+00;
68 mass[4] = 1.7220e+00;
69 mass[5] = 1.3500e+00;
70
71 width[0] = 4.7400e-02;
72 width[1] = 4.2660e-03;
73 width[2] = 0.272;
74 width[3] = 2.7000e-01;
75 width[4] = 1.3500e-01;
76 width[5] = 2.6500e-01;
77
78 modetype[0]= 0;
79 modetype[1]= 1;
80 modetype[2]= 13;
81 modetype[3]= 3;
82 modetype[4]= 4;
83 modetype[5]= 5;
84
85
86 phi[0] = 0;
87 rho[0] = 1;
88
89 //cout << "DsToKKpi :"<< endl;
90 //for (int i=0; i<6; i++) {
91 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
92 //}
93
94 mD = 1.86486;
95 mDs = 1.9683;
96 rRes = 1.5;
97 rD = 3.0;
98 mkstr = 0.89594;
99 mk0 = 0.497614;
100 mass_Kaon = 0.49368;
101 mass_Pion = 0.13957;
102 mass_Pi0 = 0.1349766;
103 mass_EtaP = 0.95778;
104 mass_Eta = 0.547862;
105 math_pi = 3.1415926;
106 afRatio= 2.04835; //BABAR // afRatio= 1.23202; //BES2
107 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
108 for (int i=0; i<4; i++) {
109 for (int j=0; j<4; j++) {
110 G[i][j] = GG[i][j];
111 }
112 }
113}
114
115void EvtDsToKKpi::initProbMax() {
116 setProbMax(35638.0);
117}
118
119void EvtDsToKKpi::decay( EvtParticle *p){
120/*
121 double maxprob = 0.0;
122 cout<<" calculate maxprob for DsToKKpi: "<<endl;
123 for(int ir=0;ir<=60000000;ir++){
124 p->initializePhaseSpace(getNDaug(),getDaugs());
125 EvtVector4R D1 = p->getDaug(0)->getP4();
126 EvtVector4R D2 = p->getDaug(1)->getP4();
127 EvtVector4R D3 = p->getDaug(2)->getP4();
128 double P1[4], P2[4], P3[4];
129 P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
130 P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
131 P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
132 double value;
133 int g0[6]={1,1,1,1,1,1};
134 int nstates=6;
135 calEvaMy( P1, P2, P3, mass, width, rho, phi, g0, modetype, nstates, value);
136 if(value>maxprob) {
137 maxprob=value;
138 cout << "Max PDF = " << ir << " " << maxprob << endl;
139 }
140 }
141 printf("MAXprob = %.10f\n",maxprob);
142*/
143 p->initializePhaseSpace(getNDaug(),getDaugs());
144 EvtVector4R D1 = p->getDaug(0)->getP4();
145 EvtVector4R D2 = p->getDaug(1)->getP4();
146 EvtVector4R D3 = p->getDaug(2)->getP4();
147 double P1[4], P2[4], P3[4];
148 P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
149 P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
150 P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
151
152 double value;
153 int g0[6]={1,1,1,1,1,1};
154 int nstates=6;
155 calEvaMy(P1, P2, P3, mass, width, rho, phi, g0, modetype, nstates, value);
156 setProb(value);
157 return ;
158}
159
160void EvtDsToKKpi::MIP_LineShape(double sa, double pro[2]){
161 double m0 = mass[2], cKK = width[2];
162 pro[0] = sqrt(1/(((m0 * m0)-sa)*((m0 * m0)-sa) + (cKK*m0* sqrt(1 - 4 * mass_Kaon*mass_Kaon/sa) )*(cKK*m0* sqrt(1 - 4 * mass_Kaon*mass_Kaon/sa) ) ));
163 pro[1] = 0;
164}
165void EvtDsToKKpi::calEvaMy(double* pKm, double* pKp, double* pPi, double *mass1, double *width1, double *amp, double *phase,int* g0, int* modetype, int nstates, double & Result)
166{
167 double pF0_980[4],pPhi1020[4], pKstr[4],pD[4],p23[4];
168 pKp[0] = sqrt( mass_Kaon*mass_Kaon + pKp[1]*pKp[1] + pKp[2]*pKp[2] + pKp[3]*pKp[3]);
169 pKm[0] = sqrt( mass_Kaon*mass_Kaon + pKm[1]*pKm[1] + pKm[2]*pKm[2] + pKm[3]*pKm[3]);
170 pPi[0] = sqrt( mass_Pion*mass_Pion + pPi[1]*pPi[1] + pPi[2]*pPi[2] + pPi[3]*pPi[3]);
171 for(int u=0; u<4;u++){
172 pF0_980[u]=pKp[u]+pKm[u];
173 pPhi1020[u]=pKp[u]+pKm[u];
174 pKstr[u]=pKm[u]+pPi[u];
175 p23[u]=pKp[u]+pPi[u];
176 pD[u]=pKp[u]+pKm[u]+pPi[u];
177 }
178 double cof[2], amp_tmp1[2], amp_tmp2[2], amp_tmp[2], amp_PDF[2], PDF[2];
179 double rrho = 3.0;
180 double ra1 = 3.0;
181 //----------------------------------------------------
182 amp_PDF[0] = 0;
183 amp_PDF[1] = 0;
184 PDF[0] = 0;
185 PDF[1] = 0;
186 double temp_PDF, tmp1, tmp2, temp[2];
187 double pro[2], pro0[2], pro1[2];
188 double t1D[4],t2D[4][4],B[3], t1Tm[4];
189 double t1kstr1[4],t1phi1020[4], t1D1[4],t1D2[4];
190
191 double sD, sF0_980, sF0_1710,sF0_1370, sKp, sKm, sPi, sKstr, sPhi1020, sKstr1430;
192 sF0_980=SCADot(pF0_980,pF0_980);
193 sF0_1710=sF0_980;
194 sF0_1370=sF0_980;
195 sKstr=SCADot(pKstr,pKstr);
196 sD=SCADot(pD,pD);
197 sKstr1430=sKstr;
198 sPhi1020=SCADot(pPhi1020,pPhi1020);
199 sKp=SCADot(pKp,pKp);
200 sKm=SCADot(pKm,pKm);
201 sPi=SCADot(pPi,pPi);
202
203 calt1(pKm,pPi,t1kstr1);
204 calt1(pKm,pKp,t1phi1020);
205 calt1(pKstr,pKp,t1D1);
206 calt1(pPhi1020,pPi,t1D2);
207
208 double mDs=sqrt(sD);
209 for(int i=0; i<nstates; i++){
210 int d=0;
211 amp_tmp[0] = 0;
212 amp_tmp[1] = 0;
213 amp_tmp1[0] = 0;
214 amp_tmp1[1] = 0;
215 amp_tmp2[0] = 0;
216 amp_tmp2[1] = 0;
217 tmp1 = 0;
218 tmp2 = 0;
219 temp_PDF = 0;
220 cof[0] = amp[i]*cos(phase[i]);
221 cof[1] = amp[i]*sin(phase[i]);
222 if(modetype[i]==13){
223 //a0(980) and f0(980) mixture
224 double amp_a0[2]={0};
225 double sa0_980=sF0_980;
226 double sKp2[2]={sKp,mass_Pion*mass_Pion};
227 double sKma2[2]={sKm,mass_Eta*mass_Eta};
228 //propagatora0980(mass1[i],sa0_980,sKp2,sKma2,pro);
229 MIP_LineShape(sa0_980, pro);
230 B[0]=barrier(0,sa0_980,sKp,sKm,rRes);
231 temp_PDF=1;
232 tmp1 = temp_PDF*B[0];
233 amp_a0[0] = tmp1*pro[0];
234 amp_a0[1] = tmp1*pro[1];
235 amp_tmp1[0] = amp_a0[0] ;
236 amp_tmp1[1] = amp_a0[1] ;
237 }
238 else if(modetype[i]==0){
239 //K*(892) K+
240 if(g0[i] == 0){
241 pro[0] = 1;
242 pro[1] = 0;
243 }
244 bool neo=true;
245 if(neo){
246 double sBC=SCADot(p23,p23);
247 //if(g0[i] == 1) propagatorRBW(mass1[i],width1[i],sKstr,sKm,sPi,rRes,1,pro);
248 if(g0[i] == 1) propagatorRBWNeoKstr892(mass1[i],width1[i],sKstr,sPi,sKm,rRes,1,pro);
249 B[0]=barrierNeo(1,sKstr,sPi,sKm,rRes,mass1[i]);
250 B[1]=barrierNeoDs(1,sD,sKstr,sKp,rD,mDs,mass1[i]);
251 temp_PDF=(sBC-sPhi1020+((sD-sKp)*(sKm-sPi)/(sKstr)));
252 }
253 else{
254 if(g0[i] == 1) propagatorRBW(mass1[i],width1[i],sKstr,sKm,sPi,rRes,1,pro);
255 //Com_Multi(pro0,pro1,pro);
256 B[0]=barrier(1,sKstr,sKm,sPi,rRes);
257 B[1]=barrier(1,sD,sKstr,sKp,rD);
258 for(int m=0; m<4; m++){
259 for(int j=0; j<4; j++){
260 temp_PDF += t1D1[m]*t1kstr1[j]*G[m][j];
261 }
262 }
263 }
264 tmp1 = temp_PDF*B[0]*B[1];
265 amp_tmp1[0] = tmp1*pro[0];
266 amp_tmp1[1] = tmp1*pro[1];
267 }
268 else if(modetype[i]==1){
269 //Phi1020 Pi v
270 if(g0[i] == 0){
271 pro[0] = 1;
272 pro[1] = 0;
273 }
274 /*if(g0[i] == 1) propagatorRBW(mass1[i],width1[i],sPhi1020,sKm,sKp,rRes,1,pro);
275 B[0]=barrier(1,sPhi1020,sKp,sKm,rRes);
276 B[1]=barrier(1,sD,sPhi1020,sPi,rD);
277 for(int i=0; i<4; i++){
278 for(int j=0; j<4; j++){
279 temp_PDF += t1D2[i]*t1phi1020[j]*G[i][j];
280 }
281 }
282 tmp1 = temp_PDF*B[0]*B[1];*/
283 bool neo=true;
284 if(neo){
285 if(g0[i] == 1) propagatorRBWNeo(mass1[i],width1[i],sPhi1020,sKm,sKp,rRes,1,pro);
286 B[0]=barrierNeo(1,sPhi1020,sKp,sKm,rRes,mass1[i]);
287 B[1]=barrierNeoDs(1,sD,sPhi1020,sPi,rD,mDs,mass1[i]);
288 double sBC=SCADot(p23,p23);
289 temp_PDF=(sKstr-sBC+((sD-sPi)*(sKp-sKm)/(sKstr)));
290 }
291 else{
292 if(g0[i] == 1) propagatorRBW(mass1[i],width1[i],sPhi1020,sKm,sKp,rRes,1,pro);
293 B[0]=barrier(1,sPhi1020,sKp,sKm,rRes);
294 B[1]=barrier(1,sD,sPhi1020,sPi,rD);
295 for(int m=0; m<4; m++){
296 for(int j=0; j<4; j++){
297 temp_PDF += t1D2[m]*t1phi1020[j]*G[m][j];
298 }
299 }
300 }
301 //Com_Multi(pro0,pro1,pro);
302
303 tmp1 = temp_PDF*B[0]*B[1];
304 amp_tmp1[0] = tmp1*pro[0];
305 amp_tmp1[1] = tmp1*pro[1];
306 }
307 else if(modetype[i]==3){
308 //Kstr1430 K S
309 /**
310 * Normal Model
311 * */
312 //if(g0[i] == 1) propagatorRBW(mass1[i],width1[i],sKstr1430,sKm,sPi,rRes,0,pro);
313 //propagatorRBWNeo(mass1[i],width1[i],sKstr1430,sKm,sPi,rRes,0,pro);
314 double sKm2[2]={sKm, mass_EtaP *mass_EtaP};
315 double sPi2[2]={sPi, mass_Kaon *mass_Kaon};
316 propagatorKstr1430(mass1[i],sKstr1430,sKm2,sPi2,pro);
317 //Com_Multi(pro0,pro1,pro);
318 B[0]=barrier(0,sPhi1020,sKp,sKm,rRes);
319 tmp1= 1* B[0];
320 amp_tmp1[0] = tmp1*pro[0];
321 amp_tmp1[1] = tmp1*pro[1];
322
323 /**
324 * KPi-S wave LASS
325 * */
326 /*amp_tmp1[0]=Amp_KPiS1[0];
327 amp_tmp1[1]=Amp_KPiS1[1];*/
328
329
330 }
331 else if(modetype[i]==4){
332 //f0(1710) Pi+ S
333 if(g0[i] == 1) propagatorRBWNeo(mass1[i],width1[i],sF0_1710,sKp,sKm,rRes,0,pro);
334 //if(g0[i] == 1) propagatorRBW(mass1[i],width1[i],sF0_1710,sKp,sKm,rRes,0,pro);
335 if(g0[i] == 0){
336 pro[0] = 1;
337 pro[1] = 0;
338 }
339 //Com_Multi(pro0,pro1,pro);
340 B[0]=barrier(0,sF0_1710,sKp,sKm,rRes);
341 temp_PDF=1;
342 tmp1= temp_PDF* B[0];
343 amp_tmp1[0] = tmp1*pro[0];
344 amp_tmp1[1] = tmp1*pro[1];
345 /*amp_tmp1[0]=Amp_KPiS1[0];
346 amp_tmp1[1]=Amp_KPiS1[1];*/
347 const bool debug89=false;
348 }
349 else if(modetype[i]==5){
350 //f0(1370) Pi+ S
351 if(g0[i] == 1) propagatorRBWNeo(mass1[i],width1[i],sF0_1370,sKp,sKm,rRes,0,pro);
352 //if(g0[i] == 1) propagatorRBW(mass1[i],width1[i],sF0_1370,sKp,sKm,rRes,0,pro);
353 if(g0[i] == 0){
354 pro[0] = 1;
355 pro[1] = 0;
356 }
357 //Com_Multi(pro0,pro1,pro);
358 B[0]=barrier(0,sF0_1370,sKp,sKm,rRes);
359 tmp1= 1* B[0];
360 amp_tmp1[0] = tmp1*pro[0];
361 amp_tmp1[1] = tmp1*pro[1];
362 }
363 amp_tmp[0] = amp_tmp1[0];
364 amp_tmp[1] = amp_tmp1[1];
365 Com_Multi(amp_tmp,cof,amp_PDF);
366 //printf("Line: %u @ %u modetype[%d]=%d amp_tmp[0]=%f amp_tmp[1]=%f cof[0]=%f cof[1]=%f \n",__LINE__,__FILE__,i,modetype[i],amp_tmp[0],amp_tmp[1],cof[0],cof[1]);
367 PDF[0] += amp_PDF[0];
368 PDF[1] += amp_PDF[1];
369 double valTmp=amp_PDF[0] * amp_PDF[0] + amp_PDF[1] * amp_PDF[1];
370 }
371
372 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
373 if(value <=0) {
374 value = 1e-20;
375 }
376 //printf("Line: %u @ %u value=%f\n",__LINE__,__FILE__,value);
377 Result = value;
378}
379
380void EvtDsToKKpi::Com_Multi(double a1[2], double a2[2], double res[2]){
381 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
382 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
383}
384void EvtDsToKKpi::Com_Divide(double a1[2], double a2[2], double res[2]){
385 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/(a2[0]*a2[0]+a2[1]*a2[1]);
386 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/(a2[0]*a2[0]+a2[1]*a2[1]);
387}
388
389double EvtDsToKKpi::SCADot(double a1[4], double a2[4])
390{
391 double _cal = 0.;
392 _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
393 return _cal;
394}
395double EvtDsToKKpi::barrier(int l, double sa, double sb, double sc, double r)
396{
397 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
398 if(q < 0) q = 1e-16;
399 double z;
400 z = q*r*r;
401 double F = 0.0;
402 if(l == 0) F = 1;
403 if(l == 1) F = sqrt(2*z/(1+z));
404 if(l == 2) F = sqrt(13*z*z/(9+3*z+z*z));
405 return F;
406}
407double EvtDsToKKpi::barrierNeo(int l, double sa, double sb, double sc, double r, double mR)
408{
409 double pAB=( (sa-sb-sc)*(sa-sb-sc)/4.0 -(sb*sc))/sa;
410 double pR =( (mR*mR-sb-sc)*(mR*mR-sb-sc)/4.0 -(sb*sc))/(mR*mR);
411 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
412 double zAB= pAB*r*r;
413 double zR = pR *r*r;
414 double F = 0.0;
415 if(l == 0) F = 1;
416 if(l == 1) F = sqrt((1+zR)/(1+zAB));
417 if(l == 2) F = sqrt((9+3*zAB+zAB*zAB)/(9+3*zAB+zAB*zAB));
418 return F;
419}
420double EvtDsToKKpi::barrierNeoDs(int l, double sa, double sb, double sc, double r, double mR,double mb)
421{
422 double pAB=( (sa-sb-sc)*(sa-sb-sc)/4.0 -(sb*sc))/sa;
423 double pR =( (sa-mb*mb-sc)*(sa-mb*mb-sc)/4.0 -(mb*mb*sc))/(mR*mR);
424 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
425 double zAB= pAB*r*r;
426 double zR = pR *r*r;
427 double F = 0.0;
428 if(l == 0) F = 1;
429 if(l == 1) F = sqrt((1+zR)/(1+zAB));
430 if(l == 2) F = sqrt((9+3*zAB+zAB*zAB)/(9+3*zAB+zAB*zAB));
431 return F;
432}
433//------------------spin-------------------------------------------
434void EvtDsToKKpi::calt1(double daug1[4], double daug2[4], double t1[4]){
435 double p, pq;
436 double pa[4], qa[4];
437 for(int i=0; i<4; i++){
438 pa[i] = daug1[i] + daug2[i];
439 qa[i] = daug1[i] - daug2[i];
440 }
441 p = SCADot(pa,pa);
442 pq = SCADot(pa,qa);
443 for(int i=0; i<4; i++){
444 t1[i] = qa[i] - pq/p*pa[i];
445 }
446}
447void EvtDsToKKpi::calt2(double daug1[4], double daug2[4], double t2[4][4]){
448 double p, r;
449 double pa[4], t1[4];
450 calt1(daug1,daug2,t1);
451 r = SCADot(t1,t1);
452 for(int i=0; i<4; i++){
453 pa[i] = daug1[i] + daug2[i];
454 }
455 p = SCADot(pa,pa);
456 for(int i=0; i<4; i++){
457 for(int j=0; j<4; j++){
458 t2[i][j] = t1[i]*t1[j] - 1.0/3*r*(G[i][j]-pa[i]*pa[j]/p);
459 }
460 }
461}
462//-------------------prop--------------------------------------------
463void EvtDsToKKpi::propagator(double mass, double width, double sx, double prop[2]){
464 double a[2], b[2];
465 a[0] = 1;
466 a[1] = 0;
467 b[0] = mass*mass-sx;
468 b[1] = -mass*width;
469 Com_Divide(a,b,prop);
470}
471double EvtDsToKKpi::wid(double mass, double sa, double sb, double sc, double r, int l){
472 double widm = 0.;
473 double q = 0.;
474 double q0 = 0.;
475 double sa0 = mass*mass;
476 double m = sqrt(sa);
477 q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
478 if(q<0) q = 1e-16;
479 q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
480 if(q0<0) q0 = 1e-16;
481 // double F = barrier(l,sa,sb,sc,r);
482 double z = q*r*r;
483 double z0 = q0*r*r;
484 double F = 0;
485 if(l == 0) F = 1;
486 if(l == 1) F = sqrt((1+z0)/(1+z));
487 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
488 double t = q/q0;
489 widm = pow(t,l+0.5)*mass/m*F*F;
490 return widm;
491}
492
493void EvtDsToKKpi::Flatte_rhoab(double sa, double sb, double sc, double rho[2])
494{
495 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
496 double b[2];
497 if(q>0) {
498 rho[0]=2* sqrt(q/sa);
499 rho[1]=0;
500 }
501 else if(q<0){
502 rho[0]=0;
503 rho[1]=2*sqrt(-q/sa);
504 }
505}
506
507void EvtDsToKKpi::propagatorFlatte(double mass, double width, double sx, double *sb, double *sc,double prop[2])
508{
509 double unit[2]={1.0};
510 double ci[2]={0,1};
511 double rho1[2];
512 Flatte_rhoab(sx,sb[0],sc[0],rho1);
513 double rho2[2];
514 Flatte_rhoab(sx,sb[1],sc[1],rho2);
515 double g1_f980=0.165, g2_f980=0.69465;
516 double tmp1[2]={g1_f980,0};
517 double tmp11[2];
518 double tmp2[2]={g2_f980,0};
519 double tmp22[2];
520 Com_Multi(tmp1,rho1,tmp11);
521 Com_Multi(tmp2,rho2,tmp22);
522 double tmp3[2]={tmp11[0]+tmp22[0],tmp11[1]+tmp22[1]};
523 double tmp31[2];
524 Com_Multi(tmp3, ci,tmp31);
525 double tmp4[2]={mass*mass-sx- tmp31[0], -1.0*tmp3[1]};
526 Com_Divide( unit,tmp4, prop);
527}
528void EvtDsToKKpi::propagator980(double mass, double sx, double *sb, double *sc, double prop[2])
529{
530 double unit[2]={1.0};
531 double ci[2]={0,1};
532 double rho1[2];
533 Flatte_rhoab(sx,sb[0],sc[0],rho1);
534 double rho2[2];
535 Flatte_rhoab(sx,sb[1],sc[1],rho2);
536 double gK_f980=0.69466, gPi_f980=0.165;
537 double tmp1[2]={gK_f980,0};
538 double tmp11[2];
539 double tmp2[2]={gPi_f980,0};
540 double tmp22[2];
541 Com_Multi(tmp1,rho1,tmp11);
542 Com_Multi(tmp2,rho2,tmp22);
543 double tmp3[2]={tmp11[0]+tmp22[0],tmp11[1]+tmp22[1]};
544 double tmp31[2];
545 Com_Multi(tmp3, ci,tmp31);
546 double tmp4[2]={mass*mass-sx-tmp31[0], -1.0*tmp31[1]};
547 Com_Divide( unit,tmp4, prop);
548}
549
550void EvtDsToKKpi::propagatora0980(double mass, double sx, double *sb, double *sc, double prop[2])
551{
552 double unit[2]={1.0};
553 double ci[2]={0,1};
554 double rho1[2];
555 Flatte_rhoab(sx,sb[0],sc[0],rho1);
556 double rho2[2];
557 Flatte_rhoab(sx,sb[1],sc[1],rho2);
558 double gKK_a980=0.892*0.341, gPiEta_a980=0.341;
559 double tmp1[2]={gKK_a980,0};
560 double tmp11[2];
561 double tmp2[2]={gPiEta_a980,0};
562 double tmp22[2];
563 Com_Multi(tmp1,rho1,tmp11);
564 Com_Multi(tmp2,rho2,tmp22);
565 double tmp3[2]={tmp11[0]+tmp22[0],tmp11[1]+tmp22[1]};
566 double tmp31[2];
567 Com_Multi(tmp3, ci,tmp31);
568 double tmp4[2]={mass*mass-sx-tmp31[0], -1.0*tmp31[1]};
569 Com_Divide( unit,tmp4, prop);
570}
571
572void EvtDsToKKpi::propagatorKstr1430(double mass, double sx, double *sb, double *sc, double prop[2])
573{
574 double unit[2]={1.0};
575 double ci[2]={0,1};
576 double rho1[2];
577 Flatte_rhoab(sx,sb[0],sc[0],rho1);
578 double rho2[2];
579 Flatte_rhoab(sx,sb[1],sc[1],rho2);
580 double gKPi_Kstr1430=0.2990, gEtaPK_Kstr1430=0.0529;
581 double tmp1[2]={gKPi_Kstr1430,0};
582 double tmp11[2];
583 double tmp2[2]={gEtaPK_Kstr1430,0};
584 double tmp22[2];
585 Com_Multi(tmp1,rho1,tmp11);
586 Com_Multi(tmp2,rho2,tmp22);
587 double tmp3[2]={tmp11[0]+tmp22[0],tmp11[1]+tmp22[1]};
588 double tmp31[2];
589 Com_Multi(tmp3, ci,tmp31);
590 double tmp4[2]={mass*mass-sx-tmp31[0], -1.0*tmp31[1]};
591 Com_Divide( unit,tmp4, prop);
592}
593
594void EvtDsToKKpi::propagatorRBW(double mass, double width, double sa, double sb, double sc, double r, int l, double prop[2]){
595 double a[2], b[2];
596 a[0] = 1;
597 a[1] = 0;
598 b[0] = mass*mass-sa;
599 b[1] = -mass*width*wid(mass,sa,sb,sc,r,l);
600 Com_Divide(a,b,prop);
601}
602
603void EvtDsToKKpi::propagatorRBWNeo(double mass, double width, double sa, double sb, double sc, double r, int l, double prop[2]){
604 double a[2], b[2];
605 a[0] = 1;
606 a[1] = 0;
607 b[0] = mass*mass-sa;
608 double tmp1=(sa-sb-sc);
609 double tmp2=sb*sc;
610 double pAB=sqrt( (tmp1*tmp1/4.0 -tmp2)/sa);
611 double pR =sqrt(( (mass*mass-sb-sc)*(mass*mass-sb-sc)/4.0 -(sb*sc))/(mass*mass));
612 double fR=sqrt(1.0+1.5*1.5*pR*pR)/sqrt(1.0+1.5*1.5*pAB*pAB);
613 double power;
614 if(!l){
615 power=1;
616 fR=1;
617 }
618 else if(l ==1){
619 power=3;
620 }
621 double gammaAB= width*pow(pAB/pR,power)*(mass/sqrt(sa))*fR*fR;
622 b[1] = -mass*gammaAB;
623 Com_Divide(a,b,prop);
624}
625
626void EvtDsToKKpi::propagatorRBWNeoKstr892(double mass, double width, double sa, double sb, double sc, double r, int l, double prop[2]){
627 double a[2], b[2];
628 a[0] = 1;
629 a[1] = 0;
630 b[0] = mass*mass-sa;
631 double tmp1=(sa-sb-sc);
632 double tmp2=sb*sc;
633 double pAB=sqrt( (tmp1*tmp1/4.0 -tmp2)/sa);
634 double pR =sqrt(( (mass*mass-sb-sc)*(mass*mass-sb-sc)/4.0 -(sb*sc))/(mass*mass));
635 double fR=sqrt(1.0+1.5*1.5*pR*pR)/sqrt(1.0+1.5*1.5*pAB*pAB);
636 double power;
637 if(!l){
638 power=1;
639 }
640 else if(l ==1){
641 power=3;
642 }
643 double gammaAB= width*pow(pAB/pR,power)*(mass/sqrt(sa))*fR*fR;
644 b[1] = -mass*gammaAB;
645 Com_Divide(a,b,prop);
646}
647
648//------------GS---used by rho----------------------------
649double EvtDsToKKpi::h(double m, double q){
650 double h = 2/math_pi*q/m*log((m+2*q)/(2*mass_Pion));
651 return h;
652}
653double EvtDsToKKpi::dh(double mass, double q0){
654 double dh = h(mass,q0)*(1.0/(8*q0*q0)-1.0/(2*mass*mass))+1.0/(2*math_pi*mass*mass);
655 return dh;
656}
657double EvtDsToKKpi::f(double mass, double sx, double q0, double q){
658 double m = sqrt(sx);
659 double f = mass*mass/(pow(q0,3))*(q*q*(h(m,q)-h(mass,q0))+(mass*mass-sx)*q0*q0*dh(mass,q0));
660 return f;
661}
662double EvtDsToKKpi::d(double mass, double q0){
663 double d = 3.0/math_pi*mass_Pion*mass_Pion/(q0*q0)*log((mass+2*q0)/(2*mass_Pion))+mass/(2*math_pi*q0)
664 - (mass_Pion*mass_Pion*mass)/(math_pi*pow(q0,3));
665 return d;
666}
667
668//rho
669void EvtDsToKKpi::propagatorGS(double mass, double width, double sa, double sb, double sc, double r, int l, double prop[2]){
670 double a[2], b[2];
671 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
672 double sa0 = mass*mass;
673 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
674 if(q<0) q = 1e-16;
675 if(q0<0) q0 = 1e-16;
676 q = sqrt(q);
677 q0 = sqrt(q0);
678 a[0] = 1+d(mass,q0)*width/mass;
679 a[1] = 0;
680 b[0] = mass*mass-sa+width*f(mass,sa,q0,q);
681 b[1] = -mass*width*wid(mass,sa,sb,sc,r,l);
682 Com_Divide(a,b,prop);
683}
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
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
EvtDsToKKpi::initProbMax
void initProbMax()
Definition EvtDsToKKpi.cc:115
EvtDsToKKpi::decay
void decay(EvtParticle *p)
Definition EvtDsToKKpi.cc:119
EvtDsToKKpi::getName
void getName(std::string &name)
Definition EvtDsToKKpi.cc:35
EvtDsToKKpi::clone
EvtDecayBase * clone()
Definition EvtDsToKKpi.cc:39
EvtDsToKKpi::~EvtDsToKKpi
virtual ~EvtDsToKKpi()
Definition EvtDsToKKpi.cc:33
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