Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-06-30 07:51:45

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