helix_plane_intersector.hpp

// 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/.
 
#pragma once
 
// Project include(s)
#include "detray/definitions/algebra.hpp"
#include "detray/definitions/detail/qualifiers.hpp"
#include "detray/definitions/math.hpp"
#include "detray/definitions/units.hpp"
#include "detray/geometry/coordinates/cartesian2D.hpp"
#include "detray/geometry/coordinates/polar2D.hpp"
#include "detray/navigation/intersection/intersection.hpp"
#include "detray/tracks/helix.hpp"
#include "detray/utils/root_finding.hpp"
 
// System include(s)
#include <cmath>
#include <limits>
#include <type_traits>
 
namespace detray {
 
template <typename frame_t, concepts::algebra algebra_t>
struct helix_intersector_impl;

template <detray::concepts::aos algebra_t>
struct helix_intersector_impl<cartesian2D<algebra_t>, algebra_t> {
 private:
  using scalar_t = dscalar<algebra_t>;
 
 public:
  using algebra_type = algebra_t;

  template <typename surface_descr_t>
  using intersection_type =
      intersection2D<surface_descr_t, algebra_t, intersection::contains_pos>;

  template <typename other_algebra_t>
  using trajectory_type = detail::helix<other_algebra_t>;
 
  // Maximum number of solutions this intersector can produce
  static constexpr std::uint8_t n_solutions{1u};

  using result_type = intersection_point_err<algebra_t>;

  DETRAY_HOST_DEVICE constexpr result_type point_of_intersection(
      const trajectory_type<algebra_t> &h, const dtransform3D<algebra_t> &trf,
      const scalar_t /*unused*/ = 0.f) const {
    using point3_t = dpoint3D<algebra_t>;
    using vector3_t = dvector3D<algebra_t>;
 
    // Surface normal
    const vector3_t sn = trf.z();
    // Surface translation
    const point3_t st = trf.translation();
 
    // For a plane parallel to the magnetic field, the projected helix is
    // a circle. Reject planes outside that circle before attempting to find
    // a root: bracket expansion cannot succeed when no intersection exists.
    const vector3_t field = h.b_field();
    const scalar_t curvature = -h.qop() * h.B();
    if (vector::dot(sn, field) == 0.f && curvature != 0.f &&
        h.radius() * h.B() >= 1e-6f) {
      const vector3_t field_axis = vector::normalize(field);
      // The transverse circle centre is offset from the initial position.
      const vector3_t center_offset =
          (1.f / curvature) * vector::cross(field_axis, h.dir());
      const scalar_t center_distance = math::fabs(
          vector::dot(sn, h.pos() - st) + vector::dot(sn, center_offset));
      const scalar_t radius = h.radius();
      // Keep near-tangent cases within the solver's tolerance, with an
      // additional allowance for floating-point roundoff.
      const scalar_t margin = convergence_tolerance +
                              64.f * std::numeric_limits<scalar_t>::epsilon() *
                                  (center_distance + radius);
      if (std::isfinite(center_distance) && std::isfinite(radius) &&
          center_distance - radius > margin) {
        constexpr scalar_t inv = detail::invalid_value<scalar_t>();
        return {inv, point3_t{inv, inv, inv}, inv};
      }
    }
 
    // Starting point on the helix for the Newton iteration
    const vector3_t dist{st - h.pos(0.f)};
    scalar_t denom{vector::dot(sn, h.dir(0.5f * vector::norm(dist)))};
    scalar_t s_ini;
    if (denom == 0.f) {
      s_ini = vector::norm(dist);
    } else {
      s_ini = vector::dot(sn, dist) / denom;
    }

    auto plane_inters_func = [&h, &st, &sn](const scalar_t x) {
      // f(s) = sn * (h.pos(s) - st) == 0
      const scalar_t f_s{vector::dot(sn, (h.pos(x) - st))};
      // f'(s) = sn * h.dir(s)
      const scalar_t df_s{vector::dot(sn, h.dir(x))};
 
      return std::make_tuple(f_s, df_s);
    };
 
    // Run the root finding algorithm
    scalar_t s;
    scalar_t ds;
    if (run_rtsafe) {
      std::tie(s, ds) =
          newton_raphson_safe(plane_inters_func, s_ini, convergence_tolerance,
                              max_n_tries, max_path);
    } else {
      std::tie(s, ds) =
          newton_raphson(plane_inters_func, s_ini, convergence_tolerance,
                         max_n_tries, max_path);
    }
 
    return {s, h.pos(s), ds};
  }

  scalar_t convergence_tolerance{1.f * unit<scalar_t>::um};
  std::size_t max_n_tries{1000u};
  scalar_t max_path{5.f * unit<scalar_t>::m};
  bool run_rtsafe{true};
};
 
template <detray::concepts::aos algebra_t>
struct helix_intersector_impl<polar2D<algebra_t>, algebra_t>
    : public helix_intersector_impl<cartesian2D<algebra_t>, algebra_t> {};
 
}  // namespace detray