[geom] Add cache layer for TGeoParallelWorld safeties in TGeoNavigator (

#16470)

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,34 @@ 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"
+      if (d_sq >= (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/src/TGeoNavigator.cxx

@@ -47,6 +47,8 @@ static Double_t gTolerance = TGeoShape::Tolerance();
 const char *kGeoOutsidePath = " ";
 const Int_t kN3 = 3 * sizeof(Double_t);
 
+Bool_t TGeoNavigator::fgUsePWSafetyCaching = kFALSE;
+
 ClassImp(TGeoNavigator);
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -56,6 +58,7 @@ TGeoNavigator::TGeoNavigator()
    : fStep(0.),
      fSafety(0.),
      fLastSafety(0.),
+     fLastPWSafety(-1.),
      fThreadId(0),
      fLevel(0),
      fNmany(0),
@@ -107,6 +110,7 @@ TGeoNavigator::TGeoNavigator(TGeoManager *geom)
    : fStep(0.),
      fSafety(0.),
      fLastSafety(0.),
+     fLastPWSafety(-1.),
      fThreadId(0),
      fLevel(0),
      fNmany(0),
@@ -1645,6 +1649,39 @@ void TGeoNavigator::ResetState()
    fIsStepEntering = fIsStepExiting = kFALSE;
 }
 
+//////////////////////////////////////////////////////////////////////////////////
+/// Wrapper for getting the safety from the parallel world. Takes care of
+/// caching mechanics + talking to the parallel world.
+
+double TGeoNavigator::GetPWSafety(Double_t cpoint[3], Double_t saf_max)
+{
+   if (!IsPWSafetyCaching()) {
+      return fGeometry->GetParallelWorld()->Safety(cpoint, saf_max);
+   }
+   const auto cached = GetPWSafetyEstimateFromCache(cpoint);
+   if (cached > 0) {
+      // if cache is valid, just use it
+      return cached;
+   }
+   // otherwise we need to evaluate it and update the cache
+   // we evaluate this with saf_max = infinity to get the best
+   // possible safety value
+   auto pw = fGeometry->GetParallelWorld();
+   const auto newsafety = pw->Safety(cpoint /*saf_max*/);
+
+   // we need to be a bit careful: A returned safety value of TGeoShape::Big()
+   // is not the actual safety and should not be cached
+   if (newsafety < TGeoShape::Big()) {
+      fLastPWSafety = newsafety;
+      fLastPWSaftyPnt[0] = cpoint[0];
+      fLastPWSaftyPnt[1] = cpoint[1];
+      fLastPWSaftyPnt[2] = cpoint[2];
+   } else {
+      fLastPWSafety = -1;
+   }
+   return newsafety;
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 /// Compute safe distance from the current point. This represent the distance
 /// from POINT to the closest boundary.
@@ -1673,7 +1710,7 @@ Double_t TGeoNavigator::Safety(Bool_t inside)
 
       // cross-check against the parallel world safety, using fSafety as limit
       if (have_PW) {
-         safpar = fGeometry->GetParallelWorld()->Safety(fPoint, fSafety);
+         safpar = GetPWSafety(fPoint, fSafety);
       }
       return TMath::Min(fSafety, safpar);
    }
@@ -1695,7 +1732,7 @@ Double_t TGeoNavigator::Safety(Bool_t inside)
    //---> Check against the parallel geometry safety
    // cross-check against the parallel world safety, using fSafety as limit
    if (have_PW) {
-      safpar = fGeometry->GetParallelWorld()->Safety(fPoint, fSafety);
+      safpar = GetPWSafety(fPoint, fSafety);
    }
    if (safpar < fSafety)
       fSafety = safpar;