test/validation/material_validation_utils.hpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 // Project include(s)
0012 #include "detray/definitions/algebra.hpp"
0013 #include "detray/material/concepts.hpp"
0014 #include "detray/material/detail/material_accessor.hpp"
0015 #include "detray/navigation/caching_navigator.hpp"
0016 #include "detray/navigation/detail/print_state.hpp"
0017 #include "detray/propagator/actors.hpp"
0018 #include "detray/propagator/line_stepper.hpp"
0019 #include "detray/propagator/propagator.hpp"
0020 #include "detray/utils/logging.hpp"
0021 #include "detray/utils/type_list.hpp"
0022
0023 // Detray IO include(s)
0024 #include "detray/io/utils/create_path.hpp"
0025 #include "detray/io/utils/file_handle.hpp"
0026
0027 // Vecmem include(s)
0028 #include <vecmem/memory/memory_resource.hpp>
0029
0030 // System include(s)
0031 #include <filesystem>
0032
0033 namespace detray::material_validator {
0034
0035 /// @brief Record the material budget per thickness or pathlength
0036 template <concepts::scalar scalar_t>
0037 struct material_record {
0038   /// Phi and eta values of the track for which the material was recorded
0039   /// @{
0040   scalar_t phi{detail::invalid_value<scalar_t>()};
0041   scalar_t eta{detail::invalid_value<scalar_t>()};
0042   /// @}
0043   /// Accumulated radiation length per pathlength through the material
0044   scalar_t sX0{0.f};
0045   /// Accumulated radiation length per thickness
0046   scalar_t tX0{0.f};
0047   /// Accumulated interaction length per pathlength through the material
0048   scalar_t sL0{0.f};
0049   /// Accumulated interaction length per thickness
0050   scalar_t tL0{0.f};
0051 };
0052
0053 /// @brief Return type that contains the material parameters and the pathlength
0054 template <concepts::scalar scalar_t>
0055 struct material_params {
0056   /// The surface the material belongs to
0057   geometry::identifier geo_id{};
0058   /// Pathlength of the track through the material
0059   scalar_t path{detail::invalid_value<scalar_t>()};
0060   /// Material thickness/radius
0061   scalar_t thickness{detail::invalid_value<scalar_t>()};
0062   /// Radiation length
0063   scalar_t mat_X0{0.f};
0064   /// Interaction length
0065   scalar_t mat_L0{0.f};
0066 };
0067
0068 /// @brief Functor to retrieve the material parameters for a given local
0069 /// position
0070 struct get_material_params {
0071   template <typename mat_group_t, typename index_t, concepts::point2D point2_t,
0072             concepts::scalar scalar_t>
0073   DETRAY_HOST_DEVICE auto operator()(
0074       [[maybe_unused]] const mat_group_t &mat_group,
0075       [[maybe_unused]] const index_t &index,
0076       [[maybe_unused]] const point2_t &loc,
0077       [[maybe_unused]] const scalar_t cos_inc_angle) const {
0078     using material_t = typename mat_group_t::value_type;
0079
0080     constexpr auto inv{detail::invalid_value<scalar_t>()};
0081     constexpr geometry::identifier inv_sf{};
0082
0083     // Access homogeneous surface material or material maps
0084     if constexpr (concepts::surface_material<material_t>) {
0085       // Slab or rod
0086       const auto mat = detail::material_accessor::get(mat_group, index, loc);
0087
0088       // Empty material can occur in material maps, skip it
0089       if (!mat) {
0090         // Set the pathlength and thickness to zero so that they
0091         // are not counted
0092         return material_params<scalar_t>{inv_sf, 0.f, 0.f, inv, inv};
0093       }
0094
0095       const scalar_t seg{mat.path_segment(cos_inc_angle, loc[0])};
0096       const scalar_t t{mat.thickness()};
0097       const scalar_t mat_X0{mat.get_material().X0()};
0098       const scalar_t mat_L0{mat.get_material().L0()};
0099
0100       return material_params<scalar_t>{inv_sf, seg, t, mat_X0, mat_L0};
0101     } else {
0102       return material_params<scalar_t>{inv_sf, inv, inv, inv, inv};
0103     }
0104   }
0105 };
0106
0107 /// @brief Actor that collects all material encountered by a track during
0108 ///        navigation
0109 ///
0110 /// The material is scaled with either the slab thickness or pathlength through
0111 /// the material.
0112 template <concepts::scalar scalar_t, template <typename...> class vector_t>
0113 struct material_tracer : public detray::base_actor {
0114   using material_record_type = material_record<scalar_t>;
0115   using material_params_type = material_params<scalar_t>;
0116
0117   struct state {
0118     friend struct material_tracer;
0119
0120     /// Construct the vector containers with a given resource
0121     /// @param resource
0122     DETRAY_HOST
0123     explicit state(vecmem::memory_resource &resource)
0124         : m_mat_steps(&resource) {}
0125
0126     /// Construct from externally provided vector for the @param steps
0127     DETRAY_HOST_DEVICE
0128     explicit state(vector_t<material_params<scalar_t>> &&steps)
0129         : m_mat_steps(std::move(steps)) {}
0130
0131     /// Access to the total recorded material along the track - const
0132     DETRAY_HOST_DEVICE
0133     const auto &get_material_record() const { return m_mat_record; }
0134
0135     /// Move the total recorded material out of the actor
0136     DETRAY_HOST
0137     auto &&release_material_record() && { return std::move(m_mat_record); }
0138
0139     /// Access to the recorded material steps along the track - const
0140     DETRAY_HOST_DEVICE
0141     const auto &get_material_steps() const { return m_mat_steps; }
0142
0143     /// Move the recorded material steps pout of the actor
0144     DETRAY_HOST_DEVICE
0145     auto &&release_material_steps() && { return std::move(m_mat_steps); }
0146
0147    private:
0148     /// Accumulated material data for the track
0149     material_record<scalar_t> m_mat_record{};
0150
0151     /// Collect material parameters for every step
0152     vector_t<material_params<scalar_t>> m_mat_steps{};
0153   };
0154
0155   /// Run as observer to the parameter transporter (covariance transport)
0156   template <typename propagator_state_t, concepts::algebra algebra_t>
0157   DETRAY_HOST_DEVICE void operator()(
0158       state &tracer, const propagator_state_t &prop_state,
0159       const detray::actor::parameter_transporter_result<algebra_t> &res) const {
0160     record_track_dir(tracer, prop_state);
0161
0162     // Only count material if navigator encountered it
0163     const auto &navigation = prop_state.navigation();
0164     if (!navigation.encountered_sf_material()) {
0165       return;
0166     }
0167
0168     // For now use default context
0169     typename propagator_state_t::detector_type::geometry_context gctx{};
0170
0171     const auto &track_param = res.destination_params();
0172     dpoint2D<algebra_t> loc_pos = track_param.bound_local();
0173
0174     record_mat_step(tracer, gctx, navigation.current_surface(), loc_pos,
0175                     track_param.dir());
0176   }
0177
0178   /// Run in a propagation chain without covariance transport
0179   template <typename propagator_state_t>
0180   DETRAY_HOST_DEVICE void operator()(
0181       state &tracer, const propagator_state_t &prop_state) const {
0182     using algebra_t = typename propagator_state_t::detector_type::algebra_type;
0183
0184     record_track_dir(tracer, prop_state);
0185
0186     // Only count material if navigator encountered it
0187     const auto &navigation = prop_state.navigation();
0188     if (!navigation.encountered_sf_material()) {
0189       return;
0190     }
0191
0192     // For now use default context
0193     typename propagator_state_t::detector_type::geometry_context gctx{};
0194
0195     // Current surface
0196     const auto sf = navigation.current_surface();
0197
0198     const auto &track_param = prop_state.stepping()();
0199     dvector3D<algebra_t> glob_dir = track_param.dir();
0200     dpoint2D<algebra_t> loc_pos =
0201         sf.global_to_bound(gctx, track_param.pos(), glob_dir);
0202
0203     record_mat_step(tracer, gctx, sf, loc_pos, glob_dir);
0204   }
0205
0206  private:
0207   /// Record the track direction
0208   template <typename propagator_state_t>
0209   DETRAY_HOST_DEVICE inline auto record_track_dir(
0210       state &tracer, const propagator_state_t &prop_state) const {
0211     using algebra_t = typename propagator_state_t::detector_type::algebra_type;
0212     using vector3_t = dvector3D<algebra_t>;
0213
0214     // Record the initial track direction
0215     vector3_t glob_dir = prop_state.stepping()().dir();
0216     if (detray::detail::is_invalid_value(tracer.m_mat_record.eta) &&
0217         detray::detail::is_invalid_value(tracer.m_mat_record.phi)) {
0218       tracer.m_mat_record.eta = vector::eta(glob_dir);
0219       tracer.m_mat_record.phi = vector::phi(glob_dir);
0220     }
0221   }
0222
0223   /// Record the data for a material step
0224   template <typename detector_t>
0225   DETRAY_HOST_DEVICE inline auto record_mat_step(
0226       state &tracer, const typename detector_t::geometry_context &gctx,
0227       const tracking_surface<detector_t> sf,
0228       const dpoint2D<typename detector_t::algebra_type> &loc_pos,
0229       const dvector3D<typename detector_t::algebra_type> &glob_dir) const {
0230     // Fetch the material parameters and pathlength through the material
0231     const auto mat_params = sf.template visit_material<get_material_params>(
0232         loc_pos, cos_angle(gctx, sf, glob_dir, loc_pos));
0233
0234     const scalar_t seg{mat_params.path};
0235     const scalar_t t{mat_params.thickness};
0236     const scalar_t mx0{mat_params.mat_X0};
0237     const scalar_t ml0{mat_params.mat_L0};
0238
0239     // Fill the material record
0240     if (mx0 > 0.f) {
0241       tracer.m_mat_record.sX0 += seg / mx0;
0242       tracer.m_mat_record.tX0 += t / mx0;
0243     }
0244     if (ml0 > 0.f) {
0245       tracer.m_mat_record.sL0 += seg / ml0;
0246       tracer.m_mat_record.tL0 += t / ml0;
0247     }
0248     if (t > 0.f) {
0249       tracer.m_mat_steps.push_back({sf.identifier(), seg, t, mx0, ml0});
0250     }
0251   }
0252 };
0253
0254 /// Run the propagation and record test data along the way
0255 template <typename detector_t>
0256 inline auto record_material(
0257     const typename detector_t::geometry_context /*gctx*/,
0258     vecmem::memory_resource *host_mr, const detector_t &det,
0259     const propagation::config &cfg,
0260     const free_track_parameters<typename detector_t::algebra_type> &track) {
0261   using algebra_t = typename detector_t::algebra_type;
0262   using scalar_t = dscalar<algebra_t>;
0263
0264   using stepper_t = line_stepper<algebra_t>;
0265   using navigator_t = caching_navigator<detector_t>;
0266
0267   // Propagator with pathlimit aborter
0268   using material_tracer_t =
0269       material_validator::material_tracer<scalar_t, vecmem::vector>;
0270   using pathlimit_aborter_t = actor::pathlimit_aborter<scalar_t>;
0271   using parameter_updater_t =
0272       actor::parameter_updater<algebra_t,
0273                                actor::pointwise_material_interactor<algebra_t>,
0274                                material_tracer_t>;
0275   using actor_chain_t = actor_chain<pathlimit_aborter_t, parameter_updater_t>;
0276   using propagator_t = propagator<stepper_t, navigator_t, actor_chain_t>;
0277
0278   // Propagator
0279   propagator_t prop{cfg};
0280
0281   // Build actor and propagator states
0282   typename pathlimit_aborter_t::state pathlimit_aborter_state{
0283       cfg.stepping.path_limit};
0284   actor::parameter_updater_state<algebra_t> updater_state{cfg};
0285   typename actor::pointwise_material_interactor<algebra_t>::state
0286       interactor_state{};
0287   typename material_tracer_t::state mat_tracer_state{*host_mr};
0288
0289   auto actor_states = detray::tie(pathlimit_aborter_state, updater_state,
0290                                   interactor_state, mat_tracer_state);
0291
0292   typename propagator_t::state propagation{track, det, cfg.context};
0293
0294   // Run the propagation
0295   bool success = prop.propagate(propagation, actor_states);
0296
0297   return std::make_tuple(success,
0298                          std::move(mat_tracer_state).release_material_record(),
0299                          std::move(mat_tracer_state).release_material_steps());
0300 }
0301
0302 /// Compare the result of two material tracers
0303 ///
0304 /// @param reference the reference trace to compare against
0305 /// @param ref_record the reference total material to compare against (s/X0)
0306 /// @param mat_trace the recoded material trace
0307 /// @param mat_record the recoded total material along the track (s/X0)
0308 /// @param trk_i the track index/identifier
0309 /// @param rel_tol the maximum allowed relative error for any parameters
0310 /// @param verbose debug output
0311 ///
0312 /// @returns if the traces contains matching material steps and if matching
0313 /// traces contain the same total material
0314 template <concepts::scalar scalar_t>
0315 inline auto compare_traces(
0316     const dvector<material_params<scalar_t>> &reference,
0317     const material_record<scalar_t> &ref_record,
0318     const dvector<material_params<scalar_t>> &mat_trace,
0319     const material_record<scalar_t> &mat_record,
0320     std::size_t trk_i = detail::invalid_value<std::size_t>(),
0321     const double rel_tol = 0.01, const bool verbose = true) {
0322   // Material traces contain records of different surfaces/material
0323   bool is_bad_comp{reference.size() != mat_trace.size()};
0324   // Overall recorded material is different (in case small errors accumulate)
0325   bool is_diff_mat{false};
0326
0327   std::stringstream debug_msg{};
0328   debug_msg << "Track No. " << trk_i << ":\n----------------" << std::endl;
0329
0330   if (is_bad_comp) {
0331     debug_msg << "-> Different no. of surfaces: " << mat_trace.size()
0332               << " (ref.: " << reference.size() << ")\n"
0333               << std::endl;
0334   } else {
0335     for (std::size_t j = 0u; j < math::min(reference.size(), mat_trace.size());
0336          ++j) {
0337       if (reference[j].geo_id != mat_trace[j].geo_id) {
0338         is_bad_comp = true;
0339         debug_msg << "-> Surfaces don't match: " << mat_trace[j].geo_id
0340                   << " (ref.: " << reference[j].geo_id << ")" << std::endl;
0341         continue;
0342       }
0343
0344       // TODO: Use approx_equal from algebra-plugins
0345       // Compare thickness of the surface material
0346       if ((reference[j].thickness - mat_trace[j].thickness) /
0347               reference[j].thickness >
0348           rel_tol) {
0349         is_bad_comp = true;
0350         debug_msg << "-> On surface: " << reference[j].geo_id << ":"
0351                   << std::endl;
0352         debug_msg << "-> thickness: " << mat_trace[j].thickness
0353                   << ", thickness ref.: " << reference[j].thickness
0354                   << std::endl;
0355       }
0356
0357       // Compare radiation length of the surface material
0358       if ((reference[j].mat_X0 - mat_trace[j].mat_X0) / reference[j].mat_X0 >
0359           rel_tol) {
0360         is_bad_comp = true;
0361         debug_msg << "-> On surface: " << reference[j].geo_id << ":"
0362                   << std::endl;
0363         debug_msg << "-> X0: " << mat_trace[j].mat_X0
0364                   << ", X0 ref.: " << reference[j].mat_X0 << std::endl;
0365       }
0366
0367       // Compare interaction length of the surface material
0368       if ((reference[j].mat_L0 - mat_trace[j].mat_L0) / reference[j].mat_L0 >
0369           rel_tol) {
0370         is_bad_comp = true;
0371         debug_msg << "-> On surface: " << reference[j].geo_id << ":"
0372                   << std::endl;
0373         debug_msg << "-> L0: " << mat_trace[j].mat_L0
0374                   << ", L0 ref.: " << reference[j].mat_L0 << std::endl;
0375       }
0376
0377       // Compare path of the track through the surface material
0378       if ((reference[j].path - mat_trace[j].path) / reference[j].path >
0379           rel_tol) {
0380         is_bad_comp = true;
0381         debug_msg << "-> On surface: " << reference[j].geo_id << ":"
0382                   << std::endl;
0383         debug_msg << "-> Mat. path: " << mat_trace[j].path / unit<scalar_t>::mm
0384                   << " mm, mat. path ref.: "
0385                   << reference[j].path / unit<scalar_t>::mm << " mm"
0386                   << std::endl;
0387       }
0388     }
0389   }
0390
0391   // Compare the total accumulated material along the track
0392   const double ref_mat_X0{static_cast<double>(ref_record.sX0)};
0393   const double mat_X0{static_cast<double>(mat_record.sX0)};
0394   const double rel_error{(ref_mat_X0 - mat_X0) / ref_mat_X0};
0395   const bool small_mat{ref_mat_X0 < rel_tol && mat_X0 < rel_tol};
0396
0397   // If almost no material was collected, the relative error can be
0398   // large, but also has little consequence for tracking
0399   if (!(small_mat || (rel_error <= rel_tol))) {
0400     // Already know that material cannot add up: Only flag the additional
0401     // cases here
0402     if (!is_bad_comp) {
0403       is_diff_mat = true;
0404     }
0405     debug_msg << "\nTotal material discrepancy of " << 100. * rel_error << "%\n"
0406               << std::endl;
0407   }
0408
0409   if (verbose && (is_bad_comp || is_diff_mat)) {
0410     DETRAY_INFO_HOST(debug_msg.str());
0411   }
0412
0413   return std::make_tuple(is_bad_comp, is_diff_mat);
0414 }
0415
0416 /// Write the accumulated material of a track from @param mat_records to a csv
0417 /// file to the path @param mat_file_name
0418 template <concepts::scalar scalar_t>
0419 auto write_material(const std::string &mat_file_name,
0420                     const dvector<material_record<scalar_t>> &mat_records) {
0421   const auto file_path = std::filesystem::path{mat_file_name};
0422   assert(file_path.extension() == ".csv");
0423
0424   // Make sure path to file exists
0425   io::create_path(file_path.parent_path());
0426
0427   detray::io::file_handle outfile{
0428       mat_file_name, std::ios::out | std::ios::binary | std::ios::trunc};
0429   *outfile << "eta,phi,mat_sX0,mat_sL0,mat_tX0,mat_tL0" << std::endl;
0430
0431   for (const auto &rec : mat_records) {
0432     *outfile << rec.eta << "," << rec.phi << "," << rec.sX0 << "," << rec.sL0
0433              << "," << rec.tX0 << "," << rec.tL0 << std::endl;
0434   }
0435 }
0436
0437 }  // namespace detray::material_validator