/*
  test.cc
  Stefan Gieseke 2005
  simple testprogram to test functionality of LHAPDFWrap class, 
  a C++ wrapper for the LHAPDF package. 
  use as a template for your program. 

  pdfName is a nice tool if you want to work with several distributions
  and don't like typing the literal paths all the time. 

  the main program can also be called with two arguments 
  ./pdfset <pdfset> <subset>
  in order to check specific sets. 
  
  Adapted for lhapdfv4 by Mike Whalley
 */

#include <iostream>
#include <string>
#include "LHAPDFWrap.h"

using namespace std;

string pdfName(const int& num) {
  // EDIT THIS PATH!
  string path("../PDFsets/");
  string file;
  switch (num) {
// proton sets
  case 1: file = "a02m_lo.LHgrid"; break;
  case 2: file = "a02m_nlo.LHgrid"; break;
  case 3: file = "a02m_nnlo.LHgrid"; break;
  case 4: file = "Alekhin_100.LHpdf"; break;
  case 5: file = "Alekhin_1000.LHpdf"; break;
  case 6: file = "Botje_100.LHpdf"; break;
  case 7: file = "Botje_1000.LHpdf"; break;
  case 8: file = "cteq61.LHpdf"; break;
  case 9: file = "cteq61.LHgrid"; break;
  case 10: file = "cteq6.LHpdf"; break;
  case 11: file = "cteq6m.LHpdf"; break;
  case 12: file = "cteq6mE.LHgrid"; break;
  case 13: file = "cteq6l.LHpdf"; break;
  case 14: file = "cteq6ll.LHpdf"; break;
  case 15: file = "cteq5m.LHgrid"; break;
  case 16: file = "cteq5m1.LHgrid"; break;
  case 17: file = "cteq5d.LHgrid"; break;
  case 18: file = "cteq5l.LHgrid"; break;
  case 19: file = "cteq4m.LHgrid"; break;
  case 20: file = "cteq4d.LHgrid"; break;
  case 21: file = "cteq4l.LHgrid"; break;
  case 22: file = "Fermi2002_100.LHpdf"; break;
  case 23: file = "Fermi2002_1000.LHpdf"; break;
  case 24: file = "GRV98lo.LHgrid"; break;
  case 25: file = "GRV98nlo.LHgrid"; break;
  case 26: file = "H12000ms.LHgrid"; break;
  case 27: file = "H12000msE.LHgrid"; break;
  case 28: file = "H12000dis.LHgrid"; break;
  case 29: file = "H12000disE.LHgrid"; break;
  case 30: file = "H12000lo.LHgrid"; break;
  case 31: file = "H12000loE.LHgrid"; break;
  case 32: file = "H12000lo2.LHgrid"; break;
  case 33: file = "H12000lo2E.LHgrid"; break;
  case 34: file = "MRST2004nlo.LHpdf"; break;
  case 35: file = "MRST2004nlo.LHgrid"; break;
  case 36: file = "MRST2004nnlo.LHgrid"; break;
  case 37: file = "MRST2003cnlo.LHgrid"; break;
  case 38: file = "MRST2003cnnlo.LHgrid"; break;
  case 39: file = "MRST2002nlo.LHpdf"; break;
  case 40: file = "MRST2002nlo.LHgrid"; break;
  case 41: file = "MRST2002nnlo.LHgrid"; break;
  case 42: file = "MRST2001E.LHpdf"; break;
  case 43: file = "MRST2001E.LHgrid"; break;
  case 44: file = "MRST2001lo.LHgrid"; break;
  case 45: file = "MRST2001nlo.LHpdf"; break;
  case 46: file = "MRST2001nlo.LHgrid"; break;
  case 47: file = "MRST2001nnlo.LHgrid"; break;
  case 48: file = "MRST98.LHpdf"; break;
  case 49: file = "MRST98lo.LHgrid"; break;
  case 50: file = "MRST98nlo.LHgrid"; break;
  case 51: file = "MRST98dis.LHgrid"; break;
  case 52: file = "MRST98ht.LHgrid"; break;
  case 53: file = "ZEUS2002_TR.LHpdf"; break;
  case 54: file = "ZEUS2002_FF.LHpdf"; break;
  case 55: file = "ZEUS2002_ZM.LHpdf"; break;
  case 56: file = "ZEUS2005_ZJ.LHpdf"; break;
// pion sets
  case 57: file = "OWPI.LHgrid"; break;
  case 58: file = "SMRSPI.LHgrid"; break;
  case 59: file = "GRVPI0.LHgrid"; break;
  case 60: file = "GRVPI1.LHgrid"; break;
  case 61: file = "ABFKWPI.LHgrid"; break;
// photon sets
  case 62: file = "DOG0.LHgrid"; break;
  case 63: file = "DOG1.LHgrid"; break;
  case 64: file = "DOGG.LHgrid"; break;
  case 65: file = "LACG.LHgrid"; break;
  case 66: file = "GSG0.LHgrid"; break;
  case 67: file = "GSG1.LHgrid"; break;
  case 68: file = "GSG960.LHgrid"; break;
  case 69: file = "GSG961.LHgrid"; break;
  case 70: file = "GRVG0.LHgrid"; break;
  case 71: file = "GRVG1.LHgrid"; break;
  case 72: file = "ACFGPG.LHgrid"; break;
  case 73: file = "WHITG.LHgrid"; break;
  case 74: file = "SASG.LHgrid"; break;
//
  default: file = "MRST2004nlo.LHgrid"; break;
  }
  path+=file;
  return path;
}

