lhapdf is hosted by Hepforge, IPPP Durham
LHAPDF  6.5.4
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

1 #! /usr/bin/env python
2 
3 ## Python LHAPDF6 usage example
4 
5 import lhapdf
6 
7 ## Getting a PDF member object
8 p = lhapdf.mkPDF("CT10nlo", 0)
9 p = lhapdf.mkPDF("CT10nlo/0")
10 
11 ## Gluon PDF querying at x=0.001, Q2=10000
12 print(p.xfxQ2(21, 1e-3, 1e4))
13 
14 ## Basic all-flavour PDF querying at x=0.01, Q=M_Z
15 for pid in p.flavors():
16  print(p.xfxQ(pid, 0.01, 91.2))
17 
18 # TODO: demonstrate looping over PDF set members
19 pset = lhapdf.getPDFSet("CT10nlo")
20 print(pset.description)
21 pcentral = pset.mkPDF(0)
22 pdfs1 = pset.mkPDFs()
23 pdfs2 = lhapdf.mkPDFs("CT10nlo") # a direct way to get all the set's PDFs
24 
25 ## Filling a numpy 2D array with xf(x,Q) samples
26 import numpy as np
27 xs = [x for x in np.logspace(-7, 0, 5)]
28 qs = [q for q in np.logspace(1, 4, 4)]
29 gluon_xfs = np.empty([len(xs), len(qs)])
30 for ix, x in enumerate(xs):
31  for iq, q in enumerate(qs):
32  gluon_xfs[ix,iq] = p.xfxQ(21, x, q)
33 print(gluon_xfs)
34 
35 ## Version info, search paths, and metadata
36 print(lhapdf.version())
37 print(lhapdf.__version__)
38 lhapdf.pathsPrepend("/path/to/extra/pdfsets")
39 print(lhapdf.paths())
40 # ...

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
// cout << "UndefFlavorAction: " << cfg.get_entry("UndefFlavorAction") << endl;
cout << "Verbosity: " << cfg.get_entry("Verbosity") << endl;
cfg.set_entry("Verbosity", 5);
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;
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:
// 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;
}