[geom] Reduce thread false sharing and improve multithreaded TGeo navigation

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.

geom/geom/inc/TGeoBoolNode.h

@@ -14,7 +14,8 @@
 
 #include "TGeoShape.h"
 
-#include <mutex>
+#include <algorithm>
+#include <atomic>
 #include <vector>
 
 // forward declarations
@@ -23,21 +24,32 @@ class TGeoMatrix;
 class TGeoHMatrix;
 
 class TGeoBoolNode : public TObject {
+   static std::atomic<UInt_t> 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,
       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.
+   /// The vector retains its high-water size until the owning thread exits.
+   ThreadData_t &GetThreadData() const
+   {
+      thread_local std::vector<ThreadData_t> tdata;
+      if (tdata.size() <= fIndex)
+         tdata.resize(std::max<size_t>(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};          ///<! number of points on the mesh
    mutable Double_t *fPoints{nullptr}; ///<! array of mesh points
 
-   mutable std::vector<ThreadData_t *> 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
-   mutable Bool_t fMeshValid{kFALSE};               ///<! Flag for mesh cache validity
-                                                    // methods
+   mutable Bool_t fMeshValid{kFALSE}; ///<! Flag for mesh cache validity
+   // methods
    Bool_t MakeBranch(const char *expr, Bool_t left);
    void AssignPoints(Int_t npoints, Double_t *points);
 

geom/geom/inc/TGeoPatternFinder.h

@@ -14,7 +14,8 @@
 
 #include "TObject.h"
 
-#include <mutex>
+#include <algorithm>
+#include <atomic>
 #include <vector>
 
 #include "TGeoVolume.h"
@@ -23,36 +24,57 @@ class TGeoMatrix;
 
 /// base finder class for patterns. A pattern is specifying a division type
 class TGeoPatternFinder : public TObject {
+   static std::atomic<UInt_t> 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<Int_t> 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().
+   /// The vector retains its high-water size until the owning thread exits.
+   ThreadData_t &GetThreadData() const
+   {
+      thread_local std::vector<ThreadData_t> tdata;
+      if (tdata.size() <= fIndex)
+         tdata.resize(std::max<size_t>(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
    Int_t fNdivisions;   // number of divisions
    Int_t fDivIndex;     // index of first div. node
    TGeoVolume *fVolume; // volume to which applies
 
-   mutable std::vector<ThreadData_t *> 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 &);

geom/geom/inc/TGeoPgon.h

@@ -14,28 +14,51 @@
 
 #include "TGeoPcon.h"
 
+#include <algorithm>
+#include <atomic>
+#include <memory>
 #include <mutex>
 #include <vector>
 
 class TGeoPgon : public TGeoPcon {
+   static std::atomic<UInt_t> 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<Int_t> 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
-
-      ThreadData_t();
-      ~ThreadData_t();
+      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 &GetThreadData() const;
+
+   /// Per-thread non-owning cache of scratch buffers 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<ThreadData_t> tdata;
+      if (tdata.size() <= fIndex)
+         tdata.resize(std::max<size_t>(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<ThreadData_t *> 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)
+   mutable std::vector<std::unique_ptr<OwnedThreadData_t>> fOwnedData; ///<! Object-owned per-thread buffers
+   mutable std::mutex fOwnedDataMutex;                                 ///<! Protects cold allocation and cleanup
 
    // internal utility methods
    Int_t GetPhiCrossList(const Double_t *point, const Double_t *dir, Int_t istart, Double_t *sphi, Int_t *iphi,

geom/geom/inc/TGeoVolume.h

@@ -21,6 +21,8 @@
 #include "TObjArray.h"
 #include "TGeoMedium.h"
 #include "TGeoShape.h"
+#include <algorithm>
+#include <atomic>
 #include <mutex>
 #include <vector>
 
@@ -315,23 +317,33 @@ class TGeoVolumeMulti : public TGeoVolume {
 ////////////////////////////////////////////////////////////////////////////
 
 class TGeoVolumeAssembly : public TGeoVolume {
+   static std::atomic<UInt_t> 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; ///<! 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.
+   /// The vector retains its high-water size until the owning thread exits.
+   ThreadData_t &GetThreadData() const
+   {
+      thread_local std::vector<ThreadData_t> tdata;
+      if (tdata.size() <= fIndex)
+         tdata.resize(std::max<size_t>(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<ThreadData_t *> 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;
@@ -349,13 +361,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
 };

geom/geom/inc/TGeoXtru.h

@@ -14,28 +14,53 @@
 
 #include "TGeoBBox.h"
 
+#include <algorithm>
+#include <atomic>
+#include <memory>
 #include <mutex>
 #include <vector>
 
 class TGeoPolygon;
 
 class TGeoXtru : public TGeoBBox {
+   static std::atomic<UInt_t> 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<Int_t> fGeneration{0};           //! bumped whenever the per-thread state must be rebuilt
+   mutable std::atomic<Bool_t> 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<ThreadData_t> tdata;
+      if (tdata.size() <= fIndex)
+         tdata.resize(std::max<size_t>(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<ThreadData_t *> fThreadData; ///<! Navigation data per thread
-   mutable Int_t fThreadSize;                       ///<! size of thread-specific array
-   mutable std::mutex fMutex;                       ///<! mutex for thread data
+   mutable std::vector<std::unique_ptr<OwnedThreadData_t>> fOwnedData; ///<! Object-owned per-thread buffers
+   mutable std::mutex fOwnedDataMutex;                                 ///<! Protects cold allocation and cleanup
 
    TGeoXtru(const TGeoXtru &) = delete;
    TGeoXtru &operator=(const TGeoXtru &) = delete;