BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToomegaenu.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// Environment:
3// This software is part of models developed at BES collaboration
4// based on the EvtGen framework. If you use all or part
5// of it, please give an appropriate acknowledgement.
6//
7// Copyright Information: See EvtGen/BesCopyright
8// Copyright (A) 2006 Ping Rong-Gang @IHEP
9//
10// Module: EvtDsToomegaenu.cc
11// the necessary file: EvtDsToomegaenu.hh
12//
13// Description: Ds+ -> omega e+ nu, omega -> pi+ pi- pi0
14//
15// Modification history:
16// Liaoyuan Dong May 29 00:41:03 2023 Module created
17//------------------------------------------------------------------------
18#include "EvtGenBase/EvtPatches.hh"
19#include "EvtGenBase/EvtParticle.hh"
20#include "EvtGenBase/EvtGenKine.hh"
21#include "EvtGenBase/EvtPDL.hh"
22#include "EvtGenBase/EvtReport.hh"
23#include "EvtGenModels/EvtDsToomegaenu.hh"
24#include "EvtGenBase/EvtComplex.hh"
25#include "EvtGenBase/EvtDecayTable.hh"
26#include "EvtGenBase/EvtVector3R.hh"
27#include <stdlib.h>
28
29EvtDsToomegaenu::~EvtDsToomegaenu() {}
30
31void EvtDsToomegaenu::getName(std::string& model_name){
32 model_name="DsToomegaenu";
33}
34
38
39void EvtDsToomegaenu::init(){
40 checkNArg(0);
41 checkNDaug(5); // Ds+ -> pi+ pi- pi0 e+ nu
42 // 0 --> 1 2 3 4 5
43 checkSpinParent(EvtSpinType::SCALAR);
44 checkSpinDaughter(0,EvtSpinType::SCALAR);
45 checkSpinDaughter(1,EvtSpinType::SCALAR);
46 checkSpinDaughter(2,EvtSpinType::SCALAR);
47
48 std::cout << "Initializing EvtDsToomegaenu" << std::endl;
49 mV = 1.81;
50 mA = 2.61;
51 V_0 = 1.411;
52 A1_0 = 1;
53 A2_0 = 0.788;
54 mV2 = mV*mV;
55 mA2 = mA*mA;
56
57 m_omega = 0.78266; // omega
58 w_omega = 0.00868;
59 m_rho = 0.77526; // rho
60 w_rho = 0.14910;
61 rBW = 3.0;
62 rBW2 = rBW*rBW;
63 mw_omega = m_omega*w_omega;
64 m2_omega = m_omega*m_omega;
65 m2_rho = m_rho*m_rho;
66
67 mDs = 1.9683;
68 m2Ds = mDs*mDs;
69 mpi = 0.13957;
70 m2pi = mpi*mpi;
71 mpi0 = 0.1349766;
72 m2pi0 = mpi0*mpi0;
73
74 Pi = atan2(0.0,-1.0);
75 root2 = sqrt(2.);
76 root2d3 = sqrt(2./3);
77 root1d2 = sqrt(0.5);
78 root3d2 = sqrt(1.5);
79}
80
81void EvtDsToomegaenu::initProbMax() {
82 setProbMax(228000.0);
83}
84
85void EvtDsToomegaenu::decay( EvtParticle *p){
86/*
87 double maxprob = 0.0;
88 for(int ir=0;ir<=60000000;ir++){
89 p->initializePhaseSpace(getNDaug(),getDaugs());
90 EvtVector4R _pip = p->getDaug(0)->getP4(); // pi+
91 EvtVector4R _pim = p->getDaug(1)->getP4(); // pi-
92 EvtVector4R _pi0 = p->getDaug(2)->getP4(); // pi0
93 EvtVector4R _e = p->getDaug(3)->getP4(); // e+
94 EvtVector4R _nu = p->getDaug(4)->getP4(); // nu
95 int pid = EvtPDL::getStdHep(p->getDaug(3)->getId());
96 int charm;
97 if(pid == -11) charm = 1;
98 else charm = -1;
99 double m2, q2, cosV, cosL, chi, s12, s13, s23, spin_V;
100 //EvtVector4R pip(0.43708151134234324298,-0.0063461065118795333476,0.411859161170534116,0.043499647450492222311);
101 //EvtVector4R pim(0.25237575453703847694,-0.078585671855406255548,-0.19138558256021315218,-0.03754444974755244413);
102 //EvtVector4R pi0(0.50598020209955174575,0.37013555928038383014,0.14129868575000212316,-0.28430903904443316499);
103 //EvtVector4R e(0.33268353968172292845,-0.18935977852196397841,0.067639189874503583,0.26503941085708321301);
104 //EvtVector4R nu(0.44017899233934343339,-0.095844002391134053287,-0.42941145423482673937,0.013314430484410192876);
105
106 KinVGen(_pip, _pim, _pi0, _e, _nu, charm, m2, q2, cosV, cosL, chi, s12, s13, s23, spin_V);
107 double _prob = calPDF(m2, q2, cosV, cosL, chi, s12, s13, s23, spin_V);
108 if(_prob>=maxprob) {
109 maxprob=_prob;
110 std::cout << ir << " prob= " << _prob << std::endl;
111 cout << "EvtVector4R pip" << _pip << ";"<< endl;
112 cout << "EvtVector4R pim" << _pim << ";"<< endl;
113 cout << "EvtVector4R pi0" << _pi0 << ";"<< endl;
114 cout << "EvtVector4R e" << _e << ";"<< endl;
115 cout << "EvtVector4R nu" << _nu << ";"<< endl;
116 }
117 }
118 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
119*/
120 p->initializePhaseSpace(getNDaug(),getDaugs());
121 EvtVector4R pip = p->getDaug(0)->getP4();
122 EvtVector4R pim = p->getDaug(1)->getP4();
123 EvtVector4R pi0 = p->getDaug(2)->getP4();
124 EvtVector4R e = p->getDaug(3)->getP4();
125 EvtVector4R nu = p->getDaug(4)->getP4();
126
127 int pid = EvtPDL::getStdHep(p->getDaug(3)->getId());
128 int charm;
129 if(pid == -11) charm = 1;
130 else charm = -1;
131 double m2, q2, cosV, cosL, chi, s12, s13, s23, spin_V;
132 KinVGen(pip, pim, pi0, e, nu, charm, m2, q2, cosV, cosL, chi, s12, s13, s23, spin_V);
133 double prob = calPDF(m2, q2, cosV, cosL, chi, s12, s13, s23, spin_V);
134 setProb(prob);
135 return;
136}
137
138void EvtDsToomegaenu::KinVGen(EvtVector4R vp4_pip, EvtVector4R vp4_pim, EvtVector4R vp4_pi0, EvtVector4R vp4_Lep, EvtVector4R vp4_Nu, int charm, double& m2, double& q2, double& cosV, double& cosL, double& chi, double& s12, double& s13, double& s23, double& spin_V)
139{
140 EvtVector4R vp4_3pi = vp4_pip + vp4_pim + vp4_pi0;
141 EvtVector4R vp4_rho0 = vp4_pip + vp4_pim;
142 EvtVector4R vp4_rhop = vp4_pip + vp4_pi0;
143 EvtVector4R vp4_rhom = vp4_pim + vp4_pi0;
144 EvtVector4R vp4_W = vp4_Lep + vp4_Nu;
145
146 m2 = vp4_3pi.mass2();
147 q2 = vp4_W.mass2();
148 s12 = vp4_rho0.mass2();
149 s13 = vp4_rhop.mass2();
150 s23 = vp4_rhom.mass2();
151
152 EvtVector4R boost;
153 boost.set(vp4_3pi.get(0), -vp4_3pi.get(1), -vp4_3pi.get(2), -vp4_3pi.get(3));
154 EvtVector4R vp4_rhob = boostTo(vp4_rho0, boost);
155 cosV = vp4_rhob.dot(vp4_3pi)/(vp4_rhob.d3mag()*vp4_3pi.d3mag());
156
157 EvtVector4R vp4_pipb = boostTo(vp4_pip, boost);
158 EvtVector4R vp4_pimb = boostTo(vp4_pim, boost);
159 EvtVector4R R = vp4_pipb.cross(vp4_pimb);
160 spin_V = R.d3mag(); // sart(|pi+ x pi-|^2)
161
162 boost.set(vp4_W.get(0), -vp4_W.get(1), -vp4_W.get(2), -vp4_W.get(3));
163 EvtVector4R vp4_Lepp = boostTo(vp4_Lep, boost);
164 cosL = vp4_Lepp.dot(vp4_W)/(vp4_Lepp.d3mag()*vp4_W.d3mag());
165
166 EvtVector4R V = vp4_3pi/vp4_3pi.d3mag();
167 EvtVector4R C = vp4_rhob.cross(V);
168 C/=C.d3mag();
169 EvtVector4R D = vp4_Lepp.cross(V);
170 D/=D.d3mag();
171 double sinx = C.cross(V).dot(D);
172 double cosx = C.dot(D);
173 chi = sinx>0? acos(cosx): -acos(cosx);
174 if(charm==-1) chi=-chi;
175}
176
177double EvtDsToomegaenu::calPDF(double m2, double q2, double cosV, double cosL, double chi, double s12, double s13, double s23, double spin_V) {
178 double m = sqrt(m2);
179 double q = sqrt(q2);
180 double m12 = sqrt(s12); // m(pi+pi-)
181 double m13 = sqrt(s13); // m(pi+pi0)
182 double m23 = sqrt(s23); // m(pi-pi0)
183
184 EvtComplex A_rhopi(0.0,0.0);
185 ResonanceV(s12, m12, s13, m13, s23, m23, A_rhopi);
186 double A2_phi = spin_V*spin_V*abs2(A_rhopi);
187
188 EvtComplex f11(0.0,0.0);
189 EvtComplex f21(0.0,0.0);
190 EvtComplex f31(0.0,0.0);
191 ResonanceP(m2, m, q2, q, s12, m12, s13, m13, s23, m23, f11, f21, f31);
192
193 //begin to calculate pdf value
194 double cosV2 = cosV*cosV;
195 double sinV = sqrt(1.0-cosV2);
196 double sinV2 = sinV*sinV;
197
198 EvtComplex F1 = f11*cosV;
199 EvtComplex F2 = f21*root1d2;
200 EvtComplex F3 = f31*root1d2;
201 //cout << "F1= " << F1 << abs2(F1) << endl;
202 //cout << "F2= " << F2 << abs2(F2) << endl;
203 //cout << "F3= " << F3 << abs2(F3) << endl;
204
205 double I1 = 0.25*( abs2(F1) + 1.5*sinV2*(abs2(F2) + abs2(F3) ) );
206 double I2 = -0.25*( abs2(F1) - 0.5*sinV2*(abs2(F2) + abs2(F3) ) );
207 double I3 = -0.25*( abs2(F2) - abs2(F3) )*sinV2;
208 double I4 = real( conj(F1)*F2 )*sinV*0.5;
209 double I5 = real( conj(F1)*F3 )*sinV;
210 double I6 = real( conj(F2)*F3 )*sinV2;
211 //cout << "Ix = " << I1 << ", " << I2 << ", " << I3 << ", " << I4 << ", " << I5 << ", " << I6 << endl;
212
213 double coschi = cos(chi);
214 double sinchi = sin(chi);
215 double sin2chi = 2.0*sinchi*coschi;
216 double cos2chi = 1.0 - 2.0*sinchi*sinchi;
217
218 double sinL = sqrt(1. - cosL*cosL);
219 double sinL2 = sinL*sinL;
220 double sin2L = 2.0*sinL*cosL;
221 double cos2L = 1.0 - 2.0*sinL2;
222
223 double I = A2_phi*(I1 + I2*cos2L + I3*sinL2*cos2chi + I4*sin2L*coschi + I5*sinL*coschi + I6*cosL);
224 return I;
225}
226
227void EvtDsToomegaenu::ResonanceV(double s12, double m12, double s13, double m13, double s23, double m23, EvtComplex& A_rhopi)
228{
229 double G_rho0 = w_rho*m2_rho/s12*pow((s12-4.0*m2pi0)/(m2_rho-4.0*m2pi0), 1.5); // Gamma(rho0) Eq.(9)
230 double G_rhop = w_rho*m2_rho/s13*pow((s13-4.0*m2pi)/( m2_rho-4.0*m2pi), 1.5); // Gamma(rho+) Eq.(9)
231 double G_rhom = w_rho*m2_rho/s23*pow((s23-4.0*m2pi)/( m2_rho-4.0*m2pi), 1.5); // Gamma(rho-) Eq.(9)
232
233 double AA0 = s12/m2_rho - 1.0;
234 double BB0 = m12*G_rho0/m2_rho;
235 //cout << "AA0 = " << AA0 << " BB0 = " << BB0 << endl;
236
237 EvtComplex amp0 = EvtComplex(1.0, 0.0)/EvtComplex(AA0,BB0);
238
239 double AAp = s13/m2_rho - 1.0;
240 double BBp = m13*G_rhop/m2_rho;
241 //cout << "AAp = " << AAp << " BBp = " << BBp << endl;
242
243 EvtComplex ampp = EvtComplex(1.0, 0.0)/EvtComplex(AAp,BBp);
244
245 double AAm = s23/m2_rho - 1.0;
246 double BBm = m23*G_rhom/m2_rho;
247 //cout << "AAm = " << AAm << " BBm = " << BBm << endl;
248
249 EvtComplex ampm = EvtComplex(1.0, 0.0)/EvtComplex(AAm,BBm);
250
251 A_rhopi = amp0 + ampp + ampm; // Eq.(7)
252}
253
254void EvtDsToomegaenu::ResonanceP(double m2, double m, double q2, double q, double s12, double m12, double s13, double m13, double s23, double m23, EvtComplex& F11, EvtComplex& F21, EvtComplex& F31)
255{
256 double p_phi = getPStar(m2Ds, m2, q2);
257 double summDsm = mDs+m;
258 double V = V_0 /(1.0 - q2/(mV2));
259 double A1 = A1_0/(1.0 - q2/(mA2));
260 double A2 = A2_0/(1.0 - q2/(mA2));
261 double A = summDsm*A1;
262 double B = 2.0*mDs*p_phi/summDsm*V;
263
264 //construct the helicity form factor
265 double H0 = 0.5/(m*q)*((m2Ds-m2-q2)*summDsm*A1-4.0*(m2Ds*p_phi*p_phi)/summDsm*A2);
266 double Hp = A-B;
267 double Hm = A+B;
268 //cout << "H0 = " << H0 << " Hp = " << Hp << " Hm = " << Hm << endl;
269
270 //calculate alpha
271 //double pStar0 = getPStar(m2_omega, s12, m2pi0);
272 double pStar0 = getPStar(m2_omega, m2_rho, m2pi0);
273 double alpha = sqrt(3.*Pi*0.892/(pStar0*w_omega));
274 //cout << "alpha = " << alpha << endl;
275
276 double AA = m2_omega-m2;
277 //cout << "AA = " << AA << endl;
278
279 //construct BW for omega -> rho0 pi0
280 //double F0 = getF1(m2, m, s12, m2pi0);
281 double F0 = getF1(m2, m, m2_rho, m2pi0);
282 //cout << "F0 = " << F0 << endl;
283 //double width0 = getWidth1(m2, m, s12, m2pi0, w_omega);
284 double width0 = getWidth1(m2, m, m2_rho, m2pi0, w_omega);
285 //cout << "width0 = " << width0 << endl;
286
287 EvtComplex C0(mw_omega*F0,0.0);
288 //cout << "C0 = " << C0 << endl;
289 double BB0 = -m_omega*width0;
290 //cout << "BB0 = " << BB0 << endl;
291 EvtComplex amp0 = C0/EvtComplex(AA,BB0);
292 //cout << "amp0 = " << amp0 << endl;
293
294 //construct BW for omega -> rho+ pi-
295 //double Fp = getF1(m2, m, s13, m2pi);
296 double Fp = getF1(m2, m, m2_rho, m2pi);
297 //cout << "Fp = " << Fp << endl;
298 //double widthp = getWidth1(m2, m, s13, m2pi, w_omega);
299 double widthp = getWidth1(m2, m, m2_rho, m2pi, w_omega);
300 //cout << "widthp = " << widthp << endl;
301
302 EvtComplex Cp(mw_omega*Fp,0.0);
303 //cout << "Cp = " << Cp << endl;
304 double BBp = -m_omega*widthp;
305 //cout << "BBp = " << BBp << endl;
306 EvtComplex ampp = Cp/EvtComplex(AA,BBp);
307 //cout << "ampp = " << ampp << endl;
308
309 //construct BW for omega -> rho- pi+
310 //double Fm = getF1(m2, m, s23, m2pi);
311 double Fm = getF1(m2, m, m2_rho, m2pi);
312 //cout << "Fm = " << Fm << endl;
313 //double widthm = getWidth1(m2, m, s23, m2pi, w_omega);
314 double widthm = getWidth1(m2, m, m2_rho, m2pi, w_omega);
315 //cout << "widthm = " << widthm << endl;
316
317 EvtComplex Cm(mw_omega*Fm,0.0);
318 //cout << "Cm = " << Cm << endl;
319 double BBm = -m_omega*widthm;
320 //cout << "BBm = " << BBm << endl;
321 EvtComplex ampm = Cm/EvtComplex(AA,BBm);
322 //cout << "ampm = " << ampm << endl;
323
324 EvtComplex amp = amp0 + ampp + ampm; // Eq. (15)
325 //cout << "amp = " << amp << endl;
326
327 double alpham2 = alpha*2.0;
328 F11 = amp*alpham2*q*H0*root2;
329 F21 = amp*alpham2*q*(Hp+Hm);
330 F31 = amp*alpham2*q*(Hp-Hm);
331}
332
333double EvtDsToomegaenu::getF1(double sa, double ma, double sb, double sc) // a -> b + c
334{
335 double pStar = getPStar(sa, sb, sc);
336 double pStar0 = getPStar(m2_omega, sb, sc);
337 double B = 1./sqrt(1.+rBW2*pStar*pStar);
338 double B0 = 1./sqrt(1.+rBW2*pStar0*pStar0);
339 double F = pStar/pStar0*B/B0;
340 return F;
341}
342
343double EvtDsToomegaenu::getWidth1(double sa, double ma, double sb, double sc, double width0)
344{
345 double pStar = getPStar(sa, sb, sc);
346 double pStar0 = getPStar(m2_omega, sb, sc);
347 double F = getF1(sa, ma, sb, sc);
348 double width = width0*pStar/pStar0*m_omega/ma*F*F; // Eq. (17)
349 return width;
350}
351
352double EvtDsToomegaenu::getPStar(double s, double s1, double s2)
353{
354 double x = s + s1 - s2;
355 double t = 0.25*x*x/s - s1;
356 double p;
357 if (t>0.0) {
358 p = sqrt(t);
359 } else {
360 //std::cout << " Hello, pstar is less than 0.0" << std::endl;
361 p = sqrt(-t);
362 //p = 0.02;
363 }
364 return p;
365}
sin
double sin(const BesAngle a)
Definition BesAngle.h:210
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
x
Double_t x[10]
Definition DataBase/tau_mode.c:57
I
const DifComplex I
conj
Evt3Rank3C conj(const Evt3Rank3C &t2)
Definition Evt3Rank3C.cc:175
abs2
double abs2(const EvtComplex &c)
Definition EvtComplex.hh:221
boostTo
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
Definition EvtDiracSpinor.cc:75
alpha
const double alpha
Definition FastVertexFit.cxx:4
s
XmlRpcServer s
Definition HelloServer.cpp:11
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
C
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations C
Definition RRes.h:29
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
EvtDsToomegaenu::clone
EvtDecayBase * clone()
Definition EvtDsToomegaenu.cc:35
EvtDsToomegaenu::~EvtDsToomegaenu
virtual ~EvtDsToomegaenu()
Definition EvtDsToomegaenu.cc:29
EvtDsToomegaenu::decay
void decay(EvtParticle *p)
Definition EvtDsToomegaenu.cc:85
EvtDsToomegaenu::getName
void getName(std::string &name)
Definition EvtDsToomegaenu.cc:31
EvtPDL::getStdHep
static int getStdHep(EvtId id)
Definition EvtPDL.hh:56
EvtParticle::getId
EvtId getId() const
Definition EvtParticle.cc:112
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::dot
double dot(const EvtVector4R &v2) const
Definition EvtVector4R.cc:199
EvtVector4R::get
double get(int i) const
Definition EvtVector4R.hh:179
EvtVector4R::cross
EvtVector4R cross(const EvtVector4R &v2)
Definition EvtVector4R.cc:171
EvtVector4R::d3mag
double d3mag() const
Definition EvtVector4R.cc:186
EvtVector4R::mass2
double mass2() const
Definition EvtVector4R.hh:116
EvtVector4R::set
void set(int i, double d)
Definition EvtVector4R.hh:183
s12
double double double double * s12
Definition qcdloop1.h:77
s23
double double double double double * s23
Definition qcdloop1.h:77
m2
double double * m2
Definition qcdloop1.h:75