RivetYODA.hh

#ifndef RIVET_RIVETYODA_HH
#define RIVET_RIVETYODA_HH
 
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/Tools/TypeTraits.hh"
#include "YODA/AnalysisObject.h"
#include "YODA/BinnedEstimate.h"
#include "YODA/Counter.h"
#include "YODA/Estimate0D.h"
#include "YODA/Histo.h"
#include "YODA/Profile.h"
#include "YODA/Scatter.h"
 
// Use execinfo for backtrace if available
#ifdef HAVE_EXECINFO_H
  #include <execinfo.h>
#endif
 
#include <map>
#include <unordered_map>
#include <valarray>
 
 
namespace YODA {
 
  template <size_t DbnN, typename... AxisT>
  using BinnedDbnPtr = std::shared_ptr<YODA::BinnedDbn<DbnN, AxisT...>>;
 
  template <typename... AxisT>
  using BinnedHistoPtr = BinnedDbnPtr<sizeof...(AxisT), AxisT...>;
 
  template <typename... AxisT>
  using BinnedProfilePtr = BinnedDbnPtr<sizeof...(AxisT) + 1, AxisT...>;
 
  template <typename... AxisT>
  using BinnedEstimatePtr = std::shared_ptr<YODA::BinnedEstimate<AxisT...>>;
 
  template <size_t N>
  using ScatterNDPtr = std::shared_ptr<YODA::ScatterND<N>>;
 
  using AnalysisObjectPtr = std::shared_ptr<YODA::AnalysisObject>;
  using CounterPtr = std::shared_ptr<YODA::Counter>;
  using Estimate0DPtr = std::shared_ptr<YODA::Estimate0D>;
  using Histo1DPtr = BinnedHistoPtr<double>;
  using Histo2DPtr = BinnedHistoPtr<double, double>;
  using Histo3DPtr = BinnedHistoPtr<double, double, double>;
  using Profile1DPtr = BinnedProfilePtr<double>;
  using Profile2DPtr = BinnedProfilePtr<double, double>;
  using Profile3DPtr = BinnedProfilePtr<double, double, double>;
  using Estimate1DPtr = BinnedEstimatePtr<double>;
  using Estimate2DPtr = BinnedEstimatePtr<double, double>;
  using Estimate3DPtr = BinnedEstimatePtr<double, double, double>;
  using Scatter1DPtr = ScatterNDPtr<1>;
  using Scatter2DPtr = ScatterNDPtr<2>;
  using Scatter3DPtr = ScatterNDPtr<3>;
 
}
 
namespace Rivet {

  template <typename T>
  bool copyAO(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst, const double scale = 1.0) {
    if (dst->hasAnnotation("Type") && src->type() != dst->type()) {
      throw YODA::LogicError("Operation requries types to be the same!");
    }
    for (const std::string& a : src->annotations()) {
      dst->setAnnotation(a, src->annotation(a));
    }
    shared_ptr<T> dstPtr = std::static_pointer_cast<T>(dst);
    *dstPtr = *std::static_pointer_cast<T>(src);
    if constexpr (isFillable<T>::value) {
      dstPtr->scaleW(scale);
    }
    return true;
  }
 

  struct TypeBaseHandle {
 
    TypeBaseHandle() = default;
 
    virtual ~TypeBaseHandle() { }
 
    virtual bool copyAO(YODA::AnalysisObjectPtr src,
                        YODA::AnalysisObjectPtr dst,
                        const double scale = 1.0) const = 0;
 
    virtual bool addAO(YODA::AnalysisObjectPtr src,
                       YODA::AnalysisObjectPtr& dst,
                       const double scale = 1.0) const = 0;
  };
 

  template <typename T>
  struct TypeHandle : public TypeBaseHandle {
 
    bool addAO(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr& dst, const double scale = 1.0) const {
      if constexpr (isFillable<T>::value) {
        std::shared_ptr<T> srcPtr = std::static_pointer_cast<T>(src);
        srcPtr->scaleW(scale);
        if (dst == nullptr) {
          dst = src;
          return true;
        }
        try {
          *std::static_pointer_cast<T>(dst) += *srcPtr;
        }
        catch (YODA::BinningError&) {
          return false;
        }
        return true;
      }
      else if (dst == nullptr) {
        dst = src;
        return true;
      }
      return false;
    }
 
    bool copyAO(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst, const double scale = 1.0) const {
      return ::Rivet::copyAO<T>(src, dst, scale);
    }
  };
 


  using Weight = double;

  template <typename T>
  using Fill = pair<typename T::FillType, Weight>;

  template <typename T>
  using Fills = vector<Fill<T>>;
 

  template <typename T>
  class FillCollector;
 

  template <>
  class FillCollector<YODA::Counter> : public YODA::Counter {
  public:
 
    using YAO = YODA::Counter;
    using Ptr = shared_ptr<FillCollector<YAO>>;
    using YAO::operator=;
 
    FillCollector()
        : YAO() { }

    FillCollector(typename YAO::Ptr yao)
        : YAO(yao->path()) { }

    int fill(const double weight = 1.0, const double fraction = 1.0) {
      (void)fraction; // suppress unused variable warning
      _fills.insert(_fills.end(), {YAO::FillType(), weight});
      return 0;
    }

    void reset() {
      _fills.clear();
    }

    const Fills<YAO>& fills() const {
      return _fills;
    }
 
  private:
 
    Fills<YAO> _fills;
  };
 

  template <size_t DbnN, typename... AxisT>
  class FillCollector<YODA::BinnedDbn<DbnN, AxisT...>> : public YODA::BinnedDbn<DbnN, AxisT...> {
  public:
 
    using YAO = YODA::BinnedDbn<DbnN, AxisT...>;
    using Ptr = shared_ptr<FillCollector<YAO>>;
    using YAO::operator=;
 
    FillCollector()
        : YAO() { }

    FillCollector(typename YAO::Ptr yao)
        : YAO(yao->binning()) {
      YAO::setPath(yao->path());
    }

