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