BinnedDbn.h

Go to the documentation of this file.

// -*- C++ -*-
//
// This file is part of YODA -- Yet more Objects for Data Analysis
// Copyright (C) 2008-2023 The YODA collaboration (see AUTHORS for details)
//
#ifndef YODA_BinnedDbn_h
#define YODA_BinnedDbn_h
 
#include "YODA/AnalysisObject.h"
#include "YODA/Fillable.h"
#include "YODA/FillableStorage.h"
#include "YODA/Dbn.h"
#include "YODA/BinnedEstimate.h"
#include "YODA/Scatter.h"
 
#include <memory>
#include <utility>
#include <iostream>
#include <iomanip>
 
namespace YODA {
 
  /*
  */
 
 
  template <size_t DbnN, typename... AxisT>
  class DbnStorage;
 
  template <size_t DbnN, typename... AxisT>
  class BinnedDbn : public DbnStorage<DbnN, AxisT...> {
  public:
    using HistoT = BinnedDbn<DbnN, AxisT...>;
    using BaseT = DbnStorage<DbnN, AxisT...>;
    using FillType = typename BaseT::FillType;
    using BinType = typename BaseT::BinT;
    using Ptr = std::shared_ptr<HistoT>;
 
    using BaseT::BaseT;
 
    BinnedDbn() = default;
    BinnedDbn(const HistoT&) = default;
    BinnedDbn(HistoT&&) = default;
    BinnedDbn& operator =(const HistoT&) = default;
    BinnedDbn& operator =(HistoT&&) = default;
    using AnalysisObject::operator =;
 
    BinnedDbn(const BaseT& other) : BaseT(other) {}
    //
    BinnedDbn(const HistoT& other, const std::string& path) : BaseT(other, path) {}
 
    BinnedDbn(BaseT&& other) : BaseT(std::move(other)) {}
    //
    BinnedDbn(HistoT&& other, const std::string& path) : BaseT(std::move(other), path) {}
 
    HistoT clone() const noexcept {
      return HistoT(*this);
    }
 
    HistoT* newclone() const noexcept {
      return new HistoT(*this);
    }
 
  };
 
 
 
  template <typename... AxisTypes>
  using BinnedHisto = BinnedDbn<sizeof...(AxisTypes), AxisTypes...>;
 
  template <typename... AxisTypes>
  using BinnedProfile = BinnedDbn<sizeof...(AxisTypes)+1, AxisTypes...>;
 
 
 
