diff --git a/README/ReleaseNotes/v642/index.md b/README/ReleaseNotes/v642/index.md index e83165d9f748d..c7e3c336ab392 100644 --- a/README/ReleaseNotes/v642/index.md +++ b/README/ReleaseNotes/v642/index.md @@ -36,6 +36,7 @@ The following people have contributed to this new version: Devajith Valaparambil Sreeramaswamy, CERN/EP-SFT,\ Vassil Vassilev, Princeton,\ Sandro Wenzel, CERN/EP-ALICE,\ + Tristan Wenzel, ETHZ,\ ## Deprecation and Removal @@ -207,6 +208,15 @@ Such file can be loaded locally in any web browser or send as attachment in emai ## Geometry +### Improved multithreaded `TGeo` navigation + +Multithreaded `TGeo` navigation is now faster and more scalable, with improved thread-local state management that avoids +false sharing, releases temporary memory during geometry cleanup, and correctly supports concurrent navigation of +multiple geometries. + +For ALICE material-budget lookup-table generation on 28 cores, these changes reduced the runtime from 139 s to 72 s +and improved scaling from 12x to 23x. + The [TGeometry](https://root.cern/doc/master/classTGeometry.html) classes (Geant 3 shapes) have been moved out of Graf3D into their own library. To link to these classes, use the cmake target `TGeometry` (preferred), `root-config --libs`, or link with `-lTGeometry`. When ROOT is configured with `-Dgeom=Off`, these classes are now off as well. diff --git a/geom/geom/inc/TGeoBoolNode.h b/geom/geom/inc/TGeoBoolNode.h index bf3f5dbc005e0..025209fcd27c8 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,6 +24,9 @@ class TGeoMatrix; class TGeoHMatrix; class TGeoBoolNode : public TObject { + static std::atomic fgInstanceCount; //! source of monotonic per-object indices + UInt_t fIndex{fgInstanceCount++}; //! non-reused index of this node into the per-thread vector + public: enum EGeoBoolType { kGeoUnion, @@ -30,14 +34,22 @@ class TGeoBoolNode : public TObject { 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. + /// The vector retains its high-water size until the owning thread exits. + 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 {} + /// No-op: this node allocates its scratch state lazily for every calling thread. + void CreateThreadData(Int_t) {} private: TGeoBoolNode(const TGeoBoolNode &) = delete; @@ -51,11 +63,8 @@ class TGeoBoolNode : public TObject { mutable Int_t fNpoints{0}; /// fThreadData; /// +#include +#include #include #include "TGeoVolume.h" @@ -23,25 +24,50 @@ class TGeoMatrix; /// base finder class for patterns. A pattern is specifying a division type class TGeoPatternFinder : public TObject { + static std::atomic fgInstanceCount; //! source of monotonic per-object indices + UInt_t fIndex{fgInstanceCount++}; //! non-reused 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; /// 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: this finder allocates its scratch state lazily for every calling thread. + void CreateThreadData(Int_t) {} protected: - enum EGeoPatternFlags { kPatternReflected = BIT(14), kPatternSpacedOut = BIT(15) }; + void InitThreadSlot(ThreadData_t &td) const; + TGeoManager *GetOwnerManager() const { return fVolume ? fVolume->GetGeoManager() : nullptr; } + TGeoMatrix *GetOwnerIdentity() const; + void RegisterMatrix(TGeoMatrix *matrix) const; + + enum EGeoPatternFlags { + kPatternReflected = BIT(14), + kPatternSpacedOut = BIT(15) + }; Double_t fStep; // division step length Double_t fStart; // starting point on divided axis Double_t fEnd; // ending point @@ -49,10 +75,6 @@ class TGeoPatternFinder : public TObject { Int_t fDivIndex; // index of first div. node TGeoVolume *fVolume; // volume to which applies - mutable std::vector fThreadData; /// +#include +#include #include #include class TGeoPgon : public TGeoPcon { + static std::atomic fgInstanceCount; //! source of monotonic per-object indices + UInt_t fIndex{fgInstanceCount++}; //! non-reused 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; /// 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; + } + /// Release object-owned scratch buffers and invalidate the non-owning TLS slots. + /// Navigation using this shape must not be active when this method is called. void ClearThreadData() const override; - void CreateThreadData(Int_t nthreads) override; + /// No-op: this shape allocates scratch data lazily for every calling thread. + void CreateThreadData(Int_t) override {} protected: + struct OwnedThreadData_t; + void InitThreadSlot(ThreadData_t &td) const; + // data members - Int_t fNedges; // number of edges (at least one) - mutable std::vector fThreadData; ///> fOwnedData; /// +#include #include #include @@ -315,23 +317,33 @@ class TGeoVolumeMulti : public TGeoVolume { //////////////////////////////////////////////////////////////////////////// class TGeoVolumeAssembly : public TGeoVolume { + static std::atomic fgInstanceCount; //! source of monotonic per-object indices + UInt_t fIndex{fgInstanceCount++}; //! non-reused index of this assembly into the per-thread vector + public: struct ThreadData_t { - Int_t fCurrent; /// 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; /// +#include +#include #include #include class TGeoPolygon; class TGeoXtru : public TGeoBBox { + static std::atomic fgInstanceCount; //! source of monotonic per-object indices + UInt_t fIndex{fgInstanceCount++}; //! non-reused 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 - - ThreadData_t(); - ~ThreadData_t(); + 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 &GetThreadData() const; + + /// Per-thread non-owning cache of scratch state indexed by this shape. + /// Hot path: a TLS read plus an indexed load; the cold rebuild lives in InitThreadSlot(). + /// The vector retains its high-water size until the owning thread exits. + 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; + } + /// Release object-owned scratch buffers and invalidate the non-owning TLS slots. + /// Navigation using this shape must not be active when this method is called. void ClearThreadData() const override; - void CreateThreadData(Int_t nthreads) override; + /// No-op: this shape allocates scratch data lazily for every calling thread. + void CreateThreadData(Int_t) override {} protected: + struct OwnedThreadData_t; + 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) @@ -46,10 +71,8 @@ class TGeoXtru : public TGeoBBox { Double_t *fScale; //[fNz] array of scale factors (for each Z) 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; ///> fOwnedData; /// 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; -} - -//////////////////////////////////////////////////////////////////////////////// -/// 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); -} - +std::atomic TGeoBoolNode::fgInstanceCount{0}; //////////////////////////////////////////////////////////////////////////////// /// Set the selected branch. @@ -131,8 +63,6 @@ TGeoBoolNode::TGeoBoolNode() fRightMat = nullptr; fNpoints = 0; fPoints = nullptr; - fThreadSize = 0; - CreateThreadData(1); } //////////////////////////////////////////////////////////////////////////////// @@ -146,8 +76,6 @@ TGeoBoolNode::TGeoBoolNode(const char *expr1, const char *expr2) fRightMat = nullptr; fNpoints = 0; fPoints = nullptr; - fThreadSize = 0; - CreateThreadData(1); if (!MakeBranch(expr1, kTRUE)) { return; } @@ -166,8 +94,6 @@ TGeoBoolNode::TGeoBoolNode(TGeoShape *left, TGeoShape *right, TGeoMatrix *lmat, fLeftMat = lmat; fNpoints = 0; fPoints = nullptr; - fThreadSize = 0; - CreateThreadData(1); if (!fLeftMat) fLeftMat = gGeoIdentity; else @@ -195,7 +121,6 @@ TGeoBoolNode::~TGeoBoolNode() { if (fPoints) delete[] fPoints; - ClearThreadData(); } //////////////////////////////////////////////////////////////////////////////// diff --git a/geom/geom/src/TGeoManager.cxx b/geom/geom/src/TGeoManager.cxx index 4de5bc656e20d..f2a2d6c14f688 100644 --- a/geom/geom/src/TGeoManager.cxx +++ b/geom/geom/src/TGeoManager.cxx @@ -976,7 +976,9 @@ void TGeoManager::RemoveNavigator(const TGeoNavigator *nav) } //////////////////////////////////////////////////////////////////////////////// -/// Set maximum number of threads for navigation. +/// Enable multi-threaded navigation for at most `nthreads` worker threads. +/// The geometry must be closed and navigation must not be active when this method is called. +/// This enables ROOT thread safety and prepares the manager and geometry objects for concurrent navigation. void TGeoManager::SetMaxThreads(Int_t nthreads) { diff --git a/geom/geom/src/TGeoPatternFinder.cxx b/geom/geom/src/TGeoPatternFinder.cxx index 9616aa2b7ad2d..882d3edc10b7e 100644 --- a/geom/geom/src/TGeoPatternFinder.cxx +++ b/geom/geom/src/TGeoPatternFinder.cxx @@ -36,55 +36,61 @@ on different axis. Implemented patterns are: #include "TGeoManager.h" #include "TMath.h" +std::atomic TGeoPatternFinder::fgInstanceCount{0}; //////////////////////////////////////////////////////////////////////////////// -/// Constructor. +/// (Re)build the per-thread scratch state for this finder into the given slot. +/// Cold path: runs once per (thread, finder, generation). -TGeoPatternFinder::ThreadData_t::ThreadData_t() : fMatrix(nullptr), fCurrent(-1), fNextIndex(-1) {} - -//////////////////////////////////////////////////////////////////////////////// -/// Destructor. - -TGeoPatternFinder::ThreadData_t::~ThreadData_t() +void TGeoPatternFinder::InitThreadSlot(ThreadData_t &td) const { - // if (fMatrix != gGeoIdentity) delete fMatrix; + TGeoManager *manager = GetOwnerManager(); + if (!manager) { + Error("InitThreadSlot", "Pattern finder has no owning geometry manager"); + return; + } + 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 for + // its owning manager's lifetime and is deliberately reused. The manager releases it during + // destruction; replacing it here would retain another matrix on every ClearThreadData(). + td.fCurrent = -1; + td.fNextIndex = -1; + td.fInitGen = fGeneration.load(std::memory_order_acquire); } //////////////////////////////////////////////////////////////////////////////// +/// Return the identity matrix owned by this finder's geometry manager. -TGeoPatternFinder::ThreadData_t &TGeoPatternFinder::GetThreadData() const +TGeoMatrix *TGeoPatternFinder::GetOwnerIdentity() const { - Int_t tid = TGeoManager::ThreadId(); - return *fThreadData[tid]; + TGeoManager *manager = GetOwnerManager(); + return manager ? static_cast(manager->GetListOfMatrices()->At(0)) : nullptr; } //////////////////////////////////////////////////////////////////////////////// +/// Register a lazily-created pattern matrix with the manager owning this finder's volume. -void TGeoPatternFinder::ClearThreadData() const +void TGeoPatternFinder::RegisterMatrix(TGeoMatrix *matrix) const { - std::lock_guard guard(fMutex); - std::vector::iterator i = fThreadData.begin(); - while (i != fThreadData.end()) { - delete *i; - ++i; + TGeoManager *manager = GetOwnerManager(); + if (!manager || !matrix) + return; + if (!matrix->IsRegistered()) { + manager->RegisterMatrix(matrix); + matrix->SetBit(TGeoMatrix::kGeoRegistered); } - fThreadData.clear(); - fThreadSize = 0; -} - -//////////////////////////////////////////////////////////////////////////////// -/// Create thread data for n threads max. - -void TGeoPatternFinder::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; - fThreadData[tid]->fMatrix = CreateMatrix(); - } + if (matrix->IsCombi()) { + TGeoRotation *rotation = static_cast(matrix)->GetRotation(); + if (rotation && rotation->IsRotation()) + RegisterMatrix(rotation); } } @@ -99,7 +105,6 @@ TGeoPatternFinder::TGeoPatternFinder() fStart = 0; fEnd = 0; fVolume = nullptr; - fThreadSize = 0; } //////////////////////////////////////////////////////////////////////////////// @@ -113,7 +118,6 @@ TGeoPatternFinder::TGeoPatternFinder(TGeoVolume *vol, Int_t ndiv) fStep = 0; fStart = 0; fEnd = 0; - fThreadSize = 0; } //////////////////////////////////////////////////////////////////////////////// @@ -300,11 +304,11 @@ TGeoMatrix *TGeoPatternX::CreateMatrix() const { if (!IsReflected()) { TGeoMatrix *matrix = new TGeoTranslation(0., 0., 0.); - matrix->RegisterYourself(); + RegisterMatrix(matrix); return matrix; } TGeoCombiTrans *combi = new TGeoCombiTrans(); - combi->RegisterYourself(); + RegisterMatrix(combi); combi->ReflectZ(kTRUE); combi->ReflectZ(kFALSE); return combi; @@ -494,11 +498,11 @@ TGeoMatrix *TGeoPatternY::CreateMatrix() const { if (!IsReflected()) { TGeoMatrix *matrix = new TGeoTranslation(0., 0., 0.); - matrix->RegisterYourself(); + RegisterMatrix(matrix); return matrix; } TGeoCombiTrans *combi = new TGeoCombiTrans(); - combi->RegisterYourself(); + RegisterMatrix(combi); combi->ReflectZ(kTRUE); combi->ReflectZ(kFALSE); return combi; @@ -684,11 +688,11 @@ TGeoMatrix *TGeoPatternZ::CreateMatrix() const { if (!IsReflected()) { TGeoMatrix *matrix = new TGeoTranslation(0., 0., 0.); - matrix->RegisterYourself(); + RegisterMatrix(matrix); return matrix; } TGeoCombiTrans *combi = new TGeoCombiTrans(); - combi->RegisterYourself(); + RegisterMatrix(combi); combi->ReflectZ(kTRUE); combi->ReflectZ(kFALSE); return combi; @@ -947,11 +951,11 @@ TGeoMatrix *TGeoPatternParaX::CreateMatrix() const { if (!IsReflected()) { TGeoMatrix *matrix = new TGeoTranslation(0., 0., 0.); - matrix->RegisterYourself(); + RegisterMatrix(matrix); return matrix; } TGeoCombiTrans *combi = new TGeoCombiTrans(); - combi->RegisterYourself(); + RegisterMatrix(combi); combi->ReflectZ(kTRUE); combi->ReflectZ(kFALSE); return combi; @@ -1125,11 +1129,11 @@ TGeoMatrix *TGeoPatternParaY::CreateMatrix() const { if (!IsReflected()) { TGeoMatrix *matrix = new TGeoTranslation(0., 0., 0.); - matrix->RegisterYourself(); + RegisterMatrix(matrix); return matrix; } TGeoCombiTrans *combi = new TGeoCombiTrans(); - combi->RegisterYourself(); + RegisterMatrix(combi); combi->ReflectZ(kTRUE); combi->ReflectZ(kFALSE); return combi; @@ -1307,11 +1311,11 @@ TGeoMatrix *TGeoPatternParaZ::CreateMatrix() const { if (!IsReflected()) { TGeoMatrix *matrix = new TGeoTranslation(0., 0., 0.); - matrix->RegisterYourself(); + RegisterMatrix(matrix); return matrix; } TGeoCombiTrans *combi = new TGeoCombiTrans(); - combi->RegisterYourself(); + RegisterMatrix(combi); combi->ReflectZ(kTRUE); combi->ReflectZ(kFALSE); return combi; @@ -1496,11 +1500,11 @@ TGeoMatrix *TGeoPatternTrapZ::CreateMatrix() const { if (!IsReflected()) { TGeoMatrix *matrix = new TGeoTranslation(0., 0., 0.); - matrix->RegisterYourself(); + RegisterMatrix(matrix); return matrix; } TGeoCombiTrans *combi = new TGeoCombiTrans(); - combi->RegisterYourself(); + RegisterMatrix(combi); combi->ReflectZ(kTRUE); combi->ReflectZ(kFALSE); return combi; @@ -1660,7 +1664,7 @@ void TGeoPatternCylR::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= TGeoMatrix *TGeoPatternCylR::CreateMatrix() const { - return gGeoIdentity; + return GetOwnerIdentity(); } //////////////////////////////////////////////////////////////////////////////// @@ -1856,11 +1860,11 @@ TGeoMatrix *TGeoPatternCylPhi::CreateMatrix() const { if (!IsReflected()) { TGeoRotation *matrix = new TGeoRotation(); - matrix->RegisterYourself(); + RegisterMatrix(matrix); return matrix; } TGeoRotation *rot = new TGeoRotation(); - rot->RegisterYourself(); + RegisterMatrix(rot); rot->ReflectZ(kTRUE); rot->ReflectZ(kFALSE); return rot; @@ -1978,7 +1982,7 @@ void TGeoPatternSphR::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= TGeoMatrix *TGeoPatternSphR::CreateMatrix() const { - return gGeoIdentity; + return GetOwnerIdentity(); } //////////////////////////////////////////////////////////////////////////////// @@ -2091,7 +2095,7 @@ void TGeoPatternSphTheta::SavePrimitive(std::ostream &out, Option_t * /*option*/ TGeoMatrix *TGeoPatternSphTheta::CreateMatrix() const { - return gGeoIdentity; + return GetOwnerIdentity(); } //////////////////////////////////////////////////////////////////////////////// @@ -2270,11 +2274,11 @@ TGeoMatrix *TGeoPatternSphPhi::CreateMatrix() const { if (!IsReflected()) { TGeoRotation *matrix = new TGeoRotation(); - matrix->RegisterYourself(); + RegisterMatrix(matrix); return matrix; } TGeoRotation *rot = new TGeoRotation(); - rot->RegisterYourself(); + RegisterMatrix(rot); rot->ReflectZ(kTRUE); rot->ReflectZ(kFALSE); return rot; @@ -2371,7 +2375,7 @@ TGeoNode *TGeoPatternHoneycomb::FindNode(Double_t * /*point*/, const Double_t * TGeoMatrix *TGeoPatternHoneycomb::CreateMatrix() const { - return gGeoIdentity; + return GetOwnerIdentity(); } //////////////////////////////////////////////////////////////////////////////// diff --git a/geom/geom/src/TGeoPgon.cxx b/geom/geom/src/TGeoPgon.cxx index 5444113c2ed40..bdfe0fe0ab05b 100644 --- a/geom/geom/src/TGeoPgon.cxx +++ b/geom/geom/src/TGeoPgon.cxx @@ -66,60 +66,40 @@ polygons, between `phi1` and `phi1+dphi.` #include "TBuffer3DTypes.h" #include "TMath.h" +std::atomic TGeoPgon::fgInstanceCount{0}; -//////////////////////////////////////////////////////////////////////////////// -/// Constructor. +struct TGeoPgon::OwnedThreadData_t { + std::unique_ptr fIntBuffer; + std::unique_ptr fDblBuffer; -TGeoPgon::ThreadData_t::ThreadData_t() : fIntBuffer(nullptr), fDblBuffer(nullptr) {} + explicit OwnedThreadData_t(std::size_t size) : fIntBuffer(new Int_t[size]), fDblBuffer(new Double_t[size]) {} +}; //////////////////////////////////////////////////////////////////////////////// -/// Destructor. +/// (Re)build the per-thread scratch buffers for this shape into the given slot. +/// Cold path: runs once per (thread, shape, generation). -TGeoPgon::ThreadData_t::~ThreadData_t() +void TGeoPgon::InitThreadSlot(ThreadData_t &td) const { - delete[] fIntBuffer; - delete[] fDblBuffer; -} - -//////////////////////////////////////////////////////////////////////////////// - -TGeoPgon::ThreadData_t &TGeoPgon::GetThreadData() const -{ - Int_t tid = TGeoManager::ThreadId(); - return *fThreadData[tid]; + auto data = std::make_unique(fNedges + 10); + Int_t *intBuffer = data->fIntBuffer.get(); + Double_t *dblBuffer = data->fDblBuffer.get(); + + std::lock_guard guard(fOwnedDataMutex); + fOwnedData.push_back(std::move(data)); + td.fIntBuffer = intBuffer; + td.fDblBuffer = dblBuffer; + td.fInitGen = fGeneration.load(std::memory_order_acquire); } //////////////////////////////////////////////////////////////////////////////// +/// Release the large scratch buffers. Navigation using this shape must not be active. void TGeoPgon::ClearThreadData() const { - std::lock_guard guard(fMutex); - std::vector::iterator i = fThreadData.begin(); - while (i != fThreadData.end()) { - delete *i; - ++i; - } - fThreadData.clear(); - fThreadSize = 0; -} - -//////////////////////////////////////////////////////////////////////////////// -/// Create thread data for n threads max. - -void TGeoPgon::CreateThreadData(Int_t nthreads) -{ - 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]; - } - } + std::lock_guard guard(fOwnedDataMutex); + fOwnedData.clear(); + fGeneration.fetch_add(1, std::memory_order_release); } //////////////////////////////////////////////////////////////////////////////// @@ -129,7 +109,6 @@ TGeoPgon::TGeoPgon() { SetShapeBit(TGeoShape::kGeoPgon); fNedges = 0; - fThreadSize = 0; } //////////////////////////////////////////////////////////////////////////////// @@ -139,8 +118,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); } //////////////////////////////////////////////////////////////////////////////// @@ -151,8 +128,6 @@ TGeoPgon::TGeoPgon(const char *name, Double_t phi, Double_t dphi, Int_t nedges, { SetShapeBit(TGeoShape::kGeoPgon); fNedges = nedges; - fThreadSize = 0; - CreateThreadData(1); } //////////////////////////////////////////////////////////////////////////////// @@ -171,8 +146,6 @@ TGeoPgon::TGeoPgon(Double_t *param) : TGeoPcon("") SetShapeBit(TGeoShape::kGeoPgon); SetDimensions(param); ComputeBBox(); - fThreadSize = 0; - CreateThreadData(1); } //////////////////////////////////////////////////////////////////////////////// @@ -466,8 +439,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; @@ -956,7 +927,7 @@ Bool_t TGeoPgon::SliceCrossingIn(const Double_t *point, const Double_t *dir, Int } ipl += incseg; } // end loop Z - } // end loop phi + } // end loop phi snext = TGeoShape::Big(); return kFALSE; } @@ -1253,8 +1224,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 d640852acaa50..34ca4fbaaf93d 100644 --- a/geom/geom/src/TGeoVolume.cxx +++ b/geom/geom/src/TGeoVolume.cxx @@ -2973,91 +2973,11 @@ void TGeoVolumeMulti::SetVisibility(Bool_t vis) } } -//////////////////////////////////////////////////////////////////////////////// -/// 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]; -} - -//////////////////////////////////////////////////////////////////////////////// - -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; -} - +std::atomic TGeoVolumeAssembly::fgInstanceCount{0}; //////////////////////////////////////////////////////////////////////////////// /// 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 @@ -3070,8 +2990,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 36f16eb751fdb..4233f96a6d939 100644 --- a/geom/geom/src/TGeoXtru.cxx +++ b/geom/geom/src/TGeoXtru.cxx @@ -102,69 +102,55 @@ Double_t y0, Double_t scale); #include "TGeoVolume.h" #include "TGeoPolygon.h" -//////////////////////////////////////////////////////////////////////////////// -/// Constructor. +std::atomic TGeoXtru::fgInstanceCount{0}; -TGeoXtru::ThreadData_t::ThreadData_t() : fSeg(0), fIz(0), fXc(nullptr), fYc(nullptr), fPoly(nullptr) {} +struct TGeoXtru::OwnedThreadData_t { + std::unique_ptr fXc; + std::unique_ptr fYc; + std::unique_ptr fPoly; -//////////////////////////////////////////////////////////////////////////////// -/// Destructor. - -TGeoXtru::ThreadData_t::~ThreadData_t() -{ - delete[] fXc; - delete[] fYc; - delete fPoly; -} + explicit OwnedThreadData_t(std::size_t size) : fXc(new Double_t[size]), fYc(new Double_t[size]) {} +}; //////////////////////////////////////////////////////////////////////////////// +/// (Re)build the per-thread scratch state for this shape into the given slot. +/// Cold path: runs once per (thread, shape, generation). -TGeoXtru::ThreadData_t &TGeoXtru::GetThreadData() const +void TGeoXtru::InitThreadSlot(ThreadData_t &td) const { - if (!fThreadSize) - ((TGeoXtru *)this)->CreateThreadData(1); - Int_t tid = TGeoManager::ThreadId(); - return *fThreadData[tid]; + auto data = std::make_unique(fNvert); + memcpy(data->fXc.get(), fX, fNvert * sizeof(Double_t)); + memcpy(data->fYc.get(), fY, fNvert * sizeof(Double_t)); + data->fPoly = std::make_unique(fNvert); + data->fPoly->SetXY(data->fXc.get(), data->fYc.get()); + data->fPoly->FinishPolygon(); + Double_t *xc = data->fXc.get(); + Double_t *yc = data->fYc.get(); + TGeoPolygon *poly = data->fPoly.get(); + + std::lock_guard guard(fOwnedDataMutex); + fOwnedData.push_back(std::move(data)); + td.fSeg = 0; + td.fIz = 0; + td.fXc = xc; + td.fYc = yc; + td.fPoly = poly; + td.fInitGen = fGeneration.load(std::memory_order_acquire); + // 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()); } //////////////////////////////////////////////////////////////////////////////// +/// Release the large scratch buffers. Navigation using this shape must not be active. void TGeoXtru::ClearThreadData() const { - std::lock_guard guard(fMutex); - std::vector::iterator i = fThreadData.begin(); - while (i != fThreadData.end()) { - delete *i; - ++i; - } - fThreadData.clear(); - fThreadSize = 0; -} - -//////////////////////////////////////////////////////////////////////////////// -/// Create thread data for n threads max. - -void TGeoXtru::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; - 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()); - } - } - } + std::lock_guard guard(fOwnedDataMutex); + fOwnedData.clear(); + fGeneration.fetch_add(1, std::memory_order_release); + fIllegalChecked.store(kFALSE, std::memory_order_relaxed); } //////////////////////////////////////////////////////////////////////////////// @@ -195,9 +181,7 @@ TGeoXtru::TGeoXtru() fZ(nullptr), fScale(nullptr), fX0(nullptr), - fY0(nullptr), - fThreadData(0), - fThreadSize(0) + fY0(nullptr) { SetShapeBit(TGeoShape::kGeoXtru); } @@ -215,9 +199,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) { @@ -251,9 +233,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); diff --git a/geom/test/CMakeLists.txt b/geom/test/CMakeLists.txt index c52ed9610d743..11e4c69f8a083 100644 --- a/geom/test/CMakeLists.txt +++ b/geom/test/CMakeLists.txt @@ -25,6 +25,10 @@ ROOT_ADD_GTEST(overlap_navigation test_overlap_navigation.cxx LIBRARIES Geom) +ROOT_ADD_GTEST(thread_navigation + test_thread_navigation.cxx + LIBRARIES Geom) + if(imt) ROOT_ADD_GTEST(manager_lifetime test_manager_lifetime.cxx diff --git a/geom/test/test_thread_navigation.cxx b/geom/test/test_thread_navigation.cxx new file mode 100644 index 0000000000000..59e5c6e288549 --- /dev/null +++ b/geom/test/test_thread_navigation.cxx @@ -0,0 +1,326 @@ +#include + +#include +#include +#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 GeometryWithFinders { + TGeoManager *manager; + TGeoPatternFinder *linearFinder; + TGeoPatternFinder *radialFinder; +}; + +/// Build two divided volumes without touching their lazily-created matrices. +GeometryWithFinders MakeDividedGeometry(const char *name) +{ + auto *manager = new TGeoManager(name, name); + auto *material = new TGeoMaterial("vacuum", 0., 0., 0.); + auto *medium = new TGeoMedium("vacuum", 1, material); + auto *top = manager->MakeBox("top", medium, 100., 100., 100.); + manager->SetTopVolume(top); + + auto *slab = manager->MakeBox("slab", medium, 20., 5., 5.); + slab->Divide("linear_slice", 1, 10, -20., 4.); + top->AddNode(slab, 1); + + auto *tube = manager->MakeTube("tube", medium, 1., 20., 10.); + tube->Divide("radial_slice", 1, 4, 1., 19. / 4.); + top->AddNode(tube, 1, new TGeoTranslation(0., 40., 0.)); + + manager->CloseGeometry(); + return {manager, slab->GetFinder(), tube->GetFinder()}; +} + +void MakeCurrent(TGeoManager *manager) +{ + gGeoManager = manager; + gGeoIdentity = static_cast(manager->GetListOfMatrices()->At(0)); +} + +class InspectablePgon : public TGeoPgon { +public: + using TGeoPgon::TGeoPgon; + + std::size_t GetOwnedThreadDataCount() const + { + std::lock_guard guard(fOwnedDataMutex); + return fOwnedData.size(); + } +}; + +class InspectableXtru : public TGeoXtru { +public: + using TGeoXtru::TGeoXtru; + + std::size_t GetOwnedThreadDataCount() const + { + std::lock_guard guard(fOwnedDataMutex); + return fOwnedData.size(); + } +}; + +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; +} + +TEST(Geometry, PatternMatricesBelongToOwningManager) +{ + auto geometryA = MakeDividedGeometry("pattern_owner_A"); + + // Keep A alive while constructing B. This is how applications that cache several + // managers switch away from the current geometry before creating another one. + gGeoManager = nullptr; + gGeoIdentity = nullptr; + auto geometryB = MakeDividedGeometry("pattern_owner_B"); + + // First-touch A while B is current. Matrix ownership must follow the divided volume, + // not the ambient globals. + TGeoMatrix *matrix = geometryA.linearFinder->GetMatrix(); + ASSERT_NE(matrix, nullptr); + EXPECT_GE(geometryA.manager->GetListOfMatrices()->IndexOf(matrix), 0); + EXPECT_LT(geometryB.manager->GetListOfMatrices()->IndexOf(matrix), 0); + + TGeoMatrix *identity = geometryA.radialFinder->GetMatrix(); + EXPECT_EQ(identity, geometryA.manager->GetListOfMatrices()->At(0)); + EXPECT_NE(identity, geometryB.manager->GetListOfMatrices()->At(0)); + + delete geometryB.manager; + MakeCurrent(geometryA.manager); + + // Under ASan this dereference also catches a matrix that was wrongly owned and deleted by B. + geometryA.linearFinder->cd(0); + EXPECT_EQ(geometryA.linearFinder->GetMatrix(), matrix); + EXPECT_TRUE(geometryA.radialFinder->GetMatrix()->IsIdentity()); + + delete geometryA.manager; +} + +TEST(Geometry, ShapeScratchDataReleasedOnClear) +{ + InspectablePgon pgon(0., 360., 64, 2); + pgon.DefineSection(0, -10., 1., 5.); + pgon.DefineSection(1, 10., 1., 5.); + + InspectableXtru xtru(2); + Double_t x[] = {-5., 5., 5., -5.}; + Double_t y[] = {-5., -5., 5., 5.}; + xtru.DefinePolygon(4, x, y); + xtru.DefineSection(0, -10.); + xtru.DefineSection(1, 10.); + + // DefineSection computes the Xtru bounding box using the main-thread slot. + pgon.ClearThreadData(); + xtru.ClearThreadData(); + + constexpr int kNThreads = 8; + std::atomic valid{true}; + std::vector threads; + threads.reserve(kNThreads); + for (int i = 0; i < kNThreads; ++i) { + threads.emplace_back([&] { + auto &pgonData = pgon.GetThreadData(); + auto &xtruData = xtru.GetThreadData(); + if (!pgonData.fIntBuffer || !pgonData.fDblBuffer || !xtruData.fXc || !xtruData.fYc || !xtruData.fPoly) + valid.store(false, std::memory_order_relaxed); + }); + } + for (auto &thread : threads) + thread.join(); + + ASSERT_TRUE(valid.load(std::memory_order_relaxed)); + EXPECT_EQ(pgon.GetOwnedThreadDataCount(), kNThreads); + EXPECT_EQ(xtru.GetOwnedThreadDataCount(), kNThreads); + + pgon.ClearThreadData(); + xtru.ClearThreadData(); + EXPECT_EQ(pgon.GetOwnedThreadDataCount(), 0u); + EXPECT_EQ(xtru.GetOwnedThreadDataCount(), 0u); + + // The main thread's stale non-owning slots must rebuild on the next access. + EXPECT_NE(pgon.GetThreadData().fIntBuffer, nullptr); + EXPECT_NE(xtru.GetThreadData().fPoly, nullptr); + EXPECT_EQ(pgon.GetOwnedThreadDataCount(), 1u); + EXPECT_EQ(xtru.GetOwnedThreadDataCount(), 1u); +}