From 650549f46bb18c6b139d9671273bc5c2eeef831d Mon Sep 17 00:00:00 2001
From: Jon Burr <jon.burr@cern.ch>
Date: Tue, 8 Dec 2020 13:34:24 +0100
Subject: [PATCH] Add iterators over combinations and first implementation of
 matching tool using them

---
 .../TriggerMatchingTool/Root/MatchingTool.cxx |   1 -
 .../Root/R3MatchingTool.cxx                   | 174 ++++++++++++++----
 .../TriggerMatchingTool/R3MatchingTool.h      |  10 +
 .../TrigCompositeUtils/CMakeLists.txt         |   3 +-
 .../TrigCompositeUtils/Root/Combinations.cxx  |  77 ++++++++
 .../TrigCompositeUtils/Root/IPartCombItr.cxx  | 173 +++++++++++++++++
 .../TrigCompositeUtils/Root/KFromNItr.cxx     | 109 +++++++++++
 .../TrigCompositeUtils/Combinations.h         |  50 +++++
 .../TrigCompositeUtils/IPartCombItr.h         | 139 ++++++++++++++
 .../TrigCompositeUtils/KFromNItr.h            |  77 ++++++++
 10 files changed, 777 insertions(+), 36 deletions(-)
 create mode 100644 Trigger/TrigSteer/TrigCompositeUtils/Root/Combinations.cxx
 create mode 100644 Trigger/TrigSteer/TrigCompositeUtils/Root/IPartCombItr.cxx
 create mode 100644 Trigger/TrigSteer/TrigCompositeUtils/Root/KFromNItr.cxx
 create mode 100644 Trigger/TrigSteer/TrigCompositeUtils/TrigCompositeUtils/Combinations.h
 create mode 100644 Trigger/TrigSteer/TrigCompositeUtils/TrigCompositeUtils/IPartCombItr.h
 create mode 100644 Trigger/TrigSteer/TrigCompositeUtils/TrigCompositeUtils/KFromNItr.h

diff --git a/Trigger/TrigAnalysis/TriggerMatchingTool/Root/MatchingTool.cxx b/Trigger/TrigAnalysis/TriggerMatchingTool/Root/MatchingTool.cxx
index 3ebdf9ee50ec..c96e4e612811 100644
--- a/Trigger/TrigAnalysis/TriggerMatchingTool/Root/MatchingTool.cxx
+++ b/Trigger/TrigAnalysis/TriggerMatchingTool/Root/MatchingTool.cxx
@@ -17,7 +17,6 @@ MatchingTool::MatchingTool(const std::string& name) :
   m_matchingThreshold(0)
 {
   declareProperty( "TrigDecisionTool", m_trigDecTool);
-
 }
 
 #ifndef XAOD_STANDALONE
diff --git a/Trigger/TrigAnalysis/TriggerMatchingTool/Root/R3MatchingTool.cxx b/Trigger/TrigAnalysis/TriggerMatchingTool/Root/R3MatchingTool.cxx
index 6ef34c6e8bbb..a7a85ba1d6a6 100644
--- a/Trigger/TrigAnalysis/TriggerMatchingTool/Root/R3MatchingTool.cxx
+++ b/Trigger/TrigAnalysis/TriggerMatchingTool/Root/R3MatchingTool.cxx
@@ -4,6 +4,10 @@
 
 #include "TriggerMatchingTool/R3MatchingTool.h"
 #include "xAODBase/IParticleContainer.h"
+#include "TrigCompositeUtils/Combinations.h"
+#include "xAODEgamma/Egamma.h"
+#include <numeric>
+#include <algorithm>
 
 // Anonymous namespace for helper functions
 namespace
@@ -11,14 +15,15 @@ namespace
   // Helper function to efficiently decide if two objects are within a drThreshold of each other
   bool fastDR(float eta1, float phi1, float eta2, float phi2, float drThreshold)
   {
-    float dEta = eta1 - eta2;
+    float dEta = std::abs(eta1 - eta2);
     if (dEta > drThreshold)
       return false;
-    float dPhi = phi1 - phi2;
+    float dPhi = std::abs(TVector2::Phi_mpi_pi(phi1 - phi2));
     if (dPhi > drThreshold)
       return false;
     return dEta * dEta + dPhi * dPhi < drThreshold * drThreshold;
   }
+
 } // namespace
 
 namespace Trig
@@ -34,6 +39,7 @@ namespace Trig
   StatusCode R3MatchingTool::initialize()
   {
     ATH_CHECK(m_trigDecTool.retrieve());
+    m_trigDecTool->ExperimentalAndExpertMethods()->enable();
     return StatusCode::SUCCESS;
   }
 
@@ -48,47 +54,108 @@ namespace Trig
       return true;
     // Make the LinkInfo type less verbose
     using IPartLinkInfo_t = TrigCompositeUtils::LinkInfo<xAOD::IParticleContainer>;
+    using VecLinkInfo_t = std::vector<IPartLinkInfo_t>;
     // TODO - detect if we're looking at run 3 data.
     // If we are, then setting rerun to true should give a warning as it no
     // longer makes sense for run 3
