Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-09-30 08:01:40

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 #include "Acts/Navigation/NavigationStream.hpp"
0010 
0011 #include "Acts/Propagator/NavigationTarget.hpp"
0012 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0013 #include "Acts/Surfaces/Surface.hpp"
0014 #include "Acts/Utilities/Enumerate.hpp"
0015 #include "Acts/Utilities/StringHelpers.hpp"
0016 
0017 #include <algorithm>
0018 
0019 namespace Acts {
0020 
0021 bool NavigationStream::initialize(const GeometryContext& gctx,
0022                                   const QueryPoint& queryPoint,
0023                                   const Logger& logger,
0024                                   const double onSurfaceTolerance,
0025                                   const bool candidatesAreUnique) {
0026   // Position and direction from the query point
0027   const Vector3& position = queryPoint.position;
0028   const Vector3& direction = queryPoint.direction;
0029 
0030   // De-duplicate by surface pointer. The first entry wins, so its tolerance is
0031   // the one used.
0032   if (!candidatesAreUnique) {
0033     ACTS_VERBOSE("De-duplicate the candidates:" << m_candidates);
0034     std::size_t writeIdx = 0;
0035     for (std::size_t readIdx = 0; readIdx < m_candidates.size(); ++readIdx) {
0036       const Surface* surface = &m_candidates[readIdx].surface();
0037       bool alreadySeen = false;
0038       for (std::size_t k = 0; k < writeIdx; ++k) {
0039         if (&m_candidates[k].surface() == surface) {
0040           alreadySeen = true;
0041           break;
0042         }
0043       }
0044       if (!alreadySeen) {
0045         if (writeIdx != readIdx) {
0046           m_candidates[writeIdx] = m_candidates[readIdx];
0047         }
0048         ++writeIdx;
0049       }
0050     }
0051     m_candidates.erase(m_candidates.begin() + writeIdx, m_candidates.end());
0052   }
0053 
0054   // Collect additional candidates for the second valid intersection. Reuse the
0055   // member scratch buffer to avoid a heap allocation on every call.
0056   std::vector<NavigationTarget>& additionalCandidates = m_additionalCandidates;
0057   additionalCandidates.clear();
0058   for (auto& candidate : m_candidates) {
0059     // Get the surface from the object intersection
0060     const Surface& surface = candidate.surface();
0061     // Intersect the surface
0062     auto multiIntersection =
0063         surface.intersect(gctx, position, direction,
0064                           candidate.boundaryTolerance(), onSurfaceTolerance);
0065 
0066     bool firstValid = multiIntersection.at(0).isValid();
0067     bool secondValid = multiIntersection.at(1).isValid();
0068     if (firstValid && !secondValid) {
0069       if (multiIntersection.at(0).pathLength() < -onSurfaceTolerance) {
0070         continue;
0071       }
0072       candidate.intersection() = multiIntersection.at(0);
0073       candidate.intersectionIndex() = 0;
0074     } else if (!firstValid && secondValid) {
0075       if (multiIntersection.at(1).pathLength() < -onSurfaceTolerance) {
0076         continue;
0077       }
0078       candidate.intersection() = multiIntersection.at(1);
0079       candidate.intersectionIndex() = 1;
0080     } else {
0081       // Split them into valid intersections, keep track of potentially
0082       // additional candidates
0083       bool originalCandidateUpdated = false;
0084       for (auto [intersectionIndex, intersection] :
0085            enumerate(multiIntersection)) {
0086         // Skip negative solutions, respecting the on surface tolerance
0087         if (intersection.pathLength() < -onSurfaceTolerance) {
0088           continue;
0089         }
0090         // Valid solution is either on surface or updates the distance
0091         if (intersection.isValid()) {
0092           if (!originalCandidateUpdated) {
0093             candidate.intersection() = intersection;
0094             candidate.intersectionIndex() = intersectionIndex;
0095             originalCandidateUpdated = true;
0096           } else {
0097             NavigationTarget additionalCandidate = candidate;
0098             additionalCandidate.intersection() = intersection;
0099             additionalCandidate.intersectionIndex() = intersectionIndex;
0100             additionalCandidates.emplace_back(additionalCandidate);
0101           }
0102         }
0103       }
0104     }
0105   }
0106 
0107   // Append the multi intersection candidates
0108   m_candidates.insert(m_candidates.end(), additionalCandidates.begin(),
0109                       additionalCandidates.end());
0110 
0111   // Sort the candidates by path length
0112   std::ranges::sort(m_candidates, NavigationTarget::pathLengthOrder);
0113   ACTS_VERBOSE("Sorted candidates:" << m_candidates);
0114 
0115   // If we have duplicates, we expect them to be close by in path length, so we
0116   // don't need to re-sort Remove duplicates on basis of the surface pointer
0117 
0118   /// But but but... What about the surfaces with multiple intersections?
0119   auto nonUniqueRange = std::ranges::unique(
0120       m_candidates.begin(), m_candidates.end(),
0121       [](const NavigationTarget& a, const NavigationTarget& b) {
0122         return &a.surface() == &b.surface();
0123       });
0124   m_candidates.erase(nonUniqueRange.begin(), nonUniqueRange.end());
0125 
0126   // The we find the first invalid candidate
0127   auto firstInvalid = std::ranges::find_if(
0128       m_candidates,
0129       [](const NavigationTarget& a) { return !a.intersection().isValid(); });
0130 
0131   // Set the range and initialize
0132   m_candidates.resize(std::distance(m_candidates.begin(), firstInvalid),
0133                       NavigationTarget::None());
0134 
0135   m_currentIndex = 0;
0136   return isValid();
0137 }
0138 
0139 bool NavigationStream::update(const GeometryContext& gctx,
0140                               const QueryPoint& queryPoint,
0141                               const Logger& logger, double onSurfaceTolerance) {
0142   ACTS_VERBOSE("Update from position " << toString(queryPoint.position)
0143                                        << " and direction "
0144                                        << toString(queryPoint.direction));
0145   // Loop over the (currently valid) candidates and update
0146   for (; m_currentIndex < m_candidates.size(); ++m_currentIndex) {
0147     // Get the candidate, and resolve the tuple
0148     NavigationTarget& candidate = currentCandidate();
0149     // Get the surface from the object intersection
0150     const Surface& surface = candidate.surface();
0151     // (re-)Intersect the surface
0152     auto multiIntersection =
0153         surface.intersect(gctx, queryPoint.position, queryPoint.direction,
0154                           candidate.boundaryTolerance(), onSurfaceTolerance);
0155     // Split them into valid intersections
0156     for (auto [intersectionIndex, intersection] :
0157          enumerate(multiIntersection)) {
0158       // Skip wrong index solution
0159       if (intersectionIndex != candidate.intersectionIndex()) {
0160         continue;
0161       }
0162       // Valid solution is either on surface or updates the distance
0163       if (intersection.isValid()) {
0164         candidate.intersection() = intersection;
0165         ACTS_VERBOSE("Updated candidate " << candidate);
0166         return true;
0167       }
0168     }
0169   }
0170   // No candidate was reachable
0171   return false;
0172 }
0173 
0174 void NavigationStream::reset() {
0175   m_candidates.clear();
0176   m_currentIndex = 0;
0177 }
0178 
0179 void NavigationStream::addSurfaceCandidate(
0180     const Surface& surface, const BoundaryTolerance& bTolerance) {
0181   m_candidates.emplace_back(Intersection3D::Invalid(), 0, surface, bTolerance);
0182 }
0183 
0184 void NavigationStream::addSurfaceCandidates(
0185     std::span<const Surface*> surfaces, const BoundaryTolerance& bTolerance) {
0186   m_candidates.reserve(m_candidates.size() + surfaces.size());
0187   std::ranges::for_each(surfaces, [&](const Surface* surface) {
0188     m_candidates.emplace_back(Intersection3D::Invalid(), 0, *surface,
0189                               bTolerance);
0190   });
0191 }
0192 
0193 void NavigationStream::addPortalCandidate(const Portal& portal) {
0194   m_candidates.emplace_back(Intersection3D::Invalid(), 0, portal,
0195                             BoundaryTolerance::None());
0196 }
0197 
0198 AppendOnlyNavigationStream::AppendOnlyNavigationStream(NavigationStream& stream)
0199     : m_stream{&stream} {}
0200 
0201 void AppendOnlyNavigationStream::addPortalCandidate(const Portal& portal) {
0202   m_stream->addPortalCandidate(portal);
0203 }
0204 
0205 void AppendOnlyNavigationStream::addSurfaceCandidate(
0206     const Surface& surface, const BoundaryTolerance& bTolerance) {
0207   m_stream->addSurfaceCandidate(surface, bTolerance);
0208 }
0209 
0210 }  // namespace Acts