Precision fixes for the cut tubes.

VecGeom/volumes/kernel/CutTubeImplementation.h

@@ -208,7 +208,7 @@ void CutTubeImplementation::DistanceToInKernel(UnplacedStruct_t const &unplaced,
 
 #if USE_CONV_WRONG_SIDE == 1
   if (vecCore::EarlyReturnAllowed()) {
-    if (vecCore::MaskFull((inside_cutplanes != EInside::kInside) && !hitplanes)) // No particles are hitting
+    if (vecCore::MaskFull((inside_cutplanes == EInside::kOutside) && !hitplanes)) // No particles are hitting
       return;
   }
 #endif
@@ -262,6 +262,9 @@ void CutTubeImplementation::DistanceToInKernel(UnplacedStruct_t const &unplaced,
   Real_v dtube = InfinityLength<Real_v>();
   TubeImplementation<TubeTypes::UniversalTube>::DistanceToInKernel<Real_v>(unplaced.GetTubeStruct(), propagated,
                                                                            direction, stepMax, dtube);
+  // A.G Propagation to cut planes can put the point inside the tube, so DistanceToIn may return -1
+  // In such case we need to set dtube to 0, otherwise we may get wrong negative answers
+  vecCore__MaskedAssignFunc(dtube, dtube < 0., Real_v(0.));
   vecCore__MaskedAssignFunc(dtube, dexit < dtube, InfinityLength<Real_v>());
   vecCore__MaskedAssignFunc(distance, !done && (dtube + dplanes) < stepMax, dtube + dplanes);
 // The line below is needed for the convention
@@ -286,8 +289,8 @@ void CutTubeImplementation::DistanceToOut(UnplacedStruct_t const &unplaced, Vect
   TubeImplementation<TubeTypes::UniversalTube>::DistanceToOut<Real_v>(unplaced.GetTubeStruct(), point, direction,
                                                                       stepMax, dtube);
   vecCore::MaskedAssign(distance, dtube < distance, dtube);
-  // The line below is needed for the start on boundary convention
-  //  vecCore__MaskedAssignFunc(distance, vecCore::math::Abs(distance) < Real_v(kTolerance), Real_v(0.));
+  // The line below is needed to avoid din=dout=0 when starting from a boundary
+  vecCore__MaskedAssignFunc(distance, distance >= Real_v(0.) && distance < Real_v(kTolerance), Real_v(kTolerance));
 }
 
 //______________________________________________________________________________

VecGeom/volumes/kernel/TubeImplementation.h

@@ -269,16 +269,18 @@ void PhiPlaneTrajectoryIntersection(Precision alongX, Precision alongY, Precisio
   // approaching phi plane from the right side?
   // this depends whether we use it for DistanceToIn or DistanceToOut
   // Note: wedge normals poing towards the wedge inside, by convention!
+  Real_v dirDotNorm = dir.x() * normalX + dir.y() * normalY;
   if (insectorCheck)
-    ok = (dir.x() * normalX + dir.y() * normalY > Real_v(0.)); // DistToIn  -- require tracks entering volume
+    ok = (dirDotNorm > Real_v(0.)); // DistToIn  -- require tracks entering volume
   else
-    ok = (dir.x() * normalX + dir.y() * normalY < Real_v(0.)); // DistToOut -- require tracks leaving volume
+    ok = (dirDotNorm < Real_v(0.)); // DistToOut -- require tracks leaving volume
 
   // if( vecCore::EarlyReturnAllowed() && vecCore::MaskEmpty(ok) ) return;
 
   Real_v dirDotXY = (dir.y() * alongX - dir.x() * alongY);
   dist            = (alongY * pos.x() - alongX * pos.y()) / NonZero(dirDotXY);
-  ok &= dist > -kHalfTolerance;
+  // A.G to check validity, we have to compare with tolerance the safety rather than the distance to plane
+  ok &= (dist * Abs(dirDotNorm)) > -kHalfTolerance;
   // if( vecCore::EarlyReturnAllowed() && vecCore::MaskEmpty(ok) ) return;
 
   if (insectorCheck) {

test/shape_tester/ShapeTester.cpp

@@ -788,6 +788,8 @@ int ShapeTester<ImplT>::TestSurfacePoint()
     distOut = CallDistanceToOut(fVolume, point, v, normal, convex);
 
     if (distIn == 0. && distOut == 0.) {
+    distIn  = fVolume->DistanceToIn(point, v);
+    distOut = CallDistanceToOut(fVolume, point, v, normal, convex);
       icount1++;
       ReportError(&nError, point, v, 0., "TS: DistanceToIn=DistanceToOut=0 for point on Surface");
     }
@@ -1086,7 +1088,7 @@ int ShapeTester<ImplT>::TestOutsidePoint()
       if (dist >= kInfLength) {
         ReportError(&nError, p, v, dist, "TO2: DistanceToOut(p,v) == kInfLength");
         continue;
-      } else if (dist < 0) {
+      } else if (dist < -fSolidTolerance) { // Not an error if distance is negative
         ReportError(&nError, p, v, dist, "TO2: DistanceToOut(p,v) < 0");
         continue;
       }

test/shape_tester/shape_testCutTube.cpp

@@ -57,10 +57,9 @@ int main(int argc, char *argv[])
   tester.setDebug(debug);
   tester.setStat(stat);
   tester.SetMaxPoints(npoints);
-  tester.SetTestBoundaryErrors(true);
-  #ifdef VECGEOM_FLOAT_PRECISION
-    tester.SetSolidTolerance(1e-4);
-  #endif
+  //tester.SetTestBoundaryErrors(true);
+  tester.SetSolidTolerance(vecgeom::kTolerance);
+
   int errCode = tester.Run(cuttube);
 
   std::cout << "Final Error count for Shape *** " << cuttube->GetName() << "*** = " << errCode << "\n";