HistoGroup.hh

// -*- C++ -*-
#ifndef RIVET_HISTOGROUP_HH
#define RIVET_HISTOGROUP_HH
 
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/Tools/RivetYODA.hh"
 
namespace Rivet {
 
  namespace {
 
    template <size_t FillDim, typename T, typename AxisT>
    typename YODA::FillableStorage<FillDim, T, AxisT>::FillAdapterT
    groupAdapter = [](auto& /* ptr */, auto&& /* coords */, double /* weight */, double /* fraction */) { };
 
  }
 
  template <typename GroupAxisT, typename... AxisT>
  class HistoGroup : public YODA::FillableStorage<sizeof...(AxisT)+1, BinnedHistoPtr<AxisT...>, GroupAxisT>,
                     public YODA::Fillable {
  public:
 
    using BaseT = YODA::FillableStorage<sizeof...(AxisT)+1, BinnedHistoPtr<AxisT...>, GroupAxisT>;
    using BinContentT = BinnedHistoPtr<AxisT...>;
    using FillDim = std::integral_constant<size_t, sizeof...(AxisT)+1>;
 
    HistoGroup() : BaseT(groupAdapter<FillDim{}, BinContentT, GroupAxisT>) { }
 
    HistoGroup(const std::vector<GroupAxisT>& edges)
         : BaseT(YODA::Axis<GroupAxisT>(edges), groupAdapter<FillDim{}, BinContentT, GroupAxisT>) { }
 
    HistoGroup(std::initializer_list<GroupAxisT>&& edges)
         : BaseT(YODA::Axis<GroupAxisT>(std::move(edges)), groupAdapter<FillDim{}, BinContentT, GroupAxisT>) { }
 
    int fill(GroupAxisT tmpCoord, AxisT... coords, const double weight = 1.0, const double fraction = 1.0) {
      auto& bin = BaseT::binAt(tmpCoord);
      if (!bin.raw())  return -1; // nullptr if bin not booked
      return bin->fill({coords...}, weight, fraction);
    }
 
 
    void reset() noexcept { BaseT::reset(); }
 
 
 
    size_t fillDim() const noexcept { return FillDim{}; }
 
    size_t dim() const noexcept { return fillDim()+1; }
 
    double numEntries(const bool includeOverflows=true) const noexcept {
      double n = 0;
      for (const auto& b : BaseT::bins(includeOverflows)) {
        if (!b.get())  continue;
        n += b->numEntries(includeOverflows);
      }
      return n;
    }
 
    double effNumEntries(const bool includeOverflows=true) const noexcept {
      double n = 0;
      for (const auto& b : BaseT::bins(includeOverflows)) {
        if (!b.get())  continue;
        n += b->effNumEntries(includeOverflows);
      }
      return n;
    }
 
    double sumW(const bool includeOverflows=true) const noexcept {
      double sumw = 0;
      for (const auto& b : BaseT::bins(includeOverflows)) {
        if (!b.get())  continue;
        sumw += b->sumW(includeOverflows);
      }
      return sumw;
    }
 
    double sumW2(const bool includeOverflows=true) const noexcept {
      double sumw2 = 0;
      for (const auto& b : BaseT::bins(includeOverflows)) {
        if (!b.get())  continue;
        sumw2 += b->sumW2(includeOverflows);
      }
      return sumw2;
    }
 
    double integral(const bool includeOverflows=true) const noexcept {
      return sumW(includeOverflows);
    }
 
    double integralError(const bool includeOverflows=true) const noexcept {
      return sqrt(sumW2(includeOverflows));
    }
 
 
 
    void divByGroupWidth() const noexcept {
      for (auto& bin : BaseT::bins(true, true)) {
        if (!bin.get())  continue;
        const double bw = bin.dVol();
        if (bw)  bin->scaleW(1.0/bw);
      }
    }
 
    void scaleW(const double scalefactor) noexcept {
      for (auto& bin : BaseT::bins(true, true)) {
        if (!bin.get())  continue;
        bin->scaleW(scalefactor);
      }
    }
 
    void scale(const size_t i, const double scalefactor) noexcept {
      for (auto& bin : BaseT::bins(true, true)) {
        if (!bin.get())  continue;
        bin->scale(i, scalefactor);
      }
    }
 
    void normalize(const double normto=1.0, const bool includeOverflows=true) {
      for (auto& bin : BaseT::bins(true, true)) {
        if (!bin.get())  continue;
        if (bin->integral(includeOverflows) != 0.)  bin->normalize(normto, includeOverflows);
      }
    }
 
    void normalizeGroup(const double normto=1.0, const bool includeOverflows=true) {
      const double oldintegral = integral(includeOverflows);
      if (oldintegral == 0) throw WeightError("Attempted to normalize a histogram group with null area");
      scaleW(normto / oldintegral);
    }
 
 
  };
 
 
  template <typename GroupAxisT, typename... AxisT>
  using HistoGroupPtr = std::shared_ptr<HistoGroup<GroupAxisT, AxisT...>>;
  using Histo1DGroupPtr = HistoGroupPtr<double,double>;
  using Histo2DGroupPtr = HistoGroupPtr<double,double,double>;
 
 
}
 
#endif