From dc944fd88aaa2d0d8f67e1f6e89981e543e6f6c2 Mon Sep 17 00:00:00 2001 From: Andrei Gheata Date: Fri, 29 Nov 2024 15:57:38 +0100 Subject: [PATCH] Better use of safety. Implementing two safety features: 1. Safety caching: the tracks hold the point where safety is computed, and the safety value in that point. When safety is needed in a different point, the method track.GetSafety(new_point) will subtract from the cached value and return the positive remainder or zero otherwise. The GetSafety method can take an accurate limit value, and when the remainder is smaller than this limit the returned safety is zero. It is the user responsibility to call track.SetSafety(pos, safety) whenever the safety is fully recomputed in a new position. 2. Safety calculation usses a new surface model feature to only compute safety accurately when close to boundaries, and only use the aligned bounding box safety when far away. The distance limit is passed as third parameter to AdePTNavigator::ComputeSafety, and must be typically equal to the discrete interaction step. --- include/AdePT/core/Track.cuh | 56 ++++++++++++++----- include/AdePT/kernels/electrons.cuh | 39 +++++++++---- .../magneticfield/fieldPropagatorConstBz.h | 5 ++ include/AdePT/navigation/SurfNavigator.h | 6 +- 4 files changed, 78 insertions(+), 28 deletions(-) diff --git a/include/AdePT/core/Track.cuh b/include/AdePT/core/Track.cuh index a1e56155..2cbe2abf 100644 --- a/include/AdePT/core/Track.cuh +++ b/include/AdePT/core/Track.cuh @@ -21,19 +21,21 @@ struct Track { using Precision = vecgeom::Precision; RanluxppDouble rngState; - double eKin; - double numIALeft[4]; - double initialRange; - double dynamicRangeFactor; - double tlimitMin; - - double globalTime{0}; - double localTime{0}; - double properTime{0}; - - vecgeom::Vector3D pos; - vecgeom::Vector3D dir; - vecgeom::NavigationState navState; + double eKin{0.}; + double numIALeft[4]{0., 0., 0., 0.}; + double initialRange{0.}; + double dynamicRangeFactor{0.}; + double tlimitMin{0.}; + + double globalTime{0.}; + double localTime{0.}; + double properTime{0.}; + + vecgeom::Vector3D pos; ///< track position + vecgeom::Vector3D dir; ///< track direction + vecgeom::Vector3D safetyPos; ///< last position where the safety was computed + float safety{0.f}; ///< last computed safety value + vecgeom::NavigationState navState; ///< current navigation state #ifdef USE_SPLIT_KERNELS // Variables used to store track info needed for scoring @@ -46,7 +48,6 @@ struct Track { // Variables used to store navigation results double geometryStepLength{0}; - double safety{0}; long hitsurfID{0}; #endif @@ -60,6 +61,29 @@ struct Track { bool stopped{false}; #endif + /// @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 const &new_pos, + float accurate_limit = 0.f) const + { + float dsafe = safety - accurate_limit; + if (dsafe <= 0.f) return 0.f; + float distSq = (vecgeom::Vector3D(new_pos) - safetyPos).Mag2(); + if (dsafe * dsafe < distSq) return 0.f; + return (safety - vecCore::math::Sqrt(distSq)); + } + + /// @brief Set Safety value computed in a new point + /// @param new_pos Position where the safety is computed + /// @param safe Safety value + __host__ __device__ VECGEOM_FORCE_INLINE void SetSafety(vecgeom::Vector3D const &new_pos, float safe) + { + safetyPos.Set(static_cast(new_pos[0]), static_cast(new_pos[1]), static_cast(new_pos[2])); + safety = vecCore::math::Max(safe, 0.f); + } + __host__ __device__ double Uniform() { return rngState.Rndm(); } __host__ __device__ void InitAsSecondary(const vecgeom::Vector3D &parentPos, @@ -77,7 +101,9 @@ struct Track { // A secondary inherits the position of its parent; the caller is responsible // to update the directions. - this->pos = parentPos; + this->pos = parentPos; + this->safetyPos.Set(0.f, 0.f, 0.f); + this->safety = 0.0f; this->navState = parentNavState; // The global time is inherited from the parent diff --git a/include/AdePT/kernels/electrons.cuh b/include/AdePT/kernels/electrons.cuh index 65762869..4a2a1f25 100644 --- a/include/AdePT/kernels/electrons.cuh +++ b/include/AdePT/kernels/electrons.cuh @@ -110,14 +110,6 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManagerSetSafety(safety); - G4HepEmRandomEngine rnge(¤tTrack.rngState); // Sample the `number-of-interaction-left` and put it into the track. @@ -132,6 +124,24 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManagerSetSafety(safety); bool restrictedPhysicalStepLength = false; if (BzFieldValue != 0) { const double momentumMag = sqrt(eKin * (eKin + 2.0 * restMass)); @@ -142,11 +152,11 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager limit ? safety : limit; - double physicalStepLength = elTrack.GetPStepLength(); if (physicalStepLength > limit) { physicalStepLength = limit; 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. @@ -201,6 +211,7 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManagerSetDirection(dir.x(), dir.y(), dir.z()); @@ -222,7 +233,7 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager kGeomMinLength2) { const double dispR = std::sqrt(dLength2); // Estimate safety by subtracting the geometrical step length. - safety -= geometryStepLength; + safety = currentTrack.GetSafety(pos); constexpr double sFact = 0.99; double reducedSafety = sFact * safety; @@ -232,7 +243,13 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager 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; diff --git a/include/AdePT/navigation/SurfNavigator.h b/include/AdePT/navigation/SurfNavigator.h index b5a10d70..d71449bc 100644 --- a/include/AdePT/navigation/SurfNavigator.h +++ b/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 safety 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()) { - auto safety = vgbrep::protonav::BVHSurfNavigator::ComputeSafety(globalpoint, state); + auto safety = vgbrep::protonav::BVHSurfNavigator::ComputeSafety(globalpoint, state, limit); return safety; }