From 1354f4782537938bda60026875aceb34e4ce8caf Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Fri, 2 Oct 2026 11:23:55 +0200 Subject: [PATCH 1/3] Let VecGeom turn large Boolean solids into MultiUnions This lets VecGeom's BooleanFactory rewrite Boolean solids of at least 8 components into MultiUnions in the kVecGeom mode. - In the ALICE geometry 170 composite shapes are rewritten, among them the MFT disk supports and the ITS end wheels. - G4.vecgeomBooleanThreshold sets the threshold; 0 keeps the solids as converted. - Other users of buildVecGeomGeometry keep the solids as converted. Co-Authored-By: Claude Opus 5.5 --- Common/SimConfig/include/SimConfig/G4Params.h | 2 ++ .../Base/include/DetectorsBase/GeometryManager.h | 5 +++-- Detectors/Base/src/GeometryManager.cxx | 13 +++++++++++-- Detectors/gconfig/src/VecGeomNavigation.cxx | 2 +- 4 files changed, 17 insertions(+), 5 deletions(-) diff --git a/Common/SimConfig/include/SimConfig/G4Params.h b/Common/SimConfig/include/SimConfig/G4Params.h index 57ae4a28da91c..25e06801560c9 100644 --- a/Common/SimConfig/include/SimConfig/G4Params.h +++ b/Common/SimConfig/include/SimConfig/G4Params.h @@ -65,6 +65,8 @@ struct G4Params : public o2::conf::ConfigurableParamHelper { // shortens steps and so changes the random history bool vecgeomFlattenAssemblies = true; // dissolve TGeo assemblies into their content when converting // to VecGeom; the Geant4 touchable keeps the assembly levels + int vecgeomBooleanThreshold = 8; // Boolean solids with at least this many components become + // VecGeom MultiUnions; 0 keeps them as converted int vecgeomCheckRays = 0; // if > 0, step this many rays out of the interaction point with // TGeo and VecGeom and report the volumes they enter differently int vecgeomCheckLocation = 0; // if > 0, locate this many random points with both and report diff --git a/Detectors/Base/include/DetectorsBase/GeometryManager.h b/Detectors/Base/include/DetectorsBase/GeometryManager.h index d325a68929dd7..f918c278dcf33 100644 --- a/Detectors/Base/include/DetectorsBase/GeometryManager.h +++ b/Detectors/Base/include/DetectorsBase/GeometryManager.h @@ -139,8 +139,9 @@ class GeometryManager : public TObject /// locator and a safety estimator to every logical volume. Does the work once per process; later /// calls, whatever they ask for, return the geometry already built, so a caller that needs a /// particular assembly treatment must come first. \param flattenAssemblies dissolves TGeo - /// assemblies into their content. - static void buildVecGeomGeometry(bool flattenAssemblies); + /// assemblies into their content. \param booleanThreshold lets VecGeom turn Boolean solids with at + /// least this many components into MultiUnions; 0 keeps them as converted. + static void buildVecGeomGeometry(bool flattenAssemblies, int booleanThreshold = 0); #else static constexpr bool isVecGeomAvailable() { return false; } #endif diff --git a/Detectors/Base/src/GeometryManager.cxx b/Detectors/Base/src/GeometryManager.cxx index a3b61217aee01..f0ea0d41e7d04 100644 --- a/Detectors/Base/src/GeometryManager.cxx +++ b/Detectors/Base/src/GeometryManager.cxx @@ -587,10 +587,10 @@ void ensureVecGeomWorldBuilt() } } // namespace -void GeometryManager::buildVecGeomGeometry(bool flattenAssemblies) +void GeometryManager::buildVecGeomGeometry(bool flattenAssemblies, int booleanThreshold) { static std::once_flag onceFlag; - std::call_once(onceFlag, [flattenAssemblies]() { + std::call_once(onceFlag, [flattenAssemblies, booleanThreshold]() { if (!gGeoManager) { LOG(fatal) << "Cannot build VecGeom geometry: no TGeo geometry loaded (call GeometryManager::loadGeometry() first)"; } @@ -598,6 +598,15 @@ void GeometryManager::buildVecGeomGeometry(bool flattenAssemblies) tgeo2vecgeom::RootGeoManager::Instance().SetMaterialConversionHook([](TGeoMaterial const* m) { return (void*)m; }); LOG(info) << "VecGeom conversion: flattenAssemblies=" << flattenAssemblies; tgeo2vecgeom::RootGeoManager::Instance().SetFlattenAssemblies(flattenAssemblies); +#if __has_include() + // Applied when the conversion closes the geometry. + LOG(info) << "VecGeom conversion: booleanThreshold=" << booleanThreshold; + vecgeom::GeoManager::Instance().SetBooleanOptimizationThreshold(booleanThreshold > 0 ? booleanThreshold : 0); +#else + if (booleanThreshold > 0) { + LOG(warning) << "VecGeom without BooleanFactory: Boolean solids are kept as converted"; + } +#endif tgeo2vecgeom::RootGeoManager::Instance().LoadRootGeometry(); // Acceleration structures must be built before the navigators/locators reference them. diff --git a/Detectors/gconfig/src/VecGeomNavigation.cxx b/Detectors/gconfig/src/VecGeomNavigation.cxx index 68dbc96f883c3..fa3da5a6e340e 100644 --- a/Detectors/gconfig/src/VecGeomNavigation.cxx +++ b/Detectors/gconfig/src/VecGeomNavigation.cxx @@ -80,7 +80,7 @@ void installVecGeomNavigator() TStopwatch timer; timer.Start(); - o2::base::GeometryManager::buildVecGeomGeometry(g4Params.vecgeomFlattenAssemblies); + o2::base::GeometryManager::buildVecGeomGeometry(g4Params.vecgeomFlattenAssemblies, g4Params.vecgeomBooleanThreshold); timer.Stop(); LOG(info) << "VecGeom geometry built in " << timer.RealTime() << " s"; From bc3b37398b7cb9ecb42106c863bc8c876a0789b4 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Fri, 2 Oct 2026 11:24:35 +0200 Subject: [PATCH 2/3] 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 --- Detectors/gconfig/src/VecGeomG4Navigator.cxx | 24 ++++++++------ Detectors/gconfig/src/VecGeomG4Navigator.h | 4 ++- .../gconfig/src/VecGeomG4NavigatorBase.cxx | 31 +++++++++++++++++++ .../gconfig/src/VecGeomG4NavigatorBase.h | 4 +++ .../src/VecGeomG4PropagatingNavigator.cxx | 20 +++++++++--- .../src/VecGeomG4PropagatingNavigator.h | 2 ++ 6 files changed, 69 insertions(+), 16 deletions(-) diff --git a/Detectors/gconfig/src/VecGeomG4Navigator.cxx b/Detectors/gconfig/src/VecGeomG4Navigator.cxx index 8644f9a29e625..d23e341beeda3 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; @@ -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; } 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; From 2beea349e2a3df17d603d16b3a981cfc7e7aba6e Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Fri, 2 Oct 2026 11:24:43 +0200 Subject: [PATCH 3/3] Clear the previous exit block when entering a daughter This fixes a stale exit block in the relocating VecGeom navigator. - Entering a daughter adopted the state of the step, which could still carry the block of the volume the previous crossing left, so that volume stayed blocked one step longer than Geant4's rule allows. - Seven of 50000 geantino rays through the ITS outer-barrel services saw up to 9 X0 too much copper because of it; they now agree with TGeo. Co-Authored-By: Claude Opus 5.5 --- Detectors/gconfig/src/VecGeomG4Navigator.cxx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Detectors/gconfig/src/VecGeomG4Navigator.cxx b/Detectors/gconfig/src/VecGeomG4Navigator.cxx index d23e341beeda3..eb7af0440b8f6 100644 --- a/Detectors/gconfig/src/VecGeomG4Navigator.cxx +++ b/Detectors/gconfig/src/VecGeomG4Navigator.cxx @@ -315,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;