#ifndef _LHAPDFWRAP_H
#define _LHAPDFWRAP_H

#include <vector>

/*
  This class is a wrapper around the LHAPDFv2/3 package for parton
  distribution functions of the proton.  
  Adapted for LHAPDFv4 by Mike Whalley
  Adapted for LHAPDFv5 by Craig Group/Mike Whalley
 */

extern "C" {
  void finitpdfset_(char&);
  void finitpdfsetbyname_(char&);
  void finitpdf_(int&);
  void fevolvepdf_(double&, double &, double*);
  void fevolvepdfp_(double&, double &, double&, int&, double*);
  void fevolvepdfa_(double&, double &, double &, double*);
  void fevolvepdfphoton_(double&, double &, double*, double&);
  void fnumberpdf_(int&);
  void falphaspdf_(double&, double &);
  void fgetorderpdf_(int&);
  void fgetorderas_(int&);
  void fgetdesc_();
  void fgetqmass_(int&, double&);
  void fgetthreshold_(int&, double&);
  void fgetnf_(int&);
  void fgetlam4_(int&, double&);
  void fgetlam5_(int&, double&);
  void fgetxmin_(int&, double&);
  void fgetxmax_(int&, double&);
  void fgetq2min_(int&, double&);
  void fgetq2max_(int&, double&);
  void fgetminmax_(int&, double&, double&, double&, double&);
  void fextrapolate_();

  //v5 subroutines for multiple set initialization
  void finitpdfsetm_(int&, char&);
  void finitpdfsetbynamem_(int&, char&);
  void finitpdfm_(int&, int&);
  void fevolvepdfm_(int&, double&, double &, double*);
  void fevolvepdfpm_(int&, double&, double &, double&, int&, double*);
  void fevolvepdfam_(int&, double&, double &, double &, double*);
  void fevolvepdfphotonm_(int&, double&, double &, double*, double&);
  void fnumberpdfm_(int&, int&);
  void falphaspdfm_(int&, double&, double &);
  void fgetorderpdfm_(int&, int&);
  void fgetorderasm_(int&, int&);
  void fgetdescm_(int&);
  void fgetqmassm_(int&, int&, double&);
  void fgetthresholdm_(int&, int&, double&);
  void fgetnfm_(int&, int&);
  void fgetlam4m_(int&, int&, double&);
  void fgetlam5m_(int&, int&, double&);
  void fgetxminm_(int&, int&, double&);
  void fgetxmaxm_(int&, int&, double&);
  void fgetq2minm_(int&, int&, double&);
  void fgetq2maxm_(int&, int&, double&);
  void fgetminmaxm_(int&, int&, double&, double&, double&, double&);

}

class LHAPDFWrap {

public:
  LHAPDFWrap();
  ~LHAPDFWrap();
  // std c'tor d'tor
  LHAPDFWrap(char *name);
  // typical constructor with pdfset 'name' 
  // 'name' is the name of the file
  // the grid or data file of the desired set.
 
  LHAPDFWrap(char *name, int memset);
  // typical constructor with pdfset 'name' and subset 'memset' 
  //   'name' is the name 
  // the grid or data file of the desired set.


  LHAPDFWrap(int nset, char *name);
  // typical constructor (when multiple PDF sets need to be initialized) 
  // with pdfset 'name' and subset 'memset' 
  //  'name' is name of
  // the grid or data file of the desired set.  
  // int nset specifies the reference number for the set to be initialized

  LHAPDFWrap(int nset, char *name, int memset);
  // typical constructor (when multiple PDF sets need to be initialized) 
  // with pdfset 'name' and subset 'memset' 
  //'name' is the name
  // the grid or data file of the desired set.  
  // int nset specifies the reference number for the set to be initialized


  std::vector< double > xfx(const double &x, const double &Q);
  // returns a vector xf(x, Q) with index 0 < i < 12.
  // 0..5 = tbar, ..., ubar, dbar; 
  // 6 = g; 
  // 7..12 = d, u, ..., t
  double xfx(const double &x, const double &Q, int fl);
  // returns xf(x, Q) for flavour fl - this time the flavour encoding
  // is as in the LHAPDF manual...
  // -6..-1 = tbar,...,ubar, dbar
  // 1..6 = duscbt
  // 0 = g
  std::vector< double > xfxp(const double &x, const double &Q, const double &P2, int ip);
  double xfxp(const double &x, const double &Q, const double &P2, int ip, int fl);

