GridPDF.h

// -*- C++ -*-
//
// This file is part of LHAPDF
// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#pragma once
#ifndef LHAPDF_GridPDF_H
#define LHAPDF_GridPDF_H

#include "LHAPDF/PDF.h"
#include "LHAPDF/Interpolator.h"
#include "LHAPDF/Extrapolator.h"
#include "LHAPDF/KnotArray.h"

namespace LHAPDF {


  class GridPDF : public PDF {
  public:



    GridPDF() {
      _mempath = "";
      _info = PDFInfo();
      _forcePos = -1;
    }

    GridPDF(const std::string& path) {
      _loadInfo(path); // Sets _mempath
      _loadData(_mempath);
      _forcePos = -1;
    }

    GridPDF(const std::string& setname, int member) {
      _loadInfo(setname, member); // Sets _mempath
      _loadData(_mempath);
      _forcePos = -1;
    }

    GridPDF(int lhaid) {
      _loadInfo(lhaid); // Sets _mempath
      _loadData(_mempath);
      _forcePos = -1;
    }

    virtual ~GridPDF() { }



  protected:

    void _loadData(const std::string& mempath);


  public:



    void setInterpolator(Interpolator* ipol) {
      _interpolator.reset(ipol);
      _interpolator->bind(this);
    }

    template <typename INTERPOLATOR>
    void setInterpolator(INTERPOLATOR ipol) {
      setInterpolator(new INTERPOLATOR(ipol));
    }

    void setInterpolator(const std::string& ipolname) {
      setInterpolator(mkInterpolator(ipolname));
    }

    bool hasInterpolator() const {
      return _interpolator.get() != 0;
    }

    const Interpolator& interpolator() const {
      if (_interpolator.get() == 0) { // Load the default interpolator lazily
        // NB. The following is equiv to set-by-name but is explicitly implemented here for const correctness
        const string ipolname = info().get_entry("Interpolator");
        Interpolator* ipol = mkInterpolator(ipolname);
        _interpolator.reset(ipol);
        _interpolator->bind(this);
      }
      return *_interpolator;
    }



    void setExtrapolator(Extrapolator* xpol) {
      _extrapolator.reset(xpol);
      _extrapolator->bind(this);
    }

    template <typename EXTRAPOLATOR>
    void setExtrapolator(EXTRAPOLATOR xpol) {
      setExtrapolator(new EXTRAPOLATOR(xpol));
    }

    void setExtrapolator(const std::string& xpolname) {
      setExtrapolator(mkExtrapolator(xpolname));
    }

    bool hasExtrapolator() const {
      return _extrapolator.get() != 0;
    }

    const Extrapolator& extrapolator() const {
      if (_extrapolator.get() == 0) { // Load the default extrapolator lazily
        // NB. The following is equiv to set-by-name but is explicitly implemented here for const correctness
        const string xpolname = info().get_entry("Extrapolator");
        Extrapolator* xpol = mkExtrapolator(xpolname);
        _extrapolator.reset(xpol);
        _extrapolator->bind(this);
      }
      return *_extrapolator;
    }



  protected:

    double _xfxQ2(int id, double x, double q2) const;


  public:



    std::map<double, KnotArrayNF>& knotarrays() {
      return _knotarrays;
    }

    const KnotArrayNF& subgrid(double q2) const {
      assert(q2 >= 0);
      assert(!q2Knots().empty());
      map<double, KnotArrayNF>::const_iterator it = _knotarrays.upper_bound(q2);
      if (it == _knotarrays.begin())
        throw GridError("Requested Q2 " + to_str(q2) + " is lower than any available Q2 subgrid (lowest Q2 = " + to_str(q2Knots().front()) + ")");
      if (it == _knotarrays.end() && q2 > q2Knots().back())
        throw GridError("Requested Q2 " + to_str(q2) + " is higher than any available Q2 subgrid (highest Q2 = " + to_str(q2Knots().back()) + ")");
      --it; // upper_bound (and lower_bound) returns the entry *above* q2: we need to decrement by one element
      // std::cout << "Using subgrid #" << std::distance(_knotarrays.begin(), it) << std::endl;
      return it->second;
    }

    const KnotArray1F& subgrid(int id, double q2) const {
      return subgrid(q2).get_pid(id);
    }

    const vector<double>& xKnots() const {
      const KnotArrayNF& subgrid1 = _knotarrays.begin()->second;
      const KnotArray1F& grid1 = subgrid1.get_first();
      return grid1.xs();
    }

    const vector<double>& q2Knots() const {
      if (_q2knots.empty()) {
        // Get the list of Q2 knots by combining all subgrids
        for (map<double, KnotArrayNF>::const_iterator isub = _knotarrays.begin(); isub != _knotarrays.end(); ++isub) {
          const KnotArrayNF& subgrid = isub->second;
          const KnotArray1F& grid1 = subgrid.get_first();
          if (grid1.q2s().empty()) continue; //< @todo This shouldn't be possible, right? Throw instead, or ditch the check?
          for (double q2 : grid1.q2s()) {
            if (_q2knots.empty() || q2 != _q2knots.back()) _q2knots.push_back(q2);
          }
        }
      }
      return _q2knots;
    }


  public:

    bool inRangeX(double x) const {
      assert(!xKnots().empty());
      if (x < xKnots().front()) return false;
      if (x > xKnots().back()) return false;
      return true;
    }

    bool inRangeQ2(double q2) const {
      assert(!q2Knots().empty());
      if (q2 < q2Knots().front()) return false;
      if (q2 > q2Knots().back()) return false;
      return true;
    }



  private:

    std::map<double, KnotArrayNF> _knotarrays;

    // /// Caching vector of x knot values
    // mutable std::vector<double> _xknots;

    mutable std::vector<double> _q2knots;

    typedef unique_ptr<Interpolator> InterpolatorPtr;

    typedef unique_ptr<Extrapolator> ExtrapolatorPtr;

    mutable InterpolatorPtr _interpolator;

    mutable ExtrapolatorPtr _extrapolator;

  };


}
#endif