diff --git a/ColliderBit/include/gambit/ColliderBit/Utils.hpp b/ColliderBit/include/gambit/ColliderBit/Utils.hpp
index da112d6803..7b0341ca7c 100644
--- a/ColliderBit/include/gambit/ColliderBit/Utils.hpp
+++ b/ColliderBit/include/gambit/ColliderBit/Utils.hpp
@@ -52,6 +52,20 @@ namespace Gambit
     /// Unit conversions (multiply to construct in standard units, divide to decode to that unit)
     static const double GeV = 1, MeV = 1e-3, TeV = 1e3;
 
+    /// Struct of different jet collection settings
+    struct jet_collection_settings
+    {
+      std::string key;
+      std::string algorithm;
+      double R;
+      std::string recombination_scheme;
+      std::string strategy;
+    };
+
+    /// Storage of different FastJet methods
+    FJNS::JetAlgorithm FJalgorithm_map(std::string);
+    FJNS::Strategy FJstrategy_map(std::string);
+    FJNS::RecombinationScheme FJRecomScheme_map(std::string);
 
     /// Use the HEPUtils Event without needing namespace qualification
     using HEPUtils::Event;
@@ -305,10 +319,10 @@ namespace Gambit
 
 
     /// Check if there's a physics object above ptmin in an annulus rmin..rmax around the given four-momentum p4
-    inline bool object_in_cone(const HEPUtils::Event& e, const HEPUtils::P4& p4, double ptmin, double rmax, double rmin=0.05) {
+    inline bool object_in_cone(const HEPUtils::Event& e, std::string jetcollection, const HEPUtils::P4& p4, double ptmin, double rmax, double rmin=0.05) {
       for (const HEPUtils::Particle* p : e.visible_particles())
         if (p->pT() > ptmin && HEPUtils::in_range(HEPUtils::deltaR_eta(p4, *p), rmin, rmax)) return true;
-      for (const HEPUtils::Jet* j : e.jets())
+      for (const HEPUtils::Jet* j : e.jets(jetcollection))
         if (j->pT() > ptmin && HEPUtils::in_range(HEPUtils::deltaR_eta(p4, *j), rmin, rmax)) return true;
       return false;
     }
diff --git a/ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp b/ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp
index 08b8588a30..d29690567a 100644
--- a/ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp
+++ b/ColliderBit/include/gambit/ColliderBit/colliders/BaseCollider.hpp
@@ -21,6 +21,8 @@
 #include <string>
 #include <vector>
 
+#include "gambit/ColliderBit/Utils.hpp"
+
 namespace Gambit
 {
   namespace ColliderBit
@@ -33,7 +35,7 @@ namespace Gambit
       public:
 
         /// Constructor
-        BaseCollider() : partonOnly(false), antiktR(0.4) {}
+        BaseCollider() : partonOnly(false), all_jet_collection_settings({}), jetcollection_taus("") {}
         /// Destructor
         virtual ~BaseCollider() {}
         /// Reset this instance for reuse, avoiding the need for "new" or "delete".
@@ -67,8 +69,12 @@ namespace Gambit
 
         /// Flag indicating if events from this collider should be processed as parton-only or full events
         bool partonOnly;
-        ///The jet radius used for the anti-kt jet clustering.
-        double antiktR;
+
+        /// Vector of different jet collection settings
+        std::vector<jet_collection_settings> all_jet_collection_settings;
+
+        /// Key for jet collection used in adding taus
+        std::string jetcollection_taus;
 
     };
 
diff --git a/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp b/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp
index 47d9daf963..3a07adfd74 100644
--- a/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp
+++ b/ColliderBit/include/gambit/ColliderBit/colliders/Pythia8/Py8EventConversions.hpp
@@ -14,12 +14,14 @@
 ///  \author Pat Scott
 ///  \author Martin White
 ///  \author Are Raklev  June 2021
+///  \author Chris Chang Nov 2023
 ///
 ///  *********************************************
 
 #pragma once
 
 #include "gambit/ColliderBit/colliders/EventConversionUtils.hpp"
+#include "gambit/ColliderBit/colliders/BaseCollider.hpp"
 
 //#define COLLIDERBIT_DEBUG
 
@@ -35,7 +37,7 @@ namespace Gambit
     /// Convert a hadron-level EventT into an unsmeared HEPUtils::Event
     /// @todo Overlap between jets and prompt containers: need some isolation in MET calculation
     template<typename EventT>
-    void convertParticleEvent(const EventT& pevt, HEPUtils::Event& result, double antiktR, double jet_pt_min)
+    void convertParticleEvent(const EventT& pevt, HEPUtils::Event& result, std::vector<jet_collection_settings> all_jet_collection_settings, str jetcollection_taus, double jet_pt_min)
     {
       result.clear();
 
@@ -113,7 +115,11 @@ namespace Gambit
           {
             // Veto bosons not decaying into quarks or gluons
             abschildID = abs(childID);
-            if (abschildID == MCUtils::PID::Z0 || abschildID == MCUtils::PID::WPLUS || abschildID == MCUtils::PID::HIGGS || abschildID == MCUtils::PID::ELECTRON || abschildID == MCUtils::PID::MUON || abschildID == MCUtils::PID::TAU || abschildID == MCUtils::PID::NU_E || abschildID == MCUtils::PID::NU_MU || abschildID == MCUtils::PID::NU_TAU || abschildID == MCUtils::PID::GAMMA)
+            if (abschildID == MCUtils::PID::Z0 || abschildID == MCUtils::PID::WPLUS ||
+                abschildID == MCUtils::PID::HIGGS || abschildID == MCUtils::PID::ELECTRON ||
+                abschildID == MCUtils::PID::MUON || abschildID == MCUtils::PID::TAU ||
+                abschildID == MCUtils::PID::NU_E || abschildID == MCUtils::PID::NU_MU ||
+                abschildID == MCUtils::PID::NU_TAU || abschildID == MCUtils::PID::GAMMA)
             {
               isGoodBoson = false;
             }
@@ -127,7 +133,8 @@ namespace Gambit
           {
             isGoodBoson = false;
           }
-          // Check that the vector bosons do not come from a Higgs boson or top quark (in which case the tagging efficiency would be different)
+          // Check that the vector bosons do not come from a Higgs boson or top quark
+          // (in which case the tagging efficiency would be different)
           int absmotherID = abs(get_unified_mother1_pid(p, pevt));
           if(absmotherID == MCUtils::PID::HIGGS || absmotherID == MCUtils::PID::TQUARK)
           {
@@ -140,7 +147,7 @@ namespace Gambit
             if(apid == MCUtils::PID::WPLUS) WCandidates.push_back(HEPUtils::Particle(p4,pid));
             if(apid == MCUtils::PID::HIGGS) hCandidates.push_back(HEPUtils::Particle(p4,pid));
           }
-        }        
+        }
 
         //We only want final state particles:
         if (!get_unified_isFinal(p)) continue;
@@ -188,91 +195,103 @@ namespace Gambit
       }
 
       /// Jet finding
-      /// @todo Choose jet algorithm via detector _settings? Run several algs?
-      const FJNS::JetDefinition jet_def(FJNS::antikt_algorithm, antiktR);
-      FJNS::ClusterSequence cseq(jetparticles, jet_def);
-      std::vector<FJNS::PseudoJet> pjets = sorted_by_pt(cseq.inclusive_jets(jet_pt_min));
-
-      /// Do jet b-tagging, etc. and add to the Event
-      /// @todo Use ghost tagging?
-      /// @note We need to _remove_ this b-tag in the detector sim if outside the tracker acceptance!
-      for (auto& pj : pjets)
+      for (jet_collection_settings jetcollection : all_jet_collection_settings)
       {
-        HEPUtils::P4 jetMom = HEPUtils::mk_p4(pj);
-        /// @todo Replace with HEPUtils::any(bhadrons, [&](const auto& pb){ pj.delta_R(pb) < 0.4 })
-        bool isB = false;
-        for (HEPUtils::Particle& pb : bpartons)
+        FJNS::JetAlgorithm jet_algorithm = FJalgorithm_map(jetcollection.algorithm);
+        FJNS::Strategy jet_strategy = FJstrategy_map(jetcollection.strategy);
+        FJNS::RecombinationScheme jet_recomscheme = FJRecomScheme_map(jetcollection.recombination_scheme);
+        const FJNS::JetDefinition jet_def(jet_algorithm, jetcollection.R, jet_strategy, jet_recomscheme);
+
+        /// Create and run a new cluster sequence for the given jet collection. 
+        /// The HEPUtils::Event instance ('result') takes ownership of the 
+        /// cluster sequence and a shared_ptr is returned here.
+        std::shared_ptr<const FJNS::ClusterSequence> CSeqBasePtr = result.emplace_clusterseq(jetparticles, jet_def, jetcollection.key);
+        /// Get the resulting pseudojets
+        std::vector<FJNS::PseudoJet> pjets = sorted_by_pt(CSeqBasePtr->inclusive_jets(jet_pt_min));
+
+        /// Do jet b-tagging, etc. and add to the Event
+        /// @todo Use ghost tagging?
+        /// @note We need to _remove_ this b-tag in the detector sim if outside the tracker acceptance!
+        for (auto& pj : pjets)
         {
-          if (jetMom.deltaR_eta(pb.mom()) < 0.4) ///< @todo Hard-coded radius!!!
+          HEPUtils::P4 jetMom = HEPUtils::mk_p4(pj);
+          /// @todo Replace with HEPUtils::any(bhadrons, [&](const auto& pb){ pj.delta_R(pb) < 0.4 })
+          bool isB = false;
+          for (HEPUtils::Particle& pb : bpartons)
           {
-            isB = true;
-            break;
+            if (jetMom.deltaR_eta(pb.mom()) < 0.4) ///< @todo Hard-coded radius!!!
+            {
+              isB = true;
+              break;
+            }
           }
-        }
 
-        bool isC = false;
-        for (HEPUtils::Particle& pc : cpartons)
-        {
-          if (jetMom.deltaR_eta(pc.mom()) < 0.4) ///< @todo Hard-coded radius!!!
+          bool isC = false;
+          for (HEPUtils::Particle& pc : cpartons)
           {
-            isC = true;
-            break;
+            if (jetMom.deltaR_eta(pc.mom()) < 0.4) ///< @todo Hard-coded radius!!!
+            {
+              isC = true;
+              break;
+            }
           }
-        }
 
-        bool isTau = false;
-        int signedTauPID = MCUtils::PID::TAU;
-        for (HEPUtils::Particle& ptau : tauCandidates)
-        {
-          if (jetMom.deltaR_eta(ptau.mom()) < 0.5) ///< @todo Hard-coded radius!!!
+          bool isTau = false;
+          int signedTauPID = MCUtils::PID::TAU;
+          for (HEPUtils::Particle& ptau : tauCandidates)
           {
-            isTau = true;
-            signedTauPID = ptau.pid();
-            break;
+            if (jetMom.deltaR_eta(ptau.mom()) < 0.5) ///< @todo Hard-coded radius!!!
+            {
+              isTau = true;
+              signedTauPID = ptau.pid();
+              break;
+            }
           }
-        }
 
-        bool isW = false;
-        for (HEPUtils::Particle& pW : WCandidates)
-        {
-          if (jetMom.deltaR_eta(pW.mom()) < 1.0) ///< @todo Hard-coded radius from ATLAS-CONF-2021-022, make selectable?
+          bool isW = false;
+          for (HEPUtils::Particle& pW : WCandidates)
           {
-            isW = true;
-            break;
+            if (jetMom.deltaR_eta(pW.mom()) < 1.0) ///< @todo Hard-coded radius from ATLAS-CONF-2021-022, make selectable?
+            {
+              isW = true;
+              break;
+            }
           }
-        }
 
-        bool isZ = false;
-        for (HEPUtils::Particle& pZ : ZCandidates)
-        {
-          if (jetMom.deltaR_eta(pZ.mom()) < 1.0) ///< @todo Hard-coded radius from ATLAS-CONF-2021-022, make selectable?
+          bool isZ = false;
+          for (HEPUtils::Particle& pZ : ZCandidates)
           {
-            isZ = true;
-            break;
+            if (jetMom.deltaR_eta(pZ.mom()) < 1.0) ///< @todo Hard-coded radius from ATLAS-CONF-2021-022, make selectable?
+            {
+              isZ = true;
+              break;
+            }
           }
-        }
 
 
-        bool ish = false;
-        for (HEPUtils::Particle& ph : hCandidates)
-        {
-          if (jetMom.deltaR_eta(ph.mom()) < 1.0) ///< @todo Hard-coded radius from ATLAS-CONF-2021-022, make selectable?
+          bool ish = false;
+          for (HEPUtils::Particle& ph : hCandidates)
           {
-            ish = true;
-            break;
+            if (jetMom.deltaR_eta(ph.mom()) < 1.0) ///< @todo Hard-coded radius from ATLAS-CONF-2021-022, make selectable?
+            {
+              ish = true;
+              break;
+            }
           }
-        }
 
-        // Add to the event (use jet momentum for tau)
-        if (isTau)
-        {
-          HEPUtils::Particle* gp = new HEPUtils::Particle(HEPUtils::mk_p4(pj), signedTauPID);
-          gp->set_prompt();
-          result.add_particle(gp);
-        }
+          // Add to the event (use jet momentum for tau)
+          // Only do this for a single jet collection
+          if (isTau && (jetcollection.key == jetcollection_taus))
+          {
+            HEPUtils::Particle* gp = new HEPUtils::Particle(HEPUtils::mk_p4(pj), signedTauPID);
+            gp->set_prompt();
+            result.add_particle(gp);
+          }
 
-        // Add jet to collection including tags
-        result.add_jet(new HEPUtils::Jet(HEPUtils::mk_p4(pj), isB, isC, isW, isZ, ish));
+          // Add jet to collection including tags and PseudoJet
+          HEPUtils::Jet::TagCounts tags{ {5,int(isB)}, {4,int(isC)}, {23,int(isZ)}, {24,int(isW)}, {25,int(ish)} };
+          result.add_jet(new HEPUtils::Jet(pj, tags), jetcollection.key);
+        }
       }
 
       /// Calculate missing momentum
@@ -299,11 +318,12 @@ namespace Gambit
 
       #ifdef COLLIDERBIT_DEBUG
         // Print event summary
+        cout << "For jet Collection:  " << jetcollection.key << endl;
         cout << "  MET  = " << result.met() << " GeV" << endl;
         cout << "  #e   = " << result.electrons().size() << endl;
         cout << "  #mu  = " << result.muons().size() << endl;
         cout << "  #tau = " << result.taus().size() << endl;
-        cout << "  #jet = " << result.jets().size() << endl;
+        cout << "  #jet = " << result.jets(jetcollection.key).size() << endl;
         cout << "  #pho  = " << result.photons().size() << endl;
         cout << endl;
       #endif
@@ -312,7 +332,7 @@ namespace Gambit
 
     /// Convert a partonic (no hadrons) EventT into an unsmeared HEPUtils::Event
     template<typename EventT>
-    void convertPartonEvent(const EventT& pevt, HEPUtils::Event& result, double antiktR, double jet_pt_min)
+    void convertPartonEvent(const EventT& pevt, HEPUtils::Event& result, std::vector<jet_collection_settings> all_jet_collection_settings, str jetcollection_taus, double jet_pt_min)
     {
       result.clear();
 
@@ -389,40 +409,46 @@ namespace Gambit
       }
 
       /// Jet finding
-      /// @todo choose jet algorithm via _settings?
-      const FJNS::JetDefinition jet_def(FJNS::antikt_algorithm, antiktR);
-      FJNS::ClusterSequence cseq(jetparticles, jet_def);
-      std::vector<FJNS::PseudoJet> pjets = sorted_by_pt(cseq.inclusive_jets(jet_pt_min));
-
-      // Add to the event, with b-tagging info"
-      for (const FJNS::PseudoJet& pj : pjets)
+      for (jet_collection_settings jetcollection : all_jet_collection_settings)
       {
-        // Do jet b-tagging, etc. by looking for b quark constituents (i.e. user index = |parton ID| = 5)
-        /// @note This b-tag is removed in the detector sim if outside the tracker acceptance!
-        const bool isB = HEPUtils::any(pj.constituents(),
-                 [](const FJNS::PseudoJet& c){ return c.user_index() == MCUtils::PID::BQUARK; });
-        const bool isC = HEPUtils::any(pj.constituents(),
-                 [](const FJNS::PseudoJet& c){ return c.user_index() == MCUtils::PID::CQUARK; });
-        result.add_jet(new HEPUtils::Jet(HEPUtils::mk_p4(pj), isB, isC, false, false, false)); // This does not currently deal with boson tagging
-
-        bool isTau=false;
-        int signedTauPID = MCUtils::PID::TAU;
-        for (auto& ptau : tauCandidates)
+        FJNS::JetAlgorithm jet_algorithm = FJalgorithm_map(jetcollection.algorithm);
+        FJNS::Strategy jet_strategy = FJstrategy_map(jetcollection.strategy);
+        FJNS::RecombinationScheme jet_recomscheme = FJRecomScheme_map(jetcollection.recombination_scheme);
+        const FJNS::JetDefinition jet_def(jet_algorithm, jetcollection.R, jet_strategy, jet_recomscheme);
+        std::shared_ptr<const FJNS::ClusterSequence> CSeqBasePtr = result.emplace_clusterseq(jetparticles, jet_def, jetcollection.key);
+        std::vector<FJNS::PseudoJet> pjets = sorted_by_pt(CSeqBasePtr->inclusive_jets(jet_pt_min));
+
+        // Add to the event, with b-tagging info"
+        for (const FJNS::PseudoJet& pj : pjets)
         {
-          HEPUtils::P4 jetMom = HEPUtils::mk_p4(pj);
-          if (jetMom.deltaR_eta(ptau.mom()) < 0.5)
+          // Do jet b-tagging, etc. by looking for b quark constituents (i.e. user index = |parton ID| = 5)
+          /// @note This b-tag is removed in the detector sim if outside the tracker acceptance!
+          const bool isB = HEPUtils::any(pj.constituents(),
+                   [](const FJNS::PseudoJet& c){ return c.user_index() == MCUtils::PID::BQUARK; });
+          const bool isC = HEPUtils::any(pj.constituents(),
+                   [](const FJNS::PseudoJet& c){ return c.user_index() == MCUtils::PID::CQUARK; });
+          result.add_jet(new HEPUtils::Jet(HEPUtils::mk_p4(pj), isB, isC), jetcollection.key); // This does not currently deal with boson tagging
+
+          bool isTau=false;
+          int signedTauPID = MCUtils::PID::TAU;
+          for (auto& ptau : tauCandidates)
           {
-            isTau=true;
-            signedTauPID = ptau.pid();
-            break;
+            HEPUtils::P4 jetMom = HEPUtils::mk_p4(pj);
+            if (jetMom.deltaR_eta(ptau.mom()) < 0.5)
+            {
+              isTau=true;
+              signedTauPID = ptau.pid();
+              break;
+            }
+          }
+          // Add to the event (use jet momentum for tau)
+          // Only do this for a single jet collection
+          if (isTau && (jetcollection.key == jetcollection_taus))
+          {
+            HEPUtils::Particle* gp = new HEPUtils::Particle(HEPUtils::mk_p4(pj), signedTauPID);
+            gp->set_prompt();
+            result.add_particle(gp);
           }
-        }
-        // Add to the event (use jet momentum for tau)
-        if (isTau)
-        {
-          HEPUtils::Particle* gp = new HEPUtils::Particle(HEPUtils::mk_p4(pj), signedTauPID);
-          gp->set_prompt();
-          result.add_particle(gp);
         }
       }
 
