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