Histogram.h

Go to the documentation of this file.

/*
 *  CepGen: a central exclusive processes event generator
 *  Copyright (C) 2021-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_Histogram_h
#define CepGen_Utils_Histogram_h
 
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_histogram2d.h>
 
#include <array>
#include <functional>
#include <memory>
#include <set>
 
#include "CepGen/Utils/Drawable.h"
 
namespace cepgen {
  class ParametersList;
  namespace utils {
    class RandomGenerator;
    class Histogram {
    public:
      Histogram() = default;
      virtual ~Histogram() = default;
 
      enum BinMode { low = 0, high, both };
 
      virtual void clear() = 0;
      virtual void scale(double) = 0;
      virtual double integral(bool include_out_of_range = false) const = 0;
      virtual double minimum() const = 0;
      virtual double maximum() const = 0;
      void normalise(double integ = 1.);
 
    protected:
      std::set<double> extractBins(BinMode mode,
                                   size_t num_bins,
                                   const std::function<Limits(size_t)>& bins_extractor) const;
    };
 
    class Hist1D : public Histogram, public Drawable {
    public:
      explicit Hist1D(const ParametersList&);
      explicit Hist1D(size_t num_bins_x, const Limits&, const std::string& name = "", const std::string& title = "");
      explicit Hist1D(const std::vector<double>& bins, const std::string& name = "", const std::string& title = "");
      Hist1D(const Hist1D&);  
 
      void clear() override;
      void fill(double x, double weight = 1.);
      void add(Hist1D, double scaling = 1.);
      void scale(double) override;
      double sample(RandomGenerator&) const;
 
      double chi2test(const Hist1D&, size_t& ndf) const;
 
      std::vector<Value> values() const;
      Value value(size_t bin) const;
      void setValue(size_t bin, Value value);
 
      axis_t axis() const;
      size_t nbins() const;
      Limits range() const;
      Limits binRange(size_t bin) const;
      std::vector<double> bins(BinMode) const;
      size_t bin(double x) const;
 
      double mean() const;
      double rms() const;
      double minimum() const override;
      double maximum() const override;
      double integral(bool = false) const override;
      size_t underflow() const { return underflow_; }
      size_t overflow() const { return overflow_; }
 
      bool isHist1D() const override final { return true; }
 
    private:
      void buildFromBins(const std::vector<double>&);
      void buildFromRange(size_t, const Limits&);
 
      struct gsl_histogram_deleter {
        void operator()(gsl_histogram* h) { gsl_histogram_free(h); }
      };
      typedef std::unique_ptr<gsl_histogram, gsl_histogram_deleter> gsl_histogram_ptr;
      gsl_histogram_ptr hist_, hist_w2_;
      size_t underflow_{0ull}, overflow_{0ull};
      struct gsl_histogram_pdf_deleter {
        void operator()(gsl_histogram_pdf* h) { gsl_histogram_pdf_free(h); }
      };
      typedef std::unique_ptr<gsl_histogram_pdf, gsl_histogram_pdf_deleter> gsl_histogram_pdf_ptr;
      mutable gsl_histogram_pdf_ptr pdf_;
    };
 
    class Hist2D : public Histogram, public Drawable {
    public:
      explicit Hist2D(const ParametersList&);
      explicit Hist2D(size_t num_bins_x,
                      const Limits& xlim,
                      size_t num_bins_y,
                      const Limits& ylim,
                      const std::string& name = "",
                      const std::string& title = "");
      explicit Hist2D(const std::vector<double>& xbins,
                      const std::vector<double>& ybins,
                      const std::string& name = "",
                      const std::string& title = "");
      Hist2D(const Hist2D&);  
 
      void clear() override;
      void fill(double x, double y, double weight = 1.);
      inline void fill(const std::pair<double, double>& xy, double weight = 1.) { fill(xy.first, xy.second, weight); }
      void add(Hist2D, double scaling = 1.);
      void scale(double) override;
      std::pair<double, double> sample(RandomGenerator&) const;
 
      Value value(size_t bin_x, size_t bin_y) const;
      void setValue(size_t bin_x, size_t bin_y, Value value);
 
      size_t nbinsX() const;
      Limits rangeX() const;
      Limits binRangeX(size_t bin) const;
      std::vector<double> binsX(BinMode) const;
      size_t nbinsY() const;
      Limits rangeY() const;
      Limits binRangeY(size_t bin) const;
      std::vector<double> binsY(BinMode) const;
      std::pair<size_t, size_t> bin(double x, double y) const;
 
      double meanX() const;
      double rmsX() const;
      double meanY() const;
      double rmsY() const;
      double minimum() const override;
      double maximum() const override;
      double integral(bool = false) const override;
 
      struct contents_t : public std::array<size_t, 8> {
        contents_t() { std::fill(begin(), end(), 0ull); }
        size_t total() const;
        friend std::ostream& operator<<(std::ostream&, const contents_t&);
        friend contents_t operator*(double, const contents_t&);
        contents_t& operator+=(const contents_t&);
        enum {
          LT_GT,
          IN_GT,
          GT_GT,
          LT_IN,
          /* INSIDE */ GT_IN,
          LT_LT,
          IN_LT,
          GT_LT,
          num_content
        };
      };
      const contents_t& outOfRange() const { return out_of_range_values_; }
 
      bool isHist2D() const override final { return true; }
 
    private:
      void buildFromBins(const std::vector<double>&, const std::vector<double>&);
      void buildFromRange(size_t, const Limits&, size_t, const Limits&);
 
      struct gsl_histogram2d_deleter {
        void operator()(gsl_histogram2d* h) { gsl_histogram2d_free(h); }
      };
      typedef std::unique_ptr<gsl_histogram2d, gsl_histogram2d_deleter> gsl_histogram2d_ptr;
      gsl_histogram2d_ptr hist_, hist_w2_;
      contents_t out_of_range_values_;
      struct gsl_histogram2d_pdf_deleter {
        void operator()(gsl_histogram2d_pdf* h) { gsl_histogram2d_pdf_free(h); }
      };
      typedef std::unique_ptr<gsl_histogram2d_pdf, gsl_histogram2d_pdf_deleter> gsl_histogram2d_pdf_ptr;
      mutable gsl_histogram2d_pdf_ptr pdf_;
    };
  }  // namespace utils
}  // namespace cepgen
 
#endif