Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-07-26 08:18:25

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/Direction.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/EventData/SpacePointContainer.hpp"
0014 #include "Acts/Utilities/Delegate.hpp"
0015 #include "Acts/Utilities/detail/ContainerIterator.hpp"
0016 
0017 #include <cstdint>
0018 #include <vector>
0019 
0020 namespace Acts {
0021 
0022 /// Container for doublets found by the doublet seed finder.
0023 ///
0024 /// This implementation uses partial AoS/SoA depending on the access pattern in
0025 /// the doublet finding process.
0026 class DoubletsForMiddleSp {
0027  public:
0028   /// Type alias for index type used in doublets container
0029   using Index = std::uint32_t;
0030   /// Type alias for range of indices in doublets container
0031   using IndexRange = std::pair<Index, Index>;
0032   /// Type alias for subset of indices in doublets container
0033   using IndexSubset = std::span<const Index>;
0034 
0035   /// Check if the doublets container is empty
0036   /// @return True if container has no doublets
0037   [[nodiscard]] bool empty() const { return m_spacePoints.empty(); }
0038   /// Get the number of doublets in container
0039   /// @return Number of doublets stored
0040   [[nodiscard]] Index size() const {
0041     return static_cast<Index>(m_spacePoints.size());
0042   }
0043 
0044   /// Clear all stored doublets and associated data
0045   void clear() {
0046     m_spacePoints.clear();
0047     m_cotTheta.clear();
0048     m_er_iDeltaR.clear();
0049     m_uv.clear();
0050     m_xy.clear();
0051   }
0052 
0053   /// Add a new doublet with associated parameters
0054   /// @param sp space point index for the doublet
0055   /// @param cotTheta Cotangent of polar angle
0056   /// @param iDeltaR Inverse delta R parameter
0057   /// @param er Error in R coordinate
0058   /// @param u U coordinate parameter
0059   /// @param v V coordinate parameter
0060   /// @param x X coordinate
0061   /// @param y Y coordinate
0062   void emplace_back(SpacePointIndex sp, float cotTheta, float iDeltaR, float er,
0063                     float u, float v, float x, float y) {
0064     m_spacePoints.push_back(sp);
0065     m_cotTheta.push_back(cotTheta);
0066     m_er_iDeltaR.push_back({er, iDeltaR});
0067     m_uv.push_back({u, v});
0068     m_xy.push_back({x, y});
0069   }
0070 
0071   /// Get reference to space point indices container
0072   /// @return Const reference to space point indices vector
0073   const std::vector<SpacePointIndex>& spacePoints() const {
0074     return m_spacePoints;
0075   }
0076   /// Get reference to cotTheta values container
0077   /// @return Const reference to cotTheta values vector
0078   const std::vector<float>& cotTheta() const { return m_cotTheta; }
0079 
0080   /// Pair of doublet index and cotTheta value.
0081   struct IndexAndCotTheta {
0082     /// Doublet index
0083     Index index{};
0084     /// Cotangent of theta
0085     float cotTheta{};
0086   };
0087 
0088   /// Type alias for subset of index and cotTheta pairs
0089   using IndexAndCotThetaSubset = std::span<const IndexAndCotTheta>;
0090 
0091   /// Sort doublets by cotTheta within given range
0092   /// @param range Index range to sort within
0093   /// @param indexAndCotTheta Output vector containing sorted index and cotTheta pairs
0094   void sortByCotTheta(const IndexRange& range,
0095                       std::vector<IndexAndCotTheta>& indexAndCotTheta) const {
0096     indexAndCotTheta.clear();
0097     indexAndCotTheta.reserve(range.second - range.first);
0098     for (Index i = range.first; i < range.second; ++i) {
0099       indexAndCotTheta.emplace_back(i, m_cotTheta[i]);
0100     }
0101     std::ranges::sort(indexAndCotTheta, {}, [](const IndexAndCotTheta& item) {
0102       return item.cotTheta;
0103     });
0104   }
0105 
0106   /// Proxy accessor for a single doublet entry.
0107   class Proxy {
0108    public:
0109     /// Constructor
0110     /// @param container Pointer to the doublet container
0111     /// @param index Index of the doublet
0112     Proxy(const DoubletsForMiddleSp* container, Index index)
0113         : m_container(container), m_index(index) {}
0114 
0115     /// Get the doublet container
0116     /// @return Reference to the container
0117     const DoubletsForMiddleSp& container() const { return *m_container; }
0118     /// Get the doublet index
0119     /// @return The index
0120     Index index() const { return m_index; }
0121 
0122     /// Get space point index pair
0123     /// @return The space point index
0124     SpacePointIndex spacePointIndex() const {
0125       return m_container->m_spacePoints[m_index];
0126     }
0127 
0128     /// Get cotangent of theta
0129     /// @return The cotTheta value
0130     float cotTheta() const { return m_container->m_cotTheta[m_index]; }
0131     /// Get er value
0132     /// @return The er value
0133     float er() const { return m_container->m_er_iDeltaR[m_index][0]; }
0134     /// Get inverse delta r
0135     /// @return The inverse delta r value
0136     float iDeltaR() const { return m_container->m_er_iDeltaR[m_index][1]; }
0137     /// Get u coordinate
0138     /// @return The u value
0139     float u() const { return m_container->m_uv[m_index][0]; }
0140     /// Get v coordinate
0141     /// @return The v value
0142     float v() const { return m_container->m_uv[m_index][1]; }
0143     /// Get x coordinate
0144     /// @return The x value
0145     float x() const { return m_container->m_xy[m_index][0]; }
0146     /// Get y coordinate
0147     /// @return The y value
0148     float y() const { return m_container->m_xy[m_index][1]; }
0149 
0150    private:
0151     const DoubletsForMiddleSp* m_container{};
0152     Index m_index{};
0153   };
0154   /// Same as `Proxy` but also contains `cotTheta`. This is useful after sorting
0155   /// doublets by `cotTheta` to avoid indirect access.
0156   class Proxy2 : public Proxy {
0157    public:
0158     /// Constructor for Proxy2 with precomputed cotTheta
0159     /// @param container Pointer to the doublets container
0160     /// @param indexAndCotTheta Index and cotTheta pair
0161     Proxy2(const DoubletsForMiddleSp* container,
0162            IndexAndCotTheta indexAndCotTheta)
0163         : Proxy(container, indexAndCotTheta.index),
0164           m_cotTheta(indexAndCotTheta.cotTheta) {}
0165 
0166     /// Get precomputed cotTheta value (avoids indirect access)
0167     /// @return Cotangent of polar angle
0168     float cotTheta() const { return m_cotTheta; }
0169 
0170    private:
0171     float m_cotTheta{};
0172   };
0173 
0174   /// Access doublet by index
0175   /// @param index Index of the doublet to access
0176   /// @return Proxy object for the doublet
0177   Proxy operator[](Index index) const { return Proxy(this, index); }
0178   /// Access doublet by index and cotTheta pair
0179   /// @param indexAndCotTheta Index and cotTheta pair for the doublet
0180   /// @return Proxy2 object for the doublet with precomputed cotTheta
0181   Proxy2 operator[](IndexAndCotTheta indexAndCotTheta) const {
0182     return Proxy2(this, indexAndCotTheta);
0183   }
0184 
0185   /// Type alias for const iterator over doublets in container
0186   using const_iterator =
0187       detail::ContainerIterator<DoubletsForMiddleSp, Proxy, Index, true>;
0188 
0189   /// Get iterator to beginning of doublets container
0190   /// @return Const iterator to first doublet
0191   const_iterator begin() const { return const_iterator(*this, 0); }
0192   /// Get iterator to end of doublets container
0193   /// @return Const iterator past the last doublet
0194   const_iterator end() const { return const_iterator(*this, size()); }
0195 
0196   /// Range view over doublets in the container.
0197   class Range : public detail::ContainerRange<Range, Range, DoubletsForMiddleSp,
0198                                               Index, true> {
0199    public:
0200     /// Base class type alias
0201     using Base =
0202         detail::ContainerRange<Range, Range, DoubletsForMiddleSp, Index, true>;
0203 
0204     using Base::Base;
0205   };
0206 
0207   /// Get range view of all doublets
0208   /// @return Range object covering all doublets
0209   Range range() const noexcept { return Range(*this, {0, size()}); }
0210   /// Get range view of doublets within specified index range
0211   /// @param range Index range to create view for
0212   /// @return Range object covering specified doublets
0213   Range range(const IndexRange& range) const noexcept {
0214     return Range(*this, range);
0215   }
0216 
0217   /// Subset view of doublets addressed by indices.
0218   class Subset
0219       : public detail::ContainerSubset<Subset, Subset, DoubletsForMiddleSp,
0220                                        Proxy, std::span<const Index>, true> {
0221    public:
0222     /// Base class type alias
0223     using Base = detail::ContainerSubset<Subset, Subset, DoubletsForMiddleSp,
0224                                          Proxy, std::span<const Index>, true>;
0225 
0226     using Base::Base;
0227   };
0228   /// Subset view of doublets addressed by index and cotTheta pairs.
0229   class Subset2 : public detail::ContainerSubset<
0230                       Subset2, Subset2, DoubletsForMiddleSp, Proxy2,
0231                       std::span<const IndexAndCotTheta>, true> {
0232    public:
0233     /// Base class type alias
0234     using Base =
0235         detail::ContainerSubset<Subset2, Subset2, DoubletsForMiddleSp, Proxy2,
0236                                 std::span<const IndexAndCotTheta>, true>;
0237 
0238     using Base::Base;
0239   };
0240 
0241   /// Create subset view from index subset
0242   /// @param subset Span of indices to include in subset
0243   /// @return Subset object for the specified indices
0244   Subset subset(const IndexSubset& subset) const noexcept {
0245     return Subset(*this, subset);
0246   }
0247   /// Create subset view from index and cotTheta subset
0248   /// @param subset Span of index and cotTheta pairs to include
0249   /// @return Subset2 object with precomputed cotTheta values
0250   Subset2 subset(const IndexAndCotThetaSubset& subset) const noexcept {
0251     return Subset2(*this, subset);
0252   }
0253 
0254  private:
0255   std::vector<SpacePointIndex> m_spacePoints;
0256 
0257   // parameters required to calculate a circle with linear equation
0258   std::vector<float> m_cotTheta;
0259   std::vector<std::array<float, 2>> m_er_iDeltaR;
0260   std::vector<std::array<float, 2>> m_uv;
0261   std::vector<std::array<float, 2>> m_xy;
0262 };
0263 
0264 /// Derived quantities for the middle space point in a doublet.
0265 struct MiddleSpInfo {
0266   /// minus one over radius of middle SP
0267   float uIP{};
0268   /// square of uIP
0269   float uIP2{};
0270   /// ratio between middle SP x position and radius
0271   float cosPhiM{};
0272   /// ratio between middle SP y position and radius
0273   float sinPhiM{};
0274 };
0275 
0276 /// Interface and a collection of standard implementations for a doublet seed
0277 /// finder. Given a starting space point and a collection of candidates, it
0278 /// finds all doublets that satisfy the selection criteria. For the standard
0279 /// implementations the criteria are given by interaction point cuts.
0280 ///
0281 /// @note The standard implementations rely on virtual function dispatch which
0282 /// did not turn out to affect the performance after measurement.
0283 class DoubletSeedFinder {
0284  public:
0285   /// Collection of configuration parameters for the doublet seed finder. This
0286   /// includes doublet cuts, steering switches, and assumptions about the space
0287   /// points.
0288   struct Config {
0289     /// Whether the input space points are sorted by radius
0290     bool spacePointsSortedByRadius = false;
0291 
0292     /// Direction of the doublet candidate space points. Either forward, also
0293     /// called top doublet, or backward, also called bottom doublet.
0294     Direction candidateDirection = Direction::Forward();
0295 
0296     /// Minimum radial distance between two doublet components
0297     float deltaRMin = 5 * UnitConstants::mm;
0298     /// Maximum radial distance between two doublet components
0299     float deltaRMax = 270 * UnitConstants::mm;
0300 
0301     /// Minimal z distance between two doublet components
0302     float deltaZMin = -std::numeric_limits<float>::infinity();
0303     /// Maximum z distance between two doublet components
0304     float deltaZMax = std::numeric_limits<float>::infinity();
0305 
0306     /// Maximum value of impact parameter estimation of the seed candidates
0307     float impactMax = 20 * UnitConstants::mm;
0308 
0309     /// Enable cut on the compatibility between interaction point and doublet,
0310     /// this is an useful approximation to speed up the seeding
0311     bool interactionPointCut = false;
0312 
0313     /// Limiting location of collision region in z-axis used to check if doublet
0314     /// origin is within reasonable bounds
0315     float collisionRegionMin = -150 * UnitConstants::mm;
0316     /// Maximum collision region boundary in z-axis for doublet origin checks
0317     float collisionRegionMax = +150 * UnitConstants::mm;
0318 
0319     /// Maximum allowed cotTheta between two space-points in doublet, used to
0320     /// check if forward angle is within bounds
0321     float cotThetaMax = 10.01788;  // equivalent to eta = 3 (pseudorapidity)
0322 
0323     /// Minimum transverse momentum (pT) used to check the r-z slope
0324     /// compatibility of triplets with maximum multiple scattering effect
0325     /// (produced by the minimum allowed pT particle) + a certain uncertainty
0326     /// term. Check the documentation for more information
0327     /// https://acts.readthedocs.io/en/latest/core/reconstruction/pattern_recognition/seeding.html
0328     float minPt = 400 * UnitConstants::MeV;
0329     /// Parameter which can loosen the tolerance of the track seed to form a
0330     /// helix. This is useful for e.g. misaligned seeding.
0331     float helixCutTolerance = 1;
0332 
0333     /// Type alias for delegate to apply experiment specific cuts during doublet
0334     /// finding
0335     using ExperimentCuts =
0336         Delegate<bool(const ConstSpacePointProxy& /*middle*/,
0337                       const ConstSpacePointProxy& /*other*/, float /*cotTheta*/,
0338                       bool /*isBottomCandidate*/)>;
0339 
0340     /// Delegate to apply experiment specific cuts during doublet finding
0341     ExperimentCuts experimentCuts;
0342   };
0343 
0344   /// Derived configuration for the doublet seed finder using a magnetic field.
0345   struct DerivedConfig : public Config {
0346     /// Constructor for derived configuration with magnetic field
0347     /// @param config Base configuration to derive from
0348     /// @param bFieldInZ Magnetic field strength in z-direction
0349     DerivedConfig(const Config& config, float bFieldInZ);
0350 
0351     /// Magnetic field strength in z-direction for helix calculation
0352     float bFieldInZ = std::numeric_limits<float>::quiet_NaN();
0353     /// Squared minimum helix diameter derived from magnetic field and minimum
0354     /// pT
0355     float minHelixDiameter2 = std::numeric_limits<float>::quiet_NaN();
0356   };
0357 
0358   /// Computes additional quantities from the middle space point which can be
0359   /// reused during doublet finding.
0360   /// @param spM Middle space point for doublet computation
0361   /// @return MiddleSpInfo structure with computed quantities
0362   static MiddleSpInfo computeMiddleSpInfo(const ConstSpacePointProxy& spM);
0363 
0364   /// Creates a new doublet seed finder instance given the configuration.
0365   /// @param config Configuration for the doublet seed finder
0366   /// @return Unique pointer to new DoubletSeedFinder instance
0367   static std::unique_ptr<DoubletSeedFinder> create(const DerivedConfig& config);
0368 
0369   virtual ~DoubletSeedFinder() = default;
0370 
0371   /// Returns the configuration of the doublet seed finder.
0372   /// @return Reference to the configuration object
0373   virtual const DerivedConfig& config() const = 0;
0374 
0375   /// Creates compatible dublets by applying a series of cuts that can be
0376   /// tested with only two SPs.
0377   ///
0378   /// @param middleSp Space point candidate to be used as middle SP in a seed
0379   /// @param middleSpInfo Information about the middle space point
0380   /// @param candidateSps Subset of space points to be used as candidates for
0381   ///   middle SP in a seed
0382   /// @param compatibleDoublets Output container for compatible doublets
0383   virtual void createDoublets(
0384       const ConstSpacePointProxy& middleSp, const MiddleSpInfo& middleSpInfo,
0385       SpacePointContainer::ConstSubset& candidateSps,
0386       DoubletsForMiddleSp& compatibleDoublets) const = 0;
0387 
0388   /// Creates compatible dublets by applying a series of cuts that can be
0389   /// tested with only two SPs.
0390   ///
0391   /// @param middleSp Space point candidate to be used as middle SP in a seed
0392   /// @param middleSpInfo Information about the middle space point
0393   /// @param candidateSps Range of space points to be used as candidates for
0394   ///   middle SP in a seed
0395   /// @param compatibleDoublets Output container for compatible doublets
0396   virtual void createDoublets(
0397       const ConstSpacePointProxy& middleSp, const MiddleSpInfo& middleSpInfo,
0398       SpacePointContainer::ConstRange& candidateSps,
0399       DoubletsForMiddleSp& compatibleDoublets) const = 0;
0400 };
0401 
0402 }  // namespace Acts