Skip to content

Commit

Permalink
Adapt pointwise interactor to handle material maps
Browse files Browse the repository at this point in the history
  • Loading branch information
niermann999 committed Sep 21, 2023
1 parent b56814c commit 946b8e8
Show file tree
Hide file tree
Showing 2 changed files with 65 additions and 43 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -39,7 +41,6 @@ struct pointwise_material_interactor : actor {
scalar_type mass{105.7f * unit<scalar_type>::MeV};
/// The particle pdg
int pdg = 13; // default muon

/// Evaluated energy loss
scalar_type e_loss{0.f};
/// Evaluated projected scattering angle
Expand Down Expand Up @@ -72,52 +73,72 @@ struct pointwise_material_interactor : actor {
const bound_track_parameters<transform3_type> &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 <class material_coll_t, class point_t,
std::enable_if_t<not detail::is_grid_v<
typename material_coll_t::value_type>,
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 <class material_coll_t, class point_t,
std::enable_if_t<
detail::is_grid_v<typename material_coll_t::value_type>,
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];
}
};

Expand Down
11 changes: 6 additions & 5 deletions core/include/detray/surface_finders/grid/grid.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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 <typename point_t>
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 <typename point_t>
DETRAY_HOST_DEVICE auto search(const point_t &p) {
return bin(m_axes.bins(p));
}

Expand Down

0 comments on commit 946b8e8

Please sign in to comment.