#include "HADRONS++/Current_Library/VA_0_PPP.H"
#include "ATOOLS/Org/Run_Parameter.H"
#include "ATOOLS/Org/MyStrStream.H"
#include "ATOOLS/Org/My_MPI.H"

using namespace HADRONS;
using namespace ATOOLS;
using namespace std;

VA_0_PPP::FF_Base::~FF_Base()
{
}

VA_0_PPP::VA_0_PPP(const ATOOLS::Flavour_Vector& flavs,
                   const std::vector<int>& indices, const std::string& name) :
  Current_Base(flavs, indices, name)
{
  int nPion_0(0), nPion_ch(0), nKaon_0(0), nKaon_ch(0), nKaon_S (0), nKaon_L(0);
  // count number of pions, kaons and calc. mass^2
  for( int i=0; i<3; i++ ) {
    if( m_flavs[p_i[i]].Kfcode() == kf_pi_plus )   nPion_ch++;
    if( m_flavs[p_i[i]].Kfcode() == kf_pi )        nPion_0++;
    if( m_flavs[p_i[i]].Kfcode() == kf_K_plus )    nKaon_ch++;
    if( m_flavs[p_i[i]].Kfcode() == kf_K )         nKaon_0++;
    if( m_flavs[p_i[i]].Kfcode() == kf_K_L )	   nKaon_L++;
    if( m_flavs[p_i[i]].Kfcode() == kf_K_S )       nKaon_S++;
    m_ms[i] = sqr( m_flavs[p_i[i]].HadMass() );
  }
  // sanity check
  if (nPion_ch+nPion_0+nKaon_ch+nKaon_L+nKaon_S+nKaon_0 != 3) {
    msg_Error()<<"ERROR in HADRONS::VA_0_PPP constructor\n"
        <<"     number of three outgoing pseudoscalars != 3 ?!.\n"
        <<"     Don't know, what to do. Will abort."<<endl;
    Abort();           
  }
  // define mode number
  m_mode = nPion_ch*1000 + nPion_0*100 + nKaon_ch*10 + (nKaon_L+nKaon_S+nKaon_0);
  m_kaon_mode = 100*nKaon_L + 10*nKaon_S + nKaon_0;
}

void VA_0_PPP::SetModelParameters( struct GeneralModel _md )
{
  m_Vud      = _md("Vud", Tools::Vud);
  m_Vus      = _md("Vus", Tools::Vus);
  double fpi = _md("fpi", 0.1307 );
  
  m_B123 = 0.;
  double A123;
  switch( m_mode ) {
    case 3000 : /* pi- pi- pi+ mode */
                A123 = m_Vud*SQRT_05;
		m_B123 = 0.;
                break;
    case 1200 : /* pi0 pi0 pi- mode */
                A123 = m_Vud*SQRT_05;
		m_B123 = 0.;
                break;
    case 1020 : /* K- pi- K+ */
                m_B123 = 2.;
                A123 = -0.5*m_Vud;
                break;
    case 1002 : /* K0 pi- K0b */
                m_B123 = -2.;
                A123 = -0.5*m_Vud;
                break;
    case  111 : /* K- pi0 K0 */
                A123 = 1.5*m_Vud*SQRT_05;
		m_B123 = -2./3.;
                break;
    case  210 : /* pi0 pi0 K- mode */
                A123 = m_Vus/4.;
		m_B123 = 1.;
		break;
    case 2010 : /* K- pi- pi+ */
                m_B123 = -1.;
                A123 = -0.5*m_Vus;
                break;
    case 1101 : /* pi- K0b pi0 */
                m_B123 = -2/3.;
                A123 = 1.5*m_Vus*SQRT_05;
                break;
    case   30 : /* K- K- K+ */
                A123 = m_Vus;
                break;
    case   12 : /* K- K0b K0 */
                A123 = -0.5*m_Vus;
                break;
    default   : msg_Error()<<"Warning in HADRONS::Tau_Decay_MEs.C in Tau_Three_Pseudo::GetA123() :\n"
                           <<"     Obviously this three pseudoscalar channel (code "<<m_mode<<")\n"
                           <<"     doesn't have a global A123. Maybe it is not implemented yet.\n"
                           <<"     Take A123=1., will continue and hope for the best."<<endl;
                A123 = 1.;
                break;
  }
  m_global   = 4./(3.*fpi)*A123;
  switch( int(_md("FORM_FACTOR", 1)) ) {
    case 2 : p_ff = new RChT(m_mode,m_kaon_mode,_md,m_ms);
    break;
    case 1 : p_ff = new KS95(m_mode,m_kaon_mode,_md,m_ms);
    break;
    case 3 : p_ff = new KS(m_mode,m_kaon_mode,_md,m_ms);
    break;
  }
}


