Skip to content

ref: Remove code duplication from helix intersectors #873

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
+469 −568
Open

ref: Remove code duplication from helix intersectors #873

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
309 changes: 106 additions & 203 deletions core/include/detray/navigation/intersection/helix_cylinder_intersector.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
/** Detray library, part of the ACTS project (R&D line)
*
* (c) 2022-2024 CERN for the benefit of the ACTS project
* (c) 2022-2025 CERN for the benefit of the ACTS project
*
* Mozilla Public License Version 2.0
*/
Expand All @@ -19,10 +19,7 @@
#include "detray/navigation/intersection/ray_cylinder_intersector.hpp"
#include "detray/tracks/helix.hpp"
#include "detray/utils/invalid_values.hpp"

// System include(s)
#include <limits>
#include <type_traits>
#include "detray/utils/root_finding.hpp"

namespace detray {

Expand All @@ -39,12 +36,11 @@ template <algebra::concepts::aos algebra_t>
struct helix_intersector_impl<cylindrical2D<algebra_t>, algebra_t>
: public ray_intersector_impl<cylindrical2D<algebra_t>, algebra_t,
intersection::contains_pos> {
private:
using scalar_t = dscalar<algebra_t>;

public:
using algebra_type = algebra_t;
using scalar_type = dscalar<algebra_t>;
using point3_type = dpoint3D<algebra_t>;
using vector3_type = dvector3D<algebra_t>;
using transform3_type = dtransform3D<algebra_t>;

template <typename surface_descr_t>
using intersection_type =
Expand All @@ -58,204 +54,117 @@ struct helix_intersector_impl<cylindrical2D<algebra_t>, algebra_t>

using result_type = darray<intersection_point_err<algebra_t>, n_solutions>;

/// Operator function to find intersections between ray and planar mask
/// Operator function to find intersections between a helix and a cylinder
/// surface
///
/// @param ray is the input ray trajectory
/// @param sf the surface handle the mask is associated with
/// @param mask is the input mask that defines the surface extent
/// @param h is the input helix trajectory
/// @param trf is the surface placement transform
/// @param mask_tolerance is the tolerance for mask edges
/// @param overstep_tol negative cutoff for the path
/// @param mask is the input mask that defines the surface extent
///
/// @return the intersection
template <typename mask_t>
DETRAY_HOST_DEVICE constexpr result_type point_of_intersection(
const trajectory_type<algebra_t> &h, const transform3_type &trf,
const mask_t &mask, const scalar_type = 0.f) const {
const trajectory_type<algebra_t> &h, const dtransform3D<algebra_t> &trf,
const mask_t &mask, const scalar_t = 0.f) const {

using point3_t = dpoint3D<algebra_t>;
using vector3_t = dvector3D<algebra_t>;

result_type ret{};

if (!run_rtsafe) {
// Get the surface placement
const auto &sm = trf.matrix();
// Cylinder z axis
const vector3_type sz = getter::vector<3>(sm, 0u, 2u);
// Cylinder centre
const point3_type sc = getter::vector<3>(sm, 0u, 3u);

// Starting point on the helix for the Newton iteration
// The mask is a cylinder -> it provides its radius as the first
// value
const scalar_type r{mask[cylinder2D::e_r]};

// Try to guess the best starting positions for the iteration

// Direction of the track at the helix origin
const auto h_dir = h.dir(0.f);
// Default starting path length for the Newton iteration (assumes
// concentric cylinder)
const scalar_type default_s{r * vector::perp(h_dir)};

// Initial helix path length parameter
darray<scalar_type, 2> paths{default_s, default_s};

// try to guess good starting path by calculating the intersection
// path of the helix tangential with the cylinder. This only has a
// chance of working for tracks with reasonably high p_T !
detail::ray<algebra_t> t{h.pos(), h.time(), h_dir, h.qop()};
const auto qe = this->solve_intersection(t, mask, trf);

// Obtain both possible solutions by looping over the (different)
// starting positions
auto n_runs{static_cast<unsigned int>(qe.solutions())};

// Note: the default path length might be smaller than either
// solution
switch (qe.solutions()) {
case 2:
paths[1] = qe.larger();
// If there are two solutions, reuse the case for a single
// solution to setup the intersection with the smaller path
// in ret[0]
[[fallthrough]];
case 1: {
paths[0] = qe.smaller();
break;
}
// Even if the ray is parallel to the cylinder, the helix
// might still hit it
default: {
n_runs = 2u;
paths[0] = r;
paths[1] = -r;
}
// Cylinder z axis
const vector3_t sz = trf.z();
// Cylinder centre
const point3_t sc = trf.translation();

// Starting point on the helix for the Newton iteration
// The mask is a cylinder -> it provides its radius as the first
// value
const scalar_t r{mask[cylinder2D::e_r]};

// Try to guess the best starting positions for the iteration

// Direction of the track at the helix origin
const auto h_dir = h.dir(0.5f * r);
// Default starting path length for the Newton iteration (assumes
// concentric cylinder)
const scalar_t default_s{r * vector::perp(h_dir)};

// Initial helix path length parameter
darray<scalar_t, 2> paths{default_s, default_s};

// try to guess good starting path by calculating the intersection
// path of the helix tangential with the cylinder. This only has a
// chance of working for tracks with reasonably high p_T !
detail::ray<algebra_t> t{h.pos(), h.time(), h_dir, h.qop()};
const auto qe = this->solve_intersection(t, mask, trf);

// Obtain both possible solutions by looping over the (different)
// starting positions
auto n_runs{static_cast<unsigned int>(qe.solutions())};

// Note: the default path length might be smaller than either
// solution
switch (qe.solutions()) {
case 2:
paths[1] = qe.larger();
// If there are two solutions, reuse the case for a single
// solution to setup the intersection with the smaller path
// in ret[0]
[[fallthrough]];
case 1: {
paths[0] = qe.smaller();
break;
}

for (unsigned int i = 0u; i < n_runs; ++i) {

scalar_type &s = paths[i];

// Path length in the previous iteration step
scalar_type s_prev{0.f};

// f(s) = ((h.pos(s) - sc) x sz)^2 - r^2 == 0
// Run the iteration on s
std::size_t n_tries{0u};
while (math::fabs(s - s_prev) > convergence_tolerance &&
n_tries < max_n_tries) {

// f'(s) = 2 * ( (h.pos(s) - sc) x sz) * (h.dir(s) x sz) )
const vector3_type crp = vector::cross(h.pos(s) - sc, sz);
const scalar_type denom{
2.f * vector::dot(crp, vector::cross(h.dir(s), sz))};

// No intersection can be found if dividing by zero
if (denom == 0.f) {
return {};
}

// x_n+1 = x_n - f(s) / f'(s)
s_prev = s;
s -= (vector::dot(crp, crp) - r * r) / denom;

++n_tries;
}
// No intersection found within max number of trials
if (n_tries == max_n_tries) {
return {};
}

ret[i] = {s, h.pos(s), s - s_prev};
default: {
n_runs = 2u;
paths[0] = r;
paths[1] = -r;
}
}

return ret;
} else {
// Cylinder z axis
const vector3_type sz = trf.z();
// Cylinder centre
const point3_type sc = trf.translation();

// Starting point on the helix for the Newton iteration
// The mask is a cylinder -> it provides its radius as the first
// value
const scalar_type r{mask[cylinder2D::e_r]};

// Try to guess the best starting positions for the iteration

// Direction of the track at the helix origin
const auto h_dir = h.dir(0.5f * r);
// Default starting path length for the Newton iteration (assumes
// concentric cylinder)
const scalar_type default_s{r * vector::perp(h_dir)};

// Initial helix path length parameter
darray<scalar_type, 2> paths{default_s, default_s};

// try to guess good starting path by calculating the intersection
// path of the helix tangential with the cylinder. This only has a
// chance of working for tracks with reasonably high p_T !
detail::ray<algebra_t> t{h.pos(), h.time(), h_dir, h.qop()};
const auto qe = this->solve_intersection(t, mask, trf);

// Obtain both possible solutions by looping over the (different)
// starting positions
auto n_runs{static_cast<unsigned int>(qe.solutions())};

// Note: the default path length might be smaller than either
// solution
switch (qe.solutions()) {
case 2:
paths[1] = qe.larger();
// If there are two solutions, reuse the case for a single
// solution to setup the intersection with the smaller path
// in ret[0]
[[fallthrough]];
case 1: {
paths[0] = qe.smaller();
break;
}
default: {
n_runs = 2u;
paths[0] = r;
paths[1] = -r;
}
}

/// Evaluate the function and its derivative at the point @param x
auto cyl_inters_func = [&h, &r, &sz, &sc](const scalar_type x) {
const vector3_type crp = vector::cross(h.pos(x) - sc, sz);

// f(s) = ((h.pos(s) - sc) x sz)^2 - r^2 == 0
const scalar_type f_s{(vector::dot(crp, crp) - r * r)};
// f'(s) = 2 * ( (h.pos(s) - sc) x sz) * (h.dir(s) x sz) )
const scalar_type df_s{
2.f * vector::dot(crp, vector::cross(h.dir(x), sz))};

return std::make_tuple(f_s, df_s);
};

for (unsigned int i = 0u; i < n_runs; ++i) {

const scalar_type &s_ini = paths[i];

// Run the root finding algorithm
const auto [s, ds] = newton_raphson_safe(cyl_inters_func, s_ini,
convergence_tolerance,
max_n_tries, max_path);

ret[i] = {s, h.pos(s), ds};
/// Evaluate the function and its derivative at the point @param x
auto cyl_inters_func = [&h, &r, &sz, &sc](const scalar_t x) {
const vector3_t crp = vector::cross(h.pos(x) - sc, sz);

// f(s) = ((h.pos(s) - sc) x sz)^2 - r^2 == 0
const scalar_t f_s{(vector::dot(crp, crp) - r * r)};
// f'(s) = 2 * ( (h.pos(s) - sc) x sz) * (h.dir(s) x sz) )
const scalar_t df_s{2.f *
vector::dot(crp, vector::cross(h.dir(x), sz))};

return std::make_tuple(f_s, df_s);
};

for (unsigned int i = 0u; i < n_runs; ++i) {

const scalar_t &s_ini = paths[i];

// Run the root finding algorithm
scalar_t s;
scalar_t ds;
if (run_rtsafe) {
std::tie(s, ds) = newton_raphson_safe(cyl_inters_func, s_ini,
convergence_tolerance,
max_n_tries, max_path);
} else {
std::tie(s, ds) = newton_raphson(cyl_inters_func, s_ini,
convergence_tolerance,
max_n_tries, max_path);
}

return ret;
ret[i] = {s, h.pos(s), ds};
}

return ret;
}

/// Tolerance for convergence
scalar_type convergence_tolerance{1.f * unit<scalar_type>::um};
scalar_t convergence_tolerance{1.f * unit<scalar_t>::um};
// Guard against inifinite loops
std::size_t max_n_tries{1000u};
// Early exit, if the intersection is too far away
scalar_type max_path{5.f * unit<scalar_type>::m};
scalar_t max_path{5.f * unit<scalar_t>::m};
// Complement the Newton algorithm with Bisection steps
bool run_rtsafe{true};
};
Expand All @@ -264,14 +173,7 @@ template <concepts::algebra algebra_t>
struct helix_intersector_impl<concentric_cylindrical2D<algebra_t>, algebra_t>
: public helix_intersector_impl<cylindrical2D<algebra_t>, algebra_t> {

using base_type =
helix_intersector_impl<cylindrical2D<algebra_t>, algebra_t>;

using algebra_type = algebra_t;
using scalar_type = dscalar<algebra_t>;
using point3_type = dpoint3D<algebra_t>;
using vector3_type = dvector3D<algebra_t>;
using transform3_type = dtransform3D<algebra_t>;

template <typename surface_descr_t>
using intersection_type =
Expand All @@ -285,24 +187,25 @@ struct helix_intersector_impl<concentric_cylindrical2D<algebra_t>, algebra_t>

using result_type = intersection_point_err<algebra_t>;

/// Operator function to find intersections between ray and planar mask
/// Operator function to find intersections between helix and a
/// concentric cylinder surface
///
/// @param ray is the input ray trajectory
/// @param sf the surface handle the mask is associated with
/// @param mask is the input mask that defines the surface extent
/// @param h is the input helix trajectory
/// @param trf is the surface placement transform
/// @param mask_tolerance is the tolerance for mask edges
/// @param overstep_tol negative cutoff for the path
/// @param mask is the input mask that defines the surface extent
///
/// @return the intersection
template <typename mask_t>
DETRAY_HOST_DEVICE constexpr result_type point_of_intersection(
const trajectory_type<algebra_t> &h, const transform3_type &trf,
const mask_t &mask, const scalar_type = 0.f) const {
const trajectory_type<algebra_t> &h, const dtransform3D<algebra_t> &trf,
const mask_t &mask, const dscalar<algebra_t> = 0.f) const {

using base_t =
helix_intersector_impl<cylindrical2D<algebra_t>, algebra_t>;

// Array of two solutions
typename base_type::result_type results =
base_type::point_of_intersection(h, trf, mask, 0.f);
typename base_t::result_type results =
base_t::point_of_intersection(h, trf, mask, 0.f);

// For portals, only take the solution along the helix direction,
// not the one in the opposite direction
Expand Down
Loading
Loading