-    if (recoObjects.size() != 1)
-    {
-      ATH_MSG_WARNING("Matching of multiple objects is not yet supported!");
-      return false;
-    }
-    // This has to assume that the last IParticle feature in each chain is the
-    // correct one - I don't *think* this will cause any issues.
-    // Unlike with TEs, each chain should build in nodes that it actually uses.
-    // This means that we don't need anything clever for the etcut triggers any more (for example)
-    std::vector<IPartLinkInfo_t> features = m_trigDecTool->features<xAOD::IParticleContainer>(
-        chain,
-        rerun ? TrigDefs::Physics | TrigDefs::allowResurrectedDecision : TrigDefs::Physics);
-    ATH_MSG_DEBUG("Found " << features.size() << " features for chain group " << chain);
-    // TODO:
-    //  Right now we are only looking at single object chains, so we only need
-    //  to find a single match for one reco-object.
-    //  Longer term we need to deal with combinations so this gets harder
-    const xAOD::IParticle &recoObject = *recoObjects.at(0);
-    for (const IPartLinkInfo_t &info : features)
+    unsigned int condition = rerun ? TrigDefs::Physics | TrigDefs::allowResurrectedDecision : TrigDefs::Physics;
+
+    // In what follows, the same comparisons between reco and trigger objects will done
+    // fairly frequently. As these include DR checks we want to minimise how often we do these
+    // Therefore we keep track of any comparisons that we've already done
+    // There is one map per input reco object, and then each map entry is the keyed on information
+    // extracted from the element link
+    std::vector<std::map<std::pair<uint32_t, uint32_t>, bool>> cachedComparisons(recoObjects.size());
+
+    // Note that the user can supply a regex pattern which matches multiple chains
+    // We should return true if any individual chain matches
+    const Trig::ChainGroup *chainGroup = m_trigDecTool->getChainGroup(chain);
+    for (const std::string &chainName : chainGroup->getListOfTriggers())
     {
-      if (!info.link.isValid())
+      // Now we have to build up combinations
+      // TODO - right now we use a filter that passes everything through.
+      // This will probably need to be fixed to something else later - at least the unique RoI filter
+      TrigCompositeUtils::Combinations combinations(TrigCompositeUtils::FilterType::All);
+      const TrigConf::HLTChain *chainInfo = m_trigDecTool->ExperimentalAndExpertMethods()->getChainConfigurationDetails(chainName);
+      std::vector<std::size_t> multiplicities = chainInfo->leg_multiplicities();
+      combinations.reserve(multiplicities.size());
+      for (std::size_t legIdx = 0; legIdx < multiplicities.size(); ++legIdx)
       {
-        // This could result in false negatives...
-        ATH_MSG_WARNING("Invalid link to trigger feature for chain " << chain);
+        HLT::Identifier legID = TrigCompositeUtils::createLegName(chainName, legIdx);
+        combinations.addLeg(
+            multiplicities.at(legIdx),
+            m_trigDecTool->features<xAOD::IParticleContainer>(legID.name(), condition));
+      }
+      // Warn once per call if one of the chain groups is too small to match anything
+      if (combinations.size() < recoObjects.size())
+      {
+        ATH_MSG_WARNING(
+            "Chain " << chainName << " (matching pattern " << chain << ") has too few objects ("
+                     << combinations.size() << ") to match the number of provided reco objects (" << recoObjects.size() << ")");
         continue;
       }
-      const xAOD::IParticle &trigObject = **info.link;
-      if (fastDR(
-              recoObject.eta(),
-              recoObject.phi(),
-              trigObject.eta(),
-              trigObject.phi(),
-              matchThreshold))
-        // Once we find one match, that is enough
-        return true;
+      // Now we iterate through the available combinations
+      for (const VecLinkInfo_t &combination : combinations)
+      {
+        // Prepare the index vector
+        std::vector<std::size_t> onlineIndices(combination.size());
+        std::iota(onlineIndices.begin(), onlineIndices.end(), 0);
+        do
+        {
+          bool match = true;
+          for (std::size_t recoIdx = 0; recoIdx < recoObjects.size(); ++recoIdx)
+            if (!matchObjects(recoObjects[recoIdx], combination[onlineIndices[recoIdx]].link, cachedComparisons[recoIdx], matchThreshold))
+            {
+              match = false;
+              break;
+            }
+          if (match)
+            return true;
+        } while (std::next_permutation(onlineIndices.begin(), onlineIndices.end()));
+      }
     }
-    // If we get here then there was no good match
+
+    // If we reach here we've tried all combinations from all chains in the group and none of them matched
     return false;
