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