Belle II Software development
EvtPHSPBMix.cc
1//--------------------------------------------------------------------------
2//
3// Environment:
4// This software is part of the EvtGen package developed jointly
5// for the BaBar and CLEO collaborations. If you use all or part
6// of it, please give an appropriate acknowledgement.
7//
8// Copyright Information: See EvtGen/COPYRIGHT
9// Copyright (C) 2002 INFN-Pisa
10//
11// Module: EvtPHSPBMix.cc
12//
13// Description:
14// Routine to decay vector-> particle particle with B0 mixing, coherent B0B0-like mixing if any.
15// EvtPHSPBBMix: handle states with two neutral B
16// EvtPHSPBMix : handles states with only one neutral B
17// Phase-space kinematics, CP conservation, deltaGamma=0, p/q=1
18//
19// Based on EvtVSSBMixCPT
20//
21// Modification history:
22//
23// R. Louvot, EPFL, 2010/03/09 Module created
24// C. MacQueen, 2016/10/03 Adapted to Belle II
25//
26//------------------------------------------------------------------------
27//
28#include <cmath>
29#include <stdlib.h>
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"
38#include <string>
39#include <sstream>
40#include "EvtGenBase/EvtRandom.hh"
41
42using std::endl;
43
45
46
48
50{
51 return "PHSP_BB_MIX";
52}
53
54
55EvtDecayBase* EvtPHSPBBMix::clone()
56{
57 return new EvtPHSPBBMix;
58}
59
61{
62 // check that there we are provided exactly one parameter
63 //One arg: Delta m
64 const unsigned short int narg(2);
65 if (getNArg() != narg) {
66 EvtGenReport(EVTGEN_ERROR, "EvtGen") << "EvtPHSPBBMix generator expected "
67 << " " << narg
68 << "argument(s) (deltam, C=-1, +1 or 0 (incoherent)) but found:"
69 << getNArg() << endl;
70 EvtGenReport(EVTGEN_ERROR, "EvtGen") << "Will terminate execution!" << endl;
71 ::abort();
72 }
73
74
75 //report(DEBUG,"EvtGen") << "EvtPHSPBBMIX::init point (2)"<<std::endl;
76 //usage: mixing P0-P0b : 2 arg. -> P0 P0b PHSP_BB_MIX dm -1;
77 // : mixing P0*-P0b + cc : 4 arg. -> P0* P0b P0*b P0 PHSP_BB_MIX dm +1;
78 // : mixing P0-P0b pi : 3 arg. -> P0 P0b pi0 PHSP_BB_MIX dm -1;
79 // : mixing P0*-P0b pi : 5 arg. -> P0* P0b pi0 P0*b P0 PHSP_BB_MIX dm +1;
80 // : mixing P0-P0b pi pi : 4 arg. -> P0 P0b pi pi PHSP_BB_MIX dm -1;
81
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;
86 ::abort();
87 }
88 // checkNDaug(2,3,4,5,6);
89 // report(DEBUG,"EvtGen") << "EvtPHSPBBMIX::init point (3)"<<std::endl;
90
91
92
93 _BBpipi = false;
94
95 if (getNDaug() == 4 && (EvtPDL::chargeConj(getDaug(0)) == getDaug(1)))
96 _BBpipi = true;
97
98
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()
106 << " and "
107 << EvtPDL::name(getDaug(3)).c_str() << endl;
108 EvtGenReport(EVTGEN_ERROR, "EvtGen") << "Will terminate execution!" << endl;
109 ::abort();
110 }
111 } else {
112
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()
117 << " and "
118 << EvtPDL::name(getDaug(1)).c_str() << endl;
119 EvtGenReport(EVTGEN_ERROR, "EvtGen") << "Will terminate execution!" << endl;
120 ::abort();
121 }
122
123 }
124 //report(DEBUG,"EvtGen") << "EvtPHSPBBMIX::init point (4)"<<std::endl;
125 // check that we are asked to decay a vector particle into a pair
126 // of scalar particles
127
128 // checkSpinParent(EvtSpinType::VECTOR);
129
130 // checkSpinDaughter(0,EvtSpinType::SCALAR);
131 // checkSpinDaughter(1,EvtSpinType::SCALAR);
132
133 // check that our daughter particles are charge conjugates of each other
134
135
136 // check that both daughter particles have the same lifetime
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
140 << " Found ctau = "
141 << EvtPDL::getctau(getDaug(0)) << " mm and "
142 << EvtPDL::getctau(getDaug(1)) << " mm" << endl;
143 EvtGenReport(EVTGEN_ERROR, "EvtGen") << "Will terminate execution!" << endl;
144 ::abort();
145 }
146 // precompute quantities that will be used to generate events
147 // and print out a summary of parameters for this decay
148 //report(DEBUG,"EvtGen") << "EvtPHSPBBMIX::init point (5)"<<std::endl;
149
150 // mixing frequency in hbar/mm
151 _freq = getArg(0) / EvtConst::c;
152
153 _C = getArg(1);
154
155 // deltaG
156 //double gamma= 1/EvtPDL::getctau(getDaug(0)); // gamma/c (1/mm)
157 // q/p
158 // decay amplitudes
159 // CPT violation in mixing
160 // some printout
161 double tau = 1e12 * EvtPDL::getctau(getDaug(0)) / EvtConst::c; // in ps
162 double dm = 1e-12 * getArg(0); // B0/anti-B0 mass difference in hbar/ps
163 double x = dm * tau;
164
165 //report(DEBUG,"EvtGen") << "EvtPHSPBBMIX::init point (6)"<<std::endl;
166
167
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()) {
173 case 3:
174 case 5:
175 sss << " + " << EvtPDL::name(getDaug(2)).c_str();
176 break;
177 case 4:
178 if (_BBpipi)
179 sss << " + " << EvtPDL::name(getDaug(2)).c_str() << " + " << EvtPDL::name(getDaug(3)).c_str();
180 break;
181 default: ;
182 }
183 EvtGenReport(EVTGEN_INFO, "EvtGen") << "PHSP_BB_MIX will generate mixing :"
184 << endl << endl
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
195 << endl;
196
197
198
199
200
201}
202
204{
205 if (getNDaug() > 3) {
206 if (_BBpipi)
207 return 4;
208 else
209 return getNDaug() - 2;
210 }
211 //else
212 return getNDaug();
213
214
215}
216
217
219{
220 // this value is ok for reasonable values of all the parameters
221 setProbMax(4.0);
222 //noProbMax();
223}
224
225void EvtPHSPBBMix::prlp(int i) const
226{
227 EvtGenReport(EVTGEN_INFO, "EvtGen") << "decay p" << i << endl;
228}
229
230void EvtPHSPBBMix::decay(EvtParticle* p)
231{
232 //p->initializePhaseSpace(getNDaug(),getDaugs());
233 //return;
234
235 if (_print_info) prlp(1);
236
237 static const EvtId B0(EvtPDL::getId("B0"));
238 static const EvtId B0B(EvtPDL::getId("anti-B0"));
239
240 // generate a final state according to phase space
241
242 // const bool ADaug(getNDaug()==4);//true if aliased daughters
243 if (_print_info) prlp(102);
244
245 const bool NCC(getNDaug() >= 4 && !_BBpipi); //true if daughters not charged-conjugated
246 // const bool BBpi(getNDaug()==
247 if (_print_info) prlp(103);
248 if (NCC) {
249 if (_print_info) prlp(2);
250 std::vector<EvtId> tempDaug(getNDaug() - 2);
251 tempDaug[0] = getDaug(0);
252 tempDaug[1] = getDaug(1);
253 if (getNDaug() == 5)
254 tempDaug[2] = getDaug(2);
255
256 p->initializePhaseSpace(getNDaug() - 2, tempDaug.data());
257
258 if (_print_info) prlp(222);
259 } else { //nominal case.
260 p->initializePhaseSpace(getNDaug(), getDaugs());
261
262 }
263 if (_print_info) prlp(3);
264 EvtParticle* s1, *s2;
265
266 s1 = p->getDaug(0);
267 s2 = p->getDaug(1);
268 //delete any daughters - if there are daughters, they
269 //are from the initialization and will be redone later
270 if (s1->getNDaug() > 0) { s1->deleteDaughters();}
271 if (s2->getNDaug() > 0) { s2->deleteDaughters();}
272
273 EvtVector4R p1 = s1->getP4();
274 EvtVector4R p2 = s2->getP4();
275
276 // throw 2 random numbers to decide whether B1 and B2 are B0 or B0B
277 const EvtId B1(EvtRandom::random() > 0.5 ? B0 : B0B);
278 const EvtId B2(EvtRandom::random() > 0.5 ? B0 : B0B);
279 if (_print_info) prlp(5);
280
281 if (NCC)
282
283 {
284 if (getNDaug() == 4) {
285 s1->init(B1 == B0 ? getDaug(0) : getDaug(2), p1);
286 s2->init(B2 == B0 ? getDaug(3) : getDaug(1), p2);
287 }
288 if (getNDaug() == 5) {
289 s1->init(B1 == B0 ? getDaug(0) : getDaug(3), p1);
290 s2->init(B2 == B0 ? getDaug(4) : getDaug(1), p2);
291 }
292 } else {
293 s1->init(B1 == B0 ? getDaug(0) : getDaug(1), p1);
294 s2->init(B2 == B0 ? getDaug(0) : getDaug(1), p2);
295 }
296
297 if (_print_info) prlp(6);
298
299 // choose a decay time for each final state particle using the
300 // lifetime (which must be the same for both particles) in pdt.table
301 // and calculate the lifetime difference for this event
302 s1->setLifetime();
303 s2->setLifetime();
304
305 const double t1(s1->getLifetime());
306 const double t2(s2->getLifetime());
307
308 // calculate the oscillation amplitude, based on whether this event is mixed or not
309 EvtComplex osc_amp(Amplitude(t1, t2, B1 == B0, B2 == B0)); //Amplitude return <B0(B)B0(B)|B1B2>
310 if (_print_info) prlp(8);
311
312
313 // store the amplitudes for each parent spin basis state
314 double norm = 1.0 / p1.d3mag();
315 if (_print_info) prlp(9);
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)));
319 if (_print_info) prlp(10);
320 return ;
321}
322
323EvtComplex EvtPHSPBBMix::Amplitude(const double& t1, const double& t2, bool B1_is_B0, bool B2_is_B0) const
324{
325 if (_C != 0 && _C != -1 && _C != 1)
326 return EvtComplex(0., 0.);
327
328 const double f(_freq);
329 if (B1_is_B0 && !B2_is_B0) {
330 if (_C == 0)
331 return EvtComplex(cos(f * t1 / 2.) * cos(f * t2 / 2.), 0.);
332 else
333 return EvtComplex(cos(f * (t2 + _C * t1) / 2.) / sqrt(2.), 0.);
334 }
335 if (!B1_is_B0 && B2_is_B0) {
336 if (_C == 0)
337 return EvtComplex(-sin(f * t1 / 2.) * sin(f * t2 / 2.), 0.);
338 else
339 return EvtComplex(cos(f * (t2 + _C * t1) / 2.) * _C / sqrt(2.), 0.);
340 }
341 if (B1_is_B0 && B2_is_B0) {
342 if (_C == 0)
343 return EvtComplex(0., -cos(f * t1 / 2.) * sin(f * t2 / 2.));
344 else
345 return EvtComplex(0., -sin(f * (t2 + _C * t1) / 2.) / sqrt(2.));
346 }
347 if (!B1_is_B0 && !B2_is_B0) {
348 if (_C == 0)
349 return EvtComplex(0., -sin(f * t1 / 2.) * cos(f * t2 / 2.));
350 else
351 return EvtComplex(0., -sin(f * (t2 + _C * t1) / 2.) * _C / sqrt(2.));
352 }
353 // no way to reach this but compiler complains without a return statement
354 return EvtComplex(0., 0.);
355}
356
357
360
361
363
364
366
368{
369 return "PHSP_B_MIX";
370}
371
372
373EvtDecayBase* EvtPHSPBMix::clone()
374{
375 return new EvtPHSPBMix;
376}
377
379{
380 // check that there we are provided exactly one parameter
381 //One arg: Delta m
382 const unsigned short int narg(1);
383 if (getNArg() != narg) {
384 EvtGenReport(EVTGEN_ERROR, "EvtGen") << "EvtPHSPBMix generator expected "
385 << " " << narg
386 << " argument(s) (delta m) but found:"
387 << getNArg() << endl;
388 EvtGenReport(EVTGEN_ERROR, "EvtGen") << "Will terminate execution!" << endl;
389 ::abort();
390 }
391
392
393 //report(DEBUG,"EvtGen") << "EvtPHSPBMIX::init point (2)"<<std::endl;
394 //usage: mixing P0 P- pi+ : 4 arg. -> P0 P- pi+ P0b PHSP_B_MIX dm;
395 // : mixing P0 P- pi+ pi : 5 arg. -> P0 P- pi+ pi0 P0b PHSP_B_MIX dm;
396
397
398 checkNDaug(4, 5);
399 // report(DEBUG,"EvtGen") << "EvtPHSPBMIX::init point (3)"<<std::endl;
400
401 bool BBpipi(getNDaug() == 5);
402
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;
409 ::abort();
410 }
411 //parameters for mixing of particle1
412
413 _freq = getArg(0) / EvtConst::c;
414
415 //double gamma= 1/EvtPDL::getctau(getDaug(0)); // gamma/c (1/mm)
416
417 double tau = 1e12 * EvtPDL::getctau(getDaug(0)) / EvtConst::c; // in ps
418 double dm = 1e-12 * getArg(0); // B0/anti-B0 mass difference in hbar/ps
419 double x = dm * tau;
420
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();
426 if (BBpipi)
427 sss << " + " << EvtPDL::name(getDaug(3)).c_str();
428
429 EvtGenReport(EVTGEN_INFO, "EvtGen") << "PHSP_B_MIX will generate mixing :"
430 << endl << endl
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
440 << endl;
441
442}
443
445{
446 return getNDaug() - 1;
447}
448
450{
451 // this value is ok for reasonable values of all the parameters
452 setProbMax(4.0);
453 //noProbMax();
454}
455
456void EvtPHSPBMix::decay(EvtParticle* p)
457{
458
459 // generate a final state according to phase space
460
461 // const bool BBpipi(getNDaug()!=4);//true if BBpipi, if not BBpi
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);
467 if (getNDaug() == 5)
468 tempDaug[3] = getDaug(3);
469
470 if (_print_info) std::cout << "decay B_MIX (1)" << std::endl;
471
472
473 p->initializePhaseSpace(getNDaug() - 1, tempDaug.data());
474
475 EvtParticle* s1;
476
477 s1 = p->getDaug(0);
478
479 //delete any daughters - if there are daughters, they
480 //are from the initialization and will be redone later
481 if (s1->getNDaug() > 0) { s1->deleteDaughters();}
482
483 EvtVector4R p1 = s1->getP4();
484
485 // throw 1 random numbers to decide whether B1 is B0 or B0B
486 const bool changed_flavor(EvtRandom::random() > 0.5); //Daug(0) becomes Daug(1)
487 if (_print_info) std::cout << "decay B_MIX (3)" << std::endl;
488 s1->init(changed_flavor ? getDaug(getNDaug() - 1) : getDaug(0), p1);
489
490 // choose a decay time for each final state particle using the
491 // lifetime (which must be the same for both particles) in pdt.table
492 // and calculate the lifetime difference for this event
493 s1->setLifetime();
494 if (_print_info) std::cout << "decay B_MIX (4)" << std::endl;
495 const double t1(s1->getLifetime());
496
497 // calculate the oscillation amplitude, based on whether this event is mixed or not
498 EvtComplex osc_amp(changed_flavor ? EvtComplex(0,
499 -sin(_freq * t1 / 2.)) : EvtComplex(cos(_freq * t1 / 2.))); //Amplitude return <B0(B)|B1>
500 if (_print_info) std::cout << "decay B_MIX (5)" << std::endl;
501 // store the amplitudes for each parent spin basis state
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;
508 return ;
509}
510
511
512
The class provides routine to decay vector-> particle particle with B0 mixing, coherent B0B0-like mix...
Definition: EvtPHSPBMix.h:36
void init()
Init function.
Definition: EvtPHSPBMix.cc:60
void prlp(int) const
Number of real daughters.
Definition: EvtPHSPBMix.cc:225
bool _print_info
Print evtgeninfo?
Definition: EvtPHSPBMix.h:82
virtual ~EvtPHSPBBMix()
Default destructor.
Definition: EvtPHSPBMix.cc:47
EvtPHSPBBMix()
Default constructor.
Definition: EvtPHSPBMix.h:41
EvtDecayBase * clone()
Clone the decay
Definition: EvtPHSPBMix.cc:55
void initProbMax()
Init maximal prob.
Definition: EvtPHSPBMix.cc:218
int nRealDaughters()
Number of real daughters.
Definition: EvtPHSPBMix.cc:203
double _C
C eigenvalue, 0= incoherent.
Definition: EvtPHSPBMix.h:76
double _freq
mixing frequency in hbar/mm
Definition: EvtPHSPBMix.h:73
EvtComplex Amplitude(const double &t1, const double &t2, bool B1_is_B0, bool B2_is_B0) const
Calculate amplitude.
Definition: EvtPHSPBMix.cc:323
std::string getName()
Get function Name
Definition: EvtPHSPBMix.cc:49
bool _BBpipi
Is BBpipi?
Definition: EvtPHSPBMix.h:79
void decay(EvtParticle *p)
Decay function.
Definition: EvtPHSPBMix.cc:230
The class provides routine to decay vector-> particle particle with B0 mixing, handles states with on...
Definition: EvtPHSPBMix.h:88
void init()
Init function.
Definition: EvtPHSPBMix.cc:378
bool _print_info
Print evtgeninfo?
Definition: EvtPHSPBMix.h:122
EvtDecayBase * clone()
Clone the decay
Definition: EvtPHSPBMix.cc:373
virtual ~EvtPHSPBMix()
Default destructor.
Definition: EvtPHSPBMix.cc:365
void initProbMax()
Init maximal prob.
Definition: EvtPHSPBMix.cc:449
int nRealDaughters()
Number of real daughters.
Definition: EvtPHSPBMix.cc:444
EvtPHSPBMix()
Default constructor.
Definition: EvtPHSPBMix.h:93
double _freq
mixing frequency in hbar/mm
Definition: EvtPHSPBMix.h:119
std::string getName()
Get function Name
Definition: EvtPHSPBMix.cc:367
void decay(EvtParticle *p)
Decay function.
Definition: EvtPHSPBMix.cc:456
#define B2_EVTGEN_REGISTER_MODEL(classname)
Class to register B2_EVTGEN_REGISTER_MODEL.