Support to use VecGeom geometry navigation in the material scan

This commit provides support to use VecGeom as an alternative
geometry backend to perform the material budget LUT creation.

It is the first time that ALICE interfaces full VecGeom navigation.
This should be a useful milestone towards adoption at detector
simulation level.

First scans indicate performance advantages of 2x over TGeo but
a precise performance investigation will be done later. This commit
will be the foundation for this and/or serve for debugging issues
or driving further development in VecGeom itself.

Note that the original method name `populateFromTGeo` remains untouched
for API stability, despite the fact that the method can now internally
use VecGeom.

Author: Tristan Wenzel

Detectors/Base/CMakeLists.txt

@@ -10,7 +10,13 @@
 # or submit itself to any jurisdiction.
 #add_compile_options(-O0 -g -fPIC)
 
+# Optional VecGeom backend for material-budget LUT filling, off by default when TGeo2VecGeom
+# isn't installed (guarded by O2_WITH_VECGEOM). Linked PRIVATE: VecGeom types never appear in
+# DetectorsBase's public headers, so consumers need neither VecGeom headers nor VecGeom itself.
+find_package(TGeo2VecGeom CONFIG QUIET)
+
 o2_add_library(DetectorsBase
+               TARGETVARNAME targetDetectorsBase
                SOURCES src/Detector.cxx
                        src/GeometryManager.cxx
                        src/MaterialManager.cxx
@@ -51,6 +57,13 @@ o2_add_library(DetectorsBase
                                      ROOT::Gdml
                              )
 
+if(TGeo2VecGeom_FOUND)
+  target_compile_definitions(${targetDetectorsBase} PRIVATE O2_WITH_VECGEOM)
+  target_link_libraries(${targetDetectorsBase} PRIVATE TGeo2VecGeom::TGeo2VecGeom)
+else()
+  message(STATUS "TGeo2VecGeom not found: DetectorsBase built without the optional VecGeom material-budget backend")
+endif()
+
 o2_target_root_dictionary(DetectorsBase
                           HEADERS include/DetectorsBase/Detector.h
                                   include/DetectorsBase/GeometryManager.h

Detectors/Base/include/DetectorsBase/GeometryManager.h

@@ -40,6 +40,12 @@ class AlignParam;
 
 namespace base
 {
+/// Backend used to compute material budget for LUT filling: ROOT/TGeo (default, always
+/// available) or VecGeom (requires O2 to be built against the optional TGeo2VecGeom
+/// package; see GeometryManager::isVecGeomAvailable()).
+enum class MatbudGeomBackend : int { ROOT = 0,
+                                     VECGEOM = 1 };
+
 /// Class for interfacing to the geometry; it also builds and manages the look-up tables for fast
 /// access to geometry and alignment information for sensitive alignable volumes:
 /// 1) the look-up table mapping unique volume ids to TGeoPNEntries. This allows to access
@@ -121,6 +127,17 @@ class GeometryManager : public TObject
     return meanMaterialBudgetExt(start.X(), start.Y(), start.Z(), end.X(), end.Y(), end.Z());
   }
 
+  /// Whether this build of O2 was configured with the optional VecGeom material-budget
+  /// backend (i.e. TGeo2VecGeom was found at CMake configure time).
+#ifdef O2_WITH_VECGEOM
+  static constexpr bool isVecGeomAvailable() { return true; }
+  /// 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);
+#else
+  static constexpr bool isVecGeomAvailable() { return false; }
+#endif
+
  private:
   /// Default constructor
   GeometryManager() = default;

Detectors/Base/include/DetectorsBase/MatLayerCyl.h

@@ -20,6 +20,9 @@
 #include <cstring>
 #endif
 #include "GPUCommonDef.h"
+#ifndef GPUCA_ALIGPUCODE
+#include "DetectorsBase/GeometryManager.h" // for MatbudGeomBackend
+#endif
 #include "FlatObject.h"
 #include "GPUCommonRtypes.h"
 #include "GPUCommonMath.h"
@@ -67,8 +70,8 @@ class MatLayerCyl : public o2::gpu::FlatObject
 
   void initSegmentation(float rMin, float rMax, float zHalfSpan, int nz, int nphi);
   void initSegmentation(float rMin, float rMax, float zHalfSpan, float dzMin, float drphiMin);
-  void populateFromTGeo(int ntrPerCell = 10);
-  void populateFromTGeo(int ip, int iz, int ntrPerCell, TGeoNavigator* nav = nullptr);
+  void populateFromTGeo(int ntrPerCell = 10, MatbudGeomBackend backend = MatbudGeomBackend::ROOT);
+  void populateFromTGeo(int ip, int iz, int ntrPerCell, TGeoNavigator* nav = nullptr, MatbudGeomBackend backend = MatbudGeomBackend::ROOT);
   void print(bool data = false) const;
 #endif // !GPUCA_ALIGPUCODE
 

Detectors/Base/include/DetectorsBase/MatLayerCylSet.h

@@ -73,9 +73,11 @@ class MatLayerCylSet : public o2::gpu::FlatObject
 #ifndef GPUCA_ALIGPUCODE // this part is unvisible on GPU version
   void print(bool data = false) const;
   void addLayer(float rmin, float rmax, float zmax, float dz, float drphi);
-  /// Populate the LUT from TGeo. nThreads > 1 fills the cells in parallel (one TGeoNavigator
-  /// per thread); nThreads < 0 takes the count from the NTHREADS_MATBUD environment variable.
-  void populateFromTGeo(int ntrPerCel = 10, int nThreads = -1);
+  /// Populate the LUT from TGeo or VecGeom. nThreads > 1 fills cells in parallel (one
+  /// TGeoNavigator per thread for ROOT; VecGeom navigation is thread-safe on its own);
+  /// nThreads < 0 takes the count from NTHREADS_MATBUD. VECGEOM requires O2 built against
+  /// TGeo2VecGeom, see GeometryManager::isVecGeomAvailable().
+  void populateFromTGeo(int ntrPerCel = 10, int nThreads = -1, MatbudGeomBackend backend = MatbudGeomBackend::ROOT);
   static int getNThreadsFromEnv();
   void optimizePhiSlices(float maxRelDiff = 0.05);
 

Detectors/Base/src/DetectorsBaseLinkDef.h

@@ -24,6 +24,7 @@
 
 #pragma link C++ class o2::base::GeometryManager + ;
 #pragma link C++ class o2::base::GeometryManager::MatBudgetExt + ;
+#pragma link C++ enum o2::base::MatbudGeomBackend;
 #pragma link C++ class o2::base::MaterialManager + ;
 #pragma link C++ class o2::MaterialManagerParam + ;
 #pragma link C++ class o2::GeometryManagerParam + ;

Detectors/Base/src/GeometryManager.cxx

@@ -30,6 +30,22 @@
 #include "CommonUtils/NameConf.h"
 #include "DetectorsBase/Aligner.h"
 
+#ifdef O2_WITH_VECGEOM
+#include "TGeo2VecGeom/RootGeoManager.h"
+#include <VecGeom/management/GeoManager.h>
+#include <VecGeom/management/ABBoxManager.h>
+#include <VecGeom/management/BVHManager.h>
+#include <VecGeom/navigation/GlobalLocator.h>
+#include <VecGeom/navigation/NavigationState.h>
+#include <VecGeom/navigation/NewSimpleNavigator.h>
+#include <VecGeom/navigation/BVHNavigator.h>
+#include <VecGeom/navigation/SimpleLevelLocator.h>
+#include <VecGeom/navigation/BVHLevelLocator.h>
+#include <VecGeom/navigation/VNavigator.h>
+#include <VecGeom/volumes/LogicalVolume.h>
+#include <mutex>
+#endif
+
 using namespace o2::detectors;
 using namespace o2::base;
 
@@ -536,3 +552,129 @@ void GeometryManager::loadGeometry(std::string_view simPrefix, bool applyMisalig
     applyMisalignent(applyMisalignment);
   }
 }
+
+#ifdef O2_WITH_VECGEOM
+
+namespace
+{
+/// Converts the currently loaded TGeo geometry to VecGeom and sets up navigators, once per
+/// 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()
+{
+  static std::once_flag onceFlag;
+  std::call_once(onceFlag, []() {
+    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);
+    tgeo2vecgeom::RootGeoManager::Instance().LoadRootGeometry();
+
+    // Acceleration structures must be built before the navigators/locators reference them.
+    vecgeom::ABBoxManager::Instance().InitABBoxesForCompleteGeometry();
+    // Builds a BVH per logical volume from the ABBoxes computed above.
+    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);
+    // volumes with very few daughters are cheaper to brute-force than to accelerate.
+    for (auto& lvol : vecgeom::GeoManager::Instance().GetLogicalVolumesMap()) {
+      auto* vol = lvol.second;
+      if (vol->GetDaughtersp()->size() <= 2) {
+        vol->SetNavigator(vecgeom::NewSimpleNavigator<>::Instance());
+        vol->SetLevelLocator(vecgeom::SimpleLevelLocator::GetInstance());
+      } else {
+        vol->SetNavigator(vecgeom::BVHNavigator<>::Instance());
+        vol->SetLevelLocator(vecgeom::BVHLevelLocator::GetInstance());
+      }
+    }
+  });
+}
+} // namespace
+
+//_____________________________________________________________________________________
+o2::base::MatBudget GeometryManager::vecGeomMaterialBudget(float x0, float y0, float z0, float x1, float y1, float z1)
+{
+  // Mean material budget between "0" and "1" via VecGeom's BVH-accelerated ray/boundary
+  // intersection, instead of TGeo.
+  ensureVecGeomWorldBuilt();
+
+  using Vector3D = vecgeom::Vector3D<vecgeom::Precision>;
+
+  double length, start[3] = {x0, y0, z0};
+  double dir[3] = {x1 - x0, y1 - y0, z1 - z0};
+  if ((length = dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2]) < TGeoShape::Tolerance() * TGeoShape::Tolerance()) {
+    return o2::base::MatBudget(); // return empty struct
+  }
+  length = std::sqrt(length);
+  double invlen = 1. / length;
+  for (int i = 3; i--;) {
+    dir[i] *= invlen;
+  }
+
+  thread_local static vecgeom::NavigationState* newnavstate = vecgeom::NavigationState::MakeInstance(vecgeom::GeoManager::Instance().getMaxDepth());
+  thread_local static vecgeom::NavigationState* currnavstate = vecgeom::NavigationState::MakeInstance(vecgeom::GeoManager::Instance().getMaxDepth());
+  thread_local static vecgeom::NavigationState* startCache = vecgeom::NavigationState::MakeInstance(vecgeom::GeoManager::Instance().getMaxDepth());
+  thread_local static bool startCacheValid = false;
+
+  Vector3D currPoint(x0, y0, z0);
+  Vector3D dirr(dir[0], dir[1], dir[2]);
+  constexpr double kPush = 1.E-6; // mimick the nudging of TGeo's FindNextBoundaryAndStep
+  auto world = vecgeom::GeoManager::Instance().GetWorld();
+  o2::base::MatBudget budTot, budStep;
+  budStep.length = length;
+
+  // Locate the starting volume, reusing the path from the previous call when still valid.
+  if (startCacheValid && !startCache->IsOutside()) {
+    startCache->CopyTo(currnavstate);
+    vecgeom::Transformation3D m;
+    currnavstate->TopMatrix(m);
+    vecgeom::GlobalLocator::RelocatePointFromPath(m.Transform(currPoint), *currnavstate);
+  } else {
+    currnavstate->Clear();
+    vecgeom::GlobalLocator::LocateGlobalPoint(world, currPoint, *currnavstate, true);
+  }
+  if (currnavstate->IsOutside() || currnavstate->Top() == nullptr) {
+    LOG(error) << "start point out of geometry: " << x0 << ':' << y0 << ':' << z0;
+    startCacheValid = false;
+    return o2::base::MatBudget();
+  }
+  currnavstate->CopyTo(startCache);
+  startCacheValid = true;
+
+  double stepTot = 0.;
+  double remainingDist = length;
+  Int_t nzero = 0;
+  while (remainingDist > 1.E-10) {
+    auto* lvol = currnavstate->Top()->GetLogicalVolume();
+    accountMaterial(static_cast<TGeoMaterial*>(lvol->GetMaterialPtr()), budStep);
+    vecgeom::VNavigator const* navigator = lvol->GetNavigator();
+    double step = static_cast<double>(navigator->ComputeStepAndPropagatedState(currPoint, dirr, remainingDist, *currnavstate, *newnavstate));
+    if (step < 2.E-10) {
+      nzero++;
+    } else {
+      nzero = 0;
+    }
+    if (nzero > 3) {
+      // This means navigation has problems on one boundary
+      LOG(warning) << "Cannot cross boundary at (" << currPoint[0] << ',' << currPoint[1] << ',' << currPoint[2] << ')';
+      budTot.meanRho /= stepTot;
+      budTot.length = stepTot;
+      return o2::base::MatBudget(budTot);
+    }
+
+    remainingDist -= step;
+    stepTot += step;
+    budTot.meanRho += step * budStep.meanRho;
+    budTot.meanX2X0 += step / budStep.meanX2X0;
+    currPoint = currPoint + (step + kPush) * dirr;
+    std::swap(currnavstate, newnavstate);
+  }
+  budTot.meanRho /= stepTot;
+  budTot.length = stepTot;
+  return o2::base::MatBudget(budTot);
+}
+
+#endif // O2_WITH_VECGEOM

