diff --git a/Phys/VertexFit/CMakeLists.txt b/Phys/VertexFit/CMakeLists.txt
index 7acec814d00931cfd9341d4ba712c71855e29c43..0bd023a3be640bf4fce5311e64b8a07479fbdd05 100644
--- a/Phys/VertexFit/CMakeLists.txt
+++ b/Phys/VertexFit/CMakeLists.txt
@@ -53,6 +53,7 @@ gaudi_add_module(VertexFit
         Rec::TrackKernel
         ROOT::MathCore
         VertexFitLib
+        Rec::SelToolsLib
 )
 
 gaudi_add_tests(QMTest)
diff --git a/Phys/VertexFit/src/ParticleVertexFitter.cpp b/Phys/VertexFit/src/ParticleVertexFitter.cpp
index 9cc8ea539fd4bec591f3401d6cca7ee928ab0f5a..3b13aee4961bf1f2dfdc7a9edafff2ba878bcb40 100644
--- a/Phys/VertexFit/src/ParticleVertexFitter.cpp
+++ b/Phys/VertexFit/src/ParticleVertexFitter.cpp
@@ -20,6 +20,7 @@
 #include "Kernel/ParticleProperty.h"
 #include "LHCbMath/MatrixTransforms.h"
 #include "MassConstrainer.h"
+#include "SelTools/DistanceCalculator.h"
 #include "TrackInterfaces/ITrackStateProvider.h"
 #include "TrackKernel/TrackTraj.h"
 #include <boost/container/static_vector.hpp>
@@ -99,7 +100,6 @@ private:
 private:
   Gaudi::Property<int>                      m_maxnumiter{ this, "MaxNumIter", 5 };
   Gaudi::Property<double>                   m_maxdchisq{ this, "MaxDeltaChi2", 0.01 };
-  Gaudi::Property<bool>                     m_extrapolateTtracks{ this, "extrapolateTtracks", false };
   PublicToolHandle<ITrackStateProvider>     m_stateprovider{ this, "StateProvider", "TrackStateProvider" };
   ToolHandle<ITrajPoca>                     m_trajpoca{ "TrajPoca" };
   const LHCb::IParticlePropertySvc*         m_ppsvc{ nullptr };
