diff --git a/Phys/SelTools/include/SelTools/DistanceCalculator.h b/Phys/SelTools/include/SelTools/DistanceCalculator.h
index ebb186e1b60896dba5210d9e27429a43720e539e..dccce9a8b86c29537c30f482e61a45df4bbf5d6d 100644
--- a/Phys/SelTools/include/SelTools/DistanceCalculator.h
+++ b/Phys/SelTools/include/SelTools/DistanceCalculator.h
@@ -319,10 +319,22 @@ namespace Sel {
       for ( size_t row = 0; row < 3; ++row ) {
         for ( size_t col = 0; col <= row; ++col ) {
           W( row, col ) = vertexCov( row, col ) + pCov( row, col ) + pCov( row + 3, col + 3 ) * a * a -
-                          ( pCov( row, col + 3 ) + pCov( col + 3, row ) ) * a;
+                          ( pCov( row + 3, col ) + pCov( col + 3, row ) ) * a;
         }
       }
 
+      // WH July 2025: This is a slightly more efficient method to compute the same thing. It was used to find a bug in
+      // the lines above.
+      // const auto posCov = posCovMatrix( particle );
+      // const auto momCov = momCovMatrix( particle );
+      // const auto momPosCov = momPosCovMatrix( particle );
+      // for ( size_t row = 0; row < 3; ++row ) {
+      //  for ( size_t col = 0; col <= row; ++col ) {
+      //    W( row, col ) = vertexCov( row, col ) + posCov( row, col ) + momCov( row, col) * a * a -
+      //    ( momPosCov( row, col) + momPosCov( col, row ) ) * a;
+      //  }
+      // }
+
       auto [success, W_Inv] = W.invChol();
       const auto nfailures  = popcount( !success );
       if ( m_dls_failure.has_value() && nfailures > 0 ) *m_dls_failure += nfailures;
@@ -331,11 +343,11 @@ namespace Sel {
         for ( size_t row = 0; row < 3; ++row ) {
           for ( size_t col = 0; col <= row; ++col ) { W( row, col ) = vertexCov( row, col ) + pCov( row, col ); }
         }
+        std::tie( success, W_Inv ) = W.invChol();
       }
-      std::tie( success, W_Inv ) = W.invChol();
-      auto halfdChi2dLam2        = similarity( dir, W_Inv );
-      auto decayLength           = dot( dir, W_Inv * flight ) / halfdChi2dLam2;
-      auto decayLengthErr        = sqrt( 1 / halfdChi2dLam2 );
+      auto halfdChi2dLam2 = similarity( dir, W_Inv );
+      auto decayLength    = dot( dir, W_Inv * flight ) / halfdChi2dLam2;
+      auto decayLengthErr = sqrt( 1 / halfdChi2dLam2 );
       if constexpr ( std::is_arithmetic_v<decltype( success )> ) {
         return success ? decayLength / decayLengthErr : NonPhysicalValue;
       } else {