Belle II Software development
EvtDToKSKpi0.cc
1// Model: EvtDToKSKpi0
2// This file is an amplitude model for D+ -> K_S0 K+ pi0.
3// The model is from the BESIII Collaboration in PRD 104, 012006 (2021). DOI:  https://doi.org/10.1103/PhysRevD.104.012006
4//
5// Permission to include these files in basf2 was generously granted by the BESIII Collaboration.
6//
7// Please cite the original reference for any public/published results where this model was used.
8
9#include <iostream>
10#include <cmath>
11#include <string>
12#include <EvtGenBase/EvtCPUtil.hh>
13#include <EvtGenBase/EvtTensor4C.hh>
14#include <EvtGenBase/EvtPatches.hh>
15#include <stdlib.h>
16#include <EvtGenBase/EvtParticle.hh>
17#include <EvtGenBase/EvtGenKine.hh>
18#include <EvtGenBase/EvtPDL.hh>
19#include <EvtGenBase/EvtVector4R.hh>
20#include <EvtGenBase/EvtReport.hh>
21
22#include <generators/evtgen/EvtGenModelRegister.h>
23#include <generators/evtgen/models/besiii/EvtDToKSKpi0.h>
24#include <string>
25
26using namespace std;
27
28namespace Belle2 {
33
35 B2_EVTGEN_REGISTER_MODEL(EvtDToKSKpi0);
36
37 EvtDToKSKpi0::~EvtDToKSKpi0() {}
38
39 std::string EvtDToKSKpi0::getName()
40 {
41 return "DToKSKpi0";
42 }
43
44 EvtDecayBase* EvtDToKSKpi0::clone()
45 {
46 return new EvtDToKSKpi0;
47 }
48
49 void EvtDToKSKpi0::init()
50 {
51 checkNArg(0);
52 checkNDaug(3);
53 checkSpinParent(EvtSpinType::SCALAR);
54
55 _nd = 3;
56 _mDp = 1.86965;
57 c_motherMass = _mDp;
58 _mDp2 = _mDp * _mDp;
59 _mDp2inv = 1. / _mDp2;
60 KsMass = 0.497611;
61 KpMass = 0.493677;
62 pi0Mass = 0.134977;
63 etamass = 0.547862;
64 pipMass = 0.13957061;
65 c_meson_radius_inter = 1.5;
66 c_meson_radius_Dp = 5;
67 }
68 void EvtDToKSKpi0::initProbMax()
69 {
70 setProbMax(1343.2);
71 }
72
73 void EvtDToKSKpi0::decay(EvtParticle* p)
74 {
75 p->initializePhaseSpace(getNDaug(), getDaugs());
76 for (int i = 0; i < _nd; i++) {
77 _p4Lab[i] = p->getDaug(i)->getP4Lab();
78 _p4CM[i] = p->getDaug(i)->getP4();
79 }
80 double prob = AmplitudeSquare();
81 setProb(prob);
82 return;
83 }
84
85 double EvtDToKSKpi0::twoBodyCMmom(double rMassSq, double d1m, double d2m)
86 {
87 double kin1 = 1 - pow(d1m + d2m, 2) / rMassSq;
88 kin1 = kin1 >= 0 ? sqrt(kin1) : 1;
89 double kin2 = 1 - pow(d1m - d2m, 2) / rMassSq;
90 kin2 = kin2 >= 0 ? sqrt(kin2) : 1;
91
92 double ret = 0.5 * sqrt(rMassSq) * kin1 * kin2;
93 return ret;
94 }
95
96 double EvtDToKSKpi0::dampingFactorSquare(const double& cmmom, const int& spin, const double& mRadius)
97 {
98 double square = mRadius * mRadius * cmmom * cmmom;
99 double dfsq = 1 + square;
100 double dfsqres = dfsq + 8 + 2 * square + square * square;
101
102 double ret = (spin == 2) ? dfsqres : dfsq;
103 return ret;
104 }
105
106 double EvtDToKSKpi0::spinFactor(int spin, double motherMass, double daug1Mass, double daug2Mass, double daug3Mass,
107 double m12, double m13, double m23)
108 {
109 if (spin == 0) return 1;
110
111 double _mA = daug1Mass;
112 double _mB = daug2Mass;
113 double _mC = daug3Mass;
114 double _mAB = m12;
115 double _mAC = m13;
116 double _mBC = m23;
117
118 double massFactor = 1.0 / _mAB;
119 double sFactor = -1;
120 sFactor *= ((_mBC - _mAC) + (massFactor * (motherMass * motherMass - _mC * _mC) * (_mA * _mA - _mB * _mB)));
121
122 if (spin == 2) {
123 sFactor *= sFactor;
124 double extraterm = ((_mAB - (2 * motherMass * motherMass) - (2 * _mC * _mC))
125 + massFactor * pow(motherMass * motherMass - _mC * _mC, 2));
126 extraterm *= ((_mAB - (2 * _mA * _mA) - (2 * _mB * _mB)) + massFactor * pow(_mA * _mA - _mB * _mB, 2));
127 extraterm /= 3;
128 sFactor -= extraterm;
129 }
130
131 return sFactor;
132 }
133
134 EvtComplex EvtDToKSKpi0::RBW(int id, double resmass, double reswidth, int spin)
135 {
136 double resmass2 = pow(resmass, 2);
137
138 EvtVector4R p1, p2, p3;
139 double mass_daug1 = 0, mass_daug2 = 0, mass_daug3 = 0;
140 if (id == 1) {
141 p1 = _pd[0];
142 mass_daug1 = pi0Mass;
143 p2 = _pd[1];
144 mass_daug2 = KpMass;
145 p3 = _pd[2];
146 mass_daug3 = KsMass;
147 }
148 if (id == 2) {
149 p1 = _pd[2];
150 mass_daug1 = KsMass;
151 p2 = _pd[0];
152 mass_daug2 = pi0Mass;
153 p3 = _pd[1];
154 mass_daug3 = KpMass;
155 }
156 if (id == 3) {
157 p1 = _pd[1];
158 mass_daug1 = KpMass;
159 p2 = _pd[2];
160 mass_daug2 = KsMass;
161 p3 = _pd[0];
162 mass_daug3 = pi0Mass;
163 }
164
165 double rMassSq = (p1 + p2).mass2();
166 double m12 = (p1 + p2).mass2();
167 double m13 = (p1 + p3).mass2();
168 double m23 = (p2 + p3).mass2();
169
170 double rMass = sqrt(rMassSq);
171 double frFactor = 1;
172 double fdFactor = 1;
173
174 double measureDaughterMoms = twoBodyCMmom(rMassSq, mass_daug1, mass_daug2);
175 double nominalDaughterMoms = twoBodyCMmom(resmass2, mass_daug1, mass_daug2);
176
177 if (spin != 0) {
178 frFactor = dampingFactorSquare(nominalDaughterMoms, spin, c_meson_radius_inter) / dampingFactorSquare(measureDaughterMoms, spin,
179 c_meson_radius_inter);
180
181 double measureDaughterMoms2 = twoBodyCMmom(c_motherMass * c_motherMass, rMass, mass_daug3);
182 double nominalDaughterMoms2 = twoBodyCMmom(c_motherMass * c_motherMass, resmass, mass_daug3);
183 fdFactor = dampingFactorSquare(nominalDaughterMoms2, spin, c_meson_radius_Dp) / dampingFactorSquare(measureDaughterMoms2, spin,
184 c_meson_radius_Dp);
185 }
186 double A = (resmass2 - rMassSq);
187 double B = resmass2 * reswidth * pow(measureDaughterMoms / nominalDaughterMoms, 2.0 * spin + 1) * frFactor / sqrt(rMassSq);
188 double C = 1.0 / (pow(A, 2) + pow(B, 2));
189
190 EvtComplex ret(A * C, B * C);
191 ret *= sqrt(frFactor * fdFactor);
192 ret *= spinFactor(spin, c_motherMass, mass_daug1, mass_daug2, mass_daug3, m12, m13, m23);
193 return ret;
194 }
195
196 EvtComplex EvtDToKSKpi0::Flatte(int id, double resmass, double g1, double rg2og1)
197 {
198
199 EvtVector4R p1, p2, p3;
200 double mass_daug1 __attribute__((unused)), mass_daug2 __attribute__((unused)), mass_daug3 __attribute__((unused));
201 if (id == 1) {
202 p1 = _pd[0];
203 mass_daug1 = pi0Mass;
204 p2 = _pd[1];
205 mass_daug2 = KpMass;
206 p3 = _pd[2];
207 mass_daug3 = KsMass;
208 }
209 if (id == 2) {
210 p1 = _pd[2];
211 mass_daug1 = KsMass;
212 p2 = _pd[0];
213 mass_daug2 = pi0Mass;
214 p3 = _pd[1];
215 mass_daug3 = KpMass;
216 }
217 if (id == 3) {
218 p1 = _pd[1];
219 mass_daug1 = KpMass;
220 p2 = _pd[2];
221 mass_daug2 = KsMass;
222 p3 = _pd[0];
223 mass_daug3 = pi0Mass;
224 }
225
226 double rMassSq __attribute__((unused)) = (p1 + p2).mass2();
227 double m12 = (p1 + p2).mass2();
228
229 double s = m12;
230
231 double rhoetapi = 2 * twoBodyCMmom(s, KsMass, KpMass) / sqrt(s);
232 double rhoKKbar = 2 * twoBodyCMmom(s, etamass, pipMass) / sqrt(s);
233 double img = rhoetapi * g1 + rhoKKbar * g1 * rg2og1;
234
235 EvtComplex ret = EvtComplex(1, 0) / EvtComplex(resmass * resmass - s, -img);
236 return ret;
237
238 }
239
240 EvtComplex EvtDToKSKpi0::LASS(int id, double resmass, double reswidth)
241 {
242 int spin = 0;
243 double resmass2 = pow(resmass, 2);
244
245 EvtVector4R p1, p2, p3;
246 double mass_daug1 = 0, mass_daug2 = 0, mass_daug3 = 0;
247 if (id == 1) {
248 p1 = _pd[0];
249 mass_daug1 = pi0Mass;
250 p2 = _pd[1];
251 mass_daug2 = KpMass;
252 p3 = _pd[2];
253 mass_daug3 = KsMass;
254 }
255 if (id == 2) {
256 p1 = _pd[2];
257 mass_daug1 = KsMass;
258 p2 = _pd[0];
259 mass_daug2 = pi0Mass;
260 p3 = _pd[1];
261 mass_daug3 = KpMass;
262 }
263 if (id == 3) {
264 p1 = _pd[1];
265 mass_daug1 = KpMass;
266 p2 = _pd[2];
267 mass_daug2 = KsMass;
268 p3 = _pd[0];
269 mass_daug3 = pi0Mass;
270 }
271
272 double rMassSq = (p1 + p2).mass2();
273 double m12 = (p1 + p2).mass2();
274 double m13 = (p1 + p3).mass2();
275 double m23 = (p2 + p3).mass2();
276
277 double rMass = sqrt(rMassSq);
278 double frFactor = 1;
279 double fdFactor = 1;
280
281 double measureDaughterMoms = twoBodyCMmom(rMassSq, mass_daug1, mass_daug2);
282 double nominalDaughterMoms = twoBodyCMmom(resmass2, mass_daug1, mass_daug2);
283 if (spin != 0) {
284 frFactor = dampingFactorSquare(nominalDaughterMoms, spin, c_meson_radius_inter) / dampingFactorSquare(measureDaughterMoms, spin,
285 c_meson_radius_inter);
286 double measureDaughterMoms2 = twoBodyCMmom(c_motherMass * c_motherMass, rMass, mass_daug3);
287 double nominalDaughterMoms2 = twoBodyCMmom(c_motherMass * c_motherMass, resmass, mass_daug3);
288 fdFactor = dampingFactorSquare(nominalDaughterMoms2, spin, c_meson_radius_Dp) / dampingFactorSquare(measureDaughterMoms2, spin,
289 c_meson_radius_Dp);
290 }
291
292 double q = measureDaughterMoms;
293 double g = reswidth * pow(measureDaughterMoms / nominalDaughterMoms, 2.0 * spin + 1) * frFactor / sqrt(rMassSq);
294 g *= resmass;
295
296 const double _a = 0.113;
297 const double _r = -33.8;
298 const double _R = 1;
299 const double _phiR = -109.7 * 3.141592653 / 180.0;
300 const double _B = 0.96;
301 const double _phiB = 0.1 * 3.141592653 / 180.0;
302
303 double cot_deltaB = (1.0 / (_a * q)) + 0.5 * _r * q;
304 double qcot_deltaB = (1.0 / _a) + 0.5 * _r * q * q;
305
306 EvtComplex expi2deltaB = EvtComplex(qcot_deltaB, q) / EvtComplex(qcot_deltaB, -q);
307 EvtComplex resT = EvtComplex(cos(_phiR + 2 * _phiB), sin(_phiR + 2 * _phiB)) * _R;
308
309 EvtComplex prop = EvtComplex(1, 0) / EvtComplex(resmass2 - rMassSq, -(resmass) * g);
310 resT *= prop * (resmass2 * reswidth / nominalDaughterMoms) * expi2deltaB;
311
312 resT += EvtComplex(cos(_phiB), sin(_phiB)) * _B * (cos(_phiB) + cot_deltaB * sin(_phiB)) * sqrt(rMassSq) / EvtComplex(qcot_deltaB,
313 -q);
314
315 resT *= sqrt(frFactor * fdFactor);
316 resT *= spinFactor(spin, c_motherMass, mass_daug1, mass_daug2, mass_daug3, m12, m13, m23);
317
318 return resT;
319 }
320
321 double EvtDToKSKpi0::AmplitudeSquare()
322 {
323 EvtVector4R dp1 = GetDaugMomLab(0);
324 EvtVector4R dp2 = GetDaugMomLab(1);
325 EvtVector4R dp3 = GetDaugMomLab(2);
326 _pd[0] = dp3;
327 _pd[1] = dp2;
328 _pd[2] = dp1;
329
330 const double K892pMass = 0.89176;
331 const double K892pWidth = 0.0503;
332
333 const double K892zeroMass = 0.89555;
334 const double K892zeroWidth = 0.0473;
335
336 const double SwaveKppi0Mass = 1.441;
337 const double SwaveKppi0Width = 0.193;
338
339 const double SwaveKspi0Mass = 1.441;
340 const double SwaveKspi0Width = 0.193;
341
342 EvtComplex temp(0.0, 0.0);
343
344 EvtComplex amp_K892p(1, 0);
345 EvtComplex amp_K892zero(-0.3903972065719, 0.1298035433874);
346 EvtComplex amp_SwaveKppi0(-1.543197997647, 1.30109134697);
347 EvtComplex amp_SwaveKspi0(-3.123793580183, -0.3449005761848);
348
349 temp += amp_K892p * (RBW(1, K892pMass, K892pWidth, 1));
350 temp += amp_K892zero * (RBW(2, K892zeroMass, K892zeroWidth, 1));
351 temp += amp_SwaveKppi0 * (LASS(1, SwaveKppi0Mass, SwaveKppi0Width));
352 temp += amp_SwaveKspi0 * (LASS(2, SwaveKspi0Mass, SwaveKspi0Width));
353
354 double ret = pow(abs(temp), 2);
355 return (ret <= 0) ? 1e-20 : ret;
356 }
357
359} // Belle2 namespace
B2_EVTGEN_REGISTER_MODEL
#define B2_EVTGEN_REGISTER_MODEL(classname)
Class to register B2_EVTGEN_REGISTER_MODEL.
Definition EvtGenModelRegister.h:75
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
Belle2
Abstract base class for different kinds of events.
Definition MillepedeAlgorithm.h:17
std
STL namespace.