diff --git a/examples/WRapidityAnalysisROOT.cc b/examples/WRapidityAnalysisROOT.cc
index 12e88d3..b57f245 100644
--- a/examples/WRapidityAnalysisROOT.cc
+++ b/examples/WRapidityAnalysisROOT.cc
@@ -1,109 +1,131 @@
 //! Simple HEJ analysis to plot y_W with ROOT
-// To compile analysis:
-// g++ `HEJ-config --cxxflags` `root-config --cflags --libs` `fastjet-config --cxxflags` -fPIC -shared -O2 -fvisibility=hidden -Wl,-soname,libWRapidityAnalysisROOT.so -o libWRapidityAnalysisROOT.so WRapidityAnalysisROOT.cc
-// To run HEJFOG:
-// LD_PRELOAD=/usr/lib64/root/libCore.so.6.26 bin/HEJFOG ../hej/FixedOrderGen/configFO.yml
+//!
+//! 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 "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;
+      const double weight = event.central().weight;
       htotal_.Fill(0.0,weight);
-      auto const & jets=event.jets(); // jets are sorted in rapidity
+      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 HEJ::is_AWZ_boson(particle); });
+      auto const & outgoing = event.outgoing();
+      auto the_boson = std::find_if(
+        outgoing.begin(), outgoing.end(),
+        [](auto const & particle) { return std::abs(particle.pid) == 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_,
-          fabs(jets.back().rapidity()-jets.front().rapidity()), weight);
+         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)
+      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.pid) == HEJ::pid::Wp; }
+      );
+      if(the_W_boson == outgoing.end())
         return false;
-      // Find the charged lepton momentum
-      const auto the_decay = event.decays().begin();
-      const auto the_charged_lepton = *std::find_if(
+      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);
         }
       );
-      // momentum: the_charged_lepton.p
-      if (fabs(the_charged_lepton.p.rapidity())>2.5)
+      if(the_charged_lepton == the_decay->second.end())
         return false;
-      return true;
+      return std::abs(the_charged_lepton->rapidity()) > 2.5;
     }
 
     void finalise() override {
       hfile_.Write();
       hfile_.Close();
       std::cout << "WRapidityAnalysisROOT bids you farewell " << 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<WRapidityAnalysisROOT>(config, heprup);
 }