Better use of safety.

examples/CMakeLists.txt

@@ -9,4 +9,4 @@ configure_file(data/ppttbar.hepmc3 ${PROJECT_BINARY_DIR}/ppttbar.hepmc3)
 # - Subprojects
 add_subdirectory(Example1)
 add_subdirectory(IntegrationBenchmark)
-add_subdirectory(AsyncExample)
+#add_subdirectory(AsyncExample)

include/AdePT/core/Track.cuh

@@ -31,9 +31,36 @@ struct Track {
   double localTime{0};
   double properTime{0};
 
-  vecgeom::Vector3D<Precision> pos;
-  vecgeom::Vector3D<Precision> dir;
-  vecgeom::NavigationState navState;
+  vecgeom::Vector3D<Precision> pos;   ///< track position
+  vecgeom::Vector3D<Precision> dir;   ///< track direction
+  vecgeom::Vector3D<float> safetyPos; ///< last position where the safety was computed
+  float safety{0};                    ///< last computed safety value
+  vecgeom::NavigationState navState;  ///< current navigation state
+
+  /// @brief Get recomputed cached safety ay a given track position
+  /// @param new_pos Track position
+  /// @param accurate_limit Only return non-zero if the recomputed safety if larger than the accurate_limit
+  /// @return Recomputed safety.
+  __host__ __device__ VECGEOM_FORCE_INLINE 
+  float GetSafety(vecgeom::Vector3D<Precision> const &new_pos, float accurate_limit = 0) const
+  {
+    float dsafe = safety - accurate_limit;
+    if (dsafe <= 0) return 0;
+    float distSq = (vecgeom::Vector3D<float>(new_pos) - safetyPos).Mag2();
+    if (dsafe * dsafe < distSq) return 0.;
+    return (safety - vecCore::math::Sqrt(distSq));
+  }
+
+  __host__ __device__ VECGEOM_FORCE_INLINE
+
+  /// @brief Set Safety value computed in a new point
+  /// @param new_pos Position where the safety is computed
+  /// @param safe Safety value
+  void SetSafety(vecgeom::Vector3D<Precision> const &new_pos, float safe)
+  {
+    safetyPos.Set(static_cast<float>(new_pos[0]), static_cast<float>(new_pos[1]), static_cast<float>(new_pos[2]));
+    safety = vecCore::math::Max(safe, 0.f);
+  }
 
 #ifdef USE_SPLIT_KERNELS
   // Variables used to store track info needed for scoring

include/AdePT/kernels/electrons.cuh

@@ -43,6 +43,8 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
                                                           Scoring *userScoring, VolAuxData const *auxDataArray)
 {
   using namespace adept_impl;
+  // #define DEBUG_NO_SECONDARY 1
+
 #ifdef VECGEOM_FLOAT_PRECISION
   const Precision kPush = 10 * vecgeom::kTolerance;
 #else
@@ -110,14 +112,6 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
     // used, this provides a fresh round of random numbers and reduces thread
     // divergence because the RNG state doesn't need to be advanced later.
     RanluxppDouble newRNG(currentTrack.rngState.BranchNoAdvance());
-
-    // Compute safety, needed for MSC step limit.
-    double safety = 0;
-    if (!navState.IsOnBoundary()) {
-      safety = AdePTNavigator::ComputeSafety(pos, navState);
-    }
-    theTrack->SetSafety(safety);
-
     G4HepEmRandomEngine rnge(&currentTrack.rngState);
 
     // Sample the `number-of-interaction-left` and put it into the track.
@@ -131,6 +125,33 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
 
     G4HepEmElectronManager::HowFarToDiscreteInteraction(&g4HepEmData, &g4HepEmPars, &elTrack);
 
+    auto physicalStepLength = elTrack.GetPStepLength();
+#ifdef DEBUG_NO_SECONDARY
+    printf("== inside: ");
+    navState.Print();
+    printf("   on_boundary = %d at pos: %g, %g, %g  and dir: %g, %g, %g  ekin = %g  phys_step = %g\n",
+           navState.IsOnBoundary(), pos[0], pos[1], pos[2], dir[0], dir[1], dir[2], eKin, physicalStepLength);
+#endif
+    // Compute safety, needed for MSC step limit. The accuracy range is physicalStepLength
+    double safety = 0.;
+    if (!navState.IsOnBoundary()) {
+      // Get the remaining safety only if larger than physicalStepLength
+      safety = currentTrack.GetSafety(pos, physicalStepLength);
+      if (safety < physicalStepLength) {
+        // Recompute safety and update it in the track.
+#ifdef ADEPT_USE_SURF
+        // Use maximum accuracy only if safety is samller than physicalStepLength
+        safety = AdePTNavigator::ComputeSafety(pos, navState, physicalStepLength);
+#else
+        safety = AdePTNavigator::ComputeSafety(pos, navState);
+#endif
+      }
+    }
+    currentTrack.SetSafety(pos, safety);
+    theTrack->SetSafety(safety);
+#ifdef DEBUG_NO_SECONDARY
+    printf("   safety = %g\n", safety);
+#endif
     bool restrictedPhysicalStepLength = false;
     if (BzFieldValue != 0) {
       const double momentumMag = sqrt(eKin * (eKin + 2.0 * restMass));
@@ -141,11 +162,14 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
       double limit                = MaxSafeLength * safeLength;
       limit                       = safety > limit ? safety : limit;
 
-      double physicalStepLength = elTrack.GetPStepLength();
       if (physicalStepLength > limit) {
-        physicalStepLength           = limit;
+        physicalStepLength = limit;
+#ifdef DEBUG_NO_SECONDARY
+        printf("   physics step resticted to %g\n", limit);
+#endif
         restrictedPhysicalStepLength = true;
         elTrack.SetPStepLength(physicalStepLength);
+
         // Note: We are limiting the true step length, which is converted to
         // a shorter geometry step length in HowFarToMSC. In that sense, the
         // limit is an over-approximation, but that is fine for our purpose.
@@ -179,6 +203,10 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
       geometryStepLength = fieldPropagatorBz.ComputeStepAndNextVolume<AdePTNavigator>(
           eKin, restMass, Charge, geometricalStepLengthFromPhysics, pos, dir, navState, nextState, hitsurf_index,
           propagated, safety);
+#ifdef DEBUG_NO_SECONDARY
+      printf("   field step = %g  to pos = (%g, %g, %g)  on_boundary = %d\n", geometryStepLength, pos[0], pos[1],
+             pos[2], nextState.IsOnBoundary());
+#endif
     } else {
 #ifdef ADEPT_USE_SURF
       geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics,
@@ -194,6 +222,7 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
     // correct information (navState = nextState only if relocated
     // in case of a boundary; see below)
     navState.SetBoundaryState(nextState.IsOnBoundary());
+    if (nextState.IsOnBoundary()) currentTrack.SetSafety(pos, 0.);
 
     // Propagate information from geometrical step to MSC.
     theTrack->SetDirection(dir.x(), dir.y(), dir.z());
@@ -215,7 +244,8 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
       if (dLength2 > kGeomMinLength2) {
         const double dispR = std::sqrt(dLength2);
         // Estimate safety by subtracting the geometrical step length.
-        safety -= geometryStepLength;
+        // safety -= geometryStepLength;
+        safety                 = currentTrack.GetSafety(pos);
         constexpr double sFact = 0.99;
         double reducedSafety   = sFact * safety;
 
@@ -225,7 +255,13 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
           pos += displacement;
         } else {
           // Recompute safety.
-          safety        = AdePTNavigator::ComputeSafety(pos, navState);
+#ifdef ADEPT_USE_SURF
+          // Use maximum accuracy only if safety is samller than physicalStepLength
+          safety = AdePTNavigator::ComputeSafety(pos, navState, dispR);
+#else
+          safety = AdePTNavigator::ComputeSafety(pos, navState);
+#endif
+          currentTrack.SetSafety(pos, safety);
           reducedSafety = sFact * safety;
 
           // 1b. Far away from geometry boundary:
@@ -243,6 +279,9 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
     // Collect the charged step length (might be changed by MSC). Collect the changes in energy and deposit.
     eKin                 = theTrack->GetEKin();
     double energyDeposit = theTrack->GetEnergyDeposit();
+#ifdef DEBUG_NO_SECONDARY
+    printf("   after MSC eKin = %g eDep = %g pos = {%g, %g, %g}\n", eKin, energyDeposit, pos[0], pos[1], pos[2]);
+#endif
 
     // Update the flight times of the particle
     // By calculating the velocity here, we assume that all the energy deposit is done at the PreStepPoint, and
@@ -262,6 +301,10 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
         AdePTNavigator::RelocateToNextVolume(pos, dir, nextState);
 #endif
       }
+#ifdef DEBUG_NO_SECONDARY
+      printf("   after relocation: ");
+      nextState.Print();
+#endif
     }
 
     if (auxData.fSensIndex >= 0)
@@ -290,6 +333,9 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
     }
 
     if (stopped) {
+#ifdef DEBUG_NO_SECONDARY
+      continue;
+#endif
       if (!IsElectron) {
         // Annihilate the stopped positron into two gammas heading to opposite
         // directions (isotropic).
@@ -388,6 +434,7 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
       double dirSecondary[3];
       G4HepEmElectronInteractionIoni::SampleDirections(eKin, deltaEkin, dirSecondary, dirPrimary, &rnge);
 
+#ifndef DEBUG_NO_SECONDARY
       Track &secondary = secondaries.electrons->NextTrack();
 
       adept_scoring::AccountProduced(userScoring, /*numElectrons*/ 1, /*numPositrons*/ 0, /*numGammas*/ 0);
@@ -397,6 +444,7 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
       secondary.rngState = newRNG;
       secondary.eKin     = deltaEkin;
       secondary.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]);
+#endif
 
       eKin -= deltaEkin;
       dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]);
