Belle II Software development
EvtD0TopipiEta.cc
1// Model: EvtD0TopipiEta
2// This file is an amplitude model for D0 -> pi- pi+ 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 <EvtGenBase/EvtDecayTable.hh>
17#include <stdlib.h>
18
19#include <generators/evtgen/EvtGenModelRegister.h>
20#include <generators/evtgen/models/besiii/EvtD0TopipiEta.h>
21
22namespace Belle2 {
27
28 using namespace std;
29
30 B2_EVTGEN_REGISTER_MODEL(EvtD0TopipiEta);
31
32 EvtD0TopipiEta::~EvtD0TopipiEta() {}
33
34 std::string EvtD0TopipiEta::getName()
35 {
36 return "D0TopipiEta";
37 }
38
39 EvtDecayBase* EvtD0TopipiEta::clone()
40 {
41 return new EvtD0TopipiEta;
42 }
43
44 void EvtD0TopipiEta::init()
45 {
46 // check that there are 0 arguments
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] = 0; rho[0] = 1; //rho eta
55 phi[1] = -0.98109; rho[1] = -0.02447; //omega eta (rho-omega mixing)
56 phi[2] = 0.71358; rho[2] = 1.0848; //a0- pi+
57 phi[3] = -0.83115; rho[3] = 2.6444; //a0+ pi-
58 phi[4] = -0.058521; rho[4] = 7.0274; //(pi+ eta)_{2+} pi-
59
60 mrho = 0.77511;
61 ma0 = 0.99;
62 Grho = 0.1491;
63 Ga0 = 0.0756;
64
65 const double mk0 = 0.497614;
66 const double mass_Kaon = 0.49368;
67 const double mass_Pion = 0.13957;
68 const double mass_Pi0 = 0.1349766;
69 const double meta = 0.547862;
70 mpi = 0.13957;
71 mD = 1.86483;
72 sD = mD * mD;
73 spi = mpi * mpi;
74 snk = mk0 * mk0;
75 sck = mass_Kaon * mass_Kaon;
76 scpi = mass_Pion * mass_Pion;
77 snpi = mass_Pi0 * mass_Pi0;
78 seta = meta * meta;
79
80 pi = 3.1415926;
81
82 ci = EvtComplex(0.0, 1.0);
83 one = EvtComplex(1.0, 0.0);
84
85 int GG[4][4] = { {1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1} };
86 for (int i = 0; i < 4; i++) {
87 for (int j = 0; j < 4; j++) {
88 G[i][j] = GG[i][j];
89 }
90 }
91
92
93 }
94
95 void EvtD0TopipiEta::initProbMax()
96 {
97 setProbMax(476.5);
98 }
99
100 void EvtD0TopipiEta::decay(EvtParticle* p)
101 {
102 /*
103 // This piece of code could in principle be used to calculate maximum
104 // probablity on fly. But as it uses high number of points and model
105 // deals with single final state, we keep hardcoded number for now rather
106 // than adapting code to work here.
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
115 double P1[4], P2[4], P3[4];
116 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
117 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
118 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
119
120 double value;
121 value = calDalEva(P1, P2, P3);
122 if(value>maxprob) {
123 maxprob=value;
124 }
125 }
126 */
127 p->initializePhaseSpace(getNDaug(), getDaugs());
128 EvtVector4R D1 = p->getDaug(0)->getP4();
129 EvtVector4R D2 = p->getDaug(1)->getP4();
130 EvtVector4R D3 = p->getDaug(2)->getP4();
131
132 double P1[4], P2[4], P3[4];
133 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
134 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
135 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
136
137 double value;
138 value = calDalEva(P1, P2, P3);
139 setProb(value);
140
141 return ;
142 }
143
144 double EvtD0TopipiEta::calDalEva(double P1[], double P2[], double P3[])
145 {
146 //pi- pi+ eta
147 //0: non-resonance
148 //1: resonance, RBW
149 //2: resonance, GS
150 //3: resonance, Flatte
151 //4: rho-omega mxing for omega
152 EvtComplex PDF[6];
153 EvtComplex cof, pdf, module;
154 double value;
155 PDF[0] = Spin_factor(P1, P2, P3, 1, 2, mrho, Grho); // rho eta
156 PDF[1] = Spin_factor(P1, P2, P3, 1, 4, mrho, Grho); // rho-omega mixing
157 PDF[2] = Spin_factor(P1, P3, P2, 0, 3, ma0, Ga0); // a0- pi+
158 PDF[3] = Spin_factor(P2, P3, P1, 0, 3, ma0, Ga0); // a0+ pi-
159 PDF[4] = Spin_factor(P2, P3, P1, 2, 0, 1.698, 0.265); // pi+ eta 2+ non-res
160
161 pdf = EvtComplex(0.0, 0.0);
162 for (int i = 0; i < 5; i++) {
163 cof = EvtComplex(rho[i] * cos(phi[i]), rho[i] * sin(phi[i]));
164 pdf = pdf + cof * PDF[i];
165 }
166 module = conj(pdf) * pdf;
167 value = real(module);
168 return (value <= 0) ? 1e-20 : value;
169 }
170
171 EvtComplex EvtD0TopipiEta::Spin_factor(double P1[], double P2[], double P3[], int spin, int flag, double mass_R, double width_R)
172 {
173 //D-> R P3, R->P1 P2, 0: non-resonance 1: resonance, RBW 2: resonance, GS 3: resonance, Flatte 4: rho-omega mxing for omega
174 double R[4], s[3], sp2, B[2];
175 double tmp;
176 for (int i = 0; i < 4; i++) {
177 R[i] = P1[i] + P2[i];
178 }
179 s[0] = dot(R, R);
180 s[1] = dot(P1, P1);
181 s[2] = dot(P2, P2);
182 sp2 = dot(P3, P3);
183
184 EvtComplex amp, prop, prop1, prop2;
185
186 //-----------for prop-------------------------
187 EvtComplex rhokk, rhopieta;
188 if (spin == 0) {
189 if (flag == 0) prop = one;
190 if (flag == 1) prop = propagatorRBW(mass_R, width_R, s[0], s[1], s[2], 3.0, 0);
191 if (flag == 3) {
192 rhokk = Flatte_rhoab(s[0], snk, sck);
193 rhopieta = Flatte_rhoab(s[0], scpi, seta);
194 prop = 1.0 / (mass_R * mass_R - s[0] - ci * (0.341 * rhopieta + 0.341 * 0.892 * rhokk));
195 }
196 amp = prop;
197 } else if (spin == 1) {
198 if (flag == 0) {
199 prop = EvtComplex(1.0, 0.0);
200 }
201 if (flag == 1) {
202 prop = propagatorRBW(mass_R, width_R, s[0], s[1], s[2], 3.0, 1);
203 }
204 if (flag == 2) {
205 prop = propagatorGS(mass_R, width_R, s[0], s[1], s[2], 3.0, 1);
206 }
207 if (flag == 4) {
208 prop1 = propagatorGS(mass_R, width_R, s[0], s[1], s[2], 3.0, 1);
209 prop2 = propagatorRBW(0.78266, 0.01358, s[0], s[1], s[2], 3.0, 1);
210 prop = prop1 * prop2;
211 }
212 double T1[4], t1[4];
213 calt1(R, P3, T1);
214 calt1(P1, P2, t1);
215 B[0] = barrier(1, s[0], s[1], s[2], 3.0, mass_R);
216 B[1] = barrier(1, sD, s[0], sp2, 5.0, mD);
217 tmp = 0.0;
218 for (int i = 0; i < 4; i++) {
219 tmp += T1[i] * t1[i] * G[i][i];
220 }
221 amp = tmp * prop * B[0] * B[1];
222 } else if (spin == 2) {
223 double T2[4][4], t2[4][4];
224 calt2(R, P3, T2);
225 calt2(P1, P2, t2);
226 B[0] = barrier(2, s[0], s[1], s[2], 3.0, mass_R);
227 B[1] = barrier(2, sD, s[0], sp2, 5.0, mD);
228 tmp = 0.0;
229 for (int i = 0; i < 4; i++) {
230 for (int j = 0; j < 4; j++) {
231 tmp += T2[i][j] * t2[j][i] * G[j][j] * G[i][i];
232 }
233 }
234 if (flag == 0) prop = one;
235 if (flag == 1) prop = propagatorRBW(mass_R, width_R, s[0], s[1], s[2], 3.0, 2);
236 amp = tmp * prop * B[0] * B[1];
237 }
238 return amp;
239 }
240
241 double EvtD0TopipiEta::dot(double* a1, double* a2)
242 {
243 double Dot = 0;
244 for (int i = 0; i != 4; i++) {
245 Dot += a1[i] * a2[i] * G[i][i];
246 }
247 return Dot;
248 }
249
250 double EvtD0TopipiEta::Qabcs(double sa, double sb, double sc)
251 {
252 double Qabcs = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
253 if (Qabcs < 0) Qabcs = 1e-16;
254 return Qabcs;
255 }
256
257 double EvtD0TopipiEta::barrier(double l, double sa, double sb, double sc, double r, double mass)
258 {
259 double sa0 = mass * mass;
260 double q0 = Qabcs(sa0, sb, sc);
261 double z0 = q0 * r * r;
262 double q = Qabcs(sa, sb, sc);
263 q = sqrt(q);
264 double z = q * r;
265 z = z * z;
266 double F = 1;
267 if (l > 2) F = 0;
268 if (l == 0) F = 1;
269 if (l == 1) F = sqrt((1 + z0) / (1 + z));
270 if (l == 2) F = sqrt((9 + 3 * z0 + z0 * z0) / (9 + 3 * z + z * z));
271 return F;
272 }
273
274 void EvtD0TopipiEta::calt1(double daug1[], double daug2[], double t1[])
275 {
276 double p, pq;
277 double pa[4], qa[4];
278 for (int i = 0; i != 4; i++) {
279 pa[i] = daug1[i] + daug2[i];
280 qa[i] = daug1[i] - daug2[i];
281 }
282 p = dot(pa, pa);
283 pq = dot(pa, qa);
284 for (int i = 0; i != 4; i++) {
285 t1[i] = qa[i] - pq / p * pa[i];
286 }
287 }
288
289 void EvtD0TopipiEta::calt2(double daug1[], double daug2[], double t2[][4])
290 {
291 double p, r;
292 double pa[4], t1[4];
293 calt1(daug1, daug2, t1);
294 r = dot(t1, t1);
295 for (int i = 0; i != 4; i++) {
296 pa[i] = daug1[i] + daug2[i];
297 }
298 p = dot(pa, pa);
299 for (int i = 0; i != 4; i++) {
300 for (int j = 0; j != 4; j++) {
301 t2[i][j] = t1[i] * t1[j] - 1.0 / 3 * r * (G[i][j] - pa[i] * pa[j] / p);
302 }
303 }
304 }
305
306 double EvtD0TopipiEta::wid(double mass, double sa, double sb, double sc, double r, int l)
307 {
308 double widm(0.), q(0.), q0(0.);
309 double sa0 = mass * mass;
310 double m = sqrt(sa);
311 q = Qabcs(sa, sb, sc);
312 q0 = Qabcs(sa0, sb, sc);
313 double z, z0;
314 z = q * r * r;
315 z0 = q0 * r * r;
316 double t = q / q0;
317 double F(0.);
318 if (l == 0) F = 1;
319 if (l == 1) F = sqrt((1 + z0) / (1 + z));
320 if (l == 2) F = sqrt((9 + 3 * z0 + z0 * z0) / (9 + 3 * z + z * z));
321 widm = pow(t, l + 0.5) * mass / m * F * F;
322 return widm;
323 }
324
325 EvtComplex EvtD0TopipiEta::propagatorRBW(double mass, double width, double sa, double sb, double sc, double r, int l)
326 {
327 EvtComplex prop = 1.0 / (mass * mass - sa - ci * mass * width * wid(mass, sa, sb, sc, r, l));
328 return prop;
329 }
330
331 double EvtD0TopipiEta::h(double m, double q)
332 {
333 double h(0.);
334 h = 2 / pi * q / m * log((m + 2 * q) / (2 * mpi));
335 return h;
336 }
337
338 double EvtD0TopipiEta::dh(double mass, double q0)
339 {
340 double dh = h(mass, q0) * (1.0 / (8 * q0 * q0) - 1.0 / (2 * mass * mass)) + 1.0 / (2 * pi * mass * mass);
341 return dh;
342 }
343
344 double EvtD0TopipiEta::f(double mass, double sx, double q0, double q)
345 {
346 double m = sqrt(sx);
347 double f = mass * mass / (pow(q0, 3)) * (q * q * (h(m, q) - h(mass, q0)) + (mass * mass - sx) * q0 * q0 * dh(mass, q0));
348 return f;
349 }
350
351 double EvtD0TopipiEta::d(double mass, double q0)
352 {
353 double d = 3.0 / pi * spi / (q0 * q0) * log((mass + 2 * q0) / (2 * mpi)) + mass / (2 * pi * q0) - (spi * mass) / (pi * pow(q0, 3));
354 return d;
355 }
356
357 EvtComplex EvtD0TopipiEta::propagatorGS(double mass, double width, double sa, double sb, double sc, double r, int l)
358 {
359 double q = Qabcs(sa, sb, sc);
360 double sa0 = mass * mass;
361 double q0 = Qabcs(sa0, sb, sc);
362 q = sqrt(q);
363 q0 = sqrt(q0);
364 EvtComplex prop = (1 + d(mass, q0) * width / mass) / (mass * mass - sa + width * f(mass, sa, q0, q) - ci * mass * width * wid(mass,
365 sa, sb, sc, r, l));
366 return prop;
367 }
368
369 EvtComplex EvtD0TopipiEta::Flatte_rhoab(double sa, double sb, double sc)
370 {
371 double q = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
372 EvtComplex rhoo;
373 if (q > 0) {
374 rhoo = 2.0 * sqrt(q / sa) * one;
375 }
376 if (q < 0) {
377 rhoo = 2.0 * sqrt(-q / sa) * ci;
378 }
379 return rhoo;
380 }
381
382 EvtComplex EvtD0TopipiEta::propagatorFlatte(double mass, double width, double sx, double* sb, double* sc)
383 {
384 (void)width;
385 const double g1sq = 0.5468 * 0.5468;
386 const double g2sq = 0.23 * 0.23;
387 EvtComplex rho1 = Flatte_rhoab(sx, sb[0], sc[0]);
388 EvtComplex rho2 = Flatte_rhoab(sx, sb[1], sc[1]);
389 EvtComplex prop = 1.0 / (mass * mass - sx - ci * (g1sq * rho1 + g2sq * rho2));
390 return prop;
391 }
393} // Belle 2 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