PDF.h

// -*- C++ -*-
//
// This file is part of LHAPDF
// Copyright (C) 2012-2023 The LHAPDF collaboration (see AUTHORS for details)
//
#pragma once
#ifndef LHAPDF_PDF_H
#define LHAPDF_PDF_H
 
#include "LHAPDF/PDFInfo.h"
#include "LHAPDF/PDFIndex.h"
#include "LHAPDF/Factories.h"
#include "LHAPDF/AlphaS.h"
#include "LHAPDF/Utils.h"
#include "LHAPDF/Paths.h"
#include "LHAPDF/Exceptions.h"
#include "LHAPDF/Version.h"
#include "LHAPDF/Config.h"
 
namespace LHAPDF {
 
 
  namespace PIDs { //< for scoping to avoid conflicts, e.g. PID::GLUON rather than just GLUON
    enum PIDCode {
      ATOP = -6, ABOTTOM = -5, ACHARM = -4, ASTRANGE = -3, AUP = -2, ADOWN = -1,
      GLUON = 0, // equivalent to 21
      DOWN = 1, UP = 2, STRANGE = 3, CHARM = 4, BOTTOM = 5, TOP = 6
    };
  }
 
 
 
  class PDF {
  protected: //< These constructors should only be called by subclasses
 
    typedef unique_ptr<AlphaS> AlphaSPtr;
 
    PDF() : _forcePos(0) { }
 
 
  public:
 
    virtual ~PDF() { }
 
 
  protected:
 
 
 
    void _loadInfo(const std::string& mempath);
 
    void _loadInfo(const std::string& setname, int member) {
      const string searchpath = findpdfmempath(setname, member);
      if (searchpath.empty())
        throw UserError("Can't find a valid PDF " + setname + "/" + to_str(member));
      _loadInfo(searchpath);
    }
 
    void _loadInfo(int lhaid) {
      const pair<string,int> setname_memid = lookupPDF(lhaid);
      if (setname_memid.second == -1)
        throw IndexError("Can't find a PDF with LHAPDF ID = " + to_str(lhaid));
      _loadInfo(setname_memid.first, setname_memid.second);
    }
 
 
 
  public:
 
 
    double xfxQ2(int id, double x, double q2) const;
 
 
    double xfxQ(int id, double x, double q) const {
      return xfxQ2(id, x, q*q);
    }
 
 
    void xfxQ2(double x, double q2, std::map<int, double>& rtn) const;
 
 
    void xfxQ(double x, double q, std::map<int, double>& rtn) const {
      xfxQ2(x, q*q, rtn);
    }
 
 
    void xfxQ2(double x, double q2, std::vector<double>& rtn) const;
 
 
    void xfxQ(double x, double q, std::vector<double>& rtn) const {
      xfxQ2(x, q*q, rtn);
    }
 
 
    std::map<int, double> xfxQ2(double x, double q2) const;
 
    std::map<int, double> xfxQ(double x, double q) const {
      return xfxQ2(x, q*q);
    }
 
 
  protected:
 
    virtual double _xfxQ2(int id, double x, double q2) const = 0;
 
    virtual void _xfxQ2(double x, double q2, std::vector<double>& ret) const = 0;
 
 
 
  public:
 
 
    virtual double xMin() {
      if (info().has_key("XMin"))
        return info().get_entry_as<double>("XMin");
      return numeric_limits<double>::epsilon();
    }
 
    virtual double xMax() {
      if (info().has_key("XMax"))
        return info().get_entry_as<double>("XMax");
      return 1.0;
    }
 
    virtual double qMin() {
      return info().get_entry_as<double>("QMin", 0);
    }
 
    virtual double qMax() {
      return info().get_entry_as<double>("QMax", numeric_limits<double>::max());
    }
 
    virtual double q2Min() {
      return sqr(this->qMin());
    }
 
