Analysis.hh

// -*- C++ -*-
#ifndef RIVET_Analysis_HH
#define RIVET_Analysis_HH
 
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/AnalysisInfo.hh"
#include "Rivet/Event.hh"
#include "Rivet/Projection.hh"
#include "Rivet/ProjectionApplier.hh"
#include "Rivet/ProjectionHandler.hh"
#include "Rivet/AnalysisLoader.hh"
#include "Rivet/Tools/Cuts.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Tools/ParticleUtils.hh"
#ifdef HAVE_H5
#include "Rivet/Tools/RivetHDF5.hh"
#endif
#include "Rivet/Tools/HistoGroup.hh"
#include "Rivet/Tools/RivetMT2.hh"
#include "Rivet/Tools/RivetPaths.hh"
#include "Rivet/Tools/RivetYODA.hh"
#include "Rivet/Tools/Percentile.hh"
#include "Rivet/Tools/Cutflow.hh"
#include "Rivet/Projections/CentralityProjection.hh"
#include "Rivet/Projections/FastJets.hh"
#include <tuple>
 

#define vetoEvent                                                       \
  do { MSG_DEBUG("Vetoing event on line " << __LINE__ << " of " << __FILE__); return; } while(0)
 
namespace Rivet {
 
 
  // Convenience for analysis writers
  using std::cout;
  using std::cerr;
  using std::endl;
  using std::tuple;
  using std::stringstream;
  using std::swap;
  using std::numeric_limits;
 
 
  // Forward declaration
  class AnalysisHandler;
 

  class Analysis : public ProjectionApplier {
  public:

    friend class AnalysisHandler;
 

    Analysis(const std::string& name);

    virtual ~Analysis() {}

    Analysis& operator = (const Analysis&) = delete;
 
 
  public:


    virtual void init() { }

    virtual void analyze(const Event& event) = 0;

    virtual void finalize() { }

 


    virtual void preInit() { }
    virtual void postInit() { }
    virtual void preAnalyze(const Event&) { }
    virtual void postAnalyze(const Event&) { }
    virtual void preFinalize() { }
    virtual void postFinalize() { }

    void syncDeclQueue() {
      this->_syncDeclQueue();
      this->markAsOwned();
    }

 
 
  public:


    void loadInfo() { info().parseInfoFile(); }

    const AnalysisInfo& info() const {
      if (!_info) throw Error("No AnalysisInfo object :-O");
      return *_info;
    }

    virtual std::string name() const {
      return ( (info().name().empty()) ? _defaultname : info().name() ) + _optstring;
    }

    std::string analysisDataPath(const std::string& extn, const std::string& suffix="") {
      string filename = name() + (suffix.empty() ? "" : "-") + suffix + "." + extn;
      return findAnalysisDataFile(filename);
    }

    virtual std::string inspireID() const {
      return info().inspireID();
    }

    virtual std::string spiresID() const {
      return info().spiresID();
    }

    virtual std::vector<std::string> authors() const {
      return info().authors();
    }

    virtual std::string summary() const {
      return info().summary();
    }

    virtual std::string description() const {
      return info().description();
    }

    virtual std::string runInfo() const {
      return info().runInfo();
    }

    virtual std::string experiment() const {
      return info().experiment();
    }

    virtual std::string collider() const {
      return info().collider();
    }

    virtual std::string year() const {
      return info().year();
    }

    virtual double luminosityfb() const {
      return info().luminosityfb();
    }

    virtual double luminosity() const {
      return info().luminosity();
    }

    double luminositypb() const { return luminosity(); }

    virtual std::vector<std::string> references() const {
      return info().references();
    }

    virtual std::string bibKey() const {
      return info().bibKey();
    }

    virtual std::string bibTeX() const {
      return info().bibTeX();
    }

    virtual std::string status() const {
      return (info().status().empty()) ? "UNVALIDATED" : info().status();
    }

    virtual std::string warning() const {
      return info().warning();
    }

    virtual std::vector<std::string> todos() const {
      return info().todos();
    }

    virtual std::vector<std::string> validation() const {
      return info().validation();
    }

    virtual bool reentrant() const {
      return info().reentrant();
    }

    virtual const std::vector<std::string>& keywords() const {
      return info().keywords();
    }

    virtual std::string refMatch() const {
      return info().refMatch();
    }

    virtual std::string refUnmatch() const {
      return info().refUnmatch();
    }

    virtual std::string writerDoublePrecision() const {
      return info().writerDoublePrecision();
    }

    virtual const std::vector<PdgIdPair>& requiredBeamIDs() const {
      return info().beamIDs();
    }

    virtual Analysis& setRequiredBeamIDs(const std::vector<PdgIdPair>& beamids) {
      info().setBeamIDs(beamids);
      return *this;
    }

    virtual const std::vector<std::pair<double, double> >& requiredBeamEnergies() const {
      return info().energies();
    }

    virtual Analysis& setRequiredBeamEnergies(const std::vector<std::pair<double, double> >& energies) {
      info().setEnergies(energies);
      return *this;
    }
 

    virtual std::string refFile() const {
      return info().refFile();
    }

    virtual std::string refDataName() const {
      return (info().getRefDataName().empty()) ? _defaultname : info().getRefDataName();
    }

    virtual void setRefDataName(const std::string& ref_data="") {
      info().setRefDataName(!ref_data.empty() ? ref_data : name());
    }
 

    AnalysisInfo& info() {
      if (!_info) throw Error("No AnalysisInfo object :-O");
      return *_info;
    }

 


    bool merging() const {
      return sqrtS() <= 0.0;
    }

    bool compatibleWithRun() const;

    void raiseBeamErrorIf(const bool condition) const;

 


    const ParticlePair& beams() const;

    PdgIdPair beamIDs() const;

    pair<double,double> beamEnergies() const;

    vector<double> allowedEnergies() const;

    double sqrtS() const;

    bool beamsMatch(const ParticlePair& beams) const;

    bool beamsMatch(PdgId beam1, PdgId beam2, double e1, double e2) const;

    bool beamsMatch(const PdgIdPair& beams, const std::pair<double,double>& energies) const;

    bool beamIDsMatch(PdgId beam1, PdgId beam2) const;

    bool beamIDsMatch(const PdgIdPair& beamids) const;

    bool beamEnergiesMatch(double e1, double e2) const;

    bool beamEnergiesMatch(const std::pair<double,double>& energies) const;

    bool beamEnergyMatch(const std::pair<double,double>& energies) const;

    bool beamEnergyMatch(double sqrts) const;

    bool isCompatibleWithSqrtS(double energy, double tolerance=1e-5) const;

 

    AnalysisHandler& handler() const { return *_analysishandler; }

    const Event& currentEvent() const {
      if (!_currentevent) throw Error("No current event set: did you try to access it in init() or finalize()?");
      return *_currentevent;
    }
 
 
  protected:

    Log& getLog() const;

    double crossSection() const;

    double crossSectionPerEvent() const;

    double crossSectionError() const;

    double crossSectionErrorPerEvent() const;

    size_t numEvents() const;

    double effNumEvents() const;

