Honour the safety bound and reuse the last safety within its sphere

This makes the stand-alone ComputeSafety of both VecGeom navigators cheaper, as Geant4's own navigator does.

- The caller's proposedMaxLength now limits the BVH search: the value is exact below it and a valid lower bound above it, as in G4VoxelSafety.
- A query at distance d from the last safety origin returns s0 - d when that covers the bound, without a computation. Every locate drops this cache.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Author: sawenzel

Detectors/gconfig/src/VecGeomG4Navigator.cxx

@@ -123,7 +123,7 @@ G4double VecGeomG4Navigator::ComputeStep(const G4ThreeVector& globalPoint, const
   // On the point a boundary locate left the track on, the safety is zero; a point seen before
   // reuses its safety. Otherwise the navigator computes it with the step.
   bool calcSafety = !mZeroSafety && !(mLocatedOnBoundary && samePoint(globalPoint, mLastLocatedPoint));
-  if (calcSafety && samePoint(globalPoint, mSafetyOrig)) {
+  if (calcSafety && samePoint(globalPoint, mSafetyOrig) && mLastSafety < mLastSafetyLimit) {
     calcSafety = false;
     newSafety = mLastSafety;
   }
@@ -154,6 +154,7 @@ G4double VecGeomG4Navigator::ComputeStep(const G4ThreeVector& globalPoint, const
     newSafety = safety * kVGToG4;
     mSafetyOrig = globalPoint;
     mLastSafety = newSafety;
+    mLastSafetyLimit = kInfinity;
   }
   const bool boundaryLimited = vgStep < limit;
   G4double step = std::max(vgStep, 0.) * kVGToG4;
@@ -416,7 +417,8 @@ void VecGeomG4Navigator::LocateGlobalPointWithinVolume(const G4ThreeVector& posi
   clearLastExited();
 }
 
-G4double VecGeomG4Navigator::ComputeSafety(const G4ThreeVector& globalPoint, const G4double, const G4bool)
+G4double VecGeomG4Navigator::ComputeSafety(const G4ThreeVector& globalPoint, const G4double proposedMaxLength,
+                                           const G4bool)
 {
   if (mZeroSafety) {
     return 0.;
@@ -427,24 +429,26 @@ G4double VecGeomG4Navigator::ComputeSafety(const G4ThreeVector& globalPoint, con
   if ((mWouldEnter || mWouldExit) && samePoint(globalPoint, mNextPoint)) {
     return 0.;
   }
-  if (samePoint(globalPoint, mSafetyOrig)) {
+  if (samePoint(globalPoint, mSafetyOrig) && (mLastSafety < mLastSafetyLimit || mLastSafety >= proposedMaxLength)) {
     return mLastSafety;
   }
-  auto const* top = topOf(mCurState);
-  if (top == nullptr) {
+  if (topOf(mCurState) == nullptr) {
     return 0.;
   }
-  auto const* estimator = top->GetLogicalVolume()->GetSafetyEstimator();
-  if (estimator == nullptr) {
-    return 0.;
+  // Every point closer to the last safety origin than its safety is in the same volume, so s0 - d is a
+  // valid safety there; it is used when it covers the caller's bound. Every locate drops the cache.
+  const double rest = mLastSafety - std::sqrt(globalPoint.diff2(mSafetyOrig));
+  if (rest >= proposedMaxLength) {
+    return rest;
   }
-  double safety = estimator->ComputeSafety(toVG(globalPoint), mCurState);
+  double safety = boundedSafety(mCurState, globalPoint, proposedMaxLength);
   if (safety < 0.) {
     ++mNegativeSafetyCount;
     safety = 0.;
   }
   mSafetyOrig = globalPoint;
-  mLastSafety = safety * kVGToG4;
+  mLastSafety = safety;
+  mLastSafetyLimit = proposedMaxLength;
   return mLastSafety;
 }
 

Detectors/gconfig/src/VecGeomG4Navigator.h

@@ -28,7 +28,8 @@ namespace o2::simsetup
 ///   The point is first pushed across the face by a small depth, and afterwards leaves every volume
 ///   it is flush with and heading out of. As in Geant4, the volume left is blocked only in the first
 ///   ComputeStep after the exit, and only while the direction points away from it.
-/// - Safety is zero only at the boundary point itself.
+/// - Safety is zero only at the boundary point itself. A bounded query is answered from the last safety
+///   computed in the same volume when that covers the bound.
 ///
 /// Unlike TG4VecGeomNavigator, the VecGeom geometry is converted from TGeo, not from Geant4, so one
 /// VecGeom placement can stand for a chain of g4root volumes (VecGeomG4Map).
@@ -110,6 +111,7 @@ class VecGeomG4Navigator : public VecGeomG4NavigatorBase
   bool mLocatedOnBoundary = false;
   G4ThreeVector mSafetyOrig{-1e8, -1e8, -1e8}; ///< the last point a safety was computed for
   double mLastSafety = 0.;                     ///< mm
+  double mLastSafetyLimit = kInfinity;         ///< mm; mLastSafety is exact below it, a lower bound above
   bool mNormalEnter = false;
   bool mNormalValid = false;
   G4ThreeVector mNormalPoint{-1e8, -1e8, -1e8};

Detectors/gconfig/src/VecGeomG4NavigatorBase.cxx

@@ -13,10 +13,17 @@
 
 #include "G4VPhysicalVolume.hh"
 
+#include <VecGeom/management/BVHManager.h>
+#include <VecGeom/navigation/BVHNavigator.h>
+#include <VecGeom/navigation/BVHSafetyEstimator.h>
+#include <VecGeom/navigation/VSafetyEstimator.h>
+#include <VecGeom/volumes/LogicalVolume.h>
 #include <VecGeom/volumes/PlacedVolume.h>
 
 #include <fairlogger/Logger.h>
 
+#include <algorithm>
+
 namespace
 {
 /// Longest run of Geant4 levels one flattened VecGeom placement can stand for.
@@ -28,6 +35,30 @@ constexpr int kMaxDepth = 64;
 namespace o2::simsetup
 {
 
+double VecGeomG4NavigatorBase::boundedSafety(vecgeom::NavigationState const& state, const G4ThreeVector& point,
+                                             double limit)
+{
+  auto const* pvol = state.Top();
+  auto const* lvol = pvol->GetLogicalVolume();
+  auto const* estimator = lvol->GetSafetyEstimator();
+  if (estimator == nullptr) {
+    return 0.;
+  }
+  if (limit < kInfinity && estimator == vecgeom::BVHSafetyEstimator::Instance() && lvol->GetDaughters().size() > 0) {
+    // The BVH estimator's own computation, with the search limited.
+    vecgeom::Transformation3D m;
+    state.TopMatrix(m);
+    const V3 local = m.Transform(toVG(point));
+    double safety = pvol->SafetyToOut(local);
+    if (safety > 0.) {
+      safety = vecgeom::BVHNavigator::ComputeBVHSafety<vecgeom::BVHSafetyEstimator>(
+        *vecgeom::BVHManager::GetBVH(lvol), local, safety, std::min(safety, limit * kG4ToVG));
+    }
+    return safety * kVGToG4;
+  }
+  return estimator->ComputeSafety(toVG(point), state) * kVGToG4;
+}
+
 G4VPhysicalVolume* VecGeomG4NavigatorBase::historyFromState(vecgeom::NavigationState const& state)
 {
   // Collect the Geant4 volumes the state stands for, then keep the levels the history already has

Detectors/gconfig/src/VecGeomG4NavigatorBase.h

@@ -50,6 +50,10 @@ class VecGeomG4NavigatorBase : public G4Navigator
   /// False if a level matches no VecGeom placement; \a state is then empty.
   bool stateFromHistory(vecgeom::NavigationState& state) const;
 
+  /// The safety (mm) at \a point in the top volume of \a state, exact below \a limit (mm) and a valid
+  /// lower bound above it: daughters beyond the limit are bounded by their boxes, as in G4VoxelSafety.
+  static double boundedSafety(vecgeom::NavigationState const& state, const G4ThreeVector& point, double limit);
+
   /// Geant4's rule for the volume a track just left (G4NormalNavigation, G4VoxelNavigation): it is
   /// blocked only while the direction points away from it, along its outward normal \a n (global).
   [[gnu::always_inline]] static bool directionLeaves(const V3& n, const V3& dir)

Detectors/gconfig/src/VecGeomG4PropagatingNavigator.cxx

@@ -24,6 +24,7 @@
 
 #include <fairlogger/Logger.h>
 
+#include <algorithm>
 #include <cmath>
 
 namespace
@@ -133,6 +134,7 @@ G4VPhysicalVolume* VecGeomG4PropagatingNavigator::ResetHierarchyAndLocate(const
   mWouldExit = false;
   mOnBoundary = false;
   mHaveNextState = false;
+  mLastSafety = -1.;
   fHistory = *history.GetHistory();
   if (!stateFromHistory(mCurState) && fHistory.GetVolume(0) != nullptr) {
     LOG(fatal) << "Geant4 handed back a touchable that matches no VecGeom path";
@@ -153,6 +155,7 @@ G4VPhysicalVolume* VecGeomG4PropagatingNavigator::LocateGlobalPointAndSetup(cons
   mPrevState = mCurState;
   mLocatedPoint = point;
   mExitBlockPending = false;
+  mLastSafety = -1.;
 
   if (!mForceReInit && relativeSearch && onBoundary && mHaveNextState) {
     // The state on the far side of the boundary was already worked out, and relocated, by the step
@@ -190,17 +193,24 @@ void VecGeomG4PropagatingNavigator::LocateGlobalPointWithinVolume(const G4ThreeV
   fExitedMother = false;
 }
 
-G4double VecGeomG4PropagatingNavigator::ComputeSafety(const G4ThreeVector& globalPoint, const G4double, const G4bool)
+G4double VecGeomG4PropagatingNavigator::ComputeSafety(const G4ThreeVector& globalPoint, const G4double proposedMaxLength,
+                                                     const G4bool)
 {
   if (mZeroSafety || mOnBoundary || mCrossed || fEnteredDaughter || fExitedMother || mWouldEnter || mWouldExit) {
     return 0.;
   }
-  auto const* top = mCurState.Top();
-  if (top == nullptr) {
+  if (mCurState.Top() == nullptr) {
     return 0.;
   }
-  const double safety = top->GetLogicalVolume()->GetSafetyEstimator()->ComputeSafety(toVG(globalPoint), mCurState);
-  return (safety > 0.) ? safety * kVGToG4 : 0.;
+  // Every point closer to the last safety origin than its safety is in the same volume, so s0 - d is a
+  // valid safety there; it is used when it covers the caller's bound. Every locate drops the cache.
+  const double rest = mLastSafety - std::sqrt(globalPoint.diff2(mSafetyOrig));
+  if (rest >= proposedMaxLength) {
+    return rest;
+  }
+  mSafetyOrig = globalPoint;
+  mLastSafety = std::max(boundedSafety(mCurState, globalPoint, proposedMaxLength), 0.);
+  return mLastSafety;
 }
 
 G4ThreeVector VecGeomG4PropagatingNavigator::GetGlobalExitNormal(const G4ThreeVector& point, G4bool* valid)

Detectors/gconfig/src/VecGeomG4PropagatingNavigator.h

@@ -65,6 +65,8 @@ class VecGeomG4PropagatingNavigator : public VecGeomG4NavigatorBase
   bool mCrossed = false;                         ///< the last locate acted on a boundary crossing
   bool mExitBlockPending = false;                ///< the next step is the first after leaving mPrevState's volume
   G4ThreeVector mLocatedPoint{-1e8, -1e8, -1e8}; ///< where the last locate put the track
+  G4ThreeVector mSafetyOrig{-1e8, -1e8, -1e8};   ///< the last point a safety was computed for
+  double mLastSafety = -1.;                      ///< mm; negative when there is none
 
   int mZeroSteps = 0;
   long mNudgedSteps = 0;