From 41925768a29fe72e25d1804a8c55eba85069af76 Mon Sep 17 00:00:00 2001
From: Andrei Gheata <Andrei.Gheata@cern.ch>
Date: Thu, 2 Sep 2021 14:58:04 +0200
Subject: [PATCH] Fix for on surface and moving out condition for grazing rays
 on sphere.

---
 VecGeom/base/Math.h               |  2 ++
 VecGeom/volumes/SphereUtilities.h | 10 +++++-----
 2 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/VecGeom/base/Math.h b/VecGeom/base/Math.h
index 78ab6f931..adbeb282e 100644
--- a/VecGeom/base/Math.h
+++ b/VecGeom/base/Math.h
@@ -17,12 +17,14 @@ inline namespace VECGEOM_IMPL_NAMESPACE {
 #ifdef VECGEOM_FLOAT_PRECISION
 using Precision                        = float;
 VECGEOM_CONST Precision kTolerance     = 1e-3;
+VECGEOM_CONST Precision kSqrtTolerance = 3.1622777e-2;
 VECGEOM_CONST Precision kAngTolerance  = 1e-2;
 VECGEOM_CONST Precision kConeTolerance = 1e-3;
 VECGEOM_CONST Precision kFarAway       = 1e5;
 #else
 using Precision                        = double;
 VECGEOM_CONST Precision kTolerance     = 1e-9;
+VECGEOM_CONST Precision kSqrtTolerance = 3.1622777e-5;
 VECGEOM_CONST Precision kAngTolerance  = 1e-9;
 VECGEOM_CONST Precision kConeTolerance = 1e-7;
 VECGEOM_CONST Precision kFarAway       = 1e10;
diff --git a/VecGeom/volumes/SphereUtilities.h b/VecGeom/volumes/SphereUtilities.h
index 598aa4d3e..7081337c2 100644
--- a/VecGeom/volumes/SphereUtilities.h
+++ b/VecGeom/volumes/SphereUtilities.h
@@ -146,18 +146,18 @@ typename vecCore::Mask_v<Real_v> IsPointOnRadialSurfaceAndMovingOut(UnplacedStru
                                                                     Vector3D<Real_v> const &point,
                                                                     Vector3D<Real_v> const &dir)
 {
-
+  // Rays from rmax+tolerance or rmin-tolerance can be moving out even if not going fully "backward"
   if (MovingOut) {
     if (ForInnerRadius) {
-      return IsPointOnInnerRadius<Real_v>(unplaced, point) && (dir.Dot(-point) > Real_v(0.));
+      return IsPointOnInnerRadius<Real_v>(unplaced, point) && (dir.Dot(-point) > Real_v(kSqrtTolerance));
     } else {
-      return IsPointOnOuterRadius<Real_v>(unplaced, point) && (dir.Dot(point) > Real_v(0.));
+      return IsPointOnOuterRadius<Real_v>(unplaced, point) && (dir.Dot(point) > Real_v(kSqrtTolerance));
     }
   } else {
     if (ForInnerRadius) {
-      return IsPointOnInnerRadius<Real_v>(unplaced, point) && (dir.Dot(-point) < Real_v(0.));
+      return IsPointOnInnerRadius<Real_v>(unplaced, point) && (dir.Dot(-point) < Real_v(-kSqrtTolerance));
     } else
-      return IsPointOnOuterRadius<Real_v>(unplaced, point) && (dir.Dot(point) < Real_v(0.));
+      return IsPointOnOuterRadius<Real_v>(unplaced, point) && (dir.Dot(point) < Real_v(-kSqrtTolerance));
   }
 }
 
-- 
GitLab