#include "EvtGenBase/EvtPatches.hh"
/*****************************************************************************
 * Project: BaBar detector at the SLAC PEP-II B-factory
 * Package: EvtGenBase
 *    File: $Id: EvtDalitzReso.cc,v 1.4 2008/11/27 16:29:30 jordix Exp $
 *
 * Description:
 *   Class to compute Dalitz amplitudes based on many models that cannot be
 *     handled with EvtResonance.
 *
 * Modification history:
 *   Jordi Garra Ticó     2008/07/03         File created
 *****************************************************************************/


#include <assert.h>
#include <cmath>
#include <iostream>

#include <stdlib.h>
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtMatrix.hh"
#include "EvtGenBase/EvtDalitzReso.hh"

#include "EvtGenBase/EvtdFunction.hh"
#include "EvtGenBase/EvtCyclic3.hh"

#define PRECISION ( 1.e-3 )

using EvtCyclic3::Index;
using EvtCyclic3::Pair;


// single Breit-Wigner
EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, Pair pairAng, Pair pairRes, 
			     EvtSpinType::spintype spin, double m0, double g0, NumType typeN) 
  : _dp(dp),
    _pairAng(pairAng),
    _pairRes(pairRes),
    _spin(spin),
    _typeN(typeN),
    _m0(m0),_g0(g0),
    _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
    _m0_mix(-1.),_g0_mix(0.),_delta_mix(0.),_amp_mix(0.,0.),
    _g1(-1.),_g2(-1.),_coupling2(Undefined),
    _kmatrix_index(-1),_fr12prod(0.,0.),_fr13prod(0.,0.),_fr14prod(0.,0.),_fr15prod(0.,0.),_s0prod(0.),
    _a(0.),_r(0.),_Blass(0.),_phiB(0.),_R(0.),_phiR(0.)
{
  _vb = EvtTwoBodyVertex(_m0,_dp.m(EvtCyclic3::other(_pairRes)),_dp.bigM(),_spin); 
  _vd = EvtTwoBodyVertex(_massFirst,_massSecond,_m0,_spin);
  _vb.set_f( 0.0 ); // Default values for Blatt-Weisskopf factors.
  _vd.set_f( 1.5 );
  assert(_typeN != K_MATRIX && _typeN != K_MATRIX_I && _typeN != K_MATRIX_II);  // single BW cannot be K-matrix
}


// Breit-Wigner with electromagnetic mass mixing
EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, Pair pairAng, Pair pairRes, 
			     EvtSpinType::spintype spin, double m0, double g0, NumType typeN,
			     double m0_mix, double g0_mix, double delta_mix, EvtComplex amp_mix) 
  : _dp(dp),
    _pairAng(pairAng),
    _pairRes(pairRes),
    _spin(spin),
    _typeN(typeN),
    _m0(m0),_g0(g0),
    _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
    _m0_mix(m0_mix),_g0_mix(g0_mix),_delta_mix(delta_mix),_amp_mix(amp_mix),
    _g1(-1.),_g2(-1.),_coupling2(Undefined),
    _kmatrix_index(-1),_fr12prod(0.,0.),_fr13prod(0.,0.),_fr14prod(0.,0.),_fr15prod(0.,0.),_s0prod(0.),
    _a(0.),_r(0.),_Blass(0.),_phiB(0.),_R(0.),_phiR(0.)
{
  _vb = EvtTwoBodyVertex(_m0,_dp.m(EvtCyclic3::other(_pairRes)),_dp.bigM(),_spin); 
  _vd = EvtTwoBodyVertex(_massFirst,_massSecond,_m0,_spin);
  _vb.set_f( 0.0 ); // Default values for Blatt-Weisskopf factors.
  _vd.set_f( 1.5 );
  // single BW (with electromagnetic mixing) cannot be K-matrix
  assert(_typeN != K_MATRIX && _typeN != K_MATRIX_I && _typeN != K_MATRIX_II);
}

