#include <TMath.h>
#include <kin_funcs.h>

double BH_cros_section( double *x, double *par )
{
  const double radian = 57.2957795130823229;
  double phi = x[0]/radian;
  double theta = x[1]/radian;
  
  double t = par[0];
  double Q2 = par[1];
  double Eg = par[2];

  double weight;
  
  const double m_e = 0.00051;
  const double M_p = 0.938272;
  const double alpha_em = 1./137.;
  const double PI = 3.14159265358979312;
  const double ammp = 2.793;
  const double ammn = -1.913;

  //int sigma = 1;

  double s = M_p*( M_p + 2*Eg );
  double beta = sqrt( 1 - (4*m_e*m_e)/Q2 );
  double r = sqrt( (s - Q2 - M_p*M_p)*(s - Q2 - M_p*M_p) - 4*Q2*M_p*M_p );
  double tau = Q2/( s - M_p*M_p );
  double cos_TH_Cm = (2*s*(t - 2*M_p*M_p) + (s + M_p*M_p)*(s + M_p*M_p - Q2))/sqrt(Lambda(s, M_p*M_p, 0)*Lambda(s, M_p*M_p, Q2));
  double sin_TH_Cm = sqrt( 1 - cos_TH_Cm*cos_TH_Cm );
  double Delta_Perp = sin_TH_Cm*r/(2*sqrt(s));
  //double Delta_Perp = sqrt((-t)*(1 - tau) - tau*tau*M_p*M_p );
  double a = beta*r*cos(theta);
  // double b = sigma*beta*sqrt( (Q2-t)*(Q2-t) - TMath::Power((2*(s-M_p*M_p)*sqrt(Q2)*Delta_Perp)/r , 2) )*cos(theta) - 
  //   beta*( (2*(s-M_p*M_p)*sqrt(Q2)*Delta_Perp)/r )*sin(theta)*cos(phi);

  double b = beta*((Q2*(s - M_p*M_p - Q2) + t*(s - M_p*M_p + Q2))/r)*cos(theta) -
    beta*( (2*(s-M_p*M_p)*sqrt(Q2)*Delta_Perp)/r )*sin(theta)*cos(phi);

  double L_BH = ((Q2 - t)*(Q2 - t) - b*b)/4.;
  double L0_BH = Q2*Q2*sin(theta)*sin(theta)/4.;
  double A_BH = (s - M_p*M_p)*(s - M_p*M_p)*Delta_Perp*Delta_Perp - t*a*(a + b) - M_p*M_p*b*b - t*(4*M_p*M_p - t)*Q2 +
    (m_e*m_e/L_BH)*(  TMath::Power( (Q2 - t)*(a + b) - (s - M_p*M_p)*b, 2) + t*(4*M_p*M_p - t)*(Q2 - t)*(Q2 - t) );
  
  double B_BH = (Q2 + t)*(Q2 + t) + b*b + 8*m_e*m_e*Q2 - 4*m_e*m_e*(t + 2*m_e*m_e)*(Q2 - t)*(Q2 - t)/L_BH;
    
  double F1p = (1./((1 - t/0.71)*(1 - t/0.71)))*(1/(1 - t/(4.*M_p*M_p) ))*( 1 - 2.79*t/( 4*M_p*M_p ) );
  double F2p = (1./((1 - t/0.71)*(1 - t/0.71)))*(1/(1 - t/(4.*M_p*M_p) ))*(ammp - 1);
  //double F1n = (1./((1 - t/0.71)*(1 - t/0.71)))*(1/(1 - t/(4.*M_p*M_p) ))*(-ammn*t/(4*M_p*M_p));
  //double F2n = (1./((1 - t/0.71)*(1 - t/0.71)))*(1/(1 - t/(4.*M_p*M_p) ))*ammn;

  if( par[3] == 1 )
    { weight = L_BH/L0_BH; }
  else
    {
      weight = 1.;
    }
  
  double crs = weight*sin(theta)*(TMath::Power(alpha_em, 3)/(4*PI*(s - M_p*M_p)*(s - M_p*M_p)) )*(beta/(-t*L_BH))*
    ( (A_BH/(-t))*(F1p*F1p - (t/(4*M_p*M_p))*F2p*F2p ) + (F1p + F2p)*(F1p + F2p)*B_BH/2 )*
    0.389379*1e9;

  return crs;
}