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/Counter.h"
#include "YODA/Histo.h"
#include "YODA/Profile.h"
#include "YODA/Estimate0D.h"
#include "YODA/BinnedEstimate.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 {
 
    using namespace std;
 
    template<typename... Args>
    double distance(const tuple<Args...>& a, const 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(get<I>(a) - get<I>(b));
        }
      };
      MetaUtils::staticFor<sizeof...(Args)>(calculateDistance);
      return rtn;
    }
 
    template <typename T>
    vector<Fills<T>> applyEmptyFillPaddingAndTranspose(const vector<typename FillCollector<T>::Ptr>& subevents) {
 
      static Fill<T> EmptyFill = { typename T::FillType{}, 0.0 };
 
      vector<Fills<T>> matched;
      matched.reserve(subevents.size());
      // First just copy subevents into vectors and find the longest vector.
      size_t maxFillLen = 0; // length of biggest vector
      size_t maxFillPos = 0; // index position of biggest vector
      for (const auto& subevt : subevents) {
        auto&& fills = subevt->fills();
        if (fills.size() > maxFillLen) {
          maxFillLen = fills.size();
          maxFillPos = matched.size();
        }
        matched.push_back(std::move(fills));
      }
 
      // Now, go through all subevents with missing fills.
      const Fills<T>& maxFill = matched[maxFillPos]; // the longest fill
      for (auto& fills : matched) {
 
        if (fills.size() == maxFillLen)  continue;
 
        while (fills.size() < maxFillLen) {
          fills.push_back(EmptyFill);
        }
 
        for (int i = maxFillLen - 1; i >= 0; --i) {
          if (fills[i] == EmptyFill)  continue;
          int j = i;
          while (j+1 < (int)maxFillLen &&
                 fills[j + 1] == EmptyFill &&
                 distance(fills[j].first,
                          maxFill[j].first) >
                 distance(fills[j].first,
                          maxFill[j+1].first)) {
            swap(fills[j], fills[j+1]);
            ++j;
          }
        }
      }
      // finally, transpose and we're done
      vector<Fills<T>> result(maxFillLen, Fills<T>(matched.size()));
      for (size_t i = 0; i < matched.size(); ++i) {
        for (size_t j = 0; j < maxFillLen; ++j) {
          result.at(j).at(i) = std::move(matched.at(i).at(j));
        }
      }
      return result;
    }
 
    template<typename T>
    struct SubwindowType;
    //
    template<typename... Args>
    struct SubwindowType<tuple<Args...>> {
      using type = YODA::Binning<std::decay_t<decltype(std::declval<YODA::Axis<Args>>())>...>;
    };
 
    template<typename YAO>
    vector<std::tuple<typename YAO::FillType, valarray<double>, double>>
    applyFillWindows(shared_ptr<YAO> ao, const Fills<YAO>& subevents,
                     const vector<valarray<double>>& weights, const double fsmear=0.5) {
 
 
      using SubwindowT = typename SubwindowType<typename YAO::FillType>::type;
      SubwindowT subwindows;
 
      const size_t nFills = subevents.size();
      vector<vector<double>> edgesLo, edgesHi;
      edgesLo.resize(YAO::FillDimension::value);
      edgesHi.resize(YAO::FillDimension::value);
 
      auto constructWindows = [&](auto I) {
 
        using FillAxisT = typename SubwindowT::template getAxisT<I>;
        using isContinuous = typename SubwindowT::template is_CAxis<I>;
 
        if constexpr(!isContinuous::value) { // discrete axes don't need smearing
          // use single-edge axis for discrete coordinates
          subwindows.template axis<I>() = FillAxisT({ std::get<I>(subevents[0].first) });
          return;
        }
        else { // continuous axes need windowing
          edgesHi[I].resize(nFills);
          edgesLo[I].resize(nFills);
 
          if constexpr(I < YAO::BinningT::Dimension::value) {
            // this fill axis is binned: window sizes will depend on bin width
            const auto& axis = ao->binning().template axis<I>();
            size_t over = 0, under = 0;
            const double edgeMax = ao->template max<I>();
            const double edgeMin = ao->template min<I>();
            const size_t binLast = axis.numBins(); // index of last visible bin
 
            for (size_t i = 0; i < nFills; ++i) {
              const double edge = get<I>(subevents[i].first);
              size_t idx = axis.index(edge);
              if (edge >= edgeMax) {
                if (edge > edgeMax) ++over;
                idx = binLast; // cut off at highest visible bin
              }
              else if (edge < edgeMin) {
                ++under;
                idx = 1; // cut off at lowest visible bin
              }
              // find index of closest neighbouring bin
              size_t ibn = idx;
              if (edge > axis.mid(idx)) {
                if (idx != binLast) ++ibn;
              }
              else {
                if (idx != 1) --ibn;
              }
 
              // construct rectangular windows of with = 2*delta
              const double ibw = axis.width(idx) < axis.width(ibn)? idx : ibn;
              if ( fsmear > 0.0 ) {
                const double delta = 0.5*fsmear*axis.width(ibw);
                edgesHi[I][i] = edge + delta;
                edgesLo[I][i] = edge - delta;
              }
              else {
                const double delta = 0.5*axis.width(ibw);
                if (edge > edgeMax) {
                  edgesHi[I][i] = max(edgeMax + 2*delta, edge + delta);
                  edgesLo[I][i] = max(edgeMax, edge - delta);
                }
                else if (edge < edgeMin) {
                  edgesHi[I][i] = min(edgeMin, edge + delta);
                  edgesLo[I][i] = min(edgeMin - 2*delta, edge - delta);
                }
                else {
                  edgesHi[I][i] = axis.max(idx);
                  edgesLo[I][i] = axis.min(idx);
                }
              }
            }
            for (size_t i = 0; i < nFills; ++i) {
              const double wsize = edgesHi[I][i] - edgesLo[I][i];
              if (over == nFills && edgesLo[I][i] < edgeMax && edgesHi[I][i] > edgeMax) {
                edgesHi[I][i] = edgeMax + wsize;
                edgesLo[I][i] = edgeMax;
              }
              else if (over == 0 && edgesLo[I][i] < edgeMax && edgesHi[I][i] > edgeMax) {
                edgesLo[I][i] = edgeMax - wsize;
                edgesHi[I][i] = edgeMax;
              }
              else if (under == nFills && edgesLo[I][i] < edgeMin && edgesHi[I][i] > edgeMin) {
                edgesLo[I][i] = edgeMin - wsize;
                edgesHi[I][i] = edgeMin;
              }
              else if (under == 0 && edgesLo[I][i] < edgeMin && edgesHi[I][i] > edgeMin) {
                edgesHi[I][i] = edgeMin + wsize;
                edgesLo[I][i] = edgeMin;
              }
            }
          } // end of constexpr-check for binned axes
          else {
            // this fill axis is unbinned (e.g. in Profiles)
            for (size_t i = 0; i < nFills; ++i) {
              const double edge = get<I>(subevents[i].first);
              const double delta = 0.1*fabs(edge);
              edgesHi[I][i] = edge + delta;
              edgesLo[I][i] = edge - delta;
            }
          }
          // create CAxis with subwindows from the set of window edges
          vector<double> windowEdges;
          std::copy(edgesLo[I].begin(), edgesLo[I].end(), std::back_inserter(windowEdges));
          std::copy(edgesHi[I].begin(), edgesHi[I].end(), std::back_inserter(windowEdges));
          std::sort(windowEdges.begin(), windowEdges.end());
          windowEdges.erase( std::unique(windowEdges.begin(), windowEdges.end()), windowEdges.end() );
          subwindows.template axis<I>() = FillAxisT(windowEdges);
        } // end of if constexpr isContinuous
      };
      // execute for each fill dimension
      MetaUtils::staticFor<YAO::FillDimension::value>(constructWindows);
 
      // Placeholder for the return type: a tuple of {FillType, multi-weights, fill fraction}
      vector<std::tuple<typename YAO::FillType, valarray<double>, double>> rtn;
 
      // Subwindows are given by visible bins of the
      // SubwindowT, so skip all under-/overflows
      const vector<size_t> overflows = subwindows.calcOverflowBinsIndices();
      const auto& itEnd = overflows.cend();
      for (size_t i = 0; i < subwindows.numBins(); ++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.;
        valarray<double> sumw(0.0, weights[0].size()); // one per multiweight
        for (size_t j = 0; j < nFills; ++j) {
          bool pass = true;
          double windowArea = 1.0;
          auto checkSubwindowOverlap = [&](auto I) {
           using isContinuous = typename SubwindowT::template is_CAxis<I>;
            if constexpr (isContinuous::value) {
              const double edge = std::get<I>(coords);
              pass &= (edgesLo[I][j] <= edge && edge <= edgesHi[I][j]);
              windowArea *= edgesHi[I][j] - edgesLo[I][j];
            }
          };
          MetaUtils::staticFor<YAO::FillDimension::value>(checkSubwindowOverlap);
          if (pass) {
            windowFrac = subwindowArea/windowArea;
            sumw += subevents[j].second * weights[j];
            ++nSubfills;
          }
        }
        if (nSubfills) {
          const double fillFrac = (double)nSubfills/(double)nFills;
          rtn.emplace_back(coords, sumw/fillFrac, fillFrac*windowFrac); // normalise to a single fill
        }
      }
      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=0.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::_fillOutcomes;
    using MultiplexedAO::_fillFractions;
 
    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=0.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 constexpr( isFillable<T>::value ) {
        if constexpr (!std::is_same<T, YODA::Counter>::value ) { // binned objects
          for (const Fills<T>& subEvents : applyEmptyFillPaddingAndTranspose<T>(_evgroup)) {
            // construct fill windows with fractional fills
            for (const auto& f : applyFillWindows(_persistent[0], subEvents, weights, nlowfrac)) {
              for ( size_t m = 0; m < _persistent.size(); ++m ) { // for each multiweight
                _persistent[m]->fill( typename T::FillType(get<0>(f)), get<1>(f)[m], get<2>(f) );
              }
            } // end of loop over fill windows
          } // end of loop over subevents
        }
        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));
            }
          }
        }
      }
    }
 
    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::Histo1D;
  using YODA::Histo2D;
  using YODA::Histo3D;
  using YODA::Profile1D;
  using YODA::Profile2D;
  using YODA::Profile3D;
  using YODA::Estimate1D;
  using YODA::Estimate2D;
  using YODA::Estimate3D;
  using YODA::Scatter1D;
  using YODA::Scatter2D;
  using YODA::Scatter3D;
  using YODA::Point1D;
  using YODA::Point2D;
  using YODA::Point3D;
 
 
 
 
  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