From 8bd7f12a0431b2ac5c4b73c676d3e72f3a39c5bb Mon Sep 17 00:00:00 2001
From: Miguel Ramos Pernas <miguel.ramos.pernas@cern.ch>
Date: Fri, 15 Oct 2021 11:57:18 +0200
Subject: [PATCH 01/26] Utilities for the ThOr combiner

---
 Phys/SelKernel/include/SelKernel/Utilities.h | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/Phys/SelKernel/include/SelKernel/Utilities.h b/Phys/SelKernel/include/SelKernel/Utilities.h
index ff73ce1a562..4d72084e099 100644
--- a/Phys/SelKernel/include/SelKernel/Utilities.h
+++ b/Phys/SelKernel/include/SelKernel/Utilities.h
@@ -194,6 +194,13 @@ namespace Sel::Utils {
   constexpr bool                    has_tracklike_API =
       details::is_proxy<T>::value&& details::unpack_and_contains<LHCb::Event::v3::Tracks, T>::value;
 
+  template <typename T>
+  constexpr bool                    has_neutrals_API =
+    details::is_proxy<T>::value&& details::unpack_and_contains<LHCb::Event::HypothesesWithOrigin, T>::value;
+
+  template <typename T>
+  constexpr bool has_basic_particle_API = (has_tracklike_API<T> || has_neutrals_API<T>);
+
   using StateLocation = LHCb::Event::v3::Tracks::StateLocation;
   template <typename Tracklike>
   auto get_kinematic_state( Tracklike track ) {
-- 
GitLab


From fde79c3533d7ab54c244d1f0cae4a3099c484ad0 Mon Sep 17 00:00:00 2001
From: Miguel Ramos Pernas <miguel.ramos.pernas@cern.ch>
Date: Fri, 15 Oct 2021 11:58:53 +0200
Subject: [PATCH 02/26] lb-format

---
 Phys/SelKernel/include/SelKernel/Utilities.h | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Phys/SelKernel/include/SelKernel/Utilities.h b/Phys/SelKernel/include/SelKernel/Utilities.h
index 4d72084e099..6a6040c72af 100644
--- a/Phys/SelKernel/include/SelKernel/Utilities.h
+++ b/Phys/SelKernel/include/SelKernel/Utilities.h
@@ -196,10 +196,10 @@ namespace Sel::Utils {
 
   template <typename T>
   constexpr bool                    has_neutrals_API =
-    details::is_proxy<T>::value&& details::unpack_and_contains<LHCb::Event::HypothesesWithOrigin, T>::value;
+      details::is_proxy<T>::value&& details::unpack_and_contains<LHCb::Event::HypothesesWithOrigin, T>::value;
 
   template <typename T>
-  constexpr bool has_basic_particle_API = (has_tracklike_API<T> || has_neutrals_API<T>);
+  constexpr bool has_basic_particle_API = ( has_tracklike_API<T> || has_neutrals_API<T> );
 
   using StateLocation = LHCb::Event::v3::Tracks::StateLocation;
   template <typename Tracklike>
-- 
GitLab


From a735b9da9ba9f61745ff1a59e1395b1085c747f1 Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Fri, 22 Oct 2021 10:11:44 +0100
Subject: [PATCH 03/26] renaming

---
 Phys/SelKernel/include/SelKernel/Utilities.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Phys/SelKernel/include/SelKernel/Utilities.h b/Phys/SelKernel/include/SelKernel/Utilities.h
index a105e7dbe37..904378150f6 100644
--- a/Phys/SelKernel/include/SelKernel/Utilities.h
+++ b/Phys/SelKernel/include/SelKernel/Utilities.h
@@ -197,7 +197,7 @@ namespace Sel::Utils {
 
   template <typename T>
   constexpr bool                    has_neutrals_API =
-      details::is_proxy<T>::value&& details::unpack_and_contains<LHCb::Event::HypothesesWithOrigin, T>::value;
+      details::is_proxy<T>::value&& details::unpack_and_contains<LHCb::Event::CaloHypothesesWithOrigin, T>::value;
 
   template <typename T>
   constexpr bool has_basic_particle_API = ( has_tracklike_API<T> || has_neutrals_API<T> );
-- 
GitLab


From 448ddfc8aa09139fd7de048c7ba2fd9017de53ed Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Fri, 5 Nov 2021 10:48:01 +0000
Subject: [PATCH 04/26] Add empty function to return a state vector from
 neutral objects

---
 Phys/SelTools/include/SelTools/Utilities.h | 17 ++++++++++++++++-
 1 file changed, 16 insertions(+), 1 deletion(-)

diff --git a/Phys/SelTools/include/SelTools/Utilities.h b/Phys/SelTools/include/SelTools/Utilities.h
index 294caba6621..53fa5b80be1 100644
--- a/Phys/SelTools/include/SelTools/Utilities.h
+++ b/Phys/SelTools/include/SelTools/Utilities.h
@@ -82,4 +82,19 @@ namespace Sel {
     }
     return s;
   }
-} // namespace Sel
\ No newline at end of file
+
+  /** @fn    stateVectorFromOther
+   @brief Convert a particle (4-momentum + position) to 4-state (at some z).
+
+   TODO: This must be properly filled
+   */
+  template <typename particle_t>
+  auto stateVectorFromOther( particle_t const& p ) {
+
+    using float_v = std::decay_t<decltype( p.pt() )>;
+
+    State4<float_v> s{};
+
+    return s;
+  }
+} // namespace Sel
-- 
GitLab


From c87bbbd6af3c1107a5010d7b51d09cd970b6819f Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Mon, 22 Nov 2021 12:37:51 +0000
Subject: [PATCH 05/26] Rename state accessor for neutrals

---
 Phys/SelTools/include/SelTools/Utilities.h | 15 +++++++++++++--
 1 file changed, 13 insertions(+), 2 deletions(-)

diff --git a/Phys/SelTools/include/SelTools/Utilities.h b/Phys/SelTools/include/SelTools/Utilities.h
index 53fa5b80be1..5f2be8da784 100644
--- a/Phys/SelTools/include/SelTools/Utilities.h
+++ b/Phys/SelTools/include/SelTools/Utilities.h
@@ -83,17 +83,28 @@ namespace Sel {
     return s;
   }
 
-  /** @fn    stateVectorFromOther
+  /** @fn    stateVectorFromNeutral
    @brief Convert a particle (4-momentum + position) to 4-state (at some z).
 
    TODO: This must be properly filled
    */
   template <typename particle_t>
-  auto stateVectorFromOther( particle_t const& p ) {
+  auto stateVectorFromNeutral( particle_t const& p ) {
 
     using float_v = std::decay_t<decltype( p.pt() )>;
 
+    auto const invpz = 1.f / p.pz();
+
     State4<float_v> s{};
+    s.x()  = p.hypothesis_position_x();
+    s.y()  = p.hypothesis_position_y();
+    s.z()  = p.hypothesis_position_z();
+    s.tx() = p.px() * invpz;
+    s.ty() = p.py() * invpz;
+
+    s.covXX() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::Mat<float_v, 2, 3>>();
+    s.covTT() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::Mat<float_v, 2, 3>>();
+    s.covXT() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::Mat<float_v, 2, 3>>();
 
     return s;
   }
-- 
GitLab


From 995f7ce19ad67d4f1a62ac4bbbbf38b0143593ad Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Wed, 24 Nov 2021 12:06:51 +0000
Subject: [PATCH 06/26] Fix setting the state as empty

---
 Phys/SelTools/include/SelTools/Utilities.h | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/Phys/SelTools/include/SelTools/Utilities.h b/Phys/SelTools/include/SelTools/Utilities.h
index 5f2be8da784..e06e304c388 100644
--- a/Phys/SelTools/include/SelTools/Utilities.h
+++ b/Phys/SelTools/include/SelTools/Utilities.h
@@ -102,9 +102,9 @@ namespace Sel {
     s.tx() = p.px() * invpz;
     s.ty() = p.py() * invpz;
 
-    s.covXX() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::Mat<float_v, 2, 3>>();
-    s.covTT() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::Mat<float_v, 2, 3>>();
-    s.covXT() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::Mat<float_v, 2, 3>>();
+    s.covXX() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::MatSym<float_v, 2>>();
+    s.covTT() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::MatSym<float_v, 2>>();
+    s.covXT() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::Mat<float_v, 2, 2>>();
 
     return s;
   }
-- 
GitLab


From c338b482a11633d007d1dde8309fabf17839338f Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Fri, 26 Nov 2021 11:44:34 +0000
Subject: [PATCH 07/26] Change the name of the container

---
 Phys/SelKernel/include/SelKernel/Utilities.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Phys/SelKernel/include/SelKernel/Utilities.h b/Phys/SelKernel/include/SelKernel/Utilities.h
index 3c7a6c48ed3..60ba84228ae 100644
--- a/Phys/SelKernel/include/SelKernel/Utilities.h
+++ b/Phys/SelKernel/include/SelKernel/Utilities.h
@@ -201,7 +201,7 @@ namespace Sel::Utils {
 
   template <typename T>
   constexpr bool                    has_neutrals_API =
-      details::is_proxy<T>::value&& details::unpack_and_contains<LHCb::Event::CaloHypothesesWithOrigin, T>::value;
+      details::is_proxy<T>::value&& details::unpack_and_contains<LHCb::Event::CaloHypothesesWithDirection, T>::value;
 
   template <typename T>
   constexpr bool has_basic_particle_API = ( has_tracklike_API<T> || has_neutrals_API<T> );
-- 
GitLab


From e8603da46a3ddb1cb27fb74159f9954b4f3c22f1 Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Thu, 2 Dec 2021 19:35:05 +0000
Subject: [PATCH 08/26] Add some interfaces for neutrals

---
 Phys/SelAlgorithms/src/MakeZip.cpp                  | 3 +++
 Phys/SelAlgorithms/src/Monitor.cpp                  | 2 ++
 Phys/SelTools/include/SelTools/DistanceCalculator.h | 2 ++
 3 files changed, 7 insertions(+)

diff --git a/Phys/SelAlgorithms/src/MakeZip.cpp b/Phys/SelAlgorithms/src/MakeZip.cpp
index 0daf3654532..e5dca6a180b 100644
--- a/Phys/SelAlgorithms/src/MakeZip.cpp
+++ b/Phys/SelAlgorithms/src/MakeZip.cpp
@@ -63,6 +63,9 @@ struct MakeZip final : public Helper<OutputType, Inputs...>::Base {
 using MakeZip__ChargedBasics_as_Particles = MakeZip<LHCb::Event::Particles, LHCb::Event::ChargedBasics>;
 DECLARE_COMPONENT_WITH_ID( MakeZip__ChargedBasics_as_Particles, "MakeZip__ChargedBasics_as_Particles" )
 
+using MakeZip__NeutralBasics_as_Particles = MakeZip<LHCb::Event::Particles, LHCb::Event::NeutralBasics>;
+DECLARE_COMPONENT_WITH_ID( MakeZip__NeutralBasics_as_Particles, "MakeZip__NeutralBasics_as_Particles" )
+
 using MakeZip__PrFittedForwardTracks__BestVertexRelations =
     MakeZip<deduce_return_type_t, LHCb::Pr::Fitted::Forward::Tracks, BestVertexRelations>;
 DECLARE_COMPONENT_WITH_ID( MakeZip__PrFittedForwardTracks__BestVertexRelations,
diff --git a/Phys/SelAlgorithms/src/Monitor.cpp b/Phys/SelAlgorithms/src/Monitor.cpp
index 8154407d841..59230867ba7 100644
--- a/Phys/SelAlgorithms/src/Monitor.cpp
+++ b/Phys/SelAlgorithms/src/Monitor.cpp
@@ -31,6 +31,8 @@ namespace SelAlgorithms {
 
   DECLARE_COMPONENT_WITH_ID( Monitor<LHCb::Event::ChargedBasics>, "Monitor__ChargedBasics" )
 
+  DECLARE_COMPONENT_WITH_ID( Monitor<LHCb::Event::NeutralBasics>, "Monitor__NeutralBasics" )
+
   DECLARE_COMPONENT_WITH_ID( Monitor<LHCb::Event::Composites>, "Monitor__Composites" )
 
   DECLARE_COMPONENT_WITH_ID( Monitor<LHCb::Particle::Range>, "Monitor__ParticleRange" )
diff --git a/Phys/SelTools/include/SelTools/DistanceCalculator.h b/Phys/SelTools/include/SelTools/DistanceCalculator.h
index 4328e55e9f3..4b0ef410d79 100644
--- a/Phys/SelTools/include/SelTools/DistanceCalculator.h
+++ b/Phys/SelTools/include/SelTools/DistanceCalculator.h
@@ -45,6 +45,8 @@ namespace Sel {
         return p.state( StateLocation::ClosestToBeam );
       } else if constexpr ( Sel::type_traits::has_closestToBeamState_v<Particle> ) {
         return p.closestToBeamState();
+      } else if constexpr ( Sel::Utils::has_neutrals_API<Particle> ) {
+        return stateVectorFromNeutral( p );
       } else {
         // Make a state vector from a particle (position + 4-momentum)
         return stateVectorFromComposite( p );
-- 
GitLab


From b0103c4af653c1632b02e0e77968cb8a1dd1ecd0 Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Fri, 3 Dec 2021 10:47:11 +0000
Subject: [PATCH 09/26] Remove the MakeZip__NeutralBasics_as_Particles
 algorithm

---
 Phys/SelAlgorithms/src/MakeZip.cpp | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/Phys/SelAlgorithms/src/MakeZip.cpp b/Phys/SelAlgorithms/src/MakeZip.cpp
index e5dca6a180b..0daf3654532 100644
--- a/Phys/SelAlgorithms/src/MakeZip.cpp
+++ b/Phys/SelAlgorithms/src/MakeZip.cpp
@@ -63,9 +63,6 @@ struct MakeZip final : public Helper<OutputType, Inputs...>::Base {
 using MakeZip__ChargedBasics_as_Particles = MakeZip<LHCb::Event::Particles, LHCb::Event::ChargedBasics>;
 DECLARE_COMPONENT_WITH_ID( MakeZip__ChargedBasics_as_Particles, "MakeZip__ChargedBasics_as_Particles" )
 
-using MakeZip__NeutralBasics_as_Particles = MakeZip<LHCb::Event::Particles, LHCb::Event::NeutralBasics>;
-DECLARE_COMPONENT_WITH_ID( MakeZip__NeutralBasics_as_Particles, "MakeZip__NeutralBasics_as_Particles" )
-
 using MakeZip__PrFittedForwardTracks__BestVertexRelations =
     MakeZip<deduce_return_type_t, LHCb::Pr::Fitted::Forward::Tracks, BestVertexRelations>;
 DECLARE_COMPONENT_WITH_ID( MakeZip__PrFittedForwardTracks__BestVertexRelations,
-- 
GitLab


From b15ab2cceefa2bdfab4d8d128f9c7cc5eaa6d0df Mon Sep 17 00:00:00 2001
From: Miguel Ramos Pernas <miguel.ramos.pernas@cern.ch>
Date: Fri, 15 Oct 2021 11:56:42 +0200
Subject: [PATCH 10/26] Extend ThOr combiner to use neutrals

---
 Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h b/Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h
index 88cc1324eac..628568c515d 100644
--- a/Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h
+++ b/Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h
@@ -165,11 +165,11 @@ namespace ThOr::detail::Combiner {
 
   template <typename T0, typename T1>
   auto have_overlap( T0 const& p0, T1 const& p1 ) {
-    using Sel::Utils::has_tracklike_API;
-    if constexpr ( Sel::Utils::has_tracklike_API<T0> && Sel::Utils::has_tracklike_API<T1> ) {
+    using Sel::Utils::has_basic_particle_API;
+    if constexpr ( Sel::Utils::has_basic_particle_API<T0> && Sel::Utils::has_basic_particle_API<T1> ) {
       // track-like/track-like
       return p0.unique_id() == p1.unique_id();
-    } else if constexpr ( Sel::Utils::has_tracklike_API<T0> && !Sel::Utils::has_tracklike_API<T1> ) {
+    } else if constexpr ( Sel::Utils::has_basic_particle_API<T0> && !Sel::Utils::has_basic_particle_API<T1> ) {
       // track-like/composite
       using mask_t = decltype( p0.unique_id() == p0.unique_id() );
       mask_t mask{false};
@@ -177,7 +177,7 @@ namespace ThOr::detail::Combiner {
         mask = mask | ( p1.descendant_unique_id( i ) == p0.unique_id() );
       }
       return mask;
-    } else if constexpr ( !Sel::Utils::has_tracklike_API<T0> && Sel::Utils::has_tracklike_API<T1> )
+    } else if constexpr ( !Sel::Utils::has_basic_particle_API<T0> && Sel::Utils::has_basic_particle_API<T1> )
       // composite/track-like
       return have_overlap( p1, p0 );
     else {
-- 
GitLab


From 8cd40d564cb311afab71c1c39b987e6f0a91d2a6 Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Fri, 5 Nov 2021 10:49:01 +0000
Subject: [PATCH 11/26] Add support for CALO objects in the vertex fitter

---
 .../include/VertexFit/ParticleVertexFitter.h  | 27 ++++++++++++-------
 1 file changed, 18 insertions(+), 9 deletions(-)

diff --git a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
index a2178deded6..1e2da8abae4 100644
--- a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
+++ b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
@@ -118,12 +118,17 @@ namespace Sel::Fitters {
     enum struct ChildType : unsigned char {
       TrackWithVelo = 0,
       TrackWithoutVelo,
+      Track,
       Composite,
       Resonance,
       Other,
       NumTypes
     };
 
+    template <ChildType C>
+    static constexpr bool has_full_state_v = ( C == ChildType::Composite || C == ChildType::Resonance ||
+                                               C == ChildType::Other );
+
     /** Indirectly stolen from TrackVertexUtils::poca
      */
     template <typename float_v, typename mask_v>
@@ -180,12 +185,12 @@ namespace Sel::Fitters {
     }
 
     // Class for daughters with a straight-line trajectory
-    template <typename particle_t, bool represents_composite>
+    template <typename particle_t, ChildType child_type>
     struct VertexTraj final {
       using simd_t  = typename particle_t::simd_t;
       using int_v   = typename simd_t::int_v;
       using float_v = typename simd_t::float_v;
-      using state_t = std::conditional_t<represents_composite, State4<float_v>, MiniState<float_v>>;
+      using state_t = std::conditional_t<has_full_state_v<child_type>, State4<float_v>, MiniState<float_v>>;
 
     private:
       particle_t            m_particle; // TODO maybe just cache some information from this?
@@ -252,10 +257,11 @@ namespace Sel::Fitters {
         // ParticleVertexFitter this was done with a full specialisation, but
         // now it would have to be a partial specialisation and it all gets a
         // bit messy.
-        if constexpr ( represents_composite ) {
+        if constexpr ( child_type == ChildType::Composite ) {
           addToFourVector_Composite( vertexpos, p4, p4cov, gainmatrix );
-        } else {
+        } else if constexpr ( child_type == ChildType::Track ) {
           addToFourVector_Track( vertexpos, p4, p4cov, gainmatrix );
+        } else {
         }
       }
 
@@ -456,8 +462,8 @@ namespace Sel::Fitters {
                                 child_t const& child ) {
         assert( type == ChildType::TrackWithVelo || type == ChildType::TrackWithoutVelo );
         using VertexTrack =
-            VertexTraj<child_t, false /* flag this represents a track not a weakly-decaying composite */>;
-        return VertexTrack{child, child.state( LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam )};
+            VertexTraj<child_t, ChildType::Track /* flag this represents a track not a weakly-decaying composite */>;
+        return VertexTrack{child, child.state( child_t::StateLocation::ClosestToBeam )};
       }
     };
 
@@ -510,7 +516,8 @@ namespace Sel::Fitters {
                                 child_t const& child ) {
         assert( type == ChildType::Composite || type == ChildType::Resonance );
         using VertexComposite =
-            VertexTraj<child_t, true /* flag this represents a weakly-decaying composite and not a track */>;
+            VertexTraj<child_t,
+                       ChildType::Composite /* flag this represents a weakly-decaying composite and not a track */>;
         return VertexComposite{child, stateVectorFromComposite( child )};
       }
     };
@@ -533,9 +540,11 @@ namespace Sel::Fitters {
       template <typename mask_v, typename float_v>
       static auto prepareChild( [[maybe_unused]] VecN<3, float_v> const& vertex_position,
                                 [[maybe_unused]] mask_v const& vertex_position_mask, [[maybe_unused]] ChildType type,
-                                [[maybe_unused]] child_t const& child ) {
+                                child_t const& child ) {
         assert( type == ChildType::Other );
-        return bool{};
+        using VertexOther =
+            VertexTraj<child_t, ChildType::Other /* flag this represents a track not a weakly-decaying composite */>;
+        return VertexOther{child, stateVectorFromOther( child )};
       }
     };
 
-- 
GitLab


From 0933d2cbed568b2534f7e29e4445cf0caa285269 Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Wed, 24 Nov 2021 11:39:54 +0000
Subject: [PATCH 12/26] Include Neutral as a ChildType in the vertex fitter

---
 .../include/VertexFit/ParticleVertexFitter.h  | 107 +++++++++++++++---
 1 file changed, 92 insertions(+), 15 deletions(-)

diff --git a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
index 1e2da8abae4..66719c52a4b 100644
--- a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
+++ b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
@@ -118,16 +118,16 @@ namespace Sel::Fitters {
     enum struct ChildType : unsigned char {
       TrackWithVelo = 0,
       TrackWithoutVelo,
-      Track,
       Composite,
       Resonance,
+      Neutral,
       Other,
       NumTypes
     };
 
     template <ChildType C>
     static constexpr bool has_full_state_v = ( C == ChildType::Composite || C == ChildType::Resonance ||
-                                               C == ChildType::Other );
+                                               C == ChildType::Neutral || C == ChildType::Other );
 
     /** Indirectly stolen from TrackVertexUtils::poca
      */
@@ -259,7 +259,7 @@ namespace Sel::Fitters {
         // bit messy.
         if constexpr ( child_type == ChildType::Composite ) {
           addToFourVector_Composite( vertexpos, p4, p4cov, gainmatrix );
-        } else if constexpr ( child_type == ChildType::Track ) {
+        } else if constexpr ( child_type == ChildType::TrackWithVelo || child_type == ChildType::TrackWithoutVelo ) {
           addToFourVector_Track( vertexpos, p4, p4cov, gainmatrix );
         } else {
         }
@@ -423,7 +423,7 @@ namespace Sel::Fitters {
      *  @todo Handle TrackWithoutVelo.
      */
     template <typename child_t>
-    struct TrackChildHelper {
+    struct TrackWithVeloChildHelper {
       /** Determine if the given `child_t` type represents a [chunk of] tracks.
        *  @todo Make this more sophisticated and follow changes in the event
        *        model.
@@ -460,9 +460,60 @@ namespace Sel::Fitters {
       static auto prepareChild( [[maybe_unused]] VecN<3, float_v> const& vertex_position,
                                 [[maybe_unused]] mask_v const& vertex_position_mask, [[maybe_unused]] ChildType type,
                                 child_t const& child ) {
-        assert( type == ChildType::TrackWithVelo || type == ChildType::TrackWithoutVelo );
+        assert( type == ChildType::TrackWithVelo );
         using VertexTrack =
-            VertexTraj<child_t, ChildType::Track /* flag this represents a track not a weakly-decaying composite */>;
+            VertexTraj<child_t,
+                       ChildType::TrackWithVelo /* flag this represents a track not a weakly-decaying composite */>;
+        return VertexTrack{child, child.state( child_t::StateLocation::ClosestToBeam )};
+      }
+    };
+
+    /** Helper type for identifying and handling TrackWithVelo and
+     *  TrackWithoutVelo children.
+     *  @todo Handle TrackWithoutVelo.
+     */
+    template <typename child_t>
+    struct TrackWithoutVeloChildHelper {
+      /** Determine if the given `child_t` type represents a [chunk of] tracks.
+       *  @todo Make this more sophisticated and follow changes in the event
+       *        model.
+       */
+      constexpr static bool value = Sel::Utils::has_tracklike_API<child_t>;
+
+      /** Given that `value` was `true`, derive (at runtime) the `ChildType` of
+       *  `child_t`. This basically means "is it a long or downstream track".
+       *  @todo Implement check and maybe return `ChildType::TrackWithoutVelo`.
+       *  @todo Clearly state what is supported when `child` is a vector type.
+       *        For example, if it is a proxy into a container that supports
+       *        heterogeneous row types (e.g. some rows are long, some rows are
+       *        downstream) then there is not a single correct return value we
+       *        can give here. Is that a runtime error? Should the return type
+       *        be changed so that can be supported?
+       */
+      static ChildType deriveType( LHCb::IParticlePropertySvc const&, [[maybe_unused]] child_t const& child ) {
+        return ChildType::TrackWithoutVelo;
+      }
+
+      /** Given that `value` was `true`, and that `type` was returned by our
+       *  `deriveType` function, return the narrowest-possible std::variant of
+       *  types implementing the API used in the vertex fit. For tracks there
+       *  should just be one return type (i.e. a variant of size 1), but there
+       *  will be some runtime branching on the `type` to figure out how to get
+       *  the relevant state.
+       *  @todo Use `vertex_position` and its validity mask
+       *        `vertex_position_mask` to determine if a sophisticated
+       *        extrapolation is needed.
+       *  @todo Implement this logic, don't just blindly call the
+       *        `state(StateLocation::ClosestToBeam)` accessor. Remove [[maybe_unused]].
+       */
+      template <typename mask_v, typename float_v>
+      static auto prepareChild( [[maybe_unused]] VecN<3, float_v> const& vertex_position,
+                                [[maybe_unused]] mask_v const& vertex_position_mask, [[maybe_unused]] ChildType type,
+                                child_t const& child ) {
+        assert( type == ChildType::TrackWithoutVelo );
+        using VertexTrack =
+            VertexTraj<child_t,
+                       ChildType::TrackWithoutVelo /* flag this represents a track not a weakly-decaying composite */>;
         return VertexTrack{child, child.state( child_t::StateLocation::ClosestToBeam )};
       }
     };
@@ -522,7 +573,33 @@ namespace Sel::Fitters {
       }
     };
 
-    /** Helper type for the default `Other` case.
+    /** Helper type for the default `Neutral` case.
+     */
+    template <typename child_t>
+    struct NeutralChildHelper {
+      /** This helper just tries generic logic, it should be tried last and it
+       *  accepts everything. If it doesn't work, a compile error is the
+       *  correct outcome.
+       */
+      constexpr static bool value = Sel::Utils::has_neutrals_API<child_t>;
+      /** Catch-all, don't expect other logic here.
+       */
+      static ChildType deriveType( LHCb::IParticlePropertySvc const&, child_t const& ) { return ChildType::Neutral; }
+      /** @todo Need to implement the helper class that handles generic
+       *        children and return an instance of it. As a 1-member variant.
+       */
+      template <typename mask_v, typename float_v>
+      static auto prepareChild( [[maybe_unused]] VecN<3, float_v> const& vertex_position,
+                                [[maybe_unused]] mask_v const& vertex_position_mask, [[maybe_unused]] ChildType type,
+                                child_t const& child ) {
+        assert( type == ChildType::Neutral );
+        using VertexNeutral =
+            VertexTraj<child_t, ChildType::Neutral /* flag this represents a track not a weakly-decaying composite */>;
+        return VertexNeutral{child, stateVectorFromNeutral( child )};
+      }
+    };
+
+    /** Helper type for the default `Neutral` case.
      */
     template <typename child_t>
     struct LeftoverHelper {
@@ -533,7 +610,7 @@ namespace Sel::Fitters {
       constexpr static bool value = true;
       /** Catch-all, don't expect other logic here.
        */
-      static ChildType deriveType( LHCb::IParticlePropertySvc const&, child_t const& ) { return ChildType::Other; }
+      static ChildType deriveType( LHCb::IParticlePropertySvc const&, child_t const& ) { return ChildType::Neutral; }
       /** @todo Need to implement the helper class that handles generic
        *        children and return an instance of it. As a 1-member variant.
        */
@@ -542,9 +619,7 @@ namespace Sel::Fitters {
                                 [[maybe_unused]] mask_v const& vertex_position_mask, [[maybe_unused]] ChildType type,
                                 child_t const& child ) {
         assert( type == ChildType::Other );
-        using VertexOther =
-            VertexTraj<child_t, ChildType::Other /* flag this represents a track not a weakly-decaying composite */>;
-        return VertexOther{child, stateVectorFromOther( child )};
+        return bool{};
       }
     };
 
@@ -553,7 +628,8 @@ namespace Sel::Fitters {
      */
     template <typename child_t>
     using all_helper_types =
-        boost::mp11::mp_list<TrackChildHelper<child_t>, CompositeChildHelper<child_t>, LeftoverHelper<child_t>>;
+        boost::mp11::mp_list<TrackWithVeloChildHelper<child_t>, TrackWithoutVeloChildHelper<child_t>,
+                             NeutralChildHelper<child_t>, CompositeChildHelper<child_t>, LeftoverHelper<child_t>>;
 
     /** Helper types that match (::value == true) for a given child type.
      */
@@ -721,7 +797,7 @@ namespace Sel::Fitters {
         oss << "Could not initialise Sel::Fitters::ParticleVertex: #TrackWithVelo = "
             << count( ChildType::TrackWithVelo ) << ", #TrackWithoutVelo = " << count( ChildType::TrackWithoutVelo )
             << ", #Composite = " << count( ChildType::Composite ) << ", #Resonance = " << count( ChildType::Resonance )
-            << ", #Other = " << count( ChildType::Other );
+            << ", #Neutral = " << count( ChildType::Neutral ) << ", #Other = " << count( ChildType::Other );
         throw exception( oss.str() );
       }
     }
@@ -757,7 +833,8 @@ namespace Sel::Fitters {
       // - TrackWithoutVelo: tracks without velo hits (only downstream, hopefully)
       // - Resonance: composites with ctau < 1micron
       // - Composite: other composites
-      // - Other: everything else, e.g. photons, pi0, jets
+      // - Neutral: neutral objects, e.g. photons, pi0
+      // - Other: everything else, e.g. jets
       // Not all of these classifications can be made at compile time based on
       // the template pack `child_ts...`. If the input children are themselves
       // variants then `deriveType` is used as a visitor, increasing the amount
@@ -780,7 +857,7 @@ namespace Sel::Fitters {
 
       // Now estimate everything that goes into the vertex, this array will be
       // filled with variants; the different variant types encode the different
-      // handling required for TrackWith[out]Velo/Composite/Resonance/Other.
+      // handling required for TrackWith[out]Velo/Composite/Resonance/Neutral/Other.
       std::tuple extrapolated_children{
           prepareChild( vertex_position, prelim_fit_mask, types[IChild], children.template get<IChild>() )...};
 
-- 
GitLab


From f6b9e8a4f4b897971a761d68f767f740ff866976 Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Fri, 26 Nov 2021 11:46:41 +0000
Subject: [PATCH 13/26] Add a NeutralBasics particle maker

---
 Phys/ParticleMaker/src/NeutralMakers.cpp | 43 ++++++++++++++++++++++++
 1 file changed, 43 insertions(+)

diff --git a/Phys/ParticleMaker/src/NeutralMakers.cpp b/Phys/ParticleMaker/src/NeutralMakers.cpp
index 43d39bdea1a..346f16657a8 100644
--- a/Phys/ParticleMaker/src/NeutralMakers.cpp
+++ b/Phys/ParticleMaker/src/NeutralMakers.cpp
@@ -14,6 +14,7 @@
 #include "GaudiAlg/Transformer.h"
 
 #include "Event/Particle.h"
+#include "Event/Particle_v2.h"
 #include "Event/ProtoParticle.h"
 #include "Event/RecVertex.h"
 
@@ -634,4 +635,46 @@ namespace LHCb::Phys::ParticleMakers {
     debug() << " Skipped : " << m_SkipPi0sCounter << endmsg;
     debug() << "--------------------" << endmsg;
   }
+
+  using neutral_basics_maker_output_t =
+      std::tuple<LHCb::Event::NeutralBasics, std::unique_ptr<LHCb::Event::CaloHypothesesWithDirection>>;
+
+  /** @class NeutralBasicsMaker
+   *
+   * @brief Create a container of neutral objects from the output of the reconstruction sequence and a set of vertices
+   *
+   *
+   */
+  class NeutralBasicsMaker : public Gaudi::Functional::Transformer<neutral_basics_maker_output_t(
+                                 LHCb::UniqueIDGenerator const&, LHCb::Event::Calo::v2::Hypotheses const&,
+                                 LHCb::Event::v2::RecVertices const& )> {
+
+  public:
+    NeutralBasicsMaker( const std::string& name, ISvcLocator* pSvcLocator )
+        : Transformer( name, pSvcLocator,
+                       {KeyValue{"InputUniqueIDGenerator", LHCb::UniqueIDGeneratorLocation::Default},
+                        KeyValue{"InputCaloObjects", LHCb::Event::Calo::v2::HypothesesLocation::Default},
+                        KeyValue{"InputPrimaryVertices", LHCb::Event::v2::RecVertexLocation::Primary}},
+                       {KeyValue{"OutputParticles", ""}} ) {}
+
+    neutral_basics_maker_output_t operator()( LHCb::UniqueIDGenerator const&           unique_id_gen,
+                                              LHCb::Event::Calo::v2::Hypotheses const& calo_hypos,
+                                              LHCb::Event::v2::RecVertices const& primary_vertices ) const override {
+
+      auto zn                        = Zipping::generateZipIdentifier();
+      auto calo_hypos_with_direction = std::make_unique<LHCb::Event::CaloHypothesesWithDirection>( unique_id_gen, zn );
+
+      calo_hypos_with_direction->reserve( calo_hypos.size() * primary_vertices.size() );
+
+      // here we are creating one hypothesis per primary vertex
+      for ( auto const& ch : calo_hypos )
+        for ( auto const& vtx : primary_vertices ) calo_hypos_with_direction->emplace_back( unique_id_gen, ch, vtx );
+
+      return {LHCb::Pr::make_zip( std::as_const( *calo_hypos_with_direction ) ),
+              std::move( calo_hypos_with_direction )};
+    }
+  };
+
+  DECLARE_COMPONENT( NeutralBasicsMaker )
+
 } // namespace LHCb::Phys::ParticleMakers
-- 
GitLab


From a1020cc4760634d9e4709ea8ebb911786090f232 Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Mon, 29 Nov 2021 14:03:25 +0000
Subject: [PATCH 14/26] Fix NeutralBasicsMaker

---
 Phys/ParticleMaker/src/NeutralMakers.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Phys/ParticleMaker/src/NeutralMakers.cpp b/Phys/ParticleMaker/src/NeutralMakers.cpp
index 346f16657a8..e918e242897 100644
--- a/Phys/ParticleMaker/src/NeutralMakers.cpp
+++ b/Phys/ParticleMaker/src/NeutralMakers.cpp
@@ -655,7 +655,7 @@ namespace LHCb::Phys::ParticleMakers {
                        {KeyValue{"InputUniqueIDGenerator", LHCb::UniqueIDGeneratorLocation::Default},
                         KeyValue{"InputCaloObjects", LHCb::Event::Calo::v2::HypothesesLocation::Default},
                         KeyValue{"InputPrimaryVertices", LHCb::Event::v2::RecVertexLocation::Primary}},
-                       {KeyValue{"OutputParticles", ""}} ) {}
+                       {KeyValue{"Particles", ""}} ) {}
 
     neutral_basics_maker_output_t operator()( LHCb::UniqueIDGenerator const&           unique_id_gen,
                                               LHCb::Event::Calo::v2::Hypotheses const& calo_hypos,
@@ -675,6 +675,6 @@ namespace LHCb::Phys::ParticleMakers {
     }
   };
 
-  DECLARE_COMPONENT( NeutralBasicsMaker )
+  DECLARE_COMPONENT( LHCb::Phys::ParticleMakers::NeutralBasicsMaker )
 
 } // namespace LHCb::Phys::ParticleMakers
-- 
GitLab


From 53ca06d58c97e7040f04a8afebd42c8555048259 Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Wed, 1 Dec 2021 16:31:43 +0000
Subject: [PATCH 15/26] Define a NeutralBasics filter using Pr::Filter

---
 Phys/ParticleCombiners/src/Filter.cpp                   | 2 ++
 Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h | 7 +++++++
 2 files changed, 9 insertions(+)

diff --git a/Phys/ParticleCombiners/src/Filter.cpp b/Phys/ParticleCombiners/src/Filter.cpp
index d87b7eee794..db12706628b 100644
--- a/Phys/ParticleCombiners/src/Filter.cpp
+++ b/Phys/ParticleCombiners/src/Filter.cpp
@@ -28,4 +28,6 @@ namespace LHCb {
 
 DECLARE_COMPONENT_WITH_ID( Pr::Filter<LHCb::Event::ChargedBasics>, "ChargedBasicsFilter" )
 
+DECLARE_COMPONENT_WITH_ID( Pr::Filter<LHCb::Event::NeutralBasics>, "NeutralBasicsFilter" )
+
 DECLARE_COMPONENT_WITH_ID( Pr::Filter<LHCb::Particle::Range>, "ParticleRangeFilter" )
diff --git a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
index 66719c52a4b..dedc95c348e 100644
--- a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
+++ b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
@@ -261,10 +261,17 @@ namespace Sel::Fitters {
           addToFourVector_Composite( vertexpos, p4, p4cov, gainmatrix );
         } else if constexpr ( child_type == ChildType::TrackWithVelo || child_type == ChildType::TrackWithoutVelo ) {
           addToFourVector_Track( vertexpos, p4, p4cov, gainmatrix );
+        } else if constexpr ( child_type == ChildType::Neutral ) {
+          addToFourVector_Neutral( vertexpos, p4, p4cov, gainmatrix );
         } else {
         }
       }
 
+      /** @todo FIXME correctly process the errors
+       */
+      void addToFourVector_Neutral( VecN<3, float_v> const&, VecN<4, float_v>&, SymNxN<4, float_v>&,
+                                    Matrix<4, 3, float_v>& ) const {}
+
       /** @todo FIXME implement logic for bremsstrahlung photons.
        */
       void addToFourVector_Track( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>& p4cov,
-- 
GitLab


From fefdc0dbd03d69e57b19b5bd61bfdb1d4eac6986 Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Thu, 2 Dec 2021 19:34:39 +0000
Subject: [PATCH 16/26] Add a basic composite + neutral particle combiner

---
 .../src/ThOrCombiner_Composites.cpp                |  4 ++++
 Phys/ParticleMaker/src/NeutralMakers.cpp           | 14 +++++++-------
 2 files changed, 11 insertions(+), 7 deletions(-)

diff --git a/Phys/ParticleCombiners/src/ThOrCombiner_Composites.cpp b/Phys/ParticleCombiners/src/ThOrCombiner_Composites.cpp
index b90e62625d7..57808bdc750 100644
--- a/Phys/ParticleCombiners/src/ThOrCombiner_Composites.cpp
+++ b/Phys/ParticleCombiners/src/ThOrCombiner_Composites.cpp
@@ -19,6 +19,10 @@ using ThOrCombiner__Composites3ChargedBasics =
     ThOr::CombinerBest<LHCb::Event::Composites, LHCb::Event::ChargedBasics, LHCb::Event::ChargedBasics,
                        LHCb::Event::ChargedBasics>;
 
+using ThOrCombiner__CompositesNeutralBasics = ThOr::CombinerBest<LHCb::Event::Composites, LHCb::Event::NeutralBasics>;
+
 DECLARE_COMPONENT_WITH_ID( ThOrCombiner__CompositesChargedBasics, "ThOrCombiner__CompositesChargedBasics" )
 DECLARE_COMPONENT_WITH_ID( ThOrCombiner__Composites2ChargedBasics, "ThOrCombiner__Composites2ChargedBasics" )
 DECLARE_COMPONENT_WITH_ID( ThOrCombiner__Composites3ChargedBasics, "ThOrCombiner__Composites3ChargedBasics" )
+
+DECLARE_COMPONENT_WITH_ID( ThOrCombiner__CompositesNeutralBasics, "ThOrCombiner__CompositesNeutralBasics" )
diff --git a/Phys/ParticleMaker/src/NeutralMakers.cpp b/Phys/ParticleMaker/src/NeutralMakers.cpp
index e918e242897..3d95931b5f1 100644
--- a/Phys/ParticleMaker/src/NeutralMakers.cpp
+++ b/Phys/ParticleMaker/src/NeutralMakers.cpp
@@ -637,7 +637,7 @@ namespace LHCb::Phys::ParticleMakers {
   }
 
   using neutral_basics_maker_output_t =
-      std::tuple<LHCb::Event::NeutralBasics, std::unique_ptr<LHCb::Event::CaloHypothesesWithDirection>>;
+      std::tuple<LHCb::Event::NeutralBasics, std::unique_ptr<LHCb::Event::CaloHypothesesWithDirection>>> ;
 
   /** @class NeutralBasicsMaker
    *
@@ -645,17 +645,17 @@ namespace LHCb::Phys::ParticleMakers {
    *
    *
    */
-  class NeutralBasicsMaker : public Gaudi::Functional::Transformer<neutral_basics_maker_output_t(
+  class NeutralBasicsMaker : public Gaudi::Functional::MultiTransformer<neutral_basics_maker_output_t(
                                  LHCb::UniqueIDGenerator const&, LHCb::Event::Calo::v2::Hypotheses const&,
                                  LHCb::Event::v2::RecVertices const& )> {
 
   public:
     NeutralBasicsMaker( const std::string& name, ISvcLocator* pSvcLocator )
-        : Transformer( name, pSvcLocator,
-                       {KeyValue{"InputUniqueIDGenerator", LHCb::UniqueIDGeneratorLocation::Default},
-                        KeyValue{"InputCaloObjects", LHCb::Event::Calo::v2::HypothesesLocation::Default},
-                        KeyValue{"InputPrimaryVertices", LHCb::Event::v2::RecVertexLocation::Primary}},
-                       {KeyValue{"Particles", ""}} ) {}
+        : MultiTransformer( name, pSvcLocator,
+                            {KeyValue{"InputUniqueIDGenerator", LHCb::UniqueIDGeneratorLocation::Default},
+                             KeyValue{"InputCaloObjects", LHCb::Event::Calo::v2::HypothesesLocation::Default},
+                             KeyValue{"InputPrimaryVertices", LHCb::Event::v2::RecVertexLocation::Primary}},
+                            {KeyValue{"Particles", ""}, KeyValue{"CaloHypotheses", ""}} ) {}
 
     neutral_basics_maker_output_t operator()( LHCb::UniqueIDGenerator const&           unique_id_gen,
                                               LHCb::Event::Calo::v2::Hypotheses const& calo_hypos,
-- 
GitLab


From cef33116009660d349867ee3f7f35b6da678c81e Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Fri, 3 Dec 2021 10:30:45 +0000
Subject: [PATCH 17/26] Use CaloHypothesesWithDirection as NeutralBasics

---
 Phys/ParticleMaker/src/NeutralMakers.cpp | 32 +++++++++++-------------
 1 file changed, 14 insertions(+), 18 deletions(-)

diff --git a/Phys/ParticleMaker/src/NeutralMakers.cpp b/Phys/ParticleMaker/src/NeutralMakers.cpp
index 3d95931b5f1..80558d58f7e 100644
--- a/Phys/ParticleMaker/src/NeutralMakers.cpp
+++ b/Phys/ParticleMaker/src/NeutralMakers.cpp
@@ -636,42 +636,38 @@ namespace LHCb::Phys::ParticleMakers {
     debug() << "--------------------" << endmsg;
   }
 
-  using neutral_basics_maker_output_t =
-      std::tuple<LHCb::Event::NeutralBasics, std::unique_ptr<LHCb::Event::CaloHypothesesWithDirection>>> ;
-
   /** @class NeutralBasicsMaker
    *
    * @brief Create a container of neutral objects from the output of the reconstruction sequence and a set of vertices
    *
    *
    */
-  class NeutralBasicsMaker : public Gaudi::Functional::MultiTransformer<neutral_basics_maker_output_t(
+  class NeutralBasicsMaker : public Gaudi::Functional::Transformer<LHCb::Event::NeutralBasics(
                                  LHCb::UniqueIDGenerator const&, LHCb::Event::Calo::v2::Hypotheses const&,
                                  LHCb::Event::v2::RecVertices const& )> {
 
   public:
     NeutralBasicsMaker( const std::string& name, ISvcLocator* pSvcLocator )
-        : MultiTransformer( name, pSvcLocator,
-                            {KeyValue{"InputUniqueIDGenerator", LHCb::UniqueIDGeneratorLocation::Default},
-                             KeyValue{"InputCaloObjects", LHCb::Event::Calo::v2::HypothesesLocation::Default},
-                             KeyValue{"InputPrimaryVertices", LHCb::Event::v2::RecVertexLocation::Primary}},
-                            {KeyValue{"Particles", ""}, KeyValue{"CaloHypotheses", ""}} ) {}
+        : Transformer( name, pSvcLocator,
+                       {KeyValue{"InputUniqueIDGenerator", LHCb::UniqueIDGeneratorLocation::Default},
+                        KeyValue{"InputCaloObjects", LHCb::Event::Calo::v2::HypothesesLocation::Default},
+                        KeyValue{"InputPrimaryVertices", LHCb::Event::v2::RecVertexLocation::Primary}},
+                       {KeyValue{"Particles", ""}} ) {}
 
-    neutral_basics_maker_output_t operator()( LHCb::UniqueIDGenerator const&           unique_id_gen,
-                                              LHCb::Event::Calo::v2::Hypotheses const& calo_hypos,
-                                              LHCb::Event::v2::RecVertices const& primary_vertices ) const override {
+    LHCb::Event::NeutralBasics operator()( LHCb::UniqueIDGenerator const&           unique_id_gen,
+                                           LHCb::Event::Calo::v2::Hypotheses const& calo_hypos,
+                                           LHCb::Event::v2::RecVertices const&      primary_vertices ) const override {
 
-      auto zn                        = Zipping::generateZipIdentifier();
-      auto calo_hypos_with_direction = std::make_unique<LHCb::Event::CaloHypothesesWithDirection>( unique_id_gen, zn );
+      auto                       zn = Zipping::generateZipIdentifier();
+      LHCb::Event::NeutralBasics neutrals( unique_id_gen, zn );
 
-      calo_hypos_with_direction->reserve( calo_hypos.size() * primary_vertices.size() );
+      neutrals.reserve( calo_hypos.size() * primary_vertices.size() );
 
       // here we are creating one hypothesis per primary vertex
       for ( auto const& ch : calo_hypos )
-        for ( auto const& vtx : primary_vertices ) calo_hypos_with_direction->emplace_back( unique_id_gen, ch, vtx );
+        for ( auto const& vtx : primary_vertices ) neutrals.emplace_back( unique_id_gen, ch, vtx );
 
-      return {LHCb::Pr::make_zip( std::as_const( *calo_hypos_with_direction ) ),
-              std::move( calo_hypos_with_direction )};
+      return neutrals;
     }
   };
 
-- 
GitLab


From 837fe1bf77951c9f36b74d9087bd2cac43d96a39 Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Mon, 6 Dec 2021 11:18:25 +0000
Subject: [PATCH 18/26] Allow to choose whether to use all the vertices or the
 average to build the neutrals

---
 Phys/ParticleMaker/src/NeutralMakers.cpp | 31 +++++++++++++++++++-----
 1 file changed, 25 insertions(+), 6 deletions(-)

diff --git a/Phys/ParticleMaker/src/NeutralMakers.cpp b/Phys/ParticleMaker/src/NeutralMakers.cpp
index 80558d58f7e..62243130782 100644
--- a/Phys/ParticleMaker/src/NeutralMakers.cpp
+++ b/Phys/ParticleMaker/src/NeutralMakers.cpp
@@ -17,6 +17,7 @@
 #include "Event/Particle_v2.h"
 #include "Event/ProtoParticle.h"
 #include "Event/RecVertex.h"
+#include "Event/RecVertex_v2.h"
 
 #include "CaloFutureUtils/CaloMomentum.h"
 #include "Kernel/IParticlePropertySvc.h"
@@ -658,17 +659,35 @@ namespace LHCb::Phys::ParticleMakers {
                                            LHCb::Event::Calo::v2::Hypotheses const& calo_hypos,
                                            LHCb::Event::v2::RecVertices const&      primary_vertices ) const override {
 
-      auto                       zn = Zipping::generateZipIdentifier();
-      LHCb::Event::NeutralBasics neutrals( unique_id_gen, zn );
+      LHCb::Event::NeutralBasics neutrals( unique_id_gen );
 
-      neutrals.reserve( calo_hypos.size() * primary_vertices.size() );
+      if ( m_AveragePVs.value() ) {
+        // Calculate the average the primary vertex positions
+        neutrals.reserve( calo_hypos.size() );
 
-      // here we are creating one hypothesis per primary vertex
-      for ( auto const& ch : calo_hypos )
-        for ( auto const& vtx : primary_vertices ) neutrals.emplace_back( unique_id_gen, ch, vtx );
+        auto avg = std::accumulate( primary_vertices.cbegin(), primary_vertices.cend(), Gaudi::XYZPoint{0.f, 0.f, 0.f},
+                                    [&primary_vertices]( auto const& avg, auto const& vtx ) -> Gaudi::XYZPoint {
+                                      return {avg.x() + vtx.position().x() / primary_vertices.size(),
+                                              avg.y() + vtx.position().y() / primary_vertices.size(),
+                                              avg.z() + vtx.position().z() / primary_vertices.size()};
+                                    } );
+
+        for ( auto const& ch : calo_hypos ) neutrals.emplace_back( unique_id_gen, ch, avg );
+      } else {
+        // Create one particle per primary vertex and ECAL hypothesis
+        neutrals.reserve( calo_hypos.size() * primary_vertices.size() );
+
+        for ( auto const& ch : calo_hypos )
+          for ( auto const& vtx : primary_vertices ) neutrals.emplace_back( unique_id_gen, ch, vtx.position() );
+      }
 
       return neutrals;
     }
+
+  private:
+    Gaudi::Property<bool> m_AveragePVs{this, "AveragePVs", true,
+                                       "If set to \"true\"; uses the average of the primary vertices in order to "
+                                       "determine the pointing direction of the neutral particles"};
   };
 
   DECLARE_COMPONENT( LHCb::Phys::ParticleMakers::NeutralBasicsMaker )
-- 
GitLab


From 3d2312d646c1230e40c1eff670345f0c15445035 Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Tue, 7 Dec 2021 16:58:01 +0000
Subject: [PATCH 19/26] Add the four-momenta of the neutral objects corrected
 by the vertex position to the result of the combination

---
 .../include/VertexFit/ParticleVertexFitter.h  | 24 ++++++++++++++++---
 1 file changed, 21 insertions(+), 3 deletions(-)

diff --git a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
index dedc95c348e..0502f8e744d 100644
--- a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
+++ b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
@@ -267,10 +267,28 @@ namespace Sel::Fitters {
         }
       }
 
-      /** @todo FIXME correctly process the errors
+      /** Add the momenta of a neutral particle to the result of a combination
+
+          Update the information of the momenta of the neutral objects based on the vertex of
+          the decaying particle. It is crucial that we add the neutral objects after
+          determining the decay vertex (i.e. at the end of the combination).
+
+          @todo FIXME correctly process the covariance matrix
        */
-      void addToFourVector_Neutral( VecN<3, float_v> const&, VecN<4, float_v>&, SymNxN<4, float_v>&,
-                                    Matrix<4, 3, float_v>& ) const {}
+      void addToFourVector_Neutral( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>&,
+                                    Matrix<4, 3, float_v>& ) const {
+        // TODO: make this a global function that modifies the neutral objects at the
+        // persistency stage
+        auto dx = m_state.x() - vertexpos( 0 );
+        auto dy = m_state.y() - vertexpos( 1 );
+        auto dz = m_state.z() - vertexpos( 2 );
+
+        auto d = sqrt( dx * dx + dy * dy + dz * dz );
+
+        auto p = m_particle.p();
+
+        p4 = VecN<4, float_v>{p * dx / d, p * dy / d, p * dz / d, m_particle.energy()};
+      }
 
       /** @todo FIXME implement logic for bremsstrahlung photons.
        */
-- 
GitLab


From fd29320e8c877f16e0e055e722b94ff098762a90 Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Mon, 17 Jan 2022 12:24:29 +0000
Subject: [PATCH 20/26] Fixes due to the new zipping infrastructure

---
 Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
index 0502f8e744d..6ca5b104046 100644
--- a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
+++ b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
@@ -489,7 +489,7 @@ namespace Sel::Fitters {
         using VertexTrack =
             VertexTraj<child_t,
                        ChildType::TrackWithVelo /* flag this represents a track not a weakly-decaying composite */>;
-        return VertexTrack{child, child.state( child_t::StateLocation::ClosestToBeam )};
+        return VertexTrack{child, child.state( LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam )};
       }
     };
 
@@ -539,7 +539,7 @@ namespace Sel::Fitters {
         using VertexTrack =
             VertexTraj<child_t,
                        ChildType::TrackWithoutVelo /* flag this represents a track not a weakly-decaying composite */>;
-        return VertexTrack{child, child.state( child_t::StateLocation::ClosestToBeam )};
+        return VertexTrack{child, child.state( LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam )};
       }
     };
 
@@ -788,8 +788,8 @@ namespace Sel::Fitters {
                       velotrajs[nvelotrajs++] =
                           StateVector4<float_v>{child.state( LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam )};
                     } else {
-                      throw exception(
-                          "TrackWithVelo type didn't have a state( child_t::StateLocation::ClosestToBeam)" );
+                      throw exception( "TrackWithVelo type didn't have a state( "
+                                       "LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam)" );
                     }
                   } else if ( types[i] == ChildType::Composite ) {
                     if constexpr ( Sel::type_traits::has_endVertex_v<child_t> &&
-- 
GitLab


From e63f366248a8ef18e9121d50e215e7d30dcf29eb Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Tue, 18 Jan 2022 14:53:35 +0000
Subject: [PATCH 21/26] Implement missing member functions for
 ParticleCombination

---
 .../SelKernel/include/SelKernel/ParticleCombination.h | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/Phys/SelKernel/include/SelKernel/ParticleCombination.h b/Phys/SelKernel/include/SelKernel/ParticleCombination.h
index 3bd96f8bafe..3b8bad9e038 100644
--- a/Phys/SelKernel/include/SelKernel/ParticleCombination.h
+++ b/Phys/SelKernel/include/SelKernel/ParticleCombination.h
@@ -120,10 +120,19 @@ namespace Sel::detail {
       return transform_reduce( []( auto const& p ) { return p.charge(); }, std::plus<>{} );
     }
 
+    auto mass2() const {
+      auto mom = this->momentum();
+      auto px  = mom( 0 );
+      auto py  = mom( 1 );
+      auto pz  = mom( 2 );
+      auto e   = mom( 3 );
+      return e * e - px * px + py * py + pz * pz;
+    }
+
     auto mass() const {
       using Sel::Utils::mass2;
       using std::sqrt;
-      return sqrt( mass2( this->fourMomentum() ) );
+      return sqrt( mass2( *this ) );
     }
 
     /** Unfortunately this is needed in order for explicit mask functionality
-- 
GitLab


From b203bf2cabee1b977c8ed381f9512b466c5ab71d Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Thu, 31 Mar 2022 10:52:32 +0100
Subject: [PATCH 22/26] Fix compilation issues by defining new members in proxy
 types

---
 .../include/CombKernel/ThOrCombiner.h         |    8 +-
 Phys/SelKernel/include/SelKernel/Utilities.h  |    6 +
 Phys/SelTools/include/SelTools/Utilities.h    |    2 +
 .../include/VertexFit/ParticleVertexFitter.h  | 1116 ++++++++---------
 4 files changed, 508 insertions(+), 624 deletions(-)

diff --git a/Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h b/Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h
index 17af73e8ca9..252db5d0eff 100644
--- a/Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h
+++ b/Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h
@@ -169,11 +169,10 @@ namespace ThOr::detail::Combiner {
 
   template <typename T0, typename T1>
   auto have_overlap( T0 const& p0, T1 const& p1 ) {
-    if constexpr ( Sel::Utils::canBeExtrapolatedDownstream_v<T0> && Sel::Utils::canBeExtrapolatedDownstream_v<T1> ) {
+    if constexpr ( Sel::Utils::isBasicParticle_v<T0> && Sel::Utils::isBasicParticle_v<T1> ) {
       // track-like/track-like
       return p0.unique_id() == p1.unique_id();
-    } else if constexpr ( Sel::Utils::canBeExtrapolatedDownstream_v<T0> &&
-                          !Sel::Utils::canBeExtrapolatedDownstream_v<T1> ) {
+    } else if constexpr ( Sel::Utils::isBasicParticle_v<T0> && !Sel::Utils::isBasicParticle_v<T1> ) {
       // track-like/composite
       using mask_t = decltype( p0.unique_id() == p0.unique_id() );
       mask_t mask{false};
@@ -181,8 +180,7 @@ namespace ThOr::detail::Combiner {
         mask = mask | ( p1.descendant_unique_id( i ) == p0.unique_id() );
       }
       return mask;
-    } else if constexpr ( !Sel::Utils::canBeExtrapolatedDownstream_v<T0> &&
-                          Sel::Utils::canBeExtrapolatedDownstream_v<T1> ) {
+    } else if constexpr ( !Sel::Utils::isBasicParticle_v<T0> && Sel::Utils::isBasicParticle_v<T1> ) {
       // composite/track-like
       return have_overlap( p1, p0 );
     } else {
diff --git a/Phys/SelKernel/include/SelKernel/Utilities.h b/Phys/SelKernel/include/SelKernel/Utilities.h
index f1c30e4f6ce..7422842d3a9 100644
--- a/Phys/SelKernel/include/SelKernel/Utilities.h
+++ b/Phys/SelKernel/include/SelKernel/Utilities.h
@@ -294,6 +294,12 @@ namespace Sel::Utils {
   template <typename T>
   constexpr auto canBeExtrapolatedDownstream_v = std::remove_pointer_t<std::decay_t<T>>::canBeExtrapolatedDownstream;
 
+  template <typename T>
+  constexpr auto isBasicParticle_v = std::remove_pointer_t<std::decay_t<T>>::isBasicParticle;
+
+  template <typename T>
+  constexpr auto hasTrack_v = std::remove_pointer_t<std::decay_t<T>>::hasTrack;
+
   template <typename T>
   auto threeMomPosCovMatrix( T const& item ) -> decltype( item.threeMomPosCovMatrix() ) {
     return item.threeMomPosCovMatrix();
diff --git a/Phys/SelTools/include/SelTools/Utilities.h b/Phys/SelTools/include/SelTools/Utilities.h
index 6075bac5262..0dd1ebe108a 100644
--- a/Phys/SelTools/include/SelTools/Utilities.h
+++ b/Phys/SelTools/include/SelTools/Utilities.h
@@ -9,10 +9,12 @@
 * or submit itself to any jurisdiction.                                       *
 \*****************************************************************************/
 #pragma once
+#include "Event/Particle_v2.h"
 #include "LHCbMath/MatVec.h"
 #include "SelKernel/State4.h"
 
 namespace Sel {
+
   /** @fn    stateVectorFromComposite
    *  @brief Convert a particle (3-momentum + position) to 4-state (at some z).
    *
diff --git a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
index b2e88434e79..400448403c8 100644
--- a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
+++ b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
@@ -22,6 +22,473 @@
 
 #include "Gaudi/Algorithm.h"
 
+namespace {
+
+  template <std::size_t N, typename T>
+  using VecN = LHCb::LinAlg::Vec<T, N>;
+
+  template <std::size_t N, typename T>
+  using SymNxN = LHCb::LinAlg::MatSym<T, N>;
+
+  template <std::size_t N, std::size_t M, typename T>
+  using Matrix = LHCb::LinAlg::Mat<T, N, M>;
+
+  template <std::size_t N, typename T>
+  using MatrixNxN = LHCb::LinAlg::Mat<T, N>;
+
+  // Fast utility function for inverting a 2x2 matrix. This just comes
+  // from ROOT/Math/Dinv, but I took out the test on the determinant
+  // since we do not need it.
+  template <typename float_v>
+  SymNxN<2, float_v> Invert2x2( SymNxN<2, float_v> const& rhs ) {
+    SymNxN<2, float_v> ret;
+    auto const         det    = rhs( 0, 0 ) * rhs( 1, 1 ) - rhs( 0, 1 ) * rhs( 0, 1 );
+    auto const         invdet = 1.f / det;
+    // if (det == T(0.)) { return false; }
+    ret( 0, 0 ) = rhs( 1, 1 ) * invdet;
+    ret( 1, 1 ) = rhs( 0, 0 ) * invdet;
+    ret( 0, 1 ) = -1.f * rhs( 0, 1 ) * invdet;
+    return ret;
+  }
+
+  // Different daughter types
+  // TODO look at using the meta enum thing
+  enum struct ChildType : unsigned char {
+    TrackWithVelo = 0,
+    TrackWithoutVelo,
+    Composite,
+    Resonance,
+    Neutral,
+    Other,
+    NumTypes
+  };
+
+  template <ChildType C>
+  static constexpr bool has_full_state_v = ( C == ChildType::Composite || C == ChildType::Resonance ||
+                                             C == ChildType::Neutral || C == ChildType::Other );
+
+  /* TODO look at merging this with the stuff in State4.h
+   */
+  template <typename F>
+  struct MiniState {
+    SymNxN<5, F> cov{};
+    auto&        x() { return m_x; }
+    auto const&  x() const { return m_x; }
+    auto&        y() { return m_y; }
+    auto const&  y() const { return m_y; }
+    auto&        z() { return m_z; }
+    auto const&  z() const { return m_z; }
+    auto&        qOverP() { return m_qOverP; }
+    auto const&  qOverP() const { return m_qOverP; }
+    auto&        tx() { return m_tx; }
+    auto const&  tx() const { return m_tx; }
+    auto&        ty() { return m_ty; }
+    auto const&  ty() const { return m_ty; }
+
+    template <typename T>
+    MiniState( T const& other )
+        : m_x{other.x()}, m_y{other.y()}, m_tx{other.tx()}, m_ty{other.ty()}, m_z{other.z()}, m_qOverP{other.qOverP()} {
+      static_assert( Sel::Utils::canBeExtrapolatedDownstream_v<T> );
+      auto const& other_cov = other.covariance();
+      LHCb::Utils::unwind<0, 5>( [this, &other_cov]( auto i ) {
+        LHCb::Utils::unwind<i, 5>( [this, i, &other_cov]( auto j ) { cov( i, j ) = other_cov( i, j ); } );
+      } );
+    }
+
+    void linearTransportTo( F const& new_z ) {
+      const auto dz  = new_z - m_z;
+      const auto dz2 = dz * dz;
+      m_x += dz * m_tx;
+      m_y += dz * m_ty;
+      cov( 0, 0 ) += dz2 * cov( 2, 2 ) + 2 * dz * cov( 2, 0 );
+      cov( 1, 0 ) += dz2 * cov( 3, 2 ) + dz * ( cov( 3, 0 ) + cov( 2, 1 ) );
+      cov( 2, 0 ) += dz * cov( 2, 2 );
+      cov( 3, 0 ) += dz * cov( 3, 2 );
+      cov( 4, 0 ) += dz * cov( 4, 2 );
+      cov( 1, 1 ) += dz2 * cov( 3, 3 ) + 2 * dz * cov( 3, 1 );
+      cov( 2, 1 ) += dz * cov( 3, 2 );
+      cov( 3, 1 ) += dz * cov( 3, 3 );
+      cov( 4, 1 ) += dz * cov( 4, 3 );
+      m_z = new_z;
+    }
+
+    SymNxN<2, F>    covXX() const { return cov.template sub<SymNxN<2, F>, 0, 0>(); }
+    MatrixNxN<2, F> covXT() const { return cov.template sub<MatrixNxN<2, F>, 2, 0>(); }
+
+  private:
+    F m_x{}, m_y{}, m_tx{}, m_ty{}, m_z{}, m_qOverP{};
+  };
+
+  // Class for daughters with a straight-line trajectory
+  template <typename particle_t, ChildType child_type>
+  struct VertexTraj final {
+    using simd_t  = typename particle_t::simd_t;
+    using int_v   = typename simd_t::int_v;
+    using float_v = typename simd_t::float_v;
+    using state_t = std::conditional_t<has_full_state_v<child_type>, Sel::State4<float_v>, MiniState<float_v>>;
+
+  private:
+    particle_t            m_particle; // TODO maybe just cache some information from this?
+    state_t               m_state;    // TODO evaluate what is really needed
+    VecN<2, float_v>      m_q;        // predicted/fitted slope (tx,ty)
+    SymNxN<2, float_v>    m_G;        // weight matrix of (x,y) of state
+    Matrix<2, 3, float_v> m_A;        // projection matrix for vertex position
+
+  public:
+    VertexTraj( particle_t p, state_t state )
+        : m_particle{std::move( p )}, m_state{std::move( state )}, m_q{m_state.tx(), m_state.ty()} {}
+
+    void project( VecN<3, float_v> const& vertexpos, VecN<3, float_v>& halfDChisqDX, SymNxN<3, float_v>& halfD2ChisqDX2,
+                  float_v& chi2, int_v& ndof ) {
+      // move the state (not sure we should do it this way: maybe better to just cache the z.)
+      m_state.linearTransportTo( vertexpos( 2 ) );
+
+      // compute the weight matrix
+      // TODO Just use invChol?
+      m_G = Invert2x2( m_state.covXX() );
+
+      // compute residual
+      VecN<2, float_v> residual{vertexpos( 0 ) - m_state.x(), vertexpos( 1 ) - m_state.y()};
+
+      // fill the projection matrix: use the fitted momentum!
+      m_A( 0, 0 ) = m_A( 1, 1 ) = 1.f;
+      m_A( 0, 2 )               = -1.f * m_q( 0 );
+      m_A( 1, 2 )               = -1.f * m_q( 1 );
+
+      // I tried to make this faster by writing it out, but H does
+      // not contain sufficiently many zeroes. Better to
+      // parallelise.
+
+      // I am not sure that the compiler realizes that A^T*G is used
+      // more than once here, so I'll substitute the code from
+      // similarity. Note that similarity does not return an expression.
+      // halfD2ChisqDX2 += ROOT::Math::Similarity( ROOT::Math::Transpose(m_A), m_G ) ;
+      // halfDChisqDX   += ( ROOT::Math::Transpose(m_A) * m_G) * res ;
+      auto const ATG = m_A.transpose() * m_G.cast_to_mat();
+      // updating relevant fit results
+      halfD2ChisqDX2 = halfD2ChisqDX2 + ( ATG * m_A ).cast_to_sym();
+      halfDChisqDX   = halfDChisqDX + ATG * residual;
+      chi2           = chi2 + LHCb::LinAlg::similarity( residual, m_G );
+      ndof           = ndof + 2;
+    }
+
+    void updateSlopes( VecN<3, float_v> const& vertexpos ) {
+      // first update the residual. (note the subtle difference with that in project!)
+      float_v const          dz = vertexpos( 2 ) - m_state.z();
+      VecN<2, float_v> const residual{vertexpos( 0 ) - ( m_state.x() + m_state.tx() * dz ),
+                                      vertexpos( 1 ) - ( m_state.y() + m_state.ty() * dz )};
+
+      // get the matrix that is the correlation of (x,y) and (tx,ty,qop)
+      // ROOT::Math::SMatrix<double,2,2> Vba = m_state.covariance().template
+      // Sub<ROOT::Math::SMatrix<double,2,2>>(2,0); compute the corresponding gain matrix (this is WBG in the BFR fit,
+      // but here it is Vba * Vaa^-1) ROOT::Math::SMatrix<double,2,2> K = Vba*m_G ; compute the momentum vector
+      m_q = VecN<2, float_v>{{m_state.tx(), m_state.ty()}} + m_state.covXT().transpose() * m_G.cast_to_mat() * residual;
+    }
+
+    void addToFourVector( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>& p4cov,
+                          Matrix<4, 3, float_v>& gainmatrix ) const {
+      // Despatch is a bit crude, but it'll do...in the original
+      // ParticleVertexFitter this was done with a full specialisation, but
+      // now it would have to be a partial specialisation and it all gets a
+      // bit messy.
+      if constexpr ( child_type == ChildType::Composite ) {
+        addToFourVector_Composite( vertexpos, p4, p4cov, gainmatrix );
+      } else if constexpr ( child_type == ChildType::TrackWithVelo || child_type == ChildType::TrackWithoutVelo ) {
+        addToFourVector_Track( vertexpos, p4, p4cov, gainmatrix );
+      } else if constexpr ( child_type == ChildType::Neutral ) {
+        addToFourVector_Neutral( vertexpos, p4, p4cov, gainmatrix );
+      } else {
+      }
+    }
+
+    /** Add the momenta of a neutral particle to the result of a combination
+
+        Update the information of the momenta of the neutral objects based on the vertex of
+        the decaying particle. It is crucial that we add the neutral objects after
+        determining the decay vertex (i.e. at the end of the combination).
+
+        @todo FIXME correctly process the covariance matrix
+     */
+    void addToFourVector_Neutral( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>&,
+                                  Matrix<4, 3, float_v>& ) const {
+      // TODO: make this a global function that modifies the neutral objects at the
+      // persistency stage
+      auto dx = m_state.x() - vertexpos( 0 );
+      auto dy = m_state.y() - vertexpos( 1 );
+      auto dz = m_state.z() - vertexpos( 2 );
+
+      auto d = sqrt( dx * dx + dy * dy + dz * dz );
+
+      auto p = m_particle.p();
+
+      p4 = VecN<4, float_v>{p * dx / d, p * dy / d, p * dz / d, m_particle.energy()};
+    }
+
+    /** @todo FIXME implement logic for bremsstrahlung photons.
+     */
+    void addToFourVector_Track( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>& p4cov,
+                                Matrix<4, 3, float_v>& gain_matrix ) const {
+      // This specialisation is for TrackWithVelo and TrackWithoutVelo
+      // we first need to update q/p. since we also need the full gain matrix, we could as well redo that part.
+      float_v const          dz = vertexpos( 2 ) - m_state.z();
+      VecN<2, float_v> const residual{vertexpos( 0 ) - ( m_state.x() + m_q( 0 ) * dz ),
+                                      vertexpos( 1 ) - ( m_state.y() + m_q( 1 ) * dz )};
+
+      auto const Vba    = m_state.cov.template sub<Matrix<3, 2, float_v>, 2, 0>();
+      auto       K      = Vba * m_G.cast_to_mat();
+      auto       tx     = m_state.tx() + ( K * residual )( 0 );
+      auto       ty     = m_state.ty() + ( K * residual )( 1 );
+      auto       qOverP = m_state.qOverP() + ( K * residual )( 2 );
+      // update 4-momentum
+      // TODO can we just store this value in the constructor? what else is needed? bremsstrahlung info?
+      auto const mass = m_particle.mass();
+      auto const p    = 1.f / abs( qOverP );
+      auto const n2   = 1.f + tx * tx + ty * ty;
+      auto const n    = sqrt( n2 );
+      auto const pz   = p / n;
+      auto const px   = tx * pz;
+      auto const py   = ty * pz;
+      auto const E    = sqrt( p * p + mass * mass );
+      p4              = p4 + VecN<4, float_v>{{px, py, pz, E}}; // TODO check whether this gives right results
+                                                                // it was E px py pz before
+
+      Matrix<4, 3, float_v> dp4dmom{}; // jacobi matrix of del(p4)/del(tx,ty,qOverP)
+      auto const            n3    = n2 * n;
+      auto const            invn3 = 1.f / n3;
+      dp4dmom( 0, 0 )             = p * ( 1.f + ty * ty ) * invn3; // dpx/dtx
+      dp4dmom( 0, 1 )             = p * tx * -ty * invn3;          // dpx/dty
+      dp4dmom( 0, 2 )             = -px / qOverP;                  // dpx/dqop
+
+      dp4dmom( 1, 0 ) = p * ty * -tx * invn3;        // dpy/dtx
+      dp4dmom( 1, 1 ) = p * ( 1 + tx * tx ) * invn3; // dpy/dty
+      dp4dmom( 1, 2 ) = -py / qOverP;                // dpy/dqop
+
+      dp4dmom( 2, 0 ) = pz * -tx / n2; // dpz/dtx
+      dp4dmom( 2, 1 ) = pz * -ty / n2; // dpz/dtx
+      dp4dmom( 2, 2 ) = -pz / qOverP;  // dpz/dqop
+
+      dp4dmom( 3, 0 ) = 0.0;                 // dE/dtx
+      dp4dmom( 3, 1 ) = 0.0;                 // dE/dty
+      dp4dmom( 3, 2 ) = p / E * -p / qOverP; // dE/dqop
+
+      auto FK = dp4dmom * K;
+
+      p4cov = p4cov + LHCb::LinAlg::similarity( dp4dmom, m_state.cov.template sub<SymNxN<3, float_v>, 2, 2>() ) -
+              LHCb::LinAlg::similarity( FK, m_state.cov.template sub<SymNxN<2, float_v>, 0, 0>() );
+
+      gain_matrix = gain_matrix + FK * m_A;
+      // // For brem-recovered electrons, we need to do something
+      // // special. So, electrons are an ugly exception here. We could
+      // // also add to q/p instead, which is cheaper, but perhaps even
+      // // more ugly.
+      // if ( particle().particleID().abspid() == 11 ) {
+      //   const double absqop               = std::abs( m_state.qOverP() );
+      //   const double momentumFromParticle = particle().momentum().P();
+      //   // is 1% a reasonable threshold?
+      //   if ( momentumFromParticle * absqop > 1.01 ) {
+      //     const double momentumFromTrack       = 1 / absqop;
+      //     const double momentumError2FromTrack = m_state.covariance()( 4, 4 ) * std::pow(
+      //     momentumFromTrack, 4 ); const double momentumError2FromParticle =
+      //         particle().momCovMatrix()( 3, 3 ) * std::pow( particle().momentum().E() /
+      //         momentumFromParticle, 2 );
+      //     const double bremEnergyCov = momentumError2FromParticle - momentumError2FromTrack;
+      //     const double bremEnergy    = momentumFromParticle - momentumFromTrack;
+      //     // if the correction is unphysical, ignore it.
+      //     if ( bremEnergyCov > 0 ) {
+      //       const auto tx = q( 0 );
+      //       const auto ty = q( 1 );
+      //       auto       t  = std::sqrt( 1 + tx * tx + ty * ty );
+      //       // we could also 'scale' the momentum, but we anyway need the components for the
+      //       jacobian Gaudi::LorentzVector p4brem{bremEnergy * tx / t, bremEnergy * ty / t,
+      //       bremEnergy / t, bremEnergy}; p4 += p4brem; ROOT::Math::SMatrix<double, 4, 1> J;
+      //       J( 0, 0 ) = tx / t;
+      //       J( 1, 0 ) = ty / t;
+      //       J( 2, 0 ) = 1 / t;
+      //       J( 3, 0 ) = 1;
+      //       p4cov += ROOT::Math::Similarity( J, Gaudi::SymMatrix1x1{bremEnergyCov} );
+      //     }
+      //   }
+      // }
+    }
+
+    void addToFourVector_Composite( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>& p4cov,
+                                    Matrix<4, 3, float_v>& gainmatrix ) const {
+      // first need to 'update' the momentum vector. but for that we
+      // first need to 'transport'.
+      using Sel::Utils::fourMomentum;
+      auto          mom = fourMomentum( m_particle ); // 4-momentum
+      float_v const pz  = Z( mom );
+      float_v const tx  = m_q( 0 );
+      float_v const ty  = m_q( 1 );
+
+      // first update the residual
+      using namespace LHCb::LinAlg;
+      using Sel::Utils::referencePoint;
+      auto             particlepos = referencePoint( m_particle );
+      float_v const    dz          = Z( vertexpos ) - particlepos.Z();
+      VecN<2, float_v> res{X( vertexpos ) - ( particlepos.X() + tx * dz ),
+                           Y( vertexpos ) - ( particlepos.Y() + ty * dz )};
+      auto const&      R    = m_state.covXX();
+      auto const&      Rinv = m_G;
+      using Sel::Utils::momCovMatrix;
+      auto const momCov = momCovMatrix( m_particle );
+      using Sel::Utils::momPosCovMatrix;
+      auto const momPosCov = momPosCovMatrix( m_particle );
+
+      // To do this right we need THREE projection matrices for the residual:
+      //    r = A*x_vertex + B*mom_dau + C*x_dau
+      // Note that C=-A.
+      // Lets call the matrix F^T = (B^T C^T). then the uncertainty on the residual is
+      //    R = F V77 F^T where V77 is the 7x7 covariance matrix of the daughter
+      // We will need R every iteration
+      // The correlation matrix between residual and the 7 parameters of the daughter is
+      //    Vr_xmom = V77 F^T
+      // Consequently, the gain matrix for the momentum-decayvertex is
+      //    K72 = V77 F^T R^-1
+      // Now, we do not need the updated parameters for the decay
+      // vertex. If we just need the momentum part, then this reduces to
+      //    K42 = ( V43 * C^T + V44 * B^T ) R^-1
+
+      // The easiest implementation is this ...
+      //  Gaudi::Matrix3x2 CT{ ROOT::Math::Transpose(m_A) } ;
+      //  ROOT::Math::SMatrix<double,4,2> BT ;
+      //  BT(0,0) = +dz/pz;
+      //  BT(2,0) = -dz/pz*tx ;
+      //  BT(1,1) = +dz/pz ;
+      //  BT(2,1) = -dz/pz*ty ;
+      //  ROOT::Math::SMatrix<double,4,2> K42 = ( momPosCov * CT + momCov * BT ) * Rinv ;
+      // but since B and C are largely empty, we better write it out:
+      float_v const         dzopz = dz / pz;
+      Matrix<4, 2, float_v> K42part; // uninitialised here
+      for ( int i = 0; i < 4; ++i ) {
+        K42part( i, 0 ) =
+            momPosCov( i, 0 ) - momPosCov( i, 2 ) * tx + momCov( i, 0 ) * dzopz - momCov( i, 2 ) * ( dzopz * tx );
+        K42part( i, 1 ) =
+            momPosCov( i, 1 ) - momPosCov( i, 2 ) * ty + momCov( i, 1 ) * dzopz - momCov( i, 2 ) * ( dzopz * ty );
+      }
+      Matrix<4, 2, float_v> const K42 = K42part * Rinv.cast_to_mat();
+
+      // update the fourvector
+      auto const deltap4 = K42 * res;
+      p4                 = p4 + mom + deltap4;
+      // to understand this, compare to what is done for the track
+      p4cov      = p4cov + momCov - similarity( K42, R );
+      gainmatrix = gainmatrix + K42 * m_A;
+    }
+  };
+
+  template <bool ParticleCondition, bool TrackCondition, bool ExtrapolationCondition>
+  struct child_type_for_conditions;
+
+  template <>
+  struct child_type_for_conditions<true, true, true> {
+    static constexpr auto value = ChildType::TrackWithVelo;
+  };
+
+  template <>
+  struct child_type_for_conditions<true, false, false> {
+    static constexpr auto value = ChildType::Neutral;
+  };
+
+  template <>
+  struct child_type_for_conditions<false, false, false> {
+    static constexpr auto value = ChildType::Composite;
+  };
+
+  template <class Proxy>
+  struct child_type_for_proxy
+      : child_type_for_conditions<Sel::Utils::isBasicParticle_v<Proxy>, Sel::Utils::hasTrack_v<Proxy>,
+                                  Sel::Utils::canBeExtrapolatedDownstream_v<Proxy>> {};
+
+  template <class Proxy>
+  static constexpr auto child_type_for_proxy_v = child_type_for_proxy<Proxy>::value;
+
+  template <ChildType>
+  struct ChildHelper;
+
+  /** Helper type for identifying and handling TrackWithVelo and
+   *  TrackWithoutVelo children.
+   *  @todo Handle TrackWithoutVelo.
+   */
+  template <>
+  struct ChildHelper<ChildType::TrackWithVelo> {
+
+    static constexpr auto child_type = ChildType::TrackWithVelo;
+
+    /** Given that `value` was `true`, and that `type` was returned by our
+     *  `deriveType` function, return the narrowest-possible std::variant of
+     *  types implementing the API used in the vertex fit. For tracks there
+     *  should just be one return type (i.e. a variant of size 1), but there
+     *  will be some runtime branching on the `type` to figure out how to get
+     *  the relevant state.
+     *  @todo Use `vertex_position` and its validity mask
+     *        `vertex_position_mask` to determine if a sophisticated
+     *        extrapolation is needed.
+     *  @todo Implement this logic, don't just blindly call the
+     *        `state(StateLocation::ClosestToBeam)` accessor. Remove [[maybe_unused]].
+     */
+    template <typename mask_v, typename float_v, typename proxy>
+    static auto prepareChild( [[maybe_unused]] LHCb::LinAlg::Vec<float_v, 3> const& vertex_position,
+                              [[maybe_unused]] mask_v const& vertex_position_mask, proxy const& child ) {
+      using VertexTrack =
+          VertexTraj<proxy, child_type /* flag this represents a track not a weakly-decaying composite */>;
+      return VertexTrack{child, child.state( LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam )};
+    }
+  };
+
+  /** Helper type for identifying and handling Composite and Resonance
+   *  children.
+   *  @todo Handle Resonance, adding a runtime check on PID and another type
+   *        to go in the variant returned by prepareChild.
+   */
+  template <>
+  struct ChildHelper<ChildType::Composite> {
+
+    static constexpr auto child_type = ChildType::Composite;
+
+    /** Given that `value` was `true`, and that `type` was returned by our
+     *  `deriveType` function, return the narrowest-possible std::variant of
+     *  types implementing the API used in the vertex fit. Because weakly
+     *  decaying particles are handled differently in the fit, the return
+     *  type here will have to be a variant of two different types.
+     *
+     *  The original ParticleVertexFitter didn't use the estimated vertex
+     *  position for `Resonance` and `Composite` children, so we don't
+     *  either.
+     *  @todo Add logic for the `ChildType::Resonance` value of `type` and
+     *        add another member to the return type; remove [[maybe_unused]]
+     */
+    template <typename mask_v, typename float_v, typename proxy>
+    static auto prepareChild( LHCb::LinAlg::Vec<float_v, 3> const&, mask_v const&, proxy const& child ) {
+      using VertexComposite =
+          VertexTraj<proxy, child_type /* flag this represents a weakly-decaying composite and not a track */>;
+      return VertexComposite{child, Sel::stateVectorFromComposite( child )};
+    }
+  };
+
+  /** Helper type for the default `Neutral` case.
+   */
+  template <>
+  struct ChildHelper<ChildType::Neutral> {
+
+    static constexpr auto child_type = ChildType::Neutral;
+
+    /** @todo Need to implement the helper class that handles generic
+     *        children and return an instance of it. As a 1-member variant.
+     */
+    template <typename mask_v, typename float_v, typename proxy>
+    static auto prepareChild( [[maybe_unused]] LHCb::LinAlg::Vec<float_v, 3> const& vertex_position,
+                              [[maybe_unused]] mask_v const& vertex_position_mask, proxy const& child ) {
+      using VertexNeutral =
+          VertexTraj<proxy, child_type /* flag this represents a track not a weakly-decaying composite */>;
+      return VertexNeutral{child, Sel::stateVectorFromNeutral( child )};
+    }
+  };
+
+  template <class Proxy>
+  using child_helper_for_proxy_t = ChildHelper<child_type_for_proxy_v<Proxy>>;
+} // namespace
+
 namespace Sel::Fitters {
 
   /** @class Sel::Fitters::ParticleVertex
@@ -44,91 +511,6 @@ namespace Sel::Fitters {
     }
 
   private:
-    template <std::size_t N, typename T>
-    using VecN = LHCb::LinAlg::Vec<T, N>;
-
-    template <std::size_t N, typename T>
-    using SymNxN = LHCb::LinAlg::MatSym<T, N>;
-
-    template <std::size_t N, std::size_t M, typename T>
-    using Matrix = LHCb::LinAlg::Mat<T, N, M>;
-
-    template <std::size_t N, typename T>
-    using MatrixNxN = LHCb::LinAlg::Mat<T, N>;
-
-    /* TODO look at merging this with the stuff in State4.h
-     */
-    template <typename F>
-    struct MiniState {
-      SymNxN<5, F> cov{};
-      auto&        x() { return m_x; }
-      auto const&  x() const { return m_x; }
-      auto&        y() { return m_y; }
-      auto const&  y() const { return m_y; }
-      auto&        z() { return m_z; }
-      auto const&  z() const { return m_z; }
-      auto&        qOverP() { return m_qOverP; }
-      auto const&  qOverP() const { return m_qOverP; }
-      auto&        tx() { return m_tx; }
-      auto const&  tx() const { return m_tx; }
-      auto&        ty() { return m_ty; }
-      auto const&  ty() const { return m_ty; }
-
-      template <typename T>
-      MiniState( T const& other )
-          : m_x{other.x()}
-          , m_y{other.y()}
-          , m_tx{other.tx()}
-          , m_ty{other.ty()}
-          , m_z{other.z()}
-          , m_qOverP{other.qOverP()} {
-        static_assert( Sel::Utils::canBeExtrapolatedDownstream_v<T> );
-        auto const& other_cov = other.covariance();
-        LHCb::Utils::unwind<0, 5>( [this, &other_cov]( auto i ) {
-          LHCb::Utils::unwind<i, 5>( [this, i, &other_cov]( auto j ) { cov( i, j ) = other_cov( i, j ); } );
-        } );
-      }
-
-      void linearTransportTo( F const& new_z ) {
-        const auto dz  = new_z - m_z;
-        const auto dz2 = dz * dz;
-        m_x += dz * m_tx;
-        m_y += dz * m_ty;
-        cov( 0, 0 ) += dz2 * cov( 2, 2 ) + 2 * dz * cov( 2, 0 );
-        cov( 1, 0 ) += dz2 * cov( 3, 2 ) + dz * ( cov( 3, 0 ) + cov( 2, 1 ) );
-        cov( 2, 0 ) += dz * cov( 2, 2 );
-        cov( 3, 0 ) += dz * cov( 3, 2 );
-        cov( 4, 0 ) += dz * cov( 4, 2 );
-        cov( 1, 1 ) += dz2 * cov( 3, 3 ) + 2 * dz * cov( 3, 1 );
-        cov( 2, 1 ) += dz * cov( 3, 2 );
-        cov( 3, 1 ) += dz * cov( 3, 3 );
-        cov( 4, 1 ) += dz * cov( 4, 3 );
-        m_z = new_z;
-      }
-
-      SymNxN<2, F>    covXX() const { return cov.template sub<SymNxN<2, F>, 0, 0>(); }
-      MatrixNxN<2, F> covXT() const { return cov.template sub<MatrixNxN<2, F>, 2, 0>(); }
-
-    private:
-      F m_x{}, m_y{}, m_tx{}, m_ty{}, m_z{}, m_qOverP{};
-    };
-
-    // Different daughter types
-    // TODO look at using the meta enum thing
-    enum struct ChildType : unsigned char {
-      TrackWithVelo = 0,
-      TrackWithoutVelo,
-      Composite,
-      Resonance,
-      Neutral,
-      Other,
-      NumTypes
-    };
-
-    template <ChildType C>
-    static constexpr bool has_full_state_v = ( C == ChildType::Composite || C == ChildType::Resonance ||
-                                               C == ChildType::Neutral || C == ChildType::Other );
-
     /** Indirectly stolen from TrackVertexUtils::poca
      */
     template <typename float_v, typename mask_v>
@@ -169,518 +551,21 @@ namespace Sel::Fitters {
       return {std::string{message}, "Sel::Fitters::ParticleVertex", StatusCode::FAILURE};
     }
 
-    // Fast utility function for inverting a 2x2 matrix. This just comes
-    // from ROOT/Math/Dinv, but I took out the test on the determinant
-    // since we do not need it.
-    template <typename float_v>
-    static SymNxN<2, float_v> Invert2x2( SymNxN<2, float_v> const& rhs ) {
-      SymNxN<2, float_v> ret;
-      auto const         det    = rhs( 0, 0 ) * rhs( 1, 1 ) - rhs( 0, 1 ) * rhs( 0, 1 );
-      auto const         invdet = 1.f / det;
-      // if (det == T(0.)) { return false; }
-      ret( 0, 0 ) = rhs( 1, 1 ) * invdet;
-      ret( 1, 1 ) = rhs( 0, 0 ) * invdet;
-      ret( 0, 1 ) = -1.f * rhs( 0, 1 ) * invdet;
-      return ret;
-    }
-
-    // Class for daughters with a straight-line trajectory
-    template <typename particle_t, ChildType child_type>
-    struct VertexTraj final {
-      using simd_t  = typename particle_t::simd_t;
-      using int_v   = typename simd_t::int_v;
-      using float_v = typename simd_t::float_v;
-      using state_t = std::conditional_t<has_full_state_v<child_type>, State4<float_v>, MiniState<float_v>>;
-
-    private:
-      particle_t            m_particle; // TODO maybe just cache some information from this?
-      state_t               m_state;    // TODO evaluate what is really needed
-      VecN<2, float_v>      m_q;        // predicted/fitted slope (tx,ty)
-      SymNxN<2, float_v>    m_G;        // weight matrix of (x,y) of state
-      Matrix<2, 3, float_v> m_A;        // projection matrix for vertex position
-
-    public:
-      VertexTraj( particle_t p, state_t state )
-          : m_particle{std::move( p )}, m_state{std::move( state )}, m_q{m_state.tx(), m_state.ty()} {}
-
-      void project( VecN<3, float_v> const& vertexpos, VecN<3, float_v>& halfDChisqDX,
-                    SymNxN<3, float_v>& halfD2ChisqDX2, float_v& chi2, int_v& ndof ) {
-        // move the state (not sure we should do it this way: maybe better to just cache the z.)
-        m_state.linearTransportTo( vertexpos( 2 ) );
-
-        // compute the weight matrix
-        // TODO Just use invChol?
-        m_G = Invert2x2( m_state.covXX() );
-
-        // compute residual
-        VecN<2, float_v> residual{vertexpos( 0 ) - m_state.x(), vertexpos( 1 ) - m_state.y()};
-
-        // fill the projection matrix: use the fitted momentum!
-        m_A( 0, 0 ) = m_A( 1, 1 ) = 1.f;
-        m_A( 0, 2 )               = -1.f * m_q( 0 );
-        m_A( 1, 2 )               = -1.f * m_q( 1 );
-
-        // I tried to make this faster by writing it out, but H does
-        // not contain sufficiently many zeroes. Better to
-        // parallelise.
-
-        // I am not sure that the compiler realizes that A^T*G is used
-        // more than once here, so I'll substitute the code from
-        // similarity. Note that similarity does not return an expression.
-        // halfD2ChisqDX2 += ROOT::Math::Similarity( ROOT::Math::Transpose(m_A), m_G ) ;
-        // halfDChisqDX   += ( ROOT::Math::Transpose(m_A) * m_G) * res ;
-        auto const ATG = m_A.transpose() * m_G.cast_to_mat();
-        // updating relevant fit results
-        halfD2ChisqDX2 = halfD2ChisqDX2 + ( ATG * m_A ).cast_to_sym();
-        halfDChisqDX   = halfDChisqDX + ATG * residual;
-        chi2           = chi2 + LHCb::LinAlg::similarity( residual, m_G );
-        ndof           = ndof + 2;
-      }
-
-      void updateSlopes( VecN<3, float_v> const& vertexpos ) {
-        // first update the residual. (note the subtle difference with that in project!)
-        float_v const          dz = vertexpos( 2 ) - m_state.z();
-        VecN<2, float_v> const residual{vertexpos( 0 ) - ( m_state.x() + m_state.tx() * dz ),
-                                        vertexpos( 1 ) - ( m_state.y() + m_state.ty() * dz )};
-
-        // get the matrix that is the correlation of (x,y) and (tx,ty,qop)
-        // ROOT::Math::SMatrix<double,2,2> Vba = m_state.covariance().template
-        // Sub<ROOT::Math::SMatrix<double,2,2>>(2,0); compute the corresponding gain matrix (this is WBG in the BFR fit,
-        // but here it is Vba * Vaa^-1) ROOT::Math::SMatrix<double,2,2> K = Vba*m_G ; compute the momentum vector
-        m_q =
-            VecN<2, float_v>{{m_state.tx(), m_state.ty()}} + m_state.covXT().transpose() * m_G.cast_to_mat() * residual;
-      }
-
-      void addToFourVector( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>& p4cov,
-                            Matrix<4, 3, float_v>& gainmatrix ) const {
-        // Despatch is a bit crude, but it'll do...in the original
-        // ParticleVertexFitter this was done with a full specialisation, but
-        // now it would have to be a partial specialisation and it all gets a
-        // bit messy.
-        if constexpr ( child_type == ChildType::Composite ) {
-          addToFourVector_Composite( vertexpos, p4, p4cov, gainmatrix );
-        } else if constexpr ( child_type == ChildType::TrackWithVelo || child_type == ChildType::TrackWithoutVelo ) {
-          addToFourVector_Track( vertexpos, p4, p4cov, gainmatrix );
-        } else if constexpr ( child_type == ChildType::Neutral ) {
-          addToFourVector_Neutral( vertexpos, p4, p4cov, gainmatrix );
-        } else {
-        }
-      }
-
-      /** Add the momenta of a neutral particle to the result of a combination
-
-          Update the information of the momenta of the neutral objects based on the vertex of
-          the decaying particle. It is crucial that we add the neutral objects after
-          determining the decay vertex (i.e. at the end of the combination).
-
-          @todo FIXME correctly process the covariance matrix
-       */
-      void addToFourVector_Neutral( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>&,
-                                    Matrix<4, 3, float_v>& ) const {
-        // TODO: make this a global function that modifies the neutral objects at the
-        // persistency stage
-        auto dx = m_state.x() - vertexpos( 0 );
-        auto dy = m_state.y() - vertexpos( 1 );
-        auto dz = m_state.z() - vertexpos( 2 );
-
-        auto d = sqrt( dx * dx + dy * dy + dz * dz );
-
-        auto p = m_particle.p();
-
-        p4 = VecN<4, float_v>{p * dx / d, p * dy / d, p * dz / d, m_particle.energy()};
-      }
-
-      /** @todo FIXME implement logic for bremsstrahlung photons.
-       */
-      void addToFourVector_Track( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>& p4cov,
-                                  Matrix<4, 3, float_v>& gain_matrix ) const {
-        // This specialisation is for TrackWithVelo and TrackWithoutVelo
-        // we first need to update q/p. since we also need the full gain matrix, we could as well redo that part.
-        float_v const          dz = vertexpos( 2 ) - m_state.z();
-        VecN<2, float_v> const residual{vertexpos( 0 ) - ( m_state.x() + m_q( 0 ) * dz ),
-                                        vertexpos( 1 ) - ( m_state.y() + m_q( 1 ) * dz )};
-
-        auto const Vba    = m_state.cov.template sub<Matrix<3, 2, float_v>, 2, 0>();
-        auto       K      = Vba * m_G.cast_to_mat();
-        auto       tx     = m_state.tx() + ( K * residual )( 0 );
-        auto       ty     = m_state.ty() + ( K * residual )( 1 );
-        auto       qOverP = m_state.qOverP() + ( K * residual )( 2 );
-        // update 4-momentum
-        // TODO can we just store this value in the constructor? what else is needed? bremsstrahlung info?
-        auto const mass = m_particle.mass();
-        auto const p    = 1.f / abs( qOverP );
-        auto const n2   = 1.f + tx * tx + ty * ty;
-        auto const n    = sqrt( n2 );
-        auto const pz   = p / n;
-        auto const px   = tx * pz;
-        auto const py   = ty * pz;
-        auto const E    = sqrt( p * p + mass * mass );
-        p4              = p4 + VecN<4, float_v>{{px, py, pz, E}}; // TODO check whether this gives right results
-                                                                  // it was E px py pz before
-
-        Matrix<4, 3, float_v> dp4dmom{}; // jacobi matrix of del(p4)/del(tx,ty,qOverP)
-        auto const            n3    = n2 * n;
-        auto const            invn3 = 1.f / n3;
-        dp4dmom( 0, 0 )             = p * ( 1.f + ty * ty ) * invn3; // dpx/dtx
-        dp4dmom( 0, 1 )             = p * tx * -ty * invn3;          // dpx/dty
-        dp4dmom( 0, 2 )             = -px / qOverP;                  // dpx/dqop
-
-        dp4dmom( 1, 0 ) = p * ty * -tx * invn3;        // dpy/dtx
-        dp4dmom( 1, 1 ) = p * ( 1 + tx * tx ) * invn3; // dpy/dty
-        dp4dmom( 1, 2 ) = -py / qOverP;                // dpy/dqop
-
-        dp4dmom( 2, 0 ) = pz * -tx / n2; // dpz/dtx
-        dp4dmom( 2, 1 ) = pz * -ty / n2; // dpz/dtx
-        dp4dmom( 2, 2 ) = -pz / qOverP;  // dpz/dqop
-
-        dp4dmom( 3, 0 ) = 0.0;                 // dE/dtx
-        dp4dmom( 3, 1 ) = 0.0;                 // dE/dty
-        dp4dmom( 3, 2 ) = p / E * -p / qOverP; // dE/dqop
-
-        auto FK = dp4dmom * K;
-
-        p4cov = p4cov + LHCb::LinAlg::similarity( dp4dmom, m_state.cov.template sub<SymNxN<3, float_v>, 2, 2>() ) -
-                LHCb::LinAlg::similarity( FK, m_state.cov.template sub<SymNxN<2, float_v>, 0, 0>() );
-
-        gain_matrix = gain_matrix + FK * m_A;
-        // // For brem-recovered electrons, we need to do something
-        // // special. So, electrons are an ugly exception here. We could
-        // // also add to q/p instead, which is cheaper, but perhaps even
-        // // more ugly.
-        // if ( particle().particleID().abspid() == 11 ) {
-        //   const double absqop               = std::abs( m_state.qOverP() );
-        //   const double momentumFromParticle = particle().momentum().P();
-        //   // is 1% a reasonable threshold?
-        //   if ( momentumFromParticle * absqop > 1.01 ) {
-        //     const double momentumFromTrack       = 1 / absqop;
-        //     const double momentumError2FromTrack = m_state.covariance()( 4, 4 ) * std::pow(
-        //     momentumFromTrack, 4 ); const double momentumError2FromParticle =
-        //         particle().momCovMatrix()( 3, 3 ) * std::pow( particle().momentum().E() /
-        //         momentumFromParticle, 2 );
-        //     const double bremEnergyCov = momentumError2FromParticle - momentumError2FromTrack;
-        //     const double bremEnergy    = momentumFromParticle - momentumFromTrack;
-        //     // if the correction is unphysical, ignore it.
-        //     if ( bremEnergyCov > 0 ) {
-        //       const auto tx = q( 0 );
-        //       const auto ty = q( 1 );
-        //       auto       t  = std::sqrt( 1 + tx * tx + ty * ty );
-        //       // we could also 'scale' the momentum, but we anyway need the components for the
-        //       jacobian Gaudi::LorentzVector p4brem{bremEnergy * tx / t, bremEnergy * ty / t,
-        //       bremEnergy / t, bremEnergy}; p4 += p4brem; ROOT::Math::SMatrix<double, 4, 1> J;
-        //       J( 0, 0 ) = tx / t;
-        //       J( 1, 0 ) = ty / t;
-        //       J( 2, 0 ) = 1 / t;
-        //       J( 3, 0 ) = 1;
-        //       p4cov += ROOT::Math::Similarity( J, Gaudi::SymMatrix1x1{bremEnergyCov} );
-        //     }
-        //   }
-        // }
-      }
-
-      void addToFourVector_Composite( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4,
-                                      SymNxN<4, float_v>& p4cov, Matrix<4, 3, float_v>& gainmatrix ) const {
-        // first need to 'update' the momentum vector. but for that we
-        // first need to 'transport'.
-        using Sel::Utils::fourMomentum;
-        auto          mom = fourMomentum( m_particle ); // 4-momentum
-        float_v const pz  = Z( mom );
-        float_v const tx  = m_q( 0 );
-        float_v const ty  = m_q( 1 );
-
-        // first update the residual
-        using namespace LHCb::LinAlg;
-        using Sel::Utils::referencePoint;
-        auto             particlepos = referencePoint( m_particle );
-        float_v const    dz          = Z( vertexpos ) - particlepos.Z();
-        VecN<2, float_v> res{X( vertexpos ) - ( particlepos.X() + tx * dz ),
-                             Y( vertexpos ) - ( particlepos.Y() + ty * dz )};
-        auto const&      R    = m_state.covXX();
-        auto const&      Rinv = m_G;
-        using Sel::Utils::momCovMatrix;
-        auto const momCov = momCovMatrix( m_particle );
-        using Sel::Utils::momPosCovMatrix;
-        auto const momPosCov = momPosCovMatrix( m_particle );
-
-        // To do this right we need THREE projection matrices for the residual:
-        //    r = A*x_vertex + B*mom_dau + C*x_dau
-        // Note that C=-A.
-        // Lets call the matrix F^T = (B^T C^T). then the uncertainty on the residual is
-        //    R = F V77 F^T where V77 is the 7x7 covariance matrix of the daughter
-        // We will need R every iteration
-        // The correlation matrix between residual and the 7 parameters of the daughter is
-        //    Vr_xmom = V77 F^T
-        // Consequently, the gain matrix for the momentum-decayvertex is
-        //    K72 = V77 F^T R^-1
-        // Now, we do not need the updated parameters for the decay
-        // vertex. If we just need the momentum part, then this reduces to
-        //    K42 = ( V43 * C^T + V44 * B^T ) R^-1
-
-        // The easiest implementation is this ...
-        //  Gaudi::Matrix3x2 CT{ ROOT::Math::Transpose(m_A) } ;
-        //  ROOT::Math::SMatrix<double,4,2> BT ;
-        //  BT(0,0) = +dz/pz;
-        //  BT(2,0) = -dz/pz*tx ;
-        //  BT(1,1) = +dz/pz ;
-        //  BT(2,1) = -dz/pz*ty ;
-        //  ROOT::Math::SMatrix<double,4,2> K42 = ( momPosCov * CT + momCov * BT ) * Rinv ;
-        // but since B and C are largely empty, we better write it out:
-        float_v const         dzopz = dz / pz;
-        Matrix<4, 2, float_v> K42part; // uninitialised here
-        for ( int i = 0; i < 4; ++i ) {
-          K42part( i, 0 ) =
-              momPosCov( i, 0 ) - momPosCov( i, 2 ) * tx + momCov( i, 0 ) * dzopz - momCov( i, 2 ) * ( dzopz * tx );
-          K42part( i, 1 ) =
-              momPosCov( i, 1 ) - momPosCov( i, 2 ) * ty + momCov( i, 1 ) * dzopz - momCov( i, 2 ) * ( dzopz * ty );
-        }
-        Matrix<4, 2, float_v> const K42 = K42part * Rinv.cast_to_mat();
-
-        // update the fourvector
-        auto const deltap4 = K42 * res;
-        p4                 = p4 + mom + deltap4;
-        // to understand this, compare to what is done for the track
-        p4cov      = p4cov + momCov - similarity( K42, R );
-        gainmatrix = gainmatrix + K42 * m_A;
-      }
-    };
-
-    /** Helper type for identifying and handling TrackWithVelo and
-     *  TrackWithoutVelo children.
-     *  @todo Handle TrackWithoutVelo.
-     */
-    template <typename child_t>
-    struct TrackWithVeloChildHelper {
-      /** Determine if the given `child_t` type represents a [chunk of] tracks.
-       *  @todo Make this more sophisticated and follow changes in the event
-       *        model.
-       */
-      constexpr static bool value = Sel::Utils::canBeExtrapolatedDownstream_v<child_t>;
-
-      /** Given that `value` was `true`, derive (at runtime) the `ChildType` of
-       *  `child_t`. This basically means "is it a long or downstream track".
-       *  @todo Implement check and maybe return `ChildType::TrackWithoutVelo`.
-       *  @todo Clearly state what is supported when `child` is a vector type.
-       *        For example, if it is a proxy into a container that supports
-       *        heterogeneous row types (e.g. some rows are long, some rows are
-       *        downstream) then there is not a single correct return value we
-       *        can give here. Is that a runtime error? Should the return type
-       *        be changed so that can be supported?
-       */
-      static ChildType deriveType( LHCb::IParticlePropertySvc const&, [[maybe_unused]] child_t const& child ) {
-        return ChildType::TrackWithVelo;
-      }
-
-      /** Given that `value` was `true`, and that `type` was returned by our
-       *  `deriveType` function, return the narrowest-possible std::variant of
-       *  types implementing the API used in the vertex fit. For tracks there
-       *  should just be one return type (i.e. a variant of size 1), but there
-       *  will be some runtime branching on the `type` to figure out how to get
-       *  the relevant state.
-       *  @todo Use `vertex_position` and its validity mask
-       *        `vertex_position_mask` to determine if a sophisticated
-       *        extrapolation is needed.
-       *  @todo Implement this logic, don't just blindly call the
-       *        `state(StateLocation::ClosestToBeam)` accessor. Remove [[maybe_unused]].
-       */
-      template <typename mask_v, typename float_v>
-      static auto prepareChild( [[maybe_unused]] VecN<3, float_v> const& vertex_position,
-                                [[maybe_unused]] mask_v const& vertex_position_mask, [[maybe_unused]] ChildType type,
-                                child_t const& child ) {
-        assert( type == ChildType::TrackWithVelo );
-        using VertexTrack =
-            VertexTraj<child_t,
-                       ChildType::TrackWithVelo /* flag this represents a track not a weakly-decaying composite */>;
-        return VertexTrack{child, child.state( LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam )};
-      }
-    };
-
-    /** Helper type for identifying and handling TrackWithVelo and
-     *  TrackWithoutVelo children.
-     *  @todo Handle TrackWithoutVelo.
-     */
-    template <typename child_t>
-    struct TrackWithoutVeloChildHelper {
-      /** Determine if the given `child_t` type represents a [chunk of] tracks.
-       *  @todo Make this more sophisticated and follow changes in the event
-       *        model.
-       */
-      constexpr static bool value = Sel::Utils::has_tracklike_API<child_t>;
-
-      /** Given that `value` was `true`, derive (at runtime) the `ChildType` of
-       *  `child_t`. This basically means "is it a long or downstream track".
-       *  @todo Implement check and maybe return `ChildType::TrackWithoutVelo`.
-       *  @todo Clearly state what is supported when `child` is a vector type.
-       *        For example, if it is a proxy into a container that supports
-       *        heterogeneous row types (e.g. some rows are long, some rows are
-       *        downstream) then there is not a single correct return value we
-       *        can give here. Is that a runtime error? Should the return type
-       *        be changed so that can be supported?
-       */
-      static ChildType deriveType( LHCb::IParticlePropertySvc const&, [[maybe_unused]] child_t const& child ) {
-        return ChildType::TrackWithoutVelo;
-      }
-
-      /** Given that `value` was `true`, and that `type` was returned by our
-       *  `deriveType` function, return the narrowest-possible std::variant of
-       *  types implementing the API used in the vertex fit. For tracks there
-       *  should just be one return type (i.e. a variant of size 1), but there
-       *  will be some runtime branching on the `type` to figure out how to get
-       *  the relevant state.
-       *  @todo Use `vertex_position` and its validity mask
-       *        `vertex_position_mask` to determine if a sophisticated
-       *        extrapolation is needed.
-       *  @todo Implement this logic, don't just blindly call the
-       *        `state(StateLocation::ClosestToBeam)` accessor. Remove [[maybe_unused]].
-       */
-      template <typename mask_v, typename float_v>
-      static auto prepareChild( [[maybe_unused]] VecN<3, float_v> const& vertex_position,
-                                [[maybe_unused]] mask_v const& vertex_position_mask, [[maybe_unused]] ChildType type,
-                                child_t const& child ) {
-        assert( type == ChildType::TrackWithoutVelo );
-        using VertexTrack =
-            VertexTraj<child_t,
-                       ChildType::TrackWithoutVelo /* flag this represents a track not a weakly-decaying composite */>;
-        return VertexTrack{child, child.state( LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam )};
-      }
-    };
-
-    /** Helper type for identifying and handling Composite and Resonance
-     *  children.
-     *  @todo Handle Resonance, adding a runtime check on PID and another type
-     *        to go in the variant returned by prepareChild.
-     */
-    template <typename child_t>
-    struct CompositeChildHelper {
-      /** Determine if the given `child_t` represents a [chunk of] composite
-       *  particles.
-       */
-      constexpr static bool value = !Sel::Utils::canBeExtrapolatedDownstream_v<child_t>;
-
-      /** Given that `value` was `true`, derive (at runtime) the `ChildType` of
-       *  `child_t`. This basically means "is it a weakly decaying particle".
-       *  @todo Implement check and maybe return `ChildType::Resonance`.
-       *  @todo Can use `pp_svc` to get the "ctau" of the given `child`. Remove
-       *        [[maybe_unused]] once we do that.
-       */
-      static ChildType deriveType( [[maybe_unused]] LHCb::IParticlePropertySvc const& pp_svc,
-                                   [[maybe_unused]] child_t const&                    child ) {
-        return ChildType::Composite;
-      }
-      /** Given that `value` was `true`, and that `type` was returned by our
-       *  `deriveType` function, return the narrowest-possible std::variant of
-       *  types implementing the API used in the vertex fit. Because weakly
-       *  decaying particles are handled differently in the fit, the return
-       *  type here will have to be a variant of two different types.
-       *
-       *  The original ParticleVertexFitter didn't use the estimated vertex
-       *  position for `Resonance` and `Composite` children, so we don't
-       *  either.
-       *  @todo Add logic for the `ChildType::Resonance` value of `type` and
-       *        add another member to the return type; remove [[maybe_unused]]
-       */
-      template <typename mask_v, typename float_v>
-      static auto prepareChild( VecN<3, float_v> const&, mask_v const&, [[maybe_unused]] ChildType type,
-                                child_t const& child ) {
-        assert( type == ChildType::Composite || type == ChildType::Resonance );
-        using VertexComposite =
-            VertexTraj<child_t,
-                       ChildType::Composite /* flag this represents a weakly-decaying composite and not a track */>;
-        return VertexComposite{child, stateVectorFromComposite( child )};
-      }
-    };
-
-    /** Helper type for the default `Neutral` case.
-     */
-    template <typename child_t>
-    struct NeutralChildHelper {
-      /** This helper just tries generic logic, it should be tried last and it
-       *  accepts everything. If it doesn't work, a compile error is the
-       *  correct outcome.
-       */
-      constexpr static bool value = Sel::Utils::has_neutrals_API<child_t>;
-      /** Catch-all, don't expect other logic here.
-       */
-      static ChildType deriveType( LHCb::IParticlePropertySvc const&, child_t const& ) { return ChildType::Neutral; }
-      /** @todo Need to implement the helper class that handles generic
-       *        children and return an instance of it. As a 1-member variant.
-       */
-      template <typename mask_v, typename float_v>
-      static auto prepareChild( [[maybe_unused]] VecN<3, float_v> const& vertex_position,
-                                [[maybe_unused]] mask_v const& vertex_position_mask, [[maybe_unused]] ChildType type,
-                                child_t const& child ) {
-        assert( type == ChildType::Neutral );
-        using VertexNeutral =
-            VertexTraj<child_t, ChildType::Neutral /* flag this represents a track not a weakly-decaying composite */>;
-        return VertexNeutral{child, stateVectorFromNeutral( child )};
-      }
-    };
-
-    /** Helper type for the default `Neutral` case.
-     */
-    template <typename child_t>
-    struct LeftoverHelper {
-      /** This helper just tries generic logic, it should be tried last and it
-       *  accepts everything. If it doesn't work, a compile error is the
-       *  correct outcome.
-       */
-      constexpr static bool value = true;
-      /** Catch-all, don't expect other logic here.
-       */
-      static ChildType deriveType( LHCb::IParticlePropertySvc const&, child_t const& ) { return ChildType::Neutral; }
-      /** @todo Need to implement the helper class that handles generic
-       *        children and return an instance of it. As a 1-member variant.
-       */
-      template <typename mask_v, typename float_v>
-      static auto prepareChild( [[maybe_unused]] VecN<3, float_v> const& vertex_position,
-                                [[maybe_unused]] mask_v const& vertex_position_mask, [[maybe_unused]] ChildType type,
-                                child_t const& child ) {
-        assert( type == ChildType::Other );
-        return bool{};
-      }
-    };
-
-    /** Master list of helper types. The order matters; the **first** matching
-     *  helper is preferred.
-     */
-    template <typename child_t>
-    using all_helper_types =
-        boost::mp11::mp_list<TrackWithVeloChildHelper<child_t>, TrackWithoutVeloChildHelper<child_t>,
-                             NeutralChildHelper<child_t>, CompositeChildHelper<child_t>, LeftoverHelper<child_t>>;
-
-    /** Helper types that match (::value == true) for a given child type.
-     */
-    template <typename child_t>
-    using matching_helper_types = boost::mp11::mp_copy_if<all_helper_types<child_t>, boost::mp11::mp_to_bool>;
-
-    /** The first matching helper type for a given child type.
-     */
-    template <typename child_t>
-    using first_matching_helper_type = boost::mp11::mp_front<matching_helper_types<child_t>>;
-
-    /** Take a [SIMD vector worth of] particles and assign a single "type"
-     *  value that will steer the fit algorithm below. It is an error
-     *  (unfortunately only a runtime one...) to give children that are of
-     *  heterogeneous `ChildType`, but it should be supported that there is
-     *  a mix of actual PDG IDs (e.g. `child` could be a mix of D0 and D~0, but
-     *  there should be a single return type (ChildType::Composite, in that
-     *  case).
-     */
-    template <typename child_t>
-    ChildType deriveType( child_t const& child ) const {
-      assert( m_pp_svc );
-      using helper_t = first_matching_helper_type<child_t>;
-      return helper_t::deriveType( *m_pp_svc, child );
-    }
-
     /** Helper for getting the return type of the prepareChild method of a
      *  helper type.
      */
-    template <typename mask_v, typename float_v>
-    struct prepareChild_helper {
-      template <typename helper_t, typename variant_member_t>
-      using fn = decltype( helper_t::prepareChild( std::declval<VecN<3, float_v>>(), std::declval<mask_v>(),
-                                                   std::declval<ChildType>(), std::declval<variant_member_t>() ) );
+    template <typename mask_v, typename float_v, typename VariantTypes>
+    struct return_types_helper;
+
+    template <typename mask_v, typename float_v, template <class...> class Variant, typename... Proxy>
+    struct return_types_helper<mask_v, float_v, Variant<Proxy...>> {
+      using type = std::variant<decltype( child_helper_for_proxy_t<Proxy>::prepareChild(
+          std::declval<VecN<3, float_v>>(), std::declval<mask_v>(), std::declval<Proxy>() ) )...>;
     };
 
+    template <typename mask_v, typename float_v, typename VariantTypes>
+    using return_types_helper_t = typename return_types_helper<mask_v, float_v, VariantTypes>::type;
+
     /** Choose the first match for the given child type in `all_helper_types`
      *  and call its `prepareChild` method. All the extra complexity is to cope
      *  with `child_t` being a variant type whose members might match different
@@ -689,28 +574,21 @@ namespace Sel::Fitters {
      *  @todo Include the `vertexinsidefoil` logic from ParticleVertexFitter.cpp
      */
     template <typename mask_v, typename float_v, typename child_t>
-    auto prepareChild( VecN<3, float_v> const& vertex_position, mask_v const& vertex_position_mask, ChildType type,
+    auto prepareChild( VecN<3, float_v> const& vertex_position, mask_v const& vertex_position_mask,
                        child_t const& possibly_variant_child ) const {
       // child_t might be variant<Ts...>, or it might be a value type itself.
       // if the former, this yields tuple<Ts...>, otherwise tuple<child_t>
       using variant_types = typename LHCb::is_variant<child_t>::tuple_type;
-      // get the [first] matching helper for each element in `variant_types`
-      using matched_helpers = boost::mp11::mp_transform<first_matching_helper_type, variant_types>;
-      // get the return types of the `prepareChild` methods of all of these
-      // helpers when they are invoked with the corresponding element of
-      // `variant_types`. e.g. this might be tuple<variant<A, B>, variant<C>>
-      using return_types =
-          boost::mp11::mp_transform_q<prepareChild_helper<mask_v, float_v>, matched_helpers, variant_types>;
-      // Turn this into a variant<...> type that can accommodate all relevant
-      // types, in the example just above `ret_t` would be variant<A, B, C>
-      using ret_t = boost::mp11::mp_rename<return_types, std::variant>;
+
+      using ret_t = return_types_helper_t<mask_v, float_v, variant_types>;
+
       static_assert( boost::mp11::mp_is_set<ret_t>::value,
                      "Got duplicate entries in derived variant type. Looks like the code needs an mp_unique<...>." );
       // Now we know what return type to impose, invoke/visit on `child`
       return LHCb::invoke_or_visit(
           [&]( auto const& child ) -> ret_t {
-            using helper_t = first_matching_helper_type<std::decay_t<decltype( child )>>;
-            return helper_t::prepareChild( vertex_position, vertex_position_mask, type, child );
+            using helper_t = child_helper_for_proxy_t<std::decay_t<decltype( child )>>;
+            return helper_t::prepareChild( vertex_position, vertex_position_mask, child );
           },
           possibly_variant_child );
     }
@@ -748,7 +626,7 @@ namespace Sel::Fitters {
             LHCb::invoke_or_visit(
                 [&]( auto const& child ) {
                   using Sel::Utils::endVertexPos;
-                  if constexpr ( !Sel::Utils::canBeExtrapolatedDownstream_v<decltype( child )> ) {
+                  if constexpr ( !Sel::Utils::isBasicParticle_v<decltype( child )> ) {
                     tmp.emplace( LHCb::LinAlg::convert( endVertexPos( child ) ), mask_v{true} );
                   } else {
                     throw exception( "Resonance type didn't have an endVertex()" );
@@ -770,18 +648,17 @@ namespace Sel::Fitters {
             // a variant type.
             LHCb::invoke_or_visit(
                 [&]( auto const& child ) {
+                  using child_type = decltype( child );
+
                   using Sel::Utils::trackState;
                   using Sel::Utils::referencePoint;
                   using Sel::Utils::threeMomentum;
-                  if ( types[i] == ChildType::TrackWithVelo ) {
-                    if constexpr ( Sel::Utils::canBeExtrapolatedDownstream_v<decltype( child )> ) {
+
+                  if constexpr ( Sel::Utils::isBasicParticle_v<child_type> ) {
+                    if constexpr ( Sel::Utils::canBeExtrapolatedDownstream_v<child_type> )
                       velotrajs[nvelotrajs++] = {trackState( child )};
-                    } else {
-                      throw exception( "TrackWithVelo type didn't have a state( "
-                                       "LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam)" );
-                    }
-                  } else if ( types[i] == ChildType::Composite ) {
-                    if constexpr ( !Sel::Utils::canBeExtrapolatedDownstream_v<decltype( child )> ) {
+                  } else {
+                    if constexpr ( !Sel::Utils::canBeExtrapolatedDownstream_v<child_type> ) {
                       velotrajs[nvelotrajs++] = {referencePoint( child ), threeMomentum( child )};
                     } else {
                       throw exception( "Composite type didn't one/both of referencePoint() and threeMomentum()" );
@@ -852,8 +729,9 @@ namespace Sel::Fitters {
       // the template pack `child_ts...`. If the input children are themselves
       // variants then `deriveType` is used as a visitor, increasing the amount
       // of runtime dispatch.
-      std::array const types{LHCb::invoke_or_visit( [this]( auto const& child ) { return this->deriveType( child ); },
-                                                    children.template get<IChild>() )...};
+      std::array const types{LHCb::invoke_or_visit(
+          [this]( auto const& child ) { return child_type_for_proxy_v<std::decay_t<decltype( child )>>; },
+          children.template get<IChild>() )...};
 
       // Obtain an estimate of the vertex position for initialization
       auto [vertex_position, vertex_position_mask] = calculateInitialPosition<simd_t>( types, children );
@@ -872,7 +750,7 @@ namespace Sel::Fitters {
       // filled with variants; the different variant types encode the different
       // handling required for TrackWith[out]Velo/Composite/Resonance/Neutral/Other.
       std::tuple extrapolated_children{
-          prepareChild( vertex_position, prelim_fit_mask, types[IChild], children.template get<IChild>() )...};
+          prepareChild( vertex_position, prelim_fit_mask, children.template get<IChild>() )...};
 
       // These are required for populating child relations at the end, but the
       // API for dealing with them should be agnostic to the child types
@@ -883,7 +761,7 @@ namespace Sel::Fitters {
 
       auto num_descendants = ( LHCb::invoke_or_visit(
                                    []( [[maybe_unused]] auto const& child ) -> std::size_t {
-                                     if constexpr ( !Sel::Utils::canBeExtrapolatedDownstream_v<decltype( child )> )
+                                     if constexpr ( !Sel::Utils::isBasicParticle_v<decltype( child )> )
                                        return child.numDescendants();
                                      else
                                        return 1;
@@ -896,7 +774,7 @@ namespace Sel::Fitters {
 
       ( LHCb::invoke_or_visit(
             [&current_unique_ids]( auto const& child ) -> void {
-              if constexpr ( !Sel::Utils::canBeExtrapolatedDownstream_v<decltype( child )> ) {
+              if constexpr ( !Sel::Utils::isBasicParticle_v<decltype( child )> ) {
                 for ( auto i = 0u; i < child.numDescendants(); ++i )
                   current_unique_ids.emplace_back( child.descendant_unique_id( i ) );
               } else {
-- 
GitLab


From 8e01508e4c0087477d9fac2efc889f2df53a952e Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Thu, 31 Mar 2022 16:35:25 +0100
Subject: [PATCH 23/26] Redefine states to allow neutrals

---
 Phys/ParticleMaker/src/NeutralMakers.cpp      |   1 -
 Phys/SelTools/include/SelTools/Utilities.h    | 109 ++++++-
 .../include/VertexFit/ParticleVertexFitter.h  | 274 ++++++++----------
 3 files changed, 221 insertions(+), 163 deletions(-)

diff --git a/Phys/ParticleMaker/src/NeutralMakers.cpp b/Phys/ParticleMaker/src/NeutralMakers.cpp
index 62243130782..1fadde1eb4f 100644
--- a/Phys/ParticleMaker/src/NeutralMakers.cpp
+++ b/Phys/ParticleMaker/src/NeutralMakers.cpp
@@ -641,7 +641,6 @@ namespace LHCb::Phys::ParticleMakers {
    *
    * @brief Create a container of neutral objects from the output of the reconstruction sequence and a set of vertices
    *
-   *
    */
   class NeutralBasicsMaker : public Gaudi::Functional::Transformer<LHCb::Event::NeutralBasics(
                                  LHCb::UniqueIDGenerator const&, LHCb::Event::Calo::v2::Hypotheses const&,
diff --git a/Phys/SelTools/include/SelTools/Utilities.h b/Phys/SelTools/include/SelTools/Utilities.h
index 0dd1ebe108a..f2f9f847707 100644
--- a/Phys/SelTools/include/SelTools/Utilities.h
+++ b/Phys/SelTools/include/SelTools/Utilities.h
@@ -15,6 +15,97 @@
 
 namespace Sel {
 
+  /** Common state type for non-composite objects with covariance matrix information
+
+     TODO look at merging this with the stuff in State4.h
+   */
+  template <typename FloatType>
+  struct MiniState {
+    LHCb::LinAlg::MatSym<FloatType, 5> cov{};
+    auto&                              x() { return m_x; }
+    auto const&                        x() const { return m_x; }
+    auto&                              y() { return m_y; }
+    auto const&                        y() const { return m_y; }
+    auto&                              z() { return m_z; }
+    auto const&                        z() const { return m_z; }
+    auto&                              qOverP() { return m_qOverP; }
+    auto const&                        qOverP() const { return m_qOverP; }
+    auto&                              tx() { return m_tx; }
+    auto const&                        tx() const { return m_tx; }
+    auto&                              ty() { return m_ty; }
+    auto const&                        ty() const { return m_ty; }
+
+    template <typename T>
+    MiniState( T const& other )
+        : m_x{other.x()}, m_y{other.y()}, m_tx{other.tx()}, m_ty{other.ty()}, m_z{other.z()}, m_qOverP{other.qOverP()} {
+      static_assert( Sel::Utils::canBeExtrapolatedDownstream_v<T> );
+      auto const& other_cov = other.covariance();
+      LHCb::Utils::unwind<0, 5>( [this, &other_cov]( auto i ) {
+        LHCb::Utils::unwind<i, 5>( [this, i, &other_cov]( auto j ) { cov( i, j ) = other_cov( i, j ); } );
+      } );
+    }
+
+    void linearTransportTo( FloatType const& new_z ) {
+      const auto dz  = new_z - m_z;
+      const auto dz2 = dz * dz;
+      m_x += dz * m_tx;
+      m_y += dz * m_ty;
+      cov( 0, 0 ) += dz2 * cov( 2, 2 ) + 2 * dz * cov( 2, 0 );
+      cov( 1, 0 ) += dz2 * cov( 3, 2 ) + dz * ( cov( 3, 0 ) + cov( 2, 1 ) );
+      cov( 2, 0 ) += dz * cov( 2, 2 );
+      cov( 3, 0 ) += dz * cov( 3, 2 );
+      cov( 4, 0 ) += dz * cov( 4, 2 );
+      cov( 1, 1 ) += dz2 * cov( 3, 3 ) + 2 * dz * cov( 3, 1 );
+      cov( 2, 1 ) += dz * cov( 3, 2 );
+      cov( 3, 1 ) += dz * cov( 3, 3 );
+      cov( 4, 1 ) += dz * cov( 4, 3 );
+      m_z = new_z;
+    }
+
+    LHCb::LinAlg::MatSym<FloatType, 2> covXX() const {
+      return cov.template sub<LHCb::LinAlg::MatSym<FloatType, 2>, 0, 0>();
+    }
+    LHCb::LinAlg::Mat<FloatType, 2> covXT() const { return cov.template sub<LHCb::LinAlg::Mat<FloatType, 2>, 2, 0>(); }
+
+  private:
+    FloatType m_x{}, m_y{}, m_tx{}, m_ty{}, m_z{}, m_qOverP{};
+  };
+
+  /** State without covariance matrix information
+   */
+  template <typename FloatType>
+  class UnfittedState {
+
+  public:
+    UnfittedState( FloatType x, FloatType y, FloatType z, FloatType tx, FloatType ty, FloatType qOverP )
+        : m_x{std::move( x )}
+        , m_y{std::move( y )}
+        , m_z{std::move( z )}
+        , m_tx{std::move( tx )}
+        , m_ty{std::move( ty )}
+        , m_qOverP{std::move( qOverP )} {}
+
+    template <class StateType>
+    UnfittedState( StateType const& other )
+        : UnfittedState( other.x(), other.y(), other.z(), other.tx(), other.ty(), other.qOverP() ) {}
+
+    auto&       x() { return m_x; }
+    auto const& x() const { return m_x; }
+    auto&       y() { return m_y; }
+    auto const& y() const { return m_y; }
+    auto&       z() { return m_z; }
+    auto const& z() const { return m_z; }
+    auto&       qOverP() { return m_qOverP; }
+    auto const& qOverP() const { return m_qOverP; }
+    auto&       tx() { return m_tx; }
+    auto const& tx() const { return m_tx; }
+    auto&       ty() { return m_ty; }
+    auto const& ty() const { return m_ty; }
+
+  private:
+    FloatType m_x, m_y, m_z, m_tx, m_ty, m_qOverP;
+  };
+
   /** @fn    stateVectorFromComposite
    *  @brief Convert a particle (3-momentum + position) to 4-state (at some z).
    *
@@ -95,8 +186,6 @@ namespace Sel {
 
   /** @fn    stateVectorFromNeutral
    @brief Convert a particle (4-momentum + position) to 4-state (at some z).
-
-   TODO: This must be properly filled
    */
   template <typename particle_t>
   auto stateVectorFromNeutral( particle_t const& p ) {
@@ -105,16 +194,12 @@ namespace Sel {
 
     auto const invpz = 1.f / p.pz();
 
-    State4<float_v> s{};
-    s.x()  = p.hypothesis_position_x();
-    s.y()  = p.hypothesis_position_y();
-    s.z()  = p.hypothesis_position_z();
-    s.tx() = p.px() * invpz;
-    s.ty() = p.py() * invpz;
-
-    s.covXX() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::MatSym<float_v, 2>>();
-    s.covTT() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::MatSym<float_v, 2>>();
-    s.covXT() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::Mat<float_v, 2, 2>>();
+    UnfittedState<float_v> s{p.hypothesis_position_x(),
+                             p.hypothesis_position_y(),
+                             p.hypothesis_position_z(),
+                             p.px() * invpz,
+                             p.py() * invpz,
+                             1. / p.p()};
 
     return s;
   }
diff --git a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
index 400448403c8..d2849f21967 100644
--- a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
+++ b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
@@ -33,9 +33,6 @@ namespace {
   template <std::size_t N, std::size_t M, typename T>
   using Matrix = LHCb::LinAlg::Mat<T, N, M>;
 
-  template <std::size_t N, typename T>
-  using MatrixNxN = LHCb::LinAlg::Mat<T, N>;
-
   // Fast utility function for inverting a 2x2 matrix. This just comes
   // from ROOT/Math/Dinv, but I took out the test on the determinant
   // since we do not need it.
@@ -63,79 +60,110 @@ namespace {
     NumTypes
   };
 
-  template <ChildType C>
-  static constexpr bool has_full_state_v = ( C == ChildType::Composite || C == ChildType::Resonance ||
-                                             C == ChildType::Neutral || C == ChildType::Other );
+  /** Define the type of child based on the properties of a proxy
 
-  /* TODO look at merging this with the stuff in State4.h
+      Note that it is hardly possible to define functions and classes based on the
+      proxy types due to the presence of zip, cv-qualifiers and the high volatility
+      of the containers in their definition. We rely on static member values that
+      determine which operations can done and what information can be accessed.
    */
-  template <typename F>
-  struct MiniState {
-    SymNxN<5, F> cov{};
-    auto&        x() { return m_x; }
-    auto const&  x() const { return m_x; }
-    auto&        y() { return m_y; }
-    auto const&  y() const { return m_y; }
-    auto&        z() { return m_z; }
-    auto const&  z() const { return m_z; }
-    auto&        qOverP() { return m_qOverP; }
-    auto const&  qOverP() const { return m_qOverP; }
-    auto&        tx() { return m_tx; }
-    auto const&  tx() const { return m_tx; }
-    auto&        ty() { return m_ty; }
-    auto const&  ty() const { return m_ty; }
+  template <bool ParticleCondition, bool TrackCondition, bool ExtrapolationCondition>
+  struct child_type_for_conditions;
 
-    template <typename T>
-    MiniState( T const& other )
-        : m_x{other.x()}, m_y{other.y()}, m_tx{other.tx()}, m_ty{other.ty()}, m_z{other.z()}, m_qOverP{other.qOverP()} {
-      static_assert( Sel::Utils::canBeExtrapolatedDownstream_v<T> );
-      auto const& other_cov = other.covariance();
-      LHCb::Utils::unwind<0, 5>( [this, &other_cov]( auto i ) {
-        LHCb::Utils::unwind<i, 5>( [this, i, &other_cov]( auto j ) { cov( i, j ) = other_cov( i, j ); } );
-      } );
-    }
+  template <>
+  struct child_type_for_conditions<true /* particle */, true /* track */, true /* extrapolation */> {
+    static constexpr auto value = ChildType::TrackWithVelo;
+  };
 
-    void linearTransportTo( F const& new_z ) {
-      const auto dz  = new_z - m_z;
-      const auto dz2 = dz * dz;
-      m_x += dz * m_tx;
-      m_y += dz * m_ty;
-      cov( 0, 0 ) += dz2 * cov( 2, 2 ) + 2 * dz * cov( 2, 0 );
-      cov( 1, 0 ) += dz2 * cov( 3, 2 ) + dz * ( cov( 3, 0 ) + cov( 2, 1 ) );
-      cov( 2, 0 ) += dz * cov( 2, 2 );
-      cov( 3, 0 ) += dz * cov( 3, 2 );
-      cov( 4, 0 ) += dz * cov( 4, 2 );
-      cov( 1, 1 ) += dz2 * cov( 3, 3 ) + 2 * dz * cov( 3, 1 );
-      cov( 2, 1 ) += dz * cov( 3, 2 );
-      cov( 3, 1 ) += dz * cov( 3, 3 );
-      cov( 4, 1 ) += dz * cov( 4, 3 );
-      m_z = new_z;
-    }
+  template <>
+  struct child_type_for_conditions<true /* particle */, false /* track */, false /* extrapolation */> {
+    static constexpr auto value = ChildType::Neutral;
+  };
+
+  template <>
+  struct child_type_for_conditions<false /* particle */, false /* track */, false /* extrapolation */> {
+    static constexpr auto value = ChildType::Composite;
+  };
+
+  /// Determine the child type for a given proxy
+  template <class Proxy>
+  struct child_type_for_proxy
+      : child_type_for_conditions<Sel::Utils::isBasicParticle_v<Proxy>, Sel::Utils::hasTrack_v<Proxy>,
+                                  Sel::Utils::canBeExtrapolatedDownstream_v<Proxy>> {};
+
+  /// Child type for a given proxy
+  template <class Proxy>
+  static constexpr auto child_type_for_proxy_v = child_type_for_proxy<Proxy>::value;
+
+  template <typename Proxy>
+  static constexpr auto is_composite_v = child_type_for_proxy_v<Proxy> == ChildType::Composite;
+
+  template <typename Proxy>
+  static constexpr bool has_covariance_v = ( child_type_for_proxy_v<Proxy> == ChildType::Composite ||
+                                             child_type_for_proxy_v<Proxy> == ChildType::Resonance ||
+                                             child_type_for_proxy_v<Proxy> == ChildType::TrackWithVelo ||
+                                             child_type_for_proxy_v<Proxy> == ChildType::TrackWithoutVelo ||
+                                             child_type_for_proxy_v<Proxy> == ChildType::Other );
 
-    SymNxN<2, F>    covXX() const { return cov.template sub<SymNxN<2, F>, 0, 0>(); }
-    MatrixNxN<2, F> covXT() const { return cov.template sub<MatrixNxN<2, F>, 2, 0>(); }
+  /** Basic proxy for a particle and a state (no covariance matrix)
+   */
+  template <typename ProxyType>
+  struct UnfittedParticle final {
+
+  public:
+    using proxy_type = ProxyType;
+    using simd_type  = typename proxy_type::simd_t;
+    using float_v    = typename simd_type::float_v;
+    using state_type = Sel::UnfittedState<float_v>;
 
   private:
-    F m_x{}, m_y{}, m_tx{}, m_ty{}, m_z{}, m_qOverP{};
+    proxy_type m_particle;
+    state_type m_state;
+
+  public:
+    UnfittedParticle( proxy_type p, state_type state ) : m_particle{std::move( p )}, m_state{std::move( state )} {}
+
+    /** Add the momenta of a particle with no covariance information to the result of a combination
+
+        Update the information of the momenta of the objects based on the vertex of
+        the decaying particle. It is crucial that we add the objects after
+        determining the decay vertex (i.e. at the end of the combination).
+     */
+    void addToFourVector( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>&,
+                          Matrix<4, 3, float_v>& ) const {
+
+      auto dx = m_state.x() - vertexpos( 0 );
+      auto dy = m_state.y() - vertexpos( 1 );
+      auto dz = m_state.z() - vertexpos( 2 );
+
+      auto d = sqrt( dx * dx + dy * dy + dz * dz );
+
+      auto p = m_particle.p();
+
+      p4 = VecN<4, float_v>{p * dx / d, p * dy / d, p * dz / d, m_particle.energy()};
+    }
   };
 
   // Class for daughters with a straight-line trajectory
-  template <typename particle_t, ChildType child_type>
-  struct VertexTraj final {
-    using simd_t  = typename particle_t::simd_t;
-    using int_v   = typename simd_t::int_v;
-    using float_v = typename simd_t::float_v;
-    using state_t = std::conditional_t<has_full_state_v<child_type>, Sel::State4<float_v>, MiniState<float_v>>;
+  template <typename ProxyType>
+  struct FittedParticle final {
+
+  public:
+    using proxy_type = ProxyType;
+    using simd_type  = typename proxy_type::simd_t;
+    using int_v      = typename simd_type::int_v;
+    using float_v    = typename simd_type::float_v;
+    using state_type = std::conditional_t<is_composite_v<proxy_type>, Sel::State4<float_v>, Sel::MiniState<float_v>>;
 
   private:
-    particle_t            m_particle; // TODO maybe just cache some information from this?
-    state_t               m_state;    // TODO evaluate what is really needed
+    proxy_type            m_particle; // TODO maybe just cache some information from this?
+    state_type            m_state;    // TODO evaluate what is really needed
     VecN<2, float_v>      m_q;        // predicted/fitted slope (tx,ty)
     SymNxN<2, float_v>    m_G;        // weight matrix of (x,y) of state
     Matrix<2, 3, float_v> m_A;        // projection matrix for vertex position
 
   public:
-    VertexTraj( particle_t p, state_t state )
+    FittedParticle( proxy_type p, state_type state )
         : m_particle{std::move( p )}, m_state{std::move( state )}, m_q{m_state.tx(), m_state.ty()} {}
 
     void project( VecN<3, float_v> const& vertexpos, VecN<3, float_v>& halfDChisqDX, SymNxN<3, float_v>& halfD2ChisqDX2,
@@ -187,6 +215,9 @@ namespace {
 
     void addToFourVector( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>& p4cov,
                           Matrix<4, 3, float_v>& gainmatrix ) const {
+
+      static constexpr auto child_type = child_type_for_proxy_v<proxy_type>;
+
       // Despatch is a bit crude, but it'll do...in the original
       // ParticleVertexFitter this was done with a full specialisation, but
       // now it would have to be a partial specialisation and it all gets a
@@ -195,35 +226,10 @@ namespace {
         addToFourVector_Composite( vertexpos, p4, p4cov, gainmatrix );
       } else if constexpr ( child_type == ChildType::TrackWithVelo || child_type == ChildType::TrackWithoutVelo ) {
         addToFourVector_Track( vertexpos, p4, p4cov, gainmatrix );
-      } else if constexpr ( child_type == ChildType::Neutral ) {
-        addToFourVector_Neutral( vertexpos, p4, p4cov, gainmatrix );
       } else {
       }
     }
 
-    /** Add the momenta of a neutral particle to the result of a combination
-
-        Update the information of the momenta of the neutral objects based on the vertex of
-        the decaying particle. It is crucial that we add the neutral objects after
-        determining the decay vertex (i.e. at the end of the combination).
-
-        @todo FIXME correctly process the covariance matrix
-     */
-    void addToFourVector_Neutral( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>&,
-                                  Matrix<4, 3, float_v>& ) const {
-      // TODO: make this a global function that modifies the neutral objects at the
-      // persistency stage
-      auto dx = m_state.x() - vertexpos( 0 );
-      auto dy = m_state.y() - vertexpos( 1 );
-      auto dz = m_state.z() - vertexpos( 2 );
-
-      auto d = sqrt( dx * dx + dy * dy + dz * dz );
-
-      auto p = m_particle.p();
-
-      p4 = VecN<4, float_v>{p * dx / d, p * dy / d, p * dz / d, m_particle.energy()};
-    }
-
     /** @todo FIXME implement logic for bremsstrahlung photons.
      */
     void addToFourVector_Track( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>& p4cov,
@@ -377,43 +383,16 @@ namespace {
     }
   };
 
-  template <bool ParticleCondition, bool TrackCondition, bool ExtrapolationCondition>
-  struct child_type_for_conditions;
-
-  template <>
-  struct child_type_for_conditions<true, true, true> {
-    static constexpr auto value = ChildType::TrackWithVelo;
-  };
-
-  template <>
-  struct child_type_for_conditions<true, false, false> {
-    static constexpr auto value = ChildType::Neutral;
-  };
-
-  template <>
-  struct child_type_for_conditions<false, false, false> {
-    static constexpr auto value = ChildType::Composite;
-  };
-
-  template <class Proxy>
-  struct child_type_for_proxy
-      : child_type_for_conditions<Sel::Utils::isBasicParticle_v<Proxy>, Sel::Utils::hasTrack_v<Proxy>,
-                                  Sel::Utils::canBeExtrapolatedDownstream_v<Proxy>> {};
-
-  template <class Proxy>
-  static constexpr auto child_type_for_proxy_v = child_type_for_proxy<Proxy>::value;
-
+  /// Define the constructor of child objects for each child type
   template <ChildType>
-  struct ChildHelper;
+  struct prepare_child_t;
 
   /** Helper type for identifying and handling TrackWithVelo and
    *  TrackWithoutVelo children.
    *  @todo Handle TrackWithoutVelo.
    */
   template <>
-  struct ChildHelper<ChildType::TrackWithVelo> {
-
-    static constexpr auto child_type = ChildType::TrackWithVelo;
+  struct prepare_child_t<ChildType::TrackWithVelo> {
 
     /** Given that `value` was `true`, and that `type` was returned by our
      *  `deriveType` function, return the narrowest-possible std::variant of
@@ -428,11 +407,9 @@ namespace {
      *        `state(StateLocation::ClosestToBeam)` accessor. Remove [[maybe_unused]].
      */
     template <typename mask_v, typename float_v, typename proxy>
-    static auto prepareChild( [[maybe_unused]] LHCb::LinAlg::Vec<float_v, 3> const& vertex_position,
-                              [[maybe_unused]] mask_v const& vertex_position_mask, proxy const& child ) {
-      using VertexTrack =
-          VertexTraj<proxy, child_type /* flag this represents a track not a weakly-decaying composite */>;
-      return VertexTrack{child, child.state( LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam )};
+    constexpr auto operator()( [[maybe_unused]] LHCb::LinAlg::Vec<float_v, 3> const& vertex_position,
+                               [[maybe_unused]] mask_v const& vertex_position_mask, proxy const& child ) const {
+      return FittedParticle<proxy>( child, child.state( LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam ) );
     }
   };
 
@@ -442,9 +419,7 @@ namespace {
    *        to go in the variant returned by prepareChild.
    */
   template <>
-  struct ChildHelper<ChildType::Composite> {
-
-    static constexpr auto child_type = ChildType::Composite;
+  struct prepare_child_t<ChildType::Composite> {
 
     /** Given that `value` was `true`, and that `type` was returned by our
      *  `deriveType` function, return the narrowest-possible std::variant of
@@ -459,34 +434,31 @@ namespace {
      *        add another member to the return type; remove [[maybe_unused]]
      */
     template <typename mask_v, typename float_v, typename proxy>
-    static auto prepareChild( LHCb::LinAlg::Vec<float_v, 3> const&, mask_v const&, proxy const& child ) {
-      using VertexComposite =
-          VertexTraj<proxy, child_type /* flag this represents a weakly-decaying composite and not a track */>;
-      return VertexComposite{child, Sel::stateVectorFromComposite( child )};
+    constexpr auto operator()( LHCb::LinAlg::Vec<float_v, 3> const&, mask_v const&, proxy const& child ) const {
+      return FittedParticle<proxy>( child, Sel::stateVectorFromComposite( child ) );
     }
   };
 
   /** Helper type for the default `Neutral` case.
    */
   template <>
-  struct ChildHelper<ChildType::Neutral> {
-
-    static constexpr auto child_type = ChildType::Neutral;
+  struct prepare_child_t<ChildType::Neutral> {
 
     /** @todo Need to implement the helper class that handles generic
      *        children and return an instance of it. As a 1-member variant.
      */
     template <typename mask_v, typename float_v, typename proxy>
-    static auto prepareChild( [[maybe_unused]] LHCb::LinAlg::Vec<float_v, 3> const& vertex_position,
-                              [[maybe_unused]] mask_v const& vertex_position_mask, proxy const& child ) {
-      using VertexNeutral =
-          VertexTraj<proxy, child_type /* flag this represents a track not a weakly-decaying composite */>;
-      return VertexNeutral{child, Sel::stateVectorFromNeutral( child )};
+    constexpr auto operator()( [[maybe_unused]] LHCb::LinAlg::Vec<float_v, 3> const& vertex_position,
+                               [[maybe_unused]] mask_v const& vertex_position_mask, proxy const& child ) const {
+      return UnfittedParticle<proxy>{child, Sel::stateVectorFromNeutral( child )};
     }
   };
 
-  template <class Proxy>
-  using child_helper_for_proxy_t = ChildHelper<child_type_for_proxy_v<Proxy>>;
+  /// Return the object that will handle the information of a particle and its state
+  template <class Proxy, typename... Args>
+  auto prepare_child( Args&&... args ) {
+    return prepare_child_t<child_type_for_proxy_v<Proxy>>{}( std::forward<Args>( args )... );
+  }
 } // namespace
 
 namespace Sel::Fitters {
@@ -559,8 +531,8 @@ namespace Sel::Fitters {
 
     template <typename mask_v, typename float_v, template <class...> class Variant, typename... Proxy>
     struct return_types_helper<mask_v, float_v, Variant<Proxy...>> {
-      using type = std::variant<decltype( child_helper_for_proxy_t<Proxy>::prepareChild(
-          std::declval<VecN<3, float_v>>(), std::declval<mask_v>(), std::declval<Proxy>() ) )...>;
+      using type = std::variant<decltype(
+          prepare_child<Proxy>( std::declval<VecN<3, float_v>>(), std::declval<mask_v>(), std::declval<Proxy>() ) )...>;
     };
 
     template <typename mask_v, typename float_v, typename VariantTypes>
@@ -580,15 +552,10 @@ namespace Sel::Fitters {
       // if the former, this yields tuple<Ts...>, otherwise tuple<child_t>
       using variant_types = typename LHCb::is_variant<child_t>::tuple_type;
 
-      using ret_t = return_types_helper_t<mask_v, float_v, variant_types>;
-
-      static_assert( boost::mp11::mp_is_set<ret_t>::value,
-                     "Got duplicate entries in derived variant type. Looks like the code needs an mp_unique<...>." );
       // Now we know what return type to impose, invoke/visit on `child`
       return LHCb::invoke_or_visit(
-          [&]( auto const& child ) -> ret_t {
-            using helper_t = child_helper_for_proxy_t<std::decay_t<decltype( child )>>;
-            return helper_t::prepareChild( vertex_position, vertex_position_mask, child );
+          [&]( auto const& child ) -> return_types_helper_t<mask_v, float_v, variant_types> {
+            return prepare_child<std::decay_t<decltype( child )>>( vertex_position, vertex_position_mask, child );
           },
           possibly_variant_child );
     }
@@ -800,9 +767,12 @@ namespace Sel::Fitters {
         auto halfDChisqDX   = LHCb::LinAlg::initialize_with_zeros<VecN<3, float_v>>();
 
         // add all particles that contribute to the vertex chi2
-        ( std::visit( [&, &vertex_position = vertex_position](
-                          auto& child ) { child.project( vertex_position, halfDChisqDX, halfD2ChisqDX2, chi2, ndof ); },
-                      std::get<IChild>( extrapolated_children ) ),
+        ( std::visit(
+              [&, &vertex_position = vertex_position]( auto& child ) {
+                if constexpr ( has_covariance_v<typename std::decay_t<decltype( child )>::proxy_type> )
+                  child.project( vertex_position, halfDChisqDX, halfD2ChisqDX2, chi2, ndof );
+              },
+              std::get<IChild>( extrapolated_children ) ),
           ... );
 
         // calculate the covariance and the change in the position
@@ -822,8 +792,12 @@ namespace Sel::Fitters {
         //      again"? Right now then if we hit max_iters we might update the
         //      slopes but not the chi2
         if ( !all( converged ) ) {
-          ( std::visit( [& vertex_position = vertex_position]( auto& child ) { child.updateSlopes( vertex_position ); },
-                        std::get<IChild>( extrapolated_children ) ),
+          ( std::visit(
+                [& vertex_position = vertex_position]( auto& child ) {
+                  if constexpr ( has_covariance_v<typename std::decay_t<decltype( child )>::proxy_type> )
+                    child.updateSlopes( vertex_position );
+                },
+                std::get<IChild>( extrapolated_children ) ),
             ... );
         }
       }
-- 
GitLab


From aac7f64b528a4f0496a0771b849161f3a820bd73 Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Fri, 1 Apr 2022 11:29:17 +0100
Subject: [PATCH 24/26] Fixes since there is no implementation of Resonance
 objects with the new event model yet

---
 .../include/VertexFit/ParticleVertexFitter.h  | 25 +++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
index d2849f21967..48aefc23619 100644
--- a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
+++ b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h
@@ -564,6 +564,9 @@ namespace Sel::Fitters {
      *  1. if there are Resonance daughters, take the vertex position of the first.
      *  2. if not, use the states of the first two TrackWithVelo or Composite
      *  3. if not, try with TrackWithoutVelo and trajpoca
+
+     TODO: Remove the case 1-bis once the differentiation between composites and
+     resonances is done
      */
     template <typename simd_t, std::size_t N, typename Combination>
     std::tuple<VecN<3, typename simd_t::float_v>, typename simd_t::mask_v>
@@ -604,6 +607,28 @@ namespace Sel::Fitters {
         } );
         assert( tmp.has_value() );
         return *tmp; // no validity check, but we have the assert just above
+      } else if ( count( ChildType::TrackWithVelo ) == 0 && count( ChildType::Composite ) == 1 ) {
+        // Case 1-bis
+        std::optional<std::tuple<VecN<3, float_v>, mask_v>> tmp{std::nullopt};
+        LHCb::Utils::unwind<0, N>( [&]( auto i ) {
+          // Find the first `Composite` child
+          if ( !tmp.has_value() && types[i] == ChildType::Composite ) {
+            // Use invoke_or_visit in case the relevant entry in `children` is
+            // a variant type.
+            LHCb::invoke_or_visit(
+                [&]( auto const& child ) {
+                  using Sel::Utils::endVertexPos;
+                  if constexpr ( !Sel::Utils::isBasicParticle_v<decltype( child )> ) {
+                    tmp.emplace( LHCb::LinAlg::convert( endVertexPos( child ) ), mask_v{true} );
+                  } else {
+                    throw exception( "Resonance type didn't have an endVertex()" );
+                  }
+                },
+                children.template get<i>() );
+          }
+        } );
+        assert( tmp.has_value() );
+        return *tmp; // no validity check, but we have the assert just above
       } else if ( count( ChildType::TrackWithVelo ) + count( ChildType::Composite ) >= 2 ) {
         // Case 2
         // Use the POCA of TrackWithVelo + closestToBeamState or Composite + position/momentum
-- 
GitLab


From cd1a067096d14ee5c9ae6c58e226ae976390e192 Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Fri, 1 Apr 2022 17:09:13 +0100
Subject: [PATCH 25/26] Make the pointing strategy an enum-based configurable

---
 Phys/ParticleMaker/src/NeutralMakers.cpp | 26 +++++++++++++++++++-----
 1 file changed, 21 insertions(+), 5 deletions(-)

diff --git a/Phys/ParticleMaker/src/NeutralMakers.cpp b/Phys/ParticleMaker/src/NeutralMakers.cpp
index 1fadde1eb4f..6f6997ebfc4 100644
--- a/Phys/ParticleMaker/src/NeutralMakers.cpp
+++ b/Phys/ParticleMaker/src/NeutralMakers.cpp
@@ -22,6 +22,7 @@
 #include "CaloFutureUtils/CaloMomentum.h"
 #include "Kernel/IParticlePropertySvc.h"
 #include "Kernel/ParticleProperty.h"
+#include "Kernel/meta_enum.h"
 
 // File containing definition AND implementation of neutral makers (photons, pi0s (...))
 /** @class PhotonMaker
@@ -637,6 +638,8 @@ namespace LHCb::Phys::ParticleMakers {
     debug() << "--------------------" << endmsg;
   }
 
+  meta_enum_class( PointingStrategy, int, Unknown = -1, AveragePVs, OneCandidatePerPV );
+
   /** @class NeutralBasicsMaker
    *
    * @brief Create a container of neutral objects from the output of the reconstruction sequence and a set of vertices
@@ -660,7 +663,10 @@ namespace LHCb::Phys::ParticleMakers {
 
       LHCb::Event::NeutralBasics neutrals( unique_id_gen );
 
-      if ( m_AveragePVs.value() ) {
+      switch ( m_PointingStrategy.value() ) {
+
+      case PointingStrategy::AveragePVs: {
+
         // Calculate the average the primary vertex positions
         neutrals.reserve( calo_hypos.size() );
 
@@ -672,21 +678,31 @@ namespace LHCb::Phys::ParticleMakers {
                                     } );
 
         for ( auto const& ch : calo_hypos ) neutrals.emplace_back( unique_id_gen, ch, avg );
-      } else {
+
+        break;
+      }
+      case PointingStrategy::OneCandidatePerPV: {
+
         // Create one particle per primary vertex and ECAL hypothesis
         neutrals.reserve( calo_hypos.size() * primary_vertices.size() );
 
         for ( auto const& ch : calo_hypos )
           for ( auto const& vtx : primary_vertices ) neutrals.emplace_back( unique_id_gen, ch, vtx.position() );
+
+        break;
+      }
+      case PointingStrategy::Unknown:
+        // we will never end-up here since Gaudi::Property will trigger an error beforehand
+        throw std::runtime_error( "Pointing strategy set to \"Unknown\" in NeutralBasicsMaker" );
       }
 
       return neutrals;
     }
 
   private:
-    Gaudi::Property<bool> m_AveragePVs{this, "AveragePVs", true,
-                                       "If set to \"true\"; uses the average of the primary vertices in order to "
-                                       "determine the pointing direction of the neutral particles"};
+    Gaudi::Property<PointingStrategy> m_PointingStrategy{
+        this, "PointingStrategy", PointingStrategy::AveragePVs,
+        "Strategy to adopt in order to determine the direction of the momenta"};
   };
 
   DECLARE_COMPONENT( LHCb::Phys::ParticleMakers::NeutralBasicsMaker )
-- 
GitLab


From d267ae0728355131f194d3630d02237c95add02b Mon Sep 17 00:00:00 2001
From: CI Runner <ci-runner@gitlab.cern>
Date: Fri, 1 Apr 2022 17:16:49 +0100
Subject: [PATCH 26/26] Add a new option to NeutralBasicsMaker

---
 Phys/ParticleMaker/src/NeutralMakers.cpp | 17 ++++++++++++++---
 1 file changed, 14 insertions(+), 3 deletions(-)

diff --git a/Phys/ParticleMaker/src/NeutralMakers.cpp b/Phys/ParticleMaker/src/NeutralMakers.cpp
index 6f6997ebfc4..032040c6ffb 100644
--- a/Phys/ParticleMaker/src/NeutralMakers.cpp
+++ b/Phys/ParticleMaker/src/NeutralMakers.cpp
@@ -638,7 +638,7 @@ namespace LHCb::Phys::ParticleMakers {
     debug() << "--------------------" << endmsg;
   }
 
-  meta_enum_class( PointingStrategy, int, Unknown = -1, AveragePVs, OneCandidatePerPV );
+  meta_enum_class( PointingStrategy, int, Unknown = -1, AveragePVs, OneCandidatePerPV, OriginAsPV );
 
   /** @class NeutralBasicsMaker
    *
@@ -667,7 +667,7 @@ namespace LHCb::Phys::ParticleMakers {
 
       case PointingStrategy::AveragePVs: {
 
-        // Calculate the average the primary vertex positions
+        // calculate the average the primary vertex positions
         neutrals.reserve( calo_hypos.size() );
 
         auto avg = std::accumulate( primary_vertices.cbegin(), primary_vertices.cend(), Gaudi::XYZPoint{0.f, 0.f, 0.f},
@@ -683,7 +683,7 @@ namespace LHCb::Phys::ParticleMakers {
       }
       case PointingStrategy::OneCandidatePerPV: {
 
-        // Create one particle per primary vertex and ECAL hypothesis
+        // create one particle per primary vertex and ECAL hypothesis
         neutrals.reserve( calo_hypos.size() * primary_vertices.size() );
 
         for ( auto const& ch : calo_hypos )
@@ -691,6 +691,17 @@ namespace LHCb::Phys::ParticleMakers {
 
         break;
       }
+      case PointingStrategy::OriginAsPV: {
+
+        // all particles have the (0, 0, 0) as the origin
+        neutrals.reserve( calo_hypos.size() );
+
+        Gaudi::XYZPoint origin{0.f, 0.f, 0.f};
+
+        for ( auto const& ch : calo_hypos ) neutrals.emplace_back( unique_id_gen, ch, origin );
+
+        break;
+      }
       case PointingStrategy::Unknown:
         // we will never end-up here since Gaudi::Property will trigger an error beforehand
         throw std::runtime_error( "Pointing strategy set to \"Unknown\" in NeutralBasicsMaker" );
-- 
GitLab