BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToEtapi2pi0.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: EvtDsToEtapi2pi0.cc
11//
12// Description: Routine to handle three-body decays of Ds+/Ds-
13//
14// Modification history:
15// Liaoyuan Dong Sun Jun 02 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/EvtDsToEtapi2pi0.hh"
27#include <string>
28#include "EvtGenBase/EvtConst.hh"
29#include "EvtGenBase/EvtFlatte.hh"
30#include "EvtGenBase/EvtDecayTable.hh"
31using std::endl;
32
33EvtDsToEtapi2pi0::~EvtDsToEtapi2pi0() {}
34
35void EvtDsToEtapi2pi0::getName(std::string& model_name){
36 model_name="DsToEtapi2pi0";
37}
38
42
43void EvtDsToEtapi2pi0::init(){
44
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 phi[0] = 0.0000e+00;
57 phi[1] = 2.2539e+00;
58 phi[2] = -4.1144e+00;
59 rho[0] = 1.0000e+00;
60 rho[1] = 3.0726e-01;
61 rho[2] = 8.8537e-01;
62 mass1[0] = 7.7526e-01;
63 mass1[1] = 9.8000e-01;
64 mass1[2] = 9.8000e-01;
65 mass2[0] = 1.2300e+00;
66 mass2[1] = 7.7526e-01;
67 mass2[2] = 1.8000e+00;
68 width1[0] = 1.4910e-01;
69 width1[1] = 7.0000e-02;
70 width1[2] = 7.0000e-02;
71 width2[0] = 3.3049e-01;
72 width2[1] = 1.4910e-01;
73 width2[2] = 3.8900e-01;
74
75 modetype[0] = 1;
76 modetype[1] = 2;
77 modetype[2] = 3;
78
79 //cout << "DsToEtapi2pi0 :" << endl;
80 //for (int i=0; i<3; i++) {
81 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << " m1= " << mass1[i] << " m2= " << mass2[i] << " w1= " << width1[i] << " w2= " << width2[i] << endl;
82 //}
83
84 rRes = 3.0;
85 rRes2 = 9.0;
86 rD = 5.0;
87 rD2 = 25.0;
88 mass_Ks = 0.497611;
89 mass_Kaon = 0.493677;
90 mass_Pion = 0.13957;
91 mass_Pi0 = 0.1349768;
92 meta = 0.547862;
93 GS1 = 0.636619783;
94 GS2 = 0.01860182466;
95 GS3 = 0.1591549458; // 1/(2*math_2pi)
96 GS4 = 0.00620060822;
97
98 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
99 for (int a=0; a<4; a++) {
100 for (int j=0; j<4; j++) {
101 G[a][j] = GG[a][j];
102 }
103 }
104}
105
106void EvtDsToEtapi2pi0::initProbMax() {
107 setProbMax(4400.0);
108}
109
110void EvtDsToEtapi2pi0::decay( EvtParticle *p){
111/*
112 double maxprob = 0.0;
113 for(int ir=0;ir<=60000000;ir++){
114 p->initializePhaseSpace(getNDaug(),getDaugs());
115 EvtVector4R pic = p->getDaug(0)->getP4();
116 EvtVector4R pi01 = p->getDaug(1)->getP4();
117 EvtVector4R pi02 = p->getDaug(2)->getP4();
118 EvtVector4R eta = p->getDaug(3)->getP4();
119
120 double Pic[4],Pi01[4],Pi02[4],Eta[4];
121 Pic[0]=pic.get(0); Pic[1]=pic.get(1); Pic[2]=pic.get(2); Pic[3]=pic.get(3);
122 Pi01[0]=pi01.get(0); Pi01[1]=pi01.get(1); Pi01[2]=pi01.get(2); Pi01[3]=pi01.get(3);
123 Pi02[0]=pi02.get(0); Pi02[1]=pi02.get(1); Pi02[2]=pi02.get(2); Pi02[3]=pi02.get(3);
124 Eta[0] = eta.get(0); Eta[1] = eta.get(1); Eta[2] = eta.get(2); Eta[3] = eta.get(3);
125
126 double prob;
127 int g0[3]={1,1,1};
128 int g1[3]={1,1,1};
129 int g2[3]={0,1,0};
130 int nstates=3;
131
132 calEvaMy(Pic, Pi01, Pi02, Eta, mass1, mass2, width1, width2, rho, phi, g0, g1, g2, modetype, nstates, prob);
133
134 if(prob>maxprob) {
135 maxprob=prob;
136 printf("nrun = %d,maxprob = %.10f,prob = %.10f\n\n",ir,maxprob,prob);
137 }
138 }
139*/
140
141 p->initializePhaseSpace(getNDaug(),getDaugs());
142 EvtVector4R pic = p->getDaug(0)->getP4();
143 EvtVector4R pi01 = p->getDaug(1)->getP4();
144 EvtVector4R pi02 = p->getDaug(2)->getP4();
145 EvtVector4R eta = p->getDaug(3)->getP4();
146
147 double Pic[4],Pi01[4],Pi02[4],Eta[4];
148 Pic[0]=pic.get(0); Pic[1]=pic.get(1); Pic[2]=pic.get(2); Pic[3]=pic.get(3);
149 Pi01[0]=pi01.get(0); Pi01[1]=pi01.get(1); Pi01[2]=pi01.get(2); Pi01[3]=pi01.get(3);
150 Pi02[0]=pi02.get(0); Pi02[1]=pi02.get(1); Pi02[2]=pi02.get(2); Pi02[3]=pi02.get(3);
151 Eta[0]=eta.get(0); Eta[1]=eta.get(1); Eta[2]=eta.get(2); Eta[3]=eta.get(3);
152
153 double prob;
154 int g0[3]={1,1,1};
155 int g1[3]={1,1,1};
156 int g2[3]={0,1,0};
157 int nstates=3;
158
159 calEvaMy(Pic, Pi01, Pi02, Eta, mass1, mass2, width1, width2, rho, phi, g0, g1, g2, modetype, nstates, prob);
160 setProb(prob);
161
162 return ;
163}
164
165void EvtDsToEtapi2pi0::calEvaMy(double* Pic, double* Pi01, double* Pi02, double* Eta, double *mass1, double *mass2, double* width1, double* width2, double *amp, double *phase,int* g0, int* g1, int* g2, int* modetype, int nstates, double & Result){
166
167 double P12[4],P13[4],P14[4],P23[4],P24[4],P34[4],P123[4],P124[4],P134[4],P234[4],PD[4];
168 for(int i=0;i!=4;i++) {
169 P12[i] = Pic[i] + Pi01[i];
170 P13[i] = Pic[i] + Pi02[i];
171 P14[i] = Pic[i] + Eta[i];
172 P23[i] = Pi01[i] + Pi02[i];
173 P24[i] = Pi01[i] + Eta[i];
174 P34[i] = Pi02[i] + Eta[i];
175 P123[i] = Pic[i] + Pi01[i] + Pi02[i];
176 P124[i] = Pic[i] + Pi01[i] + Eta[i];
177 P134[i] = Pic[i] + Pi02[i] + Eta[i];
178 P234[i] = Pi01[i] + Pi02[i] + Eta[i];
179 PD[i] = P123[i] + Eta[i];
180 }
181
182 double cof[2], amp_tmp1[2], amp_tmp2[2], amp_tmp[2], amp_PDF[2], PDF[2];
183 double s124,s134,t1f11[4],t1f14[4],t1f15[4],t1f16[4],t1f17[4];
184 double spic,spi01,spi02,seta,s12,s13,s14,s23,s24,s34,s123,s234,sD;
185
186 spic = SCADot(Pic,Pic);
187 spi01 = SCADot(Pi01,Pi01);
188 spi02 = SCADot(Pi02,Pi02);
189 seta = SCADot(Eta,Eta);
190 s12 = SCADot(P12,P12);
191 s13 = SCADot(P13,P13);
192 s14 = SCADot(P14,P14);
193 s23 = SCADot(P23,P23);
194 s24 = SCADot(P24,P24);
195 s34 = SCADot(P34,P34);
196 s123 = SCADot(P123,P123);
197 s124 = SCADot(P124,P124);
198 s134 = SCADot(P134,P134);
199 s234 = SCADot(P234,P234);
200 sD = SCADot(PD,PD);
201
202
203 double seta2[2]={mass_Ks*mass_Ks,seta};
204 double spion12[2]={mass_Kaon*mass_Kaon,spi01};
205 double spion22[2]={mass_Kaon*mass_Kaon,spi02};
206 double spim2[2]={mass_Kaon*mass_Kaon,spic};
207
208 double t1rho1[4],t1rho2[4],t1a01[4],t1a02[4],t1d01[4],t1d02[4],t1D[4],t1a1[4],t1sigma[4],t1d03[4],t2a11[4][4],t2a12[4][4];
209 double t1f12[4],t1f13[4];
210
211 calt1(Pic,Pi01,t1rho1);
212 calt1(Pic,Pi02,t1rho2);
213
214 calt1(Pi01,Eta,t1a01);
215 calt1(Pi02,Eta,t1a02);
216
217 calt1(P24,P13,t1d01);
218 calt1(P34,P12,t1d02);
219
220 calt1(P14,P23,t1d03);
221 calt1(Pi01,Pi02,t1sigma);
222
223
224 calt2(P12,Pi02,t2a11);
225 calt2(P13,Pi01,t2a12);
226//\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
227
228 PDF[0] = 0.0;
229 PDF[1] = 0.0;
230 double mass1sq, mass2sq;
231 double pro[2],pro0[2],pro1[2],pro2[2],pro3[2];
232 double B[3];
233 amp_PDF[0] = 0;
234 amp_PDF[1] = 0;
235 PDF[0] = 0;
236 PDF[1] = 0;
237
238 double temp_PDF, tmp1, tmp2;
239
240 for(int i=0; i<nstates; i++){
241
242 amp_tmp[0] = 0;
243 amp_tmp[1] = 0;
244 mass1sq = mass1[i]*mass1[i];
245 mass2sq = mass2[i]*mass2[i];
246 amp_tmp1[0] = 0;
247 amp_tmp1[1] = 0;
248 amp_tmp2[0] = 0;
249 amp_tmp2[1] = 0;
250 temp_PDF = 0;
251 tmp1 = 0;
252 tmp2 = 0;
253
254
255 cof[0] = amp[i]*cos(phase[i]);
256 cof[1] = amp[i]*sin(phase[i]);
257
258 if(modetype[i]==1)
259 {
260 if(g0[i]==1) propagatorGS(mass1sq,mass1[i],width1[i],s12,spic,spi01,rRes2,pro0);//rho
261 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; }
262 if(g1[i]==1) propagatorRBW(mass2sq,mass2[i],width2[i],s123,s12,spi02,rRes2,g2[i],pro1);//a1
263 if(g1[i]==0) { pro1[0] = 1; pro1[1] = 0; }
264 Com_Multi(pro0,pro1,pro);//Pn
265 calt1(P123,Eta,t1D);//L1(D)
266
267 B[0] = barrier(1,s12,spi01,spic,rRes2,mass1[i]);//Fn(rho)
268 B[2] = barrier(1,sD,s123,seta,rD2,1.9683);//Fn(D)
269
270 if(g2[i]==0)//a1--rhopi
271 {
272
273 for(int a=0; a<4; a++)
274 {
275 for(int j=0; j<4; j++)
276 {
277 temp_PDF += t1D[a]*(G[a][j]-P123[a]*P123[j]/s123)*t1rho1[j]*G[a][a]*G[j][j];
278 }
279 }
280 tmp1 = B[0]*B[2]*temp_PDF;//
281 }
282
283 amp_tmp1[0] = tmp1*pro[0];
284 amp_tmp1[1] = tmp1*pro[1];
285
286 temp_PDF=0;
287 if(g0[i]==1) propagatorGS(mass1sq,mass1[i],width1[i],s13,spic,spi02,rRes2,pro0);
288 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; }
289 if(g1[i]==1) propagatorRBW(mass2sq,mass2[i],width2[i],s123,s13,spi01,rRes2,g2[i],pro1);
290 if(g1[i]==0) { pro1[0] = 1; pro1[1] = 0; }
291 Com_Multi(pro0,pro1,pro);
292 calt1(P123,Eta,t1D);
293 B[0] = barrier(1,s13,spi02,spic,rRes2,mass1[i]);
294 B[2] = barrier(1,sD,s123,seta,rD2,1.9683);
295
296 if(g2[i]==0)//a1--rhopi
297 {
298 for(int a=0; a<4; a++)
299 {
300 for(int j=0; j<4; j++)
301 {
302 temp_PDF += t1D[a]*(G[a][j]-P123[a]*P123[j]/s123)*t1rho2[j]*G[a][a]*G[j][j];
303 }
304 }
305 tmp2 = B[0]*B[2]*temp_PDF;
306 }
307 amp_tmp2[0] = tmp2*pro[0];
308 amp_tmp2[1] = tmp2*pro[1];
309
310 }
311
312 if(modetype[i] == 2){
313
314 if(g0[i]==1) propagatora0980(mass1[i],s24,seta2,spion12,pro0);
315 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; }
316 if(g1[i]==1) propagatorGS(mass2sq,mass2[i],width2[i],s13,spic,spi02,rRes2,pro1);
317 if(g1[i]==0) { pro1[0] = 1; pro1[1] = 0; }
318 Com_Multi(pro0,pro1,pro);//Pn
319
320 for(int a=0; a<4; a++)
321 {
322 temp_PDF += G[a][a]*t1d01[a]*t1rho2[a];//Sn=L1(D)*L1(V:rho)
323 }
324
325 B[1] = barrier(1,s13,spi02,spic,rRes2,mass2[i]);//Fn1(rho)
326 B[2] = barrier(1,sD,s24,s13,rD2,1.9683);//Fn1(D)...Fn0(a0)==1
327
328 tmp1 = B[1]*B[2]*temp_PDF;//Fn*Sn
329
330 amp_tmp1[0] = tmp1*pro[0];//An
331 amp_tmp1[1] = tmp1*pro[1];
332
333 temp_PDF=0;
334 if(g0[i]==1) propagatora0980(mass1[i],s34,seta2,spion22,pro0);
335 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; }
336 if(g1[i]==1) propagatorGS(mass2sq,mass2[i],width2[i],s12,spic,spi01,rRes2,pro1);
337 if(g1[i]==0) { pro1[0] = 1; pro1[1] = 0; }
338 Com_Multi(pro0,pro1,pro);
339
340 for(int a=0; a<4; a++)
341 {
342 temp_PDF += G[a][a]*t1d02[a]*t1rho1[a];
343 }
344
345 B[1] = barrier(1,s12,spi01,spic,rRes2,mass2[i]);//Fn1(rho)
346 B[2] = barrier(1,sD,s34,s12,rD2,1.9683);
347
348 tmp2 = B[1]*B[2]*temp_PDF;
349 amp_tmp2[0] = tmp2*pro[0];
350 amp_tmp2[1] = tmp2*pro[1];
351
352 }
353
354 if(modetype[i]==3)
355 {
356 if(g0[i]==1) propagatora0980(mass1[i],s14,spion12,seta2,pro0);
357 if(g0[i]==0) { pro0[0] = 1; pro0[1] = 0; pro2[0]=1; pro2[1]=0;}
358 if(g1[i]==1) { pro1[0] = 1; pro1[1] = 0; pro3[0]=1; pro3[1]=0;}
359 Com_Multi(pro0,pro1,pro);
360 tmp1 = 1;
361 amp_tmp1[0] = tmp1*pro[0];
362 amp_tmp1[1] = tmp1*pro[1];
363
364 amp_tmp2[0] = amp_tmp1[0];
365 amp_tmp2[1] = amp_tmp1[1];
366
367 }
368
369 amp_tmp[0] = amp_tmp1[0]+amp_tmp2[0];
370 amp_tmp[1] = amp_tmp1[1]+amp_tmp2[1];
371 Com_Multi(amp_tmp,cof,amp_PDF);
372 PDF[0] += amp_PDF[0];
373 PDF[1] += amp_PDF[1];
374 }
375
376 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
377 Result = value;
378}
379
380void EvtDsToEtapi2pi0::Com_Multi(double a1[2], double a2[2], double res[2])
381{
382 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
383 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
384}
385void EvtDsToEtapi2pi0::Com_Divide(double a1[2], double a2[2], double res[2])
386{
387 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
388 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
389 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
390}
391double EvtDsToEtapi2pi0::SCADot(double a1[4], double a2[4])
392{
393 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
394 return _cal;
395}
396
397double EvtDsToEtapi2pi0::barrier(int l, double sa, double sb, double sc, double r, double mass)
398{
399 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
400 if(q < 0) q = 1e-16;
401 double z;
402 z = q*r;
403 double sa0;
404 sa0 = mass*mass;
405 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
406 if(q0 < 0) q0 = 1e-16;
407 //if(q0 < 0) q0 = -1.0*q0;
408 double z0 = q0*r;
409 double F = 0.0;
410 if(l == 0) F = 1;
411 if(l == 1) F = sqrt((1.0+z0)/(1.0+z));
412 if(l == 2) F = sqrt((9.0+3.0*z0+z0*z0)/(9.0+3.0*z+z*z));
413 //if(l == 1) F = sqrt((2.0*z)/(1+z));
414 //if(l == 2) F = sqrt((13.0*z*z)/(9+3*z+z*z));
415 return F;
416}
417
418void EvtDsToEtapi2pi0::calt1(double daug1[4], double daug2[4], double t1[4])
419{
420 double p, pq, tmp;
421 double pa[4], qa[4];
422 for(int i=0; i<4; i++) {
423 pa[i] = daug1[i] + daug2[i];
424 qa[i] = daug1[i] - daug2[i];
425 }
426 p = SCADot(pa,pa);
427 pq = SCADot(pa,qa);
428 tmp = pq/p;
429 for(int i=0; i<4; i++) {
430 t1[i] = qa[i] - tmp*pa[i];
431 }
432}
433
434void EvtDsToEtapi2pi0::calt2(double daug1[4], double daug2[4], double t2[4][4])
435{
436 double p, r;
437 double pa[4], t1[4];
438 calt1(daug1,daug2,t1);
439 r = SCADot(t1,t1)/3.0;
440 for(int i=0; i<4; i++) {
441 pa[i] = daug1[i] + daug2[i];
442 }
443 p = SCADot(pa,pa);
444 for(int i=0; i<4; i++) {
445 for(int j=0; j<4; j++) {
446 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
447 }
448 }
449}
450
451void EvtDsToEtapi2pi0::propagator(double mass2, double mass, double width, double sx, double prop[2])
452{
453 double a[2], b[2];
454 a[0] = 1;
455 a[1] = 0;
456 b[0] = mass2-sx;
457 b[1] = -mass*width;
458 Com_Divide(a,b,prop);
459}
460double EvtDsToEtapi2pi0::wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l)
461{
462 double widm = 0.;
463 double m = sqrt(sa);
464 double tmp = sb-sc;
465 double tmp1 = sa+tmp;
466 double q = fabs(0.25*tmp1*tmp1/sa-sb);
467 double tmp2 = mass2+tmp;
468 double q0 = fabs(0.25*tmp2*tmp2/mass2-sb);
469 double z = q*r2;
470 double z0 = q0*r2;
471 double t = q/q0;
472 if(l == 0) {widm = sqrt(t)*mass/m;}
473 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
474 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
475 return widm;
476
477}
478double EvtDsToEtapi2pi0::widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
479{
480 double widm = 0.;
481 double m = sqrt(sa);
482 double tmp = sb-sc;
483 double tmp1 = sa+tmp;
484 double q = fabs(0.25*tmp1*tmp1/sa-sb);
485 double tmp2 = mass2+tmp;
486 double q0 = fabs(0.25*tmp2*tmp2/mass2-sb);
487 double z = q*r2;
488 double z0 = q0*r2;
489 double F = (1+z0)/(1+z);
490 double t = q/q0;
491 widm = t*sqrt(t)*mass/m*F;
492 return widm;
493
494}
495void EvtDsToEtapi2pi0::propagatorRBW(double mass2, double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
496{
497 double a[2], b[2];
498 a[0] = 1;
499 a[1] = 0;
500 b[0] = mass2-sa;
501 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
502 Com_Divide(a,b,prop);
503}
504
505
506void EvtDsToEtapi2pi0::propagatorGS(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
507{
508 double a[2], b[2];
509 double tmp = sb-sc;
510 double tmp1 = sa+tmp;
511 double q2 = fabs(0.25*tmp1*tmp1/sa-sb);
512
513 double tmp2 = mass2+tmp;
514 double q02 = fabs(0.25*tmp2*tmp2/mass2-sb);
515
516 double q = sqrt(q2);
517 double q0 = sqrt(q02);
518 double m = sqrt(sa);
519 double q03 = q0*q02;
520 double tmp3 = log(mass+2*q0)+1.2760418309; // log(mass_2Pion) = 1.2760418309;
521
522 double h = GS1*q/m*(log(m+2*q)+1.2760418309);
523 double h0 = GS1*q0/mass*tmp3;
524 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
525 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
526 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
527
528 a[0] = 1.0+d*width/mass;
529 a[1] = 0.0;
530 b[0] = mass2-sa+width*f;
531 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
532 Com_Divide(a,b,prop);
533}
534
535void EvtDsToEtapi2pi0::rhoab(double sa, double sb, double sc, double res[2]) {
536 double tmp = sa+sb-sc;
537 double q = 0.25*tmp*tmp/sa-sb;
538 if(q>=0) {
539 res[0]=2.0*sqrt(q/sa);
540 res[1]=0.0;
541 } else {
542 res[0]=0.0;
543 res[1]=2.0*sqrt(-q/sa);
544 }
545}
546
547void EvtDsToEtapi2pi0::rho4Pi(double sa, double res[2]) {
548 double temp = 1.0-0.3116765584/sa; // 0.3116765584=0.13957*0.13957*16
549 if(temp>=0) {
550 res[0]=sqrt(temp)/(1.0+exp(9.8-3.5*sa));
551 res[1]=0.0;
552 } else {
553 res[0]=0.0;
554 res[1]=sqrt(-temp)/(1.0+exp(9.8-3.5*sa));
555 }
556}
557
558void EvtDsToEtapi2pi0::Flatte_rhoab(double sa, double sb, double sc, double rho[2])
559{
560 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
561 if(q>0) {
562 rho[0]=2* sqrt(q/sa);
563 rho[1]=0;
564 }
565 else if(q<0){
566 rho[0]=0;
567 rho[1]=2*sqrt(-q/sa);
568 }
569}
570void EvtDsToEtapi2pi0::propagatora0980(double mass, double sx, double *sb, double *sc, double prop[2])
571{
572 double unit[2]={1.0};
573 double ci[2]={0,1};
574 double rho1[2];
575 Flatte_rhoab(sx,sb[0],sc[0],rho1);
576 double rho2[2];
577 Flatte_rhoab(sx,sb[1],sc[1],rho2);
578
579 double gKK_a980=0.892*0.341, gPiEta_a980=0.341;
580 double tmp1[2]={gKK_a980,0};
581 double tmp11[2];
582 double tmp2[2]={gPiEta_a980,0};
583 double tmp22[2];
584 Com_Multi(tmp1,rho1,tmp11);
585 Com_Multi(tmp2,rho2,tmp22);
586 double tmp3[2]={tmp11[0]+tmp22[0],tmp11[1]+tmp22[1]};
587 double tmp31[2];
588 Com_Multi(tmp3, ci,tmp31);
589 double tmp4[2]={mass*mass-sx-tmp31[0], -1.0*tmp31[1]};
590 Com_Divide( unit,tmp4, prop);
591}
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
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:252
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
EvtDsToEtapi2pi0::decay
void decay(EvtParticle *p)
Definition EvtDsToEtapi2pi0.cc:110
EvtDsToEtapi2pi0::getName
void getName(std::string &name)
Definition EvtDsToEtapi2pi0.cc:35
EvtDsToEtapi2pi0::~EvtDsToEtapi2pi0
virtual ~EvtDsToEtapi2pi0()
Definition EvtDsToEtapi2pi0.cc:33
EvtDsToEtapi2pi0::clone
EvtDecayBase * clone()
Definition EvtDsToEtapi2pi0.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
s12
double double double double * s12
Definition qcdloop1.h:77
s23
double double double double double * s23
Definition qcdloop1.h:77
b
const double b
Definition slope.cxx:9