// coupled Breit-Wigner
EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, Pair pairAng, Pair pairRes, 
			     EvtSpinType::spintype spin, double m0, NumType typeN, double g1, double g2, CouplingType coupling2)
  : _dp(dp),
    _pairAng(pairAng),
    _pairRes(pairRes),
    _spin(spin),
    _typeN(typeN),
    _m0(m0),_g0(-1.),
    _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
    _m0_mix(-1.),_g0_mix(0.),_delta_mix(0.),_amp_mix(0.,0.),
    _g1(g1),_g2(g2),_coupling2(coupling2),
    _kmatrix_index(-1),_fr12prod(0.,0.),_fr13prod(0.,0.),_fr14prod(0.,0.),_fr15prod(0.,0.),_s0prod(0.),
    _a(0.),_r(0.),_Blass(0.),_phiB(0.),_R(0.),_phiR(0.)
{
  _vb = EvtTwoBodyVertex(_m0,_dp.m(EvtCyclic3::other(_pairRes)),_dp.bigM(),_spin);   
  _vd = EvtTwoBodyVertex(_massFirst,_massSecond,_m0,_spin);
  _vb.set_f( 0.0 ); // Default values for Blatt-Weisskopf factors.
  _vd.set_f( 1.5 );
  assert(_coupling2 != Undefined);
  assert(_typeN != K_MATRIX && _typeN != K_MATRIX_I && _typeN != K_MATRIX_II); // coupled BW cannot be K-matrix
  assert(_typeN != LASS);     // coupled BW cannot be LASS
  assert(_typeN != NBW);      // for coupled BW, only relativistic BW 
}


// K-Matrix (A&S)
EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, Pair pairRes, std::string nameIndex, NumType typeN,
			     EvtComplex fr12prod, EvtComplex fr13prod, EvtComplex fr14prod, EvtComplex fr15prod, double s0prod) 
  : _dp(dp),
    _pairRes(pairRes),
    _typeN(typeN),
    _m0(0.),_g0(0.),
    _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
    _m0_mix(-1.),_g0_mix(0.),_delta_mix(0.),_amp_mix(0.,0.),
    _g1(-1.),_g2(-1.),_coupling2(Undefined),
    _kmatrix_index(-1),_fr12prod(fr12prod),_fr13prod(fr13prod),_fr14prod(fr14prod),_fr15prod(fr15prod),_s0prod(s0prod),
    _a(0.),_r(0.),_Blass(0.),_phiB(0.),_R(0.),_phiR(0.)
{
  assert(_typeN==K_MATRIX || _typeN==K_MATRIX_I || _typeN==K_MATRIX_II);
  _spin=EvtSpinType::SCALAR;
  if (nameIndex=="Pole1") _kmatrix_index=1;
  else if (nameIndex=="Pole2") _kmatrix_index=2;
  else if (nameIndex=="Pole3") _kmatrix_index=3;
  else if (nameIndex=="Pole4") _kmatrix_index=4;
  else if (nameIndex=="Pole5") _kmatrix_index=5;
  else if (nameIndex=="f11prod") _kmatrix_index=6;
  else assert(0);
}


// LASS parameterization
EvtDalitzReso::EvtDalitzReso(const EvtDalitzPlot& dp, Pair pairRes, 
			     double m0, double g0, double a, double r, double B, double phiB, double R, double phiR) 
  : _dp(dp),
    _pairRes(pairRes),
    _typeN(LASS),
    _m0(m0),_g0(g0),
    _massFirst(dp.m(first(pairRes))),_massSecond(dp.m(second(pairRes))),
    _m0_mix(-1.),_g0_mix(0.),_delta_mix(0.),_amp_mix(0.,0.),
    _g1(-1.),_g2(-1.),_coupling2(Undefined),
    _kmatrix_index(-1),_fr12prod(0.,0.),_fr13prod(0.,0.),_fr14prod(0.,0.),_fr15prod(0.,0.),_s0prod(0.),
    _a(a),_r(r),_Blass(B),_phiB(phiB),_R(R),_phiR(phiR)
{
  _spin=EvtSpinType::SCALAR;
  _vd = EvtTwoBodyVertex(_massFirst,_massSecond,_m0,_spin);
  _vd.set_f( 1.5 ); // Default values for Blatt-Weisskopf factors.
}



EvtDalitzReso::EvtDalitzReso(const EvtDalitzReso& other) 
  : _dp(other._dp),
    _pairAng(other._pairAng),
    _pairRes(other._pairRes),
    _spin(other._spin),
    _typeN(other._typeN),
    _m0(other._m0),_g0(other._g0),
    _vb(other._vb),_vd(other._vd),
    _massFirst(other._massFirst),_massSecond(other._massSecond),
    _m0_mix(other._m0_mix),_g0_mix(other._g0_mix),_delta_mix(other._delta_mix),_amp_mix(other._amp_mix),
    _g1(other._g1),_g2(other._g2),_coupling2(other._coupling2),
    _kmatrix_index(other._kmatrix_index),
    _fr12prod(other._fr12prod),_fr13prod(other._fr13prod),_fr14prod(other._fr14prod),_fr15prod(other._fr15prod),
    _s0prod(other._s0prod),
    _a(other._a),_r(other._r),_Blass(other._Blass),_phiB(other._phiB),_R(other._R),_phiR(other._phiR)
{}


