Code examples

The following examples should help you to get to grips with using LHAPDF from C++:

Using and testing PDFs in C++

// Program to test LHAPDF6 PDF behaviour by writing out their values at lots of x and Q points
// Note: the OpenMP directives are there as an example. In fact, in this case OpenMP slows things
// down because of the need to make the stream operations critical!
 
#include "LHAPDF/LHAPDF.h"
#include <iostream>
#include <fstream>
using namespace LHAPDF;
using namespace std;
 
int main(int argc, char* argv[]) {
 
  if (argc < 3) {
    cerr << "You must specify a PDF set and member number" << endl;
    return 1;
  }
 
  const string setname = argv[1];
  const string smem = argv[2];
  const int imem = lexical_cast<int>(smem);
 
  const PDF* pdf = mkPDF(setname, imem);
  vector<int> pids = pdf->flavors();
 
  const double MINLOGX = -10;
  const double MAXLOGX =   0;
  const double DX = 0.01;
  const int NX = (int) floor((MAXLOGX - MINLOGX)/DX) + 1;
 
  const double MINLOGQ2 = 1;
  const double MAXLOGQ2 = 8;
  const double DQ2 = 0.01;
  const int NQ2 = (int) floor((MAXLOGQ2 - MINLOGQ2)/DQ2) + 1;
 
  for (int pid : pids) {
    const string spid = lexical_cast<string>(pid);
    const string filename = setname + "_" + smem + "_" + spid + ".dat";
    ofstream f(filename.c_str());
    for (int ix = 0; ix < NX; ++ix) {
      const double log10x = (MINLOGX + ix*DX < -1e-3) ? MINLOGX + ix*DX : 0;
      const double x = pow(10, log10x);
      for (int iq2 = 0; iq2 < NQ2; ++iq2) {
        const double log10q2 = MINLOGQ2 + iq2*DQ2;
        const double q2 = pow(10, log10q2);
        const double xf = pdf->xfxQ2(pid, x, q2);
        f << x << " " << q2 << " " << xf << endl;
      }
    }
    f.close();
  }
 
  for (double log10q2 = MINLOGQ2; log10q2 <= MAXLOGQ2; log10q2 += 0.2) {
    const double q2 = pow(10, log10q2);
    cout << "alpha_s(" << setprecision(1) << fixed << sqrt(q2) << " GeV) = "
         << setprecision(5) << pdf->alphasQ2(q2) << endl;
  }
 
  delete pdf;
  return 0;
}

// Program to test LHAPDF6 PDF CPU/memory performance by sampling PDFs in several ways
 
#include "LHAPDF/LHAPDF.h"
#include <iostream>
#include <cmath>
#include <ctime>
#ifdef HAVE_MPI
#include <mpi.h>
#endif
using namespace std;
 
int main(int argc, char* argv[]) {
 
#ifdef HAVE_MPI
  MPI_Init(&argc, &argv);
#endif
  const double MINLOGX = -7.5;
  const double MINLOGQ = 1;
  const double MAXLOGQ = 3;
  const double dx = 0.005;
  const double dq = 0.005;
 
  const clock_t start = clock();
 
  #if LHAPDF_MAJOR_VERSION > 5
  const LHAPDF::PDF* pdf = LHAPDF::mkPDF("CT10nlo", 0);
  #define XFS(X, Q) pdf->xfxQ(X, Q, xfs)
  #else
  LHAPDF::initPDFSetByName("CT10nlo.LHgrid");
  LHAPDF::initPDF(0);
  #define XFS(X, Q) LHAPDF::xfx(X, Q, &xfs[0])
  #endif
 
  const clock_t init = clock();
 
  vector<double> xfs; xfs.resize(13);
  for (double log10x = MINLOGX; log10x <= 0.0; log10x += dx) {
    for (double log10q = MINLOGQ; log10q <= MAXLOGQ; log10q += dq) {
      XFS(pow(10, log10x), pow(10, log10q));
    }
  }
 
  const clock_t end = clock();
 
  std::cout << "Init  = " << (init - start) << std::endl;
  std::cout << "Query = " << (end - init) << std::endl;
  std::cout << "Total = " << (end - start) << std::endl;
 
  #if LHAPDF_MAJOR_VERSION > 5
  delete pdf;
  #endif
 
#ifdef HAVE_MPI
  MPI_Finalize();
#endif
  return 0;
}