+
+    // if (recoObjects.size() != 1)
+    // {
+    //   ATH_MSG_WARNING("Matching of multiple objects is not yet supported!");
+    //   return false;
+    // }
+    // // This has to assume that the last IParticle feature in each chain is the
+    // // correct one - I don't *think* this will cause any issues.
+    // // Unlike with TEs, each chain should build in nodes that it actually uses.
+    // // This means that we don't need anything clever for the etcut triggers any more (for example)
+    // VecLinkInfo features = m_trigDecTool->features<xAOD::IParticleContainer>(
+    //     chain,
+    //     rerun ? TrigDefs::Physics | TrigDefs::allowResurrectedDecision : TrigDefs::Physics);
+    // ATH_MSG_DEBUG("Found " << features.size() << " features for chain group " << chain);
+    // // TODO:
+    // //  Right now we are only looking at single object chains, so we only need
+    // //  to find a single match for one reco-object.
+    // //  Longer term we need to deal with combinations so this gets harder
+    // const xAOD::IParticle &recoObject = *recoObjects.at(0);
+    // for (const IPartLinkInfo_t &info : features)
+    // {
+    //   if (!info.link.isValid())
+    //   {
+    //     // This could result in false negatives...
+    //     ATH_MSG_WARNING("Invalid link to trigger feature for chain " << chain);
+    //     continue;
+    //   }
+    //   const xAOD::IParticle &trigObject = **info.link;
+    //   if (fastDR(
+    //           recoObject.eta(),
+    //           recoObject.phi(),
+    //           trigObject.eta(),
+    //           trigObject.phi(),
+    //           matchThreshold))
+    //     // Once we find one match, that is enough
+    //     return true;
+    // }
+    // // If we get here then there was no good match
+    // return false;
   }
 
   bool R3MatchingTool::match(
@@ -101,4 +168,43 @@ namespace Trig
     return match(tmpVec, chain, matchThreshold, rerun);
   }
 
+  bool R3MatchingTool::matchObjects(
+      const xAOD::IParticle *reco,
+      const ElementLink<xAOD::IParticleContainer> &onlineLink,
+      std::map<std::pair<uint32_t, uint32_t>, bool> &cache,
+      double drThreshold) const
+  {
+    if (!onlineLink.isValid())
+    {
+      ATH_MSG_WARNING("Invalid element link!");
+      return false;
+    }
+    std::pair<uint32_t, uint32_t> linkIndices(onlineLink.persKey(), onlineLink.persIndex());
+    auto cacheItr = cache.find(linkIndices);
+    if (cacheItr == cache.end())
+    {
+      const xAOD::IParticle *online = *onlineLink;
+      bool match = online->type() == reco->type();
+      if (online->type() == xAOD::Type::CaloCluster && (reco->type() == xAOD::Type::Electron || reco->type() == xAOD::Type::Photon))
+      {
+        // Calo cluster is a special case - some of the egamma chains can return these (the etcut chains)
+        // In these cases we need to match this against the caloCluster object contained in electrons or photons
+        const xAOD::Egamma *egamma = dynamic_cast<const xAOD::Egamma *>(reco);
+        if (!egamma)
+          // this should never happen
+          throw std::runtime_error("Failed to cast to egamma object");
+        const xAOD::IParticle *cluster = egamma->caloCluster();
+        if (cluster)
+          reco = cluster;
+        else
+          ATH_MSG_WARNING("Cannot retrieve egamma object's primary calorimeter cluster, will match to the egamma object");
+        match = true;
+      }
+      if (match)
+        match = fastDR(reco->eta(), reco->phi(), online->eta(), online->phi(), drThreshold);
+      cacheItr = cache.insert(std::make_pair(linkIndices, match)).first;
+    }
+    return cacheItr->second;
+  }
+
 } // namespace Trig
\ No newline at end of file
diff --git a/Trigger/TrigAnalysis/TriggerMatchingTool/TriggerMatchingTool/R3MatchingTool.h b/Trigger/TrigAnalysis/TriggerMatchingTool/TriggerMatchingTool/R3MatchingTool.h
index 00ee2d004b3b..71460b103315 100644
--- a/Trigger/TrigAnalysis/TriggerMatchingTool/TriggerMatchingTool/R3MatchingTool.h
+++ b/Trigger/TrigAnalysis/TriggerMatchingTool/TriggerMatchingTool/R3MatchingTool.h
@@ -35,6 +35,16 @@ namespace Trig
 
   private:
     ToolHandle<TrigDecisionTool> m_trigDecTool;
+    bool matchObjects(
+        const xAOD::IParticle *reco,
+        const ElementLink<xAOD::IParticleContainer> &onlineLink,
+        std::map<std::pair<uint32_t, uint32_t>, bool> &cache,
+        double drThreshold) const;
+
+    // Internal functions
+    /// Inherited from the interface but does nothing
+    virtual MatchingImplementation *impl() const override { return nullptr; }
+
   }; //> end class R3MatchingTool
 } // namespace Trig
 
