12#include <EvtGenBase/EvtCPUtil.hh>
13#include <EvtGenBase/EvtTensor4C.hh>
14#include <EvtGenBase/EvtPatches.hh>
17#include <EvtGenBase/EvtParticle.hh>
18#include <EvtGenBase/EvtGenKine.hh>
19#include <EvtGenBase/EvtPDL.hh>
20#include <EvtGenBase/EvtReport.hh>
21#include <EvtGenBase/EvtResonance.hh>
22#include <EvtGenBase/EvtResonance2.hh>
24#include <EvtGenBase/EvtConst.hh>
25#include <EvtGenBase/EvtDecayTable.hh>
27#include <generators/evtgen/EvtGenModelRegister.h>
28#include <generators/evtgen/models/besiii/EvtDsToKpipi.h>
39 EvtDsToKpipi::~EvtDsToKpipi() {}
41 std::string EvtDsToKpipi::getName()
46 EvtDecayBase* EvtDsToKpipi::clone()
48 return new EvtDsToKpipi;
51 void EvtDsToKpipi::init()
55 checkSpinParent(EvtSpinType::SCALAR);
56 checkSpinDaughter(0, EvtSpinType::SCALAR);
57 checkSpinDaughter(1, EvtSpinType::SCALAR);
58 checkSpinDaughter(2, EvtSpinType::SCALAR);
78 phi[2] = -1.249864467;
79 phi[3] = -0.3067504308;
80 phi[4] = 0.9229242379;
81 phi[5] = -3.278567926;
82 phi[6] = -0.6346408622;
86 rho[2] = 0.7361782665;
118 mass[7] = 1.432787726;
129 mass_Pi0 = 0.1349766;
144 int GG[4][4] = { {1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1} };
145 for (
int i = 0; i < 4; i++) {
146 for (
int j = 0; j < 4; j++) {
152 void EvtDsToKpipi::initProbMax()
157 void EvtDsToKpipi::decay(EvtParticle* p)
159 p->initializePhaseSpace(getNDaug(), getDaugs());
160 EvtVector4R D1 = p->getDaug(0)->getP4();
161 EvtVector4R D2 = p->getDaug(1)->getP4();
162 EvtVector4R D3 = p->getDaug(2)->getP4();
164 double P1[4], P2[4], P3[4];
165 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
166 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
167 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
170 int g0[8] = {1, 1, 4, 2, 500, 1, 1, 2};
171 int spin_param[8] = {1, 1, 0, 0, 0, 1, 1, 0};
173 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin_param, modetype, nstates, value);
180 void EvtDsToKpipi::Com_Multi(
double a1[2],
double a2[2],
double res[2])
182 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
183 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
185 void EvtDsToKpipi::Com_Divide(
double a1[2],
double a2[2],
double res[2])
187 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
188 res[0] = (a1[0] * a2[0] + a1[1] * a2[1]) / tmp;
189 res[1] = (a1[1] * a2[0] - a1[0] * a2[1]) / tmp;
192 double EvtDsToKpipi::SCADot(
double a1[4],
double a2[4])
194 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
197 double EvtDsToKpipi::barrier(
int l,
double sa,
double sb,
double sc,
double r,
double mass_param)
199 double q = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
200 if (q < 0) q = 1e-16;
204 sa0 = mass_param * mass_param;
205 double q0 = (sa0 + sb - sc) * (sa0 + sb - sc) / (4 * sa0) - sb;
206 if (q0 < 0) q0 = 1e-16;
207 double z0 = q0 * r * r;
210 if (l == 1) F =
sqrt((1 + z0) / (1 + z));
211 if (l == 2) F =
sqrt((9 + 3 * z0 + z0 * z0) / (9 + 3 * z + z * z));
214 void EvtDsToKpipi::calt1(
double daug1[4],
double daug2[4],
double t1[4])
218 for (
int i = 0; i < 4; i++) {
219 pa[i] = daug1[i] + daug2[i];
220 qa[i] = daug1[i] - daug2[i];
225 for (
int i = 0; i < 4; i++) {
226 t1[i] = qa[i] - tmp * pa[i];
229 void EvtDsToKpipi::calt2(
double daug1[4],
double daug2[4],
double t2[4][4])
233 calt1(daug1, daug2, t1);
234 r = SCADot(t1, t1) / 3.0;
235 for (
int i = 0; i < 4; i++) {
236 pa[i] = daug1[i] + daug2[i];
239 for (
int i = 0; i < 4; i++) {
240 for (
int j = 0; j < 4; j++) {
241 t2[i][j] = t1[i] * t1[j] - r * (G[i][j] - pa[i] * pa[j] / p);
248 void EvtDsToKpipi::propagatorCBW(
double mass_param,
double width_param,
double sx,
double prop[2])
253 b[0] = mass_param * mass_param - sx;
254 b[1] = -mass_param * width_param;
255 Com_Divide(a, b, prop);
257 double EvtDsToKpipi::wid(
double mass2,
double mass_param,
double sa,
double sb,
double sc,
double r2,
int l)
261 double tmp = sb - sc;
262 double tmp1 = sa + tmp;
263 double q = 0.25 * tmp1 * tmp1 / sa - sb;
264 if (q < 0) q = 1e-16;
265 double tmp2 = mass2 + tmp;
266 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
267 if (q0 < 0) q0 = 1e-16;
271 if (l == 0) {widm =
sqrt(t) * mass_param / m;}
272 else if (l == 1) {widm = t *
sqrt(t) * mass_param / m * (1 + z0) / (1 + z);}
273 else if (l == 2) {widm = t * t *
sqrt(t) * mass_param / m * (9 + 3 * z0 + z0 * z0) / (9 + 3 * z + z * z);}
277 double EvtDsToKpipi::widl1(
double mass2,
double mass_param,
double sa,
double sb,
double sc,
double r2)
281 double tmp = sb - sc;
282 double tmp1 = sa + tmp;
283 double q = 0.25 * tmp1 * tmp1 / sa - sb;
284 if (q < 0) q = 1e-16;
285 double tmp2 = mass2 + tmp;
286 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
287 if (q0 < 0) q0 = 1e-16;
290 double F = (1 + z0) / (1 + z);
292 widm = t *
sqrt(t) * mass_param / m * F;
295 void EvtDsToKpipi::propagatorRBW(
double mass_param,
double width_param,
double sa,
double sb,
double sc,
double r2,
int l,
299 double mass2 = mass_param * mass_param;
304 b[1] = -mass_param * width_param * wid(mass2, mass_param, sa, sb, sc, r2, l);
305 Com_Divide(a, b, prop);
307 void EvtDsToKpipi::propagatorFlatte(
double mass_param,
double width_param __attribute__((unused)),
double sa,
double prop[2])
311 double rhoPi[2] = {0, 0}, rhoKa[2] = {0, 0};
313 q2_Pi = 0.25 * sa - mPi * mPi;
314 q2_Ka = 0.25 * sa - mKa * mKa;
317 rhoPi[0] = 2.0 *
sqrt(q2_Pi / sa);
322 rhoPi[1] = 2.0 *
sqrt(-q2_Pi / sa);
326 rhoKa[0] = 2.0 *
sqrt(q2_Ka / sa);
331 rhoKa[1] = 2.0 *
sqrt(-q2_Ka / sa);
337 b[0] = mass_param * mass_param - sa + 0.165 * rhoPi[1] + 0.69465 * rhoKa[1];
338 b[1] = - (0.165 * rhoPi[0] + 0.69465 * rhoKa[0]);
339 Com_Divide(a, b, prop);
341 void EvtDsToKpipi::propagatorGS(
double mass_param,
double width_param,
double sa,
double sb,
double sc,
double r2,
double prop[2])
344 double mass2 = mass_param * mass_param;
345 double tmp = sb - sc;
346 double tmp1 = sa + tmp;
347 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
348 if (q2 < 0) q2 = 1e-16;
350 double tmp2 = mass2 + tmp;
351 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
352 if (q02 < 0) q02 = 1e-16;
355 double q0 =
sqrt(q02);
357 double q03 = q0 * q02;
358 double tmp3 = log(mass_param + 2 * q0) + 1.2760418309;
360 double h = GS1 * q / m * (log(m + 2 * q) + 1.2760418309);
361 double h0 = GS1 * q0 / mass_param * tmp3;
362 double dh = h0 * (0.125 / q02 - 0.5 / mass2) + GS3 / mass2;
363 double d = GS2 / q02 * tmp3 + GS3 * mass_param / q0 - GS4 * mass_param / q03;
364 double f = mass2 / q03 * (q2 * (h - h0) + (mass2 - sa) * q02 * dh);
366 a[0] = 1.0 + d * width_param / mass_param;
368 b[0] = mass2 - sa + width_param * f;
369 b[1] = -mass_param * width_param * widl1(mass2, mass_param, sa, sb, sc, r2);
370 Com_Divide(a, b, prop);
373 void EvtDsToKpipi::Flatte_rhoab(
double sa,
double sb,
double sc,
double rho_param[2])
375 double q = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
377 rho_param[0] = 2 *
sqrt(q / sa);
381 rho_param[1] = 2 *
sqrt(-q / sa);
386 void EvtDsToKpipi::propagatorKstr1430(
double mass_param,
double sx,
double* sb,
double* sc,
double prop[2])
388 double unit[2] = {1.0};
389 double ci[2] = {0, 1};
391 Flatte_rhoab(sx, sb[0], sc[0], rho1);
393 Flatte_rhoab(sx, sb[1], sc[1], rho2);
394 double gKPi_Kstr1430 = 0.2990, gEtaPK_Kstr1430 = 0.0529;
395 double tmp1[2] = {gKPi_Kstr1430, 0};
397 double tmp2[2] = {gEtaPK_Kstr1430, 0};
399 Com_Multi(tmp1, rho1, tmp11);
400 Com_Multi(tmp2, rho2, tmp22);
401 double tmp3[2] = {tmp11[0] + tmp22[0], tmp11[1] + tmp22[1]};
403 Com_Multi(tmp3, ci, tmp31);
404 double tmp4[2] = {mass_param* mass_param - sx - tmp31[0], -1.0 * tmp31[1]};
405 Com_Divide(unit, tmp4, prop);
408 double EvtDsToKpipi::DDalitz(
double P1[4],
double P2[4],
double P3[4],
int Ang,
double mass_param)
411 double temp_PDF, v_re;
414 double B[2], s1, s2, s3, sR, sD;
415 for (
int i = 0; i < 4; i++) {
416 pR[i] = P1[i] + P2[i];
417 pD[i] = pR[i] + P3[i];
425 for (
int i = 0; i != 4; i++) {
426 for (
int j = 0; j != 4; j++) {
428 if (i == 0) GG[i][j] = 1;
439 B[0] = barrier(1, sR, s1, s2, 3.0, mass_param);
440 B[1] = barrier(1, sD, sR, s3, 5.0, mDsM);
445 for (
int i = 0; i < 4; i++) {
446 temp_PDF += t1[i] * T1[i] * GG[i][i];
450 B[0] = barrier(2, sR, s1, s2, 3.0, mass_param);
451 B[1] = barrier(2, sD, sR, s3, 5.0, mDsM);
452 double t2[4][4], T2[4][4];
456 for (
int i = 0; i < 4; i++) {
457 for (
int j = 0; j < 4; j++) {
458 temp_PDF += t2[i][j] * T2[j][i] * GG[i][i] * GG[j][j];
462 v_re = temp_PDF * B[0] * B[1];
467 void EvtDsToKpipi::rhoab(
double sa,
double sb,
double sc,
double res[2])
469 double tmp = sa + sb - sc;
470 double q = 0.25 * tmp * tmp / sa - sb;
472 res[0] = 2.0 *
sqrt(q / sa);
476 res[1] = 2.0 *
sqrt(-q / sa);
479 void EvtDsToKpipi::rho4Pi(
double sa,
double res[2])
481 double temp = 1.0 - 0.3116765584 / sa;
483 res[0] =
sqrt(temp) / (1.0 + exp(9.8 - 3.5 * sa));
487 res[1] =
sqrt(-temp) / (1.0 + exp(9.8 - 3.5 * sa));
491 void EvtDsToKpipi::propagatorsigma500(
double sa,
double sb,
double sc,
double prop[2])
493 double f = 0.5843 + 1.6663 * sa;
494 const double M = 0.9264;
495 const double mass2 = 0.85821696;
496 const double mpi2d2 = 0.00973989245;
497 double g1 = f * (sa - mpi2d2) / (mass2 - mpi2d2) * exp((mass2 - sa) / 1.082);
498 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
499 rhoab(sa, sb, sc, rho1s);
500 rhoab(mass2, sb, sc, rho1M);
502 rho4Pi(mass2, rho2M);
503 Com_Divide(rho1s, rho1M, rho1);
504 Com_Divide(rho2s, rho2M, rho2);
508 b[0] = mass2 - sa + M * (g1 * rho1[1] + 0.0024 * rho2[1]);
509 b[1] = -M * (g1 * rho1[0] + 0.0024 * rho2[0]);
510 Com_Divide(a, b, prop);
514 void EvtDsToKpipi::calEva(
double*
K,
double* Pi1,
double* Pi2,
double* mass1,
double* width1,
double* amp,
double* phase,
int* g0,
515 int* spin_param,
int* modetype_param,
int nstates,
double& Result)
517 double P23[4], P13[4];
518 double cof[2], amp_PDF[2], PDF[2];
520 for (
int i = 0; i < 4; i++) {
521 P13[i] =
K[i] + Pi2[i];
522 P23[i] = Pi1[i] + Pi2[i];
524 s13 = SCADot(P13, P13);
525 s23 = SCADot(P23, P23);
528 s2 = SCADot(Pi1, Pi1);
529 s3 = SCADot(Pi2, Pi2);
530 double pro[2], temp_PDF, amp_tmp[2];
538 for (
int i = 0; i < nstates; i++) {
541 cof[0] = amp[i] * cos(phase[i]);
542 cof[1] = amp[i] * sin(phase[i]);
545 if (modetype_param[i] == 13) {
546 temp_PDF = DDalitz(
K, Pi2, Pi1, spin_param[i], mass1[i]);
547 if (g0[i] == 1) propagatorRBW(mass1[i], width1[i], s13, mKa2, mPi2, rRess, spin_param[i], pro);
549 double skm2[2] = {s1, 0.95778 * 0.95778};
550 double spi2[2] = {s3, mass_Kaon * mass_Kaon};
551 propagatorKstr1430(mass1[i], s13, skm2, spi2, pro);
557 amp_tmp[0] = temp_PDF * pro[0];
558 amp_tmp[1] = temp_PDF * pro[1];
561 if (modetype_param[i] == 23) {
562 temp_PDF = DDalitz(Pi1, Pi2,
K, spin_param[i], mass1[i]);
563 if (g0[i] == 4) propagatorFlatte(mass1[i], width1[i], s23, pro);
564 if (g0[i] == 3) propagatorCBW(mass1[i], width1[i], s23, pro);
565 if (g0[i] == 2) propagatorRBW(mass1[i], width1[i], s23, mPi2, mPi2, rRess, spin_param[i], pro);
566 if (g0[i] == 1) propagatorGS(mass1[i], width1[i], s23, mPi2, mPi2, 9.0, pro);
567 if (g0[i] == 500) propagatorsigma500(s23, s2, s3, pro);
572 amp_tmp[0] = temp_PDF * pro[0];
573 amp_tmp[1] = temp_PDF * pro[1];
575 Com_Multi(amp_tmp, cof, amp_PDF);
576 PDF[0] += amp_PDF[0];
577 PDF[1] += amp_PDF[1];
580 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
581 if (value <= 0) value = 1e-20;
#define B2_EVTGEN_REGISTER_MODEL(classname)
Class to register B2_EVTGEN_REGISTER_MODEL.
double sqrt(double a)
sqrt for double
Abstract base class for different kinds of events.