test/cpu/material_scan.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/geometry/surface.hpp"
0013 #include "detray/tracks/ray.hpp"
0014 #include "detray/utils/logging.hpp"
0015
0016 // Detray IO include(s)
0017 #include "detray/io/utils/file_handle.hpp"
0018
0019 // Detray test include(s)
0020 #include "detray/test/common/track_generators.hpp"
0021 #include "detray/test/framework/fixture_base.hpp"
0022 #include "detray/test/framework/types.hpp"
0023 #include "detray/test/framework/whiteboard.hpp"
0024 #include "detray/test/validation/detector_scan_utils.hpp"
0025 #include "detray/test/validation/detector_scanner.hpp"
0026 #include "detray/test/validation/material_validation_utils.hpp"
0027
0028 // System include(s)
0029 #include <ios>
0030 #include <iostream>
0031 #include <memory>
0032 #include <string>
0033
0034 namespace detray::test {
0035
0036 /// @brief Test class that runs the material ray scan on a given detector.
0037 ///
0038 /// @note The lifetime of the detector needs to be guaranteed.
0039 template <typename detector_t>
0040 class material_scan : public test::fixture_base<> {
0041   using algebra_t = typename detector_t::algebra_type;
0042   using point2_t = dpoint2D<algebra_t>;
0043   using scalar_t = dscalar<algebra_t>;
0044   using ray_t = detail::ray<algebra_t>;
0045   using track_generator_t = uniform_track_generator<ray_t>;
0046
0047  public:
0048   using fixture_type = test::fixture_base<>;
0049
0050   struct config : public fixture_type::configuration {
0051     using trk_gen_config_t = typename track_generator_t::configuration;
0052
0053     std::string m_name{"material_scan"};
0054     trk_gen_config_t m_trk_gen_cfg{};
0055     /// Name of the output file, containing the complete ray material traces
0056     std::string m_material_file{"material_scan"};
0057     /// Perform overlaps removal (needed for detector converted from ACTS)
0058     bool m_overlaps_removal{true};
0059     /// Tolerance for overlaps
0060     float m_overlaps_tol{1e-4f * unit<float>::mm};
0061
0062     /// Getters
0063     /// @{
0064     const std::string &name() const { return m_name; }
0065     trk_gen_config_t &track_generator() { return m_trk_gen_cfg; }
0066     const trk_gen_config_t &track_generator() const { return m_trk_gen_cfg; }
0067     const std::string &material_file() const { return m_material_file; }
0068     bool overlaps_removal() const { return m_overlaps_removal; }
0069     float overlaps_tol() const { return m_overlaps_tol; }
0070     /// @}
0071
0072     /// Setters
0073     /// @{
0074     config &name(const std::string &n) {
0075       m_name = n;
0076       return *this;
0077     }
0078     config &material_file(const std::string &f) {
0079       m_material_file = f;
0080       return *this;
0081     }
0082     config &overlaps_removal(const bool o) {
0083       m_overlaps_removal = o;
0084       return *this;
0085     }
0086     config &overlaps_tol(const scalar_t tol) {
0087       m_overlaps_tol = static_cast<float>(tol);
0088       return *this;
0089     }
0090     /// @}
0091   };
0092
0093   explicit material_scan(const detector_t &det,
0094                          const typename detector_t::name_map &names,
0095                          const config &cfg = {},
0096                          std::shared_ptr<test::whiteboard> wb = nullptr,
0097                          const typename detector_t::geometry_context &gctx = {})
0098       : m_cfg{cfg},
0099         m_gctx{gctx},
0100         m_det{det},
0101         m_names{names},
0102         m_whiteboard{std::move(wb)} {
0103     if (!m_whiteboard) {
0104       throw std::invalid_argument("No white board was passed to " +
0105                                   m_cfg.name());
0106     }
0107   }
0108
0109   /// Run the ray scan
0110   void TestBody() override {
0111     using material_record_t = material_validator::material_record<scalar_t>;
0112
0113     std::size_t n_tracks{0u};
0114     auto ray_generator = track_generator_t(m_cfg.track_generator());
0115
0116     DETRAY_INFO_HOST("Running material scan on: "
0117                      << m_det.name(m_names) << "\n(" << ray_generator.size()
0118                      << " rays) ...\n");
0119
0120     // Trace material per ray
0121     dvector<free_track_parameters<algebra_t>> tracks{};
0122     tracks.reserve(ray_generator.size());
0123
0124     dvector<material_record_t> mat_records{};
0125     mat_records.reserve(ray_generator.size());
0126
0127     for (const auto ray : ray_generator) {
0128       // Record all intersections and surfaces along the ray
0129       auto intersection_record =
0130           detector_scanner::run<detray::ray_scan>(m_gctx, m_det, ray);
0131
0132       // Remove certain allowed duplications
0133       if (m_cfg.overlaps_removal()) {
0134         detector_scanner::overlaps_removal(intersection_record,
0135                                            m_cfg.overlaps_tol());
0136       }
0137
0138       if (intersection_record.empty()) {
0139         DETRAY_FATAL_HOST("Intersection trace empty for ray "
0140                           << n_tracks << "/" << ray_generator.size() << ": "
0141                           << ray);
0142         break;
0143       }
0144
0145       // Record track parameters
0146       tracks.push_back({ray.pos(), 0.f, ray.dir(), 0.f});
0147
0148       // New material record
0149       material_record_t mat_record{};
0150       mat_record.eta = vector::eta(ray.dir());
0151       mat_record.phi = vector::phi(ray.dir());
0152
0153       // Record material for this ray
0154       for (const auto &[i, record] :
0155            detray::views::enumerate(intersection_record)) {
0156         // Check whether this record has a successor
0157         if (i < intersection_record.size() - 1) {
0158           const auto &current_intr = record.intersection;
0159           const auto &next_intr = intersection_record[i + 1].intersection;
0160
0161           const bool is_same_intrs{next_intr == current_intr};
0162           const bool current_is_pt{current_intr.surface().is_portal()};
0163           const bool next_is_pt{next_intr.surface().is_portal()};
0164
0165           const bool is_dummy_record{i == 0};
0166
0167           // Prevent double counting of material on adjacent portals
0168           // (the navigator automatically skips the exit portal)
0169           if ((is_same_intrs && next_is_pt && current_is_pt) ||
0170               is_dummy_record) {
0171             continue;
0172           }
0173         }
0174
0175         // Don't count the last portal, because navigation terminates
0176         // before the material is counted
0177         if (detail::is_invalid_value(record.intersection.volume_link())) {
0178           continue;
0179         }
0180
0181         const auto sf = geometry::surface{m_det, record.intersection.surface()};
0182
0183         if (!sf.has_material()) {
0184           continue;
0185         }
0186
0187         const auto &p = record.intersection.local();
0188         const auto mat_params =
0189             sf.template visit_material<material_validator::get_material_params>(
0190                 point2_t{p[0], p[1]}, cos_angle(m_gctx, sf, ray.dir(), p));
0191
0192         const scalar_t seg{mat_params.path};
0193         const scalar_t t{mat_params.thickness};
0194         const scalar_t mx0{mat_params.mat_X0};
0195         const scalar_t ml0{mat_params.mat_L0};
0196
0197         if (mx0 > 0.f) {
0198           mat_record.sX0 += seg / mx0;
0199           mat_record.tX0 += t / mx0;
0200         } else {
0201           DETRAY_ERROR_HOST(
0202               "Encountered invalid X_0: " << mx0 << "\nOn surface: " << sf);
0203         }
0204         if (ml0 > 0.f) {
0205           mat_record.sL0 += seg / ml0;
0206           mat_record.tL0 += t / ml0;
0207         } else {
0208           DETRAY_ERROR_HOST(
0209               "Encountered invalid L_0: " << ml0 << "\nOn surface: " << sf);
0210         }
0211       }
0212
0213       if (mat_record.sX0 == 0.f || mat_record.sL0 == 0.f ||
0214           mat_record.tX0 == 0.f || mat_record.tL0 == 0.f) {
0215         DETRAY_VERBOSE_HOST("No material recorded for ray "
0216                             << n_tracks << "/" << ray_generator.size() << ": "
0217                             << ray);
0218       }
0219
0220       mat_records.push_back(mat_record);
0221
0222       ++n_tracks;
0223     }
0224
0225     std::clog << "-----------------------------------\n"
0226               << "Tested " << n_tracks << " tracks\n"
0227               << "-----------------------------------\n"
0228               << std::endl;
0229
0230     // Write recorded material to csv file
0231     const std::string det_name{m_det.name(m_names)};
0232     const auto material_path{std::filesystem::path{m_cfg.material_file()}};
0233     const auto data_path{material_path.parent_path()};
0234
0235     std::string file_name{det_name + "_" + material_path.stem().string() +
0236                           ".csv"};
0237     material_validator::write_material(data_path / file_name, mat_records);
0238
0239     // Pin data to whiteboard
0240     m_whiteboard->add(det_name + "_material_scan", std::move(mat_records));
0241     m_whiteboard->add(det_name + "_material_scan_tracks", std::move(tracks));
0242   }
0243
0244  private:
0245   /// The configuration of this test
0246   config m_cfg;
0247   /// The geometry context to scan
0248   typename detector_t::geometry_context m_gctx{};
0249   /// The detector to be checked
0250   const detector_t &m_det;
0251   /// Volume names
0252   const typename detector_t::name_map &m_names;
0253   /// Whiteboard to pin data
0254   std::shared_ptr<test::whiteboard> m_whiteboard{nullptr};
0255 };
0256
0257 }  // namespace detray::test