diff --git a/include/Rivet/Analysis.hh b/include/Rivet/Analysis.hh
--- a/include/Rivet/Analysis.hh
+++ b/include/Rivet/Analysis.hh
@@ -1,1226 +1,1226 @@
 // -*- 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"
 
 
 /// @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 {
 
 
   // 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;
     }
 
 
     /// Return true if this analysis needs to know the process cross-section.
     /// @todo Remove this and require HepMC >= 2.06
     bool needsCrossSection() const {
       return info().needsCrossSection();
     }
     /// Declare whether this analysis needs to know the process cross-section from the generator.
     /// @todo Remove this and require HepMC >= 2.06
     Analysis& setNeedsCrossSection(bool needed=true) {
       info().setNeedsCrossSection(needed);
       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;
 
     //@}
 
 
     /// @name Analysis / beam compatibility testing
     //@{
 
     /// 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;
 
     //@}
 
 
     /// Set the cross section from the generator
     Analysis& setCrossSection(double xs); //, double xserr=0.0);
 
     /// 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;
     /// Alias
     /// @deprecated Prefer the "mk" form, consistent with other "making function" names
     const std::string makeAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
       return mkAxisCode(datasetId, xAxisId, yAxisId);
     }
 
     //@}
 
 
     /// @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 = makeAxisCode(datasetId, xAxisId, yAxisId);
       return refData(hname);
     }
 
     //@}
 
 
     /// @name Counter booking
     //@{
 
     /// Book a counter.
     CounterPtr bookCounter(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 bookCounter(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 bookHisto1D(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 bookHisto1D(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 bookHisto1D(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 bookHisto1D(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 bookHisto1D(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 bookHisto1D(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 bookHisto2D(const std::string& name,
                            size_t nxbins, double xlower, double xupper,
                            size_t nybins, double ylower, double yupper,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="",
                            const std::string& ztitle="");
 
     /// Book a 2D histogram with non-uniform bins defined by the
     /// vectors of bin edges @a xbinedges and @a ybinedges.
     Histo2DPtr bookHisto2D(const std::string& name,
                            const std::vector<double>& xbinedges,
                            const std::vector<double>& ybinedges,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="",
                            const std::string& ztitle="");
 
     /// Book a 2D histogram with non-uniform bins defined by the
     /// vectors of bin edges @a xbinedges and @a ybinedges.
     Histo2DPtr bookHisto2D(const std::string& name,
                            const std::initializer_list<double>& xbinedges,
                            const std::initializer_list<double>& ybinedges,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="",
                            const std::string& ztitle="");
 
     /// Book a 2D histogram with binning from a reference scatter.
     Histo2DPtr bookHisto2D(const std::string& name,
                            const Scatter3D& refscatter,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="",
                            const std::string& ztitle="");
 
     /// Book a 2D histogram, using the binnings in the reference data histogram.
     Histo2DPtr bookHisto2D(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 bookHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="",
                            const std::string& ztitle="");
 
     //@}
 
 
     /// @name 1D profile histogram booking
     //@{
 
     /// Book a 1D profile histogram with @a nbins uniformly distributed across the range @a lower - @a upper .
     Profile1DPtr bookProfile1D(const std::string& name,
                                size_t nbins, double lower, double upper,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="");
 
     /// Book a 1D profile histogram with non-uniform bins defined by the vector of bin edges @a binedges .
     Profile1DPtr bookProfile1D(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 bookProfile1D(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 bookProfile1D(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 bookProfile1D(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 bookProfile1D(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 bookProfile2D(const std::string& name,
                                size_t nxbins, double xlower, double xupper,
                                size_t nybins, double ylower, double yupper,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="",
                                const std::string& ztitle="");
 
     /// Book a 2D profile histogram with non-uniform bins defined by the vectorx
     /// of bin edges @a xbinedges and @a ybinedges.
     Profile2DPtr bookProfile2D(const std::string& name,
                                const std::vector<double>& xbinedges,
                                const std::vector<double>& ybinedges,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="",
                                const std::string& ztitle="");
 
     /// Book a 2D profile histogram with non-uniform bins defined by the vectorx
     /// of bin edges @a xbinedges and @a ybinedges.
     Profile2DPtr bookProfile2D(const std::string& name,
                                const std::initializer_list<double>& xbinedges,
                                const std::initializer_list<double>& ybinedges,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="",
                                const std::string& ztitle="");
 
     /// Book a 2D profile histogram with binning from a reference scatter.
     Profile2DPtr bookProfile2D(const std::string& name,
                                const Scatter3D& refscatter,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="",
                                const std::string& ztitle="");
 
     /// Book a 2D profile histogram, using the binnings in the reference data histogram.
     Profile2DPtr bookProfile2D(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 bookProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="",
                                const std::string& ztitle="");
 
     //@}
 
 
     /// @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 bookScatter2D(const std::string& name,
                                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 bookScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                                bool copy_pts=false,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="");
 
     /// @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 bookScatter2D(const std::string& name,
                                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 bookScatter2D(const std::string& hname,
                                const std::vector<double>& binedges,
                                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 Pecentile 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.
     template <class T>
     Percentile<T> bookPercentile(string projName,
                                  vector<pair<float, float> > centralityBins,
                                  vector<tuple<int, int, int> > ref) {
       typedef typename ReferenceTraits<T>::RefT RefT;
       Percentile<T> pctl(this, projName);
 
       const int nCent = centralityBins.size();
       for (int iCent = 0; iCent < nCent; ++iCent) {
         const string axisCode = makeAxisCode(std::get<0>(ref[iCent]),
                                              std::get<1>(ref[iCent]),
                                              std::get<2>(ref[iCent])); 
         shared_ptr<T> ao;
         CounterPtr cnt;
         try {
           ao = getAnalysisObject<T>(axisCode);
           MSG_TRACE("Found old " << histoPath(axisCode));
         }
         catch (Exception) {
           const RefT & refscatter = refData<RefT>(axisCode);
           ao = make_shared<T>(refscatter, histoPath(axisCode));
           addAnalysisObject(ao);
           MSG_TRACE("Created new " << histoPath(axisCode));
         }
         try {
           cnt = getAnalysisObject<Counter>("TMP/COUNTER/" + axisCode);
           MSG_TRACE("Found old " << histoPath("TMP/COUNTER/" + axisCode));
         }
         catch (Exception) {
           cnt = make_shared<Counter>(histoPath("TMP/COUNTER/" + axisCode));
           addAnalysisObject(cnt);
           MSG_TRACE("Created new " << histoPath("TMP/COUNTER/" + axisCode));
         }
         pctl.add(ao, cnt, centralityBins[iCent]);
       }
       return pctl;
     }
 
     /// @brief Book Pecentile 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.
     template <class T>
     PercentileXaxis<T> bookPercentileXaxis(string projName,
                                            tuple<int, int, int> ref) {
       typedef typename ReferenceTraits<T>::RefT RefT;
       PercentileXaxis<T> pctl(this, projName);
 
       const string axisCode = makeAxisCode(std::get<0>(ref),
                                            std::get<1>(ref),
                                            std::get<2>(ref)); 
       shared_ptr<T> ao;
       CounterPtr cnt;
       try {
         ao = getAnalysisObject<T>(histoPath(axisCode));
       }
       catch (Exception) {
 	const RefT & refscatter = refData<RefT>(axisCode);
         ao = make_shared<T>(refscatter, axisCode);
         addAnalysisObject(ao);
       }
       pctl.add(proj, ao, make_shared<Counter>());
 
       return pctl;
     }
 
 
 
     /// @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, double 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, double factor) {
       for (auto& c : cnts) scale(c, factor);
     }
     /// @todo YUCK!
     template <std::size_t array_size>
     void scale(const CounterPtr (&cnts)[array_size], double factor) {
       // for (size_t i = 0; i < std::extent<decltype(cnts)>::value; ++i) scale(cnts[i], factor);
       for (auto& c : cnts) scale(c, factor);
     }
 
 
     /// Normalize the given histogram, @a histo, to area = @a norm.
     void normalize(Histo1DPtr histo, double 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, double 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], double 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, double 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, double factor) {
       for (auto& h : histos) scale(h, factor);
     }
     /// @todo YUCK!
     template <std::size_t array_size>
     void scale(const Histo1DPtr (&histos)[array_size], double factor) {
       for (auto& h : histos) scale(h, factor);
     }
 
 
     /// Normalize the given histogram, @a histo, to area = @a norm.
     void normalize(Histo2DPtr histo, double 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, double 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], double 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, double 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, double factor) {
       for (auto& h : histos) scale(h, factor);
     }
     /// @todo YUCK!
     template <std::size_t array_size>
     void scale(const Histo2DPtr (&histos)[array_size], double 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<AnalysisObjectPtr>& analysisObjects() const {
       return _analysisobjects;
     }
 
 
   protected:
 
     /// @name Data object registration, retrieval, and removal
     //@{
 
     /// Register a data object in the histogram system
     void addAnalysisObject(AnalysisObjectPtr ao);
 
     /// 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 Exception("Data object " + histoPath(name) + " not found");
+      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 Exception("Data object " + histoPath(name) + " not found");
+      throw LookupError("Data object " + histoPath(name) + " not found");
     }
 
     /// 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(AnalysisObjectPtr ao);
 
     /// Get all data object from the AnalysisHandler.
     vector<AnalysisObjectPtr> getAllData(bool includeorphans) const;
 
     /// Get a data object from another analysis (e.g. preloaded
     /// calibration histogram).
     /// Get a data object from the histogram system (non-const)
     template <typename AO=YODA::AnalysisObject>
     std::shared_ptr<AO> getAnalysisObject(const std::string & ananame,
                                           const std::string& name) {
       std::string path = "/" + ananame + "/" + name;
       for ( AnalysisObjectPtr ao : getAllData(true) ) {
         if ( ao->path() == path )
           return dynamic_pointer_cast<AO>(ao);
       }
       return std::shared_ptr<AO>();
     }
 
     /// 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<AnalysisObjectPtr> _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, 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/Tools/Exceptions.hh b/include/Rivet/Tools/Exceptions.hh
--- a/include/Rivet/Tools/Exceptions.hh
+++ b/include/Rivet/Tools/Exceptions.hh
@@ -1,62 +1,67 @@
 #ifndef RIVET_EXCEPTIONS_HH
 #define RIVET_EXCEPTIONS_HH
 
 #include <string>
 #include <exception>
 #include <stdexcept>
 
 namespace Rivet {
 
 
   /// @brief Generic runtime Rivet error.
   struct Error : public std::runtime_error {
     Error(const std::string& what) : std::runtime_error(what) {}
   };
 
 
   /// @brief Rivet::Exception is a synonym for Rivet::Error.
   typedef Error Exception;
 
 
   /// @brief Error for e.g. use of invalid bin ranges.
   struct RangeError : public Error {
     RangeError(const std::string& what) : Error(what) {}
   };
 
 
   /// @brief Error specialisation for places where alg logic has failed.
   struct LogicError : public Error {
     LogicError(const std::string& what) : Error(what) {}
   };
 
 
   /// @brief Error specialisation for failures relating to particle ID codes.
   struct PidError : public Error {
     PidError(const std::string& what) : Error(what) {}
   };
 
 
   /// @brief Error specialisation for failures relating to analysis info.
   struct InfoError : public Error {
     InfoError(const std::string& what) : Error(what) {}
   };
 
 
   /// @brief Errors relating to event/bin weights
   ///
   /// Arises in computing statistical quantities because e.g. the bin
   /// weight is zero or negative.
   struct WeightError : public Error {
     WeightError(const std::string& what) : Error(what) {}
   };
 
 
   /// @brief Error specialisation for where the problem is between the chair and the computer.
   struct UserError : public Error {
     UserError(const std::string& what) : Error(what) {}
   };
 
+  /// @brief Error relating to looking up analyis objects in the register
+  struct LookupError : public Error {
+    LookupError(const std::string& what) : Error(what) {}
+  };
+
 
 }
 
 #endif
diff --git a/src/Core/Analysis.cc b/src/Core/Analysis.cc
--- a/src/Core/Analysis.cc
+++ b/src/Core/Analysis.cc
@@ -1,1004 +1,1009 @@
 // -*- 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)
     : _crossSection(-1.0),
       _gotCrossSection(false),
       _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 {
     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::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;
 
     /// @todo Need to also check internal consistency of the analysis'
     /// beam requirements with those of the projections it uses.
   }
 
 
   ///////////////////////////////////////////
 
 
   Analysis& Analysis::setCrossSection(double xs) {
     _crossSection = xs;
     _gotCrossSection = true;
     return *this;
   }
 
   double Analysis::crossSection() const {
     if (!_gotCrossSection || std::isnan(_crossSection)) {
       string errMsg = "You did not set the cross section for the analysis " + name();
       throw Error(errMsg);
     }
     return _crossSection;
   }
 
   double Analysis::crossSectionPerEvent() const {
     const double sumW = sumOfWeights();
     assert(sumW != 0.0);
     return _crossSection / sumW;
   }
 
 
 
   ////////////////////////////////////////////////////////////
   // Histogramming
 
 
   void Analysis::_cacheRefData() const {
     if (_refdata.empty()) {
       MSG_TRACE("Getting refdata cache for paper " << name());
       _refdata = getRefData(getRefDataName());
     }
   }
 
   vector<AnalysisObjectPtr> Analysis::getAllData(bool includeorphans) const{
     return handler().getData(includeorphans);
   }
 
   CounterPtr Analysis::bookCounter(const string& cname,
                                    const string& title) {
                                    // const string& xtitle,
                                    // const string& ytitle) {
     const string path = histoPath(cname);
     CounterPtr ctr = make_shared<Counter>(path, title);
     addAnalysisObject(ctr);
     MSG_TRACE("Made counter " << cname << " for " << name());
     // hist->setAnnotation("XLabel", xtitle);
     // hist->setAnnotation("YLabel", ytitle);
     return ctr;
   }
 
 
   CounterPtr Analysis::bookCounter(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 bookCounter(axisCode, title);
   }
 
 
   Histo1DPtr Analysis::bookHisto1D(const string& hname,
                                    size_t nbins, double lower, double upper,
                                    const string& title,
                                    const string& xtitle,
                                    const string& ytitle) {
-    Histo1DPtr hist;
+    Histo1DPtr hist = make_shared<Histo1D>(nbins, lower, upper, histoPath(hname), title);
     try { // try to bind to pre-existing
+      if ( !getHisto1D(hname)->sameBinning(*hist) ) {
+        throw Exception("Histogram " + hname + " already exists with a different binning");
+      }
       // AnalysisObjectPtr ao = getAnalysisObject(path);
       // hist = dynamic_pointer_cast<Histo1D>(ao);
       hist = getHisto1D(hname);
       /// @todo Test that cast worked
       /// @todo Also test that binning is as expected?
       MSG_TRACE("Bound pre-existing histogram " << hname <<  " for " << name());
-    } catch (...) { // binding failed; make it from scratch
-      hist = make_shared<Histo1D>(nbins, lower, upper, histoPath(hname), title);
+    } catch (LookupError) { // binding failed; make it from scratch 
       addAnalysisObject(hist);
       MSG_TRACE("Made histogram " << hname <<  " for " << name());
     }
     hist->setTitle(title);
     hist->setAnnotation("XLabel", xtitle);
     hist->setAnnotation("YLabel", ytitle);
     return hist;
   }
 
 
   Histo1DPtr Analysis::bookHisto1D(const string& hname,
                                    const vector<double>& binedges,
                                    const string& title,
                                    const string& xtitle,
                                    const string& ytitle) {
-    Histo1DPtr hist;
+    Histo1DPtr hist = make_shared<Histo1D>(binedges, histoPath(hname), title);
     try { // try to bind to pre-existing
+      if ( !getHisto1D(hname)->sameBinning(*hist) ) {
+        throw Exception("Histogram " + hname + " already exists with a different binning");
+      }
       // AnalysisObjectPtr ao = getAnalysisObject(path);
       // hist = dynamic_pointer_cast<Histo1D>(ao);
       hist = getHisto1D(hname);
       /// @todo Test that cast worked
       /// @todo Also test that binning is as expected?
       MSG_TRACE("Bound pre-existing histogram " << hname <<  " for " << name());
-    } catch (...) { // binding failed; make it from scratch
-      hist = make_shared<Histo1D>(binedges, histoPath(hname), title);
+    } catch (LookupError) { // binding failed; make it from scratch 
       addAnalysisObject(hist);
       MSG_TRACE("Made histogram " << hname <<  " for " << name());
     }
     hist->setTitle(title);
     hist->setAnnotation("XLabel", xtitle);
     hist->setAnnotation("YLabel", ytitle);
     return hist;
   }
 
 
   Histo1DPtr Analysis::bookHisto1D(const string& hname,
                                    const initializer_list<double>& binedges,
                                    const string& title,
                                    const string& xtitle,
                                    const string& ytitle) {
     return bookHisto1D(hname, vector<double>{binedges}, title, xtitle, ytitle);
   }
 
 
   Histo1DPtr Analysis::bookHisto1D(const string& hname,
                                    const Scatter2D& refscatter,
                                    const string& title,
                                    const string& xtitle,
                                    const string& ytitle) {
-    Histo1DPtr hist;
+    Histo1DPtr hist = make_shared<Histo1D>(refscatter, histoPath(hname));
     try { // try to bind to pre-existing
+      if ( !getHisto1D(hname)->sameBinning(*hist) ) {
+        throw Exception("Histogram " + hname + " already exists with a different binning");
+      }
       // AnalysisObjectPtr ao = getAnalysisObject(path);
       // hist = dynamic_pointer_cast<Histo1D>(ao);
       hist = getHisto1D(hname);
       /// @todo Test that cast worked
       /// @todo Also test that binning is as expected?
       MSG_TRACE("Bound pre-existing histogram " << hname <<  " for " << name());
-    } catch (...) { // binding failed; make it from scratch
-      hist = make_shared<Histo1D>(refscatter, histoPath(hname));
-      if (hist->hasAnnotation("IsRef")) hist->rmAnnotation("IsRef");
+    } catch (LookupError) { // binding failed; make it from scratch 
       addAnalysisObject(hist);
       MSG_TRACE("Made histogram " << hname <<  " for " << name());
     }
     hist->setTitle(title);
     hist->setAnnotation("XLabel", xtitle);
     hist->setAnnotation("YLabel", ytitle);
     return hist;
   }
 
 
   Histo1DPtr Analysis::bookHisto1D(const string& hname,
                                    const string& title,
                                    const string& xtitle,
                                    const string& ytitle) {
     const Scatter2D& refdata = refData(hname);
     return bookHisto1D(hname, refdata, title, xtitle, ytitle);
   }
 
 
   Histo1DPtr Analysis::bookHisto1D(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 bookHisto1D(axisCode, title, xtitle, ytitle);
   }
 
 
   /// @todo Add booking methods which take a path, titles and *a reference Scatter from which to book*
 
 
   /////////////////
 
 
   Histo2DPtr Analysis::bookHisto2D(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);
     Histo2DPtr hist = make_shared<Histo2D>(nxbins, xlower, xupper, nybins, ylower, yupper, path, title);
     addAnalysisObject(hist);
     MSG_TRACE("Made 2D histogram " << hname <<  " for " << name());
     hist->setAnnotation("XLabel", xtitle);
     hist->setAnnotation("YLabel", ytitle);
     hist->setAnnotation("ZLabel", ztitle);
     return hist;
   }
 
 
   Histo2DPtr Analysis::bookHisto2D(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);
     Histo2DPtr hist = make_shared<Histo2D>(xbinedges, ybinedges, path, title);
     addAnalysisObject(hist);
     MSG_TRACE("Made 2D histogram " << hname <<  " for " << name());
     hist->setAnnotation("XLabel", xtitle);
     hist->setAnnotation("YLabel", ytitle);
     hist->setAnnotation("ZLabel", ztitle);
     return hist;
   }
 
 
   Histo2DPtr Analysis::bookHisto2D(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 bookHisto2D(hname, vector<double>{xbinedges}, vector<double>{ybinedges},
                        title, xtitle, ytitle, ztitle);
   }
 
 
   Histo2DPtr Analysis::bookHisto2D(const string& hname,
                                    const Scatter3D& refscatter,
                                    const string& title,
                                    const string& xtitle,
                                    const string& ytitle,
                                    const string& ztitle) {
     const string path = histoPath(hname);
     Histo2DPtr hist( new Histo2D(refscatter, path) );
     addAnalysisObject(hist);
     MSG_TRACE("Made 2D histogram " << hname <<  " for " << name());
     if (hist->hasAnnotation("IsRef")) hist->rmAnnotation("IsRef");
     hist->setTitle(title);
     hist->setAnnotation("XLabel", xtitle);
     hist->setAnnotation("YLabel", ytitle);
     hist->setAnnotation("ZLabel", ztitle);
     return hist;
   }
 
 
   Histo2DPtr Analysis::bookHisto2D(const string& hname,
                                    const string& title,
                                    const string& xtitle,
                                    const string& ytitle,
                                    const string& ztitle) {
     const Scatter3D& refdata = refData<Scatter3D>(hname);
     return bookHisto2D(hname, refdata, title, xtitle, ytitle, ztitle);
   }
 
 
   Histo2DPtr Analysis::bookHisto2D(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 bookHisto2D(axisCode, title, xtitle, ytitle, ztitle);
   }
 
 
   /////////////////
 
 
   Profile1DPtr Analysis::bookProfile1D(const string& hname,
                                        size_t nbins, double lower, double upper,
                                        const string& title,
                                        const string& xtitle,
                                        const string& ytitle) {
     const string path = histoPath(hname);
     Profile1DPtr prof = make_shared<Profile1D>(nbins, lower, upper, path, title);
     addAnalysisObject(prof);
     MSG_TRACE("Made profile histogram " << hname <<  " for " << name());
     prof->setAnnotation("XLabel", xtitle);
     prof->setAnnotation("YLabel", ytitle);
     return prof;
   }
 
 
   Profile1DPtr Analysis::bookProfile1D(const string& hname,
                                        const vector<double>& binedges,
                                        const string& title,
                                        const string& xtitle,
                                        const string& ytitle) {
     const string path = histoPath(hname);
     Profile1DPtr prof = make_shared<Profile1D>(binedges, path, title);
     addAnalysisObject(prof);
     MSG_TRACE("Made profile histogram " << hname <<  " for " << name());
     prof->setAnnotation("XLabel", xtitle);
     prof->setAnnotation("YLabel", ytitle);
     return prof;
   }
 
 
   Profile1DPtr Analysis::bookProfile1D(const string& hname,
                                        const initializer_list<double>& binedges,
                                        const string& title,
                                        const string& xtitle,
                                        const string& ytitle)
   {
     return bookProfile1D(hname, vector<double>{binedges}, title, xtitle, ytitle);
   }
 
 
   Profile1DPtr Analysis::bookProfile1D(const string& hname,
                                        const Scatter2D& refscatter,
                                        const string& title,
                                        const string& xtitle,
                                        const string& ytitle) {
     const string path = histoPath(hname);
     Profile1DPtr prof = make_shared<Profile1D>(refscatter, path);
     addAnalysisObject(prof);
     MSG_TRACE("Made profile histogram " << hname <<  " for " << name());
     if (prof->hasAnnotation("IsRef")) prof->rmAnnotation("IsRef");
     prof->setTitle(title);
     prof->setAnnotation("XLabel", xtitle);
     prof->setAnnotation("YLabel", ytitle);
     return prof;
   }
 
 
   Profile1DPtr Analysis::bookProfile1D(const string& hname,
                                        const string& title,
                                        const string& xtitle,
                                        const string& ytitle) {
     const Scatter2D& refdata = refData(hname);
     return bookProfile1D(hname, refdata, title, xtitle, ytitle);
   }
 
 
   Profile1DPtr Analysis::bookProfile1D(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 bookProfile1D(axisCode, title, xtitle, ytitle);
   }
 
 
   ///////////////////
 
 
 
   Profile2DPtr Analysis::bookProfile2D(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);
     Profile2DPtr prof = make_shared<Profile2D>(nxbins, xlower, xupper, nybins, ylower, yupper, path, title);
     addAnalysisObject(prof);
     MSG_TRACE("Made 2D profile histogram " << hname <<  " for " << name());
     prof->setAnnotation("XLabel", xtitle);
     prof->setAnnotation("YLabel", ytitle);
     prof->setAnnotation("ZLabel", ztitle);
     return prof;
   }
 
 
   Profile2DPtr Analysis::bookProfile2D(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);
     Profile2DPtr prof = make_shared<Profile2D>(xbinedges, ybinedges, path, title);
     addAnalysisObject(prof);
     MSG_TRACE("Made 2D profile histogram " << hname <<  " for " << name());
     prof->setAnnotation("XLabel", xtitle);
     prof->setAnnotation("YLabel", ytitle);
     prof->setAnnotation("ZLabel", ztitle);
     return prof;
   }
 
 
   Profile2DPtr Analysis::bookProfile2D(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 bookProfile2D(hname, vector<double>{xbinedges}, vector<double>{ybinedges},
                          title, xtitle, ytitle, ztitle);
   }
 
 
   Profile2DPtr Analysis::bookProfile2D(const string& hname,
                                        const Scatter3D& refscatter,
                                        const string& title,
                                        const string& xtitle,
                                        const string& ytitle,
                                        const string& ztitle) {
     const string path = histoPath(hname);
     Profile2DPtr prof( new Profile2D(refscatter, path) );
     addAnalysisObject(prof);
     MSG_TRACE("Made 2D profile histogram " << hname <<  " for " << name());
     if (prof->hasAnnotation("IsRef")) prof->rmAnnotation("IsRef");
     prof->setTitle(title);
     prof->setAnnotation("XLabel", xtitle);
     prof->setAnnotation("YLabel", ytitle);
     prof->setAnnotation("ZLabel", ztitle);
     return prof;
   }
 
 
   Profile2DPtr Analysis::bookProfile2D(const string& hname,
                                        const string& title,
                                        const string& xtitle,
                                        const string& ytitle,
                                        const string& ztitle) {
     const Scatter3D& refdata = refData<Scatter3D>(hname);
     return bookProfile2D(hname, refdata, title, xtitle, ytitle, ztitle);
   }
 
 
   Profile2DPtr Analysis::bookProfile2D(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 bookProfile2D(axisCode, title, xtitle, ytitle, ztitle);
   }
 
 
   /////////////////
 
 
   Scatter2DPtr Analysis::bookScatter2D(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 bookScatter2D(axisCode, copy_pts, title, xtitle, ytitle);
   }
 
 
   Scatter2DPtr Analysis::bookScatter2D(const string& hname,
                                        bool copy_pts,
                                        const string& title,
                                        const string& xtitle,
                                        const string& ytitle) {
     Scatter2DPtr s;
     const string path = histoPath(hname);
     if (copy_pts) {
       const Scatter2D& refdata = refData(hname);
       s = make_shared<Scatter2D>(refdata, path);
       for (Point2D& p : s->points()) p.setY(0, 0);
     } else {
       s = make_shared<Scatter2D>(path);
     }
     addAnalysisObject(s);
     MSG_TRACE("Made scatter " << hname <<  " for " << name());
     if (s->hasAnnotation("IsRef")) s->rmAnnotation("IsRef");
     s->setTitle(title);
     s->setAnnotation("XLabel", xtitle);
     s->setAnnotation("YLabel", ytitle);
     return s;
   }
 
 
   Scatter2DPtr Analysis::bookScatter2D(const string& hname,
                                        size_t npts, double lower, double upper,
                                        const string& title,
                                        const string& xtitle,
                                        const string& ytitle) {
     Scatter2DPtr s;
     const string path = histoPath(hname);
     try { // try to bind to pre-existing
       s = getAnalysisObject<Scatter2D>(hname);
       /// @todo Also test that binning is as expected?
       MSG_TRACE("Bound pre-existing scatter " << path <<  " for " << name());
     } catch (...) { // binding failed; make it from scratch
       s = make_shared<Scatter2D>(path);
       const double binwidth = (upper-lower)/npts;
       for (size_t pt = 0; pt < npts; ++pt) {
         const double bincentre = lower + (pt + 0.5) * binwidth;
         s->addPoint(bincentre, 0, binwidth/2.0, 0);
       }
       addAnalysisObject(s);
       MSG_TRACE("Made scatter " << hname <<  " for " << name());
     }
     s->setTitle(title);
     s->setAnnotation("XLabel", xtitle);
     s->setAnnotation("YLabel", ytitle);
     return s;
   }
 
 
   Scatter2DPtr Analysis::bookScatter2D(const string& hname,
                                        const vector<double>& binedges,
                                        const string& title,
                                        const string& xtitle,
                                        const string& ytitle) {
     const string path = histoPath(hname);
     Scatter2DPtr s = make_shared<Scatter2D>(path);
     for (size_t pt = 0; pt < binedges.size()-1; ++pt) {
       const double bincentre = (binedges[pt] + binedges[pt+1]) / 2.0;
       const double binwidth = binedges[pt+1] - binedges[pt];
       s->addPoint(bincentre, 0, binwidth/2.0, 0);
     }
     addAnalysisObject(s);
     MSG_TRACE("Made scatter " << hname <<  " for " << name());
     s->setTitle(title);
     s->setAnnotation("XLabel", xtitle);
     s->setAnnotation("YLabel", ytitle);
     return s;
   }
 
 
   /////////////////////
 
 
   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, double factor) {
     if (!cnt) {
       MSG_WARNING("Failed to scale counter=NULL in analysis " << name() << " (scale=" << factor << ")");
       return;
     }
     if (std::isnan(factor) || std::isinf(factor)) {
       MSG_WARNING("Failed to scale counter=" << cnt->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")");
       factor = 0;
     }
     MSG_TRACE("Scaling counter " << cnt->path() << " by factor " << factor);
     try {
       cnt->scaleW(factor);
     } catch (YODA::Exception& we) {
       MSG_WARNING("Could not scale counter " << cnt->path());
       return;
     }
   }
 
 
   void Analysis::normalize(Histo1DPtr histo, double norm, bool includeoverflows) {
     if (!histo) {
       MSG_WARNING("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << norm << ")");
       return;
     }
     MSG_TRACE("Normalizing histo " << histo->path() << " to " << 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, double factor) {
     if (!histo) {
       MSG_WARNING("Failed to scale histo=NULL in analysis " << name() << " (scale=" << factor << ")");
       return;
     }
     if (std::isnan(factor) || std::isinf(factor)) {
       MSG_WARNING("Failed to scale histo=" << histo->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")");
       factor = 0;
     }
     MSG_TRACE("Scaling histo " << histo->path() << " by factor " << factor);
     try {
       histo->scaleW(factor);
     } catch (YODA::Exception& we) {
       MSG_WARNING("Could not scale histo " << histo->path());
       return;
     }
   }
 
 
   void Analysis::normalize(Histo2DPtr histo, double norm, bool includeoverflows) {
     if (!histo) {
       MSG_ERROR("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << norm << ")");
       return;
     }
     MSG_TRACE("Normalizing histo " << histo->path() << " to " << 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, double factor) {
     if (!histo) {
       MSG_ERROR("Failed to scale histo=NULL in analysis " << name() << " (scale=" << factor << ")");
       return;
     }
     if (std::isnan(factor) || std::isinf(factor)) {
       MSG_ERROR("Failed to scale histo=" << histo->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")");
       factor = 0;
     }
     MSG_TRACE("Scaling histo " << histo->path() << " by factor " << 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?!?
 
 
   //////////////////////////////////
 
 
   void Analysis::addAnalysisObject(AnalysisObjectPtr ao) {
     _analysisobjects.push_back(ao);
   }
 
 
   void Analysis::removeAnalysisObject(const string& path) {
     for (vector<AnalysisObjectPtr>::iterator it = _analysisobjects.begin();  it != _analysisobjects.end(); ++it) {
       if ((*it)->path() == path) {
         _analysisobjects.erase(it);
         break;
       }
     }
   }
 
 
   void Analysis::removeAnalysisObject(AnalysisObjectPtr ao) {
     for (vector<AnalysisObjectPtr>::iterator 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;
     auto refmap = getRefData(calAnaName);
     if ( refmap.find(calHistName) != refmap.end() )
       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);
     }
   }
   else if ( sel == "GEN" ) {
     Histo1DPtr genhist;
     string histpath = "/" + calAnaName + "/" + calHistName;
     for ( AnalysisObjectPtr ao : handler().getData(true) ) {
       if ( ao->path() == histpath )
         genhist = dynamic_pointer_cast<Histo1D>(ao);
     }
     if ( !genhist || genhist->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);
     }
   }
   else if ( sel == "IMP" ) {
     Histo1DPtr imphist =
       getAnalysisObject<Histo1D>(calAnaName, calHistName + "_IMP");
     if ( !imphist || imphist->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());
       cproj.add(PercentileProjection(ImpactParameterProjection(),
                                      imphist, true), sel);
     }
   }
   else if ( sel == "USR" ) {
 #if HEPMC_VERSION_CODE >= 3000000
     Histo1DPtr usrhist =
       getAnalysisObject<Histo1D>(calAnaName, calHistName + "_USR");
     if ( !usrhist || usrhist->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);
      }
 #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);
   
 }
 
 }