lhapdf

view plotpdfs @ 1110:2ba52dbf2312

Add release instructions for new sets
author Andy Buckley <andy@insectnation.org>
date Tue, 17 Oct 2017 17:12:04 +0100
parents 799f323cdbce
children
line source
1 #! /usr/bin/env python
3 """\
4 Usage: %prog <PDF1> [<PDF2> [...]]
6 Plot PDF and alpha_s values for the named LHAPDF6 sets.
8 TODO:
9 * Allow user specification of the various PDF/parton line colours and styles
10 * Allow user specification of the plot types to be shown
11 * Allow user specification of the discrete xs, Qs, and PIDs lists
12 """
14 import optparse, os, sys
16 op = optparse.OptionParser(usage=__doc__)
17 op.add_option("--xfmin", dest="XFMIN", metavar="NUM", help="minimum xf value [default: %default]", type=float, default=1e-2)
18 op.add_option("--xmin", dest="XMIN", metavar="NUM", help="minimum x value [default: %default]", type=float, default=1e-10)
19 op.add_option("--qmax", dest="QMAX", metavar="NUM", help="maximum Q value in GeV [default: %default]", type=float, default=1e4)
20 op.add_option("-f", "--format", dest="FORMAT", metavar="F", help="plot file format, i.e. file extension pdf/png/... [default: %default]", default="pdf")
21 op.add_option("--qs", dest="QS", metavar="Q1,Q2,...", help="discrete Q values to use on plots vs. x [default: %default]", default="1,10,100,1000,10000")
22 op.add_option("--xs", dest="XS", metavar="X1,X2,...", help="discrete x values to use on plots vs. Q [default: %default]", default="1e-5,1e-3,1e-2,1e-1")
23 op.add_option("--pids", dest="PIDS", metavar="ID1,ID2,...", help="PID values to use on PDF plots [default: %default]", default="0,1,2,3,4,5,-1,-2,-3,-4,-5")
24 op.add_option("--plots", dest="PLOTS", metavar="PLOT1,PLOT2,...", help="plot types to show, default value lists all types [default: %default]",
25 default="alphas,xf_x/pid,xf_x/q,xf_q/pid,xf_q/x")
26 op.add_option("-q", "--quiet", dest="VERBOSITY", action="store_const", const=0, help="suppress non-essential messages", default=1)
27 op.add_option("-v", "--verbose", dest="VERBOSITY", action="store_const", const=2, help="output debug messages", default=1)
28 opts, args = op.parse_args()
30 opts.PLOTS = opts.PLOTS.upper().split(",")
31 if not args:
32 print __doc__
33 sys.exit(1)
35 pnames = args
37 import matplotlib.pyplot as plt
38 plt.rcParams["font.family"] = "serif"
39 # plt.rcParams["font.serif"] = "Computer Modern Roman"
40 plt.rcParams["text.usetex"] = True
42 STYLES = ["-", "--", "-.", (0, (5,2,1,2,1,2)), ":"]
43 COLORS = ["red", "blue", "darkgreen", "orange", "purple", "magenta", "gray", "cyan"]
44 PARTONS = {-5:r"$\bar{b}$",-4:r"$\bar{c}$",-3:r"$\bar{s}$",-2:r"$\bar{u}$",-1:r"$\bar{d}$",
45 1:r"$d$",2:r"$u$",3:r"$s$",4:r"$c$",5:r"$b$",0:r"$g$"}
47 def tex_str(a):
48 return a.replace("_", r"\_").replace("#", r"\#")
50 def tex_float(f):
51 float_str = "{0:.2g}".format(f)
52 if "e" in float_str:
53 mant, exp = float_str.split("e")
54 exp = int(exp)
55 return r"{0} \times 10^{{{1}}}".format(mant, exp)
56 else:
57 return float_str
60 ## Get sampling points in x,Q
61 from math import log10
62 import numpy as np
63 xs = np.logspace(log10(opts.XMIN), 0, 100)
64 qs = np.logspace(0, log10(opts.QMAX), 100)
65 xs_few = [float(x) for x in opts.XS.split(",")] #[1e-5, 1e-3, 1e-2, 1e-1]
66 qs_few = [float(q) for q in opts.QS.split(",")] #[1, 10, 100, 1000, 10000]
67 pids = [int(i) for i in opts.PIDS.split(",")] #[0] + range(1,5+1) + [-i for i in range(1,5+1)]
68 #print xs_few, qs_few, pids
70 ## Load PDFs for plotting, indexed by name
71 import lhapdf
72 lhapdf.setVerbosity(opts.VERBOSITY)
73 pdfs = { pname : lhapdf.mkPDF(pname) for pname in pnames }
74 print
77 ## Make PDF xf vs. x & Q plots for each parton flavour, and a single alpha_s vs. Q plot
78 fig = plt.figure()
79 ax = fig.add_subplot(111)
81 ## alpha_s vs Q plot
82 if "ALPHAS" in opts.PLOTS:
83 plt.cla()
84 for i, pname in enumerate(pnames):
85 color = COLORS[i % len(COLORS)]
86 as_vals = [pdfs[pname].alphasQ(q) for q in qs]
87 ax.plot(qs, as_vals, label=tex_str(pname), color=color, ls="-")
88 ax.set_xlabel("$Q$")
89 ax.set_ylabel(r"$\alpha_s(Q)$")
90 ax.set_ylim(bottom=0)
91 ax.set_xscale("log")
92 l = ax.legend(loc=0, ncol=2, frameon=False, fontsize="xx-small")
93 fname = "alpha_s.{}".format(opts.FORMAT)
94 if opts.VERBOSITY > 0:
95 print "Writing plot file", fname
96 fig.savefig(fname)
98 ## xf vs. x plots (per PID)
99 if "XF_X/PID" in opts.PLOTS:
100 for pid in pids:
101 plt.cla()
102 ax.text(0.95, 0.5, "PID={:d}".format(pid), transform=ax.transAxes, ha="right", va="top")
103 ax.set_xlabel("$x$")
104 ax.set_ylabel("$xf(x,Q)$")
105 for i, pname in enumerate(pnames):
106 for j, q in enumerate(qs_few):
107 color = COLORS[i % len(COLORS)]
108 style = STYLES[j % len(STYLES)]
109 xf_vals = [pdfs[pname].xfxQ(pid, x,q) for x in xs]
110 # title = "{}, $Q={}$ ID={:d}".format(tex_str(pname), tex_float(q), pid)
111 title = "{}, $Q={}$".format(tex_str(pname), tex_float(q))
112 # print title
113 plt.plot(xs, xf_vals, label=title, color=color, ls=style, lw=1.0)
114 # if pid != 0:
115 # xf_vals = [pdfs[pname].xfxQ(-pid, x,q) for x in xs]
116 # plt.plot(xs, xf_vals, label="{}, $Q={}$, ID={:d}".format(pname, tex_float(q), -pid), color=color, ls=style, lw=0.5)
117 ax.set_xscale("log")
118 #ax.set_ylim(bottom=0)
119 ax.set_ylim(bottom=opts.XFMIN)
120 ax.set_yscale("log")
121 l = ax.legend(loc=0, ncol=2, frameon=False, fontsize="xx-small")
122 fname = "pdf_pid{:d}_x.{}".format(pid, opts.FORMAT)
123 if opts.VERBOSITY > 0:
124 print "Writing plot file", fname
125 fig.savefig(fname)
127 ## xf vs. x plots (per Q)
128 if "XF_X/Q" in opts.PLOTS:
129 for j, q in enumerate(qs_few):
130 plt.cla()
131 ax.text(0.95, 0.5, "$Q={}$".format(tex_float(q)), transform=ax.transAxes, ha="right", va="top")
132 ax.set_xlabel("$x$")
133 ax.set_ylabel("$xf(x,Q)$")
134 for pid in pids:
135 for i, pname in enumerate(pnames):
136 color = COLORS[pid % len(COLORS)]
137 style = STYLES[i % len(STYLES)]
138 xf_vals = [pdfs[pname].xfxQ(pid, x,q) for x in xs]
139 title = "{}, ID={:d}".format(tex_str(pname), pid)
140 plt.plot(xs, xf_vals, label=title, color=color, ls=style, lw=1.0)
141 # if pid != 0:
142 # xf_vals = [pdfs[pname].xfxQ(-pid, x,q) for x in xs]
143 # plt.plot(xs, xf_vals, label="{}, $Q={}$, ID={:d}".format(pname, tex_float(q), -pid), color=color, ls=style, lw=0.5)
144 ax.set_xscale("log")
145 #ax.set_ylim(bottom=0)
146 ax.set_ylim(bottom=opts.XFMIN)
147 ax.set_yscale("log")
148 l = ax.legend(loc=0, ncol=2, frameon=False, fontsize="xx-small")
149 fname = "pdf_q{:d}_x.{}".format(int(q), opts.FORMAT)
150 if opts.VERBOSITY > 0:
151 print "Writing plot file", fname
152 fig.savefig(fname)
154 ## xf vs. Q plots (per PID)
155 if "XF_Q/PID" in opts.PLOTS:
156 for pid in pids:
157 plt.cla()
158 ax.text(0.05, 0.7, "PID={:d}".format(pid), transform=ax.transAxes, ha="left", va="center")
159 ax.set_xlabel("$Q$")
160 ax.set_ylabel("$xf(x,Q)$")
161 for i, pname in enumerate(pnames):
162 for j, x in enumerate(xs_few):
163 color = COLORS[i % len(COLORS)]
164 style = STYLES[j % len(STYLES)]
165 xf_vals = [pdfs[pname].xfxQ(pid, x,q) for q in qs]
166 #print pname, x, xf_vals
167 #title = "{}, $x={}$ ID={:d}".format(tex_str(pname), tex_float(x), pid)
168 title = "{}, $x={}$".format(tex_str(pname), tex_float(x))
169 plt.plot(qs, xf_vals, label=title, color=color, ls=style, lw=1.0)
170 # if pid != 0:
171 # xf_vals = [pdfs[pname].xfxQ(-pid, x,q) for x in xs]
172 # plt.plot(xs, xf_vals, label="{}, $Q={}$, ID={:d}".format(pname, tex_float(q), -pid), color=color, ls=style, lw=0.5)
173 ax.set_xscale("log")
174 #ax.set_ylim(bottom=0)
175 ax.set_ylim(bottom=opts.XFMIN)
176 ax.set_yscale("log")
177 l = ax.legend(loc=0, ncol=2, frameon=False, fontsize="xx-small")
178 fname = "pdf_pid{:d}_q.{}".format(pid, opts.FORMAT)
179 if opts.VERBOSITY > 0:
180 print "Writing plot file", fname
181 fig.savefig(fname)
183 ## xf vs. Q plots (per x)
184 if "XF_Q/X" in opts.PLOTS:
185 for j, x in enumerate(xs_few):
186 plt.cla()
187 ax.text(0.95, 0.5, "$x={}$".format(tex_float(x)), transform=ax.transAxes, ha="right", va="bottom")
188 ax.set_xlabel("$Q$")
189 ax.set_ylabel("$xf(x,Q)$")
190 for pid in pids:
191 for i, pname in enumerate(pnames):
192 color = COLORS[pid % len(COLORS)]
193 style = STYLES[i % len(STYLES)]
194 xf_vals = [pdfs[pname].xfxQ(pid, x,q) for q in qs]
195 title = "{}, ID={:d}".format(tex_str(pname), pid)
196 plt.plot(qs, xf_vals, label=title, color=color, ls=style, lw=1.0)
197 # if pid != 0:
198 # xf_vals = [pdfs[pname].xfxQ(-pid, x,q) for x in xs]
199 # plt.plot(xs, xf_vals, label="{}, $Q={}$, ID={:d}".format(pname, tex_float(q), -pid), color=color, ls=style, lw=0.5)
200 ax.set_xscale("log")
201 #ax.set_ylim(bottom=0)
202 ax.set_ylim(bottom=opts.XFMIN)
203 ax.set_yscale("log")
204 l = ax.legend(loc=0, ncol=2, frameon=False, fontsize="xx-small")
205 fname = "pdf_x{:0.2f}_q.{}".format(x, opts.FORMAT)
206 if opts.VERBOSITY > 0:
207 print "Writing plot file", fname
208 fig.savefig(fname)
210 plt.close(fig)
211 print