#include <iostream>
#include <fstream>
#include <math.h>
int ndata;
double xdata[500], xfdata[500];
double partonLuminosity(double tau);
void readinPDF();
double g(double x);
double I2(double x);
double alphaS(double Q);

int main()
{
double mtop = 172.0;
double ecm = 7000.;   // CM energy
double pi=3.141593;
//
readinPDF();
const double GF=1.16639e-5;
for(int i=0; i<90; i++) {
//------------------------------------------------------
double mH =100.+10.*i;
double tau = mH*mH/(ecm*ecm);
double c1 = GF*pow(alphaS(mH*mH),2)/(288.*1.414214*pi);
double crossGeVInv2 = c1*I2(mtop*mtop/(mH*mH))*tau*partonLuminosity(tau);
double crossSection = crossGeVInv2*197*197*1.e-6*100.*1.e12;
//   ****** THIS PART MUST BE CODED ******
  

//------------------------------------------------------
  std::cout <<mH<<"     "<<crossSection<<std::endl;
}
     
return 0;
}

//----------------------------------------------
//   calculate parton Luminosity by integration
//----------------------------------------------
double partonLuminosity(double tau) {
  int nint = 50000;
  double sum=0.;
  double deltax = (1.-tau)/nint;
  for (int i=0; i<nint; i++) {
    double x = tau + (i+0.5)*deltax;
    double d= deltax*g(x)*g(tau/x)/x;
    sum += d;
  }
return sum;
}

//----------------------------------------------
//   read-in gluon density from gluon_10000.txt
//----------------------------------------------
void readinPDF() {
  double dmy;
  std::ifstream fin("MSTW_10000.txt");
  ndata = 0;
  char line[255];
  fin.getline(line,sizeof(line)); // skipping the 1st dummy line
  while (fin >> xdata[ndata]>>dmy>>dmy>>dmy>>dmy
  >>dmy>>dmy>>dmy>>dmy>>xfdata[ndata]) ndata++;
fin.close();
std::cout<<"read in gluon density ended, ndata ="<< ndata+1 <<std::endl;
return;
}

//--------------------------------------------
//   calculate gluon density g(x) by interpolation
//--------------------------------------------
double g(double x){
if(x < xdata[0]) return 0.;
if(x > xdata[ndata-1] ) return 0.;
int i;
for ( i=1; i<ndata ; i++) if( x < xdata[i] ) break;
// get g(x) by linear interpolation between xdata[i-1] and xdata[i]  
double g1 = xfdata[i-1] / xdata[i-1];
double g2 = xfdata[i] / xdata[i];
double dx= xdata[i] - xdata[i-1];
return g1 * (xdata[i]-x)/dx + g2 * (x-xdata[i-1])/dx;
}

//-------------------------------------------
//   calculate QCD coupling constant alphaS()
//--------------------------------------------
double alphaS(double Q) {
return 12.*3.141593/(33-2*5)/log(Q*Q/0.1/0.1);
}

//-----------------------------------------
//   form factor |I(x)|^2
//-----------------------------------------
double I2(double x){
double pi = 3.141593;
double lambda= x*x;
double freal = -2.*pow(asin(0.5/sqrt(lambda)),2);
double fimag = 0.;
if(lambda < 0.25) {
   double etap = 0.5 + sqrt(0.25-lambda);
   double etam = 0.5 - sqrt(0.25-lambda);
   freal = 0.5*pow(log(etap/etam),2) - pi*pi/2;
   fimag = pi * log(etap/etam);
}
double Ireal = 3.*(2.*lambda + lambda*(4.*lambda-1) * freal);
double Iimag = 3.*lambda*(4.*lambda-1) * fimag;
return Ireal * Ireal + Iimag * Iimag;
}