    double sumW() const;
    double sumOfWeights() const { return sumW(); }

    double sumW2() const;
 
 
  protected:


    const std::string histoDir() const;

    const std::string histoPath(const std::string& hname) const;

    const std::string histoPath(unsigned int datasetID, unsigned int xAxisID, unsigned int yAxisID) const;

    const std::string mkAxisCode(unsigned int datasetID, unsigned int xAxisID, unsigned int yAxisID) const;

 
 
    #ifdef HAVE_H5

    H5::File auxFile() const {
      return H5::readFile(name()+".h5");
    }

    template <typename T>
    bool auxData(const string& dsname, T& rtndata) {
      return H5::readData(name()+".h5", dsname, rtndata);
    }
 
    template <typename T>
    T auxData(const string& dsname) {
      return H5::readData<T>(name()+".h5", dsname);
    }

    #endif
 


    const std::map<std::string, YODA::AnalysisObjectPtr>& refData() const {
      _cacheRefData();
      return _refdata;
    }

    bool hasRefData(const string& hname) const {
      _cacheRefData();
      return _refdata.count(hname);
    }

    bool hasRefData(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
      return hasRefData(mkAxisCode(datasetId, xAxisId, yAxisId));
    }
 

    template <typename T=YODA::Estimate1D>
    const T& refData(const string& hname) const {
      _cacheRefData();
      MSG_TRACE("Using histo bin edges for " << name() << ":" << hname);
      if (!_refdata[hname]) {
        MSG_ERROR("Can't find reference histogram " << hname);
        throw Exception("Reference data " + hname + " not found.");
      }
      try {
        return dynamic_cast<T&>(*_refdata[hname]);
      } catch (...) {
        throw Exception("Expected type " + _refdata[hname]->type()+" for reference data \"" + hname + "\".\n");
      }
    }
 

    template <typename T=YODA::Estimate1D>
    const T& refData(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
      const string hname = mkAxisCode(datasetId, xAxisId, yAxisId);
      return refData<T>(hname);
    }

 


    CounterPtr& book(CounterPtr&, const std::string& name);

    CounterPtr& book(CounterPtr&, unsigned int datasetID, unsigned int xAxisID, unsigned int yAxisID);

 


    Estimate0DPtr& book(Estimate0DPtr&, const std::string& name);

    Estimate0DPtr& book(Estimate0DPtr&, unsigned int datasetID, unsigned int xAxisID, unsigned int yAxisID);
 



    template<size_t DbnN, typename... AxisT, typename = YODA::enable_if_all_CAxisT<AxisT...>>
    BinnedDbnPtr<DbnN, AxisT...>& book(BinnedDbnPtr<DbnN, AxisT...>& ao,
                                       const std::string& name, const std::vector<size_t>& nbins,
                                       const std::vector<std::pair<double,double>>& loUpPairs) {
      if (nbins.size() != loUpPairs.size()) {
        throw RangeError("Vectors should have the same size!");
      }
      const string path = histoPath(name);
 
      YODA::BinnedDbn<DbnN, AxisT...> yao(nbins, loUpPairs, path);
      _setWriterPrecision(path, yao);
 
      return ao = registerAO(yao);
    }
 
    // Specialiation for 1D
    Histo1DPtr& book(Histo1DPtr& ao, const std::string& name,
                                     const size_t nbins, const double lower, const double upper) {
      return book(ao, name, vector<size_t>{nbins},
                  vector<pair<double,double>>{{lower,upper}});
    }
    //
    Profile1DPtr& book(Profile1DPtr& ao, const std::string& name,
                                         const size_t nbins, const double lower, const double upper) {
      return book(ao, name, vector<size_t>{nbins},
                  vector<pair<double,double>>{{lower,upper}});
    }
 
    // Specialiation for 2D
    Histo2DPtr& book(Histo2DPtr& ao, const std::string& name,
                                     const size_t nbinsX, const double lowerX, const double upperX,
                                     const size_t nbinsY, const double lowerY, const double upperY) {
      return book(ao, name, vector<size_t>{nbinsX,nbinsY},
                  vector<pair<double,double>>{{lowerX,upperX}, {lowerY,upperY}});
    }
    //
    Profile2DPtr& book(Profile2DPtr& ao, const std::string& name,
                                         const size_t nbinsX, const double lowerX, const double upperX,
                                         const size_t nbinsY, const double lowerY, const double upperY) {
      return book(ao, name, vector<size_t>{nbinsX,nbinsY},
                  vector<pair<double,double>>{{lowerX,upperX}, {lowerY,upperY}});
    }
 
    // Specialiation for 3D
    Histo3DPtr& book(Histo3DPtr& ao, const std::string& name,
                                     const size_t nbinsX, const double lowerX, const double upperX,
                                     const size_t nbinsY, const double lowerY, const double upperY,
                                     const size_t nbinsZ, const double lowerZ, const double upperZ) {
      return book(ao, name, vector<size_t>{nbinsX,nbinsY,nbinsZ},
                  vector<pair<double,double>>{{lowerX,upperX}, {lowerY,upperY}, {lowerZ,upperZ}});
    }
    //
    Profile3DPtr& book(Profile3DPtr& ao, const std::string& name,
                                         const size_t nbinsX, const double lowerX, const double upperX,
                                         const size_t nbinsY, const double lowerY, const double upperY,
                                         const size_t nbinsZ, const double lowerZ, const double upperZ) {
      return book(ao, name, vector<size_t>{nbinsX,nbinsY,nbinsZ},
                  vector<pair<double,double>>{{lowerX,upperX}, {lowerY,upperY}, {lowerZ,upperZ}});
    }

    template<size_t DbnN, typename... AxisT>
    BinnedDbnPtr<DbnN, AxisT...>& book(BinnedDbnPtr<DbnN, AxisT...>& ao, const std::string& name,
                                       const std::vector<AxisT>&... binedges) {
      const string path = histoPath(name);
      YODA::BinnedDbn<DbnN, AxisT...> yao(binedges..., path);
      _setWriterPrecision(path, yao);
 
      return ao = registerAO(yao);
    }

    template<size_t DbnN, typename... AxisT>
    BinnedDbnPtr<DbnN, AxisT...>& book(BinnedDbnPtr<DbnN, AxisT...>& ao, const std::string& name,
                                       const std::initializer_list<AxisT>&... binedges) {
      return book(ao, name, vector<AxisT>{binedges} ...);
    }

    template<size_t DbnN, typename... AxisT>
    BinnedDbnPtr<DbnN,  AxisT...>& book(BinnedDbnPtr<DbnN, AxisT...>& ao, const std::string& name,
                                        const YODA::BinnedEstimate<AxisT...>& refest) {
      const string path = histoPath(name);
 
      YODA::BinnedDbn<DbnN, AxisT...> yao(refest.binning(), path);
      for (const string& a : yao.annotations()) {
        if (a != "Path")  yao.rmAnnotation(a);
      }
      _setWriterPrecision(path, yao);
      return ao = registerAO(yao);
    }