    int fill(typename YAO::FillType&& fillCoords, const double weight = 1.0, const double fraction = 1.0) {
      (void)fraction; // suppress unused variable warning
      if (YODA::containsNan(fillCoords)) {
        _fills.insert(_fills.end(), {std::move(fillCoords), weight});
        return -1;
      }
      // Could be that DbnN > number of binned axes, so should
      // extract subset of bin coordinates to pinpoint bin
      typename YAO::BinningT::EdgeTypesTuple binCoords{};
      auto extractBinCoords = [&binCoords, &fillCoords](auto I) {
        std::get<I>(binCoords) = std::get<I>(fillCoords);
      };
      MetaUtils::staticFor<sizeof...(AxisT)>(extractBinCoords);
      _fills.insert(_fills.end(), {std::move(fillCoords), weight});
      return (int)YAO::_binning.globalIndexAt(binCoords);
    }

    void reset() noexcept {
      _fills.clear();
    }

    const Fills<YAO>& fills() const {
      return _fills;
    }
 
  private:
 
    Fills<YAO> _fills;
  };

  template <typename AxisT>
  class FillCollector<YODA::BinnedDbn<1, AxisT>> : public YODA::BinnedDbn<1, AxisT> {
  public:
 
    using YAO = YODA::BinnedDbn<1, AxisT>;
    using Ptr = shared_ptr<FillCollector<YAO>>;
    using YAO::operator=;
 
    FillCollector()
        : YAO() { }

    FillCollector(typename YAO::Ptr yao)
        : YAO(yao->binning()) {
      YAO::setPath(yao->path());
    }

    int fill(typename YAO::FillType&& fillCoords, const double weight = 1.0, const double fraction = 1.0) {
      (void)fraction; // suppress unused variable warning
      if (YODA::containsNan(fillCoords)) {
        _fills.insert(_fills.end(), {std::move(fillCoords), weight});
        return -1;
      }
      // Could be that DbnN > number of binned axes, so should
      // extract subset of bin coordinates to pinpoint bin
      typename YAO::BinningT::EdgeTypesTuple binCoords{};
      auto extractBinCoords = [&binCoords, &fillCoords](auto I) {
        std::get<I>(binCoords) = std::get<I>(fillCoords);
      };
      MetaUtils::staticFor<1>(extractBinCoords);
      _fills.insert(_fills.end(), {std::move(fillCoords), weight});
      return (int)YAO::_binning.globalIndexAt(binCoords);
    }
    //
    int fill(const AxisT x, const double weight = 1.0, const double fraction = 1.0) {
      return fill(typename YAO::FillType{x}, weight, fraction);
    }

    void reset() noexcept {
      _fills.clear();
    }

    const Fills<YAO>& fills() const {
      return _fills;
    }
 
  private:
 
    Fills<YAO> _fills;
  };

  template <typename AxisT1, typename AxisT2>
  class FillCollector<YODA::BinnedDbn<2, AxisT1, AxisT2>> : public YODA::BinnedDbn<2, AxisT1, AxisT2> {
  public:
 
    using YAO = YODA::BinnedDbn<2, AxisT1, AxisT2>;
    using Ptr = shared_ptr<FillCollector<YAO>>;
    using YAO::operator=;
 
    FillCollector()
        : YAO() { }

    FillCollector(typename YAO::Ptr yao)
        : YAO(yao->binning()) {
      YAO::setPath(yao->path());
    }

    int fill(typename YAO::FillType&& fillCoords, const double weight = 1.0, const double fraction = 1.0) {
      (void)fraction; // suppress unused variable warning
      if (YODA::containsNan(fillCoords)) {
        _fills.insert(_fills.end(), {std::move(fillCoords), weight});
        return -1;
      }
      // Could be that DbnN > number of binned axes, so should
      // extract subset of bin coordinates to pinpoint bin
      typename YAO::BinningT::EdgeTypesTuple binCoords{};
      auto extractBinCoords = [&binCoords, &fillCoords](auto I) {
        std::get<I>(binCoords) = std::get<I>(fillCoords);
      };
      MetaUtils::staticFor<1>(extractBinCoords);
      _fills.insert(_fills.end(), {std::move(fillCoords), weight});
      return (int)YAO::_binning.globalIndexAt(binCoords);
    }
    //
    int fill(const AxisT1 x, const AxisT2 y, const double weight = 1.0, const double fraction = 1.0) {
      return fill(typename YAO::FillType{x, y}, weight, fraction);
    }

    void reset() noexcept {
      _fills.clear();
    }

    const Fills<YAO>& fills() const {
      return _fills;
    }
 
  private:
 
    Fills<YAO> _fills;
  };

  template <typename AxisT1, typename AxisT2, typename AxisT3>
  class FillCollector<YODA::BinnedDbn<3, AxisT1, AxisT2, AxisT3>>
      : public YODA::BinnedDbn<3, AxisT1, AxisT2, AxisT3> {
  public:
 
    using YAO = YODA::BinnedDbn<3, AxisT1, AxisT2, AxisT3>;
    using Ptr = shared_ptr<FillCollector<YAO>>;
    using YAO::operator=;
 
    FillCollector()
        : YAO() { }

    FillCollector(typename YAO::Ptr yao)
        : YAO(yao->binning()) {
      YAO::setPath(yao->path());
    }

    int fill(typename YAO::FillType&& fillCoords, const double weight = 1.0, const double fraction = 1.0) {
      (void)fraction; // suppress unused variable warning
      if (YODA::containsNan(fillCoords)) {
        _fills.insert(_fills.end(), {std::move(fillCoords), weight});
        return -1;
      }
      // Could be that DbnN > number of binned axes, so should
      // extract subset of bin coordinates to pinpoint bin
      typename YAO::BinningT::EdgeTypesTuple binCoords{};
      auto extractBinCoords = [&binCoords, &fillCoords](auto I) {
        std::get<I>(binCoords) = std::get<I>(fillCoords);
      };
      MetaUtils::staticFor<1>(extractBinCoords);
      _fills.insert(_fills.end(), {std::move(fillCoords), weight});
      return (int)YAO::_binning.globalIndexAt(binCoords);
    }
    //
    int fill(const AxisT1 x,
             const AxisT2 y,
             const AxisT3 z,
             const double weight = 1.0,
             const double fraction = 1.0) {
      return fill(typename YAO::FillType{x, y, z}, weight, fraction);
    }

    void reset() noexcept {
      _fills.clear();
    }

    const Fills<YAO>& fills() const {
      return _fills;
    }
 
  private:
 
    Fills<YAO> _fills;
  };