    virtual double q2Max() {
      // Explicitly re-access this from the info, to avoid an overflow from squaring double_max
      return (info().has_key("QMax")) ? sqr(info().get_entry_as<double>("QMax")) : numeric_limits<double>::max();
    }
 
 
    int forcePositive() const {
      if (_forcePos < 0) //< Caching
        _forcePos = info().get_entry_as<unsigned int>("ForcePositive", 0);
      return _forcePos;
    }
    void setForcePositive(int mode) {
      _forcePos = mode;
    }
 
 
    bool inPhysicalRangeX(double x) const {
      return x >= 0.0 && x <= 1.0;
    }
 
    bool inPhysicalRangeQ2(double q2) const {
      return q2 >= 0.0;
    }
 
    bool inPhysicalRangeQ(double q) const {
      return inPhysicalRangeQ2(q*q);
    }
 
    bool inPhysicalRangeXQ2(double x, double q2) const {
      return inPhysicalRangeX(x) && inPhysicalRangeQ2(q2);
    }
 
    bool inPhysicalRangeXQ(double x, double q) const {
      return inPhysicalRangeX(x) && inPhysicalRangeQ(q);
    }
 
    virtual bool inRangeQ(double q) const {
      return inRangeQ2(q*q);
    }
 
    virtual bool inRangeQ2(double q2) const = 0;
 
    virtual bool inRangeX(double x) const = 0;
 
    virtual bool inRangeXQ(double x, double q) const {
      return inRangeX(x) && inRangeQ(q);
    }
 
    bool inRangeXQ2(double x, double q2) const {
      return inRangeX(x) && inRangeQ2(q2);
    }
 
 
 
 
    PDFInfo& info() { return _info; }
 
    const PDFInfo& info() const { return _info; }
 
    PDFSet& set() const {
      return getPDFSet(_setname());
    }
 
 
 
 
    int memberID() const {
      const string memname = file_stem(_mempath);
      assert(memname.length() > 5); // There must be more to the filename stem than just the _nnnn suffix
      const int memid = lexical_cast<int>(memname.substr(memname.length()-4)); //< Last 4 chars should be the member number
      return memid;
    }
 
    int lhapdfID() const;
 
    std::string description() const {
      return info().get_entry("MemDesc", info().get_entry("PdfDesc", ""));
    }
 
    int dataversion() const {
      return info().get_entry_as<int>("DataVersion", -1);
    }
 
    std::string type() const {
      return to_lower(info().get_entry("MemType", info().get_entry("PdfType")));
    }
 
 
 
    void print(std::ostream& os=std::cout, int verbosity=1) const;
 
 
 
    virtual const std::vector<int>& flavors() const {
      if (_flavors.empty()) {
        _flavors = info().get_entry_as< vector<int> >("Flavors");
        sort(_flavors.begin(), _flavors.end());
      }
      return _flavors;
    }
 
    void setFlavors(std::vector<int> const& flavors) {
      _flavors = flavors;
      sort(_flavors.begin(), _flavors.end());
    }
 
    bool hasFlavor(int id) const;
 
    int orderQCD() const {
      return info().get_entry_as<int>("OrderQCD");
    }
    int qcdOrder() const { return orderQCD(); }
 
    double quarkMass(int id) const;
 
    double quarkThreshold(int id) const;
 
 
 
 
    void setAlphaS(AlphaS* alphas) {
      _alphas.reset(alphas);
    }
 
    void setAlphaS(AlphaSPtr alphas) {
      _alphas = std::move(alphas);
    }
 
    bool hasAlphaS() const {
      return bool(_alphas);
    }
 
    AlphaS& alphaS() {
      return *_alphas;
    }
 
    const AlphaS& alphaS() const {
      return *_alphas;
    }
 
    double alphasQ(double q) const {
      return alphasQ2(q*q);
    }
 
    double alphasQ2(double q2) const {
      if (!hasAlphaS()) throw Exception("No AlphaS pointer has been set");
      return _alphas->alphasQ2(q2);
    }
 
 
 
  protected:
 
    void _loadAlphaS() {
      _alphas.reset( mkAlphaS(info()) );
    }
 
    std::string _setname() const {
      return basename(dirname(_mempath));
    }
 
    std::string _mempath;
 
    PDFInfo _info;
 
    mutable vector<int> _flavors;
 
    mutable AlphaSPtr _alphas;
 
    mutable int _forcePos;
 
  };
 
 
}
#endif