Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-04-21 07:46:50

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/Definitions/Units.hpp"
0013 #include "Acts/Geometry/GeometryContext.hpp"
0014 #include "Acts/Geometry/ProtoLayer.hpp"
0015 #include "Acts/Surfaces/RegularSurface.hpp"
0016 #include "Acts/Surfaces/Surface.hpp"
0017 #include "Acts/Surfaces/SurfaceArray.hpp"
0018 #include "Acts/Utilities/Axis.hpp"
0019 #include "Acts/Utilities/AxisDefinitions.hpp"
0020 #include "Acts/Utilities/BinningType.hpp"
0021 #include "Acts/Utilities/Logger.hpp"
0022 
0023 #include <algorithm>
0024 #include <cmath>
0025 #include <cstddef>
0026 #include <functional>
0027 #include <iterator>
0028 #include <memory>
0029 #include <numbers>
0030 #include <optional>
0031 #include <utility>
0032 #include <vector>
0033 
0034 namespace ActsTests {
0035 struct SurfaceArrayCreatorFixture;
0036 }
0037 
0038 namespace Acts {
0039 
0040 /// Function type for comparing two surfaces in a given geometry context and
0041 /// axis direction.
0042 /// @param gctx The geometry context for the comparison
0043 /// @param dir The axis direction to consider
0044 /// @param s1 First surface to compare
0045 /// @param s2 Second surface to compare
0046 /// @return True if the surfaces are considered equivalent for binning purposes
0047 using SurfaceMatcher =
0048     std::function<bool(const GeometryContext& gctx, AxisDirection,
0049                        const Surface*, const Surface*)>;
0050 
0051 /// Vector of pointers to constant Surface objects
0052 using SurfaceVector = std::vector<const Surface*>;
0053 /// Matrix (2D vector) of pointers to constant Surface objects
0054 using SurfaceMatrix = std::vector<SurfaceVector>;
0055 
0056 /// @typedef V3Vector
0057 /// Vector of 3D vectors, used for storing collections of 3D points.
0058 using V3Vector = std::vector<Vector3>;
0059 
0060 /// @typedef V3Matrix
0061 /// Matrix (2D vector) of 3D vectors, used for storing grid-like collections of
0062 /// 3D points.
0063 using V3Matrix = std::vector<V3Vector>;
0064 
0065 /// @brief Scalar type used for axis values in surface array binning
0066 using AxisScalar = Vector3::Scalar;
0067 
0068 /// @class SurfaceArrayCreator
0069 ///
0070 /// It is designed create sub surface arrays to be ordered on Surfaces
0071 ///
0072 /// @todo write more documentation on how this is done
0073 class SurfaceArrayCreator {
0074  public:
0075   friend struct ActsTests::SurfaceArrayCreatorFixture;
0076   friend class SurfaceArray;
0077 
0078   /// Prototype axis definition for surface binning.
0079   struct ProtoAxis {
0080     /// Binning type (equidistant or variable)
0081     BinningType bType = BinningType::equidistant;
0082     /// Axis direction for binning
0083     AxisDirection axisDir = AxisDirection::AxisX;
0084     /// Number of bins
0085     std::size_t nBins = 0;
0086     /// Minimum value of the axis
0087     AxisScalar min = 0;
0088     /// Maximum value of the axis
0089     AxisScalar max = 0;
0090     /// Bin edges for variable binning
0091     std::vector<AxisScalar> binEdges;
0092 
0093     /// Get the bin index for a given value
0094     /// @param x The value to find the bin for
0095     /// @return The bin index
0096     std::size_t getBin(AxisScalar x) const {
0097       if (binEdges.empty()) {
0098         // equidistant
0099         AxisScalar w = (max - min) / nBins;
0100         return static_cast<std::size_t>(std::floor((x - min) / w));
0101       } else {
0102         // variable
0103         const auto it = std::ranges::upper_bound(binEdges, x);
0104         return static_cast<std::size_t>(
0105                    std::ranges::distance(binEdges.begin(), it)) -
0106                1;
0107       }
0108     }
0109   };
0110 
0111   /// Configuration options for the surface array creator.
0112   struct Config {
0113     /// Type-erased function which determines whether two surfaces are
0114     /// supposed to be considered equivalent in terms of the binning
0115     SurfaceMatcher surfaceMatcher = SurfaceArrayCreator::isSurfaceEquivalent;
0116 
0117     /// Optimize the binning in phi for disc layers. Reduces the number
0118     /// of bins to the lowest number of non-equivalent phi surfaces
0119     /// of all r-bins. If false, this step is skipped.
0120     bool doPhiBinningOptimization = true;
0121   };
0122 
0123   /// Constructor with default config
0124   ///
0125   /// @param logger logging instance
0126   explicit SurfaceArrayCreator(std::unique_ptr<const Logger> logger =
0127                                    getDefaultLogger("SurfaceArrayCreator",
0128                                                     Logging::INFO))
0129       : m_cfg(Config()), m_logger(std::move(logger)) {}
0130   /// Constructor with explicit config
0131   ///
0132   /// @param cfg Explicit config struct
0133   /// @param logger logging instance
0134   explicit SurfaceArrayCreator(const Config& cfg,
0135                                std::unique_ptr<const Logger> logger =
0136                                    getDefaultLogger("SurfaceArrayCreator",
0137                                                     Logging::INFO))
0138       : m_cfg(cfg), m_logger(std::move(logger)) {}
0139 
0140   /// Destructor
0141   virtual ~SurfaceArrayCreator() = default;
0142 
0143   /// SurfaceArrayCreator interface method
0144   ///
0145   /// - create an array in a cylinder, binned in phi, z when extrema and
0146   /// bin numbers are known
0147   /// @warning This function requires the cylinder aligned with the z-axis
0148   /// @param surfaces is the vector of pointers to sensitive surfaces
0149   /// to be ordered on the cylinder
0150   /// @pre the pointers to the sensitive surfaces in the surfaces vectors all
0151   /// need to be valid, since no check is performed
0152   /// @param [in] gctx The geometry context for this building call
0153   /// @param protoLayerOpt The proto layer containing the layer size
0154   /// @param binsPhi is the number of bins in phi for the surfaces
0155   /// @param binsZ is the number of bin in Z for the surfaces
0156   /// @param transform is the (optional) additional transform applied
0157   /// @param maxNeighborDistance Maximum next neighbor distance to be included in neighbor lookups
0158   ///
0159   /// @return a unique pointer to a new SurfaceArray
0160   std::unique_ptr<SurfaceArray> surfaceArrayOnCylinder(
0161       const GeometryContext& gctx,
0162       std::vector<std::shared_ptr<const Surface>> surfaces, std::size_t binsPhi,
0163       std::size_t binsZ, std::optional<ProtoLayer> protoLayerOpt = std::nullopt,
0164       const Transform3& transform = Transform3::Identity(),
0165       std::uint8_t maxNeighborDistance = 1) const;
0166 
0167   /// SurfaceArrayCreator interface method
0168   ///
0169   /// - create an array in a cylinder, binned in phi, z when extrema and bin
0170   /// numbers are unknown - this method goes through the surfaces and finds
0171   /// out the needed information
0172   /// @warning This function requires the cylinder aligned with the z-axis
0173   /// @param surfaces is the vector of pointers to sensitive surfaces
0174   /// to be ordered on the cylinder
0175   /// @pre the pointers to the sensitive surfaces in the surfaces vectors all
0176   /// need to be valid, since no check is performed
0177   /// @param [in] gctx The geometry context for this building call
0178   /// @param protoLayerOpt The proto layer containing the layer size
0179   /// @param bTypePhi the binning type in phi direction (equidistant/arbitrary)
0180   /// @param bTypeZ the binning type in z direction (equidistant/arbitrary)
0181   /// @param transform is the (optional) additional transform applied
0182   /// @param maxNeighborDistance Maximum next neighbor distance to be included in neighbor lookups
0183   ///
0184   /// @return a unique pointer a new SurfaceArray
0185   std::unique_ptr<SurfaceArray> surfaceArrayOnCylinder(
0186       const GeometryContext& gctx,
0187       std::vector<std::shared_ptr<const Surface>> surfaces,
0188       BinningType bTypePhi = equidistant, BinningType bTypeZ = equidistant,
0189       std::optional<ProtoLayer> protoLayerOpt = std::nullopt,
0190       const Transform3& transform = Transform3::Identity(),
0191       std::uint8_t maxNeighborDistance = 1) const;
0192 
0193   /// SurfaceArrayCreator interface method
0194   /// - create an array on a disc, binned in r, phi when extrema and
0195   /// bin numbers are known
0196   ///
0197   /// @param surfaces is the vector of pointers to sensitive surfaces
0198   /// to be ordered on the disc
0199   /// @pre the pointers to the sensitive surfaces in the surfaces vectors all
0200   /// need to be valid, since no check is performed
0201   /// @warning This function requires the disc aligned with the z-axis
0202   /// @param [in] gctx The geometry context for this building call
0203   /// @param protoLayerOpt The proto layer containing the layer size
0204   /// @param binsPhi is the number of bins in phi for the surfaces
0205   /// @param binsR is the number of bin in R for the surfaces
0206   /// @param transform is the (optional) additional transform applied
0207   /// @param maxNeighborDistance Maximum next neighbor distance to be included in neighbor lookups
0208   ///
0209   /// @return a unique pointer a new SurfaceArray
0210   std::unique_ptr<SurfaceArray> surfaceArrayOnDisc(
0211       const GeometryContext& gctx,
0212       std::vector<std::shared_ptr<const Surface>> surfaces, std::size_t binsR,
0213       std::size_t binsPhi,
0214       std::optional<ProtoLayer> protoLayerOpt = std::nullopt,
0215       const Transform3& transform = Transform3::Identity(),
0216       std::uint8_t maxNeighborDistance = 1) const;
0217 
0218   /// SurfaceArrayCreator interface method
0219   ///
0220   /// - create an array in a cylinder, binned in phi, r when extrema and bin
0221   /// numbers are unknown - this method goes through the surfaces and finds
0222   /// out the needed information
0223   /// @param surfaces is the vector of pointers to sensitive surfaces
0224   /// to be ordered on the disc
0225   /// @pre the pointers to the sensitive surfaces in the surfaces vectors all
0226   /// need to be valid, since no check is performed
0227   /// @warning This function requires the disc aligned with the z-axis
0228   /// @param [in] gctx The geometry context for this building call
0229   /// @param protoLayerOpt The proto layer containing the layer size
0230   /// @param bTypeR the binning type in r direction (equidistant/arbitrary)
0231   /// @param bTypePhi the binning type in phi direction (equidistant/arbitrary)
0232   /// @param transform is the (optional) additional transform applied
0233   /// @param maxNeighborDistance Maximum next neighbor distance to be included in neighbor lookups
0234   ///
0235   /// @return a unique pointer a new SurfaceArray
0236   ///
0237   /// @note If there is more than on R-Ring, number of phi bins
0238   ///       will be set to lowest number of surfaces of any R-ring.
0239   ///       This ignores bTypePhi and produces equidistant binning in phi
0240   std::unique_ptr<SurfaceArray> surfaceArrayOnDisc(
0241       const GeometryContext& gctx,
0242       std::vector<std::shared_ptr<const Surface>> surfaces, BinningType bTypeR,
0243       BinningType bTypePhi,
0244       std::optional<ProtoLayer> protoLayerOpt = std::nullopt,
0245       const Transform3& transform = Transform3::Identity(),
0246       std::uint8_t maxNeighborDistance = 1) const;
0247 
0248   /// SurfaceArrayCreator interface method
0249   /// - create an array on a plane
0250   ///
0251   /// @param [in] gctx The geometry context for this building call
0252   /// @param [in] surfaces is the vector of pointers to sensitive surfaces
0253   /// to be ordered on the plane
0254   /// @pre the pointers to the sensitive surfaces in the surfaces vectors all
0255   /// need to be valid, since no check is performed
0256   /// @warning This function requires the plane aligned with either the x-, y-
0257   /// or z-axis
0258   /// @param [in] bins1 is the number of bins in the orthogonal direction to @p
0259   /// aDir
0260   /// @param [in] bins2 is the number of bins in the orthogonal direction to @p
0261   /// aDir
0262   /// @param [in] aDir Direction of the aligned surfaces
0263   /// @param [in] protoLayerOpt Optional @c ProtoLayer instance
0264   /// @param [in] transform is the (optional) additional transform applied
0265   /// @param maxNeighborDistance Maximum next neighbor distance to be included in neighbor lookups
0266   ///
0267   /// @return a unique pointer a new SurfaceArray
0268   std::unique_ptr<SurfaceArray> surfaceArrayOnPlane(
0269       const GeometryContext& gctx,
0270       std::vector<std::shared_ptr<const Surface>> surfaces, std::size_t bins1,
0271       std::size_t bins2, AxisDirection aDir,
0272       std::optional<ProtoLayer> protoLayerOpt = std::nullopt,
0273       const Transform3& transform = Transform3::Identity(),
0274       std::uint8_t maxNeighborDistance = 1) const;
0275 
0276   /// Static check function for surface equivalent
0277   ///
0278   /// @param [in] gctx the geometry context for this check
0279   /// @param aDir the axis direction for the binning
0280   /// @param a first surface for checking
0281   /// @param b second surface for checking
0282   /// @return true if surfaces are equivalent for binning purposes
0283   static bool isSurfaceEquivalent(const GeometryContext& gctx,
0284                                   AxisDirection aDir, const Surface* a,
0285                                   const Surface* b) {
0286     using namespace UnitLiterals;
0287     using VectorHelpers::perp;
0288 
0289     if (aDir == AxisDirection::AxisPhi) {
0290       // Take the two binning positions
0291       Vector3 pos1 = a->referencePosition(gctx, AxisDirection::AxisR);
0292       Vector3 pos2 = b->referencePosition(gctx, AxisDirection::AxisR);
0293 
0294       // Project them on the (x, y) plane, where Phi angles are calculated
0295       auto proj1 = pos1.head<2>(), proj2 = pos2.head<2>();
0296 
0297       // Basic dot and cross products identities give us the cosine and sine
0298       // of these angles, time the squared vector norm
0299       auto cos_dPhi_n2 = proj1.dot(proj2);
0300       auto sin_dPhi_n2 = proj1.x() * proj2.y() - proj2.x() * proj1.y();
0301 
0302       // ...so by injecting them into atan2, we get the angle between them
0303       auto dPhi = std::atan2(sin_dPhi_n2, cos_dPhi_n2);
0304       return std::abs(dPhi) < std::numbers::pi / 180.;
0305     }
0306 
0307     if (aDir == AxisDirection::AxisZ) {
0308       return (std::abs(a->referencePosition(gctx, AxisDirection::AxisR).z() -
0309                        b->referencePosition(gctx, AxisDirection::AxisR).z()) <
0310               1_um);
0311     }
0312 
0313     if (aDir == AxisDirection::AxisR) {
0314       return (std::abs(perp(a->referencePosition(gctx, AxisDirection::AxisR)) -
0315                        perp(b->referencePosition(gctx, AxisDirection::AxisR))) <
0316               1_um);
0317     }
0318 
0319     return false;
0320   }
0321 
0322   /// Set logging instance
0323   /// @param logger is the logging instance to be set
0324   void setLogger(std::unique_ptr<const Logger> logger) {
0325     m_logger = std::move(logger);
0326   }
0327 
0328  private:
0329   /// configuration object
0330   Config m_cfg;
0331 
0332   /// Private access to logger
0333   const Logger& logger() const { return *m_logger; }
0334 
0335   std::vector<const Surface*> findKeySurfaces(
0336       const std::vector<const Surface*>& surfaces,
0337       const std::function<bool(const Surface*, const Surface*)>& equal) const;
0338 
0339   std::size_t determineBinCount(const GeometryContext& gctx,
0340                                 const std::vector<const Surface*>& surfaces,
0341                                 AxisDirection aDir) const;
0342 
0343   /// SurfaceArrayCreator internal method
0344   /// Creates a variable @c ProtoAxis from a vector of (unsorted) surfaces with
0345   /// PlanarBounds
0346   /// It loops through the surfaces and finds out the needed information
0347   /// First the surfaces are sorted in the binning direction and the so called
0348   /// "key" surfaces (surfaces with different positions in the binning
0349   /// direction) are extracted. The boundary value between two surfaces is the
0350   /// mean value of the two center position in the binning direction. The
0351   /// first and the last boundaries are calculated from the vertices of the
0352   /// first and last surface.
0353   /// @note currently implemented for phi, r and z bining
0354   /// @todo implement for x,y binning
0355   /// @param [in] gctx the geometry context for this call
0356   /// @param surfaces are the sensitive surfaces to be
0357   /// @param aDir the AxisDirection in which direction should be binned
0358   /// (currently possible: AxisPhi, AxisR, AxisZ)
0359   /// @param protoLayer Instance of @c ProtoLayer holding generic layer info
0360   /// @param transform is the (optional) additional transform applied
0361   /// @return Instance of @c ProtoAxis containing determined properties
0362   /// @note This only creates the @c ProtoAxis, this needs to be turned
0363   ///       into an actual @c Axis object to be used
0364   ProtoAxis createVariableAxis(const GeometryContext& gctx,
0365                                const std::vector<const Surface*>& surfaces,
0366                                AxisDirection aDir, const ProtoLayer& protoLayer,
0367                                Transform3& transform) const;
0368 
0369   /// SurfaceArrayCreator internal method
0370   /// Creates a equidistant @c ProtoAxis when the extrema and the bin number
0371   /// are
0372   /// It loops through the surfaces and finds out the needed information
0373   /// First the surfaces are sorted in the binning direction and the so called
0374   /// "key" surfaces (surfaces with different positions in the binning
0375   /// direction) are extracted. The number of key surfaces equals the number
0376   /// of bins. Afterwards the minimum and maximum are calculated by
0377   /// subtracting/adding half of a bin size to the center position (in the
0378   /// binning direction) to the first/last surface.
0379   /// @note currently implemented for phi, r and z bining
0380   /// @todo implement for x,y binning
0381   /// @param [in] gctx the geometry context for this call
0382   /// @param surfaces are the sensitive surfaces to be
0383   /// @param aDir the AxisDirection in which direction should be binned
0384   /// (currently possible: AxisPhi, AxisR, AxisZ)
0385   /// @param protoLayer Instance of @c ProtoLayer holding generic layer info
0386   /// @param transform is the (optional) additional transform applied
0387   /// @param nBins Number of bins to use, 0 means determine automatically
0388   /// @return Instance of @c ProtoAxis containing determined properties
0389   /// @note This only creates the @c ProtoAxis, this needs to be turned
0390   ///       into an actual @c Axis object to be used
0391   ProtoAxis createEquidistantAxis(const GeometryContext& gctx,
0392                                   const std::vector<const Surface*>& surfaces,
0393                                   AxisDirection aDir,
0394                                   const ProtoLayer& protoLayer,
0395                                   Transform3& transform,
0396                                   std::size_t nBins = 0) const;
0397 
0398   /// SurfaceArrayCreator internal method
0399   /// @brief Creates a SurfaceGridLookup instance within an any
0400   /// This is essentially a factory which absorbs some if/else logic
0401   /// that is required by the templating.
0402   /// @tparam bdtA AxisBoundaryType of axis A
0403   /// @tparam bdtB AxisBoundaryType of axis B
0404   /// @param surface the surface of the grid
0405   /// @param layerTolerance the layer tolerance
0406   /// @param pAxisA ProtoAxis object for axis A
0407   /// @param pAxisB ProtoAxis object for axis B
0408   /// @param maxNeighborDistance the maximum neighbor distance for the grid lookup
0409   template <AxisBoundaryType bdtA, AxisBoundaryType bdtB>
0410   static std::unique_ptr<SurfaceArray::ISurfaceGridLookup>
0411   makeSurfaceGridLookup2D(std::shared_ptr<RegularSurface> surface,
0412                           double layerTolerance, const ProtoAxis& pAxisA,
0413                           const ProtoAxis& pAxisB,
0414                           std::uint8_t maxNeighborDistance) {
0415     using ISGL = SurfaceArray::ISurfaceGridLookup;
0416     std::unique_ptr<ISGL> ptr;
0417 
0418     if (pAxisA.bType == equidistant && pAxisB.bType == equidistant) {
0419       Axis<AxisType::Equidistant, bdtA> axisA(pAxisA.min, pAxisA.max,
0420                                               pAxisA.nBins);
0421       Axis<AxisType::Equidistant, bdtB> axisB(pAxisB.min, pAxisB.max,
0422                                               pAxisB.nBins);
0423 
0424       ptr = SurfaceArray::makeSurfaceGridLookup(
0425           std::move(surface), layerTolerance, std::tuple{axisA, axisB},
0426           maxNeighborDistance);
0427 
0428     } else if (pAxisA.bType == equidistant && pAxisB.bType == arbitrary) {
0429       Axis<AxisType::Equidistant, bdtA> axisA(pAxisA.min, pAxisA.max,
0430                                               pAxisA.nBins);
0431       Axis<AxisType::Variable, bdtB> axisB(pAxisB.binEdges);
0432 
0433       ptr = SurfaceArray::makeSurfaceGridLookup(
0434           std::move(surface), layerTolerance, std::tuple{axisA, axisB},
0435           maxNeighborDistance);
0436 
0437     } else if (pAxisA.bType == arbitrary && pAxisB.bType == equidistant) {
0438       Axis<AxisType::Variable, bdtA> axisA(pAxisA.binEdges);
0439       Axis<AxisType::Equidistant, bdtB> axisB(pAxisB.min, pAxisB.max,
0440                                               pAxisB.nBins);
0441 
0442       ptr = SurfaceArray::makeSurfaceGridLookup(
0443           std::move(surface), layerTolerance, std::tuple{axisA, axisB},
0444           maxNeighborDistance);
0445 
0446     } else /*if (pAxisA.bType == arbitrary && pAxisB.bType == arbitrary)*/ {
0447       Axis<AxisType::Variable, bdtA> axisA(pAxisA.binEdges);
0448       Axis<AxisType::Variable, bdtB> axisB(pAxisB.binEdges);
0449 
0450       ptr = SurfaceArray::makeSurfaceGridLookup(
0451           std::move(surface), layerTolerance, std::tuple{axisA, axisB},
0452           maxNeighborDistance);
0453     }
0454 
0455     return ptr;
0456   }
0457 
0458   /// logging instance
0459   std::unique_ptr<const Logger> m_logger;
0460 
0461   /// Private helper method to transform the  vertices of surface bounds into
0462   /// global coordinates
0463   /// @param [in] gctx the geometry context for this call
0464   /// @param surface the surface associated with the given vertices
0465   /// @param locVertices a vector of the vertices in local coordinates
0466   /// @return a vector of the vertices in global coordinates
0467   std::vector<Vector3> makeGlobalVertices(
0468       const GeometryContext& gctx, const Surface& surface,
0469       const std::vector<Vector2>& locVertices) const;
0470 };
0471 
0472 }  // namespace Acts