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/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; 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";