Fixes for cut tube DistanceToIn.

VecGeom/volumes/CutPlanes.h

@@ -72,7 +72,6 @@ public:
   VECGEOM_FORCE_INLINE
   VECCORE_ATT_HOST_DEVICE
   void SafetyToOut(Vector3D<Real_v> const &point, Real_v &distance) const;
-
 }; // class CutPlanes
 
 std::ostream &operator<<(std::ostream &os, CutPlanes const &planes);
@@ -132,7 +131,11 @@ void CutPlanes::DistanceToIn(Vector3D<Real_v> const &point, Vector3D<Real_v> con
   // has to be taken
   Real_v d0, d1;
   fCutPlanes[0].DistanceToIn(point, direction, d0);
+  vecCore::MaskedAssign(d0, direction.Dot(Vector3D<Real_v>(fCutPlanes[0].GetNormal())) > Real_v(0.),
+                        Real_v(-kInfLength));
   fCutPlanes[1].DistanceToIn(point, direction, d1);
+  vecCore::MaskedAssign(d1, direction.Dot(Vector3D<Real_v>(fCutPlanes[1].GetNormal())) > Real_v(0.),
+                        Real_v(-kInfLength));
   distance = vecCore::math::Max(d0, d1);
 }
 

VecGeom/volumes/kernel/CutTubeImplementation.h

