[VECGEOM-564] Adapt FillUncontainedPoints for "full" volumes, add 2nd...

[VECGEOM-564] Adapt FillUncontainedPoints for "full" volumes, add 2nd CheckDaughterIntersections in NewSimpleNavigator

CMakeLists.txt

@@ -249,7 +249,7 @@ endif()
 
 ################################################################################
 # Minimum version of VecCore we need.
-set(VecCore_VERSION "0.5.2")
+set(VecCore_VERSION "0.8.0")
 
 # Let's see if we can find VecCore
 if(NOT BUILTIN_VECCORE)
@@ -468,7 +468,6 @@ endif()
 if(VECGEOM_ROOT)
   # ROOT install may be relying on ROOTSYS
   list(APPEND CMAKE_PREFIX_PATH $ENV{ROOTSYS})
-
   find_package(ROOT REQUIRED COMPONENTS Geom Graf3d)
 
   #---setup ROOT include + lib dirs

ReleaseNotes.txt

+
+Version 1.1.18       5th November 2021
+======================================
+Overview: 
+- Further implementation of functions needed for inclusion of VecGeom Navigation into Geant4.
+- Moved included version of VecCore to 0.8.0
+- Use precision-dependent push rather than fixed one.
+- Fix for GeometryTest in CUDA + vector mode
+- More improvements and fixes to multiple shapes for improved test results in single precision mode
+- Fixes for single-precision mode on GPU.
+
+Migration / configuration / behaviour change(s):
+- Added arb8 to gdml reader
+- Export required flags for VecGeom as compile options
+- Disable usage of fast math for nvcc unless requested via FAST_MATH option.
+
+Details:
+- VNavigators (HybridNavigator2, NewSimpleNavigator) improved
+- VolumeUtilities::FillUncontainedPoints revised to cope with no uncontained volume
+- Improvements for different compilers: Intel, LLVM, Clang 13.0
+- Type consistency fixes.
+- More robust type casting for utility methods.
+- Compilation fixes for single-precision vector mode.
+- Use vecCore::math function rather than standard ones, for consistent host/device behaviour.
+- Improved cone performance: avoid a few redundant IsOnConicalSurface(...) function calls
+- Add Normal() CPU implementation also to GPU
+- Added doxygen docs for ApproachNextDaughter
+- BVH: Declare stacks as unsigned int and reduced code duplication of LevelLocate
+
+Version 1.1.17    3rd September 2021
+===================================
+
+Overview:
+- Improvements and fixes to multiple shapes to enable tests to pass in single precision mode
+- Fixes for BVH
+- Added accessors for navigation state index.
+
+Migration - configuration or interface changes:
+- Removed navigation index setter since the state should rather be recomputed for consistency.
+
+Details:
+- Fix for on surface and moving out condition for grazing rays on sphere.
+- Fix for BVH::ApproachNextDaughter: point/dir conversion to checked daughter frame.
+- Fixes for cut tube DistanceToIn. Precision fixes for the cut tubes.
+- Fixed cone shape tester. Consistent Inside tolerance check for cone.
+- Fixed approaching bbox and removed leftover debug checks.
+- Fixed convention checking for points on edges.
+- Fixed formula for MakePlus/MinuTolerantSquare.
+- More precision fixes in the ShapeTester, improved use of ApproachSolid
+- Tolerance-related fixes for Tet to pass the ShapeTester.
+
+
+- Introduce standalone kFarAway in ShapeTester and further reduce errors
+- Fix types in ConventionChecker and ApproxEqual
+
+- Unify construction of Quadrilaterals
+- Consistently mark functions & destructors with VECCORE_ATT_HOST_DEVICE
+
+- Use VECGEOM_CONST for constants of kTolerance
+- Rename kernel function for ConstructManyOnGpu
+- UnplacedTrd now exposes interface 'normal' on GPU. Value is computed even if point not exactly on surface..
+- Transformation3D::Set copies from memory only if has_trans or has_rot are set.
+- Added new method to approach the next daughter to its bounding box.
+- Fixes for navigation index table in single precision mode.
+- Add GPU transfer testing for cms2018 geometry
+- (upstream/fix_bboxes_gpu) Fix bounding boxes on GPU by copying the values calculated on host.
+- Fix per solid tolerances and change types from double to Precision
+
+
+- Add bounding box tests to ctest
+- Implement ApproachSolid utility and store BBox in UnplacedVolume
+- Fix Polyhedron::Contains() around phi tolerance
+- Remove unused ScalarSegmentContainsKernel
+- Activate bulk copying of PlacedVolume and Transformation3D to the GPU.
+- Implement functions for bulk copying Transformation3D instances to GPU.
+- Implement functions for bulk-copying of placed volumes to GPU.
+- Improve const correctness in CudaManager.
+- Move verbose printouts to cout instead of cerr in CudaManager.
+- Add "Testing/" to .gitignore, which is created by ctest.
+
+- Reduce HOST_DEVICE attribute for PlacedVolume::GetName().
+- Add geometry test that compares CPU and GPU geometries.
+- [CMake] Remove -G from RelWithDebInfo.
+- [clang-format] Prevent clang-format from breaking doxygen comments.
+- [git] Extend git ignore to "build*" instead of build.
+- [NFC] Make detailed BVH description show up in doxygen.
+- [NFC] Add docstring for VECGEOM_DEVICE_INST_PLACED_VOLUME.
+- Fix failing test by creating new utility function
+- Fix failing tests in single precision
+- Change kTolerance to kHalfTolerance
+- [VECGEOM-578] Fix an FP comparison that was failing on the GPU
+- Fix a typo in debugging printout
+- [CMake] Fix doxygen search path.
+- Fix TestOrb unit test when using single precision

