Belle II Software development
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 {
35
38
39 EvtDsToKpipi::~EvtDsToKpipi() {}
40
41 std::string EvtDsToKpipi::getName()
42 {
43 return "DsToKpipi";
44 }
45
46 EvtDecayBase* EvtDsToKpipi::clone()
47 {
48 return new EvtDsToKpipi;
49 }
50
51 void EvtDsToKpipi::init()
52 {
53 checkNArg(0);
54 checkNDaug(3);
55 checkSpinParent(EvtSpinType::SCALAR);
56 checkSpinDaughter(0, EvtSpinType::SCALAR);
57 checkSpinDaughter(1, EvtSpinType::SCALAR);
58 checkSpinDaughter(2, EvtSpinType::SCALAR);
59
60 phi[0] = 0;
61 rho[0] = 1;
62 phi[1] = 0;
63 rho[1] = 0;
64 phi[2] = 0;
65 rho[2] = 0;
66 phi[3] = 0;
67 rho[3] = 0;
68 phi[4] = 0;
69 rho[4] = 0;
70 phi[5] = 0;
71 rho[5] = 0;
72 phi[6] = 0;
73 rho[6] = 0;
74 phi[7] = 0;
75 rho[7] = 0;
76
77 phi[1] = -3.47995752;
78 phi[2] = -1.249864467;
79 phi[3] = -0.3067504308;
80 phi[4] = 0.9229242379;
81 phi[5] = -3.278567926;
82 phi[6] = -0.6346408622;
83 phi[7] = 1.762377475;
84
85 rho[1] = 2.463898984;
86 rho[2] = 0.7361782665;
87 rho[3] = 1.90216812;
88 rho[4] = 2.140687169;
89 rho[5] = 0.914684056;
90 rho[6] = 1.116206881;
91 rho[7] = 1.483440555;
92
93 modetype[0] = 23;
94 modetype[1] = 23;
95 modetype[2] = 23;
96 modetype[3] = 23;
97 modetype[4] = 23;
98 modetype[5] = 13;
99 modetype[6] = 13;
100 modetype[7] = 13;
101
102 width[0] = 0.1478;
103 width[1] = 0.400;
104 width[2] = 0.100;
105 width[3] = 0.265;
106 width[4] = 0.270;
107 width[5] = 0.0473;
108 width[6] = 0.232;
109 width[7] = 0.270;
110
111 mass[0] = 0.77526;
112 mass[1] = 1.465;
113 mass[2] = 0.965;
114 mass[3] = 1.35;
115 mass[4] = 1.425;
116 mass[5] = 0.89555;
117 mass[6] = 1.414;
118 mass[7] = 1.432787726;
119
120 mDsM = 1.9683;
121 mD = 1.86486;
122 mDs = 1.9683;
123 rD = 25.0;
124 metap = 0.95778;
125 mkstr = 0.89594;
126 mk0 = 0.497614;
127 mass_Kaon = 0.49368;
128 mass_Pion = 0.13957;
129 mass_Pi0 = 0.1349766;
130 math_pi = 3.1415926;
131 mKa2 = 0.24371994;
132 mPi2 = 0.01947978;
133 mPi = 0.13957;
134 mKa = 0.49368;
135 ma0 = 0.99;
136 Ga0 = 0.0756;
137 meta = 0.547862;
138
139 GS1 = 0.636619783;
140 GS2 = 0.01860182466;
141 GS3 = 0.1591549458;
142 GS4 = 0.00620060822;
143
144 int GG[4][4] = { {1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1} };
145 for (int i = 0; i < 4; i++) {
146 for (int j = 0; j < 4; j++) {
147 G[i][j] = GG[i][j];
148 }
149 }
150 }
151
152 void EvtDsToKpipi::initProbMax()
153 {
154 setProbMax(825.0);
155 }
156
157 void EvtDsToKpipi::decay(EvtParticle* p)
158 {
159 p->initializePhaseSpace(getNDaug(), getDaugs());
160 EvtVector4R D1 = p->getDaug(0)->getP4();
161 EvtVector4R D2 = p->getDaug(1)->getP4();
162 EvtVector4R D3 = p->getDaug(2)->getP4();
163
164 double P1[4], P2[4], P3[4];
165 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
166 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
167 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
168
169 double value;
170 int g0[8] = {1, 1, 4, 2, 500, 1, 1, 2};
171 int spin_param[8] = {1, 1, 0, 0, 0, 1, 1, 0};
172 int nstates = 8;
173 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin_param, modetype, nstates, value);
174
175 setProb(value);
176
177 return ;
178 }
179
180 void EvtDsToKpipi::Com_Multi(double a1[2], double a2[2], double res[2])
181 {
182 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
183 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
184 }
185 void EvtDsToKpipi::Com_Divide(double a1[2], double a2[2], double res[2])
186 {
187 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
188 res[0] = (a1[0] * a2[0] + a1[1] * a2[1]) / tmp;
189 res[1] = (a1[1] * a2[0] - a1[0] * a2[1]) / tmp;
190 }
191
192 double EvtDsToKpipi::SCADot(double a1[4], double a2[4])
193 {
194 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
195 return _cal;
196 }
197 double EvtDsToKpipi::barrier(int l, double sa, double sb, double sc, double r, double mass_param)
198 {
199 double q = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
200 if (q < 0) q = 1e-16;
201 double z;
202 z = q * r * r;
203 double sa0;
204 sa0 = mass_param * mass_param;
205 double q0 = (sa0 + sb - sc) * (sa0 + sb - sc) / (4 * sa0) - sb;
206 if (q0 < 0) q0 = 1e-16;
207 double z0 = q0 * r * r;
208 double F = 0.0;
209 if (l == 0) F = 1;
210 if (l == 1) F = sqrt((1 + z0) / (1 + z));
211 if (l == 2) F = sqrt((9 + 3 * z0 + z0 * z0) / (9 + 3 * z + z * z));
212 return F;
213 }
214 void EvtDsToKpipi::calt1(double daug1[4], double daug2[4], double t1[4])
215 {
216 double p, pq, tmp;
217 double pa[4], qa[4];
218 for (int i = 0; i < 4; i++) {
219 pa[i] = daug1[i] + daug2[i];
220 qa[i] = daug1[i] - daug2[i];
221 }
222 p = SCADot(pa, pa);
223 pq = SCADot(pa, qa);
224 tmp = pq / p;
225 for (int i = 0; i < 4; i++) {
226 t1[i] = qa[i] - tmp * pa[i];
227 }
228 }
229 void EvtDsToKpipi::calt2(double daug1[4], double daug2[4], double t2[4][4])
230 {
231 double p, r;
232 double pa[4], t1[4];
233 calt1(daug1, daug2, t1);
234 r = SCADot(t1, t1) / 3.0;
235 for (int i = 0; i < 4; i++) {
236 pa[i] = daug1[i] + daug2[i];
237 }
238 p = SCADot(pa, pa);
239 for (int i = 0; i < 4; i++) {
240 for (int j = 0; j < 4; j++) {
241 t2[i][j] = t1[i] * t1[j] - r * (G[i][j] - pa[i] * pa[j] / p);
242 }
243 }
244 }
245
246//-------------------prop--------------------------------------------
247
248 void EvtDsToKpipi::propagatorCBW(double mass_param, double width_param, double sx, double prop[2])
249 {
250 double a[2], b[2];
251 a[0] = 1;
252 a[1] = 0;
253 b[0] = mass_param * mass_param - sx;
254 b[1] = -mass_param * width_param;
255 Com_Divide(a, b, prop);
256 }
257 double EvtDsToKpipi::wid(double mass2, double mass_param, double sa, double sb, double sc, double r2, int l)
258 {
259 double widm = 0.;
260 double m = sqrt(sa);
261 double tmp = sb - sc;
262 double tmp1 = sa + tmp;
263 double q = 0.25 * tmp1 * tmp1 / sa - sb;
264 if (q < 0) q = 1e-16;
265 double tmp2 = mass2 + tmp;
266 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
267 if (q0 < 0) q0 = 1e-16;
268 double z = q * r2;
269 double z0 = q0 * r2;
270 double t = q / q0;
271 if (l == 0) {widm = sqrt(t) * mass_param / m;}
272 else if (l == 1) {widm = t * sqrt(t) * mass_param / m * (1 + z0) / (1 + z);}
273 else if (l == 2) {widm = t * t * sqrt(t) * mass_param / m * (9 + 3 * z0 + z0 * z0) / (9 + 3 * z + z * z);}
274 return widm;
275 }
276
277 double EvtDsToKpipi::widl1(double mass2, double mass_param, double sa, double sb, double sc, double r2)
278 {
279 double widm = 0.;
280 double m = sqrt(sa);
281 double tmp = sb - sc;
282 double tmp1 = sa + tmp;
283 double q = 0.25 * tmp1 * tmp1 / sa - sb;
284 if (q < 0) q = 1e-16;
285 double tmp2 = mass2 + tmp;
286 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
287 if (q0 < 0) q0 = 1e-16;
288 double z = q * r2;
289 double z0 = q0 * r2;
290 double F = (1 + z0) / (1 + z);
291 double t = q / q0;
292 widm = t * sqrt(t) * mass_param / m * F;
293 return widm;
294 }
295 void EvtDsToKpipi::propagatorRBW(double mass_param, double width_param, double sa, double sb, double sc, double r2, int l,
296 double prop[2])
297 {
298 double a[2], b[2];
299 double mass2 = mass_param * mass_param;
300
301 a[0] = 1;
302 a[1] = 0;
303 b[0] = mass2 - sa;
304 b[1] = -mass_param * width_param * wid(mass2, mass_param, sa, sb, sc, r2, l);
305 Com_Divide(a, b, prop);
306 }
307 void EvtDsToKpipi::propagatorFlatte(double mass_param, double width_param __attribute__((unused)), double sa, double prop[2])
308 {
309
310 double q2_Pi, q2_Ka;
311 double rhoPi[2] = {0, 0}, rhoKa[2] = {0, 0};
312
313 q2_Pi = 0.25 * sa - mPi * mPi;
314 q2_Ka = 0.25 * sa - mKa * mKa;
315
316 if (q2_Pi > 0) {
317 rhoPi[0] = 2.0 * sqrt(q2_Pi / sa);
318 rhoPi[1] = 0.0;
319 }
320 if (q2_Pi <= 0) {
321 rhoPi[0] = 0.0;
322 rhoPi[1] = 2.0 * sqrt(-q2_Pi / sa);
323 }
324
325 if (q2_Ka > 0) {
326 rhoKa[0] = 2.0 * sqrt(q2_Ka / sa);
327 rhoKa[1] = 0.0;
328 }
329 if (q2_Ka <= 0) {
330 rhoKa[0] = 0.0;
331 rhoKa[1] = 2.0 * sqrt(-q2_Ka / sa);
332 }
333
334 double a[2], b[2];
335 a[0] = 1;
336 a[1] = 0;
337 b[0] = mass_param * mass_param - sa + 0.165 * rhoPi[1] + 0.69465 * rhoKa[1];
338 b[1] = - (0.165 * rhoPi[0] + 0.69465 * rhoKa[0]);
339 Com_Divide(a, b, prop);
340 }
341 void EvtDsToKpipi::propagatorGS(double mass_param, double width_param, double sa, double sb, double sc, double r2, double prop[2])
342 {
343 double a[2], b[2];
344 double mass2 = mass_param * mass_param;
345 double tmp = sb - sc;
346 double tmp1 = sa + tmp;
347 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
348 if (q2 < 0) q2 = 1e-16;
349
350 double tmp2 = mass2 + tmp;
351 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
352 if (q02 < 0) q02 = 1e-16;
353
354 double q = sqrt(q2);
355 double q0 = sqrt(q02);
356 double m = sqrt(sa);
357 double q03 = q0 * q02;
358 double tmp3 = log(mass_param + 2 * q0) + 1.2760418309; // log(mass_2Pion) = 1.2760418309;
359
360 double h = GS1 * q / m * (log(m + 2 * q) + 1.2760418309);
361 double h0 = GS1 * q0 / mass_param * tmp3;
362 double dh = h0 * (0.125 / q02 - 0.5 / mass2) + GS3 / mass2;
363 double d = GS2 / q02 * tmp3 + GS3 * mass_param / q0 - GS4 * mass_param / q03;
364 double f = mass2 / q03 * (q2 * (h - h0) + (mass2 - sa) * q02 * dh);
365
366 a[0] = 1.0 + d * width_param / mass_param;
367 a[1] = 0.0;
368 b[0] = mass2 - sa + width_param * f;
369 b[1] = -mass_param * width_param * widl1(mass2, mass_param, sa, sb, sc, r2);
370 Com_Divide(a, b, prop);
371 }
372
373 void EvtDsToKpipi::Flatte_rhoab(double sa, double sb, double sc, double rho_param[2])
374 {
375 double q = (sa + sb - sc) * (sa + sb - sc) / (4 * sa) - sb;
376 if (q > 0) {
377 rho_param[0] = 2 * sqrt(q / sa);
378 rho_param[1] = 0;
379 } else if (q < 0) {
380 rho_param[0] = 0;
381 rho_param[1] = 2 * sqrt(-q / sa);
382 }
383 }
384
385
386 void EvtDsToKpipi::propagatorKstr1430(double mass_param, double sx, double* sb, double* sc, double prop[2]) //K*1430 Flatte
387 {
388 double unit[2] = {1.0};
389 double ci[2] = {0, 1};
390 double rho1[2];
391 Flatte_rhoab(sx, sb[0], sc[0], rho1);
392 double rho2[2];
393 Flatte_rhoab(sx, sb[1], sc[1], rho2);
394 double gKPi_Kstr1430 = 0.2990, gEtaPK_Kstr1430 = 0.0529;
395 double tmp1[2] = {gKPi_Kstr1430, 0};
396 double tmp11[2];
397 double tmp2[2] = {gEtaPK_Kstr1430, 0};
398 double tmp22[2];
399 Com_Multi(tmp1, rho1, tmp11);
400 Com_Multi(tmp2, rho2, tmp22);
401 double tmp3[2] = {tmp11[0] + tmp22[0], tmp11[1] + tmp22[1]};
402 double tmp31[2];
403 Com_Multi(tmp3, ci, tmp31);
404 double tmp4[2] = {mass_param* mass_param - sx - tmp31[0], -1.0 * tmp31[1]};
405 Com_Divide(unit, tmp4, prop);
406 }
407
408 double EvtDsToKpipi::DDalitz(double P1[4], double P2[4], double P3[4], int Ang, double mass_param)
409 {
410 double pR[4], pD[4];
411 double temp_PDF, v_re;
412 temp_PDF = 0.0;
413 v_re = 0.0;
414 double B[2], s1, s2, s3, sR, sD;
415 for (int i = 0; i < 4; i++) {
416 pR[i] = P1[i] + P2[i];
417 pD[i] = pR[i] + P3[i];
418 }
419 s1 = SCADot(P1, P1);
420 s2 = SCADot(P2, P2);
421 s3 = SCADot(P3, P3);
422 sR = SCADot(pR, pR);
423 sD = SCADot(pD, pD);
424 int GG[4][4];
425 for (int i = 0; i != 4; i++) {
426 for (int j = 0; j != 4; j++) {
427 if (i == j) {
428 if (i == 0) GG[i][j] = 1;
429 else GG[i][j] = -1;
430 } else GG[i][j] = 0;
431 }
432 }
433 if (Ang == 0) {
434 B[0] = 1;
435 B[1] = 1;
436 temp_PDF = 1;
437 }
438 if (Ang == 1) {
439 B[0] = barrier(1, sR, s1, s2, 3.0, mass_param);
440 B[1] = barrier(1, sD, sR, s3, 5.0, mDsM);
441 double t1[4], T1[4];
442 calt1(P1, P2, t1);
443 calt1(pR, P3, T1);
444 temp_PDF = 0;
445 for (int i = 0; i < 4; i++) {
446 temp_PDF += t1[i] * T1[i] * GG[i][i];
447 }
448 }
449 if (Ang == 2) {
450 B[0] = barrier(2, sR, s1, s2, 3.0, mass_param);
451 B[1] = barrier(2, sD, sR, s3, 5.0, mDsM);
452 double t2[4][4], T2[4][4];
453 calt2(P1, P2, t2);
454 calt2(pR, P3, T2);
455 temp_PDF = 0;
456 for (int i = 0; i < 4; i++) {
457 for (int j = 0; j < 4; j++) {
458 temp_PDF += t2[i][j] * T2[j][i] * GG[i][i] * GG[j][j];
459 }
460 }
461 }
462 v_re = temp_PDF * B[0] * B[1];
463 return v_re;
464 }
465
466
467 void EvtDsToKpipi::rhoab(double sa, double sb, double sc, double res[2])
468 {
469 double tmp = sa + sb - sc;
470 double q = 0.25 * tmp * tmp / sa - sb;
471 if (q >= 0) {
472 res[0] = 2.0 * sqrt(q / sa);
473 res[1] = 0.0;
474 } else {
475 res[0] = 0.0;
476 res[1] = 2.0 * sqrt(-q / sa);
477 }
478 }
479 void EvtDsToKpipi::rho4Pi(double sa, double res[2])
480 {
481 double temp = 1.0 - 0.3116765584 / sa;
482 if (temp >= 0) {
483 res[0] = sqrt(temp) / (1.0 + exp(9.8 - 3.5 * sa));
484 res[1] = 0.0;
485 } else {
486 res[0] = 0.0;
487 res[1] = sqrt(-temp) / (1.0 + exp(9.8 - 3.5 * sa));
488 }
489 }
490
491 void EvtDsToKpipi::propagatorsigma500(double sa, double sb, double sc, double prop[2])
492 {
493 double f = 0.5843 + 1.6663 * sa;
494 const double M = 0.9264;
495 const double mass2 = 0.85821696;
496 const double mpi2d2 = 0.00973989245;
497 double g1 = f * (sa - mpi2d2) / (mass2 - mpi2d2) * exp((mass2 - sa) / 1.082);
498 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
499 rhoab(sa, sb, sc, rho1s);
500 rhoab(mass2, sb, sc, rho1M);
501 rho4Pi(sa, rho2s);
502 rho4Pi(mass2, rho2M);
503 Com_Divide(rho1s, rho1M, rho1);
504 Com_Divide(rho2s, rho2M, rho2);
505 double a[2], b[2];
506 a[0] = 1.0;
507 a[1] = 0.0;
508 b[0] = mass2 - sa + M * (g1 * rho1[1] + 0.0024 * rho2[1]);
509 b[1] = -M * (g1 * rho1[0] + 0.0024 * rho2[0]);
510 Com_Divide(a, b, prop);
511 }
512
513
514 void EvtDsToKpipi::calEva(double* K, double* Pi1, double* Pi2, double* mass1, double* width1, double* amp, double* phase, int* g0,
515 int* spin_param, int* modetype_param, int nstates, double& Result)
516 {
517 double P23[4], P13[4];
518 double cof[2], amp_PDF[2], PDF[2];
519 double s13, s23;
520 for (int i = 0; i < 4; i++) {
521 P13[i] = K[i] + Pi2[i];
522 P23[i] = Pi1[i] + Pi2[i];
523 }
524 s13 = SCADot(P13, P13);
525 s23 = SCADot(P23, P23);
526 double s1, s2, s3;
527 s1 = SCADot(K, K);
528 s2 = SCADot(Pi1, Pi1);
529 s3 = SCADot(Pi2, Pi2);
530 double pro[2], temp_PDF, amp_tmp[2];
531 amp_PDF[0] = 0;
532 amp_PDF[1] = 0;
533 PDF[0] = 0;
534 PDF[1] = 0;
535 amp_tmp[0] = 0;
536 amp_tmp[1] = 0;
537 double rRess = 9.0;
538 for (int i = 0; i < nstates; i++) {
539 amp_tmp[0] = 0;
540 amp_tmp[1] = 0;
541 cof[0] = amp[i] * cos(phase[i]);
542 cof[1] = amp[i] * sin(phase[i]);
543 temp_PDF = 0;
544
545 if (modetype_param[i] == 13) {
546 temp_PDF = DDalitz(K, Pi2, Pi1, spin_param[i], mass1[i]);
547 if (g0[i] == 1) propagatorRBW(mass1[i], width1[i], s13, mKa2, mPi2, rRess, spin_param[i], pro);
548 if (g0[i] == 2) { // K*1430 Flatte
549 double skm2[2] = {s1, 0.95778 * 0.95778};
550 double spi2[2] = {s3, mass_Kaon * mass_Kaon};
551 propagatorKstr1430(mass1[i], s13, skm2, spi2, pro);
552 }
553 if (g0[i] == 0) {
554 pro[0] = 1;
555 pro[1] = 0;
556 }
557 amp_tmp[0] = temp_PDF * pro[0];
558 amp_tmp[1] = temp_PDF * pro[1];
559 }
560
561 if (modetype_param[i] == 23) {
562 temp_PDF = DDalitz(Pi1, Pi2, K, spin_param[i], mass1[i]);
563 if (g0[i] == 4) propagatorFlatte(mass1[i], width1[i], s23, pro); // Only for f0(980)
564 if (g0[i] == 3) propagatorCBW(mass1[i], width1[i], s23, pro);
565 if (g0[i] == 2) propagatorRBW(mass1[i], width1[i], s23, mPi2, mPi2, rRess, spin_param[i], pro);
566 if (g0[i] == 1) propagatorGS(mass1[i], width1[i], s23, mPi2, mPi2, 9.0, pro);
567 if (g0[i] == 500) propagatorsigma500(s23, s2, s3, pro);
568 if (g0[i] == 0) {
569 pro[0] = 1;
570 pro[1] = 0;
571 }
572 amp_tmp[0] = temp_PDF * pro[0];
573 amp_tmp[1] = temp_PDF * pro[1];
574 }
575 Com_Multi(amp_tmp, cof, amp_PDF);
576 PDF[0] += amp_PDF[0];
577 PDF[1] += amp_PDF[1];
578 }
579
580 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
581 if (value <= 0) value = 1e-20;
582 Result = value;
583 }
584
586} // Belle2 namespace
#define K(x)
macro autogenerated by FFTW
#define B2_EVTGEN_REGISTER_MODEL(classname)
Class to register B2_EVTGEN_REGISTER_MODEL.
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
Abstract base class for different kinds of events.