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"
#include "Rivet/Tools/BinnedHistogram.hh"
#include "Rivet/Tools/RivetMT2.hh"
#include "Rivet/Tools/RivetYODA.hh"
#include "Rivet/Tools/Percentile.hh"
#include "Rivet/Projections/CentralityProjection.hh"


#define vetoEvent                                                       \
  do { MSG_DEBUG("Vetoing event on line " << __LINE__ << " of " << __FILE__); return; } while(0)


namespace Rivet {


  // Forward declaration
  class AnalysisHandler;

  class Analysis : public ProjectionApplier {

    friend class AnalysisHandler;

  public:



    // /// The default constructor.
    // Analysis();

    Analysis(const std::string& name);

    virtual ~Analysis() {}



  public:



    virtual void init() { }

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

    virtual void finalize() { }



  public:



    const AnalysisInfo& info() const {
      assert(_info && "No AnalysisInfo object :O");
      return *_info;
    }

    virtual std::string name() const {
      return  ( (info().name().empty()) ? _defaultname : info().name() ) + _optstring;
    }
    // get name of reference data file, which could be different from plugin name
    virtual std::string getRefDataName() const {
      return (info().getRefDataName().empty()) ? _defaultname : info().getRefDataName();
    }

    // set name of reference data file, which could be different from plugin name
    virtual void setRefDataName(const std::string& ref_data="") {
      info().setRefDataName(!ref_data.empty() ? ref_data : name());
    }

    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 std::string luminosityfb() const {
      return info().luminosityfb();
    }

    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::vector<std::string> todos() const {
      return info().todos();
    }


    virtual const std::vector<PdgIdPair>& requiredBeams() const {
      return info().beams();
    }
    virtual Analysis& setRequiredBeams(const std::vector<PdgIdPair>& requiredBeams) {
      info().setBeams(requiredBeams);
      return *this;
    }


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

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

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


    bool needsCrossSection() const {
      return info().needsCrossSection();
    }
    Analysis& setNeedsCrossSection(bool needed=true) {
      info().setNeedsCrossSection(needed);
      return *this;
    }





    AnalysisInfo& info() {
      assert(_info && "No AnalysisInfo object :O");
      return *_info;
    }





    const ParticlePair& beams() const;

    const PdgIdPair beamIds() const;

    double sqrtS() const;

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





    bool isCompatible(const ParticlePair& beams) const;

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

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



    Analysis& setCrossSection(double xs); //, double xserr=0.0);

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



  protected:

    Log& getLog() const;

    double crossSection() const;

    double crossSectionPerEvent() const;

    size_t numEvents() 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;
    const std::string makeAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
      return mkAxisCode(datasetId, xAxisId, yAxisId);
    }





