diff --git a/Phys/ParticleCombiners/CombKernel/ParticleCombination.h b/Phys/ParticleCombiners/CombKernel/ParticleCombination.h
index 3a331a674392d700180807f4467f62989fb5ba86..1ae462cfe1fc2b3d08ab7877b8e234073a5da187 100644
--- a/Phys/ParticleCombiners/CombKernel/ParticleCombination.h
+++ b/Phys/ParticleCombiners/CombKernel/ParticleCombination.h
@@ -9,27 +9,58 @@
 * or submit itself to any jurisdiction.                                       *
 \*****************************************************************************/
 #pragma once
-#include "ParticleCombiner.h"
 #include <array>
+#include <cmath>
 #include <functional>
 #include <limits>
 
+#include "Kernel/HeaderMapping.h"
+#include "Kernel/IDistanceCalculator.h"
+
 namespace Combination::details {
   using fp = float;
 }
 
-template <typename Particle, std::size_t N>
+struct IDistanceCalculator;
+
+/** @brief A list of particle objects.
+ *
+ * A 'combination' is simply an ordered list of particle objects. This
+ * combination object holds such a list and has several member functions for
+ * computing kinematic and geometric 'properties' of the list.
+ *
+ * The 'four momentum' of a combination is defined as the sum of the four
+ * momentum of the particles held by it. Most other intrinsic properties can be
+ * derived from this, such as the 'invariant mass' of a combination.
+ *
+ * Some other methods exist to make this object compatible with the API expected
+ * by the ThOr functor framework, allowing an instance of ParticleCombination to
+ * be passed in as an argument to a functor.
+ */
+template <typename... InputTypes>
 struct ParticleCombination {
+  static constexpr auto N = sizeof...( InputTypes );
+  static_assert( N >= 1 );
+  // Assert that all members of the parameter pack have the same type
+  using Particle_t = std::tuple_element_t<0, std::tuple<InputTypes...>>;
+  static_assert( std::conjunction_v<std::is_same<Particle_t, InputTypes>...> );
+
+  ParticleCombination( IDistanceCalculator const* distance_calculator, IGeometryInfo const& geometry )
+      : m_distance_calculator( distance_calculator ), m_geometry{geometry} {}
 
-  using child_array = std::array<Particle const*, N>;
-  child_array m_children{};
+  /// Pointers to each of the N child objects in this combination
+  std::array<Particle_t const*, N> m_children{};
 
-  NBodyParticleCombinerBase const* m_legacy_accessor;
-  ParticleCombination( NBodyParticleCombinerBase const* legacy_accessor ) : m_legacy_accessor( legacy_accessor ) {}
+  /// Distance calculator used to compute distance-related quantities
+  IDistanceCalculator const* m_distance_calculator;
 
+  /// Geometry information passed to the distance calculator
+  IGeometryInfo const& m_geometry;
+
+  /** @brief Squared-magnitude of the four momentum sum of the children.
+   */
   Combination::details::fp mom2() const {
     using Combination::details::fp;
-    using std::sqrt;
     double px = std::numeric_limits<fp>::min();
     double py = std::numeric_limits<fp>::min();
     double pz = std::numeric_limits<fp>::min();
@@ -41,6 +72,29 @@ struct ParticleCombination {
     return px * px + py * py + pz * pz;
   }
 
+  /** @brief Magnitude of the four momentum sum of the children.
+   */
+  Combination::details::fp p() const {
+    using std::sqrt;
+    return sqrt( mom2() );
+  }
+
+  /** @brief Transverse momentum component of the four momentum sum of the children.
+   */
+  Combination::details::fp pt() const {
+    using Combination::details::fp;
+    using std::sqrt;
+    double px = std::numeric_limits<fp>::min();
+    double py = std::numeric_limits<fp>::min();
+    for ( auto* i : m_children ) {
+      px += i->momentum().px();
+      py += i->momentum().py();
+    }
+    return sqrt( px * px + py * py );
+  }
+
+  /** @brief Invariant mass of the four momentum sum of the children.
+   */
   Combination::details::fp mass() const {
     using Combination::details::fp;
     using std::sqrt;
@@ -49,37 +103,47 @@ struct ParticleCombination {
     return sqrt( energy * energy - mom2() );
   }
 
+  /** @brief Distance of closest approach between two children.
+   */
   template <std::size_t I1, std::size_t I2, typename DistanceCalculator>
   auto doca( DistanceCalculator const& ) const {
     static_assert( I1 < N && I2 < N );
     double current_distance{};
-    m_legacy_accessor->m_dist_calc->distance( m_children[I1], m_children[I2], current_distance ).ignore();
+    m_distance_calculator->distance( m_children[I1], m_children[I2], current_distance, m_geometry ).ignore();
     return current_distance;
   }
 
+  /** @brief Return true if the maximum distance of closest approach across all
+   *         possible pairs of children is less than the threshold.
+   */
   template <typename DistanceCalculator>
   bool maxdocacut( DistanceCalculator const&, Combination::details::fp thresh ) const {
     using Combination::details::fp;
     for ( auto i1 = m_children.begin(); i1 != m_children.end(); ++i1 ) {
       for ( auto i2 = std::next( i1 ); i2 != m_children.end(); ++i2 ) {
         double current_distance = std::numeric_limits<fp>::min();
-        m_legacy_accessor->m_dist_calc->distance( *i1, *i2, current_distance, m_legacy_accessor->geometry() );
+        m_distance_calculator->distance( *i1, *i2, current_distance, m_geometry ).ignore();
+        // Short circuit if we find a DOCA above the threshold
         if ( current_distance > thresh ) { return false; };
       }
     }
     return true;
   }
 
+  /** @brief Distance of closest approach chi^2 between two children.
+   */
   template <std::size_t I1, std::size_t I2, typename DistanceCalculator>
   auto docachi2( DistanceCalculator const& ) const {
     static_assert( I1 < N && I2 < N );
     double current_chi2{}, current_distance{};
-    m_legacy_accessor->m_dist_calc
-        ->distance( m_children[I1], m_children[I2], current_distance, current_chi2, m_legacy_accessor->geometry() )
+    m_distance_calculator->distance( m_children[I1], m_children[I2], current_distance, current_chi2, m_geometry )
         .ignore();
     return current_chi2;
   }
 
+  /** @brief Maximum distance of closest approach chi^2 across all possible pairs
+   *         of children is less.
+   */
   template <typename DistanceCalculator>
   Combination::details::fp maxdocachi2( DistanceCalculator const& ) const {
     using Combination::details::fp;
@@ -88,13 +152,16 @@ struct ParticleCombination {
       for ( auto i2 = std::next( i1 ); i2 != m_children.end(); ++i2 ) {
         double current_chi2     = std::numeric_limits<fp>::min();
         double current_distance = std::numeric_limits<fp>::min();
-        m_legacy_accessor->m_dist_calc->distance( *i1, *i2, current_distance, current_chi2 );
+        m_distance_calculator->distance( *i1, *i2, current_distance, current_chi2, m_geometry ).ignore();
         chi2 = std::max( chi2, current_chi2 );
       }
     }
     return chi2;
   }
 
+  /** @brief Return true if the maximum distance of closest approach chi^2 across
+   *         all possible pairs of children is less than the threshold.
+   */
   template <typename DistanceCalculator>
   bool maxdocachi2cut( DistanceCalculator const&, Combination::details::fp thresh ) const {
     using Combination::details::fp;
@@ -102,9 +169,8 @@ struct ParticleCombination {
       for ( auto i2 = std::next( i1 ); i2 != m_children.end(); ++i2 ) {
         double current_chi2     = std::numeric_limits<fp>::min();
         double current_distance = std::numeric_limits<fp>::min();
-        m_legacy_accessor->m_dist_calc
-            ->distance( *i1, *i2, current_distance, current_chi2, m_legacy_accessor->geometry() )
-            .ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ );
+        m_distance_calculator->distance( *i1, *i2, current_distance, current_chi2, m_geometry ).ignore();
+        // Short circuit if we find a DOCA chi^2 above the threshold
         if ( current_chi2 > thresh ) { return false; };
       }
     }
@@ -112,12 +178,18 @@ struct ParticleCombination {
   }
 };
 
-template <typename Particle, std::size_t N>
-auto begin( ParticleCombination<Particle, N> const& combination ) {
+template <typename... InputTypes>
+auto begin( ParticleCombination<InputTypes...> const& combination ) {
   return combination.m_children.begin();
 }
 
-template <typename Particle, std::size_t N>
-auto end( ParticleCombination<Particle, N> const& combination ) {
+template <typename... InputTypes>
+auto end( ParticleCombination<InputTypes...> const& combination ) {
   return combination.m_children.end();
 }
+
+// Register headers
+template <typename... child_ts>
+struct LHCb::header_map<ParticleCombination<child_ts...>> {
+  static constexpr auto value = ( LHCb::header_map_v<child_ts> + ... ) + "CombKernel/ParticleCombination.h";
+};
diff --git a/Phys/ParticleCombiners/CombKernel/ParticleCombiner.h b/Phys/ParticleCombiners/CombKernel/ParticleCombiner.h
index 18b2fce0ffb2efbb0c2fa64c1723dc8c9e4b9912..043b95a2686cb0f2cb618794fee5f780560867fc 100644
--- a/Phys/ParticleCombiners/CombKernel/ParticleCombiner.h
+++ b/Phys/ParticleCombiners/CombKernel/ParticleCombiner.h
@@ -1,5 +1,5 @@
 /*****************************************************************************\
-* (c) Copyright 2000-2019 CERN for the benefit of the LHCb Collaboration      *
+* (c) Copyright 2000-2021 CERN for the benefit of the LHCb Collaboration      *
 *                                                                             *
 * This software is distributed under the terms of the GNU General Public      *
 * Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".   *
@@ -8,262 +8,633 @@
 * granted to it by virtue of its status as an Intergovernmental Organization  *
 * or submit itself to any jurisdiction.                                       *
 \*****************************************************************************/
-
 #pragma once
+#include <algorithm>
 
-// legacy stuff
-#include "Kernel/IDistanceCalculator.h"
+#include <boost/mp11.hpp>
 
-#include "DetDesc/DetectorElement.h"
-#include "DetDesc/GenericConditionAccessorHolder.h"
+#include "CombKernel/ParticleCombination.h"
 #include "Event/Particle.h"
-#include "Event/RecVertex.h"
-#include "Event/Track_v1.h"
-#include "Functors/Core.h"
-#include "Functors/Function.h"
-#include "Functors/IFactory.h"
-#include "Functors/with_functor_maps.h"
-#include "Functors/with_functors.h"
 #include "Kernel/GetDecay.h"
 #include "Kernel/ICheckOverlap.h"
 #include "Kernel/IDecodeSimpleDecayString.h"
+#include "Kernel/IDistanceCalculator.h"
 #include "Kernel/IParticleCombiner.h"
 #include "Kernel/IParticlePropertySvc.h"
-#include "Kernel/ParticleProperty.h"
-#include "PrKernel/PrSelection.h"
 
-#include "GaudiAlg/MergingTransformer.h"
-#include "GaudiAlg/Transformer.h"
+#include "DetDesc/DetectorElement.h"
+#include "DetDesc/GenericConditionAccessorHolder.h"
 
-#include <boost/container/small_vector.hpp>
+#include "GaudiAlg/Transformer.h"
 
-template <typename particle, std::size_t N>
-struct ParticleCombination;
+#include "Functors/with_functors.h"
 
 namespace Combiner {
   namespace constants {
-    static constexpr auto max_particlecontainer_size    = 500;   // TODO optimize
-    static constexpr auto max_outcontainer_size         = 200;   // TODO optimize
-    static constexpr auto max_combinationcontainer_size = 10000; // TODO optimize
-
-    constexpr auto make_combiner_map = []( auto i ) {
-      if constexpr ( i == 12 )
-        return "Comb12Cut";
-      else if constexpr ( i == 13 )
-        return "Comb13Cut";
-      else if constexpr ( i == 14 )
-        return "Comb14Cut";
-      else if constexpr ( i == 123 )
-        return "Comb123Cut";
-      else
-        return "Comb1234Cut";
-    };
-
-    template <std::size_t i>
-    constexpr static auto combiner_name = make_combiner_map( std::integral_constant<int, i>{} );
-
+    // TODO optimise based on real-world usage
+    static constexpr auto max_outcontainer_size = 200;
   } // namespace constants
 
-  template <typename T, std::size_t N, std::size_t i>
+  // Get the name of the Mth combination cut of an N-particle combiner
+  // Very ugly, but not worth the effort to do something beautiful and constexpr
+  template <std::size_t N, std::size_t M>
+  std::string combinationCutName() {
+    static_assert( N >= 2 );
+    static_assert( M + 1 < N );
+    std::string name{"Combination"};
+    if ( M + 2 < N ) {
+      // M = 0 => "12"
+      // M = 1 => "123"
+      // etc.
+      name += "12";
+      for ( auto i = 0ul; i < M; ++i ) { name += std::to_string( i + 3 ); }
+    }
+    name += "Cut";
+    return name;
+  }
+
+  // Tag types for combination functors
+  template <typename... Ts>
   struct CombTypes {
-    using Combination             = ParticleCombination<T, N>;
-    constexpr static auto Headers = std::array{"CombKernel/ParticleCombination.h"};
-    struct Cut {
-      using Signature                    = bool( Combination const& ); // TODO mask_v instead of bool
-      constexpr static auto PropertyName = constants::combiner_name<i>;
-      constexpr static auto ExtraHeaders = Headers;
-    };
-  };
+    static constexpr auto N = sizeof...( Ts );
+    static_assert( N >= 2 );
 
-  template <typename T>
-  struct MotherTypes {
-    constexpr static auto Headers = std::array{"Event/Particle.h"};
+    // Adapted from https://stackoverflow.com/a/45172773/596068
+    template <template <typename...> class C, std::size_t M, typename... T, std::size_t... I>
+    static std::enable_if_t<( M <= sizeof...( T ) ), C<std::tuple_element_t<I, std::tuple<T...>>...>>
+    sub( std::index_sequence<I...> );
+
+    template <template <typename...> class C, std::size_t M, typename... T>
+    using subpack = decltype( sub<C, M, T...>( std::make_index_sequence<M>{} ) );
+
+    // Helper to return a ParticleCombination with the first M elements of the
+    // parameter pack `Ts`
+    template <std::size_t M>
+    using Combination = subpack<ParticleCombination, M, Ts...>;
+
+    // Combination12Cut (M = 0), Combination123Cut (M = 1), ..., CombinationCut (M = N-1)
+    template <std::size_t M>
     struct Cut {
-      using Out                          = bool;
-      using Signature                    = Out( T const& ); // TODO mask_v instead of bool
-      constexpr static auto PropertyName = "MotherCut";
-      constexpr static auto ExtraHeaders = Headers;
+      static_assert( M < N );
+      using Signature                    = bool( Combination<M + 2> const& );
+      inline static auto    PropertyName = combinationCutName<N, M>();
+      constexpr static auto ExtraHeaders = LHCb::header_map_v<Combination<M + 2>>;
     };
   };
 
+  // Tag types for composite functors (acting on post-vertex-fit particles)
   template <typename T>
-  struct DauTypes {
-    constexpr static auto Headers = std::array{"Event/Particle.h"};
+  struct CompositeTypes {
     struct Cut {
-      using Out                          = bool;
-      using Signature                    = Out( T const& );
-      constexpr static auto PropertyName = "DaughterCuts";
-      constexpr static auto ExtraHeaders = Headers;
+      using Signature                    = bool( T const& );
+      constexpr static auto PropertyName = "CompositeCut";
+      constexpr static auto ExtraHeaders = LHCb::header_map_v<T>;
     };
   };
 
   namespace types {
-    using Vertex_t       = LHCb::Vertex;
-    using Particle_t     = LHCb::Particle;
-    using Particle_ptr_t = LHCb::Particle const*;
-    using Track_t        = LHCb::Event::v1::Track;
-    using Combination4_t = ParticleCombination<LHCb::Particle, 4>;
-    using Combination3_t = ParticleCombination<LHCb::Particle, 3>;
-    using Combination2_t = ParticleCombination<LHCb::Particle, 2>;
-
-    using ParticleCut_t     = Functors::Functor<bool( Particle_t const& )>;
-    using DaughterCuts_t    = std::map<int, ParticleCut_t const*>;
-    using Combination4Cut_t = Functors::Functor<bool( Combination4_t const& )>;
-    using Combination3Cut_t = Functors::Functor<bool( Combination3_t const& )>;
-    using Combination2Cut_t = Functors::Functor<bool( Combination2_t const& )>;
-
-    using ParticleContainer_t =
-        boost::container::small_vector<Particle_t const*, constants::max_particlecontainer_size>;
-    using InputParticleContainer_t = Particle_t::Range;
-    using In_t                     = Gaudi::Functional::details::vector_of_const_<InputParticleContainer_t>;
-    using OutVertices_t            = Vertex_t::Container;
-    using OutParticles_t           = Particle_t::Container;
-    using Out_t                    = std::tuple<OutParticles_t, OutVertices_t>;
-    using FilterOut_t              = std::tuple<bool, OutParticles_t, OutVertices_t>;
-    using functional_base_t        = Gaudi::Functional::MergingMultiTransformerFilter<
-        Out_t( In_t const& ),
-        // Needed for handling of LHCb ConditionHandle.
-        // FIXME should be integrated in functional algorithm, but MergingTransformers
-        // do not allow extra (non vector) inputs
-        Gaudi::Functional::Traits::BaseClass_t<LHCb::DetDesc::ConditionAccessorHolder<FixTESPath<Gaudi::Algorithm>>>>;
-    using DauType    = DauTypes<LHCb::Particle>;
-    using MotherType = MotherTypes<LHCb::Particle>;
-    template <std::size_t N, std::size_t i>
-    using CombType = CombTypes<LHCb::Particle, N, i>;
-
-    using base_t = with_functors<with_functor_maps<functional_base_t, DauType::Cut>, MotherType::Cut>;
+    using Particle_t = LHCb::Particle;
+    using Vertex_t   = std::remove_pointer<decltype( LHCb::Particle{}.endVertex() )>::type;
+
+    using OutVertices_t  = Vertex_t::Container;
+    using OutParticles_t = Particle_t::Container;
+    using Output_t       = std::tuple<OutParticles_t, OutVertices_t>;
+    using FilterOutput_t = std::tuple<bool, OutParticles_t, OutVertices_t>;
+    using CompositeType  = CompositeTypes<Particle_t>;
+
+    template <typename T>
+    using ToRange = typename T::Range;
+
+    template <typename... Ts>
+    struct base_helper {
+      // Assert that all members of the parameter pack have the same type and
+      // that that type is LHCb::Particle
+      using first_t = std::tuple_element_t<0, std::tuple<Ts...>>;
+      static_assert( std::conjunction_v<std::is_same<first_t, Ts>...> );
+      static_assert( std::conjunction_v<std::is_same<first_t, Particle_t>> );
+
+      // Shorthand
+      using FilterTransform = Gaudi::Functional::MultiTransformerFilter<
+          Output_t( ToRange<Ts> const&..., DetectorElement const& ),
+          LHCb::DetDesc::usesBaseAndConditions<GaudiAlgorithm, DetectorElement>>;
+
+      // Bake T and N into aliases that we can pass to mp_product.
+      template <std::size_t M>
+      using CombCut = typename CombTypes<Ts...>::template Cut<M>;
+      template <typename M>
+      using _CombCut = CombCut<M::value>;
+      // Make with_functors<FilterTransform<Ts...>,
+      //                    [cartesian expansion of CombinationCut<A>...],
+      //                    CompositeCut>
+      // Where A is an integer from 0 to N-2
+      using type = boost::mp11::mp_append<
+          with_functors<FilterTransform>,
+          boost::mp11::mp_product<_CombCut,
+                                  boost::mp11::mp_from_sequence<std::make_index_sequence<sizeof...( Ts ) - 1>>>,
+          boost::mp11::mp_list<CompositeType::Cut>>;
+    };
 
+    template <typename... Ts>
+    using base_t = typename base_helper<Ts...>::type;
   } // namespace types
 
 } // namespace Combiner
 
-struct NBodyParticleCombinerBase : Combiner::types::base_t {
-  NBodyParticleCombinerBase( const std::string& name, ISvcLocator* pSvcLocator )
-      : Combiner::types::base_t( name, pSvcLocator, {"Inputs", {}},
-                                 {KeyValue{"OutputParticles", ""}, KeyValue{"OutputVertices", ""}} ) {}
-
-  using base = Combiner::types::base_t;
+template <typename... InputTypes>
+class NBodyParticleCombiner final : public Combiner::types::base_t<InputTypes...> {
+public:
+  /// Number of children in the decay tree we'll be reconstructing
+  static constexpr auto N = sizeof...( InputTypes );
+  using base_t            = typename Combiner::types::base_t<InputTypes...>;
+  using KeyValue          = typename base_t::KeyValue;
+  template <typename T>
+  using InputContainer = typename Combiner::types::ToRange<T>;
+  template <std::size_t M>
+  using CombCut      = typename Combiner::types::base_helper<InputTypes...>::template CombCut<M>;
+  using CompositeCut = typename Combiner::types::CompositeType::Cut;
+  // Container type used to store indices into the input containers
+  using index_vector_t = std::vector<int>;
+
+  NBodyParticleCombiner( const std::string& name, ISvcLocator* pSvcLocator )
+      : base_t{name,
+               pSvcLocator,
+               input_names( std::index_sequence_for<InputTypes...>{} ),
+               {KeyValue{"OutputParticles", ""}, KeyValue{"OutputVertices", ""}}} {}
+
+  /** @brief Return an array of names for our N input properties.
+   *
+   * Returns Input1, ..., InputN, starting from 1 to match the indices shown in
+   * the combination cuts (e.g. Combination12Cut).
+   */
+  template <std::size_t... Is>
+  static auto input_names( std::index_sequence<Is...> ) {
+    return std::array{KeyValue{"Input" + std::to_string( Is + 1 ), ""}...,
+                      KeyValue{"StandardGeometryTop", "/dd/Structure/LHCb"}};
+  }
 
   StatusCode initialize() override {
-    auto sc = base::initialize();
-
-    m_decoder.retrieve().ignore(); // need this?
-    m_overlapTool.retrieve().ignore();
-    m_vertexFitter.retrieve().ignore();
-    m_dist_calc.retrieve().ignore();
-    m_decays = Decays::decays( m_decaydescriptors, &*m_decoder );
-
-    auto const& daughter_cuts = this->template getFunctorMap<Combiner::types::DauType::Cut>();
-
-    std::map<int, LHCb::ParticleProperty const*> children_pids;
-    // TODO only works for one decay now!
+    StatusCode sc = base_t::initialize();
+    if ( sc.isFailure() ) { return sc; }
+
+    // We will parse at most two decays to reconstruct: one if the descriptor is
+    // just `...`, and two if its `[...]cc`
+    m_decays = Decays::decays( m_decay_descriptor, &*m_decay_descriptor_decoder );
+    assert( 0 < m_decays.size() && m_decays.size() <= 2 );
+
+    // For each decay to reconstruct, validate that:
+    // 1. Number of children in descriptor == `N`
+    // 2. Identical PIDs appear consecutively
+    auto decay_index = 0;
     for ( auto const& decay : m_decays ) {
-      for ( auto const& child : decay.children() ) { children_pids.emplace( child.pid().pid(), child.pp() ); }
-    }
-
-    // put daughter cut into the map of <pid, daughter cuts>
-
-    for ( auto const& child : children_pids ) {
-      auto  child_pid = std::get<0>( child );
-      auto* partProp  = std::get<1>( child );
-      // try particle
-      auto iter = std::find_if( daughter_cuts.begin(), daughter_cuts.end(), [&]( auto const& dau ) {
-        return this->m_particlePropSvc->find( dau.first )->pid().pid() == child_pid;
-      } );
-      // try antiparticle
-      if ( iter == daughter_cuts.end() ) {
-        iter = std::find_if( daughter_cuts.begin(), daughter_cuts.end(), [&]( auto const& dau ) {
-          return this->m_particlePropSvc->find( dau.first )->anti()->pid().pid() == child_pid;
-        } );
+      auto const& children = decay.children();
+      if ( children.size() != N ) {
+        throw std::runtime_error( "Cannot combine " + std::to_string( children.size() ) + " children in a " +
+                                  std::to_string( N ) + "-body combiner" );
       }
-      if ( iter == daughter_cuts.end() ) {
-        // if its not there, emplace nullptr to indicate default
-        m_daughter_cuts.emplace( child_pid, nullptr );
-      } else {
-        m_daughter_cuts.emplace( child_pid, &iter->second );
+
+      std::set<int> pids{children[0].pid().pid()};
+      for ( uint i = 1; i < N; ++i ) {
+        int pid = children[i].pid().pid();
+        if ( pids.count( pid ) && children[i - 1].pid().pid() != pid ) {
+          throw std::runtime_error( std::string{"Children with idential PID "} +
+                                    "must appear consecutively in the decay descriptor; e.g.:\n" +
+                                    "YES: B+ -> pi+ pi+ pi-\n" + "NO:  B+ -> pi+ pi- pi+\n" );
+        }
+        pids.emplace( pid );
       }
-      m_daughters_passed.emplace( std::piecewise_construct, std::tuple{child_pid},
-                                  std::tuple{this, "# passed " + partProp->name() + " cut"} );
+      ++decay_index;
     }
-    // prepare the cuts
     return sc;
   }
 
-  IGeometryInfo const& geometry() const {
-    /// FIXME should be integrated into functional algorithm, but MultiMergingTranformer does not allow it
-    return *m_lhcb.get( getConditionContext( Gaudi::Hive::currentContext() ) ).geometry();
+  Combiner::types::FilterOutput_t operator()( InputContainer<InputTypes> const&... particles,
+                                              DetectorElement const& lhcb ) const override {
+    // Create a parameter pack for manipulating the N-1 combination cuts
+    return combine( particles..., *lhcb.geometry(), std::make_index_sequence<N - 1>{} );
   }
 
-  Combiner::types::DaughterCuts_t                               m_daughter_cuts;
-  mutable std::map<int, Gaudi::Accumulators::BinomialCounter<>> m_daughters_passed{};
-  mutable Gaudi::Accumulators::BinomialCounter<>                m_mother_passed{this, "# passed mother cut"};
+private:
+  template <std::size_t... ICombinationCuts>
+  Combiner::types::FilterOutput_t combine( InputContainer<InputTypes> const&... particles,
+                                           IGeometryInfo const& geometry,
+                                           std::index_sequence<ICombinationCuts...> ) const {
+    // Retrieve the combination and composite cuts
+    auto const& unprepared_combination_cuts = this->template getFunctors<CombCut<ICombinationCuts>...>();
+    auto const& unprepared_composite_cut    = this->template getFunctor<CompositeCut>();
+
+    // Create the output, reserving enough space for the 'typical' number of
+    // passing candidates
+    Combiner::types::Output_t out;
+    std::get<0>( out ).reserve( Combiner::constants::max_outcontainer_size );
+    std::get<1>( out ).reserve( Combiner::constants::max_outcontainer_size );
+
+    process_decays( particles..., out,
+                    std::tuple{std::get<ICombinationCuts>( unprepared_combination_cuts ).prepare()...},
+                    unprepared_composite_cut.prepare(), geometry,
+                    // Buffer counters to reduce atomic operations
+                    std::tuple{m_npassed_combination_cuts[ICombinationCuts].buffer()...}, m_npassed_vertex_fit.buffer(),
+                    m_npassed_composite_cut.buffer(),
+                    // Generate some integer packs for convenience
+                    std::index_sequence<ICombinationCuts...>{}, std::index_sequence_for<InputTypes...>{} );
+
+    const auto filter_passed = !std::get<0>( out ).empty();
+    m_npassed += filter_passed;
+    return std::tuple_cat( std::make_tuple( filter_passed ), std::move( out ) );
+  }
 
-  /// FIXME should be integrated in functional algorithm, but MergingTransformers
-  /// do not allow extra (non vector) inputs
-  LHCb::DetDesc::ConditionAccessor<DetectorElement> m_lhcb{this, "StandardGeometryTop", "/dd/Structure/LHCb",
-                                                           "Access to the LHCb geometry"};
+  template <typename CombinationCuts, typename CompositeCut, typename Counter, typename Counters,
+            std::size_t... ICombinationCuts, std::size_t... IInput>
+  void process_decays( InputContainer<InputTypes> const&... particles, Combiner::types::Output_t& storage,
+                       CombinationCuts combination_cuts, CompositeCut composite_cut, IGeometryInfo const& geometry,
+                       Counters npassed_combination_cuts, Counter npassed_vertex_fit, Counter npassed_composite_cut,
+                       std::index_sequence<ICombinationCuts...>, std::index_sequence<IInput...> ) const {
+    // Cannot use a structured binding (auto& [ps, vs] = storage) because Clang
+    // cannot compile the vertex fit lambda when it uses a variable from a
+    // structured binding (see LLVM bug #39963)
+    auto& output_particles = std::get<0>( storage );
+    auto& output_vertices  = std::get<1>( storage );
+
+    if ( UNLIKELY( this->msgLevel( MSG::DEBUG ) ) ) {
+      ( ( this->debug() << "Container #" << IInput << ": " << &( particles )
+                        << "; first element: " << *particles.begin() << endmsg ),
+        ... );
+    }
+    // Record input container sizes
+    ( ( m_ninput_particles[IInput] += particles.size() ), ... );
+
+    // Lambda that handles N-body combinations once we've got that far
+    // `indices`: indices into the input containers of the members of this
+    //            N-body combination
+    // `decay_index`: indices into `m_decays` of the decay descriptor
+    //                       this N-body combination was produced from
+    auto do_vertex_fit = [&]( std::array<int, N> const& indices, const int decay_index ) {
+      // Prepare for the vertex fit
+      const auto& pid      = m_decays[decay_index].mother().pid();
+      auto        particle = std::make_unique<Combiner::types::Particle_t>( pid );
+      auto        vertex   = std::make_unique<Combiner::types::Vertex_t>();
+      if ( UNLIKELY( this->msgLevel( MSG::DEBUG ) ) ) {
+        this->debug() << "Promoting " << N << "-body combination to composite, "
+                      << "indices = " << indices << ", decay index = " << decay_index << ", parent PID = " << pid.pid()
+                      << endmsg;
+      }
 
-  // legacy stuff
-  ToolHandle<IDistanceCalculator> m_dist_calc{this, "DistanceCalculator", "LoKi::DistanceCalculator"};
-  // get particle pids from names:
-  ServiceHandle<LHCb::IParticlePropertySvc> m_particlePropSvc{this, "ParticlePropertySvc", "LHCb::ParticlePropertySvc"};
+      // Load the children from the input containers
+      LHCb::Particle::ConstVector children{[&]( auto const& container, auto const& index ) {
+        return container[index];
+      }( particles, std::get<IInput>( indices ) )...};
+      // Run the fit
+      auto const fit_success = ( m_vertex_fitter->combine( children, *particle, *vertex, geometry ) ).isSuccess();
+      // Record the fit success statistics.
+      npassed_vertex_fit += fit_success;
+      if ( UNLIKELY( this->msgLevel( MSG::DEBUG ) ) ) {
+        this->debug() << "Vertex fit " << ( fit_success ? "passed" : "failed" ) << endmsg;
+      }
 
-  // decay descriptors
-  ToolHandle<IDecodeSimpleDecayString> m_decoder{this, "DecayDecoder", "DecodeSimpleDecayString"};
-  ToolHandle<ICheckOverlap>            m_overlapTool{this, "OverlapTool", "CheckOverlap"};
-  ToolHandle<IParticleCombiner>        m_vertexFitter{this, "ParticleCombiner", "LoKi::VertexFitter"};
+      if ( fit_success ) {
+        auto const composite_cut_success = composite_cut( *particle );
+        npassed_composite_cut += composite_cut_success;
+        if ( UNLIKELY( this->msgLevel( MSG::DEBUG ) ) ) {
+          this->debug() << "Composite cut " << ( fit_success ? "passed" : "failed" ) << endmsg;
+        }
+        if ( composite_cut_success ) {
+          output_particles.add( particle.release() );
+          output_vertices.add( vertex.release() );
+        }
+      }
+    };
 
-  Gaudi::Property<std::vector<std::string>> m_decaydescriptors{
-      this, "DecayDescriptors", {}, "Please provide a list of decay descriptors"};
-  std::vector<Decays::Decay> m_decays;
-};
+    for ( int n_decay = 0; n_decay < int( m_decays.size() ); ++n_decay ) {
+      auto const& decay = m_decays[n_decay];
+      // Collect indices of particles in the input containers that satisfy
+      // the constraints imposed by the decay descriptor. e.g. if the 0th
+      // input container is a mix of D0 and D~0 then get the indices of the
+      // D0s. If the Nth input is the same as the N-1th input then returns a
+      // null std::optional indicating that a "triangle" loop should be used
+      // TODO(AP) could hoist this outside the loop for the first iteration, and on
+      // the second iteration we can re-use the indices in the case where the
+      // PID is the same (the child is self-conjugate)
+      std::array const matching_indices_storage{get_indices<IInput>( particles..., decay )...};
+      // See if we can shortcut -- if we got zero entries for any compulsory
+      // input then we can abort processing here
+      if ( std::any_of( matching_indices_storage.begin(), matching_indices_storage.end(),
+                        []( auto const& opt_indices ) { return opt_indices.has_value() && opt_indices->empty(); } ) ) {
+        continue;
+      }
+      // Cache which list of indices we will actually use in each slot
+      std::array const matching_indices{[&]( std::size_t InputIndex ) {
+        // get_indices indicated that we can use the previous (InputIndex - 1) container
+        bool const do_triangle_loop = !matching_indices_storage[InputIndex].has_value();
+        if ( do_triangle_loop ) {
+          // Look in indices below InputIndex for a populated entry in matching_indices_storage
+          assert( InputIndex > 0 );
+          for ( std::size_t index_to_try = InputIndex - 1;; --index_to_try ) {
+            if ( matching_indices_storage[index_to_try].has_value() ) {
+              // Found the match
+              return std::pair{do_triangle_loop, std::cref( matching_indices_storage[index_to_try].value() )};
+            }
+            assert( index_to_try != 0 );
+          }
+          assert( false ); // we should never get here
+        }
+        return std::pair{do_triangle_loop, std::cref( matching_indices_storage[InputIndex].value() )};
+      }( IInput )...};
+
+      // Defer to another function to actually generate the combinations
+      process_combinations( particles..., n_decay, matching_indices, do_vertex_fit, combination_cuts,
+                            npassed_combination_cuts, geometry, std::index_sequence<ICombinationCuts...>{},
+                            std::make_index_sequence<sizeof...( ICombinationCuts ) - 1>{} );
+    }
+  }
 
-template <std::size_t N>
-class NBodyParticleCombiner;
+  template <typename ArrayOfIndexVectors, typename CombinationCuts, typename CombinationCutCounters, typename VertexFit,
+            std::size_t... ICombinationCuts, std::size_t... ICombinationCutsMinusOne>
+  void process_combinations( InputContainer<InputTypes> const&... particles, int n_decay,
+                             ArrayOfIndexVectors const& matching_indices, VertexFit& do_vertex_fit,
+                             CombinationCuts const& combination_cuts, CombinationCutCounters& npassed_combination_cuts,
+                             IGeometryInfo const& geometry, std::index_sequence<ICombinationCuts...>,
+                             std::index_sequence<ICombinationCutsMinusOne...> ) const {
+    // Indicies into the first container
+    auto const& indices0 = matching_indices[0].second.get();
+    // Indicies into the second container
+    auto const& indices1 = matching_indices[1].second.get();
+    // Whether we should do a 'triangle loop' to generate all 2-body combinations of the above indicies
+    bool const do_triangle_loop = matching_indices[1].first;
+    // The logic in get_indices should never mark the first container for the triangle loop (it has no preceeding
+    // container)
+    assert( !matching_indices[0].first );
+    if ( UNLIKELY( this->msgLevel( MSG::DEBUG ) ) ) {
+      this->debug() << "indices0.size() == " << indices0.size() << ", indices1.size() == " << indices1.size()
+                    << ", do_triangle_loop == " << do_triangle_loop << endmsg;
+    }
 
-template <>
-class NBodyParticleCombiner<2>
-    : public with_functors<NBodyParticleCombinerBase, Combiner::types::CombType<2, 12>::Cut> {
-public:
-  using base = with_functors<NBodyParticleCombinerBase, Combiner::types::CombType<2, 12>::Cut>;
-  using base::base;
+    // We will need these objects for all combinations passing the overlap cut,
+    // which is likely to be most combinations, so create them once outside the
+    // loop
+    auto const& comb_cut         = std::get<0>( combination_cuts );
+    auto&       npassed_comb_cut = std::get<0>( npassed_combination_cuts );
+    auto        combination      = particle_combination<2>( geometry );
+
+    auto const inputs = std::tuple{particles...};
+    for ( auto i0 = 0; i0 < int( indices0.size() ); ++i0 ) {
+      // Load the particles matching the first child of the current decay
+      // descriptor from the first container
+      auto const p0_index = indices0[i0];
+      auto const p0       = std::get<0>( inputs )[p0_index];
+      for ( auto i1 = do_triangle_loop ? i0 + 1 : 0; i1 < int( indices1.size() ); ++i1 ) {
+        auto const       p1_index = indices1[i1];
+        auto const       p1       = std::get<1>( inputs )[p1_index];
+        std::array const combination_indices{p0_index, p1_index};
+        if ( UNLIKELY( this->msgLevel( MSG::DEBUG ) ) ) {
+          this->debug() << "Processing 2-body combination; indices = " << combination_indices << endmsg;
+        }
+
+        auto const has_overlap = m_check_overlap->foundOverlap( p0, p1 );
+        if ( UNLIKELY( this->msgLevel( MSG::DEBUG ) ) ) {
+          this->debug() << "Overlap check " << ( has_overlap ? "failed" : "passed" ) << endmsg;
+        }
+
+        combination.m_children     = {p0, p1};
+        auto const passed_comb_cut = !has_overlap && comb_cut( combination );
+        npassed_comb_cut += passed_comb_cut;
+        if ( UNLIKELY( this->msgLevel( MSG::DEBUG ) ) ) {
+          this->debug() << "Combination12Cut " << ( passed_comb_cut ? "passed" : "failed" ) << endmsg;
+        }
+        if ( !passed_comb_cut ) { continue; }
+
+        if constexpr ( N == 2 ) {
+          // We now have a fully selected N-body combination, so run the vertex fit
+          std::invoke( do_vertex_fit, combination_indices, n_decay );
+        } else {
+          // Otherwise, we need to add another particle to these combinations
+          extend_combination( particles..., n_decay, matching_indices, combination_indices, do_vertex_fit,
+                              combination_cuts, npassed_combination_cuts, geometry, std::make_index_sequence<2>{} );
+        }
+      }
+    }
+  }
 
-  NBodyParticleCombiner( const std::string& name, ISvcLocator* pSvcLocator ) : base( name, pSvcLocator ) {}
+  /** Add one more child to an M-body combination.
+   *
+   * The `IParticle` index has the size of the M-body combinations which have
+   * already been built and survived the corresponding combination cuts.
+   *
+   * `NBodies` is the size of the incoming combinations, i.e. the B-body
+   * combinations that we will try to promote to an M+1-body combination.  If
+   * the promotion is successful and M == N the combination will enter the
+   * vertex fit, otherwise this method recursively extends the new combination
+   * until that condition is reached.
+   */
+  template <std::size_t... IParticle, typename CombinationCuts, typename VertexFit, typename CombinationCutCounters,
+            std::size_t NBodies, typename ArrayOfIndexVectors>
+  void extend_combination( InputContainer<InputTypes> const&... particles, int n_decay,
+                           ArrayOfIndexVectors const&      matching_indices,
+                           std::array<int, NBodies> const& current_combination_indices, VertexFit const& do_vertex_fit,
+                           CombinationCuts const& combination_cuts, CombinationCutCounters& npassed_combination_cuts,
+                           IGeometryInfo const& geometry, std::index_sequence<IParticle...> ) const {
+    // Combination indices are the indices into input containers of the elements
+    // of the current combination
+    // NBodies is the size of the combination we've been given and that we will
+    // try to extend
+    static_assert( NBodies == sizeof...( IParticle ) );
+    static_assert( NBodies >= 2 );
+    static_assert( NBodies < N );
+    static_assert( std::tuple_size_v<CombinationCuts> == N - 1 );
+    static_assert( std::tuple_size_v<CombinationCutCounters> == N - 1 );
+    // Size of the combination we will build in this call
+    constexpr std::size_t NBodiesExtended = NBodies + 1;
+    // Get the index of the combination cut that will act on this type
+    constexpr std::size_t CombCutIndex = NBodies + 1 - 2;
+
+    if ( UNLIKELY( this->msgLevel( MSG::DEBUG ) ) ) {
+      this->debug() << "Promoting " << NBodies << "-body combination, indices = " << current_combination_indices
+                    << endmsg;
+    }
+    bool const            do_triangle_loop = std::get<0>( matching_indices[NBodies] );
+    index_vector_t const& new_indices      = std::get<1>( matching_indices[NBodies] );
+    // If the previous container is identical (i.e. do_triangle_loop is true)
+    // then we only need to consider candidates after the last candidate in the
+    // current combination
+    int const first_index = [&]() -> int {
+      if ( do_triangle_loop ) {
+        return std::distance( new_indices.begin(), std::upper_bound( new_indices.begin(), new_indices.end(),
+                                                                     current_combination_indices.back() ) );
+      } else {
+        return 0;
+      }
+    }();
+    if ( UNLIKELY( this->msgLevel( MSG::DEBUG ) ) ) {
+      this->debug() << "indices" << NBodies << ".size() = " << new_indices.size()
+                    << ", do_triangle_loop = " << do_triangle_loop << ", starting index = " << first_index << endmsg;
+    }
 
-  Combiner::types::FilterOut_t operator()( Combiner::types::In_t const& particles ) const override;
+    // We will need these objects for all combinations passing the overlap cut,
+    // which is likely to be most combinations, so create them once outside the
+    // loop
+    auto const& comb_cut         = std::get<CombCutIndex>( combination_cuts );
+    auto&       npassed_comb_cut = std::get<CombCutIndex>( npassed_combination_cuts );
+    auto        combination      = particle_combination<NBodiesExtended>( geometry );
+
+    // Loop over the new additions to the input combinations
+    auto const inputs = std::tuple{particles...};
+    for ( auto i = first_index; i < int( new_indices.size() ); ++i ) {
+      auto const                             pN1_index = new_indices[i];
+      auto const                             pN1       = std::get<NBodies>( inputs )[pN1_index];
+      std::array<int, NBodiesExtended> const extended_combination_indices{
+          []( auto const index ) { return index; }( std::get<IParticle>( current_combination_indices ) )..., pN1_index};
+      if ( UNLIKELY( this->msgLevel( MSG::DEBUG ) ) ) {
+        this->debug() << "Processing " << NBodiesExtended
+                      << "-body combination; indices = " << extended_combination_indices << endmsg;
+      }
+      // Update the combination children with the current set of particles
+      combination.m_children = {
+          ( std::get<IParticle>( inputs )[std::get<IParticle>( extended_combination_indices )] )..., pN1};
+
+      // We only need to check the overlap between each child in the current
+      // N-body object and the current particle (i.e. we don't need to check
+      // the overlap *within* the N-body object; that has already been done)
+      auto const has_overlap =
+          ( m_check_overlap->foundOverlap( std::get<IParticle>( combination.m_children ), pN1 ) || ... );
+      if ( UNLIKELY( this->msgLevel( MSG::DEBUG ) ) ) {
+        this->debug() << "Overlap check " << ( has_overlap ? "failed" : "passed" ) << endmsg;
+      }
 
-protected:
-  mutable Gaudi::Accumulators::BinomialCounter<> m_combination12_passed{this, "# passed comb12 cut"};
-};
+      auto const passed_comb_cut = !has_overlap && comb_cut( combination );
+      npassed_comb_cut += passed_comb_cut;
+      if ( UNLIKELY( this->msgLevel( MSG::DEBUG ) ) ) {
+        this->debug() << Combiner::combinationCutName<N, NBodies - 2>() << " cut "
+                      << ( passed_comb_cut ? "passed" : "failed" ) << endmsg;
+      }
+      if ( !passed_comb_cut ) { continue; }
 
-template <>
-class NBodyParticleCombiner<3>
-    : public with_functors<NBodyParticleCombiner<2>, typename Combiner::types::CombType<3, 123>::Cut,
-                           typename Combiner::types::CombType<2, 13>::Cut> {
-public:
-  using base = with_functors<NBodyParticleCombiner<2>, typename Combiner::types::CombType<3, 123>::Cut,
-                             typename Combiner::types::CombType<2, 13>::Cut>;
-  using base::base;
+      if constexpr ( NBodiesExtended == N ) {
+        // We now have a fully selected N-body combination, so run the vertex fit
+        std::invoke( do_vertex_fit, extended_combination_indices, n_decay );
+      } else {
+        // We need to add another particle to these combinations
+        extend_combination( particles..., n_decay, matching_indices, extended_combination_indices, do_vertex_fit,
+                            combination_cuts, npassed_combination_cuts, geometry,
+                            std::make_index_sequence<NBodiesExtended>{} );
+      }
+    }
+  }
 
-  Combiner::types::FilterOut_t operator()( Combiner::types::In_t const& particles ) const override;
+  // Collect relevant indices from the input containers, based on pid+charge as defined in the decay descriptor
+  template <std::size_t IInput>
+  auto get_indices( InputContainer<InputTypes> const&... particles, Decays::Decay const& decay ) const {
+    using IterableInputTuple             = std::tuple<InputContainer<InputTypes>...>;
+    auto const                    inputs = std::tuple{particles...};
+    std::optional<index_vector_t> opt_indices{std::nullopt};
+    auto const&                   input = std::get<IInput>( inputs );
+    if constexpr ( IInput > 0 ) {
+      // Not the first input, we need to figure out if this might be one of a
+      // contiguous series of decay descriptor children of the same type. In
+      // that case we have to generate a different kind of loop, and we have
+      // to check that we were configured with a matching sequence of
+      // identical inputs (mild EWWW :()
+      if constexpr ( std::is_same_v<std::tuple_element_t<IInput - 1, IterableInputTuple>,
+                                    std::tuple_element_t<IInput, IterableInputTuple>> ) {
+        if ( decay.children()[IInput].pid() == decay.children()[IInput - 1].pid() ) {
+          // Yes, in this case we should insist that the `IInput`th and
+          // `IInput-1`th inputs were, in fact, the same...
+          // FIXME(AP): we cannot compare containers because the inputs are Range
+          // objects which get instantiated when this algorithm's operator() is
+          // called; the same underlying container passed in N times will result
+          // in N Range objects with different memory addresses. Instead we use
+          // the memory address of the first element of the container; this is
+          // very shaky and we could come up with a better definition of
+          // 'identical containers'
+          if ( UNLIKELY( *( std::get<IInput>( inputs ).begin() ) != *( std::get<IInput - 1>( inputs ).begin() ) ) ) {
+            throw GaudiException{"Got contiguous children with the same PIDs, but the corresponding algorithm inputs "
+                                 "were not the same. This is unsupported.",
+                                 "NBodyCombiner<name>", StatusCode::FAILURE};
+          }
+          // In this case there is no need to come up with a new list of
+          // indices, we will make a "triangular" (need a better name) loop
+          // using the indices from IInput - 1.
+          if ( UNLIKELY( this->msgLevel( MSG::DEBUG ) ) ) {
+            this->debug() << "Not looking at container " << &input << " as the previous container is identical"
+                          << endmsg;
+          }
+          return opt_indices;
+        }
+      }
+    }
+    // We didn't fall into the special case, so we just loop through and
+    // store indices. In this case we need an actual vector to return them in
+    opt_indices.emplace();
+    opt_indices->reserve( input.size() );
+    auto const target_pid = decay.children()[IInput].pid().pid();
+    auto       index      = 0;
+    if ( UNLIKELY( this->msgLevel( MSG::DEBUG ) ) ) {
+      this->debug() << "Looking at container = " << &input << ", size = " << input.size() << " for PID = " << target_pid
+                    << endmsg;
+    }
+    for ( auto const& particle : input ) {
+      if ( UNLIKELY( this->msgLevel( MSG::VERBOSE ) ) ) {
+        this->debug() << "Looking at container " << &input << " element " << &particle << " index " << index << "... ";
+      }
+      if ( particle->particleID().pid() == target_pid ) {
+        opt_indices->push_back( index );
+        if ( UNLIKELY( this->msgLevel( MSG::VERBOSE ) ) ) { this->debug() << "Matched" << endmsg; }
+      } else {
+        if ( UNLIKELY( this->msgLevel( MSG::VERBOSE ) ) ) { this->debug() << "Did not match" << endmsg; }
+      }
+      ++index;
+    }
+    if ( UNLIKELY( this->msgLevel( MSG::DEBUG ) ) ) {
+      this->debug() << "Container " << &input << " has " << opt_indices->size() << " matches: " << *opt_indices
+                    << endmsg;
+    }
+    return opt_indices;
+  }
 
-protected:
-  mutable Gaudi::Accumulators::BinomialCounter<> m_combination123_passed{this, "# passed comb123 cut"};
-  mutable Gaudi::Accumulators::BinomialCounter<> m_combination13_passed{this, "# passed comb13 cut"};
-};
+  /// Retrieve a ParticleCombination object suitable for an M-body combination
+  template <std::size_t M>
+  auto particle_combination( IGeometryInfo const& geometry ) const {
+    using ParticleSubCombination = typename Combiner::CombTypes<InputTypes...>::template Combination<M>;
+    return ParticleSubCombination{this->m_distance_calculator.get(), geometry};
+  }
 
-template <>
-class NBodyParticleCombiner<4>
-    : public with_functors<NBodyParticleCombiner<3>, typename Combiner::types::CombType<4, 1234>::Cut,
-                           typename Combiner::types::CombType<2, 14>::Cut> {
-public:
-  using base = with_functors<NBodyParticleCombiner<3>, typename Combiner::types::CombType<4, 1234>::Cut,
-                             typename Combiner::types::CombType<2, 14>::Cut>;
-  using base::base;
+  template <typename C, std::size_t... Ms>
+  static std::array<C, sizeof...( Ms )> make_input_particle_counters( NBodyParticleCombiner* ptr,
+                                                                      std::index_sequence<Ms...> ) {
+    // No space in InputN to match the property name
+    return {C{ptr, "Input" + std::to_string( Ms + 1 ) + " size"}...};
+  }
+
+  template <typename C, std::size_t... Ms>
+  static std::array<C, sizeof...( Ms )> make_combination_cut_counters( NBodyParticleCombiner* ptr,
+                                                                       std::index_sequence<Ms...> ) {
+    return {C{ptr, "# passed " + Combiner::combinationCutName<N, Ms>()}...};
+  }
 
-  Combiner::types::FilterOut_t operator()( Combiner::types::In_t const& particles ) const override;
+  /// Distance calculator passed to the ParticleCombination instance
+  ToolHandle<IDistanceCalculator>      m_distance_calculator{this, "DistanceCalculator", "LoKi::DistanceCalculator"};
+  ToolHandle<IDecodeSimpleDecayString> m_decay_descriptor_decoder{this, "DecayDecoder", "DecodeSimpleDecayString"};
+  ToolHandle<ICheckOverlap>            m_check_overlap{this, "OverlapTool", "CheckOverlap"};
+  ToolHandle<IParticleCombiner>        m_vertex_fitter{this, "ParticleCombiner", "ParticleVertexFitter"};
+
+  // get particle pids from names:
+  ServiceHandle<LHCb::IParticlePropertySvc> m_particlePropSvc{this, "ParticlePropertySvc", "LHCb::ParticlePropertySvc"};
+
+  Gaudi::Property<std::string> m_decay_descriptor{this, "DecayDescriptor", {}, "The decay topology to reconstruct."};
+  /// Store decay objects parsed from the descriptor (this will be a maximum of
+  /// two, in the case of a `[]cc` decay descriptor)
+  std::vector<Decays::Decay> m_decays;
 
-protected:
-  mutable Gaudi::Accumulators::BinomialCounter<> m_combination1234_passed{this, "# passed comb1234 cut"};
-  mutable Gaudi::Accumulators::BinomialCounter<> m_combination14_passed{this, "# passed comb14 cut"};
+  using counter_t = Gaudi::Accumulators::BinomialCounter<>;
+  /// Number of objects in each input container
+  mutable std::array<Gaudi::Accumulators::SummingCounter<unsigned int>, N> m_ninput_particles{
+      make_input_particle_counters<Gaudi::Accumulators::SummingCounter<unsigned int>>( this,
+                                                                                       std::make_index_sequence<N>{} )};
+  /// Number of combinations passing the combination cut; note that the overlap
+  /// check *is counted* as a 'combination cut'
+  mutable std::array<counter_t, N - 1> m_npassed_combination_cuts{
+      make_combination_cut_counters<counter_t>( this, std::make_index_sequence<N - 1>{} )};
+  /// Number of successfully fitted combinations
+  mutable counter_t m_npassed_vertex_fit{this, "# passed vertex fit"};
+  /// Number of post-vertex-fit particles passing the composite cut
+  mutable counter_t m_npassed_composite_cut{this, "# passed CompositeCut"};
+  /// Number of operator() calls with at least one candidate produced
+  mutable counter_t m_npassed{this, "# passed"};
 };
diff --git a/Phys/ParticleCombiners/src/ParticleCombiner.cpp b/Phys/ParticleCombiners/src/ParticleCombiner.cpp
index eab6a86e425a1023e2643fcc9715d1e4dc38cf56..fc2e2620e413d74c82271c8ae5fe618e98fae184 100644
--- a/Phys/ParticleCombiners/src/ParticleCombiner.cpp
+++ b/Phys/ParticleCombiners/src/ParticleCombiner.cpp
@@ -1,5 +1,5 @@
 /*****************************************************************************\
-* (c) Copyright 2000-2019 CERN for the benefit of the LHCb Collaboration      *
+* (c) Copyright 2000-2021 CERN for the benefit of the LHCb Collaboration      *
 *                                                                             *
 * This software is distributed under the terms of the GNU General Public      *
 * Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".   *
@@ -8,437 +8,21 @@
 * granted to it by virtue of its status as an Intergovernmental Organization  *
 * or submit itself to any jurisdiction.                                       *
 \*****************************************************************************/
-
-#include <functional>
-
-#include <boost/algorithm/string/join.hpp>
-
-#include "CombKernel/ParticleCombination.h"
 #include "CombKernel/ParticleCombiner.h"
-#include "LoKi/Combiner.h"
-#include "LoKi/ParticleProperties.h"
-#include "PrKernel/PrSelection.h"
-#include "TrackKernel/TrackStateVertex.h"
-#include "TrackKernel/TrackVertexUtils.h"
-
-/*
-TODO features:
-  1. Functors
-  4. deduplicate children
-
-TODO optimizations:
-  1. the overlap tool asks about the several properties every time that are the same for entire containers (like
-isBasicParticle)
-  2. we make too many comparisons with the overlaptool, it would sometimes be better to call it with more than 2
-particles
-  2. the default combination cut is ALL, in which case selecting through entire containers just to get n* true is
-useless.
-  3. vectorize cuts/selection/fits
-  4. check at what point in the loop we should invoke cuts (for instance: in 3 body combinations, upfront comb13 cut or
-in the inner loop?)
-  5. improve the vertex fit in docachi2
-  6. if p2 and p3 are of the same type and the cut12/cut13 are the same, the selections need not be executed twice!
-*/
-
-DECLARE_COMPONENT_WITH_ID( NBodyParticleCombiner<2>, "TwoBodyCombiner" )
-DECLARE_COMPONENT_WITH_ID( NBodyParticleCombiner<3>, "ThreeBodyCombiner" )
-DECLARE_COMPONENT_WITH_ID( NBodyParticleCombiner<4>, "FourBodyCombiner" )
-
-using namespace Combiner::types;
-
-namespace {
-
-  auto prepare_daughter_cut = []( int pid, Combiner::types::ParticleCut_t const& cut, auto& ctr ) {
-    return [dau_cut = cut.prepare(), _ctr = ctr.buffer(), pid]( auto const* p ) mutable {
-      if ( not p ) return false; // TODO hopefully to be removed?
-      if ( p->particleID().pid() != pid ) return false;
-      bool passed = dau_cut( *p );
-      _ctr += passed;
-      return passed;
-    };
-  };
-
-  auto prepare_mother_cut = []( Combiner::types::ParticleCut_t const& cut, auto& ctr ) {
-    return [mother_cut = cut.prepare(), _ctr = ctr.buffer()]( auto const& p ) mutable {
-      bool passed = mother_cut( p );
-      _ctr += passed;
-      return passed;
-    };
-  };
 
-  auto prepare_combination_cut = []( auto& comb, auto const& cut, auto& ctr ) {
-    return [comb_cut = cut.prepare(), _ctr = ctr.buffer(), &comb]( auto const*... ps ) mutable {
-      comb.m_children = {ps...};
-      bool passed     = comb_cut( comb );
-      _ctr += passed;
-      return passed;
-    };
+// TODO this should live in Event/Particle.h, otherwise every ThOr functor
+// algorithm using LHCb::Particle has to define this header map themselves
+namespace LHCb {
+  template <>
+  struct header_map<LHCb::Particle> {
+    constexpr static auto value = LHCb::string_array{"Event/Particle.h"};
   };
+} // namespace LHCb
 
-  template <typename Pred_t>
-  ParticleContainer_t select_particles( ParticleContainer_t const& pcontainer, Pred_t&& pred ) {
-    ParticleContainer_t out{};
-    // out.reserve( pcontainer.size() ); //need this only in case of normal vectors
-    // TODO vectorize this pls, requires more soa like particle container
-    std::copy_if( pcontainer.begin(), pcontainer.end(), std::back_inserter( out ), std::forward<Pred_t>( pred ) );
-    return out;
-  }
-
-  auto select_daughters = []( auto const& particles, auto const& daughter_cuts, auto& cntrs ) {
-    std::map<int, ParticleContainer_t> selected_daughters{};
-    for ( auto const& daughter_cut : daughter_cuts ) {
-      if ( daughter_cut.second == nullptr ) {
-        selected_daughters.emplace( daughter_cut.first, select_particles( particles, [&]( Particle_ptr_t p ) {
-                                      if ( not p ) return false; // TODO hopefully to be removed?
-                                      bool passed = p->particleID().pid() == daughter_cut.first;
-                                      cntrs[daughter_cut.first] += passed;
-                                      return passed;
-                                    } ) );
-      } else {
-        selected_daughters.emplace(
-            daughter_cut.first,
-            select_particles( particles, prepare_daughter_cut( daughter_cut.first, *daughter_cut.second,
-                                                               cntrs[daughter_cut.first] ) ) );
-      }
-    }
-    return selected_daughters;
-  };
-
-  auto fill_candidates = []( auto& candidates, auto const& decay, auto const& selected_daughters ) {
-    const Decays::Decay::Items& items = decay.children();
-    for ( auto i = 0u; i < candidates.size(); ++i ) {
-      auto it_to_selected_daughter = selected_daughters.find( items[i].pid().pid() );
-      if ( it_to_selected_daughter != selected_daughters.end() ) {
-        candidates[i] = &it_to_selected_daughter->second;
-      } else {
-        throw std::runtime_error( "no entry in the daughter cut list for PID " +
-                                  std::to_string( items[i].pid().pid() ) );
-      }
-    }
-  };
-
-  template <std::size_t n>
-  using comb_t = std::array<Particle_ptr_t, n>;
-
-  // check order if both particles are of same type
-  template <bool same_type = false>
-  bool check_order( Particle_ptr_t p1, Particle_ptr_t p2 ) {
-    if constexpr ( not same_type ) {
-      return true;
-    } else {
-      // if the pid is the same, we merged all and thus only accept 1 order to "deduplicate"
-      auto const& p1mom = p1->momentum();
-      auto const& p2mom = p2->momentum();
-      if ( p1mom.perp2() < p2mom.perp2() ) {
-        return true;
-      } else if ( p2mom.perp2() < p1mom.perp2() ) {
-        return false;
-      } // RETURN
-      // compare by E
-      const double e1 = p1mom.E();
-      const double e2 = p2mom.E();
-      return e1 < e2;
-    }
-  }
-
-  // transform a combination cut into a cut that only takes one particle as input (the other one is captured
-  template <typename Cut, typename O>
-  auto make_cut_predicate( Cut&& cut, O&& ordered, ICheckOverlap const* overlapCheck,
-                           Combiner::types::Particle_ptr_t p1 ) { // a predicate for selection
-    return [p1, &cut, ordered = std::forward<O>( ordered ), overlapCheck]( Particle_ptr_t p2 ) -> bool {
-      return ordered( p1, p2 ) and not overlapCheck->foundOverlap( p1, p2 ) and cut( p1, p2 );
-    };
-  }
-
-  // we compare pids from every particle in the container, to see whether we have to order in the combination
-  //(see check_order)
-  template <typename T, std::size_t n>
-  auto compare_pids( std::array<T, n> const& container ) {
-    auto                              get_pid = []( auto const& i ) { return i[0]->particleID().pid(); };
-    std::array<bool, n*( n - 1 ) / 2> ret;
-    int                               idx = 0;
-    for ( auto i = 0u; i < n; ++i ) {
-      for ( auto j = i + 1; j < n; ++j ) { ret[idx++] = get_pid( *container[i] ) == get_pid( *container[j] ); }
-    }
-    return ret;
-  }
-
-  template <typename Out_t, typename Candidates_t, typename cut12_t>
-  Out_t combine_2( Candidates_t const& candidates, bool const pid_comparison, ICheckOverlap const* overlapChecker,
-                   cut12_t&& cut12 ) {
-    Out_t combinations;
-    auto  is_ordered12 = pid_comparison ? check_order<true> : check_order<false>;
-
-    for ( Particle_ptr_t p1 : *candidates[0] ) {
-      for ( Particle_ptr_t p2 :
-            select_particles( *candidates[1], make_cut_predicate( cut12, is_ordered12, overlapChecker, p1 ) ) ) {
-        combinations.emplace_back( comb_t<2>{p1, p2} );
-      }
-    }
-    return combinations;
-  }
-
-  template <typename Out_t, typename Candidates_t, typename cut12_t, typename cut13_t, typename cut123_t,
-            typename pid_comparison_t>
-  Out_t combine_3( Candidates_t const& candidates, pid_comparison_t const& pid_comparison,
-                   ICheckOverlap const* overlapChecker, cut12_t&& cut12, cut13_t&& cut13, cut123_t&& cut123 ) {
-    Out_t combinations;
-    auto  is_ordered12 = pid_comparison[0] ? check_order<true> : check_order<false>;
-    auto  is_ordered13 = pid_comparison[1] ? check_order<true> : check_order<false>;
-    auto  is_ordered23 = pid_comparison[2] ? check_order<true> : check_order<false>;
-
-    for ( Particle_ptr_t p1 : *( candidates[0] ) ) {
-      auto p2s = select_particles( *( candidates[1] ), make_cut_predicate( cut12, is_ordered12, overlapChecker, p1 ) );
-      if ( p2s.empty() ) continue;
-      auto p3s = select_particles( *( candidates[2] ), make_cut_predicate( cut13, is_ordered13, overlapChecker, p1 ) );
-      if ( p3s.empty() ) continue;
-      for ( Particle_ptr_t p2 : p2s ) {
-        for ( Particle_ptr_t p3 : p3s ) {
-          if ( is_ordered23( p2, p3 ) and not overlapChecker->foundOverlap( p2, p3 ) and cut123( p1, p2, p3 ) ) {
-            combinations.emplace_back( comb_t<3>{p1, p2, p3} );
-          }
-        }
-      }
-    }
-    return combinations;
-  }
-
-  template <typename Out_t, typename Candidates_t, typename cut12_t, typename cut13_t, typename cut14_t,
-            typename cut123_t, typename cut1234_t, typename pid_comparison_t>
-  Out_t combine_4( Candidates_t const& candidates, pid_comparison_t const& pid_comparison,
-                   ICheckOverlap const* overlapChecker, cut12_t&& cut12, cut13_t&& cut13, cut14_t&& cut14,
-                   cut123_t&& cut123, cut1234_t&& cut1234 ) {
-    Out_t combinations{};
-    auto  is_ordered12 = pid_comparison[0] ? check_order<true> : check_order<false>;
-    auto  is_ordered13 = pid_comparison[1] ? check_order<true> : check_order<false>;
-    auto  is_ordered14 = pid_comparison[2] ? check_order<true> : check_order<false>;
-    auto  is_ordered23 = pid_comparison[3] ? check_order<true> : check_order<false>;
-    auto  is_ordered24 = pid_comparison[4] ? check_order<true> : check_order<false>;
-    auto  is_ordered34 = pid_comparison[5] ? check_order<true> : check_order<false>;
-
-    for ( Particle_ptr_t p1 : *( candidates[0] ) ) {
-      auto p2s = select_particles( *( candidates[1] ), make_cut_predicate( cut12, is_ordered12, overlapChecker, p1 ) );
-      if ( p2s.empty() ) continue;
-      auto p3s = select_particles( *( candidates[2] ), make_cut_predicate( cut13, is_ordered13, overlapChecker, p1 ) );
-      if ( p3s.empty() ) continue;
-      auto p4s = select_particles( *( candidates[3] ), make_cut_predicate( cut14, is_ordered14, overlapChecker, p1 ) );
-      if ( p4s.empty() ) continue;
-      for ( Particle_ptr_t p2 : p2s ) {
-        for ( Particle_ptr_t p3 : p3s ) {
-          if ( is_ordered23( p2, p3 ) and cut123( p1, p2, p3 ) ) {
-            for ( Particle_ptr_t p4 : p4s ) {
-              if ( is_ordered24( p2, p4 ) and is_ordered34( p3, p4 ) and
-                   not overlapChecker->foundOverlap( p2, p3, p4 ) and cut1234( p1, p2, p3, p4 ) ) {
-                combinations.emplace_back( comb_t<4>{p1, p2, p3, p4} );
-              }
-            }
-          }
-        }
-      }
-    }
-    return combinations;
-  }
-
-  template <typename Cut_t, typename Comb_t>
-  void fit_and_append_mothers( Combiner::types::Out_t& ret, const LHCb::ParticleID& pid, const Comb_t& combinations,
-                               const IParticleCombiner* combiner, Cut_t&& mother_cut, IGeometryInfo const& geometry ) {
-    auto& [mothers, vertices] = ret;
-    for ( auto const& combination : combinations ) {
-      auto mother = new LHCb::Particle( pid );
-      auto vertex = new LHCb::Vertex();
-      // use default particle combiner to create the composed particle
-      // TODO remove this hack, and in general a better combiner
-      LHCb::Particle::ConstVector c;
-      c.reserve( combination.size() );
-      std::transform( combination.begin(), combination.end(), std::back_inserter( c ),
-                      []( auto const* p ) { return p; } );
-      if ( ( combiner->combine( c, *mother, *vertex, geometry ) ).isFailure() or not mother_cut( *mother ) ) {
-        delete mother; // TODO remove this ugly design, rather only build the mother if we mothercut passes
-        delete vertex;
-      } else {
-        // TODO relate bestpv // legacy: const LHCb::Particle* particle = this->markTree( &mother ) ;
-        mothers.add( mother );
-        vertices.add( vertex );
-      }
-    }
-  }
-
-  /** @brief Return true if a combiner producing `out` should set its filter to 'passed'.
-   */
-  bool filter_passed( const Out_t& out ) {
-    // Filter passes if at least one candidate passes all cuts
-    return !std::get<0>( out ).empty();
-  }
-
-} // end of namespace
-
-Combiner::types::FilterOut_t NBodyParticleCombiner<2>::
-                             operator()( Combiner::types::In_t const& all_particles_soc ) const {
-
-  // create the output type (particles and vertices)
-  Out_t out;
-  std::get<0>( out ).reserve( Combiner::constants::max_outcontainer_size );
-  std::get<1>( out ).reserve( Combiner::constants::max_outcontainer_size );
-
-  using namespace Combiner::types;
-  constexpr int nbodies = 2;
-
-  ParticleContainer_t particles{};
-  auto                all_sizes = 0u;
-  for ( auto const& container : all_particles_soc ) all_sizes += container.size();
-  particles.reserve( all_sizes );
-
-  for ( auto const& container : all_particles_soc )
-    for ( Particle_ptr_t particle : container ) particles.push_back( particle );
-
-  auto selected_daughters = select_daughters( particles, m_daughter_cuts, m_daughters_passed );
-
-  Combination2_t two_combination{this}; // this not reinstantiated every time, only children are
-                                        // changed
-  auto cut12 = prepare_combination_cut( two_combination, this->getFunctor<Combiner::types::CombType<2, 12>::Cut>(),
-                                        m_combination12_passed );
-
-  std::array<ParticleContainer_t const*, nbodies> candidates;
-
-  for ( auto const& decay : m_decays ) {
-    fill_candidates( candidates, decay, selected_daughters );
-
-    if ( std::any_of( candidates.begin(), candidates.end(), []( auto const* c ) { return c->empty(); } ) ) continue;
-
-    auto same_pid = compare_pids( candidates );
-
-    auto combinations =
-        combine_2<boost::container::small_vector<comb_t<nbodies>, Combiner::constants::max_combinationcontainer_size>>(
-            candidates, same_pid[0], m_overlapTool.get(), cut12 );
-    fit_and_append_mothers(
-        out, decay.mother().pid(), combinations, m_vertexFitter.get(),
-        prepare_mother_cut( static_cast<base::base const*>( this )->getFunctor<Combiner::types::MotherType::Cut>(),
-                            m_mother_passed ),
-        geometry() );
-  }
-
-  return std::tuple_cat( std::make_tuple( filter_passed( out ) ), std::move( out ) );
-}
-
-Combiner::types::FilterOut_t NBodyParticleCombiner<3>::
-                             operator()( Combiner::types::In_t const& all_particles_soc ) const {
-
-  // create the output type (particles and vertices)
-  Out_t out;
-  std::get<0>( out ).reserve( Combiner::constants::max_outcontainer_size );
-  std::get<1>( out ).reserve( Combiner::constants::max_outcontainer_size );
-
-  using namespace Combiner::types;
-  constexpr int nbodies = 3;
-
-  ParticleContainer_t particles{};
-  auto                all_sizes = 0u;
-  for ( auto const& container : all_particles_soc ) all_sizes += container.size();
-  particles.reserve( all_sizes );
-
-  for ( Combiner::types::InputParticleContainer_t const& container : all_particles_soc )
-    for ( Particle_ptr_t particle : container ) particles.push_back( particle );
-
-  auto selected_daughters = select_daughters( particles, m_daughter_cuts, m_daughters_passed );
-
-  Combination3_t three_combination{this}; // this not reinstantiated every time, only children are
-                                          // changed
-  Combination2_t two_combination{this};   // this not reinstantiated every time, only children are
-                                          // changed
-  auto cut12 = prepare_combination_cut(
-      two_combination, static_cast<base::base const*>( this )->getFunctor<Combiner::types::CombType<2, 12>::Cut>(),
-      m_combination12_passed );
-  auto cut13  = prepare_combination_cut( two_combination, this->getFunctor<Combiner::types::CombType<2, 13>::Cut>(),
-                                        m_combination13_passed );
-  auto cut123 = prepare_combination_cut( three_combination, this->getFunctor<Combiner::types::CombType<3, 123>::Cut>(),
-                                         m_combination123_passed );
-
-  std::array<ParticleContainer_t const*, nbodies> candidates;
-
-  for ( auto const& decay : m_decays ) {
-    fill_candidates( candidates, decay, selected_daughters );
-
-    if ( std::any_of( candidates.begin(), candidates.end(), []( auto const* c ) { return c->empty(); } ) ) continue;
-
-    auto pid_comparison = compare_pids( candidates );
-
-    auto combinations =
-        combine_3<boost::container::small_vector<comb_t<nbodies>, Combiner::constants::max_combinationcontainer_size>>(
-            candidates, pid_comparison, m_overlapTool.get(), cut12, cut13, cut123 );
-
-    fit_and_append_mothers(
-        out, decay.mother().pid(), combinations, m_vertexFitter.get(),
-        prepare_mother_cut(
-            static_cast<base::base::base const*>( this )->getFunctor<Combiner::types::MotherType::Cut>(),
-            m_mother_passed ),
-        geometry() );
-  }
-
-  return std::tuple_cat( std::make_tuple( filter_passed( out ) ), std::move( out ) );
-}
-
-Combiner::types::FilterOut_t NBodyParticleCombiner<4>::
-                             operator()( Combiner::types::In_t const& all_particles_soc ) const {
-
-  // create the output type (particles and vertices)
-  Out_t out;
-  std::get<0>( out ).reserve( Combiner::constants::max_outcontainer_size );
-  std::get<1>( out ).reserve( Combiner::constants::max_outcontainer_size );
-
-  using namespace Combiner::types;
-  constexpr int nbodies = 4;
-
-  ParticleContainer_t particles{};
-  auto                all_sizes = 0u;
-  for ( auto const& container : all_particles_soc ) all_sizes += container.size();
-  particles.reserve( all_sizes );
-
-  for ( auto const& container : all_particles_soc )
-    for ( Particle_ptr_t particle : container ) particles.push_back( particle );
-
-  auto selected_daughters = select_daughters( particles, m_daughter_cuts, m_daughters_passed );
-
-  Combination4_t four_combination{this};  // this not reinstantiated every time, only children are
-                                          // changed
-  Combination3_t three_combination{this}; // this not reinstantiated every time, only children are
-                                          // changed
-  Combination2_t two_combination{this};   // this not reinstantiated every time, only children are
-                                          // changed
-  auto cut12 = prepare_combination_cut(
-      two_combination,
-      static_cast<base::base::base const*>( this )->getFunctor<Combiner::types::CombType<2, 12>::Cut>(),
-      m_combination12_passed );
-  auto cut13 = prepare_combination_cut(
-      two_combination, static_cast<base::base const*>( this )->getFunctor<Combiner::types::CombType<2, 13>::Cut>(),
-      m_combination13_passed );
-  auto cut14  = prepare_combination_cut( two_combination, this->getFunctor<Combiner::types::CombType<2, 14>::Cut>(),
-                                        m_combination14_passed );
-  auto cut123 = prepare_combination_cut(
-      three_combination, static_cast<base::base const*>( this )->getFunctor<Combiner::types::CombType<3, 123>::Cut>(),
-      m_combination123_passed );
-  auto cut1234 = prepare_combination_cut( four_combination, this->getFunctor<Combiner::types::CombType<4, 1234>::Cut>(),
-                                          m_combination1234_passed );
-
-  std::array<ParticleContainer_t const*, nbodies> candidates;
-
-  for ( auto const& decay : m_decays ) {
-    fill_candidates( candidates, decay, selected_daughters );
-
-    if ( std::any_of( candidates.begin(), candidates.end(), []( auto const* c ) { return c->empty(); } ) ) continue;
-
-    auto pid_comparison = compare_pids( candidates );
-
-    auto combinations =
-        combine_4<boost::container::small_vector<comb_t<nbodies>, Combiner::constants::max_combinationcontainer_size>>(
-            candidates, pid_comparison, m_overlapTool.get(), cut12, cut13, cut14, cut123, cut1234 );
-
-    fit_and_append_mothers(
-        out, decay.mother().pid(), combinations, m_vertexFitter.get(),
-        prepare_mother_cut(
-            static_cast<base::base::base::base const*>( this )->getFunctor<Combiner::types::MotherType::Cut>(),
-            m_mother_passed ),
-        geometry() );
-  }
+using TwoBodyCombiner   = NBodyParticleCombiner<LHCb::Particle, LHCb::Particle>;
+using ThreeBodyCombiner = NBodyParticleCombiner<LHCb::Particle, LHCb::Particle, LHCb::Particle>;
+using FourBodyCombiner  = NBodyParticleCombiner<LHCb::Particle, LHCb::Particle, LHCb::Particle, LHCb::Particle>;
 
-  return std::tuple_cat( std::make_tuple( filter_passed( out ) ), std::move( out ) );
-}
+DECLARE_COMPONENT_WITH_ID( TwoBodyCombiner, "TwoBodyCombiner" )
+DECLARE_COMPONENT_WITH_ID( ThreeBodyCombiner, "ThreeBodyCombiner" )
+DECLARE_COMPONENT_WITH_ID( FourBodyCombiner, "FourBodyCombiner" )
diff --git a/Phys/ParticleCombiners/src/ThOrCombiner.h b/Phys/ParticleCombiners/src/ThOrCombiner.h
index 1e6e08190961477db3db53fa64dfc3ce2b1f2b53..83c5d73aeb180d159e358747773a3489da314f7d 100644
--- a/Phys/ParticleCombiners/src/ThOrCombiner.h
+++ b/Phys/ParticleCombiners/src/ThOrCombiner.h
@@ -121,13 +121,13 @@ namespace ThOr::detail::Combiner {
   };
 
   template <SIMDWrapper::InstructionSet Backend = SIMDWrapper::Scalar>
-  struct VertexCut {
+  struct CompositeCut {
     using particle_t = LHCb::v2::Event::ZipCompat::proxy_t<Backend, LHCb::Pr::ProxyBehaviour::Contiguous, Output_t>;
     using Signature  = mask_b_v<Backend>( Functors::mask_arg_t, mask_b_v<Backend> const&, particle_t const& );
     constexpr static bool is_simd = Backend != SIMDWrapper::Scalar;
     constexpr static auto Compilation =
         is_simd ? Functors::IFactory::QuietCacheOnly : Functors::IFactory::DefaultCompilationBehaviour;
-    constexpr static auto PropertyName = "VertexCut";
+    constexpr static auto PropertyName = "CompositeCut";
     constexpr static auto ExtraHeaders = LHCb::header_map_v<particle_t>;
   };
 
@@ -146,16 +146,16 @@ namespace ThOr::detail::Combiner {
     template <typename M, typename Backend>
     using CombCut = typename Comb<Backend::value, Ts...>::template Cut<M::value>;
     template <typename Backend>
-    using VtxCut = VertexCut<Backend::value>;
+    using CompCut = CompositeCut<Backend::value>;
     // Make with_functors<FilterTransform<Ts...>,
     //                    [cartesian expansion of CombCut<A, B>...],
-    //                    VtxCut<B>...>
+    //                    CompCut<B>...>
     // Where A is an integer from 0 to N-2 and B is the set of SIMD backends
     using type = boost::mp11::mp_append<
         with_functors<FilterTransform<Ts...>>,
         boost::mp11::mp_product<CombCut, boost::mp11::mp_from_sequence<std::make_index_sequence<sizeof...( Ts ) - 1>>,
                                 backend_ts>,
-        boost::mp11::mp_product<VtxCut, backend_ts>>;
+        boost::mp11::mp_product<CompCut, backend_ts>>;
   };
 
   template <SIMDWrapper::InstructionSet SIMDBackend, typename... Ts>
@@ -349,14 +349,14 @@ namespace ThOr {
    *    Combination123Cut -> cut applied to 3-particle tuples (only for N >=4)
    *    ...
    *    CombinationCut    -> cut applied to N-particle tuples (KS0 K+ pi-)
-   *    VertexCut         -> cut applied to fitted N-particle composite
+   *    CompositeCut      -> cut applied to fitted N-particle composite
    *                         particles (D0)
-   *  and N input locations: Input0 ... Input{N-1}.
+   *  and N input locations: Input1 ... Input{N}.
    *
    *  If (see above) an input container is used for multiple children, as in
    *    D+ -> pi+ pi+ K-
-   *  then the same container (of charged pions) must be given both as Input0
-   *  and Input1, with Input2 pointing to a container of charged kaons.
+   *  then the same container (of charged pions) must be given both as Input1
+   *  and Input2, with Input3 pointing to a container of charged kaons.
    *
    *  The algorithm supports two "modes of operation" for the input container
    *  types. It can be instantiated explicitly:
@@ -404,7 +404,7 @@ namespace ThOr {
 
     template <std::size_t... Is>
     static auto input_names( std::index_sequence<Is...> ) {
-      return std::array{KeyValue{"Input" + std::to_string( Is ), ""}...};
+      return std::array{KeyValue{"Input" + std::to_string( Is + 1 ), ""}...};
     }
 
     Combiner( const std::string& name, ISvcLocator* pSvcLocator )
@@ -446,13 +446,13 @@ namespace ThOr {
     template <std::size_t... ICombinationCuts>
     std::tuple<bool, Output_t> dispatch( EventContext const& evtCtx, InputTypes const&... raw_in,
                                          std::index_sequence<ICombinationCuts...> ) const {
-      // Try to get the vectorised versions of the combination and vertex cuts
-      auto const  vector_comb_cuts  = this->template getFunctors<CombCut<ICombinationCuts, SIMDBackend>...>();
-      auto const& vector_vertex_cut = this->template getFunctor<detail::Combiner::VertexCut<SIMDBackend>>();
+      // Try to get the vectorised versions of the combination and composite cuts
+      auto const  vector_comb_cuts     = this->template getFunctors<CombCut<ICombinationCuts, SIMDBackend>...>();
+      auto const& vector_composite_cut = this->template getFunctor<detail::Combiner::CompositeCut<SIMDBackend>>();
       // See if we got **all** the vectorised versions
-      if ( LIKELY( ( std::get<ICombinationCuts>( vector_comb_cuts ) && ... ) && vector_vertex_cut ) ) {
+      if ( LIKELY( ( std::get<ICombinationCuts>( vector_comb_cuts ) && ... ) && vector_composite_cut ) ) {
         // Yes, do all the hard work using the vectorised versions
-        return do_the_work<SIMDBackend>( evtCtx, raw_in..., vector_comb_cuts, vector_vertex_cut,
+        return do_the_work<SIMDBackend>( evtCtx, raw_in..., vector_comb_cuts, vector_composite_cut,
                                          std::index_sequence<ICombinationCuts...>{} );
       } else {
         // Presumably there was a cache miss, despatch to the implementation with JIT-able cuts
@@ -461,7 +461,7 @@ namespace ThOr {
         }
         return do_the_work<SIMDWrapper::Scalar>(
             evtCtx, raw_in..., this->template getFunctors<CombCut<ICombinationCuts, SIMDWrapper::Scalar>...>(),
-            this->template getFunctor<detail::Combiner::VertexCut<SIMDWrapper::Scalar>>(),
+            this->template getFunctor<detail::Combiner::CompositeCut<SIMDWrapper::Scalar>>(),
             std::index_sequence<ICombinationCuts...>{} );
       }
     }
@@ -473,8 +473,8 @@ namespace ThOr {
     std::tuple<bool, Output_t> do_the_work(
         EventContext const& evtCtx, InputTypes const&... raw_in,
         std::tuple<Functors::Functor<typename CombCut<ICombinationCuts, Backend>::Signature> const&...>
-                                                                                           unprepared_combination_cuts,
-        Functors::Functor<typename detail::Combiner::VertexCut<Backend>::Signature> const& unprepared_vertex_cut,
+                                                                                              unprepared_combination_cuts,
+        Functors::Functor<typename detail::Combiner::CompositeCut<Backend>::Signature> const& unprepared_composite_cut,
         std::index_sequence<ICombinationCuts...> ) const {
       using simd_t = typename SIMDWrapper::type_map_t<Backend>;
       // Get an iterable version of the input if need be -- this should be a noop if the input is already a zip
@@ -496,7 +496,7 @@ namespace ThOr {
         // Number of events processed so far
         auto const num_events = std::max( 10ul, m_output_container_reallocated.nFalseEntries() );
         // Number of accepted composite particles so far
-        auto const num_parts = std::max( 100ul, m_npassed_vertex_cut.nTrueEntries() );
+        auto const num_parts = std::max( 100ul, m_npassed_composite_cut.nTrueEntries() );
         // Average number of accepted composites in an event, and a fudge...
         constexpr auto fudge_factor = 4;
         storage.reserve( fudge_factor * num_parts / num_events );
@@ -526,13 +526,13 @@ namespace ThOr {
           storage, tmp_storage, in, algLocalMemResource,
           // Prepare functors, in case that's expensive
           std::tuple{std::get<ICombinationCuts>( unprepared_combination_cuts ).prepare( evtCtx )...},
-          unprepared_vertex_cut.prepare( evtCtx ),
+          unprepared_composite_cut.prepare( evtCtx ),
           // Buffer counters to reduce atomic operations
           std::tuple{m_buffer_shuffles[ICombinationCuts].buffer()...},
           std::tuple{m_combination_cut_loads[ICombinationCuts].buffer()...},
-          std::tuple{m_combination_cut_ideal_loads[ICombinationCuts].buffer()...}, m_vertex_cut_load.buffer(),
+          std::tuple{m_combination_cut_ideal_loads[ICombinationCuts].buffer()...}, m_composite_cut_load.buffer(),
           std::tuple{m_npassed_combination_cuts[ICombinationCuts].buffer()...}, m_npassed_vertex_fit.buffer(),
-          m_npassed_vertex_cut.buffer(),
+          m_npassed_composite_cut.buffer(),
           // Generate some integer packs for convenience
           std::index_sequence<ICombinationCuts...>{}, std::index_sequence_for<InputTypes...>{} );
 
@@ -604,13 +604,13 @@ namespace ThOr {
      *  of the combination generation.
      */
     template <SIMDWrapper::InstructionSet Backend, typename IterableInputTuple, typename CombinationCuts,
-              typename VertexCut, typename Counter, typename Counters, typename ShuffleCounters,
+              typename CompositeCut, typename Counter, typename Counters, typename ShuffleCounters,
               std::size_t... ICombinationCuts, std::size_t... IInput>
     void process_decays( Output_t& storage, Output_t& tmp_storage, IterableInputTuple const& in,
                          LHCb::Allocators::MemoryResource* algLocalMemResource, CombinationCuts combination_cuts,
-                         VertexCut vertex_cut, ShuffleCounters buffer_shuffles, Counters combination_cut_loads,
-                         Counters combination_cut_ideal_loads, Counter vertex_cut_load,
-                         Counters npassed_combination_cuts, Counter npassed_vertex_fit, Counter npassed_vertex_cut,
+                         CompositeCut composite_cut, ShuffleCounters buffer_shuffles, Counters combination_cut_loads,
+                         Counters combination_cut_ideal_loads, Counter composite_cut_load,
+                         Counters npassed_combination_cuts, Counter npassed_vertex_fit, Counter npassed_composite_cut,
                          std::index_sequence<ICombinationCuts...>, std::index_sequence<IInput...> ) const {
       using simd_t = typename SIMDWrapper::type_map_t<Backend>;
       using int_v  = typename simd_t::int_v;
@@ -651,7 +651,7 @@ namespace ThOr {
                                 int_v const& descriptor_indices ) {
         // Tasks:
         // - run a vertex fit
-        // - apply the vertex cut to surviving entries
+        // - apply the composite cut to surviving entries
         // - store surviving entries in the output container
 
         // Get the PIDs that we should assign to the new composites
@@ -688,20 +688,20 @@ namespace ThOr {
         // be quite well-populated. The exception might be the last call in each
         // event, where we are tidying up the leftovers.
         if ( LIKELY( fit_success_count ) ) {
-          // Last fit_count elements of tmp_storage are valid particles that we can apply the vertex cut to
+          // Last fit_count elements of tmp_storage are valid particles that we can apply the composite cut to
           auto const& particle_chunk = *tmp_storage_view.begin();
-          // Record how well we're SIMDifying the vertex cut
-          if constexpr ( Backend != SIMDWrapper::Scalar ) { vertex_cut_load += {fit_success_count, simd_t::size}; }
-          // Apply the vertex cut to those particles
-          auto const vertex_cut_mask =
-              fit_success_mask && vertex_cut( Functors::mask_arg, fit_success_mask, particle_chunk );
-          npassed_vertex_cut += {std::size_t( popcount( vertex_cut_mask ) ), fit_success_count};
+          // Record how well we're SIMDifying the composite cut
+          if constexpr ( Backend != SIMDWrapper::Scalar ) { composite_cut_load += {fit_success_count, simd_t::size}; }
+          // Apply the composite cut to those particles
+          auto const composite_cut_mask =
+              fit_success_mask && composite_cut( Functors::mask_arg, fit_success_mask, particle_chunk );
+          npassed_composite_cut += {std::size_t( popcount( composite_cut_mask ) ), fit_success_count};
           if ( UNLIKELY( this->msgLevel( MSG::DEBUG ) ) ) {
-            this->debug() << popcount( vertex_cut_mask ) << '/' << fit_success_count << " passed vertex cut ("
-                          << vertex_cut_mask << ')' << endmsg;
+            this->debug() << popcount( composite_cut_mask ) << '/' << fit_success_count << " passed composite cut ("
+                          << composite_cut_mask << ')' << endmsg;
           }
           // Copy the particles that matched into the actual storage
-          storage.template copy_back<simd_t>( tmp_storage, 0, vertex_cut_mask );
+          storage.template copy_back<simd_t>( tmp_storage, 0, composite_cut_mask );
         }
         // Clear the temporary storage again
         tmp_storage.clear();
@@ -713,7 +713,7 @@ namespace ThOr {
       // or 2 descriptors, with the second one being a top-level charge
       // conjugation of the first one.
       // The exception to this rule is that combination_indices are outside the
-      // loop, and we allow the vertexcut to execute with a mix of candidates
+      // loop, and we allow the composite cut to execute with a mix of candidates
       // from the two descriptors, which should allow the SIMD utilisation to be
       // increased slightly.
       for ( int n_decay = 0; n_decay < int( m_decays.size() ); ++n_decay ) {
@@ -778,7 +778,7 @@ namespace ThOr {
           ... );
       }
 
-      this->debug() << npassed_vertex_cut.nTrueEntries() << '/' << npassed_vertex_cut.nEntries()
+      this->debug() << npassed_composite_cut.nTrueEntries() << '/' << npassed_composite_cut.nEntries()
                     << " composites passed cut" << endmsg;
     }
 
@@ -932,7 +932,7 @@ namespace ThOr {
     /** Process a SIMD vector width worth of NBodies-body combinations defined
      *  by indices stored in 'combination_indices'. If NBodies < N this means
      *  creating (NBodies+1)-body combinations. If NBodies == N this means
-     *  running the vertex fit and applying the vertex cut.
+     *  running the vertex fit and applying the composite cut.
      *
      *  Only the first two template parameters need to be template parameters,
      *  everything else **could** be written out explicitly...
@@ -1097,12 +1097,12 @@ namespace ThOr {
         makeBufferShuffleCounters<shuffle_counter_t>( this, std::make_index_sequence<N - 1>{} )};
     /// Number of successfully fitted combinations
     mutable counter_t m_npassed_vertex_fit{this, "# passed vertex fit"};
-    /// Number of fitted vertices passing vertex cut
-    mutable counter_t m_npassed_vertex_cut{this, "# passed VertexCut"};
+    /// Number of fitted vertices passing composite cut
+    mutable counter_t m_npassed_composite_cut{this, "# passed CompositeCut"};
     /// How often we had to reallocate
     mutable counter_t m_output_container_reallocated{this, "# output reallocations"};
-    /// SIMD utilisation of the vertex cuts
-    mutable counter_t m_vertex_cut_load{this, "VertexCut SIMD utilisation"};
+    /// SIMD utilisation of the composite cuts
+    mutable counter_t m_composite_cut_load{this, "CompositeCut SIMD utilisation"};
     // decay descriptors
     ToolHandle<IDecodeSimpleDecayString> m_decoder{this, "DecayDecoder", "DecodeSimpleDecayString"};
     // decays
diff --git a/Phys/ParticleCombiners/src/TrackCombiner.cpp b/Phys/ParticleCombiners/src/TrackCombiner.cpp
deleted file mode 100644
index 1e14d793f6e76a6b56236967eb60292283a3a8d2..0000000000000000000000000000000000000000
--- a/Phys/ParticleCombiners/src/TrackCombiner.cpp
+++ /dev/null
@@ -1,137 +0,0 @@
-/*****************************************************************************\
-* (c) Copyright 2000-2019 CERN for the benefit of the LHCb Collaboration      *
-*                                                                             *
-* This software is distributed under the terms of the GNU General Public      *
-* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".   *
-*                                                                             *
-* In applying this licence, CERN does not waive the privileges and immunities *
-* granted to it by virtue of its status as an Intergovernmental Organization  *
-* or submit itself to any jurisdiction.                                       *
-\*****************************************************************************/
-
-#include <assert.h>
-
-#include <functional>
-
-#include <boost/algorithm/string/join.hpp>
-
-#include "LoKi/Combiner.h"
-#include "LoKi/VertexCuts.h"
-#include "TrackKernel/TrackStateVertex.h"
-#include "TrackKernel/TrackVertexUtils.h"
-
-#include "TrackCombiner.h"
-
-/// Vertex fitters used by the class
-namespace {
-  /// Fit an N-body vertex using TrackStateVertex
-  void fit_vertex_n_body( LHCb::RecVertex& vertex ) {
-    auto const&                     tracks = vertex.tracks();
-    std::vector<const LHCb::State*> states;
-    for ( auto const& track : tracks ) {
-      auto state = track->stateAt( LHCb::State::ClosestToBeam );
-      assert( state != nullptr );
-      states.push_back( state );
-    }
-
-    LHCb::TrackStateVertex fitter( states );
-
-    // Fill the vertex with the fit result
-    vertex.setCovMatrix( fitter.covMatrix() );
-    vertex.setPosition( fitter.position() );
-    vertex.setChi2AndDoF( fitter.chi2(), fitter.nDoF() );
-  }
-} // namespace
-
-DECLARE_COMPONENT( TrackCombiner )
-
-TrackCombiner::TrackCombiner( const std::string& name, ISvcLocator* pSvcLocator )
-    : Transformer( name, pSvcLocator, {KeyValue{"InputTracks", LHCb::TrackLocation::Default}},
-                   {KeyValue{"OutputVertices", ""}} )
-    , m_combination_cut( LoKi::Cuts::VFALSE )
-    , m_vertex_cut( LoKi::Cuts::VFALSE ) {}
-
-StatusCode TrackCombiner::initialize() {
-  auto sc = Transformer::initialize();
-  if ( sc.isFailure() ) { return sc; }
-
-  // Can't fit a one-body vertex
-  if ( m_nbodies < 2 ) { return Error( "Cannot fit a N-body vertex with N < 2" ); }
-
-  // Form the preamble by concatenating the list items with newlines
-  m_preamble = boost::algorithm::join( m_preamble_def.value(), "\n" );
-
-  // Initialise the pre-fit selection object
-  sc = m_factory->get( m_combination_cut_def, m_combination_cut, m_preamble );
-  if ( sc.isFailure() ) { return sc; }
-
-  // Initialise the post-fit selection object
-  sc = m_factory->get( m_vertex_cut_def, m_vertex_cut, m_preamble );
-  if ( sc.isFailure() ) { return sc; }
-
-  return StatusCode::SUCCESS;
-}
-
-std::vector<LHCb::RecVertex> TrackCombiner::operator()( Pr::Selection<LHCb::Event::v1::Track> const& tracks ) const {
-  using Tracks        = std::vector<const LHCb::Event::v1::Track*>;
-  using Combiner      = LoKi::Combiner_<Tracks>;
-  using CombinerRange = typename Combiner::Range;
-
-  // Create vector of pointers from vector of objects,
-  // as needed for the combinatorics engine
-  // TODO This is a hack; the combinatorics engine should accept Pr::Selection
-  Tracks ptracks;
-  ptracks.reserve( tracks.size() );
-  for ( auto const& track : tracks ) { ptracks.push_back( &track ); }
-
-  // Create the combinatorics iterator
-  Combiner      combiner;
-  CombinerRange range( ptracks.begin(), ptracks.end() );
-  for ( int i = 0; i < m_nbodies; ++i ) { combiner.add( range ); }
-
-  // Create vertices from each combination
-  std::vector<LHCb::RecVertex> vertices;
-  auto                         ncombinations = 0;
-  // A comparator to impose pair-wise order within each combination
-  std::less<> comparator;
-  for ( ; combiner.valid(); combiner.next() ) {
-    // Filter out combinations that are not ordered
-    // Because only one permutation of a given set of elements will satisfy the
-    // ordering criteria, we will skip all other permutations (which is good,
-    // because they would result in the same vertex)
-    // This criteria also means we skip combinations with duplicate elements
-    if ( !combiner.unique( comparator ) ) { continue; }
-
-    // Get the current combination
-    Tracks combination( combiner.dim() );
-    combiner.current( combination.begin() );
-    ++ncombinations;
-
-    // Create an (effectively uninitialised) vertex filled with the tracks
-    LHCb::RecVertex vertex;
-    for ( auto const* track : combination ) { vertex.addToTracks( track ); }
-
-    // Apply the pre-fit selection
-    auto passed_combination = m_combination_cut( &vertex );
-    m_npassed_combination_cut += passed_combination;
-    if ( !passed_combination ) { continue; }
-
-    // Perform the fit
-    fit_vertex_n_body( vertex );
-    auto fit_chi2   = vertex.chi2();
-    auto passed_fit = fit_chi2 > 0;
-    m_npassed_vertex_fit += passed_fit;
-    if ( !passed_fit ) { continue; }
-
-    // Apply the post-fit selection
-    auto passed_vertex = m_vertex_cut( &vertex );
-    m_npassed_vertex_cut += passed_vertex;
-    if ( !passed_vertex ) { continue; }
-
-    vertices.push_back( std::move( vertex ) );
-  }
-
-  m_ncombinations += ncombinations;
-
-  return vertices;
-}
diff --git a/Phys/ParticleCombiners/src/TrackCombiner.h b/Phys/ParticleCombiners/src/TrackCombiner.h
deleted file mode 100644
index e77fb2d60dac0c69aa8048c3f18586040426d832..0000000000000000000000000000000000000000
--- a/Phys/ParticleCombiners/src/TrackCombiner.h
+++ /dev/null
@@ -1,95 +0,0 @@
-/*****************************************************************************\
-* (c) Copyright 2000-2019 CERN for the benefit of the LHCb Collaboration      *
-*                                                                             *
-* This software is distributed under the terms of the GNU General Public      *
-* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".   *
-*                                                                             *
-* In applying this licence, CERN does not waive the privileges and immunities *
-* granted to it by virtue of its status as an Intergovernmental Organization  *
-* or submit itself to any jurisdiction.                                       *
-\*****************************************************************************/
-
-#pragma once
-
-#include "GaudiAlg/Transformer.h"
-
-#include "Event/RecVertex.h"
-#include "Event/Track_v1.h"
-#include "LoKi/IHybridFactory.h"
-#include "PrKernel/PrSelection.h"
-
-/** @class TrackCombiner TrackCombiner.h
- *
- * @brief N-body combiner for LHCb::Event::v1::Track objects.
- *
- * All input Track objects will be taken from a single input container. For
- * each combination of Track objects, a vertex is fitted from the set of
- * LHCb::State::ClosestToBeam states.
- *
- * @todo Consider adding an overlap checking tool, e.g. LoKi::CheckOverlap, for
- * usage before the combination cut
- * @todo Reserve the output vertex vector once we have an understanding of the
- * number of vertices a typical selection creates
- * @todo The comparator used for ensuring unique track object permutations does
- * not guarantee ordering from one application run to the next (as it uses
- * pointer values), neither for the Track vector inserted into the vertex nor
- * the vertex creation order itself. Should consider ordering these for
- * reproducibility (ideally after as much selection as possible).
- * @todo The RecVertex creation/destruction, one per combination, adds up to a
- * huge overall cost (30-40%). Should consider avoiding this, e.g. optimising
- * RecVertex or using a (as-yet non-existent) family of ATr functors that can
- * act on Track vectors directly.
- *
- * Example usage:
- *
- * @code {.py}
- * combiner = TrackCombiner("TwoBodyCombiner")
- * combiner.NBodies = 2
- * combiner.Preamble = [
- *     "from GaudiKernel.SystemOfUnits import MeV, mm",
- *     "from LoKiPhys.decorators import RV_MASS, VCHI2PDOF, VZ"
- * ]
- * combiner.CombinationCut = "RV_MASS('mu-', 'mu+') < 1000*MeV"
- * combiner.VertexCut = "(VCHI2PDOF < 8) & (VZ > -100*mm)"
- * @endcode
- *
- * @see CombineParticles Similar class for LHCb::Particle objects
- *
- */
-class TrackCombiner : public Gaudi::Functional::Transformer<std::vector<LHCb::RecVertex>(
-                          Pr::Selection<LHCb::Event::v1::Track> const& )> {
-
-public:
-  TrackCombiner( const std::string& name, ISvcLocator* pSvcLocator );
-
-  StatusCode initialize() override;
-
-  std::vector<LHCb::RecVertex> operator()( Pr::Selection<LHCb::Event::v1::Track> const& tracks ) const override;
-
-private:
-  ToolHandle<LoKi::IHybridFactory> m_factory = {"LoKi::Hybrid::Tool/HybridFactory:PUBLIC", this};
-
-  Gaudi::Property<int> m_nbodies{this, "NBodies", 2, "Number of tracks to include in each combination"};
-  Gaudi::Property<std::vector<std::string>> m_preamble_def{this, "Preamble", {}, "List of preamble strings"};
-  Gaudi::Property<std::string>              m_combination_cut_def{this, "CombinationCut", "DEFINE_ME",
-                                                     "Cut string applied to the N-body combination"};
-  Gaudi::Property<std::string>              m_vertex_cut_def{this, "VertexCut", "DEFINE_ME",
-                                                "Cut string applied to the fitted vertex"};
-
-  /// Newline-delimited preamble
-  std::string m_preamble;
-
-  /// Cut object for the pre-fit combination selection
-  LoKi::Types::VCut m_combination_cut;
-  /// Cut object for the fitted vertex selection
-  LoKi::Types::VCut m_vertex_cut;
-
-  /// Number of N-body Track combinations created
-  mutable Gaudi::Accumulators::StatCounter<unsigned int> m_ncombinations{this, "Nb created combinations"};
-  /// Number of combinations passing the combination cut
-  mutable Gaudi::Accumulators::BinomialCounter<> m_npassed_combination_cut{this, "# passed combination cut"};
-  /// Number of successfully fitted combinations
-  mutable Gaudi::Accumulators::BinomialCounter<> m_npassed_vertex_fit{this, "# passed vertex fit"};
-  /// Number of fitted vertices passing vertex cut
-  mutable Gaudi::Accumulators::BinomialCounter<> m_npassed_vertex_cut{this, "# passed vertex cut"};
-};
diff --git a/Phys/ParticleCombiners/tests/options/test_trackcombiner.py b/Phys/ParticleCombiners/tests/options/test_trackcombiner.py
deleted file mode 100644
index 647695ea1a59f1ebd6b4fecee0ef17438aca7b2e..0000000000000000000000000000000000000000
--- a/Phys/ParticleCombiners/tests/options/test_trackcombiner.py
+++ /dev/null
@@ -1,113 +0,0 @@
-from __future__ import print_function
-###############################################################################
-# (c) Copyright 2000-2019 CERN for the benefit of the LHCb Collaboration      #
-#                                                                             #
-# This software is distributed under the terms of the GNU General Public      #
-# Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".   #
-#                                                                             #
-# In applying this licence, CERN does not waive the privileges and immunities #
-# granted to it by virtue of its status as an Intergovernmental Organization  #
-# or submit itself to any jurisdiction.                                       #
-###############################################################################
-from builtins import range
-from Configurables import (
-    ApplicationMgr, LHCb__Converters__Track__v1__fromLHCbTrackVector, LHCbApp,
-    MakeSelection__Track_v1, TrackCombiner, UnpackTrack)
-import GaudiPython as GP
-from PRConfig import TestFileDB
-
-
-def remove_davinci_hybrid_modules():
-    """Remove Hybrid modules that require a DVAlgorithm.
-
-    Some modules import functors that depend on the DVAlgorithm context being
-    available. The LoKi::Hybrid::Tool tool loads these modules by default, but
-    we don't want nor need them, so remove them from the default list.
-    """
-    # List of modules we will delete from the default list
-    dv_modules = ["LoKiPhys.decorators", "LoKiArrayFunctors.decorators"]
-
-    from Configurables import LoKi__Hybrid__Tool
-    # Fetch the public tool instance used by the TrackCombiner algorithm
-    factory = LoKi__Hybrid__Tool("HybridFactory")
-    for m in dv_modules:
-        del factory.Modules[factory.Modules.index(m)]
-
-
-# Unpack pRec/Track/Best
-unpacker = UnpackTrack()
-
-# Convert the KeyedContainer to a std::vector
-converter = LHCb__Converters__Track__v1__fromLHCbTrackVector(
-    "ConvertTrackVector")
-converter.InputTracksName = unpacker.getDefaultProperty("OutputName")
-converter.OutputTracksName = str(converter.InputTracksName) + "Vector"
-
-# Apply a dummy selection (selects all tracks)
-selection = MakeSelection__Track_v1("SelectionMaker")
-selection.Input = converter.OutputTracksName
-selection.Output = str(selection.Input) + "Selected"
-
-combiner = TrackCombiner("TwoBodyCombiner")
-combiner.NBodies = 2
-combiner.Preamble = [
-    "from GaudiKernel.SystemOfUnits import MeV, GeV",
-    "from LoKiPhys.decorators import RV_MASS, RV_PT, RV_TrNUM, VCHI2PDOF",
-    "from LoKiTrack.decorators import TrLONG"
-]
-# Form high-pT phi(1020) candidates
-combiner.CombinationCut = "(RV_TrNUM(TrLONG) == 2) & (RV_PT > 3*GeV) & (RV_MASS('K-', 'K+') < 1040*MeV)"
-combiner.VertexCut = "(VCHI2PDOF < 8)"
-combiner.InputTracks = selection.Output
-combiner.OutputVertices = 'Rec/Vertex/Phi'
-remove_davinci_hybrid_modules()
-
-# Define the sequence
-ApplicationMgr().TopAlg = [unpacker, converter, selection, combiner]
-
-app = LHCbApp()
-app.EvtMax = 100
-
-# Pick a file that has the reconstruction available
-f = TestFileDB.test_file_db[
-    "2017HLTCommissioning_BnoC_MC2015_Bs2PhiPhi_L0Processed"]
-f.setqualifiers(configurable=app)
-f.run(configurable=app, withDB=True)
-
-appMgr = GP.AppMgr()
-TES = appMgr.evtsvc()
-
-# FIXME
-# without this calling `vertex.chi2PerDoF()` causes a warning on stderr
-# input_line_495:10:7: warning: ignoring return value of function declared with 'nodiscard' attribute [-Wunused-result]
-#       ((const LHCb::VertexBase*)obj)->chi2PerDoF();
-# This is due to a bug in ROOT, now reported. Whenever it's fixed, this should be removed
-import cppyy
-cppyy.gbl.gInterpreter.Declare(
-    '''#pragma clang diagnostic ignored "-Wunused-result"''')
-
-event = 0
-
-while event < app.EvtMax:
-    anydata = TES[str(combiner.OutputVertices)]
-    if not anydata:
-        event += 1
-        appMgr.run(1)
-        continue
-
-    container = anydata.getData()
-    for i in range(container.size()):
-        vertex = container.at(i)
-        tracks = vertex.tracks()
-        assert len(tracks) == 2, (i, vertex, len(tracks))
-
-        chi2pdof = vertex.chi2PerDoF()
-        pt = sum(t.pT().value() for t in tracks)
-        all_long = all(t.type() == t.Long for t in tracks)
-
-        # Check that our cuts have been applied
-        assert chi2pdof < 8, (i, vertex, chi2pdof)
-        assert pt > 3000, (i, vertex, pt)
-
-    event += 1
-    appMgr.run(1)
diff --git a/Phys/ParticleCombiners/tests/qmtest/particlecombiners.qms/trackcombiner.qmt b/Phys/ParticleCombiners/tests/qmtest/particlecombiners.qms/trackcombiner.qmt
deleted file mode 100644
index e643ee6e429495cec38bbcacdf11014207d18d27..0000000000000000000000000000000000000000
--- a/Phys/ParticleCombiners/tests/qmtest/particlecombiners.qms/trackcombiner.qmt
+++ /dev/null
@@ -1,26 +0,0 @@
-<?xml version="1.0" ?><!DOCTYPE extension  PUBLIC '-//QM/2.3/Extension//EN'  'http://www.codesourcery.com/qm/dtds/2.3/-//qm/2.3/extension//en.dtd'>
-<!--
-    (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration
-
-    This software is distributed under the terms of the GNU General Public
-    Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".
-
-    In applying this licence, CERN does not waive the privileges and immunities
-    granted to it by virtue of its status as an Intergovernmental Organization
-    or submit itself to any jurisdiction.
--->
-<extension class="GaudiTest.GaudiExeTest" kind="test">
-  <argument name="program"><text>python</text></argument>
-  <argument name="args"><set>
-      <text>../options/test_trackcombiner.py</text>
-  </set></argument>
-  <argument name="reference"><text>../refs/trackcombiner.ref</text></argument>
-  <argument name="validator">
-    <text>
-from GaudiConf.QMTest.LHCbExclusions import preprocessor
-pp = preprocessor + LineSkipper(["|-ThreadPoolSize ="])
-validateWithReference(preproc = pp)
-# TODO ideally, check only counters of `TwoBodyCombiner`
-    </text>
-  </argument>
-</extension>