diff --git a/ColliderBit/include/gambit/ColliderBit/generateEventPy8Collider.hpp b/ColliderBit/include/gambit/ColliderBit/generateEventPy8Collider.hpp
index a347b1212b..236bd9fd67 100644
--- a/ColliderBit/include/gambit/ColliderBit/generateEventPy8Collider.hpp
+++ b/ColliderBit/include/gambit/ColliderBit/generateEventPy8Collider.hpp
@@ -185,9 +185,9 @@ namespace Gambit
       try
       {
         if (HardScatteringSim.partonOnly)
-          convertPartonEvent(pythia_event, event, HardScatteringSim.antiktR, jet_pt_min);
+          convertPartonEvent(pythia_event, event, HardScatteringSim.all_jet_collection_settings, HardScatteringSim.jetcollection_taus, jet_pt_min);
         else
-          convertParticleEvent(pythia_event, event, HardScatteringSim.antiktR, jet_pt_min);
+          convertParticleEvent(pythia_event, event, HardScatteringSim.all_jet_collection_settings, HardScatteringSim.jetcollection_taus, jet_pt_min);
       }
       // No good.
       catch (Gambit::exception& e)
diff --git a/ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp b/ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp
index 5de42a7dce..04ddc385a6 100644
--- a/ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp
+++ b/ColliderBit/include/gambit/ColliderBit/getPy8Collider.hpp
@@ -98,14 +98,44 @@ namespace Gambit
         // Get options from yaml file.
         const double xsec_veto_default = 0.0;
         const bool partonOnly_default = false;
-        const double antiktR_default = 0.4;
         if (runOptions.hasKey(RunMC.current_collider()))
         {
           YAML::Node colNode = runOptions.getValue<YAML::Node>(RunMC.current_collider());
           Options colOptions(colNode);
           xsec_veto_fb = colOptions.getValueOrDef<double>(xsec_veto_default, "xsec_veto");
           result.partonOnly = colOptions.getValueOrDef<bool>(partonOnly_default, "partonOnly");
-          result.antiktR = colOptions.getValueOrDef<double>(antiktR_default, "antiktR");
+
+          // Fill the jet collection settings
+          if (colOptions.hasKey("jet_collections"))
+          {
+            YAML::Node all_jetcollections_node = colOptions.getValue<YAML::Node>("jet_collections");
+            Options all_jetcollection_options(all_jetcollections_node);
+            std::vector<str> jetcollection_names = all_jetcollection_options.getNames();
+
+            for (str key : jetcollection_names)
+            {
+              YAML::Node current_jc_node = all_jetcollection_options.getValue<YAML::Node>(key);
+              Options current_jc_options(current_jc_node);
+
+              str algorithm = current_jc_options.getValue<str>("algorithm");
+              double R = current_jc_options.getValue<double>("R");
+              str recombination_scheme = current_jc_options.getValue<str>("recombination_scheme");
+              str strategy = current_jc_options.getValue<str>("strategy");
+
+              (result.all_jet_collection_settings).push_back({key, algorithm, R, recombination_scheme, strategy});
+            }
+
+            result.jetcollection_taus = colOptions.getValue<str>("jet_collection_taus");
+            // Throw an error if the "jet_collection_taus" setting does not match an entry in "jet_collections".
+            if (std::find(jetcollection_names.begin(), jetcollection_names.end(), result.jetcollection_taus) == jetcollection_names.end())
+            {
+              ColliderBit_error().raise(LOCAL_INFO,"Please provide the jet_collection_taus setting for jet collections.");
+            }
+          }
+          else
+          {
+            ColliderBit_error().raise(LOCAL_INFO,"Could not find jet_collections option for collider " + RunMC.current_collider() + ". Please provide this in the YAML file.");
+          }
           if (colOptions.hasKey("pythia_settings"))
           {
             std::vector<str> addPythiaOptions = colNode["pythia_settings"].as<std::vector<str> >();
@@ -114,9 +144,7 @@ namespace Gambit
         }
         else
         {
-          xsec_veto_fb = xsec_veto_default;
-          result.partonOnly = partonOnly_default;
-          result.antiktR = antiktR_default;
+          ColliderBit_error().raise(LOCAL_INFO,"Could not find runOptions for collider " + RunMC.current_collider() + ".");
         }
 
         // We need showProcesses for the xsec veto.
diff --git a/ColliderBit/include/gambit/ColliderBit/lhef2heputils.hpp b/ColliderBit/include/gambit/ColliderBit/lhef2heputils.hpp
index 82d2dd98b0..934241ec5b 100644
--- a/ColliderBit/include/gambit/ColliderBit/lhef2heputils.hpp
+++ b/ColliderBit/include/gambit/ColliderBit/lhef2heputils.hpp
@@ -24,6 +24,7 @@
 ///  *********************************************
 
 #include "gambit/cmake/cmake_variables.hpp"
+#include "gambit/ColliderBit/Utils.hpp"
 
 #ifndef EXCLUDE_HEPMC
 
@@ -32,7 +33,17 @@
 /// Forward declaration to cut down on includes
 namespace LHEF { class Reader; }
 
-/// Extract an LHE event as a HEPUtils::Event
-void get_HEPUtils_event(const LHEF::Reader&, HEPUtils::Event&, double);
+namespace Gambit
+{
 
-#endif
\ No newline at end of file
+  namespace ColliderBit
+  {
+
+    /// Extract an LHE event as a HEPUtils::Event
+    void get_HEPUtils_event(const LHEF::Reader&, HEPUtils::Event&, double, std::vector<Gambit::ColliderBit::jet_collection_settings>);
+
+  }
+
+}
+
+#endif
diff --git a/ColliderBit/src/ColliderBit_measurements.cpp b/ColliderBit/src/ColliderBit_measurements.cpp
index 9e632c5101..1964fc12f2 100644
--- a/ColliderBit/src/ColliderBit_measurements.cpp
+++ b/ColliderBit/src/ColliderBit_measurements.cpp
@@ -484,7 +484,7 @@ namespace Gambit
           summary_line << "LHC Contur LogLikes per pool: ";
           result = (*Dep::LHC_measurements).pool_LLR;
 
-          for( auto const& entry : result)
+          for (auto const& entry : result)
           {
             summary_line << entry.first << ":" << entry.second << ", ";
           }
@@ -506,7 +506,7 @@ namespace Gambit
               result[pool_LLR_entry.first + "_" + contur_output_instance.first] = pool_LLR_entry.second;
             }
           }
-          for( auto const& entry : result)
+          for (auto const& entry : result)
           {
             summary_line << entry.first << ":" << entry.second << ", ";
           }
@@ -522,7 +522,7 @@ namespace Gambit
           summary_line << "LHC Contur LogLikes per pool: ";
           result = (*Dep::LHC_measurements).pool_tags;
 
-          for( auto const& entry : result)
+          for (auto const& entry : result)
           {
             summary_line << entry.first << ":" << entry.second << ", ";
           }
@@ -547,7 +547,7 @@ namespace Gambit
               result[pool_LLR_entry.first + "_" + contur_output_instance.first] = pool_LLR_entry.second;
             }
           }
-          for( auto const& entry : result)
+          for (auto const& entry : result)
           {
             summary_line << entry.first << ":" << entry.second << ", ";
           }
diff --git a/ColliderBit/src/Utils.cpp b/ColliderBit/src/Utils.cpp
index a4e93a2946..f8250bc34e 100644
--- a/ColliderBit/src/Utils.cpp
+++ b/ColliderBit/src/Utils.cpp
@@ -31,6 +31,7 @@
 ///  *********************************************
 
 #include "gambit/ColliderBit/Utils.hpp"
+#include "gambit/ColliderBit/ColliderBit_eventloop.hpp"
 #include "gambit/Utils/threadsafe_rng.hpp"
 #include <iostream>
 using namespace std;