    template<size_t DbnN, typename... AxisT>
    BinnedDbnPtr<DbnN, AxisT...>& book(BinnedDbnPtr<DbnN, AxisT...>& ao, const std::string& name) {
      return book(ao, name, refData<YODA::BinnedEstimate<AxisT...>>(name));
    }

    template<size_t DbnN, typename... AxisT>
    BinnedDbnPtr<DbnN, AxisT...>& book(BinnedDbnPtr<DbnN, AxisT...>& ao, const unsigned int datasetID,
                                       const unsigned int xAxisID, const unsigned int yAxisID) {
      const string name = mkAxisCode(datasetID, xAxisID, yAxisID);
      return book(ao, name);
    }


 
    template <typename GroupAxisT, typename... AxisT>
    HistoGroupPtr<GroupAxisT, AxisT...>& book(HistoGroupPtr<GroupAxisT, AxisT...>& ao,
                                              const std::vector<GroupAxisT>& edges,
                                              const std::vector<std::string>& names) {
      ao = make_shared<HistoGroup<GroupAxisT, AxisT...>>(edges);
      if (ao->numBins() != names.size()) {
        throw RangeError("Binning and reference-data names don't match!");
      }
      for (auto& b : ao->bins()) {
        const string& refname = names[b.index()-1];
        book(b, refname, refData<YODA::BinnedEstimate<AxisT...>>(refname));
      }
      return ao;
    }
 
    template <typename GroupAxisT, typename... AxisT>
    HistoGroupPtr<GroupAxisT, AxisT...>& book(HistoGroupPtr<GroupAxisT, AxisT...>& ao,
                                              const std::vector<GroupAxisT>& edges) {
      return ao = make_shared<HistoGroup<GroupAxisT, AxisT...>>(edges);
    }
 
    template <typename GroupAxisT, typename... AxisT>
    HistoGroupPtr<GroupAxisT, AxisT...>& book(HistoGroupPtr<GroupAxisT, AxisT...>& ao,
                                              std::initializer_list<GroupAxisT>&& edges) {
      return ao = make_shared<HistoGroup<GroupAxisT, AxisT...>>(std::move(edges));
    }



    template<typename... AxisT, typename = YODA::enable_if_all_CAxisT<AxisT...>>
    BinnedEstimatePtr<AxisT...>& book(BinnedEstimatePtr<AxisT...>& ao,
                                      const std::string& name, const std::vector<size_t>& nbins,
                                      const std::vector<std::pair<double,double>>& loUpPairs) {
      if (nbins.size() != loUpPairs.size()) {
        throw RangeError("Vectors should have the same size!");
      }
      const string path = histoPath(name);
 
      YODA::BinnedEstimate<AxisT...> yao(nbins, loUpPairs, path);
      _setWriterPrecision(path, yao);
 
      return ao = registerAO(yao);
    }
 
 
    // Specialiation for 1D
    Estimate1DPtr& book(Estimate1DPtr& ao, const std::string& name,
                                           const size_t nbins, const double lower, const double upper) {
      return book(ao, name, vector<size_t>{nbins},
                  vector<pair<double,double>>{{lower,upper}});
    }
 
    // Specialiation for 2D
    Estimate2DPtr& book(Estimate2DPtr& ao, const std::string& name,
                                           const size_t nbinsX, const double lowerX, const double upperX,
                                           const size_t nbinsY, const double lowerY, const double upperY) {
      return book(ao, name, vector<size_t>{nbinsX,nbinsY},
                  vector<pair<double,double>>{{lowerX,upperX}, {lowerY,upperY}});
    }
 
    // Specialiation for 3D
    Estimate3DPtr& book(Estimate3DPtr& ao, const std::string& name,
                                           const size_t nbinsX, const double lowerX, const double upperX,
                                           const size_t nbinsY, const double lowerY, const double upperY,
                                           const size_t nbinsZ, const double lowerZ, const double upperZ) {
      return book(ao, name, vector<size_t>{nbinsX,nbinsY,nbinsZ},
                  vector<pair<double,double>>{{lowerX,upperX}, {lowerY,upperY}, {lowerZ,upperZ}});
    }

    template<typename... AxisT>
    BinnedEstimatePtr<AxisT...>& book(BinnedEstimatePtr<AxisT...>& ao, const std::string& name,
                                      const std::vector<AxisT>&... binedges) {
      const string path = histoPath(name);
      YODA::BinnedEstimate<AxisT...> yao(binedges..., path);
      _setWriterPrecision(path, yao);
 
      return ao = registerAO(yao);
    }

    template<typename... AxisT>
    BinnedEstimatePtr<AxisT...>& book(BinnedEstimatePtr<AxisT...>& ao, const std::string& name,
                                      const std::initializer_list<AxisT>&... binedges) {
      return book(ao, name, vector<AxisT>{binedges} ...);
    }

    template<typename... AxisT>
    BinnedEstimatePtr<AxisT...>& book(BinnedEstimatePtr<AxisT...>& ao, const std::string& name) {
 
      const string path = histoPath(name);
      YODA::BinnedEstimate<AxisT...> yao;
      try {
        yao = YODA::BinnedEstimate<AxisT...>(refData<YODA::BinnedEstimate<AxisT...>>(name).binning());
      } catch (...) {
        MSG_DEBUG("Couldn't retrieve reference binning, continue with nullary AO constructor.");
      }
      yao.setPath(path);
      _setWriterPrecision(path, yao);
      return ao = registerAO(yao);
    }

    template<typename... AxisT>
    BinnedEstimatePtr<AxisT...>& book(BinnedEstimatePtr<AxisT...>& ao, const unsigned int datasetID,
                                      const unsigned int xAxisID, const unsigned int yAxisID) {
      const string name = mkAxisCode(datasetID, xAxisID, yAxisID);
      return book(ao, name);
    }

 


    template<size_t N>
    ScatterNDPtr<N>& book(ScatterNDPtr<N>& snd, const string& name, const bool copy_pts = false) {
      const string path = histoPath(name);
      YODA::ScatterND<N> scat(path);
      if (copy_pts) {
        const YODA::ScatterND<N> ref = refData<YODA::EstimateND<N-1>>(name).mkScatter();
        for (YODA::PointND<N> p : ref.points()) {
          p.setVal(N-1, 0.0);
          p.setErr(N-1, 0.0);
          scat.addPoint(std::move(p));
        }
      }
      _setWriterPrecision(path, scat);
      return snd = registerAO(scat);
    }

    template<size_t N>
    ScatterNDPtr<N>& book(ScatterNDPtr<N>& snd, const unsigned int datasetID, const unsigned int xAxisID,
                                                const unsigned int yAxisID, const bool copy_pts = false) {
      const string axisCode = mkAxisCode(datasetID, xAxisID, yAxisID);
      return book(snd, axisCode, copy_pts);
    }

