AliceO2Group · sawenzel · Jul 31, 2026
@@ -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

@@ -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;

@@ -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
 

@@ -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);
 

@@ -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 + ;

@@ -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
@@ -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