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 {
 
    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::_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=-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::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