  template <typename AxisT>
  class FillCollector<YODA::BinnedDbn<2, AxisT>> : public YODA::BinnedDbn<2, AxisT> {
  public:
 
    using YAO = YODA::BinnedDbn<2, AxisT>;
    using Ptr = shared_ptr<FillCollector<YAO>>;
    using YAO::operator=;
 
    FillCollector()
        : YAO() { }

    FillCollector(typename YAO::Ptr yao)
        : YAO(yao->binning()) {
      YAO::setPath(yao->path());
    }

    int fill(typename YAO::FillType&& fillCoords, const double weight = 1.0, const double fraction = 1.0) {
      (void)fraction; // suppress unused variable warning
      if (YODA::containsNan(fillCoords)) {
        _fills.insert(_fills.end(), {std::move(fillCoords), weight});
        return -1;
      }
      // Could be that DbnN > number of binned axes, so should
      // extract subset of bin coordinates to pinpoint bin
      typename YAO::BinningT::EdgeTypesTuple binCoords{};
      auto extractBinCoords = [&binCoords, &fillCoords](auto I) {
        std::get<I>(binCoords) = std::get<I>(fillCoords);
      };
      MetaUtils::staticFor<1>(extractBinCoords);
      _fills.insert(_fills.end(), {std::move(fillCoords), weight});
      return (int)YAO::_binning.globalIndexAt(binCoords);
    }
    //
    int fill(const AxisT x, const double y, const double weight = 1.0, const double fraction = 1.0) {
      return fill(typename YAO::FillType{x, y}, weight, fraction);
    }

    void reset() noexcept {
      _fills.clear();
    }

    const Fills<YAO>& fills() const {
      return _fills;
    }
 
  private:
 
    Fills<YAO> _fills;
  };

  template <typename AxisT1, typename AxisT2>
  class FillCollector<YODA::BinnedDbn<3, AxisT1, AxisT2>> : public YODA::BinnedDbn<3, AxisT1, AxisT2> {
  public:
 
    using YAO = YODA::BinnedDbn<3, AxisT1, AxisT2>;
    using Ptr = shared_ptr<FillCollector<YAO>>;
    using YAO::operator=;
 
    FillCollector()
        : YAO() { }

    FillCollector(typename YAO::Ptr yao)
        : YAO(yao->binning()) {
      YAO::setPath(yao->path());
    }

    int fill(typename YAO::FillType&& fillCoords, const double weight = 1.0, const double fraction = 1.0) {
      (void)fraction; // suppress unused variable warning
      if (YODA::containsNan(fillCoords)) {
        _fills.insert(_fills.end(), {std::move(fillCoords), weight});
        return -1;
      }
      // Could be that DbnN > number of binned axes, so should
      // extract subset of bin coordinates to pinpoint bin
      typename YAO::BinningT::EdgeTypesTuple binCoords{};
      auto extractBinCoords = [&binCoords, &fillCoords](auto I) {
        std::get<I>(binCoords) = std::get<I>(fillCoords);
      };
      MetaUtils::staticFor<1>(extractBinCoords);
      _fills.insert(_fills.end(), {std::move(fillCoords), weight});
      return (int)YAO::_binning.globalIndexAt(binCoords);
    }
    //
    int fill(const AxisT1 x,
             const AxisT2 y,
             const double z,
             const double weight = 1.0,
             const double fraction = 1.0) {
      return fill(typename YAO::FillType{x, y, z}, weight, fraction);
    }

    void reset() noexcept {
      _fills.clear();
    }

    const Fills<YAO>& fills() const {
      return _fills;
    }
 
  private:
 
    Fills<YAO> _fills;
  };

  template <typename AxisT1, typename AxisT2, typename AxisT3>
  class FillCollector<YODA::BinnedDbn<4, AxisT1, AxisT2, AxisT3>>
      : public YODA::BinnedDbn<4, AxisT1, AxisT2, AxisT3> {
  public:
 
    using YAO = YODA::BinnedDbn<4, AxisT1, AxisT2, AxisT3>;
    using Ptr = shared_ptr<FillCollector<YAO>>;
    using YAO::operator=;
 
    FillCollector()
        : YAO() { }

    FillCollector(typename YAO::Ptr yao)
        : YAO(yao->binning()) {
      YAO::setPath(yao->path());
    }

    int fill(typename YAO::FillType&& fillCoords, const double weight = 1.0, const double fraction = 1.0) {
      (void)fraction; // suppress unused variable warning
      if (YODA::containsNan(fillCoords)) {
        _fills.insert(_fills.end(), {std::move(fillCoords), weight});
        return -1;
      }
      // Could be that DbnN > number of binned axes, so should
      // extract subset of bin coordinates to pinpoint bin
      typename YAO::BinningT::EdgeTypesTuple binCoords{};
      auto extractBinCoords = [&binCoords, &fillCoords](auto I) {
        std::get<I>(binCoords) = std::get<I>(fillCoords);
      };
      MetaUtils::staticFor<1>(extractBinCoords);
      _fills.insert(_fills.end(), {std::move(fillCoords), weight});
      return (int)YAO::_binning.globalIndexAt(binCoords);
    }
    //
    int fill(const AxisT1 x,
             const AxisT2 y,
             const AxisT3 z,
             const double zPlus,
             const double weight = 1.0,
             const double fraction = 1.0) {
      return fill(typename YAO::FillType{x, y, z, zPlus}, weight, fraction);
    }

    void reset() noexcept {
      _fills.clear();
    }

    const Fills<YAO>& fills() const {
      return _fills;
    }
 
  private:
 
    Fills<YAO> _fills;
  };
 

