#include "HADRONS++/Current_Library/VA_P_P.H"
#include <cmath>
// BGL parametrization -> arXiv:hep-ph/9508211

namespace HADRONS {
namespace VA_P_P_FFs {

class BGL : public FormFactor_Base {
 int case_number; // case_number selects the used phi- and f_plus function
 //Variables used by case 0 arXiv:1208.1253 [hep-ex] and case 2, arXiv:1510.03657 [hep-ex]
 // a0 is for m_fplus, a0_0 is for m_f0
 double mstar ; double a0; double a1; double a2; double a3; double a0_0; double a1_0; double a2_0; double a3_0; double q02;
 double fplus_0; double r1; double r2; // Variables used by case 1
 double zCalc(double q2, double t_0, double t_plus, double t_minus);

public:
  BGL(GeneralModel model, double* masses, const Flavour_Vector& flavs,
        std::vector<int>& indices);
  void CalcFFs(ATOOLS::Vec4D p0, ATOOLS::Vec4D p1);
};

BGL::BGL(GeneralModel model, double* masses, const Flavour_Vector& flavs,
             std::vector<int>& indices) :
  FormFactor_Base(model, masses, flavs, indices)
{
  kf_code kf0=m_flavs[p_i[0]].Kfcode();
  kf_code kf1=m_flavs[p_i[1]].Kfcode();
  // B-Decay:
  if (kf0==kf_B) {
    if (kf1==kf_pi_plus) {
    case_number = 0;
    const double mminus = m_m0-m_m1;
    mstar = 5.32483;
    a0 = 0.026;
    a1 = -1.15*a0;
    a2 = -4.52*a0;
    a3 = 0;
    q02 = 0.65*mminus*mminus;  
    }
    else if (kf1==kf_D_plus) {
      case_number = 2;
      a0 = 0.0126;
      a1 = -0.094;
      a2 = 0.34;
      a3 = -0.1;

      a0_0  = 0.0115;
      a1_0  = -0.057;
      a2_0  = 0.12;
      a3_0  = 0.4;
     }
  }
  // B+-Decay:
   else if (kf0==kf_B_plus) {
     if (kf1==kf_pi) {
       const double mminus = m_m0-m_m1;
       case_number = 0;
       mstar = 5.32483;
       a0 = 0.026;
       a1 = -0.63*a0;
       a2 = -5.80*a0;
       a3 = 0;
       q02 = 0.65*mminus*mminus;
    }
    else if (kf1==kf_eta) {
      case_number = 0;
      mstar=5.32483;
      a0 = 0.0031;
      a1 = -0.005301;
      q02 = 14.14;
    }
    else if (kf1==kf_D) {
      case_number = 2;
      a0 = 0.0126;
      a1 = -0.094;
      a2 = 0.34;
      a3 = -0.1;

      a0_0 = 0.0115;
      a1_0 = -0.057;
      a2_0 = 0.12;
      a3_0 = 0.4;
     }
  }
  // D0-Decay:
  else if (kf0==kf_D) {
    if (kf1==kf_pi_plus) {
      case_number = 1;
      fplus_0 = 0.6748956576;
      r1 = -2.0;
      r2 = -1.6;
    }
    else if (kf1==kf_K_plus) {
      case_number = 1;
      fplus_0 = 0.7653349496;  
      r1 = -2.4;
      r2 = 15.6;
    }
  }
  // D+-Decay
  else if (kf0==kf_D_plus) {
    if (kf1==kf_pi) {
      case_number = 1;
      fplus_0 = 0.6340467099; 
      r1 = -1.95;
      r2 = -0.11;
    }
    else if (kf1==kf_K) {
      case_number = 1;
      fplus_0 = 0.738;
      r1 = -1.9;
      r2 = 16.6;
    }
  }
  
  mstar = model("mstar",mstar);
  a0 = model("a0",a0);
  a1 = model("a1",a1);
  a2 = model("a2",a2);
  a3 = model("a3",a3); 
  a0_0 = model("a0_0",a0_0);
  a1_0 = model("a1_0",a1_0);
  a2_0 = model("a2_0",a2_0);
  a3_0 = model("a3_0",a3_0); 
  fplus_0 = model("fplus_0",fplus_0);
  r1 = model("r1",r1);
  r2 = model("r2",r2);
  case_number = model("case_number",case_number);
}

void BGL::CalcFFs( Vec4D p0, Vec4D p1 )
{
  if (case_number==0) {
    const double mplus = m_m0+m_m1;
    const double mminus = m_m0-m_m1;
    const double chi = pow(6.889,-4);
    double q2 = (p0-p1).Abs2();
    double z = (sqrt(mplus*mplus-q2)-sqrt(mplus*mplus-q02))/(sqrt(mplus*mplus-q2)+sqrt(mplus*mplus-q02));
    // Blasche factor P
    double P = (sqrt(mplus*mplus-q2)-sqrt(mplus*mplus-mstar*mstar ))/(sqrt(mplus*mplus-q2)+sqrt(mplus*mplus-q02));
    double phi = sqrt(1/(32*M_PI*chi))*(sqrt(mplus*mplus-q2)+sqrt(mplus*mplus-q02))*pow(sqrt(mplus*mplus-q2)+sqrt(mplus*mplus-(mminus*mminus)),1.5)*pow(sqrt(mplus*mplus-q2)+mplus,-5)*(mplus*mplus-q2)/pow(mplus*mplus-q02,0.25);  
    m_fplus  = 1/(P*phi) * (a0 + a1*z + a2*z*z + a3*z*z*z);
    m_f0     = 0.0;
    m_calced = true;
  }

  else if (case_number==1) {
    const double m_Dstar_plus = 2.1121;
    const double zero = 0.0;
    const double m_c = 1.29;

    kf_code kf0=m_flavs[p_i[0]].Kfcode();
    kf_code kf1=m_flavs[p_i[1]].Kfcode();

    double q2 = (p0 - p1).Abs2();
    double m2_c = pow(m_c,2);
    double t_plus = pow(m_m0 + m_m1,2);
    double t_minus = pow(m_m0 - m_m1,2);
    double t_0 = t_plus * (1 - sqrt(1 - t_minus/t_plus));
    double m2_Dstar_plus = pow(m_Dstar_plus,2);
    double c = sqrt(M_PI * m2_c/3);
    // c fuer Parameter von Paper: http://arxiv.org/pdf/1412.5502v1.pdf
    //   double c = sqrt(4 * M_PI * m2_c/9);
    double P_q2;
    double P_0;

    if (kf1==kf_pi || kf1==kf_pi_plus) { P_q2 = 1; }
    else if (kf1==kf_K || kf1==kf_K_plus) { P_q2 = zCalc(q2, m2_Dstar_plus, t_plus, t_minus); }

    if (kf1==kf_pi || kf1==kf_pi_plus) { P_0 = 1; }
    else if (kf1==kf_K || kf1==kf_K_plus) { P_0 = zCalc(zero, m2_Dstar_plus, t_plus, t_minus); }
   
    // Berechnung von Phi:
    double phi = c * pow(zCalc(q2, zero, t_plus, t_minus)/(-1 * q2),2.5) *
      pow(zCalc(q2, t_0, t_plus, t_minus)/(t_0 - q2),-0.5) *
      pow(zCalc(q2, t_minus, t_plus, t_minus)/(t_minus - q2),-0.75) *
      (t_plus - q2)/(pow(t_plus - t_0,0.25));
    // Berechnung von Phi(0):
    double phi_0 = c * pow(zCalc(zero, t_0, t_plus, t_minus)/(t_0 - zero),-0.5) *
      pow(zCalc(zero, t_minus, t_plus, t_minus)/(t_minus - zero),-0.75) *
      (t_plus - zero)/(pow(t_plus - t_0,0.25));
    // Berechnung von z:
    double z = zCalc(q2, t_0, t_plus, t_minus);
    // Berechnung von z(0)
    double z_0 = zCalc(zero, t_0, t_plus, t_minus);

    m_fplus = (fplus_0 * P_0 * phi_0 * (1 + r1*z + r2*z*z)) / (P_q2 * phi * (1 + r1*z_0 + r2*z_0*z_0));
    m_f0 = 0.0;
    m_calced = true;
  }

  else if (case_number==2) {
    double r = m_m1/m_m0;
    double q2 = (p0-p1).Abs2();
    double w = (m_m0*m_m0+m_m1*m_m1-q2)/(2*m_m0*m_m1);
    double z = (sqrt(w+1)-sqrt(2))/(sqrt(w+1)+sqrt(2));
    double phi_plus = 1.1213*pow(1+z,2)*pow(1-z,0.5)*pow((1+r)*(1-z)+2*sqrt(r)*(1+z),-5);  
    double phi_0 = 0.5299*(1+z)*pow(1-z,1.5)*pow((1+r)*(1-z)+2*sqrt(r)*(1+z),-4); 
    m_fplus  = 1/phi_plus*(a0+a1*z+a2*z*z+a3*z*z*z);
    m_f0     = 1/phi_0*(a0_0+a1_0*z+a2_0*z*z+a3_0*z*z*z);
    m_calced = true;
  }
}

double BGL::zCalc(double q2, double t_0, double t_plus, double t_minus)
{
  double z = (sqrt(t_plus - q2) - sqrt(t_plus - t_0))/(sqrt(t_plus - q2) + sqrt(t_plus - t_0));
  return z;
}

} // namespace VA_P_V
} // namespace HADRONS