VecGeom/navigation/HybridNavigator2.h

@@ -242,35 +242,39 @@ public:
                                  [&](HybridManager2::BoxIdDistancePair_t hitbox) {
                                    // only consider those hitboxes which are within potential reach of this step
                                    if (!(step < hitbox.second)) {
+                                     vecgeom::Vector3D<double>
+                                         normal; // To reuse in printing below - else move it into 'if'
                                      VPlacedVolume const *candidate = LookupDaughter(lvol, hitbox.first);
                                      if (candidate == blocked) {
-                                       return false; // return early and go on in the looper
+                                       // return false; // return early and go on in the looper
+                                       candidate->Normal(localpoint, normal);
+                                       if (normal.Dot(localdir) >= 0.0) {
+                                         std::cerr << "HybridNav2> blocked " << candidate
+                                                   << " has normal.dir = " << normal.Dot(localdir) << " and distToIn = "
+                                                   << candidate->DistanceToIn(localpoint, localdir, step) << "\n";
+                                       }
                                      }
                                      const Precision ddistance = candidate->DistanceToIn(localpoint, localdir, step);
-                                     const auto valid          = !IsInf(ddistance) && ddistance < step;
+                                     const auto valid          = !IsInf(ddistance) && ddistance < step &&
+                                                        !((ddistance <= 0.) &&
+                                                          blocked == candidate); // && normal.Dot(localdir) > 0.0);
                                      hitcandidate              = valid ? candidate : hitcandidate;
                                      step                      = valid ? ddistance : step;
 #if 0 // enable for debugging
         if ( ddistance<=0 ) {
-           std::cerr << "negative distance found for " << candidate->GetName() << "\n"; 
+           std::cerr << "HybridNav2> negative distance found for " << candidate->GetName() << "\n"; 
            auto inside = candidate->Inside(localpoint);
-           if (inside == kSurface) {
-             std::cerr << "on surface\n";
-           }
-           if (inside == kOutside) {
-             std::cerr << "outside\n";
-           }
-           if (inside == kInside) {
-             std::cerr << "inside\n";
-           }
+           static std::string InsideCode[4] = { "N/A", "Inside", "Surface", "Outside" } ;
+           std::cerr << InsideCode[inside];
            const auto transf = candidate->GetTransformation();
            const auto unpl = candidate->GetUnplacedVolume();
-           Vector3D<Precision> normal;
+           Vector3D<Precision> normalDg;
            const auto testdaughterlocal = transf->Transform(localpoint);      
-           auto valid = unpl->Normal(testdaughterlocal, normal);
+           auto valid = unpl->Normal(testdaughterlocal, normalDg);
            
            const auto directiondaughterlocal = transf->TransformDirection(localdir);
-           const auto dot = normal.Dot(directiondaughterlocal);
+           const auto dot = normalDg.Dot(directiondaughterlocal);
+           std::cerr << " normal.dir = " << dot;
            if (dot >= 0) {
              std::cerr << " exiting " << valid << "\n";
            }

VecGeom/navigation/NewSimpleNavigator.h

@@ -55,6 +55,10 @@ public:
         step         = valid ? ddistance : step;
 #ifdef CHECKCONTAINS
       } else {
+        std::cerr << " INDA "
+                  << " contained in daughter " << daughter << " - inside = " << daughter->Inside(localpoint)
+                  << " , distToIn(p,v,s) = " << daughter->DistanceToIn(localpoint, localdir, step) << " \n";
+
         std::cerr << " INDA ";
         step         = -1.;
         hitcandidate = daughter;
@@ -67,10 +71,54 @@ public:
 
   VECGEOM_FORCE_INLINE
   VECCORE_ATT_HOST_DEVICE
-  virtual bool CheckDaughterIntersections(LogicalVolume const * /*lvol*/, Vector3D<Precision> const & /*localpoint*/,
-                                          Vector3D<Precision> const & /*localdir*/, VPlacedVolume const * /*blocked*/,
-                                          Precision & /*step*/, VPlacedVolume const *& /*hitcandidate*/) const override
+  virtual bool CheckDaughterIntersections(LogicalVolume const *lvol, Vector3D<Precision> const &localpoint,
+                                          Vector3D<Precision> const &localdir, VPlacedVolume const *blocked,
+                                          Precision &step, VPlacedVolume const *&hitcandidate) const override
   {
+    //  New Implementation JA 2021.03.18
+    static const double kMinExitingCos = 1.e-3;
+    VPlacedVolume const *excludedVol   = nullptr;
+    if (blocked) {
+      UVector3 normal;
+      blocked->Normal(localpoint, normal);
+      if (normal.Dot(localdir) >= kMinExitingCos) {
+        excludedVol = blocked;
+      }
+    }
+
+    // iterate over all daughters
+    auto *daughters = lvol->GetDaughtersp();
+    auto ndaughters = daughters->size();
+    for (decltype(ndaughters) d = 0; d < ndaughters; ++d) {
+      auto daughter = daughters->operator[](d);
+//    previous distance becomes step estimate, distance to daughter returned in workspace
+// SW: this makes the navigation more robust and it appears that I have to
+// put this at the moment since not all shapes respond yet with a negative distance if
+// the point is actually inside the daughter
+#ifdef CHECKCONTAINS
+      bool contains = daughter->Contains(localpoint);
+      if (!contains) {
+#endif
+        if (daughter != excludedVol) {
+          Precision ddistance = daughter->DistanceToIn(localpoint, localdir, step);
+
+          // if distance is negative; we are inside that daughter and should relocate
+          // unless distance is minus infinity
+          const bool valid = (ddistance < step && !IsInf(ddistance));
+          hitcandidate     = valid ? daughter : hitcandidate;
+          step             = valid ? ddistance : step;
+        }
+#ifdef CHECKCONTAINS
+      } else {
+        std::cerr << " INDA: contained in daughter " << daughter << " - inside = " << daughter->Inside(localpoint)
+                  << " , distToIn(p,v,s) = " << daughter->DistanceToIn(localpoint, localdir, step) << " \n";
+        step         = -1.;
+        hitcandidate = daughter;
+        break;
+      }
+#endif
+    }
+    // assert(false); --- Was not implemented before
     return false;
   }
 

VecGeom/navigation/SimpleABBoxNavigator.h

@@ -148,6 +148,75 @@ public:
     return false;
   }
 
+  virtual bool CheckDaughterIntersections(LogicalVolume const *lvol, Vector3D<Precision> const &localpoint,
+                                          Vector3D<Precision> const &localdir, VPlacedVolume const *blocked,
+                                          Precision &step, VPlacedVolume const *&hitcandidate) const override
+  {
+    static size_t counter = 0;
+    counter++;
+    // The following construct reserves stackspace for objects
+    // of type IdDistPair_t WITHOUT initializing those objects
+    using IdDistPair_t = ABBoxManager::BoxIdDistancePair_t;
+    char stackspace[VECGEOM_MAXDAUGHTERS * sizeof(IdDistPair_t)];
+    IdDistPair_t *hitlist = reinterpret_cast<IdDistPair_t *>(&stackspace);
+
+    if (lvol->GetDaughtersp()->size() == 0) return false;
+
+    int size;
+    ABBoxManager::ABBoxContainer_v bboxes = fABBoxManager.GetABBoxes_v(lvol, size);
+    auto ncandidates                      = GetHitCandidates_v(lvol, localpoint, localdir, bboxes, size, hitlist);
+
+    // sort candidates according to their bounding volume hit distance
+    insertionsort(hitlist, ncandidates);
+
+    for (size_t index = 0; index < ncandidates; ++index) {
+      auto &hitbox                   = hitlist[index];
+      VPlacedVolume const *candidate = LookupDaughter(lvol, hitbox.first);
+
+      // only consider those hitboxes which are within potential reach of this step
+      if (hitbox.second <= step) { // !(step < hitbox.second)) {
+#ifdef VERBOSE
+        std::cerr << "SimpleABBoxNavigator: checking id " << hitbox.first << " at box distance " << hitbox.second
+                  << "\n";
+#endif
+        //           if( hitbox.second < 0 ){
+        //            //   std::cerr << "hit dist < 0. \n";
+        //              bool checkindaughter = candidate->Contains( localpoint );
+        //              if( checkindaughter == true ){
+        //                  std::cerr << " -- unexpected : in daughter vol \n";
+        //
+        //                  // need to relocate
+        //                  step = 0;
+        //                  hitcandidate = candidate;
+        //                  // THE ALTERNATIVE WOULD BE TO PUSH THE CURRENT STATE AND RETURN DIRECTLY
+        //                  break;
+        //              }
+        //          }
+        Precision ddistance = candidate->DistanceToIn(localpoint, localdir, step);
+        vecgeom::Vector3D<double> normal; // To reuse in printing below - else move it into 'if'
+#ifdef VERBOSE
+        std::cerr << " . Distance to " << candidate->GetLabel() << " is " << ddistance;
+#endif
+        if (ddistance <= 0. /* && candidate == blocked */) {
+          candidate->Normal(localpoint, normal);
+#ifdef VERBOSE
+          std::cerr << "  SimpleABBoxNavigator:  Negative daughter dist = " << ddistance
+                    << (blocked == candidate ? " Is " : " Not ") << "Blocked. "
+                    << " normal= " << normal << "  normal.dir= " << normal.Dot(localdir) << "\n";
+#endif
+        }
+        const auto valid = !IsInf(ddistance) && ddistance < step &&
+                           !((ddistance <= 0.) && (blocked == candidate || normal.Dot(localdir) > 0.0));
+
+        hitcandidate = valid ? candidate : hitcandidate;
+        step         = valid ? ddistance : step;
+      } else {
+        break;
+      }
+    }
+    return false;
+  }
+
   static VNavigator *Instance()
   {
     static SimpleABBoxNavigator instance;

VecGeom/navigation/VNavigator.h

@@ -120,6 +120,7 @@ public:
                                           Vector3D<Precision> const & /*localdir*/, VPlacedVolume const * /*blocked*/,
                                           Precision & /*step*/, VPlacedVolume const *& /*hitcandidate*/) const
   {
+    assert(false); // Not implemented --- notify of failure !!
     return false;
   }
 

VecGeom/navigation/VoxelSafetyEstimator.h

@@ -103,7 +103,15 @@ public:
       const int *safetycandidates{nullptr};
 
       int length{0};
-      safetycandidates = structure->fVoxelToCandidate->getProperties(lp, length);
+      auto voxelHashMap = structure->fVoxelToCandidate;
+      if (voxelHashMap == nullptr) {
+        std::cerr << "ERROR> VoxelSafetyEstimator::ComputeSafetyForLocalPoint call# " << counter
+                  << " no structure VoxelToCandidate structure found for volume "
+                  << " phys: " << pvol->GetName() << "  physvol id = " << pvol->id() << " log : " << lvol->GetName()
+                  << "  log vol id = " << lvol->id() << std::endl;
+        return 0.0;
+      }
+      safetycandidates = voxelHashMap->getProperties(lp, length);
       if (length > 0) {
         const bool needmother        = true; //(safetycandidates[0] == -1);
         const Precision safetymother = needmother ? lvol->GetUnplacedVolume()->SafetyToOut(localpoint) : 0.;
@@ -146,14 +154,16 @@ public:
         }
 #endif
 
-        if (returnvalue < 0) {
-          std::cerr << "returning negative value " << safetyToDaughters << " " << safetymother << "\n";
+        if (returnvalue < -1.0e-10) {
+          std::cerr << "VoxelSafetyEstimator: returning negative value.  saf-to-daughters= " << safetyToDaughters
+                    << " saf-to-mother = " << safetymother << "\n";
         }
 
         return returnvalue;
       } else {
         // no information for this voxel present
-        std::cerr << counter << " no information for this voxel present " << localpoint << " at key "
+        std::cerr << "Warning> ComputeSafetyForLocalPoint call# " << counter
+                  << " no information for this voxel present " << localpoint << " at key "
                   << structure->fVoxelToCandidate->getKey(lp.x(), lp.y(), lp.z()) << " \n ";
         return 0.;
       }

VecGeom/volumes/utilities/VolumeUtilities.h

@@ -17,6 +17,9 @@
 #include "VecGeom/navigation/VNavigator.h"
 #include "VecGeom/navigation/GlobalLocator.h"
 #include "VecGeom/management/GeoManager.h"
+
+#include "VecGeom/management/GeoManager.h"
+
 #include <cstdio>
 #include <random>
 #include <vector>
@@ -309,19 +312,18 @@ Precision UncontainedCapacity(VPlacedVolume const &volume)
  * @param points is the output container, provided by the caller.
  */
 template <typename TrackContainer>
-VECGEOM_FORCE_INLINE
-void FillUncontainedPoints(VPlacedVolume const &volume, TrackContainer &points)
+VECGEOM_FORCE_INLINE bool FillUncontainedPoints(VPlacedVolume const &volume, TrackContainer &points)
 {
   static double lastUncontCap = 0.0;
   double uncontainedCapacity  = UncontainedCapacity(volume);
   if (uncontainedCapacity != lastUncontCap) {
-    printf("Uncontained capacity for %s: %g units\n", volume.GetLabel().c_str(), uncontainedCapacity);
+    std::cout << "Uncontained capacity for " << volume.GetLabel() << ":" << uncontainedCapacity << " units\n";
     lastUncontCap = uncontainedCapacity;
   }
   if (uncontainedCapacity <= 1000 * kTolerance) {
-    std::cout << "\nVolUtil: FillUncontPts: ERROR: Volume provided <" << volume.GetLabel()
-              << "> does not have uncontained capacity!  Aborting.\n";
-    assert(false);
+    std::cout << "\nVolUtil: FillUncontPts: WARNING: Volume provided <" << volume.GetLabel()
+              << "> does not have uncontained capacity!  Method returns false.\n";
+    return false;
   }
 
   const int size = points.capacity();
@@ -332,18 +334,25 @@ void FillUncontainedPoints(VPlacedVolume const &volume, TrackContainer &points)
   offset                        = 0.5 * (upper + lower);
   const Vector3D<Precision> dim = 0.5 * (upper - lower);
 
-  int tries = 0;
+  int totaltries = 0;
   for (int i = 0; i < size; ++i) {
     bool contained;
     Vector3D<Precision> point;
-    tries = 0;
+    totaltries = 0;
     do {
       // ensure that point is contained in mother volume
+      int onego = 0;
       do {
-        ++tries;
-        if (tries % 1000000 == 0) {
+        ++totaltries;
+        ++onego;
+        if (totaltries % 10000 == 0) {
           printf("%s line %i: Warning: %i tries to find uncontained points... volume=%s.  Please check.\n", __FILE__,
-                 __LINE__, tries, volume.GetLabel().c_str());
+                 __LINE__, totaltries, volume.GetLabel().c_str());
+        }
+        if (totaltries % 5000000 == 0) {
+          double ratio = 1.0 * i / totaltries;
+          printf("Progress : %i tries ( succeeded = %i , ratio %f %% ) to find uncontained points... volume=%s.\n",
+                 totaltries, i, 100. * ratio, volume.GetLabel().c_str());
         }
 
         point = offset + SamplePoint(dim);
@@ -361,15 +370,17 @@ void FillUncontainedPoints(VPlacedVolume const &volume, TrackContainer &points)
       }
     } while (contained);
   }
+  return true;
 }
 
 template <typename TrackContainer>
-VECGEOM_FORCE_INLINE
-void FillUncontainedPoints(LogicalVolume const &volume, TrackContainer &points)
+VECGEOM_FORCE_INLINE bool FillUncontainedPoints(LogicalVolume const &volume, TrackContainer &points)
 {
   VPlacedVolume const *const placed = volume.Place();
-  FillUncontainedPoints(*placed, points);
+  bool good                         = FillUncontainedPoints(*placed, points);
   delete placed;
+
+  return good;
 }
 
 // *** The following functions allow to give an external generator
@@ -384,8 +395,8 @@ void FillUncontainedPoints(LogicalVolume const &volume, TrackContainer &points)
  * @param points is the output container, provided by the caller.
  */
 template <typename RandomEngine, typename TrackContainer>
-VECGEOM_FORCE_INLINE
-void FillUncontainedPoints(VPlacedVolume const &volume, RandomEngine &rngengine, TrackContainer &points)
+VECGEOM_FORCE_INLINE bool FillUncontainedPoints(VPlacedVolume const &volume, RandomEngine &rngengine,
+                                                TrackContainer &points)
 {
   static double lastUncontCap = 0.0;
   double uncontainedCapacity  = UncontainedCapacity(volume);
@@ -393,10 +404,22 @@ void FillUncontainedPoints(VPlacedVolume const &volume, RandomEngine &rngengine,
     printf("Uncontained capacity for %s: %g units\n", volume.GetLabel().c_str(), uncontainedCapacity);
     lastUncontCap = uncontainedCapacity;
   }
+  double totalcapacity = const_cast<VPlacedVolume &>(volume).Capacity();
+
+  std::cout << "\nVolUtil: FillUncontPts: Volume <" << volume.GetLabel() << "  capacities: total =  " << totalcapacity
+            << " uncontained = " << uncontainedCapacity << "\n";
+
   if (uncontainedCapacity <= 1000 * kTolerance) {
+    // double checkUC= UncontainedCapacity(volume); // Rerun - for debugging ...
     std::cout << "\nVolUtil: FillUncontPts: ERROR: Volume provided <" << volume.GetLabel()
-              << "> does not have uncontained capacity!  Aborting.\n";
-    assert(false);
+              << "> does not have uncontained capacity!  "
+              << "    Value = " << uncontainedCapacity << " \n"
+              << "      contained = " << totalcapacity
+        // << "    check = " << checkUC << " \n"
+        ;
+    // if( checkUC < 0 ) { assert(false); }
+    return false;
+    // TODO --- try to find points anyway, and decide if real points were found
   }
 
   const int size = points.capacity();
@@ -407,18 +430,28 @@ void FillUncontainedPoints(VPlacedVolume const &volume, RandomEngine &rngengine,
   offset                        = 0.5 * (upper + lower);
   const Vector3D<Precision> dim = 0.5 * (upper - lower);
 
-  int tries = 0;
-  for (int i = 0; i < size; ++i) {
+  const int maxtries = 100 * 1000 * 1000;
+
+  int tries = 0; // count total trials ...
+  int i;
+  for (i = 0; i < size; ++i) {
     bool contained;
     Vector3D<Precision> point;
-    tries = 0;
     do {
       // ensure that point is contained in mother volume
+      int onego = 0;
       do {
         ++tries;
-        if (tries % 1000000 == 0) {
-          printf("%s line %i: Warning: %i tries to find uncontained points... volume=%s.  Please check.\n", __FILE__,
-                 __LINE__, tries, volume.GetLabel().c_str());
+        onego++;
+        if (onego % 100000 == 0) {
+          printf("%s line %i: Warning: %i tries ( success = %i ) to find uncontained points... volume=%s.  Please "
+                 "check.\n",
+                 __FILE__, __LINE__, tries, i, volume.GetLabel().c_str());
+        }
+        if (tries % 5000000 == 0) {
+          double ratio = 1.0 * i / tries;
+          printf("Progress : %i tries ( succeeded = %i , ratio %f %% ) to find uncontained points... volume=%s.\n",
+                 tries, i, 100.0 * ratio, volume.GetLabel().c_str());
         }
 
         point = offset + SamplePoint(dim, rngengine);
@@ -434,17 +467,24 @@ void FillUncontainedPoints(VPlacedVolume const &volume, RandomEngine &rngengine,
           break;
         }
       }
-    } while (contained);
+    } while (contained && tries < maxtries);
+
+    if (tries >= maxtries) break;
   }
+  std::cout << " FillUncontained:  trials " << tries << " for num points = " << i << " ( out of " << size
+            << " requested - "  << " success ratio = " << (i * 1.0) / tries << "\n";
+  return (i > 0);
 }
 
 template <typename RandomEngine, typename TrackContainer>
-VECGEOM_FORCE_INLINE
-void FillUncontainedPoints(LogicalVolume const &volume, RandomEngine &rngengine, TrackContainer &points)
+VECGEOM_FORCE_INLINE bool FillUncontainedPoints(LogicalVolume const &volume, RandomEngine &rngengine,
+                                                TrackContainer &points)
 {
   VPlacedVolume const *const placed = volume.Place();
-  FillUncontainedPoints(*placed, rngengine, points);
+  bool good                         = FillUncontainedPoints(*placed, rngengine, points);
   delete placed;
+
+  return good;
 }
 
 /**
@@ -571,16 +611,16 @@ void FillContainedPoints(VPlacedVolume const &volume, const double bias, TrackCo
 
   // remove contained points to reduce bias as needed
   int i     = 0;
-  int tries = 0;
+  int totaltries = 0;
   while (static_cast<double>(insideCount) / static_cast<double>(size) > bias) {
     while (!insideVector[i])
       ++i;
     bool contained;
     do {
-      ++tries;
-      if (tries % 1000000 == 0) {
-        printf("%s line %i: Warning: %i tries to reduce bias... volume=%s.  Please check.\n", __FILE__, __LINE__, tries,
-               volume.GetLabel().c_str());
+      ++totaltries;
+      if (totaltries % 1000000 == 0) {
+        printf("%s line %i: Warning: %i totaltries to reduce bias... volume=%s.  Please check.\n", __FILE__, __LINE__,
+               totaltries, volume.GetLabel().c_str());
       }
 
       points.set(i, offset + SamplePoint(dim));
@@ -595,11 +635,12 @@ void FillContainedPoints(VPlacedVolume const &volume, const double bias, TrackCo
       }
     } while (contained);
     insideVector[i] = false;
-    tries           = 0;
+    // tries           = 0;
     --insideCount;
     ++i;
   }
 
+  int tries;
   // add contained points to increase bias as needed
   i     = 0;
   tries = 0;
@@ -747,7 +788,12 @@ inline void FillGlobalPointsAndDirectionsForLogicalVolume(LogicalVolume const *l
     VPlacedVolume const *pvol = allpaths.front()->Top();
 
     // generate points which are in lvol but not in its daughters
-    FillUncontainedPoints(*pvol, localpoints);
+    bool good = FillUncontainedPoints(*pvol, localpoints);
+    // assert(good);
+    if (!good) {
+      std::cerr << "FATAL ERROR> FillUncontainedPoints failed for volume " << pvol->GetName() << std::endl;
+      exit(1);
+    }
 
     // now have the points in the local reference frame of the logical volume
     FillBiasedDirections(*lvol, localpoints, fraction, directions);
@@ -857,7 +903,12 @@ inline void FillGlobalPointsForLogicalVolume(LogicalVolume const *lvol, TrackCon
       FillContainedPoints(*pvol, localpoints);
     } else {
       // generate points which are in lvol but not in its daughters
-      FillUncontainedPoints(*pvol, localpoints);
+      bool good = FillUncontainedPoints(*pvol, localpoints);
+      // assert(good);