diff --git a/include/Rivet/Tools/RivetHepMC.hh b/include/Rivet/Tools/RivetHepMC.hh
--- a/include/Rivet/Tools/RivetHepMC.hh
+++ b/include/Rivet/Tools/RivetHepMC.hh
@@ -1,96 +1,96 @@
 // -*- C++ -*-
 #ifndef RIVET_RivetHepMC_HH
 #define RIVET_RivetHepMC_HH
 
 #ifdef ENABLE_HEPMC_3
 #include "HepMC3/HepMC3.h"
 #include "HepMC3/Relatives.h"
 #include "HepMC3/Reader.h"
 
 namespace Rivet{
   namespace RivetHepMC = HepMC3;
   using RivetHepMC::ConstGenParticlePtr;
   using RivetHepMC::ConstGenVertexPtr;
   using RivetHepMC::Relatives;
   using RivetHepMC::ConstGenHeavyIonPtr;
   
   using HepMC_IO_type = RivetHepMC::Reader;
 
   using PdfInfo = RivetHepMC::GenPdfInfo;
 }
 
 #else
 #include "HepMC/GenEvent.h"
 #include "HepMC/GenParticle.h"
 #include "HepMC/GenVertex.h"
 #include "HepMC/Version.h"
 #include "HepMC/GenRanges.h"
 #include "HepMC/IO_GenEvent.h"
 
 
 namespace Rivet{
   namespace RivetHepMC = HepMC;
   
   // HepMC 2.07 provides its own #defines
   #define ConstGenParticlePtr const HepMC::GenParticle*
   #define ConstGenVertexPtr const HepMC::GenVertex*
   #define ConstGenHeavyIonPtr const HepMC::HeavyIon*
   
   /// @brief Replicated the HepMC3 Relatives syntax using HepMC2 IteratorRanges
   /// This is necessary mainly because of capitalisation differences
   class Relatives{
     
     public:
     
     constexpr Relatives(HepMC::IteratorRange relo): _internal(relo){}
     
     constexpr HepMC::IteratorRange operator()() const {return _internal;}
     operator HepMC::IteratorRange() const {return _internal;}
     
     const static Relatives PARENTS;
     const static Relatives CHILDREN;
     const static Relatives ANCESTORS;
     const static Relatives DESCENDANTS;
     
     private:
     const HepMC::IteratorRange _internal;
     
   };
   
   using HepMC_IO_type = HepMC::IO_GenEvent;
   using PdfInfo = RivetHepMC::PdfInfo;
 
 }
   
 #endif
 
 #include "Rivet/Tools/RivetSTL.hh"
 #include "Rivet/Tools/Exceptions.hh"
 
 
 namespace Rivet {
 
   using RivetHepMC::GenEvent;
   using ConstGenEventPtr = std::shared_ptr<const GenEvent>;
   /// @todo Use mcutils?
 
   namespace HepMCUtils{
   
     std::vector<ConstGenParticlePtr> particles(ConstGenEventPtr ge);
     std::vector<ConstGenParticlePtr> particles(const GenEvent *ge);
     std::vector<ConstGenVertexPtr>   vertices(ConstGenEventPtr ge);
     std::vector<ConstGenVertexPtr>   vertices(const GenEvent *ge);
     std::vector<ConstGenParticlePtr> particles(ConstGenVertexPtr gv, const Relatives &relo);
     std::vector<ConstGenParticlePtr> particles(ConstGenParticlePtr gp, const Relatives &relo);
     int uniqueId(ConstGenParticlePtr gp);
     int particles_size(ConstGenEventPtr ge);
     int particles_size(const GenEvent *ge);
-    std::vector<ConstGenParticlePtr> beams(const GenEvent *ge);
+    std::pair<ConstGenParticlePtr,ConstGenParticlePtr> beams(const GenEvent *ge);
     std::shared_ptr<HepMC_IO_type> makeReader(std::istream &istr,
                                               std::string * errm = 0);
     bool readEvent(std::shared_ptr<HepMC_IO_type> io,
                    std::shared_ptr<GenEvent> evt);
   }
 }
 
 #endif
