detray/builders/grid_factory.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/definitions/algebra.hpp"
0013 #include "detray/definitions/containers.hpp"
0014 #include "detray/definitions/units.hpp"
0015 #include "detray/geometry/concepts.hpp"
0016 #include "detray/geometry/mask.hpp"
0017 #include "detray/geometry/shapes.hpp"
0018 #include "detray/utils/grid/axis.hpp"
0019 #include "detray/utils/grid/concepts.hpp"
0020 #include "detray/utils/grid/detail/axis_helpers.hpp"
0021 #include "detray/utils/grid/grid.hpp"
0022 #include "detray/utils/grid/grid_collection.hpp"
0023 #include "detray/utils/grid/populators.hpp"
0024 #include "detray/utils/grid/serializers.hpp"
0025
0026 // Vecmem include(s)
0027 #include <vecmem/memory/memory_resource.hpp>
0028
0029 // System include(s)
0030 #include <cassert>
0031 #include <iostream>
0032 #include <string>
0033 #include <type_traits>
0034
0035 namespace detray {
0036
0037 /// @brief Provides functionality to instantiate grids and grid collections.
0038 ///
0039 /// @tparam bin_t type bins in the grid
0040 /// @tparam serialzier_t  type of the serializer to the storage representations
0041 /// @tparam algebra_t the matrix/vector/point types to use
0042 /// @tparam container_t the container types to use
0043 ///
0044 /// throughout:
0045 /// @note that bin_edges is only used for variable binning
0046 /// @note that if non-zero axis_spans are provided the values of the
0047 /// mask is overwritten
0048 template <typename bin_t, template <std::size_t> class serializer_t,
0049           concepts::algebra algebra_t>
0050 class grid_factory {
0051  public:
0052   // All grids are owning since they are used to fill the data
0053   static constexpr bool is_owning = true;
0054
0055   using bin_type = bin_t;
0056   template <typename grid_shape_t>
0057   using grid_type = grid<algebra_t, axes<grid_shape_t>, bin_type, serializer_t>;
0058   template <typename grid_shape_t>
0059   using loc_bin_index = typename grid_type<grid_shape_t>::loc_bin_index;
0060
0061   using scalar_type = dscalar<algebra_t>;
0062   template <typename T>
0063   using vector_type = host_container_types::template vector_type<T>;
0064   using algebra_type = algebra_t;
0065
0066   grid_factory() = default;
0067
0068   /// Takes the resource of the detector to allocate memory correctly
0069   explicit grid_factory(vecmem::memory_resource &resource)
0070       : m_resource(&resource) {}
0071
0072   //
0073   // annulus 2D
0074   //
0075   template <
0076       typename r_bounds = axis::closed<axis::label::e_r>,
0077       typename phi_bounds = axis::circular<>,
0078       typename r_binning = axis::regular<scalar_type, host_container_types>,
0079       typename phi_binning = axis::regular<scalar_type, host_container_types>>
0080     requires std::is_enum_v<decltype(r_bounds::label)>
0081   auto new_grid(const mask<annulus2D, algebra_type> &grid_bounds,
0082                 const darray<std::size_t, 2UL> n_bins,
0083                 const std::vector<std::pair<loc_bin_index<annulus2D>, dindex>>
0084                     &bin_capacities = {},
0085                 const darray<std::vector<scalar_type>, 2UL> &bin_edges =
0086                     darray<std::vector<scalar_type>, 2UL>(),
0087                 const darray<std::vector<scalar_type>, 2UL> &axis_spans =
0088                     darray<std::vector<scalar_type>, 2UL>()) const {
0089     static_assert(
0090         std::is_same_v<phi_bounds, axis::circular<>>,
0091         "Phi axis bounds need to be circular for stereo annulus shape");
0092
0093     // Axes boundaries and local indices
0094     using boundary = annulus2D::boundaries;
0095     using axes_t = axes<annulus2D>::template type<algebra_t>;
0096     using local_frame = typename axes_t::template frame<algebra_t>;
0097
0098     constexpr auto e_r_axis = static_cast<dindex>(axes_t::label0);
0099     constexpr auto e_phi_axis = static_cast<dindex>(axes_t::label1);
0100
0101     auto b_values = grid_bounds.values();
0102     // Overwrite the mask values if axis spans are provided
0103     if (!axis_spans[0UL].empty()) {
0104       assert(axis_spans[0UL].size() == 2UL);
0105       b_values[boundary::e_min_r] = axis_spans[0UL].at(0UL);
0106       b_values[boundary::e_max_r] = axis_spans[0UL].at(1UL);
0107     }
0108     scalar_type min_phi =
0109         b_values[boundary::e_average_phi] - b_values[boundary::e_min_phi_rel];
0110     scalar_type max_phi =
0111         b_values[boundary::e_average_phi] + b_values[boundary::e_max_phi_rel];
0112     if (!axis_spans[1UL].empty()) {
0113       assert(axis_spans[1UL].size() == 2UL);
0114       min_phi = axis_spans[1UL].at(0UL);
0115       max_phi = axis_spans[1UL].at(1UL);
0116     }
0117
0118     return new_grid<local_frame>(
0119         {b_values[boundary::e_min_r], b_values[boundary::e_max_r], min_phi,
0120          max_phi},
0121         {n_bins[e_r_axis], n_bins[e_phi_axis]}, bin_capacities,
0122         {bin_edges[e_r_axis], bin_edges[e_phi_axis]},
0123         types::list<r_bounds, phi_bounds>{},
0124         types::list<r_binning, phi_binning>{});
0125   }
0126
0127   //
0128   // cuboid 3D
0129   //
0130   template <
0131       typename x_bounds = axis::closed<axis::label::e_x>,
0132       typename y_bounds = axis::closed<axis::label::e_y>,
0133       typename z_bounds = axis::closed<axis::label::e_z>,
0134       typename x_binning = axis::regular<scalar_type, host_container_types>,
0135       typename y_binning = axis::regular<scalar_type, host_container_types>,
0136       typename z_binning = axis::regular<scalar_type, host_container_types>>
0137     requires std::is_enum_v<decltype(x_bounds::label)>
0138   auto new_grid(const mask<cuboid3D, algebra_type> &grid_bounds,
0139                 const darray<std::size_t, 3UL> n_bins,
0140                 const std::vector<std::pair<loc_bin_index<cuboid3D>, dindex>>
0141                     &bin_capacities = {},
0142                 const darray<std::vector<scalar_type>, 3UL> &bin_edges =
0143                     darray<std::vector<scalar_type>, 3UL>(),
0144                 const darray<std::vector<scalar_type>, 3UL> &axis_spans =
0145                     darray<std::vector<scalar_type>, 3UL>()) const {
0146     // Axes boundaries and local indices
0147     using boundary = cuboid3D::boundaries;
0148     using axes_t = axes<cuboid3D>::template type<algebra_t>;
0149     using local_frame = typename axes_t::template frame<algebra_t>;
0150
0151     constexpr auto e_x_axis = static_cast<dindex>(axes_t::label0);
0152     constexpr auto e_y_axis = static_cast<dindex>(axes_t::label1);
0153     constexpr auto e_z_axis = static_cast<dindex>(axes_t::label2);
0154
0155     // Overwrite the mask values if axis spans are provided
0156     auto b_values = grid_bounds.values();
0157     if (!axis_spans[0UL].empty()) {
0158       assert(axis_spans[0UL].size() == 2UL);
0159       b_values[boundary::e_min_x] = axis_spans[0UL].at(0UL);
0160       b_values[boundary::e_max_x] = axis_spans[0UL].at(1UL);
0161     }
0162     if (!axis_spans[1UL].empty()) {
0163       assert(axis_spans[1UL].size() == 2UL);
0164       b_values[boundary::e_min_y] = axis_spans[1UL].at(0UL);
0165       b_values[boundary::e_max_y] = axis_spans[1UL].at(1UL);
0166     }
0167     if (!axis_spans[2UL].empty()) {
0168       assert(axis_spans[2UL].size() == 2UL);
0169       b_values[boundary::e_min_z] = axis_spans[2UL].at(0UL);
0170       b_values[boundary::e_max_z] = axis_spans[2UL].at(1UL);
0171     }
0172
0173     return new_grid<local_frame>(
0174         {b_values[boundary::e_min_x], b_values[boundary::e_max_x],
0175          b_values[boundary::e_min_y], b_values[boundary::e_max_y],
0176          b_values[boundary::e_min_z], b_values[boundary::e_max_z]},
0177         {n_bins[e_x_axis], n_bins[e_y_axis], n_bins[e_z_axis]}, bin_capacities,
0178         {bin_edges[e_x_axis], bin_edges[e_y_axis], bin_edges[e_z_axis]},
0179         types::list<x_bounds, y_bounds, z_bounds>{},
0180         types::list<x_binning, y_binning, z_binning>{});
0181   }
0182
0183   //
0184   // cylinder 2D
0185   //
0186   template <
0187       typename rphi_bounds = axis::circular<axis::label::e_rphi>,
0188       typename z_bounds = axis::closed<axis::label::e_cyl_z>,
0189       typename rphi_binning = axis::regular<scalar_type, host_container_types>,
0190       typename z_binning = axis::regular<scalar_type, host_container_types>>
0191     requires std::is_enum_v<decltype(rphi_bounds::label)>
0192   auto new_grid(const mask<cylinder2D, algebra_type> &grid_bounds,
0193                 const darray<std::size_t, 2UL> n_bins,
0194                 const std::vector<std::pair<loc_bin_index<cylinder2D>, dindex>>
0195                     &bin_capacities = {},
0196                 const darray<std::vector<scalar_type>, 2UL> &bin_edges =
0197                     darray<std::vector<scalar_type>, 2UL>(),
0198                 const darray<std::vector<scalar_type>, 2UL> &axis_spans =
0199                     darray<std::vector<scalar_type>, 2UL>()) const {
0200     static_assert(
0201         std::is_same_v<rphi_bounds, axis::circular<axis::label::e_rphi>>,
0202         "Phi axis bounds need to be circular for cylinder2D shape");
0203
0204     // Axes boundaries and local indices
0205     using boundary = cylinder2D::boundaries;
0206     using axes_t = axes<cylinder2D>::template type<algebra_t>;
0207     using local_frame = typename axes_t::template frame<algebra_t>;
0208
0209     constexpr auto e_rphi_axis = static_cast<dindex>(axes_t::label0);
0210     constexpr auto e_z_axis = static_cast<dindex>(axes_t::label1);
0211
0212     auto b_values = grid_bounds.values();
0213     if (!axis_spans[1UL].empty()) {
0214       assert(axis_spans[1UL].size() == 2UL);
0215       b_values[boundary::e_lower_z] = axis_spans[1UL].at(0UL);
0216       b_values[boundary::e_upper_z] = axis_spans[1UL].at(1UL);
0217     }
0218     return new_grid<local_frame>(
0219         {-constant<scalar_type>::pi, constant<scalar_type>::pi,
0220          b_values[boundary::e_lower_z], b_values[boundary::e_upper_z]},
0221         {n_bins[e_rphi_axis], n_bins[e_z_axis]}, bin_capacities,
0222         {bin_edges[e_rphi_axis], bin_edges[e_z_axis]},
0223         types::list<rphi_bounds, z_bounds>{},
0224         types::list<rphi_binning, z_binning>{});
0225   }
0226
0227   //
0228   // concentric cylinder 2D
0229   //
0230   template <
0231       typename rphi_bounds = axis::circular<axis::label::e_rphi>,
0232       typename z_bounds = axis::closed<axis::label::e_cyl_z>,
0233       typename rphi_binning = axis::regular<scalar_type, host_container_types>,
0234       typename z_binning = axis::regular<scalar_type, host_container_types>>
0235     requires std::is_enum_v<decltype(rphi_bounds::label)>
0236   auto new_grid(
0237       const mask<concentric_cylinder2D, algebra_type> &grid_bounds,
0238       const darray<std::size_t, 2UL> n_bins,
0239       const std::vector<std::pair<loc_bin_index<concentric_cylinder2D>, dindex>>
0240           &bin_capacities = {},
0241       const darray<std::vector<scalar_type>, 2UL> &bin_edges =
0242           darray<std::vector<scalar_type>, 2UL>(),
0243       const darray<std::vector<scalar_type>, 2UL> &axis_spans =
0244           darray<std::vector<scalar_type>, 2UL>()) const {
0245     static_assert(
0246         std::is_same_v<rphi_bounds, axis::circular<axis::label::e_rphi>>,
0247         "Phi axis bounds need to be circular for cylinder2D portal shape");
0248
0249     // Axes boundaries and local indices
0250     using boundary = concentric_cylinder2D::boundaries;
0251     using axes_t = axes<concentric_cylinder2D>::template type<algebra_t>;
0252     using local_frame = typename axes_t::template frame<algebra_t>;
0253
0254     constexpr auto e_rphi_axis = static_cast<dindex>(axes_t::label0);
0255     constexpr auto e_z_axis = static_cast<dindex>(axes_t::label1);
0256
0257     auto b_values = grid_bounds.values();
0258     if (!axis_spans[1UL].empty()) {
0259       assert(axis_spans[1UL].size() == 2UL);
0260       b_values[boundary::e_lower_z] = axis_spans[1UL].at(0UL);
0261       b_values[boundary::e_upper_z] = axis_spans[1UL].at(1UL);
0262     }
0263
0264     return new_grid<local_frame>(
0265         {-constant<scalar_type>::pi, constant<scalar_type>::pi,
0266          b_values[boundary::e_lower_z], b_values[boundary::e_upper_z]},
0267         {n_bins[e_rphi_axis], n_bins[e_z_axis]}, bin_capacities,
0268         {bin_edges[e_rphi_axis], bin_edges[e_z_axis]},
0269         types::list<rphi_bounds, z_bounds>{},
0270         types::list<rphi_binning, z_binning>{});
0271   }
0272
0273   //
0274   // cylinder 3D
0275   //
0276   template <
0277       typename r_bounds = axis::closed<axis::label::e_r>,
0278       typename phi_bounds = axis::circular<>,
0279       typename z_bounds = axis::closed<axis::label::e_z>,
0280       typename r_binning = axis::regular<scalar_type, host_container_types>,
0281       typename phi_binning = axis::regular<scalar_type, host_container_types>,
0282       typename z_binning = axis::regular<scalar_type, host_container_types>>
0283     requires std::is_enum_v<decltype(r_bounds::label)>
0284   auto new_grid(const mask<cylinder3D, algebra_type> &grid_bounds,
0285                 const darray<std::size_t, 3UL> n_bins,
0286                 const std::vector<std::pair<loc_bin_index<cylinder3D>, dindex>>
0287                     &bin_capacities = {},
0288                 const darray<std::vector<scalar_type>, 3UL> &bin_edges =
0289                     darray<std::vector<scalar_type>, 3UL>(),
0290                 const darray<std::vector<scalar_type>, 3UL> &axis_spans =
0291                     darray<std::vector<scalar_type>, 3UL>()) const {
0292     static_assert(std::is_same_v<phi_bounds, axis::circular<>>,
0293                   "Phi axis bounds need to be circular for cylinder3D shape");
0294
0295     // Axes boundaries and local indices
0296     using boundary = cylinder3D::boundaries;
0297     using axes_t = axes<cylinder3D>::template type<algebra_t>;
0298     using local_frame = typename axes_t::template frame<algebra_t>;
0299
0300     constexpr auto e_r_axis = static_cast<dindex>(axes_t::label0);
0301     constexpr auto e_phi_axis = static_cast<dindex>(axes_t::label1);
0302     constexpr auto e_z_axis = static_cast<dindex>(axes_t::label2);
0303
0304     auto b_values = grid_bounds.values();
0305
0306     // Overwrite the mask values if axis spans are provided
0307     if (!axis_spans[0UL].empty()) {
0308       assert(axis_spans[0UL].size() == 2UL);
0309       b_values[boundary::e_min_r] = axis_spans[0UL].at(0UL);
0310       b_values[boundary::e_max_r] = axis_spans[0UL].at(1UL);
0311     }
0312     scalar_type min_phi = -constant<scalar_type>::pi;
0313     scalar_type max_phi = constant<scalar_type>::pi;
0314     if (!axis_spans[1UL].empty()) {
0315       assert(axis_spans[1UL].size() == 2UL);
0316       min_phi = axis_spans[1UL].at(0UL);
0317       max_phi = axis_spans[1UL].at(1UL);
0318     }
0319     if (!axis_spans[2UL].empty()) {
0320       assert(axis_spans[2UL].size() == 2UL);
0321       b_values[boundary::e_min_z] = axis_spans[2UL].at(0UL);
0322       b_values[boundary::e_max_z] = axis_spans[2UL].at(1UL);
0323     }
0324
0325     return new_grid<local_frame>(
0326         {b_values[boundary::e_min_r], b_values[boundary::e_max_r], min_phi,
0327          max_phi, -b_values[boundary::e_min_z], b_values[boundary::e_max_z]},
0328         {n_bins[e_r_axis], n_bins[e_phi_axis], n_bins[e_z_axis]},
0329         bin_capacities,
0330         {bin_edges[e_r_axis], bin_edges[e_phi_axis], bin_edges[e_z_axis]},
0331         types::list<r_bounds, phi_bounds, z_bounds>{},
0332         types::list<r_binning, phi_binning, z_binning>{});
0333   }
0334
0335   //
0336   // polar 2D
0337   //
0338   template <
0339       typename r_bounds = axis::closed<axis::label::e_r>,
0340       typename phi_bounds = axis::circular<>,
0341       typename r_binning = axis::regular<scalar_type, host_container_types>,
0342       typename phi_binning = axis::regular<scalar_type, host_container_types>>
0343     requires std::is_enum_v<decltype(r_bounds::label)>
0344   auto new_grid(const mask<ring2D, algebra_type> &grid_bounds,
0345                 const darray<std::size_t, 2UL> n_bins,
0346                 const std::vector<std::pair<loc_bin_index<ring2D>, dindex>>
0347                     &bin_capacities = {},
0348                 const darray<std::vector<scalar_type>, 2UL> &bin_edges =
0349                     darray<std::vector<scalar_type>, 2UL>(),
0350                 const darray<std::vector<scalar_type>, 2UL> &axis_spans =
0351                     darray<std::vector<scalar_type>, 2UL>()) const {
0352     static_assert(std::is_same_v<phi_bounds, axis::circular<>>,
0353                   "Phi axis bounds need to be circular for ring shape");
0354
0355     // Axes boundaries and local indices
0356     using boundary = ring2D::boundaries;
0357     using axes_t = axes<ring2D>::template type<algebra_t>;
0358     using local_frame = typename axes_t::template frame<algebra_t>;
0359
0360     constexpr auto e_r_axis = static_cast<dindex>(axes_t::label0);
0361     constexpr auto e_phi_axis = static_cast<dindex>(axes_t::label1);
0362
0363     auto b_values = grid_bounds.values();
0364     // Overwrite the mask values if axis spans are provided
0365     if (!axis_spans[0UL].empty()) {
0366       assert(axis_spans[0UL].size() == 2UL);
0367       b_values[boundary::e_inner_r] = axis_spans[0UL].at(0UL);
0368       b_values[boundary::e_outer_r] = axis_spans[0UL].at(1UL);
0369     }
0370
0371     return new_grid<local_frame>(
0372         {b_values[boundary::e_inner_r], b_values[boundary::e_outer_r],
0373          -constant<scalar_type>::pi, constant<scalar_type>::pi},
0374         {n_bins[e_r_axis], n_bins[e_phi_axis]}, bin_capacities,
0375         {bin_edges[e_r_axis], bin_edges[e_phi_axis]},
0376         types::list<r_bounds, phi_bounds>{},
0377         types::list<r_binning, phi_binning>{});
0378   }
0379
0380   //
0381   // rectangle 2D
0382   //
0383   template <
0384       typename x_bounds = axis::closed<axis::label::e_x>,
0385       typename y_bounds = axis::closed<axis::label::e_y>,
0386       typename x_binning = axis::regular<scalar_type, host_container_types>,
0387       typename y_binning = axis::regular<scalar_type, host_container_types>>
0388     requires std::is_enum_v<decltype(x_bounds::label)>
0389   auto new_grid(const mask<rectangle2D, algebra_type> &grid_bounds,
0390                 const darray<std::size_t, 2UL> n_bins,
0391                 const std::vector<std::pair<loc_bin_index<rectangle2D>, dindex>>
0392                     &bin_capacities = {},
0393                 const darray<std::vector<scalar_type>, 2UL> &bin_edges =
0394                     darray<std::vector<scalar_type>, 2UL>(),
0395                 const darray<std::vector<scalar_type>, 2UL> &axis_spans =
0396                     darray<std::vector<scalar_type>, 2UL>()) const {
0397     // Axes boundaries and local indices
0398     using boundary = rectangle2D::boundaries;
0399     using axes_t = axes<rectangle2D>::template type<algebra_t>;
0400     using local_frame = typename axes_t::template frame<algebra_t>;
0401
0402     constexpr auto e_x_axis = static_cast<dindex>(axes_t::label0);
0403     constexpr auto e_y_axis = static_cast<dindex>(axes_t::label1);
0404
0405     auto b_values = grid_bounds.values();
0406     // Overwrite the mask values if axis spans are provided
0407     if (!axis_spans[0UL].empty()) {
0408       assert(axis_spans[0UL].size() == 2UL);
0409       b_values[boundary::e_half_x] = axis_spans[0UL].at(1UL);
0410     }
0411     if (!axis_spans[1UL].empty()) {
0412       assert(axis_spans[1UL].size() == 2UL);
0413       b_values[boundary::e_half_y] = axis_spans[1UL].at(1UL);
0414     }
0415
0416     return new_grid<local_frame>(
0417         {-b_values[boundary::e_half_x], b_values[boundary::e_half_x],
0418          -b_values[boundary::e_half_y], b_values[boundary::e_half_y]},
0419         {n_bins[e_x_axis], n_bins[e_y_axis]}, bin_capacities,
0420         {bin_edges[e_x_axis], bin_edges[e_y_axis]},
0421         types::list<x_bounds, y_bounds>{}, types::list<x_binning, y_binning>{});
0422   }
0423
0424   //
0425   // trapezoid 2D
0426   //
0427   template <
0428       typename x_bounds = axis::closed<axis::label::e_x>,
0429       typename y_bounds = axis::closed<axis::label::e_y>,
0430       typename x_binning = axis::regular<scalar_type, host_container_types>,
0431       typename y_binning = axis::regular<scalar_type, host_container_types>>
0432     requires std::is_enum_v<decltype(x_bounds::label)>
0433   auto new_grid(const mask<trapezoid2D, algebra_type> &grid_bounds,
0434                 const darray<std::size_t, 2UL> n_bins,
0435                 const std::vector<std::pair<loc_bin_index<trapezoid2D>, dindex>>
0436                     &bin_capacities = {},
0437                 const darray<std::vector<scalar_type>, 2UL> &bin_edges =
0438                     darray<std::vector<scalar_type>, 2UL>(),
0439                 const darray<std::vector<scalar_type>, 2UL> &axis_spans =
0440                     darray<std::vector<scalar_type>, 2UL>()) const {
0441     // Axes boundaries and local indices
0442     using boundary = trapezoid2D::boundaries;
0443     using axes_t = axes<trapezoid2D>::template type<algebra_t>;
0444     using local_frame = typename axes_t::template frame<algebra_t>;
0445
0446     constexpr auto e_x_axis = static_cast<dindex>(axes_t::label0);
0447     constexpr auto e_y_axis = static_cast<dindex>(axes_t::label1);
0448
0449     auto b_values = grid_bounds.values();
0450     // Overwrite the mask values if axis spans are provided
0451     if (!axis_spans[0UL].empty()) {
0452       assert(axis_spans[0UL].size() == 2UL);
0453       b_values[boundary::e_half_length_1] = axis_spans[0UL].at(1UL);
0454     }
0455     if (!axis_spans[1UL].empty()) {
0456       assert(axis_spans[1UL].size() == 2UL);
0457       b_values[boundary::e_half_length_2] = axis_spans[1UL].at(1UL);
0458     }
0459
0460     return new_grid<local_frame>(
0461         {-b_values[boundary::e_half_length_1],
0462          b_values[boundary::e_half_length_1],
0463          -b_values[boundary::e_half_length_2],
0464          b_values[boundary::e_half_length_2]},
0465         {n_bins[e_x_axis], n_bins[e_y_axis]}, bin_capacities,
0466         {bin_edges[e_x_axis], bin_edges[e_y_axis]},
0467         types::list<x_bounds, y_bounds>{}, types::list<x_binning, y_binning>{});
0468   }
0469
0470   /// @brief Create and empty grid with fully initialized axes.
0471   ///
0472   /// @tparam grid_shape_t the shape of the resulting grid
0473   ///         (e.g. cylinder2D).
0474   /// @tparam e_bounds the bounds of the regular axes
0475   ///         (open vs. closed bounds).
0476   /// @tparam binnings the binning types of the axes
0477   ///         (regular vs. irregular)
0478   ///
0479   /// @param spans the span of the axis values for regular axes, otherwise
0480   ///              ignored.
0481   /// @param n_bins the number of bins for regular axes, otherwise ignored
0482   /// @param ax_bin_edges the explicit bin edges for irregular axes
0483   ///                     (lower bin edges + the the upper edge of the
0484   ///                     last bin), otherwise ignored.
0485   template <concepts::coordinate_frame grid_frame_t, typename... bound_ts,
0486             typename... binning_ts>
0487   auto new_grid(
0488       const std::vector<scalar_type> &spans,
0489       const std::vector<std::size_t> &n_bins,
0490       const std::vector<std::pair<axis::multi_bin<sizeof...(bound_ts)>, dindex>>
0491           &bin_capacities = {},
0492       const std::vector<std::vector<scalar_type>> &ax_bin_edges = {},
0493       const types::list<bound_ts...> & /*unused*/ = {},
0494       const types::list<binning_ts...> & /*unused*/ = {}) const {
0495     static_assert(sizeof...(bound_ts) == sizeof...(binning_ts),
0496                   "number of axis bounds and binning types has to match");
0497
0498     // Build the coordinate axes and the grid
0499     using axes_t = axis::multi_axis<is_owning, grid_frame_t,
0500                                     axis::single_axis<bound_ts, binning_ts>...>;
0501     using grid_t = grid<algebra_t, axes_t, bin_t, serializer_t>;
0502
0503     return new_grid<grid_t>(spans, n_bins, bin_capacities, ax_bin_edges);
0504   }
0505
0506   /// Helper to build grid from shape plus binning and bounds types
0507   template <typename grid_shape_t, typename... bound_ts, typename... binning_ts>
0508     requires concepts::shape<grid_shape_t, algebra_t>
0509   auto new_grid(
0510       const std::vector<scalar_type> &spans,
0511       const std::vector<std::size_t> &n_bins,
0512       const std::vector<std::pair<axis::multi_bin<sizeof...(bound_ts)>, dindex>>
0513           &bin_capacities = {},
0514       const std::vector<std::vector<scalar_type>> &ax_bin_edges = {},
0515       const types::list<bound_ts...> &bounds = {},
0516       const types::list<binning_ts...> &binnings = {}) const {
0517     return new_grid<
0518         typename grid_shape_t::template local_frame_type<algebra_t>>(
0519         spans, n_bins, bin_capacities, ax_bin_edges, bounds, binnings);
0520   }
0521
0522   /// Helper overload for grid builder: Build from mask and resolve bounds
0523   /// and binnings from concrete grid type
0524   template <concepts::grid grid_t, typename grid_shape_t>
0525     requires concepts::shape<grid_shape_t, algebra_t>
0526   auto new_grid(
0527       const mask<grid_shape_t, algebra_type> &m,
0528       const darray<std::size_t, grid_t::dim> &n_bins,
0529       const std::vector<std::pair<typename grid_t::loc_bin_index, dindex>>
0530           &bin_capacities = {},
0531       const darray<std::vector<scalar_type>, grid_t::dim> &ax_bin_edges = {})
0532       const {
0533     return new_grid(m, n_bins, bin_capacities, ax_bin_edges,
0534                     typename grid_t::axes_type::bounds{},
0535                     typename grid_t::axes_type::binnings{});
0536   }
0537
0538   /// Helper overload for grid builder: Build from mask and resolve bounds
0539   /// and binnings
0540   template <typename grid_shape_t, typename... bound_ts, typename... binning_ts>
0541     requires concepts::shape<grid_shape_t, algebra_t>
0542   auto new_grid(
0543       const mask<grid_shape_t, algebra_type> &m,
0544       const darray<std::size_t, grid_shape_t::dim> &n_bins,
0545       const std::vector<std::pair<axis::multi_bin<sizeof...(bound_ts)>, dindex>>
0546           &bin_capacities = {},
0547       const darray<std::vector<scalar_type>, grid_shape_t::dim> &ax_bin_edges =
0548           {},
0549       const types::list<bound_ts...> & /*unused*/ = {},
0550       const types::list<binning_ts...> & /*unused*/ = {}) const {
0551     return new_grid<bound_ts..., binning_ts...>(m, n_bins, bin_capacities,
0552                                                 ax_bin_edges);
0553   }
0554
0555   /// @brief Create and empty grid with fully initialized axes.
0556   ///
0557   /// @tparam grid_t the type of the resulting grid
0558   ///
0559   /// @param spans the span of the axis values for regular axes, otherwise
0560   ///              ignored.
0561   /// @param n_bins the number of bins for regular axes, otherwise ignored
0562   /// @param ax_bin_edges the explicit bin edges for irregular axes
0563   ///                     (lower bin edges + the the upper edge of the
0564   ///                     last bin), otherwise ignored.
0565   template <concepts::grid grid_t>
0566   auto new_grid(
0567       const std::vector<scalar_type> &spans,
0568       const std::vector<std::size_t> &n_bins,
0569       [[maybe_unused]] const std::vector<
0570           std::pair<typename grid_t::loc_bin_index, dindex>> &bin_capacities =
0571           {},
0572       const std::vector<std::vector<scalar_type>> &ax_bin_edges = {}) const {
0573     using owning_grid_t = typename grid_t::template type<true>;
0574     using axes_t = typename owning_grid_t::axes_type;
0575
0576     // Prepare data
0577     vector_type<dsized_index_range> axes_data{};
0578     vector_type<scalar_type> bin_edges{};
0579
0580     // Call init for every axis
0581     unroll_axis_init<typename axes_t::binnings>(
0582         spans, n_bins, ax_bin_edges, axes_data, bin_edges,
0583         std::make_index_sequence<axes_t::dim>{});
0584
0585     // Assemble the grid and return it
0586     axes_t axes(std::move(axes_data), std::move(bin_edges));
0587
0588     typename owning_grid_t::bin_container_type bin_data{};
0589
0590     // Bins with dynamic capacity and "index grids" need different treatment
0591     if constexpr (std::is_same_v<bin_t, bins::dynamic_array<
0592                                             typename grid_t::value_type>>) {
0593       // Bin and bin entries vector are separate containers
0594       bin_data.bins.resize(axes.nbins());
0595       // Set correct bin capacities, if they were provided
0596       if (!bin_capacities.empty()) {
0597         // Set memory offsets correctly
0598         typename grid_t::template serializer_type<grid_t::dim> serializer{};
0599         dindex total_cap{0u};
0600         for (const auto &capacity : bin_capacities) {
0601           // Get the empty bin by its global index
0602           auto &data = bin_data.bins[serializer(axes, capacity.first)];
0603
0604           data.offset = total_cap;
0605           data.capacity = capacity.second;
0606           data.size = 0u;
0607
0608           total_cap += data.capacity;
0609         }
0610         assert(total_cap > 0);
0611         bin_data.entries.resize(
0612             total_cap, detail::invalid_value<typename grid_t::value_type>());
0613       }
0614     } else {
0615       // Bin entries are contained in the bins directly
0616       bin_data.resize(axes.nbins(), bin_type{});
0617     }
0618
0619     return owning_grid_t(std::move(bin_data), std::move(axes));
0620   }
0621
0622  private:
0623   /// Initialize a single axis (either regular or irregular)
0624   /// @note change to template lambda as soon as it becomes available.
0625   template <std::size_t I, typename binnings>
0626   auto axis_init([[maybe_unused]] const std::vector<scalar_type> &spans,
0627                  [[maybe_unused]] const std::vector<std::size_t> &n_bins,
0628                  [[maybe_unused]] const std::vector<std::vector<scalar_type>>
0629                      &ax_bin_edges,
0630                  vector_type<dsized_index_range> &axes_data,
0631                  vector_type<scalar_type> &bin_edges) const {
0632     if constexpr (std::is_same_v<
0633                       types::at<binnings, I>,
0634                       axis::regular<scalar_type, host_container_types>>) {
0635       axes_data.push_back({static_cast<dindex>(bin_edges.size()),
0636                            static_cast<dindex>(n_bins.at(I))});
0637       bin_edges.push_back(spans.at(I * 2u));
0638       bin_edges.push_back(spans.at(I * 2u + 1u));
0639     } else {
0640       const auto &bin_edges_loc = ax_bin_edges.at(I);
0641       axes_data.push_back({static_cast<dindex>(bin_edges.size()),
0642                            static_cast<dindex>(bin_edges_loc.size() - 1u)});
0643       bin_edges.insert(bin_edges.end(), bin_edges_loc.begin(),
0644                        bin_edges_loc.end());
0645     }
0646   }
0647
0648   /// Call axis init for every dimension
0649   template <typename binnings, std::size_t... I>
0650   auto unroll_axis_init(
0651       const std::vector<scalar_type> &spans,
0652       const std::vector<std::size_t> &n_bins,
0653       const std::vector<std::vector<scalar_type>> &ax_bin_edges,
0654       vector_type<dsized_index_range> &axes_data,
0655       vector_type<scalar_type> &bin_edges,
0656       std::index_sequence<I...> /*ids*/) const {
0657     (axis_init<I, binnings>(spans, n_bins, ax_bin_edges, axes_data, bin_edges),
0658      ...);
0659   }
0660
0661   vecmem::memory_resource *m_resource{};
0662 };
0663
0664 // Infer a grid factory type from an already completely assembled grid type
0665 template <concepts::grid grid_t>
0666 using grid_factory_type =
0667     grid_factory<typename grid_t::bin_type,
0668                  simple_serializer /*grid_t::template serializer_type*/,
0669                  typename grid_t::algebra_type>;
0670
0671 }  // namespace detray