#include "sigma_zz.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;
}