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
      _loadPlugins();
      _loadData(_mempath);
      _forcePos = -1;
    }

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

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

    virtual ~GridPDF() { }



  protected:

    void _loadInterpolator();

    void _loadExtrapolator();

    void _loadPlugins() {
      _loadAlphaS();
      _loadInterpolator();
      _loadExtrapolator();
    }

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


  public:



    void setInterpolator(Interpolator* ipol);

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

    void setInterpolator(const std::string& ipolname);

    bool hasInterpolator() const { return bool(_interpolator); }

    const Interpolator& interpolator() const;


    void setExtrapolator(Extrapolator* xpol);

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

    void setExtrapolator(const std::string& xpolname);

    bool hasExtrapolator() const { return bool(_extrapolator); }

    const Extrapolator& extrapolator() const;



  protected:

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


  public:



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

    const KnotArrayNF& subgrid(double q2) const;

    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;

  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