detray/builders/grid_builder.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/builders/bin_fillers.hpp"
0013 #include "detray/builders/detail/radius_getter.hpp"
0014 #include "detray/builders/grid_factory.hpp"
0015 #include "detray/builders/surface_factory_interface.hpp"
0016 #include "detray/builders/volume_builder.hpp"
0017 #include "detray/builders/volume_builder_interface.hpp"
0018 #include "detray/geometry/tracking_volume.hpp"
0019 #include "detray/navigation/accelerators/concepts.hpp"
0020 #include "detray/navigation/accelerators/spatial_grid.hpp"
0021 #include "detray/utils/grid/concepts.hpp"
0022 #include "detray/utils/logging.hpp"
0023
0024 // System include(s)
0025 #include <algorithm>
0026 #include <cassert>
0027 #include <memory>
0028 #include <vector>
0029
0030 namespace detray {
0031
0032 /// @brief Build a grid of a certain shape.
0033 ///
0034 /// Decorator class to a volume builder that adds a grid as the volumes
0035 /// geometry accelerator structure.
0036 template <typename detector_t, concepts::grid grid_t,
0037           typename bin_filler_t = fill_by_pos,
0038           typename grid_factory_t = grid_factory_type<grid_t>>
0039 class grid_builder : public volume_decorator<detector_t> {
0040   using link_id_t = typename detector_t::volume_type::object_id;
0041
0042   // Is the given grid type already an acceleration structure?
0043   using spatial_grid_t =
0044       std::conditional_t<concepts::surface_accelerator<grid_t>, grid_t,
0045                          spatial_grid_impl<grid_t>>;
0046   using spatial_grid_owning_t = typename spatial_grid_t::template type<true>;
0047
0048  public:
0049   using detector_type = detector_t;
0050   using algebra_type = typename detector_t::algebra_type;
0051   using value_type = typename detector_type::surface_type;
0052   using scalar_type = dscalar<algebra_type>;
0053
0054   /// Decorate a volume with a surface accelerator grid
0055   DETRAY_HOST
0056   explicit grid_builder(
0057       std::unique_ptr<volume_builder_interface<detector_t>> vol_builder)
0058       : volume_decorator<detector_t>(std::move(vol_builder)) {
0059     DETRAY_VERBOSE_HOST("Add grid builder to volume: " << this->name());
0060
0061     // The grid builder provides an acceleration structure to the
0062     // volume, so don't add sensitive surfaces to the brute force method
0063     if (this->get_builder()) {
0064       this->has_accel(true);
0065     }
0066   }
0067
0068   /// Should the passive surfaces be added to the grid ?
0069   void set_add_passives(bool is_add_passive = false) {
0070     m_add_passives = is_add_passive;
0071
0072     DETRAY_VERBOSE_HOST("Add passive surfaces to grid: " << std::boolalpha
0073                                                          << m_add_passives
0074                                                          << std::noboolalpha);
0075   }
0076
0077   /// Set the surface category this grid should contain (type id in the
0078   /// accelerator link in the volume descriptor)
0079   void set_type(std::size_t sf_id) { set_type(static_cast<link_id_t>(sf_id)); }
0080
0081   /// Set the surface category this grid should contain (type id in the
0082   /// accelerator link in the volume descriptor)
0083   void set_type(link_id_t sf_id) {
0084     // Exclude zero, it is reserved for the brute force method
0085     assert(static_cast<int>(sf_id) > 0);
0086     // Make sure the id fits in the volume accelerator link
0087     assert(sf_id < link_id_t::e_size);
0088
0089     m_id = sf_id;
0090   }
0091
0092   /// Delegate init call depending on @param span type
0093   template <typename grid_shape_t>
0094   DETRAY_HOST void init_grid(
0095       const mask<grid_shape_t, algebra_type> &m,
0096       const darray<std::size_t, grid_t::dim> &n_bins,
0097       const std::vector<std::pair<typename grid_t::loc_bin_index, dindex>>
0098           &bin_capacities = {},
0099       const darray<std::vector<scalar_type>, grid_t::dim> &ax_bin_edges =
0100           darray<std::vector<scalar_type>, grid_t::dim>()) {
0101     DETRAY_VERBOSE_HOST("Initialize surface grid...");
0102
0103     static_assert(
0104         std::is_same_v<typename grid_shape_t::template local_frame_type<
0105                            typename detector_t::algebra_type>,
0106                        typename grid_t::local_frame_type>,
0107         "Mask has incorrect shape");
0108
0109     // Build spatial grid from grid utility and (so far empty) mask
0110     m_grid = spatial_grid_owning_t{m_factory.template new_grid<grid_t>(
0111                                        m, n_bins, bin_capacities, ax_bin_edges),
0112                                    typename spatial_grid_t::mask_type{}};
0113
0114     DETRAY_VERBOSE_HOST("Created empty grid:\n"
0115                         << DETRAY_TYPENAME(grid_t) << "\n"
0116                         << m_grid.axes());
0117   }
0118
0119   /// Build the empty grid from axis parameters
0120   DETRAY_HOST void init_grid(
0121       const std::vector<scalar_type> &spans,
0122       const std::vector<std::size_t> &n_bins,
0123       const std::vector<std::pair<typename grid_t::loc_bin_index, dindex>>
0124           &bin_capacities = {},
0125       const std::vector<std::vector<scalar_type>> &ax_bin_edges =
0126           std::vector<std::vector<scalar_type>>()) {
0127     DETRAY_VERBOSE_HOST("Initialize surface grid...");
0128
0129     // Build spatial grid from grid utility and (so far empty) mask
0130     m_grid =
0131         spatial_grid_owning_t{m_factory.template new_grid<grid_t>(
0132                                   spans, n_bins, bin_capacities, ax_bin_edges),
0133                               typename spatial_grid_t::mask_type{}};
0134
0135     DETRAY_VERBOSE_HOST("Created empty grid:\n"
0136                         << DETRAY_TYPENAME(grid_t) << "\n"
0137                         << m_grid.axes());
0138   }
0139
0140   /// Fill grid from existing volume using a bin filling strategy
0141   /// This can also be called without a volume builder
0142   template <typename volume_type, typename... Args>
0143   DETRAY_HOST void fill_grid(
0144       const detector_t &det, const volume_type &vol,
0145       const typename detector_t::geometry_context ctx = {},
0146       const bin_filler_t bin_filler = {}, Args &&...args) {
0147     bin_filler(m_grid, det, vol, ctx, args...);
0148   }
0149
0150   /// Fill grid from externally provided surfaces - temporary solution until
0151   /// the volume builders can be deployed in the toy detector
0152   template <typename volume_type, typename surface_container_t,
0153             typename transform_container_t, typename mask_container_t,
0154             typename... Args>
0155   DETRAY_HOST void fill_grid(
0156       const volume_type &vol, const surface_container_t &surfaces,
0157       const transform_container_t &transforms, const mask_container_t &masks,
0158       const typename detector_t::geometry_context ctx = {},
0159       const bin_filler_t bin_filler = {}, Args &&...args) {
0160     bin_filler(m_grid, vol, surfaces, transforms, masks, ctx, args...);
0161   }
0162
0163   /// Add the volume and the grid to the detector @param det
0164   DETRAY_HOST
0165   auto build(detector_t &det, typename detector_t::geometry_context ctx = {}) ->
0166       typename detector_t::volume_type * override {
0167     DETRAY_VERBOSE_HOST("Build surface grid...");
0168
0169     DETRAY_VERBOSE_HOST(
0170         " -> Defer to other builders to get complete surface descriptors "
0171         "first:");
0172
0173     using surface_desc_t = typename detector_t::surface_type;
0174
0175     // Add the surfaces (portals and/or passives) that are owned by the vol
0176     typename detector_t::volume_type *vol_ptr =
0177         volume_decorator<detector_t>::build(det, ctx);
0178
0179     DETRAY_VERBOSE_HOST("Resume building with updated surface descriptors");
0180
0181     // Find the surfaces that should be filled into the grid
0182     const auto vol = tracking_volume{det, vol_ptr->index()};
0183
0184     // Grid has not been filled previously, fill it automatically
0185     if (m_grid.size() == 0u) {
0186       DETRAY_DEBUG_HOST("-> Filling grid automatically...");
0187
0188       std::vector<surface_desc_t> surfaces{};
0189       for (auto &sf_desc : vol.surfaces()) {
0190         if (sf_desc.is_sensitive() ||
0191             (m_add_passives && sf_desc.is_passive())) {
0192           surfaces.push_back(sf_desc);
0193         }
0194       }
0195
0196       this->fill_grid(
0197           tracking_volume{det, volume_decorator<detector_t>::operator()()},
0198           surfaces, det.transform_store(), det.mask_store(), ctx);
0199     } else {
0200       DETRAY_DEBUG_HOST("-> Grid is prefilled...");
0201
0202       // The grid is prefilled with surface descriptors that contain the
0203       // correct LOCAL surface indices per bin (e.g. from file IO).
0204       // Now add the rest of the linking information, which is only
0205       // available after the volume builder ran
0206       for (surface_desc_t &sf_desc : m_grid.all()) {
0207         assert(!detail::is_invalid_value(sf_desc.index()));
0208
0209         dindex glob_idx{vol_ptr->to_global_sf_index(sf_desc.index())};
0210         const auto &new_sf_desc = det.surface(glob_idx);
0211
0212         assert(new_sf_desc.index() == glob_idx);
0213         assert(!new_sf_desc.identifier().is_invalid());
0214
0215         sf_desc = new_sf_desc;
0216       }
0217     }
0218
0219     // For cylinder grids, additional mask information is needed
0220     if constexpr (concepts::cylindrical<typename grid_t::local_frame_type>) {
0221       assert(vol_ptr->id() == volume_id::e_cylinder);
0222
0223       constexpr auto inv_vol_link{
0224           detray::detail::invalid_value<std::uint8_t>()};
0225       constexpr auto inf{std::numeric_limits<scalar_type>::max()};
0226
0227       DETRAY_DEBUG_HOST("Find radius for cylinder grid...");
0228
0229       // Passive surfaces could be in the brute force finder, but no
0230       // sensitive surfaces, since the volume has a grid. Their radii are,
0231       // however, always within the interval of the portal radii
0232       std::vector<scalar_type> radii{};
0233       for (auto pt_desc : vol.portals()) {
0234         auto r = det.mask_store().template visit<detail::outer_radius_getter>(
0235             pt_desc.mask());
0236
0237         if (r.has_value()) {
0238           DETRAY_DEBUG_HOST("Found radius " << *r
0239                                             << " mm of portal: " << pt_desc);
0240           radii.push_back(*r);
0241         }
0242       }
0243
0244       scalar_type grid_r{0.f};
0245       if (!radii.empty()) {
0246         const scalar_type inner_r{*std::ranges::min_element(radii)};
0247         const scalar_type outer_r{*std::ranges::max_element(radii)};
0248
0249         grid_r = 0.5f * (inner_r + outer_r);
0250       }
0251
0252       typename spatial_grid_t::mask_type cyl_mask{inv_vol_link, grid_r, -inf,
0253                                                   inf};
0254
0255       DETRAY_DEBUG_HOST("Setting mask for cylinder grid: " << cyl_mask);
0256
0257       m_grid.mask(cyl_mask);
0258     }
0259
0260     constexpr auto gid{types::id<typename detector_t::accel, spatial_grid_t>};
0261     const dindex grid_idx{det.accelerator_store().template size<gid>()};
0262
0263     DETRAY_DEBUG_HOST("Adding grid to volume in detector. Surface type: "
0264                       << m_id << ", grid: " << gid
0265                       << ", grid idx: " << grid_idx);
0266
0267     // Set: contained surface type, grid type, grid instance index
0268     vol_ptr->set_accel_link(m_id, gid, grid_idx);
0269
0270     // Add to detector
0271     det._accelerators.template push_back<gid>(m_grid);
0272
0273     DETRAY_DEBUG_HOST("Accelerator link: " << vol_ptr->accel_link());
0274     DETRAY_VERBOSE_HOST("Successfully built "
0275                         << gid << " for volume: " << this->name());
0276
0277     DETRAY_DEBUG_HOST("Finished grid: " << m_grid);
0278
0279     return vol_ptr;
0280   }
0281
0282   /// @returns access to the new grid
0283   DETRAY_HOST
0284   auto &get() { return m_grid; }
0285
0286  private:
0287   link_id_t m_id{link_id_t::e_sensitive};
0288   grid_factory_t m_factory{};
0289   // Data owning grid type, so that surface data can be filled into memory
0290   spatial_grid_owning_t m_grid{};
0291   bin_filler_t m_bin_filler{};
0292   bool m_add_passives{false};
0293 };
0294
0295 /// Grid builder from single components
0296 template <typename detector_t,
0297           template <class, template <std::size_t> class,
0298                     typename> class grid_factory_t,
0299           typename grid_shape_t, typename bin_t,
0300           template <std::size_t> class serializer_t,
0301           axis::bounds e_bounds = axis::bounds::e_closed,
0302           template <typename, typename> class... binning_ts>
0303 using grid_builder_type = grid_builder<
0304     detector_t,
0305     typename grid_factory_t<bin_t, serializer_t,
0306                             typename detector_t::algebra_type>::
0307         template grid_type<axes<grid_shape_t, e_bounds, binning_ts...>>,
0308     grid_factory_t<bin_t, serializer_t, typename detector_t::algebra_type>>;
0309
0310 }  // namespace detray