Using PDFs in Fortran

program example1
  implicit double precision (a-h,o-z)
  character name*64
  double precision f(-6:6)
  character*20 lparm
  logical has_photon
  dimension z(10), xx(10)
  Data (z(i), i=1,10) /.05, .1, .2, .3, .4, .5, .6, .7, .8, .9/
 
  Do i = 1, 10
     xx(i) = z(i) **3
  EndDo
 
  name='CT10nlo.LHgrid'
  ! name='cteq6.LHgrid'
  ! name='cteq65.LHgrid'
  ! name='cteq66.LHgrid'
  ! name='MRST2006nnlo.LHgrid'
  ! name='MRST2001E.LHgrid'
  ! name='H12000ms.LHgrid'
  call initpdfsetbyname(name)
 
  qmz=91.18d0
  write(*,*)
  call numberpdf(n)
  print *,'There are ',n,' PDF sets'
  do i=0,n
     write(*,*) '---------------------------------------------'
     call initpdf(i)
     write(*,*) 'PDF set ',i
     call getxmin(i,xmin)
     call getxmax(i,xmax)
     call getq2min(i,q2min)
     call getq2max(i,q2max)
     print *,'xmin=',xmin,' xmax=',xmax,' Q2min=',q2min,' Q2max=',q2max
     call getminmax(i,xmin,xmax,q2min,q2max)
     print *,'xmin=',xmin,' xmax=',xmax,' Q2min=',q2min,' Q2max=',q2max
     call setlhaparm('EXTRAPOLATE') 
     call getlhaparm(18,lparm) 
     print *,'lhaparm(18)=',lparm
     write(*,*)
     a=alphaspdf(qmz)
     write(*,*) 'alpha_S(M_Z) = ',a
     call getlam4m(1,i,xlam4)
     call getlam5m(1,i,xlam5)
     print *,' lambda5: ',xlam5, ' lambda4: ',xlam4
     write(*,*)
     write(*,*) 'x*up'
     write(*,*) '   x     Q=10 GeV     Q=100 GeV    Q=1000 GeV'
     ! q2 = 10.0d0
     ! q = dsqrt(q2)
     q = 50.0d0
     print *,q
     do ix=1,10
        ! x = (ix-0.5d0)/10.0d0
        x = xx(ix)
        ! x = z(ix)
        if(has_photon()) then
          print *,"This set has a photon"
          call evolvepdfphoton(x,q,f,photon)
        else
          call evolvepdf(x,q,f)
        endif
        g = f(0)
        u = f(2)
        d = f(1)
        s = f(3)
        c = f(4)
        b = f(5)
        ubar = f(-2)
        dbar = f(-1)
        sbar = f(-3)
        cbar = f(-4)
        bbar = f(-5)
        write(*,'(F7.4,13(1pE10.3))') x,u,d,ubar,dbar,s,sbar,c,cbar,b,bbar,g,photon
     enddo
  enddo
 
end program example1

Using and testing PDFs in Python

#! /usr/bin/env python
 
 
 
import lhapdf
 
 
p = lhapdf.mkPDF("CT10nlo", 0)
p = lhapdf.mkPDF("CT10nlo/0")
 
 
print(p.xfxQ2(21, 1e-3, 1e4))
 
 
for pid in p.flavors():
    print(p.xfxQ(pid, 0.01, 91.2))
 
# TODO: demonstrate looping over PDF set members
pset = lhapdf.getPDFSet("CT10nlo")
print(pset.description)
pcentral = pset.mkPDF(0)
pdfs1 = pset.mkPDFs()
pdfs2 = lhapdf.mkPDFs("CT10nlo") # a direct way to get all the set's PDFs
 
 
import numpy as np
xs = [x for x in np.logspace(-7, 0, 5)]
qs = [q for q in np.logspace(1, 4, 4)]
gluon_xfs = np.empty([len(xs), len(qs)])
for ix, x in enumerate(xs):
    for iq, q in enumerate(qs):
        gluon_xfs[ix,iq] = p.xfxQ(21, x, q)
