From 6c36a3dacfe4f7083c2fc9c40331ad004d82caa5 Mon Sep 17 00:00:00 2001 From: kjplows Date: Thu, 27 Nov 2025 12:09:27 -0600 Subject: [PATCH 1/4] Fixes two segfaults in CAF processing of SBND v10_06_00_09. First fix: Added guard in G4 weight calculator against trajectory points outside world volume. Second fix: Added guard in RecoUtils.cc for when a trajectory does not intersect detector walls. --- sbncode/CAFMaker/RecoUtils/RecoUtils.cc | 22 ++++++++++++++++++- .../Calculators/Geant4/Geant4WeightCalc.cxx | 3 +++ 2 files changed, 24 insertions(+), 1 deletion(-) diff --git a/sbncode/CAFMaker/RecoUtils/RecoUtils.cc b/sbncode/CAFMaker/RecoUtils/RecoUtils.cc index a2cdc8922..edb7ab358 100644 --- a/sbncode/CAFMaker/RecoUtils/RecoUtils.cc +++ b/sbncode/CAFMaker/RecoUtils/RecoUtils.cc @@ -143,7 +143,27 @@ caf::Wall_t sbn::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); + //assert(intersections.size() == 2); + /* + * There are either infinity, 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 + */ + switch( intersections.size() ) { + case 0: // Set kWallNone and return + return caf::kWallNone; + break; + case 1: // Set the second intersection to be the same as the first, and fall through + intersections.emplace_back(intersections.front()); + [[fallthrough]]; + case 2: + break; + default: // There are infinity points to consider. Just use the first two. + intersections.resize(2); + break; + } // get the intersection point closer to p0 int intersection_i = ((intersections[0] - p0).Mag() < (intersections[1] - p0).Mag()) ? 0 : 1; diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx index ec9d8cc6e..e7e76ecd4 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 || testmaterial1 == nullptr ) { + 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" ) ){ From 239646ec7687f45d4d7c7f6f895ab609cc8644d1 Mon Sep 17 00:00:00 2001 From: kjplows Date: Wed, 3 Dec 2025 06:45:27 -0600 Subject: [PATCH 2/4] Compute closest intersection in a loop --- sbncode/CAFMaker/RecoUtils/RecoUtils.cc | 49 ++++++++++--------------- 1 file changed, 20 insertions(+), 29 deletions(-) diff --git a/sbncode/CAFMaker/RecoUtils/RecoUtils.cc b/sbncode/CAFMaker/RecoUtils/RecoUtils.cc index edb7ab358..a67e7d0fe 100644 --- a/sbncode/CAFMaker/RecoUtils/RecoUtils.cc +++ b/sbncode/CAFMaker/RecoUtils/RecoUtils.cc @@ -143,58 +143,49 @@ caf::Wall_t sbn::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 infinity, two, or one, or zero intersection points. + * 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 */ - switch( intersections.size() ) { - case 0: // Set kWallNone and return - return caf::kWallNone; - break; - case 1: // Set the second intersection to be the same as the first, and fall through - intersections.emplace_back(intersections.front()); - [[fallthrough]]; - case 2: - break; - default: // There are infinity points to consider. Just use the first two. - intersections.resize(2); - break; - } + + if( intersections.empty() ) return caf::kWallNone; + // ensure intersections has two points. No op if already 2 + intersections.resize( 2, intersections.front() ); // 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 From a70eb00214d3707e07c63f6032dc63e0e03672d6 Mon Sep 17 00:00:00 2001 From: kjplows Date: Wed, 3 Dec 2025 06:46:52 -0600 Subject: [PATCH 3/4] Remove duplicate check on testmaterial1 --- sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx index e7e76ecd4..5f45db672 100644 --- a/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx +++ b/sbncode/SBNEventWeight/Calculators/Geant4/Geant4WeightCalc.cxx @@ -297,7 +297,7 @@ 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 || testmaterial1 == nullptr ) { + 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 From ebafba001e3f0388b78c8d0242de35858db16471 Mon Sep 17 00:00:00 2001 From: kjplows Date: Wed, 3 Dec 2025 10:36:46 -0600 Subject: [PATCH 4/4] no resizing needed --- sbncode/CAFMaker/RecoUtils/RecoUtils.cc | 2 -- 1 file changed, 2 deletions(-) diff --git a/sbncode/CAFMaker/RecoUtils/RecoUtils.cc b/sbncode/CAFMaker/RecoUtils/RecoUtils.cc index a67e7d0fe..0c70d25df 100644 --- a/sbncode/CAFMaker/RecoUtils/RecoUtils.cc +++ b/sbncode/CAFMaker/RecoUtils/RecoUtils.cc @@ -152,8 +152,6 @@ caf::Wall_t sbn::GetWallCross(const geo::BoxBoundedGeo &volume, const TVector3 p */ if( intersections.empty() ) return caf::kWallNone; - // ensure intersections has two points. No op if already 2 - intersections.resize( 2, intersections.front() ); // get the intersection point closer to p0 TVector3 closestIntersection;