void VA_0_PPP::Calc(const ATOOLS::Vec4D_Vector& moms, bool m_anti)
{
  Vec4D p1( moms[p_i[0]] ),
        p2( moms[p_i[1]] ),
        p3( moms[p_i[2]] );
  Vec4D Q( p1+p2+p3 );
  double s = (p1+p3).Abs2(),
         t = (p2+p3).Abs2();
  double Q2 = Q.Abs2();
  double dot1 = Q*(p1-p3),
  dot2 = Q*(p2-p3);
  double d1 = dot1/Q2,
  d2 = dot2/Q2;
  Complex F1 = FormFactor( 1, Q2, s, t );
  Complex F2 = FormFactor( 2, Q2, s, t );
  Complex FS = FormFactor( 3, Q2, s, t );
  Complex FV = m_B123*FormFactor( 4, Q2, s, t );
//  F1 = Complex(0.,0.);
//  F2 = Complex(0.,0.);
//  FS = Complex(0.,0.);
  
  Insert( m_global*( FS+F1*(1.-d1)-F2*d2 )*Vec4C(p1)
	+ m_global*( FS-F1*d1+F2*(1.-d2) )*Vec4C(p2)
	+ m_global*( FS-F1*(1.+d1)-F2*(1.+d2) )*Vec4C(p3)
	+ m_global*Complex(0.,1.)*FV*Vec4C(cross(p1,p2,p3)) , 0);
//   for( int h=0; h<4; h++ ) {        // helicity comb. (nutau,tau)
//     _ampls_tensor->push_back( 
//         ( F.X( m_nutau, m_pseudo_1, 0, h, m_cR, m_cL ) * ( FS + F1*(1.-d1) - F2*d2 )
//         + F.X( m_nutau, m_pseudo_2, 0, h, m_cR, m_cL ) * ( FS - F1*d1 + F2*(1.-d2) )
//         + F.X( m_nutau, m_pseudo_3, 0, h, m_cR, m_cL ) * ( FS - F1*(1.+d1) - F2*(1.+d2) ) 
//         + Complex(0.,1.)*F.X( m_nutau, cross(p1,p2,p3), 0, h, m_cR, m_cL ) * FV
//         ) * m_global
//                             );
//   }
}


Complex VA_0_PPP::FormFactor( int j, double Q2, double s, double t )
{
  return( p_ff->FormFactor(j,Q2,s,t) );     
}


VA_0_PPP::FF_Base::FF_Base(int mode, int kaon_mode, GeneralModel _md, double * _ms)
  : m_mode (mode), m_kaon_mode (kaon_mode)
{
  // set resonances and parameters
  m_ms[0] = _ms[0];
  m_ms[1] = _ms[1];
  m_ms[2] = _ms[2];
  m_deltas = 0;
  kf_code resA, resAA, resV[2], resVV[2];
  switch( m_mode ) {
    case 3000 : /* pi- pi- pi+ mode */
      resV[0] = kf_rho_770;
      resV[1] = kf_rho_770;
      break;
    case 1200 : /* pi0 pi0 pi- mode */
      resV[0] = kf_rho_770_plus;
      resV[1] = kf_rho_770_plus;
      break;
    case 1020 : /* K- pi- K+ */
      resV[0] = kf_rho_770;
      resV[1] = kf_K_star_892; 
      break;
    case 1002 : /* K0 pi- K0b */
      resV[0] = kf_rho_770;
      resV[1] = kf_K_star_892_plus; 
      break;
    case  111 : /* K- pi0 K0 */
      resV[0] = kf_rho_770_plus;
      resV[1] = kf_K_star_892;
      break;
    case  210 : /* pi0 pi0 K- mode */
      m_deltas= 1;
      resV[0] = kf_K_star_892_plus;
      resV[1] = kf_K_star_892_plus;
      break;
    case 2010 : /* K- pi- pi+ */
      m_deltas= 1;
      resV[0] = kf_K_star_892; 
      resV[1] = kf_rho_770;
      break;
    case 1101 : /* pi- K0b pi0 */
      m_deltas= 1;
      resV[0] = kf_rho_770_plus;
      resV[1] = kf_K_star_892;
      break;
    case   30 : /* K- K- K+ */
      m_deltas= 1;
      resV[0] = kf_rho_770;
      resV[1] = kf_rho_770;
      break;
    case   12 : /* K- K0b K0 */
      m_deltas= 1;
      resV[0] = kf_rho_770_plus;
      resV[1] = kf_rho_770_plus;
      break;
    default:
      THROW(fatal_error, "Three pseudoscalar mode not recognized: "+ToString(m_mode));
      break;
  }

  // higher resonances
  for( int i=0; i<2; ++i ) {
    if( resV[i] == kf_rho_770 )          resVV[i] = kf_rho_1450;
    if( resV[i] == kf_rho_770_plus )     resVV[i] = kf_rho_1450_plus;
    if( resV[i] == kf_K_star_892 )       resVV[i] = kf_K_star_1410;
    if( resV[i] == kf_K_star_892_plus )  resVV[i] = kf_K_star_1410_plus;
  }

  // axial resonance
  int running    = int( _md("RUNNING_WIDTH", 3 ) );     // running width
  double mA, mAA, wA, wAA;
  if( m_deltas ) {
	running = 0;
	resA  = kf_K_1_1400_plus;
    resAA = kf_K_1_1270_plus;
	mA      = _md("Mass_"+Flavour(resA).IDName(),     1.402 );          // mass of axial resonance
	mAA     = _md("Mass_"+Flavour(resAA).IDName(),    1.270 );         // mass of axial resonance'
	wA      = _md("Width_"+Flavour(resA).IDName(),     0.174 );          // width of axial resonance
	wAA     = _md("Width_"+Flavour(resAA).IDName(),    0.090 );         // width of axial resonance'
	m_alpha        = _md("alpha_"+Flavour(resAA).IDName(),  1.3 );                              // weight factor for A'
  }
  else {
	resA  = kf_a_1_1260_plus;
	resAA = kf_a_1_1260_plus;
	mA      = _md("Mass_"+Flavour(resA).IDName(),     1.251 );          // mass of axial resonance
	mAA     = _md("Mass_"+Flavour(resAA).IDName(),    1.251 );         // mass of axial resonance'
	wA      = _md("Width_"+Flavour(resA).IDName(),     0.599 );          // width of axial resonance
	wAA     = _md("Width_"+Flavour(resAA).IDName(),    0.599 );         // width of axial resonance'
	m_alpha        = _md("alpha_"+Flavour(resAA).IDName(),  0. );                              // weight factor for A'
  }
  m_A = ResonanceFlavour( resA, mA, wA, running&1);
  m_AA = ResonanceFlavour( resAA, mAA, wAA, running&1);
  ResonanceFlavour help_rho = ResonanceFlavour(
      kf_rho_770_plus, 
      _md("Mass_rho(770)+", _md("Mass_rho(770)", 0.773 )),
      _md("Width_rho(770)+", _md("Width_rho(770)", 0.145 )),
      running&2 
      );
  ResonanceFlavour help_rhoP = ResonanceFlavour(
      kf_rho_1450_plus, 
      _md("Mass_rho(1450)+", _md("Mass_rho(1450)", 1.370 )),
      _md("Width_rho(1450)+", _md("Width_rho(1450)", 0.510 )),
      running&2 
      );
  double help_beta = _md("beta_rho(1450)+", _md("beta_rho(1450)", -0.145) );
  m_A.InitialiseThreeBodyResonance(help_rho, help_rhoP, help_beta);
  m_AA.InitialiseThreeBodyResonance(help_rho, help_rhoP, help_beta);

  double mV[2], mVV[2], wV[2], wVV[2];
  for( int i=0; i<2; ++i ) {
	if( resV[i] == kf_rho_770 || resV[i] == kf_rho_770_plus ) {
	  mV[i]       = _md("Mass_"+Flavour(resV[i]).IDName(),  0.773 );       // mass^2 of vector resonance ij
	  wV[i]       = _md("Width_"+Flavour(resV[i]).IDName(), 0.145 );        // width^2 of vector resonance ij
	  mVV[i]      = _md("Mass_"+Flavour(resVV[i]).IDName(), 1.370 );      // mass^2 of vector resonance' ij
	  wVV[i]      = _md("Width_"+Flavour(resVV[i]).IDName(),0.510 );       // width^2 of vector resonance' ij
	  m_Beta[i]   = _md("beta_"+Flavour(resVV[i]).IDName(),  -0.145 );                            // weight factor for Vij'
	}
	else if( resV[i] == kf_K_star_892 || resV[i] == kf_K_star_892_plus ) {
	  mV[i]       = _md("Mass_"+Flavour(resV[i]).IDName(),  0.8921 );       // mass^2 of vector resonance ij
	  wV[i]       = _md("Width_"+Flavour(resV[i]).IDName(), 0.0513 );        // width^2 of vector resonance ij
	  mVV[i]      = _md("Mass_"+Flavour(resVV[i]).IDName(), 1.412 );      // mass^2 of vector resonance' ij
	  wVV[i]      = _md("Width_"+Flavour(resVV[i]).IDName(),0.227 );       // width^2 of vector resonance' ij
	  m_Beta[i]   = _md("beta_"+Flavour(resVV[i]).IDName(),  -0.135 );                            // weight factor for Vij'
	}
	else {
	  mV[i]       = _md("Mass_"+Flavour(resV[i]).IDName(),  Flavour(resV[i]).HadMass());       // mass^2 of vector resonance ij
	  wV[i]       = _md("Width_"+Flavour(resV[i]).IDName(), Flavour(resV[i]).Width());        // width^2 of vector resonance ij
	  mVV[i]      = _md("Mass_"+Flavour(resVV[i]).IDName(), Flavour(resVV[i]).HadMass());      // mass^2 of vector resonance' ij
	  wVV[i]      = _md("Width_"+Flavour(resVV[i]).IDName(),Flavour(resVV[i]).Width());       // width^2 of vector resonance' ij
	  m_Beta[i]   = _md("beta_"+Flavour(resVV[i]).IDName(),  0. );                            // weight factor for Vij'
	}
	m_V[i]      = ResonanceFlavour( resV[i], mV[i], wV[i], running&2 );
	m_VV[i]     = ResonanceFlavour( resVV[i], mVV[i], wVV[i], running&2 );
  }

  m_fpi2      = sqr( _md("fpi", 0.1307) );          // pion decay constant
}


VA_0_PPP::RChT::RChT(int mode, int kaon_mode, GeneralModel _md, double * _ms)
  : FF_Base(mode,kaon_mode,_md,_ms)
{
  // set resonances and parameters
  if( m_mode != 1200 && 
      m_mode != 3000 && 
      m_mode != 1020 ) {
    msg_Error()<<"Error: The mode "<<m_mode<<endl
	       <<"     hasn't been implemented yet (RChT). Please use KS model."
	       <<"     Don't know what to do. Will abort"<<endl;
    Abort();
  }
  
  m_MO     = _md("Mass_omega(782)", Flavour(kf_omega_782).HadMass());
  m_MO2    = sqr(m_MO);
  m_GO     = _md("Width_omega(782)", Flavour(kf_omega_782).Width());
  
  m_l0     = _md("lambda0", 1.);                    // fit parameter lambda0
  m_gammaR = _md("gamma_rho(770)", 1.);              // global factor for rho width
  m_m      = Flavour( kf_pi_plus ).HadMass();         // pion mass
  m_m2     = sqr(m_m);                              // pion mass^2
  m_mK2    = sqr( Flavour( kf_K_plus ).HadMass() );   // Kaon mass^2
  m_exp_alpha = _md("exp_alpha", 2.45);             // exponent in off-shell GA
  m_l1     = _md("lambda1", 0.5);                   // fit parameter        
  m_l2     = _md("lambda2", 0.);                    // fit parameter
  m_lsum   = m_l1 + m_l2; 
  // constraints due to short-distance behaviour
  m_FV2    = 2*m_fpi2;                              // vector coupling
  m_FV     = sqrt(m_FV2);
  m_GV     = m_FV/2.;
  m_FA2    = m_FV2 - m_fpi2;                        // axial coupling
  m_FA     = sqrt(m_FA2);
}

