diff --git a/Detectors/Base/CMakeLists.txt b/Detectors/Base/CMakeLists.txt index 83a9193274e4f..74d2d02a9246c 100644 --- a/Detectors/Base/CMakeLists.txt +++ b/Detectors/Base/CMakeLists.txt @@ -92,6 +92,7 @@ if(BUILD_SIMULATION) endif() install(FILES test/buildMatBudLUT.C + test/compareMatBudLUT.C test/extractLUTLayers.C test/rescaleLUT.C DESTINATION share/macro/) @@ -100,6 +101,10 @@ o2_add_test_root_macro(test/buildMatBudLUT.C PUBLIC_LINK_LIBRARIES O2::DetectorsBase LABELS detectorsbase) +o2_add_test_root_macro(test/compareMatBudLUT.C + PUBLIC_LINK_LIBRARIES O2::DetectorsBase + LABELS detectorsbase) + o2_add_test_root_macro(test/extractLUTLayers.C PUBLIC_LINK_LIBRARIES O2::DetectorsBase LABELS detectorsbase) diff --git a/Detectors/Base/include/DetectorsBase/GeometryManager.h b/Detectors/Base/include/DetectorsBase/GeometryManager.h index ad1d77b10f49a..5e296a4045984 100644 --- a/Detectors/Base/include/DetectorsBase/GeometryManager.h +++ b/Detectors/Base/include/DetectorsBase/GeometryManager.h @@ -29,6 +29,7 @@ #include class TGeoHMatrix; // lines 11-11 class TGeoManager; // lines 9-9 +class TGeoNavigator; namespace o2 { @@ -96,14 +97,18 @@ class GeometryManager : public TObject ClassDefNV(MatBudgetExt, 1); }; - static o2::base::MatBudget meanMaterialBudget(float x0, float y0, float z0, float x1, float y1, float z1); - static o2::base::MatBudget meanMaterialBudget(const math_utils::Point3D& start, const math_utils::Point3D& end) + /// Mean material budget between two points. Pass a navigator owned by the calling thread to + /// run lock-free from several threads; with nav = nullptr the shared navigator is used under + /// a mutex, as before. + static o2::base::MatBudget meanMaterialBudget(float x0, float y0, float z0, float x1, float y1, float z1, + TGeoNavigator* nav = nullptr); + static o2::base::MatBudget meanMaterialBudget(const math_utils::Point3D& start, const math_utils::Point3D& end, TGeoNavigator* nav = nullptr) { - return meanMaterialBudget(start.X(), start.Y(), start.Z(), end.X(), end.Y(), end.Z()); + return meanMaterialBudget(start.X(), start.Y(), start.Z(), end.X(), end.Y(), end.Z(), nav); } - static o2::base::MatBudget meanMaterialBudget(const math_utils::Point3D& start, const math_utils::Point3D& end) + static o2::base::MatBudget meanMaterialBudget(const math_utils::Point3D& start, const math_utils::Point3D& end, TGeoNavigator* nav = nullptr) { - return meanMaterialBudget(start.X(), start.Y(), start.Z(), end.X(), end.Y(), end.Z()); + return meanMaterialBudget(start.X(), start.Y(), start.Z(), end.X(), end.Y(), end.Z(), nav); } static MatBudgetExt meanMaterialBudgetExt(float x0, float y0, float z0, float x1, float y1, float z1); diff --git a/Detectors/Base/include/DetectorsBase/MatLayerCyl.h b/Detectors/Base/include/DetectorsBase/MatLayerCyl.h index e63de51e0a6ca..d65ea86ad2a83 100644 --- a/Detectors/Base/include/DetectorsBase/MatLayerCyl.h +++ b/Detectors/Base/include/DetectorsBase/MatLayerCyl.h @@ -25,6 +25,8 @@ #include "GPUCommonMath.h" #include "DetectorsBase/MatCell.h" +class TGeoNavigator; + namespace o2 { namespace base @@ -66,7 +68,7 @@ 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); + void populateFromTGeo(int ip, int iz, int ntrPerCell, TGeoNavigator* nav = nullptr); 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 cba6e5cebcfc8..4408261dc740a 100644 --- a/Detectors/Base/include/DetectorsBase/MatLayerCylSet.h +++ b/Detectors/Base/include/DetectorsBase/MatLayerCylSet.h @@ -73,7 +73,10 @@ 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); - void populateFromTGeo(int ntrPerCel = 10); + /// 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); + static int getNThreadsFromEnv(); void optimizePhiSlices(float maxRelDiff = 0.05); void dumpToTree(const std::string& outName = "matbudTree.root") const; diff --git a/Detectors/Base/src/GeometryManager.cxx b/Detectors/Base/src/GeometryManager.cxx index a067767752a69..8aaf902aa95c5 100644 --- a/Detectors/Base/src/GeometryManager.cxx +++ b/Detectors/Base/src/GeometryManager.cxx @@ -16,6 +16,7 @@ #include // for TIter #include #include // for TGeoHMatrix +#include // for TGeoNavigator #include // for TGeoNode #include // for TGeoPhysicalNode, TGeoPNEntry #include @@ -398,7 +399,8 @@ GeometryManager::MatBudgetExt GeometryManager::meanMaterialBudgetExt(float x0, f } //_____________________________________________________________________________________ -o2::base::MatBudget GeometryManager::meanMaterialBudget(float x0, float y0, float z0, float x1, float y1, float z1) +o2::base::MatBudget GeometryManager::meanMaterialBudget(float x0, float y0, float z0, float x1, float y1, float z1, + TGeoNavigator* nav) { // // Calculate mean material budget and material properties between @@ -414,6 +416,8 @@ o2::base::MatBudget GeometryManager::meanMaterialBudget(float x0, float y0, floa // // Ported to O2: ruben.shahoyan@cern.ch // + // Multi-threaded execution: pass a navigator owned by the calling thread. + // double length, startD[3] = {x0, y0, z0}; double dir[3] = {x1 - x0, y1 - y0, z1 - z0}; @@ -425,9 +429,17 @@ o2::base::MatBudget GeometryManager::meanMaterialBudget(float x0, float y0, floa for (int i = 3; i--;) { dir[i] *= invlen; } - std::lock_guard guard(sTGMutex); + // A caller that passes its own navigator owns it exclusively, so no locking is needed. A caller + // that passes none shares gGeoManager's current navigator and must still serialize. Deciding + // this from the argument keeps the choice local: it does not depend on -- and cannot be broken + // by -- process-global state such as TGeoManager::GetMaxThreads(). + std::unique_lock guard(sTGMutex, std::defer_lock); + if (!nav) { + guard.lock(); + nav = gGeoManager->GetCurrentNavigator(); + } // Initialize start point and direction - TGeoNode* currentnode = gGeoManager->InitTrack(startD, dir); + TGeoNode* currentnode = nav->InitTrack(startD, dir); if (!currentnode) { LOG(error) << "start point out of geometry: " << x0 << ':' << y0 << ':' << z0; return o2::base::MatBudget(); // return empty struct @@ -439,11 +451,11 @@ o2::base::MatBudget GeometryManager::meanMaterialBudget(float x0, float y0, floa // Locate next boundary within length without computing safety. // Propagate either with length (if no boundary found) or just cross boundary - gGeoManager->FindNextBoundaryAndStep(length, kFALSE); + nav->FindNextBoundaryAndStep(length, kFALSE); Double_t stepTot = 0.0; // Step made - Double_t step = gGeoManager->GetStep(); + Double_t step = nav->GetStep(); // If no boundary within proposed length, return current step data - if (!gGeoManager->IsOnBoundary()) { + if (!nav->IsOnBoundary()) { budStep.meanX2X0 = budStep.length / budStep.meanX2X0; return o2::base::MatBudget(budStep); } @@ -458,7 +470,7 @@ o2::base::MatBudget GeometryManager::meanMaterialBudget(float x0, float y0, floa if (nzero > 3) { // This means navigation has problems on one boundary // Try to cross by making a small step - const double* curPos = gGeoManager->GetCurrentPoint(); + const double* curPos = nav->GetCurrentPoint(); LOG(warning) << "Cannot cross boundary at (" << curPos[0] << ',' << curPos[1] << ',' << curPos[2] << ')'; budTotal.meanRho /= stepTot; budTotal.length = stepTot; @@ -472,14 +484,14 @@ o2::base::MatBudget GeometryManager::meanMaterialBudget(float x0, float y0, floa if (step >= length) { break; } - currentnode = gGeoManager->GetCurrentNode(); + currentnode = nav->GetCurrentNode(); if (!currentnode) { break; } length -= step; accountMaterial(currentnode->GetVolume()->GetMedium()->GetMaterial(), budStep); - gGeoManager->FindNextBoundaryAndStep(length, kFALSE); - step = gGeoManager->GetStep(); + nav->FindNextBoundaryAndStep(length, kFALSE); + step = nav->GetStep(); } budTotal.meanRho /= stepTot; budTotal.length = stepTot; diff --git a/Detectors/Base/src/MatLayerCyl.cxx b/Detectors/Base/src/MatLayerCyl.cxx index 2efe60235b895..dee179a84ca06 100644 --- a/Detectors/Base/src/MatLayerCyl.cxx +++ b/Detectors/Base/src/MatLayerCyl.cxx @@ -123,7 +123,7 @@ void MatLayerCyl::populateFromTGeo(int ntrPerCell) } //________________________________________________________________________________ -void MatLayerCyl::populateFromTGeo(int ip, int iz, int ntrPerCell) +void MatLayerCyl::populateFromTGeo(int ip, int iz, int ntrPerCell, TGeoNavigator* nav) { /// populate cell with info extracted from TGeometry, using ntrPerCell test tracks per cell @@ -136,7 +136,7 @@ void MatLayerCyl::populateFromTGeo(int ip, int iz, int ntrPerCell) 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); + auto bud = o2::base::GeometryManager::meanMaterialBudget(rMin * cs, rMin * sn, zs - dzt, rMax * cs, rMax * sn, zs + dzt, nav); 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 c390c8d617326..2db82b03e7cc3 100644 --- a/Detectors/Base/src/MatLayerCylSet.cxx +++ b/Detectors/Base/src/MatLayerCylSet.cxx @@ -17,7 +17,16 @@ #ifndef GPUCA_ALIGPUCODE // this part is unvisible on GPU version #include "GPUCommonLogger.h" #include +#include #include "CommonUtils/TreeStreamRedirector.h" +#include +#include +#include +#include +#include +#include +#include +#include //#define _DBG_LOC_ // for local debugging only #endif // !GPUCA_ALIGPUCODE @@ -69,9 +78,26 @@ void MatLayerCylSet::addLayer(float rmin, float rmax, float zmax, float dz, floa } //________________________________________________________________________________ -void MatLayerCylSet::populateFromTGeo(int ntrPerCell) +int MatLayerCylSet::getNThreadsFromEnv() { - ///< populate layers, using ntrPerCell test tracks per cell + ///< number of threads requested via NTHREADS_MATBUD, or 1 if unset/invalid + const char* env = std::getenv("NTHREADS_MATBUD"); + if (!env) { + return 1; + } + int n = std::atoi(env); + if (n < 1) { + LOG(warning) << "Ignoring invalid NTHREADS_MATBUD=" << env; + return 1; + } + return n; +} + +//________________________________________________________________________________ +void MatLayerCylSet::populateFromTGeo(int ntrPerCell, int nThreads) +{ + ///< populate layers, using ntrPerCell test tracks per cell. + ///< nThreads < 0 takes the number of threads from the NTHREADS_MATBUD environment variable. assert(mConstructionMask == InProgress); int nlr = getNLayers(); @@ -83,12 +109,86 @@ void MatLayerCylSet::populateFromTGeo(int ntrPerCell) LOG(error) << "The LUT is already populated"; return; } + + if (nThreads < 0) { + nThreads = getNThreadsFromEnv(); + } + + using Clock = std::chrono::steady_clock; + auto seconds = [](Clock::time_point a, Clock::time_point b) { + return std::chrono::duration(b - a).count(); + }; + + if (nThreads <= 1) { + for (int i = 0; i < nlr; i++) { + LOG(info) << "Populating with " << ntrPerCell << " trials Lr " << i; + get()->mLayers[i].print(); + } + const auto tFillStart = Clock::now(); + for (int i = 0; i < nlr; i++) { + get()->mLayers[i].populateFromTGeo(ntrPerCell); + } + const auto tFillEnd = Clock::now(); + finalizeStructures(); + LOG(info) << "LUT fill: 1 thread, cells " << seconds(tFillStart, tFillEnd) << " s"; + return; + } + + // Cells of all layers form one flat index range so that the load is balanced across + // threads even though layers differ a lot in cell count. layerOffsets[i] is the first + // flat index of layer i; a binary search maps a flat index back to (layer, iz, iphi). + std::vector layerOffsets(nlr + 1, 0); for (int i = 0; i < nlr; i++) { - printf("Populating with %d trials Lr %3d ", ntrPerCell, i); + LOG(info) << "Queuing " << ntrPerCell << " trials Lr " << i; get()->mLayers[i].print(); - get()->mLayers[i].populateFromTGeo(ntrPerCell); + const auto& lr = get()->mLayers[i]; + layerOffsets[i + 1] = layerOffsets[i] + size_t(lr.getNZBins()) * lr.getNPhiBins(); } + const size_t totalCells = layerOffsets[nlr]; + + 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); + } + finalizeStructures(); + const auto tEnd = Clock::now(); + + // Reported separately because only the middle term scales: the setup walks every volume + // in the geometry (TGeoManager::SetMaxThreads) and the teardown is serial by nature. + LOG(info) << "LUT fill: " << nThreads << " threads, setup " << seconds(tSetupStart, tFillStart) + << " s, cells " << seconds(tFillStart, tFillEnd) + << " s, finalize " << seconds(tFillEnd, tEnd) << " s"; } //________________________________________________________________________________ diff --git a/Detectors/Base/test/README.md b/Detectors/Base/test/README.md index f5f9fd4c04b29..97e8c7f569fd1 100644 --- a/Detectors/Base/test/README.md +++ b/Detectors/Base/test/README.md @@ -13,6 +13,16 @@ 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: +``` +export NTHREADS_MATBUD=16 +``` +The result does not depend on the number of threads. Scaling beyond a few threads needs +ROOT >= v6-36-10-alice3, which removes the per-query thread-id lookup and the false sharing +between the per-thread scratch buffers of TGeo shapes; with older ROOT the parallel path is +still correct, just slower. + 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 860fcbd5da940..44019198685f8 100644 --- a/Detectors/Base/test/buildMatBudLUT.C +++ b/Detectors/Base/test/buildMatBudLUT.C @@ -27,7 +27,10 @@ o2::base::MatLayerCylSet mbLUT; bool testMBLUT(const std::string& lutFile = "matbud.root"); -bool buildMatBudLUT(int nTst = 60, int maxLr = -1, const std::string& outFile = "matbud.root", const std::string& geomName = "o2sim_geometry-aligned.root"); +/// Build the material budget LUT. nThreads < 0 takes the thread count from NTHREADS_MATBUD. +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); struct LrData { float rMin = 0.f; @@ -42,7 +45,8 @@ struct LrData { std::vector lrData; void configLayers(); -bool buildMatBudLUT(int nTst, int maxLr, const std::string& outFile, const std::string& geomNamePrefix, const std::string& opts) +bool buildMatBudLUT(int nTst, int maxLr, const std::string& outFile, const std::string& geomNamePrefix, + const std::string& opts, int nThreads) { auto geomName = o2::base::NameConf::getGeomFileName(geomNamePrefix); if (gSystem->AccessPathName(geomName.c_str())) { // if needed, create geometry @@ -67,7 +71,7 @@ bool buildMatBudLUT(int nTst, int maxLr, const std::string& outFile, const std:: } TStopwatch sw; - mbLUT.populateFromTGeo(nTst); + mbLUT.populateFromTGeo(nTst, nThreads); mbLUT.optimizePhiSlices(); // move to populateFromTGeo mbLUT.flatten(); // move to populateFromTGeo diff --git a/Detectors/Base/test/compareMatBudLUT.C b/Detectors/Base/test/compareMatBudLUT.C new file mode 100644 index 0000000000000..7abaa9bd7182a --- /dev/null +++ b/Detectors/Base/test/compareMatBudLUT.C @@ -0,0 +1,93 @@ +// 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. + +/// \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: +/// +/// root -b -q 'compareMatBudLUT.C("matbud_serial.root","matbud_parallel.root")' + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include "DetectorsBase/MatLayerCylSet.h" +#include "GPUCommonLogger.h" +#include +#include +#endif + +/// Returns true if the two LUTs agree everywhere within tol (relative). +bool compareMatBudLUT(const std::string& fileA = "matbud_serial.root", + const std::string& fileB = "matbud_parallel.root", + float tol = 0.f) +{ + auto* lutA = o2::base::MatLayerCylSet::loadFromFile(fileA); + auto* lutB = o2::base::MatLayerCylSet::loadFromFile(fileB); + if (!lutA) { + LOG(error) << "Failed to load LUT from " << fileA; + return false; + } + if (!lutB) { + LOG(error) << "Failed to load LUT from " << fileB; + return false; + } + + if (lutA->getNLayers() != lutB->getNLayers()) { + LOG(error) << "Layer count differs: " << lutA->getNLayers() << " vs " << lutB->getNLayers(); + return false; + } + + size_t nCells = 0, nBad = 0; + double maxRelRho = 0., maxRelX2X0 = 0.; + + for (int il = 0; il < lutA->getNLayers(); il++) { + const auto& la = lutA->getLayer(il); + const auto& lb = lutB->getLayer(il); + if (la.getNZBins() != lb.getNZBins() || la.getNPhiBins() != lb.getNPhiBins()) { + LOG(error) << "Layer " << il << " segmentation differs: " + << la.getNZBins() << "x" << la.getNPhiBins() << " vs " + << lb.getNZBins() << "x" << lb.getNPhiBins(); + return false; + } + for (int iz = 0; iz < la.getNZBins(); iz++) { + for (int ip = 0; ip < la.getNPhiBins(); ip++) { + const auto& ca = la.getCellPhiBin(ip, iz); + const auto& cb = lb.getCellPhiBin(ip, iz); + nCells++; + + auto rel = [](float a, float b) { + const float den = std::max(std::abs(a), std::abs(b)); + return den > 0.f ? std::abs(a - b) / den : 0.f; + }; + const double rRho = rel(ca.meanRho, cb.meanRho); + const double rX = rel(ca.meanX2X0, cb.meanX2X0); + maxRelRho = std::max(maxRelRho, rRho); + maxRelX2X0 = std::max(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); + } + nBad++; + } + } + } + } + + printf("Compared %zu cells over %d layers\n", nCells, lutA->getNLayers()); + printf("Max relative difference: meanRho %.3g, meanX2X0 %.3g (tolerance %.3g)\n", maxRelRho, maxRelX2X0, tol); + if (nBad) { + LOG(error) << nBad << " cells differ beyond tolerance"; + return false; + } + LOG(info) << "LUTs agree"; + return true; +}