int main(int argc, char **argv) {
  int pdfset = 73, subset = 0;
  if (argc == 3) { 
    pdfset = strtol(argv[1], 0, 10);
    subset = strtol(argv[2], 0, 10);
  }

//  for (int pdfset=1; pdfset<57; pdfset++){
  
  
  char name[160];
  strcpy(name, pdfName(pdfset).c_str());
  LHAPDFWrap pdf = LHAPDFWrap(name, 0);
  pdf.initPDF(subset);
  double x = 0.01, Q = 10.0, Q2 = 100.0, P2 = 1.0;
  int ip = 0;

    if(pdf.numberPDF() > 10){ cout <<"yyyy "<<pdf.numberPDF();}
    if(pdf.numberPDF() < 10){ cout <<"yyyy  "<<pdf.numberPDF();}
    cout <<" "<<pdf.alphasPDF(91.2);
    cout <<" "<<pdf.getLam4(subset);
    cout <<" "<<pdf.getLam5(subset);
    cout <<" "<<pdf.getOrderPDF();
    cout <<" "<<pdf.getOrderAlphaS();
    cout <<" "<<pdf.getNf();
    cout <<" "<<name;
    cout <<endl;

  cout << "----------------------------------------"
      << "----------------------------------------" << endl
       << " test: checking class LHAPDFWrap with package " << endl
       << name << endl
       << " set number " << subset << "/" << pdf.numberPDF() << endl;
  cout << " - pdf.getdescription() gives -----------"
      << "----------------------------------------" << endl;
  pdf.getDescription();
  cout << " ----------------------------------------"
       << "----------------------------------------" << endl;
  cout << " #pdf sets in this distribution = " << pdf.numberPDF() << endl
       << " as(91.2) = " << pdf.alphasPDF(91.2) << endl
       << " lam4 for   set = " << pdf.getLam4(subset) << endl
       << " lam5 for   set = " << pdf.getLam5(subset) << endl
       << " O(pdf) = " << pdf.getOrderPDF() << endl
       << " O(as) = " << pdf.getOrderAlphaS() << endl
       << " Nf in this set = " << pdf.getNf() << endl;
  int nf = pdf.getNf();
  cout << "----------------------------------------"
       << "---------------------------------------" << endl;
  for (int f=-nf; f<nf+1; f++) {
    cout << "f = " << f << "; mf = " << pdf.getQMass(f)
       << ", mth = " << pdf.getThreshold(f) << endl;
  }
  cout << "xxxx ----------------------------------------"
       << "----------------------------------------" << endl;
  cout << "x = " << x << ", Q = " << Q << endl;
// next for(){} is for protons pions
  for (int f=-nf; f<nf+1; f++) {
    cout << "f(x, Q, f = " << f << ") = " << pdf.xfx(x, Q, f) << endl;
  }
// next for(){} is for protons photons
  for (int f=-nf; f<nf+1; f++) {
    cout << "f(x, Q2, P2, IP, f = " << f << ") = " << pdf.xfxp(x, Q2, P2, ip, f) << endl;
  }

//  }
  
  return 0;
}