EvtDalitzReso::~EvtDalitzReso()
{}


EvtComplex EvtDalitzReso::evaluate(const EvtDalitzPoint& x) 
{
  double m = sqrt(x.q(_pairRes));

  // do use always hash table (speed up fitting)
  if (_typeN==K_MATRIX || _typeN==K_MATRIX_I || _typeN==K_MATRIX_II)
    return Fvector( m*m, _kmatrix_index );

  if (_typeN==LASS)
    return lass(m*m);

  EvtComplex amp(1.0,0.0);

  if (_dp.bigM() != x.bigM()) _vb = EvtTwoBodyVertex(_m0,_dp.m(EvtCyclic3::other(_pairRes)),x.bigM(),_spin); 
  EvtTwoBodyKine vb(m,x.m(EvtCyclic3::other(_pairRes)),x.bigM());
  EvtTwoBodyKine vd(_massFirst,_massSecond,m);   

  EvtComplex prop(0,0);
  if (_typeN==NBW) {
    prop = propBreitWigner(_m0,_g0,m);
  } else if (_typeN==GAUSS_CLEO || _typeN==GAUSS_CLEO_ZEMACH) {
    prop = propGauss(_m0,_g0,m);
  } else {
    if (_coupling2==Undefined) {  
      // single BW
      double g = (_g0<=0. || _vd.pD()<=0.)? -_g0 : _g0*_vd.widthFactor(vd);  // running width
      if (_typeN==GS_CLEO || _typeN==GS_CLEO_ZEMACH) {
	// Gounaris-Sakurai (GS)
	assert(_massFirst==_massSecond);
	prop = propGounarisSakurai(_m0,fabs(_g0),_vd.pD(),m,g,vd.p());
      } else {
	// standard relativistic BW
	prop = propBreitWignerRel(_m0,g,m);
      }
    } else {    
      // coupled width BW
      EvtComplex G1,G2;
      switch (_coupling2) { 
      case PicPic: {
	G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
	static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
	G2 = _g2*_g2*psFactor(mPic,mPic,m);
	break;
      }
      case PizPiz: {
	G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
	static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
	G2 = _g2*_g2*psFactor(mPiz,mPiz,m);
	break;
      }
      case PiPi: {
	G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
	static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
	static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
	G2 = _g2*_g2*psFactor(mPic,mPic,mPiz,mPiz,m);
	break;
      }
      case KcKc: {
	G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
	static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
	G2 = _g2*_g2*psFactor(mKc,mKc,m);
	break;
      }
      case KzKz: {
	G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
	static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
	G2 = _g2*_g2*psFactor(mKz,mKz,m);
	break;
      }
      case KK: {
	G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
	static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
	static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
	G2 = _g2*_g2*psFactor(mKc,mKc,mKz,mKz,m);
	break;
      }
      case EtaPic: {
	G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
	static double mEta = EvtPDL::getMass( EvtPDL::getId( "eta" ) );
	static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
	G2 = _g2*_g2*psFactor(mEta,mPic,m);
	break;
      }
      case EtaPiz: {
	G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
	static double mEta = EvtPDL::getMass( EvtPDL::getId( "eta" ) );
	static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
	G2 = _g2*_g2*psFactor(mEta,mPiz,m);
	break;
      }
      case PicPicKK: {
	static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
	//G1 = _g1*_g1*psFactor(mPic,mPic,m);
	G1 = _g1*psFactor(mPic,mPic,m);
	static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
	static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
	//G2 = _g2*_g2*psFactor(mKc,mKc,mKz,mKz,m);
	G2 = _g2*psFactor(mKc,mKc,mKz,mKz,m);
	break;
      }
      default:
	std::cout << "EvtDalitzReso:evaluate(): PANIC, wrong coupling2 state." << std::endl;
	assert(0);
	break;
      }
      // calculate standard couple BW propagator
      if (_coupling2 != WA76)
	prop = _g1*propBreitWignerRelCoupled(_m0,G1,G2,m);
    } 
  }
  amp *= prop;

  // Compute form-factors (Blatt-Weisskopf penetration factor)
  amp *= _vb.formFactor(vb);  
  amp *= _vd.formFactor(vd);  

  // Compute numerator (angular distribution)
  amp *= numerator(x,vb,vd);  

  // Compute electromagnetic mass mixing factor
  if (_m0_mix>0.) {
    EvtComplex prop_mix;
    if (_typeN==NBW) {
      prop_mix = propBreitWigner(_m0_mix,_g0_mix,m);
    } else {
      assert(_g1<0.); // running width only
      double g_mix = _g0_mix*_vd.widthFactor(vd);
      prop_mix = propBreitWignerRel(_m0_mix,g_mix,m);
    }
    amp *= mixFactor(prop,prop_mix);
  }

  return amp;
}


