o2-sim: VecGeom navigation mode for Geant4

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