double VA_0_PPP::RChT::MassWidthVector( int a, double s )
{ 
  if(m_V[0].Running()) {
    switch(m_mode) {
    case 3000:
    case 1200: return MassWidthVector( s );
    case 1020: return (a==0)? MassWidthVector(s) : m_V[1].MassWidth();
    }
    // default:
    msg_Error()<<"Warning: this form factor (RChT) for the three-pseudoe mode "<<m_mode<<"\n"
	       <<"     hasn't been implemented yet. Please use KS model."<<endl;
  }
  return( m_V[a].MassWidth() );
}

double VA_0_PPP::RChT::MassWidthVector( double s )
{ 
  double MVGV (0.);
  if( s>4.*m_m2 )  MVGV += pow( 1.-4.*m_m2/s, 1.5 );
  if( s>4.*m_mK2 ) MVGV += pow( 1.-4.*m_mK2/s, 1.5 ) / 2.;
  MVGV *= m_gammaR*m_V[0].Mass2()*s/(96.*M_PI*m_fpi2); 
  return MVGV;
}

double VA_0_PPP::RChT::MassWidthAxial( double Q2 )
{
  if( !m_deltas && m_A.Running() )
    return(  m_A.OffShellMassWidth(Q2)*pow(m_A.Mass2()/Q2,m_exp_alpha-2.) );
  return m_A.MassWidth();
}

Complex VA_0_PPP::RChT::FormFactor( int j, double Q2, double s, double t )
{
  switch( m_mode ) {
  case 1200: 
  case 3000: { // 3pion mode
    if (j==1 || j==2) {        // axial contributions
      double u = Q2-s-t+Mass2(1-1)+Mass2(2-1)+Mass2(3-1);
      double x = (j==1)? s : t;
      double y = (j==1)? t : s;
      double MVGV_x = MassWidthVector(x);
      double MVGV_y = MassWidthVector(y);
      double MAGA = MassWidthAxial(Q2);
      double F_Q2_x = x/2./Q2 - m_l0*m_m2/Q2;
      double F_Q2_y = y/2./Q2 - m_l0*m_m2/Q2;
      double MV2 = m_V[0].Mass2();
      Complex alpha = 1. - 3./2.*x/Complex(x-MV2, MVGV_x);
      Complex beta =  -3./2.*x/Complex(x-MV2, MVGV_x)
	+ F_Q2_x*(2.*Q2+x-u)/Complex(x-MV2, MVGV_x)
	+ F_Q2_y*(u-x)/Complex(y-MV2, MVGV_y);
      return alpha - Q2/Complex(Q2-m_A.Mass2(),MAGA)*beta;
    }
    else {                     // pseudoscalar, vector
      return Complex(0.,0.);
    }
  }
  case 1020: { // K- pi- K+ mode
    double u = Q2-s-t+Mass2(1-1)+Mass2(2-1)+Mass2(3-1);
    switch(j) {
    case 1:
    case 2: {                // axial contributions
      double x = (j==1)? s : t;
      double y = (j==1)? t : s;
      int    a = j-1;                          // resonance assoc. with x
      int    b = (a==1)? 0 : 1;                 // resonance assoc. with y
      Complex FAchi = Complex(1.,0.);
      Complex FA1r  = 0.5 * m_FV*m_GV/m_fpi2 * ( 
						1./Complex(m_V[a].Mass2()-x,-1.*MassWidthVector(a,x))*
						( 3.*x + Mass2(2)-Mass2(a) + (1.-2.*m_GV/m_FV)*(2.*Q2-2.*x-u+Mass2(a)-Mass2(b)) )
						+1./Complex(m_V[b].Mass2()-y,-1.*MassWidthVector(b,y))*
                  ( 2.*(Mass2(b)-Mass2(2)) + (1.-2.*m_GV/m_FV)*(u-x+Mass2(2)-Mass2(b)) )    
                                                       );
              Complex FA2r  = -sqrt(2.) * m_FA*m_GV/m_fpi2 * Q2/Complex(m_A.Mass2()-Q2,-1.*MassWidthAxial(Q2)) * (
                  1./Complex(m_V[a].Mass2()-x,-1.*MassWidthVector(a,x))*
                      ( m_lsum*(-3.*x+Mass2(a)-Mass2(2)) + FFunc(x,Mass2(b),Q2)*(2.*Q2+x-u+Mass2(2)-Mass2(b)) )
                      +1./Complex(m_V[b].Mass2()-y,-1.*MassWidthVector(b,y))*
                      ( 2.*m_lsum*(Mass2(2)+Mass2(b)) + FFunc(y,Mass2(a),Q2)*(u-x+Mass2(b)-Mass2(2)) )
                                                                                );
              return FAchi + FA1r + FA2r;
            }
            case 3: {                // pseudoscalar contribution
              Complex FSchi = 3./2.* Mass2(1)/(Q2-Mass2(1)) *
                  ( 1. + (Mass2(2)-u)/Q2 );
              Complex FS1r  = 3./2.*sqr(m_GV)/m_fpi2 * Mass2(1)/(Q2*(Q2-Mass2(1))) *
                  (   s*(t-u)/Complex(m_V[0].Mass2()-s,-1.*MassWidthVector(0,s))
                      + (t*(s-u)+(Q2-Mass2(0))*(Mass2(1)-Mass2(2)))/Complex(m_V[1].Mass2()-t,-1.*MassWidthVector(1,t)) 
                  );
              return FSchi + FS1r;
            }
            case 4: {                // vector contribution
              double c1c2c5 = 0.;
//              double c1c2c52c6 = -3./(32.*sqr(M_PI)) * m_V[0].Mass()/(sqrt(2.)*m_FV);
              double c1c2c52c6 = -3./(96.*sqr(M_PI)) * m_V[0].Mass()/(sqrt(2.)*m_GV);
              double c1c28c3c5 = 0.;
              double c4 = 0.;
//              double d3 = -3./(64.*sqr(M_PI)) * m_V[0].Mass2()/m_FV2 + m_fpi2/(8.*m_FV2);
              double d3 = -3./(192.*sqr(M_PI)) * m_V[0].Mass2()/(m_FV*m_GV);
              double d18d2d3 = m_fpi2/(8.*m_FV2);
              double g12g2g3 = 0.;
              double g2 = 3./(192.*sqr(M_PI)) * m_V[0].Mass()/(sqrt(2.)*m_FV);
             
              Complex FVchi  = 3./(2.*sqr(M_PI)*m_fpi2); 
              Complex FV1r1  = 6.*m_GV/(sqrt(2.)*m_fpi2*m_V[0].Mass()) * (
                  1./Complex(m_MO2-s,-1.*m_MO*m_GO) * (
                    c1c2c5*Q2 - c1c2c52c6*s + c1c28c3c5*Mass2(1) ) +
                  1./Complex(m_V[1].Mass2()-t,-1.*MassWidthVector(1,t)) * (
                    c1c2c5*Q2 - c1c2c52c6*t + c1c28c3c5*Mass2(0) + 8.*c4*(Mass2(0)-Mass2(1)) )
                  );
              Complex FV1r2  = -12.*m_FV/(sqrt(2.)*m_V[0].Mass()*m_fpi2) * 
                1./Complex(m_V[0].Mass2()-Q2,-1*MassWidthVector(0,Q2)) * (
                  g12g2g3*(s+t)-2.*g2*Q2 );
              Complex FV2r   = -6.*m_FV*m_GV/m_fpi2 * 
                1./Complex(m_V[0].Mass2()-Q2,-1.*MassWidthVector(0,Q2)) * (
                    1./Complex(m_MO2-s,-1.*m_MO*m_GO) * (
                      d3*(Q2+s) + d18d2d3*Mass2(1) ) +
                    1./Complex(m_V[1].Mass2()-t,-1.*MassWidthVector(1,t)) * (
                      d3*(Q2+t) + d18d2d3*Mass2(0) )
                    );
              return FVchi + FV1r1 + FV1r2 + FV2r;
            }
        }
      }
  }
  return  (j>3)? Complex(0.,0.) : Complex(1.,0.); 
}