@@ -678,6 +678,10 @@ StatusCode ParticleVertexFitter::fit( const LHCb::Particle::ConstVector& origdau
   // 3. if not, try with downstream tracks and trajpoca
   bool posinitialized{ false };
 
+  const bool ttracks_only = std::all_of( daughters.begin(), daughters.end(), []( const LHCb::Particle* p ) {
+    return p && p->proto() && p->proto()->track() && !p->proto()->track()->hasVelo() && !p->proto()->track()->hasUT();
+  } );
+
   Gaudi::XYZPoint position;
   if ( counttypes[Resonance] > 0 ) {
     // try 1
@@ -717,43 +721,19 @@ StatusCode ParticleVertexFitter::fit( const LHCb::Particle::ConstVector& origdau
       }
     }
     bool parallelTracks = false;
-    if ( ntrajs == 2 && m_extrapolateTtracks ) {
-      // this is NOT default behaviour
-      // For T-tracks it does not make sense to use 3-D POCA starting position
-      // Trajectory in YZ plane is approx straight line
-      // calculate the yz intersection z and extrapolate there
-
-      const float ty0 = daughters[0]->proto()->track()->firstState().ty();
-      const float ty1 = daughters[1]->proto()->track()->firstState().ty();
-
-      if ( LHCb::essentiallyEqual( ty0, ty1 ) ) {
-        debug() << "tracks are parallel in yz plane. Cannot calculate yz intersection. Default to 3-D POCA starting "
-                   "position. "
-                << endmsg;
-        parallelTracks = true;
-      } else {
-        const float z0 = daughters[0]->proto()->track()->firstState().z();
-        const float z1 = daughters[1]->proto()->track()->firstState().z();
-
-        const float y0 = daughters[0]->proto()->track()->firstState().y();
-        const float y1 = daughters[1]->proto()->track()->firstState().y();
-
-        // calculate the y-axis intersection of each track
-        const float c0 = y0 - ty0 * z0;
-        const float c1 = y1 - ty1 * z1;
-
-        // calculate the point in the yz-plane where the tracks intersect
-        const float yzIntersectionZ = ( c1 - c0 ) / ( ty0 - ty1 );
-        const float yzIntersectionY = ty0 * yzIntersectionZ + c1;
-
-        position.SetZ( yzIntersectionZ );
-        position.SetY( yzIntersectionY );
-        position.SetX( 0 );
+    if ( ntrajs == 2 && ttracks_only ) {
+      // Calculate the yz intersection and use it as a seed for POCA
+      // If tracks are essentially parallel in y-z plane (|dty| < m_ttracks_poca_min_dty), use a constant starting point
+      // (m_ttracks_poca_initial_z)
+
+      const auto poca = Sel::NonlinearDistanceCalculator::findPOCA( *( trajs[0] ), *( trajs[1] ), *m_trajpoca );
+      if ( poca.has_value() ) {
+        position       = poca->poca;
         posinitialized = true;
       }
     }
 
-    if ( ( ntrajs == 2 && !m_extrapolateTtracks ) || ( ntrajs == 2 && parallelTracks ) ) {
+    if ( ( ntrajs == 2 && !ttracks_only ) || ( ntrajs == 2 && parallelTracks ) ) {
       // this is the default behaviour
       double           mu0( 0 ), mu1( 0 );
       Gaudi::XYZVector deltaX;
@@ -799,7 +779,7 @@ StatusCode ParticleVertexFitter::fit( const LHCb::Particle::ConstVector& origdau
       }
     } break;
     case TrackWithoutVelo: {
-      if ( !m_extrapolateTtracks ) {
+      if ( !ttracks_only ) {
         // this is the default behaviour
         const LHCb::Track*     track     = p->proto()->track();
         const LHCb::TrackTraj* tracktraj = m_stateprovider->trajectory( *( p->proto()->track() ), geometry );
diff --git a/Tr/TrackExtrapolators/src/TrackStateProvider.cpp b/Tr/TrackExtrapolators/src/TrackStateProvider.cpp
index 209d391cfb5be422683bbc356c67e81672d825b4..06d2442ba24063c6c513e23364c81326290816e9 100644
--- a/Tr/TrackExtrapolators/src/TrackStateProvider.cpp
+++ b/Tr/TrackExtrapolators/src/TrackStateProvider.cpp
@@ -181,14 +181,11 @@ private:
 
   /// Get the track cache from the event store
   TrackCaches& trackcache() const {
-    // auto* obj = m_caches.getIfExists();
-    // return const_cast<TrackCaches&>( obj ? *obj : *( m_caches.put( TrackCaches{} ) ) );
-    static constexpr auto loc = "TrackStateProviderCache";
-    using CacheTES            = AnyDataWrapper<TrackCaches>;
-    auto d                    = getIfExists<CacheTES>( loc );
+    using CacheTES = AnyDataWrapper<TrackCaches>;
+    auto d         = getIfExists<CacheTES>( m_cachelocation );
     if ( !d ) {
       d = new CacheTES( TrackCaches{} );
-      put( d, loc );
+      put( d, m_cachelocation );
     }
     return d->getData();
   }
@@ -203,6 +200,10 @@ private:
   // mutable DataObjectHandle<AnyDataWrapper<TrackCaches>> m_caches{this, Gaudi::DataHandle::Writer, "CacheLocation",
   //                                                               "TrackStateProviderCache"};
 
+  // Depending on TrackStateProvider configuration, the resulting track extrapolation might differ.
+  // In this case we can't simply reuse the cache created by standard configuration (as it's stored in TES).
+  std::string m_cachelocation = std::string( "TrackStateProviderCache_" ) + name();
+
   ToolHandle<ITrackExtrapolator> m_extrapolator{ "TrackMasterExtrapolator", this };
   ToolHandle<ITrackInterpolator> m_interpolator{ "TrackInterpolator", this };
 
@@ -382,6 +383,17 @@ TrackCache TrackStateProvider::createCacheEntry( TkCacheKey key, const LHCb::Tra
     }
   }
 
+  // the same for T-Tracks
+  if ( track.type() == LHCb::Track::Types::Ttrack ) {
+    for ( const auto& loc : { std::pair{ LHCb::State::Location::EndMagnet, StateParameters::ZEndMag },
+                              std::pair{ LHCb::State::Location::MidMagnet, StateParameters::ZMidMag },
+                              std::pair{ LHCb::State::Location::BegMagnet, StateParameters::ZBegMag },
+                              std::pair{ LHCb::State::Location::EndUT, StateParameters::ZEndUT },
+                              std::pair{ LHCb::State::Location::EndVelo, StateParameters::ZEndVelo } } ) {
+      if ( !track.stateAt( loc.first ) ) { addState( tc, loc.second, geometry, loc.first ); }
+    }
+  }
+
   // make sure all tracks (incl. Downstream) get assigned a state at
   // the beamline. this is useful for the trajectory approximation.
   if ( ( track.hasVelo() || track.hasUT() ) && track.firstState().location() != LHCb::State::Location::ClosestToBeam &&