detray/utils/consistency_checker.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/geometry/tracking_volume.hpp"
0014 #include "detray/material/concepts.hpp"
0015 #include "detray/material/predefined_materials.hpp"
0016 #include "detray/utils/logging.hpp"
0017 #include "detray/utils/ranges.hpp"
0018 #include "detray/utils/type_registry.hpp"
0019
0020 // System include(s)
0021 #include <iostream>
0022 #include <sstream>
0023 #include <stdexcept>
0024 #include <utility>
0025
0026 namespace detray::detail {
0027
0028 /// Checks every collection in a multistore to be empty and prints a warning
0029 template <typename store_t, std::size_t... I>
0030 void report_empty(const store_t &store,
0031                   [[maybe_unused]] const std::string &store_name,
0032                   std::index_sequence<I...> /*seq*/) {
0033   ((store.template empty<types::id_cast<typename store_t::value_types, I>>()
0034         ?
0035 #if DETRAY_LOG_LVL < 0
0036         std::clog << ""
0037 #else
0038         // The log macro does not compile here...
0039         std::clog << "DETRAY WARNING (HOST): " << __FILENAME__ << ":"
0040                   << __LINE__ << " " << store_name
0041                   << " has empty collection no. " << I << std::endl
0042 #endif
0043         : std::clog << ""),
0044    ...);
0045 }
0046
0047 /// A functor that checks the surface descriptor and correct volume index in
0048 /// every acceleration data structure for a given volume
0049 struct surface_checker {
0050   /// Test the contained surfaces for consistency
0051   template <typename detector_t>
0052   DETRAY_HOST_DEVICE void operator()(
0053       const typename detector_t::surface_type &sf_descr, const detector_t &det,
0054       const dindex vol_idx, const typename detector_t::name_map &names) const {
0055     const auto sf = geometry::surface{det, sf_descr};
0056     const auto vol = tracking_volume{det, vol_idx};
0057     std::stringstream err_stream{};
0058     err_stream << "Volume \"" << vol.name(names) << "\":\n";
0059
0060     if (!sf.self_check(err_stream)) {
0061       throw std::invalid_argument(err_stream.str());
0062     }
0063
0064     if (sf.volume() != vol_idx) {
0065       err_stream << "Incorrect volume index on surface: vol. index " << vol_idx
0066                  << ", sf: " << sf;
0067
0068       DETRAY_FATAL_HOST(err_stream.str());
0069       throw std::invalid_argument(err_stream.str());
0070     }
0071
0072     // Does the mask link to an existing volume?
0073     const auto vol_links = sf.volume_links();
0074     for (const auto vol_link : vol_links) {
0075       if (!detail::is_invalid_value(vol_link) &&
0076           (vol_link >= det.volumes().size())) {
0077         err_stream << "Incorrect volume link to non-existent volume "
0078                    << vol_link << " on surface " << sf;
0079
0080         DETRAY_FATAL_HOST(err_stream.str());
0081         throw std::invalid_argument(err_stream.str());
0082       }
0083     }
0084
0085     // A passive surface should have material, if the detector was
0086     // constructed with material
0087     if (!det.material_store().all_empty() &&
0088         (sf.is_passive() && !sf.has_material())) {
0089       DETRAY_WARN_HOST(err_stream.str());
0090       DETRAY_WARN_HOST("Passive surface without material: " << sf_descr
0091                                                             << "\n");
0092     }
0093
0094     // Check that the same surface is registered in the detector surface
0095     // lookup
0096     const auto sf_from_lkp =
0097         geometry::surface{det, det.surface(sf.identifier())};
0098     if (sf_from_lkp != sf) {
0099       err_stream << "Surfaces in volume and detector lookups "
0100                  << "differ:\n In volume: " << vol
0101                  << "\nFound in detector surface lookup: " << sf_from_lkp;
0102
0103       DETRAY_FATAL_HOST(err_stream.str());
0104       throw std::runtime_error(err_stream.str());
0105     }
0106   }
0107
0108   /// Test whether a given surface @param check_desc is properly registered at
0109   /// least once in one of the volume acceleration data structures
0110   ///
0111   /// @param ref_descr one of the surfaces in the volumes acceleration data
0112   /// @param check_descr surface that we are searching for
0113   /// @param success communicate success to the outside
0114   template <typename sf_descr_t, typename sf_source_descr_t>
0115   DETRAY_HOST_DEVICE void operator()(const sf_descr_t &ref_descr,
0116                                      const sf_source_descr_t &check_descr,
0117                                      bool &success) const {
0118     // Check that the other surfaces in the acceleration structure belong
0119     // there The volume index of the check_descr must be checked to be
0120     // correct beforehand, e.g. by the call operator above
0121     if (ref_descr.volume() != check_descr.volume()) {
0122       std::stringstream err_stream{};
0123       err_stream << "Inconsistent volume index in accel. structure: "
0124                  << "Expected: " << check_descr.volume()
0125                  << ". Got surface in accel. structure: " << ref_descr;
0126
0127       DETRAY_FATAL_HOST(err_stream.str());
0128       throw std::invalid_argument(err_stream.str());
0129     }
0130
0131     // Check if it is the surface we are looking for
0132     if (ref_descr == check_descr) {
0133       success = true;
0134     }
0135   }
0136 };
0137
0138 /// A functor that checks the material parametrization for a surface/volume
0139 struct material_checker {
0140   /// Error message for material consistency check
0141   template <typename material_t>
0142   void throw_material_error(const std::string &type, const dindex idx,
0143                             const material_t &mat) const {
0144     std::stringstream err_stream{};
0145     err_stream << "Invalid material found in: " << type << " at index " << idx
0146                << ": " << mat;
0147     DETRAY_FATAL_HOST(err_stream.str());
0148     throw std::invalid_argument(err_stream.str());
0149   }
0150
0151   /// Test whether a given material map contains invalid material
0152   ///
0153   /// @param material_coll collection of material grids
0154   /// @param idx the specific grid to be checked
0155   /// @param id type id of the material grid collection
0156   template <typename material_coll_t, typename index_t, typename id_t>
0157     requires concepts::material_map<typename material_coll_t::value_type>
0158   DETRAY_HOST_DEVICE void operator()(const material_coll_t &material_coll,
0159                                      const index_t idx, const id_t id) const {
0160     try {
0161       const auto mat_map = material_coll.at(idx);
0162
0163       // Check whether there are any entries in the bins
0164       if (mat_map.size() == 0u) {
0165         std::stringstream err_stream{};
0166         err_stream << "Empty material grid: " << static_cast<int>(id)
0167                    << " at index " << idx;
0168         DETRAY_FATAL_HOST(err_stream.str());
0169         throw std::invalid_argument(err_stream.str());
0170       } else {
0171         for (const auto &bin : mat_map.bins()) {
0172           if (bin.size() == 0u) {
0173             std::stringstream err_stream{};
0174             err_stream << "Empty material bin: " << static_cast<int>(id)
0175                        << " at index " << idx;
0176             DETRAY_FATAL_HOST(err_stream.str());
0177             throw std::invalid_argument(err_stream.str());
0178           }
0179         }
0180       }
0181     } catch (std::out_of_range & /*unused*/) {
0182       std::stringstream err_stream{};
0183       err_stream << "Out of range material access in: "
0184                  << "binned material collection at index " << idx;
0185       DETRAY_FATAL_HOST(err_stream.str());
0186       throw std::invalid_argument(err_stream.str());
0187     }
0188   }
0189
0190   /// Test whether a given collection of material contains invalid material
0191   ///
0192   /// @param material_coll collection of material slabs/rods/raw mat
0193   /// @param idx the specific instance to be checked
0194   template <typename material_coll_t, typename index_t, typename id_t>
0195     requires concepts::homogeneous_material<
0196         typename material_coll_t::value_type>
0197   DETRAY_HOST_DEVICE void operator()(const material_coll_t &material_coll,
0198                                      const index_t idx,
0199                                      const id_t /*unused*/) const {
0200     using material_t = typename material_coll_t::value_type;
0201     using scalar_t = typename material_t::scalar_type;
0202
0203     try {
0204       const material_t &mat = material_coll.at(idx);
0205
0206       // Homogeneous volume material
0207       if constexpr (concepts::material_params<material_t>) {
0208         if (mat == detray::vacuum<scalar_t>{}) {
0209           throw_material_error("homogeneous volume material", idx, mat);
0210         }
0211
0212       } else {
0213         // Material slabs and rods
0214         if (!mat) {
0215           throw_material_error("homogeneous surface material", idx, mat);
0216         }
0217       }
0218     } catch (std::out_of_range & /*unused*/) {
0219       std::stringstream err_stream{};
0220       err_stream << "Out of range material access in: "
0221                  << "homogeneous material collection at index " << idx;
0222       DETRAY_FATAL_HOST(err_stream.str());
0223       throw std::invalid_argument(err_stream.str());
0224     }
0225   }
0226 };
0227
0228 /// @brief Checks whether the data containers of a detector are empty
0229 ///
0230 /// In case the default metadata is used, the unused containers are allowed to
0231 /// be empty.
0232 template <typename detector_t>
0233 inline void check_empty(const detector_t &det, const bool verbose) {
0234   // Check if there is at least one portal in the detector
0235   auto find_portals = [&det]() {
0236     if (det.portals().empty()) {
0237       return false;
0238     }
0239     // In the brute force finder, also other surfaces can be contained, e.g.
0240     // passive surfaces (depends on the detector)
0241     return std::ranges::any_of(
0242         det.portals(), [](auto pt_desc) { return pt_desc.is_portal(); });
0243   };
0244
0245   // TODO: Check for empty volume acceleration structures
0246
0247   // Fatal errors
0248   if (det.volumes().empty()) {
0249     std::string err_str{"No volumes in detector"};
0250     DETRAY_FATAL_HOST(err_str);
0251     throw std::runtime_error(err_str);
0252   }
0253   if (det.surfaces().empty()) {
0254     std::string err_str{"No surfaces found"};
0255     DETRAY_FATAL_HOST(err_str);
0256     throw std::runtime_error(err_str);
0257   }
0258   if (det.transform_store().empty()) {
0259     std::string err_str{"No transforms in detector"};
0260     DETRAY_FATAL_HOST(err_str);
0261     throw std::runtime_error(err_str);
0262   }
0263   if (det.mask_store().all_empty()) {
0264     std::string err_str{"No masks in detector"};
0265     DETRAY_FATAL_HOST(err_str);
0266     throw std::runtime_error(err_str);
0267   }
0268   if (!find_portals()) {
0269     std::string err_str{"No portals in detector"};
0270     DETRAY_FATAL_HOST(err_str);
0271     throw std::runtime_error(err_str);
0272   }
0273
0274   // Warnings
0275
0276   // Check the material description
0277   if (det.material_store().all_empty()) {
0278     DETRAY_WARN_HOST("No material in detector");
0279   } else if (verbose) {
0280     // Check for empty material collections
0281     detail::report_empty(
0282         det.material_store(), "material store",
0283         std::make_index_sequence<detector_t::material::n_types>{});
0284   }
0285
0286   // Check for empty acceleration data structure collections (e.g. grids)
0287   if (verbose) {
0288     // Check for empty mask collections
0289     detail::report_empty(
0290         det.mask_store(), "mask store",
0291         std::make_index_sequence<detector_t::masks::n_types>{});
0292
0293     detail::report_empty(
0294         det.accelerator_store(), "acceleration data structures store",
0295         std::make_index_sequence<detector_t::accel::n_types>{});
0296   }
0297
0298   // TODO: Implement volume acceleration structure check
0299 }
0300
0301 /// @brief Checks the internal consistency of a detector
0302 template <typename detector_t>
0303 inline bool check_consistency(const detector_t &det, const bool verbose = false,
0304                               const typename detector_t::name_map &names = {}) {
0305   DETRAY_INFO_HOST("Checking detector consistency...");
0306
0307   check_empty(det, verbose);
0308
0309   // Check the volumes
0310   constexpr bool portal_eq_passive{
0311       detector_t::volume_type::object_id::e_portal ==
0312       detector_t::volume_type::object_id::e_passive};
0313
0314   for (const auto &[idx, vol_desc] : detray::views::enumerate(det.volumes())) {
0315     const auto vol = tracking_volume{det, vol_desc};
0316
0317     std::stringstream err_stream{};
0318     err_stream << "Volume \"" << vol.name(names) << "\":\n";
0319
0320     // Check that nothing is obviously broken
0321     if (!vol.self_check(err_stream)) {
0322       DETRAY_FATAL_HOST(err_stream.str());
0323       throw std::invalid_argument(err_stream.str());
0324     }
0325
0326     // Check consistency in the context of the owning detector
0327     if (vol.index() != idx) {
0328       err_stream << "Incorrect volume index! Found volume:\n"
0329                  << vol << "\nat index " << idx;
0330
0331       DETRAY_FATAL_HOST(err_stream.str());
0332       throw std::invalid_argument(err_stream.str());
0333     }
0334
0335     // Check that only surfaces belonging to this volume are in the surface
0336     // range of the detector
0337
0338     // Check volume indices
0339     for (const auto sf_desc : vol.template surfaces<surface_id::e_all>()) {
0340       if (sf_desc.volume() != vol.index()) {
0341         err_stream << "Surface of other volume detected:\n\nvolume " << vol
0342                    << "\n\nsurface " << sf_desc;
0343
0344         DETRAY_FATAL_HOST(err_stream.str());
0345         throw std::invalid_argument(err_stream.str());
0346       }
0347     }
0348     // Check surface type: portal or passive
0349     for (const auto sf_desc : vol.portals()) {
0350       if (!portal_eq_passive && !sf_desc.is_portal()) {
0351         err_stream << "Non-portal surface discovered in portal "
0352                       "range: "
0353                    << vol_desc.template sf_link<surface_id::e_portal>()
0354                    << "\n\nvolume " << vol << "\n\nsurface " << sf_desc;
0355
0356         DETRAY_FATAL_HOST(err_stream.str());
0357         throw std::invalid_argument(err_stream.str());
0358       }
0359
0360       if (portal_eq_passive && !(sf_desc.is_portal() || sf_desc.is_passive())) {
0361         err_stream << "Sensitive surface discovered in portal/passive "
0362                       "range: "
0363                    << vol_desc.template sf_link<surface_id::e_portal>()
0364                    << "\n\nvolume " << vol << "\n\nsurface " << sf_desc;
0365
0366         DETRAY_FATAL_HOST(err_stream.str());
0367         throw std::invalid_argument(err_stream.str());
0368       }
0369     }
0370     // Check surface type: passive or portal
0371     for (const auto sf_desc : vol.template surfaces<surface_id::e_passive>()) {
0372       if (!portal_eq_passive && !sf_desc.is_passive()) {
0373         err_stream << "Non-passive surface discovered in passive "
0374                       "range: "
0375                    << vol_desc.template sf_link<surface_id::e_passive>()
0376                    << "\n\nvolume " << vol << "\n\nsurface " << sf_desc;
0377
0378         DETRAY_FATAL_HOST(err_stream.str());
0379         throw std::invalid_argument(err_stream.str());
0380       }
0381
0382       if (portal_eq_passive && !(sf_desc.is_portal() || sf_desc.is_passive())) {
0383         err_stream << "Sensitive surface discovered in portal/passive "
0384                       "range: "
0385                    << vol_desc.template sf_link<surface_id::e_passive>()
0386                    << "\n\nvolume " << vol << "\n\nsurface " << sf_desc;
0387
0388         DETRAY_FATAL_HOST(err_stream.str());
0389         throw std::invalid_argument(err_stream.str());
0390       }
0391     }
0392     // Check surface type: sensitive
0393     for (const auto sf_desc :
0394          vol.template surfaces<surface_id::e_sensitive>()) {
0395       if (!sf_desc.is_sensitive()) {
0396         err_stream << "Non-sensitive surface discovered in sensitive "
0397                       "range: "
0398                    << vol_desc.template sf_link<surface_id::e_sensitive>()
0399                    << "\n\nvolume " << vol << "\n\nsurface " << sf_desc;
0400
0401         DETRAY_FATAL_HOST(err_stream.str());
0402         throw std::invalid_argument(err_stream.str());
0403       }
0404     }
0405
0406     // Go through the acceleration data structures and check the surfaces
0407     vol.template visit_surfaces<surface_id::e_all, detail::surface_checker>(
0408         det, vol.index(), names);
0409
0410     // Check the volume material, if present
0411     if (vol.has_material()) {
0412       vol.template visit_material<detail::material_checker>(
0413           vol_desc.material().id());
0414     }
0415   }
0416
0417   // Check the surfaces in the detector's surface lookup
0418   for (const auto &[idx, sf_desc] : detray::views::enumerate(det.surfaces())) {
0419     const auto sf = geometry::surface{det, sf_desc};
0420     const auto vol = tracking_volume{det, sf.volume()};
0421
0422     std::stringstream err_stream{};
0423     err_stream << "Volume \"" << vol.name(names) << "\":\n";
0424
0425     // Check that nothing is obviously broken
0426     if (!sf.self_check(err_stream)) {
0427       err_stream << "\nat surface no. " << std::to_string(idx);
0428       throw std::invalid_argument(err_stream.str());
0429     }
0430
0431     // Check consistency in the context of the owning detector
0432     if (sf.index() != idx) {
0433       err_stream << "Incorrect surface index! Found surface:\n"
0434                  << sf << "\nat index " << idx;
0435
0436       DETRAY_FATAL_HOST(err_stream.str());
0437       throw std::invalid_argument(err_stream.str());
0438     }
0439
0440     // Check that the surface can be found in its volume's acceleration
0441     // data structures (if there are no grids, must at least be in the
0442     // brute force method)
0443     bool is_registered = false;
0444     vol.template visit_surfaces<surface_id::e_all, detail::surface_checker>(
0445         sf_desc, is_registered);
0446
0447     if (!is_registered) {
0448       err_stream << "Found surface that is not part of its "
0449                  << "volume's navigation acceleration data structures:\n"
0450                  << "Surface: " << sf;
0451
0452       DETRAY_FATAL_HOST(err_stream.str());
0453       throw std::invalid_argument(err_stream.str());
0454     }
0455
0456     // Check the surface material, if present
0457     if (sf.has_material()) {
0458       DETRAY_DEBUG_HOST("Checking surface sf=" << sf);
0459       sf.template visit_material<detail::material_checker>(
0460           sf_desc.material().id());
0461     }
0462   }
0463
0464   DETRAY_INFO_HOST("Detector consistency check: OK");
0465
0466   return true;
0467 }
0468
0469 }  // namespace detray::detail