LHAPDF  6.5.5
PDF.h
1 // -*- C++ -*-
2 //
3 // This file is part of LHAPDF
4 // Copyright (C) 2012-2024 The LHAPDF collaboration (see AUTHORS for details)
5 //
6 #pragma once
7 #ifndef LHAPDF_PDF_H
8 #define LHAPDF_PDF_H
9 
10 #include "LHAPDF/PDFInfo.h"
11 #include "LHAPDF/PDFIndex.h"
12 #include "LHAPDF/Factories.h"
13 #include "LHAPDF/AlphaS.h"
14 #include "LHAPDF/Utils.h"
15 #include "LHAPDF/Paths.h"
16 #include "LHAPDF/Exceptions.h"
17 #include "LHAPDF/Version.h"
18 #include "LHAPDF/Config.h"
19 
20 namespace LHAPDF {
21 
22 
23  /// @defgroup partonids Parton ID codes
24  /// @{
25  /// Convenience/clarity enum for specifying parton IDs
26  namespace PIDs { //< for scoping to avoid conflicts, e.g. PID::GLUON rather than just GLUON
27  enum PIDCode {
28  ATOP = -6, ABOTTOM = -5, ACHARM = -4, ASTRANGE = -3, AUP = -2, ADOWN = -1,
29  GLUON = 0, // equivalent to 21
30  DOWN = 1, UP = 2, STRANGE = 3, CHARM = 4, BOTTOM = 5, TOP = 6
31  };
32  }
33  /// @}
34 
35 
36 
37  /// @brief PDF is the general interface for access to parton density information.
38  ///
39  /// The PDF interface declares the general form of all PDF types, such as Grid based or analytic.
40  class PDF {
41  protected: //< These constructors should only be called by subclasses
42 
43  /// Internal convenience typedef for the AlphaS object handle
45 
46  /// Force initialization of the only non-class member.
47  PDF() : _forcePos(0) { }
48 
49 
50  public:
51 
52  /// Virtual destructor, to allow unfettered inheritance
53  virtual ~PDF() { }
54 
55 
56  protected:
57 
58 
59  /// @name Helper methods for info loading / path setting, used by derived types
60  ///@{
61 
62  void _loadInfo(const std::string& mempath);
63 
64  void _loadInfo(const std::string& setname, int member) {
65  const string searchpath = findpdfmempath(setname, member);
66  if (searchpath.empty())
67  throw UserError("Can't find a valid PDF " + setname + "/" + to_str(member));
68  _loadInfo(searchpath);
69  }
70 
71  void _loadInfo(int lhaid) {
72  const pair<string,int> setname_memid = lookupPDF(lhaid);
73  if (setname_memid.second == -1)
74  throw IndexError("Can't find a PDF with LHAPDF ID = " + to_str(lhaid));
75  _loadInfo(setname_memid.first, setname_memid.second);
76  }
77 
78  ///@}
79 
80 
81  public:
82 
83  /// @name PDF values
84  ///@{
85 
86  /// @brief Get the PDF xf(x) value at (x,q2) for the given PID.
87  ///
88  /// All grids are defined in Q2 rather than Q since the natural value
89  /// in MC programs is squared, so we typically avoid an expensive sqrt() call.
90  ///
91  /// @param id PDG parton ID
92  /// @param x Momentum fraction
93  /// @param q2 Squared energy (renormalization) scale
94  /// @return The value of xf(x,q2)
95  double xfxQ2(int id, double x, double q2) const;
96 
97 
98  /// @brief Get the PDF xf(x) value at (x,q) for the given PID.
99  ///
100  /// xfxQ will square the given q and return the value from xfxQ2.
101  /// All grids are defined in q2 rather than q since the natural value
102  /// in MC programs is squared, so we typically avoid an expensive sqrt() call.
103  ///
104  /// @param id PDG parton ID
105  /// @param x Momentum fraction
106  /// @param q Energy (renormalization) scale
107  /// @return The value of xf(x,q2)
108  double xfxQ(int id, double x, double q) const {
109  return xfxQ2(id, x, q*q);
110  }
111 
112 
113  /// @brief Get the PDF xf(x) value at (x,q2) for all supported PIDs.
114  ///
115  /// This version fills a user-supplied map to avoid container construction
116  /// costs on every call.
117  ///
118  /// @param x Momentum fraction
119  /// @param q2 Squared energy (renormalization) scale
120  /// @param rtn Map of PDF xf(x,q2) values, to be filled
121  void xfxQ2(double x, double q2, std::map<int, double>& rtn) const;
122 
123 
124  /// @brief Get the PDF xf(x) value at (x,q) for all supported PIDs.
125  ///
126  /// This version fills a user-supplied map to avoid container construction
127  /// costs on every call.
128  ///
129  /// @param x Momentum fraction
130  /// @param q Energy (renormalization) scale
131  /// @param rtn Map of PDF xf(x,q) values, to be filled
132  void xfxQ(double x, double q, std::map<int, double>& rtn) const {
133  xfxQ2(x, q*q, rtn);
134  }
135 
136 
137  /// @brief Get the PDF xf(x) value at (x,q2) for "standard" PIDs.
138  ///
139  /// This version fills a user-supplied vector to avoid container
140  /// construction costs on every call.
141  ///
142  /// The filled vector follows the LHAPDF5 convention, with 13 entries
143  /// running in the PDF ID order [-6, -5, ..., -1, 21, 1, ... 5, 6], i.e.
144  /// quark PDF values will be at vector index pid+6 and the gluon at index 6.
145  ///
146  /// @param x Momentum fraction
147  /// @param q2 Squared energy (renormalization) scale
148  /// @param rtn Vector of PDF xf(x,q2) values, to be filled
149  void xfxQ2(double x, double q2, std::vector<double>& rtn) const;
150 
151 
152  /// @brief Get the PDF xf(x) value at (x,q) for "standard" PIDs.
153  ///
154  /// This version fills a user-supplied vector to avoid container
155  /// construction costs on every call.
156  ///
157  /// The filled vector follows the LHAPDF5 convention, with 13 entries
158  /// running in the PDF ID order [-6, -5, ..., -1, 21, 1, ... 5, 6], i.e.
159  /// quark PDF values will be at vector index pid+6 and the gluon at index 6.
160  ///
161  /// @param x Momentum fraction
162  /// @param q Energy (renormalization) scale
163  /// @param rtn Vector of PDF xf(x,q) values, to be filled
164  void xfxQ(double x, double q, std::vector<double>& rtn) const {
165  xfxQ2(x, q*q, rtn);
166  }
167 
168 
169  /// @brief Get the PDF xf(x) value at (x,q2) for all supported PIDs.
170  ///
171  /// This version creates a new map on every call: prefer to use the
172  /// fill-in-place version with a user-supplied map for many calls.
173  ///
174  /// @param x Momentum fraction
175  /// @param q2 Squared energy (renormalization) scale
176  /// @return A map of PDF xf(x,q2) values
177  std::map<int, double> xfxQ2(double x, double q2) const;
178 
179  /// @brief Get the PDF xf(x) value at (x,q) for all supported PIDs.
180  ///
181  /// This version creates a new map on every call: prefer to use the
182  /// fill-in-place version with a user-supplied map for many calls.
183  ///
184  /// xfxQ will square the given q and return the value from xfxQ2.
185  /// All grids are defined in q2 rather than q since the natural value
186  /// in MC programs is squared, so we typically avoid an expensive sqrt() call.
187  ///
188  /// @param x Momentum fraction
189  /// @param q Energy (renormalization) scale
190  /// @return A map of PDF xf(x,q) values
191  std::map<int, double> xfxQ(double x, double q) const {
192  return xfxQ2(x, q*q);
193  }
194 
195 
196  protected:
197 
198  /// @brief Calculate the PDF xf(x) value at (x,q2) for the given PID.
199  ///
200  /// This is the key function to be overridden in concrete PDF types, since
201  /// it actually does the calculation of xf(x,Q2) by analytic, interpolation,
202  /// or other means. The user-called xfxQ2 method exists so that the physical
203  /// range and PID checks need only be done in one place, rather than need to
204  /// be re-implemented in each concrete implementation.
205  ///
206  /// @param id Parton ID in the PDG scheme
207  /// @param x Momentum fraction
208  /// @param q2 Squared energy (renormalization) scale
209  /// @return the value of xf(x,q2)
210  virtual double _xfxQ2(int id, double x, double q2) const = 0;
211 
212  virtual void _xfxQ2(double x, double q2, std::vector<double>& ret) const = 0;
213 
214  ///@}
215 
216 
217  public:
218 
219  /// @name Ranges of validity
220  ///@{
221 
222  /// Minimum valid x value for this PDF.
223  virtual double xMin() {
224  if (info().has_key("XMin"))
225  return info().get_entry_as<double>("XMin");
226  return numeric_limits<double>::epsilon();
227  }
228 
229  /// Maximum valid x value for this PDF.
230  virtual double xMax() {
231  if (info().has_key("XMax"))
232  return info().get_entry_as<double>("XMax");
233  return 1.0;
234  }
235 
236  /// Minimum valid Q value for this PDF (in GeV).
237  /// @note This function calls sqrt(q2Min()). For better CPU efficiency and accuracy use q2Min() directly.
238  virtual double qMin() {
239  return info().get_entry_as<double>("QMin", 0);
240  }
241 
242  /// @brief Maximum valid Q value for this PDF (in GeV).
243  /// @note This function calls sqrt(q2Max()). For better CPU efficiency and accuracy use q2Max() directly.
244  virtual double qMax() {
245  return info().get_entry_as<double>("QMax", numeric_limits<double>::max());
246  }
247 
248  /// Minimum valid Q2 value for this PDF (in GeV2).
249  virtual double q2Min() {
250  return sqr(this->qMin());
251  }
252 
253  /// Maximum valid Q2 value for this PDF (in GeV2).
254  virtual double q2Max() {
255  // Explicitly re-access this from the info, to avoid an overflow from squaring double_max
256  return (info().has_key("QMax")) ? sqr(info().get_entry_as<double>("QMax")) : numeric_limits<double>::max();
257  }
258 
259 
260  /// @brief Check whether the PDF is set to only return positive (definite) values or not.
261  ///
262  /// This is to avoid overshooting in to negative values when
263  /// interpolating/extrapolating PDFs that sharply decrease towards zero.
264  /// 0 = unforced, 1 = forced positive, 2 = forced positive definite (>= 1e-10)
265  int forcePositive() const {
266  if (_forcePos < 0) //< Caching
267  _forcePos = info().get_entry_as<unsigned int>("ForcePositive", 0);
268  return _forcePos;
269  }
270  /// @brief Set whether the PDF will only return positive (definite) values or not.
271  void setForcePositive(int mode) {
272  _forcePos = mode;
273  }
274 
275 
276  /// @brief Check whether the given x is physically valid
277  ///
278  /// Returns false for x less than 0 or greater than 1, since it
279  /// is a momentum fraction and not valid outside those values.
280  bool inPhysicalRangeX(double x) const {
281  return x >= 0.0 && x <= 1.0;
282  }
283 
284  /// @brief Check whether the given Q2 is physically valid
285  ///
286  /// Returns false for Q2 less than 0 (Q must be real-valued).
287  bool inPhysicalRangeQ2(double q2) const {
288  return q2 >= 0.0;
289  }
290 
291  /// @brief Check whether the given Q is physically valid
292  ///
293  /// Returns false for Q less than 0 (Q must be positive).
294  bool inPhysicalRangeQ(double q) const {
295  return inPhysicalRangeQ2(q*q);
296  }
297 
298  /// Check whether the given (x,Q2) is physically valid
299  bool inPhysicalRangeXQ2(double x, double q2) const {
301  }
302 
303  /// Check whether the given (x,Q) is physically valid
304  bool inPhysicalRangeXQ(double x, double q) const {
306  }
307 
308  /// @brief Grid range check for Q
309  ///
310  /// Return true when given Q is in the coverage range of this PDF.
311  /// It actually squares the given Q and returns value from inRangeQ2.
312  ///
313  /// @param q Energy scale
314  /// @return Whether q is in range
315  virtual bool inRangeQ(double q) const {
316  return inRangeQ2(q*q);
317  }
318 
319  /// @brief Grid range check for Q2
320  ///
321  /// Return true when given Q2 is in the coverage range of this PDF.
322  ///
323  /// @param q2 Squared energy scale
324  /// @return Whether q2 is in range
325  virtual bool inRangeQ2(double q2) const = 0;
326 
327  /// @brief Grid range check for x
328  ///
329  /// Return true when given x is in the coverage range of this PDF.
330  ///
331  /// @param x Momentum fraction
332  /// @return Whether x is in range
333  virtual bool inRangeX(double x) const = 0;
334 
335  /// Combined range check for x and Q
336  virtual bool inRangeXQ(double x, double q) const {
337  return inRangeX(x) && inRangeQ(q);
338  }
339 
340  /// Combined range check for x and Q2
341  bool inRangeXQ2(double x, double q2) const {
342  return inRangeX(x) && inRangeQ2(q2);
343  }
344 
345  ///@}
346 
347 
348  /// @name Generic member-level metadata (including cascaded metadata from set & config level)
349  ///@{
350 
351  /// Get the info class that actually stores and handles the metadata
352  PDFInfo& info() { return _info; }
353 
354  /// Get the info class that actually stores and handles the metadata (const version)
355  const PDFInfo& info() const { return _info; }
356 
357  /// @brief Get the PDF set of which this is a member
358  ///
359  /// Obtained from the member file path, not Info-based metadata.
360  PDFSet& set() const {
361  return getPDFSet(_setname());
362  }
363 
364  ///@}
365 
366 
367  /// @name Member-level metadata
368  ///@{
369 
370  /// @brief PDF member local ID number
371  ///
372  /// Obtained from the member file path, not Info-based metadata.
373  int memberID() const {
374  const string memname = file_stem(_mempath);
375  assert(memname.length() > 5); // There must be more to the filename stem than just the _nnnn suffix
376  const int memid = lexical_cast<int>(memname.substr(memname.length()-4)); //< Last 4 chars should be the member number
377  return memid;
378  }
379 
380  /// @brief PDF member global LHAPDF ID number
381  ///
382  /// Obtained from the member ID and the set's LHAPDF ID index
383  int lhapdfID() const;
384 
385  /// Description of this PDF member
386  std::string description() const {
387  return info().get_entry("MemDesc", info().get_entry("PdfDesc", ""));
388  }
389 
390  /// Version of this PDF's data file
391  int dataversion() const {
392  return info().get_entry_as<int>("DataVersion", -1);
393  }
394 
395  /// Get the type of PDF member that this object represents (central, error)
396  std::string type() const {
397  return to_lower(info().get_entry("MemType", info().get_entry("PdfType")));
398  }
399 
400  ///@}
401 
402 
403  /// Summary printout
404  void print(std::ostream& os=std::cout, int verbosity=1) const;
405 
406 
407  /// @name Parton content and QCD parameters
408  ///@{
409 
410  /// @brief List of flavours defined by this PDF set.
411  ///
412  /// This list is stored locally after its initial read from the Info object
413  /// to avoid unnecessary lookups and string decoding, since e.g. it is
414  /// looked at by every call to the GridPDF's Interpolator and Extrapolator
415  /// classes.
416  ///
417  /// @todo Make virtual for AnalyticPDF? Or allow manual setting of the Info?
418  virtual const std::vector<int>& flavors() const {
419  if (_flavors.empty()) {
420  _flavors = info().get_entry_as< vector<int> >("Flavors");
421  sort(_flavors.begin(), _flavors.end());
422  }
423  return _flavors;
424  }
425 
426  /// @brief Manually set/override the list of flavours defined by this PDF set.
427  ///
428  /// @note The provided vector will be internally sorted into ascending numeric order.
429  void setFlavors(std::vector<int> const& flavors) {
430  _flavors = flavors;
431  sort(_flavors.begin(), _flavors.end());
432  }
433 
434  /// Checks whether @a id is a valid parton for this PDF.
435  bool hasFlavor(int id) const;
436 
437  /// @brief Order of QCD at which this PDF has been constructed
438  ///
439  /// "Order" is defined here and throughout LHAPDF as the maximum number of
440  /// loops included in the matrix elements, in order to have an integer value
441  /// for easy use in comparisons, as opposed to "LO", "NLO", etc. strings.
442  ///
443  /// @todo Provide a setter function?
444  int orderQCD() const {
445  return info().get_entry_as<int>("OrderQCD");
446  }
447  /// @deprecated Use orderQCD instead
448  int qcdOrder() const { return orderQCD(); }
449 
450  /// @brief Get a quark mass in GeV by PDG code (|PID| = 1-6 only)
451  ///
452  /// Convenience interface to the Mass* info keywords.
453  /// Returns -1 for an undefined PID.
454  ///
455  /// @todo Provide a setter function?
456  double quarkMass(int id) const;
457 
458  /// @brief Get a flavor scale threshold in GeV by PDG code (|PID| = 1-6 only)
459  /// Convenience interface to the Mass* and Threshold* info keywords.
460  /// Returns -1 for an undefined PID.
461  ///
462  /// @todo Provide a setter function?
463  double quarkThreshold(int id) const;
464 
465  ///@}
466 
467 
468  /// @name QCD running coupling calculation
469  ///@{
470 
471  /// @brief Set the AlphaS calculator by pointer
472  ///
473  /// The provided AlphaS must have been new'd, as it will not be copied and
474  /// ownership passes to this GridPDF: it will be deleted when this PDF goes
475  /// out of scope or another setAlphaS call is made.
476  void setAlphaS(AlphaS* alphas) {
477  _alphas.reset(alphas);
478  }
479 
480  /// @brief Set the AlphaS calculator by smart pointer
481  void setAlphaS(AlphaSPtr alphas) {
482  _alphas = std::move(alphas);
483  }
484 
485  /// @brief Check if an AlphaS calculator is set
486  bool hasAlphaS() const {
487  return bool(_alphas);
488  }
489 
490  /// @brief Retrieve the AlphaS object for this PDF
492  return *_alphas;
493  }
494 
495  /// @brief Retrieve the AlphaS object for this PDF (const)
496  const AlphaS& alphaS() const {
497  return *_alphas;
498  }
499 
500  /// @brief Value of alpha_s(Q2) used by this PDF
501  ///
502  /// Calculated numerically, analytically, or interpolated according to
503  /// metadata, using the AlphaS classes.
504  double alphasQ(double q) const {
505  return alphasQ2(q*q);
506  }
507 
508  /// @brief Value of alpha_s(Q2) used by this PDF
509  ///
510  /// Calculated numerically, analytically, or interpolated according to
511  /// metadata, using the AlphaS classes.
512  double alphasQ2(double q2) const {
513  if (!hasAlphaS()) throw Exception("No AlphaS pointer has been set");
514  return _alphas->alphasQ2(q2);
515  }
516 
517  ///@}
518 
519 
520  protected:
521 
522  void _loadAlphaS() {
523  _alphas.reset( mkAlphaS(info()) );
524  }
525 
526  /// Get the set name from the member data file path (for internal use only)
527  std::string _setname() const {
528  return basename(dirname(_mempath));
529  }
530 
531  /// Member data file path
533 
534  /// Metadata container
536 
537  /// Locally cached list of supported PIDs (mutable for laziness/caching)
538  mutable vector<int> _flavors;
539 
540  /// Optionally loaded AlphaS object (mutable for laziness/caching)
542 
543  /// @brief Cached flag for whether to return only positive (or positive definite) PDF values
544  ///
545  /// A negative value indicates that the flag has not been set. 0 = no
546  /// forcing, 1 = force positive (i.e. 0 is permitted, negative values are
547  /// not), 2 = force positive definite (i.e. no values less than 1e-10).
548  mutable int _forcePos;
549 
550  };
551 
552 
553 }
554 #endif