Detectors/Base/src/MatLayerCyl.cxx

@@ -109,21 +109,21 @@ void MatLayerCyl::initSegmentation(float rMin, float rMax, float zHalfSpan, int
 }
 
 //________________________________________________________________________________
-void MatLayerCyl::populateFromTGeo(int ntrPerCell)
+void MatLayerCyl::populateFromTGeo(int ntrPerCell, MatbudGeomBackend backend)
 {
   /// populate layer with info extracted from TGeometry, using ntrPerCell test tracks per cell
   assert(mConstructionMask != Constructed);
   mConstructionMask = InProgress;
   ntrPerCell = ntrPerCell > 1 ? ntrPerCell : 1;
   for (int iz = getNZBins(); iz--;) {
     for (int ip = getNPhiBins(); ip--;) {
-      populateFromTGeo(ip, iz, ntrPerCell);
+      populateFromTGeo(ip, iz, ntrPerCell, nullptr, backend);
     }
   }
 }
 
 //________________________________________________________________________________
-void MatLayerCyl::populateFromTGeo(int ip, int iz, int ntrPerCell, TGeoNavigator* nav)
+void MatLayerCyl::populateFromTGeo(int ip, int iz, int ntrPerCell, TGeoNavigator* nav, MatbudGeomBackend backend)
 {
   /// populate cell with info extracted from TGeometry, using ntrPerCell test tracks per cell
 
@@ -136,7 +136,16 @@ void MatLayerCyl::populateFromTGeo(int ip, int iz, int ntrPerCell, TGeoNavigator
     float dzt = zs > 0.f ? 0.25 * dz : -0.25 * dz; // to avoid 90 degree polar angle
     for (int isp = ntrPerCell; isp--;) {
       o2::math_utils::sincos(phmn + (isp + 0.5) * getDPhi() / ntrPerCell, sn, cs);
-      auto bud = o2::base::GeometryManager::meanMaterialBudget(rMin * cs, rMin * sn, zs - dzt, rMax * cs, rMax * sn, zs + dzt, nav);
+      o2::base::MatBudget bud;
+      if (backend == MatbudGeomBackend::ROOT) {
+        bud = o2::base::GeometryManager::meanMaterialBudget(rMin * cs, rMin * sn, zs - dzt, rMax * cs, rMax * sn, zs + dzt, nav);
+      } else {
+#ifdef O2_WITH_VECGEOM
+        bud = o2::base::GeometryManager::vecGeomMaterialBudget(rMin * cs, rMin * sn, zs - dzt, rMax * cs, rMax * sn, zs + dzt);
+#else
+        LOG(fatal) << "MatbudGeomBackend::VECGEOM requested but O2 was built without VecGeom support (TGeo2VecGeom not found at configure time)";
+#endif
+      }
       if (bud.length > 0.) {
         meanRho += bud.length * bud.meanRho;
         meanX2X0 += bud.meanX2X0; // we store actually not X2X0 but 1./X0