EvtComplex EvtDalitzReso::psFactor(double & ma, double & mb, double& m)
{
  if (m>(ma+mb)) {
    EvtTwoBodyKine vd(ma,mb,m);
    return EvtComplex(0,2*vd.p()/m);
  } else { 
    // analytical continuation
    double s = m*m;
    double phaseFactor_analyticalCont = -0.5*(sqrt(4*ma*ma/s-1)+sqrt(4*mb*mb/s-1)); 
    return EvtComplex(phaseFactor_analyticalCont,0);
  }
}


EvtComplex EvtDalitzReso::psFactor(double & ma1,double & mb1, double & ma2, double & mb2, double& m)
{
  return 0.5*(psFactor(ma1,mb1,m)+psFactor(ma2,mb2,m));
}


EvtComplex EvtDalitzReso::propGauss(const double& m0, const double& s0, const double& m) 
{
  // Gaussian
  double gauss = 1./sqrt(EvtConst::twoPi)/s0*exp(-(m-m0)*(m-m0)/2./(s0*s0));
  return EvtComplex(gauss,0.);
}


EvtComplex EvtDalitzReso::propBreitWigner(const double& m0, const double& g0, const double& m) 
{
  // non-relativistic BW
  return sqrt(g0/EvtConst::twoPi)/(m-m0-EvtComplex(0.0,g0/2.));
}


EvtComplex EvtDalitzReso::propBreitWignerRel(const double& m0, const double& g0, const double& m) 
{
  // relativistic BW with real width
  return 1./(m0*m0-m*m-EvtComplex(0.,m0*g0));
}



EvtComplex EvtDalitzReso::propBreitWignerRel(const double& m0, const EvtComplex& g0, const double& m) 
{
  // relativistic BW with complex width
  return 1./(m0*m0-m*m-EvtComplex(0.,m0)*g0);
}


EvtComplex EvtDalitzReso::propBreitWignerRelCoupled(const double& m0, const EvtComplex& g1, const EvtComplex& g2, const double& m)
{
  // relativistic coupled BW
  return 1./(m0*m0-m*m-(g1+g2));
}

EvtComplex EvtDalitzReso::propGounarisSakurai(const double& m0, const double& g0, const double& k0,
					    const double& m, const double& g, const double& k) 
{
  // Gounaris-Sakurai parameterization of pi+pi- P wave. PRD, Vol61, 112002. PRL, Vol21, 244.
  // Expressions taken from BAD637v4, after fixing the imaginary part of the BW denominator: i M_R Gamma_R(s) --> i sqrt(s) Gamma_R(s) 
  return (1.+GS_d(m0,k0)*g0/m0)/(m0*m0-m*m-EvtComplex(0.,m*g)+GS_f(m0,g0,k0,m,k));
}


inline double EvtDalitzReso::GS_f(const double& m0, const double& g0, const double& k0, const double& m, const double& k) 
{
  // m: sqrt(s)
  // m0: nominal resonance mass
  // k: momentum of pion in resonance rest frame (at m)
  // k0: momentum of pion in resonance rest frame (at nominal resonance mass)
  return g0*m0*m0/(k0*k0*k0)*( k*k*(GS_h(m,k)-GS_h(m0,k0)) + (m0*m0-m*m)*k0*k0*GS_dhods(m0,k0) );
}

inline double EvtDalitzReso::GS_h(const double& m, const double& k) 
{return 2./EvtConst::pi*k/m*log((m+2.*k)/(2.*_massFirst)) ;}

