diff --git a/examples/WRapidityAnalysisROOT.cc b/examples/WRapidityAnalysisROOT.cc
index 5668bbd..43ee6e2 100644
--- a/examples/WRapidityAnalysisROOT.cc
+++ b/examples/WRapidityAnalysisROOT.cc
@@ -1,131 +1,157 @@
 //! 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:
     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 {
+      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);
+        } 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 " << std::endl; // be polite
+      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);
 }