Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-05-29 07:46:32

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/ReferenceGenerators.hpp"
0013 #include "Acts/Surfaces/Surface.hpp"
0014 #include "Acts/Utilities/AxisDefinitions.hpp"
0015 #include "Acts/Utilities/KDTree.hpp"
0016 
0017 #include <array>
0018 #include <ranges>
0019 #include <stdexcept>
0020 #include <tuple>
0021 #include <vector>
0022 
0023 namespace Acts {
0024 
0025 /// @brief A wrapper class around a KDTree of surfaces
0026 ///
0027 /// It also deals with the conversion from global query to
0028 /// KDTree lookup positions
0029 ///
0030 template <std::size_t kDIM = 2u, std::size_t bSize = 100u>
0031 class KdtSurfaces {
0032  public:
0033   /// Broadcast the surface KDT type
0034   using KDTS =
0035       KDTree<kDIM, std::shared_ptr<Surface>, double, std::array, bSize>;
0036 
0037   /// Broadcast the query definition
0038   using Query = std::array<double, kDIM>;
0039 
0040   /// Broadcast the entry
0041   using Entry = std::pair<Query, std::shared_ptr<Surface>>;
0042 
0043   /// Constructor from a vector of surfaces
0044   ///
0045   /// @param gctx the geometry context of this call
0046   /// @param surfaces the surfaces to be filled into the tree
0047   /// @param casts the cast list from global position into kdtree local
0048   /// @param rgen the reference point generator
0049   KdtSurfaces(const GeometryContext& gctx,
0050               const std::vector<std::shared_ptr<Surface>>& surfaces,
0051               const std::array<AxisDirection, kDIM>& casts,
0052               const std::shared_ptr<IReferenceGenerator>& rgen =
0053                   std::make_shared<PolyhedronReferenceGenerator>())
0054       : m_kdt(nullptr), m_casts(casts), m_rGenerator(rgen) {
0055     // Simple check if the dimension is correct
0056     if (kDIM == 0u) {
0057       throw std::invalid_argument(
0058           "KdtSurfaces: dimension and/or cast rules are incorrect.");
0059     }
0060     // Fill the tree from surfaces
0061     std::vector<Entry> kdtEntries;
0062     kdtEntries.reserve(surfaces.size());
0063     for (auto& s : surfaces) {
0064       // Generate the references and the center of gravity from it
0065       const auto references = m_rGenerator->references(gctx, *s);
0066       std::vector<Query> castedReferences;
0067       castedReferences.reserve(references.size());
0068       for (const auto& r : references) {
0069         //  Now cast into the correct fill position
0070         Query rc = {};
0071         fillCasts(r, rc, std::make_integer_sequence<std::size_t, kDIM>{});
0072         castedReferences.push_back(rc);
0073       }
0074       // Calculate the center of gravity in casted frame
0075       kdtEntries.push_back({cog(castedReferences), s});
0076     }
0077     // Create the KDTree
0078     m_kdt = std::make_unique<KDTS>(std::move(kdtEntries));
0079   }
0080 
0081   /// Query with a Range object
0082   ///
0083   /// @param range is the range to be queried
0084   ///
0085   /// @return the matching surfaces from the KDT structure
0086   std::vector<std::shared_ptr<Surface>> surfaces(
0087       const RangeXD<kDIM, double>& range) const {
0088     // Strip the surfaces
0089     std::vector<std::shared_ptr<Surface>> surfacePtrs;
0090     auto surfaceQuery = m_kdt->rangeSearchWithKey(range);
0091     std::ranges::for_each(
0092         surfaceQuery, [&](auto& surf) { surfacePtrs.push_back(surf.second); });
0093     return surfacePtrs;
0094   }
0095 
0096   /// Query with an Extent object
0097   ///
0098   /// @param extent is the range Extent to be queried
0099   ///
0100   /// @return the matching surfaces fpulled from the KDT structure
0101   std::vector<std::shared_ptr<Surface>> surfaces(const Extent& extent) const {
0102     RangeXD<kDIM, double> qRange;
0103     for (auto [ibv, v] : enumerate(m_casts)) {
0104       qRange[ibv] = extent.range(v);
0105     }
0106     return surfaces(qRange);
0107   }
0108 
0109  private:
0110   /// The KDTree as single source for the surfaces (maybe shared)
0111   std::unique_ptr<KDTS> m_kdt = nullptr;
0112 
0113   /// Cast values that turn a global position to lookup position
0114   std::array<AxisDirection, kDIM> m_casts = {};
0115 
0116   /// Helper to generate reference points for filling
0117   std::shared_ptr<IReferenceGenerator> m_rGenerator{nullptr};
0118 
0119   /// Unroll the cast loop
0120   /// @param position is the position of the update call
0121   /// @param a is the array to be filled
0122   template <typename Array, std::size_t... idx>
0123   void fillCasts(const Vector3& position, Array& a,
0124                  std::index_sequence<idx...> /*indices*/) const {
0125     ((a[idx] = VectorHelpers::cast(position, m_casts[idx])), ...);
0126   }
0127 
0128   /// Helper method to calculate the center of gravity in the
0129   /// casted frame (i.e. query frame)
0130   ///
0131   /// @param cQueries are the casted query positions
0132   /// @note will do nothing if vector size is equal to 1
0133   ///
0134   /// @note no checking on qQueries.empty() is done as the
0135   /// positions are to be provided by a generator which itself
0136   /// is tested for consistency
0137   ///
0138   /// @return the center of gravity as a query object
0139   Query cog(const std::vector<Query>& cQueries) const {
0140     // If there is only one position, return it
0141     if (cQueries.size() == 1) {
0142       return cQueries.front();
0143     }
0144     // Build the center of gravity of the n positions
0145     Query c{};
0146     float weight = 1. / cQueries.size();
0147     for (auto& q : cQueries) {
0148       std::transform(c.begin(), c.end(), q.begin(), c.begin(),
0149                      std::plus<double>());
0150     }
0151     std::ranges::for_each(c, [&](auto& v) { v *= weight; });
0152     return c;
0153   }
0154 };
0155 
0156 }  // namespace Acts