double VA_0_PPP::RChT::FFunc( double a, double b, double c)
{
  return m_l2 + m_l1*a/c - m_l0*b/c;
}


// Parameterisation
// DECKER, FINKEMEIER, MIRKES hep-ph/9310270

VA_0_PPP::KS::KS(int mode, int kaon_mode, GeneralModel _md, double * _ms)
  : FF_Base(mode,kaon_mode,_md,_ms)
{
  bool anomaly (false);
   
  // special constants for this particular parameterisation
  switch( m_mode ) {
    case 3000 : /* pi- pi- pi+ mode */
      m_X123  = 2.*Mass2(0);
      m_ms123 = Mass2(0);
      m_G123  = 1;
      break;
    case 1200 : /* pi0 pi0 pi- mode */
      m_X123  = Mass2(2);
      m_ms123 = Mass2(2);
      m_G123  = 1;
      break;
    case 1020 : /* K- pi- K+ */
      m_X123  = Mass2(1) + Mass2(0);
      m_ms123 = Mass2(1);
      m_G123  = 1;
      anomaly = true;
      break;
    case 1002 : /* K0 pi- K0b */
      m_X123  = Mass2(1) + Mass2(0);
      m_ms123 = Mass2(1);
      m_G123  = 1;
      anomaly = true;
      break;
    case  111 : /* K- pi0 K0 */
      m_X123  = 0.;
      m_ms123 = Mass2(1);
      m_G123  = 0;
      break;
    case  210 : /* pi0 pi0 K- mode */
      m_X123  = -2.*(Mass2(1)+Mass2(2));
      m_ms123 = Mass2(2);
      m_G123  = 1;
      break;
    case 2010 : /* K- pi- pi+ */
      m_X123  = Mass2(1)+Mass2(0);
      m_ms123 = Mass2(0);
      m_G123  = 1;
      anomaly = true;
      break;
    case 1101 : /* pi- K0b pi0 */
      m_X123  = 0.;
      m_ms123 = Mass2(1);
      m_G123  = 0;
      anomaly = true;
      break;
    case   30 : /* K- K- K+ */
      m_X123  = 2.*Mass2(2);
      m_ms123 = Mass2(2);
      m_G123  = 1;
      break;
    case   12 : /* K- K0b K0 */
      m_X123  = 2.*Mass2(0);
      m_ms123 = Mass2(0);
      m_G123  = 1;
      break;
    default:
      msg_Error()<<METHOD<<" Warning: Three pseudoscalar mode not recognized. "
        <<"mode="<<m_mode<<endl;
      break;
  }
   
  kf_code resAnoV   = (m_deltas)? kf_K_star_892_plus : kf_rho_770_plus;
  kf_code resAnoVV  = (m_deltas)? kf_K_star_1410_plus : kf_rho_1450_plus;
  kf_code resAnoVVV = (m_deltas)? kf_K_star_1680_plus : kf_rho_1700_plus;

  if( anomaly ) {
	int running    = int( _md("RUNNING_WIDTH", 3 ) );     // running width
    double MV   = _md("Mass_anomaly_"+Flavour(resAnoV).IDName(),   Flavour(resAnoV).HadMass()   );   // mass V
    double MVV  = _md("Mass_anomaly_"+Flavour(resAnoVV).IDName(),  Flavour(resAnoVV).HadMass()   );  // mass V'
    double MVVV = _md("Mass_anomaly_"+Flavour(resAnoVVV).IDName(), Flavour(resAnoVVV).HadMass()   ); // mass V'
    double GV   = _md("Width_anomaly_"+Flavour(resAnoV).IDName(),   Flavour(resAnoV).Width()   );   // width V
    double GVV  = _md("Width_anomaly_"+Flavour(resAnoVV).IDName(),  Flavour(resAnoVV).Width()   );  // width V'
    double GVVV = _md("Width_anomaly_"+Flavour(resAnoVVV).IDName(), Flavour(resAnoVVV).Width()   ); // width V'
    m_AnoV      = ResonanceFlavour( resAnoV, MV, GV, running&2 );
    m_AnoVV     = ResonanceFlavour( resAnoVV, MVV, GVV, running&2 );
    m_AnoVVV    = ResonanceFlavour( resAnoVVV, MVVV, GVVV, running&2 );
    m_AlphaV    = _md("alpha_anomaly_K*(892)", _md("alpha_anomaly_K*(892)+", 0. ));   // alpha for K*
    m_BetaV[0]  = _md("beta_anomaly_"+Flavour(resAnoVV).IDName(), 0. );        // beta for V'
    m_BetaV[1]  = _md("gamma_anomaly_"+Flavour(resAnoVVV).IDName(), 0. );      // gamma for V''
  }
  else {
    m_AnoV      = ResonanceFlavour( kf_rho_770_plus, 0., 0., 0 );
    m_AnoVV     = ResonanceFlavour( kf_rho_770_plus, 0., 0., 0 );
    m_AnoVVV    = ResonanceFlavour( kf_rho_770_plus, 0., 0., 0 );
    m_AlphaV    = 0.;
    m_BetaV[0]  = 0.;
    m_BetaV[1]  = 0.;
  }

} 
 
