Belle II Software development
EvtB0toK0K0K0.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <generators/evtgen/EvtGenModelRegister.h>
10
11#include <stdlib.h>
12#include <cmath>
13#include <cstring>
14#include "EvtGenBase/EvtRandom.hh"
15#include "EvtGenBase/EvtParticle.hh"
16#include "EvtGenBase/EvtGenKine.hh"
17#include "EvtGenBase/EvtPDL.hh"
18#include "EvtGenBase/EvtReport.hh"
19#include "EvtGenBase/EvtConst.hh"
20#include "EvtGenBase/EvtId.hh"
21
22#include "generators/evtgen/models/EvtB0toK0K0K0.h"
23
24using std::endl;
25using namespace std::complex_literals;
26
27namespace Belle2 {
34 B2_EVTGEN_REGISTER_MODEL(EvtB0toK0K0K0);
35
36 EvtB0toK0K0K0::~EvtB0toK0K0K0() {}
37
38 std::string EvtB0toK0K0K0::getName()
39 {
40 return "B0TOK0K0K0";
41 }
42
43 EvtDecayBase* EvtB0toK0K0K0::clone()
44 {
45 return new EvtB0toK0K0K0;
46 }
47
48 void EvtB0toK0K0K0::decay(EvtParticle* p)
49 {
50 // follow PhysRevD.85.054023
51
52 int Ntry = 0;
53 const int max_Ntry = 10000;
54
55 while (true) {
56 p->initializePhaseSpace(getNDaug(), getDaugs());
57
58 EvtParticle* NeutralKaon_1 = p->getDaug(0);
59 EvtParticle* NeutralKaon_2 = p->getDaug(1);
60 EvtParticle* NeutralKaon_3 = p->getDaug(2);
61
62 EvtVector4R mom_NeutralKaon_1 = NeutralKaon_1->getP4();
63 EvtVector4R mom_NeutralKaon_2 = NeutralKaon_2->getP4();
64 EvtVector4R mom_NeutralKaon_3 = NeutralKaon_3->getP4();
65
66 EvtVector4R mom_sum_12 = mom_NeutralKaon_1 + mom_NeutralKaon_2;
67 EvtVector4R mom_sum_13 = mom_NeutralKaon_1 + mom_NeutralKaon_3;
68 EvtVector4R mom_sum_23 = mom_NeutralKaon_2 + mom_NeutralKaon_3;
69
70 double s12 = mom_sum_12.mass2();
71 double s13 = mom_sum_13.mass2();
72 double s23 = mom_sum_23.mass2();
73
74 // follow the index of PhysRevD.85.054023
75 double smax = std::max(std::max(s12, s13), s23);
76 double smin = std::min(std::min(s12, s13), s23);
77
78 // the maximum value of probability is obtained by brute fource method. scan all region and find.
79 const double Probability_max = 0.409212;
80 double Probability_value = Probability(smax, smin);
81 double Probability_ratio = Probability_value / Probability_max;
82
83 double xbox = EvtRandom::Flat(0.0, 1.0);
84
85 if (xbox < Probability_ratio) break;
86 Ntry++;
87 if (Ntry > max_Ntry) {
88 std::cout << "try to set the kinematics more than 10000 times. abort.\n";
89 ::abort();
90 }
91 }
92
93
94 return;
95
96 }
97
98
99 void EvtB0toK0K0K0::initProbMax()
100 {
101 noProbMax();
102 }
103
104
105 void EvtB0toK0K0K0::init()
106 {
107
108 // check that there are no arguments
109 checkNArg(0);
110
111 // there should be three daughters: K+ KL0 KL0
112 checkNDaug(3);
113
114 // Check that the daughters are correct
115 EvtId KaonZeroType_1 = getDaug(0);
116 EvtId KaonZeroType_2 = getDaug(1);
117 EvtId KaonZeroType_3 = getDaug(2);
118
119 int KaonZerotyp = 0;
120
121 if ((KaonZeroType_1 == EvtPDL::getId("K0")) ||
122 (KaonZeroType_1 == EvtPDL::getId("anti-K0")) ||
123 (KaonZeroType_1 == EvtPDL::getId("K_S0")) ||
124 (KaonZeroType_1 == EvtPDL::getId("K_L0"))) KaonZerotyp++;
125 if ((KaonZeroType_2 == EvtPDL::getId("K0")) ||
126 (KaonZeroType_2 == EvtPDL::getId("anti-K0")) ||
127 (KaonZeroType_2 == EvtPDL::getId("K_S0")) ||
128 (KaonZeroType_2 == EvtPDL::getId("K_L0"))) KaonZerotyp++;
129 if ((KaonZeroType_3 == EvtPDL::getId("K0")) ||
130 (KaonZeroType_3 == EvtPDL::getId("anti-K0")) ||
131 (KaonZeroType_3 == EvtPDL::getId("K_S0")) ||
132 (KaonZeroType_3 == EvtPDL::getId("K_L0"))) KaonZerotyp++;
133
134 if (KaonZerotyp != 3) {
135
136 std::cout << "All daughters should be K0, anti-K0, K_L0, or K_S0.\n";
137 ::abort();
138 }
139
140 // initialize q0 and pstar0
141 GetMasses();
142 GetZeros();
143
144 }
145
146 double EvtB0toK0K0K0::Probability(double s13, double s23)
147 {
148 std::complex<double> total_amplitude = Amplitude(s13, s23, "f980") + Amplitude(s13, s23, "f1710") + Amplitude(s13, s23,
149 "f2010") + Amplitude(s13, s23, "chic0") + Amplitude(s13, s23, "NR");
150 return std::abs(total_amplitude) * std::abs(total_amplitude);
151 }
152
153
154 std::complex<double> EvtB0toK0K0K0::Amplitude(double s13, double s23, const char* resonance)
155 {
156
157 // this factors scale each resonance. Target FFs are written in the paper
158 double MagicNumber_f980 = std::sqrt(0.44 / 30.222805945);
159 double MagicNumber_f1710 = std::sqrt(0.07 / 136.555962030);
160 double MagicNumber_f2010 = std::sqrt(0.09 / 310762.546712675);
161 double MagicNumber_chic0 = std::sqrt(0.07 / 677.282967269);
162 double MagicNumber_NR = std::sqrt(2.16 / 51.262701105);
163
164 double s12 = mB0 * mB0 + mKS0 * mKS0 + mKL0 * mKL0 + mKL0 * mKL0 - s13 - s23;
165
166 if (strcmp(resonance, "f980") == 0) {
167 std::complex<double> a;
168 a = c_f980 * std::exp(1i * phi_f980);
169 return MagicNumber_f980 * a * (DynamicalAmplitude(s13, s23, resonance) + DynamicalAmplitude(s12, s23,
170 resonance) + DynamicalAmplitude(s12, s13, resonance));
171 } else if (strcmp(resonance, "f1710") == 0) {
172 std::complex<double> a;
173 a = c_f1710 * std::exp(1i * phi_f1710);
174 return MagicNumber_f1710 * a * (DynamicalAmplitude(s13, s23, resonance) + DynamicalAmplitude(s12, s23,
175 resonance) + DynamicalAmplitude(s12, s13, resonance));
176 } else if (strcmp(resonance, "f2010") == 0) {
177 std::complex<double> a;
178 a = c_f2010 * std::exp(1i * phi_f2010);
179 return MagicNumber_f2010 * a * (DynamicalAmplitude(s13, s23, resonance) + DynamicalAmplitude(s12, s23,
180 resonance) + DynamicalAmplitude(s12, s13, resonance));
181 } else if (strcmp(resonance, "chic0") == 0) {
182 std::complex<double> a;
183 a = c_chic0 * std::exp(1i * phi_chic0);
184 return MagicNumber_chic0 * a * (DynamicalAmplitude(s13, s23, resonance) + DynamicalAmplitude(s12, s23,
185 resonance) + DynamicalAmplitude(s12, s13, resonance));
186 } else if (strcmp(resonance, "NR") == 0) {
187
188 std::complex<double> a;
189 a = c_NR * std::exp(1i * phi_NR);
190
191 return MagicNumber_NR * a * (std::exp(alpha_NR * s13) + std::exp(alpha_NR * s12) + std::exp(alpha_NR * s23));
192 } else {
193 printf("[Amplitude] unsupported resonance\n");
194 exit(1);
195 }
196
197 }
198
199 std::complex<double> EvtB0toK0K0K0::DynamicalAmplitude(double s13, double s23, const char* resonance) // resonance: 1+2
200 {
201
202 if (strcmp(resonance, "f980") == 0) {
203 return Flatte(s13, s23, resonance) * BlattWeisskopf_pstarrprime(s13, s23, resonance) * BlattWeisskopf_qr(s13, s23,
204 resonance) * Zemach(s13, s23, resonance);
205 }
206 if (strcmp(resonance, "f1710") == 0) {
207 return RBW(s13, s23, resonance) * BlattWeisskopf_pstarrprime(s13, s23, resonance) * BlattWeisskopf_qr(s13, s23,
208 resonance) * Zemach(s13, s23, resonance);
209 } else if (strcmp(resonance, "f2010") == 0) {
210 return RBW(s13, s23, resonance) * BlattWeisskopf_pstarrprime(s13, s23, resonance) * BlattWeisskopf_qr(s13, s23,
211 resonance) * Zemach(s13, s23, resonance);
212 } else if (strcmp(resonance, "chic0") == 0) {
213 return RBW(s13, s23, resonance) * BlattWeisskopf_pstarrprime(s13, s23, resonance) * BlattWeisskopf_qr(s13, s23,
214 resonance) * Zemach(s13, s23, resonance);
215 } else {
216 printf("[Amplitude] unsupported resonance\n");
217 exit(1);
218 }
219
220 }
221
222 std::complex<double> EvtB0toK0K0K0::Flatte(double s13, double s23, const char* resonance)
223 {
224 double m0;
225 double gpi;
226 double gK;
227
228 if (strcmp(resonance, "f980") == 0) {
229 m0 = m0_f980;
230 gpi = gpi_f980;
231 gK = gK_f980;
232 } else {
233 printf("[Flatte] unsupported resonance\n");
234 exit(1);
235 }
236
237 double m = Calculate_m(s13, s23);
238
239 double Gammapipi = gpi * ((1.0 / 3.0) * std::sqrt(1 - 4.0 * mpi0 * mpi0 / (m * m)) + (2.0 / 3.0) * std::sqrt(
240 1 - 4.0 * mpic * mpic / (m * m)));
241 double GammaKK = gK * (0.5 * std::sqrt(1 - 4.0 * mKp * mKp / (m * m)) + 0.5 * std::sqrt(1 - 4.0 * mK0 * mK0 / (m * m)));
242
243 return 1.0 / ((m0 * m0 - m * m) - 1i * (Gammapipi + GammaKK));
244 }
245
246 std::complex<double> EvtB0toK0K0K0::RBW(double s13, double s23, const char* resonance)
247 {
248 double m0;
249
250 double m = Calculate_m(s13, s23);
251
252 if (strcmp(resonance, "f1710") == 0) {
253 m0 = m0_f1710;
254 } else if (strcmp(resonance, "f2010") == 0) {
255 m0 = m0_f2010;
256 } else if (strcmp(resonance, "chic0") == 0) {
257 m0 = m0_chic0;
258 } else {
259 printf("[RBW] unsupported resonance\n");
260 exit(1);
261 }
262
263 return 1.0 / (m0 * m0 - m * m - 1i * m0 * MassDepWidth(s13, s23, resonance));
264 }
265
266 double EvtB0toK0K0K0::MassDepWidth(double s13, double s23, const char* resonance)
267 {
268 // Gamma(m) when s12 and s23
269 int spin = -1;
270 double m0;
271 double q0;
272 double Gamma0;
273
274 double q_mag = Calculate_q_mag(s13, s23);
275 double m = Calculate_m(s13, s23);
276
277 if (strcmp(resonance, "f1710") == 0) {
278 spin = spin_f1710;
279 m0 = m0_f1710;
280 q0 = q0_f1710;
281 Gamma0 = Gamma0_f1710;
282 } else if (strcmp(resonance, "f2010") == 0) {
283 spin = spin_f2010;
284 m0 = m0_f2010;
285 q0 = q0_f2010;
286 Gamma0 = Gamma0_f2010;
287 } else if (strcmp(resonance, "chic0") == 0) {
288 spin = spin_chic0;
289 m0 = m0_chic0;
290 q0 = q0_chic0;
291 Gamma0 = Gamma0_chic0;
292 } else {
293 printf("[MassDepWidth] unsupported resonance\n");
294 exit(1);
295 }
296
297 return Gamma0 * std::pow(q_mag / q0, 2 * spin + 1) * (m0 / m) * std::pow(BlattWeisskopf_qr(s13, s23, resonance), 2);
298 }
299
300 double EvtB0toK0K0K0::BlattWeisskopf_pstarrprime(double s13, double s23, const char* resonance)
301 {
302 // X_L(|p*|r') when s12 and s23
303 int spin = -1;
304 double z0 = -1;
305 double z = Calculate_pstar_mag(s13, s23) * rprime;
306
307 if (strcmp(resonance, "f980") == 0) {
308 spin = spin_f980;
309 z0 = -1;
310 } else if (strcmp(resonance, "f1710") == 0) {
311 spin = spin_f1710;
312 z0 = pstar0_f1710 * rprime;
313 } else if (strcmp(resonance, "f2010") == 0) {
314 spin = spin_f2010;
315 z0 = pstar0_f2010 * rprime;
316 } else if (strcmp(resonance, "chic0") == 0) {
317 spin = spin_chic0;
318 z0 = pstar0_chic0 * rprime;
319 } else {
320 printf("[BlattWeisskopf_pstarrprime] unsupported resonance\n");
321 exit(1);
322 }
323
324 if (spin == 0) return 1;
325 else if (spin == 1) return std::sqrt((1 + z0 * z0) / (1 + z * z));
326 else if (spin == 2) return std::sqrt((9 + 3 * z0 * z0 + z0 * z0 * z0 * z0) / (9 + 3 * z * z + z * z * z * z));
327 else {
328 printf("[BlattWeisskopf_pstarrprime] unsupported spin\n");
329 exit(1);
330 }
331
332 }
333
334 double EvtB0toK0K0K0::BlattWeisskopf_qr(double s13, double s23, const char* resonance)
335 {
336 // X_L(|q|r) when s12 and s23
337 int spin = -1;
338 double z0 = -1;
339 double z = Calculate_q_mag(s13, s23) * r;
340
341 if (strcmp(resonance, "f980") == 0) {
342 spin = spin_f980;
343 z0 = -1;
344 } else if (strcmp(resonance, "f1710") == 0) {
345 spin = spin_f1710;
346 z0 = q0_f1710 * r;
347 } else if (strcmp(resonance, "f2010") == 0) {
348 spin = spin_f2010;
349 z0 = q0_f2010 * r;
350 } else if (strcmp(resonance, "chic0") == 0) {
351 spin = spin_chic0;
352 z0 = q0_chic0 * r;
353 } else {
354 printf("[BlattWeisskopf_qr] unsupported resonance\n");
355 exit(1);
356 }
357
358 if (spin == 0) return 1;
359 else if (spin == 1) return std::sqrt((1 + z0 * z0) / (1 + z * z));
360 else if (spin == 2) return std::sqrt((9 + 3 * z0 * z0 + z0 * z0 * z0 * z0) / (9 + 3 * z * z + z * z * z * z));
361 else {
362 printf("[BlattWeisskopf_qr] unsupported spin\n");
363 exit(1);
364 }
365
366 }
367
368 double EvtB0toK0K0K0::Zemach(double s13, double s23, const char* resonance) // resonance: 1+2
369 {
370
371 int spin;
372 double p_dot_q = Calculate_q_dot_p_mag(s13, s23);
373 double p_mag = Calculate_p_mag(s13, s23);
374 double q_mag = Calculate_q_mag(s13, s23);
375
376 if (strcmp(resonance, "f980") == 0) {
377 spin = spin_f980;
378 } else if (strcmp(resonance, "f1710") == 0) {
379 spin = spin_f1710;
380 } else if (strcmp(resonance, "f2010") == 0) {
381 spin = spin_f2010;
382 } else if (strcmp(resonance, "chic0") == 0) {
383 spin = spin_chic0;
384 } else {
385 printf("[Zemach] unsupported resonance\n");
386 exit(1);
387 }
388
389 if (spin == 0) return 1.0;
390 else if (spin == 2) return (8.0 / 3.0) * (3 * p_dot_q * p_dot_q - p_mag * p_mag * q_mag * q_mag);
391 else {
392 printf("[Zemach] unsupported spin\n");
393 exit(1);
394 }
395
396 }
397
398 double EvtB0toK0K0K0::Calculate_m(double s13, double s23) // resonance: 1+2
399 {
400 return std::sqrt(mB0 * mB0 + mKS0 * mKS0 + mKL0 * mKL0 + mKL0 * mKL0 - s13 - s23);
401 }
402
403 double EvtB0toK0K0K0::Calculate_q_mag(double s13, double s23) // resonance: 1+2
404 {
405 double m = Calculate_m(s13, s23);
406 return std::sqrt(m * m / 4.0 - mKL0 * mKL0);
407 }
408
409 double EvtB0toK0K0K0::Calculate_pstar_mag(double s13, double s23) // resonance: 1+2
410 {
411 double m = Calculate_m(s13, s23);
412 return std::sqrt(std::pow(mB0 * mB0 - m * m - mKS0 * mKS0, 2) / 4.0 - m * m * mKS0 * mKS0) / mB0;
413 }
414
415 double EvtB0toK0K0K0::Calculate_p_mag(double s13, double s23) // resonance: 1+2
416 {
417 double m = Calculate_m(s13, s23);
418 double s12 = m * m;
419 return std::sqrt(std::pow(mB0 * mB0 - s12 - mKS0 * mKS0, 2) / (4 * s12) - mKS0 * mKS0);
420 }
421
422 double EvtB0toK0K0K0::Calculate_q_dot_p_mag(double s13, double s23) // resonance: 1+2
423 {
424 return std::abs((s13 - s23) / 4.0);
425 }
426
427 void EvtB0toK0K0K0::GetZeros()
428 {
429
430 // get q0 and pstar0 for f1710
431 q0_f1710 = std::sqrt((m0_f1710 * m0_f1710) / 4.0 - mKL0 * mKL0);
432 pstar0_f1710 = std::sqrt(std::pow(mB0 * mB0 - m0_f1710 * m0_f1710 - mKS0 * mKS0,
433 2) / 4.0 - m0_f1710 * m0_f1710 * mKS0 * mKS0) / mB0;
434
435 // get q0 and pstar0 for f2010
436 q0_f2010 = std::sqrt((m0_f2010 * m0_f2010) / 4.0 - mKL0 * mKL0);
437 pstar0_f2010 = std::sqrt(std::pow(mB0 * mB0 - m0_f2010 * m0_f2010 - mKS0 * mKS0,
438 2) / 4.0 - m0_f2010 * m0_f2010 * mKS0 * mKS0) / mB0;
439
440 // get q0 and pstar0 for chic0
441 q0_chic0 = std::sqrt((m0_chic0 * m0_chic0) / 4.0 - mKL0 * mKL0);
442 pstar0_chic0 = std::sqrt(std::pow(mB0 * mB0 - m0_chic0 * m0_chic0 - mKS0 * mKS0,
443 2) / 4.0 - m0_chic0 * m0_chic0 * mKS0 * mKS0) / mB0;
444
445 }
446
447 void EvtB0toK0K0K0::GetMasses()
448 {
449 mB0 = EvtPDL::getMass(EvtPDL::getId("B0"));
450 mKp = EvtPDL::getMass(EvtPDL::getId("K+"));
451 mKL0 = EvtPDL::getMass(EvtPDL::getId("K_L0"));
452 mKS0 = EvtPDL::getMass(EvtPDL::getId("K_S0"));
453 mK0 = EvtPDL::getMass(EvtPDL::getId("K0"));
454 mpic = EvtPDL::getMass(EvtPDL::getId("pi+"));
455 mpi0 = EvtPDL::getMass(EvtPDL::getId("pi0"));
456 }
457
459}
460
B2_EVTGEN_REGISTER_MODEL
#define B2_EVTGEN_REGISTER_MODEL(classname)
Class to register B2_EVTGEN_REGISTER_MODEL.
Definition: EvtGenModelRegister.h:75
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17