Belle II Software development
EvtBtoKK0K0.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/EvtBtoKK0K0.h"
23
24using std::endl;
25using namespace std::complex_literals;
26
27namespace Belle2 {
34 B2_EVTGEN_REGISTER_MODEL(EvtBtoKK0K0);
35
36 EvtBtoKK0K0::~EvtBtoKK0K0() {}
37
38 std::string EvtBtoKK0K0::getName()
39 {
40 return "BTOKK0K0";
41 }
42
43 EvtDecayBase* EvtBtoKK0K0::clone()
44 {
45 return new EvtBtoKK0K0;
46 }
47
48 void EvtBtoKK0K0::decay(EvtParticle* p)
49 {
50 // follow PhysRevD.85.112010
51
52 int Ntry = 0;
53 const int max_Ntry = 10000;
54
55 while (true) {
56 p->initializePhaseSpace(getNDaug(), getDaugs());
57
58 EvtParticle* ChargedKaon = p->getDaug(0);
59 EvtParticle* NeutralKaon_1 = p->getDaug(1);
60 EvtParticle* NeutralKaon_2 = p->getDaug(2);
61
62 EvtVector4R mom_ChargedKaon = ChargedKaon->getP4();
63 EvtVector4R mom_NeutralKaon_1 = NeutralKaon_1->getP4();
64 EvtVector4R mom_NeutralKaon_2 = NeutralKaon_2->getP4();
65
66 EvtVector4R mom_sum_1 = mom_ChargedKaon + mom_NeutralKaon_1;
67 EvtVector4R mom_sum_2 = mom_ChargedKaon + mom_NeutralKaon_2;
68
69 double s_1 = mom_sum_1.mass2();
70 double s_2 = mom_sum_2.mass2();
71
72 // follow the index of PhysRevD.85.112010
73 double s13 = -1;
74 double s23 = -1;
75 if (s_1 > s_2) {
76 s23 = s_1;
77 s13 = s_2;
78 } else {
79 s23 = s_2;
80 s13 = s_1;
81 }
82
83 // the maximum value of probability is obtained by brute fource method. scan all region and find.
84 const double Probability_max = 871.390583;
85 double Probability_value = Probability(s13, s23);
86 double Probability_ratio = Probability_value / Probability_max;
87
88 double xbox = EvtRandom::Flat(0.0, 1.0);
89
90 if (xbox < Probability_ratio) break;
91 Ntry++;
92 if (Ntry > max_Ntry) {
93 std::cout << "try to set the kinematics more than 10000 times. abort.\n";
94 ::abort();
95 }
96 }
97
98
99 return;
100
101 }
102
103
104 void EvtBtoKK0K0::initProbMax()
105 {
106 noProbMax();
107 }
108
109
110 void EvtBtoKK0K0::init()
111 {
112
113 // check that there are no arguments
114 checkNArg(0);
115
116 // there should be three daughters: K+ KL0 KL0
117 checkNDaug(3);
118
119 // Check that the daughters are correct
120 EvtId KaonPlusType = getDaug(0);
121 EvtId KaonZeroType_1 = getDaug(1);
122 EvtId KaonZeroType_2 = getDaug(2);
123
124 int KaonPlustyp = 0;
125 int KaonZerotyp = 0;
126
127 if ((KaonPlusType == EvtPDL::getId("K+")) || (KaonPlusType == EvtPDL::getId("K-"))) KaonPlustyp++;
128 if ((KaonZeroType_1 == EvtPDL::getId("K0")) ||
129 (KaonZeroType_1 == EvtPDL::getId("anti-K0")) ||
130 (KaonZeroType_1 == EvtPDL::getId("K_S0")) ||
131 (KaonZeroType_1 == EvtPDL::getId("K_L0"))) KaonZerotyp++;
132 if ((KaonZeroType_2 == EvtPDL::getId("K0")) ||
133 (KaonZeroType_2 == EvtPDL::getId("anti-K0")) ||
134 (KaonZeroType_2 == EvtPDL::getId("K_S0")) ||
135 (KaonZeroType_2 == EvtPDL::getId("K_L0"))) KaonZerotyp++;
136
137 if ((KaonPlustyp != 1) || (KaonZerotyp != 2)) {
138
139 std::cout << "The first dauther should be charged Kaon. The second and third daughters should be K0, anti-K0, K_L0, or K_S0.\n";
140 ::abort();
141 }
142
143 // initialize q0 and pstar0
144 GetMasses();
145 GetZeros();
146
147 }
148
149 double EvtBtoKK0K0::Probability(double s13, double s23)
150 {
151 std::complex<double> total_amplitude = Amplitude(s13, s23, "f980", false) + Amplitude(s13, s23, "f1500", false) + Amplitude(s13,
152 s23, "f1525", false) + Amplitude(s13, s23, "f1710", false) + Amplitude(s13, s23, "chic0", false) + Amplitude(s13, s23, "NR", false);
153 std::complex<double> total_amplitude_isobar = Amplitude(s13, s23, "f980", true) + Amplitude(s13, s23, "f1500",
154 true) + Amplitude(s13, s23, "f1525", true) + Amplitude(s13, s23, "f1710", true) + Amplitude(s13, s23, "chic0",
155 true) + Amplitude(s13, s23, "NR", true);
156 return std::abs(total_amplitude) * std::abs(total_amplitude) + std::abs(total_amplitude_isobar) * std::abs(total_amplitude_isobar);
157 }
158
159
160 std::complex<double> EvtBtoKK0K0::Amplitude(double s13, double s23, const char* resonance, bool isobar = false)
161 {
162
163 if (strcmp(resonance, "f980") == 0) {
164 std::complex<double> a;
165 if (isobar == false) a = c_f980 * std::exp(1i * DegreeToRadian(phi_f980));
166 else a = c_f980 * std::exp(1i * DegreeToRadian(phi_f980));
167 return a * (DynamicalAmplitude(s13, s23, resonance) + DynamicalAmplitude(s23, s13, resonance));
168 } else if (strcmp(resonance, "f1500") == 0) {
169 std::complex<double> a;
170 if (isobar == false) a = c_f1500 * std::exp(1i * DegreeToRadian(phi_f1500));
171 else a = c_f1500 * std::exp(1i * DegreeToRadian(phi_f1500));
172 return a * (DynamicalAmplitude(s13, s23, resonance) + DynamicalAmplitude(s23, s13, resonance));
173 } else if (strcmp(resonance, "f1525") == 0) {
174 std::complex<double> a;
175 if (isobar == false) a = c_f1525 * std::exp(1i * DegreeToRadian(phi_f1525));
176 else a = c_f1525 * std::exp(1i * DegreeToRadian(phi_f1525));
177 return a * (DynamicalAmplitude(s13, s23, resonance) + DynamicalAmplitude(s23, s13, resonance));
178 } else if (strcmp(resonance, "f1710") == 0) {
179 std::complex<double> a;
180 if (isobar == false) a = c_f1710 * std::exp(1i * DegreeToRadian(phi_f1710));
181 else a = c_f1710 * std::exp(1i * DegreeToRadian(phi_f1710));
182 return a * (DynamicalAmplitude(s13, s23, resonance) + DynamicalAmplitude(s23, s13, resonance));
183 } else if (strcmp(resonance, "chic0") == 0) {
184 std::complex<double> a;
185 if (isobar == false) a = c_chic0 * std::exp(1i * (DegreeToRadian(phi_chic0) + DegreeToRadian(delta_chic0)));
186 else a = c_chic0 * std::exp(1i * (DegreeToRadian(phi_chic0) - DegreeToRadian(delta_chic0)));
187 return a * (DynamicalAmplitude(s13, s23, resonance) + DynamicalAmplitude(s23, s13, resonance));
188 } else if (strcmp(resonance, "NR") == 0) {
189
190 double Omega = 0.5 * (mB + (1.0 / 3.0) * (mKp + mKL0 + mKL0));
191 double m12 = std::sqrt(mB * mB + mKp * mKp + mKL0 * mKL0 + mKL0 * mKL0 - s13 - s23); // sqrt(s12)
192 double x = m12 - Omega;
193
194 std::complex<double> aS0;
195 std::complex<double> aS1;
196 std::complex<double> aS2;
197
198 if (isobar == false) {
199 aS0 = aS0_c_NR * (1 + b_NR) * std::exp(1i * DegreeToRadian(aS0_phi_NR));
200 aS1 = aS1_c_NR * (1 + b_NR) * std::exp(1i * DegreeToRadian(aS1_phi_NR));
201 aS2 = aS2_c_NR * (1 + b_NR) * std::exp(1i * DegreeToRadian(aS2_phi_NR));
202 } else {
203 aS0 = aS0_c_NR * (1 - b_NR) * std::exp(1i * DegreeToRadian(aS0_phi_NR));
204 aS1 = aS1_c_NR * (1 - b_NR) * std::exp(1i * DegreeToRadian(aS1_phi_NR));
205 aS2 = aS2_c_NR * (1 - b_NR) * std::exp(1i * DegreeToRadian(aS2_phi_NR));
206 }
207
208 return 2.0 * (aS0 + aS1 * x + aS2 * x * x);
209 } else {
210 printf("[Amplitude] unsupported resonance\n");
211 ::abort();
212 }
213
214 }
215
216 std::complex<double> EvtBtoKK0K0::DynamicalAmplitude(double s13, double s23, const char* resonance)
217 {
218
219 if (strcmp(resonance, "f980") == 0) {
220 return Flatte(s13, s23, resonance) * BlattWeisskopf_pstarrprime(s13, s23, resonance) * BlattWeisskopf_qr(s13, s23,
221 resonance) * Zemach(s13, s23, resonance);
222 } else if (strcmp(resonance, "f1500") == 0) {
223 return RBW(s13, s23, resonance) * BlattWeisskopf_pstarrprime(s13, s23, resonance) * BlattWeisskopf_qr(s13, s23,
224 resonance) * Zemach(s13, s23, resonance);
225 } else if (strcmp(resonance, "f1525") == 0) {
226 return RBW(s13, s23, resonance) * BlattWeisskopf_pstarrprime(s13, s23, resonance) * BlattWeisskopf_qr(s13, s23,
227 resonance) * Zemach(s13, s23, resonance);
228 } else if (strcmp(resonance, "f1710") == 0) {
229 return RBW(s13, s23, resonance) * BlattWeisskopf_pstarrprime(s13, s23, resonance) * BlattWeisskopf_qr(s13, s23,
230 resonance) * Zemach(s13, s23, resonance);
231 } else if (strcmp(resonance, "chic0") == 0) {
232 return RBW(s13, s23, resonance) * BlattWeisskopf_pstarrprime(s13, s23, resonance) * BlattWeisskopf_qr(s13, s23,
233 resonance) * Zemach(s13, s23, resonance);
234 } else {
235 printf("[Amplitude] unsupported resonance\n");
236 ::abort();
237 }
238
239 }
240
241 std::complex<double> EvtBtoKK0K0::Flatte(double s13, double s23, const char* resonance)
242 {
243 double m0;
244 double gpi;
245 double gK;
246
247 if (strcmp(resonance, "f980") == 0) {
248 m0 = m0_f980;
249 gpi = gpi_f980;
250 gK = gKgpi_f980 * gpi_f980;
251 } else {
252 printf("[Flatte] unsupported resonance\n");
253 ::abort();
254 }
255
256 double m = Calculate_m(s13, s23);
257
258 double rho_pipi = std::sqrt(1 - 4 * mpic * mpic / (m * m));
259 double rho_KK = std::sqrt(1 - 4 * mK * mK / (m * m));
260
261 return 1.0 / ((m0 * m0 - m * m) - 1i * (gpi * rho_pipi + gK * rho_KK));
262 }
263
264 std::complex<double> EvtBtoKK0K0::RBW(double s13, double s23, const char* resonance)
265 {
266 double m0;
267
268 double m = Calculate_m(s13, s23);
269
270 if (strcmp(resonance, "f1500") == 0) {
271 m0 = m0_f1500;
272 } else if (strcmp(resonance, "f1525") == 0) {
273 m0 = m0_f1525;
274 } else if (strcmp(resonance, "f1710") == 0) {
275 m0 = m0_f1710;
276 } else if (strcmp(resonance, "chic0") == 0) {
277 m0 = m0_chic0;
278 } else {
279 printf("[RBW] unsupported resonance\n");
280 ::abort();
281 }
282
283 return 1.0 / (m0 * m0 - m * m - 1i * m0 * MassDepWidth(s13, s23, resonance));
284 }
285
286 double EvtBtoKK0K0::MassDepWidth(double s13, double s23, const char* resonance)
287 {
288 // Gamma(m) when s12 and s23
289 int spin = -1;
290 double m0;
291 double q0;
292 double Gamma0;
293
294 double q_mag = Calculate_q_mag(s13, s23);
295 double m = Calculate_m(s13, s23);
296
297 if (strcmp(resonance, "f1500") == 0) {
298 spin = spin_f1500;
299 m0 = m0_f1500;
300 q0 = q0_f1500;
301 Gamma0 = Gamma0_f1500;
302 } else if (strcmp(resonance, "f1525") == 0) {
303 spin = spin_f1525;
304 m0 = m0_f1525;
305 q0 = q0_f1525;
306 Gamma0 = Gamma0_f1525;
307 } else if (strcmp(resonance, "f1710") == 0) {
308 spin = spin_f1710;
309 m0 = m0_f1710;
310 q0 = q0_f1710;
311 Gamma0 = Gamma0_f1710;
312 } else if (strcmp(resonance, "chic0") == 0) {
313 spin = spin_chic0;
314 m0 = m0_chic0;
315 q0 = q0_chic0;
316 Gamma0 = Gamma0_chic0;
317 } else {
318 printf("[MassDepWidth] unsupported resonance\n");
319 ::abort();
320 }
321
322 return Gamma0 * std::pow(q_mag / q0, 2 * spin + 1) * (m0 / m) * std::pow(BlattWeisskopf_qr(s13, s23, resonance), 2);
323 }
324
325 double EvtBtoKK0K0::BlattWeisskopf_pstarrprime(double s13, double s23, const char* resonance)
326 {
327 // X_L(|p*|r') when s12 and s23
328 int spin = -1;
329 double z0 = -1;
330 double z = Calculate_pstar_mag(s13, s23) * rprime;
331
332 if (strcmp(resonance, "f980") == 0) {
333 spin = spin_f980;
334 z0 = -1;
335 } else if (strcmp(resonance, "f1500") == 0) {
336 spin = spin_f1500;
337 z0 = pstar0_f1500 * rprime;
338 } else if (strcmp(resonance, "f1525") == 0) {
339 spin = spin_f1525;
340 z0 = pstar0_f1525 * rprime;
341 } else if (strcmp(resonance, "f1710") == 0) {
342 spin = spin_f1710;
343 z0 = pstar0_f1710 * rprime;
344 } else if (strcmp(resonance, "chic0") == 0) {
345 spin = spin_chic0;
346 z0 = pstar0_chic0 * rprime;
347 } else {
348 printf("[BlattWeisskopf_pstarrprime] unsupported resonance\n");
349 ::abort();
350 }
351
352 if (spin == 0) return 1;
353 else if (spin == 1) return std::sqrt((1 + z0 * z0) / (1 + z * z));
354 else if (spin == 2) return std::sqrt((9 + 3 * z0 * z0 + z0 * z0 * z0 * z0) / (9 + 3 * z * z + z * z * z * z));
355 else {
356 printf("[BlattWeisskopf_pstarrprime] unsupported spin\n");
357 ::abort();
358 }
359
360 }
361
362 double EvtBtoKK0K0::BlattWeisskopf_qr(double s13, double s23, const char* resonance)
363 {
364 // X_L(|q|r) when s12 and s23
365 int spin = -1;
366 double z0 = -1;
367 double z = Calculate_q_mag(s13, s23) * r;
368
369 if (strcmp(resonance, "f980") == 0) {
370 spin = spin_f980;
371 z0 = -1;
372 } else if (strcmp(resonance, "f1500") == 0) {
373 spin = spin_f1500;
374 z0 = q0_f1500 * r;
375 } else if (strcmp(resonance, "f1525") == 0) {
376 spin = spin_f1525;
377 z0 = q0_f1525 * r;
378 } else if (strcmp(resonance, "f1710") == 0) {
379 spin = spin_f1710;
380 z0 = q0_f1710 * r;
381 } else if (strcmp(resonance, "chic0") == 0) {
382 spin = spin_chic0;
383 z0 = q0_chic0 * r;
384 } else {
385 printf("[BlattWeisskopf_qr] unsupported resonance\n");
386 ::abort();
387 }
388
389 if (spin == 0) return 1;
390 else if (spin == 1) return std::sqrt((1 + z0 * z0) / (1 + z * z));
391 else if (spin == 2) return std::sqrt((9 + 3 * z0 * z0 + z0 * z0 * z0 * z0) / (9 + 3 * z * z + z * z * z * z));
392 else {
393 printf("[BlattWeisskopf_qr] unsupported spin\n");
394 ::abort();
395 }
396
397 }
398
399 double EvtBtoKK0K0::Zemach(double s13, double s23, const char* resonance)
400 {
401
402 int spin;
403 double p_dot_q = Calculate_q_dot_p_mag(s13, s23);
404 double p_mag = Calculate_p_mag(s13, s23);
405 double q_mag = Calculate_q_mag(s13, s23);
406
407 if (strcmp(resonance, "f980") == 0) {
408 spin = spin_f980;
409 } else if (strcmp(resonance, "f1500") == 0) {
410 spin = spin_f1500;
411 } else if (strcmp(resonance, "f1525") == 0) {
412 spin = spin_f1525;
413 } else if (strcmp(resonance, "f1710") == 0) {
414 spin = spin_f1710;
415 } else if (strcmp(resonance, "chic0") == 0) {
416 spin = spin_chic0;
417 } else {
418 printf("[Zemach] unsupported resonance\n");
419 ::abort();
420 }
421
422 if (spin == 0) return 1.0;
423 else if (spin == 1) return 4 * p_dot_q;
424 else if (spin == 2) return (16.0 / 3.0) * (3 * p_dot_q * p_dot_q - p_mag * p_mag * q_mag * q_mag);
425 else {
426 printf("[Zemach] unsupported spin\n");
427 ::abort();
428 }
429
430 }
431
432 double EvtBtoKK0K0::Calculate_m(double s13, double s23)
433 {
434 return std::sqrt(mB * mB + mKp * mKp + mKL0 * mKL0 + mKL0 * mKL0 - s13 - s23);
435 }
436
437 double EvtBtoKK0K0::Calculate_q_mag(double s13, double s23)
438 {
439 double m = Calculate_m(s13, s23);
440 return std::sqrt(m * m / 4.0 - mKL0 * mKL0);
441 }
442
443 double EvtBtoKK0K0::Calculate_pstar_mag(double s13, double s23)
444 {
445 double m = Calculate_m(s13, s23);
446 return std::sqrt(std::pow(mB * mB - m * m - mKp * mKp, 2) / 4.0 - m * m * mKp * mKp) / mB;
447 }
448
449 double EvtBtoKK0K0::Calculate_p_mag(double s13, double s23)
450 {
451 double m = Calculate_m(s13, s23);
452 double s12 = m * m;
453 return std::sqrt(std::pow(mB * mB - s12 - mKp * mKp, 2) / (4 * s12) - mKp * mKp);
454 }
455
456 double EvtBtoKK0K0::Calculate_q_dot_p_mag(double s13, double s23)
457 {
458 return std::abs((s13 - s23) / 4.0);
459 }
460
461 double EvtBtoKK0K0::DegreeToRadian(double degree)
462 {
463 return (3.141592 * degree) / 180.0;
464 }
465
466 void EvtBtoKK0K0::GetZeros()
467 {
468
469 // get q0 and pstar0 for f1500
470 q0_f1500 = std::sqrt((m0_f1500 * m0_f1500) / 4.0 - mKL0 * mKL0);
471 pstar0_f1500 = std::sqrt(std::pow(mB * mB - m0_f1500 * m0_f1500 - mKp * mKp, 2) / 4.0 - m0_f1500 * m0_f1500 * mKp * mKp) / mB;
472
473 // get q0 and pstar0 for f1525
474 q0_f1525 = std::sqrt((m0_f1525 * m0_f1525) / 4.0 - mKL0 * mKL0);
475 pstar0_f1525 = std::sqrt(std::pow(mB * mB - m0_f1525 * m0_f1525 - mKp * mKp, 2) / 4.0 - m0_f1525 * m0_f1525 * mKp * mKp) / mB;
476
477 // get q0 and pstar0 for f1710
478 q0_f1710 = std::sqrt((m0_f1710 * m0_f1710) / 4.0 - mKL0 * mKL0);
479 pstar0_f1710 = std::sqrt(std::pow(mB * mB - m0_f1710 * m0_f1710 - mKp * mKp, 2) / 4.0 - m0_f1710 * m0_f1710 * mKp * mKp) / mB;
480
481 // get q0 and pstar0 for chic0
482 q0_chic0 = std::sqrt((m0_chic0 * m0_chic0) / 4.0 - mKL0 * mKL0);
483 pstar0_chic0 = std::sqrt(std::pow(mB * mB - m0_chic0 * m0_chic0 - mKp * mKp, 2) / 4.0 - m0_chic0 * m0_chic0 * mKp * mKp) / mB;
484
485 }
486
487 void EvtBtoKK0K0::GetMasses()
488 {
489 mB = EvtPDL::getMass(EvtPDL::getId("B+"));
490 mKp = EvtPDL::getMass(EvtPDL::getId("K+"));
491 mKL0 = EvtPDL::getMass(EvtPDL::getId("K_L0"));
492 mK = (EvtPDL::getMass(EvtPDL::getId("K+")) + EvtPDL::getMass(EvtPDL::getId("K0"))) / 2.0;
493 mpic = EvtPDL::getMass(EvtPDL::getId("pi+"));
494 }
495
497}
498
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