diff --git a/sbncode/CAFMaker/FillTrue.cxx b/sbncode/CAFMaker/FillTrue.cxx index bbc936e88..ef1b4cc0c 100644 --- a/sbncode/CAFMaker/FillTrue.cxx +++ b/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 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::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 diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx index ec9d8cc6e..896eaa2ef 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx +++ b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx @@ -297,6 +297,9 @@ std::vector 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" ) ){