Belle II Software development
EvtDTopipi0Eta.cc
1// Model: EvtDTopipi0Eta
2// This file is an amplitude model for D+ -> pi+ pi0 eta.
3// The model is from the BESIII Collaboration in Phys. Rev. D 110, L111102 (2024). DOI:  https://doi.org/10.1103/PhysRevD.110.L111102
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 <EvtGenBase/EvtPatches.hh>
10#include <EvtGenBase/EvtParticle.hh>
11#include <EvtGenBase/EvtGenKine.hh>
12#include <EvtGenBase/EvtPDL.hh>
13#include <EvtGenBase/EvtReport.hh>
14#include <EvtGenBase/EvtVector4R.hh>
15#include <EvtGenBase/EvtComplex.hh>
16#include <generators/evtgen/models/besiii/EvtDTopipi0Eta.h>
17#include <EvtGenBase/EvtDecayTable.hh>
18#include <stdlib.h>
19
20#include <generators/evtgen/EvtGenModelRegister.h>
21
22using namespace std;
23
24namespace Belle2 {
29
31 B2_EVTGEN_REGISTER_MODEL(EvtDTopipi0Eta);
32
33 EvtDTopipi0Eta::~EvtDTopipi0Eta() {}
34
35 std::string EvtDTopipi0Eta::getName()
36 {
37 return "DTopipi0Eta";
38 }
39
40 EvtDecayBase* EvtDTopipi0Eta::clone()
41 {
42 return new EvtDTopipi0Eta;
43 }
44
45 void EvtDTopipi0Eta::init()
46 {
47 checkNArg(0);
48 checkNDaug(3);
49 checkSpinParent(EvtSpinType::SCALAR);
50 checkSpinDaughter(0, EvtSpinType::SCALAR);
51 checkSpinDaughter(1, EvtSpinType::SCALAR);
52 checkSpinDaughter(2, EvtSpinType::SCALAR);
53
54 phi[0] = -3.3276; rho[0] = 0.31478; //rho eta
55 phi[1] = 0.0; rho[1] = 1.0; //a0+ pi0
56 phi[2] = 0.0; rho[2] = -1.0; //a00 pi+
57
58 mrho = 0.77511;
59 ma0 = 0.99;
60 Grho = 0.1491;
61 Ga0 = 0.0756;
62
63 const double mk0 = 0.497614;
64 const double mass_Kaon = 0.49368;
65 const double mass_Pion = 0.13957;
66 const double mass_Pi0 = 0.1349766;
67 const double meta = 0.547862;
68 mpi = 0.13957;
69 mD = 1.86966;
70 sD = mD * mD;
71 spi = mpi * mpi;
72 snk = mk0 * mk0;
73 sck = mass_Kaon * mass_Kaon;
74 scpi = mass_Pion * mass_Pion;
75 snpi = mass_Pi0 * mass_Pi0;
76 seta = meta * meta;
77
78 pi = 3.1415926;
79
80 ci = EvtComplex(0.0, 1.0);
81 one = EvtComplex(1.0, 0.0);
82
83 int GG[4][4] = { {1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1} };
84 for (int i = 0; i < 4; i++) {
85 for (int j = 0; j < 4; j++) {
86 G[i][j] = GG[i][j];
87 }
88 }
89
90 }
91
92 void EvtDTopipi0Eta::initProbMax()
93 {
94 setProbMax(20.0);
95 }
96
97 void EvtDTopipi0Eta::decay(EvtParticle* p)
98 {
99 p->initializePhaseSpace(getNDaug(), getDaugs());
100 EvtVector4R D1 = p->getDaug(0)->getP4();
101 EvtVector4R D2 = p->getDaug(1)->getP4();
102 EvtVector4R D3 = p->getDaug(2)->getP4();
103
104 double P1[4], P2[4], P3[4];
105 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
106 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
107 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
108
109 double value;
110 value = calDalEva(P1, P2, P3);
111 setProb(value);
112 return;
113
114 }
115
116 double EvtDTopipi0Eta::calDalEva(double P1[], double P2[], double P3[])
117 {
118 //pi+ pi0 eta
119 //0: non-resonance
120 //1: resonance, RBW
121 //2: resonance, GS
122 //3: resonance, Flatte
123 //4: rho-omega mxing for omega
124 EvtComplex PDF[3];
125 EvtComplex cof, pdf, module;
126 double value;
127 PDF[0] = Spin_factor(P1, P2, P3, 1, 2, mrho, Grho); // rho+ eta
128 PDF[1] = Spin_factor(P1, P3, P2, 0, 30, ma0, Ga0); // a0+ pi0
129 PDF[2] = Spin_factor(P2, P3, P1, 0, 31, ma0, Ga0); // a00 pi+
130
131 pdf = EvtComplex(0.0, 0.0);
132 for (int i = 0; i < 3; i++) {
133 cof = EvtComplex(rho[i] * cos(phi[i]), rho[i] * sin(phi[i]));
134 pdf = pdf + cof * PDF[i];
135 }
136 module = conj(pdf) * pdf;
137 value = real(module);
138 return (value <= 0) ? 1e-20 : value;
139 }
140
141 EvtComplex EvtDTopipi0Eta::Spin_factor(double P1[], double P2[], double P3[], int spin, int flag, double mass_R, double width_R)
142 {
143 //D-> R P3, R->P1 P2, 0: non-resonance 1: resonance, RBW 2: resonance, GS 3: resonance, Flatte 4: rho-omega mxing for omega
144 double R[4], s[3], sp2, B[2];
145 double tmp;
146 for (int i = 0; i < 4; i++) {
147 R[i] = P1[i] + P2[i];
148 }
149 s[0] = dot(R, R);
150 s[1] = dot(P1, P1);
151 s[2] = dot(P2, P2);
152 sp2 = dot(P3, P3);
153
154 EvtComplex amp, prop, prop1, prop2;
155 EvtComplex rhokk, rhopieta;
156 if (spin == 0) {
157 if (flag == 0) prop = one;
158 if (flag == 1) prop = propagatorRBW(mass_R, width_R, s[0], s[1], s[2], 3.0, 0);
159 if (flag == 30) {
160 rhokk = Flatte_rhoab(s[0], snk, sck);
161 rhopieta = Flatte_rhoab(s[0], scpi, seta);
162 prop = 1.0 / (mass_R * mass_R - s[0] - ci * (0.341 * rhopieta + 0.341 * 0.892 * rhokk));
163 }
164 if (flag == 31) {
165 double qKsK;
166 qKsK = 0.25 * (s[0] + 3.899750596e-03) * (s[0] + 3.899750596e-03) / s[0] - 0.497614 * 0.497614;
167 if (qKsK > 0) rhokk = 2.0 * sqrt(qKsK / s[0]) * one;
168 if (qKsK < 0) rhokk = 2.0 * sqrt(qKsK / s[0]) * ci;
169 rhopieta = Flatte_rhoab(s[0], snpi, seta);
170 prop = 1.0 / (mass_R * mass_R - s[0] - ci * (0.341 * rhopieta + 0.341 * 0.892 * rhokk));
171 }
172 amp = prop;
173 } else if (spin == 1) {
174 if (flag == 0) {
175 prop = EvtComplex(1.0, 0.0);
176 }
177 if (flag == 1) {
178 prop = propagatorRBW(mass_R, width_R, s[0], s[1], s[2], 3.0, 1);
179 }
180 if (flag == 2) {
181 prop = propagatorGS(mass_R, width_R, s[0], s[1], s[2], 3.0, 1);
182 }
183 if (flag == 4) {
184 prop1 = propagatorGS(mass_R, width_R, s[0], s[1], s[2], 3.0, 1);
185 prop2 = propagatorRBW(0.78266, 0.01358, s[0], s[1], s[2], 3.0, 1);
186 prop = prop1 * prop2;
187 }
188 double T1[4], t1[4];
189 calt1(R, P3, T1);
190 calt1(P1, P2, t1);
191 B[0] = barrier(1, s[0], s[1], s[2], 3.0, mass_R);
192 B[1] = barrier(1, sD, s[0], sp2, 5.0, mD);
193 tmp = 0.0;
194 for (int i = 0; i < 4; i++) {
195 tmp += T1[i] * t1[i] * G[i][i];
196 }
197 amp = tmp * prop * B[0] * B[1];
198 } else if (spin == 2) {
199 double T2[4][4], t2[4][4];
200 calt2(R, P3, T2);
201 calt2(P1, P2, t2);
202 B[0] = barrier(2, s[0], s[1], s[2], 3.0, mass_R);
203 B[1] = barrier(2, sD, s[0], sp2, 5.0, mD);
204 tmp = 0.0;
205 for (int i = 0; i < 4; i++) {
206 for (int j = 0; j < 4; j++) {
207 tmp += T2[i][j] * t2[j][i] * G[j][j] * G[i][i];
208 }
209 }
210 if (flag == 0) prop = one;
211 if (flag == 1) prop = propagatorRBW(mass_R, width_R, s[0], s[1], s[2], 3.0, 2);
212 amp = tmp * prop * B[0] * B[1];
213 }
214 return amp;
215 }
216
217 double EvtDTopipi0Eta::dot(double* a1, double* a2)
218 {
219 double Dot = 0;
220 for (int i = 0; i != 4; i++) {
221 Dot += a1[i] * a2[i] * G[i][i];
222 }
223 return Dot;
224 }
225
226 double EvtDTopipi0Eta::Qabcs(double sa, double sb, double sc)
227 {
228 double Qabcs = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
229 if (Qabcs < 0) Qabcs = 1e-16;
230 return Qabcs;
231 }
232
233 double EvtDTopipi0Eta::barrier(double l, double sa, double sb, double sc, double r, double mass)
234 {
235 double sa0 = mass * mass;
236 double q0 = Qabcs(sa0, sb, sc);
237 double z0 = q0 * r * r;
238 double q = Qabcs(sa, sb, sc);
239 q = sqrt(q);
240 double z = q * r;
241 z = z * z;
242 double F = 1;
243 if (l > 2) F = 0;
244 if (l == 0) F = 1;
245 if (l == 1) F = sqrt((1 + z0) / (1 + z));
246 if (l == 2) F = sqrt((9 + 3 * z0 + z0 * z0) / (9 + 3 * z + z * z));
247 return F;
248 }
249
250 void EvtDTopipi0Eta::calt1(double daug1[], double daug2[], double t1[])
251 {
252 double p, pq;
253 double pa[4], qa[4];
254 for (int i = 0; i != 4; i++) {
255 pa[i] = daug1[i] + daug2[i];
256 qa[i] = daug1[i] - daug2[i];
257 }
258 p = dot(pa, pa);
259 pq = dot(pa, qa);
260 for (int i = 0; i != 4; i++) {
261 t1[i] = qa[i] - pq / p * pa[i];
262 }
263 }
264
265 void EvtDTopipi0Eta::calt2(double daug1[], double daug2[], double t2[][4])
266 {
267 double p, r;
268 double pa[4], t1[4];
269 calt1(daug1, daug2, t1);
270 r = dot(t1, t1);
271 for (int i = 0; i != 4; i++) {
272 pa[i] = daug1[i] + daug2[i];
273 }
274 p = dot(pa, pa);
275 for (int i = 0; i != 4; i++) {
276 for (int j = 0; j != 4; j++) {
277 t2[i][j] = t1[i] * t1[j] - 1.0 / 3 * r * (G[i][j] - pa[i] * pa[j] / p);
278 }
279 }
280 }
281
282 double EvtDTopipi0Eta::wid(double mass, double sa, double sb, double sc, double r, int l)
283 {
284 double widm(0.), q(0.), q0(0.);
285 double sa0 = mass * mass;
286 double m = sqrt(sa);
287 q = Qabcs(sa, sb, sc);
288 q0 = Qabcs(sa0, sb, sc);
289 double z, z0;
290 z = q * r * r;
291 z0 = q0 * r * r;
292 double t = q / q0;
293 double F(0.);
294 if (l == 0) F = 1;
295 if (l == 1) F = sqrt((1 + z0) / (1 + z));
296 if (l == 2) F = sqrt((9 + 3 * z0 + z0 * z0) / (9 + 3 * z + z * z));
297 widm = pow(t, l + 0.5) * mass / m * F * F;
298 return widm;
299 }
300
301 EvtComplex EvtDTopipi0Eta::propagatorRBW(double mass, double width, double sa, double sb, double sc, double r, int l)
302 {
303 EvtComplex prop = 1.0 / (mass * mass - sa - ci * mass * width * wid(mass, sa, sb, sc, r, l));
304 return prop;
305 }
306
307 double EvtDTopipi0Eta::h(double m, double q)
308 {
309 double h = 2.0 / pi * q / m * log((m + 2 * q) / (0.13957 + 0.134976));
310 return h;
311 }
312
313 double EvtDTopipi0Eta::dh(double mass, double q0)
314 {
315 double dh = h(mass, q0) * (1.0 / (8 * q0 * q0) - 1.0 / (2 * mass * mass)) + 1.0 / (2 * pi * mass * mass);
316 return dh;
317 }
318
319 double EvtDTopipi0Eta::f(double mass, double sx, double q0, double q)
320 {
321 double m = sqrt(sx);
322 double f = mass * mass / (pow(q0, 3)) * (q * q * (h(m, q) - h(mass, q0)) + (mass * mass - sx) * q0 * q0 * dh(mass, q0));
323 return f;
324 }
325
326 double EvtDTopipi0Eta::d(double mass, double q0)
327 {
328 double cmpi = 0.5 * (0.13957 + 0.134976);
329 double mpi2 = cmpi * cmpi;
330 double d = 3.0 / pi * mpi2 / (q0 * q0) * log((mass + 2 * q0) / (2 * cmpi)) + mass / (2 * pi * q0) - (mpi2 * mass) / (pi * pow(q0,
331 3));
332 return d;
333 }
334
335 EvtComplex EvtDTopipi0Eta::propagatorGS(double mass, double width, double sa, double sb, double sc, double r, int l)
336 {
337 double q = Qabcs(sa, sb, sc);
338 double sa0 = mass * mass;
339 double q0 = Qabcs(sa0, sb, sc);
340 q = sqrt(q);
341 q0 = sqrt(q0);
342 EvtComplex prop = (1 + d(mass, q0) * width / mass) / (mass * mass - sa + width * f(mass, sa, q0, q) - ci * mass * width * wid(mass,
343 sa, sb, sc, r, l));
344 return prop;
345 }
346
347 EvtComplex EvtDTopipi0Eta::Flatte_rhoab(double sa, double sb, double sc)
348 {
349 double q = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
350 EvtComplex rho_val;
351 if (q > 0) {
352 rho_val = 2.0 * sqrt(q / sa) * one;
353 }
354 if (q < 0) {
355 rho_val = 2.0 * sqrt(-q / sa) * ci;
356 }
357 return rho_val;
358 }
359
360 EvtComplex EvtDTopipi0Eta::propagatorFlatte(double mass, double width __attribute__((unused)), double sx, double* sb, double* sc)
361 {
362 const double g1sq = 0.5468 * 0.5468;
363 const double g2sq = 0.23 * 0.23;
364 EvtComplex rho1 = Flatte_rhoab(sx, sb[0], sc[0]);
365 EvtComplex rho2 = Flatte_rhoab(sx, sb[1], sc[1]);
366 EvtComplex prop = 1.0 / (mass * mass - sx - ci * (g1sq * rho1 + g2sq * rho2));
367 return prop;
368 }
369
371} // Belle2 namespace
R
double R
typedef autogenerated by FFTW
Definition DiscreteCosineTransform_31points.cc:35
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.