print(gluon_xfs)
 
 
print(lhapdf.version())
print(lhapdf.__version__)
lhapdf.pathsPrepend("/path/to/extra/pdfsets")
print(lhapdf.paths())
# ...

Handling LHAPDF v5/v6 compatibility in C++

// -*- C++ -*-
// LHAPDFv5/v6 compatibility example
 
#include "LHAPDF/LHAPDF.h"
#include <iostream>
 
 
int main() {
  const double x = 1e-3, Q = 200;
 
  #if LHAPDF_MAJOR_VERSION == 6
  LHAPDF::PDF* pdf = LHAPDF::mkPDF("CT10nlo", 0);
  std::cout << "xf_g = " << pdf->xfxQ(21, x, Q) << std::endl;
  delete pdf;
  #else
  LHAPDF::initPDFSet("CT10nlo", LHAPDF::LHGRID, 0);
  std::cout << "xf_g = " << LHAPDF::xfx(x, Q, 0) << std::endl;
  #endif
 
  return 0;
}

Using and testing PDF grids in C++

// Example program to test PDF grid format reading and interpolation
 
#include "LHAPDF/GridPDF.h"
#include <iostream>
#include <fstream>
#ifdef HAVE_MPI
#include <mpi.h>
#endif
using namespace std;
 
 
void safeprint(const LHAPDF::PDF& pdf, const string& key) {
  if (pdf.info().has_key(key))
    cout << key << " = " << pdf.info().get_entry(key) << endl;
}
 
 
int main(int argc, char* argv[]) {
#ifdef HAVE_MPI
  MPI_Init(&argc, &argv);
#endif
 
  if (argc < 2) {
    cout << "Usage: testgrid <PDFNAME=CT10nlo>" << endl;
    //exit(1);
  }
  const string setname = (argc < 2) ? "CT10nlo" : argv[1];
 
  const LHAPDF::PDF* basepdf = LHAPDF::mkPDF(setname);
  const LHAPDF::GridPDF& pdf = * dynamic_cast<const LHAPDF::GridPDF*>(basepdf);
 
  for (const string& p : LHAPDF::paths()) cout << p << " : ";
  cout << endl;
 
  safeprint(pdf, "Verbosity");
  safeprint(pdf, "PdfDesc");
  safeprint(pdf, "SetDesc");
  cout << "Flavors (str) = " << pdf.info().get_entry("Flavors") << endl;
  vector<int> pids = pdf.info().get_entry_as< vector<int> >("Flavors");
  cout << "Flavors (ints) = ";
  for (int f : pids) cout << f << " ";
  cout << endl;
  cout << "Flavors (vec<int>) = " << LHAPDF::to_str(pids) << endl;
 
  cout << "x0, Q0 = " << pdf.knotarray().xf(21, 0, 0) << endl;
  cout << "x1, Q0 = " << pdf.knotarray().xf(21, 1, 0) << endl;
  cout << "x0, Q1 = " << pdf.knotarray().xf(21, 0, 1) << endl;
  cout << "x1, Q1 = " << pdf.knotarray().xf(21, 1, 1) << endl;
 
  cout << pdf.xfxQ(21, 0.7, 10.0) << endl;
  cout << pdf.xfxQ(21, 0.2, 126) << endl;
  for (int pid : pdf.flavors()) {
    cout << pid << " = " << pdf.xfxQ(pid, 0.2, 124) << endl;
  }
 
  ofstream f("pdf.dat");
  for (double x = 0; x <= 1; x += 0.02) {
    for (double log10q2 = 1; log10q2 < 5; log10q2 += 0.05) {
      f << x << " " << log10q2 << " " << pdf.xfxQ2(21, x, pow(10, log10q2)) << endl;
    }
  }
  f.close();
  cout << endl;
#ifdef HAVE_MPI
  MPI_Finalize();
#endif
 
  return 0;
}

Testing the Info system in C++

// Example program for testing the info system
 
#include "LHAPDF/Info.h"
#include "LHAPDF/Config.h"
#include "LHAPDF/PDFInfo.h"
#include "LHAPDF/PDFSet.h"
#include "LHAPDF/Factories.h"
#include <iostream>
#ifdef HAVE_MPI
#include <mpi.h>
#endif
using namespace std;
 