  template <>
  class FillCollector<YODA::Estimate0D> : public YODA::Estimate0D {
  public:
 
    using YAO = YODA::Estimate0D;
    using Ptr = shared_ptr<FillCollector<YAO>>;
    using YAO::operator=;
 
    FillCollector()
        : YAO() { }

    FillCollector(typename YAO::Ptr yao)
        : YAO(yao->path()) { }
  };

  template <typename... AxisT>
  class FillCollector<YODA::BinnedEstimate<AxisT...>> : public YODA::BinnedEstimate<AxisT...> {
  public:
 
    using YAO = YODA::BinnedEstimate<AxisT...>;
    using Ptr = shared_ptr<FillCollector<YAO>>;
    using YAO::operator=;
 
    FillCollector()
        : YAO() { }

    FillCollector(typename YAO::Ptr yao)
        : YAO(yao->path()) { }
  };

  template <size_t N>
  class FillCollector<YODA::ScatterND<N>> : public YODA::ScatterND<N> {
  public:
 
    using YAO = YODA::ScatterND<N>;
    using Ptr = shared_ptr<FillCollector<YAO>>;
    using YAO::operator=;
 
    FillCollector()
        : YAO() { }

    FillCollector(typename YAO::Ptr yao)
        : YAO(yao->path()) { }
  };


  namespace {
 
    template <typename... Args>
    double distance(const std::tuple<Args...>& a, const std::tuple<Args...>& b) {
      double rtn = 0;
      auto calculateDistance = [&](auto I) {
        if constexpr (std::is_floating_point<std::tuple_element_t<I, std::tuple<Args...>>>::value) {
          rtn += Rivet::sqr(std::get<I>(a) - std::get<I>(b));
        }
      };
      MetaUtils::staticFor<sizeof...(Args)>(calculateDistance);
      return rtn;
    }

    template <typename T>
    struct SubwindowType;
    //
    template <typename... Args>
    struct SubwindowType<std::tuple<Args...>> {
      using type = YODA::Binning<std::decay_t<decltype(std::declval<YODA::Axis<Args>>())>...>;
    };

    template <typename T>
    struct WindowType;
    //
    template <typename... Args>
    struct WindowType<std::tuple<Args...>> {
      using type = std::tuple<std::decay_t<decltype(std::declval<vector<Args>>())>...>;
    };

