// 09/08/2012 by Graeme Watt <Graeme.Watt(at)cern.ch>.
// ROOT script to convert MSTW .LHgrid files from Hessian to random.
// Input: best-fit and 2*20 eigenvector PDF sets in one .LHgrid file.
// Output: average and 40 random PDF sets in a single new .LHgrid file.
// Pass random seed as argument.  If no argument is given, then a
// different seed is used each time, so the results are not reproducible.
// (This may be an advantage in generating many independent .LHgrid files.)

void conversion(int seed=0)
{

  gROOT->Reset();
  random = new TRandom3(seed); // random number generator

  string oldprefix = "MSTW2008nlo68cl"; // prefix of existing .LHgrid file
  string newprefix = oldprefix+"_rand"; // prefix for new random .LHgrid file

  // Options to print out Hessian and random predictions for each grid point.
  bool lPrintout=true; // print out for all 9*63*48 or 11*63*48 grid points
  //bool lPrintout=false; // suppress printing out of predictions

  // Options to symmetrise (or not) uncertainties and random predictions.
  bool lSymmetrise=true; // symmetric (then average equals best-fit)
  //bool lSymmetrise=false; // asymmetric (then average differs from best-fit)
  
  // Options to correct from 90% C.L. eigenvector PDF sets to 68% C.L.
  // It is recommended to only use the 68% C.L. sets from MSTW.
  bool l90cl=false; // 68% C.L.
  //bool l90cl=true; // 90% C.L.


  // Open input and output .LHgrid files.
  ifstream infile; ofstream outfile;
  infile.open((oldprefix+".LHgrid").c_str());
  outfile.open((newprefix+".LHgrid").c_str());
  if (!infile.good()) {
    cerr << "Error reading " << (oldprefix+".LHgrid") << endl;
    exit(-1);
  }
  if (!outfile.good()) {
    cerr << "Error reading " << (newprefix+".LHgrid") << endl;
    exit(-1);
  }

  // Copy header mostly unchanged, other than two lines of description.
  string str,line; int found;
  while (infile && line!="40,0") {
    getline(infile,line);
    str="central value";
    found=line.find(str);
    if (found!=string::npos)
      line.replace(found,str.length(),"average value over 40 random sets");
    str="20 eigenvector sets (+/- directions)";
    found=line.find(str);
    if (found!=string::npos)
      line.replace(found,str.length(),"40 random sets");
    outfile << line << endl;
    //cout << line << endl; // for testing
  }

  // Read best-fit and eigenvector PDF sets from input .LHgrid file.
  cout << "Reading grids from " << (oldprefix+".LHgrid") << " ..." << endl;
  // Define size of grid: number of members, flavours, x and q points.
  const int nmem=40, np=12, nx=64, nq=48;
  double f[np+1][nx+1][nq+1][nmem+1]; // value of PDF at grid points
  char buffer[256],ch; int alphaSorder;
  date = new TDatime(); // store the current date and time
  string **header = new string*[nmem+1]; // string array to store header
  bool lNNLO; // true for NNLO .LHgrid files (two extra columns)
  for (int imem=0;imem<=nmem;imem++) { // loop over members
    header[imem] = new string[12]; // string array to store header
    // Copy header mostly unchanged, other than three lines of description.
    for (int i=0;i<=10+(imem!=0);i++) {
      getline(infile,line);
      if (i == 0+(imem!=0)) // replace the date and time
	line = " # " + string(date->AsString());
      str="central";
      found=line.find(str);
      if (found!=string::npos)
	line.replace(found,line.length(),"average over 40 random sets");
      str="eigenvector";
      found=line.find(str);
      if (found!=string::npos) {
	sprintf(buffer,"random set number %d",imem);
	line.replace(found,line.length(),buffer);
      }
      // Distance and tolerance values are not relevant
      // for random PDF sets, so just replace by zeros.
      str="tolerance";
      found=line.find(str);
      if (found!=string::npos)
	line.replace(found,line.length(),"tolerance(T) =   0.  0.");
      header[imem][i] = line; // store to write later in new .LHgrid file
      //cout << i << " " << header[imem][i] << endl; // for testing
      if (i == 7+(imem!=0)) {
	// Use value of "alphaSorder" to determine if NNLO.
	stringstream ss; ss << line;
	ss >> ch >> str >> ch >> alphaSorder;
	if (alphaSorder==2) lNNLO = true;
	else lNNLO = false;
      }
    }
    // Read in the grid for this member set.
    for (int ix=1;ix<=nx-1;ix++) { // loop over x values
      for (int iq=1;iq<=nq;iq++) { // loop over q values
	for (int ip=1;ip<=9;ip++) { // loop over flavours
	  infile >> f[ip][ix][iq][imem]; // value of PDF
	}
	if (lNNLO) { // only at NNLO
	  infile >> f[10][ix][iq][imem]; // = chm-cbar
	  infile >> f[11][ix][iq][imem]; // = bot-bbar
	}
      } // ix
    } // iq
  } // imem
  cout << "... DONE!" << endl;

  // Now generate the random predictions.
  cout << "Calculating random predictions ..." << endl;
  const int Npdf = nmem; // = 40 (number of random PDF sets)
  const int neigen = nmem/2; // = 20 (number of eigenvectors)
  double frand[np+1][nx+1][nq+1][Npdf+1]; // value of random PDF
  double *rand = new double[neigen*Npdf+1]; // array of random numbers
  for (int ir=0;ir<=neigen*Npdf;ir++) { // 20*40+1 random numbers 
    rand[ir] = random->Gaus(0.0,1.0); // Gaussian random number
  }
  for (int ip=1;ip<=9+2*lNNLO;ip++) { // loop over flavours
    for (int ix=1;ix<=nx-1;ix++) { // loop over x values
      for (int iq=1;iq<=nq;iq++) { // loop over q values
   	// Get all Npdf random predictions for one parton flavour f(x,q).
	double *frandtemp = getRandomPredictions(neigen,f[ip][ix][iq],rand,Npdf,lSymmetrise,l90cl);
	if (lPrintout) {
	  // Calculate the Hessian PDF uncertainty (only needed to print out).
	  double *unchess = getPDFuncertainty(nmem,f[ip][ix][iq],false,l90cl);
	  cout << ip << " " << ix << " " << iq << ", Hess: " << unchess[0] << "+/-" << unchess[3] << " (+" << unchess[1] << ", -" << unchess[2] << ")" << endl;
	}
  	// Calculate the average and s.d. over the Npdf random predictions.
	double *uncrand = getPDFuncertainty(Npdf,frandtemp,true,l90cl);
  	if (lPrintout) cout << ip << " " << ix << " " << iq << ", Rand: " << uncrand[0] << "+/-" << uncrand[3] << endl;
  	// Store the average as the zeroth member set.
  	frand[ip][ix][iq][0] = uncrand[0];
	for (int ipdf=1;ipdf<=Npdf;ipdf++) {
	  frand[ip][ix][iq][ipdf] = frandtemp[ipdf];
	}
      } // iq
    } // ix
  } // ip
  cout << "... DONE!" << endl;

  // Write average and random PDF sets to new .LHgrid file.
  cout << "Writing grids to " << (newprefix+".LHgrid") << " ..." << endl;
  for (int imem=0;imem<=nmem;imem++) { // loop over members
    for (int i=0+(imem!=0);i<=10+(imem!=0);i++) {
      // Copy header previously stored in a string array.
      outfile << header[imem][i] << endl;
    }
    // Write the grid for this member set.
    for (int ix=1;ix<=nx-1;ix++) { // loop over x values
      for (int iq=1;iq<=nq;iq++) { // loop over q values
  	for (int ip=1;ip<=9;ip++) { // loop over flavours
  	  sprintf(buffer,"%12.4E",frand[ip][ix][iq][imem]);
  	  outfile << buffer;
  	  //cout << buffer; // for testing
  	}
  	if (lNNLO) { // only at NNLO (chm-cbar and bot-bbar)
	  sprintf(buffer,"%12.4E%12.4E",frand[10][ix][iq][imem],frand[11][ix][iq][imem]);
	  outfile << buffer;
  	  //cout << buffer; // for testing
  	}
  	outfile << endl;
  	//cout << endl; // for testing
      } // ix
    } // iq
  } // imem

  // Copy final line of input .LHgrid file.
  getline(infile,line); getline(infile,line); outfile << line << endl;
  //cout << line << endl; // for testing (should be 'End:')

  // Check that end of input .LHgrid file has been reached.
  getline(infile,line);
  if (!infile.eof()) {
    cerr << "Error reading " << (oldprefix+".LHgrid") << endl;
    exit(-1);
  }
  cout << "... DONE!" << endl;

  // Close input and output .LHgrid files.
  infile.close();
  outfile.close();

}
// End of main function "void conversion(int seed=0)".


