GridPDF.h

// -*- C++ -*-
//
// This file is part of LHAPDF
// Copyright (C) 2012-2022 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
      _loadData(_mempath);
      _loadPlugins();
      _forcePos = -1;
    }

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

    virtual ~GridPDF() { }


  protected:

    void _loadInterpolator();

    void _loadExtrapolator();

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

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

    void _computePolynomialCoefficients(bool logspace);


  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;

    void   _xfxQ2(double x, double q2, std::vector<double>& ret) const;


  public:


    // accerror to new data structure
    const KnotArray& knotarray() const {
      return data;
    }

    KnotArray& Data() {
      return data;
    }

    const vector<double>& xKnots() const;

    const vector<double>& q2Knots() const;

  public:

    bool inRangeX(double x) const {
      return data.inRangeX(x);
    }

    bool inRangeQ2(double q2) const {
      return data.inRangeQ2(q2);
    }

  private:
    // *new* memory object, to handle basically everything
    // name?
    KnotArray data;
        
    typedef unique_ptr<Interpolator> InterpolatorPtr;

    typedef unique_ptr<Extrapolator> ExtrapolatorPtr;

    mutable InterpolatorPtr _interpolator;

    mutable ExtrapolatorPtr _extrapolator;

  };


}
#endif