PDFxTMDLib  1.0.0
PDFSet.h
Go to the documentation of this file.
1 #pragma once
15 #include <PDFxTMDLib/Factory.h>
25 #include <memory>
26 #include <mutex>
27 #include <string>
28 #include <type_traits>
29 #include <vector>
30 
31 namespace PDFxTMD
32 {
35 template <typename T> struct PDFType;
36 
38 template <> struct PDFType<TMDPDFTag>
39 {
40  using type = ITMD;
41 };
42 
44 template <> struct PDFType<CollinearPDFTag>
45 {
46  using type = ICPDF;
47 };
48 
59 template <typename Tag> class PDFSet
60 {
61  public:
63  using PDF_t = typename PDFType<Tag>::type;
64 
70  explicit PDFSet(std::string pdfSetName, bool alternativeReplicaUncertainty = false)
71  : m_pdfSetName(std::move(pdfSetName)),
72  m_alternativeReplicaUncertainty(alternativeReplicaUncertainty),
73  m_uncertaintyStrategy_(NullUncertaintyStrategy()),
74  m_qcdCoupling(CouplingFactory().mkCoupling(m_pdfSetName))
75  {
76  Initialize();
78  }
79 
81  PDFSet():m_uncertaintyStrategy_(NullUncertaintyStrategy()), m_qcdCoupling(NullQCDCoupling()){};
82 
88  double alphasQ(double q) const
89  {
90  return alphasQ2(q * q);
91  }
92 
98  double alphasQ2(double q2) const
99  {
100  return m_qcdCoupling.AlphaQCDMu2(q2);
101  }
102 
108  PDF_t *operator[](int member)
109  {
110  std::lock_guard<std::mutex> lock(m_pdfSetMtx);
111  if (m_PDFSet_.find(member) == m_PDFSet_.end())
112  {
113  CreatePDFSet(member);
114  }
115  return m_PDFSet_[member].get();
116  }
117 
123  PDF_t *operator[](int member) const
124  {
125  if (m_PDFSet_.find(member) == m_PDFSet_.end())
126  return nullptr;
127  return m_PDFSet_.at(member).get();
128  }
129 
139  template <typename T = Tag, typename = std::enable_if_t<std::is_same_v<T, TMDPDFTag>>>
140  void Uncertainty(PartonFlavor flavor, double x, double kt2, double mu2, double cl,
141  PDFUncertainty &resUncertainty)
142  {
143  const auto pdfs = CalculatePDFValues(flavor, x, kt2, mu2);
144  PDFUncertaintyInternalEvaluation(pdfs, cl, resUncertainty);
145  }
146 
155  template <typename T = Tag, typename = std::enable_if_t<std::is_same_v<T, CollinearPDFTag>>>
156  void Uncertainty(PartonFlavor flavor, double x, double mu2, double cl,
157  PDFUncertainty &resUncertainty)
158  {
159  const auto pdfs = CalculatePDFValues(flavor, x, mu2);
160  PDFUncertaintyInternalEvaluation(pdfs, cl, resUncertainty);
161  }
162 
172  template <typename T = Tag, typename = std::enable_if_t<std::is_same_v<T, TMDPDFTag>>>
173  PDFUncertainty Uncertainty(PartonFlavor flavor, double x, double kt2, double mu2,
174  double cl = NO_REQUESTED_CONFIDENCE_LEVEL)
175  {
176  PDFUncertainty resUncertainty;
177  const auto pdfs = CalculatePDFValues(flavor, x, kt2, mu2);
178  PDFUncertaintyInternalEvaluation(pdfs, cl, resUncertainty);
179  return resUncertainty;
180  }
181 
190  template <typename T = Tag, typename = std::enable_if_t<std::is_same_v<T, CollinearPDFTag>>>
191  PDFUncertainty Uncertainty(PartonFlavor flavor, double x, double mu2,
192  double cl = NO_REQUESTED_CONFIDENCE_LEVEL)
193  {
194  PDFUncertainty resUncertainty;
195  const auto pdfs = CalculatePDFValues(flavor, x, mu2);
196  PDFUncertaintyInternalEvaluation(pdfs, cl, resUncertainty);
197  return resUncertainty;
198  }
199 
206  void Uncertainty(const std::vector<double> &values, double cl, PDFUncertainty &resUncertainty)
207  {
208  if (values.size() != m_pdfSetErrorInfo.size)
209  throw InvalidInputError("Error in PDFxTMD::PDFSet::Uncertainty. Input vector must "
210  "contain values for all PDF members.");
211  PDFUncertaintyInternalEvaluation(values, cl, resUncertainty);
212  }
213 
220  PDFUncertainty Uncertainty(const std::vector<double> &values,
221  double cl = NO_REQUESTED_CONFIDENCE_LEVEL)
222  {
223  if (values.size() != m_pdfSetErrorInfo.size)
224  throw InvalidInputError("Error in PDFxTMD::PDFSet::Uncertainty. Input vector must "
225  "contain values for all PDF members.");
226  PDFUncertainty resUncertainty;
227  PDFUncertaintyInternalEvaluation(values, cl, resUncertainty);
228  return resUncertainty;
229  }
230 
241  template <typename T = Tag, typename = std::enable_if_t<std::is_same_v<T, CollinearPDFTag>>>
242  double Correlation(PartonFlavor flavorA, double xA, double mu2A, PartonFlavor flavorB,
243  double xB, double mu2B)
244  {
245  const auto pdfsA = CalculatePDFValues(flavorA, xA, mu2A);
246  const auto pdfsB = CalculatePDFValues(flavorB, xB, mu2B);
247  return m_uncertaintyStrategy_.Correlation(pdfsA, pdfsB, m_pdfErrInfo.nmemCore());
248  }
249 
262  template <typename T = Tag, typename = std::enable_if_t<std::is_same_v<T, TMDPDFTag>>>
263  double Correlation(PartonFlavor flavorA, double xA, double kt2A, double mu2A,
264  PartonFlavor flavorB, double xB, double kt2B, double mu2B)
265  {
266  const auto pdfsA = CalculatePDFValues(flavorA, xA, kt2A, mu2A);
267  const auto pdfsB = CalculatePDFValues(flavorB, xB, kt2B, mu2B);
268  return m_uncertaintyStrategy_.Correlation(pdfsA, pdfsB, m_pdfErrInfo.nmemCore());
269  }
270 
277  double Correlation(const std::vector<double> &valuesA, const std::vector<double> &valuesB) const
278  {
279  if (valuesA.size() != m_pdfSetErrorInfo.size || valuesB.size() != m_pdfSetErrorInfo.size)
280  throw InvalidInputError("Error in PDFxTMD::PDFSet::Correlation. Input vectors must "
281  "contain values for all PDF members.");
282  return m_uncertaintyStrategy_.Correlation(valuesA, valuesB, m_pdfErrInfo.nmemCore());
283  }
284 
289  void CreatePDFSet(unsigned int setMember)
290  {
291  {
292  if constexpr (std::is_same_v<Tag, TMDPDFTag>)
293  {
294  m_PDFSet_.insert_or_assign(
295  setMember, std::make_unique<PDF_t>(
296  PDFxTMD::GenericTMDFactory().mkTMD(m_pdfSetName, setMember)));
297  }
298  else if constexpr (std::is_same_v<Tag, CollinearPDFTag>)
299  {
300  m_PDFSet_.insert_or_assign(
301  setMember, std::make_unique<PDF_t>(
302  PDFxTMD::GenericCPDFFactory().mkCPDF(m_pdfSetName, setMember)));
303  }
304  else
305  {
306  static_assert(!std::is_same_v<Tag, Tag>, "Unsupported Tag");
307  }
308  }
309  }
310 
315  {
316  for (int i = 0; i < m_pdfSetStdInfo.NumMembers; ++i)
317  {
318  if (m_PDFSet_.find(i) == m_PDFSet_.end()) {
319  CreatePDFSet(i);
320  }
321  }
322  }
323 
328  void InitailizePDFSetName(std::string pdfSetName)
329  {
330  m_pdfSetName = std::move(pdfSetName);
331  Initialize();
332  }
333 
338  size_t size() const
339  {
340  return m_pdfSetStdInfo.NumMembers;
341  }
342 
348  {
349  return m_pdfSetStdInfo;
350  }
351 
357  {
358  return m_pdfSetErrorInfo;
359  }
360 
366  {
367  return m_pdfSetInfo;
368  }
369 
370  // Rule of Five: disable copying, allow moving.
371  PDFSet(const PDFSet &) = delete;
372  PDFSet &operator=(const PDFSet &) = delete;
373  PDFSet(PDFSet &&) = default;
374  PDFSet &operator=(PDFSet &&) = default;
375 
376  private:
383  inline void PDFUncertaintyInternalEvaluation(const std::vector<double> &pdfs, double cl,
384  PDFUncertainty &resUncertainty)
385  {
386  const double reqCL = ValidateAndGetCL(cl);
387 
388  m_uncertaintyStrategy_.Uncertainty(pdfs, m_pdfErrInfo.nmemCore(), reqCL, resUncertainty);
389  resUncertainty.central = pdfs[0];
390 
391  ApplyConfidenceLevelScaling(resUncertainty, reqCL);
392  StoreCoreVariationErros(resUncertainty);
393  CalculateParameterVariationErrors(resUncertainty, pdfs);
394  }
395 
397  void Initialize()
398  {
399  ValidatePDFSetName();
400  LoadYamlInfo();
401  InitializeUncertaintyStrategy();
402  }
403 
405  void InitializeQCDCoupling()
406  {
407  m_qcdCoupling = CouplingFactory().mkCoupling(m_pdfSetName);
408  }
409 
411  void ValidatePDFSetName()
412  {
413  if (m_pdfSetName.empty())
414  {
415  throw InvalidInputError(
416  "PDFxTMD::PDFSet::ValidatePDFSetName: PDF set name cannot be empty.");
417  }
418  }
419 
421  void LoadYamlInfo()
422  {
423  auto infoPathPair = StandardInfoFilePath(m_pdfSetName);
424  if (infoPathPair.second != ErrorType::None)
425  throw FileLoadException(
426  "PDFxTMD::PDFSet::LoadYamlInfo: Unable to find info file of PDF set " +
427  m_pdfSetName);
428 
429  m_isValid = true;
430  auto yamlInfoErrPair = YamlErrorInfoReader(*infoPathPair.first);
431  if (yamlInfoErrPair.second != ErrorType::None || !yamlInfoErrPair.first.has_value())
432  {
433  m_isValid = false;
434  }
435  else
436  {
437  m_pdfSetErrorInfo = *yamlInfoErrPair.first;
438  }
439 
440  auto yamlStdInfoErrPair = YamlStandardPDFInfoReader(*infoPathPair.first);
441  if (yamlStdInfoErrPair.second != ErrorType::None || !yamlStdInfoErrPair.first.has_value())
442  {
443  m_isValid = false;
444  }
445  else
446  {
447  m_pdfSetStdInfo = *yamlStdInfoErrPair.first;
448  }
449 
450  if (!m_isValid)
451  {
452  throw InvalidInputError("PDFSet: Invalid PDF set '" + m_pdfSetName +
453  "'. Please check the YAML file.");
454  }
455  if (!m_pdfSetInfo.loadFromFile(*infoPathPair.first, ConfigWrapper::Format::YAML))
456  {
457  throw FileLoadException("path: " + *infoPathPair.first + "not found.");
458  }
459  }
460 
462  void InitializeUncertaintyStrategy()
463  {
464  m_pdfErrInfo = PDFErrInfo::CalculateErrorInfo(m_pdfSetErrorInfo);
465  if (m_pdfErrInfo.nmemCore() <= 0)
466  {
467  m_uncertaintyStrategy_ = IUncertainty(NullUncertaintyStrategy());
468  PDFxTMDLOG << "Error in PDFxTMD::PDFSet::InitializeUncertaintyStrategy. PDF "
469  "set must contain more than just the central value.";
470  return;
471  }
472 
473  const auto &coreType = m_pdfErrInfo.coreType();
474  if (coreType == "replicas")
475  {
476  if (m_alternativeReplicaUncertainty)
477  {
478  m_uncertaintyStrategy_ = IUncertainty(ReplicasPercentileStrategy());
479  }
480  else
481  {
482  m_uncertaintyStrategy_ = IUncertainty(ReplicasStdDevStrategy());
483  }
484  }
485  else if (coreType == "hessian")
486  {
487  m_uncertaintyStrategy_ = IUncertainty(HessianStrategy());
488  }
489  else if (coreType == "symmhessian" || coreType == "symm-hessian")
490  {
491  m_uncertaintyStrategy_ = IUncertainty(SymmHessianStrategy());
492  }
493  else if (coreType == "unknown")
494  {
495  m_uncertaintyStrategy_ = IUncertainty(NullUncertaintyStrategy());
496  }
497  else
498  {
499  throw NotSupportError(
500  "PDFxTMD::PDFSet::InitializeUncertaintyStrategy: Unsupported error type '" +
501  coreType + "' for PDF set " + m_pdfSetName + ".");
502  }
503 
504  m_setCL =
505  (coreType != "replicas") ? m_pdfSetErrorInfo.ErrorConfLevel / 100.0 : CL1SIGMA / 100.0;
506  }
507 
509  template <typename... Args>
510  std::vector<double> CalculatePDFValues(PartonFlavor flavor, Args... args) const
511  {
512  std::vector<double> pdfs;
513  pdfs.reserve(m_pdfSetStdInfo.NumMembers);
514  for (size_t i = 0; i < size(); i++)
515  {
516  if constexpr (sizeof...(args) == 3)
517  { // TMD case
518  pdfs.emplace_back(operator[](i)->tmd(flavor, args...));
519  }
520  else
521  { // Collinear case
522  pdfs.emplace_back(operator[](i)->pdf(flavor, args...));
523  }
524  }
525  return pdfs;
526  }
527 
529  double ValidateAndGetCL(double cl) const
530  {
531  const double reqCL = (cl >= 0) ? cl / 100.0 : m_setCL;
532  if (!in_range(reqCL, 0, 1) || !in_range(m_setCL, 0, 1))
533  {
534  throw InvalidInputError(
535  "Error in PDFxTMD::PDFSet::ValidateAndGetCL. Requested or PDF set "
536  "confidence level outside [0,1] range.");
537  }
538  return reqCL;
539  }
540 
542  void ApplyConfidenceLevelScaling(PDFUncertainty &rtn, double reqCL) const
543  {
544  if (m_setCL != reqCL && !m_alternativeReplicaUncertainty)
545  {
546  double qsetCL = chisquared_quantile(m_setCL, 1);
547  double qreqCL = chisquared_quantile(reqCL, 1);
548  const double scale = std::sqrt(qreqCL / qsetCL);
549 
550  rtn.scale = scale;
551  rtn.errplus *= scale;
552  rtn.errminus *= scale;
553  rtn.errsymm *= scale;
554  }
555  }
556 
558  void StoreCoreVariationErros(PDFUncertainty &resUncertainty)
559  {
560  // Store core variation uncertainties
561  resUncertainty.errplus_pdf = resUncertainty.errplus;
562  resUncertainty.errminus_pdf = resUncertainty.errminus;
563  resUncertainty.errsymm_pdf = resUncertainty.errsymm;
564  resUncertainty.errparts.push_back(
565  {resUncertainty.errplus_pdf,
566  resUncertainty.errminus_pdf});
567  }
568 
570  void CalculateParameterVariationErrors(PDFUncertainty &rtn,
571  const std::vector<double> &values) const
572  {
573  double errsq_par_plus = 0.0;
574  double errsq_par_minus = 0.0;
575  size_t index = m_pdfErrInfo.nmemCore();
576 
577  for (size_t iq = 1; iq < m_pdfErrInfo.qparts.size(); ++iq)
578  {
579  const auto &eparts = m_pdfErrInfo.qparts[iq];
580  double vmin = rtn.central;
581  double vmax = rtn.central;
582 
583  for (const auto &epart : eparts)
584  {
585  const bool symm = StartsWith(epart.first, "$");
586  for (size_t ie = 0; ie < epart.second; ++ie)
587  {
588  index++;
589  if (!symm)
590  {
591  vmin = std::min(values[index], vmin);
592  vmax = std::max(values[index], vmax);
593  }
594  else
595  {
596  const double delta = values[index] - rtn.central;
597  vmin = std::min({values[index], rtn.central - delta, vmin});
598  vmax = std::max({values[index], rtn.central - delta, vmax});
599  }
600  }
601  }
602 
603  const double eplus = vmax - rtn.central;
604  const double eminus = rtn.central - vmin;
605  rtn.errparts.push_back({eplus, eminus});
606  errsq_par_plus += SQR(eplus);
607  errsq_par_minus += SQR(eminus);
608  }
609  }
610  std::string m_pdfSetName;
611  ConfigWrapper m_pdfSetInfo;
612  YamlErrorInfo m_pdfSetErrorInfo;
613  YamlStandardTMDInfo m_pdfSetStdInfo;
614  PDFErrInfo m_pdfErrInfo;
615  std::map<unsigned, std::unique_ptr<PDF_t>> m_PDFSet_;
616  std::vector<unsigned int> m_createdPDFSetsMember;
617  IUncertainty m_uncertaintyStrategy_;
618  IQCDCoupling m_qcdCoupling;
619  double m_setCL;
620  bool m_alternativeReplicaUncertainty;
621  bool m_isValid = false;
622  std::mutex m_pdfSetMtx;
623 };
624 
625 } // namespace PDFxTMD
#define PDFxTMDLOG
Definition: Logger.h:60
#define NO_REQUESTED_CONFIDENCE_LEVEL
Definition: PartonUtils.h:26
#define SQR(x)
Definition: PartonUtils.h:25
Definition: ConfigWrapper.h:29
bool loadFromFile(const std::filesystem::path &filepath, Format format)
Factory class for creating QCD coupling objects.
Definition: Factory.h:17
Factory class for creating collinear PDF objects.
Definition: Factory.h:84
Factory class for creating TMD (Transverse Momentum Dependent) PDF objects.
Definition: Factory.h:51
Interface for Collinear Parton Distribution Functions (CPDFs).
Definition: ICPDF.h:21
double AlphaQCDMu2(double mu2) const
Definition: IQCDCoupling.h:36
Interface for Transverse Momentum Dependent (TMD) parton distribution functions.
Definition: ITMD.h:23
void Uncertainty(const std::vector< double > &values, const int numCoreErrMember, const double cl, PDFUncertainty &uncertainty) const
Definition: IUncertainty.h:41
double Correlation(const std::vector< double > &valuesA, const std::vector< double > &valuesB, const int numCoreErrMember) const
Definition: IUncertainty.h:47
Definition: Exception.h:72
null uncertainty for pdf set that do not impelement uncertainty.
Definition: NullQCDCoupling.h:10
null uncertainty for pdf set that do not impelement uncertainty.
Definition: NullUncertaintyStrategy.h:9
Manages a set of Parton Distribution Functions (PDFs), providing tools for uncertainty and correlatio...
Definition: PDFSet.h:60
PDFUncertainty Uncertainty(PartonFlavor flavor, double x, double mu2, double cl=NO_REQUESTED_CONFIDENCE_LEVEL)
Calculate the collinear PDF uncertainty and return the result.
Definition: PDFSet.h:191
size_t size() const
Get the total number of members in this PDF set.
Definition: PDFSet.h:338
PDFSet(const PDFSet &)=delete
YamlStandardTMDInfo getStdPDFInfo() const
Get the standard metadata for the PDF set.
Definition: PDFSet.h:347
PDFSet()
Default constructor.
Definition: PDFSet.h:81
typename PDFType< Tag >::type PDF_t
The specific PDF interface type (ITMD or ICPDF) determined by the Tag.
Definition: PDFSet.h:63
double alphasQ2(double q2) const
Get the strong coupling constant alpha_s at a given squared scale Q^2.
Definition: PDFSet.h:98
void InitailizePDFSetName(std::string pdfSetName)
Re-initializes the PDFSet with a new PDF set name.
Definition: PDFSet.h:328
void Uncertainty(PartonFlavor flavor, double x, double mu2, double cl, PDFUncertainty &resUncertainty)
Calculate the collinear PDF uncertainty.
Definition: PDFSet.h:156
YamlErrorInfo getPDFErrorInfo() const
Get the error metadata for the PDF set.
Definition: PDFSet.h:356
double Correlation(PartonFlavor flavorA, double xA, double kt2A, double mu2A, PartonFlavor flavorB, double xB, double kt2B, double mu2B)
Calculate the correlation for TMDs.
Definition: PDFSet.h:263
ConfigWrapper info()
Get the full configuration metadata object for the set.
Definition: PDFSet.h:365
double alphasQ(double q) const
Get the strong coupling constant alpha_s at a given scale Q.
Definition: PDFSet.h:88
PDFSet & operator=(const PDFSet &)=delete
PDFSet & operator=(PDFSet &&)=default
PDF_t * operator[](int member) const
Access a specific PDF member from the set (const version).
Definition: PDFSet.h:123
PDF_t * operator[](int member)
Access a specific PDF member from the set.
Definition: PDFSet.h:108
void Uncertainty(const std::vector< double > &values, double cl, PDFUncertainty &resUncertainty)
Calculate uncertainty from a pre-computed vector of PDF values.
Definition: PDFSet.h:206
PDFUncertainty Uncertainty(const std::vector< double > &values, double cl=NO_REQUESTED_CONFIDENCE_LEVEL)
Calculate uncertainty from a pre-computed vector of PDF values and return the result.
Definition: PDFSet.h:220
void CreateAllPDFSets()
Pre-loads all PDF members in the set.
Definition: PDFSet.h:314
double Correlation(PartonFlavor flavorA, double xA, double mu2A, PartonFlavor flavorB, double xB, double mu2B)
Calculate the correlation for collinear PDFs.
Definition: PDFSet.h:242
PDFSet(std::string pdfSetName, bool alternativeReplicaUncertainty=false)
Constructs a PDFSet for a given PDF name.
Definition: PDFSet.h:70
void CreatePDFSet(unsigned int setMember)
Explicitly creates a single PDF member.
Definition: PDFSet.h:289
PDFUncertainty Uncertainty(PartonFlavor flavor, double x, double kt2, double mu2, double cl=NO_REQUESTED_CONFIDENCE_LEVEL)
Calculate the TMD uncertainty and return the result.
Definition: PDFSet.h:173
void Uncertainty(PartonFlavor flavor, double x, double kt2, double mu2, double cl, PDFUncertainty &resUncertainty)
Calculate the TMD uncertainty.
Definition: PDFSet.h:140
double Correlation(const std::vector< double > &valuesA, const std::vector< double > &valuesB) const
Calculate correlation from two pre-computed vectors of PDF values.
Definition: PDFSet.h:277
PDFSet(PDFSet &&)=default
Definition: AllFlavorsShape.h:14
std::pair< std::optional< YamlStandardTMDInfo >, ErrorType > YamlStandardPDFInfoReader(const std::string &yamlInfoPath)
std::pair< std::optional< YamlErrorInfo >, ErrorType > YamlErrorInfoReader(const std::string &yamlInfoPath)
PartonFlavor
Definition: PartonUtils.h:58
double chisquared_quantile(double p, double ndf)
Compute quantiles of the chi-squared probability distribution function.
int in_range(double x, double low, double high)
Check if a number is in a range (closed-open) (from lhapdf)
Definition: PartonUtils.h:111
std::pair< std::optional< std::string >, ErrorType > StandardInfoFilePath(const std::string &pdfSetName)
bool StartsWith(const std::string &str, const std::string &prefix)
Definition: StringUtils.h:58
int mu2
Definition: pdfset_tutorial.py:14
int kt2
Definition: pdfset_tutorial.py:15
Definition: GenericPDF.h:30
size_t nmemCore() const
Number of core-set members.
QuadParts qparts
Error-set quadrature parts.
Definition: PDFErrInfo.h:27
std::string coreType() const
Calculated name of a quadrature part.
Definition: PDFErrInfo.h:36
static PDFErrInfo CalculateErrorInfo(const YamlErrorInfo &yamlErrInfo)
Parse extended error type syntax.
A type trait to map a tag to a specific PDF interface type.
Definition: PDFSet.h:35
Structure for storage of uncertainty info calculated over a PDF error set.
Definition: Uncertainty.h:10
double central
Variables for the central value, +ve, -ve & symmetrised errors, and a CL scalefactor.
Definition: Uncertainty.h:25
Definition: GenericPDF.h:27
Definition: YamlErrorInfo.h:9
size_t size
Definition: YamlErrorInfo.h:13
double ErrorConfLevel
Definition: YamlErrorInfo.h:10
int NumMembers
Definition: YamlStandardPDFInfo.h:18
Definition: YamlStandardPDFInfo.h:26