    template <typename YAO>
    vector<std::tuple<typename YAO::FillType, std::valarray<double>, double>> mkFillWindows(
        const vector<typename FillCollector<YAO>::Ptr>& subevents,
        const std::string& path,
        const vector<std::valarray<double>>& weights,
        const double fsmear = -1.0) {
 
      // Placeholder for the return type: a tuple of {FillType, multi-weights, fill fraction}
      vector<std::tuple<typename YAO::FillType, std::valarray<double>, double>> rtn;
 
      using WindowT = typename WindowType<typename YAO::FillType>::type;
      using SubwindowT = typename SubwindowType<typename YAO::FillType>::type;
      SubwindowT subwindows;
 
      // calculate the maximum number of fills
      size_t fidx = 0;
      const size_t nSubevents = subevents.size();
      for (size_t sidx = 0; sidx < nSubevents; ++sidx) {
        const size_t nFills = subevents[sidx]->fills().size();
        if (nFills == 0) continue;
        fidx = std::max(fidx, nFills);
      }
 
      while (fidx--) { // loop over fills, back to front
 
        WindowT windows;
 
        auto constructWindows = [&](auto I) {
          using FillAxisT = typename SubwindowT::template getAxisT<I>;
          using isContinuous = typename SubwindowT::template is_CAxis<I>;
          using EdgeT = typename FillAxisT::EdgeT;
 
          if constexpr (!isContinuous::value) { // discrete axes don't need smearing
            vector<EdgeT> edges;
            edges.resize(nSubevents);
            for (const auto& subevt : subevents) {
              if (fidx >= subevt->fills().size()) continue;
              edges.push_back(std::get<I>(subevt->fills()[fidx].first));
            }
            subwindows.template axis<I>() = FillAxisT(edges);
            std::get<I>(windows) = std::move(edges);
            return;
          }
          else { // continuous axes need windowing
            vector<EdgeT>& edges = std::get<I>(windows);
            edges.reserve(2 * nSubevents);
 
            if constexpr (I < YAO::BinningT::Dimension::value) {
              // this fill axis is binned: window sizes will depend on bin width
              const auto& axis = subevents[0]->binning().template axis<I>();
              size_t over = 0, under = 0;
              const EdgeT edgeMax = subevents[0]->template max<I>();
              const EdgeT edgeMin = subevents[0]->template min<I>();
              const size_t binLast = axis.numBins(); // index of last visible bin
 
              for (const auto& subevt : subevents) {
                if (fidx >= subevt->fills().size()) continue;
                const EdgeT coord = std::get<I>(subevt->fills()[fidx].first);
                size_t idx = axis.index(coord);
                if (coord >= edgeMax) {
                  if (coord > edgeMax) ++over;
                  idx = binLast; // cut off at highest visible bin
                }
                else if (coord < edgeMin) {
                  ++under;
                  idx = 1; // cut off at lowest visible bin
                }
                // find index of closest neighbouring bin
                size_t ibn = idx;
                if (coord > axis.mid(idx)) {
                  if (idx != binLast) ++ibn;
                }
                else {
                  if (idx != 1) --ibn;
                }
 
                // construct rectangular windows of width = 2*delta
                EdgeT lo, hi;
                const EdgeT ibw = axis.width(idx) < axis.width(ibn) ? idx : ibn;
                if (fsmear > 0.0) {
                  const EdgeT delta = 0.5 * fsmear * axis.width(ibw);
                  lo = coord - delta;
                  hi = coord + delta;
                }
                else {
                  const EdgeT delta = 0.5 * axis.width(ibw);
                  if (coord > edgeMax) {
                    lo = max(edgeMax, coord - delta);
                    hi = max(edgeMax + 2 * delta, coord + delta);
                  }
                  else if (coord < edgeMin) {
                    lo = min(edgeMin - 2 * delta, coord - delta);
                    hi = min(edgeMin, coord + delta);
                  }
                  else {
                    lo = axis.min(idx);
                    hi = axis.max(idx);
                  }
                }
                edges.push_back(lo);
                edges.push_back(hi);
              }
              for (size_t i = 0; i < edges.size(); i += 2) {
                const EdgeT wsize = edges[i + 1] - edges[i];
                if (over == nSubevents && edges[i] < edgeMax && edges[i + 1] > edgeMax) {
                  edges[i] = edgeMax;
                  edges[i + 1] = edgeMax + wsize;
                }
                else if (over == 0 && edges[i] < edgeMax && edges[i + 1] > edgeMax) {
                  edges[i] = edgeMax - wsize;
                  edges[i + 1] = edgeMax;
                }
                else if (under == nSubevents && edges[i] < edgeMin && edges[i + 1] > edgeMin) {
                  edges[i] = edgeMin - wsize;
                  edges[i + 1] = edgeMin;
                }
                else if (under == 0 && edges[i] < edgeMin && edges[i + 1] > edgeMin) {
                  edges[i] = edgeMin;
                  edges[i + 1] = edgeMin + wsize;
                }
              }
            } // end of constexpr-check for binned axes
            else {
              // this fill axis is unbinned (e.g. in Profiles)
              for (const auto& subevt : subevents) {
                if (fidx >= subevt->fills().size()) continue;
                const EdgeT coord = std::get<I>(subevt->fills()[fidx].first);
                const EdgeT delta = 0.1 * std::abs(coord);
                edges.push_back(coord - delta);
                edges.push_back(coord + delta);
              }
            }
            // create CAxis with subwindows from the set of window edges
            vector<EdgeT> subwindowEdges = edges;
            std::sort(subwindowEdges.begin(), subwindowEdges.end());
            subwindowEdges.erase(std::unique(subwindowEdges.begin(), subwindowEdges.end()),
                                 subwindowEdges.end());
            subwindows.template axis<I>() = FillAxisT(std::move(subwindowEdges));
          } // end of if constexpr isContinuous
        };
        // execute for each fill dimension
        MetaUtils::staticFor<YAO::FillDimension::value>(constructWindows);
 
        // Subwindows are given by visible bins of the SubwindowT,
        // so we can skip all under-/overflows
        const vector<size_t> overflows = subwindows.calcOverflowBinsIndices();
        const auto& itEnd = overflows.cend();
        size_t nSubwindows = subwindows.numBins();
        for (size_t i = 0; i < nSubwindows; ++i) {
          if (std::find(overflows.cbegin(), itEnd, i) != itEnd) continue;
 
          const auto coords = subwindows.edgeTuple(i);
          const double subwindowArea = subwindows.dVol(i);
          size_t nSubfills = 0;
          double windowFrac = 0.;
          std::valarray<double> sumw(0.0, weights[0].size());
          for (size_t sidx = 0, eidx = 0; sidx < nSubevents; ++sidx) {
            if (fidx >= subevents[sidx]->fills().size()) continue;
            bool pass = true;
            double windowArea = 1.0;
            auto checkSubwindowOverlap = [&](auto I) {
              using isContinuous = typename SubwindowT::template is_CAxis<I>;
              using EdgeT = typename SubwindowT::template getAxisT<I>::EdgeT;
              const EdgeT coord = std::get<I>(coords);
              const vector<EdgeT>& edges = std::get<I>(windows);
              if constexpr (isContinuous::value) {
                pass &= (edges[2 * eidx] <= coord && coord <= edges[2 * eidx + 1]);
                windowArea *= edges[2 * eidx + 1] - edges[2 * eidx];
              }
              else {
                pass &= (coord == edges[eidx]);
              }
            };
            MetaUtils::staticFor<YAO::FillDimension::value>(checkSubwindowOverlap);
            if (pass) {
              windowFrac = subwindowArea / windowArea;
              sumw += subevents[sidx]->fills()[fidx].second * weights[sidx];
              ++nSubfills;
            }
            ++eidx;
          }
          if (nSubfills) {
            const double fillFrac = (double)nSubfills / (double)nSubevents;
            rtn.emplace_back(coords, sumw / fillFrac, fillFrac * windowFrac); // normalise to a single fill
          }
        }
      } // end of loop over fills
      return rtn;
    } // end of applyFillWindows
 
  } // end of anonymous name space
 


  class MultiplexedAO {
  public:
 
    virtual ~MultiplexedAO() { }

    using Inner = YODA::AnalysisObject;

    virtual void newSubEvent() = 0;

    virtual void collapseEventGroup(const vector<std::valarray<double>>& weight,
                                    const double nlowfrac = -1.0) = 0;

    virtual void pushToFinal() = 0;

    virtual YODA::AnalysisObjectPtr activeAO() const = 0;

    virtual string basePath() const = 0;

    virtual YODA::AnalysisObject* operator->() = 0;

    virtual YODA::AnalysisObject* operator->() const = 0;

    virtual const YODA::AnalysisObject& operator*() const = 0;

    virtual void setActiveWeightIdx(size_t iWeight) = 0;

    virtual void setActiveFinalWeightIdx(size_t iWeight) = 0;

    virtual void unsetActiveWeight() = 0;

    virtual void initBootstrap() = 0;

    bool operator==(const MultiplexedAO& p) {
      return (this == &p);
    }

    bool operator!=(const MultiplexedAO& p) {
      return (this != &p);
    }
 
    const vector<bool>& fillOutcomes() const {
      return _fillOutcomes;
    }
 
