diff --git a/analyses/pluginMC/TEST.cc b/analyses/pluginMC/TEST.cc
new file mode 100644
--- /dev/null
+++ b/analyses/pluginMC/TEST.cc
@@ -0,0 +1,80 @@
+// -*- C++ -*-
+#include "Rivet/Analysis.hh"
+#include "Rivet/Projections/ChargedFinalState.hh"
+#include "Rivet/Projections/PrimaryParticles.hh"
+#include "Rivet/Tools/Correlators.hh"
+
+
+namespace Rivet {
+
+
+  class TEST : public CumulantAnalysis {
+  public:
+
+    /// @name Constructors etc.
+    //@{
+
+    /// Constructor
+    TEST() : CumulantAnalysis("TEST") {
+    }
+    //@}
+
+
+  public:
+
+    /// @name Analysis methods
+    //@{
+    /// Book histograms and initialise projections before the run
+    void init() {
+
+      ChargedFinalState cfs(-1.0, 1.0);
+      declare(cfs, "CFS");
+      ChargedFinalState pp(Cuts::abseta < 2.0);
+      declare(pp, "PP");
+      h_c22 = bookScatter2D("c22",120,0,120);
+      h_c23 = bookScatter2D("c23",120,0,120);
+      ec22 = bookECorrelator<2,2>("ec22",h_c22);
+      ec23 = bookECorrelator<3,2>("ec32",h_c22);
+      pair<int, int> max = getMaxValues(); 
+      // Declare correlator projections.
+      declare(Correlators(pp, max.first, max.second),"CRS");
+    }
+
+    /// Perform the per-event analysis
+    void analyze(const Event& event) {
+      ec22->fill(apply<ChargedFinalState>(event,"CFS").particles().size(), 
+        apply<Correlators>(event,"CRS"), event.weight());
+      ec23->fill(apply<ChargedFinalState>(event,"CFS").particles().size(), 
+        apply<Correlators>(event,"CRS"), event.weight());
+    }
+
+
+    /// Normalise histograms etc., after the run
+    void finalize() {
+      CumulantAnalysis::finalize();
+      cnTwoInt(h_c22,ec22);
+    }
+
+
+    //@}
+
+	  private:
+
+    /// @name Histograms
+    //@{
+
+    Scatter2DPtr h_c22;
+    ECorrPtr ec22;
+    Scatter2DPtr h_c23;
+    ECorrPtr ec23;
+    //@}
+
+
+  };
+
+
+
+  // The hook for the plugin system
+  DECLARE_RIVET_PLUGIN(TEST);
+
+}
diff --git a/analyses/pluginMC/TEST.info b/analyses/pluginMC/TEST.info
new file mode 100644
--- /dev/null
+++ b/analyses/pluginMC/TEST.info
@@ -0,0 +1,28 @@
+Name: TEST
+Year: 2010
+Summary: Pseudorapidities at three energies, charged multiplicity at 7 TeV
+Experiment: ALICE
+RunInfo:
+  Diffractive events need to be enabled.
+NumEvents: 1000000
+Beams: [p+, p+]
+Energies: [900]
+Description:
+  'This is an ALICE publication with pseudorapities for 0.9, 2.36 and $7\;\TeV$
+  and the charged multiplicity at $7\;\TeV$. The analysis requires at least on
+  charged particle in the event. Only the INEL distributions are considered here'
+BibKey: Aamodt:2010pp
+BibTeX: '@article{Aamodt:2010pp,
+      author         = "Aamodt, K. and others",
+      title          = "{Charged-particle multiplicity measurement in
+                        proton-proton collisions at $\sqrt{s} = 7$ TeV with ALICE at LHC}",
+      collaboration  = "ALICE",
+      journal        = "Eur.Phys.J.",
+      volume         = "C68",
+      pages          = "345-354",
+      doi            = "10.1140/epjc/s10052-010-1350-2",
+      year           = "2010",
+      eprint         = "1004.3514",
+      archivePrefix  = "arXiv",
+      primaryClass   = "hep-ex",
+}'
diff --git a/analyses/pluginMC/TEST.plot b/analyses/pluginMC/TEST.plot
new file mode 100644
--- /dev/null
+++ b/analyses/pluginMC/TEST.plot
@@ -0,0 +1,6 @@
+# BEGIN PLOT /TEST/d04-x01-y01
+Title=Pseudorapidity $\sqrt(s)=0.9$ TeV, INEL $>0$
+XLabel=$\eta$
+YLabel=$\text{d}N/\text{d}\eta$
+LogY=0
+# END PLOT
diff --git a/analyses/pluginMC/TEST.yoda b/analyses/pluginMC/TEST.yoda
new file mode 100644
--- /dev/null
+++ b/analyses/pluginMC/TEST.yoda
@@ -0,0 +1,18 @@
+BEGIN YODA_SCATTER2D /REF/TEST/d04-x01-y01
+IsRef=1
+Path=/REF/TEST/d04-x01-y01
+Title=
+Type=Scatter2D
+# xval	 xerr-	 xerr+	 yval	 yerr-	 yerr+
+-9.000000e-01	1.000000e-01	1.000000e-01	4.000000e+00	7.615773e-02	7.615773e-02
+-7.000000e-01	1.000000e-01	1.000000e-01	3.870000e+00	7.615773e-02	7.615773e-02
+-5.000000e-01	1.000000e-01	1.000000e-01	3.800000e+00	7.615773e-02	7.615773e-02
+-3.000000e-01	1.000000e-01	1.000000e-01	3.700000e+00	6.324555e-02	6.324555e-02
+-1.000000e-01	1.000000e-01	1.000000e-01	3.670000e+00	6.324555e-02	6.324555e-02
+1.000000e-01	1.000000e-01	1.000000e-01	3.730000e+00	6.324555e-02	6.324555e-02
+3.000000e-01	1.000000e-01	1.000000e-01	3.720000e+00	6.708204e-02	6.708204e-02
+5.000000e-01	1.000000e-01	1.000000e-01	3.770000e+00	7.615773e-02	7.615773e-02
+7.000000e-01	1.000000e-01	1.000000e-01	3.920000e+00	7.615773e-02	7.615773e-02
+9.000000e-01	1.000000e-01	1.000000e-01	4.010000e+00	7.615773e-02	7.615773e-02
+END YODA_SCATTER2D
+
diff --git a/include/Rivet/Analysis.hh b/include/Rivet/Analysis.hh
--- a/include/Rivet/Analysis.hh
+++ b/include/Rivet/Analysis.hh
@@ -1,1226 +1,1225 @@
 // -*- C++ -*-
 #ifndef RIVET_Analysis_HH
 #define RIVET_Analysis_HH
 
 #include "Rivet/Config/RivetCommon.hh"
 #include "Rivet/AnalysisInfo.hh"
 #include "Rivet/Event.hh"
 #include "Rivet/Projection.hh"
 #include "Rivet/ProjectionApplier.hh"
 #include "Rivet/ProjectionHandler.hh"
 #include "Rivet/AnalysisLoader.hh"
 #include "Rivet/Tools/Cuts.hh"
 #include "Rivet/Tools/Logging.hh"
 #include "Rivet/Tools/ParticleUtils.hh"
 #include "Rivet/Tools/BinnedHistogram.hh"
 #include "Rivet/Tools/RivetMT2.hh"
 #include "Rivet/Tools/RivetYODA.hh"
 #include "Rivet/Tools/Percentile.hh"
 #include "Rivet/Projections/CentralityProjection.hh"
 
 
 /// @def vetoEvent
 /// Preprocessor define for vetoing events, including the log message and return.
 #define vetoEvent                                                       \
   do { MSG_DEBUG("Vetoing event on line " << __LINE__ << " of " << __FILE__); return; } while(0)
 
 
 namespace Rivet {
 
 
   // Forward declaration
   class AnalysisHandler;
 
   /// @brief This is the base class of all analysis classes in Rivet.
   ///
   /// There are
   /// three virtual functions which should be implemented in base classes:
   ///
   /// void init() is called by Rivet before a run is started. Here the
   /// analysis class should book necessary histograms. The needed
   /// projections should probably rather be constructed in the
   /// constructor.
   ///
   /// void analyze(const Event&) is called once for each event. Here the
   /// analysis class should apply the necessary Projections and fill the
   /// histograms.
   ///
   /// void finalize() is called after a run is finished. Here the analysis
   /// class should do whatever manipulations are necessary on the
   /// histograms. Writing the histograms to a file is, however, done by
   /// the Rivet class.
   class Analysis : public ProjectionApplier {
 
     /// The AnalysisHandler is a friend.
     friend class AnalysisHandler;
 
-
   public:
 
     /// @name Standard constructors and destructors.
     //@{
 
     // /// The default constructor.
     // Analysis();
 
     /// Constructor
     Analysis(const std::string& name);
 
     /// The destructor.
     virtual ~Analysis() {}
 
     //@}
 
 
   public:
 
     /// @name Main analysis methods
     //@{
 
     /// Initialize this analysis object. A concrete class should here
     /// book all necessary histograms. An overridden function must make
     /// sure it first calls the base class function.
     virtual void init() { }
 
     /// Analyze one event. A concrete class should here apply the
     /// necessary projections on the \a event and fill the relevant
     /// histograms. An overridden function must make sure it first calls
     /// the base class function.
     virtual void analyze(const Event& event) = 0;
 
     /// Finalize this analysis object. A concrete class should here make
     /// all necessary operations on the histograms. Writing the
     /// histograms to a file is, however, done by the Rivet class. An
     /// overridden function must make sure it first calls the base class
     /// function.
     virtual void finalize() { }
 
     //@}
 
 
   public:
 
     /// @name Metadata
     /// Metadata is used for querying from the command line and also for
     /// building web pages and the analysis pages in the Rivet manual.
     //@{
 
     /// Get the actual AnalysisInfo object in which all this metadata is stored.
     const AnalysisInfo& info() const {
       assert(_info && "No AnalysisInfo object :O");
       return *_info;
     }
 
     /// @brief Get the name of the analysis.
     ///
     /// By default this is computed by combining the results of the
     /// experiment, year and Spires ID metadata methods and you should
     /// only override it if there's a good reason why those won't
     /// work. If options has been set for this instance, a
     /// corresponding string is appended at the end.
     virtual std::string name() const {
       return  ( (info().name().empty()) ? _defaultname : info().name() ) + _optstring;
     }
     // get name of reference data file, which could be different from plugin name
     virtual std::string getRefDataName() const {
       return (info().getRefDataName().empty()) ? _defaultname : info().getRefDataName();
     }
 
     // set name of reference data file, which could be different from plugin name
     virtual void setRefDataName(const std::string& ref_data="") {
       info().setRefDataName(!ref_data.empty() ? ref_data : name());
     }
 
     /// Get the Inspire ID code for this analysis.
     virtual std::string inspireId() const {
       return info().inspireId();
     }
 
     /// Get the SPIRES ID code for this analysis (~deprecated).
     virtual std::string spiresId() const {
       return info().spiresId();
     }
 
     /// @brief Names & emails of paper/analysis authors.
     ///
     /// Names and email of authors in 'NAME \<EMAIL\>' format. The first
     /// name in the list should be the primary contact person.
     virtual std::vector<std::string> authors() const {
       return info().authors();
     }
 
     /// @brief Get a short description of the analysis.
     ///
     /// Short (one sentence) description used as an index entry.
     /// Use @a description() to provide full descriptive paragraphs
     /// of analysis details.
     virtual std::string summary() const {
       return info().summary();
     }
 
     /// @brief Get a full description of the analysis.
     ///
     /// Full textual description of this analysis, what it is useful for,
     /// what experimental techniques are applied, etc. Should be treated
     /// as a chunk of restructuredText (http://docutils.sourceforge.net/rst.html),
     /// with equations to be rendered as LaTeX with amsmath operators.
     virtual std::string description() const {
       return info().description();
     }
 
     /// @brief Information about the events needed as input for this analysis.
     ///
     /// Event types, energies, kinematic cuts, particles to be considered
     /// stable, etc. etc. Should be treated as a restructuredText bullet list
     /// (http://docutils.sourceforge.net/rst.html)
     virtual std::string runInfo() const {
       return info().runInfo();
     }
 
     /// Experiment which performed and published this analysis.
     virtual std::string experiment() const {
       return info().experiment();
     }
 
     /// Collider on which the experiment ran.
     virtual std::string collider() const {
       return info().collider();
     }
 
     /// When the original experimental analysis was published.
     virtual std::string year() const {
       return info().year();
     }
 
     /// The luminosity in inverse femtobarn
     virtual std::string luminosityfb() const {
       return info().luminosityfb();
     }
 
     /// Journal, and preprint references.
     virtual std::vector<std::string> references() const {
       return info().references();
     }
 
     /// BibTeX citation key for this article.
     virtual std::string bibKey() const {
       return info().bibKey();
     }
 
     /// BibTeX citation entry for this article.
     virtual std::string bibTeX() const {
       return info().bibTeX();
     }
 
     /// Whether this analysis is trusted (in any way!)
     virtual std::string status() const {
       return (info().status().empty()) ? "UNVALIDATED" : info().status();
     }
 
     /// Any work to be done on this analysis.
     virtual std::vector<std::string> todos() const {
       return info().todos();
     }
 
 
     /// Return the allowed pairs of incoming beams required by this analysis.
     virtual const std::vector<PdgIdPair>& requiredBeams() const {
       return info().beams();
     }
     /// Declare the allowed pairs of incoming beams required by this analysis.
     virtual Analysis& setRequiredBeams(const std::vector<PdgIdPair>& requiredBeams) {
       info().setBeams(requiredBeams);
       return *this;
     }
 
 
     /// Sets of valid beam energy pairs, in GeV
     virtual const std::vector<std::pair<double, double> >& requiredEnergies() const {
       return info().energies();
     }
 
     /// Get vector of analysis keywords
     virtual const std::vector<std::string> & keywords() const {
       return info().keywords();
     }
 
     /// Declare the list of valid beam energy pairs, in GeV
     virtual Analysis& setRequiredEnergies(const std::vector<std::pair<double, double> >& requiredEnergies) {
       info().setEnergies(requiredEnergies);
       return *this;
     }
 
 
     /// Return true if this analysis needs to know the process cross-section.
     /// @todo Remove this and require HepMC >= 2.06
     bool needsCrossSection() const {
       return info().needsCrossSection();
     }
     /// Declare whether this analysis needs to know the process cross-section from the generator.
     /// @todo Remove this and require HepMC >= 2.06
     Analysis& setNeedsCrossSection(bool needed=true) {
       info().setNeedsCrossSection(needed);
       return *this;
     }
 
     //@}
 
 
     /// @name Internal metadata modifying methods
     //@{
 
     /// Get the actual AnalysisInfo object in which all this metadata is stored (non-const).
     AnalysisInfo& info() {
       assert(_info && "No AnalysisInfo object :O");
       return *_info;
     }
 
     //@}
 
 
     /// @name Run conditions
     //@{
 
     /// Incoming beams for this run
     const ParticlePair& beams() const;
 
     /// Incoming beam IDs for this run
     const PdgIdPair beamIds() const;
 
     /// Centre of mass energy for this run
     double sqrtS() const;
 
     //@}
 
 
     /// @name Analysis / beam compatibility testing
     //@{
 
     /// Check if analysis is compatible with the provided beam particle IDs and energies
     bool isCompatible(const ParticlePair& beams) const;
 
     /// Check if analysis is compatible with the provided beam particle IDs and energies
     bool isCompatible(PdgId beam1, PdgId beam2, double e1, double e2) const;
 
     /// Check if analysis is compatible with the provided beam particle IDs and energies
     bool isCompatible(const PdgIdPair& beams, const std::pair<double,double>& energies) const;
 
     //@}
 
 
     /// Set the cross section from the generator
     Analysis& setCrossSection(double xs); //, double xserr=0.0);
 
     /// Access the controlling AnalysisHandler object.
     AnalysisHandler& handler() const { return *_analysishandler; }
 
 
 
   protected:
 
     /// Get a Log object based on the name() property of the calling analysis object.
     Log& getLog() const;
 
     /// Get the process cross-section in pb. Throws if this hasn't been set.
     double crossSection() const;
 
     /// Get the process cross-section per generated event in pb. Throws if this
     /// hasn't been set.
     double crossSectionPerEvent() const;
 
     /// @brief Get the number of events seen (via the analysis handler).
     ///
     /// @note Use in the finalize phase only.
     size_t numEvents() const;
 
     /// @brief Get the sum of event weights seen (via the analysis handler).
     ///
     /// @note Use in the finalize phase only.
     double sumW() const;
     /// Alias
     double sumOfWeights() const { return sumW(); }
 
     /// @brief Get the sum of squared event weights seen (via the analysis handler).
     ///
     /// @note Use in the finalize phase only.
     double sumW2() const;
 
 
   protected:
 
     /// @name Histogram paths
     //@{
 
     /// Get the canonical histogram "directory" path for this analysis.
     const std::string histoDir() const;
 
     /// Get the canonical histogram path for the named histogram in this analysis.
     const std::string histoPath(const std::string& hname) const;
 
     /// Get the canonical histogram path for the numbered histogram in this analysis.
     const std::string histoPath(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const;
 
     /// Get the internal histogram name for given d, x and y (cf. HepData)
     const std::string mkAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const;
     /// Alias
     /// @deprecated Prefer the "mk" form, consistent with other "making function" names
     const std::string makeAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
       return mkAxisCode(datasetId, xAxisId, yAxisId);
     }
 
     //@}
 
 
     /// @name Histogram reference data
     //@{
 
     /// Get reference data for a named histo
     /// @todo SFINAE to ensure that the type inherits from YODA::AnalysisObject?
     template <typename T=YODA::Scatter2D>
     const T& refData(const string& hname) const {
       _cacheRefData();
       MSG_TRACE("Using histo bin edges for " << name() << ":" << hname);
       if (!_refdata[hname]) {
         MSG_ERROR("Can't find reference histogram " << hname);
         throw Exception("Reference data " + hname + " not found.");
       }
       return dynamic_cast<T&>(*_refdata[hname]);
     }
 
     /// Get reference data for a numbered histo
     /// @todo SFINAE to ensure that the type inherits from YODA::AnalysisObject?
     template <typename T=YODA::Scatter2D>
     const T& refData(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const {
       const string hname = makeAxisCode(datasetId, xAxisId, yAxisId);
       return refData(hname);
     }
 
     //@}
 
 
     /// @name Counter booking
     //@{
 
     /// Book a counter.
     CounterPtr bookCounter(const std::string& name,
                            const std::string& title="");
                            // const std::string& valtitle=""
 
     /// Book a counter, using a path generated from the dataset and axis ID codes
     ///
     /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
     CounterPtr bookCounter(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                            const std::string& title="");
                            // const std::string& valtitle=""
 
     //@}
 
 
     /// @name 1D histogram booking
     //@{
 
     /// Book a 1D histogram with @a nbins uniformly distributed across the range @a lower - @a upper .
     Histo1DPtr bookHisto1D(const std::string& name,
                            size_t nbins, double lower, double upper,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="");
 
     /// Book a 1D histogram with non-uniform bins defined by the vector of bin edges @a binedges .
     Histo1DPtr bookHisto1D(const std::string& name,
                            const std::vector<double>& binedges,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="");
 
     /// Book a 1D histogram with non-uniform bins defined by the vector of bin edges @a binedges .
     Histo1DPtr bookHisto1D(const std::string& name,
                            const std::initializer_list<double>& binedges,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="");
 
         /// Book a 1D histogram with binning from a reference scatter.
     Histo1DPtr bookHisto1D(const std::string& name,
                            const Scatter2D& refscatter,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="");
 
     /// Book a 1D histogram, using the binnings in the reference data histogram.
     Histo1DPtr bookHisto1D(const std::string& name,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="");
 
     /// Book a 1D histogram, using the binnings in the reference data histogram.
     ///
     /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
     Histo1DPtr bookHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="");
 
     //@}
 
 
     /// @name 2D histogram booking
     //@{
 
     /// Book a 2D histogram with @a nxbins and @a nybins uniformly
     /// distributed across the ranges @a xlower - @a xupper and @a
     /// ylower - @a yupper respectively along the x- and y-axis.
     Histo2DPtr bookHisto2D(const std::string& name,
                            size_t nxbins, double xlower, double xupper,
                            size_t nybins, double ylower, double yupper,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="",
                            const std::string& ztitle="");
 
     /// Book a 2D histogram with non-uniform bins defined by the
     /// vectors of bin edges @a xbinedges and @a ybinedges.
     Histo2DPtr bookHisto2D(const std::string& name,
                            const std::vector<double>& xbinedges,
                            const std::vector<double>& ybinedges,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="",
                            const std::string& ztitle="");
 
     /// Book a 2D histogram with non-uniform bins defined by the
     /// vectors of bin edges @a xbinedges and @a ybinedges.
     Histo2DPtr bookHisto2D(const std::string& name,
                            const std::initializer_list<double>& xbinedges,
                            const std::initializer_list<double>& ybinedges,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="",
                            const std::string& ztitle="");
 
     /// Book a 2D histogram with binning from a reference scatter.
     Histo2DPtr bookHisto2D(const std::string& name,
                            const Scatter3D& refscatter,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="",
                            const std::string& ztitle="");
 
     /// Book a 2D histogram, using the binnings in the reference data histogram.
     Histo2DPtr bookHisto2D(const std::string& name,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="",
                            const std::string& ztitle="");
 
     /// Book a 2D histogram, using the binnings in the reference data histogram.
     ///
     /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
     Histo2DPtr bookHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                            const std::string& title="",
                            const std::string& xtitle="",
                            const std::string& ytitle="",
                            const std::string& ztitle="");
 
     //@}
 
 
     /// @name 1D profile histogram booking
     //@{
 
     /// Book a 1D profile histogram with @a nbins uniformly distributed across the range @a lower - @a upper .
     Profile1DPtr bookProfile1D(const std::string& name,
                                size_t nbins, double lower, double upper,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="");
 
     /// Book a 1D profile histogram with non-uniform bins defined by the vector of bin edges @a binedges .
     Profile1DPtr bookProfile1D(const std::string& name,
                                const std::vector<double>& binedges,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="");
 
     /// Book a 1D profile histogram with non-uniform bins defined by the vector of bin edges @a binedges .
     Profile1DPtr bookProfile1D(const std::string& name,
                                const std::initializer_list<double>& binedges,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="");
 
     /// Book a 1D profile histogram with binning from a reference scatter.
     Profile1DPtr bookProfile1D(const std::string& name,
                                const Scatter2D& refscatter,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="");
 
     /// Book a 1D profile histogram, using the binnings in the reference data histogram.
     Profile1DPtr bookProfile1D(const std::string& name,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="");
 
     /// Book a 1D profile histogram, using the binnings in the reference data histogram.
     ///
     /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
     Profile1DPtr bookProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="");
 
     //@}
 
 
     /// @name 2D profile histogram booking
     //@{
 
     /// Book a 2D profile histogram with @a nxbins and @a nybins uniformly
     /// distributed across the ranges @a xlower - @a xupper and @a ylower - @a
     /// yupper respectively along the x- and y-axis.
     Profile2DPtr bookProfile2D(const std::string& name,
                                size_t nxbins, double xlower, double xupper,
                                size_t nybins, double ylower, double yupper,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="",
                                const std::string& ztitle="");
 
     /// Book a 2D profile histogram with non-uniform bins defined by the vectorx
     /// of bin edges @a xbinedges and @a ybinedges.
     Profile2DPtr bookProfile2D(const std::string& name,
                                const std::vector<double>& xbinedges,
                                const std::vector<double>& ybinedges,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="",
                                const std::string& ztitle="");
 
     /// Book a 2D profile histogram with non-uniform bins defined by the vectorx
     /// of bin edges @a xbinedges and @a ybinedges.
     Profile2DPtr bookProfile2D(const std::string& name,
                                const std::initializer_list<double>& xbinedges,
                                const std::initializer_list<double>& ybinedges,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="",
                                const std::string& ztitle="");
 
     /// Book a 2D profile histogram with binning from a reference scatter.
     Profile2DPtr bookProfile2D(const std::string& name,
                                const Scatter3D& refscatter,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="",
                                const std::string& ztitle="");
 
     /// Book a 2D profile histogram, using the binnings in the reference data histogram.
     Profile2DPtr bookProfile2D(const std::string& name,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="",
                                const std::string& ztitle="");
 
     /// Book a 2D profile histogram, using the binnings in the reference data histogram.
     ///
     /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
     Profile2DPtr bookProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="",
                                const std::string& ztitle="");
 
     //@}
 
 
     /// @name 2D scatter booking
     //@{
 
     /// @brief Book a 2-dimensional data point set with the given name.
     ///
     /// @note Unlike histogram booking, scatter booking by default makes no
     /// attempt to use reference data to pre-fill the data object. If you want
     /// this, which is sometimes useful e.g. when the x-position is not really
     /// meaningful and can't be extracted from the data, then set the @a
     /// copy_pts parameter to true. This creates points to match the reference
     /// data's x values and errors, but with the y values and errors zeroed...
     /// assuming that there is a reference histo with the same name: if there
     /// isn't, an exception will be thrown.
     Scatter2DPtr bookScatter2D(const std::string& name,
                                bool copy_pts=false,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="");
 
     /// @brief Book a 2-dimensional data point set, using the binnings in the reference data histogram.
     ///
     /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way.
     ///
     /// @note Unlike histogram booking, scatter booking by default makes no
     /// attempt to use reference data to pre-fill the data object. If you want
     /// this, which is sometimes useful e.g. when the x-position is not really
     /// meaningful and can't be extracted from the data, then set the @a
     /// copy_pts parameter to true. This creates points to match the reference
     /// data's x values and errors, but with the y values and errors zeroed.
     Scatter2DPtr bookScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId,
                                bool copy_pts=false,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="");
 
     /// @brief Book a 2-dimensional data point set with equally spaced x-points in a range.
     ///
     /// The y values and errors will be set to 0.
     Scatter2DPtr bookScatter2D(const std::string& name,
                                size_t npts, double lower, double upper,
                                const std::string& title="",
                                const std::string& xtitle="",
                                const std::string& ytitle="");
 
     /// @brief Book a 2-dimensional data point set based on provided contiguous "bin edges".
     ///
     /// The y values and errors will be set to 0.
     Scatter2DPtr bookScatter2D(const std::string& hname,
                                const std::vector<double>& binedges,
                                const std::string& title,
                                const std::string& xtitle,
                                const std::string& ytitle);
 
     //@}
 
 
   public:
 
     /// @name accessing options for this Analysis instance.
     //@{
 
     /// Get an option for this analysis instance as a string.
     std::string getOption(std::string optname) {
       if ( _options.find(optname) != _options.end() )
         return _options.find(optname)->second;
       return "";
     }
 
     /// Get an option for this analysis instance converted to a
     /// specific type (given by the specified @a def value).
     template<typename T>
     T getOption(std::string optname, T def) {
       if (_options.find(optname) == _options.end()) return def;
       std::stringstream ss;
       ss << _options.find(optname)->second;
       T ret;
       ss >> ret;
       return ret;
     }
 
     //@}
     /// @brief Book a CentralityProjection
     ///
     /// Using a SingleValueProjection, @a proj, giving the value of an
     /// experimental observable to be used as a centrality estimator,
     /// book a CentralityProjection based on the experimentally
     /// measured pecentiles of this observable (as given by the
     /// reference data for the @a calHistName histogram in the @a
     /// calAnaName analysis. If a preloaded file with the output of a
     /// run using the @a calAnaName analysis contains a valid
     /// generated @a calHistName histogram, it will be used as an
     /// optional percentile binning. Also if this preloaded file
     /// contains a histogram with the name @a calHistName with an
     /// appended "_IMP" This histogram will be used to add an optional
     /// centrality percentile based on the generated impact
     /// parameter. If @increasing is true, a low (high) value of @proj
     /// is assumed to correspond to a more peripheral (central) event.
     const CentralityProjection&
     declareCentrality(const SingleValueProjection &proj,
                       string calAnaName, string calHistName,
                       const string projName, bool increasing = false);
 
     /// @brief Book a Pecentile wrapper around AnalysisObjects.
     ///
     /// Based on a previously registered CentralityProjection named @a
     /// projName book one AnalysisObject for each @a centralityBin and
     /// name them according to the corresponding code in the @a ref
     /// vector.
     template <class T>
     Percentile<T> bookPercentile(string projName,
                                  vector<pair<float, float> > centralityBins,
                                  vector<tuple<int, int, int> > ref) {
       typedef typename ReferenceTraits<T>::RefT RefT;
       Percentile<T> pctl(this, projName);
 
       const int nCent = centralityBins.size();
       for (int iCent = 0; iCent < nCent; ++iCent) {
         const string axisCode = makeAxisCode(std::get<0>(ref[iCent]),
                                              std::get<1>(ref[iCent]),
                                              std::get<2>(ref[iCent])); 
         shared_ptr<T> ao;
         CounterPtr cnt;
         try {
           ao = getAnalysisObject<T>(axisCode);
           MSG_TRACE("Found old " << histoPath(axisCode));
         }
         catch (Exception) {
           const RefT & refscatter = refData<RefT>(axisCode);
           ao = make_shared<T>(refscatter, histoPath(axisCode));
           addAnalysisObject(ao);
           MSG_TRACE("Created new " << histoPath(axisCode));
         }
         try {
           cnt = getAnalysisObject<Counter>("TMP/COUNTER/" + axisCode);
           MSG_TRACE("Found old " << histoPath("TMP/COUNTER/" + axisCode));
         }
         catch (Exception) {
           cnt = make_shared<Counter>(histoPath("TMP/COUNTER/" + axisCode));
           addAnalysisObject(cnt);
           MSG_TRACE("Created new " << histoPath("TMP/COUNTER/" + axisCode));
         }
         pctl.add(ao, cnt, centralityBins[iCent]);
       }
       return pctl;
     }
 
     /// @brief Book Pecentile wrappers around AnalysisObjects.
     ///
     /// Based on a previously registered CentralityProjection named @a
     /// projName book one (or several) AnalysisObject(s) named
     /// according to @a ref where the x-axis will be filled according
     /// to the percentile output(s) of the @projName.
     template <class T>
     PercentileXaxis<T> bookPercentileXaxis(string projName,
                                            tuple<int, int, int> ref) {
       typedef typename ReferenceTraits<T>::RefT RefT;
       PercentileXaxis<T> pctl(this, projName);
 
       const string axisCode = makeAxisCode(std::get<0>(ref),
                                            std::get<1>(ref),
                                            std::get<2>(ref)); 
       shared_ptr<T> ao;
       CounterPtr cnt;
       try {
         ao = getAnalysisObject<T>(histoPath(axisCode));
       }
       catch (Exception) {
 	const RefT & refscatter = refData<RefT>(axisCode);
         ao = make_shared<T>(refscatter, axisCode);
         addAnalysisObject(ao);
       }
       pctl.add(proj, ao, make_shared<Counter>());
 
       return pctl;
     }
 
 
 
     /// @name Analysis object manipulation
     /// @todo Should really be protected: only public to keep BinnedHistogram happy for now...
     //@{
 
     /// Multiplicatively scale the given counter, @a cnt, by factor @s factor.
     void scale(CounterPtr cnt, double factor);
 
     /// Multiplicatively scale the given counters, @a cnts, by factor @s factor.
     /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh)
     /// @todo Use SFINAE for a generic iterable of CounterPtrs
     void scale(const std::vector<CounterPtr>& cnts, double factor) {
       for (auto& c : cnts) scale(c, factor);
     }
     /// @todo YUCK!
     template <std::size_t array_size>
     void scale(const CounterPtr (&cnts)[array_size], double factor) {
       // for (size_t i = 0; i < std::extent<decltype(cnts)>::value; ++i) scale(cnts[i], factor);
       for (auto& c : cnts) scale(c, factor);
     }
 
 
     /// Normalize the given histogram, @a histo, to area = @a norm.
     void normalize(Histo1DPtr histo, double norm=1.0, bool includeoverflows=true);
 
     /// Normalize the given histograms, @a histos, to area = @a norm.
     /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh)
     /// @todo Use SFINAE for a generic iterable of Histo1DPtrs
     void normalize(const std::vector<Histo1DPtr>& histos, double norm=1.0, bool includeoverflows=true) {
       for (auto& h : histos) normalize(h, norm, includeoverflows);
     }
     /// @todo YUCK!
     template <std::size_t array_size>
     void normalize(const Histo1DPtr (&histos)[array_size], double norm=1.0, bool includeoverflows=true) {
       for (auto& h : histos) normalize(h, norm, includeoverflows);
     }
 
     /// Multiplicatively scale the given histogram, @a histo, by factor @s factor.
     void scale(Histo1DPtr histo, double factor);
 
     /// Multiplicatively scale the given histograms, @a histos, by factor @s factor.
     /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh)
     /// @todo Use SFINAE for a generic iterable of Histo1DPtrs
     void scale(const std::vector<Histo1DPtr>& histos, double factor) {
       for (auto& h : histos) scale(h, factor);
     }
     /// @todo YUCK!
     template <std::size_t array_size>
     void scale(const Histo1DPtr (&histos)[array_size], double factor) {
       for (auto& h : histos) scale(h, factor);
     }
 
 
     /// Normalize the given histogram, @a histo, to area = @a norm.
     void normalize(Histo2DPtr histo, double norm=1.0, bool includeoverflows=true);
 
     /// Normalize the given histograms, @a histos, to area = @a norm.
     /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh)
     /// @todo Use SFINAE for a generic iterable of Histo2DPtrs
     void normalize(const std::vector<Histo2DPtr>& histos, double norm=1.0, bool includeoverflows=true) {
       for (auto& h : histos) normalize(h, norm, includeoverflows);
     }
     /// @todo YUCK!
     template <std::size_t array_size>
     void normalize(const Histo2DPtr (&histos)[array_size], double norm=1.0, bool includeoverflows=true) {
       for (auto& h : histos) normalize(h, norm, includeoverflows);
     }
 
     /// Multiplicatively scale the given histogram, @a histo, by factor @s factor.
     void scale(Histo2DPtr histo, double factor);
 
     /// Multiplicatively scale the given histograms, @a histos, by factor @s factor.
     /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh)
     /// @todo Use SFINAE for a generic iterable of Histo2DPtrs
     void scale(const std::vector<Histo2DPtr>& histos, double factor) {
       for (auto& h : histos) scale(h, factor);
     }
     /// @todo YUCK!
     template <std::size_t array_size>
     void scale(const Histo2DPtr (&histos)[array_size], double factor) {
       for (auto& h : histos) scale(h, factor);
     }
 
 
     /// Helper for counter division.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void divide(CounterPtr c1, CounterPtr c2, Scatter1DPtr s) const;
 
     /// Helper for histogram division with raw YODA objects.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void divide(const YODA::Counter& c1, const YODA::Counter& c2, Scatter1DPtr s) const;
 
 
     /// Helper for histogram division.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void divide(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const;
 
     /// Helper for histogram division with raw YODA objects.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void divide(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const;
 
 
     /// Helper for profile histogram division.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void divide(Profile1DPtr p1, Profile1DPtr p2, Scatter2DPtr s) const;
 
     /// Helper for profile histogram division with raw YODA objects.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void divide(const YODA::Profile1D& p1, const YODA::Profile1D& p2, Scatter2DPtr s) const;
 
 
     /// Helper for 2D histogram division.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void divide(Histo2DPtr h1, Histo2DPtr h2, Scatter3DPtr s) const;
 
     /// Helper for 2D histogram division with raw YODA objects.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void divide(const YODA::Histo2D& h1, const YODA::Histo2D& h2, Scatter3DPtr s) const;
 
 
     /// Helper for 2D profile histogram division.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void divide(Profile2DPtr p1, Profile2DPtr p2, Scatter3DPtr s) const;
 
     /// Helper for 2D profile histogram division with raw YODA objects
     ///
     /// @note Assigns to the (already registered) output scatter, @a s.  Preserves the path information of the target.
     void divide(const YODA::Profile2D& p1, const YODA::Profile2D& p2, Scatter3DPtr s) const;
 
 
     /// Helper for histogram efficiency calculation.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void efficiency(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const;
 
     /// Helper for histogram efficiency calculation.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void efficiency(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const;
 
 
     /// Helper for histogram asymmetry calculation.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void asymm(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const;
 
     /// Helper for histogram asymmetry calculation.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void asymm(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const;
 
 
     /// Helper for converting a differential histo to an integral one.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void integrate(Histo1DPtr h, Scatter2DPtr s) const;
 
     /// Helper for converting a differential histo to an integral one.
     ///
     /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target.
     void integrate(const Histo1D& h, Scatter2DPtr s) const;
 
     //@}
 
 
   public:
 
     /// List of registered analysis data objects
     const vector<AnalysisObjectPtr>& analysisObjects() const {
       return _analysisobjects;
     }
 
 
   protected:
 
     /// @name Data object registration, retrieval, and removal
     //@{
 
     /// Register a data object in the histogram system
     void addAnalysisObject(AnalysisObjectPtr ao);
 
     /// Get a data object from the histogram system
     template <typename AO=YODA::AnalysisObject>
     const std::shared_ptr<AO> getAnalysisObject(const std::string& name) const {
       foreach (const AnalysisObjectPtr& ao, analysisObjects()) {
         if (ao->path() == histoPath(name)) return dynamic_pointer_cast<AO>(ao);
       }
       throw 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/Makefile.am b/include/Rivet/Makefile.am
--- a/include/Rivet/Makefile.am
+++ b/include/Rivet/Makefile.am
@@ -1,185 +1,186 @@
 ## Internal headers - not to be installed
 nobase_dist_noinst_HEADERS =
 ## Public headers - to be installed
 nobase_pkginclude_HEADERS =
 
 
 ## Rivet interface
 nobase_pkginclude_HEADERS += \
   Rivet.hh \
   Run.hh \
   Event.hh \
   ParticleBase.hh \
   Particle.fhh Particle.hh \
   Jet.fhh Jet.hh \
   Projection.fhh Projection.hh \
   ProjectionApplier.hh \
   ProjectionHandler.hh \
   Analysis.hh \
   AnalysisHandler.hh \
   AnalysisInfo.hh \
   AnalysisBuilder.hh \
   AnalysisLoader.hh
 
 
 ## Build config stuff
 nobase_pkginclude_HEADERS += \
   Config/RivetCommon.hh \
   Config/RivetConfig.hh \
   Config/BuildOptions.hh
 
 
 ## Projections
 nobase_pkginclude_HEADERS += \
   Projections/AliceCommon.hh \
   Projections/AxesDefinition.hh \
   Projections/Beam.hh \
   Projections/BeamThrust.hh \
   Projections/CentralEtHCM.hh \
   Projections/CentralityProjection.hh \
   Projections/ChargedFinalState.hh \
   Projections/ChargedLeptons.hh \
   Projections/ConstLossyFinalState.hh \
   Projections/DirectFinalState.hh \
   Projections/DISFinalState.hh \
   Projections/DISKinematics.hh \
   Projections/DISLepton.hh \
   Projections/DressedLeptons.hh \
   Projections/FastJets.hh \
   Projections/FinalPartons.hh \
   Projections/FinalState.hh \
   Projections/FoxWolframMoments.hh \
   Projections/FParameter.hh \
   Projections/GeneratedPercentileProjection.hh \
   Projections/HadronicFinalState.hh \
   Projections/HeavyHadrons.hh \
   Projections/Hemispheres.hh \
   Projections/IdentifiedFinalState.hh \
   Projections/ImpactParameterProjection.hh \
   Projections/IndirectFinalState.hh \
   Projections/InitialQuarks.hh \
   Projections/InvMassFinalState.hh \
   Projections/JetAlg.hh \
   Projections/JetShape.hh \
   Projections/LeadingParticlesFinalState.hh \
   Projections/LossyFinalState.hh \
   Projections/MergedFinalState.hh \
   Projections/MissingMomentum.hh \
   Projections/MixedFinalState.hh \
   Projections/NeutralFinalState.hh \
   Projections/NonHadronicFinalState.hh \
   Projections/NonPromptFinalState.hh \
   Projections/ParisiTensor.hh \
   Projections/ParticleFinder.hh \
   Projections/PartonicTops.hh \
   Projections/PercentileProjection.hh \
   Projections/PrimaryHadrons.hh \
   Projections/PrimaryParticles.hh \
   Projections/PromptFinalState.hh \
   Projections/SingleValueProjection.hh \
   Projections/SmearedParticles.hh \
   Projections/SmearedJets.hh \
   Projections/SmearedMET.hh \
   Projections/Sphericity.hh \
   Projections/Spherocity.hh \
   Projections/TauFinder.hh \
   Projections/Thrust.hh \
   Projections/TriggerCDFRun0Run1.hh \
   Projections/TriggerCDFRun2.hh \
   Projections/TriggerProjection.hh \
   Projections/TriggerUA5.hh \
   Projections/UnstableFinalState.hh \
   Projections/UserCentEstimate.hh \
   Projections/VetoedFinalState.hh \
   Projections/VisibleFinalState.hh \
   Projections/WFinder.hh \
   Projections/ZFinder.hh
 ## Meta-projection convenience headers
 nobase_pkginclude_HEADERS +=   \
   Projections/FinalStates.hh \
   Projections/Smearing.hh
 
 
 ## Analysis base class headers
 # TODO: Move to Rivet/AnalysisTools header dir?
 nobase_pkginclude_HEADERS += \
   Analyses/MC_Cent_pPb.hh \
   Analyses/MC_ParticleAnalysis.hh \
   Analyses/MC_JetAnalysis.hh \
   Analyses/MC_JetSplittings.hh
 
 
 ## Tools
 nobase_pkginclude_HEADERS += \
   Tools/AliceCommon.hh \
   Tools/AtlasCommon.hh \
   Tools/BeamConstraint.hh \
   Tools/BinnedHistogram.hh \
   Tools/CentralityBinner.hh \
   Tools/Cmp.fhh \
   Tools/Cmp.hh \
+  Tools/Correlators.hh \
   Tools/Cutflow.hh \
   Tools/Cuts.fhh \
   Tools/Cuts.hh \
   Tools/Exceptions.hh \
   Tools/JetUtils.hh \
   Tools/Logging.hh  \
   Tools/Random.hh  \
   Tools/ParticleBaseUtils.hh \
   Tools/ParticleIdUtils.hh \
   Tools/ParticleUtils.hh \
   Tools/ParticleName.hh \
   Tools/Percentile.hh \
   Tools/PrettyPrint.hh \
   Tools/RivetPaths.hh \
   Tools/RivetSTL.hh \
   Tools/RivetFastJet.hh \
   Tools/RivetHepMC.hh \
   Tools/RivetYODA.hh \
   Tools/RivetMT2.hh \
   Tools/SmearingFunctions.hh \
   Tools/MomentumSmearingFunctions.hh \
   Tools/ParticleSmearingFunctions.hh \
   Tools/JetSmearingFunctions.hh \
   Tools/TypeTraits.hh \
   Tools/Utils.hh \
   Tools/Nsubjettiness/AxesDefinition.hh \
   Tools/Nsubjettiness/ExtraRecombiners.hh \
   Tools/Nsubjettiness/MeasureDefinition.hh \
   Tools/Nsubjettiness/Njettiness.hh \
   Tools/Nsubjettiness/NjettinessPlugin.hh \
   Tools/Nsubjettiness/Nsubjettiness.hh \
   Tools/Nsubjettiness/TauComponents.hh \
   Tools/Nsubjettiness/XConePlugin.hh
 nobase_dist_noinst_HEADERS += \
   Tools/osdir.hh
 
 
 ## Maths
 nobase_pkginclude_HEADERS += \
   Math/Matrices.hh \
   Math/Vector3.hh \
   Math/VectorN.hh \
   Math/MatrixN.hh \
   Math/MatrixDiag.hh \
   Math/MathHeader.hh \
   Math/Vectors.hh \
   Math/LorentzTrans.hh \
   Math/Matrix3.hh \
   Math/MathUtils.hh \
   Math/Vector4.hh \
   Math/Math.hh \
   Math/Units.hh \
   Math/Constants.hh \
   Math/eigen/util.h \
   Math/eigen/regressioninternal.h \
   Math/eigen/regression.h \
   Math/eigen/vector.h \
   Math/eigen/ludecompositionbase.h \
   Math/eigen/ludecomposition.h \
   Math/eigen/linearsolver.h \
   Math/eigen/linearsolverbase.h \
   Math/eigen/matrix.h \
   Math/eigen/vectorbase.h \
   Math/eigen/projective.h \
   Math/eigen/matrixbase.h
diff --git a/include/Rivet/Tools/Correlators.hh b/include/Rivet/Tools/Correlators.hh
new file mode 100644
--- /dev/null
+++ b/include/Rivet/Tools/Correlators.hh
@@ -0,0 +1,1084 @@
+// -*- C++ -*-
+#ifndef RIVET_Correlators_HH
+#define RIVET_Correlators_HH
+
+#include "Rivet/Projection.hh"
+#include "Rivet/Projections/ParticleFinder.hh"
+#include <complex>
+#include "YODA/Scatter2D.h"
+#include "Rivet/Analysis.hh"
+/* File contains tools for calculating flow coefficients using 
+ * correlators.
+ * Classes:
+ * Correlators: Calculates single event correlators of a given harmonic.
+ * Cumulants: An additional base class for flow analyses 
+ * (use as: class MyAnalysis : public Analysis, Cumulants {};)
+ * Includes a framework for calculating cumulants and flow coefficients
+ * from single event correlators, including automatic handling of statistical
+ * errors. Contains multiple internal sub-classes:
+ * CorBinBase: Base class for correlators binned in event or particle observables.
+ * CorSingleBin: A simple bin for correlators.
+ * CorBin: Has the interface of a simple bin, but does automatic calculation
+ * of statistical errors by a bootstrap method.
+ * ECorrelator: Data type for event averaged correlators. 
+ * @author Vytautas Vislavicius, Christine O. Rasmussen, Christian Bierlich.
+ */
+
+namespace Rivet {
+
+  /* @brief Projection for calculating correlators for flow measurements.
+  *  
+  *   A projection which calculates Q-vectors and P-vectors, and projects
+  *   them out as correlators. Implementation follows the description of 
+  *   the ''Generic Framework'' 
+  *   Phys. Rev. C 83 (2011) 044913, arXiv: 1010.0233
+  *   Phys. Rev. C 89 (2014) 064904, arXiv: 1312.3572
+  *   
+  */
+
+  class Correlators : public Projection {
+  
+  public:
+
+    // Constructor. Takes parameters @parm fsp, the Final State
+    // projection correlators should be constructed from, @parm nMaxIn, 
+    // the maximal sum of harmonics, eg. for 
+    // c_2{2} = {2,-2} = 2 + 2 = 4
+    // c_2{4} = {2,2,-2,-2} = 2 + 2 + 2 + 2 = 8
+    // c_4{2} = {4,-4} = 4 + 4 = 8
+    // c_4{4} = {4,4,-4,-4} = 4 + 4 + 4 + 4 = 16. 
+    // @parm pMaxIn is the maximal number of particles
+    // you want to correlate and @parm pTbinEdgesIn is the (lower) 
+    // edges of pT bins, the last one the upper edge of the final bin.
+    Correlators(const ParticleFinder& fsp, int nMaxIn = 2, 
+      int pMaxIn = 0, vector<double> pTbinEdgesIn = {});
+    
+    // Constructor which takes a Scatter2D to estimate bin edges.
+    Correlators(const ParticleFinder& fsp, int nMaxIn, 
+      int pMaxIn, const Scatter2DPtr hIn);
+    
+    /// @brief Integrated correlator of @parm n harmonic, with the
+    /// number of powers being the size of @parm n.
+    /// Eg. @parm n should be:
+    /// <<2>>_2 => n = {2, -2}
+    /// <<4>>_2 => n = {2, 2, -2, -2}
+    /// <<2>>_4 => n = {4, -4}
+    /// <<4>>_4 => n = {4, 4, -4, 4} and so on.
+    const pair<double,double> intCorrelator(vector<int> n) const;
+
+    /// @brief pT differential correlator of @parm n harmonic, with the
+    /// number of powers being the size of @parm n.
+    /// The method can include overflow/underflow bins in the 
+    /// beginning/end of the returned vector, by toggling
+    /// @parm overflow = true.
+    const vector<pair<double,double>> pTBinnedCorrelators(vector<int> n,
+      bool overflow = false) const;
+
+    /// @brief Integrated correlator of @parm n1 harmonic, with the
+    /// number of powers being the size of @parm n1. This method
+    /// imposes an eta gap, correlating with another phase space, 
+    /// where another Correlators projection @parm other should be 
+    /// defined. The harmonics of the other phase space is given 
+    /// as @parm n2.
+    /// To get eg. integrated <<4>>_2, @parm n1 should be:
+    /// n1 = {2, 2} and n2 = {-2, -2}
+    const pair<double,double> intCorrelatorGap(const Correlators& other, 
+      vector<int> n1, vector<int> n2) const;
+
+    /// @brief pT differential correlators of @parm n1 harmonic, with the
+    /// number of powers being the size of @parm n1. This method
+    /// imposes an eta gap, correlating with another phase space, 
+    /// where another Correlators projection @parm other should be 
+    /// defined. The harmonics of the other phase space is given 
+    /// as @parm n2.
+    /// To get eg. differential <<4'>_2, @parm n1 should be:
+    /// n1 = {2, 2} and @parm n2: n2 = {-2, -2}.
+    /// To get eg. differential <<2'>>_4, @parm n1 should be:
+    /// n1 = {4} and @parm n2: n2 = {-4}.
+    /// The method can include overflow/underflow
+    /// bins in the beginning/end of the returned vector, by toggling
+    /// @parm overflow = true.
+    const vector<pair<double,double>> pTBinnedCorrelatorsGap(
+      const Correlators& other, vector<int> n1, vector<int> n2, 
+      bool oveflow = false) const; 
+    
+    /// @brief Construct a harmonic vectors from @parm n harmonics
+    /// and @parm m number of particles.
+    /// TODO: In C++14 this can be done much nicer with TMP.
+    static vector<int> hVec(int n, int m) {
+      if (m%2 != 0) {
+        cout << "Harmonic Vector: Number of particles "
+	        "must be an even number." << endl;
+        return {};
+      }
+      vector<int> ret = {};
+      for (int i = 0; i < m; ++i) {
+        if (i < m/2) ret.push_back(n);
+        else ret.push_back(-n);
+      }
+      return ret;
+    }
+
+    /// @brief Return the maximal values for n, p to be used in the
+    /// constructor of Correlators(xxx, nMax, pMax, xxxx)
+    static pair<int,int> getMaxValues(vector< vector<int> >& hList){
+      int nMax = 0, pMax = 0;
+      for (vector<int> h : hList) {
+        int tmpN = 0, tmpP = 0;
+        for ( int i = 0; i < int(h.size()); ++i) {
+          tmpN += abs(h[i]);
+          ++tmpP;
+        }
+        if (tmpN > nMax) nMax = tmpN;
+        if (tmpP > pMax) pMax = tmpP;
+      }
+      return make_pair(nMax,pMax);
+    }
+    // Clone on the heap.
+    DEFAULT_RIVET_PROJ_CLONE(Correlators);
+ 
+  protected:
+    
+    // @brief Project function. Loops over array and calculates Q vectors
+    // and P vectors if needed 
+    void project(const Event& e); 
+    
+    // @brief Compare against other projection. Test if harmonics, pT-bins
+    // and underlying final state are similar.
+    int compare(const Projection& p) const {
+      const Correlators* other = dynamic_cast<const Correlators*>(&p);
+      if (nMax != other->nMax) return UNDEFINED;
+      if (pMax != other->pMax) return UNDEFINED;
+      if (pTbinEdges != other->pTbinEdges) return UNDEFINED;
+      return mkPCmp(p, "FS");
+    }; 
+
+    // @brief Calculate correlators from one particle. 
+    void fillCorrelators(const Particle& p, const double& weight);   
+    
+    // @brief Return a Q-vector.
+    const complex<double> getQ(int n, int p) const {
+      bool isNeg = (n < 0);
+      if (isNeg) return conj( qVec[abs(n)][p] );
+      else       return qVec[n][p];
+    };
+    
+    // @brief Return a P-vector. 
+    const complex<double> getP(int n, int p, double pT = 0.) const {
+      bool isNeg = (n < 0);
+      map<double, Vec2D>::const_iterator pTitr = pVec.lower_bound(pT);
+      if (pTitr == pVec.end()) return DBL_NAN;
+      if (isNeg) return conj( pTitr->second[abs(n)][p] );
+      else       return pTitr->second[n][p];
+    };
+
+  private:
+    // Find correlators by recursion. Order = M (# of particles), 
+    // n's are harmonics, p's are the powers of the weights
+    const complex<double> recCorr(int order, vector<int> n,  
+      vector<int> p, bool useP, double pT = 0.) const;  
+
+    // Two-particle correlator Eq. (19) p. 6
+    // Flag if p-vectors or q-vectors should be used to
+    // calculate the correlator.
+    const complex<double> twoPartCorr(int n1, int n2, int p1 = 1, 
+      int p2 = 1, double pT = 0., bool useP = false) const;
+
+    // Set elements in vectors to zero.
+    void setToZero();
+    
+    // Shorthands for setting and comparing to zero.
+    const complex<double> _ZERO = {0., 0.};
+    const double _TINY = 1e-10;
+
+    // Shorthand typedefs for vec<vec<complex>>.
+    typedef vector< vector<complex<double>> > Vec2D;
+
+    // Define Q-vectors and p-vectors
+    Vec2D qVec; // Q[n][p] 
+    map<double, Vec2D> pVec; // p[pT][n][p]
+   
+    // The max values of n and p to be calculated.
+    int nMax, pMax;
+    // pT bin-edges
+    vector<double> pTbinEdges;
+    bool isPtDiff;
+ 
+  /// End class Correlators
+  };
+
+
+  /// @brief Tools for flow analyses.
+  /// The following are helper classes to construct event averaged correlators
+  /// as well as cummulants and flow coefficents from the basic event 
+  //  correlators defined above. They are all encapsulated in a Cumulants 
+  //  class, which can be used as a(nother) base class for flow analyses, 
+  //  to ensure access.
+ 
+  class CumulantAnalysis : public Analysis {
+  private:
+
+  // Number of bins used for bootstrap calculation of statistical 
+  // uncertainties. It is hard coded, and shout NOT be changed unless there
+  // are good reasons to do so. 
+  static const int BOOT_BINS = 9;
+
+  // Enum for choosing the method of error calculation.
+  enum Error {
+    VARIANCE,
+    ENVELOPE
+  };
+  
+  // The desired error method. Can be changed in the analysis constructor
+  // by setting it appropriately. 
+  Error errorMethod;
+
+  /// @brief Base class for correlator bins.
+  class CorBinBase {
+  public:
+    // Derived class should have fill and mean defined.
+    virtual void fill(const pair<double, double>& cor, const double& weight) = 0;
+    virtual const double mean() const = 0;
+  };
+
+  /// @brief The CorSingleBin is the basic quantity filled in an ECorrelator.
+  /// It is a simple counter with an even simpler structure than normal
+  /// YODA type DBNs, but added functionality to test for out of
+  /// bounds correlators.
+  class CorSingleBin : public CorBinBase {
+  
+  public:
+    /// @brief The default constructor.
+    CorSingleBin() : _sumWX(0.), _sumW(0.), _sumW2(0.), _numEntries(0.) {}
+
+    /// @brief Fill a correlator bin with the return type from a 
+    /// Correlator (a pair giving numerator and denominator of <M>_event).
+    void fill(const pair<double, double>& cor, const double& weight) {
+      _numEntries += 1.;
+      // Test if denominator for the single event average is zero.
+      if (cor.second < 1e-10) return;
+      // The single event average <M> is then cor.first / cor.second.
+      // With weights this becomes just:
+      _sumWX += cor.first * weight;
+      _sumW += weight * cor.second;
+      _sumW2 += weight * weight * cor.second * cor.second;
+    }
+
+    const double mean() const {
+      if (_sumW < 1e-10) return 0;
+      return _sumWX / _sumW;
+    }
+
+    // @brief Sum of weights.
+    const double sumW() const {
+      return _sumW;
+    }
+
+    const double sumW2() const {
+      return _sumW2;
+    }
+
+    // @brief Sum of weight * X.
+    const double sumWX() const {
+      return _sumWX;
+    }
+
+    // @brief Number of entries.
+    const double numEntries() const {
+      return _numEntries;
+    }
+    
+  private:
+    double _sumWX, _sumW, _sumW2, _numEntries;
+  
+  }; // End of CorSingleBin sub-class.
+
+  /// @brief The CorBin is the basic bin quantity in ECorrelators.
+  /// It consists of several CorSingleBins, to facilitate 
+  /// bootstrapping calculation of statistical uncertainties.
+  class CorBin : public CorBinBase {
+  public:
+    // @brief The constructor. nBins signifies the period of the bootstrap
+    // calculation, and should never be changed here, but in its definition
+    // above -- and only if there are good reasons to do so.
+    CorBin() : binIndex(0), nBins(BOOT_BINS) {
+      for(size_t i = 0; i < nBins; ++i) bins.push_back(CorSingleBin()); 
+    }
+
+    // @brief Fill the correct underlying bin and take a step.
+    void fill(const pair<double, double>& cor, const double& weight) {
+      // Test if denominator for the single event average is zero.
+      if (cor.second < 1e-10) return;
+      // Fill the correct bin.
+      bins[binIndex].fill(cor, weight);
+      if (binIndex == nBins - 1) binIndex = 0;
+      else ++binIndex;
+    }
+
+    // @brief Calculate the total sample mean with all
+    // available statistics.
+    const double mean() const {
+      double sow = 0;
+      double sowx = 0;
+      for(auto b : bins) {
+	if (b.sumW() < 1e-10) continue;
+	sow += b.sumW();
+	sowx += b.sumWX();
+      }
+      return sowx / sow;
+    }
+
+    // @brief Return a copy of the bins.
+    const vector<CorSingleBin> getBins() const {
+      return bins;
+    }
+
+    // @brief Return the bins as pointers to the base class.
+    template<class T=CorBinBase>
+    const vector<T*> getBinPtrs() {
+      vector<T*> ret(bins.size());
+      transform(bins.begin(), bins.end(), ret.begin(),
+        [](CorSingleBin& b) {return &b;});
+      return ret;
+    }
+
+  private:
+    vector<CorSingleBin> bins;
+    size_t binIndex;
+    size_t nBins; 
+
+  }; // End of CorBin sub-class.
+  
+  public:
+  /// @brief The ECorrelator is a helper class to calculate all event 
+  /// averages of correlators, in order to construct cumulants.
+  /// It can be binned in any variable.
+  class ECorrelator {
+  
+  public:
+    /// @brief Constructor. Takes as argument the desired harmonic and number
+    /// of correlated particles as a generic framework style vector, eg,
+    /// {2, -2} for <<2>>_2, no binning.
+    /// TODO: Implement functionality for this if needed. 
+    //ECorrelator(vector<int> h)  : h1(h), h2({}),
+    //  binX(0), binContent(0), reference() {    
+    //};
+
+    /// @brief Constructor. Takes as argument the desired harmonic and number
+    /// of correlated particles as a generic framework style vector, eg,
+    /// {2, -2} for <<2>>_2 and binning.
+    ECorrelator(vector<int> h, vector<double> binIn) : h1(h), h2({}),
+      binX(binIn), binContent(binIn.size() - 1), reference() {};
+
+    /// @brief Constructor for gapped correlator. Takes as argument the 
+    /// desired harmonics for the two final states, and binning. 
+    ECorrelator(vector<int> h1In, vector<int> h2In, vector<double> binIn) : 
+      h1(h1In), h2(h2In), binX(binIn), binContent(binIn.size() - 1), 
+      reference() {};
+    
+    /// @brief Fill the appropriate bin given an input (per event)
+    /// observable, eg. centrality.
+    void fill(const double& obs, const Correlators& c, 
+      const double weight = 1.0) {
+      int index = getBinIndex(obs);
+      if (index < 0) return;
+      binContent[index].fill(c.intCorrelator(h1), weight);
+    }
+
+    /// @brief Fill the appropriate bin given an input (per event)
+    /// observable, eg. centrality. Using a rapidity gap between
+    /// two Correlators.
+    void fill(const double& obs, const Correlators& c1,
+      const Correlators& c2, const double weight = 1.0) {
+      if (!h2.size()) {
+        cout << "Trying to fill gapped correlator, but harmonics behind "
+        "the gap (h2) are not given!" << endl;
+      	return;
+      }
+      int index = getBinIndex(obs);
+      if (index < 0) return;
+      binContent[index].fill(c1.intCorrelatorGap(c2, h1, h2), weight);
+    }
+
+    /// @brief Fill the bins with the appropriate correlator, taking the 
+    /// binning directly from the Correlators object, and filling also the 
+    /// reference flow.
+    void fill(const Correlators& c, const double& weight = 1.0) {
+      vector< pair<double, double> > diffCorr = c.pTBinnedCorrelators(h1);
+      // We always skip overflow when calculating the all event average.
+      if (diffCorr.size() != binX.size() - 1)
+        cout << "Tried to fill event with wrong binning (ungapped)" << endl;
+      for (size_t i = 0; i < diffCorr.size(); ++i) {
+        int index = getBinIndex(binX[i]);
+        if (index < 0) return;
+        binContent[index].fill(diffCorr[i], weight);
+      }
+      reference.fill(c.intCorrelator(h1), weight);
+    }
+
+    /// @brief Fill bins with the appropriate correlator, taking the binning
+    /// directly from the Correlators object, and also the reference flow.
+    /// Using a rapidity gap between two Correlators.
+    void fill(const Correlators& c1, const Correlators& c2, 
+      const double& weight = 1.0) {
+      if (!h2.size()) {
+        cout << "Trying to fill gapped correlator, but harmonics behind "
+        "the gap (h2) are not given!" << endl;
+        return;
+      }
+      vector< pair<double, double> > diffCorr = c1.pTBinnedCorrelatorsGap(c2, h1, h2);
+      // We always skip overflow when calculating the all event average.
+      if (diffCorr.size() != binX.size() - 1)
+        cout << "Tried to fill event with wrong binning (gapped)" << endl;
+      for (size_t i = 0; i < diffCorr.size(); ++i) {
+	int index = getBinIndex(binX[i]);
+	if (index < 0) return;
+        binContent[index].fill(diffCorr[i], weight);
+      }
+      reference.fill(c1.intCorrelatorGap(c2, h1, h2), weight);
+    }
+    
+    /// @brief Get a copy of the bin contents.
+    const vector<CorBin> getBins() const {
+      return binContent;
+    }
+
+    // @brief Return the bins as pointers to the base class.
+    const vector<CorBinBase*> getBinPtrs() {
+      vector<CorBinBase*> ret(binContent.size());
+      transform(binContent.begin(), binContent.end(), ret.begin(),
+        [](CorBin& b) {return &b;});
+      return ret;
+    }
+
+    /// @brief Get a copy of the bin x-values.
+    const vector<double> getBinX() const {
+      return binX;
+    }
+
+    /// @brief Get a copy of the @ret h1 harmonic vector.
+    const vector<int> getH1() const {
+      return h1;
+    }
+
+    /// @brief Get a copy of the @ret h2 harmonic vector.
+    const vector<int> getH2() const {
+      return h2;
+    }
+
+    /// @brief Replace reference flow bin with another reference
+    /// flow bin, eg. calculated in another phase space or with
+    /// other pid.
+    void setReference(CorBin refIn) {
+      reference = refIn;
+    }
+
+    /// @brief Extract the reference flow from a differential event
+    /// averaged correlator.
+    const CorBin getReference() const {
+      if (reference.mean() < 1e-10)
+        cout << "Warning: ECorrelator, reference bin is zero." << endl;
+      return reference;
+    }
+    /// @brief set the @parm prIn list of profile histograms associated
+    /// with the internal bins. Is called automatically when booking, no
+    /// need to call it yourself.
+    void setProfs(list<Profile1DPtr> prIn) {
+      profs = prIn;
+      cout << "TEST: " << (*prIn.begin())->effNumEntries() << endl;
+    }
+
+    /// @brief begin() iterator for the list of associated profile histograms.
+    list<Profile1DPtr>::iterator profBegin() {
+      return profs.begin();
+    } 
+    
+    /// @brief end() iterator for the list of associated profile histograms.
+    list<Profile1DPtr>::iterator profEnd() {
+      return profs.end();
+    } 
+
+  private:
+    /// @brief Get correct bin index for a given @parm obs value.
+    const int getBinIndex(const double& obs) const {
+      // Find the correct index of binContent.
+      // If we are in overflow, just skip.
+      if (obs >= binX.back()) return -1;
+      // If we are in underflow, ditto.
+      if (obs < binX[0]) return -1;
+      int index = 0;
+      for (int i = 0, N = binX.size() - 1; i < N; ++i, ++index)
+        if (obs >= binX[i] && obs < binX[i + 1]) break;
+      return index;
+    }
+
+    // The harmonics vectors.
+    vector<int> h1;
+    vector<int> h2;
+    // The bins.
+    vector<double> binX;
+    vector<CorBin> binContent;
+    // The reference flow.
+    CorBin reference;
+    // The profile histograms associated with the CorBins for streaming.
+    list<Profile1DPtr> profs;
+
+  }; // End of ECorrelator sub-class.
+
+  // @brief Get the correct max N and max P for the set of 
+  // booked correlators.
+  const pair<int, int> getMaxValues() const {
+    vector< vector<int>> harmVecs;
+    for ( auto eItr = eCorrPtrs.begin(); eItr != eCorrPtrs.end(); ++eItr) {
+      vector<int> h1 = (*eItr)->getH1();
+      vector<int> h2 = (*eItr)->getH2();
+      if (h1.size() > 0) harmVecs.push_back(h1);
+      if (h2.size() > 0) harmVecs.push_back(h2);
+    }
+    if (harmVecs.size() == 0) {
+      cout << "Warning: You tried to extract max values from harmonic "
+	      "vectors, but have not booked any." << endl;
+      return pair<int,int>();
+    }
+    return Correlators::getMaxValues(harmVecs);
+  }
+
+  // Typedeffing shared pointer to ECorrelator.
+  typedef shared_ptr<ECorrelator> ECorrPtr;
+  
+  // @brief Book an ECorrelator in the same way as a histogram.
+  ECorrPtr bookECorrelator(const string name, const vector<int>& h, const Scatter2DPtr hIn) {
+    vector<double> binIn;
+    for (auto b : hIn->points()) binIn.push_back(b.xMin());
+    binIn.push_back(hIn->points().back().xMax());
+    ECorrPtr ecPtr = ECorrPtr(new ECorrelator(h, binIn));
+    list<Profile1DPtr> eCorrProfs;
+    for (int i = 0; i < BOOT_BINS; ++i)
+      eCorrProfs.push_back(bookProfile1D(name+"-e"+to_string(i),*hIn));
+    ecPtr->setProfs(eCorrProfs);
+    eCorrPtrs.push_back(ecPtr);
+    return ecPtr;
+  }
+
+  // @brief Book a gapped ECorrelator with two harmonic vectors.
+  ECorrPtr bookECorrelator(const string name, const vector<int>& h1, 
+    const vector<int>& h2, const Scatter2DPtr hIn ) {
+    vector<double> binIn;
+    for (auto b : hIn->points()) binIn.push_back(b.xMin());
+    binIn.push_back(hIn->points().back().xMax());
+    ECorrPtr ecPtr = ECorrPtr(new ECorrelator(h1, h2, binIn));
+    list<Profile1DPtr> eCorrProfs;
+    for (int i = 0; i < BOOT_BINS; ++i)
+      eCorrProfs.push_back(bookProfile1D(name+"-e"+to_string(i),*hIn));
+    ecPtr->setProfs(eCorrProfs);
+    eCorrPtrs.push_back(ecPtr);
+    return ecPtr;
+  }
+  
+  // @brief Short hand for gapped correlators which splits the harmonic 
+  // vector into negative and positive components.
+  ECorrPtr bookECorrelatorGap (const string name, const vector<int>& h, 
+    const Scatter2DPtr hIn) {
+    const vector<int> h1(h.begin(), h.begin() + h.size() / 2);
+    const vector<int> h2(h.begin() + h.size() / 2, h.end());
+    return bookECorrelator(name, h1, h2, hIn);
+  }
+
+  // @brief Templated version of correlator booking which takes
+  // @parm N desired harmonic and @parm M number of particles.
+  template<unsigned int N, unsigned int M>
+  ECorrPtr bookECorrelator(const string name, const Scatter2DPtr hIn) {
+    return bookECorrelator(name, Correlators::hVec(N, M), hIn);
+  }
+
+  // @brief Templated version of gapped correlator booking which takes
+  // @parm N desired harmonic and @parm M number of particles.
+  template<unsigned int N, unsigned int M>
+  ECorrPtr bookECorrelatorGap(const string name, const Scatter2DPtr hIn) {
+    return bookECorrelatorGap(name, Correlators::hVec(N, M), hIn);
+  }
+  
+  // @brief Finalize MUST be explicitly called for this base class with a 
+  // CumulantsAnalysis::finalize() call as the first thing in each analysis,
+  // in order to stream the contents of ECorrelators to a yoda file for 
+  // reentrant finalize.
+  void finalize() {
+    for (auto ecItr = eCorrPtrs.begin(); ecItr != eCorrPtrs.end(); ++ecItr)
+      corrPlot(list<Profile1DPtr>((*ecItr)->profBegin(),
+        (*ecItr)->profEnd()), *ecItr);
+  }
+  private:
+  
+  /// Bookkeeping of the event averaged correlators. 
+  list<ECorrPtr> eCorrPtrs;
+  
+  public:
+  
+    // @brief Constructor. Use CumulantAnalysis as base class for the
+    // analysis to have access to functionality.
+    CumulantAnalysis (string n) : Analysis(n), errorMethod(ENVELOPE) {};
+    // @brief Helper method for turning correlators into Scatter2Ds.
+    // Takes @parm h a pointer to the resulting Scatter2D, @parm binx
+    // the x-bins and a function @parm func defining the transformation.
+    // Makes no attempt at statistical errors.
+    // See usage in the methods below.
+    // Can also be used directly in the analysis if a user wants to 
+    // perform an unforseen transformation from correlators to 
+    // Scatter2D.
+    template<typename T>
+    static void fillScatter(Scatter2DPtr h, vector<double>& binx, T func) {
+      vector<YODA::Point2D> points;
+      for (int i = 0, N = binx.size() - 1; i < N; ++i) {
+        double xMid = (binx[i] + binx[i + 1]) / 2.0;
+        double xeMin = fabs(xMid - binx[i]);
+        double xeMax = fabs(xMid - binx[i + 1]);
+        double yVal = func(i);
+        if (isnan(yVal)) yVal = 0.;
+        double yErr = 0;
+        points.push_back(YODA::Point2D(xMid, yVal, xeMin, xeMax, yErr, yErr));
+        }
+       h->reset();
+       h->points().clear();
+       for (int i = 0, N = points.size(); i < N; ++i) h->addPoint(points[i]);
+    }
+
+    // @brief Helper method for turning correlators into Scatter2Ds
+    // with error estimates.
+    // Takes @parm h a pointer to the resulting Scatter2D, @parm binx
+    // the x-bins, a function @parm func defining the transformation 
+    // and a vector of errors @parm err.
+    // See usage in the methods below.
+    // Can also be used directly in the analysis if a user wants to 
+    // perform an unforseen transformation from correlators to 
+    // Scatter2D.
+    template<typename F>
+    const void fillScatter(Scatter2DPtr h, vector<double>& binx, F func,
+      vector<pair<double, double> >& yErr) const {
+      vector<YODA::Point2D> points;
+      for (int i = 0, N = binx.size() - 1; i < N; ++i) {
+        double xMid = (binx[i] + binx[i + 1]) / 2.0;
+        double xeMin = fabs(xMid - binx[i]);
+        double xeMax = fabs(xMid - binx[i + 1]);
+        double yVal = func(i);
+        if (isnan(yVal))
+          points.push_back(YODA::Point2D(xMid, 0., xeMin, xeMax,0., 0.));
+	else
+          points.push_back(YODA::Point2D(xMid, yVal, xeMin, xeMax,
+	    yErr[i].first, yErr[i].second));
+        }
+       h->reset();
+       h->points().clear();
+       for (int i = 0, N = points.size(); i < N; ++i) 
+         h->addPoint(points[i]);
+    }
+
+
+    /// @brief Take the @parm n th power of all points in @parm hIn,
+    /// and put the result in @parm hOut. Optionally put a  
+    /// @parm k constant below the root.
+    static void nthPow(Scatter2DPtr hOut, const Scatter2DPtr hIn,
+      const double& n, const double& k = 1.0) {
+      if (n == 0 || n == 1) {
+        cout << "Error: Do not take the 0th or 1st power of a Scatter2D,"
+	       " use scale instead." << endl;
+	return;
+      }
+      if (hIn->points().size() != hOut->points().size()) {
+        cout << "nthRoot: Scatterplots: " << hIn->name() << " and " <<
+          hOut->name() << " not initialized with same length." << endl;
+        return;
+      }
+      vector<YODA::Point2D> points;
+      // The error pre-factor is k^(1/n) / n by Taylors formula.
+      double eFac = pow(k,1./n)/n;
+      for (auto b : hIn->points()) {
+        double yVal = pow(k * b.y(),n);
+        if (isnan(yVal))
+          points.push_back(YODA::Point2D(b.x(), 0., b.xErrMinus(),
+	    b.xErrPlus(), 0, 0 ));
+	else {
+	  double yemin = abs(eFac * pow(yVal,1./(n - 1.))) * b.yErrMinus();
+	  if (isnan(yemin)) yemin = b.yErrMinus();
+	  double yemax = abs(eFac * pow(yVal,1./(n - 1.))) * b.yErrPlus();
+	  if (isnan(yemax)) yemax = b.yErrPlus();
+	  points.push_back(YODA::Point2D(b.x(), yVal, b.xErrMinus(), 
+	    b.xErrPlus(), yemin, yemax ));
+	}  
+      }
+      hOut->reset();
+      hOut->points().clear();
+      for (int i = 0, N = points.size(); i < N; ++i) 
+        hOut->addPoint(points[i]);
+    }
+    
+    /// @brief Take the @parm n th power of all points in @parm h,
+    /// and put the result back in the same Scatter2D. Optionally put a  
+    /// @parm k constant below the root.
+    static void nthPow(Scatter2DPtr h, const double& n, 
+      const double& k = 1.0) {
+      if (n == 0 || n == 1) {
+        cout << "Error: Do not take the 0th or 1st power of a Scatter2D,"
+	       " use scale instead." << endl;
+	return;
+      }
+      vector<YODA::Point2D> points;
+      vector<YODA::Point2D> pIn = h->points();
+      // The error pre-factor is k^(1/n) / n by Taylors formula.
+      double eFac = pow(k,1./n)/n;
+      for (auto b : pIn) {
+        double yVal =  pow(k * b.y(),n);
+        if (isnan(yVal))
+          points.push_back(YODA::Point2D(b.x(), 0., b.xErrMinus(), 
+	    b.xErrPlus(), 0, 0 ));
+	else {
+	  double yemin = abs(eFac * pow(yVal,1./(n - 1.))) * b.yErrMinus();
+	  if (isnan(yemin)) yemin = b.yErrMinus();
+	  double yemax = abs(eFac * pow(yVal,1./(n - 1.))) * b.yErrPlus();
+	  if (isnan(yemax)) yemax = b.yErrPlus();
+	  points.push_back(YODA::Point2D(b.x(), yVal, b.xErrMinus(), 
+	    b.xErrPlus(), yemin, yemax ));
+	}  
+      }
+      h->reset();
+      h->points().clear();
+      for (int i = 0, N = points.size(); i < N; ++i) h->addPoint(points[i]);
+    }
+   
+    // @brief Calculate the bootstrapped sample variance for the observable
+    // given by correlators and the transformation @parm func. 
+    template<typename T>
+    static pair<double, double> sampleVariance(T func) {
+      // First we calculate the mean (two pass calculation).
+      double avg = 0.;
+      for (int i = 0; i < BOOT_BINS; ++i) avg += func(i);
+      avg /= double(BOOT_BINS);
+      // Then we find the variance.
+      double var = 0.;
+      for (int i = 0; i < BOOT_BINS; ++i) var += pow(func(i) - avg, 2.);
+      var /= (double(BOOT_BINS) - 1);
+      return pair<double, double>(var, var);
+    }
+
+    // @brief Calculate the bootstrapped sample envelope for the observable
+    // given by correlators and the transformation @parm func. 
+    template<typename T>
+    static pair<double, double> sampleEnvelope(T func) {
+      // First we calculate the mean.
+      double avg = 0.;
+      for (int i = 0; i < BOOT_BINS; ++i) avg += func(i);
+      avg /= double(BOOT_BINS);
+      double yMax = avg;
+      double yMin = avg;
+      // Then we find the envelope using the mean as initial value.
+      for (int i = 0; i < BOOT_BINS; ++i) {
+        double yVal = func(i);
+        if (yMin > yVal) yMin = yVal;
+        else if (yMax < yVal) yMax = yVal;
+      }
+      return pair<double, double>(fabs(avg - yMin), fabs(yMax - avg));
+    }
+
+    // @brief Selection method for which sample error to use, given
+    // in the constructor.
+    template<typename T>
+    const pair<double, double> sampleError(T func) const {
+      if (errorMethod == VARIANCE) return sampleVariance(func);
+      else if (errorMethod == ENVELOPE) return sampleEnvelope(func);
+      else
+      cout << "Error: Error method not found!" << endl;
+      return pair<double, double>(0.,0.);
+    }    
+
+    // @brief Two particle integrated cn.
+    const void cnTwoInt(Scatter2DPtr h, ECorrPtr e2) const {
+      vector<CorBin> bins = e2->getBins();
+      vector<double> binx = e2->getBinX();
+      // Assert bin size.
+      if (binx.size() - 1 != bins.size()){
+        cout << "cnTwoInt: Bin size (x,y) differs!" << endl;
+        return;
+      }
+      vector<CorBinBase*> binPtrs;
+      // The mean value of the cumulant.
+      auto cn = [&] (int i) { return binPtrs[i]->mean(); };
+      // Error calculation.
+      vector<pair<double, double> > yErr;
+      for (int j = 0, N = bins.size(); j < N; ++j) {
+        binPtrs = bins[j].getBinPtrs();
+	yErr.push_back(sampleError(cn));
+      }
+      binPtrs = e2->getBinPtrs();
+      fillScatter(h, binx, cn, yErr);
+    }
+    
+    // @brief Two particle integrated vn.
+    const void vnTwoInt(Scatter2DPtr h, ECorrPtr e2) const {
+      cnTwoInt(h, e2);
+      nthPow(h, 0.5);
+    }
+
+    // @brief Put an event averaged correlator into a Scatter2D.
+    // Reduces to cnTwoInt, but better with a proper name.
+    const void corrPlot(Scatter2DPtr h, ECorrPtr e) const {
+      cnTwoInt(h, e);
+    }
+
+    // @brief Put an event averaged correlator into Profile1Ds, one
+    // for each bootstrapping bin.
+    const void corrPlot(list<Profile1DPtr> profs, ECorrPtr e) const {
+      vector<CorBin> corBins = e->getBins();
+      vector<double> binx = e->getBinX();
+      // Assert bin size.
+      if (binx.size() - 1 != corBins.size()){
+        cout << "corrPlot: Bin size (x,y) differs!" << endl;
+        return;
+      }
+      list<Profile1DPtr>::iterator hItr = profs.begin();
+      // Loop over the boostrapped correlators.
+      for (size_t i = 0; i < profs.size(); ++i, ++hItr) {
+        vector<YODA::ProfileBin1D> profBins;
+	for (size_t j = 0, N = binx.size() - 1; j < N; ++j) {
+	  vector<CorSingleBin*> binPtrs = 
+	    corBins[j].getBinPtrs<CorSingleBin>();
+	  // Construct bin of the profiled quantities. We have no information
+	  // (and no desire to add it) of sumWX of the profile, so really 
+	  // we should use a Dbn1D - but that does not work for Profile1D's.
+          profBins.push_back( YODA::ProfileBin1D((*hItr)->bin(j).xEdges(),
+	    YODA::Dbn2D( binPtrs[i]->numEntries(), binPtrs[i]->sumW(), 
+	    binPtrs[i]->sumW2(), 0., 0., binPtrs[i]->sumWX(), 0, 0)));
+        }
+	// Reset the bins of the profiles.
+	(*hItr)->reset();
+	(*hItr)->bins().clear();
+	// Add our new bins.
+	for (int j = 0, N = profBins.size(); j < N; ++j) {
+	  (*hItr)->addBin(profBins[j]);
+	}
+      } // End loop of bootstrapped correlators.
+    }
+
+    // @brief Four particle integrated cn.
+    const void cnFourInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4) const {
+      auto e2bins = e2->getBins();
+      auto e4bins = e4->getBins();
+      auto binx = e2->getBinX();
+      if (binx.size() - 1 != e2bins.size()){
+        cout << "cnFourInt: Bin size (x,y) differs!" << endl;
+        return;
+      }
+      if (binx != e4->getBinX()) {
+        cout << "Error in cnFourInt: Correlator x-binning differs!" << endl;
+        return;
+      }
+      vector<CorBinBase*> e2binPtrs;
+      vector<CorBinBase*> e4binPtrs;
+      auto cn = [&] (int i) {
+        double e22 = e2binPtrs[i]->mean() * e2binPtrs[i]->mean();
+        return e4binPtrs[i]->mean() - 2. * e22;
+      };
+      // Error calculation.
+      vector<pair<double, double> > yErr;
+      for (int j = 0, N = e2bins.size(); j < N; ++j) {
+        e2binPtrs = e2bins[j].getBinPtrs();
+        e4binPtrs = e4bins[j].getBinPtrs();
+	yErr.push_back(sampleError(cn));
+      }
+      // Put the bin ptrs back in place.
+      e2binPtrs = e2->getBinPtrs();
+      e4binPtrs = e4->getBinPtrs();
+      fillScatter(h, binx, cn, yErr);  
+    }
+
+    // @brief Four particle integrated vn
+    const void vnFourInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4) const {
+      cnFourInt(h, e2, e4);
+      nthPow(h, 0.25, -1.0);
+    }
+    
+    // @brief Six particle integrated cn.
+    const void cnSixInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4, 
+      ECorrPtr e6) const {
+      auto e2bins = e2->getBins();
+      auto e4bins = e4->getBins();
+      auto e6bins = e6->getBins();
+      auto binx = e2->getBinX();
+      if (binx.size() - 1 != e2bins.size()){
+        cout << "cnSixInt: Bin size (x,y) differs!" << endl;
+        return;
+      }
+      if (binx != e4->getBinX() || binx != e6->getBinX()) {
+        cout << "Error in cnSixInt: Correlator x-binning differs!" << endl;
+        return;
+      }
+      vector<CorBinBase*> e2binPtrs;
+      vector<CorBinBase*> e4binPtrs;
+      vector<CorBinBase*> e6binPtrs;
+      auto cn = [&] (int i) {
+        double e23 = pow(e2binPtrs[i]->mean(), 3.0); 
+        return e6binPtrs[i]->mean() - 9.*e2binPtrs[i]->mean()*e4binPtrs[i]->mean() + 
+          12.*e23;
+      };
+      // Error calculation.
+      vector<pair<double, double> > yErr;
+      for (int j = 0, N = e2bins.size(); j < N; ++j) {
+        e2binPtrs = e2bins[j].getBinPtrs();
+        e4binPtrs = e4bins[j].getBinPtrs();
+        e6binPtrs = e6bins[j].getBinPtrs();
+	yErr.push_back(sampleError(cn));
+      }
+      // Put the bin ptrs back in place.
+      e2binPtrs = e2->getBinPtrs();
+      e4binPtrs = e4->getBinPtrs();
+      e6binPtrs = e6->getBinPtrs();
+      fillScatter(h, binx, cn, yErr);
+    }
+  
+    // @brief Six particle integrated vn
+    const void vnSixInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4,
+      ECorrPtr e6) const {
+      cnSixInt(h, e2, e4, e6);
+      nthPow(h, 1./6., 0.25);
+    }
+    
+    // @brief Eight particle integrated cn.
+    const void cnEightInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4,
+      ECorrPtr e6, ECorrPtr e8) const {
+      auto e2bins = e2->getBins();
+      auto e4bins = e4->getBins();
+      auto e6bins = e6->getBins();
+      auto e8bins = e8->getBins();
+      auto binx = e2->getBinX();
+      if (binx.size() - 1 != e2bins.size()){
+        cout << "cnEightInt: Bin size (x,y) differs!" << endl;
+        return;
+      }
+      if (binx != e4->getBinX() || binx != e6->getBinX() || 
+        binx != e8->getBinX()) {
+        cout << "Error in cnEightInt: Correlator x-binning differs!" << endl;
+        return;
+      }
+      vector<CorBinBase*> e2binPtrs;
+      vector<CorBinBase*> e4binPtrs;
+      vector<CorBinBase*> e6binPtrs;
+      vector<CorBinBase*> e8binPtrs;
+      auto cn = [&] (int i ) {
+        double e22 = e2binPtrs[i]->mean() * e2binPtrs[i]->mean();
+        double e24 = e22 * e22;
+        double e42 = e4binPtrs[i]->mean() * e4binPtrs[i]->mean();
+        return e8binPtrs[i]->mean() - 16. * e6binPtrs[i]->mean() * 
+          e2binPtrs[i]->mean() - 18. * e42 + 144. * e4binPtrs[i]->mean()*e22
+	  - 144. * e24; 
+      };
+      // Error calculation.
+      vector<pair<double, double> > yErr;
+      for (int j = 0, N = e2bins.size(); j < N; ++j) {
+        e2binPtrs = e2bins[j].getBinPtrs();
+        e4binPtrs = e4bins[j].getBinPtrs();
+        e6binPtrs = e6bins[j].getBinPtrs();
+        e8binPtrs = e8bins[j].getBinPtrs();
+	yErr.push_back(sampleError(cn));
+      }
+      // Put the bin ptrs back in place.
+      e2binPtrs = e2->getBinPtrs();
+      e4binPtrs = e4->getBinPtrs();
+      e6binPtrs = e6->getBinPtrs();
+      e8binPtrs = e8->getBinPtrs();
+      fillScatter(h, binx, cn, yErr);
+    }
+
+    // @brief Eight particle integrated vn
+    const void vnEightInt(Scatter2DPtr h, ECorrPtr e2, ECorrPtr e4,
+      ECorrPtr e6, ECorrPtr e8) const {
+      cnEightInt(h, e2, e4, e6, e8);
+      nthPow(h, 1./8., -1./33.);
+    }
+
+    // @brief Two particle differential vn.
+    const void vnTwoDiff(Scatter2DPtr h, ECorrPtr e2Dif) const {
+      auto e2bins = e2Dif->getBins();
+      auto ref = e2Dif->getReference();
+      auto binx = e2Dif->getBinX();
+      if (binx.size() -1 != e2bins.size()) {
+        cout << "vnTwoDif: Bin size (x,y) differs!" << endl;
+        return;
+      }
+      vector<CorBinBase*> e2binPtrs;
+      vector<CorBinBase*> refPtrs;
+      auto vn = [&] (int i) {
+        // Test reference flow.
+        if (ref.mean() <= 0) return 0.;
+        return e2binPtrs[i]->mean() / sqrt(ref.mean());
+      };
+      // We need here a seperate error function, as we don't
+      // iterate over the reference flow.
+      auto vnerr = [&] (int i) {
+        // Test reference flow.
+	if (refPtrs[i]->mean() <=0) return 0.;
+	return e2binPtrs[i]->mean() / sqrt(refPtrs[i]->mean());
+      };
+      // Error calculation.
+      vector<pair<double, double> > yErr;
+      refPtrs = ref.getBinPtrs();
+      for (int j = 0, N = e2bins.size(); j < N; ++j) {
+        e2binPtrs = e2bins[j].getBinPtrs();
+	yErr.push_back(sampleError(vnerr));
+      }
+      // Put the e2binPtrs back in place.
+      e2binPtrs = e2Dif->getBinPtrs();
+      fillScatter(h, binx, vn);
+    }
+
+    // @brief Four particle differential vn.
+    const void vnFourDiff(Scatter2DPtr h, ECorrPtr e2Dif, 
+      ECorrPtr e4Dif) const {
+      auto e2bins = e2Dif->getBins();
+      auto e4bins = e4Dif->getBins();
+      auto ref2 = e2Dif->getReference();
+      auto ref4 = e4Dif->getReference();
+      auto binx = e2Dif->getBinX();
+      if (binx.size() - 1 != e2bins.size()){
+        cout << "vnFourDif: Bin size (x,y) differs!" << endl;
+        return;
+      }
+      if (binx != e4Dif->getBinX()) {
+        cout << "Error in vnFourDif: Correlator x-binning differs!" << endl;
+        return;
+      }
+      vector<CorBinBase*> e2binPtrs;
+      vector<CorBinBase*> e4binPtrs;
+      vector<CorBinBase*> ref2Ptrs;
+      vector<CorBinBase*> ref4Ptrs;
+      double denom = 2 * ref2.mean() * ref2.mean() - ref4.mean();
+      auto vn = [&] (int i) {
+        // Test denominator.
+        if (denom <= 0 ) return 0.;
+        return ((2 * ref2.mean() * e2bins[i].mean() - e4bins[i].mean()) /
+	  pow(denom, 0.75));
+      };
+      // We need here a seperate error function, as we don't
+      // iterate over the reference flow.
+      auto vnerr = [&] (int i) {
+        double denom2 = 2 * ref2Ptrs[i]->mean() * ref2Ptrs[i]->mean() -
+	  ref4Ptrs[i]->mean();
+	// Test denominator.
+	if (denom2 <= 0) return 0.;
+	return ((2 * ref2Ptrs[i]->mean() * e2binPtrs[i]->mean() - 
+	  e4binPtrs[i]->mean()) / pow(denom2, 0.75));
+      };
+      // Error calculation.
+      vector<pair<double, double> > yErr;
+      ref2Ptrs = ref2.getBinPtrs();
+      ref4Ptrs = ref4.getBinPtrs();
+      for (int j = 0, N = e2bins.size(); j < N; ++j) {
+        e2binPtrs = e2bins[j].getBinPtrs();
+	e4binPtrs = e4bins[j].getBinPtrs();
+	yErr.push_back(sampleError(vnerr));
+      }
+      // Put the binPtrs back in place.
+      e2binPtrs = e2Dif->getBinPtrs();
+      e4binPtrs = e4Dif->getBinPtrs();
+      fillScatter(h, binx, vn, yErr);
+    }
+  }; // End class CumulantAnalysis.
+} // End Rivet namespace.
+#endif
diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc
--- a/src/Core/AnalysisHandler.cc
+++ b/src/Core/AnalysisHandler.cc
@@ -1,556 +1,553 @@
 // -*- 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::
   mergeYodas(const vector<string> & aofiles, bool equiv) {
     vector< vector<AnalysisObjectPtr> > aosv;
     vector<double> xsecs;
     vector<double> xsecerrs;
     vector<CounterPtr> sows;
     set<string> ananames;
 
     // 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;
       _eventcounter.reset();
       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 {
             string ananame = split(ao->path(), "/")[0];
             if ( ananames.insert(ananame).second ) addAnalysis(ananame);
             aos.push_back(ao);
           }
         }
         if ( !xsec || !sow ) {
           MSG_ERROR( "Error in AnalysisHandler::mergeYodas: The file " << file
                      << " did not contain weights and cross section info.");
           exit(1);
         }
         xsecs.push_back(xsec->point(0).x());
         xsecerrs.push_back(sqr(xsec->point(0).xErrAvg()));
         _eventcounter += *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;
         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;
   }
 
 
 }