    Scatter2DPtr& book(Scatter2DPtr& snd, const string& name,
                       const size_t npts, const double lower, const double upper) {
      const string path = histoPath(name);
 
      Scatter2D scat(path);
      const double binwidth = (upper-lower)/npts;
      for (size_t pt = 0; pt < npts; ++pt) {
        const double bincentre = lower + (pt + 0.5) * binwidth;
        scat.addPoint(bincentre, 0, binwidth/2.0, 0);
      }
      _setWriterPrecision(path, scat);
 
      return snd = registerAO(scat);
    }
    //
    Scatter3DPtr& book(Scatter3DPtr& snd, const string& name,
                       const size_t nptsX, const double lowerX, const double upperX,
                       const size_t nptsY, const double lowerY, const double upperY) {
      const string path = histoPath(name);
 
      Scatter3D scat(path);
      const double xbinwidth = (upperX-lowerX)/nptsX;
      const double ybinwidth = (upperY-lowerY)/nptsY;
      for (size_t xpt = 0; xpt < nptsX; ++xpt) {
        const double xbincentre = lowerX + (xpt + 0.5) * xbinwidth;
        for (size_t ypt = 0; ypt < nptsY; ++ypt) {
          const double ybincentre = lowerY + (ypt + 0.5) * ybinwidth;
          scat.addPoint(xbincentre, ybincentre, 0, 0.5*xbinwidth, 0.5*ybinwidth, 0);
        }
      }
      _setWriterPrecision(path, scat);
 
      return snd = registerAO(scat);
    }

    Scatter2DPtr& book(Scatter2DPtr& snd, const string& name,
                       const std::vector<double>& binedges) {
      const string path = histoPath(name);
 
      Scatter2D scat(path);
      for (size_t pt = 0; pt < binedges.size()-1; ++pt) {
        const double bincentre = (binedges[pt] + binedges[pt+1]) / 2.0;
        const double binwidth = binedges[pt+1] - binedges[pt];
        scat.addPoint(bincentre, 0, binwidth/2.0, 0);
      }
      _setWriterPrecision(path, scat);
 
      return snd = registerAO(scat);
    }
    //
    Scatter3DPtr& book(Scatter3DPtr& snd, const string& name,
                       const std::vector<double>& binedgesX,
                       const std::vector<double>& binedgesY) {
      const string path = histoPath(name);
 
      Scatter3D scat(path);
      for (size_t xpt = 0; xpt < binedgesX.size()-1; ++xpt) {
        const double xbincentre = (binedgesX[xpt] + binedgesX[xpt+1]) / 2.0;
        const double xbinwidth = binedgesX[xpt+1] - binedgesX[xpt];
        for (size_t ypt = 0; ypt < binedgesY.size()-1; ++ypt) {
          const double ybincentre = (binedgesY[ypt] + binedgesY[ypt+1]) / 2.0;
          const double ybinwidth = binedgesY[ypt+1] - binedgesY[ypt];
          scat.addPoint(xbincentre, ybincentre, 0, 0.5*xbinwidth, 0.5*ybinwidth, 0);
        }
      }
      _setWriterPrecision(path, scat);
 
      return snd = registerAO(scat);
    }

    template<size_t N>
    ScatterNDPtr<N>& book(ScatterNDPtr<N>& snd, const string& name, const YODA::ScatterND<N>& refscatter) {
      const string path = histoPath(name);
 
      YODA::ScatterND<N> scat(refscatter, path);
      for (const string& a : scat.annotations()) {
        if (a != "Path")  scat.rmAnnotation(a);
      }
      _setWriterPrecision(path, scat);
 
      return snd = registerAO(scat);
    }



    CutflowPtr& book(CutflowPtr& ao, const string& name, const std::vector<std::string>& edges) {
      const string path = histoPath(name);
      Cutflow yao(edges, path);
      _setWriterPrecision(path, yao);
      return ao = registerAO(yao);
    }

    CutflowPtr& book(CutflowPtr& ao, const string& name, const std::initializer_list<std::string>& edges) {
      return book(ao, name, vector<std::string>{edges});
    }


 
    CutflowsPtr& book(CutflowsPtr& ao, const std::vector<std::string>& edges,
                                           const std::vector<std::vector<std::string>>& innerEdges) {
      ao = make_shared<Cutflows>(edges);
      if (ao->numBins() !=innerEdges.size()) {
        throw RangeError("Outer and Inner edges don't match");
      }
      for (auto& b : ao->bins()) {
        book(b, b.xEdge(), innerEdges[b.index()-1]);
      }
      return ao;
    }
 
    CutflowsPtr& book(CutflowsPtr& ao, const std::vector<std::string>& edges) {
      return ao = make_shared<Cutflows>(edges);
    }
 
    CutflowsPtr& book(CutflowsPtr& ao, std::initializer_list<std::string>&& edges) {
      return ao = make_shared<Cutflows>(std::move(edges));
    }

 

    virtual void rawHookIn(YODA::AnalysisObjectPtr yao) {
      (void) yao; // suppress unused variable warning
    }

    virtual void rawHookOut(const vector<MultiplexAOPtr>& raos, size_t iW) {
      (void) raos; // suppress unused variable warning
      (void) iW; // suppress unused variable warning
    }
 
 
  public:


    const std::map<std::string,std::string>& options() const {
      return _options;
    }

    std::string getOption(std::string optname, string def="") const {
      if ( _options.find(optname) != _options.end() )
        return _options.find(optname)->second;
      return def;
    }

    std::string getOption(std::string optname, const char* def) {
      return getOption<std::string>(optname, def);
    }

    template<typename T>
    T getOption(std::string optname, T def) const {
      if (_options.find(optname) == _options.end()) return def;
      std::stringstream ss;
      ss.exceptions(std::ios::failbit);
      T ret;
      ss << _options.find(optname)->second;
      try {
        ss >> ret;
      } catch (...) {
        throw ReadError("Could not read user-provided option into requested type");
      }
      return ret;
    }

    // template<>
    // bool getOption<bool>(std::string optname, bool def) const {
    bool getOption(std::string optname, bool def) const {
      if (_options.find(optname) == _options.end()) return def;
      const std::string val = getOption(optname);
      const std::string lval = toLower(val);
      if (lval.empty()) return false;
      if (lval == "true" || lval == "yes" || lval == "on") return true;
      if (lval == "false" || lval == "no" || lval == "off") return false;
      return bool(getOption<int>(optname, 0));
    }

 


    const CentralityProjection&
    declareCentrality(const SingleValueProjection &proj,
                      string calAnaName, string calHistName,
                      const string projName,
              PercentileOrder pctorder=PercentileOrder::DECREASING);
 

    template<typename T>
    Percentile<T> book(const string& projName,
                       const vector<pair<double, double>>& centralityBins,
                       const vector<tuple<size_t, size_t, size_t>>& ref) {
 
      using RefT = typename ReferenceTraits<T>::RefT;
      using WrapT = MultiplexPtr<Multiplexer<T>>;
 
      Percentile<T> pctl(this, projName);
 
      const size_t nCent = centralityBins.size();
      for (size_t iCent = 0; iCent < nCent; ++iCent) {
        const string axisCode = mkAxisCode(std::get<0>(ref[iCent]),
                                           std::get<1>(ref[iCent]),
                                           std::get<2>(ref[iCent]));
        const RefT& refscatter = refData<RefT>(axisCode);
 
        WrapT wtf(_weightNames(), T(refscatter, histoPath(axisCode)));
        wtf = addAnalysisObject(wtf);
 
        CounterPtr cnt(_weightNames(), Counter(histoPath("TMP/COUNTER/" + axisCode)));
        cnt = addAnalysisObject(cnt);
 
        pctl.add(wtf, cnt, centralityBins[iCent]);
      }
      return pctl;
    }
 
 
    // /// @brief Book Percentile Multiplexers around AnalysisObjects.
    // ///
    // /// Based on a previously registered CentralityProjection named @a
    // /// projName book one (or several) AnalysisObject(s) named
    // /// according to @a ref where the x-axis will be filled according
    // /// to the percentile output(s) of the @projName.
    // ///
    // /// @todo Convert to just be called book() cf. others
    // template <class T>
    // PercentileXaxis<T> bookPercentileXaxis(string projName,
    //                                        tuple<int, int, int> ref) {
 
