feat(detray): perigee extrapolator

Detray/codegen/detray-sympy/requirements.txt

@@ -0,0 +1,10 @@
+# This file was autogenerated by uv via the following command:
+#    uv pip compile --universal --python-version 3.11 Detray/codegen/detray-sympy/pyproject.toml -o Detray/codegen/detray-sympy/requirements.txt
+mpmath==1.3.0
+    # via sympy
+numpy==2.4.6 ; python_full_version < '3.12'
+    # via detray-sympy (Detray/codegen/detray-sympy/pyproject.toml)
+numpy==2.5.3 ; python_full_version >= '3.12'
+    # via detray-sympy (Detray/codegen/detray-sympy/pyproject.toml)
+sympy==1.14.0
+    # via detray-sympy (Detray/codegen/detray-sympy/pyproject.toml)

Detray/core/include/detray/navigation/direct_navigator.hpp

@@ -30,7 +30,7 @@
 
 namespace detray {
 
-template <concepts::detector detector_t, typename surface_t = void>
+template <typename detector_t, typename surface_sequence_t = void>
 class direct_navigator {
   using algebra_t = typename detector_t::algebra_type;
   using scalar_t = dscalar<algebra_t>;
@@ -39,10 +39,18 @@ class direct_navigator {
   using detector_type = detector_t;
   using context_type = detector_type::geometry_context;
 
-  // Use an external surface type, if one was provided
-  using surface_type =
-      std::conditional_t<std::same_as<surface_t, void>,
-                         typename detector_t::surface_type, surface_t>;
+  // Use an external surface sequence type (rather than a single surface type),
+  // perigee/external surface flows instantiate direct_navigator with
+  // single_view-like range wrappers.
+  using surface_sequence_type = std::conditional_t<
+      std::same_as<surface_sequence_t, void>,
+      vecmem::device_vector<typename detector_t::surface_type>,
+      surface_sequence_t>;
+
+  static_assert(detray::ranges::range<surface_sequence_type>,
+                "Surface sequence must be iterable");
+
+  using surface_type = detray::ranges::range_value_t<surface_sequence_type>;
   using intersection_type =
       intersection2D<surface_type, algebra_t, !intersection::contains_pos>;
   using inspector_type = navigation::void_inspector;
@@ -64,10 +72,16 @@ class direct_navigator {
 
    public:
     using value_type = intersection_type;
-    using sequence_type = vecmem::device_vector<surface_type>;
+    using sequence_type = surface_sequence_type;
 
-    using view_type = dvector_view<surface_type>;
-    using const_view_type = dvector_view<const surface_type>;
+    using view_type =
+        std::conditional_t<std::same_as<surface_sequence_t, void>,
+                           detail::get_view_t<surface_sequence_type>,
+                           surface_sequence_type>;
+    using const_view_type =
+        std::conditional_t<std::same_as<surface_sequence_t, void>,
+                           detail::get_view_t<const surface_sequence_type>,
+                           surface_sequence_type>;
 
     /// Constructor using the detector and an externally provided sequence
     /// of detector surfaces
@@ -247,9 +261,7 @@ class direct_navigator {
 
       // Update the current target. If it cannot be reached, direct
       // navigation is broken
-      if (!navigation::update_candidate(navigation.target(), tangential, det,
-                                        intr_cfg, navigation.external_tol(),
-                                        ctx)) {
+      if (!update_candidate(tangential, det, navigation, intr_cfg, ctx)) {
         navigation.abort("Could not reach current target");
         return !is_init;
       }
@@ -273,9 +285,7 @@ class direct_navigator {
 
       // Otherwise, track is on surface: Update the next target
       if (navigation.has_next_external() &&
-          !navigation::update_candidate(navigation.target(), tangential, det,
-                                        intr_cfg, navigation.external_tol(),
-                                        ctx)) {
+          !update_candidate(tangential, det, navigation, intr_cfg, ctx)) {
         navigation.abort("Could not find new target after surface was reached");
         return !is_init;
       }
@@ -301,6 +311,59 @@ class direct_navigator {
 
     return !is_init;
   }
+
+ private:
+  /// Update the next candidate: Either an externally specified detector
+  /// surface or a standalone surface.
+  DETRAY_HOST_DEVICE DETRAY_INLINE constexpr bool update_candidate(
+      const detray::detail::ray<algebra_t> &tangential, const detector_t &det,
+      state &navigation, const intersection::config &intr_cfg,
+      [[maybe_unused]] const context_type &ctx) const {
+    if constexpr (std::same_as<surface_sequence_t, void>) {
+      // Update a detector surface.
+      return navigation::update_candidate(navigation.target(), tangential, det,
+                                          intr_cfg, navigation.external_tol(),
+                                          ctx);
+    } else {
+      // Update a standalone surface.
+      return update_external_candidate(
+          navigation.direction(), navigation.target(), tangential,
+          navigation.next_external(), intr_cfg, navigation.external_tol());
+    }
+  }
+
+  /// Update a standalone surface.
+  DETRAY_HOST_DEVICE DETRAY_INLINE constexpr bool update_external_candidate(
+      const navigation::direction nav_dir,
+      typename state::value_type &candidate,
+      const detray::detail::ray<algebra_t> &tangential, const surface_type &sf,
+      const intersection::config &cfg,
+      const scalar_t external_mask_tolerance) const {
+    constexpr ray_intersector<typename surface_type::mask_type::shape,
+                              algebra_t, !intersection::contains_pos>
+        intersector{};
+
+    // Intersect the standalone surface.
+    typename decltype(intersector)::result_type result{};
+    if constexpr (concepts::cylindrical<typename surface_type::mask_type>) {
+      result = intersector.point_of_intersection(
+          tangential, sf.transform(), sf.mask(), cfg.overstep_tolerance);
+    } else {
+      result = intersector.point_of_intersection(tangential, sf.transform(),
+                                                 cfg.overstep_tolerance);
+    }
+
+    // Build resulting intersection candidate(s).
+    if constexpr (decltype(intersector)::n_solutions > 1) {
+      resolve_mask(candidate, tangential, result[0], candidate.surface(),
+                   sf.mask(), sf.transform(), cfg, external_mask_tolerance);
+    } else {
+      resolve_mask(candidate, tangential, result, candidate.surface(),
+                   sf.mask(), sf.transform(), cfg, external_mask_tolerance);
+    }
+
+    return candidate.is_probably_inside();
+  }
 };
 
 }  // namespace detray

Detray/core/include/detray/navigation/external_surface.hpp

@@ -0,0 +1,194 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/.
+
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2025 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Project include(s)
+#include "detray/definitions/detail/qualifiers.hpp"
+#include "detray/definitions/geometry.hpp"
+#include "detray/definitions/indexing.hpp"
+#include "detray/definitions/units.hpp"
+#include "detray/geometry/identifier.hpp"
+#include "detray/geometry/mask.hpp"
+
+namespace detray {
+
+template <typename algebra_t, typename shape_t, typename material_t,
+          typename nav_link_t>
+class external_surface {
+ public:
+  using transform3_type = dtransform3D<algebra_t>;
+  using mask_type = detray::mask<shape_t, algebra_t, nav_link_t>;
+  using material_type = material_t;
+  /// Link type of the mask to a volume.
+  using navigation_link = nav_link_t;
+
+  /// Default constructor
+  constexpr external_surface() = default;
+
+  /// Constructor with full arguments
+  ///
+  /// @param trf the transform for positioning and 3D local frame
+  /// @param mask the surface mask
+  /// @param material the surface material
+  /// @param volume the volume this surface belongs to
+  /// @param sf_id remember whether this is a portal or not
+  DETRAY_HOST_DEVICE
+  constexpr external_surface(const transform3_type& trf, const mask_type& mask,
+                             const material_type& material, const dindex volume,
+                             const surface_id sf_id)
+      : m_transform{trf},
+        m_mask(mask),
+        m_material(material),
+        m_barcode{geometry::identifier{}.set_volume(volume).set_id(sf_id)} {}
+
+  /// Constructor with full arguments - move semantics
+  ///
+  /// @param trf the transform for positioning and 3D local frame
+  /// @param mask the surface mask
+  /// @param material the surface material
+  /// @param volume the volume this surface belongs to
+  /// @param sf_id remember whether this is a portal or not
+  DETRAY_HOST_DEVICE
+  constexpr external_surface(transform3_type&& trf, mask_type&& mask,
+                             material_type&& material, const dindex volume,
+                             const surface_id sf_id)
+      : m_transform{std::move(trf)},
+        m_mask(std::move(mask)),
+        m_material(std::move(material)),
+        m_barcode{geometry::identifier{}.set_volume(volume).set_id(sf_id)} {}
+
+  /// Equality operator
+  ///
+  /// @param rhs is the right hand side to be compared to
+  DETRAY_HOST_DEVICE
+  constexpr auto operator==(const external_surface& rhs) const -> bool {
+    return (m_barcode == rhs.m_barcode);
+  }
+
+  /// Sets a new surface barcode
+  DETRAY_HOST_DEVICE
+  auto set_barcode(const geometry::identifier bcd) -> void { m_barcode = bcd; }
+
+  /// @returns the surface barcode
+  DETRAY_HOST_DEVICE
+  constexpr auto barcode() const -> geometry::identifier { return m_barcode; }
+
+  /// @returns the surface identifier (identifier-era naming)
+  DETRAY_HOST_DEVICE
+  constexpr auto identifier() const -> geometry::identifier {
+    return m_barcode;
+  }
+
+  /// Implicit conversion for APIs that take geometry::identifier.
+  DETRAY_HOST_DEVICE
+  constexpr operator geometry::identifier() const { return m_barcode; }
+
+  /// Sets a new surface id (portal/passive/sensitive)
+  DETRAY_HOST_DEVICE
+  auto set_id(const surface_id new_id) -> void { m_barcode.set_id(new_id); }
+
+  /// @returns the surface id (sensitive, passive or portal)
+  DETRAY_HOST_DEVICE
+  constexpr auto id() const -> surface_id { return m_barcode.id(); }
+
+  /// Sets a new volume link (index in volume collection of detector)
+  DETRAY_HOST_DEVICE
+  auto set_volume(const dindex new_idx) -> void {
+    m_barcode.set_volume(new_idx);
+  }
+
+  /// @returns the surface id (sensitive, passive or portal)
+  DETRAY_HOST_DEVICE
+  constexpr auto volume() const -> dindex { return m_barcode.volume(); }
+
+  /// Sets a new surface index (index in surface collection of surface store)
+  DETRAY_HOST_DEVICE
+  auto set_index(const dindex new_idx) -> void { m_barcode.set_index(new_idx); }
+
+  /// @returns the surface id (sensitive, passive or portal)
+  DETRAY_HOST_DEVICE
+  constexpr auto index() const -> dindex { return m_barcode.index(); }
+
+  /// Update the transform index
+  ///
+  /// @param offset update the position when move into new collection
+  DETRAY_HOST_DEVICE
+  void set_transform(const transform3_type& trf) { m_transform = trf; }
+
+  /// @return the transform index
+  DETRAY_HOST_DEVICE
+  constexpr auto transform() const -> const transform3_type& {
+    return m_transform;
+  }
+
+  /// Update the mask link
+  ///
+  /// @param offset update the position when move into new collection
+  DETRAY_HOST_DEVICE
+  void set_mask(const mask_type m) { m_mask = m; }
+
+  /// @return the mask - const
+  DETRAY_HOST_DEVICE
+  constexpr auto mask() const -> const mask_type& { return m_mask; }
+
+  /// Update the material link
+  ///
+  /// @param offset update the position when move into new collection
+  DETRAY_HOST_DEVICE
+  void set_material(const material_type& offset) { m_material.shift(offset); }
+
+  /// Access to the material
+  DETRAY_HOST_DEVICE
+  constexpr auto material() -> material_type& { return m_material; }
+
+  /// @return the material link
+  DETRAY_HOST_DEVICE
+  constexpr auto material() const -> const material_type& { return m_material; }
+
+  /// @returns true if the surface descriptor has a valid material link
+  DETRAY_HOST_DEVICE
+  constexpr auto has_material() const -> bool { return false; }
+
+  /// @returns true if the surface is a sensitive detector module.
+  DETRAY_HOST_DEVICE
+  constexpr auto is_sensitive() const -> bool {
+    return m_barcode.id() == surface_id::e_sensitive;
+  }
+
+  /// @returns true if the surface is a portal.
+  DETRAY_HOST_DEVICE
+  constexpr auto is_portal() const -> bool {
+    return m_barcode.id() == surface_id::e_portal;
+  }
+
+  /// @returns true if the surface is a passive detector element.
+  DETRAY_HOST_DEVICE
+  constexpr auto is_passive() const -> bool {
+    return m_barcode.id() == surface_id::e_passive;
+  }
+
+ private:
+  /// Surface placement transform and its inverse
+  transform3_type m_transform{};
+  /// Surface mask
+  mask_type m_mask{};
+  /// Surface material
+  material_type m_material{};
+  /// Surface hash
+  geometry::identifier m_barcode{};
+};
+
+}  // namespace detray