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

int main()
{
readinPDF();
//
for(int i=0; i<30; i++) {
  double tau = 0.00001*i ;
  std::cout<<"tau="<<tau<<"  L(tau)="<<partonLuminosity(tau)<<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;
ndata = 0;
std::ifstream fin("gluon_10000.txt");
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 <<std::endl;
return;
}
 
//--------------------------------------------
//   calculate gluon density g(x)
//--------------------------------------------
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;
}