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;
238 static const EvtId B0(EvtPDL::getId(
"B0"));
239 static const EvtId B0B(EvtPDL::getId(
"anti-B0"));
246 const bool NCC(getNDaug() >= 4 && !
_BBpipi);
251 std::vector<EvtId> tempDaug(getNDaug() - 2);
252 tempDaug[0] = getDaug(0);
253 tempDaug[1] = getDaug(1);
255 tempDaug[2] = getDaug(2);
257 p->initializePhaseSpace(getNDaug() - 2, tempDaug.data());
261 p->initializePhaseSpace(getNDaug(), getDaugs());
265 EvtParticle* s1, *s2;
271 if (s1->getNDaug() > 0) { s1->deleteDaughters();}
272 if (s2->getNDaug() > 0) { s2->deleteDaughters();}
274 EvtVector4R p1 = s1->getP4();
275 EvtVector4R p2 = s2->getP4();
278 const EvtId B1(EvtRandom::random() > 0.5 ? B0 : B0B);
279 const EvtId B2(EvtRandom::random() > 0.5 ? B0 : B0B);
285 if (getNDaug() == 4) {
286 s1->init(B1 == B0 ? getDaug(0) : getDaug(2), p1);
287 s2->init(B2 == B0 ? getDaug(3) : getDaug(1), p2);
289 if (getNDaug() == 5) {
290 s1->init(B1 == B0 ? getDaug(0) : getDaug(3), p1);
291 s2->init(B2 == B0 ? getDaug(4) : getDaug(1), p2);
294 s1->init(B1 == B0 ? getDaug(0) : getDaug(1), p1);
295 s2->init(B2 == B0 ? getDaug(0) : getDaug(1), p2);
306 const double t1(s1->getLifetime());
307 const double t2(s2->getLifetime());
310 EvtComplex osc_amp(
Amplitude(t1, t2, B1 == B0, B2 == B0));
315 double norm = 1.0 / p1.d3mag();
317 vertex(0, norm * osc_amp * p1 * (p->eps(0)));
318 vertex(1, norm * osc_amp * p1 * (p->eps(1)));
319 vertex(2, norm * osc_amp * p1 * (p->eps(2)));
326 if (
_C != 0 &&
_C != -1 &&
_C != 1)
327 return EvtComplex(0., 0.);
329 const double f(
_freq);
330 if (B1_is_B0 && !B2_is_B0) {
332 return EvtComplex(cos(f * t1 / 2.) * cos(f * t2 / 2.), 0.);
334 return EvtComplex(cos(f * (t2 +
_C * t1) / 2.) / sqrt(2.), 0.);
336 if (!B1_is_B0 && B2_is_B0) {
338 return EvtComplex(-sin(f * t1 / 2.) * sin(f * t2 / 2.), 0.);
340 return EvtComplex(cos(f * (t2 +
_C * t1) / 2.) *
_C / sqrt(2.), 0.);
342 if (B1_is_B0 && B2_is_B0) {
344 return EvtComplex(0., -cos(f * t1 / 2.) * sin(f * t2 / 2.));
346 return EvtComplex(0., -sin(f * (t2 +
_C * t1) / 2.) / sqrt(2.));
348 if (!B1_is_B0 && !B2_is_B0) {
350 return EvtComplex(0., -sin(f * t1 / 2.) * cos(f * t2 / 2.));
352 return EvtComplex(0., -sin(f * (t2 +
_C * t1) / 2.) *
_C / sqrt(2.));
355 return EvtComplex(0., 0.);
383 const unsigned short int narg(1);
384 if (getNArg() != narg) {
385 EvtGenReport(EVTGEN_ERROR,
"EvtGen") <<
"EvtPHSPBMix generator expected "
387 <<
" argument(s) (delta m) but found:"
388 << getNArg() << endl;
389 EvtGenReport(EVTGEN_ERROR,
"EvtGen") <<
"Will terminate execution!" << endl;
402 bool BBpipi(getNDaug() == 5);
404 if (!(EvtPDL::chargeConj(getDaug(0)) == getDaug(getNDaug() - 1))) {
405 EvtGenReport(EVTGEN_ERROR,
"EvtGen") <<
"EvtPHSPBMix generator expected two first particles "
406 <<
"to be charge conjugate." << endl
407 <<
" Found " << EvtPDL::name(getDaug(0)).c_str()
408 <<
"," << EvtPDL::name(getDaug(getNDaug() - 1)).c_str() << endl;
409 EvtGenReport(EVTGEN_ERROR,
"EvtGen") <<
"Will terminate execution!" << endl;
414 _freq = getArg(0) / EvtConst::c;
418 double tau = 1e12 * EvtPDL::getctau(getDaug(0)) / EvtConst::c;
419 double dm = 1e-12 * getArg(0);
422 std::ostringstream sss;
423 sss <<
" " << EvtPDL::name(getParentId()).c_str() <<
" --> "
424 << EvtPDL::name(getDaug(0)).c_str() <<
" + "
425 << EvtPDL::name(getDaug(1)).c_str();
426 sss <<
" + " << EvtPDL::name(getDaug(2)).c_str();
428 sss <<
" + " << EvtPDL::name(getDaug(3)).c_str();
430 EvtGenReport(EVTGEN_INFO,
"EvtGen") <<
"PHSP_B_MIX will generate mixing :"
432 << sss.str() << endl << endl
433 <<
"using parameters:" << endl << endl
434 <<
" delta(m) = " << dm <<
" hbar/ps" << endl
435 <<
" _freq = " <<
_freq <<
" hbar/mm" << endl
436 <<
" dgog = " << 0 << endl
437 <<
" dGamma = " << 0 <<
" hbar/mm" << endl
438 <<
" q/p = " << 1 << endl
439 <<
" tau = " << tau <<
" ps" << endl
440 <<
" x = " << x << endl
447 return getNDaug() - 1;
464 if (pr) std::cout <<
"decay B_MIX (0)" << std::endl;
465 std::vector<EvtId> tempDaug(getNDaug() - 1);
466 tempDaug[0] = getDaug(0);
467 tempDaug[1] = getDaug(1);
468 tempDaug[2] = getDaug(2);
470 tempDaug[3] = getDaug(3);
472 if (pr) std::cout <<
"decay B_MIX (1)" << std::endl;
475 p->initializePhaseSpace(getNDaug() - 1, tempDaug.data());
483 if (s1->getNDaug() > 0) { s1->deleteDaughters();}
485 EvtVector4R p1 = s1->getP4();
488 const bool changed_flavor(EvtRandom::random() > 0.5);
489 if (pr) std::cout <<
"decay B_MIX (3)" << std::endl;
490 s1->init(changed_flavor ? getDaug(getNDaug() - 1) : getDaug(0), p1);
496 if (pr) std::cout <<
"decay B_MIX (4)" << std::endl;
497 const double t1(s1->getLifetime());
500 EvtComplex osc_amp(changed_flavor ? EvtComplex(0,
501 -sin(
_freq * t1 / 2.)) : EvtComplex(cos(
_freq * t1 / 2.)));
502 if (pr) std::cout <<
"decay B_MIX (5)" << std::endl;
504 double norm = 1.0 / p1.d3mag();
505 if (pr) std::cout <<
"decay B_MIX (6)" << std::endl;
506 vertex(0, norm * osc_amp * p1 * (p->eps(0)));
507 vertex(1, norm * osc_amp * p1 * (p->eps(1)));
508 vertex(2, norm * osc_amp * p1 * (p->eps(2)));
509 if (pr) std::cout <<
"decay B_MIX (7)" << std::endl;