BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtD0ToKSKK.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2//
3// Environment:
4// This software is part of the EvtGen package developed jointly
5// for the BaBar and CLEO collaborations. If you use all or part
6// of it, please give an appropriate acknowledgement.
7//
8// Copyright Information: See EvtGen/COPYRIGHT
9// Copyright (C) 1998 Caltech, UCSB
10//
11// Module: EvtD0ToKSKK.cc
12//
13// Description: Routine to handle three-body decays of D0/D0_bar or D+/D- (BESIII
14// arXiv:2006.02800)
15//
16// Modification history:
17//
18// Liaoyuan Dong Dec 11 05:11:33 2022 Module created
19//
20//------------------------------------------------------------------------
21#include "EvtD0ToKSKK.hh"
22#include "../EvtGenBase/EvtConst.hh"
23#include "../EvtGenBase/EvtDecayTable.hh"
24#include "../EvtGenBase/EvtGenKine.hh"
25#include "../EvtGenBase/EvtPDL.hh"
26#include "../EvtGenBase/EvtParticle.hh"
27#include "../EvtGenBase/EvtPatches.hh"
28#include "../EvtGenBase/EvtReport.hh"
29#include <cmath>
30#include <complex>
31#include <stdlib.h>
32using std::endl;
33
34EvtD0ToKSKK::~EvtD0ToKSKK() {}
35
36void EvtD0ToKSKK::getName( std::string& model_name ) { model_name = "D0ToKSKK"; }
37
38EvtDecayBase* EvtD0ToKSKK::clone() { return new EvtD0ToKSKK; }
39
40void EvtD0ToKSKK::init() {
41 checkNArg( 0 );
42 checkNDaug( 3 );
43 checkSpinParent( EvtSpinType::SCALAR );
44 checkSpinDaughter( 0, EvtSpinType::SCALAR );
45 checkSpinDaughter( 1, EvtSpinType::SCALAR );
46 checkSpinDaughter( 2, EvtSpinType::SCALAR );
47
48 // cout << "Initializing D0ToKSKK" << endl;
49
50 MD = 1.86484; // 0 mass D0
51 m1 = 0.4976; // 1 mass KS
52 m2 = 0.49368; // 2 mass K-
53 m3 = 0.49368; // 3 mass K+
54
55 mag_a00 = 1.0;
56 phase_a00 = 0.0;
57 spin_a00 = 0;
58 mass_a00 = 0.994; // 23 -> 0:m23sq 8:cosTheta23_12
59 gA_a00 = 3.80179;
60 gB_a00 = 2.66;
61 gC_a00 = 3.80179;
62 massB1_a00 = 0.547862;
63 massB2_a00 = 0.134977;
64 massC1_a00 = 0.497614;
65 massC2_a00 = 0.497614;
66
67 mag_a0p = 0.591716;
68 phase_a0p = 2.95792;
69 spin_a0p = 0;
70 mass_a0p = 0.994; // 13 -> 1:m13sq 5:cosTheta13_12
71 gA_a0p = 3.80179;
72 gB_a0p = 2.66;
73 massB1_a0p = 0.547862;
74 massB2_a0p = 0.13957;
75
76 mag_phi = 0.713488;
77 phase_phi = 1.66557;
78 spin_phi = 1;
79 mass_phi = 1.01946; // 23 -> 0 8
80 width_phi = 0.004266;
81
82 mag_a2p = 0.126298;
83 phase_a2p = -2.9734;
84 spin_a2p = 2;
85 mass_a2p = 1.3181; // 13 -> 1 5
86 width_a2p = 0.1098;
87
88 mag_a2m = 0.0965778;
89 phase_a2m = -0.138591;
90 spin_a2m = 2;
91 mass_a2m = 1.3181; // 12 -> 2:m12sq 4:cosTheta12_23
92 width_a2m = 0.1098;
93
94 mag_a0_1450m = 0.210031;
95 phase_a0_1450m = -0.23412;
96 spin_a0_1450m = 0;
97 mass_a0_1450m = 1.474; // 12 -> 2 4
98 width_a0_1450m = 0.265;
99
100 mesonRadius = 1.5;
101 motherRadius = 1.0;
102}
103
104void EvtD0ToKSKK::initProbMax() { setProbMax( 143.0 ); }
105
106void EvtD0ToKSKK::decay( EvtParticle* p ) {
107 /*
108 double maxprob = 0.0;
109 for(int ir=0;ir<=60000000;ir++){
110 p->initializePhaseSpace(getNDaug(),getDaugs());
111 EvtVector4R D1 = p->getDaug(0)->getP4();
112 EvtVector4R D2 = p->getDaug(1)->getP4();
113 EvtVector4R D3 = p->getDaug(2)->getP4();
114 double value1 = EvaluateAmp(D1, D2, D3);
115 if(value1>maxprob) {
116 maxprob=value1;
117 cout << "ir = " << ir << " maxprob = " << value1 << endl;
118 }
119 }
120 cout << "MAXprob = " << maxprob << endl;
121 */
122 p->initializePhaseSpace( getNDaug(), getDaugs() );
123 EvtVector4R D1 = p->getDaug( 0 )->getP4(); // KS0
124 EvtVector4R D2 = p->getDaug( 1 )->getP4(); // K-
125 EvtVector4R D3 = p->getDaug( 2 )->getP4(); // K+
126 double value = EvaluateAmp( D1, D2, D3 );
127 setProb( value );
128 return;
129}
130
131double EvtD0ToKSKK::EvaluateAmp( EvtVector4R& p4_d1, EvtVector4R& p4_d2, EvtVector4R& p4_d3 ) {
132 double m23sq = ( p4_d2 + p4_d3 ).mass2();
133 double m13sq = ( p4_d1 + p4_d3 ).mass2();
134 double m12sq = ( p4_d1 + p4_d2 ).mass2();
135 double cosTheta23_1v2 =
136 helicityAngle( MD, m2, m3, m1, m23sq, m12sq ); // Angle in RF12 between 3 and 2
137 double cosTheta13_1v2 =
138 helicityAngle( MD, m1, m3, m2, m13sq, m12sq ); // Angle in RF13 between 2 and 1
139 double cosTheta12_2v3 =
140 helicityAngle( MD, m2, m1, m3, m12sq, m23sq ); // Angle in RF23 between 1 and 2
141 // res[0] = m23sq; res[1] = m13sq; res[2] = m12sq; res[3] = cosTheta23_1v2; res[4] =
142 // cosTheta13_1v2; res[5] = cosTheta12_2v3;
143
144 std::complex<double> Amp_a00 =
145 GetCoefficient( mag_a00, phase_a00 ) *
146 AmpFlatteRes( m23sq, mass_a00, m2, m3, gA_a00, massB1_a00, massB2_a00, gB_a00,
147 massC1_a00, massC2_a00, gC_a00, spin_a00, mesonRadius ) /
148 10.194701039;
149 std::complex<double> Amp_a0p = GetCoefficient( mag_a0p, phase_a0p ) *
150 AmpFlatteRes( m13sq, mass_a0p, m1, m3, gA_a0p, massB1_a0p,
151 massB2_a0p, gB_a0p, spin_a0p, mesonRadius ) /
152 12.861154457;
153 std::complex<double> Amp_phi =
154 GetCoefficient( mag_phi, phase_phi ) *
155 AmpRelBreitWignerRes( m23sq, mass_phi, m2, m3, width_phi, spin_phi, mesonRadius ) *
156 Wigner_d( spin_phi, 0, 0, acos( cosTheta23_1v2 ) ) / 715.1406308838;
157 std::complex<double> Amp_a2p =
158 GetCoefficient( mag_a2p, phase_a2p ) *
159 AmpRelBreitWignerRes( m13sq, mass_a2p, m1, m3, width_a2p, spin_a2p, mesonRadius ) *
160 Wigner_d( spin_a2p, 0, 0, acos( cosTheta13_1v2 ) ) / 340.70301058;
161 std::complex<double> Amp_a2m =
162 GetCoefficient( mag_a2m, phase_a2m ) *
163 AmpRelBreitWignerRes( m12sq, mass_a2m, m1, m2, width_a2m, spin_a2m, mesonRadius ) *
164 Wigner_d( spin_a2m, 0, 0, acos( cosTheta12_2v3 ) ) / 340.737568318;
165 std::complex<double> Amp_a0_1450m =
166 GetCoefficient( mag_a0_1450m, phase_a0_1450m ) *
167 AmpRelBreitWignerRes( m12sq, mass_a0_1450m, m1, m2, width_a0_1450m, spin_a0_1450m,
168 mesonRadius ) /
169 8.9320553176;
170 // res[6] = std::norm(Amp_a00); res[7] = std::norm(Amp_a0p); res[8] = std::norm(Amp_phi);
171 // res[9] = std::norm(Amp_a2p); res[10] = std::norm(Amp_a2m); res[11] =
172 // std::norm(Amp_a0_1450m);
173
174 std::complex<double> Amp = Amp_a00 + Amp_a0p + Amp_phi + Amp_a2p + Amp_a2m + Amp_a0_1450m;
175 double value = std::norm( Amp ); // res[12] = std::norm(Amp);
176 return value;
177}
178
179double EvtD0ToKSKK::helicityAngle( double M, double m, double m2, double mSpec,
180 double invMassSqA, double invMassSqB ) {
181 // Calculate energy and momentum of m1/m2 in the invMassSqA rest frame
182 double eCms = ( invMassSqA + m * m - m2 * m2 ) / ( 2 * sqrt( invMassSqA ) );
183 double pCms = eCms * eCms - m * m;
184 // Calculate energy and momentum of mSpec in the invMassSqA rest frame
185 double eSpecCms = ( M * M - mSpec * mSpec - invMassSqA ) / ( 2 * sqrt( invMassSqA ) );
186 double pSpecCms = eSpecCms * eSpecCms - mSpec * mSpec;
187 double cosAngle = -( invMassSqB - m * m - mSpec * mSpec - 2 * eCms * eSpecCms ) /
188 ( 2 * sqrt( pCms * pSpecCms ) );
189 return cosAngle;
190}
191
192std::complex<double> EvtD0ToKSKK::AmpRelBreitWignerRes( double mSq, double mR, double ma,
193 double mb, double width,
194 unsigned int J, double mesonRadius ) {
195 // J Angular momentum between A&B. In case A&B have spin 0 this is the spin for the resonace.
196 // mSq Invariant mass, mR Resonance mass, ma Mass particle A, mb Mass particle B, width Width
197 // of resonance, mesonRadius Scale of interaction range.
198 std::complex<double> i( 0, 1 );
199 double sqrtS = sqrt( mSq );
200 double barrier =
201 FormFactor( sqrtS, ma, mb, J, mesonRadius ) / FormFactor( mR, ma, mb, J, mesonRadius );
202 std::complex<double> qTerm =
203 pow( ( phspFactor( sqrtS, ma, mb ) / phspFactor( mR, ma, mb ) ) * sqrtS / mR,
204 ( 2 * J + 1 ) ) *
205 mR / sqrtS;
206 std::complex<double> g_final = widthToCoupling( mSq, mR, width, ma, mb, J, mesonRadius );
207 std::complex<double> denom = std::complex<double>( mR * mR - mSq, 0 ) +
208 ( -1.0 ) * i * sqrtS * ( width * qTerm * barrier );
209 std::complex<double> result = g_final / denom;
210 return result;
211}
212
213std::complex<double> EvtD0ToKSKK::widthToCoupling( double mSq, double mR, double width,
214 double ma, double mb, double spin,
215 double mesonRadius ) {
216 double sqrtS = sqrt( mSq );
217 std::complex<double> gammaA( 1, 0 ); // spin==0
218 if ( spin > 0 )
219 {
220 std::complex<double> q = qValue( mR, ma, mb );
221 double ffR = FormFactor( mR, ma, mb, spin, mesonRadius );
222 std::complex<double> qR = pow( q, spin );
223 gammaA = ffR * qR;
224 }
225 std::complex<double> rho = phspFactor( sqrtS, ma, mb );
226 std::complex<double> denom = gammaA * sqrt( rho );
227 std::complex<double> res = std::complex<double>( sqrt( mR * width ), 0 ) / denom;
228 return res;
229}
230
231std::complex<double> EvtD0ToKSKK::AmpFlatteRes( double mSq, double mR, double massA1,
232 double massA2, double gA, double massB1,
233 double massB2, double couplingB,
234 unsigned int J, double mesonRadius ) {
235 std::complex<double> i( 0, 1 );
236 double sqrtS = sqrt( mSq );
237 std::complex<double> termA =
238 gA * gA * phspFactor( sqrtS, massA1, massA2 ); // channel A - signal channel
239 std::complex<double> termB =
240 couplingB * couplingB *
241 phspFactor( sqrtS, massB1, massB2 ); // channel B - hidden channel
242 std::complex<double> denom =
243 std::complex<double>( mR * mR - mSq, 0 ) + ( -1.0 ) * i * ( termA + termB );
244 std::complex<double> result = std::complex<double>( gA, 0 ) / denom;
245 return result;
246}
247std::complex<double> EvtD0ToKSKK::AmpFlatteRes( double mSq, double mR, double massA1,
248 double massA2, double gA, double massB1,
249 double massB2, double couplingB, double massC1,
250 double massC2, double couplingC,
251 unsigned int J, double mesonRadius ) {
252 std::complex<double> i( 0, 1 );
253 double sqrtS = sqrt( mSq );
254 std::complex<double> termA =
255 gA * gA * phspFactor( sqrtS, massA1, massA2 ); // channel A - signal channel
256 std::complex<double> termB =
257 couplingB * couplingB *
258 phspFactor( sqrtS, massB1, massB2 ); // channel B - hidden channel
259 std::complex<double> termC =
260 couplingC * couplingC *
261 phspFactor( sqrtS, massC1, massC2 ); // channel C - hidden channel
262 std::complex<double> denom =
263 std::complex<double>( mR * mR - mSq, 0 ) + ( -1.0 ) * i * ( termA + termB + termC );
264 std::complex<double> result = std::complex<double>( gA, 0 ) / denom;
265 return result;
266}
267
268double EvtD0ToKSKK::FormFactor( double sqrtS, double ma, double mb, double spin,
269 double mesonRadius ) {
270 if ( spin == 0 ) { return 1.0; }
271 std::complex<double> q = qValue( sqrtS, ma, mb );
272 // Blatt-Weisskopt form factors with normalization F(x=mR) = 1. Reference: S.U.Chung Annalen
273 // der Physik 4(1995) 404-430 z = q / (interaction range). For the interaction range we
274 // assume 1/mesonRadius
275 double qSq = std::norm( q );
276 double z = qSq * mesonRadius * mesonRadius;
277 z = fabs( z );
278 if ( spin == 1 ) { return ( sqrt( 2 * z / ( z + 1 ) ) ); }
279 else if ( spin == 2 ) { return ( sqrt( 13 * z * z / ( ( z - 3 ) * ( z - 3 ) + 9 * z ) ) ); }
280 return 0;
281}
282
283std::complex<double> EvtD0ToKSKK::phspFactor( double sqrtS, double ma, double mb ) {
284 // Two types of analytic continuation
285 // 1) Complex sqrt
286 // std::complex<double> rhoOld;
287 // rhoOld = sqrt(std::complex<double>(qSqValue(sqrtS,ma,mb))) / (8*M_PI*sqrtS); //PDG
288 // definition rhoOld = sqrt(std::complex<double>(qSqValue(sqrtS,ma,mb))) / (0.5*sqrtS);
289 // //BaBar definition return rhoOld; //complex sqrt
290
291 // 2) Correct analytic continuation
292 // proper analytic continuation (PDG 2014 - Resonances (47.2.2))
293 // I'm not sure of this is correct for the case of two different masses ma and mb.
294 // Furthermore we divide by the factor 16*Pi*Sqrt[s]). This is more or less arbitrary
295 // and not mentioned in the reference, but it leads to a good agreement between both
296 // approaches.
297 double s = sqrtS * sqrtS;
298 std::complex<double> i( 0, 1 );
299 std::complex<double> rho;
300 double q = sqrt( fabs( qSqValue( sqrtS, ma, mb ) * 4 / s ) );
301 if ( s <= 0 )
302 { // below 0
303 rho = -q / M_PI * log( fabs( ( 1 + q ) / ( 1 - q ) ) );
304 }
305 else if ( 0 < s && s <= ( ma + mb ) * ( ma + mb ) )
306 { // below threshold
307 rho = ( -2.0 * q / M_PI * atan( 1 / q ) ) / ( i * 16.0 * M_PI * sqrtS );
308 }
309 else if ( ( ma + mb ) * ( ma + mb ) < s )
310 { // above threshold
311 rho = ( -q / M_PI * log( fabs( ( 1 + q ) / ( 1 - q ) ) ) + i * q ) /
312 ( i * 16.0 * M_PI * sqrtS );
313 }
314 return rho;
315}
316
317std::complex<double> EvtD0ToKSKK::qValue( double sqrtS, double ma, double mb ) {
318 return phspFactor( sqrtS, ma, mb ) * 8.0 * M_PI * sqrtS;
319}
320
321double EvtD0ToKSKK::qSqValue( double sqrtS, double ma, double mb ) {
322 double mapb = ma + mb;
323 double mamb = ma - mb;
324 double xSq = sqrtS * sqrtS;
325 double t1 = xSq - mapb * mapb;
326 double t2 = xSq - mamb * mamb;
327 return ( t1 * t2 / ( 4 * xSq ) );
328}
329
330std::complex<double> EvtD0ToKSKK::GetCoefficient( double mag, double phase ) const {
331 return std::complex<double>( fabs( mag ) * cos( phase ), fabs( mag ) * sin( phase ) );
332}
333
334double EvtD0ToKSKK::Wigner_d( unsigned int __j, unsigned int __m, unsigned int __n,
335 double __beta ) {
336 int J = (int)( 2. * __j );
337 int M = (int)( 2. * __m );
338 int N = (int)( 2. * __n );
339 int temp_M, k, k_low, k_hi;
340 double const_term = 0.0, sum_term = 0.0, d = 1.0;
341 int m_p_n, j_p_m, j_p_n, j_m_m, j_m_n;
342 int kmn1, kmn2, jmnk, jmk, jnk;
343 double kk;
344
345 if ( J < 0 || abs( M ) > J || abs( N ) > J )
346 {
347 cout << endl;
348 cout << "d: you have entered an illegal number for J, M, N." << endl;
349 cout << "Must follow these rules: J >= 0, abs(M) <= J, and abs(N) <= J." << endl;
350 cout << "J = " << J << " M = " << M << " N = " << N << endl;
351 return 0.;
352 }
353
354 if ( __beta < 0 )
355 {
356 __beta = fabs( __beta );
357 temp_M = M;
358 M = N;
359 N = temp_M;
360 }
361
362 m_p_n = ( M + N ) / 2;
363 j_p_m = ( J + M ) / 2;
364 j_m_m = ( J - M ) / 2;
365 j_p_n = ( J + N ) / 2;
366 j_m_n = ( J - N ) / 2;
367
368 kk = (double)factorial( j_p_m ) * (double)factorial( j_m_m ) * (double)factorial( j_p_n ) *
369 (double)factorial( j_m_n );
370 const_term = pow( ( -1.0 ), ( j_p_m ) ) * sqrt( kk );
371
372 k_low = ( ( m_p_n > 0 ) ? m_p_n : 0 );
373 k_hi = ( ( j_p_m < j_p_n ) ? j_p_m : j_p_n );
374
375 for ( k = k_low; k <= k_hi; k++ )
376 {
377 kmn1 = 2 * k - ( M + N ) / 2;
378 jmnk = J + ( M + N ) / 2 - 2 * k;
379 jmk = ( J + M ) / 2 - k;
380 jnk = ( J + N ) / 2 - k;
381 kmn2 = k - ( M + N ) / 2;
382 sum_term +=
383 pow( ( -1.0 ), ( k ) ) *
384 ( ( pow( cos( __beta / 2.0 ), kmn1 ) ) * ( pow( sin( __beta / 2.0 ), jmnk ) ) ) /
385 ( factorial( k ) * factorial( jmk ) * factorial( jnk ) * factorial( kmn2 ) );
386 }
387
388 d = const_term * sum_term * sqrt( 2.0 * __j + 1.0 );
389 return d;
390}
391
392int EvtD0ToKSKK::factorial( int i ) {
393 int f = 1;
394 if ( ( i == 0 ) || ( i == 1 ) ) f = 1;
395 else
396 {
397 while ( i > 0 )
398 {
399 f = f * i;
400 i--;
401 }
402 }
403 return f;
404}
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
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
M_PI
#define M_PI
Definition TConstant.h:4
EvtD0ToKSKK::initProbMax
void initProbMax()
Definition EvtD0ToKSKK.cc:104
EvtD0ToKSKK::getName
void getName(std::string &name)
Definition EvtD0ToKSKK.cc:36
EvtD0ToKSKK::decay
void decay(EvtParticle *p)
Definition EvtD0ToKSKK.cc:106
EvtD0ToKSKK::~EvtD0ToKSKK
virtual ~EvtD0ToKSKK()
Definition EvtD0ToKSKK.cc:34
EvtD0ToKSKK::clone
EvtDecayBase * clone()
Definition EvtD0ToKSKK.cc:38
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
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
m2
double double * m2
Definition qcdloop1.h:83
m1
double * m1
Definition qcdloop1.h:83