diff --git a/Detectors/gconfig/src/VecGeomG4Navigator.cxx b/Detectors/gconfig/src/VecGeomG4Navigator.cxx index 8644f9a29e625..eb7af0440b8f6 100644 --- a/Detectors/gconfig/src/VecGeomG4Navigator.cxx +++ b/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; @@ -314,6 +315,9 @@ G4VPhysicalVolume* VecGeomG4Navigator::LocateGlobalPointAndSetup(const G4ThreeVe // Into the daughter the step hit, then down inside it. mReloScratch = mCurState; mCurState = mNextState; + // mNextState is a copy of the step state and can still carry the block of the volume the previous + // crossing left, which was meant for the first step after that exit only. + mCurState.SetLastExited(mEmptyState.GetLastExitedState()); auto const* daughter = mCurState.Top(); mCurState.Pop(); vecgeom::Transformation3D m; @@ -416,7 +420,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 +432,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; } diff --git a/Detectors/gconfig/src/VecGeomG4Navigator.h b/Detectors/gconfig/src/VecGeomG4Navigator.h index 12e281c040016..300c318090f10 100644 --- a/Detectors/gconfig/src/VecGeomG4Navigator.h +++ b/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}; diff --git a/Detectors/gconfig/src/VecGeomG4NavigatorBase.cxx b/Detectors/gconfig/src/VecGeomG4NavigatorBase.cxx index 7c6ebb22e98f1..cffe9374bb182 100644 --- a/Detectors/gconfig/src/VecGeomG4NavigatorBase.cxx +++ b/Detectors/gconfig/src/VecGeomG4NavigatorBase.cxx @@ -13,10 +13,17 @@ #include "G4VPhysicalVolume.hh" +#include +#include +#include +#include +#include #include #include +#include + 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::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 diff --git a/Detectors/gconfig/src/VecGeomG4NavigatorBase.h b/Detectors/gconfig/src/VecGeomG4NavigatorBase.h index 6f86368a2b221..f3ebeb1dcc9d7 100644 --- a/Detectors/gconfig/src/VecGeomG4NavigatorBase.h +++ b/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) diff --git a/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.cxx b/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.cxx index 4057c9c65dd3e..b03fa14b0b959 100644 --- a/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.cxx +++ b/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.cxx @@ -24,6 +24,7 @@ #include +#include #include 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) diff --git a/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.h b/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.h index 4ccf469966304..2a5133d6d4514 100644 --- a/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.h +++ b/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;