diff --git a/include/Rivet/Analysis.hh b/include/Rivet/Analysis.hh
--- a/include/Rivet/Analysis.hh
+++ b/include/Rivet/Analysis.hh
@@ -1,1333 +1,1335 @@
 // -*- 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"
 #include <tuple>
 
 
 /// @def vetoEvent
 /// Preprocessor define for vetoing events, including the log message and return.
 #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;
 
 
   /// @brief This is the base class of all analysis classes in Rivet.
   ///
   /// There are
   /// three virtual functions which should be implemented in base classes:
   ///
   /// void init() is called by Rivet before a run is started. Here the
   /// analysis class should book necessary histograms. The needed
   /// projections should probably rather be constructed in the
   /// constructor.
   ///
   /// void analyze(const Event&) is called once for each event. Here the
   /// analysis class should apply the necessary Projections and fill the
   /// histograms.
   ///
   /// void finalize() is called after a run is finished. Here the analysis
   /// class should do whatever manipulations are necessary on the
   /// histograms. Writing the histograms to a file is, however, done by
   /// the Rivet class.
   class Analysis : public ProjectionApplier {
 
     /// The AnalysisHandler is a friend.
     friend class AnalysisHandler;
 
 
   public:
 
     /// @name Standard constructors and destructors.
     //@{
 
     // /// The default constructor.
     // Analysis();
 
     /// Constructor
     Analysis(const std::string& name);
 
     /// The destructor.
     virtual ~Analysis() {}
 
     //@}
 
 
   public:
 
     /// @name Main analysis methods
     //@{
 
     /// Initialize this analysis object. A concrete class should here
     /// book all necessary histograms. An overridden function must make
     /// sure it first calls the base class function.
     virtual void init() { }
 
     /// Analyze one event. A concrete class should here apply the
     /// necessary projections on the \a event and fill the relevant
     /// histograms. An overridden function must make sure it first calls
     /// the base class function.
     virtual void analyze(const Event& event) = 0;
 
     /// Finalize this analysis object. A concrete class should here make
     /// all necessary operations on the histograms. Writing the
     /// histograms to a file is, however, done by the Rivet class. An
     /// overridden function must make sure it first calls the base class
     /// function.
     virtual void finalize() { }
 
     //@}
 
 
   public:
 
     /// @name Metadata
     /// Metadata is used for querying from the command line and also for
     /// building web pages and the analysis pages in the Rivet manual.
     //@{
 
     /// Get the actual AnalysisInfo object in which all this metadata is stored.
     const AnalysisInfo& info() const {
       assert(_info && "No AnalysisInfo object :O");
       return *_info;
     }
 
     /// @brief Get the name of the analysis.
     ///
     /// By default this is computed by combining the results of the
     /// experiment, year and Spires ID metadata methods and you should
     /// only override it if there's a good reason why those won't
     /// work. If options has been set for this instance, a
     /// corresponding string is appended at the end.
     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());
     }
 
     /// Get the Inspire ID code for this analysis.
     virtual std::string inspireId() const {
       return info().inspireId();
     }
 
     /// Get the SPIRES ID code for this analysis (~deprecated).
     virtual std::string spiresId() const {
       return info().spiresId();
     }
 
     /// @brief Names & emails of paper/analysis authors.
     ///
     /// Names and email of authors in 'NAME \<EMAIL\>' format. The first
     /// name in the list should be the primary contact person.
     virtual std::vector<std::string> authors() const {
       return info().authors();
     }
 
     /// @brief Get a short description of the analysis.
     ///
     /// Short (one sentence) description used as an index entry.
     /// Use @a description() to provide full descriptive paragraphs
     /// of analysis details.
     virtual std::string summary() const {
       return info().summary();
     }
 
     /// @brief Get a full description of the analysis.
     ///
     /// Full textual description of this analysis, what it is useful for,
     /// what experimental techniques are applied, etc. Should be treated
     /// as a chunk of restructuredText (http://docutils.sourceforge.net/rst.html),
     /// with equations to be rendered as LaTeX with amsmath operators.
     virtual std::string description() const {
       return info().description();
     }
 
     /// @brief Information about the events needed as input for this analysis.
     ///
     /// Event types, energies, kinematic cuts, particles to be considered
     /// stable, etc. etc. Should be treated as a restructuredText bullet list
     /// (http://docutils.sourceforge.net/rst.html)
     virtual std::string runInfo() const {
       return info().runInfo();
     }
 
     /// Experiment which performed and published this analysis.
     virtual std::string experiment() const {
       return info().experiment();
     }
 
     /// Collider on which the experiment ran.
     virtual std::string collider() const {
       return info().collider();
     }
 
     /// When the original experimental analysis was published.
     virtual std::string year() const {
       return info().year();
     }
 
     /// The luminosity in inverse femtobarn
     virtual std::string luminosityfb() const {
       return info().luminosityfb();
     }
 
     /// Journal, and preprint references.
     virtual std::vector<std::string> references() const {
       return info().references();
     }
 
     /// BibTeX citation key for this article.
     virtual std::string bibKey() const {
       return info().bibKey();
     }
 
     /// BibTeX citation entry for this article.
     virtual std::string bibTeX() const {
       return info().bibTeX();
     }
 
     /// Whether this analysis is trusted (in any way!)
     virtual std::string status() const {
       return (info().status().empty()) ? "UNVALIDATED" : info().status();
     }
 
     /// Any work to be done on this analysis.
     virtual std::vector<std::string> todos() const {
       return info().todos();
     }
 
 
     /// Return the allowed pairs of incoming beams required by this analysis.
     virtual const std::vector<PdgIdPair>& requiredBeams() const {
       return info().beams();
     }
     /// Declare the allowed pairs of incoming beams required by this analysis.
     virtual Analysis& setRequiredBeams(const std::vector<PdgIdPair>& requiredBeams) {
       info().setBeams(requiredBeams);
       return *this;
     }
 
 
     /// Sets of valid beam energy pairs, in GeV
     virtual const std::vector<std::pair<double, double> >& requiredEnergies() const {
       return info().energies();
     }
 
     /// Get vector of analysis keywords
     virtual const std::vector<std::string> & keywords() const {
       return info().keywords();
     }
 
     /// Declare the list of valid beam energy pairs, in GeV
     virtual Analysis& setRequiredEnergies(const std::vector<std::pair<double, double> >& requiredEnergies) {
       info().setEnergies(requiredEnergies);
       return *this;
     }
 
     //@}
 
 
     /// @name Internal metadata modifying methods
     //@{
 
     /// Get the actual AnalysisInfo object in which all this metadata is stored (non-const).
     AnalysisInfo& info() {
       assert(_info && "No AnalysisInfo object :O");
       return *_info;
     }
 
     //@}
 
 
     /// @name Run conditions
     //@{
 
     /// Incoming beams for this run
     const ParticlePair& beams() const;
 
     /// Incoming beam IDs for this run
     const PdgIdPair beamIds() const;
 
     /// Centre of mass energy for this run
     double sqrtS() const;
 
     /// Check if we are running rivet-merge.
     bool merging() const {
       return sqrtS() <= 0.0;
     }
 
     //@}
 
 
     /// @name Analysis / beam compatibility testing
     /// @todo Replace with beamsCompatible() with no args (calling beams() function internally)
     /// @todo Add beamsMatch() methods with same (shared-code?) tolerance as in beamsCompatible()
     //@{
 
     /// Check if analysis is compatible with the provided beam particle IDs and energies
     bool isCompatible(const ParticlePair& beams) const;
 
     /// Check if analysis is compatible with the provided beam particle IDs and energies
     bool isCompatible(PdgId beam1, PdgId beam2, double e1, double e2) const;
 
     /// Check if analysis is compatible with the provided beam particle IDs and energies
     bool isCompatible(const PdgIdPair& beams, const std::pair<double,double>& energies) const;
 
     //@}
 
     /// Access the controlling AnalysisHandler object.
     AnalysisHandler& handler() const { return *_analysishandler; }
 
 
   protected:
 
     /// Get a Log object based on the name() property of the calling analysis object.
     Log& getLog() const;
 
     /// Get the process cross-section in pb. Throws if this hasn't been set.
     double crossSection() const;
 
     /// Get the process cross-section per generated event in pb. Throws if this
     /// hasn't been set.
     double crossSectionPerEvent() const;
 
     /// @brief Get the number of events seen (via the analysis handler).
     ///
     /// @note Use in the finalize phase only.
     size_t numEvents() const;
 
     /// @brief Get the sum of event weights seen (via the analysis handler).
     ///
     /// @note Use in the finalize phase only.
     double sumW() const;
     /// Alias
     double sumOfWeights() const { return sumW(); }
 
     /// @brief Get the sum of squared event weights seen (via the analysis handler).
     ///
     /// @note Use in the finalize phase only.
     double sumW2() const;
 
 
   protected:
 
     /// @name Histogram paths
     //@{
 
     /// Get the canonical histogram "directory" path for this analysis.
     const std::string histoDir() const;
 
     /// Get the canonical histogram path for the named histogram in this analysis.
     const std::string histoPath(const std::string& hname) const;
 
     /// Get the canonical histogram path for the numbered histogram in this analysis.
     const std::string histoPath(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const;
 
     /// Get the internal histogram name for given d, x and y (cf. HepData)
     const std::string mkAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const;
 
     //@}
 
 
     /// @name Histogram reference data
     //@{
 
     /// Get reference data for a named histo
     /// @todo SFINAE to ensure that the type inherits from YODA::AnalysisObject?
     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]);
     }
 
     /// Get reference data for a numbered histo
     /// @todo SFINAE to ensure that the type inherits from YODA::AnalysisObject?
     template <typename T=YODA::Scatter2D>
     const T& refData(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
       const string hname = mkAxisCode(datasetId, xAxisId, yAxisId);
       return refData(hname);
     }
 
     //@}
 
 
     /// @name Counter booking
     //@{
 
     /// Book a counter.
     CounterPtr & book(CounterPtr &, const std::string& name,
                            const std::string& title="");
                            // const std::string& valtitle=""
 
     /// Book a counter, using a path generated from the dataset and axis ID codes
     ///
     /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
     CounterPtr & book(CounterPtr &, unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                            const std::string& title="");
                            // const std::string& valtitle=""
 
     //@}
 
 
     /// @name 1D histogram booking
     //@{
 
     /// Book a 1D histogram with @a nbins uniformly distributed across the range @a lower - @a upper .
     Histo1DPtr & book(Histo1DPtr &,const std::string& name,
                            size_t nbins, double lower, double upper,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="");
 
     /// Book a 1D histogram with non-uniform bins defined by the vector of bin edges @a binedges .
     Histo1DPtr & book(Histo1DPtr &,const std::string& name,
                            const std::vector<double>& binedges,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="");
 
     /// Book a 1D histogram with non-uniform bins defined by the vector of bin edges @a binedges .
     Histo1DPtr & book(Histo1DPtr &,const std::string& name,
                            const std::initializer_list<double>& binedges,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="");
 
     /// Book a 1D histogram with binning from a reference scatter.
     Histo1DPtr & book(Histo1DPtr &,const std::string& name,
                            const Scatter2D& refscatter,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="");
 
     /// Book a 1D histogram, using the binnings in the reference data histogram.
     Histo1DPtr & book(Histo1DPtr &,const std::string& name,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="");
 
     /// Book a 1D histogram, using the binnings in the reference data histogram.
     ///
     /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
     Histo1DPtr & book(Histo1DPtr &,unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="");
 
     //@}
 
 
     /// @name 2D histogram booking
     //@{
 
     /// Book a 2D histogram with @a nxbins and @a nybins uniformly
     /// distributed across the ranges @a xlower - @a xupper and @a
     /// ylower - @a yupper respectively along the x- and y-axis.
     Histo2DPtr & book(Histo2DPtr &,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="");
 
     /// Book a 2D histogram with non-uniform bins defined by the
     /// vectors of bin edges @a xbinedges and @a ybinedges.
     Histo2DPtr & book(Histo2DPtr &,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="");
 
     /// Book a 2D histogram with non-uniform bins defined by the
     /// vectors of bin edges @a xbinedges and @a ybinedges.
     Histo2DPtr & book(Histo2DPtr &,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="");
 
     /// Book a 2D histogram with binning from a reference scatter.
     Histo2DPtr & book(Histo2DPtr &,const std::string& name,
                            const Scatter3D& refscatter,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="",
                            const std::string& ztitle="");
 
     /// Book a 2D histogram, using the binnings in the reference data histogram.
     Histo2DPtr & book(Histo2DPtr &,const std::string& name,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="",
                            const std::string& ztitle="");
 
     /// Book a 2D histogram, using the binnings in the reference data histogram.
     ///
     /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
     Histo2DPtr & book(Histo2DPtr &,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="");
 
     //@}
 
 
     /// @name 1D profile histogram booking
     //@{
 
     /// Book a 1D profile histogram with @a nbins uniformly distributed across the range @a lower - @a upper .
     Profile1DPtr & book(Profile1DPtr &,  const std::string& name,
                                size_t nbins, double lower, double upper,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="");
 
     /// Book a 1D profile histogram with non-uniform bins defined by the vector of bin edges @a binedges .
     Profile1DPtr & book(Profile1DPtr &,  const std::string& name,
                                const std::vector<double>& binedges,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="");
 
     /// Book a 1D profile histogram with non-uniform bins defined by the vector of bin edges @a binedges .
     Profile1DPtr & book(Profile1DPtr &,  const std::string& name,
                                const std::initializer_list<double>& binedges,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="");
 
     /// Book a 1D profile histogram with binning from a reference scatter.
     Profile1DPtr & book(Profile1DPtr &,  const std::string& name,
                                const Scatter2D& refscatter,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="");
 
     /// Book a 1D profile histogram, using the binnings in the reference data histogram.
     Profile1DPtr & book(Profile1DPtr &,  const std::string& name,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="");
 
     /// Book a 1D profile histogram, using the binnings in the reference data histogram.
     ///
     /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
     Profile1DPtr & book(Profile1DPtr &,  unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="");
 
     //@}
 
 
     /// @name 2D profile histogram booking
     //@{
 
     /// Book a 2D profile histogram with @a nxbins and @a nybins uniformly
     /// distributed across the ranges @a xlower - @a xupper and @a ylower - @a
     /// yupper respectively along the x- and y-axis.
     Profile2DPtr & book(Profile2DPtr &,  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="");
 
     /// Book a 2D profile histogram with non-uniform bins defined by the vectorx
     /// of bin edges @a xbinedges and @a ybinedges.
     Profile2DPtr & book(Profile2DPtr &,  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="");
 
     /// Book a 2D profile histogram with non-uniform bins defined by the vectorx
     /// of bin edges @a xbinedges and @a ybinedges.
     Profile2DPtr & book(Profile2DPtr &,  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="");
 
     /// @todo REINSTATE
 
     // /// Book a 2D profile histogram with binning from a reference scatter.
     // Profile2DPtr & book(const Profile2DPtr &, const std::string& name,
     //                            const Scatter3D& refscatter,
     //                            const std::string& title="",
     //                            const std::string& xtitle="",
     //                            const std::string& ytitle="",
     //                            const std::string& ztitle="");
 
     // /// Book a 2D profile histogram, using the binnings in the reference data histogram.
     // Profile2DPtr & book(const Profile2DPtr &, const std::string& name,
     //                            const std::string& title="",
     //                            const std::string& xtitle="",
     //                            const std::string& ytitle="",
     //                            const std::string& ztitle="");
 
     // /// Book a 2D profile histogram, using the binnings in the reference data histogram.
     // ///
     // /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
     // Profile2DPtr & book(const Profile2DPtr &, 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="");
 
     //@}
 
 
     /// @name 2D scatter booking
     //@{
 
     /// @brief Book a 2-dimensional data point set with the given name.
     ///
     /// @note Unlike histogram booking, scatter booking by default makes no
     /// attempt to use reference data to pre-fill the data object. If you want
     /// this, which is sometimes useful e.g. when the x-position is not really
     /// meaningful and can't be extracted from the data, then set the @a
     /// copy_pts parameter to true. This creates points to match the reference
     /// data's x values and errors, but with the y values and errors zeroed...
     /// assuming that there is a reference histo with the same name: if there
     /// isn't, an exception will be thrown.
     Scatter2DPtr & book(Scatter2DPtr & s2d, const string& hname,
                                bool copy_pts=false,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="");
 
     /// @brief Book a 2-dimensional data point set, using the binnings in the reference data histogram.
     ///
     /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
     ///
     /// @note Unlike histogram booking, scatter booking by default makes no
     /// attempt to use reference data to pre-fill the data object. If you want
     /// this, which is sometimes useful e.g. when the x-position is not really
     /// meaningful and can't be extracted from the data, then set the @a
     /// copy_pts parameter to true. This creates points to match the reference
     /// data's x values and errors, but with the y values and errors zeroed.
     Scatter2DPtr & book(Scatter2DPtr & s2d, 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="");
 
     /// @brief Book a 2-dimensional data point set with equally spaced x-points in a range.
     ///
     /// The y values and errors will be set to 0.
     Scatter2DPtr & book(Scatter2DPtr & s2d, const string& hname,
                                size_t npts, double lower, double upper,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="");
 
     /// @brief Book a 2-dimensional data point set based on provided contiguous "bin edges".
     ///
     /// The y values and errors will be set to 0.
     Scatter2DPtr & book(Scatter2DPtr & s2d, const string& hname,
                                const std::vector<double>& binedges,
                                const std::string& title,
                                const std::string& xtitle,
                                const std::string& ytitle);
 
     /// Book a 2-dimensional data point set with x-points from an existing scatter and a new path.
     Scatter2DPtr book(Scatter2DPtr & s2d, const Scatter2DPtr & scPtr,
                                const std::string& path,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="");
 
     //@}
 
 
   public:
 
     /// @name Accessing options for this Analysis instance.
     //@{
 
     /// Get an option for this analysis instance as a string.
     std::string getOption(std::string optname) {
       if ( _options.find(optname) != _options.end() )
         return _options.find(optname)->second;
       return "";
     }
 
     /// Get an option for this analysis instance converted to a
     /// specific type (given by the specified @a def value).
     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;
     }
 
     //@}
     /// @brief Book a CentralityProjection
     ///
     /// Using a SingleValueProjection, @a proj, giving the value of an
     /// experimental observable to be used as a centrality estimator,
     /// book a CentralityProjection based on the experimentally
     /// measured pecentiles of this observable (as given by the
     /// reference data for the @a calHistName histogram in the @a
     /// calAnaName analysis. If a preloaded file with the output of a
     /// run using the @a calAnaName analysis contains a valid
     /// generated @a calHistName histogram, it will be used as an
     /// optional percentile binning. Also if this preloaded file
     /// contains a histogram with the name @a calHistName with an
     /// appended "_IMP" This histogram will be used to add an optional
     /// centrality percentile based on the generated impact
     /// parameter. If @increasing is true, a low (high) value of @proj
     /// is assumed to correspond to a more peripheral (central) event.
     const CentralityProjection&
     declareCentrality(const SingleValueProjection &proj,
                       string calAnaName, string calHistName,
                       const string projName, bool increasing = false);
 
     /// @brief Book a Percentile wrapper around AnalysisObjects.
     ///
     /// Based on a previously registered CentralityProjection named @a
     /// projName book one AnalysisObject for each @a centralityBin and
     /// name them according to the corresponding code in the @a ref
     /// vector.
     ///
     /// @todo Convert to just be called book() cf. others
     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;
       typedef rivet_shared_ptr<Wrapper<T>> WrapT;
 
       Percentile<T> pctl(this, projName);
 
       const int nCent = centralityBins.size();
       for (int 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 wrappers 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 rivet_shared_ptr<Wrapper<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
 
     /// Get the list of weight names from the handler
     vector<string> _weightNames() const;
 
     /// Get the default/nominal weight index
     size_t _defaultWeightIndex() const;
 
     /// Get an AO from another analysis
     MultiweightAOPtr _getOtherAnalysisObject(const std::string & ananame, const std::string& name);
 
     /// Check that analysis objects aren't being booked/registered outside the init stage
     void _checkBookInit() const;
 
 
   private:
 
     /// To be used in finalize context only:
     class CounterAdapter {
     public:
       CounterAdapter(double          x) : x_(x                ) {}
 
       CounterAdapter(const YODA::Counter   & c) : x_(c.val()          ) {}
 
       // CounterAdapter(CounterPtr     cp) : x_(cp->val()        ) {}
 
       CounterAdapter(const YODA::Scatter1D & s) : x_(s.points()[0].x()) {
         assert( s.numPoints() == 1 || "Can only scale by a single value.");
       }
 
       // CounterAdapter(Scatter1DPtr   sp) : x_(sp->points()[0].x()) {
       // 	assert( sp->numPoints() == 1 || "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::Scatter1D & s) {
       assert( s.numPoints() == 1 );
       return s.points()[0].x();
     }
 
     /// @name Analysis object manipulation
     /// @todo Should really be protected: only public to keep BinnedHistogram happy for now...
     //@{
 
     /// Multiplicatively scale the given counter, @a cnt, by factor @s factor.
     void scale(CounterPtr cnt, CounterAdapter factor);
 
     /// Multiplicatively scale the given counters, @a cnts, by factor @s factor.
     /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh)
     /// @todo Use SFINAE for a generic iterable of CounterPtrs
     void scale(const std::vector<CounterPtr>& cnts, CounterAdapter factor) {
       for (auto& c : cnts) scale(c, factor);
     }
     /// @todo YUCK!
     template <std::size_t array_size>
     void scale(const CounterPtr (&cnts)[array_size], CounterAdapter factor) {
       // for (size_t i = 0; i < std::extent<decltype(cnts)>::value; ++i) scale(cnts[i], factor);
       for (auto& c : cnts) scale(c, factor);
     }
 
 
     /// Normalize the given histogram, @a histo, to area = @a norm.
     void normalize(Histo1DPtr histo, CounterAdapter norm=1.0, bool includeoverflows=true);
 
     /// Normalize the given histograms, @a histos, to area = @a norm.
     /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh)
     /// @todo Use SFINAE for a generic iterable of Histo1DPtrs
     void normalize(const std::vector<Histo1DPtr>& histos, CounterAdapter norm=1.0, bool includeoverflows=true) {
       for (auto& h : histos) normalize(h, norm, includeoverflows);
     }
     /// @todo YUCK!
     template <std::size_t array_size>
     void normalize(const Histo1DPtr (&histos)[array_size], CounterAdapter norm=1.0, bool includeoverflows=true) {
       for (auto& h : histos) normalize(h, norm, includeoverflows);
     }
 
     /// Multiplicatively scale the given histogram, @a histo, by factor @s factor.
     void scale(Histo1DPtr histo, CounterAdapter factor);
 
     /// Multiplicatively scale the given histograms, @a histos, by factor @s factor.
     /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh)
     /// @todo Use SFINAE for a generic iterable of Histo1DPtrs
     void scale(const std::vector<Histo1DPtr>& histos, CounterAdapter factor) {
       for (auto& h : histos) scale(h, factor);
     }
     /// @todo YUCK!
     template <std::size_t array_size>
     void scale(const Histo1DPtr (&histos)[array_size], CounterAdapter factor) {
       for (auto& h : histos) scale(h, factor);
     }
 
 
     /// Normalize the given histogram, @a histo, to area = @a norm.
     void normalize(Histo2DPtr histo, CounterAdapter norm=1.0, bool includeoverflows=true);
 
     /// Normalize the given histograms, @a histos, to area = @a norm.
     /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh)
     /// @todo Use SFINAE for a generic iterable of Histo2DPtrs
     void normalize(const std::vector<Histo2DPtr>& histos, CounterAdapter norm=1.0, bool includeoverflows=true) {
       for (auto& h : histos) normalize(h, norm, includeoverflows);
     }
     /// @todo YUCK!
     template <std::size_t array_size>
     void normalize(const Histo2DPtr (&histos)[array_size], CounterAdapter norm=1.0, bool includeoverflows=true) {
       for (auto& h : histos) normalize(h, norm, includeoverflows);
     }
 
     /// Multiplicatively scale the given histogram, @a histo, by factor @s factor.
     void scale(Histo2DPtr histo, CounterAdapter factor);
 
     /// Multiplicatively scale the given histograms, @a histos, by factor @s factor.
     /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh)
     /// @todo Use SFINAE for a generic iterable of Histo2DPtrs
     void scale(const std::vector<Histo2DPtr>& histos, CounterAdapter factor) {
       for (auto& h : histos) scale(h, factor);
     }
     /// @todo YUCK!
     template <std::size_t array_size>
     void scale(const Histo2DPtr (&histos)[array_size], CounterAdapter factor) {
       for (auto& h : histos) scale(h, factor);
     }
 
 
     /// Helper for counter division.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void divide(CounterPtr c1, CounterPtr c2, Scatter1DPtr s) const;
 
     /// Helper for histogram division with raw YODA objects.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void divide(const YODA::Counter& c1, const YODA::Counter& c2, Scatter1DPtr s) const;
 
 
     /// Helper for histogram division.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void divide(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const;
 
     /// Helper for histogram division with raw YODA objects.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void divide(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const;
 
 
     /// Helper for profile histogram division.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void divide(Profile1DPtr p1, Profile1DPtr p2, Scatter2DPtr s) const;
 
     /// Helper for profile histogram division with raw YODA objects.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void divide(const YODA::Profile1D& p1, const YODA::Profile1D& p2, Scatter2DPtr s) const;
 
 
     /// Helper for 2D histogram division.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void divide(Histo2DPtr h1, Histo2DPtr h2, Scatter3DPtr s) const;
 
     /// Helper for 2D histogram division with raw YODA objects.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void divide(const YODA::Histo2D& h1, const YODA::Histo2D& h2, Scatter3DPtr s) const;
 
 
     /// Helper for 2D profile histogram division.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void divide(Profile2DPtr p1, Profile2DPtr p2, Scatter3DPtr s) const;
 
     /// Helper for 2D profile histogram division with raw YODA objects
     ///
     /// @note Assigns to the (already registered) output scatter, @a s.  Preserves the path information of the target.
     void divide(const YODA::Profile2D& p1, const YODA::Profile2D& p2, Scatter3DPtr s) const;
 
 
     /// Helper for histogram efficiency calculation.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void efficiency(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const;
 
     /// Helper for histogram efficiency calculation.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void efficiency(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const;
 
 
     /// Helper for histogram asymmetry calculation.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void asymm(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const;
 
     /// Helper for histogram asymmetry calculation.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void asymm(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const;
 
 
     /// Helper for converting a differential histo to an integral one.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void integrate(Histo1DPtr h, Scatter2DPtr s) const;
 
     /// Helper for converting a differential histo to an integral one.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void integrate(const Histo1D& h, Scatter2DPtr s) const;
 
     //@}
 
 
   public:
 
     /// List of registered analysis data objects
     const vector<MultiweightAOPtr>& analysisObjects() const {
       return _analysisobjects;
     }
 
 
   protected:
 
     /// @name Data object registration, retrieval, and removal
     //@{
 
     /// Register a data object in the histogram system
     template <typename AO=MultiweightAOPtr>
     AO addAnalysisObject(const AO & aonew) {
       _checkBookInit();
 
       for (const MultiweightAOPtr& 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 rivet_shared_ptr directly
-        AO aoold = AO(dynamic_pointer_cast<typename AO::value_type::Inner>(ao.get())); //< OMG
+        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;
     }
 
     /// Unregister a data object from the histogram system (by name)
     void removeAnalysisObject(const std::string& path);
 
     /// Unregister a data object from the histogram system (by pointer)
     void removeAnalysisObject(const MultiweightAOPtr & ao);
 
     // /// Get all data objects, for all analyses, from the AnalysisHandler
     // /// @todo Can we remove this? Why not call handler().getData()?
     // vector<YODA::AnalysisObjectPtr> getAllData(bool includeorphans) const;
 
 
     /// Get a Rivet data object from the histogram system
     template <typename AO=MultiweightAOPtr>
-    const AO getAnalysisObject(const std::string& name) const {
+    const AO getAnalysisObject(const std::string& aoname) const {
       for (const MultiweightAOPtr& ao : analysisObjects()) {
         ao.get()->setActiveWeightIdx(_defaultWeightIndex());
-        if (ao->path() == histoPath(name)) {
-          return dynamic_pointer_cast<AO>(ao);
+        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(name) + " not found");
+      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");
     // }
 
 
     /// Get a data object from another analysis (e.g. preloaded
     /// calibration histogram).
     template <typename AO=MultiweightAOPtr>
-    AO getAnalysisObject(const std::string & ananame,
-                         const std::string& name) {
-      MultiweightAOPtr ao = _getOtherAnalysisObject(ananame, name);
-      return dynamic_pointer_cast<AO>(ao);
+    AO getAnalysisObject(const std::string& ananame,
+                         const std::string& aoname) {
+      MultiweightAOPtr ao = _getOtherAnalysisObject(ananame, aoname);
+      // return dynamic_pointer_cast<AO>(ao);
+      return AO(dynamic_pointer_cast<typename AO::value_type>(ao.get()));
     }
 
 
     // /// Get a named Histo1D object from the histogram system
     // const Histo1DPtr getHisto1D(const std::string& name) const {
     //   return getAnalysisObject<Histo1D>(name);
     // }
 
     // /// Get a named Histo1D object from the histogram system (non-const)
     // Histo1DPtr getHisto1D(const std::string& name) {
     //   return getAnalysisObject<Histo1D>(name);
     // }
 
     // /// Get a Histo1D object from the histogram system by axis ID codes (non-const)
     // const Histo1DPtr getHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
     //   return getAnalysisObject<Histo1D>(makeAxisCode(datasetId, xAxisId, yAxisId));
     // }
 
     // /// Get a Histo1D object from the histogram system by axis ID codes (non-const)
     // Histo1DPtr getHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) {
     //   return getAnalysisObject<Histo1D>(makeAxisCode(datasetId, xAxisId, yAxisId));
     // }
 
 
     // /// Get a named Histo2D object from the histogram system
     // const Histo2DPtr getHisto2D(const std::string& name) const {
     //   return getAnalysisObject<Histo2D>(name);
     // }
 
     // /// Get a named Histo2D object from the histogram system (non-const)
     // Histo2DPtr getHisto2D(const std::string& name) {
     //   return getAnalysisObject<Histo2D>(name);
     // }
 
     // /// Get a Histo2D object from the histogram system by axis ID codes (non-const)
     // const Histo2DPtr getHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
     //   return getAnalysisObject<Histo2D>(makeAxisCode(datasetId, xAxisId, yAxisId));
     // }
 
     // /// Get a Histo2D object from the histogram system by axis ID codes (non-const)
     // Histo2DPtr getHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) {
     //   return getAnalysisObject<Histo2D>(makeAxisCode(datasetId, xAxisId, yAxisId));
     // }
 
 
     // /// Get a named Profile1D object from the histogram system
     // const Profile1DPtr getProfile1D(const std::string& name) const {
     //   return getAnalysisObject<Profile1D>(name);
     // }
 
     // /// Get a named Profile1D object from the histogram system (non-const)
     // Profile1DPtr getProfile1D(const std::string& name) {
     //   return getAnalysisObject<Profile1D>(name);
     // }
 
     // /// Get a Profile1D object from the histogram system by axis ID codes (non-const)
     // const Profile1DPtr getProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
     //   return getAnalysisObject<Profile1D>(makeAxisCode(datasetId, xAxisId, yAxisId));
     // }
 
     // /// Get a Profile1D object from the histogram system by axis ID codes (non-const)
     // Profile1DPtr getProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) {
     //   return getAnalysisObject<Profile1D>(makeAxisCode(datasetId, xAxisId, yAxisId));
     // }
 
 
     // /// Get a named Profile2D object from the histogram system
     // const Profile2DPtr getProfile2D(const std::string& name) const {
     //   return getAnalysisObject<Profile2D>(name);
     // }
 
     // /// Get a named Profile2D object from the histogram system (non-const)
     // Profile2DPtr getProfile2D(const std::string& name) {
     //   return getAnalysisObject<Profile2D>(name);
     // }
 
     // /// Get a Profile2D object from the histogram system by axis ID codes (non-const)
     // const Profile2DPtr getProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
     //   return getAnalysisObject<Profile2D>(makeAxisCode(datasetId, xAxisId, yAxisId));
     // }
 
     // /// Get a Profile2D object from the histogram system by axis ID codes (non-const)
     // Profile2DPtr getProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) {
     //   return getAnalysisObject<Profile2D>(makeAxisCode(datasetId, xAxisId, yAxisId));
     // }
 
 
     // /// Get a named Scatter2D object from the histogram system
     // const Scatter2DPtr getScatter2D(const std::string& name) const {
     //   return getAnalysisObject<Scatter2D>(name);
     // }
 
     // /// Get a named Scatter2D object from the histogram system (non-const)
     // Scatter2DPtr getScatter2D(const std::string& name) {
     //   return getAnalysisObject<Scatter2D>(name);
     // }
 
     // /// Get a Scatter2D object from the histogram system by axis ID codes (non-const)
     // const Scatter2DPtr getScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
     //   return getAnalysisObject<Scatter2D>(makeAxisCode(datasetId, xAxisId, yAxisId));
     // }
 
     // /// Get a Scatter2D object from the histogram system by axis ID codes (non-const)
     // Scatter2DPtr getScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) {
     //   return getAnalysisObject<Scatter2D>(makeAxisCode(datasetId, xAxisId, yAxisId));
     // }
 
     //@}
 
 
   private:
 
     /// Name passed to constructor (used to find .info analysis data file, and as a fallback)
     string _defaultname;
 
     /// Pointer to analysis metadata object
     unique_ptr<AnalysisInfo> _info;
 
     /// Storage of all plot objects
     /// @todo Make this a map for fast lookup by path?
     vector<MultiweightAOPtr> _analysisobjects;
 
     /// @name Cross-section variables
     //@{
     double _crossSection;
     bool _gotCrossSection;
     //@}
 
     /// The controlling AnalysisHandler object.
     AnalysisHandler* _analysishandler;
 
     /// Collection of cached refdata to speed up many autobookings: the
     /// reference data file should only be read once.
     mutable std::map<std::string, YODA::AnalysisObjectPtr> _refdata;
 
      /// Options the (this instance of) the analysis
      map<string, string> _options;
 
     /// The string of options.
     string _optstring;
 
   private:
 
     /// @name Utility functions
     //@{
 
     /// Get the reference data for this paper and cache it.
     void _cacheRefData() const;
 
     //@}
 
 
     /// The assignment operator is private and must never be called.
     /// In fact, it should not even be implemented.
     Analysis& operator=(const Analysis&);
 
   };
 
 
 }
 
 
 // Include definition of analysis plugin system so that analyses automatically see it when including Analysis.hh
 #include "Rivet/AnalysisBuilder.hh"
 
 /// @def DECLARE_RIVET_PLUGIN
 /// Preprocessor define to prettify the global-object plugin hook mechanism.
 #define DECLARE_RIVET_PLUGIN(clsname) Rivet::AnalysisBuilder<clsname> plugin_ ## clsname
 
 /// @def DECLARE_ALIASED_RIVET_PLUGIN
 /// Preprocessor define to prettify the global-object plugin hook mechanism, with an extra alias name for this analysis.
 // #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 )
 
 /// @def DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR
 /// Preprocessor define to prettify the manky constructor with name string argument
 #define DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR(clsname) clsname() : Analysis(# clsname) {}
 
 /// @def DEFAULT_RIVET_ANALYSIS_CTOR
 /// Slight abbreviation for DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR
 #define DEFAULT_RIVET_ANALYSIS_CTOR(clsname) DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR(clsname)
 
 
 
 #endif
diff --git a/include/Rivet/AnalysisHandler.hh b/include/Rivet/AnalysisHandler.hh
--- a/include/Rivet/AnalysisHandler.hh
+++ b/include/Rivet/AnalysisHandler.hh
@@ -1,342 +1,358 @@
 // -*- C++ -*-
 #ifndef RIVET_RivetHandler_HH
 #define RIVET_RivetHandler_HH
 
 #include "Rivet/Config/RivetCommon.hh"
 #include "Rivet/Particle.hh"
 #include "Rivet/AnalysisLoader.hh"
 #include "Rivet/Tools/RivetYODA.hh"
 
 namespace Rivet {
 
 
   // Forward declaration and smart pointer for Analysis
   class Analysis;
   typedef std::shared_ptr<Analysis> AnaHandle;
 
 
   // Needed to make smart pointers compare equivalent in the STL set
   struct CmpAnaHandle {
     bool operator() (const AnaHandle& a, const AnaHandle& b) const {
       return a.get() < b.get();
     }
   };
 
 
   /// A class which handles a number of analysis objects to be applied to
   /// generated events. An {@link Analysis}' AnalysisHandler is also responsible
   /// for handling the final writing-out of histograms.
   class AnalysisHandler {
   public:
 
     /// @name Constructors and destructors. */
     //@{
 
     /// Preferred constructor, with optional run name.
     AnalysisHandler(const string& runname="");
 
     /// @brief Destructor
     /// The destructor is not virtual, as this class should not be inherited from.
     ~AnalysisHandler();
 
     //@}
 
 
   private:
 
     /// Get a logger object.
     Log& getLog() const;
 
 
   public:
 
     /// @name Run properties
     //@{
 
     /// Get the name of this run.
     string runName() const;
 
 
     /// Get the number of events seen. Should only really be used by external
     /// steering code or analyses in the finalize phase.
     size_t numEvents() const;
 
     /// @brief Access the sum of the event weights seen
     ///
     /// This is the weighted equivalent of the number of events. It should only
     /// be used by external steering code or analyses in the finalize phase.
     double sumW() const { return _eventCounter->sumW(); }
     /// Access to the sum of squared-weights
     double sumW2() const { return _eventCounter->sumW2(); }
 
     /// Names of event weight categories
     const vector<string>& weightNames() const {
         return _weightNames;
     }
 
 
     /// Get the index of the nominal weight-stream
     size_t defaultWeightIndex() const {
       return _defaultWeightIdx;
     }
 
 
     /// Set the cross-section for the process being generated.
     AnalysisHandler& setCrossSection(double xs, double xserr);
 
     /// Get the cross-section known to the handler.
     Scatter1DPtr crossSection() const { return _xs; }
 
     /// Get the nominal cross-section
     double nominalCrossSection() const {
       _xs.get()->setActiveWeightIdx(_defaultWeightIdx);
       const YODA::Scatter1D::Points& ps = _xs->points();
       if (ps.size() != 1) {
         string errMsg = "cross section missing when requesting nominal cross section";
         throw Error(errMsg);
       }
       double xs = ps[0].x();
       _xs.get()->unsetActiveWeight();
       return xs;
     }
 
 
     /// Set the beam particles for this run
     AnalysisHandler& setRunBeams(const ParticlePair& beams) {
       _beams = beams;
       MSG_DEBUG("Setting run beams = " << beams << " @ " << sqrtS()/GeV << " GeV");
       return *this;
     }
 
     /// Get the beam particles for this run, usually determined from the first event.
     const ParticlePair& beams() const { return _beams; }
 
     /// Get beam IDs for this run, usually determined from the first event.
     /// @deprecated Use standalone beamIds(ah.beams()), to clean AH interface
     PdgIdPair beamIds() const;
 
     /// Get energy for this run, usually determined from the first event.
     /// @deprecated Use standalone sqrtS(ah.beams()), to clean AH interface
     double sqrtS() const;
 
     /// Setter for _ignoreBeams
     void setIgnoreBeams(bool ignore=true);
 
     //@}
 
 
     /// @name Handle analyses
     //@{
 
     /// Get a list of the currently registered analyses' names.
     std::vector<std::string> analysisNames() const;
 
     /// Get the collection of currently registered analyses.
     const std::set<AnaHandle, CmpAnaHandle>& analyses() const {
       return _analyses;
     }
 
     /// Get a registered analysis by name.
     const AnaHandle analysis(const std::string& analysisname) const;
 
 
     /// Add an analysis to the run list by object
     AnalysisHandler& addAnalysis(Analysis* analysis);
 
     /// @brief Add an analysis to the run list using its name.
     ///
     /// The actual Analysis to be used will be obtained via
     /// AnalysisLoader::getAnalysis(string).  If no matching analysis is found,
     /// no analysis is added (i.e. the null pointer is checked and discarded.
     AnalysisHandler& addAnalysis(const std::string& analysisname);
 
     /// @brief Add an analysis with a map of analysis options.
     AnalysisHandler& addAnalysis(const std::string& analysisname, std::map<string, string> pars);
 
     /// @brief Add analyses to the run list using their names.
     ///
     /// The actual {@link Analysis}' to be used will be obtained via
     /// AnalysisHandler::addAnalysis(string), which in turn uses
     /// AnalysisLoader::getAnalysis(string). If no matching analysis is found
     /// for a given name, no analysis is added, but also no error is thrown.
     AnalysisHandler& addAnalyses(const std::vector<std::string>& analysisnames);
 
 
     /// Remove an analysis from the run list using its name.
     AnalysisHandler& removeAnalysis(const std::string& analysisname);
 
     /// Remove analyses from the run list using their names.
     AnalysisHandler& removeAnalyses(const std::vector<std::string>& analysisnames);
 
     //@}
 
 
     /// @name Main init/execute/finalise
     //@{
 
     /// Initialize a run, with the run beams taken from the example event.
     void init(const GenEvent& event);
 
     /// @brief Analyze the given \a event by reference.
     ///
     /// This function will call the AnalysisBase::analyze() function of all
     /// included analysis objects.
     void analyze(const GenEvent& event);
 
     /// @brief Analyze the given \a event by pointer.
     ///
     /// This function will call the AnalysisBase::analyze() function of all
     /// included analysis objects, after checking the event pointer validity.
     void analyze(const GenEvent* event);
 
     /// Finalize a run. This function calls the AnalysisBase::finalize()
     /// functions of all included analysis objects.
     void finalize();
 
     //@}
 
 
     /// @name Histogram / data object access
     //@{
 
     /// Add a vector of analysis objects to the current state.
     void addData(const std::vector<YODA::AnalysisObjectPtr>& aos);
 
     /// Read analysis plots into the histo collection (via addData) from the named file.
     void readData(const std::string& filename);
 
     /// Get all multi-weight Rivet analysis object wrappers
     vector<MultiweightAOPtr> getRivetAOs() const;
 
     /// Get single-weight YODA analysis objects
     vector<YODA::AnalysisObjectPtr> getYodaAOs(bool includeorphans,
                                                bool includetmps,
                                                bool usefinalized) const;
 
     /// Get all analyses' plots as a vector of analysis objects.
     std::vector<YODA::AnalysisObjectPtr> getData(bool includeorphans = false,
                                                  bool includetmps = false,
                                                  bool usefinalized = true) const;
 
     /// Write all analyses' plots (via getData) to the named file.
     void writeData(const std::string& filename) const;
 
-
     /// Tell Rivet to dump intermediate result to a file named @a
     /// dumpfile every @a period'th event. If @period is not positive,
     /// no dumping will be done.
     void dump(string dumpfile, int period) {
       _dumpPeriod = period;
       _dumpFile = dumpfile;
     }
 
     /// Take the vector of yoda files and merge them together using
     /// the cross section and weight information provided in each
     /// file. Each file in @a aofiles is assumed to have been produced
     /// by Rivet. By default the files are assumed to contain
     /// different processes (or the same processs but mutually
     /// exclusive cuts), but if @a equiv if ture, the files are
     /// assumed to contain output of completely equivalent (but
     /// statistically independent) Rivet runs. The corresponding
     /// analyses will be loaded and their analysis objects will be
     /// filled with the merged result. finalize() will be run on each
     /// relevant analysis. The resulting YODA file can then be rwitten
     /// out by writeData(). If delopts is non-empty, it is assumed to
     /// contain names different options to be merged into the same
     /// analysis objects.
     void mergeYodas(const vector<string> & aofiles,
                     const vector<string> & delopts = vector<string>(),
                     bool equiv = false);
 
     /// Helper function to strip specific options from data object paths.
     void stripOptions(YODA::AnalysisObjectPtr ao,
                       const vector<string> & delopts) const;
 
     //@}
 
+
+    /// @name Event weight access
+    //@{
+
+    /// Are any of the weights non-numeric?
+    size_t numWeights() { return _weightNames.size(); }
+
+    /// Are any of the weights non-numeric?
+    bool haveNamedWeights();
+
+    /// Set the weight names from a GenEvent
+    void setWeightNames(const GenEvent& ge);
+
+    //@}
+
+
     /// Indicate which Rivet stage we're in.
     /// At the moment, only INIT is used to enable booking.
     enum class Stage { OTHER, INIT };
 
     /// Which stage are we in?
     Stage stage() const { return _stage; }
 
 
   private:
 
     /// Current handler stage
     Stage _stage = Stage::OTHER;
 
     /// The collection of Analysis objects to be used.
     set<AnaHandle, CmpAnaHandle> _analyses;
 
     /// A vector of pre-loaded object which do not have a valid
     /// Analysis plugged in.
     vector<YODA::AnalysisObjectPtr> _orphanedPreloads;
 
     /// A vector containing copies of analysis objects after
     /// finalize() has been run.
     vector<YODA::AnalysisObjectPtr> _finalizedAOs;
 
+
     /// @name Run properties
     //@{
 
     /// Weight names
     std::vector<std::string> _weightNames;
     std::vector<std::valarray<double> > _subEventWeights;
     size_t _numWeightTypes; // always == WeightVector.size()
 
     /// Run name
     std::string _runname;
 
     /// Event counter
     mutable CounterPtr _eventCounter;
 
     /// Cross-section known to AH
     Scatter1DPtr _xs;
 
     /// Beams used by this run.
     ParticlePair _beams;
 
     /// Flag to check if init has been called
     bool _initialised;
 
     /// Flag whether input event beams should be ignored in compatibility check
     bool _ignoreBeams;
 
     /// Current event number
     int _eventNumber;
 
     /// The index in the weight vector for the nominal weight stream
     size_t _defaultWeightIdx;
 
     /// Determines how often Rivet runs finalize() and writes the
     /// result to a YODA file.
     int _dumpPeriod;
 
     /// The name of a YODA file to which Rivet periodically dumps
     /// results.
     string _dumpFile;
 
     /// Flag to indicate periodic dumping is in progress
     bool _dumping;
 
     //@}
 
 
   private:
 
     /// The assignment operator is private and must never be called.
     /// In fact, it should not even be implemented.
     AnalysisHandler& operator=(const AnalysisHandler&);
 
     /// The copy constructor is private and must never be called.  In
     /// fact, it should not even be implemented.
     AnalysisHandler(const AnalysisHandler&);
 
   };
 
 
 }
 
 #endif
diff --git a/include/Rivet/Projections/PercentileProjection.hh b/include/Rivet/Projections/PercentileProjection.hh
--- a/include/Rivet/Projections/PercentileProjection.hh
+++ b/include/Rivet/Projections/PercentileProjection.hh
@@ -1,138 +1,139 @@
 // -*- C++ -*-
 #ifndef RIVET_PERCENTILEPROJECTION_HH
 #define RIVET_PERCENTILEPROJECTION_HH
 
 #include "Rivet/Projections/SingleValueProjection.hh"
 #include "Rivet/Tools/RivetYODA.hh"
 #include <map>
 
 namespace Rivet {
 
 /**
    @brief class for projections that reports the percentile for a
    given SingleValueProjection when initialized with a Histo1D of the
    distribution in the SingleValueProjection.
 
    @author Leif Lönnblad
 
 */
 class PercentileProjection : public SingleValueProjection {
-
 public:
 
-  // Constructor taking a SingleValueProjection and a calibration
-  // histogram. If increasing it means that low values corresponds to
-  // lower percentiles.
-  PercentileProjection(const SingleValueProjection & sv, Histo1DPtr calhist,
+  /// Constructor taking a SingleValueProjection and a calibration
+  /// histogram. If increasing it means that low values corresponds to
+  /// lower percentiles.
+  ///
+  /// @todo Use mkScatter to pass this to the Scatter2D-calibrated version?
+  PercentileProjection(const SingleValueProjection & sv, const Histo1D& calhist,
                   bool increasing = false)
     : _calhist("EMPTY"), _increasing(increasing) {
     declare(sv, "OBSERVABLE");
-    if ( !calhist ) return;
-    MSG_INFO("Constructing PercentileProjection from " << calhist->path());
-    _calhist = calhist->path();
-    int N = calhist->numBins();
-    double sum = calhist->sumW();
+    //if ( !calhist ) return;
+    MSG_INFO("Constructing PercentileProjection from " << calhist.path());
+    _calhist = calhist.path();
+    int N = calhist.numBins();
+    double sum = calhist.sumW();
 
     if ( increasing ) {
-      double acc = calhist->underflow().sumW();
-      _table.insert(make_pair(calhist->bin(0).xEdges().first, 100.0*acc/sum));
+      double acc = calhist.underflow().sumW();
+      _table.insert(make_pair(calhist.bin(0).xEdges().first, 100.0*acc/sum));
       for ( int i = 0; i < N; ++i ) {
-        acc += calhist->bin(i).sumW();
-        _table.insert(make_pair(calhist->bin(i).xEdges().second, 100.0*acc/sum));
+        acc += calhist.bin(i).sumW();
+        _table.insert(make_pair(calhist.bin(i).xEdges().second, 100.0*acc/sum));
       }
     } else {
-      double acc = calhist->overflow().sumW();
-      _table.insert(make_pair(calhist->bin(N - 1).xEdges().second, 100.0*acc/sum));
+      double acc = calhist.overflow().sumW();
+      _table.insert(make_pair(calhist.bin(N - 1).xEdges().second, 100.0*acc/sum));
       for ( int i = N - 1; i >= 0; --i ) {
-        acc += calhist->bin(i).sumW();
-        _table.insert(make_pair(calhist->bin(i).xEdges().first, 100.0*acc/sum));
+        acc += calhist.bin(i).sumW();
+        _table.insert(make_pair(calhist.bin(i).xEdges().first, 100.0*acc/sum));
       }
     }
   }
 
   // Constructor taking a SingleValueProjection and a calibration
   // histogram. If increasing it means that low values corresponds to
   // lower percentiles.
-  PercentileProjection(const SingleValueProjection & sv, Scatter2DPtr calscat,
+  PercentileProjection(const SingleValueProjection & sv, const Scatter2D& calscat,
                   bool increasing = false)
     : _calhist("EMPTY"), _increasing(increasing) {
     declare(sv, "OBSERVABLE");
 
-    if ( !calscat ) return;
-    MSG_INFO("Constructing PercentileProjection from " << calscat->path());
-    _calhist = calscat->path();
-    int N = calscat->numPoints();
+    //if ( !calscat ) return;
+    MSG_INFO("Constructing PercentileProjection from " << calscat.path());
+    _calhist = calscat.path();
+    int N = calscat.numPoints();
     double sum = 0.0;
-    for ( const auto & p : calscat->points() ) sum += p.y();
+    for ( const auto & p : calscat.points() ) sum += p.y();
 
     double acc = 0.0;
     if ( increasing ) {
-      _table.insert(make_pair(calscat->point(0).xMin(), 100.0*acc/sum));
+      _table.insert(make_pair(calscat.point(0).xMin(), 100.0*acc/sum));
       for ( int i = 0; i < N; ++i ) {
-        acc += calscat->point(i).y();
-        _table.insert(make_pair(calscat->point(i).xMax(), 100.0*acc/sum));
+        acc += calscat.point(i).y();
+        _table.insert(make_pair(calscat.point(i).xMax(), 100.0*acc/sum));
       }
     } else {
-      _table.insert(make_pair(calscat->point(N - 1).xMax(), 100.0*acc/sum));
+      _table.insert(make_pair(calscat.point(N - 1).xMax(), 100.0*acc/sum));
       for ( int i = N - 1; i >= 0; --i ) {
-        acc += calscat->point(i).y();
-        _table.insert(make_pair(calscat->point(i).xMin(), 100.0*acc/sum));
+        acc += calscat.point(i).y();
+        _table.insert(make_pair(calscat.point(i).xMin(), 100.0*acc/sum));
       }
     }
   }
 
   DEFAULT_RIVET_PROJ_CLONE(PercentileProjection);
 
   // The projection function takes the assigned SingeValueProjection
   // and sets the value of this projection to the corresponding
   // percentile. If no calibration has been provided, -1 will be
   // returned. If values are outside of the calibration histogram, 0
   // or 100 will be returned.
   void project(const Event& e) {
     clear();
     if ( _table.empty() ) return;
     double obs = apply<SingleValueProjection>(e, "OBSERVABLE")();
     double pcnt = lookup(obs);
     if ( pcnt >= 0.0 ) set(pcnt);
   }
 
 
   // Standard comparison function.
   CmpState compare(const Projection& p) const {
     const PercentileProjection pp = dynamic_cast<const PercentileProjection&>(p);
     return mkNamedPCmp(p, "OBSERVABLE") ||
       cmp(_increasing, pp._increasing) ||
       cmp(_calhist, pp._calhist);
   }
 
 private:
 
   // The (interpolated) lookup table
   double lookup(double obs) const {
     auto low = _table.upper_bound(obs);
     if ( low == _table.end() ) return _increasing? 100.0: 0.0;
     if ( low == _table.begin() ) return _increasing? 0.0: 100.0;
     auto high = low--;
     return low->second + (obs - low->first)*(high->second - low->second)/
       (high->first - low->first);
   }
 
   // Astring identifying the calibration histogram.
   string _calhist;
 
   // A lookup table to find (by interpolation) the percentile given
   // the value of the underlying SingleValueProjection.
   map<double,double> _table;
 
   // A flag to say whether the distribution should be integrated from
   // below or above.
   bool _increasing;
 
 };
 
 
 
 
 }
 
 #endif
diff --git a/include/Rivet/Tools/RivetYODA.hh b/include/Rivet/Tools/RivetYODA.hh
--- a/include/Rivet/Tools/RivetYODA.hh
+++ b/include/Rivet/Tools/RivetYODA.hh
@@ -1,566 +1,578 @@
 #ifndef RIVET_RIVETYODA_HH
 #define RIVET_RIVETYODA_HH
 
 #include "Rivet/Config/RivetCommon.hh"
 #include "YODA/AnalysisObject.h"
 #include "YODA/Counter.h"
 #include "YODA/Histo1D.h"
 #include "YODA/Histo2D.h"
 #include "YODA/Profile1D.h"
 #include "YODA/Profile2D.h"
 #include "YODA/Scatter1D.h"
 #include "YODA/Scatter2D.h"
 #include "YODA/Scatter3D.h"
 
 #include <map>
 #include <valarray>
 
 namespace YODA {
   typedef std::shared_ptr<YODA::AnalysisObject> AnalysisObjectPtr;
+
+  typedef std::shared_ptr<YODA::Counter> CounterPtr;
+  typedef std::shared_ptr<YODA::Histo1D> Histo1DPtr;
+  typedef std::shared_ptr<YODA::Histo2D> Histo2DPtr;
+  typedef std::shared_ptr<YODA::Profile1D> Profile1DPtr;
+  typedef std::shared_ptr<YODA::Profile2D> Profile2DPtr;
+  typedef std::shared_ptr<YODA::Scatter1D> Scatter1DPtr;
+  typedef std::shared_ptr<YODA::Scatter2D> Scatter2DPtr;
+  typedef std::shared_ptr<YODA::Scatter3D> Scatter3DPtr;
+
 }
 
 
 namespace Rivet {
 
 
   class AnalysisObjectWrapper {
   public:
     virtual ~AnalysisObjectWrapper() {}
 
     virtual YODA::AnalysisObject* operator->() = 0;
     virtual YODA::AnalysisObject* operator->() const = 0;
     virtual const YODA::AnalysisObject & operator*() const = 0;
 
     /// @todo Rename to setActive(idx)
     virtual void setActiveWeightIdx(unsigned int iWeight) = 0;
 
     bool operator ==(const AnalysisObjectWrapper& p) { return (this == &p); }
     bool operator !=(const AnalysisObjectWrapper& p) { return (this != &p); }
 
   protected:
     /// @todo do we need this?
     // virtual void reset() = 0;
   };
 
   /// @todo
   /// implement scatter1dptr and scatter2dptr here
   /// these need to be multi-weighted eventually.
   /*
     class Scatter1DPtr : public AnalysisObjectPtr {
     public:
     Scatter1DPtr() : _persistent() { }
 
     Scatter1DPtr(size_t len_of_weightvec, const YODA::Scatter1D& p) {
     for (size_t m = 0; m < len_of_weightvec; ++m)
     _persistent.push_back(make_shared<YODA::Scatter1D>(p));
     }
 
     bool operator!() const { return !_persistent; }
     explicit operator bool() const { return bool(_persistent); }
 
     YODA::Scatter1D* operator->() { return _persistent.get(); }
 
     YODA::Scatter1D* operator->() const { return _persistent.get(); }
 
     YODA::Scatter1D & operator*() { return *_persistent; }
 
     const YODA::Scatter1D & operator*() const { return *_persistent; }
 
     protected:
     vector<YODA::Scatter1DPtr> _persistent;
     };
 
     class Scatter2DPtr : public AnalysisObjectPtr {
     public:
     Scatter2DPtr(size_t len_of_weightvec, const YODA::Scatter2D& p) {
     for (size_t m = 0; m < len_of_weightvec; ++m)
     _persistent.push_back(make_shared<YODA::Scatter2D>(p));
     }
 
     Scatter2DPtr() : _persistent() { }
 
     bool operator!() { return !_persistent; }
     explicit operator bool() { return bool(_persistent); }
 
     YODA::Scatter2D* operator->() { return _persistent.get(); }
 
     YODA::Scatter2D* operator->() const { return _persistent.get(); }
 
     YODA::Scatter2D & operator*() { return *_persistent; }
 
     const YODA::Scatter2D & operator*() const { return *_persistent; }
 
     protected:
     vector<YODA::Scatter2DPtr> _persistent;
     };
 
     class Scatter3DPtr : public AnalysisObjectPtr {
     public:
     Scatter3DPtr(size_t len_of_weightvec, const YODA::Scatter3D& p) {
     for (size_t m = 0; m < len_of_weightvec; ++m)
     _persistent.push_back(make_shared<YODA::Scatter3D>(p));
     }
 
     Scatter3DPtr() : _persistent() { }
 
     bool operator!() { return !_persistent; }
     explicit operator bool() { return bool(_persistent); }
 
     YODA::Scatter3D* operator->() { return _persistent.get(); }
 
     YODA::Scatter3D* operator->() const { return _persistent.get(); }
 
     YODA::Scatter3D & operator*() { return *_persistent; }
 
     const YODA::Scatter3D & operator*() const { return *_persistent; }
 
     protected:
     vector<YODA::Scatter3DPtr> _persistent;
     };
   */
 
 
   class MultiweightAOWrapper : public AnalysisObjectWrapper {
 
   public:
     using Inner = YODA::AnalysisObject;
 
     virtual void newSubEvent() = 0;
 
     virtual void pushToPersistent(const vector<std::valarray<double> >& weight) = 0;
 
     virtual YODA::AnalysisObjectPtr activeYODAPtr() const = 0;
   };
 
 
   using Weight = double;
 
   template <class T>
   using Fill = pair<typename T::FillType, Weight>;
 
   template <class T>
   using Fills = multiset<Fill<T>>;
 
 
   // TODO TODO TODO
   // need to override the old fill method too!
   // otherwise we can't intercept existing fill calls in analysis code
   // TODO TODO TODO
 
 
   /// Wrappers for analysis objects to store all fills unaggregated, until collapsed
   template <class T>
   class TupleWrapper;
 
   template<>
   class TupleWrapper<YODA::Counter> : public YODA::Counter {
   public:
     typedef shared_ptr<TupleWrapper<YODA::Counter>> Ptr;
     TupleWrapper(const YODA::Counter & h) : YODA::Counter(h) {}
     // todo: do we need to deal with users using fractions directly?
     void fill( double weight=1.0, double fraction=1.0 ) {
       fills_.insert( {YODA::Counter::FillType(),weight} );
     }
     void reset() { fills_.clear(); }
     const Fills<YODA::Counter> & fills() const { return fills_; }
   private:
     // x / weight pairs
     Fills<YODA::Counter> fills_;
   };
 
   template<>
   class TupleWrapper<YODA::Histo1D> : public YODA::Histo1D {
   public:
     typedef shared_ptr<TupleWrapper<YODA::Histo1D>> Ptr;
     TupleWrapper(const YODA::Histo1D & h) : YODA::Histo1D(h) {}
     // todo: do we need to deal with users using fractions directly?
     void fill( double x, double weight=1.0, double fraction=1.0 ) {
       if ( std::isnan(x) ) throw YODA::RangeError("X is NaN");
       fills_.insert( { x , weight } );
     }
     void reset() { fills_.clear(); }
     const Fills<YODA::Histo1D> & fills() const { return fills_; }
   private:
     // x / weight pairs
     Fills<YODA::Histo1D> fills_;
   };
 
   template<>
   class TupleWrapper<YODA::Profile1D> : public YODA::Profile1D {
   public:
     typedef shared_ptr<TupleWrapper<YODA::Profile1D>> Ptr;
     TupleWrapper(const YODA::Profile1D & h) : YODA::Profile1D(h) {}
     // todo: do we need to deal with users using fractions directly?
     void fill( double x, double y, double weight=1.0, double fraction=1.0 ) {
       if ( std::isnan(x) ) throw YODA::RangeError("X is NaN");
       if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN");
       fills_.insert( { YODA::Profile1D::FillType{x,y}, weight } );
     }
     void reset() { fills_.clear(); }
     const Fills<YODA::Profile1D> & fills() const { return fills_; }
   private:
     // x / weight pairs
     Fills<YODA::Profile1D> fills_;
   };
 
 
   template<>
   class TupleWrapper<YODA::Histo2D> : public YODA::Histo2D {
   public:
     typedef shared_ptr<TupleWrapper<YODA::Histo2D>> Ptr;
     TupleWrapper(const YODA::Histo2D & h) : YODA::Histo2D(h) {}
     // todo: do we need to deal with users using fractions directly?
     void fill( double x, double y, double weight=1.0, double fraction=1.0 ) {
       if ( std::isnan(x) ) throw YODA::RangeError("X is NaN");
       if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN");
       fills_.insert( { YODA::Histo2D::FillType{x,y}, weight } );
     }
     void reset() { fills_.clear(); }
     const Fills<YODA::Histo2D> & fills() const { return fills_; }
   private:
     // x / weight pairs
     Fills<YODA::Histo2D> fills_;
   };
 
   template<>
   class TupleWrapper<YODA::Profile2D> : public YODA::Profile2D {
   public:
     typedef shared_ptr<TupleWrapper<YODA::Profile2D>> Ptr;
     TupleWrapper(const YODA::Profile2D & h) : YODA::Profile2D(h) {}
     // todo: do we need to deal with users using fractions directly?
     void fill( double x, double y, double z, double weight=1.0, double fraction=1.0 ) {
       if ( std::isnan(x) ) throw YODA::RangeError("X is NaN");
       if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN");
       if ( std::isnan(z) ) throw YODA::RangeError("Z is NaN");
       fills_.insert( { YODA::Profile2D::FillType{x,y,z}, weight } );
     }
     void reset() { fills_.clear(); }
     const Fills<YODA::Profile2D> & fills() const { return fills_; }
   private:
     // x / weight pairs
     Fills<YODA::Profile2D> fills_;
   };
 
   template<>
   class TupleWrapper<YODA::Scatter1D> : public YODA::Scatter1D {
   public:
     typedef shared_ptr<TupleWrapper<YODA::Scatter1D>> Ptr;
     TupleWrapper(const YODA::Scatter1D & h) : YODA::Scatter1D(h) {}
   };
 
   template<>
   class TupleWrapper<YODA::Scatter2D> : public YODA::Scatter2D {
   public:
     typedef shared_ptr<TupleWrapper<YODA::Scatter2D>> Ptr;
     TupleWrapper(const YODA::Scatter2D & h) : YODA::Scatter2D(h) {}
   };
 
   template<>
   class TupleWrapper<YODA::Scatter3D> : public YODA::Scatter3D {
   public:
     typedef shared_ptr<TupleWrapper<YODA::Scatter3D>> Ptr;
     TupleWrapper(const YODA::Scatter3D & h) : YODA::Scatter3D(h) {}
   };
 
 
 
   template <class T>
   class Wrapper : public MultiweightAOWrapper {
     friend class Analysis;
 
   public:
 
     using Inner = T;
     /* @todo
      * some things are not really well-defined here
      * for instance: fill() in the finalize() method and integral() in
      * the analyze() method.
      */
 
     Wrapper() = default;
 
     Wrapper(const vector<string>& weightnames, const T & p);
 
     ~Wrapper();
 
     typename T::Ptr active() const;
 
     /* @todo this probably need to loop over all? */
     bool operator!() const { return !_active; } // Don't use active() here, assert will catch
 
     explicit operator bool() const { return static_cast<bool>(_active); } // Don't use active() here, assert will catch
 
     T * operator->() { return active().get(); }
 
     T * operator->() const { return active().get(); }
 
     T & operator*() { return *active(); }
 
     const T & operator*() const { return *active(); }
     /* @todo
      * these need to be re-thought out.
 
      void reset() { active()->reset(); }
     */
 
     /* @todo
      * these probably need to loop over all?
      * do we even want to provide equality?
      */
     /* @todo
      * how about no.
      friend bool operator==(Wrapper a, Wrapper 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;
      }
 
      friend bool operator!=(Wrapper a, Wrapper b){
      return !(a == b);
      }
 
 
      friend bool operator<(Wrapper a, Wrapper 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;
      }
     */
 
 
   private:
     void setActiveWeightIdx(unsigned int iWeight) {
       _active = _persistent.at(iWeight);
     }
 
     /* this is for dev only---we shouldn't need this in real runs. */
     void unsetActiveWeight() { _active.reset(); }
 
     void newSubEvent();
 
     virtual YODA::AnalysisObjectPtr activeYODAPtr() const { return _active; }
 
     const vector<typename T::Ptr> & persistent() const { return _persistent; }
 
     /* to be implemented for each type */
     void pushToPersistent(const vector<std::valarray<double> >& weight);
 
     /* M of these, one for each weight */
     vector<typename T::Ptr> _persistent;
 
     /* N of these, one for each event in evgroup */
     vector<typename TupleWrapper<T>::Ptr> _evgroup;
 
     typename T::Ptr _active;
 
     // do we need implicit cast?
     // operator typename T::Ptr () {
     //     return _active;
     // }
 
     friend class AnalysisHandler;
   };
 
 
   /// We need our own shared_ptr class, so we can dispatch -> and *
   /// all the way down to the inner YODA analysis objects
   ///
   /// TODO: provide remaining functionality that shared_ptr has (not needed right now)
   ///
   template <typename T>
   class rivet_shared_ptr {
   public:
     typedef T value_type;
 
     rivet_shared_ptr() = default;
 
     rivet_shared_ptr(decltype(nullptr)) : _p(nullptr) {}
 
     /// Convenience constructor, pass through to the Wrapper constructor
     rivet_shared_ptr(const vector<string>& weightNames, const typename T::Inner & p)
       : _p( make_shared<T>(weightNames, p) )
     {}
 
     template <typename U>
     rivet_shared_ptr(const shared_ptr<U> & p)
       : _p(p)
     {}
 
     template <typename U>
     rivet_shared_ptr(const rivet_shared_ptr<U> & p)
       : _p(p.get())
     {}
 
     // Goes right through to the active YODA object's members
     T & operator->()                            { return  *_p; }
     const T & operator->() const                { return  *_p; }
 
     // The active YODA object
     typename T::Inner & operator*()             { return **_p; }
     const typename T::Inner & operator*() const { return **_p; }
 
     bool operator!() const { return !_p || !(*_p);   }
     explicit operator bool()  const { return _p && bool(*_p); }
 
     template <typename U>
     bool operator==(const rivet_shared_ptr<U> & other) const {
       return _p == other._p;
     }
 
     template <typename U>
     bool operator!=(const rivet_shared_ptr<U> & other) const {
       return _p != other._p;
     }
 
     template <typename U>
     bool operator<(const rivet_shared_ptr<U> & other) const {
       return _p < other._p;
     }
 
     template <typename U>
     bool operator>(const rivet_shared_ptr<U> & other) const {
       return _p > other._p;
     }
 
     template <typename U>
     bool operator<=(const rivet_shared_ptr<U> & other) const {
       return _p <= other._p;
     }
 
     template <typename U>
     bool operator>=(const rivet_shared_ptr<U> & other) const {
       return _p >= other._p;
     }
 
     shared_ptr<T> get() const { return _p; }
+
   private:
     shared_ptr<T> _p;
+
   };
 
 
 
   // every object listed here needs a virtual fill method in YODA,
   // otherwise the Tuple fakery won't work.
 
   using MultiweightAOPtr = rivet_shared_ptr<MultiweightAOWrapper>;
 
   using Histo1DPtr   = rivet_shared_ptr<Wrapper<YODA::Histo1D>>;
   using Histo2DPtr   = rivet_shared_ptr<Wrapper<YODA::Histo2D>>;
   using Profile1DPtr = rivet_shared_ptr<Wrapper<YODA::Profile1D>>;
   using Profile2DPtr = rivet_shared_ptr<Wrapper<YODA::Profile2D>>;
   using CounterPtr   = rivet_shared_ptr<Wrapper<YODA::Counter>>;
   using Scatter1DPtr = rivet_shared_ptr<Wrapper<YODA::Scatter1D>>;
   using Scatter2DPtr = rivet_shared_ptr<Wrapper<YODA::Scatter2D>>;
   using Scatter3DPtr = rivet_shared_ptr<Wrapper<YODA::Scatter3D>>;
 
   using YODA::Counter;
   using YODA::Histo1D;
   using YODA::HistoBin1D;
   using YODA::Histo2D;
   using YODA::HistoBin2D;
   using YODA::Profile1D;
   using YODA::ProfileBin1D;
   using YODA::Profile2D;
   using YODA::ProfileBin2D;
   using YODA::Scatter1D;
   using YODA::Point1D;
   using YODA::Scatter2D;
   using YODA::Point2D;
   using YODA::Scatter3D;
   using YODA::Point3D;
 
   /// Function to get a map of all the refdata in a paper with the
   /// given @a papername.
   map<string, YODA::AnalysisObjectPtr> getRefData(const string& papername);
 
   /// @todo Also provide a Scatter3D getRefData() version?
 
   /// Get the file system path to the reference file for this paper.
   string getDatafilePath(const string& papername);
 
 
   /// Traits class to access the type of the AnalysisObject in the
   /// reference files.
   template<typename T> struct ReferenceTraits {};
   template<> struct ReferenceTraits<Counter> { typedef Counter RefT; };
   template<> struct ReferenceTraits<Scatter1D> { typedef Scatter1D RefT; };
   template<> struct ReferenceTraits<Histo1D> { typedef Scatter2D RefT; };
   template<> struct ReferenceTraits<Profile1D> { typedef Scatter2D RefT; };
   template<> struct ReferenceTraits<Scatter2D> { typedef Scatter2D RefT; };
   template<> struct ReferenceTraits<Histo2D> { typedef Scatter3D RefT; };
   template<> struct ReferenceTraits<Profile2D> { typedef Scatter3D RefT; };
   template<> struct ReferenceTraits<Scatter3D> { typedef Scatter3D RefT; };
 
 
   /// If @a dst and @a src both are of same subclass T, copy the
   /// contents of @a src into @a dst and return true. Otherwise return
   /// false.
   template <typename T>
   inline bool aocopy(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst) {
     shared_ptr<T> tsrc = dynamic_pointer_cast<T>(src);
     if ( !tsrc ) return false;
     shared_ptr<T> tdst = dynamic_pointer_cast<T>(dst);
     if ( !tdst ) return false;
     *tdst = *tsrc;
     return true;
   }
 
   /// If @a dst and @a src both are of same subclass T, add the
   /// contents of @a src into @a dst and return true. Otherwise return
   /// false.
   template <typename T>
   inline bool aoadd(YODA::AnalysisObjectPtr dst, YODA::AnalysisObjectPtr src, double scale) {
     shared_ptr<T> tsrc = dynamic_pointer_cast<T>(src);
     if ( !tsrc ) return false;
     shared_ptr<T> tdst = dynamic_pointer_cast<T>(dst);
     if ( !tdst ) return false;
     tsrc->scaleW(scale);
     *tdst += *tsrc;
     return true;
   }
 
   /// If @a dst is the same subclass as @a src, copy the contents of @a
   /// src into @a dst and return true. Otherwise return false.
   bool copyao(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst);
 
   /// If @a dst is the same subclass as @a src, scale the contents of
   /// @a src with @a scale and add it to @a dst and return true. Otherwise
   /// return false.
   bool addaos(YODA::AnalysisObjectPtr dst, YODA::AnalysisObjectPtr src, double scale);
 
   /// Check if two analysis objects have the same binning or, if not
   /// binned, are in other ways compatible.
   template <typename TPtr>
   inline bool bookingCompatible(TPtr a, TPtr b) {
     return a->sameBinning(*b);
   }
   inline bool bookingCompatible(CounterPtr, CounterPtr) {
     return true;
   }
   inline bool bookingCompatible(Scatter1DPtr a, Scatter1DPtr b) {
     return a->numPoints() == b->numPoints();
   }
   inline bool bookingCompatible(Scatter2DPtr a, Scatter2DPtr b) {
     return a->numPoints() == b->numPoints();
   }
   inline bool bookingCompatible(Scatter3DPtr a, Scatter3DPtr b) {
     return a->numPoints() == b->numPoints();
   }
 
 }
 
 #endif
diff --git a/src/Core/Analysis.cc b/src/Core/Analysis.cc
--- a/src/Core/Analysis.cc
+++ b/src/Core/Analysis.cc
@@ -1,1048 +1,1050 @@
 // -*- C++ -*-
 #include "Rivet/Config/RivetCommon.hh"
 #include "Rivet/Analysis.hh"
 #include "Rivet/AnalysisHandler.hh"
 #include "Rivet/AnalysisInfo.hh"
 #include "Rivet/Tools/BeamConstraint.hh"
 #include "Rivet/Projections/ImpactParameterProjection.hh"
 #include "Rivet/Projections/GeneratedPercentileProjection.hh"
 #include "Rivet/Projections/UserCentEstimate.hh"
 #include "Rivet/Projections/CentralityProjection.hh"
 
 namespace Rivet {
 
 
   Analysis::Analysis(const string& name)
       : _analysishandler(NULL)
   {
     ProjectionApplier::_allowProjReg = false;
     _defaultname = name;
 
     unique_ptr<AnalysisInfo> ai = AnalysisInfo::make(name);
     assert(ai);
     _info = move(ai);
     assert(_info);
   }
 
   double Analysis::sqrtS() const {
     return handler().sqrtS();
   }
 
   const ParticlePair& Analysis::beams() const {
     return handler().beams();
   }
 
   const PdgIdPair Analysis::beamIds() const {
     return handler().beamIds();
   }
 
 
   const string Analysis::histoDir() const {
     /// @todo Cache in a member variable
     string _histoDir;
     if (_histoDir.empty()) {
       _histoDir = "/" + name();
       if (handler().runName().length() > 0) {
         _histoDir = "/" + handler().runName() + _histoDir;
       }
       replace_all(_histoDir, "//", "/"); //< iterates until none
     }
     return _histoDir;
   }
 
 
   const string Analysis::histoPath(const string& hname) const {
     const string path = histoDir() + "/" + hname;
     return path;
   }
 
 
   const string Analysis::histoPath(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
     return histoDir() + "/" + mkAxisCode(datasetId, xAxisId, yAxisId);
   }
 
 
   const string Analysis::mkAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
     std::stringstream axisCode;
     axisCode << "d";
     if (datasetId < 10) axisCode << 0;
     axisCode << datasetId;
     axisCode << "-x";
     if (xAxisId < 10) axisCode << 0;
     axisCode << xAxisId;
     axisCode << "-y";
     if (yAxisId < 10) axisCode << 0;
     axisCode << yAxisId;
     return axisCode.str();
   }
 
 
   Log& Analysis::getLog() const {
     string logname = "Rivet.Analysis." + name();
     return Log::getLog(logname);
   }
 
 
   ///////////////////////////////////////////
 
 
   size_t Analysis::numEvents() const {
     return handler().numEvents();
   }
 
   double Analysis::sumW() const {
     return handler().sumW();
   }
 
   double Analysis::sumW2() const {
     return handler().sumW2();
   }
 
 
   ///////////////////////////////////////////
 
 
   bool Analysis::isCompatible(const ParticlePair& beams) const {
     return isCompatible(beams.first.pid(),  beams.second.pid(),
                         beams.first.energy(), beams.second.energy());
   }
 
 
   bool Analysis::isCompatible(PdgId beam1, PdgId beam2, double e1, double e2) const {
     PdgIdPair beams(beam1, beam2);
     pair<double,double> energies(e1, e2);
     return isCompatible(beams, energies);
   }
 
 
   // bool Analysis::beamIDsCompatible(const PdgIdPair& beams) const {
   //   bool beamIdsOk = false;
   //   for (const PdgIdPair& bp : requiredBeams()) {
   //     if (compatible(beams, bp)) {
   //       beamIdsOk =  true;
   //       break;
   //     }
   //   }
   //   return beamIdsOk;
   // }
 
 
   // /// Check that the energies are compatible (within 1% or 1 GeV, whichever is larger, for a bit of UI forgiveness)
   // bool Analysis::beamEnergiesCompatible(const pair<double,double>& energies) const {
   //   /// @todo Use some sort of standard ordering to improve comparisons, esp. when the two beams are different particles
   //   bool beamEnergiesOk = requiredEnergies().size() > 0 ? false : true;
   //   typedef pair<double,double> DoublePair;
   //   for (const DoublePair& ep : requiredEnergies()) {
   //     if ((fuzzyEquals(ep.first, energies.first, 0.01) && fuzzyEquals(ep.second, energies.second, 0.01)) ||
   //         (fuzzyEquals(ep.first, energies.second, 0.01) && fuzzyEquals(ep.second, energies.first, 0.01)) ||
   //         (abs(ep.first - energies.first) < 1*GeV && abs(ep.second - energies.second) < 1*GeV) ||
   //         (abs(ep.first - energies.second) < 1*GeV && abs(ep.second - energies.first) < 1*GeV)) {
   //       beamEnergiesOk =  true;
   //       break;
   //     }
   //   }
   //   return beamEnergiesOk;
   // }
 
 
   // bool Analysis::beamsCompatible(const PdgIdPair& beams, const pair<double,double>& energies) const {
   bool Analysis::isCompatible(const PdgIdPair& beams, const pair<double,double>& energies) const {
     // First check the beam IDs
     bool beamIdsOk = false;
     for (const PdgIdPair& bp : requiredBeams()) {
       if (compatible(beams, bp)) {
         beamIdsOk =  true;
         break;
       }
     }
     if (!beamIdsOk) return false;
 
     // Next check that the energies are compatible (within 1% or 1 GeV, whichever is larger, for a bit of UI forgiveness)
 
     /// @todo Use some sort of standard ordering to improve comparisons, esp. when the two beams are different particles
     bool beamEnergiesOk = requiredEnergies().size() > 0 ? false : true;
     typedef pair<double,double> DoublePair;
     for (const DoublePair& ep : requiredEnergies()) {
       if ((fuzzyEquals(ep.first, energies.first, 0.01) && fuzzyEquals(ep.second, energies.second, 0.01)) ||
           (fuzzyEquals(ep.first, energies.second, 0.01) && fuzzyEquals(ep.second, energies.first, 0.01)) ||
           (abs(ep.first - energies.first) < 1*GeV && abs(ep.second - energies.second) < 1*GeV) ||
           (abs(ep.first - energies.second) < 1*GeV && abs(ep.second - energies.first) < 1*GeV)) {
         beamEnergiesOk =  true;
         break;
       }
     }
     return beamEnergiesOk;
   }
 
 
   ///////////////////////////////////////////
 
   double Analysis::crossSection() const {
     const YODA::Scatter1D::Points& ps = handler().crossSection()->points();
     if (ps.size() != 1) {
       string errMsg = "cross section missing for analysis " + name();
       throw Error(errMsg);
     }
     return ps[0].x();
   }
 
   double Analysis::crossSectionPerEvent() const {
     return crossSection()/sumW();
   }
 
 
 
   ////////////////////////////////////////////////////////////
   // Histogramming
 
 
   void Analysis::_cacheRefData() const {
     if (_refdata.empty()) {
       MSG_TRACE("Getting refdata cache for paper " << name());
       _refdata = getRefData(getRefDataName());
     }
   }
 
   // vector<YODA::AnalysisObjectPtr> Analysis::getAllData(bool includeorphans) const{
   //   return handler().getData(includeorphans);
   // }
 
 
   CounterPtr & Analysis::book(CounterPtr & ctr,
                       const string& cname,
                       const string& title) {
     const string path = histoPath(cname);
     ctr = CounterPtr(handler().weightNames(), Counter(path, title));
     ctr = addAnalysisObject(ctr);
     return ctr;
   }
 
 
   CounterPtr & Analysis::book(CounterPtr & ctr, unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                               const string& title) {
     const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId);
     return book(ctr, axisCode, title);
   }
 
 
 
 
   Histo1DPtr & Analysis::book(Histo1DPtr & histo, const string& hname,
                               size_t nbins, double lower, double upper,
                               const string& title,
                               const string& xtitle,
                               const string& ytitle) {
     const string path = histoPath(hname);
 
     Histo1D hist = Histo1D(nbins, lower, upper, path, title);
     hist.setAnnotation("XLabel", xtitle);
     hist.setAnnotation("YLabel", ytitle);
 
     histo = Histo1DPtr(handler().weightNames(), hist);
     histo = addAnalysisObject(histo);
     return histo;
   }
 
   Histo1DPtr & Analysis::book(Histo1DPtr & histo, const string& hname,
                               const initializer_list<double>& binedges,
                               const string& title,
                               const string& xtitle,
                               const string& ytitle) {
   	return book(histo, hname, vector<double>{binedges}, title, xtitle, ytitle);
   }
 
   Histo1DPtr & Analysis::book(Histo1DPtr & histo, const string& hname,
                               const vector<double>& binedges,
                               const string& title,
                               const string& xtitle,
                               const string& ytitle) {
     const string path = histoPath(hname);
 
     Histo1D hist = Histo1D(binedges, path, title);
     hist.setAnnotation("XLabel", xtitle);
     hist.setAnnotation("YLabel", ytitle);
 
     histo = Histo1DPtr(handler().weightNames(), hist);
     histo = addAnalysisObject(histo);
     return histo;
   }
 
   Histo1DPtr & Analysis::book(Histo1DPtr & histo, const string& hname,
                               const string& title,
                               const string& xtitle,
                               const string& ytitle) {
     const Scatter2D& refdata = refData(hname);
     return book(histo, hname, refdata, title, xtitle, ytitle);
   }
 
 
   Histo1DPtr & Analysis::book(Histo1DPtr & histo, unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                               const string& title,
                               const string& xtitle,
                               const string& ytitle) {
     const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId);
     return book(histo, axisCode, title, xtitle, ytitle);
   }
 
   Histo1DPtr & Analysis::book(Histo1DPtr& histo,
                               const string& hname,
                               const Scatter2D& refscatter,
                               const string& title,
                               const string& xtitle,
                               const string& ytitle) {
     const string path = histoPath(hname);
 
     Histo1D hist = Histo1D(refscatter, path);
     hist.setTitle(title);
     hist.setAnnotation("XLabel", xtitle);
     hist.setAnnotation("YLabel", ytitle);
     if (hist.hasAnnotation("IsRef")) hist.rmAnnotation("IsRef");
 
     histo = Histo1DPtr(handler().weightNames(), hist);
     histo = addAnalysisObject(histo);
     return histo;
   }
 
 
   /////////////////
 
 
   Histo2DPtr & Analysis::book(Histo2DPtr & h2d,const string& hname,
                                    size_t nxbins, double xlower, double xupper,
                                    size_t nybins, double ylower, double yupper,
                                    const string& title,
                                    const string& xtitle,
                                    const string& ytitle,
                                    const string& ztitle)
   {
     const string path = histoPath(hname);
 
     Histo2D hist(nxbins, xlower, xupper, nybins, ylower, yupper, path, title);
     hist.setAnnotation("XLabel", xtitle);
     hist.setAnnotation("YLabel", ytitle);
     hist.setAnnotation("ZLabel", ztitle);
 
     h2d = Histo2DPtr(handler().weightNames(), hist);
     h2d = addAnalysisObject(h2d);
     return h2d;
   }
 
   Histo2DPtr & Analysis::book(Histo2DPtr & h2d,const string& hname,
                                    const initializer_list<double>& xbinedges,
                                    const initializer_list<double>& ybinedges,
                                    const string& title,
                                    const string& xtitle,
                                    const string& ytitle,
                                    const string& ztitle)
   {
   	return book(h2d, hname, vector<double>{xbinedges}, vector<double>{ybinedges}, title, xtitle, ytitle, ztitle);
   }
 
   Histo2DPtr & Analysis::book(Histo2DPtr & h2d,const string& hname,
                                    const vector<double>& xbinedges,
                                    const vector<double>& ybinedges,
                                    const string& title,
                                    const string& xtitle,
                                    const string& ytitle,
                                    const string& ztitle)
   {
     const string path = histoPath(hname);
 
     Histo2D hist(xbinedges, ybinedges, path, title);
     hist.setAnnotation("XLabel", xtitle);
     hist.setAnnotation("YLabel", ytitle);
     hist.setAnnotation("ZLabel", ztitle);
 
     h2d = Histo2DPtr(handler().weightNames(), hist);
     h2d = addAnalysisObject(h2d);
     return h2d;
   }
 
 
   Histo2DPtr & Analysis::book(Histo2DPtr & histo, const string& hname,
                                    const Scatter3D& refscatter,
                                    const string& title,
                                    const string& xtitle,
                                    const string& ytitle,
                                    const string& ztitle) {
     const string path = histoPath(hname);
 
     Histo2D hist = Histo2D(refscatter, path);
     hist.setTitle(title);
     hist.setAnnotation("XLabel", xtitle);
     hist.setAnnotation("YLabel", ytitle);
     hist.setAnnotation("ZLabel", ztitle);
     if (hist.hasAnnotation("IsRef")) hist.rmAnnotation("IsRef");
 
     histo = Histo2DPtr(handler().weightNames(), hist);
     histo = addAnalysisObject(histo);
     return histo;
   }
 
 
   Histo2DPtr & Analysis::book(Histo2DPtr & histo, const string& hname,
                                    const string& title,
                                    const string& xtitle,
                                    const string& ytitle,
                                    const string& ztitle) {
     const Scatter3D& refdata = refData<Scatter3D>(hname);
     return book(histo, hname, refdata, title, xtitle, ytitle, ztitle);
   }
 
 
   Histo2DPtr & Analysis::book(Histo2DPtr & histo, unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                                    const string& title,
                                    const string& xtitle,
                                    const string& ytitle,
                                    const string& ztitle) {
     const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId);
     return book(histo, axisCode, title, xtitle, ytitle, ztitle);
   }
 
 
   /////////////////
 
 
   Profile1DPtr & Analysis::book(Profile1DPtr & p1d,const string& hname,
                                        size_t nbins, double lower, double upper,
                                        const string& title,
                                        const string& xtitle,
                                        const string& ytitle) {
     const string path = histoPath(hname);
 
     Profile1D prof(nbins, lower, upper, path, title);
     prof.setAnnotation("XLabel", xtitle);
     prof.setAnnotation("YLabel", ytitle);
 
     p1d = Profile1DPtr(handler().weightNames(), prof);
     p1d = addAnalysisObject(p1d);
     return p1d;
   }
 
 
   Profile1DPtr & Analysis::book(Profile1DPtr & p1d,const string& hname,
                                        const initializer_list<double>& binedges,
                                        const string& title,
                                        const string& xtitle,
                                        const string& ytitle) {
   	return book(p1d, hname, vector<double>{binedges}, title, xtitle, ytitle);
   }
 
   Profile1DPtr & Analysis::book(Profile1DPtr & p1d, const string& hname,
                                        const vector<double>& binedges,
                                        const string& title,
                                        const string& xtitle,
                                        const string& ytitle) {
     const string path = histoPath(hname);
 
     Profile1D prof(binedges, path, title);
     prof.setAnnotation("XLabel", xtitle);
     prof.setAnnotation("YLabel", ytitle);
 
     p1d = Profile1DPtr(handler().weightNames(), prof);
     p1d = addAnalysisObject(p1d);
     return p1d;
   }
 
   Profile1DPtr & Analysis::book(Profile1DPtr & p1d, const string& hname,
                                        const Scatter2D& refscatter,
                                        const string& title,
                                        const string& xtitle,
                                        const string& ytitle) {
     const string path = histoPath(hname);
 
     Profile1D prof(refscatter, path);
     prof.setTitle(title);
     prof.setAnnotation("XLabel", xtitle);
     prof.setAnnotation("YLabel", ytitle);
     if (prof.hasAnnotation("IsRef")) prof.rmAnnotation("IsRef");
 
     p1d = Profile1DPtr(handler().weightNames(), prof);
     p1d = addAnalysisObject(p1d);
     return p1d;
   }
 
 
   Profile1DPtr & Analysis::book(Profile1DPtr & p1d,const string& hname,
                                        const string& title,
                                        const string& xtitle,
                                        const string& ytitle) {
     const Scatter2D& refdata = refData(hname);
     book(p1d, hname, refdata, title, xtitle, ytitle);
     return p1d;
   }
 
 
   Profile1DPtr & Analysis::book(Profile1DPtr & p1d,unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                                        const string& title,
                                        const string& xtitle,
                                        const string& ytitle) {
     const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId);
     return book(p1d, axisCode, title, xtitle, ytitle);
   }
 
 
   ///////////////////
 
 
   Profile2DPtr & Analysis::book(Profile2DPtr & p2d, const string& hname,
                                    size_t nxbins, double xlower, double xupper,
                                    size_t nybins, double ylower, double yupper,
                                    const string& title,
                                    const string& xtitle,
                                    const string& ytitle,
                                    const string& ztitle)
   {
     const string path = histoPath(hname);
 
     Profile2D prof(nxbins, xlower, xupper, nybins, ylower, yupper, path, title);
     prof.setAnnotation("XLabel", xtitle);
     prof.setAnnotation("YLabel", ytitle);
     prof.setAnnotation("ZLabel", ztitle);
 
     p2d = Profile2DPtr(handler().weightNames(), prof);
     p2d = addAnalysisObject(p2d);
     return p2d;
   }
 
 
   Profile2DPtr & Analysis::book(Profile2DPtr & p2d, const string& hname,
                                    const initializer_list<double>& xbinedges,
                                    const initializer_list<double>& ybinedges,
                                    const string& title,
                                    const string& xtitle,
                                    const string& ytitle,
                                    const string& ztitle)
   {
   	return book(p2d, hname, vector<double>{xbinedges}, vector<double>{ybinedges}, title, xtitle, ytitle, ztitle);
   }
 
 
   Profile2DPtr & Analysis::book(Profile2DPtr & p2d, const string& hname,
                                    const vector<double>& xbinedges,
                                    const vector<double>& ybinedges,
                                    const string& title,
                                    const string& xtitle,
                                    const string& ytitle,
                                    const string& ztitle)
   {
     const string path = histoPath(hname);
 
     Profile2D prof(xbinedges, ybinedges, path, title);
     prof.setAnnotation("XLabel", xtitle);
     prof.setAnnotation("YLabel", ytitle);
     prof.setAnnotation("ZLabel", ztitle);
 
     p2d = Profile2DPtr(handler().weightNames(), prof);
     p2d = addAnalysisObject(p2d);
     return p2d;
   }
 
 
   /// @todo REINSTATE
 
   // Profile2DPtr Analysis::book(Profile2DPtr& prof,const string& hname,
   //                                      const Scatter3D& refscatter,
   //                                      const string& title,
   //                                      const string& xtitle,
   //                                      const string& ytitle,
   //                                      const string& ztitle) {
   //   const string path = histoPath(hname);
 
   //   /// @todo Add no-metadata argument to YODA copy constructors
   //   Profile2D prof(refscatter, path);
   //   prof.setTitle(title);
   //   prof.setAnnotation("XLabel", xtitle);
   //   prof.setAnnotation("YLabel", ytitle);
   //   prof.setAnnotation("ZLabel", ztitle);
   //   if (prof.hasAnnotation("IsRef")) prof.rmAnnotation("IsRef");
 
   //   p2d = Profile2DPtr(handler().weightNames(), prof);
   //   p2d = addAnalysisObject(p2d);
   //   return p2d;
   // }
 
 
   // Profile2DPtr Analysis::book(Profile2DPtr& prof, const string& hname,
   //                                      const string& title,
   //                                      const string& xtitle,
   //                                      const string& ytitle,
   //                                      const string& ztitle) {
   //   const Scatter3D& refdata = refData<Scatter3D>(hname);
   //   return book(prof, hname, refdata, title, xtitle, ytitle, ztitle);
   // }
 
 
   ///////////////
 
 
   /// @todo Should be able to book Scatter1Ds
 
 
   ///////////////
 
 
   Scatter2DPtr & Analysis::book(Scatter2DPtr & s2d, unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                                        bool copy_pts,
                                        const string& title,
                                        const string& xtitle,
                                        const string& ytitle) {
     const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId);
     return book(s2d, axisCode, copy_pts, title, xtitle, ytitle);
   }
 
 
   Scatter2DPtr & Analysis::book(Scatter2DPtr & s2d, const string& hname,
                                        bool copy_pts,
                                        const string& title,
                                        const string& xtitle,
                                        const string& ytitle) {
     const string path = histoPath(hname);
 
     Scatter2D scat;
     if (copy_pts) {
       const Scatter2D& refdata = refData(hname);
       scat = Scatter2D(refdata, path);
       for (Point2D& p : scat.points()) p.setY(0, 0);
     } else {
       scat = Scatter2D(path);
     }
 
     scat.setTitle(title);
     scat.setAnnotation("XLabel", xtitle);
     scat.setAnnotation("YLabel", ytitle);
     if (scat.hasAnnotation("IsRef")) scat.rmAnnotation("IsRef");
 
     s2d = Scatter2DPtr(handler().weightNames(), scat);
     s2d = addAnalysisObject(s2d);
     return s2d;
   }
 
 
   Scatter2DPtr & Analysis::book(Scatter2DPtr & s2d, const string& hname,
                                        size_t npts, double lower, double upper,
                                        const string& title,
                                        const string& xtitle,
                                        const string& ytitle) {
     const string path = histoPath(hname);
 
     Scatter2D scat;
     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);
     }
 
     scat.setTitle(title);
     scat.setAnnotation("XLabel", xtitle);
     scat.setAnnotation("YLabel", ytitle);
 
     s2d = Scatter2DPtr(handler().weightNames(), scat);
     s2d = addAnalysisObject(s2d);
     return s2d;
   }
 
   Scatter2DPtr & Analysis::book(Scatter2DPtr & s2d, const string& hname,
                                        const vector<double>& binedges,
                                        const string& title,
                                        const string& xtitle,
                                        const string& ytitle) {
     const string path = histoPath(hname);
 
     Scatter2D scat;
     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);
     }
 
     scat.setTitle(title);
     scat.setAnnotation("XLabel", xtitle);
     scat.setAnnotation("YLabel", ytitle);
 
     s2d = Scatter2DPtr(handler().weightNames(), scat);
     s2d = addAnalysisObject(s2d);
     return s2d;
   }
 
 
   /// @todo Book Scatter3Ds?
 
 
   /////////////////////
 
 
   void Analysis::divide(CounterPtr c1, CounterPtr c2, Scatter1DPtr s) const {
     const string path = s->path();
     *s = *c1 / *c2;
     s->setPath(path);
   }
 
   void Analysis::divide(const Counter& c1, const Counter& c2, Scatter1DPtr s) const {
     const string path = s->path();
     *s = c1 / c2;
     s->setPath(path);
   }
 
 
   void Analysis::divide(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const {
     const string path = s->path();
     *s = *h1 / *h2;
     s->setPath(path);
   }
 
   void Analysis::divide(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const {
     const string path = s->path();
     *s = h1 / h2;
     s->setPath(path);
   }
 
 
   void Analysis::divide(Profile1DPtr p1, Profile1DPtr p2, Scatter2DPtr s) const {
     const string path = s->path();
     *s = *p1 / *p2;
     s->setPath(path);
   }
 
   void Analysis::divide(const Profile1D& p1, const Profile1D& p2, Scatter2DPtr s) const {
     const string path = s->path();
     *s = p1 / p2;
     s->setPath(path);
   }
 
 
   void Analysis::divide(Histo2DPtr h1, Histo2DPtr h2, Scatter3DPtr s) const {
     const string path = s->path();
     *s = *h1 / *h2;
     s->setPath(path);
   }
 
   void Analysis::divide(const Histo2D& h1, const Histo2D& h2, Scatter3DPtr s) const {
     const string path = s->path();
     *s = h1 / h2;
     s->setPath(path);
   }
 
 
   void Analysis::divide(Profile2DPtr p1, Profile2DPtr p2, Scatter3DPtr s) const {
     const string path = s->path();
     *s = *p1 / *p2;
     s->setPath(path);
   }
 
   void Analysis::divide(const Profile2D& p1, const Profile2D& p2, Scatter3DPtr s) const {
     const string path = s->path();
     *s = p1 / p2;
     s->setPath(path);
   }
 
 
   /// @todo Counter and Histo2D efficiencies and asymms
 
 
   void Analysis::efficiency(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const {
     const string path = s->path();
     *s = YODA::efficiency(*h1, *h2);
     s->setPath(path);
   }
 
   void Analysis::efficiency(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const {
     const string path = s->path();
     *s = YODA::efficiency(h1, h2);
     s->setPath(path);
   }
 
 
   void Analysis::asymm(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const {
     const string path = s->path();
     *s = YODA::asymm(*h1, *h2);
     s->setPath(path);
   }
 
   void Analysis::asymm(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const {
     const string path = s->path();
     *s = YODA::asymm(h1, h2);
     s->setPath(path);
   }
 
 
   void Analysis::scale(CounterPtr cnt, Analysis::CounterAdapter factor) {
     if (!cnt) {
       MSG_WARNING("Failed to scale counter=NULL in analysis " << name() << " (scale=" << double(factor) << ")");
       return;
     }
     if (std::isnan(double(factor)) || std::isinf(double(factor))) {
       MSG_WARNING("Failed to scale counter=" << cnt->path() << " in analysis: " << name() << " (invalid scale factor = " << double(factor) << ")");
       factor = 0;
     }
     MSG_TRACE("Scaling counter " << cnt->path() << " by factor " << double(factor));
     try {
       cnt->scaleW(factor);
     } catch (YODA::Exception& we) {
       MSG_WARNING("Could not scale counter " << cnt->path());
       return;
     }
   }
 
 
   void Analysis::normalize(Histo1DPtr histo, Analysis::CounterAdapter norm, bool includeoverflows) {
     if (!histo) {
       MSG_WARNING("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << double(norm) << ")");
       return;
     }
     MSG_TRACE("Normalizing histo " << histo->path() << " to " << double(norm));
     try {
       const double hint = histo->integral(includeoverflows);
       if (hint == 0)  MSG_WARNING("Skipping histo with null area " << histo->path());
       else            histo->normalize(norm, includeoverflows);
     } catch (YODA::Exception& we) {
       MSG_WARNING("Could not normalize histo " << histo->path());
       return;
     }
   }
 
 
   void Analysis::scale(Histo1DPtr histo, Analysis::CounterAdapter factor) {
     if (!histo) {
       MSG_WARNING("Failed to scale histo=NULL in analysis " << name() << " (scale=" << double(factor) << ")");
       return;
     }
     if (std::isnan(double(factor)) || std::isinf(double(factor))) {
       MSG_WARNING("Failed to scale histo=" << histo->path() << " in analysis: " << name() << " (invalid scale factor = " << double(factor) << ")");
       factor = 0;
     }
     MSG_TRACE("Scaling histo " << histo->path() << " by factor " << double(factor));
     try {
       histo->scaleW(factor);
     } catch (YODA::Exception& we) {
       MSG_WARNING("Could not scale histo " << histo->path());
       return;
     }
   }
 
 
   void Analysis::normalize(Histo2DPtr histo, Analysis::CounterAdapter norm, bool includeoverflows) {
     if (!histo) {
       MSG_ERROR("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << double(norm) << ")");
       return;
     }
     MSG_TRACE("Normalizing histo " << histo->path() << " to " << double(norm));
     try {
       const double hint = histo->integral(includeoverflows);
       if (hint == 0)  MSG_WARNING("Skipping histo with null area " << histo->path());
       else            histo->normalize(norm, includeoverflows);
     } catch (YODA::Exception& we) {
       MSG_WARNING("Could not normalize histo " << histo->path());
       return;
     }
   }
 
 
   void Analysis::scale(Histo2DPtr histo, Analysis::CounterAdapter factor) {
     if (!histo) {
       MSG_ERROR("Failed to scale histo=NULL in analysis " << name() << " (scale=" << double(factor) << ")");
       return;
     }
     if (std::isnan(double(factor)) || std::isinf(double(factor))) {
       MSG_ERROR("Failed to scale histo=" << histo->path() << " in analysis: " << name() << " (invalid scale factor = " << double(factor) << ")");
       factor = 0;
     }
     MSG_TRACE("Scaling histo " << histo->path() << " by factor " << double(factor));
     try {
       histo->scaleW(factor);
     } catch (YODA::Exception& we) {
       MSG_WARNING("Could not scale histo " << histo->path());
       return;
     }
   }
 
 
   void Analysis::integrate(Histo1DPtr h, Scatter2DPtr s) const {
     // preserve the path info
     const string path = s->path();
     *s = toIntegralHisto(*h);
     s->setPath(path);
   }
 
   void Analysis::integrate(const Histo1D& h, Scatter2DPtr s) const {
     // preserve the path info
     const string path = s->path();
     *s = toIntegralHisto(h);
     s->setPath(path);
   }
 
 }
   /// @todo 2D versions of integrate... defined how, exactly?!?
 
 
   //////////////////////////////////
 
-namespace {
-  void errormsg(std::string name) {
-// #ifdef HAVE_BACKTRACE
-//      void * buffer[4];
-//      backtrace(buffer, 4);
-//      backtrace_symbols_fd(buffer, 4 , 1);
-// #endif
-    std::cerr << name << ": Can't book objects outside of init().\n";
-    assert(false);
-  }
-}
+// namespace {
+//   void errormsg(std::string name) {
+//     // #ifdef HAVE_BACKTRACE
+//     //      void * buffer[4];
+//     //      backtrace(buffer, 4);
+//     //      backtrace_symbols_fd(buffer, 4 , 1);
+//     // #endif
+//     std::cerr << name << ": Can't book objects outside of init().\n";
+//     assert(false);
+//   }
+// }
 
 
 namespace Rivet {
 
 
   // void Analysis::addAnalysisObject(const MultiweightAOPtr & ao) {
   //   if (handler().stage() == AnalysisHandler::Stage::INIT) {
   //     _analysisobjects.push_back(ao);
   //   }
   //   else {
   //     errormsg(name());
   //   }
   // }
 
 
   void Analysis::removeAnalysisObject(const string& path) {
     for (auto it = _analysisobjects.begin();
          it != _analysisobjects.end(); ++it) {
       if ((*it)->path() == path) {
         _analysisobjects.erase(it);
         break;
       }
     }
   }
 
   void Analysis::removeAnalysisObject(const MultiweightAOPtr & ao) {
     for (auto it = _analysisobjects.begin();  it != _analysisobjects.end(); ++it) {
       if ((*it) == ao) {
         _analysisobjects.erase(it);
         break;
       }
     }
  }
 
 const CentralityProjection &
 Analysis::declareCentrality(const SingleValueProjection &proj,
                             string calAnaName, string calHistName,
                             const string projName, bool increasing) {
 
   CentralityProjection cproj;
 
   // Select the centrality variable from option. Use REF as default.
   // Other selections are "GEN", "IMP" and "USR" (USR only in HEPMC 3).
   string sel = getOption<string>("cent","REF");
   set<string> done;
 
   if ( sel == "REF" ) {
-    Scatter2DPtr refscat;
+    YODA::Scatter2DPtr refscat;
     auto refmap = getRefData(calAnaName);
     if ( refmap.find(calHistName) != refmap.end() )
-      refscat =
-        dynamic_pointer_cast<Scatter2D>(refmap.find(calHistName)->second);
+      refscat = dynamic_pointer_cast<Scatter2D>(refmap.find(calHistName)->second);
 
     if ( !refscat ) {
       MSG_WARNING("No reference calibration histogram for " <<
                   "CentralityProjection " << projName << " found " <<
                   "(requested histogram " << calHistName << " in " <<
                   calAnaName << ")");
     }
     else {
       MSG_INFO("Found calibration histogram " << sel << " " << refscat->path());
-      cproj.add(PercentileProjection(proj, refscat, increasing), sel);
+      cproj.add(PercentileProjection(proj, *refscat, increasing), sel);
     }
   }
   else if ( sel == "GEN" ) {
-    Histo1DPtr genhist;
-    string histpath = "/" + calAnaName + "/" + calHistName;
-    for ( YODA::AnalysisObjectPtr ao : handler().getData(true) ) {
-      if ( ao->path() == histpath )
-        genhist = dynamic_pointer_cast<Histo1D>(ao);
-    }
-    if ( !genhist || genhist->numEntries() <= 1 ) {
+    Histo1DPtr genhists =
+      getAnalysisObject<Histo1DPtr>(calAnaName, calHistName + "_IMP");
+    // for ( YODA::AnalysisObjectPtr ao : handler().getData(true) ) {
+    //   if ( ao->path() == histpath )
+    //     genhist = dynamic_pointer_cast<Histo1D>(ao);
+    // }
+    if ( !genhists || genhists->numEntries() <= 1 ) {
       MSG_WARNING("No generated calibration histogram for " <<
                "CentralityProjection " << projName << " found " <<
                "(requested histogram " << calHistName << " in " <<
                calAnaName << ")");
     }
     else {
-      MSG_INFO("Found calibration histogram " << sel << " " << genhist->path());
-      cproj.add(PercentileProjection(proj, genhist, increasing), sel);
+      MSG_INFO("Found calibration histogram " << sel << " " << genhists->path());
+      genhists.get()->setActiveWeightIdx(_defaultWeightIndex());
+      cproj.add(PercentileProjection(proj, *genhists, increasing), sel);
     }
   }
   else if ( sel == "IMP" ) {
-    Histo1DPtr imphist =
+    Histo1DPtr imphists =
       getAnalysisObject<Histo1DPtr>(calAnaName, calHistName + "_IMP");
-    if ( !imphist || imphist->numEntries() <= 1 ) {
+    if ( !imphists || imphists->numEntries() <= 1 ) {
       MSG_WARNING("No impact parameter calibration histogram for " <<
                "CentralityProjection " << projName << " found " <<
                "(requested histogram " << calHistName << "_IMP in " <<
                calAnaName << ")");
     }
     else {
-      MSG_INFO("Found calibration histogram " << sel << " " << imphist->path());
+      MSG_INFO("Found calibration histogram " << sel << " " << imphists->path());
+      imphists.get()->setActiveWeightIdx(_defaultWeightIndex());
       cproj.add(PercentileProjection(ImpactParameterProjection(),
-                                     imphist, true), sel);
+                                     *imphists, true), sel);
     }
   }
   else if ( sel == "USR" ) {
 #if HEPMC_VERSION_CODE >= 3000000
-    Histo1DPtr usrhist =
+    Histo1DPtr usrhists =
       getAnalysisObject<Histo1DPtr>(calAnaName, calHistName + "_USR");
-    if ( !usrhist || usrhist->numEntries() <= 1 ) {
+    if ( !usrhists || usrhists->numEntries() <= 1 ) {
       MSG_WARNING("No user-defined calibration histogram for " <<
                "CentralityProjection " << projName << " found " <<
                "(requested histogram " << calHistName << "_USR in " <<
                   calAnaName << ")");
       continue;
     }
     else {
-      MSG_INFO("Found calibration histogram " << sel << " " << usrhist->path());
-      cproj.add((UserCentEstimate(), usrhist, true), sel);
+      MSG_INFO("Found calibration histogram " << sel << " " << usrhists->path());
+      usrhists.get()->setActiveWeightIdx(_defaultWeightIndex());
+      cproj.add((UserCentEstimate(), usrhists*, true), sel);
      }
 #else
       MSG_WARNING("UserCentEstimate is only available with HepMC3.");
 #endif
     }
   else if ( sel == "RAW" ) {
 #if HEPMC_VERSION_CODE >= 3000000
     cproj.add(GeneratedCentrality(), sel);
 #else
     MSG_WARNING("GeneratedCentrality is only available with HepMC3.");
 #endif
   }
     else
       MSG_WARNING("'" << sel << "' is not a valid PercentileProjection tag.");
 
   if ( cproj.empty() )
     MSG_WARNING("CentralityProjection " << projName
                 << " did not contain any valid PercentileProjections.");
 
   return declare(cproj, projName);
 
 }
 
 
   vector<string> Analysis::_weightNames() const {
     return handler().weightNames();
   }
 
   size_t Analysis::_defaultWeightIndex() const {
     return handler().defaultWeightIndex();
   }
 
   MultiweightAOPtr Analysis::_getOtherAnalysisObject(const std::string & ananame, const std::string& name) {
     std::string path = "/" + ananame + "/" + name;
     const auto& ana = handler().analysis(ananame);
     return ana->getAnalysisObject(name); //< @todo includeorphans check??
   }
 
   void Analysis::_checkBookInit() const {
     if (handler().stage() != AnalysisHandler::Stage::INIT) {
       MSG_ERROR("Can't book objects outside of init()");
       throw UserError(name() + ": Can't book objects outside of init().");
     }
   }
 
 
 }
diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc
--- a/src/Core/AnalysisHandler.cc
+++ b/src/Core/AnalysisHandler.cc
@@ -1,763 +1,755 @@
 // -*- C++ -*-
 #include "Rivet/Config/RivetCommon.hh"
 #include "Rivet/AnalysisHandler.hh"
 #include "Rivet/Analysis.hh"
 #include "Rivet/Tools/ParticleName.hh"
 #include "Rivet/Tools/BeamConstraint.hh"
 #include "Rivet/Tools/Logging.hh"
 #include "Rivet/Projections/Beam.hh"
-#include "YODA/ReaderYODA.h"
-#include "YODA/WriterYODA.h"
+#include "YODA/IO.h"
 #include <regex>
 #include <iostream>
 
 using std::cout;
 using std::cerr;
 
 namespace {
 
   inline std::vector<std::string> split(const std::string& input, const std::string& regex) {
     // passing -1 as the submatch index parameter performs splitting
     std::regex re(regex);
     std::sregex_token_iterator
       first{input.begin(), input.end(), re, -1},
       last;
       return {first, last};
   }
 
 }
 
 
 namespace Rivet {
 
 
   AnalysisHandler::AnalysisHandler(const string& runname)
     : _runname(runname),
       _eventCounter(vector<string>(), Counter()), _xs(vector<string>(), Scatter1D()),
       _initialised(false), _ignoreBeams(false),
       _defaultWeightIdx(0)
   {}
 
 
   AnalysisHandler::~AnalysisHandler()
   {  }
 
 
   Log& AnalysisHandler::getLog() const {
     return Log::getLog("Rivet.Analysis.Handler");
   }
 
 
   /// http://stackoverflow.com/questions/4654636/how-to-determine-if-a-string-is-a-number-with-c
-  bool is_number(const std::string& s)
-  {
+  bool is_number(const std::string& s) {
     std::string::const_iterator it = s.begin();
     while (it != s.end() && std::isdigit(*it)) ++it;
     return !s.empty() && it == s.end();
   }
 
   /// Check if any of the weightnames is not a number
   bool AnalysisHandler::haveNamedWeights() {
     bool dec=false;
     for (unsigned int i=0;i<_weightNames.size();++i) {
       string s = _weightNames[i];
       if (!is_number(s)) {
         dec=true;
         break;
       }
     }
     return dec;
   }
 
 
   void AnalysisHandler::init(const GenEvent& ge) {
     if (_initialised)
       throw UserError("AnalysisHandler::init has already been called: cannot re-initialize!");
 
     /// @todo Should the Rivet analysis objects know about weight names?
 
     setRunBeams(Rivet::beams(ge));
     MSG_DEBUG("Initialising the analysis handler");
     _eventNumber = ge.event_number();
 
     setWeightNames(ge);
     if (haveNamedWeights())
         MSG_INFO("Using named weights");
     else
         MSG_INFO("NOT using named weights. Using first weight as nominal weight");
 
     _eventCounter = CounterPtr(weightNames(), Counter("_EVTCOUNT"));
 
     // Set the cross section based on what is reported by this event.
     if (ge.cross_section()) {
       MSG_TRACE("Getting cross section.");
       double xs = ge.cross_section()->cross_section();
       double xserr = ge.cross_section()->cross_section_error();
       setCrossSection(xs, xserr);
     }
 
     // Check that analyses are beam-compatible, and remove those that aren't
     const size_t num_anas_requested = analysisNames().size();
     vector<string> anamestodelete;
     for (const AnaHandle a : _analyses) {
       if (!_ignoreBeams && !a->isCompatible(beams())) {
         //MSG_DEBUG(a->name() << " requires beams " << a->requiredBeams() << " @ " << a->requiredEnergies() << " GeV");
         anamestodelete.push_back(a->name());
       }
     }
     for (const string& aname : anamestodelete) {
       MSG_WARNING("Analysis '" << aname << "' is incompatible with the provided beams: removing");
       removeAnalysis(aname);
     }
     if (num_anas_requested > 0 && analysisNames().empty()) {
       cerr << "All analyses were incompatible with the first event's beams\n"
            << "Exiting, since this probably wasn't intentional!" << endl;
       exit(1);
     }
 
     // Warn if any analysis' status is not unblemished
     for (const AnaHandle a : analyses()) {
       if (toUpper(a->status()) == "PRELIMINARY") {
         MSG_WARNING("Analysis '" << a->name() << "' is preliminary: be careful, it may change and/or be renamed!");
       } else if (toUpper(a->status()) == "OBSOLETE") {
         MSG_WARNING("Analysis '" << a->name() << "' is obsolete: please update!");
       } else if (toUpper(a->status()).find("UNVALIDATED") != string::npos) {
         MSG_WARNING("Analysis '" << a->name() << "' is unvalidated: be careful, it may be broken!");
       }
     }
 
     // Initialize the remaining analyses
     _stage = Stage::INIT;
     for (AnaHandle a : _analyses) {
       MSG_DEBUG("Initialising analysis: " << a->name());
       try {
         // Allow projection registration in the init phase onwards
         a->_allowProjReg = true;
         a->init();
         //MSG_DEBUG("Checking consistency of analysis: " << a->name());
         //a->checkConsistency();
       } catch (const Error& err) {
         cerr << "Error in " << a->name() << "::init method: " << err.what() << endl;
         exit(1);
       }
       MSG_DEBUG("Done initialising analysis: " << a->name());
     }
     _stage = Stage::OTHER;
     _initialised = true;
     MSG_DEBUG("Analysis handler initialised");
   }
 
 
   void AnalysisHandler::setWeightNames(const GenEvent& ge) {
     /// reroute the print output to a std::stringstream and process
     /// The iteration is done over a map in hepmc2 so this is safe
     std::ostringstream stream;
     ge.weights().print(stream);  // Super lame, I know
     string str =  stream.str();
 
     std::regex re("(([^()]+))"); // Regex for stuff enclosed by parentheses ()
     size_t idx = 0;
     for (std::sregex_iterator i = std::sregex_iterator(str.begin(), str.end(), re); i != std::sregex_iterator(); ++i ) {
       std::smatch m = *i;
       vector<string> temp = ::split(m.str(), "[,]");
       if (temp.size() ==2) {
         MSG_DEBUG("Name of weight #" << _weightNames.size() << ": " << temp[0]);
 
         // store the default weight based on weight names
         if (temp[0] == "Weight" || temp[0] == "0" || temp[0] == "Default") {
           MSG_DEBUG(_weightNames.size() << " is being used as the nominal.");
           _weightNames.push_back("");
           _defaultWeightIdx = idx;
         } else
           _weightNames.push_back(temp[0]);
 
 
         idx++;
       }
     }
   }
 
 
   void AnalysisHandler::analyze(const GenEvent& ge) {
     // Call init with event as template if not already initialised
     if (!_initialised) init(ge);
     assert(_initialised);
 
     // Ensure that beam details match those from the first event (if we're checking beams)
     if ( !_ignoreBeams ) {
       const PdgIdPair beams = Rivet::beamIds(ge);
       const double sqrts = Rivet::sqrtS(ge);
       if (!compatible(beams, _beams) || !fuzzyEquals(sqrts, sqrtS())) {
         cerr << "Event beams mismatch: "
              << PID::toBeamsString(beams) << " @ " << sqrts/GeV << " GeV" << " vs. first beams "
              << this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << endl;
         exit(1);
       }
     }
 
-
     // Create the Rivet event wrapper
     /// @todo Filter/normalize the event here
     Event event(ge);
 
-    // set the cross section based on what is reported by this event.
-    // if no cross section
-    MSG_TRACE("getting cross section.");
+    // Set the cross section based on what is reported by this event.
+    MSG_TRACE("Getting cross-section.");
     if (ge.cross_section()) {
-      MSG_TRACE("getting cross section from GenEvent.");
+      MSG_TRACE("Getting cross-section from GenEvent.");
       double xs = ge.cross_section()->cross_section();
       double xserr = ge.cross_section()->cross_section_error();
       setCrossSection(xs, xserr);
     }
-    #endif
 
-    // won't happen for first event because _eventNumber is set in
-    // init()
+    // Won't happen for first event because _eventNumber is set in init()
     if (_eventNumber != ge.event_number()) {
-        /// @todo
-        /// can we get away with not passing a matrix?
+        /// @todo Can we get away with not passing a matrix?
 
         MSG_TRACE("AnalysisHandler::analyze(): Pushing _eventCounter to persistent.");
         _eventCounter.get()->pushToPersistent(_subEventWeights);
         // if this is indeed a new event, push the temporary
         // histograms and reset
         for (const AnaHandle& a : _analyses) {
             for (auto ao : a->analysisObjects()) {
                 MSG_TRACE("AnalysisHandler::analyze(): Pushing " << a->name() << "'s " << ao->name() << " to persistent.");
                 ao.get()->pushToPersistent(_subEventWeights);
             }
             MSG_TRACE("AnalysisHandler::analyze(): finished pushing " << a->name() << "'s objects to persistent.");
         }
 
         _eventNumber = ge.event_number();
 
         MSG_DEBUG("nominal event # " << _eventCounter.get()->_persistent[0]->numEntries());
         MSG_DEBUG("nominal sum of weights: " << _eventCounter.get()->_persistent[0]->sumW());
         MSG_DEBUG("Event has " << _subEventWeights.size() << " sub events.");
         _subEventWeights.clear();
     }
 
 
     MSG_TRACE("starting new sub event");
     _eventCounter.get()->newSubEvent();
 
     for (const AnaHandle& a : _analyses) {
         for (auto ao : a->analysisObjects()) {
             ao.get()->newSubEvent();
         }
     }
 
     _subEventWeights.push_back(event.weights());
     MSG_DEBUG("Analyzing subevent #" << _subEventWeights.size() - 1 << ".");
 
     _eventCounter->fill();
     // Run the analyses
     for (AnaHandle a : _analyses) {
       MSG_TRACE("About to run analysis " << a->name());
       try {
         a->analyze(event);
       } catch (const Error& err) {
         cerr << "Error in " << a->name() << "::analyze method: " << err.what() << endl;
         exit(1);
       }
       MSG_TRACE("Finished running analysis " << a->name());
     }
 
     if ( _dumpPeriod > 0 && numEvents()%_dumpPeriod == 0 ) {
       MSG_INFO("Dumping intermediate results to " << _dumpFile << ".");
       _dumping = true;
       finalize();
       _dumping = false;
       writeData(_dumpFile);
     }
 
   }
 
 
   void AnalysisHandler::analyze(const GenEvent* ge) {
     if (ge == nullptr) {
       MSG_ERROR("AnalysisHandler received null pointer to GenEvent");
       //throw Error("AnalysisHandler received null pointer to GenEvent");
     }
     analyze(*ge);
   }
 
 
   void AnalysisHandler::finalize() {
     if (!_initialised) return;
     MSG_INFO("Finalising analyses");
 
     // First push all analyses' objects to persistent
     MSG_TRACE("AnalysisHandler::finalize(): Pushing analysis objects to persistent.");
     _eventCounter.get()->pushToPersistent(_subEventWeights);
     for (const AnaHandle& a : _analyses) {
       for (auto ao : a->analysisObjects())
         ao.get()->pushToPersistent(_subEventWeights);
     }
 
     // First we make copies of all analysis objects.
-    map<string,AnalysisObjectPtr> backupAOs;
+    map<string,YODA::AnalysisObjectPtr> backupAOs;
     for (auto ao : getYodaAOs(false, true, false) )
-      backupAOs[ao->path()] = AnalysisObjectPtr(ao->newclone());
+      backupAOs[ao->path()] = YODA::AnalysisObjectPtr(ao->newclone());
 
     // Finalize all the histograms
     for (const AnaHandle& a : _analyses) {
       // a->setCrossSection(_xs);
       for (size_t iW = 0; iW < numWeights(); iW++) {
         _eventCounter.get()->setActiveWeightIdx(iW);
         _xs.get()->setActiveWeightIdx(iW);
         for (auto ao : a->analysisObjects())
           ao.get()->setActiveWeightIdx(iW);
 
         MSG_TRACE("Running " << a->name() << "::finalize() for weight " << iW << ".");
 
         try {
           if ( !_dumping || a->info().reentrant() )  a->finalize();
           else if ( _dumping == 1 && iW == 0 )
             MSG_INFO("Skipping periodic dump of " << a->name()
                      << " as it is not declared reentrant.");
         } catch (const Error& err) {
           cerr << "Error in " << a->name() << "::finalize method: " << err.what() << endl;
           exit(1);
         }
 
       }
     }
 
     // Now we copy all analysis objects to the list of finalized
     // ones, and restore the value to their original ones.
     _finalizedAOs.clear();
     for ( auto ao : getYodaAOs(false, false, false) )
-      _finalizedAOs.push_back(AnalysisObjectPtr(ao->newclone()));
+      _finalizedAOs.push_back(YODA::AnalysisObjectPtr(ao->newclone()));
     for ( auto ao : getYodaAOs(false, true, false) ) {
       // TODO: This should be possible to do in a nicer way, with a flag etc.
       if (ao->path().find("/FINAL") != std::string::npos) continue;
       auto aoit = backupAOs.find(ao->path());
       if ( aoit == backupAOs.end() ) {
         AnaHandle ana = analysis(split(ao->path(), "/")[0]);
         if ( ana ) ana->removeAnalysisObject(ao->path());
       } else
         copyao(aoit->second, ao);
     }
 
     // Print out number of events processed
     const int nevts = numEvents();
     MSG_INFO("Processed " << nevts << " event" << (nevts != 1 ? "s" : ""));
 
     // // Delete analyses
     // MSG_DEBUG("Deleting analyses");
     // _analyses.clear();
 
     // Print out MCnet boilerplate
     cout << '\n';
     cout << "The MCnet usage guidelines apply to Rivet: see http://www.montecarlonet.org/GUIDELINES" << '\n';
     cout << "Please acknowledge plots made with Rivet analyses, and cite arXiv:1003.0694 (http://arxiv.org/abs/1003.0694)" << endl;
   }
 
   AnalysisHandler& AnalysisHandler::addAnalysis(const string& analysisname, std::map<string, string> pars) {
      // Make an option handle.
     std::string parHandle = "";
     for (map<string, string>::iterator par = pars.begin(); par != pars.end(); ++par) {
       parHandle +=":";
       parHandle += par->first + "=" + par->second;
     }
     return addAnalysis(analysisname + parHandle);
 
   }
 
   AnalysisHandler& AnalysisHandler::addAnalysis(const string& analysisname) {
     // Check for a duplicate analysis
     /// @todo Might we want to be able to run an analysis twice, with different params?
     ///       Requires avoiding histo tree clashes, i.e. storing the histos on the analysis objects.
     string ananame = analysisname;
     vector<string> anaopt = split(analysisname, ":");
     if ( anaopt.size() > 1 ) ananame = anaopt[0];
     AnaHandle analysis( AnalysisLoader::getAnalysis(ananame) );
     if (analysis.get() != 0) { // < Check for null analysis.
       MSG_DEBUG("Adding analysis '" << analysisname << "'");
       map<string,string> opts;
       for ( int i = 1, N = anaopt.size(); i < N; ++i ) {
         vector<string> opt = split(anaopt[i], "=");
         if ( opt.size() != 2 ) {
           MSG_WARNING("Error in option specification. Skipping analysis "
                       << analysisname);
           return *this;
         }
         if ( !analysis->info().validOption(opt[0], opt[1]) ) {
           MSG_WARNING("Cannot set option '" << opt[0] << "' to '" << opt[1]
                       << "'. Skipping analysis " << analysisname);
           return *this;
         }
         opts[opt[0]] = opt[1];
       }
       for ( auto opt: opts) {
         analysis->_options[opt.first] = opt.second;
         analysis->_optstring += ":" + opt.first + "=" + opt.second;
       }
       for (const AnaHandle& a : _analyses) {
         if (a->name() == analysis->name() ) {
           MSG_WARNING("Analysis '" << analysisname << "' already registered: skipping duplicate");
           return *this;
         }
       }
       analysis->_analysishandler = this;
       _analyses.insert(analysis);
     } else {
       MSG_WARNING("Analysis '" << analysisname << "' not found.");
     }
     // MSG_WARNING(_analyses.size());
     // for (const AnaHandle& a : _analyses) MSG_WARNING(a->name());
     return *this;
   }
 
 
   AnalysisHandler& AnalysisHandler::removeAnalysis(const string& analysisname) {
     std::shared_ptr<Analysis> toremove;
     for (const AnaHandle a : _analyses) {
       if (a->name() == analysisname) {
         toremove = a;
         break;
       }
     }
     if (toremove.get() != 0) {
       MSG_DEBUG("Removing analysis '" << analysisname << "'");
       _analyses.erase(toremove);
     }
     return *this;
   }
 
 
   void AnalysisHandler::addData(const std::vector<YODA::AnalysisObjectPtr>& aos) {
     for (const YODA::AnalysisObjectPtr ao : aos) {
       string path = ao->path();
       if ( path.substr(0, 5) != "/RAW/" ) {
         _orphanedPreloads.push_back(ao);
         continue;
       }
 
       path = path.substr(4);
       ao->setPath(path);
       if (path.size() > 1) { // path > "/"
         try {
           const string ananame =  ::split(path, "/")[0];
           AnaHandle a = analysis(ananame);
+          /// @todo FIXXXXX
           //MultiweightAOPtr mao = ????; /// @todo generate right Multiweight object from ao
-          a->addAnalysisObject(mao); /// @todo Need to statistically merge...
+          //a->addAnalysisObject(mao); /// @todo Need to statistically merge...
         } catch (const Error& e) {
           MSG_TRACE("Adding analysis object " << path <<
                     " to the list of orphans.");
           _orphanedPreloads.push_back(ao);
         }
       }
     }
   }
 
-  void AnalysisHandler::stripOptions(AnalysisObjectPtr ao,
+
+  void AnalysisHandler::stripOptions(YODA::AnalysisObjectPtr ao,
                                      const vector<string> & delopts) const {
     string path = ao->path();
     string ananame = split(path, "/")[0];
     vector<string> anaopts = split(ananame, ":");
     for ( int i = 1, N = anaopts.size(); i < N; ++i )
       for ( auto opt : delopts )
         if ( opt == "*" || anaopts[i].find(opt + "=") == 0 )
           path.replace(path.find(":" + anaopts[i]), (":" + anaopts[i]).length(), "");
     ao->setPath(path);
   }
 
 
-
-
   void AnalysisHandler::
   mergeYodas(const vector<string> & aofiles, const vector<string> & delopts, bool equiv) {
-    vector< vector<AnalysisObjectPtr> > aosv;
+    vector< vector<YODA::AnalysisObjectPtr> > aosv;
     vector<double> xsecs;
     vector<double> xsecerrs;
     vector<CounterPtr> sows;
     set<string> ananames;
-     _eventcounter.reset();
+    _eventCounter->reset();
 
     // First scan all files and extract analysis objects and add the
-    // corresponding anayses..
-    for ( auto file : aofiles ) {
+    // corresponding analyses..
+    for (const string& file : aofiles ) {
       Scatter1DPtr xsec;
       CounterPtr sow;
 
       // For each file make sure that cross section and sum-of-weights
       // objects are present and stor all RAW ones in a vector;
-      vector<AnalysisObjectPtr> aos;
+      vector<YODA::AnalysisObjectPtr> aos;
       try {
         /// @todo Use new YODA SFINAE to fill the smart ptr vector directly
         vector<YODA::AnalysisObject*> aos_raw;
         YODA::read(file, aos_raw);
-        for (AnalysisObject* aor : aos_raw) {
-          AnalysisObjectPtr ao = AnalysisObjectPtr(aor);
+        for (YODA::AnalysisObject* aor : aos_raw) {
+          YODA::AnalysisObjectPtr ao(aor);
           if ( ao->path().substr(0, 5) != "/RAW/" ) continue;
           ao->setPath(ao->path().substr(4));
           if ( ao->path() == "/_XSEC" )
             xsec = dynamic_pointer_cast<Scatter1D>(ao);
           else if ( ao->path() == "/_EVTCOUNT" )
             sow = dynamic_pointer_cast<Counter>(ao);
           else {
             stripOptions(ao, delopts);
             string ananame = split(ao->path(), "/")[0];
             if ( ananames.insert(ananame).second ) addAnalysis(ananame);
             aos.push_back(ao);
           }
         }
         if ( !xsec || !sow ) {
           MSG_ERROR( "Error in AnalysisHandler::mergeYodas: The file " << file
                      << " did not contain weights and cross section info.");
           exit(1);
         }
         xsecs.push_back(xsec->point(0).x());
         sows.push_back(sow);
-	xsecerrs.push_back(sqr(xsec->point(0).xErrAvg()));
-        _eventcounter += *sow;
+        xsecerrs.push_back(sqr(xsec->point(0).xErrAvg()));
+        _eventCounter->operator+=(*sow); //< HAHAHAHAHA!
         sows.push_back(sow);
         aosv.push_back(aos);
       } catch (...) { //< YODA::ReadError&
         throw UserError("Unexpected error in reading file: " + file);
       }
     }
 
     // Now calculate the scale to be applied for all bins in a file
     // and get the common cross section and sum of weights.
     _xs = _xserr = 0.0;
     for ( int i = 0, N = sows.size(); i < N; ++i ) {
       double effnent = sows[i]->effNumEntries();
       _xs += (equiv? effnent: 1.0)*xsecs[i];
       _xserr += (equiv? sqr(effnent): 1.0)*xsecerrs[i];
     }
 
     vector<double> scales(sows.size(), 1.0);
     if ( equiv ) {
-      _xs /= _eventcounter.effNumEntries();
-      _xserr = sqrt(_xserr)/_eventcounter.effNumEntries();
+      _xs /= _eventCounter.effNumEntries();
+      _xserr = sqrt(_xserr)/_eventCounter.effNumEntries();
     } else {
       _xserr = sqrt(_xserr);
       for ( int i = 0, N = sows.size(); i < N; ++i )
-        scales[i] = (_eventcounter.sumW()/sows[i]->sumW())*(xsecs[i]/_xs);
+        scales[i] = (_eventCounter.sumW()/sows[i]->sumW())*(xsecs[i]/_xs);
     }
 
     // Initialize the analyses allowing them to book analysis objects.
     for (AnaHandle a : _analyses) {
       MSG_DEBUG("Initialising analysis: " << a->name());
       if ( !a->info().reentrant() )
         MSG_WARNING("Analysis " << a->name() << " has not been validated to have "
                     << "a reentrant finalize method. The result is unpredictable.");
       try {
         // Allow projection registration in the init phase onwards
         a->_allowProjReg = true;
         cerr << "sqrtS " << sqrtS() << endl;
         a->init();
         //MSG_DEBUG("Checking consistency of analysis: " << a->name());
         //a->checkConsistency();
       } catch (const Error& err) {
         cerr << "Error in " << a->name() << "::init method: " << err.what() << endl;
         exit(1);
       }
       MSG_DEBUG("Done initialising analysis: " << a->name());
     }
     _initialised = true;
     // Get a list of all anaysis objects to handle.
     map<string,AnalysisObjectPtr> current;
     for ( auto ao : getData(false, true) ) current[ao->path()] = ao;
     // Go through all objects to be merged and add them to current
     // after appropriate scaling.
     for ( int i = 0, N = aosv.size(); i < N; ++i)
       for ( auto ao : aosv[i] ) {
         if ( ao->path() == "/_XSEC" || ao->path() == "_EVTCOUNT" ) continue;
-	auto aoit = current.find(ao->path());
+        auto aoit = current.find(ao->path());
         if ( aoit == current.end() ) {
           MSG_WARNING("" << ao->path() << " was not properly booked.");
           continue;
         }
         if ( !addaos(aoit->second, ao, scales[i]) )
           MSG_WARNING("Cannot merge objects with path " << ao->path()
                       <<" of type " << ao->annotation("Type") );
       }
     // Now we can simply finalize() the analysis, leaving the
     // controlling program to write it out some yoda-file.
     finalize();
-
   }
 
 
   void AnalysisHandler::readData(const string& filename) {
     vector<YODA::AnalysisObjectPtr> aos;
     try {
       /// @todo Use new YODA SFINAE to fill the smart ptr vector directly
       vector<YODA::AnalysisObject*> aos_raw;
       YODA::ReaderYODA::read(filename, aos_raw);
       for (YODA::AnalysisObject* aor : aos_raw) aos.push_back(YODA::AnalysisObjectPtr(aor));
     //} catch (const YODA::ReadError & e) {
     } catch (...) { //< YODA::ReadError&
       throw UserError("Unexpected error in reading file: " + filename);
     }
     if (!aos.empty()) addData(aos);
   }
 
 
   vector<MultiweightAOPtr> AnalysisHandler::getRivetAOs() const {
       vector<MultiweightAOPtr> rtn;
 
       for (AnaHandle a : _analyses) {
           for (const auto & ao : a->analysisObjects()) {
               rtn.push_back(ao);
           }
       }
 
       rtn.push_back(_eventCounter);
       rtn.push_back(_xs);
 
       return rtn;
   }
 
   vector<YODA::AnalysisObjectPtr> AnalysisHandler::getYodaAOs(bool includeorphans,
                                                               bool includetmps,
                                                               bool usefinalized) const {
       vector<YODA::AnalysisObjectPtr> rtn;
       if (usefinalized)
         rtn = _finalizedAOs;
       else {
         for (auto rao : getRivetAOs()) {
           // need to set the index
           // before we can search the PATH
           rao.get()->setActiveWeightIdx(_defaultWeightIdx);
           // Exclude paths from final write-out if they contain a "TMP" layer (i.e. matching "/TMP/")
           if (!includetmps && rao->path().find("/TMP/") != string::npos)
             continue;
 
           for (size_t iW = 0; iW < numWeights(); iW++) {
             rao.get()->setActiveWeightIdx(iW);
             rtn.push_back(rao.get()->activeYODAPtr());
           }
         }
       }
 
       // Sort histograms alphanumerically by path before write-out
       sort(rtn.begin(), rtn.end(),
            [](YODA::AnalysisObjectPtr a, YODA::AnalysisObjectPtr b) {
                 return a->path() < b->path();
             }
           );
 
       return rtn;
   }
 
 
   vector<YODA::AnalysisObjectPtr> AnalysisHandler::getData(bool includeorphans,
                                                            bool includetmps,
                                                            bool usefinalized) const {
     return getYodaAOs(includeorphans, includetmps, usefinalized);
   }
 
 
   void AnalysisHandler::writeData(const string& filename) const {
     vector<YODA::AnalysisObjectPtr> out = _finalizedAOs;
     set<string> finalana;
     for ( auto ao : out) finalana.insert(ao->path());
     out.reserve(2*out.size());
     vector<YODA::AnalysisObjectPtr> aos = getData(false, true, false);
 
     if ( _dumping ) {
       for ( auto ao : aos ) {
         if ( finalana.find(ao->path()) == finalana.end() )
           out.push_back(AnalysisObjectPtr(ao->newclone()));
       }
     }
 
     for ( auto ao : aos ) {
       ao = YODA::AnalysisObjectPtr(ao->newclone());
       ao->setPath("/RAW" + ao->path());
       out.push_back(ao);
     }
 
     try {
       YODA::WriterYODA::write(filename, aos.begin(), aos.end());
     } catch (...) { //< YODA::WriteError&
       throw UserError("Unexpected error in writing file: " + filename);
     }
   }
 
 
   string AnalysisHandler::runName() const { return _runname; }
   size_t AnalysisHandler::numEvents() const { return _eventCounter->numEntries(); }
 
 
   std::vector<std::string> AnalysisHandler::analysisNames() const {
     std::vector<std::string> rtn;
     for (AnaHandle a : _analyses) {
       rtn.push_back(a->name());
     }
     return rtn;
   }
 
 
   const AnaHandle AnalysisHandler::analysis(const std::string& analysisname) const {
     for (const AnaHandle a : analyses())
       if (a->name() == analysisname) return a;
     throw Error("No analysis named '" + analysisname + "' registered in AnalysisHandler");
   }
 
 
   AnalysisHandler& AnalysisHandler::addAnalyses(const std::vector<std::string>& analysisnames) {
     for (const string& aname : analysisnames) {
       //MSG_DEBUG("Adding analysis '" << aname << "'");
       addAnalysis(aname);
     }
     return *this;
   }
 
 
   AnalysisHandler& AnalysisHandler::removeAnalyses(const std::vector<std::string>& analysisnames) {
     for (const string& aname : analysisnames) {
       removeAnalysis(aname);
     }
     return *this;
   }
 
 
   AnalysisHandler& AnalysisHandler::setCrossSection(double xs, double xserr) {
     _xs = Scatter1DPtr(weightNames(), Scatter1D("_XSEC"));
     _eventCounter.get()->setActiveWeightIdx(_defaultWeightIdx);
     double nomwgt = sumOfWeights();
 
     // The cross section of each weight variation is the nominal cross section
     // times the sumOfWeights(variation) / sumOfWeights(nominal).
     // This way the cross section will work correctly
     for (size_t iW = 0; iW < numWeights(); iW++) {
       _eventCounter.get()->setActiveWeightIdx(iW);
       double s = sumOfWeights() / nomwgt;
       _xs.get()->setActiveWeightIdx(iW);
       _xs->addPoint(xs*s, xserr*s);
     }
 
     _eventCounter.get()->unsetActiveWeight();
     _xs.get()->unsetActiveWeight();
     return *this;
   }
 
 
   AnalysisHandler& AnalysisHandler::addAnalysis(Analysis* analysis) {
     analysis->_analysishandler = this;
     _analyses.insert(AnaHandle(analysis));
     return *this;
   }
 
 
   PdgIdPair AnalysisHandler::beamIds() const {
     return Rivet::beamIds(beams());
   }
 
 
   double AnalysisHandler::sqrtS() const {
     return Rivet::sqrtS(beams());
   }
 
   void AnalysisHandler::setIgnoreBeams(bool ignore) {
     _ignoreBeams=ignore;
   }
 
 
 }