BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToEtappipi0.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: EvtDsToEtappipi0.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/EvtDsToEtappipi0.hh"
27#include <string>
28#include "EvtGenBase/EvtConst.hh"
29#include "EvtGenBase/EvtDecayTable.hh"
30using std::endl;
31
32EvtDsToEtappipi0::~EvtDsToEtappipi0() {}
33
34void EvtDsToEtappipi0::getName(std::string& model_name){
35 model_name="DsToEtappipi0";
36}
37
41
42void EvtDsToEtappipi0::init(){
43 // check that there are 0 arguments
44 checkNArg(0);
45 checkNDaug(3);
46 checkSpinParent(EvtSpinType::SCALAR);
47 checkSpinDaughter(0,EvtSpinType::SCALAR);
48 checkSpinDaughter(1,EvtSpinType::SCALAR);
49 checkSpinDaughter(2,EvtSpinType::SCALAR);
50
51 phi[0] = 0;
52 rho[0] = 1;
53 phi[1] = 0;
54 rho[1] = 0;
55 phi[2] = 0;
56 rho[2] = 0;
57 phi[3] = 0;
58 rho[3] = 0;
59 phi[4] = 0;
60 rho[4] = 0;
61 phi[5] = 0;
62 rho[5] = 0;
63
64 phi[0] = 0.0;
65 rho[0] = 1.0;
66
67 modetype[0]= 23;
68
69 //cout << "DsToEtappipi0 :" << endl;
70 //for (int i=0; i<6; i++) {
71 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
72 //}
73
74 width[0] = 0.1491;
75 mass[0] = 0.77526;
76
77 mD = 1.86486;
78 mDs = 1.9683;
79 rRes = 9.0;
80 rD = 5.0;
81 metap = 0.95778;
82 mkstr = 0.89594;
83 mk0 = 0.497614;
84 mass_Kaon = 0.49368;
85 mass_Pion = 0.13957;
86 mass_Pi0 = 0.1349766;
87 math_pi = 3.1415926;
88 ma0 = 0.99;
89 Ga0 = 0.0756;
90 meta = 0.547862;
91
92 GS1 = 0.636619783;
93 GS2 = 0.01860182466;
94 GS3 = 0.1591549458; // 1/(2*math_2pi)
95 GS4 = 0.00620060822; // mass_Pion2/math_pi
96
97 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
98 for (int i=0; i<4; i++) {
99 for (int j=0; j<4; j++) {
100 G[i][j] = GG[i][j];
101 }
102 }
103}
104
108
109void EvtDsToEtappipi0::decay( EvtParticle *p){
110 /*
111 double maxprob = 0.0;
112 for(int ir=0;ir<=40000000;ir++){
113 p->initializePhaseSpace(getNDaug(),getDaugs());
114 EvtVector4R D1 = p->getDaug(0)->getP4();
115 EvtVector4R D2 = p->getDaug(1)->getP4();
116 EvtVector4R D3 = p->getDaug(2)->getP4();
117
118 double P1[4], P2[4], P3[4];
119 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
120 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
121 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
122
123 double value;
124 //value = calDalEva(P1, P2, P3);
125 int g0[6]={1};
126 int spin[6]={1};
127 int nstates=1;
128 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
129 if (value<0) continue;
130 if(value>maxprob) {
131 maxprob=value;
132 cout << "ir= " << ir << endl;
133 cout << "double P1[4] = {" << P1[0] <<","<< P1[1] <<","<< P1[2] <<","<< P1[3] <<"};"<< endl;
134 cout << "double P2[4] = {" << P2[0] <<","<< P2[1] <<","<< P2[2] <<","<< P2[3] <<"};"<< endl;
135 cout << "double P3[4] = {" << P3[0] <<","<< P3[1] <<","<< P3[2] <<","<< P3[3] <<"};"<< endl;
136 cout << "MAX====> " << maxprob << endl;
137 }
138 }
139 printf("MAXprob = %.10f\n",maxprob);
140 */
141 p->initializePhaseSpace(getNDaug(),getDaugs());
142 EvtVector4R D1 = p->getDaug(0)->getP4();
143 EvtVector4R D2 = p->getDaug(1)->getP4();
144 EvtVector4R D3 = p->getDaug(2)->getP4();
145
146 double P1[4], P2[4], P3[4];
147 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
148 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
149 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
150
151 double value;
152 int g0[6]={1};
153 int spin[6]={1};
154 int nstates=1;
155
156 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
157
158 setProb(value);
159 return ;
160
161}
162
163void EvtDsToEtappipi0::Com_Multi(double a1[2], double a2[2], double res[2])
164{
165 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
166 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
167}
168void EvtDsToEtappipi0::Com_Divide(double a1[2], double a2[2], double res[2])
169{
170 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
171 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
172 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
173}
174//------------base---------------------------------
175double EvtDsToEtappipi0::SCADot(double a1[4], double a2[4])
176{
177 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
178 return _cal;
179}
180double EvtDsToEtappipi0::barrier(int l, double sa, double sb, double sc, double r, double mass)
181{
182 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
183 if(q < 0) q = 1e-16;
184 double z;
185 z = q*r*r;
186 double sa0;
187 sa0 = mass*mass;
188 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
189 if(q0 < 0) q0 = 1e-16;
190 double z0 = q0*r*r;
191 double F = 0.0;
192 if(l == 0) F = 1;
193 if(l == 1) F = sqrt((1+z0)/(1+z));
194 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
195 return F;
196}
197
198void EvtDsToEtappipi0::calt1(double daug1[4], double daug2[4], double t1[4])
199{
200 double p, pq, tmp;
201 double pa[4], qa[4];
202 for(int i=0; i<4; i++) {
203 pa[i] = daug1[i] + daug2[i];
204 qa[i] = daug1[i] - daug2[i];
205 }
206 p = SCADot(pa,pa);
207 pq = SCADot(pa,qa);
208 tmp = pq/p;
209 for(int i=0; i<4; i++) {
210 t1[i] = qa[i] - tmp*pa[i];
211 }
212}
213void EvtDsToEtappipi0::calt2(double daug1[4], double daug2[4], double t2[4][4])
214{
215 double p, r;
216 double pa[4], t1[4];
217 calt1(daug1,daug2,t1);
218 r = SCADot(t1,t1)/3.0;
219 for(int i=0; i<4; i++) {
220 pa[i] = daug1[i] + daug2[i];
221 }
222 p = SCADot(pa,pa);
223 for(int i=0; i<4; i++) {
224 for(int j=0; j<4; j++) {
225 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
226 }
227 }
228}
229//-------------------prop--------------------------------------------
230
231double EvtDsToEtappipi0::wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l)
232{
233 double widm = 0.;
234 double m = sqrt(sa);
235 double tmp = sb-sc;
236 double tmp1 = sa+tmp;
237 double q = 0.25*tmp1*tmp1/sa-sb;
238 if(q<0) q = 1e-16;
239 double tmp2 = mass2+tmp;
240 double q0 = 0.25*tmp2*tmp2/mass2-sb;
241 if(q0<0) q0 = 1e-16;
242 double z = q*r2;
243 double z0 = q0*r2;
244 double t = q/q0;
245 if(l == 0) {widm = sqrt(t)*mass/m;}
246 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
247 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
248 return widm;
249}
250double EvtDsToEtappipi0::widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
251{
252 double widm = 0.;
253 double m = sqrt(sa);
254 double tmp = sb-sc;
255 double tmp1 = sa+tmp;
256 double q = 0.25*tmp1*tmp1/sa-sb;
257 if(q<0) q = 1e-16;
258 double tmp2 = mass2+tmp;
259 double q0 = 0.25*tmp2*tmp2/mass2-sb;
260 if(q0<0) q0 = 1e-16;
261 double z = q*r2;
262 double z0 = q0*r2;
263 double F = (1+z0)/(1+z);
264 double t = q/q0;
265 widm = t*sqrt(t)*mass/m*F;
266 return widm;
267}
268void EvtDsToEtappipi0::propagatorRBW(double mass2, double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
269{
270 double a[2], b[2];
271 a[0] = 1;
272 a[1] = 0;
273 b[0] = mass2-sa;
274 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
275 Com_Divide(a,b,prop);
276}
277
278void EvtDsToEtappipi0::propagatorGS(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
279{
280 double a[2], b[2];
281 double tmp = sb-sc;
282 double tmp1 = sa+tmp;
283 double q2 = 0.25*tmp1*tmp1/sa-sb;
284 if(q2<0) q2 = 1e-16;
285
286 double tmp2 = mass2+tmp;
287 double q02 = 0.25*tmp2*tmp2/mass2-sb;
288 if(q02<0) q02 = 1e-16;
289
290 double q = sqrt(q2);
291 double q0 = sqrt(q02);
292 double m = sqrt(sa);
293 double q03 = q0*q02;
294 double tmp3 = log(mass+2*q0)+1.2926305904; // log(mpi0+mpip) = -1.2926305904;
295
296 double h = GS1*q/m*(log(m+2*q)+1.2926305904);
297 double h0 = GS1*q0/mass*tmp3;
298 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
299 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
300 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
301
302 a[0] = 1.0+d*width/mass;
303 a[1] = 0.0;
304 b[0] = mass2-sa+width*f;
305 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
306 Com_Divide(a,b,prop);
307}
308
309void EvtDsToEtappipi0::PiPiSWASS(double sa, double sb, double sc, double prop[2]) {
310 double tmp = sb-sc;
311 double tmp2 = sa+tmp;
312 double qs = 0.25*tmp2*tmp2/sa-sb;
313 double q = sqrt(qs);
314 double a0 = -0.11/mass_Pion;
315 prop[0] = 1/(1+a0*a0*q*q);
316 prop[1] = a0*q/(1+a0*a0*q*q);
317}
318
319void EvtDsToEtappipi0::Flatte_rhoab(double sa, double sb, double rho[2])
320{
321 double q = 1.0-(4*sb/sa);
322
323 if(q>0) {
324 rho[0]=sqrt(q);
325 rho[1]=0;
326 }
327 else if(q<0){
328 rho[0]=0;
329 rho[1]=sqrt(-q);
330 }
331}
332
333void EvtDsToEtappipi0:: propagator980(double mass, double sx, double *sb, double prop[2])
334{
335
336 double gpipi1 = 2.0/3.0*0.165;
337 double gpipi2 = 1.0/3.0*0.165;
338 double gKK1 = 0.5*0.69466;
339 double gKK2 = 0.5*0.69466;
340
341 double unit[2]={1.0};
342 double ci[2]={0,1};
343 double rho1[2];
344 Flatte_rhoab(sx,sb[0],rho1);
345 double rho2[2];
346 Flatte_rhoab(sx,sb[1],rho2);
347 double rho3[2];
348 Flatte_rhoab(sx,sb[2],rho3);
349 double rho4[2];
350 Flatte_rhoab(sx,sb[3],rho4);
351
352 double tmp1[2]={gpipi1,0};
353 double tmp11[2];
354 double tmp2[2]={gpipi2,0};
355 double tmp22[2];
356 double tmp3[2]={gKK1,0};
357 double tmp33[2];
358 double tmp4[2]={gKK2,0};
359 double tmp44[2];
360
361 Com_Multi(tmp1,rho1,tmp11);
362 Com_Multi(tmp2,rho2,tmp22);
363 Com_Multi(tmp3,rho3,tmp33);
364 Com_Multi(tmp4,rho4,tmp44);
365
366 double tmp5[2]={tmp11[0]+tmp22[0]+tmp33[0]+tmp44[0],tmp11[1]+tmp22[1]+tmp33[1]+tmp44[1]};
367 double tmp51[2];
368 Com_Multi(tmp5, ci,tmp51);
369 double tmp6[2]={mass*mass-sx-tmp51[0], -1.0*tmp51[1]};
370
371 Com_Divide( unit,tmp6, prop);
372}
373
374
375
376double EvtDsToEtappipi0::DDalitz(double P1[4], double P2[4], double P3[4], int Ang, double mass){
377 double pR[4], pD[4];
378 double temp_PDF, v_re;
379 temp_PDF = 0.0;
380 v_re = 0.0;
381 double B[2], s1, s2, s3, sR, sD;
382 for(int i=0; i<4; i++){
383 pR[i] = P1[i] + P2[i];
384 pD[i] = pR[i] + P3[i];
385 }
386 s1 = SCADot(P1,P1);
387 s2 = SCADot(P2,P2);
388 s3 = SCADot(P3,P3);
389 sR = SCADot(pR,pR);
390 sD = SCADot(pD,pD);
391 //int G[4][4];
392 //for(int i=0; i!=4; i++){
393 // for(int j=0; j!=4; j++){
394 // if(i==j){
395 // if(i==0) G[i][j] = 1;
396 // else G[i][j] = -1;
397 // }
398 // else G[i][j] = 0;
399 // }
400 //}
401 if(Ang == 0){
402 B[0] = 1;
403 B[1] = 1;
404 temp_PDF = 1;
405 }
406 if(Ang == 1){
407 B[0] = barrier(1,sR,s1,s2,3.0,mass);
408 B[1] = barrier(1,sD,sR,s3,5.0,1.9683);
409 double t1[4], T1[4];
410 calt1(P1,P2,t1);
411 calt1(pR,P3,T1);
412 temp_PDF = 0;
413 for(int i=0; i<4; i++){
414 temp_PDF += t1[i]*T1[i]*G[i][i];
415 }
416 }
417 if(Ang == 2){
418 B[0] = barrier(2,sR,s1,s2,3.0,mass);
419 B[1] = barrier(2,sD,sR,s3,5.0,1.9683);
420 double t2[4][4], T2[4][4];
421 calt2(P1,P2,t2);
422 calt2(pR,P3,T2);
423 temp_PDF = 0;
424 for(int i=0; i<4; i++){
425 for(int j=0; j<4; j++){
426 temp_PDF += t2[i][j]*T2[j][i]*G[i][i]*G[j][j];
427 }
428 }
429 }
430 v_re = temp_PDF*B[0]*B[1];
431 return v_re;
432}
433
434
435void EvtDsToEtappipi0::calEva(double* EtaP, double* Pi, double* Pi0, double *mass1, double *width1, double *amp, double *phase,int* g0,int* spin, int* modetype, int nstates, double & Result)
436{
437 double P12[4], P23[4], P13[4];
438 double cof[2], amp_PDF[2], PDF[2];
439 double spi, spi0, setap;
440 double s12, s13, s23;
441 for(int i=0; i<4; i++){
442 P23[i] = Pi[i] + Pi0[i];
443 P12[i] = EtaP[i] + Pi[i];
444 P13[i] = EtaP[i] + Pi0[i];
445 }
446 setap = SCADot(EtaP,EtaP);
447 spi = SCADot(Pi,Pi);
448 spi0 = SCADot(Pi0,Pi0);
449 s12 = SCADot(P12,P12);
450 s13 = SCADot(P13,P13);
451 s23 = SCADot(P23,P23);
452
453 double pro[2], temp_PDF, amp_tmp[2];
454 double mass1sq;
455 amp_PDF[0] = 0;
456 amp_PDF[1] = 0;
457 PDF[0] = 0;
458 PDF[1] = 0;
459 amp_tmp[0] = 0;
460 amp_tmp[1] = 0;
461 for(int i=0; i<nstates; i++) {
462 amp_tmp[0] = 0;
463 amp_tmp[1] = 0;
464 mass1sq = mass1[i]*mass1[i];
465 cof[0] = amp[i]*cos(phase[i]);
466 cof[1] = amp[i]*sin(phase[i]);
467 temp_PDF = 0;
468
469 if(modetype[i] == 23){
470 temp_PDF = DDalitz(Pi, Pi0, EtaP, spin[i], mass1[i]);
471 if(g0[i]==1) propagatorGS(mass1sq,mass1[i],width1[i],s23,spi,spi0,9.0,pro);
472 if(g0[i]==12) PiPiSWASS(s23,spi,spi0,pro);
473 if(g0[i]==0){
474 pro[0] = 1;
475 pro[1] = 0;
476 }
477 // printf("cof1 = %.15f \n",cof[1]);
478
479 amp_tmp[0] = temp_PDF*pro[0];
480 amp_tmp[1] = temp_PDF*pro[1];
481
482 }
483 Com_Multi(amp_tmp,cof,amp_PDF);
484 PDF[0] += amp_PDF[0];
485 PDF[1] += amp_PDF[1];
486 }
487 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
488 //printf("value = %.15f\n",value);
489 Result = value;
490}
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")
a0
character *LEPTONflag integer iresonances real zeta5 real a0
Definition Ekhara/Ekhara-00-01-02/src/ekhara/common.h:4
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
EvtDsToEtappipi0::clone
EvtDecayBase * clone()
Definition EvtDsToEtappipi0.cc:38
EvtDsToEtappipi0::~EvtDsToEtappipi0
virtual ~EvtDsToEtappipi0()
Definition EvtDsToEtappipi0.cc:32
EvtDsToEtappipi0::decay
void decay(EvtParticle *p)
Definition EvtDsToEtappipi0.cc:109
EvtDsToEtappipi0::getName
void getName(std::string &name)
Definition EvtDsToEtappipi0.cc:34
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