Belle II Software development
EvtDsToKKpi.cc
1// Model: EvtDsToKKpi
2// This file is an amplitude model for Ds+ -> K- K+ pi+.
3// The model is from the BESIII Collaboration in PRD 104, 012016 (2021). DOI: https://doi.org/10.1103/PhysRevD.104.012016
4//
5// Permission to use 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 <iomanip>
10#include <cmath>
11#include <string>
12#include <EvtGenBase/EvtCPUtil.hh>
13#include <EvtGenBase/EvtTensor4C.hh>
14
15#include <EvtGenBase/EvtPatches.hh>
16#include <fstream>
17#include <stdlib.h>
18#include <EvtGenBase/EvtParticle.hh>
19#include <EvtGenBase/EvtGenKine.hh>
20#include <EvtGenBase/EvtPDL.hh>
21#include <EvtGenBase/EvtReport.hh>
22#include <EvtGenBase/EvtResonance.hh>
23#include <EvtGenBase/EvtResonance2.hh>
24#include <string>
25#include <EvtGenBase/EvtConst.hh>
26#include <EvtGenBase/EvtFlatte.hh>
27#include <EvtGenBase/EvtDecayTable.hh>
28
29#include <generators/evtgen/EvtGenModelRegister.h>
30#include <generators/evtgen/models/besiii/EvtDsToKKpi.h>
31
32namespace Belle2 {
37
39 B2_EVTGEN_REGISTER_MODEL(EvtDsToKKpi);
40
41 EvtDsToKKpi::~EvtDsToKKpi() {}
42
43 std::string EvtDsToKKpi::getName()
44 {
45 return "DsToKKpi";
46 }
47
48 EvtDecayBase* EvtDsToKKpi::clone()
49 {
50 return new EvtDsToKKpi;
51 }
52
53 void EvtDsToKKpi::init()
54 {
55 checkNArg(0);
56 checkNDaug(3);
57 checkSpinParent(EvtSpinType::SCALAR);
58 checkSpinDaughter(0, EvtSpinType::SCALAR);
59 checkSpinDaughter(1, EvtSpinType::SCALAR);
60 checkSpinDaughter(2, EvtSpinType::SCALAR);
61
62
63 phi[0] = 0;
64 phi[1] = 6.1944e+00;
65 phi[2] = 4.7398e+00;
66 phi[3] = 2.9047e+00;
67 phi[4] = 1.0068e+00;
68 phi[5] = 5.8035e-01;
69
70 rho[0] = 1;
71 rho[1] = 1.0963e+00;
72 rho[2] = 2.7818e+00;
73 rho[3] = 1.2570e+00;
74 rho[4] = 7.7351e-01;
75 rho[5] = 5.6277e-01;
76
77 mass[0] = 8.9581e-01;
78 mass[1] = 1.019461e+00;
79 mass[2] = 0.919;
80 mass[3] = 1.4712e+00;
81 mass[4] = 1.7220e+00;
82 mass[5] = 1.3500e+00;
83
84 width[0] = 4.7400e-02;
85 width[1] = 4.2660e-03;
86 width[2] = 0.272;
87 width[3] = 2.7000e-01;
88 width[4] = 1.3500e-01;
89 width[5] = 2.6500e-01;
90
91 modetype[0] = 0;
92 modetype[1] = 1;
93 modetype[2] = 13;
94 modetype[3] = 3;
95 modetype[4] = 4;
96 modetype[5] = 5;
97
98
99 mD = 1.86484;
100 mDs = 1.96835;
101 rRes = 1.5;
102 rD = 3.0;
103 mkstr = 0.89555;
104 mk0 = 0.497611;
105 mass_Kaon = 0.49368;
106 mass_Pion = 0.13957;
107 mass_Pi0 = 0.1349768;
108 mass_EtaP = 0.95778;
109 mass_Eta = 0.547862;
110 math_pi = 3.1415926;
111 afRatio = 2.04835; //BABAR // afRatio= 1.23202; //BES2
112 int GG[4][4] = { {1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1} };
113 for (int i = 0; i < 4; i++) {
114 for (int j = 0; j < 4; j++) {
115 G[i][j] = GG[i][j];
116 }
117 }
118 }
119
120 void EvtDsToKKpi::initProbMax()
121 {
122 setProbMax(35638.0);
123 }
124
125 void EvtDsToKKpi::decay(EvtParticle* p)
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 double P1[4], P2[4], P3[4];
132 P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
133 P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
134 P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
135
136 double value;
137 int g0[6] = {1, 1, 1, 1, 1, 1};
138 int nstates = 6;
139 calEvaMy(P1, P2, P3, mass, width, rho, phi, g0, modetype, nstates, value);
140 setProb(value);
141 return ;
142 }
143
144 void EvtDsToKKpi::MIP_LineShape(double sa, double pro[2])
145 {
146 double m0 = mass[2], cKK = width[2];
147 pro[0] = sqrt(1 / (((m0 * m0) - sa) * ((m0 * m0) - sa) + (cKK * m0 * sqrt(1 - 4 * mass_Kaon * mass_Kaon / sa)) * (cKK * m0 * sqrt(
148 1 - 4 * mass_Kaon * mass_Kaon / sa))));
149 pro[1] = 0;
150 }
151 void EvtDsToKKpi::calEvaMy(double* pKm, double* pKp, double* pPi, double* mass1, double* width1, double* amp, double* phase,
152 int* g0, int* modetype_param, int nstates, double& Result) // Renamed modetype to modetype_param to avoid shadowing member variable
153 {
154 double pF0_980[4], pPhi1020[4], pKstr[4], pD[4], p23[4];
155 pKp[0] = sqrt(mass_Kaon * mass_Kaon + pKp[1] * pKp[1] + pKp[2] * pKp[2] + pKp[3] * pKp[3]);
156 pKm[0] = sqrt(mass_Kaon * mass_Kaon + pKm[1] * pKm[1] + pKm[2] * pKm[2] + pKm[3] * pKm[3]);
157 pPi[0] = sqrt(mass_Pion * mass_Pion + pPi[1] * pPi[1] + pPi[2] * pPi[2] + pPi[3] * pPi[3]);
158 for (int u = 0; u < 4; u++) {
159 pF0_980[u] = pKp[u] + pKm[u];
160 pPhi1020[u] = pKp[u] + pKm[u];
161 pKstr[u] = pKm[u] + pPi[u];
162 p23[u] = pKp[u] + pPi[u];
163 pD[u] = pKp[u] + pKm[u] + pPi[u];
164 }
165 double cof[2], amp_tmp1[2], amp_tmp[2], amp_PDF[2], PDF[2];
166 amp_PDF[0] = 0;
167 amp_PDF[1] = 0;
168 PDF[0] = 0;
169 PDF[1] = 0;
170 double temp_PDF, tmp1;
171 double pro[2], B[3];
172 double t1kstr1[4], t1phi1020[4], t1D1[4], t1D2[4];
173
174 double sD, sF0_980, sF0_1710, sF0_1370, sKp, sKm, sPi, sKstr, sPhi1020, sKstr1430;
175 sF0_980 = SCADot(pF0_980, pF0_980);
176 sF0_1710 = sF0_980;
177 sF0_1370 = sF0_980;
178 sKstr = SCADot(pKstr, pKstr);
179 sD = SCADot(pD, pD);
180 sKstr1430 = sKstr;
181 sPhi1020 = SCADot(pPhi1020, pPhi1020);
182 sKp = SCADot(pKp, pKp);
183 sKm = SCADot(pKm, pKm);
184 sPi = SCADot(pPi, pPi);
185
186 calt1(pKm, pPi, t1kstr1);
187 calt1(pKm, pKp, t1phi1020);
188 calt1(pKstr, pKp, t1D1);
189 calt1(pPhi1020, pPi, t1D2);
190
191 for (int i = 0; i < nstates; i++) {
192 amp_tmp[0] = 0;
193 amp_tmp[1] = 0;
194 amp_tmp1[0] = 0;
195 amp_tmp1[1] = 0;
196 tmp1 = 0;
197 temp_PDF = 0;
198 cof[0] = amp[i] * cos(phase[i]);
199 cof[1] = amp[i] * sin(phase[i]);
200 if (modetype_param[i] == 13) {
201 //a0(980) and f0(980) mixture
202 double amp_a0[2] = {0};
203 double sa0_980 = sF0_980;
204 MIP_LineShape(sa0_980, pro);
205 B[0] = barrier(0, sa0_980, sKp, sKm, rRes);
206 temp_PDF = 1;
207 tmp1 = temp_PDF * B[0];
208 amp_a0[0] = tmp1 * pro[0];
209 amp_a0[1] = tmp1 * pro[1];
210 amp_tmp1[0] = amp_a0[0] ;
211 amp_tmp1[1] = amp_a0[1] ;
212 } else if (modetype_param[i] == 0) {
213 //K*(892) K+
214 if (g0[i] == 0) {
215 pro[0] = 1;
216 pro[1] = 0;
217 }
218 bool neo = true;
219 if (neo) {
220 double sBC = SCADot(p23, p23);
221 if (g0[i] == 1) propagatorRBWNeoKstr892(mass1[i], width1[i], sKstr, sPi, sKm, rRes, 1, pro);
222 B[0] = barrierNeo(1, sKstr, sPi, sKm, rRes, mass1[i]);
223 B[1] = barrierNeoDs(1, sD, sKstr, sKp, rD, mDs, mass1[i]);
224 temp_PDF = (sBC - sPhi1020 + ((sD - sKp) * (sKm - sPi) / (sKstr)));
225 } else {
226 if (g0[i] == 1) propagatorRBW(mass1[i], width1[i], sKstr, sKm, sPi, rRes, 1, pro);
227 B[0] = barrier(1, sKstr, sKm, sPi, rRes);
228 B[1] = barrier(1, sD, sKstr, sKp, rD);
229 for (int m = 0; m < 4; m++) {
230 for (int j = 0; j < 4; j++) {
231 temp_PDF += t1D1[m] * t1kstr1[j] * G[m][j];
232 }
233 }
234 }
235 tmp1 = temp_PDF * B[0] * B[1];
236 amp_tmp1[0] = tmp1 * pro[0];
237 amp_tmp1[1] = tmp1 * pro[1];
238 } else if (modetype_param[i] == 1) {
239 if (g0[i] == 0) {
240 pro[0] = 1;
241 pro[1] = 0;
242 }
243 bool neo = true;
244 if (neo) {
245 if (g0[i] == 1) propagatorRBWNeo(mass1[i], width1[i], sPhi1020, sKm, sKp, rRes, 1, pro);
246 B[0] = barrierNeo(1, sPhi1020, sKp, sKm, rRes, mass1[i]);
247 B[1] = barrierNeoDs(1, sD, sPhi1020, sPi, rD, mDs, mass1[i]);
248 double sBC = SCADot(p23, p23);
249 temp_PDF = (sKstr - sBC + ((sD - sPi) * (sKp - sKm) / (sKstr)));
250 } else {
251 if (g0[i] == 1) propagatorRBW(mass1[i], width1[i], sPhi1020, sKm, sKp, rRes, 1, pro);
252 B[0] = barrier(1, sPhi1020, sKp, sKm, rRes);
253 B[1] = barrier(1, sD, sPhi1020, sPi, rD);
254 for (int m = 0; m < 4; m++) {
255 for (int j = 0; j < 4; j++) {
256 temp_PDF += t1D2[m] * t1phi1020[j] * G[m][j];
257 }
258 }
259 }
260
261 tmp1 = temp_PDF * B[0] * B[1];
262 amp_tmp1[0] = tmp1 * pro[0];
263 amp_tmp1[1] = tmp1 * pro[1];
264 } else if (modetype_param[i] == 3) {
265 //Kstr1430 K S
266 double sKm2[2] = {sKm, mass_EtaP * mass_EtaP};
267 double sPi2[2] = {sPi, mass_Kaon * mass_Kaon};
268 propagatorKstr1430(mass1[i], sKstr1430, sKm2, sPi2, pro);
269 B[0] = barrier(0, sPhi1020, sKp, sKm, rRes);
270 tmp1 = 1 * B[0];
271 amp_tmp1[0] = tmp1 * pro[0];
272 amp_tmp1[1] = tmp1 * pro[1];
273
274 } else if (modetype_param[i] == 4) {
275 if (g0[i] == 1) propagatorRBWNeo(mass1[i], width1[i], sF0_1710, sKp, sKm, rRes, 0, pro);
276 if (g0[i] == 0) {
277 pro[0] = 1;
278 pro[1] = 0;
279 }
280 B[0] = barrier(0, sF0_1710, sKp, sKm, rRes);
281 temp_PDF = 1;
282 tmp1 = temp_PDF * B[0];
283 amp_tmp1[0] = tmp1 * pro[0];
284 amp_tmp1[1] = tmp1 * pro[1];
285 } else if (modetype_param[i] == 5) {
286 //f0(1370) Pi+ S
287 if (g0[i] == 1) propagatorRBWNeo(mass1[i], width1[i], sF0_1370, sKp, sKm, rRes, 0, pro);
288 if (g0[i] == 0) {
289 pro[0] = 1;
290 pro[1] = 0;
291 }
292 B[0] = barrier(0, sF0_1370, sKp, sKm, rRes);
293 tmp1 = 1 * B[0];
294 amp_tmp1[0] = tmp1 * pro[0];
295 amp_tmp1[1] = tmp1 * pro[1];
296 }
297 amp_tmp[0] = amp_tmp1[0];
298 amp_tmp[1] = amp_tmp1[1];
299 Com_Multi(amp_tmp, cof, amp_PDF);
300 PDF[0] += amp_PDF[0];
301 PDF[1] += amp_PDF[1];
302 }
303
304 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
305 if (value <= 0) {
306 value = 1e-20;
307 }
308 Result = value;
309 }
310
311 void EvtDsToKKpi::Com_Multi(double a1[2], double a2[2], double res[2])
312 {
313 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
314 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
315 }
316 void EvtDsToKKpi::Com_Divide(double a1[2], double a2[2], double res[2])
317 {
318 res[0] = (a1[0] * a2[0] + a1[1] * a2[1]) / (a2[0] * a2[0] + a2[1] * a2[1]);
319 res[1] = (a1[1] * a2[0] - a1[0] * a2[1]) / (a2[0] * a2[0] + a2[1] * a2[1]);
320 }
321
322 double EvtDsToKKpi::SCADot(double a1[4], double a2[4])
323 {
324 double _cal = 0.;
325 _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
326 return _cal;
327 }
328 double EvtDsToKKpi::barrier(int l, double sa, double sb, double sc, double r)
329 {
330 double q = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
331 if (q < 0) q = 1e-16;
332 double z;
333 z = q * r * r;
334 double F = 0.0;
335 if (l == 0) F = 1;
336 if (l == 1) F = sqrt(2 * z / (1 + z));
337 if (l == 2) F = sqrt(13 * z * z / (9 + 3 * z + z * z));
338 return F;
339 }
340 double EvtDsToKKpi::barrierNeo(int l, double sa, double sb, double sc, double r, double mR)
341 {
342 double pAB = ((sa - sb - sc) * (sa - sb - sc) / 4.0 - (sb * sc)) / sa;
343 double pR = ((mR * mR - sb - sc) * (mR * mR - sb - sc) / 4.0 - (sb * sc)) / (mR * mR);
344 double zAB = pAB * r * r;
345 double zR = pR * r * r;
346 double F = 0.0;
347 if (l == 0) F = 1;
348 if (l == 1) F = sqrt((1 + zR) / (1 + zAB));
349 if (l == 2) F = sqrt((9 + 3 * zAB + zAB * zAB) / (9 + 3 * zAB + zAB * zAB));
350 return F;
351 }
352 double EvtDsToKKpi::barrierNeoDs(int l, double sa, double sb, double sc, double r, double mR, double mb)
353 {
354 double pAB = ((sa - sb - sc) * (sa - sb - sc) / 4.0 - (sb * sc)) / sa;
355 double pR = ((sa - mb * mb - sc) * (sa - mb * mb - sc) / 4.0 - (mb * mb * sc)) / (mR * mR);
356 double zAB = pAB * r * r;
357 double zR = pR * r * r;
358 double F = 0.0;
359 if (l == 0) F = 1;
360 if (l == 1) F = sqrt((1 + zR) / (1 + zAB));
361 if (l == 2) F = sqrt((9 + 3 * zAB + zAB * zAB) / (9 + 3 * zAB + zAB * zAB));
362 return F;
363 }
364//------------------spin-------------------------------------------
365 void EvtDsToKKpi::calt1(double daug1[4], double daug2[4], double t1[4])
366 {
367 double p, pq;
368 double pa[4], qa[4];
369 for (int i = 0; i < 4; i++) {
370 pa[i] = daug1[i] + daug2[i];
371 qa[i] = daug1[i] - daug2[i];
372 }
373 p = SCADot(pa, pa);
374 pq = SCADot(pa, qa);
375 for (int i = 0; i < 4; i++) {
376 t1[i] = qa[i] - pq / p * pa[i];
377 }
378 }
379 void EvtDsToKKpi::calt2(double daug1[4], double daug2[4], double t2[4][4])
380 {
381 double p, r;
382 double pa[4], t1[4];
383 calt1(daug1, daug2, t1);
384 r = SCADot(t1, t1);
385 for (int i = 0; i < 4; i++) {
386 pa[i] = daug1[i] + daug2[i];
387 }
388 p = SCADot(pa, pa);
389 for (int i = 0; i < 4; i++) {
390 for (int j = 0; j < 4; j++) {
391 t2[i][j] = t1[i] * t1[j] - 1.0 / 3 * r * (G[i][j] - pa[i] * pa[j] / p);
392 }
393 }
394 }
395//-------------------prop--------------------------------------------
396 void EvtDsToKKpi::propagator(double mass_param, double width_param, double sx, double prop[2])
397 {
398 double a[2], b[2];
399 a[0] = 1;
400 a[1] = 0;
401 b[0] = mass_param * mass_param - sx;
402 b[1] = -mass_param * width_param;
403 Com_Divide(a, b, prop);
404 }
405 double EvtDsToKKpi::wid(double mass_param, double sa, double sb, double sc, double r, int l)
406 {
407 double widm = 0.;
408 double q = 0.;
409 double q0 = 0.;
410 double sa0 = mass_param * mass_param;
411 double m = sqrt(sa);
412 q = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
413 if (q < 0) q = 1e-16;
414 q0 = (sa0 + sb - sc) * (sa0 + sb - sc) / (4 * sa0) - sb;
415 if (q0 < 0) q0 = 1e-16;
416 double z = q * r * r;
417 double z0 = q0 * r * r;
418 double F = 0;
419 if (l == 0) F = 1;
420 if (l == 1) F = sqrt((1 + z0) / (1 + z));
421 if (l == 2) F = sqrt((9 + 3 * z0 + z0 * z0) / (9 + 3 * z + z * z));
422 double t = q / q0;
423 widm = pow(t, l + 0.5) * mass_param / m * F * F;
424 return widm;
425 }
426
427 void EvtDsToKKpi::Flatte_rhoab(double sa, double sb, double sc, double rho_param[2])
428 {
429 double q = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
430 if (q > 0) {
431 rho_param[0] = 2 * sqrt(q / sa);
432 rho_param[1] = 0;
433 } else if (q < 0) {
434 rho_param[0] = 0;
435 rho_param[1] = 2 * sqrt(-q / sa);
436 }
437 }
438
439 void EvtDsToKKpi::propagatorFlatte(double mass_param, double width_param __attribute__((unused)), double sx, double* sb, double* sc,
440 double prop[2])
441 {
442 double unit[2] = {1.0};
443 double ci[2] = {0, 1};
444 double rho1[2];
445 Flatte_rhoab(sx, sb[0], sc[0], rho1);
446 double rho2[2];
447 Flatte_rhoab(sx, sb[1], sc[1], rho2);
448 double g1_f980 = 0.165, g2_f980 = 0.69465;
449 double tmp1[2] = {g1_f980, 0};
450 double tmp11[2];
451 double tmp2[2] = {g2_f980, 0};
452 double tmp22[2];
453 Com_Multi(tmp1, rho1, tmp11);
454 Com_Multi(tmp2, rho2, tmp22);
455 double tmp3[2] = {tmp11[0] + tmp22[0], tmp11[1] + tmp22[1]};
456 double tmp31[2];
457 Com_Multi(tmp3, ci, tmp31);
458 double tmp4[2] = {mass_param* mass_param - sx - tmp31[0], -1.0 * tmp3[1]};
459 Com_Divide(unit, tmp4, prop);
460 }
461 void EvtDsToKKpi::propagator980(double mass_param, double sx, double* sb, double* sc, double prop[2])
462 {
463 double unit[2] = {1.0};
464 double ci[2] = {0, 1};
465 double rho1[2];
466 Flatte_rhoab(sx, sb[0], sc[0], rho1);
467 double rho2[2];
468 Flatte_rhoab(sx, sb[1], sc[1], rho2);
469 double gK_f980 = 0.69466, gPi_f980 = 0.165;
470 double tmp1[2] = {gK_f980, 0};
471 double tmp11[2];
472 double tmp2[2] = {gPi_f980, 0};
473 double tmp22[2];
474 Com_Multi(tmp1, rho1, tmp11);
475 Com_Multi(tmp2, rho2, tmp22);
476 double tmp3[2] = {tmp11[0] + tmp22[0], tmp11[1] + tmp22[1]};
477 double tmp31[2];
478 Com_Multi(tmp3, ci, tmp31);
479 double tmp4[2] = {mass_param* mass_param - sx - tmp31[0], -1.0 * tmp31[1]};
480 Com_Divide(unit, tmp4, prop);
481 }
482
483 void EvtDsToKKpi::propagatora0980(double mass_param, double sx, double* sb, double* sc, double prop[2])
484 {
485 double unit[2] = {1.0};
486 double ci[2] = {0, 1};
487 double rho1[2];
488 Flatte_rhoab(sx, sb[0], sc[0], rho1);
489 double rho2[2];
490 Flatte_rhoab(sx, sb[1], sc[1], rho2);
491 double gKK_a980 = 0.892 * 0.341, gPiEta_a980 = 0.341;
492 double tmp1[2] = {gKK_a980, 0};
493 double tmp11[2];
494 double tmp2[2] = {gPiEta_a980, 0};
495 double tmp22[2];
496 Com_Multi(tmp1, rho1, tmp11);
497 Com_Multi(tmp2, rho2, tmp22);
498 double tmp3[2] = {tmp11[0] + tmp22[0], tmp11[1] + tmp22[1]};
499 double tmp31[2];
500 Com_Multi(tmp3, ci, tmp31);
501 double tmp4[2] = {mass_param* mass_param - sx - tmp31[0], -1.0 * tmp31[1]};
502 Com_Divide(unit, tmp4, prop);
503 }
504
505 void EvtDsToKKpi::propagatorKstr1430(double mass_param, double sx, double* sb, double* sc, double prop[2])
506 {
507 double unit[2] = {1.0};
508 double ci[2] = {0, 1};
509 double rho1[2];
510 Flatte_rhoab(sx, sb[0], sc[0], rho1);
511 double rho2[2];
512 Flatte_rhoab(sx, sb[1], sc[1], rho2);
513 double gKPi_Kstr1430 = 0.2990, gEtaPK_Kstr1430 = 0.0529;
514 double tmp1[2] = {gKPi_Kstr1430, 0};
515 double tmp11[2];
516 double tmp2[2] = {gEtaPK_Kstr1430, 0};
517 double tmp22[2];
518 Com_Multi(tmp1, rho1, tmp11);
519 Com_Multi(tmp2, rho2, tmp22);
520 double tmp3[2] = {tmp11[0] + tmp22[0], tmp11[1] + tmp22[1]};
521 double tmp31[2];
522 Com_Multi(tmp3, ci, tmp31);
523 double tmp4[2] = {mass_param* mass_param - sx - tmp31[0], -1.0 * tmp31[1]};
524 Com_Divide(unit, tmp4, prop);
525 }
526
527 void EvtDsToKKpi::propagatorRBW(double mass_param, double width_param, double sa, double sb, double sc, double r, int l,
528 double prop[2])
529 {
530 double a[2], b[2];
531 a[0] = 1;
532 a[1] = 0;
533 b[0] = mass_param * mass_param - sa;
534 b[1] = -mass_param * width_param * wid(mass_param, sa, sb, sc, r, l);
535 Com_Divide(a, b, prop);
536 }
537
538 void EvtDsToKKpi::propagatorRBWNeo(double mass_param, double width_param, double sa, double sb, double sc,
539 double r __attribute__((unused)), int l, double prop[2])
540 {
541 double a[2], b[2];
542 a[0] = 1;
543 a[1] = 0;
544 b[0] = mass_param * mass_param - sa;
545 double tmp1 = (sa - sb - sc);
546 double tmp2 = sb * sc;
547 double pAB = sqrt((tmp1 * tmp1 / 4.0 - tmp2) / sa);
548 double pR = sqrt(((mass_param * mass_param - sb - sc) * (mass_param * mass_param - sb - sc) / 4.0 - (sb * sc)) /
549 (mass_param * mass_param));
550 double fR = sqrt(1.0 + 1.5 * 1.5 * pR * pR) / sqrt(1.0 + 1.5 * 1.5 * pAB * pAB);
551 double power = 1;
552 if (!l) {
553 power = 1;
554 fR = 1;
555 } else if (l == 1) {
556 power = 3;
557 }
558 double gammaAB = width_param * pow(pAB / pR, power) * (mass_param / sqrt(sa)) * fR * fR;
559 b[1] = -mass_param * gammaAB;
560 Com_Divide(a, b, prop);
561 }
562
563 void EvtDsToKKpi::propagatorRBWNeoKstr892(double mass_param, double width_param, double sa, double sb, double sc,
564 double r __attribute__((unused)), int l, double prop[2])
565 {
566 double a[2], b[2];
567 a[0] = 1;
568 a[1] = 0;
569 b[0] = mass_param * mass_param - sa;
570 double tmp1 = (sa - sb - sc);
571 double tmp2 = sb * sc;
572 double pAB = sqrt((tmp1 * tmp1 / 4.0 - tmp2) / sa);
573 double pR = sqrt(((mass_param * mass_param - sb - sc) * (mass_param * mass_param - sb - sc) / 4.0 - (sb * sc)) /
574 (mass_param * mass_param));
575 double fR = sqrt(1.0 + 1.5 * 1.5 * pR * pR) / sqrt(1.0 + 1.5 * 1.5 * pAB * pAB);
576 double power = 1;
577 if (!l) {
578 power = 1;
579 } else if (l == 1) {
580 power = 3;
581 }
582 double gammaAB = width_param * pow(pAB / pR, power) * (mass_param / sqrt(sa)) * fR * fR;
583 b[1] = -mass_param * gammaAB;
584 Com_Divide(a, b, prop);
585 }
586
587//------------GS---used by rho----------------------------
588 double EvtDsToKKpi::h(double m, double q)
589 {
590 double h = 2 / math_pi * q / m * log((m + 2 * q) / (2 * mass_Pion));
591 return h;
592 }
593 double EvtDsToKKpi::dh(double mass_param, double q0)
594 {
595 double dh = h(mass_param, q0) * (1.0 / (8 * q0 * q0) - 1.0 / (2 * mass_param * mass_param)) + 1.0 /
596 (2 * math_pi * mass_param * mass_param);
597 return dh;
598 }
599 double EvtDsToKKpi::f(double mass_param, double sx, double q0, double q)
600 {
601 double m = sqrt(sx);
602 double f = mass_param * mass_param / (pow(q0, 3)) * (q * q * (h(m, q) - h(mass_param,
603 q0)) + (mass_param * mass_param - sx) * q0 * q0 * dh(mass_param, q0));
604 return f;
605 }
606 double EvtDsToKKpi::d(double mass_param, double q0)
607 {
608 double d = 3.0 / math_pi * mass_Pion * mass_Pion / (q0 * q0) * log((mass_param + 2 * q0) / (2 * mass_Pion)) + mass_param /
609 (2 * math_pi * q0)
610 - (mass_Pion * mass_Pion * mass_param) / (math_pi * pow(q0, 3));
611 return d;
612 }
613
614//rho
615 void EvtDsToKKpi::propagatorGS(double mass_param, double width_param, double sa, double sb, double sc, double r, int l,
616 double prop[2])
617 {
618 double a[2], b[2];
619 double q = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
620 double sa0 = mass_param * mass_param;
621 double q0 = (sa0 + sb - sc) * (sa0 + sb - sc) / (4 * sa0) - sb;
622 if (q < 0) q = 1e-16;
623 if (q0 < 0) q0 = 1e-16;
624 q = sqrt(q);
625 q0 = sqrt(q0);
626 a[0] = 1 + d(mass_param, q0) * width_param / mass_param;
627 a[1] = 0;
628 b[0] = mass_param * mass_param - sa + width_param * f(mass_param, sa, q0, q);
629 b[1] = -mass_param * width_param * wid(mass_param, sa, sb, sc, r, l);
630 Com_Divide(a, b, prop);
631 }
632
633
634
636} // 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