diff --git a/Core/include/Acts/Propagator/TryAllNavigator.hpp b/Core/include/Acts/Propagator/TryAllNavigator.hpp index 6f9b5df3e12..7f4a2157923 100644 --- a/Core/include/Acts/Propagator/TryAllNavigator.hpp +++ b/Core/include/Acts/Propagator/TryAllNavigator.hpp @@ -27,16 +27,8 @@ namespace Acts { -/// @brief Alternative @c Navigator which tries all possible intersections -/// -/// See @c Navigator for more general information about the Navigator concept. -/// -/// This Navigator tries all possible intersections with all surfaces in the -/// current volume. It does not use any information about the geometry to -/// optimize the search. It is therefore very slow, but can be used as a -/// reference implementation. -/// -class TryAllNavigator { +/// @brief Captures the common functionality of the `TryAllNavigator`s +class TryAllNavigatorBase { public: /// @brief Configuration for this Navigator struct Config { @@ -63,7 +55,6 @@ class TryAllNavigator { // while initialization. NOTE: This information is mostly used by actors to // check if we are on the starting surface (e.g. MaterialInteraction). const Surface* startSurface = nullptr; - const TrackingVolume* startVolume = nullptr; // Target geometry information of the navigation which should only be set // while initialization. NOTE: This information is mostly used by actors to @@ -90,19 +81,9 @@ class TryAllNavigator { /// /// @param cfg The navigator configuration /// @param _logger a logger instance - TryAllNavigator(Config cfg, - std::unique_ptr _logger = - getDefaultLogger("TryAllNavigator", Logging::INFO)) + TryAllNavigatorBase(Config cfg, std::unique_ptr _logger) : m_cfg(std::move(cfg)), m_logger{std::move(_logger)} {} - State makeState(const Surface* startSurface, - const Surface* targetSurface) const { - State result; - result.startSurface = startSurface; - result.targetSurface = targetSurface; - return result; - } - const Surface* currentSurface(const State& state) const { return state.currentSurface; } @@ -195,6 +176,100 @@ class TryAllNavigator { ACTS_VERBOSE(volInfo(state) << "No start surface set."); } } + } + + protected: + /// Helper method to initialize navigation candidates for the current volume. + template + void initializeVolumeCandidates(propagator_state_t& state) const { + const TrackingVolume* volume = state.navigation.currentVolume; + ACTS_VERBOSE(volInfo(state) << "Initialize volume"); + + if (volume == nullptr) { + state.navigation.navigationBreak = true; + ACTS_VERBOSE(volInfo(state) << "No volume set. Good luck."); + return; + } + + emplaceAllVolumeCandidates(state.navigation.navigationCandidates, *volume, + m_cfg.resolveSensitive, m_cfg.resolveMaterial, + m_cfg.resolvePassive, + m_cfg.boundaryCheckSurfaceApproach, logger()); + } + + template + std::string volInfo(const propagator_state_t& state) const { + return (state.navigation.currentVolume != nullptr + ? state.navigation.currentVolume->volumeName() + : "No Volume") + + " | "; + } + + const Logger& logger() const { return *m_logger; } + + Config m_cfg; + + std::unique_ptr m_logger; +}; + +/// @brief Alternative @c Navigator which tries all possible intersections +/// +/// See @c Navigator for more general information about the Navigator concept. +/// +/// This Navigator tries all possible intersections with all surfaces in the +/// current volume. It does not use any information about the geometry to +/// optimize the search. It is therefore very slow, but can be used as a +/// reference implementation. +/// +class TryAllNavigator : public TryAllNavigatorBase { + public: + /// @brief Configuration for this Navigator + struct Config : public TryAllNavigatorBase::Config {}; + + /// @brief Nested State struct + /// + /// It acts as an internal state which is created for every propagation and + /// meant to keep thread-local navigation information. + struct State : public TryAllNavigatorBase::State {}; + + /// Constructor with configuration object + /// + /// @param cfg The navigator configuration + /// @param logger a logger instance + TryAllNavigator(Config cfg, + std::unique_ptr logger = + getDefaultLogger("TryAllNavigator", Logging::INFO)) + : TryAllNavigatorBase(std::move(cfg), std::move(logger)) {} + + State makeState(const Surface* startSurface, + const Surface* targetSurface) const { + State result; + result.startSurface = startSurface; + result.targetSurface = targetSurface; + return result; + } + + using TryAllNavigatorBase::currentSurface; + using TryAllNavigatorBase::currentVolume; + using TryAllNavigatorBase::currentVolumeMaterial; + using TryAllNavigatorBase::endOfWorldReached; + using TryAllNavigatorBase::navigationBreak; + using TryAllNavigatorBase::startSurface; + using TryAllNavigatorBase::targetReached; + using TryAllNavigatorBase::targetSurface; + + using TryAllNavigatorBase::initialize; + + /// @brief Initialize call - start of navigation + /// + /// @tparam propagator_state_t The state type of the propagator + /// @tparam stepper_t The type of stepper used for the propagation + /// + /// @param [in,out] state is the propagation state object + /// @param [in] stepper Stepper in use + template + void initialize(propagator_state_t& state, const stepper_t& stepper) const { + TryAllNavigatorBase::initialize(state, stepper); // Initialize navigation candidates for the start volume reinitializeCandidates(state); @@ -430,46 +505,404 @@ class TryAllNavigator { } private: - /// Helper method to initialize navigation candidates for the current volume. + /// Helper method to reset and reinitialize the navigation candidates. template - void initializeVolumeCandidates(propagator_state_t& state) const { - const TrackingVolume* volume = state.navigation.currentVolume; - ACTS_VERBOSE(volInfo(state) << "Initialize volume"); + void reinitializeCandidates(propagator_state_t& state) const { + state.navigation.navigationCandidates.clear(); + state.navigation.intersectionCandidates.clear(); + + initializeVolumeCandidates(state); + } +}; + +/// @brief Alternative @c Navigator which tries all possible intersections +/// +/// See @c Navigator for more general information about the Navigator concept. +/// +/// This Navigator tries all possible intersections with all surfaces in the +/// current volume. It does not use any information about the geometry to +/// optimize the search. It is therefore very slow, but can be used as a +/// reference implementation. +/// +/// Different to @c TryAllNavigator, this Navigator discovers intersections by +/// stepping forward blindly and then checking for intersections with all +/// surfaces in the current volume. This is slower, but more robust against +/// bent tracks. +/// +class TryAllOverstepNavigator : public TryAllNavigatorBase { + public: + /// @brief Configuration for this Navigator + struct Config : public TryAllNavigatorBase::Config {}; + + /// @brief Nested State struct + /// + /// It acts as an internal state which is created for every propagation and + /// meant to keep thread-local navigation information. + struct State : public TryAllNavigatorBase::State { + /// The vector of navigation candidates to work through + std::vector navigationCandidates; + /// The vector of active intersection candidates to work through + std::vector activeCandidates; + /// The current active candidate index of the navigation state + std::size_t activeCandidateIndex = 0; + + /// The position before the last step + std::optional lastPosition; + /// The last intersection used to avoid rehitting the same surface + std::optional lastIntersection; + + /// Provides easy access to the active intersection candidate + const detail::IntersectionCandidate& activeCandidate() const { + return activeCandidates.at(activeCandidateIndex); + } + }; + + /// Constructor with configuration object + /// + /// @param cfg The navigator configuration + /// @param logger a logger instance + TryAllOverstepNavigator(Config cfg, + std::unique_ptr logger = + getDefaultLogger("TryAllOverstepNavigator", + Logging::INFO)) + : TryAllNavigatorBase(std::move(cfg), std::move(logger)) {} + + State makeState(const Surface* startSurface, + const Surface* targetSurface) const { + State result; + result.startSurface = startSurface; + result.targetSurface = targetSurface; + return result; + } + + const Surface* currentSurface(const State& state) const { + return state.currentSurface; + } + + const TrackingVolume* currentVolume(const State& state) const { + return state.currentVolume; + } + + const IVolumeMaterial* currentVolumeMaterial(const State& state) const { + if (state.currentVolume == nullptr) { + return nullptr; + } + return state.currentVolume->volumeMaterial(); + } + + const Surface* startSurface(const State& state) const { + return state.startSurface; + } + + const Surface* targetSurface(const State& state) const { + return state.targetSurface; + } + + bool targetReached(const State& state) const { return state.targetReached; } + + bool endOfWorldReached(State& state) const { + return state.currentVolume == nullptr; + } + + bool navigationBreak(const State& state) const { + return state.navigationBreak; + } + + void currentSurface(State& state, const Surface* surface) const { + state.currentSurface = surface; + } + + void targetReached(State& state, bool targetReached) const { + state.targetReached = targetReached; + } + + void navigationBreak(State& state, bool navigationBreak) const { + state.navigationBreak = navigationBreak; + } + + /// @brief Initialize call - start of navigation + /// + /// @tparam propagator_state_t The state type of the propagator + /// @tparam stepper_t The type of stepper used for the propagation + /// + /// @param [in,out] state is the propagation state object + /// @param [in] stepper Stepper in use + template + void initialize(propagator_state_t& state, const stepper_t& stepper) const { + TryAllNavigatorBase::initialize(state, stepper); + + // Initialize navigation candidates for the start volume + reinitializeCandidates(state); + + state.navigation.lastPosition.reset(); + state.navigation.lastIntersection.reset(); + } + + /// @brief Navigator pre step call + /// + /// This determines the next surface to be targeted and sets the step length + /// accordingly. + /// + /// @tparam propagator_state_t is the type of Propagatgor state + /// @tparam stepper_t is the used type of the Stepper by the Propagator + /// + /// @param [in,out] state is the mutable propagator state object + /// @param [in] stepper Stepper in use + template + void preStep(propagator_state_t& state, const stepper_t& stepper) const { + ACTS_VERBOSE(volInfo(state) << "pre step"); + + // Navigator preStep always resets the current surface + state.navigation.currentSurface = nullptr; + + ACTS_VERBOSE(volInfo(state) << "handle active candidates"); + + // Check next navigation candidate + while (state.navigation.activeCandidateIndex != + state.navigation.activeCandidates.size()) { + // Screen output how much is left to try + ACTS_VERBOSE(volInfo(state) + << (state.navigation.activeCandidates.size() - + state.navigation.activeCandidateIndex) + << " out of " << state.navigation.activeCandidates.size() + << " surfaces remain to try."); + + const auto& candidate = state.navigation.activeCandidate(); + const auto& intersection = candidate.intersection; + const Surface& surface = *intersection.object(); + BoundaryCheck boundaryCheck = candidate.boundaryCheck; + + // Screen output which surface you are on + ACTS_VERBOSE(volInfo(state) << "Next surface candidate will be " + << surface.geometryId()); + + // Estimate the surface status + auto surfaceStatus = stepper.updateSurfaceStatus( + state.stepping, surface, intersection.index(), + state.options.direction, boundaryCheck, + state.options.surfaceTolerance, logger()); + + if (surfaceStatus == IntersectionStatus::onSurface) { + ACTS_ERROR(volInfo(state) + << "We are on surface " << surface.geometryId() + << " before trying to reach it. This should not happen. " + "Good luck."); + ++state.navigation.activeCandidateIndex; + continue; + } + + if (surfaceStatus == IntersectionStatus::reachable) { + ACTS_VERBOSE(volInfo(state) + << "Surface reachable, step size updated to " + << stepper.outputStepSize(state.stepping)); + break; + } + + ACTS_WARNING(volInfo(state) << "Surface unreachable, skip."); + ++state.navigation.activeCandidateIndex; + } + + if (state.navigation.activeCandidateIndex == + state.navigation.activeCandidates.size()) { + state.navigation.lastPosition = stepper.position(state.stepping); + + stepper.releaseStepSize(state.stepping, ConstrainedStep::actor); + + ACTS_VERBOSE(volInfo(state) + << "blindly step forwards. step size updated to " + << stepper.outputStepSize(state.stepping)); - if (volume == nullptr) { - state.navigation.navigationBreak = true; - ACTS_VERBOSE(volInfo(state) << "No volume set. Good luck."); return; } + } - emplaceAllVolumeCandidates(state.navigation.navigationCandidates, *volume, - m_cfg.resolveSensitive, m_cfg.resolveMaterial, - m_cfg.resolvePassive, - m_cfg.boundaryCheckSurfaceApproach, logger()); + /// @brief Navigator post step call + /// + /// This determines if we hit the next navigation candidate and deals with it + /// accordingly. It sets the current surface, enters layers and changes + /// volumes. + /// + /// @tparam propagator_state_t is the type of Propagatgor state + /// @tparam stepper_t is the used type of the Stepper by the Propagator + /// + /// @param [in,out] state is the mutable propagator state object + /// @param [in] stepper Stepper in use + /// + /// @return Boolean to indicate if we continue with the actors and + /// aborters or if we should target again. + template + void postStep(propagator_state_t& state, const stepper_t& stepper) const { + ACTS_VERBOSE(volInfo(state) << "post step"); + + assert(state.navigation.currentSurface == nullptr && + "Current surface must be reset."); + + if (state.navigation.activeCandidateIndex == + state.navigation.activeCandidates.size()) { + ACTS_VERBOSE(volInfo(state) << "evaluate blind step"); + + state.navigation.activeCandidates.clear(); + + assert(state.navigation.lastPosition.has_value() && + "last position not set"); + + Vector3 start = state.navigation.lastPosition.value(); + Vector3 end = stepper.position(state.stepping); + Vector3 step = end - start; + double distance = step.norm(); + Vector3 direction = step.normalized(); + + double nearLimit = -distance + state.options.surfaceTolerance; + double farLimit = state.options.surfaceTolerance; + + // Find intersections with all candidates + for (const auto& candidate : state.navigation.navigationCandidates) { + auto intersections = candidate.intersect( + state.geoContext, end, direction, state.options.surfaceTolerance); + for (const auto& intersection : intersections.first.split()) { + // exclude invalid intersections + if (!intersection || + !detail::checkIntersection(intersection, nearLimit, farLimit)) { + continue; + } + // exclude last candidate + if (state.navigation.lastIntersection.has_value() && + state.navigation.lastIntersection->object() == + intersection.object() && + state.navigation.lastIntersection->index() == + intersection.index()) { + continue; + } + // store candidate + state.navigation.activeCandidates.emplace_back( + intersection, intersections.second, candidate.boundaryCheck); + } + } + + std::sort(state.navigation.activeCandidates.begin(), + state.navigation.activeCandidates.end(), + detail::IntersectionCandidate::forwardOrder); + + state.navigation.activeCandidateIndex = 0; + + ACTS_VERBOSE(volInfo(state) + << "Found " << state.navigation.activeCandidates.size() + << " intersections"); + } + + if (state.navigation.activeCandidateIndex != + state.navigation.activeCandidates.size()) { + ACTS_VERBOSE(volInfo(state) << "handle active candidates"); + + std::vector hitCandidates; + + while (state.navigation.activeCandidateIndex != + state.navigation.activeCandidates.size()) { + const auto& candidate = state.navigation.activeCandidate(); + const auto& intersection = candidate.intersection; + const Surface& surface = *intersection.object(); + + Intersection3D::Status surfaceStatus = stepper.updateSurfaceStatus( + state.stepping, surface, intersection.index(), + state.options.direction, BoundaryCheck(false), + state.options.surfaceTolerance, logger()); + + if (surfaceStatus != IntersectionStatus::onSurface) { + break; + } + + hitCandidates.emplace_back(candidate); + + ++state.navigation.activeCandidateIndex; + } + + if (hitCandidates.empty()) { + ACTS_VERBOSE(volInfo(state) << "Staying focussed on surface."); + return; + } + + state.navigation.lastIntersection.reset(); + + std::vector trueHitCandidates; + + for (const auto& candidate : hitCandidates) { + const auto& intersection = candidate.intersection; + const Surface& surface = *intersection.object(); + + Intersection3D::Status surfaceStatus = stepper.updateSurfaceStatus( + state.stepping, surface, intersection.index(), + state.options.direction, BoundaryCheck(true), + state.options.surfaceTolerance, logger()); + + if (surfaceStatus != IntersectionStatus::onSurface) { + continue; + } + + trueHitCandidates.emplace_back(candidate); + } + + ACTS_VERBOSE(volInfo(state) + << "Found " << trueHitCandidates.size() + << " intersections on surface with bounds check."); + + if (trueHitCandidates.empty()) { + ACTS_VERBOSE(volInfo(state) + << "Surface successfully hit, but outside bounds."); + return; + } + + if (trueHitCandidates.size() > 1) { + ACTS_VERBOSE(volInfo(state) + << "Only using first intersection within bounds."); + } + + const auto& candidate = trueHitCandidates.front(); + const auto& intersection = candidate.intersection; + const Surface& surface = *intersection.object(); + + state.navigation.lastIntersection = intersection; + + ACTS_VERBOSE(volInfo(state) << "Surface successfully hit, storing it."); + // Set in navigation state, so actors and aborters can access it + state.navigation.currentSurface = &surface; + if (state.navigation.currentSurface) { + ACTS_VERBOSE(volInfo(state) << "Current surface set to surface " + << surface.geometryId()); + } + + if (candidate.template checkType()) { + ACTS_VERBOSE(volInfo(state) << "This is a surface"); + } else if (candidate.template checkType()) { + ACTS_VERBOSE(volInfo(state) << "This is a layer"); + } else if (candidate.template checkType()) { + ACTS_VERBOSE(volInfo(state) + << "This is a boundary. Reinitialize navigation"); + + const auto& boundary = *candidate.template object(); + + state.navigation.currentVolume = boundary.attachedVolume( + state.geoContext, stepper.position(state.stepping), + state.options.direction * stepper.direction(state.stepping)); + + ACTS_VERBOSE(volInfo(state) << "Switched volume"); + + reinitializeCandidates(state); + } else { + ACTS_ERROR(volInfo(state) << "Unknown intersection type"); + } + } } + private: /// Helper method to reset and reinitialize the navigation candidates. template void reinitializeCandidates(propagator_state_t& state) const { state.navigation.navigationCandidates.clear(); - state.navigation.intersectionCandidates.clear(); + state.navigation.activeCandidates.clear(); + state.navigation.activeCandidateIndex = 0; initializeVolumeCandidates(state); } - - template - std::string volInfo(const propagator_state_t& state) const { - return (state.navigation.currentVolume != nullptr - ? state.navigation.currentVolume->volumeName() - : "No Volume") + - " | "; - } - - const Logger& logger() const { return *m_logger; } - - Config m_cfg; - - std::unique_ptr m_logger; }; } // namespace Acts diff --git a/Tests/UnitTests/Core/Propagator/NavigatorTests.cpp b/Tests/UnitTests/Core/Propagator/NavigatorTests.cpp index 7591c1548d0..cfda6d847fb 100644 --- a/Tests/UnitTests/Core/Propagator/NavigatorTests.cpp +++ b/Tests/UnitTests/Core/Propagator/NavigatorTests.cpp @@ -1001,6 +1001,10 @@ using StraightLinePropagator = Propagator; using Reference1EigenPropagator = Propagator; using Reference1StraightLinePropagator = Propagator; +using Reference2EigenPropagator = + Propagator; +using Reference2StraightLinePropagator = + Propagator; EigenStepper estepper(bField); StraightLineStepper slstepper; @@ -1026,6 +1030,18 @@ Reference1StraightLinePropagator refslpropagator1( getDefaultLogger("ref1_sl_nav", Logging::INFO)), getDefaultLogger("ref1_sl_prop", Logging::INFO)); +Reference2EigenPropagator refepropagator2( + estepper, + TryAllOverstepNavigator({tGeometry, true, true, false, + BoundaryCheck(false)}, + getDefaultLogger("ref2_e_nav", Logging::INFO)), + getDefaultLogger("ref2_e_prop", Logging::INFO)); +Reference2StraightLinePropagator refslpropagator2( + slstepper, + TryAllOverstepNavigator({tGeometry, true, true, false}, + getDefaultLogger("ref2_sl_nav", Logging::INFO)), + getDefaultLogger("ref2_sl_prop", Logging::INFO)); + BOOST_DATA_TEST_CASE( Navigator_random, bdata::random((bdata::engine = std::mt19937(), bdata::seed = 20, @@ -1070,10 +1086,22 @@ BOOST_DATA_TEST_CASE( runSelfConsistencyTest(slpropagator, start, debugMode); if (debugMode) { - std::cout << ">>> Test consistency epropagator vs refepropagator1" - << std::endl; + std::cout << ">>> Test reference 1 consistency epropagator" << std::endl; } runConsistencyTest(epropagator, refepropagator1, start, debugMode); + if (debugMode) { + std::cout << ">>> Test reference 1 consistency slpropagator" << std::endl; + } + runConsistencyTest(slpropagator, refslpropagator1, start, debugMode); + + if (debugMode) { + std::cout << ">>> Test reference 2 consistency epropagator" << std::endl; + } + runConsistencyTest(epropagator, refepropagator2, start, debugMode); + if (debugMode) { + std::cout << ">>> Test reference 2 consistency slpropagator" << std::endl; + } + runConsistencyTest(slpropagator, refslpropagator2, start, debugMode); } } // namespace Acts::Test