@@ -200,15 +200,19 @@ void CutTubeImplementation::DistanceToInKernel(UnplacedStruct_t const &unplaced,
 #endif
 
   Bool_v inside;
-  Real_v dplanes;
+  Real_v dplanes, safplanes;
   // Compute distance to cut planes
   unplaced.GetCutPlanes().DistanceToIn<Real_v>(point, direction, dplanes);
+  // We need the safety to cut out points already inside the cut planes.
+  // Unfortunately there is no easy way to get the right normal in vector mode to avoid safety calculation
+  unplaced.GetCutPlanes().SafetyToIn<Real_v>(point, safplanes);
+
   // Mark tracks hitting the planes
-  Bool_v hitplanes = dplanes < stepMax && dplanes > Real_v(-kTolerance);
+  Bool_v hitplanes = vecCore::math::Abs(dplanes) < stepMax && safplanes > Real_v(-kTolerance);
 
 #if USE_CONV_WRONG_SIDE == 1
   if (vecCore::EarlyReturnAllowed()) {
-    if (vecCore::MaskFull((inside_cutplanes == EInside::kOutside) && !hitplanes)) // No particles are hitting
+    if (vecCore::MaskFull((inside_cutplanes != EInside::kInside) && !hitplanes)) // No particles are hitting
       return;
   }
 #endif
@@ -241,7 +245,7 @@ void CutTubeImplementation::DistanceToInKernel(UnplacedStruct_t const &unplaced,
   // All propagated points not yet marked as done should now lie between the cut planes
   // Check now which of the propagated points already entering the solid
   TubeImplementation<TubeTypes::UniversalTube>::Inside<Real_v>(unplaced.GetTubeStruct(), propagated, instart);
-  inside = instart == EInside::kInside;
+  inside = hitplanes && (instart != EInside::kOutside);
   vecCore::MaskedAssign(distance, inside && !done, dplanes);
   done |= inside;
 

VecGeom/volumes/kernel/TubeImplementation.h

@@ -288,7 +288,7 @@ void PhiPlaneTrajectoryIntersection(Precision alongX, Precision alongY, Precisio
     Real_v hity = pos.y() + dist * dir.y();
     Real_v hitz = pos.z() + dist * dir.z();
     Real_v r2   = hitx * hitx + hity * hity;
-    ok &= Abs(hitz) <= tube.fTolIz && (r2 >= tube.fTolIrmin2) && (r2 <= tube.fTolIrmax2);
+    ok &= Abs(hitz) <= tube.fTolOz && (r2 >= tube.fTolOrmin2) && (r2 <= tube.fTolOrmax2);
 
     // GL: tested with this if(PosDirPhiVec) around if(insector), so
     // if(insector){} requires PosDirPhiVec==true to run
@@ -312,11 +312,9 @@ typename vecCore::Mask_v<Real_v> IsOnTubeSurface(UnplacedStruct_t const &tube, V
 {
   const Real_v rho = point.Perp2();
   if (ForInnerSurface) {
-    return (rho >= (tube.fRmin2 - kTolerance * tube.fRmin)) && (rho <= (tube.fRmin2 + kTolerance * tube.fRmin)) &&
-           (Abs(point.z()) < (tube.fZ + kTolerance));
+    return (rho >= tube.fTolOrmin2) && (rho <= tube.fTolIrmin2) && (Abs(point.z()) < (tube.fZ + kTolerance));
   } else {
-    return (rho >= (tube.fRmax2 - kTolerance * tube.fRmax)) && (rho <= (tube.fRmax2 + kTolerance * tube.fRmax)) &&
-           (Abs(point.z()) < (tube.fZ + kTolerance));
+    return (rho >= tube.fTolIrmax2) && (rho <= tube.fTolOrmax2) && (Abs(point.z()) < (tube.fZ + kTolerance));
   }
 }
 
@@ -709,14 +707,14 @@ struct TubeImplementation {
     Real_v crmin = rsq;
 
     // if outside of Rmax, return -1
-    done |= crmax > kTolerance * tube.fRmax;
+    done |= crmax > 2. * kTolerance * tube.fRmax;
     if (vecCore::EarlyReturnAllowed() && vecCore::MaskFull(done)) return;
 
     if (checkRminTreatment<tubeTypeT>(tube)) {
       // if point is within inner-hole of a hollow tube, it is outside of the
       // tube --> return -1
       crmin -= tube.fRmin2; // avoid a division for now
-      done |= crmin < -kTolerance * tube.fRmin;
+      done |= crmin < -2. * kTolerance * tube.fRmin;
       if (vecCore::EarlyReturnAllowed() && vecCore::MaskFull(done)) return;
     }
 

test/unit_tests/TestTube.cpp

@@ -321,33 +321,37 @@ bool TestTubs()
   // std::cout<<"Dist=t5.DistanceToIn((30.0,-20.0,0),vxy) = "<<Dist<<std::endl;
   assert(ApproxEqual<Precision>(Dist, 28.284271));
 
+  // This ray is passing through an edge, may or may not hit depending on rounding
   Dist = t5.DistanceToIn(Vec_t(30.0, -70.0, 0), vxy);
   // std::cout<<"Dist=t5.DistanceToIn((30.0,-70.0,0),vxy) = "<<Dist<<std::endl;
-  assert(ApproxEqual<Precision>(Dist, kInfLength));
+  assert(ApproxEqual<Precision>(Dist, kInfLength) || ApproxEqual<Precision>(Dist, 70. * std::sqrt(2.)));
 
   Dist = t5.DistanceToIn(Vec_t(30.0, -20.0, 0), vmxmy);
   //  std::cout<<"Dist=t5.DistanceToIn((30.0,-20.0,0),vmxmy) = "<<Dist<<std::endl;
   assert(ApproxEqual<Precision>(Dist, 42.426407));
 
+  // This ray is passing through an edge, may or may not hit depending on rounding
   Dist = t5.DistanceToIn(Vec_t(30.0, -70.0, 0), vmxmy);
   // std::cout<<"Dist=t5.DistanceToIn((30.0,-70.0,0),vmxmy) = "<<Dist<<std::endl;
-  assert(ApproxEqual<Precision>(Dist, kInfLength));
+  assert(ApproxEqual<Precision>(Dist, kInfLength) || ApproxEqual<Precision>(Dist, 30. * std::sqrt(2.)));
 
   Dist = t5.DistanceToIn(Vec_t(50.0, -20.0, 0), vy);
   // std::cout<<"Dist=t5.DistanceToIn((50.0,-20.0,0),vy) = "<<Dist<<std::endl;
   assert(ApproxEqual<Precision>(Dist, 20));
 
+  // This ray is passing through an edge, may or may not hit depending on rounding
   Dist = t5.DistanceToIn(Vec_t(100.0, -20.0, 0), vy);
   // std::cout<<"Dist=t5.DistanceToIn((100.0,-20.0,0),vy) = "<<Dist<<std::endl;
-  assert(ApproxEqual<Precision>(Dist, kInfLength));
+  assert(ApproxEqual<Precision>(Dist, kInfLength) || ApproxEqual<Precision>(Dist, 20.));
 
   Dist = t5.DistanceToIn(Vec_t(30.0, -50.0, 0), vmx);
   //  std::cout<<"Dist=t5.DistanceToIn((30.0,-50.0,0),vmx) = "<<Dist<<std::endl;
   assert(ApproxEqual<Precision>(Dist, 30));
 
+  // This ray is passing through an edge, may or may not hit depending on rounding
   Dist = t5.DistanceToIn(Vec_t(30.0, -100.0, 0), vmx);
   //  std::cout<<"Dist=t5.DistanceToIn((30.0,-100.0,0),vmx) = "<<Dist<<std::endl;
-  assert(ApproxEqual<Precision>(Dist, kInfLength));
+  assert(ApproxEqual<Precision>(Dist, kInfLength) || ApproxEqual<Precision>(Dist, 30.));
 
   // ********************************