    //   typedef typename ReferenceTraits<T>::RefT RefT;
    //   typedef MultiplexPtr<Multiplexer<T>> WrapT;
 
    //   PercentileXaxis<T> pctl(this, projName);
 
    //   const string axisCode = mkAxisCode(std::get<0>(ref),
    //                                      std::get<1>(ref),
    //                                      std::get<2>(ref));
    //   const RefT & refscatter = refData<RefT>(axisCode);
 
    //   WrapT wtf(_weightNames(), T(refscatter, histoPath(axisCode)));
    //   wtf = addAnalysisObject(wtf);
 
    //   CounterPtr cnt(_weightNames(), Counter());
    //   cnt = addAnalysisObject(cnt);
 
    //   pctl.add(wtf, cnt);
    //   return pctl;
    // }

 
 
  private:
 
    // Functions that have to be defined in the .cc file to avoid circular #includes

    vector<string> _weightNames() const;

    YODA::AnalysisObjectPtr _getPreload (const string& name) const;

    MultiplexAOPtr _getOtherAnalysisObject(const std::string & ananame, const std::string& name);

    void _checkBookInit() const;

    bool _inInit() const;

    bool _inFinalize() const;

    template <typename YODAT>
    void _setWriterPrecision(const string& path, YODAT& yao) {
      const string re = _info->writerDoublePrecision();
      if (re != "") {
        std::smatch match;
        const bool needsDP = std::regex_search(path, match, std::regex(re));
        if (needsDP)  yao.setAnnotation("WriterDoublePrecision", "1");
      }
    }
 
 
  private:

    class CounterAdapter {
    public:
 
      CounterAdapter(double x) : x_(x) {}
 
      CounterAdapter(const YODA::Counter& c) : x_(c.val()) {}
 
      CounterAdapter(const YODA::Estimate& e) : x_(e.val()) {}
 
      CounterAdapter(const YODA::Scatter1D& s) : x_(s.points()[0].x()) {
        if (s.numPoints() != 1)  throw RangeError("Can only scale by a single value.");
      }
 
      operator double() const { return x_; }
 
    private:
      double x_;
 
    };
 
 
  public:
 
    double dbl(double x) { return x; }
    double dbl(const YODA::Counter& c) { return c.val(); }
    double dbl(const YODA::Estimate0D& e) { return e.val(); }
    double dbl(const YODA::Scatter1D& s) {
      if ( s.numPoints() != 1 ) throw RangeError("Only scatter with single value supported.");
      return s.points()[0].x();
    }
 
 
  protected:


    template<typename T>
    void scale(MultiplexPtr<Multiplexer<T>>& ao, CounterAdapter factor) {
      if (!ao) {
        MSG_WARNING("Failed to scale AnalysisObject=NULL in analysis "
                    << name() << " (scale=" << double(factor) << ")");
        return;
      }
      if (std::isnan(double(factor)) || std::isinf(double(factor))) {
        MSG_WARNING("Failed to scale AnalysisObject=" << ao->path() << " in analysis: "
                    << name() << " (invalid scale factor = " << double(factor) << ")");
        factor = 0;
      }
      MSG_TRACE("Scaling AnalysisObject " << ao->path() << " by factor " << double(factor));
      try {
        if constexpr( isFillable<T>::value ) {
          ao->scaleW(factor);
        }
        else {
          ao->scale(factor);
        }
      }
      catch (YODA::Exception& we) {
        MSG_WARNING("Could not scale AnalysisObject " << ao->path());
        return;
      }
    }

    template<typename GroupAxisT, typename... AxisT>
    void scale(HistoGroupPtr<GroupAxisT, AxisT...>& group, CounterAdapter factor) {
      if (!group) {
        MSG_WARNING("Failed to scale AnalysisObject=NULL in analysis "
                    << name() << " (scale=" << double(factor) << ")");
        return;
      }
      if (std::isnan(double(factor)) || std::isinf(double(factor))) {
        MSG_WARNING("Failed to scale histo group in analysis: "
                    << name() << " (invalid scale factor = " << double(factor) << ")");
        factor = 0;
      }
      MSG_TRACE("Scaling histo group by factor " << double(factor));
      try {
        group->scaleW(factor);
      }
      catch (YODA::Exception& we) {
        MSG_WARNING("Could not scale histo group.");
        return;
      }
    }

    template<typename GroupAxisT, typename... AxisT>
    void scale(HistoGroupPtr<GroupAxisT, AxisT...>& group, const vector<double>& factors) {
      if (!group) {
        MSG_WARNING("Failed to scale AnalysisObject=NULL in analysis " << name());
        return;
      }
      if (group->numBins(true) != factors.size()) {
        throw RangeError(name() + ": Number of scale factors does not match group binning");
        return;
      }
      for (auto& b : group->bins(true)) {
        if (!b.get())  continue;
        double factor = factors[b.index()];
        if (std::isnan(factor) || std::isinf(factor)) {
          MSG_WARNING("Failed to scale componment of histo group in analysis: "
                      << name() << " (invalid scale factor = " << factor << ")");
          factor = 0;
        }
        MSG_TRACE("Scaling histo group element by factor " << factor);
        try {
          b->scaleW(factor);
        }
        catch (YODA::Exception& we) {
          MSG_WARNING("Could not scale component of histo group.");
        }
      }
    }

    void scale(CutflowsPtr& group, CounterAdapter factor) {
      if (!group) {
        MSG_WARNING("Failed to scale AnalysisObject=NULL in analysis "
                    << name() << " (scale=" << double(factor) << ")");
        return;
      }
      if (std::isnan(double(factor)) || std::isinf(double(factor))) {
        MSG_WARNING("Failed to scale histo group in analysis: "
                    << name() << " (invalid scale factor = " << double(factor) << ")");
        factor = 0;
      }
      MSG_TRACE("Scaling histo group by factor " << double(factor));
      try {
        group->scale(factor);
      }
      catch (YODA::Exception& we) {
        MSG_WARNING("Could not scale histo group.");
        return;
      }
    }

    template<typename T, typename U>
    void scale(std::map<T, U>& aos, CounterAdapter factor) {
      for (auto& item : aos)  scale(item.second, factor);
    }

    template <typename AORange, typename = std::enable_if_t<YODA::isIterable<AORange>>>
    void scale(AORange& aos, CounterAdapter factor) {
      for (auto& ao : aos)  scale(ao, factor);
    }