inline double EvtDalitzReso::GS_dhods(const double& m0, const double& k0)  
{return GS_h(m0,k0)*( 0.125/(k0*k0) - 0.5/(m0*m0) ) + 0.5/(EvtConst::pi*m0*m0) ;}

inline double EvtDalitzReso::GS_d(const double& m0, const double& k0) 
{return 3./EvtConst::pi*_massFirst*_massFirst/(k0*k0)*log((m0+2.*k0)/(2.*_massFirst)) + 
   m0/(2.*EvtConst::pi*k0) - _massFirst*_massFirst*m0/(EvtConst::pi*k0*k0*k0) ;}


EvtComplex EvtDalitzReso::numerator(const EvtDalitzPoint& x, const EvtTwoBodyKine& vb, const EvtTwoBodyKine& vd) 
{
  EvtComplex ret(0.,0.);

  // Non-relativistic Breit-Wigner
  if(NBW == _typeN) {
    ret = angDep(x);
  }

  // Standard relativistic Zemach propagator
  else if(RBW_ZEMACH == _typeN) {
    ret = _vd.phaseSpaceFactor(vd,EvtTwoBodyKine::AB)*angDep(x);
  }

  // Standard relativistic Zemach propagator
  else if(RBW_ZEMACH2 == _typeN) {
    ret = _vd.phaseSpaceFactor(vd,EvtTwoBodyKine::AB)*_vb.phaseSpaceFactor(vb,EvtTwoBodyKine::AB)*angDep(x);
    if(_spin == EvtSpinType::VECTOR) {
      ret *= -4.;
    } else if(_spin == EvtSpinType::TENSOR) {
      ret *= 16./3.;
    } else if(_spin != EvtSpinType::SCALAR)
      assert(0);
  }

  // Kuehn-Santamaria normalization:
  else if(RBW_KUEHN == _typeN) {
    ret = _m0*_m0 * angDep(x);
  }  

  // CLEO amplitude 
  else if( ( RBW_CLEO        == _typeN ) || ( GS_CLEO           == _typeN ) ||
	   ( RBW_CLEO_ZEMACH == _typeN ) || ( GS_CLEO_ZEMACH    == _typeN ) ||
	   ( GAUSS_CLEO      == _typeN ) || ( GAUSS_CLEO_ZEMACH == _typeN)) {

    Index iA = other(_pairAng);           // A = other(BC)
    Index iB = common(_pairRes,_pairAng); // B = common(AB,BC)
    Index iC = other(_pairRes);           // C = other(AB)
    
    double M = x.bigM();
    double mA = x.m(iA);
    double mB = x.m(iB);
    double mC = x.m(iC);
    double qAB = x.q(combine(iA,iB));
    double qBC = x.q(combine(iB,iC));
    double qCA = x.q(combine(iC,iA));

    double M2 = M*M;
    double m02 = ((RBW_CLEO_ZEMACH == _typeN)||(GS_CLEO_ZEMACH == _typeN)||(GAUSS_CLEO_ZEMACH == _typeN))?  qAB : _m0*_m0;
    double mA2 = mA*mA;
    double mB2 = mB*mB;
    double mC2 = mC*mC;
    
    if (_spin == EvtSpinType::SCALAR) ret = EvtComplex(1.,0.);
    else if(_spin == EvtSpinType::VECTOR) {
      ret = qCA - qBC + (M2 - mC2)*(mB2 - mA2)/m02;
    } else if(_spin == EvtSpinType::TENSOR) {
      double x1 = qBC - qCA + (M2 - mC2)*(mA2 - mB2)/m02;       
      double x2 = M2 - mC2;      
      double x3 = qAB - 2*M2 - 2*mC2 + x2*x2/m02;      
      double x4 = mA2 - mB2;
      double x5 = qAB - 2*mB2 - 2*mA2 + x4*x4/m02;
      ret = x1*x1 - x3*x5/3.;
    } else assert(0);
  }
  
  return ret;
}


double EvtDalitzReso::angDep(const EvtDalitzPoint& x)  
{ 
  // Angular dependece for factorizable amplitudes  
  // unphysical cosines indicate we are in big trouble
  double cosTh = x.cosTh(_pairAng,_pairRes);  // angle between common(reso,ang) and other(reso)
  if(fabs(cosTh) > 1.) {
    report(INFO,"EvtGen") << "cosTh " << cosTh << std::endl; 
    assert(0);
  }
  
  // in units of half-spin
  return EvtdFunction::d(EvtSpinType::getSpin2(_spin),2*0,2*0,acos(cosTh));
}


