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;