    template <typename T>
    void scale(std::initializer_list<T> aos, CounterAdapter factor) {
      for (auto& ao : std::vector<T>{aos})  scale(ao, factor);
    }

    template<typename T, typename U>
    void scale(std::map<T, U>& aos, const vector<double>& factors) {
      for (auto& item : aos)  scale(item.second, factors);
    }

    template <typename AORange, typename = std::enable_if_t<YODA::isIterable<AORange>>>
    void scale(AORange& aos, const vector<double>& factors) {
      for (auto& ao : aos)  scale(ao, factors);
    }

    template <typename T>
    void scale(std::initializer_list<T> aos, const vector<double>& factors) {
      for (auto& ao : std::vector<T>{aos})  scale(ao, factors);
    }

    template<typename GroupAxisT, typename... AxisT>
    void divByGroupWidth(HistoGroupPtr<GroupAxisT, AxisT...>& group) {
      if (!group) {
        MSG_WARNING("Failed to scale HistoGroup=NULL in analysis "
                    << name() << " by group axis width");
        return;
      }
      group->divByGroupWidth();
    }

    template<typename T, typename U>
    void divByGroupWidth(std::map<T, U>& aos) {
      for (auto& item : aos)  divByGroupWidth(item.second);
    }

    template <typename AORange, typename = std::enable_if_t<YODA::isIterable<AORange>>>
    void divByGroupWidth(AORange& aos) {
      for (auto& ao : aos)  divByGroupWidth(ao);
    }

    template <typename T>
    void divByGroupWidth(std::initializer_list<T> aos) {
      for (auto& ao : std::vector<T>{aos})  divByGroupWidth(ao);
    }
 

    template <size_t DbnN, typename... AxisT>
    void normalize(BinnedDbnPtr<DbnN, AxisT...> ao, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
      if (!ao) {
        MSG_WARNING("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << double(norm) << ")");
        return;
      }
      MSG_TRACE("Normalizing histo " << ao->path() << " to " << double(norm));
      try {
        const double hint = ao->integral(includeoverflows);
        if (hint == 0)  MSG_DEBUG("Skipping histo with null area " << ao->path());
        else            ao->normalize(norm, includeoverflows);
      }
      catch (YODA::Exception& we) {
        MSG_WARNING("Could not normalize histo " << ao->path());
        return;
      }
    }

    template <typename GroupAxisT, typename... AxisT>
    void normalize(HistoGroupPtr<GroupAxisT, AxisT...> group, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
      if (!group) {
        MSG_WARNING("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << double(norm) << ")");
        return;
      }
      MSG_TRACE("Normalizing histo group  to " << double(norm));
      try {
        const double hint = group->integral(includeoverflows);
        if (hint == 0)  MSG_DEBUG("Skipping histo group with null area.");
        else            group->normalize(norm, includeoverflows);
      }
      catch (YODA::Exception& we) {
        MSG_WARNING("Could not normalize histo group.");
        return;
      }
    }
 

    template <typename AORange, typename = std::enable_if_t<YODA::isIterable<AORange>>>
    void normalize(AORange& aos, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
      for (auto& ao : aos)  normalize(ao, norm, includeoverflows);
    }

    template<typename T>
    void normalize(std::initializer_list<T>&& aos, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
      for (auto& ao : aos)  normalize(ao, norm, includeoverflows);
    }

    template<typename T, typename U>
    void normalize(std::map<T, U>& aos, //BinnedDbnPtr<DbnN, AxisT...>>& aos,
                   const CounterAdapter norm=1.0, const bool includeoverflows=true) {
      for (auto& item : aos)  normalize(item.second, norm, includeoverflows);
    }

    template <typename GroupAxisT, typename... AxisT>
    void normalizeGroup(HistoGroupPtr<GroupAxisT, AxisT...> group, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
      if (!group) {
        MSG_WARNING("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << double(norm) << ")");
        return;
      }
      MSG_TRACE("Normalizing histo group  to " << double(norm));
      try {
        const double hint = group->integral(includeoverflows);
        if (hint == 0)  MSG_DEBUG("Skipping histo group with null area.");
        else            group->normalizeGroup(norm, includeoverflows);
      }
      catch (YODA::Exception& we) {
        MSG_WARNING("Could not normalize histo group.");
        return;
      }
    }

    template <typename AORange, typename = std::enable_if_t<YODA::isIterable<AORange>>>
    void normalizeGroup(AORange& aos, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
      for (auto& ao : aos)  normalizeGroup(ao, norm, includeoverflows);
    }

    template<typename T>
    void normalizeGroup(std::initializer_list<T>&& aos, const CounterAdapter norm=1.0, const bool includeoverflows=true) {
      for (auto& ao : aos)  normalizeGroup(ao, norm, includeoverflows);
    }

    template<typename T, typename U>
    void normalizeGroup(std::map<T, U>& aos, //BinnedDbnPtr<DbnN, AxisT...>>& aos,
                   const CounterAdapter norm=1.0, const bool includeoverflows=true) {
      for (auto& item : aos)  normalizeGroup(item.second, norm, includeoverflows);
    }
 

    template<size_t DbnN, typename... AxisT>
    void barchart(BinnedDbnPtr<DbnN, AxisT...> ao, BinnedEstimatePtr<AxisT...> est) const {
      const string path = est->path();
      *est = ao->mkEstimate(path, "stats", false); //< do NOT divide by bin area cf. a differential dsigma/dX histogram
    }

    void divide(CounterPtr c1, CounterPtr c2, Estimate0DPtr est) const;

    void divide(const YODA::Counter& c1, const YODA::Counter& c2, Estimate0DPtr est) const;

    void divide(Estimate0DPtr e1, Estimate0DPtr e2, Estimate0DPtr est) const;

    void divide(const YODA::Estimate0D& e1, const YODA::Estimate0D& e2, Estimate0DPtr est) const;
 

    template<size_t DbnN, typename... AxisT>
    void divide(const YODA::BinnedDbn<DbnN, AxisT...>& h1, const YODA::BinnedDbn<DbnN, AxisT...>& h2,
                BinnedEstimatePtr<AxisT...> est) const {
      const string path = est->path();
      *est = h1 / h2;
      est->setPath(path);
    }
    //
    template<size_t DbnN, typename... AxisT>
    void divide(BinnedDbnPtr<DbnN, AxisT...> h1, BinnedDbnPtr<DbnN, AxisT...> h2,
                             BinnedEstimatePtr<AxisT...> est) const {
      return divide(*h1, *h2, est);
    }

    template<typename... AxisT>
    void divide(const YODA::BinnedEstimate<AxisT...>& e1, const YODA::BinnedEstimate<AxisT...>& e2,
                BinnedEstimatePtr<AxisT...> est) const {
      const string path = est->path();
      *est = e1 / e2;
      est->setPath(path);
    }
    //
    template<typename... AxisT>
    void divide(BinnedEstimatePtr<AxisT...> e1, BinnedEstimatePtr<AxisT...> e2,
                             BinnedEstimatePtr<AxisT...> est) const {
      return divide(*e1, *e2, est);
    }
 
 

