HistoGroup.hh

// -*- C++ -*-
#ifndef RIVET_HISTOGROUP_HH
#define RIVET_HISTOGROUP_HH
 
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/Tools/Logging.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;
    }

    vector<double> sumWGroup(const bool includeOverflows = true) const noexcept {
      vector<double> rtn;
      rtn.reserve(BaseT::numBins(true));
      for (const auto& b : BaseT::bins(true)) {
        if (!b.get()) {
          rtn.push_back(0.0);
          continue;
        }
        rtn.push_back(b->sumW(includeOverflows));
      }
      return rtn;
    }

    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) {
        MSG_DEBUG("Attempted to normalize a histogram group with null area; skipping.");
        return;
      }
      scaleW(normto / oldintegral);
    }

 
  protected:

    Log& getLog() const {
      string logname = "Rivet.HistoGroup";
      return Log::getLog(logname);
    }
  };

 
  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