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