diff --git a/Trigger/TrigSteer/TrigCompositeUtils/CMakeLists.txt b/Trigger/TrigSteer/TrigCompositeUtils/CMakeLists.txt
index dca99bc423be..981108688e02 100644
--- a/Trigger/TrigSteer/TrigCompositeUtils/CMakeLists.txt
+++ b/Trigger/TrigSteer/TrigCompositeUtils/CMakeLists.txt
@@ -13,7 +13,8 @@ endif()
 atlas_add_library( TrigCompositeUtilsLib
   TrigCompositeUtils/*.h TrigCompositeUtils/*.icc Root/*.cxx
   PUBLIC_HEADERS TrigCompositeUtils
-  LINK_LIBRARIES TrigConfHLTUtilsLib CxxUtils AsgMessagingLib AsgDataHandlesLib AsgTools TrigDecisionInterface xAODBase xAODTrigger ${extra_libs}  )
+  LINK_LIBRARIES TrigConfHLTUtilsLib CxxUtils AsgMessagingLib AsgDataHandlesLib AsgTools TrigDecisionInterface xAODBase xAODTrigger ${extra_libs}
+  PRIVATE_LINK_LIBRARIES TrigSteeringEvent )
 
 if( NOT XAOD_STANDALONE )
   atlas_add_component( TrigCompositeUtils
diff --git a/Trigger/TrigSteer/TrigCompositeUtils/Root/Combinations.cxx b/Trigger/TrigSteer/TrigCompositeUtils/Root/Combinations.cxx
new file mode 100644
index 000000000000..e081c942ecd2
--- /dev/null
+++ b/Trigger/TrigSteer/TrigCompositeUtils/Root/Combinations.cxx
@@ -0,0 +1,77 @@
+#include "TrigCompositeUtils/Combinations.h"
+#include <tuple>
+#include <algorithm>
+
+namespace TrigCompositeUtils
+{
+  Combinations::Combinations(
+      const std::vector<std::size_t> &legMultiplicities,
+      const std::vector<VecLInfo_t> &legFeatures,
+      const std::function<bool(const VecLInfo_t &)> &filter)
+      : m_filter(filter),
+        m_legMultiplicities(legMultiplicities),
+        m_legFeatures(legFeatures)
+  {
+    if (legMultiplicities.size() != legFeatures.size())
+      throw std::invalid_argument("Different numbers of multiplicities and features provided");
+  }
+
+  Combinations::Combinations(
+      const std::vector<std::size_t> &legMultiplicities,
+      const std::vector<VecLInfo_t> &legFeatures,
+      FilterType filter)
+      : Combinations(legMultiplicities, legFeatures, getFilter(filter))
+  {
+  }
+
+  Combinations::Combinations(const std::function<bool(const VecLInfo_t &)> &filter)
+      : Combinations(std::vector<std::size_t>{}, std::vector<VecLInfo_t>{}, filter)
+  {
+  }
+
+  Combinations::Combinations(FilterType filter) : Combinations(getFilter(filter)) {}
+
+  void Combinations::reserve(std::size_t capacity)
+  {
+    m_legMultiplicities.reserve(capacity);
+    m_legFeatures.reserve(capacity);
+  }
+
+  bool Combinations::empty() const
+  {
+    return begin() == end();
+  }
+
+  std::size_t Combinations::size() const
+  {
+    return std::accumulate(m_legMultiplicities.begin(), m_legMultiplicities.end(), std::size_t(0));
+  }
+
+  void Combinations::addLeg(std::size_t multiplicity, const VecLInfo_t &features)
+  {
+    m_legMultiplicities.push_back(multiplicity);
+    m_legFeatures.push_back(features);
+  }
+
+  void Combinations::addLeg(const VecLInfo_t &features)
+  {
+    addLeg(1, features);
+  }
+
+  IPartCombItr Combinations::begin() const
+  {
+    // Build up the constructor arguments
+    std::vector<std::tuple<std::size_t, VecLInfo_t::const_iterator, VecLInfo_t::const_iterator>> args;
+    for (std::size_t ii = 0; ii < m_legFeatures.size(); ++ii)
+      args.push_back(std::make_tuple(
+          m_legMultiplicities.at(ii),
+          m_legFeatures.at(ii).begin(),
+          m_legFeatures.at(ii).end()));
+    return IPartCombItr(args, m_filter);
+  }
+
+  IPartCombItr Combinations::end() const
+  {
+    return IPartCombItr();
+  }
+} // namespace TrigCompositeUtils
\ No newline at end of file
diff --git a/Trigger/TrigSteer/TrigCompositeUtils/Root/IPartCombItr.cxx b/Trigger/TrigSteer/TrigCompositeUtils/Root/IPartCombItr.cxx
new file mode 100644
index 000000000000..232df7007c3e
--- /dev/null
+++ b/Trigger/TrigSteer/TrigCompositeUtils/Root/IPartCombItr.cxx
@@ -0,0 +1,173 @@
+#include "TrigCompositeUtils/IPartCombItr.h"
+#include <set>
+#include "TrigSteeringEvent/TrigRoiDescriptorCollection.h"
+
+namespace TrigCompositeUtils
+{
+  bool uniqueRoIs(const std::vector<LinkInfo<xAOD::IParticleContainer>> &links)
+  {
+    std::set<std::pair<uint32_t, uint32_t>> seen;
+    for (const auto &info : links)
+    {
+      LinkInfo<TrigRoiDescriptorCollection> roi = findLink<TrigRoiDescriptorCollection>(info.source, "initialRoI");
+      if (!seen.insert(std::make_pair(roi.link.persKey(), roi.link.persIndex())).second)
+        // Insert returns false if that item already exists in it
+        return false;
+    }
+    return true;
+  }
+
+  std::function<bool(const std::vector<LinkInfo<xAOD::IParticleContainer>> &)> getFilter(FilterType filter)
+  {
+    switch (filter)
+    {
+    case FilterType::All:
+      return [](const std::vector<LinkInfo<xAOD::IParticleContainer>> &) { return true; };
+    case FilterType::UniqueRoIs:
+      return uniqueRoIs;
+    default:
+      throw std::runtime_error("Unhandled FilterType enum value!");
+      return {};
+    }
+  }
+
+  IPartCombItr::IPartCombItr(
+      const std::vector<std::tuple<std::size_t, LInfoItr_t, LInfoItr_t>> &pieces,
+      const std::function<bool(const VecLInfo_t &)> filter)
+      : m_filter(filter)
+  {
+    m_itrs.reserve(pieces.size());
+    auto currentItr = std::back_inserter(m_current);
+    for (const auto &tup : pieces)
+    {
+      LInfoItr_t begin = std::get<1>(tup);
+      LInfoItr_t end = std::get<2>(tup);
+      m_itrs.push_back(std::make_pair(
+          KFromNItr(std::get<0>(tup), std::distance(begin, end)), begin));
+      const KFromNItr &idxItr = std::get<0>(m_itrs.back());
+      if (idxItr.exhausted())
+        for (std::size_t ii = 0; ii < idxItr.size(); ++ii)
+          *(currentItr++) = LinkInfo<xAOD::IParticleContainer>{};
+      else
+        std::transform(
+            idxItr->begin(), idxItr->end(), currentItr,
+            [begin](std::size_t idx) { return *(begin + idx); });
+    }
+    // make sure that the starting set makes sense
+    if (!exhausted() && !m_filter(m_current))
+      this->operator++();
+  }
+
+  IPartCombItr::IPartCombItr(
+      const std::vector<std::tuple<std::size_t, LInfoItr_t, LInfoItr_t>> &pieces,
+      FilterType filter)
+      : IPartCombItr(pieces, getFilter(filter))
+  {
+  }
+
+  IPartCombItr::IPartCombItr(const std::function<bool(const VecLInfo_t &)> filter)
+      : m_filter(filter) {}
+
+  IPartCombItr::IPartCombItr(FilterType filter)
+      : m_filter(getFilter(filter)) {}
+
+  void IPartCombItr::reset()
+  {
+    // Reset each individual iterator and set our current value accordingly
+    auto outItr = m_current.begin();
+    for (auto &itrPair : m_itrs)
+    {
+      KFromNItr &idxItr = itrPair.first;
+      idxItr.reset();
+      if (idxItr.exhausted())
+        for (std::size_t ii = 0; ii < idxItr.size(); ++ii)
+          *(outItr++) = LinkInfo<xAOD::IParticleContainer>{};
+      else
+        std::transform(
+            idxItr->begin(), idxItr->end(), outItr,
+            [begin = itrPair.second](std::size_t idx) { return *(begin + idx); });
+    }
+    // make sure that the starting set makes sense
+    if (!exhausted() && !m_filter(m_current))
+      this->operator++();
+  }
+
+  bool IPartCombItr::exhausted() const
+  {
+    return m_itrs.size() == 0 ||
+           std::any_of(m_itrs.begin(), m_itrs.end(),
+                       [](const auto &itrPair) { return itrPair.first.exhausted(); });
+  }
+
+  IPartCombItr::reference IPartCombItr::operator*() const
+  {
+    if (exhausted())
+      throw std::runtime_error("Dereferencing past-the-end iterator");
+    return m_current;
+  }
+
+  IPartCombItr::pointer IPartCombItr::operator->() const
+  {
+    if (exhausted())
+      throw std::runtime_error("Dereferencing past-the-end iterator");
+    return &m_current;
+  }
+
+  IPartCombItr &IPartCombItr::operator++()
+  {
+    if (exhausted())
+      // Don't iterate an iterator that is already past the end
+      return *this;
+    auto backItr = m_itrs.rbegin();
+    std::size_t step = 0;
+    while (backItr != m_itrs.rend())
+    {
+      KFromNItr &idxItr = backItr->first;
+      step += idxItr.size();
+      if (!backItr++->first.exhausted())
+      {
+        // This is the starting point for a good combination
+        // We need to update the current value
+        auto outItr = m_current.end() - step;
+        for (std::size_t idx : *idxItr)
+          *(outItr) = *(backItr->second + idx);
+
+        // Any iterators we passed by up to this point were exhausted so we have
+        // to reset them before we use their values
+        for (auto itr = backItr.base(); itr != m_itrs.end(); ++itr)
+        {
+          itr->first.reset();
+          for (std::size_t idx : *itr->first)
+            *(outItr++) = *(itr->second + idx);
+        }
+        break;
+      }
+      ++backItr;
+    }
+    if (!m_filter(m_current))
+      // If we fail the filter condition, advance the iterator again
+      this->operator++();
+    return *this;
+  }
+
+  IPartCombItr IPartCombItr::operator++(int)
+  {
+    IPartCombItr ret(*this);
+    this->operator++();
+    return ret;
+  }
+
+  bool IPartCombItr::operator==(const IPartCombItr &other) const
+  {
+    // All past-the-end iterators compare equal
+    if (exhausted() && other.exhausted())
+      return true;
+    return m_itrs == other.m_itrs;
+  }
+
+  bool IPartCombItr::operator!=(const IPartCombItr &other) const
+  {
+    return !(*this == other);
+  }
+
+} // namespace TrigCompositeUtils
\ No newline at end of file
diff --git a/Trigger/TrigSteer/TrigCompositeUtils/Root/KFromNItr.cxx b/Trigger/TrigSteer/TrigCompositeUtils/Root/KFromNItr.cxx
new file mode 100644
index 000000000000..e33a17c1862b
--- /dev/null
+++ b/Trigger/TrigSteer/TrigCompositeUtils/Root/KFromNItr.cxx
@@ -0,0 +1,109 @@
+#include "TrigCompositeUtils/KFromNItr.h"
+#include <numeric>
+#include <stdexcept>
+
+namespace TrigCompositeUtils {
+  KFromNItr::KFromNItr(std::size_t k, std::size_t N)
+      : m_N(N),
+        m_current(k, 0)
+  {
+    // Fill the iterator with the first combination
+    reset();
+  }
+
+  void KFromNItr::reset()
+  {
+    // The first combination is just an ascending sequence starting from 0
+    std::iota(m_current.begin(), m_current.end(), 0);
+  }
+
+  bool KFromNItr::exhausted() const
+  {
+    return size() == 0 || m_current.back() >= m_N;
+  }
+
+  KFromNItr::reference KFromNItr::operator*() const
+  {
+    if (exhausted())
+      throw std::runtime_error("Dereferencing past-the-end iterator!");
+    return m_current;
+  }
+
+  KFromNItr::pointer KFromNItr::operator->() const
+  {
+    if (exhausted())
+      throw std::runtime_error("Dereferencing past-the-end iterator!");
+    return &m_current;
+  }
+
+  KFromNItr &KFromNItr::operator++()
+  {
+    if (exhausted())
+      // Don't iterate an iterator that is already past the end
+      return *this;
+    // All index combinations pointed to by this iterator are in ascending order.
+    // In order to generate the next combination, go through the following steps,
+    // starting with the highest (last) index
+    //
+    // 1. Increment the current index
+    // 2. If it is less than its end point then go to step 4
+    // 3. If it is equal to its end point then shift the current index to be the
+    //    before this one and go to step 1.
+    // 4. Set each index after the current one to larger than its preceding
+    //    index's value
+    //
+    // For an example, consider picking 3 values from the range 0 - 4.
+    // The maximum value here is '5'
+    // The combination starts at  (0, 1, 2). For the first two increments only the
+    // third index changes as the comparison in step 2 is always true and the
+    // iterator yields the sequence (0, 1, 3), (0, 1, 4).
+    // After this, the third index increments to 5 (the max value) and so we
+    // increment the second index. As this is now less than its maximum value (4)
+    // we reset all following indices in an ascending sequence, yielding (0, 2, 3)
+
+    // Start from the last index
+    auto backItr = m_current.rbegin();
+    // the end point decreases as we go back through the sequence
+    std::size_t end = m_N;
+    while (backItr != m_current.rend())
+    {
+      // This statement increments the value pointed to by backItr, then checks
+      // that it's still below its maximum value
+      if (++(*backItr) < end)
+      {
+        // step 4 - this value is now OK. We have to set each following element
+        // in increasing sequence
+        // backItr.base() returns an iterator pointing to the element *after* this
+        // one
+        std::iota(backItr.base(), m_current.end(), *backItr + 1);
+        break;
+      }
+      ++backItr;
+      --end;
+    }
+    return *this;
+  }
+
+  KFromNItr KFromNItr::operator++(int)
+  {
+    /// Copy our current state
+    KFromNItr ret(*this);
+    // step us past this point
+    this->operator++();
+    return ret;
+  }
+
+  bool KFromNItr::operator==(const KFromNItr &other) const
+  {
+    // All past-the-end iterators compare equal
+    if (exhausted() && other.exhausted())
+      return true;
+    // Otherwise make sure that the iterators point to the same combination
+    return m_current == other.m_current;
+  }
+
+  bool KFromNItr::operator!=(const KFromNItr &other) const
+  {
+    return !(*this == other);
+  }
+} //> end namespace TrigCompositeUtils
\ No newline at end of file
diff --git a/Trigger/TrigSteer/TrigCompositeUtils/TrigCompositeUtils/Combinations.h b/Trigger/TrigSteer/TrigCompositeUtils/TrigCompositeUtils/Combinations.h
new file mode 100644
index 000000000000..66d87d09c6e2
--- /dev/null
+++ b/Trigger/TrigSteer/TrigCompositeUtils/TrigCompositeUtils/Combinations.h
@@ -0,0 +1,50 @@
+#ifndef TRIGCOMPOSITEUTILS_COMBINATIONS_H
+#define TRIGCOMPOSITEUTILS_COMBINATIONS_H
+
+#include "TrigCompositeUtils/IPartCombItr.h"
+#include "TrigCompositeUtils/TrigCompositeUtils.h"
+#include "xAODBase/IParticleContainer.h"
+#include <functional>
+#include <vector>
+
+namespace TrigCompositeUtils
+{
+
+  class Combinations
+  {
+  public:
+    using VecLInfo_t = std::vector<LinkInfo<xAOD::IParticleContainer>>;
+    Combinations(
+        const std::vector<std::size_t> &legMultiplicities,
+        const std::vector<VecLInfo_t> &legFeatures,
+        const std::function<bool(const VecLInfo_t &)> &filter);
+
+    Combinations(
+        const std::vector<std::size_t> &legMultiplicities,
+        const std::vector<VecLInfo_t> &legFeatures,
+        FilterType filter = FilterType::All);
+
+    Combinations(const std::function<bool(const VecLInfo_t &)> &filter);
+
+    Combinations(FilterType filter = FilterType::All);
+
+    void reserve(std::size_t capacity);
+
+    bool empty() const;
+
+    std::size_t size() const;
+
+    void addLeg(std::size_t multiplicity, const VecLInfo_t &features);
+    void addLeg(const VecLInfo_t &features);
+
+    IPartCombItr begin() const;
+    IPartCombItr end() const;
+
+  private:
+    std::function<bool(const VecLInfo_t &)> m_filter;
+    std::vector<std::size_t> m_legMultiplicities;
+    std::vector<VecLInfo_t> m_legFeatures;
+  }; //> end class Combinations
+} // namespace TrigCompositeUtils
+
+#endif //> TRIGCOMPOSITEUTILS_COMBINATIONS_H
\ No newline at end of file
diff --git a/Trigger/TrigSteer/TrigCompositeUtils/TrigCompositeUtils/IPartCombItr.h b/Trigger/TrigSteer/TrigCompositeUtils/TrigCompositeUtils/IPartCombItr.h
new file mode 100644
index 000000000000..aa36ccee3cf6
--- /dev/null
+++ b/Trigger/TrigSteer/TrigCompositeUtils/TrigCompositeUtils/IPartCombItr.h
@@ -0,0 +1,139 @@
+#ifndef TRIGCOMPOSITEUTILS_IPARTCOMBITR_H
+#define TRIGCOMPOSITEUTILS_IPARTCOMBITR_H
+
+#include <iterator>
+#include <vector>
+#include <tuple>
+#include <utility>
+#include <algorithm>
+#include <functional>
+#include "xAODBase/IParticleContainer.h"
+#include "TrigCompositeUtils/KFromNItr.h"
+#include "TrigCompositeUtils/TrigCompositeUtils.h"
+
+namespace TrigCompositeUtils
+{
+  enum class FilterType
+  {
+    /// Allow all combinations
+    All,
+    /// Do not allow any two objects to share an RoI
+    UniqueRoIs,
+  };
+  /// Helper function that returns true if no objects share an initial RoI
+  bool uniqueRoIs(const std::vector<LinkInfo<xAOD::IParticleContainer>> &links);
+
+  /// Get a lambda corresponding to the specified FilterType enum.
+  std::function<bool(const std::vector<LinkInfo<xAOD::IParticleContainer>> &)> getFilter(FilterType filter);
+  class IPartCombItr
+  {
+  public:
+    using VecLInfo_t = std::vector<LinkInfo<xAOD::IParticleContainer>>;
+    using LInfoItr_t = VecLInfo_t::const_iterator;
+
+    using iterator_category = std::input_iterator_tag;
+    using value_type = VecLInfo_t;
+    using reference = const value_type &;
+    using pointer = const value_type *;
+    using difference_type = std::ptrdiff_t;
+
+    /**
+     * @brief The direct constructor
+     * 
+     * Takes a vector of tuples.
+     * The last two elements of each tuple are start and end iterators
+     * describing ranges of IParticleContainer LinkInfos.
+     * The first element describes how many particles from this range are used
+     * in the output combinations
+     */
+    IPartCombItr(
+        const std::vector<std::tuple<std::size_t, LInfoItr_t, LInfoItr_t>> &pieces,
+        const std::function<bool(const VecLInfo_t &)> filter);
+
+    /**
+     * @brief The direct constructor
+     * 
+     * Takes a vector of tuples.
+     * The last two elements of each tuple are start and end iterators
+     * describing ranges of IParticleContainer LinkInfos.
+     * The first element describes how many particles from this range are used
+     * in the output combinations
+     */
+    IPartCombItr(
+        const std::vector<std::tuple<std::size_t, LInfoItr_t, LInfoItr_t>> &pieces,
+        FilterType filter = FilterType::All);
+
+    /// Base case constructor for the variadic constructors
+    IPartCombItr(const std::function<bool(const VecLInfo_t &)> filter);
+
+    /// Base case constructor for the variadict constructors
+    IPartCombItr(FilterType filter = FilterType::All);
+
+    template <typename... Ts>
+    IPartCombItr(std::size_t k, const LInfoItr_t &begin, const LInfoItr_t &end, Ts &&... args)
+        : IPartCombItr(std::forward<Ts>(args)...)
+    {
+      m_itrs.insert(m_itrs.begin(), std::make_pair(KFromNItr(k, std::distance(begin, end)), begin));
+      m_current.insert(m_current.begin(), k, {});
+      const KFromNItr &idxItr = *std::get<0>(m_itrs.front());
+      if (!idxItr.exhausted())
+        std::transform(idxItr->begin(), idxItr->end(), m_current.begin(),
+                       [begin](std::size_t idx) { return *(begin + idx); });
+    }
+
+    template <typename... Ts>
+    IPartCombItr(const LInfoItr_t &begin, const LInfoItr_t &end, Ts &&... args)
+        : IPartCombItr(1, begin, end, std::forward<Ts>(args)...)
+    {
+    }
+
+    template <typename... Ts>
+    IPartCombItr(std::size_t k, const VecLInfo_t &linfos, Ts &&... args)
+        : IPartCombItr(k, linfos.begin(), linfos.end(), std::forward<Ts>(args)...)
+    {
+    }
+
+    template <typename... Ts>
+    IPartCombItr(const VecLInfo_t &linfos, Ts &&... args)
+        : IPartCombItr(1, linfos.begin(), linfos.end(), std::forward<Ts>(args)...)
+    {
+    }
+
+    /// The size of each combination
+    std::size_t size() const { return m_current.size(); }
+
+    /**
+     * @brief Reset the iterator to its starting point
+     * 
+     * Note that this resets all component iterators, so it will not reset this
+     * object to its starting point if any of the iterators passed to it were not
+     * initially in their starting positions
+     */
+    void reset();
+
+    /// True if this iterator is past the end
+    bool exhausted() const;
+
+    /// Dereference
+    reference operator*() const;
+    pointer operator->() const;
+
+    /// Pre-increment operator
+    IPartCombItr &operator++();
+
+    /// Post-increment operator
+    IPartCombItr operator++(int);
+
+    /// Iterator comparison functions
+    bool operator==(const IPartCombItr &other) const;
+    bool operator!=(const IPartCombItr &other) const;
+
+  private:
+    std::function<bool(const VecLInfo_t &)> m_filter;
+    std::vector<std::pair<KFromNItr, LInfoItr_t>> m_itrs;
+    VecLInfo_t m_current;
+
+  }; //> end class IPartCombItr
+} // namespace TrigCompositeUtils
+
+#endif //> !TRIGCOMPOSITEUTILS_IPARTCOMBITR_H
\ No newline at end of file
diff --git a/Trigger/TrigSteer/TrigCompositeUtils/TrigCompositeUtils/KFromNItr.h b/Trigger/TrigSteer/TrigCompositeUtils/TrigCompositeUtils/KFromNItr.h
new file mode 100644
index 000000000000..2aa61f3c4851
--- /dev/null
+++ b/Trigger/TrigSteer/TrigCompositeUtils/TrigCompositeUtils/KFromNItr.h
@@ -0,0 +1,77 @@
+#ifndef TRIGCOMPOSITEUTILS_KFROMNITR_H
+#define TRIGCOMPOSITEUTILS_KFROMNITR_H
+
+#include <iterator>
+#include <vector>
+
+namespace TrigCompositeUtils {
+  /**
+  * @brief Iterates over all combinations of k values chosen from a range n
+  * 
+  * Generates all distinct combinations of k distinct values less than N.
+  * 
+  * This is an input iterator as it's impossible to satisfy the multi-pass
+  * guarantee in the forward iterator category where the return value is
+  * generated on the fly.
+  * However, this still follows the forward iterator convention that all
+  * past-the-end iterators compare equal, and that a default constructed
+  * iterator counts as a past-the-end iterator.
+  * 
+  * Combinations are generated in ascending order, with the indices in
+  * ascending order. The highest index is always incremented first if possible.
+  * Therefore the combinations generated by an iterator constructed as
+  * KFromNItr(3, 4);
+  * would be
+  * (0, 1, 2)
+  * (0, 1, 3)
+  * (0, 2, 3)
+  * (1, 2, 3)
+  */
+  class KFromNItr
+  {
+  public:
+    /// Iterator traits
+    using iterator_category = std::input_iterator_tag;
+    using value_type = std::vector<std::size_t>;
+    using reference = const value_type &;
+    using pointer = const value_type *;
+    using difference_type = std::ptrdiff_t;
+
+    /// Default constructor creates a generic past-the-end iterator
+    KFromNItr() = default;
+
+    /// Construct the iterator choosing k from N
+    KFromNItr(std::size_t k, std::size_t N);
+
+    /// The size of each combination (k)
+    std::size_t size() const { return m_current.size(); }
+
+    /// Reset the iterator to its start position
+    void reset();
+
+    /// True if this iterator is past the end
+    bool exhausted() const;
+
+    /// Dereference
+    reference operator*() const;
+    pointer operator->() const;
+
+    /// Pre-increment operator
+    KFromNItr &operator++();
+
+    /// Post-increment operator
+    KFromNItr operator++(int);
+
+    /// Iterator comparison functions
+    bool operator==(const KFromNItr &other) const;
+    bool operator!=(const KFromNItr &other) const;
+
+  private:
+    /// The number of indices
+    std::size_t m_N;
+    /// The current combination
+    std::vector<std::size_t> m_current;
+  }; //> end class KFromNItr
+} //> end namespace TrigCompositeUtils
+
+#endif //> !TRIGCOMPOSITEUTILS_KFROMNITR_H
\ No newline at end of file
-- 
GitLab