diff --git a/include/Rivet/AnalysisHandler.hh b/include/Rivet/AnalysisHandler.hh
--- a/include/Rivet/AnalysisHandler.hh
+++ b/include/Rivet/AnalysisHandler.hh
@@ -1,315 +1,318 @@
 // -*- C++ -*-
 #ifndef RIVET_RivetHandler_HH
 #define RIVET_RivetHandler_HH
 
 #include "Rivet/Config/RivetCommon.hh"
 #include "Rivet/Particle.hh"
 #include "Rivet/AnalysisLoader.hh"
 #include "Rivet/Tools/RivetYODA.hh"
 
 namespace Rivet {
 
 
   // Forward declaration and smart pointer for Analysis
   class Analysis;
   typedef std::shared_ptr<Analysis> AnaHandle;
 
 
   // Needed to make smart pointers compare equivalent in the STL set
   struct CmpAnaHandle {
     bool operator() (const AnaHandle& a, const AnaHandle& b) const {
       return a.get() < b.get();
     }
   };
 
 
   /// A class which handles a number of analysis objects to be applied to
   /// generated events. An {@link Analysis}' AnalysisHandler is also responsible
   /// for handling the final writing-out of histograms.
   class AnalysisHandler {
   public:
 
     /// @name Constructors and destructors. */
     //@{
 
     /// Preferred constructor, with optional run name.
     AnalysisHandler(const string& runname="");
 
     /// @brief Destructor
     /// The destructor is not virtual, as this class should not be inherited from.
     ~AnalysisHandler();
 
     //@}
 
 
   private:
 
     /// Get a logger object.
     Log& getLog() const;
 
 
   public:
 
     /// @name Run properties
     //@{
 
     /// Get the name of this run.
     string runName() const { return _runname; }
 
 
     /// Get the number of events seen. Should only really be used by external
     /// steering code or analyses in the finalize phase.
     size_t numEvents() const { return _eventcounter.numEntries(); }
 
     /// @brief Access the sum of the event weights seen
     ///
     /// This is the weighted equivalent of the number of events. It should only
     /// be used by external steering code or analyses in the finalize phase.
     double sumW() const { return _eventcounter.sumW(); }
     /// Access to the sum of squared-weights
     double sumW2() const { return _eventcounter.sumW2(); }
 
     /// @brief Compatibility alias for sumOfWeights
     ///
     /// @deprecated Prefer sumW
     double sumOfWeights() const { return sumW(); }
 
     /// @brief Set the sum of weights
     ///
     /// This is useful if Rivet is steered externally and
     /// the analyses are run for a sub-contribution of the events
     /// (but of course have to be normalised to the total sum of weights)
     ///
     /// @todo What about the sumW2 term? That needs to be set coherently. Need a
     /// new version, with all three N,sumW,sumW2 numbers (or a counter)
     /// supplied.
     ///
     /// @deprecated Weight sums are no longer tracked this way...
     void setSumOfWeights(const double& sum) {
       //_sumOfWeights = sum;
       throw Error("Can't set sum of weights independently, since it's now tracked by a Counter. "
                   "Please contact the Rivet authors if you need this.");
     }
 
 
     /// Is cross-section information required by at least one child analysis?
     /// @deprecated Should no-longer be an issue: does any generator not write the cross-section?
     bool needCrossSection() const;
     /// Whether the handler knows about a cross-section.
     /// @deprecated Should no-longer be an issue: does any generator not write the cross-section?
     bool hasCrossSection() const;
 
     /// Get the cross-section known to the handler.
     double crossSection() const { return _xs; }
 
     /// Set the cross-section for the process being generated.
     AnalysisHandler& setCrossSection(double xs, double xserr=0);
 
 
     /// Set the beam particles for this run
     AnalysisHandler& setRunBeams(const ParticlePair& beams) {
       _beams = beams;
       MSG_DEBUG("Setting run beams = " << beams << " @ " << sqrtS()/GeV << " GeV");
       return *this;
     }
 
     /// Get the beam particles for this run, usually determined from the first event.
     const ParticlePair& beams() const { return _beams; }
 
     /// Get beam IDs for this run, usually determined from the first event.
     /// @deprecated Use standalone beamIds(ah.beams()), to clean AH interface
     PdgIdPair beamIds() const;
 
     /// Get energy for this run, usually determined from the first event.
     /// @deprecated Use standalone sqrtS(ah.beams()), to clean AH interface
     double sqrtS() const;
 
     /// Setter for _ignoreBeams
     void setIgnoreBeams(bool ignore=true);
 
     //@}
 
 
     /// @name Handle analyses
     //@{
 
     /// Get a list of the currently registered analyses' names.
     std::vector<std::string> analysisNames() const;
 
     /// Get the collection of currently registered analyses.
     const std::set<AnaHandle, CmpAnaHandle>& analyses() const {
       return _analyses;
     }
 
     /// Get a registered analysis by name.
     const AnaHandle analysis(const std::string& analysisname) const;
 
 
     /// Add an analysis to the run list by object
     AnalysisHandler& addAnalysis(Analysis* analysis);
 
     /// @brief Add an analysis to the run list using its name.
     ///
     /// The actual Analysis to be used will be obtained via
     /// AnalysisLoader::getAnalysis(string).  If no matching analysis is found,
     /// no analysis is added (i.e. the null pointer is checked and discarded.
     AnalysisHandler& addAnalysis(const std::string& analysisname);
 
     /// @brief Add an analysis with a map of analysis options.
     AnalysisHandler& addAnalysis(const std::string& analysisname, std::map<string, string> pars);
     
     /// @brief Add analyses to the run list using their names.
     ///
     /// The actual {@link Analysis}' to be used will be obtained via
     /// AnalysisHandler::addAnalysis(string), which in turn uses
     /// AnalysisLoader::getAnalysis(string). If no matching analysis is found
     /// for a given name, no analysis is added, but also no error is thrown.
     AnalysisHandler& addAnalyses(const std::vector<std::string>& analysisnames);
 
 
     /// Remove an analysis from the run list using its name.
     AnalysisHandler& removeAnalysis(const std::string& analysisname);
 
     /// Remove analyses from the run list using their names.
     AnalysisHandler& removeAnalyses(const std::vector<std::string>& analysisnames);
 
     //@}
 
 
     /// @name Main init/execute/finalise
     //@{
 
     /// Initialize a run, with the run beams taken from the example event.
     void init(const GenEvent& event);
 
     /// @brief Analyze the given \a event by reference.
     ///
     /// This function will call the AnalysisBase::analyze() function of all
     /// included analysis objects.
     void analyze(const GenEvent& event);
 
     /// @brief Analyze the given \a event by pointer.
     ///
     /// This function will call the AnalysisBase::analyze() function of all
     /// included analysis objects, after checking the event pointer validity.
     void analyze(const GenEvent* event);
 
     /// Finalize a run. This function calls the AnalysisBase::finalize()
     /// functions of all included analysis objects.
     void finalize();
 
     //@}
 
 
     /// @name Histogram / data object access
     //@{
 
     /// Add a vector of analysis objects to the current state.
     void addData(const std::vector<AnalysisObjectPtr>& aos);
 
     /// Read analysis plots into the histo collection (via addData) from the named file.
     void readData(const std::string& filename);
 
     /// Get all analyses' plots as a vector of analysis objects.
     std::vector<AnalysisObjectPtr> getData(bool includeorphans = false,
                                            bool includetmps = false) const;
-
+    
+    /// Get all finalized analyses' plots as a vector of analysis objects.
+    std::vector<AnalysisObjectPtr> getFinalizedData(bool includeorphans = false) const;
+    
     /// Write all analyses' plots (via getData) to the named file.
     void writeData(const std::string& filename) const;
 
     /// Tell Rivet to dump intermediate result to a file named @a
     /// dumpfile every @a period'th event. If @period is not positive,
     /// no dumping will be done.
     void dump(string dumpfile, int period) {
       _dumpPeriod = period;
       _dumpFile = dumpfile;
     }
 
     /// Take the vector of yoda files and merge them together using
     /// the cross section and weight information provided in each
     /// file. Each file in @a aofiles is assumed to have been produced
     /// by Rivet. By default the files are assumed to contain
     /// different processes (or the same processs but mutually
     /// exclusive cuts), but if @a equiv if ture, the files are
     /// assumed to contain output of completely equivalent (but
     /// statistically independent) Rivet runs. The corresponding
     /// analyses will be loaded and their analysis objects will be
     /// filled with the merged result. finalize() will be run on each
     /// relevant anslysis. The resulting YODA file can then be rwitten
     /// out by writeData(). If delopts is non-empty, it is assumed to
     /// contain names different options to be merged into the same
     /// analysis objects.
     void mergeYodas(const vector<string> & aofiles,
                     const vector<string> & delopts = vector<string>(),
                     bool equiv = false);
 
     /// Helper function to strip specific options from data object paths.
     void stripOptions(AnalysisObjectPtr ao,
                       const vector<string> & delopts) const;
 
     //@}
 
   private:
 
     /// The collection of Analysis objects to be used.
     set<AnaHandle, CmpAnaHandle> _analyses;
 
     /// A vector of pre-loaded object which do not have a valid
     /// Analysis plugged in.
     vector<AnalysisObjectPtr> _orphanedPreloads;
 
     /// A vector containing copies of analysis objects after
     /// finalize() has been run.
     vector<AnalysisObjectPtr> _finalizedAOs;
 
     /// @name Run properties
     //@{
 
     /// Run name
     std::string _runname;
 
     /// Event counter
     Counter _eventcounter;
 
     /// Cross-section known to AH
     double _xs, _xserr;
 
     /// Beams used by this run.
     ParticlePair _beams;
 
     /// Flag to check if init has been called
     bool _initialised;
 
     /// Flag whether input event beams should be ignored in compatibility check
     bool _ignoreBeams;
 
     /// Determines how often Rivet runs finalize() and writes the
     /// result to a YODA file.
     int _dumpPeriod;
 
     /// The name of a YODA file to which Rivet periodically dumps
     /// results.
     string _dumpFile;
 
     /// Flag to indicate periodic dumping is in progress
     bool _dumping;
 
     //@}
 
 
   private:
 
     /// The assignment operator is private and must never be called.
     /// In fact, it should not even be implemented.
     AnalysisHandler& operator=(const AnalysisHandler&);
 
     /// The copy constructor is private and must never be called.  In
     /// fact, it should not even be implemented.
     AnalysisHandler(const AnalysisHandler&);
 
   };
 
 
 }
 
 #endif
diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc
--- a/src/Core/AnalysisHandler.cc
+++ b/src/Core/AnalysisHandler.cc
@@ -1,575 +1,591 @@
 // -*- 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) ) {
       // TODO: This should be possible to do in a nicer way, with a flag etc.
       if (ao->path().find("/FINAL") != std::string::npos) continue;
       auto aoit = backupAOs.find(ao->path());
       if ( aoit == backupAOs.end() ) {
         AnaHandle ana = analysis(split(ao->path(), "/")[0]);
         if ( ana ) ana->removeAnalysisObject(ao->path());
       } else
         copyao(aoit->second, ao);
     }
 
     // Print out number of events processed
     const int nevts = _eventcounter.numEntries();
     MSG_INFO("Processed " << nevts << " event" << (nevts != 1 ? "s" : ""));
 
     // // Delete analyses
     // MSG_DEBUG("Deleting analyses");
     // _analyses.clear();
 
     // Print out MCnet boilerplate
     cout << endl;
     cout << "The MCnet usage guidelines apply to Rivet: see http://www.montecarlonet.org/GUIDELINES" << endl;
     cout << "Please acknowledge plots made with Rivet analyses, and cite arXiv:1003.0694 (http://arxiv.org/abs/1003.0694)" << endl;
   }
 
   AnalysisHandler& AnalysisHandler::addAnalysis(const string& analysisname, std::map<string, string> pars) {
      // Make an option handle.
     std::string parHandle = "";
     for (map<string, string>::iterator par = pars.begin(); par != pars.end(); ++par) {
       parHandle +=":";
       parHandle += par->first + "=" + par->second;
     }
     return addAnalysis(analysisname + parHandle);
 
   }
 
   AnalysisHandler& AnalysisHandler::addAnalysis(const string& analysisname) {
     // Check for a duplicate analysis
     /// @todo Might we want to be able to run an analysis twice, with different params?
     ///       Requires avoiding histo tree clashes, i.e. storing the histos on the analysis objects.
     string ananame = analysisname;
     vector<string> anaopt = split(analysisname, ":");
     if ( anaopt.size() > 1 ) ananame = anaopt[0];
     AnaHandle analysis( AnalysisLoader::getAnalysis(ananame) );
     if (analysis.get() != 0) { // < Check for null analysis.
       MSG_DEBUG("Adding analysis '" << analysisname << "'");
       map<string,string> opts;
       for ( int i = 1, N = anaopt.size(); i < N; ++i ) {
         vector<string> opt = split(anaopt[i], "=");
         if ( opt.size() != 2 ) {
           MSG_WARNING("Error in option specification. Skipping analysis "
                       << analysisname);
           return *this;
         }
         if ( !analysis->info().validOption(opt[0], opt[1]) ) {
           MSG_WARNING("Cannot set option '" << opt[0] << "' to '" << opt[1]
                       << "'. Skipping analysis " << analysisname);
           return *this;
         }
         opts[opt[0]] = opt[1];
       }
       for ( auto opt: opts) {
         analysis->_options[opt.first] = opt.second;
         analysis->_optstring += ":" + opt.first + "=" + opt.second;
       }
       for (const AnaHandle& a : _analyses) {
         if (a->name() == analysis->name() ) {
           MSG_WARNING("Analysis '" << analysisname << "' already registered: skipping duplicate");
           return *this;
         }
       }
       analysis->_analysishandler = this;
       _analyses.insert(analysis);
     } else {
       MSG_WARNING("Analysis '" << analysisname << "' not found.");
     }
     // MSG_WARNING(_analyses.size());
     // for (const AnaHandle& a : _analyses) MSG_WARNING(a->name());
     return *this;
   }
 
 
   AnalysisHandler& AnalysisHandler::removeAnalysis(const string& analysisname) {
     std::shared_ptr<Analysis> toremove;
     for (const AnaHandle a : _analyses) {
       if (a->name() == analysisname) {
         toremove = a;
         break;
       }
     }
     if (toremove.get() != 0) {
       MSG_DEBUG("Removing analysis '" << analysisname << "'");
       _analyses.erase(toremove);
     }
     return *this;
   }
 
 
   /////////////////////////////
 
 
   void AnalysisHandler::addData(const std::vector<AnalysisObjectPtr>& aos) {
     for (const AnalysisObjectPtr ao : aos) {
       string path = ao->path();
       if ( path.substr(0, 5) != "/RAW/" ) {
         _orphanedPreloads.push_back(ao);
         continue;
       }
  
       path = path.substr(4);
       ao->setPath(path);
       if (path.size() > 1) { // path > "/"
         try {
           const string ananame =  split(path, "/")[0];
           AnaHandle a = analysis(ananame);
           a->addAnalysisObject(ao); /// @todo Need to statistically merge...
         } catch (const Error& e) {
           MSG_TRACE("Adding analysis object " << path <<
                     " to the list of orphans.");
           _orphanedPreloads.push_back(ao);
         }
       }
     }
   }
 
   void AnalysisHandler::stripOptions(AnalysisObjectPtr ao,
                                      const vector<string> & delopts) const {
     string path = ao->path();
     string ananame = split(path, "/")[0];
     vector<string> anaopts = split(ananame, ":");
     for ( int i = 1, N = anaopts.size(); i < N; ++i )
       for ( auto opt : delopts )
         if ( opt == "*" || anaopts[i].find(opt + "=") == 0 )
           path.replace(path.find(":" + anaopts[i]), (":" + anaopts[i]).length(), "");
     ao->setPath(path);
   }
    
 
 
 
   void AnalysisHandler::
   mergeYodas(const vector<string> & aofiles, const vector<string> & delopts, bool equiv) {
     vector< vector<AnalysisObjectPtr> > aosv;
     vector<double> xsecs;
     vector<double> xsecerrs;
     vector<CounterPtr> sows;
     set<string> ananames;
      _eventcounter.reset();
  
     // First scan all files and extract analysis objects and add the
     // corresponding anayses..
     for ( auto file : aofiles ) {
       Scatter1DPtr xsec;
       CounterPtr sow;
 
       // For each file make sure that cross section and sum-of-weights
       // objects are present and stor all RAW ones in a vector;
       vector<AnalysisObjectPtr> aos;
       try {
         /// @todo Use new YODA SFINAE to fill the smart ptr vector directly
         vector<YODA::AnalysisObject*> aos_raw;
         YODA::read(file, aos_raw);
         for (AnalysisObject* aor : aos_raw) {
           AnalysisObjectPtr ao = AnalysisObjectPtr(aor);
           if ( ao->path().substr(0, 5) != "/RAW/" ) continue;
           ao->setPath(ao->path().substr(4));
           if ( ao->path() == "/_XSEC" )
             xsec = dynamic_pointer_cast<Scatter1D>(ao);
           else if ( ao->path() == "/_EVTCOUNT" )
             sow = dynamic_pointer_cast<Counter>(ao);
           else {
             stripOptions(ao, delopts);
             string ananame = split(ao->path(), "/")[0];
             if ( ananames.insert(ananame).second ) addAnalysis(ananame);
             aos.push_back(ao);
           }
         }
         if ( !xsec || !sow ) {
           MSG_ERROR( "Error in AnalysisHandler::mergeYodas: The file " << file
                      << " did not contain weights and cross section info.");
           exit(1);
         }
         xsecs.push_back(xsec->point(0).x());
         sows.push_back(sow);
 	xsecerrs.push_back(sqr(xsec->point(0).xErrAvg()));
         _eventcounter += *sow;
         sows.push_back(sow);
         aosv.push_back(aos);
       } catch (...) { //< YODA::ReadError&
         throw UserError("Unexpected error in reading file: " + file);
       }
     }
 
     // Now calculate the scale to be applied for all bins in a file
     // and get the common cross section and sum of weights.
     _xs = _xserr = 0.0;
     for ( int i = 0, N = sows.size(); i < N; ++i ) {
       double effnent = sows[i]->effNumEntries();
       _xs += (equiv? effnent: 1.0)*xsecs[i];
       _xserr += (equiv? sqr(effnent): 1.0)*xsecerrs[i];
     }
 
     vector<double> scales(sows.size(), 1.0);
     if ( equiv ) {
       _xs /= _eventcounter.effNumEntries();
       _xserr = sqrt(_xserr)/_eventcounter.effNumEntries();
     } else {
       _xserr = sqrt(_xserr);
       for ( int i = 0, N = sows.size(); i < N; ++i ) 
         scales[i] = (_eventcounter.sumW()/sows[i]->sumW())*(xsecs[i]/_xs);
     }
 
     // Initialize the analyses allowing them to book analysis objects.
     for (AnaHandle a : _analyses) {
       MSG_DEBUG("Initialising analysis: " << a->name());
       if ( !a->info().reentrant() )
         MSG_WARNING("Analysis " << a->name() << " has not been validated to have "
                     << "a reentrant finalize method. The result is unpredictable.");
       try {
         // Allow projection registration in the init phase onwards
         a->_allowProjReg = true;
         cerr << "sqrtS " << sqrtS() << endl;
         a->init();
         //MSG_DEBUG("Checking consistency of analysis: " << a->name());
         //a->checkConsistency();
       } catch (const Error& err) {
         cerr << "Error in " << a->name() << "::init method: " << err.what() << endl;
         exit(1);
       }
       MSG_DEBUG("Done initialising analysis: " << a->name());
     }
     _initialised = true;
     // Get a list of all anaysis objects to handle.
     map<string,AnalysisObjectPtr> current;
     for ( auto ao : getData(false, true) ) current[ao->path()] = ao;
     // Go through all objects to be merged and add them to current
     // after appropriate scaling.
     for ( int i = 0, N = aosv.size(); i < N; ++i)
       for ( auto ao : aosv[i] ) {
         if ( ao->path() == "/_XSEC" || ao->path() == "_EVTCOUNT" ) continue;
 	auto aoit = current.find(ao->path());
         if ( aoit == current.end() ) {
           MSG_WARNING("" << ao->path() << " was not properly booked.");
           continue;
         }
         if ( !addaos(aoit->second, ao, scales[i]) )
           MSG_WARNING("Cannot merge objects with path " << ao->path()
                       <<" of type " << ao->annotation("Type") );
       }
     // Now we can simply finalize() the analysis, leaving the
     // controlling program to write it out some yoda-file.
     finalize();
 
   }
 
 
   void AnalysisHandler::readData(const string& filename) {
     vector<AnalysisObjectPtr> aos;
     try {
       /// @todo Use new YODA SFINAE to fill the smart ptr vector directly
       vector<YODA::AnalysisObject*> aos_raw;
       YODA::read(filename, aos_raw);
       for (AnalysisObject* aor : aos_raw) aos.push_back(AnalysisObjectPtr(aor));
     } catch (...) { //< YODA::ReadError&
       throw UserError("Unexpected error in reading file: " + filename);
     }
     if (!aos.empty()) addData(aos);
   }
 
 
   vector<AnalysisObjectPtr> AnalysisHandler::
   getData(bool includeorphans, bool includetmps) const {
     vector<AnalysisObjectPtr> rtn;
     // Event counter
     rtn.push_back( make_shared<Counter>(_eventcounter) );
     // Cross-section + err as scatter
     YODA::Scatter1D::Points pts; pts.insert(YODA::Point1D(_xs, _xserr));
     rtn.push_back( make_shared<Scatter1D>(pts, "/_XSEC") );
     // Analysis histograms
     for (const AnaHandle a : analyses()) {
       vector<AnalysisObjectPtr> aos = a->analysisObjects();
       // MSG_WARNING(a->name() << " " << aos.size());
       for (const AnalysisObjectPtr ao : aos) {
         // Exclude paths from final write-out if they contain a "TMP" layer (i.e. matching "/TMP/")
         /// @todo This needs to be much more nuanced for re-entrant histogramming
         if ( !includetmps && ao->path().find("/TMP/" ) != string::npos) continue;
         rtn.push_back(ao);
       }
     }
     // Sort histograms alphanumerically by path before write-out
     sort(rtn.begin(), rtn.end(), [](AnalysisObjectPtr a, AnalysisObjectPtr b) {return a->path() < b->path();});
     if ( includeorphans )
       rtn.insert(rtn.end(), _orphanedPreloads.begin(), _orphanedPreloads.end());
     return rtn;
   }
+  
+  vector<AnalysisObjectPtr> AnalysisHandler::
+  getFinalizedData(bool includeorphans) const {
+    // Finalized analysis histograms
+    vector<AnalysisObjectPtr> rtn = _finalizedAOs;
+    // 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") );
+    // 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, double xserr) {
     _xs = xs;
     _xserr = xserr;
     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;
   }
 
 
 }