int main(int argc, char* argv[]) {
#ifdef HAVE_MPI
  MPI_Init(&argc, &argv);
#endif
 
  LHAPDF::Info& cfg = LHAPDF::getConfig();
  // cout << "UndefFlavorAction: " << cfg.get_entry("UndefFlavorAction") << endl;
  cout << "Verbosity: " << cfg.get_entry("Verbosity") << endl;
  cfg.set_entry("Verbosity", 5);
  const LHAPDF::Info& cfg2 = LHAPDF::getConfig();
  cout << "New Verbosity from second Config: " << cfg2.get_entry("Verbosity") << endl;
 
  const LHAPDF::PDFSet set("CT10nlo");
  cout << "SetDesc: " << set.get_entry("SetDesc") << endl;
  cout << "Verbosity from set: " << set.get_entry("Verbosity") << endl;
 
  const LHAPDF::PDFInfo info("CT10nlo", 2);
  if (info.has_key("PdfDesc")) cout << "PdfDesc: " << info.get_entry("PdfDesc") << endl;
  cout << "PdfType: " << info.get_entry("PdfType") << endl;
  cout << "Verbosity from PDF: " << info.get_entry("Verbosity") << endl;
  vector<int> pids = info.get_entry_as< vector<int> >("Flavors");
  cout << "PIDs (1): "; for (int f : pids) { cout << f << " "; } cout << endl;
  cout << "PIDs (2): " << LHAPDF::to_str(pids) << endl;
 
  // Now test loading of all central PDFs
  for (const string& name : LHAPDF::availablePDFSets()) {
    cout << "Testing PDFInfo for " << name << endl;
    LHAPDF::PDFInfo* i = LHAPDF::mkPDFInfo(name, 0);
    i->has_key("Foo"); // < Force loading of all info levels
    delete i;
  }
#ifdef HAVE_MPI
  MPI_Finalize();
#endif
 
  return 0;
}

Testing the LHAPDF ID indexing system in C++

// Example program to test LHAPDF index lookup and loading
 
#include "LHAPDF/PDFIndex.h"
#include <iostream>
#ifdef HAVE_MPI
#include <mpi.h>
#endif
using namespace LHAPDF;
using namespace std;
 
void lookup(int id) {
  pair<string, int> set_id = lookupPDF(id);
  cout << "ID=" << id << " -> set=" << set_id.first << ", mem=" << set_id.second << endl;
}
 
 
int main(int argc, char* argv[]) {
#ifdef HAVE_MPI
  MPI_Init(&argc, &argv);
#endif
  lookup(10800);
  lookup(10801);
  lookup(10042);
  lookup(10041);
  lookup(10799);
  lookup(12346);
#ifdef HAVE_MPI
  MPI_Finalize();
#endif
  return 0;
}

Testing the path searching system in C++

// Test program for path searching machinery
 
#include "LHAPDF/Paths.h"
#include <iostream>
using namespace std;
 
#ifdef HAVE_MPI
#include <mpi.h>
#endif
 
int main(int argc, char* argv[]) {
 
#ifdef HAVE_MPI
  MPI_Init(&argc, &argv);
#endif
 
  for (const string& p : LHAPDF::paths())
    cout << p << endl;
 
  cout << "@" << LHAPDF::findFile("lhapdf.conf") << "@" << endl;
 
  cout << "List of available PDFs:" << endl;
  for (const string& s : LHAPDF::availablePDFSets())
    cout << " " << s << endl;
 
#ifdef HAVE_MPI
  MPI_Finalize();
#endif
  return 0;
}

Testing the AlphaS classes in C++

 // Example program to test alpha_s calculation
 
#include "LHAPDF/LHAPDF.h"
#include <iostream>
#include <fstream>
#include <limits>
 
#ifdef HAVE_MPI
#include <mpi.h>
#endif
using namespace LHAPDF;
using namespace std;
 
