TEvtGen/EvtGenModels/EvtD0mixDalitz.cxx

   1 /*****************************************************************************
   2  * Project: BaBar detector at the SLAC PEP-II B-factory
   3  * Package: EvtGenModels
   4  *    File: $Id: EvtD0mixDalitz.cc,v 1.9 2008/11/27 16:37:24 jordix 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       {
  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
 117 void EvtD0mixDalitz::decay( EvtParticle* part )
 118 {
 119   // Same structure for all of these decays.
 120   part->initializePhaseSpace( getNDaug(), getDaugs() );
 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
 176 void 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
 229 EvtComplex 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
 321 EvtComplex 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 ).
 357 EvtComplex EvtD0mixDalitz::h1( const double& gt ) const
 358 {
 359   return exp( - EvtComplex( _y, _x ) * gt / 2. );
 360 }
 361
 362
 363 EvtComplex EvtD0mixDalitz::h2( const double& gt ) const
 364 {
 365   return exp(   EvtComplex( _y, _x ) * gt / 2. );
 366 }
 367