Belle II Software development
EvtD0ToKpipipi.cc
1// Model: EvtD0ToKpipipi
2// This file is an amplitude model for D0 -> K- pi+ pi+ pi-.
3// The model is from the BESIII Collaboration in PRD 95, 072010 (2017). DOI:  https://doi.org/10.1103/PhysRevD.95.072010
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/EvtD0ToKpipipi.h>
20
21namespace Belle2 {
26
27 B2_EVTGEN_REGISTER_MODEL(EvtD0ToKpipipi);
28
29 EvtD0ToKpipipi::~EvtD0ToKpipipi() {}
30
31 std::string EvtD0ToKpipipi::getName()
32 {
33 return "D0ToKpipipi";
34 }
35
36 EvtDecayBase* EvtD0ToKpipipi::clone()
37 {
38 return new EvtD0ToKpipipi;
39 }
40
41 void EvtD0ToKpipipi::init()
42 {
43 checkNArg(0);
44 checkNDaug(4);
45 checkSpinParent(EvtSpinType::SCALAR);
46 /*
47 checkSpinDaughter(0,EvtSpinType::SCALAR);
48 checkSpinDaughter(1,EvtSpinType::SCALAR);
49 checkSpinDaughter(3,EvtSpinType::SCALAR);
50 checkSpinDaughter(4,EvtSpinType::SCALAR);
51 */
52
53 width[0] = 0.09;
54 width[1] = 0.044183653178315;
55 width[2] = 0.541879469380012;
56 width[3] = 0.148423336450619;
57 mass[0] = 1.272;
58 mass[1] = 0.894781734682169;
59 mass[2] = 1.3622013558915;
60 mass[3] = 0.779143408171384;
61
62 phi[0] = 2.34794687054858;
63 rho[0] = 0.0759345115620669;
64 phi[1] = -2.24641399153466;
65 rho[1] = 0.0383327604903577;
66 phi[2] = 2.48955684856045;
67 rho[2] = 0.0931445480476023;
68
69 phi[3] = 0;
70 rho[3] = 1;
71
72 phi[4] = -2.10558220063012;
73 rho[4] = 0.347041869435286;
74 phi[5] = 1.47445088061872;
75 phi[6] = 3.00243265559304;
76 rho[5] = 0.00965088341753795;
77 rho[6] = 0.120536507325731;
78 phi[7] = -2.45477499325158;
79 rho[7] = 0.101419048440676;
80 phi[8] = -1.35809992343491;
81 rho[8] = 4.28149643321317;
82 phi[9] = -2.45149221243198;
83 rho[9] = 0.339492272598394;
84 phi[10] = -0.17419389225461;
85 rho[10] = -0.143619437541254;
86
87 phi[11] = -2.08744386934208;
88 rho[11] = 0.296286583716349;
89 phi[12] = 0.;
90 rho[12] = 0.;
91 phi[13] = -0.432190571560873;
92 rho[13] = 0.657344690733276;
93 phi[14] = -1.39790294886865;
94 rho[14] = 1.71208007006123;
95 phi[15] = 1.58945300476228;
96 rho[15] = 3.58248347683687;
97 phi[16] = 2.58249107256307;
98 rho[16] = -1.10728829503506;
99 phi[17] = -0.163623135170955;
100 rho[17] = 1.70863070178363;
101 phi[18] = -0.134699023080211;
102 rho[18] = 0.567531283682344;
103 phi[19] = -2.12670610368279;
104 rho[19] = 0.276571752504914;
105 phi[20] = -1.3352622107357;
106 rho[20] = 0.416634203151278;
107
108 phi[21] = -2.91571684221842;
109 rho[21] = 0.423062298489176;
110 phi[22] = 2.4544220004327;
111 rho[22] = 1.4017194038459;
112 phi[23] = -2.23388390670423;
113 rho[23] = 4.11110400629068;
114
115 mD = 1.86486;
116 rRes = 3.0;
117 rD = 5.0;
118 metap = 0.95778;
119 mkstr = 0.89594;
120 mk0 = 0.497614;
121 mass_Kaon = 0.49368;
122 mass_Pion = 0.13957;
123 mass_Pi0 = 0.1349766;
124 math_pi = 3.1415926;
125
126 pi = 3.1415926;
127 mpi = 0.13957;
128 g1 = 0.5468;
129 g2 = 0.23;
130
131 int GG[4][4] = { {1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1} };
132 int EE[4][4][4][4] = {
133 { {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0} },
134 {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 1}, {0, 0, -1, 0}},
135 {{0, 0, 0, 0}, {0, 0, 0, -1}, {0, 0, 0, 0}, {0, 1, 0, 0} },
136 {{0, 0, 0, 0}, {0, 0, 1, 0}, {0, -1, 0, 0}, {0, 0, 0, 0} }
137 },
138 { {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, -1}, {0, 0, 1, 0} },
139 {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0} },
140 {{0, 0, 0, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}, {-1, 0, 0, 0}},
141 {{0, 0, -1, 0}, {0, 0, 0, 0}, {1, 0, 0, 0}, {0, 0, 0, 0} }
142 },
143 { {{0, 0, 0, 0}, {0, 0, 0, 1}, {0, 0, 0, 0}, {0, -1, 0, 0}},
144 {{0, 0, 0, -1}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 0, 0, 0} },
145 {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0} },
146 {{0, 1, 0, 0}, {-1, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0} }
147 },
148 { {{0, 0, 0, 0}, {0, 0, -1, 0}, {0, 1, 0, 0}, {0, 0, 0, 0} },
149 {{0, 0, 1, 0}, {0, 0, 0, 0}, {-1, 0, 0, 0}, {0, 0, 0, 0} },
150 {{0, -1, 0, 0}, {1, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0} },
151 {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0} }
152 }
153 };
154 for (int i = 0; i < 4; i++) {
155 for (int j = 0; j < 4; j++) {
156 G[i][j] = GG[i][j];
157 for (int k = 0; k < 4; k++) {
158 for (int l = 0; l < 4; l++) {
159 E[i][j][k][l] = EE[i][j][k][l];
160 }
161 }
162 }
163 }
164 }
165
166 void EvtD0ToKpipipi::initProbMax()
167 {
168 setProbMax(720.0);
169 }
170
171 void EvtD0ToKpipipi::decay(EvtParticle* p)
172 {
173 /*
174 double maxprob = 0.0;
175 for(int ir=0;ir<=60000000;ir++){
176 p->initializePhaseSpace(getNDaug(),getDaugs());
177 EvtVector4R Km0 = p->getDaug(0)->getP4();
178 EvtVector4R pi1 = p->getDaug(1)->getP4();
179 EvtVector4R pi2 = p->getDaug(2)->getP4();
180 EvtVector4R pi3 = p->getDaug(3)->getP4();
181 double Km[4],Pip1[4],Pip2[4],Pim[4];
182 Km[0] = Km0.get(0); Pip1[0] = pi1.get(0); Pip2[0] = pi2.get(0); Pim[0] = pi3.get(0);
183 Km[1] = Km0.get(1); Pip1[1] = pi1.get(1); Pip2[1] = pi2.get(1); Pim[1] = pi3.get(1);
184 Km[2] = Km0.get(2); Pip1[2] = pi1.get(2); Pip2[2] = pi2.get(2); Pim[2] = pi3.get(2);
185 Km[3] = Km0.get(3); Pip1[3] = pi1.get(3); Pip2[3] = pi2.get(3); Pim[3] = pi3.get(3);
186 double Prob = calPDF(Km, Pip1, Pip2, Pim);
187 if(Prob>maxprob) {
188 maxprob=Prob;
189 }
190 }
191 */
192 p->initializePhaseSpace(getNDaug(), getDaugs());
193 EvtVector4R Km0 = p->getDaug(0)->getP4();
194 EvtVector4R pi1 = p->getDaug(1)->getP4();
195 EvtVector4R pi2 = p->getDaug(2)->getP4();
196 EvtVector4R pi3 = p->getDaug(3)->getP4();
197
198 double Km[4], Pip1[4], Pip2[4], Pim[4];
199 Km[0] = Km0.get(0); Pip1[0] = pi1.get(0); Pip2[0] = pi2.get(0); Pim[0] = pi3.get(0);
200 Km[1] = Km0.get(1); Pip1[1] = pi1.get(1); Pip2[1] = pi2.get(1); Pim[1] = pi3.get(1);
201 Km[2] = Km0.get(2); Pip1[2] = pi1.get(2); Pip2[2] = pi2.get(2); Pim[2] = pi3.get(2);
202 Km[3] = Km0.get(3); Pip1[3] = pi1.get(3); Pip2[3] = pi2.get(3); Pim[3] = pi3.get(3);
203 double prob = calPDF(Km, Pip1, Pip2, Pim);
204 setProb(prob);
205 return;
206 }
207
208 double EvtD0ToKpipipi::calPDF(double* Km, double* Pip1, double* Pip2, double* Pim)
209 {
210 Km[0] = sqrt(mass_Kaon * mass_Kaon + Km[1] * Km[1] + Km[2] * Km[2] + Km[3] * Km[3]);
211 Pip1[0] = sqrt(mass_Pion * mass_Pion + Pip1[1] * Pip1[1] + Pip1[2] * Pip1[2] + Pip1[3] * Pip1[3]);
212 Pip2[0] = sqrt(mass_Pion * mass_Pion + Pip2[1] * Pip2[1] + Pip2[2] * Pip2[2] + Pip2[3] * Pip2[3]);
213 Pim[0] = sqrt(mass_Pion * mass_Pion + Pim[1] * Pim[1] + Pim[2] * Pim[2] + Pim[3] * Pim[3]);
214
215 EvtComplex PDF[24];
216 int g[3];
217 g[0] = 1; g[1] = 1;
218
219 //-----------D->K*rho--------
220 g[2] = 0;
221 PDF[0] = D2VV(Km, Pip1, Pip2, Pim, g) + D2VV(Km, Pip2, Pip1, Pim, g);
222 g[2] = 1;
223 PDF[1] = D2VV(Km, Pip1, Pip2, Pim, g) + D2VV(Km, Pip2, Pip1, Pim, g);
224 g[2] = 2;
225 PDF[2] = D2VV(Km, Pip1, Pip2, Pim, g) + D2VV(Km, Pip2, Pip1, Pim, g);
226 //-----------D->K*rho finish--
227 //----------D->a1K------------
228 g[2] = 0;
229 PDF[3] = D2AP_A2VP(Km, Pip2, Pip1, Pim, g, 2) + D2AP_A2VP(Km, Pip1, Pip2, Pim, g, 2);
230 g[2] = 2;
231 PDF[4] = D2AP_A2VP(Km, Pip2, Pip1, Pim, g, 2) + D2AP_A2VP(Km, Pip1, Pip2, Pim, g, 2);
232 //----------D->a1K finish-----
233 //--D->K1Pi--K1->K*Pi---------
234 g[2] = 0;
235 PDF[5] = D2AP_A2VP(Pip2, Pim, Km, Pip1, g, 1) + D2AP_A2VP(Pip1, Pim, Km, Pip2, g, 1);
236 g[2] = 2;
237 PDF[6] = D2AP_A2VP(Pip2, Pim, Km, Pip1, g, 1) + D2AP_A2VP(Pip1, Pim, Km, Pip2, g, 1);
238 //--D->K1Pi--K1->K*Pi-finish--
239 //--D->K1Pi--K1->rhoK---------
240 g[2] = 0;
241 PDF[7] = D2AP_A2VP(Pip2, Km, Pip1, Pim, g, 3) + D2AP_A2VP(Pip1, Km, Pip2, Pim, g, 3);
242 //--D->K1Pi--K1->rhoK-finish--
243 //--D->K1Pi--K1->K*0(1430)Pi--
244 PDF[8] = D2AP_A2SP(Pip2, Pim, Km, Pip1, 1) + D2AP_A2SP(Pip1, Pim, Km, Pip2, 1);
245 //--------res finish----------
246 //----------------non-res------------------
247 //--------KPiSrho------------------
248 PDF[9] = D2VS(Pip1, Pim, Km, Pip2, 1, 2) + D2VS(Pip2, Pim, Km, Pip1, 1, 2);
249 //--------K*PiPiS-----------------
250 PDF[10] = D2VS(Km, Pip1, Pip2, Pim, 1, 1) + D2VS(Km, Pip2, Pip1, Pim, 1, 1);
251 //--------K*PiP-------------------
252 PDF[11] = D2PP_P2VP(Pip2, Pim, Km, Pip1, 1) + D2PP_P2VP(Pip1, Pim, Km, Pip2, 1);
253 //--------rhoKA--------------------
254 g[0] = 1; g[1] = 0;
255 g[2] = 0;
256 PDF[12] = 0;
257 g[2] = 2;
258 PDF[13] = D2AP_A2VP(Pip2, Km, Pip1, Pim, g, 3) + D2AP_A2VP(Pip1, Km, Pip2, Pim, g, 3);
259 //-------PHSP-----------------------
260 PDF[14] = PHSP(Km, Pip1) + PHSP(Km, Pip2);
261 //------KPiVPiPiV-----------------------
262 g[0] = 0; g[1] = 0; g[2] = 0;
263 PDF[15] = D2VV(Km, Pip1, Pip2, Pim, g) + D2VV(Km, Pip2, Pip1, Pim, g);
264 //------KPiVPiPiS----------------------
265 PDF[16] = D2VS(Km, Pip1, Pip2, Pim, 0, 1) + D2VS(Km, Pip2, Pip1, Pim, 0, 1);
266 //------KPiSPiPiV----------------------
267 PDF[17] = D2VS(Pip1, Pim, Km, Pip2, 0, 2) + D2VS(Pip2, Pim, Km, Pip1, 0, 2);
268 //-------------------------------------
269 PDF[18] = D2PP_P2VP(Pip2, Km, Pip1, Pim, 2) + D2PP_P2VP(Pip1, Km, Pip2, Pim, 2);
270 //-------------------------------------
271 PDF[19] = D2VP_V2VP(Pip2, Pim, Km, Pip1, 1) + D2VP_V2VP(Pip1, Pim, Km, Pip2, 1);
272 //-------------------------------------
273 PDF[20] = D2VP_V2VP(Pip2, Km, Pip1, Pim, 2) + D2VP_V2VP(Pip1, Km, Pip2, Pim, 2);
274 //-------------------------------------
275 PDF[21] = D2TS(Km, Pip1, Pip2, Pim, 1) + D2TS(Km, Pip2, Pip1, Pim, 1);
276 PDF[22] = D2TS(Pip1, Pim, Km, Pip2, 2) + D2TS(Pip2, Pim, Km, Pip1, 2);
277 PDF[23] = D2AP_A2SP(Km, Pip2, Pip1, Pim, 2) + D2AP_A2SP(Km, Pip1, Pip2, Pim, 2);
278//------------------------------------------
279 EvtComplex cof;
280 EvtComplex pdf, module;
281 pdf = EvtComplex(0, 0);
282 for (int i = 0; i != 24; i++) {
283 cof = EvtComplex(rho[i] * cos(phi[i]), rho[i] * sin(phi[i]));
284 pdf = pdf + cof * PDF[i];
285 }
286 module = conj(pdf) * pdf;
287 double value;
288 value = real(module);
289 return (value <= 0) ? 1e-20 : value;
290 }
291
292 EvtComplex EvtD0ToKpipipi::KPiSFormfactor(double sa, double sb, double sc, double r)
293 {
294 (void)r;
295 double m1430 = 1.463;
296 double sa0 = m1430 * m1430;
297 double w1430 = 0.233;
298 double q0 = (sa0 + sb - sc) * (sa0 + sb - sc) / (4 * sa0) - sb;
299 if (q0 < 0) q0 = 1e-16;
300 double qs = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
301 double q = sqrt(qs);
302 double Width = w1430 * q * m1430 / sqrt(sa * q0);
303 double temp_R = atan(m1430 * Width / (sa0 - sa));
304 if (temp_R < 0) temp_R += math_pi;
305 double deltaR = -5.31 + temp_R;
306 double temp_F = atan(2 * 1.07 * q / (2 + (-1.8) * 1.07 * qs));
307 if (temp_F < 0) temp_F += math_pi;
308 double deltaF = 2.33 + temp_F;
309 EvtComplex cR(cos(deltaR), sin(deltaR));
310 EvtComplex cF(cos(deltaF), sin(deltaF));
311 EvtComplex amp = 0.8 * sin(deltaF) * cF + sin(deltaR) * cR * cF * cF;
312 return amp;
313 }
314
315 EvtComplex EvtD0ToKpipipi::D2VV(const double* P1, const double* P2, const double* P3, const double* P4, int* g)
316 {
317 double t1V1[4], t1V2[4], t1D[4], t2D[4][4];
318 double temp_PDF = 0;
319 EvtComplex amp_PDF(0, 0);
320 EvtComplex pro[2];
321 //---------use g[0],g[1] to define res or non-res, g[2] represent S, P or D wave
322 double sa[3], sb[3], sc[3], B[3];
323 double pV1[4], pV2[4], pD[4];
324 for (int i = 0; i != 4; i++) {
325 pV1[i] = P1[i] + P2[i];
326 pV2[i] = P3[i] + P4[i];
327 pD[i] = pV1[i] + pV2[i];
328 }
329 sa[0] = dot(pV1, pV1);
330 sb[0] = dot(P1, P1);
331 sc[0] = dot(P2, P2);
332 sa[1] = dot(pV2, pV2);
333 sb[1] = dot(P3, P3);
334 sc[1] = dot(P4, P4);
335 sa[2] = dot(pD, pD);
336 sb[2] = sa[0];
337 sc[2] = sa[1];
338 if (g[1] == 1) {
339 pro[1] = propagatorGS(mass[3], width[3], sa[1], sb[1], sc[1], rRes, 1);
340 }
341 if (g[0] == 1) {
342 pro[0] = propagatorRBW(mass[1], width[1], sa[0], sb[0], sc[0], rRes, 1);
343 }
344 if (g[0] == 0) pro[0] = 1;
345 if (g[1] == 0) pro[1] = 1;
346 B[0] = barrier(1, sa[0], sb[0], sc[0], rRes);
347 B[1] = barrier(1, sa[1], sb[1], sc[1], rRes);
348 calt1(P1, P2, t1V1);
349 calt1(P3, P4, t1V2);
350 if (g[2] == 0) {
351 for (int i = 0; i != 4; i++) {
352 temp_PDF += (G[i][i]) * t1V1[i] * t1V2[i];
353 }
354 B[2] = 1;
355 }
356 if (g[2] == 1) {
357 calt1(pV1, pV2, t1D);
358 for (int i = 0; i != 4; i++) {
359 for (int j = 0; j != 4; j++) {
360 for (int k = 0; k != 4; k++) {
361 for (int l = 0; l != 4; l++) {
362 temp_PDF += E[i][j][k][l] * pD[i] * t1D[j] * t1V1[k] * t1V2[l] * (G[i][i]) * (G[j][j]) * (G[l][l]) * (G[k][k]);
363 }
364 }
365 }
366 }
367 B[2] = barrier(1, sa[2], sb[2], sc[2], rD);
368 }
369 if (g[2] == 2) {
370 calt2(pV1, pV2, t2D);
371 for (int i = 0; i != 4; i++) {
372 for (int j = 0; j != 4; j++) {
373 temp_PDF += t2D[i][j] * t1V1[i] * t1V2[j] * (G[i][i]) * (G[j][j]);
374 }
375 }
376 B[2] = barrier(2, sa[2], sb[2], sc[2], rD);
377 }
378 amp_PDF = temp_PDF * B[0] * B[1] * B[2] * pro[0] * pro[1];
379 return amp_PDF;
380 }
381
382 EvtComplex EvtD0ToKpipipi::D2AP_A2VP(const double* P1, const double* P2, const double* P3, const double* P4, int* g, const int flag)
383 {
384 //flag = 1, V = K*, A = K1(1270); flag = 2, V = rho, A = a1(1260);
385 //flag = 3, V = rho, A = K1(1270);
386 double temp_PDF = 0;
387 EvtComplex amp_PDF(0, 0);
388 EvtComplex pro[2];
389 double t1V[4], t1D[4], t2A[4][4];
390 double sa[3], sb[3], sc[3], B[3];
391 double pV[4], pA[4], pD[4];
392 for (int i = 0; i != 4; i++) {
393 pV[i] = P3[i] + P4[i];
394 pA[i] = pV[i] + P2[i];
395 pD[i] = pA[i] + P1[i];
396 }
397 sa[0] = dot(pV, pV);
398 sb[0] = dot(P3, P3);
399 sc[0] = dot(P4, P4);
400 sa[1] = dot(pA, pA);
401 sb[1] = sa[0];
402 sc[1] = dot(P2, P2);
403 sa[2] = dot(pD, pD);
404 sb[2] = sa[1];
405 sc[2] = dot(P1, P1);
406 if (g[0] == 1) {
407 if (flag == 1) pro[0] = propagatorRBW(mass[1], width[1], sa[0], sb[0], sc[0], rRes, 1);
408 if (flag == 2 || flag == 3) pro[0] = propagatorGS(mass[3], width[3], sa[0], sb[0], sc[0], rRes, 1);
409 }
410 if (g[1] == 1) {
411 if (flag == 1) pro[1] = propogator(mass[0], width[0], sa[1]);
412 if (flag == 2) pro[1] = propagatorRBW(mass[2], width[2], sa[1], sb[1], sc[1], rRes, g[2]);
413 if (flag == 3) pro[1] = propogator(mass[0], width[0], sa[1]);
414 }
415 if (g[0] == 0) pro[0] = 1;
416 if (g[1] == 0) pro[1] = 1;
417 B[0] = barrier(1, sa[0], sb[0], sc[0], rRes);
418 B[2] = barrier(1, sa[2], sb[2], sc[2], rD);
419 calt1(P3, P4, t1V);
420 calt1(pA, P1, t1D);
421 if (g[2] == 0) {
422 for (int i = 0; i != 4; i++) {
423 for (int j = 0; j != 4; j++) {
424 temp_PDF += t1D[i] * (G[i][j] - pA[i] * pA[j] / sa[1]) * t1V[j] * (G[i][i]) * (G[j][j]);
425 }
426 }
427 B[1] = 1;
428 }
429 if (g[2] == 2) {
430 calt2(pV, P2, t2A);
431 for (int i = 0; i != 4; i++) {
432 for (int j = 0; j != 4; j++) {
433 temp_PDF += t1D[i] * t2A[i][j] * t1V[j] * (G[i][i]) * (G[j][j]);
434 }
435 }
436 B[1] = barrier(2, sa[1], sb[1], sc[1], rRes);
437 }
438 amp_PDF = temp_PDF * B[0] * B[1] * B[2] * pro[0] * pro[1];
439 return amp_PDF;
440 }
441 EvtComplex EvtD0ToKpipipi::D2AP_A2SP(const double* P1, const double* P2, const double* P3, const double* P4, const int flag)
442 {
443 //flag = 1, S = K*; flag = 2, S = rho
444 double temp_PDF = 0;
445 EvtComplex amp_PDF(0, 0);
446 EvtComplex pro;
447 double sa[3], sb[3], sc[3], B[3];
448 double t1D[4], t1A[4];
449 double pS[4], pA[4], pD[4];
450 for (int i = 0; i != 4; i++) {
451 pS[i] = P3[i] + P4[i];
452 pA[i] = pS[i] + P2[i];
453 pD[i] = pA[i] + P1[i];
454 }
455 sa[0] = dot(pS, pS);
456 sb[0] = dot(P3, P3);
457 sc[0] = dot(P4, P4);
458 sa[1] = dot(pA, pA);
459 sb[1] = sa[0];
460 sc[1] = dot(P2, P2);
461 sa[2] = dot(pD, pD);
462 sb[2] = sa[1];
463 sc[2] = dot(P1, P1);
464 B[1] = barrier(1, sa[1], sb[1], sc[1], rRes);
465 B[2] = barrier(1, sa[2], sb[2], sc[2], rD);
466 calt1(pA, P1, t1D);
467 calt1(pS, P2, t1A);
468 for (int i = 0; i != 4; i++) {
469 temp_PDF += t1D[i] * t1A[i] * (G[i][i]);
470 }
471 amp_PDF = temp_PDF * B[1] * B[2];
472 if (flag == 1) amp_PDF *= KPiSFormfactor(sa[0], sb[0], sc[0], rRes);
473 return amp_PDF;
474 }
475 EvtComplex EvtD0ToKpipipi::D2PP_P2VP(const double* P1, const double* P2, const double* P3, const double* P4, const int flag)
476 {
477 //flag = 1, V = K*; flag = 2, V = rho
478 double temp_PDF = 0;
479 EvtComplex amp_PDF(0, 0);
480 EvtComplex pro;
481 double sa[3], sb[3], sc[3], B[3];
482 double t1V[4];
483 double pV[4], pP[4], pD[4];
484 for (int i = 0; i != 4; i++) {
485 pV[i] = P3[i] + P4[i];
486 pP[i] = pV[i] + P2[i];
487 pD[i] = pP[i] + P1[i];
488 }
489 sa[0] = dot(pV, pV);
490 sb[0] = dot(P3, P3);
491 sc[0] = dot(P4, P4);
492 sa[1] = dot(pP, pP);
493 sb[1] = sa[0];
494 sc[1] = dot(P2, P2);
495 sa[2] = dot(pD, pD);
496 sb[2] = sa[1];
497 sc[2] = dot(P1, P1);
498 B[0] = barrier(1, sa[0], sb[0], sc[0], rRes);
499 B[1] = barrier(1, sa[1], sb[1], sc[1], rRes);
500 if (flag == 1) pro = propagatorRBW(mass[1], width[1], sa[0], sb[0], sc[0], rRes, 1);
501 if (flag == 2) pro = propagatorGS(mass[3], width[3], sa[0], sb[0], sc[0], rRes, 1);
502 calt1(P3, P4, t1V);
503 for (int i = 0; i != 4; i++) {
504 temp_PDF += P2[i] * t1V[i] * (G[i][i]);
505 }
506 amp_PDF = temp_PDF * B[0] * B[1] * pro;
507 return amp_PDF;
508 }
509 EvtComplex EvtD0ToKpipipi::D2VP_V2VP(const double* P1, const double* P2, const double* P3, const double* P4, const int flag)
510 {
511 //flag = 1, (K*Pi)V; flag = 2, (rhoK)V
512 double temp_PDF = 0;
513 EvtComplex amp_PDF(0, 0);
514 EvtComplex pro;
515 double sa[3], sb[3], sc[3], B[3];
516 double pV1[4], pV2[4], qV1[4], qV2[4], pD[4];
517 for (int i = 0; i != 4; i++) {
518 pV2[i] = P3[i] + P4[i];
519 qV2[i] = P3[i] - P4[i];
520 pV1[i] = pV2[i] + P2[i];
521 qV1[i] = pV2[i] - P2[i];
522 pD[i] = pV1[i] + P1[i];
523 }
524 for (int i = 0; i != 4; i++) {
525 for (int j = 0; j != 4; j++) {
526 for (int k = 0; k != 4; k++) {
527 for (int l = 0; l != 4; l++) {
528 temp_PDF += E[i][j][k][l] * pV1[i] * qV1[j] * P1[k] * qV2[l] * (G[i][i]) * (G[j][j]) * (G[k][k]) * (G[l][l]);
529 }
530 }
531 }
532 }
533 sa[0] = dot(pV2, pV2);
534 sb[0] = dot(P3, P3);
535 sc[0] = dot(P4, P4);
536 sa[1] = dot(pV1, pV1);
537 sb[1] = sa[0];
538 sc[1] = dot(P2, P2);
539 sa[2] = dot(pD, pD);
540 sb[2] = sa[1];
541 sc[2] = dot(P1, P1);
542 if (flag == 1) pro = propagatorRBW(mass[1], width[1], sa[0], sb[0], sc[0], rRes, 1);
543 if (flag == 2) pro = propagatorGS(mass[3], width[3], sa[0], sb[0], sc[0], rRes, 1);
544 B[0] = barrier(1, sa[0], sb[0], sc[0], rRes);
545 B[1] = barrier(1, sa[1], sb[1], sc[1], rRes);
546 B[2] = barrier(1, sa[2], sb[2], sc[2], rD);
547 amp_PDF = temp_PDF * B[0] * B[1] * B[2] * pro;
548 return amp_PDF;
549 }
550 EvtComplex EvtD0ToKpipipi::D2VS(const double* P1, const double* P2, const double* P3, const double* P4, int g, const int flag)
551 {
552 //flag = 1, V = K*; flag = 2, V = rho
553 double temp_PDF = 0;
554 EvtComplex amp_PDF(0, 0);
555 EvtComplex pro;
556 double sa[3], sb[3], sc[3], B[3];
557 double t1D[4], t1V[4];
558 double pS[4], pV[4], pD[4];
559 for (int i = 0; i != 4; i++) {
560 pS[i] = P3[i] + P4[i];
561 pV[i] = P1[i] + P2[i];
562 pD[i] = pS[i] + pV[i];
563 }
564 sa[0] = dot(pS, pS);
565 sb[0] = dot(P3, P3);
566 sc[0] = dot(P4, P4);
567 sa[1] = dot(pV, pV);
568 sb[1] = dot(P1, P1);
569 sc[1] = dot(P2, P2);
570 sa[2] = dot(pD, pD);
571 sb[2] = sa[0];
572 sc[2] = sa[1];
573 if (g == 1) {
574 if (flag == 2) pro = propagatorGS(mass[3], width[3], sa[1], sb[1], sc[1], rRes, 1);
575 if (flag == 1) pro = propagatorRBW(mass[1], width[1], sa[1], sb[1], sc[1], rRes, 1);
576 }
577 if (g == 0) pro = 1;
578 B[1] = barrier(1, sa[1], sb[1], sc[1], rRes);
579 B[2] = barrier(1, sa[2], sb[2], sc[2], rD);
580 calt1(P1, P2, t1V);
581 calt1(pS, pV, t1D);
582 for (int i = 0; i != 4; i++) {
583 temp_PDF += G[i][i] * t1D[i] * t1V[i];
584 }
585 amp_PDF = temp_PDF * B[1] * B[2] * pro;
586 if (flag == 2) amp_PDF *= KPiSFormfactor(sa[0], sb[0], sc[0], rRes);
587 return amp_PDF;
588 }
589 EvtComplex EvtD0ToKpipipi::D2TS(const double* P1, const double* P2, const double* P3, const double* P4, const int flag)
590 {
591 //flag = 1, T = K*; flag = 2, T = rho
592 double temp_PDF = 0;
593 EvtComplex amp_PDF(0, 0);
594 double sa[3], sb[3], sc[3], B[3];
595 double t2D[4][4], t2T[4][4];
596 double pS[4], pT[4], pD[4];
597 for (int i = 0; i != 4; i++) {
598 pS[i] = P3[i] + P4[i];
599 pT[i] = P1[i] + P2[i];
600 pD[i] = pT[i] + pS[i];
601 }
602 sa[0] = dot(pT, pT);
603 sb[0] = dot(P1, P1);
604 sc[0] = dot(P2, P2);
605 sa[1] = dot(pS, pS);
606 sb[1] = dot(P3, P3);
607 sc[1] = dot(P4, P4);
608 sa[2] = dot(pD, pD);
609 sb[2] = sa[0];
610 sc[2] = sa[1];
611 B[0] = barrier(2, sa[0], sb[0], sc[0], rRes);
612 B[2] = barrier(2, sa[2], sb[2], sc[2], rD);
613 calt2(P1, P2, t2T);
614 calt2(pT, pS, t2D);
615 for (int i = 0; i != 4; i++) {
616 for (int j = 0; j != 4; j++) {
617 temp_PDF += t2D[i][j] * t2T[j][i] * (G[i][i]) * (G[j][j]);
618 }
619 }
620 amp_PDF = temp_PDF * B[0] * B[2];
621 if (flag == 2) amp_PDF *= KPiSFormfactor(sa[1], sb[1], sc[1], rRes);
622 return amp_PDF;
623 }
624 EvtComplex EvtD0ToKpipipi::PHSP(double* Km, double* Pip)
625 {
626 EvtComplex amp_PDF(0, 0);
627 double sa, sb, sc;
628 double KPi[4];
629 for (int i = 0; i != 4; i++) {
630 KPi[i] = Km[i] + Pip[i];
631 }
632 sa = dot(KPi, KPi);
633 sb = dot(Km, Km);
634 sc = dot(Pip, Pip);
635 amp_PDF = KPiSFormfactor(sa, sb, sc, rRes);
636 return amp_PDF;
637 }
638 EvtComplex EvtD0ToKpipipi::propogator(double Mass, double Width, double sx)const
639 {
640 EvtComplex ci(0, 1);
641 EvtComplex prop = 1.0 / (Mass * Mass - sx - ci * Mass * Width);
642 return prop;
643 }
644 double EvtD0ToKpipipi::wid(double Mass, double sa, double sb, double sc, double r, int l)const
645 {
646 double widm(0.), q(0.), q0(0.);
647 double sa0 = Mass * Mass;
648 double m = sqrt(sa);
649 q = Qabcs(sa, sb, sc);
650 q0 = Qabcs(sa0, sb, sc);
651 double z, z0;
652 z = q * r * r;
653 z0 = q0 * r * r;
654 double t = q / q0;
655 double F(0.);
656 if (l == 0) F = 1;
657 if (l == 1) F = sqrt((1 + z0) / (1 + z));
658 if (l == 2) F = sqrt((9 + 3 * z0 + z0 * z0) / (9 + 3 * z + z * z));
659 widm = pow(t, l + 0.5) * Mass / m * F * F;
660 return widm;
661 }
662 double EvtD0ToKpipipi::h(double m, double q)const
663 {
664 double h(0.);
665 h = 2 / pi * q / m * log((m + 2 * q) / (2 * mpi));
666 return h;
667 }
668 double EvtD0ToKpipipi::dh(double Mass, double q0)const
669 {
670 double dh = h(Mass, q0) * (1.0 / (8 * q0 * q0) - 1.0 / (2 * Mass * Mass)) + 1.0 / (2 * pi * Mass * Mass);
671 return dh;
672 }
673 double EvtD0ToKpipipi::f(double Mass, double sx, double q0, double q)const
674 {
675 double m = sqrt(sx);
676 double f = Mass * Mass / (pow(q0, 3)) * (q * q * (h(m, q) - h(Mass, q0)) + (Mass * Mass - sx) * q0 * q0 * dh(Mass, q0));
677 return f;
678 }
679 double EvtD0ToKpipipi::d(double Mass, double q0)const
680 {
681 double d = 3.0 / pi * mpi * mpi / (q0 * q0) * log((Mass + 2 * q0) / (2 * mpi)) + Mass / (2 * pi * q0) - (mpi * mpi * Mass) /
682 (pi * pow(q0, 3));
683 return d;
684 }
685 EvtComplex EvtD0ToKpipipi::propagatorRBW(double Mass, double Width, double sa, double sb, double sc, double r, int l)const
686 {
687 EvtComplex ci(0, 1);
688 EvtComplex prop = 1.0 / (Mass * Mass - sa - ci * Mass * Width * wid(Mass, sa, sb, sc, r, l));
689 return prop;
690 }
691 EvtComplex EvtD0ToKpipipi::propagatorGS(double Mass, double Width, double sa, double sb, double sc, double r, int l)const
692 {
693 EvtComplex ci(0, 1);
694 double q = Qabcs(sa, sb, sc);
695 double sa0 = Mass * Mass;
696 double q0 = Qabcs(sa0, sb, sc);
697 q = sqrt(q);
698 q0 = sqrt(q0);
699 EvtComplex prop = (1 + d(Mass, q0) * Width / Mass) / (Mass * Mass - sa + Width * f(Mass, sa, q0, q) - ci * Mass * Width * wid(Mass,
700 sa, sb, sc, r, l));
701 return prop;
702 }
703 double EvtD0ToKpipipi::dot(const double* a1, const double* a2)const
704 {
705 double dot = 0;
706 for (int i = 0; i != 4; i++) {
707 dot += a1[i] * a2[i] * G[i][i];
708 }
709 return dot;
710 }
711 double EvtD0ToKpipipi::Qabcs(double sa, double sb, double sc)const
712 {
713 double Qabcs = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
714 if (Qabcs < 0) Qabcs = 1e-16;
715 return Qabcs;
716 }
717 double EvtD0ToKpipipi::barrier(double l, double sa, double sb, double sc, double r)const
718 {
719 double q = Qabcs(sa, sb, sc);
720 q = sqrt(q);
721 double z = q * r;
722 z = z * z;
723 double F = 1;
724 if (l > 2) F = 0;
725 if (l == 0)
726 F = 1;
727 if (l == 1) {
728 F = sqrt((2 * z) / (1 + z));
729 }
730 if (l == 2) {
731 F = sqrt((13 * z * z) / (9 + 3 * z + z * z));
732 }
733 return F;
734 }
735 void EvtD0ToKpipipi::calt1(const double* daug1, const double* daug2, double* t1) const
736 {
737 double p, pq;
738 double pa[4], qa[4];
739 for (int i = 0; i != 4; i++) {
740 pa[i] = daug1[i] + daug2[i];
741 qa[i] = daug1[i] - daug2[i];
742 }
743 p = dot(pa, pa);
744 pq = dot(pa, qa);
745 for (int i = 0; i != 4; i++) {
746 t1[i] = qa[i] - pq / p * pa[i];
747 }
748 }
749 void EvtD0ToKpipipi::calt2(const double* daug1, const double* daug2, double (*t2)[4]) const
750 {
751 double p, r;
752 double pa[4], t1[4];
753 calt1(daug1, daug2, t1);
754 r = dot(t1, t1);
755 for (int i = 0; i != 4; i++) {
756 pa[i] = daug1[i] + daug2[i];
757 }
758 p = dot(pa, pa);
759 for (int i = 0; i != 4; i++) {
760 for (int j = 0; j != 4; j++) {
761 t2[i][j] = t1[i] * t1[j] - 1.0 / 3 * r * (G[i][j] - pa[i] * pa[j] / p);
762 }
763 }
764 }
765
767} // Belle 2 namespace
Belle2::EvtD0ToKpipipi::mD
double mD
Fixed parameters.
Definition EvtD0ToKpipipi.h:92
Belle2::EvtD0ToKpipipi::d
double d(double mass, const double q0) const
d function in Gounaris-Sakurai lineshape
Definition EvtD0ToKpipipi.cc:679
Belle2::EvtD0ToKpipipi::propagatorRBW
EvtComplex propagatorRBW(double mass, double width, const double sa, const double sb, const double sc, const double r, const int l) const
Relativistic Breit-Wigner Lineshape Function.
Definition EvtD0ToKpipipi.cc:685
Belle2::EvtD0ToKpipipi::dh
double dh(double mass, const double q0) const
derivative h function in Gounaris-Sakurai lineshape
Definition EvtD0ToKpipipi.cc:668
Belle2::EvtD0ToKpipipi::calPDF
double calPDF(double *Km, double *Pip1, double *Pip2, double *Pim)
Probability distribution function of the decay.
Definition EvtD0ToKpipipi.cc:208
Belle2::EvtD0ToKpipipi::D2VV
EvtComplex D2VV(const double *P1, const double *P2, const double *P3, const double *P4, int *g)
Amplitude modes.
Definition EvtD0ToKpipipi.cc:315
Belle2::EvtD0ToKpipipi::f
double f(double mass, const double sx, const double q0, const double q) const
f function in Gounaris-Sakurai lineshape
Definition EvtD0ToKpipipi.cc:673
Belle2::EvtD0ToKpipipi::propagatorGS
EvtComplex propagatorGS(double mass, double width, const double sa, const double sb, const double sc, const double r, const int l) const
Gounaris-Sakurai lineshape Function.
Definition EvtD0ToKpipipi.cc:691
Belle2::EvtD0ToKpipipi::dot
double dot(const double *a1, const double *a2) const
Four-Vector Scalar Product.
Definition EvtD0ToKpipipi.cc:703
Belle2::EvtD0ToKpipipi::calt1
void calt1(const double *daug1, const double *daug2, double *t1) const
Covariant Spin-1 Projector.
Definition EvtD0ToKpipipi.cc:735
Belle2::EvtD0ToKpipipi::calt2
void calt2(const double *daug1, const double *daug2, double(*t2)[4]) const
Covariant Spin-2 Projector.
Definition EvtD0ToKpipipi.cc:749
Belle2::EvtD0ToKpipipi::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 EvtD0ToKpipipi.cc:711
Belle2::EvtD0ToKpipipi::KPiSFormfactor
EvtComplex KPiSFormfactor(const double sa, const double sb, const double sc, const double r)
K pi S-wave form factor.
Definition EvtD0ToKpipipi.cc:292
Belle2::EvtD0ToKpipipi::barrier
double barrier(const double l, const double sa, const double sb, const double sc, const double r) const
Blatt-Weisskopf barrier factors.
Definition EvtD0ToKpipipi.cc:717
Belle2::atan
double atan(double a)
atan for double
Definition beamHelpers.h:34
B2_EVTGEN_REGISTER_MODEL
#define B2_EVTGEN_REGISTER_MODEL(classname)
Class to register B2_EVTGEN_REGISTER_MODEL.
Definition EvtGenModelRegister.h:75
Belle2::EvtD0ToKpipipi::wid
double wid(double mass, const double sa, const double sb, const double sc, const double r, const int l) const
Energy dependent width.
Definition EvtD0ToKpipipi.cc:644
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
Belle2::EvtD0ToKpipipi::propogator
EvtComplex propogator(double mass, double width, const double sx) const
Relativistic Breit-Wigner Lineshape Function (Fixed Width)
Definition EvtD0ToKpipipi.cc:638
Belle2::EvtD0ToKpipipi::h
double h(const double m, const double q) const
h function in Gounaris-Sakurai lineshape
Definition EvtD0ToKpipipi.cc:662
Belle2
Abstract base class for different kinds of events.
Definition MillepedeAlgorithm.h:17