test/cpu/material_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/core/detector.hpp"
0013 #include "detray/tracks/tracks.hpp"
0014 #include "detray/utils/logging.hpp"
0015 #include "detray/utils/ranges.hpp"
0016
0017 // Detray test include(s)
0018 #include "detray/test/framework/fixture_base.hpp"
0019 #include "detray/test/framework/whiteboard.hpp"
0020 #include "detray/test/validation/material_validation_config.hpp"
0021 #include "detray/test/validation/material_validation_utils.hpp"
0022
0023 // Vecmem include(s)
0024 #include <vecmem/memory/host_memory_resource.hpp>
0025
0026 // System include(s)
0027 #include <iostream>
0028 #include <limits>
0029 #include <memory>
0030 #include <stdexcept>
0031 #include <string>
0032 #include <string_view>
0033
0034 namespace detray::test {
0035
0036 /// Run the material validation on the host
0037 struct run_material_validation {
0038   static constexpr std::string_view name{"cpu"};
0039
0040   template <typename detector_t>
0041   auto operator()(
0042       vecmem::memory_resource *host_mr, vecmem::memory_resource * /*mr*/,
0043       const detector_t &det, const propagation::config &cfg,
0044       const dvector<free_track_parameters<typename detector_t::algebra_type>>
0045           &tracks,
0046       const std::vector<std::size_t> & /*mask*/ = {}) {
0047     using scalar_t = dscalar<typename detector_t::algebra_type>;
0048
0049     typename detector_t::geometry_context gctx{};
0050
0051     dvector<material_validator::material_record<scalar_t>> mat_records{host_mr};
0052     mat_records.reserve(tracks.size());
0053
0054     dvector<dvector<material_validator::material_params<scalar_t>>>
0055         mat_steps_vec{host_mr};
0056
0057     mat_steps_vec.reserve(tracks.size());
0058
0059     for (const auto &[i, track] : detray::views::enumerate(tracks)) {
0060       auto [success, mat_record, mat_steps] =
0061           detray::material_validator::record_material(gctx, host_mr, det, cfg,
0062                                                       track);
0063       mat_records.push_back(mat_record);
0064       mat_steps_vec.push_back(std::move(mat_steps));
0065
0066       if (!success) {
0067         DETRAY_ERROR_HOST("Propagation failed for track "
0068                           << i << ": "
0069                           << "Material record may be incomplete!");
0070       }
0071     }
0072
0073     return std::make_tuple(std::move(mat_records), std::move(mat_steps_vec));
0074   }
0075 };
0076
0077 /// @brief Test class that runs the material validation for a given detector.
0078 ///
0079 /// @note The lifetime of the detector needs to be guaranteed outside this class
0080 template <typename detector_t, typename material_validator_t>
0081 class material_validation_impl : public test::fixture_base<> {
0082   using algebra_t = typename detector_t::algebra_type;
0083   using scalar_t = dscalar<algebra_t>;
0084   using free_track_parameters_t = free_track_parameters<algebra_t>;
0085   using material_record_t = material_validator::material_record<scalar_t>;
0086
0087  public:
0088   using fixture_type = test::fixture_base<>;
0089   using config = detray::test::material_validation_config<algebra_t>;
0090
0091   explicit material_validation_impl(
0092       const detector_t &det, const typename detector_t::name_map &names,
0093       const config &cfg = {}, std::shared_ptr<test::whiteboard> wb = nullptr,
0094       const typename detector_t::geometry_context gctx = {})
0095       : m_cfg{cfg},
0096         m_gctx{gctx},
0097         m_det{det},
0098         m_names{names},
0099         m_whiteboard{std::move(wb)} {
0100     if (!m_whiteboard) {
0101       throw std::invalid_argument("No white board was passed to " +
0102                                   m_cfg.name() + " test");
0103     }
0104
0105     // Name of the material scan data collection
0106     m_scan_data_name = m_det.name(m_names) + "_material_scan";
0107     m_track_data_name = m_det.name(m_names) + "_material_scan_tracks";
0108
0109     // Check that data is available in memory
0110     if (!m_whiteboard->exists(m_scan_data_name)) {
0111       throw std::invalid_argument(
0112           "Material validation: Could not find scan data on whiteboard."
0113           "Please run material scan first.");
0114     }
0115     if (!m_whiteboard->exists(m_track_data_name)) {
0116       throw std::invalid_argument(
0117           "Material validation: Could not find track data on whiteboard."
0118           "Please run material scan first.");
0119     }
0120   }
0121
0122   /// Run the check
0123   void TestBody() override {
0124     using namespace detray;
0125
0126     // Fetch the input data
0127     const auto &tracks =
0128         m_whiteboard->template get<dvector<free_track_parameters_t>>(
0129             m_track_data_name);
0130
0131     const auto &truth_mat_records =
0132         m_whiteboard->template get<dvector<material_record_t>>(
0133             m_scan_data_name);
0134
0135     DETRAY_INFO_HOST("Running material validation on: " << m_det.name(m_names)
0136                                                         << "...\n");
0137
0138     // only needed for device material steps allocations
0139     // @TODO: For now, guess how many surface might be encountered
0140     std::vector<std::size_t> capacities(tracks.size(), 80u);
0141
0142     // Run the propagation on device and record the accumulated material
0143     auto [mat_records, mat_steps] =
0144         material_validator_t{}(&m_host_mr, m_cfg.device_mr(), m_det,
0145                                m_cfg.propagation(), tracks, capacities);
0146
0147     // One material record per track
0148     ASSERT_EQ(tracks.size(), mat_records.size());
0149
0150     // Collect some statistics
0151     std::size_t n_tracks{0u};
0152     const scalar_t rel_error{m_cfg.relative_error()};
0153     for (std::size_t i = 0u; i < mat_records.size(); ++i) {
0154       if (n_tracks >= m_cfg.n_tracks()) {
0155         break;
0156       }
0157
0158       const auto &truth_mat = truth_mat_records[i];
0159       const auto &recorded_mat = mat_records[i];
0160
0161       auto get_rel_error = [](const scalar_t truth, const scalar_t rec) {
0162         constexpr scalar_t e{std::numeric_limits<scalar_t>::epsilon()};
0163
0164         if (truth <= e && rec <= e) {
0165           // No material for this ray => valid
0166           return scalar_t{0.f};
0167         } else if (truth <= e) {
0168           // Material found where none should be
0169           return detail::invalid_value<scalar_t>();
0170         } else {
0171           return math::fabs(truth - rec) / truth;
0172         }
0173       };
0174
0175       EXPECT_LT(get_rel_error(truth_mat.sX0, recorded_mat.sX0), rel_error)
0176           << "Track " << n_tracks << " (X0 / path): Truth " << truth_mat.sX0
0177           << ", Nav. " << recorded_mat.sX0;
0178       EXPECT_LT(get_rel_error(truth_mat.tX0, recorded_mat.tX0), rel_error)
0179           << "Track " << n_tracks << " (X0 / thickness): Truth "
0180           << truth_mat.tX0 << ", Nav. " << recorded_mat.tX0;
0181       EXPECT_LT(get_rel_error(truth_mat.sL0, recorded_mat.sL0), rel_error)
0182           << "Track " << n_tracks << " (L0 / path): Truth " << truth_mat.sL0
0183           << ", Nav. " << recorded_mat.sL0;
0184       EXPECT_LT(get_rel_error(truth_mat.tL0, recorded_mat.tL0), rel_error)
0185           << "Track " << n_tracks << " (L0 / thickness): Truth "
0186           << truth_mat.tL0 << ", Nav. " << recorded_mat.tL0;
0187
0188       ++n_tracks;
0189     }
0190
0191     std::clog << "-----------------------------------\n"
0192               << "Tested " << n_tracks << " tracks\n"
0193               << "-----------------------------------\n"
0194               << std::endl;
0195
0196     // Print accumulated material per track
0197     std::filesystem::path mat_path{m_cfg.material_file()};
0198     const auto data_path{mat_path.parent_path()};
0199
0200     auto file_name{data_path /
0201                    (m_det.name(m_names) + "_" + mat_path.stem().string() + "_" +
0202                     std::string(material_validator_t::name) + ".csv")};
0203
0204     material_validator::write_material(file_name.string(), mat_records);
0205   }
0206
0207  private:
0208   /// The configuration of this test
0209   config m_cfg;
0210   /// Vecmem memory resource for the host allocations
0211   vecmem::host_memory_resource m_host_mr{};
0212   /// Name of the input data collections
0213   std::string m_scan_data_name{""};
0214   std::string m_track_data_name{""};
0215   /// The geometry context to check
0216   typename detector_t::geometry_context m_gctx{};
0217   /// The detector to be checked
0218   const detector_t &m_det;
0219   /// Volume names
0220   const typename detector_t::name_map &m_names;
0221   /// Whiteboard to pin data
0222   std::shared_ptr<test::whiteboard> m_whiteboard{nullptr};
0223 };
0224
0225 template <typename detector_t>
0226 using material_validation =
0227     material_validation_impl<detector_t, run_material_validation>;
0228
0229 }  // namespace detray::test