TGeoParallelWorld improvements

geom/geom/inc/TGeoNavigator.h

@@ -49,6 +49,8 @@ class TGeoNavigator : public TObject {
    Double_t fPoint[3];           //! current point
    Double_t fDirection[3];       //! current direction
    Double_t fLastPoint[3];       //! last point for which safety was computed
+   Double_t fLastPWSaftyPnt[3];  //! last point for which parallel world safety was "evaluated"
+   Double_t fLastPWSafety{-1};   //! last safety returned from parallel world (negative if invalid)
    Int_t fThreadId;              //! thread id for this navigator
    Int_t fLevel;                 //! current geometry level;
    Int_t fNmany;                 //! number of overlapping nodes on current branch
@@ -81,6 +83,8 @@ class TGeoNavigator : public TObject {
    TGeoHMatrix *fDivMatrix;      //! current local matrix of the selected division cell
    TString fPath;                //! path to current node
 
+   static Bool_t fgUsePWSafetyCaching; //! global mode is caching enabled for parallel world safety calls
+
 public:
    TGeoNavigator();
    TGeoNavigator(TGeoManager *geom);
@@ -199,6 +203,35 @@ class TGeoNavigator : public TObject {
       fLastPoint[1] = y, fLastPoint[2] = z;
    }
 
+   // Check if we have a cached safety value from parallel world, and if this can still be used.
+   // Return negative value if no cache available.
+   Double_t GetPWSafetyEstimateFromCache(Double_t cpoint[3]) const
+   {
+      // disregard too small or invalid safeties
+      if (fLastPWSafety < TGeoShape::Tolerance()) {
+         return -1.;
+      }
+      const auto d0 = fLastPWSaftyPnt[0] - cpoint[0];
+      const auto d1 = fLastPWSaftyPnt[1] - cpoint[1];
+      const auto d2 = fLastPWSaftyPnt[2] - cpoint[2];
+      const auto d_sq = d0 * d0 + d1 * d1 + d2 * d2;
+      // if we have moved too much return -1 as "invalid"
+      static constexpr double factor = 1.; // factor between 0 and 1 that "invalidates cache"
+      if (d_sq >= factor * (fLastPWSafety * fLastPWSafety)) {
+         return -1.;
+      }
+      // or return a reasonable cache estimate for safety
+      return fLastPWSafety - std::sqrt(d_sq);
+   }
+
+   // Wrapper for getting the safety from the parallel world.
+   // Takes care of caching mechanics and talking to the Safety function of parallel world.
+   Double_t GetPWSafety(Double_t cpoint[3], Double_t saf_max);
+
+   // enable/disable parallel world safety caching
+   static void SetPWSafetyCaching(Bool_t b) { fgUsePWSafetyCaching = b; }
+   static Bool_t IsPWSafetyCaching() { return fgUsePWSafetyCaching; }
+
    //--- point/vector reference frame conversion
    void LocalToMaster(const Double_t *local, Double_t *master) const { fCache->LocalToMaster(local, master); }
    void LocalToMasterVect(const Double_t *local, Double_t *master) const { fCache->LocalToMasterVect(local, master); }

geom/geom/inc/TGeoParallelWorld.h

@@ -6,19 +6,37 @@
  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
  *************************************************************************/
 
-// Author: Andrei Gheata   30/06/14
+// Authors: Andrei Gheata   30/06/14
+//          Sandro Wenzel   01/09/24
 
 #ifndef ROOT_TGeoParallelWorld
 #define ROOT_TGeoParallelWorld
 
 #include "TNamed.h"
+#include "TGeoVoxelGrid.h"
 
 // forward declarations
 class TGeoManager;
 class TGeoPhysicalNode;
 class TGeoVolume;
 
 class TGeoParallelWorld : public TNamed {
+
+public:
+   // internal enum letting choose between
+   // VoxelFinder or BVH-based algorithms
+   enum class AccelerationMode { kVoxelFinder, kBVH };
+
+   // a structure for safety evaluation (caching) purpose
+   // to be stored per 3D grid voxel
+   struct SafetyVoxelInfo {
+      float min_safety{-1.f}; // the minimum safety from the mid-point of this voxel to any leaf bounding box
+      int idx_start{-1}; // the index into an external storage, where candidate bounding boxes to search for this voxel
+                         // are stored (if -1) --> VoxelInfo not yet initialized
+      unsigned int num_candidates{0}; // the number of candidate bounding boxes to search
+      bool isInitialized() const { return idx_start >= 0; }
+   };
+
 protected:
    TGeoManager *fGeoManager;     // base geometry
    TObjArray *fPaths;            // array of paths
@@ -28,6 +46,14 @@ class TGeoParallelWorld : public TNamed {
    TGeoPhysicalNode *fLastState; //! Last PN touched
    TObjArray *fPhysical;         //! array of physical nodes
 
+   void *fBVH = nullptr; //! BVH helper structure for safety and navigation
+   TGeoVoxelGrid<SafetyVoxelInfo> *fSafetyVoxelCache =
+      nullptr;                                        //! A regular 3D cache layer for fast point-based safety lookups
+   std::vector<unsigned int> fSafetyCandidateStore{}; //! stores bounding boxes serving a quick safety candidates (to be
+                                                      //! used with the VoxelGrid and SafetyVoxelInfo)
+   void *fBoundingBoxes = nullptr;                    //! to keep the vector of primitive axis aligned bounding boxes
+   AccelerationMode fAccMode = AccelerationMode::kVoxelFinder; //! switch between different algorithm implementations
+
    TGeoParallelWorld(const TGeoParallelWorld &) = delete;
    TGeoParallelWorld &operator=(const TGeoParallelWorld &) = delete;
 
@@ -65,10 +91,52 @@ class TGeoParallelWorld : public TNamed {
    // Refresh structures in case of re-alignment
    void RefreshPhysicalNodes();
 
-   // Navigation interface
-   TGeoPhysicalNode *FindNode(Double_t point[3]);
-   TGeoPhysicalNode *FindNextBoundary(Double_t point[3], Double_t dir[3], Double_t &step, Double_t stepmax = 1.E30);
-   Double_t Safety(Double_t point[3], Double_t safmax = 1.E30);
+   // ability to choose algorithm implementation; should be called before CloseGeometry
+   void SetAccelerationMode(AccelerationMode const &mode) { fAccMode = mode; }
+   AccelerationMode const &GetAccelerationMode() const { return fAccMode; }
+
+   // BVH related functions for building, printing, checking
+   void BuildBVH();
+   void PrintBVH() const;
+   bool CheckBVH(void *, size_t) const;
+
+   // --- main navigation interfaces ----
+
+   // FindNode
+   TGeoPhysicalNode *FindNode(Double_t point[3])
+   {
+      switch (fAccMode) {
+      case AccelerationMode::kVoxelFinder: return FindNodeOrig(point);
+      case AccelerationMode::kBVH:
+         return FindNodeBVH(point);
+         // case AccelerationMode::kLoop : return FindNodeLoop(point);
+      }
+      return nullptr;
+   }
+
+   // main interface for Safety
+   Double_t Safety(Double_t point[3], Double_t safmax = 1.E30)
+   {
+      switch (fAccMode) {
+      case AccelerationMode::kVoxelFinder: return SafetyOrig(point, safmax);
+      case AccelerationMode::kBVH:
+         return VoxelSafety(point, safmax);
+         // case AccelerationMode::kLoop : return SafetyLoop(point, safmax);
+      }
+      return 0;
+   }
+
+   // main interface for FindNextBoundary
+   TGeoPhysicalNode *FindNextBoundary(Double_t point[3], Double_t dir[3], Double_t &step, Double_t stepmax = 1.E30)
+   {
+      switch (fAccMode) {
+      case AccelerationMode::kVoxelFinder: return FindNextBoundaryOrig(point, dir, step, stepmax);
+      case AccelerationMode::kBVH:
+         return FindNextBoundaryBVH(point, dir, step, stepmax);
+         // case AccelerationMode::kLoop : return FindNextBoundaryLoop(point, dir, step, stepmax);
+      }
+      return nullptr;
+   }
 
    // Getters
    TGeoManager *GetGeometry() const { return fGeoManager; }
@@ -79,6 +147,40 @@ class TGeoParallelWorld : public TNamed {
    void CheckOverlaps(Double_t ovlp = 0.001); // default 10 microns
    void Draw(Option_t *option) override;
 
+private:
+   // various implementations for FindNextBoundary
+   TGeoPhysicalNode *FindNextBoundaryLoop(Double_t point[3], Double_t dir[3], Double_t &step, Double_t stepmax = 1.E30);
+   TGeoPhysicalNode *FindNextBoundaryOrig(Double_t point[3], Double_t dir[3], Double_t &step, Double_t stepmax = 1.E30);
+   TGeoPhysicalNode *FindNextBoundaryBVH(Double_t point[3], Double_t dir[3], Double_t &step, Double_t stepmax = 1.E30);
+
+   // various implementations for FindNode
+   TGeoPhysicalNode *FindNodeLoop(Double_t point[3]);
+   TGeoPhysicalNode *FindNodeOrig(Double_t point[3]);
+   TGeoPhysicalNode *FindNodeBVH(Double_t point[3]);
+
+   // various implementations for Safety
+   Double_t SafetyLoop(Double_t point[3], Double_t safmax = 1.E30);
+   Double_t SafetyBVH(Double_t point[3], Double_t safmax = 1.E30);
+   Double_t SafetyOrig(Double_t point[3], Double_t safmax = 1.E30);
+   Double_t VoxelSafety(Double_t point[3], Double_t safmax = 1.E30);
+
+   // helper functions related to local safety caching (3D voxel grid)
+
+   // Given a 3D point, returns the
+   // a) the min radius R such that there is at least one leaf bounding box fully enclosed
+   //    in the sphere of radius R around point
+   // b) the set of leaf bounding boxes who partially lie within radius + margin
+
+   // ... using BVH
+   std::pair<double, double>
+   GetBVHSafetyCandidates(double point[3], std::vector<int> &candidates, double margin = 0.) const;
+   // .... same with a simpler, slower algorithm
+   std::pair<double, double>
+   GetLoopSafetyCandidates(double point[3], std::vector<int> &candidates, double margin = 0.) const;
+
+   void InitSafetyVoxel(TGeoVoxelGridIndex const &);
+   void TestVoxelGrid(); // a debug method to play with the voxel grid
+
    ClassDefOverride(TGeoParallelWorld, 3) // parallel world base class
 };
 

geom/geom/inc/TGeoVoxelGrid.h

@@ -0,0 +1,176 @@
+/// \author Sandro Wenzel <sandro.wenzel@cern.ch>
+/// \date 2024-02-22
+
+/*************************************************************************
+ * Copyright (C) 1995-2024, 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_TGeoVoxelGrid
+#define ROOT_TGeoVoxelGrid
+
+// a simple structure to encode voxel indices, to address
+// individual voxels in the 3D grid.
+struct TGeoVoxelGridIndex {
+   int ix{-1};
+   int iy{-1};
+   int iz{-1};
+   size_t idx{std::numeric_limits<size_t>::max()};
+   bool isValid() const { return idx != std::numeric_limits<size_t>::max(); }
+};
+
+/// A finite 3D grid structure, mapping/binning arbitrary 3D cartesian points
+/// onto discrete "voxels". Each such voxel can store an object of type T.
+/// The precision of the voxel binning is done with S (float or double).
+template <typename T, typename S = float>
+class TGeoVoxelGrid {
+public:
+   TGeoVoxelGrid(S xmin, S ymin, S zmin, S xmax, S ymax, S zmax, S Lx_, S Ly_, S Lz_)
+      : fMinBound{xmin, ymin, zmin}, fMaxBound{xmax, ymax, zmax}, fLx(Lx_), fLy(Ly_), fLz(Lz_)
+   {
+
+      // Calculate the number of voxels in each dimension
+      fNx = static_cast<int>((fMaxBound[0] - fMinBound[0]) / fLx);
+      fNy = static_cast<int>((fMaxBound[1] - fMinBound[1]) / fLy);
+      fNz = static_cast<int>((fMaxBound[2] - fMinBound[2]) / fLz);
+
+      finvLx = 1. / fLx;
+      finvLy = 1. / fLy;
+      finvLz = 1. / fLz;
+
+      fHalfDiag = std::sqrt(fLx / 2. * fLx / 2. + fLy / 2. * fLy / 2. + fLz / 2. * fLz / 2.);
+
+      // Resize the grid to hold the voxels
+      fGrid.resize(fNx * fNy * fNz);
+   }
+
+   T &at(int i, int j, int k) { return fGrid[index(i, j, k)]; }
+
+   // check if point is covered by voxel structure
+   bool inside(std::array<S, 3> const &p) const
+   {
+      for (int i = 0; i < 3; ++i) {
+         if (p[i] < fMinBound[i] || p[i] > fMaxBound[i]) {
+            return false;
+         }
+      }
+      return true;
+   }
+
+   // Access a voxel given a 3D point P
+   T &at(std::array<S, 3> const &P)
+   {
+      int i, j, k;
+      pointToVoxelIndex(P, i, j, k); // Convert point to voxel index
+      return fGrid[index(i, j, k)];  // Return reference to voxel's data
+   }
+
+   T *at(TGeoVoxelGridIndex const &vi)
+   {
+      if (!vi.isValid()) {
+         return nullptr;
+      }
+      return &fGrid[vi.idx];
+   }
+
+   // Set the data of a voxel at point P
+   void set(std::array<S, 3> const &p, const T &value)
+   {
+      int i, j, k;
+      pointToVoxelIndex(p, i, j, k); // Convert point to voxel index
+      fGrid[index(i, j, k)] = value; // Set the value at the voxel
+   }
+
+   // Set the data of a voxel at point P
+   void set(int i, int j, int k, const T &value)
+   {
+      fGrid[index(i, j, k)] = value; // Set the value at the voxel
+   }
+
+   void set(TGeoVoxelGridIndex const &vi, const T &value) { fGrid[vi.idx] = value; }
+
+   // Get voxel dimensions
+   int getVoxelCountX() const { return fNx; }
+   int getVoxelCountY() const { return fNy; }
+   int getVoxelCountZ() const { return fNz; }
+
+   // returns the cartesian mid-point coordinates of a voxel given by a VoxelIndex
+   std::array<S, 3> getVoxelMidpoint(TGeoVoxelGridIndex const &vi) const
+   {
+      const S midX = fMinBound[0] + (vi.ix + 0.5) * fLx;
+      const S midY = fMinBound[1] + (vi.iy + 0.5) * fLy;
+      const S midZ = fMinBound[2] + (vi.iz + 0.5) * fLz;
+
+      return {midX, midY, midZ};
+   }
+
+   S getDiagonalLength() const { return fHalfDiag; }
+
+   // Convert a point p(x, y, z) to voxel indices (i, j, k)
+   // if point is outside set indices i,j,k to -1
+   void pointToVoxelIndex(std::array<S, 3> const &p, int &i, int &j, int &k) const
+   {
+      if (!inside(p)) {
+         i = -1;
+         j = -1;
+         k = -1;
+      }
+
+      i = static_cast<int>((p[0] - fMinBound[0]) * finvLx);
+      j = static_cast<int>((p[1] - fMinBound[1]) * finvLy);
+      k = static_cast<int>((p[2] - fMinBound[2]) * finvLz);
+
+      // Clamp the indices to valid ranges
+      i = std::min(i, fNx - 1);
+      j = std::min(j, fNy - 1);
+      k = std::min(k, fNz - 1);
+   }
+
+   // Convert a point p(x, y, z) to voxel index object
+   // if outside, an invalid index object will be returned
+   TGeoVoxelGridIndex pointToVoxelIndex(std::array<S, 3> const &p) const
+   {
+      if (!inside(p)) {
+         return TGeoVoxelGridIndex(); // invalid voxel index
+      }
+
+      int i = static_cast<int>((p[0] - fMinBound[0]) * finvLx);
+      int j = static_cast<int>((p[1] - fMinBound[1]) * finvLy);
+      int k = static_cast<int>((p[2] - fMinBound[2]) * finvLz);
+
+      // Clamp the indices to valid ranges
+      i = std::min(i, fNz - 1);
+      j = std::min(j, fNy - 1);
+      k = std::min(k, fNz - 1);
+
+      return TGeoVoxelGridIndex{i, j, k, index(i, j, k)};
+   }
+
+   TGeoVoxelGridIndex pointToVoxelIndex(S x, S y, S z) const { return pointToVoxelIndex(std::array<S, 3>{x, y, z}); }
+
+   // Convert voxel indices (i, j, k) to a linear index in the grid array
+   size_t index(int i, int j, int k) const { return i + fNx * (j + fNy * k); }
+
+   void indexToIndices(size_t idx, int &i, int &j, int &k) const
+   {
+      k = idx % fNz;
+      j = (idx / fNz) % fNy;
+      i = idx / (fNy * fNz);
+   }
+
+   // member data
+
+   std::array<S, 3> fMinBound;
+   std::array<S, 3> fMaxBound; // 3D bounds for grid structure
+   S fLx, fLy, fLz;            // Voxel dimensions
+   S finvLx, finvLy, finvLz;   // inverse voxel dimensions
+   S fHalfDiag;                // cached value for voxel half diagonal length
+
+   int fNx, fNy, fNz;    // Number of voxels in each dimension
+   std::vector<T> fGrid; // The actual voxel grid data container
+};
+
+#endif