Skip to content

Commit 91a3606

Browse files
niermann999niermann999
committed
Remove raw newton solver in favour of rtsafe
1 parent 26680b6 commit 91a3606

File tree

5 files changed

+163
-462
lines changed

5 files changed

+163
-462
lines changed

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

Lines changed: 70 additions & 188 deletions
Original file line numberDiff line numberDiff line change
@@ -71,202 +71,86 @@ struct helix_intersector_impl<cylindrical2D<algebra_t>, algebra_t>
7171

7272
std::array<intersection_type<surface_descr_t>, 2> ret{};
7373

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

187-
return ret;
188-
} else {
189-
// Cylinder z axis
190-
const vector3_type sz = trf.z();
191-
// Cylinder centre
192-
const point3_type sc = trf.translation();
193-
194-
// Starting point on the helix for the Newton iteration
195-
// The mask is a cylinder -> it provides its radius as the first
196-
// value
197-
const scalar_type r{mask[cylinder2D::e_r]};
198-
199-
// Try to guess the best starting positions for the iteration
200-
201-
// Direction of the track at the helix origin
202-
const auto h_dir = h.dir(0.5f * r);
203-
// Default starting path length for the Newton iteration (assumes
204-
// concentric cylinder)
205-
const scalar_type default_s{r * getter::perp(h_dir)};
206-
207-
// Initial helix path length parameter
208-
std::array<scalar_type, 2> paths{default_s, default_s};
209-
210-
// try to guess good starting path by calculating the intersection
211-
// path of the helix tangential with the cylinder. This only has a
212-
// chance of working for tracks with reasonably high p_T !
213-
detail::ray<algebra_t> t{h.pos(), h.time(), h_dir, h.qop()};
214-
const auto qe = this->solve_intersection(t, mask, trf);
215-
216-
// Obtain both possible solutions by looping over the (different)
217-
// starting positions
218-
auto n_runs{static_cast<unsigned int>(qe.solutions())};
219-
220-
// Note: the default path length might be smaller than either
221-
// solution
222-
switch (qe.solutions()) {
223-
case 2:
224-
paths[1] = qe.larger();
225-
// If there are two solutions, reuse the case for a single
226-
// solution to setup the intersection with the smaller path
227-
// in ret[0]
228-
[[fallthrough]];
229-
case 1: {
230-
paths[0] = qe.smaller();
231-
break;
232-
}
233-
default: {
234-
n_runs = 2u;
235-
paths[0] = r;
236-
paths[1] = -r;
237-
}
238-
}
239-
240-
/// Evaluate the function and its derivative at the point @param x
241-
auto cyl_inters_func = [&h, &r, &sz, &sc](const scalar_type x) {
242-
const vector3_type crp = vector::cross(h.pos(x) - sc, sz);
243-
244-
// f(s) = ((h.pos(s) - sc) x sz)^2 - r^2 == 0
245-
const scalar_type f_s{(vector::dot(crp, crp) - r * r)};
246-
// f'(s) = 2 * ( (h.pos(s) - sc) x sz) * (h.dir(s) x sz) )
247-
const scalar_type df_s{
248-
2.f * vector::dot(crp, vector::cross(h.dir(x), sz))};
125+
/// Evaluate the function and its derivative at the point @param x
126+
auto cyl_inters_func = [&h, &r, &sz, &sc](const scalar_type x) {
127+
const vector3_type crp = vector::cross(h.pos(x) - sc, sz);
249128

250-
return std::make_tuple(f_s, df_s);
251-
};
129+
// f(s) = ((h.pos(s) - sc) x sz)^2 - r^2 == 0
130+
const scalar_type f_s{(vector::dot(crp, crp) - r * r)};
131+
// f'(s) = 2 * ( (h.pos(s) - sc) x sz) * (h.dir(s) x sz) )
132+
const scalar_type df_s{
133+
2.f * vector::dot(crp, vector::cross(h.dir(x), sz))};
252134

253-
for (unsigned int i = 0u; i < n_runs; ++i) {
135+
return std::make_tuple(f_s, df_s);
136+
};
254137

255-
const scalar_type &s_ini = paths[i];
256-
intersection_type<surface_descr_t> &sfi = ret[i];
138+
for (unsigned int i = 0u; i < n_runs; ++i) {
257139

258-
// Run the root finding algorithm
259-
const auto [s, ds] = newton_raphson_safe(cyl_inters_func, s_ini,
260-
convergence_tolerance,
261-
max_n_tries, max_path);
140+
const scalar_type &s_ini = paths[i];
141+
intersection_type<surface_descr_t> &sfi = ret[i];
262142

263-
// Build intersection struct from the root
264-
build_intersection(h, sfi, s, ds, sf_desc, mask, trf,
265-
mask_tolerance);
266-
}
143+
// Run the root finding algorithm
144+
const auto [s, ds] = newton_raphson_safe(cyl_inters_func, s_ini,
145+
convergence_tolerance,
146+
max_n_tries, max_path);
267147

268-
return ret;
148+
// Build intersection struct from the root
149+
build_intersection(h, sfi, s, ds, sf_desc, mask, trf,
150+
mask_tolerance);
269151
}
152+
153+
return ret;
270154
}
271155

272156
/// Interface to use fixed mask tolerance
@@ -286,8 +170,6 @@ struct helix_intersector_impl<cylindrical2D<algebra_t>, algebra_t>
286170
std::size_t max_n_tries{1000u};
287171
// Early exit, if the intersection is too far away
288172
scalar_type max_path{5.f * unit<scalar_type>::m};
289-
// Complement the Newton algorithm with Bisection steps
290-
bool run_rtsafe{true};
291173
};
292174

293175
template <typename algebra_t>

0 commit comments

Comments
 (0)