    const vector<double>& fillFractions() const {
      return _fillFractions;
    }
 
  protected:
 
    vector<bool> _fillOutcomes;
 
    vector<double> _fillFractions;
  };
 

  template <typename T>
  class Multiplexer : public MultiplexedAO {
 
  public:
 
    friend class Analysis;
    //friend class AnalysisHandler;

    using Inner = T;
    using MultiplexedAO::_fillFractions;
    using MultiplexedAO::_fillOutcomes;
 
    Multiplexer() = default;
 
    Multiplexer(const vector<string>& weightNames, const T& p) {
      _basePath = p.path();
      _baseName = p.name();
      for (const string& weightname : weightNames) {
        _persistent.push_back(make_shared<T>(p));
        _final.push_back(make_shared<T>(p));
 
        typename T::Ptr obj = _persistent.back();
        obj->setPath("/RAW" + obj->path());
        typename T::Ptr final = _final.back();
        if (weightname != "") {
          obj->setPath(obj->path() + "[" + weightname + "]");
          final->setPath(final->path() + "[" + weightname + "]");
        }
      }
    }
 
    ~Multiplexer() = default;

    typename T::Ptr active() const {
      if (!_active) {
#ifdef HAVE_BACKTRACE
        void* buffer[4];
        backtrace(buffer, 4);
        backtrace_symbols_fd(buffer, 4, 1);
#endif
        assert(false && "No active pointer set. Was this object booked in init()?");
      }
      return _active;
    }
 
    template <typename U = T>
    auto binning() const -> std::enable_if_t<hasBinning<T>::value, const typename U::BinningT&> {
      return _persistent.back()->binning();
    }

    string basePath() const {
      return _basePath;
    }

    string baseName() const {
      return _baseName;
    }
 

    explicit operator bool() const {
      return static_cast<bool>(_active);
    }

    bool operator!() const {
      return !_active;
    }
 

    T* operator->() {
      return active().get();
    }

    T* operator->() const {
      return active().get();
    }

    T& operator*() {
      return *active();
    }

    const T& operator*() const {
      return *active();
    }
 

    friend bool operator==(const Multiplexer& a, const Multiplexer& b) {
      if (a._persistent.size() != b._persistent.size()) {
        return false;
      }
      // also test for binning compatibility
      for (size_t i = 0; i < a._persistent.size(); ++i) {
        if (a._persistent.at(i) != b._persistent.at(i)) {
          return false;
        }
      }
      return true;
    }

    friend bool operator!=(const Multiplexer& a, const Multiplexer& b) {
      return !(a == b);
    }

    friend bool operator<(const Multiplexer a, const Multiplexer& b) {
      if (a._persistent.size() >= b._persistent.size()) {
        return false;
      }
      for (size_t i = 0; i < a._persistent.size(); ++i) {
        if (*(a._persistent.at(i)) >= *(b._persistent.at(i))) {
          return false;
        }
      }
      return true;
    }
 



    void setActiveWeightIdx(size_t iWeight) {
      _active = _persistent.at(iWeight);
    }

    void setActiveFinalWeightIdx(size_t iWeight) {
      _active = _final.at(iWeight);
    }

    void unsetActiveWeight() {
      _active.reset();
    }

    void reset() {
      active()->reset();
    }
 

    void newSubEvent() {
      _evgroup.emplace_back(new FillCollector<T>(_persistent[0]));
      _active = _evgroup.back();
      assert(_active);
    }

    void collapseEventGroup(const vector<std::valarray<double>>& weights, const double nlowfrac = -1.0) {

      if constexpr (isFillable<T>::value) {
 
        // Should have as many subevent fills as subevent weights
        assert(_evgroup.size() == weights.size());
 
        if (_evgroup.size() == 1) { // Have we had subevents at all?
          // ensure vector length matches size of multiplexed AO
          if (_fillOutcomes.empty()) initBootstrap();
          // reset fill positions
          std::fill(_fillOutcomes.begin(), _fillOutcomes.end(), false);
          std::fill(_fillFractions.begin(), _fillFractions.end(), 0.0);
 
          // Simple replay of all collected fills:
          // each fill is inserted into every persistent AO
          for (auto f : _evgroup[0]->fills()) {
            int pos = -1;
            double frac = 0.0;
            for (size_t m = 0; m < _persistent.size(); ++m) { //< m is the variation index
              if (!m) frac = f.second;
              pos = _persistent[m]->fill(std::move(f.first), std::move(f.second) * weights[0][m]);
            }
            if (pos >= 0) {
              _fillOutcomes[pos] = true;
              _fillFractions[pos] += frac;
            }
          }
        }
        else {
          collapseSubevents(weights, nlowfrac);
        }
      }
 
      _evgroup.clear();
      _active.reset();
    }

    void collapseSubevents(const vector<std::valarray<double>>& weights, const double nlowfrac) {
      if (_evgroup.empty()) return;
      if constexpr (isFillable<T>::value) {
        if constexpr (!std::is_same<T, YODA::Counter>::value) { // binned objects
          for (auto&& fw : mkFillWindows<T>(_evgroup, _persistent[0]->path(), weights, nlowfrac)) {
            for (size_t m = 0; m < _persistent.size(); ++m) {                        // for each multiweight
              _persistent[m]->fill(typename T::FillType(std::move(std::get<0>(fw))), //< coords
                                   std::move(std::get<1>(fw)[m]),                    //< weight
                                   std::move(std::get<2>(fw)));                      //< fraction
            }
          } // end of loop over fill windows
        }
        else {
          for (size_t m = 0; m < _persistent.size(); ++m) { //< m is the variation index
            vector<double> sumfw{0.0};                      // final fill weights (one per fill)
            for (size_t n = 0; n < _evgroup.size(); ++n) {  //< n is the correlated sub-event index
              const auto& fills = _evgroup[n]->fills();
              // resize if this subevent has an
              // even larger number of fills
              if (fills.size() > sumfw.size()) {
                sumfw.resize(fills.size(), 0.0);
              }
              size_t fi = 0;
              for (const auto& f : fills) { // collapse sub-events and aggregate final fill weights
                sumfw[fi++] += f.second * weights[n][m]; // f.second is optional user-supplied scaling
              }
            }
            for (double fw : sumfw) { // fill persistent Counters
              _persistent[m]->fill(std::move(fw));
            }
          } // end of loop over weights
        }
      } // end of isFillable<T>::value
    }

