BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDToKSKpi0.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2//
3// Environment:
4// This software is part of models developed at BES collaboration
5// based on the EvtGen framework. If you use all or part
6// of it, please give an appropriate acknowledgement.
7//
8// Copyright Information: See EvtGen/BesCopyright
9// Copyright (A) 2006 Ping Rong-Gang @IHEP
10//
11// Module: EvtDToKSKpi0.cc
12//
13// Description: Phys. Rev. D 104, 012006 (2021)
14//
15// Modification history:
16//
17// Liaoyuan Dong Aug 11, 2022 Module created
18//
19//------------------------------------------------------------------------
20#include "EvtDToKSKpi0.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/EvtVector4R.hh"
27#include <stdlib.h>
28#include <string>
29using namespace std; //::endl;
30
31EvtDToKSKpi0::~EvtDToKSKpi0() {}
32
33void EvtDToKSKpi0::getName( std::string& model_name ) { model_name = "DToKSKpi0"; }
34
35EvtDecayBase* EvtDToKSKpi0::clone() { return new EvtDToKSKpi0; }
36
37void EvtDToKSKpi0::init() {
38 checkNArg( 0 );
39 checkNDaug( 3 );
40 checkSpinParent( EvtSpinType::SCALAR );
41 // std::cout << "Initializing EvtDToKSKpi0" << std::endl;
42
43 _nd = 3;
44 _mDp = 1.86965;
45 c_motherMass = _mDp;
46 _mDp2 = _mDp * _mDp;
47 _mDp2inv = 1. / _mDp2;
48 KsMass = 0.497611;
49 KpMass = 0.493677;
50 pi0Mass = 0.134977;
51 etamass = 0.547862;
52 pipMass = 0.13957061;
53 c_meson_radius_inter = 1.5;
54 c_meson_radius_Dp = 5;
55}
56void EvtDToKSKpi0::initProbMax() { setProbMax( 1343.2 ); }
57
58void EvtDToKSKpi0::decay( EvtParticle* p ) {
59 /*
60 double maxprob = 0.0;
61 for(int ir=0;ir<=60000000;ir++){
62 p->initializePhaseSpace(getNDaug(),getDaugs());
63 for(int i=0; i<_nd; i++){
64 _p4Lab[i]=p->getDaug(i)->getP4Lab();
65 _p4CM[i]=p->getDaug(i)->getP4();
66 }
67 double Prob = AmplitudeSquare();
68 if(Prob>maxprob) {
69 maxprob=Prob;
70 std::cout << "Max PDF = " << ir << " prob= " << Prob << std::endl;
71 }
72 }
73 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
74 */
75
76 p->initializePhaseSpace( getNDaug(), getDaugs() );
77 for ( int i = 0; i < _nd; i++ )
78 {
79 _p4Lab[i] = p->getDaug( i )->getP4Lab();
80 _p4CM[i] = p->getDaug( i )->getP4();
81 }
82 double prob = AmplitudeSquare();
83 setProb( prob );
84 return;
85}
86
87double EvtDToKSKpi0::twoBodyCMmom( double rMassSq, double d1m, double d2m ) {
88 double kin1 = 1 - pow( d1m + d2m, 2 ) / rMassSq;
89 kin1 = kin1 >= 0 ? sqrt( kin1 ) : 1;
90 double kin2 = 1 - pow( d1m - d2m, 2 ) / rMassSq;
91 kin2 = kin2 >= 0 ? sqrt( kin2 ) : 1;
92
93 double ret = 0.5 * sqrt( rMassSq ) * kin1 * kin2;
94 return ret;
95}
96
97double EvtDToKSKpi0::dampingFactorSquare( const double& cmmom, const int& spin,
98 const double& mRadius ) {
99 double square = mRadius * mRadius * cmmom * cmmom;
100 double dfsq = 1 + square; // This accounts for spin 1
101 // if (2 == spin) dfsq += 8 + 2*square + square*square; // Coefficients are 9, 3, 1.
102 double dfsqres = dfsq + 8 + 2 * square + square * square;
103
104 // Spin 3 and up not accounted for.
105 double ret = ( spin == 2 ) ? dfsqres : dfsq;
106 return ret;
107}
108
109// here, m12 is the resonance
110double EvtDToKSKpi0::spinFactor( int spin, double motherMass, double daug1Mass,
111 double daug2Mass, double daug3Mass, double m12, double m13,
112 double m23 ) {
113 if ( spin == 0 ) return 1;
114
115 double _mA = daug1Mass;
116 double _mB = daug2Mass;
117 double _mC = daug3Mass;
118 double _mAB = m12;
119 double _mAC = m13;
120 double _mBC = m23;
121
122 double massFactor = 1.0 / _mAB;
123 double sFactor = -1;
124 sFactor *= ( ( _mBC - _mAC ) + ( massFactor * ( motherMass * motherMass - _mC * _mC ) *
125 ( _mA * _mA - _mB * _mB ) ) );
126
127 if ( spin == 2 )
128 {
129 sFactor *= sFactor;
130 double extraterm = ( ( _mAB - ( 2 * motherMass * motherMass ) - ( 2 * _mC * _mC ) ) +
131 massFactor * pow( motherMass * motherMass - _mC * _mC, 2 ) );
132 extraterm *= ( ( _mAB - ( 2 * _mA * _mA ) - ( 2 * _mB * _mB ) ) +
133 massFactor * pow( _mA * _mA - _mB * _mB, 2 ) );
134 extraterm /= 3;
135 sFactor -= extraterm;
136 }
137
138 return sFactor;
139}
140
141// here, id == 1 -> pi0 K+ resonance, 2 -> pi0 Ks resonance, 3 -> K+ Ks resonance
142EvtComplex EvtDToKSKpi0::RBW( int id, double resmass, double reswidth, int spin ) {
143 double resmass2 = pow( resmass, 2 );
144
145 EvtVector4R p1, p2, p3;
146 double mass_daug1, mass_daug2, mass_daug3;
147 if ( id == 1 )
148 {
149 p1 = _pd[0];
150 mass_daug1 = pi0Mass;
151 p2 = _pd[1];
152 mass_daug2 = KpMass;
153 p3 = _pd[2];
154 mass_daug3 = KsMass;
155 }
156 if ( id == 2 )
157 {
158 p1 = _pd[2];
159 mass_daug1 = KsMass;
160 p2 = _pd[0];
161 mass_daug2 = pi0Mass;
162 p3 = _pd[1];
163 mass_daug3 = KpMass;
164 }
165 if ( id == 3 )
166 {
167 p1 = _pd[1];
168 mass_daug1 = KpMass;
169 p2 = _pd[2];
170 mass_daug2 = KsMass;
171 p3 = _pd[0];
172 mass_daug3 = pi0Mass;
173 }
174
175 double rMassSq = ( p1 + p2 ).mass2();
176 double m12 = ( p1 + p2 ).mass2();
177 double m13 = ( p1 + p3 ).mass2();
178 double m23 = ( p2 + p3 ).mass2();
179
180 double rMass = sqrt( rMassSq );
181 double frFactor = 1;
182 double fdFactor = 1;
183
184 // Calculate momentum of the two daughters in the resonance rest frame
185 double measureDaughterMoms = twoBodyCMmom( rMassSq, mass_daug1, mass_daug2 );
186 double nominalDaughterMoms = twoBodyCMmom( resmass2, mass_daug1, mass_daug2 );
187
188 if ( spin != 0 )
189 {
190 frFactor = dampingFactorSquare( nominalDaughterMoms, spin, c_meson_radius_inter ) /
191 dampingFactorSquare( measureDaughterMoms, spin, c_meson_radius_inter );
192
193 double measureDaughterMoms2 =
194 twoBodyCMmom( c_motherMass * c_motherMass, rMass, mass_daug3 );
195 double nominalDaughterMoms2 =
196 twoBodyCMmom( c_motherMass * c_motherMass, resmass, mass_daug3 );
197 fdFactor = dampingFactorSquare( nominalDaughterMoms2, spin, c_meson_radius_Dp ) /
198 dampingFactorSquare( measureDaughterMoms2, spin, c_meson_radius_Dp );
199 }
200 double A = ( resmass2 - rMassSq );
201 double B = resmass2 * reswidth *
202 pow( measureDaughterMoms / nominalDaughterMoms, 2.0 * spin + 1 ) * frFactor /
203 sqrt( rMassSq );
204 double C = 1.0 / ( pow( A, 2 ) + pow( B, 2 ) );
205
206 EvtComplex ret( A * C, B * C );
207 ret *= sqrt( frFactor * fdFactor );
208 ret *= spinFactor( spin, c_motherMass, mass_daug1, mass_daug2, mass_daug3, m12, m13, m23 );
209
210 // if(id == 2)
211 // cout << "m13 resonance spinFactor = " << spinFactor(spin, c_motherMass, mass_daug1,
212 // mass_daug2, mass_daug3, m12, m13, m23) << endl; if(id == 1) cout << "m12 resonance
213 // spinFactor = " << spinFactor(spin, c_motherMass, mass_daug1, mass_daug2, mass_daug3, m12,
214 // m13, m23) << endl;
215
216 return ret;
217}
218
219EvtComplex EvtDToKSKpi0::Flatte( int id, double resmass, double g1, double rg2og1 ) {
220 double resmass2 = pow( resmass, 2 );
221
222 EvtVector4R p1, p2, p3;
223 double mass_daug1, mass_daug2, mass_daug3;
224 if ( id == 1 )
225 {
226 p1 = _pd[0];
227 mass_daug1 = pi0Mass;
228 p2 = _pd[1];
229 mass_daug2 = KpMass;
230 p3 = _pd[2];
231 mass_daug3 = KsMass;
232 }
233 if ( id == 2 )
234 {
235 p1 = _pd[2];
236 mass_daug1 = KsMass;
237 p2 = _pd[0];
238 mass_daug2 = pi0Mass;
239 p3 = _pd[1];
240 mass_daug3 = KpMass;
241 }
242 if ( id == 3 )
243 {
244 p1 = _pd[1];
245 mass_daug1 = KpMass;
246 p2 = _pd[2];
247 mass_daug2 = KsMass;
248 p3 = _pd[0];
249 mass_daug3 = pi0Mass;
250 }
251
252 double rMassSq = ( p1 + p2 ).mass2();
253 double m12 = ( p1 + p2 ).mass2();
254 double m13 = ( p1 + p3 ).mass2();
255 double m23 = ( p2 + p3 ).mass2();
256 double rMass = sqrt( rMassSq );
257
258 double s = m12;
259
260 double rhoetapi = 2 * twoBodyCMmom( s, KsMass, KpMass ) / sqrt( s );
261 double rhoKKbar = 2 * twoBodyCMmom( s, etamass, pipMass ) / sqrt( s );
262 double img = rhoetapi * g1 + rhoKKbar * g1 * rg2og1;
263
264 EvtComplex ret = EvtComplex( 1, 0 ) / EvtComplex( resmass * resmass - s, -img );
265 return ret;
266}
267
268// here, id == 1 -> pi0 K+ resonance, 2 -> pi0 Ks resonance, 3 -> K+ Ks resonance
269EvtComplex EvtDToKSKpi0::LASS( int id, double resmass, double reswidth ) {
270 int spin = 0;
271 double resmass2 = pow( resmass, 2 );
272
273 EvtVector4R p1, p2, p3;
274 double mass_daug1, mass_daug2, mass_daug3;
275 if ( id == 1 )
276 {
277 p1 = _pd[0];
278 mass_daug1 = pi0Mass;
279 p2 = _pd[1];
280 mass_daug2 = KpMass;
281 p3 = _pd[2];
282 mass_daug3 = KsMass;
283 }
284 if ( id == 2 )
285 {
286 p1 = _pd[2];
287 mass_daug1 = KsMass;
288 p2 = _pd[0];
289 mass_daug2 = pi0Mass;
290 p3 = _pd[1];
291 mass_daug3 = KpMass;
292 }
293 if ( id == 3 )
294 {
295 p1 = _pd[1];
296 mass_daug1 = KpMass;
297 p2 = _pd[2];
298 mass_daug2 = KsMass;
299 p3 = _pd[0];
300 mass_daug3 = pi0Mass;
301 }
302
303 double rMassSq = ( p1 + p2 ).mass2();
304 double m12 = ( p1 + p2 ).mass2();
305 double m13 = ( p1 + p3 ).mass2();
306 double m23 = ( p2 + p3 ).mass2();
307
308 double rMass = sqrt( rMassSq );
309 double frFactor = 1;
310 double fdFactor = 1;
311
312 // Calculate momentum of the two daughters in the resonance rest frame
313 double measureDaughterMoms = twoBodyCMmom( rMassSq, mass_daug1, mass_daug2 );
314 double nominalDaughterMoms = twoBodyCMmom( resmass2, mass_daug1, mass_daug2 );
315 if ( spin != 0 )
316 {
317 frFactor = dampingFactorSquare( nominalDaughterMoms, spin, c_meson_radius_inter ) /
318 dampingFactorSquare( measureDaughterMoms, spin, c_meson_radius_inter );
319 double measureDaughterMoms2 =
320 twoBodyCMmom( c_motherMass * c_motherMass, rMass, mass_daug3 );
321 double nominalDaughterMoms2 =
322 twoBodyCMmom( c_motherMass * c_motherMass, resmass, mass_daug3 );
323 fdFactor = dampingFactorSquare( nominalDaughterMoms2, spin, c_meson_radius_Dp ) /
324 dampingFactorSquare( measureDaughterMoms2, spin, c_meson_radius_Dp );
325 }
326
327 double q = measureDaughterMoms;
328 double g = reswidth * pow( measureDaughterMoms / nominalDaughterMoms, 2.0 * spin + 1 ) *
329 frFactor / sqrt( rMassSq );
330 g *= resmass;
331
332 const double _a = 0.113;
333 const double _r = -33.8;
334 const double _R = 1;
335 const double _phiR = -109.7 * 3.141592653 / 180.0;
336 const double _B = 0.96;
337 const double _phiB = 0.1 * 3.141592653 / 180.0;
338
339 // background phase motion
340 double cot_deltaB = ( 1.0 / ( _a * q ) ) + 0.5 * _r * q;
341 double qcot_deltaB = ( 1.0 / _a ) + 0.5 * _r * q * q;
342
343 // calculate resonant part
344 EvtComplex expi2deltaB = EvtComplex( qcot_deltaB, q ) / EvtComplex( qcot_deltaB, -q );
345 EvtComplex resT = EvtComplex( cos( _phiR + 2 * _phiB ), sin( _phiR + 2 * _phiB ) ) * _R;
346
347 EvtComplex prop = EvtComplex( 1, 0 ) / EvtComplex( resmass2 - rMassSq, -(resmass)*g );
348 resT *= prop * ( resmass2 * reswidth / nominalDaughterMoms ) * expi2deltaB;
349
350 // calculate bkg part
351 resT += EvtComplex( cos( _phiB ), sin( _phiB ) ) * _B *
352 ( cos( _phiB ) + cot_deltaB * sin( _phiB ) ) * sqrt( rMassSq ) /
353 EvtComplex( qcot_deltaB, -q );
354
355 resT *= sqrt( frFactor * fdFactor );
356 resT *= spinFactor( spin, c_motherMass, mass_daug1, mass_daug2, mass_daug3, m12, m13, m23 );
357
358 return resT;
359}
360
361double EvtDToKSKpi0::AmplitudeSquare() {
362 EvtVector4R dp1 = GetDaugMomLab( 0 ); // KS
363 EvtVector4R dp2 = GetDaugMomLab( 1 ); // K+
364 EvtVector4R dp3 = GetDaugMomLab( 2 ); // pi0
365 _pd[0] = dp3;
366 _pd[1] = dp2;
367 _pd[2] = dp1;
368
369 const double K892pMass = 0.89176;
370 const double K892pWidth = 0.0503;
371
372 const double K892zeroMass = 0.89555;
373 const double K892zeroWidth = 0.0473;
374
375 const double a980pMass = 0.980;
376 const double c_g1 = 0.341;
377 const double c_g2og1 = 0.892;
378
379 const double K1410zeroMass = 1.421;
380 const double K1410zeroWidth = 0.236;
381
382 const double a1450pMass = 1.474;
383 const double a1450pWidth = 0.265;
384
385 const double SwaveKppi0Mass = 1.441;
386 const double SwaveKppi0Width = 0.193;
387
388 const double SwaveKspi0Mass = 1.441;
389 const double SwaveKspi0Width = 0.193;
390
391 EvtComplex temp( 0.0, 0.0 );
392
393 EvtComplex amp_K892p( 1, 0 );
394 EvtComplex amp_K892zero( -0.3903972065719, 0.1298035433874 );
395 EvtComplex amp_SwaveKppi0( -1.543197997647, 1.30109134697 );
396 EvtComplex amp_SwaveKspi0( -3.123793580183, -0.3449005761848 );
397
398 temp += amp_K892p * ( RBW( 1, K892pMass, K892pWidth, 1 ) );
399 temp += amp_K892zero * ( RBW( 2, K892zeroMass, K892zeroWidth, 1 ) );
400 temp += amp_SwaveKppi0 * ( LASS( 1, SwaveKppi0Mass, SwaveKppi0Width ) );
401 temp += amp_SwaveKspi0 * ( LASS( 2, SwaveKspi0Mass, SwaveKspi0Width ) );
402
403 double ret = pow( abs( temp ), 2 );
404 return ( ret <= 0 ) ? 1e-20 : ret;
405}
p2
double p2[4]
Definition BabayagaNLO.cxx:41
p1
double p1[4]
Definition BabayagaNLO.cxx:41
g1
TF1 * g1
Definition DataBase/tau_mode.c:58
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
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
EvtDToKSKpi0::clone
EvtDecayBase * clone()
Definition EvtDToKSKpi0.cc:35
EvtDToKSKpi0::~EvtDToKSKpi0
virtual ~EvtDToKSKpi0()
Definition EvtDToKSKpi0.cc:31
EvtDToKSKpi0::initProbMax
void initProbMax()
Definition EvtDToKSKpi0.cc:56
EvtDToKSKpi0::getName
void getName(std::string &name)
Definition EvtDToKSKpi0.cc:33
EvtDToKSKpi0::decay
void decay(EvtParticle *p)
Definition EvtDToKSKpi0.cc:58
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
EvtParticle::getP4Lab
EvtVector4R getP4Lab()
Definition EvtParticle.cc:694
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