From 946b8e8437f965fe95494b491989a5881719cbca Mon Sep 17 00:00:00 2001 From: Joana Niermann Date: Thu, 21 Sep 2023 16:30:54 +0200 Subject: [PATCH] Adapt pointwise interactor to handle material maps --- .../actors/pointwise_material_interactor.hpp | 97 +++++++++++-------- .../detray/surface_finders/grid/grid.hpp | 11 ++- 2 files changed, 65 insertions(+), 43 deletions(-) diff --git a/core/include/detray/propagator/actors/pointwise_material_interactor.hpp b/core/include/detray/propagator/actors/pointwise_material_interactor.hpp index 1e1e150ca1..ef9711908a 100644 --- a/core/include/detray/propagator/actors/pointwise_material_interactor.hpp +++ b/core/include/detray/propagator/actors/pointwise_material_interactor.hpp @@ -12,6 +12,8 @@ #include "detray/definitions/track_parametrization.hpp" #include "detray/geometry/surface.hpp" #include "detray/materials/interaction.hpp" +#include "detray/materials/material_rod.hpp" +#include "detray/materials/material_slab.hpp" #include "detray/propagator/base_actor.hpp" #include "detray/tracks/bound_track_parameters.hpp" #include "detray/utils/ranges.hpp" @@ -39,7 +41,6 @@ struct pointwise_material_interactor : actor { scalar_type mass{105.7f * unit::MeV}; /// The particle pdg int pdg = 13; // default muon - /// Evaluated energy loss scalar_type e_loss{0.f}; /// Evaluated projected scattering angle @@ -72,52 +73,72 @@ struct pointwise_material_interactor : actor { const bound_track_parameters &bound_params, const scalar_type cos_inc_angle, const scalar_type approach) const { - const scalar qop = bound_params.qop(); - const scalar charge = bound_params.charge(); + const auto &mat = this->get_material(material_group, mat_index, + bound_params.bound_local()); - bool success = false; - - if constexpr (detail::is_grid_v< - typename material_group_t::value_type>) { + // return early in case of vacuum or zero thickness + if (not mat) { return false; - } else { - const auto &mat = material_group[mat_index]; + } - // return early in case of vacuum or zero thickness - if (not mat) { - return false; - } - success = true; + const scalar qop = bound_params.qop(); + const scalar charge = bound_params.charge(); - const scalar_type path_segment{ - mat.path_segment(cos_inc_angle, approach)}; + const scalar_type path_segment{ + mat.path_segment(cos_inc_angle, approach)}; + + // Energy Loss + if (s.do_energy_loss) { + s.e_loss = interaction_type().compute_energy_loss_bethe( + path_segment, mat, s.pdg, s.mass, qop, charge); + } - // Energy Loss - if (s.do_energy_loss) { - s.e_loss = interaction_type().compute_energy_loss_bethe( + // @todo: include the radiative loss (Bremsstrahlung) + if (s.do_energy_loss && s.do_covariance_transport) { + s.sigma_qop = + interaction_type().compute_energy_loss_landau_sigma_QOverP( path_segment, mat, s.pdg, s.mass, qop, charge); - } + } - // @todo: include the radiative loss (Bremsstrahlung) - if (s.do_energy_loss && s.do_covariance_transport) { - s.sigma_qop = - interaction_type() - .compute_energy_loss_landau_sigma_QOverP( - path_segment, mat, s.pdg, s.mass, qop, charge); - } + // Covariance update + if (s.do_multiple_scattering) { + // @todo: use momentum before or after energy loss in + // backward mode? + s.projected_scattering_angle = + interaction_type().compute_multiple_scattering_theta0( + mat.path_segment_in_X0(cos_inc_angle, approach), s.pdg, + s.mass, qop, charge); + } - // Covariance update - if (s.do_multiple_scattering) { - // @todo: use momentum before or after energy loss in - // backward mode? - s.projected_scattering_angle = - interaction_type().compute_multiple_scattering_theta0( - mat.path_segment_in_X0(cos_inc_angle, approach), - s.pdg, s.mass, qop, charge); - } + return true; + } - return success; - } + private: + /// Access to material slabs or rods in a homogeneous material + /// description + template , + bool> = true> + inline constexpr auto &get_material( + const material_coll_t &material_coll, const dindex idx, + const point_t &) const noexcept { + return material_coll[idx]; + } + + /// Access to material slabs in a material map + template , + bool> = true> + inline constexpr auto &get_material( + const material_coll_t &material_coll, const dindex idx, + const point_t &loc_point) const noexcept { + // Get the material grid instance + const auto &mat_grid = material_coll[idx]; + + // Find the material slab (only one entry in bin) + return (mat_grid.search(loc_point))[0]; } }; diff --git a/core/include/detray/surface_finders/grid/grid.hpp b/core/include/detray/surface_finders/grid/grid.hpp index 17939a75d8..b3643c95d2 100644 --- a/core/include/detray/surface_finders/grid/grid.hpp +++ b/core/include/detray/surface_finders/grid/grid.hpp @@ -284,20 +284,21 @@ class grid { /// Find the value of a single bin - const /// - /// @param p is point in the local frame + /// @param p is point in the local (bound) frame /// /// @return the iterable view of the bin content - DETRAY_HOST_DEVICE auto search( - const typename local_frame_type::point3 &p) const { + template + DETRAY_HOST_DEVICE auto search(const point_t &p) const { return bin(m_axes.bins(p)); } /// Find the value of a single bin /// - /// @param p is point in the local frame + /// @param p is point in the local (bound) frame /// /// @return the iterable view of the bin content - DETRAY_HOST_DEVICE auto search(const typename local_frame_type::point3 &p) { + template + DETRAY_HOST_DEVICE auto search(const point_t &p) { return bin(m_axes.bins(p)); }