test/validation/propagation_validation.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/navigation/intersection/intersection.hpp"
0013 #include "detray/propagator/rk_stepper.hpp"
0014 #include "detray/utils/logging.hpp"
0015
0016 // Detray test include(s)
0017 #include "detray/test/utils/perigee_stopper.hpp"
0018 #include "detray/test/validation/material_validation_utils.hpp"
0019 #include "detray/test/validation/navigation_validation_utils.hpp"
0020
0021 // VecMem include(s).
0022 #include <vecmem/memory/host_memory_resource.hpp>
0023
0024 // System include(s).
0025 #include <iostream>
0026 #include <memory>
0027 #include <optional>
0028 #include <string>
0029
0030 namespace detray {
0031
0032 template <typename detector_t>
0033 using intersection_type =
0034     intersection2D<typename detector_t::surface_type,
0035                    typename detector_t::algebra_type, true>;
0036
0037 template <typename detector_t>
0038 using candidate_type =
0039     navigation::detail::candidate_record<intersection_type<detector_t>>;
0040
0041 /// Configure the Kalman Filter comparison test
0042 struct propagation_validation_config {
0043   propagation::config propagation;
0044
0045   pdg_particle<float> particle = detray::muon<float>();
0046
0047   bool do_multiple_scattering{true};
0048   bool do_energy_loss{true};
0049
0050   bool display_svg{true};
0051   float max_percent_missed{1.f};
0052   float max_percent_additional{15.f};
0053 };
0054
0055 /// Run a detray propagation and compare the surfaces that were encountered
0056 /// against a reference of truth tracks givben as intersection traces
0057 ///
0058 /// @param det the detector
0059 /// @param names the detector and volume names
0060 /// @param bfield the magnetic field representation
0061 /// @param cfg the full test configuration
0062 /// @param tracks the initial particle vertices as free track parameters
0063 /// @param truth_traces_fw the reference data for each encountered module
0064 ///
0065 /// @returns whether the validation was successful
0066 template <typename detector_t, typename bfield_view_t>
0067 bool propagation_validation(
0068     const detector_t& det, const typename detector_t::name_map& names,
0069     const std::optional<bfield_view_t>& bfield,
0070     const propagation_validation_config& cfg,
0071     const std::vector<free_track_parameters<typename detector_t::algebra_type>>&
0072         tracks,
0073     std::vector<vecmem::vector<candidate_type<detector_t>>>& truth_traces_fw) {
0074   // 'false' if any failures were detected
0075   bool test_successful{true};
0076
0077   using algebra_t = typename detector_t::algebra_type;
0078   using scalar_t = dscalar<algebra_t>;
0079   using stepper_t =
0080       rk_stepper<bfield_view_t, algebra_t, constrained_step<scalar_t>,
0081                  stepper_rk_policy<scalar_t>, stepping::print_inspector>;
0082   using sf_candidate_t =
0083       traccc::propagation_validator::candidate_type<detector_t>;
0084
0085   // Host memory resource
0086   vecmem::host_memory_resource host_mr;
0087
0088   // Geometry context
0089   const typename detector_t::geometry_context ctx{};
0090
0091   // Configure the test
0092   test::navigation_validation_config<algebra_t> test_cfg{};
0093   test_cfg.n_tracks(tracks.size());
0094   test_cfg.ptc_hypothesis(cfg.particle);
0095   test_cfg.collect_sensitives_only(true);
0096   test_cfg.fail_on_diff(false);
0097   test_cfg.display_svg(cfg.display_svg);
0098   test_cfg.display_only_missed(true);
0099   test_cfg.verbose(false);
0100
0101   // Create the backward truth traces by reverting the order
0102   std::vector<vecmem::vector<sf_candidate_t>> truth_traces_bw;
0103   truth_traces_bw.reserve(tracks.capacity());
0104
0105   for (const auto& truth_trace_fw : truth_traces_fw) {
0106     // Revert the forward trace for the backward propagation
0107     vecmem::vector<sf_candidate_t> truth_trace_bw(truth_trace_fw.size());
0108     std::ranges::reverse_copy(truth_trace_fw, truth_trace_bw.begin());
0109
0110     assert(!truth_trace_bw.empty());
0111
0112     truth_traces_bw.push_back(std::move(truth_trace_bw));
0113   }
0114
0115   // Check truth data
0116   if (truth_traces_fw.empty() || truth_traces_bw.empty()) {
0117     DETRAY_ERROR_HOST("Propagation truth data empty");
0118     return false;
0119   }
0120
0121   // Make a tuple of references from a tuple
0122   auto setup_actor_states = []<typename... T>(detray::dtuple<T...>& t) {
0123     return detray::tie(detray::detail::get<T>(t)...);
0124   };
0125
0126   // Define the actors
0127   using interactor = actor::pointwise_material_interactor<algebra_t>;
0128   using perigee_stopper = perigee_stopper<algebra_t>;
0129
0130   // Reusable actor states
0131   typename interactor::state interactor_state{};
0132   interactor_state.do_multiple_scattering = cfg.do_multiple_scattering;
0133   interactor_state.do_energy_loss = cfg.do_energy_loss;
0134
0135   typename perigee_stopper::state stopper_state{};
0136
0137   std::cout << "-----------------------------------"
0138             << "\nFORWARD - No KF" << std::endl
0139             << "-----------------------------------\n";
0140
0141   using parameter_updater = actor::parameter_updater<algebra_t, interactor>;
0142
0143   using actor_chain_t = actor_chain<parameter_updater>;
0144
0145   typename parameter_updater::state updater_state{};
0146   updater_state.noise_estimation_cfg().n_stddev =
0147       cfg.propagation.navigation.n_scattering_stddev;
0148   updater_state.noise_estimation_cfg().accumulated_error =
0149       cfg.propagation.navigation.accumulated_error;
0150   updater_state.noise_estimation_cfg().estimate_scattering_noise =
0151       cfg.propagation.navigation.estimate_scattering_noise;
0152
0153   // Prepare the actor states for every track
0154   vecmem::vector<typename actor_chain_t::state_tuple> state_tuples{};
0155   vecmem::vector<typename actor_chain_t::state_ref_tuple> state_ref_tuples{};
0156   state_tuples.reserve(tracks.size());
0157   state_ref_tuples.reserve(tracks.size());
0158
0159   for (std::size_t i = 0u; i < tracks.size(); ++i) {
0160     // Define the initial covariance
0161     updater_state.init(tracks.at(i));
0162
0163     // Copy the configured updater state for this track
0164     state_tuples.push_back(detray::make_tuple(updater_state, interactor_state));
0165     state_ref_tuples.push_back(setup_actor_states(state_tuples.back()));
0166   }
0167
0168   // Forward navigation
0169   test_cfg.name(det.name(names) + "_GeV_fw");
0170   test_cfg.navigation_direction(navigation::direction::e_forward);
0171   const auto [trk_stats_fw, n_surfaces_fw, n_miss_nav_fw, n_miss_truth_fw,
0172               step_traces_fw, mat_traces_fw, mat_records_fw] =
0173       navigation_validator::compare_to_navigation<stepper_t, parameter_updater>(
0174           test_cfg, host_mr, det, names, ctx, bfield, cfg.propagation,
0175           truth_traces_fw, tracks, state_ref_tuples);
0176
0177   // Check stats
0178   auto n_tracks{static_cast<double>(trk_stats_fw.n_tracks)};
0179   if (static_cast<double>(trk_stats_fw.n_tracks_w_holes) / n_tracks >
0180       cfg.max_percent_missed / 100.) {
0181     DETRAY_ERROR_HOST("Too many tracks with missing surfaces");
0182     test_successful = false;
0183   }
0184   if (static_cast<double>(trk_stats_fw.n_tracks_w_extra) / n_tracks >
0185       cfg.max_percent_additional / 100.) {
0186     DETRAY_ERROR_HOST("Too many tracks with additional surfaces");
0187     test_successful = false;
0188   }
0189
0190   std::cout << "\n-----------------------------------\n"
0191             << "BACKWARD - No KF" << std::endl
0192             << "-----------------------------------\n";
0193
0194   using actor_chain_bw_t = actor_chain<parameter_updater, perigee_stopper>;
0195
0196   // Prepare the actor states for every track
0197   vecmem::vector<typename actor_chain_bw_t::state_tuple> bw_state_tuples{};
0198   vecmem::vector<typename actor_chain_bw_t::state_ref_tuple>
0199       bw_state_ref_tuples{};
0200   bw_state_tuples.reserve(tracks.size());
0201   bw_state_ref_tuples.reserve(tracks.size());
0202
0203   for (std::size_t i = 0u; i < tracks.size(); ++i) {
0204     // Define the initial covariance
0205     updater_state.init(tracks.at(i));
0206
0207     // Copy the configured updater state for this track
0208     bw_state_tuples.push_back(
0209         detray::make_tuple(stopper_state, interactor_state, updater_state));
0210     bw_state_ref_tuples.push_back(setup_actor_states(bw_state_tuples.back()));
0211   }
0212
0213   // Backward navigation
0214   test_cfg.name(det.name(names) + "_GeV_bw");
0215   test_cfg.navigation_direction(navigation::direction::e_backward);
0216   const auto [trk_stats_bw, n_surfaces_bw, n_miss_nav_bw, n_miss_truth_bw,
0217               step_traces_bw, mat_traces_bw, mat_records_bw] =
0218       navigation_validator::compare_to_navigation<stepper_t, parameter_updater,
0219                                                   perigee_stopper>(
0220           test_cfg, host_mr, det, names, ctx, bfield, cfg.propagation,
0221           truth_traces_bw, tracks, bw_state_ref_tuples);
0222
0223   // Make sure some data was collected
0224   assert(trk_stats_fw.n_tracks > 0u);
0225   assert(n_surfaces_fw.n_total() > 0u);
0226   assert(trk_stats_bw.n_tracks > 0u);
0227   assert(n_surfaces_bw.n_total() > 0u);
0228
0229   assert(trk_stats_fw.n_tracks == trk_stats_bw.n_tracks);
0230
0231   // Check, the amount of collected material between forward and backward
0232   assert(step_traces_fw.size() == trk_stats_fw.n_tracks);
0233   assert(mat_traces_fw.size() == trk_stats_fw.n_tracks);
0234   assert(mat_records_fw.size() == trk_stats_fw.n_tracks);
0235   assert(mat_records_fw.size() == mat_records_bw.size());
0236   assert(step_traces_fw.size() == step_traces_bw.size());
0237   assert(mat_traces_fw.size() == mat_traces_bw.size());
0238
0239   // Check stats
0240   n_tracks = static_cast<double>(trk_stats_bw.n_tracks);
0241   if (static_cast<double>(trk_stats_bw.n_tracks_w_holes) / n_tracks >
0242       cfg.max_percent_missed / 100.) {
0243     DETRAY_ERROR_HOST("Too many tracks with missing surfaces");
0244     test_successful = false;
0245   }
0246   if (static_cast<double>(trk_stats_bw.n_tracks_w_extra) / n_tracks >
0247       cfg.max_percent_additional / 100.) {
0248     DETRAY_ERROR_HOST("Too many tracks with additional surfaces");
0249     test_successful = false;
0250   }
0251
0252   std::cout << "\n-----------------------------------\n"
0253             << "MATERIAL TRACE - No KF" << std::endl
0254             << "-----------------------------------\n";
0255
0256   constexpr double rel_mat_error{0.01};
0257
0258   // Material traces contain different surfaces
0259   std::size_t n_bad_comp{0u};
0260   // Overall integrated material differs while surface seq. is identical
0261   std::size_t n_diff_mat{0u};
0262
0263   // Loop over tracks
0264   for (std::size_t i = 0u; i < mat_records_fw.size(); ++i) {
0265     // No material on that track (e.g. detector model without material)
0266     if (mat_traces_fw[i].empty() && mat_traces_bw[i].empty()) {
0267       continue;
0268     }
0269
0270     std::remove_cvref_t<decltype(mat_traces_bw[i])> inv_mat_trace_bw{};
0271     if (!mat_traces_bw[i].empty()) {
0272       inv_mat_trace_bw.resize(mat_traces_bw[i].size());
0273
0274       // Revert the backward trace to compare to the forward trace
0275       std::ranges::reverse_copy(mat_traces_bw[i], inv_mat_trace_bw.begin());
0276
0277       assert(mat_traces_bw[i].size() == inv_mat_trace_bw.size());
0278       assert(mat_traces_bw[i].front().bcd == inv_mat_trace_bw.back().bcd);
0279     }
0280
0281     // Compare the material traces and total integrated material per trk
0282     const auto [is_bad_comp, is_diff_mat] = material_validator::compare_traces(
0283         mat_traces_fw[i], mat_records_fw[i], inv_mat_trace_bw,
0284         mat_records_bw[i], i, rel_mat_error, test_cfg.verbose());
0285
0286     if (is_bad_comp) {
0287       n_bad_comp++;
0288     }
0289     if (is_diff_mat) {
0290       n_diff_mat++;
0291     }
0292   }
0293
0294   std::cout << "Total no. tracks with diff. material: "
0295             << (n_bad_comp + n_diff_mat) << " ("
0296             << 100. * static_cast<double>(n_bad_comp + n_diff_mat) / n_tracks
0297             << "%)" << std::endl;
0298
0299   std::cout << "No. identical tracks with diff. material: " << n_diff_mat
0300             << " (" << 100. * static_cast<double>(n_diff_mat) / n_tracks << "%)"
0301             << std::endl;
0302   std::cout << "-----------------------------------\n" << std::endl;
0303
0304   // Trigger test failures
0305   if (n_diff_mat != 0) {
0306     DETRAY_ERROR_HOST("" << n_diff_mat << " tracks have differing material");
0307     test_successful = false;
0308   }
0309
0310   return test_successful;
0311 }
0312
0313 }  // namespace detray