Belle II Software development
EvtB0toKsKK.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 <iomanip>
10#include <cmath>
11#include <string>
12
13#include "EvtGenBase/EvtCPUtil.hh"
14#include "EvtGenBase/EvtPDL.hh"
15#include "EvtGenBase/EvtTensor4C.hh"
16
17#include <generators/evtgen/EvtGenModelRegister.h>
18#include "generators/evtgen/models/EvtB0toKsKK.h"
19
20namespace Belle2 {
28
29 EvtB0toKsKK::~EvtB0toKsKK() {}
30
32 {
33 return "B0toKsKK";
34 }
35
36 EvtDecayBase* EvtB0toKsKK::clone()
37 {
38 return new EvtB0toKsKK;
39 }
40
42 {
43 // Check number of arguments
44 checkNArg(32);
45
46 // Check number of daughters
47 checkNDaug(3);
48
49 // Check decay chain
50 if ((getParentId() != EvtPDL::getId("B0") &&
51 getParentId() != EvtPDL::getId("anti-B0")) ||
52 (getDaug(0) != EvtPDL::getId("K_S0")) ||
53 (getDaug(1) != EvtPDL::getId("K+") &&
54 getDaug(1) != EvtPDL::getId("K-")) ||
55 (getDaug(2) != EvtPDL::getId("K+") &&
56 getDaug(2) != EvtPDL::getId("K-"))) {
57 std::cout << "ERROR: Invalid decay" << std::endl;
58 std::cout << "USAGE: K_S0 K+ K-" << std::endl;
59 exit(1);
60 }
61
62 a_f0ks_ =
63 (1.0 + getArg(3)) * EvtComplex(getArg(1) * cos(getArg(2) * M_PI / 180.0),
64 getArg(1) * sin(getArg(2) * M_PI / 180.0));
65 a_phiks_ =
66 (1.0 + getArg(7)) * EvtComplex(getArg(5) * cos(getArg(6) * M_PI / 180.0),
67 getArg(5) * sin(getArg(6) * M_PI / 180.0));
68 a_fxks_ =
69 (1.0 + getArg(11)) * EvtComplex(getArg(9) * cos(getArg(10) * M_PI / 180.0),
70 getArg(9) * sin(getArg(10) * M_PI / 180.0));
72 (1.0 + getArg(15)) * EvtComplex(getArg(13) * cos(getArg(14) * M_PI / 180.0),
73 getArg(13) * sin(getArg(14) * M_PI / 180.0));
74 a_kpkmnr_ =
75 (1.0 + getArg(19)) * EvtComplex(getArg(17) * cos(getArg(18) * M_PI / 180.0),
76 getArg(17) * sin(getArg(18) * M_PI / 180.0));
77 a_kskpnr_ =
78 (1.0 + getArg(24)) * EvtComplex(getArg(22) * cos(getArg(23) * M_PI / 180.0),
79 getArg(22) * sin(getArg(23) * M_PI / 180.0));
80 a_kskmnr_ =
81 (1.0 + getArg(29)) * EvtComplex(getArg(27) * cos(getArg(28) * M_PI / 180.0),
82 getArg(27) * sin(getArg(28) * M_PI / 180.0));
83
85 (1.0 - getArg(3)) * EvtComplex(getArg(1) * cos(getArg(2) * M_PI / 180.0),
86 getArg(1) * sin(getArg(2) * M_PI / 180.0));
88 (1.0 - getArg(7)) * EvtComplex(getArg(5) * cos(getArg(6) * M_PI / 180.0),
89 getArg(5) * sin(getArg(6) * M_PI / 180.0));
91 (1.0 - getArg(11)) * EvtComplex(getArg(9) * cos(getArg(10) * M_PI / 180.0),
92 getArg(9) * sin(getArg(10) * M_PI / 180.0));
94 (1.0 - getArg(15)) * EvtComplex(getArg(13) * cos(getArg(14) * M_PI / 180.0),
95 getArg(13) * sin(getArg(14) * M_PI / 180.0));
97 (1.0 - getArg(19)) * EvtComplex(getArg(17) * cos(getArg(18) * M_PI / 180.0),
98 getArg(17) * sin(getArg(18) * M_PI / 180.0));
100 (1.0 - getArg(24)) * EvtComplex(getArg(22) * cos(getArg(23) * M_PI / 180.0),
101 getArg(22) * sin(getArg(23) * M_PI / 180.0));
103 (1.0 - getArg(29)) * EvtComplex(getArg(27) * cos(getArg(28) * M_PI / 180.0),
104 getArg(27) * sin(getArg(28) * M_PI / 180.0));
105
106 alpha_kpkmnr = getArg(21);
107 alpha_kskpnr = getArg(26);
108 alpha_kskmnr = getArg(31);
109
110 std::cout << setiosflags(std::ios::left) << std::endl
111 << "B0 Channel " << std::setw(20) << "Relative amplitude" << std::setw(20) << "Relative phase" << std::endl
112 << "f0ks " << std::setw(20) << real(a_f0ks_) << std::setw(20) << imag(a_f0ks_) << std::endl
113 << "phiks " << std::setw(20) << real(a_phiks_) << std::setw(20) << imag(a_phiks_) << std::endl
114 << "fxks " << std::setw(20) << real(a_fxks_) << std::setw(20) << imag(a_fxks_) << std::endl
115 << "chic0ks " << std::setw(20) << real(a_chic0ks_) << std::setw(20) << imag(a_chic0ks_) << std::endl
116 << "kpkmnr " << std::setw(20) << real(a_kpkmnr_) << std::setw(20) << imag(a_kpkmnr_) << std::endl
117 << "kskpnr " << std::setw(20) << real(a_kskpnr_) << std::setw(20) << imag(a_kskpnr_) << std::endl
118 << "kskmnr " << std::setw(20) << real(a_kskmnr_) << std::setw(20) << imag(a_kskmnr_) << std::endl
119 << std::endl;
120
121 std::cout << setiosflags(std::ios::left) << std::endl
122 << "B0B Channel " << std::setw(20) << "Relative amplitude" << std::setw(20) << "Relative phase" << std::endl
123 << "f0ks " << std::setw(20) << real(abar_f0ks_) << std::setw(20) << imag(abar_f0ks_) << std::endl
124 << "phiks " << std::setw(20) << real(abar_phiks_) << std::setw(20) << imag(abar_phiks_) << std::endl
125 << "fxks " << std::setw(20) << real(abar_fxks_) << std::setw(20) << imag(abar_fxks_) << std::endl
126 << "chic0ks " << std::setw(20) << real(abar_chic0ks_) << std::setw(20) << imag(abar_chic0ks_) << std::endl
127 << "kpkmnr " << std::setw(20) << real(abar_kpkmnr_) << std::setw(20) << imag(abar_kpkmnr_) << std::endl
128 << "kskpnr " << std::setw(20) << real(abar_kskpnr_) << std::setw(20) << imag(abar_kskpnr_) << std::endl
129 << "kskmnr " << std::setw(20) << real(abar_kskmnr_) << std::setw(20) << imag(abar_kskmnr_) << std::endl
130 << std::endl;
131
132 // Open debugging file
133 debugfile_.open("debug.dat", std::ios::out);
134 }
135
137 {
138 //setProbMax(50000.0);
139 //setProbMax(100000.0);
140 //setProbMax(200000.0);
141 //setProbMax(400000.0);
142 //setProbMax(600000.0);
143 setProbMax(1100000.0);
144 }
145
146 void EvtB0toKsKK::decay(EvtParticle* p)
147 {
148 // Btag
149 static EvtId B0 = EvtPDL::getId("B0");
150 static EvtId B0B = EvtPDL::getId("anti-B0");
151
152 double t;
153 EvtId other_b;
154
155 //std::cout << EvtCPUtil::getInstance()->getMixingType() << std::endl;
156 //EvtCPUtil::getInstance()->setMixingType(0);
157 EvtCPUtil::getInstance()->OtherB(p, t, other_b, 0.5);
158 //EvtCPUtil::getInstance()->OtherB(p,t,other_b,0.4);
159 //EvtCPUtil::getInstance()->OtherB(p,t,other_b);
160
161 // Brec
162 p->initializePhaseSpace(getNDaug(), getDaugs());
163 EvtVector4R p4ks, p4kp, p4km;
164 if (p->getId() == B0) {
165 p4ks = p->getDaug(0)->getP4();
166 p4kp = p->getDaug(1)->getP4();
167 p4km = p->getDaug(2)->getP4();
170 /*p4ks = p->getDaug(0)->getP4();
171 p4kp = p->getDaug(2)->getP4();
172 p4km = p->getDaug(1)->getP4();*/
173 } else {
174 p4ks = p->getDaug(0)->getP4();
175 p4kp = p->getDaug(2)->getP4();
176 p4km = p->getDaug(1)->getP4();
178 /*p4ks = p->getDaug(0)->getP4();
179 p4kp = p->getDaug(1)->getP4();
180 p4km = p->getDaug(2)->getP4();*/
181 }
182
183 /*std::cout << (p4ks + p4kp).mass() << " " << (p4ks + p4km).mass() << " "
184 << (p4kp + p4km).mass() << std::endl;
185 if( p->getId() == other_b )
186 std::cout << "same flavour" << std::endl;
187 std::cout << p->getId() << " --> " << p->getDaug(0)->getId() << " "
188 << p->getDaug(1)->getId() << " " << p->getDaug(2)->getId()
189 << std::endl;
190 std::cout << B0 << std::endl;
191 std::cout << B0B << std::endl;
192 std::cout << other_b << std::endl;*/
193 /*std::cout << p->getLifetime() << std::endl;
194 std::cout << p->get4Pos() << std::endl;
195 std::cout << p->getP4() << std::endl;
196 std::cout << p->getP4Lab() << std::endl;
197 EvtParticle *parent=p->getParent();
198 if (parent->getDaug(0)!=p){
199 std::cout << parent->getDaug(0)->getLifetime() << std::endl;
200 std::cout << parent->getDaug(0)->get4Pos() << std::endl;
201 std::cout << parent->getDaug(0)->getP4() << std::endl;
202 std::cout << parent->getDaug(0)->getP4Lab() << std::endl;
203 std::cout << parent->getDaug(0)->getP4Lab()+p->getP4Lab() << std::endl;
204 }
205 else{
206 std::cout << parent->getDaug(1)->getLifetime() << std::endl;
207 std::cout << parent->getDaug(1)->get4Pos() << std::endl;
208 std::cout << parent->getDaug(1)->getP4() << std::endl;
209 std::cout << parent->getDaug(1)->getP4Lab() << std::endl;
210 std::cout << parent->getDaug(1)->getP4Lab()+p->getP4Lab() << std::endl;
211 }
212 std::cout << parent->getP4() << std::endl;
213 std::cout << t << std::endl;*/
214
215 // Relative amplides and phases with direct CP violation
216 a_f0ks_ =
217 (1.0 + getArg(3)) * EvtComplex(getArg(1) * cos(getArg(2) * M_PI / 180.0),
218 getArg(1) * sin(getArg(2) * M_PI / 180.0));
219 a_phiks_ =
220 (1.0 + getArg(7)) * EvtComplex(getArg(5) * cos(getArg(6) * M_PI / 180.0),
221 getArg(5) * sin(getArg(6) * M_PI / 180.0));
222 a_fxks_ =
223 (1.0 + getArg(11)) * EvtComplex(getArg(9) * cos(getArg(10) * M_PI / 180.0),
224 getArg(9) * sin(getArg(10) * M_PI / 180.0));
225 a_chic0ks_ =
226 (1.0 + getArg(15)) * EvtComplex(getArg(13) * cos(getArg(14) * M_PI / 180.0),
227 getArg(13) * sin(getArg(14) * M_PI / 180.0));
228 a_kpkmnr_ =
229 (1.0 + getArg(19)) * EvtComplex(getArg(17) * cos(getArg(18) * M_PI / 180.0),
230 getArg(17) * sin(getArg(18) * M_PI / 180.0));
231 a_kskpnr_ =
232 (1.0 + getArg(24)) * EvtComplex(getArg(22) * cos(getArg(23) * M_PI / 180.0),
233 getArg(22) * sin(getArg(23) * M_PI / 180.0));
234 a_kskmnr_ =
235 (1.0 + getArg(29)) * EvtComplex(getArg(27) * cos(getArg(28) * M_PI / 180.0),
236 getArg(27) * sin(getArg(28) * M_PI / 180.0));
237
238 abar_f0ks_ =
239 (1.0 - getArg(3)) * EvtComplex(getArg(1) * cos(getArg(2) * M_PI / 180.0),
240 getArg(1) * sin(getArg(2) * M_PI / 180.0));
242 (1.0 - getArg(7)) * EvtComplex(getArg(5) * cos(getArg(6) * M_PI / 180.0),
243 getArg(5) * sin(getArg(6) * M_PI / 180.0));
244 abar_fxks_ =
245 (1.0 - getArg(11)) * EvtComplex(getArg(9) * cos(getArg(10) * M_PI / 180.0),
246 getArg(9) * sin(getArg(10) * M_PI / 180.0));
248 (1.0 - getArg(15)) * EvtComplex(getArg(13) * cos(getArg(14) * M_PI / 180.0),
249 getArg(13) * sin(getArg(14) * M_PI / 180.0));
251 (1.0 - getArg(19)) * EvtComplex(getArg(17) * cos(getArg(18) * M_PI / 180.0),
252 getArg(17) * sin(getArg(18) * M_PI / 180.0));
254 (1.0 - getArg(24)) * EvtComplex(getArg(22) * cos(getArg(23) * M_PI / 180.0),
255 getArg(22) * sin(getArg(23) * M_PI / 180.0));
257 (1.0 - getArg(29)) * EvtComplex(getArg(27) * cos(getArg(28) * M_PI / 180.0),
258 getArg(27) * sin(getArg(28) * M_PI / 180.0));
259
260 // Mixing-induced CP asymmetry
261 const double pCP_f0ks = getArg(4) * M_PI / 180.0;
262 const double pCP_phiks = getArg(8) * M_PI / 180.0;
263 const double pCP_fxks = getArg(12) * M_PI / 180.0;
264 const double pCP_chic0ks = getArg(16) * M_PI / 180.0;
265 const double pCP_kpkmnr = getArg(20) * M_PI / 180.0;
266 const double pCP_kskpnr = getArg(25) * M_PI / 180.0;
267 const double pCP_kskmnr = getArg(30) * M_PI / 180.0;
268
269 if (other_b == B0) {
270 a_f0ks_ *=
271 EvtComplex(cos(+2.0 * pCP_f0ks), sin(+2.0 * pCP_f0ks));
272 a_phiks_ *=
273 EvtComplex(cos(+2.0 * pCP_phiks), sin(+2.0 * pCP_phiks));
274 a_fxks_ *=
275 EvtComplex(cos(+2.0 * pCP_fxks), sin(+2.0 * pCP_fxks));
276 a_chic0ks_ *=
277 EvtComplex(cos(+2.0 * pCP_chic0ks), sin(+2.0 * pCP_chic0ks));
278 a_kpkmnr_ *=
279 EvtComplex(cos(+2.0 * pCP_kpkmnr), sin(+2.0 * pCP_kpkmnr));
280 a_kskpnr_ *=
281 EvtComplex(cos(+2.0 * pCP_kskpnr), sin(+2.0 * pCP_kskpnr));
282 a_kskmnr_ *=
283 EvtComplex(cos(+2.0 * pCP_kskmnr), sin(+2.0 * pCP_kskmnr));
284 }
285 if (other_b == B0B) {
286 abar_f0ks_ *=
287 EvtComplex(cos(-2.0 * pCP_f0ks), sin(-2.0 * pCP_f0ks));
288 abar_phiks_ *=
289 EvtComplex(cos(-2.0 * pCP_phiks), sin(-2.0 * pCP_phiks));
290 abar_fxks_ *=
291 EvtComplex(cos(-2.0 * pCP_fxks), sin(-2.0 * pCP_fxks));
293 EvtComplex(cos(-2.0 * pCP_chic0ks), sin(-2.0 * pCP_chic0ks));
294 abar_kpkmnr_ *=
295 EvtComplex(cos(-2.0 * pCP_kpkmnr), sin(-2.0 * pCP_kpkmnr));
296 abar_kskpnr_ *=
297 EvtComplex(cos(-2.0 * pCP_kskpnr), sin(-2.0 * pCP_kskpnr));
298 abar_kskmnr_ *=
299 EvtComplex(cos(-2.0 * pCP_kskmnr), sin(-2.0 * pCP_kskmnr));
300 }
301
302 //Form Factors
303 EvtComplex Amp_f0ks = A_f0ks(p4ks, p4kp, p4km);
304 EvtComplex Amp_phiks = A_phiks(p4ks, p4kp, p4km);
305 EvtComplex Amp_fxks = A_fxks(p4ks, p4kp, p4km);
306 EvtComplex Amp_chic0ks = A_chic0ks(p4ks, p4kp, p4km);
307 EvtComplex Amp_kpkmnr = A_kknr(p4kp, p4km, alpha_kpkmnr);
308 EvtComplex Amp_kskpnr = A_kknr(p4ks, p4kp, alpha_kskpnr);
309 EvtComplex Amp_kskmnr = A_kknr(p4ks, p4km, alpha_kskmnr);
310
311 EvtComplex Ampbar_f0ks = A_f0ks(p4ks, p4km, p4kp);
312 EvtComplex Ampbar_phiks = A_phiks(p4ks, p4km, p4kp);
313 EvtComplex Ampbar_fxks = A_fxks(p4ks, p4km, p4kp);
314 EvtComplex Ampbar_chic0ks = A_chic0ks(p4ks, p4km, p4kp);
315 EvtComplex Ampbar_kpkmnr = A_kknr(p4km, p4kp, alpha_kpkmnr);
316 EvtComplex Ampbar_kskpnr = A_kknr(p4ks, p4km, alpha_kskpnr);
317 EvtComplex Ampbar_kskmnr = A_kknr(p4ks, p4kp, alpha_kskmnr);
318
319 const EvtComplex A_B0toKsKK =
320 (a_f0ks_ * Amp_f0ks) +
321 (a_phiks_ * Amp_phiks) +
322 (a_fxks_ * Amp_fxks) +
323 (a_chic0ks_ * Amp_chic0ks) +
324 (a_kpkmnr_ * Amp_kpkmnr) +
325 (a_kskpnr_ * Amp_kskpnr) +
326 (a_kskmnr_ * Amp_kskmnr);
327
328 const EvtComplex Abar_B0toKsKK =
329 (abar_f0ks_ * Ampbar_f0ks) +
330 (abar_phiks_ * Ampbar_phiks) +
331 (abar_fxks_ * Ampbar_fxks) +
332 (abar_chic0ks_ * Ampbar_chic0ks) +
333 (abar_kpkmnr_ * Ampbar_kpkmnr) +
334 (abar_kskpnr_ * Ampbar_kskpnr) +
335 (abar_kskmnr_ * Ampbar_kskmnr);
336
337 // CP asymmetry
338 const double dm = getArg(0);
339
340 EvtComplex amp;
341 if (other_b == B0B) {
342 amp =
343 A_B0toKsKK * cos(dm * t / (2.0 * EvtConst::c)) +
344 EvtComplex(0.0, 1.0) * Abar_B0toKsKK * sin(dm * t / (2.0 * EvtConst::c));
345 }
346 if (other_b == B0) {
347 amp =
348 A_B0toKsKK *
349 EvtComplex(0.0, 1.0) * sin(dm * t / (2.0 * EvtConst::c)) +
350 Abar_B0toKsKK * cos(dm * t / (2.0 * EvtConst::c));
351 }
352
353 if (abs2(amp) > 50000.0)
354 debugfile_ << abs2(amp) << std::endl;
355 vertex(amp);
356
357 return;
358 }
359
360 EvtVector4R EvtB0toKsKK::umu(const EvtVector4R& p4a, const EvtVector4R& p4b,
361 const EvtVector4R& p4c)
362 {
363 const double s = (p4a + p4b + p4c) * (p4a + p4b + p4c);
364
365 const EvtVector4R umu = (p4a + p4b + p4c) / sqrt(s);
366
367 return umu;
368 }
369
370 EvtVector4R EvtB0toKsKK::Smu(const EvtVector4R& p4a, const EvtVector4R& p4b,
371 const EvtVector4R& /*p4c*/)
372 {
373 const double ma = p4a.mass();
374 const double mb = p4b.mass();
375 const double mres = (p4a + p4b).mass();
376 const double ma2 = ma * ma;
377 const double mb2 = mb * mb;
378 const double mres2 = mres * mres;
379
380 const double N1 =
381 mres /
382 (sqrt(mres2 - ((ma + mb) * (ma + mb))) * sqrt(mres2 - ((ma - mb) * (ma - mb))));
383
384 const EvtVector4R Smu =
385 N1 * ((p4a - p4b) - ((p4a + p4b) * (ma2 - mb2) / mres2));
386
387 return Smu;
388 }
389
390 EvtVector4R EvtB0toKsKK::Lmu(const EvtVector4R& p4a, const EvtVector4R& p4b,
391 const EvtVector4R& p4c)
392 {
393 const double mc = p4c.mass();
394 const double mres = (p4a + p4b).mass();
395 const double mc2 = mc * mc;
396 const double mres2 = mres * mres;
397 const double s = (p4a + p4b + p4c) * (p4a + p4b + p4c);
398
399 const double N2 =
400 sqrt(s) /
401 (sqrt(s - ((mres + mc) * (mres + mc))) * sqrt(s - ((mres - mc) * (mres - mc))));
402
403 const EvtVector4R Lmu =
404 N2 * ((p4c - (p4a + p4b)) - ((p4c + (p4a + p4b)) * (mc2 - mres2) / s));
405
406 return Lmu;
407 }
408
409 EvtTensor4C EvtB0toKsKK::gmunu_tilde(const EvtVector4R& p4a,
410 const EvtVector4R& p4b,
411 const EvtVector4R& p4c)
412 {
413 const double s = (p4a + p4b + p4c) * (p4a + p4b + p4c);
414
415 const EvtVector4R umu = (p4a + p4b + p4c) / sqrt(s);
416
417 const EvtTensor4C gmunu_tilde =
418 EvtTensor4C::g() - EvtGenFunctions::directProd(umu, umu);
419
420 return gmunu_tilde;
421 }
422
423 EvtTensor4C EvtB0toKsKK::Tmunu(const EvtVector4R& p4a, const EvtVector4R& p4b,
424 const EvtVector4R& p4c)
425 {
426 const double mres = (p4a + p4b).mass();
427 const EvtVector4R wmu = (p4a + p4b) / mres;
428
429 const EvtTensor4C Tmunu =
430 sqrt(3.0 / 2.0) *
431 (EvtGenFunctions::directProd(Smu(p4a, p4b, p4c), Smu(p4a, p4b, p4c)) +
432 (1.0 / 3.0 * (EvtTensor4C::g() - EvtGenFunctions::directProd(wmu, wmu))));
433
434 return Tmunu;
435 }
436
437 EvtTensor4C EvtB0toKsKK::Multiply(const EvtTensor4C& t1,
438 const EvtTensor4C& t2)
439 {
440 EvtTensor4C t;
441 for (unsigned int i = 0; i < 4; i++)
442 for (unsigned int j = 0; j < 4; j++) {
443 const EvtComplex element =
444 t1.get(i, 0) * t2.get(0, j) +
445 t1.get(i, 1) * t2.get(1, j) +
446 t1.get(i, 2) * t2.get(2, j) +
447 t1.get(i, 3) * t2.get(3, j);
448 t.set(i, j, element);
449 }
450 return t;
451 }
452
453 EvtTensor4C EvtB0toKsKK::RaiseIndices(const EvtTensor4C& t)
454 {
455 const EvtTensor4C tmp = Multiply(t, EvtTensor4C::g());
456 const EvtTensor4C t_raised = Multiply(EvtTensor4C::g(), tmp);
457
458 return t_raised;
459 }
460
461 void EvtB0toKsKK::RaiseIndex(EvtVector4R& vector)
462 {
463 vector.set(1, -vector.get(1));
464 vector.set(2, -vector.get(2));
465 vector.set(3, -vector.get(3));
466 }
467
468 EvtTensor4C EvtB0toKsKK::Mmunu(const EvtVector4R& p4a, const EvtVector4R& p4b,
469 const EvtVector4R& p4c)
470 {
471 const EvtTensor4C Mmunu =
472 sqrt(3.0 / 2.0) *
473 (EvtGenFunctions::directProd(Lmu(p4a, p4b, p4c), Lmu(p4a, p4b, p4c)) +
474 (1.0 / 3.0 * gmunu_tilde(p4a, p4b, p4c)));
475
476 return Mmunu;
477 }
478
479 double EvtB0toKsKK::BWBF(const double& q, const unsigned int& L)
480 {
481 //Meson radius
482 const double d = 1.0;
483
484 const double z = q * q * d * d;
485
486 double bwbf = 1.0;
487
488 if (L == 1)
489 bwbf = sqrt(2.0 * z / (1.0 + z));
490 if (L == 2)
491 bwbf = sqrt(13.0 * z * z / (((z - 3.0) * (z - 3.0)) + (9.0 * z)));
492 if (L > 2) {
493 std::cout << "ERROR: BWBF not implemented for L>2" << std::endl;
494 exit(1);
495 }
496 return bwbf;
497 }
498
499 double EvtB0toKsKK::BWBF(const double& q, const double& q0,
500 const unsigned int& L)
501 {
502 //Meson radius
503 const double d = 1.0;
504
505 const double z = q * q * d * d;
506 const double z0 = q0 * q0 * d * d;
507
508 double bwbf = 1.0;
509
510 if (L == 1)
511 bwbf = sqrt((1.0 + z0) / (1.0 + z));
512 if (L == 2)
513 bwbf = sqrt((((z0 - 3.0) * (z0 - 3.0)) + (9.0 * z0)) / (((z - 3.0) * (z - 3.0)) + (9.0 * z)));
514 if (L > 2) {
515 std::cout << "ERROR: BWBF not implemented for L>2" << std::endl;
516 exit(1);
517 }
518 return bwbf;
519 }
520
521 EvtComplex EvtB0toKsKK::BreitWigner(const double& m, const double& m0,
522 const double& Gamma0,
523 const double& q, const double& q0,
524 const unsigned int& L)
525 {
526 const double s = m * m;
527 const double s0 = m0 * m0;
528
529 const double Gamma =
530 Gamma0 * m0 / m * pow(q / q0, (2 * L) + 1) *
531 BWBF(q, q0, L) * BWBF(q, q0, L);
532
533 const EvtComplex BreitWigner = 1.0 / EvtComplex(s0 - s, -m0 * Gamma);
534
535 return BreitWigner;
536 }
537
538 EvtVector4R EvtB0toKsKK::Boost(const EvtVector4R& p4,
539 const EvtVector4R& boost)
540 {
541 const double bx = boost.get(1) / boost.get(0);
542 const double by = boost.get(2) / boost.get(0);
543 const double bz = boost.get(3) / boost.get(0);
544
545 const double bx2 = bx * bx;
546 const double by2 = by * by;
547 const double bz2 = bz * bz;
548
549 const double b2 = bx2 + by2 + bz2;
550 if (b2 == 0.0)
551 return p4;
552 assert(b2 < 1.0);
553
554 const double gamma = 1.0 / sqrt(1 - b2);
555
556 const double gb2 = (gamma - 1.0) / b2;
557
558 const double gb2xy = gb2 * bx * by;
559 const double gb2xz = gb2 * bx * bz;
560 const double gb2yz = gb2 * by * bz;
561
562 const double gbx = gamma * bx;
563 const double gby = gamma * by;
564 const double gbz = gamma * bz;
565
566 const double e_b =
567 (gamma * p4.get(0)) - (gbx * p4.get(1)) - (gby * p4.get(2)) - (gbz * p4.get(3));
568 const double x_b =
569 (-gbx * p4.get(0)) + ((1.0 + (gb2 * bx2)) * p4.get(1)) +
570 (gb2xy * p4.get(2)) + (gb2xz * p4.get(3));
571 const double y_b =
572 (-gby * p4.get(0)) + (gb2xy * p4.get(1)) +
573 ((1.0 + (gb2 * by2)) * p4.get(2)) + (gb2yz * p4.get(3));
574 const double z_b =
575 (-gbz * p4.get(0)) + (gb2xz * p4.get(1)) +
576 (gb2yz * p4.get(2)) + ((1.0 + (gb2 * bz2)) * p4.get(3));
577
578 const EvtVector4R p4_b(e_b, x_b, y_b, z_b);
579
580 return p4_b;
581 }
582
583 double EvtB0toKsKK::p(const double& mab, const double& M, const double& mc)
584 {
585 const double sab = mab * mab;
586 const double M_p_mc2 = (M + mc) * (M + mc);
587 const double M_m_mc2 = (M - mc) * (M - mc);
588
589 const double p = sqrt((sab - M_p_mc2) * (sab - M_m_mc2)) / (2.0 * mab);
590
591 return p;
592 }
593
594 double EvtB0toKsKK::q(const double& mab, const double& ma, const double& mb)
595 {
596 const double mab2 = mab * mab;
597 const double ma_p_mb2 = (ma + mb) * (ma + mb);
598 const double ma_m_mb2 = (ma - mb) * (ma - mb);
599
600 const double q = sqrt((mab2 - ma_p_mb2) * (mab2 - ma_m_mb2)) / (2.0 * mab);
601
602 return q;
603 }
604
605 EvtComplex EvtB0toKsKK::Flatte_k(const double& s, const double& m_h)
606 {
607 const double k2 = 1.0 - (4.0 * m_h * m_h / s);
608 EvtComplex k;
609 if (k2 < 0.0)
610 k = EvtComplex(0.0, sqrt(fabs(k2)));
611 else
612 k = sqrt(k2);
613
614 return k;
615 }
616
617 EvtComplex EvtB0toKsKK::Flatte(const double& m, const double& m0)
618 {
619 const double g_pipi = 0.165;
620 const double g_kk = 4.21 * g_pipi;
621
622 static EvtId pip = EvtPDL::getId("pi+");
623 static EvtId kp = EvtPDL::getId("K+");
624
625 const double s = m * m;
626 const double s0 = m0 * m0;
627
628 const EvtComplex rho_pipi = 2.0 * Flatte_k(s, EvtPDL::getMeanMass(pip)) / m;
629 const EvtComplex rho_kk = 2.0 * Flatte_k(s, EvtPDL::getMeanMass(kp)) / m;
630
631 const EvtComplex Gamma =
632 EvtComplex(0.0, 1.0) * ((g_pipi * rho_pipi) + (g_kk * rho_kk));
633
634 const EvtComplex Flatte = 1.0 / (s0 - s - Gamma);
635
636 return Flatte;
637 }
638
639 EvtComplex EvtB0toKsKK::A_f0ks(const EvtVector4R& /*p4ks*/,
640 const EvtVector4R& p4kp,
641 const EvtVector4R& p4km)
642 {
643 static EvtId f0 = EvtPDL::getId("f_0");
644
645 const double f0_m = EvtPDL::getMeanMass(f0);
646
647 const EvtVector4R p4kpkm = p4kp + p4km;
648 const double mkpkm = p4kpkm.mass();
649
650 //Angular distribution
651 const EvtComplex H_f0ks = 1.0;
652
653 //Barrier factors
654 const EvtComplex D_f0ks = 1.0;
655
656 //Line shape
657 const EvtComplex L_f0ks = Flatte(mkpkm, f0_m);
658
659 //Amplitude
660 const EvtComplex A_f0ks = D_f0ks * H_f0ks * L_f0ks;
661
662 return A_f0ks;
663 }
664
666 double s13_min(const double& s12, const double& M,
667 const double& m1, const double& m2, const double& m3)
668 {
669 const double E2star = (s12 - (m2 * m2) + (m1 * m1)) / 2.0 / sqrt(s12);
670 const double E3star = (M * M - s12 - (m3 * m3)) / 2.0 / sqrt(s12);
671
672 const double s23_min =
673 ((E2star + E3star) * (E2star + E3star)) -
674 ((sqrt((E2star * E2star) - (m1 * m1)) + sqrt((E3star * E3star) - (m3 * m3))) *
675 (sqrt((E2star * E2star) - (m1 * m1)) + sqrt((E3star * E3star) - (m3 * m3))));
676
677 return s23_min;
678 }
679
681 double s13_max(const double& s12, const double& M,
682 const double& m1, const double& m2, const double& m3)
683 {
684 const double E2star = (s12 - (m2 * m2) + (m1 * m1)) / 2.0 / sqrt(s12);
685 const double E3star = (M * M - s12 - (m3 * m3)) / 2.0 / sqrt(s12);
686
687 const double s23_max =
688 ((E2star + E3star) * (E2star + E3star)) -
689 ((sqrt((E2star * E2star) - (m1 * m1)) - sqrt((E3star * E3star) - (m3 * m3))) *
690 (sqrt((E2star * E2star) - (m1 * m1)) - sqrt((E3star * E3star) - (m3 * m3))));
691
692 return s23_max;
693 }
694
695 EvtComplex EvtB0toKsKK::A_phiks(const EvtVector4R& p4ks,
696 const EvtVector4R& p4kp,
697 const EvtVector4R& p4km)
698 {
699 static EvtId phi = EvtPDL::getId("phi");
700
701 const double phi_m = EvtPDL::getMeanMass(phi);
702 const double phi_w = EvtPDL::getWidth(phi);
703
704 const EvtVector4R p4kpkm = p4kp + p4km;
705 const double mkpkm = p4kpkm.mass();
706
707 //Angular distribution
708 const EvtVector4R S_mu = Smu(p4kp, p4km, p4ks);
709 const EvtVector4R L_mu = Lmu(p4kp, p4km, p4ks);
710
711 const EvtTensor4C g_munu_tilde = gmunu_tilde(p4kp, p4km, p4ks);
712 const EvtTensor4C g__munu_tilde = RaiseIndices(g_munu_tilde);
713
714 EvtComplex H_phiks = 0.0;
715 for (unsigned int i = 0; i < 4; i++)
716 for (unsigned int j = 0; j < 4; j++)
717 H_phiks += g__munu_tilde.get(i, j) * S_mu.get(i) * L_mu.get(j);
718
719 //Barrier factors
720 const EvtVector4R p4kp_kpkm = Boost(p4kp, p4kpkm);
721
722 const EvtComplex D_phiks = BWBF(p4kp_kpkm.d3mag(), 1);
723
724 //Line shape
725 const double q0_phi = q(phi_m, p4kp.mass(), p4km.mass());
726 const EvtComplex L_phiks =
727 BreitWigner(mkpkm, phi_m, phi_w, p4kp_kpkm.d3mag(), q0_phi, 1);
728
729 //Amplitude
730 const EvtComplex A_phiks = D_phiks * H_phiks * L_phiks;
731
733 /*const EvtVector4R p4ks_kpkm = Boost(p4ks, p4kpkm);
734 const double p4ks_kpkm_alt = p(mkpkm, (p4ks+p4kp+p4km).mass(), p4ks.mass());
735 const double p4kp_kpkm_alt = q(mkpkm, p4kp.mass(), p4km.mass());
736 //std::cout << p4ks_kpkm.d3mag() << std::endl;
737 //std::cout << p4ks_kpkm_alt << std::endl;
738 //std::cout << p4kp_kpkm.d3mag() << std::endl;
739 //std::cout << p4kp_kpkm_alt << std::endl;
740 const double costheta =
741 p4ks_kpkm.dot(p4kp_kpkm)/p4ks_kpkm.d3mag()/p4kp_kpkm.d3mag();
742 //const double z = p4ks_kpkm.d3mag()/(p4ks+p4kp+p4km).mass();
743 //const double H_phiks_alt = -sqrt(1.0+(z*z))*costheta;
744 const double sp_min =
745 s13_min( mkpkm*mkpkm, (p4ks+p4kp+p4km).mass(),
746 p4kp.mass(), p4km.mass(), p4ks.mass() );
747 const double sp_max =
748 s13_max( mkpkm*mkpkm, (p4ks+p4kp+p4km).mass(),
749 p4kp.mass(), p4km.mass(), p4ks.mass() );
750 const double sp = (p4ks+p4kp).mass2();
751 const double costheta_alt = (sp_max + sp_min - (2.0*sp))/(sp_max - sp_min);
752 //std::cout << costheta << std::endl;
753 //std::cout << costheta_alt << std::endl;
754 //exit(1);
755 const double z_alt = p4ks_kpkm_alt/(p4ks+p4kp+p4km).mass();
756 const double H_phiks_alt = -sqrt(1.0+(z_alt*z_alt))*costheta_alt;
757 //std::cout << H_phiks << std::endl;
758 //std::cout << H_phiks_alt << std::endl;
759 if( real(H_phiks) != H_phiks_alt )
760 {
761 std::cout << H_phiks << std::endl;
762 std::cout << H_phiks_alt << std::endl;
763 }
764 const EvtComplex A_phiks = D_phiks*H_phiks_alt*L_phiks;*/
766
767 return A_phiks;
768 }
769
770 EvtComplex EvtB0toKsKK::A_fxks(const EvtVector4R& /*p4ks*/,
771 const EvtVector4R& p4kp,
772 const EvtVector4R& p4km)
773 {
774 static EvtId fx = EvtPDL::getId("f_0(1500)");
775
776 const double fx_m = EvtPDL::getMeanMass(fx);
777 const double fx_w = EvtPDL::getWidth(fx);
778
779 const EvtVector4R p4kpkm = p4kp + p4km;
780 const double mkpkm = p4kpkm.mass();
781
782 //Angular distribution
783 const EvtComplex H_fxks = 1.0;
784
785 //Barrier factors
786 const EvtComplex D_fxks = 1.0;
787
788 //Line shape
789 const EvtVector4R p4kp_kpkm = Boost(p4kp, p4kpkm);
790 const double q0_fx = q(fx_m, p4kp.mass(), p4km.mass());
791 const EvtComplex L_fxks =
792 BreitWigner(mkpkm, fx_m, fx_w, p4kp_kpkm.d3mag(), q0_fx, 0);
793
794 //Amplitude
795 const EvtComplex A_fxks = D_fxks * H_fxks * L_fxks;
796
797 return A_fxks;
798 }
799
800 EvtComplex EvtB0toKsKK::A_chic0ks(const EvtVector4R& /*p4ks*/,
801 const EvtVector4R& p4kp,
802 const EvtVector4R& p4km)
803 {
804 static EvtId chic0 = EvtPDL::getId("chi_c0");
805
806 const double chic0_m = EvtPDL::getMeanMass(chic0);
807 const double chic0_w = EvtPDL::getWidth(chic0);
808
809 const EvtVector4R p4kpkm = p4kp + p4km;
810 const double mkpkm = p4kpkm.mass();
811
812 //Angular distribution
813 const EvtComplex H_chic0ks = 1.0;
814
815 //Barrier factors
816 const EvtComplex D_chic0ks = 1.0;
817
818 //Line shape
819 const EvtVector4R p4kp_kpkm = Boost(p4kp, p4kpkm);
820 const double q0_chic0 = q(chic0_m, p4kp.mass(), p4km.mass());
821 const EvtComplex L_chic0ks =
822 BreitWigner(mkpkm, chic0_m, chic0_w, p4kp_kpkm.d3mag(), q0_chic0, 0);
823
824 //Amplitude
825 const EvtComplex A_chic0ks = D_chic0ks * H_chic0ks * L_chic0ks;
826
827 return A_chic0ks;
828 }
829
830 EvtComplex EvtB0toKsKK::A_kknr(const EvtVector4R& p4k1,
831 const EvtVector4R& p4k2,
832 const double& alpha_kk)
833 {
834 const EvtVector4R p4kk = p4k1 + p4k2;
835 const double m2kk = p4kk.mass2();
836
837 const EvtComplex A_kknr = exp(-alpha_kk * m2kk);
838
839 return A_kknr;
840 }
841
843} // Belle 2 Namespace
Register Decay model EvtB0toKsKK.
Definition: EvtB0toKsKK.h:21
EvtComplex abar_kskpnr_
Variable member abar_kskpnr_
Definition: EvtB0toKsKK.h:100
double alpha_kskmnr
Variable member alpha_kskmnr.
Definition: EvtB0toKsKK.h:105
EvtComplex a_chic0ks_
Variable member a_chic0ks_.
Definition: EvtB0toKsKK.h:90
EvtComplex a_fxks_
Variable member a_fxks_
Definition: EvtB0toKsKK.h:89
std::ofstream debugfile_
debuging stream
Definition: EvtB0toKsKK.h:107
EvtComplex a_f0ks_
<Variable names for form factors
Definition: EvtB0toKsKK.h:87
double alpha_kpkmnr
Variable member alpha_kpkmnr.
Definition: EvtB0toKsKK.h:103
EvtComplex abar_fxks_
Variable member abar_fxks_
Definition: EvtB0toKsKK.h:97
EvtComplex abar_f0ks_
Variable member abar_f0ks_
Definition: EvtB0toKsKK.h:95
double alpha_kskpnr
Variable member alpha_kskpnr.
Definition: EvtB0toKsKK.h:104
EvtComplex abar_chic0ks_
Variable member abar_chic0ks_.
Definition: EvtB0toKsKK.h:98
EvtComplex a_kskmnr_
Variable member a_kskmnr_.
Definition: EvtB0toKsKK.h:93
EvtComplex abar_kskmnr_
Variable member abar_kskmnr_
Definition: EvtB0toKsKK.h:101
EvtComplex a_kpkmnr_
Variable member a_kpkmnr_.
Definition: EvtB0toKsKK.h:91
EvtComplex abar_phiks_
Variable member abar_phiks_.
Definition: EvtB0toKsKK.h:96
EvtComplex a_kskpnr_
Variable member a_kskpnr_.
Definition: EvtB0toKsKK.h:92
EvtComplex a_phiks_
Variable member a_phiks_
Definition: EvtB0toKsKK.h:88
EvtComplex abar_kpkmnr_
Variable member abar_kpkmnr_
Definition: EvtB0toKsKK.h:99
void init()
Initialize standard stream objects
Definition: EvtB0toKsKK.cc:41
EvtVector4R Smu(const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
Function 4Vector Smu.
Definition: EvtB0toKsKK.cc:370
EvtComplex A_f0ks(const EvtVector4R &p4ks, const EvtVector4R &p4kp, const EvtVector4R &p4km)
A_f0ks is amplitude of f0.
Definition: EvtB0toKsKK.cc:639
EvtComplex A_phiks(const EvtVector4R &p4ks, const EvtVector4R &p4kp, const EvtVector4R &p4km)
A_phiks is amplitude of phi.
Definition: EvtB0toKsKK.cc:695
EvtVector4R umu(const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
Function 4Vector umu.
Definition: EvtB0toKsKK.cc:360
EvtComplex Flatte(const double &m, const double &m0)
Constant Flatte.
Definition: EvtB0toKsKK.cc:617
double s13_max(const double &s12, const double &M, const double &m1, const double &m2, const double &m3)
maximum s13
Definition: EvtB0toKsKK.cc:681
double q(const double &mab, const double &ma, const double &mb)
Constants q.
Definition: EvtB0toKsKK.cc:594
EvtTensor4C gmunu_tilde(const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
Function Tensor gmunu
Definition: EvtB0toKsKK.cc:409
EvtComplex BreitWigner(const double &m, const double &m0, const double &Gamma0, const double &q, const double &q0, const unsigned int &L)
BreitWigner Shape.
Definition: EvtB0toKsKK.cc:521
EvtTensor4C Mmunu(const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
Function Tensor Mmunu.
Definition: EvtB0toKsKK.cc:468
EvtTensor4C RaiseIndices(const EvtTensor4C &t)
Function RaiseIndices
Definition: EvtB0toKsKK.cc:453
void RaiseIndex(EvtVector4R &vector)
Member function RaiseIndices.
Definition: EvtB0toKsKK.cc:461
EvtDecayBase * clone()
Clone the decay of B0toKsKK.
Definition: EvtB0toKsKK.cc:36
EvtVector4R Lmu(const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
Function 4Vector Lmu.
Definition: EvtB0toKsKK.cc:390
EvtComplex A_fxks(const EvtVector4R &p4ks, const EvtVector4R &p4kp, const EvtVector4R &p4km)
A_fxks is amplitude of fxks.
Definition: EvtB0toKsKK.cc:770
double p(const double &mab, const double &M, const double &mc)
Constants p
Definition: EvtB0toKsKK.cc:583
EvtComplex Flatte_k(const double &s, const double &m_h)
Constant Flatte_k.
Definition: EvtB0toKsKK.cc:605
EvtComplex A_kknr(const EvtVector4R &p4k1, const EvtVector4R &p4k2, const double &alpha_kk)
A_kknr is amplitude of kknr.
Definition: EvtB0toKsKK.cc:830
void initProbMax()
Initialize standard stream objects for probability function
Definition: EvtB0toKsKK.cc:136
double BWBF(const double &q, const unsigned int &L)
Meson radius
Definition: EvtB0toKsKK.cc:479
EvtTensor4C Tmunu(const EvtVector4R &p4a, const EvtVector4R &p4b, const EvtVector4R &p4c)
Function Tensor Tmunu
Definition: EvtB0toKsKK.cc:423
EvtVector4R Boost(const EvtVector4R &p4, const EvtVector4R &boost)
Parameter for boost frame
Definition: EvtB0toKsKK.cc:538
std::string getName()
Get function Name
Definition: EvtB0toKsKK.cc:31
#define B2_EVTGEN_REGISTER_MODEL(classname)
Class to register B2_EVTGEN_REGISTER_MODEL.
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
void decay(EvtParticle *p)
Member of particle in EvtGen.
Definition: EvtB0toKsKK.cc:146
EvtTensor4C Multiply(const EvtTensor4C &t1, const EvtTensor4C &t2)
Function Tensor Multiply
Definition: EvtB0toKsKK.cc:437
double s13_min(const double &s12, const double &M, const double &m1, const double &m2, const double &m3)
minimum s13
Definition: EvtB0toKsKK.cc:666
EvtComplex A_chic0ks(const EvtVector4R &p4ks, const EvtVector4R &p4kp, const EvtVector4R &p4km)
A_chic0ks is amplitude of chic0ks.
Definition: EvtB0toKsKK.cc:800
Abstract base class for different kinds of events.