Belle II Software development
EvtD0ToKpipi0pi0.cc
1// Model: EvtD0ToKpipi0pi0
2// This file is an amplitude model for D0 -> K- pi+ pi0 pi0.
3// The model is from the BESIII Collaboration in PRD 99, 092008 (2019). DOI:  https://doi.org/10.1103/PhysRevD.99.092008
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/EvtD0ToKpipi0pi0.h>
20
21namespace Belle2 {
26
28 B2_EVTGEN_REGISTER_MODEL(EvtD0ToKpipi0pi0);
29
30 EvtD0ToKpipi0pi0::~EvtD0ToKpipi0pi0() {}
31
32 std::string EvtD0ToKpipi0pi0::getName()
33 {
34 return "D0ToKpipi0pi0";
35 }
36
37 EvtDecayBase* EvtD0ToKpipi0pi0::clone()
38 {
39 return new EvtD0ToKpipi0pi0;
40 }
41
42 void EvtD0ToKpipi0pi0::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 mod[0] = 1; rho[0] = 2.02; phi[0] = -0.75;
55 mod[1] = 1; rho[1] = 1.66; phi[1] = -2.90;
56 mod[2] = 0; rho[2] = 0; phi[2] = 0;
57 mod[3] = 0; rho[3] = 0; phi[3] = 0;
58 mod[4] = 0; rho[4] = 0; phi[4] = 0;
59 mod[5] = 0; rho[5] = 0; phi[5] = 0;
60 mod[6] = 0; rho[6] = 0; phi[6] = 0;
61 mod[7] = 1; rho[7] = 1; phi[7] = 0;
62 mod[8] = 1; rho[8] = 0.842; phi[8] = -2.05;
63 mod[9] = 1; rho[9] = 0.0218; phi[9] = 1.84;
64 mod[10] = 0; rho[10] = 0; phi[10] = 0;
65 mod[11] = 0; rho[11] = 0; phi[11] = 0;
66 mod[12] = 0; rho[12] = 0; phi[12] = 0;
67 mod[13] = 1; rho[13] = 0.0336; phi[13] = -1.55;
68 mod[14] = 1; rho[14] = 0.109; phi[14] = -1.35;
69 mod[15] = 1; rho[15] = 0.196; phi[15] = -2.07;
70 mod[16] = 0; rho[16] = 0; phi[16] = 0;
71 mod[17] = 0; rho[17] = 0; phi[17] = 0;
72 mod[18] = 0; rho[18] = 0; phi[18] = 0;
73 mod[19] = 1; rho[19] = 0.363; phi[19] = 1.93;
74 mod[20] = 0; rho[20] = 0; phi[20] = 0;
75 mod[21] = 0; rho[21] = 0; phi[21] = 0;
76 mod[22] = 0; rho[22] = 0; phi[22] = 0;
77 mod[23] = 1; rho[23] = 0.555; phi[23] = 0.44;
78 mod[24] = 1; rho[24] = 0.526; phi[24] = -1.84;
79 mod[25] = 0; rho[25] = 0; phi[25] = 0;
80 mod[26] = 1; rho[26] = 1; phi[26] = 0.64;
81 mod[27] = 0; rho[27] = 0; phi[27] = 0;
82 mod[28] = 0; rho[28] = 0; phi[28] = 0;
83 mod[29] = 1; rho[29] = 3.34; phi[29] = -0.02;
84 mod[30] = 0; rho[30] = 0; phi[30] = 0;
85 mod[31] = 0; rho[31] = 0; phi[31] = 0;
86 mod[32] = 0; rho[32] = 0; phi[32] = 0;
87 mod[33] = 0; rho[33] = 0; phi[33] = 0;
88 mod[34] = 1; rho[34] = 1.76; phi[34] = -2.39;
89 mod[35] = 1; rho[35] = 0.175; phi[35] = 1.59;
90 mod[36] = 1; rho[36] = 0.397; phi[36] = 1.45;
91 mod[37] = 0; rho[37] = 0; phi[37] = 0;
92 mod[38] = 0; rho[38] = 0; phi[38] = 0;
93 mod[39] = 0; rho[39] = 0; phi[39] = 0;
94 mod[40] = 0; rho[40] = 0; phi[40] = 0;
95 mod[41] = 1; rho[41] = 1.02; phi[41] = 0.52;
96 mod[42] = 0; rho[42] = 0; phi[42] = 0;
97 mod[43] = 0; rho[43] = 0; phi[43] = 0;
98 mod[44] = 0; rho[44] = 0; phi[44] = 0;
99 mod[45] = 0; rho[45] = 0; phi[45] = 0;
100 mod[46] = 1; rho[46] = 0.146; phi[46] = 1.24;
101 mod[47] = 1; rho[47] = 0.0978; phi[47] = -2.89;
102 mod[48] = 1; rho[48] = 0.233; phi[48] = 2.41;
103 mod[49] = 0; rho[49] = 0; phi[49] = 0;
104 mod[50] = 1; rho[50] = 0.424; phi[50] = -0.94;
105 mod[51] = 1; rho[51] = 1.03; phi[51] = -1.93;
106 mod[52] = 0; rho[52] = 0; phi[52] = 0;
107 mod[53] = 0; rho[53] = 0; phi[53] = 0;
108 mod[54] = 1; rho[54] = 0.474; phi[54] = -1.17;
109 mod[55] = 0; rho[55] = 0; phi[55] = 0;
110 mod[56] = 0; rho[56] = 0; phi[56] = 0;
111 mod[57] = 0; rho[57] = 0; phi[57] = 0;
112 mod[58] = 0; rho[58] = 0; phi[58] = 0;
113 mod[59] = 0; rho[59] = 0; phi[59] = 0;
114 mod[60] = 0; rho[60] = 0; phi[60] = 0;
115 mod[61] = 1; rho[61] = 6.74; phi[61] = -1.74;
116 mod[62] = 0; rho[62] = 0; phi[62] = 0;
117 mod[63] = 0; rho[63] = 0; phi[63] = 0;
118 mod[64] = 0; rho[64] = 0; phi[64] = 0;
119 mod[65] = 0; rho[65] = 0; phi[65] = 0;
120 mod[66] = 1; rho[66] = 1.54; phi[66] = -2.93;
121 mod[67] = 1; rho[67] = 1.36; phi[67] = 2.23;
122
123 mass[0] = 1.23; width[0] = 0.50204;
124 mass[1] = 1.2723; width[1] = 0.09;
125 mass[2] = 0.89166; width[2] = 0.0508;
126 mass[3] = 0.89581; width[3] = 0.0474;
127 mass[4] = 0.77511; width[4] = 0.1491;
128
129 mD = 1.86486;
130 rRes = 3.0;
131 rD = 5.0;
132 metap = 0.95778;
133 mk0 = 0.497614;
134 mass_Kaon = 0.49368;
135 mass_Pion = 0.13957;
136 math_pi = 3.1415926;
137 mass_Pi0 = 0.1349766;
138 mkstrm = 0.89594;
139 mkstr0 = 0.89594;
140
141 pi = 3.1415926;
142 mpi = 0.13957;
143 g1 = 0.5468;
144 g2 = 0.23;
145
146 int GG[4][4] = { {1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1} };
147 int EE[4][4][4][4] = {
148 { {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0} },
149 {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 1}, {0, 0, -1, 0}},
150 {{0, 0, 0, 0}, {0, 0, 0, -1}, {0, 0, 0, 0}, {0, 1, 0, 0} },
151 {{0, 0, 0, 0}, {0, 0, 1, 0}, {0, -1, 0, 0}, {0, 0, 0, 0} }
152 },
153 { {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, -1}, {0, 0, 1, 0} },
154 {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0} },
155 {{0, 0, 0, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}, {-1, 0, 0, 0}},
156 {{0, 0, -1, 0}, {0, 0, 0, 0}, {1, 0, 0, 0}, {0, 0, 0, 0} }
157 },
158 { {{0, 0, 0, 0}, {0, 0, 0, 1}, {0, 0, 0, 0}, {0, -1, 0, 0}},
159 {{0, 0, 0, -1}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 0, 0, 0} },
160 {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0} },
161 {{0, 1, 0, 0}, {-1, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0} }
162 },
163 { {{0, 0, 0, 0}, {0, 0, -1, 0}, {0, 1, 0, 0}, {0, 0, 0, 0} },
164 {{0, 0, 1, 0}, {0, 0, 0, 0}, {-1, 0, 0, 0}, {0, 0, 0, 0} },
165 {{0, -1, 0, 0}, {1, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0} },
166 {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0} }
167 }
168 };
169 for (int i = 0; i < 4; i++) {
170 for (int j = 0; j < 4; j++) {
171 G[i][j] = GG[i][j];
172 for (int k = 0; k < 4; k++) {
173 for (int l = 0; l < 4; l++) {
174 E[i][j][k][l] = EE[i][j][k][l];
175 }
176 }
177 }
178 }
179 }
180
181 void EvtD0ToKpipi0pi0::initProbMax()
182 {
183 setProbMax(8060.0);
184 }
185
186 void EvtD0ToKpipi0pi0::decay(EvtParticle* p)
187 {
188 /*
189 double maxprob = 0.0;
190 for(int ir=0;ir<=60000000;ir++){
191 p->initializePhaseSpace(getNDaug(),getDaugs());
192 EvtVector4R Km0 = p->getDaug(0)->getP4();
193 EvtVector4R pi1 = p->getDaug(1)->getP4();
194 EvtVector4R pi2 = p->getDaug(2)->getP4();
195 EvtVector4R pi3 = p->getDaug(3)->getP4();
196 double Km[4],Pip[4],Pi01[4],Pi02[4];
197 Km[0] = Km0.get(0); Pip[0] = pi1.get(0); Pi01[0] = pi2.get(0); Pi02[0] = pi3.get(0);
198 Km[1] = Km0.get(1); Pip[1] = pi1.get(1); Pi01[1] = pi2.get(1); Pi02[1] = pi3.get(1);
199 Km[2] = Km0.get(2); Pip[2] = pi1.get(2); Pi01[2] = pi2.get(2); Pi02[2] = pi3.get(2);
200 Km[3] = Km0.get(3); Pip[3] = pi1.get(3); Pi01[3] = pi2.get(3); Pi02[3] = pi3.get(3);
201 double Prob = calPDF(Km, Pip, Pi01, Pi02);
202 if(Prob>maxprob) {
203 maxprob=Prob;
204 }
205 }
206 */
207 p->initializePhaseSpace(getNDaug(), getDaugs());
208 EvtVector4R Km0 = p->getDaug(0)->getP4();
209 EvtVector4R pi1 = p->getDaug(1)->getP4();
210 EvtVector4R pi2 = p->getDaug(2)->getP4();
211 EvtVector4R pi3 = p->getDaug(3)->getP4();
212
213 double Km[4], Pip[4], Pi01[4], Pi02[4];
214 Km[0] = Km0.get(0); Pip[0] = pi1.get(0); Pi01[0] = pi2.get(0); Pi02[0] = pi3.get(0);
215 Km[1] = Km0.get(1); Pip[1] = pi1.get(1); Pi01[1] = pi2.get(1); Pi02[1] = pi3.get(1);
216 Km[2] = Km0.get(2); Pip[2] = pi1.get(2); Pi01[2] = pi2.get(2); Pi02[2] = pi3.get(2);
217 Km[3] = Km0.get(3); Pip[3] = pi1.get(3); Pi01[3] = pi2.get(3); Pi02[3] = pi3.get(3);
218 double prob = PDF(Km, Pip, Pi01, Pi02);
219 setProb(prob);
220 return;
221 }
222
223 double EvtD0ToKpipi0pi0::PDF(double* Km, double* Pip, double* Pi01, double* Pi02)
224 {
225 Km[0] = sqrt(mass_Kaon * mass_Kaon + Km[1] * Km[1] + Km[2] * Km[2] + Km[3] * Km[3]);
226 Pip[0] = sqrt(mass_Pion * mass_Pion + Pip[1] * Pip[1] + Pip[2] * Pip[2] + Pip[3] * Pip[3]);
227 Pi01[0] = sqrt(mass_Pi0 * mass_Pi0 + Pi01[1] * Pi01[1] + Pi01[2] * Pi01[2] + Pi01[3] * Pi01[3]);
228 Pi02[0] = sqrt(mass_Pi0 * mass_Pi0 + Pi02[1] * Pi02[1] + Pi02[2] * Pi02[2] + Pi02[3] * Pi02[3]);
229
230 EvtComplex PDF[68];
231 int g[3];
232 //-------PHSP---------
233 PDF[0] = PHSP(Km, Pip) + PHSP(Km, Pip);
234 PDF[1] = PHSP(Km, Pi01) + PHSP(Km, Pi02);
235 //-------D2PP,P2VP------
236// PDF[2] = D2PP_P2VP(Pi01,Pi02,Km,Pip,0) + D2PP_P2VP(Pi02,Pi01,Km,Pip,0);
237// PDF[3] = D2PP_P2VP(Pi01,Pip,Km,Pi02,10) + D2PP_P2VP(Pi02,Pip,Km,Pi01,10);
238// PDF[4] = D2PP_P2VP(Pip,Pi01,Km,Pi02,20) + D2PP_P2VP(Pip,Pi02,Km,Pi01,20);
239// PDF[5] = D2PP_P2VP(Pi01,Km,Pip,Pi02,1) + D2PP_P2VP(Pi02,Km,Pip,Pi01,1);
240// PDF[6] = D2PP_P2VP(Km,Pi01,Pip,Pi02,11) + D2PP_P2VP(Km,Pi02,Pip,Pi01,11);
241 //----------D2AP,A2VP--------------
242 g[0] = 1; g[1] = 1; g[2] = 0;
243 PDF[7] = D2AP_A2VP(Km, Pi01, Pip, Pi02, g, 0) + D2AP_A2VP(Km, Pi02, Pip, Pi01, g, 0);
244 g[2] = 2;
245 PDF[8] = D2AP_A2VP(Km, Pi01, Pip, Pi02, g, 0) + D2AP_A2VP(Km, Pi02, Pip, Pi01, g, 0);
246 g[2] = 0;
247 PDF[9] = D2AP_A2VP(Pip, Pi01, Km, Pi02, g, 1) + D2AP_A2VP(Pip, Pi02, Km, Pi01, g, 1);
248// g[2] = 2;
249// PDF[10] = D2AP_A2VP(Pip,Pi01,Km,Pi02,g,1) + D2AP_A2VP(Pip,Pi02,Km,Pi01,g,1);
250// g[2] = 0;
251// PDF[11] = D2AP_A2VP(Pi01,Pip,Km,Pi02,g,21) + D2AP_A2VP(Pi02,Pip,Km,Pi01,g,21);
252// g[2] = 2;
253// PDF[12] = D2AP_A2VP(Pi01,Pip,Km,Pi02,g,21) + D2AP_A2VP(Pi02,Pip,Km,Pi01,g,21);
254 g[2] = 0;
255 PDF[13] = D2AP_A2VP(Pi01, Pi02, Km, Pip, g, 31) + D2AP_A2VP(Pi02, Pi01, Km, Pip, g, 31);
256 g[2] = 2;
257 PDF[14] = D2AP_A2VP(Pi01, Pi02, Km, Pip, g, 31) + D2AP_A2VP(Pi02, Pi01, Km, Pip, g, 31);
258 g[2] = 0;
259 PDF[15] = D2AP_A2VP(Pi01, Km, Pip, Pi02, g, 3) + D2AP_A2VP(Pi02, Km, Pip, Pi01, g, 3);
260// g[2] = 2;
261// PDF[16] = D2AP_A2VP(Pi01,Km,Pip,Pi02,g,3) + D2AP_A2VP(Pi02,Km,Pip,Pi01,g,3);
262 g[0] = 1; g[1] = 0; g[2] = 0;
263// PDF[17] = D2AP_A2VP(Km,Pi01,Pip,Pi02,g,0) + D2AP_A2VP(Km,Pi02,Pip,Pi01,g,0);
264// g[2] = 2;
265// PDF[18] = D2AP_A2VP(Km,Pi01,Pip,Pi02,g,0) + D2AP_A2VP(Km,Pi02,Pip,Pi01,g,0);
266 g[2] = 0;
267 PDF[19] = D2AP_A2VP(Pip, Pi01, Km, Pi02, g, 1) + D2AP_A2VP(Pip, Pi02, Km, Pi01, g, 1);
268// g[2] = 2;
269// PDF[20] = D2AP_A2VP(Pip,Pi01,Km,Pi02,g,1) + D2AP_A2VP(Pip,Pi02,Km,Pi01,g,1);
270// g[2] = 0;
271// PDF[21] = D2AP_A2VP(Pi01,Pip,Km,Pi02,g,21) + D2AP_A2VP(Pi02,Pip,Km,Pi01,g,21);
272// g[2] = 2;
273// PDF[22] = D2AP_A2VP(Pi01,Pip,Km,Pi02,g,21) + D2AP_A2VP(Pi02,Pip,Km,Pi01,g,21);
274 g[2] = 0;
275 PDF[23] = D2AP_A2VP(Pi01, Pi02, Km, Pip, g, 31) + D2AP_A2VP(Pi02, Pi01, Km, Pip, g, 31);
276 g[2] = 2;
277 PDF[24] = D2AP_A2VP(Pi01, Pi02, Km, Pip, g, 31) + D2AP_A2VP(Pi02, Pi01, Km, Pip, g, 31);
278// g[2] = 0;
279// PDF[25] = D2AP_A2VP(Pi01,Km,Pip,Pi02,g,3) + D2AP_A2VP(Pi02,Km,Pip,Pi01,g,3);
280 g[2] = 2;
281 PDF[26] = D2AP_A2VP(Pi01, Km, Pip, Pi02, g, 3) + D2AP_A2VP(Pi02, Km, Pip, Pi01, g, 3);
282 //--------D2AP,A2SP-----------------------------------
283// PDF[27] = D2AP_A2SP(Km,Pi01,Pip,Pi02,0) + D2AP_A2SP(Km,Pi02,Pip,Pi01,0);
284// PDF[28] = D2AP_A2SP(Km,Pip,Pi01,Pi02,10) + D2AP_A2SP(Km,Pip,Pi02,Pi01,10);
285 PDF[29] = D2AP_A2SP(Pi01, Pi02, Km, Pip, 1) + D2AP_A2SP(Pi02, Pi01, Km, Pip, 1);
286// PDF[30] = D2AP_A2SP(Pi01,Pip,Km,Pi02,11) + D2AP_A2SP(Pi02,Pip,Km,Pi01,11);
287// PDF[31] = D2AP_A2SP(Pip,Pi01,Km,Pi02,21) + D2AP_A2SP(Pip,Pi02,Km,Pi01,21);
288// PDF[32] = D2AP_A2SP(Pi01,Km,Pip,Pi02,31) + D2AP_A2SP(Pi02,Km,Pip,Pi01,31);
289// PDF[33] = D2AP_A2SP(Pip,Km,Pi01,Pi02,41) + D2AP_A2SP(Pip,Km,Pi02,Pi01,41);
290 //--------D2VS----------------------------
291 PDF[34] = D2VS(Pip, Pi01, Km, Pi02, 1, 0) + D2VS(Pip, Pi02, Km, Pi01, 1, 0);
292 PDF[35] = D2VS(Km, Pi01, Pip, Pi02, 1, 1) + D2VS(Km, Pi02, Pip, Pi01, 1, 1);
293 PDF[36] = D2VS(Km, Pip, Pi01, Pi02, 1, 11) + D2VS(Km, Pip, Pi02, Pi01, 1, 11);
294// PDF[37] = D2VS(Pi01,Pi02,Km,Pip,0,10) + D2VS(Pi02,Pi01,Km,Pip,0,10);
295// PDF[38] = D2VS(Pip,Pi01,Km,Pi02,0,0) + D2VS(Pip,Pi02,Km,Pi01,0,0);
296// PDF[39] = D2VS(Km,Pip,Pi01,Pi02,0,11) + D2VS(Km,Pip,Pi02,Pi01,0,11);
297// PDF[40] = D2VS(Km,Pi01,Pip,Pi02,0,1) + D2VS(Km,Pi02,Pip,Pi01,0,1);
298 //---------D2VP,V2VP----------------------
299 PDF[41] = D2VP_V2VP(Pi01, Pip, Km, Pi02, 0) + D2VP_V2VP(Pi02, Pip, Km, Pi01, 0);
300// PDF[42] = D2VP_V2VP(Pip,Pi01,Km,Pi02,10) + D2VP_V2VP(Pip,Pi02,Km,Pi01,10);
301// PDF[43] = D2VP_V2VP(Pi01,Pi02,Km,Pip,20) + D2VP_V2VP(Pi02,Pi01,Km,Pip,20);
302// PDF[44] = D2VP_V2VP(Pi01,Km,Pip,Pi02,1) + D2VP_V2VP(Pi02,Km,Pip,Pi01,1);
303// PDF[45] = D2VP_V2VP(Km,Pi01,Pip,Pi02,2) + D2VP_V2VP(Km,Pi02,Pip,Pi01,2);
304 //-----------D2VV--------------------------
305 g[0] = 1; g[1] = 1; g[2] = 0;
306 PDF[46] = D2VV(Km, Pi01, Pip, Pi02, g, 0) + D2VV(Km, Pi02, Pip, Pi01, g, 0);
307 g[2] = 1;
308 PDF[47] = D2VV(Km, Pi01, Pip, Pi02, g, 0) + D2VV(Km, Pi02, Pip, Pi01, g, 0);
309 g[2] = 2;
310 PDF[48] = D2VV(Km, Pi01, Pip, Pi02, g, 0) + D2VV(Km, Pi02, Pip, Pi01, g, 0);
311 g[0] = 0; g[1] = 1; g[2] = 0;
312// PDF[49] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
313 g[2] = 1;
314 PDF[50] = D2VV(Km, Pi01, Pip, Pi02, g, 0) + D2VV(Km, Pi02, Pip, Pi01, g, 0);
315 g[2] = 2;
316 PDF[51] = D2VV(Km, Pi01, Pip, Pi02, g, 0) + D2VV(Km, Pi02, Pip, Pi01, g, 0);
317 g[0] = 1; g[1] = 0; g[2] = 0;
318// PDF[52] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
319// g[2] = 1;
320// PDF[53] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
321 g[2] = 2;
322 PDF[54] = D2VV(Km, Pi01, Pip, Pi02, g, 0) + D2VV(Km, Pi02, Pip, Pi01, g, 0);
323 g[0] = 1; g[1] = 0; g[2] = 0;
324// PDF[55] = D2VV(Km,Pip,Pi01,Pi01,g,1) + D2VV(Km,Pip,Pi02,Pi01,g,1);
325// g[2] = 1;
326// PDF[56] = D2VV(Km,Pip,Pi01,Pi01,g,1) + D2VV(Km,Pip,Pi02,Pi01,g,1);
327// g[2] = 2;
328// PDF[57] = D2VV(Km,Pip,Pi01,Pi01,g,1) + D2VV(Km,Pip,Pi02,Pi01,g,1);
329 g[0] = 0; g[1] = 0; g[2] = 0;
330// PDF[58] = D2VV(Km,Pip,Pi01,Pi01,g,1) + D2VV(Km,Pip,Pi02,Pi01,g,1);
331// g[2] = 1;
332// PDF[59] = D2VV(Km,Pip,Pi01,Pi01,g,1) + D2VV(Km,Pip,Pi02,Pi01,g,1);
333// g[2] = 2;
334// PDF[60] = D2VV(Km,Pip,Pi01,Pi01,g,1) + D2VV(Km,Pip,Pi02,Pi01,g,1);
335 g[0] = 0; g[1] = 0; g[2] = 0;
336 PDF[61] = D2VV(Km, Pi01, Pip, Pi02, g, 0) + D2VV(Km, Pi02, Pip, Pi01, g, 0);
337 g[2] = 1;
338// PDF[62] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
339 g[2] = 2;
340// PDF[63] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
341 //----------D2TS--------------------
342// PDF[64] = D2TS(Km,Pip,Pi01,Pi02,0) + D2TS(Km,Pip,Pi02,Pi01,0);
343// PDF[65] = D2TS(Km,Pi01,Pip,Pi02,10) + D2TS(Km,Pi02,Pip,Pi01,10);
344 PDF[66] = D2TS(Pi02, Pi01, Km, Pip, 1) + D2TS(Pi01, Pi02, Km, Pip, 1);
345 PDF[67] = D2TS(Pip, Pi01, Km, Pi02, 11) + D2TS(Pip, Pi02, Km, Pi01, 11);
346
347//------------------------------------------
348 EvtComplex cof;
349 EvtComplex pdf, module;
350 pdf = EvtComplex(0, 0);
351 for (int i = 0; i < 68; i++) {
352 if (mod[i] == 0) continue;
353 cof = EvtComplex(rho[i] * cos(phi[i]), rho[i] * sin(phi[i]));
354 pdf = pdf + cof * PDF[i];
355 }
356 module = conj(pdf) * pdf;
357 double value;
358 value = real(module);
359 return (value <= 0) ? 1e-20 : value;
360 }
361 EvtComplex EvtD0ToKpipi0pi0::KPiSFormfactor(double sa, double sb, double sc, double r)
362 {
363 (void)r;
364 double m1430 = 1.463;
365 double sa0 = m1430 * m1430;
366 double w1430 = 0.233;
367 double q0 = (sa0 + sb - sc) * (sa0 + sb - sc) / (4 * sa0) - sb;
368 if (q0 < 0) q0 = 1e-16;
369 double qs = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
370 double q = sqrt(qs);
371 double Width = w1430 * q * m1430 / sqrt(sa * q0) * r / r;
372 double temp_R = atan(m1430 * Width / (sa0 - sa));
373 if (temp_R < 0) temp_R += math_pi;
374 double deltaR = -5.31 + temp_R;
375 double temp_F = atan(2 * 1.07 * q / (2 + (-1.8) * 1.07 * qs));
376 if (temp_F < 0) temp_F += math_pi;
377 double deltaF = 2.33 + temp_F;
378 EvtComplex cR(cos(deltaR), sin(deltaR));
379 EvtComplex cF(cos(deltaF), sin(deltaF));
380 EvtComplex amp = 0.8 * sin(deltaF) * cF + sin(deltaR) * cR * cF * cF;
381 return amp;
382 }
383 EvtComplex EvtD0ToKpipi0pi0::D2VV(const double* P1, const double* P2, const double* P3, const double* P4, int* g, const int flag)
384 {
385 double t1V1[4], t1V2[4], t1D[4], t2D[4][4];
386 double temp_PDF = 0;
387 EvtComplex amp_PDF(0, 0);
388 EvtComplex pro[2];
389
390 double sa[3], sb[3], sc[3], B[3];
391 double pV1[4], pV2[4], pD[4];
392 for (int i = 0; i != 4; i++) {
393 pV1[i] = P1[i] + P2[i];
394 pV2[i] = P3[i] + P4[i];
395 pD[i] = pV1[i] + pV2[i];
396 }
397 sa[0] = LorentzDotProduct(pV1, pV1);
398 sb[0] = LorentzDotProduct(P1, P1);
399 sc[0] = LorentzDotProduct(P2, P2);
400 sa[1] = LorentzDotProduct(pV2, pV2);
401 sb[1] = LorentzDotProduct(P3, P3);
402 sc[1] = LorentzDotProduct(P4, P4);
403 sa[2] = LorentzDotProduct(pD, pD);
404 sb[2] = sa[0];
405 sc[2] = sa[1];
406 if (g[0] == 1) {
407 if (flag == 0)pro[0] = propagatorRBW(mass[2], width[2], sa[0], sb[0], sc[0], rRes, 1);
408 if (flag == 1)pro[0] = propagatorRBW(mass[3], width[3], sa[0], sb[0], sc[0], rRes, 1);
409 }
410 if (g[1] == 1) {
411 if (flag == 0) pro[1] = propagatorGS(mass[4], width[4], sa[1], sb[1], sc[1], rRes, 1); //rho
412 if (flag == 1) pro[1] = 1;
413 }
414 if (g[0] == 0) pro[0] = 1;
415 if (g[1] == 0) pro[1] = 1;
416 B[0] = BWBarrierFactor(1, sa[0], sb[0], sc[0], rRes);
417 B[1] = BWBarrierFactor(1, sa[1], sb[1], sc[1], rRes);
418 covariantTensor1(P1, P2, t1V1);
419 covariantTensor1(P3, P4, t1V2);
420 if (g[2] == 0) {
421 for (int i = 0; i != 4; i++) {
422 temp_PDF += (G[i][i]) * t1V1[i] * t1V2[i];
423 }
424 B[2] = 1;
425 }
426 if (g[2] == 1) {
427 covariantTensor1(pV1, pV2, t1D);
428 for (int i = 0; i != 4; i++) {
429 for (int j = 0; j != 4; j++) {
430 for (int k = 0; k != 4; k++) {
431 for (int l = 0; l != 4; l++) {
432 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]);
433 }
434 }
435 }
436 }
437 B[2] = BWBarrierFactor(1, sa[2], sb[2], sc[2], rD);
438 }
439 if (g[2] == 2) {
440 covariantTensor2(pV1, pV2, t2D);
441 for (int i = 0; i != 4; i++) {
442 for (int j = 0; j != 4; j++) {
443 temp_PDF += t2D[i][j] * t1V1[i] * t1V2[j] * (G[i][i]) * (G[j][j]);
444 }
445 }
446 B[2] = BWBarrierFactor(2, sa[2], sb[2], sc[2], rD);
447 }
448 amp_PDF = temp_PDF * B[0] * B[1] * B[2] * pro[0] * pro[1];
449 return amp_PDF;
450 }
451
452 EvtComplex EvtD0ToKpipi0pi0::D2AP_A2VP(const double* P1, const double* P2, const double* P3, const double* P4, int* g,
453 const int flag)
454 {
455 double temp_PDF = 0;
456 EvtComplex amp_PDF(0, 0);
457 EvtComplex pro[2];
458 double t1V[4], t1D[4], t2A[4][4];
459 double sa[3], sb[3], sc[3], B[3] = {0, 0, 0};
460 double pV[4], pA[4], pD[4];
461 for (int i = 0; i != 4; i++) {
462 pV[i] = P3[i] + P4[i];
463 pA[i] = pV[i] + P2[i];
464 pD[i] = pA[i] + P1[i];
465 }
466 sa[0] = LorentzDotProduct(pV, pV);
467 sb[0] = LorentzDotProduct(P3, P3);
468 sc[0] = LorentzDotProduct(P4, P4);
469 sa[1] = LorentzDotProduct(pA, pA);
470 sb[1] = sa[0];
471 sc[1] = LorentzDotProduct(P2, P2);
472 sa[2] = LorentzDotProduct(pD, pD);
473 sb[2] = sa[1];
474 sc[2] = LorentzDotProduct(P1, P1);
475 if (g[0] == 1) {
476 if (flag == 0 || flag == 3) pro[0] = propagatorGS(mass[4], width[4], sa[0], sb[0], sc[0], rRes, 1);
477 else if (flag == 1 || flag == 21) pro[0] = propagatorRBW(mass[2], width[2], sa[0], sb[0], sc[0], rRes, 1);
478 else if (flag == 31) pro[0] = propagatorRBW(mass[3], width[3], sa[0], sb[0], sc[0], rRes, 1);
479 } else if (g[0] == 0) pro[0] = 1;
480 if (g[1] == 1) {
481 if (flag == 0) pro[1] = propagatorRBW(mass[0], width[0], sa[1], sb[1], sc[1], rRes, g[2]);
482 if (flag == 1 || flag == 21 || flag == 31 || flag == 3) pro[1] = propagatorRBW(mass[1], width[1], sa[1], sb[1], sc[1], rRes, g[2]);
483 } else if (g[1] == 0) pro[1] = 1;
484 B[0] = BWBarrierFactor(1, sa[0], sb[0], sc[0], rRes);
485 B[2] = BWBarrierFactor(1, sa[2], sb[2], sc[2], rD);
486 covariantTensor1(P3, P4, t1V);
487 covariantTensor1(pA, P1, t1D);
488 if (g[2] == 0) {
489 for (int i = 0; i != 4; i++) {
490 for (int j = 0; j != 4; j++) {
491 temp_PDF += t1D[i] * (G[i][j] - pA[i] * pA[j] / sa[1]) * t1V[j] * (G[i][i]) * (G[j][j]);
492 }
493 }
494 B[1] = 1;
495 } else if (g[2] == 2) {
496 covariantTensor2(pV, P2, t2A);
497 for (int i = 0; i != 4; i++) {
498 for (int j = 0; j != 4; j++) {
499 temp_PDF += t1D[i] * t2A[i][j] * t1V[j] * (G[i][i]) * (G[j][j]);
500 }
501 }
502 B[1] = BWBarrierFactor(2, sa[1], sb[1], sc[1], rRes);
503 }
504 amp_PDF = temp_PDF * B[0] * B[1] * B[2] * pro[0] * pro[1];
505 return amp_PDF;
506 }
507
508 EvtComplex EvtD0ToKpipi0pi0::D2AP_A2SP(const double* P1, const double* P2, const double* P3, const double* P4, const int flag)
509 {
510 //flag = 0, S = rho; flag = 1, S = K*
511 double temp_PDF = 0;
512 EvtComplex amp_PDF(0, 0);
513 EvtComplex pro;
514 double sa[3], sb[3], sc[3], B[3];
515 double t1D[4], t1A[4];
516 double pS[4], pA[4], pD[4];
517 for (int i = 0; i != 4; i++) {
518 pS[i] = P3[i] + P4[i];
519 pA[i] = pS[i] + P2[i];
520 pD[i] = pA[i] + P1[i];
521 }
522 sa[0] = LorentzDotProduct(pS, pS);
523 sb[0] = LorentzDotProduct(P3, P3);
524 sc[0] = LorentzDotProduct(P4, P4);
525 sa[1] = LorentzDotProduct(pA, pA);
526 sb[1] = sa[0];
527 sc[1] = LorentzDotProduct(P2, P2);
528 sa[2] = LorentzDotProduct(pD, pD);
529 sb[2] = sa[1];
530 sc[2] = LorentzDotProduct(P1, P1);
531 B[1] = BWBarrierFactor(1, sa[1], sb[1], sc[1], rRes);
532 B[2] = BWBarrierFactor(1, sa[2], sb[2], sc[2], rD);
533 covariantTensor1(pA, P1, t1D);
534 covariantTensor1(pS, P2, t1A);
535 for (int i = 0; i != 4; i++) {
536 temp_PDF += t1D[i] * t1A[i] * (G[i][i]);
537 }
538 amp_PDF = temp_PDF * B[1] * B[2];
539 if (flag == 1 || flag == 11 || flag == 21) amp_PDF = amp_PDF * KPiSFormfactor(sa[0], sb[0], sc[0], rRes);
540 return amp_PDF;
541 }
542
543 EvtComplex EvtD0ToKpipi0pi0::D2PP_P2VP(const double* P1, const double* P2, const double* P3, const double* P4, const int flag)
544 {
545 //modeidx = 0 :(K*0 pi0)pi0
546 //modeidx = 10:(K*- pi+)pi0
547 //modeidx = 20:(K*- pi0)pi+
548 //modeidx = 1 :(K- rho+)pi0
549 //modeidx = 11:K-(rho+ pi0)
550 double temp_PDF = 0;
551 EvtComplex amp(0, 0);
552 EvtComplex prop;
553 double sa[3], sb[3], sc[3], B[3];
554 double t1V[4];
555 double pV[4], pP[4], pD[4];
556 for (int i = 0; i != 4; i++) {
557 pV[i] = P3[i] + P4[i];
558 pP[i] = pV[i] + P2[i];
559 pD[i] = pP[i] + P1[i];
560 }
561 sa[0] = LorentzDotProduct(pV, pV);
562 sb[0] = LorentzDotProduct(P3, P3);
563 sc[0] = LorentzDotProduct(P4, P4);
564 sa[1] = LorentzDotProduct(pP, pP);
565 sb[1] = sa[0];
566 sc[1] = LorentzDotProduct(P2, P2);
567 sa[2] = LorentzDotProduct(pD, pD);
568 sb[2] = sa[1];
569 sc[2] = LorentzDotProduct(P1, P1);
570 B[0] = BWBarrierFactor(1, sa[0], sb[0], sc[0], rRes);
571 B[1] = BWBarrierFactor(1, sa[1], sb[1], sc[1], rRes);
572 if (flag == 0) prop = propagatorRBW(mass[3], width[3], sa[0], sb[0], sc[0], rRes, 1);
573 else if (flag == 10 || 20) prop = propagatorRBW(mass[2], width[2], sa[0], sb[0], sc[0], rRes, 1);
574 else if (flag == 1 || 11) prop = propagatorGS(mass[4], width[4], sa[0], sb[0], sc[0], rRes, 1);
575 covariantTensor1(P3, P4, t1V);
576 for (int i = 0; i != 4; i++) {
577 temp_PDF += P2[i] * t1V[i] * (G[i][i]);
578 }
579 amp = temp_PDF * B[0] * B[1] * prop;
580 return amp;
581 }
582
583 EvtComplex EvtD0ToKpipi0pi0::D2VP_V2VP(const double* P1, const double* P2, const double* P3, const double* P4, const int flag)
584 {
585 double temp_PDF = 0;
586 EvtComplex amp_PDF(0, 0);
587 EvtComplex pro;
588 double sa[3], sb[3], sc[3], B[3];
589 double pV1[4], pV2[4], qV1[4], qV2[4], pD[4];
590 for (int i = 0; i != 4; i++) {
591 pV2[i] = P3[i] + P4[i];
592 qV2[i] = P3[i] - P4[i];
593 pV1[i] = pV2[i] + P2[i];
594 qV1[i] = pV2[i] - P2[i];
595 pD[i] = pV1[i] + P1[i];
596 }
597 for (int i = 0; i != 4; i++) {
598 for (int j = 0; j != 4; j++) {
599 for (int k = 0; k != 4; k++) {
600 for (int l = 0; l != 4; l++) {
601 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]);
602 }
603 }
604 }
605 }
606 sa[0] = LorentzDotProduct(pV2, pV2);
607 sb[0] = LorentzDotProduct(P3, P3);
608 sc[0] = LorentzDotProduct(P4, P4);
609 sa[1] = LorentzDotProduct(pV1, pV1);
610 sb[1] = sa[0];
611 sc[1] = LorentzDotProduct(P2, P2);
612 sa[2] = LorentzDotProduct(pD, pD);
613 sb[2] = sa[1];
614 sc[2] = LorentzDotProduct(P1, P1);
615 if (flag == 0 || flag == 10) pro = propagatorRBW(mass[2], width[2], sa[0], sb[0], sc[0], rRes, 1);
616 else if (flag == 20) pro = propagatorRBW(mass[3], width[3], sa[0], sb[0], sc[0], rRes, 1);
617 else if (flag == 1 || flag == 2) pro = propagatorGS(mass[4], width[4], sa[0], sb[0], sc[0], rRes, 1);
618 B[0] = BWBarrierFactor(1, sa[0], sb[0], sc[0], rRes);
619 B[1] = BWBarrierFactor(1, sa[1], sb[1], sc[1], rRes);
620 B[2] = BWBarrierFactor(1, sa[2], sb[2], sc[2], rD);
621 amp_PDF = temp_PDF * B[0] * B[1] * B[2] * pro;
622 return amp_PDF;
623 }
624
625 EvtComplex EvtD0ToKpipi0pi0::D2VS(const double* P1, const double* P2, const double* P3, const double* P4, int g, const int flag)
626 {
627 double temp_PDF = 0;
628 EvtComplex amp_PDF(0, 0);
629 EvtComplex pro;
630 double sa[3], sb[3], sc[3], B[3];
631 double t1D[4], t1V[4];
632 double pS[4], pV[4], pD[4];
633 for (int i = 0; i != 4; i++) {
634 pS[i] = P3[i] + P4[i];
635 pV[i] = P1[i] + P2[i];
636 pD[i] = pS[i] + pV[i];
637 }
638 sa[0] = LorentzDotProduct(pS, pS);
639 sb[0] = LorentzDotProduct(P3, P3);
640 sc[0] = LorentzDotProduct(P4, P4);
641 sa[1] = LorentzDotProduct(pV, pV);
642 sb[1] = LorentzDotProduct(P1, P1);
643 sc[1] = LorentzDotProduct(P2, P2);
644 sa[2] = LorentzDotProduct(pD, pD);
645 sb[2] = sa[0];
646 sc[2] = sa[1];
647 if (g == 1) {
648 if (flag == 0) pro = propagatorGS(mass[4], width[4], sa[1], sb[1], sc[1], rRes, 1);
649 else if (flag == 1) pro = propagatorRBW(mass[2], width[2], sa[1], sb[1], sc[1], rRes, 1);
650 else if (flag == 11) pro = propagatorRBW(mass[3], width[3], sa[1], sb[1], sc[1], rRes, 1);
651 else if (flag == 10) pro = 1;
652 } else if (g == 0) pro = 1;
653 B[1] = BWBarrierFactor(1, sa[1], sb[1], sc[1], rRes);
654 B[2] = BWBarrierFactor(1, sa[2], sb[2], sc[2], rD);
655 covariantTensor1(P1, P2, t1V);
656 covariantTensor1(pS, pV, t1D);
657 for (int i = 0; i != 4; i++) {
658 temp_PDF += G[i][i] * t1D[i] * t1V[i];
659 }
660 amp_PDF = temp_PDF * B[1] * B[2] * pro;
661 if (flag == 0 || flag == 10) amp_PDF *= KPiSFormfactor(sa[0], sb[0], sc[0], rRes);
662 // if(modeidx == 1 || modeidx == 11) amp_PDF *= -1.0; /// why ???????
663 return amp_PDF;
664 }
665
666 EvtComplex EvtD0ToKpipi0pi0::D2TS(const double* P1, const double* P2, const double* P3, const double* P4, const int flag)
667 {
668 // flag == 0 KPiT. 1 PiPiT
669 double temp_PDF = 0;
670 EvtComplex amp_PDF(0, 0);
671 double sa[3], sb[3], sc[3], B[3];
672 double t2D[4][4], t2T[4][4];
673 double pS[4], pT[4], pD[4];
674 for (int i = 0; i != 4; i++) {
675 pS[i] = P3[i] + P4[i];
676 pT[i] = P1[i] + P2[i];
677 pD[i] = pT[i] + pS[i];
678 }
679 sa[0] = LorentzDotProduct(pT, pT);
680 sb[0] = LorentzDotProduct(P1, P1);
681 sc[0] = LorentzDotProduct(P2, P2);
682 sa[1] = LorentzDotProduct(pS, pS);
683 sb[1] = LorentzDotProduct(P3, P3);
684 sc[1] = LorentzDotProduct(P4, P4);
685 sa[2] = LorentzDotProduct(pD, pD);
686 sb[2] = sa[0];
687 sc[2] = sa[1];
688 B[0] = BWBarrierFactor(2, sa[0], sb[0], sc[0], rRes);
689 B[2] = BWBarrierFactor(2, sa[2], sb[2], sc[2], rD);
690 covariantTensor2(P1, P2, t2T);
691 covariantTensor2(pT, pS, t2D);
692 for (int i = 0; i != 4; i++) {
693 for (int j = 0; j != 4; j++) {
694 temp_PDF += t2D[i][j] * t2T[j][i] * (G[i][i]) * (G[j][j]);
695 }
696 }
697 amp_PDF = temp_PDF * B[0] * B[2];
698 if (flag == 1 || flag == 11) { amp_PDF = amp_PDF * KPiSFormfactor(sa[1], sb[1], sc[1], rRes);}
699 return amp_PDF;
700 }
701
702 EvtComplex EvtD0ToKpipi0pi0::PHSP(double* P1, double* P2)
703 {
704 EvtComplex amp_PDF(0, 0);
705 double sa, sb, sc;
706 double KPi[4];
707 for (int i = 0; i != 4; i++) {
708 KPi[i] = P1[i] + P2[i];
709 }
710 sa = LorentzDotProduct(KPi, KPi);
711 sb = LorentzDotProduct(P1, P1);
712 sc = LorentzDotProduct(P2, P2);
713 amp_PDF = KPiSFormfactor(sa, sb, sc, rRes);
714 return amp_PDF;
715 }
716
717 double EvtD0ToKpipi0pi0::energyDependentWidth(const double Mass, double sa, double sb, double sc, double r, int l)const
718 {
719 double widm(0.), q(0.), q0(0.);
720 double sa0 = Mass * Mass;
721 double m = sqrt(sa);
722 q = Qabcs(sa, sb, sc);
723 q0 = Qabcs(sa0, sb, sc);
724 double z, z0;
725 z = q * r * r;
726 z0 = q0 * r * r;
727 double t = q / q0;
728 double F(0.);
729 if (l == 0) F = 1;
730 if (l == 1) F = sqrt((1 + z0) / (1 + z));
731 if (l == 2) F = sqrt((9 + 3 * z0 + z0 * z0) / (9 + 3 * z + z * z));
732 widm = pow(t, l + 0.5) * Mass / m * F * F;
733 return widm;
734 }
735 double EvtD0ToKpipi0pi0::h(double m, double q)const
736 {
737 double h(0.);
738 h = 2 / pi * q / m * log((m + 2 * q) / (2 * mpi));
739 return h;
740 }
741 double EvtD0ToKpipi0pi0::dh(double Mass, double q0)const
742 {
743 double dh = h(Mass, q0) * (1.0 / (8 * q0 * q0) - 1.0 / (2 * Mass * Mass)) + 1.0 / (2 * pi * Mass * Mass);
744 return dh;
745 }
746 double EvtD0ToKpipi0pi0::f(double Mass, double sx, double q0, double q)const
747 {
748 double m = sqrt(sx);
749 double f = Mass * Mass / (pow(q0, 3)) * (q * q * (h(m, q) - h(Mass, q0)) + (Mass * Mass - sx) * q0 * q0 * dh(Mass,
750 q0));
751 return f;
752 }
753 double EvtD0ToKpipi0pi0::d(double Mass, double q0)const
754 {
755 double d = 3.0 / pi * mpi * mpi / (q0 * q0) * log((Mass + 2 * q0) / (2 * mpi)) + Mass / (2 * pi * q0) - (mpi * mpi * Mass) /
756 (pi * pow(q0, 3));
757 return d;
758 }
759 EvtComplex EvtD0ToKpipi0pi0::propagatorRBW(double Mass, double Width, double sa, double sb, double sc, double r,
760 int l)const
761 {
762 EvtComplex ci(0, 1);
763 EvtComplex prop = 1.0 / (Mass * Mass - sa - ci * Mass * Width * energyDependentWidth(Mass, sa, sb, sc, r,
764 l));
765 return prop;
766 }
767 EvtComplex EvtD0ToKpipi0pi0::propagatorGS(double Mass, double Width, double sa, double sb, double sc, double r,
768 int l)const
769 {
770 EvtComplex ci(0, 1);
771 double q = Qabcs(sa, sb, sc);
772 double sa0 = Mass * Mass;
773 double q0 = Qabcs(sa0, sb, sc);
774 q = sqrt(q);
775 q0 = sqrt(q0);
776 EvtComplex prop = (1 + d(Mass, q0) * Width / Mass) / (Mass * Mass - sa + Width * f(Mass, sa, q0,
777 q) - ci * Mass * Width * energyDependentWidth(Mass,
778 sa, sb, sc, r, l));
779 return prop;
780 }
781 double EvtD0ToKpipi0pi0::LorentzDotProduct(const double* a1, const double* a2)const
782 {
783 double dot = 0;
784 for (int i = 0; i != 4; i++) {
785 dot += a1[i] * a2[i] * G[i][i];
786 }
787 return dot;
788 }
789 double EvtD0ToKpipi0pi0::Qabcs(double sa, double sb, double sc)const
790 {
791 double Qabcs = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
792 if (Qabcs < 0) Qabcs = 1e-16;
793 return Qabcs;
794 }
795 double EvtD0ToKpipi0pi0::BWBarrierFactor(double l, double sa, double sb, double sc, double r)const
796 {
797 double q = Qabcs(sa, sb, sc);
798 q = sqrt(q);
799 double z = q * r;
800 z = z * z;
801 double F = 1;
802 if (l > 2) F = 0;
803 if (l == 0)
804 F = 1;
805 if (l == 1) {
806 F = sqrt((2 * z) / (1 + z));
807 }
808 if (l == 2) {
809 F = sqrt((13 * z * z) / (9 + 3 * z + z * z));
810 }
811 return F;
812 }
813 void EvtD0ToKpipi0pi0::covariantTensor1(const double* daug1, const double* daug2, double* t1) const
814 {
815 double p, pq;
816 double pa[4], qa[4];
817 for (int i = 0; i != 4; i++) {
818 pa[i] = daug1[i] + daug2[i];
819 qa[i] = daug1[i] - daug2[i];
820 }
821 p = LorentzDotProduct(pa, pa);
822 pq = LorentzDotProduct(pa, qa);
823 for (int i = 0; i != 4; i++) {
824 t1[i] = qa[i] - pq / p * pa[i];
825 }
826 }
827 void EvtD0ToKpipi0pi0::covariantTensor2(const double* daug1, const double* daug2, double (*t2)[4]) const
828 {
829 double p, r;
830 double pa[4], t1[4];
831 covariantTensor1(daug1, daug2, t1);
832 r = LorentzDotProduct(t1, t1);
833 for (int i = 0; i != 4; i++) {
834 pa[i] = daug1[i] + daug2[i];
835 }
836 p = LorentzDotProduct(pa, pa);
837 for (int i = 0; i != 4; i++) {
838 for (int j = 0; j != 4; j++) {
839 t2[i][j] = t1[i] * t1[j] - 1.0 / 3 * r * (G[i][j] - pa[i] * pa[j] / p);
840 }
841 }
842 }
843
845} // Belle 2 namespace
Belle2::EvtD0ToKpipi0pi0::mD
double mD
Fixed parameters.
Definition EvtD0ToKpipi0pi0.h:91
Belle2::EvtD0ToKpipi0pi0::d
double d(double mass, const double q0) const
d function in Gounaris-Sakurai lineshape
Definition EvtD0ToKpipi0pi0.cc:753
Belle2::EvtD0ToKpipi0pi0::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 EvtD0ToKpipi0pi0.cc:759
Belle2::EvtD0ToKpipi0pi0::dh
double dh(double mass, const double q0) const
derivative h function in Gounaris-Sakurai lineshape
Definition EvtD0ToKpipi0pi0.cc:741
Belle2::EvtD0ToKpipi0pi0::BWBarrierFactor
double BWBarrierFactor(const double l, const double sa, const double sb, const double sc, const double r) const
Blatt-Weisskopf barrier factors.
Definition EvtD0ToKpipi0pi0.cc:795
Belle2::EvtD0ToKpipi0pi0::f
double f(double mass, const double sx, const double q0, const double q) const
f function in Gounaris-Sakurai lineshape
Definition EvtD0ToKpipi0pi0.cc:746
Belle2::EvtD0ToKpipi0pi0::covariantTensor1
void covariantTensor1(const double *daug1, const double *daug2, double *t1) const
Covariant Spin-1 Projector.
Definition EvtD0ToKpipi0pi0.cc:813
Belle2::EvtD0ToKpipi0pi0::LorentzDotProduct
double LorentzDotProduct(const double *a1, const double *a2) const
Four-vector dot product.
Definition EvtD0ToKpipi0pi0.cc:781
Belle2::EvtD0ToKpipi0pi0::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 EvtD0ToKpipi0pi0.cc:767
Belle2::EvtD0ToKpipi0pi0::D2VV
EvtComplex D2VV(const double *P1, const double *P2, const double *P3, const double *P4, int *g, const int flag)
Amplitude modes.
Definition EvtD0ToKpipi0pi0.cc:383
Belle2::EvtD0ToKpipi0pi0::PDF
double PDF(double *Km, double *Pip, double *Pi01, double *Pi02)
Probability distribution function of the decay.
Definition EvtD0ToKpipi0pi0.cc:223
Belle2::EvtD0ToKpipi0pi0::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 EvtD0ToKpipi0pi0.cc:789
Belle2::EvtD0ToKpipi0pi0::KPiSFormfactor
EvtComplex KPiSFormfactor(const double sa, const double sb, const double sc, const double r)
K pi S-wave form factor.
Definition EvtD0ToKpipi0pi0.cc:361
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::EvtD0ToKpipi0pi0::covariantTensor2
void covariantTensor2(const double *daug1, const double *daug2, double(*t2)[4]) const
Covariant Spin-2 Projector.
Definition EvtD0ToKpipi0pi0.cc:827
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
Belle2::EvtD0ToKpipi0pi0::h
double h(const double m, const double q) const
h function in Gounaris-Sakurai lineshape
Definition EvtD0ToKpipi0pi0.cc:735
Belle2::dot
T dot(GeneralVector< T > a, GeneralVector< T > b)
dot product of two general vectors
Definition beamHelpers.h:163
Belle2::EvtD0ToKpipi0pi0::energyDependentWidth
double energyDependentWidth(const double mass, const double sa, const double sb, const double sc, const double r, const int l) const
Energy dependent width.
Definition EvtD0ToKpipi0pi0.cc:717
Belle2
Abstract base class for different kinds of events.
Definition MillepedeAlgorithm.h:17