lhapdf
is hosted by
Hepforge
,
IPPP Durham
LHAPDF
6.5.4
Main page
PDF sets
Class hierarchy
Functions
Examples
More...
All
Classes
Namespaces
Files
Functions
Variables
Typedefs
Enumerations
Groups
Pages
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
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;
}
Generated by
1.8.5