GridHandler.h

Go to the documentation of this file.

/*
 *  CepGen: a central exclusive processes event generator
 *  Copyright (C) 2018-2024  Laurent Forthomme
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
#ifndef CepGen_Utils_GridHandler_h
#define CepGen_Utils_GridHandler_h
 
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_version.h>
#if defined(GSL_MAJOR_VERSION) && (GSL_MAJOR_VERSION > 2 || (GSL_MAJOR_VERSION == 2 && GSL_MINOR_VERSION >= 1))
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>
#define GSL_VERSION_ABOVE_2_1 1
#endif
 
#include <array>
#include <map>
#include <memory>
 
#include "CepGen/Utils/Limits.h"
 
namespace cepgen {
  enum struct GridType { linear, logarithmic, square };
  template <size_t D, size_t N = 1>
  class GridHandler {
  public:
    typedef std::vector<double> coord_t;     
    typedef std::array<double, N> values_t;  
 
  public:
    explicit GridHandler(const GridType& grid_type);  
    virtual ~GridHandler() = default;
 
    values_t eval(coord_t in_coords) const;  
 
    void insert(coord_t coord, values_t value);                                
    inline std::map<coord_t, values_t> values() const { return values_raw_; }  
 
    void initialise();                         
    std::array<Limits, D> boundaries() const;  
    std::array<double, D> min() const;         
    std::array<double, D> max() const;         
 
  protected:
    const GridType grid_type_;                
    std::map<coord_t, values_t> values_raw_;  
    std::vector<std::unique_ptr<gsl_interp_accel, void (*)(gsl_interp_accel*)> > accel_;
    std::vector<std::unique_ptr<gsl_spline, void (*)(gsl_spline*)> > splines_1d_;  
#ifdef GSL_VERSION_ABOVE_2_1
    std::vector<std::unique_ptr<gsl_spline2d, void (*)(gsl_spline2d*)> > splines_2d_;
#endif
    std::array<coord_t, D> coords_;                    
    std::array<std::unique_ptr<double[]>, N> values_;  
 
  private:
    void findIndices(const coord_t& coord, coord_t& min, coord_t& max) const;  
    struct gridpoint_t : values_t {
      gridpoint_t(const values_t& arr) : values_t(arr) {}
      gridpoint_t operator*(double c) const;
      gridpoint_t operator+(const gridpoint_t& rhs) const;
    };
    bool init_{false};  
  };
}  // namespace cepgen
 
#endif