diff --git a/src/Tools/Correlators.cc b/src/Tools/Correlators.cc
new file mode 100644
--- /dev/null
+++ b/src/Tools/Correlators.cc
@@ -0,0 +1,233 @@
+// -*- C++ -*-
+#include "Rivet/Tools/Correlators.hh"
+
+namespace Rivet {
+  
+  // Constructor
+  Correlators::Correlators(const ParticleFinder& fsp, int nMaxIn, 
+    int pMaxIn, vector<double> pTbinEdgesIn) : 
+      nMax(nMaxIn + 1), pMax(pMaxIn + 1), pTbinEdges(pTbinEdgesIn) {
+    setName("Correlators");
+    declareProjection(fsp, "FS");
+    isPtDiff   = !pTbinEdges.empty();
+    if (isPtDiff) {
+      vector<double>::iterator underflow = pTbinEdges.begin();
+      pTbinEdges.insert(underflow,pTbinEdges[0]-1);
+    } 
+    setToZero();
+  }
+
+  // Alternative constructor.
+  Correlators::Correlators(const ParticleFinder& fsp, int nMaxIn, 
+    int pMaxIn, const Scatter2DPtr hIn) : nMax(nMaxIn + 1), pMax(pMaxIn + 1) {
+    for (auto b : hIn->points()) pTbinEdges.push_back(b.xMin());
+    pTbinEdges.push_back(hIn->points().back().xMax());
+    setName("Correlators");
+    declareProjection(fsp, "FS");
+    isPtDiff   = !pTbinEdges.empty();
+    if (isPtDiff) {
+      vector<double>::iterator underflow = pTbinEdges.begin();
+      pTbinEdges.insert(underflow,pTbinEdges[0]-1);
+    } 
+    setToZero();
+  }
+   
+
+  // Set all elements in vectors to zero
+  void Correlators::setToZero(){
+    vector< complex<double> > pTmp(pMax, _ZERO);
+    Vec2D qTmp(nMax, pTmp);
+    qVec = qTmp;
+    if (isPtDiff) {
+      pVec.erase(pVec.begin(), pVec.end());
+      for (double pT : pTbinEdges)
+         pVec.insert(pair<double, Vec2D>(pT, qVec));
+    }
+  }
+
+  // Functions for output:
+  const pair<double,double> Correlators::intCorrelator(vector<int> n) const {
+    // Create vector of zeros for normalisation and vector of initial
+    // powers
+    int m = n.size();
+    vector<int> powers(m, 1);
+    vector<int> zeros(m, 0);
+    complex<double> num = recCorr(m, n, powers, false);
+    complex<double> den = recCorr(m, zeros, powers, false);
+    pair<double, double> ret;
+    ret.second = (den.real() < _TINY) ? 0. : den.real();
+    ret.first = num.real();
+    return ret;
+  } 
+  
+  const vector<pair<double,double>> Correlators::pTBinnedCorrelators(vector<int> n,
+    bool overflow) const {
+    // Create vector of zeros for normalisation and vector of initial
+    // powers
+    if (!isPtDiff) 
+      cout << "You must book the correlator with a binning if you want to"
+	      " extract binned correlators! Failing." << endl;
+    int m = n.size();
+    vector<int> powers(m, 1);
+    vector<int> zeros(m, 0);
+    vector<pair<double,double>> ret;
+    for (double pT : pTbinEdges){
+      complex<double> num = recCorr(m, n, powers, true, pT);
+      complex<double> den = recCorr(m, zeros, powers, true, pT);
+      pair<double, double> tmp;
+      tmp.second = (den.real() < _TINY) ? 0. : den.real();
+      tmp.first = num.real();
+      ret.push_back(tmp);
+    }
+    if (!overflow) 
+      return vector<pair<double, double> > (ret.begin() + 1, ret.end() - 1);
+    return ret;
+  } 
+  
+  // M-particle correlation with eta-gap
+  const pair<double,double> Correlators::intCorrelatorGap(const Correlators& other,
+    vector<int> n1, vector<int> n2) const {
+    // Create vectors of zeros for normalisation and vectors of initial
+    // powers
+    int m1 = n1.size();
+    int m2 = n2.size();
+    // Return if too few particles in event.
+    vector<int> zero1(m1, 0);
+    vector<int> zero2(m2, 0);
+    vector<int> p1(m1, 1);
+    vector<int> p2(m2, 1);
+    complex<double> num1 = recCorr(m1, n1, p1, false);
+    complex<double> den1 = recCorr(m1, zero1, p1, false);
+    complex<double> num2 = other.recCorr(m2, n2, p2, false);
+    complex<double> den2 = other.recCorr(m2, zero2, p2, false);
+    complex<double> num  = num1 * num2;
+    complex<double> den  = den1 * den2;
+    pair<double, double> ret;
+    ret.second = (den1.real() < _TINY || den2.real() < _TINY)  ? 0. :
+	   den.real();
+    ret.first = num.real();
+    return ret;
+  } 
+  
+  // M-particle correlation with eta-gap
+  const vector<pair<double,double>> Correlators::pTBinnedCorrelatorsGap(
+    const Correlators& other, vector<int> n1, vector<int> n2,
+    bool overflow) const {
+    if (!isPtDiff) 
+      cout << "You must book the correlator with a binning if you want to"
+	      " extract binned correlators! Failing." << endl;
+    // Create vectors of zeros for normalisation and vectors of initial
+    // powers
+    int m1 = n1.size();
+    int m2 = n2.size();
+    vector<int> zero1(m1, 0);
+    vector<int> zero2(m2, 0);
+    vector<int> p1(m1, 1);
+    vector<int> p2(m2, 1);
+    vector<pair<double,double>> ret;
+    for (double pT : pTbinEdges) {
+      complex<double> num1 = recCorr(m1, n1, p1, true, pT);
+      complex<double> den1 = recCorr(m1, zero1, p1, true, pT);
+      complex<double> num2 = other.recCorr(m2, n2, p2, false);
+      complex<double> den2 = other.recCorr(m2, zero2, p2, false);
+      complex<double> num  = num1 * num2;
+      complex<double> den  = den1 * den2;
+      pair<double, double> tmp;
+      tmp.second = (den1.real() < _TINY || den2.real() < _TINY) 
+        ? 0. : den.real();
+      tmp.first = num.real();
+      ret.push_back(tmp);
+    }
+    //If we don't want to include underflow/overflow, remove them here
+    if (!overflow) 
+      return vector<pair<double, double> > (ret.begin() + 1, ret.end() - 1);
+    return ret;
+  } 
+
+
+  // Project function. Loops over array and calculates Q vectors
+  void Correlators::project(const Event& e) {
+    setToZero();
+    // @TODO: Weight could be implemented to account for non-uniform
+    // acceptance if detector simulation is needed. If not, the weight
+    // can be unity. Note that this weight is not the MC event weight, which
+    // should be used when filling histograms (as usual).
+    const double w = 1.0;;
+    const Particles& parts = applyProjection<ParticleFinder>(e, "FS").particles();
+    // Check that we have at least two particles in the event
+    if (parts.size() > 2)  {
+      for(const Particle& p : parts)
+        fillCorrelators(p, w);
+    }
+  } 
+  
+  // Calculate correlators from one particle
+  void Correlators::fillCorrelators(const Particle& p, const double& weight = 1.) {
+    for (int iN = 0; iN < nMax; ++iN) 
+    for (int iP = 0; iP < pMax; ++iP) {
+      double real = cos(iN * p.phi());
+      double imag = sin(iN * p.phi());
+      complex<double> expi(real, imag);
+      complex<double> tmp = pow(weight, iP) * expi;
+      qVec[iN][iP] += tmp;
+      if (isPtDiff) {
+        map<double, Vec2D>::iterator pTitr = pVec.lower_bound(p.pT());
+        // Move to the correct bin.
+        if (pTitr != pVec.begin()) pTitr--; 
+        pTitr->second[iN][iP] += tmp;
+      }
+    } 
+  } 
+  
+  // Two-particle correlator Eq. (19) p. 6 in Generic Fr. paper.
+  const complex<double> Correlators::twoPartCorr(int n1, int n2, int p1, 
+    int p2, double pT, bool useP) const {
+    complex<double> tmp1 = (!useP) ? getQ(n1, p1) : 
+                           getP(n1, p1, pT);
+    complex<double> tmp2 = getQ(n2, p2);
+    complex<double> tmp3 = (!useP) ? getQ(n1+n2, p1+p2) : 
+                           getP(n1+n2, p1+p2, pT); 
+    complex<double> sum  = tmp1 * tmp2 - tmp3;
+    return sum;
+  }
+
+  // Find correlators by recursion. Order = M (# of particles), 
+  // n's are harmonics, p's are the powers of the weights
+  const complex<double> Correlators::recCorr(int order, vector<int> n, 
+    vector<int> p, bool useP, double pT) const {
+    // Sanity checks
+    int nUsed = 0;
+    for (int i = 0, N = n.size(); i < N; ++i) nUsed += n[i]; 
+    if (nMax < nUsed)
+      cout <<"Requested n = " << nUsed << ", nMax = " << nMax << endl;
+    if (int(p.size()) > pMax)
+      cout << "Requested p = " << p.size() << ", pMax = " << pMax << endl;
+    // If order is 1, then return Q/p vector (important when dealing 
+    // with gaps and one side has only one particle
+    if (order < 2) 
+      return (!useP) ? getQ(n[0], p[0]) : getP(n[0], p[0], pT);
+    // Return 2-p correlator explicitly.
+    if ( order < 3 )
+      return twoPartCorr(n[0], n[1], p[0], p[1], pT, useP);      
+
+    // Else find nth order harmonics by recursion 
+    // at order M - 1
+    int orderNow           = order - 1;
+    int nNow               = n[orderNow];
+    int pNow               = p[orderNow];
+    complex<double> recNow = getQ(n[orderNow], p[orderNow]);
+    recNow                *= recCorr(orderNow, n, p, useP, pT); 
+    for (int i = 0; i < orderNow; ++i){
+      vector<int> tmpN, tmpP;
+      for (int j = 0; j < orderNow; ++j){
+        tmpN.push_back(n[j]);
+        tmpP.push_back(p[j]);
+      }
+      tmpN[i] += nNow;
+      tmpP[i] += pNow;
+      recNow -= recCorr(orderNow, tmpN, tmpP, useP, pT);
+    }
+    return recNow;
+  }
+
+}
diff --git a/src/Tools/Makefile.am b/src/Tools/Makefile.am
--- a/src/Tools/Makefile.am
+++ b/src/Tools/Makefile.am
@@ -1,22 +1,23 @@
 SUBDIRS = Nsubjettiness
 
 noinst_LTLIBRARIES = libRivetTools.la
 
 libRivetTools_la_SOURCES = \
   BinnedHistogram.cc \
+  Correlators.cc \
   Cuts.cc \
   JetUtils.cc \
   Random.cc \
   Logging.cc \
   ParticleUtils.cc \
   ParticleName.cc \
   Percentile.cc \
   RivetYODA.cc \
   RivetMT2.cc \
   RivetPaths.cc \
   Utils.cc \
   binreloc.c
 
 dist_noinst_HEADERS = binreloc.h lester_mt2_bisect.hh
 
 libRivetTools_la_CPPFLAGS = $(AM_CPPFLAGS) -DENABLE_BINRELOC -DDEFAULTDATADIR=\"$(datadir)\" -DDEFAULTLIBDIR=\"$(libdir)\"