TEvtGen/EvtGenModels/EvtD0gammaDalitz.cpp

   1 //--------------------------------------------------------------------------
   2 //
   3 // Module: EvtD0gammaDalitz.cc
   4 //
   5 // Modification history:
   6 //
   7 //    JGT     February 13, 2012         Module created
   8 //
   9 //------------------------------------------------------------------------
  10
  11 #include <cstdlib>
  12 #include <string>
  13 #include <cmath>
  14
  15 #include "EvtGenBase/EvtPatches.hh"
  16
  17 #include "EvtGenBase/EvtParticle.hh"
  18 #include "EvtGenBase/EvtGenKine.hh"
  19 #include "EvtGenBase/EvtPDL.hh"
  20 #include "EvtGenBase/EvtReport.hh"
  21 #include "EvtGenBase/EvtResonance.hh"
  22 #include "EvtGenBase/EvtResonance2.hh"
  23 #include "EvtGenModels/EvtD0gammaDalitz.hh"
  24 #include "EvtGenBase/EvtConst.hh"
  25 #include "EvtGenBase/EvtFlatte.hh"
  26 #include "EvtGenBase/EvtDecayTable.hh"
  27
  28
  29 // Initialize the static variables.
  30 const EvtSpinType::spintype& EvtD0gammaDalitz::_SCALAR = EvtSpinType::SCALAR;
  31 const EvtSpinType::spintype& EvtD0gammaDalitz::_VECTOR = EvtSpinType::VECTOR;
  32 const EvtSpinType::spintype& EvtD0gammaDalitz::_TENSOR = EvtSpinType::TENSOR;
  33
  34 const EvtDalitzReso::CouplingType& EvtD0gammaDalitz::_EtaPic   = EvtDalitzReso::EtaPic;
  35 const EvtDalitzReso::CouplingType& EvtD0gammaDalitz::_PicPicKK = EvtDalitzReso::PicPicKK;
  36
  37 const EvtDalitzReso::NumType& EvtD0gammaDalitz::_RBW   = EvtDalitzReso::RBW_CLEO_ZEMACH;
  38 const EvtDalitzReso::NumType& EvtD0gammaDalitz::_GS    = EvtDalitzReso::GS_CLEO_ZEMACH;
  39 const EvtDalitzReso::NumType& EvtD0gammaDalitz::_KMAT  = EvtDalitzReso::K_MATRIX;
  40
  41 const EvtCyclic3::Pair& EvtD0gammaDalitz::_AB = EvtCyclic3::AB;
  42 const EvtCyclic3::Pair& EvtD0gammaDalitz::_AC = EvtCyclic3::AC;
  43 const EvtCyclic3::Pair& EvtD0gammaDalitz::_BC = EvtCyclic3::BC;
  44
  45
  46 EvtD0gammaDalitz::EvtD0gammaDalitz()
  47 {
  48   /* Empty constructor. */
  49 }
  50
  51
  52 EvtD0gammaDalitz::~EvtD0gammaDalitz()
  53 {
  54   /* Empty destructor. */
  55 }
  56
  57
  58 std::string EvtD0gammaDalitz::getName()
  59 {
  60   return "D0GAMMADALITZ";
  61 }
  62
  63
  64 EvtDecayBase* EvtD0gammaDalitz::clone()
  65 {
  66   return new EvtD0gammaDalitz;
  67 }
  68
  69
  70 void EvtD0gammaDalitz::init()
  71 {
  72   // check that there are 0 arguments
  73   checkNArg(0);
  74
  75   // Check that this model is valid for the specified decay.
  76   checkNDaug( 3 );
  77   checkSpinParent  (    _SCALAR );
  78   checkSpinDaughter( 0, _SCALAR );
  79   checkSpinDaughter( 1, _SCALAR );
  80   checkSpinDaughter( 2, _SCALAR );
  81
  82   // Get the values of the EvtId objects from the data files.
  83   readPDGValues();
  84
  85   // Get the EvtId of the D0 and its 3 daughters.
  86   getParentId();
  87
  88   EvtId dau[ 3 ];
  89   for ( int index = 0; index < 3; index++ )
  90   {
  91     dau[ index ] = getDaug( index );
  92   }
  93
  94   // Look for K0bar h+ h-. The order will be K[0SL] h+ h-
  95   for ( int index = 0; index < 3; index++ )
  96   {
  97     if      ( ( dau[ index ] == _K0B ) || ( dau[ index ] == _KS ) || ( dau[ index ] == _KL ) )
  98     {
  99       _d1 = index;
 100     }
 101     else if ( ( dau[ index ] == _PIP ) || ( dau[ index ] == _KP ) )
 102     {
 103       _d2 = index;
 104     }
 105     else if ( ( dau[ index ] == _PIM ) || ( dau[ index ] == _KM ) )
 106     {
 107       _d3 = index;
 108     }
 109     else
 110     {
 111       reportInvalidAndExit();
 112     }
 113   }
 114
 115   // Check if we're dealing with Ks pi pi or with Ks K K.
 116   _isKsPiPi = false;
 117   if ( dau[ _d2 ] == _PIP || dau[ _d2 ] == _PIM )
 118   {
 119     _isKsPiPi = true;
 120   }
 121
 122 }
 123
 124
 125 void EvtD0gammaDalitz::initProbMax()
 126 {
 127   setProbMax( 5200. );
 128 }
 129
 130
 131 void EvtD0gammaDalitz::decay( EvtParticle* part )
 132 {
 133   // Check if the D is from a B+- -> D0 K+- decay with the appropriate model.
 134   EvtParticle* parent = part->getParent(); // If there are no mistakes, should be B+ or B-.
 135   if (parent != 0 && EvtDecayTable::getInstance()->getDecayFunc( parent )->getName() == "BTODDALITZCPK" )
 136   {
 137     EvtId parId = parent->getId();
 138     if ( ( parId == _BP ) || ( parId == _BM ) ||
 139          ( parId == _B0 ) || ( parId == _B0B) )
 140     {
 141       _bFlavor = parId;
 142     }
 143     else
 144     {
 145       reportInvalidAndExit();
 146     }
 147   }
 148   else
 149   {
 150     reportInvalidAndExit();
 151   }
 152
 153   // Read the D decay parameters from the B decay model.
 154   // Gamma angle in rad.
 155   double gamma = EvtDecayTable::getInstance()->getDecayFunc( parent )->getArg( 0 );
 156   // Strong phase in rad.
 157   double delta = EvtDecayTable::getInstance()->getDecayFunc( parent )->getArg( 1 );
 158   // Ratio between B->D0K and B->D0barK
 159   double rB    = EvtDecayTable::getInstance()->getDecayFunc( parent )->getArg( 2 );
 160
 161   // Same structure for all of these decays.
 162   part->initializePhaseSpace( getNDaug(), getDaugs() );
 163   EvtVector4R pA = part->getDaug( _d1 )->getP4();
 164   EvtVector4R pB = part->getDaug( _d2 )->getP4();
 165   EvtVector4R pC = part->getDaug( _d3 )->getP4();
 166
 167   // Squared invariant masses.
 168   double mSqAB = ( pA + pB ).mass2();
 169   double mSqAC = ( pA + pC ).mass2();
 170   double mSqBC = ( pB + pC ).mass2();
 171
 172   EvtComplex amp( 1.0, 0.0 );
 173
 174   // Direct and conjugated amplitudes.
 175   EvtComplex ampDir;
 176   EvtComplex ampCnj;
 177
 178   if ( _isKsPiPi )
 179   {
 180     // Direct and conjugated Dalitz points.
 181     EvtDalitzPoint pointDir( _mKs, _mPi, _mPi, mSqAB, mSqBC, mSqAC );
 182     EvtDalitzPoint pointCnj( _mKs, _mPi, _mPi, mSqAC, mSqBC, mSqAB );
 183
 184     // Direct and conjugated amplitudes.
 185     ampDir = dalitzKsPiPi( pointDir );
 186     ampCnj = dalitzKsPiPi( pointCnj );
 187   }
 188   else
 189   {
 190     // Direct and conjugated Dalitz points.
 191     EvtDalitzPoint pointDir( _mKs, _mK, _mK, mSqAB, mSqBC, mSqAC );
 192     EvtDalitzPoint pointCnj( _mKs, _mK, _mK, mSqAC, mSqBC, mSqAB );
 193
 194     // Direct and conjugated amplitudes.
 195     ampDir = dalitzKsKK( pointDir );
 196     ampCnj = dalitzKsKK( pointCnj );
 197   }
 198
 199   if ( _bFlavor == _BP || _bFlavor == _B0 )
 200   {
 201     amp = ampCnj + rB * exp( EvtComplex( 0., delta + gamma ) ) * ampDir;
 202   }
 203   else
 204   {
 205     amp = ampDir + rB * exp( EvtComplex( 0., delta - gamma ) ) * ampCnj;
 206   }
 207
 208   vertex( amp );
 209
 210   return;
 211
 212 }
 213
 214
 215 EvtComplex EvtD0gammaDalitz::dalitzKsPiPi( const EvtDalitzPoint& point ) const
 216 {
 217   static const EvtDalitzPlot plot( _mKs, _mPi, _mPi, _mD0 );
 218
 219   EvtComplex amp = 0.;
 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   return amp;
 256 }
 257
 258
 259 EvtComplex EvtD0gammaDalitz::dalitzKsKK( const EvtDalitzPoint& point ) const
 260 {
 261   static const EvtDalitzPlot plot( _mKs, _mK, _mK, _mD0 );
 262
 263   // Defining resonances.
 264   static EvtDalitzReso a00_980 ( plot, _AC, _BC, _SCALAR, 0.999  , _RBW, .550173, .324, _EtaPic   );
 265   static EvtDalitzReso phi     ( plot, _AC, _BC, _VECTOR, 1.01943,       .00459319    , _RBW      );
 266   static EvtDalitzReso a0p_980 ( plot, _AC, _AB, _SCALAR, 0.999  , _RBW, .550173, .324, _EtaPic   );
 267   static EvtDalitzReso f0_1370 ( plot, _AC, _BC, _SCALAR, 1.350  ,       .265         , _RBW      );
 268   static EvtDalitzReso a0m_980 ( plot, _AB, _AC, _SCALAR, 0.999  , _RBW, .550173, .324, _EtaPic   );
 269   static EvtDalitzReso f0_980  ( plot, _AC, _BC, _SCALAR, 0.965  , _RBW, .695   , .165, _PicPicKK );
 270   static EvtDalitzReso f2_1270 ( plot, _AC, _BC, _TENSOR, 1.2754 ,       .1851        , _RBW      );
 271   static EvtDalitzReso a00_1450( plot, _AC, _BC, _SCALAR, 1.474  ,       .265         , _RBW      );
 272   static EvtDalitzReso a0p_1450( plot, _AC, _AB, _SCALAR, 1.474  ,       .265         , _RBW      );
 273   static EvtDalitzReso a0m_1450( plot, _AB, _AC, _SCALAR, 1.474  ,       .265         , _RBW      );
 274
 275   // Adding terms to the amplitude with their corresponding amplitude and phase terms.
 276   EvtComplex amp( 0., 0. ); // Phase space amplitude.
 277   amp += EvtComplex( 1.0          , 0.0        ) * a00_980 .evaluate( point );
 278   amp += EvtComplex( -.126314     ,  .188701   ) * phi     .evaluate( point );
 279   amp += EvtComplex( -.561428     ,  .0135338  ) * a0p_980 .evaluate( point );
 280   amp += EvtComplex(  .035        , -.00110488 ) * f0_1370 .evaluate( point );
 281   amp += EvtComplex( -.0872735    ,  .0791190  ) * a0m_980 .evaluate( point );
 282   amp += EvtComplex( 0.           , 0.         ) * f0_980  .evaluate( point );
 283   amp += EvtComplex(  .257341     , -.0408343  ) * f2_1270 .evaluate( point );
 284   amp += EvtComplex( -.0614342    , -.649930   ) * a00_1450.evaluate( point );
 285   amp += EvtComplex( -.104629     ,  .830120   ) * a0p_1450.evaluate( point );
 286   amp += EvtComplex( 0.           , 0.         ) * a0m_1450.evaluate( point );
 287
 288   return 2.8 * amp; // Multiply by 2.8 in order to reuse the same probmax as Ks pi pi.
 289 }
 290
 291
 292 void EvtD0gammaDalitz::readPDGValues()
 293 {
 294   // Define the EvtIds.
 295   _BP  = EvtPDL::getId( "B+"      );
 296   _BM  = EvtPDL::getId( "B-"      );
 297   _B0  = EvtPDL::getId( "B0"      );
 298   _B0B = EvtPDL::getId( "anti-B0" );
 299   _D0  = EvtPDL::getId( "D0"      );
 300   _D0B = EvtPDL::getId( "anti-D0" );
 301   _KM  = EvtPDL::getId( "K-"      );
 302   _KP  = EvtPDL::getId( "K+"      );
 303   _K0  = EvtPDL::getId( "K0"      );
 304   _K0B = EvtPDL::getId( "anti-K0" );
 305   _KL  = EvtPDL::getId( "K_L0"    );
 306   _KS  = EvtPDL::getId( "K_S0"    );
 307   _PIM = EvtPDL::getId( "pi-"     );
 308   _PIP = EvtPDL::getId( "pi+"     );
 309
 310   // Read the relevant masses.
 311   _mD0 = EvtPDL::getMass( _D0  );
 312   _mKs = EvtPDL::getMass( _KS  );
 313   _mPi = EvtPDL::getMass( _PIP );
 314   _mK  = EvtPDL::getMass( _KP  );
 315 }
 316
 317
 318 void EvtD0gammaDalitz::reportInvalidAndExit() const
 319 {
 320   report( ERROR, "EvtD0gammaDalitz" ) << "EvtD0gammaDalitz: Invalid mode." << std::endl;
 321   exit( 1 );
 322 }