Skip to content

Commit aabd5d3

Browse files
niermann999niermann999
committed
Deduplaucate the newton methods in the helix intersectors
1 parent d7c0aab commit aabd5d3

File tree

8 files changed

+469
-568
lines changed

8 files changed

+469
-568
lines changed

‎core/include/detray/navigation/intersection/helix_cylinder_intersector.hpp‎

Lines changed: 106 additions & 203 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
/** Detray library, part of the ACTS project (R&D line)
22
*
3-
* (c) 2022-2024 CERN for the benefit of the ACTS project
3+
* (c) 2022-2025 CERN for the benefit of the ACTS project
44
*
55
* Mozilla Public License Version 2.0
66
*/
@@ -19,10 +19,7 @@
1919
#include "detray/navigation/intersection/ray_cylinder_intersector.hpp"
2020
#include "detray/tracks/helix.hpp"
2121
#include "detray/utils/invalid_values.hpp"
22-
23-
// System include(s)
24-
#include <limits>
25-
#include <type_traits>
22+
#include "detray/utils/root_finding.hpp"
2623

2724
namespace detray {
2825

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

42+
public:
4343
using algebra_type = algebra_t;
44-
using scalar_type = dscalar<algebra_t>;
45-
using point3_type = dpoint3D<algebra_t>;
46-
using vector3_type = dvector3D<algebra_t>;
47-
using transform3_type = dtransform3D<algebra_t>;
4844

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

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

61-
/// Operator function to find intersections between ray and planar mask
57+
/// Operator function to find intersections between a helix and a cylinder
58+
/// surface
6259
///
63-
/// @param ray is the input ray trajectory
64-
/// @param sf the surface handle the mask is associated with
65-
/// @param mask is the input mask that defines the surface extent
60+
/// @param h is the input helix trajectory
6661
/// @param trf is the surface placement transform
67-
/// @param mask_tolerance is the tolerance for mask edges
68-
/// @param overstep_tol negative cutoff for the path
62+
/// @param mask is the input mask that defines the surface extent
6963
///
7064
/// @return the intersection
7165
template <typename mask_t>
7266
DETRAY_HOST_DEVICE constexpr result_type point_of_intersection(
73-
const trajectory_type<algebra_t> &h, const transform3_type &trf,
74-
const mask_t &mask, const scalar_type = 0.f) const {
67+
const trajectory_type<algebra_t> &h, const dtransform3D<algebra_t> &trf,
68+
const mask_t &mask, const scalar_t = 0.f) const {
69+
70+
using point3_t = dpoint3D<algebra_t>;
71+
using vector3_t = dvector3D<algebra_t>;
7572

7673
result_type ret{};
7774

78-
if (!run_rtsafe) {
79-
// Get the surface placement
80-
const auto &sm = trf.matrix();
81-
// Cylinder z axis
82-
const vector3_type sz = getter::vector<3>(sm, 0u, 2u);
83-
// Cylinder centre
84-
const point3_type sc = getter::vector<3>(sm, 0u, 3u);
85-
86-
// Starting point on the helix for the Newton iteration
87-
// The mask is a cylinder -> it provides its radius as the first
88-
// value
89-
const scalar_type r{mask[cylinder2D::e_r]};
90-
91-
// Try to guess the best starting positions for the iteration
92-
93-
// Direction of the track at the helix origin
94-
const auto h_dir = h.dir(0.f);
95-
// Default starting path length for the Newton iteration (assumes
96-
// concentric cylinder)
97-
const scalar_type default_s{r * vector::perp(h_dir)};
98-
99-
// Initial helix path length parameter
100-
darray<scalar_type, 2> paths{default_s, default_s};
101-
102-
// try to guess good starting path by calculating the intersection
103-
// path of the helix tangential with the cylinder. This only has a
104-
// chance of working for tracks with reasonably high p_T !
105-
detail::ray<algebra_t> t{h.pos(), h.time(), h_dir, h.qop()};
106-
const auto qe = this->solve_intersection(t, mask, trf);
107-
108-
// Obtain both possible solutions by looping over the (different)
109-
// starting positions
110-
auto n_runs{static_cast<unsigned int>(qe.solutions())};
111-
112-
// Note: the default path length might be smaller than either
113-
// solution
114-
switch (qe.solutions()) {
115-
case 2:
116-
paths[1] = qe.larger();
117-
// If there are two solutions, reuse the case for a single
118-
// solution to setup the intersection with the smaller path
119-
// in ret[0]
120-
[[fallthrough]];
121-
case 1: {
122-
paths[0] = qe.smaller();
123-
break;
124-
}
125-
// Even if the ray is parallel to the cylinder, the helix
126-
// might still hit it
127-
default: {
128-
n_runs = 2u;
129-
paths[0] = r;
130-
paths[1] = -r;
131-
}
75+
// Cylinder z axis
76+
const vector3_t sz = trf.z();
77+
// Cylinder centre
78+
const point3_t sc = trf.translation();
79+
80+
// Starting point on the helix for the Newton iteration
81+
// The mask is a cylinder -> it provides its radius as the first
82+
// value
83+
const scalar_t r{mask[cylinder2D::e_r]};
84+
85+
// Try to guess the best starting positions for the iteration
86+
87+
// Direction of the track at the helix origin
88+
const auto h_dir = h.dir(0.5f * r);
89+
// Default starting path length for the Newton iteration (assumes
90+
// concentric cylinder)
91+
const scalar_t default_s{r * vector::perp(h_dir)};
92+
93+
// Initial helix path length parameter
94+
darray<scalar_t, 2> paths{default_s, default_s};
95+
96+
// try to guess good starting path by calculating the intersection
97+
// path of the helix tangential with the cylinder. This only has a
98+
// chance of working for tracks with reasonably high p_T !
99+
detail::ray<algebra_t> t{h.pos(), h.time(), h_dir, h.qop()};
100+
const auto qe = this->solve_intersection(t, mask, trf);
101+
102+
// Obtain both possible solutions by looping over the (different)
103+
// starting positions
104+
auto n_runs{static_cast<unsigned int>(qe.solutions())};
105+
106+
// Note: the default path length might be smaller than either
107+
// solution
108+
switch (qe.solutions()) {
109+
case 2:
110+
paths[1] = qe.larger();
111+
// If there are two solutions, reuse the case for a single
112+
// solution to setup the intersection with the smaller path
113+
// in ret[0]
114+
[[fallthrough]];
115+
case 1: {
116+
paths[0] = qe.smaller();
117+
break;
132118
}
133-
134-
for (unsigned int i = 0u; i < n_runs; ++i) {
135-
136-
scalar_type &s = paths[i];
137-
138-
// Path length in the previous iteration step
139-
scalar_type s_prev{0.f};
140-
141-
// f(s) = ((h.pos(s) - sc) x sz)^2 - r^2 == 0
142-
// Run the iteration on s
143-
std::size_t n_tries{0u};
144-
while (math::fabs(s - s_prev) > convergence_tolerance &&
145-
n_tries < max_n_tries) {
146-
147-
// f'(s) = 2 * ( (h.pos(s) - sc) x sz) * (h.dir(s) x sz) )
148-
const vector3_type crp = vector::cross(h.pos(s) - sc, sz);
149-
const scalar_type denom{
150-
2.f * vector::dot(crp, vector::cross(h.dir(s), sz))};
151-
152-
// No intersection can be found if dividing by zero
153-
if (denom == 0.f) {
154-
return {};
155-
}
156-
157-
// x_n+1 = x_n - f(s) / f'(s)
158-
s_prev = s;
159-
s -= (vector::dot(crp, crp) - r * r) / denom;
160-
161-
++n_tries;
162-
}
163-
// No intersection found within max number of trials
164-
if (n_tries == max_n_tries) {
165-
return {};
166-
}
167-
168-
ret[i] = {s, h.pos(s), s - s_prev};
119+
default: {
120+
n_runs = 2u;
121+
paths[0] = r;
122+
paths[1] = -r;
169123
}
124+
}
170125

171-
return ret;
172-
} else {
173-
// Cylinder z axis
174-
const vector3_type sz = trf.z();
175-
// Cylinder centre
176-
const point3_type sc = trf.translation();
177-
178-
// Starting point on the helix for the Newton iteration
179-
// The mask is a cylinder -> it provides its radius as the first
180-
// value
181-
const scalar_type r{mask[cylinder2D::e_r]};
182-
183-
// Try to guess the best starting positions for the iteration
184-
185-
// Direction of the track at the helix origin
186-
const auto h_dir = h.dir(0.5f * r);
187-
// Default starting path length for the Newton iteration (assumes
188-
// concentric cylinder)
189-
const scalar_type default_s{r * vector::perp(h_dir)};
190-
191-
// Initial helix path length parameter
192-
darray<scalar_type, 2> paths{default_s, default_s};
193-
194-
// try to guess good starting path by calculating the intersection
195-
// path of the helix tangential with the cylinder. This only has a
196-
// chance of working for tracks with reasonably high p_T !
197-
detail::ray<algebra_t> t{h.pos(), h.time(), h_dir, h.qop()};
198-
const auto qe = this->solve_intersection(t, mask, trf);
199-
200-
// Obtain both possible solutions by looping over the (different)
201-
// starting positions
202-
auto n_runs{static_cast<unsigned int>(qe.solutions())};
203-
204-
// Note: the default path length might be smaller than either
205-
// solution
206-
switch (qe.solutions()) {
207-
case 2:
208-
paths[1] = qe.larger();
209-
// If there are two solutions, reuse the case for a single
210-
// solution to setup the intersection with the smaller path
211-
// in ret[0]
212-
[[fallthrough]];
213-
case 1: {
214-
paths[0] = qe.smaller();
215-
break;
216-
}
217-
default: {
218-
n_runs = 2u;
219-
paths[0] = r;
220-
paths[1] = -r;
221-
}
222-
}
223-
224-
/// Evaluate the function and its derivative at the point @param x
225-
auto cyl_inters_func = [&h, &r, &sz, &sc](const scalar_type x) {
226-
const vector3_type crp = vector::cross(h.pos(x) - sc, sz);
227-
228-
// f(s) = ((h.pos(s) - sc) x sz)^2 - r^2 == 0
229-
const scalar_type f_s{(vector::dot(crp, crp) - r * r)};
230-
// f'(s) = 2 * ( (h.pos(s) - sc) x sz) * (h.dir(s) x sz) )
231-
const scalar_type df_s{
232-
2.f * vector::dot(crp, vector::cross(h.dir(x), sz))};
233-
234-
return std::make_tuple(f_s, df_s);
235-
};
236-
237-
for (unsigned int i = 0u; i < n_runs; ++i) {
238-
239-
const scalar_type &s_ini = paths[i];
240-
241-
// Run the root finding algorithm
242-
const auto [s, ds] = newton_raphson_safe(cyl_inters_func, s_ini,
243-
convergence_tolerance,
244-
max_n_tries, max_path);
245-
246-
ret[i] = {s, h.pos(s), ds};
126+
/// Evaluate the function and its derivative at the point @param x
127+
auto cyl_inters_func = [&h, &r, &sz, &sc](const scalar_t x) {
128+
const vector3_t crp = vector::cross(h.pos(x) - sc, sz);
129+
130+
// f(s) = ((h.pos(s) - sc) x sz)^2 - r^2 == 0
131+
const scalar_t f_s{(vector::dot(crp, crp) - r * r)};
132+
// f'(s) = 2 * ( (h.pos(s) - sc) x sz) * (h.dir(s) x sz) )
133+
const scalar_t df_s{2.f *
134+
vector::dot(crp, vector::cross(h.dir(x), sz))};
135+
136+
return std::make_tuple(f_s, df_s);
137+
};
138+
139+
for (unsigned int i = 0u; i < n_runs; ++i) {
140+
141+
const scalar_t &s_ini = paths[i];
142+
143+
// Run the root finding algorithm
144+
scalar_t s;
145+
scalar_t ds;
146+
if (run_rtsafe) {
147+
std::tie(s, ds) = newton_raphson_safe(cyl_inters_func, s_ini,
148+
convergence_tolerance,
149+
max_n_tries, max_path);
150+
} else {
151+
std::tie(s, ds) = newton_raphson(cyl_inters_func, s_ini,
152+
convergence_tolerance,
153+
max_n_tries, max_path);
247154
}
248155

249-
return ret;
156+
ret[i] = {s, h.pos(s), ds};
250157
}
158+
159+
return ret;
251160
}
252161

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

267-
using base_type =
268-
helix_intersector_impl<cylindrical2D<algebra_t>, algebra_t>;
269-
270176
using algebra_type = algebra_t;
271-
using scalar_type = dscalar<algebra_t>;
272-
using point3_type = dpoint3D<algebra_t>;
273-
using vector3_type = dvector3D<algebra_t>;
274-
using transform3_type = dtransform3D<algebra_t>;
275177

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

286188
using result_type = intersection_point_err<algebra_t>;
287189

288-
/// Operator function to find intersections between ray and planar mask
190+
/// Operator function to find intersections between helix and a
191+
/// concentric cylinder surface
289192
///
290-
/// @param ray is the input ray trajectory
291-
/// @param sf the surface handle the mask is associated with
292-
/// @param mask is the input mask that defines the surface extent
193+
/// @param h is the input helix trajectory
293194
/// @param trf is the surface placement transform
294-
/// @param mask_tolerance is the tolerance for mask edges
295-
/// @param overstep_tol negative cutoff for the path
195+
/// @param mask is the input mask that defines the surface extent
296196
///
297197
/// @return the intersection
298198
template <typename mask_t>
299199
DETRAY_HOST_DEVICE constexpr result_type point_of_intersection(
300-
const trajectory_type<algebra_t> &h, const transform3_type &trf,
301-
const mask_t &mask, const scalar_type = 0.f) const {
200+
const trajectory_type<algebra_t> &h, const dtransform3D<algebra_t> &trf,
201+
const mask_t &mask, const dscalar<algebra_t> = 0.f) const {
202+
203+
using base_t =
204+
helix_intersector_impl<cylindrical2D<algebra_t>, algebra_t>;
302205

303206
// Array of two solutions
304-
typename base_type::result_type results =
305-
base_type::point_of_intersection(h, trf, mask, 0.f);
207+
typename base_t::result_type results =
208+
base_t::point_of_intersection(h, trf, mask, 0.f);
306209

307210
// For portals, only take the solution along the helix direction,
308211
// not the one in the opposite direction

0 commit comments

Comments
 (0)