/usr/share/doc/python-lhapdf/examples/pdf-plot.py is in python-lhapdf 5.9.1-5.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 | #! /usr/bin/env python
import math, numpy
import lhapdf
import matplotlib.pyplot as plt
q = math.sqrt(6400.0)
pdfsets = ["cteq6ll.LHpdf", "MSTW2008lo68cl.LHgrid", "MRST2001lo.LHgrid", "MRST2007lomod.LHgrid", "MRSTMCal.LHgrid"]
#pdfsets = ["cteq6ll.LHpdf", "MRST2001lo.LHgrid", "MRST2007lomod.LHgrid", "MRSTMCal.LHgrid"]
partons = { 0 : "gluon", 2 : "up" }
NPOINTS = 1000
xs = numpy.logspace(-4, -0.001, NPOINTS)
plt.figure(figsize=(13,7))
for n, parton in enumerate(sorted(partons.keys())):
plt.subplot(1, len(partons), n+1)
lines = []
for pdfset in pdfsets:
lhapdf.initPDFSetByName(pdfset)
lhapdf.initPDF(0)
xfxs = numpy.zeros([NPOINTS])
for i, x in enumerate(xs):
xfx = lhapdf.xfx(x, q, parton)
xfxs[i] = xfx
l = plt.plot(xs, xfxs)
lines.append(l)
plt.xscale("log")
#plt.ylim(0.5, 3)
plt.legend(lines, pdfsets)
plt.title(partons[parton])
plt.xlabel("$x$")
if n == 0:
plt.ylabel("$x f(x, Q^2)$")
plt.show()
|