Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-12-03 09:14:34

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