Merge branch 'develop' into feature/mdeltutt_fluxweight_g4bnb

Author: kjplows

sbncode/CAFMaker/FillTrue.cxx

@@ -1260,38 +1260,47 @@ caf::Wall_t caf::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/LArG4/CMakeLists.txt

@@ -1,15 +1,5 @@
 add_subdirectory(PhysicsLists)
 
-cet_build_plugin(G4InfoReducer art::module
-                       LIBRARIES
-                           larsim::Simulation
-                           larcore::Geometry_Geometry_service
-                           lardata::DetectorInfoServices_DetectorClocksServiceStandard_service
-                           art::Persistency_Common
-                           art::Persistency_Provenance
-                           art::Utilities canvas::canvas
-                           BASENAME_ONLY)
-
 set(
   MODULE_LIBRARIES
     art::Utilities
@@ -61,14 +51,5 @@ cet_build_plugin(MergeSimSourcesSBN art::EDProducer
   BASENAME_ONLY
 )
 
-cet_build_plugin(SimpleMerge art::EDProducer
-  LIBRARIES
-  lardataobj::MCBase
-  larcorealg::headers
-  larcorealg::CoreUtils
-  messagefacility::MF_MessageLogger
-)
-
 install_headers()
-install_fhicl()
 install_source()

sbncode/LArG4/g4inforeducer.fcl

@@ -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" ) ){