int main(int argc, char* argv[]) {
#ifdef HAVE_MPI
  MPI_Init(&argc, &argv);
#endif
 
  // Set up three standalone AlphaS solvers:
 
 
  AlphaS_Analytic as_ana;
 
  // Can set order of QCD (up to 4)
  // 0 returns a constant value -- needs to be set by as_ana.setAlphaSMZ(double value);
  as_ana.setOrderQCD(4);
 
  // Set quark masses for evolution (both transition points and gradients by default)
  as_ana.setQuarkMass(1, 0.0017);
  as_ana.setQuarkMass(2, 0.0041);
  as_ana.setQuarkMass(3, 0.1);
  as_ana.setQuarkMass(4, 1.29);
  as_ana.setQuarkMass(5, 4.1);
  as_ana.setQuarkMass(6, 172.5);
  as_ana.setLambda(3, 0.339);
  as_ana.setLambda(4, 0.296);
  as_ana.setLambda(5, 0.213);
 
  // Can override quark masses for Nf transitions thresholds by setting them explicitly
  // You can't mix the two: if you set one flavor threshold explicitly you need to set all of them.
  //as_ode.setQuarkThreshold(6, 650);
  //as_ode.setQuarkThreshold(5, 10);
  //as_ode.setQuarkThreshold(4, 2);
  //as_ode.setQuarkThreshold(3, 0.3);
  //as_ode.setQuarkThreshold(2, 0.1);
  //as_ode.setQuarkThreshold(1, 0.08);
 
  // Uncomment to use fixed flavor scheme for analytic solver
  // as_ana.setFlavorScheme(AlphaS::FIXED, 5);
 
 
 
  AlphaS_ODE as_ode;
 
  // As above: order = 0 returns
  // constant value set by
  // as_ode.setAlphaSMZ(double value);
  as_ode.setOrderQCD(5);
 
  as_ode.setMZ(91);
  as_ode.setAlphaSMZ(0.118);
  // as_ode.setMassReference(4.1);
  // as_ode.setAlphaSReference(0.21);
  as_ode.setQuarkMass(1, 0.0017);
  as_ode.setQuarkMass(2, 0.0041);
  as_ode.setQuarkMass(3, 0.1);
  as_ode.setQuarkMass(4, 1.29);
  as_ode.setQuarkMass(5, 4.1);
  as_ode.setQuarkMass(6, 172.5);
 
  // Can override quark masses for Nf transitions thresholds by setting them explicitly
  // You can't mix the two: if you set one flavor threshold explicitly you need to set all of them.
  //  as_ode.setQuarkThreshold(6, 650);
  //  as_ode.setQuarkThreshold(5, 10);
  //  as_ode.setQuarkThreshold(4, 2);
  //  as_ode.setQuarkThreshold(3, 0.3);
  //  as_ode.setQuarkThreshold(2, 0.1);
  //  as_ode.setQuarkThreshold(1, 0.08);
 
  // Uncomment to use fixed flavor scheme for ODE solver
  // as_ode.setFlavorScheme(AlphaS::FIXED, 4);
 
 
 
  AlphaS_Ipol as_ipol;
  vector<double> qs = { 1.300000e+00, 1.300000e+00, 1.560453e+00, 1.873087e+00, 2.248357e+00, 2.698811e+00, 3.239513e+00, 3.888544e+00, 4.667607e+00, 5.602754e+00, 6.725257e+00, 8.072650e+00, 9.689992e+00, 1.163137e+01, 1.396169e+01, 1.675889e+01, 2.011651e+01, 2.414681e+01, 2.898459e+01, 3.479160e+01, 4.176203e+01, 5.012899e+01, 6.017224e+01, 7.222765e+01, 8.669834e+01, 1.040682e+02, 1.249181e+02, 1.499452e+02, 1.799865e+02, 2.160465e+02, 2.593310e+02, 3.112875e+02, 3.736534e+02, 4.485143e+02, 5.383733e+02, 6.462355e+02, 7.757077e+02, 9.311194e+02, 1.117668e+03, 1.341590e+03, 1.610376e+03, 1.933012e+03, 2.320287e+03, 2.785153e+03, 3.343154e+03, 4.012949e+03, 4.816936e+03, 4.816936e+03 };
  vector<double> alphas = { 4.189466e-01, 4.189466e-01, 3.803532e-01, 3.482705e-01, 3.211791e-01, 2.979983e-01, 2.779383e-01, 2.604087e-01, 2.451922e-01, 2.324903e-01, 2.210397e-01, 2.106640e-01, 2.012187e-01, 1.925841e-01, 1.846600e-01, 1.773623e-01, 1.706194e-01, 1.643704e-01, 1.585630e-01, 1.531520e-01, 1.480981e-01, 1.433671e-01, 1.389290e-01, 1.347574e-01, 1.308290e-01, 1.271232e-01, 1.236215e-01, 1.203076e-01, 1.171667e-01, 1.141857e-01, 1.113525e-01, 1.086566e-01, 1.060881e-01, 1.036382e-01, 1.012990e-01, 9.906295e-02, 9.692353e-02, 9.487457e-02, 9.291044e-02, 9.102599e-02, 8.921646e-02, 8.747747e-02, 8.580498e-02, 8.419524e-02, 8.264479e-02, 8.115040e-02, 7.970910e-02, 7.970910e-02 };
  as_ipol.setQValues(qs);
  as_ipol.setAlphaSValues(alphas);
 
  // Can interpolate ODE with given knots in Q
  // as_ode.setQValues(qs);
 
 
  // Test these solvers and the CT10nlo PDF's default behaviours:
 
 
  PDF* pdf = mkPDF("CT10", 0);
  const double inf = numeric_limits<double>::infinity();
  ofstream fa("alphas_ana.dat"), fo("alphas_ode.dat"), fi("alphas_ipol.dat"), fc("alphas_ct10nlo.dat");
  cout << endl;
  for (double log10q = -0.5; log10q < 3; log10q += 0.05) {
    const double q = pow(10, log10q);
    const double as_ana_q = as_ana.alphasQ(q);
    cout << fixed;
    cout << "Q = " << setprecision(3) << q << " GeV" << endl;
    cout << "Analytical solution:           " << setprecision(3) << setw(6)  << ( (as_ana_q > 2) ? inf : as_ana_q )
         << "    num flavs = " << as_ana.numFlavorsQ(q) << endl;
    fa << q << " " << as_ana_q << endl;
 
    const double as_ode_q = as_ode.alphasQ(q);
    cout << "ODE solution:                  " << setprecision(3) << setw(6)  << ( (as_ode_q > 2) ? inf : as_ode_q )
         << "    num flavs = " << as_ode.numFlavorsQ(q) << endl;
    fo << q << " " << as_ode_q << endl;
    // const double as_ipol_q = as_ipol.alphasQ(q);
    // cout << "Interpolated solution:         " << setprecision(3) << setw(6)  << ( (as_ipol_q > 2) ? inf : as_ipol_q ) << endl;
    // fi << q << " " << as_ipol_q << endl;
    // const double as_ct10_q = pdf->alphasQ(q);
    // cout << "CT10 solution:                 " << setprecision(3) << setw(6) << as_ct10_q << endl;
    // fc << q << " " << as_ct10_q << endl;
    // const double as_ct10_q_2 = pdf->alphaS().alphasQ(q);
    // cout << "CT10 AlphaS solution:          " << setprecision(3) << setw(6)  << as_ct10_q_2
    //      << "    agrees = " << boolalpha << (as_ct10_q == as_ct10_q_2) << endl;
    // cout << endl;
  }
  fa.close(); fo.close(); fi.close(); fc.close();
 
  delete pdf;
#ifdef HAVE_MPI
  MPI_Finalize();
#endif  
  return 0;
}