// methods for axial and scalar FF
 
Complex VA_0_PPP::KS::BW_V( int a, double s )
{
  if( !m_G123 && a==1 ) return Complex(0.,0.);          // BW_V23 = 0 if G123=0
  return m_V[a].BreitWigner(s);
}
 
Complex VA_0_PPP::KS::BW_VV( int a, double s )
{
  if( !m_G123 && a==1 ) return Complex(0.,0.);          // BW_V23' = 0 if G123=0
  return m_VV[a].BreitWigner(s);
}

Complex VA_0_PPP::KS::BW_A( double s )
{
  if (!IsZero(m_alpha)) {
    return( ( m_A.BreitWigner(s)+m_alpha*m_AA.BreitWigner(s) ) / (1.+m_alpha) );
  }
  return m_A.BreitWigner(s);
}


Complex VA_0_PPP::KS::Tvector1( int a, int b, double x )
{
  Complex ret =           
                BW_V(a,x) *( 1.-(Mass2(a)-Mass2(b))/(3.*m_V[a].Mass2()) )
  + m_Beta[a]*( BW_VV(a,x)*( 1.-(Mass2(a)-Mass2(b))/(3.*m_VV[a].Mass2()) ) );
  return ret/(1.+m_Beta[a]);
}

Complex VA_0_PPP::KS::Tvector2( int a, int b, double x )
{
  Complex ret = BW_V(a,x)/m_V[a].Mass2() + m_Beta[a]*BW_VV(a,x)/m_V[a].Mass2();
  return ret*( 2.*(Mass2(a)-Mass2(b)) )/( 3.*(1.+m_Beta[a]) );    
}

Complex VA_0_PPP::KS::TSvector( int a, int b, int c, double Q2, double s, double t )
{
  Complex ret =
                BW_V(a,s)  * (   m_ms123*( Q2-2.*t-s+2.*Mass2(a)+Mass2(c) )
                               - (Mass2(a)-Mass2(b))/m_V[a].Mass2()*( m_ms123*(Q2+s-Mass2(c))-Q2*(s-m_V[a].Mass2()) ) )
    + m_Beta[a]*BW_VV(a,s) * (   m_ms123*( Q2-2.*t-s+2.*Mass2(a)+Mass2(c) ) 
                               - (Mass2(a)-Mass2(b))/m_VV[a].Mass2()*( m_ms123*(Q2+s-Mass2(c))-Q2*(s-m_VV[a].Mass2()) ) );
  return ret/(1.+m_Beta[a]);
}

Complex VA_0_PPP::KS::Tgen( int a, int b, int c, double s, double t)
{
  return( Tvector1(a-1,c-1,s) + Tvector2(b-1,c-1,t) );
}

// methods for vector FF
 
Complex VA_0_PPP::KS::T_V( double q2 )
{
  Complex ret(0.,0.);
  ret  =           m_AnoV.BreitWigner(q2)
    + m_BetaV[0] * m_AnoVV.BreitWigner(q2)
    + m_BetaV[1] * m_AnoVVV.BreitWigner(q2);
  ret /= 1.+m_BetaV[0]+m_BetaV[1];
  return ret;
}

Complex VA_0_PPP::KS::T_V13V23( double s, double t )
{
  Complex ret(0.,0.);
  int n_rho (1-1), n_k (2-1);                       // what is the rho and K* resonance
  double x (s), y (t);
  if( m_mode==2010 ) { 
    n_rho = 2-1; n_k = 1-1;                         // swap them in Kpipi channel
    x=t; y=s;
  }    
  ret  = ( BW_V(n_rho,x) + m_Beta[n_rho]*BW_VV(n_rho,x) )/( 1.+m_Beta[n_rho] );
  ret += m_AlphaV * BW_V(n_k,y);
  ret /= 1.+m_AlphaV;
  return ret;
}