  template <size_t DbnN, typename... AxisT>
  class DbnStorage : public FillableStorage<DbnN, Dbn<DbnN>, AxisT...>,
                     public AnalysisObject, public Fillable {
  public:
 
    using BaseT = FillableStorage<DbnN, Dbn<DbnN>, AxisT...>;
    using BinningT = typename BaseT::BinningT;
    using BinT = typename BaseT::BinT;
    using BinType = typename BaseT::BinT;
    using FillType = typename BaseT::FillType;
    using AnalysisObject::operator =;
 
 
    DbnStorage() : BaseT(), AnalysisObject(mkTypeString<DbnN, AxisT...>(), "") { }
 
    DbnStorage(std::vector<AxisT>&&... binsEdges,
               const std::string& path = "", const std::string& title = "")
         : BaseT(Axis<AxisT>(std::move(binsEdges))...),
           AnalysisObject(mkTypeString<DbnN, AxisT...>(), path, title) { }
 
    DbnStorage(const std::vector<AxisT>&... binsEdges,
               const std::string& path = "", const std::string& title = "")
         : BaseT(Axis<AxisT>(binsEdges)...),
           AnalysisObject(mkTypeString<DbnN, AxisT...>(), path, title) { }
 
    DbnStorage(std::initializer_list<AxisT>&&... binsEdges,
               const std::string& path = "", const std::string& title = "")
         : BaseT(Axis<AxisT>(std::move(binsEdges))...),
           AnalysisObject(mkTypeString<DbnN, AxisT...>(), path, title) { }
 
    template <typename EdgeT = double, typename = enable_if_all_CAxisT<EdgeT, AxisT...>>
    DbnStorage(const std::vector<size_t>& nBins, const std::vector<std::pair<EdgeT, EdgeT>>& limitsLowUp,
               const std::string& path = "", const std::string& title = "")
         : BaseT( _mkBinning(nBins, limitsLowUp, std::make_index_sequence<sizeof...(AxisT)>{}) ),
           AnalysisObject(mkTypeString<DbnN, AxisT...>(), path, title) { }
 
    DbnStorage(const Axis<AxisT>&... axes, const std::string& path = "", const std::string& title = "")
         : BaseT(axes...), AnalysisObject(mkTypeString<DbnN, AxisT...>(), path, title) { }
 
    DbnStorage(Axis<AxisT>&&... axes, const std::string& path = "", const std::string& title = "")
         : BaseT(std::move(axes)...), AnalysisObject(mkTypeString<DbnN, AxisT...>(), path, title) { }
 
    DbnStorage(const BinningT& binning, const std::string& path = "", const std::string& title = "")
         : BaseT(binning), AnalysisObject(mkTypeString<DbnN, AxisT...>(), path, title) { }
 
    DbnStorage(BinningT&& binning, const std::string& path = "", const std::string& title = "")
         : BaseT(std::move(binning)), AnalysisObject(mkTypeString<DbnN, AxisT...>(), path, title) { }
 
    template <typename EdgeT = double, typename = enable_if_all_CAxisT<EdgeT, AxisT...>>
    DbnStorage(const ScatterND<sizeof...(AxisT)+1>& s, const std::string& path = "", const std::string& title = "")
         : BaseT(_mkBinning(s, std::make_index_sequence<sizeof...(AxisT)>{})),
           AnalysisObject(mkTypeString<DbnN, AxisT...>(), path, title) { }
 
    DbnStorage(const DbnStorage& other, const std::string& path = "") : BaseT(other),
         AnalysisObject(mkTypeString<DbnN, AxisT...>(), path!=""? path : other.path(), other, other.title()) { }
 
    DbnStorage(DbnStorage&& other, const std::string& path = "") : BaseT(std::move(other)),
         AnalysisObject(mkTypeString<DbnN, AxisT...>(), path!=""? path : other.path(), other, other.title()) {  }
 
    DbnStorage clone() const noexcept {
      return DbnStorage(*this);
    }
 
    DbnStorage* newclone() const noexcept {
      return new DbnStorage(*this);
    }
 
 
 
 
    virtual int fill(FillType&& coords, const double weight = 1.0, const double fraction = 1.0) {
      return BaseT::fill(std::move(coords), std::make_index_sequence<sizeof...(AxisT)>{}, weight, fraction);
    }
 
    void scaleW(const double scalefactor) noexcept {
      setAnnotation("ScaledBy", annotation<double>("ScaledBy", 1.0) * scalefactor);
      for (auto& bin : BaseT::_bins) {
        bin.scaleW(scalefactor);
      }
    }
 
    void scale(const size_t i, const double scalefactor) noexcept {
      setAnnotation("ScaledBy", annotation<double>("ScaledBy", 1.0) * scalefactor);
      for (auto& bin : BaseT::_bins) {
        bin.scale(i, scalefactor);
      }
    }
 
 
    void normalize(const double normto=1.0, const bool includeOverflows=true) {
      const double oldintegral = integral(includeOverflows);
      if (oldintegral == 0) throw WeightError("Attempted to normalize a histogram with null area");
      scaleW(normto / oldintegral);
    }
 
 
    template <size_t axisN>
    void rebinBy(unsigned int n, size_t begin=1, size_t end=UINT_MAX) {
      if (n < 1)  throw UserError("Rebinning requested in groups of 0!");
      if (!begin) throw UserError("Visible bins start with index 1!");
      if (end > BaseT::numBinsAt(axisN)+1)  end = BaseT::numBinsAt(axisN) + 1;
      for (size_t m = begin; m < end; ++m) {
        if (m > BaseT::numBinsAt(axisN)+1) break; // nothing to be done
        const size_t myend = (m+n-1 < BaseT::numBinsAt(axisN)+1) ? m+n-1 : BaseT::numBinsAt(axisN);
        if (myend > m) {
          BaseT::template mergeBins<axisN>({m, myend});
          end -= myend-m; //< reduce upper index by the number of removed bins
        }
      }
    }
 
    template <size_t axisN>
    void rebin(unsigned int n, size_t begin=1, size_t end=UINT_MAX) {
      rebinBy<axisN>(n, begin, end);
    }
 
    template <size_t axisN>
    void rebinTo(const std::vector<typename BinningT::template getAxisT<axisN>::EdgeT>& newedges) {
      if (newedges.size() < 2)
        throw UserError("Requested rebinning to an edge list which defines no bins");
      using thisAxisT = typename BinningT::template getAxisT<axisN>;
      using thisEdgeT = typename thisAxisT::EdgeT;
      // get list of shared edges
      thisAxisT& oldAxis = BaseT::_binning.template axis<axisN>();
      const thisAxisT newAxis(newedges);
      const std::vector<thisEdgeT> eshared = oldAxis.sharedEdges(newAxis);
      if (eshared.size() != newAxis.edges().size())
        throw BinningError("Requested rebinning to incompatible edges");
      // loop over new lower bin edges (= first bin index of merge range)
      for (size_t begin = 0; begin < eshared.size() - 1; ++begin) {
        // find index of upper edge along old axis
        // (subtracting 1 gives index of last bin to be merged)
        size_t end = oldAxis.index(eshared[begin+1]) - 1;
        // if the current edge is the last visible edge before the overflow
        // merge the remaining bins into the overflow
        if (begin == newAxis.numBins(true)-1)  end = oldAxis.numBins(true)-1;
        // merge this range
        if (end > begin)  BaseT::template mergeBins<axisN>({begin, end});
        if (eshared.size() == oldAxis.edges().size())  break; // we're done
      }
    }
 
    template <size_t axisN>
    void rebin(const std::vector<typename BinningT::template getAxisT<axisN>::EdgeT>& newedges) {
      rebinTo<axisN>(newedges);
    }
 
    DbnStorage& operator = (const DbnStorage& dbn) noexcept {
      if (this != &dbn) {
        AnalysisObject::operator = (dbn);
        BaseT::operator = (dbn);
      }
      return *this;
    }
 
    DbnStorage& operator = (DbnStorage&& dbn) noexcept {
      if (this != &dbn) {
        AnalysisObject::operator = (dbn);
        BaseT::operator = (std::move(dbn));
      }
      return *this;
    }
 
 
    DbnStorage& operator += (const DbnStorage& dbn) {
      if (*this != dbn)
        throw BinningError("Arithmetic operation requires compatible binning!");
      if (AO::hasAnnotation("ScaledBy")) AO::rmAnnotation("ScaledBy");
      for (size_t i = 0; i < BaseT::numBins(true, true); ++i) {
        BaseT::bin(i) += dbn.bin(i);
      }
      BaseT::maskBins(dbn.maskedBins(), true);
      return *this;
    }
    //
    DbnStorage& operator += (DbnStorage&& dbn) {
      if (*this != dbn)
        throw BinningError("Arithmetic operation requires compatible binning!");
      if (AO::hasAnnotation("ScaledBy")) AO::rmAnnotation("ScaledBy");
      for (size_t i = 0; i < BaseT::numBins(true, true); ++i) {
        BaseT::bin(i) += std::move(dbn.bin(i));
      }
      BaseT::maskBins(dbn.maskedBins(), true);
      return *this;
    }
 
 
    DbnStorage& operator -= (const DbnStorage& dbn) {
      if (*this != dbn)
        throw BinningError("Arithmetic operation requires compatible binning!");
      if (AO::hasAnnotation("ScaledBy")) AO::rmAnnotation("ScaledBy");
      for (size_t i = 0; i < BaseT::numBins(true, true); ++i) {
        BaseT::bin(i) -= dbn.bin(i);
      }
      BaseT::maskBins(dbn.maskedBins(), true);
      return *this;
    }
    //
    DbnStorage& operator -= (DbnStorage&& dbn) {
      if (*this != dbn)
        throw BinningError("Arithmetic operation requires compatible binning!");
      if (AO::hasAnnotation("ScaledBy")) AO::rmAnnotation("ScaledBy");
      for (size_t i = 0; i < BaseT::numBins(true, true); ++i) {
        BaseT::bin(i) -= std::move(dbn.bin(i));
      }
      BaseT::maskBins(dbn.maskedBins(), true);
      return *this;
    }
 
 
 
    void reset() noexcept { BaseT::reset(); }
 
 
 
 
    size_t fillDim() const noexcept { return BaseT::fillDim(); }
 
    size_t dim() const noexcept { return sizeof...(AxisT) + 1; }
 
    std::string _config() const noexcept { return mkAxisConfig<AxisT...>(); }
 
    template <size_t I, typename E = typename BinningT::template getEdgeT<I>>
    std::vector<E> edges(const bool includeOverflows = false) const noexcept {
      return BaseT::_binning.template edges<I>(includeOverflows);
    }
 
    template <size_t I, typename E = typename BinningT::template getEdgeT<I>>
    std::enable_if_t<std::is_floating_point<E>::value, std::vector<E>>
    widths(const bool includeOverflows = false) const noexcept {
      return BaseT::_binning.template widths<I>(includeOverflows);
    }
 
    template <size_t I, typename E = typename BinningT::template getEdgeT<I>>
    std::enable_if_t<std::is_floating_point<E>::value, E> min() const noexcept {
      return BaseT::_binning.template min<I>();
    }
 
    template <size_t I, typename E = typename BinningT::template getEdgeT<I>>
    std::enable_if_t<std::is_floating_point<E>::value, E> max() const noexcept {
      return BaseT::_binning.template max<I>();
    }
 
 
 
 
    double integral(const bool includeOverflows=true) const noexcept {
      return sumW(includeOverflows);
    }
 
    double integralError(const bool includeOverflows=true) const noexcept {
      return sqrt(sumW2(includeOverflows));
    }
 
    double integralTo(const size_t binIndex) const noexcept {
      return integralRange(0, binIndex);
    }
 
    double integralRange(const size_t binIndex1, size_t binIndex2) const {
      assert(binIndex2 >= binIndex1);
      if (binIndex1 >= BaseT::numBins(true)) throw RangeError("binindex1 is out of range");
      if (binIndex2 >= BaseT::numBins(true)) throw RangeError("binindex2 is out of range");
      double sumw = 0;
      for (size_t idx = binIndex1; idx <= binIndex2; ++idx) {
        if (BaseT::bin(idx).isMasked())  continue;
        sumw += BaseT::bin(idx).sumW();
      }
      return sumw;
    }
 
    double integralRangeError(const size_t binIndex1, size_t binIndex2) const {
      assert(binIndex2 >= binIndex1);
      if (binIndex1 >= BaseT::numBins(true)) throw RangeError("binindex1 is out of range");
      if (binIndex2 >= BaseT::numBins(true)) throw RangeError("binindex2 is out of range");
      double sumw2 = 0;
      for (size_t idx = binIndex1; idx <= binIndex2; ++idx) {
        if (BaseT::bin(idx).isMasked())  continue;
        sumw2 += BaseT::bin(idx).sumW2();
      }
      return sumw2;
    }
 
    double numEntries(const bool includeOverflows=true) const noexcept {
      double n = 0;
      for (const auto& b : BaseT::bins(includeOverflows))
        n += b.numEntries();
      return n;
    }
 
    double effNumEntries(const bool includeOverflows=true) const noexcept {
      double n = 0;
      for (const auto& b : BaseT::bins(includeOverflows))
        n += b.effNumEntries();
      return n;
    }
 
    double sumW(const bool includeOverflows=true) const noexcept {
      double sumw = 0;
      for (const auto& b : BaseT::bins(includeOverflows))
        sumw += b.sumW();
      return sumw;
    }
 
    double sumW2(const bool includeOverflows=true) const noexcept {
      double sumw2 = 0;
      for (const auto& b : BaseT::bins(includeOverflows))
        sumw2 += b.sumW2();
      return sumw2;
    }
 
    double sumWA(const size_t dim, const bool includeOverflows=true) const {
      if (dim >= DbnN)  throw RangeError("Invalid axis int, must be in range 0..dim-1");
      double sumwa = 0;
      for (const auto& b : BaseT::bins(includeOverflows))
        sumwa += b.sumW(dim+1);
      return sumwa;
    }
 
    double sumWA2(const size_t dim, const bool includeOverflows=true) const {
      if (dim >= DbnN)  throw RangeError("Invalid axis int, must be in range 0..dim-1");
      double sumwa2 = 0;
      for (const auto& b : BaseT::bins(includeOverflows))
        sumwa2 += b.sumW2(dim+1);
      return sumwa2;
    }
 
    template<size_t dim = DbnN, typename = std::enable_if_t<(dim >= 2)>>
    double crossTerm(const size_t A1, const size_t A2, const bool includeOverflows=true) const {
      if (A1 >= DbnN || A2 >= DbnN)  throw RangeError("Invalid axis int, must be in range 0..dim-1");
      if (A1 >= A2)  throw RangeError("Indices need to be different for cross term");
      double sumw = 0;
      for (const auto& b : BaseT::bins(includeOverflows))
        sumw += b.crossTerm(A1, A2);
      return sumw;
    }
 
    double mean(size_t axisN, const bool includeOverflows=true) const noexcept {
      Dbn<DbnN> dbn;
      for (const auto& b : BaseT::bins(includeOverflows))
        dbn += b;
      return dbn.mean(axisN+1);
    }
 
    double variance(size_t axisN, const bool includeOverflows=true) const noexcept {
      Dbn<DbnN> dbn;
      for (const auto& b : BaseT::bins(includeOverflows))
        dbn += b;
      return dbn.variance(axisN+1);
    }
 
    double stdDev(size_t axisN, const bool includeOverflows=true) const noexcept {
      return std::sqrt(variance(axisN, includeOverflows));
    }
 
    double stdErr(size_t axisN, const bool includeOverflows=true) const noexcept {
      Dbn<DbnN> dbn;
      for (const auto& b : BaseT::bins(includeOverflows))
        dbn += b;
      return dbn.stdErr(axisN+1);
    }
 
    double rms(size_t axisN, const bool includeOverflows=true) const noexcept {
      Dbn<DbnN> dbn;
      for (const auto& b : BaseT::bins(includeOverflows))
        dbn += b;
      return dbn.RMS(axisN+1);
    }
 
    double dVol(const bool includeOverflows=true) const noexcept {
      double vol = 0.0;
      for (const auto& b : BaseT::bins(includeOverflows))
        vol += b.dVol();
      return vol;
    }
 
    double density(const bool includeOverflows=true) const noexcept {
      const double vol = dVol(includeOverflows);
      if (vol)  return integral(includeOverflows) / vol;
      return std::numeric_limits<double>::quiet_NaN();
    }
 
    double densityError(const bool includeOverflows=true) const noexcept {
      const double vol = dVol(includeOverflows);
      if (vol)  return integralError(includeOverflows) / vol;
      return std::numeric_limits<double>::quiet_NaN();
    }
 
    double densitySum(const bool includeOverflows=true) const noexcept {
      double rho = 0.0;
      for (const auto& b : BaseT::bins(includeOverflows))
        rho += b.sumW() / b.dVol();
      return rho;
    }
 
    double maxDensity(const bool includeOverflows=true) const noexcept {
      std::vector<double> vals;
      for (auto& b : BaseT::bins(includeOverflows))
        vals.emplace_back(b.sumW() / b.dVol());
      return *max_element(vals.begin(), vals.end());
    }
 
 
 
  private:
 
    // @brief Render information about this AO (private implementation)
    template<size_t... Is>
    void _renderYODA_aux(std::ostream& os, const int width, std::index_sequence<Is...>) const noexcept {
 
      // YODA1-style metadata
      if ( effNumEntries(true) > 0 ) {
        os << "# Mean: ";
        if (DbnN > 1) {  os << "("; }
        (( os <<  std::string(Is? ", " : "") << mean(Is, true)), ...);
        if (DbnN > 1) {  os << ")"; }
        os << "\n# Integral: " << integral(true) << "\n";
      }
 
      // render bin edges
      BaseT::_binning._renderYODA(os);
 
      // column header: content types
      os << std::setw(width) << std::left << "# sumW" << "\t";
      os << std::setw(width) << std::left << "sumW2" << "\t";
      (( os << std::setw(width) << std::left << ("sumW(A"  + std::to_string(Is+1) + ")") << "\t"
            << std::setw(width) << std::left << ("sumW2(A" + std::to_string(Is+1) + ")") << "\t"), ...);
      if constexpr (DbnN >= 2) {
        for (size_t i = 0; i < (DbnN-1); ++i) {
          for (size_t j = i+1; j < DbnN; ++j) {
            const std::string scross = "sumW(A" + std::to_string(i+1) + ",A" + std::to_string(j+1) + ")";
            os << std::setw(width) << std::left << scross << "\t";
          }
        }
      }
      os << "numEntries\n";
      // now write one row per bin
      for (const auto& b : BaseT::bins(true, true)) {
        os << std::setw(width) << std::left << b.sumW() << "\t"; // renders sumW
        os << std::setw(width) << std::left << b.sumW2() << "\t"; // renders sumW2
        ((os << std::setw(width) << std::left << b.sumW( Is+1) << "\t"
             << std::setw(width) << std::left << b.sumW2(Is+1) << "\t"), ...); // renders first moments
        if constexpr (DbnN >= 2) {
          for (size_t i = 0; i < (DbnN-1); ++i) {
            for (size_t j = i+1; j < DbnN; ++j) {
              os << std::setw(width) << std::left << b.crossTerm(i,j) << "\t";
            }
          }
        }
        os << std::setw(width) << std::left << b.numEntries() << "\n"; // renders raw event count
      }
    }
 
  public:
 
    // @brief Render information about this AO (public API)
    void _renderYODA(std::ostream& os, const int width = 13) const noexcept {
      _renderYODA_aux(os, width, std::make_index_sequence<DbnN>{});
    }
 
    // @brief Render scatter-like information about this AO
    void _renderFLAT(std::ostream& os, const int width = 13) const noexcept {
      const ScatterND<sizeof...(AxisT)+1> tmp = mkScatter();
      tmp._renderYODA(os, width);
    }
 
 
 
 
    size_t lengthContent(bool = false) const noexcept {
      return BaseT::numBins(true, true) * Dbn<DbnN>::DataSize::value;
    }
 
    std::vector<double> serializeContent(bool = false) const noexcept {
      std::vector<double> rtn;
      const size_t nBins = BaseT::numBins(true, true);
      rtn.reserve(nBins * Dbn<DbnN>::DataSize::value);
      for (size_t i = 0; i < nBins; ++i) {
        std::vector<double> bdata = BaseT::bin(i)._serializeContent();
        rtn.insert(std::end(rtn),
                   std::make_move_iterator(std::begin(bdata)),
                   std::make_move_iterator(std::end(bdata)));
      }
      return rtn;
    }
 
    void deserializeContent(const std::vector<double>& data) {
 
      constexpr size_t dbnSize = Dbn<DbnN>::DataSize::value;
      const size_t nBins = BaseT::numBins(true, true);
      if (data.size() != nBins * dbnSize)
        throw UserError("Length of serialized data should be "
                        + std::to_string(nBins * dbnSize)+"!");
 
      const auto itr = data.cbegin();
      for (size_t i = 0; i < nBins; ++i) {
        auto first = itr + i*dbnSize;
        auto last = first + dbnSize;
        BaseT::bin(i)._deserializeContent(std::vector<double>{first, last});
      }
 
    }
 
    // @}
 
 
    BinnedEstimate<AxisT...> mkEstimate(const std::string& path="",
                                        const std::string& source = "",
                                        const bool divbyvol=true) const {
 
      // @todo Should we check BaseT::nanCount() and report?
      BinnedEstimate<AxisT...> rtn(BaseT::_binning);
      for (const std::string& a : annotations()) {
        if (a != "Type")  rtn.setAnnotation(a, annotation(a));
      }
      rtn.setAnnotation("Path", path);
 
      if (BaseT::nanCount()) {
        const double nanc = BaseT::nanCount();
        const double nanw = BaseT::nanSumW();
        const double frac = nanc / (nanc + numEntries());
        const double wtot = nanw + effNumEntries();
        rtn.setAnnotation("NanFraction", frac);
        if (wtot)  rtn.setAnnotation("WeightedNanFraction", nanw/wtot);
      }
 
      for (const auto& b : BaseT::bins(true, true)) {
        if (!b.isVisible() && b.numEntries() == 0)  continue;
        const double scale = divbyvol? b.dVol() : 1.0;
        double v, e;
        if constexpr(DbnN > sizeof...(AxisT)) {
          v = b.mean(DbnN) / scale;
          e = b.stdErr(DbnN) / scale;
        }
        else {
          v = b.sumW() / scale;
          e = b.errW() / scale;
        }
        rtn.bin(b.index()).set(v, e, source);
      }
 
      return rtn;
    }
 
    template<size_t axisN, typename = std::enable_if_t< (axisN < sizeof...(AxisT)) >>
    auto mkEstimates(const std::string& path="", const std::string source = "",
                     const bool divbyvol=true, const bool includeOverflows=false) {
 
      BinnedEstimate<AxisT...> est = mkEstimate(path, source, divbyvol);
      return est.template mkEstimates<axisN>(path, includeOverflows);
    }
 
 
    auto mkScatter(const std::string& path="", const bool divbyvol=true,
                                               const bool usefocus=false,
                                               const bool includeOverflows=false,
                                               const bool includeMaskedBins=false) const {
      const BinnedEstimate<AxisT...> est = mkEstimate("", "", divbyvol);
      ScatterND<sizeof...(AxisT)+1> rtn = est.mkScatter(path, includeOverflows, includeMaskedBins);
      if (usefocus) {
        size_t idx = 0;
        for (const auto& b : BaseT::bins(includeOverflows, includeMaskedBins)) {
          auto shiftIfContinuous = [&rtn, &b, &idx](auto I) {
            using isContinuous = typename BinningT::template is_CAxis<I>;
            if constexpr (isContinuous::value) {
              const double oldMax = rtn.point(idx).max(I);
              const double oldMin = rtn.point(idx).min(I);
              const double newVal = b.mean(I+1);
              rtn.point(idx).set(I, newVal, newVal - oldMin, oldMax - newVal);
            }
          };
          MetaUtils::staticFor<BinningT::Dimension::value>(shiftIfContinuous);
          ++idx;
        }
      }
      return rtn;
    }
 
    template<size_t N = DbnN, typename = std::enable_if_t< (N == sizeof...(AxisT)+1) >>
    BinnedHisto<AxisT...> mkHisto(const std::string& path="") const {
 
      BinnedHisto<AxisT...> rtn(BaseT::_binning);
      rtn.setNanLog(BaseT::nanCount(), BaseT::nanSumW(), BaseT::nanSumW2());
      for (const std::string& a : annotations()) {
        if (a != "Type")  rtn.setAnnotation(a, annotation(a));
      }
      rtn.setAnnotation("Path", path);
 
      for (const auto& b : BaseT::bins(true)) {
        rtn.bin(b.index()) += b.template reduce<N-1>();
      }
 
      return rtn;
    }
 
    template<size_t axisN, typename = std::enable_if_t< (axisN < sizeof...(AxisT)) >>
    auto mkMarginalProfile(const std::string& path="") const {
 
      auto rtn = BaseT::template _mkBinnedT<BinnedProfile>(BaseT::_binning.template _getAxesExcept<axisN>());
      rtn.setNanLog(BaseT::nanCount(), BaseT::nanSumW(), BaseT::nanSumW2());
      for (const std::string& a : annotations()) {
        if (a != "Type")  rtn.setAnnotation(a, annotation(a));
      }
      rtn.setAnnotation("Path", path);
 
      auto collapseStorageBins =
        [&oldBinning = BaseT::_binning, &oldBins = BaseT::_bins, &rtn](auto I, auto dbnRed) {
 
        auto collapse = [&oldBins, &rtn](const auto& binsIndicesToMerge, auto axis) {
          assert(rtn.numBins(true) == binsIndicesToMerge.size());
 
          // for any given pivot, add the content
          // from the old slice to the new slice
          for (size_t i = 0; i < rtn.numBins(true); ++i) {
            auto& pivotBin = rtn.bin(i);
            auto& binToAppend = oldBins[binsIndicesToMerge[i]];
            pivotBin += binToAppend.template reduce<axis>();
          }
        };
 
        // get bin slice for any given bin along the axis that is to be
        // collapsed, then copy the values into the new binning
        ssize_t nBinRowsToBeMerged = oldBinning.numBinsAt(I);
        while (nBinRowsToBeMerged--) {
          collapse(oldBinning.sliceIndices(I, nBinRowsToBeMerged), dbnRed);
        }
      };
      auto dbnRed = std::integral_constant<size_t, (sizeof...(AxisT) == DbnN)? DbnN : axisN>();
      (void)collapseStorageBins(std::integral_constant<std::size_t, axisN>(), dbnRed);
 
      return rtn;
    }
 
    template<size_t axisN, typename = std::enable_if_t< (axisN < sizeof...(AxisT)) >>
    auto mkMarginalHisto(const std::string& path="") const {
 
      if constexpr (DbnN != sizeof...(AxisT)) {
        // Case 1: BP(N+1) -> BH(N+1) -> BHN
        return mkHisto().template mkMarginalHisto<axisN>(path);
      }
      else {
        // Case 2: BH(N+1) -> BHN
 
        auto rtn = BaseT::template _mkBinnedT<BinnedHisto>(BaseT::_binning.template _getAxesExcept<axisN>());
        rtn.setNanLog(BaseT::nanCount(), BaseT::nanSumW(), BaseT::nanSumW2());
        for (const std::string& a : annotations()) {
          if (a != "Type")  rtn.setAnnotation(a, annotation(a));
        }
        rtn.setAnnotation("Path", path);
 
        auto collapseStorageBins =
          [&oldBinning = BaseT::_binning, &oldBins = BaseT::_bins, &rtn](auto I, auto dbnRed) {
 
          auto collapse = [&oldBins, &rtn](const auto& binsIndicesToMerge, auto axis) {
            assert(rtn.numBins(true) == binsIndicesToMerge.size());
 
            // for any given pivot, add the content
            // from the old slice to the new slice
            for (size_t i = 0; i < rtn.numBins(true); ++i) {
              auto& pivotBin = rtn.bin(i);
              auto& binToAppend = oldBins[binsIndicesToMerge[i]];
              pivotBin += binToAppend.template reduce<axis>();
            }
          };
 
          // get bin slice for any given bin along the axis that is to be
          // collapsed, then copy the values into the new binning
          ssize_t nBinRowsToBeMerged = oldBinning.numBinsAt(I);
          while (nBinRowsToBeMerged--) {
            collapse(oldBinning.sliceIndices(I, nBinRowsToBeMerged), dbnRed);
          }
        };
        // collapse Dbn along axisN
        auto dbnRed = std::integral_constant<size_t, axisN>();
        (void)collapseStorageBins(std::integral_constant<std::size_t, axisN>(), dbnRed);
 
        return rtn;
      }
    }
 
 
    template<size_t axisN, typename = std::enable_if_t< (axisN < sizeof...(AxisT) &&
                                                         sizeof...(AxisT)>=2 &&
                                                         DbnN > sizeof...(AxisT)) >>
    auto mkProfiles(const std::string& path="", const bool includeOverflows=false) const {
 
      // Need to provide a prescription for how to add the two bin contents
      auto how2add = [](auto& pivot, const BinType& toCopy) { pivot = toCopy.template reduce<axisN>(); };
      auto rtn = BaseT::template mkBinnedSlices<axisN, BinnedProfile>(how2add, includeOverflows);
      for (const std::string& a : annotations()) {
        if (a == "Type")  continue;
        for (size_t i = 0; i < rtn.size(); ++i) {
          rtn[i].setAnnotation(a, annotation(a));
        }
      }
      for (size_t i = 0; i < rtn.size(); ++i) {
        rtn[i].setAnnotation("Path", path);
      }
      return rtn;
    }
 
 
    template<size_t axisN, typename = std::enable_if_t< (axisN < sizeof...(AxisT) && sizeof...(AxisT)>=2) >>
    auto mkHistos(const std::string& path="", const bool includeOverflows=false) const {
 
      if constexpr (DbnN != sizeof...(AxisT)) {
        // Case 1: BP(N+1) -> BH(N+1) -> BHN
        return mkHisto().template mkHistos<axisN>(path, includeOverflows);
      }
      else {
        // Case 2: BH(N+1) -> BHN
 
        // Need to provide a prescription for how to add the two bin contents
        auto how2add = [](auto& pivot, const BinType& toCopy) { pivot = toCopy.template reduce<axisN>(); };
        auto rtn = BaseT::template mkBinnedSlices<axisN,BinnedHisto>(how2add, includeOverflows);
        for (const std::string& a : annotations()) {
          if (a == "Type")  continue;
          for (size_t i = 0; i < rtn.size(); ++i) {
            rtn[i].setAnnotation(a, annotation(a));
          }
        }
        for (size_t i = 0; i < rtn.size(); ++i) {
          rtn[i].setAnnotation("Path", path);
        }
        return rtn;
      }
    }
 
 
    AnalysisObject* mkInert(const std::string& path = "",
                            const std::string& source = "") const noexcept {
      return mkEstimate(path, source).newclone();
    }
 
 
    private:
 
    template<size_t... Is>
    BinningT _mkBinning(const std::vector<size_t>& nBins,
                        const std::vector<std::pair<double, double>>& limitsLowUp,
                        std::index_sequence<Is...>) const {
      return BinningT({((void)Is, Axis<AxisT>(nBins[Is], limitsLowUp[Is].first, limitsLowUp[Is].second))...});
    }
 
    template<size_t... Is>
    BinningT _mkBinning(const ScatterND<sizeof...(AxisT)+1>& s, std::index_sequence<Is...>) const {
      return BinningT(Axis<AxisT>(s.edges(Is))...);
    }
 
  };
 
 
 
 
  template<size_t DbnN, typename... AxisT>
  inline BinnedDbn<DbnN, AxisT...>
  operator + (BinnedDbn<DbnN, AxisT...> first, BinnedDbn<DbnN, AxisT...>&& second) {
    first += std::move(second);
    return first;
  }
  //
  template <size_t DbnN, typename... AxisT>
  inline BinnedDbn<DbnN, AxisT...>
  operator + (BinnedDbn<DbnN, AxisT...> first, const BinnedDbn<DbnN, AxisT...>& second) {
    first += second;
    return first;
  }
 
 
  template <size_t DbnN, typename... AxisT>
  inline BinnedDbn<DbnN, AxisT...>
  operator - (BinnedDbn<DbnN, AxisT...> first, BinnedDbn<DbnN, AxisT...>&& second) {
    first -= std::move(second);
    return first;
  }
  //
  template <size_t DbnN, typename... AxisT>
  inline BinnedDbn<DbnN, AxisT...>
  operator - (BinnedDbn<DbnN, AxisT...> first, const BinnedDbn<DbnN, AxisT...>& second) {
    first -= second;
    return first;
  }
 
 
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  divide(const BinnedDbn<DbnN, AxisT...>& numer, const BinnedDbn<DbnN, AxisT...>& denom) {
 
    if (numer != denom) {
      throw BinningError("Arithmetic operation requires compatible binning!");
    }
 
    BinnedEstimate<AxisT...> rtn = numer.mkEstimate();
    if (numer.path() == denom.path())  rtn.setPath(numer.path());
    if (rtn.hasAnnotation("ScaledBy")) rtn.rmAnnotation("ScaledBy");
 
    for (const auto& b_num : numer.bins(true, true)) {
      const size_t idx = b_num.index();
      const auto& b_den = denom.bin(idx);
      double v, e;
      if (!b_den.effNumEntries()) {
        v = std::numeric_limits<double>::quiet_NaN();
        e = std::numeric_limits<double>::quiet_NaN();
      }
      else {
        if constexpr(DbnN > sizeof...(AxisT)) {
          v = b_num.mean(DbnN) / b_den.mean(DbnN);
          const double e_num = b_num.effNumEntries()? b_num.relStdErr(DbnN) : 0;
          const double e_den = b_den.effNumEntries()? b_den.relStdErr(DbnN) : 0;
          e = fabs(v) * sqrt(sqr(e_num) + sqr(e_den));
        }
        else {
          v = b_num.sumW() / b_den.sumW();
          const double e_num = b_num.sumW()? b_num.relErrW() : 0;
          const double e_den = b_den.sumW()? b_den.relErrW() : 0;
          e = fabs(v) * sqrt(sqr(e_num) + sqr(e_den));
        }
      }
      rtn.bin(idx).set(v, {-e, e}); // @todo put "stats" as source?
    }
    rtn.maskBins(denom.maskedBins(), true);
 
    return rtn;
  }
  //
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  operator / (const BinnedDbn<DbnN, AxisT...>& numer, const BinnedDbn<DbnN, AxisT...>& denom) {
    return divide(numer, denom);
  }
  //
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  operator / (const BinnedDbn<DbnN, AxisT...>& numer, BinnedDbn<DbnN, AxisT...>&& denom) {
    return divide(numer, std::move(denom));
  }
  //
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  operator / (BinnedDbn<DbnN, AxisT...>&& numer, const BinnedDbn<DbnN, AxisT...>& denom) {
    return divide(std::move(numer), denom);
  }
  //
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  operator / (BinnedDbn<DbnN, AxisT...>&& numer, BinnedDbn<DbnN, AxisT...>&& denom) {
    return divide(std::move(numer), std::move(denom));
  }
 
 
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  efficiency(const BinnedDbn<DbnN, AxisT...>& accepted, const BinnedDbn<DbnN, AxisT...>& total) {
 
    if (accepted != total) {
      throw BinningError("Arithmetic operation requires compatible binning!");
    }
 
    BinnedEstimate<AxisT...> rtn = divide(accepted, total);
 
    for (const auto& b_acc : accepted.bins(true, true)) {
      const auto& b_tot = total.bin(b_acc.index());
      auto& b_rtn = rtn.bin(b_acc.index());
 
      // Check that the numerator is consistent with being a subset of the denominator
      if (b_acc.numEntries() > b_tot.numEntries())
        throw UserError("Attempt to calculate an efficiency when the numerator is not a subset of the denominator: "
                        + Utils::toStr(b_acc.numEntries()) + " entries / " + Utils::toStr(b_tot.numEntries()) + " entries");
 
      // If no entries on the denominator, set eff = err = 0 and move to the next bin
      double eff = std::numeric_limits<double>::quiet_NaN();
      double err = std::numeric_limits<double>::quiet_NaN();
      try {
        if (b_tot.sumW()) {
          eff = b_rtn.val();
          err = sqrt(abs( ((1-2*eff)*b_acc.sumW2() + sqr(eff)*b_tot.sumW2()) / sqr(b_tot.sumW()) ));
        }
      } catch (const LowStatsError& e) {
        //
      }
 
      b_rtn.setErr({-err, err}); // @todo put "stats" as source?
    }
    return rtn;
  }
 
 
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  asymm(const BinnedDbn<DbnN, AxisT...>& a, const BinnedDbn<DbnN, AxisT...>& b) {
    return (a-b) / (a+b);
  }
 
 
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  mkIntegral(const BinnedDbn<DbnN, AxisT...>& histo, const bool includeOverflows = true) {
 
    BinnedEstimate<AxisT...> rtn = histo.mkEstimate();
 
    double sumW = 0.0, sumW2 = 0.0;
    for (const auto& b : histo.bins(includeOverflows)) {
      sumW  += b.sumW();
      sumW2 += b.sumW2();
      const double e = sqrt(sumW2);
      rtn.bin(b.index()).set(sumW, {-e, e});
    }
 
    return rtn;
  }
 
 
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  mkIntegralEff(const BinnedDbn<DbnN, AxisT...>& histo, const bool includeOverflows = true) {
 
    BinnedEstimate<AxisT...> rtn = mkIntegral(histo, includeOverflows);
    const double integral = histo.integral(includeOverflows);
 
    // If the integral is empty, the (integrated) efficiency values may as well all be zero, so return here
    if (!integral) return rtn;
 
    const double integral_err = histo.integralError(includeOverflows);
    for (const auto& b : rtn.bins(includeOverflows)) {
      const double eff = b.val() / integral;
      const double err = sqrt(std::abs( ((1-2*eff)*sqr(b.relTotalErrAvg()) + sqr(eff)*sqr(integral_err)) / sqr(integral) ));
      b.set(eff, {-err,err});
    }
 
    return rtn;
  }
 
 
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  add(const BinnedDbn<DbnN, AxisT...>& dbn, const BinnedEstimate<AxisT...>& est) {
    return dbn.mkEstimate() + est;
  }
  //
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  operator + (const BinnedDbn<DbnN, AxisT...>& dbn, const BinnedEstimate<AxisT...>& est) {
    return add(dbn, est);
  }
  //
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  operator + (BinnedDbn<DbnN, AxisT...>&& dbn, const BinnedEstimate<AxisT...>& est) {
    return add(std::move(dbn), est);
  }
  //
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  operator + (const BinnedDbn<DbnN, AxisT...>& dbn, BinnedEstimate<AxisT...>&& est) {
    return add(dbn, std::move(est));
  }
  //
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  operator + (BinnedDbn<DbnN, AxisT...>&& dbn, BinnedEstimate<AxisT...>&& est) {
    return add(std::move(dbn), std::move(est));
  }
 
 
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  subtract(const BinnedDbn<DbnN, AxisT...>& dbn, const BinnedEstimate<AxisT...>& est) {
    return dbn.mkEstimate() - est;
  }
  //
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  operator - (const BinnedDbn<DbnN, AxisT...>& dbn, const BinnedEstimate<AxisT...>& est) {
    return subtract(dbn, est);
  }
  //
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  operator - (BinnedDbn<DbnN, AxisT...>&& dbn, const BinnedEstimate<AxisT...>& est) {
    return subtract(std::move(dbn), est);
  }
  //
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  operator - (const BinnedDbn<DbnN, AxisT...>& dbn, BinnedEstimate<AxisT...>&& est) {
    return subtract(dbn, std::move(est));
  }
  //
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  operator - (BinnedDbn<DbnN, AxisT...>&& dbn, BinnedEstimate<AxisT...>&& est) {
    return subtract(std::move(dbn), std::move(est));
  }
 
 
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  divide(const BinnedDbn<DbnN, AxisT...>& dbn, const BinnedEstimate<AxisT...>& est) {
    return dbn.mkEstimate() / est;
  }
  //
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  operator / (const BinnedDbn<DbnN, AxisT...>& dbn, const BinnedEstimate<AxisT...>& est) {
    return divide(dbn, est);
  }
  //
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  operator / (BinnedDbn<DbnN, AxisT...>&& dbn, const BinnedEstimate<AxisT...>& est) {
    return divide(std::move(dbn), est);
  }
  //
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  operator / (const BinnedDbn<DbnN, AxisT...>& dbn, BinnedEstimate<AxisT...>&& est) {
    return divide(dbn, std::move(est));
  }
  //
  template <size_t DbnN, typename... AxisT>
  inline BinnedEstimate<AxisT...>
  operator / (BinnedDbn<DbnN, AxisT...>&& dbn, BinnedEstimate<AxisT...>&& est) {
    return divide(std::move(dbn), std::move(est));
  }
 
 
}
 
#endif