Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
39 changes: 24 additions & 15 deletions sbncode/CAFMaker/RecoUtils/RecoUtils.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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" ) ){
Expand Down