    template <typename T=YODA::Scatter2D>
    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.");
      }
      return dynamic_cast<T&>(*_refdata[hname]);
    }

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





    CounterPtr bookCounter(const std::string& name,
                           const std::string& title="");
                           // const std::string& valtitle=""

    CounterPtr bookCounter(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                           const std::string& title="");
                           // const std::string& valtitle=""





    Histo1DPtr bookHisto1D(const std::string& name,
                           size_t nbins, double lower, double upper,
                           const std::string& title="",
                           const std::string& xtitle="",
                           const std::string& ytitle="");

    Histo1DPtr bookHisto1D(const std::string& name,
                           const std::vector<double>& binedges,
                           const std::string& title="",
                           const std::string& xtitle="",
                           const std::string& ytitle="");

    Histo1DPtr bookHisto1D(const std::string& name,
                           const std::initializer_list<double>& binedges,
                           const std::string& title="",
                           const std::string& xtitle="",
                           const std::string& ytitle="");

    Histo1DPtr bookHisto1D(const std::string& name,
                           const Scatter2D& refscatter,
                           const std::string& title="",
                           const std::string& xtitle="",
                           const std::string& ytitle="");

    Histo1DPtr bookHisto1D(const std::string& name,
                           const std::string& title="",
                           const std::string& xtitle="",
                           const std::string& ytitle="");

    Histo1DPtr bookHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                           const std::string& title="",
                           const std::string& xtitle="",
                           const std::string& ytitle="");





    Histo2DPtr bookHisto2D(const std::string& name,
                           size_t nxbins, double xlower, double xupper,
                           size_t nybins, double ylower, double yupper,
                           const std::string& title="",
                           const std::string& xtitle="",
                           const std::string& ytitle="",
                           const std::string& ztitle="");

    Histo2DPtr bookHisto2D(const std::string& name,
                           const std::vector<double>& xbinedges,
                           const std::vector<double>& ybinedges,
                           const std::string& title="",
                           const std::string& xtitle="",
                           const std::string& ytitle="",
                           const std::string& ztitle="");

    Histo2DPtr bookHisto2D(const std::string& name,
                           const std::initializer_list<double>& xbinedges,
                           const std::initializer_list<double>& ybinedges,
                           const std::string& title="",
                           const std::string& xtitle="",
                           const std::string& ytitle="",
                           const std::string& ztitle="");

    Histo2DPtr bookHisto2D(const std::string& name,
                           const Scatter3D& refscatter,
                           const std::string& title="",
                           const std::string& xtitle="",
                           const std::string& ytitle="",
                           const std::string& ztitle="");

    Histo2DPtr bookHisto2D(const std::string& name,
                           const std::string& title="",
                           const std::string& xtitle="",
                           const std::string& ytitle="",
                           const std::string& ztitle="");

    Histo2DPtr bookHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                           const std::string& title="",
                           const std::string& xtitle="",
                           const std::string& ytitle="",
                           const std::string& ztitle="");





    Profile1DPtr bookProfile1D(const std::string& name,
                               size_t nbins, double lower, double upper,
                               const std::string& title="",
                               const std::string& xtitle="",
                               const std::string& ytitle="");

    Profile1DPtr bookProfile1D(const std::string& name,
                               const std::vector<double>& binedges,
                               const std::string& title="",
                               const std::string& xtitle="",
                               const std::string& ytitle="");

    Profile1DPtr bookProfile1D(const std::string& name,
                               const std::initializer_list<double>& binedges,
                               const std::string& title="",
                               const std::string& xtitle="",
                               const std::string& ytitle="");

    Profile1DPtr bookProfile1D(const std::string& name,
                               const Scatter2D& refscatter,
                               const std::string& title="",
                               const std::string& xtitle="",
                               const std::string& ytitle="");

    Profile1DPtr bookProfile1D(const std::string& name,
                               const std::string& title="",
                               const std::string& xtitle="",
                               const std::string& ytitle="");

    Profile1DPtr bookProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                               const std::string& title="",
                               const std::string& xtitle="",
                               const std::string& ytitle="");





    Profile2DPtr bookProfile2D(const std::string& name,
                               size_t nxbins, double xlower, double xupper,
                               size_t nybins, double ylower, double yupper,
                               const std::string& title="",
                               const std::string& xtitle="",
                               const std::string& ytitle="",
                               const std::string& ztitle="");

    Profile2DPtr bookProfile2D(const std::string& name,
                               const std::vector<double>& xbinedges,
                               const std::vector<double>& ybinedges,
                               const std::string& title="",
                               const std::string& xtitle="",
                               const std::string& ytitle="",
                               const std::string& ztitle="");

    Profile2DPtr bookProfile2D(const std::string& name,
                               const std::initializer_list<double>& xbinedges,
                               const std::initializer_list<double>& ybinedges,
                               const std::string& title="",
                               const std::string& xtitle="",
                               const std::string& ytitle="",
                               const std::string& ztitle="");

    Profile2DPtr bookProfile2D(const std::string& name,
                               const Scatter3D& refscatter,
                               const std::string& title="",
                               const std::string& xtitle="",
                               const std::string& ytitle="",
                               const std::string& ztitle="");

    Profile2DPtr bookProfile2D(const std::string& name,
                               const std::string& title="",
                               const std::string& xtitle="",
                               const std::string& ytitle="",
                               const std::string& ztitle="");

    Profile2DPtr bookProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                               const std::string& title="",
                               const std::string& xtitle="",
                               const std::string& ytitle="",
                               const std::string& ztitle="");





    Scatter2DPtr bookScatter2D(const std::string& name,
                               bool copy_pts=false,
                               const std::string& title="",
                               const std::string& xtitle="",
                               const std::string& ytitle="");

    Scatter2DPtr bookScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                               bool copy_pts=false,
                               const std::string& title="",
                               const std::string& xtitle="",
                               const std::string& ytitle="");

    Scatter2DPtr bookScatter2D(const std::string& name,
                               size_t npts, double lower, double upper,
                               const std::string& title="",
                               const std::string& xtitle="",
                               const std::string& ytitle="");

    Scatter2DPtr bookScatter2D(const std::string& hname,
                               const std::vector<double>& binedges,
                               const std::string& title,
                               const std::string& xtitle,
                               const std::string& ytitle);

    Scatter2DPtr bookScatter2D(const Scatter2DPtr scPtr, 
    const std::string& path, const std::string& title = "", 
    const std::string& xtitle = "", const std::string& ytitle = "" );
    


  public:



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

    template<typename T>
    T getOption(std::string optname, T def) {
      if (_options.find(optname) == _options.end()) return def;
      std::stringstream ss;
      ss << _options.find(optname)->second;
      T ret;
      ss >> ret;
      return ret;
    }

    const CentralityProjection&
    declareCentrality(const SingleValueProjection &proj,
                      string calAnaName, string calHistName,
                      const string projName, bool increasing = false);

    template <class T>
    Percentile<T> bookPercentile(string projName,
                                 vector<pair<float, float> > centralityBins,
                                 vector<tuple<int, int, int> > ref) {
      typedef typename ReferenceTraits<T>::RefT RefT;
      Percentile<T> pctl(this, projName);

      const int nCent = centralityBins.size();
      for (int iCent = 0; iCent < nCent; ++iCent) {
        const string axisCode = makeAxisCode(std::get<0>(ref[iCent]),
                                             std::get<1>(ref[iCent]),
                                             std::get<2>(ref[iCent])); 
    const RefT & refscatter = refData<RefT>(axisCode);
        shared_ptr<T> ao = addOrGetCompatAO(make_shared<T>(refscatter, histoPath(axisCode)));
        CounterPtr cnt =
          addOrGetCompatAO(make_shared<Counter>(histoPath("TMP/COUNTER/" + axisCode)));
        pctl.add(ao, cnt, centralityBins[iCent]);
      }
      return pctl;
    }

    template <class T>
    PercentileXaxis<T> bookPercentileXaxis(string projName,
                                           tuple<int, int, int> ref) {
      typedef typename ReferenceTraits<T>::RefT RefT;
      PercentileXaxis<T> pctl(this, projName);

      const string axisCode = makeAxisCode(std::get<0>(ref),
                                           std::get<1>(ref),
                                           std::get<2>(ref)); 
      const RefT & refscatter = refData<RefT>(axisCode);
      shared_ptr<T> ao = addOrGetCompatAO(make_shared<T>(refscatter, histoPath(axisCode)));
      pctl.add(proj, ao, make_shared<Counter>());

      return pctl;
    }





    void scale(CounterPtr cnt, double factor);

    void scale(const std::vector<CounterPtr>& cnts, double factor) {
      for (auto& c : cnts) scale(c, factor);
    }
    template <std::size_t array_size>
    void scale(const CounterPtr (&cnts)[array_size], double factor) {
      // for (size_t i = 0; i < std::extent<decltype(cnts)>::value; ++i) scale(cnts[i], factor);
      for (auto& c : cnts) scale(c, factor);
    }


    void normalize(Histo1DPtr histo, double norm=1.0, bool includeoverflows=true);

    void normalize(const std::vector<Histo1DPtr>& histos, double norm=1.0, bool includeoverflows=true) {
      for (auto& h : histos) normalize(h, norm, includeoverflows);
    }
    template <std::size_t array_size>
    void normalize(const Histo1DPtr (&histos)[array_size], double norm=1.0, bool includeoverflows=true) {
      for (auto& h : histos) normalize(h, norm, includeoverflows);
    }

    void scale(Histo1DPtr histo, double factor);

    void scale(const std::vector<Histo1DPtr>& histos, double factor) {
      for (auto& h : histos) scale(h, factor);
    }
    template <std::size_t array_size>
    void scale(const Histo1DPtr (&histos)[array_size], double factor) {
      for (auto& h : histos) scale(h, factor);
    }


    void normalize(Histo2DPtr histo, double norm=1.0, bool includeoverflows=true);

    void normalize(const std::vector<Histo2DPtr>& histos, double norm=1.0, bool includeoverflows=true) {
      for (auto& h : histos) normalize(h, norm, includeoverflows);
    }
    template <std::size_t array_size>
    void normalize(const Histo2DPtr (&histos)[array_size], double norm=1.0, bool includeoverflows=true) {
      for (auto& h : histos) normalize(h, norm, includeoverflows);
    }

    void scale(Histo2DPtr histo, double factor);

    void scale(const std::vector<Histo2DPtr>& histos, double factor) {
      for (auto& h : histos) scale(h, factor);
    }
    template <std::size_t array_size>
    void scale(const Histo2DPtr (&histos)[array_size], double factor) {
      for (auto& h : histos) scale(h, factor);
    }


    void divide(CounterPtr c1, CounterPtr c2, Scatter1DPtr s) const;

    void divide(const YODA::Counter& c1, const YODA::Counter& c2, Scatter1DPtr s) const;


    void divide(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const;

    void divide(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const;


    void divide(Profile1DPtr p1, Profile1DPtr p2, Scatter2DPtr s) const;

    void divide(const YODA::Profile1D& p1, const YODA::Profile1D& p2, Scatter2DPtr s) const;


    void divide(Histo2DPtr h1, Histo2DPtr h2, Scatter3DPtr s) const;

    void divide(const YODA::Histo2D& h1, const YODA::Histo2D& h2, Scatter3DPtr s) const;


    void divide(Profile2DPtr p1, Profile2DPtr p2, Scatter3DPtr s) const;

    void divide(const YODA::Profile2D& p1, const YODA::Profile2D& p2, Scatter3DPtr s) const;


    void efficiency(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const;

    void efficiency(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const;


    void asymm(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const;

    void asymm(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const;


    void integrate(Histo1DPtr h, Scatter2DPtr s) const;

    void integrate(const Histo1D& h, Scatter2DPtr s) const;



  public:

    const vector<AnalysisObjectPtr>& analysisObjects() const {
      return _analysisobjects;
    }


  protected:



    void addAnalysisObject(AnalysisObjectPtr ao);

    template <typename AO=YODA::AnalysisObject>
    std::shared_ptr<AO> addOrGetCompatAO(std::shared_ptr<AO> aonew) {
      foreach (const AnalysisObjectPtr& ao, analysisObjects()) {
        if ( ao->path() == aonew->path() ) {
           std::shared_ptr<AO> aoold = dynamic_pointer_cast<AO>(ao);
           if ( aoold && bookingCompatible(aonew, aoold) ) {
             MSG_TRACE("Bound pre-existing data object " << aonew->path()
                       <<  " for " << name());
             return aoold;
           } else {
             MSG_WARNING("Found incompatible pre-existing data object with same path " 
                         << aonew->path() <<  " for " << name());
           }
        }
      }
      MSG_TRACE("Registered " << aonew->annotation("Type") << " " << aonew->path()
                <<  " for " << name());
      addAnalysisObject(aonew);
      return aonew;
    }

    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");
    }

    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");
    }

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

    void removeAnalysisObject(AnalysisObjectPtr ao);

    vector<AnalysisObjectPtr> getAllData(bool includeorphans) const;

    template <typename AO=YODA::AnalysisObject>
    std::shared_ptr<AO> getAnalysisObject(const std::string & ananame,
                                          const std::string& name) {
      std::string path = "/" + ananame + "/" + name;
      for ( AnalysisObjectPtr ao : getAllData(true) ) {
        if ( ao->path() == path )
          return dynamic_pointer_cast<AO>(ao);
      }
      return std::shared_ptr<AO>();
    }

    const Histo1DPtr getHisto1D(const std::string& name) const {
      return getAnalysisObject<Histo1D>(name);
    }

    Histo1DPtr getHisto1D(const std::string& name) {
      return getAnalysisObject<Histo1D>(name);
    }

    const Histo1DPtr getHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
      return getAnalysisObject<Histo1D>(makeAxisCode(datasetId, xAxisId, yAxisId));
    }

    Histo1DPtr getHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) {
      return getAnalysisObject<Histo1D>(makeAxisCode(datasetId, xAxisId, yAxisId));
    }


    const Histo2DPtr getHisto2D(const std::string& name) const {
      return getAnalysisObject<Histo2D>(name);
    }

    Histo2DPtr getHisto2D(const std::string& name) {
      return getAnalysisObject<Histo2D>(name);
    }

    const Histo2DPtr getHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
      return getAnalysisObject<Histo2D>(makeAxisCode(datasetId, xAxisId, yAxisId));
    }

    Histo2DPtr getHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) {
      return getAnalysisObject<Histo2D>(makeAxisCode(datasetId, xAxisId, yAxisId));
    }


    const Profile1DPtr getProfile1D(const std::string& name) const {
      return getAnalysisObject<Profile1D>(name);
    }

    Profile1DPtr getProfile1D(const std::string& name) {
      return getAnalysisObject<Profile1D>(name);
    }

    const Profile1DPtr getProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
      return getAnalysisObject<Profile1D>(makeAxisCode(datasetId, xAxisId, yAxisId));
    }

    Profile1DPtr getProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) {
      return getAnalysisObject<Profile1D>(makeAxisCode(datasetId, xAxisId, yAxisId));
    }


    const Profile2DPtr getProfile2D(const std::string& name) const {
      return getAnalysisObject<Profile2D>(name);
    }

    Profile2DPtr getProfile2D(const std::string& name) {
      return getAnalysisObject<Profile2D>(name);
    }

    const Profile2DPtr getProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
      return getAnalysisObject<Profile2D>(makeAxisCode(datasetId, xAxisId, yAxisId));
    }

    Profile2DPtr getProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) {
      return getAnalysisObject<Profile2D>(makeAxisCode(datasetId, xAxisId, yAxisId));
    }


    const Scatter2DPtr getScatter2D(const std::string& name) const {
      return getAnalysisObject<Scatter2D>(name);
    }

    Scatter2DPtr getScatter2D(const std::string& name) {
      return getAnalysisObject<Scatter2D>(name);
    }

    const Scatter2DPtr getScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
      return getAnalysisObject<Scatter2D>(makeAxisCode(datasetId, xAxisId, yAxisId));
    }

    Scatter2DPtr getScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) {
      return getAnalysisObject<Scatter2D>(makeAxisCode(datasetId, xAxisId, yAxisId));
    }



  private:

    string _defaultname;

    unique_ptr<AnalysisInfo> _info;

    vector<AnalysisObjectPtr> _analysisobjects;
    

    double _crossSection;
    bool _gotCrossSection;

    AnalysisHandler* _analysishandler;

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

     map<string, string> _options;

    string _optstring;

  private:



    void _cacheRefData() const;



    Analysis& operator=(const Analysis&);

  };


}


// Include definition of analysis plugin system so that analyses automatically see it when including Analysis.hh
#include "Rivet/AnalysisBuilder.hh"

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

// #define DECLARE_ALIASED_RIVET_PLUGIN(clsname, alias) Rivet::AnalysisBuilder<clsname> plugin_ ## clsname ## ( ## #alias ## )
#define DECLARE_ALIASED_RIVET_PLUGIN(clsname, alias) DECLARE_RIVET_PLUGIN(clsname)( #alias )

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

#define DEFAULT_RIVET_ANALYSIS_CTOR(clsname) DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR(clsname)



#endif