Merge pull request #609 from SBNSoftware/bugfix/kplows_fix_springProcessingCAFs

Fixes two segfaults in CAF processing of SBND v10_06_00_09.

Author: kjplows

sbncode/CAFMaker/RecoUtils/RecoUtils.cc

@@ -143,38 +143,47 @@ caf::Wall_t sbn::GetWallCross(const geo::BoxBoundedGeo &volume, const TVector3 p
   TVector3 direction = (p1 - p0) * ( 1. / (p1 - p0).Mag());
   std::vector<TVector3> intersections = volume.GetIntersections(p0, direction);
 
-  assert(intersections.size() == 2);
+  /*
+   * There are either two, or one, or zero intersection points.
+   * As per larcorealg/Geometry/BoxBoundedGeo.h documentation:
+   * If the return std::vector is empty the trajectory does not intersect with the box.
+   * Normally the return value should have one (if the trajectory originates in the box) or two (else) entries.
+   * If the return value has two entries the first represents the entry point and the second the exit point
+   */
+
+  if( intersections.empty() ) return caf::kWallNone;
 
   // get the intersection point closer to p0
-  int intersection_i = ((intersections[0] - p0).Mag() < (intersections[1] - p0).Mag()) ? 0 : 1;
+  TVector3 closestIntersection;
+  double minDistance2 = std::numeric_limits<double>::max();
+
+  for( TVector3 const & point : intersections ) {
+    const double d2 = (point - p0).Mag2();
+    if( d2 > minDistance2 ) continue;
+    minDistance2 = d2;
+    closestIntersection = point;
+  }
 
   double eps = 1e-3;
-  if (abs(intersections[intersection_i].X() - volume.MinX()) < eps) {
-    //std::cout << "Left\n";
+  if (abs(closestIntersection.X() - volume.MinX()) < eps) {
     return caf::kWallLeft;
   }
-  else if (abs(intersections[intersection_i].X() - volume.MaxX()) < eps) {
-    //std::cout << "Right\n";
+  else if (abs(closestIntersection.X() - volume.MaxX()) < eps) {
     return caf::kWallRight;
   }
-  else if (abs(intersections[intersection_i].Y() - volume.MinY()) < eps) {
-    //std::cout << "Bottom\n";
+  else if (abs(closestIntersection.Y() - volume.MinY()) < eps) {
     return caf::kWallBottom;
   }
-  else if (abs(intersections[intersection_i].Y() - volume.MaxY()) < eps) {
-    //std::cout << "Top\n";
+  else if (abs(closestIntersection.Y() - volume.MaxY()) < eps) {
     return caf::kWallTop;
   }
-  else if (abs(intersections[intersection_i].Z() - volume.MinZ()) < eps) {
-    //std::cout << "Front\n";
+  else if (abs(closestIntersection.Z() - volume.MinZ()) < eps) {
     return caf::kWallFront;
   }
-  else if (abs(intersections[intersection_i].Z() - volume.MaxZ()) < eps) {
-    //std::cout << "Back\n";
+  else if (abs(closestIntersection.Z() - volume.MaxZ()) < eps) {
     return caf::kWallBack;
   }
   else assert(false);
-  //std::cout << "None\n";
 
   return caf::kWallNone;
 }//GetWallCross

sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx

@@ -297,6 +297,9 @@ std::vector<float> Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) {
         double Z = p.Position(i).Z();
         geo::Point_t testpoint1 { X, Y, Z };
         const TGeoMaterial* testmaterial1 = fGeometryService->Material( testpoint1 );
+	if( ! testmaterial1 ) {  
+	  break; // The trajectory point is outside the world, or the geometry service is unable to handle the point
+	}
         //For now, just going to reweight the points within the LAr of the TPC
         // TODO check if this is right
         if ( !strcmp( testmaterial1->GetName(), "LAr" ) ){