diff --git a/core/include/detray/builders/detector_builder.hpp b/core/include/detray/builders/detector_builder.hpp index 1f80a379f..96340e06e 100644 --- a/core/include/detray/builders/detector_builder.hpp +++ b/core/include/detray/builders/detector_builder.hpp @@ -119,8 +119,8 @@ class detector_builder { detector_type det{resource}; - DETRAY_INFO_HOST("Have " << m_volumes.size() - << " configured volume builders"); + DETRAY_VERBOSE_HOST("Have " << m_volumes.size() + << " configured volume builders"); DETRAY_VERBOSE_HOST("Start building the volumes..."); for (auto& vol_builder : m_volumes) { @@ -132,7 +132,6 @@ class detector_builder { // TODO: Add sorting, data deduplication etc. here later... - DETRAY_INFO_HOST("Detector building complete: " << name()); DETRAY_INFO_HOST("-> Built " << det.volumes().size() << " volumes"); DETRAY_INFO_HOST("-> Built " << det.surfaces().size() << " surfaces:"); DETRAY_INFO_HOST("--> portals: " << detray::n_portals(det)); @@ -155,6 +154,7 @@ class detector_builder { DETRAY_INFO_HOST("--> slabs: " << detray::n_material_slabs(det)); DETRAY_INFO_HOST("--> rods: " << detray::n_material_rods(det)); } + DETRAY_INFO_HOST("Detector building complete: " << name()); return det; } diff --git a/core/include/detray/core/detail/multi_store.hpp b/core/include/detray/core/detail/multi_store.hpp index d7fbf3f02..b7c5e5880 100644 --- a/core/include/detray/core/detail/multi_store.hpp +++ b/core/include/detray/core/detail/multi_store.hpp @@ -88,7 +88,7 @@ class multi_store { /// (host-side only) template > - requires(std::is_same_v> && + requires(std::is_same_v> && std::derived_from) DETRAY_HOST explicit multi_store(allocator_t &resource, const Ts &...args) : m_tuple_container(resource, args...) {} diff --git a/core/include/detray/definitions/actor.hpp b/core/include/detray/definitions/actor.hpp new file mode 100644 index 000000000..092739168 --- /dev/null +++ b/core/include/detray/definitions/actor.hpp @@ -0,0 +1,45 @@ +/** Detray library, part of the ACTS project (R&D line) + * + * (c) 2026 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" + +// System include(s) +#include +#include + +namespace detray::actor { + +enum class status : std::uint8_t { + e_notify = 0u, + e_success = 1u, + e_unknown = 2u, + e_failure = 3u, +}; + +// Print the values of an enum by identifier +#define ENUM_PRINT(x) \ + case x: \ + os << #x; \ + break + +DETRAY_HOST inline std::ostream& operator<<(std::ostream& os, status s) { + + switch (s) { + using enum status; + ENUM_PRINT(e_notify); + ENUM_PRINT(e_success); + ENUM_PRINT(e_unknown); + ENUM_PRINT(e_failure); + } + return os; +} + +#undef ENUM_PRINT +} // namespace detray::actor diff --git a/core/include/detray/materials/material_slab.hpp b/core/include/detray/materials/material_slab.hpp index 9069c8cfd..20b19c948 100644 --- a/core/include/detray/materials/material_slab.hpp +++ b/core/include/detray/materials/material_slab.hpp @@ -102,7 +102,7 @@ struct material_slab { const material_slab& mat) { os << "slab: "; os << mat.get_material(); - os << " | thickness = " << mat.thickness() << "mm"; + os << " | thickness = " << mat.thickness() << " mm"; return os; } diff --git a/core/include/detray/navigation/policies.hpp b/core/include/detray/navigation/policies.hpp index 72ae5f304..88b39b5c5 100644 --- a/core/include/detray/navigation/policies.hpp +++ b/core/include/detray/navigation/policies.hpp @@ -21,7 +21,7 @@ namespace detray { /// Struct that represents the most conservative navigation policy: alway re- /// initialize the current volume -struct always_init : actor { +struct always_init : public base_actor { struct state {}; @@ -38,7 +38,7 @@ struct always_init : actor { /// During guided navigation only the next surface should be re-evaluated. This /// maps to the 'high trust' level in the navigator -struct guided_navigation : actor { +struct guided_navigation : public base_actor { struct state {}; @@ -58,7 +58,7 @@ struct guided_navigation : actor { /// The reasoning is, that the track state might have changed much when a /// constraint was triggered. template -struct stepper_default_policy : actor { +struct stepper_default_policy : public base_actor { struct state { scalar_t tol{std::numeric_limits::epsilon()}; @@ -96,7 +96,7 @@ struct stepper_default_policy : actor { /// amount of step size correction as a measure for the change in direction of /// the track state. template -struct stepper_rk_policy : actor { +struct stepper_rk_policy : public base_actor { struct state { scalar_t m_threshold_fair_trust{0.05f}; diff --git a/core/include/detray/propagator/actor_chain.hpp b/core/include/detray/propagator/actor_chain.hpp index 0b91c1989..a54899996 100644 --- a/core/include/detray/propagator/actor_chain.hpp +++ b/core/include/detray/propagator/actor_chain.hpp @@ -38,11 +38,12 @@ class actor_chain { using actor_tuple = dtuple; // Tuple of actor states (including states of observing actors, if present) - using state_tuple = detail::tuple_cat_t...>; + using state_tuple = detail::unique_t< + detail::tuple_cat_t...>>; // Tuple of state references that is used in the propagator - using state_ref_tuple = - detail::tuple_cat_t...>; + using state_ref_tuple = detail::unique_t< + detail::tuple_cat_t...>>; /// Call all actors in the chain. /// @@ -63,10 +64,11 @@ class actor_chain { /// @returns a tuple of default constructible actor states DETRAY_HOST_DEVICE static constexpr auto make_default_actor_states() { - // Only possible if each state is default initializable - if constexpr ((std::default_initializable && - ...)) { - return state_tuple{}; + // Only possible if each state is default initializable (including + // obsevers to the actors in actors_t) + if constexpr (detail::tuple_all_v) { + return state_tuple(); } else { return std::nullopt; } @@ -93,7 +95,7 @@ class actor_chain { if constexpr (!typename actor_t::is_comp_actor()) { // No actor state defined (empty) if constexpr (std::same_as) { + detray::base_actor::state>) { actr(p_state); } else { actr(detail::get(states), p_state); diff --git a/core/include/detray/propagator/actors.hpp b/core/include/detray/propagator/actors.hpp index 235fb8219..3b88df088 100644 --- a/core/include/detray/propagator/actors.hpp +++ b/core/include/detray/propagator/actors.hpp @@ -10,8 +10,7 @@ // Include all core actors #include "detray/propagator/actor_chain.hpp" #include "detray/propagator/actors/aborters.hpp" -#include "detray/propagator/actors/parameter_resetter.hpp" -#include "detray/propagator/actors/parameter_transporter.hpp" +#include "detray/propagator/actors/parameter_updater.hpp" #include "detray/propagator/actors/pointwise_material_interactor.hpp" #include "detray/propagator/actors/surface_sequencer.hpp" #include "detray/propagator/concepts.hpp" diff --git a/core/include/detray/propagator/actors/aborters.hpp b/core/include/detray/propagator/actors/aborters.hpp index 079dc290e..bfcecfb8a 100644 --- a/core/include/detray/propagator/actors/aborters.hpp +++ b/core/include/detray/propagator/actors/aborters.hpp @@ -17,11 +17,11 @@ // System include(s) #include -namespace detray { +namespace detray::actor { /// Aborter that checks whether the track has exceeded its pathlimit template -struct pathlimit_aborter : actor { +struct pathlimit_aborter : public base_actor { /// Pathlimit for a single propagation workflow struct state { @@ -77,7 +77,7 @@ struct pathlimit_aborter : actor { /// Aborter that checks whether the track fell below a minimum momentum template -struct momentum_aborter : actor { +struct momentum_aborter : public base_actor { struct state { /// @returns the momentum limit. @@ -102,7 +102,7 @@ struct momentum_aborter : actor { scalar_t m_min_pT = 10.f * unit::MeV; }; - /// Enforces a minimum momentum magnitude + /// Actor interface: Enforces a minimum momentum magnitude /// /// @param abrt_state contains the momentum limit /// @param prop_state state of the propagation @@ -143,7 +143,7 @@ struct momentum_aborter : actor { }; /// Aborter checks whether a specific surface was reached -struct target_aborter : actor { +struct target_aborter : public base_actor { /// Keeps the index for the target surface struct state { @@ -174,4 +174,4 @@ struct target_aborter : actor { } }; -} // namespace detray +} // namespace detray::actor diff --git a/core/include/detray/propagator/actors/parameter_resetter.hpp b/core/include/detray/propagator/actors/parameter_resetter.hpp deleted file mode 100644 index ce24290f9..000000000 --- a/core/include/detray/propagator/actors/parameter_resetter.hpp +++ /dev/null @@ -1,83 +0,0 @@ -/** Detray library, part of the ACTS project (R&D line) - * - * (c) 2022-2025 CERN for the benefit of the ACTS project - * - * Mozilla Public License Version 2.0 - */ - -#pragma once - -// Project include(s). -#include "detray/definitions/algebra.hpp" -#include "detray/definitions/detail/qualifiers.hpp" -#include "detray/definitions/track_parametrization.hpp" -#include "detray/propagator/base_actor.hpp" -#include "detray/propagator/detail/jacobian_engine.hpp" -#include "detray/propagator/detail/noise_estimation.hpp" -#include "detray/propagator/propagation_config.hpp" -#include "detray/utils/logging.hpp" - -namespace detray { - -template -struct parameter_resetter : actor { - - struct state { - - /// Default construction - state() = default; - - /// Build from propagation configuration set by the user - DETRAY_HOST_DEVICE - explicit constexpr state(const propagation::config& cfg) - : accumulated_error{cfg.navigation.accumulated_error}, - n_stddev{cfg.navigation.n_scattering_stddev}, - estimate_scattering_noise{ - cfg.navigation.estimate_scattering_noise} {} - - /// Percentage of total track path to assume as accumulated error - float accumulated_error{0.001f}; - /// Number of standard deviations to assume to model the scattering - /// noise - int n_stddev{2}; - /// Estimate mask tolerance for navigation to for compensate scattering - bool estimate_scattering_noise{true}; - }; - - template - DETRAY_HOST_DEVICE void operator()(const state& resetter_state, - propagator_state_t& propagation) const { - - using scalar_t = dscalar; - - auto& navigation = propagation.navigation(); - auto& stepping = propagation.stepping(); - - // Do covariance transport when the track is on surface - if (!(navigation.is_on_sensitive() || - navigation.encountered_sf_material())) { - return; - } - DETRAY_VERBOSE_HOST_DEVICE("Actor: Update the free track parameters"); - - // Update free params after bound params were changed by actors - const auto sf = navigation.current_surface(); - const auto& bound_params = stepping.bound_params(); - stepping() = - sf.bound_to_free_vector(propagation.context(), bound_params); - assert(!stepping().is_invalid()); - - // Reset jacobian transport to identity matrix - stepping.reset_transport_jacobian(); - - // Track pos/dir is not known precisely: adjust navigation tolerances - if (resetter_state.estimate_scattering_noise) { - detail::estimate_external_mask_tolerance( - bound_params, propagation, - static_cast(resetter_state.n_stddev), - resetter_state.accumulated_error); - } - } -}; - -} // namespace detray diff --git a/core/include/detray/propagator/actors/parameter_transporter.hpp b/core/include/detray/propagator/actors/parameter_transporter.hpp deleted file mode 100644 index 2c5c4a30e..000000000 --- a/core/include/detray/propagator/actors/parameter_transporter.hpp +++ /dev/null @@ -1,300 +0,0 @@ -/** Detray library, part of the ACTS project (R&D line) - * - * (c) 2022-2024 CERN for the benefit of the ACTS project - * - * Mozilla Public License Version 2.0 - */ - -#pragma once - -// Project include(s). -#include "./codegen/covariance_transport.hpp" -#include "./codegen/full_jacobian.hpp" -#include "detray/algebra/utils/approximately_equal.hpp" -#include "detray/definitions/algebra.hpp" -#include "detray/definitions/detail/qualifiers.hpp" -#include "detray/definitions/track_parametrization.hpp" -#include "detray/geometry/tracking_surface.hpp" -#include "detray/propagator/base_actor.hpp" -#include "detray/propagator/detail/jacobian_engine.hpp" -#include "detray/utils/logging.hpp" -#include "detray/utils/type_registry.hpp" - -namespace detray { - -namespace detail { - -/// Filter the masks of a detector according to the local frame type -struct select_frame { - template - using type = typename mask_t::local_frame; -}; - -} // namespace detail - -template -struct parameter_transporter : actor { - - /// @name Type definitions for the struct - /// @{ - using scalar_type = dscalar; - // Transformation matching this struct - using transform3_type = dtransform3D; - // bound matrix type - using bound_matrix_t = bound_matrix; - // free matrix type - using free_matrix_t = free_matrix; - // Matrix type for bound to free jacobian - using bound_to_free_matrix_t = bound_to_free_matrix; - // Matrix type for free to bound jacobian - using free_to_bound_matrix_t = free_to_bound_matrix; - /// @} - - struct state { - bound_matrix_t* _full_jacobian_ptr = nullptr; - - DETRAY_HOST_DEVICE - explicit state(bound_matrix_t* full_jac = nullptr) - : _full_jacobian_ptr(full_jac) {} - }; - - struct get_bound_to_free_dpos_dloc_visitor { - template - DETRAY_HOST_DEVICE constexpr dmatrix operator()( - const frame_t& /*frame*/, const transform3_type& trf3, - const free_track_parameters& params) const { - - return detail::jacobian_engine:: - template bound_to_free_jacobian_submatrix_dpos_dloc( - trf3, params.pos(), params.dir()); - } - }; - - struct get_bound_to_free_dpos_dangle_visitor { - template - DETRAY_HOST_DEVICE constexpr dmatrix operator()( - const frame_t& /*frame*/, const transform3_type& trf3, - const free_track_parameters& params, - const dmatrix& ddir_dangle) const { - - return detail::jacobian_engine:: - template bound_to_free_jacobian_submatrix_dpos_dangle( - trf3, params.pos(), params.dir(), ddir_dangle); - } - }; - - struct get_free_to_bound_dloc_dpos_visitor { - template - DETRAY_HOST_DEVICE constexpr dmatrix operator()( - const frame_t& /*frame*/, const transform3_type& trf3, - const stepper_state_t& stepping) const { - - return detail::jacobian_engine:: - template free_to_bound_jacobian_submatrix_dloc_dpos( - trf3, stepping().pos(), stepping().dir()); - } - }; - - template - DETRAY_HOST_DEVICE void operator()(state& actor_state, - propagator_state_t& propagation) const { - auto& stepping = propagation.stepping(); - const auto& navigation = propagation.navigation(); - - // Do covariance transport when the track is on surface - if (!(navigation.is_on_sensitive() || - navigation.encountered_sf_material())) { - return; - } - DETRAY_VERBOSE_HOST_DEVICE( - "Actor: Transport track parameters to current surface"); - - // Geometry context for this track - const auto& gctx = propagation.context(); - - // Current Surface - const auto sf = navigation.current_surface(); - - // Bound track params of departure surface - auto& bound_params = stepping.bound_params(); - - // Covariance is transported only when the previous surface is an - // actual tracking surface. (i.e. This disables the covariance transport - // from curvilinear frame) - if (!bound_params.surface_link().is_invalid()) { - const auto full_jacobian = get_full_jacobian(propagation); - const bound_matrix_t old_cov = stepping.bound_params().covariance(); - - detail::transport_covariance_to_bound_impl( - old_cov, full_jacobian, stepping.bound_params().covariance()); - - if (actor_state._full_jacobian_ptr != nullptr) { - const auto aggregate_full_jacobian = - full_jacobian * (*(actor_state._full_jacobian_ptr)); - (*(actor_state._full_jacobian_ptr)) = aggregate_full_jacobian; - } - } - - // Convert free to bound vector - bound_params.set_parameter_vector( - sf.free_to_bound_vector(gctx, stepping())); - - // Set surface link - bound_params.set_surface_link(sf.barcode()); - - assert(!bound_params.is_invalid()); - } - - template - DETRAY_HOST_DEVICE constexpr bound_matrix_t get_full_jacobian( - propagator_state_t& propagation) const { - - // Map the surface shapes of the detector down to the common frames - using detector_t = typename propagator_state_t::detector_type; - using frame_registry_t = - types::mapped_registry; - - const auto& stepping = propagation.stepping(); - const auto& navigation = propagation.navigation(); - - // Geometry context for this track - const auto& gctx = propagation.context(); - - // Our goal here is to compute the full Jacobian, which is given as: - // - // J_full = J_F2B * (D + I) * J_transport * J_B2F - // - // The transport Jacobian is given, but the free-to-bound and - // bound-to-free matrices as well as the derivative matrix $D$ need to - // be computed still. In order to avoid full-rank matrix - // multiplications, we use the known substructure of the - // aforementioned matrices in order to perform fewer operations. We - // describe in detail the mathematics performed throughout the - // remainder of this function. - - // Current Surface - const auto sf = navigation.current_surface(); - - // Bound track params of departure surface - auto& bound_params = stepping.bound_params(); - - // Previous surface - tracking_surface prev_sf{navigation.detector(), - bound_params.surface_link()}; - - // Free track params of departure surface - const free_track_parameters free_params = - prev_sf.bound_to_free_vector(gctx, bound_params); - - // First, we will compute the bound-to-free Jacobian, which is given - // by three sub-Jacobians. In particular, the matrix has an 8x6 shape - // and looks like this: - // - // l0 l1 phi th q/p t - // px [[ A, A, B, B, 0, 0], - // py [ A, A, B, B, 0, 0], - // pz [ A, A, B, B, 0, 0], - // t [ 0, 0, 0, 0, 0, 1], - // dx [ 0, 0, C, C, 0, 0], - // dy [ 0, 0, C, C, 0, 0], - // dz [ 0, 0, 0, C, 0, 0], - // q/p [ 0, 0, 0, 0, 1, 0]] - // - // Note that we thus only have 17 out of 48 non-trivial matrix - // elements. - // - // In this matrix, A represents the d(pos)/d(loc) submatrix, B is the - // d(pos)/d(angle) submatrix, and C is the d(dir)/d(angle) submatrix, - // all of which are 3x2 in size. Also, A and B depend on the frame - // type while C is computed in the same way for all frames. Finally, - // note that submatrix B is the zero matrix for most frame types. - const auto& prev_trf3 = prev_sf.transform(gctx); - const dmatrix b2f_dpos_dloc = - types::visit( - prev_sf.shape_id(), prev_trf3, free_params); - - const dmatrix b2f_ddir_dangle = - detail::jacobian_engine:: - bound_to_free_jacobian_submatrix_ddir_dangle(bound_params); - - const dmatrix b2f_dpos_dangle = - types::visit( - prev_sf.shape_id(), prev_trf3, free_params, b2f_ddir_dangle); - - // Next, we compute the derivative which is defined as the outer - // product of the two 8x1 vectors representing the path to free - // derivative and the free to path derivative. Thus, it is the product - // of the following two vectors: - // - // px [[ A], [[ B],^T - // py [ A], [ B], - // pz [ A], [ B], - // t [ 0]] [ 0], - // dx [ A], [ C], - // dy [ A], [ C], - // dz [ A], [ C], - // q/p [ A], [ 0]] - // - // Where A is frame independent and non-zero, B is frame-dependent and - // non-zero, while C is frame-dependent and non-zero only for line - // frames. - auto vol = navigation.current_volume(); - const auto vol_mat_ptr = vol.has_material() - ? vol.material_parameters(stepping().pos()) - : nullptr; - - const auto path_to_free_derivative = - detail::jacobian_engine::path_to_free_derivative( - stepping().dir(), stepping.dtds(), - stepping.dqopds(vol_mat_ptr)); - - const auto free_to_path_derivative = sf.free_to_path_derivative( - gctx, stepping().pos(), stepping().dir(), stepping.dtds()); - - // Now, we compute the free-to-bound Jacobian which is of size 6x8 and - // has the following structure: - // - // px py pz t dx dy dz q/p - // l0 [[ A, A, A, 0, 0, 0, 0, 0], - // l1 [ A, A, A, 0, 0, 0, 0, 0], - // phi [ 0, 0, 0, 0, B, B, 0, 0], - // th [ 0, 0, 0, 0, B, B, B, 0], - // q/p [ 0, 0, 0, 0, 0, 0, 0, 1], - // t [ 0, 0, 0, 1, 0, 0, 0, 0]] - // - // Thus, the number of non-trivial elements is only 11 out of the 48 - // matrix elements. Also, submatrix A depends on the frame type while - // submatrix B is the same for all frame types. - const dmatrix f2b_dloc_dpos = - types::visit( - sf.shape_id(), sf.transform(gctx), propagation.stepping()); - - const dmatrix f2b_dangle_ddir = - detail::jacobian_engine:: - free_to_bound_jacobian_submatrix_dangle_ddir(stepping().dir()); - - // Finally, we can use our Sympy-generated full Jacobian computation - // and return its result. - bound_matrix_t full_jacobian; - - if constexpr (std::decay_t::stepper_uses_gradient) { - detail::update_full_jacobian_with_gradient_impl( - stepping.internal_transport_jacobian(), b2f_dpos_dloc, - b2f_ddir_dangle, b2f_dpos_dangle, path_to_free_derivative, - free_to_path_derivative, f2b_dloc_dpos, f2b_dangle_ddir, - full_jacobian); - } else { - detail::update_full_jacobian_without_gradient_impl( - stepping.internal_transport_jacobian(), b2f_dpos_dloc, - b2f_ddir_dangle, b2f_dpos_dangle, path_to_free_derivative, - free_to_path_derivative, f2b_dloc_dpos, f2b_dangle_ddir, - full_jacobian); - } - - return full_jacobian; - } -}; - -} // namespace detray diff --git a/core/include/detray/propagator/actors/parameter_updater.hpp b/core/include/detray/propagator/actors/parameter_updater.hpp new file mode 100644 index 000000000..4a4359f7a --- /dev/null +++ b/core/include/detray/propagator/actors/parameter_updater.hpp @@ -0,0 +1,636 @@ +/** Detray library, part of the ACTS project (R&D line) + * + * (c) 2026 CERN for the benefit of the ACTS project + * + * Mozilla Public License Version 2.0 + */ + +#pragma once + +// Project include(s). +#include "./codegen/covariance_transport.hpp" +#include "./codegen/full_jacobian.hpp" +#include "detray/definitions/algebra.hpp" +#include "detray/definitions/detail/qualifiers.hpp" +#include "detray/definitions/track_parametrization.hpp" +#include "detray/geometry/tracking_surface.hpp" +#include "detray/propagator/base_actor.hpp" +#include "detray/propagator/composite_actor.hpp" +#include "detray/propagator/detail/jacobian_engine.hpp" +#include "detray/propagator/detail/noise_estimation.hpp" +#include "detray/propagator/propagation_config.hpp" +#include "detray/utils/curvilinear_frame.hpp" +#include "detray/utils/logging.hpp" +#include "detray/utils/type_registry.hpp" + +namespace detray::actor { + +template +struct parameter_transporter; + +template +struct parameter_setter; + +/// Configuration of noise estimation +struct noise_cfg { + /// Percentage of total track path to assume as accumulated error + float accumulated_error{0.001f}; + /// Number of standard deviations to assume to model the scattering noise + int n_stddev{2}; + /// Estimate mask tolerance for navigation to for compensate scattering + bool estimate_scattering_noise{true}; +}; + +/// State of the parameter updater. Used by two distinct steps: Transporter and +/// setter +template +struct parameter_updater_state { + + // Allow only the corresponding actors to modify the parameter data + friend struct parameter_transporter; + friend struct parameter_setter; + + constexpr parameter_updater_state() = default; + + /// Start without explicit track parameters + DETRAY_HOST_DEVICE + explicit constexpr parameter_updater_state( + const propagation::config& cfg, + bound_matrix* full_jac = nullptr) + : m_full_jacobian(full_jac), + m_cfg{cfg.navigation.accumulated_error, + cfg.navigation.n_scattering_stddev, + cfg.navigation.estimate_scattering_noise} {} + + /// Start from free track parameters + DETRAY_HOST_DEVICE + constexpr parameter_updater_state( + const propagation::config& cfg, + const free_track_parameters& free_params, + bound_matrix* full_jac = nullptr) + : m_full_jacobian(full_jac), + m_cfg{cfg.navigation.accumulated_error, + cfg.navigation.n_scattering_stddev, + cfg.navigation.estimate_scattering_noise} { + // Set bound track parameters + curvilinear_frame cf(free_params); + m_bound_params.set_parameter_vector(cf.m_bound_vec); + } + + /// Start from bound track parameters + DETRAY_HOST_DEVICE + constexpr parameter_updater_state( + const propagation::config& cfg, + const bound_track_parameters& bound_params, + bound_matrix* full_jac = nullptr) + : m_bound_params{bound_params}, + m_full_jacobian(full_jac), + m_cfg{cfg.navigation.accumulated_error, + cfg.navigation.n_scattering_stddev, + cfg.navigation.estimate_scattering_noise} {} + + /// Initialize the state from free track parameters + DETRAY_HOST_DEVICE + constexpr void init(const free_track_parameters& free_params) { + // Set bound track parameters + curvilinear_frame cf(free_params); + m_bound_params.set_parameter_vector(cf.m_bound_vec); + + // A dummy covariance - should not be used + m_bound_params.set_covariance( + matrix::identity>()); + + // An invalid barcode - should not be used + m_bound_params.set_surface_link(geometry::barcode{}); + + assert(!m_bound_params.is_invalid()); + } + + /// Initialize the state from bound track parameters + DETRAY_HOST_DEVICE + constexpr void init(const bound_track_parameters& bound_params) { + assert(!bound_params.is_invalid()); + m_bound_params = bound_params; + } + + /// Always update the track parameters in the stepper state. Otherwise, + /// observing actors decide if the paramters should be updated. + DETRAY_HOST_DEVICE + void always_update(const bool do_update = true) { + m_always_update = do_update; + } + + /// Notify the observing actors on the initial propagation surface (i.e. + /// before the stepper advanced the state from the seed parameters) + DETRAY_HOST_DEVICE + void notify_on_initial(const bool do_notify = true) { + m_notify_on_initial = do_notify; + } + + /// @return access to the noise estimation configuration - const + DETRAY_HOST_DEVICE + const noise_cfg& noise_estimation_cfg() const { return m_cfg; } + + /// @return access to the noise estimation configuration + DETRAY_HOST_DEVICE + noise_cfg& noise_estimation_cfg() { return m_cfg; } + + /// @returns bound track parameters - const + /// @note Only meaningful, if the navigation is at the corresponding surface + DETRAY_HOST_DEVICE + constexpr const bound_track_parameters& bound_params() const { + return m_bound_params; + } + + /// @returns bound track parameters + /// @note Only meaningful, if the navigation is at the corresponding surface + DETRAY_HOST_DEVICE + constexpr bound_track_parameters& bound_params() { + return m_bound_params; + } + + /// @returns true if the full Jacobian matrix should be assembled. + DETRAY_HOST_DEVICE + constexpr bool has_full_jacobian() const { + return m_full_jacobian != nullptr; + } + + /// @returns the current full Jacbian. + DETRAY_HOST_DEVICE + constexpr const bound_matrix& full_jacobian() const { + assert(has_full_jacobian()); + return *m_full_jacobian; + } + + /// Set new full Jacbian. + DETRAY_HOST_DEVICE + constexpr void set_full_jacobian(const bound_matrix& jac) { + assert(has_full_jacobian()); + *m_full_jacobian = jac; + } + + /// Set new full Jacbian. + DETRAY_HOST_DEVICE + constexpr void set_full_jacobian(bound_matrix* jac_ptr) { + assert(jac_ptr); + m_full_jacobian = jac_ptr; + } + + private: + /// Bound track parameters and covariance + bound_track_parameters m_bound_params{}; + /// Full jacobian for up to the current destination surface + bound_matrix* m_full_jacobian{nullptr}; + /// Configuration for the noise estimation + noise_cfg m_cfg{}; + /// Always update the free track parameters, even if nothing changed + bool m_always_update{false}; + /// Notify observing actors on initial surface (propagation init) + bool m_notify_on_initial{true}; +}; + +/// Result of the param. transporter: bound track parameters at dest. sf. +template +struct parameter_transporter_result : public actor::result { + + constexpr parameter_transporter_result() = default; + + /// Parameterized constructor using a pointer to the destination parameters + DETRAY_HOST_DEVICE + constexpr parameter_transporter_result( + actor::status stat, bound_track_parameters* params_ptr) + : actor::result{stat}, m_destination_params_ptr{params_ptr} { + assert(m_destination_params_ptr); + } + + /// Parameterized constructor using destination parameters + DETRAY_HOST_DEVICE + constexpr parameter_transporter_result( + actor::status stat, bound_track_parameters& params) + : actor::result{stat}, m_destination_params_ptr{¶ms} { + assert(m_destination_params_ptr); + } + + /// @returns access to the destination bound track parameters - const + DETRAY_HOST_DEVICE + constexpr const bound_track_parameters& destination_params() + const { + assert(this->status == actor::status::e_unknown || + m_destination_params_ptr); + return *m_destination_params_ptr; + } + + /// @returns access to the destination bound track parameters + DETRAY_HOST_DEVICE + constexpr bound_track_parameters& destination_params() { + assert(this->status == actor::status::e_unknown || + m_destination_params_ptr); + return *m_destination_params_ptr; + } + + private: + /// @returns a string stream that prints the transporter result details + DETRAY_HOST + friend std::ostream& operator<<( + std::ostream& os, const parameter_transporter_result& res) { + os << static_cast(res) << std::endl; + os << "destination params:\n" << res.destination_params() << std::endl; + return os; + } + + /// Bound track parameters of destination surface + bound_track_parameters* m_destination_params_ptr{nullptr}; +}; + +/// Transport the free track parameters at the destination surface +/// and the covariance at the departure surface to bound track parameters +/// at the destination surface (current sensitive/material surface) +template +struct parameter_transporter : base_actor { + + /// @name Type definitions for the struct + /// @{ + // Transformation matching this struct + using transform3_type = dtransform3D; + // The current free track parameters in the stepper algorithm + using free_track_parameters_type = free_track_parameters; + // The track parameters bound to the current sensitive/material surface + using bound_track_parameters_type = bound_track_parameters; + // Bound matrix type (bound covariance) + using bound_matrix_type = bound_matrix; + /// @} + + /// Use the parameter updater state + using state = parameter_updater_state; + using result = parameter_transporter_result; + + /// Filter the masks of a detector according to the local frame type + struct select_frame { + template + using type = typename mask_t::local_frame; + }; + + /// Visitors to the surface local coordinate frame + /// @{ + struct get_bound_to_free_dpos_dloc_visitor { + template + DETRAY_HOST_DEVICE constexpr dmatrix operator()( + const frame_t& /*frame*/, const transform3_type& trf3, + const free_track_parameters_type& params) const { + + return detail::jacobian_engine:: + template bound_to_free_jacobian_submatrix_dpos_dloc( + trf3, params.pos(), params.dir()); + } + }; + + struct get_bound_to_free_dpos_dangle_visitor { + template + DETRAY_HOST_DEVICE constexpr dmatrix operator()( + const frame_t& /*frame*/, const transform3_type& trf3, + const free_track_parameters_type& params, + const dmatrix& ddir_dangle) const { + + return detail::jacobian_engine:: + template bound_to_free_jacobian_submatrix_dpos_dangle( + trf3, params.pos(), params.dir(), ddir_dangle); + } + }; + + struct get_free_to_bound_dloc_dpos_visitor { + template + DETRAY_HOST_DEVICE constexpr dmatrix operator()( + const frame_t& /*frame*/, const transform3_type& trf3, + const stepper_state_t& stepping) const { + + return detail::jacobian_engine:: + template free_to_bound_jacobian_submatrix_dloc_dpos( + trf3, stepping().pos(), stepping().dir()); + } + }; + /// @} + + /// Actor interface + template + DETRAY_HOST_DEVICE result + operator()(state& updater_state, propagator_state_t& propagation) const { + + const auto& navigation = propagation.navigation(); + + // Do covariance transport only when the track is on surface + if (!(navigation.is_on_sensitive() || + navigation.encountered_sf_material())) { + return {}; + } + assert(navigation.is_on_surface()); + + // Current track position and transport jacobian + auto& stepping = propagation.stepping(); + + // Destination surface (current) + const auto dest_sf = navigation.current_surface(); + + // Bound track params of departure surface (not yet updated) + auto& departure_params = updater_state.bound_params(); + + // Result to be passed on to observing actors + result res{actor::status::e_notify, departure_params}; + + // Covariance is transported only when the departure surface is an + // actual tracking surface (i.e. this disables the covariance + // transport from curvilinear frames). + if (!departure_params.surface_link().is_invalid()) { + + DETRAY_DEBUG_HOST("Actor: Departure surface: " + << departure_params.surface_link()); + + // There is no need to transport the identical/initial parameters + assert(!dest_sf.barcode().is_invalid()); + if (departure_params.surface_link() == dest_sf.barcode()) { + DETRAY_VERBOSE_HOST_DEVICE( + "Actor: On initial surface (%d), parameter transport not " + "required", + dest_sf.index()); + + // Notify observers? + if (!updater_state.m_notify_on_initial) { + res.status = actor::status::e_unknown; + } + + return res; + } + + DETRAY_DEBUG_HOST( + "Actor: Destination surface: " << dest_sf.barcode()); + DETRAY_VERBOSE_HOST_DEVICE( + "Actor: Transport track covariance to surface %d", + dest_sf.index()); + + // Transport the covariance + const bound_matrix propagation_step_jacobian = + get_full_jacobian(propagation, departure_params); + + // Update the full Jacobian, if required + if (math::fabs(stepping.path_length()) > 0.f) { + + if (updater_state.has_full_jacobian()) { + updater_state.set_full_jacobian( + propagation_step_jacobian * + updater_state.full_jacobian()); + } + + // Reset transport Jacobian to identity matrix + stepping.reset_transport_jacobian(); + } + + // Copy old covariance in order to update new covariance in place + const bound_matrix_type old_cov = departure_params.covariance(); + bound_matrix_type& new_cov = departure_params.covariance(); + + detray::detail::transport_covariance_to_bound_impl( + old_cov, propagation_step_jacobian, new_cov); + } else { + DETRAY_VERBOSE_HOST_DEVICE( + "Actor: Departure surface link invalid, setting covariance to " + "identity"); + + // Dummy covariance that allows to multiply the transport jacobian + departure_params.set_covariance( + matrix::identity()); + } + + DETRAY_VERBOSE_HOST_DEVICE( + "Actor: Convert track parameter vector to local for surface %d", + dest_sf.index()); + + // Get the bound parameters at the destination surface + res.destination_params().set_parameter_vector( + dest_sf.free_to_bound_vector(propagation.context(), stepping())); + + // Set new surface link + res.destination_params().set_surface_link(dest_sf.barcode()); + + assert(!res.destination_params().is_invalid()); + + DETRAY_DEBUG_HOST("Actor: Transported bound param.:\n" + << res.destination_params()); + + return res; + } + + /// @returns the full Jacobian between departure and destination surfaces + template + DETRAY_HOST_DEVICE constexpr bound_matrix_type get_full_jacobian( + propagator_state_t& propagation, + const bound_track_parameters_type& departure_params) const { + + // Map the surface shapes of the detector down to the common frames + using detector_t = typename propagator_state_t::detector_type; + using frame_registry_t = + types::mapped_registry; + + const auto& gctx = propagation.context(); + const auto& stepping = propagation.stepping(); + const auto& navigation = propagation.navigation(); + + // Our goal here is to compute the full Jacobian, which is given as: + // + // J_full = J_F2B * (D + I) * J_transport * J_B2F + // + // The transport Jacobian is given, but the free-to-bound and + // bound-to-free matrices as well as the derivative matrix $D$ need + // to be computed still. In order to avoid full-rank matrix + // multiplications, we use the known substructure of the + // aforementioned matrices in order to perform fewer operations. We + // describe in detail the mathematics performed throughout the + // remainder of this function. + + // Departure surface + const tracking_surface dep_sf{navigation.detector(), + departure_params.surface_link()}; + + // Destination Surface + const tracking_surface dest_sf = navigation.current_surface(); + + // Free track params of departure surface + const free_track_parameters dep_free_params = + dep_sf.bound_to_free_vector(gctx, departure_params); + + // Free track params of destination surface + const free_track_parameters& dest_free_params = stepping(); + + // First, we will compute the bound-to-free Jacobian, which is given + // by three sub-Jacobians. In particular, the matrix has an 8x6 + // shape and looks like this: + // + // l0 l1 phi th q/p t + // px [[ A, A, B, B, 0, 0], + // py [ A, A, B, B, 0, 0], + // pz [ A, A, B, B, 0, 0], + // t [ 0, 0, 0, 0, 0, 1], + // dx [ 0, 0, C, C, 0, 0], + // dy [ 0, 0, C, C, 0, 0], + // dz [ 0, 0, 0, C, 0, 0], + // q/p [ 0, 0, 0, 0, 1, 0]] + // + // Note that we thus only have 17 out of 48 non-trivial matrix + // elements. + // + // In this matrix, A represents the d(pos)/d(loc) submatrix, B is + // the d(pos)/d(angle) submatrix, and C is the d(dir)/d(angle) + // submatrix, all of which are 3x2 in size. Also, A and B depend on + // the frame type while C is computed in the same way for all + // frames. Finally, note that submatrix B is the zero matrix for + // most frame types. + const auto& dep_trf3 = dep_sf.transform(gctx); + const dmatrix b2f_dpos_dloc = + types::visit( + dep_sf.shape_id(), dep_trf3, dep_free_params); + + const dmatrix b2f_ddir_dangle = + detail::jacobian_engine:: + bound_to_free_jacobian_submatrix_ddir_dangle(departure_params); + + const dmatrix b2f_dpos_dangle = + types::visit( + dep_sf.shape_id(), dep_trf3, dep_free_params, b2f_ddir_dangle); + + // Next, we compute the derivative which is defined as the outer + // product of the two 8x1 vectors representing the path to free + // derivative and the free to path derivative. Thus, it is the + // product of the following two vectors: + // + // px [[ A], [[ B],^T + // py [ A], [ B], + // pz [ A], [ B], + // t [ 0]] [ 0], + // dx [ A], [ C], + // dy [ A], [ C], + // dz [ A], [ C], + // q/p [ A], [ 0]] + // + // Where A is frame independent and non-zero, B is frame-dependent + // and non-zero, while C is frame-dependent and non-zero only for + // line frames. + const dpoint3D dest_glob_pos{dest_free_params.pos()}; + const dvector3D dest_glob_dir{dest_free_params.dir()}; + + auto vol = navigation.current_volume(); + const auto vol_mat_ptr = vol.has_material() + ? vol.material_parameters(dest_glob_pos) + : nullptr; + + const auto path_to_free_derivative = + detail::jacobian_engine::path_to_free_derivative( + dest_glob_dir, stepping.dtds(), stepping.dqopds(vol_mat_ptr)); + + const auto free_to_path_derivative = dest_sf.free_to_path_derivative( + gctx, dest_glob_pos, dest_glob_dir, stepping.dtds()); + + // Now, we compute the free-to-bound Jacobian which is of size 6x8 + // and has the following structure: + // + // px py pz t dx dy dz q/p + // l0 [[ A, A, A, 0, 0, 0, 0, 0], + // l1 [ A, A, A, 0, 0, 0, 0, 0], + // phi [ 0, 0, 0, 0, B, B, 0, 0], + // th [ 0, 0, 0, 0, B, B, B, 0], + // q/p [ 0, 0, 0, 0, 0, 0, 0, 1], + // t [ 0, 0, 0, 1, 0, 0, 0, 0]] + // + // Thus, the number of non-trivial elements is only 11 out of the 48 + // matrix elements. Also, submatrix A depends on the frame type + // while submatrix B is the same for all frame types. + const dmatrix f2b_dloc_dpos = + types::visit( + dest_sf.shape_id(), dest_sf.transform(gctx), stepping); + + const dmatrix f2b_dangle_ddir = + detail::jacobian_engine:: + free_to_bound_jacobian_submatrix_dangle_ddir(dest_glob_dir); + + // Finally, we can use our Sympy-generated full Jacobian computation + // and return its result. + bound_matrix_type full_jacobian; + + if constexpr (std::decay_t::stepper_uses_gradient) { + detail::update_full_jacobian_with_gradient_impl( + stepping.internal_transport_jacobian(), b2f_dpos_dloc, + b2f_ddir_dangle, b2f_dpos_dangle, path_to_free_derivative, + free_to_path_derivative, f2b_dloc_dpos, f2b_dangle_ddir, + full_jacobian); + } else { + detail::update_full_jacobian_without_gradient_impl( + stepping.internal_transport_jacobian(), b2f_dpos_dloc, + b2f_ddir_dangle, b2f_dpos_dangle, path_to_free_derivative, + free_to_path_derivative, f2b_dloc_dpos, f2b_dangle_ddir, + full_jacobian); + } + + return full_jacobian; + } +}; + +/// Set the free track parameters in the stepper state, in case +/// observing actors (e.g. the material interactor) changed the bound parameters +template +struct parameter_setter : base_actor { + + /// Acces the same parameter updater state as the parameter transporter + using state = parameter_updater_state; + + /// Actor interface: Observer to the parameter transporter and its observers + template + DETRAY_HOST_DEVICE void operator()( + state& updater_state, propagator_state_t& propagation, + const parameter_transporter_result& res) const { + + // Update the track parameters only if necessary + if (res.status != actor::status::e_success && + !updater_state.m_always_update) { + return; + } + + const auto& navigation = propagation.navigation(); + auto& stepping = propagation.stepping(); + + // One of the transporter observers might have exited the navigation + assert(navigation.is_on_surface() || !navigation.is_alive()); + + DETRAY_VERBOSE_HOST_DEVICE("Actor: Update the track parameters"); + + // Updated bound track parameters + const bound_track_parameters& bound_params = + res.destination_params(); + + // Update free params after bound params were changed by actors + const tracking_surface dest_sf{navigation.detector(), + bound_params.surface_link()}; + stepping() = + dest_sf.bound_to_free_vector(propagation.context(), bound_params); + + assert(!stepping().is_invalid()); + + // Track pos/dir is not always known precisely: adjust navigation + // tolerances according to bound covariance + if (const noise_cfg& cfg = updater_state.noise_estimation_cfg(); + cfg.estimate_scattering_noise) { + detail::estimate_external_mask_tolerance( + bound_params, propagation, + static_cast>(cfg.n_stddev), + cfg.accumulated_error); + } + + DETRAY_DEBUG_HOST("Actor: Updated bound param.:\n" << bound_params); + } +}; + +/// Call actors that depend on the bound track parameters safely together +/// with the parameter transporter and parameter setter +template +using parameter_updater = + composite_actor, transporter_observers..., + parameter_setter>; + +} // namespace detray::actor diff --git a/core/include/detray/propagator/actors/pointwise_material_interactor.hpp b/core/include/detray/propagator/actors/pointwise_material_interactor.hpp index 0f9c9db7e..09fdf5086 100644 --- a/core/include/detray/propagator/actors/pointwise_material_interactor.hpp +++ b/core/include/detray/propagator/actors/pointwise_material_interactor.hpp @@ -8,21 +8,23 @@ #pragma once // Project include(s). +#include "detray/definitions/actor.hpp" #include "detray/definitions/algebra.hpp" #include "detray/definitions/detail/qualifiers.hpp" #include "detray/definitions/track_parametrization.hpp" #include "detray/materials/concepts.hpp" #include "detray/materials/detail/material_accessor.hpp" #include "detray/materials/interaction.hpp" +#include "detray/propagator/actors/parameter_updater.hpp" #include "detray/propagator/base_actor.hpp" #include "detray/tracks/bound_track_parameters.hpp" #include "detray/utils/geometry_utils.hpp" #include "detray/utils/logging.hpp" -namespace detray { +namespace detray::actor { template -struct pointwise_material_interactor : actor { +struct pointwise_material_interactor : public base_actor { using algebra_type = algebra_t; using scalar_type = dscalar; @@ -132,24 +134,29 @@ struct pointwise_material_interactor : actor { template DETRAY_HOST_DEVICE inline void operator()( - state &interactor_state, propagator_state_t &prop_state) const { - - interactor_state.reset(); + state &interactor_state, propagator_state_t &prop_state, + parameter_transporter_result &res) const { const auto &navigation = prop_state.navigation(); // Do material interaction when the track is on material surface - if (navigation.encountered_sf_material()) { + if (!navigation.encountered_sf_material()) { + return; + } - DETRAY_VERBOSE_HOST_DEVICE("Actor: Resolve material effects:"); + DETRAY_VERBOSE_HOST_DEVICE("Actor: Resolve material effects:"); + + interactor_state.reset(); - auto &stepping = prop_state.stepping(); + const auto &stepping = prop_state.stepping(); + const bool success = this->update(prop_state.context(), stepping.particle_hypothesis(), - stepping.bound_params(), interactor_state, + res.destination_params(), interactor_state, static_cast(navigation.direction()), navigation.current_surface()); - } + + res.status = success ? actor::status::e_success : res.status; } /// @brief Update the bound track parameter @@ -159,7 +166,7 @@ struct pointwise_material_interactor : actor { /// @param[in] nav_dir navigation direction /// @param[in] sf the surface template - DETRAY_HOST_DEVICE inline void update( + DETRAY_HOST_DEVICE inline bool update( const context_t gctx, const pdg_particle &ptc, bound_track_parameters &bound_params, state &interactor_state, const int nav_dir, const surface_t &sf) const { @@ -170,10 +177,10 @@ struct pointwise_material_interactor : actor { const scalar_type cos_inc_angle{cos_angle(gctx, sf, bound_params.dir(), bound_params.bound_local())}; - const bool succeed = sf.template visit_material( + const bool success = sf.template visit_material( interactor_state, ptc, bound_params, cos_inc_angle, approach); - if (succeed) { + if (success) { auto &covariance = bound_params.covariance(); @@ -196,6 +203,8 @@ struct pointwise_material_interactor : actor { } assert(!bound_params.is_invalid()); + + return success; } /// @brief Update the q over p of bound track parameter @@ -219,7 +228,7 @@ struct pointwise_material_interactor : actor { math::sqrt(m * m + p * p) - math::copysign(e_loss, static_cast(sign))}; - DETRAY_DEBUG_HOST_DEVICE("-> new energy: %f", next_e); + DETRAY_DEBUG_HOST_DEVICE("-> new energy: %f GeV", next_e); // Put particle at rest if energy loss is too large const scalar_type next_p{ @@ -230,10 +239,11 @@ struct pointwise_material_interactor : actor { const scalar_type next_qop{(q != 0.f) ? q / next_p : 1.f / next_p}; vector.set_qop((next_p == 0.f) ? inv : next_qop); - DETRAY_DEBUG_HOST_DEVICE("-> new abs. momentum: %f", next_p); + DETRAY_DEBUG_HOST_DEVICE("-> new abs. momentum: %f GeV", next_p); - DETRAY_DEBUG_HOST_DEVICE("-> Update qop: before: %f, after: %f", q / p, - vector.qop()); + DETRAY_DEBUG_HOST_DEVICE( + "-> Update qop: before: %f e/GeV, after: %f e/GeV", q / p, + vector.qop()); } /// @brief Update the variance of q over p of bound track parameter @@ -246,13 +256,14 @@ struct pointwise_material_interactor : actor { const scalar_type variance_qop{sigma_qop * sigma_qop}; - DETRAY_DEBUG_HOST_DEVICE("-> qop variance: %f", variance_qop); + DETRAY_DEBUG_HOST_DEVICE("-> qop variance: %f (e/GeV)^2", variance_qop); getter::element(covariance, e_bound_qoverp, e_bound_qoverp) += variance_qop; DETRAY_DEBUG_HOST_DEVICE( - "-> Update qop variance: before: %f, after: %f", + "-> Update qop variance: before: %f (e/GeV)^2, after: %f " + "(e/GeV)^2", getter::element(covariance, e_bound_qoverp, e_bound_qoverp) - variance_qop, getter::element(covariance, e_bound_qoverp, e_bound_qoverp)); @@ -271,14 +282,14 @@ struct pointwise_material_interactor : actor { const scalar_type var_scattering_angle{projected_scattering_angle * projected_scattering_angle}; - DETRAY_DEBUG_HOST_DEVICE("-> Scattering angle variance: %f", + DETRAY_DEBUG_HOST_DEVICE("-> Scattering angle variance: %f rad^2", var_scattering_angle); constexpr auto inv{detail::invalid_value()}; DETRAY_DEBUG_HOST_DEVICE("-> Update phi variance:"); DETRAY_DEBUG_HOST_DEVICE( - "--> before: %f", + "--> before: %f rad^2", getter::element(covariance, e_bound_phi, e_bound_phi)); getter::element(covariance, e_bound_phi, e_bound_phi) += @@ -286,21 +297,21 @@ struct pointwise_material_interactor : actor { : var_scattering_angle / (1.f - dir[2] * dir[2]); DETRAY_DEBUG_HOST_DEVICE( - "--> after: %f", + "--> after: %f rad^2", getter::element(covariance, e_bound_phi, e_bound_phi)); DETRAY_DEBUG_HOST_DEVICE("-> Update theta variance:"); DETRAY_DEBUG_HOST_DEVICE( - "--> before: %f", + "--> before: %f rad^2", getter::element(covariance, e_bound_theta, e_bound_theta)); getter::element(covariance, e_bound_theta, e_bound_theta) += var_scattering_angle; DETRAY_DEBUG_HOST_DEVICE( - "--> after: %f", + "--> after: %f rad^2", getter::element(covariance, e_bound_theta, e_bound_theta)); } }; -} // namespace detray +} // namespace detray::actor diff --git a/core/include/detray/propagator/actors/surface_sequencer.hpp b/core/include/detray/propagator/actors/surface_sequencer.hpp index 4fd605b26..ba3bcedc6 100644 --- a/core/include/detray/propagator/actors/surface_sequencer.hpp +++ b/core/include/detray/propagator/actors/surface_sequencer.hpp @@ -1,6 +1,6 @@ /** Detray library, part of the ACTS project (R&D line) * - * (c) 2025 CERN for the benefit of the ACTS project + * (c) 2025-2026 CERN for the benefit of the ACTS project * * Mozilla Public License Version 2.0 */ @@ -16,10 +16,13 @@ // Vecmem include(s) #include -namespace detray { +// System includes +#include + +namespace detray::actor { template -struct surface_sequencer : actor { +struct surface_sequencer : public base_actor { struct state { using surface_type = sf_descriptor_t; @@ -102,4 +105,4 @@ struct surface_sequencer : actor { } }; -} // namespace detray +} // namespace detray::actor diff --git a/core/include/detray/propagator/base_actor.hpp b/core/include/detray/propagator/base_actor.hpp index f239b1f51..0d1731f3e 100644 --- a/core/include/detray/propagator/base_actor.hpp +++ b/core/include/detray/propagator/base_actor.hpp @@ -1,19 +1,22 @@ /** Detray library, part of the ACTS project (R&D line) * - * (c) 2022-2025 CERN for the benefit of the ACTS project + * (c) 2022-2026 CERN for the benefit of the ACTS project * * Mozilla Public License Version 2.0 */ #pragma once +// Project include(s) +#include "detray/definitions/actor.hpp" + // System include(s) #include namespace detray { /// Base class actor implementation -struct actor { +struct base_actor { /// Tag whether this is a composite type struct is_comp_actor : public std::false_type {}; @@ -21,4 +24,23 @@ struct actor { struct state {}; }; +namespace actor { + +/// Dummy result type that signals that no result is present +struct empty_result {}; + +/// Result of the principal actor to be passed to the obsevers +struct result { + actor::status status{actor::status::e_unknown}; + + /// @returns a string stream that prints the transporter result details + DETRAY_HOST + friend std::ostream &operator<<(std::ostream &os, const result &res) { + os << "status: " << res.status << std::endl; + return os; + } +}; + +} // namespace actor + } // namespace detray diff --git a/core/include/detray/propagator/base_stepper.hpp b/core/include/detray/propagator/base_stepper.hpp index eab95cfcf..6adc00143 100644 --- a/core/include/detray/propagator/base_stepper.hpp +++ b/core/include/detray/propagator/base_stepper.hpp @@ -68,20 +68,7 @@ class base_stepper { DETRAY_HOST_DEVICE explicit state(const free_track_parameters_type &free_params) : m_track(free_params) { - - curvilinear_frame cf(free_params); - - // Set bound track parameters - m_bound_params.set_parameter_vector(cf.m_bound_vec); - - // A dummy covariance - should not be used - m_bound_params.set_covariance( - matrix::identity()); - - // An invalid barcode - should not be used - m_bound_params.set_surface_link(geometry::barcode{}); - - assert(!m_bound_params.is_invalid()); + assert(!m_track.is_invalid()); // HACK: When the overload resolution for the transport Jacobian // type is resolved, turn this into a default member @@ -99,11 +86,10 @@ class base_stepper { DETRAY_HOST_DEVICE state( const bound_track_parameters_type &bound_params, const detector_t &det, - const typename detector_t::geometry_context &ctx) - : m_bound_params(bound_params) { + const typename detector_t::geometry_context &ctx) { - assert(!m_bound_params.is_invalid()); - assert(!m_bound_params.surface_link().is_invalid()); + assert(!bound_params.is_invalid()); + assert(!bound_params.surface_link().is_invalid()); // Departure surface const auto sf = tracking_surface{det, bound_params.surface_link()}; @@ -132,16 +118,6 @@ class base_stepper { DETRAY_HOST_DEVICE const free_track_parameters_type &operator()() const { return m_track; } - /// @returns bound track parameters - const access - DETRAY_HOST_DEVICE - bound_track_parameters_type &bound_params() { return m_bound_params; } - - /// @returns bound track parameters - non-const access - DETRAY_HOST_DEVICE - const bound_track_parameters_type &bound_params() const { - return m_bound_params; - } - /// Get stepping direction DETRAY_HOST_DEVICE inline step::direction direction() const { @@ -288,9 +264,6 @@ class base_stepper { /// Jacobian transport matrix internal_jacobian_type m_jac_transport; - /// Bound covariance - bound_track_parameters_type m_bound_params; - /// Free track parameters free_track_parameters_type m_track; diff --git a/core/include/detray/propagator/composite_actor.hpp b/core/include/detray/propagator/composite_actor.hpp index ab2048116..a89808470 100644 --- a/core/include/detray/propagator/composite_actor.hpp +++ b/core/include/detray/propagator/composite_actor.hpp @@ -1,13 +1,14 @@ /** Detray library, part of the ACTS project (R&D line) * - * (c) 2022-2025 CERN for the benefit of the ACTS project + * (c) 2022-2026 CERN for the benefit of the ACTS project * * Mozilla Public License Version 2.0 */ #pragma once -// Propagate include(s) +// Project include(s) +#include "detray/definitions/actor.hpp" #include "detray/definitions/containers.hpp" #include "detray/definitions/detail/qualifiers.hpp" #include "detray/propagator/base_actor.hpp" @@ -21,6 +22,7 @@ #include namespace detray { + /// Composition of actors /// /// The composition represents an actor together with its observers. In @@ -29,7 +31,7 @@ namespace detray { /// @tparam principal_actor_t the actor the compositions implements itself. /// @tparam observers a pack of observing actors that get called on the updated /// actor state of the compositions actor implementation. -template class composite_actor final : public principal_actor_t { @@ -40,6 +42,10 @@ class composite_actor final : public principal_actor_t { /// The composite is an actor in itself. using actor_type = principal_actor_t; using state = typename actor_type::state; + using result = typename actor_type::result; + + // Make sure the result type contains a status code + static_assert(std::is_base_of_v); /// Tuple of states of observing actors using observer_states = @@ -55,13 +61,13 @@ class composite_actor final : public principal_actor_t { /// /// @param states the states of all actors in the chain /// @param p_state the state of the propagator (stepper and navigator) - /// @param subject_state the state of the actor this actor observes. Uses + /// @param subject_res the result of the actor this actor observes. Uses /// a dummy type if this is not an observing actor. template - DETRAY_HOST_DEVICE void operator()( - actor_states_t &states, propagator_state_t &p_state, - subj_state_t &&subject_state = {}) const { + typename subj_result_t = typename actor::empty_result> + DETRAY_HOST_DEVICE void operator()(actor_states_t &states, + propagator_state_t &p_state, + subj_result_t &&subject_res = {}) const { // State of the primary actor that is implement by this composite actor auto &actor_state = detail::get(states); @@ -69,16 +75,19 @@ class composite_actor final : public principal_actor_t { // Do your own work ... // Two cases: This is a simple actor or observing actor (pass on its // subject's state) - if constexpr (std::same_as) { - actor_type::operator()(actor_state, p_state); + result res{}; + if constexpr (std::same_as) { + res = actor_type::operator()(actor_state, p_state); } else { - actor_type::operator()(actor_state, p_state, - std::forward(subject_state)); + res = actor_type::operator()( + actor_state, p_state, std::forward(subject_res)); } - // ... then run the observers on the updated state - notify(m_observers, states, actor_state, p_state, - std::make_index_sequence{}); + // ... then run the observers on the new result + if (res.status == actor::status::e_notify) { + notify(states, p_state, res, + std::make_index_sequence{}); + } } private: @@ -93,20 +102,20 @@ class composite_actor final : public principal_actor_t { typename propagator_state_t> DETRAY_HOST_DEVICE inline void notify(const observer_t &observer, actor_states_t &states, - state &actor_state, - propagator_state_t &p_state) const { + propagator_state_t &p_state, + result &res) const { // Two cases: observer is a simple actor or a composite actor if constexpr (!concepts::composite_actor) { // No actor state defined (empty) if constexpr (std::same_as) { - observer(actor_state, p_state); + detray::base_actor::state>) { + observer(p_state, res); } else { observer(detail::get(states), - actor_state, p_state); + p_state, res); } } else { - observer(states, actor_state, p_state); + observer(states, p_state, res); } } @@ -122,13 +131,10 @@ class composite_actor final : public principal_actor_t { template DETRAY_HOST_DEVICE inline void notify( - const dtuple &observer_list, actor_states_t &states, - state &actor_state, propagator_state_t &p_state, + actor_states_t &states, propagator_state_t &p_state, result &res, std::index_sequence /*ids*/) const { - (notify(detail::get(observer_list), states, actor_state, - p_state), - ...); + (notify(detail::get(m_observers), states, p_state, res), ...); } /// Keep the observers (might be composites again) diff --git a/core/include/detray/propagator/concepts.hpp b/core/include/detray/propagator/concepts.hpp index 7c5b618fd..5e1b09295 100644 --- a/core/include/detray/propagator/concepts.hpp +++ b/core/include/detray/propagator/concepts.hpp @@ -21,7 +21,7 @@ namespace detray::concepts { /// Concept for a simple actor template -concept actor = std::derived_from && +concept actor = std::derived_from && requires(const A a) { typename A::state; }; /// Concept for an actor including observing actors diff --git a/core/include/detray/propagator/detail/noise_estimation.hpp b/core/include/detray/propagator/detail/noise_estimation.hpp index 44d96f625..0ca618055 100644 --- a/core/include/detray/propagator/detail/noise_estimation.hpp +++ b/core/include/detray/propagator/detail/noise_estimation.hpp @@ -93,14 +93,12 @@ DETRAY_HOST_DEVICE constexpr void estimate_external_mask_tolerance( DETRAY_DEBUG_HOST_DEVICE("-> accumulated noise: %f", accumulated_noise); - navigation.set_external_tol(math::sqrt( + const scalar_t ext_tol{math::sqrt( n_stddev * n_stddev * (var_loc0 + var_loc1) + - vector::dot(displ, displ) + accumulated_noise * accumulated_noise)); + vector::dot(displ, displ) + accumulated_noise * accumulated_noise)}; // Clip to 5mm if the covariances are very large - navigation.set_external_tol(navigation.external_tol() > max_tol - ? max_tol - : navigation.external_tol()); + navigation.set_external_tol(ext_tol > max_tol ? max_tol : ext_tol); assert(std::isfinite(navigation.external_tol())); diff --git a/core/include/detray/propagator/detail/type_traits.hpp b/core/include/detray/propagator/detail/type_traits.hpp index 98c0759b4..34e12a9b4 100644 --- a/core/include/detray/propagator/detail/type_traits.hpp +++ b/core/include/detray/propagator/detail/type_traits.hpp @@ -1,6 +1,6 @@ /** Detray library, part of the ACTS project (R&D line) * - * (c) 2025 CERN for the benefit of the ACTS project + * (c) 2025-2026 CERN for the benefit of the ACTS project * * Mozilla Public License Version 2.0 */ @@ -26,10 +26,11 @@ struct get_state_tuple { using state_t = typename actor_t::state; // Remove empty default state of base actor type from tuple - using principal = std::conditional_t, - dtuple<>, dtuple>; + using principal = + std::conditional_t, dtuple<>, + dtuple>; using principal_ref = - std::conditional_t, dtuple<>, + std::conditional_t, dtuple<>, dtuple>; public: diff --git a/core/include/detray/propagator/propagator.hpp b/core/include/detray/propagator/propagator.hpp index 92a55e3f3..e76c6ddde 100644 --- a/core/include/detray/propagator/propagator.hpp +++ b/core/include/detray/propagator/propagator.hpp @@ -48,8 +48,8 @@ struct propagator { const propagation::config &m_cfg; - stepper_t m_stepper; - navigator_t m_navigator; + stepper_t m_stepper{}; + navigator_t m_navigator{}; /// Register the actor types const actor_chain_type run_actors{}; @@ -80,7 +80,7 @@ struct propagator { const free_track_parameters_type &free_params, const detector_type &det, const context_type &ctx) requires(!states_as_reference) - : _stepping(free_params), _navigation(det), _context(ctx) {} + : m_stepping(free_params), m_navigation(det), m_context(ctx) {} /// Construct the propagation state with free parameter template @@ -90,9 +90,9 @@ struct propagator { const free_track_parameters_type &free_params, const field_t &magnetic_field, const detector_type &det, const context_type &ctx = {}) - : _stepping(free_params, magnetic_field), - _navigation(det), - _context(ctx) {} + : m_stepping(free_params, magnetic_field), + m_navigation(det), + m_context(ctx) {} /// Construct the propagation state from the navigator state view DETRAY_HOST_DEVICE state_base( @@ -101,9 +101,9 @@ struct propagator { typename navigator_type::state::view_type nav_view, const context_type &ctx = {}) requires(!states_as_reference) - : _stepping(free_params), - _navigation(det, nav_view), - _context(ctx) {} + : m_stepping(free_params), + m_navigation(det, nav_view), + m_context(ctx) {} /// Construct the propagation state from the navigator state view template @@ -114,17 +114,17 @@ struct propagator { const field_t &magnetic_field, const detector_type &det, typename navigator_type::state::view_type nav_view, const context_type &ctx = {}) - : _stepping(free_params, magnetic_field), - _navigation(det, nav_view), - _context(ctx) {} + : m_stepping(free_params, magnetic_field), + m_navigation(det, nav_view), + m_context(ctx) {} /// Construct the propagation state with bound parameter DETRAY_HOST_DEVICE state_base(const bound_track_parameters_type ¶m, const detector_type &det, const context_type &ctx = {}) requires(!states_as_reference) - : _stepping(param, det, ctx), _navigation(det), _context(ctx) { - _navigation.set_volume(param.surface_link().volume()); + : m_stepping(param, det, ctx), m_navigation(det), m_context(ctx) { + m_navigation.set_volume(param.surface_link().volume()); } /// Construct the propagation state with propagation and navigation @@ -133,7 +133,7 @@ struct propagator { navigator_state_type &navigation, const context_type &ctx = {}) requires(states_as_reference) - : _stepping(stepping), _navigation(navigation), _context(ctx) {} + : m_stepping(stepping), m_navigation(navigation), m_context(ctx) {} /// Construct the propagation state with bound parameter template @@ -142,10 +142,10 @@ struct propagator { const field_t &magnetic_field, const detector_type &det, const context_type &ctx = {}) - : _stepping(param, magnetic_field, det, ctx), - _navigation(det), - _context(ctx) { - _navigation.set_volume(param.surface_link().volume()); + : m_stepping(param, magnetic_field, det, ctx), + m_navigation(det), + m_context(ctx) { + m_navigation.set_volume(param.surface_link().volume()); } /// Construct the propagation state with bound parameter and @@ -157,67 +157,67 @@ struct propagator { const field_t &magnetic_field, const detector_type &det, typename navigator_type::state::view_type nav_view, const context_type &ctx = {}) - : _stepping(param, magnetic_field, det, ctx), - _navigation(det, nav_view), - _context(ctx) { - _navigation.set_volume(param.surface_link().volume()); + : m_stepping(param, magnetic_field, det, ctx), + m_navigation(det, nav_view), + m_context(ctx) { + m_navigation.set_volume(param.surface_link().volume()); } /// Set the particle hypothesis DETRAY_HOST_DEVICE void set_particle(const pdg_particle &ptc) { - _stepping.set_particle(ptc); + m_stepping.set_particle(ptc); } /// @returns the propagation heartbeat DETRAY_HOST_DEVICE - bool is_alive() const { return _heartbeat; } + bool is_alive() const { return m_heartbeat; } /// @returns the propagation heartbeat DETRAY_HOST_DEVICE - bool heartbeat() const { return _heartbeat; } + bool heartbeat() const { return m_heartbeat; } /// @returns the propagation heartbeat DETRAY_HOST_DEVICE - void heartbeat(bool heartbeat) { _heartbeat = heartbeat; } + void heartbeat(bool heartbeat) { m_heartbeat = heartbeat; } DETRAY_HOST_DEVICE - const stepper_state_type &stepping() const { return _stepping; } + const stepper_state_type &stepping() const { return m_stepping; } DETRAY_HOST_DEVICE - stepper_state_type &stepping() { return _stepping; } + stepper_state_type &stepping() { return m_stepping; } DETRAY_HOST_DEVICE - const navigator_state_type &navigation() const { return _navigation; } + const navigator_state_type &navigation() const { return m_navigation; } DETRAY_HOST_DEVICE - navigator_state_type &navigation() { return _navigation; } + navigator_state_type &navigation() { return m_navigation; } DETRAY_HOST_DEVICE - const context_type &context() const { return _context; } + const context_type &context() const { return m_context; } DETRAY_HOST_DEVICE - context_type &context() { return _context; } + context_type &context() { return m_context; } DETRAY_HOST_DEVICE - bool debug() const { return do_debug; } + bool debug() const { return m_do_debug; } DETRAY_HOST_DEVICE - void debug(bool b) { do_debug = b; } + void debug(bool b) { m_do_debug = b; } private: - // Is the propagation still alive? - bool _heartbeat = false; - std::conditional_t - _stepping; + m_stepping; std::conditional_t - _navigation; - context_type _context; + m_navigation; + + context_type m_context; - bool do_debug = false; + // Is the propagation still alive? + bool m_heartbeat = false; + bool m_do_debug = false; }; using state = state_base; @@ -271,15 +271,6 @@ struct propagator { DETRAY_VERBOSE_HOST("Starting propagation for track:\n" << track); - // Open the navigation area according to uncertainties in initital - // track params - if (m_cfg.navigation.estimate_scattering_noise && - !stepping.bound_params().is_invalid()) { - detail::estimate_external_mask_tolerance( - stepping.bound_params(), propagation, - static_cast(m_cfg.navigation.n_scattering_stddev)); - } - // Initialize the navigation DETRAY_VERBOSE_HOST("Initialize navigation..."); m_navigator.init(track, navigation, m_cfg.navigation, context); @@ -419,7 +410,7 @@ struct propagator { std::stringstream debug_stream{}; debug_stream << std::left << std::setw(10); - debug_stream << "status: " << navigation.status() << std::endl; + debug_stream << "\nstatus: " << navigation.status() << std::endl; debug_stream << "volume: " << std::setw(10); if (detail::is_invalid_value(navigation.volume())) { @@ -432,18 +423,20 @@ struct propagator { debug_stream << "navigation:" << std::endl; if (navigation.is_on_surface()) { debug_stream << std::setw(10) - << " ->on surface: " << navigation.barcode(); + << " -> on surface: " << navigation.barcode(); } else { - debug_stream << std::setw(10) << " ->target: " + debug_stream << std::setw(10) << " -> target: " << navigation.target().surface().barcode(); } debug_stream << std::endl; - debug_stream << " ->path: " << navigation() << "mm" << std::endl; + debug_stream << std::setw(10) << " -> path: " << navigation() << " mm" + << std::endl; debug_stream << "stepping:" << std::endl; - debug_stream << " -> step size: " << std::setw(10) - << stepping.step_size() << "mm" << std::endl; - debug_stream << " ->" << detail::ray(stepping()); + debug_stream << std::setw(10) + << " -> step size: " << stepping.step_size() << " mm" + << std::endl; + debug_stream << " -> " << detail::ray(stepping()); return debug_stream.str(); } diff --git a/core/include/detray/tracks/bound_track_parameters.hpp b/core/include/detray/tracks/bound_track_parameters.hpp index ab82e3222..2e26965ac 100644 --- a/core/include/detray/tracks/bound_track_parameters.hpp +++ b/core/include/detray/tracks/bound_track_parameters.hpp @@ -378,7 +378,9 @@ struct bound_track_parameters : public bound_parameters_vector { return out_stream; } + /// Bound covaraicne matrix of the track parameters covariance_type m_covariance = matrix::zero(); + /// Identifier of the surface the track parameters are associated with geometry::barcode m_barcode{}; }; diff --git a/core/include/detray/utils/print_detector.hpp b/core/include/detray/utils/print_detector.hpp index 65bb2df75..2a5b5d2ec 100644 --- a/core/include/detray/utils/print_detector.hpp +++ b/core/include/detray/utils/print_detector.hpp @@ -49,7 +49,7 @@ struct material_printer { DETRAY_HOST void operator()(const material_coll_t &material_coll, const index_t idx, std::stringstream &os) const { - os << material_coll[idx]; + os << material_coll[idx] << std::endl; } }; @@ -95,7 +95,7 @@ DETRAY_HOST inline std::string print_detector( // Check the surface material, if present if (sf.has_material()) { - debug_stream << "[>>>>] Surface material:" << std::endl; + debug_stream << "[>>>>] Surface material: "; sf.template visit_material( debug_stream); } diff --git a/core/include/detray/utils/ranges/cartesian_product.hpp b/core/include/detray/utils/ranges/cartesian_product.hpp index 6c8ae3386..49eb9ed96 100644 --- a/core/include/detray/utils/ranges/cartesian_product.hpp +++ b/core/include/detray/utils/ranges/cartesian_product.hpp @@ -34,12 +34,11 @@ template struct cartesian_product_view : public detray::ranges::view_interface< cartesian_product_view> { - using iterator_coll_t = - detray::tuple...>; + using iterator_coll_t = dtuple...>; using iterator_t = detray::ranges::detail::cartesian_product_iterator< detray::ranges::iterator_t...>; - using value_type = detray::tuple< - std::iter_reference_t>...>; + using value_type = + dtuple>...>; /// Default constructor constexpr cartesian_product_view() = default; @@ -132,8 +131,8 @@ struct cartesian_product_iterator { /// Construct from a collection of @param begin and @param end positions DETRAY_HOST_DEVICE - constexpr cartesian_product_iterator(detray::tuple begins, - detray::tuple ends) + constexpr cartesian_product_iterator(dtuple begins, + dtuple ends) : m_begins(begins), m_ends(ends), m_itrs(begins) {} /// @returns true if the last range iterators are equal. @@ -225,10 +224,10 @@ struct cartesian_product_iterator { } /// Global range collection of begin and end iterators - detray::tuple m_begins; - detray::tuple m_ends; + dtuple m_begins; + dtuple m_ends; /// Current iterator states - detray::tuple m_itrs; + dtuple m_itrs; }; } // namespace detail diff --git a/core/include/detray/utils/tuple_helpers.hpp b/core/include/detray/utils/tuple_helpers.hpp index 38d08110e..f6ec486d6 100644 --- a/core/include/detray/utils/tuple_helpers.hpp +++ b/core/include/detray/utils/tuple_helpers.hpp @@ -29,25 +29,25 @@ using std::get; template DETRAY_HOST_DEVICE constexpr decltype(auto) get( - const ::detray::tuple& tuple) noexcept { + const ::detray::dtuple& tuple) noexcept { return ::detray::get(tuple); } template DETRAY_HOST_DEVICE constexpr decltype(auto) get( - ::detray::tuple& tuple) noexcept { + ::detray::dtuple& tuple) noexcept { return ::detray::get(tuple); } template DETRAY_HOST_DEVICE constexpr decltype(auto) get( - const ::detray::tuple& tuple) noexcept { + const ::detray::dtuple& tuple) noexcept { return ::detray::get>(tuple); } template DETRAY_HOST_DEVICE constexpr decltype(auto) get( - ::detray::tuple& tuple) noexcept { + ::detray::dtuple& tuple) noexcept { return ::detray::get>(tuple); } /// @} @@ -65,11 +65,11 @@ template struct tuple_element> : std::tuple_element> {}; -// detray::tuple +// detray::dtuple template -struct tuple_element> { +struct tuple_element> { using type = std::decay_t( - std::declval>()))>; + std::declval>()))>; }; template @@ -89,9 +89,9 @@ template struct tuple_size> : std::tuple_size> {}; -// detray::tuple +// detray::dtuple template -struct tuple_size<::detray::tuple> { +struct tuple_size<::detray::dtuple> { static constexpr std::size_t value = sizeof...(value_types); }; @@ -126,13 +126,13 @@ DETRAY_HOST constexpr std::tuple...> make_tuple( return std::make_tuple(std::forward(args)...); } -// make_tuple for detray::tuple +// make_tuple for detray::dtuple template