Making a new PDF subclass in C++

#include "LHAPDF/PDF.h"
#include <iostream>
 
using namespace std;
 
 
struct AnalyticPDF : public LHAPDF::PDF {
 
  AnalyticPDF() {
    info().set_entry("Flavors", "-5,-4,-3,-2,-1,21,1,2,3,4,5");
  }
 
  double _xfxQ2(int id, double x, double q2) const {
    if (abs(id) > 5 && id != 21) return 0;
    return 0.15 * sin(20.0*x) * sin(20.0*q2);
  }
 
  void _xfxQ2(double x, double q2, std::vector<double>& ret) const {
    for (int id(-5); id<5; ++id)
      _xfxQ2(id,x,q2);
  }
 
  bool inRangeX(double x) const { return true; }
  bool inRangeQ2(double q2) const { return true; }
 
};
 
 
int main(int argc, const char* argv[]) {
  AnalyticPDF apdf;
  LHAPDF::PDF& pdf = apdf;
  for (double x = 0; x < 1.0; x += 0.1) {
    for (double logq2 = 1; logq2 < 6; logq2 += 0.5) {
      const double q2 = pow(10, logq2);
      cout << x << " " << q2 << " " << pdf.xfxQ2(21, x, q2) << endl;
    }
  }
  return 0;
}

