evtgen is hosted by Hepforge, IPPP Durham
EvtGen  2.0.0
Monte Carlo generator of particle decays, in particular the weak decays of heavy flavour particles such as B mesons.
EvtD0mixDalitz.cpp
Go to the documentation of this file.
1 
2 /***********************************************************************
3 * Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
4 * *
5 * This file is part of EvtGen. *
6 * *
7 * EvtGen is free software: you can redistribute it and/or modify *
8 * it under the terms of the GNU General Public License as published by *
9 * the Free Software Foundation, either version 3 of the License, or *
10 * (at your option) any later version. *
11 * *
12 * EvtGen is distributed in the hope that it will be useful, *
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15 * GNU General Public License for more details. *
16 * *
17 * You should have received a copy of the GNU General Public License *
18 * along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
19 ***********************************************************************/
20 
22 
25 #include "EvtGenBase/EvtPDL.hh"
27 #include "EvtGenBase/EvtPatches.hh"
28 #include "EvtGenBase/EvtRandom.hh"
30 
31 #include <cmath> // for std::fabs
32 
33 // Initialize the static variables.
37 
41 
45 
49 
51 {
52  // check that there are 0 arguments
53  checkNDaug( 3 );
54 
55  if ( getNArg() ) {
56  if ( getNArg() == 2 ) {
57  _x = getArg( 0 );
58  _y = getArg( 1 );
59  } else if ( getNArg() == 4 ) {
60  _x = getArg( 0 );
61  _y = getArg( 1 );
62  _qp = EvtComplex( getArg( 2 ), getArg( 3 ) );
63  } else if ( getNArg() == 5 ) {
64  _x = getArg( 0 );
65  _y = getArg( 1 );
66  _qp = EvtComplex( getArg( 2 ), getArg( 3 ) );
68  4 ); // RBW by default. If arg4 is set, do K-matrix.
69  } else {
70  EvtGenReport( EVTGEN_ERROR, "EvtD0mixDalitz" )
71  << "Number of arguments for this model must be 0, 2, 4 or 5:"
72  << std::endl
73  << "[ x y ][ qp.re qp.im ][ doK-matrix ]" << std::endl
74  << "Check your dec file." << std::endl;
75  exit( 1 );
76  }
77  }
78 
83 
84  readPDGValues();
85 
86  // Get the EvtId of the D0 and its (3) daughters.
87  EvtId parId = getParentId();
88 
89  EvtId dau[3];
90  for ( int index = 0; index < 3; index++ )
91  dau[index] = getDaug( index );
92 
93  if ( parId == _D0 ) // Look for K0bar h+ h-. The order must be K[0SL] h+ h-
94  for ( int index = 0; index < 3; index++ )
95  if ( ( dau[index] == _K0B ) || ( dau[index] == _KS ) ||
96  ( dau[index] == _KL ) )
97  _d1 = index;
98  else if ( ( dau[index] == _PIP ) || ( dau[index] == _KP ) )
99  _d2 = index;
100  else if ( ( dau[index] == _PIM ) || ( dau[index] == _KM ) )
101  _d3 = index;
102  else
104  else if ( parId == _D0B ) // Look for K0 h+ h-. The order must be K[0SL] h- h+
105  for ( int index = 0; index < 3; index++ )
106  if ( ( dau[index] == _K0 ) || ( dau[index] == _KS ) ||
107  ( dau[index] == _KL ) )
108  _d1 = index;
109  else if ( ( dau[index] == _PIM ) || ( dau[index] == _KM ) )
110  _d2 = index;
111  else if ( ( dau[index] == _PIP ) || ( dau[index] == _KP ) )
112  _d3 = index;
113  else
115  else
117 
118  // If the D meson is a D0bar, the expressions should use p/q instead of q/p.
119  if ( parId == _D0B )
120  _qp = 1.0 / _qp;
121 
122  // At this point, if parId is D0bar, the amplitude is the D0bar amplitude, the conjugated amplitude
123  // is the amplitude of the D0 decay, and _qp means p/q, so it is like changing the meaning of
124  // A <-> Abar, and p <-> q. It is just a trick so after this point the code for D0bar can be the
125  // same as the code for D0.
126 
127  // Check if we're dealing with Ks pi pi or with Ks K K.
128  _isKsPiPi = false;
129  if ( dau[_d2] == _PIP || dau[_d2] == _PIM )
130  _isKsPiPi = true;
131 }
132 
134 {
135  // Same structure for all of these decays.
136  part->initializePhaseSpace( getNDaug(), getDaugs() );
137  EvtVector4R pA = part->getDaug( _d1 )->getP4();
138  EvtVector4R pB = part->getDaug( _d2 )->getP4();
139  EvtVector4R pC = part->getDaug( _d3 )->getP4();
140 
141  // Squared invariant masses.
142  double m2AB = ( pA + pB ).mass2();
143  double m2AC = ( pA + pC ).mass2();
144  double m2BC = ( pB + pC ).mass2();
145 
146  // Dalitz amplitudes of the decay of the particle and that of the antiparticle.
147  EvtComplex ampDalitz;
148  EvtComplex ampAntiDalitz;
149 
150  if ( _isKsPiPi ) { // For Ks pi pi
151  EvtDalitzPoint point( _mKs, _mPi, _mPi, m2AB, m2BC, m2AC );
152  EvtDalitzPoint antiPoint( _mKs, _mPi, _mPi, m2AC, m2BC, m2AB );
153 
154  ampDalitz = dalitzKsPiPi( point );
155  ampAntiDalitz = dalitzKsPiPi( antiPoint );
156  } else { // For Ks K K
157  EvtDalitzPoint point( _mKs, _mK, _mK, m2AB, m2BC, m2AC );
158  EvtDalitzPoint antiPoint( _mKs, _mK, _mK, m2AC, m2BC, m2AB );
159 
160  ampDalitz = dalitzKsKK( point );
161  ampAntiDalitz = dalitzKsKK( antiPoint );
162  }
163 
164  // Assume there's no direct CP violation.
165  EvtComplex barAOverA = ampAntiDalitz / ampDalitz;
166 
167  // CP violation in the interference. _qp implements CP violation in the mixing.
168  EvtComplex chi = _qp * barAOverA;
169 
170  // Generate a negative exponential life time. p( gt ) = ( 1 - y ) * e^{ - ( 1 - y ) gt }
171  double gt = -log( EvtRandom::Flat() ) / ( 1.0 - std::fabs( _y ) );
172  part->setLifetime( gt / _gamma );
173 
174  // Compute time dependent amplitude.
175  EvtComplex amp = 0.5 * ampDalitz * exp( -std::fabs( _y ) * gt / 2.0 ) *
176  ( ( 1.0 + chi ) * h1( gt ) + ( 1.0 - chi ) * h2( gt ) );
177 
178  vertex( amp );
179 
180  return;
181 }
182 
184 {
185  // Define the EvtIds.
186  _D0 = EvtPDL::getId( "D0" );
187  _D0B = EvtPDL::getId( "anti-D0" );
188  _KM = EvtPDL::getId( "K-" );
189  _KP = EvtPDL::getId( "K+" );
190  _K0 = EvtPDL::getId( "K0" );
191  _K0B = EvtPDL::getId( "anti-K0" );
192  _KL = EvtPDL::getId( "K_L0" );
193  _KS = EvtPDL::getId( "K_S0" );
194  _PIM = EvtPDL::getId( "pi-" );
195  _PIP = EvtPDL::getId( "pi+" );
196 
197  // Read the relevant masses.
198  _mD0 = EvtPDL::getMass( _D0 );
199  _mKs = EvtPDL::getMass( _KS );
201  _mK = EvtPDL::getMass( _KP );
202 
203  // Compute the decay rate from the parameter in the evt.pdl file.
204  _ctau = EvtPDL::getctau( EvtPDL::getId( "D0" ) );
205 
206  _gamma = 1.0 / _ctau; // ALERT: Gamma is not 1 / tau.
207 }
208 
210 {
211  static const EvtDalitzPlot plot( _mKs, _mPi, _mPi, _mD0 );
212 
213  EvtComplex amp = 0.;
214 
215  if ( _isRBWmodel ) {
216  // This corresponds to relativistic Breit-Wigner distributions. Not K-matrix.
217  // Defining resonances.
218  static EvtDalitzReso KStarm( plot, _BC, _AC, _VECTOR, 0.893606,
219  0.0463407, _RBW );
220  static EvtDalitzReso KStarp( plot, _BC, _AB, _VECTOR, 0.893606,
221  0.0463407, _RBW );
222  static EvtDalitzReso rho0( plot, _AC, _BC, _VECTOR, 0.7758, 0.1464, _GS );
223  static EvtDalitzReso omega( plot, _AC, _BC, _VECTOR, 0.78259, 0.00849,
224  _RBW );
225  static EvtDalitzReso f0_980( plot, _AC, _BC, _SCALAR, 0.975, 0.044, _RBW );
226  static EvtDalitzReso f0_1370( plot, _AC, _BC, _SCALAR, 1.434, 0.173,
227  _RBW );
228  static EvtDalitzReso f2_1270( plot, _AC, _BC, _TENSOR, 1.2754, 0.1851,
229  _RBW );
230  static EvtDalitzReso K0Starm_1430( plot, _BC, _AC, _SCALAR, 1.459,
231  0.175, _RBW );
232  static EvtDalitzReso K0Starp_1430( plot, _BC, _AB, _SCALAR, 1.459,
233  0.175, _RBW );
234  static EvtDalitzReso K2Starm_1430( plot, _BC, _AC, _TENSOR, 1.4256,
235  0.0985, _RBW );
236  static EvtDalitzReso K2Starp_1430( plot, _BC, _AB, _TENSOR, 1.4256,
237  0.0985, _RBW );
238  static EvtDalitzReso sigma( plot, _AC, _BC, _SCALAR, 0.527699, 0.511861,
239  _RBW );
240  static EvtDalitzReso sigma2( plot, _AC, _BC, _SCALAR, 1.03327,
241  0.0987890, _RBW );
242  static EvtDalitzReso KStarm_1680( plot, _BC, _AC, _VECTOR, 1.677, 0.205,
243  _RBW );
244 
245  // Adding terms to the amplitude with their corresponding amplitude and phase terms.
246  amp += EvtComplex( 0.848984, 0.893618 );
247  amp += EvtComplex( -1.16356, 1.19933 ) * KStarm.evaluate( point );
248  amp += EvtComplex( 0.106051, -0.118513 ) * KStarp.evaluate( point );
249  amp += EvtComplex( 1.0, 0.0 ) * rho0.evaluate( point );
250  amp += EvtComplex( -0.0249569, 0.0388072 ) * omega.evaluate( point );
251  amp += EvtComplex( -0.423586, -0.236099 ) * f0_980.evaluate( point );
252  amp += EvtComplex( -2.16486, 3.62385 ) * f0_1370.evaluate( point );
253  amp += EvtComplex( 0.217748, -0.133327 ) * f2_1270.evaluate( point );
254  amp += EvtComplex( 1.62128, 1.06816 ) * K0Starm_1430.evaluate( point );
255  amp += EvtComplex( 0.148802, 0.0897144 ) * K0Starp_1430.evaluate( point );
256  amp += EvtComplex( 1.15489, -0.773363 ) * K2Starm_1430.evaluate( point );
257  amp += EvtComplex( 0.140865, -0.165378 ) * K2Starp_1430.evaluate( point );
258  amp += EvtComplex( -1.55556, -0.931685 ) * sigma.evaluate( point );
259  amp += EvtComplex( -0.273791, -0.0535596 ) * sigma2.evaluate( point );
260  amp += EvtComplex( -1.69720, 0.128038 ) * KStarm_1680.evaluate( point );
261  } else {
262  // This corresponds to the complete model (RBW, GS, LASS and K-matrix).
263  // Defining resonances.
264  static EvtDalitzReso KStarm( plot, _BC, _AC, _VECTOR, 0.893619,
265  0.0466508, _RBW );
266  static EvtDalitzReso KStarp( plot, _BC, _AB, _VECTOR, 0.893619,
267  0.0466508, _RBW );
268  static EvtDalitzReso rho0( plot, _AC, _BC, _VECTOR, 0.7758, 0.1464, _GS );
269  static EvtDalitzReso omega( plot, _AC, _BC, _VECTOR, 0.78259, 0.00849,
270  _RBW );
271  static EvtDalitzReso f2_1270( plot, _AC, _BC, _TENSOR, 1.2754, 0.1851,
272  _RBW );
273  static EvtDalitzReso K0Starm_1430( plot, _AC, 1.46312, 0.232393, 1.0746,
274  -1.83214, .803516, 2.32788, 1.0,
275  -5.31306 ); // LASS
276  static EvtDalitzReso K0Starp_1430( plot, _AB, 1.46312, 0.232393, 1.0746,
277  -1.83214, .803516, 2.32788, 1.0,
278  -5.31306 ); // LASS
279  static EvtDalitzReso K2Starm_1430( plot, _BC, _AC, _TENSOR, 1.4256,
280  0.0985, _RBW );
281  static EvtDalitzReso K2Starp_1430( plot, _BC, _AB, _TENSOR, 1.4256,
282  0.0985, _RBW );
283  static EvtDalitzReso KStarm_1680( plot, _BC, _AC, _VECTOR, 1.677, 0.205,
284  _RBW );
285 
286  // Defining K-matrix.
287  static EvtComplex fr12( 1.87981, -0.628378 );
288  static EvtComplex fr13( 4.3242, 2.75019 );
289  static EvtComplex fr14( 3.22336, 0.271048 );
290  static EvtComplex fr15( 0.0, 0.0 );
291  static EvtDalitzReso Pole1( plot, _BC, "Pole1", _KMAT, fr12, fr13, fr14,
292  fr15, -0.0694725 );
293  static EvtDalitzReso Pole2( plot, _BC, "Pole2", _KMAT, fr12, fr13, fr14,
294  fr15, -0.0694725 );
295  static EvtDalitzReso Pole3( plot, _BC, "Pole3", _KMAT, fr12, fr13, fr14,
296  fr15, -0.0694725 );
297  static EvtDalitzReso Pole4( plot, _BC, "Pole4", _KMAT, fr12, fr13, fr14,
298  fr15, -0.0694725 );
299  static EvtDalitzReso kmatrix( plot, _BC, "f11prod", _KMAT, fr12, fr13,
300  fr14, fr15, -0.0694725 );
301 
302  // Adding terms to the amplitude with their corresponding amplitude and phase terms.
303  amp += EvtComplex( -1.31394, 1.14072 ) * KStarm.evaluate( point );
304  amp += EvtComplex( 0.116239, -0.107287 ) * KStarp.evaluate( point );
305  amp += EvtComplex( 1.0, 0.0 ) * rho0.evaluate( point );
306  amp += EvtComplex( -0.0313343, 0.0424013 ) * omega.evaluate( point );
307  amp += EvtComplex( 0.559412, -0.232336 ) * f2_1270.evaluate( point );
308  amp += EvtComplex( 7.35400, -3.67637 ) * K0Starm_1430.evaluate( point );
309  amp += EvtComplex( 0.255913, -0.190459 ) * K0Starp_1430.evaluate( point );
310  amp += EvtComplex( 1.05397, -0.936297 ) * K2Starm_1430.evaluate( point );
311  amp += EvtComplex( -0.00760136, -0.0908624 ) *
312  K2Starp_1430.evaluate( point );
313  amp += EvtComplex( -1.45336, -0.164494 ) * KStarm_1680.evaluate( point );
314  amp += EvtComplex( -1.81830, 9.10680 ) * Pole1.evaluate( point );
315  amp += EvtComplex( 10.1751, 3.87961 ) * Pole2.evaluate( point );
316  amp += EvtComplex( 23.6569, -4.94551 ) * Pole3.evaluate( point );
317  amp += EvtComplex( 0.0725431, -9.16264 ) * Pole4.evaluate( point );
318  amp += EvtComplex( -2.19449, -7.62666 ) * kmatrix.evaluate( point );
319 
320  amp *= 0.97; // Multiply by a constant in order to use the same maximum as RBW model.
321  }
322 
323  return amp;
324 }
325 
327 {
328  static const EvtDalitzPlot plot( _mKs, _mK, _mK, _mD0 );
329 
330  // Defining resonances.
331  static EvtDalitzReso a00_980( plot, _AC, _BC, _SCALAR, 0.999, _RBW,
332  0.550173, 0.324, _EtaPic );
333  static EvtDalitzReso phi( plot, _AC, _BC, _VECTOR, 1.01943, 0.00459319, _RBW );
334  static EvtDalitzReso a0p_980( plot, _AC, _AB, _SCALAR, 0.999, _RBW,
335  0.550173, 0.324, _EtaPic );
336  static EvtDalitzReso f0_1370( plot, _AC, _BC, _SCALAR, 1.350, 0.265, _RBW );
337  static EvtDalitzReso a0m_980( plot, _AB, _AC, _SCALAR, 0.999, _RBW,
338  0.550173, 0.324, _EtaPic );
339  static EvtDalitzReso f0_980( plot, _AC, _BC, _SCALAR, 0.965, _RBW, 0.695,
340  0.165, _PicPicKK );
341  static EvtDalitzReso f2_1270( plot, _AC, _BC, _TENSOR, 1.2754, 0.1851, _RBW );
342  static EvtDalitzReso a00_1450( plot, _AC, _BC, _SCALAR, 1.474, 0.265, _RBW );
343  static EvtDalitzReso a0p_1450( plot, _AC, _AB, _SCALAR, 1.474, 0.265, _RBW );
344  static EvtDalitzReso a0m_1450( plot, _AB, _AC, _SCALAR, 1.474, 0.265, _RBW );
345 
346  // Adding terms to the amplitude with their corresponding amplitude and phase terms.
347  EvtComplex amp( 0., 0. ); // Phase space amplitude.
348  amp += EvtComplex( 1.0, 0.0 ) * a00_980.evaluate( point );
349  amp += EvtComplex( -0.126314, 0.188701 ) * phi.evaluate( point );
350  amp += EvtComplex( -0.561428, 0.0135338 ) * a0p_980.evaluate( point );
351  amp += EvtComplex( 0.035, -0.00110488 ) * f0_1370.evaluate( point );
352  amp += EvtComplex( -0.0872735, 0.0791190 ) * a0m_980.evaluate( point );
353  amp += EvtComplex( 0.0, 0.0 ) * f0_980.evaluate( point );
354  amp += EvtComplex( 0.257341, -0.0408343 ) * f2_1270.evaluate( point );
355  amp += EvtComplex( -0.0614342, -0.649930 ) * a00_1450.evaluate( point );
356  amp += EvtComplex( -0.104629, 0.830120 ) * a0p_1450.evaluate( point );
357  amp += EvtComplex( 0.0, 0.0 ) * a0m_1450.evaluate( point );
358 
359  return 2.8 *
360  amp; // Multiply by 2.8 in order to reuse the same probmax as Ks pi pi.
361 }
362 
363 // < f | H | D^0 (t) > = 1/2 * [ ( 1 + \chi_f ) * A_f * e_1(gt) + ( 1 - \chi_f ) * A_f * e_2(gt) ]
364 // < f | H | D^0 (t) > = 1/2 * exp( -gamma t / 2 ) * [ ( 1 + \chi_f ) * A_f * h_1(t) + ( 1 - \chi_f ) * A_f * h_2(t) ]
365 // e{1,2}( gt ) = exp( -gt / 2 ) * h{1,2}( gt ).
366 EvtComplex EvtD0mixDalitz::h1( const double& gt ) const
367 {
368  return exp( -EvtComplex( _y, _x ) * gt / 2. );
369 }
370 
371 EvtComplex EvtD0mixDalitz::h2( const double& gt ) const
372 {
373  return exp( EvtComplex( _y, _x ) * gt / 2. );
374 }
static const EvtCyclic3::Pair & _BC
double getArg(unsigned int j)
EvtComplex evaluate(const EvtDalitzPoint &p)
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=0)
Definition: EvtReport.cpp:33
static const EvtDalitzReso::NumType & _GS
static const EvtSpinType::spintype & _SCALAR
EvtComplex dalitzKsKK(const EvtDalitzPoint &point)
static double getctau(EvtId i)
Definition: EvtPDL.cpp:357
static const EvtSpinType::spintype & _VECTOR
void decay(EvtParticle *p) override
void init() override
void setLifetime(double tau)
EvtId * getDaugs()
Definition: EvtDecayBase.hh:66
EvtComplex h1(const double &ct) const
Definition: EvtId.hh:27
static const EvtDalitzReso::CouplingType & _EtaPic
EvtId getParentId() const
Definition: EvtDecayBase.hh:61
void vertex(const EvtComplex &amp)
Definition: EvtDecayAmp.hh:37
double initializePhaseSpace(unsigned int numdaughter, EvtId *daughters, bool forceResetMasses=false, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
EvtComplex h2(const double &ct) const
void checkNDaug(int d1, int d2=-1)
static const EvtCyclic3::Pair & _AC
void checkSpinParent(EvtSpinType::spintype sp)
static double Flat()
Definition: EvtRandom.cpp:72
static const EvtDalitzReso::CouplingType & _PicPicKK
static EvtId getId(const std::string &name)
Definition: EvtPDL.cpp:287
const EvtVector4R & getP4() const
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
int getNDaug() const
Definition: EvtDecayBase.hh:65
static const EvtDalitzReso::NumType & _RBW
EvtComplex exp(const EvtComplex &c)
Definition: EvtComplex.hh:240
EvtComplex dalitzKsPiPi(const EvtDalitzPoint &point)
EvtParticle * getDaug(int i)
Definition: EvtParticle.cpp:91
EvtComplex _qp
void reportInvalidAndExit() const
static const EvtCyclic3::Pair & _AB
static const EvtSpinType::spintype & _TENSOR
static double getMass(EvtId i)
Definition: EvtPDL.cpp:319
static const EvtDalitzReso::NumType & _KMAT
int getNArg() const
Definition: EvtDecayBase.hh:68
EvtId getDaug(int i) const
Definition: EvtDecayBase.hh:67