Parallel overlap checking with staged workflow

geom/geom/CMakeLists.txt

@@ -124,6 +124,8 @@ ROOT_STANDARD_LIBRARY_PACKAGE(Geom
     RIO
     MathCore
     Hist
+    Physics
+    Imt
 )
 
 # GCC has bugs with -O3 or -Ofast that break Geom

geom/geom/inc/TGeoArb8.h

@@ -72,6 +72,7 @@ class TGeoArb8 : public TGeoBBox {
    static void GetPlaneNormal(Double_t *p1, Double_t *p2, Double_t *p3, Double_t *norm);
    Double_t *GetVertices() { return &fXY[0][0]; }
    Double_t GetTwist(Int_t iseg) const;
+   Bool_t IsConvex() const override { return !IsTwisted(); }
    Bool_t IsCylType() const override { return kFALSE; }
    static Bool_t IsSamePoint(const Double_t *p1, const Double_t *p2)
    {
@@ -140,6 +141,7 @@ class TGeoTrap : public TGeoArb8 {
    Double_t GetTl2() const { return fTl2; }
    Double_t GetAlpha2() const { return fAlpha2; }
    TGeoShape *GetMakeRuntimeShape(TGeoShape *mother, TGeoMatrix *mat) const override;
+   Bool_t IsConvex() const override { return kTRUE; }
    void SetDimensions(Double_t *param) override;
    Double_t Safety(const Double_t *point, Bool_t in = kTRUE) const override;
    void Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const override;
@@ -171,6 +173,7 @@ class TGeoGtra : public TGeoTrap {
                           Double_t *step) const override;
    TGeoShape *GetMakeRuntimeShape(TGeoShape *mother, TGeoMatrix *mat) const override;
    Double_t GetTwistAngle() const { return fTwistAngle; }
+   Bool_t IsConvex() const final { return kFALSE; }
    Double_t Safety(const Double_t *point, Bool_t in = kTRUE) const override;
    void Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const override;
    void SetDimensions(Double_t *param) override;

geom/geom/inc/TGeoBBox.h

@@ -13,6 +13,7 @@
 #define ROOT_TGeoBBox
 
 #include "TGeoShape.h"
+#include "TVector3.h"
 
 class TGeoBBox : public TGeoShape {
 protected:
@@ -36,7 +37,7 @@ class TGeoBBox : public TGeoShape {
    // destructor
    ~TGeoBBox() override;
    // methods
-   
+
    Double_t Capacity() const override;
    void ComputeBBox() override;
    void ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm) const override;
@@ -65,6 +66,7 @@ class TGeoBBox : public TGeoShape {
    void GetBoundingCylinder(Double_t *param) const override;
    const TBuffer3D &GetBuffer3D(Int_t reqSections, Bool_t localFrame) const override;
    Int_t GetByteCount() const override { return 36; }
+   inline TVector3 GetDimensions() const { return TVector3(fDX, fDY, fDZ); }
    virtual Double_t GetFacetArea(Int_t index = 0) const;
    virtual Bool_t GetPointsOnFacet(Int_t index, Int_t npoints, Double_t *array) const;
    Bool_t GetPointsOnSegments(Int_t npoints, Double_t *array) const override;
@@ -77,11 +79,17 @@ class TGeoBBox : public TGeoShape {
    virtual Double_t GetDY() const { return fDY; }
    virtual Double_t GetDZ() const { return fDZ; }
    virtual const Double_t *GetOrigin() const { return fOrigin; }
+   TVector3 GetWorldCenter(const TGeoMatrix *m) const;
    void InspectShape() const override;
+   Bool_t IsConvex() const override { return kTRUE; }
    Bool_t IsCylType() const override { return kFALSE; }
    Bool_t IsValidBox() const override { return ((fDX < 0) || (fDY < 0) || (fDZ < 0)) ? kFALSE : kTRUE; }
    virtual Bool_t IsNullBox() const { return ((fDX < 1.E-16) && (fDY < 1.E-16) && (fDZ < 1.E-16)) ? kTRUE : kFALSE; }
+   static Bool_t IsSeparatingAxis(const TVector3 &L, const TVector3 &D, const TVector3 &Ax, const TVector3 &Ay,
+                                  const TVector3 &Az, const TVector3 &Bx, const TVector3 &By, const TVector3 &Bz,
+                                  const TVector3 &dA, const TVector3 &dB, Double_t tol);
    TBuffer3D *MakeBuffer3D() const override;
+   static bool MayIntersect(const TGeoBBox *boxA, const TGeoMatrix *mA, const TGeoBBox *boxB, const TGeoMatrix *mB);
    Double_t Safety(const Double_t *point, Bool_t in = kTRUE) const override;
    void Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const override;
    void SavePrimitive(std::ostream &out, Option_t *option = "") override;

geom/geom/inc/TGeoBoolNode.h

@@ -24,7 +24,11 @@ class TGeoHMatrix;
 
 class TGeoBoolNode : public TObject {
 public:
-   enum EGeoBoolType { kGeoUnion, kGeoIntersection, kGeoSubtraction };
+   enum EGeoBoolType {
+      kGeoUnion,
+      kGeoIntersection,
+      kGeoSubtraction
+   };
    struct ThreadData_t {
       Int_t fSelected; // ! selected branch
 
@@ -40,16 +44,17 @@ class TGeoBoolNode : public TObject {
    TGeoBoolNode &operator=(const TGeoBoolNode &) = delete;
 
 protected:
-   TGeoShape *fLeft{nullptr};      // shape on the left branch
-   TGeoShape *fRight{nullptr};     // shape on the right branch
-   TGeoMatrix *fLeftMat{nullptr};  // transformation that applies to the left branch
-   TGeoMatrix *fRightMat{nullptr}; // transformation that applies to the right branch
-   Int_t fNpoints{0};              //! number of points on the mesh
-   Double_t *fPoints{nullptr};     //! array of mesh points
+   TGeoShape *fLeft{nullptr};          // shape on the left branch
+   TGeoShape *fRight{nullptr};         // shape on the right branch
+   TGeoMatrix *fLeftMat{nullptr};      // transformation that applies to the left branch
+   TGeoMatrix *fRightMat{nullptr};     // transformation that applies to the right branch
+   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
    Bool_t MakeBranch(const char *expr, Bool_t left);
    void AssignPoints(Int_t npoints, Double_t *points);
@@ -78,6 +83,7 @@ class TGeoBoolNode : public TObject {
    TGeoShape *GetLeftShape() const { return fLeft; }
    TGeoShape *GetRightShape() const { return fRight; }
    TGeoShape::EInside Inside(const Double_t *point) const;
+   void InvalidateMeshCaches();
    virtual TGeoBoolNode *MakeClone() const = 0;
    void Paint(Option_t *option) override;
    void RegisterMatrices();

geom/geom/inc/TGeoCompositeShape.h

@@ -64,7 +64,9 @@ class TGeoCompositeShape : public TGeoBBox {
    Bool_t GetPointsOnSegments(Int_t /*npoints*/, Double_t * /*array*/) const override { return kFALSE; }
    void InspectShape() const override;
    Bool_t IsComposite() const override { return kTRUE; }
+   Bool_t IsConvex() const final { return kFALSE; }
    Bool_t IsCylType() const override { return kFALSE; }
+   void InvalidateMeshCaches();
    void MakeNode(const char *expression);
    virtual Bool_t PaintComposite(Option_t *option = "") const;
    void RegisterYourself();

geom/geom/inc/TGeoCone.h

@@ -77,10 +77,10 @@ class TGeoCone : public TGeoBBox {
    virtual Double_t GetRmax1() const { return fRmax1; }
    virtual Double_t GetRmin2() const { return fRmin2; }
    virtual Double_t GetRmax2() const { return fRmax2; }
-
+   inline Bool_t HasRmin() const { return (fRmin1 > 0. || fRmin2 > 0.); }
    void InspectShape() const override;
+   Bool_t IsConvex() const override { return !HasRmin(); }
    Bool_t IsCylType() const override { return kTRUE; }
-   TBuffer3D *MakeBuffer3D() const override;
    Double_t Safety(const Double_t *point, Bool_t in = kTRUE) const override;
    void Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const override;
    static Double_t SafetyS(const Double_t *point, Bool_t in, Double_t dz, Double_t rmin1, Double_t rmax1,
@@ -165,8 +165,9 @@ class TGeoConeSeg : public TGeoCone {
    Bool_t GetPointsOnSegments(Int_t npoints, Double_t *array) const override;
    Double_t GetPhi1() const { return fPhi1; }
    Double_t GetPhi2() const { return fPhi2; }
+   inline Bool_t HasDphi() const { return (fPhi2 - fPhi1) < 360.; }
    void InspectShape() const override;
-   TBuffer3D *MakeBuffer3D() const override;
+   Bool_t IsConvex() const final { return !HasRmin() && !HasDphi(); }
    Double_t Safety(const Double_t *point, Bool_t in = kTRUE) const override;
    void Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const override;
    static Double_t SafetyS(const Double_t *point, Bool_t in, Double_t dz, Double_t rmin1, Double_t rmax1,

geom/geom/inc/TGeoEltu.h

@@ -53,6 +53,7 @@ class TGeoEltu : public TGeoTube {
    Int_t GetNmeshVertices() const override;
    Bool_t GetPointsOnSegments(Int_t /*npoints*/, Double_t * /*array*/) const override { return kFALSE; }
    void InspectShape() const override;
+   Bool_t IsConvex() const final { return kTRUE; }
    Bool_t IsCylType() const override { return kTRUE; }
    Double_t Safety(const Double_t *point, Bool_t in = kTRUE) const override;
    void Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const override;

geom/geom/inc/TGeoHype.h

@@ -75,6 +75,7 @@ class TGeoHype : public TGeoTube {
    Double_t RadiusHypeSq(Double_t z, Bool_t inner) const;
    Double_t ZHypeSq(Double_t r, Bool_t inner) const;
    void InspectShape() const override;
+   Bool_t IsConvex() const final { return kFALSE; }
    Bool_t IsCylType() const override { return kTRUE; }
    TBuffer3D *MakeBuffer3D() const override;
    // virtual void          Paint(Option_t *option);

geom/geom/inc/TGeoManager.h

@@ -44,7 +44,10 @@ class TGeoBorderSurface;
 
 class TGeoManager : public TNamed {
 public:
-   enum EDefaultUnits { kG4Units = 0, kRootUnits = 1 };
+   enum EDefaultUnits {
+      kG4Units = 0,
+      kRootUnits = 1
+   };
 
 protected:
    static std::mutex fgMutex;           //! mutex for navigator booking in MT mode
@@ -221,12 +224,14 @@ class TGeoManager : public TNamed {
    Double_t GetVisDensity() const { return fVisDensity; }
    Int_t GetVisLevel() const;
    Int_t GetVisOption() const;
+   void InvalidateMeshCaches();
    Bool_t IsInPhiRange() const;
    Bool_t IsDrawingExtra() const { return fDrawExtra; }
    Bool_t IsNodeSelectable() const { return fIsNodeSelectable; }
    Bool_t IsVisLeaves() const { return fVisOption == 1; }
    void ModifiedPad() const;
-   void OptimizeVoxels(const char *filename = "tgeovox.C");                                     // *MENU*
+   void OptimizeVoxels(const char *filename = "tgeovox.C"); // *MENU*
+   void RebuildVoxels();
    void SetClipping(Bool_t flag = kTRUE) { SetClippingShape(flag ? fClippingShape : nullptr); } // *MENU*
    void SetClippingShape(TGeoShape *clip);
    void SetExplodedView(Int_t iopt = 0); // *MENU*
@@ -255,8 +260,10 @@ class TGeoManager : public TNamed {
    void CheckGeometryFull(Int_t ntracks = 1000000, Double_t vx = 0., Double_t vy = 0., Double_t vz = 0.,
                           Option_t *option = "ob"); // *MENU*
    void CheckGeometry(Option_t *option = "");
-   void CheckOverlaps(Double_t ovlp = 0.1, Option_t *option = "");                         // *MENU*
-   void CheckPoint(Double_t x = 0, Double_t y = 0, Double_t z = 0, Option_t *option = "", Double_t safety = 0.); // *MENU*
+   void CheckOverlaps(Double_t ovlp = 0.1, Option_t *option = ""); // *MENU*
+   void CheckOverlapsBySampling(Double_t ovlp, Int_t npoints);
+   void CheckPoint(Double_t x = 0, Double_t y = 0, Double_t z = 0, Option_t *option = "",
+                   Double_t safety = 0.); // *MENU*
    void CheckShape(TGeoShape *shape, Int_t testNo, Int_t nsamples, Option_t *option);
    void ConvertReflections();
    void DrawCurrentPoint(Int_t color = 2); // *MENU*

geom/geom/inc/TGeoMatrix.h

@@ -19,6 +19,7 @@
  *************************************************************************/
 
 #include "TNamed.h"
+#include "TVector3.h"
 
 //--- globals
 const Double_t kNullVector[3] = {0.0, 0.0, 0.0};
@@ -73,7 +74,7 @@ class TGeoMatrix : public TNamed {
    Bool_t IsRotAboutZ() const;
    void GetHomogenousMatrix(Double_t *hmat) const;
    const char *GetPointerName() const;
-
+   void GetWorldAxes(TVector3 &ax, TVector3 &ay, TVector3 &az) const;
    virtual Int_t GetByteCount() const;
    virtual const Double_t *GetTranslation() const = 0;
    virtual const Double_t *GetRotationMatrix() const = 0;
@@ -463,7 +464,7 @@ class TGeoHMatrix : public TGeoMatrix {
 
 public:
    TGeoHMatrix();
-   TGeoHMatrix(const TGeoHMatrix &other) : TGeoHMatrix((TGeoMatrix &)other) {}
+   TGeoHMatrix(const TGeoHMatrix &other);
    TGeoHMatrix(const TGeoMatrix &matrix);
    TGeoHMatrix(const char *name);
    ~TGeoHMatrix() override;

geom/geom/inc/TGeoNode.h

@@ -70,6 +70,7 @@ class TGeoNode : public TNamed, public TGeoAtt {
    void Browse(TBrowser *b) override;
    virtual void cd() const {}
    void CheckOverlaps(Double_t ovlp = 0.1, Option_t *option = ""); // *MENU*
+   void CheckOverlapsBySampling(Double_t ovlp = 0.1, Int_t npoints = 1000000); // *MENU*
    void CheckShapes();
    Int_t CountDaughters(Bool_t unique_volumes = kFALSE);
    Int_t DistancetoPrimitive(Int_t px, Int_t py) override;

geom/geom/inc/TGeoOverlapCandidate.h

@@ -0,0 +1,44 @@
+// @(#)root/geom:$Id$
+// Author: Andrei Gheata   05/01/26
+
+/*************************************************************************
+ * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
+ * All rights reserved.                                                  *
+ *                                                                       *
+ * For the licensing terms see $ROOTSYS/LICENSE.                         *
+ * For the list of contributors see $ROOTSYS/README/CREDITS.             *
+ *************************************************************************/
+
+#ifndef ROOT_TGeoOverlapCandidate
+#define ROOT_TGeoOverlapCandidate
+
+#include "TBuffer3D.h"
+#include "TString.h"
+#include "TGeoMatrix.h"
+class TGeoShape;
+class TGeoVolume;
+
+// Lightweight description of a single overlap/extrusion check to perform
+struct TGeoOverlapCandidate {
+   TString fName;               ///< display name
+   TGeoVolume *fVol1 = nullptr; ///< first volume
+   TGeoVolume *fVol2 = nullptr; ///< second volume
+   TGeoHMatrix fMat1;           ///< matrix for first volume
+   TGeoHMatrix fMat2;           ///< matrix for second volume
+   Bool_t fIsOverlap = kTRUE;   ///< kTRUE=overlap, kFALSE=extrusion
+   Double_t fOvlp = 0.0;        ///< threshold for "illegal"
+};
+
+// Output of the numerical check (no ROOT object allocation inside)
+struct TGeoOverlapResult {
+   TString fName;                                ///< display name
+   TGeoVolume *fVol1 = nullptr;                  ///< first volume
+   TGeoVolume *fVol2 = nullptr;                  ///< second volume
+   TGeoHMatrix fMat1;                            /// matrix for first volume
+   TGeoHMatrix fMat2;                            /// matrix for second volume
+   Bool_t fIsOverlap = kTRUE;                    ///< kTRUE=overlap, kFALSE=extrusion
+   Double_t fMaxOverlap = 0.0;                   ///< overlap distance found
+   std::vector<std::array<Double_t, 3>> fPoints; ///< up to N points (e.g. 100)
+};
+
+#endif // ROOT_TGeoOverlapWorkState

geom/geom/inc/TGeoPara.h

@@ -73,6 +73,7 @@ class TGeoPara : public TGeoBBox {
    Double_t GetTxz() const { return fTxz; }
    Double_t GetTyz() const { return fTyz; }
    void InspectShape() const override;
+   Bool_t IsConvex() const final { return kTRUE; }
    Bool_t IsCylType() const override { return kFALSE; }
    Double_t Safety(const Double_t *point, Bool_t in = kTRUE) const override;
    void Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const override;

geom/geom/inc/TGeoParaboloid.h

@@ -63,6 +63,7 @@ class TGeoParaboloid : public TGeoBBox {
    Int_t GetNmeshVertices() const override;
    Bool_t GetPointsOnSegments(Int_t /*npoints*/, Double_t * /*array*/) const override { return kFALSE; }
    void InspectShape() const override;
+   Bool_t IsConvex() const final { return kTRUE; }
    Bool_t IsCylType() const override { return kTRUE; }
    TBuffer3D *MakeBuffer3D() const override;
    Double_t Safety(const Double_t *point, Bool_t in = kTRUE) const override;

geom/geom/inc/TGeoPcon.h

@@ -85,9 +85,12 @@ class TGeoPcon : public TGeoBBox {
    Double_t GetZ(Int_t ipl) const;
    TGeoShape *GetMakeRuntimeShape(TGeoShape * /*mother*/, TGeoMatrix * /*mat*/) const override { return nullptr; }
    Int_t GetNmeshVertices() const override;
-   Bool_t GetPointsOnSegments(Int_t /*npoints*/, Double_t * /*array*/) const override { return kFALSE; }
+   Bool_t GetPointsOnSegments(Int_t npoints, Double_t *array) const override;
    void GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const override;
+   inline Bool_t HasDphi() const { return (TGeoShape::IsSameWithinTolerance(fDphi, 360.)); }
+   inline Bool_t HasRmin(Int_t isec) const { return (fRmin[isec] == 0. && fRmin[isec + 1] == 0.); }
    void InspectShape() const override;
+   Bool_t IsConvex() const final { return (fNz == 2 && !HasDphi() && !HasRmin(0)); }
    Bool_t IsCylType() const override { return kTRUE; }
    TBuffer3D *MakeBuffer3D() const override;
    Double_t &Phi1() { return fPhi1; }

geom/geom/inc/TGeoScaledShape.h

@@ -56,6 +56,7 @@ class TGeoScaledShape : public TGeoBBox {
    TGeoScale *GetScale() const { return fScale; }
    void InspectShape() const override;
    Bool_t IsAssembly() const override;
+   Bool_t IsConvex() const final { return fShape->IsConvex(); }
    Bool_t IsCylType() const override { return fShape->IsCylType(); }
    Bool_t IsReflected() const override;
    TBuffer3D *MakeBuffer3D() const override;

geom/geom/inc/TGeoShape.h

@@ -98,7 +98,7 @@ class TGeoShape : public TNamed {
    static Double_t Tolerance() { return 1.E-10; }
    static Double_t ComputeEpsMch();
    static Double_t EpsMch();
-   virtual void AfterStreamer(){};
+   virtual void AfterStreamer() {};
    virtual Double_t Capacity() const = 0;
    void CheckShape(Int_t testNo, Int_t nsamples = 10000, Option_t *option = "");
    virtual void ComputeBBox() = 0;
@@ -137,6 +137,7 @@ class TGeoShape : public TNamed {
    virtual EInside Inside(const Double_t *point) const;
    virtual Bool_t IsAssembly() const { return kFALSE; }
    virtual Bool_t IsComposite() const { return kFALSE; }
+   virtual Bool_t IsConvex() const { return kFALSE; }
    virtual Bool_t IsCylType() const = 0;
    static Bool_t
    IsCloseToPhi(Double_t epsil, const Double_t *point, Double_t c1, Double_t s1, Double_t c2, Double_t s2);

geom/geom/inc/TGeoShapeAssembly.h

@@ -51,6 +51,7 @@ class TGeoShapeAssembly : public TGeoBBox {
    Int_t GetNmeshVertices() const override { return 0; }
    void InspectShape() const override;
    Bool_t IsAssembly() const override { return kTRUE; }
+   Bool_t IsConvex() const final { return kFALSE; }
    Bool_t IsCylType() const override { return kFALSE; }
    void NeedsBBoxRecompute() { fBBoxOK = kFALSE; }
    void RecomputeBoxLast();

geom/geom/inc/TGeoSphere.h

@@ -75,6 +75,10 @@ class TGeoSphere : public TGeoBBox {
    Double_t GetTheta2() const { return fTheta2; }
    Double_t GetPhi1() const { return fPhi1; }
    Double_t GetPhi2() const { return fPhi2; }
+   Bool_t IsConvex() const override
+   {
+      return !TGeoShape::TestBit(kGeoRSeg) && !TGeoShape::TestBit(kGeoPhiSeg) && !TGeoShape::TestBit(kGeoThetaSeg);
+   }
    void InspectShape() const override;
    Bool_t IsCylType() const override { return kFALSE; }
    Int_t IsOnBoundary(const Double_t *point) const;