diff --git a/src/Core/Particle.cc b/src/Core/Particle.cc
--- a/src/Core/Particle.cc
+++ b/src/Core/Particle.cc
@@ -1,299 +1,296 @@
 #include "Rivet/Particle.hh"
 #include "Rivet/Tools/Cuts.hh"
 #include "Rivet/Tools/ParticleIdUtils.hh"
 #include "Rivet/Tools/ParticleUtils.hh"
 
 namespace Rivet {
 
 
   void Particle::setConstituents(const vector<Particle>& cs, bool setmom) {
     _constituents = cs;
     if (setmom) _momentum = sum(cs, p4, FourMomentum());
   }
 
 
   void Particle::addConstituent(const Particle& c, bool addmom) {
     _constituents += c;
     if (addmom) _momentum += c;
   }
 
 
   void Particle::addConstituents(const vector<Particle>& cs, bool addmom) {
     _constituents += cs;
     if (addmom)
       for (const Particle& c : cs)
         _momentum += c;
   }
 
 
   vector<Particle> Particle::rawConstituents() const {
     if (!isComposite()) return Particles{*this};
     vector<Particle> rtn;
     for (const Particle& p : constituents()) rtn += p.rawConstituents();
     return rtn;
   }
 
 
   Particle& Particle::transformBy(const LorentzTransform& lt) {
     _momentum = lt.transform(_momentum);
     return *this;
   }
 
 
   bool Particle::isVisible() const {
     // Charged particles are visible
     if ( PID::threeCharge(pid()) != 0 ) return true;
     // Neutral hadrons are visible
     if ( PID::isHadron(pid()) ) return true;
     // Photons are visible
     if ( pid() == PID::PHOTON ) return true;
     // Gluons are visible (for parton level analyses)
     if ( pid() == PID::GLUON ) return true;
     // Everything else is invisible
     return false;
   }
 
 
   bool Particle::isStable() const {
     return genParticle() != NULL &&
       genParticle()->status() == 1 &&
       genParticle()->end_vertex() == NULL;
   }
 
 
   vector<Particle> Particle::ancestors(const Cut& c, bool physical_only) const {
     
     /// @todo: use HepMC FindParticles when it actually works properly for const objects
     //HepMC::FindParticles searcher(genParticle(), HepMC::ANCESTORS, HepMC::STATUS==1);
     
     vector<Particle> rtn;
 
     // this case needed protecting against (at least for the latest Herwig... not sure why
     // it didn't show up earlier
     if (genParticle() == nullptr) return rtn;
 
     ConstGenVertexPtr gv = genParticle()->production_vertex();
 
     if (gv == nullptr) return rtn;
     
     vector<ConstGenParticlePtr> ancestors = HepMCUtils::particles(genParticle(), Relatives::ANCESTORS);
     
     for(const auto &a: ancestors){
       if(physical_only && a->status() != 1 && a->status() != 2) continue;
       const Particle p(a);
       if(c != Cuts::OPEN && !c->accept(p)) continue;
       rtn += p;
     }
 
     return rtn;
   }
 
 
   vector<Particle> Particle::parents(const Cut& c) const {
     vector<Particle> rtn;
     ConstGenVertexPtr gv = genParticle()->production_vertex();
     if (gv == nullptr) return rtn;
     for(ConstGenParticlePtr it: HepMCUtils::particles(gv, Relatives::PARENTS)){
     const Particle p(it);
       if (c != Cuts::OPEN && !c->accept(p)) continue;
       rtn += p;
     }
     return rtn;
   }
 
 
   vector<Particle> Particle::children(const Cut& c) const {
     vector<Particle> rtn;
     if (isStable()) return rtn;
     ConstGenVertexPtr gv = genParticle()->end_vertex();
     if (gv == nullptr) return rtn;
     /// @todo Would like to do this, but the range objects are broken
     // foreach (const GenParticlePtr gp, gv->particles(HepMC::children))
     //   rtn += Particle(gp);
     for(ConstGenParticlePtr it: HepMCUtils::particles(gv, Relatives::CHILDREN)){
       const Particle p(it);
       if (c != Cuts::OPEN && !c->accept(p)) continue;
       rtn += p;
     }
     return rtn;
   }
 
 
   /// @todo Insist that the current particle is post-hadronization, otherwise throw an exception?
   /// @todo Use recursion through replica-avoiding functions to avoid bookkeeping duplicates
   vector<Particle> Particle::allDescendants(const Cut& c, bool remove_duplicates) const {
     vector<Particle> rtn;
     if (isStable()) return rtn;
 
     ConstGenVertexPtr gv = genParticle()->end_vertex();
     if (gv == nullptr) return rtn;
     /// @todo Would like to do this, but the range objects are broken
     // foreach (const GenParticlePtr gp, gv->particles(HepMC::descendants))
     for(ConstGenParticlePtr it: HepMCUtils::particles(gv, Relatives::DESCENDANTS)){
       const Particle p(it);
       if (c != Cuts::OPEN && !c->accept(p)) continue;
       if (remove_duplicates && it->end_vertex() != NULL) {
         // size_t n = 0; ///< @todo Only remove 1-to-1 duplicates?
         bool dup = false;
         /// @todo Yuck, HepMC
         for(ConstGenParticlePtr it2: HepMCUtils::particles(it->end_vertex(), Relatives::CHILDREN)){
           // n += 1; if (n > 1) break;
           if (it->pdg_id() == it2->pdg_id()) { dup = true; break; }
         }
         if (dup) continue;
       }
       rtn += p;
     }
     return rtn;
   }
 
 
   /// @todo Insist that the current particle is post-hadronization, otherwise throw an exception?
   vector<Particle> Particle::stableDescendants(const Cut& c) const {
     vector<Particle> rtn;
     if (isStable()) return rtn;
     ConstGenVertexPtr gv = genParticle()->end_vertex();
     if (gv == nullptr) return rtn;
     /// @todo Would like to do this, but the range objects are broken
     // foreach (const GenParticlePtr gp, gv->particles(HepMC::descendants))
     for(ConstGenParticlePtr it: HepMCUtils::particles(gv, Relatives::DESCENDANTS)){
     //for (GenVertex::particle_iterator it = gv->particles_begin(HepMC::descendants); it != gv->particles_end(HepMC::descendants); ++it) {
       // if ((*it)->status() != 1 || (*it)->end_vertex() != NULL) continue;
       const Particle p(*it);
       if (!p.isStable()) continue;
       if (c != Cuts::OPEN && !c->accept(p)) continue;
       rtn += p;
     }
     return rtn;
   }
 
 
   double Particle::flightLength() const {
     if (isStable()) return -1;
     if (genParticle() == NULL) return 0;
     if (genParticle()->production_vertex() == NULL) return 0;
     const RivetHepMC::FourVector v1 = genParticle()->production_vertex()->position();
     const RivetHepMC::FourVector v2 = genParticle()->end_vertex()->position();
     return sqrt(sqr(v2.x()-v1.x()) + sqr(v2.y()-v1.y()) + sqr(v2.z()-v1.z()));
   }
 
 
   bool Particle::hasParent(PdgId pid) const {
     return hasParentWith(hasPID(pid));
   }
 
   bool Particle::hasParentWith(const Cut& c) const {
     return hasParentWith([&](const Particle& p){return c->accept(p);});
   }
 
 
   bool Particle::hasAncestor(PdgId pid, bool only_physical) const {
     return hasAncestorWith(hasPID(pid), only_physical);
   }
 
   bool Particle::hasAncestorWith(const Cut& c, bool only_physical) const {
     return hasAncestorWith([&](const Particle& p){return c->accept(p);}, only_physical);
   }
 
 
   bool Particle::hasChildWith(const Cut& c) const {
     return hasChildWith([&](const Particle& p){return c->accept(p);});
   }
 
 
   bool Particle::hasDescendantWith(const Cut& c, bool remove_duplicates) const {
     return hasDescendantWith([&](const Particle& p){return c->accept(p);}, remove_duplicates);
   }
 
   bool Particle::hasStableDescendantWith(const Cut& c) const {
     return hasStableDescendantWith([&](const Particle& p){return c->accept(p);});
   }
 
 
 
   bool Particle::fromBottom() const {
     return hasAncestorWith([](const Particle& p){
         return p.genParticle()->status() == 2 && p.isHadron() && p.hasBottom();
       });
   }
 
   bool Particle::fromCharm() const {
     return hasAncestorWith([](const Particle& p){
         return p.genParticle()->status() == 2 && p.isHadron() && p.hasCharm();
       });
   }
 
   bool Particle::fromHadron() const {
     return hasAncestorWith([](const Particle& p){
         return p.genParticle()->status() == 2 && p.isHadron();
       });
   }
 
   bool Particle::fromTau(bool prompt_taus_only) const {
     if (prompt_taus_only && fromHadron()) return false;
     return hasAncestorWith([](const Particle& p){
         return p.genParticle()->status() == 2 && isTau(p);
       });
   }
 
   bool Particle::fromHadronicTau(bool prompt_taus_only) const {
     return hasAncestorWith([&](const Particle& p){
         return p.genParticle()->status() == 2 && isTau(p) && (!prompt_taus_only || p.isPrompt()) && hasHadronicDecay(p);
       });
   }
 
 
 
 
   bool Particle::isDirect(bool allow_from_direct_tau, bool allow_from_direct_mu) const {
     if (genParticle() == nullptr) return false; // no HepMC connection, give up! Throw UserError exception?
     ConstGenVertexPtr prodVtx = genParticle()->production_vertex();
     if (prodVtx == nullptr) return false; // orphaned particle, has to be assume false
 
     /// @todo Would be nicer to be able to write this recursively up the chain, exiting as soon as a parton or string/cluster is seen
     for (ConstGenParticlePtr ancestor : HepMCUtils::particles(prodVtx, Relatives::ANCESTORS)) {
       const PdgId pid = ancestor->pdg_id();
       if (ancestor->status() != 2) continue; // no non-standard statuses or beams to be used in decision making
-      bool isBeam = false;
-      for(ConstGenParticlePtr b: HepMCUtils::beams(prodVtx->parent_event())){
-        if(ancestor == b){
-          isBeam = true;
-          break;
-        }
-      }
-      if(isBeam) continue; // PYTHIA6 uses status 2 for beams, I think... (sigh)
+
+      std::pair<ConstGenParticlePtr,ConstGenParticlePtr> thebeams = HepMCUtils::beams(prodVtx->parent_event());
+      // PYTHIA6 uses status 2 for beams, I think... (sigh)
+      if ( ancestor == thebeams.first || ancestor == thebeams.second ) continue;
+
       if (PID::isParton(pid)) continue; // PYTHIA6 also uses status 2 for some partons, I think... (sigh)
       if (PID::isHadron(pid)) return false; // direct particles can't be from hadron decays
       if (abs(pid) == PID::TAU && abspid() != PID::TAU && !allow_from_direct_tau) return false; // allow or ban particles from tau decays (permitting tau copies)
       if (abs(pid) == PID::MUON && abspid() != PID::MUON && !allow_from_direct_mu) return false; // allow or ban particles from muon decays (permitting muon copies)
     }
     return true;
   }
 
 
 
   ///////////////////////
 
 
   /// Allow a Particle to be passed to an ostream.
   std::ostream& operator << (std::ostream& os, const Particle& p) {
     string pname;
     try {
       pname = PID::toParticleName(p.pid());
     } catch (...) {
       pname = "PID=" + to_str(p.pid());
     }
     os << "Particle<" << pname << " @ " << p.momentum()/GeV << " GeV>";
     return os;
   }
 
 
   /// Allow ParticlePair to be passed to an ostream.
   std::ostream& operator << (std::ostream& os, const ParticlePair& pp) {
     os << "[" << pp.first << ", " << pp.second << "]";
     return os;
   }
 
 
 
 }
