BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtD0mixDalitz.cc
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: BaBar detector at the SLAC PEP-II B-factory
3 * Package: EvtGenModels
4 * File: $Id: EvtD0mixDalitz.cc,v 1.1 2009/05/08 01:59:56 pingrg Exp $
5 *
6 * Description:
7 * The D0mixDalitz model, with many resonances and mixing implemented.
8 *
9 * Modification history:
10 * Jordi Garra Tic� 2008/07/03 File created
11 *****************************************************************************/
12
13#include "EvtD0mixDalitz.hh"
21
22// Initialize the static variables.
23const EvtSpinType::spintype& EvtD0mixDalitz::_SCALAR = EvtSpinType::SCALAR;
24const EvtSpinType::spintype& EvtD0mixDalitz::_VECTOR = EvtSpinType::VECTOR;
25const EvtSpinType::spintype& EvtD0mixDalitz::_TENSOR = EvtSpinType::TENSOR;
26
27const EvtDalitzReso::CouplingType& EvtD0mixDalitz::_EtaPic = EvtDalitzReso::EtaPic;
28const EvtDalitzReso::CouplingType& EvtD0mixDalitz::_PicPicKK = EvtDalitzReso::PicPicKK;
29
30const EvtDalitzReso::NumType& EvtD0mixDalitz::_RBW = EvtDalitzReso::RBW_CLEO_ZEMACH;
31const EvtDalitzReso::NumType& EvtD0mixDalitz::_GS = EvtDalitzReso::GS_CLEO_ZEMACH;
32const EvtDalitzReso::NumType& EvtD0mixDalitz::_KMAT = EvtDalitzReso::K_MATRIX;
33
34const EvtCyclic3::Pair& EvtD0mixDalitz::_AB = EvtCyclic3::AB;
35const EvtCyclic3::Pair& EvtD0mixDalitz::_AC = EvtCyclic3::AC;
36const EvtCyclic3::Pair& EvtD0mixDalitz::_BC = EvtCyclic3::BC;
37
39 // check that there are 0 arguments
40 checkNDaug( 3 );
41
42 if ( getNArg() )
43 if ( getNArg() == 2 )
44 {
45 _x = getArg( 0 );
46 _y = getArg( 1 );
47 }
48 else if ( getNArg() == 4 )
49 {
50 _x = getArg( 0 );
51 _y = getArg( 1 );
52 _qp = EvtComplex( getArg( 2 ), getArg( 3 ) );
53 }
54 else if ( getNArg() == 5 )
55 {
56 _x = getArg( 0 );
57 _y = getArg( 1 );
58 _qp = EvtComplex( getArg( 2 ), getArg( 3 ) );
59 _isRBWmodel = !getArg( 4 ); // RBW by default. If arg4 is set, do K-matrix.
60 }
61 else
62 {
63 report( ERROR, "EvtD0mixDalitz" )
64 << "Number of arguments for this model must be 0, 2, 4 or 5:" << std::endl
65 << "[ x y ][ qp.re qp.im ][ doK-matrix ]" << std::endl
66 << "Check your dec file." << std::endl;
67 exit( 1 );
68 }
69
70 checkSpinParent( _SCALAR );
71 checkSpinDaughter( 0, _SCALAR );
72 checkSpinDaughter( 1, _SCALAR );
73 checkSpinDaughter( 2, _SCALAR );
74
75 readPDGValues();
76
77 // Get the EvtId of the D0 and its (3) daughters.
78 EvtId parId = getParentId();
79
80 EvtId dau[3];
81 for ( int index = 0; index < 3; index++ ) dau[index] = getDaug( index );
82
83 if ( parId == _D0 ) // Look for K0bar h+ h-. The order must be K[0SL] h+ h-
84 for ( int index = 0; index < 3; index++ )
85 if ( ( dau[index] == _K0B ) || ( dau[index] == _KS ) || ( dau[index] == _KL ) )
86 _d1 = index;
87 else if ( ( dau[index] == _PIP ) || ( dau[index] == _KP ) ) _d2 = index;
88 else if ( ( dau[index] == _PIM ) || ( dau[index] == _KM ) ) _d3 = index;
89 else reportInvalidAndExit();
90 else if ( parId == _D0B ) // Look for K0 h+ h-. The order must be K[0SL] h- h+
91 for ( int index = 0; index < 3; index++ )
92 if ( ( dau[index] == _K0 ) || ( dau[index] == _KS ) || ( dau[index] == _KL ) )
93 _d1 = index;
94 else if ( ( dau[index] == _PIM ) || ( dau[index] == _KM ) ) _d2 = index;
95 else if ( ( dau[index] == _PIP ) || ( dau[index] == _KP ) ) _d3 = index;
96 else reportInvalidAndExit();
97 else reportInvalidAndExit();
98
99 // Check if we're dealing with Ks pi pi or with Ks K K.
100 _isKsPiPi = false;
101 if ( dau[_d2] == _PIP || dau[_d2] == _PIM ) _isKsPiPi = true;
102}
103
105 // Same structure for all of these decays.
107 EvtVector4R pA = part->getDaug( _d1 )->getP4();
108 EvtVector4R pB = part->getDaug( _d2 )->getP4();
109 EvtVector4R pC = part->getDaug( _d3 )->getP4();
110
111 // Squared invariant masses.
112 double m2AB = ( pA + pB ).mass2();
113 double m2AC = ( pA + pC ).mass2();
114 double m2BC = ( pB + pC ).mass2();
115
116 // Dalitz amplitudes of the decay of the particle and that of the antiparticle.
117 EvtComplex ampDalitz;
118 EvtComplex ampAntiDalitz;
119
120 if ( _isKsPiPi )
121 { // For Ks pi pi
122 EvtDalitzPoint point( _mKs, _mPi, _mPi, m2AB, m2BC, m2AC );
123 EvtDalitzPoint antiPoint( _mKs, _mPi, _mPi, m2AC, m2BC, m2AB );
124
125 ampDalitz = dalitzKsPiPi( point );
126 ampAntiDalitz = dalitzKsPiPi( antiPoint );
127 }
128 else
129 { // For Ks K K
130 EvtDalitzPoint point( _mKs, _mK, _mK, m2AB, m2BC, m2AC );
131 EvtDalitzPoint antiPoint( _mKs, _mK, _mK, m2AC, m2BC, m2AB );
132
133 ampDalitz = dalitzKsKK( point );
134 ampAntiDalitz = dalitzKsKK( antiPoint );
135 }
136
137 //_i1 += ampDalitz * conj( ampDalitz ) / 1.e8;
138 //_iChi += ampAntiDalitz * conj( ampDalitz ) / 1.e8;
139 //_iChi2 += ampAntiDalitz * conj( ampAntiDalitz ) / 1.e8;
140
141 // std::cout << "INTEGRALS: " << _i1 << " " << _iChi << " " << _iChi2 << " " << _iChi / _i1
142 // << " " << _iChi2 / _i1 << std::endl;
143
144 // Assume there's no direct CP violation.
145 EvtComplex barAOverA = ampAntiDalitz / ampDalitz;
146
147 // CP violation in the interference. _qp implements CP violation in the mixing.
148 EvtComplex chi = _qp * barAOverA;
149
150 // Generate a negative exponential life time. p( gt ) = ( 1 - y ) * e^{ - ( 1 - y ) gt }
151 double gt = -log( EvtRandom::Flat() ) / ( 1. - _y );
152 part->setLifetime( gt / _gamma );
153
154 // Compute time dependent amplitude.
155 EvtComplex amp = .5 * ampDalitz * exp( -_y * gt / 2. ) *
156 ( ( 1. + chi ) * h1( gt ) + ( 1. - chi ) * h2( gt ) );
157
158 vertex( amp );
159
160 return;
161}
162
163void EvtD0mixDalitz::readPDGValues() {
164 // Define the EvtIds.
165 _D0 = EvtPDL::getId( "D0" );
166 _D0B = EvtPDL::getId( "anti-D0" );
167 _KM = EvtPDL::getId( "K-" );
168 _KP = EvtPDL::getId( "K+" );
169 _K0 = EvtPDL::getId( "K0" );
170 _K0B = EvtPDL::getId( "anti-K0" );
171 _KL = EvtPDL::getId( "K_L0" );
172 _KS = EvtPDL::getId( "K_S0" );
173 _PIM = EvtPDL::getId( "pi-" );
174 _PIP = EvtPDL::getId( "pi+" );
175
176 // Read the relevant masses.
177 _mD0 = EvtPDL::getMass( _D0 );
178 _mKs = EvtPDL::getMass( _KS );
179 _mPi = EvtPDL::getMass( _PIP );
180 _mK = EvtPDL::getMass( _KP );
181
182 // Compute the decay rate from the parameter in the evt.pdl file.
183 _ctau = EvtPDL::getctau( EvtPDL::getId( "D0" ) );
184
185 //_iChi = _qp * EvtComplex( 0.089723 , 0.0004776 ); // All resonances RBW, also Rho0.
186
187 //_iChi = _qp * EvtComplex( 0.0481807, 0.0003043 ); // KStarm only
188 //_iChi = _qp * EvtComplex( 0.0594099, 0.00023803 ); // All resonances RBW but GS Rho
189 //_iChi = _qp * EvtComplex( 0.0780186, 0.000417646 ); // All resonances for KsKK
190 //_iChi2 = _qp * 1.;
191
192 /*
193 // Compute the gamma correction factor avgBeta = Gamma tau.
194 // Compute the norm of the unnormalized p(\beta).
195 double factorY = ( 1. + abs( _iChi2 ) ) / 2. - _y * real( _iChi );
196 double factorX = ( 1. - abs( _iChi2 ) ) / 2. + _x * imag( _iChi );
197 double norm = factorY / ( 1. - pow( _y, 2 ) ) + factorX / ( 1. + pow( _x, 2 ) );
198
199 // Compute the integral of p(\beta) \beta d\beta.
200 double termY = ( 1. + abs( _iChi2 ) ) / 2. - 2. * _y / ( 1. + pow( _y, 2 ) ) * real( _iChi );
201 double termX = ( 1. - abs( _iChi2 ) ) / 2. + 2. * _x / ( 1. - pow( _x, 2 ) ) * imag( _iChi );
202 double quotientY = ( 1. + pow( _y, 2 ) ) / pow( 1. - pow( _y, 2 ), 2 );
203 double quotientX = ( 1. - pow( _x, 2 ) ) / pow( 1. + pow( _x, 2 ), 2 );
204 double normTimesAvg = termY * quotientY + termX * quotientX;
205
206 double avgBeta = normTimesAvg / norm;
207
208 _gamma = avgBeta / _ctau;
209 */
210
211 _gamma = 1. / _ctau; // ALERT: Gamma is not 1 / tau.
212}
213
214EvtComplex EvtD0mixDalitz::dalitzKsPiPi( const EvtDalitzPoint& point ) {
215 static const EvtDalitzPlot plot( _mKs, _mPi, _mPi, _mD0 );
216
217 EvtComplex amp = 0.;
218
219 if ( _isRBWmodel )
220 {
221 // This corresponds to relativistic Breit-Wigner distributions. Not K-matrix.
222 // Defining resonances.
223 static EvtDalitzReso KStarm( plot, _BC, _AC, _VECTOR, 0.893606, 0.0463407, _RBW );
224 static EvtDalitzReso KStarp( plot, _BC, _AB, _VECTOR, 0.893606, 0.0463407, _RBW );
225 static EvtDalitzReso rho0( plot, _AC, _BC, _VECTOR, 0.7758, 0.1464, _GS );
226 static EvtDalitzReso omega( plot, _AC, _BC, _VECTOR, 0.78259, 0.00849, _RBW );
227 static EvtDalitzReso f0_980( plot, _AC, _BC, _SCALAR, 0.975, 0.044, _RBW );
228 static EvtDalitzReso f0_1370( plot, _AC, _BC, _SCALAR, 1.434, 0.173, _RBW );
229 static EvtDalitzReso f2_1270( plot, _AC, _BC, _TENSOR, 1.2754, 0.1851, _RBW );
230 static EvtDalitzReso K0Starm_1430( plot, _BC, _AC, _SCALAR, 1.459, 0.175, _RBW );
231 static EvtDalitzReso K0Starp_1430( plot, _BC, _AB, _SCALAR, 1.459, 0.175, _RBW );
232 static EvtDalitzReso K2Starm_1430( plot, _BC, _AC, _TENSOR, 1.4256, 0.0985, _RBW );
233 static EvtDalitzReso K2Starp_1430( plot, _BC, _AB, _TENSOR, 1.4256, 0.0985, _RBW );
234 static EvtDalitzReso sigma( plot, _AC, _BC, _SCALAR, 0.527699, 0.511861, _RBW );
235 static EvtDalitzReso sigma2( plot, _AC, _BC, _SCALAR, 1.03327, 0.0987890, _RBW );
236 static EvtDalitzReso KStarm_1680( plot, _BC, _AC, _VECTOR, 1.677, 0.205, _RBW );
237
238 // Adding terms to the amplitude with their corresponding amplitude and phase terms.
239 amp += EvtComplex( .848984, .893618 );
240 amp += EvtComplex( -1.16356, 1.19933 ) * KStarm.evaluate( point );
241 amp += EvtComplex( .106051, -.118513 ) * KStarp.evaluate( point );
242 amp += EvtComplex( 1.0, 0.0 ) * rho0.evaluate( point );
243 amp += EvtComplex( -.0249569, .0388072 ) * omega.evaluate( point );
244 amp += EvtComplex( -.423586, -.236099 ) * f0_980.evaluate( point );
245 amp += EvtComplex( -2.16486, 3.62385 ) * f0_1370.evaluate( point );
246 amp += EvtComplex( .217748, -.133327 ) * f2_1270.evaluate( point );
247 amp += EvtComplex( 1.62128, 1.06816 ) * K0Starm_1430.evaluate( point );
248 amp += EvtComplex( .148802, .0897144 ) * K0Starp_1430.evaluate( point );
249 amp += EvtComplex( 1.15489, -.773363 ) * K2Starm_1430.evaluate( point );
250 amp += EvtComplex( .140865, -.165378 ) * K2Starp_1430.evaluate( point );
251 amp += EvtComplex( -1.55556, -.931685 ) * sigma.evaluate( point );
252 amp += EvtComplex( -.273791, -.0535596 ) * sigma2.evaluate( point );
253 amp += EvtComplex( -1.69720, .128038 ) * KStarm_1680.evaluate( point );
254 }
255 else
256 {
257 // This corresponds to the complete model (RBW, GS, LASS and K-matrix).
258 // Defining resonances.
259 static EvtDalitzReso KStarm( plot, _BC, _AC, _VECTOR, 0.893619, 0.0466508, _RBW );
260 static EvtDalitzReso KStarp( plot, _BC, _AB, _VECTOR, 0.893619, 0.0466508, _RBW );
261 static EvtDalitzReso rho0( plot, _AC, _BC, _VECTOR, 0.7758, 0.1464, _GS );
262 static EvtDalitzReso omega( plot, _AC, _BC, _VECTOR, 0.78259, 0.00849, _RBW );
263 static EvtDalitzReso f2_1270( plot, _AC, _BC, _TENSOR, 1.2754, 0.1851, _RBW );
264 static EvtDalitzReso K0Starm_1430( plot, _AC, 1.46312, 0.232393, 1.0746, -1.83214, .803516,
265 2.32788, 1., -5.31306 ); // LASS
266 static EvtDalitzReso K0Starp_1430( plot, _AB, 1.46312, 0.232393, 1.0746, -1.83214, .803516,
267 2.32788, 1., -5.31306 ); // LASS
268 static EvtDalitzReso K2Starm_1430( plot, _BC, _AC, _TENSOR, 1.4256, 0.0985, _RBW );
269 static EvtDalitzReso K2Starp_1430( plot, _BC, _AB, _TENSOR, 1.4256, 0.0985, _RBW );
270 static EvtDalitzReso KStarm_1680( plot, _BC, _AC, _VECTOR, 1.677, 0.205, _RBW );
271
272 // Defining K-matrix.
273 static EvtComplex fr12( 1.87981, -.628378 );
274 static EvtComplex fr13( 4.3242, 2.75019 );
275 static EvtComplex fr14( 3.22336, .271048 );
276 static EvtComplex fr15( .0, .0 );
277 static EvtDalitzReso Pole1( plot, _BC, "Pole1", _KMAT, fr12, fr13, fr14, fr15, -.0694725 );
278 static EvtDalitzReso Pole2( plot, _BC, "Pole2", _KMAT, fr12, fr13, fr14, fr15, -.0694725 );
279 static EvtDalitzReso Pole3( plot, _BC, "Pole3", _KMAT, fr12, fr13, fr14, fr15, -.0694725 );
280 static EvtDalitzReso Pole4( plot, _BC, "Pole4", _KMAT, fr12, fr13, fr14, fr15, -.0694725 );
281 static EvtDalitzReso kmatrix( plot, _BC, "f11prod", _KMAT, fr12, fr13, fr14, fr15,
282 -.0694725 );
283
284 // Adding terms to the amplitude with their corresponding amplitude and phase terms.
285 amp += EvtComplex( -1.31394, 1.14072 ) * KStarm.evaluate( point );
286 amp += EvtComplex( .116239, -.107287 ) * KStarp.evaluate( point );
287 amp += EvtComplex( 1.0, 0.0 ) * rho0.evaluate( point );
288 amp += EvtComplex( -.0313343, .0424013 ) * omega.evaluate( point );
289 amp += EvtComplex( .559412, -.232336 ) * f2_1270.evaluate( point );
290 amp += EvtComplex( 7.35400, -3.67637 ) * K0Starm_1430.evaluate( point );
291 amp += EvtComplex( .255913, -.190459 ) * K0Starp_1430.evaluate( point );
292 amp += EvtComplex( 1.05397, -.936297 ) * K2Starm_1430.evaluate( point );
293 amp += EvtComplex( -.00760136, -.0908624 ) * K2Starp_1430.evaluate( point );
294 amp += EvtComplex( -1.45336, -.164494 ) * KStarm_1680.evaluate( point );
295 amp += EvtComplex( -1.81830, 9.10680 ) * Pole1.evaluate( point );
296 amp += EvtComplex( 10.1751, 3.87961 ) * Pole2.evaluate( point );
297 amp += EvtComplex( 23.6569, -4.94551 ) * Pole3.evaluate( point );
298 amp += EvtComplex( .0725431, -9.16264 ) * Pole4.evaluate( point );
299 amp += EvtComplex( -2.19449, -7.62666 ) * kmatrix.evaluate( point );
300
301 amp *= .97; // Multiply by a constant in order to use the same maximum as RBW model.
302 }
303
304 return amp;
305}
306
307EvtComplex EvtD0mixDalitz::dalitzKsKK( const EvtDalitzPoint& point ) {
308 static const EvtDalitzPlot plot( _mKs, _mK, _mK, _mD0 );
309
310 // Defining resonances.
311 static EvtDalitzReso a00_980( plot, _AC, _BC, _SCALAR, 0.999, _RBW, .550173, .324, _EtaPic );
312 static EvtDalitzReso phi( plot, _AC, _BC, _VECTOR, 1.01943, .00459319, _RBW );
313 static EvtDalitzReso a0p_980( plot, _AC, _AB, _SCALAR, 0.999, _RBW, .550173, .324, _EtaPic );
314 static EvtDalitzReso f0_1370( plot, _AC, _BC, _SCALAR, 1.350, .265, _RBW );
315 static EvtDalitzReso a0m_980( plot, _AB, _AC, _SCALAR, 0.999, _RBW, .550173, .324, _EtaPic );
316 static EvtDalitzReso f0_980( plot, _AC, _BC, _SCALAR, 0.965, _RBW, .695, .165, _PicPicKK );
317 static EvtDalitzReso f2_1270( plot, _AC, _BC, _TENSOR, 1.2754, .1851, _RBW );
318 static EvtDalitzReso a00_1450( plot, _AC, _BC, _SCALAR, 1.474, .265, _RBW );
319 static EvtDalitzReso a0p_1450( plot, _AC, _AB, _SCALAR, 1.474, .265, _RBW );
320 static EvtDalitzReso a0m_1450( plot, _AB, _AC, _SCALAR, 1.474, .265, _RBW );
321
322 // Adding terms to the amplitude with their corresponding amplitude and phase terms.
323 EvtComplex amp( 0., 0. ); // Phase space amplitude.
324 amp += EvtComplex( 1.0, 0.0 ) * a00_980.evaluate( point );
325 amp += EvtComplex( -.126314, .188701 ) * phi.evaluate( point );
326 amp += EvtComplex( -.561428, .0135338 ) * a0p_980.evaluate( point );
327 amp += EvtComplex( .035, -.00110488 ) * f0_1370.evaluate( point );
328 amp += EvtComplex( -.0872735, .0791190 ) * a0m_980.evaluate( point );
329 amp += EvtComplex( 0., 0. ) * f0_980.evaluate( point );
330 amp += EvtComplex( .257341, -.0408343 ) * f2_1270.evaluate( point );
331 amp += EvtComplex( -.0614342, -.649930 ) * a00_1450.evaluate( point );
332 amp += EvtComplex( -.104629, .830120 ) * a0p_1450.evaluate( point );
333 amp += EvtComplex( 0., 0. ) * a0m_1450.evaluate( point );
334
335 return 2.8 * amp; // Multiply by 2.8 in order to reuse the same probmax as Ks pi pi.
336}
337
338// < f | H | D^0 (t) > = 1/2 * [ ( 1 + \chi_f ) * A_f * e_1(gt) + ( 1 - \chi_f ) * A_f *
339// e_2(gt) ] < f | H | D^0 (t) > = 1/2 * exp( -gamma t / 2 ) * [ ( 1 + \chi_f ) * A_f * h_1(t)
340// + ( 1 - \chi_f ) * A_f * h_2(t) ] e{1,2}( gt ) = exp( -gt / 2 ) * h{1,2}( gt ).
341EvtComplex EvtD0mixDalitz::h1( const double& gt ) const {
342 return exp( -EvtComplex( _y, _x ) * gt / 2. );
343}
344
345EvtComplex EvtD0mixDalitz::h2( const double& gt ) const {
346 return exp( EvtComplex( _y, _x ) * gt / 2. );
347}
EvtComplex exp(const EvtComplex &c)
ostream & report(Severity severity, const char *facility)
Definition EvtReport.cc:34
@ ERROR
Definition EvtReport.hh:49
void decay(EvtParticle *p)
void vertex(const EvtComplex &amp)
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
void checkSpinParent(EvtSpinType::spintype sp)
double getArg(int j)
EvtId getParentId()
void checkNDaug(int d1, int d2=-1)
EvtId * getDaugs()
EvtId getDaug(int i)
Definition EvtId.hh:27
static double getMass(EvtId i)
Definition EvtPDL.hh:44
static EvtId getId(const std::string &name)
Definition EvtPDL.cc:272
static double getctau(EvtId i)
Definition EvtPDL.hh:60
const EvtVector4R & getP4() const
void setLifetime(double tau)
EvtParticle * getDaug(int i)
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
static double Flat(double min, double max)
Definition EvtRandom.cc:55