// handling of form factors

Complex VA_0_PPP::KS::FormFactor( int j, double Q2, double s, double t )
{
  Complex FF(0.,0.);
  switch( j ) {
    case 1 : { FF = BW_A(Q2)*Tgen(1,2,3,s,t);                   // axial
             break; }
    case 2 : { FF = BW_A(Q2)*Tgen(2,1,3,t,s);                   // axial
             break; }
    case 3 : { FF = m_X123 + BW_A(Q2)*(Q2-m_A.Mass2())/(m_A.Mass2()*Q2)     // scalar
                         *( TSvector(1-1,3-1,2-1,Q2,s,t) + TSvector(2-1,3-1,1-1,Q2,t,s) );
             FF /= 2.*(Q2-m_ms123);          
             break; }
    case 4 : { FF = 3./(8.*sqr(M_PI)*m_fpi2)*T_V(Q2)*T_V13V23(s,t); // vector
               break; }
  }
  return FF;
}

// Parameterisation
// FINKEMEIER, MIRKES hep-ph/9503474

VA_0_PPP::KS95::KS95(int mode, int kaon_mode, GeneralModel _md, double * _ms)
  : FF_Base(mode,kaon_mode,_md,_ms)
{
  kf_code resV, resVV, resVVV;
  int running    = int( _md("RUNNING_WIDTH", 3 ) );     // running width
  double MV, MVV, MVVV, GV, GVV, GVVV;
  if( m_deltas==0 ) {
	resV   = kf_rho_770_plus;
	resVV  = kf_rho_1450_plus;
	resVVV = kf_rho_1700_plus;
	MV   = _md("q2_Mass_"+Flavour(resV).IDName(),   0.7769 );   // mass V
	MVV  = _md("q2_Mass_"+Flavour(resVV).IDName(),  1.363 );  // mass V'
	MVVV = _md("q2_Mass_"+Flavour(resVVV).IDName(), 1.700 ); // mass V"
	GV   = _md("q2_Width_"+Flavour(resV).IDName(),   0.1490 );   // width V
	GVV  = _md("q2_Width_"+Flavour(resVV).IDName(),  0.310 );  // width V'
	GVVV = _md("q2_Width_"+Flavour(resVVV).IDName(), 0.235 ); // width V"
	m_beta_VectorFF   = _md("q2_beta_"+Flavour(resVV).IDName(), -0.108 );        // beta for V'
	m_gamma_VectorFF  = _md("q2_gamma_"+Flavour(resVVV).IDName(), 0.020 );      // gamma for V"
  } 
  else {
	resV   = kf_K_star_892_plus;
	resVV  = kf_K_star_1410_plus;
	resVVV = kf_K_star_1680_plus;
	MV   = _md("q2_Mass_"+Flavour(resV).IDName(),   0.892    );   // mass V
	MVV  = _md("q2_Mass_"+Flavour(resVV).IDName(),  1.412    );  // mass V'
	MVVV = _md("q2_Mass_"+Flavour(resVVV).IDName(), 1.714    ); // mass V"
	GV   = _md("q2_Width_"+Flavour(resV).IDName(),   0.051   );   // width V
	GVV  = _md("q2_Width_"+Flavour(resVV).IDName(),  0.227   );  // width V'
	GVVV = _md("q2_Width_"+Flavour(resVVV).IDName(), 0.323   ); // width V"
	m_beta_VectorFF   = _md("q2_beta_"+Flavour(resVV).IDName(), -0.25 );        // beta for V'
	m_gamma_VectorFF  = _md("q2_gamma_"+Flavour(resVVV).IDName(), 0.038 );      // gamma for V"
  }

  m_V_VectorFF      = ResonanceFlavour( resV, MV, GV, running&2 );
  m_VV_VectorFF     = ResonanceFlavour( resVV, MVV, GVV, running&2 );
  m_VVV_VectorFF    = ResonanceFlavour( resVVV, MVVV, GVVV, running&2 );

  // omega - phi - resonances
  kf_code resO = kf_omega_782;
  kf_code resP = kf_phi_1020;
  double MO   = _md("G3_Mass_"+Flavour(resO).IDName(),   0.782   );    // mass omega
  double MP   = _md("G3_Mass_"+Flavour(resP).IDName(),   1.020   );    // mass phi
  double GO   = _md("G3_Width_"+Flavour(resO).IDName(),   0.00843   );   // mass omega
  double GP   = _md("G3_Width_"+Flavour(resP).IDName(),   0.00443   );   // mass phi
  m_Omega_VectorFF	= ResonanceFlavour( resO, MO, GO, 0 );
  m_Phi_VectorFF	= ResonanceFlavour( resP, MP, GP, 0 );
  m_eps_VectorFF    = _md("G3_epsilon_"+Flavour(resP).IDName(), 0.05 );        			// epsilon for phi

} 
 
Complex VA_0_PPP::KS95::Axial(double s)
{
  if (!IsZero(m_alpha)) {
    return( ( m_A.BreitWigner(s)+m_alpha*m_AA.BreitWigner(s) ) / (1.+m_alpha) );
  }
  return m_A.BreitWigner(s);
}
 
Complex VA_0_PPP::KS95::TVector(int i, double s)
{
  return (m_V[i].BreitWigner(s)+m_Beta[i]*m_VV[i].BreitWigner(s))/(1.+m_Beta[i]);
}
 