@@ -416,6 +464,7 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
       double dirSecondary[3];
       G4HepEmElectronInteractionBrem::SampleDirections(eKin, deltaEkin, dirSecondary, dirPrimary, &rnge);
 
+#ifndef DEBUG_NO_SECONDARY
       Track &gamma = secondaries.gammas->NextTrack();
       adept_scoring::AccountProduced(userScoring, /*numElectrons*/ 0, /*numPositrons*/ 0, /*numGammas*/ 1);
 
@@ -424,6 +473,7 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
       gamma.rngState = newRNG;
       gamma.eKin     = deltaEkin;
       gamma.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]);
+#endif
 
       eKin -= deltaEkin;
       dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]);

include/AdePT/magneticfield/fieldPropagatorConstBz.h

@@ -142,7 +142,12 @@ __host__ __device__ Precision fieldPropagatorConstBz::ComputeStepAndNextVolume(
     } else {
       Precision newSafety = 0;
       if (stepDone > 0) {
+#ifdef ADEPT_USE_SURF
+        // Use maximum accuracy only if safety is samller than physicalStepLength
+        newSafety = Navigator::ComputeSafety(position, current_state, remains);
+#else
         newSafety = Navigator::ComputeSafety(position, current_state);
+#endif
       }
       if (newSafety > chordLen) {
         move         = chordLen;

include/AdePT/navigation/SurfNavigator.h

@@ -49,10 +49,12 @@ class SurfNavigator {
   /// @brief Computes the isotropic safety from the globalpoint.
   /// @param globalpoint Point in global coordinates
   /// @param state Path where to compute safety
+  /// @param limit Limit to which safty should be computed accurately
   /// @return Isotropic safe distance
-  __host__ __device__ static Precision ComputeSafety(Vector3D const &globalpoint, vecgeom::NavigationState const &state)
+  __host__ __device__ static Precision ComputeSafety(Vector3D const &globalpoint, vecgeom::NavigationState const &state,
+                                                     Precision limit = vecgeom::InfinityLength<Real_t>())
   {
-    auto safety = vgbrep::protonav::BVHSurfNavigator<Real_t>::ComputeSafety(globalpoint, state);
+    auto safety = vgbrep::protonav::BVHSurfNavigator<Real_t>::ComputeSafety(globalpoint, state, limit);
     return safety;
   }