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