refactor: New DetectorNavigator tests and some fixes (acts-project#2988)

Adding new tests for the `DetectorNavigator` and introducing some fixes.
 
In the `DetectorNavigator` itself:

1) Adding `DetectorVolume`, `CandidateSurfaces` and `NavigationState`
initialization into `initialize`. Previously, if you ran with the
`endOfWorld` aborter, Propagation didn't enter the stepping loop,
because the current volume was `nullptr` and there were no candidates to
try.

2) Adding implicit `endOfWorld` check after the `Portal` update.
Previously, the `initialize` was called and volume was set to the one
the Navigator is currently exiting (volume finders treat the boundary as
a part of the volume). This resulted in the stepping loop continuing
until the limit in the number of steps is reached.

3) Adding navigation breaks. Some of the aborters (related to the KF and
its tests) were not exiting properly before.

4) Removed `initialize` call in the post step to insure the proper
abortion.

5) Added iteration over the intersection list to store candidates.
Previously, if the surface was intersected more than once, only one of
the intersections was taken into account, even if more of them were
valid.

6) Took direction (e.g. `Acts::Forward`) into account during the
`NavigationState` filling. Previously, the `Backward` direction
propagation didn't work.

7) Small changes in the material-related files to resolve the
`currentVolume` conflicts in the compilation.

In the tests:

1) Added initialization tests, analogous to the ones for the `Navigator`

2) Navigation through a simple cubic geometry in the `Forward` and
`Backward` directions. Consistency checks, and things like that.

3) Tests that ensure the proper behaviour when the multiple valid
intersections are encountered. Like the double crossing of the
Cylindrical surface.

---------

Co-authored-by: Andreas Stefl <stefl.andreas@gmail.com>
Co-authored-by: Andreas Salzburger <asalzburger@gmail.com>

AUTHORS

@@ -20,3 +20,4 @@ The following people have contributed to the project (in alphabetical order):
 - Bastian Schlag, CERN, JGU Mainz
 - Andreas Stefl, CERN, TU Wien
 - Stephen Nicholas Swatman, CERN, University of Amsterdam
+- Roman Urmanov, Weizmann Institute of Science

Core/include/Acts/Detector/GeometryCompatibilityConcept.hpp

@@ -0,0 +1,53 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2023 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/Detector/DetectorVolume.hpp"
+#include "Acts/Geometry/TrackingVolume.hpp"
+#include "Acts/Utilities/TypeTraits.hpp"
+
+namespace Acts {
+
+namespace Concepts {
+
+// Types to check compatibility of
+template <typename propagator_state_t, typename navigator_t>
+using ReturnTypeCurrentVolume =
+    decltype(std::declval<navigator_t>().currentVolume(
+        std::declval<propagator_state_t>().navigation));
+
+/// @brief Concept ensuring compatibility TrackingGeometry
+/// and Detector navigation interfaces with the client code
+/// @tparam propagator_state_t Type of the object for navigation state
+/// @tparam navigator_t Type of the navigator object
+template <typename propagator_state_t, typename navigator_t>
+struct NavigationCompatibilityConceptImpl {
+  /// @brief Ensure that the currentVolume method
+  /// returns one of the known volume types
+  constexpr static bool isCurrentVolumePtr =
+      (Acts::Concepts::identical_to<const TrackingVolume*,
+                                    ReturnTypeCurrentVolume, propagator_state_t,
+                                    navigator_t> ||
+       Acts::Concepts::identical_to<const Acts::Experimental::DetectorVolume*,
+                                    ReturnTypeCurrentVolume, propagator_state_t,
+                                    navigator_t>);
+
+  static_assert(isCurrentVolumePtr,
+                "Return type is not a known volume pointer type");
+
+  constexpr static bool value = Acts::Concepts::require<isCurrentVolumePtr>;
+};
+
+template <typename propagator_state_t, typename navigator_t>
+constexpr bool NavigationCompatibilityConcept =
+    NavigationCompatibilityConceptImpl<propagator_state_t, navigator_t>::value;
+
+}  // namespace Concepts
+
+}  // namespace Acts

