Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-07-08 07:50:54

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 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Geometry/GeometryContext.hpp"
0013 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0014 #include "Acts/Surfaces/CylinderBounds.hpp"
0015 #include "Acts/Surfaces/RegularSurface.hpp"
0016 #include "Acts/Surfaces/Surface.hpp"
0017 #include "Acts/Utilities/AnyGridView.hpp"
0018 #include "Acts/Utilities/AxisDefinitions.hpp"
0019 #include "Acts/Utilities/Grid.hpp"
0020 #include "Acts/Utilities/IAxis.hpp"
0021 #include "Acts/Utilities/Intersection.hpp"
0022 
0023 #include <iostream>
0024 #include <limits>
0025 #include <vector>
0026 
0027 namespace Acts {
0028 
0029 using SurfaceVector = std::vector<const Surface*>;
0030 
0031 /// @brief Provides Surface binning in 2 dimensions
0032 ///
0033 /// Uses @c Grid under the hood to implement the storage and lookup
0034 /// Contains a lookup struct which talks to the @c Grid
0035 /// and performs utility actions. This struct needs to be initialised
0036 /// externally and passed to @c SurfaceArray on construction.
0037 class SurfaceArray {
0038  public:
0039   /// @brief Base interface for all surface lookups.
0040   struct ISurfaceGridLookup {
0041     /// @brief Fill provided surfaces into the contained @c Grid.
0042     /// @param gctx The current geometry context object, e.g. alignment
0043     /// @param surfaces Input surface pointers
0044     virtual void fill(const GeometryContext& gctx,
0045                       const SurfaceVector& surfaces) = 0;
0046 
0047     /// @brief Performs lookup at @c pos and returns bin content as const
0048     /// reference
0049     /// @param position Lookup position
0050     /// @param direction Lookup direction
0051     /// @return @c SurfaceVector at given bin
0052     virtual const SurfaceVector& lookup(const Vector3& position,
0053                                         const Vector3& direction) const = 0;
0054 
0055     /// @brief Performs lookup at global bin and returns bin content as
0056     /// reference
0057     /// @param bin Global lookup bin
0058     /// @return @c SurfaceVector at given bin
0059     virtual SurfaceVector& lookup(std::size_t bin) = 0;
0060 
0061     /// @brief Performs lookup at global bin and returns bin content as const
0062     /// reference
0063     /// @param bin Global lookup bin
0064     /// @return @c SurfaceVector at given bin
0065     virtual const SurfaceVector& lookup(std::size_t bin) const = 0;
0066 
0067     /// @brief Performs a lookup at @c pos, but returns neighbors as well
0068     ///
0069     /// @param position Lookup position
0070     /// @param direction Lookup direction
0071     /// @return @c SurfaceVector at given bin. Copy of all bins selected
0072     virtual const SurfaceVector& neighbors(const Vector3& position,
0073                                            const Vector3& direction) const = 0;
0074 
0075     /// @brief Returns the total size of the grid (including under/overflow
0076     /// bins)
0077     /// @return Size of the grid data structure
0078     virtual std::size_t size() const = 0;
0079 
0080     /// @brief Gets the center position of bin @c bin in global coordinates
0081     /// @param bin the global bin index
0082     /// @return The bin center
0083     virtual Vector3 getBinCenter(std::size_t bin) const = 0;
0084 
0085     /// @brief Returns copies of the axes used in the grid as @c AnyAxis
0086     /// @return The axes
0087     /// @note This returns copies. Use for introspection and querying.
0088     virtual std::vector<const IAxis*> getAxes() const = 0;
0089 
0090     /// @brief Get a view of the grid for inspection
0091     /// @return Optional grid view containing surface vectors
0092     virtual std::optional<AnyGridConstView<SurfaceVector>> getGridView()
0093         const = 0;
0094 
0095     /// @brief Get the representative surface used for this lookup
0096     /// @return Surface pointer
0097     virtual const Surface* surfaceRepresentation() const = 0;
0098 
0099     /// @brief Checks if global bin is valid
0100     /// @param bin the global bin index
0101     /// @return bool if the bin is valid
0102     /// @note Valid means that the index points to a bin which is not a under
0103     ///       or overflow bin or out of range in any axis.
0104     virtual bool isValidBin(std::size_t bin) const = 0;
0105 
0106     /// @brief The binning values described by this surface grid lookup
0107     /// They are in order of the axes (optional) and empty for eingle lookups
0108     /// @return Vector of axis directions for binning
0109     virtual std::vector<AxisDirection> binningValues() const { return {}; };
0110 
0111     /// Pure virtual destructor
0112     virtual ~ISurfaceGridLookup() = 0;
0113   };
0114 
0115   /// @brief Lookup helper which encapsulates a @c Grid
0116   /// @tparam Axes The axes used for the grid
0117   template <class Axis1, class Axis2>
0118   struct SurfaceGridLookup : ISurfaceGridLookup {
0119     /// Grid type storing surface vectors with two axes
0120     using Grid_t = Grid<SurfaceVector, Axis1, Axis2>;
0121 
0122     /// Construct a surface grid lookup
0123     /// @param representative The surface which is used as representative
0124     /// @param tolerance The tolerance used for intersection checks
0125     /// @param axes The axes used for the grid
0126     /// @param bValues Optional vector of axis directions for binning
0127     SurfaceGridLookup(std::shared_ptr<RegularSurface> representative,
0128                       double tolerance, std::tuple<Axis1, Axis2> axes,
0129                       std::vector<AxisDirection> bValues = {})
0130         : m_representative(std::move(representative)),
0131           m_tolerance(tolerance),
0132           m_grid(std::move(axes)),
0133           m_binValues(std::move(bValues)) {
0134       m_neighborMap.resize(m_grid.size());
0135     }
0136 
0137     /// @brief Fill provided surfaces into the contained @c Grid.
0138     ///
0139     /// This is done by iterating, accessing the referencePosition, lookup
0140     /// and append.
0141     /// Also populates the neighbor map by combining the filled bins of
0142     /// all bins around a given one.
0143     ///
0144     /// @param gctx The current geometry context object, e.g. alignment
0145     /// @param surfaces Input surface pointers
0146     void fill(const GeometryContext& gctx,
0147               const SurfaceVector& surfaces) override {
0148       for (const Surface* surface : surfaces) {
0149         const std::size_t globalBin = fillSurfaceToBinMapping(gctx, *surface);
0150         if (globalBin == 0) {
0151           continue;
0152         }
0153 
0154         fillBinToSurfaceMapping(gctx, *surface, globalBin);
0155       }
0156 
0157       populateNeighborCache();
0158     }
0159 
0160     const SurfaceVector& lookup(const Vector3& position,
0161                                 const Vector3& direction) const override {
0162       return m_grid.at(findGlobalBin(position, direction,
0163                                      std::numeric_limits<double>::infinity()));
0164     }
0165 
0166     /// @brief Performs lookup at global bin and returns bin content as
0167     /// reference
0168     /// @param bin Global lookup bin
0169     /// @return @c SurfaceVector at given bin
0170     SurfaceVector& lookup(std::size_t bin) override { return m_grid.at(bin); }
0171 
0172     /// @brief Performs lookup at global bin and returns bin content as const
0173     /// reference
0174     /// @param bin Global lookup bin
0175     /// @return @c SurfaceVector at given bin
0176     const SurfaceVector& lookup(std::size_t bin) const override {
0177       return m_grid.at(bin);
0178     }
0179 
0180     /// @brief Performs a lookup at @c pos, but returns neighbors as well
0181     ///
0182     /// @param position Lookup position
0183     /// @param direction Lookup direction
0184     /// @return @c SurfaceVector at given bin. Copy of all bins selected
0185     const SurfaceVector& neighbors(const Vector3& position,
0186                                    const Vector3& direction) const override {
0187       return m_neighborMap.at(findGlobalBin(
0188           position, direction, std::numeric_limits<double>::infinity()));
0189     }
0190 
0191     /// @brief Returns the total size of the grid (including under/overflow
0192     /// bins)
0193     /// @return Size of the grid data structure
0194     std::size_t size() const override { return m_grid.size(); }
0195 
0196     /// @brief The binning values described by this surface grid lookup
0197     /// They are in order of the axes
0198     /// @return Vector of axis directions for binning
0199     std::vector<AxisDirection> binningValues() const override {
0200       return m_binValues;
0201     }
0202 
0203     /// @brief Gets the center position of bin @c bin in global coordinates
0204     /// @param bin the global bin index
0205     /// @return The bin center
0206     Vector3 getBinCenter(std::size_t bin) const override {
0207       auto gctx = GeometryContext::dangerouslyDefaultConstruct();
0208       return getBinCenterImpl(gctx, bin);
0209     }
0210 
0211     /// @brief Returns copies of the axes used in the grid as @c AnyAxis
0212     /// @return The axes
0213     /// @note This returns copies. Use for introspection and querying.
0214     std::vector<const IAxis*> getAxes() const override {
0215       auto arr = m_grid.axes();
0216       return std::vector<const IAxis*>(arr.begin(), arr.end());
0217     }
0218 
0219     std::optional<AnyGridConstView<SurfaceVector>> getGridView()
0220         const override {
0221       return AnyGridConstView<SurfaceVector>{m_grid};
0222     }
0223 
0224     const Surface* surfaceRepresentation() const override {
0225       return m_representative.get();
0226     }
0227 
0228     /// @brief Checks if global bin is valid
0229     /// @param bin the global bin index
0230     /// @return bool if the bin is valid
0231     /// @note Valid means that the index points to a bin which is not a under
0232     ///       or overflow bin or out of range in any axis.
0233     bool isValidBin(std::size_t bin) const override {
0234       std::array<std::size_t, 2> indices = m_grid.localBinsFromGlobalBin(bin);
0235       std::array<std::size_t, 2> nBins = m_grid.numLocalBins();
0236       for (std::size_t i = 0; i < indices.size(); ++i) {
0237         std::size_t idx = indices.at(i);
0238         if (idx <= 0 || idx >= nBins.at(i) + 1) {
0239           return false;
0240         }
0241       }
0242       return true;
0243     }
0244 
0245    private:
0246     /// map surface center to grid
0247     std::size_t fillSurfaceToBinMapping(const GeometryContext& gctx,
0248                                         const Surface& surface) {
0249       const Vector3 pos = surface.referencePosition(gctx, AxisDirection::AxisR);
0250       const Vector3 normal = m_representative->normal(gctx, pos);
0251       const std::size_t globalBin = findGlobalBin(pos, normal, m_tolerance);
0252       if (globalBin != 0) {
0253         m_grid.at(globalBin).push_back(&surface);
0254       }
0255       return globalBin;
0256     };
0257 
0258     /// flood fill neighboring bins given a starting bin
0259     void fillBinToSurfaceMapping(const GeometryContext& gctx,
0260                                  const Surface& surface, std::size_t startBin) {
0261       const std::array<std::size_t, 2> startIndices =
0262           m_grid.localBinsFromGlobalBin(startBin);
0263       const auto startNeighborIndices =
0264           m_grid.neighborHoodIndices(startIndices, 1u);
0265 
0266       std::set<std::size_t> visited({startBin});
0267       std::vector<std::size_t> queue(startNeighborIndices.begin(),
0268                                      startNeighborIndices.end());
0269 
0270       while (!queue.empty()) {
0271         const std::size_t current = queue.back();
0272         queue.pop_back();
0273         if (visited.contains(current)) {
0274           continue;
0275         }
0276 
0277         const std::array<std::size_t, 2> currentIndices =
0278             m_grid.localBinsFromGlobalBin(current);
0279         visited.insert(current);
0280 
0281         const std::array<double, 2> gridLocal =
0282             m_grid.binCenter(currentIndices);
0283         const Vector2 surfaceLocal = gridToSurfaceLocal(gridLocal);
0284         const Vector3 normal = m_representative->normal(gctx, surfaceLocal);
0285         const Vector3 global =
0286             m_representative->localToGlobal(gctx, surfaceLocal, normal);
0287 
0288         const Intersection3D intersection =
0289             surface.intersect(gctx, global, normal, BoundaryTolerance::None())
0290                 .closest();
0291         if (!intersection.isValid() ||
0292             std::abs(intersection.pathLength()) > m_tolerance) {
0293           continue;
0294         }
0295         m_grid.at(current).push_back(&surface);
0296 
0297         const auto neighborIndices =
0298             m_grid.neighborHoodIndices(currentIndices, 1u);
0299         queue.insert(queue.end(), neighborIndices.begin(),
0300                      neighborIndices.end());
0301       }
0302     };
0303 
0304     void populateNeighborCache() {
0305       // calculate neighbors for every bin and store in map
0306       for (std::size_t i = 0; i < m_grid.size(); i++) {
0307         if (!isValidBin(i)) {
0308           continue;
0309         }
0310         const std::array<std::size_t, 2> indices =
0311             m_grid.localBinsFromGlobalBin(i);
0312         std::vector<const Surface*>& neighbors = m_neighborMap.at(i);
0313         neighbors.clear();
0314 
0315         for (std::size_t idx : m_grid.neighborHoodIndices(indices, 1u)) {
0316           const std::vector<const Surface*>& binContent = m_grid.at(idx);
0317           std::copy(binContent.begin(), binContent.end(),
0318                     std::back_inserter(neighbors));
0319         }
0320 
0321         std::ranges::sort(neighbors);
0322         auto last = std::ranges::unique(neighbors);
0323         neighbors.erase(last.begin(), last.end());
0324         neighbors.shrink_to_fit();
0325       }
0326     }
0327 
0328     Vector3 getBinCenterImpl(const GeometryContext& gctx,
0329                              std::size_t bin) const {
0330       const std::array<double, 2> gridLocal =
0331           m_grid.binCenter(m_grid.localBinsFromGlobalBin(bin));
0332       const Vector2 surfaceLocal = gridToSurfaceLocal(gridLocal);
0333       return m_representative->localToGlobal(gctx, surfaceLocal);
0334     }
0335 
0336     const CylinderBounds* getCylinderBounds() const {
0337       return dynamic_cast<const CylinderBounds*>(&m_representative->bounds());
0338     }
0339 
0340     Vector2 gridToSurfaceLocal(std::array<double, 2> gridLocal) const {
0341       Vector2 surfaceLocal = Eigen::Map<Vector2>(gridLocal.data());
0342       if (const CylinderBounds* bounds = getCylinderBounds();
0343           bounds != nullptr) {
0344         surfaceLocal[0] *= bounds->get(CylinderBounds::eR);
0345       }
0346       return surfaceLocal;
0347     }
0348     std::array<double, 2> surfaceToGridLocal(Vector2 local) const {
0349       std::array<double, 2> gridLocal = {local[0], local[1]};
0350       if (const CylinderBounds* bounds = getCylinderBounds();
0351           bounds != nullptr) {
0352         gridLocal[0] /= bounds->get(CylinderBounds::eR);
0353       }
0354       return gridLocal;
0355     }
0356 
0357     std::size_t findGlobalBin(const Vector3& position, const Vector3& direction,
0358                               double tolerance) const {
0359       auto gctx = GeometryContext::dangerouslyDefaultConstruct();
0360 
0361       const Intersection3D intersection =
0362           m_representative
0363               ->intersect(gctx, position, direction,
0364                           BoundaryTolerance::Infinite())
0365               .closest();
0366       if (!intersection.isValid() ||
0367           std::abs(intersection.pathLength()) > tolerance) {
0368         return 0;  // overflow bin
0369       }
0370       const Vector2 surfaceLocal =
0371           m_representative
0372               ->globalToLocal(gctx, intersection.position(), direction)
0373               .value();
0374       const std::array<double, 2> gridLocal = surfaceToGridLocal(surfaceLocal);
0375       return m_grid.globalBinFromPosition(gridLocal);
0376     }
0377 
0378     std::shared_ptr<RegularSurface> m_representative;
0379     double m_tolerance{};
0380     Grid_t m_grid;
0381     std::vector<AxisDirection> m_binValues;
0382     std::vector<SurfaceVector> m_neighborMap;
0383   };
0384 
0385   /// @brief Lookup implementation which wraps one element and always returns
0386   ///        this element when lookup is called
0387   struct SingleElementLookup : ISurfaceGridLookup {
0388     /// @brief Default constructor.
0389     /// @param element the one and only element.
0390     explicit SingleElementLookup(SurfaceVector::value_type element)
0391         : m_element({element}) {}
0392 
0393     /// @brief Default constructor.
0394     /// @param elements the surfaces that are provided through a single lookup
0395     explicit SingleElementLookup(const SurfaceVector& elements)
0396         : m_element(elements) {}
0397 
0398     /// @brief Lookup, always returns @c element
0399     /// @return reference to vector containing only @c element
0400     const SurfaceVector& lookup(const Vector3& /*position*/,
0401                                 const Vector3& /*direction*/) const override {
0402       return m_element;
0403     }
0404 
0405     /// @brief Lookup, always returns @c element
0406     /// @return reference to vector containing only @c element
0407     SurfaceVector& lookup(std::size_t /*bin*/) override { return m_element; }
0408 
0409     /// @brief Lookup, always returns @c element
0410     /// @return reference to vector containing only @c element
0411     const SurfaceVector& lookup(std::size_t /*bin*/) const override {
0412       return m_element;
0413     }
0414 
0415     /// @brief Lookup, always returns @c element
0416     /// @return reference to vector containing only @c element
0417     const SurfaceVector& neighbors(
0418         const Vector3& /*position*/,
0419         const Vector3& /*direction*/) const override {
0420       return m_element;
0421     }
0422 
0423     /// @brief returns 1
0424     /// @return 1
0425     std::size_t size() const override { return 1; }
0426 
0427     /// @brief Gets the bin center, but always returns (0, 0, 0)
0428     /// @return (0, 0, 0)
0429     Vector3 getBinCenter(std::size_t /*bin*/) const override {
0430       return Vector3(0, 0, 0);
0431     }
0432 
0433     /// @brief Returns an empty vector of @c AnyAxis
0434     /// @return empty vector
0435     std::vector<const IAxis*> getAxes() const override { return {}; }
0436 
0437     std::optional<AnyGridConstView<SurfaceVector>> getGridView()
0438         const override {
0439       return std::nullopt;
0440     }
0441 
0442     const Surface* surfaceRepresentation() const override { return nullptr; }
0443 
0444     /// @brief Comply with concept and provide fill method
0445     /// @note Does nothing
0446     void fill(const GeometryContext& /*gctx*/,
0447               const SurfaceVector& /*surfaces*/) override {}
0448 
0449     /// @brief Returns if the bin is valid (it is)
0450     /// @return always true
0451     bool isValidBin(std::size_t /*bin*/) const override { return true; }
0452 
0453    private:
0454     SurfaceVector m_element;
0455   };
0456 
0457   /// @brief Default constructor which takes a @c SurfaceLookup and a vector of
0458   /// surfaces
0459   /// @param gridLookup The grid storage. @c SurfaceArray does not fill it on
0460   /// its own
0461   /// @param surfaces The input vector of surfaces. This is only for
0462   /// bookkeeping, so we can ask
0463   /// @param transform Optional additional transform for this SurfaceArray
0464   explicit SurfaceArray(std::unique_ptr<ISurfaceGridLookup> gridLookup,
0465                         std::vector<std::shared_ptr<const Surface>> surfaces,
0466                         const Transform3& transform = Transform3::Identity());
0467 
0468   /// @brief Constructor with a single surface
0469   /// @param srf The one and only surface
0470   explicit SurfaceArray(std::shared_ptr<const Surface> srf);
0471 
0472   /// @brief Get all surfaces in bin given by position @p pos.
0473   /// @param position the lookup position
0474   /// @param direction the lookup direction
0475   /// @return const reference to @c SurfaceVector contained in bin at that
0476   /// position
0477   const SurfaceVector& at(const Vector3& position,
0478                           const Vector3& direction) const {
0479     return p_gridLookup->lookup(position, direction);
0480   }
0481 
0482   /// @brief Get all surfaces in bin given by global bin index @p bin.
0483   /// @param bin the global bin index
0484   /// @return reference to @c SurfaceVector contained in bin
0485   SurfaceVector& at(std::size_t bin) { return p_gridLookup->lookup(bin); }
0486 
0487   /// @brief Get all surfaces in bin given by global bin index.
0488   /// @param bin the global bin index
0489   /// @return const reference to @c SurfaceVector contained in bin
0490   const SurfaceVector& at(std::size_t bin) const {
0491     return p_gridLookup->lookup(bin);
0492   }
0493 
0494   /// @brief Get all surfaces in bin at @p pos and its neighbors
0495   /// @param position The position to lookup
0496   /// @param direction The direction to lookup
0497   /// @return Merged @c SurfaceVector of neighbors and nominal
0498   /// @note The @c SurfaceVector will be combined. For technical reasons, the
0499   ///       different bin content vectors have to be copied, so the resulting
0500   ///       vector contains copies.
0501   const SurfaceVector& neighbors(const Vector3& position,
0502                                  const Vector3& direction) const {
0503     return p_gridLookup->neighbors(position, direction);
0504   }
0505 
0506   /// @brief Get the size of the underlying grid structure including
0507   /// under/overflow bins
0508   /// @return the size
0509   std::size_t size() const { return p_gridLookup->size(); }
0510 
0511   /// @brief Get the center of the bin identified by global bin index @p bin
0512   /// @param bin the global bin index
0513   /// @return Center position of the bin in global coordinates
0514   Vector3 getBinCenter(std::size_t bin) const {
0515     return p_gridLookup->getBinCenter(bin);
0516   }
0517 
0518   /// @brief Get all surfaces attached to this @c SurfaceArray
0519   /// @return Reference to @c SurfaceVector containing all surfaces
0520   /// @note This does not reflect the actual state of the grid. It only
0521   ///       returns what was given in the constructor, without any checks
0522   ///       if that is actually what's in the grid.
0523   const SurfaceVector& surfaces() const { return m_surfacesRawPointers; }
0524 
0525   /// @brief Get vector of axes spanning the grid as @c AnyAxis
0526   /// @return vector of @c AnyAxis
0527   /// @note The axes in the vector are copies. Only use for introspection and
0528   ///       querying.
0529   std::vector<const IAxis*> getAxes() const { return p_gridLookup->getAxes(); }
0530 
0531   /// @brief Checks if global bin is valid
0532   /// @param bin the global bin index
0533   /// @return bool if the bin is valid
0534   /// @note Valid means that the index points to a bin which is not a under
0535   ///       or overflow bin or out of range in any axis.
0536   bool isValidBin(std::size_t bin) const {
0537     return p_gridLookup->isValidBin(bin);
0538   }
0539 
0540   /// Get the transform of this surface array.
0541   /// @return Reference to the transformation matrix
0542   const Transform3& transform() const { return m_transform; }
0543 
0544   /// @brief The binning values described by this surface grid lookup
0545   /// They are in order of the axes
0546   /// @return Vector of axis directions for binning
0547   std::vector<AxisDirection> binningValues() const {
0548     return p_gridLookup->binningValues();
0549   };
0550 
0551   /// @brief String representation of this @c SurfaceArray
0552   /// @param gctx The current geometry context object, e.g. alignment
0553   /// @param sl Output stream to write to
0554   /// @return the output stream given as @p sl
0555   std::ostream& toStream(const GeometryContext& gctx, std::ostream& sl) const;
0556 
0557   /// Return the lookup object
0558   /// @return Reference to the surface grid lookup interface
0559   const ISurfaceGridLookup& gridLookup() const { return *p_gridLookup; }
0560 
0561  private:
0562   std::unique_ptr<ISurfaceGridLookup> p_gridLookup;
0563   // this vector makes sure we have shared ownership over the surfaces
0564   std::vector<std::shared_ptr<const Surface>> m_surfaces;
0565   // this vector is returned, so that (expensive) copying of the shared_ptr
0566   // vector does not happen by default
0567   SurfaceVector m_surfacesRawPointers;
0568   // this is only used to keep info on transform applied
0569   // by l2g and g2l
0570   Transform3 m_transform;
0571 };
0572 
0573 }  // namespace Acts