    void pushToFinal() {
      for (size_t m = 0; m < _persistent.size(); ++m) { //< variation weight index
        _final.at(m)->clearAnnotations();               // in case this is a repeat call
        copyAO<T>(_persistent.at(m), _final.at(m));
        // Remove the /RAW prefix, if there is one, from the final copy
        if (_final[m]->path().substr(0, 4) == "/RAW") _final[m]->setPath(_final[m]->path().substr(4));
      }
    }

    const vector<typename T::Ptr>& persistent() const {
      return _persistent;
    }

    typename T::Ptr persistent(const size_t iWeight) {
      return _persistent.at(iWeight);
    }

    const vector<typename T::Ptr>& final() const {
      return _final;
    }

    typename T::Ptr final(const size_t iWeight) {
      return _final.at(iWeight);
    }

    YODA::AnalysisObjectPtr activeAO() const {
      return _active;
    }

    void initBootstrap() {
      if constexpr (isFillable<T>::value) {
        size_t nPos = 1;                                        // Counter only has a single fill position
        if constexpr (!std::is_same<T, YODA::Counter>::value) { // binned objects
          nPos = _persistent.back()->numBins(true, true);
        }
        _fillOutcomes.resize(nPos);
        _fillFractions.resize(nPos);
      }
    }

 
  private:


    vector<typename T::Ptr> _persistent;

    vector<typename T::Ptr> _final;

    vector<typename FillCollector<T>::Ptr> _evgroup;

    typename T::Ptr _active;

    string _basePath;

    string _baseName;

  };
 

  template <typename T>
  class MultiplexPtr {
 
  public:
 
    using value_type = T;
 
    MultiplexPtr() = default;
 
    MultiplexPtr(decltype(nullptr))
        : _p(nullptr) { }

    MultiplexPtr(const vector<string>& weightNames, const typename T::Inner& p)
        : _p(make_shared<T>(weightNames, p)) { }
 
    // Ensure a shared_ptr<T> can be instantiated from a shared_ptr<U>
    template <typename U, typename = decltype(shared_ptr<T>(shared_ptr<U>{}))>
    MultiplexPtr(const shared_ptr<U>& p)
        : _p(p) { }
 
    // Ensure a shared_ptr<T> can be instantiated from a shared_ptr<U>
    template <typename U, typename = decltype(shared_ptr<T>(shared_ptr<U>{}))>
    MultiplexPtr(const MultiplexPtr<U>& p)
        : _p(p.get()) { }

    T& operator->() {
      if (_p == nullptr) {
        throw Error("Dereferencing null AnalysisObject pointer. Is there an unbooked histogram variable?");
      }
      return *_p;
    }
 
    template <typename U = T>
    auto binning() const
        -> std::enable_if_t<hasBinning<typename T::Inner>::value, const typename U::Inner::BinningT&> {
      if (_p == nullptr) {
        throw Error("Dereferencing null AnalysisObject pointer. Is there an unbooked histogram variable?");
      }
      return _p->binning();
    }
 

    const T& operator->() const {
      if (_p == nullptr) {
        throw Error("Dereferencing null AnalysisObject pointer. Is there an unbooked histogram variable?");
      }
      return *_p;
    }

    typename T::Inner& operator*() {
      return **_p;
    }
    const typename T::Inner& operator*() const {
      return **_p;
    }

    explicit operator bool() const {
      return _p && bool(*_p);
    }

    bool operator!() const {
      return !_p || !(*_p);
    }

    bool operator==(const MultiplexPtr& other) const {
      return _p == other._p;
    }

    bool operator!=(const MultiplexPtr& other) const {
      return _p != other._p;
    }

    bool operator<(const MultiplexPtr& other) const {
      return _p < other._p;
    }

    bool operator>(const MultiplexPtr other) const {
      return _p > other._p;
    }

    bool operator<=(const MultiplexPtr& other) const {
      return _p <= other._p;
    }

    bool operator>=(const MultiplexPtr& other) const {
      return _p >= other._p;
    }

    shared_ptr<T> get() const {
      return _p;
    }
 
  private:

    shared_ptr<T> _p;
  };

 

 
  using MultiplexAOPtr = MultiplexPtr<MultiplexedAO>;
 
  template <size_t DbnN, typename... AxisT>
  using BinnedDbnPtr = MultiplexPtr<Multiplexer<YODA::BinnedDbn<DbnN, AxisT...>>>;
 
  template <typename... AxisT>
  using BinnedHistoPtr = BinnedDbnPtr<sizeof...(AxisT), AxisT...>;
 
  template <typename... AxisT>
  using BinnedProfilePtr = BinnedDbnPtr<sizeof...(AxisT) + 1, AxisT...>;
 
  template <typename... AxisT>
  using BinnedEstimatePtr = MultiplexPtr<Multiplexer<YODA::BinnedEstimate<AxisT...>>>;
 
  template <size_t N>
  using ScatterNDPtr = MultiplexPtr<Multiplexer<YODA::ScatterND<N>>>;
 
  using CounterPtr = MultiplexPtr<Multiplexer<YODA::Counter>>;
  using Estimate0DPtr = MultiplexPtr<Multiplexer<YODA::Estimate0D>>;
  using Histo1DPtr = BinnedHistoPtr<double>;
  using Histo2DPtr = BinnedHistoPtr<double, double>;
  using Histo3DPtr = BinnedHistoPtr<double, double, double>;
  using Profile1DPtr = BinnedProfilePtr<double>;
  using Profile2DPtr = BinnedProfilePtr<double, double>;
  using Profile3DPtr = BinnedProfilePtr<double, double, double>;
  using Estimate1DPtr = BinnedEstimatePtr<double>;
  using Estimate2DPtr = BinnedEstimatePtr<double, double>;
  using Estimate3DPtr = BinnedEstimatePtr<double, double, double>;
  using Scatter1DPtr = ScatterNDPtr<1>;
  using Scatter2DPtr = ScatterNDPtr<2>;
  using Scatter3DPtr = ScatterNDPtr<3>;
 