Calculating PDF uncertainties in C++

// Program to test automatic calculation of PDF uncertainties.
//
// Written March 2014 by G. Watt <Graeme.Watt(at)durham.ac.uk>
// Extended April 2022 by A. Buckley <andy.buckley(at)cern.ch>
//
// Use formulae for PDF uncertainties and correlations in:
//   G. Watt, JHEP 1109 (2011) 069 [arXiv:1106.5788 [hep-ph]].
// Also see Section 6 of LHAPDF6 paper [arXiv:1412.7420 [hep-ph]].
// Extended in July 2015 for ErrorType values ending in "+as".
// Modified in September 2015 for more general ErrorType values:
// number of parameter variations determined by counting "+" symbols.
 
#include "LHAPDF/LHAPDF.h"
#include <random>
using namespace std;
 
 
// Helper function for computing and printing errors from the given vals vector
void printUncs(const LHAPDF::PDFSet& set, const vector<double>& vals, double cl, const string& varname, bool alt=false) {
  if (cl == 0) cl = LHAPDF::CL1SIGMA;
  cout << "PDF uncertainties on " << varname << " computed with " << set.name() << " at CL=" << cl << "%" << endl;
  const LHAPDF::PDFErrInfo errinfo = set.errorInfo();
  const LHAPDF::PDFUncertainty err = set.uncertainty(vals, cl, alt);
  if (cl >= 0) cout << "Scaled PDF uncertainties using scale = " << err.scale << endl;
  // Print summary numbers
  cout.precision(5);
  cout << varname << " = " << err.central << " +" << err.errplus << " -" << err.errminus << " (+-" << err.errsymm << ")" << endl;
  // Break down into quadrature-combined uncertainty components
  for (size_t i = 0; i < errinfo.qparts.size(); ++i) {
    //cout << "  " << errinfo.qpartName(i) << endl;
    cout << "  " << setw(12) << err.errparts[i].first << setw(12) << err.errparts[i].second << "  " << errinfo.qpartName(i) << endl;
  }
}
 
 
 
