PDF.h

// -*- C++ -*-
//
// This file is part of LHAPDF
// Copyright (C) 2012-2022 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("PdfDesc", "");
    }

    int dataversion() const {
      return info().get_entry_as<int>("DataVersion", -1);
    }

    std::string type() const {
      return to_lower(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