Bug Summary

File:generators/evtgen/models/besiii/src/EvtDsToKpipi.cc
Warning:line 178, column 20
The right operand of '*' is a garbage value

Annotated Source Code

Press '?' to see keyboard shortcuts

clang -cc1 -cc1 -triple x86_64-unknown-linux-gnu -O3 -analyze -disable-free -clear-ast-before-backend -disable-llvm-verifier -discard-value-names -main-file-name EvtDsToKpipi.cc -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=cplusplus -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -setup-static-analyzer -mrelocation-model pic -pic-level 2 -fhalf-no-semantic-interposition -mframe-pointer=none -fmath-errno -ffp-contract=on -fno-rounding-math -mconstructor-aliases -funwind-tables=2 -target-cpu x86-64 -tune-cpu generic -debugger-tuning=gdb -fdebug-compilation-dir=/data/b2soft/buildbot/development/build -fcoverage-compilation-dir=/data/b2soft/buildbot/development/build -resource-dir /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/lib/clang/21 -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/c++ -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/c++/x86_64-redhat-linux -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/c++/backward -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/include -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/python3.12 -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/include/CLHEP -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/Geant4 -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/include/root -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/include/belle_legacy -I include/ -D _PACKAGE_="generators" -D G4UI_USE_TCSH -D RaveDllExport= -D HAS_SQLITE -D HAS_CALLGRIND -I include -I /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/libxml2 -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/bin/../lib64/gcc/x86_64-redhat-linux/15.2.0/../../../../include/c++ -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/bin/../lib64/gcc/x86_64-redhat-linux/15.2.0/../../../../include/c++/x86_64-redhat-linux -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/bin/../lib64/gcc/x86_64-redhat-linux/15.2.0/../../../../include/c++/backward -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/lib/clang/21/include -internal-isystem /usr/local/include -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/bin/../lib64/gcc/x86_64-redhat-linux/15.2.0/../../../../x86_64-redhat-linux/include -internal-externc-isystem /include -internal-externc-isystem /usr/include -Wno-missing-braces -Wno-unused-command-line-argument -std=c++20 -fdeprecated-macro -ferror-limit 19 -fgnuc-version=4.2.1 -fno-implicit-modules -fskip-odr-check-in-gmf -fcxx-exceptions -fexceptions -vectorize-loops -vectorize-slp -analyzer-output=html -faddrsig -D__GCC_HAVE_DWARF2_CFI_ASM=1 -o /scan_build/2026-05-31-004316-385593-1 -x c++ generators/evtgen/models/besiii/src/EvtDsToKpipi.cc
1// Model: EvtDsToKpipi
2// This file is an amplitude model for Ds+ -> K- pi+ pi+.
3// The model is from the BESIII Collaboration in JHEP08 (2022) 196. DOI: https://doi.org/10.1007/JHEP08(2022)196
4//
5// Permission to use these files in basf2 was generously granted by the BESIII Collaboration.
6//
7// Please cite the original reference for any public/published results where this model was used.
8
9#include <iostream>
10#include <cmath>
11#include <string>
12#include <EvtGenBase/EvtCPUtil.hh>
13#include <EvtGenBase/EvtTensor4C.hh>
14#include <EvtGenBase/EvtPatches.hh>
15#include <fstream>
16#include <stdlib.h>
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>
23#include <string>
24#include <EvtGenBase/EvtConst.hh>
25#include <EvtGenBase/EvtDecayTable.hh>
26
27#include <generators/evtgen/EvtGenModelRegister.h>
28#include <generators/evtgen/models/besiii/EvtDsToKpipi.h>
29
30namespace Belle2 {
31
32 /** register the model in EvtGen */
33 B2_EVTGEN_REGISTER_MODEL(EvtDsToKpipi)namespace { Belle2::EvtGenModelRegister::Factory<EvtDsToKpipi
> EvtGenModelFactory_EvtDsToKpipi; };
34
35 EvtDsToKpipi::~EvtDsToKpipi() {}
36
37 std::string EvtDsToKpipi::getName()
38 {
39 return "DsToKpipi";
40 }
41
42 EvtDecayBase* EvtDsToKpipi::clone()
43 {
44 return new EvtDsToKpipi;
45 }
46
47 void EvtDsToKpipi::init()
48 {
49 checkNArg(0);
50 checkNDaug(3);
51 checkSpinParent(EvtSpinType::SCALAR);
52 checkSpinDaughter(0, EvtSpinType::SCALAR);
53 checkSpinDaughter(1, EvtSpinType::SCALAR);
54 checkSpinDaughter(2, EvtSpinType::SCALAR);
55
56 phi[0] = 0;
57 rho[0] = 1;
58 phi[1] = 0;
59 rho[1] = 0;
60 phi[2] = 0;
61 rho[2] = 0;
62 phi[3] = 0;
63 rho[3] = 0;
64 phi[4] = 0;
65 rho[4] = 0;
66 phi[5] = 0;
67 rho[5] = 0;
68 phi[6] = 0;
69 rho[6] = 0;
70 phi[7] = 0;
71 rho[7] = 0;
72
73 phi[1] = -3.47995752;
74 phi[2] = -1.249864467;
75 phi[3] = -0.3067504308;
76 phi[4] = 0.9229242379;
77 phi[5] = -3.278567926;
78 phi[6] = -0.6346408622;
79 phi[7] = 1.762377475;
80
81 rho[1] = 2.463898984;
82 rho[2] = 0.7361782665;
83 rho[3] = 1.90216812;
84 rho[4] = 2.140687169;
85 rho[5] = 0.914684056;
86 rho[6] = 1.116206881;
87 rho[7] = 1.483440555;
88
89 modetype[0] = 23;
90 modetype[1] = 23;
91 modetype[2] = 23;
92 modetype[3] = 23;
93 modetype[4] = 23;
94 modetype[5] = 13;
95 modetype[6] = 13;
96 modetype[7] = 13;
97
98 width[0] = 0.1478;
99 width[1] = 0.400;
100 width[2] = 0.100;
101 width[3] = 0.265;
102 width[4] = 0.270;
103 width[5] = 0.0473;
104 width[6] = 0.232;
105 width[7] = 0.270;
106
107 mass[0] = 0.77526;
108 mass[1] = 1.465;
109 mass[2] = 0.965;
110 mass[3] = 1.35;
111 mass[4] = 1.425;
112 mass[5] = 0.89555;
113 mass[6] = 1.414;
114 mass[7] = 1.432787726;
115
116 mDsM = 1.9683;
117 mD = 1.86486;
118 mDs = 1.9683;
119 rD = 25.0;
120 metap = 0.95778;
121 mkstr = 0.89594;
122 mk0 = 0.497614;
123 mass_Kaon = 0.49368;
124 mass_Pion = 0.13957;
125 mass_Pi0 = 0.1349766;
126 math_pi = 3.1415926;
127 mKa2 = 0.24371994;
128 mPi2 = 0.01947978;
129 mPi = 0.13957;
130 mKa = 0.49368;
131 ma0 = 0.99;
132 Ga0 = 0.0756;
133 meta = 0.547862;
134
135 GS1 = 0.636619783;
136 GS2 = 0.01860182466;
137 GS3 = 0.1591549458;
138 GS4 = 0.00620060822;
139
140 int GG[4][4] = { {1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1} };
141 for (int i = 0; i < 4; i++) {
142 for (int j = 0; j < 4; j++) {
143 G[i][j] = GG[i][j];
144 }
145 }
146 }
147
148 void EvtDsToKpipi::initProbMax()
149 {
150 setProbMax(825.0);
151 }
152
153 void EvtDsToKpipi::decay(EvtParticle* p)
154 {
155 p->initializePhaseSpace(getNDaug(), getDaugs());
156 EvtVector4R D1 = p->getDaug(0)->getP4();
157 EvtVector4R D2 = p->getDaug(1)->getP4();
158 EvtVector4R D3 = p->getDaug(2)->getP4();
159
160 double P1[4], P2[4], P3[4];
161 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
162 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
163 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
164
165 double value;
166 int g0[8] = {1, 1, 4, 2, 500, 1, 1, 2};
167 int spin_param[8] = {1, 1, 0, 0, 0, 1, 1, 0};
168 int nstates = 8;
169 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin_param, modetype, nstates, value);
170
171 setProb(value);
172
173 return ;
174 }
175
176 void EvtDsToKpipi::Com_Multi(double a1[2], double a2[2], double res[2])
177 {
178 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
9
←
The right operand of '*' is a garbage value
179 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
180 }
181 void EvtDsToKpipi::Com_Divide(double a1[2], double a2[2], double res[2])
182 {
183 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
184 res[0] = (a1[0] * a2[0] + a1[1] * a2[1]) / tmp;
185 res[1] = (a1[1] * a2[0] - a1[0] * a2[1]) / tmp;
186 }
187
188 double EvtDsToKpipi::SCADot(double a1[4], double a2[4])
189 {
190 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
191 return _cal;
192 }
193 double EvtDsToKpipi::barrier(int l, double sa, double sb, double sc, double r, double mass_param)
194 {
195 double q = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
196 if (q < 0) q = 1e-16;
197 double z;
198 z = q * r * r;
199 double sa0;
200 sa0 = mass_param * mass_param;
201 double q0 = (sa0 + sb - sc) * (sa0 + sb - sc) / (4 * sa0) - sb;
202 if (q0 < 0) q0 = 1e-16;
203 double z0 = q0 * r * r;
204 double F = 0.0;
205 if (l == 0) F = 1;
206 if (l == 1) F = sqrt((1 + z0) / (1 + z));
207 if (l == 2) F = sqrt((9 + 3 * z0 + z0 * z0) / (9 + 3 * z + z * z));
208 return F;
209 }
210 void EvtDsToKpipi::calt1(double daug1[4], double daug2[4], double t1[4])
211 {
212 double p, pq, tmp;
213 double pa[4], qa[4];
214 for (int i = 0; i < 4; i++) {
215 pa[i] = daug1[i] + daug2[i];
216 qa[i] = daug1[i] - daug2[i];
217 }
218 p = SCADot(pa, pa);
219 pq = SCADot(pa, qa);
220 tmp = pq / p;
221 for (int i = 0; i < 4; i++) {
222 t1[i] = qa[i] - tmp * pa[i];
223 }
224 }
225 void EvtDsToKpipi::calt2(double daug1[4], double daug2[4], double t2[4][4])
226 {
227 double p, r;
228 double pa[4], t1[4];
229 calt1(daug1, daug2, t1);
230 r = SCADot(t1, t1) / 3.0;
231 for (int i = 0; i < 4; i++) {
232 pa[i] = daug1[i] + daug2[i];
233 }
234 p = SCADot(pa, pa);
235 for (int i = 0; i < 4; i++) {
236 for (int j = 0; j < 4; j++) {
237 t2[i][j] = t1[i] * t1[j] - r * (G[i][j] - pa[i] * pa[j] / p);
238 }
239 }
240 }
241
242//-------------------prop--------------------------------------------
243
244 void EvtDsToKpipi::propagatorCBW(double mass_param, double width_param, double sx, double prop[2])
245 {
246 double a[2], b[2];
247 a[0] = 1;
248 a[1] = 0;
249 b[0] = mass_param * mass_param - sx;
250 b[1] = -mass_param * width_param;
251 Com_Divide(a, b, prop);
252 }
253 double EvtDsToKpipi::wid(double mass2, double mass_param, double sa, double sb, double sc, double r2, int l)
254 {
255 double widm = 0.;
256 double m = sqrt(sa);
257 double tmp = sb - sc;
258 double tmp1 = sa + tmp;
259 double q = 0.25 * tmp1 * tmp1 / sa - sb;
260 if (q < 0) q = 1e-16;
261 double tmp2 = mass2 + tmp;
262 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
263 if (q0 < 0) q0 = 1e-16;
264 double z = q * r2;
265 double z0 = q0 * r2;
266 double t = q / q0;
267 if (l == 0) {widm = sqrt(t) * mass_param / m;}
268 else if (l == 1) {widm = t * sqrt(t) * mass_param / m * (1 + z0) / (1 + z);}
269 else if (l == 2) {widm = t * t * sqrt(t) * mass_param / m * (9 + 3 * z0 + z0 * z0) / (9 + 3 * z + z * z);}
270 return widm;
271 }
272
273 double EvtDsToKpipi::widl1(double mass2, double mass_param, double sa, double sb, double sc, double r2)
274 {
275 double widm = 0.;
276 double m = sqrt(sa);
277 double tmp = sb - sc;
278 double tmp1 = sa + tmp;
279 double q = 0.25 * tmp1 * tmp1 / sa - sb;
280 if (q < 0) q = 1e-16;
281 double tmp2 = mass2 + tmp;
282 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
283 if (q0 < 0) q0 = 1e-16;
284 double z = q * r2;
285 double z0 = q0 * r2;
286 double F = (1 + z0) / (1 + z);
287 double t = q / q0;
288 widm = t * sqrt(t) * mass_param / m * F;
289 return widm;
290 }
291 void EvtDsToKpipi::propagatorRBW(double mass_param, double width_param, double sa, double sb, double sc, double r2, int l,
292 double prop[2])
293 {
294 double a[2], b[2];
295 double mass2 = mass_param * mass_param;
296
297 a[0] = 1;
298 a[1] = 0;
299 b[0] = mass2 - sa;
300 b[1] = -mass_param * width_param * wid(mass2, mass_param, sa, sb, sc, r2, l);
301 Com_Divide(a, b, prop);
302 }
303 void EvtDsToKpipi::propagatorFlatte(double mass_param, double width_param __attribute__((unused)), double sa, double prop[2])
304 {
305
306 double q2_Pi, q2_Ka;
307 double rhoPi[2] = {0, 0}, rhoKa[2] = {0, 0};
308
309 q2_Pi = 0.25 * sa - mPi * mPi;
310 q2_Ka = 0.25 * sa - mKa * mKa;
311
312 if (q2_Pi > 0) {
313 rhoPi[0] = 2.0 * sqrt(q2_Pi / sa);
314 rhoPi[1] = 0.0;
315 }
316 if (q2_Pi <= 0) {
317 rhoPi[0] = 0.0;
318 rhoPi[1] = 2.0 * sqrt(-q2_Pi / sa);
319 }
320
321 if (q2_Ka > 0) {
322 rhoKa[0] = 2.0 * sqrt(q2_Ka / sa);
323 rhoKa[1] = 0.0;
324 }
325 if (q2_Ka <= 0) {
326 rhoKa[0] = 0.0;
327 rhoKa[1] = 2.0 * sqrt(-q2_Ka / sa);
328 }
329
330 double a[2], b[2];
331 a[0] = 1;
332 a[1] = 0;
333 b[0] = mass_param * mass_param - sa + 0.165 * rhoPi[1] + 0.69465 * rhoKa[1];
334 b[1] = - (0.165 * rhoPi[0] + 0.69465 * rhoKa[0]);
335 Com_Divide(a, b, prop);
336 }
337 void EvtDsToKpipi::propagatorGS(double mass_param, double width_param, double sa, double sb, double sc, double r2, double prop[2])
338 {
339 double a[2], b[2];
340 double mass2 = mass_param * mass_param;
341 double tmp = sb - sc;
342 double tmp1 = sa + tmp;
343 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
344 if (q2 < 0) q2 = 1e-16;
345
346 double tmp2 = mass2 + tmp;
347 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
348 if (q02 < 0) q02 = 1e-16;
349
350 double q = sqrt(q2);
351 double q0 = sqrt(q02);
352 double m = sqrt(sa);
353 double q03 = q0 * q02;
354 double tmp3 = log(mass_param + 2 * q0) + 1.2760418309; // log(mass_2Pion) = 1.2760418309;
355
356 double h = GS1 * q / m * (log(m + 2 * q) + 1.2760418309);
357 double h0 = GS1 * q0 / mass_param * tmp3;
358 double dh = h0 * (0.125 / q02 - 0.5 / mass2) + GS3 / mass2;
359 double d = GS2 / q02 * tmp3 + GS3 * mass_param / q0 - GS4 * mass_param / q03;
360 double f = mass2 / q03 * (q2 * (h - h0) + (mass2 - sa) * q02 * dh);
361
362 a[0] = 1.0 + d * width_param / mass_param;
363 a[1] = 0.0;
364 b[0] = mass2 - sa + width_param * f;
365 b[1] = -mass_param * width_param * widl1(mass2, mass_param, sa, sb, sc, r2);
366 Com_Divide(a, b, prop);
367 }
368
369 void EvtDsToKpipi::Flatte_rhoab(double sa, double sb, double sc, double rho_param[2])
370 {
371 double q = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
372 if (q > 0) {
2
←
Assuming 'q' is <= 0
→
3
←
Taking false branch
→
373 rho_param[0] = 2 * sqrt(q / sa);
374 rho_param[1] = 0;
375 } else if (q < 0) {
4
←
Assuming 'q' is >= 0
→
5
←
Taking false branch
→
376 rho_param[0] = 0;
377 rho_param[1] = 2 * sqrt(-q / sa);
378 }
379 }
6
←
Returning without writing to '*rho_param'
→
380
381
382 void EvtDsToKpipi::propagatorKstr1430(double mass_param, double sx, double* sb, double* sc, double prop[2]) //K*1430 Flatte
383 {
384 double unit[2] = {1.0};
385 double ci[2] = {0, 1};
386 double rho1[2];
387 Flatte_rhoab(sx, sb[0], sc[0], rho1);
1
Calling 'EvtDsToKpipi::Flatte_rhoab'
→
7
←
Returning from 'EvtDsToKpipi::Flatte_rhoab'
→
388 double rho2[2];
389 Flatte_rhoab(sx, sb[1], sc[1], rho2);
390 double gKPi_Kstr1430 = 0.2990, gEtaPK_Kstr1430 = 0.0529;
391 double tmp1[2] = {gKPi_Kstr1430, 0};
392 double tmp11[2];
393 double tmp2[2] = {gEtaPK_Kstr1430, 0};
394 double tmp22[2];
395 Com_Multi(tmp1, rho1, tmp11);
8
←
Calling 'EvtDsToKpipi::Com_Multi'
→
396 Com_Multi(tmp2, rho2, tmp22);
397 double tmp3[2] = {tmp11[0] + tmp22[0], tmp11[1] + tmp22[1]};
398 double tmp31[2];
399 Com_Multi(tmp3, ci, tmp31);
400 double tmp4[2] = {mass_param* mass_param - sx - tmp31[0], -1.0 * tmp31[1]};
401 Com_Divide(unit, tmp4, prop);
402 }
403
404 double EvtDsToKpipi::DDalitz(double P1[4], double P2[4], double P3[4], int Ang, double mass_param)
405 {
406 double pR[4], pD[4];
407 double temp_PDF, v_re;
408 temp_PDF = 0.0;
409 v_re = 0.0;
410 double B[2], s1, s2, s3, sR, sD;
411 for (int i = 0; i < 4; i++) {
412 pR[i] = P1[i] + P2[i];
413 pD[i] = pR[i] + P3[i];
414 }
415 s1 = SCADot(P1, P1);
416 s2 = SCADot(P2, P2);
417 s3 = SCADot(P3, P3);
418 sR = SCADot(pR, pR);
419 sD = SCADot(pD, pD);
420 int GG[4][4];
421 for (int i = 0; i != 4; i++) {
422 for (int j = 0; j != 4; j++) {
423 if (i == j) {
424 if (i == 0) GG[i][j] = 1;
425 else GG[i][j] = -1;
426 } else GG[i][j] = 0;
427 }
428 }
429 if (Ang == 0) {
430 B[0] = 1;
431 B[1] = 1;
432 temp_PDF = 1;
433 }
434 if (Ang == 1) {
435 B[0] = barrier(1, sR, s1, s2, 3.0, mass_param);
436 B[1] = barrier(1, sD, sR, s3, 5.0, mDsM);
437 double t1[4], T1[4];
438 calt1(P1, P2, t1);
439 calt1(pR, P3, T1);
440 temp_PDF = 0;
441 for (int i = 0; i < 4; i++) {
442 temp_PDF += t1[i] * T1[i] * GG[i][i];
443 }
444 }
445 if (Ang == 2) {
446 B[0] = barrier(2, sR, s1, s2, 3.0, mass_param);
447 B[1] = barrier(2, sD, sR, s3, 5.0, mDsM);
448 double t2[4][4], T2[4][4];
449 calt2(P1, P2, t2);
450 calt2(pR, P3, T2);
451 temp_PDF = 0;
452 for (int i = 0; i < 4; i++) {
453 for (int j = 0; j < 4; j++) {
454 temp_PDF += t2[i][j] * T2[j][i] * GG[i][i] * GG[j][j];
455 }
456 }
457 }
458 v_re = temp_PDF * B[0] * B[1];
459 return v_re;
460 }
461
462
463 void EvtDsToKpipi::rhoab(double sa, double sb, double sc, double res[2])
464 {
465 double tmp = sa + sb - sc;
466 double q = 0.25 * tmp * tmp / sa - sb;
467 if (q >= 0) {
468 res[0] = 2.0 * sqrt(q / sa);
469 res[1] = 0.0;
470 } else {
471 res[0] = 0.0;
472 res[1] = 2.0 * sqrt(-q / sa);
473 }
474 }
475 void EvtDsToKpipi::rho4Pi(double sa, double res[2])
476 {
477 double temp = 1.0 - 0.3116765584 / sa;
478 if (temp >= 0) {
479 res[0] = sqrt(temp) / (1.0 + exp(9.8 - 3.5 * sa));
480 res[1] = 0.0;
481 } else {
482 res[0] = 0.0;
483 res[1] = sqrt(-temp) / (1.0 + exp(9.8 - 3.5 * sa));
484 }
485 }
486
487 void EvtDsToKpipi::propagatorsigma500(double sa, double sb, double sc, double prop[2])
488 {
489 double f = 0.5843 + 1.6663 * sa;
490 const double M = 0.9264;
491 const double mass2 = 0.85821696;
492 const double mpi2d2 = 0.00973989245;
493 double g1 = f * (sa - mpi2d2) / (mass2 - mpi2d2) * exp((mass2 - sa) / 1.082);
494 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
495 rhoab(sa, sb, sc, rho1s);
496 rhoab(mass2, sb, sc, rho1M);
497 rho4Pi(sa, rho2s);
498 rho4Pi(mass2, rho2M);
499 Com_Divide(rho1s, rho1M, rho1);
500 Com_Divide(rho2s, rho2M, rho2);
501 double a[2], b[2];
502 a[0] = 1.0;
503 a[1] = 0.0;
504 b[0] = mass2 - sa + M * (g1 * rho1[1] + 0.0024 * rho2[1]);
505 b[1] = -M * (g1 * rho1[0] + 0.0024 * rho2[0]);
506 Com_Divide(a, b, prop);
507 }
508
509
510 void EvtDsToKpipi::calEva(double* K, double* Pi1, double* Pi2, double* mass1, double* width1, double* amp, double* phase, int* g0,
511 int* spin_param, int* modetype_param, int nstates, double& Result)
512 {
513 double P23[4], P13[4];
514 double cof[2], amp_PDF[2], PDF[2];
515 double s13, s23;
516 for (int i = 0; i < 4; i++) {
517 P13[i] = K[i] + Pi2[i];
518 P23[i] = Pi1[i] + Pi2[i];
519 }
520 s13 = SCADot(P13, P13);
521 s23 = SCADot(P23, P23);
522 double s1, s2, s3;
523 s1 = SCADot(K, K);
524 s2 = SCADot(Pi1, Pi1);
525 s3 = SCADot(Pi2, Pi2);
526 double pro[2], temp_PDF, amp_tmp[2];
527 amp_PDF[0] = 0;
528 amp_PDF[1] = 0;
529 PDF[0] = 0;
530 PDF[1] = 0;
531 amp_tmp[0] = 0;
532 amp_tmp[1] = 0;
533 double rRess = 9.0;
534 for (int i = 0; i < nstates; i++) {
535 amp_tmp[0] = 0;
536 amp_tmp[1] = 0;
537 cof[0] = amp[i] * cos(phase[i]);
538 cof[1] = amp[i] * sin(phase[i]);
539 temp_PDF = 0;
540
541 if (modetype_param[i] == 13) {
542 temp_PDF = DDalitz(K, Pi2, Pi1, spin_param[i], mass1[i]);
543 if (g0[i] == 1) propagatorRBW(mass1[i], width1[i], s13, mKa2, mPi2, rRess, spin_param[i], pro);
544 if (g0[i] == 2) { // K*1430 Flatte
545 double skm2[2] = {s1, 0.95778 * 0.95778};
546 double spi2[2] = {s3, mass_Kaon * mass_Kaon};
547 propagatorKstr1430(mass1[i], s13, skm2, spi2, pro);
548 }
549 if (g0[i] == 0) {
550 pro[0] = 1;
551 pro[1] = 0;
552 }
553 amp_tmp[0] = temp_PDF * pro[0];
554 amp_tmp[1] = temp_PDF * pro[1];
555 }
556
557 if (modetype_param[i] == 23) {
558 temp_PDF = DDalitz(Pi1, Pi2, K, spin_param[i], mass1[i]);
559 if (g0[i] == 4) propagatorFlatte(mass1[i], width1[i], s23, pro); // Only for f0(980)
560 if (g0[i] == 3) propagatorCBW(mass1[i], width1[i], s23, pro);
561 if (g0[i] == 2) propagatorRBW(mass1[i], width1[i], s23, mPi2, mPi2, rRess, spin_param[i], pro);
562 if (g0[i] == 1) propagatorGS(mass1[i], width1[i], s23, mPi2, mPi2, 9.0, pro);
563 if (g0[i] == 500) propagatorsigma500(s23, s2, s3, pro);
564 if (g0[i] == 0) {
565 pro[0] = 1;
566 pro[1] = 0;
567 }
568 amp_tmp[0] = temp_PDF * pro[0];
569 amp_tmp[1] = temp_PDF * pro[1];
570 }
571 Com_Multi(amp_tmp, cof, amp_PDF);
572 PDF[0] += amp_PDF[0];
573 PDF[1] += amp_PDF[1];
574 }
575
576 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
577 if (value <= 0) value = 1e-20;
578 Result = value;
579 }
580
581} // Belle2 namespace