Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:11:04

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