EvtComplex EvtDalitzReso::mixFactor(EvtComplex prop, EvtComplex prop_mix) 
{
  double Delta = _delta_mix*(_m0+_m0_mix);
  return 1/(1-Delta*Delta*prop*prop_mix)*(1+_amp_mix*Delta*prop_mix);
}



EvtComplex EvtDalitzReso::Fvector( double s, int index )
{
  assert(index>=1 && index<=6);

  //Define the complex coupling constant
  //The convection is as follow
  //i=0 --> pi+ pi-
  //i=1 --> KK
  //i=2 --> 4pi
  //i=3 --> eta eta
  //i=4 --> eta eta'
  //The first index is the resonace-pole index
      
  double g[5][5]; // Coupling constants. The first index is the pole index. The second index is the decay channel
  double ma[5];   // Pole masses. The unit is in GeV

  int solution = (_typeN==K_MATRIX)? 3 : (    (_typeN==K_MATRIX_I)? 1 : ( (_typeN==K_MATRIX_II)? 2 : 0 )    ) ;
  if (solution==0) { std::cout << "EvtDalitzReso::Fvector() error. Kmatrix solution incorrectly chosen ! " << std::endl; abort(); } 

  if (solution == 3 ) {

    // coupling constants
    //pi+pi- channel
    g[0][0]=0.22889;
    g[1][0]=0.94128;
    g[2][0]=0.36856;
    g[3][0]=0.33650;
    g[4][0]=0.18171;
    //K+K- channel
    g[0][1]=-0.55377;
    g[1][1]=0.55095;
    g[2][1]=0.23888;
    g[3][1]=0.40907;
    g[4][1]=-0.17558;
    //4pi channel
    g[0][2]=0;
    g[1][2]=0;
    g[2][2]=0.55639;
    g[3][2]=0.85679;
    g[4][2]=-0.79658;
    //eta eta channel
    g[0][3]=-0.39899;
    g[1][3]=0.39065;
    g[2][3]=0.18340;
    g[3][3]=0.19906;
    g[4][3]=-0.00355;
    //eta eta' channel
    g[0][4]=-0.34639;
    g[1][4]=0.31503;
    g[2][4]=0.18681;
    g[3][4]=-0.00984;
    g[4][4]=0.22358;

    // Pole masses
    ma[0]=0.651;      
    ma[1]=1.20360;
    ma[2]=1.55817;
    ma[3]=1.21000;
    ma[4]=1.82206;

  } else if (solution == 1) { // solnI.txt 
    
    // coupling constants
    //pi+pi- channel
    g[0][0]=0.31896;
    g[1][0]=0.85963;
    g[2][0]=0.47993;
    g[3][0]=0.45121;
    g[4][0]=0.39391;
    //K+K- channel
    g[0][1]=-0.49998;
    g[1][1]=0.52402;
    g[2][1]=0.40254;
    g[3][1]=0.42769;
    g[4][1]=-0.30860;
    //4pi channel
    g[0][2]=0;
    g[1][2]=0;
    g[2][2]=1.0;
    g[3][2]=1.15088;
    g[4][2]=0.33999;
    //eta eta channel
    g[0][3]=-0.21554;
    g[1][3]=0.38093;
    g[2][3]=0.21811;
    g[3][3]=0.22925;
    g[4][3]=0.06919;
    //eta eta' channel
    g[0][4]=-0.18294;
    g[1][4]=0.23788;
    g[2][4]=0.05454;
    g[3][4]=0.06444;
    g[4][4]=0.32620;

    // Pole masses
    ma[0]=0.7369;
    ma[1]=1.24347;
    ma[2]=1.62681;
    ma[3]=1.21900;
    ma[4]=1.74932;

  } else if (solution == 2) { // solnIIa.txt 
    
    // coupling constants
    //pi+pi- channel
    g[0][0]=0.26014;
    g[1][0]=0.95289;
    g[2][0]=0.46244;
    g[3][0]=0.41848;
    g[4][0]=0.01804;
    //K+K- channel
    g[0][1]=-0.57849;
    g[1][1]=0.55887;
    g[2][1]=0.31712;
    g[3][1]=0.49910;
    g[4][1]=-0.28430;
    //4pi channel
    g[0][2]=0;
    g[1][2]=0;
    g[2][2]=0.70340;
    g[3][2]=0.96819;
    g[4][2]=-0.90100;
    //eta eta channel
    g[0][3]=-0.32936;
    g[1][3]=0.39910;
    g[2][3]=0.22963;
    g[3][3]=0.24415;
    g[4][3]=-0.07252;
    //eta eta' channel
    g[0][4]=-0.30906;
    g[1][4]=0.31143;
    g[2][4]=0.19802;
    g[3][4]=-0.00522;
    g[4][4]=0.17097;

    // Pole masses
    ma[0]=0.67460;
    ma[1]=1.21094;
    ma[2]=1.57896;
    ma[3]=1.21900;
    ma[4]=1.86602;
  } 

  //Now define the K-matrix pole
  double  rho1sq,rho2sq,rho4sq,rho5sq;
  EvtComplex rho[5];
  double f[5][5];

  //Initalize the mass of the resonance
  double mpi=0.13957;
  double mK=0.493677;     //using charged K value
  double meta=0.54775;    //using PDG value
  double metap=0.95778;   //using PDG value
    
  //Initialize the matrix to value zero
  EvtComplex K[5][5];
  for(int i=0;i<5;i++) { 
    for(int j=0;j<5;j++) {
      K[i][j]=EvtComplex(0,0);
      f[i][j]=0;
    }
  }

  //Input the _f[i][j] scattering data
  double s_scatt=0.0 ; 
  if (solution == 3) 
    s_scatt=-3.92637; 
  else if (solution == 1) 
    s_scatt= -5.0 ;
  else if (solution == 2) 
    s_scatt= -5.0 ; 
  double sa=1.0;
  double sa_0=-0.15;
  if (solution == 3) {
    f[0][0]=0.23399;  // f^scatt
    f[0][1]=0.15044;
    f[0][2]=-0.20545;
    f[0][3]=0.32825;
    f[0][4]=0.35412;
  }else if (solution == 1) {
    f[0][0]=0.04214;  // f^scatt
    f[0][1]=0.19865;
    f[0][2]=-0.63764;
    f[0][3]=0.44063;
    f[0][4]=0.36717;
  }else if (solution == 2) {
    f[0][0]=0.26447;  // f^scatt
    f[0][1]=0.10400;
    f[0][2]=-0.35445;
    f[0][3]=0.31596;
    f[0][4]=0.42483;
  }   
  f[1][0]=f[0][1];
  f[2][0]=f[0][2];
  f[3][0]=f[0][3];
  f[4][0]=f[0][4];

  //Now construct the phase-space factor
  //For eta-eta' there is no difference term
  rho1sq = 1. - pow( mpi + mpi, 2 ) / s;   //pi+ pi- phase factor
  if( rho1sq >= 0 )
    rho[ 0 ] = EvtComplex( sqrt( rho1sq ), 0 );
  else
    rho[ 0 ] = EvtComplex( 0, sqrt( -rho1sq ) );  

  rho2sq = 1. - pow( mK + mK, 2 ) / s;
  if( rho2sq >= 0 )
    rho[ 1 ] = EvtComplex( sqrt( rho2sq ), 0 );
  else
    rho[ 1 ] = EvtComplex( 0, sqrt( -rho2sq ) );

  //using the A&S 4pi phase space Factor:
  //Shit, not continue
  if( s <= 1 )
    {
      double real   = 1.2274 + .00370909 / ( s * s ) - .111203 / s - 6.39017 * s + 16.8358*s*s - 21.8845*s*s*s + 11.3153*s*s*s*s;
      double cont32 = sqrt(1.0-(16.0*mpi*mpi));
      rho[ 2 ] = EvtComplex( cont32 * real, 0 );
    }
  else
    rho[ 2 ] = EvtComplex( sqrt( 1. - 16. * mpi * mpi / s ), 0 );

  rho4sq = 1. - pow( meta + meta, 2 ) / s;
  if( rho4sq >= 0 )
    rho[ 3 ] = EvtComplex( sqrt( rho4sq ), 0 );
  else
    rho[ 3 ] = EvtComplex( 0, sqrt( -rho4sq ) );

  rho5sq = 1. - pow( meta + metap, 2 ) / s;
  if( rho5sq >= 0 )
    rho[ 4 ] = EvtComplex( sqrt( rho5sq ), 0 );
  else
    rho[ 4 ] = EvtComplex( 0, sqrt( -rho5sq ) );

  double smallTerm = 1; // Factor to prevent divergences.

  // Check if some pole may arise problems.
  for ( int pole = 0; pole < 5; pole++ )
    if ( fabs( pow( ma[ pole ], 2 ) - s ) < PRECISION )
      smallTerm = pow( ma[ pole ], 2 ) - s;

  //now sum all the pole
  //equation (3) in the E791 K-matrix paper
  for(int i=0;i<5;i++) { 
    for(int j=0;j<5;j++) {  
      for (int pole_index=0;pole_index<5;pole_index++) {
	double A=g[pole_index][i]*g[pole_index][j];
	double B=ma[pole_index]*ma[pole_index]-s;

	if ( fabs( B ) < PRECISION )
	  K[ i ][ j ] += EvtComplex( A    , 0 );
	else
	  K[ i ][ j ] += EvtComplex( A / B, 0 ) * smallTerm;
      }
    }
  }

  //now add the SVT part
  for(int i=0;i<5;i++) { 
    for(int j=0;j<5;j++) {
      double C=f[i][j]*(1.0-s_scatt);
      double D=(s-s_scatt);
      K[ i ][ j ] += EvtComplex( C / D, 0 ) * smallTerm;
    }
  }

  //Fix the bug in the FOCUS paper
  //Include the Alder zero term:
  for(int i=0;i<5;i++) { 
    for(int j=0;j<5;j++) {
      double E=(s-(sa*mpi*mpi*0.5))*(1.0-sa_0);
      double F=(s-sa_0);    
      K[ i ][ j ] *= EvtComplex(E/F,0);
    }
  }

  //This is not correct!
  //(1-ipK) != (1-iKp)
  static EvtMatrix< EvtComplex > mat;
  mat.setRange( 5 ); // Try to do in only the first time. DEFINE ALLOCATION IN CONSTRUCTOR.

  for ( int row = 0; row < 5; row++ )
    for ( int col = 0; col < 5; col++ )
      mat( row, col ) = ( row == col ) * smallTerm - EvtComplex( 0., 1. ) * K[ row ][ col ] * rho[ col ];


  EvtMatrix< EvtComplex >* matInverse = mat.inverse();  //The 1st row of the inverse matrix. This matrix is {(I-iKp)^-1}_0j
  vector< EvtComplex > U1j;
  for ( int j = 0; j < 5; j++ )
    U1j.push_back( (*matInverse)[ 0 ][ j ] );

  delete matInverse;

  //this calculates final F0 factor
  EvtComplex value( 0, 0 );
  if (index<=5) {
    //this calculates the beta_idx Factors
    for(int j=0;j<5;j++) {        // sum for 5 channel
      EvtComplex top    = U1j[j]*g[index-1][j];
      double     bottom = ma[index-1]*ma[index-1]-s;

      if ( fabs( bottom ) < PRECISION )
	value += top;
      else
	value += top / bottom * smallTerm;
    }
  } else {
    //this calculates fprod Factors
    value += U1j[0];
    value += U1j[1]*_fr12prod;
    value += U1j[2]*_fr13prod;
    value += U1j[3]*_fr14prod;
    value += U1j[4]*_fr15prod;

    value *= (1-_s0prod)/(s-_s0prod) * smallTerm;
  }

  return value;
}


//replace Breit-Wigner with LASS
EvtComplex EvtDalitzReso::lass(double s)
{
  EvtTwoBodyKine vd(_massFirst,_massSecond, sqrt(s));
  double q = vd.p();
  double GammaM = _g0*_vd.widthFactor(vd);  // running width;

  //calculate the background phase motion
  double cot_deltaB = 1.0/(_a*q) + 0.5*_r*q;
  double deltaB = atan( 1.0/cot_deltaB);
  double totalB = deltaB + _phiB ;

  //calculate the resonant phase motion
  double deltaR = atan((_m0*GammaM/(_m0*_m0 - s)));
  double totalR = deltaR + _phiR ;

  //sum them up
  EvtComplex  bkgB,resT;
  bkgB = EvtComplex(_Blass*sin(totalB),0)*EvtComplex(cos(totalB),sin(totalB));
  resT = EvtComplex(_R*sin(deltaR),0)*EvtComplex(cos(totalR),sin(totalR))*EvtComplex(cos(2*totalB),sin(2*totalB));
  EvtComplex T = bkgB + resT;      

  return T;
}