// Simple test program to demonstrate the PDFSet member functions.
//   set.errorInfo();
//   set.uncertainty(values, cl=68.268949..., alternative=false);
//   set.correlation(valuesA, valuesB);
//   set.randomValueFromHessian(values, randoms, symmetrise=true);
int main(int argc, char* argv[]) {
 
  if (argc < 2) {
    cerr << "You must specify a PDF set:  ./testpdfunc setname" << endl;
    return 1;
  }
 
  const string setname = argv[1];
  const LHAPDF::PDFSet set(setname);
  const size_t nmem = set.size()-1;
 
  double x = 0.1; // momentum fraction
  double Q = 100.0; // factorisation scale in GeV
 
  // Fill vectors xgAll and xuAll using all PDF members.
  // Could replace xg, xu by cross section, acceptance etc.
  const vector<LHAPDF::PDF*> pdfs = set.mkPDFs();
  vector<double> xgAll, xuAll;
  vector<string> pdftypes;
  for (size_t imem = 0; imem <= nmem; imem++) {
    xgAll.push_back(pdfs[imem]->xfxQ(21,x,Q)); // gluon distribution
    xuAll.push_back(pdfs[imem]->xfxQ(2,x,Q)); // up-quark distribution
    pdftypes.push_back(pdfs[imem]->type()); // PdfType of member
  }
 
  cout << endl;
  const LHAPDF::PDFErrInfo errinfo = set.errorInfo();
  cout << "ErrorType: " << errinfo.errtype << endl;
  cout << "ErrorConfLevel: " << errinfo.conflevel << endl;
 
  // Count number of parameter variations = number of '+' characters.
  const size_t npar = errinfo.nmemPar();
  if (npar > 0) cout << "Last " << npar << " members are parameter variations" << endl;
  cout << endl;
 
  // Check that the PdfType of each member matches the ErrorType of the set.
  // NB. "Hidden" expert-only functionality -- API may change
  set._checkPdfType(pdftypes);
 
 
  // Calculate PDF uncertainty on gluon distribution.
  cout << "Gluon distribution at Q = " << Q << " GeV (normal uncertainties)" << endl;
  printUncs(set, xgAll, -1, "xg"); //< -1 => same C.L. as set
  cout << endl;
 
  // Calculate PDF uncertainty on up-quark distribution.
  cout << "Up-quark distribution at Q = " << Q << " GeV (normal uncertainties)" << endl;
  printUncs(set, xuAll, -1, "xu"); //< -1 => same C.L. as set
  cout << endl;
 
  // Calculate sanity-check PDF self-correlation between gluon and gluon.
  const double autocorr = set.correlation(xgAll, xgAll);
  cout << "Self-correlation of xg = " << autocorr << endl;
  cout << endl;
 
  // Calculate PDF correlation between gluon and up-quark.
  // (This is the PDF-only correlation if npar > 0.)
  const double corr = set.correlation(xgAll, xuAll);
  cout << "Correlation between xg and xu = " << corr << endl;
  cout << endl;
 
 
  // Calculate gluon PDF uncertainty scaled to 90% C.L.
  cout << "Gluon distribution at Q = " << Q << " GeV (scaled uncertainties)" << endl;
  printUncs(set, xgAll, 90, "xg"); //< -1 => same C.L. as set
  cout << endl;
 
  // Calculate up-quark PDF uncertainty scaled to 1-sigma.
  cout << "Up-quark distribution at Q = " << Q << " GeV (scaled uncertainties)" << endl;
  printUncs(set, xuAll, 0, "xu"); //< -1 => same C.L. as set
  cout << endl;
 
 
  if (LHAPDF::startswith(set.errorType(), "replicas")) {
 
    // Calculate gluon PDF as median and 90% C.L. of replica probability distribution.
    cout << "Gluon distribution at Q = " << Q << " GeV (median and 90% C.L.)" << endl;
    printUncs(set, xgAll, 90, "xg", true);
    cout << endl;
 
    // Calculate up-quark PDF as median and 68% C.L. of replica probability distribution.
    cout << "Up-quark distribution at Q = " << Q << " GeV (median and 68% C.L.)" << endl;
    printUncs(set, xuAll, 68, "xu", true);
    cout << endl;
 
  } else if (LHAPDF::startswith(set.errorType(), "hessian") || LHAPDF::startswith(set.errorType(), "symmhessian")) {
 
    // Generate random values from Hessian best-fit and eigenvector values.
    // See: G. Watt and R.S. Thorne, JHEP 1208 (2012) 052 [arXiv:1205.4024 [hep-ph]].
 
    // If npar > 0 exclude the last 2*npar members (parameter variations).
    const int npdfmem = errinfo.nmemCore();
    const int neigen = (LHAPDF::startswith(set.errorType(), "hessian")) ? npdfmem/2 : npdfmem;
    const unsigned seed = 1234;
    default_random_engine generator(seed);
    normal_distribution<double> distribution; //< mean 0.0, s.d. = 1.0
    const int nrand = 5; // generate nrand random values
    for (int irand = 1; irand <= nrand; irand++) {
      // Fill vector "randoms" with neigen Gaussian random numbers.
      vector<double> randoms;
      for (int ieigen=1; ieigen <= neigen; ieigen++) {
        double r = distribution(generator);
        randoms.push_back(r);
      }
      // const bool symmetrise = false; // average differs from best-fit
      const bool symmetrise = true; // default: average tends to best-fit
      double xgrand = set.randomValueFromHessian(xgAll, randoms, symmetrise);
      // Pass *same* random numbers to preserve correlations between xg and xu.
      double xurand = set.randomValueFromHessian(xuAll, randoms, symmetrise);
      cout << "Random " << irand << ": xg = " << xgrand << ", xu = " << xurand << endl;
    }
    // Random values generated in this way can subsequently be used for
    // applications such as Bayesian reweighting or combining predictions
    // from different groups (as an alternative to taking the envelope).
    // See, for example, material at http://mstwpdf.hepforge.org/random/.
    // The "randomValueFromHessian" function is the basis of the program
    // "examples/hessian2replicas.cc" to convert a Hessian set to replicas.
    cout << endl;
 
  }
 
  return 0;
 
}