  using YODA::Counter;
  using YODA::Estimate0D;
  using YODA::Estimate1D;
  using YODA::Estimate2D;
  using YODA::Estimate3D;
  using YODA::Histo1D;
  using YODA::Histo2D;
  using YODA::Histo3D;
  using YODA::Point1D;
  using YODA::Point2D;
  using YODA::Point3D;
  using YODA::Profile1D;
  using YODA::Profile2D;
  using YODA::Profile3D;
  using YODA::Scatter1D;
  using YODA::Scatter2D;
  using YODA::Scatter3D;

 

 
  inline bool isTmpPath(const std::string& path, const bool tmp_only = false) {
    if (tmp_only) return path.find("/TMP/") != string::npos;
    return path.find("/TMP/") != string::npos || path.find("/_") != string::npos;
  }

  map<string, YODA::AnalysisObjectPtr> getRefData(const string& papername);


  string getDatafilePath(const string& papername);
 

  template <typename T>
  struct ReferenceTraits { };
 
  template <>
  struct ReferenceTraits<YODA::Counter> {
    using RefT = YODA::Estimate0D;
  };
 
  template <>
  struct ReferenceTraits<YODA::Estimate0D> {
    using RefT = YODA::Estimate0D;
  };
 
  template <typename... AxisT>
  struct ReferenceTraits<YODA::BinnedEstimate<AxisT...>> {
    using RefT = YODA::BinnedEstimate<AxisT...>;
  };
 
  template <size_t DbnN, typename... AxisT>
  struct ReferenceTraits<YODA::BinnedDbn<DbnN, AxisT...>> {
    using RefT = YODA::ScatterND<sizeof...(AxisT) + 1>;
  };
 
  template <size_t N>
  struct ReferenceTraits<YODA::ScatterND<N>> {
    using RefT = YODA::ScatterND<N>;
  };

  template <typename TPtr>
  inline bool bookingCompatible(TPtr a, TPtr b) {
    return *a == *b;
  }
  //
  inline bool bookingCompatible(CounterPtr, CounterPtr) {
    return true;
  }
  //
  inline bool bookingCompatible(YODA::CounterPtr, YODA::CounterPtr) {
    return true;
  }
  //
  inline bool bookingCompatible(Estimate0DPtr, Estimate0DPtr) {
    return true;
  }
  //
  inline bool bookingCompatible(YODA::Estimate0DPtr, YODA::Estimate0DPtr) {
    return true;
  }
  //
  template <size_t N>
  inline bool bookingCompatible(ScatterNDPtr<N> a, ScatterNDPtr<N> b) {
    return a->numPoints() == b->numPoints();
  }
  //
  template <size_t N>
  inline bool bookingCompatible(YODA::ScatterNDPtr<N> a, YODA::ScatterNDPtr<N> b) {
    return a->numPoints() == b->numPoints();
  }
 
  inline bool beamInfoCompatible(YODA::AnalysisObjectPtr a, YODA::AnalysisObjectPtr b) {
    YODA::BinnedEstimatePtr<string> beamsA = std::dynamic_pointer_cast<YODA::BinnedEstimate<string>>(a);
    YODA::BinnedEstimatePtr<string> beamsB = std::dynamic_pointer_cast<YODA::BinnedEstimate<string>>(b);
    return beamsA && beamsB && (*beamsA == *beamsB) && beamsA->numBins() == 2
        && fuzzyEquals(beamsA->bin(1).val(), beamsB->bin(1).val())
        && fuzzyEquals(beamsA->bin(2).val(), beamsB->bin(2).val());
  }

 

  class AOPath {
  public:

    AOPath(string fullpath)
        : _valid(false), _path(fullpath), _raw(false), _tmp(false), _ref(false) {
      _valid = init(fullpath);
    }

    string path() const {
      return _path;
    }

    string analysis() const {
      return _analysis;
    }

    string analysisWithOptions() const {
      return _analysis + _optionstring;
    }

    string name() const {
      return _name;
    }

    string weight() const {
      return _weight;
    }

    string weightComponent() const {
      if (_weight == "") return _weight;
      return "[" + _weight + "]";
    }

    bool isRaw() const {
      return _raw;
    }
 
    // Is this a temporary (filling) object?
    bool isTmp() const {
      return _tmp;
    }

    bool isRef() const {
      return _ref;
    }

    string optionString() const {
      return _optionstring;
    }

    bool hasOptions() const {
      return !_options.empty();
    }

    void removeOption(string opt) {
      _options.erase(opt);
      fixOptionString();
    }

    void setOption(string opt, string val) {
      _options[opt] = val;
      fixOptionString();
    }

    bool hasOption(string opt) const {
      return _options.find(opt) != _options.end();
    }

    string getOption(string opt) const {
      auto it = _options.find(opt);
      if (it != _options.end()) return it->second;
      return "";
    }

    void fixOptionString();

    string mkPath() const;
    string setPath() {
      return _path = mkPath();
    }

    void debug() const;

    bool operator<(const AOPath& other) const {
      return _path < other._path;
    }

    bool valid() const {
      return _valid;
    };
    bool operator!() const {
      return !valid();
    }
 
  private:

    bool init(string fullpath);
    bool chopweight(string& fullpath);
    bool chopoptions(string& anal);
 
    bool _valid;
    string _path;
    string _analysis;
    string _optionstring;
    string _name;
    string _weight;
    bool _raw;
    bool _tmp;
    bool _ref;
    map<string, string> _options;
  };
 
}
 
#endif