    void efficiency(CounterPtr c1, CounterPtr c2, Estimate0DPtr est) const {
      efficiency(*c1, *c2, est);
    }

    void efficiency(const YODA::Counter& c1, const YODA::Counter& c2, Estimate0DPtr est) const {
      const string path = est->path();
      *est = YODA::efficiency(c1, c2);
      est->setPath(path);
    }
 
 

    template<size_t DbnN, typename... AxisT>
    void efficiency(const YODA::BinnedDbn<DbnN, AxisT...>& h1, const YODA::BinnedDbn<DbnN, AxisT...>& h2,
                    BinnedEstimatePtr<AxisT...> est) const {
      const string path = est->path();
      *est = YODA::efficiency(h1, h2);
      est->setPath(path);
    }
    //
    template<size_t DbnN, typename... AxisT>
    void efficiency(BinnedDbnPtr<DbnN, AxisT...> h1, BinnedDbnPtr<DbnN, AxisT...> h2,
                    BinnedEstimatePtr<AxisT...> est) const {
      efficiency(*h1, *h2, est);
    }
 

    template<typename... AxisT>
    void efficiency(const YODA::BinnedEstimate<AxisT...>& e1, const YODA::BinnedEstimate<AxisT...>& e2,
                    BinnedEstimatePtr<AxisT...> est) const {
      const string path = est->path();
      *est = YODA::efficiency(e1, e2);
      est->setPath(path);
    }
    //
    template<typename... AxisT>
    void efficiency(BinnedEstimatePtr<AxisT...> e1, BinnedEstimatePtr<AxisT...> e2,
                    BinnedEstimatePtr<AxisT...> est) const {
      efficiency(*e1, *e2, est);
    }
 

    template<size_t DbnN, typename... AxisT>
    void asymm(const YODA::BinnedDbn<DbnN, AxisT...>& h1, const YODA::BinnedDbn<DbnN, AxisT...>& h2,
               BinnedEstimatePtr<AxisT...> est) const {
      const string path = est->path();
      *est = YODA::asymm(h1, h2);
      est->setPath(path);
    }
    //
    template<size_t DbnN, typename... AxisT>
    void asymm(BinnedDbnPtr<DbnN, AxisT...> h1, BinnedDbnPtr<DbnN, AxisT...> h2,
               BinnedEstimatePtr<AxisT...> est) const {
      asymm(*h1, *h2, est);
    }

    template<typename... AxisT>
    void asymm(const YODA::BinnedEstimate<AxisT...>& e1, const YODA::BinnedEstimate<AxisT...>& e2,
               BinnedEstimatePtr<AxisT...> est) const {
      const string path = est->path();
      *est = YODA::asymm(e1, e2);
      est->setPath(path);
    }
    //
    template<typename... AxisT>
    void asymm(BinnedEstimatePtr<AxisT...> e1, BinnedEstimatePtr<AxisT...> e2,
               BinnedEstimatePtr<AxisT...> est) const {
      asymm(*e1, *e2, est);
    }

    template<size_t DbnN, typename... AxisT>
    void integrate(const YODA::BinnedDbn<DbnN, AxisT...>& h, BinnedEstimatePtr<AxisT...> est) const {
      const string path = est->path();
      *est = mkIntegral(h);
      est->setPath(path);
    }
    //
    template<size_t DbnN, typename... AxisT>
    void integrate(BinnedDbnPtr<DbnN, AxisT...>& h, BinnedEstimatePtr<AxisT...> est) const {
      integrate(*h, est);
    }

 
 
  public:

    const vector<MultiplexAOPtr>& analysisObjects() const {
      return _analysisobjects;
    }
 
 
  protected:


    size_t defaultWeightIndex() const;

    template <typename YODAT>
    shared_ptr<YODAT> getPreload(const string& path) const {
      return dynamic_pointer_cast<YODAT>(_getPreload(path));
    }
 

    template <typename YODAT>
    MultiplexPtr< Multiplexer<YODAT> > registerAO(const YODAT& yao) {
      using MultiplexerT = Multiplexer<YODAT>;
      using YODAPtrT = shared_ptr<YODAT>;
      using RAOT = MultiplexPtr<MultiplexerT>;
 
      if ( !_inInit() && !_inFinalize() ) {
        MSG_ERROR("Can't book objects outside of init() or finalize()");
        throw UserError(name() + ": Can't book objects outside of init() or finalize().");
      }
 
      // First check that we haven't booked this before.
      // This is allowed when booking in finalize: just warn in that case.
      // If in init(), throw an exception: it's 99.9% never going to be intentional.
      for (auto& waold : analysisObjects()) {
        if ( yao.path() == waold.get()->basePath() ) {
          const string msg = "Found double-booking of " + yao.path() + " in " + name();
          if ( _inInit() ) {
            MSG_ERROR(msg);
            throw LookupError(msg);
          } else {
            MSG_WARNING(msg + ". Keeping previous booking");
          }
          return RAOT(dynamic_pointer_cast<MultiplexerT>(waold.get()));
        }
      }
 
      shared_ptr<MultiplexerT> wao = make_shared<MultiplexerT>();
      wao->_basePath = yao.path();
      YODAPtrT yaop = make_shared<YODAT>(yao);
 
      for (const string& weightname : _weightNames()) {
        // Create two YODA objects for each weight. Copy from
        // preloaded YODAs if present. First the finalized yoda:
        string finalpath = yao.path();
        if ( weightname != "" ) finalpath +=  "[" + weightname + "]";
        YODAPtrT preload = getPreload<YODAT>(finalpath);
        if ( preload ) {
          if ( !bookingCompatible(preload, yaop) ) {
            MSG_WARNING("Found incompatible pre-existing data object with same base path "
                        << finalpath <<  " for " << name());
            preload = nullptr;
          } else {
            MSG_TRACE("Using preloaded " << finalpath << " in " <<name());
            wao->_final.push_back(make_shared<YODAT>(*preload));
          }
        }
        else {
          wao->_final.push_back(make_shared<YODAT>(yao));
          wao->_final.back()->setPath(finalpath);
        }
 
        // Then the raw filling yodas.
        string rawpath = "/RAW" + finalpath;
        preload = getPreload<YODAT>(rawpath);
        if ( preload ) {
          if ( !bookingCompatible(preload, yaop) ) {
            MSG_WARNING("Found incompatible pre-existing data object with same base path "
                        << rawpath <<  " for " << name());
            preload = nullptr;
          } else {
            MSG_TRACE("Using preloaded " << rawpath << " in " <<name());
            wao->_persistent.push_back(make_shared<YODAT>(*preload));
          }
        }
        else {
          wao->_persistent.push_back(make_shared<YODAT>(yao));
          wao->_persistent.back()->setPath(rawpath);
        }
      }
      MultiplexPtr<MultiplexerT> ret(wao);
 
      ret.get()->unsetActiveWeight();
      if ( _inFinalize() ) {
        // If booked in finalize() we assume it is the first time
        // finalize is run.
        ret.get()->pushToFinal();
        ret.get()->setActiveFinalWeightIdx(0);
      }
      _analysisobjects.push_back(ret);
 
      return ret;
    }
 

