diff --git a/Core/include/Acts/Propagator/NavigationTarget.hpp b/Core/include/Acts/Propagator/NavigationTarget.hpp index b95e5a7f074..47404508b18 100644 --- a/Core/include/Acts/Propagator/NavigationTarget.hpp +++ b/Core/include/Acts/Propagator/NavigationTarget.hpp @@ -189,7 +189,11 @@ class NavigationTarget { constexpr const BoundaryTolerance& boundaryTolerance() const noexcept { return m_boundaryTolerance; } - + /// Returns whether the target has been reached by the propagator. + /// @return the target passed flag + constexpr bool isReached() const noexcept { return m_reached; } + /// Marks the target as reached + constexpr void targetReached() noexcept { m_reached = true; } /// Returns whether the intersection was successful or not /// @return true if the intersection is valid constexpr bool isValid() const noexcept { return m_intersection.isValid(); } @@ -279,11 +283,13 @@ class NavigationTarget { const Surface* m_surfaceRepresentation = nullptr; /// The boundary tolerance used for this intersection BoundaryTolerance m_boundaryTolerance = BoundaryTolerance::None(); + /// State toggling whether the target has been reached by the propgator + bool m_reached{false}; /// Default constructor creating a none target constexpr NavigationTarget() = default; - /// @brief print method + /// Print method /// @param ostr: Stream to which the object is printed void print(std::ostream& ostr) const; }; diff --git a/Core/include/Acts/Propagator/Navigator.hpp b/Core/include/Acts/Propagator/Navigator.hpp index 633d46aee1a..e8e2e44d160 100644 --- a/Core/include/Acts/Propagator/Navigator.hpp +++ b/Core/include/Acts/Propagator/Navigator.hpp @@ -241,7 +241,6 @@ class Navigator final { navLayerIndex.reset(); navBoundaries.clear(); navBoundaryIndex.reset(); - navCandidates.clear(); navCandidateIndex.reset(); currentLayer = nullptr; diff --git a/Core/include/Acts/Propagator/NavigatorOptions.hpp b/Core/include/Acts/Propagator/NavigatorOptions.hpp index 31de4953eb3..28382e2b577 100644 --- a/Core/include/Acts/Propagator/NavigatorOptions.hpp +++ b/Core/include/Acts/Propagator/NavigatorOptions.hpp @@ -44,6 +44,10 @@ struct NavigatorPlainOptions { /// The far limit to resolve surfaces double farLimit = std::numeric_limits::max(); + /// Switch deciding whether surfaces without boundary + /// are cleared after a volume switch although they've not been reached yet + bool eraseUnboundVolChange{true}; + /// Delegate to decide whether free surfaces are appended to the navigation /// stream given the current volume and the track coordinates. If the /// delegate is set, it is called in each candidate resolution step diff --git a/Core/include/Acts/Surfaces/BoundaryTolerance.hpp b/Core/include/Acts/Surfaces/BoundaryTolerance.hpp index 3bdaed60176..c86abf1e090 100644 --- a/Core/include/Acts/Surfaces/BoundaryTolerance.hpp +++ b/Core/include/Acts/Surfaces/BoundaryTolerance.hpp @@ -224,7 +224,20 @@ class BoundaryTolerance { bool isTolerated(const Vector2& boundDelta, const SquareMatrix2& boundToCartesian) const; + /// Define the ostream operator to print the object + /// @param ostr: Reference to the ostream + /// @param tolerance: Reference to the tolerance to print + friend std::ostream& operator<<(std::ostream& ostr, + const BoundaryTolerance& tolerance) { + tolerance.print(ostr); + return ostr; + } + private: + /// Print method + /// @param ostr: Stream to which the object is printed + void print(std::ostream& ostr) const; + Variant m_variant; /// Check if the boundary check is of a specific type. diff --git a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp index b0f2e5ec93f..5efaa3c8c0e 100644 --- a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp +++ b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp @@ -30,6 +30,7 @@ #include "Acts/TrackFitting/detail/VoidFitterComponents.hpp" #include "Acts/Utilities/CalibrationContext.hpp" #include "Acts/Utilities/Delegate.hpp" +#include "Acts/Utilities/Helpers.hpp" #include "Acts/Utilities/Logger.hpp" #include "Acts/Utilities/Result.hpp" #include "Acts/Utilities/TrackHelpers.hpp" @@ -37,6 +38,7 @@ #include #include #include +#include #include #include @@ -185,6 +187,10 @@ struct Gx2FitterOptions { /// Max number of iterations during the fit (abort condition) std::size_t nUpdateMax = 5; + /// Number of iterations before missed measurements are tried to + /// be added back to the track + std::size_t nCoolDownMissed = 10; + /// Check for convergence (abort condition). Set to 0 to skip. double relChi2changeCutOff = 1e-7; }; @@ -231,9 +237,8 @@ struct Gx2FitterResult { /// Measurement surfaces without hits std::vector missedActiveSurfaces; - /// Measurement surfaces handled in both forward and - /// backward filtering - std::vector passedAgainSurfaces; + /// Hit measurement surfaces + std::vector hitActiveSurfaces; /// Count how many surfaces have been hit std::size_t surfaceCount = 0; @@ -756,6 +761,8 @@ class Gx2Fitter { /// Allows retrieving measurements for a surface const std::unordered_map* inputMeasurements{}; + /// Allows to skip measurements that cannot be considered as active + std::span vetoedActiveSurfaces{}; /// Whether to consider multiple scattering. bool multipleScattering = false; @@ -900,8 +907,10 @@ class Gx2Fitter { // Here we handle all measurements if (auto sourceLinkIt = inputMeasurements->find(surface); - sourceLinkIt != inputMeasurements->end()) { + sourceLinkIt != inputMeasurements->end() && + !Acts::rangeContainsValue(vetoedActiveSurfaces, surface)) { ACTS_DEBUG(" The surface contains a measurement."); + result.hitActiveSurfaces.push_back(surface); // Transport the covariance to the surface stepper.transportCovarianceToBound(state.stepping, *surface, @@ -1210,6 +1219,9 @@ class Gx2Fitter { ACTS_VERBOSE("Preparing " << std::distance(it, end) << " input measurements"); std::unordered_map inputMeasurements{}; + // Keep track of the hits that were not hit during the propagation + std::vector missedMeasurements{}; + std::size_t lastMissedIter{0u}; for (; it != end; ++it) { inputMeasurements.try_emplace(gx2fOptions.extensions.surfaceAccessor(*it), @@ -1266,6 +1278,12 @@ class Gx2Fitter { for (nUpdate = 0; nUpdate < gx2fOptions.nUpdateMax; nUpdate++) { ACTS_DEBUG("nUpdate = " << nUpdate + 1 << "/" << gx2fOptions.nUpdateMax); + if (!missedMeasurements.empty() && + lastMissedIter + gx2fOptions.nCoolDownMissed <= nUpdate) { + ACTS_DEBUG("Try to add back " << missedMeasurements.size() + << " measurements"); + missedMeasurements.clear(); + } // set up propagator and co PropagatorOptions propagatorOptions{gx2fOptions.propagatorPlainOptions}; @@ -1277,6 +1295,7 @@ class Gx2Fitter { auto& gx2fActor = propagatorOptions.actorList.template get(); gx2fActor.inputMeasurements = &inputMeasurements; + gx2fActor.vetoedActiveSurfaces = missedMeasurements; gx2fActor.multipleScattering = false; gx2fActor.extensions = gx2fOptions.extensions; gx2fActor.calibrationContext = &gx2fOptions.calibrationContext.get(); @@ -1296,6 +1315,7 @@ class Gx2Fitter { auto& r = propagatorState.template get>(); r.fittedStates = &trajectoryTempBackend; + r.hitActiveSurfaces.reserve(inputMeasurements.size()); // Clear the track container. It could be more performant to update the // existing states, but this needs some more thinking. @@ -1316,11 +1336,31 @@ class Gx2Fitter { // TODO Improve Propagator + Actor [allocate before loop], rewrite // makeMeasurements auto& propRes = *result; + GX2FResult gx2fResult = std::move(propRes.template get()); auto track = trackContainerTemp.makeTrack(); tipIndex = gx2fResult.lastMeasurementIndex; + if (gx2fResult.hitActiveSurfaces.size() + missedMeasurements.size() != + inputMeasurements.size()) { + ACTS_DEBUG("Propagation only hit " + << gx2fResult.hitActiveSurfaces.size() << "/" + << (inputMeasurements.size() - missedMeasurements.size()) + << " surfaces. Remove non hit surfaces from the list"); + for (const auto& [surface, _] : inputMeasurements) { + if (!Acts::rangeContainsValue(gx2fResult.hitActiveSurfaces, + surface)) { + missedMeasurements.push_back(surface); + ACTS_VERBOSE(" -- Remove hit associated to " + << surface->geometryId() << "."); + } + } + } + ACTS_VERBOSE("Found " << gx2fResult.hitActiveSurfaces.size() + << " measurements, " + << ", " << gx2fResult.measurementHoles << " holes"); + // It could happen, that no measurements were found. Then the track would // be empty and the following operations would be invalid. Usually, this // only happens during the first iteration, due to bad initial parameters. @@ -1372,6 +1412,9 @@ class Gx2Fitter { ACTS_INFO("Not enough measurements. Require " << extendedSystem.findRequiredNdf() + 1 << ", but only " << extendedSystem.ndf() << " could be used."); + for (const Surface* surface : gx2fResult.hitActiveSurfaces) { + ACTS_INFO(" --- " << surface->geometryId()); + } return Experimental::GlobalChiSquareFitterError::NotEnoughMeasurements; } @@ -1543,7 +1586,8 @@ class Gx2Fitter { ACTS_VERBOSE("aMatrix:\n" << extendedSystem.aMatrix() << "\n" - << "bVector:\n" + << "determinant: " << extendedSystem.aMatrix().determinant() + << ", bVector:\n" << extendedSystem.bVector() << "\n" << "deltaParamsExtended:\n" << deltaParamsExtended << "\n" diff --git a/Core/src/Navigation/NavigationStream.cpp b/Core/src/Navigation/NavigationStream.cpp index 20bfbe6d14e..eeda4cf66a6 100644 --- a/Core/src/Navigation/NavigationStream.cpp +++ b/Core/src/Navigation/NavigationStream.cpp @@ -84,8 +84,9 @@ bool NavigationStream::initialize(const GeometryContext& gctx, } // Append the multi intersection candidates - m_candidates.insert(m_candidates.end(), additionalCandidates.begin(), - additionalCandidates.end()); + m_candidates.insert(m_candidates.end(), + std::make_move_iterator(additionalCandidates.begin()), + std::make_move_iterator(additionalCandidates.end())); // Sort the candidates by path length std::ranges::sort(m_candidates, NavigationTarget::pathLengthOrder); diff --git a/Core/src/Navigation/Navigator.cpp b/Core/src/Navigation/Navigator.cpp index bc9de8a6785..d701c130689 100644 --- a/Core/src/Navigation/Navigator.cpp +++ b/Core/src/Navigation/Navigator.cpp @@ -313,6 +313,8 @@ void Navigator::handleSurfaceReached(State& state, const Vector3& position, // handling portals in gen3 configuration if (m_geometryVersion == GeometryVersion::Gen3) { + state.navCandidate().targetReached(); + if (state.navCandidate().isPortalTarget() && &state.navCandidate().surface() == &surface) { ACTS_VERBOSE(volInfo(state) << "Handling portal status."); @@ -581,12 +583,6 @@ void Navigator::resolveCandidates(State& state, const Vector3& position, "Navigator: No navigation policy found for current volume."); } - auto policyState = state.policyStateManager.currentState(); - state.currentVolume->initializeNavigationCandidates( - state.options.geoContext, args, policyState, appendOnly, logger()); - - ACTS_VERBOSE(volInfo(state) << "Found " << state.stream.candidates().size() - << " navigation candidates."); for (const Surface* surface : state.options.externalSurfaces) { const GeometryIdentifier geoId = surface->geometryId(); // Don't add any surface which is not in the same volume (volume bits) @@ -598,6 +594,28 @@ void Navigator::resolveCandidates(State& state, const Vector3& position, << " surface " << geoId); appendOnly.addSurfaceCandidate(*surface, BoundaryTolerance::Infinite()); } + + if (!state.options.eraseUnboundVolChange) { + // unbound targets are copied over to the next surface if they've not been + // yet reached by the propagator + for (const NavigationTarget& target : state.navCandidates) { + if (target.isSurfaceTarget() && target.boundaryTolerance().isInfinite() && + !target.isReached()) { + appendOnly.addSurfaceCandidate(target.surface(), + BoundaryTolerance::Infinite()); + ACTS_VERBOSE(volInfo(state) << "Target " << target + << " not been reached yet. Try now"); + } + } + } + + auto policyState = state.policyStateManager.currentState(); + state.currentVolume->initializeNavigationCandidates( + state.options.geoContext, args, policyState, appendOnly, logger()); + + ACTS_VERBOSE(volInfo(state) << "Found " << state.stream.candidates().size() + << " navigation candidates."); + bool pruneFreeCand{false}; if (!state.freeCandidates.empty()) { for (const auto& [surface, wasReached] : state.freeCandidates) { @@ -615,8 +633,11 @@ void Navigator::resolveCandidates(State& state, const Vector3& position, appendOnly.addSurfaceCandidate(*surface, BoundaryTolerance::Infinite()); pruneFreeCand = !state.options.freeSurfaceSelector.connected(); } - }; + } } + + state.navCandidates.clear(); + state.stream.initialize(state.options.geoContext, {position, direction}, BoundaryTolerance::None(), state.options.surfaceTolerance); @@ -626,8 +647,6 @@ void Navigator::resolveCandidates(State& state, const Vector3& position, << " navigation candidates after initialization.\n" << state.stream.candidates()); - state.navCandidates.clear(); - double farLimit = state.options.farLimit; // If the user has not provided the selection delegate, then // just apply a simple candidate pruning. Constrain the maximum diff --git a/Core/src/Propagator/NavigationTarget.cpp b/Core/src/Propagator/NavigationTarget.cpp index ce04ceb4eda..25bc7995cc6 100644 --- a/Core/src/Propagator/NavigationTarget.cpp +++ b/Core/src/Propagator/NavigationTarget.cpp @@ -32,6 +32,8 @@ void NavigationTarget::print(std::ostream& ostr) const { } }, m_target); - ostr << ", path length: " << pathLength(); + ostr << ", path length: " << pathLength() + << ", reached: " << (isReached() ? "yes" : "no") << ", " + << boundaryTolerance(); } } // namespace Acts diff --git a/Core/src/Surfaces/BoundaryTolerance.cpp b/Core/src/Surfaces/BoundaryTolerance.cpp index e467340ed3a..44308d35804 100644 --- a/Core/src/Surfaces/BoundaryTolerance.cpp +++ b/Core/src/Surfaces/BoundaryTolerance.cpp @@ -97,4 +97,20 @@ bool BoundaryTolerance::isTolerated( throw std::logic_error("Unsupported tolerance type"); } +void BoundaryTolerance::print(std::ostream& ostr) const { + if (isInfinite()) { + ostr << "no boundary limit applied"; + } else if (isNone()) { + ostr << "strict boundary limit"; + } else if (hasAbsoluteEuclidean()) { + ostr << "boundary limit dX < " << asAbsoluteEuclidean().tolerance; + } else if (hasChi2Bound()) { + ostr << "local chi2 " << asChi2Bound().maxChi2 << " with weights:\n" + << asChi2Bound().weightMatrix() << "\n"; + } else if (hasChi2Cartesian()) { + ostr << "local chi2 " << asChi2Cartesian().maxChi2 << " with weights:\n" + << asChi2Cartesian().weightMatrix() << "\n"; + } +} + } // namespace Acts