"""
Module containing an LHAPDF class compatible with validphys
using the official lhapdf python interface.
The ``.members`` and ``.central_member`` of the ``LHAPDFSet`` are
LHAPDF objects (the typical output from ``mkPDFs``) and can be used normally.
Examples
--------
>>> from validphys.lhapdfset import LHAPDFSet
>>> pdf = LHAPDFSet("NNPDF40_nnlo_as_01180", "replicas")
>>> len(pdf.members)
101
>>> pdf.central_member.alphasQ(91.19)
0.11800
>>> pdf.members[0].xfxQ2(0.5, 15625)
{-5: 6.983360500601136e-05,
-4: 0.0021818063617227604,
-3: 0.00172453472243952,
-2: 0.0010906577230485718,
-1: 0.0022049272225017286,
1: 0.020051104853608722,
2: 0.0954139944889494,
3: 0.004116641378803191,
4: 0.002180124185625795,
5: 6.922722705177504e-05,
21: 0.007604124516892057}
"""
import logging
import numpy as np
from validphys.lhapdf_compatibility import make_pdf
log = logging.getLogger(__name__)
[docs]
class LHAPDFSet:
"""Wrapper for the lhapdf python interface.
Once instantiated this class will load the PDF set from LHAPDF.
If it is a T0 set only the CV will be loaded.
"""
def __init__(self, name, error_type):
self._name = name
self._error_type = error_type
if self.is_t0:
# If at this point we already know this is a T0 set, load only the CV
self._lhapdf_set = make_pdf(name, 0)
else:
self._lhapdf_set = make_pdf(name)
self._flavors = None
@property
def is_t0(self):
"""Check whether we are in t0 mode"""
return self._error_type == "t0"
@property
def n_members(self):
"""Return the number of active members in the PDF set"""
return len(self.members)
@property
def members(self):
"""Return the members of the set
the special error type t0 returns only member 0
"""
if self.is_t0:
return self._lhapdf_set[0:1]
return self._lhapdf_set
@property
def central_member(self):
"""Returns a reference to member 0 of the PDF list"""
return self._lhapdf_set[0]
[docs]
def xfxQ(self, x, Q, n, fl):
"""Return the PDF value for one single point for one single member
If the flavour is not included in the PDF (for instance top/antitop) return 0.0
"""
member_pdf = self.members[n]
res = member_pdf.xfxQ(x, Q)
return res.get(fl, 0.0)
@property
def flavors(self):
"""Returns the list of accepted flavors by the LHAPDF set"""
if self._flavors is None:
self._flavors = self.members[0].flavors()
return self._flavors
[docs]
def grid_values(self, flavors: np.ndarray, xgrid: np.ndarray, qgrid: np.ndarray):
"""Returns the PDF values for every member for the required
flavours, points in x and pointx in q
The return shape is
(members, flavors, xgrid, qgrid)
Return
------
ndarray of shape (members, flavors, xgrid, qgrid)
Examples
--------
>>> import numpy as np
>>> from validphys.lhapdfset import LHAPDFSet
>>> pdf = LHAPDFSet("NNPDF40_nnlo_as_01180", "replicas")
>>> xgrid = np.random.rand(10)
>>> qgrid = np.random.rand(3)
>>> flavs = np.arange(-4,4)
>>> flavs[4] = 21
>>> results = pdf.grid_values(flavs, xgrid, qgrid)
"""
# Create an array of x and q of equal length for LHAPDF
xarr, qarr = (g.ravel() for g in np.meshgrid(xgrid, qgrid))
# Ask LHAPDF for the values and swap the flavours and xgrid-qgrid axes
raw = np.array([member.xfxQ(flavors, xarr, qarr) for member in self.members]).swapaxes(1, 2)
# Unroll the xgrid-qgrid axes
return raw.reshape(self.n_members, len(flavors), len(xgrid), len(qgrid))