diff --git a/include/Rivet/Tools/RivetYODA.hh b/include/Rivet/Tools/RivetYODA.hh
--- a/include/Rivet/Tools/RivetYODA.hh
+++ b/include/Rivet/Tools/RivetYODA.hh
@@ -1,595 +1,595 @@
 #ifndef RIVET_RIVETYODA_HH
 #define RIVET_RIVETYODA_HH
 
 #include "Rivet/Config/RivetCommon.hh"
 #include "YODA/AnalysisObject.h"
 #include "YODA/Counter.h"
 #include "YODA/Histo1D.h"
 #include "YODA/Histo2D.h"
 #include "YODA/Profile1D.h"
 #include "YODA/Profile2D.h"
 #include "YODA/Scatter1D.h"
 #include "YODA/Scatter2D.h"
 #include "YODA/Scatter3D.h"
 
 #include <map>
 #include <valarray>
 
 namespace YODA {
 
   typedef std::shared_ptr<YODA::AnalysisObject> AnalysisObjectPtr;
 
   typedef std::shared_ptr<YODA::Counter> CounterPtr;
   typedef std::shared_ptr<YODA::Histo1D> Histo1DPtr;
   typedef std::shared_ptr<YODA::Histo2D> Histo2DPtr;
   typedef std::shared_ptr<YODA::Profile1D> Profile1DPtr;
   typedef std::shared_ptr<YODA::Profile2D> Profile2DPtr;
   typedef std::shared_ptr<YODA::Scatter1D> Scatter1DPtr;
   typedef std::shared_ptr<YODA::Scatter2D> Scatter2DPtr;
   typedef std::shared_ptr<YODA::Scatter3D> Scatter3DPtr;
 
 }
 
 
 namespace Rivet {
 
 
   class AnalysisObjectWrapper {
   public:
     virtual ~AnalysisObjectWrapper() {}
 
     virtual YODA::AnalysisObject* operator->() = 0;
     virtual YODA::AnalysisObject* operator->() const = 0;
     virtual const YODA::AnalysisObject & operator*() const = 0;
 
     /// @todo Rename to setActive(idx)
     virtual void setActiveWeightIdx(unsigned int iWeight) = 0;
 
     /// @todo Set active object for finalize
     virtual void setActiveFinalWeightIdx(unsigned int iWeight) = 0;
 
     bool operator ==(const AnalysisObjectWrapper& p) { return (this == &p); }
     bool operator !=(const AnalysisObjectWrapper& p) { return (this != &p); }
 
   protected:
     /// @todo do we need this?
     // virtual void reset() = 0;
   };
 
   /// @todo
   /// implement scatter1dptr and scatter2dptr here
   /// these need to be multi-weighted eventually.
   /*
     class Scatter1DPtr : public AnalysisObjectPtr {
     public:
     Scatter1DPtr() : _persistent() { }
 
     Scatter1DPtr(size_t len_of_weightvec, const YODA::Scatter1D& p) {
     for (size_t m = 0; m < len_of_weightvec; ++m)
     _persistent.push_back(make_shared<YODA::Scatter1D>(p));
     }
 
     bool operator!() const { return !_persistent; }
     explicit operator bool() const { return bool(_persistent); }
 
     YODA::Scatter1D* operator->() { return _persistent.get(); }
 
     YODA::Scatter1D* operator->() const { return _persistent.get(); }
 
     YODA::Scatter1D & operator*() { return *_persistent; }
 
     const YODA::Scatter1D & operator*() const { return *_persistent; }
 
     protected:
     vector<YODA::Scatter1DPtr> _persistent;
     };
 
     class Scatter2DPtr : public AnalysisObjectPtr {
     public:
     Scatter2DPtr(size_t len_of_weightvec, const YODA::Scatter2D& p) {
     for (size_t m = 0; m < len_of_weightvec; ++m)
     _persistent.push_back(make_shared<YODA::Scatter2D>(p));
     }
 
     Scatter2DPtr() : _persistent() { }
 
     bool operator!() { return !_persistent; }
     explicit operator bool() { return bool(_persistent); }
 
     YODA::Scatter2D* operator->() { return _persistent.get(); }
 
     YODA::Scatter2D* operator->() const { return _persistent.get(); }
 
     YODA::Scatter2D & operator*() { return *_persistent; }
 
     const YODA::Scatter2D & operator*() const { return *_persistent; }
 
     protected:
     vector<YODA::Scatter2DPtr> _persistent;
     };
 
     class Scatter3DPtr : public AnalysisObjectPtr {
     public:
     Scatter3DPtr(size_t len_of_weightvec, const YODA::Scatter3D& p) {
     for (size_t m = 0; m < len_of_weightvec; ++m)
     _persistent.push_back(make_shared<YODA::Scatter3D>(p));
     }
 
     Scatter3DPtr() : _persistent() { }
 
     bool operator!() { return !_persistent; }
     explicit operator bool() { return bool(_persistent); }
 
     YODA::Scatter3D* operator->() { return _persistent.get(); }
 
     YODA::Scatter3D* operator->() const { return _persistent.get(); }
 
     YODA::Scatter3D & operator*() { return *_persistent; }
 
     const YODA::Scatter3D & operator*() const { return *_persistent; }
 
     protected:
     vector<YODA::Scatter3DPtr> _persistent;
     };
   */
 
 
   class MultiweightAOWrapper : public AnalysisObjectWrapper {
 
   public:
     using Inner = YODA::AnalysisObject;
 
     virtual void newSubEvent() = 0;
 
     virtual void pushToPersistent(const vector<std::valarray<double> >& weight) = 0;
-    virtual void pushToFinal(const vector<std::valarray<double> >& weight) = 0;
+    virtual void pushToFinal() = 0;
 
     virtual YODA::AnalysisObjectPtr activeYODAPtr() const = 0;
   };
 
 
   using Weight = double;
 
   template <class T>
   using Fill = pair<typename T::FillType, Weight>;
 
   template <class T>
   using Fills = multiset<Fill<T>>;
 
 
   // TODO TODO TODO
   // need to override the old fill method too!
   // otherwise we can't intercept existing fill calls in analysis code
   // TODO TODO TODO
 
 
   /// Wrappers for analysis objects to store all fills unaggregated, until collapsed
   template <class T>
   class TupleWrapper;
 
   template<>
   class TupleWrapper<YODA::Counter> : public YODA::Counter {
   public:
     typedef shared_ptr<TupleWrapper<YODA::Counter>> Ptr;
     TupleWrapper(const YODA::Counter & h) : YODA::Counter(h) {}
     // todo: do we need to deal with users using fractions directly?
     void fill( double weight=1.0, double fraction=1.0 ) {
       fills_.insert( {YODA::Counter::FillType(),weight} );
     }
     void reset() { fills_.clear(); }
     const Fills<YODA::Counter> & fills() const { return fills_; }
   private:
     // x / weight pairs
     Fills<YODA::Counter> fills_;
   };
 
   template<>
   class TupleWrapper<YODA::Histo1D> : public YODA::Histo1D {
   public:
     typedef shared_ptr<TupleWrapper<YODA::Histo1D>> Ptr;
     TupleWrapper(const YODA::Histo1D & h) : YODA::Histo1D(h) {}
     // todo: do we need to deal with users using fractions directly?
     void fill( double x, double weight=1.0, double fraction=1.0 ) {
       if ( std::isnan(x) ) throw YODA::RangeError("X is NaN");
       fills_.insert( { x , weight } );
     }
     void reset() { fills_.clear(); }
     const Fills<YODA::Histo1D> & fills() const { return fills_; }
   private:
     // x / weight pairs
     Fills<YODA::Histo1D> fills_;
   };
 
   template<>
   class TupleWrapper<YODA::Profile1D> : public YODA::Profile1D {
   public:
     typedef shared_ptr<TupleWrapper<YODA::Profile1D>> Ptr;
     TupleWrapper(const YODA::Profile1D & h) : YODA::Profile1D(h) {}
     // todo: do we need to deal with users using fractions directly?
     void fill( double x, double y, double weight=1.0, double fraction=1.0 ) {
       if ( std::isnan(x) ) throw YODA::RangeError("X is NaN");
       if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN");
       fills_.insert( { YODA::Profile1D::FillType{x,y}, weight } );
     }
     void reset() { fills_.clear(); }
     const Fills<YODA::Profile1D> & fills() const { return fills_; }
   private:
     // x / weight pairs
     Fills<YODA::Profile1D> fills_;
   };
 
 
   template<>
   class TupleWrapper<YODA::Histo2D> : public YODA::Histo2D {
   public:
     typedef shared_ptr<TupleWrapper<YODA::Histo2D>> Ptr;
     TupleWrapper(const YODA::Histo2D & h) : YODA::Histo2D(h) {}
     // todo: do we need to deal with users using fractions directly?
     void fill( double x, double y, double weight=1.0, double fraction=1.0 ) {
       if ( std::isnan(x) ) throw YODA::RangeError("X is NaN");
       if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN");
       fills_.insert( { YODA::Histo2D::FillType{x,y}, weight } );
     }
     void reset() { fills_.clear(); }
     const Fills<YODA::Histo2D> & fills() const { return fills_; }
   private:
     // x / weight pairs
     Fills<YODA::Histo2D> fills_;
   };
 
   template<>
   class TupleWrapper<YODA::Profile2D> : public YODA::Profile2D {
   public:
     typedef shared_ptr<TupleWrapper<YODA::Profile2D>> Ptr;
     TupleWrapper(const YODA::Profile2D & h) : YODA::Profile2D(h) {}
     // todo: do we need to deal with users using fractions directly?
     void fill( double x, double y, double z, double weight=1.0, double fraction=1.0 ) {
       if ( std::isnan(x) ) throw YODA::RangeError("X is NaN");
       if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN");
       if ( std::isnan(z) ) throw YODA::RangeError("Z is NaN");
       fills_.insert( { YODA::Profile2D::FillType{x,y,z}, weight } );
     }
     void reset() { fills_.clear(); }
     const Fills<YODA::Profile2D> & fills() const { return fills_; }
   private:
     // x / weight pairs
     Fills<YODA::Profile2D> fills_;
   };
 
   template<>
   class TupleWrapper<YODA::Scatter1D> : public YODA::Scatter1D {
   public:
     typedef shared_ptr<TupleWrapper<YODA::Scatter1D>> Ptr;
     TupleWrapper(const YODA::Scatter1D & h) : YODA::Scatter1D(h) {}
   };
 
   template<>
   class TupleWrapper<YODA::Scatter2D> : public YODA::Scatter2D {
   public:
     typedef shared_ptr<TupleWrapper<YODA::Scatter2D>> Ptr;
     TupleWrapper(const YODA::Scatter2D & h) : YODA::Scatter2D(h) {}
   };
 
   template<>
   class TupleWrapper<YODA::Scatter3D> : public YODA::Scatter3D {
   public:
     typedef shared_ptr<TupleWrapper<YODA::Scatter3D>> Ptr;
     TupleWrapper(const YODA::Scatter3D & h) : YODA::Scatter3D(h) {}
   };
 
 
 
   template <class T>
   class Wrapper : public MultiweightAOWrapper {
     friend class Analysis;
 
   public:
 
     using Inner = T;
     /* @todo
      * some things are not really well-defined here
      * for instance: fill() in the finalize() method and integral() in
      * the analyze() method.
      */
 
     Wrapper() = default;
 
     Wrapper(const vector<string>& weightnames, const T & p);
 
     ~Wrapper();
 
     typename T::Ptr active() const;
 
     /* @todo this probably need to loop over all? */
     bool operator!() const { return !_active; } // Don't use active() here, assert will catch
 
     explicit operator bool() const { return static_cast<bool>(_active); } // Don't use active() here, assert will catch
 
     T * operator->() { return active().get(); }
 
     T * operator->() const { return active().get(); }
 
     T & operator*() { return *active(); }
 
     const T & operator*() const { return *active(); }
     /* @todo
      * these need to be re-thought out.
 
      void reset() { active()->reset(); }
     */
 
     /* @todo
      * these probably need to loop over all?
      * do we even want to provide equality?
      */
     /* @todo
      * how about no.
      friend bool operator==(Wrapper a, Wrapper b){
      if (a._persistent.size() != b._persistent.size())
      return false;
 
      for (size_t i = 0; i < a._persistent.size(); i++) {
      if (a._persistent.at(i) != b._persistent.at(i)) {
      return false;
      }
      }
 
      return true;
      }
 
      friend bool operator!=(Wrapper a, Wrapper b){
      return !(a == b);
      }
 
 
      friend bool operator<(Wrapper a, Wrapper b){
      if (a._persistent.size() >= b._persistent.size())
      return false;
 
      for (size_t i = 0; i < a._persistent.size(); i++) {
      if (*(a._persistent.at(i)) >= *(b._persistent.at(i))) {
      return false;
      }
      }
 
      return true;
      }
     */
 
 
   private:
     void setActiveWeightIdx(unsigned int iWeight) {
       _active = _persistent.at(iWeight);
     }
 
     void setActiveFinalWeightIdx(unsigned int iWeight) {
       _active = _final.at(iWeight);
     }
 
 
     /* this is for dev only---we shouldn't need this in real runs. */
     void unsetActiveWeight() { _active.reset(); }
 
     void newSubEvent();
 
     virtual YODA::AnalysisObjectPtr activeYODAPtr() const { return _active; }
 
     const vector<typename T::Ptr> & persistent() const { return _persistent; }
 
     const vector<typename T::Ptr> & final() const { return _final; }
 
     /* to be implemented for each type */
     void pushToPersistent(const vector<std::valarray<double> >& weight);
-    void pushToFinal(const vector<std::valarray<double> >& weight);
+    void pushToFinal();
 
 
     /* M of these, one for each weight */
     vector<typename T::Ptr> _persistent;
 
     /* This is the copy of _persistent that will be passed to finalize(). */
     vector<typename T::Ptr> _final;
 
     /* N of these, one for each event in evgroup */
     vector<typename TupleWrapper<T>::Ptr> _evgroup;
 
     typename T::Ptr _active;
 
     // do we need implicit cast?
     // operator typename T::Ptr () {
     //     return _active;
     // }
 
     friend class AnalysisHandler;
   };
 
 
   /// We need our own shared_ptr class, so we can dispatch -> and *
   /// all the way down to the inner YODA analysis objects
   ///
   /// TODO: provide remaining functionality that shared_ptr has (not needed right now)
   ///
   template <typename T>
   class rivet_shared_ptr {
   public:
     typedef T value_type;
 
     rivet_shared_ptr() = default;
 
     rivet_shared_ptr(decltype(nullptr)) : _p(nullptr) {}
 
     /// Convenience constructor, pass through to the Wrapper constructor
     rivet_shared_ptr(const vector<string>& weightNames, const typename T::Inner & p)
       : _p( make_shared<T>(weightNames, p) )
     {}
 
     template <typename U>
     rivet_shared_ptr(const shared_ptr<U> & p)
       : _p(p)
     {}
 
     template <typename U>
     rivet_shared_ptr(const rivet_shared_ptr<U> & p)
       : _p(p.get())
     {}
 
     // Goes right through to the active YODA object's members
     T & operator->()                            { return  *_p; }
     const T & operator->() const                { return  *_p; }
 
     // The active YODA object
     typename T::Inner & operator*()             { return **_p; }
     const typename T::Inner & operator*() const { return **_p; }
 
     bool operator!() const { return !_p || !(*_p);   }
     explicit operator bool()  const { return _p && bool(*_p); }
 
     template <typename U>
     bool operator==(const rivet_shared_ptr<U> & other) const {
       return _p == other._p;
     }
 
     template <typename U>
     bool operator!=(const rivet_shared_ptr<U> & other) const {
       return _p != other._p;
     }
 
     template <typename U>
     bool operator<(const rivet_shared_ptr<U> & other) const {
       return _p < other._p;
     }
 
     template <typename U>
     bool operator>(const rivet_shared_ptr<U> & other) const {
       return _p > other._p;
     }
 
     template <typename U>
     bool operator<=(const rivet_shared_ptr<U> & other) const {
       return _p <= other._p;
     }
 
     template <typename U>
     bool operator>=(const rivet_shared_ptr<U> & other) const {
       return _p >= other._p;
     }
 
     shared_ptr<T> get() const { return _p; }
 
   private:
     shared_ptr<T> _p;
 
   };
 
 
 
   // every object listed here needs a virtual fill method in YODA,
   // otherwise the Tuple fakery won't work.
 
   using MultiweightAOPtr = rivet_shared_ptr<MultiweightAOWrapper>;
 
   using Histo1DPtr   = rivet_shared_ptr<Wrapper<YODA::Histo1D>>;
   using Histo2DPtr   = rivet_shared_ptr<Wrapper<YODA::Histo2D>>;
   using Profile1DPtr = rivet_shared_ptr<Wrapper<YODA::Profile1D>>;
   using Profile2DPtr = rivet_shared_ptr<Wrapper<YODA::Profile2D>>;
   using CounterPtr   = rivet_shared_ptr<Wrapper<YODA::Counter>>;
   using Scatter1DPtr = rivet_shared_ptr<Wrapper<YODA::Scatter1D>>;
   using Scatter2DPtr = rivet_shared_ptr<Wrapper<YODA::Scatter2D>>;
   using Scatter3DPtr = rivet_shared_ptr<Wrapper<YODA::Scatter3D>>;
 
   using YODA::Counter;
   using YODA::Histo1D;
   using YODA::HistoBin1D;
   using YODA::Histo2D;
   using YODA::HistoBin2D;
   using YODA::Profile1D;
   using YODA::ProfileBin1D;
   using YODA::Profile2D;
   using YODA::ProfileBin2D;
   using YODA::Scatter1D;
   using YODA::Point1D;
   using YODA::Scatter2D;
   using YODA::Point2D;
   using YODA::Scatter3D;
   using YODA::Point3D;
 
   /// Function to get a map of all the refdata in a paper with the
   /// given @a papername.
   map<string, YODA::AnalysisObjectPtr> getRefData(const string& papername);
 
   /// @todo Also provide a Scatter3D getRefData() version?
 
   /// Get the file system path to the reference file for this paper.
   string getDatafilePath(const string& papername);
 
 
   /// Traits class to access the type of the AnalysisObject in the
   /// reference files.
   template<typename T> struct ReferenceTraits {};
   template<> struct ReferenceTraits<Counter> { typedef Counter RefT; };
   template<> struct ReferenceTraits<Scatter1D> { typedef Scatter1D RefT; };
   template<> struct ReferenceTraits<Histo1D> { typedef Scatter2D RefT; };
   template<> struct ReferenceTraits<Profile1D> { typedef Scatter2D RefT; };
   template<> struct ReferenceTraits<Scatter2D> { typedef Scatter2D RefT; };
   template<> struct ReferenceTraits<Histo2D> { typedef Scatter3D RefT; };
   template<> struct ReferenceTraits<Profile2D> { typedef Scatter3D RefT; };
   template<> struct ReferenceTraits<Scatter3D> { typedef Scatter3D RefT; };
 
 
   /// If @a dst and @a src both are of same subclass T, copy the
   /// contents of @a src into @a dst and return true. Otherwise return
   /// false.
   template <typename T>
   inline bool aocopy(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst) {
     shared_ptr<T> tsrc = dynamic_pointer_cast<T>(src);
     if ( !tsrc ) return false;
     shared_ptr<T> tdst = dynamic_pointer_cast<T>(dst);
     if ( !tdst ) return false;
     *tdst = *tsrc;
     return true;
   }
 
   /// If @a dst and @a src both are of same subclass T, add the
   /// contents of @a src into @a dst and return true. Otherwise return
   /// false.
   template <typename T>
   inline bool aoadd(YODA::AnalysisObjectPtr dst, YODA::AnalysisObjectPtr src, double scale) {
     shared_ptr<T> tsrc = dynamic_pointer_cast<T>(src);
     if ( !tsrc ) return false;
     shared_ptr<T> tdst = dynamic_pointer_cast<T>(dst);
     if ( !tdst ) return false;
     tsrc->scaleW(scale);
     *tdst += *tsrc;
     return true;
   }
 
   /// If @a dst is the same subclass as @a src, copy the contents of @a
   /// src into @a dst and return true. Otherwise return false.
   bool copyao(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst);
 
   /// If @a dst is the same subclass as @a src, scale the contents of
   /// @a src with @a scale and add it to @a dst and return true. Otherwise
   /// return false.
   bool addaos(YODA::AnalysisObjectPtr dst, YODA::AnalysisObjectPtr src, double scale);
 
   /// Check if two analysis objects have the same binning or, if not
   /// binned, are in other ways compatible.
   template <typename TPtr>
   inline bool bookingCompatible(TPtr a, TPtr b) {
     return a->sameBinning(*b);
   }
   inline bool bookingCompatible(CounterPtr, CounterPtr) {
     return true;
   }
   inline bool bookingCompatible(Scatter1DPtr a, Scatter1DPtr b) {
     return a->numPoints() == b->numPoints();
   }
   inline bool bookingCompatible(Scatter2DPtr a, Scatter2DPtr b) {
     return a->numPoints() == b->numPoints();
   }
   inline bool bookingCompatible(Scatter3DPtr a, Scatter3DPtr b) {
     return a->numPoints() == b->numPoints();
   }
 
 }
 
 #endif
diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc
--- a/src/Core/AnalysisHandler.cc
+++ b/src/Core/AnalysisHandler.cc
@@ -1,752 +1,739 @@
 // -*- 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"
 #include <regex>
 #include <iostream>
 
 using std::cout;
 using std::cerr;
 
 namespace {
 
   inline std::vector<std::string> split(const std::string& input, const std::string& regex) {
     // passing -1 as the submatch index parameter performs splitting
     std::regex re(regex);
     std::sregex_token_iterator
       first{input.begin(), input.end(), re, -1},
       last;
       return {first, last};
   }
 
 }
 
 
 namespace Rivet {
 
 
   AnalysisHandler::AnalysisHandler(const string& runname)
     : _runname(runname),
       // *** LEIF *** temporarily removed this
       // _eventcounter("/_EVTCOUNT"),
       // _xs(NAN), _xserr(NAN),
       _initialised(false), _ignoreBeams(false),
       _defaultWeightIdx(0), _dumpPeriod(0), _dumping(false)
   {  }
 
 
   AnalysisHandler::~AnalysisHandler()
   {  }
 
 
   Log& AnalysisHandler::getLog() const {
     return Log::getLog("Rivet.Analysis.Handler");
   }
 
 
   /// http://stackoverflow.com/questions/4654636/how-to-determine-if-a-string-is-a-number-with-c
   namespace {
     bool is_number(const std::string& s) {
       std::string::const_iterator it = s.begin();
       while (it != s.end() && std::isdigit(*it)) ++it;
       return !s.empty() && it == s.end();
     }
   }
 
   /// Check if any of the weightnames is not a number
   bool AnalysisHandler::haveNamedWeights() const {
     bool dec=false;
     for (unsigned int i=0;i<_weightNames.size();++i) {
       string s = _weightNames[i];
       if (!is_number(s)) {
         dec=true;
         break;
       }
     }
     return dec;
   }
 
 
   void AnalysisHandler::init(const GenEvent& ge) {
     if (_initialised)
       throw UserError("AnalysisHandler::init has already been called: cannot re-initialize!");
 
     /// @todo Should the Rivet analysis objects know about weight names?
 
     setRunBeams(Rivet::beams(ge));
     MSG_DEBUG("Initialising the analysis handler");
     _eventNumber = ge.event_number();
 
     setWeightNames(ge);
     if (haveNamedWeights())
         MSG_INFO("Using named weights");
     else
         MSG_INFO("NOT using named weights. Using first weight as nominal weight");
 
     _eventCounter = CounterPtr(weightNames(), Counter("_EVTCOUNT"));
 
     // Set the cross section based on what is reported by this event.
     if (ge.cross_section()) {
       MSG_TRACE("Getting cross section.");
       double xs = ge.cross_section()->cross_section();
       double xserr = ge.cross_section()->cross_section_error();
       setCrossSection(xs, xserr);
     }
 
     // Check that analyses are beam-compatible, and remove those that aren't
     const size_t num_anas_requested = analysisNames().size();
     vector<string> anamestodelete;
     for (const AnaHandle a : analyses()) {
       if (!_ignoreBeams && !a->isCompatible(beams())) {
         //MSG_DEBUG(a->name() << " requires beams " << a->requiredBeams() << " @ " << a->requiredEnergies() << " GeV");
         anamestodelete.push_back(a->name());
       }
     }
     for (const string& aname : anamestodelete) {
       MSG_WARNING("Analysis '" << aname << "' is incompatible with the provided beams: removing");
       removeAnalysis(aname);
     }
     if (num_anas_requested > 0 && analysisNames().empty()) {
       cerr << "All analyses were incompatible with the first event's beams\n"
            << "Exiting, since this probably wasn't intentional!" << endl;
       exit(1);
     }
 
     // Warn if any analysis' status is not unblemished
     for (const AnaHandle a : analyses()) {
       if ( a->info().preliminary() ) {
         MSG_WARNING("Analysis '" << a->name() << "' is preliminary: be careful, it may change and/or be renamed!");
       } else if ( a->info().obsolete() ) {
         MSG_WARNING("Analysis '" << a->name() << "' is obsolete: please update!");
       } else if (( a->info().unvalidated() ) ) {
         MSG_WARNING("Analysis '" << a->name() << "' is unvalidated: be careful, it may be broken!");
       }
     }
 
     // Initialize the remaining analyses
     _stage = Stage::INIT;
     for (AnaHandle a : analyses()) {
       MSG_DEBUG("Initialising analysis: " << a->name());
       try {
         // Allow projection registration in the init phase onwards
         a->_allowProjReg = true;
         a->init();
         //MSG_DEBUG("Checking consistency of analysis: " << a->name());
         //a->checkConsistency();
       } catch (const Error& err) {
         cerr << "Error in " << a->name() << "::init method: " << err.what() << endl;
         exit(1);
       }
       MSG_DEBUG("Done initialising analysis: " << a->name());
     }
     _stage = Stage::OTHER;
     _initialised = true;
     MSG_DEBUG("Analysis handler initialised");
   }
 
 // *** LEIF *** Reinstated this.
   void AnalysisHandler::setWeightNames(const GenEvent& ge) {
     /// reroute the print output to a std::stringstream and process
     /// The iteration is done over a map in hepmc2 so this is safe
     std::ostringstream stream;
     ge.weights().print(stream);  // Super lame, I know
     string str =  stream.str();
 
     std::regex re("(([^()]+))"); // Regex for stuff enclosed by parentheses ()
     size_t idx = 0;
     for (std::sregex_iterator i = std::sregex_iterator(str.begin(), str.end(), re); i != std::sregex_iterator(); ++i ) {
       std::smatch m = *i;
       vector<string> temp = ::split(m.str(), "[,]");
       if (temp.size() ==2) {
         MSG_DEBUG("Name of weight #" << _weightNames.size() << ": " << temp[0]);
 
         // store the default weight based on weight names
         if (temp[0] == "Weight" || temp[0] == "0" || temp[0] == "Default") {
           MSG_DEBUG(_weightNames.size() << " is being used as the nominal.");
           _weightNames.push_back("");
           _defaultWeightIdx = idx;
         } else
           _weightNames.push_back(temp[0]);
 
 
         idx++;
       }
     }
   }
 
 
 
   void AnalysisHandler::analyze(const GenEvent& ge) {
     // Call init with event as template if not already initialised
     if (!_initialised) init(ge);
     assert(_initialised);
 
     // Ensure that beam details match those from the first event (if we're checking beams)
     if ( !_ignoreBeams ) {
       const PdgIdPair beams = Rivet::beamIds(ge);
       const double sqrts = Rivet::sqrtS(ge);
       if (!compatible(beams, _beams) || !fuzzyEquals(sqrts, sqrtS())) {
         cerr << "Event beams mismatch: "
              << PID::toBeamsString(beams) << " @ " << sqrts/GeV << " GeV" << " vs. first beams "
              << this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << endl;
         exit(1);
       }
     }
 
     // Create the Rivet event wrapper
     /// @todo Filter/normalize the event here
     bool strip = ( getEnvParam("RIVET_STRIP_HEPMC", string("NOOOO") ) != "NOOOO" );
     Event event(ge, strip);
 
     // Weights
     /// @todo Drop this / just report first weight when we support multiweight events
     // *** LEIF *** temporarily removed this
     // _eventcounter.fill(event.weight());
     // MSG_DEBUG("Event #" << _eventcounter.numEntries() << " weight = " << event.weight());
 
     // *** LEIF *** temporarily removed this
     // // Cross-section
     // #if defined ENABLE_HEPMC_3
     // if (ge.cross_section()) {
     //   //@todo HepMC3::GenCrossSection methods aren't const accessible :(
     //   RivetHepMC::GenCrossSection gcs = *(event.genEvent()->cross_section());
     //   _xs = gcs.xsec();
     //   _xserr = gcs.xsec_err();
     // }
     // #elif defined HEPMC_HAS_CROSS_SECTION
     // if (ge.cross_section()) {
     //   _xs = ge.cross_section()->cross_section();
     //   _xserr = ge.cross_section()->cross_section_error();
     // }
     // #endif
 
     // Won't happen for first event because _eventNumber is set in init()
     if (_eventNumber != ge.event_number()) {
         /// @todo Can we get away with not passing a matrix?
 
         MSG_TRACE("AnalysisHandler::analyze(): Pushing _eventCounter to persistent.");
         _eventCounter.get()->pushToPersistent(_subEventWeights);
         // if this is indeed a new event, push the temporary
         // histograms and reset
         for (const AnaHandle& a : analyses()) {
             for (auto ao : a->analysisObjects()) {
                 MSG_TRACE("AnalysisHandler::analyze(): Pushing " << a->name() << "'s " << ao->name() << " to persistent.");
                 ao.get()->pushToPersistent(_subEventWeights);
             }
             MSG_TRACE("AnalysisHandler::analyze(): finished pushing " << a->name() << "'s objects to persistent.");
         }
 
         _eventNumber = ge.event_number();
 
         MSG_DEBUG("nominal event # " << _eventCounter.get()->_persistent[0]->numEntries());
         MSG_DEBUG("nominal sum of weights: " << _eventCounter.get()->_persistent[0]->sumW());
         MSG_DEBUG("Event has " << _subEventWeights.size() << " sub events.");
         _subEventWeights.clear();
     }
 
 
     MSG_TRACE("starting new sub event");
     _eventCounter.get()->newSubEvent();
 
     for (const AnaHandle& a : analyses()) {
         for (auto ao : a->analysisObjects()) {
             ao.get()->newSubEvent();
         }
     }
 
     _subEventWeights.push_back(event.weights());
     MSG_DEBUG("Analyzing subevent #" << _subEventWeights.size() - 1 << ".");
 
     _eventCounter->fill();
     // Run the analyses
     for (AnaHandle a : analyses()) {
       MSG_TRACE("About to run analysis " << a->name());
       try {
         a->analyze(event);
       } catch (const Error& err) {
         cerr << "Error in " << a->name() << "::analyze method: " << err.what() << endl;
         exit(1);
       }
       MSG_TRACE("Finished running analysis " << a->name());
     }
 
     if ( _dumpPeriod > 0 && numEvents() > 0 && numEvents()%_dumpPeriod == 0 ) {
       MSG_INFO("Dumping intermediate results to " << _dumpFile << ".");
       _dumping = numEvents()/_dumpPeriod;
       finalize();
       writeData(_dumpFile);
       _dumping = 0;
     }
 
   }
 
 
   void AnalysisHandler::analyze(const GenEvent* ge) {
     if (ge == nullptr) {
       MSG_ERROR("AnalysisHandler received null pointer to GenEvent");
       //throw Error("AnalysisHandler received null pointer to GenEvent");
     }
     analyze(*ge);
   }
 
 
   void AnalysisHandler::finalize() {
     if (!_initialised) return;
     MSG_INFO("Finalising analyses");
 
-    // First push all analyses' objects to persistent
+    // First push all analyses' objects to persistent and final
     MSG_TRACE("AnalysisHandler::finalize(): Pushing analysis objects to persistent.");
     _eventCounter.get()->pushToPersistent(_subEventWeights);
+    _eventCounter.get()->pushToFinal();
+    _xs.get()->pushToFinal();
     for (const AnaHandle& a : analyses()) {
-      for (auto ao : a->analysisObjects())
+      for (auto ao : a->analysisObjects()) {
         ao.get()->pushToPersistent(_subEventWeights);
-    }
-
-    // First we make copies of all analysis objects.
-    map<string,YODA::AnalysisObjectPtr> backupAOs;
-    for (auto ao : getData(false, true, false) )
-      backupAOs[ao->path()] = YODA::AnalysisObjectPtr(ao->newclone());
-
-    // Now we run the (re-entrant) finalize() functions for all analyses.
-    MSG_INFO("Finalising analyses");
-    for (AnaHandle a : analyses()) {
-    // *** LEIF *** temporarily removed this
-      // a->setCrossSection(_xs);
-      try {
-        if ( !_dumping || a->info().reentrant() )  a->finalize();
-        else if ( _dumping == 1 )
-          MSG_INFO("Skipping finalize in 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);
+        ao.get()->pushToFinal();
       }
     }
 
-    // Now we copy all analysis objects to the list of finalized
-    // ones, and restore the value to their original ones.
-    _finalizedAOs.clear();
-    for ( auto ao : getYodaAOs(false, false, false) )
-      _finalizedAOs.push_back(YODA::AnalysisObjectPtr(ao->newclone()));
-    for ( auto ao : getYodaAOs(false, true, false) ) {
-      // TODO: This should be possible to do in a nicer way, with a flag etc.
-      if (ao->path().find("/FINAL") != std::string::npos) continue;
-      auto aoit = backupAOs.find(ao->path());
-      if ( aoit == backupAOs.end() ) {
-        AnaHandle ana = analysis(split(ao->path(), "/")[0]);
-        if ( ana ) ana->removeAnalysisObject(ao->path());
-      } else
-        copyao(aoit->second, ao);
+    for (AnaHandle a : analyses()) {
+      if ( _dumping && !a->info().reentrant() )  {
+        if ( _dumping == 1 )
+          MSG_INFO("Skipping finalize in periodic dump of " << a->name()
+                   << " as it is not declared reentrant.");
+        continue;
+      }
+      for (size_t iW = 0; iW < numWeights(); iW++) {
+        _eventCounter.get()->setActiveFinalWeightIdx(iW);
+        _xs.get()->setActiveFinalWeightIdx(iW);
+        for (auto ao : a->analysisObjects())
+          ao.get()->setActiveFinalWeightIdx(iW);
+        try {
+          MSG_TRACE("running " << a->name() << "::finalize() for weight " << iW << ".");
+          a->finalize();
+        } catch (const Error& err) {
+          cerr << "Error in " << a->name() << "::finalize method: " << err.what() << '\n';
+          exit(1);
+        }
+      }
     }
 
     // Print out number of events processed
     const int nevts = numEvents();
     MSG_INFO("Processed " << nevts << " event" << (nevts != 1 ? "s" : ""));
 
-    // // Delete analyses
-    // MSG_DEBUG("Deleting analyses");
-    // _analyses.clear();
-
+    if ( _dumping ) return;
+ 
     // 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[analysisname] = 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) {
     MSG_DEBUG("Removing analysis '" << analysisname << "'");
     if (_analyses.find(analysisname) != _analyses.end()) _analyses.erase(analysisname);
     // }
     return *this;
   }
 
 
   void AnalysisHandler::addData(const std::vector<YODA::AnalysisObjectPtr>& aos) {
     for (const YODA::AnalysisObjectPtr ao : aos) {
       string path = ao->path();
       if ( path.substr(0, 5) != "/RAW/" ) {
         _orphanedPreloads.push_back(ao);
         continue;
       }
 
       path = path.substr(4);
       ao->setPath(path);
       if (path.size() > 1) { // path > "/"
         try {
           const string ananame =  ::split(path, "/")[0];
           AnaHandle a = analysis(ananame);
           /// @todo FIXXXXX
           //MultiweightAOPtr mao = ????; /// @todo generate right Multiweight object from ao
           //a->addAnalysisObject(mao); /// @todo Need to statistically merge...
         } catch (const Error& e) {
           MSG_TRACE("Adding analysis object " << path <<
                     " to the list of orphans.");
           _orphanedPreloads.push_back(ao);
         }
       }
     }
   }
 
 
   void AnalysisHandler::stripOptions(YODA::AnalysisObjectPtr ao,
                                      const vector<string> & delopts) const {
     string path = ao->path();
     string ananame = split(path, "/")[0];
     vector<string> anaopts = split(ananame, ":");
     for ( int i = 1, N = anaopts.size(); i < N; ++i )
       for ( auto opt : delopts )
         if ( opt == "*" || anaopts[i].find(opt + "=") == 0 )
           path.replace(path.find(":" + anaopts[i]), (":" + anaopts[i]).length(), "");
     ao->setPath(path);
   }
 
 
   void AnalysisHandler::mergeYodas(const vector<string> & ,
                                    const vector<string> & , bool ) {}
   // void AnalysisHandler::mergeYodas(const vector<string> & aofiles,
   //                                  const vector<string> & delopts, bool equiv) {
   //   vector< vector<YODA::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 analyses..
   //   for (const string& file : aofiles ) {
   //     Scatter1DPtr xsec;
   //     CounterPtr sow;
 
   //     // For each file make sure that cross section and sum-of-weights
   //     // objects are present and stor all RAW ones in a vector;
   //     vector<YODA::AnalysisObjectPtr> aos;
   //     try {
   //       /// @todo Use new YODA SFINAE to fill the smart ptr vector directly
   //       vector<YODA::AnalysisObject*> aos_raw;
   //       YODA::read(file, aos_raw);
   //       for (YODA::AnalysisObject* aor : aos_raw) {
   //         YODA::AnalysisObjectPtr ao(aor);
   //         if ( ao->path().substr(0, 5) != "/RAW/" ) continue;
   //         ao->setPath(ao->path().substr(4));
   //         if ( ao->path() == "/_XSEC" )
   //           xsec = dynamic_pointer_cast<Scatter1D>(ao);
   //         else if ( ao->path() == "/_EVTCOUNT" )
   //           sow = dynamic_pointer_cast<Counter>(ao);
   //         else {
   //           stripOptions(ao, delopts);
   //           string ananame = split(ao->path(), "/")[0];
   //           if ( ananames.insert(ananame).second ) addAnalysis(ananame);
   //           aos.push_back(ao);
   //         }
   //       }
   //       if ( !xsec || !sow ) {
   //         MSG_ERROR( "Error in AnalysisHandler::mergeYodas: The file " << file
   //                    << " did not contain weights and cross section info.");
   //         exit(1);
   //       }
   //       xsecs.push_back(xsec->point(0).x());
   //       sows.push_back(sow);
   //       xsecerrs.push_back(sqr(xsec->point(0).xErrAvg()));
   //       _eventCounter->operator+=(*sow); //< HAHAHAHAHA!
   //       sows.push_back(sow);
   //       aosv.push_back(aos);
   //     } catch (...) { //< YODA::ReadError&
   //       throw UserError("Unexpected error in reading file: " + file);
   //     }
   //   }
 
   //   // Now calculate the scale to be applied for all bins in a file
   //   // and get the common cross section and sum of weights.
   //   _xs = _xserr = 0.0;
   //   for ( int i = 0, N = sows.size(); i < N; ++i ) {
   //     double effnent = sows[i]->effNumEntries();
   //     _xs += (equiv? effnent: 1.0)*xsecs[i];
   //     _xserr += (equiv? sqr(effnent): 1.0)*xsecerrs[i];
   //   }
 
   //   vector<double> scales(sows.size(), 1.0);
   //   if ( equiv ) {
   //     _xs /= _eventCounter.effNumEntries();
   //     _xserr = sqrt(_xserr)/_eventCounter.effNumEntries();
   //   } 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<YODA::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 (YODA::AnalysisObject* aor : aos_raw) aos.push_back(YODA::AnalysisObjectPtr(aor));
     //} catch (const YODA::ReadError & e) {
     } catch (...) { //< YODA::ReadError&
       throw UserError("Unexpected error in reading file: " + filename);
     }
     if (!aos.empty()) addData(aos);
   }
 
 
   vector<MultiweightAOPtr> AnalysisHandler::getRivetAOs() const {
       vector<MultiweightAOPtr> rtn;
 
       for (AnaHandle a : analyses()) {
           for (const auto & ao : a->analysisObjects()) {
               rtn.push_back(ao);
           }
       }
 
       rtn.push_back(_eventCounter);
       rtn.push_back(_xs);
 
       return rtn;
   }
 
   vector<YODA::AnalysisObjectPtr> AnalysisHandler::getYodaAOs(bool includeorphans,
                                                               bool includetmps,
                                                               bool usefinalized) const {
       vector<YODA::AnalysisObjectPtr> rtn;
       if (usefinalized)
         rtn = _finalizedAOs;
       else {
         for (auto rao : getRivetAOs()) {
           // need to set the index
           // before we can search the PATH
           rao.get()->setActiveWeightIdx(_defaultWeightIdx);
           // Exclude paths from final write-out if they contain a "TMP" layer (i.e. matching "/TMP/")
           if (!includetmps && rao->path().find("/TMP/") != string::npos)
             continue;
 
           for (size_t iW = 0; iW < numWeights(); iW++) {
             rao.get()->setActiveWeightIdx(iW);
             rtn.push_back(rao.get()->activeYODAPtr());
           }
         }
       }
 
       // Sort histograms alphanumerically by path before write-out
       sort(rtn.begin(), rtn.end(),
            [](YODA::AnalysisObjectPtr a, YODA::AnalysisObjectPtr b) {
                 return a->path() < b->path();
             }
           );
 
       return rtn;
   }
 
 
   vector<YODA::AnalysisObjectPtr> AnalysisHandler::getData(bool includeorphans,
                                                            bool includetmps,
                                                            bool usefinalized) const {
     return getYodaAOs(includeorphans, includetmps, usefinalized);
   }
 
 
   void AnalysisHandler::writeData(const string& filename) const {
     vector<YODA::AnalysisObjectPtr> out = _finalizedAOs;
     set<string> finalana;
     for ( auto ao : out) finalana.insert(ao->path());
     out.reserve(2*out.size());
     vector<YODA::AnalysisObjectPtr> aos = getData(false, true);
     for ( auto ao : aos ) {
       ao = YODA::AnalysisObjectPtr(ao->newclone());
       ao->setPath("/RAW" + ao->path());
       out.push_back(ao);
     }
 
     try {
       YODA::write(filename, aos.begin(), aos.end());
     } catch (...) { //< YODA::WriteError&
       throw UserError("Unexpected error in writing file: " + filename);
     }
   }
 
 
   string AnalysisHandler::runName() const { return _runname; }
   size_t AnalysisHandler::numEvents() const { return _eventCounter->numEntries(); }
 
 
   std::vector<std::string> AnalysisHandler::analysisNames() const {
     std::vector<std::string> rtn;
     for (AnaHandle a : analyses()) {
       rtn.push_back(a->name());
     }
     return rtn;
   }
 
 
   AnalysisHandler& AnalysisHandler::addAnalyses(const std::vector<std::string>& analysisnames) {
     for (const string& aname : analysisnames) {
       //MSG_DEBUG("Adding analysis '" << aname << "'");
       addAnalysis(aname);
     }
     return *this;
   }
 
 
   AnalysisHandler& AnalysisHandler::removeAnalyses(const std::vector<std::string>& analysisnames) {
     for (const string& aname : analysisnames) {
       removeAnalysis(aname);
     }
     return *this;
   }
 
 
   AnalysisHandler& AnalysisHandler::setCrossSection(double xs, double xserr) {
     _xs = Scatter1DPtr(weightNames(), Scatter1D("_XSEC"));
     _eventCounter.get()->setActiveWeightIdx(_defaultWeightIdx);
     double nomwgt = sumW();
 
     // The cross section of each weight variation is the nominal cross section
     // times the sumW(variation) / sumW(nominal).
     // This way the cross section will work correctly
     for (size_t iW = 0; iW < numWeights(); iW++) {
       _eventCounter.get()->setActiveWeightIdx(iW);
       double s = sumW() / nomwgt;
       _xs.get()->setActiveWeightIdx(iW);
       _xs->addPoint(xs*s, xserr*s);
     }
 
     _eventCounter.get()->unsetActiveWeight();
     _xs.get()->unsetActiveWeight();
     return *this;
   }
 
 
   // *** LEIF *** temporarily removed this
   // bool AnalysisHandler::hasCrossSection() const {
   //   return (!std::isnan(crossSection()));
   // }
 
 
   AnalysisHandler& AnalysisHandler::addAnalysis(Analysis* analysis) {
     analysis->_analysishandler = this;
     // _analyses.insert(AnaHandle(analysis));
     _analyses[analysis->name()] = 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/RivetYODA.cc b/src/Tools/RivetYODA.cc
--- a/src/Tools/RivetYODA.cc
+++ b/src/Tools/RivetYODA.cc
@@ -1,428 +1,428 @@
 #include "Rivet/Config/RivetCommon.hh"
 #include "Rivet/Tools/RivetYODA.hh"
 #include "Rivet/Tools/RivetPaths.hh"
 #include "YODA/ReaderYODA.h"
 #include "YODA/ReaderAIDA.h"
 
 // use execinfo for backtrace if available
 #include "DummyConfig.hh"
 #ifdef HAVE_EXECINFO_H
 #include <execinfo.h>
 #endif
 
 using namespace std;
 
 namespace Rivet {
 
 
 template <class T>
 Wrapper<T>::~Wrapper() {}
 
 
 template <class T>
 Wrapper<T>::Wrapper(const vector<string>& weightNames, const T & p)
 {
   for (const string& weightname : weightNames) {
     _persistent.push_back(make_shared<T>(p));
     _final.push_back(make_shared<T>(p));
 
     auto obj = _persistent.back();
     auto final = _final.back();
     if (weightname != "") {
       obj->setPath("/RAW" + obj->path() + "[" + weightname + "]");
-      final->setPath(obj->path() + "[" + weightname + "]");
+      final->setPath(final->path() + "[" + weightname + "]");
     }
   }
 }
 
 
 template <class T>
 typename T::Ptr Wrapper<T>::active() const {
     if ( !_active ) {
 #ifdef HAVE_BACKTRACE
       void * buffer[4];
       backtrace(buffer, 4);
       backtrace_symbols_fd(buffer, 4 , 1);
 #endif
       assert(false && "No active pointer set. Was this object booked in init()?");
     }
     return _active;
 }
 
 
 template <class T>
 void Wrapper<T>::newSubEvent() {
   typename TupleWrapper<T>::Ptr tmp
       = make_shared<TupleWrapper<T>>(_persistent[0]->clone());
   tmp->reset();
   _evgroup.push_back( tmp );
   _active = _evgroup.back();
   assert(_active);
 }
 
 
 
   string getDatafilePath(const string& papername) {
     /// Try to find YODA otherwise fall back to try AIDA
     const string path1 = findAnalysisRefFile(papername + ".yoda");
     if (!path1.empty()) return path1;
     const string path2 = findAnalysisRefFile(papername + ".aida");
     if (!path2.empty()) return path2;
     throw Rivet::Error("Couldn't find ref data file '" + papername + ".yoda" +
                        " in data path, '" + getRivetDataPath() + "', or '.'");
   }
 
 
   map<string, YODA::AnalysisObjectPtr> getRefData(const string& papername) {
     const string datafile = getDatafilePath(papername);
 
     // Make an appropriate data file reader and read the data objects
     /// @todo Remove AIDA support some day...
     YODA::Reader& reader = (datafile.find(".yoda") != string::npos) ?   \
       YODA::ReaderYODA::create() : YODA::ReaderAIDA::create();
     vector<YODA::AnalysisObject *> aovec;
     reader.read(datafile, aovec);
 
     // Return value, to be populated
     map<string, YODA::AnalysisObjectPtr> rtn;
     for ( YODA::AnalysisObject* ao : aovec ) {
       YODA::AnalysisObjectPtr refdata(ao);
       if (!refdata) continue;
       const string plotpath = refdata->path();
       // Split path at "/" and only return the last field, i.e. the histogram ID
       const size_t slashpos = plotpath.rfind("/");
       const string plotname = (slashpos+1 < plotpath.size()) ? plotpath.substr(slashpos+1) : "";
       rtn[plotname] = refdata;
     }
     return rtn;
   }
 
 
 }
 
 
 namespace {
 
   using Rivet::Fill;
   using Rivet::Fills;
   using Rivet::TupleWrapper;
 
     template <class T>
     double get_window_size(const typename T::Ptr & histo,
                            typename T::BinType x) {
         // the bin index we fall in
         const auto binidx = histo->binIndexAt(x);
         // gaps, overflow, underflow don't contribute
         if ( binidx == -1 )
             return 0;
 
 
         const auto & b = histo->bin(binidx);
 
         // if we don't have a valid neighbouring bin,
         // we use infinite width
         typename T::Bin b1(-1.0/0.0, 1.0/0.0);
 
         // points in the top half compare to the upper neighbour
         if ( x > b.xMid() ) {
             size_t nextidx = binidx + 1;
             if ( nextidx < histo->bins().size() )
                 b1 = histo->bin(nextidx);
         }
         else { // compare to the lower neighbour
             int nextidx = binidx - 1;
             if ( nextidx >= 0 )
                 b1 = histo->bin(nextidx);
         }
         // the factor 2 is arbitrary, could poss. be smaller
         return min( b.width(), b1.width() ) / 2.0;
     }
 
     template <class T>
      typename T::BinType
      fillT2binT(typename T::FillType a) {
        return a;
     }
 
     template <>
     YODA::Profile1D::BinType
     fillT2binT<YODA::Profile1D>(YODA::Profile1D::FillType a) {
       return get<0>(a);
     }
 
     template <>
     YODA::Profile2D::BinType
     fillT2binT<YODA::Profile2D>(YODA::Profile2D::FillType a) {
       return YODA::Profile2D::BinType{ get<0>(a), get<1>(a) };
     }
 
 
 
   template <class T>
   void commit(vector<typename T::Ptr> & persistent,
               const vector< vector<Fill<T>> > & tuple,
               const vector<valarray<double>> & weights ) {
 
     // TODO check if all the xs are in the same bin anyway!
     // Then no windowing needed
 
     assert(persistent.size() == weights[0].size());
 
     for ( const auto & x : tuple ) {
       double maxwindow = 0.0;
       for ( const auto & xi : x ) {
         // TODO check for NOFILL here
         // persistent[0] has the same binning as all the other persistent objects
         double window = get_window_size<T>(persistent[0], fillT2binT<T>(xi.first));
         if ( window > maxwindow )
           maxwindow = window;
       }
       const double wsize = maxwindow;
       // all windows have same size
 
       set<double> edgeset;
       // bin edges need to be in here!
       for ( const auto & xi : x ) {
         edgeset.insert(fillT2binT<T>(xi.first) - wsize);
         edgeset.insert(fillT2binT<T>(xi.first) + wsize);
       }
 
           vector< std::tuple<double,valarray<double>,double> > hfill;
           double sumf = 0.0;
           auto edgit = edgeset.begin();
           double ehi = *edgit;
           while ( ++edgit != edgeset.end() ) {
             double elo = ehi;
             ehi = *edgit;
             valarray<double> sumw(0.0, persistent.size()); // need m copies of this
             bool gap = true; // Check for gaps between the sub-windows.
             for ( size_t i = 0; i < x.size(); ++i  ) {
               // check equals comparisons here!
               if ( fillT2binT<T>(x[i].first) + wsize >= ehi
                    &&
                    fillT2binT<T>(x[i].first) - wsize <= elo ) {
                 sumw += x[i].second * weights[i];
                 gap = false;
               }
             }
             if ( gap ) continue;
             hfill.push_back( make_tuple( (ehi + elo)/2.0, sumw, (ehi - elo) ) );
             sumf += ehi - elo;
           }
 
           for ( auto f : hfill )
             for ( size_t m = 0; m < persistent.size(); ++m )
               persistent[m]->fill( get<0>(f), get<1>(f)[m], get<2>(f)/sumf );
             // Note the scaling to one single fill
 
     }
 
   }
 
 
   template<>
   void commit<YODA::Histo2D>(vector<YODA::Histo2D::Ptr> & persistent,
               const vector< vector<Fill<YODA::Histo2D>> > & tuple,
               const vector<valarray<double>> & weights)
   {}
 
   template<>
   void commit<YODA::Profile2D>(vector<YODA::Profile2D::Ptr> & persistent,
               const vector< vector<Fill<YODA::Profile2D>> > & tuple,
               const vector<valarray<double>> & weights)
   {}
 
     template <class T>
     double distance(T a, T b) {
       return abs(a - b);
     }
 
     template <>
     double distance<tuple<double,double> >(tuple<double,double> a, tuple<double,double> b) {
       return Rivet::sqr(get<0>(a) - get<0>(b)) + Rivet::sqr(get<1>(a) - get<1>(b));
     }
 
   }
 
 
   namespace Rivet {
 
     bool copyao(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst) {
       for (const std::string& a : src->annotations())
         dst->setAnnotation(a, src->annotation(a));
       if ( aocopy<Counter>(src,dst) ) return true;
       if ( aocopy<Histo1D>(src,dst) ) return true;
       if ( aocopy<Histo2D>(src,dst) ) return true;
       if ( aocopy<Profile1D>(src,dst) ) return true;
       if ( aocopy<Profile2D>(src,dst) ) return true;
       if ( aocopy<Scatter1D>(src,dst) ) return true;
       if ( aocopy<Scatter2D>(src,dst) ) return true;
       if ( aocopy<Scatter3D>(src,dst) ) return true;
       return false;
     }
 
     bool addaos(YODA::AnalysisObjectPtr dst, YODA::AnalysisObjectPtr src, double scale) {
       if ( aoadd<Counter>(dst,src,scale) ) return true;
       if ( aoadd<Histo1D>(dst,src,scale) ) return true;
       if ( aoadd<Histo2D>(dst,src,scale) ) return true;
       if ( aoadd<Profile1D>(dst,src,scale) ) return true;
       if ( aoadd<Profile2D>(dst,src,scale) ) return true;
       return false;
     }
 
   }
 
 
 namespace {
 
 /// fills is a vector of sub-event with an ordered set of x-values of
 /// the fills in each sub-event. NOFILL should be an "impossible"
 /// value for this histogram. Returns a vector of sub-events with
 /// an ordered vector of fills (including NOFILLs) for each sub-event.
 template <class T>
 vector< vector<Fill<T> > >
 match_fills(const vector<typename TupleWrapper<T>::Ptr> & evgroup, const Fill<T> & NOFILL)
 {
   vector< vector<Fill<T> > > matched;
   // First just copy subevents into vectors and find the longest vector.
   unsigned int maxfill = 0; // length of biggest vector
   int imax = 0; // index position of biggest vector
   for ( const auto & it : evgroup ) {
     const auto & subev = it->fills();
     if ( subev.size() > maxfill ) {
       maxfill = subev.size();
       imax = matched.size();
     }
     matched.push_back(vector<Fill<T> >(subev.begin(), subev.end()));
   }
   // Now, go through all subevents with missing fills.
   const vector<Fill<T>> & full = matched[imax]; // the longest one
   for ( auto & subev : matched ) {
     if ( subev.size() == maxfill ) continue;
 
     // Add NOFILLs to the end;
     while ( subev.size() < maxfill ) subev.push_back(NOFILL);
 
     // Iterate from the back and shift all fill values backwards by
     // swapping with NOFILLs so that they better match the full
     // subevent.
     for ( int i = maxfill - 1; i >= 0; --i ) {
       if ( subev[i] == NOFILL ) continue;
       size_t j = i;
       while ( j + 1 < maxfill && subev[j + 1] == NOFILL &&
               distance(fillT2binT<T>(subev[j].first),
                        fillT2binT<T>(full[j].first))
               >
               distance(fillT2binT<T>(subev[j].first),
                        fillT2binT<T>(full[j + 1].first)) )
       {
             swap(subev[j], subev[j + 1]);
             ++j;
       }
     }
   }
   // transpose
   vector<vector<Fill<T>>> result(maxfill,vector<Fill<T>>(matched.size()));
   for (size_t i = 0; i < matched.size(); ++i)
       for (size_t j = 0; j < maxfill; ++j)
           result.at(j).at(i) = matched.at(i).at(j);
   return result;
 }
 
 }
 
 
 
 
 namespace Rivet {
 
   template <class T>
   void Wrapper<T>::pushToPersistent(const vector<valarray<double> >& weight) {
       assert( _evgroup.size() == weight.size() );
 
       // have we had subevents at all?
       const bool have_subevents = _evgroup.size() > 1;
       if ( ! have_subevents ) {
 
 
 
         // simple replay of all tuple entries
         // each recorded fill is inserted into all persistent weightname histos
         for ( size_t m = 0; m < _persistent.size(); ++m )
             for ( const auto & f : _evgroup[0]->fills() )
                 _persistent[m]->fill( f.first, f.second * weight[0][m] );
 
 
 
       } else {
 
 
 
         // outer index is subevent, inner index is jets in the event
         vector<vector<Fill<T>>> linedUpXs
             = match_fills<T>(_evgroup, {typename T::FillType(), 0.0});
         commit<T>( _persistent, linedUpXs, weight );
 
 
 
       }
       _evgroup.clear();
       _active.reset();
   }
 
   template <class T>
-  void Wrapper<T>::pushToFinal(const vector<valarray<double> >& weight) {
+  void Wrapper<T>::pushToFinal() {
     for ( size_t m = 0; m < _persistent.size(); ++m ) {
       copyao(_persistent.at(m), _final.at(m));
     }
   }
 
 
 
   template <>
   void Wrapper<YODA::Counter>::pushToPersistent(const vector<valarray<double> >& weight) {
     for ( size_t m = 0; m < _persistent.size(); ++m ) {
       for ( size_t n = 0; n < _evgroup.size(); ++n ) {
         for ( const auto & f : _evgroup[n]->fills() ) {
           _persistent[m]->fill( f.second * weight[n][m] );
         }
       }
     }
 
     _evgroup.clear();
     _active.reset();
   }
 
   template <>
   void Wrapper<YODA::Scatter1D>::pushToPersistent(const vector<valarray<double> >& weight) {
 
     _evgroup.clear();
     _active.reset();
   }
 
   template <>
   void Wrapper<YODA::Scatter2D>::pushToPersistent(const vector<valarray<double> >& weight) {
 
     _evgroup.clear();
     _active.reset();
   }
 
   template <>
   void Wrapper<YODA::Scatter3D>::pushToPersistent(const vector<valarray<double> >& weight) {
 
     _evgroup.clear();
     _active.reset();
   }
 
   // explicitly instantiate all wrappers
 
   template class Wrapper<YODA::Histo1D>;
   template class Wrapper<YODA::Histo2D>;
   template class Wrapper<YODA::Profile1D>;
   template class Wrapper<YODA::Profile2D>;
   template class Wrapper<YODA::Counter>;
   template class Wrapper<YODA::Scatter1D>;
   template class Wrapper<YODA::Scatter2D>;
   template class Wrapper<YODA::Scatter3D>;
 
 }