refactor: Assert to be on surface for surface functions with free p…

…aram input (acts-project#2932)

This popped up in acts-project#2811 - as this can be a pitfall we are asserting if we are really on surface. This is done via an `assert` which is not evaluated in standard release build

This also fixes a bug in the covariance engine. We used to first do the cov transport and then check if we are actually on surface. I turned these checks around.

Core/src/Propagator/detail/CovarianceEngine.cpp

@@ -41,6 +41,13 @@ Result<BoundState> detail::boundState(
     const ParticleHypothesis& particleHypothesis, bool covTransport,
     double accumulatedPath,
     const FreeToBoundCorrection& freeToBoundCorrection) {
+  // Create the bound parameters
+  Result<BoundVector> bv = detail::transformFreeToBoundParameters(
+      freeParameters, surface, geoContext);
+  if (!bv.ok()) {
+    return bv.error();
+  }
+
   // Covariance transport
   std::optional<BoundSquareMatrix> cov = std::nullopt;
   if (covTransport) {
@@ -58,12 +65,6 @@ Result<BoundState> detail::boundState(
     cov = boundCovariance;
   }
 
-  // Create the bound parameters
-  Result<BoundVector> bv = detail::transformFreeToBoundParameters(
-      freeParameters, surface, geoContext);
-  if (!bv.ok()) {
-    return bv.error();
-  }
   // Create the bound state
   return std::make_tuple(
       BoundTrackParameters(surface.getSharedPtr(), *bv, std::move(cov),

Core/src/Surfaces/ConeSurface.cpp

@@ -336,6 +336,9 @@ Acts::AlignmentToPathMatrix Acts::ConeSurface::alignmentToPathDerivative(
   const auto position = parameters.segment<3>(eFreePos0);
   // The direction
   const auto direction = parameters.segment<3>(eFreeDir0);
+
+  assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
+
   // The vector between position and center
   const auto pcRowVec = (position - center(gctx)).transpose().eval();
   // The rotation

Core/src/Surfaces/CylinderSurface.cpp

@@ -291,6 +291,9 @@ Acts::AlignmentToPathMatrix Acts::CylinderSurface::alignmentToPathDerivative(
   const auto position = parameters.segment<3>(eFreePos0);
   // The direction
   const auto direction = parameters.segment<3>(eFreeDir0);
+
+  assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
+
   // The vector between position and center
   const auto pcRowVec = (position - center(gctx)).transpose().eval();
   // The rotation

Core/src/Surfaces/DiscSurface.cpp

@@ -219,6 +219,9 @@ Acts::BoundToFreeMatrix Acts::DiscSurface::boundToFreeJacobian(
   const Vector3 position = freeParams.segment<3>(eFreePos0);
   // The direction
   const Vector3 direction = freeParams.segment<3>(eFreeDir0);
+
+  assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
+
   // special polar coordinates for the Disc
   double lrad = boundParams[eBoundLoc0];
   double lphi = boundParams[eBoundLoc1];

Core/src/Surfaces/LineSurface.cpp

@@ -208,6 +208,8 @@ Acts::BoundToFreeMatrix Acts::LineSurface::boundToFreeJacobian(
   // retrieve the reference frame
   auto rframe = referenceFrame(gctx, position, direction);
 
+  assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
+
   // Initialize the jacobian from local to global
   BoundToFreeMatrix jacToGlobal = BoundToFreeMatrix::Zero();
 
@@ -251,6 +253,9 @@ Acts::FreeToPathMatrix Acts::LineSurface::freeToPathDerivative(
   Vector3 position = parameters.segment<3>(eFreePos0);
   // The direction
   Vector3 direction = parameters.segment<3>(eFreeDir0);
+
+  assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
+
   // The vector between position and center
   Vector3 pcRowVec = position - center(gctx);
   // The local frame z axis
@@ -281,6 +286,9 @@ Acts::AlignmentToPathMatrix Acts::LineSurface::alignmentToPathDerivative(
   Vector3 position = parameters.segment<3>(eFreePos0);
   // The direction
   Vector3 direction = parameters.segment<3>(eFreeDir0);
+
+  assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
+
   // The vector between position and center
   Vector3 pcRowVec = position - center(gctx);
   // The local frame z axis

Core/src/Surfaces/Surface.cpp

@@ -59,6 +59,15 @@ bool Acts::Surface::isOnSurface(const GeometryContext& gctx,
 Acts::AlignmentToBoundMatrix Acts::Surface::alignmentToBoundDerivative(
     const GeometryContext& gctx, const FreeVector& parameters,
     const FreeVector& pathDerivative) const {
+  // The global posiiton
+  const auto position = parameters.segment<3>(eFreePos0);
+  // The direction
+  const auto direction = parameters.segment<3>(eFreeDir0);
+
+  (void)position;
+  (void)direction;
+  assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
+
   // 1) Calculate the derivative of bound parameter local position w.r.t.
   // alignment parameters without path length correction
   const auto alignToBoundWithoutCorrection =
@@ -81,6 +90,12 @@ Acts::Surface::alignmentToBoundDerivativeWithoutCorrection(
     const GeometryContext& gctx, const FreeVector& parameters) const {
   // The global posiiton
   const auto position = parameters.segment<3>(eFreePos0);
+  // The direction
+  const auto direction = parameters.segment<3>(eFreeDir0);
+
+  (void)direction;
+  assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
+
   // The vector between position and center
   const auto pcRowVec = (position - center(gctx)).transpose().eval();
   // The local frame rotation
@@ -122,6 +137,9 @@ Acts::AlignmentToPathMatrix Acts::Surface::alignmentToPathDerivative(
   const auto position = parameters.segment<3>(eFreePos0);
   // The direction
   const auto direction = parameters.segment<3>(eFreeDir0);
+
+  assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
+
   // The vector between position and center
   const auto pcRowVec = (position - center(gctx)).transpose().eval();
   // The local frame rotation
@@ -264,6 +282,9 @@ Acts::BoundToFreeMatrix Acts::Surface::boundToFreeJacobian(
   const Vector3 direction = freeParams.segment<3>(eFreeDir0);
   // retrieve the reference frame
   const auto rframe = referenceFrame(gctx, position, direction);
+
+  assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
+
   // Initialize the jacobian from local to global
   BoundToFreeMatrix jacToGlobal = BoundToFreeMatrix::Zero();
   // the local error components - given by reference frame
@@ -286,6 +307,9 @@ Acts::FreeToBoundMatrix Acts::Surface::freeToBoundJacobian(
   // The measurement frame of the surface
   RotationMatrix3 rframeT =
       referenceFrame(gctx, position, direction).transpose();
+
+  assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
+
   // Initialize the jacobian from global to local
   FreeToBoundMatrix jacToLocal = FreeToBoundMatrix::Zero();
   // Local position component given by the reference frame
@@ -305,6 +329,9 @@ Acts::FreeToPathMatrix Acts::Surface::freeToPathDerivative(
   const auto position = parameters.segment<3>(eFreePos0);
   // The direction
   const auto direction = parameters.segment<3>(eFreeDir0);
+
+  assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
+
   // The measurement frame of the surface
   const RotationMatrix3 rframe = referenceFrame(gctx, position, direction);
   // The measurement frame z axis

Tests/UnitTests/Core/EventData/CorrectedTransformFreeToBoundTests.cpp

@@ -77,8 +77,7 @@ BOOST_AUTO_TEST_CASE(CorrectedFreeToBoundTrackParameters) {
 
   // the intersection of the track with the end surface
   SurfaceIntersection intersection =
-      eSurface
-          ->intersect(geoCtx, Vector3(0, 0, 0), dir, Acts::BoundaryCheck(true))
+      eSurface->intersect(geoCtx, Vector3(0, 0, 0), dir, BoundaryCheck(true))
           .closest();
   Vector3 tpos = intersection.position();
   auto s = intersection.pathLength();
@@ -87,9 +86,9 @@ BOOST_AUTO_TEST_CASE(CorrectedFreeToBoundTrackParameters) {
 
   // construct the free parameters vector
   FreeVector eFreeParams = FreeVector::Zero();
-  eFreeParams.segment<3>(0) = tpos;
+  eFreeParams.segment<3>(eFreePos0) = tpos;
   eFreeParams[eFreeTime] = t;
-  eFreeParams.segment<3>(4) = dir;
+  eFreeParams.segment<3>(eFreeDir0) = dir;
   eFreeParams[eFreeQOverP] = qOverP;
 
   // the jacobian from local to global at the starting position

Tests/UnitTests/Core/Propagator/EigenStepperTests.cpp

@@ -505,17 +505,14 @@ BOOST_AUTO_TEST_CASE(eigen_stepper_test) {
   // Update in context of a surface
   freeParams = detail::transformBoundToFreeParameters(
       bp.referenceSurface(), tgContext, bp.parameters());
-  freeParams.segment<3>(eFreePos0) *= 2;
-  freeParams[eFreeTime] *= 2;
-  freeParams[eFreeQOverP] *= -0.5;
 
   es.update(esState, freeParams, bp.parameters(), 2 * (*bp.covariance()),
             *plane);
-  CHECK_CLOSE_OR_SMALL(es.position(esState), 2. * pos, eps, eps);
+  CHECK_CLOSE_OR_SMALL(es.position(esState), pos, eps, eps);
   CHECK_CLOSE_OR_SMALL(es.direction(esState), dir, eps, eps);
-  CHECK_CLOSE_REL(es.absoluteMomentum(esState), 2 * absMom, eps);
-  BOOST_CHECK_EQUAL(es.charge(esState), -1. * charge);
-  CHECK_CLOSE_OR_SMALL(es.time(esState), 2. * time, eps, eps);
+  CHECK_CLOSE_REL(es.absoluteMomentum(esState), absMom, eps);
+  BOOST_CHECK_EQUAL(es.charge(esState), charge);
+  CHECK_CLOSE_OR_SMALL(es.time(esState), time, eps, eps);
   CHECK_CLOSE_COVARIANCE(esState.cov, Covariance(2. * cov), eps);
 
   // Test a case where no step size adjustment is required

Tests/UnitTests/Core/Propagator/StraightLineStepperTests.cpp

@@ -382,16 +382,15 @@ BOOST_AUTO_TEST_CASE(straight_line_stepper_test) {
   // Update in context of a surface
   freeParams = detail::transformBoundToFreeParameters(
       bp.referenceSurface(), tgContext, bp.parameters());
-  freeParams.segment<3>(eFreePos0) *= 2;
-  freeParams[eFreeTime] *= 2;
 
   BOOST_CHECK(bp.covariance().has_value());
   sls.update(slsState, freeParams, bp.parameters(), 2 * (*bp.covariance()),
              *plane);
-  CHECK_CLOSE_OR_SMALL(sls.position(slsState), 2. * pos, eps, eps);
-  BOOST_CHECK_EQUAL(sls.charge(slsState), 1. * charge);
-  CHECK_CLOSE_OR_SMALL(sls.time(slsState), 2. * time, eps, eps);
+  CHECK_CLOSE_OR_SMALL(sls.position(slsState), pos, eps, eps);
+  BOOST_CHECK_EQUAL(sls.charge(slsState), charge);
+  CHECK_CLOSE_OR_SMALL(sls.time(slsState), time, eps, eps);
   CHECK_CLOSE_COVARIANCE(slsState.cov, Covariance(2. * cov), 1e-6);
 }
+
 }  // namespace Test
 }  // namespace Acts