[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtD0gammaDalitz.cpp

//--------------------------------------------------------------------------
//
// Module: EvtD0gammaDalitz.cc
//
// Modification history:
//
//    JGT     February 13, 2012         Module created
//
//------------------------------------------------------------------------

#include <cstdlib>
#include <string>
#include <cmath>

#include "EvtGenBase/EvtPatches.hh"

#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtResonance.hh"
#include "EvtGenBase/EvtResonance2.hh"
#include "EvtGenModels/EvtD0gammaDalitz.hh"
#include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtFlatte.hh"
#include "EvtGenBase/EvtDecayTable.hh"


// Initialize the static variables.
const EvtSpinType::spintype& EvtD0gammaDalitz::_SCALAR = EvtSpinType::SCALAR;
const EvtSpinType::spintype& EvtD0gammaDalitz::_VECTOR = EvtSpinType::VECTOR;
const EvtSpinType::spintype& EvtD0gammaDalitz::_TENSOR = EvtSpinType::TENSOR;

const EvtDalitzReso::CouplingType& EvtD0gammaDalitz::_EtaPic   = EvtDalitzReso::EtaPic;
const EvtDalitzReso::CouplingType& EvtD0gammaDalitz::_PicPicKK = EvtDalitzReso::PicPicKK;

const EvtDalitzReso::NumType& EvtD0gammaDalitz::_RBW   = EvtDalitzReso::RBW_CLEO_ZEMACH;
const EvtDalitzReso::NumType& EvtD0gammaDalitz::_GS    = EvtDalitzReso::GS_CLEO_ZEMACH;
const EvtDalitzReso::NumType& EvtD0gammaDalitz::_KMAT  = EvtDalitzReso::K_MATRIX;

const EvtCyclic3::Pair& EvtD0gammaDalitz::_AB = EvtCyclic3::AB;
const EvtCyclic3::Pair& EvtD0gammaDalitz::_AC = EvtCyclic3::AC;
const EvtCyclic3::Pair& EvtD0gammaDalitz::_BC = EvtCyclic3::BC;


EvtD0gammaDalitz::EvtD0gammaDalitz()
{
  /* Empty constructor. */
}


EvtD0gammaDalitz::~EvtD0gammaDalitz()
{
  /* Empty destructor. */
}


std::string EvtD0gammaDalitz::getName()
{
  return "D0GAMMADALITZ";
}


EvtDecayBase* EvtD0gammaDalitz::clone()
{
  return new EvtD0gammaDalitz;
}


void EvtD0gammaDalitz::init()
{
  // check that there are 0 arguments
  checkNArg(0);

  // Check that this model is valid for the specified decay.
  checkNDaug( 3 );
  checkSpinParent  (    _SCALAR );
  checkSpinDaughter( 0, _SCALAR );
  checkSpinDaughter( 1, _SCALAR );
  checkSpinDaughter( 2, _SCALAR );

  // Get the values of the EvtId objects from the data files.
  readPDGValues();

  // Get the EvtId of the D0 and its 3 daughters.
  getParentId();

  EvtId dau[ 3 ];
  for ( int index = 0; index < 3; index++ )
  {
    dau[ index ] = getDaug( index );
  }

  // Look for K0bar h+ h-. The order will be K[0SL] h+ h-
  for ( int index = 0; index < 3; index++ )
  {
    if      ( ( dau[ index ] == _K0B ) || ( dau[ index ] == _KS ) || ( dau[ index ] == _KL ) )
    {
      _d1 = index;
    }
    else if ( ( dau[ index ] == _PIP ) || ( dau[ index ] == _KP ) )
    {
      _d2 = index;
    }
    else if ( ( dau[ index ] == _PIM ) || ( dau[ index ] == _KM ) )
    {
      _d3 = index;
    }
    else
    {
      reportInvalidAndExit();
    }
  }

  // Check if we're dealing with Ks pi pi or with Ks K K.
  _isKsPiPi = false;
  if ( dau[ _d2 ] == _PIP || dau[ _d2 ] == _PIM )
  {
    _isKsPiPi = true;
  }

}


void EvtD0gammaDalitz::initProbMax()
{
  setProbMax( 5200. );
}


void EvtD0gammaDalitz::decay( EvtParticle* part )
{
  // Check if the D is from a B+- -> D0 K+- decay with the appropriate model.
  EvtParticle* parent = part->getParent(); // If there are no mistakes, should be B+ or B-.
  if (parent != 0 && EvtDecayTable::getInstance()->getDecayFunc( parent )->getName() == "BTODDALITZCPK" )
  {
    EvtId parId = parent->getId();
    if ( ( parId == _BP ) || ( parId == _BM ) ||
         ( parId == _B0 ) || ( parId == _B0B) )
    {
      _bFlavor = parId;
    }
    else
    {
      reportInvalidAndExit();
    }
  }
  else
  {
    reportInvalidAndExit();
  }

  // Read the D decay parameters from the B decay model.
  // Gamma angle in rad.
  double gamma = EvtDecayTable::getInstance()->getDecayFunc( parent )->getArg( 0 );
  // Strong phase in rad.
  double delta = EvtDecayTable::getInstance()->getDecayFunc( parent )->getArg( 1 );
  // Ratio between B->D0K and B->D0barK
  double rB    = EvtDecayTable::getInstance()->getDecayFunc( parent )->getArg( 2 );

  // Same structure for all of these decays.
  part->initializePhaseSpace( getNDaug(), getDaugs() );
  EvtVector4R pA = part->getDaug( _d1 )->getP4();
  EvtVector4R pB = part->getDaug( _d2 )->getP4();
  EvtVector4R pC = part->getDaug( _d3 )->getP4();

  // Squared invariant masses.
  double mSqAB = ( pA + pB ).mass2();
  double mSqAC = ( pA + pC ).mass2();
  double mSqBC = ( pB + pC ).mass2();

  EvtComplex amp( 1.0, 0.0 );

  // Direct and conjugated amplitudes.
  EvtComplex ampDir;
  EvtComplex ampCnj;

  if ( _isKsPiPi )
  {
    // Direct and conjugated Dalitz points.
    EvtDalitzPoint pointDir( _mKs, _mPi, _mPi, mSqAB, mSqBC, mSqAC );
    EvtDalitzPoint pointCnj( _mKs, _mPi, _mPi, mSqAC, mSqBC, mSqAB );

    // Direct and conjugated amplitudes.
    ampDir = dalitzKsPiPi( pointDir );
    ampCnj = dalitzKsPiPi( pointCnj );
  }
  else
  {
    // Direct and conjugated Dalitz points.
    EvtDalitzPoint pointDir( _mKs, _mK, _mK, mSqAB, mSqBC, mSqAC );
    EvtDalitzPoint pointCnj( _mKs, _mK, _mK, mSqAC, mSqBC, mSqAB );

    // Direct and conjugated amplitudes.
    ampDir = dalitzKsKK( pointDir );
    ampCnj = dalitzKsKK( pointCnj );
  }

  if ( _bFlavor == _BP || _bFlavor == _B0 )
  {
    amp = ampCnj + rB * exp( EvtComplex( 0., delta + gamma ) ) * ampDir;
  }
  else
  {
    amp = ampDir + rB * exp( EvtComplex( 0., delta - gamma ) ) * ampCnj;
  }

  vertex( amp );

  return;

}


EvtComplex EvtD0gammaDalitz::dalitzKsPiPi( const EvtDalitzPoint& point ) const
{
  static const EvtDalitzPlot plot( _mKs, _mPi, _mPi, _mD0 );

  EvtComplex amp = 0.;

  // This corresponds to relativistic Breit-Wigner distributions. Not K-matrix.
  // Defining resonances.
  static EvtDalitzReso KStarm      ( plot, _BC, _AC, _VECTOR, 0.893606, 0.0463407, _RBW );
  static EvtDalitzReso KStarp      ( plot, _BC, _AB, _VECTOR, 0.893606, 0.0463407, _RBW );
  static EvtDalitzReso rho0        ( plot, _AC, _BC, _VECTOR, 0.7758  , 0.1464   , _GS  );
  static EvtDalitzReso omega       ( plot, _AC, _BC, _VECTOR, 0.78259 , 0.00849  , _RBW );
  static EvtDalitzReso f0_980      ( plot, _AC, _BC, _SCALAR, 0.975   , 0.044    , _RBW );
  static EvtDalitzReso f0_1370     ( plot, _AC, _BC, _SCALAR, 1.434   , 0.173    , _RBW );
  static EvtDalitzReso f2_1270     ( plot, _AC, _BC, _TENSOR, 1.2754  , 0.1851   , _RBW );
  static EvtDalitzReso K0Starm_1430( plot, _BC, _AC, _SCALAR, 1.459   , 0.175    , _RBW );
  static EvtDalitzReso K0Starp_1430( plot, _BC, _AB, _SCALAR, 1.459   , 0.175    , _RBW );
  static EvtDalitzReso K2Starm_1430( plot, _BC, _AC, _TENSOR, 1.4256  , 0.0985   , _RBW );
  static EvtDalitzReso K2Starp_1430( plot, _BC, _AB, _TENSOR, 1.4256  , 0.0985   , _RBW );
  static EvtDalitzReso sigma       ( plot, _AC, _BC, _SCALAR, 0.527699, 0.511861 , _RBW );
  static EvtDalitzReso sigma2      ( plot, _AC, _BC, _SCALAR, 1.03327 , 0.0987890, _RBW );
  static EvtDalitzReso KStarm_1680 ( plot, _BC, _AC, _VECTOR, 1.677   , 0.205    , _RBW );

  // Adding terms to the amplitude with their corresponding amplitude and phase terms.
  amp += EvtComplex(   .848984 ,   .893618  );
  amp += EvtComplex( -1.16356  ,  1.19933   ) * KStarm      .evaluate( point );
  amp += EvtComplex(   .106051 , - .118513  ) * KStarp      .evaluate( point );
  amp += EvtComplex(  1.0      ,  0.0       ) * rho0        .evaluate( point );
  amp += EvtComplex( - .0249569,   .0388072 ) * omega       .evaluate( point );
  amp += EvtComplex( - .423586 , - .236099  ) * f0_980      .evaluate( point );
  amp += EvtComplex( -2.16486  ,  3.62385   ) * f0_1370     .evaluate( point );
  amp += EvtComplex(   .217748 , - .133327  ) * f2_1270     .evaluate( point );
  amp += EvtComplex(  1.62128  ,  1.06816   ) * K0Starm_1430.evaluate( point );
  amp += EvtComplex(   .148802 ,   .0897144 ) * K0Starp_1430.evaluate( point );
  amp += EvtComplex(  1.15489  , - .773363  ) * K2Starm_1430.evaluate( point );
  amp += EvtComplex(   .140865 , - .165378  ) * K2Starp_1430.evaluate( point );
  amp += EvtComplex( -1.55556  , - .931685  ) * sigma       .evaluate( point );
  amp += EvtComplex( - .273791 , - .0535596 ) * sigma2      .evaluate( point );
  amp += EvtComplex( -1.69720  ,   .128038  ) * KStarm_1680 .evaluate( point );

  return amp;
}


EvtComplex EvtD0gammaDalitz::dalitzKsKK( const EvtDalitzPoint& point ) const
{
  static const EvtDalitzPlot plot( _mKs, _mK, _mK, _mD0 );

  // Defining resonances.
  static EvtDalitzReso a00_980 ( plot, _AC, _BC, _SCALAR, 0.999  , _RBW, .550173, .324, _EtaPic   );
  static EvtDalitzReso phi     ( plot, _AC, _BC, _VECTOR, 1.01943,       .00459319    , _RBW      );
  static EvtDalitzReso a0p_980 ( plot, _AC, _AB, _SCALAR, 0.999  , _RBW, .550173, .324, _EtaPic   );
  static EvtDalitzReso f0_1370 ( plot, _AC, _BC, _SCALAR, 1.350  ,       .265         , _RBW      );
  static EvtDalitzReso a0m_980 ( plot, _AB, _AC, _SCALAR, 0.999  , _RBW, .550173, .324, _EtaPic   );
  static EvtDalitzReso f0_980  ( plot, _AC, _BC, _SCALAR, 0.965  , _RBW, .695   , .165, _PicPicKK );
  static EvtDalitzReso f2_1270 ( plot, _AC, _BC, _TENSOR, 1.2754 ,       .1851        , _RBW      );
  static EvtDalitzReso a00_1450( plot, _AC, _BC, _SCALAR, 1.474  ,       .265         , _RBW      );
  static EvtDalitzReso a0p_1450( plot, _AC, _AB, _SCALAR, 1.474  ,       .265         , _RBW      );
  static EvtDalitzReso a0m_1450( plot, _AB, _AC, _SCALAR, 1.474  ,       .265         , _RBW      );

  // Adding terms to the amplitude with their corresponding amplitude and phase terms.
  EvtComplex amp( 0., 0. ); // Phase space amplitude.
  amp += EvtComplex( 1.0          , 0.0        ) * a00_980 .evaluate( point );
  amp += EvtComplex( -.126314     ,  .188701   ) * phi     .evaluate( point );
  amp += EvtComplex( -.561428     ,  .0135338  ) * a0p_980 .evaluate( point );
  amp += EvtComplex(  .035        , -.00110488 ) * f0_1370 .evaluate( point );
  amp += EvtComplex( -.0872735    ,  .0791190  ) * a0m_980 .evaluate( point );
  amp += EvtComplex( 0.           , 0.         ) * f0_980  .evaluate( point );
  amp += EvtComplex(  .257341     , -.0408343  ) * f2_1270 .evaluate( point );
  amp += EvtComplex( -.0614342    , -.649930   ) * a00_1450.evaluate( point );
  amp += EvtComplex( -.104629     ,  .830120   ) * a0p_1450.evaluate( point );
  amp += EvtComplex( 0.           , 0.         ) * a0m_1450.evaluate( point );

  return 2.8 * amp; // Multiply by 2.8 in order to reuse the same probmax as Ks pi pi.
}


void EvtD0gammaDalitz::readPDGValues()
{
  // Define the EvtIds.
  _BP  = EvtPDL::getId( "B+"      );
  _BM  = EvtPDL::getId( "B-"      );
  _B0  = EvtPDL::getId( "B0"      );
  _B0B = EvtPDL::getId( "anti-B0" );
  _D0  = EvtPDL::getId( "D0"      );
  _D0B = EvtPDL::getId( "anti-D0" );
  _KM  = EvtPDL::getId( "K-"      );
  _KP  = EvtPDL::getId( "K+"      );
  _K0  = EvtPDL::getId( "K0"      );
  _K0B = EvtPDL::getId( "anti-K0" );
  _KL  = EvtPDL::getId( "K_L0"    );
  _KS  = EvtPDL::getId( "K_S0"    );
  _PIM = EvtPDL::getId( "pi-"     );
  _PIP = EvtPDL::getId( "pi+"     );

  // Read the relevant masses.
  _mD0 = EvtPDL::getMass( _D0  );
  _mKs = EvtPDL::getMass( _KS  );
  _mPi = EvtPDL::getMass( _PIP );
  _mK  = EvtPDL::getMass( _KP  );
}


void EvtD0gammaDalitz::reportInvalidAndExit() const
{
  report( ERROR, "EvtD0gammaDalitz" ) << "EvtD0gammaDalitz: Invalid mode." << std::endl;
  exit( 1 );
}
Commit	Line	Data
0ca57c2f	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	}