diff --git a/examples/AnalysisPrint.cc b/examples/AnalysisPrint.cc
index 522fe2b..02dbc56 100644
--- a/examples/AnalysisPrint.cc
+++ b/examples/AnalysisPrint.cc
@@ -1,78 +1,80 @@
 //! HEJ analysis to output the cross section to a file
 
 #include <cmath>
 #include <fstream>
 #include <iostream>
 #include <memory>
 #include <string>
 
 #include "HEJ/Analysis.hh"
 #include "HEJ/Event.hh"
 
 #include "yaml-cpp/yaml.h"
 
 #include "LHEF/LHEF.h"
 
 namespace {
   class AnalysisPrint: public HEJ::Analysis {
   public:
     AnalysisPrint(YAML::Node const & config, LHEF::HEPRUP const & heprup):
-      xsection_{0.}, xsection_error_{0.},
+      wt_sum_{0.}, wt2_sum_{0.},
       outfile_{config["output"].as<std::string>()},
       generators_{heprup.generators}
     {}
 
     void fill(
       HEJ::Event const & event,
       HEJ::Event const & /* FO_event */
     ) override {
       const double wt = event.central().weight;
-      xsection_ += wt;
-      xsection_error_ += wt*wt; // this error estimate is too small
+      wt_sum_ += wt;
+      wt2_sum_ += wt*wt; // this error estimate is too small
     }
 
     bool pass_cuts(
       HEJ::Event const & /* event */,
       HEJ::Event const & /* FO_event */
     ) override {
       return true;
     }
 
-    void scale(double xsscale) override {
-      // used to scale the cross sections and histograms before .finalise in case the
-      // normalisation is unknown until the end of the run
-      xsection_*=xsscale;
-      xsection_error_*=xsscale*xsscale;
+    void set_xs_scale(double xsscale) override {
+      xs_scale_ = xsscale;
     }
 
     void finalise() override {
+      // compute cross section from sum of weights and scaling factor
+      const double xsection = wt_sum_ * xs_scale_;
+      const double xsection_error = std::sqrt(wt2_sum_) * xs_scale_;
+
       // print to screen
       std::cout << "Generated with:\n";
       for(auto const & generator: generators_)
         std::cout << generator.name << " " << generator.version << "\n";
-      std::cout << "cross section: " << xsection_ << " +- "
-        << std::sqrt(xsection_error_) << std::endl;
+      std::cout << "cross section: " << xsection << " +- "
+        << xsection_error << std::endl;
 
       // print to file
       std::ofstream fout{outfile_};
       fout << "Generated with\n";
       for(auto const & generator: generators_)
         fout << generator.name << " " << generator.version << "\n";
-      fout << "cross section: " << xsection_ << " +- "
-        << std::sqrt(xsection_error_) << std::endl;
+      fout << "cross section: " << xsection << " +- "
+        << xsection_error << std::endl;
     }
 
     private:
-    double xsection_, xsection_error_;
+    double wt_sum_, wt2_sum_;
+    double xs_scale_{1.0};
     std::string outfile_;
     std::vector<LHEF::Generator> generators_;
   };
 }
 
 extern "C"
 __attribute__((visibility("default")))
 std::unique_ptr<HEJ::Analysis> make_analysis(
     YAML::Node const & config, LHEF::HEPRUP const & heprup
 ){
   return std::make_unique<AnalysisPrint>(config, heprup);
 }
diff --git a/examples/AnalysisROOT.cc b/examples/AnalysisROOT.cc
index 2419823..f73ab03 100644
--- a/examples/AnalysisROOT.cc
+++ b/examples/AnalysisROOT.cc
@@ -1,112 +1,118 @@
 //! Simple HEJ analysis to plot Delta y_{fb} with ROOT
 
 #include <memory>     // for std::unique_ptr
 #include <iostream>   // for std::cout
 #include <stdexcept>
 
 #include "TFile.h"
 #include "TKey.h"
 #include "TH1.h"
 
 #include "HEJ/Analysis.hh"
 #include "HEJ/Event.hh"
 #include "HEJ/Version.hh"
 
 namespace YAML {
   class Node;
 }
 
 namespace LHEF {
   class HEPRUP;
 }
 namespace {
   void Fill2DHistogram(TH1D & h,double x,double wt)
   {
     const int binnumber=h.GetXaxis()->FindBin(x);
     const double binwidth=h.GetXaxis()->GetBinWidth(binnumber);
     h.Fill(x,wt/binwidth);
   }
 
   class AnalysisROOT: public HEJ::Analysis {
   private:
+    double xs_scale_{1.0};
     TFile hfile_;
     TH1D htotal_,hydif_,hyf_,hyb_;
   public:
     AnalysisROOT(
       YAML::Node const & /* config */, LHEF::HEPRUP const & /* heprup */
     ):hfile_{"HEJ.root","RECREATE",
              ("Output from "+HEJ::Version::package_name_full()).c_str()},
       htotal_{"total","total",1,-0.5,.5},
       hydif_{"ydif","ydif",40,0.,10.},
       hyf_{"yf","yf",40,-5.,5.},
       hyb_{"yb","yb",40,-5.,5.}
     {
       // Setup histograms
       htotal_.Sumw2();
       hydif_.Sumw2();
       hyf_.Sumw2();
       hyb_.Sumw2();
     }
 
     void fill(
       HEJ::Event const & event,
       HEJ::Event const & /* FO_event */
     ) override {
       const double weight=event.central().weight;
       htotal_.Fill(0.0,weight);
       auto const & jets=event.jets(); // jets are sorted in rapidity
       Fill2DHistogram(hydif_,
           fabs(jets.back().rapidity()-jets.front().rapidity()), weight);
       Fill2DHistogram(hyf_, jets.front().rapidity(), weight);
       Fill2DHistogram(hyb_, jets.back().rapidity(), weight);
     }
 
     bool pass_cuts(
       HEJ::Event const & event,
       HEJ::Event const & /* FO_event */
     ) override {
       if (event.jets().size()<2)
         return false;
       return true;
     }
 
-    void scale(const double xsscale) override {
+    void set_xs_scale(double xsscale) override {
+      xs_scale_ = xsscale;
+    }
+
+    void finalise() override {
+      scale();
+      hfile_.Write();
+      hfile_.Close();
+      std::cout << "AnalysisROOT bids you farewell\n"; // be polite
+    }
+
+  private:
+    void scale() {
       std::cout << "AnalysisROOT scales the histograms\n";
 
       TList* list = hfile_.GetListOfKeys() ;
       if(!list) throw std::logic_error{"<E> No keys defined\n"};
       TIter next(list) ;
       TKey* key ;
       TObject* obj ;
 
       while ( (key = (TKey*)next()) ) {
         obj = key->ReadObj() ;
         auto hist = dynamic_cast<TH1*>(obj);
         if(hist) {
           std::cerr << "Histo name: " << hist->GetName()
                     << " title: " << hist->GetTitle()
                     << "will be scaled\n";
-          hist->Scale(xsscale);
+          hist->Scale(xs_scale_);
         } else {
-           std::cerr << "<W> Object " << obj->GetName()
-                     << " is not 1D or 2D histogram : will not be scaled\n";
+          std::cerr << "<W> Object " << obj->GetName()
+                    << " is not 1D or 2D histogram : will not be scaled\n";
         }
       }
     }
-
-    void finalise() override {
-      hfile_.Write();
-      hfile_.Close();
-      std::cout << "AnalysisROOT bids you farewell\n"; // be polite
-    }
-
   };
 }
 
 extern "C"
 __attribute__((visibility("default")))
 std::unique_ptr<HEJ::Analysis> make_analysis(
     YAML::Node const & config, LHEF::HEPRUP const & heprup
 ){
   return std::make_unique<AnalysisROOT>(config, heprup);
 }
diff --git a/examples/AnalysisTemplate.cc b/examples/AnalysisTemplate.cc
index 0269c93..96bb5c6 100644
--- a/examples/AnalysisTemplate.cc
+++ b/examples/AnalysisTemplate.cc
@@ -1,52 +1,52 @@
 //! simple HEJ analysis that doesn't do anything
 
 #include <memory>     // for std::unique_ptr
 #include <iostream>
 
 #include "HEJ/Analysis.hh"
 
 namespace YAML {
   class Node;
 }
 
 namespace LHEF {
   class HEPRUP;
 }
 
 namespace {
   class AnalysisTemplate: public HEJ::Analysis {
     public:
       AnalysisTemplate(
         YAML::Node const & /* config */, LHEF::HEPRUP const & /* heprup */
       ) {}
 
       void fill(
         HEJ::Event const & /* event */,
         HEJ::Event const & /* FO_event */
       ) override {
       }
 
       bool pass_cuts(
         HEJ::Event const & /* event */,
         HEJ::Event const & /* FO_event */
       ) override {
         return true;
       }
 
-      void scale(double /* xsscale */) override {
-        std::cerr <<"AnalysisTemplate Warning:: Function scale does nothing\n";
+      void set_xs_scale(double /* xsscale */) override {
+        // typically store this information for later
       }
 
       void finalise() override {
         std::cout << "bye" << std::endl; // be polite
       }
     };
 }
 
 extern "C"
 __attribute__((visibility("default")))
 std::unique_ptr<HEJ::Analysis> make_analysis(
     YAML::Node const & config, LHEF::HEPRUP const & heprup
 ){
   return std::make_unique<AnalysisTemplate>(config, heprup);
 }
diff --git a/examples/WRapidityAnalysisROOT.cc b/examples/WRapidityAnalysisROOT.cc
index 43ee6e2..1ebdcd0 100644
--- a/examples/WRapidityAnalysisROOT.cc
+++ b/examples/WRapidityAnalysisROOT.cc
@@ -1,157 +1,164 @@
 //! Simple HEJ analysis to plot y_W with ROOT
 //!
 //! See README.md for installation and running instructions.
 //!
 #include <cassert>
 #include <cmath>
 #include <iterator>
 #include <memory>     // for std::unique_ptr
 #include <iostream>   // for std::cout
 
 #include "HEJ/PDG_codes.hh"
 
 #include "TFile.h"
 #include "TKey.h"
 #include "TH1.h"
 
 #include "HEJ/Analysis.hh"
 #include "HEJ/Event.hh"
 #include "HEJ/Version.hh"
 
 namespace YAML {
   class Node;
 }
 
 namespace LHEF {
   class HEPRUP;
 }
 namespace {
   void Fill2DHistogram(TH1D & h,double x,double wt)
   {
     const int binnumber=h.GetXaxis()->FindBin(x);
     const double binwidth=h.GetXaxis()->GetBinWidth(binnumber);
     h.Fill(x,wt/binwidth);
   }
 
   class WRapidityAnalysisROOT: public HEJ::Analysis {
   private:
+    double xs_scale_{1.0};
     TFile hfile_;
     TH1D htotal_,hydif_,hyf_,hyb_,hyW_,hyWsampling_;
   public:
     WRapidityAnalysisROOT(
       YAML::Node const & /* config */, LHEF::HEPRUP const & /* heprup */
     ):hfile_{"HEJ.root","RECREATE",
              ("Output from "+HEJ::Version::package_name_full()).c_str()},
       htotal_{"total","total",1,-0.5,.5},
       hydif_{"ydif","ydif",40,0.,10.},
       hyf_{"yf","yf",40,-5.,5.},
       hyb_{"yb","yb",40,-5.,5.},
       hyW_{"yW","yW",40,-5.,5.},
       hyWsampling_{"yWsampling","yWsampling",40,-5.,5.}
     {
       // Setup histograms
       htotal_.Sumw2();
       hydif_.Sumw2();
       hyf_.Sumw2();
       hyb_.Sumw2();
       hyW_.Sumw2();
     }
 
     void fill(
       HEJ::Event const & event,
       HEJ::Event const & /* FO_event */
     ) override {
       const double weight = event.central().weight;
       htotal_.Fill(0.0,weight);
       auto const & jets = event.jets(); // jets are sorted in rapidity
       // find the W momentum
       auto const & outgoing = event.outgoing();
       auto the_boson = std::find_if(
         outgoing.begin(), outgoing.end(),
         [](auto const & particle) { return std::abs(particle.type) == HEJ::pid::Wp; }
       );
       // we already checked in `pass_cuts` that there is a boson,
       // so the following assert cannot fail
       assert(the_boson != outgoing.end());
       Fill2DHistogram(hydif_,
          std::abs(jets.back().rapidity()-jets.front().rapidity()), weight);
       Fill2DHistogram(hyf_, jets.front().rapidity(), weight);
       Fill2DHistogram(hyb_, jets.back().rapidity(), weight);
       Fill2DHistogram(hyW_, the_boson->rapidity(), weight);
       Fill2DHistogram(hyWsampling_, the_boson->rapidity(),1);
     }
 
     bool pass_cuts(
       HEJ::Event const & event,
       HEJ::Event const & /* FO_event */
     ) override {
       if(event.jets().size() < 2)
         return false;
       auto const & outgoing = event.outgoing();
 
       // Find the W boson
       auto the_W_boson = std::find_if(
         outgoing.begin(), outgoing.end(),
         [](auto const & particle) { return std::abs(particle.type) == HEJ::pid::Wp; }
       );
       if(the_W_boson == outgoing.end())
         return false;
       const auto W_pos = std::distance(outgoing.begin(), the_W_boson);
 
       // Find the decay products
       const auto the_decay = event.decays().find(W_pos);
       if(the_decay == event.decays().end())
         return false;
 
       // Find the charged lepton
       const auto the_charged_lepton = std::find_if(
         the_decay->second.begin(), the_decay->second.end(),
         [](HEJ::Particle const & p) {
           return HEJ::is_anylepton(p) && !HEJ::is_anyneutrino(p);
         }
       );
       if(the_charged_lepton == the_decay->second.end())
         return false;
       return std::abs(the_charged_lepton->rapidity()) > 2.5;
     }
 
-    void scale(const double xsscale) override {
+    void set_xs_scale(double xsscale) override {
+      xs_scale_ = xsscale;
+    }
+
+    void finalise() override {
+      scale();
+      hfile_.Write();
+      hfile_.Close();
+      std::cout << "WRapidityAnalysisROOT bids you farewell\n"; // be polite
+    }
+
+  private:
+    void scale() {
       std::cout << "AnalysisROOT scales the histograms\n";
 
       TList* list = hfile_.GetListOfKeys() ;
       if(!list) throw std::logic_error{"<E> No keys defined\n"};
       TIter next(list) ;
       TKey* key ;
       TObject* obj ;
 
       while ( (key = (TKey*)next()) ) {
         obj = key->ReadObj() ;
         auto hist = dynamic_cast<TH1*>(obj);
         if(hist) {
           std::cerr << "Histo name: " << hist->GetName()
                     << " title: " << hist->GetTitle()
                     << "will be scaled\n";
-          hist->Scale(xsscale);
+          hist->Scale(xs_scale_);
         } else {
            std::cerr << "<W> Object " << obj->GetName()
                      << " is not 1D or 2D histogram : will not be scaled\n";
         }
       }
     }
 
-    void finalise() override {
-      hfile_.Write();
-      hfile_.Close();
-      std::cout << "WRapidityAnalysisROOT bids you farewell\n"; // be polite
-    }
-
   };
 }
 
 extern "C"
 __attribute__((visibility("default")))
 std::unique_ptr<HEJ::Analysis> make_analysis(
     YAML::Node const & config, LHEF::HEPRUP const & heprup
 ){
   return std::make_unique<WRapidityAnalysisROOT>(config, heprup);
 }
diff --git a/include/HEJ/Analysis.hh b/include/HEJ/Analysis.hh
index 3391569..3f6b45c 100644
--- a/include/HEJ/Analysis.hh
+++ b/include/HEJ/Analysis.hh
@@ -1,57 +1,53 @@
 /** \file
  *  \brief Header file for the Analysis interface
  *
  * This header contains declarations that facilitate creating custom analyses
  * to be used with HEJ 2.
  * \todo link to user documentation
  *
  *  \authors   The HEJ collaboration (see AUTHORS for details)
- *  \date      2022
+ *  \date      2022-2023
  *  \copyright GPLv2 or later
  */
 #pragma once
 
 namespace LHEF {
   class HEPRUP;
 }
 
 //! Main HEJ 2 Namespace
 namespace HEJ {
   class Event;
 
   //! Analysis base class
   /**
    *  This is the interface that all analyses should implement,
    *  i.e. all custom analyses have to be derived from this struct.
    */
   struct Analysis {
     /**
      * @param res_event     The event in resummation phase space
      * @param FO_event      The original fixed-order event
      */
     virtual void fill(Event const & res_event, Event const & FO_event) = 0;
     //! Decide whether an event passes the cuts
     /**
      * @param res_event     The event in resummation phase space
      * @param FO_event      The original fixed-order event
      * @returns             Whether the event passes all cuts
      */
     virtual bool pass_cuts(Event const & res_event, Event const & FO_event) = 0;
     //! Finalise analysis
     /**
      * This function is called after all events have been processed and
      * can be used for example to print out or save the results.
      */
     virtual void finalise() = 0;
 
-    /**
-     * This function is called after all events have been processed and
-     * can be used to impose the correct normalisation of weights, in case
-     this is unknown until the very end (e.g. division by number of events).
-     */
-    virtual void scale(double xsscale) = 0;
+    //! Set the ratio (cross section) / (sum of event weights)
+    virtual void set_xs_scale(double scale) = 0;
 
     virtual ~Analysis() = default;
   };
 
 } // namespace HEJ
diff --git a/include/HEJ/RivetAnalysis.hh b/include/HEJ/RivetAnalysis.hh
index 674e39e..64c3955 100644
--- a/include/HEJ/RivetAnalysis.hh
+++ b/include/HEJ/RivetAnalysis.hh
@@ -1,71 +1,71 @@
 /** \file
  *  \brief HEJ 2 interface to rivet analyses
  *
  *  \authors   The HEJ collaboration (see AUTHORS for details)
- *  \date      2019-2022
+ *  \date      2019-2023
  *  \copyright GPLv2 or later
  */
 #pragma once
 
 #include <memory>
 #include <string>
 #include <vector>
 
 #include "LHEF/LHEF.h"
 
 #include "HEJ/Analysis.hh"
 
 namespace YAML {
   class Node;
 }
 
 namespace HEJ {
   /**
     * @brief Class representing a Rivet analysis
     *
     * This class inherits from Analysis and can therefore be used
     * like any other HEJ 2 analysis.
     */
   class RivetAnalysis: public Analysis {
   public:
     //! Create RivetAnalysis
     static std::unique_ptr<Analysis> create(
         YAML::Node const & config, LHEF::HEPRUP const & heprup);
 
     //! Constructor
     /**
      *  @param config    Configuration parameters
      *  @param heprup    General run informations
      *
      * config["rivet"] should be the name of a single Rivet analysis or
      * a list of Rivet analyses. config["output"] is the prefix for
      * the .yoda output files.
      */
     RivetAnalysis(YAML::Node const & config, LHEF::HEPRUP heprup);
     ~RivetAnalysis() override;
     //! Pass an event to the underlying Rivet analysis
     void fill(Event const & event, Event const & /*unused*/) override;
     //! no additional cuts are applied
     bool pass_cuts(Event const & /*unused*/, Event const & /*unused*/
     ) override { return true; }
-    void scale(double scale) override;
+    void set_xs_scale(double scale) override;
     void finalise() override;
 
   private:
     std::vector<std::string> analyses_names_;
     std::string output_name_;
     LHEF::HEPRUP heprup_;
 
     /// struct to organise the infos per rivet run/scale setting
     struct rivet_info;
     std::vector<rivet_info> rivet_runs_;
 
     /**
      *  \internal
      * @brief Calculates the scale variation from the first event for the output
      *        file
      */
     void init(Event const & event);
     bool first_event_;
   };
 } // namespace HEJ
diff --git a/include/HEJ/detail/EmptyAnalysis.hh b/include/HEJ/detail/EmptyAnalysis.hh
index ddc83b6..08563cc 100644
--- a/include/HEJ/detail/EmptyAnalysis.hh
+++ b/include/HEJ/detail/EmptyAnalysis.hh
@@ -1,52 +1,52 @@
 /** \file
  *  \brief Internal header for empty analysis, to be removed in HEJ 2.3
  *
  *  \authors   The HEJ collaboration (see AUTHORS for details)
- *  \date      2022
+ *  \date      2022-2023
  *  \copyright GPLv2 or later
  */
 
 #pragma once
 
 #include <memory>
 
 #include "HEJ/Analysis.hh"
 
 namespace YAML {
   class Node;
 }
 namespace LHEF {
   class HEPRUP;
 }
 
 namespace HEJ {
   class Event;
 
   namespace detail {
     struct EmptyAnalysis: Analysis{
       //! Create EmptyAnalysis
       static std::unique_ptr<Analysis> create(
         YAML::Node const & parameters, LHEF::HEPRUP const & /*unused*/);
 
       //! Fill event into analysis (e.g. to histograms)
       /**
        *  This function does nothing
        */
       void fill(Event const & /*res_event*/, Event const & /*FO_event*/) override;
       //! Whether a resummation event passes all cuts
       /**
        *  There are no cuts, so all events pass
        */
       bool pass_cuts(Event const & /*res_event*/, Event const & /*FO_event*/) override;
+      //! Set the ratio (cross section) / (sum of event weights)
+      void set_xs_scale(double scale) override;
       //! Finalise analysis
       /**
        * This function does nothing
        */
       void finalise() override;
 
-      void scale(double /*xsscale*/) override;
-
       ~EmptyAnalysis() override = default;
     };
   }
 } // namespace HEJ
diff --git a/src/EmptyAnalysis.cc b/src/EmptyAnalysis.cc
index 1d5806a..a95c26c 100644
--- a/src/EmptyAnalysis.cc
+++ b/src/EmptyAnalysis.cc
@@ -1,86 +1,86 @@
 /**
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019-2020
  *  \copyright GPLv2 or later
  */
 #include "HEJ/EmptyAnalysis.hh"
 
 #include <string>
 #include <vector>
 #include <iostream>
 
 #include "yaml-cpp/yaml.h"
 
 #include "HEJ/exceptions.hh"
 
 namespace HEJ::detail {
   namespace {
     std::vector<std::string> param_as_strings(YAML::Node const & parameters){
       using YAML::NodeType;
       switch(parameters.Type()){
       case NodeType::Null:
       case NodeType::Undefined:
         return {};
       case NodeType::Scalar:
         return {parameters.as<std::string>()};
       case NodeType::Sequence: {
         std::vector<std::string> param_strings;
         for(auto && param: parameters){
           param_strings.emplace_back(param.as<std::string>());
         }
         return param_strings;
       }
       case NodeType::Map: {
         std::vector<std::string> param_strings;
         for(auto && param: parameters){
           param_strings.emplace_back(param.first.as<std::string>());
         }
         return param_strings;
       }
       default:;
       }
       throw std::logic_error{"unreachable"};
     }
   } // namespace
 
   std::unique_ptr<Analysis> EmptyAnalysis::create(
       YAML::Node const & parameters, LHEF::HEPRUP const & /*unused*/
   ){
     const auto param_strings = param_as_strings(parameters);
     if(! param_strings.empty()){
       std::string error{"Unknown analysis parameter(s):"};
       for(auto && p: param_strings) error += " " + p;
       throw unknown_option{error};
     }
     return std::make_unique<EmptyAnalysis>();
   }
 
   void EmptyAnalysis::fill(
       Event const & /*res_event*/, Event const & /*FO_event*/
   ){
     // do nothing
   }
 
   bool EmptyAnalysis::pass_cuts(
       Event const & /*res_event*/, Event const & /*FO_event*/
   ){
     return true;
   }
 
-  void EmptyAnalysis::scale(double /*xsscale*/){
+  void EmptyAnalysis::set_xs_scale(double /*xsscale*/){
     // do nothing
   }
 
   void EmptyAnalysis::finalise(){
     // do nothing
   }
 
 } // namespace HEJ::detail
 
 namespace HEJ {
   std::unique_ptr<Analysis> EmptyAnalysis::create(
     YAML::Node const & parameters, LHEF::HEPRUP const & heprup
   ){
     return detail::EmptyAnalysis::create(parameters, heprup);
   }
 }
diff --git a/src/RivetAnalysis.cc b/src/RivetAnalysis.cc
index 13e7ad3..6e148ec 100644
--- a/src/RivetAnalysis.cc
+++ b/src/RivetAnalysis.cc
@@ -1,208 +1,199 @@
 /**
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019-2020
  *  \copyright GPLv2 or later
  */
 #include "HEJ/RivetAnalysis.hh"
 
 #include "HEJ/ConfigFlags.hh"
 #include "HEJ/utility.hh"
 
 #ifdef HEJ_BUILD_WITH_RIVET
 
 #include <cstddef>
 #include <ostream>
 #include <utility>
 
 #include "yaml-cpp/yaml.h"
 
 #include "Rivet/AnalysisHandler.hh"
 #include "Rivet/Config/RivetConfig.hh"
 
 #ifdef RIVET_ENABLE_HEPMC_3
 
 #include "HepMC3/GenEvent.h"
 #include "HEJ/HepMC3Interface.hh"
 
 #else // rivet with HepMC 2
 
 #include "HEJ/HepMC2Interface.hh"
 #include "HepMC/GenEvent.h"
 
 #endif // RIVET_ENABLE_HEPMC_3
 
 #include "HEJ/Event.hh"
 #include "HEJ/exceptions.hh"
 
 #endif
 
 namespace HEJ {
   std::unique_ptr<Analysis> RivetAnalysis::create(
       YAML::Node const & config, LHEF::HEPRUP const & heprup
   ){
     return std::make_unique<RivetAnalysis>(config, heprup);
   }
 }
 
 #ifdef HEJ_BUILD_WITH_RIVET
 
 namespace HEJ {
 #ifdef RIVET_ENABLE_HEPMC_3
   using HepMCInterface=HepMC3Interface;
 #else
   using HepMCInterface=HepMC2Interface;
 #endif
 
   struct RivetAnalysis::rivet_info {
     rivet_info(std::unique_ptr<Rivet::AnalysisHandler> && h,
         std::string && s,  HepMCInterface && m):
       handler{std::move(h)}, name{std::move(s)}, hepmc{std::move(m)}
     {}
     // AnalysisHandler doesn't allow a copy constructor -> use ptr
     std::unique_ptr<Rivet::AnalysisHandler> handler;
     std::string name;
     HepMCInterface hepmc;
   };
 
   RivetAnalysis::RivetAnalysis(YAML::Node const & config,
                                LHEF::HEPRUP heprup
   ):
     heprup_{std::move(heprup)},
     first_event_(true)
   {
 
     // assert that only supported options are provided
     if(!config.IsMap()) throw invalid_type{"Rivet config has to be a map!"};
     for(auto const & entry: config){
       const auto name = entry.first.as<std::string>();
       if(name != "output" && name != "rivet")
         throw unknown_option{"Unknown option \'"+name+"\' provided to rivet."};
     }
 
     // get output name
     if(!config["output"])
       throw std::invalid_argument{
         "No output was provided to rivet. "
         "Either give an output or deactivate rivet."
       };
     if(!config["output"].IsScalar())
       throw invalid_type{
         "Output for rivet must be a scalar."
       };
     output_name_ = config["output"].as<std::string>();
 
     // read in analyses
     auto const & name_node = config["rivet"];
     switch(name_node.Type()){
     case YAML::NodeType::Scalar:
       analyses_names_.push_back(name_node.as<std::string>());
       break;
     case YAML::NodeType::Sequence:
       for(YAML::const_iterator it = name_node.begin(); it != name_node.end(); ++it){
         analyses_names_.push_back(it->as<std::string>());
       }
       break;
     default:
       throw std::invalid_argument{
         "No analysis was provided to rivet. "
         "Either give an analysis or deactivate rivet."
       };
     }
   }
 
   // it is a bit ugly that we can't do this directly in `initialise`
   void RivetAnalysis::init(Event const & event){
 #ifdef HEJ_USE_RIVET2
     rivet_runs_.reserve(event.variations().size()+1);
 #else
     ignore(event); // shut up compiler
 #endif
     rivet_runs_.emplace_back( std::make_unique<Rivet::AnalysisHandler>(),
       std::string(""), HepMCInterface(heprup_)
     );
     rivet_runs_.back().handler->addAnalyses(analyses_names_);
 
 #ifdef HEJ_USE_RIVET2
     //! scale variation in rivet 2 through multiple analyses
     if( !event.variations().empty() ){
       for(auto const & vari : event.variations()){
         rivet_runs_.emplace_back( std::make_unique<Rivet::AnalysisHandler>(),
           "."+to_simple_string(*vari.description), HepMCInterface(heprup_)
         );
         rivet_runs_.back().handler->addAnalyses(analyses_names_);
       }
     }
 #endif
   }
 
   void RivetAnalysis::fill(Event const & event, Event const & /*unused*/){
     if(first_event_){
       first_event_=false;
       init(event);
     }
     auto hepmc_event(rivet_runs_.front().hepmc(event));
     rivet_runs_.front().handler->analyze(hepmc_event);
 
 #ifdef HEJ_USE_RIVET2
     for(std::size_t i = 1; i < rivet_runs_.size(); ++i){
       auto & run = rivet_runs_[i];
       run.hepmc.set_central(hepmc_event, event, i-1);
       run.handler->analyze(hepmc_event);
     }
 #endif
   }
 
-  void RivetAnalysis::scale(double xsscale){
-    // used to scale the cross sections and histograms before .finalise in case the
-    // normalisation is unknown until the end of the run
-    for(auto & run: rivet_runs_){
-      // get yoda analysis objects:
-      std::vector<YODA::AnalysisObjectPtr> yodaAOS=run.handler->getYodaAOs(true);
-      for (auto& t : yodaAOS) {
-        auto * obj = dynamic_cast<YODA::Fillable*>(t.get());
-        if(obj) {
-          obj->scaleW(xsscale);
-        }
-      }
+  void RivetAnalysis::set_xs_scale(double scale) {
+    for(auto & run: rivet_runs_) {
+      run.hepmc.set_xs_scale(scale);
     }
   }
 
   void RivetAnalysis::finalise(){
     for(auto & run: rivet_runs_){
       run.handler->finalize();
       run.handler->writeData(output_name_+run.name+std::string(".yoda"));
     }
   }
 } // namespace HEJ
 
 #else // no rivet => create empty analysis
 
 namespace Rivet {
   class AnalysisHandler {};
 }
 
 namespace HEJ {
   struct RivetAnalysis::rivet_info{};
 
   RivetAnalysis::RivetAnalysis(
     YAML::Node const & /*config*/, LHEF::HEPRUP /*heprup*/
   ){
     ignore(first_event_);
     throw std::invalid_argument(
         "Failed to create RivetAnalysis: "
         "HEJ 2 was built without rivet support"
     );
   }
 
   void RivetAnalysis::init(Event const & /*event*/){}
   void RivetAnalysis::fill(Event const & /*event*/, Event const & /*unused*/){}
-  void RivetAnalysis::scale(double /*xsscale*/){}
+  void RivetAnalysis::set_xs_scale(double /*xsscale*/){}
 
   void RivetAnalysis::finalise(){}
 } // namespace HEJ
 
 #endif
 
 namespace HEJ {
   RivetAnalysis::~RivetAnalysis() = default;
 }
diff --git a/src/bin/HEJ.cc b/src/bin/HEJ.cc
index 383320f..365816d 100644
--- a/src/bin/HEJ.cc
+++ b/src/bin/HEJ.cc
@@ -1,448 +1,444 @@
 /**
  *  \authors   The HEJ collaboration (see AUTHORS for details)
  *  \date      2019-2022
  *  \copyright GPLv2 or later
  */
 #include <array>
 #include <chrono>
 #include <iostream>
 #include <limits>
 #include <memory>
 #include <numeric>
 #include <optional>
 
 
 #include "yaml-cpp/yaml.h"
 
 #include "fastjet/ClusterSequence.hh"
 
 #include "HEJ/CombinedEventWriter.hh"
 #include "HEJ/Config.hh"
 #include "HEJ/CrossSectionAccumulator.hh"
 #include "HEJ/Event.hh"
 #include "HEJ/EventReader.hh"
 #include "HEJ/BufferedEventReader.hh"
 #include "HEJ/EventReweighter.hh"
 #include "HEJ/filesystem.hh"
 #include "HEJ/get_analysis.hh"
 #include "HEJ/make_RNG.hh"
 #include "HEJ/ProgressBar.hh"
 #include "HEJ/stream.hh"
 #include "HEJ/Unweighter.hh"
 #include "HEJ/Version.hh"
 #include "HEJ/YAMLreader.hh"
 
 HEJ::Config load_config(char const * filename){
   try{
     return HEJ::load_config(filename);
   }
   catch(std::exception const & exc){
     std::cerr << "Error: " << exc.what() << '\n';
     std::exit(EXIT_FAILURE);
   }
 }
 
 std::vector<std::unique_ptr<HEJ::Analysis>> get_analyses(
     std::vector<YAML::Node> const & parameters, LHEF::HEPRUP const & heprup
 ){
   try{
     return HEJ::get_analyses(parameters, heprup);
   }
   catch(std::exception const & exc){
     std::cerr << "Failed to load analysis: " << exc.what() << '\n';
     std::exit(EXIT_FAILURE);
   }
 }
 
 std::string time_to_string(const time_t time){
   char s[30];
   struct tm * p = localtime(&time);
   strftime(s, 30, "%a %b %d %Y %H:%M:%S", p);
   return s;
 }
 
 void check_one_parton_per_jet(HEJ::Event const & ev) {
   const std::size_t npartons = std::count_if(
     ev.outgoing().begin(), ev.outgoing().end(),
     [](HEJ::Particle const & p) { return is_parton(p); }
   );
   if(ev.jets().size() < npartons) {
     throw std::invalid_argument{
       "Number of partons (" + std::to_string(npartons) + ") in input event"
         " does not match number of jets (" + std::to_string(ev.jets().size())+  ")"
       };
   }
 }
 
 HEJ::Event to_event(
     LHEF::HEPEUP const & hepeup,
     HEJ::JetParameters const & fixed_order_jets,
     HEJ::EWConstants const & ew_parameters,
     const double off_shell_tolerance
 ) {
     HEJ::Event::EventData event_data{hepeup};
     event_data.reconstruct_intermediate(ew_parameters);
     event_data.repair_momenta(off_shell_tolerance);
 
     HEJ::Event ev{
       std::move(event_data).cluster(
           fixed_order_jets.def, fixed_order_jets.min_pt
       )
     };
     check_one_parton_per_jet(ev);
     return ev;
 }
 
 void unweight(
   HEJ::Unweighter & unweighter,
   HEJ::WeightType weight_type,
   std::vector<HEJ::Event> & events,
   HEJ::RNG & ran
 ) {
     if(weight_type == HEJ::WeightType::unweighted_resum){
       unweighter.set_cut_to_maxwt(events);
     }
     events.erase(
       unweighter.unweight(begin(events), end(events), ran),
       end(events)
     );
 }
 
 // peek up to nevents events from reader
 std::vector<LHEF::HEPEUP> peek_events(
     HEJ::BufferedEventReader & reader,
     const int nevents
 ) {
   std::vector<LHEF::HEPEUP> events;
   while(
       static_cast<int>(events.size()) < nevents
       && reader.read_event()
   ) {
     events.emplace_back(reader.hepeup());
   }
   // put everything back into the reader
   for(auto it = rbegin(events); it != rend(events); ++it) {
     reader.emplace(*it);
   }
   return events;
 }
 
 void append_resummed_events(
     std::vector<HEJ::Event> & resummation_events,
     HEJ::EventReweighter & reweighter,
     LHEF::HEPEUP const & hepeup,
     const size_t trials,
     HEJ::JetParameters const & fixed_order_jets,
     HEJ::EWConstants const & ew_parameters,
     const double off_shell_tolerance
 ) {
   const HEJ::Event FO_event = to_event(
     hepeup,
     fixed_order_jets,
     ew_parameters,
     off_shell_tolerance
   );
   if(reweighter.treatment(FO_event.type()) != HEJ::EventTreatment::reweight) {
     return;
   }
   const auto resummed = reweighter.reweight(FO_event, trials);
   resummation_events.insert(
       end(resummation_events),
       begin(resummed), end(resummed)
   );
 }
 
 std::optional<HEJ::ProgressBar<std::size_t>> make_progress_bar(
   std::optional<std::size_t> n_input_events,
   char const * event_file
 ) {
   if(n_input_events) {
     return HEJ::ProgressBar{std::cout, *n_input_events};
   }
   if(HEJ::is_regular(event_file)) {
       auto t_reader = HEJ::make_reader(event_file);
       std::size_t n_input_events = 0;
       while(t_reader->read_event()) ++n_input_events;
       return HEJ::ProgressBar{std::cout, n_input_events};
   }
   return {};
 }
 
 void train(
     HEJ::Unweighter & unweighter,
     HEJ::BufferedEventReader & reader,
     HEJ::EventReweighter & reweighter,
     const size_t total_trials,
     const double max_dev,
     double reweight_factor,
     HEJ::JetParameters const & fixed_order_jets,
     HEJ::EWConstants const & ew_parameters,
     const double off_shell_tolerance
 ) {
   std::cout << "Reading up to " << total_trials << " training events...\n";
   auto FO_events = peek_events(reader, total_trials);
   if(FO_events.empty()) {
     throw std::runtime_error{
       "No events generated to calibrate the unweighting weight!"
       "Please increase the number \"trials\" or deactivate the unweighting."
     };
   }
   const size_t trials = total_trials/FO_events.size();
   // adjust reweight factor so that the overall normalisation
   // is the same as in the full run
   reweight_factor *= trials;
   for(auto & hepeup: FO_events) {
     hepeup.XWGTUP *= reweight_factor;
   }
   std::cout << "Training unweighter with "
             << trials << '*' << FO_events.size() << " events\n";
 
   auto progress = HEJ::ProgressBar<int>{
      std::cout, static_cast<int>(FO_events.size())
   };
   std::vector<HEJ::Event> resummation_events;
   for(auto const & hepeup: FO_events) {
     append_resummed_events(
         resummation_events,
         reweighter, hepeup, trials, fixed_order_jets,
         ew_parameters,
         off_shell_tolerance
     );
     ++progress;
   }
 
   unweighter.set_cut_to_peakwt(resummation_events, max_dev);
   std::cout << "\nUnweighting events with weight up to "
             << unweighter.get_cut() << '\n';
 }
 
 int main(int argn, char** argv) {
   using clock = std::chrono::system_clock;
 
   if (argn != 3) {
     std::cerr << "\n# Usage:\n."<< argv[0] <<" config_file input_file\n\n";
     return EXIT_FAILURE;
   }
 
   const auto start_time = clock::now();
   {
     std::cout << "Starting " << HEJ::Version::package_name_full()
       << ", revision " << HEJ::Version::revision() << " ("
       << time_to_string(clock::to_time_t(start_time)) << ")" << std::endl;
   }
   fastjet::ClusterSequence::print_banner();
 
   // read configuration
   const HEJ::Config config = load_config(argv[1]);
   auto reader = HEJ::make_reader(argv[2]);
   assert(reader);
 
   auto heprup{ reader->heprup() };
   heprup.generators.emplace_back(LHEF::XMLTag{});
   heprup.generators.back().name = HEJ::Version::package_name();
   heprup.generators.back().version = HEJ::Version::String();
 
   auto analyses = get_analyses( config.analyses_parameters, heprup );
   assert(analyses.empty() || analyses.front() != nullptr);
 
   HEJ::CombinedEventWriter writer{config.output, std::move(heprup)};
 
   double global_reweight = 1.;
   auto const & max_events = config.max_events;
 
   // if we need the event number:
   if(std::abs(heprup.IDWTUP) == 4 || std::abs(heprup.IDWTUP) == 1 || max_events){
     // try to read from LHE head
     auto input_events{reader->number_events()};
     if(!input_events) {
       // else count manually
       auto t_reader = HEJ::make_reader(argv[2]);
       input_events = 0;
       while(t_reader->read_event()) ++(*input_events);
     }
     if(std::abs(heprup.IDWTUP) == 4 || std::abs(heprup.IDWTUP) == 1){
       // IDWTUP 4 or 1 assume average(weight)=xs, but we need sum(weights)=xs
       std::cout << "Found IDWTUP " << heprup.IDWTUP << ": "
         << "assuming \"cross section = average weight\".\n"
         << "converting to \"cross section = sum of weights\" ";
       global_reweight /= *input_events;
     }
     if(max_events && (*input_events > *max_events)){
       // maximal number of events given
       global_reweight *= *input_events/static_cast<double>(*max_events);
       std::cout << "Processing " << *max_events
                 << " out of " << *input_events << " events\n";
     }
   }
 
   HEJ::ScaleGenerator scale_gen{
     config.scales.base,
     config.scales.factors,
     config.scales.max_ratio
   };
   std::shared_ptr<HEJ::RNG> ran{
     HEJ::make_RNG(config.rng.name, config.rng.seed)};
   assert(ran != nullptr);
   HEJ::EventReweighter hej{
     reader->heprup(),
     std::move(scale_gen),
     to_EventReweighterConfig(config),
     ran
   };
 
   std::optional<HEJ::Unweighter> unweighter{};
   if(config.weight_type != HEJ::WeightType::weighted) {
     unweighter = HEJ::Unweighter{};
   }
   if(config.weight_type == HEJ::WeightType::partially_unweighted) {
     HEJ::BufferedEventReader buffered_reader{std::move(reader)};
     assert(config.unweight_config);
     train(
         *unweighter,
         buffered_reader,
         hej,
         config.unweight_config->trials,
         config.unweight_config->max_dev,
         global_reweight/config.trials,
         config.fixed_order_jets,
         config.ew_parameters,
         config.off_shell_tolerance
     );
     reader = std::make_unique<HEJ::BufferedEventReader>(
         std::move(buffered_reader)
     );
   }
 
   // status infos & eye candy
   if (config.nlo.enabled){
     std::cout << "HEJ@NLO Truncation Enabled for NLO order: "
     << config.nlo.nj << std::endl;
   }
   size_t nevent = 0;
   std::array<int, HEJ::event_type::last_type + 1>
     nevent_type{0}, nfailed_type{0};
 
   auto progress = make_progress_bar(reader->number_events(), argv[2]);
   if(!progress) {
     std::cout
       << "Cannot determine number of input events. "
       "Disabling progress bar"
       << std::endl;
   }
 
   HEJ::CrossSectionAccumulator xs;
   std::map<HEJ::StatusCode, int> status_counter;
   size_t total_trials = 0;
   size_t total_resum = 0;
 
   // Loop over the events in the input file
   while(reader->read_event() && (!max_events || nevent < *max_events) ){
     ++nevent;
 
     // reweight events so that the total cross section is conserved
     auto hepeup = reader->hepeup();
     hepeup.XWGTUP *= global_reweight;
 
     const auto FO_event = to_event(
       hepeup,
       config.fixed_order_jets,
       config.ew_parameters,
       config.off_shell_tolerance
     );
     if(FO_event.central().weight == 0) {
       static const bool warned_once = [argv,nevent](){
         std::cerr
           << "WARNING: event number " << nevent
           << " in " << argv[2] << " has zero weight. "
           "Ignoring this and all further events with vanishing weight.\n";
         return true;
       }();
       (void) warned_once; // shut up compiler warnings
       continue;
     }
 
     auto resummed_events{ hej.reweight(FO_event, config.trials) };
 
     // some bookkeeping
     for(auto const & s: hej.status())
       ++status_counter[s];
     total_trials+=hej.status().size();
     ++nevent_type[FO_event.type()];
 
     if(resummed_events.empty()) ++nfailed_type[FO_event.type()];
 
     if(unweighter) {
       unweight(*unweighter, config.weight_type, resummed_events, *ran);
     }
 
     // analysis
     for(auto & ev: resummed_events){
       //TODO: move pass_cuts to after phase space point generation
       bool passed = analyses.empty();
       for(auto const & analysis: analyses){
+        analysis->set_xs_scale(reader->scalefactor());
         if(analysis->pass_cuts(ev, FO_event)){
           passed = true;
           analysis->fill(ev, FO_event);
         };
       }
       if(passed){
         writer.set_xs_scale(reader->scalefactor());
         writer.write(ev);
       } else {
         ev.parameters()*=0; // do not use discarded events afterwards
       }
     }
     xs.fill_correlated(resummed_events);
     total_resum += resummed_events.size();
     if(progress) ++*progress;
   } // main event loop
   std::cout << '\n';
   // scale to account for num_trials of sherpa
-  if(const double scalefactor = reader->scalefactor(); scalefactor != 1.) {
-    for(auto const & analysis: analyses){
-      analysis->scale(scalefactor);
-    }
-    xs.scale(scalefactor);
-  }
+  xs.scale(reader->scalefactor());
 
   for(auto const & analysis: analyses){
     analysis->finalise();
   }
   writer.finish();
 
   using namespace HEJ::event_type;
   std::cout<< "Events processed: " << nevent << " (" << total_resum << " resummed)"<< '\n';
   std::cout << '\t' << name(EventType::first_type) << ": "
             << nevent_type[EventType::first_type]
             << ", failed to reconstruct " << nfailed_type[EventType::first_type]
             << '\n';
   for(auto i=EventType::first_type+1; i<=EventType::last_type; i*=2){
     std::cout << '\t' << name(static_cast<EventType>(i)) << ": "
               << nevent_type[i]
               << ", failed to reconstruct " << nfailed_type[i]
               << '\n';
   }
 
   std::cout << '\n' << xs << '\n';
 
   std::cout << "Generation statistic: "
     << status_counter[HEJ::StatusCode::good] << "/" << total_trials
     << " trials successful.\n";
   for(auto && entry: status_counter){
     const double fraction = static_cast<double>(entry.second)/total_trials;
     const int percent = std::round(100*fraction);
     std::cout << std::left << std::setw(17) << (to_string(entry.first) + ":")
               << " [";
     for(int i = 0; i < percent/2; ++i) std::cout << '#';
     for(int i = percent/2; i < 50; ++i) std::cout << ' ';
     std::cout << "] " <<std::setw(2)<<std::right<< percent << "%\n";
   }
 
   std::chrono::duration<double> run_time = (clock::now() - start_time);
   std::cout << "\nFinished " << HEJ::Version::package_name() << " at "
     << time_to_string(clock::to_time_t(clock::now()))
     << "\n=> Runtime: " << run_time.count() << " sec ("
     << nevent/run_time.count() << " Events/sec).\n";
 
   return EXIT_SUCCESS;
 }