Acts/Navigation/DetectorNavigator.hpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2023 CERN for the benefit of the Acts project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 #include "Acts/Definitions/Units.hpp"
0012 #include "Acts/Detector/Detector.hpp"
0013 #include "Acts/Detector/DetectorVolume.hpp"
0014 #include "Acts/Detector/Portal.hpp"
0015 #include "Acts/Geometry/BoundarySurfaceT.hpp"
0016 #include "Acts/Geometry/GeometryIdentifier.hpp"
0017 #include "Acts/Geometry/Layer.hpp"
0018 #include "Acts/Navigation/NavigationState.hpp"
0019 #include "Acts/Propagator/NavigatorOptions.hpp"
0020 #include "Acts/Propagator/Propagator.hpp"
0021 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0022 #include "Acts/Surfaces/Surface.hpp"
0023 #include "Acts/Utilities/Logger.hpp"
0024
0025 #include <iomanip>
0026 #include <iterator>
0027 #include <sstream>
0028 #include <string>
0029
0030 #include <boost/algorithm/string.hpp>
0031 #include <boost/container/small_vector.hpp>
0032
0033 namespace Acts::Experimental {
0034
0035 class DetectorNavigator {
0036  public:
0037   struct Config {
0038     /// Detector for this Navigation
0039     const Detector* detector = nullptr;
0040
0041     /// Configuration for this Navigator
0042     /// stop at every sensitive surface (whether it has material or not)
0043     bool resolveSensitive = true;
0044     /// stop at every material surface (whether it is passive or not)
0045     bool resolveMaterial = true;
0046     /// stop at every surface regardless what it is
0047     bool resolvePassive = false;
0048   };
0049
0050   struct Options : public NavigatorPlainOptions {
0051     void setPlainOptions(const NavigatorPlainOptions& options) {
0052       static_cast<NavigatorPlainOptions&>(*this) = options;
0053     }
0054   };
0055
0056   /// Nested State struct
0057   ///
0058   /// It acts as an internal state which is
0059   /// created for every propagation/extrapolation step
0060   /// and keep thread-local navigation information
0061   struct State : public NavigationState {
0062     Options options;
0063
0064     /// Navigation state - external state: the current surface
0065     const Surface* currentSurface = nullptr;
0066     /// Indicator if the target is reached
0067     bool targetReached = false;
0068     /// Navigation state : a break has been detected
0069     bool navigationBreak = false;
0070   };
0071
0072   /// Constructor with configuration object
0073   ///
0074   /// @param cfg The navigator configuration
0075   /// @param _logger a logger instance
0076   explicit DetectorNavigator(Config cfg,
0077                              std::shared_ptr<const Logger> _logger =
0078                                  getDefaultLogger("DetectorNavigator",
0079                                                   Logging::Level::INFO))
0080       : m_cfg{cfg}, m_logger{std::move(_logger)} {}
0081
0082   State makeState(const Options& options) const {
0083     State state;
0084     state.options = options;
0085     return state;
0086   }
0087
0088   const Surface* currentSurface(const State& state) const {
0089     return state.currentSurface;
0090   }
0091
0092   const DetectorVolume* currentVolume(const State& state) const {
0093     return state.currentVolume;
0094   }
0095
0096   const IVolumeMaterial* currentVolumeMaterial(const State& state) const {
0097     return state.currentVolume->volumeMaterial();
0098   }
0099
0100   const Surface* startSurface(const State& state) const {
0101     return state.options.startSurface;
0102   }
0103
0104   const Surface* targetSurface(const State& state) const {
0105     return state.options.targetSurface;
0106   }
0107
0108   bool targetReached(const State& state) const { return state.targetReached; }
0109
0110   bool endOfWorldReached(State& state) const {
0111     return state.currentVolume == nullptr;
0112   }
0113
0114   bool navigationBreak(const State& state) const {
0115     return state.navigationBreak;
0116   }
0117
0118   void currentSurface(State& state, const Surface* surface) const {
0119     state.currentSurface = surface;
0120   }
0121
0122   void targetReached(State& state, bool targetReached) const {
0123     state.targetReached = targetReached;
0124   }
0125
0126   void navigationBreak(State& state, bool navigationBreak) const {
0127     state.navigationBreak = navigationBreak;
0128   }
0129
0130   /// Initialize call - start of propagation
0131   ///
0132   /// @tparam propagator_state_t The state type of the propagator
0133   /// @tparam stepper_t The type of stepper used for the propagation
0134   ///
0135   /// @param [in,out] state is the propagation state object
0136   /// @param [in] stepper Stepper in use
0137   ///
0138   /// @return boolean return triggers exit to stepper
0139   template <typename propagator_state_t, typename stepper_t>
0140   void initialize(propagator_state_t& state, const stepper_t& stepper) const {
0141     ACTS_VERBOSE(volInfo(state) << posInfo(state, stepper) << "initialize");
0142
0143     auto& nState = state.navigation;
0144
0145     if (nState.currentDetector == nullptr) {
0146       ACTS_VERBOSE("Assigning detector from the config.");
0147       nState.currentDetector = m_cfg.detector;
0148     }
0149     if (nState.currentDetector == nullptr) {
0150       throw std::invalid_argument("DetectorNavigator: no detector assigned");
0151     }
0152
0153     fillNavigationState(state, stepper, nState);
0154     if (nState.currentVolume == nullptr) {
0155       nState.currentVolume = nState.currentDetector->findDetectorVolume(
0156           state.geoContext, nState.position);
0157     }
0158     if (nState.currentVolume == nullptr) {
0159       throw std::invalid_argument("DetectorNavigator: no current volume found");
0160     }
0161     updateCandidateSurfaces(state, stepper);
0162   }
0163
0164   /// @brief Navigator pre step call
0165   ///
0166   /// This will invalid the current surface and current portal in order
0167   /// to navigate to the next ones.
0168   ///
0169   /// @tparam propagator_state_t is the type of Propagatgor state
0170   /// @tparam stepper_t is the used type of the Stepper by the Propagator
0171   ///
0172   /// @param [in,out] state is the mutable propagator state object
0173   /// @param [in] stepper Stepper in use
0174   template <typename propagator_state_t, typename stepper_t>
0175   void preStep(propagator_state_t& state, const stepper_t& stepper) const {
0176     ACTS_VERBOSE(volInfo(state)
0177                  << posInfo(state, stepper) << "Entering navigator::preStep.");
0178
0179     if (inactive()) {
0180       ACTS_VERBOSE(volInfo(state)
0181                    << posInfo(state, stepper) << "navigator inactive");
0182       return;
0183     }
0184
0185     auto& nState = state.navigation;
0186     fillNavigationState(state, stepper, nState);
0187
0188     if (nState.currentSurface != nullptr) {
0189       ACTS_VERBOSE(volInfo(state)
0190                    << posInfo(state, stepper) << "stepping through surface");
0191     }
0192
0193     for (; nState.surfaceCandidateIndex != nState.surfaceCandidates.size();
0194          ++nState.surfaceCandidateIndex) {
0195       // Screen output how much is left to try
0196       ACTS_VERBOSE(volInfo(state)
0197                    << posInfo(state, stepper)
0198                    << (nState.surfaceCandidates.size() -
0199                        nState.surfaceCandidateIndex)
0200                    << " out of " << nState.surfaceCandidates.size()
0201                    << " surfaces remain to try.");
0202       // Take the surface
0203       const auto& c = nState.surfaceCandidate();
0204       const auto& surface =
0205           (c.surface != nullptr) ? (*c.surface) : (c.portal->surface());
0206       // Screen output which surface you are on
0207       ACTS_VERBOSE(volInfo(state)
0208                    << posInfo(state, stepper)
0209                    << "next surface candidate will be " << surface.geometryId()
0210                    << " (" << surface.center(state.geoContext).transpose()
0211                    << ")");
0212       // Estimate the surface status
0213       auto surfaceStatus = stepper.updateSurfaceStatus(
0214           state.stepping, surface, c.objectIntersection.index(),
0215           state.options.direction, c.boundaryTolerance,
0216           state.options.surfaceTolerance, logger());
0217
0218       ACTS_VERBOSE(volInfo(state) << posInfo(state, stepper)
0219                                   << "surface status is " << surfaceStatus);
0220
0221       if (surfaceStatus == Intersection3D::Status::reachable) {
0222         ACTS_VERBOSE(volInfo(state)
0223                      << posInfo(state, stepper) << "surface "
0224                      << surface.center(state.geoContext).transpose()
0225                      << " is reachable, step size updated to "
0226                      << stepper.outputStepSize(state.stepping));
0227         break;
0228       }
0229     }
0230
0231     nState.currentSurface = nullptr;
0232     nState.currentPortal = nullptr;
0233   }
0234
0235   /// @brief Navigator post step call
0236   ///
0237   /// @tparam propagator_state_t is the type of Propagatgor state
0238   /// @tparam stepper_t is the used type of the Stepper by the Propagator
0239   ///
0240   /// @param [in,out] state is the mutable propagator state object
0241   /// @param [in] stepper Stepper in use
0242   template <typename propagator_state_t, typename stepper_t>
0243   void postStep(propagator_state_t& state, const stepper_t& stepper) const {
0244     ACTS_VERBOSE(volInfo(state)
0245                  << posInfo(state, stepper) << "Entering navigator::postStep.");
0246
0247     auto& nState = state.navigation;
0248     fillNavigationState(state, stepper, nState);
0249
0250     if (inactive()) {
0251       ACTS_VERBOSE(volInfo(state)
0252                    << posInfo(state, stepper) << "navigator inactive");
0253       return;
0254     }
0255
0256     if (nState.surfaceCandidateIndex == nState.surfaceCandidates.size()) {
0257       ACTS_VERBOSE(volInfo(state)
0258                    << posInfo(state, stepper)
0259                    << "no surface candidates - waiting for target call");
0260       return;
0261     }
0262
0263     const Portal* nextPortal = nullptr;
0264     const Surface* nextSurface = nullptr;
0265     bool isPortal = false;
0266     BoundaryTolerance boundaryTolerance =
0267         nState.surfaceCandidate().boundaryTolerance;
0268
0269     if (nState.surfaceCandidate().surface != nullptr) {
0270       nextSurface = nState.surfaceCandidate().surface;
0271     } else if (nState.surfaceCandidate().portal != nullptr) {
0272       nextPortal = nState.surfaceCandidate().portal;
0273       nextSurface = &nextPortal->surface();
0274       isPortal = true;
0275     } else {
0276       std::string msg = "DetectorNavigator: " + volInfo(state) +
0277                         posInfo(state, stepper) +
0278                         "panic: not a surface not a portal - what is it?";
0279       throw std::runtime_error(msg);
0280     }
0281
0282     // TODO not sure about the boundary check
0283     auto surfaceStatus = stepper.updateSurfaceStatus(
0284         state.stepping, *nextSurface,
0285         nState.surfaceCandidate().objectIntersection.index(),
0286         state.options.direction, boundaryTolerance,
0287         state.options.surfaceTolerance, logger());
0288
0289     // Check if we are at a surface
0290     if (surfaceStatus == Intersection3D::Status::onSurface) {
0291       ACTS_VERBOSE(volInfo(state)
0292                    << posInfo(state, stepper) << "landed on surface");
0293
0294       if (isPortal) {
0295         ACTS_VERBOSE(volInfo(state)
0296                      << posInfo(state, stepper)
0297                      << "this is a portal, updating to new volume.");
0298         nState.currentPortal = nextPortal;
0299         nState.currentSurface = &nextPortal->surface();
0300         nState.surfaceCandidates.clear();
0301         nState.surfaceCandidateIndex = 0;
0302
0303         nState.currentPortal->updateDetectorVolume(state.geoContext, nState);
0304
0305         // If no Volume is found, we are at the end of the world
0306         if (nState.currentVolume == nullptr) {
0307           ACTS_VERBOSE(volInfo(state)
0308                        << posInfo(state, stepper)
0309                        << "no volume after Portal update, end of world.");
0310           nState.navigationBreak = true;
0311           return;
0312         }
0313
0314         // Switched to a new volume
0315         // Update candidate surfaces
0316         updateCandidateSurfaces(state, stepper);
0317
0318         ACTS_VERBOSE(volInfo(state)
0319                      << posInfo(state, stepper) << "current portal set to "
0320                      << nState.currentPortal->surface().geometryId());
0321
0322       } else {
0323         ACTS_VERBOSE(volInfo(state) << posInfo(state, stepper)
0324                                     << "this is a surface, storing it.");
0325
0326         // If we are on the surface pointed at by the iterator, we can make
0327         // it the current one to pass it to the other actors
0328         nState.currentSurface = nextSurface;
0329         ACTS_VERBOSE(volInfo(state)
0330                      << posInfo(state, stepper) << "current surface set to "
0331                      << nState.currentSurface->geometryId());
0332         ++nState.surfaceCandidateIndex;
0333       }
0334     }
0335   }
0336
0337  private:
0338   Config m_cfg;
0339
0340   std::shared_ptr<const Logger> m_logger;
0341
0342   template <typename propagator_state_t>
0343   std::string volInfo(const propagator_state_t& state) const {
0344     auto& nState = state.navigation;
0345
0346     return (nState.currentVolume ? nState.currentVolume->name() : "No Volume") +
0347            " | ";
0348   }
0349
0350   template <typename propagator_state_t, typename stepper_t>
0351   std::string posInfo(const propagator_state_t& state,
0352                       const stepper_t& stepper) const {
0353     std::stringstream ss;
0354     ss << stepper.position(state.stepping).transpose();
0355     ss << " | ";
0356     return ss.str();
0357   }
0358
0359   const Logger& logger() const { return *m_logger; }
0360
0361   /// This checks if a navigation break had been triggered or navigator
0362   /// is misconfigured
0363   ///
0364   /// boolean return triggers exit to stepper
0365   bool inactive() const {
0366     if (m_cfg.detector == nullptr) {
0367       return true;
0368     }
0369
0370     if (!m_cfg.resolveSensitive && !m_cfg.resolveMaterial &&
0371         !m_cfg.resolvePassive) {
0372       return true;
0373     }
0374
0375     return false;
0376   }
0377
0378   /// @brief Navigation (re-)initialisation for the target
0379   ///
0380   /// @note This is only called a few times every propagation/extrapolation
0381   ///
0382   /// As a straight line estimate can lead you to the wrong destination
0383   /// Volume, this will be called at:
0384   /// - initialization
0385   /// - attempted volume switch
0386   /// Target finding by association will not be done again
0387   ///
0388   /// @tparam propagator_state_t The state type of the propagator
0389   /// @tparam stepper_t The type of stepper used for the propagation
0390   ///
0391   /// @param [in,out] state is the propagation state object
0392   /// @param [in] stepper Stepper in use
0393   ///
0394   /// boolean return triggers exit to stepper
0395   template <typename propagator_state_t, typename stepper_t>
0396   void updateCandidateSurfaces(propagator_state_t& state,
0397                                const stepper_t& stepper) const {
0398     ACTS_VERBOSE(volInfo(state)
0399                  << posInfo(state, stepper) << "initialize target");
0400
0401     auto& nState = state.navigation;
0402
0403     // Here we get the candidate surfaces
0404     nState.currentVolume->updateNavigationState(state.geoContext, nState);
0405
0406     // Sort properly the surface candidates
0407     auto& nCandidates = nState.surfaceCandidates;
0408     std::sort(nCandidates.begin(), nCandidates.end(),
0409               [&](const auto& a, const auto& b) {
0410                 // The two path lengths
0411                 ActsScalar pathToA = a.objectIntersection.pathLength();
0412                 ActsScalar pathToB = b.objectIntersection.pathLength();
0413                 return pathToA < pathToB;
0414               });
0415     // Set the surface candidate
0416     nState.surfaceCandidateIndex = 0;
0417   }
0418
0419   template <typename propagator_state_t, typename stepper_t>
0420   void fillNavigationState(propagator_state_t& state, const stepper_t& stepper,
0421                            NavigationState& nState) const {
0422     nState.position = stepper.position(state.stepping);
0423     nState.direction =
0424         state.options.direction * stepper.direction(state.stepping);
0425     nState.absMomentum = stepper.absoluteMomentum(state.stepping);
0426     auto fieldResult = stepper.getField(state.stepping, nState.position);
0427     if (!fieldResult.ok()) {
0428       std::string msg = "DetectorNavigator: " + volInfo(state) +
0429                         posInfo(state, stepper) +
0430                         "could not read from the magnetic field";
0431       throw std::runtime_error(msg);
0432     }
0433     nState.magneticField = *fieldResult;
0434   }
0435 };
0436
0437 }  // namespace Acts::Experimental