Belle II Software development
EvtDToKSpipipi.cc
1// Model: EvtDToKSpipipi
2// This file is an amplitude model for D+ -> K_S0 pi+ pi+ pi-.
3// The model is from the BESIII Collaboration in PRD 100, 072008 (2019). DOI:  https://doi.org/10.1103/PhysRevD.100.072008
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/EvtComplex.hh>
15#include <EvtGenBase/EvtDecayTable.hh>
16#include <stdlib.h>
17
18#include <generators/evtgen/EvtGenModelRegister.h>
19#include <generators/evtgen/models/besiii/EvtDToKSpipipi.h>
20
21namespace Belle2 {
26
28 B2_EVTGEN_REGISTER_MODEL(EvtDToKSpipipi);
29
30 EvtDToKSpipipi::~EvtDToKSpipipi() {}
31
32 std::string EvtDToKSpipipi::getName()
33 {
34 return "DToKSpipipi";
35 }
36
37 EvtDecayBase* EvtDToKSpipipi::clone()
38 {
39 return new EvtDToKSpipipi;
40 }
41
42 void EvtDToKSpipipi::init()
43 {
44 checkNArg(0);
45 checkNDaug(4);
46 checkSpinParent(EvtSpinType::SCALAR);
47 /*
48 checkSpinDaughter(0,EvtSpinType::SCALAR);
49 checkSpinDaughter(1,EvtSpinType::SCALAR);
50 checkSpinDaughter(3,EvtSpinType::SCALAR);
51 checkSpinDaughter(4,EvtSpinType::SCALAR);
52 */
53
54 mK1270 = 1.272; mK1400 = 1.403;
55 GK1270 = 0.09; GK1400 = 0.174;
56 mKstr = 0.89166; mrho = 0.77549;
57 GKstr = 0.0462; Grho = 0.1491;
58 msigma = 0.472;
59 Gsigma = 0.542;
60 phi_omega = -0.02;
61 mK1650 = 1.65;
62 GK1650 = 0.158;
63 rho[0] = 1.0;
64 phi[0] = 0.0;
65
66 ma1 = 1.22;
67 Ga1 = 0.4282;
68 mK1460 = 1.4152;
69 GK1460 = 0.2485;
70 rho_omega = 0.00294;
71
72 phi[1] = -1.55;
73 rho[1] = 0.5843;
74
75 phi[2] = -1.8223;
76 rho[2] = 2.0974;
77
78 phi[3] = -2.6751;
79 rho[3] = 0.46642;
80
81 phi[4] = -2.2429;
82 rho[4] = 0.33334;
83
84 phi[5] = -0.55888;
85 rho[5] = 0.15549;
86
87 phi[6] = -1.8778;
88 rho[6] = 0.94452;
89
90 phi[7] = 2.7724;
91 rho[7] = 0.99411;
92
93 phi[8] = -2.6461;
94 rho[8] = 0.21231;
95
96 phi[9] = -0.95137;
97 rho[9] = 0.29895;
98
99 phi[10] = -3.0843;
100 rho[10] = 3.6361;
101
102 phi[11] = 2.0954;
103 rho[11] = 0.96472;
104
105 phi[12] = -2.4965;
106 rho[12] = 0.48470;
107
108 mD = 1.86486;
109 rD = 5;
110 metap = 0.95778;
111 mkstr = 0.89594;
112 mk0 = 0.497614;
113 mass_Kaon = 0.49368;
114 mass_Pion = 0.13957;
115 math_pi = 3.1415926;
116
117 pi = 3.1415926;
118 mpi = 0.13957;
119 g1 = 0.5468;
120 g2 = 0.23;
121
122 int GG[4][4] = { {1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1} };
123 int EE[4][4][4][4] = {
124 { {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0} },
125 {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 1}, {0, 0, -1, 0}},
126 {{0, 0, 0, 0}, {0, 0, 0, -1}, {0, 0, 0, 0}, {0, 1, 0, 0} },
127 {{0, 0, 0, 0}, {0, 0, 1, 0}, {0, -1, 0, 0}, {0, 0, 0, 0} }
128 },
129 { {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, -1}, {0, 0, 1, 0} },
130 {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0} },
131 {{0, 0, 0, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}, {-1, 0, 0, 0}},
132 {{0, 0, -1, 0}, {0, 0, 0, 0}, {1, 0, 0, 0}, {0, 0, 0, 0} }
133 },
134 { {{0, 0, 0, 0}, {0, 0, 0, 1}, {0, 0, 0, 0}, {0, -1, 0, 0}},
135 {{0, 0, 0, -1}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 0, 0, 0} },
136 {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0} },
137 {{0, 1, 0, 0}, {-1, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0} }
138 },
139 { {{0, 0, 0, 0}, {0, 0, -1, 0}, {0, 1, 0, 0}, {0, 0, 0, 0} },
140 {{0, 0, 1, 0}, {0, 0, 0, 0}, {-1, 0, 0, 0}, {0, 0, 0, 0} },
141 {{0, -1, 0, 0}, {1, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0} },
142 {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0} }
143 }
144 };
145 for (int i = 0; i < 4; i++) {
146 for (int j = 0; j < 4; j++) {
147 G[i][j] = GG[i][j];
148 for (int k = 0; k < 4; k++) {
149 for (int l = 0; l < 4; l++) {
150 E[i][j][k][l] = EE[i][j][k][l];
151 }
152 }
153 }
154 }
155 }
156
157 void EvtDToKSpipipi::initProbMax()
158 {
159 setProbMax(3700.0);
160 }
161
162 void EvtDToKSpipipi::decay(EvtParticle* p)
163 {
164 /*
165 double maxprob = 0.0;
166 for(int ir=0;ir<=60000000;ir++){
167 p->initializePhaseSpace(getNDaug(),getDaugs());
168 EvtVector4R Ks0 = p->getDaug(0)->getP4();
169 EvtVector4R pi1 = p->getDaug(1)->getP4();
170 EvtVector4R pi2 = p->getDaug(2)->getP4();
171 EvtVector4R pi3 = p->getDaug(3)->getP4();
172 double Ks[4],Pip1[4],Pip2[4],Pim[4];
173 Ks[0] = Ks0.get(0); Pip1[0] = pi1.get(0); Pip2[0] = pi2.get(0); Pim[0] = pi3.get(0);
174 Ks[1] = Ks0.get(1); Pip1[1] = pi1.get(1); Pip2[1] = pi2.get(1); Pim[1] = pi3.get(1);
175 Ks[2] = Ks0.get(2); Pip1[2] = pi1.get(2); Pip2[2] = pi2.get(2); Pim[2] = pi3.get(2);
176 Ks[3] = Ks0.get(3); Pip1[3] = pi1.get(3); Pip2[3] = pi2.get(3); Pim[3] = pi3.get(3);
177 double Prob = calPDF(Ks, Pip1, Pip2, Pim);
178 if(Prob>maxprob) {
179 maxprob=Prob;
180 }
181 }
182 */
183 p->initializePhaseSpace(getNDaug(), getDaugs());
184 EvtVector4R Ks0 = p->getDaug(0)->getP4();
185 EvtVector4R pi1 = p->getDaug(1)->getP4();
186 EvtVector4R pi2 = p->getDaug(2)->getP4();
187 EvtVector4R pi3 = p->getDaug(3)->getP4();
188
189 double Ks[4], Pip1[4], Pip2[4], Pim[4];
190 Ks[0] = Ks0.get(0); Pip1[0] = pi1.get(0); Pip2[0] = pi2.get(0); Pim[0] = pi3.get(0);
191 Ks[1] = Ks0.get(1); Pip1[1] = pi1.get(1); Pip2[1] = pi2.get(1); Pim[1] = pi3.get(1);
192 Ks[2] = Ks0.get(2); Pip1[2] = pi1.get(2); Pip2[2] = pi2.get(2); Pim[2] = pi3.get(2);
193 Ks[3] = Ks0.get(3); Pip1[3] = pi1.get(3); Pip2[3] = pi2.get(3); Pim[3] = pi3.get(3);
194 double prob = calPDF(Ks, Pip1, Pip2, Pim);
195 setProb(prob);
196 return;
197 }
198
199 double EvtDToKSpipipi::calPDF(double Ks[], double Pip1[], double Pip2[], double Pim[])
200 {
201 EvtComplex PDF[100];
202 double P14[4], P24[4], P34[4];
203 for (int i = 0; i < 4; i++) {
204 P14[i] = Ks[i] + Pim[i];
205 P24[i] = Pip1[i] + Pim[i];
206 P34[i] = Pip2[i] + Pim[i];
207 }
208 //----------D->a1Ks--------------
209 //----------a1->rhoPi------------
210 PDF[0] = D2AP_A2VP(Ks, Pip2, Pip1, Pim, 0) * getprop(P24, Pip2, ma1, Ga1, 1, 0) *
211 getprop(Pip1, Pim, mrho, Grho, 2, 1) +
212 D2AP_A2VP(Ks, Pip1, Pip2, Pim, 0) * getprop(P34, Pip1, ma1, Ga1, 1, 0) *
213 getprop(Pip2, Pim, mrho, Grho, 2, 1);
214 PDF[1] = D2AP_A2VP(Ks, Pip2, Pip1, Pim, 2) * getprop(P24, Pip2, ma1, Ga1, 1, 2) *
215 getprop(Pip1, Pim, mrho, Grho, 2, 1) +
216 D2AP_A2VP(Ks, Pip1, Pip2, Pim, 2) * getprop(P34, Pip1, ma1, Ga1, 1, 2) *
217 getprop(Pip2, Pim, mrho, Grho, 2, 1);
218 //----------a1->sigma pi---------
219 PDF[2] = D2AP_A2SP(Ks, Pip2, Pip1, Pim) * getprop(P24, Pip2, ma1, Ga1, 1, 1) *
220 getprop(Pip1, Pim, msigma, Gsigma, 4, 0) +
221 D2AP_A2SP(Ks, Pip1, Pip2, Pim) * getprop(P34, Pip1, ma1, Ga1, 1, 1) *
222 getprop(Pip2, Pim, msigma, Gsigma, 4, 0);
223 //----------D->a1K finish-----
224 //---------D->K1(1400) pi-----
225 //K1400[S]->K* pi
226 PDF[3] = D2AP_A2VP(Pip2, Pip1, Ks, Pim, 0) * getprop(P14, Pip1, mK1400, GK1400, 1, 0) *
227 getprop(Ks, Pim, mKstr, GKstr, 1, 1) +
228 D2AP_A2VP(Pip1, Pip2, Ks, Pim, 0) * getprop(P14, Pip2, mK1400, GK1400, 1, 0) *
229 getprop(Ks, Pim, mKstr, GKstr, 1, 1);
230 //K1400[D]->K* pi
231 PDF[4] = D2AP_A2VP(Pip2, Pip1, Ks, Pim, 2) * getprop(P14, Pip1, mK1400, GK1400, 1, 2) *
232 getprop(Ks, Pim, mKstr, GKstr, 1, 1) +
233 D2AP_A2VP(Pip1, Pip2, Ks, Pim, 2) * getprop(P14, Pip2, mK1400, GK1400, 1, 2) *
234 getprop(Ks, Pim, mKstr, GKstr, 1, 1);
235 //-----------------------------
236 //-------K1270[S]->Ksrho-------
237 PDF[5] = D2AP_A2VP(Pip2, Ks, Pip1, Pim, 0) * getprop(P24, Ks, mK1270, GK1270, 0, 0) *
238 getprop(Pip1, Pim, mrho, Grho, 2, 1) +
239 D2AP_A2VP(Pip1, Ks, Pip2, Pim, 0) * getprop(P34, Ks, mK1270, GK1270, 0, 0) *
240 getprop(Pip2, Pim, mrho, Grho, 2, 1);
241 //-------D->rhoKAD------------
242 PDF[6] = D2AP_A2VP(Pip2, Ks, Pip1, Pim, 0) * getprop(Pip1, Pim, mrho, Grho, 2, 1) +
243 D2AP_A2VP(Pip1, Ks, Pip2, Pim, 0) * getprop(Pip2, Pim, mrho, Grho, 2, 1);
244 PDF[7] = D2AP_A2VP(Pip2, Ks, Pip1, Pim, 2) * getprop(Pip1, Pim, mrho, Grho, 2, 1) +
245 D2AP_A2VP(Pip1, Ks, Pip2, Pim, 2) * getprop(Pip2, Pim, mrho, Grho, 2, 1);
246 //-------D->K1460, K1460->Ks rho---------
247 PDF[8] = D2PP_P2VP(Pip2, Ks, Pip1, Pim) * getprop(P24, Ks, mK1460, GK1460, 1, 1) *
248 getprop(Pip1, Pim, mrho, Grho, 2, 1) +
249 D2PP_P2VP(Pip1, Ks, Pip2, Pim) * getprop(P34, Ks, mK1460, GK1460, 1, 1) *
250 getprop(Pip2, Pim, mrho, Grho, 2, 1);
251 //--------K*PiA (K1650)---------------------
252 PDF[9] = D2AP_A2VP(Pip2, Pip1, Ks, Pim, 0) * getprop(P14, Pip1, mK1650, GK1650, 1, 0) *
253 getprop(Ks, Pim, mKstr, GKstr, 1, 1) +
254 D2AP_A2VP(Pip1, Pip2, Ks, Pim, 0) * getprop(P14, Pip2, mK1650, GK1650, 1, 0) *
255 getprop(Ks, Pim, mKstr, GKstr, 1, 1);
256 //-------KsPiPiSPiA-----------------
257 PDF[10] = D2AP_A2SP(Pip2, Ks, Pip1, Pim) + D2AP_A2SP(Pip1, Ks, Pip2, Pim);
258 //-------KPiS wave------------------
259 EvtComplex corr(2, 0);
260 PDF[11] = corr * PHSP(Ks, Pim);
261 //-------D->K1460pi, K1460->K*-pi+--------
262 PDF[12] = D2PP_P2VP(Pip2, Pip1, Ks, Pim) * getprop(P14, Pip1, mK1460, GK1460, 1, 1) *
263 getprop(Ks, Pim, mKstr, GKstr, 1, 1) +
264 D2PP_P2VP(Pip1, Pip2, Ks, Pim) * getprop(P14, Pip2, mK1460, GK1460, 1, 1) *
265 getprop(Ks, Pim, mKstr, GKstr, 1, 1);
266 //-------------------------------------------
267 EvtComplex cof;
268 EvtComplex pdf, module;
269 pdf = EvtComplex(0, 0);
270 for (int i = 0; i < 13; i++) {
271 cof = EvtComplex(rho[i] * cos(phi[i]), rho[i] * sin(phi[i]));
272 pdf = pdf + cof * PDF[i];
273 }
274 module = conj(pdf) * pdf;
275 double value;
276 value = real(module);
277 return (value <= 0) ? 1e-20 : value;
278 }
279
280 EvtComplex EvtDToKSpipipi::KPiSFormfactor(const double sa, const double sb, const double sc, const double r)
281 {
282 (void)r;
283 double m1430 = 1.463;
284 double sa0 = m1430 * m1430;
285 double w1430 = 0.233;
286 double q0 = (sa0 + sb - sc) * (sa0 + sb - sc) / (4 * sa0) - sb;
287 if (q0 < 0) q0 = 1e-16;
288 double qs = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
289 double q = sqrt(qs);
290 double width_ = w1430 * q * m1430 / sqrt(sa * q0);
291 double temp_R = atan(m1430 * width_ / (sa0 - sa));
292 if (temp_R < 0) temp_R += math_pi;
293 double deltaR = -5.31 + temp_R;
294 double temp_F = atan(2 * 1.07 * q / (2 + (-1.8) * 1.07 * qs));
295 if (temp_F < 0) temp_F += math_pi;
296 double deltaF = 2.33 + temp_F;
297 EvtComplex cR(cos(deltaR), sin(deltaR));
298 EvtComplex cF(cos(deltaF), sin(deltaF));
299 EvtComplex amp = 0.8 * sin(deltaF) * cF + sin(deltaR) * cR * cF * cF;
300 return amp;
301 }
302 EvtComplex EvtDToKSpipipi::D2AP_A2VP(double P1[], double P2[], double P3[], double P4[], int L)
303 {
304 double temp_PDF = 0;
305 EvtComplex amp_PDF(0, 0);
306 double t1V[4], t1D[4], t2A[4][4];
307 double sa[3], sb[3], sc[3], B[3];
308 double pV[4], pA[4], pD[4];
309 for (int i = 0; i != 4; i++) {
310 pV[i] = P3[i] + P4[i];
311 pA[i] = pV[i] + P2[i];
312 pD[i] = pA[i] + P1[i];
313 }
314 sa[0] = dot(pV, pV);
315 sb[0] = dot(P3, P3);
316 sc[0] = dot(P4, P4);
317 sa[1] = dot(pA, pA);
318 sb[1] = sa[0];
319 sc[1] = dot(P2, P2);
320 sa[2] = dot(pD, pD);
321 sb[2] = sa[1];
322 sc[2] = dot(P1, P1);
323 B[0] = barrier(1, sa[0], sb[0], sc[0], 3);
324 B[2] = barrier(1, sa[2], sb[2], sc[2], rD);
325 calt1(P3, P4, t1V);
326 calt1(pA, P1, t1D);
327 if (L == 0) {
328 for (int i = 0; i != 4; i++) {
329 for (int j = 0; j != 4; j++) {
330 temp_PDF += t1D[i] * (G[i][j] - pA[i] * pA[j] / sa[1]) * t1V[j] * (G[i][i]) * (G[j][j]);
331 }
332 }
333 B[1] = 1;
334 }
335 if (L == 2) {
336 calt2(pV, P2, t2A);
337 for (int i = 0; i != 4; i++) {
338 for (int j = 0; j != 4; j++) {
339 temp_PDF += t1D[i] * t2A[i][j] * t1V[j] * (G[i][i]) * (G[j][j]);
340 }
341 }
342 B[1] = barrier(2, sa[1], sb[1], sc[1], 3);
343 }
344 amp_PDF = temp_PDF * B[0] * B[1] * B[2];
345 return amp_PDF;
346 }
347 EvtComplex EvtDToKSpipipi::D2AP_A2SP(double P1[], double P2[], double P3[], double P4[])
348 {
349 double temp_PDF = 0;
350 EvtComplex amp_PDF(0, 0);
351 EvtComplex pro;
352 double sa[3], sb[3], sc[3], B[3];
353 double t1D[4], t1A[4];
354 double pS[4], pA[4], pD[4];
355 for (int i = 0; i != 4; i++) {
356 pS[i] = P3[i] + P4[i];
357 pA[i] = pS[i] + P2[i];
358 pD[i] = pA[i] + P1[i];
359 }
360 sa[0] = dot(pS, pS);
361 sb[0] = dot(P3, P3);
362 sc[0] = dot(P4, P4);
363 sa[1] = dot(pA, pA);
364 sb[1] = sa[0];
365 sc[1] = dot(P2, P2);
366 sa[2] = dot(pD, pD);
367 sb[2] = sa[1];
368 sc[2] = dot(P1, P1);
369 B[1] = barrier(1, sa[1], sb[1], sc[1], 3);
370 B[2] = barrier(1, sa[2], sb[2], sc[2], rD);
371 calt1(pA, P1, t1D);
372 calt1(pS, P2, t1A);
373 for (int i = 0; i != 4; i++) {
374 temp_PDF += t1D[i] * t1A[i] * (G[i][i]);
375 }
376 amp_PDF = temp_PDF * B[1] * B[2];
377 return amp_PDF;
378 }
379 EvtComplex EvtDToKSpipipi::D2PP_P2VP(double P1[], double P2[], double P3[], double P4[])
380 {
381 double temp_PDF = 0;
382 EvtComplex amp_PDF(0, 0);
383 EvtComplex pro;
384 double sa[3], sb[3], sc[3], B[3];
385 double t1V[4];
386 double pV[4], pP[4], pD[4];
387 for (int i = 0; i != 4; i++) {
388 pV[i] = P3[i] + P4[i];
389 pP[i] = pV[i] + P2[i];
390 pD[i] = pP[i] + P1[i];
391 }
392 sa[0] = dot(pV, pV);
393 sb[0] = dot(P3, P3);
394 sc[0] = dot(P4, P4);
395 sa[1] = dot(pP, pP);
396 sb[1] = sa[0];
397 sc[1] = dot(P2, P2);
398 sa[2] = dot(pD, pD);
399 sb[2] = sa[1];
400 sc[2] = dot(P1, P1);
401 B[0] = barrier(1, sa[0], sb[0], sc[0], 3);
402 B[1] = barrier(1, sa[1], sb[1], sc[1], 3);
403 calt1(P3, P4, t1V);
404 for (int i = 0; i != 4; i++) {
405 temp_PDF += P2[i] * t1V[i] * (G[i][i]);
406 }
407 amp_PDF = temp_PDF * B[0] * B[1];
408 return amp_PDF;
409 }
410 EvtComplex EvtDToKSpipipi::PHSP(double Km[], double Pip[])
411 {
412 EvtComplex amp_PDF(0, 0);
413 double sa, sb, sc;
414 double KPi[4];
415 for (int i = 0; i != 4; i++) {
416 KPi[i] = Km[i] + Pip[i];
417 }
418 sa = dot(KPi, KPi);
419 sb = dot(Km, Km);
420 sc = dot(Pip, Pip);
421 amp_PDF = KPiSFormfactor(sa, sb, sc, 3.0);
422 return amp_PDF;
423 }
424 EvtComplex EvtDToKSpipipi::getprop(double daug1[], double daug2[], double mass_, double width_, int flag, int L)
425 {
426 //flag = 0, RBW with constant width
427 //flag = 1, RBW with width depends on momentum
428 //flag = 2, GS convolute with RBW_omega
429 //flag = 3, flatte for f_0(980)/a_0(980)
430 //flag = 4, Bugg formula for sigma;
431 //effect radii == 3.0 GeV^-1
432 EvtComplex prop1, prop2, prop;
433 EvtComplex one(1, 0);
434 double pR[4];
435 for (int i = 0; i < 4; i++) {
436 pR[i] = daug1[i] + daug2[i];
437 }
438 double sa, sb, sc;
439 sa = dot(pR, pR);
440 sb = dot(daug1, daug1);
441 sc = dot(daug2, daug2);
442 if (flag == 0) return propogator(mass_, width_, sa);
443 if (flag == 1) return propagatorRBW(mass_, width_, sa, sb, sc, 3.0, L);
444 if (flag == 2) {
445 prop1 = propagatorGS(mass_, width_, sa, sb, sc, 3.0, L);
446 prop2 = propagatorRBW(0.783, 0.008, sa, sb, sc, 3.0, L);
447 EvtComplex coef_omega(rho_omega * cos(phi_omega), rho_omega * sin(phi_omega));
448 prop = prop1 * (one + 0.783 * 0.783 * coef_omega * prop2);
449 return prop;
450 }
451 if (flag == 3) {
452 //Not need for D+ -> Ks 3pi
453 }
454 if (flag == 4) {
455 EvtComplex ci(0, 1);
456 double f = 0.5843 + 1.6663 * sa;
457 double M = 0.9264;
458 double mpi2 = mass_Pion * mass_Pion;
459 double mass2 = M * M;
460 double g1_ = f * (sa - mpi2 / 2) / (mass2 - mpi2 / 2) * exp((mass2 - sa) / 1.082);
461 EvtComplex rho1s = rhoab(sa, sb, sc);
462 EvtComplex rho1M = rhoab(mass2, sb, sc);
463 EvtComplex rho2s = rho4Pi(sa);
464 EvtComplex rho2M = rho4Pi(mass2);
465 prop = 1.0 / (M * M - sa - ci * M * (g1_ * rho1s / rho1M + 0.0024 * rho2s / rho2M));
466 return prop;
467 }
468 return one;
469 }
470 EvtComplex EvtDToKSpipipi::rhoab(const double sa, const double sb, const double sc)
471 {
472 EvtComplex one(1, 0);
473 double q = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
474 EvtComplex rho_;
475 EvtComplex ci(0, 1);
476 if (q > 0) rho_ = one * sqrt(q / sa);
477 if (q < 0) rho_ = ci * sqrt(-q / sa);
478 rho_ = 2.0 * rho_;
479 return rho_;
480 }
481 EvtComplex EvtDToKSpipipi::rho4Pi(const double sa)
482 {
483 double mpi_ = 0.13957;
484 EvtComplex one(1, 0);
485 EvtComplex rho_(0, 0);
486 EvtComplex ci(0, 1);
487 double temp = 1 - 16 * mpi_ * mpi_ / sa;
488 if (temp > 0) rho_ = one * sqrt(temp) / (1 + exp(9.8 - 3.5 * sa));
489 if (temp < 0) rho_ = ci * sqrt(-temp) / (1 + exp(9.8 - 3.5 * sa));
490 return rho_;
491 }
492
493 EvtComplex EvtDToKSpipipi::propogator(const double mass_, const double width_, const double sx)const
494 {
495 EvtComplex ci(0, 1);
496 EvtComplex prop = 1.0 / (mass_ * mass_ - sx - ci * mass_ * width_);
497 return prop;
498 }
499 double EvtDToKSpipipi::wid(double mass_, double sa, double sb, double sc, double r, int l)const
500 {
501 double widm(0.), q(0.), q0(0.);
502 double sa0 = mass_ * mass_;
503 double m = sqrt(sa);
504 q = Qabcs(sa, sb, sc);
505 q0 = Qabcs(sa0, sb, sc);
506 double z, z0;
507 z = q * r * r;
508 z0 = q0 * r * r;
509 double t = q / q0;
510 double F(0.);
511 if (l == 0) F = 1;
512 if (l == 1) F = sqrt((1 + z0) / (1 + z));
513 if (l == 2) F = sqrt((9 + 3 * z0 + z0 * z0) / (9 + 3 * z + z * z));
514 widm = pow(t, l + 0.5) * mass_ / m * F * F;
515 return widm;
516 }
517 double EvtDToKSpipipi::h(const double m, const double q)const
518 {
519 double h(0.);
520 h = 2 / pi * q / m * log((m + 2 * q) / (2 * mpi));
521 return h;
522 }
523 double EvtDToKSpipipi::dh(const double mass_, const double q0)const
524 {
525 double dh = h(mass_, q0) * (1.0 / (8 * q0 * q0) - 1.0 / (2 * mass_ * mass_)) + 1.0 / (2 * pi * mass_ * mass_);
526 return dh;
527 }
528 double EvtDToKSpipipi::f(const double mass_, const double sx, const double q0, const double q)const
529 {
530 double m = sqrt(sx);
531 double f = mass_ * mass_ / (pow(q0, 3)) * (q * q * (h(m, q) - h(mass_, q0)) + (mass_ * mass_ - sx) * q0 * q0 * dh(mass_, q0));
532 return f;
533 }
534 double EvtDToKSpipipi::d(const double mass_, const double q0)const
535 {
536 double d = 3.0 / pi * mpi * mpi / (q0 * q0) * log((mass_ + 2 * q0) / (2 * mpi)) + mass_ / (2 * pi * q0) - (mpi * mpi * mass_) /
537 (pi * pow(q0, 3));
538 return d;
539 }
540 EvtComplex EvtDToKSpipipi::propagatorRBW(const double mass_, const double width_, const double sa, const double sb,
541 const double sc, const double r, const int l)const
542 {
543 EvtComplex ci(0, 1);
544 EvtComplex prop = 1.0 / (mass_ * mass_ - sa - ci * mass_ * width_ * wid(mass_, sa, sb, sc, r, l));
545 return prop;
546 }
547 EvtComplex EvtDToKSpipipi::propagatorGS(const double mass_, const double width_, const double sa, const double sb, const double sc,
548 const double r, const int l)const
549 {
550 EvtComplex ci(0, 1);
551 double q = Qabcs(sa, sb, sc);
552 double sa0 = mass_ * mass_;
553 double q0 = Qabcs(sa0, sb, sc);
554 q = sqrt(q);
555 q0 = sqrt(q0);
556 EvtComplex prop = (1 + d(mass_, q0) * width_ / mass_) / (mass_ * mass_ - sa + width_ * f(mass_, sa, q0,
557 q) - ci * mass_ * width_ * wid(mass_, sa, sb, sc, r, l));
558 return prop;
559 }
560 double EvtDToKSpipipi::dot(double* a1, double* a2)const
561 {
562 double dot = 0;
563 for (int i = 0; i != 4; i++) {
564 dot += a1[i] * a2[i] * G[i][i];
565 }
566 return dot;
567 }
568 double EvtDToKSpipipi::Qabcs(const double sa, const double sb, const double sc)const
569 {
570 double Qabcs = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
571 if (Qabcs < 0) Qabcs = 1e-16;
572 return Qabcs;
573 }
574 double EvtDToKSpipipi::barrier(const double l, const double sa, const double sb, const double sc, const double r)const
575 {
576 double q = Qabcs(sa, sb, sc);
577 q = sqrt(q);
578 double z = q * r;
579 z = z * z;
580 double F = 1;
581 if (l > 2) F = 0;
582 if (l == 0)
583 F = 1;
584 if (l == 1) {
585 F = sqrt((2 * z) / (1 + z));
586 }
587 if (l == 2) {
588 F = sqrt((13 * z * z) / (9 + 3 * z + z * z));
589 }
590 return F;
591 }
592 void EvtDToKSpipipi::calt1(double daug1[], double daug2[], double t1[]) const
593 {
594 double p, pq;
595 double pa[4], qa[4];
596 for (int i = 0; i != 4; i++) {
597 pa[i] = daug1[i] + daug2[i];
598 qa[i] = daug1[i] - daug2[i];
599 }
600 p = dot(pa, pa);
601 pq = dot(pa, qa);
602 for (int i = 0; i != 4; i++) {
603 t1[i] = qa[i] - pq / p * pa[i];
604 }
605 }
606 void EvtDToKSpipipi::calt2(double daug1[], double daug2[], double t2[][4]) const
607 {
608 double p, r;
609 double pa[4], t1[4];
610 calt1(daug1, daug2, t1);
611 r = dot(t1, t1);
612 for (int i = 0; i != 4; i++) {
613 pa[i] = daug1[i] + daug2[i];
614 }
615 p = dot(pa, pa);
616 for (int i = 0; i != 4; i++) {
617 for (int j = 0; j != 4; j++) {
618 t2[i][j] = t1[i] * t1[j] - 1.0 / 3 * r * (G[i][j] - pa[i] * pa[j] / p);
619 }
620 }
621 }
622
624} // Belle 2 namespace
Belle2::EvtDToKSpipipi::mD
double mD
Fixed parameters.
Definition EvtDToKSpipipi.h:91
Belle2::EvtDToKSpipipi::getprop
EvtComplex getprop(double daug1[], double daug2[], double mass, double width, int flag, int L)
Propagator Lineshapes.
Definition EvtDToKSpipipi.cc:424
Belle2::EvtDToKSpipipi::f
double f(const double mass, const double sx, const double q0, const double q) const
f function in Gounaris-Sakurai lineshape
Definition EvtDToKSpipipi.cc:528
Belle2::EvtDToKSpipipi::D2AP_A2VP
EvtComplex D2AP_A2VP(double P1[], double P2[], double P3[], double P4[], int L)
Amplitude modes.
Definition EvtDToKSpipipi.cc:302
Belle2::EvtDToKSpipipi::propagatorGS
EvtComplex propagatorGS(const double mass, const double width, const double sa, const double sb, const double sc, const double r, const int l) const
Gounaris-Sakurai lineshape Function.
Definition EvtDToKSpipipi.cc:547
Belle2::EvtDToKSpipipi::calPDF
double calPDF(double Km[], double Pip1[], double Pip2[], double Pim[])
Probability distribution function of the decay.
Definition EvtDToKSpipipi.cc:199
Belle2::EvtDToKSpipipi::Qabcs
double Qabcs(const double sa, const double sb, const double sc) const
Magnitudes of daughter particle momenta in the rest system of the mother particle.
Definition EvtDToKSpipipi.cc:568
Belle2::EvtDToKSpipipi::calt1
void calt1(double daug1[], double daug2[], double t1[]) const
Covariant Spin-1 Projector.
Definition EvtDToKSpipipi.cc:592
Belle2::EvtDToKSpipipi::propogator
EvtComplex propogator(const double mass, const double width, const double sx) const
Relativistic Breit-Wigner Lineshape Function (Fixed Width)
Definition EvtDToKSpipipi.cc:493
Belle2::EvtDToKSpipipi::KPiSFormfactor
EvtComplex KPiSFormfactor(const double sa, const double sb, const double sc, const double r)
K pi S-wave form factor.
Definition EvtDToKSpipipi.cc:280
Belle2::EvtDToKSpipipi::propagatorRBW
EvtComplex propagatorRBW(const double mass, const double width, const double sa, const double sb, const double sc, const double r, const int l) const
Relativistic Breit-Wigner Lineshape Function.
Definition EvtDToKSpipipi.cc:540
Belle2::EvtDToKSpipipi::dh
double dh(const double mass, const double q0) const
derivative h function in Gounaris-Sakurai lineshape
Definition EvtDToKSpipipi.cc:523
Belle2::EvtDToKSpipipi::wid
double wid(const double mass, const double sa, const double sb, const double sc, const double r, const int l) const
Energy Dependent Width.
Definition EvtDToKSpipipi.cc:499
Belle2::EvtDToKSpipipi::barrier
double barrier(const double l, const double sa, const double sb, const double sc, const double r) const
Blatt-Weisskopf barrier factors.
Definition EvtDToKSpipipi.cc:574
Belle2::EvtDToKSpipipi::d
double d(const double mass, const double q0) const
d function in Gounaris-Sakurai lineshape
Definition EvtDToKSpipipi.cc:534
Belle2::atan
double atan(double a)
atan for double
Definition beamHelpers.h:34
Belle2::EvtDToKSpipipi::rho4Pi
EvtComplex rho4Pi(const double sa)
Two-body Phase-space Function (Two Pions)
Definition EvtDToKSpipipi.cc:481
Belle2::EvtDToKSpipipi::rhoab
EvtComplex rhoab(const double sa, const double sb, const double sc)
Two-body Phase-space Function.
Definition EvtDToKSpipipi.cc:470
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::EvtDToKSpipipi::calt2
void calt2(double daug1[], double daug2[], double t2[][4]) const
Covariant Spin-2 Projector.
Definition EvtDToKSpipipi.cc:606
Belle2::EvtDToKSpipipi::h
double h(const double m, const double q) const
h function in Gounaris-Sakurai lineshape
Definition EvtDToKSpipipi.cc:517
Belle2::EvtDToKSpipipi::dot
double dot(double *a1, double *a2) const
Four-Vector Scalar Product.
Definition EvtDToKSpipipi.cc:560
Belle2
Abstract base class for different kinds of events.
Definition MillepedeAlgorithm.h:17