Merge branch 'main' into save-vertex-seed

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