diff --git a/CI/setup_cvmfs_lcg.sh b/CI/setup_cvmfs_lcg.sh
index 41fc7573dff..de711e0548e 100644
--- a/CI/setup_cvmfs_lcg.sh
+++ b/CI/setup_cvmfs_lcg.sh
@@ -1,4 +1,4 @@
-# setup lowest, supported LCG release via cvmfs
+# setup latest, supported LCG release via cvmfs
 
 if test -n "$BASH_SOURCE"; then
   this_script=$BASH_SOURCE
@@ -11,4 +11,4 @@ else
 fi
 
 dir="$( cd "$( dirname "${this_script}" )" && pwd )"
-source $dir/setup_cvmfs_lcg101.sh
+source $dir/setup_cvmfs_lcg105.sh
diff --git a/CI/setup_cvmfs_lcg100.sh b/CI/setup_cvmfs_lcg100.sh
deleted file mode 100644
index a7bede7f296..00000000000
--- a/CI/setup_cvmfs_lcg100.sh
+++ /dev/null
@@ -1,19 +0,0 @@
-# setup appropriate LCG 100 release via cvmfs
-
-if test -e /etc/centos-release && grep 'CentOS Linux release 7' /etc/centos-release; then
-  lcg_os=centos7
-elif test -e /etc/centos-release && grep 'CentOS Stream release 8' /etc/centos-release; then
-  lcg_os=centos8
-else
-  echo "Unsupported system" 1>&2
-  return
-fi
-
-lcg_release=LCG_100
-lcg_compiler=gcc10-opt
-lcg_platform=x86_64-${lcg_os}-${lcg_compiler}
-lcg_view=/cvmfs/sft.cern.ch/lcg/views/${lcg_release}/${lcg_platform}
-
-source ${lcg_view}/setup.sh
-# extra variables required to build acts
-export DD4hep_DIR=${lcg_view}
diff --git a/CI/setup_cvmfs_lcg101.sh b/CI/setup_cvmfs_lcg101.sh
deleted file mode 100644
index 20b816e2fe0..00000000000
--- a/CI/setup_cvmfs_lcg101.sh
+++ /dev/null
@@ -1,19 +0,0 @@
-# setup appropriate LCG 101 release via cvmfs
-
-if test -e /etc/centos-release && grep 'CentOS Linux release 7' /etc/centos-release; then
-  lcg_os=centos7
-elif test -e /etc/centos-release && grep 'CentOS Stream release 8' /etc/centos-release; then
-  lcg_os=centos8
-else
-  echo "Unsupported system" 1>&2
-  return
-fi
-
-lcg_release=LCG_101
-lcg_compiler=gcc11-opt
-lcg_platform=x86_64-${lcg_os}-${lcg_compiler}
-lcg_view=/cvmfs/sft.cern.ch/lcg/views/${lcg_release}/${lcg_platform}
-
-source ${lcg_view}/setup.sh
-# extra variables required to build acts
-export DD4hep_DIR=${lcg_view}
diff --git a/CI/setup_cvmfs_lcg105.sh b/CI/setup_cvmfs_lcg105.sh
index c45bba4493e..531d8c827fa 100644
--- a/CI/setup_cvmfs_lcg105.sh
+++ b/CI/setup_cvmfs_lcg105.sh
@@ -13,11 +13,10 @@ else
 fi
 
 lcg_release=LCG_105
-lcg_compiler=gcc11-opt
-lcg_platform=x86_64-${lcg_os}-${lcg_compiler}
+lcg_compiler=gcc13
+lcg_platform=x86_64-${lcg_os}-${lcg_compiler}-opt
 lcg_view=/cvmfs/sft.cern.ch/lcg/views/${lcg_release}/${lcg_platform}
 
 source ${lcg_view}/setup.sh
 # extra variables required to build acts
 export DD4hep_DIR=${lcg_view}
-
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 714a9521804..f2bedbc0b40 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -222,10 +222,17 @@ set(_acts_boost_recommended_version 1.78.0)
 # Include the sources for the external dependencies.
 include(ActsExternSources)
 
+# Controls behavior of DOWNLOAD_EXTRACT_TIMESTAMP
+if(POLICY CMP0135)
+  cmake_policy(SET CMP0135 NEW)
+endif()
+
 # required packages
 if (ACTS_SETUP_BOOST)
   if (ACTS_USE_SYSTEM_BOOST)
-    # NOTE FindBoost.cmake looks for BoostConfig.cmake first, before running it's own logic.
+    if(POLICY CMP0167)
+      cmake_policy(SET CMP0167 NEW)
+    endif()
 
     # Enable both program_options and unit_test_framework to reduce complexity
     # Also Cuda tests seem to use program_options
