Belle II Software  release-08-01-10
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 
42 using std::endl;
43 
45 
46 
48 
49 std::string EvtPHSPBBMix::getName()
50 {
51  return "PHSP_BB_MIX";
52 }
53 
54 
55 EvtDecayBase* 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 
225 void EvtPHSPBBMix::prlp(int i) const
226 {
227  EvtGenReport(EVTGEN_INFO, "EvtGen") << "decay p" << i << endl;
228 }
229 
230 void 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 
323 EvtComplex 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 
367 std::string EvtPHSPBMix::getName()
368 {
369  return "PHSP_B_MIX";
370 }
371 
372 
373 EvtDecayBase* 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 
456 void 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
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
#define B2_EVTGEN_REGISTER_MODEL(classname)
Class to register B2_EVTGEN_REGISTER_MODEL.