No OneTemporary
Actions

Size

138 KB

Subscribers

None

View Options

	diff --git a/bin/rivet-merge b/bin/rivet-merge
	--- a/bin/rivet-merge
	+++ b/bin/rivet-merge
	@@ -1,66 +1,70 @@
	#! /usr/bin/env python

	"""\
	Merge together yoda files produced by Rivet.

	Examples:

	%prog [options] <yodafile> [<yodafile2> ...]

	ENVIRONMENT:
	* RIVET_ANALYSIS_PATH: list of paths to be searched for plugin
	analysis libraries at runtime
	* RIVET_DATA_PATH: list of paths to be searched for data files
	"""

	import os, sys

	## Load the rivet module
	try:
	import rivet
	except:
	## If rivet loading failed, try to bootstrap the Python path!
	try:
	# TODO: Is this a good idea? Maybe just notify the user that their PYTHONPATH is wrong?
	import commands
	modname = sys.modules[__name__].__file__
	binpath = os.path.dirname(modname)
	rivetconfigpath = os.path.join(binpath, "rivet-config")
	rivetpypath = commands.getoutput(rivetconfigpath + " --pythonpath")
	sys.path.append(rivetpypath)
	import rivet
	except:
	sys.stderr.write("The rivet Python module could not be loaded: is your PYTHONPATH set correctly?\n")
	sys.exit(5)

	rivet.util.check_python_version()
	rivet.util.set_process_name("rivet-merge")

	import time, datetime, logging, signal

	## Parse command line options
	from optparse import OptionParser, OptionGroup
	parser = OptionParser(usage=__doc__, version="rivet v%s" % rivet.version())

	extragroup = OptionGroup(parser, "Run settings")
	extragroup.add_option("-o", "--output-file", dest="OUTPUTFILE",
	default="Rivet.yoda", help="specify the output histo file path (default = %default)")
	extragroup.add_option("-e", "--equiv", dest="EQUIV", action="store_true", default=False,
	help="assume that the yoda files are equivalent but statistically independent (default= assume that different files contains different processes)")
	+
	+extragroup.add_option("-O", "--merge-option", dest="MERGEOPTIONS", action="append",
	+ default=[], help="specify an analysis option name where different options should be merged into the default analysis.")
	+
	parser.add_option_group(extragroup)


	opts, args = parser.parse_args()


	############################
	## Actual analysis runs



	## Set up analysis handler
	ah = rivet.AnalysisHandler("Merging")
	-ah.mergeYodas(args, opts.EQUIV)
	+ah.mergeYodas(args, opts.MERGEOPTIONS, opts.EQUIV)
	ah.writeData(opts.OUTPUTFILE)


	diff --git a/include/Rivet/Analysis.hh b/include/Rivet/Analysis.hh
	--- a/include/Rivet/Analysis.hh
	+++ b/include/Rivet/Analysis.hh
	@@ -1,1225 +1,1232 @@
	// -- 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;

	+ /// Check if we are running rivet-merge.
	+ bool merging() const {
	+ return sqrtS() <= 0.0;
	+ }
	+
	//@}


	/// @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));
	- }
	+ const RefT & refscatter = refData<RefT>(axisCode);
	+ shared_ptr<T> ao = addOrGetCompatAO(make_shared<T>(refscatter, histoPath(axisCode)));
	+ CounterPtr cnt =
	+ addOrGetCompatAO(make_shared<Counter>(histoPath("TMP/COUNTER/" + axisCode)));
	pctl.add(ao, cnt, centralityBins[iCent]);
	}
	return pctl;
	}

	/// @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);
	- }
	+ const RefT & refscatter = refData<RefT>(axisCode);
	+ shared_ptr<T> ao = addOrGetCompatAO(make_shared<T>(refscatter, histoPath(axisCode)));
	pctl.add(proj, ao, make_shared<Counter>());

	return pctl;
	}



	/// @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);

	+ /// Register a data object in the system and return its pointer,
	+ /// or, if an object of the same path is already there, check if it
	+ /// is compatible (eg. same type and same binning) and return that
	+ /// object instead. Emits a warning if an incompatible object with
	+ /// the same name is found and replcaces that with the given data
	+ /// object.
	+ template <typename AO=YODA::AnalysisObject>
	+ std::shared_ptr<AO> addOrGetCompatAO(std::shared_ptr<AO> aonew) {
	+ foreach (const AnalysisObjectPtr& ao, analysisObjects()) {
	+ if ( ao->path() == aonew->path() ) {
	+ std::shared_ptr<AO> aoold = dynamic_pointer_cast<AO>(ao);
	+ if ( aoold && bookingCompatible(aonew, aoold) ) {
	+ MSG_TRACE("Bound pre-existing data object " << aonew->path()
	+ << " for " << name());
	+ return aoold;
	+ } else {
	+ MSG_WARNING("Found incompatible pre-existing data object with same path "
	+ << aonew->path() << " for " << name());
	+ }
	+ }
	+ }
	+ MSG_TRACE("Registered " << aonew->annotation("Type") << " " << aonew->path()
	+ << " for " << name());
	+ addAnalysisObject(aonew);
	+ return aonew;
	+ }
	+
	/// 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");
	}

	/// 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/AnalysisHandler.hh b/include/Rivet/AnalysisHandler.hh
	--- a/include/Rivet/AnalysisHandler.hh
	+++ b/include/Rivet/AnalysisHandler.hh
	@@ -1,308 +1,316 @@
	// -- 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 { return _runname; }


	/// 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 { return _eventcounter.numEntries(); }

	/// @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(); }

	/// @brief Compatibility alias for sumOfWeights
	///
	/// @deprecated Prefer sumW
	double sumOfWeights() const { return sumW(); }

	/// @brief Set the sum of weights
	///
	/// This is useful if Rivet is steered externally and
	/// the analyses are run for a sub-contribution of the events
	/// (but of course have to be normalised to the total sum of weights)
	///
	/// @todo What about the sumW2 term? That needs to be set coherently. Need a
	/// new version, with all three N,sumW,sumW2 numbers (or a counter)
	/// supplied.
	///
	/// @deprecated Weight sums are no longer tracked this way...
	void setSumOfWeights(const double& sum) {
	//_sumOfWeights = sum;
	throw Error("Can't set sum of weights independently, since it's now tracked by a Counter. "
	"Please contact the Rivet authors if you need this.");
	}


	/// Is cross-section information required by at least one child analysis?
	/// @deprecated Should no-longer be an issue: does any generator not write the cross-section?
	bool needCrossSection() const;
	/// Whether the handler knows about a cross-section.
	/// @deprecated Should no-longer be an issue: does any generator not write the cross-section?
	bool hasCrossSection() const;

	/// Get the cross-section known to the handler.
	double crossSection() const { return _xs; }

	/// Set the cross-section for the process being generated.
	/// @todo What about the xsec uncertainty? Add a second, optional argument?
	AnalysisHandler& setCrossSection(double 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<AnalysisObjectPtr>& aos);

	/// Read analysis plots into the histo collection (via addData) from the named file.
	void readData(const std::string& filename);

	/// Get all analyses' plots as a vector of analysis objects.
	std::vector<AnalysisObjectPtr> getData(bool includeorphans = false,
	bool includetmps = false) 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 anslysis. The resulting YODA file can then be rwitten
	- /// out by writeData().
	- void mergeYodas(const vector<string> & aofiles, bool equiv = false);
	+ /// 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(AnalysisObjectPtr ao,
	+ const vector<string> & delopts) const;

	//@}

	private:

	/// 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<AnalysisObjectPtr> _orphanedPreloads;

	/// A vector containing copies of analysis objects after
	/// finalize() has been run.
	vector<AnalysisObjectPtr> _finalizedAOs;

	/// @name Run properties
	//@{

	/// Run name
	std::string _runname;

	/// Event counter
	Counter _eventcounter;

	/// Cross-section known to AH
	double _xs, _xserr;

	/// 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;

	/// 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/Tools/RivetYODA.hh b/include/Rivet/Tools/RivetYODA.hh
	--- a/include/Rivet/Tools/RivetYODA.hh
	+++ b/include/Rivet/Tools/RivetYODA.hh
	@@ -1,103 +1,125 @@
	#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"

	namespace Rivet {

	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;

	using YODA::AnalysisObject;
	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, AnalysisObjectPtr> getRefData(const string& papername);

	/// 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.
	+ /// 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>
	- bool aocopy(AnalysisObjectPtr src, AnalysisObjectPtr dst) {
	+ inline bool aocopy(AnalysisObjectPtr src, 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.
	+ /// 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>
	- bool aoadd(AnalysisObjectPtr dst, AnalysisObjectPtr src, double scale) {
	+ inline bool aoadd(AnalysisObjectPtr dst, 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.
	+ /// 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(AnalysisObjectPtr src, 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.
	+ /// 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(AnalysisObjectPtr dst, AnalysisObjectPtr src, double scale);

	+ /// Check if two analysis objects have the same binning or, if not
	+ /// binned, are in other ways compatible.
	+ // inline bool bookingCompatible(CounterPtr, CounterPtr) {
	+ // return true;
	+ // }
	+ 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/pyext/rivet/core.pyx b/pyext/rivet/core.pyx
	--- a/pyext/rivet/core.pyx
	+++ b/pyext/rivet/core.pyx
	@@ -1,238 +1,238 @@
	# distutils: language = c++

	cimport rivet as c
	from cython.operator cimport dereference as deref
	# Need to be careful with memory management -- perhaps use the base object that
	# we used in YODA?

	cdef extern from "<utility>" namespace "std" nogil:
	cdef c.unique_ptr[c.Analysis] move(c.unique_ptr[c.Analysis])

	cdef class AnalysisHandler:
	cdef c.AnalysisHandler *_ptr

	def __cinit__(self):
	self._ptr = new c.AnalysisHandler()

	def __del__(self):
	del self._ptr

	def setIgnoreBeams(self, ignore=True):
	self._ptr.setIgnoreBeams(ignore)

	def addAnalysis(self, name):
	self._ptr.addAnalysis(name.encode('utf-8'))
	return self

	def analysisNames(self):
	anames = self._ptr.analysisNames()
	return [ a.decode('utf-8') for a in anames ]

	# def analysis(self, aname):
	# cdef c.Analysis* ptr = self._ptr.analysis(aname)
	# cdef Analysis pyobj = Analysis.__new__(Analysis)
	# if not ptr:
	# return None
	# pyobj._ptr = ptr
	# return pyobj

	def readData(self, name):
	self._ptr.readData(name.encode('utf-8'))

	def writeData(self, name):
	self._ptr.writeData(name.encode('utf-8'))

	def crossSection(self):
	return self._ptr.crossSection()

	def finalize(self):
	self._ptr.finalize()

	def dump(self, file, period):
	self._ptr.dump(file, period)

	- def mergeYodas(self, filelist, equiv):
	- self._ptr.mergeYodas(filelist, equiv)
	+ def mergeYodas(self, filelist, delopts, equiv):
	+ self._ptr.mergeYodas(filelist, delopts, equiv)


	cdef class Run:
	cdef c.Run *_ptr

	def __cinit__(self, AnalysisHandler h):
	self._ptr = new c.Run(h._ptr[0])

	def __del__(self):
	del self._ptr

	def setCrossSection(self, double x):
	self._ptr.setCrossSection(x)
	return self

	def setListAnalyses(self, choice):
	self._ptr.setListAnalyses(choice)
	return self

	def init(self, name, weight=1.0):
	return self._ptr.init(name.encode('utf-8'), weight)

	def openFile(self, name, weight=1.0):
	return self._ptr.openFile(name.encode('utf-8'), weight)

	def readEvent(self):
	return self._ptr.readEvent()

	def skipEvent(self):
	return self._ptr.skipEvent()

	def processEvent(self):
	return self._ptr.processEvent()

	def finalize(self):
	return self._ptr.finalize()


	cdef class Analysis:
	cdef c.unique_ptr[c.Analysis] _ptr

	def __init__(self):
	raise RuntimeError('This class cannot be instantiated')

	def requiredBeams(self):
	return deref(self._ptr).requiredBeams()

	def requiredEnergies(self):
	return deref(self._ptr).requiredEnergies()

	def keywords(self):
	kws = deref(self._ptr).keywords()
	return [ k.decode('utf-8') for k in kws ]

	def authors(self):
	auths = deref(self._ptr).authors()
	return [ a.decode('utf-8') for a in auths ]

	def bibKey(self):
	return deref(self._ptr).bibKey().decode('utf-8')

	def name(self):
	return deref(self._ptr).name().decode('utf-8')

	def bibTeX(self):
	return deref(self._ptr).bibTeX().decode('utf-8')

	def references(self):
	refs = deref(self._ptr).references()
	return [ r.decode('utf-8') for r in refs ]

	def collider(self):
	return deref(self._ptr).collider().decode('utf-8')

	def description(self):
	return deref(self._ptr).description().decode('utf-8')

	def experiment(self):
	return deref(self._ptr).experiment().decode('utf-8')

	def inspireId(self):
	return deref(self._ptr).inspireId().decode('utf-8')

	def spiresId(self):
	return deref(self._ptr).spiresId().decode('utf-8')

	def runInfo(self):
	return deref(self._ptr).runInfo().decode('utf-8')

	def status(self):
	return deref(self._ptr).status().decode('utf-8')

	def summary(self):
	return deref(self._ptr).summary().decode('utf-8')

	def year(self):
	return deref(self._ptr).year().decode('utf-8')

	def luminosityfb(self):
	return deref(self._ptr).luminosityfb().decode('utf-8')

	#cdef object
	LEVELS = dict(TRACE = 0, DEBUG = 10, INFO = 20,
	WARN = 30, WARNING = 30, ERROR = 40,
	CRITICAL = 50, ALWAYS = 50)


	cdef class AnalysisLoader:
	@staticmethod
	def analysisNames():
	names = c.AnalysisLoader_analysisNames()
	return [ n.decode('utf-8') for n in names ]


	@staticmethod
	def getAnalysis(name):
	name = name.encode('utf-8')
	cdef c.unique_ptr[c.Analysis] ptr = c.AnalysisLoader_getAnalysis(name)
	cdef Analysis pyobj = Analysis.__new__(Analysis)
	if not ptr:
	return None
	pyobj._ptr = move(ptr)
	# Create python object
	return pyobj


	def getAnalysisLibPaths():
	ps = c.getAnalysisLibPaths()
	return [ p.decode('utf-8') for p in ps ]

	def setAnalysisLibPaths(xs):
	bs = [ x.encode('utf-8') for x in xs ]
	c.setAnalysisLibPaths(bs)

	def addAnalysisLibPath(path):
	c.addAnalysisLibPath(path.encode('utf-8'))


	def setAnalysisDataPaths(xs):
	bs = [ x.encode('utf-8') for x in xs ]
	c.setAnalysisDataPaths(bs)

	def addAnalysisDataPath(path):
	c.addAnalysisDataPath(path.encode('utf-8'))

	def getAnalysisDataPaths():
	ps = c.getAnalysisDataPaths()
	return [ p.decode('utf-8') for p in ps ]

	def findAnalysisDataFile(q):
	f = c.findAnalysisDataFile(q.encode('utf-8'))
	return f.decode('utf-8')

	def getAnalysisRefPaths():
	ps = c.getAnalysisRefPaths()
	return [ p.decode('utf-8') for p in ps ]

	def findAnalysisRefFile(q):
	f = c.findAnalysisRefFile(q.encode('utf-8'))
	return f.decode('utf-8')


	def getAnalysisInfoPaths():
	ps = c.getAnalysisInfoPaths()
	return [ p.decode('utf-8') for p in ps ]

	def findAnalysisInfoFile(q):
	f = c.findAnalysisInfoFile(q.encode('utf-8'))
	return f.decode('utf-8')

	def getAnalysisPlotPaths():
	ps = c.getAnalysisPlotPaths()
	return [ p.decode('utf-8') for p in ps ]

	def findAnalysisPlotFile(q):
	f = c.findAnalysisPlotFile(q.encode('utf-8'))
	return f.decode('utf-8')

	def version():
	return c.version().decode('utf-8')

	def setLogLevel(name, level):
	c.setLogLevel(name.encode('utf-8'), level)
	diff --git a/pyext/rivet/rivet.pxd b/pyext/rivet/rivet.pxd
	--- a/pyext/rivet/rivet.pxd
	+++ b/pyext/rivet/rivet.pxd
	@@ -1,94 +1,94 @@
	from libcpp.map cimport map
	from libcpp.pair cimport pair
	from libcpp.vector cimport vector
	from libcpp cimport bool
	from libcpp.string cimport string
	from libcpp.memory cimport unique_ptr

	ctypedef int PdgId
	ctypedef pair[PdgId,PdgId] PdgIdPair

	cdef extern from "Rivet/AnalysisHandler.hh" namespace "Rivet":
	cdef cppclass AnalysisHandler:
	void setIgnoreBeams(bool)
	AnalysisHandler& addAnalysis(string)
	vector[string] analysisNames() const
	# Analysis* analysis(string)
	void writeData(string&)
	void readData(string&)
	double crossSection()
	void finalize()
	void dump(string, int)
	- void mergeYodas(vector[string], bool)
	+ void mergeYodas(vector[string], vector[string], bool)

	cdef extern from "Rivet/Run.hh" namespace "Rivet":
	cdef cppclass Run:
	Run(AnalysisHandler)
	Run& setCrossSection(double) # For chaining?
	Run& setListAnalyses(bool)
	bool init(string, double) except + # $2=1.0
	bool openFile(string, double) except + # $2=1.0
	bool readEvent() except +
	bool skipEvent() except +
	bool processEvent() except +
	bool finalize() except +

	cdef extern from "Rivet/Analysis.hh" namespace "Rivet":
	cdef cppclass Analysis:
	vector[PdgIdPair]& requiredBeams()
	vector[pair[double, double]] requiredEnergies()
	vector[string] authors()
	vector[string] references()
	vector[string] keywords()
	string name()
	string bibTeX()
	string bibKey()
	string collider()
	string description()
	string experiment()
	string inspireId()
	string spiresId()
	string runInfo()
	string status()
	string summary()
	string year()
	string luminosityfb()

	# Might need to translate the following errors, although I believe 'what' is now
	# preserved. But often, we need the exception class name.
	#Error
	#RangeError
	#LogicError
	#PidError
	#InfoError
	#WeightError
	#UserError

	cdef extern from "Rivet/AnalysisLoader.hh":
	vector[string] AnalysisLoader_analysisNames "Rivet::AnalysisLoader::analysisNames" ()
	unique_ptr[Analysis] AnalysisLoader_getAnalysis "Rivet::AnalysisLoader::getAnalysis" (string)

	cdef extern from "Rivet/Tools/RivetPaths.hh" namespace "Rivet":
	vector[string] getAnalysisLibPaths()
	void setAnalysisLibPaths(vector[string])
	void addAnalysisLibPath(string)

	vector[string] getAnalysisDataPaths()
	void setAnalysisDataPaths(vector[string])
	void addAnalysisDataPath(string)
	string findAnalysisDataFile(string)

	vector[string] getAnalysisRefPaths()
	string findAnalysisRefFile(string)

	vector[string] getAnalysisInfoPaths()
	string findAnalysisInfoFile(string)

	vector[string] getAnalysisPlotPaths()
	string findAnalysisPlotFile(string)

	cdef extern from "Rivet/Rivet.hh" namespace "Rivet":
	string version()

	cdef extern from "Rivet/Tools/Logging.hh":
	void setLogLevel "Rivet::Log::setLevel" (string, int)
	diff --git a/src/Core/Analysis.cc b/src/Core/Analysis.cc
	--- a/src/Core/Analysis.cc
	+++ b/src/Core/Analysis.cc
	@@ -1,1009 +1,925 @@
	// -- 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) < 1GeV && abs(ep.second - energies.second) < 1GeV) \|\|
	(abs(ep.first - energies.second) < 1GeV && abs(ep.second - energies.first) < 1GeV)) {
	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;
	+ return addOrGetCompatAO(make_shared<Counter>(histoPath(cname), title));
	}


	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);
	+ return bookCounter(mkAxisCode(datasetId, xAxisId, yAxisId), 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 = 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 (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;
	+ return addOrGetCompatAO(hist);
	}


	Histo1DPtr Analysis::bookHisto1D(const string& hname,
	const vector<double>& binedges,
	const string& title,
	const string& xtitle,
	const string& ytitle) {
	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 (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;
	+ return addOrGetCompatAO(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 = 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 (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;
	+ return addOrGetCompatAO(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;
	+ return addOrGetCompatAO(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;
	+ return addOrGetCompatAO(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;
	+ return addOrGetCompatAO(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;
	+ return addOrGetCompatAO(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;
	+ return addOrGetCompatAO(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;
	+ return addOrGetCompatAO(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;
	+ return addOrGetCompatAO(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;
	+ return addOrGetCompatAO(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;
	+ return addOrGetCompatAO(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;
	+ return addOrGetCompatAO(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());
	+ Scatter2DPtr 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);
	}
	s->setTitle(title);
	s->setAnnotation("XLabel", xtitle);
	s->setAnnotation("YLabel", ytitle);
	- return s;
	+ return addOrGetCompatAO(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;
	+ return addOrGetCompatAO(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);

	}

	}
	diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc
	--- a/src/Core/AnalysisHandler.cc
	+++ b/src/Core/AnalysisHandler.cc
	@@ -1,566 +1,571 @@
	// -- 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/IO.h"

	namespace Rivet {


	AnalysisHandler::AnalysisHandler(const string& runname)
	: _runname(runname),
	_eventcounter("/_EVTCOUNT"),
	_xs(NAN), _xserr(NAN),
	_initialised(false), _ignoreBeams(false), _dumpPeriod(0), _dumping(false)
	{ }


	AnalysisHandler::~AnalysisHandler()
	{ }


	Log& AnalysisHandler::getLog() const {
	return Log::getLog("Rivet.Analysis.Handler");
	}


	void AnalysisHandler::init(const GenEvent& ge) {
	if (_initialised)
	throw UserError("AnalysisHandler::init has already been called: cannot re-initialize!");

	setRunBeams(Rivet::beams(ge));
	MSG_DEBUG("Initialising the analysis handler");
	_eventcounter.reset();

	// 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
	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());
	}
	_initialised = true;
	MSG_DEBUG("Analysis handler initialised");
	}


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

	// Weights
	/// @todo Drop this / just report first weight when we support multiweight events
	_eventcounter.fill(event.weight());
	MSG_DEBUG("Event #" << _eventcounter.numEntries() << " weight = " << event.weight());

	// Cross-section
	#ifdef HEPMC_HAS_CROSS_SECTION
	if (ge.cross_section()) {
	_xs = ge.cross_section()->cross_section();
	_xserr = ge.cross_section()->cross_section_error();
	}
	#endif

	// 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;

	// First we make copies of all analysis objects.
	map<string,AnalysisObjectPtr> backupAOs;
	for (auto ao : getData(false, true) )
	backupAOs[ao->path()] = AnalysisObjectPtr(ao->newclone());

	// Now we run the (re-entrant) finalize() functions for all analyses.
	MSG_INFO("Finalising analyses");
	for (AnaHandle a : _analyses) {
	a->setCrossSection(_xs);
	try {
	if ( !_dumping \|\| a->info().reentrant() ) a->finalize();
	else if ( _dumping )
	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 : getData() )
	_finalizedAOs.push_back(AnalysisObjectPtr(ao->newclone()));
	for ( auto ao : getData(false, true) ) {
	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 = _eventcounter.numEntries();
	MSG_INFO("Processed " << nevts << " event" << (nevts != 1 ? "s" : ""));

	// // Delete analyses
	// MSG_DEBUG("Deleting analyses");
	// _analyses.clear();

	// Print out MCnet boilerplate
	cout << endl;
	cout << "The MCnet usage guidelines apply to Rivet: see http://www.montecarlonet.org/GUIDELINES" << endl;
	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<AnalysisObjectPtr>& aos) {
	for (const 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);
	a->addAnalysisObject(ao); /// @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,
	+ 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, bool equiv) {
	+ mergeYodas(const vector<string> & aofiles, const vector<string> & delopts, bool equiv) {
	vector< vector<AnalysisObjectPtr> > aosv;
	vector<double> xsecs;
	vector<double> xsecerrs;
	vector<CounterPtr> sows;
	set<string> ananames;
	_eventcounter.reset();

	// First scan all files and extract analysis objects and add the
	// corresponding anayses..
	for ( auto 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;
	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);
	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];
	- // HERE we shoud handle merged options, if any.
	- // vector<string> anaopts = split(ananame, ":");
	- // ananame = anaopts[0];
	- // for (int i = 1, N = anaopts.size(); i < N; ++i )
	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());
	xsecerrs.push_back(sqr(xsec->point(0).xErrAvg()));
	- std::cerr << _eventcounter.numEntries() << std::endl;
	- std::cerr << _eventcounter.effNumEntries() << std::endl;
	- std::cerr << _eventcounter.sumW() << std::endl;
	- std::cerr << _eventcounter.sumW2() << std::endl;
	_eventcounter += *sow;
	- std::cerr << _eventcounter.numEntries() << std::endl;
	- std::cerr << _eventcounter.effNumEntries() << std::endl;
	- std::cerr << _eventcounter.sumW() << std::endl;
	- std::cerr << _eventcounter.sumW2() << std::endl;
	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();
	} else {
	_xserr = sqrt(_xserr);
	for ( int i = 0, N = sows.size(); i < N; ++i )
	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());
	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<AnalysisObjectPtr> aos;
	try {
	/// @todo Use new YODA SFINAE to fill the smart ptr vector directly
	vector<YODA::AnalysisObject*> aos_raw;
	YODA::read(filename, aos_raw);
	for (AnalysisObject* aor : aos_raw) aos.push_back(AnalysisObjectPtr(aor));
	} catch (...) { //< YODA::ReadError&
	throw UserError("Unexpected error in reading file: " + filename);
	}
	if (!aos.empty()) addData(aos);
	}


	vector<AnalysisObjectPtr> AnalysisHandler::
	getData(bool includeorphans, bool includetmps) const {
	vector<AnalysisObjectPtr> rtn;
	// Event counter
	rtn.push_back( make_shared<Counter>(_eventcounter) );
	// Cross-section + err as scatter
	YODA::Scatter1D::Points pts; pts.insert(YODA::Point1D(_xs, _xserr));
	rtn.push_back( make_shared<Scatter1D>(pts, "/_XSEC") );
	// Analysis histograms
	for (const AnaHandle a : analyses()) {
	vector<AnalysisObjectPtr> aos = a->analysisObjects();
	// MSG_WARNING(a->name() << " " << aos.size());
	for (const AnalysisObjectPtr ao : aos) {
	// Exclude paths from final write-out if they contain a "TMP" layer (i.e. matching "/TMP/")
	/// @todo This needs to be much more nuanced for re-entrant histogramming
	if ( !includetmps && ao->path().find("/TMP/" ) != string::npos) continue;
	rtn.push_back(ao);
	}
	}
	// Sort histograms alphanumerically by path before write-out
	sort(rtn.begin(), rtn.end(), [](AnalysisObjectPtr a, AnalysisObjectPtr b) {return a->path() < b->path();});
	if ( includeorphans )
	rtn.insert(rtn.end(), _orphanedPreloads.begin(), _orphanedPreloads.end());
	return rtn;
	}


	void AnalysisHandler::writeData(const string& filename) const {
	vector<AnalysisObjectPtr> out = _finalizedAOs;
	out.reserve(2*out.size());
	vector<AnalysisObjectPtr> aos = getData(false, true);
	for ( auto ao : aos ) {
	ao = AnalysisObjectPtr(ao->newclone());
	ao->setPath("/RAW" + ao->path());
	out.push_back(ao);
	}

	try {
	YODA::write(filename, out.begin(), out.end());
	} catch (...) { //< YODA::WriteError&
	throw UserError("Unexpected error in writing file: " + filename);
	}
	}


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


	bool AnalysisHandler::needCrossSection() const {
	bool rtn = false;
	for (const AnaHandle a : _analyses) {
	if (!rtn) rtn = a->needsCrossSection();
	if (rtn) break;
	}
	return rtn;
	}


	AnalysisHandler& AnalysisHandler::setCrossSection(double xs) {
	_xs = xs;
	return *this;
	}


	bool AnalysisHandler::hasCrossSection() const {
	return (!std::isnan(crossSection()));
	}


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


	}

File Metadata

Mime Type: text/x-diff
Expires: Sun, Feb 23, 2:45 PM (22 h, 21 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4486693
Default Alt Text: (138 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions