builders/detail/bin_association.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/detail/associator.hpp"
0013 #include "detray/definitions/algebra.hpp"
0014 #include "detray/definitions/units.hpp"
0015 #include "detray/geometry/coordinates/concentric_cylindrical2D.hpp"
0016 #include "detray/geometry/coordinates/cylindrical2D.hpp"
0017 #include "detray/geometry/coordinates/polar2D.hpp"
0018 #include "detray/geometry/detail/vertexer.hpp"
0019 #include "detray/navigation/accelerators/concepts.hpp"
0020 #include "detray/utils/grid/populators.hpp"
0021 #include "detray/utils/ranges.hpp"
0022
0023 // System include(s)
0024 #include <vector>
0025
0026 namespace detray::detail {
0027
0028 /// Run the bin association of surfaces (via their contour) to a given 2D grid.
0029 ///
0030 /// @param context is the context to win which the association is done
0031 /// @param surfaces a range of detector surfaces
0032 /// @param transforms the transforms that belong to the surfaces
0033 /// @param surface_masks the masks that belong to the surfaces
0034 /// @param grid either a cylinder or disc grid to be filled
0035 /// @param tolerance is the bin_tolerance in the two local coordinates
0036 /// @param absolute_tolerance is an indicator if the tolerance is to be
0037 ///        taken absolute or relative
0038 template <typename context_t, typename surface_container_t,
0039           typename transform_container_t, typename mask_container_t,
0040           concepts::surface_grid grid_t, concepts::scalar scalar_t>
0041 static inline void bin_association(const context_t & /*context*/,
0042                                    const surface_container_t &surfaces,
0043                                    const transform_container_t &transforms,
0044                                    const mask_container_t &surface_masks,
0045                                    grid_t &grid,
0046                                    const darray<scalar_t, 2> &bin_tolerance,
0047                                    bool absolute_tolerance = true) {
0048   using algebra_t = typename grid_t::local_frame_type::algebra_type;
0049   using point2_t = dpoint2D<algebra_t>;
0050   using point3_t = dpoint3D<algebra_t>;
0051
0052   const auto &axis_0 = grid.template get_axis<0>();
0053   const auto &axis_1 = grid.template get_axis<1>();
0054
0055   // Disk type bin association
0056   if constexpr (std::is_same_v<typename grid_t::local_frame_type,
0057                                polar2D<algebra_t>>) {
0058     // Run with two different associators: center of gravity and edge
0059     // intersection
0060     center_of_gravity_generic<algebra_t> cgs_assoc;
0061     edges_intersect_generic<algebra_t> edges_assoc;
0062
0063     // Loop over all bins and associate the surfaces
0064     for (unsigned int bin_0 = 0u; bin_0 < axis_0.nbins(); ++bin_0) {
0065       for (unsigned int bin_1 = 0u; bin_1 < axis_1.nbins(); ++bin_1) {
0066         auto r_borders = axis_0.bin_edges(bin_0);
0067         auto phi_borders = axis_1.bin_edges(bin_1);
0068
0069         scalar_t r_add = absolute_tolerance
0070                              ? bin_tolerance[0]
0071                              : bin_tolerance[0] * (r_borders[1] - r_borders[0]);
0072         scalar_t phi_add =
0073             absolute_tolerance
0074                 ? bin_tolerance[1]
0075                 : bin_tolerance[1] * (phi_borders[1] - phi_borders[0]);
0076
0077         // Create a contour for the bin
0078         std::vector<point2_t> bin_contour =
0079             detail::r_phi_polygon<scalar_t, point2_t>(
0080                 r_borders[0] - r_add, r_borders[1] + r_add,
0081                 phi_borders[0] - phi_add, phi_borders[1] + phi_add);
0082
0083         // Run through the surfaces and associate them by contour
0084         for (auto sf : surfaces) {
0085           // Add only sensitive surfaces to the grid
0086           if (sf.is_portal()) {
0087             continue;
0088           }
0089
0090           // Unroll the mask container and generate vertices
0091           const auto &transform = transforms[sf.transform()];
0092
0093           auto vertices_per_masks =
0094               surface_masks
0095                   .template visit<detail::vertexer<point2_t, point3_t>>(
0096                       sf.mask());
0097
0098           // Usually one mask per surface, but design allows - a
0099           // single association  is sufficient though
0100           for (auto &vertices : vertices_per_masks) {
0101             if (vertices.empty()) {
0102               continue;
0103             }
0104             // Create a surface contour
0105             std::vector<point2_t> surface_contour;
0106             surface_contour.reserve(vertices.size());
0107             for (const auto &v : vertices) {
0108               auto vg = transform.point_to_global(v);
0109               surface_contour.push_back({vg[0], vg[1]});
0110             }
0111             // The association has worked
0112             if (cgs_assoc(bin_contour, surface_contour) ||
0113                 edges_assoc(bin_contour, surface_contour)) {
0114               grid.template populate<attach<>>({bin_0, bin_1}, sf);
0115               break;
0116             }
0117           }
0118         }
0119       }
0120     }
0121   } else if constexpr (std::is_same_v<typename grid_t::local_frame_type,
0122                                       cylindrical2D<algebra_t>> ||
0123                        std::is_same_v<typename grid_t::local_frame_type,
0124                                       concentric_cylindrical2D<algebra_t>>) {
0125     center_of_gravity_rectangle<algebra_t> cgs_assoc;
0126     edges_intersect_generic<algebra_t> edges_assoc;
0127
0128     // Loop over all bins and associate the surfaces
0129     for (unsigned int bin_0 = 0u; bin_0 < axis_0.nbins(); ++bin_0) {
0130       for (unsigned int bin_1 = 0u; bin_1 < axis_1.nbins(); ++bin_1) {
0131         auto z_borders = axis_0.bin_edges(bin_0);
0132         auto phi_borders = axis_1.bin_edges(bin_1);
0133
0134         scalar_t z_add = absolute_tolerance
0135                              ? bin_tolerance[0]
0136                              : bin_tolerance[0] * (z_borders[1] - z_borders[0]);
0137         scalar_t phi_add =
0138             absolute_tolerance
0139                 ? bin_tolerance[1]
0140                 : bin_tolerance[1] * (phi_borders[1] - phi_borders[0]);
0141
0142         scalar_t z_min = z_borders[0];
0143         scalar_t z_max = z_borders[1];
0144         scalar_t phi_min_rep = phi_borders[0];
0145         scalar_t phi_max_rep = phi_borders[1];
0146
0147         point2_t p0_bin = {z_min - z_add, phi_min_rep - phi_add};
0148         point2_t p1_bin = {z_min - z_add, phi_max_rep + phi_add};
0149         point2_t p2_bin = {z_max + z_add, phi_max_rep + phi_add};
0150         point2_t p3_bin = {z_max + z_add, phi_min_rep - phi_add};
0151
0152         std::vector<point2_t> bin_contour = {p0_bin, p1_bin, p2_bin, p3_bin};
0153
0154         // Loop over the surfaces within a volume
0155         for (auto sf : surfaces) {
0156           // Add only sensitive surfaces to the grid
0157           if (sf.is_portal()) {
0158             continue;
0159           }
0160
0161           // Unroll the mask container and generate vertices
0162           const auto &transform = transforms.at(sf.transform());
0163
0164           auto vertices_per_masks =
0165               surface_masks
0166                   .template visit<detail::vertexer<point2_t, point3_t>>(
0167                       sf.mask());
0168
0169           for (auto &vertices : vertices_per_masks) {
0170             if (!vertices.empty()) {
0171               // Create a surface contour
0172               std::vector<point2_t> surface_contour{};
0173               surface_contour.reserve(vertices.size());
0174               scalar_t phi_min = std::numeric_limits<scalar_t>::max();
0175               scalar_t phi_max = -std::numeric_limits<scalar_t>::max();
0176               // We potentially need the split vertices
0177               std::vector<point2_t> s_c_neg{};
0178               std::vector<point2_t> s_c_pos{};
0179               scalar_t z_min_neg = std::numeric_limits<scalar_t>::max();
0180               scalar_t z_max_neg = -std::numeric_limits<scalar_t>::max();
0181               scalar_t z_min_pos = std::numeric_limits<scalar_t>::max();
0182               scalar_t z_max_pos = -std::numeric_limits<scalar_t>::max();
0183
0184               for (const auto &v : vertices) {
0185                 const point3_t vg = transform.point_to_global(v);
0186                 scalar_t phi = math::atan2(vg[1], vg[0]);
0187                 phi_min = math::min(phi, phi_min);
0188                 phi_max = math::max(phi, phi_max);
0189                 surface_contour.push_back({vg[2], phi});
0190                 if (phi < 0.f) {
0191                   s_c_neg.push_back({vg[2], phi});
0192                   z_min_neg = math::min(vg[2], z_min_neg);
0193                   z_max_neg = math::max(vg[2], z_max_neg);
0194                 } else {
0195                   s_c_pos.push_back({vg[2], phi});
0196                   z_min_pos = math::min(vg[2], z_min_pos);
0197                   z_max_pos = math::max(vg[2], z_max_pos);
0198                 }
0199               }
0200               // Check for phi wrapping
0201               std::vector<std::vector<point2_t>> surface_contours{};
0202               if (phi_max - phi_min > constant<scalar_t>::pi &&
0203                   phi_max * phi_min < 0.f) {
0204                 s_c_neg.push_back({z_max_neg, -constant<scalar_t>::pi});
0205                 s_c_neg.push_back({z_min_neg, -constant<scalar_t>::pi});
0206                 s_c_pos.push_back({z_max_pos, constant<scalar_t>::pi});
0207                 s_c_pos.push_back({z_min_pos, constant<scalar_t>::pi});
0208                 surface_contours.insert(surface_contours.end(),
0209                                         {s_c_neg, s_c_pos});
0210               } else {
0211                 surface_contours.push_back(surface_contour);
0212               }
0213
0214               // Check the association (with potential splits)
0215               bool associated = false;
0216               for (const auto &s_c : surface_contours) {
0217                 if (cgs_assoc(bin_contour, s_c) ||
0218                     edges_assoc(bin_contour, s_c)) {
0219                   associated = true;
0220                   break;
0221                 }
0222               }
0223
0224               // Register if associated
0225               if (associated) {
0226                 typename grid_t::loc_bin_index mbin{bin_0, bin_1};
0227                 grid.template populate<attach<>>(mbin, sf);
0228                 break;
0229               }
0230             }
0231           }
0232         }
0233       }
0234     }
0235   }
0236 }
0237
0238 }  // namespace detray::detail