    template <typename AO=MultiplexAOPtr>
    AO addAnalysisObject(const AO& aonew) {
      _checkBookInit();
 
      for (const MultiplexAOPtr& ao : analysisObjects()) {
 
        // Check AO base-name first
        ao.get()->setActiveWeightIdx(defaultWeightIndex());
        aonew.get()->setActiveWeightIdx(defaultWeightIndex());
        if (ao->path() != aonew->path()) continue;
 
        // If base-name matches, check compatibility
        // NB. This evil is because dynamic_ptr_cast can't work on MultiplexPtr directly
        AO aoold = AO(dynamic_pointer_cast<typename AO::value_type>(ao.get())); //< OMG
        if ( !aoold || !bookingCompatible(aonew, aoold) ) {
          MSG_WARNING("Found incompatible pre-existing data object with same base path "
                      << aonew->path() <<  " for " << name());
          throw LookupError("Found incompatible pre-existing data object with same base path during AO booking");
        }
 
        // Finally, check all weight variations
        for (size_t weightIdx = 0; weightIdx < _weightNames().size(); ++weightIdx) {
          aoold.get()->setActiveWeightIdx(weightIdx);
          aonew.get()->setActiveWeightIdx(weightIdx);
          if (aoold->path() != aonew->path()) {
            MSG_WARNING("Found incompatible pre-existing data object with different weight-path "
                        << aonew->path() <<  " for " << name());
            throw LookupError("Found incompatible pre-existing data object with same weight-path during AO booking");
          }
        }
 
        // They're fully compatible: bind and return
        aoold.get()->unsetActiveWeight();
        MSG_TRACE("Bound pre-existing data object " << aoold->path() <<  " for " << name());
        return aoold;
      }
 
      // No equivalent found
      MSG_TRACE("Registered " << aonew->annotation("Type") << " " << aonew->path() <<  " for " << name());
      aonew.get()->unsetActiveWeight();
 
      _analysisobjects.push_back(aonew);
      return aonew;
    }

    void removeAnalysisObject(const std::string& path);

    void removeAnalysisObject(const MultiplexAOPtr& ao);

    template <typename AO=MultiplexAOPtr>
    const AO getAnalysisObject(const std::string& aoname) const {
      for (const MultiplexAOPtr& ao : analysisObjects()) {
        ao.get()->setActiveWeightIdx(defaultWeightIndex());
        if (ao->path() == histoPath(aoname)) {
          // return dynamic_pointer_cast<AO>(ao);
          return AO(dynamic_pointer_cast<typename AO::value_type>(ao.get()));
        }
      }
      throw LookupError("Data object " + histoPath(aoname) + " not found");
    }
 
 
    // /// Get a data object from the histogram system
    // template <typename AO=YODA::AnalysisObject>
    // const std::shared_ptr<AO> getAnalysisObject(const std::string& name) const {
    //   foreach (const AnalysisObjectPtr& ao, analysisObjects()) {
    //     if (ao->path() == histoPath(name)) return dynamic_pointer_cast<AO>(ao);
    //   }
    //   throw LookupError("Data object " + histoPath(name) + " not found");
    // }
 
    // /// Get a data object from the histogram system (non-const)
    // template <typename AO=YODA::AnalysisObject>
    // std::shared_ptr<AO> getAnalysisObject(const std::string& name) {
    //   foreach (const AnalysisObjectPtr& ao, analysisObjects()) {
    //     if (ao->path() == histoPath(name)) return dynamic_pointer_cast<AO>(ao);
    //   }
    //   throw LookupError("Data object " + histoPath(name) + " not found");
    // }
 

    template <typename AO=MultiplexAOPtr>
    AO getAnalysisObject(const std::string& ananame,
                         const std::string& aoname) {
      MultiplexAOPtr ao = _getOtherAnalysisObject(ananame, aoname);
      // return dynamic_pointer_cast<AO>(ao);
      return AO(dynamic_pointer_cast<typename AO::value_type>(ao.get()));
    }



    template <
      typename... Args, typename CONTAINER,
      typename = std::enable_if_t<
        is_citerable_v<CONTAINER>,
        Jet
      >
    >
    static CONTAINER reclusterJets(const CONTAINER &jetsIn, Args&&... args){
      return FastJets::reclusterJets(jetsIn, std::forward<Args>(args)...);
    }
 
    template <typename T, typename U, typename... Args>
    static std::map<T, U> reclusterJets(const std::map<T, U> &jetsMap, Args&&... args){
      return FastJets::reclusterJets(jetsMap, std::forward<Args>(args)...);
    }
 
    template <
      JetAlg JETALG, typename... Args, typename CONTAINER,
      typename = std::enable_if_t<
        is_citerable_v<CONTAINER>,
        Jet
      >
    >
    static CONTAINER reclusterJets(const CONTAINER &jetsIn, Args&&... args){
      return FastJets::reclusterJets<JETALG>(jetsIn, std::forward<Args>(args)...);
    }
 
    template <JetAlg JETALG, typename T, typename U, typename... Args>
    static std::map<T, U> reclusterJets(const std::map<T, U> &jetsMap, Args&&... args){
      return FastJets::reclusterJets<JETALG>(jetsMap, std::forward<Args>(args)...);
    }

 
 
  private:

    string _defaultname;

    unique_ptr<AnalysisInfo> _info;

    vector<MultiplexAOPtr> _analysisobjects;

    double _crossSection;
    bool _gotCrossSection;

    AnalysisHandler* _analysishandler;

    const Event* _currentevent = nullptr;

    mutable std::map<std::string, YODA::AnalysisObjectPtr> _refdata;

    map<string, string> _options;

    string _optstring;
 
 
  private:


    void _cacheRefData() const;

 
  };
 
 
  // // Template specialisation for literal character strings (which don't play well with stringstream)
  // template<>
  // inline std::string Analysis::getOption(std::string optname, const char* def) {
  //   return getOption<std::string>(optname, def); //.c_str();
  // }
 
 
}
 
 
// Include definition of analysis plugin system so that analyses automatically see it when including Analysis.hh
#include "Rivet/AnalysisBuilder.hh"
 


#define RIVET_DECLARE_PLUGIN(clsname) ::Rivet::AnalysisBuilder<clsname> plugin_ ## clsname

#define RIVET_DECLARE_ALIASED_PLUGIN(clsname, alias) RIVET_DECLARE_PLUGIN(clsname)( #alias )

#define RIVET_DEFAULT_ANALYSIS_CTOR(clsname) clsname() : Analysis(# clsname) {}

#define RIVET_REGISTER_TYPE(...) handler().registerType<__VA_ARGS__>()

#define RIVET_REGISTER_BINNED_SET(...) { \
    RIVET_REGISTER_TYPE(YODA::BinnedHisto<__VA_ARGS__>); \
    RIVET_REGISTER_TYPE(YODA::BinnedProfile<__VA_ARGS__>); \
    RIVET_REGISTER_TYPE(YODA::BinnedEstimate<__VA_ARGS__>); }

 
 
#endif