o2-sim: VecGeom navigation mode for Geant4 (#15865)

* Set up the VecGeom BVH navigator and safety estimators for every volume

This lets GeometryManager set up VecGeom for Geant4 navigation as well as for the material budget.

- buildVecGeomGeometry() converts once and takes the assembly flattening as a parameter; the
  material budget keeps flattening.
- With a VecGeom that has BVHNavigatorV, volumes with more than two daughters get the BVH navigator
  and level locator, and every volume gets an explicit safety estimator.
- Without it the VecGeom v2 setup is unchanged.

https://gitlab.cern.ch/VecGeom/VecGeom/-/merge_requests/1547

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

* Add a VecGeom navigation mode to the Geant4 engine

This adds G4.navmode=kVecGeom, in which Geant4 answers every navigation query from VecGeom.

- The geometry, materials and touchables stay the ones g4root builds from TGeo.
- VecGeomG4Map pairs each VecGeom placement with the chain of g4root volumes it stands for, so
  volume ids, copy numbers and CurrentVolOffID are unchanged.
- VecGeomG4Navigator relocates at the boundary locate, with the volume just left blocked for the
  next step and the point pushed across the face by a small depth, as G4VecGeomNav's
  TG4VecGeomNavigator does.
- G4.vecgeomCheckRays, vecgeomCheckLocation and vecgeomCheckVolumes compare VecGeom with TGeo
  before transport.
- It is built only when TGeo2VecGeom and a VecGeom with BVHNavigatorV are found.

https://gitlab.cern.ch/VecGeom/g4vecgeomnav/-/merge_requests/25

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

* Add a lighter VecGeom navigator that adopts the propagated state

This adds G4.vecgeomNavigator=kPropagated, a second VecGeom navigator kept for comparison.

- It adopts the state VecGeom propagated during the step instead of relocating at the boundary.
- It returns zero safety for the step after a crossing and nudges stuck steps forward.
- The default navigator becomes G4.vecgeomNavigator=kRelocating.

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

* Block the exited volume only as Geant4 does

This fixes tracks entering the volume they just left without seeing its boundary.

- Both VecGeom navigators kept the volume just left blocked for every following step, whatever the direction.
- A track turned back into it at the exit point, as by multiple scattering, crossed it unseen in the mother.
- In the hollow RB24 copper tubes this let shower electrons pass the wall and gave about 10 % more FDD-A hits.
- The volume is now blocked only in the first step after the exit, and only while the direction points away
  from it, as in G4NormalNavigation. This follows the same fix in G4VecGeomNav.

https://gitlab.cern.ch/swenzel/g4vecgeomnav/-/tree/fix/block-and-reflection

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

* Please consider the following formatting changes

---------

Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
Co-authored-by: ALICE Action Bot <alibuild@cern.ch>

Author: 3 people

Common/SimConfig/include/SimConfig/G4Params.h

@@ -36,8 +36,15 @@ enum class EG4Physics {
 // enumerating possible geometry navigation modes
 // (understanding that geometry description is always done with TGeo)
 enum class EG4Nav {
-  kTGeo = 0, /* navigate with TGeo */
-  kG4 = 1    /* navigate with G4 native geometry */
+  kTGeo = 0,   /* navigate with TGeo */
+  kG4 = 1,     /* navigate with G4 native geometry */
+  kVecGeom = 2 /* navigate with VecGeom, on the G4 geometry built from TGeo */
+};
+
+// the Geant4 navigator used with navmode kVecGeom
+enum class EVecGeomNav {
+  kRelocating = 0, /* relocates at the boundary locate, blocking the volume just left (default) */
+  kPropagated = 1  /* adopts the state VecGeom propagated during the step; less work per crossing */
 };
 
 // parameters to influence the G4 engine
@@ -49,6 +56,22 @@ struct G4Params : public o2::conf::ConfigurableParamHelper<G4Params> {
 
   EG4Nav navmode = EG4Nav::kTGeo; // geometry navigation mode (default TGeo)
 
+  // Settings for navmode == kVecGeom; ignored otherwise.
+  // which of the two VecGeom navigators
+  EVecGeomNav vecgeomNavigator = EVecGeomNav::kRelocating;
+  double vecgeomPushDepth = 1.e-9;      // cm; how far past a face, measured across it, a boundary
+                                        // point is pushed before it is located
+  bool vecgeomZeroSafety = false;       // answer zero to every safety query; conservative, but it
+                                        // 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 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
+                                        // the volumes they disagree on
+  std::string vecgeomCheckVolumes = ""; // comma-separated volumes to cross-check by sampling inside
+                                        // their placements
+
   std::string fluenceWeightFile = ""; // file containing the scoring weights (pdg, ekin, weight)
   std::string const& getPhysicsConfigString() const;
 

Detectors/Base/include/DetectorsBase/GeometryManager.h

@@ -135,6 +135,12 @@ class GeometryManager : public TObject
   /// Mean material budget between two points, using the VecGeom backend. On first call,
   /// lazily converts the currently loaded TGeo geometry to VecGeom (once per process).
   static o2::base::MatBudget vecGeomMaterialBudget(float x0, float y0, float z0, float x1, float y1, float z1);
+  /// Converts the currently loaded TGeo geometry to VecGeom and assigns a navigator, a level
+  /// 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);
 #else
   static constexpr bool isVecGeomAvailable() { return false; }
 #endif

Detectors/Base/src/GeometryManager.cxx

@@ -49,6 +49,14 @@
 #include <VecGeom/navigation/VNavigator.h>
 #include <VecGeom/volumes/LogicalVolume.h>
 #include <mutex>
+// The BVH navigator of the VNavigator family, which Geant4 navigation needs on every volume.
+#if __has_include(<VecGeom/navigation/BVHNavigatorV.h>)
+#define O2_VECGEOM_HAS_BVH_VNAVIGATOR
+#include <VecGeom/navigation/BVHNavigatorV.h>
+#include <VecGeom/navigation/BVHLevelLocator.h>
+#include <VecGeom/navigation/BVHSafetyEstimator.h>
+#include <VecGeom/navigation/SimpleSafetyEstimator.h>
+#endif
 #endif
 
 using namespace o2::detectors;
@@ -574,15 +582,22 @@ bool usesBvhAcceleration(vecgeom::LogicalVolume const* vol)
 /// process, the first time the VecGeom backend is requested. Not part of loadGeometry(),
 /// which every job calls regardless of whether it ever uses the VecGeom backend.
 void ensureVecGeomWorldBuilt()
+{
+  GeometryManager::buildVecGeomGeometry(true);
+}
+} // namespace
+
+void GeometryManager::buildVecGeomGeometry(bool flattenAssemblies)
 {
   static std::once_flag onceFlag;
-  std::call_once(onceFlag, []() {
+  std::call_once(onceFlag, [flattenAssemblies]() {
     if (!gGeoManager) {
       LOG(fatal) << "Cannot build VecGeom geometry: no TGeo geometry loaded (call GeometryManager::loadGeometry() first)";
     }
     // Translate geometry and material pointers, then build acceleration structures.
     tgeo2vecgeom::RootGeoManager::Instance().SetMaterialConversionHook([](TGeoMaterial const* m) { return (void*)m; });
-    tgeo2vecgeom::RootGeoManager::Instance().SetFlattenAssemblies(true);
+    LOG(info) << "VecGeom conversion: flattenAssemblies=" << flattenAssemblies;
+    tgeo2vecgeom::RootGeoManager::Instance().SetFlattenAssemblies(flattenAssemblies);
     tgeo2vecgeom::RootGeoManager::Instance().LoadRootGeometry();
 
     // Acceleration structures must be built before the navigators/locators reference them.
@@ -593,15 +608,28 @@ void ensureVecGeomWorldBuilt()
     // Builds a BVH per logical volume.
     vecgeom::BVHManager::Init();
 
-    // For each logical volume, set both a navigator (used for ComputeStep) and a matched
-    // level locator (used for point relocation after a boundary crossing via GlobalLocator).
+    // For each logical volume, set a navigator (used for ComputeStep), a matched level locator
+    // (used for point relocation after a boundary crossing via GlobalLocator) and, where the
+    // VNavigator family is complete, the safety estimator LogicalVolume::GetSafetyEstimator()
+    // hands out, which is separate from the one a navigator uses internally.
     for (auto& lvol : vecgeom::GeoManager::Instance().GetLogicalVolumesMap()) {
       auto* vol = lvol.second;
       if (!usesBvhAcceleration(vol)) {
         vol->SetNavigator(vecgeom::NewSimpleNavigator<>::Instance());
+#ifdef O2_VECGEOM_HAS_BVH_VNAVIGATOR
+        vol->SetLevelLocator(vol->ContainsAssembly() ? vecgeom::SimpleAssemblyLevelLocator::GetInstance()
+                                                     : vecgeom::SimpleLevelLocator::GetInstance());
+        vol->SetSafetyEstimator(vecgeom::SimpleSafetyEstimator::Instance());
+#else
         vol->SetLevelLocator(vecgeom::SimpleLevelLocator::GetInstance());
+#endif
       } else {
-#if VECGEOM_VERSION >= 0x020000
+#if defined(O2_VECGEOM_HAS_BVH_VNAVIGATOR)
+        vol->SetNavigator(vecgeom::BVHNavigatorV<>::Instance());
+        vol->SetLevelLocator(vol->ContainsAssembly() ? vecgeom::BVHAssemblyAwareLevelLocator::GetInstance()
+                                                     : vecgeom::BVHLevelLocator::GetInstance());
+        vol->SetSafetyEstimator(vecgeom::BVHSafetyEstimator::Instance());
+#elif VECGEOM_VERSION >= 0x020000
         // VecGeom 2 turned BVHNavigator into a plain class with static entry points instead of a
         // VNavigator singleton, so there is nothing to attach: vecGeomMaterialBudget() calls it
         // directly.
@@ -627,7 +655,6 @@ void ensureVecGeomWorldBuilt()
     }
   });
 }
-} // namespace
 
 //_____________________________________________________________________________________
 o2::base::MatBudget GeometryManager::vecGeomMaterialBudget(float x0, float y0, float z0, float x1, float y1, float z1)

Detectors/gconfig/CMakeLists.txt

@@ -14,11 +14,40 @@ o2_add_library(G3Setup
                PUBLIC_LINK_LIBRARIES MC::Geant3 FairRoot::Base O2::SimulationDataFormat O2::Generators O2::SimSetup
 )
 
+# Optional VecGeom navigation for Geant4 (G4.navmode=kVecGeom). It needs TGeo2VecGeom and a VecGeom
+# with the BVH navigator of the VNavigator family (BVHNavigatorV). Linked PRIVATE: VecGeom types
+# never appear in G4Setup's public headers.
+find_package(TGeo2VecGeom CONFIG QUIET)
+set(G4SETUP_WITH_VECGEOM OFF)
+if(TGeo2VecGeom_FOUND)
+  find_path(O2_VECGEOM_BVHNAVIGATORV_INCLUDE VecGeom/navigation/BVHNavigatorV.h
+            HINTS ${VecGeom_INCLUDE_DIRS} ${VecGeom_DIR}/../../../include $ENV{VECGEOM_ROOT}/include)
+  if(O2_VECGEOM_BVHNAVIGATORV_INCLUDE)
+    set(G4SETUP_WITH_VECGEOM ON)
+  endif()
+endif()
+
+set(G4SETUP_SOURCES src/G4Config.cxx src/G4RunConfiguration.cxx src/G4LocalFieldConstruction.cxx
+                    src/VecGeomNavigation.cxx)
+if(G4SETUP_WITH_VECGEOM)
+  list(APPEND G4SETUP_SOURCES src/VecGeomG4Map.cxx src/VecGeomChecks.cxx src/VecGeomG4NavigatorBase.cxx
+                              src/VecGeomG4Navigator.cxx
+                              src/VecGeomG4PropagatingNavigator.cxx)
+endif()
+
 o2_add_library(G4Setup
-               SOURCES src/G4Config.cxx src/G4RunConfiguration.cxx src/G4LocalFieldConstruction.cxx
-               PUBLIC_LINK_LIBRARIES MC::Geant4VMC MC::Geant4 FairRoot::Base O2::SimulationDataFormat O2::Generators O2::SimSetup O2::FastSim
+               TARGETVARNAME targetG4Setup
+               SOURCES ${G4SETUP_SOURCES}
+               PUBLIC_LINK_LIBRARIES MC::Geant4VMC MC::Geant4 FairRoot::Base O2::SimulationDataFormat O2::Generators O2::SimSetup O2::FastSim O2::DetectorsBase
 )
 
+if(G4SETUP_WITH_VECGEOM)
+  target_compile_definitions(${targetG4Setup} PRIVATE O2_WITH_VECGEOM)
+  target_link_libraries(${targetG4Setup} PRIVATE TGeo2VecGeom::TGeo2VecGeom)
+else()
+  message(STATUS "G4.navmode=kVecGeom not built: it needs TGeo2VecGeom and a VecGeom with BVHNavigatorV")
+endif()
+
 o2_add_library(FLUKASetup
                SOURCES src/FlukaConfig.cxx
                PUBLIC_LINK_LIBRARIES FairRoot::Base O2::SimulationDataFormat O2::Generators O2::SimSetup

Detectors/gconfig/g4Config.C

@@ -66,6 +66,7 @@ R__LOAD_LIBRARY(libgeant4vmc)
 #include "G4VScoringMesh.hh"
 #include <unistd.h>
 #include "SimSetup/G4RunConfiguration.h"
+#include "SimSetup/VecGeomNavigation.h"
 #endif
 #include "commonConfig.C"
 
@@ -115,6 +116,10 @@ void Config()
     geomNavStr = "geomRoot";
   } else if (g4Params.navmode == o2::conf::EG4Nav::kG4) {
     geomNavStr = "geomVMC+RootToGeant4";
+  } else if (g4Params.navmode == o2::conf::EG4Nav::kVecGeom) {
+    // The geometry, its materials and the touchable stay the ones g4root builds from TGeo;
+    // only the navigator is swapped, once the engine below has built that hierarchy.
+    geomNavStr = "geomRoot";
   } else {
     LOG(fatal) << "Unsupported geometry navigation mode";
   }
@@ -137,6 +142,10 @@ void Config()
   TGeant4* geant4 = new TGeant4("TGeant4", "The Geant4 Monte Carlo", runConfiguration);
   std::cout << "Geant4 has been created." << std::endl;
 
+  if (g4Params.navmode == o2::conf::EG4Nav::kVecGeom) {
+    o2::simsetup::installVecGeomNavigator();
+  }
+
   // setup the stack
   stackSetup(geant4, FairRunSim::Instance());
 

Detectors/gconfig/include/SimSetup/VecGeomNavigation.h

@@ -0,0 +1,32 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+#ifndef O2_SIMSETUP_VECGEOMNAVIGATION_H_
+#define O2_SIMSETUP_VECGEOMNAVIGATION_H_
+
+namespace o2::simsetup
+{
+
+/// Whether this build of O2 has the VecGeom navigation backend, i.e. whether TGeo2VecGeom and a
+/// VecGeom with BVHNavigatorV were found when O2 was configured.
+bool isVecGeomNavigationAvailable();
+
+/// Replaces Geant4's tracking navigator by one that answers every navigation query from
+/// VecGeom. The Geant4 geometry, its materials and the touchable the scoring code reads stay
+/// the ones g4root built from TGeo, so only navigation changes.
+///
+/// Call after the TGeant4 engine has been constructed: the Geant4 hierarchy this needs to map
+/// onto is built while TG4RunManager configures itself. Aborts if the backend is missing.
+void installVecGeomNavigator();
+
+} // namespace o2::simsetup
+
+#endif