From f8085f6171e09c12aed7da7d7ffea0119252bba0 Mon Sep 17 00:00:00 2001 From: Tristan Wenzel Date: Tue, 28 Jul 2026 10:13:57 +0200 Subject: [PATCH 1/3] [geom] Lock the navigator map in all readers, not just the writer This commit improves multithreaded TGeo navigation. The changes come from profiling a parallel geometry scan in ALICE. AddNavigator() inserts into fNavigators under fgMutex, but GetCurrentNavigator(), GetListOfNavigators() and SetCurrentNavigator() called find() on the same std::map with no lock. A thread booking its navigator rebalances the tree underneath a concurrent reader, returning a garbage TGeoNavigatorArray* or crashing outright. The lookup is factored into FindNavigatorArray(), which takes fgMutex in MT mode and skips it otherwise, matching what AddNavigator() already does. The TLS cache in GetCurrentNavigator() is untouched, so steady-state navigation stays lock-free. Assisted-by: Claude Code (review, hardening and benchmarking) Supervised-by: Sandro Wenzel --- geom/geom/inc/TGeoManager.h | 1 + geom/geom/src/TGeoManager.cxx | 43 ++++++++++++++++++++++------------- 2 files changed, 28 insertions(+), 16 deletions(-) diff --git a/geom/geom/inc/TGeoManager.h b/geom/geom/inc/TGeoManager.h index f34ce7802aa66..faffa96b929ae 100644 --- a/geom/geom/inc/TGeoManager.h +++ b/geom/geom/inc/TGeoManager.h @@ -148,6 +148,7 @@ class TGeoManager : public TNamed { TGeoParallelWorld *fParallelWorld; // Parallel world ConstPropMap_t fProperties; // Map of user-defined constant properties //--- private methods + TGeoNavigatorArray *FindNavigatorArray() const; Bool_t IsLoopingVolumes() const { return fLoopVolumes; } void Init(); Bool_t InitArrayPNE() const; diff --git a/geom/geom/src/TGeoManager.cxx b/geom/geom/src/TGeoManager.cxx index 7fb21fd64d633..e672a8acc6acf 100644 --- a/geom/geom/src/TGeoManager.cxx +++ b/geom/geom/src/TGeoManager.cxx @@ -809,6 +809,25 @@ TGeoNavigator *TGeoManager::AddNavigator() return nav; } +//////////////////////////////////////////////////////////////////////////////// +/// Look up the navigator array booked for the calling thread, or nullptr if the +/// thread has none. +/// +/// AddNavigator() inserts into fNavigators under fgMutex, so in MT mode every +/// reader has to take the same lock: an unsynchronized find() racing an insert +/// walks a tree that is being rebalanced and returns garbage or crashes. In +/// single-threaded mode the map is never mutated concurrently and the lock is +/// skipped, as elsewhere in this class. + +TGeoNavigatorArray *TGeoManager::FindNavigatorArray() const +{ + std::unique_lock guard(fgMutex, std::defer_lock); + if (fMultiThread) + guard.lock(); + NavigatorsMap_t::const_iterator it = fNavigators.find(std::this_thread::get_id()); + return (it == fNavigators.end()) ? nullptr : it->second; +} + //////////////////////////////////////////////////////////////////////////////// /// Returns current navigator for the calling thread. @@ -820,11 +839,10 @@ TGeoNavigator *TGeoManager::GetCurrentNavigator() const TGeoNavigator *nav = tnav; // TTHREAD_TLS_GET(TGeoNavigator*,tnav); if (nav) return nav; - std::thread::id threadId = std::this_thread::get_id(); - NavigatorsMap_t::const_iterator it = fNavigators.find(threadId); - if (it == fNavigators.end()) + // The TLS cache above means the lookup below happens at most once per thread. + TGeoNavigatorArray *array = FindNavigatorArray(); + if (!array) return nullptr; - TGeoNavigatorArray *array = it->second; nav = array->GetCurrentNavigator(); tnav = nav; // TTHREAD_TLS_SET(TGeoNavigator*,tnav,nav); return nav; @@ -835,12 +853,7 @@ TGeoNavigator *TGeoManager::GetCurrentNavigator() const TGeoNavigatorArray *TGeoManager::GetListOfNavigators() const { - std::thread::id threadId = std::this_thread::get_id(); - NavigatorsMap_t::const_iterator it = fNavigators.find(threadId); - if (it == fNavigators.end()) - return nullptr; - TGeoNavigatorArray *array = it->second; - return array; + return FindNavigatorArray(); } //////////////////////////////////////////////////////////////////////////////// @@ -848,18 +861,16 @@ TGeoNavigatorArray *TGeoManager::GetListOfNavigators() const Bool_t TGeoManager::SetCurrentNavigator(Int_t index) { - std::thread::id threadId = std::this_thread::get_id(); - NavigatorsMap_t::const_iterator it = fNavigators.find(threadId); - if (it == fNavigators.end()) { + TGeoNavigatorArray *array = FindNavigatorArray(); + if (!array) { Error("SetCurrentNavigator", "No navigator defined for this thread\n"); - std::cout << " thread id: " << threadId << std::endl; + std::cout << " thread id: " << std::this_thread::get_id() << std::endl; return kFALSE; } - TGeoNavigatorArray *array = it->second; TGeoNavigator *nav = array->SetCurrentNavigator(index); if (!nav) { Error("SetCurrentNavigator", "Navigator %d not existing for this thread\n", index); - std::cout << " thread id: " << threadId << std::endl; + std::cout << " thread id: " << std::this_thread::get_id() << std::endl; return kFALSE; } if (!fMultiThread) From a4e0a956ac713845c35265aa9315b3087b83a08f Mon Sep 17 00:00:00 2001 From: Tristan Wenzel Date: Tue, 28 Jul 2026 10:13:57 +0200 Subject: [PATCH 2/3] [geom] Index per-thread shape data by object instead of by thread id This commit improves multithreaded TGeo navigation. The changes come from profiling a parallel geometry scan in ALICE: filling the material budget LUT on 28 cores goes from 12x to 23x speedup (139 s -> 72 s). Two costs dominate. Every per-thread scratch lookup went through TGeoManager::ThreadId(), a non-inlined cross-library __tls_get_addr call, paid on every boolean, section or division query. And the per-object ThreadData_t blocks are small allocations made back to back, so slots of different threads share a cache line and invalidate each other on every write -- false sharing, worst across sockets. Each object now gets a dense index and the per-thread state lives in one thread_local vector indexed by it. GetThreadData() becomes a header-inline TLS read plus an indexed load, and each thread owns its whole vector. TGeoBoolNode, TGeoVolumeAssembly, TGeoPatternFinder, TGeoPgon, TGeoXtru TGeoPgon/TGeoXtru: noexcept move ctor on ThreadData_t (owns heap buffers, and for Xtru a TGeoPolygon aliasing them) so a slot survives vector growth TGeoPatternFinder: matrix reused across generations -- it is owned by the geometry manager and never released, so a fresh one would leak per (thread, finder) Provisioning goes away: ClearThreadData() bumps a generation counter and each thread rebuilds its slot lazily, so CreateThreadData()/SetMaxThreads() are no longer needed to make a given thread count work -- any thread count works out of the box. The counters are atomic because ClearThreadData() is const and can run concurrently. Assisted-by: Claude Code (review, hardening and benchmarking) Supervised-by: Sandro Wenzel --- geom/geom/inc/TGeoBoolNode.h | 31 ++++---- geom/geom/inc/TGeoPatternFinder.h | 52 +++++++++----- geom/geom/inc/TGeoPgon.h | 48 ++++++++++--- geom/geom/inc/TGeoVolume.h | 39 +++++----- geom/geom/inc/TGeoXtru.h | 57 +++++++++++---- geom/geom/src/TGeoBoolNode.cxx | 78 +------------------- geom/geom/src/TGeoPatternFinder.cxx | 66 ++++++----------- geom/geom/src/TGeoPgon.cxx | 71 ++++++++----------- geom/geom/src/TGeoVolume.cxx | 87 +---------------------- geom/geom/src/TGeoXtru.cxx | 106 +++++++++++++++------------- 10 files changed, 267 insertions(+), 368 deletions(-) diff --git a/geom/geom/inc/TGeoBoolNode.h b/geom/geom/inc/TGeoBoolNode.h index 92ffb79e6d6bd..3bf6d776a489f 100644 --- a/geom/geom/inc/TGeoBoolNode.h +++ b/geom/geom/inc/TGeoBoolNode.h @@ -14,7 +14,8 @@ #include "TGeoShape.h" -#include +#include +#include #include // forward declarations @@ -23,17 +24,26 @@ class TGeoMatrix; class TGeoHMatrix; class TGeoBoolNode : public TObject { + static std::atomic fgInstanceCount; //! source of dense per-object indices + UInt_t fIndex{fgInstanceCount++}; //! dense index of this node into the per-thread vector + public: enum EGeoBoolType { kGeoUnion, kGeoIntersection, kGeoSubtraction }; struct ThreadData_t { - Int_t fSelected; // ! selected branch - - ThreadData_t(); - ~ThreadData_t(); + Int_t fSelected{0}; //! selected branch }; - ThreadData_t &GetThreadData() const; - void ClearThreadData() const; - void CreateThreadData(Int_t nthreads); + + /// Per-thread scratch state, owned by the calling thread and indexed by this node. + /// Each thread owns its whole vector, so no two threads ever write the same cache line. + ThreadData_t &GetThreadData() const + { + thread_local std::vector tdata; + if (tdata.size() <= fIndex) + tdata.resize(std::max(fgInstanceCount.load(std::memory_order_relaxed), fIndex + 1)); + return tdata[fIndex]; + } + void ClearThreadData() const {} + void CreateThreadData(Int_t) {} private: TGeoBoolNode(const TGeoBoolNode &) = delete; @@ -47,10 +57,7 @@ class TGeoBoolNode : public TObject { Int_t fNpoints{0}; //! number of points on the mesh Double_t *fPoints{nullptr}; //! array of mesh points - mutable std::vector fThreadData; //! Navigation data per thread - mutable Int_t fThreadSize{0}; //! Size for the navigation data array - mutable std::mutex fMutex; //! Mutex for thread data access - // methods + // methods Bool_t MakeBranch(const char *expr, Bool_t left); void AssignPoints(Int_t npoints, Double_t *points); diff --git a/geom/geom/inc/TGeoPatternFinder.h b/geom/geom/inc/TGeoPatternFinder.h index 7e3c441aae02d..3eafabf7ad07f 100644 --- a/geom/geom/inc/TGeoPatternFinder.h +++ b/geom/geom/inc/TGeoPatternFinder.h @@ -14,7 +14,8 @@ #include "TObject.h" -#include +#include +#include #include #include "TGeoVolume.h" @@ -23,24 +24,43 @@ class TGeoMatrix; /// base finder class for patterns. A pattern is specifying a division type class TGeoPatternFinder : public TObject { + static std::atomic fgInstanceCount; //! source of dense per-object indices + UInt_t fIndex{fgInstanceCount++}; //! dense index of this finder into the per-thread vector + mutable std::atomic fGeneration{0}; //! bumped whenever the per-thread state must be rebuilt + public: struct ThreadData_t { - TGeoMatrix *fMatrix; //! generic matrix - Int_t fCurrent; //! current division element - Int_t fNextIndex; //! index of next node - - ThreadData_t(); - ~ThreadData_t(); - - private: - ThreadData_t(const ThreadData_t &) = delete; - ThreadData_t &operator=(const ThreadData_t &) = delete; + // fMatrix is produced by the virtual CreateMatrix() and registered with (owned by) the + // geometry manager, so this struct does not own it. All members are trivially movable, + // which lets the slot live in a resizable thread_local vector. + TGeoMatrix *fMatrix{nullptr}; //! generic matrix (owned by TGeoManager) + Int_t fCurrent{-1}; //! current division element + Int_t fNextIndex{-1}; //! index of next node + Int_t fInitGen{-1}; //! generation this slot was last initialized for }; - ThreadData_t &GetThreadData() const; - void ClearThreadData() const; - void CreateThreadData(Int_t nthreads); + + /// Per-thread scratch state, owned by the calling thread and indexed by this finder. + /// Hot path: a TLS read plus an indexed load; the cold rebuild lives in InitThreadSlot(). + ThreadData_t &GetThreadData() const + { + thread_local std::vector tdata; + if (tdata.size() <= fIndex) + tdata.resize(std::max(fgInstanceCount.load(std::memory_order_relaxed), fIndex + 1)); + ThreadData_t &td = tdata[fIndex]; + if (td.fInitGen != fGeneration.load(std::memory_order_acquire)) + InitThreadSlot(td); + return td; + } + /// Invalidate the per-thread data. Each thread rebuilds its own slot lazily on next access, + /// so no cross-thread reach-in is needed. + void ClearThreadData() const { fGeneration.fetch_add(1, std::memory_order_release); } + /// No-op: per-thread data is allocated lazily, so no provisioning for a fixed thread count + /// is required and any number of threads works. + void CreateThreadData(Int_t) {} protected: + void InitThreadSlot(ThreadData_t &td) const; + enum EGeoPatternFlags { kPatternReflected = BIT(14), kPatternSpacedOut = BIT(15) }; Double_t fStep; // division step length Double_t fStart; // starting point on divided axis @@ -49,10 +69,6 @@ class TGeoPatternFinder : public TObject { Int_t fDivIndex; // index of first div. node TGeoVolume *fVolume; // volume to which applies - mutable std::vector fThreadData; //! Vector of thread private transient data - mutable Int_t fThreadSize; //! Size of the thread vector - mutable std::mutex fMutex; //! Mutex for thread data - protected: TGeoPatternFinder(const TGeoPatternFinder &); TGeoPatternFinder &operator=(const TGeoPatternFinder &); diff --git a/geom/geom/inc/TGeoPgon.h b/geom/geom/inc/TGeoPgon.h index 9a31281250317..dcc0badc3a976 100644 --- a/geom/geom/inc/TGeoPgon.h +++ b/geom/geom/inc/TGeoPgon.h @@ -14,28 +14,54 @@ #include "TGeoPcon.h" -#include +#include +#include #include class TGeoPgon : public TGeoPcon { + static std::atomic fgInstanceCount; //! source of dense per-object indices + UInt_t fIndex{fgInstanceCount++}; //! dense index of this shape into the per-thread vector + mutable std::atomic fGeneration{0}; //! bumped whenever the per-thread state must be rebuilt + public: struct ThreadData_t { - Int_t *fIntBuffer; //![fNedges+4] temporary int buffer array - Double_t *fDblBuffer; //![fNedges+4] temporary double buffer array + Int_t *fIntBuffer{nullptr}; //![fNedges+4] temporary int buffer array + Double_t *fDblBuffer{nullptr}; //![fNedges+4] temporary double buffer array + Int_t fInitGen{-1}; //! generation this slot was last initialized for - ThreadData_t(); + ThreadData_t() = default; ~ThreadData_t(); + // Owns the two scratch buffers: movable so the slot can live in a resizable vector, + // but not copyable. + ThreadData_t(ThreadData_t &&other) noexcept; + ThreadData_t &operator=(ThreadData_t &&other) noexcept; + ThreadData_t(const ThreadData_t &) = delete; + ThreadData_t &operator=(const ThreadData_t &) = delete; }; - ThreadData_t &GetThreadData() const; - void ClearThreadData() const override; - void CreateThreadData(Int_t nthreads) override; + + /// Per-thread scratch buffers, owned by the calling thread and indexed by this shape. + /// Hot path: a TLS read plus an indexed load; the cold rebuild lives in InitThreadSlot(). + ThreadData_t &GetThreadData() const + { + thread_local std::vector tdata; + if (tdata.size() <= fIndex) + tdata.resize(std::max(fgInstanceCount.load(std::memory_order_relaxed), fIndex + 1)); + ThreadData_t &td = tdata[fIndex]; + if (td.fInitGen != fGeneration.load(std::memory_order_acquire)) + InitThreadSlot(td); + return td; + } + /// Invalidate the per-thread data. Each thread rebuilds its own slot lazily on next access. + void ClearThreadData() const override { fGeneration.fetch_add(1, std::memory_order_release); } + /// No-op: per-thread data is allocated lazily, so no provisioning for a fixed thread count + /// is required and any number of threads works. + void CreateThreadData(Int_t) override {} protected: + void InitThreadSlot(ThreadData_t &td) const; + // data members - Int_t fNedges; // number of edges (at least one) - mutable std::vector fThreadData; //! Navigation data per thread - mutable Int_t fThreadSize; //! Size for the navigation data array - mutable std::mutex fMutex; //! Mutex for thread data + Int_t fNedges; // number of edges (at least one) // internal utility methods Int_t GetPhiCrossList(const Double_t *point, const Double_t *dir, Int_t istart, Double_t *sphi, Int_t *iphi, diff --git a/geom/geom/inc/TGeoVolume.h b/geom/geom/inc/TGeoVolume.h index 17b6716e59d71..b716b01f2e12f 100644 --- a/geom/geom/inc/TGeoVolume.h +++ b/geom/geom/inc/TGeoVolume.h @@ -21,6 +21,8 @@ #include "TObjArray.h" #include "TGeoMedium.h" #include "TGeoShape.h" +#include +#include #include #include @@ -314,23 +316,32 @@ class TGeoVolumeMulti : public TGeoVolume { //////////////////////////////////////////////////////////////////////////// class TGeoVolumeAssembly : public TGeoVolume { + static std::atomic fgInstanceCount; //! source of dense per-object indices + UInt_t fIndex{fgInstanceCount++}; //! dense index of this assembly into the per-thread vector + public: struct ThreadData_t { - Int_t fCurrent; //! index of current selected node - Int_t fNext; //! index of next node to be entered - - ThreadData_t(); - ~ThreadData_t(); + Int_t fCurrent{-1}; //! index of current selected node + Int_t fNext{-1}; //! index of next node to be entered }; - ThreadData_t &GetThreadData() const; - void ClearThreadData() const override; - void CreateThreadData(Int_t nthreads) override; + /// Per-thread scratch state, owned by the calling thread and indexed by this assembly. + /// Each thread owns its whole vector, so no two threads ever write the same cache line. + ThreadData_t &GetThreadData() const + { + thread_local std::vector tdata; + if (tdata.size() <= fIndex) + tdata.resize(std::max(fgInstanceCount.load(std::memory_order_relaxed), fIndex + 1)); + return tdata[fIndex]; + } + // ClearThreadData()/CreateThreadData() are deliberately not overridden: the assembly's own + // per-thread state is allocated lazily, and TGeoVolume's implementation still has to run so + // the shape and the division finder get their generation bumped. -protected: - mutable std::vector fThreadData; //! Thread specific data vector - mutable Int_t fThreadSize; //! Thread vector size - mutable std::mutex fMutex; //! Mutex for concurrent operations + Int_t GetCurrentNodeIndex() const override { return GetThreadData().fCurrent; } + Int_t GetNextNodeIndex() const override { return GetThreadData().fNext; } + void SetCurrentNodeIndex(Int_t index) { GetThreadData().fCurrent = index; } + void SetNextNodeIndex(Int_t index) { GetThreadData().fNext = index; } private: TGeoVolumeAssembly(const TGeoVolumeAssembly &) = delete; @@ -348,13 +359,9 @@ class TGeoVolumeAssembly : public TGeoVolume { Option_t *option = "") override; TGeoVolume *Divide(TGeoVolume *cell, TGeoPatternFinder *pattern, Option_t *option = "spacedout"); void DrawOnly(Option_t *) override {} - Int_t GetCurrentNodeIndex() const override; - Int_t GetNextNodeIndex() const override; Bool_t IsAssembly() const override { return kTRUE; } Bool_t IsVisible() const override { return kFALSE; } static TGeoVolumeAssembly *MakeAssemblyFromVolume(TGeoVolume *vol); - void SetCurrentNodeIndex(Int_t index); - void SetNextNodeIndex(Int_t index); ClassDefOverride(TGeoVolumeAssembly, 2) // an assembly of volumes }; diff --git a/geom/geom/inc/TGeoXtru.h b/geom/geom/inc/TGeoXtru.h index dde7c5073866d..ae414b2cecf46 100644 --- a/geom/geom/inc/TGeoXtru.h +++ b/geom/geom/inc/TGeoXtru.h @@ -14,28 +14,61 @@ #include "TGeoBBox.h" -#include +#include +#include #include class TGeoPolygon; class TGeoXtru : public TGeoBBox { + static std::atomic fgInstanceCount; //! source of dense per-object indices + UInt_t fIndex{fgInstanceCount++}; //! dense index of this shape into the per-thread vector + mutable std::atomic fGeneration{0}; //! bumped whenever the per-thread state must be rebuilt + mutable std::atomic fIllegalChecked{kFALSE}; //! illegal-polygon warning already emitted + public: struct ThreadData_t { - Int_t fSeg; // !current segment [0,fNvert-1] - Int_t fIz; // !current z plane [0,fNz-1] - Double_t *fXc; // ![fNvert] current X positions for polygon vertices - Double_t *fYc; // ![fNvert] current Y positions for polygon vertices - TGeoPolygon *fPoly; // !polygon defining section shape + Int_t fSeg{0}; //! current segment [0,fNvert-1] + Int_t fIz{0}; //! current z plane [0,fNz-1] + Double_t *fXc{nullptr}; //![fNvert] current X positions for polygon vertices + Double_t *fYc{nullptr}; //![fNvert] current Y positions for polygon vertices + TGeoPolygon *fPoly{nullptr}; //! polygon defining section shape + Int_t fInitGen{-1}; //! generation this slot was last initialized for - ThreadData_t(); + ThreadData_t() = default; ~ThreadData_t(); + // Owns fXc/fYc/fPoly: movable so the slot can live in a resizable vector, not copyable. + ThreadData_t(ThreadData_t &&other) noexcept; + ThreadData_t &operator=(ThreadData_t &&other) noexcept; + ThreadData_t(const ThreadData_t &) = delete; + ThreadData_t &operator=(const ThreadData_t &) = delete; }; - ThreadData_t &GetThreadData() const; - void ClearThreadData() const override; - void CreateThreadData(Int_t nthreads) override; + + /// Per-thread scratch state, owned by the calling thread and indexed by this shape. + /// Hot path: a TLS read plus an indexed load; the cold rebuild lives in InitThreadSlot(). + ThreadData_t &GetThreadData() const + { + thread_local std::vector tdata; + if (tdata.size() <= fIndex) + tdata.resize(std::max(fgInstanceCount.load(std::memory_order_relaxed), fIndex + 1)); + ThreadData_t &td = tdata[fIndex]; + if (td.fInitGen != fGeneration.load(std::memory_order_acquire)) + InitThreadSlot(td); + return td; + } + /// Invalidate the per-thread data. Each thread rebuilds its own slot lazily on next access. + void ClearThreadData() const override + { + fGeneration.fetch_add(1, std::memory_order_release); + fIllegalChecked.store(kFALSE, std::memory_order_relaxed); + } + /// No-op: per-thread data is allocated lazily, so no provisioning for a fixed thread count + /// is required and any number of threads works. + void CreateThreadData(Int_t) override {} protected: + void InitThreadSlot(ThreadData_t &td) const; + // data members Int_t fNvert; // number of vertices of the 2D polygon (at least 3) Int_t fNz; // number of z planes (at least two) @@ -47,10 +80,6 @@ class TGeoXtru : public TGeoBBox { Double_t *fX0; //[fNz] array of X offsets (for each Z) Double_t *fY0; //[fNz] array of Y offsets (for each Z) - mutable std::vector fThreadData; //! Navigation data per thread - mutable Int_t fThreadSize; //! size of thread-specific array - mutable std::mutex fMutex; //! mutex for thread data - TGeoXtru(const TGeoXtru &) = delete; TGeoXtru &operator=(const TGeoXtru &) = delete; diff --git a/geom/geom/src/TGeoBoolNode.cxx b/geom/geom/src/TGeoBoolNode.cxx index c84b4ab581e15..0cdaa7e5936d9 100644 --- a/geom/geom/src/TGeoBoolNode.cxx +++ b/geom/geom/src/TGeoBoolNode.cxx @@ -43,76 +43,9 @@ implementations for Boolean nodes are: - TGeoIntersection - representing the Boolean intersection of two positioned shapes */ -ClassImp(TGeoBoolNode); - -//////////////////////////////////////////////////////////////////////////////// -/// Constructor. - -TGeoBoolNode::ThreadData_t::ThreadData_t() : fSelected(0) {} - -//////////////////////////////////////////////////////////////////////////////// -/// Destructor. - -TGeoBoolNode::ThreadData_t::~ThreadData_t() {} - -//////////////////////////////////////////////////////////////////////////////// - -TGeoBoolNode::ThreadData_t &TGeoBoolNode::GetThreadData() const -{ - Int_t tid = TGeoManager::ThreadId(); - /* - std::lock_guard guard(fMutex); - if (tid >= fThreadSize) { - Error("GetThreadData", "Thread id=%d bigger than maximum declared thread number %d. \nUse - TGeoManager::SetMaxThreads properly !!!", tid, fThreadSize); - } - if (tid >= fThreadSize) - { - fThreadData.resize(tid + 1); - fThreadSize = tid + 1; - } - if (fThreadData[tid] == 0) - { - if (fThreadData[tid] == 0) - fThreadData[tid] = new ThreadData_t; - } - */ - return *fThreadData[tid]; -} - -//////////////////////////////////////////////////////////////////////////////// - -void TGeoBoolNode::ClearThreadData() const -{ - std::lock_guard guard(fMutex); - std::vector::iterator i = fThreadData.begin(); - while (i != fThreadData.end()) { - delete *i; - ++i; - } - fThreadData.clear(); - fThreadSize = 0; -} +std::atomic TGeoBoolNode::fgInstanceCount{0}; -//////////////////////////////////////////////////////////////////////////////// -/// Create thread data for n threads max. - -void TGeoBoolNode::CreateThreadData(Int_t nthreads) -{ - std::lock_guard guard(fMutex); - fThreadData.resize(nthreads); - fThreadSize = nthreads; - for (Int_t tid = 0; tid < nthreads; tid++) { - if (fThreadData[tid] == nullptr) { - fThreadData[tid] = new ThreadData_t; - } - } - // Propagate to components - if (fLeft) - fLeft->CreateThreadData(nthreads); - if (fRight) - fRight->CreateThreadData(nthreads); -} +ClassImp(TGeoBoolNode); //////////////////////////////////////////////////////////////////////////////// /// Set the selected branch. @@ -133,8 +66,6 @@ TGeoBoolNode::TGeoBoolNode() fRightMat = nullptr; fNpoints = 0; fPoints = nullptr; - fThreadSize = 0; - CreateThreadData(1); } //////////////////////////////////////////////////////////////////////////////// @@ -148,8 +79,6 @@ TGeoBoolNode::TGeoBoolNode(const char *expr1, const char *expr2) fRightMat = nullptr; fNpoints = 0; fPoints = nullptr; - fThreadSize = 0; - CreateThreadData(1); if (!MakeBranch(expr1, kTRUE)) { return; } @@ -168,8 +97,6 @@ TGeoBoolNode::TGeoBoolNode(TGeoShape *left, TGeoShape *right, TGeoMatrix *lmat, fLeftMat = lmat; fNpoints = 0; fPoints = nullptr; - fThreadSize = 0; - CreateThreadData(1); if (!fLeftMat) fLeftMat = gGeoIdentity; else @@ -197,7 +124,6 @@ TGeoBoolNode::~TGeoBoolNode() { if (fPoints) delete[] fPoints; - ClearThreadData(); } //////////////////////////////////////////////////////////////////////////////// diff --git a/geom/geom/src/TGeoPatternFinder.cxx b/geom/geom/src/TGeoPatternFinder.cxx index aecf1d8e3cfac..4dfe6e92ae836 100644 --- a/geom/geom/src/TGeoPatternFinder.cxx +++ b/geom/geom/src/TGeoPatternFinder.cxx @@ -36,6 +36,8 @@ on different axis. Implemented patterns are: #include "TGeoManager.h" #include "TMath.h" +#include + ClassImp(TGeoPatternFinder); ClassImp(TGeoPatternX); ClassImp(TGeoPatternY); @@ -51,55 +53,29 @@ ClassImp(TGeoPatternSphTheta); ClassImp(TGeoPatternSphPhi); ClassImp(TGeoPatternHoneycomb); -//////////////////////////////////////////////////////////////////////////////// -/// Constructor. - -TGeoPatternFinder::ThreadData_t::ThreadData_t() : fMatrix(nullptr), fCurrent(-1), fNextIndex(-1) {} - -//////////////////////////////////////////////////////////////////////////////// -/// Destructor. - -TGeoPatternFinder::ThreadData_t::~ThreadData_t() -{ - // if (fMatrix != gGeoIdentity) delete fMatrix; -} - -//////////////////////////////////////////////////////////////////////////////// - -TGeoPatternFinder::ThreadData_t &TGeoPatternFinder::GetThreadData() const -{ - Int_t tid = TGeoManager::ThreadId(); - return *fThreadData[tid]; -} - -//////////////////////////////////////////////////////////////////////////////// - -void TGeoPatternFinder::ClearThreadData() const -{ - std::lock_guard guard(fMutex); - std::vector::iterator i = fThreadData.begin(); - while (i != fThreadData.end()) { - delete *i; - ++i; - } - fThreadData.clear(); - fThreadSize = 0; -} +std::atomic TGeoPatternFinder::fgInstanceCount{0}; //////////////////////////////////////////////////////////////////////////////// -/// Create thread data for n threads max. +/// (Re)build the per-thread scratch state for this finder into the given slot. +/// Cold path: runs once per (thread, finder, generation). -void TGeoPatternFinder::CreateThreadData(Int_t nthreads) +void TGeoPatternFinder::InitThreadSlot(ThreadData_t &td) const { - std::lock_guard guard(fMutex); - fThreadData.resize(nthreads); - fThreadSize = nthreads; - for (Int_t tid = 0; tid < nthreads; tid++) { - if (fThreadData[tid] == nullptr) { - fThreadData[tid] = new ThreadData_t; - fThreadData[tid]->fMatrix = CreateMatrix(); - } + if (!td.fMatrix) { + // CreateMatrix() registers the new matrix with the geometry manager, which mutates a + // shared, unlocked TObjArray. Lazy initialization means several threads can reach this + // on first touch concurrently, so serialize the registration. Taken once per + // (thread, finder); steady-state navigation is lock-free. + static std::mutex sInitMutex; + std::lock_guard guard(sInitMutex); + td.fMatrix = CreateMatrix(); } + // A generation bump only invalidates the cached division indices. The matrix stays valid and + // is deliberately reused: it is owned by the geometry manager and never released, so creating + // a fresh one here would leak one matrix per (thread, finder) on every ClearThreadData(). + td.fCurrent = -1; + td.fNextIndex = -1; + td.fInitGen = fGeneration.load(std::memory_order_acquire); } //////////////////////////////////////////////////////////////////////////////// @@ -113,7 +89,6 @@ TGeoPatternFinder::TGeoPatternFinder() fStart = 0; fEnd = 0; fVolume = nullptr; - fThreadSize = 0; } //////////////////////////////////////////////////////////////////////////////// @@ -127,7 +102,6 @@ TGeoPatternFinder::TGeoPatternFinder(TGeoVolume *vol, Int_t ndiv) fStep = 0; fStart = 0; fEnd = 0; - fThreadSize = 0; } //////////////////////////////////////////////////////////////////////////////// diff --git a/geom/geom/src/TGeoPgon.cxx b/geom/geom/src/TGeoPgon.cxx index d7a864aab7335..59bf9abd357af 100644 --- a/geom/geom/src/TGeoPgon.cxx +++ b/geom/geom/src/TGeoPgon.cxx @@ -66,12 +66,9 @@ polygons, between `phi1` and `phi1+dphi.` #include "TBuffer3DTypes.h" #include "TMath.h" -ClassImp(TGeoPgon); - -//////////////////////////////////////////////////////////////////////////////// -/// Constructor. +std::atomic TGeoPgon::fgInstanceCount{0}; -TGeoPgon::ThreadData_t::ThreadData_t() : fIntBuffer(nullptr), fDblBuffer(nullptr) {} +ClassImp(TGeoPgon); //////////////////////////////////////////////////////////////////////////////// /// Destructor. @@ -83,44 +80,45 @@ TGeoPgon::ThreadData_t::~ThreadData_t() } //////////////////////////////////////////////////////////////////////////////// +/// Move constructor. Steals the owned buffers; the moved-from slot is left empty +/// and uninitialized (fInitGen = -1). -TGeoPgon::ThreadData_t &TGeoPgon::GetThreadData() const +TGeoPgon::ThreadData_t::ThreadData_t(ThreadData_t &&other) noexcept + : fIntBuffer(other.fIntBuffer), fDblBuffer(other.fDblBuffer), fInitGen(other.fInitGen) { - Int_t tid = TGeoManager::ThreadId(); - return *fThreadData[tid]; + other.fIntBuffer = nullptr; + other.fDblBuffer = nullptr; + other.fInitGen = -1; } //////////////////////////////////////////////////////////////////////////////// +/// Move assignment. Releases the current buffers before stealing the source's. -void TGeoPgon::ClearThreadData() const +TGeoPgon::ThreadData_t &TGeoPgon::ThreadData_t::operator=(ThreadData_t &&other) noexcept { - std::lock_guard guard(fMutex); - std::vector::iterator i = fThreadData.begin(); - while (i != fThreadData.end()) { - delete *i; - ++i; - } - fThreadData.clear(); - fThreadSize = 0; + if (this != &other) { + delete[] fIntBuffer; + delete[] fDblBuffer; + fIntBuffer = other.fIntBuffer; + fDblBuffer = other.fDblBuffer; + fInitGen = other.fInitGen; + other.fIntBuffer = nullptr; + other.fDblBuffer = nullptr; + other.fInitGen = -1; + } + return *this; } //////////////////////////////////////////////////////////////////////////////// -/// Create thread data for n threads max. +/// (Re)build the per-thread scratch buffers for this shape into the given slot. +/// Cold path: runs once per (thread, shape, generation). -void TGeoPgon::CreateThreadData(Int_t nthreads) +void TGeoPgon::InitThreadSlot(ThreadData_t &td) const { - if (fThreadSize) - ClearThreadData(); - std::lock_guard guard(fMutex); - fThreadData.resize(nthreads); - fThreadSize = nthreads; - for (Int_t tid = 0; tid < nthreads; tid++) { - if (fThreadData[tid] == nullptr) { - fThreadData[tid] = new ThreadData_t; - fThreadData[tid]->fIntBuffer = new Int_t[fNedges + 10]; - fThreadData[tid]->fDblBuffer = new Double_t[fNedges + 10]; - } - } + td = ThreadData_t{}; // release any buffers left from a previous generation + td.fIntBuffer = new Int_t[fNedges + 10]; + td.fDblBuffer = new Double_t[fNedges + 10]; + td.fInitGen = fGeneration.load(std::memory_order_acquire); } //////////////////////////////////////////////////////////////////////////////// @@ -130,7 +128,6 @@ TGeoPgon::TGeoPgon() { SetShapeBit(TGeoShape::kGeoPgon); fNedges = 0; - fThreadSize = 0; } //////////////////////////////////////////////////////////////////////////////// @@ -140,8 +137,6 @@ TGeoPgon::TGeoPgon(Double_t phi, Double_t dphi, Int_t nedges, Int_t nz) : TGeoPc { SetShapeBit(TGeoShape::kGeoPgon); fNedges = nedges; - fThreadSize = 0; - CreateThreadData(1); } //////////////////////////////////////////////////////////////////////////////// @@ -152,8 +147,6 @@ TGeoPgon::TGeoPgon(const char *name, Double_t phi, Double_t dphi, Int_t nedges, { SetShapeBit(TGeoShape::kGeoPgon); fNedges = nedges; - fThreadSize = 0; - CreateThreadData(1); } //////////////////////////////////////////////////////////////////////////////// @@ -172,8 +165,6 @@ TGeoPgon::TGeoPgon(Double_t *param) : TGeoPcon("") SetShapeBit(TGeoShape::kGeoPgon); SetDimensions(param); ComputeBBox(); - fThreadSize = 0; - CreateThreadData(1); } //////////////////////////////////////////////////////////////////////////////// @@ -467,8 +458,6 @@ TGeoPgon::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact, ipl++; } Double_t stepmax = step; - if (!fThreadSize) - ((TGeoPgon *)this)->CreateThreadData(1); ThreadData_t &td = GetThreadData(); Double_t *sph = td.fDblBuffer; Int_t *iph = td.fIntBuffer; @@ -1254,8 +1243,6 @@ TGeoPgon::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact } } } - if (!fThreadSize) - ((TGeoPgon *)this)->CreateThreadData(1); ThreadData_t &td = GetThreadData(); Double_t *sph = td.fDblBuffer; Int_t *iph = td.fIntBuffer; diff --git a/geom/geom/src/TGeoVolume.cxx b/geom/geom/src/TGeoVolume.cxx index dfb0eecf7e14c..4ab65882aafdf 100644 --- a/geom/geom/src/TGeoVolume.cxx +++ b/geom/geom/src/TGeoVolume.cxx @@ -2859,93 +2859,14 @@ void TGeoVolumeMulti::SetVisibility(Bool_t vis) } } -ClassImp(TGeoVolumeAssembly); - -//////////////////////////////////////////////////////////////////////////////// -/// Constructor. - -TGeoVolumeAssembly::ThreadData_t::ThreadData_t() : fCurrent(-1), fNext(-1) {} - -//////////////////////////////////////////////////////////////////////////////// -/// Destructor. - -TGeoVolumeAssembly::ThreadData_t::~ThreadData_t() {} - -//////////////////////////////////////////////////////////////////////////////// - -TGeoVolumeAssembly::ThreadData_t &TGeoVolumeAssembly::GetThreadData() const -{ - Int_t tid = TGeoManager::ThreadId(); - return *fThreadData[tid]; -} +std::atomic TGeoVolumeAssembly::fgInstanceCount{0}; -//////////////////////////////////////////////////////////////////////////////// - -void TGeoVolumeAssembly::ClearThreadData() const -{ - std::lock_guard guard(fMutex); - TGeoVolume::ClearThreadData(); - std::vector::iterator i = fThreadData.begin(); - while (i != fThreadData.end()) { - delete *i; - ++i; - } - fThreadData.clear(); - fThreadSize = 0; -} - -//////////////////////////////////////////////////////////////////////////////// - -void TGeoVolumeAssembly::CreateThreadData(Int_t nthreads) -{ - std::lock_guard guard(fMutex); - // Create assembly thread data here - fThreadData.resize(nthreads); - fThreadSize = nthreads; - for (Int_t tid = 0; tid < nthreads; tid++) { - if (fThreadData[tid] == nullptr) { - fThreadData[tid] = new ThreadData_t; - } - } - TGeoVolume::CreateThreadData(nthreads); -} - -//////////////////////////////////////////////////////////////////////////////// - -Int_t TGeoVolumeAssembly::GetCurrentNodeIndex() const -{ - return fThreadData[TGeoManager::ThreadId()]->fCurrent; -} - -//////////////////////////////////////////////////////////////////////////////// - -Int_t TGeoVolumeAssembly::GetNextNodeIndex() const -{ - return fThreadData[TGeoManager::ThreadId()]->fNext; -} - -//////////////////////////////////////////////////////////////////////////////// - -void TGeoVolumeAssembly::SetCurrentNodeIndex(Int_t index) -{ - fThreadData[TGeoManager::ThreadId()]->fCurrent = index; -} - -//////////////////////////////////////////////////////////////////////////////// - -void TGeoVolumeAssembly::SetNextNodeIndex(Int_t index) -{ - fThreadData[TGeoManager::ThreadId()]->fNext = index; -} +ClassImp(TGeoVolumeAssembly); //////////////////////////////////////////////////////////////////////////////// /// Default constructor -TGeoVolumeAssembly::TGeoVolumeAssembly() : TGeoVolume() -{ - fThreadSize = 0; - CreateThreadData(1); -} +TGeoVolumeAssembly::TGeoVolumeAssembly() : TGeoVolume() {} //////////////////////////////////////////////////////////////////////////////// /// Constructor. Just the name has to be provided. Assemblies does not have their own @@ -2958,8 +2879,6 @@ TGeoVolumeAssembly::TGeoVolumeAssembly(const char *name) : TGeoVolume() fShape = new TGeoShapeAssembly(this); if (fGeoManager) fNumber = fGeoManager->AddVolume(this); - fThreadSize = 0; - CreateThreadData(1); } //////////////////////////////////////////////////////////////////////////////// diff --git a/geom/geom/src/TGeoXtru.cxx b/geom/geom/src/TGeoXtru.cxx index 9cce080f0b562..ad19cd639e199 100644 --- a/geom/geom/src/TGeoXtru.cxx +++ b/geom/geom/src/TGeoXtru.cxx @@ -101,13 +101,11 @@ Double_t y0, Double_t scale); #include "TGeoManager.h" #include "TGeoVolume.h" #include "TGeoPolygon.h" +#include "TROOT.h" -ClassImp(TGeoXtru); - -//////////////////////////////////////////////////////////////////////////////// -/// Constructor. +std::atomic TGeoXtru::fgInstanceCount{0}; -TGeoXtru::ThreadData_t::ThreadData_t() : fSeg(0), fIz(0), fXc(nullptr), fYc(nullptr), fPoly(nullptr) {} +ClassImp(TGeoXtru); //////////////////////////////////////////////////////////////////////////////// /// Destructor. @@ -116,57 +114,73 @@ TGeoXtru::ThreadData_t::~ThreadData_t() { delete[] fXc; delete[] fYc; - delete fPoly; + // fPoly is a TObject, so deleting it walks ROOT's cleanup machinery (RecursiveRemove, + // TObjArray::Delete). This slot lives in a thread_local vector and can therefore be + // destroyed at thread/process exit, possibly after ROOT's globals have been torn down, + // where that machinery reads freed state. Only delete while ROOT is still alive -- at + // shutdown the OS reclaims this small per-thread buffer anyway. The in-run rebuild path + // (move-assignment in InitThreadSlot) always runs while ROOT is up, so this guard only + // ever skips the very last teardown, never steady-state churn. + // Same "is ROOT still there" test that TDirectory and TGenericClassInfo use. + if (fPoly && ROOT::Internal::gROOTLocal) + delete fPoly; } //////////////////////////////////////////////////////////////////////////////// +/// Move constructor. Steals the owned resources; the moved-from slot is left empty and +/// uninitialized (fInitGen = -1). fPoly keeps pointing at fXc/fYc, which are not relocated. -TGeoXtru::ThreadData_t &TGeoXtru::GetThreadData() const +TGeoXtru::ThreadData_t::ThreadData_t(ThreadData_t &&other) noexcept + : fSeg(other.fSeg), fIz(other.fIz), fXc(other.fXc), fYc(other.fYc), fPoly(other.fPoly), fInitGen(other.fInitGen) { - if (!fThreadSize) - ((TGeoXtru *)this)->CreateThreadData(1); - Int_t tid = TGeoManager::ThreadId(); - return *fThreadData[tid]; + other.fXc = nullptr; + other.fYc = nullptr; + other.fPoly = nullptr; + other.fInitGen = -1; } //////////////////////////////////////////////////////////////////////////////// +/// Move assignment. Releases the current resources before stealing the source's. -void TGeoXtru::ClearThreadData() const +TGeoXtru::ThreadData_t &TGeoXtru::ThreadData_t::operator=(ThreadData_t &&other) noexcept { - std::lock_guard guard(fMutex); - std::vector::iterator i = fThreadData.begin(); - while (i != fThreadData.end()) { - delete *i; - ++i; + if (this != &other) { + delete[] fXc; + delete[] fYc; + delete fPoly; + fSeg = other.fSeg; + fIz = other.fIz; + fXc = other.fXc; + fYc = other.fYc; + fPoly = other.fPoly; + fInitGen = other.fInitGen; + other.fXc = nullptr; + other.fYc = nullptr; + other.fPoly = nullptr; + other.fInitGen = -1; } - fThreadData.clear(); - fThreadSize = 0; + return *this; } //////////////////////////////////////////////////////////////////////////////// -/// Create thread data for n threads max. +/// (Re)build the per-thread scratch state for this shape into the given slot. +/// Cold path: runs once per (thread, shape, generation). -void TGeoXtru::CreateThreadData(Int_t nthreads) +void TGeoXtru::InitThreadSlot(ThreadData_t &td) const { - std::lock_guard guard(fMutex); - fThreadData.resize(nthreads); - fThreadSize = nthreads; - for (Int_t tid = 0; tid < nthreads; tid++) { - if (fThreadData[tid] == nullptr) { - fThreadData[tid] = new ThreadData_t; - ThreadData_t &td = *fThreadData[tid]; - td.fXc = new Double_t[fNvert]; - td.fYc = new Double_t[fNvert]; - memcpy(td.fXc, fX, fNvert * sizeof(Double_t)); - memcpy(td.fYc, fY, fNvert * sizeof(Double_t)); - td.fPoly = new TGeoPolygon(fNvert); - td.fPoly->SetXY(td.fXc, td.fYc); // initialize with current coordinates - td.fPoly->FinishPolygon(); - if (tid == 0 && td.fPoly->IsIllegalCheck()) { - Error("DefinePolygon", "Shape %s of type XTRU has an illegal polygon.", GetName()); - } - } - } + td = ThreadData_t{}; // release anything left from a previous generation + td.fXc = new Double_t[fNvert]; + td.fYc = new Double_t[fNvert]; + memcpy(td.fXc, fX, fNvert * sizeof(Double_t)); + memcpy(td.fYc, fY, fNvert * sizeof(Double_t)); + td.fPoly = new TGeoPolygon(fNvert); + td.fPoly->SetXY(td.fXc, td.fYc); // initialize with current coordinates + td.fPoly->FinishPolygon(); + // The polygon is identical in every thread, so report an illegal one exactly once + // instead of once per thread (previously: only for thread id 0). + if (td.fPoly->IsIllegalCheck() && !fIllegalChecked.exchange(kTRUE, std::memory_order_relaxed)) + Error("DefinePolygon", "Shape %s of type XTRU has an illegal polygon.", GetName()); + td.fInitGen = fGeneration.load(std::memory_order_acquire); } //////////////////////////////////////////////////////////////////////////////// @@ -197,9 +211,7 @@ TGeoXtru::TGeoXtru() fZ(nullptr), fScale(nullptr), fX0(nullptr), - fY0(nullptr), - fThreadData(0), - fThreadSize(0) + fY0(nullptr) { SetShapeBit(TGeoShape::kGeoXtru); } @@ -217,9 +229,7 @@ TGeoXtru::TGeoXtru(Int_t nz) fZ(new Double_t[nz]), fScale(new Double_t[nz]), fX0(new Double_t[nz]), - fY0(new Double_t[nz]), - fThreadData(0), - fThreadSize(0) + fY0(new Double_t[nz]) { SetShapeBit(TGeoShape::kGeoXtru); if (nz < 2) { @@ -253,9 +263,7 @@ TGeoXtru::TGeoXtru(Double_t *param) fZ(nullptr), fScale(nullptr), fX0(nullptr), - fY0(nullptr), - fThreadData(0), - fThreadSize(0) + fY0(nullptr) { SetShapeBit(TGeoShape::kGeoXtru); SetDimensions(param); From 5fe298d8fda047c66f3d2d6c538074639b64f0e3 Mon Sep 17 00:00:00 2001 From: Tristan Wenzel Date: Tue, 28 Jul 2026 10:13:57 +0200 Subject: [PATCH 3/3] [geom] Add a multi-threaded navigation regression test This commit improves multithreaded TGeo navigation. The changes come from profiling a parallel geometry scan in ALICE. Eight threads navigating the same geometry must give exactly the single-threaded answer. Covers TGeoXtru, TGeoPgon, TGeoVolumeAssembly, TGeoBoolNode (composite shape) and TGeoPatternFinder (divided volume). The threads book their navigators lazily, so it also exercises AddNavigator() against the navigator-map readers. Assisted-by: Claude Code (review, hardening and benchmarking) Supervised-by: Sandro Wenzel --- geom/test/CMakeLists.txt | 1 + geom/test/test_thread_navigation.cxx | 188 +++++++++++++++++++++++++++ 2 files changed, 189 insertions(+) create mode 100644 geom/test/test_thread_navigation.cxx diff --git a/geom/test/CMakeLists.txt b/geom/test/CMakeLists.txt index 1a935837a2c43..79131b795108c 100644 --- a/geom/test/CMakeLists.txt +++ b/geom/test/CMakeLists.txt @@ -8,4 +8,5 @@ configure_file(no_extrusion.gdml no_extrusion.gdml COPYONLY) ROOT_ADD_GTEST(geomTests test_material_units.cxx test_boolean_extrusion.cxx + test_thread_navigation.cxx LIBRARIES Geom) diff --git a/geom/test/test_thread_navigation.cxx b/geom/test/test_thread_navigation.cxx new file mode 100644 index 0000000000000..d01691526de78 --- /dev/null +++ b/geom/test/test_thread_navigation.cxx @@ -0,0 +1,188 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +/** + Navigating the same geometry from several threads must give exactly the same answer as + navigating it from one. This exercises every class that keeps per-thread scratch state: + TGeoXtru, TGeoPgon, TGeoVolumeAssembly, TGeoBoolNode (via a composite shape) and + TGeoPatternFinder (via a divided volume). + + Worth running under ThreadSanitizer: the threads book their navigators lazily, so this + also covers concurrent TGeoManager::AddNavigator() against navigator-map readers. +*/ + +namespace { + +/// Geometry containing one instance of each shape family that owns per-thread data. +TGeoManager *MakeGeometry() +{ + auto *geom = new TGeoManager("mt_nav_geom", "geometry for MT navigation test"); + + auto *matVac = new TGeoMaterial("Vacuum", 0, 0, 0); + auto *matAl = new TGeoMaterial("Al", 26.98, 13, 2.7); + auto *vac = new TGeoMedium("Vacuum", 1, matVac); + auto *alu = new TGeoMedium("Aluminium", 2, matAl); + + TGeoVolume *top = geom->MakeBox("TOP", vac, 100., 100., 100.); + geom->SetTopVolume(top); + + // --- TGeoXtru: convex, simple polygon extruded over two sections + Double_t xv[5] = {-10., -5., 5., 10., 0.}; + Double_t yv[5] = {-6., -10., -10., -6., 10.}; + auto *xtru = new TGeoXtru(2); + xtru->DefinePolygon(5, xv, yv); + xtru->DefineSection(0, -20., 0., 0., 1.); + xtru->DefineSection(1, 20., 0., 0., 1.); + top->AddNode(new TGeoVolume("XTRU", xtru, alu), 1, new TGeoTranslation(-45., 0., 0.)); + + // --- TGeoPgon + auto *pgon = new TGeoPgon("pgon", 0., 360., 8, 2); + pgon->DefineSection(0, -20., 5., 12.); + pgon->DefineSection(1, 20., 5., 12.); + top->AddNode(new TGeoVolume("PGON", pgon, alu), 1, new TGeoTranslation(45., 0., 0.)); + + // --- TGeoBoolNode, through a composite shape + new TGeoBBox("cbox", 12., 12., 12.); + new TGeoTube("ctub", 0., 6., 20.); + auto *comp = new TGeoCompositeShape("comp", "cbox - ctub"); + top->AddNode(new TGeoVolume("COMP", comp, alu), 1, new TGeoTranslation(0., 45., 0.)); + + // --- TGeoPatternFinder, through a divided volume + TGeoVolume *slab = geom->MakeBox("SLAB", alu, 20., 5., 5.); + slab->Divide("SLABDIV", 1, 10, -20., 4.); + top->AddNode(slab, 1, new TGeoTranslation(0., -45., 0.)); + + // --- TGeoVolumeAssembly + auto *assembly = new TGeoVolumeAssembly("ASSEMBLY"); + TGeoVolume *brick = geom->MakeBox("BRICK", alu, 3., 3., 3.); + for (int i = 0; i < 6; ++i) + assembly->AddNode(brick, i + 1, new TGeoTranslation(0., 0., -25. + 10. * i)); + top->AddNode(assembly, 1, new TGeoTranslation(0., 0., 0.)); + + geom->CloseGeometry(); + return geom; +} + +struct Ray { + Double_t point[3]; + Double_t dir[3]; +}; + +/// Deterministic fan of rays aimed at each shape, so every class with per-thread state is +/// actually traversed. Start points sit 40 cm from the target, well inside the world volume. +std::vector MakeRays() +{ + const Double_t targets[][3] = { + {-45., 0., 0.}, {-45., 0., 10.}, {-45., 3., -10.}, // XTRU + {45., 0., 0.}, {45., 0., 10.}, {45., 8., -10.}, // PGON + {0., 45., 0.}, {0., 45., 8.}, {8., 45., 0.}, // composite (box minus tube) + {0., -45., 0.}, {5., -45., 0.}, {-5., -45., 2.}, // divided slab + {0., 0., -25.}, {0., 0., -5.}, {0., 0., 15.}, // assembly bricks + }; + std::vector rays; + for (const auto &t : targets) { + for (int i = 0; i < 12; ++i) { + const double phi = 2. * M_PI * i / 12.; + const double theta = 0.4 + 0.15 * (i % 5); + const double d[3] = {std::sin(theta) * std::cos(phi), std::sin(theta) * std::sin(phi), std::cos(theta)}; + Ray ray; + for (int k = 0; k < 3; ++k) { + ray.point[k] = t[k] + 40. * d[k]; + ray.dir[k] = -d[k]; + } + rays.push_back(ray); + } + } + return rays; +} + +struct RayResult { + Int_t nsteps{0}; + Double_t pathlen{0.}; + Long64_t checksum{0}; // sequence of volumes traversed + + bool operator==(const RayResult &o) const + { + return nsteps == o.nsteps && checksum == o.checksum && std::abs(pathlen - o.pathlen) < 1e-9; + } +}; + +RayResult ShootRay(TGeoNavigator *nav, const Ray &ray) +{ + RayResult res; + nav->InitTrack(ray.point, ray.dir); + while (!nav->IsOutside() && res.nsteps < 500) { + TGeoNode *node = nav->GetCurrentNode(); + res.checksum = res.checksum * 31 + (node ? node->GetVolume()->GetNumber() : -1); + nav->FindNextBoundaryAndStep(1.e6, kFALSE); + res.pathlen += nav->GetStep(); + ++res.nsteps; + if (!nav->IsOnBoundary()) + break; + } + return res; +} + +} // namespace + +TEST(Geometry, MultiThreadedNavigationMatchesSerial) +{ + TGeoManager *geom = MakeGeometry(); + ASSERT_NE(geom, nullptr); + const std::vector rays = MakeRays(); + + // Reference: single-threaded, default navigator. + std::vector reference; + reference.reserve(rays.size()); + for (const Ray &ray : rays) + reference.push_back(ShootRay(geom->GetCurrentNavigator(), ray)); + + // Rays that miss every object legitimately cross a single boundary (the world). Guard against + // a vacuous comparison by requiring that a good share of them actually traverse the shapes. + const size_t nTraversing = + std::count_if(reference.begin(), reference.end(), [](const RayResult &r) { return r.nsteps > 2; }); + ASSERT_GT(nTraversing, reference.size() / 2); + + constexpr int kNThreads = 8; + geom->SetMaxThreads(kNThreads); + + std::vector> perThread(kNThreads); + std::vector threads; + threads.reserve(kNThreads); + for (int t = 0; t < kNThreads; ++t) { + threads.emplace_back([&, t] { + // Booked lazily and concurrently, exactly as a task-parallel workload would. + TGeoNavigator *nav = geom->AddNavigator(); + perThread[t].reserve(rays.size()); + for (const Ray &ray : rays) + perThread[t].push_back(ShootRay(nav, ray)); + }); + } + for (std::thread &th : threads) + th.join(); + + for (int t = 0; t < kNThreads; ++t) { + ASSERT_EQ(perThread[t].size(), reference.size()) << "thread " << t; + for (size_t i = 0; i < reference.size(); ++i) + EXPECT_TRUE(perThread[t][i] == reference[i]) << "thread " << t << ", ray " << i; + } + + delete geom; +}