  std::vector< double > xfxa(const double &x, const double &Q, const double &a);
  double xfxa(const double &x, const double &Q, const double &a, int fl);

  std::vector< double > xfxphoton(const double &x, const double &Q);
  double xfxphoton(const double &x, const double &Q, int fl);

  inline void initPDFSet(char *name) {finitpdfset_(*name);}
  // the pdfset by name, see subdir 'PDFset' of LHAPDFv2 for choices
  inline void initPDFSetByName(char *name) {finitpdfsetbyname_(*name);}
  // the pdfset by name, see subdir 'PDFset' of LHAPDFv2 for choices
  inline void initPDF(int memset) {finitpdf_(memset);}
  // the choice of pdf subset out of one distribution

  inline void getDescription() {fgetdesc_();}
  // prints a brief description of the current pdf set to stdout
  int numberPDF();
  // number of subsets available in the current distribution.
  double alphasPDF(double Q);
  // alphas used by the current pdf.
  int getOrderPDF();
  int getOrderAlphaS();
  // perturbative order of parton evolution and alpha_S resp.
  double getQMass(int f);
  // quark mass used for flavour f.
  double getThreshold(int f);
  // threshold for flavour f.
  int getNf();
  // number of flavours used in the current pdf set.
  double getLam4(int m);
  // value of qcd lambda4 for member m 
  double getLam5(int m);
  // value of qcd lambda5 for member m 
  double getXmin(int m);
  double getXmax(int m);
  double getQ2min(int m);
  double getQ2max(int m);

  inline void extrapolate(){ fextrapolate_();};


  //Additional functions for when more than 1 PDF set is being stored in memory


  std::vector< double > xfxM(int nset, const double &x, const double &Q);
  // returns a vector xf(x, Q) with index 0 < i < 12.
  // 0..5 = tbar, ..., ubar, dbar; 
  // 6 = g; 
  // 7..12 = d, u, ..., t
  double xfxM(int nset, const double &x, const double &Q, int fl);
  // returns xf(x, Q) for flavour fl - this time the flavour encoding
  // is as in the LHAPDF manual...
  // -6..-1 = tbar,...,ubar, dbar
  // 1..6 = duscbt
  // 0 = g
  std::vector< double > xfxpM(int nset, const double &x, const double &Q, const double &P2, int ip);
  double xfxpM(int nset, const double &x, const double &Q, const double &P2, int ip, int fl);

  std::vector< double > xfxaM(int nset, const double &x, const double &Q, const double &a);
  double xfxaM(int nset, const double &x, const double &Q, const double &a, int fl);

  std::vector< double > xfxphotonM(int nset, const double &x, const double &Q);
  double xfxphotonM(int nset, const double &x, const double &Q, int fl);

  inline void initPDFSetM(int nset, char *name) {finitpdfsetm_(nset, *name);}
  // the pdfset by name, see subdir 'PDFset' of LHAPDFv2 for choices
  inline void initPDFSetByNameM(int nset, char *name) {finitpdfsetbynamem_(nset, *name);}
  // the pdfset by name, see subdir 'PDFset' of LHAPDFv2 for choices
  inline void initPDFM(int nset, int memset) {finitpdfm_(nset, memset);}
  // the choice of pdf subset out of one distribution


  inline void getDescriptionM(int nset) {fgetdescm_(nset);}
  // prints a brief description of the current pdf set to stdout
  int numberPDFM(int nset);
  // number of subsets available in the current distribution.
  double alphasPDFM(int nset, double Q);
  // alphas used by the current pdf.
  int getOrderPDFM(int nset);
  int getOrderAlphaSM(int nset);
  // perturbative order of parton evolution and alpha_S resp.
  double getQMassM(int nset, int f);
  // quark mass used for flavour f.
  double getThresholdM(int nset, int f);
  // threshold for flavour f.
  int getNfM(int nset);
  // number of flavours used in the current pdf set.
  double getLam4M(int nset, int m);
  // value of qcd lambda4 for member m 
  double getLam5M(int nset, int m);
  // value of qcd lambda5 for member m 

  double getXminM(int nset, int m);
  double getXmaxM(int nset, int m);
  double getQ2minM(int nset, int m);
  double getQ2maxM(int nset, int m);

};

#endif