Complex VA_0_PPP::KS95::Vector(int i, double s, double u)
{
  switch(m_mode) {
	case 1101 : { 	// pi- K0b pi0
				  if(m_kaon_mode==3) {			// K0 K0
					switch(i) {
					  case 0 : return TVector(0,s);
					  case 1 : return TVector(1,s);
					}
				  }
				  else if(m_kaon_mode==20) { 	// KS KS
					switch(i) {
					  case 0 : return 0.5*TVector(1,u);
					  case 1 : return -0.5*(TVector(1,s)+TVector(1,u));
					}
				  }
				  else if(m_kaon_mode==200) {	// KL KL
					switch(i) {
					  case 0 : return 0.5*TVector(1,u);
					  case 1 : return -0.5*(TVector(1,s)+TVector(1,u));
					}
					switch(i) {
					  case 0 : return -0.5*TVector(1,u);
					  case 1 : return 0.5*(TVector(1,s)+TVector(1,u));
					}
				  }
				  else if(m_kaon_mode==110) {	// KS KL
					switch(i) {
					  case 0 : return 0.5*(2.*TVector(0,s)+TVector(1,u));
					  case 1 : return 0.5*(TVector(1,s)-TVector(1,u));
					}
				  }
				}
	case  111 : {	// K- pi0 K0
				  switch(i) {
					case 0 : return 2./3.*TVector(0,s) + 1./3.*TVector(1,u);
					case 1 : return 1./3.*TVector(1,s) - 1./3.*TVector(1,u);
				  }
				}
  }
  return TVector(i,s);
}
 
Complex VA_0_PPP::KS95::TOmega(double s)
{
  return (m_Omega_VectorFF.BreitWigner(s) + m_eps_VectorFF*m_Phi_VectorFF.BreitWigner(s)) / (1.+m_eps_VectorFF);
}
 
Complex VA_0_PPP::KS95::G3(double q2, double s, double t, double u )
{
  Complex vector1 
	=                     m_V_VectorFF.BreitWigner(q2)
	   +  m_beta_VectorFF*m_VV_VectorFF.BreitWigner(q2)
	   + m_gamma_VectorFF*m_VVV_VectorFF.BreitWigner(q2);
  vector1 /= 1.+m_beta_VectorFF+m_gamma_VectorFF;	   
	    
  Complex vector(0.,0.);
  switch(m_mode) {
	case 3000 :
	case 1200 : return Complex(0.,0.);
	case  210 : // pi0 pi0 K-
				{
				  vector  = TVector(2-1,t) - TVector(1-1,s);
				}
	case 2010 : // K- pi- pi+
				{
				  vector  = TVector(2-1,t) + TVector(1-1,s);
				}
	case 1101 : // pi- K0b pi0
				{
				  vector  = 2.*TVector(1-1,s) + TVector(2-1,t) + TVector(2-1,u);
				  vector /= 3.;
				}
	case 1020 : // K- pi- K+
				{
				  vector  = sqrt(2.)*TOmega(s) + TVector(2-1,t);
				  vector *= (sqrt(2.)-1.);
				}
	case 1002 : // K0 pi- K0b
				{
				  if(m_kaon_mode==2) {	// K0 K0
					vector  = sqrt(2.)*TOmega(s) + TVector(2-1,t);
					vector *= (sqrt(2.)-1.);
				  }
				  else if(m_kaon_mode==20) {	// KS KS
					vector  = TVector(2-1,t)-TVector(2-1,u);
					vector *= -0.5*(sqrt(2.)-1.);
				  }
				  else if(m_kaon_mode==200) {	// KL KL
					vector  = TVector(2-1,t)-TVector(2-1,u);
					vector *= -0.5*(sqrt(2.)-1.);
				  }
				  else if(m_kaon_mode==110)	{	// KS KL
					vector  = 2.*sqrt(2.)*TOmega(s)+TVector(2-1,t)+TVector(2-1,u);
					vector *= 0.5*(sqrt(2.)-1.);
				  }
				}
	case  111 : // K- pi0 K0
				{
				  vector  = TVector(1-1,s) - TVector(2-1,t);
				  vector *= (sqrt(2.)-1.);
				}
  }
  return vector1*vector;
}

// handling of form factors

Complex VA_0_PPP::KS95::FormFactor( int j, double Q2, double s, double t )
{
  double u = Q2 + Mass2(1-1) + Mass2(2-1) + Mass2(3-1) - s - t;
  switch( j ) {
    case 1 : 	// axial
		return Axial(Q2)*Vector(1-1,s,u);
    case 2 : 	// axial
		return Axial(Q2)*Vector(2-1,t,u);
    case 3 : 	// scalar
		break;
    case 4 :  	// vector
		return 6./(8.*sqr(M_PI)*m_fpi2)*G3(Q2,s,t,u);
		break; 
  }
  Complex FF(0.,0.);
  return FF;
}

DEFINE_CURRENT_GETTER(HADRONS::VA_0_PPP,"VA_0_PPP")

void ATOOLS::Getter<HADRONS::Current_Base,HADRONS::ME_Parameters,HADRONS::VA_0_PPP>::
PrintInfo(std::ostream &st,const size_t width) const {
  st<<"Example: $ 0 \\rightarrow \\pi \\pi K $ \n\n"
    <<"Order: special, see tau decays. \n\n"
    <<"Available form factors: \n "
    <<"  \\begin{itemize} \n"
    <<"    \\item {\\tt FORM\\_FACTOR = 1 :} Kuehn-Santamaria \n"
    <<"    \\item {\\tt FORM\\_FACTOR = 2 :} Resonance Chiral Theory \n"
    <<"  \\end{itemize} \n"
    <<"Reference: https://sherpa.hepforge.org/olddokuwiki/data/media/publications/theses/diplom\\__laubrich.pdf \n"
    <<std::endl;
}