From 030b4c6d343622c76c385fe7a9303960699bb2a5 Mon Sep 17 00:00:00 2001 From: Tristan Wenzel Date: Fri, 31 Jul 2026 15:33:36 +0200 Subject: [PATCH] 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. --- Detectors/Base/CMakeLists.txt | 13 ++ .../include/DetectorsBase/GeometryManager.h | 17 +++ .../Base/include/DetectorsBase/MatLayerCyl.h | 7 +- .../include/DetectorsBase/MatLayerCylSet.h | 8 +- Detectors/Base/src/DetectorsBaseLinkDef.h | 1 + Detectors/Base/src/GeometryManager.cxx | 142 ++++++++++++++++++ Detectors/Base/src/MatLayerCyl.cxx | 17 ++- Detectors/Base/src/MatLayerCylSet.cxx | 101 ++++++++----- Detectors/Base/test/README.md | 11 +- Detectors/Base/test/buildMatBudLUT.C | 27 +++- Detectors/Base/test/compareMatBudLUT.C | 72 ++++++++- 11 files changed, 363 insertions(+), 53 deletions(-) diff --git a/Detectors/Base/CMakeLists.txt b/Detectors/Base/CMakeLists.txt index 74d2d02a9246c..9830d0c1175ce 100644 --- a/Detectors/Base/CMakeLists.txt +++ b/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 diff --git a/Detectors/Base/include/DetectorsBase/GeometryManager.h b/Detectors/Base/include/DetectorsBase/GeometryManager.h index 5e296a4045984..1d62eff126d09 100644 --- a/Detectors/Base/include/DetectorsBase/GeometryManager.h +++ b/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; diff --git a/Detectors/Base/include/DetectorsBase/MatLayerCyl.h b/Detectors/Base/include/DetectorsBase/MatLayerCyl.h index 26de279468477..04aefed06f010 100644 --- a/Detectors/Base/include/DetectorsBase/MatLayerCyl.h +++ b/Detectors/Base/include/DetectorsBase/MatLayerCyl.h @@ -20,6 +20,9 @@ #include #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 diff --git a/Detectors/Base/include/DetectorsBase/MatLayerCylSet.h b/Detectors/Base/include/DetectorsBase/MatLayerCylSet.h index 4408261dc740a..f2e8937ae8448 100644 --- a/Detectors/Base/include/DetectorsBase/MatLayerCylSet.h +++ b/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); diff --git a/Detectors/Base/src/DetectorsBaseLinkDef.h b/Detectors/Base/src/DetectorsBaseLinkDef.h index 8255c143ebb4a..2da1d5bdfad15 100644 --- a/Detectors/Base/src/DetectorsBaseLinkDef.h +++ b/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 + ; diff --git a/Detectors/Base/src/GeometryManager.cxx b/Detectors/Base/src/GeometryManager.cxx index 8aaf902aa95c5..225f21d8239a1 100644 --- a/Detectors/Base/src/GeometryManager.cxx +++ b/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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#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; + + 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(lvol->GetMaterialPtr()), budStep); + vecgeom::VNavigator const* navigator = lvol->GetNavigator(); + double step = static_cast(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 diff --git a/Detectors/Base/src/MatLayerCyl.cxx b/Detectors/Base/src/MatLayerCyl.cxx index dee179a84ca06..c35293d07c10e 100644 --- a/Detectors/Base/src/MatLayerCyl.cxx +++ b/Detectors/Base/src/MatLayerCyl.cxx @@ -109,7 +109,7 @@ 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); @@ -117,13 +117,13 @@ void MatLayerCyl::populateFromTGeo(int ntrPerCell) 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 diff --git a/Detectors/Base/src/MatLayerCylSet.cxx b/Detectors/Base/src/MatLayerCylSet.cxx index b85e3042cc379..b65a7472fbcf6 100644 --- a/Detectors/Base/src/MatLayerCylSet.cxx +++ b/Detectors/Base/src/MatLayerCylSet.cxx @@ -94,12 +94,16 @@ int MatLayerCylSet::getNThreadsFromEnv() } //________________________________________________________________________________ -void MatLayerCylSet::populateFromTGeo(int ntrPerCell, int nThreads) +void MatLayerCylSet::populateFromTGeo(int ntrPerCell, int nThreads, MatbudGeomBackend backend) { ///< populate layers, using ntrPerCell test tracks per cell. ///< nThreads < 0 takes the number of threads from the NTHREADS_MATBUD environment variable. assert(mConstructionMask == InProgress); + if (backend == MatbudGeomBackend::VECGEOM && !GeometryManager::isVecGeomAvailable()) { + LOG(fatal) << "MatbudGeomBackend::VECGEOM requested but O2 was built without VecGeom support (TGeo2VecGeom not found at configure time)"; + } + int nlr = getNLayers(); if (!nlr) { LOG(error) << "The LUT is not yet initialized"; @@ -124,13 +128,22 @@ void MatLayerCylSet::populateFromTGeo(int ntrPerCell, int nThreads) LOG(info) << "Populating with " << ntrPerCell << " trials Lr " << i; get()->mLayers[i].print(); } + const auto tSetupStart = Clock::now(); +#ifdef O2_WITH_VECGEOM + if (backend == MatbudGeomBackend::VECGEOM) { + // Trigger the lazy VecGeom world build/BVH init here so it counts as "setup" below, + // not as fill time for whichever cell happens first. + GeometryManager::vecGeomMaterialBudget(0.f, 0.f, 0.f, 0.f, 0.f, 1.f); + } +#endif const auto tFillStart = Clock::now(); for (int i = 0; i < nlr; i++) { - get()->mLayers[i].populateFromTGeo(ntrPerCell); + get()->mLayers[i].populateFromTGeo(ntrPerCell, backend); } const auto tFillEnd = Clock::now(); finalizeStructures(); - LOG(info) << "LUT fill: 1 thread, cells " << seconds(tFillStart, tFillEnd) << " s"; + LOG(info) << "LUT fill: 1 thread, setup " << seconds(tSetupStart, tFillStart) + << " s, cells " << seconds(tFillStart, tFillEnd) << " s"; return; } @@ -148,37 +161,57 @@ void MatLayerCylSet::populateFromTGeo(int ntrPerCell, int nThreads) const auto tSetupStart = Clock::now(); - // TGeo has to be told that several threads will navigate it, and each thread needs its own - // navigator. SetMaxThreads() is one-way -- TGeoManager has no API to return to - // single-threaded mode -- so we do not pretend to restore it; that is harmless because - // meanMaterialBudget() decides whether to lock from its own argument, not from this global. - // The navigators we book are ours, though, so those we do give back. - gGeoManager->SetMaxThreads(nThreads); - - tbb::enumerable_thread_specific threadNavigators( - []() { return gGeoManager->AddNavigator(); }); - - const auto tFillStart = Clock::now(); - { - tbb::global_control threadControl(tbb::global_control::max_allowed_parallelism, nThreads); - tbb::parallel_for(tbb::blocked_range(0, totalCells), - [this, ntrPerCell, &layerOffsets, &threadNavigators](const tbb::blocked_range& range) { - TGeoNavigator* nav = threadNavigators.local(); - for (size_t idx = range.begin(); idx != range.end(); ++idx) { - auto it = std::upper_bound(layerOffsets.begin(), layerOffsets.end(), idx); - const int layerIdx = int(std::distance(layerOffsets.begin(), it)) - 1; - const size_t cellInLayer = idx - layerOffsets[layerIdx]; - auto& layer = this->get()->mLayers[layerIdx]; - const int nphi = layer.getNPhiBins(); - layer.populateFromTGeo(int(cellInLayer % nphi), int(cellInLayer / nphi), ntrPerCell, nav); - } - }); - } - - const auto tFillEnd = Clock::now(); - - for (TGeoNavigator* nav : threadNavigators) { - gGeoManager->RemoveNavigator(nav); + auto fillRange = [this, ntrPerCell, backend, &layerOffsets](const tbb::blocked_range& range, TGeoNavigator* nav) { + for (size_t idx = range.begin(); idx != range.end(); ++idx) { + auto it = std::upper_bound(layerOffsets.begin(), layerOffsets.end(), idx); + const int layerIdx = int(std::distance(layerOffsets.begin(), it)) - 1; + const size_t cellInLayer = idx - layerOffsets[layerIdx]; + auto& layer = this->get()->mLayers[layerIdx]; + const int nphi = layer.getNPhiBins(); + layer.populateFromTGeo(int(cellInLayer % nphi), int(cellInLayer / nphi), ntrPerCell, nav, backend); + } + }; + + Clock::time_point tFillStart, tFillEnd; + if (backend == MatbudGeomBackend::ROOT) { + // TGeo has to be told that several threads will navigate it, and each thread needs its own + // navigator. SetMaxThreads() is one-way -- TGeoManager has no API to return to + // single-threaded mode -- so we do not pretend to restore it; that is harmless because + // meanMaterialBudget() decides whether to lock from its own argument, not from this global. + // The navigators we book are ours, though, so those we do give back. + gGeoManager->SetMaxThreads(nThreads); + + tbb::enumerable_thread_specific threadNavigators( + []() { return gGeoManager->AddNavigator(); }); + + tFillStart = Clock::now(); + { + tbb::global_control threadControl(tbb::global_control::max_allowed_parallelism, nThreads); + tbb::parallel_for(tbb::blocked_range(0, totalCells), + [&fillRange, &threadNavigators](const tbb::blocked_range& range) { + fillRange(range, threadNavigators.local()); + }); + } + tFillEnd = Clock::now(); + + for (TGeoNavigator* nav : threadNavigators) { + gGeoManager->RemoveNavigator(nav); + } + } else { + // VecGeom navigation needs no per-thread navigator bookkeeping. Trigger the lazy world + // build/BVH init before tFillStart so it counts as "setup", not fill time. +#ifdef O2_WITH_VECGEOM + GeometryManager::vecGeomMaterialBudget(0.f, 0.f, 0.f, 0.f, 0.f, 1.f); +#endif + tFillStart = Clock::now(); + { + tbb::global_control threadControl(tbb::global_control::max_allowed_parallelism, nThreads); + tbb::parallel_for(tbb::blocked_range(0, totalCells), + [&fillRange](const tbb::blocked_range& range) { + fillRange(range, nullptr); + }); + } + tFillEnd = Clock::now(); } finalizeStructures(); diff --git a/Detectors/Base/test/README.md b/Detectors/Base/test/README.md index 97e8c7f569fd1..ca7fea02949ba 100644 --- a/Detectors/Base/test/README.md +++ b/Detectors/Base/test/README.md @@ -14,7 +14,7 @@ root -b -q O2/Detectors/Base/test/buildMatBudLUT.C+ The generation is quite time consuming (may take ~30 min). It can be filled in parallel, one `TGeoNavigator` per thread, by passing a thread count as the -5th argument of `buildMatBudLUT` or by setting the environment variable: +7th argument of `buildMatBudLUT` or by setting the environment variable: ``` export NTHREADS_MATBUD=16 ``` @@ -23,6 +23,15 @@ ROOT >= v6-36-10-alice3, which removes the per-query thread-id lookup and the fa between the per-thread scratch buffers of TGeo shapes; with older ROOT the parallel path is still correct, just slower. +An alternative VecGeom geometry backend can be selected via the 8th argument (`"ROOT"` +or `"VECGEOM"`), e.g. +``` +root -b -q 'O2/Detectors/Base/test/buildMatBudLUT.C(60, -1, "matbud.root", "o2sim", "", 16, "VECGEOM")' +``` +This requires O2 to have been built against the optional `TGeo2VecGeom` package +(`o2::base::GeometryManager::isVecGeomAvailable()`); it is otherwise a build-time no-op that +does not affect the default ROOT/TGeo path in any way. + The optimized LUT will be stored in the matbud.root file. Load it as: diff --git a/Detectors/Base/test/buildMatBudLUT.C b/Detectors/Base/test/buildMatBudLUT.C index 44019198685f8..2b371b90effa3 100644 --- a/Detectors/Base/test/buildMatBudLUT.C +++ b/Detectors/Base/test/buildMatBudLUT.C @@ -23,14 +23,18 @@ #include #endif +using MatbudGeomBackend = o2::base::MatbudGeomBackend; + o2::base::MatLayerCylSet mbLUT; bool testMBLUT(const std::string& lutFile = "matbud.root"); +MatbudGeomBackend parseBackend(const std::string& s); /// Build the material budget LUT. nThreads < 0 takes the thread count from NTHREADS_MATBUD. +/// geomBackend is "ROOT" (default) or "VECGEOM" (requires O2 built against TGeo2VecGeom). bool buildMatBudLUT(int nTst = 60, int maxLr = -1, const std::string& outFile = "matbud.root", const std::string& geomNamePrefix = "o2sim", const std::string& opts = "", - int nThreads = -1); + int nThreads = -1, const std::string& geomBackend = "ROOT"); struct LrData { float rMin = 0.f; @@ -46,8 +50,9 @@ std::vector lrData; void configLayers(); bool buildMatBudLUT(int nTst, int maxLr, const std::string& outFile, const std::string& geomNamePrefix, - const std::string& opts, int nThreads) + const std::string& opts, int nThreads, const std::string& geomBackend) { + MatbudGeomBackend backend = parseBackend(geomBackend); auto geomName = o2::base::NameConf::getGeomFileName(geomNamePrefix); if (gSystem->AccessPathName(geomName.c_str())) { // if needed, create geometry std::cout << geomName << " does not exist. Will create it on the fly\n"; @@ -71,7 +76,7 @@ bool buildMatBudLUT(int nTst, int maxLr, const std::string& outFile, const std:: } TStopwatch sw; - mbLUT.populateFromTGeo(nTst, nThreads); + mbLUT.populateFromTGeo(nTst, nThreads, backend); mbLUT.optimizePhiSlices(); // move to populateFromTGeo mbLUT.flatten(); // move to populateFromTGeo @@ -401,3 +406,19 @@ void configLayers() lrData.emplace_back(LrData(lrData.back().rMax, lrData.back().rMax + drStep, zSpanH, zBin, rphiBin)); } while (lrData.back().rMax < 500); } + +//_______________________________________________________________________ +MatbudGeomBackend parseBackend(const std::string& s) +{ + if (s == "ROOT") { + return MatbudGeomBackend::ROOT; + } + if (s == "VECGEOM") { + if (!o2::base::GeometryManager::isVecGeomAvailable()) { + LOG(fatal) << "geomBackend=VECGEOM requested but O2 was built without VecGeom support (TGeo2VecGeom not found at configure time)"; + } + return MatbudGeomBackend::VECGEOM; + } + LOG(fatal) << "Unknown geomBackend '" << s << "', expected ROOT or VECGEOM"; + return MatbudGeomBackend::ROOT; +} diff --git a/Detectors/Base/test/compareMatBudLUT.C b/Detectors/Base/test/compareMatBudLUT.C index 7abaa9bd7182a..a58695d05d22c 100644 --- a/Detectors/Base/test/compareMatBudLUT.C +++ b/Detectors/Base/test/compareMatBudLUT.C @@ -12,21 +12,47 @@ /// \file compareMatBudLUT.C /// \brief Compare two material budget LUTs cell by cell /// -/// Used to check that filling the LUT in parallel gives the same result as filling it serially: +/// Used to check that filling the LUT in parallel gives the same result as filling it serially, +/// or that the VecGeom and ROOT geometry backends agree within a given tolerance: /// /// root -b -q 'compareMatBudLUT.C("matbud_serial.root","matbud_parallel.root")' +/// root -b -q 'compareMatBudLUT.C("matbud_ROOT.root","matbud_VECGEOM.root", 0.01, 20, "sweep.csv")' #if !defined(__CLING__) || defined(__ROOTCLING__) #include "DetectorsBase/MatLayerCylSet.h" #include "GPUCommonLogger.h" +#include #include #include +#include +#include #endif +namespace +{ +struct CellDiff { + int layer, iz, ip; + float rhoA, rhoB, x2x0A, x2x0B; + double rRho, rX; + double score() const { return std::max(rRho, rX); } +}; + +struct LayerStat { + float rMin = 0.f, rMax = 0.f; + size_t nBad = 0; + double maxRelRho = 0., maxRelX2X0 = 0.; +}; +} // namespace + /// Returns true if the two LUTs agree everywhere within tol (relative). +/// nWorst: number of worst-offending cells to print, ranked by max(relRho, relX2X0) over the +/// whole comparison, not just the first ones found in scan order. +/// csvSummary: if non-empty, append one summary row to this CSV file (header written once). bool compareMatBudLUT(const std::string& fileA = "matbud_serial.root", const std::string& fileB = "matbud_parallel.root", - float tol = 0.f) + float tol = 0.f, + int nWorst = 10, + const std::string& csvSummary = "") { auto* lutA = o2::base::MatLayerCylSet::loadFromFile(fileA); auto* lutB = o2::base::MatLayerCylSet::loadFromFile(fileB); @@ -46,6 +72,8 @@ bool compareMatBudLUT(const std::string& fileA = "matbud_serial.root", size_t nCells = 0, nBad = 0; double maxRelRho = 0., maxRelX2X0 = 0.; + std::vector layerStats(lutA->getNLayers()); + std::vector allCells; for (int il = 0; il < lutA->getNLayers(); il++) { const auto& la = lutA->getLayer(il); @@ -56,6 +84,10 @@ bool compareMatBudLUT(const std::string& fileA = "matbud_serial.root", << lb.getNZBins() << "x" << lb.getNPhiBins(); return false; } + auto& ls = layerStats[il]; + ls.rMin = la.getRMin(); + ls.rMax = la.getRMax(); + for (int iz = 0; iz < la.getNZBins(); iz++) { for (int ip = 0; ip < la.getNPhiBins(); ip++) { const auto& ca = la.getCellPhiBin(ip, iz); @@ -70,20 +102,48 @@ bool compareMatBudLUT(const std::string& fileA = "matbud_serial.root", const double rX = rel(ca.meanX2X0, cb.meanX2X0); maxRelRho = std::max(maxRelRho, rRho); maxRelX2X0 = std::max(maxRelX2X0, rX); + ls.maxRelRho = std::max(ls.maxRelRho, rRho); + ls.maxRelX2X0 = std::max(ls.maxRelX2X0, rX); if (rRho > tol || rX > tol) { - if (nBad < 10) { - printf("Lr %3d iz %4d ip %4d : rho %.9g vs %.9g (rel %.3g) | x2x0 %.9g vs %.9g (rel %.3g)\n", - il, iz, ip, ca.meanRho, cb.meanRho, rRho, ca.meanX2X0, cb.meanX2X0, rX); - } + ls.nBad++; nBad++; } + allCells.push_back({il, iz, ip, ca.meanRho, cb.meanRho, ca.meanX2X0, cb.meanX2X0, rRho, rX}); } } } + const int nw = std::min(nWorst, (int)allCells.size()); + std::partial_sort(allCells.begin(), allCells.begin() + nw, allCells.end(), + [](const CellDiff& a, const CellDiff& b) { return a.score() > b.score(); }); + printf("--- %d worst cells (by max relative deviation) ---\n", nw); + for (int i = 0; i < nw; i++) { + const auto& c = allCells[i]; + printf("Lr %3d iz %4d ip %4d : rho %.9g vs %.9g (rel %.3g) | x2x0 %.9g vs %.9g (rel %.3g)\n", + c.layer, c.iz, c.ip, c.rhoA, c.rhoB, c.rRho, c.x2x0A, c.x2x0B, c.rX); + } + + printf("--- per-layer summary (%d layers) ---\n", lutA->getNLayers()); + for (int il = 0; il < lutA->getNLayers(); il++) { + const auto& ls = layerStats[il]; + printf("Lr %3d %8.3fgetNLayers()); printf("Max relative difference: meanRho %.3g, meanX2X0 %.3g (tolerance %.3g)\n", maxRelRho, maxRelX2X0, tol); + + if (!csvSummary.empty()) { + const bool writeHeader = gSystem->AccessPathName(csvSummary.c_str()); // true if it does NOT exist + std::ofstream csv(csvSummary, std::ios::app); + if (writeHeader) { + csv << "fileA,fileB,nLayers,nCells,nBad,maxRelRho,maxRelX2X0,tol\n"; + } + csv << fileA << "," << fileB << "," << lutA->getNLayers() << "," << nCells << "," << nBad << "," + << maxRelRho << "," << maxRelX2X0 << "," << tol << "\n"; + } + if (nBad) { LOG(error) << nBad << " cells differ beyond tolerance"; return false;