diff --git a/Core/include/Acts/AmbiguityResolution/ScoreBasedAmbiguityResolution.hpp b/Core/include/Acts/AmbiguityResolution/ScoreBasedAmbiguityResolution.hpp
index 2bbb812848b..afc19083117 100644
--- a/Core/include/Acts/AmbiguityResolution/ScoreBasedAmbiguityResolution.hpp
+++ b/Core/include/Acts/AmbiguityResolution/ScoreBasedAmbiguityResolution.hpp
@@ -117,7 +117,7 @@ class ScoreBasedAmbiguityResolution {
   /// good enough, based on some criteria. Users are free to add their own cuts
   /// with the help of this struct.
   template <typename track_container_t, typename traj_t,
-            template <typename> class holder_t, bool ReadOnly>
+            template <typename> class holder_t, bool ReadOnly = true>
   struct OptionalCuts {
     using OptionalFilter =
         std::function<bool(const Acts::TrackProxy<track_container_t, traj_t,
@@ -162,7 +162,7 @@ class ScoreBasedAmbiguityResolution {
   /// @param optionalCuts is the user defined optional cuts to be applied.
   /// @return a vector of scores for each track
   template <typename track_container_t, typename traj_t,
-            template <typename> class holder_t, bool ReadOnly>
+            template <typename> class holder_t, bool ReadOnly = true>
   std::vector<double> simpleScore(
       const TrackContainer<track_container_t, traj_t, holder_t>& tracks,
       const std::vector<std::vector<TrackFeatures>>& trackFeaturesVectors,
@@ -176,7 +176,7 @@ class ScoreBasedAmbiguityResolution {
   /// @param optionalCuts is the user defined optional cuts to be applied.
   /// @return a vector of scores for each track
   template <typename track_container_t, typename traj_t,
-            template <typename> class holder_t, bool ReadOnly>
+            template <typename> class holder_t, bool ReadOnly = true>
   std::vector<double> ambiguityScore(
       const TrackContainer<track_container_t, traj_t, holder_t>& tracks,
       const std::vector<std::vector<TrackFeatures>>& trackFeaturesVectors,
@@ -206,7 +206,7 @@ class ScoreBasedAmbiguityResolution {
   /// @param optionalCuts is the optional cuts to be applied
   /// @return a vector of IDs of the tracks we want to keep
   template <typename track_container_t, typename traj_t,
-            template <typename> class holder_t, bool ReadOnly>
+            template <typename> class holder_t, bool ReadOnly = true>
   std::vector<int> solveAmbiguity(
       const TrackContainer<track_container_t, traj_t, holder_t>& tracks,
       const std::vector<std::vector<MeasurementInfo>>& measurementsPerTrack,
diff --git a/Core/include/Acts/EventData/MultiTrajectory.hpp b/Core/include/Acts/EventData/MultiTrajectory.hpp
index d880177c57c..9a148efeb9d 100644
--- a/Core/include/Acts/EventData/MultiTrajectory.hpp
+++ b/Core/include/Acts/EventData/MultiTrajectory.hpp
@@ -226,8 +226,7 @@ class MultiTrajectory {
   /// @note Only available if the MultiTrajectory is not read-only
   /// @param istate The index to access
   /// @return Read-write proxy to the stored track state
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  TrackStateProxy getTrackState(IndexType istate) {
+  TrackStateProxy getTrackState(IndexType istate) requires(!ReadOnly) {
     return {*this, istate};
   }
 
@@ -238,9 +237,8 @@ class MultiTrajectory {
   ///       which to leave invalid
   /// @param iprevious index of the previous state, kInvalid if first
   /// @return Index of the newly added track state
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
   IndexType addTrackState(TrackStatePropMask mask = TrackStatePropMask::All,
-                          IndexType iprevious = kInvalid) {
+                          IndexType iprevious = kInvalid) requires(!ReadOnly) {
     return self().addTrackState_impl(mask, iprevious);
   }
 
@@ -248,10 +246,9 @@ class MultiTrajectory {
   /// This effectively calls @c addTrackState and @c getTrackState
   /// @note Only available if the track state container is not read-only
   /// @return a track state proxy to the newly added track state
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
   TrackStateProxy makeTrackState(
       TrackStatePropMask mask = TrackStatePropMask::All,
-      IndexType iprevious = kInvalid) {
+      IndexType iprevious = kInvalid) requires(!ReadOnly) {
     return getTrackState(addTrackState(mask, iprevious));
   }
 
@@ -276,8 +273,8 @@ class MultiTrajectory {
   /// @warning If the trajectory contains multiple components with common
   ///          points, this can have an impact on the other components.
   /// @note Only available if the MultiTrajectory is not read-only
-  template <typename F, bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  void applyBackwards(IndexType iendpoint, F&& callable) {
+  template <typename F>
+  void applyBackwards(IndexType iendpoint, F&& callable) requires(!ReadOnly) {
     static_assert(detail_lt::VisitorConcept<F, TrackStateProxy>,
                   "Callable needs to satisfy VisitorConcept");
 
@@ -327,8 +324,7 @@ class MultiTrajectory {
   /// @param iendpoint Trajectory entry point to start from
   /// @return Iterator pair to iterate over
   /// @note Mutable version
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  auto reverseTrackStateRange(IndexType iendpoint) {
+  auto reverseTrackStateRange(IndexType iendpoint) requires(!ReadOnly) {
     using range_t =
         detail_lt::TrackStateRange<true, Derived, MeasurementSizeMax, false>;
     if (iendpoint == kInvalid) {
@@ -360,8 +356,7 @@ class MultiTrajectory {
   /// @param istartpoint Trajectory state index for the innermost track
   ///        state to start from
   /// @return Iterator pair to iterate over
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  auto forwardTrackStateRange(IndexType istartpoint) {
+  auto forwardTrackStateRange(IndexType istartpoint) requires(!ReadOnly) {
     using range_t =
         detail_lt::TrackStateRange<false, Derived, MeasurementSizeMax, false>;
     if (istartpoint == kInvalid) {
@@ -386,8 +381,8 @@ class MultiTrajectory {
   /// @note This takes a string argument rather than a hashed string to maintain
   ///       compatibility with backends.
   /// @note Only available if the MultiTrajectory is not read-only
-  template <typename T, bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  void addColumn(std::string_view key) {
+  template <typename T>
+  void addColumn(std::string_view key) requires(!ReadOnly) {
     self().template addColumn_impl<T>(key);
   }
 
@@ -400,10 +395,7 @@ class MultiTrajectory {
 
   /// Clear the @c MultiTrajectory. Leaves the underlying storage untouched
   /// @note Only available if the MultiTrajectory is not read-only
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  void clear() {
-    self().clear_impl();
-  }
+  void clear() requires(!ReadOnly) { self().clear_impl(); }
 
   /// Returns the number of track states contained
   /// @return The number of track states
@@ -432,8 +424,8 @@ class MultiTrajectory {
   /// Retrieve a parameter proxy instance for parameters at a given index
   /// @param parIdx Index into the parameter column
   /// @return Mutable proxy
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  typename TrackStateProxy::Parameters parameters(IndexType parIdx) {
+  typename TrackStateProxy::Parameters parameters(IndexType parIdx) requires(
+      !ReadOnly) {
     return self().parameters_impl(parIdx);
   }
 
@@ -447,8 +439,8 @@ class MultiTrajectory {
   /// Retrieve a covariance proxy instance for a covariance at a given index
   /// @param covIdx Index into the covariance column
   /// @return Mutable proxy
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  typename TrackStateProxy::Covariance covariance(IndexType covIdx) {
+  typename TrackStateProxy::Covariance covariance(IndexType covIdx) requires(
+      !ReadOnly) {
     return self().covariance_impl(covIdx);
   }
 
@@ -462,8 +454,8 @@ class MultiTrajectory {
   /// Retrieve a jacobian proxy instance for a jacobian at a given index
   /// @param istate The track state
   /// @return Mutable proxy
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  typename TrackStateProxy::Covariance jacobian(IndexType istate) {
+  typename TrackStateProxy::Covariance jacobian(IndexType istate) requires(
+      !ReadOnly) {
     return self().jacobian_impl(istate);
   }
 
@@ -479,10 +471,9 @@ class MultiTrajectory {
   /// @tparam measdim the measurement dimension
   /// @param istate The track state
   /// @return Mutable proxy
-  template <std::size_t measdim, bool RO = ReadOnly,
-            typename = std::enable_if_t<!RO>>
+  template <std::size_t measdim>
   typename TrackStateProxy::template Calibrated<measdim> calibrated(
-      IndexType istate) {
+      IndexType istate) requires(!ReadOnly) {
     return self().template calibrated_impl<measdim>(istate);
   }
 
@@ -502,16 +493,14 @@ class MultiTrajectory {
   /// @tparam measdim the measurement dimension
   /// @param istate The track state
   /// @return Mutable proxy
-  template <std::size_t measdim, bool RO = ReadOnly,
-            typename = std::enable_if_t<!RO>>
+  template <std::size_t measdim>
   typename TrackStateProxy::template CalibratedCovariance<measdim>
-  calibratedCovariance(IndexType istate) {
+  calibratedCovariance(IndexType istate) requires(!ReadOnly) {
     return self().template calibratedCovariance_impl<measdim>(istate);
   }
 
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
   typename TrackStateProxy::EffectiveCalibrated effectiveCalibrated(
-      IndexType istate) {
+      IndexType istate) requires(!ReadOnly) {
     // This abuses an incorrectly sized vector / matrix to access the
     // data pointer! This works (don't use the matrix as is!), but be
     // careful!
@@ -528,9 +517,8 @@ class MultiTrajectory {
         calibrated<eBoundSize>(istate).data(), calibratedSize(istate)};
   }
 
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
   typename TrackStateProxy::EffectiveCalibratedCovariance
-  effectiveCalibratedCovariance(IndexType istate) {
+  effectiveCalibratedCovariance(IndexType istate) requires(!ReadOnly) {
     // This abuses an incorrectly sized vector / matrix to access the
     // data pointer! This works (don't use the matrix as is!), but be
     // careful!
@@ -576,18 +564,16 @@ class MultiTrajectory {
   ///       or projector. See @c TrackStatePropMask.
   /// @note The track states both need to be stored in the
   ///       same @c MultiTrajectory instance
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
   void shareFrom(IndexType iself, IndexType iother,
                  TrackStatePropMask shareSource,
-                 TrackStatePropMask shareTarget) {
+                 TrackStatePropMask shareTarget) requires(!ReadOnly) {
     self().shareFrom_impl(iself, iother, shareSource, shareTarget);
   }
 
   /// Unset an optional track state component
   /// @param target The component to unset
   /// @param istate The track state index to operate on
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  void unset(TrackStatePropMask target, IndexType istate) {
+  void unset(TrackStatePropMask target, IndexType istate) requires(!ReadOnly) {
     self().unset_impl(target, istate);
   }
 
@@ -595,8 +581,8 @@ class MultiTrajectory {
   /// @note Only available if the track state container is not read-only
   /// @param istate The track state index to alter
   /// @param mask The bitmask that instructs which components to allocate
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  void addTrackStateComponents(IndexType istate, TrackStatePropMask mask) {
+  void addTrackStateComponents(IndexType istate,
+                               TrackStatePropMask mask) requires(!ReadOnly) {
     self().addTrackStateComponents_impl(istate, mask);
   }
 
@@ -605,9 +591,8 @@ class MultiTrajectory {
   /// @tparam key String key for the component to access
   /// @param istate The track state index to operate on
   /// @return Mutable reference to the component given by @p key
-  template <typename T, HashedString key, bool RO = ReadOnly,
-            typename = std::enable_if_t<!RO>>
-  T& component(IndexType istate) {
+  template <typename T, HashedString key>
+  T& component(IndexType istate) requires(!ReadOnly) {
     assert(checkOptional(key, istate));
     return *std::any_cast<T*>(self().component_impl(key, istate));
   }
@@ -617,8 +602,8 @@ class MultiTrajectory {
   /// @param key String key for the component to access
   /// @param istate The track state index to operate on
   /// @return Mutable reference to the component given by @p key
-  template <typename T, bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  T& component(HashedString key, IndexType istate) {
+  template <typename T>
+  T& component(HashedString key, IndexType istate) requires(!ReadOnly) {
     assert(checkOptional(key, istate));
     return *std::any_cast<T*>(self().component_impl(key, istate));
   }
@@ -657,8 +642,8 @@ class MultiTrajectory {
     self().allocateCalibrated_impl(istate, measdim);
   }
 
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  void setUncalibratedSourceLink(IndexType istate, SourceLink sourceLink) {
+  void setUncalibratedSourceLink(IndexType istate,
+                                 SourceLink sourceLink) requires(!ReadOnly) {
     self().setUncalibratedSourceLink_impl(istate, std::move(sourceLink));
   }
 
@@ -670,15 +655,16 @@ class MultiTrajectory {
     return self().referenceSurface_impl(istate);
   }
 
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  void setReferenceSurface(IndexType istate,
-                           std::shared_ptr<const Surface> surface) {
+  void setReferenceSurface(
+      IndexType istate,
+      std::shared_ptr<const Surface> surface) requires(!ReadOnly) {
     self().setReferenceSurface_impl(istate, std::move(surface));
   }
 
  private:
-  template <typename T, bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  void copyDynamicFrom(IndexType dstIdx, const T& src, IndexType srcIdx) {
+  template <typename T>
+  void copyDynamicFrom(IndexType dstIdx, const T& src,
+                       IndexType srcIdx) requires(!ReadOnly) {
     const auto& dynamicKeys = src.self().dynamicKeys_impl();
     for (const auto key : dynamicKeys) {
       std::any srcPtr = src.self().component_impl(key, srcIdx);
diff --git a/Core/include/Acts/EventData/TrackContainer.hpp b/Core/include/Acts/EventData/TrackContainer.hpp
index 7c42a9d18f0..d2cf6c24ec2 100644
--- a/Core/include/Acts/EventData/TrackContainer.hpp
+++ b/Core/include/Acts/EventData/TrackContainer.hpp
@@ -12,6 +12,7 @@
 #include "Acts/Definitions/TrackParametrization.hpp"
 #include "Acts/Definitions/Units.hpp"
 #include "Acts/EventData/MultiTrajectory.hpp"
+#include "Acts/EventData/MultiTrajectoryBackendConcept.hpp"
 #include "Acts/EventData/TrackContainerBackendConcept.hpp"
 #include "Acts/EventData/TrackProxy.hpp"
 #include "Acts/EventData/Types.hpp"
@@ -36,7 +37,8 @@ struct IsReadOnlyTrackContainer;
 /// @tparam track_container_t the track container backend
 /// @tparam traj_t the track state container backend
 /// @tparam holder_t ownership management class for the backend
-template <TrackContainerBackend track_container_t, typename traj_t,
+template <TrackContainerBackend track_container_t,
+          CommonMultiTrajectoryBackend traj_t,
           template <typename> class holder_t = detail::RefHolder>
 class TrackContainer {
  public:
@@ -105,10 +107,8 @@ class TrackContainer {
   ///       You need to ensure suitable lifetime
   /// @param container the track container backend
   /// @param traj the track state container backend
-  template <template <typename> class H = holder_t,
-            typename = std::enable_if_t<
-                detail::is_same_template<H, detail::RefHolder>::value>>
-  TrackContainer(track_container_t& container, traj_t& traj)
+  TrackContainer(auto& container, auto& traj) requires(
+      detail::is_same_template<holder_t, detail::RefHolder>::value)
       : m_container{&container}, m_traj{&traj} {}
 
   /// Constructor from const references to a track container backend and to a
@@ -117,13 +117,9 @@ class TrackContainer {
   ///       You need to ensure suitable lifetime
   /// @param container the track container backend
   /// @param traj the track state container backend
-  template <
-      template <typename> class H = holder_t,
-      bool RO = (IsReadOnlyTrackContainer<track_container_t>::value &&
-                 IsReadOnlyMultiTrajectory<traj_t>::value),
-      typename = std::enable_if_t<
-          detail::is_same_template<H, detail::ConstRefHolder>::value && RO>>
-  TrackContainer(const track_container_t& container, const traj_t& traj)
+  TrackContainer(const auto& container, const auto& traj) requires(
+      detail::is_same_template<holder_t, detail::ConstRefHolder>::value&&
+          ReadOnly&& TrackStateReadOnly)
       : m_container{&container}, m_traj{&traj} {}
 
   /// @}
@@ -147,8 +143,7 @@ class TrackContainer {
   /// @note Only available if the track container is not read-only
   /// @param itrack the track index in the container
   /// @return A mutable track proxy for the index
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  TrackProxy getTrack(IndexType itrack) {
+  TrackProxy getTrack(IndexType itrack) requires(!ReadOnly) {
     return {*this, itrack};
   }
 
@@ -156,8 +151,7 @@ class TrackContainer {
   /// allocates memory. You can combine this with @c getTrack to obtain a track proxy
   /// @note Only available if the track container is not read-only
   /// @return the index to the newly added track
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  IndexType addTrack() {
+  IndexType addTrack() requires(!ReadOnly) {
     auto track = getTrack(m_container->addTrack_impl());
     track.tipIndex() = kInvalid;
     return track.index();
@@ -167,8 +161,7 @@ class TrackContainer {
   /// This effectively calls @c addTrack and @c getTrack
   /// @note Only available if the track container is not read-only
   /// @return a track proxy to the newly added track
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  TrackProxy makeTrack() {
+  TrackProxy makeTrack() requires(!ReadOnly) {
     return getTrack(addTrack());
   }
 
@@ -177,16 +170,14 @@ class TrackContainer {
   /// @note This invalidates track proxies that point to tracks with larger
   ///       indices than @p itrack!
   /// @param itrack The index of the track to remove
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  void removeTrack(IndexType itrack) {
+  void removeTrack(IndexType itrack) requires(!ReadOnly) {
     m_container->removeTrack_impl(itrack);
   }
 
   /// Get a mutable iterator to the first track in the container
   /// @note Only available if the track container is not read-only
   /// @return a mutable iterator to the first track
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  auto begin() {
+  auto begin() requires(!ReadOnly) {
     return detail_tc::TrackProxyIterator<std::decay_t<decltype(*this)>,
                                          TrackProxy, false>{*this, 0};
   }
@@ -194,8 +185,7 @@ class TrackContainer {
   /// Get a past-the-end iterator for this container
   /// @note Only available if the track container is not read-only
   /// @return a past-the-end iterator
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  auto end() {
+  auto end() requires(!ReadOnly) {
     return detail_tc::TrackProxyIterator<std::decay_t<decltype(*this)>,
                                          TrackProxy, false>{*this, size()};
   }
@@ -226,8 +216,8 @@ class TrackContainer {
   /// Add a dymanic column to the track container
   /// @note Only available if the track container is not read-only
   /// @param key the name of the column to be added
-  template <typename T, bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  constexpr void addColumn(std::string_view key) {
+  template <typename T>
+  constexpr void addColumn(std::string_view key) requires(!ReadOnly) {
     m_container->template addColumn_impl<T>(key);
   }
 
@@ -250,9 +240,9 @@ class TrackContainer {
   /// @note Only available if the track container is not read-only
   /// @tparam other_track_container_t Type of the other track container
   /// @param other The other track container
-  template <typename other_track_container_t, bool RO = ReadOnly,
-            typename = std::enable_if_t<!RO>>
-  void ensureDynamicColumns(const other_track_container_t& other) {
+  template <typename other_track_container_t>
+  void ensureDynamicColumns(const other_track_container_t& other) requires(
+      !ReadOnly) {
     container().ensureDynamicColumns_impl(other.container());
   }
 
@@ -267,8 +257,7 @@ class TrackContainer {
   /// Get a mutable reference to the track container backend
   /// @note Only available if the track container is not read-only
   /// @return a mutable reference to the backend
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  auto& container() {
+  auto& container() requires(!ReadOnly) {
     return *m_container;
   }
 
@@ -281,16 +270,14 @@ class TrackContainer {
   /// Get a mutable reference to the track state container backend
   /// @note Only available if the track container is not read-only
   /// @return a mutable reference to the backend
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  auto& trackStateContainer() {
+  auto& trackStateContainer() requires(!ReadOnly) {
     return *m_traj;
   }
 
   /// Retrieve the holder of the track state container
   /// @return The track state container including it's holder
   /// @note Only available if the track container is not read-only
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  auto& trackStateContainerHolder() {
+  auto& trackStateContainerHolder() requires(!ReadOnly) {
     return m_traj;
   }
 
@@ -316,21 +303,20 @@ class TrackContainer {
 
   /// Clear the content of the track container
   /// @note Only available if the track container is not read-only
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  void clear() {
+  void clear() requires(!ReadOnly) {
     m_container->clear();
     m_traj->clear();
   }
 
  protected:
-  template <typename T, HashedString key, bool RO = ReadOnly,
-            typename = std::enable_if_t<!RO>>
-  constexpr T& component(IndexType itrack) {
+  template <typename T, HashedString key>
+  constexpr T& component(IndexType itrack) requires(!ReadOnly) {
     return *std::any_cast<T*>(container().component_impl(key, itrack));
   }
 
-  template <typename T, bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  constexpr T& component(HashedString key, IndexType itrack) {
+  template <typename T>
+  constexpr T& component(HashedString key,
+                         IndexType itrack) requires(!ReadOnly) {
     return *std::any_cast<T*>(container().component_impl(key, itrack));
   }
 
@@ -344,8 +330,8 @@ class TrackContainer {
     return *std::any_cast<const T*>(container().component_impl(key, itrack));
   }
 
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  constexpr typename TrackProxy::Parameters parameters(IndexType itrack) {
+  constexpr typename TrackProxy::Parameters parameters(
+      IndexType itrack) requires(!ReadOnly) {
     return container().parameters(itrack);
   }
 
@@ -354,8 +340,8 @@ class TrackContainer {
     return container().parameters(itrack);
   }
 
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  constexpr typename TrackProxy::Covariance covariance(IndexType itrack) {
+  constexpr typename TrackProxy::Covariance covariance(
+      IndexType itrack) requires(!ReadOnly) {
     return container().covariance(itrack);
   }
 
@@ -364,8 +350,7 @@ class TrackContainer {
     return container().covariance(itrack);
   }
 
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  auto reverseTrackStateRange(IndexType itrack) {
+  auto reverseTrackStateRange(IndexType itrack) requires(!ReadOnly) {
     auto tip = component<IndexType, hashString("tipIndex")>(itrack);
     return m_traj->reverseTrackStateRange(tip);
   }
@@ -375,8 +360,7 @@ class TrackContainer {
     return m_traj->reverseTrackStateRange(tip);
   }
 
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  auto forwardTrackStateRange(IndexType itrack) {
+  auto forwardTrackStateRange(IndexType itrack) requires(!ReadOnly) {
     auto stem = component<IndexType, hashString("stemIndex")>(itrack);
     if (stem == kInvalid) {
       throw std::invalid_argument{"Track has no stem index"};
@@ -393,8 +377,9 @@ class TrackContainer {
   }
 
  private:
-  template <typename T, bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  void copyDynamicFrom(IndexType dstIdx, const T& src, IndexType srcIdx) {
+  template <typename T>
+  void copyDynamicFrom(IndexType dstIdx, const T& src,
+                       IndexType srcIdx) requires(!ReadOnly) {
     const auto& dynamicKeys = src.dynamicKeys_impl();
     for (const auto key : dynamicKeys) {
       std::any srcPtr = src.component_impl(key, srcIdx);
diff --git a/Core/include/Acts/EventData/TrackProxy.hpp b/Core/include/Acts/EventData/TrackProxy.hpp
index 09c511c0522..e7b9e56f89b 100644
--- a/Core/include/Acts/EventData/TrackProxy.hpp
+++ b/Core/include/Acts/EventData/TrackProxy.hpp
@@ -10,6 +10,7 @@
 
 #include "Acts/Definitions/TrackParametrization.hpp"
 #include "Acts/EventData/MultiTrajectory.hpp"
+#include "Acts/EventData/MultiTrajectoryBackendConcept.hpp"
 #include "Acts/EventData/ParticleHypothesis.hpp"
 #include "Acts/EventData/TrackContainerBackendConcept.hpp"
 #include "Acts/EventData/TrackParameters.hpp"
@@ -22,7 +23,8 @@
 
 namespace Acts {
 
-template <Acts::TrackContainerBackend track_container_t, typename traj_t,
+template <TrackContainerBackend track_container_t,
+          CommonMultiTrajectoryBackend traj_t,
           template <typename> class holder_t>
 class TrackContainer;
 
@@ -44,13 +46,14 @@ class TrackProxyIterator {
   using pointer = void;
   using reference = void;
 
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  TrackProxyIterator(container_t& container, IndexType itrack)
+  TrackProxyIterator(container_t& container,
+                     IndexType itrack) requires(!ReadOnly)
       : m_container(&container), m_itrack(itrack) {}
 
-  template <bool RO = ReadOnly, typename = std::enable_if_t<RO>>
-  TrackProxyIterator(const container_t& container, IndexType itrack)
-      : m_container(&container), m_itrack(itrack) {}
+  TrackProxyIterator(const container_t& container,
+                     IndexType itrack) requires ReadOnly
+      : m_container(&container),
+        m_itrack(itrack) {}
 
   TrackProxyIterator& operator++() {
     m_itrack++;
@@ -87,8 +90,7 @@ class TrackProxyIterator {
 
   ProxyType operator*() const { return m_container->getTrack(m_itrack); }
 
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  ProxyType operator*() {
+  ProxyType operator*() requires(!ReadOnly) {
     return m_container->getTrack(m_itrack);
   }
 
@@ -279,8 +281,7 @@ class TrackProxy {
   /// track container
   /// @note Only available if the track proxy is not read-only
   /// @return mutable reference to the tip index
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  IndexType& tipIndex() {
+  IndexType& tipIndex() requires(!ReadOnly) {
     return component<IndexType>(hashString("tipIndex"));
   }
 
@@ -289,8 +290,7 @@ class TrackProxy {
   /// forward-linked.
   /// @note Only available if the track proxy is not read-only
   /// @return mutable reference to the stem index
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  IndexType& stemIndex() {
+  IndexType& stemIndex() requires(!ReadOnly) {
     return component<IndexType>(hashString("stemIndex"));
   }
 
@@ -304,8 +304,8 @@ class TrackProxy {
   // looks like a false-positive. clang-tidy believes `srf` is not movable.
   /// Set a new reference surface for this track
   /// @param srf The surface to set
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  void setReferenceSurface(std::shared_ptr<const Surface> srf) {
+  void setReferenceSurface(std::shared_ptr<const Surface> srf) requires(
+      !ReadOnly) {
     m_container->container().setReferenceSurface_impl(m_index, std::move(srf));
   }
   // NOLINTEND(performance-unnecessary-value-param)
@@ -335,8 +335,7 @@ class TrackProxy {
   /// Mutable version
   /// @note Only available if the track proxy is not read-only
   /// @return Proxy vector for the parameters
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  Parameters parameters() {
+  Parameters parameters() requires(!ReadOnly) {
     return m_container->parameters(m_index);
   }
 
@@ -344,8 +343,7 @@ class TrackProxy {
   /// Mutable version
   /// @note Only available if the track proxy is not read-only
   /// @return Proxy matrix for the covariance
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  Covariance covariance() {
+  Covariance covariance() requires(!ReadOnly) {
     return m_container->covariance(m_index);
   }
 
@@ -394,8 +392,8 @@ class TrackProxy {
   /// Set a new particle hypothesis for this track
   /// @note Only available if the track proxy is not read-only
   /// @param particleHypothesis The particle hypothesis to set
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  void setParticleHypothesis(const ParticleHypothesis& particleHypothesis) {
+  void setParticleHypothesis(
+      const ParticleHypothesis& particleHypothesis) requires(!ReadOnly) {
     m_container->container().setParticleHypothesis_impl(m_index,
                                                         particleHypothesis);
   }
@@ -450,8 +448,7 @@ class TrackProxy {
   /// Return the number of measurements for the track. Const version
   /// @note Only available if the track proxy is not read-only
   /// @return The number of measurements
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  unsigned int& nMeasurements() {
+  unsigned int& nMeasurements() requires(!ReadOnly) {
     return component<unsigned int, hashString("nMeasurements")>();
   }
 
@@ -466,8 +463,7 @@ class TrackProxy {
   /// Mutable version
   /// @note Only available if the track proxy is not read-only
   /// @return The number of holes
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  unsigned int& nHoles() {
+  unsigned int& nHoles() requires(!ReadOnly) {
     return component<unsigned int, hashString("nHoles")>();
   }
 
@@ -481,8 +477,7 @@ class TrackProxy {
   /// Mutable version
   /// @note Only available if the track proxy is not read-only
   /// @return The number of outliers
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  unsigned int& nOutliers() {
+  unsigned int& nOutliers() requires(!ReadOnly) {
     return component<unsigned int, hashString("nOutliers")>();
   }
 
@@ -496,8 +491,7 @@ class TrackProxy {
   /// Mutable version
   /// @note Only available if the track proxy is not read-only
   /// @return The number of shared hits
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  unsigned int& nSharedHits() {
+  unsigned int& nSharedHits() requires(!ReadOnly) {
     return component<unsigned int, hashString("nSharedHits")>();
   }
 
@@ -511,8 +505,7 @@ class TrackProxy {
   /// Mutable version
   /// @note Only available if the track proxy is not read-only
   /// @return The chi squared
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  float& chi2() {
+  float& chi2() requires(!ReadOnly) {
     return component<float, hashString("chi2")>();
   }
 
@@ -526,8 +519,7 @@ class TrackProxy {
   /// track. Mutable version
   /// @note Only available if the track proxy is not read-only
   /// @return The number of degrees of freedom
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  unsigned int& nDoF() {
+  unsigned int& nDoF() requires(!ReadOnly) {
     return component<unsigned int, hashString("ndf")>();
   }
 
@@ -559,8 +551,7 @@ class TrackProxy {
 
   /// Return a mutable track state proxy to the outermost track state
   /// @return The outermost track state proxy
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  TrackStateProxy outermostTrackState() const {
+  TrackStateProxy outermostTrackState() requires(!ReadOnly) {
     return m_container->trackStateContainer().getTrackState(tipIndex());
   }
 
@@ -584,8 +575,7 @@ class TrackProxy {
   /// @note This is only available, if the track is forward linked
   /// @note Only available if the track proxy is not read-only
   /// @return The innermost track state proxy
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  auto innermostTrackState() {
+  auto innermostTrackState() requires(!ReadOnly) {
     using proxy_t = decltype(m_container->trackStateContainer().getTrackState(
         std::declval<IndexType>()));
 
@@ -611,8 +601,7 @@ class TrackProxy {
   /// @note Only available if the track proxy is not read-only
   /// @note This range is from the outside inwards!
   /// @return Track state range to iterate over
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  auto trackStatesReversed() {
+  auto trackStatesReversed() requires(!ReadOnly) {
     return m_container->reverseTrackStateRange(m_index);
   }
 
@@ -637,8 +626,7 @@ class TrackProxy {
   /// @warning This access direction is only possible if the track states are
   ///          **forward-linked**.
   /// @return Track state range to iterate over
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  auto trackStates() {
+  auto trackStates() requires(!ReadOnly) {
     return m_container->forwardTrackStateRange(m_index);
   }
 
@@ -652,8 +640,7 @@ class TrackProxy {
   /// Forward connect a track.
   /// This means setting indices from the inside out on all track states.
   /// @note Only available if the track proxy is not read-only
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  void linkForward() {
+  void linkForward() requires(!ReadOnly) {
     IndexType last = kInvalid;
     for (auto ts : trackStatesReversed()) {
       ts.template component<IndexType>(hashString("next")) = last;
@@ -668,8 +655,8 @@ class TrackProxy {
   /// @note Only available if the track proxy is not read-only
   /// @param mask The allocation prop mask for the new track state
   /// @return The newly added track state
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  auto appendTrackState(TrackStatePropMask mask = TrackStatePropMask::All) {
+  auto appendTrackState(
+      TrackStatePropMask mask = TrackStatePropMask::All) requires(!ReadOnly) {
     auto& tsc = m_container->trackStateContainer();
     auto ts = tsc.makeTrackState(mask, tipIndex());
     tipIndex() = ts.index();
@@ -681,9 +668,9 @@ class TrackProxy {
   /// @tparam track_proxy_t the other track proxy's type
   /// @param other The track proxy
   /// @param copyTrackStates Copy the track state sequence from @p other
-  template <typename track_proxy_t, bool RO = ReadOnly,
-            typename = std::enable_if_t<!RO>>
-  void copyFrom(const track_proxy_t& other, bool copyTrackStates = true) {
+  template <typename track_proxy_t>
+  void copyFrom(const track_proxy_t& other,
+                bool copyTrackStates = true) requires(!ReadOnly) {
     // @TODO: Add constraint on which track proxies are allowed,
     // this is only implicit right now
 
@@ -702,12 +689,14 @@ class TrackProxy {
       reverseTrackStates();
     }
 
-    parameters() = other.parameters();
-    covariance() = other.covariance();
     setParticleHypothesis(other.particleHypothesis());
+
     if (other.hasReferenceSurface()) {
       setReferenceSurface(other.referenceSurface().getSharedPtr());
+      parameters() = other.parameters();
+      covariance() = other.covariance();
     }
+
     nMeasurements() = other.nMeasurements();
     nHoles() = other.nHoles();
     nOutliers() = other.nOutliers();
@@ -719,6 +708,17 @@ class TrackProxy {
                                  other.m_index);
   }
 
+  /// Creates  a *shallow copy* of the track. Track states are not copied, but
+  /// the resulting track points at the same track states as the original.
+  /// @note Only available if the track proxy is not read-only
+  TrackProxy shallowCopy() requires(!ReadOnly) {
+    auto ts = container().makeTrack();
+    ts.copyFrom(*this, false);
+    ts.tipIndex() = tipIndex();
+    ts.stemIndex() = stemIndex();
+    return ts;
+  }
+
   /// Reverse the ordering of track states for this track
   /// Afterwards, the previous endpoint of the track state sequence will be the
   /// "innermost" track state
@@ -726,8 +726,7 @@ class TrackProxy {
   /// @note This is dangerous with branching track state sequences, as it will break them
   /// @note This also automatically forward-links the track!
   /// @param invertJacobians Whether to invert the Jacobians of the track states
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  void reverseTrackStates(bool invertJacobians = false) {
+  void reverseTrackStates(bool invertJacobians = false) requires(!ReadOnly) {
     IndexType current = tipIndex();
     IndexType next = kInvalid;
     IndexType prev = kInvalid;
@@ -779,9 +778,8 @@ class TrackProxy {
   /// @tparam T The type of the component to access
   /// @tparam key String key for the component to access
   /// @return Mutable reference to the component given by @p key
-  template <typename T, HashedString key, bool RO = ReadOnly,
-            typename = std::enable_if_t<!RO>>
-  constexpr T& component() {
+  template <typename T, HashedString key>
+  constexpr T& component() requires(!ReadOnly) {
     return m_container->template component<T, key>(m_index);
   }
 
@@ -789,8 +787,8 @@ class TrackProxy {
   /// @tparam T The type of the component to access
   /// @param key String key for the component to access
   /// @return Mutable reference to the component given by @p key
-  template <typename T, bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  constexpr T& component(HashedString key) {
+  template <typename T>
+  constexpr T& component(HashedString key) requires(!ReadOnly) {
     return m_container->template component<T>(key, m_index);
   }
 
@@ -799,8 +797,8 @@ class TrackProxy {
   /// @param key String key for the component to access
   /// @note This might hash the @p key at runtime instead of compile-time
   /// @return Mutable reference to the component given by @p key
-  template <typename T, bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  constexpr T& component(std::string_view key) {
+  template <typename T>
+  constexpr T& component(std::string_view key) requires(!ReadOnly) {
     return m_container->template component<T>(hashString(key), m_index);
   }
 
@@ -855,8 +853,7 @@ class TrackProxy {
   /// Return a reference to the track container backend, mutable version.
   /// @note Only available if the track proxy is not read-only
   /// @return reference to the track container backend
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  auto& container() {
+  auto& container() requires(!ReadOnly) {
     return *m_container;
   }
 
diff --git a/Core/include/Acts/EventData/TrackStateProxy.hpp b/Core/include/Acts/EventData/TrackStateProxy.hpp
index 505a43cf246..f0ebe42cdd9 100644
--- a/Core/include/Acts/EventData/TrackStateProxy.hpp
+++ b/Core/include/Acts/EventData/TrackStateProxy.hpp
@@ -312,8 +312,7 @@ class TrackStateProxy {
   /// the track sequence
   /// @note Only available if the track state proxy is not read-only
   /// @return The index of the previous track state.
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  IndexType& previous() {
+  IndexType& previous() requires(!ReadOnly) {
     return component<IndexType, hashString("previous")>();
   }
 
@@ -331,16 +330,14 @@ class TrackStateProxy {
   /// Unset an optional track state component
   /// @note Only available if the track state proxy is not read-only
   /// @param target The component to unset
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  void unset(TrackStatePropMask target) {
+  void unset(TrackStatePropMask target) requires(!ReadOnly) {
     m_traj->self().unset(target, m_istate);
   }
 
   /// Add additional components to the track state
   /// @note Only available if the track state proxy is not read-only
   /// @param mask The bitmask that instructs which components to allocate
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  void addComponents(TrackStatePropMask mask) {
+  void addComponents(TrackStatePropMask mask) requires(!ReadOnly) {
     m_traj->self().addTrackStateComponents_impl(m_istate, mask);
   }
 
@@ -364,8 +361,8 @@ class TrackStateProxy {
   /// Set the reference surface to a given value
   /// @param srf Shared pointer to the surface to set
   /// @note This overload is only present in case @c ReadOnly is false.
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  void setReferenceSurface(std::shared_ptr<const Surface> srf) {
+  void setReferenceSurface(std::shared_ptr<const Surface> srf) requires(
+      !ReadOnly) {
     m_traj->setReferenceSurface(m_istate, std::move(srf));
   }
   // NOLINTEND(performance-unnecessary-value-param)
@@ -375,8 +372,7 @@ class TrackStateProxy {
   /// value directly into the backing store.
   /// @note this overload is only enabled in case the proxy is not read-only
   /// @return Mutable reference to the chi2 value
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  float& chi2() {
+  float& chi2() requires(!ReadOnly) {
     return component<float, hashString("chi2")>();
   }
 
@@ -405,8 +401,7 @@ class TrackStateProxy {
   /// This overloaded is only enabled if not read-only, and returns a mutable
   /// reference.
   /// @return reference to the type flags.
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  TrackStateType typeFlags() {
+  TrackStateType typeFlags() requires(!ReadOnly) {
     return TrackStateType{
         component<TrackStateType::raw_type, hashString("typeFlags")>()};
   }
@@ -460,8 +455,7 @@ class TrackStateProxy {
         component<IndexType, hashString("predicted")>());
   }
 
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  Covariance predictedCovariance() {
+  Covariance predictedCovariance() requires(!ReadOnly) {
     assert(has<hashString("predicted")>());
     return m_traj->self().covariance(
         component<IndexType, hashString("predicted")>());
@@ -483,8 +477,7 @@ class TrackStateProxy {
   /// Filtered track parameters vector
   /// @return The filtered parameters
   /// @note Mutable version
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  Parameters filtered() {
+  Parameters filtered() requires(!ReadOnly) {
     assert(has<hashString("filtered")>());
     return m_traj->self().parameters(
         component<IndexType, hashString("filtered")>());
@@ -502,8 +495,7 @@ class TrackStateProxy {
   /// Filtered track parameters covariance matrix
   /// @return The filtered parameters covariance
   /// @note Mutable version
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  Covariance filteredCovariance() {
+  Covariance filteredCovariance() requires(!ReadOnly) {
     assert(has<hashString("filtered")>());
     return m_traj->self().covariance(
         component<IndexType, hashString("filtered")>());
@@ -525,8 +517,7 @@ class TrackStateProxy {
   /// Smoothed track parameters vector
   /// @return The smoothed parameters
   /// @note Mutable version
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  Parameters smoothed() {
+  Parameters smoothed() requires(!ReadOnly) {
     assert(has<hashString("smoothed")>());
     return m_traj->self().parameters(
         component<IndexType, hashString("smoothed")>());
@@ -544,8 +535,7 @@ class TrackStateProxy {
   /// Smoothed track parameters covariance matrix
   /// @return the parameter covariance matrix
   /// @note Mutable version
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  Covariance smoothedCovariance() {
+  Covariance smoothedCovariance() requires(!ReadOnly) {
     assert(has<hashString("smoothed")>());
     return m_traj->self().covariance(
         component<IndexType, hashString("smoothed")>());
@@ -566,8 +556,7 @@ class TrackStateProxy {
   /// Returns the jacobian from the previous trackstate to this one
   /// @return The jacobian matrix
   /// @note Mutable version
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  Covariance jacobian() {
+  Covariance jacobian() requires(!ReadOnly) {
     assert(has<hashString("jacobian")>());
     return m_traj->self().jacobian(m_istate);
   }
@@ -633,9 +622,9 @@ class TrackStateProxy {
   /// and store it.
   /// @param projector The projector in the form of a dense matrix
   /// @note @p projector is assumed to only have 0s or 1s as components.
-  template <typename Derived, bool RO = ReadOnly,
-            typename = std::enable_if_t<!RO>>
-  void setProjector(const Eigen::MatrixBase<Derived>& projector) {
+  template <typename Derived>
+  void setProjector(const Eigen::MatrixBase<Derived>& projector) requires(
+      !ReadOnly) {
     constexpr int rows = Eigen::MatrixBase<Derived>::RowsAtCompileTime;
     constexpr int cols = Eigen::MatrixBase<Derived>::ColsAtCompileTime;
 
@@ -678,8 +667,7 @@ class TrackStateProxy {
   /// @note This is mainly to copy explicitly a projector from one state
   ///       to another. If you have a projection matrix, set it with
   ///       `setProjector`.
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  void setProjectorBitset(ProjectorBitset proj) {
+  void setProjectorBitset(ProjectorBitset proj) requires(!ReadOnly) {
     assert(has<hashString("projector")>());
     component<ProjectorBitset, hashString("projector")>() = proj;
   }
@@ -690,8 +678,7 @@ class TrackStateProxy {
 
   /// Set an uncalibrated source link
   /// @param sourceLink The uncalibrated source link to set
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  void setUncalibratedSourceLink(SourceLink sourceLink) {
+  void setUncalibratedSourceLink(SourceLink sourceLink) requires(!ReadOnly) {
     m_traj->setUncalibratedSourceLink(m_istate, std::move(sourceLink));
   }
 
@@ -719,9 +706,8 @@ class TrackStateProxy {
   /// dimensions.
   /// @return The measurement vector
   /// @note Mutable version
-  template <std::size_t measdim, bool RO = ReadOnly,
-            typename = std::enable_if_t<!RO>>
-  Calibrated<measdim> calibrated() {
+  template <std::size_t measdim>
+  Calibrated<measdim> calibrated() requires(!ReadOnly) {
     assert(has<hashString("calibrated")>());
     return m_traj->self().template calibrated<measdim>(m_istate);
   }
@@ -738,9 +724,8 @@ class TrackStateProxy {
   /// Mutable full calibrated measurement covariance matrix. The effective
   /// covariance is located in the top left corner, everything else is zeroed.
   /// @return The measurement covariance matrix
-  template <std::size_t measdim, bool RO = ReadOnly,
-            typename = std::enable_if_t<!RO>>
-  CalibratedCovariance<measdim> calibratedCovariance() {
+  template <std::size_t measdim>
+  CalibratedCovariance<measdim> calibratedCovariance() requires(!ReadOnly) {
     assert(has<hashString("calibratedCov")>());
     return m_traj->self().template calibratedCovariance<measdim>(m_istate);
   }
@@ -748,8 +733,7 @@ class TrackStateProxy {
   /// Mutable dynamic measurement vector with only the valid dimensions.
   /// @warning The dynamic vector has a runtime overhead!
   /// @return The effective calibrated measurement vector
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  EffectiveCalibrated effectiveCalibrated() {
+  EffectiveCalibrated effectiveCalibrated() requires(!ReadOnly) {
     assert(has<hashString("calibrated")>());
     return m_traj->self().effectiveCalibrated(m_istate);
   }
@@ -814,9 +798,8 @@ class TrackStateProxy {
   /// @param shareSource Which component to share from
   /// @param shareTarget Which component to share as. This should be different from
   ///                    as @p shareSource, e.g. predicted can be shared as filtered.
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
   void shareFrom(TrackStatePropMask shareSource,
-                 TrackStatePropMask shareTarget) {
+                 TrackStatePropMask shareTarget) requires(!ReadOnly) {
     shareFrom(*this, shareSource, shareTarget);
   }
 
@@ -825,10 +808,9 @@ class TrackStateProxy {
   /// @param component Which component to share.
   /// @note The track states both need to be stored in the
   ///       same @c MultiTrajectory instance
-  template <bool RO = ReadOnly, bool ReadOnlyOther,
-            typename = std::enable_if_t<!RO>>
+  template <bool ReadOnlyOther>
   void shareFrom(const TrackStateProxy<Trajectory, M, ReadOnlyOther>& other,
-                 TrackStatePropMask component) {
+                 TrackStatePropMask component) requires(!ReadOnly) {
     shareFrom(other, component, component);
   }
 
@@ -839,11 +821,10 @@ class TrackStateProxy {
   ///                    as @p shareSource, e.g. predicted can be shared as filtered.
   /// @note Shareable components are predicted, filtered, smoothed, calibrated, jacobian,
   ///       or projector. See @c TrackStatePropMask.
-  template <bool RO = ReadOnly, bool ReadOnlyOther,
-            typename = std::enable_if_t<!RO>>
+  template <bool ReadOnlyOther>
   void shareFrom(const TrackStateProxy<Trajectory, M, ReadOnlyOther>& other,
                  TrackStatePropMask shareSource,
-                 TrackStatePropMask shareTarget) {
+                 TrackStatePropMask shareTarget) requires(!ReadOnly) {
     assert(m_traj == other.m_traj &&
            "Cannot share components across MultiTrajectories");
 
@@ -863,11 +844,10 @@ class TrackStateProxy {
   ///       an exception is thrown.
   /// @note The mask parameter will not cause a copy of components that are
   ///       not allocated in the source track state proxy.
-  template <TrackStateProxyConcept track_state_proxy_t, bool RO = ReadOnly,
-            typename = std::enable_if_t<!RO>>
+  template <TrackStateProxyConcept track_state_proxy_t>
   void copyFrom(const track_state_proxy_t& other,
                 TrackStatePropMask mask = TrackStatePropMask::All,
-                bool onlyAllocated = true) {
+                bool onlyAllocated = true) requires(!ReadOnly) {
     using PM = TrackStatePropMask;
 
     if (onlyAllocated) {
@@ -1022,9 +1002,8 @@ class TrackStateProxy {
   /// @tparam T The type of the component to access
   /// @tparam key String key for the component to access
   /// @return Mutable reference to the component given by @p key
-  template <typename T, HashedString key, bool RO = ReadOnly,
-            typename = std::enable_if_t<!RO>>
-  constexpr T& component() {
+  template <typename T, HashedString key>
+  constexpr T& component() requires(!ReadOnly) {
     return m_traj->template component<T, key>(m_istate);
   }
 
@@ -1032,8 +1011,8 @@ class TrackStateProxy {
   /// @tparam T The type of the component to access
   /// @param key String key for the component to access
   /// @return Mutable reference to the component given by @p key
-  template <typename T, bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  constexpr T& component(HashedString key) {
+  template <typename T>
+  constexpr T& component(HashedString key) requires(!ReadOnly) {
     return m_traj->template component<T>(key, m_istate);
   }
 
@@ -1042,8 +1021,8 @@ class TrackStateProxy {
   /// @param key String key for the component to access
   /// @note This might hash the @p key at runtime instead of compile-time
   /// @return Mutable reference to the component given by @p key
-  template <typename T, bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  constexpr T& component(std::string_view key) {
+  template <typename T>
+  constexpr T& component(std::string_view key) requires(!ReadOnly) {
     return m_traj->template component<T>(hashString(key), m_istate);
   }
 
@@ -1079,8 +1058,7 @@ class TrackStateProxy {
 
   /// Return a mutable reference to the underlying backend container
   /// @return A reference to the backend container
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  MultiTrajectory<Trajectory>& trajectory() {
+  MultiTrajectory<Trajectory>& trajectory() requires(!ReadOnly) {
     return *m_traj;
   }
 
@@ -1090,10 +1068,7 @@ class TrackStateProxy {
 
   /// Get a mutable reference to the track state container backend
   /// @return a mutable reference to the backend
-  template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>>
-  auto& container() {
-    return *m_traj;
-  }
+  auto& container() requires(!ReadOnly) { return *m_traj; }
 
   /// Get a const reference to the track state container backend
   /// @return a const reference to the backend
diff --git a/Core/include/Acts/Navigation/DetectorNavigator.hpp b/Core/include/Acts/Navigation/DetectorNavigator.hpp
index e9fdbd58e0e..02522888f14 100644
--- a/Core/include/Acts/Navigation/DetectorNavigator.hpp
+++ b/Core/include/Acts/Navigation/DetectorNavigator.hpp
@@ -296,7 +296,7 @@ class DetectorNavigator {
                      << posInfo(state, stepper)
                      << "this is a portal, updating to new volume.");
         nState.currentPortal = nextPortal;
-
+        nState.currentSurface = &nextPortal->surface();
         nState.surfaceCandidates.clear();
         nState.surfaceCandidateIndex = 0;
 
diff --git a/Core/include/Acts/Navigation/NavigationStateUpdaters.hpp b/Core/include/Acts/Navigation/NavigationStateUpdaters.hpp
index d1eed8b948b..8bd00d6a406 100644
--- a/Core/include/Acts/Navigation/NavigationStateUpdaters.hpp
+++ b/Core/include/Acts/Navigation/NavigationStateUpdaters.hpp
@@ -56,7 +56,7 @@ inline void updateCandidates(const GeometryContext& gctx,
         gctx, position, direction, c.boundaryTolerance, s_onSurfaceTolerance);
     for (auto& si : sIntersection.split()) {
       c.objectIntersection = si;
-      if (c.objectIntersection &&
+      if (c.objectIntersection.isValid() &&
           c.objectIntersection.pathLength() > overstepTolerance) {
         nextSurfaceCandidates.emplace_back(c);
       }
diff --git a/Core/include/Acts/Propagator/DirectNavigator.hpp b/Core/include/Acts/Propagator/DirectNavigator.hpp
index 20f15ad8221..d389c36e7e0 100644
--- a/Core/include/Acts/Propagator/DirectNavigator.hpp
+++ b/Core/include/Acts/Propagator/DirectNavigator.hpp
@@ -254,7 +254,7 @@ class DirectNavigator {
                    << state.navigation.currentSurface->geometryId());
       // Move the sequence to the next surface
       ++state.navigation.surfaceIndex;
-      if (state.navigation.surfaceIndex >=
+      if (state.navigation.surfaceIndex <
           state.navigation.options.surfaces.size()) {
         ACTS_VERBOSE("Next surface candidate is  "
                      << state.navigation.options.surfaces
@@ -277,8 +277,8 @@ class DirectNavigator {
                                            boundaryTolerance, tolerance);
 
     for (auto& intersection : intersections.split()) {
-      if (detail::checkIntersection(intersection, nearLimit, farLimit,
-                                    logger())) {
+      if (detail::checkPathLength(intersection.pathLength(), nearLimit,
+                                  farLimit, logger())) {
         return intersection;
       }
     }
diff --git a/Core/include/Acts/Propagator/Navigator.hpp b/Core/include/Acts/Propagator/Navigator.hpp
index ab390fa4cb1..9d8efc93428 100644
--- a/Core/include/Acts/Propagator/Navigator.hpp
+++ b/Core/include/Acts/Propagator/Navigator.hpp
@@ -989,7 +989,7 @@ class Navigator {
                   state.options.direction * stepper.direction(state.stepping),
                   BoundaryTolerance::Infinite(), state.options.surfaceTolerance)
               .closest();
-      if (targetIntersection) {
+      if (targetIntersection.isValid()) {
         ACTS_VERBOSE(volInfo(state)
                      << "Target estimate position ("
                      << targetIntersection.position().x() << ", "
diff --git a/Core/include/Acts/Propagator/StandardAborters.hpp b/Core/include/Acts/Propagator/StandardAborters.hpp
index 144edfda1df..23c045eac7e 100644
--- a/Core/include/Acts/Propagator/StandardAborters.hpp
+++ b/Core/include/Acts/Propagator/StandardAborters.hpp
@@ -130,9 +130,9 @@ struct SurfaceReached {
     bool intersectionFound = false;
 
     for (const auto& intersection : sIntersection.split()) {
-      if (intersection &&
-          detail::checkIntersection(intersection.intersection(), nearLimit,
-                                    farLimit, logger)) {
+      if (intersection.isValid() &&
+          detail::checkPathLength(intersection.pathLength(), nearLimit,
+                                  farLimit, logger)) {
         stepper.updateStepSize(state.stepping, intersection.pathLength(),
                                ConstrainedStep::aborter, false);
         ACTS_VERBOSE(
diff --git a/Core/include/Acts/Propagator/TryAllNavigator.hpp b/Core/include/Acts/Propagator/TryAllNavigator.hpp
index 744e393ab29..a516a3292d2 100644
--- a/Core/include/Acts/Propagator/TryAllNavigator.hpp
+++ b/Core/include/Acts/Propagator/TryAllNavigator.hpp
@@ -19,6 +19,7 @@
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/Intersection.hpp"
 #include "Acts/Utilities/Logger.hpp"
+#include "Acts/Utilities/StringHelpers.hpp"
 
 #include <algorithm>
 #include <cstdint>
@@ -342,7 +343,9 @@ class TryAllNavigator : public TryAllNavigatorBase {
                               state.options.surfaceTolerance);
       for (const auto& intersection : intersections.first.split()) {
         // exclude invalid intersections
-        if (!detail::checkIntersection(intersection, nearLimit, farLimit)) {
+        if (!intersection.isValid() ||
+            !detail::checkPathLength(intersection.pathLength(), nearLimit,
+                                     farLimit)) {
           continue;
         }
         // store candidate
@@ -357,6 +360,8 @@ class TryAllNavigator : public TryAllNavigatorBase {
     ACTS_VERBOSE(volInfo(state) << "found " << intersectionCandidates.size()
                                 << " intersections");
 
+    bool intersectionFound = false;
+
     for (const auto& candidate : intersectionCandidates) {
       const auto& intersection = candidate.intersection;
       const Surface& surface = *intersection.object();
@@ -375,13 +380,16 @@ class TryAllNavigator : public TryAllNavigatorBase {
         continue;
       }
 
-      ACTS_VERBOSE(volInfo(state) << "aiming at surface "
-                                  << surface.geometryId() << ". step size is "
-                                  << stepper.outputStepSize(state.stepping));
-      break;
+      if (surfaceStatus == IntersectionStatus::reachable) {
+        ACTS_VERBOSE(volInfo(state)
+                     << "Surface reachable, step size updated to "
+                     << stepper.outputStepSize(state.stepping));
+        intersectionFound = true;
+        break;
+      }
     }
 
-    if (intersectionCandidates.empty()) {
+    if (!intersectionFound) {
       stepper.releaseStepSize(state.stepping, ConstrainedStep::actor);
 
       ACTS_VERBOSE(volInfo(state) << "no intersections found. advance without "
@@ -675,7 +683,8 @@ class TryAllOverstepNavigator : public TryAllNavigatorBase {
         break;
       }
 
-      ACTS_WARNING(volInfo(state) << "Surface unreachable, skip.");
+      ACTS_VERBOSE(volInfo(state) << "Surface " << surface.geometryId()
+                                  << " unreachable, skip.");
       ++state.navigation.activeCandidateIndex;
     }
 
@@ -738,8 +747,9 @@ class TryAllOverstepNavigator : public TryAllNavigatorBase {
             state.geoContext, end, direction, state.options.surfaceTolerance);
         for (const auto& intersection : intersections.first.split()) {
           // exclude invalid intersections
-          if (!intersection ||
-              !detail::checkIntersection(intersection, nearLimit, farLimit)) {
+          if (!intersection.isValid() ||
+              !detail::checkPathLength(intersection.pathLength(), nearLimit,
+                                       farLimit)) {
             continue;
           }
           // exclude last candidate
diff --git a/Core/include/Acts/Propagator/detail/SteppingHelper.hpp b/Core/include/Acts/Propagator/detail/SteppingHelper.hpp
index 1ccd139fdb3..01db7f6c8fd 100644
--- a/Core/include/Acts/Propagator/detail/SteppingHelper.hpp
+++ b/Core/include/Acts/Propagator/detail/SteppingHelper.hpp
@@ -56,8 +56,9 @@ Acts::Intersection3D::Status updateSingleSurfaceStatus(
   const double nearLimit = std::numeric_limits<double>::lowest();
   const double farLimit = state.stepSize.value(ConstrainedStep::aborter);
 
-  if (sIntersection && detail::checkIntersection(sIntersection.intersection(),
-                                                 nearLimit, farLimit, logger)) {
+  if (sIntersection.isValid() &&
+      detail::checkPathLength(sIntersection.pathLength(), nearLimit, farLimit,
+                              logger)) {
     ACTS_VERBOSE("Surface is reachable");
     stepper.updateStepSize(state, sIntersection.pathLength(),
                            ConstrainedStep::actor);
diff --git a/Core/include/Acts/Seeding/CandidatesForMiddleSp.hpp b/Core/include/Acts/Seeding/CandidatesForMiddleSp.hpp
index be50a0546ed..9affac3af92 100644
--- a/Core/include/Acts/Seeding/CandidatesForMiddleSp.hpp
+++ b/Core/include/Acts/Seeding/CandidatesForMiddleSp.hpp
@@ -8,10 +8,13 @@
 
 #pragma once
 
+#include "Acts/Definitions/Algebra.hpp"
+
 #include <algorithm>
 #include <limits>
 #include <memory>
 #include <tuple>
+#include <utility>
 #include <vector>
 
 namespace Acts {
@@ -41,8 +44,23 @@ struct TripletCandidate {
   TripletCandidate& operator=(const TripletCandidate&) = default;
 
   /// @brief Move operations
-  TripletCandidate(TripletCandidate&&) = default;
-  TripletCandidate& operator=(TripletCandidate&&) = default;
+  TripletCandidate(TripletCandidate&& other) noexcept
+      : bottom(std::exchange(other.bottom, nullptr)),
+        middle(std::exchange(other.middle, nullptr)),
+        top(std::exchange(other.top, nullptr)),
+        weight(other.weight),
+        zOrigin(other.zOrigin),
+        isQuality(other.isQuality) {}
+
+  TripletCandidate& operator=(TripletCandidate&& other) noexcept {
+    bottom = std::exchange(other.bottom, nullptr);
+    middle = std::exchange(other.middle, nullptr);
+    top = std::exchange(other.top, nullptr);
+    weight = other.weight;
+    zOrigin = other.zOrigin;
+    isQuality = other.isQuality;
+    return *this;
+  }
 
   external_space_point_t* bottom{nullptr};
   external_space_point_t* middle{nullptr};
@@ -60,19 +78,21 @@ struct TripletCandidate {
 /// @tparam external_space_point_t The external spacepoint type.
 
 template <typename external_space_point_t>
+concept SatisfyCandidateConcept = requires(external_space_point_t spacePoint) {
+  { spacePoint.x() } -> std::convertible_to<float>;
+  { spacePoint.y() } -> std::convertible_to<float>;
+  { spacePoint.z() } -> std::convertible_to<float>;
+};
+
+template <SatisfyCandidateConcept external_space_point_t>
 class CandidatesForMiddleSp {
  public:
   using value_type = TripletCandidate<external_space_point_t>;
 
-  /// @brief constructor
-  CandidatesForMiddleSp() = default;
-  /// @brief Destructor
-  ~CandidatesForMiddleSp() = default;
-
   /// @brief Setting maximum number of candidates to keep
-  /// @param n_low Maximum number of candidates in the low-quality collection
-  /// @param n_high Maximum number of candidates in the high-quality collection
-  void setMaxElements(std::size_t n_low, std::size_t n_high);
+  /// @param nLow Maximum number of candidates in the low-quality collection
+  /// @param nHigh Maximum number of candidates in the high-quality collection
+  void setMaxElements(std::size_t nLow, std::size_t nHigh);
 
   /// @brief Retrieve the triplet candidates, the resulting vector is already sorted,
   /// elements with higher quality first
@@ -107,12 +127,20 @@ class CandidatesForMiddleSp {
   /// @returns The comparison result
   static bool ascendingByQuality(const value_type& i1, const value_type& i2);
 
+  /// @brief Retrieve the number of Low quality candidates
+  /// @returns The number of Low quality candidates
+  std::size_t nLowQualityCandidates() const;
+
+  /// @brief Retrieve the number of High quality candidates
+  /// @returns The number of High quality candidates
+  std::size_t nHighQualityCandidates() const;
+
  private:
   /// @brief dding a new triplet candidate to the collection, should it satisfy the
   /// selection criteria
   /// @param indices The collection into which the candidate should be stored
   /// @param n The current number of stored elements in the container
-  /// @param n_max The maximum number of elements that can be stored in the container
+  /// @param nMax The maximum number of elements that can be stored in the container
   /// @param SpB The bottom space point
   /// @param SpM The middle space point
   /// @param SpT The top space point
@@ -121,25 +149,25 @@ class CandidatesForMiddleSp {
   /// @param isQuality Whether the triplet candidate is high or low quality
   /// @returns whether the triplet candidate has been added or not to the collection
   bool push(std::vector<std::size_t>& indices, std::size_t& n,
-            const std::size_t n_max, external_space_point_t& SpB,
+            const std::size_t nMax, external_space_point_t& SpB,
             external_space_point_t& SpM, external_space_point_t& SpT,
             float weight, float zOrigin, bool isQuality);
 
   /// @brief Check if an element exists in the collection. The element to be checked
   /// is supposed to be in the n position of the collection.
   /// @param n Index of the requested element
-  /// @param max_size Number of elements currently stored in the collection
+  /// @param maxSize Number of elements currently stored in the collection
   /// @returns Whether the element exists
-  bool exists(const std::size_t n, const std::size_t max_size) const;
+  bool exists(const std::size_t n, const std::size_t maxSize) const;
 
   /// @brief Pop an element from a collection. The removal of the element from the collection
   /// does not imply its destruction. In fact, the number of stored elements is
   /// simply diminished by 1. The popped element is technically still available
   /// at the end of the collection.
   /// @param indices The collection
-  /// @param current_size The current number of element stored in the collection. The function will
+  /// @param currentSize The current number of element stored in the collection. The function will
   /// diminish this value by 1
-  void pop(std::vector<std::size_t>& indices, std::size_t& current_size);
+  void pop(std::vector<std::size_t>& indices, std::size_t& currentSize);
 
   /// @brief Return the weight for a candidate
   /// @param indices The collection in which the element is stored
@@ -160,29 +188,29 @@ class CandidatesForMiddleSp {
   /// is in the correct position on the tree
   /// @param indices The collection
   /// @param n The index of the element to place in the correct position
-  /// @param actual_size The current number of elements stored in the collection
+  /// @param actualSize The current number of elements stored in the collection
   void bubbledw(std::vector<std::size_t>& indices, std::size_t n,
-                std::size_t actual_size);
+                std::size_t actualSize);
 
   /// @brief Adding a new triplet candidate to the collection. The function is called after the candidate has satisfied
   /// all the selection criteria
   /// @param indices The collection
   /// @param n Current number of stored elements in the collection
-  /// @param n_max The maximum number of elements that can be stored in the collection
+  /// @param nMax The maximum number of elements that can be stored in the collection
   /// @param element The element that must be added to the collection
   void addToCollection(std::vector<std::size_t>& indices, std::size_t& n,
-                       const std::size_t n_max, value_type&& element);
+                       const std::size_t nMax, value_type&& element);
 
  private:
   // sizes
-  // m_max_size_* is the maximum size of the indices collections. These values
+  // m_maxSize* is the maximum size of the indices collections. These values
   // are set by the user once
-  std::size_t m_max_size_high{0};
-  std::size_t m_max_size_low{0};
-  // m_n_* is the current size of the indices collections [0, m_max_size_*).
+  std::size_t m_maxSizeHigh{0};
+  std::size_t m_maxSizeLow{0};
+  // m_n_* is the current size of the indices collections [0, m_maxSize*).
   // These values are set internally by the class
-  std::size_t m_n_high{0};
-  std::size_t m_n_low{0};
+  std::size_t m_nHigh{0};
+  std::size_t m_nLow{0};
 
   // storage contains the collection of the candidates
   std::vector<value_type> m_storage{};
@@ -197,9 +225,9 @@ class CandidatesForMiddleSp {
   // and std::priority_queue were tried and were found to be slower.
 
   // list of indexes of candidates with high quality in the storage
-  std::vector<std::size_t> m_indices_high{};
+  std::vector<std::size_t> m_indicesHigh{};
   // list of indexes of candidates with low quality in the storage
-  std::vector<std::size_t> m_indices_low{};
+  std::vector<std::size_t> m_indicesLow{};
 };
 
 }  // namespace Acts
diff --git a/Core/include/Acts/Seeding/CandidatesForMiddleSp.ipp b/Core/include/Acts/Seeding/CandidatesForMiddleSp.ipp
index 883e827868c..7fbaa006a59 100644
--- a/Core/include/Acts/Seeding/CandidatesForMiddleSp.ipp
+++ b/Core/include/Acts/Seeding/CandidatesForMiddleSp.ipp
@@ -8,117 +8,116 @@
 
 namespace Acts {
 
-template <typename external_space_point_t>
+template <SatisfyCandidateConcept external_space_point_t>
 inline void CandidatesForMiddleSp<external_space_point_t>::setMaxElements(
-    std::size_t n_low, std::size_t n_high) {
-  m_max_size_high = n_high;
-  m_max_size_low = n_low;
+    std::size_t nLow, std::size_t nHigh) {
+  m_maxSizeHigh = nHigh;
+  m_maxSizeLow = nLow;
 
   // protection against default numbers
   // it may cause std::bad_alloc if we don't protect
-  if (n_high == std::numeric_limits<std::size_t>::max() ||
-      n_low == std::numeric_limits<std::size_t>::max()) {
+  if (nHigh == std::numeric_limits<std::size_t>::max() ||
+      nLow == std::numeric_limits<std::size_t>::max()) {
     return;
   }
 
   // Reserve enough memory for all collections
-  m_storage.reserve(n_low + n_high);
-  m_indices_high.reserve(n_high);
-  m_indices_low.reserve(n_low);
+  m_storage.reserve(nLow + nHigh);
+  m_indicesHigh.reserve(nHigh);
+  m_indicesLow.reserve(nLow);
 }
 
-template <typename external_space_point_t>
+template <SatisfyCandidateConcept external_space_point_t>
 inline void CandidatesForMiddleSp<external_space_point_t>::pop(
-    std::vector<std::size_t>& indices, std::size_t& current_size) {
+    std::vector<std::size_t>& indices, std::size_t& currentSize) {
   // Remove the candidate with the lowest weight in the collection
   // By construction, this element is always the first element in the
   // collection.
   // The removal works this way: the first element of the collection
   // is overwritten with the last element of the collection.
-  // The number of stored elements (current_size) is lowered by 1
+  // The number of stored elements (currentSize) is lowered by 1
   // The new first element is moved down in the tree to its correct position
-  std::swap(indices[0], indices[current_size - 1]);
-  bubbledw(indices, 0, --current_size);
+  std::swap(indices[0], indices[currentSize - 1]);
+  bubbledw(indices, 0, --currentSize);
 }
 
-template <typename external_space_point_t>
+template <SatisfyCandidateConcept external_space_point_t>
 inline bool CandidatesForMiddleSp<external_space_point_t>::exists(
-    const std::size_t n, const std::size_t max_size) const {
+    const std::size_t n, const std::size_t maxSize) const {
   // If the element exists, its index is lower than the current number
   // of stored elements
-  return n < max_size;
+  return n < maxSize;
 }
 
-template <typename external_space_point_t>
+template <SatisfyCandidateConcept external_space_point_t>
 inline float CandidatesForMiddleSp<external_space_point_t>::weight(
     const std::vector<std::size_t>& indices, std::size_t n) const {
   // Get the weight of the n-th element
   return m_storage[indices[n]].weight;
 }
 
-template <typename external_space_point_t>
+template <SatisfyCandidateConcept external_space_point_t>
 inline void CandidatesForMiddleSp<external_space_point_t>::clear() {
   // do not clear max size, this is set only once
   // reset to 0 the number of stored elements
-  m_n_high = 0;
-  m_n_low = 0;
+  m_nHigh = 0;
+  m_nLow = 0;
   // Reset all values to default
   m_storage.clear();
-  m_indices_high.clear();
-  m_indices_low.clear();
+  m_indicesHigh.clear();
+  m_indicesLow.clear();
 }
 
-template <typename external_space_point_t>
+template <SatisfyCandidateConcept external_space_point_t>
 bool CandidatesForMiddleSp<external_space_point_t>::push(
     external_space_point_t& SpB, external_space_point_t& SpM,
     external_space_point_t& SpT, float weight, float zOrigin, bool isQuality) {
   // Decide in which collection this candidate may be added to according to the
   // isQuality boolean
   if (isQuality) {
-    return push(m_indices_high, m_n_high, m_max_size_high, SpB, SpM, SpT,
-                weight, zOrigin, isQuality);
+    return push(m_indicesHigh, m_nHigh, m_maxSizeHigh, SpB, SpM, SpT, weight,
+                zOrigin, isQuality);
   }
-  return push(m_indices_low, m_n_low, m_max_size_low, SpB, SpM, SpT, weight,
+  return push(m_indicesLow, m_nLow, m_maxSizeLow, SpB, SpM, SpT, weight,
               zOrigin, isQuality);
 }
 
-template <typename external_space_point_t>
+template <SatisfyCandidateConcept external_space_point_t>
 bool CandidatesForMiddleSp<external_space_point_t>::push(
-    std::vector<std::size_t>& indices, std::size_t& n, const std::size_t n_max,
+    std::vector<std::size_t>& indices, std::size_t& n, const std::size_t nMax,
     external_space_point_t& SpB, external_space_point_t& SpM,
     external_space_point_t& SpT, float weight, float zOrigin, bool isQuality) {
   // If we do not want to store candidates, returns
-  if (n_max == 0) {
+  if (nMax == 0) {
     return false;
   }
 
   // if there is still space, add anything
-  if (n < n_max) {
-    addToCollection(indices, n, n_max,
+  if (n < nMax) {
+    addToCollection(indices, n, nMax,
                     value_type(SpB, SpM, SpT, weight, zOrigin, isQuality));
     return true;
   }
 
   // if no space, replace one if quality is enough
   // compare to element with lowest weight
-  const auto& lowest_weight = this->weight(indices, 0);
-  if (weight <= lowest_weight) {
+  if (float lowestWeight = this->weight(indices, 0); weight <= lowestWeight) {
     return false;
   }
 
   // remove element with lower weight and add this one
   pop(indices, n);
-  addToCollection(indices, n, n_max,
+  addToCollection(indices, n, nMax,
                   value_type(SpB, SpM, SpT, weight, zOrigin, isQuality));
   return true;
 }
 
-template <typename external_space_point_t>
+template <SatisfyCandidateConcept external_space_point_t>
 void CandidatesForMiddleSp<external_space_point_t>::addToCollection(
-    std::vector<std::size_t>& indices, std::size_t& n, const std::size_t n_max,
+    std::vector<std::size_t>& indices, std::size_t& n, const std::size_t nMax,
     value_type&& element) {
   // adds elements to the end of the collection
-  if (indices.size() == n_max) {
+  if (indices.size() == nMax) {
     m_storage[indices[n]] = std::move(element);
   } else {
     m_storage.push_back(std::move(element));
@@ -128,16 +127,16 @@ void CandidatesForMiddleSp<external_space_point_t>::addToCollection(
   bubbleup(indices, n++);
 }
 
-template <typename external_space_point_t>
+template <SatisfyCandidateConcept external_space_point_t>
 void CandidatesForMiddleSp<external_space_point_t>::bubbledw(
-    std::vector<std::size_t>& indices, std::size_t n, std::size_t actual_size) {
-  while (n < actual_size) {
+    std::vector<std::size_t>& indices, std::size_t n, std::size_t actualSize) {
+  while (n < actualSize) {
     // The collection of indexes are sorted as min heap trees
     // left child : 2 * n + 1
     // right child: 2 * n + 2
     float current = weight(indices, n);
-    std::size_t left_child = 2 * n + 1;
-    std::size_t right_child = 2 * n + 2;
+    std::size_t leftChild = 2 * n + 1;
+    std::size_t rightChild = 2 * n + 2;
 
     // We have to move the current node down the tree to its correct position.
     // This is done by comparing its weight with the weights of its two
@@ -150,7 +149,7 @@ void CandidatesForMiddleSp<external_space_point_t>::bubbledw(
 
     // if there is no left child, that also means no right child is present.
     // We do nothing
-    if (!exists(left_child, actual_size)) {
+    if (!exists(leftChild, actualSize)) {
       break;
     }
 
@@ -158,96 +157,95 @@ void CandidatesForMiddleSp<external_space_point_t>::bubbledw(
     // have to check. We take the lowest weight of the children. By default this
     // is the weight of the left child, and we then check for the right child
 
-    float weight_left_child = weight(indices, left_child);
+    float weightLeftChild = weight(indices, leftChild);
 
-    std::size_t selected_child = left_child;
-    float selected_weight = weight_left_child;
+    std::size_t selectedChild = leftChild;
+    float selectedWeight = weightLeftChild;
 
     // Check which child has the lower weight
-    if (exists(right_child, actual_size)) {
-      float weight_right_child = weight(indices, right_child);
-      if (weight_right_child <= weight_left_child) {
-        selected_child = right_child;
-        selected_weight = weight_right_child;
+    if (exists(rightChild, actualSize)) {
+      float weightRightChild = weight(indices, rightChild);
+      if (weightRightChild <= weightLeftChild) {
+        selectedChild = rightChild;
+        selectedWeight = weightRightChild;
       }
     }
 
     // At this point we have the minimum weight of the children
     // We can compare this to the current weight
     // If weight of the children is higher we stop
-    if (selected_weight >= current) {
+    if (selectedWeight >= current) {
       break;
     }
 
     // swap and repeat the process
-    std::swap(indices[n], indices[selected_child]);
-    n = selected_child;
+    std::swap(indices[n], indices[selectedChild]);
+    n = selectedChild;
   }  // while loop
 }
 
-template <typename external_space_point_t>
+template <SatisfyCandidateConcept external_space_point_t>
 void CandidatesForMiddleSp<external_space_point_t>::bubbleup(
     std::vector<std::size_t>& indices, std::size_t n) {
   while (n != 0) {
     // The collection of indexes are sorted as min heap trees
     // parent: (n - 1) / 2;
     // this works because it is an integer operation
-    std::size_t parent_idx = (n - 1) / 2;
-
-    float weight_current = weight(indices, n);
-    float weight_parent = weight(indices, parent_idx);
+    std::size_t parentIdx = (n - 1) / 2;
 
+    float weightCurrent = weight(indices, n);
     // If weight of the parent is lower than this one, we stop
-    if (weight_parent <= weight_current) {
+    if (float weightParent = weight(indices, parentIdx);
+        weightParent <= weightCurrent) {
       break;
     }
 
     // swap and repeat the process
-    std::swap(indices[n], indices[parent_idx]);
-    n = parent_idx;
+    std::swap(indices[n], indices[parentIdx]);
+    n = parentIdx;
   }
 }
 
-template <typename external_space_point_t>
+template <SatisfyCandidateConcept external_space_point_t>
 std::vector<typename CandidatesForMiddleSp<external_space_point_t>::value_type>
 CandidatesForMiddleSp<external_space_point_t>::storage() {
   // this will retrieve the entire storage
   // the resulting vector is already sorted from high to low quality
-  std::vector<value_type> output(m_n_high + m_n_low);
-  std::size_t out_idx = output.size() - 1;
+  std::vector<value_type> output(m_nHigh + m_nLow);
+  std::size_t outIdx = output.size() - 1;
 
-  // rely on the fact that m_indices_* are both min heap trees
+  // rely on the fact that m_indices* are both min heap trees
   // Sorting comes naturally by popping elements one by one and
   // placing this element at the end of the output vector
-  while (m_n_high != 0 || m_n_low != 0) {
+  while (m_nHigh != 0 || m_nLow != 0) {
     // no entries in collection high, we attach the entire low collection
-    if (m_n_high == 0) {
-      std::size_t idx = m_n_low;
+    if (m_nHigh == 0) {
+      std::size_t idx = m_nLow;
       for (std::size_t i(0); i < idx; i++) {
-        output[out_idx--] = std::move(m_storage[m_indices_low[0]]);
-        pop(m_indices_low, m_n_low);
+        output[outIdx--] = std::move(m_storage[m_indicesLow[0]]);
+        pop(m_indicesLow, m_nLow);
       }
       break;
     }
 
     // no entries in collection low, we attach the entire high collection
-    if (m_n_low == 0) {
-      std::size_t idx = m_n_high;
+    if (m_nLow == 0) {
+      std::size_t idx = m_nHigh;
       for (std::size_t i(0); i < idx; i++) {
-        output[out_idx--] = std::move(m_storage[m_indices_high[0]]);
-        pop(m_indices_high, m_n_high);
+        output[outIdx--] = std::move(m_storage[m_indicesHigh[0]]);
+        pop(m_indicesHigh, m_nHigh);
       }
       break;
     }
 
     // Both have entries, get the minimum
-    if (descendingByQuality(m_storage[m_indices_low[0]],
-                            m_storage[m_indices_high[0]])) {
-      output[out_idx--] = std::move(m_storage[m_indices_high[0]]);
-      pop(m_indices_high, m_n_high);
+    if (descendingByQuality(m_storage[m_indicesLow[0]],
+                            m_storage[m_indicesHigh[0]])) {
+      output[outIdx--] = std::move(m_storage[m_indicesHigh[0]]);
+      pop(m_indicesHigh, m_nHigh);
     } else {
-      output[out_idx--] = std::move(m_storage[m_indices_low[0]]);
-      pop(m_indices_low, m_n_low);
+      output[outIdx--] = std::move(m_storage[m_indicesLow[0]]);
+      pop(m_indicesLow, m_nLow);
     }
 
   }  // while loop
@@ -256,7 +254,7 @@ CandidatesForMiddleSp<external_space_point_t>::storage() {
   return output;
 }
 
-template <typename external_space_point_t>
+template <SatisfyCandidateConcept external_space_point_t>
 bool CandidatesForMiddleSp<external_space_point_t>::descendingByQuality(
     const value_type& i1, const value_type& i2) {
   if (i1.weight != i2.weight) {
@@ -266,36 +264,44 @@ bool CandidatesForMiddleSp<external_space_point_t>::descendingByQuality(
   // This is for the case when the weights from different seeds
   // are same. This makes cpu & cuda results same
 
-  const auto& bottom_l1 = i1.bottom;
-  const auto& middle_l1 = i1.middle;
-  const auto& top_l1 = i1.top;
+  const auto& bottomL1 = i1.bottom;
+  const auto& middleL1 = i1.middle;
+  const auto& topL1 = i1.top;
 
-  const auto& bottom_l2 = i2.bottom;
-  const auto& middle_l2 = i2.middle;
-  const auto& top_l2 = i2.top;
+  const auto& bottomL2 = i2.bottom;
+  const auto& middleL2 = i2.middle;
+  const auto& topL2 = i2.top;
 
   float seed1_sum = 0.;
   float seed2_sum = 0.;
 
-  seed1_sum +=
-      bottom_l1->y() * bottom_l1->y() + bottom_l1->z() * bottom_l1->z();
-  seed1_sum +=
-      middle_l1->y() * middle_l1->y() + middle_l1->z() * middle_l1->z();
-  seed1_sum += top_l1->y() * top_l1->y() + top_l1->z() * top_l1->z();
+  seed1_sum += bottomL1->y() * bottomL1->y() + bottomL1->z() * bottomL1->z();
+  seed1_sum += middleL1->y() * middleL1->y() + middleL1->z() * middleL1->z();
+  seed1_sum += topL1->y() * topL1->y() + topL1->z() * topL1->z();
 
-  seed2_sum +=
-      bottom_l2->y() * bottom_l2->y() + bottom_l2->z() * bottom_l2->z();
-  seed2_sum +=
-      middle_l2->y() * middle_l2->y() + middle_l2->z() * middle_l2->z();
-  seed2_sum += top_l2->y() * top_l2->y() + top_l2->z() * top_l2->z();
+  seed2_sum += bottomL2->y() * bottomL2->y() + bottomL2->z() * bottomL2->z();
+  seed2_sum += middleL2->y() * middleL2->y() + middleL2->z() * middleL2->z();
+  seed2_sum += topL2->y() * topL2->y() + topL2->z() * topL2->z();
 
   return seed1_sum > seed2_sum;
 }
 
-template <typename external_space_point_t>
+template <SatisfyCandidateConcept external_space_point_t>
 bool CandidatesForMiddleSp<external_space_point_t>::ascendingByQuality(
     const value_type& i1, const value_type& i2) {
   return !descendingByQuality(i1, i2);
 }
 
+template <SatisfyCandidateConcept external_space_point_t>
+std::size_t
+CandidatesForMiddleSp<external_space_point_t>::nLowQualityCandidates() const {
+  return m_nLow;
+}
+
+template <SatisfyCandidateConcept external_space_point_t>
+std::size_t
+CandidatesForMiddleSp<external_space_point_t>::nHighQualityCandidates() const {
+  return m_nHigh;
+}
+
 }  // namespace Acts
diff --git a/Core/include/Acts/Seeding/SeedFilter.hpp b/Core/include/Acts/Seeding/SeedFilter.hpp
index c1b8ab3542d..95a7b81aa0e 100644
--- a/Core/include/Acts/Seeding/SeedFilter.hpp
+++ b/Core/include/Acts/Seeding/SeedFilter.hpp
@@ -31,12 +31,6 @@ struct SeedFilterState {
   // compatible top required
   float rMaxSeedConf =
       std::numeric_limits<float>::max();  // Acts::UnitConstants::mm
-  // number of high quality seeds in seed confirmation
-  std::size_t numQualitySeeds = 0;
-  // number of seeds that did not pass the quality confirmation but were still
-  // accepted, if quality confirmation is not used this is the total number of
-  // seeds
-  std::size_t numSeeds = 0;
 };
 
 /// Filter seeds at various stages with the currently
@@ -76,7 +70,6 @@ class SeedFilter final {
   /// Filter seeds once all seeds for one middle space point have been created
   /// @param spacePointData Auxiliary variables used by the seeding
   /// @param candidates_collector collection of seed candidates
-  /// @param numQualitySeeds number of high quality seeds in seed confirmation
   /// @param outputCollection Output container for the seeds
   /// for all seeds with the same middle space point
   template <typename collection_t>
@@ -84,7 +77,7 @@ class SeedFilter final {
       Acts::SpacePointData& spacePointData,
       CandidatesForMiddleSp<const InternalSpacePoint<external_spacepoint_t>>&
           candidates_collector,
-      const std::size_t numQualitySeeds, collection_t& outputCollection) const;
+      collection_t& outputCollection) const;
 
   /// Filter seeds once all seeds for one middle space point have been created
   /// @param spacePointData Auxiliary variables used by the seeding
diff --git a/Core/include/Acts/Seeding/SeedFilter.ipp b/Core/include/Acts/Seeding/SeedFilter.ipp
index a282637654e..71d2a6250ce 100644
--- a/Core/include/Acts/Seeding/SeedFilter.ipp
+++ b/Core/include/Acts/Seeding/SeedFilter.ipp
@@ -172,7 +172,8 @@ void SeedFilter<external_spacepoint_t>::filterSeeds_2SpFixed(
       int deltaSeedConf =
           compatibleSeedR.size() + 1 - seedFilterState.nTopSeedConf;
       if (deltaSeedConf < 0 ||
-          (seedFilterState.numQualitySeeds != 0 && deltaSeedConf == 0)) {
+          (candidates_collector.nHighQualityCandidates() != 0 &&
+           deltaSeedConf == 0)) {
         continue;
       }
       bool seedRangeCuts =
@@ -204,11 +205,6 @@ void SeedFilter<external_spacepoint_t>::filterSeeds_2SpFixed(
         // If this is reached, we remove the seed with the lowest weight.
         candidates_collector.push(bottomSP, middleSP, *topSpVec[topSPIndex],
                                   weight, zOrigin, true);
-        if (seedFilterState.numQualitySeeds < m_cfg.maxQualitySeedsPerSpMConf) {
-          // fill high quality seed
-          seedFilterState.numQualitySeeds++;
-        }
-
       } else if (weight > weightMax) {
         // store weight and index of the best "lower quality" seed
         weightMax = weight;
@@ -222,25 +218,17 @@ void SeedFilter<external_spacepoint_t>::filterSeeds_2SpFixed(
 
       candidates_collector.push(bottomSP, middleSP, *topSpVec[topSPIndex],
                                 weight, zOrigin, false);
-      if (seedFilterState.numSeeds < m_cfg.maxSeedsPerSpMConf) {
-        // fill seed
-        seedFilterState.numSeeds++;
-      }
     }
   }  // loop on tops
   // if no high quality seed was found for a certain middle+bottom SP pair,
   // lower quality seeds can be accepted
   if (m_cfg.seedConfirmation && maxWeightSeed &&
-      seedFilterState.numQualitySeeds == 0) {
+      candidates_collector.nHighQualityCandidates() == 0) {
     // if we have not yet reached our max number of seeds we add the new seed to
     // outCont
 
     candidates_collector.push(bottomSP, middleSP, *topSpVec[maxWeightSeedIndex],
                               weightMax, zOrigin, false);
-    if (seedFilterState.numSeeds < m_cfg.maxSeedsPerSpMConf) {
-      // fill seed
-      seedFilterState.numSeeds++;
-    }
   }
 }
 
@@ -252,10 +240,11 @@ void SeedFilter<external_spacepoint_t>::filterSeeds_1SpFixed(
     Acts::SpacePointData& spacePointData,
     CandidatesForMiddleSp<const InternalSpacePoint<external_spacepoint_t>>&
         candidates_collector,
-    const std::size_t numQualitySeeds, collection_t& outputCollection) const {
+    collection_t& outputCollection) const {
   // retrieve all candidates
   // this collection is already sorted
   // higher weights first
+  std::size_t numQualitySeeds = candidates_collector.nHighQualityCandidates();
   auto extended_collection = candidates_collector.storage();
   filterSeeds_1SpFixed(spacePointData, extended_collection, numQualitySeeds,
                        outputCollection);
diff --git a/Core/include/Acts/Seeding/SeedFinder.ipp b/Core/include/Acts/Seeding/SeedFinder.ipp
index f699820e91a..ba36b260909 100644
--- a/Core/include/Acts/Seeding/SeedFinder.ipp
+++ b/Core/include/Acts/Seeding/SeedFinder.ipp
@@ -199,8 +199,7 @@ void SeedFinder<external_spacepoint_t, grid_t, platform_t>::createSeedsForGroup(
     }
 
     m_config.seedFilter->filterSeeds_1SpFixed(
-        state.spacePointData, state.candidates_collector,
-        seedFilterState.numQualitySeeds, outputCollection);
+        state.spacePointData, state.candidates_collector, outputCollection);
 
   }  // loop on mediums
 }
@@ -572,7 +571,8 @@ SeedFinder<external_spacepoint_t, grid_t, platform_t>::filterCandidates(
         state.compatBottomSP[b]->radius() > seedFilterState.rMaxSeedConf) {
       minCompatibleTopSPs = 1;
     }
-    if (m_config.seedConfirmation && seedFilterState.numQualitySeeds) {
+    if (m_config.seedConfirmation &&
+        state.candidates_collector.nHighQualityCandidates()) {
       minCompatibleTopSPs++;
     }
 
diff --git a/Core/include/Acts/Seeding/SeedFinderOrthogonal.ipp b/Core/include/Acts/Seeding/SeedFinderOrthogonal.ipp
index 7b4663dac72..31757e3e72a 100644
--- a/Core/include/Acts/Seeding/SeedFinderOrthogonal.ipp
+++ b/Core/include/Acts/Seeding/SeedFinderOrthogonal.ipp
@@ -363,7 +363,8 @@ void SeedFinderOrthogonal<external_spacepoint_t>::filterCandidates(
         bottom[b]->radius() > seedFilterState.rMaxSeedConf) {
       minCompatibleTopSPs = 1;
     }
-    if (m_config.seedConfirmation && seedFilterState.numQualitySeeds) {
+    if (m_config.seedConfirmation &&
+        candidates_collector.nHighQualityCandidates()) {
       minCompatibleTopSPs++;
     }
 
@@ -676,9 +677,8 @@ void SeedFinderOrthogonal<external_spacepoint_t>::processFromMiddleSP(
    */
   if ((!bottom_lh_v.empty() && !top_lh_v.empty()) ||
       (!bottom_hl_v.empty() && !top_hl_v.empty())) {
-    m_config.seedFilter->filterSeeds_1SpFixed(
-        spacePointData, candidates_collector, seedFilterState.numQualitySeeds,
-        out_cont);
+    m_config.seedFilter->filterSeeds_1SpFixed(spacePointData,
+                                              candidates_collector, out_cont);
   }
 }
 
diff --git a/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp b/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp
index ae6f5a55209..974a534ae0e 100644
--- a/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp
+++ b/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp
@@ -272,9 +272,8 @@ struct CombinatorialKalmanFilterResult {
 /// The Sequencer has to be part of the Navigator of the Propagator in order to
 /// initialize and provide the measurement surfaces.
 ///
-/// The Actor is part of the Propagation call and does the Kalman update and
-/// eventually the smoothing. Updater and Calibrator are given to the Actor for
-/// further use:
+/// The Actor is part of the Propagation call and does the Kalman update.
+/// Updater and Calibrator are given to the Actor for further use:
 /// - The Updater is the implemented kalman updater formalism, it
 ///   runs via a visitor pattern through the measurements.
 ///
@@ -671,8 +670,7 @@ class CombinatorialKalmanFilter {
                          state.options.stepping.maxStepSize);
 
       // Reset the navigation state
-      // Set targetSurface to nullptr for forward filtering; it's only needed
-      // after smoothing
+      // Set targetSurface to nullptr for forward filtering
       auto navigationOptions = state.navigation.options;
       navigationOptions.startSurface = &currentState.referenceSurface();
       navigationOptions.targetSurface = nullptr;
diff --git a/Core/include/Acts/TrackFitting/GainMatrixSmoother.hpp b/Core/include/Acts/TrackFitting/GainMatrixSmoother.hpp
index 7083cfa35cc..4847364b663 100644
--- a/Core/include/Acts/TrackFitting/GainMatrixSmoother.hpp
+++ b/Core/include/Acts/TrackFitting/GainMatrixSmoother.hpp
@@ -34,13 +34,16 @@ class GainMatrixSmoother {
 
   /// Run the Kalman smoothing for one trajectory.
   ///
+  /// @param[in] gctx The geometry context to be used
   /// @param[in,out] trajectory The trajectory to be smoothed
   /// @param[in] entryIndex The index of state to start the smoothing
   /// @param[in] logger Where to write logging information to
   template <typename traj_t>
-  Result<void> operator()(const GeometryContext& /*gctx*/, traj_t& trajectory,
+  Result<void> operator()(const GeometryContext& gctx, traj_t& trajectory,
                           std::size_t entryIndex,
                           const Logger& logger = getDummyLogger()) const {
+    (void)gctx;
+
     using TrackStateProxy = typename traj_t::TrackStateProxy;
 
     GetParameters filtered;
diff --git a/Core/include/Acts/TrackFitting/GainMatrixUpdater.hpp b/Core/include/Acts/TrackFitting/GainMatrixUpdater.hpp
index c45065799ff..86b6035bad8 100644
--- a/Core/include/Acts/TrackFitting/GainMatrixUpdater.hpp
+++ b/Core/include/Acts/TrackFitting/GainMatrixUpdater.hpp
@@ -24,6 +24,13 @@ namespace Acts {
 /// Kalman update step using the gain matrix formalism.
 class GainMatrixUpdater {
   struct InternalTrackState {
+    unsigned int calibratedSize;
+    // This is used to build a covariance matrix view in the .cpp file
+    const double* calibrated;
+    const double* calibratedCovariance;
+    TrackStateTraits<MultiTrajectoryTraits::MeasurementSizeMax,
+                     false>::Projector projector;
+
     TrackStateTraits<MultiTrajectoryTraits::MeasurementSizeMax,
                      false>::Parameters predicted;
     TrackStateTraits<MultiTrajectoryTraits::MeasurementSizeMax,
@@ -32,12 +39,6 @@ class GainMatrixUpdater {
                      false>::Parameters filtered;
     TrackStateTraits<MultiTrajectoryTraits::MeasurementSizeMax,
                      false>::Covariance filteredCovariance;
-    // This is used to build a covariance matrix view in the .cpp file
-    double* calibrated;
-    double* calibratedCovariance;
-    TrackStateTraits<MultiTrajectoryTraits::MeasurementSizeMax,
-                     false>::Projector projector;
-    unsigned int calibratedSize;
   };
 
  public:
@@ -74,23 +75,16 @@ class GainMatrixUpdater {
 
     auto [chi2, error] = visitMeasurement(
         InternalTrackState{
+            trackState.calibratedSize(),
+            // Note that we pass raw pointers here which are used in the correct
+            // shape later
+            trackState.effectiveCalibrated().data(),
+            trackState.effectiveCalibratedCovariance().data(),
+            trackState.projector(),
             trackState.predicted(),
             trackState.predictedCovariance(),
             trackState.filtered(),
             trackState.filteredCovariance(),
-            // This abuses an incorrectly sized vector / matrix to access the
-            // data pointer! This works (don't use the matrix as is!), but be
-            // careful!
-            trackState
-                .template calibrated<
-                    MultiTrajectoryTraits::MeasurementSizeMax>()
-                .data(),
-            trackState
-                .template calibratedCovariance<
-                    MultiTrajectoryTraits::MeasurementSizeMax>()
-                .data(),
-            trackState.projector(),
-            trackState.calibratedSize(),
         },
         logger);
 
diff --git a/Core/include/Acts/TrackFitting/MbfSmoother.hpp b/Core/include/Acts/TrackFitting/MbfSmoother.hpp
new file mode 100644
index 00000000000..61f28e1d7b4
--- /dev/null
+++ b/Core/include/Acts/TrackFitting/MbfSmoother.hpp
@@ -0,0 +1,160 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2024 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Definitions/TrackParametrization.hpp"
+#include "Acts/EventData/MultiTrajectory.hpp"
+#include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Utilities/Logger.hpp"
+#include "Acts/Utilities/Result.hpp"
+
+#include <cstddef>
+#include <optional>
+
+namespace Acts {
+
+/// Kalman trajectory smoother based on the Modified Bryson–Frazier (mBF)
+/// smoother.
+///
+/// The benefit of the mBF smoother is that it does not require the inverse of
+/// the full covariance matrix, but only the inverse of the residual covariance
+/// matrix which can be cached by the filter step. The same holds for the
+/// Kalman gain matrix.
+///
+/// This implements not a single smoothing step, but the full backwards
+/// smoothing procedure for a filtered, forward trajectory using the stored
+/// linearization.
+///
+/// See
+/// [Wikipedia](https://en.wikipedia.org/wiki/Kalman_filter#Modified_Bryson%E2%80%93Frazier_smoother)
+/// for more information.
+class MbfSmoother {
+ public:
+  /// Run the Kalman smoothing for one trajectory.
+  ///
+  /// @param[in] gctx The geometry context to be used
+  /// @param[in,out] trajectory The trajectory to be smoothed
+  /// @param[in] entryIndex The index of state to start the smoothing
+  /// @param[in] logger Where to write logging information to
+  template <typename traj_t>
+  Result<void> operator()(const GeometryContext& gctx, traj_t& trajectory,
+                          std::size_t entryIndex,
+                          const Logger& logger = getDummyLogger()) const {
+    (void)gctx;
+    (void)logger;
+
+    using TrackStateProxy = typename traj_t::TrackStateProxy;
+
+    TrackStateProxy start_ts = trajectory.getTrackState(entryIndex);
+
+    // Notation consistent with the Wikipedia article
+    // https://en.wikipedia.org/wiki/Kalman_filter
+    BoundMatrix big_lambda_hat = BoundMatrix::Zero();
+    BoundVector small_lambda_hat = BoundVector::Zero();
+
+    trajectory.applyBackwards(start_ts.index(), [&](TrackStateProxy ts) {
+      // ensure the track state has a smoothed component
+      ts.addComponents(TrackStatePropMask::Smoothed);
+
+      InternalTrackState internalTrackState(ts);
+
+      // Smoothe the current state
+      calculateSmoothed(internalTrackState, big_lambda_hat, small_lambda_hat);
+
+      // We smoothed the last state - no need to update the lambdas
+      if (!ts.hasPrevious()) {
+        return;
+      }
+
+      // Update the lambdas depending on the type of track state
+      if (ts.typeFlags().test(TrackStateFlag::MeasurementFlag)) {
+        visitMeasurement(internalTrackState, big_lambda_hat, small_lambda_hat);
+      } else {
+        visitNonMeasurement(internalTrackState, big_lambda_hat,
+                            small_lambda_hat);
+      }
+    });
+
+    return Result<void>::success();
+  }
+
+ private:
+  /// Internal track state representation for the smoother.
+  /// @note This allows us to move parts of the implementation into the .cpp
+  struct InternalTrackState final {
+    using Projector =
+        typename TrackStateTraits<MultiTrajectoryTraits::MeasurementSizeMax,
+                                  false>::Projector;
+    using Jacobian =
+        typename TrackStateTraits<MultiTrajectoryTraits::MeasurementSizeMax,
+                                  false>::Covariance;
+    using Parameters =
+        typename TrackStateTraits<MultiTrajectoryTraits::MeasurementSizeMax,
+                                  false>::Parameters;
+    using Covariance =
+        typename TrackStateTraits<MultiTrajectoryTraits::MeasurementSizeMax,
+                                  false>::Covariance;
+
+    struct Measurement final {
+      unsigned int calibratedSize{0};
+      // This is used to build a covariance matrix view in the .cpp file
+      const double* calibrated{nullptr};
+      const double* calibratedCovariance{nullptr};
+      Projector projector;
+
+      template <typename TrackStateProxy>
+      explicit Measurement(TrackStateProxy ts)
+          : calibratedSize(ts.calibratedSize()),
+            calibrated(ts.effectiveCalibrated().data()),
+            calibratedCovariance(ts.effectiveCalibratedCovariance().data()),
+            projector(ts.projector()) {}
+    };
+
+    Jacobian jacobian;
+
+    Parameters predicted;
+    Covariance predictedCovariance;
+    Parameters filtered;
+    Covariance filteredCovariance;
+    Parameters smoothed;
+    Covariance smoothedCovariance;
+
+    std::optional<Measurement> measurement;
+
+    template <typename TrackStateProxy>
+    explicit InternalTrackState(TrackStateProxy ts)
+        : jacobian(ts.jacobian()),
+          predicted(ts.predicted()),
+          predictedCovariance(ts.predictedCovariance()),
+          filtered(ts.filtered()),
+          filteredCovariance(ts.filteredCovariance()),
+          smoothed(ts.smoothed()),
+          smoothedCovariance(ts.smoothedCovariance()),
+          measurement(ts.typeFlags().test(TrackStateFlag::MeasurementFlag)
+                          ? std::optional<Measurement>(ts)
+                          : std::nullopt) {}
+  };
+
+  /// Calculate the smoothed parameters and covariance.
+  void calculateSmoothed(InternalTrackState& ts,
+                         const BoundMatrix& big_lambda_hat,
+                         const BoundVector& small_lambda_hat) const;
+
+  /// Visit a non-measurement track state and update the lambdas.
+  void visitNonMeasurement(const InternalTrackState& ts,
+                           BoundMatrix& big_lambda_hat,
+                           BoundVector& small_lambda_hat) const;
+
+  /// Visit a measurement track state and update the lambdas.
+  void visitMeasurement(const InternalTrackState& ts,
+                        BoundMatrix& big_lambda_hat,
+                        BoundVector& small_lambda_hat) const;
+};
+
+}  // namespace Acts
diff --git a/Core/include/Acts/Utilities/AnnealingUtility.hpp b/Core/include/Acts/Utilities/AnnealingUtility.hpp
index d7075355e4f..204f6c59008 100644
--- a/Core/include/Acts/Utilities/AnnealingUtility.hpp
+++ b/Core/include/Acts/Utilities/AnnealingUtility.hpp
@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2019 CERN for the benefit of the Acts project
+// Copyright (C) 2019-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -32,6 +32,9 @@ class AnnealingUtility {
   struct Config {
     Config();
 
+    Config(double cutOff_, std::vector<double> setOfTemperatures_)
+        : cutOff(cutOff_), setOfTemperatures(std::move(setOfTemperatures_)) {}
+
     // Insensitivity of calculated weight at cutoff
     double cutOff{9.};
 
diff --git a/Core/include/Acts/Utilities/Intersection.hpp b/Core/include/Acts/Utilities/Intersection.hpp
index f30da706c07..49dc23517b7 100644
--- a/Core/include/Acts/Utilities/Intersection.hpp
+++ b/Core/include/Acts/Utilities/Intersection.hpp
@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2023 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -60,10 +60,15 @@ class Intersection {
       : m_position(position), m_pathLength(pathLength), m_status(status) {}
 
   /// Returns whether the intersection was successful or not
-  constexpr explicit operator bool() const {
-    return m_status != Status::missed;
+  /// @deprecated
+  [[deprecated("Use isValid() instead")]] constexpr explicit operator bool()
+      const {
+    return isValid();
   }
 
+  /// Returns whether the intersection was successful or not
+  constexpr bool isValid() const { return m_status != Status::missed; }
+
   constexpr const Position& position() const { return m_position; }
 
   constexpr ActsScalar pathLength() const { return m_pathLength; }
@@ -141,10 +146,15 @@ class ObjectIntersection {
       : m_intersection(intersection), m_object(object), m_index(index) {}
 
   /// Returns whether the intersection was successful or not
-  constexpr explicit operator bool() const {
-    return m_intersection.operator bool();
+  /// @deprecated
+  [[deprecated("Use isValid() instead")]] constexpr explicit operator bool()
+      const {
+    return isValid();
   }
 
+  /// Returns whether the intersection was successful or not
+  constexpr bool isValid() const { return m_intersection.isValid(); }
+
   constexpr const Intersection3D& intersection() const {
     return m_intersection;
   }
@@ -256,50 +266,16 @@ class ObjectMultiIntersection {
 
 namespace detail {
 
-/// This function checks if an intersection is valid for the specified
-/// path-limit and overstep-limit
-///
-/// @tparam intersection_t Type of the intersection object
+/// This function checks if an intersection path length is valid for the
+/// specified near-limit and far-limit
 ///
-/// @param intersection The intersection to check
-/// @param nearLimit The minimum distance for an intersection to be considered
-/// @param farLimit The maximum distance for an intersection to be considered
+/// @param pathLength The path length of the intersection
+/// @param nearLimit The minimum path length for an intersection to be considered
+/// @param farLimit The maximum path length for an intersection to be considered
 /// @param logger A optionally supplied logger which prints out a lot of infos
 ///               at VERBOSE level
-template <typename intersection_t>
-bool checkIntersection(const intersection_t& intersection, double nearLimit,
-                       double farLimit,
-                       const Logger& logger = getDummyLogger()) {
-  const double distance = intersection.pathLength();
-  // TODO why?
-  const double tolerance = s_onSurfaceTolerance;
-
-  ACTS_VERBOSE(" -> near limit, far limit, distance: "
-               << nearLimit << ", " << farLimit << ", " << distance);
-
-  const bool coCriterion = distance > nearLimit;
-  const bool cpCriterion = distance < farLimit + tolerance;
-
-  const bool accept = coCriterion && cpCriterion;
-
-  if (accept) {
-    ACTS_VERBOSE("Intersection is WITHIN limit");
-  } else {
-    ACTS_VERBOSE("Intersection is OUTSIDE limit because: ");
-    if (!coCriterion) {
-      ACTS_VERBOSE("- intersection path length "
-                   << distance << " <= near limit " << nearLimit);
-    }
-    if (!cpCriterion) {
-      ACTS_VERBOSE("- intersection path length "
-                   << distance << " is over the far limit "
-                   << (farLimit + tolerance) << " (including tolerance of "
-                   << tolerance << ")");
-    }
-  }
-
-  return accept;
-}
+bool checkPathLength(double pathLength, double nearLimit, double farLimit,
+                     const Logger& logger = getDummyLogger());
 
 }  // namespace detail
 
diff --git a/Core/include/Acts/Utilities/TrackHelpers.hpp b/Core/include/Acts/Utilities/TrackHelpers.hpp
index 46e9b757857..f9883bfeb63 100644
--- a/Core/include/Acts/Utilities/TrackHelpers.hpp
+++ b/Core/include/Acts/Utilities/TrackHelpers.hpp
@@ -85,18 +85,19 @@ Result<typename track_proxy_t::ConstTrackStateProxy> findLastMeasurementState(
 /// @brief Smooth a track using the gain matrix smoother
 ///
 /// @tparam track_proxy_t The track proxy type
+/// @tparam smoother_t The smoother type
 ///
 /// @param geoContext The geometry context
 /// @param track The track to smooth
 /// @param logger The logger
+/// @param smoother The smoother
 ///
 /// @return The result of the smoothing
-template <typename track_proxy_t>
+template <typename track_proxy_t, typename smoother_t = GainMatrixSmoother>
 Result<void> smoothTrack(
     const GeometryContext &geoContext, track_proxy_t &track,
-    const Logger &logger = *getDefaultLogger("TrackSmoother", Logging::INFO)) {
-  Acts::GainMatrixSmoother smoother;
-
+    const Logger &logger = *getDefaultLogger("TrackSmoother", Logging::INFO),
+    smoother_t smoother = GainMatrixSmoother()) {
   auto &trackContainer = track.container();
   auto &trackStateContainer = trackContainer.trackStateContainer();
 
@@ -194,7 +195,7 @@ findTrackStateForExtrapolation(
       }
 
       SurfaceIntersection intersection = intersect(*first);
-      if (!intersection) {
+      if (!intersection.isValid()) {
         ACTS_ERROR("no intersection found");
         return Result<std::pair<TrackStateProxy, double>>::failure(
             TrackExtrapolationError::ReferenceSurfaceUnreachable);
@@ -214,7 +215,7 @@ findTrackStateForExtrapolation(
       }
 
       SurfaceIntersection intersection = intersect(*last);
-      if (!intersection) {
+      if (!intersection.isValid()) {
         ACTS_ERROR("no intersection found");
         return Result<std::pair<TrackStateProxy, double>>::failure(
             TrackExtrapolationError::ReferenceSurfaceUnreachable);
@@ -245,13 +246,13 @@ findTrackStateForExtrapolation(
       double absDistanceFirst = std::abs(intersectionFirst.pathLength());
       double absDistanceLast = std::abs(intersectionLast.pathLength());
 
-      if (intersectionFirst && absDistanceFirst <= absDistanceLast) {
+      if (intersectionFirst.isValid() && absDistanceFirst <= absDistanceLast) {
         ACTS_VERBOSE("using first track state with intersection at "
                      << intersectionFirst.pathLength());
         return std::make_pair(*first, intersectionFirst.pathLength());
       }
 
-      if (intersectionLast && absDistanceLast <= absDistanceFirst) {
+      if (intersectionLast.isValid() && absDistanceLast <= absDistanceFirst) {
         ACTS_VERBOSE("using last track state with intersection at "
                      << intersectionLast.pathLength());
         return std::make_pair(*last, intersectionLast.pathLength());
diff --git a/Core/include/Acts/Utilities/detail/ContextType.hpp b/Core/include/Acts/Utilities/detail/ContextType.hpp
index 029769e292e..486e07d425d 100644
--- a/Core/include/Acts/Utilities/detail/ContextType.hpp
+++ b/Core/include/Acts/Utilities/detail/ContextType.hpp
@@ -46,7 +46,7 @@ class ContextType {
   /// @return ContextType&
   template <typename T>
   ContextType& operator=(T&& value) {
-    m_data = std::move(value);
+    m_data = std::forward<T>(value);
     return *this;
   }
 
diff --git a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.hpp b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.hpp
index 0bf05086b67..818d2c363be 100644
--- a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.hpp
+++ b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.hpp
@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2020 CERN for the benefit of the Acts project
+// Copyright (C) 2020-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -25,6 +25,7 @@
 #include <type_traits>
 
 namespace Acts {
+
 /// @brief Implements an iterative vertex finder
 class AdaptiveMultiVertexFinder final : public IVertexFinder {
   using VertexFitter = AdaptiveMultiVertexFitter;
@@ -47,113 +48,120 @@ class AdaptiveMultiVertexFinder final : public IVertexFinder {
           ipEstimator(std::move(ipEst)),
           bField{std::move(bIn)} {}
 
-    // Vertex fitter
+    /// Vertex fitter
     VertexFitter vertexFitter;
 
-    // Vertex seed finder
+    /// Vertex seed finder
     std::shared_ptr<const IVertexFinder> seedFinder;
 
-    // ImpactPointEstimator
+    /// ImpactPointEstimator
     ImpactPointEstimator ipEstimator;
 
     std::shared_ptr<const MagneticFieldProvider> bField;
 
-    // Max z interval used for adding tracks to fit:
-    // When adding a new vertex to the multi vertex fit,
-    // only the tracks whose z at PCA is closer
-    // to the seeded vertex than tracksMaxZinterval
-    // are added to this new vertex.
-    //
-    // Note: If you cut too hard, you cut out
-    // the good cases where the seed finder is not
-    // reliable, but the fit would be still able to converge
-    // towards the right vertex. If you cut too soft, you
-    // consider a lot of tracks which just slow down the fit.
+    /// Max z interval used for adding tracks to fit:
+    /// When adding a new vertex to the multi vertex fit,
+    /// only the tracks whose z at PCA is closer
+    /// to the seeded vertex than tracksMaxZinterval
+    /// are added to this new vertex.
+    ///
+    /// Note: If you cut too hard, you cut out
+    /// the good cases where the seed finder is not
+    /// reliable, but the fit would be still able to converge
+    /// towards the right vertex. If you cut too soft, you
+    /// consider a lot of tracks which just slow down the fit.
     double tracksMaxZinterval = 3. * Acts::UnitConstants::mm;
 
-    // Maximum allowed significance of track position to vertex seed to consider
-    // track as compatible to vertex. If useTime is set to true, the time
-    // coordinate also contributes to the significance and tracksMaxSignificance
-    // needs to be increased.
-    // 5 corresponds to a p-value of ~0.92 using `chi2(x=5,ndf=2)`
+    /// Maximum allowed significance of track position to vertex seed to
+    /// consider track as compatible to vertex. If useTime is set to true, the
+    /// time coordinate also contributes to the significance and
+    /// tracksMaxSignificance needs to be increased. 5 corresponds to a p-value
+    /// of ~0.92 using `chi2(x=5,ndf=2)`
     double tracksMaxSignificance = 5.;
 
-    // Max chi2 value for which tracks are considered compatible with
-    // the fitted vertex. These tracks are removed from the seedTracks
-    // after the fit has been performed.
+    /// Max chi2 value for which tracks are considered compatible with
+    /// the fitted vertex. These tracks are removed from the seedTracks
+    /// after the fit has been performed.
     double maxVertexChi2 = 18.42;
 
-    // Perform a 'real' multi-vertex fit as intended by the algorithm.
-    // If switched to true, always all (!) tracks are considered to be
-    // added to the new vertex candidate after seeding. If switched to
-    // false, only the seedTracks, i.e. all tracks that are considered
-    // as outliers of previously fitted vertices, are used.
+    /// Perform a 'real' multi-vertex fit as intended by the algorithm.
+    /// If switched to true, always all (!) tracks are considered to be
+    /// added to the new vertex candidate after seeding. If switched to
+    /// false, only the seedTracks, i.e. all tracks that are considered
+    /// as outliers of previously fitted vertices, are used.
     bool doRealMultiVertex = true;
 
-    // Decides if you want to use the ```vertexCompatibility``` of the
-    //  track (set to true) or the ```chi2Track``` (set to false) as an
-    // estimate for a track being an outlier or not.
-    // In case the track refitting is switched on in the AMVFitter, you
-    // may want to use the refitted ```chi2Track```.
+    /// Decides if you want to use the ```vertexCompatibility``` of the
+    ///  track (set to true) or the ```chi2Track``` (set to false) as an
+    /// estimate for a track being an outlier or not.
+    /// In case the track refitting is switched on in the AMVFitter, you
+    /// may want to use the refitted ```chi2Track```.
     bool useFastCompatibility = true;
 
-    // Maximum significance on the distance between two vertices
-    // to allow merging of two vertices.
-    // 3 corresponds to a p-value of ~0.92 using `chi2(x=3,ndf=1)`
+    /// Maximum significance on the distance between two vertices
+    /// to allow merging of two vertices.
+    /// 3 corresponds to a p-value of ~0.92 using `chi2(x=3,ndf=1)`
     double maxMergeVertexSignificance = 3.;
 
-    // Minimum weight a track has to have to be considered a compatible
-    // track with a vertex candidate.
-    //
-    // Note: This value has to be the same as the one in the AMVFitter.
+    /// Minimum weight a track has to have to be considered a compatible
+    /// track with a vertex candidate.
+    ///
+    /// Note: This value has to be the same as the one in the AMVFitter.
     double minWeight = 0.0001;
 
-    // Maximal number of iterations in the finding procedure
+    /// Maximal number of iterations in the finding procedure
     int maxIterations = 100;
 
-    // Include also single track vertices
+    /// Include also single track vertices
     bool addSingleTrackVertices = false;
 
-    // If doFullSplitting == true, we check the 3D distance (if useTime ==
-    // false) or the 4D distance (if useTime == true) of the vertices to
-    // determine whether they are merged.
-    // If doFullSplitting == false, we check the z distance (if useTime ==
-    // false) or the z-t distance (if useTime == true) of the vertices to
-    // determine whether they are merged.
+    /// If doFullSplitting == true, we check the 3D distance (if useTime ==
+    /// false) or the 4D distance (if useTime == true) of the vertices to
+    /// determine whether they are merged.
+    /// If doFullSplitting == false, we check the z distance (if useTime ==
+    /// false) or the z-t distance (if useTime == true) of the vertices to
+    /// determine whether they are merged.
     bool doFullSplitting = false;
 
-    // Maximum vertex contamination value
+    /// Maximum vertex contamination value
     double maximumVertexContamination = 0.5;
 
-    // Use seed vertex as a constraint for the fit
+    /// Use seed vertex as a constraint for the fit
     bool useSeedConstraint = true;
 
-    // Variances of the 4D vertex position before the vertex fit if no beamspot
-    // constraint is provided
+    /// Variances of the 4D vertex position before the vertex fit if no beamspot
+    /// constraint is provided
     Vector4 initialVariances = Vector4::Constant(1e+8);
 
-    // Default fitQuality for constraint vertex in case no beamspot
-    // constraint is provided
+    /// Default fitQuality for constraint vertex in case no beamspot
+    /// constraint is provided
     std::pair<double, double> defaultConstrFitQuality{0., -3.};
 
-    // Use the full available vertex covariance information after
-    // seeding for the IP estimation. In original implementation
-    // this is not (!) done, however, this is probably not correct.
-    // So definitely consider setting this to true.
+    /// Use the full available vertex covariance information after
+    /// seeding for the IP estimation. In original implementation
+    /// this is not (!) done, however, this is probably not correct.
+    /// So definitely consider setting this to true.
     bool useVertexCovForIPEstimation = false;
 
-    // Use time information when assigning tracks to vertices. If this is set to
-    // true, useTime of the vertex fitter configuration should also be set to
-    // true, and time seeding should be enabled.
+    /// Use time information when assigning tracks to vertices. If this is set
+    /// to true, useTime of the vertex fitter configuration should also be set
+    /// to true, and time seeding should be enabled.
     bool useTime = false;
 
-    // Function to extract parameters from InputTrack
+    /// If set to true, the vertex finder will not break the finding loop.
+    /// Some seeders are not able to cope with this therefore this is
+    /// disabled by default.
+    bool doNotBreakWhileSeeding = false;
+
+    /// Function to extract parameters from InputTrack
     InputTrack::Extractor extractParameters;
-  };  // Config struct
+  };
 
   /// State struct for fulfilling interface
   struct State {
     std::reference_wrapper<const MagneticFieldContext> magContext;
+
+    IVertexFinder::State seedFinderState;
   };
 
   /// @brief Constructor for user-defined InputTrack_t type !=
@@ -197,7 +205,8 @@ class AdaptiveMultiVertexFinder final : public IVertexFinder {
 
   IVertexFinder::State makeState(
       const Acts::MagneticFieldContext& mctx) const override {
-    return IVertexFinder::State{State{mctx}};
+    return IVertexFinder::State{
+        State{mctx, IVertexFinder::State{m_cfg.seedFinder->makeState(mctx)}}};
   }
 
   void setTracksToRemove(
diff --git a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.hpp b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.hpp
index cc6ebfb008d..2bd3827faa7 100644
--- a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.hpp
+++ b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.hpp
@@ -41,6 +41,7 @@ class AdaptiveMultiVertexFitter {
           const Acts::MagneticFieldContext& magContext)
         : ipState{field.makeCache(magContext)},
           fieldCache(field.makeCache(magContext)) {}
+
     // Vertex collection to be fitted
     std::vector<Vertex*> vertexCollection;
 
diff --git a/Core/include/Acts/Vertexing/Vertex.hpp b/Core/include/Acts/Vertexing/Vertex.hpp
index c0e6c62d839..1ce6e0ae7c5 100644
--- a/Core/include/Acts/Vertexing/Vertex.hpp
+++ b/Core/include/Acts/Vertexing/Vertex.hpp
@@ -73,11 +73,19 @@ class Vertex {
   /// @return Returns pair of (chi2, numberDoF)
   std::pair<double, double> fitQuality() const;
 
+  /// @brief Set position
+  ///
+  /// @param position Vertex position
+  void setPosition(const Vector3& position);
+
   /// @brief Set position and time
   ///
+  /// @deprecated Use setFullPosition instead
+  ///
   /// @param position Vertex position
   /// @param time The time
-  void setPosition(const Vector3& position, ActsScalar time = 0);
+  [[deprecated("Use setFullPosition instead")]] void setPosition(
+      const Vector3& position, ActsScalar time);
 
   /// @brief Set position and time
   ///
diff --git a/Core/src/EventData/VectorTrackContainer.cpp b/Core/src/EventData/VectorTrackContainer.cpp
index 27f7425182a..663634af5f2 100644
--- a/Core/src/EventData/VectorTrackContainer.cpp
+++ b/Core/src/EventData/VectorTrackContainer.cpp
@@ -80,6 +80,8 @@ void VectorTrackContainer::removeTrack_impl(IndexType itrack) {
   erase(m_tipIndex);
   erase(m_stemIndex);
 
+  erase(m_particleHypothesis);
+
   erase(m_params);
   erase(m_cov);
   erase(m_referenceSurfaces);
diff --git a/Core/src/Geometry/GenericApproachDescriptor.cpp b/Core/src/Geometry/GenericApproachDescriptor.cpp
index cc1320f19e1..e532e2f5eb2 100644
--- a/Core/src/Geometry/GenericApproachDescriptor.cpp
+++ b/Core/src/Geometry/GenericApproachDescriptor.cpp
@@ -35,8 +35,9 @@ Acts::SurfaceIntersection Acts::GenericApproachDescriptor::approachSurface(
     auto sfIntersection =
         sf->intersect(gctx, position, direction, boundaryTolerance);
     for (const auto& intersection : sfIntersection.split()) {
-      if (intersection &&
-          detail::checkIntersection(intersection, nearLimit, farLimit)) {
+      if (intersection.isValid() &&
+          detail::checkPathLength(intersection.pathLength(), nearLimit,
+                                  farLimit)) {
         sIntersections.push_back(intersection);
       }
     }
diff --git a/Core/src/Geometry/Layer.cpp b/Core/src/Geometry/Layer.cpp
index 0f1dee70659..c88bfccb301 100644
--- a/Core/src/Geometry/Layer.cpp
+++ b/Core/src/Geometry/Layer.cpp
@@ -14,6 +14,7 @@
 #include "Acts/Propagator/Navigator.hpp"
 #include "Acts/Surfaces/BoundaryTolerance.hpp"
 #include "Acts/Surfaces/Surface.hpp"
+#include "Acts/Utilities/Intersection.hpp"
 
 #include <algorithm>
 #include <functional>
@@ -138,7 +139,7 @@ Acts::Layer::compatibleSurfaces(
     // indicates wrong direction or faulty setup
     // -> do not return compatible surfaces since they may lead you on a wrong
     // navigation path
-    if (endInter) {
+    if (endInter.isValid()) {
       farLimit = endInter.pathLength();
     } else {
       return sIntersections;
@@ -196,8 +197,8 @@ Acts::Layer::compatibleSurfaces(
     // the surface intersection
     SurfaceIntersection sfi =
         sf.intersect(gctx, position, direction, boundaryTolerance).closest();
-    if (sfi &&
-        detail::checkIntersection(sfi.intersection(), nearLimit, farLimit) &&
+    if (sfi.isValid() &&
+        detail::checkPathLength(sfi.pathLength(), nearLimit, farLimit) &&
         isUnique(sfi)) {
       sIntersections.push_back(sfi);
     }
@@ -267,8 +268,8 @@ Acts::SurfaceIntersection Acts::Layer::surfaceOnApproach(
   auto findValidIntersection =
       [&](const SurfaceMultiIntersection& sfmi) -> SurfaceIntersection {
     for (const auto& sfi : sfmi.split()) {
-      if (sfi &&
-          detail::checkIntersection(sfi.intersection(), nearLimit, farLimit)) {
+      if (sfi.isValid() &&
+          detail::checkPathLength(sfi.pathLength(), nearLimit, farLimit)) {
         return sfi;
       }
     }
diff --git a/Core/src/Geometry/TrackingVolume.cpp b/Core/src/Geometry/TrackingVolume.cpp
index 8b34abf912e..9cc4afc2a37 100644
--- a/Core/src/Geometry/TrackingVolume.cpp
+++ b/Core/src/Geometry/TrackingVolume.cpp
@@ -19,6 +19,7 @@
 #include "Acts/Surfaces/RegularSurface.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/BinningType.hpp"
+#include "Acts/Utilities/Intersection.hpp"
 
 #include <algorithm>
 #include <memory>
@@ -440,7 +441,7 @@ TrackingVolume::compatibleBoundaries(const GeometryContext& gctx,
       [&](SurfaceMultiIntersection& candidates,
           const BoundarySurface* boundary) -> BoundaryIntersection {
     for (const auto& intersection : candidates.split()) {
-      if (!intersection) {
+      if (!intersection.isValid()) {
         continue;
       }
 
@@ -464,8 +465,8 @@ TrackingVolume::compatibleBoundaries(const GeometryContext& gctx,
 
       ACTS_VERBOSE("Check intersection with surface "
                    << boundary->surfaceRepresentation().geometryId());
-      if (detail::checkIntersection(intersection.intersection(), nearLimit,
-                                    farLimit, logger)) {
+      if (detail::checkPathLength(intersection.pathLength(), nearLimit,
+                                  farLimit, logger)) {
         return BoundaryIntersection(intersection, boundary);
       }
     }
@@ -495,7 +496,7 @@ TrackingVolume::compatibleBoundaries(const GeometryContext& gctx,
                                           options.boundaryTolerance);
       // Intersect and continue
       auto intersection = checkIntersection(candidates, boundary.get());
-      if (intersection.first) {
+      if (intersection.first.isValid()) {
         ACTS_VERBOSE(" - Proceed with surface");
         intersections.push_back(intersection);
       } else {
@@ -550,7 +551,7 @@ TrackingVolume::compatibleLayers(
       auto atIntersection =
           tLayer->surfaceOnApproach(gctx, position, direction, options);
       // Intersection is ok - take it (move to surface on approach)
-      if (atIntersection) {
+      if (atIntersection.isValid()) {
         // create a layer intersection
         lIntersections.push_back(LayerIntersection(atIntersection, tLayer));
       }
diff --git a/Core/src/TrackFitting/CMakeLists.txt b/Core/src/TrackFitting/CMakeLists.txt
index c370dc9aa01..922636b6d64 100644
--- a/Core/src/TrackFitting/CMakeLists.txt
+++ b/Core/src/TrackFitting/CMakeLists.txt
@@ -10,4 +10,5 @@ target_sources(
     BetheHeitlerApprox.cpp
     GsfMixtureReduction.cpp
     GlobalChiSquareFitter.cpp
+    MbfSmoother.cpp
 )
diff --git a/Core/src/TrackFitting/GainMatrixSmoother.cpp b/Core/src/TrackFitting/GainMatrixSmoother.cpp
index 464fecb20c6..d5a79f57d59 100644
--- a/Core/src/TrackFitting/GainMatrixSmoother.cpp
+++ b/Core/src/TrackFitting/GainMatrixSmoother.cpp
@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2021 CERN for the benefit of the Acts project
+// Copyright (C) 2021-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
diff --git a/Core/src/TrackFitting/MbfSmoother.cpp b/Core/src/TrackFitting/MbfSmoother.cpp
new file mode 100644
index 00000000000..fe668eb41d8
--- /dev/null
+++ b/Core/src/TrackFitting/MbfSmoother.cpp
@@ -0,0 +1,79 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2024 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "Acts/TrackFitting/MbfSmoother.hpp"
+
+namespace Acts {
+
+void MbfSmoother::calculateSmoothed(InternalTrackState& ts,
+                                    const BoundMatrix& big_lambda_hat,
+                                    const BoundVector& small_lambda_hat) const {
+  ts.smoothedCovariance = ts.filteredCovariance - ts.filteredCovariance *
+                                                      big_lambda_hat *
+                                                      ts.filteredCovariance;
+  ts.smoothed = ts.filtered - ts.filteredCovariance * small_lambda_hat;
+}
+
+void MbfSmoother::visitNonMeasurement(const InternalTrackState& ts,
+                                      BoundMatrix& big_lambda_hat,
+                                      BoundVector& small_lambda_hat) const {
+  const auto F = ts.jacobian;
+
+  big_lambda_hat = F.transpose() * big_lambda_hat * F;
+  small_lambda_hat = F.transpose() * small_lambda_hat;
+}
+
+void MbfSmoother::visitMeasurement(const InternalTrackState& ts,
+                                   BoundMatrix& big_lambda_hat,
+                                   BoundVector& small_lambda_hat) const {
+  assert(ts.measurement.has_value());
+
+  const auto& measurement = ts.measurement.value();
+
+  visit_measurement(measurement.calibratedSize, [&](auto N) -> void {
+    constexpr std::size_t kMeasurementSize = decltype(N)::value;
+
+    typename TrackStateTraits<kMeasurementSize, true>::Calibrated calibrated{
+        measurement.calibrated};
+    typename TrackStateTraits<kMeasurementSize, true>::CalibratedCovariance
+        calibratedCovariance{measurement.calibratedCovariance};
+
+    // Measurement matrix
+    const auto H = measurement.projector
+                       .template topLeftCorner<kMeasurementSize, eBoundSize>()
+                       .eval();
+
+    // Residual covariance
+    const auto S =
+        (H * ts.predictedCovariance * H.transpose() + calibratedCovariance)
+            .eval();
+    // TODO Sinv could be cached by the filter step
+    const auto S_inv = S.inverse().eval();
+
+    // Kalman gain
+    // TODO K could be cached by the filter step
+    const auto K = (ts.predictedCovariance * H.transpose() * S_inv).eval();
+
+    const auto C_hat = (BoundMatrix::Identity() - K * H).eval();
+    const auto y = (calibrated - H * ts.predicted).eval();
+
+    const auto big_lambda_tilde =
+        (H.transpose() * S_inv * H + C_hat.transpose() * big_lambda_hat * C_hat)
+            .eval();
+    const auto small_lambda_tilde =
+        (-H.transpose() * S_inv * y + C_hat.transpose() * small_lambda_hat)
+            .eval();
+
+    const auto F = ts.jacobian;
+
+    big_lambda_hat = F.transpose() * big_lambda_tilde * F;
+    small_lambda_hat = F.transpose() * small_lambda_tilde;
+  });
+}
+
+}  // namespace Acts
diff --git a/Core/src/Utilities/CMakeLists.txt b/Core/src/Utilities/CMakeLists.txt
index 64ede8dcd89..17d0492d51e 100644
--- a/Core/src/Utilities/CMakeLists.txt
+++ b/Core/src/Utilities/CMakeLists.txt
@@ -7,4 +7,5 @@ target_sources(
     SpacePointUtility.cpp
     TrackHelpers.cpp
     BinningType.cpp
+    Intersection.cpp
 )
diff --git a/Core/src/Utilities/Intersection.cpp b/Core/src/Utilities/Intersection.cpp
new file mode 100644
index 00000000000..2c9531bdfdf
--- /dev/null
+++ b/Core/src/Utilities/Intersection.cpp
@@ -0,0 +1,45 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2024 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "Acts/Utilities/Intersection.hpp"
+
+namespace Acts {
+
+bool detail::checkPathLength(double pathLength, double nearLimit,
+                             double farLimit, const Logger& logger) {
+  // TODO why?
+  const double tolerance = s_onSurfaceTolerance;
+
+  ACTS_VERBOSE(" -> near limit, far limit, distance: "
+               << nearLimit << ", " << farLimit << ", " << pathLength);
+
+  const bool coCriterion = pathLength > nearLimit;
+  const bool cpCriterion = pathLength < farLimit + tolerance;
+
+  const bool accept = coCriterion && cpCriterion;
+
+  if (accept) {
+    ACTS_VERBOSE("Intersection is WITHIN limit");
+  } else {
+    ACTS_VERBOSE("Intersection is OUTSIDE limit because: ");
+    if (!coCriterion) {
+      ACTS_VERBOSE("- intersection path length "
+                   << pathLength << " <= near limit " << nearLimit);
+    }
+    if (!cpCriterion) {
+      ACTS_VERBOSE("- intersection path length "
+                   << pathLength << " is over the far limit "
+                   << (farLimit + tolerance) << " (including tolerance of "
+                   << tolerance << ")");
+    }
+  }
+
+  return accept;
+}
+
+}  // namespace Acts
diff --git a/Core/src/Vertexing/AdaptiveMultiVertexFinder.cpp b/Core/src/Vertexing/AdaptiveMultiVertexFinder.cpp
index e410bc2fb01..c6308c17152 100644
--- a/Core/src/Vertexing/AdaptiveMultiVertexFinder.cpp
+++ b/Core/src/Vertexing/AdaptiveMultiVertexFinder.cpp
@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2020-2023 CERN for the benefit of the Acts project
+// Copyright (C) 2020-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -9,11 +9,12 @@
 #include "Acts/Vertexing/AdaptiveMultiVertexFinder.hpp"
 
 #include "Acts/Utilities/AlgebraHelpers.hpp"
+#include "Acts/Vertexing/IVertexFinder.hpp"
 #include "Acts/Vertexing/VertexingError.hpp"
 
 namespace Acts {
 
-Acts::Result<std::vector<Acts::Vertex>> AdaptiveMultiVertexFinder::find(
+Result<std::vector<Vertex>> AdaptiveMultiVertexFinder::find(
     const std::vector<InputTrack>& allTracks,
     const VertexingOptions& vertexingOptions,
     IVertexFinder::State& anyState) const {
@@ -23,19 +24,13 @@ Acts::Result<std::vector<Acts::Vertex>> AdaptiveMultiVertexFinder::find(
   }
 
   State& state = anyState.template as<State>();
-
-  // Original tracks
-  const std::vector<InputTrack>& origTracks = allTracks;
-
-  // Seed tracks
-  std::vector<InputTrack> seedTracks = allTracks;
-
+  IVertexFinder::State& seedFinderState = state.seedFinderState;
   VertexFitterState fitterState(*m_cfg.bField,
                                 vertexingOptions.magFieldContext);
-  auto seedFinderState = m_cfg.seedFinder->makeState(state.magContext);
 
+  const std::vector<InputTrack>& origTracks = allTracks;
+  std::vector<InputTrack> seedTracks = allTracks;
   std::vector<std::unique_ptr<Vertex>> allVertices;
-
   std::vector<Vertex*> allVerticesPtr;
 
   int iteration = 0;
@@ -91,7 +86,11 @@ Acts::Result<std::vector<Acts::Vertex>> AdaptiveMultiVertexFinder::find(
           "Could not prepare for fit. Discarding the vertex candindate.");
       allVertices.pop_back();
       allVerticesPtr.pop_back();
-      break;
+      if (m_cfg.doNotBreakWhileSeeding) {
+        continue;
+      } else {
+        break;
+      }
     }
     // Update fitter state with all vertices
     fitterState.addVertexToMultiMap(vtxCandidate);
@@ -201,7 +200,7 @@ void AdaptiveMultiVertexFinder::setConstraintAfterSeeding(
   }
 }
 
-Acts::Result<double> AdaptiveMultiVertexFinder::getIPSignificance(
+Result<double> AdaptiveMultiVertexFinder::getIPSignificance(
     const InputTrack& track, const Vertex& vtx,
     const VertexingOptions& vertexingOptions) const {
   // TODO: In original implementation the covariance of the given vertex is set
@@ -240,7 +239,7 @@ Acts::Result<double> AdaptiveMultiVertexFinder::getIPSignificance(
   return significance;
 }
 
-Acts::Result<void> AdaptiveMultiVertexFinder::addCompatibleTracksToVertex(
+Result<void> AdaptiveMultiVertexFinder::addCompatibleTracksToVertex(
     const std::vector<InputTrack>& tracks, Vertex& vtx,
     VertexFitterState& fitterState,
     const VertexingOptions& vertexingOptions) const {
@@ -269,7 +268,7 @@ Acts::Result<void> AdaptiveMultiVertexFinder::addCompatibleTracksToVertex(
   return {};
 }
 
-Acts::Result<bool> AdaptiveMultiVertexFinder::canRecoverFromNoCompatibleTracks(
+Result<bool> AdaptiveMultiVertexFinder::canRecoverFromNoCompatibleTracks(
     const std::vector<InputTrack>& allTracks,
     const std::vector<InputTrack>& seedTracks, Vertex& vtx,
     const Vertex& currentConstraint, VertexFitterState& fitterState,
@@ -323,7 +322,7 @@ Acts::Result<bool> AdaptiveMultiVertexFinder::canRecoverFromNoCompatibleTracks(
   return Result<bool>::success(true);
 }
 
-Acts::Result<bool> AdaptiveMultiVertexFinder::canPrepareVertexForFit(
+Result<bool> AdaptiveMultiVertexFinder::canPrepareVertexForFit(
     const std::vector<InputTrack>& allTracks,
     const std::vector<InputTrack>& seedTracks, Vertex& vtx,
     const Vertex& currentConstraint, VertexFitterState& fitterState,
@@ -577,7 +576,7 @@ Result<bool> AdaptiveMultiVertexFinder::isMergedVertex(
   return Result<bool>::success(false);
 }
 
-Acts::Result<void> AdaptiveMultiVertexFinder::deleteLastVertex(
+Result<void> AdaptiveMultiVertexFinder::deleteLastVertex(
     Vertex& vtx, std::vector<std::unique_ptr<Vertex>>& allVertices,
     std::vector<Vertex*>& allVerticesPtr, VertexFitterState& fitterState,
     const VertexingOptions& vertexingOptions) const {
@@ -616,8 +615,7 @@ Acts::Result<void> AdaptiveMultiVertexFinder::deleteLastVertex(
   return {};
 }
 
-Acts::Result<std::vector<Acts::Vertex>>
-AdaptiveMultiVertexFinder::getVertexOutputList(
+Result<std::vector<Vertex>> AdaptiveMultiVertexFinder::getVertexOutputList(
     const std::vector<Vertex*>& allVerticesPtr,
     VertexFitterState& fitterState) const {
   std::vector<Vertex> outputVec;
diff --git a/Core/src/Vertexing/Vertex.cpp b/Core/src/Vertexing/Vertex.cpp
index 637ecc07eaf..f8373ba774a 100644
--- a/Core/src/Vertexing/Vertex.cpp
+++ b/Core/src/Vertexing/Vertex.cpp
@@ -77,6 +77,10 @@ std::pair<double, double> Vertex::fitQuality() const {
   return std::pair<double, double>(m_chiSquared, m_numberDoF);
 }
 
+void Vertex::setPosition(const Vector3& position) {
+  m_position.head<3>() = position;
+}
+
 void Vertex::setPosition(const Vector3& position, ActsScalar time) {
   m_position.head<3>() = position;
   m_position[eTime] = time;
diff --git a/Examples/Algorithms/MaterialMapping/include/ActsExamples/MaterialMapping/MaterialMapping.hpp b/Examples/Algorithms/MaterialMapping/include/ActsExamples/MaterialMapping/MaterialMapping.hpp
index 69e9cecb9fb..2f1e1e7334e 100644
--- a/Examples/Algorithms/MaterialMapping/include/ActsExamples/MaterialMapping/MaterialMapping.hpp
+++ b/Examples/Algorithms/MaterialMapping/include/ActsExamples/MaterialMapping/MaterialMapping.hpp
@@ -26,7 +26,6 @@
 #include <functional>
 #include <map>
 #include <memory>
-#include <mutex>
 #include <string>
 #include <unordered_map>
 #include <utility>
@@ -102,16 +101,15 @@ class MaterialMapping : public IAlgorithm {
   MaterialMapping(const Config& cfg,
                   Acts::Logging::Level level = Acts::Logging::INFO);
 
-  /// Destructor
-  /// - it also writes out the file
-  ~MaterialMapping() override;
-
   /// Framework execute method
   ///
   /// @param context The algorithm context for event consistency
   ActsExamples::ProcessCode execute(
       const AlgorithmContext& context) const override;
 
+  // Write out the file
+  ProcessCode finalize() override;
+
   /// Return the parameters to optimised the material map for a given surface
   /// Those parameters are the variance and the number of track for each bin
   ///
diff --git a/Examples/Algorithms/MaterialMapping/src/MaterialMapping.cpp b/Examples/Algorithms/MaterialMapping/src/MaterialMapping.cpp
index c673f36ba60..6e40e27c5b9 100644
--- a/Examples/Algorithms/MaterialMapping/src/MaterialMapping.cpp
+++ b/Examples/Algorithms/MaterialMapping/src/MaterialMapping.cpp
@@ -12,7 +12,6 @@
 #include "Acts/Material/AccumulatedSurfaceMaterial.hpp"
 #include "ActsExamples/MaterialMapping/IMaterialWriter.hpp"
 
-#include <iostream>
 #include <stdexcept>
 #include <unordered_map>
 
@@ -48,7 +47,8 @@ MaterialMapping::MaterialMapping(const MaterialMapping::Config& cfg,
   }
 }
 
-MaterialMapping::~MaterialMapping() {
+ProcessCode MaterialMapping::finalize() {
+  ACTS_INFO("Finalizing material mappig output");
   Acts::DetectorMaterialMaps detectorMaterial;
 
   if (m_cfg.materialSurfaceMapper && m_cfg.materialVolumeMapper) {
@@ -93,6 +93,8 @@ MaterialMapping::~MaterialMapping() {
   for (auto& imw : m_cfg.materialWriters) {
     imw->writeMaterial(detectorMaterial);
   }
+
+  return ProcessCode::SUCCESS;
 }
 
 ProcessCode MaterialMapping::execute(const AlgorithmContext& context) const {
diff --git a/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp
index 1b748e263f5..fb9e7d55560 100644
--- a/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp
+++ b/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp
@@ -28,7 +28,6 @@
 #include "Acts/Surfaces/PerigeeSurface.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp"
-#include "Acts/TrackFitting/GainMatrixSmoother.hpp"
 #include "Acts/TrackFitting/GainMatrixUpdater.hpp"
 #include "Acts/TrackFitting/KalmanFitter.hpp"
 #include "Acts/Utilities/Delegate.hpp"
diff --git a/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithmFunction.cpp b/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithmFunction.cpp
index 7dfd713bb75..1e91b33e217 100644
--- a/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithmFunction.cpp
+++ b/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithmFunction.cpp
@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2021 CERN for the benefit of the Acts project
+// Copyright (C) 2021-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -17,7 +17,6 @@
 #include "Acts/Propagator/Propagator.hpp"
 #include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp"
 #include "Acts/TrackFitting/GainMatrixSmoother.hpp"
-#include "Acts/TrackFitting/GainMatrixUpdater.hpp"
 #include "Acts/Utilities/Logger.hpp"
 #include "ActsExamples/EventData/Track.hpp"
 #include "ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp"
@@ -35,9 +34,6 @@ class TrackingGeometry;
 
 namespace {
 
-using Updater = Acts::GainMatrixUpdater;
-using Smoother = Acts::GainMatrixSmoother;
-
 using Stepper = Acts::EigenStepper<>;
 using Navigator = Acts::Navigator;
 using Propagator = Acts::Propagator<Stepper, Navigator>;
diff --git a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/RefittingCalibrator.hpp b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/RefittingCalibrator.hpp
index a4f2b6b6f2c..d9a0c80d553 100644
--- a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/RefittingCalibrator.hpp
+++ b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/RefittingCalibrator.hpp
@@ -29,12 +29,14 @@ struct RefittingCalibrator {
 
   struct RefittingSourceLink {
     ConstProxy state;
-
-    Acts::GeometryIdentifier geometryId() const {
-      return state.referenceSurface().geometryId();
-    }
   };
 
+  static const Acts::Surface* accessSurface(
+      const Acts::SourceLink& sourceLink) {
+    const auto& refittingSl = sourceLink.get<RefittingSourceLink>();
+    return &refittingSl.state.referenceSurface();
+  }
+
   void calibrate(const Acts::GeometryContext& gctx,
                  const Acts::CalibrationContext& cctx,
                  const Acts::SourceLink& sourceLink, Proxy trackState) const;
diff --git a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/TrackFitterFunction.hpp b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/TrackFitterFunction.hpp
index 3785be1d0c0..8e00a47b4a7 100644
--- a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/TrackFitterFunction.hpp
+++ b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/TrackFitterFunction.hpp
@@ -40,6 +40,7 @@ class TrackFitterFunction {
     std::reference_wrapper<const Acts::CalibrationContext> calibrationContext;
     const Acts::Surface* referenceSurface = nullptr;
     Acts::PropagatorPlainOptions propOptions;
+    bool doRefit = false;
   };
 
   virtual ~TrackFitterFunction() = default;
diff --git a/Examples/Algorithms/TrackFitting/src/GlobalChiSquareFitterFunction.cpp b/Examples/Algorithms/TrackFitting/src/GlobalChiSquareFitterFunction.cpp
index d6e32055667..47ff781da06 100644
--- a/Examples/Algorithms/TrackFitting/src/GlobalChiSquareFitterFunction.cpp
+++ b/Examples/Algorithms/TrackFitting/src/GlobalChiSquareFitterFunction.cpp
@@ -86,9 +86,13 @@ struct GlobalChiSquareFitterFunctionImpl final : public TrackFitterFunction {
         extensions;
     extensions.calibrator.connect<&calibrator_t::calibrate>(&calibrator);
 
-    extensions.surfaceAccessor
-        .connect<&IndexSourceLink::SurfaceAccessor::operator()>(
-            &m_slSurfaceAccessor);
+    if (options.doRefit) {
+      extensions.surfaceAccessor.connect<&RefittingCalibrator::accessSurface>();
+    } else {
+      extensions.surfaceAccessor
+          .connect<&IndexSourceLink::SurfaceAccessor::operator()>(
+              &m_slSurfaceAccessor);
+    }
 
     const Acts::Experimental::Gx2FitterOptions gx2fOptions(
         options.geoContext, options.magFieldContext, options.calibrationContext,
diff --git a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp
index e92a1f1072f..ca0af32d204 100644
--- a/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp
+++ b/Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp
@@ -116,9 +116,15 @@ struct GsfFitterFunctionImpl final : public ActsExamples::TrackFitterFunction {
 
     gsfOptions.extensions.calibrator.connect<&calibrator_t::calibrate>(
         &calibrator);
-    gsfOptions.extensions.surfaceAccessor
-        .connect<&IndexSourceLink::SurfaceAccessor::operator()>(
-            &m_slSurfaceAccessor);
+
+    if (options.doRefit) {
+      gsfOptions.extensions.surfaceAccessor
+          .connect<&RefittingCalibrator::accessSurface>();
+    } else {
+      gsfOptions.extensions.surfaceAccessor
+          .connect<&IndexSourceLink::SurfaceAccessor::operator()>(
+              &m_slSurfaceAccessor);
+    }
     switch (reductionAlg) {
       case MixtureReductionAlgorithm::weightCut: {
         gsfOptions.extensions.mixtureReducer
diff --git a/Examples/Algorithms/TrackFitting/src/KalmanFitterFunction.cpp b/Examples/Algorithms/TrackFitting/src/KalmanFitterFunction.cpp
index b5425862285..58a0d23aef6 100644
--- a/Examples/Algorithms/TrackFitting/src/KalmanFitterFunction.cpp
+++ b/Examples/Algorithms/TrackFitting/src/KalmanFitterFunction.cpp
@@ -114,9 +114,15 @@ struct KalmanFitterFunctionImpl final : public TrackFitterFunction {
     kfOptions.freeToBoundCorrection = freeToBoundCorrection;
     kfOptions.extensions.calibrator.connect<&calibrator_t::calibrate>(
         &calibrator);
-    kfOptions.extensions.surfaceAccessor
-        .connect<&IndexSourceLink::SurfaceAccessor::operator()>(
-            &slSurfaceAccessor);
+
+    if (options.doRefit) {
+      kfOptions.extensions.surfaceAccessor
+          .connect<&RefittingCalibrator::accessSurface>();
+    } else {
+      kfOptions.extensions.surfaceAccessor
+          .connect<&IndexSourceLink::SurfaceAccessor::operator()>(
+              &slSurfaceAccessor);
+    }
 
     return kfOptions;
   }
diff --git a/Examples/Algorithms/TrackFitting/src/RefittingAlgorithm.cpp b/Examples/Algorithms/TrackFitting/src/RefittingAlgorithm.cpp
index e4642b23c3d..8b86fb8aea9 100644
--- a/Examples/Algorithms/TrackFitting/src/RefittingAlgorithm.cpp
+++ b/Examples/Algorithms/TrackFitting/src/RefittingAlgorithm.cpp
@@ -71,6 +71,7 @@ ActsExamples::ProcessCode ActsExamples::RefittingAlgorithm::execute(
         ctx.geoContext, ctx.magFieldContext, ctx.calibContext,
         &track.referenceSurface(),
         Acts::PropagatorPlainOptions(ctx.geoContext, ctx.magFieldContext)};
+    options.doRefit = true;
 
     const Acts::BoundTrackParameters initialParams(
         track.referenceSurface().getSharedPtr(), track.parameters(),
diff --git a/Examples/Algorithms/Vertexing/CMakeLists.txt b/Examples/Algorithms/Vertexing/CMakeLists.txt
index f35ce7f559d..b982e10bc26 100644
--- a/Examples/Algorithms/Vertexing/CMakeLists.txt
+++ b/Examples/Algorithms/Vertexing/CMakeLists.txt
@@ -4,6 +4,7 @@ add_library(
   src/IterativeVertexFinderAlgorithm.cpp
   src/VertexFitterAlgorithm.cpp
   src/SingleSeedVertexFinderAlgorithm.cpp
+  src/TruthVertexSeeder.cpp
 )
 target_include_directories(
   ActsExamplesVertexing
diff --git a/Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/AdaptiveMultiVertexFinderAlgorithm.hpp b/Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/AdaptiveMultiVertexFinderAlgorithm.hpp
index 697a646ed50..735edcd1b0c 100644
--- a/Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/AdaptiveMultiVertexFinderAlgorithm.hpp
+++ b/Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/AdaptiveMultiVertexFinderAlgorithm.hpp
@@ -22,10 +22,13 @@
 #include "Acts/Vertexing/AdaptiveMultiVertexFitter.hpp"
 #include "Acts/Vertexing/GaussianTrackDensity.hpp"
 #include "Acts/Vertexing/HelicalTrackLinearizer.hpp"
+#include "Acts/Vertexing/IVertexFinder.hpp"
 #include "Acts/Vertexing/ImpactPointEstimator.hpp"
 #include "Acts/Vertexing/TrackDensityVertexFinder.hpp"
 #include "Acts/Vertexing/Vertex.hpp"
 #include "ActsExamples/EventData/ProtoVertex.hpp"
+#include "ActsExamples/EventData/SimParticle.hpp"
+#include "ActsExamples/EventData/SimVertex.hpp"
 #include "ActsExamples/EventData/Track.hpp"
 #include "ActsExamples/EventData/Trajectories.hpp"
 #include "ActsExamples/Framework/DataHandle.hpp"
@@ -58,11 +61,16 @@ class AdaptiveMultiVertexFinderAlgorithm final : public IAlgorithm {
 
   using VertexCollection = std::vector<Acts::Vertex>;
 
-  enum class SeedFinder { GaussianSeeder, AdaptiveGridSeeder };
+  enum class SeedFinder { TruthSeeder, GaussianSeeder, AdaptiveGridSeeder };
 
   struct Config {
     /// Input track parameters collection
     std::string inputTrackParameters;
+    /// Optional
+    std::string inputTruthParticles;
+    /// Optional: Input truth vertex collection. This will only be used if
+    /// `seedFinder == SeedFinder::TruthSeeder`.
+    std::string inputTruthVertices;
     /// Output proto vertex collection
     std::string outputProtoVertices;
     /// Output vertex collection
@@ -71,10 +79,31 @@ class AdaptiveMultiVertexFinderAlgorithm final : public IAlgorithm {
     /// The magnetic field
     std::shared_ptr<Acts::MagneticFieldProvider> bField;
 
-    /// Enum member determining the choice of the vertex seed finder
-    SeedFinder seedFinder;
+    /// For more information look at `AdaptiveMultiVertexFitter.hpp`
+    Acts::AnnealingUtility::Config annealingConfig{9., {1.0}};
+
+    /// For more information look at `AdaptiveMultiVertexFitter.hpp`
+    double minWeight = 0.001;
+    /// For more information look at `AdaptiveMultiVertexFitter.hpp`
+    bool doSmoothing = true;
+
+    /// Maximum number of iterations for the vertex finding
+    int maxIterations = 1000;
     /// Use time information in vertex seeder, finder, and fitter
     bool useTime = false;
+    /// For more information look at `AdaptiveMultiVertexFinder.hpp`
+    double tracksMaxZinterval = 1. * Acts::UnitConstants::mm;
+    /// For more information look at `AdaptiveMultiVertexFinder.hpp`
+    Acts::Vector4 initialVariances = Acts::Vector4{1e+2, 1e+2, 1e+2, 1e+8};
+    /// For more information look at `AdaptiveMultiVertexFinder.hpp`
+    bool doFullSplitting = false;
+    /// For more information look at `AdaptiveMultiVertexFinder.hpp`
+    std::optional<double> tracksMaxSignificance;
+    /// For more information look at `AdaptiveMultiVertexFinder.hpp`
+    std::optional<double> maxMergeVertexSignificance;
+
+    /// Enum member determining the choice of the vertex seed finder
+    SeedFinder seedFinder;
     /// Bin extent in z-direction which is only used with `AdaptiveGridSeeder`
     double spatialBinExtent = 15. * Acts::UnitConstants::um;
     /// Bin extent in t-direction which is only used with `AdaptiveGridSeeder`
@@ -96,21 +125,26 @@ class AdaptiveMultiVertexFinderAlgorithm final : public IAlgorithm {
   const Config& config() const { return m_cfg; }
 
  private:
-  Acts::AdaptiveMultiVertexFinder makeVertexFinder() const;
+  std::unique_ptr<Acts::IVertexFinder> makeVertexSeeder() const;
+  Acts::AdaptiveMultiVertexFinder makeVertexFinder(
+      std::shared_ptr<const Acts::IVertexFinder> seedFinder) const;
 
   Config m_cfg;
 
   std::shared_ptr<const Acts::BasePropagator> m_propagator;
   Acts::ImpactPointEstimator m_ipEstimator;
   Linearizer m_linearizer;
+  std::shared_ptr<Acts::IVertexFinder> m_vertexSeeder;
   Acts::AdaptiveMultiVertexFinder m_vertexFinder;
 
   ReadDataHandle<TrackParametersContainer> m_inputTrackParameters{
       this, "InputTrackParameters"};
-
+  ReadDataHandle<SimParticleContainer> m_inputTruthParticles{
+      this, "InputTruthParticles"};
+  ReadDataHandle<SimVertexContainer> m_inputTruthVertices{this,
+                                                          "InputTruthVertices"};
   WriteDataHandle<ProtoVertexContainer> m_outputProtoVertices{
       this, "OutputProtoVertices"};
-
   WriteDataHandle<VertexCollection> m_outputVertices{this, "OutputVertices"};
 };
 
diff --git a/Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/IterativeVertexFinderAlgorithm.hpp b/Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/IterativeVertexFinderAlgorithm.hpp
index 384c16051db..15cbe4cbecf 100644
--- a/Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/IterativeVertexFinderAlgorithm.hpp
+++ b/Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/IterativeVertexFinderAlgorithm.hpp
@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2021 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -67,8 +67,12 @@ class IterativeVertexFinderAlgorithm final : public IAlgorithm {
     std::string outputProtoVertices;
     /// Output vertex collection
     std::string outputVertices = "vertices";
+
     /// The magnetic field
     std::shared_ptr<Acts::MagneticFieldProvider> bField;
+
+    /// Maximum number of iterations for the vertex finding
+    int maxIterations = 1000;
   };
 
   IterativeVertexFinderAlgorithm(const Config& config,
diff --git a/Examples/Algorithms/Vertexing/src/AdaptiveMultiVertexFinderAlgorithm.cpp b/Examples/Algorithms/Vertexing/src/AdaptiveMultiVertexFinderAlgorithm.cpp
index 91a8e5c81c3..ec933e739f1 100644
--- a/Examples/Algorithms/Vertexing/src/AdaptiveMultiVertexFinderAlgorithm.cpp
+++ b/Examples/Algorithms/Vertexing/src/AdaptiveMultiVertexFinderAlgorithm.cpp
@@ -25,21 +25,28 @@
 #include "Acts/Vertexing/Vertex.hpp"
 #include "Acts/Vertexing/VertexingOptions.hpp"
 #include "ActsExamples/EventData/ProtoVertex.hpp"
+#include "ActsExamples/EventData/SimParticle.hpp"
+#include "ActsExamples/EventData/SimVertex.hpp"
 #include "ActsExamples/Framework/AlgorithmContext.hpp"
 #include "ActsExamples/Framework/ProcessCode.hpp"
+#include "ActsExamples/TruthTracking/TruthVertexFinder.hpp"
 
 #include <memory>
 #include <optional>
 #include <ostream>
 #include <stdexcept>
 #include <system_error>
+#include <unordered_map>
+#include <utility>
 
+#include "TruthVertexSeeder.hpp"
 #include "VertexingHelpers.hpp"
 
-ActsExamples::AdaptiveMultiVertexFinderAlgorithm::
-    AdaptiveMultiVertexFinderAlgorithm(const Config& config,
-                                       Acts::Logging::Level level)
-    : ActsExamples::IAlgorithm("AdaptiveMultiVertexFinder", level),
+namespace ActsExamples {
+
+AdaptiveMultiVertexFinderAlgorithm::AdaptiveMultiVertexFinderAlgorithm(
+    const Config& config, Acts::Logging::Level level)
+    : IAlgorithm("AdaptiveMultiVertexFinder", level),
       m_cfg(config),
       m_propagator{[&]() {
         // Set up EigenStepper
@@ -63,7 +70,8 @@ ActsExamples::AdaptiveMultiVertexFinderAlgorithm::
         return Linearizer(ltConfig,
                           logger().cloneWithSuffix("HelicalTrackLinearizer"));
       }()},
-      m_vertexFinder{makeVertexFinder()} {
+      m_vertexSeeder{makeVertexSeeder()},
+      m_vertexFinder{makeVertexFinder(m_vertexSeeder)} {
   if (m_cfg.inputTrackParameters.empty()) {
     throw std::invalid_argument("Missing input track parameter collection");
   }
@@ -73,22 +81,46 @@ ActsExamples::AdaptiveMultiVertexFinderAlgorithm::
   if (m_cfg.outputVertices.empty()) {
     throw std::invalid_argument("Missing output vertices collection");
   }
+  if (m_cfg.seedFinder == SeedFinder::TruthSeeder &&
+      m_cfg.inputTruthVertices.empty()) {
+    throw std::invalid_argument("Missing input truth vertex collection");
+  }
+
+  // Sanitize the configuration
+  if (m_cfg.seedFinder != SeedFinder::TruthSeeder &&
+      (!m_cfg.inputTruthParticles.empty() ||
+       !m_cfg.inputTruthVertices.empty())) {
+    ACTS_INFO("Ignoring truth input as seed finder is not TruthSeeder");
+    m_cfg.inputTruthVertices.clear();
+    m_cfg.inputTruthVertices.clear();
+  }
 
   m_inputTrackParameters.initialize(m_cfg.inputTrackParameters);
+  m_inputTruthParticles.maybeInitialize(m_cfg.inputTruthParticles);
+  m_inputTruthVertices.maybeInitialize(m_cfg.inputTruthVertices);
   m_outputProtoVertices.initialize(m_cfg.outputProtoVertices);
   m_outputVertices.initialize(m_cfg.outputVertices);
 }
 
-auto ActsExamples::AdaptiveMultiVertexFinderAlgorithm::makeVertexFinder() const
-    -> Acts::AdaptiveMultiVertexFinder {
-  std::shared_ptr<const Acts::IVertexFinder> seedFinder;
+std::unique_ptr<Acts::IVertexFinder>
+AdaptiveMultiVertexFinderAlgorithm::makeVertexSeeder() const {
+  if (m_cfg.seedFinder == SeedFinder::TruthSeeder) {
+    using Seeder = ActsExamples::TruthVertexSeeder;
+    Seeder::Config seederConfig;
+    seederConfig.useXY = false;
+    seederConfig.useTime = m_cfg.useTime;
+    return std::make_unique<Seeder>(seederConfig);
+  }
+
   if (m_cfg.seedFinder == SeedFinder::GaussianSeeder) {
     using Seeder = Acts::TrackDensityVertexFinder;
     Acts::GaussianTrackDensity::Config trkDensityCfg;
     trkDensityCfg.extractParameters
         .connect<&Acts::InputTrack::extractParameters>();
-    seedFinder = std::make_shared<Seeder>(Seeder::Config{trkDensityCfg});
-  } else if (m_cfg.seedFinder == SeedFinder::AdaptiveGridSeeder) {
+    return std::make_unique<Seeder>(Seeder::Config{trkDensityCfg});
+  }
+
+  if (m_cfg.seedFinder == SeedFinder::AdaptiveGridSeeder) {
     // Set up track density used during vertex seeding
     Acts::AdaptiveGridTrackDensity::Config trkDensityCfg;
     // Bin extent in z-direction
@@ -103,21 +135,23 @@ auto ActsExamples::AdaptiveMultiVertexFinderAlgorithm::makeVertexFinder() const
     Seeder::Config seederConfig(trkDensity);
     seederConfig.extractParameters
         .connect<&Acts::InputTrack::extractParameters>();
-    seedFinder = std::make_shared<Seeder>(seederConfig);
-  } else {
-    throw std::invalid_argument("Unknown seed finder");
+    return std::make_unique<Seeder>(seederConfig);
   }
 
+  throw std::invalid_argument("Unknown seed finder");
+}
+
+Acts::AdaptiveMultiVertexFinder
+AdaptiveMultiVertexFinderAlgorithm::makeVertexFinder(
+    std::shared_ptr<const Acts::IVertexFinder> seedFinder) const {
   // Set up deterministic annealing with user-defined temperatures
-  Acts::AnnealingUtility::Config annealingConfig;
-  annealingConfig.setOfTemperatures = {1.};
-  Acts::AnnealingUtility annealingUtility(annealingConfig);
+  Acts::AnnealingUtility annealingUtility(m_cfg.annealingConfig);
 
   // Set up the vertex fitter with user-defined annealing
   Fitter::Config fitterCfg(m_ipEstimator);
   fitterCfg.annealingTool = annealingUtility;
-  fitterCfg.minWeight = 0.001;
-  fitterCfg.doSmoothing = true;
+  fitterCfg.minWeight = m_cfg.minWeight;
+  fitterCfg.doSmoothing = m_cfg.doSmoothing;
   fitterCfg.useTime = m_cfg.useTime;
   fitterCfg.extractParameters.connect<&Acts::InputTrack::extractParameters>();
   fitterCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&m_linearizer);
@@ -125,18 +159,18 @@ auto ActsExamples::AdaptiveMultiVertexFinderAlgorithm::makeVertexFinder() const
                 logger().cloneWithSuffix("AdaptiveMultiVertexFitter"));
 
   Acts::AdaptiveMultiVertexFinder::Config finderConfig(
-      std::move(fitter), seedFinder, m_ipEstimator, m_cfg.bField);
+      std::move(fitter), std::move(seedFinder), m_ipEstimator, m_cfg.bField);
   // Set the initial variance of the 4D vertex position. Since time is on a
   // numerical scale, we have to provide a greater value in the corresponding
   // dimension.
-  finderConfig.initialVariances << 1e+2, 1e+2, 1e+2, 1e+8;
-  finderConfig.tracksMaxZinterval = 1. * Acts::UnitConstants::mm;
-  finderConfig.maxIterations = 200;
+  finderConfig.initialVariances = m_cfg.initialVariances;
+  finderConfig.tracksMaxZinterval = m_cfg.tracksMaxZinterval;
+  finderConfig.maxIterations = m_cfg.maxIterations;
   finderConfig.useTime = m_cfg.useTime;
   // 5 corresponds to a p-value of ~0.92 using `chi2(x=5,ndf=2)`
   finderConfig.tracksMaxSignificance = 5;
   // This should be used consistently with and without time
-  finderConfig.doFullSplitting = false;
+  finderConfig.doFullSplitting = m_cfg.doFullSplitting;
   // 3 corresponds to a p-value of ~0.92 using `chi2(x=3,ndf=1)`
   finderConfig.maxMergeVertexSignificance = 3;
   if (m_cfg.useTime) {
@@ -150,21 +184,29 @@ auto ActsExamples::AdaptiveMultiVertexFinderAlgorithm::makeVertexFinder() const
     // 5 corresponds to a p-value of ~0.92 using `chi2(x=5,ndf=2)`
     finderConfig.maxMergeVertexSignificance = 5;
   }
+
   finderConfig.extractParameters
       .template connect<&Acts::InputTrack::extractParameters>();
 
+  if (m_cfg.seedFinder == SeedFinder::TruthSeeder) {
+    finderConfig.doNotBreakWhileSeeding = true;
+  }
+
+  finderConfig.tracksMaxSignificance =
+      m_cfg.tracksMaxSignificance.value_or(finderConfig.tracksMaxSignificance);
+  finderConfig.maxMergeVertexSignificance =
+      m_cfg.maxMergeVertexSignificance.value_or(
+          finderConfig.maxMergeVertexSignificance);
+
   // Instantiate the finder
   return Acts::AdaptiveMultiVertexFinder(std::move(finderConfig),
                                          logger().clone());
 }
 
-ActsExamples::ProcessCode
-ActsExamples::AdaptiveMultiVertexFinderAlgorithm::execute(
-    const ActsExamples::AlgorithmContext& ctx) const {
-  // retrieve input tracks and convert into the expected format
-
+ProcessCode AdaptiveMultiVertexFinderAlgorithm::execute(
+    const AlgorithmContext& ctx) const {
   const auto& inputTrackParameters = m_inputTrackParameters(ctx);
-  // TODO change this from pointers to tracks parameters to actual tracks
+
   auto inputTracks = makeInputTracks(inputTrackParameters);
 
   if (inputTrackParameters.size() != inputTracks.size()) {
@@ -181,13 +223,50 @@ ActsExamples::AdaptiveMultiVertexFinderAlgorithm::execute(
     }
   }
 
-  //////////////////////////////////////////////
-  /* Full tutorial example code for reference */
-  //////////////////////////////////////////////
-
   // The vertex finder state
   auto state = m_vertexFinder.makeState(ctx.magFieldContext);
 
+  // In case of the truth seeder, we need to wire the truth vertices into the
+  // vertex finder
+  if (m_cfg.seedFinder == SeedFinder::TruthSeeder) {
+    const auto& truthParticles = m_inputTruthParticles(ctx);
+    const auto& truthVertices = m_inputTruthVertices(ctx);
+
+    auto& vertexSeederState =
+        state.as<Acts::AdaptiveMultiVertexFinder::State>()
+            .seedFinderState.as<TruthVertexSeeder::State>();
+
+    std::map<SimVertexBarcode, std::size_t> vertexParticleCount;
+
+    for (const auto& truthVertex : truthVertices) {
+      // Skip secondary vertices
+      if (truthVertex.vertexId().vertexSecondary() != 0) {
+        continue;
+      }
+      vertexSeederState.truthVertices.push_back(truthVertex);
+
+      // Count the number of particles associated with each vertex
+      std::size_t particleCount = 0;
+      for (const auto& particle : truthParticles) {
+        if (particle.particleId().vertexId() == truthVertex.vertexId()) {
+          ++particleCount;
+        }
+      }
+      vertexParticleCount[truthVertex.vertexId()] = particleCount;
+    }
+
+    // sort by number of particles
+    std::sort(vertexSeederState.truthVertices.begin(),
+              vertexSeederState.truthVertices.end(),
+              [&vertexParticleCount](const auto& lhs, const auto& rhs) {
+                return vertexParticleCount[lhs.vertexId()] >
+                       vertexParticleCount[rhs.vertexId()];
+              });
+
+    ACTS_INFO("Got " << truthVertices.size() << " truth vertices and selected "
+                     << vertexSeederState.truthVertices.size() << " in event");
+  }
+
   // Default vertexing options, this is where e.g. a constraint could be set
   Options finderOpts(ctx.geoContext, ctx.magFieldContext);
 
@@ -221,5 +300,7 @@ ActsExamples::AdaptiveMultiVertexFinderAlgorithm::execute(
   // store found vertices
   m_outputVertices(ctx, std::move(vertices));
 
-  return ActsExamples::ProcessCode::SUCCESS;
+  return ProcessCode::SUCCESS;
 }
+
+}  // namespace ActsExamples
diff --git a/Examples/Algorithms/Vertexing/src/IterativeVertexFinderAlgorithm.cpp b/Examples/Algorithms/Vertexing/src/IterativeVertexFinderAlgorithm.cpp
index a812b249625..237a68010a9 100644
--- a/Examples/Algorithms/Vertexing/src/IterativeVertexFinderAlgorithm.cpp
+++ b/Examples/Algorithms/Vertexing/src/IterativeVertexFinderAlgorithm.cpp
@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2019-2021 CERN for the benefit of the Acts project
+// Copyright (C) 2019-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -26,9 +26,11 @@
 
 #include "VertexingHelpers.hpp"
 
-ActsExamples::IterativeVertexFinderAlgorithm::IterativeVertexFinderAlgorithm(
+namespace ActsExamples {
+
+IterativeVertexFinderAlgorithm::IterativeVertexFinderAlgorithm(
     const Config& config, Acts::Logging::Level level)
-    : ActsExamples::IAlgorithm("IterativeVertexFinder", level), m_cfg(config) {
+    : IAlgorithm("IterativeVertexFinder", level), m_cfg(config) {
   if (m_cfg.inputTrackParameters.empty()) {
     throw std::invalid_argument("Missing input track parameter collection");
   }
@@ -44,8 +46,8 @@ ActsExamples::IterativeVertexFinderAlgorithm::IterativeVertexFinderAlgorithm(
   m_outputVertices.initialize(m_cfg.outputVertices);
 }
 
-ActsExamples::ProcessCode ActsExamples::IterativeVertexFinderAlgorithm::execute(
-    const ActsExamples::AlgorithmContext& ctx) const {
+ProcessCode IterativeVertexFinderAlgorithm::execute(
+    const AlgorithmContext& ctx) const {
   // retrieve input tracks and convert into the expected format
 
   const auto& inputTrackParameters = m_inputTrackParameters(ctx);
@@ -100,7 +102,7 @@ ActsExamples::ProcessCode ActsExamples::IterativeVertexFinderAlgorithm::execute(
   Finder::Config finderCfg(std::move(vertexFitter), seeder, ipEst);
   finderCfg.trackLinearizer.connect<&Linearizer::linearizeTrack>(&linearizer);
 
-  finderCfg.maxVertices = 200;
+  finderCfg.maxVertices = m_cfg.maxIterations;
   finderCfg.reassignTracksAfterFirstFit = false;
   finderCfg.extractParameters.connect<&Acts::InputTrack::extractParameters>();
   finderCfg.field = m_cfg.bField;
@@ -132,5 +134,7 @@ ActsExamples::ProcessCode ActsExamples::IterativeVertexFinderAlgorithm::execute(
   // store found vertices
   m_outputVertices(ctx, std::move(vertices));
 
-  return ActsExamples::ProcessCode::SUCCESS;
+  return ProcessCode::SUCCESS;
 }
+
+}  // namespace ActsExamples
diff --git a/Examples/Algorithms/Vertexing/src/TruthVertexSeeder.cpp b/Examples/Algorithms/Vertexing/src/TruthVertexSeeder.cpp
new file mode 100644
index 00000000000..89a3ba90d91
--- /dev/null
+++ b/Examples/Algorithms/Vertexing/src/TruthVertexSeeder.cpp
@@ -0,0 +1,57 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2024 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "TruthVertexSeeder.hpp"
+
+#include "VertexingHelpers.hpp"
+
+namespace ActsExamples {
+
+TruthVertexSeeder::TruthVertexSeeder(const Config &cfg) : m_cfg(cfg) {}
+
+Acts::Result<std::vector<Acts::Vertex>> TruthVertexSeeder::find(
+    const std::vector<Acts::InputTrack> & /*trackVector*/,
+    const Acts::VertexingOptions &vertexingOptions,
+    Acts::IVertexFinder::State &anyState) const {
+  auto &state = anyState.template as<State>();
+
+  if (state.nextVertexIndex >= state.truthVertices.size()) {
+    return std::vector<Acts::Vertex>();
+  }
+
+  const auto &truthVertex = state.truthVertices[state.nextVertexIndex];
+  ++state.nextVertexIndex;
+
+  Acts::Vertex converted;
+  converted.fullPosition().z() = truthVertex.position().z();
+  if (m_cfg.useXY) {
+    converted.fullPosition().x() = truthVertex.position().x();
+    converted.fullPosition().y() = truthVertex.position().y();
+  }
+  if (m_cfg.useTime) {
+    converted.setTime(truthVertex.time());
+  }
+
+  Acts::SquareMatrix4 seedCov = vertexingOptions.constraint.fullCovariance();
+  converted.setFullCovariance(seedCov);
+
+  return std::vector<Acts::Vertex>{converted};
+}
+
+Acts::IVertexFinder::State TruthVertexSeeder::makeState(
+    const Acts::MagneticFieldContext & /*mctx*/) const {
+  return Acts::IVertexFinder::State{State{}};
+}
+
+void TruthVertexSeeder::setTracksToRemove(
+    Acts::IVertexFinder::State & /*anyState*/,
+    const std::vector<Acts::InputTrack> & /*removedTracks*/) const {
+  // nothing to do
+}
+
+}  // namespace ActsExamples
diff --git a/Examples/Algorithms/Vertexing/src/TruthVertexSeeder.hpp b/Examples/Algorithms/Vertexing/src/TruthVertexSeeder.hpp
new file mode 100644
index 00000000000..27a27425c4b
--- /dev/null
+++ b/Examples/Algorithms/Vertexing/src/TruthVertexSeeder.hpp
@@ -0,0 +1,51 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2024 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/MagneticField/MagneticFieldContext.hpp"
+#include "Acts/Vertexing/IVertexFinder.hpp"
+#include "Acts/Vertexing/TrackAtVertex.hpp"
+#include "Acts/Vertexing/Vertex.hpp"
+#include "Acts/Vertexing/VertexingOptions.hpp"
+#include "ActsExamples/EventData/SimVertex.hpp"
+
+namespace ActsExamples {
+
+class TruthVertexSeeder final : public Acts::IVertexFinder {
+ public:
+  struct Config {
+    bool useXY = false;
+    bool useTime = false;
+  };
+
+  struct State {
+    std::vector<SimVertex> truthVertices;
+
+    std::size_t nextVertexIndex = 0;
+  };
+
+  explicit TruthVertexSeeder(const Config& cfg);
+
+  Acts::Result<std::vector<Acts::Vertex>> find(
+      const std::vector<Acts::InputTrack>& trackVector,
+      const Acts::VertexingOptions& vertexingOptions,
+      IVertexFinder::State& state) const final;
+
+  IVertexFinder::State makeState(
+      const Acts::MagneticFieldContext& mctx) const final;
+
+  void setTracksToRemove(
+      IVertexFinder::State& anyState,
+      const std::vector<Acts::InputTrack>& removedTracks) const final;
+
+ private:
+  Config m_cfg;
+};
+
+}  // namespace ActsExamples
diff --git a/Examples/Detectors/DD4hepDetector/CMakeLists.txt b/Examples/Detectors/DD4hepDetector/CMakeLists.txt
index bafd0878276..88192fdbd4e 100644
--- a/Examples/Detectors/DD4hepDetector/CMakeLists.txt
+++ b/Examples/Detectors/DD4hepDetector/CMakeLists.txt
@@ -1,8 +1,7 @@
 add_library(
   ActsExamplesDetectorDD4hep SHARED
   src/DD4hepDetector.cpp
-  src/DD4hepGeometryService.cpp
-  src/DetUtils.cpp)
+  src/DD4hepGeometryService.cpp)
 
 
 target_include_directories(
diff --git a/Examples/Detectors/DD4hepDetector/src/DetUtils.cpp b/Examples/Detectors/DD4hepDetector/src/DetUtils.cpp
deleted file mode 100644
index 007a539be2a..00000000000
--- a/Examples/Detectors/DD4hepDetector/src/DetUtils.cpp
+++ /dev/null
@@ -1,40 +0,0 @@
-// This file is part of the Acts project.
-//
-// Copyright (C) 2017 CERN for the benefit of the Acts project
-//
-// This Source Code Form is subject to the terms of the Mozilla Public
-// License, v. 2.0. If a copy of the MPL was not distributed with this
-// file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#include "DetUtils.h"
-
-#include <cstddef>
-
-#include <XML/XMLTags.h>
-
-namespace det::utils {
-
-dd4hep::xml::Component getNodeByStrAttr(const dd4hep::xml::Handle_t& mother,
-                                        const std::string& nodeName,
-                                        const std::string& attrName,
-                                        const std::string& attrValue) {
-  for (dd4hep::xml::Collection_t xCompColl(mother, nodeName.c_str());
-       nullptr != xCompColl; ++xCompColl) {
-    if (xCompColl.attr<std::string>(attrName.c_str()) == attrValue) {
-      return static_cast<dd4hep::xml::Component>(xCompColl);
-    }
-  }
-  // in case there was no xml daughter with matching name
-  return dd4hep::xml::Component(nullptr);
-}
-
-double getAttrValueWithFallback(const dd4hep::xml::Component& node,
-                                const std::string& attrName,
-                                const double& defaultValue) {
-  if (node.hasAttr(_Unicode(attrName.c_str()))) {
-    return node.attr<double>(attrName.c_str());
-  } else {
-    return defaultValue;
-  }
-}
-}  // namespace det::utils
diff --git a/Examples/Detectors/DD4hepDetector/src/DetUtils.h b/Examples/Detectors/DD4hepDetector/src/DetUtils.h
deleted file mode 100644
index b5418c19801..00000000000
--- a/Examples/Detectors/DD4hepDetector/src/DetUtils.h
+++ /dev/null
@@ -1,45 +0,0 @@
-#pragma once
-
-#include <memory>
-#include <string>
-
-// DD4hep
-#include "DD4hep/DetFactoryHelper.h"
-#include "DD4hep/Detector.h"
-#include "DD4hep/Objects.h"
-#include "DD4hep/Segmentations.h"
-#include "DDSegmentation/BitField64.h"
-#include "DDSegmentation/CartesianGridXY.h"
-#include "DDSegmentation/CartesianGridXYZ.h"
-#include "DDSegmentation/PolarGridRPhi.h"
-#include "XML/XMLDetector.h"
-#include "XML/XMLElements.h"
-
-namespace Acts {
-class DigitizationModule;
-}  // namespace Acts
-
-/** Given a xml element with several daughters with the same name, e.g.
- <detector> <layer name="1" /> <layer name="2"> </detector>
- this method returns the first daughter of type nodeName whose attribute has a
- given value
- e.g. returns <layer name="2"/> when called with (detector, "layer", "name",
- "1") */
-namespace det::utils {
-
-/// @brief Retrieves the node component from a mother by the string names
-/// @param mother The Handle to the mother volume
-/// @param nodeName The name of the note
-/// @param attrName The name of the Attribute
-/// @param attrValue The attribute value
-dd4hep::xml::Component getNodeByStrAttr(const dd4hep::xml::Handle_t& mother,
-                                        const std::string& nodeName,
-                                        const std::string& attrName,
-                                        const std::string& attrValue);
-
-/// try to get attribute with double value, return defaultValue if attribute
-/// not found
-double getAttrValueWithFallback(const dd4hep::xml::Component& node,
-                                const std::string& attrName,
-                                const double& defaultValue);
-}  // namespace det::utils
diff --git a/Examples/Framework/include/ActsExamples/Framework/DataHandle.hpp b/Examples/Framework/include/ActsExamples/Framework/DataHandle.hpp
index 044d0fe665d..6e8edbf07fc 100644
--- a/Examples/Framework/include/ActsExamples/Framework/DataHandle.hpp
+++ b/Examples/Framework/include/ActsExamples/Framework/DataHandle.hpp
@@ -66,7 +66,7 @@ class WriteDataHandle final : public DataHandleBase {
   }
 
   void operator()(const AlgorithmContext& ctx, T&& value) const {
-    (*this)(ctx.eventStore, std::forward<T>(value));
+    (*this)(ctx.eventStore, std::move(value));
   }
 
   void operator()(WhiteBoard& wb, T&& value) const {
diff --git a/Examples/Io/Root/include/ActsExamples/Io/Root/RootMeasurementWriter.hpp b/Examples/Io/Root/include/ActsExamples/Io/Root/RootMeasurementWriter.hpp
index a4205e13001..ccf1be2a901 100644
--- a/Examples/Io/Root/include/ActsExamples/Io/Root/RootMeasurementWriter.hpp
+++ b/Examples/Io/Root/include/ActsExamples/Io/Root/RootMeasurementWriter.hpp
@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2020 CERN for the benefit of the Acts project
+// Copyright (C) 2020-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -8,13 +8,10 @@
 
 #pragma once
 
-#include "Acts/Definitions/Algebra.hpp"
-#include "Acts/Definitions/Common.hpp"
 #include "Acts/Definitions/TrackParametrization.hpp"
 #include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/Geometry/GeometryHierarchyMap.hpp"
 #include "Acts/Geometry/GeometryIdentifier.hpp"
-#include "Acts/Geometry/TrackingGeometry.hpp"
 #include "Acts/Utilities/Logger.hpp"
 #include "ActsExamples/EventData/Cluster.hpp"
 #include "ActsExamples/EventData/Index.hpp"
@@ -23,17 +20,12 @@
 #include "ActsExamples/Framework/DataHandle.hpp"
 #include "ActsExamples/Framework/ProcessCode.hpp"
 #include "ActsExamples/Framework/WriterT.hpp"
-#include "ActsFatras/Digitization/Channelizer.hpp"
 
-#include <array>
 #include <memory>
-#include <mutex>
 #include <string>
 #include <unordered_map>
 #include <vector>
 
-#include <TTree.h>
-
 class TFile;
 class TTree;
 namespace Acts {
@@ -66,8 +58,12 @@ class RootMeasurementWriter final : public WriterT<MeasurementContainer> {
     std::string inputSimHits;
     /// Input collection to map measured hits to simulated hits.
     std::string inputMeasurementSimHitsMap;
-    std::string filePath = "";          ///< path of the output file
-    std::string fileMode = "RECREATE";  ///< file access mode
+
+    /// path of the output file
+    std::string filePath = "";
+    /// file access mode
+    std::string fileMode = "RECREATE";
+
     /// The indices for this digitization configurations
     Acts::GeometryHierarchyMap<std::vector<Acts::BoundIndices>> boundIndices;
     /// Map of the geometry identifier to the surface
@@ -75,186 +71,6 @@ class RootMeasurementWriter final : public WriterT<MeasurementContainer> {
         surfaceByIdentifier;
   };
 
-  struct DigitizationTree {
-    const std::array<std::string, Acts::eBoundSize> bNames = {
-        "loc0", "loc1", "phi", "theta", "qop", "time"};
-
-    TTree* tree = nullptr;
-    // Identification parameters
-    int eventNr = 0;
-    int volumeID = 0;
-    int layerID = 0;
-    int surfaceID = 0;
-
-    /// Type 0 - free, 1 - bound
-    int measType = 1;
-
-    /// Truth parameters
-    float trueBound[Acts::eBoundSize] = {};
-    float trueGx = 0.;
-    float trueGy = 0.;
-    float trueGz = 0.;
-    float incidentPhi = 0.;
-    float incidentTheta = 0.;
-
-    /// Reconstruction information
-    float recBound[Acts::eBoundSize] = {};
-    float varBound[Acts::eBoundSize] = {};
-
-    /// Cluster information comprised of
-    /// nch :  number of channels
-    /// cSize : cluster size in loc0 and loc1
-    /// chId : channel identification
-    /// chValue: value/activation of the channel
-    int nch = 0;
-    int cSize[2] = {};
-    std::array<std::vector<int>*, 2> chId = {nullptr, nullptr};
-    std::vector<float>* chValue = nullptr;
-
-    /// Setup helper to create the tree and
-    /// register the branches
-    ///
-    /// @param treeName the name of the tree to be registered
-    void setupTree(const std::string& treeName) {
-      tree = new TTree(treeName.c_str(), treeName.c_str());
-      // Declare the branches
-      tree->Branch("event_nr", &eventNr);
-      tree->Branch("volume_id", &volumeID);
-      tree->Branch("layer_id", &layerID);
-      tree->Branch("surface_id", &surfaceID);
-      tree->Branch("measurement_type", &measType);
-      for (unsigned int ib = 0; ib < Acts::eBoundSize; ++ib) {
-        if (ib != Acts::eBoundQOverP) {
-          tree->Branch(("true_" + bNames[ib]).c_str(), &trueBound[ib]);
-        }
-      }
-      tree->Branch("true_x", &trueGx);
-      tree->Branch("true_y", &trueGy);
-      tree->Branch("true_z", &trueGz);
-      tree->Branch("true_incident_phi", &incidentPhi);
-      tree->Branch("true_incident_theta", &incidentTheta);
-    }
-
-    /// Constructor from GeometryIdentifier
-    DigitizationTree(Acts::GeometryIdentifier geoID) {
-      auto vID = geoID.volume();
-      auto lID = geoID.layer();
-      auto mID = geoID.sensitive();
-      std::string treeName = "vol" + std::to_string(vID);
-      if (lID > 0) {
-        treeName += "_lay" + std::to_string(lID);
-      }
-      if (mID > 0) {
-        treeName += "_mod" + std::to_string(mID);
-      }
-      setupTree(treeName);
-    }
-
-    /// Non-trivial destructor for memory cleanup
-    ~DigitizationTree() {
-      delete chId[0];
-      delete chId[1];
-      delete chValue;
-    }
-
-    /// Setup the dimension depended branches
-    ///
-    /// @param i the bound index in question
-    void setupBoundRecBranch(Acts::BoundIndices i) {
-      tree->Branch(("rec_" + bNames[i]).c_str(), &recBound[i]);
-      tree->Branch(("var_" + bNames[i]).c_str(), &varBound[i]);
-    }
-
-    /// Setup the cluster related branch
-    ///
-    /// @param bIndices the bound indices to be written
-    void setupClusterBranch(const std::vector<Acts::BoundIndices>& bIndices) {
-      chValue = new std::vector<float>;
-      tree->Branch("clus_size", &nch);
-      tree->Branch("channel_value", &chValue);
-      // Both are allocated, but only relevant ones are set
-      chId[0] = new std::vector<int>;
-      chId[1] = new std::vector<int>;
-      for (const auto& ib : bIndices) {
-        if (static_cast<unsigned int>(ib) < 2) {
-          tree->Branch(("channel_" + bNames[ib]).c_str(), &chId[ib]);
-          tree->Branch(("clus_size_" + bNames[ib]).c_str(), &cSize[ib]);
-        }
-      }
-    }
-
-    /// Convenience function to register idenfication
-    ///
-    /// @param eventNr The event number
-    /// @param geoID The geometry identifier of the measurement
-    void fillIdentification(int evnt, Acts::GeometryIdentifier geoId) {
-      eventNr = evnt;
-      volumeID = geoId.volume();
-      layerID = geoId.layer();
-      surfaceID = geoId.sensitive();
-    }
-
-    /// Convenience function to register the truth parameters
-    ///
-    /// @param lp The true local position
-    /// @param xt The true 4D global position
-    /// @param dir The true particle direction
-    void fillTruthParameters(const Acts::Vector2& lp, const Acts::Vector4& xt,
-                             const Acts::Vector3& dir,
-                             const std::pair<double, double> angles) {
-      trueBound[Acts::eBoundLoc0] = lp[Acts::eBoundLoc0];
-      trueBound[Acts::eBoundLoc1] = lp[Acts::eBoundLoc1];
-      trueBound[Acts::eBoundPhi] = Acts::VectorHelpers::phi(dir);
-      trueBound[Acts::eBoundTheta] = Acts::VectorHelpers::theta(dir);
-      trueBound[Acts::eBoundTime] = xt[Acts::eTime];
-
-      trueGx = xt[Acts::ePos0];
-      trueGy = xt[Acts::ePos1];
-      trueGz = xt[Acts::ePos2];
-
-      incidentPhi = angles.first;
-      incidentTheta = angles.second;
-    }
-
-    /// Convenience function to fill bound parameters
-    ///
-    /// @tparam measurement_t Type of the parameter set
-    ///
-    /// @param m The measurement set
-    template <typename measurement_t>
-    void fillBoundMeasurement(const measurement_t& m) {
-      Acts::BoundVector fullVect = m.expander() * m.parameters();
-      recBound[Acts::eBoundLoc0] = fullVect[Acts::eBoundLoc0];
-      recBound[Acts::eBoundLoc1] = fullVect[Acts::eBoundLoc1];
-      recBound[Acts::eBoundPhi] = fullVect[Acts::eBoundPhi];
-      recBound[Acts::eBoundTheta] = fullVect[Acts::eBoundTheta];
-      recBound[Acts::eBoundTime] = fullVect[Acts::eBoundTime];
-
-      Acts::BoundSquareMatrix fullVar =
-          m.expander() * m.covariance() * m.expander().transpose();
-      varBound[Acts::eBoundLoc0] = fullVar(Acts::eBoundLoc0, Acts::eBoundLoc0);
-      varBound[Acts::eBoundLoc1] = fullVar(Acts::eBoundLoc1, Acts::eBoundLoc1);
-      varBound[Acts::eBoundPhi] = fullVar(Acts::eBoundPhi, Acts::eBoundPhi);
-      varBound[Acts::eBoundTheta] =
-          fullVar(Acts::eBoundTheta, Acts::eBoundTheta);
-      varBound[Acts::eBoundTime] = fullVar(Acts::eBoundTime, Acts::eBoundTime);
-    }
-
-    /// Convenience function to fill the cluster information
-    ///
-    /// @param c The cluster
-    void fillCluster(const Cluster& c) {
-      nch = static_cast<int>(c.channels.size());
-      cSize[0] = static_cast<int>(c.sizeLoc0);
-      cSize[1] = static_cast<int>(c.sizeLoc1);
-      for (auto ch : c.channels) {
-        chId[0]->push_back(static_cast<int>(ch.bin[0]));
-        chId[1]->push_back(static_cast<int>(ch.bin[1]));
-        chValue->push_back(static_cast<float>(ch.activation));
-      }
-    }
-  };
-
   /// Constructor with
   /// @param cfg configuration struct
   /// @param output logging level
@@ -279,18 +95,20 @@ class RootMeasurementWriter final : public WriterT<MeasurementContainer> {
                      const MeasurementContainer& measurements) override;
 
  private:
+  struct DigitizationTree;
+
   Config m_cfg;
-  std::mutex m_writeMutex;  ///< protect multi-threaded writes
-  TFile* m_outputFile;      ///< the output file
-  Acts::GeometryHierarchyMap<std::unique_ptr<DigitizationTree>>
-      m_outputTrees;  ///< the output trees
-  std::unordered_map<Acts::GeometryIdentifier, const Acts::Surface*>
-      m_dSurfaces;  ///< All surfaces that could carry measurements
+  /// protect multi-threaded writes
+  std::mutex m_writeMutex;
+  /// the output file
+  TFile* m_outputFile;
+  /// the output tree
+  std::unique_ptr<DigitizationTree> m_outputTree;
 
+  ReadDataHandle<ClusterContainer> m_inputClusters{this, "InputClusters"};
   ReadDataHandle<SimHitContainer> m_inputSimHits{this, "InputSimHits"};
   ReadDataHandle<IndexMultimap<Index>> m_inputMeasurementSimHitsMap{
       this, "InputMeasurementSimHitsMap"};
-  ReadDataHandle<ClusterContainer> m_inputClusters{this, "InputClusters"};
 };
 
 }  // namespace ActsExamples
diff --git a/Examples/Io/Root/src/RootMeasurementWriter.cpp b/Examples/Io/Root/src/RootMeasurementWriter.cpp
index 4ab08e7744e..9cc695596e4 100644
--- a/Examples/Io/Root/src/RootMeasurementWriter.cpp
+++ b/Examples/Io/Root/src/RootMeasurementWriter.cpp
@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2017 CERN for the benefit of the Acts project
+// Copyright (C) 2017-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -8,6 +8,8 @@
 
 #include "ActsExamples/Io/Root/RootMeasurementWriter.hpp"
 
+#include "Acts/Definitions/TrackParametrization.hpp"
+#include "Acts/Utilities/Enumerate.hpp"
 #include "ActsExamples/EventData/AverageSimHits.hpp"
 #include "ActsExamples/EventData/Index.hpp"
 #include "ActsExamples/EventData/IndexSourceLink.hpp"
@@ -16,19 +18,195 @@
 
 #include <cstddef>
 #include <ios>
+#include <limits>
+#include <memory>
 #include <stdexcept>
 #include <utility>
 #include <variant>
 
 #include <TFile.h>
+#include <TTree.h>
 
-namespace Acts {
-class Surface;
-}  // namespace Acts
+namespace ActsExamples {
 
-ActsExamples::RootMeasurementWriter::RootMeasurementWriter(
-    const ActsExamples::RootMeasurementWriter::Config& config,
-    Acts::Logging::Level level)
+struct RootMeasurementWriter::DigitizationTree {
+  const std::array<std::string, Acts::eBoundSize> bNames = {
+      "loc0", "loc1", "phi", "theta", "qop", "time"};
+
+  TTree* tree = nullptr;
+
+  // Identification parameters
+  int eventNr = 0;
+  int volumeID = 0;
+  int layerID = 0;
+  int surfaceID = 0;
+
+  // Reconstruction information
+  float recBound[Acts::eBoundSize] = {};
+  float varBound[Acts::eBoundSize] = {};
+
+  // Truth parameters
+  float trueBound[Acts::eBoundSize] = {};
+  float trueGx = 0.;
+  float trueGy = 0.;
+  float trueGz = 0.;
+  float incidentPhi = 0.;
+  float incidentTheta = 0.;
+
+  // Residuals and pulls
+  float residual[Acts::eBoundSize] = {};
+  float pull[Acts::eBoundSize] = {};
+
+  // Cluster information comprised of
+  // nch :  number of channels
+  // cSize : cluster size in loc0 and loc1
+  // chId : channel identification
+  // chValue: value/activation of the channel
+  int nch = 0;
+  int cSize[2] = {};
+  std::array<std::vector<int>, 2> chId;
+  std::vector<float> chValue;
+
+  /// Constructor from tree name
+  DigitizationTree(const std::string& treeName,
+                   const std::vector<Acts::BoundIndices>& recoIndices,
+                   const std::vector<Acts::BoundIndices>& clusterIndices) {
+    tree = new TTree(treeName.c_str(), treeName.c_str());
+
+    tree->Branch("event_nr", &eventNr);
+    tree->Branch("volume_id", &volumeID);
+    tree->Branch("layer_id", &layerID);
+    tree->Branch("surface_id", &surfaceID);
+
+    for (auto ib : recoIndices) {
+      tree->Branch(("rec_" + bNames[ib]).c_str(), &recBound[ib]);
+    }
+    for (auto ib : recoIndices) {
+      tree->Branch(("var_" + bNames[ib]).c_str(), &varBound[ib]);
+    }
+
+    tree->Branch("clus_size", &nch);
+    tree->Branch("channel_value", &chValue);
+    // Both are allocated, but only relevant ones are set
+    for (auto ib : clusterIndices) {
+      if (static_cast<unsigned int>(ib) < 2) {
+        tree->Branch(("channel_" + bNames[ib]).c_str(), &chId[ib]);
+        tree->Branch(("clus_size_" + bNames[ib]).c_str(), &cSize[ib]);
+      }
+    }
+
+    for (unsigned int ib = 0; ib < Acts::eBoundSize; ++ib) {
+      tree->Branch(("true_" + bNames[ib]).c_str(), &trueBound[ib]);
+    }
+    tree->Branch("true_x", &trueGx);
+    tree->Branch("true_y", &trueGy);
+    tree->Branch("true_z", &trueGz);
+    tree->Branch("true_incident_phi", &incidentPhi);
+    tree->Branch("true_incident_theta", &incidentTheta);
+
+    for (auto ib : recoIndices) {
+      tree->Branch(("residual_" + bNames[ib]).c_str(), &residual[ib]);
+    }
+    for (auto ib : recoIndices) {
+      tree->Branch(("pull_" + bNames[ib]).c_str(), &pull[ib]);
+    }
+
+    clear();
+  }
+
+  /// Convenience function to register idenfication
+  ///
+  /// @param eventNr The event number
+  /// @param geoID The geometry identifier of the measurement
+  void fillIdentification(int evnt, Acts::GeometryIdentifier geoId) {
+    eventNr = evnt;
+    volumeID = geoId.volume();
+    layerID = geoId.layer();
+    surfaceID = geoId.sensitive();
+  }
+
+  /// Convenience function to register the truth parameters
+  ///
+  /// @param lp The true local position
+  /// @param xt The true 4D global position
+  /// @param dir The true particle direction
+  /// @param angles The incident angles
+  /// @param qop The true q/p
+  void fillTruthParameters(const Acts::Vector2& lp, const Acts::Vector4& xt,
+                           const Acts::Vector3& dir,
+                           const std::pair<double, double> angles) {
+    trueBound[Acts::eBoundLoc0] = lp[Acts::eBoundLoc0];
+    trueBound[Acts::eBoundLoc1] = lp[Acts::eBoundLoc1];
+    trueBound[Acts::eBoundPhi] = Acts::VectorHelpers::phi(dir);
+    trueBound[Acts::eBoundTheta] = Acts::VectorHelpers::theta(dir);
+    trueBound[Acts::eBoundTime] = xt[Acts::eTime];
+
+    trueGx = xt[Acts::ePos0];
+    trueGy = xt[Acts::ePos1];
+    trueGz = xt[Acts::ePos2];
+
+    incidentPhi = angles.first;
+    incidentTheta = angles.second;
+  }
+
+  /// Convenience function to fill bound parameters
+  ///
+  /// @tparam measurement_t Type of the parameter set
+  ///
+  /// @param m The measurement set
+  template <typename measurement_t>
+  void fillBoundMeasurement(const measurement_t& m) {
+    for (auto [i, ib] : Acts::enumerate(m.indices())) {
+      recBound[ib] = m.parameters()[i];
+      varBound[ib] = m.covariance()(i, i);
+
+      residual[ib] = recBound[ib] - trueBound[ib];
+      pull[ib] = residual[ib] / std::sqrt(varBound[ib]);
+    }
+  }
+
+  /// Convenience function to fill the cluster information
+  ///
+  /// @param c The cluster
+  void fillCluster(const Cluster& c) {
+    nch = static_cast<int>(c.channels.size());
+    cSize[0] = static_cast<int>(c.sizeLoc0);
+    cSize[1] = static_cast<int>(c.sizeLoc1);
+    for (auto ch : c.channels) {
+      chId[0].push_back(static_cast<int>(ch.bin[0]));
+      chId[1].push_back(static_cast<int>(ch.bin[1]));
+      chValue.push_back(static_cast<float>(ch.activation));
+    }
+  }
+
+  /// Fill the tree
+  void fill() { tree->Fill(); }
+
+  /// Clear the tree
+  void clear() {
+    for (unsigned int ib = 0; ib < Acts::eBoundSize; ++ib) {
+      trueBound[ib] = std::numeric_limits<float>::quiet_NaN();
+      recBound[ib] = std::numeric_limits<float>::quiet_NaN();
+      varBound[ib] = std::numeric_limits<float>::quiet_NaN();
+      residual[ib] = std::numeric_limits<float>::quiet_NaN();
+      pull[ib] = std::numeric_limits<float>::quiet_NaN();
+    }
+    trueGx = std::numeric_limits<float>::quiet_NaN();
+    trueGy = std::numeric_limits<float>::quiet_NaN();
+    trueGz = std::numeric_limits<float>::quiet_NaN();
+    incidentPhi = std::numeric_limits<float>::quiet_NaN();
+    incidentTheta = std::numeric_limits<float>::quiet_NaN();
+    nch = 0;
+    cSize[0] = 0;
+    cSize[1] = 0;
+    chId[0].clear();
+    chId[1].clear();
+    chValue.clear();
+  }
+};
+
+RootMeasurementWriter::RootMeasurementWriter(
+    const RootMeasurementWriter::Config& config, Acts::Logging::Level level)
     : WriterT(config.inputMeasurements, "RootMeasurementWriter", level),
       m_cfg(config) {
   // Input container for measurements is already checked by base constructor
@@ -40,9 +218,9 @@ ActsExamples::RootMeasurementWriter::RootMeasurementWriter(
         "Missing hit-to-simulated-hits map input collection");
   }
 
+  m_inputClusters.maybeInitialize(m_cfg.inputClusters);
   m_inputSimHits.initialize(m_cfg.inputSimHits);
   m_inputMeasurementSimHitsMap.initialize(m_cfg.inputMeasurementSimHitsMap);
-  m_inputClusters.maybeInitialize(m_cfg.inputClusters);
 
   if (m_cfg.surfaceByIdentifier.empty()) {
     throw std::invalid_argument("Missing Surface-GeoID association map");
@@ -55,59 +233,34 @@ ActsExamples::RootMeasurementWriter::RootMeasurementWriter(
 
   m_outputFile->cd();
 
-  // Analyze the smearers
-  std::vector<
-      std::pair<Acts::GeometryIdentifier, std::unique_ptr<DigitizationTree>>>
-      dTrees;
-  if (!m_cfg.boundIndices.empty()) {
-    ACTS_DEBUG("Bound indices are declared, preparing trees.");
-    for (std::size_t ikv = 0; ikv < m_cfg.boundIndices.size(); ++ikv) {
-      auto geoID = m_cfg.boundIndices.idAt(ikv);
-      auto bIndices = m_cfg.boundIndices.valueAt(ikv);
-      auto dTree = std::make_unique<DigitizationTree>(geoID);
-      for (const auto& bIndex : bIndices) {
-        ACTS_VERBOSE("- setup branch for index: " << bIndex);
-        dTree->setupBoundRecBranch(bIndex);
-      }
-      if (!m_cfg.inputClusters.empty()) {
-        dTree->setupClusterBranch(bIndices);
-      }
-      dTrees.push_back({geoID, std::move(dTree)});
-    }
-  } else {
-    ACTS_DEBUG("Bound indices are not declared, no reco setup.");
-  }
-
-  m_outputTrees = Acts::GeometryHierarchyMap<std::unique_ptr<DigitizationTree>>(
-      std::move(dTrees));
+  std::vector bIndices = {Acts::eBoundLoc0, Acts::eBoundLoc1, Acts::eBoundTime};
+  m_outputTree =
+      std::make_unique<DigitizationTree>("measurements", bIndices, bIndices);
 }
 
-ActsExamples::RootMeasurementWriter::~RootMeasurementWriter() {
+RootMeasurementWriter::~RootMeasurementWriter() {
   if (m_outputFile != nullptr) {
     m_outputFile->Close();
   }
 }
 
-ActsExamples::ProcessCode ActsExamples::RootMeasurementWriter::finalize() {
+ProcessCode RootMeasurementWriter::finalize() {
   /// Close the file if it's yours
   m_outputFile->cd();
-  for (auto dTree = m_outputTrees.begin(); dTree != m_outputTrees.end();
-       ++dTree) {
-    (*dTree)->tree->Write();
-  }
+  m_outputTree->tree->Write();
   m_outputFile->Close();
 
   return ProcessCode::SUCCESS;
 }
 
-ActsExamples::ProcessCode ActsExamples::RootMeasurementWriter::writeT(
+ProcessCode RootMeasurementWriter::writeT(
     const AlgorithmContext& ctx, const MeasurementContainer& measurements) {
   const auto& simHits = m_inputSimHits(ctx);
   const auto& hitSimHitsMap = m_inputMeasurementSimHitsMap(ctx);
 
-  ClusterContainer clusters;
+  const ClusterContainer* clusters = nullptr;
   if (!m_cfg.inputClusters.empty()) {
-    clusters = m_inputClusters(ctx);
+    clusters = &m_inputClusters(ctx);
   }
 
   // Exclusive access to the tree while writing
@@ -126,15 +279,9 @@ ActsExamples::ProcessCode ActsExamples::RootMeasurementWriter::writeT(
             return;
           }
           const Acts::Surface& surface = *(surfaceItr->second);
-          // find the corresponding output tree
-          auto dTreeItr = m_outputTrees.find(geoId);
-          if (dTreeItr == m_outputTrees.end()) {
-            return;
-          }
-          auto& dTree = *dTreeItr;
 
           // Fill the identification
-          dTree->fillIdentification(ctx.eventNumber, geoId);
+          m_outputTree->fillIdentification(ctx.eventNumber, geoId);
 
           // Find the contributing simulated hits
           auto indices = makeRange(hitSimHitsMap.equal_range(hitIdx));
@@ -148,25 +295,21 @@ ActsExamples::ProcessCode ActsExamples::RootMeasurementWriter::writeT(
                   .inverse();
           std::pair<double, double> angles =
               Acts::VectorHelpers::incidentAngles(dir, rot);
-          dTree->fillTruthParameters(local, pos4, dir, angles);
-          dTree->fillBoundMeasurement(m);
-          if (!clusters.empty()) {
-            const auto& c = clusters[hitIdx];
-            dTree->fillCluster(c);
-          }
-          dTree->tree->Fill();
-          if (dTree->chValue != nullptr) {
-            dTree->chValue->clear();
-          }
-          if (dTree->chId[0] != nullptr) {
-            dTree->chId[0]->clear();
-          }
-          if (dTree->chId[1] != nullptr) {
-            dTree->chId[1]->clear();
+
+          m_outputTree->fillTruthParameters(local, pos4, dir, angles);
+          m_outputTree->fillBoundMeasurement(m);
+          if (clusters != nullptr) {
+            const auto& c = (*clusters)[hitIdx];
+            m_outputTree->fillCluster(c);
           }
+
+          m_outputTree->fill();
+          m_outputTree->clear();
         },
         meas);
   }
 
-  return ActsExamples::ProcessCode::SUCCESS;
+  return ProcessCode::SUCCESS;
 }
+
+}  // namespace ActsExamples
diff --git a/Examples/Python/python/acts/examples/odd.py b/Examples/Python/python/acts/examples/odd.py
index c66eaf979fc..ec8cbcb9438 100644
--- a/Examples/Python/python/acts/examples/odd.py
+++ b/Examples/Python/python/acts/examples/odd.py
@@ -20,6 +20,13 @@ def getOpenDataDetector(
     odd_dir: Optional[Path] = None,
     logLevel=acts.logging.INFO,
 ):
+    """This function sets up the open data detector. Requires DD4hep.
+    Parameters
+    ----------
+    mdecorator: Material Decorator, take RootMaterialDecorator if non is given
+    odd_dir: if not given, try to get via ODD_PATH environment variable
+    logLevel: logging level
+    """
     import acts.examples.dd4hep
 
     customLogLevel = acts.examples.defaultLogging(logLevel=logLevel)
@@ -86,11 +93,9 @@ def geoid_hook(geoid, surface):
     )
     detector = acts.examples.dd4hep.DD4hepDetector()
 
-    config = acts.MaterialMapJsonConverter.Config()
     if mdecorator is None:
-        mdecorator = acts.JsonMaterialDecorator(
-            rConfig=config,
-            jFileName=str(odd_dir / "config/odd-material-mapping-config.json"),
+        mdecorator = acts.examples.RootMaterialDecorator(
+            fileName=str(odd_dir / "data/odd-material-maps.root"),
             level=customLogLevel(minLevel=acts.logging.WARNING),
         )
 
diff --git a/Examples/Python/python/acts/examples/reconstruction.py b/Examples/Python/python/acts/examples/reconstruction.py
index de2fefa4eaf..cff5eea0cff 100644
--- a/Examples/Python/python/acts/examples/reconstruction.py
+++ b/Examples/Python/python/acts/examples/reconstruction.py
@@ -478,7 +478,7 @@ def addTruthSmearedSeeding(
     selectedParticles: str,
     particleSmearingSigmas: ParticleSmearingSigmas,
     initialSigmas: Optional[List[float]],
-    initialVarInflation: List[float],
+    initialVarInflation: Optional[List[float]],
     particleHypothesis: Optional[acts.ParticleHypothesis],
     logLevel: acts.logging.Level = None,
 ):
@@ -1882,9 +1882,11 @@ def addVertexFitting(
     trackParameters: Optional[str] = None,
     outputProtoVertices: str = "protovertices",
     outputVertices: str = "fittedVertices",
-    seeder: Optional[acts.VertexSeedFinder] = acts.VertexSeedFinder.GaussianSeeder,
     vertexFinder: VertexFinder = VertexFinder.Truth,
+    maxIterations: Optional[int] = None,
+    initialVariances: Optional[List[float]] = None,
     useTime: Optional[bool] = False,
+    seeder: Optional[acts.VertexSeedFinder] = acts.VertexSeedFinder.GaussianSeeder,
     spatialBinExtent: Optional[float] = None,
     temporalBinExtent: Optional[float] = None,
     trackSelectorConfig: Optional[TrackSelectorConfig] = None,
@@ -1901,11 +1903,11 @@ def addVertexFitting(
     field : magnetic field
     outputDirRoot : Path|str, path, None
         the output folder for the Root output, None triggers no output
+    vertexFinder : VertexFinder, Truth
+        vertexFinder algorithm: one of Truth, AMVF, Iterative
     seeder : enum member
         determines vertex seeder for AMVF, can be acts.seeder.GaussianSeeder or
         acts.seeder.AdaptiveGridSeeder
-    vertexFinder : VertexFinder, Truth
-        vertexFinder algorithm: one of Truth, AMVF, Iterative
     useTime : bool, False
         determines whether time information is used in vertex seeder, finder,
         and fitter
@@ -1949,7 +1951,7 @@ def addVertexFitting(
     tracks = tracks if tracks is not None else ""
     inputParticles = "particles_input"
     selectedParticles = "particles_selected"
-    inputVertices = "vertices_input"
+    inputTruthVertices = "vertices_input"
 
     if vertexFinder == VertexFinder.Truth:
         findVertices = TruthVertexFinder(
@@ -1963,30 +1965,34 @@ def addVertexFitting(
         s.addAlgorithm(findVertices)
         fitVertices = VertexFitterAlgorithm(
             level=customLogLevel(),
-            bField=field,
             inputTrackParameters=trackParameters,
             inputProtoVertices=findVertices.config.outputProtoVertices,
             outputVertices=outputVertices,
+            bField=field,
         )
         s.addAlgorithm(fitVertices)
     elif vertexFinder == VertexFinder.Iterative:
         findVertices = IterativeVertexFinderAlgorithm(
             level=customLogLevel(),
-            bField=field,
             inputTrackParameters=trackParameters,
             outputProtoVertices=outputProtoVertices,
             outputVertices=outputVertices,
+            bField=field,
         )
         s.addAlgorithm(findVertices)
     elif vertexFinder == VertexFinder.AMVF:
         findVertices = AdaptiveMultiVertexFinderAlgorithm(
             level=customLogLevel(),
-            seedFinder=seeder,
-            bField=field,
             inputTrackParameters=trackParameters,
+            inputTruthParticles=selectedParticles,
+            inputTruthVertices=inputTruthVertices,
             outputProtoVertices=outputProtoVertices,
             outputVertices=outputVertices,
+            bField=field,
+            seedFinder=seeder,
             **acts.examples.defaultKWArgs(
+                maxIterations=maxIterations,
+                initialVariances=initialVariances,
                 useTime=useTime,
                 spatialBinExtent=spatialBinExtent,
                 temporalBinExtent=temporalBinExtent,
@@ -2005,7 +2011,7 @@ def addVertexFitting(
                 level=customLogLevel(),
                 inputVertices=outputVertices,
                 inputTracks=tracks,
-                inputTruthVertices=inputVertices,
+                inputTruthVertices=inputTruthVertices,
                 inputParticles=inputParticles,
                 inputSelectedParticles=selectedParticles,
                 inputTrackParticleMatching="track_particle_matching",
diff --git a/Examples/Python/python/acts/examples/simulation.py b/Examples/Python/python/acts/examples/simulation.py
index 76cf0685aa8..5b1ddccee26 100644
--- a/Examples/Python/python/acts/examples/simulation.py
+++ b/Examples/Python/python/acts/examples/simulation.py
@@ -452,21 +452,22 @@ def addFatras(
 
     # Selector
     if preSelectParticles is not None:
-        particles_selected = "particles_selected"
+        particlesSelected = "fatras_particles_preselected"
         addParticleSelection(
             s,
             preSelectParticles,
             inputParticles=inputParticles,
-            outputParticles=particles_selected,
+            outputParticles=particlesSelected,
         )
+        s.addWhiteboardAlias("particles_selected", particlesSelected)
     else:
-        particles_selected = inputParticles
+        particlesSelected = inputParticles
 
     # Simulation
     alg = acts.examples.FatrasSimulation(
         **acts.examples.defaultKWArgs(
             level=customLogLevel(),
-            inputParticles=particles_selected,
+            inputParticles=particlesSelected,
             outputParticlesInitial=outputParticlesInitial,
             outputParticlesFinal=outputParticlesFinal,
             outputSimHits=outputSimHits,
@@ -502,6 +503,7 @@ def addFatras(
             inputParticles=alg.config.outputParticlesFinal,
             outputParticles=particlesFinal,
         )
+        s.addWhiteboardAlias("particles_selected", particlesFinal)
     else:
         particlesInitial = alg.config.outputParticlesInitial
         particlesFinal = alg.config.outputParticlesFinal
@@ -673,15 +675,16 @@ def addGeant4(
 
     # Selector
     if preSelectParticles is not None:
-        particles_selected = "particles_selected"
+        particlesSelected = "geant4_particles_preselected"
         addParticleSelection(
             s,
             preSelectParticles,
             inputParticles=inputParticles,
-            outputParticles=particles_selected,
+            outputParticles=particlesSelected,
         )
+        s.addWhiteboardAlias("particles_selected", particlesSelected)
     else:
-        particles_selected = inputParticles
+        particlesSelected = inputParticles
 
     if g4DetectorConstructionFactory is None:
         if detector is None:
@@ -705,7 +708,7 @@ def addGeant4(
         geant4Handle=__geant4Handle,
         detectorConstructionFactory=g4DetectorConstructionFactory,
         randomNumbers=rnd,
-        inputParticles=particles_selected,
+        inputParticles=particlesSelected,
         outputParticlesInitial=outputParticlesInitial,
         outputParticlesFinal=outputParticlesFinal,
         outputSimHits=outputSimHits,
@@ -728,7 +731,7 @@ def addGeant4(
 
     # Selector
     if postSelectParticles is not None:
-        particlesInitial = "geant4_particles_initial_selected"
+        particlesInitial = "geant4_particles_initial_postselected"
         addParticleSelection(
             s,
             postSelectParticles,
@@ -736,13 +739,14 @@ def addGeant4(
             outputParticles=particlesInitial,
         )
 
-        particlesFinal = "geant4_particles_final_selected"
+        particlesFinal = "geant4_particles_final_postselected"
         addParticleSelection(
             s,
             postSelectParticles,
             inputParticles=alg.config.outputParticlesFinal,
             outputParticles=particlesFinal,
         )
+        s.addWhiteboardAlias("particles_selected", particlesFinal)
     else:
         particlesInitial = alg.config.outputParticlesInitial
         particlesFinal = alg.config.outputParticlesFinal
diff --git a/Examples/Python/src/Vertexing.cpp b/Examples/Python/src/Vertexing.cpp
index 79c60021ec4..44f2ebb8433 100644
--- a/Examples/Python/src/Vertexing.cpp
+++ b/Examples/Python/src/Vertexing.cpp
@@ -32,19 +32,23 @@ void addVertexing(Context& ctx) {
   auto& m = ctx.get("main");
 
   py::enum_<Seeder>(m, "VertexSeedFinder")
+      .value("TruthSeeder", Seeder::TruthSeeder)
       .value("GaussianSeeder", Seeder::GaussianSeeder)
       .value("AdaptiveGridSeeder", Seeder::AdaptiveGridSeeder);
 
   ACTS_PYTHON_DECLARE_ALGORITHM(
       ActsExamples::AdaptiveMultiVertexFinderAlgorithm, mex,
       "AdaptiveMultiVertexFinderAlgorithm", inputTrackParameters,
-      outputProtoVertices, outputVertices, seedFinder, bField, useTime,
-      spatialBinExtent, temporalBinExtent);
+      inputTruthParticles, inputTruthVertices, outputProtoVertices,
+      outputVertices, seedFinder, bField, minWeight, doSmoothing, maxIterations,
+      useTime, tracksMaxZinterval, initialVariances, doFullSplitting,
+      tracksMaxSignificance, maxMergeVertexSignificance, spatialBinExtent,
+      temporalBinExtent);
 
   ACTS_PYTHON_DECLARE_ALGORITHM(ActsExamples::IterativeVertexFinderAlgorithm,
                                 mex, "IterativeVertexFinderAlgorithm",
                                 inputTrackParameters, outputProtoVertices,
-                                outputVertices, bField);
+                                outputVertices, bField, maxIterations);
 
   ACTS_PYTHON_DECLARE_ALGORITHM(ActsExamples::VertexFitterAlgorithm, mex,
                                 "VertexFitterAlgorithm", inputTrackParameters,
diff --git a/Examples/Python/tests/root_file_hashes.txt b/Examples/Python/tests/root_file_hashes.txt
index bced4603466..082d537dc98 100644
--- a/Examples/Python/tests/root_file_hashes.txt
+++ b/Examples/Python/tests/root_file_hashes.txt
@@ -38,12 +38,12 @@ test_material_mapping__material-map_tracks.root: 938b1a855369e9304401cb10d2751df
 test_material_mapping__propagation-material.root: d081c2c54403941bc7876fa552425eef29b496883dcdaab01113462895d0905e
 test_volume_material_mapping__material-map-volume_tracks.root: 98e212d32ca054fa3d01af4167c1f49755a139d43b82c57908197f5985e0a4ff
 test_volume_material_mapping__propagation-volume-material.root: cc3a2a6df914526be7e4a0c33db615a9589441b6bfb763271ab7ce7ea7f8789e
-test_digitization_example[smeared]__measurements.root: 03184f20b7b8e69f40424f76e2fd0620fb9f319d23ad69aa4bdd25c1b5cd482e
-test_digitization_example[geometric]__measurements.root: bc5d158478d7fdc4e09728cc9b4609afaa62be9dfcc1162685d03160f4d89495
+test_digitization_example[smeared]__measurements.root: 2b583b886b76f94786c6c2832aa84d176059cbbc35882fb34b4fd1f22dd4c006
+test_digitization_example[geometric]__measurements.root: 85efa861d14207fd7d1798dab093edcd0a2685bc189b4fc9a96b07d1001013a0
 test_digitization_example_input[smeared]__particles.root: 5fe7dda2933ee6b9615b064d192322fe07831133cd998e5ed99a3b992b713a10
-test_digitization_example_input[smeared]__measurements.root: 02a12e9949da9ba40b0120462d822f89e55562956fe3811981bdfd8648b7d946
+test_digitization_example_input[smeared]__measurements.root: 243c2f69b7b0db9dbeaa7494d4ea0f3dd1691dc90f16e10df6c0491ff4dc7d62
 test_digitization_example_input[geometric]__particles.root: 5fe7dda2933ee6b9615b064d192322fe07831133cd998e5ed99a3b992b713a10
-test_digitization_example_input[geometric]__measurements.root: 919abc40e61811bb1c591b6b5125594b922fd5a4bfd65d6b9f653e6a9292f662
+test_digitization_example_input[geometric]__measurements.root: 393e5d0278dceedb9be0efccafc5fceb024e2ef8ed919f351910c60606399f88
 test_ckf_tracks_example[generic-full_seeding]__trackstates_ckf.root: fa3460d2bed49028a9b6c0b3553486b8fddc979d6e508a1813780b4ea1a04f02
 test_ckf_tracks_example[generic-full_seeding]__tracksummary_ckf.root: b8ddfc25b013ec181e2a9f77925b86e57528140ccb6250696a9869272e1afa10
 test_ckf_tracks_example[generic-full_seeding]__performance_seeding_trees.root: 0e0676ffafdb27112fbda50d1cf627859fa745760f98073261dcf6db3f2f991e
@@ -73,7 +73,7 @@ test_root_prop_step_writer[kwargsConstructor]__prop_steps.root: 02460dea7b0b56b5
 test_root_particle_writer[configPosConstructor]__particles.root: 759f7c0225a30afc28c3dda6010295207254ebb0194b88f692d538b8c65d4326
 test_root_particle_writer[configKwConstructor]__particles.root: 759f7c0225a30afc28c3dda6010295207254ebb0194b88f692d538b8c65d4326
 test_root_particle_writer[kwargsConstructor]__particles.root: 759f7c0225a30afc28c3dda6010295207254ebb0194b88f692d538b8c65d4326
-test_root_meas_writer__meas.root: c92ef6936a7a20d3eb09722e9fab6aa91286439ff57b5c89c9f75c1fd156684e
+test_root_meas_writer__meas.root: 0f7e51c3c6e0133c6c1fc0da6d936e16f59ca2b579d15108e6c7b82bdf0a1e73
 test_root_simhits_writer[configPosConstructor]__meas.root: 66fee65a672bbf28098aac7af59ef181c2ac8331f87a379132873da86bdb1ba5
 test_root_simhits_writer[configKwConstructor]__meas.root: 66fee65a672bbf28098aac7af59ef181c2ac8331f87a379132873da86bdb1ba5
 test_root_simhits_writer[kwargsConstructor]__meas.root: 66fee65a672bbf28098aac7af59ef181c2ac8331f87a379132873da86bdb1ba5
diff --git a/Examples/Python/tests/test_examples.py b/Examples/Python/tests/test_examples.py
index d92867a770c..032ee122138 100644
--- a/Examples/Python/tests/test_examples.py
+++ b/Examples/Python/tests/test_examples.py
@@ -28,7 +28,7 @@
     GenericDetector,
     AlignedDetector,
 )
-from acts.examples.odd import getOpenDataDetector
+from acts.examples.odd import getOpenDataDetector, getOpenDataDetectorDirectory
 
 
 u = acts.UnitConstants
@@ -628,6 +628,23 @@ def test_truth_tracking_gsf(tmp_path, assert_root_hash, detector_config):
             assert_root_hash(fn, fp)
 
 
+def test_refitting(tmp_path, detector_config):
+    from truth_tracking_gsf_refitting import runRefittingGsf
+
+    field = acts.ConstantBField(acts.Vector3(0, 0, 2 * u.T))
+
+    seq = Sequencer(events=1, numThreads=1)
+
+    # Only check if it runs without errors right known
+    # Changes in fitter behaviour should be caught by other tests
+    runRefittingGsf(
+        trackingGeometry=detector_config.trackingGeometry,
+        field=field,
+        outputDir=tmp_path,
+        s=seq,
+    ).run()
+
+
 def test_particle_gun(tmp_path, assert_root_hash):
     from particle_gun import runParticleGun
 
@@ -661,7 +678,14 @@ def test_material_mapping(material_recording, tmp_path, assert_root_hash):
 
     s = Sequencer(numThreads=1)
 
-    detector, trackingGeometry, decorators = getOpenDataDetector()
+    odd_dir = getOpenDataDetectorDirectory()
+    config = acts.MaterialMapJsonConverter.Config()
+    mdecorator = acts.JsonMaterialDecorator(
+        level=acts.logging.INFO,
+        rConfig=config,
+        jFileName=str(odd_dir / "config/odd-material-mapping-config.json"),
+    )
+    detector, trackingGeometry, decorators = getOpenDataDetector(mdecorator)
 
     from material_mapping import runMaterialMapping
 
@@ -916,25 +940,6 @@ def test_digitization_example(trk_geo, tmp_path, assert_root_hash, digi_config_f
     assert len(list(csv_dir.iterdir())) == 3 * s.config.events
     assert all(f.stat().st_size > 50 for f in csv_dir.iterdir())
 
-    assert_entries(root_file, "vol9", 0)
-    assert_entries(root_file, "vol14", 0)
-
-    if "smearing" in digi_config_file.name:
-        filled_entries = [f"vol{tn}" for tn in (8, 12, 13, 16, 17, 18)]
-    else:
-        # fmt: off
-        filled_entries = [
-            'vol8', 'vol8_lay2', 'vol12_lay8_mod117', 'vol12_lay10', 'vol12_lay10_mod154',
-            'vol12_lay10_mod163', 'vol12_lay12', 'vol12_lay12_mod150', 'vol13',
-            'vol13_lay2', 'vol16_lay2_mod53', 'vol16_lay4', 'vol16_lay6', 'vol16_lay8',
-            'vol16_lay10', 'vol16_lay12', 'vol17', 'vol17_lay2', 'vol18_lay2',
-            'vol18_lay2_mod1', 'vol18_lay2_mod49', 'vol18_lay2_mod86', 'vol18_lay4',
-        ]
-        # fmt: on
-
-    for entry in filled_entries:
-        assert_has_entries(root_file, entry)
-
     assert_root_hash(root_file.name, root_file)
 
 
@@ -990,29 +995,6 @@ def test_digitization_example_input(
     assert len(list(csv_dir.iterdir())) == 3 * pgs.config.events
     assert all(f.stat().st_size > 50 for f in csv_dir.iterdir())
 
-    assert_entries(root_file, "vol7", 0)
-    assert_entries(root_file, "vol9", 0)
-
-    if "smearing" in digi_config_file.name:
-        filled_entries = [f"vol{tn}" for tn in (8, 12, 13, 16, 17, 18)]
-    else:
-        # fmt: off
-        filled_entries = [
-            "vol8", "vol8_lay2", "vol12_lay8_mod120", "vol12_lay10_mod120",
-            "vol12_lay10_mod144", "vol12_lay12", "vol12_lay12_mod111",
-            "vol12_lay12_mod137", "vol12_lay12_mod170", "vol13", "vol13_lay2",
-            "vol14_lay2_mod93", "vol14_lay2_mod102", "vol14_lay2_mod112",
-            "vol14_lay2_mod118", "vol14_lay4_mod112", "vol14_lay4_mod118",
-            "vol14_lay4_mod152", "vol14_lay4_mod161", "vol16_lay4", "vol16_lay6",
-            "vol16_lay8", "vol16_lay10", "vol16_lay12", "vol17", "vol17_lay2",
-            "vol18_lay2", "vol18_lay2_mod71", "vol18_lay4", "vol18_lay6",
-            "vol18_lay8", "vol18_lay10"
-        ]
-        # fmt: on
-
-    for entry in filled_entries:
-        assert_has_entries(root_file, entry)
-
     assert_root_hash(root_file.name, root_file)
 
 
diff --git a/Examples/Scripts/MaterialMapping/material_comparison.py b/Examples/Scripts/MaterialMapping/material_comparison.py
index 74ee6281b75..0c980af3287 100644
--- a/Examples/Scripts/MaterialMapping/material_comparison.py
+++ b/Examples/Scripts/MaterialMapping/material_comparison.py
@@ -1,239 +1,125 @@
-import ROOT
-import argparse
+import matplotlib.pyplot as plt
+import numpy as np
+import uproot
 import math
+from collections import namedtuple
 
-from ROOT import (
-    TCanvas,
-    TFile,
-    TTree,
-    gDirectory,
-    gStyle,
-    TH1D,
-    TH2F,
-    TProfile,
-    TRatioPlot,
+FileRecord = namedtuple("FileRecord", ["name", "tree", "label", "color", "marker_size"])
+PlotRecord = namedtuple(
+    "PlotRecord", ["x_axis", "y_axis", "x_range", "x_bins", "x_label", "saveAs"]
 )
 
-
-def TH1D_from_TProf(tprof):
-    h = TH1D(
-        tprof.GetName() + "_th1",
-        tprof.GetTitle(),
-        tprof.GetNbinsX(),
-        tprof.GetXaxis().GetXmin(),
-        tprof.GetXaxis().GetXmax(),
-    )
-    for i in range(tprof.GetNbinsX()):
-        if tprof.GetBinContent(i + 1) == 0.0:
-            h.SetBinContent(i + 1, 0.0)
-            h.SetBinError(i + 1, 1000.0)
-            continue
-        h.SetBinContent(i + 1, tprof.GetBinContent(i + 1))
-        h.SetBinError(i + 1, tprof.GetBinError(i + 1))
-    return h
-
-
-if "__main__" == __name__:
-    p = argparse.ArgumentParser()
-
-    p.add_argument(
-        "-e", "--entries", type=int, default=100000, help="Number of events to process"
-    )
-    p.add_argument(
-        "-i",
-        "--input",
-        type=str,
-        nargs="+",
-        default="",
-        help="Input files with material tracks",
-    )
-    p.add_argument(
-        "-o",
-        "--output",
-        type=str,
-        default="",
-        help="Output file with produced material overview plots",
-    )
-    p.add_argument(
-        "-l",
-        "--labels",
-        type=str,
-        nargs="+",
-        default="",
-        help="The labels for the input files",
-    )
-    p.add_argument(
-        "-v",
-        "--variables",
-        type=str,
-        nargs="+",
-        default=["t_X0", "t_L0"],
-        help="The variables to be plotted",
-    )
-    p.add_argument(
-        "-m",
-        "--max",
-        type=float,
-        nargs="+",
-        default=[4.5, 1.5],
-        help="The variables to be plotted",
-    )
-    p.add_argument(
-        "-a",
-        "--axes",
-        type=str,
-        nargs="+",
-        default=["v_eta", "v_phi"],
-        help="The axes versus which the variables to be plotted",
-    )
-    p.add_argument(
-        "--eta",
-        type=float,
-        nargs=2,
-        default=[-4.0, 4.0],
-        help="Eta range for the plotting",
-    )
-    p.add_argument("--eta-bins", type=int, default=60, help="Eta bins for the plotting")
-    p.add_argument(
-        "--phi",
-        type=float,
-        nargs=2,
-        default=[-math.pi, math.pi],
-        help="Phi range for the plotting",
-    )
-    p.add_argument("--phi-bins", type=int, default=72, help="Phi bins for the plotting")
-    p.add_argument(
-        "-d",
-        "--detray-input",
-        type=str,
-        default="",
-        help="Optional detray input csv file",
-    )
-
-    args = p.parse_args()
-
-    ttrees = []
-
-    if len(args.input) != len(args.labels):
-        print("** ERROR ** The number of input files and labels must match")
-        exit(1)
-
-    # The histogram map
-    h_dict = {}
-    tfiles = []
-
-    ofile = ROOT.TFile.Open(args.output, "RECREATE") if args.output else None
-
-    input_files = args.input
-
-    # Detray csv file
-    if args.detray_input:
-        input_files.append(args.detray_input)
-
-    # Loop over the files and create the comparison histograms
-    for fi, ifile in enumerate(input_files):
-        # Special treament for detray input
-        if ifile == args.detray_input:
-            print("Reading detray input from file: ", args.detray_input)
-            ttree = TTree("t", "material read from detray")
-            ttree.ReadFile(args.detray_input, "v_eta/D:v_phi:t_X0:t_L0:t_X0:t_L0")
-            lbl = "detray"
-        else:
-            # Get the file
-            tfile = ROOT.TFile.Open(ifile)
-            tfiles.append(tfile)
-            ttree = tfile.Get("material-tracks")
-            ttrees.append(ttree)
-            # The label
-            lbl = args.labels[fi]
-
-        # Loop over the variables and axes
-        for iv, var in enumerate(args.variables):
-            for ax in args.axes:
-                if ax == "v_eta":
-                    bins = args.eta_bins
-                    low = args.eta[0]
-                    high = args.eta[1]
-                    cut = f"v_phi > {args.phi[0]} && v_phi < {args.phi[1]}"
-                elif ax == "v_phi":
-                    bins = args.phi_bins
-                    low = args.phi[0]
-                    high = args.phi[1]
-                    cut = f"v_eta > {args.eta[0]} && v_eta < {args.eta[1]}"
-                # Some naming magic
-                hcmmd = f"{var}:{ax}>>"
-                hname = f"{var}_vs_{ax}"
-                hfname = f"{lbl}_{var}_vs_{ax}"
-                hrange = f"({bins},{low},{high},100,0,{args.max[iv]})"
-                htitle = f"{var} vs {ax}"
-                ttree.Draw(hcmmd + hfname + hrange, cut, "")
-                h = ROOT.gDirectory.Get(hfname)
-                # Fill into comparison histogram
-                if h_dict.get(hname):
-                    h_dict[hname].append(h)
-                else:
-                    h_dict[hname] = [h]
-
-                # Write to file
-                if ofile:
-                    ofile.cd()
-                    h.Write()
-
-    colors = [
-        ROOT.kBlack,
-        ROOT.kRed + 2,
-        ROOT.kBlue - 4,
-        ROOT.kGreen + 1,
-        ROOT.kYellow - 2,
-    ]
-    markers = [
-        ROOT.kFullCircle,
-        ROOT.kFullCircle,
-        ROOT.kFullTriangleUp,
-        ROOT.kFullTriangleDown,
-        ROOT.kFullDiamond,
-    ]
-
-    # Now create the comparison histograms
-    c = ROOT.TCanvas("Comparison", "Comparison", 1200, 800)
-    c.Divide(2, 2)
-
-    # Remove the stat box
-    gStyle.SetOptStat(0)
-
-    # Memory garbage collection, thanks ROOT
-    hist_memory_pool = []
-
-    ic = 0
-    for hname in h_dict:
-        ic += 1
-        c.cd(ic)
-        h_list = h_dict[hname]
-        h_ref = None
-        h_prof_ref = None
-        for ih, h in enumerate(h_list):
-            h_prof = TH1D_from_TProf(h.ProfileX())
-            hist_memory_pool.append(h_prof)
-            h_prof.SetObjectStat(0)
-            h_prof.SetLineColor(colors[ih])
-            h_prof.SetMarkerColor(colors[ih])
-            h_prof.SetMarkerSize(0.5)
-            h_prof.SetMarkerStyle(markers[ih])
-            h_prof.GetYaxis().SetRangeUser(0.0, 1.3 * h_prof.GetMaximum())
-            if ih == 0:
-                h_ref = h
-                h_prof_ref = h_prof
-            else:
-                h_ratio = ROOT.TRatioPlot(h_prof_ref, h_prof)
-                h_ratio.SetGraphDrawOpt("pe")
-                h_ratio.SetSeparationMargin(0.005)
-                drawOption = "e,same" if ih > 1 else "e"
-                h_ratio.Draw(drawOption)
-                h_ratio.GetLowerRefGraph().SetLineColor(colors[ih])
-                h_ratio.GetLowerRefGraph().SetMarkerColor(colors[ih])
-                h_ratio.GetLowerRefGraph().SetMarkerSize(0.5)
-                h_ratio.GetLowerRefGraph().SetMarkerStyle(markers[ih])
-                c.Update()
-                hist_memory_pool.append(h_ratio)
-        h_prof.Draw(("same" if ih > 0 else ""))
-        c.SaveAs(f"{hname}.png")
+# The file records
+fileRecords = [
+    FileRecord("geant4_material_tracks.root", "material-tracks", "Geant4", "blue", 3),
+    FileRecord("acts_material_tracks.root", "material-tracks", "Acts", "orange", 4),
+]
+
+
+# The plot records
+plotRecords = [
+    PlotRecord("v_eta", "t_X0", (-4.0, 4.0), 80, "η", "tX0_vs_eta.svg"),
+    PlotRecord("v_phi", "t_X0", (-math.pi, math.pi), 72, "φ", "tX0_vs_phi.svg"),
+]
+
+# Different plot records
+for pr in plotRecords:
+
+    fig, axs = plt.subplots(2, 1, sharex=True)
+    fig.subplots_adjust(hspace=0.05)
+
+    # Prepare limit & ratios
+    y_lim = 0
+    y_ratio_values = []
+    y_ratio_errors = [0.0 for i in range(pr.x_bins)]
+
+    # Loop over the file records
+    for ifr, fr in enumerate(fileRecords):
+
+        # Load the three
+        tree = uproot.open(fr.name + ":" + fr.tree)
+
+        x_arr = tree[pr.x_axis].array(library="np")
+        y_arr = tree[pr.y_axis].array(library="np")
+        y_max = y_arr.max()
+        y_lim = y_max if y_max > y_lim else y_lim
+
+        # Generate the central bin values
+        x_step = (pr.x_range[1] - pr.x_range[0]) / pr.x_bins
+        x_vals = [pr.x_range[0] + (ix + 0.5) * x_step for ix in range(pr.x_bins)]
+
+        # Prepare the min /max
+        y_min_vals = [1000.0] * pr.x_bins
+        y_max_vals = [0.0] * pr.x_bins
+        y_vals_sorted = [np.array([])] * pr.x_bins
+
+        for iv in range(len(x_arr)):
+            x_b = int((x_arr[iv] - pr.x_range[0]) / x_step)
+            y_v = y_arr[iv]
+            # Take min / max
+            y_min_vals[x_b] = y_v if y_v < y_min_vals[x_b] else y_min_vals[x_b]
+            y_max_vals[x_b] = y_v if y_v > y_max_vals[x_b] else y_max_vals[x_b]
+            # Regulate the x value
+            y_vals_sorted[x_b] = np.append(y_vals_sorted[x_b], y_v)
+
+        axs[0].fill_between(
+            x=x_vals,
+            y1=y_min_vals,
+            y2=y_max_vals,
+            alpha=0.1,
+            label=fr.label + " spread",
+            color=fr.color,
+        )
+        axs[0].grid(axis="x")
+        y_vals_mean = [y_bin_vals.mean() for y_bin_vals in y_vals_sorted]
+
+        y_ratio_values += [y_vals_mean]
+        y_vals_mse = [
+            y_bin_vals.std() ** 2 / len(y_bin_vals) for y_bin_vals in y_vals_sorted
+        ]
+
+        axs[0].errorbar(
+            x=x_vals,
+            y=y_vals_mean,
+            yerr=y_vals_mse,
+            markersize=fr.marker_size,
+            marker="o",
+            mfc=fr.color if ifr == 0 else "none",
+            linestyle="none",
+            label=fr.label + " mean",
+            color=fr.color,
+        )
+
+        if ifr > 0:
+            y_ratios = [
+                y_ratio_values[ifr][ib] / y_ratio_values[0][ib]
+                for ib in range(pr.x_bins)
+            ]
+            axs[1].errorbar(
+                x=x_vals,
+                y=y_ratios,
+                yerr=y_ratio_errors,
+                markersize=fr.marker_size,
+                marker="o",
+                mfc="none",
+                linestyle="none",
+                color=fr.color,
+                label=fr.label,
+            )
+            axs[1].set_ylabel("Ratio to " + fileRecords[0].label)
+
+    # Some final cosmetics
+    axs[0].set_ylim(0.0, y_lim)
+    axs[0].grid(axis="x", linestyle="dotted")
+
+    axs[1].set_ylim(0.9, 1.1)
+    axs[1].grid(axis="x", linestyle="dotted")
+    axs[1].axhline(y=1.0, color="black", linestyle="-")
+
+    axs[0].legend(loc="upper center")
+    axs[1].legend(loc="upper center")
+
+    # Set the range of x-axis
+    plt.xlabel(pr.x_label)
+    plt.show()
+    fig.savefig(pr.saveAs)
diff --git a/Examples/Scripts/Python/full_chain_odd.py b/Examples/Scripts/Python/full_chain_odd.py
index b5012022bbd..eba903e7d2b 100755
--- a/Examples/Scripts/Python/full_chain_odd.py
+++ b/Examples/Scripts/Python/full_chain_odd.py
@@ -407,7 +407,7 @@
     addVertexFitting(
         s,
         field,
-        vertexFinder=VertexFinder.Iterative,
+        vertexFinder=VertexFinder.AMVF,
         outputDirRoot=outputDir if args.output_root else None,
     )
 
diff --git a/Examples/Scripts/Python/truth_tracking_gsf_refitting.py b/Examples/Scripts/Python/truth_tracking_gsf_refitting.py
index 3937143444e..e1a3003fc26 100755
--- a/Examples/Scripts/Python/truth_tracking_gsf_refitting.py
+++ b/Examples/Scripts/Python/truth_tracking_gsf_refitting.py
@@ -9,49 +9,63 @@
 
 u = acts.UnitConstants
 
-srcdir = Path(__file__).resolve().parent.parent.parent.parent
-outputDir = Path.cwd()
 
-# detector, trackingGeometry, _ = getOpenDataDetector()
-detector, trackingGeometry, decorators = acts.examples.GenericDetector.create()
-field = acts.ConstantBField(acts.Vector3(0, 0, 2 * u.T))
-
-s = runTruthTrackingKalman(
+def runRefittingGsf(
     trackingGeometry,
     field,
-    digiConfigFile=srcdir
-    / "Examples/Algorithms/Digitization/share/default-smearing-config-generic.json",
-    # "thirdparty/OpenDataDetector/config/odd-digi-smearing-config.json",
-    outputDir=outputDir,
-)
-
-gsfOptions = {
-    "betheHeitlerApprox": acts.examples.AtlasBetheHeitlerApprox.makeDefault(),
-    "maxComponents": 4,
-    "abortOnError": False,
-    "disableAllMaterialHandling": False,
-    "finalReductionMethod": acts.examples.FinalReductionMethod.maxWeight,
-    "weightCutoff": 1.0e-4,
-    "level": acts.logging.INFO,
-}
-
-s.addAlgorithm(
-    acts.examples.RefittingAlgorithm(
-        acts.logging.INFO,
-        inputTracks="kf_tracks",
-        outputTracks="gsf_tracks",
-        fit=acts.examples.makeGsfFitterFunction(trackingGeometry, field, **gsfOptions),
+    outputDir,
+    s: acts.examples.Sequencer = None,
+):
+    srcdir = Path(__file__).resolve().parent.parent.parent.parent
+
+    s = runTruthTrackingKalman(
+        trackingGeometry,
+        field,
+        digiConfigFile=srcdir
+        / "Examples/Algorithms/Digitization/share/default-smearing-config-generic.json",
+        # "thirdparty/OpenDataDetector/config/odd-digi-smearing-config.json",
+        outputDir=outputDir,
+        s=s,
+    )
+
+    gsfOptions = {
+        "betheHeitlerApprox": acts.examples.AtlasBetheHeitlerApprox.makeDefault(),
+        "maxComponents": 12,
+        "componentMergeMethod": acts.examples.ComponentMergeMethod.maxWeight,
+        "mixtureReductionAlgorithm": acts.examples.MixtureReductionAlgorithm.KLDistance,
+        "weightCutoff": 1.0e-4,
+        "level": acts.logging.INFO,
+    }
+
+    s.addAlgorithm(
+        acts.examples.RefittingAlgorithm(
+            acts.logging.INFO,
+            inputTracks="kf_tracks",
+            outputTracks="gsf_tracks",
+            fit=acts.examples.makeGsfFitterFunction(
+                trackingGeometry, field, **gsfOptions
+            ),
+        )
     )
-)
-
-s.addWriter(
-    acts.examples.TrackFitterPerformanceWriter(
-        level=acts.logging.INFO,
-        inputTracks="tracks",
-        inputParticles="truth_seeds_selected",
-        inputTrackParticleMatching="track_particle_matching",
-        filePath=str(outputDir / "performance_refitter.root"),
+
+    s.addWriter(
+        acts.examples.TrackFitterPerformanceWriter(
+            level=acts.logging.INFO,
+            inputTracks="tracks",
+            inputParticles="truth_seeds_selected",
+            inputTrackParticleMatching="track_particle_matching",
+            filePath=str(outputDir / "performance_refitter.root"),
+        )
     )
-)
 
-s.run()
+    return s
+
+
+if __name__ == "__main__":
+    outputDir = Path.cwd()
+
+    # detector, trackingGeometry, _ = getOpenDataDetector()
+    detector, trackingGeometry, decorators = acts.examples.GenericDetector.create()
+    field = acts.ConstantBField(acts.Vector3(0, 0, 2 * u.T))
+
+    runRefittingGsf(trackingGeometry, field, outputDir).run()
diff --git a/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp b/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp
index dd6b0affd69..991b8bc8f89 100644
--- a/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp
+++ b/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp
@@ -127,8 +127,8 @@ struct SimulationActor {
         (result.properTimeLimit - result.particle.properTime() <
          result.properTimeLimit * properTimeRelativeTolerance)) {
       auto descendants = decay.run(generator, result.particle);
-      for (auto &&descendant : descendants) {
-        result.generatedParticles.emplace_back(std::move(descendant));
+      for (const auto &descendant : descendants) {
+        result.generatedParticles.emplace_back(descendant);
       }
       result.isAlive = false;
       return;
diff --git a/Fatras/src/Digitization/PlanarSurfaceMask.cpp b/Fatras/src/Digitization/PlanarSurfaceMask.cpp
index badf058c611..f660082bfad 100644
--- a/Fatras/src/Digitization/PlanarSurfaceMask.cpp
+++ b/Fatras/src/Digitization/PlanarSurfaceMask.cpp
@@ -38,7 +38,7 @@ namespace {
 /// @param sLength The segment length, maximal allowed length
 void checkIntersection(std::vector<Acts::Intersection2D>& intersections,
                        const Acts::Intersection2D& candidate, double sLength) {
-  if (candidate && candidate.pathLength() > 0 &&
+  if (candidate.isValid() && candidate.pathLength() > 0 &&
       candidate.pathLength() < sLength) {
     intersections.push_back(candidate);
   }
@@ -282,7 +282,7 @@ ActsFatras::PlanarSurfaceMask::annulusMask(const Acts::AnnulusBounds& aBounds,
         Acts::VectorHelpers::phi(vertices[phii[iarc * 2]] - moduleOrigin),
         Acts::VectorHelpers::phi(vertices[phii[iarc * 2 + 1]] - moduleOrigin),
         segment[0] - moduleOrigin, sDir);
-    if (intersection) {
+    if (intersection.isValid()) {
       checkIntersection(intersections,
                         Acts::Intersection2D(
                             intersection.position() + moduleOrigin,
diff --git a/Plugins/Json/src/DetectorVolumeFinderJsonConverter.cpp b/Plugins/Json/src/DetectorVolumeFinderJsonConverter.cpp
index 45d89529c94..9b4d93532f2 100644
--- a/Plugins/Json/src/DetectorVolumeFinderJsonConverter.cpp
+++ b/Plugins/Json/src/DetectorVolumeFinderJsonConverter.cpp
@@ -46,8 +46,8 @@ struct IndexedVolumesGenerator {
             Acts::Experimental::DetectorVolumeFiller>;
 
     auto indexedDetectorVolumeImpl =
-        std::make_unique<const IndexedDetectorVolumesImpl>(std::move(grid), bv,
-                                                           transform);
+        std::make_unique<const IndexedDetectorVolumesImpl>(
+            std::forward<grid_type>(grid), bv, transform);
 
     // Create the delegate and connect it
     Acts::Experimental::ExternalNavigationDelegate vFinder;
diff --git a/Plugins/Json/src/IndexedSurfacesJsonConverter.cpp b/Plugins/Json/src/IndexedSurfacesJsonConverter.cpp
index d53b15dd5cc..a0e027c9bf0 100644
--- a/Plugins/Json/src/IndexedSurfacesJsonConverter.cpp
+++ b/Plugins/Json/src/IndexedSurfacesJsonConverter.cpp
@@ -40,7 +40,7 @@ struct IndexedSurfacesGenerator {
       const std::array<Acts::BinningValue, grid_type::DIM>& bv,
       const Acts::Transform3& transform) {
     Acts::Experimental::IndexedSurfacesNavigation<grid_type> indexedSurfaces(
-        std::move(grid), bv, transform);
+        std::forward<grid_type>(grid), bv, transform);
 
     // The portal delegate
     Acts::Experimental::AllPortalsNavigation allPortals;
diff --git a/Tests/UnitTests/Core/EventData/TrackTests.cpp b/Tests/UnitTests/Core/EventData/TrackTests.cpp
index b6fc3bb1174..9b549630151 100644
--- a/Tests/UnitTests/Core/EventData/TrackTests.cpp
+++ b/Tests/UnitTests/Core/EventData/TrackTests.cpp
@@ -504,4 +504,37 @@ BOOST_AUTO_TEST_CASE(CalculateQuantities) {
   BOOST_CHECK_EQUAL(t.nSharedHits(), 2);
 }
 
+BOOST_AUTO_TEST_CASE(ShallowCopy) {
+  TrackContainer tc{VectorTrackContainer{}, VectorMultiTrajectory{}};
+  auto t = tc.makeTrack();
+
+  auto ts1 = t.appendTrackState();
+  ts1.predicted().setRandom();
+  auto ts2 = t.appendTrackState();
+  ts2.predicted().setRandom();
+  auto ts3 = t.appendTrackState();
+  ts3.predicted().setRandom();
+
+  auto t2 = t.shallowCopy();
+
+  BOOST_CHECK_NE(t.index(), t2.index());
+  BOOST_CHECK_EQUAL(t.nTrackStates(), t2.nTrackStates());
+
+  std::vector<decltype(tc)::TrackStateProxy> trackStates;
+  for (const auto& ts : t.trackStatesReversed()) {
+    trackStates.insert(trackStates.begin(), ts);
+  }
+
+  auto t2_ts1 = trackStates.at(0);
+  auto t2_ts2 = trackStates.at(1);
+  auto t2_ts3 = trackStates.at(2);
+
+  BOOST_CHECK_EQUAL(t2_ts1.predicted(), ts1.predicted());
+  BOOST_CHECK_EQUAL(t2_ts1.index(), ts1.index());
+  BOOST_CHECK_EQUAL(t2_ts2.predicted(), ts2.predicted());
+  BOOST_CHECK_EQUAL(t2_ts2.index(), ts2.index());
+  BOOST_CHECK_EQUAL(t2_ts3.predicted(), ts3.predicted());
+  BOOST_CHECK_EQUAL(t2_ts3.index(), ts3.index());
+}
+
 BOOST_AUTO_TEST_SUITE_END()
diff --git a/Tests/UnitTests/Core/Navigation/DetectorNavigatorTests.cpp b/Tests/UnitTests/Core/Navigation/DetectorNavigatorTests.cpp
index a6a098ad248..d9d8dac01e6 100644
--- a/Tests/UnitTests/Core/Navigation/DetectorNavigatorTests.cpp
+++ b/Tests/UnitTests/Core/Navigation/DetectorNavigatorTests.cpp
@@ -329,7 +329,8 @@ BOOST_AUTO_TEST_CASE(DetectorNavigatorTestsForwardBackward) {
 
   // Step to the volume1|volume2 boundary (portal has witched volume id)
   BOOST_CHECK_EQUAL(statesFwd[2].currentVolume->geometryId(), 2);
-  BOOST_CHECK_EQUAL(statesFwd[2].currentSurface, nullptr);
+  BOOST_CHECK_EQUAL(statesFwd[2].currentSurface,
+                    &(statesFwd[2].currentPortal->surface()));
   BOOST_CHECK_EQUAL(statesFwd[2].currentPortal->surface().geometryId(), 7);
 
   // Step to the surface inside volume2
@@ -339,7 +340,8 @@ BOOST_AUTO_TEST_CASE(DetectorNavigatorTestsForwardBackward) {
 
   // Step to the volume2|volume3 boundary - volume has switched
   BOOST_CHECK_EQUAL(statesFwd[4].currentVolume->geometryId(), 3);
-  BOOST_CHECK_EQUAL(statesFwd[4].currentSurface, nullptr);
+  BOOST_CHECK_EQUAL(statesFwd[4].currentSurface,
+                    &(statesFwd[4].currentPortal->surface()));
   BOOST_CHECK_EQUAL(statesFwd[4].currentPortal->surface().geometryId(), 10);
 
   // Step to the surface inside volume3
@@ -349,7 +351,8 @@ BOOST_AUTO_TEST_CASE(DetectorNavigatorTestsForwardBackward) {
 
   // Step to the volume3|endOfWorld boundary
   BOOST_CHECK_EQUAL(statesFwd[6].currentVolume, nullptr);
-  BOOST_CHECK_EQUAL(statesFwd[6].currentSurface, nullptr);
+  BOOST_CHECK_EQUAL(statesFwd[6].currentSurface,
+                    &(statesFwd[6].currentPortal->surface()));
   BOOST_CHECK_EQUAL(statesFwd[6].currentPortal->surface().geometryId(), 11);
 
   // Step to the end of world
@@ -359,7 +362,8 @@ BOOST_AUTO_TEST_CASE(DetectorNavigatorTestsForwardBackward) {
   // Action list call before the first step
   // Starting in the volume3
   BOOST_CHECK_EQUAL(statesBwd[6].currentVolume, nullptr);
-  BOOST_CHECK_EQUAL(statesBwd[6].currentSurface, nullptr);
+  BOOST_CHECK_EQUAL(statesBwd[6].currentSurface,
+                    &(statesBwd[6].currentPortal->surface()));
   BOOST_CHECK_EQUAL(statesBwd[6].currentPortal->surface().geometryId(), 6);
 
   // Step to the surface inside volume1
@@ -369,7 +373,8 @@ BOOST_AUTO_TEST_CASE(DetectorNavigatorTestsForwardBackward) {
 
   // Step to the volume1|volume2 boundary / preStep not yet set
   BOOST_CHECK_EQUAL(statesBwd[4].currentVolume->geometryId(), 1);
-  BOOST_CHECK_EQUAL(statesBwd[4].currentSurface, nullptr);
+  BOOST_CHECK_EQUAL(statesBwd[4].currentSurface,
+                    &(statesBwd[4].currentPortal->surface()));
   BOOST_CHECK_EQUAL(statesBwd[4].currentPortal->surface().geometryId(), 7);
 
   // Step to the surface inside volume2
@@ -379,7 +384,8 @@ BOOST_AUTO_TEST_CASE(DetectorNavigatorTestsForwardBackward) {
 
   // Step to the volume2|volume3 boundary / pre-step not yet set
   BOOST_CHECK_EQUAL(statesBwd[2].currentVolume->geometryId(), 2);
-  BOOST_CHECK_EQUAL(statesBwd[2].currentSurface, nullptr);
+  BOOST_CHECK_EQUAL(statesBwd[2].currentSurface,
+                    &(statesBwd[2].currentPortal->surface()));
   BOOST_CHECK_EQUAL(statesBwd[2].currentPortal->surface().geometryId(), 10);
   BOOST_CHECK_EQUAL(statesBwd[2].surfaceCandidates.size(), 2u);
 
@@ -483,7 +489,8 @@ BOOST_AUTO_TEST_CASE(DetectorNavigatorTestsAmbiguity) {
 
   // Step to the volume|endOfWorld boundary
   BOOST_CHECK_EQUAL(statesFwd[2].currentVolume, nullptr);
-  BOOST_CHECK_EQUAL(statesFwd[2].currentSurface, nullptr);
+  BOOST_CHECK_EQUAL(statesFwd[2].currentSurface,
+                    &(statesFwd[2].currentPortal->surface()));
   BOOST_CHECK_EQUAL(statesFwd[2].currentPortal->surface().geometryId(), 6);
 
   // Step to the end of world
@@ -492,7 +499,8 @@ BOOST_AUTO_TEST_CASE(DetectorNavigatorTestsAmbiguity) {
 
   // Step to the endOfWorld|volume boundary
   BOOST_CHECK_EQUAL(statesBwd[2].currentVolume, nullptr);
-  BOOST_CHECK_EQUAL(statesBwd[2].currentSurface, nullptr);
+  BOOST_CHECK_EQUAL(statesBwd[2].currentSurface,
+                    &(statesBwd[2].currentPortal->surface()));
   BOOST_CHECK_EQUAL(statesBwd[2].currentPortal->surface().geometryId(), 5);
 
   // Step to the cylindrical surface
@@ -606,7 +614,8 @@ BOOST_AUTO_TEST_CASE(DetectorNavigatorTestsMultipleIntersection) {
 
   // Step to the volume|endOfWorld boundary
   BOOST_CHECK_EQUAL(statesFwd[3].currentVolume, nullptr);
-  BOOST_CHECK_EQUAL(statesFwd[3].currentSurface, nullptr);
+  BOOST_CHECK_EQUAL(statesFwd[3].currentSurface,
+                    &(statesFwd[3].currentPortal->surface()));
   BOOST_CHECK_EQUAL(statesFwd[3].currentPortal->surface().geometryId(), 6);
 
   // Step to the end of world
@@ -615,7 +624,8 @@ BOOST_AUTO_TEST_CASE(DetectorNavigatorTestsMultipleIntersection) {
 
   // Step to the endOfWorld|volume boundary
   BOOST_CHECK_EQUAL(statesBwd[3].currentVolume, nullptr);
-  BOOST_CHECK_EQUAL(statesBwd[3].currentSurface, nullptr);
+  BOOST_CHECK_EQUAL(statesBwd[3].currentSurface,
+                    &(statesBwd[3].currentPortal->surface()));
   BOOST_CHECK_EQUAL(statesBwd[3].currentPortal->surface().geometryId(), 5);
 
   // Second intersection of the cylindrical surface
diff --git a/Tests/UnitTests/Core/Seeding/CMakeLists.txt b/Tests/UnitTests/Core/Seeding/CMakeLists.txt
index c3187405826..5f59cf614ba 100644
--- a/Tests/UnitTests/Core/Seeding/CMakeLists.txt
+++ b/Tests/UnitTests/Core/Seeding/CMakeLists.txt
@@ -5,3 +5,4 @@ add_unittest(EstimateTrackParamsFromSeed EstimateTrackParamsFromSeedTest.cpp)
 add_unittest(BinnedGroupTest BinnedGroupTest.cpp)
 add_unittest(HoughTransformTest HoughTransformTest.cpp)
 add_unittest(UtilityFunctions UtilityFunctionsTests.cpp)
+add_unittest(CandidatesForMiddleSp CandidatesForMiddleSpTests.cpp)
diff --git a/Tests/UnitTests/Core/Seeding/CandidatesForMiddleSpTests.cpp b/Tests/UnitTests/Core/Seeding/CandidatesForMiddleSpTests.cpp
new file mode 100644
index 00000000000..3d8b672cb45
--- /dev/null
+++ b/Tests/UnitTests/Core/Seeding/CandidatesForMiddleSpTests.cpp
@@ -0,0 +1,177 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2024 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include <boost/test/unit_test.hpp>
+
+#include "Acts/Seeding/CandidatesForMiddleSp.hpp"
+
+#include <limits>
+#include <vector>
+
+#include "SpacePoint.hpp"
+
+namespace Acts::Test {
+
+BOOST_AUTO_TEST_CASE(TripletCandidateObject) {
+  using UnitTestSpacePoint = ::SpacePoint;
+  std::vector<UnitTestSpacePoint> spacePoints(3);
+
+  // Default Constructor
+  Acts::TripletCandidate<UnitTestSpacePoint> defaultCandidate;
+  BOOST_CHECK_EQUAL(defaultCandidate.bottom, nullptr);
+  BOOST_CHECK_EQUAL(defaultCandidate.middle, nullptr);
+  BOOST_CHECK_EQUAL(defaultCandidate.top, nullptr);
+  BOOST_CHECK_EQUAL(defaultCandidate.weight, 0.);
+  BOOST_CHECK_EQUAL(defaultCandidate.zOrigin, 0.);
+  BOOST_CHECK_EQUAL(defaultCandidate.isQuality, false);
+
+  // Constructor
+  Acts::TripletCandidate<UnitTestSpacePoint> constructedCandidate(
+      spacePoints[0], spacePoints[1], spacePoints[2], 2.4f, 1.1f, true);
+  BOOST_CHECK_EQUAL(constructedCandidate.bottom, &spacePoints[0]);
+  BOOST_CHECK_EQUAL(constructedCandidate.middle, &spacePoints[1]);
+  BOOST_CHECK_EQUAL(constructedCandidate.top, &spacePoints[2]);
+  BOOST_CHECK_EQUAL(constructedCandidate.weight, 2.4f);
+  BOOST_CHECK_EQUAL(constructedCandidate.zOrigin, 1.1f);
+  BOOST_CHECK_EQUAL(constructedCandidate.isQuality, true);
+
+  // Copy Constructor
+  Acts::TripletCandidate<UnitTestSpacePoint> copiedConstructedCandidate(
+      constructedCandidate);
+  BOOST_CHECK_EQUAL(copiedConstructedCandidate.bottom, &spacePoints[0]);
+  BOOST_CHECK_EQUAL(copiedConstructedCandidate.middle, &spacePoints[1]);
+  BOOST_CHECK_EQUAL(copiedConstructedCandidate.top, &spacePoints[2]);
+  BOOST_CHECK_EQUAL(copiedConstructedCandidate.weight, 2.4f);
+  BOOST_CHECK_EQUAL(copiedConstructedCandidate.zOrigin, 1.1f);
+  BOOST_CHECK_EQUAL(copiedConstructedCandidate.isQuality, true);
+
+  // Copy Assign
+  Acts::TripletCandidate<UnitTestSpacePoint> copiedAssignCandidate =
+      constructedCandidate;
+  BOOST_CHECK_EQUAL(copiedAssignCandidate.bottom, &spacePoints[0]);
+  BOOST_CHECK_EQUAL(copiedAssignCandidate.middle, &spacePoints[1]);
+  BOOST_CHECK_EQUAL(copiedAssignCandidate.top, &spacePoints[2]);
+  BOOST_CHECK_EQUAL(copiedAssignCandidate.weight, 2.4f);
+  BOOST_CHECK_EQUAL(copiedAssignCandidate.zOrigin, 1.1f);
+  BOOST_CHECK_EQUAL(copiedAssignCandidate.isQuality, true);
+
+  // Move Constructor
+  Acts::TripletCandidate<UnitTestSpacePoint> movedConstructedCandidate(
+      std::move(constructedCandidate));
+  BOOST_CHECK_EQUAL(movedConstructedCandidate.bottom, &spacePoints[0]);
+  BOOST_CHECK_EQUAL(movedConstructedCandidate.middle, &spacePoints[1]);
+  BOOST_CHECK_EQUAL(movedConstructedCandidate.top, &spacePoints[2]);
+  BOOST_CHECK_EQUAL(movedConstructedCandidate.weight, 2.4f);
+  BOOST_CHECK_EQUAL(movedConstructedCandidate.zOrigin, 1.1f);
+  BOOST_CHECK_EQUAL(movedConstructedCandidate.isQuality, true);
+  BOOST_CHECK_EQUAL(constructedCandidate.bottom, nullptr);
+  BOOST_CHECK_EQUAL(constructedCandidate.middle, nullptr);
+  BOOST_CHECK_EQUAL(constructedCandidate.top, nullptr);
+
+  // Move Assign
+  Acts::TripletCandidate<UnitTestSpacePoint> movedAssignCandidate =
+      std::move(copiedAssignCandidate);
+  BOOST_CHECK_EQUAL(movedAssignCandidate.bottom, &spacePoints[0]);
+  BOOST_CHECK_EQUAL(movedAssignCandidate.middle, &spacePoints[1]);
+  BOOST_CHECK_EQUAL(movedAssignCandidate.top, &spacePoints[2]);
+  BOOST_CHECK_EQUAL(movedAssignCandidate.weight, 2.4f);
+  BOOST_CHECK_EQUAL(movedAssignCandidate.zOrigin, 1.1f);
+  BOOST_CHECK_EQUAL(movedAssignCandidate.isQuality, true);
+  BOOST_CHECK_EQUAL(copiedAssignCandidate.bottom, nullptr);
+  BOOST_CHECK_EQUAL(copiedAssignCandidate.middle, nullptr);
+  BOOST_CHECK_EQUAL(copiedAssignCandidate.top, nullptr);
+}
+
+BOOST_AUTO_TEST_CASE(CandidatesForMiddleSpObject) {
+  using UnitTestSpacePoint = ::SpacePoint;
+  using value_t =
+      typename Acts::CandidatesForMiddleSp<UnitTestSpacePoint>::value_type;
+  UnitTestSpacePoint spacePoint;
+
+  Acts::CandidatesForMiddleSp<UnitTestSpacePoint> container;
+  container.setMaxElements(std::numeric_limits<std::size_t>::max(),
+                           std::numeric_limits<std::size_t>::max());
+  BOOST_CHECK_EQUAL(container.nLowQualityCandidates(), 0);
+  BOOST_CHECK_EQUAL(container.nHighQualityCandidates(), 0);
+  for (int i(0); i < 20; ++i) {
+    container.push(spacePoint, spacePoint, spacePoint, 1, 2.1, false);
+  }
+  BOOST_CHECK_EQUAL(container.nLowQualityCandidates(), 20);
+  BOOST_CHECK_EQUAL(container.nHighQualityCandidates(), 0);
+  container.clear();
+
+  container.setMaxElements(5, 3);
+  BOOST_CHECK_EQUAL(container.nLowQualityCandidates(), 0);
+  BOOST_CHECK_EQUAL(container.nHighQualityCandidates(), 0);
+
+  std::vector<value_t> emptyStorage = container.storage();
+  BOOST_CHECK_EQUAL(emptyStorage.size(), 0);
+  BOOST_CHECK_EQUAL(container.nLowQualityCandidates(), 0);
+  BOOST_CHECK_EQUAL(container.nHighQualityCandidates(), 0);
+
+  // push low quality
+  for (int i(0); i < 2; ++i) {
+    container.push(spacePoint, spacePoint, spacePoint, i, 2.1, false);
+  }
+  BOOST_CHECK_EQUAL(container.nLowQualityCandidates(), 2);
+  BOOST_CHECK_EQUAL(container.nHighQualityCandidates(), 0);
+
+  for (int i(0); i < 7; ++i) {
+    container.push(spacePoint, spacePoint, spacePoint, 2.01, 2.15, false);
+  }
+  BOOST_CHECK_EQUAL(container.nLowQualityCandidates(), 5);
+  BOOST_CHECK_EQUAL(container.nHighQualityCandidates(), 0);
+
+  // push high quality
+  for (int i(0); i < 5; ++i) {
+    container.push(spacePoint, spacePoint, spacePoint, 0.5f + i, 2.1, true);
+  }
+  BOOST_CHECK_EQUAL(container.nLowQualityCandidates(), 5);
+  BOOST_CHECK_EQUAL(container.nHighQualityCandidates(), 3);
+
+  std::vector<value_t> storagedValues = container.storage();
+  // check size array
+  BOOST_CHECK_EQUAL(storagedValues.size(), 5 + 3);
+  BOOST_CHECK_EQUAL(container.nLowQualityCandidates(), 0);
+  BOOST_CHECK_EQUAL(container.nHighQualityCandidates(), 0);
+
+  // check elements are sorted
+  for (std::size_t i(0); i < storagedValues.size() - 1; ++i) {
+    BOOST_CHECK(storagedValues[i].weight >= storagedValues[i + 1].weight);
+  }
+
+  std::sort(
+      storagedValues.begin(), storagedValues.end(),
+      Acts::CandidatesForMiddleSp<UnitTestSpacePoint>::ascendingByQuality);
+  // check values are sorted properly
+  for (std::size_t i(0); i < storagedValues.size() - 1; ++i) {
+    BOOST_CHECK(storagedValues[i].weight <= storagedValues[i + 1].weight);
+  }
+
+  std::sort(
+      storagedValues.begin(), storagedValues.end(),
+      Acts::CandidatesForMiddleSp<UnitTestSpacePoint>::descendingByQuality);
+  // check values are sorted properly
+  for (std::size_t i(0); i < storagedValues.size() - 1; ++i) {
+    BOOST_CHECK(storagedValues[i].weight >= storagedValues[i + 1].weight);
+  }
+  // push again and check size
+  for (int i(0); i < 7; ++i) {
+    container.push(spacePoint, spacePoint, spacePoint, i, 2.15, false);
+  }
+  for (int i(0); i < 7; ++i) {
+    container.push(spacePoint, spacePoint, spacePoint, i, 2.15, true);
+  }
+  BOOST_CHECK_EQUAL(container.nLowQualityCandidates(), 5);
+  BOOST_CHECK_EQUAL(container.nHighQualityCandidates(), 3);
+  container.clear();
+  BOOST_CHECK_EQUAL(container.nLowQualityCandidates(), 0);
+  BOOST_CHECK_EQUAL(container.nHighQualityCandidates(), 0);
+}
+
+}  // namespace Acts::Test
diff --git a/Tests/UnitTests/Core/Surfaces/DiscSurfaceTests.cpp b/Tests/UnitTests/Core/Surfaces/DiscSurfaceTests.cpp
index c2570b82c86..b50a8cfb81c 100644
--- a/Tests/UnitTests/Core/Surfaces/DiscSurfaceTests.cpp
+++ b/Tests/UnitTests/Core/Surfaces/DiscSurfaceTests.cpp
@@ -210,7 +210,7 @@ BOOST_AUTO_TEST_CASE(DiscSurfaceProperties) {
                             .closest();
   Intersection3D expectedIntersect{Vector3{1.2, 0., 0.}, 10.,
                                    Intersection3D::Status::reachable};
-  BOOST_CHECK(bool(sfIntersection));
+  BOOST_CHECK(sfIntersection.isValid());
   CHECK_CLOSE_ABS(sfIntersection.position(), expectedIntersect.position(),
                   1e-9);
   CHECK_CLOSE_ABS(sfIntersection.pathLength(), expectedIntersect.pathLength(),
diff --git a/Tests/UnitTests/Core/Surfaces/IntersectionHelper2DTests.cpp b/Tests/UnitTests/Core/Surfaces/IntersectionHelper2DTests.cpp
index a76ec58d6ca..67c3d2a3555 100644
--- a/Tests/UnitTests/Core/Surfaces/IntersectionHelper2DTests.cpp
+++ b/Tests/UnitTests/Core/Surfaces/IntersectionHelper2DTests.cpp
@@ -32,8 +32,8 @@ void basicChecks(bool circleCase = false) {
                                 rY, start, direction)
                           : detail::IntersectionHelper2D::intersectEllipse(
                                 rX, rY, start, direction);
-  BOOST_CHECK(!nosol[0]);
-  BOOST_CHECK(!nosol[1]);
+  BOOST_CHECK(!nosol[0].isValid());
+  BOOST_CHECK(!nosol[1].isValid());
 
   start = Vector2(4., -4.);
   auto twosol = circleCase ? detail::IntersectionHelper2D::intersectCircle(
@@ -41,8 +41,8 @@ void basicChecks(bool circleCase = false) {
                            : detail::IntersectionHelper2D::intersectEllipse(
                                  rX, rY, start, direction);
 
-  BOOST_CHECK(twosol[0]);
-  BOOST_CHECK(twosol[1]);
+  BOOST_CHECK(twosol[0].isValid());
+  BOOST_CHECK(twosol[1].isValid());
 
   start = Vector2(-4., 10.);
   direction = Vector2(1., 0.);
@@ -52,9 +52,9 @@ void basicChecks(bool circleCase = false) {
                             : detail::IntersectionHelper2D::intersectEllipse(
                                   rX, rY, start, direction);
 
-  BOOST_CHECK(onesolY[0]);
+  BOOST_CHECK(onesolY[0].isValid());
   CHECK_CLOSE_ABS(onesolY[0].position().x(), 0., s_epsilon);
-  BOOST_CHECK(!onesolY[1]);
+  BOOST_CHECK(!onesolY[1].isValid());
 
   start = Vector2(rX, -4);
   direction = Vector2(0., 1.);
@@ -64,9 +64,9 @@ void basicChecks(bool circleCase = false) {
                             : detail::IntersectionHelper2D::intersectEllipse(
                                   rX, rY, start, direction);
 
-  BOOST_CHECK(onesolX[0]);
+  BOOST_CHECK(onesolX[0].isValid());
   CHECK_CLOSE_ABS(onesolX[0].position().y(), 0., s_epsilon);
-  BOOST_CHECK(!onesolX[1]);
+  BOOST_CHECK(!onesolX[1].isValid());
 }
 
 /// Unit test for creating Ellipse intersection
@@ -78,13 +78,13 @@ BOOST_AUTO_TEST_CASE(LineLineIntersection) {
   auto solution = detail::IntersectionHelper2D::intersectSegment(
       Vector2(5., 3.), Vector2(6., 4), start, dir.normalized());
 
-  BOOST_CHECK(!solution);
+  BOOST_CHECK(!solution.isValid());
 
   // Possible
   solution = detail::IntersectionHelper2D::intersectSegment(
       Vector2(5., 3.), Vector2(3., -1.), start, dir.normalized());
 
-  BOOST_CHECK(solution);
+  BOOST_CHECK(solution.isValid());
 
   // In principle possible, but out of bound
   start = Vector2(2, 3);
@@ -93,12 +93,12 @@ BOOST_AUTO_TEST_CASE(LineLineIntersection) {
   solution = detail::IntersectionHelper2D::intersectSegment(
       Vector2(-1., -2.5), Vector2(3., 2.5), start, dir);
 
-  BOOST_CHECK(solution);
+  BOOST_CHECK(solution.isValid());
 
   solution = detail::IntersectionHelper2D::intersectSegment(
       Vector2(-1., -2.5), Vector2(3., 2.5), start, dir, true);
 
-  BOOST_CHECK(!solution);
+  BOOST_CHECK(!solution.isValid());
 }
 
 /// Unit test for creating Ellipse intersection
diff --git a/Tests/UnitTests/Core/Surfaces/SurfaceIntersectionTests.cpp b/Tests/UnitTests/Core/Surfaces/SurfaceIntersectionTests.cpp
index 8a06a86db04..eb7671e6733 100644
--- a/Tests/UnitTests/Core/Surfaces/SurfaceIntersectionTests.cpp
+++ b/Tests/UnitTests/Core/Surfaces/SurfaceIntersectionTests.cpp
@@ -121,7 +121,7 @@ BOOST_AUTO_TEST_CASE(CylinderIntersectionTests) {
                                               BoundaryTolerance::Infinite());
 
     // Check the validity of the intersection
-    BOOST_CHECK(!iIntersection[0]);
+    BOOST_CHECK(!iIntersection[0].isValid());
 
     // From edge tests - wo boundary test
     auto eIntersection = aCylinder->intersect(tgContext, atEdge, transTZ,
@@ -151,7 +151,7 @@ BOOST_AUTO_TEST_CASE(CylinderIntersectionTests) {
     BOOST_CHECK_EQUAL(eIntersection[0].status(),
                       Intersection3D::Status::reachable);
     // There MUST be a second solution
-    BOOST_CHECK(!eIntersection[1]);
+    BOOST_CHECK(!eIntersection[1].isValid());
     // The other intersection MUST NOT be reachable
     BOOST_CHECK_EQUAL(eIntersection[1].status(),
                       Intersection3D::Status::missed);
@@ -209,7 +209,7 @@ BOOST_AUTO_TEST_CASE(ConeIntersectionTest) {
                                           BoundaryTolerance::Infinite());
 
     // Check the validity of the intersection
-    BOOST_CHECK(!iIntersection[0]);
+    BOOST_CHECK(!iIntersection[0].isValid());
   };
 
   // In a nominal world
@@ -257,7 +257,7 @@ BOOST_AUTO_TEST_CASE(PlanarIntersectionTest) {
     // The path length MUST be positive
     BOOST_CHECK_GT(fIntersection[0].pathLength(), 0.);
     // The intersection MUST be unique
-    BOOST_CHECK(!fIntersection[1]);
+    BOOST_CHECK(!fIntersection[1].isValid());
 
     // On surface intersection
     auto oIntersection = aPlane->intersect(tgContext, onit, direction,
@@ -271,7 +271,7 @@ BOOST_AUTO_TEST_CASE(PlanarIntersectionTest) {
     BOOST_CHECK_LT(std::abs(oIntersection[0].pathLength()),
                    s_onSurfaceTolerance);
     // The intersection MUST be unique
-    BOOST_CHECK(!oIntersection[1]);
+    BOOST_CHECK(!oIntersection[1].isValid());
 
     // Intersect backwards
     auto bIntersection = aPlane->intersect(tgContext, after, direction,
@@ -284,31 +284,31 @@ BOOST_AUTO_TEST_CASE(PlanarIntersectionTest) {
     // The path length MUST be negative
     BOOST_CHECK_LT(bIntersection[0].pathLength(), 0.);
     // The intersection MUST be unique
-    BOOST_CHECK(!bIntersection[1]);
+    BOOST_CHECK(!bIntersection[1].isValid());
 
     // An out of bounds attempt: missed
     auto mIntersection = aPlane->intersect(tgContext, outside, direction,
                                            BoundaryTolerance::None());
     // The intersection MUST NOT be valid
-    BOOST_CHECK(!mIntersection[0]);
+    BOOST_CHECK(!mIntersection[0].isValid());
     // The intersection MUST be reachable
     BOOST_CHECK_EQUAL(mIntersection[0].status(),
                       Intersection3D::Status::missed);
     // The path length MUST be negative
     BOOST_CHECK_GT(mIntersection[0].pathLength(), 0.);
     // The intersection MUST be unique
-    BOOST_CHECK(!mIntersection[1]);
+    BOOST_CHECK(!mIntersection[1].isValid());
 
     // An invalid attempt
     auto iIntersection = aPlane->intersect(tgContext, before, parallel,
                                            BoundaryTolerance::None());
     // The intersection MUST NOT be valid
-    BOOST_CHECK(!iIntersection[0]);
+    BOOST_CHECK(!iIntersection[0].isValid());
     // The intersection MUST be reachable
     BOOST_CHECK_EQUAL(iIntersection[0].status(),
                       Intersection3D::Status::unreachable);
     // The intersection MUST be unique
-    BOOST_CHECK(!iIntersection[1]);
+    BOOST_CHECK(!iIntersection[1].isValid());
   };
 
   // In a nominal world
@@ -355,7 +355,7 @@ BOOST_AUTO_TEST_CASE(LineIntersectionTest) {
     // The path length MUST be positive
     BOOST_CHECK_GT(fIntersection[0].pathLength(), 0.);
     // The intersection MUST be unique
-    BOOST_CHECK(!fIntersection[1]);
+    BOOST_CHECK(!fIntersection[1].isValid());
 
     // On surface intersection - on the straw with random direction
     auto oIntersection = aLine->intersect(tgContext, onit1, direction,
@@ -369,7 +369,7 @@ BOOST_AUTO_TEST_CASE(LineIntersectionTest) {
     BOOST_CHECK_LT(std::abs(oIntersection[0].pathLength()),
                    s_onSurfaceTolerance);
     // The intersection MUST be unique
-    BOOST_CHECK(!oIntersection[1]);
+    BOOST_CHECK(!oIntersection[1].isValid());
 
     // On surface intersecion - on the surface with normal vector
     oIntersection =
@@ -383,7 +383,7 @@ BOOST_AUTO_TEST_CASE(LineIntersectionTest) {
     BOOST_CHECK_LT(std::abs(oIntersection[0].pathLength()),
                    s_onSurfaceTolerance);
     // The intersection MUST be unique
-    BOOST_CHECK(!oIntersection[1]);
+    BOOST_CHECK(!oIntersection[1].isValid());
 
     // Intersect backwards
     auto bIntersection = aLine->intersect(tgContext, after, direction,
@@ -396,31 +396,31 @@ BOOST_AUTO_TEST_CASE(LineIntersectionTest) {
     // The path length MUST be negative
     BOOST_CHECK_LT(bIntersection[0].pathLength(), 0.);
     // The intersection MUST be unique
-    BOOST_CHECK(!bIntersection[1]);
+    BOOST_CHECK(!bIntersection[1].isValid());
 
     // An out of bounds attempt: missed
     auto mIntersection = aLine->intersect(tgContext, outside, direction,
                                           BoundaryTolerance::None());
     // The intersection MUST NOT be valid
-    BOOST_CHECK(!mIntersection[0]);
+    BOOST_CHECK(!mIntersection[0].isValid());
     // The intersection MUST be reachable
     BOOST_CHECK_EQUAL(mIntersection[0].status(),
                       Intersection3D::Status::missed);
     // The path length MUST be negative
     BOOST_CHECK_LT(mIntersection[0].pathLength(), 0.);
     // The intersection MUST be unique
-    BOOST_CHECK(!mIntersection[1]);
+    BOOST_CHECK(!mIntersection[1].isValid());
 
     // An invalid attempt
     auto iIntersection = aLine->intersect(tgContext, before, parallel,
                                           BoundaryTolerance::None());
     // The intersection MUST NOT be valid
-    BOOST_CHECK(!iIntersection[0]);
+    BOOST_CHECK(!iIntersection[0].isValid());
     // The intersection MUST be reachable
     BOOST_CHECK_EQUAL(iIntersection[0].status(),
                       Intersection3D::Status::unreachable);
     // The intersection MUST be unique
-    BOOST_CHECK(!iIntersection[1]);
+    BOOST_CHECK(!iIntersection[1].isValid());
   };
 
   // In a nominal world
diff --git a/Tests/UnitTests/Core/TrackFitting/CMakeLists.txt b/Tests/UnitTests/Core/TrackFitting/CMakeLists.txt
index c415bbb571c..fb18a09b8a7 100644
--- a/Tests/UnitTests/Core/TrackFitting/CMakeLists.txt
+++ b/Tests/UnitTests/Core/TrackFitting/CMakeLists.txt
@@ -5,3 +5,4 @@ add_unittest(Gsf GsfTests.cpp)
 add_unittest(GsfComponentMerging GsfComponentMergingTests.cpp)
 add_unittest(GsfMixtureReduction GsfMixtureReductionTests.cpp)
 add_unittest(Gx2f Gx2fTests.cpp)
+add_unittest(MbfSmoother MbfSmootherTests.cpp)
diff --git a/Tests/UnitTests/Core/TrackFitting/MbfSmootherTests.cpp b/Tests/UnitTests/Core/TrackFitting/MbfSmootherTests.cpp
new file mode 100644
index 00000000000..bdd6d3acab2
--- /dev/null
+++ b/Tests/UnitTests/Core/TrackFitting/MbfSmootherTests.cpp
@@ -0,0 +1,143 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2024 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include <boost/test/unit_test.hpp>
+
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Definitions/TrackParametrization.hpp"
+#include "Acts/EventData/MultiTrajectory.hpp"
+#include "Acts/EventData/TrackStatePropMask.hpp"
+#include "Acts/EventData/TrackStateType.hpp"
+#include "Acts/EventData/VectorMultiTrajectory.hpp"
+#include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
+#include "Acts/TrackFitting/MbfSmoother.hpp"
+#include "Acts/Utilities/Result.hpp"
+#include "Acts/Utilities/detail/Subspace.hpp"
+
+#include <cmath>
+#include <cstddef>
+
+namespace {
+
+using namespace Acts;
+using namespace Acts::Test;
+
+using ParametersVector = Acts::BoundVector;
+using CovarianceMatrix = Acts::BoundSquareMatrix;
+using Jacobian = Acts::BoundMatrix;
+
+const Acts::GeometryContext tgContext;
+
+}  // namespace
+
+BOOST_AUTO_TEST_SUITE(TrackFittingMbfSmoother)
+
+BOOST_AUTO_TEST_CASE(Smooth) {
+  VectorMultiTrajectory traj;
+
+  auto projector = detail::FixedSizeSubspace<eBoundSize, 2>(
+                       std::array{eBoundLoc0, eBoundLoc1})
+                       .projector<double>();
+
+  // Make dummy track parameter
+  CovarianceMatrix covTrk;
+  covTrk.setIdentity();
+  covTrk.diagonal() << 0.08, 0.3, 1, 1, 1, 1;
+  BoundVector parValues;
+
+  std::size_t ts_idx = traj.addTrackState(TrackStatePropMask::All);
+  auto ts = traj.getTrackState(ts_idx);
+  ts.typeFlags().set(TrackStateFlag::MeasurementFlag);
+
+  ts.predicted() << 0.3, 0.5, 0.5 * M_PI, 0., 1 / 100., 0.;
+  ts.predictedCovariance() = covTrk;
+
+  ts.setProjector(projector);
+  ts.allocateCalibrated(2);
+  ts.calibrated<2>() << 0.351, 0.473;
+  ts.calibratedCovariance<2>() << 1e+8, 0., 0., 1e+8;
+
+  ts.filtered() << 0.301, 0.503, 0.5 * M_PI, 0., 1 / 100., 0.;
+  ts.filteredCovariance() = covTrk;
+  ts.pathLength() = 1.;
+  ts.jacobian().setIdentity();
+
+  ts_idx = traj.addTrackState(TrackStatePropMask::All, ts_idx);
+  ts = traj.getTrackState(ts_idx);
+  ts.typeFlags().set(TrackStateFlag::MeasurementFlag);
+
+  ts.predicted() << 0.2, 0.5, 0.5 * M_PI, 0., 1 / 100., 0.;
+  ts.predictedCovariance() = covTrk;
+
+  ts.setProjector(projector);
+  ts.allocateCalibrated(2);
+  ts.calibrated<2>() << 0.351, 0.473;
+  ts.calibratedCovariance<2>() << 1e+8, 0., 0., 1e+8;
+
+  ts.filtered() << 0.27, 0.53, 0.5 * M_PI, 0., 1 / 100., 0.;
+  ts.filteredCovariance() = covTrk;
+  ts.pathLength() = 2.;
+  ts.jacobian().setIdentity();
+
+  ts_idx = traj.addTrackState(TrackStatePropMask::All, ts_idx);
+  ts = traj.getTrackState(ts_idx);
+  ts.typeFlags().set(TrackStateFlag::MeasurementFlag);
+
+  ts.predicted() << 0.35, 0.49, 0.5 * M_PI, 0., 1 / 100., 0.;
+  ts.predictedCovariance() = covTrk;
+
+  ts.setProjector(projector);
+  ts.allocateCalibrated(2);
+  ts.calibrated<2>() << 0.351, 0.473;
+  ts.calibratedCovariance<2>() << 1e+8, 0., 0., 1e+8;
+
+  ts.filtered() << 0.33, 0.43, 0.5 * M_PI, 0., 1 / 100., 0.;
+  ts.filteredCovariance() = covTrk;
+  ts.pathLength() = 3.;
+  ts.jacobian().setIdentity();
+
+  // "smooth" these three track states
+  BOOST_CHECK(MbfSmoother()(tgContext, traj, ts_idx).ok());
+
+  // Regression tests, only tests very basic correctness of the math, but tests
+  // for regressions in the result.
+
+  auto ts1 = traj.getTrackState(0);
+  BOOST_CHECK(ts1.hasSmoothed());
+  BOOST_CHECK_NE(ts1.filtered(), ts1.smoothed());
+
+  double tol = 1e-6;
+
+  ParametersVector expPars;
+  expPars << 0.301, 0.503, 1.5707963, 0.0, 0.01, 0.0;
+  CovarianceMatrix expCov;
+  expCov.setIdentity();
+  expCov.diagonal() << 0.08, 0.3, 1.0, 1.0, 1.0, 1.0;
+  CHECK_CLOSE_ABS(ts1.smoothed(), expPars, tol);
+  CHECK_CLOSE_ABS(ts1.smoothedCovariance(), expCov, tol);
+
+  auto ts2 = traj.getTrackState(1);
+  BOOST_CHECK(ts2.hasSmoothed());
+  BOOST_CHECK_NE(ts2.filtered(), ts2.smoothed());
+
+  expPars << 0.27, 0.53, 1.5707963, 0.0, 0.01, 0.0;
+  CHECK_CLOSE_ABS(ts2.smoothed(), expPars, tol);
+  CHECK_CLOSE_ABS(ts2.smoothedCovariance(), expCov, tol);
+
+  auto ts3 = traj.getTrackState(2);
+  BOOST_CHECK(ts3.hasSmoothed());
+  // last one, smoothed == filtered
+  BOOST_CHECK_EQUAL(ts3.filtered(), ts3.smoothed());
+
+  expPars << 0.33, 0.43, 1.5707963, 0.0, 0.01, 0.0;
+  CHECK_CLOSE_ABS(ts3.smoothed(), expPars, tol);
+  CHECK_CLOSE_ABS(ts3.smoothedCovariance(), expCov, tol);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/Tests/UnitTests/Core/Utilities/IntersectionTests.cpp b/Tests/UnitTests/Core/Utilities/IntersectionTests.cpp
index c2f9afe75e2..076abbfe088 100644
--- a/Tests/UnitTests/Core/Utilities/IntersectionTests.cpp
+++ b/Tests/UnitTests/Core/Utilities/IntersectionTests.cpp
@@ -37,14 +37,14 @@ BOOST_AUTO_TEST_CASE(IntersectionTest) {
                      Intersection3D::Status::reachable);
   Intersection3D tIp(Vector3(0., 3., 0.), 3.,
                      Intersection3D::Status::reachable);
-  BOOST_CHECK(bool(fIp));
-  BOOST_CHECK(bool(sIp));
-  BOOST_CHECK(bool(tIp));
+  BOOST_CHECK(fIp.isValid());
+  BOOST_CHECK(sIp.isValid());
+  BOOST_CHECK(tIp.isValid());
 
   // a non-valid intersection
   Intersection3D nIp(Vector3(3., 3., 0.), 3.,
                      Intersection3D::Status::unreachable);
-  BOOST_CHECK(!bool(nIp));
+  BOOST_CHECK(!nIp.isValid());
 
   std::vector<Intersection3D> fstpIntersections = {fIp, sIp, tIp};
   std::vector<Intersection3D> tsfpIntersections = {tIp, sIp, fIp};
diff --git a/cmake/FindFilesystem.cmake b/cmake/FindFilesystem.cmake
index c42d585c007..c34988752dd 100644
--- a/cmake/FindFilesystem.cmake
+++ b/cmake/FindFilesystem.cmake
@@ -104,7 +104,6 @@ if(TARGET std::filesystem)
     return()
 endif()
 
-cmake_minimum_required(VERSION 3.10)
 
 include(CMakePushCheckState)
 include(CheckIncludeFileCXX)
diff --git a/docs/core/reconstruction/pattern_recognition/seeding.md b/docs/core/reconstruction/pattern_recognition/seeding.md
index 4ee17e98c6d..67fe50c6748 100644
--- a/docs/core/reconstruction/pattern_recognition/seeding.md
+++ b/docs/core/reconstruction/pattern_recognition/seeding.md
@@ -223,7 +223,7 @@ The seed confirmation also sets a limit on the number of seeds produced for each
 which retains only the higher quality seeds. If this limit is exceeded, the algorithm
 checks if there is any low-quality seed in the seed container of this middle SP that can be removed.
 
-:::{doxygenfunction} Acts::SeedFilter::filterSeeds_1SpFixed (Acts::SpacePointData &spacePointData, CandidatesForMiddleSp<const InternalSpacePoint<external_spacepoint_t>>&, const std::size_t, collection_t&) const
+:::{doxygenfunction} Acts::SeedFilter::filterSeeds_1SpFixed (Acts::SpacePointData &spacePointData, CandidatesForMiddleSp<const InternalSpacePoint<external_spacepoint_t>>&, collection_t&) const
 :outline:
 :::
 
diff --git a/thirdparty/FRNN/CMakeLists.txt b/thirdparty/FRNN/CMakeLists.txt
index 30161ebb276..3399ab3289b 100644
--- a/thirdparty/FRNN/CMakeLists.txt
+++ b/thirdparty/FRNN/CMakeLists.txt
@@ -5,8 +5,6 @@
 # This Source Code Form is subject to the terms of the Mozilla Public
 # License, v. 2.0. If a copy of the MPL was not distributed with this
 # file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-cmake_minimum_required(VERSION 3.11)
 include(FetchContent)
 
 message(STATUS "Building FRNN as part of the ACTS project")
@@ -22,5 +20,3 @@ if(NOT frnncontent_POPULATED)
   configure_file(CMakeLists.txt.in "${frnncontent_SOURCE_DIR}/CMakeLists.txt" COPYONLY)
   add_subdirectory(${frnncontent_SOURCE_DIR} ${frnncontent_BINARY_DIR})
 endif()
-
-
diff --git a/thirdparty/GeoModel/CMakeLists.txt b/thirdparty/GeoModel/CMakeLists.txt
index d61f948be2e..0adb3aede35 100644
--- a/thirdparty/GeoModel/CMakeLists.txt
+++ b/thirdparty/GeoModel/CMakeLists.txt
@@ -7,7 +7,6 @@
 # file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 # CMake include(s).
-cmake_minimum_required( VERSION 3.11 )
 include( FetchContent )
 
 # Tell the user what's happening.
diff --git a/thirdparty/actsvg/CMakeLists.txt b/thirdparty/actsvg/CMakeLists.txt
index 2da1b0c7aa5..5b1622d3290 100644
--- a/thirdparty/actsvg/CMakeLists.txt
+++ b/thirdparty/actsvg/CMakeLists.txt
@@ -7,7 +7,6 @@
 # file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 # CMake include(s).
-cmake_minimum_required( VERSION 3.11 )
 include( FetchContent )
 
 # Tell the user what's happening.
diff --git a/thirdparty/algebra-plugins/CMakeLists.txt b/thirdparty/algebra-plugins/CMakeLists.txt
index a0126cece7f..e625715417f 100644
--- a/thirdparty/algebra-plugins/CMakeLists.txt
+++ b/thirdparty/algebra-plugins/CMakeLists.txt
@@ -7,7 +7,6 @@
 # file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 # CMake include(s).
-cmake_minimum_required(VERSION 3.14)
 include(FetchContent)
 
 # Tell the user what's happening.
@@ -44,4 +43,4 @@ set(ALGEBRA_PLUGINS_SETUP_GOOGLETEST FALSE CACHE BOOL
 set(ALGEBRA_PLUGINS_SETUP_BENCHMARK FALSE CACHE BOOL "Do not setup Algebra Plugins benchmark")
 
 # Now set up its build.
-FetchContent_MakeAvailable( AlgebraPlugins )
\ No newline at end of file
+FetchContent_MakeAvailable( AlgebraPlugins )
diff --git a/thirdparty/covfie/CMakeLists.txt b/thirdparty/covfie/CMakeLists.txt
index b74710c587b..bfd6da88215 100644
--- a/thirdparty/covfie/CMakeLists.txt
+++ b/thirdparty/covfie/CMakeLists.txt
@@ -7,7 +7,6 @@
 # file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 # CMake include(s).
-cmake_minimum_required(VERSION 3.14)
 include(FetchContent)
 
 # Tell the user what's happening.
@@ -30,4 +29,4 @@ set( COVFIE_REQUIRE_CXX20 OFF CACHE BOOL "Enable covfie C++20 requirement")
 set( COVFIE_QUIET ON CACHE BOOL "Quiet covfie feature warnings")
 
 # Now set up its build.
-FetchContent_MakeAvailable( covfie )
\ No newline at end of file
+FetchContent_MakeAvailable( covfie )
diff --git a/thirdparty/detray/CMakeLists.txt b/thirdparty/detray/CMakeLists.txt
index 961057df933..8ec164e2440 100644
--- a/thirdparty/detray/CMakeLists.txt
+++ b/thirdparty/detray/CMakeLists.txt
@@ -7,7 +7,6 @@
 # file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 # CMake include(s).
-cmake_minimum_required( VERSION 3.11 )
 include( FetchContent )
 
 # Tell the user what's happening.
@@ -56,4 +55,4 @@ set( DETRAY_SETUP_DFELIBS OFF CACHE BOOL
 
 
 #Now set up its build.
-FetchContent_MakeAvailable(Detray)
\ No newline at end of file
+FetchContent_MakeAvailable(Detray)
diff --git a/thirdparty/dfelibs/CMakeLists.txt b/thirdparty/dfelibs/CMakeLists.txt
index d95941b8406..c3dc15c53d0 100644
--- a/thirdparty/dfelibs/CMakeLists.txt
+++ b/thirdparty/dfelibs/CMakeLists.txt
@@ -1,4 +1,3 @@
-cmake_minimum_required( VERSION 3.11 )
 include( FetchContent )
 
 # Tell the user what's happening.
diff --git a/thirdparty/nlohmann_json/CMakeLists.txt b/thirdparty/nlohmann_json/CMakeLists.txt
index 39aedcb2d5c..f54653556b0 100644
--- a/thirdparty/nlohmann_json/CMakeLists.txt
+++ b/thirdparty/nlohmann_json/CMakeLists.txt
@@ -1,13 +1,8 @@
-cmake_minimum_required( VERSION 3.11 )
 include( FetchContent )
 
 # Tell the user what's happening.
 message( STATUS "Building nlohmann_json as part of the ACTS project" )
 
-if(POLICY CMP0135)
-   cmake_policy(SET CMP0135 NEW)
-endif()
-
 # Declare where to get nlohmann json from.
 FetchContent_Declare( nlohmann_json ${ACTS_NLOHMANNJSON_SOURCE} )
 
diff --git a/thirdparty/pybind11/CMakeLists.txt b/thirdparty/pybind11/CMakeLists.txt
index 60b95573dfc..ac31675e057 100644
--- a/thirdparty/pybind11/CMakeLists.txt
+++ b/thirdparty/pybind11/CMakeLists.txt
@@ -1,4 +1,3 @@
-cmake_minimum_required( VERSION 3.11 )
 include( FetchContent )
 
 # Tell the user what's happening.
diff --git a/thirdparty/vecmem/CMakeLists.txt b/thirdparty/vecmem/CMakeLists.txt
index b6d74615471..fca3f49572f 100644
--- a/thirdparty/vecmem/CMakeLists.txt
+++ b/thirdparty/vecmem/CMakeLists.txt
@@ -7,7 +7,6 @@
 # file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 # CMake include(s).
-cmake_minimum_required( VERSION 3.11 )
 include( FetchContent )
 
 # Tell the user what's happening.
@@ -22,7 +21,7 @@ FetchContent_Declare( VecMem ${ACTS_VECMEM_SOURCE} )
 set( VECMEM_BUILD_TESTING FALSE CACHE BOOL
    "Turn off the build of the VecMem unit tests" )
 set( VECMEM_BUILD_BENCHMARKING FALSE CACHE BOOL
-   "Turn off the build of VecMem benchmarking" )   
+   "Turn off the build of VecMem benchmarking" )
 
 # Now set up its build.
 FetchContent_MakeAvailable( VecMem )