src/Detector/LayerStructureBuilder.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2023 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 http://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/Detector/LayerStructureBuilder.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Detector/ProtoBinning.hpp"
0013 #include "Acts/Detector/detail/IndexedSurfacesGenerator.hpp"
0014 #include "Acts/Detector/detail/ReferenceGenerators.hpp"
0015 #include "Acts/Detector/detail/SupportSurfacesHelper.hpp"
0016 #include "Acts/Geometry/Extent.hpp"
0017 #include "Acts/Geometry/Polyhedron.hpp"
0018 #include "Acts/Navigation/DetectorVolumeFinders.hpp"
0019 #include "Acts/Navigation/NavigationDelegates.hpp"
0020 #include "Acts/Surfaces/Surface.hpp"
0021 #include "Acts/Utilities/BinningData.hpp"
0022 #include "Acts/Utilities/Enumerate.hpp"
0023 #include "Acts/Utilities/Grid.hpp"
0024 #include "Acts/Utilities/GridAxisGenerators.hpp"
0025 #include "Acts/Utilities/detail/AxisFwd.hpp"
0026
0027 #include <cmath>
0028 #include <cstddef>
0029 #include <ostream>
0030 #include <set>
0031 #include <stdexcept>
0032 #include <utility>
0033
0034 namespace Acts::Experimental {
0035 class DetectorVolume;
0036 }  // namespace Acts::Experimental
0037
0038 namespace {
0039
0040 /// Check autorange for a given binning
0041 ///
0042 /// @param pBinning the proto binning
0043 /// @param extent the extent from which the range is taken
0044 ///
0045 void adaptBinningRange(std::vector<Acts::Experimental::ProtoBinning>& pBinning,
0046                        const Acts::Extent& extent) {
0047   for (auto& pb : pBinning) {
0048     // Starting values
0049     Acts::ActsScalar vmin = pb.edges.front();
0050     Acts::ActsScalar vmax = pb.edges.back();
0051     // Get the number of bins
0052     std::size_t nBins = pb.bins();
0053     // Check if extent overwrites that
0054     if (extent.constrains(pb.binValue)) {
0055       const auto& range = extent.range(pb.binValue);
0056       // Patch the edges values from the range
0057       vmin = range.min();
0058       vmax = range.max();
0059     }
0060     // Possibly update the edges
0061     if (pb.axisType == Acts::detail::AxisType::Equidistant) {
0062       Acts::ActsScalar binWidth = (vmax - vmin) / nBins;
0063       // Fill the edges
0064       pb.edges = {vmin};
0065       pb.edges.resize(nBins + 1);
0066       for (std::size_t ib = 0; ib <= nBins; ++ib) {
0067         pb.edges[ib] = vmin + ib * binWidth;
0068       }
0069     } else {
0070       pb.edges.front() = vmin;
0071       pb.edges.back() = vmax;
0072     }
0073   }
0074 }
0075
0076 /// Helper for 1-dimensional generators
0077 ///
0078 /// @tparam aType is the axis boundary type: closed or bound
0079 ///
0080 /// @param gctx the geometry context
0081 /// @param lSurfaces the surfaces of the layer
0082 /// @param assignToAll the indices assigned to all
0083 /// @param binning the binning struct
0084 ///
0085 /// @return a configured surface candidate updators
0086 template <Acts::detail::AxisBoundaryType aType>
0087 Acts::Experimental::SurfaceCandidatesUpdater createUpdater(
0088     const Acts::GeometryContext& gctx,
0089     std::vector<std::shared_ptr<Acts::Surface>> lSurfaces,
0090     std::vector<std::size_t> assignToAll,
0091     const Acts::Experimental::ProtoBinning& binning) {
0092   // The surface candidate updator & a generator for polyhedrons
0093   Acts::Experimental::SurfaceCandidatesUpdater sfCandidates;
0094   Acts::Experimental::detail::PolyhedronReferenceGenerator rGenerator;
0095   // Indexed Surface generator for this case
0096   Acts::Experimental::detail::IndexedSurfacesGenerator<
0097       decltype(lSurfaces), Acts::Experimental::IndexedSurfacesImpl>
0098       isg{std::move(lSurfaces),
0099           std::move(assignToAll),
0100           {binning.binValue},
0101           {binning.expansion}};
0102   if (binning.axisType == Acts::detail::AxisType::Equidistant) {
0103     // Equidistant
0104     Acts::GridAxisGenerators::Eq<aType> aGenerator{
0105         {binning.edges.front(), binning.edges.back()}, binning.bins()};
0106     sfCandidates = isg(gctx, aGenerator, rGenerator);
0107   } else {
0108     // Variable
0109     Acts::GridAxisGenerators::Var<aType> aGenerator{binning.edges};
0110     sfCandidates = isg(gctx, aGenerator, rGenerator);
0111   }
0112   return sfCandidates;
0113 }
0114
0115 /// Helper for 2-dimensional generators
0116 ///
0117 /// @tparam aType is the axis boundary type, axis a: closed or bound
0118 /// @tparam bType is the axis boundary type, axis b: closed or bound
0119 ///
0120 /// @param gctx the geometry context
0121 /// @param lSurfaces the surfaces of the layer
0122 /// @param assignToAll the indices assigned to all
0123 /// @param aBinning the binning struct of axis a
0124 /// @param bBinning the binning struct of axis b
0125 ///
0126 /// @return a configured surface candidate updators
0127 template <Acts::detail::AxisBoundaryType aType,
0128           Acts::detail::AxisBoundaryType bType>
0129 Acts::Experimental::SurfaceCandidatesUpdater createUpdater(
0130     const Acts::GeometryContext& gctx,
0131     const std::vector<std::shared_ptr<Acts::Surface>>& lSurfaces,
0132     const std::vector<std::size_t>& assignToAll,
0133     const Acts::Experimental::ProtoBinning& aBinning,
0134     const Acts::Experimental::ProtoBinning& bBinning) {
0135   // The surface candidate updator & a generator for polyhedrons
0136   Acts::Experimental::SurfaceCandidatesUpdater sfCandidates;
0137   Acts::Experimental::detail::PolyhedronReferenceGenerator rGenerator;
0138   // Indexed Surface generator for this case
0139   Acts::Experimental::detail::IndexedSurfacesGenerator<
0140       decltype(lSurfaces), Acts::Experimental::IndexedSurfacesImpl>
0141       isg{lSurfaces,
0142           assignToAll,
0143           {aBinning.binValue, bBinning.binValue},
0144           {aBinning.expansion, bBinning.expansion}};
0145   // Run through the cases
0146   if (aBinning.axisType == Acts::detail::AxisType::Equidistant &&
0147       bBinning.axisType == Acts::detail::AxisType::Equidistant) {
0148     // Equidistant-Equidistant
0149     Acts::GridAxisGenerators::EqEq<aType, bType> aGenerator{
0150         {aBinning.edges.front(), aBinning.edges.back()},
0151         aBinning.bins(),
0152         {bBinning.edges.front(), bBinning.edges.back()},
0153         bBinning.bins()};
0154     sfCandidates = isg(gctx, aGenerator, rGenerator);
0155   } else if (bBinning.axisType == Acts::detail::AxisType::Equidistant) {
0156     // Variable-Equidistant
0157     Acts::GridAxisGenerators::VarEq<aType, bType> aGenerator{
0158         aBinning.edges,
0159         {bBinning.edges.front(), bBinning.edges.back()},
0160         bBinning.bins()};
0161     sfCandidates = isg(gctx, aGenerator, rGenerator);
0162   } else if (aBinning.axisType == Acts::detail::AxisType::Equidistant) {
0163     // Equidistant-Variable
0164     Acts::GridAxisGenerators::EqVar<aType, bType> aGenerator{
0165         {aBinning.edges.front(), aBinning.edges.back()},
0166         aBinning.bins(),
0167         bBinning.edges};
0168     sfCandidates = isg(gctx, aGenerator, rGenerator);
0169   } else {
0170     // Variable-Variable
0171     Acts::GridAxisGenerators::VarVar<aType, bType> aGenerator{aBinning.edges,
0172                                                               bBinning.edges};
0173     sfCandidates = isg(gctx, aGenerator, rGenerator);
0174   }
0175   // Return the candidates
0176   return sfCandidates;
0177 }
0178
0179 }  // namespace
0180
0181 Acts::Experimental::LayerStructureBuilder::LayerStructureBuilder(
0182     const Acts::Experimental::LayerStructureBuilder::Config& cfg,
0183     std::unique_ptr<const Acts::Logger> logger)
0184     : IInternalStructureBuilder(), m_cfg(cfg), m_logger(std::move(logger)) {
0185   if (m_cfg.surfacesProvider == nullptr) {
0186     throw std::invalid_argument(
0187         "LayerStructureBuilder: surfaces provider is nullptr.");
0188   }
0189 }
0190
0191 Acts::Experimental::InternalStructure
0192 Acts::Experimental::LayerStructureBuilder::construct(
0193     const Acts::GeometryContext& gctx) const {
0194   // Trivialities first: internal volumes
0195   std::vector<std::shared_ptr<DetectorVolume>> internalVolumes = {};
0196   DetectorVolumeUpdater internalVolumeUpdater = tryNoVolumes();
0197
0198   // Print the auxiliary information
0199   if (!m_cfg.auxiliary.empty()) {
0200     ACTS_DEBUG(m_cfg.auxiliary);
0201   }
0202
0203   // Retrieve the layer surfaces
0204   SurfaceCandidatesUpdater internalCandidatesUpdater =
0205       tryAllPortalsAndSurfaces();
0206   auto internalSurfaces = m_cfg.surfacesProvider->surfaces(gctx);
0207   ACTS_DEBUG("Building internal layer structure from "
0208              << internalSurfaces.size() << " provided surfaces.");
0209
0210   // Check whether support structure is scheduled to be built, and if so
0211   // collect those that should be assigned to all bins
0212   std::vector<std::size_t> assignToAll = {};
0213   if (!m_cfg.supports.empty()) {
0214     ACTS_DEBUG("Adding " << m_cfg.supports.size() << " support structures.")
0215     // The surface candidate updator
0216     for (const auto& support : m_cfg.supports) {
0217       // Check if the supportsurface has already been built
0218       if (support.surface != nullptr) {
0219         ACTS_VERBOSE("- Use provided support surface directly.");
0220         if (support.assignToAll) {
0221           assignToAll.push_back(internalSurfaces.size());
0222           ACTS_VERBOSE("  Support surface is assigned to all bins.");
0223         }
0224         internalSurfaces.push_back(support.surface);
0225         continue;
0226       }
0227
0228       // Throw an exception is misconfigured
0229       if (support.type == Surface::SurfaceType::Other) {
0230         throw std::invalid_argument(
0231             "LayerStructureBuilder: support surface type not specified.");
0232       }
0233       ACTS_VERBOSE("- Build support of type '"
0234                    << Acts::Surface::s_surfaceTypeNames[support.type] << "'.");
0235       if (support.splits > 1u) {
0236         ACTS_VERBOSE("  Support surface is modelled with " << support.splits
0237                                                            << " planes.");
0238       }
0239
0240       // The support extent
0241       Extent supportExtent;
0242       // Let us start with an eventually existing volume extent, but only pick
0243       // the binning value that are not constrained by the internal surfaces
0244       for (const auto& bv : s_binningValues) {
0245         if (support.volumeExtent.constrains(bv) &&
0246             std::find(support.internalConstraints.begin(),
0247                       support.internalConstraints.end(),
0248                       bv) == support.internalConstraints.end()) {
0249           ACTS_VERBOSE("  Support surface is constrained by volume extent in "
0250                        << binningValueNames()[bv]);
0251           supportExtent.set(bv, support.volumeExtent.min(bv),
0252                             support.volumeExtent.max(bv));
0253         }
0254       }
0255
0256       // Now add the internal constraints
0257       if (!support.internalConstraints.empty()) {
0258         // Estimate the extent from the surfaces
0259         for (const auto& s : internalSurfaces) {
0260           auto sPolyhedron = s->polyhedronRepresentation(gctx, m_cfg.nSegments);
0261           supportExtent.extend(sPolyhedron.extent(),
0262                                support.internalConstraints);
0263         }
0264       }
0265
0266       // Add cylindrical support
0267       if (support.type == Surface::SurfaceType::Cylinder) {
0268         detail::SupportSurfacesHelper::CylindricalSupport cSupport{
0269             support.offset, support.volumeClearance[binZ],
0270             support.volumeClearance[binPhi]};
0271         detail::SupportSurfacesHelper::addSupport(internalSurfaces, assignToAll,
0272                                                   supportExtent, cSupport,
0273                                                   support.splits);
0274       } else if (support.type == Surface::SurfaceType::Disc) {
0275         // Add disc support
0276         detail::SupportSurfacesHelper::DiscSupport dSupport{
0277             support.offset, support.volumeClearance[binR],
0278             support.volumeClearance[binPhi]};
0279         detail::SupportSurfacesHelper::addSupport(internalSurfaces, assignToAll,
0280                                                   supportExtent, dSupport,
0281                                                   support.splits);
0282       } else if (support.type == Surface::SurfaceType::Plane) {
0283         // Set the local coordinates - cyclic permutation
0284         std::array<BinningValue, 2> locals = {binX, binY};
0285         if (support.pPlacement == binX) {
0286           locals = {binY, binZ};
0287         } else if (support.pPlacement == binY) {
0288           locals = {binZ, binX};
0289         }
0290         // Add rectangular support
0291         detail::SupportSurfacesHelper::RectangularSupport rSupport{
0292             support.pPlacement, support.offset,
0293             support.volumeClearance[locals[0u]],
0294             support.volumeClearance[locals[1u]]};
0295         detail::SupportSurfacesHelper::addSupport(internalSurfaces, assignToAll,
0296                                                   supportExtent, rSupport);
0297       }
0298
0299       else {
0300         throw std::invalid_argument(
0301             "LayerStructureBuilder: support surface type not supported.");
0302       }
0303     }
0304   }
0305
0306   if (internalSurfaces.size() >= m_cfg.nMinimalSurfaces) {
0307     // Copy as we might patch it with the surface extent
0308     auto binnings = m_cfg.binnings;
0309
0310     if (binnings.empty()) {
0311       ACTS_DEBUG(
0312           "No surface binning provided, navigation will be 'tryAll' "
0313           "(potentially slow).");
0314     } else if (binnings.size() == 1u) {
0315       // Check if autorange for binning applies
0316       if (m_cfg.extent.has_value()) {
0317         ACTS_DEBUG("- adapting the proto binning range to the surface extent.");
0318         adaptBinningRange(binnings, m_cfg.extent.value());
0319       }
0320       ACTS_DEBUG("- 1-dimensional surface binning detected.");
0321       // Capture the binning
0322       auto binning = binnings[0u];
0323       if (binning.boundaryType == Acts::detail::AxisBoundaryType::Closed) {
0324         ACTS_VERBOSE("-- closed binning option.");
0325         internalCandidatesUpdater =
0326             createUpdater<Acts::detail::AxisBoundaryType::Closed>(
0327                 gctx, internalSurfaces, assignToAll, binning);
0328       } else {
0329         ACTS_VERBOSE("-- bound binning option.");
0330         internalCandidatesUpdater =
0331             createUpdater<Acts::detail::AxisBoundaryType::Bound>(
0332                 gctx, internalSurfaces, assignToAll, binning);
0333       }
0334     } else if (binnings.size() == 2u) {
0335       // Check if autorange for binning applies
0336       if (m_cfg.extent.has_value()) {
0337         ACTS_DEBUG(
0338             "- adapting the proto binning range(s) to the surface extent.");
0339         adaptBinningRange(binnings, m_cfg.extent.value());
0340       }
0341       // Sort the binning for conventions
0342       std::sort(binnings.begin(), binnings.end(),
0343                 [](const ProtoBinning& a, const ProtoBinning& b) {
0344                   return a.binValue < b.binValue;
0345                 });
0346
0347       ACTS_DEBUG("- 2-dimensional surface binning detected.");
0348       // Capture the binnings
0349       const auto& binning0 = binnings[0u];
0350       const auto& binning1 = binnings[1u];
0351
0352       if (binning0.boundaryType == Acts::detail::AxisBoundaryType::Closed) {
0353         ACTS_VERBOSE("-- closed/bound binning option.");
0354         internalCandidatesUpdater =
0355             createUpdater<Acts::detail::AxisBoundaryType::Closed,
0356                           Acts::detail::AxisBoundaryType::Bound>(
0357                 gctx, internalSurfaces, assignToAll, binning0, binning1);
0358       } else if (binning1.boundaryType ==
0359                  Acts::detail::AxisBoundaryType::Closed) {
0360         ACTS_VERBOSE("-- bound/closed binning option.");
0361         internalCandidatesUpdater =
0362             createUpdater<Acts::detail::AxisBoundaryType::Bound,
0363                           Acts::detail::AxisBoundaryType::Closed>(
0364                 gctx, internalSurfaces, assignToAll, binning0, binning1);
0365       } else {
0366         ACTS_VERBOSE("-- bound/bound binning option.");
0367         internalCandidatesUpdater =
0368             createUpdater<Acts::detail::AxisBoundaryType::Bound,
0369                           Acts::detail::AxisBoundaryType::Bound>(
0370                 gctx, internalSurfaces, assignToAll, binning0, binning1);
0371       }
0372     }
0373   } else {
0374     ACTS_DEBUG("Only " << internalSurfaces.size() << " surfaces provided, "
0375                        << "navigation will be 'tryAll'");
0376     ACTS_DEBUG("Per configuration " << m_cfg.nMinimalSurfaces
0377                                     << " surfaces are "
0378                                     << "required to use the surface binning.");
0379   }
0380
0381   // Check if everything went ok
0382   if (!internalCandidatesUpdater.connected()) {
0383     throw std::runtime_error(
0384         "LayerStructureBuilder: could not connect surface candidate updator.");
0385   }
0386
0387   // Return the internal structure
0388   return InternalStructure{internalSurfaces, internalVolumes,
0389                            std::move(internalCandidatesUpdater),
0390                            std::move(internalVolumeUpdater)};
0391 }