diff --git a/include/Rivet/Projections/LeadingParticlesFinalState.hh b/include/Rivet/Projections/LeadingParticlesFinalState.hh
--- a/include/Rivet/Projections/LeadingParticlesFinalState.hh
+++ b/include/Rivet/Projections/LeadingParticlesFinalState.hh
@@ -1,76 +1,83 @@
 // -*- C++ -*-
 #ifndef RIVET_LeadingParticlesFinalState_HH
 #define RIVET_LeadingParticlesFinalState_HH
 
 #include "Rivet/Event.hh"
 #include "Rivet/Projection.hh"
 #include "Rivet/Projections/FinalState.hh"
 
 namespace Rivet {
 
 
   /// @brief Get the highest-pT occurrences of FS particles with the specified PDG IDs.
   class LeadingParticlesFinalState : public FinalState {
-
   public:
 
     /// Constructor: the supplied FinalState projection is assumed to live through the run.
     LeadingParticlesFinalState(const FinalState& fsp)
       :  FinalState(), _leading_only(false)
     {
       setName("LeadingParticlesFinalState");
       addProjection(fsp, "FS");
     }
 
     /// Clone on the heap.
     DEFAULT_RIVET_PROJ_CLONE(LeadingParticlesFinalState);
 
+
     /// Add a particle ID to the list of leading particles selected
-    LeadingParticlesFinalState& addParticleId(long id) {
+    LeadingParticlesFinalState& addParticleId(PdgId id) {
       _ids.insert(id);
       return *this;
     }
 
     /// Add a particle ID to the list of leading particles selected
-    LeadingParticlesFinalState& addParticleIdPair(long id) {
+    LeadingParticlesFinalState& addParticleIds(vector<PdgId> ids) {
+      for (PdgId id : ids) _ids.insert(id);
+      return *this;
+    }
+
+    /// Add a particle ID to the list of leading particles selected
+    LeadingParticlesFinalState& addParticleIdPair(PdgId id) {
       _ids.insert(id);
       _ids.insert(-id);
       return *this;
     }
 
     /// Toggle whether to keep track only of the leading particle of any ID,
     /// or the leading particle of all IDs separately
     /// Default is the latter (=false)
-    void setLeadingOnly(const bool& leadingonly) {
+    void setLeadingOnly(bool leadingonly) {
       _leading_only = leadingonly;
     }
 
     // /// Check if a particle of a particular ID was found in the current event
     // bool hasParticleId(const PdgId pid) const;
 
     // /// Get a particle of a particular ID (check it exists first)
     // bool get(const PdgId pid) const;
 
 
   protected:
 
     /// Apply the projection on the supplied event.
     void project(const Event& e);
 
     /// Compare projections.
     int compare(const Projection& p) const;
 
-    /// Check if the particle's ID is in the list
-    bool inList(const Particle& particle) const;
 
   private:
 
+    /// Check if the particle's ID is in the list
+    bool _inList(const Particle& particle) const;
+
     /// IDs of the leading particles to be selected
     std::set<long>_ids;
     bool _leading_only;
 
   };
 
 }
 
 #endif
diff --git a/src/Projections/InitialQuarks.cc b/src/Projections/InitialQuarks.cc
--- a/src/Projections/InitialQuarks.cc
+++ b/src/Projections/InitialQuarks.cc
@@ -1,65 +1,65 @@
 // -*- C++ -*-
 #include "Rivet/Projections/InitialQuarks.hh"
 
 
 #define IS_PARTON_PDGID(id) ( abs(id) <= 100 && abs(id) != 22 && (abs(id) < 11 || abs(id) > 18) )
 
 
 namespace Rivet {
 
 
   int InitialQuarks::compare(const Projection& p) const {
     return EQUIVALENT;
   }
 
 
   void InitialQuarks::project(const Event& e) {
     _theParticles.clear();
 
-    for (const GenParticle* p : Rivet::particles(e.genEvent())) {
-      const GenVertex* pv = p->production_vertex();
-      const GenVertex* dv = p->end_vertex();
+    for (const GenParticlePtr p : Rivet::particles(e.genEvent())) {
+      const GenVertexPtr pv = p->production_vertex();
+      const GenVertexPtr dv = p->end_vertex();
       const PdgId pid = abs(p->pdg_id());
       bool passed = inRange((long)pid, 1, 6);
       if (passed) {
         if (pv != 0) {
-          for (const GenParticle* pp : particles_in(pv)) {
+          for (const GenParticlePtr pp : particles_in(pv)) {
             // Only accept if parent is electron or Z0
             const PdgId pid = abs(pp->pdg_id());
             passed = (pid == PID::ELECTRON || abs(pp->pdg_id()) == PID::ZBOSON || abs(pp->pdg_id()) == PID::GAMMA);
           }
         } else {
           passed = false;
         }
       }
 
       if (getLog().isActive(Log::TRACE)) {
         const int st = p->status();
         const double pT = p->momentum().perp();
         const double eta = p->momentum().eta();
         MSG_TRACE(std::boolalpha
                   << "ID = " << p->pdg_id() << ", status = " << st << ", pT = " << pT
                   << ", eta = " << eta << ": result = " << passed);
         if (pv != 0) {
-          for (const GenParticle* pp : particles_in(pv)) {
+          for (const GenParticlePtr pp : particles_in(pv)) {
             MSG_TRACE(std::boolalpha << " parent ID = " << pp->pdg_id());
           }
         }
         if (dv != 0) {
-          for (const GenParticle* pp : particles_out(dv)) {
+          for (const GenParticlePtr pp : particles_out(dv)) {
             MSG_TRACE(std::boolalpha << " child ID  = " << pp->pdg_id());
           }
         }
       }
-      if (passed) _theParticles.push_back(Particle(*p));
+      if (passed) _theParticles.push_back(Particle(p));
     }
     MSG_DEBUG("Number of initial quarks = " << _theParticles.size());
     if (!_theParticles.empty()) {
       for (size_t i = 0; i < _theParticles.size(); i++) {
         MSG_DEBUG("Initial quark[" << i << "] = " << _theParticles[i].pid());
       }
     }
   }
 
 
 }
diff --git a/src/Projections/InvMassFinalState.cc b/src/Projections/InvMassFinalState.cc
--- a/src/Projections/InvMassFinalState.cc
+++ b/src/Projections/InvMassFinalState.cc
@@ -1,190 +1,190 @@
 // -*- C++ -*-
 #include "Rivet/Projections/InvMassFinalState.hh"
 
 namespace Rivet {
 
 
   InvMassFinalState::InvMassFinalState(const FinalState& fsp,
                                        const pair<PdgId, PdgId>& idpair, // pair of decay products
                                        double minmass, // min inv mass
                                        double maxmass, // max inv mass
                                        double masstarget)
     : _minmass(minmass), _maxmass(maxmass), _masstarget(masstarget), _useTransverseMass(false)
   {
     setName("InvMassFinalState");
     addProjection(fsp, "FS");
     _decayids.push_back(idpair);
   }
 
 
   InvMassFinalState::InvMassFinalState(const FinalState& fsp,
                                        const vector<pair<PdgId, PdgId> >& idpairs,  // vector of pairs of decay products
                                        double minmass, // min inv mass
                                        double maxmass, // max inv mass
                                        double masstarget)
     : _decayids(idpairs), _minmass(minmass), _maxmass(maxmass), _masstarget(masstarget), _useTransverseMass(false)
   {
     setName("InvMassFinalState");
     addProjection(fsp, "FS");
   }
 
 
   InvMassFinalState::InvMassFinalState(const pair<PdgId, PdgId>& idpair, // pair of decay products
                                        double minmass, // min inv mass
                                        double maxmass, // max inv mass
                                        double masstarget)
     : _minmass(minmass), _maxmass(maxmass), _masstarget(masstarget), _useTransverseMass(false)
   {
     setName("InvMassFinalState");
     _decayids.push_back(idpair);
   }
 
 
   InvMassFinalState::InvMassFinalState(const vector<pair<PdgId, PdgId> >& idpairs,  // vector of pairs of decay products
                                        double minmass, // min inv mass
                                        double maxmass, // max inv mass
                                        double masstarget)
     : _decayids(idpairs), _minmass(minmass), _maxmass(maxmass), _masstarget(masstarget), _useTransverseMass(false)
   {
     setName("InvMassFinalState");
   }
 
 
   int InvMassFinalState::compare(const Projection& p) const {
     // First compare the final states we are running on
     int fscmp = mkNamedPCmp(p, "FS");
     if (fscmp != EQUIVALENT) return fscmp;
 
     // Then compare the two as final states
     const InvMassFinalState& other = dynamic_cast<const InvMassFinalState&>(p);
     fscmp = FinalState::compare(other);
     if (fscmp != EQUIVALENT) return fscmp;
 
     // Compare the mass limits
     int masstypecmp = cmp(_useTransverseMass, other._useTransverseMass);
     if (masstypecmp != EQUIVALENT) return masstypecmp;
     int massllimcmp = cmp(_minmass, other._minmass);
     if (massllimcmp != EQUIVALENT) return massllimcmp;
     int masshlimcmp = cmp(_maxmass, other._maxmass);
     if (masshlimcmp != EQUIVALENT) return masshlimcmp;
 
     // Compare the decay species
     int decaycmp = cmp(_decayids, other._decayids);
     if (decaycmp != EQUIVALENT) return decaycmp;
 
     // Finally compare them as final states
     return FinalState::compare(other);
   }
 
 
   void InvMassFinalState::project(const Event& e) {
     const FinalState& fs = applyProjection<FinalState>(e, "FS");
     calc(fs.particles());
   }
 
 
   void InvMassFinalState::calc(const Particles& inparticles) {
     _theParticles.clear();
     _particlePairs.clear();
 
     // Containers for the particles of type specified in the pair
     vector<const Particle*> type1;
     vector<const Particle*> type2;
     // Get all the particles of the type specified in the pair from the particle list
     foreach (const Particle& ipart, inparticles) {
       // Loop around possible particle pairs (typedef needed to keep BOOST_FOREACH happy)
       foreach (const PdgIdPair& ipair, _decayids) {
         if (ipart.pid() == ipair.first) {
           if (accept(ipart)) {
             type1 += &ipart;
           }
         } else if (ipart.pid() == ipair.second) {
           if (accept(ipart)) {
             type2 += &ipart;
           }
         }
       }
     }
     if (type1.empty() || type2.empty()) return;
 
     // Temporary container of selected particles iterators
     // Useful to compare iterators and avoid double occurrences of the same
     // particle in case it matches with more than another particle
     vector<const Particle*> tmp;
 
     // Now calculate the inv mass
     pair<double, pair<Particle, Particle> > closestPair;
     closestPair.first = 1e30;
     foreach (const Particle* i1, type1) {
       foreach (const Particle* i2, type2) {
         // Check this is actually a pair
         // (if more than one pair in vector particles can be unrelated)
         bool found = false;
         foreach (const PdgIdPair& ipair, _decayids) {
           if (i1->pid() == ipair.first && i2->pid() == ipair.second) {
             found = true;
             break;
           }
         }
         if (!found) continue;
 
         FourMomentum v4 = i1->momentum() + i2->momentum();
         if (v4.mass2() < 0) {
           MSG_DEBUG("Constructed negative inv mass2: skipping!");
           continue;
         }
         bool passedMassCut = false;
         if (_useTransverseMass) {
           passedMassCut = inRange(massT(i1->momentum(), i2->momentum()), _minmass, _maxmass);
         } else {
           passedMassCut = inRange(v4.mass(), _minmass, _maxmass);
         }
 
         if (passedMassCut) {
           MSG_DEBUG("Selecting particles with IDs " << i1->pid() << " & " << i2->pid()
                     << " and mass = " << v4.mass()/GeV << " GeV");
           // Store accepted particles, avoiding duplicates
           if (find(tmp.begin(), tmp.end(), i1) == tmp.end()) {
             tmp.push_back(i1);
             _theParticles += *i1;
           }
           if (find(tmp.begin(), tmp.end(), i2) == tmp.end()) {
             tmp.push_back(i2);
             _theParticles += *i2;
           }
           // Store accepted particle pairs
           _particlePairs += make_pair(*i1, *i2);
           if (_masstarget>0.0) {
             double diff=fabs(v4.mass()-_masstarget);
             if (diff<closestPair.first) {
               closestPair.first=diff;
               closestPair.second=make_pair(*i1, *i2);
             }
           }
         }
       }
     }
     if (_masstarget > 0.0 && closestPair.first < 1e30) {
       _theParticles.clear();
       _particlePairs.clear();
       _theParticles += closestPair.second.first;
       _theParticles += closestPair.second.second;
       _particlePairs += closestPair.second;
     }
 
     MSG_DEBUG("Selected " << _theParticles.size() << " particles " << "(" << _particlePairs.size() << " pairs)");
     if (getLog().isActive(Log::TRACE)) {
       foreach (const Particle& p, _theParticles) {
-        MSG_TRACE("ID: " << p.pid() << ", barcode: " << p.genParticle()->barcode());
+        MSG_TRACE("ID: " << p.pid() << ", barcode: " << p.genParticle()->id());
       }
     }
   }
 
 
   /// Constituent pairs
   const std::vector<std::pair<Particle, Particle> >& InvMassFinalState::particlePairs() const {
     return _particlePairs;
   }
 
 
 }
diff --git a/src/Projections/LeadingParticlesFinalState.cc b/src/Projections/LeadingParticlesFinalState.cc
--- a/src/Projections/LeadingParticlesFinalState.cc
+++ b/src/Projections/LeadingParticlesFinalState.cc
@@ -1,85 +1,65 @@
 #include "Rivet/Projections/LeadingParticlesFinalState.hh"
 #include "Rivet/Particle.hh"
 
 namespace Rivet {
 
 
   int LeadingParticlesFinalState::compare(const Projection& p) const {
     // First compare the final states we are running on
     int fscmp = mkNamedPCmp(p, "FS");
     if (fscmp != EQUIVALENT) return fscmp;
 
     // Then compare the two as final states
     const LeadingParticlesFinalState& other = dynamic_cast<const LeadingParticlesFinalState&>(p);
     fscmp = FinalState::compare(other);
     if (fscmp != EQUIVALENT) return fscmp;
 
     int locmp = cmp(_leading_only, other._leading_only);
     if (locmp != EQUIVALENT) return locmp;
 
     // Finally compare the IDs
     if (_ids < other._ids) return ORDERED;
     else if (other._ids < _ids) return UNORDERED;
     return EQUIVALENT;
   }
 
 
   void LeadingParticlesFinalState::project(const Event & e) {
     _theParticles.clear();
     const FinalState& fs = applyProjection<FinalState>(e, "FS");
-
-    // Temporary container
-    map<long, Particles::const_iterator> tmp;
-
     const Particles& particles = fs.particles();
     MSG_DEBUG("Original final state particles size " << particles.size());
-    Particles::const_iterator ifs;
-    for (ifs = particles.begin(); ifs != particles.end(); ++ifs) {
-      if (inList(*ifs) && FinalState::accept(*ifs->genParticle())) {
+
+    map<PdgId, const Particle*> tmp;
+    for (const Particle& p : particles) {
+      const PdgId pid = p.pid();
+      // If it's a PID we're looking for, and passes the cuts
+      if (_ids.find(pid) != _ids.end() && FinalState::accept(p.genParticle())) {
         // Look for an existing particle in tmp container
-        map < long, Particles::const_iterator >::const_iterator itmp = tmp.find(ifs->pid());
-        if (itmp != tmp.end()) {  // if a particle with this type has been already selected
+        if (tmp.find(pid) != tmp.end()) { // if a particle with this type has been already selected
+          const Particle& p2 = *tmp.find(pid)->second;
           // If the new pT is higher than the previous one, then substitute...
-          if (ifs->pT() > itmp->second->pT()) {
-            tmp[ifs->pid()] = ifs;
-          }
+          if (p.pT() > p2.pT()) tmp[pid] = &p;
+        } else {
           // ...otherwise insert in the container
-        } else {
-          tmp[ifs->pid()] = ifs;
+          tmp[pid] = &p;
         }
       }
     }
 
     // Loop on the tmp container and fill _theParticles
-    map<long, Particles::const_iterator>::const_iterator i;
-    for (i = tmp.begin(); i != tmp.end(); ++i) {
-      MSG_DEBUG("LeadingParticlesFinalState is accepting particle ID " << i->second->pid()
-               << " with momentum " << i->second->momentum());
-      _theParticles.push_back(*(i->second));
+    for (const pair<PdgId,const Particle*>& id_p : tmp) {
+      MSG_DEBUG("Accepting particle ID " << id_p.first << " with momentum " << id_p.second->momentum());
+      _theParticles.push_back(*id_p.second);
     }
 
-    if (_leading_only) {
-      double ptmax=0.0;
-      Particle pmax;
-
-      foreach (const Particle& p, _theParticles) {
-        if (p.pT() > ptmax) {
-          ptmax = p.pT();
-          pmax = p;
-        }
-      }
-
+    if (_leading_only && !_theParticles.empty()) {
+      sortByPt(_theParticles);
+      const Particle pmax = _theParticles.front();
       _theParticles.clear();
       _theParticles.push_back(pmax);
     }
   }
 
 
-  bool LeadingParticlesFinalState::inList(const Particle & particle) const {
-    std::set < long >::const_iterator ilist = _ids.find(particle.pid());
-    if (ilist != _ids.end()) return true;
-    return false;
-  }
-
-
 }