#include "sigma_tot.h"
#include <math.h>

float sigma_zz(float Ecm, float mz){

  const float pi = 3.14159265;
  const float hbarc = 0.1973;           // GeV-fm
  const float sw2 = 0.23;
  const float G_mu =  1.16639e-5;

  float sigma;
  float ae = 1.;
  float ve = 1. - 4.*sw2;
  float s = pow(Ecm,2);
  float s1 = pow(mz,2);
  float s2 = pow(mz,2);
  float lambda = pow(s,2) + pow(s1,2) + pow(s2,2)
    - 2.*s*s1 - 2.*s1*s2 - 2.*s2*s;
  if ( lambda < 0 ) {
    sigma = 0.;
  }
  else {
    float L4 = (1./sqrt(lambda))  * log ( 
         ( s - s1 - s2 + sqrt(lambda) ) / 
         ( s - s1 - s2 - sqrt(lambda) ) );
    float G4_tpu = ( sqrt(lambda) / s ) * ( 
             ( ( pow(s,2) + pow((s1+s2),2) ) / (s - s1 - s2 ) ) * L4
             - 2. );
    sigma = ( (pow(G_mu,2) * pow(mz,4)) / (64.*pi*s) ) * 
              (pow(ve,4) + 6.*pow((ve*ae),2) + pow(ae,4)) * G4_tpu ;
    sigma = sigma * pow(hbarc,2) * 1.e10;    // convert to pb
  }
  return sigma;
}


// Cross section for e+e- -> HZ from Barger et al., Phys Rev D
// 49 (1994) 79, equation (2).

float sigma_hz(float Ecm, float mh, float mz){

  const float pi = 3.14159265;
  const float hbarc = 0.1973;           // GeV-fm
  const float sw2 = 0.23;
  const float G_mu =  1.16639e-5;

  float sigma;
  float ae = 1.;
  float ve = 1. - 4.*sw2;
  float s = pow(Ecm,2);
  float mz2 = pow(mz,2);
  float mh2 = pow(mh,2);
  float beta2 = ( 1 - pow((mh+mz),2)/s ) *
                ( 1 - pow((mh-mz),2)/s );
  float beta;
  if ( beta2 <= 0. ) {
    sigma = 0.;
  }
  else {
    sigma = ( ( pow(G_mu,2)*pow(mz,4) ) / ( 96.* pi * s ) ) *  
      ( pow(ae,2) + pow(ve,2) ) * sqrt(beta2) * 
      ( beta2 + 12.*mz2/s ) / pow((1.-mz2/s), 2);;
  }
  sigma = sigma * pow(hbarc,2) * 1.e10;    // convert to pb
  return sigma;
}