LHAPDF is hosted by Hepforge, IPPP Durham
LHAPDF  6.2.1
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>
using namespace std;
int main(int argc, char* argv[]) {
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
return 0;
}

Using and testing PDFs in Python

#! /usr/bin/env python
import lhapdf
p = lhapdf.mkPDF("CT10", 0)
p = lhapdf.mkPDF("CT10/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>
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[]) {
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.subgrid(21, 100).xf(0, 0) << endl;
cout << "x1, Q0 = " << pdf.subgrid(21, 100).xf(1, 0) << endl;
cout << "x0, Q1 = " << pdf.subgrid(21, 100).xf(0, 1) << endl;
cout << "x1, Q1 = " << pdf.subgrid(21, 100).xf(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;
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>
using namespace std;
int main() {
// 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;
}
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>
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() {
lookup(10800);
lookup(10801);
lookup(10042);
lookup(10041);
lookup(10799);
lookup(12346);
return 0;
}

Testing the path searching system in C++

// Test program for path searching machinery
#include "LHAPDF/Paths.h"
#include <iostream>
using namespace std;
int main() {
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;
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>
using namespace LHAPDF;
using namespace std;
int main() {
// 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;
return 0;
}

Making a new PDF subclass in C++

Calculating PDF uncertainties in C++

// Program to test LHAPDF6 automatic calculation of PDF uncertainties.
// Written in March 2014 by G. Watt <Graeme.Watt(at)durham.ac.uk>.
// 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;
// Simple test program to demonstrate the three PDFSet member functions.
// 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;
cout << "ErrorType: " << set.errorType() << endl;
cout << "ErrorConfLevel: " << set.errorConfLevel() << endl;
// Count number of parameter variations = number of '+' characters.
const size_t npar = LHAPDF::countchar(set.errorType(), '+');
if (npar > 0 )
cout << "Last " << 2*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);
// Define formats for printing labels and numbers in output.
string labformat = "%2s%10s%12s%12s%12s%12s\n";
string numformat = "%12.4e%12.4e%12.4e%12.4e%12.4e\n";
// Calculate PDF uncertainty on gluon distribution.
cout << "Gluon distribution at Q = " << Q << " GeV (normal uncertainties)" << endl;
printf(labformat.c_str()," #","x","xg","error+","error-","error");
const LHAPDF::PDFUncertainty xgErr = set.uncertainty(xgAll, -1); // -1 => same C.L. as set
printf(numformat.c_str(), x, xgErr.central, xgErr.errplus, xgErr.errminus, xgErr.errsymm);
if (npar > 0) {
// The last 2*npar members correspond to parameter variations (such as alphaS, etc.).
// In this case, the errors above are combined PDF+parameter uncertainties, obtained by
// adding in quadrature the PDF and parameter uncertainties, shown separately below.
printf(labformat.c_str()," #","x","error_par","error_pdf+","error_pdf-","error_pdf");
printf(numformat.c_str(), x, xgErr.err_par, xgErr.errplus_pdf, xgErr.errminus_pdf, xgErr.errsymm_pdf);
}
cout << endl;
// Calculate PDF uncertainty on up-quark distribution.
cout << "Up-quark distribution at Q = " << Q << " GeV (normal uncertainties)" << endl;
printf(labformat.c_str()," #","x","xu","error+","error-","error");
const LHAPDF::PDFUncertainty xuErr = set.uncertainty(xuAll, -1); // -1 => same C.L. as set
printf(numformat.c_str(), x, xuErr.central, xuErr.errplus, xuErr.errminus, xuErr.errsymm);
if (npar > 0) {
printf(labformat.c_str()," #","x","error_par","error_pdf+","error_pdf-","error_pdf");
printf(numformat.c_str(), x, xuErr.err_par, xuErr.errplus_pdf, xuErr.errminus_pdf, xuErr.errsymm_pdf);
}
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;
printf(labformat.c_str()," #","x","xg","error+","error-","error");
const LHAPDF::PDFUncertainty xgErr90 = set.uncertainty(xgAll, 90); // scale to 90% C.L.
printf(numformat.c_str(), x, xgErr90.central, xgErr90.errplus, xgErr90.errminus, xgErr90.errsymm);
cout << "Scaled PDF uncertainties to 90% C.L. using scale = " << xgErr90.scale << endl;
if (npar > 0) {
// If npar > 0 the same scale factor is applied to the parameter variation uncertainty.
printf(labformat.c_str()," #","x","error_par","error_pdf+","error_pdf-","error_pdf");
printf(numformat.c_str(), x, xgErr90.err_par, xgErr90.errplus_pdf, xgErr90.errminus_pdf, xgErr90.errsymm_pdf);
}
cout << endl;
// Calculate up-quark PDF uncertainty scaled to 1-sigma.
cout << "Up-quark distribution at Q = " << Q << " GeV (scaled uncertainties)" << endl;
printf(labformat.c_str()," #","x","xu","error+","error-","error");
const LHAPDF::PDFUncertainty xuErr1s = set.uncertainty(xuAll); // scale to 1-sigma (default)
printf(numformat.c_str(), x, xuErr1s.central, xuErr1s.errplus, xuErr1s.errminus, xuErr1s.errsymm);
cout << "Scaled PDF uncertainties to 1-sigma using scale = " << xuErr1s.scale << endl;
if (npar > 0) {
printf(labformat.c_str()," #","x","error_par","error_pdf+","error_pdf-","error_pdf");
printf(numformat.c_str(), x, xuErr1s.err_par, xuErr1s.errplus_pdf, xuErr1s.errminus_pdf, xuErr1s.errsymm_pdf);
}
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;
printf(labformat.c_str()," #","x","xg","error+","error-","error");
const LHAPDF::PDFUncertainty xgErr = set.uncertainty(xgAll, 90, true);
printf(numformat.c_str(), x, xgErr.central, xgErr.errplus, xgErr.errminus, xgErr.errsymm);
if (npar > 0) {
// If npar > 0 the parameter variation uncertainty is scaled to the requested C.L.
printf(labformat.c_str()," #","x","error_par","error_pdf+","error_pdf-","error_pdf");
printf(numformat.c_str(), x, xgErr.err_par, xgErr.errplus_pdf, xgErr.errminus_pdf, xgErr.errsymm_pdf);
}
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;
printf(labformat.c_str()," #","x","xu","error+","error-","error");
const LHAPDF::PDFUncertainty xuErr = set.uncertainty(xuAll, 68, true);
printf(numformat.c_str(), x, xuErr.central, xuErr.errplus, xuErr.errminus, xuErr.errsymm);
if (npar > 0) {
printf(labformat.c_str()," #","x","error_par","error_pdf+","error_pdf-","error_pdf");
printf(numformat.c_str(), x, xuErr.err_par, xuErr.errplus_pdf, xuErr.errminus_pdf, xuErr.errsymm_pdf);
}
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 = nmem - 2*npar;
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;
}