@@ -41,6 +42,47 @@ namespace Gambit
   {
 
 
+    /// Storage of different FastJet methods
+    FJNS::JetAlgorithm FJalgorithm_map(str algorithm)
+    {
+      FJNS::JetAlgorithm result;
+      if (algorithm == "antikt") {result = FJNS::antikt_algorithm;}
+      else if (algorithm == "cambridge") {result = FJNS::cambridge_algorithm;}
+      else if (algorithm == "kt") {result = FJNS::kt_algorithm;}
+      else if (algorithm == "genkt") {result = FJNS::genkt_algorithm;}
+      else if (algorithm == "cambridge_for_passive") {result = FJNS::cambridge_for_passive_algorithm;}
+      else
+      {
+        ColliderBit_error().raise(LOCAL_INFO, "Could not find jet algorithm in list available. Please add the missing option to the FJalgorithm_map function in ColliderBit/src/Utils.cpp.");
+      }
+      return result;
+    }
+
+    FJNS::Strategy FJstrategy_map(str strategy)
+    {
+      FJNS::Strategy result;
+      if (strategy == "Best") {result = FJNS::Best;}
+      else if (strategy == "NlnN") {result = FJNS::NlnN;}
+      else
+      {
+        ColliderBit_error().raise(LOCAL_INFO, "Could not find jet strategy in list available. Please add the missing option to the FJstrategy_map function in ColliderBit/src/Utils.cpp.");
+      }
+      return result;
+    }
+
+    FJNS::RecombinationScheme FJRecomScheme_map(str reco_scheme)
+    {
+      FJNS::RecombinationScheme result;
+      if (reco_scheme == "E_scheme") {result = FJNS::E_scheme;}
+      else if (reco_scheme == "pt_scheme") {result = FJNS::pt_scheme;}
+      else if (reco_scheme == "pt2_scheme") {result = FJNS::pt2_scheme;}
+      else
+      {
+        ColliderBit_error().raise(LOCAL_INFO, "Could not find jet recombination scheme in list available. Please add the missing option to the FJRecomScheme_map function in ColliderBit/src/Utils.cpp.");
+      }
+      return result;
+    }
+
     bool random_bool(double eff)
     {
       /// @todo Handle out-of-range eff values
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEPStop_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEPStop_36invfb.cpp
index d0b6028697..a870ad6e9c 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEPStop_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEPStop_36invfb.cpp
@@ -279,7 +279,7 @@ namespace Gambit {
         const std::vector<double>  b = {0,10000.};
         const std::vector<double> c = {0.77}; // set b-tag efficiency to 77%
         HEPUtils::BinnedFn2D<double> _eff2d(a,b,c);
-        for (const HEPUtils::Jet* jet : event->jets())
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
         {
           bool hasTag=has_tag(_eff2d, fabs(jet->eta()), jet->pT());
           if (jet->pT() > 20. && fabs(jet->eta()) < 2.8)
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_139invfb.cpp
index 2bb39009e9..eb91ebffd0 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_139invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_139invfb.cpp
@@ -101,7 +101,7 @@ namespace Gambit
           // Get baseline jets
           /// @todo Drop b-tag if pT < 50 GeV or |eta| > 2.5?
           vector<const Jet*> baselineJets;
-          for (const Jet* jet : event->jets())
+          for (const Jet* jet : event->jets("antikt_R04"))
           {
             if (jet->pT() > 20. && jet->abseta() < 2.8)
             {
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_13invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_13invfb.cpp
index 0b4206c9c0..61896dd221 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_13invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_13invfb.cpp
@@ -88,7 +88,7 @@ namespace Gambit
           // Get baseline jets
           /// @todo Drop b-tag if pT < 50 GeV or |eta| > 2.5?
           vector<const Jet*> baselineJets;
-          for (const Jet* jet : event->jets())
+          for (const Jet* jet : event->jets("antikt_R04"))
             if (jet->pT() > 20. && jet->abseta() < 2.8)
             {
               baselineJets.push_back(jet);
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_36invfb.cpp
index 406a083714..dfee9f7788 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_0LEP_36invfb.cpp
@@ -115,7 +115,7 @@ namespace Gambit {
         // Get baseline jets
         /// @todo Drop b-tag if pT < 50 GeV or |eta| > 2.5?
         vector<const Jet*> baselineJets;
-        for (const Jet* jet : event->jets())
+        for (const Jet* jet : event->jets("antikt_R04"))
           if (jet->pT() > 20. && jet->abseta() < 2.8) {
             baselineJets.push_back(jet);
           }
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_1LEPStop_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_1LEPStop_36invfb.cpp
index 20b46b48e1..d40bda7b5c 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_1LEPStop_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_1LEPStop_36invfb.cpp
@@ -437,7 +437,7 @@ namespace Gambit {
         const std::vector<double>  b = {0,10000.};
         const std::vector<double> c = {0.77}; // set b-tag efficiency to 77%
         HEPUtils::BinnedFn2D<double> _eff2d(a,b,c);
-        for (const HEPUtils::Jet* jet : event->jets())
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
         {
           bool hasTag=has_tag(_eff2d, fabs(jet->eta()), jet->pT());
           if (jet->pT() > 20. && fabs(jet->eta()) < 4.9)
@@ -1423,4 +1423,4 @@ namespace Gambit {
   }
 }
 
-#endif
\ No newline at end of file
+#endif
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_1Lep2b_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_1Lep2b_139invfb.cpp
index 468b221f80..70227f866b 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_1Lep2b_139invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_1Lep2b_139invfb.cpp
@@ -180,7 +180,7 @@ namespace Gambit
 
 
         vector<const Jet*> baselineJets;
-        for (const HEPUtils::Jet* jet : event->jets())
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
         {
           if (jet->pT() > 20.0 && jet->abseta() < 4.5)
           {
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2BoostedBosons_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2BoostedBosons_139invfb.cpp
index 90cb601b05..9e490760cc 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2BoostedBosons_139invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2BoostedBosons_139invfb.cpp
@@ -54,6 +54,7 @@ namespace Gambit
     {
 
       protected:
+
         // Signal region map
         std::map<string, EventCounter> _counters = {
           // Exclusion regions
@@ -166,9 +167,9 @@ namespace Gambit
 
           // Look at jets to see if they fulfil criteria for fat jets
           vector<const HEPUtils::Jet*> fatJets;
-          for (const HEPUtils::Jet* jet : event->jets())
+          for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
           {
-            //  cout  << jet->pT() << " " << jet->mass() << " Z-tag " <<  jet->Ztag() << " W-tag " << jet->Wtag() << " " << endl;
+            //  cout  << jet->pT() << " " << jet->mass() << " Z-tag " <<  jet->tagged(23) << " W-tag " << jet->tagged(24) << " " << endl;
             if (jet->pT() > 200. && fabs(jet->eta()) < 2.0 && jet->mass() > 40.)
             {
               fatJets.push_back(jet);
@@ -194,13 +195,13 @@ namespace Gambit
           for (const HEPUtils::Jet* jet : fatJets)
           {
             // Tag W
-            if( jet->Wtag() && random_bool( _eff1dW.get_at( jet->pT() ) ) ) nW++;
+            if( jet->tagged(24) && random_bool( _eff1dW.get_at( jet->pT() ) ) ) nW++;
             // Tag Z
-            if( jet->Ztag() && random_bool( _eff1dZ.get_at( jet->pT() ) ) ) nZ++;
+            if( jet->tagged(23) && random_bool( _eff1dZ.get_at( jet->pT() ) ) ) nZ++;
             // Tag SM Higgs
-            if( jet->htag() && random_bool( _eff1dH.get_at( jet->pT() )  ) ) nH++;
+            if( jet->tagged(25) && random_bool( _eff1dH.get_at( jet->pT() )  ) ) nH++;
             // Misstag as Z or W
-            if( !jet->Wtag() && !jet->Ztag() )
+            if( !jet->tagged(23) && !jet->tagged(24) )
             {
               if( random_bool( _eff1dZmiss.get_at( jet->pT() ) )  ) nZ++;
               if( random_bool( _eff1dWmiss.get_at( jet->pT() ) )  ) nW++;
@@ -212,17 +213,17 @@ namespace Gambit
           //if(nH > 0 ) cout << "nZ " << nZ << " nW " << nW << " nV " << nV << " nH " << nH << endl;
 
           // b-jet tagging
-          /* There is a difference here wrt the actual analysis where small
-           sliding radius track jets are used, and the number of such b-jets are
-           counted. This means that the rejection for b-jets in the 4Q SRs has to
-           be changed. We use the conservative choice of rejecting all events with
-           a b-labeled large radius jet and mis-tagging large radius non-b-jets
-           according to the mis-tag probabilities of the small radius track jets.
-          */
-  //        double btag = 0.83;
+          // There is a difference here wrt the actual analysis where small
+          // variable-radius track jets are used, and the number of such b-jets are
+          // counted. This means that the rejection for b-jets in the 4Q SRs has to
+          // be changed. We use the conservative choice of rejecting all events with
+          // a b-labeled large radius jet and mis-tagging large radius non-b-jets
+          // according to the mis-tag probabilities of the small radius track jets.
+
+          // double btag = 0.83;
           double cmisstag = 1/3.; double misstag = 1./33.;
           int nb = 0;
-          for ( const HEPUtils::Jet* jet : event->jets() )
+          for ( const HEPUtils::Jet* jet : event->jets("antikt_R04") )
           {
             // Tag b-jet
             if( jet->btag() ) nb++;
@@ -370,7 +371,7 @@ namespace Gambit
 
           // Compare meff spectrum
           cout << "Meff SR-4Q-VV\t" << "GAMBIT\t" << "ATLAS " << endl;
-          for( size_t j = 0; j < _meff_4QVV.size(); j++){
+          for (size_t j = 0; j < _meff_4QVV.size(); j++){
             cout << "[" << _meff_bins[j] << ", " << _meff_bins[j+1] << "]\t" << _meff_4QVV[j]*_scale << "\t" << _meff_4QVV_model[j] << endl;
           }
 
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2LEPStop_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2LEPStop_139invfb.cpp
index c2b0f1f15d..201f9b0f24 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2LEPStop_139invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2LEPStop_139invfb.cpp
@@ -163,7 +163,7 @@ namespace Gambit
             //   primary vertex (see paper)
 
             vector<const HEPUtils::Jet*> baselineJets;
-            for (const HEPUtils::Jet* jet : event->jets()) 
+            for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) 
             {
                 if (jet->pT() > 20. && jet->abseta() < 2.8) 
                 {
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2LEPStop_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2LEPStop_36invfb.cpp
index f309b7dacf..974df46699 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2LEPStop_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2LEPStop_36invfb.cpp
@@ -175,7 +175,7 @@ namespace Gambit {
             // Jets
             vector<const HEPUtils::Jet*> blJets;          // Used for SR-2body and SR-3body
             vector<const HEPUtils::Jet*> baselineJets;    // Used for SR-4body
-            for (const HEPUtils::Jet* jet : event->jets()) {
+            for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
                 if (jet->pT() > 20. && fabs(jet->eta()) < 2.8) blJets.push_back(jet);
                 if (jet->pT() > 20. && fabs(jet->eta()) < 2.8) baselineJets.push_back(jet);
             }
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_Z_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_Z_139invfb.cpp
index a11ff77657..f275327166 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_Z_139invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_Z_139invfb.cpp
@@ -148,7 +148,7 @@ namespace Gambit
         // Jets with pT < 120 GeV and |η| < 2.8 have an efficiency of 90%
         // Mising:  cut based on detector noise and non-collision backgrounds
         double jet_eff = 0.9;
-        for (const HEPUtils::Jet* jet : event->jets())
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
         {
           if (jet->pT()>20. && jet->abseta()<2.8)
             if( (jet->pT() >= 120. || jet->abseta() >= 2.5) || random_bool(jet_eff) ) baselineJets.push_back(jet);
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_chargino_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_chargino_139invfb.cpp
index bee304bfe5..af0c953e8b 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_chargino_139invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_chargino_139invfb.cpp
@@ -213,7 +213,7 @@ namespace Gambit
 
         // Jets
         vector<const HEPUtils::Jet*> candJets;
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT() > 20. && fabs(jet->eta()) < 2.4)
             candJets.push_back(jet);
         }
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_chargino_80invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_chargino_80invfb.cpp
index 6426698cf9..f7cb98872d 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_chargino_80invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2OSLEP_chargino_80invfb.cpp
@@ -199,7 +199,7 @@ namespace Gambit
 
         // Jets
         vector<const HEPUtils::Jet*> candJets;
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT() > 20. && fabs(jet->eta()) < 2.5)
             candJets.push_back(jet);
         }
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2bMET_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2bMET_36invfb.cpp
index 93d5939deb..63e54cfcb7 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2bMET_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_2bMET_36invfb.cpp
@@ -357,7 +357,7 @@ namespace Gambit
           const std::vector<double>  b = {0,10000.};
           const std::vector<double> c = {0.77}; // set b-tag efficiency to 77%
           HEPUtils::BinnedFn2D<double> _eff2d(a,b,c);
-          for (const HEPUtils::Jet* jet : event->jets())
+          for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
           {
             bool hasTag=has_tag(_eff2d, fabs(jet->eta()), jet->pT());
             if (jet->pT() > 20. && fabs(jet->eta()) < 4.8)
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3LEP_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3LEP_139invfb.cpp
index e68af59ba9..d5368bc4d4 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3LEP_139invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3LEP_139invfb.cpp
@@ -219,7 +219,7 @@ namespace Gambit
 
           // Baseline jets
           vector<const HEPUtils::Jet*> baselineJets;
-          for (const HEPUtils::Jet* jet : event->jets())
+          for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
           {
             if (jet->pT() > 20. && fabs(jet->eta()) < 4.5 )
               baselineJets.push_back(jet);
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3b_24invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3b_24invfb.cpp
index 957738b8d2..4fc7e6def0 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3b_24invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3b_24invfb.cpp
@@ -197,7 +197,7 @@ namespace Gambit {
         ATLAS::applyMuonEff(muons);
 
         vector<const HEPUtils::Jet*> candJets;
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT() > 25. && fabs(jet->eta()) < 2.5)
             candJets.push_back(jet);
         }
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3b_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3b_36invfb.cpp
index d82199d666..43b157e0a0 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3b_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_3b_36invfb.cpp
@@ -162,7 +162,7 @@ namespace Gambit {
         ATLAS::applyMuonEff(muons);
 
         vector<const HEPUtils::Jet*> candJets;
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT() > 20. && fabs(jet->eta()) < 2.8)
             candJets.push_back(jet);
         }
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_4LEP_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_4LEP_139invfb.cpp
index 54e610942a..9950a14d55 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_4LEP_139invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_4LEP_139invfb.cpp
@@ -326,7 +326,7 @@ namespace Gambit
         // Since tau efficiencies are not applied as part of the BuckFast ATLAS sim we apply it here
         ATLAS::applyTauEfficiencyR2(baselineTaus);
 
-        for (const HEPUtils::Jet* jet : event->jets())
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
         {
           if (jet->pT()>20. && jet->abseta()<2.8) baselineJets.push_back(jet);
         }
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_4LEP_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_4LEP_36invfb.cpp
index a40e20e089..4150d66976 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_4LEP_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_4LEP_36invfb.cpp
@@ -224,7 +224,7 @@ namespace Gambit
         // Since tau efficiencies are not applied as part of the BuckFast ATLAS sim we apply it here
         ATLAS::applyTauEfficiencyR2(baselineTaus);
 
-        for (const HEPUtils::Jet* jet : event->jets())
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
         {
           if (jet->pT()>20. && jet->abseta()<2.8) baselineJets.push_back(jet);
         }
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MONOJET_139infb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MONOJET_139infb.cpp
index 0a1750e5c9..46e599ea3c 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MONOJET_139infb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MONOJET_139infb.cpp
@@ -141,7 +141,7 @@ namespace Gambit
 
           // Get jets (0.9 is to emulate the requirement of coming from a primary vertex)
           vector<const Jet*> baselineJets;
-          for (const Jet* jet : event->jets())
+          for (const Jet* jet : event->jets("antikt_R04"))
           {
             if ((jet->pT() > 30) && (jet->abseta() < 2.8))
             {
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_36invfb.cpp
index 1e6d422dd1..b4acd611f9 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_36invfb.cpp
@@ -187,7 +187,7 @@ namespace Gambit {
         ATLAS::applyMuonEff(baselineMuons);
 
         vector<const HEPUtils::Jet*> baselineJets;
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT()>20. && jet->abseta()<4.5)baselineJets.push_back(jet);
         }
 
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_confnote_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_confnote_36invfb.cpp
index 261b9b45fa..5326185a2e 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_confnote_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_confnote_36invfb.cpp
@@ -176,7 +176,7 @@ namespace Gambit {
         ATLAS::applyMuonEff(baselineMuons);
 
         vector<const HEPUtils::Jet*> baselineJets;
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT()>20. && jet->abseta()<4.5)baselineJets.push_back(jet);
         }
 
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_strong_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_strong_139invfb.cpp
index 7644e2fdbd..84c9489362 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_strong_139invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_MultiLEP_strong_139invfb.cpp
@@ -179,7 +179,7 @@ namespace Gambit {
 
         // Get baseline jets
         /// @todo Drop b-tag if |eta| > 2.5?
-        for (const Jet* jet : event->jets()) 
+        for (const Jet* jet : event->jets("antikt_R04")) 
         {
           if (jet->pT() > 20. && jet->abseta() < 2.8) 
           {
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_PhotonGGM_1Photon_139invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_PhotonGGM_1Photon_139invfb.cpp
index e75cb51172..648f0d7691 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_PhotonGGM_1Photon_139invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_PhotonGGM_1Photon_139invfb.cpp
@@ -120,7 +120,7 @@ namespace Gambit
         // - pT > 30
         // - |eta| < 2.5
         vector<const HEPUtils::Jet*> signalJets;
-        for (const HEPUtils::Jet* jet : event->jets())
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
         {
           if (jet->pT() > 30. && fabs(jet->eta()) < 2.5)
           {
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_PhotonGGM_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_PhotonGGM_36invfb.cpp
index 5f54059dc5..47de7698a0 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_PhotonGGM_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_PhotonGGM_36invfb.cpp
@@ -183,7 +183,7 @@ namespace Gambit
           // Jets
           vector<const HEPUtils::Jet*> jets28;
           vector<const HEPUtils::Jet*> jets28_nophooverlap;
-          for (const HEPUtils::Jet* jet : event->jets())
+          for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
           {
             if (jet->pT() > 30. && fabs(jet->eta()) < 2.8)
             {
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb.cpp
index 14b3c7f742..5dd2c42d55 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_RJ3L_lowmass_36invfb.cpp
@@ -645,7 +645,7 @@ namespace Gambit {
         const std::vector<double>  b = {0,10000.};
         const std::vector<double> c = {0.77}; // set b-tag efficiency to 77%
         HEPUtils::BinnedFn2D<double> _eff2d(a,b,c);
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           bool hasTag=has_tag(_eff2d, jet->abseta(), jet->pT());
           if (jet->pT() > 20. && fabs(jet->eta()) < 2.4) {
             if(jet->btag() && hasTag){
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_ZGammaGrav_CONFNOTE_80invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_ZGammaGrav_CONFNOTE_80invfb.cpp
index d5a0dce055..f521963413 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_ZGammaGrav_CONFNOTE_80invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_13TeV_ZGammaGrav_CONFNOTE_80invfb.cpp
@@ -67,7 +67,7 @@ namespace Gambit {
 
         // Jets
         JetPtrs jets;
-        for (const Jet* j : event->jets())
+        for (const Jet* j : event->jets("antikt_R04"))
           if (j->pT() > 20. && j->absrap() < 4.4)
             jets.push_back(j);
 
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_7TeV_1OR2LEPStop_4_7invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_7TeV_1OR2LEPStop_4_7invfb.cpp
index 0d8dbe10f6..ec043c17ee 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_7TeV_1OR2LEPStop_4_7invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_7TeV_1OR2LEPStop_4_7invfb.cpp
@@ -211,7 +211,7 @@ Analysis_ATLAS_7TeV_1OR2LEPStop_4_7invfb()
 	incrementCut(Total_events);
 	std::vector<const Particle*> electrons = event->electrons();
 	std::vector<const Particle*> muons = event->muons();
-	std::vector<const Jet*> jets = event->jets();
+	std::vector<const Jet*> jets = event->jets("antikt_R04");
 
 	electrons = AnalysisUtil::filterPtEta(electrons, 20, 2.47);
 	muons = AnalysisUtil::filterPtEta(muons, 10, 2.4);
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_7TeV_2LEPStop_4_7invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_7TeV_2LEPStop_4_7invfb.cpp
index f48ed86d47..fc63bbbeb9 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_7TeV_2LEPStop_4_7invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_7TeV_2LEPStop_4_7invfb.cpp
@@ -49,7 +49,7 @@ namespace Gambit
         // get the jets and leptons filtered by their pt and eta requirements
         electrons = AnalysisUtil::filterPtEta(electrons, 10, 2.47);
         muons = AnalysisUtil::filterPtEta(muons, 10, 2.4);
-        std::vector<const Jet*> jets = AnalysisUtil::filterPtEta(event->jets(), 20, 4.5);
+        std::vector<const Jet*> jets = AnalysisUtil::filterPtEta(event->jets("antikt_R04"), 20, 4.5);
 
         // check if any of the triggers were triggered
         bool eeTrigger = AnalysisUtil::isMultipleParticleTriggered(electrons, {17, 17});
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_0LEPStop_20invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_0LEPStop_20invfb.cpp
index d8de98095c..83bb54ca65 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_0LEPStop_20invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_0LEPStop_20invfb.cpp
@@ -157,7 +157,7 @@ namespace Gambit {
         const std::vector<double>  b = {0,10000.};
         const std::vector<double> c = {0.7};
         HEPUtils::BinnedFn2D<double> _eff2d(a,b,c);
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT() > 20. && fabs(jet->eta()) < 4.5) baselineJets.push_back(jet);
           bool hasTag=has_tag(_eff2d, fabs(jet->eta()), jet->pT());
           if (jet->pT() > 20. && fabs(jet->eta()) < 4.5) {
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_0LEP_20invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_0LEP_20invfb.cpp
index 0192565cd6..2109e92358 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_0LEP_20invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_0LEP_20invfb.cpp
@@ -77,7 +77,7 @@ namespace Gambit {
         ATLAS::applyMuonEff(baselineMuons);
 
         vector<const HEPUtils::Jet*> baselineJets;
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT() > 20. && fabs(jet->eta()) < 4.5) baselineJets.push_back(jet);
         }
 
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_1LEPStop_20invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_1LEPStop_20invfb.cpp
index 7ceba10beb..b74a546d35 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_1LEPStop_20invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_1LEPStop_20invfb.cpp
@@ -206,6 +206,8 @@ namespace Gambit {
 
       void run(const HEPUtils::Event* event) {
 
+        static const std::string jet_collection_name = "antikt_R04";
+
         // Missing energy
         HEPUtils::P4 ptot = event->missingmom();
         double met = event->met();
@@ -214,7 +216,7 @@ namespace Gambit {
         vector<const HEPUtils::Particle*> baselineElectrons;
         for (const HEPUtils::Particle* electron : event->electrons()) {
           if (electron->pT() > 10. && electron->abseta() < 2.47 &&
-              !object_in_cone(*event, *electron, 0.1*electron->pT(), 0.2)) baselineElectrons.push_back(electron);
+              !object_in_cone(*event, jet_collection_name, *electron, 0.1*electron->pT(), 0.2)) baselineElectrons.push_back(electron);
         }
 
         // Apply electron efficiency
@@ -224,7 +226,7 @@ namespace Gambit {
         vector<const HEPUtils::Particle*> baselineMuons;
         for (const HEPUtils::Particle* muon : event->muons()) {
           if (muon->pT() > 10. && muon->abseta() < 2.4 &&
-              !object_in_cone(*event, *muon, 1.8, 0.2)) baselineMuons.push_back(muon);
+              !object_in_cone(*event, jet_collection_name, *muon, 1.8, 0.2)) baselineMuons.push_back(muon);
         }
 
         // Apply muon efficiency
@@ -232,7 +234,7 @@ namespace Gambit {
 
         // Get b jets with efficiency and mistag (fake) rates
         vector<const HEPUtils::Jet*> baselineJets, bJets; // trueBJets; //for debugging
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets(jet_collection_name)) {
           if (jet->pT() > 20. && jet->abseta() < 10.0) baselineJets.push_back(jet);
           if (jet->abseta() < 2.5 && jet->pT() > 25.) {
             if ((jet->btag() && Random::draw() < 0.75) || (!jet->btag() && Random::draw() < 0.02)) bJets.push_back(jet);
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_1LEPbb_20invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_1LEPbb_20invfb.cpp
index 0c593b452a..dcb7d864b7 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_1LEPbb_20invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_1LEPbb_20invfb.cpp
@@ -109,7 +109,7 @@ namespace Gambit {
         ATLAS::applyMuonEff(baselineMuons);
 
         vector<const HEPUtils::Jet*> baselineJets;
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT()>20. && fabs(jet->eta())<4.5) baselineJets.push_back(jet);
         }
 
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2LEPEW_20invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2LEPEW_20invfb.cpp
index dbe5c070b4..cb9c2e0ad2 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2LEPEW_20invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2LEPEW_20invfb.cpp
@@ -287,7 +287,7 @@ namespace Gambit
           ATLAS::applyMuonEff(signalMuons);
 
           vector<const HEPUtils::Jet*> signalJets;
-          for (const HEPUtils::Jet* jet : event->jets())
+          for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
           {
             if (jet->pT() > 20. && fabs(jet->eta()) < 4.5) signalJets.push_back(jet);
             //if(jet->btag() && fabs(jet->eta()) < 2.5 && jet->pT() > 20.) bJets.push_back(jet);
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2LEPStop_20invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2LEPStop_20invfb.cpp
index 2c681d4b2d..886bde501b 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2LEPStop_20invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2LEPStop_20invfb.cpp
@@ -100,7 +100,7 @@ namespace Gambit {
         vector<const HEPUtils::Jet*> baselineJets;
         vector<const HEPUtils::Jet*> bJets;
         vector<const HEPUtils::Jet*> trueBJets; //for debugging
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT() > 20. && fabs(jet->eta()) < 2.5) baselineJets.push_back(jet);
           if(jet->btag() && fabs(jet->eta()) < 2.5 && jet->pT() > 20.) bJets.push_back(jet);
         }
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2bStop_20invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2bStop_20invfb.cpp
index 18b5eaba81..367131cf37 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2bStop_20invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_2bStop_20invfb.cpp
@@ -123,7 +123,7 @@ namespace Gambit
           vector<const HEPUtils::Jet*> bJets;
           vector<const HEPUtils::Jet*> trueBJets; //for debugging
 
-          for (const HEPUtils::Jet* jet : event->jets())
+          for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
           {
             if (jet->pT() > 20. && fabs(jet->eta()) < 4.9) baselineJets.push_back(jet);
           }
diff --git a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_3LEPEW_20invfb.cpp b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_3LEPEW_20invfb.cpp
index 12f006036e..58192f8470 100644
--- a/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_3LEPEW_20invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_ATLAS_8TeV_3LEPEW_20invfb.cpp
@@ -215,7 +215,7 @@ namespace Gambit {
         vector<const HEPUtils::Jet*> signalJets;
         vector<const HEPUtils::Jet*> bJets;
 
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT() > 20. && fabs(jet->eta()) < 2.5) signalJets.push_back(jet);
           //if(jet->btag() && fabs(jet->eta()) < 2.5 && jet->pT() > 20.) bJets.push_back(jet);
         }
diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_137invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_137invfb.cpp
index 265258625e..9b150cd6d9 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_137invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_137invfb.cpp
@@ -65,7 +65,7 @@ namespace Gambit
 
           // Get jets
           vector<const Jet*> jets24, jets50;
-          for (const Jet* jet : event->jets())
+          for (const Jet* jet : event->jets("antikt_R04"))
           {
             if (jet->pT() < 30) continue;
             if (jet->abseta() < 2.4) jets24.push_back(jet);
diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_13invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_13invfb.cpp
index 985c2a5bac..6f0d7ca8ef 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_13invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_13invfb.cpp
@@ -57,7 +57,7 @@ namespace Gambit {
 
         // Get baseline jets
         vector<const Jet*> jets24, jets50;
-        for (const Jet* jet : event->jets()) {
+        for (const Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT() < 30) continue;
           if (jet->abseta() < 2.4) jets24.push_back(jet);
           if (jet->abseta() < 5.0) jets50.push_back(jet);
diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_36invfb.cpp
index 2f79a8fd7a..7d0b853edf 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_0LEP_36invfb.cpp
@@ -57,7 +57,7 @@ namespace Gambit {
 
         // Get baseline jets
         vector<const Jet*> jets24, jets50;
-        for (const Jet* jet : event->jets()) {
+        for (const Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT() < 30) continue;
           if (jet->abseta() < 2.4) jets24.push_back(jet);
           if (jet->abseta() < 5.0) jets50.push_back(jet);
diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_1LEPStop_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_1LEPStop_36invfb.cpp
index 265e290fe5..b097fe873b 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_1LEPStop_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_1LEPStop_36invfb.cpp
@@ -154,7 +154,7 @@ namespace Gambit {
             // Jets
             vector<const HEPUtils::Jet*> baselineJets;
             vector<const HEPUtils::Jet*> fullJets;
-            for (const HEPUtils::Jet* jet : event->jets()) {
+            for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
                 if (jet->pT() > 30. && jet->abseta() < 2.4) baselineJets.push_back(jet);
                 if (jet->abseta() < 5.0) fullJets.push_back(jet);
             }
@@ -294,7 +294,7 @@ namespace Gambit {
 
             // Mlb
             double deltaRlb=9999.;
-            double Mlb;
+            double Mlb = -1; ///< what if no real value is found???
             for (const HEPUtils::Jet* bj :mediumbJets) {
                 if (deltaRlb > bj->mom().deltaR_eta(Leptons.at(0)->mom())){
                     deltaRlb = bj->mom().deltaR_eta(Leptons.at(0)->mom());
@@ -395,7 +395,7 @@ namespace Gambit {
             if (baselineJets.size()>=4 and tmod>10 and Mlb<=175 and met>=450) _counters.at("aggregateSR2").add_event(event);
             if (baselineJets.size()>=4 and tmod<=0 and Mlb> 175 and met>=450) _counters.at("aggregateSR3").add_event(event);
             if (baselineJets.size()>=4 and tmod> 0 and Mlb> 175 and met>=450) _counters.at("aggregateSR4").add_event(event);
-            if(baselineJets.size()>=5 and leadjet_nob and deltaPhi_j12 >0.5 and Leptons.at(0)->pT() < 150 and Leptons.at(0)->mom().deltaPhi(ptot)<2. ){
+            if (baselineJets.size()>=5 and leadjet_nob and deltaPhi_j12 >0.5 and Leptons.at(0)->pT() < 150 and Leptons.at(0)->mom().deltaPhi(ptot)<2. ){
                 if( met>=450 ) _counters.at("aggregateSR5").add_event(event);
             }
         return;
diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_1LEPbb_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_1LEPbb_36invfb.cpp
index 462e5af12a..4f1e309e23 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_1LEPbb_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_1LEPbb_36invfb.cpp
@@ -119,7 +119,7 @@ namespace Gambit
         }
 
         vector<const HEPUtils::Jet*> baselineJets;
-        for (const HEPUtils::Jet* jet : event->jets())
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
         {
           if (jet->pT()>25. &&fabs(jet->eta())<2.4)baselineJets.push_back(jet);
         }
diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_1Photon1Lepton_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_1Photon1Lepton_36invfb.cpp
index 614fd45578..922e185e62 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_1Photon1Lepton_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_1Photon1Lepton_36invfb.cpp
@@ -270,7 +270,7 @@ namespace Gambit {
 
         // Jets
         vector<const HEPUtils::Jet*> jets;
-        for (const HEPUtils::Jet* jet : event->jets())
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
         {
           if (jet->pT()>30. && jet->abseta()<2.5) jets.push_back(jet);
         }
diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2LEPStop_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2LEPStop_36invfb.cpp
index a5062c83fa..4f295d8588 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2LEPStop_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2LEPStop_36invfb.cpp
@@ -178,7 +178,7 @@ namespace Gambit {
             ATLAS::applyLooseIDElectronSelectionR2(baselineElectrons);
             // Jets
             vector<const HEPUtils::Jet*> baselineJets;
-            for (const HEPUtils::Jet* jet : event->jets()) {
+            for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
                 if (jet->pT() > 30. && fabs(jet->eta()) < 2.4) baselineJets.push_back(jet);
             }
 
diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2LEPsoft_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2LEPsoft_36invfb.cpp
index 011e00d72d..57705b9ba8 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2LEPsoft_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2LEPsoft_36invfb.cpp
@@ -168,7 +168,7 @@ namespace Gambit {
           else if (met > 300. && muon->pT()>3.5 && muon->pT()<30. && fabs(muon->eta())<2.4 && isMu) signalMuons.push_back(muon);
         }
 
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT()>25. && fabs(jet->eta())<2.4) {
            signalJets.push_back(jet);
            if (jet->btag())signalBJets.push_back(jet);
diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_137invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_137invfb.cpp
index 826caeaed0..8875551ee2 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_137invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_137invfb.cpp
@@ -200,7 +200,7 @@ namespace Gambit
 
           // Baseline jets
           vector<const HEPUtils::Jet*> baselineJets;
-          for (const HEPUtils::Jet* jet : event->jets())
+          for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
           {
             if (jet->pT()>25. && fabs(jet->eta())<2.4) baselineJets.push_back(jet);
           }
diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_36invfb.cpp
index 87181c3812..753a2541b7 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_36invfb.cpp
@@ -139,7 +139,7 @@ namespace Gambit {
 
         // Baseline jets
         vector<const HEPUtils::Jet*> baselineJets;
-        for (const HEPUtils::Jet* jet : event->jets())
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
         {
           // We use 25 GeV rather than 35 GeV
           // if (jet->pT()>35. &&fabs(jet->eta())<2.4) baselineJets.push_back(jet);
diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_chargino_stop_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_chargino_stop_36invfb.cpp
index 66aa38e9e5..604b4d6274 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_chargino_stop_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_chargino_stop_36invfb.cpp
@@ -265,7 +265,7 @@ namespace Gambit {
 
         // Baseline jets
         vector<const HEPUtils::Jet*> baselineJets;
-        for (const HEPUtils::Jet* jet : event->jets())
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
         {
           if (jet->pT()>20. &&fabs(jet->eta())<2.4) baselineJets.push_back(jet);
         }
diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_confnote_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_confnote_36invfb.cpp
index fb3e049aaf..62a9a7a4b7 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_confnote_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2OSLEP_confnote_36invfb.cpp
@@ -97,7 +97,7 @@ namespace Gambit {
         }
 
         vector<const HEPUtils::Jet*> baselineJets;
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT()>35. &&fabs(jet->eta())<2.4)baselineJets.push_back(jet);
         }
 
diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2Photon_GMSB_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2Photon_GMSB_36invfb.cpp
index 1f4a282f56..4549ce53c5 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2Photon_GMSB_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2Photon_GMSB_36invfb.cpp
@@ -207,7 +207,7 @@ namespace Gambit {
 
         // Jets
         vector<const HEPUtils::Jet*> jets;
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           // No info on baseline jet cuts in the paper, so for now we'll
           // apply an|eta| cut for HCAL coverage and a loose jet pT cut
           if (jet->pT()>10. && jet->abseta()<3.0) jets.push_back(jet);
diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_137invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_137invfb.cpp
index 9479ffe8e6..e74cab7cbb 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_137invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_137invfb.cpp
@@ -293,10 +293,10 @@ namespace Gambit {
                     muons.push_back(muon);
             }
 
-            double HT = 0.;
             // Jets
+            double HT = 0.;
             vector<const HEPUtils::Jet*> candJets;
-            for (const HEPUtils::Jet* jet : event->jets()) {
+            for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
                 if (jet->pT() > 25. && fabs(jet->eta()) < 2.4){
                     HT += jet->pT();
                     candJets.push_back(jet);
diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_36invfb.cpp
index e6786a68c2..a7e4f8ce1f 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_2SSLEP_Stop_36invfb.cpp
@@ -290,10 +290,10 @@ namespace Gambit {
                     muons.push_back(muon);
             }
 
-            double HT = 0.;
             // Jets
+            double HT = 0.;
             vector<const HEPUtils::Jet*> candJets;
-            for (const HEPUtils::Jet* jet : event->jets()) {
+            for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
                 if (jet->pT() > 25. && fabs(jet->eta()) < 2.4){
                     HT += jet->pT();
                     candJets.push_back(jet);
diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_MONOJET_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_MONOJET_36invfb.cpp
index 1a9388fde1..650820373a 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_MONOJET_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_MONOJET_36invfb.cpp
@@ -91,7 +91,7 @@ namespace Gambit {
 
         // Get jets
         vector<const Jet*> jets4;
-        for (const Jet* jet : event->jets())
+        for (const Jet* jet : event->jets("antikt_R04"))
           if (jet->pT() > 20) jets4.push_back(jet);
 
         // Veto if there are any b-tagged jets (reduce top background)
diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_137invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_137invfb.cpp
index fa8d6a9b73..ab3d3f7bc7 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_137invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_137invfb.cpp
@@ -184,7 +184,7 @@ namespace Gambit
 
           // Baseline jets
           std::vector<const HEPUtils::Jet*> baselineJets;
-          for (const HEPUtils::Jet* jet : event->jets())
+          for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
           {
             if (jet->pT()>25. &&fabs(jet->eta())<2.4)
               baselineJets.push_back(jet);
diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_36invfb.cpp
index 1c9351564a..50eb14f3a0 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_36invfb.cpp
@@ -193,7 +193,7 @@ namespace Gambit
           }
 
           vector<const HEPUtils::Jet*> baselineJets;
-          for (const HEPUtils::Jet* jet : event->jets())
+          for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
           {
             if (jet->pT()>25. &&fabs(jet->eta())<2.4)baselineJets.push_back(jet);
           }
diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_Full_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_Full_36invfb.cpp
index 24b5f195ef..1885609fc5 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_Full_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_MultiLEP_Full_36invfb.cpp
@@ -312,7 +312,7 @@ namespace Gambit {
         }
 
         vector<const HEPUtils::Jet*> baselineJets;
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT()>25. &&fabs(jet->eta())<2.4)baselineJets.push_back(jet);
         }
 
diff --git a/ColliderBit/src/analyses/Analysis_CMS_13TeV_Photon_GMSB_36invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_13TeV_Photon_GMSB_36invfb.cpp
index cd1cfcf68a..3b5e209a19 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_13TeV_Photon_GMSB_36invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_13TeV_Photon_GMSB_36invfb.cpp
@@ -87,7 +87,7 @@ namespace Gambit {
 
         // jets
         vector<const HEPUtils::Jet*> Jets;
-        for (const HEPUtils::Jet* jet : event->jets())
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04"))
         {
           if (jet->pT()>30. &&fabs(jet->eta())<3.0) Jets.push_back(jet);
         }
diff --git a/ColliderBit/src/analyses/Analysis_CMS_8TeV_1LEPDMTOP_20invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_8TeV_1LEPDMTOP_20invfb.cpp
index b10a199891..92fbd513f9 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_8TeV_1LEPDMTOP_20invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_8TeV_1LEPDMTOP_20invfb.cpp
@@ -114,7 +114,7 @@ namespace Gambit {
         const std::vector<double> c = {0.60};
         HEPUtils::BinnedFn2D<double> _eff2d(a,b,c);
 
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT() > 30. && fabs(jet->eta()) < 4.0) {
             baselineJets.push_back(jet);
             //LorentzVector j1 (jet->mom().px(),jet->mom().py(),jet->mom().pz(),jet->mom().E()) ;
diff --git a/ColliderBit/src/analyses/Analysis_CMS_8TeV_2LEPDMTOP_20invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_8TeV_2LEPDMTOP_20invfb.cpp
index 56871097d5..77619a6910 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_8TeV_2LEPDMTOP_20invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_8TeV_2LEPDMTOP_20invfb.cpp
@@ -111,7 +111,7 @@ namespace Gambit {
         const std::vector<double> c = {0.60};
         HEPUtils::BinnedFn2D<double> _eff2d(a,b,c);
 
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT() > 30. && fabs(jet->eta()) < 5.0) {
             baselineJets.push_back(jet);
             //LorentzVector j1 (jet->mom().px(),jet->mom().py(),jet->mom().pz(),jet->mom().E()) ;
diff --git a/ColliderBit/src/analyses/Analysis_CMS_8TeV_MONOJET_20invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_8TeV_MONOJET_20invfb.cpp
index 4b3fe973ca..b577d8d2b7 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_8TeV_MONOJET_20invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_8TeV_MONOJET_20invfb.cpp
@@ -110,7 +110,7 @@ namespace Gambit {
         vector<const HEPUtils::Jet*> baselineJets;
         vector<HEPUtils::P4> jets;
 
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT() > 30. && fabs(jet->eta()) < 4.5) {
             baselineJets.push_back(jet);
           }
diff --git a/ColliderBit/src/analyses/Analysis_CMS_8TeV_MultiLEP_20invfb.cpp b/ColliderBit/src/analyses/Analysis_CMS_8TeV_MultiLEP_20invfb.cpp
index 3c1e31b638..0780ee2dd2 100644
--- a/ColliderBit/src/analyses/Analysis_CMS_8TeV_MultiLEP_20invfb.cpp
+++ b/ColliderBit/src/analyses/Analysis_CMS_8TeV_MultiLEP_20invfb.cpp
@@ -477,7 +477,7 @@ namespace Gambit {
         // - jets
         vector<const HEPUtils::Jet*> signalJets;
         vector<const HEPUtils::Jet*> signalBjets;
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT() > 30. && fabs(jet->eta()) < 2.5) signalJets.push_back(jet);
           if(jet->btag() && fabs(jet->eta()) < 2.5 && jet->pT() > 30.) signalBjets.push_back(jet);
         }
diff --git a/ColliderBit/src/analyses/Analysis_Minimum.cpp b/ColliderBit/src/analyses/Analysis_Minimum.cpp
index 4044c022a3..683332bfa4 100644
--- a/ColliderBit/src/analyses/Analysis_Minimum.cpp
+++ b/ColliderBit/src/analyses/Analysis_Minimum.cpp
@@ -56,7 +56,7 @@ namespace Gambit {
 
         // Baseline jets
         vector<const HEPUtils::Jet*> baselineJets;
-        for (const HEPUtils::Jet* jet : event->jets()) {
+        for (const HEPUtils::Jet* jet : event->jets("antikt_R04")) {
           if (jet->pT() > 20. && fabs(jet->eta()) < 4.5) baselineJets.push_back(jet);
         }
 
diff --git a/ColliderBit/src/detectors/BuckFast.cpp b/ColliderBit/src/detectors/BuckFast.cpp
index 2d2f12dd7c..70a945b2c7 100644
--- a/ColliderBit/src/detectors/BuckFast.cpp
+++ b/ColliderBit/src/detectors/BuckFast.cpp
@@ -39,12 +39,21 @@ namespace Gambit
       if (smearTaus != NULL) smearTaus(event.taus());
 
       // Smear jet momenta
-      if (smearJets != NULL) smearJets(event.jets());
+      if (smearJets != NULL)
+      {
+        for (std::string jetcollection : event.jet_collections())
+        {
+          smearJets(event.jets(jetcollection));
+        }
+      }
 
       // Unset b-tags outside |eta|=2.5
-      for (HEPUtils::Jet* j : event.jets())
+      for (std::string jetcollection : event.jet_collections())
       {
-        if (j->abseta() > 2.5) j->set_btag(false);
+        for (HEPUtils::Jet* j : event.jets(jetcollection))
+        {
+          if (j->abseta() > 2.5) j->set_btag(false);
+        }
       }
     }
 
diff --git a/ColliderBit/src/getHepMCEvent.cpp b/ColliderBit/src/getHepMCEvent.cpp
index 7018ef65e1..67a92c0429 100644
--- a/ColliderBit/src/getHepMCEvent.cpp
+++ b/ColliderBit/src/getHepMCEvent.cpp
@@ -190,6 +190,44 @@ namespace Gambit
 
     }
 
+
+    /// A helper function for collecting the jet_collections yaml settings. 
+    /// Used by functions getHepMCEvent_HEPUtils and convertHepMCEvent_HEPUtils.
+    void read_jet_collections_settings(const Options& runOptions, std::vector<jet_collection_settings>& all_jet_collection_settings, str& jetcollection_taus)
+    {
+      if (runOptions.hasKey("jet_collections"))
+      {
+        YAML::Node all_jetcollections_node = runOptions.getValue<YAML::Node>("jet_collections");
+        Options all_jetcollection_options(all_jetcollections_node);
+        std::vector<str> jetcollection_names = all_jetcollection_options.getNames();
+
+        for (str key : jetcollection_names)
+        {
+          YAML::Node current_jc_node = all_jetcollection_options.getValue<YAML::Node>(key);
+          Options current_jc_options(current_jc_node);
+
+          str algorithm = current_jc_options.getValue<str>("algorithm");
+          double R = current_jc_options.getValue<double>("R");
+          str recombination_scheme = current_jc_options.getValue<str>("recombination_scheme");
+          str strategy = current_jc_options.getValue<str>("strategy");
+
+          all_jet_collection_settings.push_back({key, algorithm, R, recombination_scheme, strategy});
+        }
+
+        jetcollection_taus = runOptions.getValue<str>("jet_collection_taus");
+        // Throw an error if the "jet_collection_taus" setting does not match an entry in "jet_collections".
+        if (std::find(jetcollection_names.begin(), jetcollection_names.end(), jetcollection_taus) == jetcollection_names.end())
+        {
+          ColliderBit_error().raise(LOCAL_INFO,"Please provide the jet_collection_taus setting for jet collections.");
+        }
+      }
+      else
+      {
+        ColliderBit_error().raise(LOCAL_INFO,"Could not find jet_collections options. Please provide this in the YAML file.");
+      }
+    }
+
+
     /// A nested function that reads in HepMC event files and converts them to HEPUtils::Event format
     void getHepMCEvent_HEPUtils(HEPUtils::Event &result)
     {
@@ -197,8 +235,10 @@ namespace Gambit
 
       // Get yaml options
       const static str HepMC_filename = runOptions->getValueOrDef<str>("", "hepmc_filename");
-      const static double antiktR = runOptions->getValueOrDef<double>(0.4, "antiktR");
       const static double jet_pt_min = runOptions->getValueOrDef<double>(10.0, "jet_pt_min");
+      std::vector<jet_collection_settings> all_jet_collection_settings = {};
+      str jetcollection_taus;
+      read_jet_collections_settings(*runOptions, all_jet_collection_settings, jetcollection_taus);
 
       // Get the HepMC event
       //HepMC3::GenEvent ge = *Dep::HardScatteringEvent;
@@ -212,7 +252,7 @@ namespace Gambit
       result.set_weight(ge.weight());
 
       //Translate to HEPUtils event by calling the unified HEPMC/Pythia event converter:
-      Gambit::ColliderBit::convertParticleEvent(ge.particles(), result, antiktR, jet_pt_min);
+      Gambit::ColliderBit::convertParticleEvent(ge.particles(), result, all_jet_collection_settings, jetcollection_taus, jet_pt_min);
 
     }
 
@@ -227,14 +267,16 @@ namespace Gambit
       HepMC3::GenEvent ge = *Dep::HardScatteringEvent;
 
       //Get yaml options required for conversion
-      const static double antiktR = runOptions->getValueOrDef<double>(0.4, "antiktR");
       const static double jet_pt_min = runOptions->getValueOrDef<double>(10.0, "jet_pt_min");
+      std::vector<jet_collection_settings> all_jet_collection_settings = {};
+      str jetcollection_taus;
+      read_jet_collections_settings(*runOptions, all_jet_collection_settings, jetcollection_taus);
 
       //Set the weight
       result.set_weight(ge.weight());
 
       //Translate to HEPUtils event by calling the unified HEPMC/Pythia event converter:
-      Gambit::ColliderBit::convertParticleEvent(ge.particles(), result, antiktR, jet_pt_min);
+      Gambit::ColliderBit::convertParticleEvent(ge.particles(), result, all_jet_collection_settings, jetcollection_taus, jet_pt_min);
     }
 
   }
diff --git a/ColliderBit/src/getLHEvent.cpp b/ColliderBit/src/getLHEvent.cpp
index c21dc95d30..b62c87e362 100644
--- a/ColliderBit/src/getLHEvent.cpp
+++ b/ColliderBit/src/getLHEvent.cpp
@@ -55,6 +55,54 @@ namespace Gambit
       }
       static LHEF::Reader lhe(lhef_filename);
 
+      // Get all jet collection settings
+      str jetcollection_taus;
+      std::vector<jet_collection_settings> all_jet_collection_settings = {};
+      if (runOptions->hasKey((*Dep::RunMC).current_collider()))
+      {
+        YAML::Node colNode = runOptions->getValue<YAML::Node>((*Dep::RunMC).current_collider());
+        Options colOptions(colNode);
+
+        // Fill the jet collection settings
+        if (colOptions.hasKey("jet_collections"))
+        {
+          YAML::Node all_jetcollections_node = colOptions.getValue<YAML::Node>("jet_collections");
+          Options all_jetcollection_options(all_jetcollections_node);
+          std::vector<str> jetcollection_names = all_jetcollection_options.getNames();
+
+          for (str key : jetcollection_names)
+          {
+            YAML::Node current_jc_node = all_jetcollection_options.getValue<YAML::Node>(key);
+            Options current_jc_options(current_jc_node);
+
+            str algorithm = current_jc_options.getValue<str>("algorithm");
+            double R = current_jc_options.getValue<double>("R");
+            str recombination_scheme = current_jc_options.getValue<str>("recombination_scheme");
+            str strategy = current_jc_options.getValue<str>("strategy");
+
+            all_jet_collection_settings.push_back({key, algorithm, R, recombination_scheme, strategy});
+          }
+
+          // Note that jetcollection_taus is not used here as get_HEPUtils_event(...) has a much simpler jet definition than in Py8Conversions.hpp
+          jetcollection_taus = colOptions.getValue<str>("jet_collection_taus");
+          // Throw an error if the "jet_collection_taus" setting does not match an entry in "jet_collections".
+          if (std::find(jetcollection_names.begin(), jetcollection_names.end(), jetcollection_taus) == jetcollection_names.end())
+          {
+            ColliderBit_error().raise(LOCAL_INFO,"Please provide the jet_collection_taus setting for jet collections.");
+          }
+        }
+        else
+        {
+          str current_collider = (*Dep::RunMC).current_collider();
+          ColliderBit_error().raise(LOCAL_INFO,"Could not find jet_collections option for collider " + current_collider +  ". Please provide this in the YAML file.");
+        }
+      }
+      else
+      {
+        str current_collider = (*Dep::RunMC).current_collider();
+        ColliderBit_error().raise(LOCAL_INFO,"Could not find runOptions for collider " + current_collider + ".");
+      }
+
       // Don't do anything during special iterations
       if (*Loop::iteration < 0) return;
 
@@ -62,7 +110,7 @@ namespace Gambit
       bool event_retrieved = true;
       #pragma omp critical (reading_LHEvent)
       {
-        if (lhe.readEvent()) get_HEPUtils_event(lhe, result, jet_pt_min);
+        if (lhe.readEvent()) get_HEPUtils_event(lhe, result, jet_pt_min, all_jet_collection_settings);
         else event_retrieved = false;
       }
       if (not event_retrieved)
diff --git a/ColliderBit/src/lhef2heputils.cpp b/ColliderBit/src/lhef2heputils.cpp
index b267d279db..cbbd44d6cb 100644
--- a/ColliderBit/src/lhef2heputils.cpp
+++ b/ColliderBit/src/lhef2heputils.cpp
@@ -44,86 +44,101 @@ using namespace std;
 using namespace HEPUtils;
 using namespace FJNS;
 
-/// Extract an LHE event as a HEPUtils::Event
-void get_HEPUtils_event(const LHEF::Reader& lhe, Event& evt, double jet_pt_min)
+namespace Gambit
 {
 
-  P4 vmet;
-  vector<PseudoJet> jetparticles;
-
-  evt.set_weight(lhe.hepeup.weight());
-
-  // Loop over all particles in the event
-  for (int i = 0; i < lhe.hepeup.NUP; ++i)
+  namespace ColliderBit
   {
-    // Get status and PID code
-    const int st = lhe.hepeup.ISTUP[i];
-    const int pid = lhe.hepeup.IDUP[i];
-    const int apid = fabs(pid);
-
-    // Use LHE-stable particles only
-    if (st != 1) continue;
-
-    // Get 4-momentum
-    const P4 p4 = P4::mkXYZM(lhe.hepeup.PUP[i][0], lhe.hepeup.PUP[i][1], lhe.hepeup.PUP[i][2], lhe.hepeup.PUP[i][4]);
 
-    // Store interacting prompt particles
-    /// @todo Dress leptons?
-    if (apid == 22 || apid == 11 || apid == 13 || apid == 15)
+    /// Extract an LHE event as a HEPUtils::Event
+    void get_HEPUtils_event(const LHEF::Reader& lhe, Event& evt, double jet_pt_min, std::vector<jet_collection_settings> all_jet_collection_settings)
     {
-      Particle* p = new Particle(p4, pid); // the event will take ownership of this pointer
-      p->set_prompt(true);
-      evt.add_particle(p);
-    }
 
-    // Aggregate missing ET
-    else if (apid == 12 || apid == 14 || apid == 16 || apid == 1000022 || apid == 1000039)
-    {
-      Particle* p = new Particle(p4, pid); // the event will take ownership of this pointer
-      p->set_prompt(true);
-      evt.add_particle(p);
-      vmet += p4;
-    }
+      P4 vmet;
+      vector<PseudoJet> jetparticles;
+
+      evt.set_weight(lhe.hepeup.weight());
+
+      // Loop over all particles in the event
+      for (int i = 0; i < lhe.hepeup.NUP; ++i)
+      {
+        // Get status and PID code
+        const int st = lhe.hepeup.ISTUP[i];
+        const int pid = lhe.hepeup.IDUP[i];
+        const int apid = fabs(pid);
+
+        // Use LHE-stable particles only
+        if (st != 1) continue;
+
+        // Get 4-momentum
+        const P4 p4 = P4::mkXYZM(lhe.hepeup.PUP[i][0], lhe.hepeup.PUP[i][1], lhe.hepeup.PUP[i][2], lhe.hepeup.PUP[i][4]);
+
+        // Store interacting prompt particles
+        /// @todo Dress leptons?
+        if (apid == 22 || apid == 11 || apid == 13 || apid == 15)
+        {
+          Particle* p = new Particle(p4, pid); // the event will take ownership of this pointer
+          p->set_prompt(true);
+          evt.add_particle(p);
+        }
+
+        // Aggregate missing ET
+        else if (apid == 12 || apid == 14 || apid == 16 || apid == 1000022 || apid == 1000039)
+        {
+          Particle* p = new Particle(p4, pid); // the event will take ownership of this pointer
+          p->set_prompt(true);
+          evt.add_particle(p);
+          vmet += p4;
+        }
+
+        // Store non-prompt momenta for Jet building
+        else
+        {
+          PseudoJet pj = mk_pj(p4);
+          pj.set_user_index(apid);
+          jetparticles.push_back(pj);
+        }
+      }
+
+      // MET
+      vmet.setPz(0);
+      vmet.setM(0);
+      evt.set_missingmom(vmet);
+
+      // Jet Finding
+      for (jet_collection_settings jetcollection : all_jet_collection_settings)
+      {
+
+        // @todo get_jets function could accept a more general jet definition
+        vector<PseudoJet> jets = HEPUtils::get_jets(jetparticles, 0.4, jet_pt_min, FJalgorithm_map(jetcollection.algorithm));
+
+        for (const PseudoJet& pj : jets)
+        {
+          bool hasC = false, hasB = false;
+          /// @todo Bug in HEPUtils::get_jets means that constituent info is lost for now...
+          // for (const PseudoJet& c : pj.constituents()) {
+          //   if (c.user_index() == 4) hasC = true;
+          //   if (c.user_index() == 5) hasB = true;
+          // }
+          evt.add_jet(new Jet(mk_p4(pj), hasB, hasC), jetcollection.key);
+        }
+
+        #ifdef COLLIDERBIT_DEBUG
+          // Print event summary
+          cout << "  MET  = " << evt.met() << " GeV" << endl;
+          cout << "  #e   = " << evt.electrons().size() << endl;
+          cout << "  #mu  = " << evt.muons().size() << endl;
+          cout << "  #tau = " << evt.taus().size() << endl;
+          cout << "  #jet = " << evt.jets().size() << endl;
+          cout << "  #pho  = " << evt.photons().size() << endl;
+          cout << endl;
+        #endif
+      }
 
-    // Store non-prompt momenta for Jet building
-    else
-    {
-      PseudoJet pj = mk_pj(p4);
-      pj.set_user_index(apid);
-      jetparticles.push_back(pj);
     }
-  }
-
-  // MET
-  vmet.setPz(0);
-  vmet.setM(0);
-  evt.set_missingmom(vmet);
 
-  // Jets
-  vector<PseudoJet> jets = get_jets(jetparticles, 0.4, jet_pt_min);
-
-  for (const PseudoJet& pj : jets)
-  {
-    bool hasC = false, hasB = false;
-    /// @todo Bug in HEPUtils::get_jets means that constituent info is lost for now...
-    // for (const PseudoJet& c : pj.constituents()) {
-    //   if (c.user_index() == 4) hasC = true;
-    //   if (c.user_index() == 5) hasB = true;
-    // }
-    evt.add_jet(new Jet(mk_p4(pj), hasB, hasC));
   }
 
-  #ifdef COLLIDERBIT_DEBUG
-    // Print event summary
-    cout << "  MET  = " << evt.met() << " GeV" << endl;
-    cout << "  #e   = " << evt.electrons().size() << endl;
-    cout << "  #mu  = " << evt.muons().size() << endl;
-    cout << "  #tau = " << evt.taus().size() << endl;
-    cout << "  #jet = " << evt.jets().size() << endl;
-    cout << "  #pho  = " << evt.photons().size() << endl;
-    cout << endl;
-  #endif
-
 }
 
-#endif
\ No newline at end of file
+#endif
diff --git a/contrib/heputils/include/HEPUtils/BinnedFn.h b/contrib/heputils/include/HEPUtils/BinnedFn.h
index 6e52ee473b..72558e1051 100644
--- a/contrib/heputils/include/HEPUtils/BinnedFn.h
+++ b/contrib/heputils/include/HEPUtils/BinnedFn.h
@@ -1,7 +1,7 @@
 // -*- C++ -*-
 //
-// This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils
-// Copyright (C) 2013-2018 Andy Buckley <andy.buckley@cern.ch>
+// This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/
+// Copyright (C) 2013-2022 Andy Buckley <andy.buckley@cern.ch>
 //
 // Embedding of HEPUtils code in other projects is permitted provided this
 // notice is retained and the HEPUtils namespace and include path are changed.
diff --git a/contrib/heputils/include/HEPUtils/Event.h b/contrib/heputils/include/HEPUtils/Event.h
index 987e11832b..acc29dfa23 100644
--- a/contrib/heputils/include/HEPUtils/Event.h
+++ b/contrib/heputils/include/HEPUtils/Event.h
@@ -1,7 +1,7 @@
 // -*- C++ -*-
 //
-// This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils
-// Copyright (C) 2013-2018 Andy Buckley <andy.buckley@cern.ch>
+// This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/
+// Copyright (C) 2013-2023 Andy Buckley <andy.buckley@cern.ch>
 //
 // Embedding of HEPUtils code in other projects is permitted provided this
 // notice is retained and the HEPUtils namespace and include path are changed.
@@ -11,6 +11,8 @@
 #include "HEPUtils/Particle.h"
 #include "HEPUtils/Jet.h"
 #include <algorithm>
+#include <memory>
+#include <map>
 
 namespace HEPUtils {
 
@@ -20,26 +22,34 @@ namespace HEPUtils {
   private:
 
     /// @name Internal particle / vector containers
-    //@{
+    /// @{
 
     /// Event weights
     std::vector<double> _weights;
     std::vector<double> _weight_errs;
 
     /// @name Separate particle collections
-    //@{
-    std::vector<Particle*> _photons, _electrons, _muons, _taus, _invisibles;
-    std::vector<const Particle*> _cphotons, _celectrons, _cmuons, _ctaus, _cinvisibles;
-    //@}
-
-    /// Jets collection (mutable to allow sorting)
-    mutable std::vector<Jet*> _jets;
-    mutable std::vector<const Jet*> _cjets;
-
+    /// @{
+    std::vector<const Particle*> _allparticles;
+    mutable std::vector<const Particle*> _visibles, _invisibles, _photons, _electrons, _muons, _taus;
+    /// @}
+
+    /// Jets collection(s) (mutable to allow sorting)
+    mutable std::map<std::string, std::vector<const Jet*>> _jets;
+
+    /// Typedef for the generic cluster-sequence type
+    using CSeqBase = FJNS::ClusterSequence;
+    /// Typedef for a smart ptr to the generic cluster-sequence type
+    using CSeqBasePtr = std::shared_ptr<const FJNS::ClusterSequence>;
+    
+    /// Hold the cluster sequences corresponding to jets, to keep them alive
+    std::map<const std::string, CSeqBasePtr> _cseqs;
+    
     /// Missing momentum vector
     P4 _pmiss;
 
-    //@}
+    /// @}
+
 
   private:
 
@@ -49,27 +59,23 @@ namespace HEPUtils {
       clear(); //< Delete current particles
       _weights = e._weights;
       _weight_errs = e._weight_errs;
+      //
+      _allparticles = e._allparticles;
+      _visibles = e._visibles;
+      _invisibles = e._invisibles;
       _photons = e._photons;
-      _cphotons = e._cphotons;
       _electrons = e._electrons;
-      _celectrons = e._celectrons;
       _muons = e._muons;
-      _cmuons = e._cmuons;
       _taus = e._taus;
-      _ctaus = e._ctaus;
-      _invisibles = e._invisibles;
-      _cinvisibles = e._cinvisibles;
+      //
       _jets = e._jets;
-      _cjets = e._cjets;
+      _cseqs = e._cseqs;
       _pmiss = e._pmiss;
     }
 
 
   public:
 
-    /// @name Constructors
-    //@{
-
     /// Default constructor
     Event() { clear(); }
 
@@ -89,7 +95,8 @@ namespace HEPUtils {
     }
 
 
-  public:
+    /// @name Cloning (= deep copy)
+    /// @{
 
     /// Clone a copy on the heap
     Event* clone() const {
@@ -112,59 +119,58 @@ namespace HEPUtils {
       for (size_t i = 0; i < ps.size(); ++i) {
         e.add_particle(new Particle(*ps[i]));
       }
-      const std::vector<const Jet*> js = jets();
-      for (size_t i = 0; i < js.size(); ++i) {
-        e.add_jet(new Jet(*js[i]));
+      for (const auto& kv : _jets ) {
+        const std::vector<const Jet*> js = jets(kv.first);
+        for (size_t i = 0; i < js.size(); ++i) {
+          e.add_jet(new Jet(*js[i]), kv.first);
+        }
       }
       e._pmiss = _pmiss;
+      e._cseqs = _cseqs;
     }
 
-    //@}
+    /// @}
 
 
-    /// Empty the event's particle, jet and MET collections
+    /// Empty the event's weight, particle, jet, and MET collections
     void clear() {
+      // Weights
       _weights.clear();
       _weight_errs.clear();
-      // TODO: indexed loop -> for (Particle* p : particles()) delete p;
-      #define DELCLEAR(v) do { if (!v.empty()) for (size_t i = 0; i < v.size(); ++i) delete v[i]; v.clear(); } while (0)
-      DELCLEAR(_photons);
-      DELCLEAR(_electrons);
-      DELCLEAR(_muons);
-      DELCLEAR(_taus);
-      DELCLEAR(_invisibles);
-      DELCLEAR(_jets);
-      #undef DELCLEAR
 
+      // Particles, first the canonical collection
+      for (const Particle* p : _allparticles) delete p;
+      _allparticles.clear();
+      // Then the caches
+      _visibles.clear();
+      _invisibles.clear();
       _photons.clear();
-      _cphotons.clear();
       _electrons.clear();
-      _celectrons.clear();
       _muons.clear();
-      _cmuons.clear();
       _taus.clear();
-      _ctaus.clear();
-      _invisibles.clear();
-      _cinvisibles.clear();
+
+      // Jets
+      for (const std::string& jc : jet_collections()) clear_jets(jc);
       _jets.clear();
-      _cjets.clear();
+      _cseqs.clear();
 
+      // MET
       _pmiss.clear();
     }
 
 
-    ///////////////////////
-
+    /// @name Weights
+    /// @{
 
     /// Set the event weights (also possible directly via non-const reference)
     void set_weights(const std::vector<double>& ws) {
       _weights = ws;
     }
 
+    /// Set the event weight errors (also possible directly via non-const reference)
     void set_weight_errs(const std::vector<double>& werrs) {
       _weight_errs = werrs;
     }
-
     /// Set the event weights to the single given weight
     void set_weight(double w) {
       _weights.clear();
@@ -196,7 +202,6 @@ namespace HEPUtils {
     std::vector<double>& weight_errs() {
       return _weight_errs;
     }
-
     /// Get a single event weight -- the nominal, by default
     double weight(size_t i=0) const {
       if (_weights.empty()) {
@@ -215,11 +220,13 @@ namespace HEPUtils {
       return _weight_errs[i];
     }
 
+    /// @}
 
-    /////////////////
 
+    /// @name Particles
+    /// @{
 
-    /// Add a particle to the event
+    /// @brief Add a particle to the event
     ///
     /// Supplied particle should be new'd, and Event will take ownership.
     ///
@@ -228,190 +235,265 @@ namespace HEPUtils {
     /// immediately deleted. Accordingly, the pointer passed by user code
     /// must be considered potentially invalid from the moment this function is called.
     ///
+    /// @note pT-sorting can be disabled, so
+    ///
     /// @todo "Lock" at some point so that jet finding etc. only get done once
-    void add_particle(Particle* p) {
-      if (!p->is_prompt())
-        delete p;
-      else if (p->pid() == 22)
-      {
-        _photons.push_back(p);
-        _cphotons.push_back(p);
-      }
-      else if (p->abspid() == 11)
-      {
-        _electrons.push_back(p);
-        _celectrons.push_back(p);
+    void add_particle(Particle* p, bool ptsort=true) {
+
+      // All particles (canonical collection)
+      _allparticles.push_back(p);
+
+      // Caching collections
+      if (!p->is_visible()) {
+        if (p->is_prompt()) _invisibles.push_back(p);
+      } else {
+        _visibles.push_back(p);
+        if (p->is_prompt()) {
+          if (p->pid() == 22) _photons.push_back(p);
+          else if (p->abspid() == 11) _electrons.push_back(p);
+          else if (p->abspid() == 13) _muons.push_back(p);
+          else if (p->abspid() == 15) _taus.push_back(p);
+        }
       }
-      else if (p->abspid() == 13)
-      {
-        _muons.push_back(p);
-        _cmuons.push_back(p);
-      }
-      else if (p->abspid() == 15)
-      {
-        _taus.push_back(p);
-        _ctaus.push_back(p);
-      }
-      else if (p->abspid() == 12 || p->abspid() == 14 || p->abspid() == 16 ||
-               p->pid() == 1000022 || p->pid() == 1000039 ||
-               in_range(p->abspid(), 50, 60)) //< invert definition to specify all *visibles*?
-      {
-        _invisibles.push_back(p);
-        _cinvisibles.push_back(p);
-      }
-      else
-        delete p;
+
+      // Sort the collections
+      if (ptsort) _sort_particles();
     }
 
 
     /// Add a collection of final state particles to the event
     ///
-    /// Supplied particles should be new'd, and Event will take ownership.
-    void add_particles(const std::vector<Particle*>& ps) {
-      for (size_t i = 0; i < ps.size(); ++i) add_particle(ps[i]);
-    }
-
-
-    /// Get all final state particles
-    /// @todo Note the return by value: it's not efficient yet!
-    /// @note Overlap of taus and e/mu
-    std::vector<const Particle*> particles() const {
-      // Add together all the vectors of the different particle types
-      std::vector<const Particle*> rtn;
-      // rtn.reserve(visible_particles().size() + _cinvisibles.size());
-      rtn.reserve(_cphotons.size() + _celectrons.size() + _cmuons.size() + _ctaus.size() + _cinvisibles.size());
-      #define APPEND_VEC(vec) rtn.insert(rtn.end(), vec.begin(), vec.end())
-      // APPEND_VEC(visible_particles());
-      APPEND_VEC(_cphotons);
-      APPEND_VEC(_celectrons);
-      APPEND_VEC(_cmuons);
-      APPEND_VEC(_ctaus);
-      APPEND_VEC(_cinvisibles);
-      #undef APPEND_VEC
-      return rtn;
-      /// @todo Or use Boost range join to iterate over separate containers transparently... I like this
-      /// @todo Cache, or otherwise think about efficiency since this gets called by the destructor
-    }
-
-
-    /// Get visible state particles
-    /// @todo Note the return by value: it's not efficient yet!
-    /// @note Overlap of taus and e/mu
-    std::vector<const Particle*> visible_particles() const {
-      // Add together all the vectors of the different particle types
-      std::vector<const Particle*> rtn;
-      rtn.reserve(_cphotons.size() + _celectrons.size() + _cmuons.size() + _ctaus.size());
-      #define APPEND_VEC(vec) rtn.insert(rtn.end(), vec.begin(), vec.end() )
-      APPEND_VEC(_cphotons);
-      APPEND_VEC(_celectrons);
-      APPEND_VEC(_cmuons);
-      APPEND_VEC(_ctaus);
-      #undef APPEND_VEC
-      return rtn;
-      /// @todo Add together all the vectors of the different particle types
-      /// @todo Or use Boost range join to iterate over separate containers transparently... I like this
+    /// @warning Supplied particles should be new'd, and Event will take ownership.
+    void add_particles(const std::vector<Particle*>& ps, bool ptsort=true) {
+      // Add each particle, without sorting each time
+      for (const Particle* p : ps) add_particle(new Particle(*p), false);
+
+      // Finally sort the collections, once all new particles are in place
+      if (ptsort) _sort_particles();
+    }
+
+
+    /// A mostly-internal function to sort the particle-vector caches
+    void _sort_particles() {
+      std::sort(_allparticles.begin(), _allparticles.end(), _cmpPtDescPtr<Particle>);
+      std::sort(_invisibles.begin(), _invisibles.end(), _cmpPtDescPtr<Particle>);
+      std::sort(_photons.begin(), _photons.end(), _cmpPtDescPtr<Particle>);
+      std::sort(_electrons.begin(), _electrons.end(), _cmpPtDescPtr<Particle>);
+      std::sort(_muons.begin(), _muons.end(), _cmpPtDescPtr<Particle>);
+      std::sort(_taus.begin(), _taus.end(), _cmpPtDescPtr<Particle>);
+    }
+
+
+    /// @brief Get all final state particles
+    ///
+    /// @note Small overlap of taus and e/mu?
+    const std::vector<const Particle*>& particles() const {
+      return _allparticles;
+    }
+
+
+    /// @brief Get visible state particles
+    ///
+    /// @note Small overlap of taus and e/mu?
+    const std::vector<const Particle*>& visible_particles() const {
+      return _visibles;
     }
 
 
-    /// Get invisible final state particles
+    /// @brief Get invisible final state particles
+    ///
+    /// @note Both prompt and non-prompt... correct?
     const std::vector<const Particle*>& invisible_particles() const {
-      return _cinvisibles;
+      return _invisibles;
     }
     /// Get invisible final state particles (non-const)
     std::vector<Particle*>& invisible_particles() {
-      return _invisibles;
+      return mkunconst(_invisibles);
     }
 
 
     /// Get prompt electrons
     const std::vector<const Particle*>& electrons() const {
-      return _celectrons;
+      return _electrons;
     }
     /// Get prompt electrons (non-const)
     std::vector<Particle*>& electrons() {
-      return _electrons;
+      return mkunconst(_electrons);
     }
 
 
     /// Get prompt muons
     const std::vector<const Particle*>& muons() const {
-      return _cmuons;
+      return _muons;
     }
     /// Get prompt muons (non-const)
     std::vector<Particle*>& muons() {
-      return _muons;
+      return mkunconst(_muons);
     }
 
 
     /// Get prompt (hadronic) taus
     const std::vector<const Particle*>& taus() const {
-      return _ctaus;
+      return _taus;
     }
     /// Get prompt (hadronic) taus (non-const)
     std::vector<Particle*>& taus() {
-      return _taus;
+      return mkunconst(_taus);
     }
 
 
     /// Get prompt photons
     const std::vector<const Particle*>& photons() const {
-      return _cphotons;
+      return _photons;
     }
     /// Get prompt photons (non-const)
     std::vector<Particle*>& photons() {
-      return _photons;
+      return mkunconst(_photons);
     }
 
+    /// @}
+
 
     /// @name Jets
-    /// @todo Why the new'ing? Can we use references more?
-    //@{
+    ///
+    /// @{
 
-    /// @brief Get anti-kT 0.4 jets (not including charged leptons or photons)
-    const std::vector<const Jet*>& jets() const {
-      std::sort(_cjets.begin(), _cjets.end(), _cmpPtDesc);
-      return _cjets;
+    // Implementation function, to avoid const/non-const duplication below
+    std::vector<const Jet*>& _get_jets(const std::string& key) const {
+      std::vector<const Jet*>& rtn = _jets[key];
+      // std::sort(rtn.begin(), rtn.end(), _cmpPtDescPtr<Jet>); //< not thread-safe!
+      return rtn;
+    }
+
+    /// @brief Get a jet collection (not including charged leptons or photons)
+    const std::vector<const Jet*>& jets(const std::string& key) const {
+      return _get_jets(key);
     }
 
-    /// @brief Get anti-kT 0.4 jets (not including charged leptons or photons) (non-const)
-    std::vector<Jet*>& jets()  {
-      std::sort(_jets.begin(), _jets.end(), _cmpPtDesc);
-      return _jets;
+    /// @brief Get a jet collection (not including charged leptons or photons) (non-const)
+    std::vector<Jet*>& jets(const std::string& key) {
+      return mkunconst(_get_jets(key));
     }
 
-    /// @brief Set the jets collection
+
+    /// Get the list of jet-collection names
+    std::vector<std::string> jet_collections() {
+      std::vector<std::string> rtn;
+      for (const auto& kv : _jets) rtn.push_back(kv.first);
+      return rtn;
+    }
+
+    
+    /// @brief Set a jet collection
+    ///
+    /// @warning The Jets should be new'd; Event will take ownership.
+    ///
+    /// @todo "Lock" at some point so that jet finding etc. only get done once
     ///
-    /// The Jets should be new'd; Event will take ownership.
-    void set_jets(const std::vector<Jet*>& jets) {
-      _jets = jets;
-      _cjets.clear();
-      for (Jet* j : jets ) _cjets.push_back(j);
+    /// @note This resets the cluster sequence, but as a shared_ptr is used for storage,
+    ///   any existing shared_ptr links to the previous one will keep it alive.
+    void set_jets(const std::vector<const Jet*>& jets, const std::string& key) {
+      _jets[key] = jets;
+      std::sort(_jets[key].begin(), _jets[key].end(), _cmpPtDescPtr<Jet>);
+      set_clusterseq(jets.front()->clusterseq(), key);
     }
+    // /// @brief Set the jets collection (non-const input)
+    // void set_jets(const std::vector<Jet*>& jets, const std::string& key) {
+    //   set_jets(mkconst(jets), key);
+    // }
+
 
-    /// @brief Add a jet to the jets collection
+    /// @brief Clear a jet collection
     ///
-    /// The Jet should be new'd; Event will take ownership.
-    void add_jet(Jet* j) {
-      _jets.push_back(j);
-      _cjets.push_back(j);
+    /// @note This resets the cluster sequence, but as a shared_ptr is used for storage,
+    ///   any existing shared_ptr links to the previous one will keep it alive.
+    void clear_jets(const std::string& key) {
+      for (const Jet* j : jets(key)) delete j;
+      _jets.erase(key);
+      _cseqs.erase(key);
     }
+    
 
-    //@}
+    /// @brief Add a jet to a jet collection
+    ///
+    /// @warning The Jet should be new'd; Event will take ownership.
+    ///
+    /// @todo "Lock" at some point so that jet finding etc. only get done once
+    void add_jet(const Jet* j, const std::string& key) {
+      _jets[key].push_back(j);
+      std::sort(_jets[key].begin(), _jets[key].end(), _cmpPtDescPtr<Jet>);
+      // Check that the CSeq on the added jet is consistent with this collection
+      /// @todo Needs more care that the cseq pointer is live
+      // if (_cseqs.find(key) != _cseqs.end() && _cseqs.at(key) &&
+      //     !_jets.at(key).empty() && !_jets.at(key).front()->clusterseq() &&
+      //     _cseqs.at(key) != _jets.at(key).front()->clusterseq()) {
+      //   throw std::runtime_error("Event::add_jet() received a jet whose cluster sequence mismatched the active one for that collection");
+      // }
+    }
+    // /// @brief Add a jet to the jets collection (non-const input)
+    // void add_jet(Jet* j, const std::string& key) {
+    //   add_jet(mkconst(j), key);
+    // }
+
+    
+    /// @brief Access the jets' ClusterSequence object if possible (can be null)
+    ///
+    /// Optional template arg can be used to cast to a specific derived CS type if wanted.
+    template <typename CS=FJNS::ClusterSequence>
+    typename std::shared_ptr<const CS> clusterseq(const std::string& key) const {
+      return std::dynamic_pointer_cast<const CS>(_cseqs.find(key)->second);
+    }
 
+    // /// @brief Non-const access to the jets' ClusterSequence object if possible (can be null)
+    // ///
+    // /// Optional template arg can be used to cast to a specific derived CS type if wanted.
+    // template <typename CS=FJNS::ClusterSequence>
+    // std::shared_ptr<CS> clusterseq(const std::string& key) {
+    //   return std::dynamic_pointer_cast<CS>(_cseq->find(key).second);
+    // }
 
-    /// @name Missing energy
-    //@{
+    /// Set a specific cluster sequence for the named jet collection (which must currently be empty)
+    ///
+    /// @warning The CS should be new'd; Event will take ownership via a shared_ptr
+    template <typename CS=FJNS::ClusterSequence>
+    void set_clusterseq(std::shared_ptr<const CS> cseq, const std::string& key) {
+      if (_cseqs.find(key) != _cseqs.end() && !_cseqs.empty()) {
+	throw std::runtime_error("Event::set_clusterseq() called for a non-empty jet collection");
+      }
+      _cseqs[key] =  cseq;
+    }
 
-    /// @brief Get the missing momentum vector
+    
+    /// Create and run a new cluster sequence for the named jet collection
+    ///
+    /// @warning The CS will be new'd, and Event will take ownership via a shared_ptr
+    ///
+    /// @note The resulting pseudojets from this still need to be manually set as Jets.
     ///
-    /// Not _necessarily_ the sum over momenta of final state invisibles
+    /// @todo How to run a more advanced CS like the active- or Voronoi-area ones?
+    template <typename CS=FJNS::ClusterSequence>
+    CSeqBasePtr emplace_clusterseq(std::vector<FJNS::PseudoJet>& jetparticles, const FJNS::JetDefinition& jetdef, const std::string& key) {
+      if (_cseqs.find(key) != _cseqs.end() && !_cseqs.empty()) {
+	throw std::runtime_error("Event::emplace_clusterseq() called for a non-empty jet collection");
+      }
+      _cseqs[key] = std::make_shared<CS>(jetparticles, jetdef);
+      return _cseqs[key];
+    }
+   
+    /// @}
+
+
+    /// @name Missing momentum
+    /// @{
+
+    /// Get the missing momentum vector
+    /// @note Not _necessarily_ the sum over momenta of final state invisibles
     const P4& missingmom() const {
       return _pmiss;
     }
 
-    /// @brief Set the missing momentum vector
-    ///
-    /// Not _necessarily_ the sum over momenta of final state invisibles
+    /// Set the missing momentum vector
+    /// @todo Not _necessarily_ the sum over momenta of final state invisibles
     void set_missingmom(const P4& pmiss) {
       _pmiss = pmiss;
     }
@@ -421,7 +503,7 @@ namespace HEPUtils {
       return missingmom().pT();
     }
 
-    //@}
+    /// @}
 
 
   };
diff --git a/contrib/heputils/include/HEPUtils/FastJet.h b/contrib/heputils/include/HEPUtils/FastJet.h
index c9b7952416..d902feaa73 100644
--- a/contrib/heputils/include/HEPUtils/FastJet.h
+++ b/contrib/heputils/include/HEPUtils/FastJet.h
@@ -1,7 +1,7 @@
 // -*- C++ -*-
 //
-// This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils
-// Copyright (C) 2013-2018 Andy Buckley <andy.buckley@cern.ch>
+// This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/
+// Copyright (C) 2013-2022 Andy Buckley <andy.buckley@cern.ch>
 //
 // Embedding of HEPUtils code in other projects is permitted provided this
 // notice is retained and the HEPUtils namespace and include path are changed.
@@ -35,8 +35,8 @@
 namespace HEPUtils {
 
 
-  /// @name Converters between HEPUtils and FastJet momentum types
-  //@{
+  /// @defgroup fastjet_cnv Converters between HEPUtils and FastJet momentum types
+  /// @{
 
   /// @todo Enable... conditionally on FJ version?
   // /// For attaching the GenParticle provenance info to a PseudoJet
@@ -64,11 +64,11 @@ namespace HEPUtils {
     return P4::mkXYZM(p.px(), p.py(), p.pz(), (m >= 0) ? m : 0);
   }
 
-  //@}
+  /// @}
 
 
-  /// @name Jet builders
-  //@{
+  /// @defgroup fastjet_mk Jet builders
+  /// @{
 
   /// Construct pT-sorted jets using the @a alg measure with jet @a R parameter, and min pT @a ptmin (in MeV)
   inline std::vector<FJNS::PseudoJet> get_jets(const std::vector<FJNS::PseudoJet>& particles, double R, double ptmin,
@@ -79,7 +79,7 @@ namespace HEPUtils {
     return sorted_by_pt(cseq.inclusive_jets(ptmin));
   }
 
-  //@}
+  /// @}
 
 
 }
diff --git a/contrib/heputils/include/HEPUtils/Jet.h b/contrib/heputils/include/HEPUtils/Jet.h
index dfe3ae73d2..9255f848e8 100644
--- a/contrib/heputils/include/HEPUtils/Jet.h
+++ b/contrib/heputils/include/HEPUtils/Jet.h
@@ -1,13 +1,14 @@
 // -*- C++ -*-
 //
-// This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils
-// Copyright (C) 2013-2017 Andy Buckley <andy.buckley@cern.ch>
+// This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/
+// Copyright (C) 2013-2023 Andy Buckley <andy.buckley@cern.ch>
 //
 // Embedding of HEPUtils code in other projects is permitted provided this
 // notice is retained and the HEPUtils namespace and include path are changed.
 //
 #pragma once
 
+#include "HEPUtils/FastJet.h"
 #include "HEPUtils/MathUtils.h"
 #include "HEPUtils/Vectors.h"
 
@@ -15,46 +16,81 @@ namespace HEPUtils {
 
 
   /// Simple jet class, encapsulating a momentum 4-vector and with some extra b-tag info
+  ///
   /// @todo Derive from a PhysObj base class to centralise the momentum handling
   class Jet {
+  public:
+
+    /// Typedef for tag PID -> counts dictionary
+    using TagCounts = std::map<int, int>;
+
+
+  private:
 
     /// @name Storage
-    //@{
+    /// @{
+
     /// Momentum vector
     P4 _p4;
-    /// Tags
-    bool _isB, _isC, _isW, _isZ, _ish;
-    //@}
+
+    /// Tagging counts
+    TagCounts _tags;
+
+    /// Optional FastJet PJ (contains link to ClusterSeq)
+    /// @todo Use std::optional when C++17 allowed
+    FJNS::PseudoJet _pj;
+
+    /// @}
 
 
   public:
 
     /// @name Constructors
-    //@{
+    /// @{
 
-    /// Constructor for jet with tags
-    Jet(const P4& mom, bool isB=false, bool isC=false, bool isW=false, bool isZ=false, bool ish=false)
-      : _p4(mom), _isB(isB), _isC(isC), _isW(isW), _isZ(isZ), _ish(ish) {  }
+    /// Constructor for a light jet without explicit constituents
+    Jet(const P4& mom, bool isB=false, bool isC=false)
+      : _p4(mom)
+    { set_btag(int(isB)); set_ctag(int(isC)); }
 
     /// "Cartesian" constructor
     Jet(double px, double py, double pz, double E, bool isB=false, bool isC=false)
-      : _p4(px, py, pz, E), _isB(isB), _isC(isC) {  }
+      : _p4(px, py, pz, E)
+    { set_btag(int(isB)); set_ctag(int(isC)); }
 
-    //@}
+    /// "PseudoJet" constructor
+    Jet(const FJNS::PseudoJet& pj, bool isB=false, bool isC=false)
+      : _p4(mk_p4(pj)), _pj(pj)
+    { set_btag(int(isB)); set_ctag(int(isC)); }
+
+
+    /// Constructor for a light jet without explicit constituents, with a general tags map
+    Jet(const P4& mom, const TagCounts& tags)
+      : _p4(mom), _tags(tags) {  }
+
+    /// "Cartesian" constructor, with a general tags map
+    Jet(double px, double py, double pz, double E, const TagCounts& tags)
+      : _p4(px, py, pz, E), _tags(tags) {  }
+
+    /// "PseudoJet" constructor, with a general tags map
+    Jet(const FJNS::PseudoJet& pj, const TagCounts& tags)
+      : _p4(mk_p4(pj)), _tags(tags), _pj(pj) {  }
+
+    /// @}
 
 
     /// @name Implicit casts
-    //@{
+    /// @{
 
     operator const P4& () const { return mom(); }
 
     operator const P4* () const { return &mom(); }
 
-    //@}
+    /// @}
 
 
     /// @name Momentum
-    //@{
+    /// @{
 
     /// Get the 4 vector
     const P4& mom() const { return _p4; }
@@ -92,43 +128,106 @@ namespace HEPUtils {
     /// Get the squared transverse momentum
     double pT() const { return mom().pT(); }
 
-    //@}
+    /// @}
 
 
     /// @name Tagging
-    //@{
+    /// @{
+
+    /// Get the tags map (const)
+    const std::map<int,int>& tags() const { return _tags; }
+    /// Get the tags map (const)
+    std::map<int,int>& tags() { return _tags; }
+
+    /// Get the number of tags for the given PDG ID
+    int ntags(int pdgid) const { auto it = _tags.find(pdgid); return (it == _tags.end()) ? 0 : it->second; }
+    /// Get whether there is a non-zero number of tags for the given PDG ID
+    bool tagged(int pdgid) const { return (ntags(pdgid) > 0); }
+    /// Set the number of tags for the given PDG ID
+    void set_ntags(int pdgid, int ntag) { _tags[pdgid] = ntag; }
 
     /// Is this particle tagged as a b?
-    bool btag() const { return _isB; }
-    /// Set BTag value
-    void set_btag(bool isb) { _isB = isb; }
+    bool btag() const { return tagged(5); }
+    /// Set b-tag value
+    void set_btag(int ntag=1) { set_ntags(5, ntag); }
 
     /// Is this particle tagged as a c?
-    /// @note Can be simultaneously btag()'d -- analyses should probably only use if fallback from b-tag.
-    bool ctag() const { return /* !btag() && */ _isC; }
-    /// Set CTag value
-    void set_ctag(bool isc) { _isC = isc; }
+    ///
+    /// @note Can be simultaneously b-tagged unless the optional bool is set true
+    bool ctag(bool not_if_b=false) const {
+      if (not_if_b) return tagged(5) ? false : tagged(4);
+      return tagged(4);
+    }
+    /// Set c-tag value
+    void set_ctag(int ntag=1) { set_ntags(4, ntag); }
 
-    /// Is this particle tagged as a W?
-    bool Wtag()  const { return _isW; }
-    
-    /// Is this particle tagged as a Z?
-    bool Ztag()  const { return _isZ; }
+    /// @}
 
-    /// Is this particle tagged as a SM Higgs?
-    bool htag()  const { return _ish; }
 
-    /// @todo Generalize for charm tags, tau tags, multiple tags of a single type?
+    /// @name FastJet information
+    /// @{
 
-    //@}
+    /// Get the contained PseudoJet object (const)
+    const FJNS::PseudoJet& pseudojet() const { return _pj; }
+
+    /// Get the contained PseudoJet object
+    FJNS::PseudoJet& pseudojet() { return _pj; }
+
+    /// Set the contained PseudoJet object
+    void set_pseudojet(const FJNS::PseudoJet& pj) { _pj = pj; }
+
+    /// @brief Access the ClusterSequence object if possible (can be null)
+    ///
+    /// Optional template arg can be used to cast to a specific derived CS type if wanted.
+    template <typename CS=FJNS::ClusterSequence>
+    typename std::shared_ptr<const CS> clusterseq() const {  
+      return std::shared_ptr<const CS>(clusterseq_raw<CS>());
+    }
+    
+    /// @brief Access the ClusterSequence object if possible (can be null)
+    ///
+    /// Optional template arg can be used to cast to a specific derived CS type if wanted.
+    template <typename CS=FJNS::ClusterSequence>
+    const CS* clusterseq_raw() const {
+      return dynamic_cast<const CS*>(_pj.associated_cs());
+    }
+
+    /// @}
 
   };
 
 
-  /// Function/functor for container<const Jet*> sorting (cf. std::less)
-  inline bool _cmpPtDesc(const Jet* a, const Jet* b) {
-    return a->pT2() >= b->pT2();
+  /// @defgroup jet_const Jet constness conversion
+  /// @{
+
+  /// Convenience Jet cast to const
+  inline const Jet* mkconst(Jet* jet) {
+    return const_cast<const Jet*>(jet);
+  }
+
+  /// Convenience Jet cast to non-const
+  inline Jet* mkunconst(const Jet* cjet) {
+    return const_cast<Jet*>(cjet);
+  }
+
+  /// Get a reference to a vector of Jets, with each member const
+  inline std::vector<const Jet*>& mkconst(const std::vector<Jet*>& jets) {
+    return * (std::vector<const Jet*>*)(void*) (&jets);
+  }
+
+  /// Get a reference to a vector of Jets, with each member non-const
+  inline std::vector<Jet*>& mkunconst(const std::vector<const Jet*>& cjets) {
+    return * (std::vector<Jet*>*) (void*) (&cjets);
   }
 
+  /// @}
+
+
+
+  // /// Function/functor for container<const Jet*> sorting (cf. std::less)
+  // inline bool _cmpPtDescJet(const Jet* a, const Jet* b) {
+  //   return a->pT2() >= b->pT2();
+  // }
+
 
 }
diff --git a/contrib/heputils/include/HEPUtils/MathUtils.h b/contrib/heputils/include/HEPUtils/MathUtils.h
index 1a10e76745..850285ff46 100644
--- a/contrib/heputils/include/HEPUtils/MathUtils.h
+++ b/contrib/heputils/include/HEPUtils/MathUtils.h
@@ -1,7 +1,7 @@
 // -*- C++ -*-
 //
-// This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils
-// Copyright (C) 2013-2018 Andy Buckley <andy.buckley@cern.ch>
+// This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/
+// Copyright (C) 2013-2022 Andy Buckley <andy.buckley@cern.ch>
 //
 // Embedding of HEPUtils code in other projects is permitted provided this
 // notice is retained and the HEPUtils namespace and include path are changed.
@@ -25,8 +25,11 @@
 namespace HEPUtils {
 
 
-  /// @name Numerical helpers
-  //@{
+  /// @defgroup mathutils Mathematical utilities
+  /// @{
+
+  /// @defgroup mathutils_numhelp Numerical helpers
+  /// @{
 
   /// Convenience function for getting the sign of a number
   template <typename N1>
@@ -58,11 +61,11 @@ namespace HEPUtils {
     return sqrt(a*a + b*b + c*c);
   }
 
-  //@}
+  /// @}
 
 
-  /// @name Range helpers
-  //@{
+  /// @defgroup mathutils_rangehelp Range helpers
+  /// @{
 
   /// @brief Boolean function to determine if @a value is within the given range
   ///
@@ -128,11 +131,11 @@ namespace HEPUtils {
     return rtn;
   }
 
-  //@}
+  /// @}
 
 
-  /// @name Physics maths utils
-  //@{
+  /// @defgroup mathutils_physics Physics maths utils
+  /// @{
 
   inline double deltaphi(double a, double b) {
     double rtn = a - b;
@@ -147,17 +150,17 @@ namespace HEPUtils {
     return rtn;
   }
 
-  //@}
+  /// @}
 
 
-  /// @name Random numbers and sampling
-  //@{
+  /// @defgroup mathutils_random Random numbers and sampling
+  /// @{
 
   inline double rand01() {
     return rand() / (double)RAND_MAX;
   }
 
-  //@}
+  /// @}
 
 
 }
diff --git a/contrib/heputils/include/HEPUtils/Particle.h b/contrib/heputils/include/HEPUtils/Particle.h
index fd7c412a67..dc70972b56 100644
--- a/contrib/heputils/include/HEPUtils/Particle.h
+++ b/contrib/heputils/include/HEPUtils/Particle.h
@@ -1,7 +1,7 @@
 // -*- C++ -*-
 //
-// This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils
-// Copyright (C) 2013-2018 Andy Buckley <andy.buckley@cern.ch>
+// This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/
+// Copyright (C) 2013-2022 Andy Buckley <andy.buckley@cern.ch>
 //
 // Embedding of HEPUtils code in other projects is permitted provided this
 // notice is retained and the HEPUtils namespace and include path are changed.
@@ -14,81 +14,77 @@
 namespace HEPUtils {
 
 
-  /// Simple particle class, encapsulating a momentum 4-vector and adding some extra ID info
-  /// @todo Derive from a PhysObj base class to centralise the momentum handling
-  /// @todo Provide cast operators to P4 and P4*
+  /// @brief Simple particle class, encapsulating a momentum 4-vector and adding some extra ID info
+  ///
+  /// @todo Derive from a PhysObj base class to centralise the momentum handling?
   class Particle {
   private:
 
     /// @name Storage
-    //@{
+    /// @{
+
     /// Momentum vector
     P4 _p4;
+
     /// PDG ID code
-    int _pdgId;
+    int _pid;
+
     /// Promptness flag
     bool _prompt;
-    // /// Isolation value
-    // double _isol4;
-    //@}
+
+    /// @}
 
 
   public:
 
     /// @name Constructors
-    //@{
+    /// @{
 
     /// Default constructor
     Particle()
-      : _pdgId(0), _prompt(false) { } //, _isol4(0.0) {  }
+      : _pid(0), _prompt(false) { }
 
     /// "Cartesian" constructor
     Particle(double px, double py, double pz, double E, int pdgid)
-      : _p4(px, py, pz, E), _pdgId(pdgid), _prompt(false) { } //, _isol4(0.0) {  }
-
-    /// "Cartesian" constructor for massless particles - or close enough
-    /// @todo AB: WTF?
-    // Particle(double px, double py, double pz, int pdgid)
-    //   : _p4(px, py, pz), _pdgId(pdgid), _prompt(false), _isol4(0.0) {  }
+      : _p4(px, py, pz, E), _pid(pdgid), _prompt(false) { }
 
     /// 4-mom + PDG ID constructor
     Particle(const P4& mom, int pdgid)
-      : _p4(mom), _pdgId(pdgid), _prompt(false) { } //, _isol4(0.0) {  }
+      : _p4(mom), _pid(pdgid), _prompt(false) { }
 
     /// Copy constructor
     Particle(const Particle& p)
-      : _p4(p.mom()), _pdgId(p.pid()), _prompt(p.is_prompt()) { } //, _isol4(p.isol()) {  }
+      : _p4(p.mom()), _pid(p.pid()), _prompt(p.is_prompt()) { }
 
     /// Copy constructor from a pointer
     Particle(const Particle* p)
-      : _p4(p->mom()), _pdgId(p->pid()), _prompt(p->is_prompt()) { } //, _isol4(p->isol()) {  }
+      : _p4(p->mom()), _pid(p->pid()), _prompt(p->is_prompt()) { }
 
     /// Copy assignment operator
-    Particle& operator=(const Particle& p) {
+    Particle& operator = (const Particle& p) {
       _p4 = p.mom();
-      _pdgId = p.pid();
+      _pid = p.pid();
       _prompt = p.is_prompt();
-      //_isol4 = p.isol();
       return *this;
     }
 
-    //@}
+    /// @}
 
 
     /// @name Implicit casts
-    //@{
+    /// @{
 
     operator const P4& () const { return mom(); }
 
     operator const P4* () const { return &mom(); }
 
-    //@}
+    /// @}
 
 
     /// @name Momentum
     ///
     /// Access to the P4 object, plus convenience mapping of a few popular properties
-    //@{
+    /// @{
 
     /// Get the 4 vector
     const P4& mom() const { return _p4; }
@@ -126,44 +122,88 @@ namespace HEPUtils {
     /// Get the squared transverse momentum
     double pT() const { return mom().pT(); }
 
-    //@}
+    /// @}
 
 
     /// @name Promptness
-    //@{
+    /// @{
 
     /// Is this particle connected to the hard process or from a hadron/tau decay?
     bool is_prompt() const { return _prompt; }
     /// Set promptness
     void set_prompt(bool isprompt=true) { _prompt = isprompt; }
 
-    //@}
+    /// @}
 
 
     /// @name Particle ID
-    //@{
+    /// @{
 
     /// Get PDG ID code
-    int pid() const { return _pdgId; }
+    int pid() const { return _pid; }
     /// Get abs PDG ID code
-    int abspid() const { return abs(_pdgId); }
+    int abspid() const { return abs(_pid); }
     /// Set PDG ID code
-    void set_pid(int pid) { _pdgId = pid; }
+    void set_pid(int pdgid) { _pid = pdgid; }
+
+    /// Is this particle usually visible in a detector?
+    bool is_visible() {
+      if (abspid() == 12 || abspid() == 14 || abspid() == 16) return false;
+      if (pid() == 1000022 || pid() == 1000039) return false;
+      if (in_range(abspid(), 50, 60)) return false; //< abspid zealousness since some -ve DM PIDs seen
+      return true;
+    }
 
-    //@}
+    /// Is this particle usually invisible in a detector?
+    bool is_invisible() {
+      return !is_visible();
+    }
+    /// @}
 
 
     // /// @name Isolation of particle
-    // //@{
+    // /// @{
 
     // /// Get isolation
     // double isol() const { return _isol4;}
     // void set_isol(double isol) { _isol4 = isol;}
 
-    // //@}
+    // /// @}
 
 
   };
 
 
+  /// @defgroup particle_const Particle constness conversion
+  /// @{
+
+  /// Convenience Particle cast to const
+  inline const Particle* mkconst(Particle* particle) {
+    return const_cast<const Particle*>(particle);
+  }
+
+  /// Convenience Particle cast to non-const
+  inline Particle* mkunconst(const Particle* cparticle) {
+    return const_cast<Particle*>(cparticle);
+  }
+
+  /// Get a reference to a vector of Particles, with each member const
+  inline std::vector<const Particle*>& mkconst(const std::vector<Particle*>& particles) {
+    return * (std::vector<const Particle*>*)(void*) (&particles);
+  }
+
+  /// Get a reference to a vector of Particles, with each member non-const
+  inline std::vector<Particle*>& mkunconst(const std::vector<const Particle*>& cparticles) {
+    return * (std::vector<Particle*>*) (void*) (&cparticles);
+  }
+
+  /// @}
+
+
+  // /// Function/functor for container<const Particle*> sorting (cf. std::less)
+  // inline bool _cmpPtDescParticle(const Particle* a, const Particle* b) {
+  //   return a->pT2() >= b->pT2();
+  // }
+
+
 }
diff --git a/contrib/heputils/include/HEPUtils/PhaseSpace.h b/contrib/heputils/include/HEPUtils/PhaseSpace.h
index c68888a43f..a1501facf2 100644
--- a/contrib/heputils/include/HEPUtils/PhaseSpace.h
+++ b/contrib/heputils/include/HEPUtils/PhaseSpace.h
@@ -1,7 +1,7 @@
 // -*- C++ -*-
 //
-// This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils
-// Copyright (C) 2013-2018 Andy Buckley <andy.buckley@cern.ch>
+// This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/
+// Copyright (C) 2013-2022 Andy Buckley <andy.buckley@cern.ch>
 //
 // Embedding of HEPUtils code in other projects is permitted provided this
 // notice is retained and the HEPUtils namespace and include path are changed.
@@ -20,8 +20,12 @@
 namespace HEPUtils {
 
 
+  /// @defgroup phasespace Phase-space helpers
+  /// @{
+
+
   /// Integral of acos(x)
-  constexpr double iacos(double x) { return sqrt(1 - x*x) - x*acos(x); }
+  constexpr double iacos(double x) { return x*acos(x) - sqrt(1 - x*x); }
 
   /// Inverse secant
   constexpr double asec(double x) { return acos(1/x); }
@@ -74,5 +78,7 @@ namespace HEPUtils {
     return psPhiIntegral;
   }
 
+  /// @}
+
 
 }
diff --git a/contrib/heputils/include/HEPUtils/Utils.h b/contrib/heputils/include/HEPUtils/Utils.h
index b40f5497b2..c6f762e466 100644
--- a/contrib/heputils/include/HEPUtils/Utils.h
+++ b/contrib/heputils/include/HEPUtils/Utils.h
@@ -1,28 +1,26 @@
 // -*- C++ -*-
 //
-// This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils
-// Copyright (C) 2013-2018 Andy Buckley <andy.buckley@cern.ch>
+// This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/
+// Copyright (C) 2013-2022 Andy Buckley <andy.buckley@cern.ch>
 //
 // Embedding of HEPUtils code in other projects is permitted provided this
 // notice is retained and the HEPUtils namespace and include path are changed.
 //
 #pragma once
 
-#if defined(__cplusplus) && __cplusplus < 201103L
-// #define XSTR(x) STR(x)
-// #define STR(x) #x
-// #pragma message STR(__cplusplus)
-#error "This library needs at least a C++11 compliant compiler"
+#if  __cplusplus < 201103L
+#pragma message "This library needs at least a C++11 compliant compiler"
 #endif
 
+
 /// @file Utility functions
 /// @author Andy Buckley <andy.buckley@cern.ch>
 
 namespace HEPUtils {
 
 
-  /// @name Container utils
-  //@{
+  /// @defgroup utils_container Container utils
+  /// @{
 
   /// Return true if f(x) is true for any x in container c, otherwise false.
   template <typename CONTAINER, typename FN>
@@ -52,7 +50,7 @@ namespace HEPUtils {
   }
 
 
-  //@}
+  /// @}
 
 
 }
diff --git a/contrib/heputils/include/HEPUtils/Vectors.h b/contrib/heputils/include/HEPUtils/Vectors.h
index 77f32701f2..f632aea81f 100644
--- a/contrib/heputils/include/HEPUtils/Vectors.h
+++ b/contrib/heputils/include/HEPUtils/Vectors.h
@@ -1,7 +1,7 @@
 // -*- C++ -*-
 //
-// This file is part of HEPUtils -- https://bitbucket.org/andybuckley/heputils
-// Copyright (C) 2013-2018 Andy Buckley <andy.buckley@cern.ch>
+// This file is part of HEPUtils -- https://gitlab.com/hepcedar/heputils/
+// Copyright (C) 2013-2023 Andy Buckley <andy.buckley@cern.ch>
 //
 // Embedding of HEPUtils code in other projects is permitted provided this
 // notice is retained and the HEPUtils namespace and include path are changed.
@@ -367,16 +367,20 @@ namespace HEPUtils {
     /// Get the transverse momentum (same as rho)
     double pT() const { return rho(); }
 
-    /// Get the spatial phi
-    double phi() const { if (rho2() == 0) return 0; else return atan2(py(),px()); }
-    /// Get the spatial theta
-    double theta() const { if (p2() == 0) return 0; else if (pz() == 0) return M_PI; else return atan2(rho(),pz()); }
-    /// Get the spatial vector pseudorapidity
-    double eta() const { return -log(tan( 0.5 * theta() )); } //< Optimise with a trig reln on tan(x/2) to avoid tan(atan(..)/2)?
-    /// Get the spatial vector absolute pseudorapidity
+    /// Get the spatial phi (in the range -pi .. pi)
+    double phi() const { if (rho2() == 0) return 0; else return atan2(py(), px()); }
+    /// Get the spatial phi (in the range 0 .. 2pi)
+    double phi_02pi() const { if (rho2() == 0) return 0; else return phi() + M_PI; }
+
+    /// Get the spatial theta (in the range 0 .. pi)
+    double theta() const { if (p2() == 0) return 0; else
+	if (pz() == 0) return M_PI/2; else return atan2(rho(), pz()); } //< atan2(+ve, z) is +ve
+    /// Get the spatial-vector pseudorapidity
+    double eta() const { return std::copysign(log((p() + fabs(pz())) / pT()), pz()); }
+    /// Get the spatial-vector absolute pseudorapidity
     double abseta() const { return fabs(eta()); }
     /// Get the 4-momentum rapidity
-    double rap() const { return 0.5 * (E() + pz()) / (E() - pz()); }
+    double rap() const { return 0.5 * log((E() + pz()) / (E() - pz())); }
     /// Get the 4-momentum absolute rapidity
     double absrap() const { return fabs(rap()); }
 
@@ -505,4 +509,17 @@ namespace HEPUtils {
   //@}
 
 
+  /// Function/functor for container<const P4> sorting (cf. std::less)
+  template <typename T>
+  inline bool _cmpPtDesc(const T& a, const T& b) {
+    return a.pT2() >= b.pT2();
+  }
+
+  /// Function/functor for container<const P4*> sorting (cf. std::less)
+  template <typename T>
+  inline bool _cmpPtDescPtr(const T* a, const T* b) {
+     return _cmpPtDesc(*a, *b);
+  }
+
+
 }
diff --git a/gum/gum.py b/gum/gum.py
index 563a92de8c..2c4b0a367e 100644
--- a/gum/gum.py
+++ b/gum/gum.py
@@ -813,10 +813,10 @@
                            num+10, reset_contents)
 
             # write all invisible particles in the model to Event header
-            num = find_string("heputils/include/HEPUtils/Event.h", "contrib",
-                                  "        _cinvisibles.push_back(p);")[1]
-            amend_file("heputils/include/HEPUtils/Event.h", "contrib", get_invisibles(gum.invisibles_pdg),
-                           num+1, reset_contents)
+            num = find_string("heputils/include/HEPUtils/Particle.h", "contrib",
+                                  "      if (pid() == 1000022 || pid() == 1000039) return false;")[1]
+            amend_file("heputils/include/HEPUtils/Particle.h", "contrib", get_invisibles(gum.invisibles_pdg),
+                           num, reset_contents)
 
         # HiggsBounds interface
         if output_opts.spheno:
diff --git a/gum/src/particledb.py b/gum/src/particledb.py
index 4f1cd30ffc..ade2ef16ae 100644
--- a/gum/src/particledb.py
+++ b/gum/src/particledb.py
@@ -394,8 +394,8 @@ def get_higgses(partlist):
     
 def get_invisibles(invisible_pdgs):
     """
-    Get the PDG codes of invisibles particles to write to 'contrib/heputils/include/HEPUtils/Event.h'
-    Will not write for photons, leptons or existing invisibles in Event.h
+    Get the PDG codes of invisibles particles to write to 'contrib/heputils/include/HEPUtils/Particle.h'
+    Will not write for photons, leptons or existing invisibles
     """
     
     # Create list of existing particle
@@ -407,15 +407,7 @@ def get_invisibles(invisible_pdgs):
     to_write = ""
 
     if (len(new_invisibles) != 0):
-        pdg_string = "p->pid() == {0}".format(new_invisibles[0])
         for i,pdg in enumerate(new_invisibles):
-            if i != 0:
-                pdg_string = pdg_string + "|| p->pid() == {0} ".format(pdg)
-    
-        to_write = ("      else if (" + pdg_string + ") \n"
-        "      { \n"
-        "        _invisibles.push_back(p); \n"
-        "        _cinvisibles.push_back(p); \n"
-        "      } \n")
+            to_write = to_write + ("      if (pid() == {0}) return false;\n".format(pdg))
     
     return to_write
diff --git a/yaml_files/CMSSM.yaml b/yaml_files/CMSSM.yaml
index 46f2fc16f6..9f34496e15 100644
--- a/yaml_files/CMSSM.yaml
+++ b/yaml_files/CMSSM.yaml
@@ -708,7 +708,13 @@ Rules:
         LHC_13TeV:
           # 0.028 fb corresponds to ~1 expected event at L = 36 fb^-1.
           xsec_veto: 0.028
-          antiktR: 0.4
+          jet_collections:
+            antikt_R04:
+              algorithm: antikt
+              R: 0.4
+              recombination_scheme: E_scheme
+              strategy: Best
+          jet_collection_taus: antikt_R04
           partonOnly: false
           pythia_settings:
             - Print:quiet = on
diff --git a/yaml_files/ColliderBit_CMSSM.yaml b/yaml_files/ColliderBit_CMSSM.yaml
index 313f2a1fc1..ef6835e9c0 100644
--- a/yaml_files/ColliderBit_CMSSM.yaml
+++ b/yaml_files/ColliderBit_CMSSM.yaml
@@ -310,7 +310,13 @@ Rules:
       LHC_8TeV:
         xsec_veto: 0.00
         partonOnly: false
-        antiktR: 0.4
+        jet_collections:
+          antikt_R04:
+            algorithm: antikt
+            R: 0.4
+            recombination_scheme: E_scheme
+            strategy: Best
+        jet_collection_taus: antikt_R04
         pythia_settings:
           - Beams:eCM = 8000
           - Next:numberShowProcess = 1
@@ -329,7 +335,18 @@ Rules:
         # 0.028 fb corresponds to ~1 expected event at L = 36 fb^-1.
         xsec_veto: 0.028
         partonOnly: false
-        antiktR: 0.4
+        jet_collections:
+          antikt_R04:
+            algorithm: antikt  # antikt, cambridge, kt, genkt, cambridge_for_passive
+            R: 0.4
+            recombination_scheme: E_scheme   # E_scheme, pt_scheme, pt2_scheme 
+            strategy: Best     # Best, NlnN
+          # antikt_R08:
+          #   algorithm: antikt
+          #   R: 0.8
+          #   recombination_scheme: E_scheme
+          #   strategy: Best
+        jet_collection_taus: antikt_R04
         pythia_settings:
           - Beams:eCM = 13000
           - Next:numberShowProcess = 1
diff --git a/yaml_files/MSSM7.yaml b/yaml_files/MSSM7.yaml
index 2e629db6b1..a4b32aa393 100644
--- a/yaml_files/MSSM7.yaml
+++ b/yaml_files/MSSM7.yaml
@@ -616,7 +616,13 @@ Rules:
         # 0.028 fb corresponds to ~1 expected event at L = 36 fb^-1.
         xsec_veto: 0.028
         partonOnly: false
-        antiktR: 0.4
+        jet_collections:
+          antikt_R04:
+            algorithm: antikt
+            R: 0.4
+            recombination_scheme: E_scheme
+            strategy: Best
+        jet_collection_taus: antikt_R04
         pythia_settings:
           - Print:quiet = on
           - Next:numberCount = 0
diff --git a/yaml_files/MSSM9.yaml b/yaml_files/MSSM9.yaml
index fe6de4ccfd..d18401677d 100644
--- a/yaml_files/MSSM9.yaml
+++ b/yaml_files/MSSM9.yaml
@@ -624,7 +624,13 @@ Rules:
         # 0.028 fb corresponds to ~1 expected event at L = 36 fb^-1.
         xsec_veto: 0.028
         partonOnly: false
-        antiktR: 0.4
+        jet_collections:
+          antikt_R04:
+            algorithm: antikt
+            R: 0.4
+            recombination_scheme: E_scheme
+            strategy: Best
+        jet_collection_taus: antikt_R04
         pythia_settings:
           - Print:quiet = on
           - Next:numberCount = 0
diff --git a/yaml_files/NUHM1.yaml b/yaml_files/NUHM1.yaml
index 8a160a4ddb..1693527f74 100644
--- a/yaml_files/NUHM1.yaml
+++ b/yaml_files/NUHM1.yaml
@@ -608,7 +608,13 @@ Rules:
         # 0.028 fb corresponds to ~1 expected event at L = 36 fb^-1.
         xsec_veto: 0.028
         partonOnly: false
-        antiktR: 0.4
+        jet_collections:
+          antikt_R04:
+            algorithm: antikt
+            R: 0.4
+            recombination_scheme: E_scheme
+            strategy: Best
+        jet_collection_taus: antikt_R04
         pythia_settings:
           - Print:quiet = on
           - Next:numberCount = 0
diff --git a/yaml_files/NUHM2.yaml b/yaml_files/NUHM2.yaml
index f5a53e5440..2d3af567c6 100644
--- a/yaml_files/NUHM2.yaml
+++ b/yaml_files/NUHM2.yaml
@@ -611,7 +611,13 @@ Rules:
         # 0.028 fb corresponds to ~1 expected event at L = 36 fb^-1.
         xsec_veto: 0.028
         partonOnly: false
-        antiktR: 0.4
+        jet_collections:
+          antikt_R04:
+            algorithm: antikt
+            R: 0.4
+            recombination_scheme: E_scheme
+            strategy: Best
+        jet_collection_taus: antikt_R04
         pythia_settings:
           - Print:quiet = on
           - Next:numberCount = 0