// Function to calculate either Monte Carlo or Hessian PDF uncertainties.
// Pass an array f[imem] (imem=0,...,Nmem).  l90cl to correct for 90% C.L.
double *getPDFuncertainty(int nmem, double *f, bool lMonteCarlo, bool l90cl)
{
  double *ret = new double[4];
  double central=0.0,errsym=0.0,errplus=0.0,errminus=0.0;
  if (lMonteCarlo) {
    // Calculate the average and standard deviation
    // using Eqs. (2.13) and (2.14) of arXiv:1205.4024v2.
    double fav = 0.0, fsd = 0.0;
    for (int imem=1;imem<=nmem;imem++) {
      fav += f[imem];
      fsd += f[imem]**2;
    }
    fav /= nmem; fsd /= nmem;
    fsd = nmem/(nmem-1.0)*(fsd-fav**2);
    if (fsd>0.0) fsd = sqrt(fsd);
    else fsd = 0.0;    
    central = fav;
    errplus = errminus = errsym = fsd;
  }
  else {
    // Calculate the symmetric and asymmetric Hessian uncertainties
    // using Eqs. (2.7), (2.8) and (2.9) of arXiv:1205.4024v2.
    central = f[0];
    for (int ieigen=1;ieigen<=nmem/2;ieigen++) {
      errplus += (TMath::Max(TMath::Max(f[2*ieigen-1]-f[0],f[2*ieigen]-f[0]),0.0))**2;
      errminus += (TMath::Max(TMath::Max(f[0]-f[2*ieigen-1],f[0]-f[2*ieigen]),0.0))**2;
      errsym += (f[2*ieigen-1]-f[2*ieigen])**2;
    }
    errsym = 0.5*sqrt(errsym);
    errplus = sqrt(errplus);
    errminus = sqrt(errminus);
    if (l90cl) { // correct from 90% C.L. to 68% C.L. (for CT10)
      errsym /= 1.64485; errplus /= 1.64485; errminus /= 1.64485;
    }
  }
  // Return central value and asymmetric/symmetric PDF uncertainties.
  ret[0] = central;
  ret[1] = errplus;
  ret[2] = errminus;
  ret[3] = errsym;
  return ret;
}


// Function to generate Npdf random predictions using
// either Eq. (6.4) or Eq. (6.5) of arXiv:1205.4024v2.
// Pass an array f[imem] (imem=0,...,2*neigen).
double *getRandomPredictions(int neigen, double *f, double *rand, int Npdf, bool lSymmetrise, bool l90cl)
{
  double *frand = new double[Npdf+1];
  double scale=1.0;
  if (l90cl) scale=1.64485; // correct from 90% C.L. to 68% C.L.
  for (int ipdf=1;ipdf<=Npdf;ipdf++) {
    frand[ipdf] = f[0];
    // Loop over number of eigenvectors.
    for (int ieigen=1;ieigen<=neigen;ieigen++) {
      double r = rand[(ipdf-1)*neigen+ieigen]; // Gaussian random number
      if (lSymmetrise) frand[ipdf] +=
			 0.5*r*TMath::Abs(f[2*ieigen-1]-f[2*ieigen])/scale;
      else { // not symmetrised
	if (r<0.0) frand[ipdf] -=
		     r*(f[2*ieigen]-f[0])/scale; // negative direction
	else frand[ipdf] +=
	       r*(f[2*ieigen-1]-f[0])/scale; // positive direction
      }
    }
  }
  return frand;
}