Core/include/Acts/Navigation/DetectorNavigator.hpp

@@ -87,8 +87,8 @@ class DetectorNavigator {
     return state.currentSurface;
   }
 
-  const TrackingVolume* currentVolume(const State& /*state*/) const {
-    return nullptr;  // TODO we do not have a tracking volume
+  const DetectorVolume* currentVolume(const State& state) const {
+    return state.currentVolume;
   }
 
   const IVolumeMaterial* currentVolumeMaterial(const State& state) const {
@@ -146,14 +146,22 @@ class DetectorNavigator {
     auto& nState = state.navigation;
 
     if (nState.currentDetector == nullptr) {
-      ACTS_VERBOSE("assigning detector from the config.");
+      ACTS_VERBOSE("Assigning detector from the config.");
       nState.currentDetector = m_cfg.detector;
     }
-
     if (nState.currentDetector == nullptr) {
-      ACTS_ERROR("panic: no detector");
-      return;
+      throw std::invalid_argument("DetectorNavigator: no detector assigned");
+    }
+
+    fillNavigationState(state, stepper, nState);
+    if (nState.currentVolume == nullptr) {
+      nState.currentVolume = nState.currentDetector->findDetectorVolume(
+          state.geoContext, nState.position);
+    }
+    if (nState.currentVolume == nullptr) {
+      throw std::invalid_argument("DetectorNavigator: no current volume found");
     }
+    updateCandidateSurfaces(state, stepper);
   }
 
   /// @brief Navigator pre step call
@@ -171,18 +179,14 @@ class DetectorNavigator {
     ACTS_VERBOSE(volInfo(state)
                  << posInfo(state, stepper) << "Entering navigator::preStep.");
 
-    auto& nState = state.navigation;
-    fillNavigationState(state, stepper, nState);
-
     if (inactive()) {
       ACTS_VERBOSE(volInfo(state)
                    << posInfo(state, stepper) << "navigator inactive");
       return;
     }
 
-    if (nState.currentVolume == nullptr) {
-      initializeTarget(state, stepper);
-    }
+    auto& nState = state.navigation;
+    fillNavigationState(state, stepper, nState);
 
     if (nState.currentSurface != nullptr) {
       ACTS_VERBOSE(volInfo(state)
@@ -196,9 +200,18 @@ class DetectorNavigator {
 
       nState.currentPortal->updateDetectorVolume(state.geoContext, nState);
 
-      initializeTarget(state, stepper);
-    }
+      // If no Volume is found, we are at the end of the world
+      if (nState.currentVolume == nullptr) {
+        ACTS_VERBOSE(volInfo(state) << posInfo(state, stepper)
+                                    << "no volume after Portal update");
+        nState.navigationBreak = true;
+        return;
+      }
 
+      // Switched to a new volume
+      // Update candidate surfaces
+      updateCandidateSurfaces(state, stepper);
+    }
     for (; nState.surfaceCandidateIndex != nState.surfaceCandidates.size();
          ++nState.surfaceCandidateIndex) {
       // Screen output how much is left to try
@@ -213,19 +226,26 @@ class DetectorNavigator {
       const auto& surface =
           (c.surface != nullptr) ? (*c.surface) : (c.portal->surface());
       // Screen output which surface you are on
-      ACTS_VERBOSE(volInfo(state) << posInfo(state, stepper)
-                                  << "next surface candidate will be "
-                                  << surface.geometryId());
+      ACTS_VERBOSE(volInfo(state)
+                   << posInfo(state, stepper)
+                   << "next surface candidate will be " << surface.geometryId()
+                   << " (" << surface.center(state.geoContext).transpose()
+                   << ")");
       // Estimate the surface status
       bool boundaryCheck = c.boundaryCheck.isEnabled();
       auto surfaceStatus = stepper.updateSurfaceStatus(
           state.stepping, surface, c.objectIntersection.index(),
           state.options.direction, BoundaryCheck(boundaryCheck),
           state.options.surfaceTolerance, logger());
+
+      ACTS_VERBOSE(volInfo(state) << posInfo(state, stepper)
+                                  << "surface status is " << surfaceStatus);
+
       if (surfaceStatus == Intersection3D::Status::reachable) {
         ACTS_VERBOSE(volInfo(state)
-                     << posInfo(state, stepper)
-                     << "surface reachable, step size updated to "
+                     << posInfo(state, stepper) << "surface "
+                     << surface.center(state.geoContext).transpose()
+                     << " is reachable, step size updated to "
                      << stepper.outputStepSize(state.stepping));
         break;
       }
@@ -256,11 +276,6 @@ class DetectorNavigator {
       return;
     }
 
-    if (nState.currentDetector == nullptr) {
-      initialize(state, stepper);
-      return;
-    }
-
     if (nState.surfaceCandidateIndex == nState.surfaceCandidates.size()) {
       ACTS_VERBOSE(volInfo(state)
                    << posInfo(state, stepper)
@@ -280,10 +295,10 @@ class DetectorNavigator {
       nextSurface = &nextPortal->surface();
       isPortal = true;
     } else {
-      ACTS_ERROR(volInfo(state)
-                 << posInfo(state, stepper)
-                 << "panic: not a surface not a portal - what is it?");
-      return;
+      std::string msg = "DetectorNavigator: " + volInfo(state) +
+                        posInfo(state, stepper) +
+                        "panic: not a surface not a portal - what is it?";
+      throw std::runtime_error(msg);
     }
 
     // TODO not sure about the boundary check
@@ -381,31 +396,18 @@ class DetectorNavigator {
   ///
   /// boolean return triggers exit to stepper
   template <typename propagator_state_t, typename stepper_t>
-  void initializeTarget(propagator_state_t& state,
-                        const stepper_t& stepper) const {
+  void updateCandidateSurfaces(propagator_state_t& state,
+                               const stepper_t& stepper) const {
     ACTS_VERBOSE(volInfo(state)
                  << posInfo(state, stepper) << "initialize target");
 
     auto& nState = state.navigation;
 
-    if (nState.currentVolume == nullptr) {
-      nState.currentVolume = nState.currentDetector->findDetectorVolume(
-          state.geoContext, nState.position);
-
-      if (nState.currentVolume != nullptr) {
-        ACTS_VERBOSE(volInfo(state)
-                     << posInfo(state, stepper) << "switched detector volume");
-      }
-    }
-
-    if (nState.currentVolume == nullptr) {
-      ACTS_ERROR(volInfo(state)
-                 << posInfo(state, stepper) << "panic: no current volume");
-      return;
-    }
-
+    // Here we get the candidate surfaces
     nState.currentVolume->updateNavigationState(state.geoContext, nState);
 
+    ACTS_VERBOSE("SURFACE CANDIDATES: " << nState.surfaceCandidates.size());
+
     // Sort properly the surface candidates
     auto& nCandidates = nState.surfaceCandidates;
     std::sort(nCandidates.begin(), nCandidates.end(),
@@ -423,12 +425,15 @@ class DetectorNavigator {
   void fillNavigationState(propagator_state_t& state, const stepper_t& stepper,
                            NavigationState& nState) const {
     nState.position = stepper.position(state.stepping);
-    nState.direction = stepper.direction(state.stepping);
+    nState.direction =
+        state.options.direction * stepper.direction(state.stepping);
     nState.absMomentum = stepper.absoluteMomentum(state.stepping);
     auto fieldResult = stepper.getField(state.stepping, nState.position);
     if (!fieldResult.ok()) {
-      ACTS_ERROR(volInfo(state) << posInfo(state, stepper)
-                                << "could not read from the magnetic field");
+      std::string msg = "DetectorNavigator: " + volInfo(state) +
+                        posInfo(state, stepper) +
+                        "could not read from the magnetic field";
+      throw std::runtime_error(msg);
     }
     nState.magneticField = *fieldResult;
   }

Core/include/Acts/Navigation/NavigationStateUpdaters.hpp

@@ -44,19 +44,20 @@ inline void updateCandidates(const GeometryContext& gctx,
   NavigationState::SurfaceCandidates nextSurfaceCandidates;
 
   for (NavigationState::SurfaceCandidate c : nState.surfaceCandidates) {
-    // Get the surface representation: either native surfcae of portal
+    // Get the surface representation: either native surface of portal
     const Surface& sRep =
         c.surface != nullptr ? *c.surface : c.portal->surface();
 
     // Get the intersection @todo make a templated intersector
     // TODO surface tolerance
     auto sIntersection = sRep.intersect(gctx, position, direction,
                                         c.boundaryCheck, s_onSurfaceTolerance);
-    c.objectIntersection = sIntersection[c.objectIntersection.index()];
-
-    if (c.objectIntersection &&
-        c.objectIntersection.pathLength() > nState.overstepTolerance) {
-      nextSurfaceCandidates.emplace_back(std::move(c));
+    for (auto& si : sIntersection.split()) {
+      c.objectIntersection = si;
+      if (c.objectIntersection &&
+          c.objectIntersection.pathLength() > nState.overstepTolerance) {
+        nextSurfaceCandidates.emplace_back(c);
+      }
     }
   }
 

Core/include/Acts/Propagator/MaterialInteractor.hpp

@@ -9,6 +9,7 @@
 #pragma once
 
 #include "Acts/Definitions/Units.hpp"
+#include "Acts/Detector/GeometryCompatibilityConcept.hpp"
 #include "Acts/Geometry/TrackingVolume.hpp"
 #include "Acts/Material/MaterialInteraction.hpp"
 #include "Acts/Material/MaterialSlab.hpp"
@@ -67,6 +68,11 @@ struct MaterialInteractor {
       return;
     }
 
+    static_assert(
+        Acts::Concepts::NavigationCompatibilityConcept<propagator_state_t,
+                                                       navigator_t>,
+        "Navigation does not fulfill geometry compatibility concept");
+
     // Handle surface material
 
     // Note that start and target surface conditions are handled in the
@@ -119,7 +125,7 @@ struct MaterialInteractor {
       updateResult(state, stepper, result);
     }
 
-    const TrackingVolume* volume = navigator.currentVolume(state.navigation);
+    auto volume = navigator.currentVolume(state.navigation);
 
     // We only have material interactions if there is potential material
     if (volume && volume->volumeMaterial()) {

Core/include/Acts/Propagator/detail/VolumeMaterialInteraction.hpp

@@ -73,6 +73,30 @@ struct VolumeMaterialInteraction {
         performCovarianceTransport(state.stepping.covTransport),
         navDir(state.options.direction) {}
 
+  /// @brief Constructor
+  ///
+  /// @tparam propagator_state_t Type of the propagator state
+  /// @tparam stepper_t Type of the stepper
+  ///
+  /// @param [in] vVolume The current volume
+  /// @param [in] state State of the propagation
+  /// @param [in] stepper Stepper in use
+  template <typename propagator_state_t, typename stepper_t>
+  VolumeMaterialInteraction(const Acts::Experimental::DetectorVolume* vVolume,
+                            const propagator_state_t& state,
+                            const stepper_t& stepper)
+      : volume(vVolume),
+        pos(stepper.position(state.stepping)),
+        time(stepper.time(state.stepping)),
+        dir(stepper.direction(state.stepping)),
+        qOverP(stepper.qOverP(state.stepping)),
+        absQ(stepper.particleHypothesis(state.stepping).absoluteCharge()),
+        momentum(stepper.absoluteMomentum(state.stepping)),
+        mass(stepper.particleHypothesis(state.stepping).mass()),
+        absPdg(stepper.particleHypothesis(state.stepping).absolutePdg()),
+        performCovarianceTransport(state.stepping.covTransport),
+        navDir(state.options.direction) {}
+
   /// @brief This function evaluates the material properties to interact with
   ///
   /// @tparam propagator_state_t Type of the propagator state