30#include "EvtGenBase/EvtConst.hh"
31#include "EvtGenBase/EvtParticle.hh"
32#include "EvtGenBase/EvtPDL.hh"
33#include "EvtGenBase/EvtReport.hh"
34#include "EvtGenBase/EvtVector4C.hh"
35#include "generators/evtgen/models/EvtPHSPBMix.h"
36#include <generators/evtgen/EvtGenModelRegister.h>
37#include "EvtGenBase/EvtId.hh"
40#include "EvtGenBase/EvtRandom.hh"
64 const unsigned short int narg(2);
65 if (getNArg() != narg) {
66 EvtGenReport(EVTGEN_ERROR,
"EvtGen") <<
"EvtPHSPBBMix generator expected "
68 <<
"argument(s) (deltam, C=-1, +1 or 0 (incoherent)) but found:"
70 EvtGenReport(EVTGEN_ERROR,
"EvtGen") <<
"Will terminate execution!" << endl;
82 if (getNDaug() < 2 || getNDaug() > 5) {
83 EvtGenReport(EVTGEN_ERROR,
"EvtGen") <<
"PHSP_BB_MIX n daughters not ok :"
84 << getNDaug() << endl;
85 EvtGenReport(EVTGEN_ERROR,
"EvtGen") <<
"Will terminate execution!" << endl;
95 if (getNDaug() == 4 && (EvtPDL::chargeConj(getDaug(0)) == getDaug(1)))
99 if (!
_BBpipi && getNDaug() > 3) {
100 if (!(EvtPDL::chargeConj(getDaug(0)) == getDaug(getNDaug() - 2) && EvtPDL::chargeConj(getDaug(1)) == getDaug(getNDaug() - 1))) {
101 EvtGenReport(EVTGEN_ERROR,
"EvtGen") <<
"EvtPHSPBBMix generator expected daughters "
102 <<
"to be charge conjugate." << endl
103 <<
" Found " << EvtPDL::name(getDaug(0)).c_str()
104 <<
"," << EvtPDL::name(getDaug(1)).c_str()
105 <<
"," << EvtPDL::name(getDaug(2)).c_str()
107 << EvtPDL::name(getDaug(3)).c_str() << endl;
108 EvtGenReport(EVTGEN_ERROR,
"EvtGen") <<
"Will terminate execution!" << endl;
113 if (!(EvtPDL::chargeConj(getDaug(0)) == getDaug(1))) {
114 EvtGenReport(EVTGEN_ERROR,
"EvtGen") <<
"EvtPHSPBBMix generator expected daughters "
115 <<
"to be charge conjugate." << endl
116 <<
" Found " << EvtPDL::name(getDaug(0)).c_str()
118 << EvtPDL::name(getDaug(1)).c_str() << endl;
119 EvtGenReport(EVTGEN_ERROR,
"EvtGen") <<
"Will terminate execution!" << endl;
137 if (EvtPDL::getctau(getDaug(0)) != EvtPDL::getctau(getDaug(1))) {
138 EvtGenReport(EVTGEN_ERROR,
"EvtGen") <<
"EvtPHSPBBMix generator expected daughters "
139 <<
"to have the same lifetime." << endl
141 << EvtPDL::getctau(getDaug(0)) <<
" mm and "
142 << EvtPDL::getctau(getDaug(1)) <<
" mm" << endl;
143 EvtGenReport(EVTGEN_ERROR,
"EvtGen") <<
"Will terminate execution!" << endl;
151 _freq = getArg(0) / EvtConst::c;
161 double tau = 1e12 * EvtPDL::getctau(getDaug(0)) / EvtConst::c;
162 double dm = 1e-12 * getArg(0);
168 std::ostringstream sss;
169 sss <<
" " << EvtPDL::name(getParentId()).c_str() <<
" --> "
170 << EvtPDL::name(getDaug(0)).c_str() <<
" + "
171 << EvtPDL::name(getDaug(1)).c_str();
172 switch (getNDaug()) {
175 sss <<
" + " << EvtPDL::name(getDaug(2)).c_str();
179 sss <<
" + " << EvtPDL::name(getDaug(2)).c_str() <<
" + " << EvtPDL::name(getDaug(3)).c_str();
183 EvtGenReport(EVTGEN_INFO,
"EvtGen") <<
"PHSP_BB_MIX will generate mixing :"
185 << sss.str() << endl << endl
186 <<
"using parameters:" << endl << endl
187 <<
" delta(m) = " << dm <<
" hbar/ps" << endl
188 <<
" C (B0-B0b) = " <<
_C << endl
189 <<
" _freq = " <<
_freq <<
" hbar/mm" << endl
190 <<
" dgog = " << 0 << endl
191 <<
" dGamma = " << 0 <<
" hbar/mm" << endl
192 <<
" q/p = " << 1 << endl
193 <<
" tau = " << tau <<
" ps" << endl
194 <<
" x = " << x << endl
205 if (getNDaug() > 3) {
209 return getNDaug() - 2;
227 EvtGenReport(EVTGEN_INFO,
"EvtGen") <<
"decay p" << i << endl;
237 static const EvtId B0(EvtPDL::getId(
"B0"));
238 static const EvtId B0B(EvtPDL::getId(
"anti-B0"));
245 const bool NCC(getNDaug() >= 4 && !
_BBpipi);
250 std::vector<EvtId> tempDaug(getNDaug() - 2);
251 tempDaug[0] = getDaug(0);
252 tempDaug[1] = getDaug(1);
254 tempDaug[2] = getDaug(2);
256 p->initializePhaseSpace(getNDaug() - 2, tempDaug.data());
260 p->initializePhaseSpace(getNDaug(), getDaugs());
264 EvtParticle* s1, *s2;
270 if (s1->getNDaug() > 0) { s1->deleteDaughters();}
271 if (s2->getNDaug() > 0) { s2->deleteDaughters();}
273 EvtVector4R p1 = s1->getP4();
274 EvtVector4R p2 = s2->getP4();
277 const EvtId B1(EvtRandom::random() > 0.5 ? B0 : B0B);
278 const EvtId B2(EvtRandom::random() > 0.5 ? B0 : B0B);
284 if (getNDaug() == 4) {
285 s1->init(B1 == B0 ? getDaug(0) : getDaug(2), p1);
286 s2->init(B2 == B0 ? getDaug(3) : getDaug(1), p2);
288 if (getNDaug() == 5) {
289 s1->init(B1 == B0 ? getDaug(0) : getDaug(3), p1);
290 s2->init(B2 == B0 ? getDaug(4) : getDaug(1), p2);
293 s1->init(B1 == B0 ? getDaug(0) : getDaug(1), p1);
294 s2->init(B2 == B0 ? getDaug(0) : getDaug(1), p2);
305 const double t1(s1->getLifetime());
306 const double t2(s2->getLifetime());
309 EvtComplex osc_amp(
Amplitude(t1, t2, B1 == B0, B2 == B0));
314 double norm = 1.0 / p1.d3mag();
316 vertex(0, norm * osc_amp * p1 * (p->eps(0)));
317 vertex(1, norm * osc_amp * p1 * (p->eps(1)));
318 vertex(2, norm * osc_amp * p1 * (p->eps(2)));
325 if (
_C != 0 &&
_C != -1 &&
_C != 1)
326 return EvtComplex(0., 0.);
328 const double f(
_freq);
329 if (B1_is_B0 && !B2_is_B0) {
331 return EvtComplex(cos(f * t1 / 2.) * cos(f * t2 / 2.), 0.);
333 return EvtComplex(cos(f * (t2 +
_C * t1) / 2.) / sqrt(2.), 0.);
335 if (!B1_is_B0 && B2_is_B0) {
337 return EvtComplex(-sin(f * t1 / 2.) * sin(f * t2 / 2.), 0.);
339 return EvtComplex(cos(f * (t2 +
_C * t1) / 2.) *
_C / sqrt(2.), 0.);
341 if (B1_is_B0 && B2_is_B0) {
343 return EvtComplex(0., -cos(f * t1 / 2.) * sin(f * t2 / 2.));
345 return EvtComplex(0., -sin(f * (t2 +
_C * t1) / 2.) / sqrt(2.));
347 if (!B1_is_B0 && !B2_is_B0) {
349 return EvtComplex(0., -sin(f * t1 / 2.) * cos(f * t2 / 2.));
351 return EvtComplex(0., -sin(f * (t2 +
_C * t1) / 2.) *
_C / sqrt(2.));
354 return EvtComplex(0., 0.);
382 const unsigned short int narg(1);
383 if (getNArg() != narg) {
384 EvtGenReport(EVTGEN_ERROR,
"EvtGen") <<
"EvtPHSPBMix generator expected "
386 <<
" argument(s) (delta m) but found:"
387 << getNArg() << endl;
388 EvtGenReport(EVTGEN_ERROR,
"EvtGen") <<
"Will terminate execution!" << endl;
401 bool BBpipi(getNDaug() == 5);
403 if (!(EvtPDL::chargeConj(getDaug(0)) == getDaug(getNDaug() - 1))) {
404 EvtGenReport(EVTGEN_ERROR,
"EvtGen") <<
"EvtPHSPBMix generator expected two first particles "
405 <<
"to be charge conjugate." << endl
406 <<
" Found " << EvtPDL::name(getDaug(0)).c_str()
407 <<
"," << EvtPDL::name(getDaug(getNDaug() - 1)).c_str() << endl;
408 EvtGenReport(EVTGEN_ERROR,
"EvtGen") <<
"Will terminate execution!" << endl;
413 _freq = getArg(0) / EvtConst::c;
417 double tau = 1e12 * EvtPDL::getctau(getDaug(0)) / EvtConst::c;
418 double dm = 1e-12 * getArg(0);
421 std::ostringstream sss;
422 sss <<
" " << EvtPDL::name(getParentId()).c_str() <<
" --> "
423 << EvtPDL::name(getDaug(0)).c_str() <<
" + "
424 << EvtPDL::name(getDaug(1)).c_str();
425 sss <<
" + " << EvtPDL::name(getDaug(2)).c_str();
427 sss <<
" + " << EvtPDL::name(getDaug(3)).c_str();
429 EvtGenReport(EVTGEN_INFO,
"EvtGen") <<
"PHSP_B_MIX will generate mixing :"
431 << sss.str() << endl << endl
432 <<
"using parameters:" << endl << endl
433 <<
" delta(m) = " << dm <<
" hbar/ps" << endl
434 <<
" _freq = " <<
_freq <<
" hbar/mm" << endl
435 <<
" dgog = " << 0 << endl
436 <<
" dGamma = " << 0 <<
" hbar/mm" << endl
437 <<
" q/p = " << 1 << endl
438 <<
" tau = " << tau <<
" ps" << endl
439 <<
" x = " << x << endl
446 return getNDaug() - 1;
462 if (
_print_info) std::cout <<
"decay B_MIX (0)" << std::endl;
463 std::vector<EvtId> tempDaug(getNDaug() - 1);
464 tempDaug[0] = getDaug(0);
465 tempDaug[1] = getDaug(1);
466 tempDaug[2] = getDaug(2);
468 tempDaug[3] = getDaug(3);
470 if (
_print_info) std::cout <<
"decay B_MIX (1)" << std::endl;
473 p->initializePhaseSpace(getNDaug() - 1, tempDaug.data());
481 if (s1->getNDaug() > 0) { s1->deleteDaughters();}
483 EvtVector4R p1 = s1->getP4();
486 const bool changed_flavor(EvtRandom::random() > 0.5);
487 if (
_print_info) std::cout <<
"decay B_MIX (3)" << std::endl;
488 s1->init(changed_flavor ? getDaug(getNDaug() - 1) : getDaug(0), p1);
494 if (
_print_info) std::cout <<
"decay B_MIX (4)" << std::endl;
495 const double t1(s1->getLifetime());
498 EvtComplex osc_amp(changed_flavor ? EvtComplex(0,
499 -sin(
_freq * t1 / 2.)) : EvtComplex(cos(
_freq * t1 / 2.)));
500 if (
_print_info) std::cout <<
"decay B_MIX (5)" << std::endl;
502 double norm = 1.0 / p1.d3mag();
503 if (
_print_info) std::cout <<
"decay B_MIX (6)" << std::endl;
504 vertex(0, norm * osc_amp * p1 * (p->eps(0)));
505 vertex(1, norm * osc_amp * p1 * (p->eps(1)));
506 vertex(2, norm * osc_amp * p1 * (p->eps(2)));
507 if (
_print_info) std::cout <<
"decay B_MIX (7)" << std::endl;
The class provides routine to decay vector-> particle particle with B0 mixing, coherent B0B0-like mix...
void init()
Init function.
void prlp(int) const
Number of real daughters.
bool _print_info
Print evtgeninfo?
virtual ~EvtPHSPBBMix()
Default destructor.
EvtPHSPBBMix()
Default constructor.
EvtDecayBase * clone()
Clone the decay
void initProbMax()
Init maximal prob.
int nRealDaughters()
Number of real daughters.
double _C
C eigenvalue, 0= incoherent.
double _freq
mixing frequency in hbar/mm
EvtComplex Amplitude(const double &t1, const double &t2, bool B1_is_B0, bool B2_is_B0) const
Calculate amplitude.
std::string getName()
Get function Name
void decay(EvtParticle *p)
Decay function.
The class provides routine to decay vector-> particle particle with B0 mixing, handles states with on...
void init()
Init function.
bool _print_info
Print evtgeninfo?
EvtDecayBase * clone()
Clone the decay
virtual ~EvtPHSPBMix()
Default destructor.
void initProbMax()
Init maximal prob.
int nRealDaughters()
Number of real daughters.
EvtPHSPBMix()
Default constructor.
double _freq
mixing frequency in hbar/mm
std::string getName()
Get function Name
void decay(EvtParticle *p)
Decay function.
#define B2_EVTGEN_REGISTER_MODEL(classname)
Class to register B2_EVTGEN_REGISTER_MODEL.