diff --git a/src/Projections/Beam.cc b/src/Projections/Beam.cc
--- a/src/Projections/Beam.cc
+++ b/src/Projections/Beam.cc
@@ -1,147 +1,147 @@
 // -*- C++ -*-
 #include "Rivet/Config/RivetCommon.hh"
 #include "Rivet/Tools/Logging.hh"
 #include "Rivet/Projections/Beam.hh"
 
 namespace Rivet {
 
 
   ParticlePair beams(const Event& e) {
 
     // First try the official way: ask the GenEvent for the beam pointers
     assert(HepMCUtils::particles_size(e.genEvent()) >= 2);
-    vector<ConstGenParticlePtr> beams = HepMCUtils::beams(e.genEvent());
-    if (beams.size() == 2 && beams[0] && beams[1]) {
-      return ParticlePair{beams[0], beams[1]};
+    std::pair<ConstGenParticlePtr,ConstGenParticlePtr> thebeams = HepMCUtils::beams(e.genEvent());
+    if ( thebeams.first && thebeams.second ) {
+      return ParticlePair{thebeams.first, thebeams.second};
     }
     // Ok, that failed: let's find the status = 4 particles by hand
     const vector<Particle> pstat4s = e.allParticles([](const Particle& p){ return p.genParticle()->status() == 4; });
     if (pstat4s.size() >= 2) {
       return ParticlePair{pstat4s[0], pstat4s[1]};
     }
     
     /// There are no barcodes in HepMC3
     /// @todo implement some other fallback rubric?
     #ifndef ENABLE_HEPMC_3
     
     // Hmm, this sucks. Last guess is that barcodes 1 and 2 are the beams
     if (e.genEvent()->barcode_to_particle(1) && e.genEvent()->barcode_to_particle(2)) {
       return ParticlePair{e.genEvent()->barcode_to_particle(1), e.genEvent()->barcode_to_particle(2)};
     }
 
     #endif
 
     // Give up: return null beams
     return ParticlePair{Particle(), Particle()};
   }
 
 
   double sqrtS(const FourMomentum& pa, const FourMomentum& pb) {
     const double mom1 = pa.pz();
     const double e1 = pa.E();
     const double mom2 = pb.pz();
     const double e2 = pb.E();
     double sqrts = sqrt( sqr(e1+e2) - sqr(mom1+mom2) );
     return sqrts;
   }
 
 
   double asqrtS(const FourMomentum& pa, const FourMomentum& pb) {
     const static double MNUCLEON = 939*MeV; //< nominal nucleon mass
     return sqrtS(pa/(pa.mass()/MNUCLEON), pb/(pb.mass()/MNUCLEON));
   }
 
   double asqrtS(const ParticlePair& beams) {
     return sqrtS(beams.first.mom()/nuclA(beams.first), beams.second.mom()/nuclA(beams.second));
   }
 
 
   FourMomentum acmsBoostVec(const FourMomentum& pa, const FourMomentum& pb) {
     const static double MNUCLEON = 939*MeV; //< nominal nucleon mass
     const double na = pa.mass()/MNUCLEON, nb = pb.mass()/MNUCLEON;
     return cmsBoostVec(pa/na, pb/nb);
   }
 
   FourMomentum acmsBoostVec(const ParticlePair& beams) {
     return cmsBoostVec(beams.first.mom()/nuclA(beams.first), beams.second.mom()/nuclA(beams.second));
   }
 
 
   Vector3 cmsBetaVec(const FourMomentum& pa, const FourMomentum& pb) {
     // const Vector3 rtn = (pa.p3() + pb.p3()) / (pa.E() + pb.E());
     const Vector3 rtn = (pa + pb).betaVec();
     return rtn;
   }
 
   Vector3 acmsBetaVec(const FourMomentum& pa, const FourMomentum& pb) {
     const static double MNUCLEON = 939*MeV; //< nominal nucleon mass
     const Vector3 rtn = cmsBetaVec(pa/(pa.mass()/MNUCLEON), pb/(pb.mass()/MNUCLEON));
     return rtn;
   }
 
   Vector3 acmsBetaVec(const ParticlePair& beams) {
     const Vector3 rtn = cmsBetaVec(beams.first.mom()/nuclA(beams.first), beams.second.mom()/nuclA(beams.second));
     return rtn;
   }
 
 
   Vector3 cmsGammaVec(const FourMomentum& pa, const FourMomentum& pb) {
     // const Vector3 rtn = (pa + pb).gammaVec();
     const double gamma = (pa.E() + pb.E()) / sqrt( sqr(pa.mass()) + sqr(pb.mass()) + 2*(pa.E()*pb.E() - dot(pa.p3(), pb.p3())) );
     const Vector3 rtn = gamma * (pa.p3() + pb.p3()).unit();
     return rtn;
   }
 
   Vector3 acmsGammaVec(const FourMomentum& pa, const FourMomentum& pb) {
     const static double MNUCLEON = 939*MeV; //< nominal nucleon mass
     Vector3 rtn = cmsGammaVec(pa/(pa.mass()/MNUCLEON), pb/(pb.mass()/MNUCLEON));
     return rtn;
   }
 
   Vector3 acmsGammaVec(const ParticlePair& beams) {
     Vector3 rtn = cmsGammaVec(beams.first.mom()/nuclA(beams.first), beams.second.mom()/nuclA(beams.second));
     return rtn;
   }
 
 
   LorentzTransform cmsTransform(const FourMomentum& pa, const FourMomentum& pb) {
     /// @todo Automatically choose to construct from beta or gamma according to which is more precise?
     return LorentzTransform::mkFrameTransformFromGamma(cmsGammaVec(pa, pb));
   }
 
   LorentzTransform acmsTransform(const FourMomentum& pa, const FourMomentum& pb) {
     /// @todo Automatically choose to construct from beta or gamma according to which is more precise?
     return LorentzTransform::mkFrameTransformFromGamma(acmsGammaVec(pa, pb));
   }
 
   LorentzTransform acmsTransform(const ParticlePair& beams) {
     return LorentzTransform::mkFrameTransformFromGamma(acmsGammaVec(beams));
   }
 
 
 
   /////////////////////////////////////////////
 
 
 
   void Beam::project(const Event& e) {
     _theBeams = Rivet::beams(e);
     MSG_DEBUG("Beam particles = " << _theBeams << " => sqrt(s) = " << sqrtS()/GeV << " GeV");
   }
 
 
   FourVector Beam::pv() const {
     RivetHepMC::FourVector v1, v2;
     const ParticlePair bpair = beams();
     if (bpair.first.genParticle() && bpair.first.genParticle()->end_vertex())
       v1 = bpair.first.genParticle()->end_vertex()->position();
     if (bpair.second.genParticle() && bpair.second.genParticle()->end_vertex())
       v2 = bpair.second.genParticle()->end_vertex()->position();
     const FourVector rtn = (v1 == v2) ? FourVector(v1.t(), v1.x(), v1.y(), v1.z()) : FourVector();
     MSG_DEBUG("Beam PV 4-position = " << rtn);
     return rtn;
   }
 
 
 
 }
diff --git a/src/Tools/RivetHepMC_2.cc b/src/Tools/RivetHepMC_2.cc
--- a/src/Tools/RivetHepMC_2.cc
+++ b/src/Tools/RivetHepMC_2.cc
@@ -1,117 +1,116 @@
 // -*- C++ -*-
 
 #include "Rivet/Tools/RivetHepMC.hh"
 
 namespace Rivet{
   
   const Relatives Relatives::PARENTS     = HepMC::parents;
   const Relatives Relatives::CHILDREN    = HepMC::children;
   const Relatives Relatives::ANCESTORS   = HepMC::ancestors;
   const Relatives Relatives::DESCENDANTS = HepMC::descendants;
   
   namespace HepMCUtils{
     
     std::vector<ConstGenParticlePtr> particles(ConstGenEventPtr ge){
       std::vector<ConstGenParticlePtr> result;
       for(GenEvent::particle_const_iterator pi = ge->particles_begin(); pi != ge->particles_end(); ++pi){
         result.push_back(*pi);
       }
       return result;
     }
     
     std::vector<ConstGenParticlePtr> particles(const GenEvent *ge){
       std::vector<ConstGenParticlePtr> result;
       for(GenEvent::particle_const_iterator pi = ge->particles_begin(); pi != ge->particles_end(); ++pi){
         result.push_back(*pi);
       }
       return result;
     }
    
     std::vector<ConstGenVertexPtr> vertices(ConstGenEventPtr ge){
       std::vector<ConstGenVertexPtr> result;
       for(GenEvent::vertex_const_iterator vi = ge->vertices_begin(); vi != ge->vertices_end(); ++vi){
         result.push_back(*vi);
       }
       return result;
     }
     
     std::vector<ConstGenVertexPtr> vertices(const GenEvent *ge){
       std::vector<ConstGenVertexPtr> result;
       for(GenEvent::vertex_const_iterator vi = ge->vertices_begin(); vi != ge->vertices_end(); ++vi){
         result.push_back(*vi);
       }
       return result;
     }
     
     std::vector<ConstGenParticlePtr> particles(ConstGenVertexPtr gv, const Relatives &relo){
       std::vector<ConstGenParticlePtr> result;
       /// @todo A particle_const_iterator on GenVertex would be nice...
       // Before HepMC 2.7.0 there were no GV::particles_const_iterators and constness consistency was all screwed up :-/
 #if HEPMC_VERSION_CODE >= 2007000
       for (HepMC::GenVertex::particle_iterator pi = gv->particles_begin(relo); pi != gv->particles_end(relo); ++pi)
       result.push_back(*pi);
 #else
       HepMC::GenVertex* gv2 = const_cast<HepMC::GenVertex*>(gv);
       for (HepMC::GenVertex::particle_iterator pi = gv2->particles_begin(relo); pi != gv2->particles_end(relo); ++pi)
       result.push_back(const_cast<ConstGenParticlePtr>(*pi));
 #endif
       return result;
     }
 
     std::vector<ConstGenParticlePtr> particles(ConstGenParticlePtr gp, const Relatives &relo){
       ConstGenVertexPtr vtx;
       
       switch(relo){
         case HepMC::parents:
         
         case HepMC::ancestors:
           vtx = gp->production_vertex();
         break;
         
         case HepMC::children:
         
         case HepMC::descendants:
           vtx = gp->end_vertex();
         break;
         
         default:
         
         throw std::runtime_error("Not implemented!");
         break;
       }
       
       return particles(vtx, relo);
     }
 
     
     
     int uniqueId(ConstGenParticlePtr gp){
       return gp->barcode();
     }
     
     int particles_size(ConstGenEventPtr ge){
       return ge->particles_size();
     }
 
     int particles_size(const GenEvent *ge){
       return ge->particles_size();
     }
     
-    std::vector<ConstGenParticlePtr> beams(const GenEvent *ge){
-      pair<HepMC::GenParticle*, HepMC::GenParticle*> beams = ge->beam_particles();
-      return std::vector<ConstGenParticlePtr>{beams.first, beams.second};
+    std::pair<ConstGenParticlePtr,ConstGenParticlePtr> beams(const GenEvent *ge){
+      return ge->beam_particles();
     }
 
     std::shared_ptr<HepMC::IO_GenEvent> makeReader(std::istream &istr,
                                                    std::string *) {
       return make_shared<HepMC::IO_GenEvent>(istr);
     }
    
     bool readEvent(std::shared_ptr<HepMC::IO_GenEvent> io, std::shared_ptr<GenEvent> evt){
       if(io->rdstate() != 0) return false;
       if(!io->fill_next_event(evt.get())) return false;
       return true;
     }
     
   }
 }
diff --git a/src/Tools/RivetHepMC_3.cc b/src/Tools/RivetHepMC_3.cc
--- a/src/Tools/RivetHepMC_3.cc
+++ b/src/Tools/RivetHepMC_3.cc
@@ -1,106 +1,108 @@
 // -*- C++ -*-
 
 #include "Rivet/Tools/RivetHepMC.hh"
 #include "HepMC3/ReaderAscii.h"
 #include "HepMC3/ReaderAsciiHepMC2.h"
 
 namespace Rivet{
   
   namespace HepMCUtils{
     
     std::vector<ConstGenParticlePtr> particles(ConstGenEventPtr ge){
       return ge->particles();
     }
     
     std::vector<ConstGenParticlePtr> particles(const GenEvent *ge){
       assert(ge);
       return ge->particles();
     }
     
     std::vector<ConstGenVertexPtr> vertices(ConstGenEventPtr ge){
       return ge->vertices();
     }
     
     std::vector<ConstGenVertexPtr> vertices(const GenEvent *ge){
       assert(ge);
       return ge->vertices();
     }
     
     std::vector<ConstGenParticlePtr> particles(ConstGenVertexPtr gv, const Relatives &relo){
       return relo(gv);
     }
     
     std::vector<ConstGenParticlePtr> particles(ConstGenParticlePtr gp, const Relatives &relo){
       return relo(gp);
     }
     
     int particles_size(ConstGenEventPtr ge){
       return particles(ge).size();
     }
 
     int particles_size(const GenEvent *ge){
       return particles(ge).size();
     }
     
     int uniqueId(ConstGenParticlePtr gp){
       return gp->id();
     }
     
-    std::vector<ConstGenParticlePtr> beams(const GenEvent *ge){
-      return ge->beams();
+    std::pair<ConstGenParticlePtr,ConstGenParticlePtr> beams(const GenEvent *ge) {
+      std::vector<ConstGenParticlePtr> beamlist = ge->beams();
+      if ( beamlist.size() < 2 ) return std::pair<ConstGenParticlePtr,ConstGenParticlePtr>();
+      return std::make_pair(beamlist[0], beamlist[1]);
     }
     
     bool readEvent(std::shared_ptr<HepMC_IO_type> io, std::shared_ptr<GenEvent> evt){
       io->read_event(*evt);
       return !io->failed();
     }
     
     shared_ptr<HepMC_IO_type> makeReader(std::istream & istr,
                                          std::string * errm) {
       shared_ptr<HepMC_IO_type> ret;
       
       // First scan forward and check if there is some hint as to what
       // kind of file we are looking att.
       int nchar = 256;
       std::string header;
       while ( nchar-- && !istr.eof() )
         header += char(istr.get());
       // If this stream was too short to contain any reasonable number
       // of events, just give up.
       if ( !istr ) {
         if ( errm ) *errm = "Could not find HepMC header information";
         return shared_ptr<HepMC_IO_type>();
       }
 
       // Now reset the stream to its original state ...
       for ( int i = header.length() - 1; i >= 0; --i )
         istr.putback(header[i]);
 
       // ... and check which kind of format it was and create the
       // corresponding reader. First try the HepM3 ascii format.
       if ( header.find("HepMC::Asciiv3-START_EVENT_LISTING") !=
            std::string::npos )
         ret = make_shared<RivetHepMC::ReaderAscii>(istr);
 
       // Check if the file is written by WriterRoot or WriterRootTree.
       else if ( header.substr(0, 4) == "root" ) {
         if ( errm ) *errm = "Rivet cancurrently not read HepMC root files.";
         return ret;
       }
 
       // The default is to assume it is a good old HepMC2 ascii file.
       else {
         ret = make_shared<RivetHepMC::ReaderAsciiHepMC2>(istr);
       }
 
       // Check that everything was ok.
       if ( ret->failed() ) {
         if ( errm ) *errm = "Problems reading from HepMC file.";
         ret = shared_ptr<HepMC_IO_type>();
       }
 
       return ret;
     }
 
   }
 }
diff --git a/test/testApi.cc b/test/testApi.cc
--- a/test/testApi.cc
+++ b/test/testApi.cc
@@ -1,35 +1,31 @@
 #include "Rivet/AnalysisHandler.hh"
 #include "HepMC/GenEvent.h"
 #include "Rivet/Tools/RivetHepMC.hh"
 
 using namespace std;
 
 int main() {
   Rivet::AnalysisHandler ah;
   Rivet::Log::setLevel("Rivet", Rivet::Log::TRACE);
 
   // Specify the analyses to be used
   ah.addAnalysis("EXAMPLE");
   ah.addAnalyses({{ "MC_JETS", "EXAMPLE_CUTS", "EXAMPLE_SMEAR" }});
 
   std::ifstream file("testApi.hepmc");
   shared_ptr<Rivet::HepMC_IO_type> reader = Rivet::HepMCUtils::makeReader(file);
-  std::shared_ptr<Rivet::GenEvent> evt;
-  Rivet::HepMCUtils::readEvent(reader, evt);
+  std::shared_ptr<Rivet::GenEvent> evt = make_shared<Rivet::GenEvent>();
   double sum_of_weights = 0.0;
-  while (evt) {
+  while ( Rivet::HepMCUtils::readEvent(reader, evt) ) {
     // Analyse current event
     ah.analyze(*evt);
     sum_of_weights += evt->weights()[0];
-
-    // Clean up and get next event
-    Rivet::HepMCUtils::readEvent(reader, evt);
   }
   file.close();
 
   ah.setCrossSection(1.0);
   ah.finalize();
   ah.writeData("out.yoda");
 
   return 0;
 }
diff --git a/test/testNaN.cc b/test/testNaN.cc
--- a/test/testNaN.cc
+++ b/test/testNaN.cc
@@ -1,76 +1,72 @@
 #include "Rivet/AnalysisHandler.hh"
 #include "Rivet/Analysis.hh"
 #include "Rivet/Tools/RivetYODA.hh"
 #include <limits>
 #include <cmath>
 #include <iostream>
 
 using namespace std;
 
 
 class Test : public Rivet::Analysis {
 public:
   Test() : Analysis("Test") {}
 
   void init() {
     _h_test = bookHisto1D("test", 50, 66.0, 116.0);
   }
 
   void analyze(const Rivet::Event & e) {
     cout << "Normal fill" << endl;
     _h_test->fill(90., 1.);
 
     cout << "Underflow fill" << endl;
     _h_test->fill(30.,1.);
 
     cout << "Overflow fill" << endl;
     _h_test->fill(130.,1.);
 
      cout << "Inf fill" << endl;
     try {
       _h_test->fill(numeric_limits<double>::infinity(), 1.);
     } catch (YODA::RangeError e) {
       cerr << e.what() << '\n';
       if ( string(e.what()) != string("X is Inf") ) throw;
     }
 
     cout << "NaN fill" << endl;
     try {
       _h_test->fill(numeric_limits<double>::quiet_NaN(), 1.);
     } catch (YODA::RangeError e) {
       cerr << e.what() << '\n';
       if ( string(e.what()) != string("X is NaN") ) throw;
     }
   }
 
 private:
   Rivet::Histo1DPtr _h_test;
 };
 
 DECLARE_RIVET_PLUGIN(Test);
 
 int main() {
   Rivet::AnalysisHandler rivet;
   rivet.addAnalysis("Test");
 
   std::ifstream file("testApi.hepmc");
   shared_ptr<Rivet::HepMC_IO_type> reader = Rivet::HepMCUtils::makeReader(file);
-  std::shared_ptr<Rivet::GenEvent> evt;
-  Rivet::HepMCUtils::readEvent(reader, evt);
+  std::shared_ptr<Rivet::GenEvent> evt = make_shared<Rivet::GenEvent>();
   double sum_of_weights = 0.0;
-  while (evt) {
+  while ( Rivet::HepMCUtils::readEvent(reader, evt) ) {
     // Analyse current event
     rivet.analyze(*evt);
     sum_of_weights += evt->weights()[0];
-
-    // Clean up and get next event
-    Rivet::HepMCUtils::readEvent(reader, evt);
   }
-   file.close();
+  file.close();
 
   rivet.setCrossSection(1.0);
   rivet.finalize();
   rivet.writeData("NaN.aida");
 
   return 0;
 }