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