Skip to content

Fix bugs in v10 upgrade. #501

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ibsafa merged 4 commits into from
Feb 14, 2025
Merged

Fix bugs in v10 upgrade. #501

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
46db6fb
Fix bugs in v10 upgrade.
gputnam Feb 2, 2025
8d355d0
Use correct half-length calculation.
gputnam Feb 3, 2025
05734b6
Merge branch 'develop' into feature/gputnam-fix-v10
ibsafa Feb 14, 2025
b7a6e6b
Merge branch 'develop' into feature/gputnam-fix-v10
ibsafa Feb 14, 2025
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
2 changes: 1 addition & 1 deletion sbncode/CAFMaker/FillTrue.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -300,7 +300,7 @@ namespace caf {
double cosgamma = std::abs(std::sin(angletovert)*dir.Y() + std::cos(angletovert)*dir.Z());
double pitch;
if (cosgamma) {
pitch = plane.View()/cosgamma;
pitch = plane.WirePitch()/cosgamma;
}
else {
pitch = 0.;
Expand Down
4 changes: 2 additions & 2 deletions sbncode/Cluster3D/SnippetHit3DBuilderSBN_tool.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1633,8 +1633,8 @@ geo::WireID SnippetHit3DBuilderSBN::NearestWireID(const Eigen::Vector3f& positio
mf::LogWarning("Cluster3D") << "Exception caught finding nearest wire, position - " << exc.what() << std::endl;

// Assume extremum for wire number depending on z coordinate
if (position[2] < 0.5 * m_geometry->TPC({0,0}).Length()) wire = 0;
else wire = m_wireReadoutAlg->Nwires(wireIDIn.asPlaneID()) - 1;
if (position[2] < m_geometry->TPC({0,0}).ActiveLength()*0.5) wire = 0;
else wire = m_wireReadoutAlg->Nwires(wireIDIn.asPlaneID()) - 1;
}

geo::WireID wireID(wireIDIn.Cryostat,wireIDIn.TPC,wireIDIn.Plane,wire);
Expand Down
6 changes: 3 additions & 3 deletions sbncode/TPCReco/NuVertexChargeTree_module.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -949,9 +949,9 @@ void sbn::NuVertexChargeTree::FillNeutrino(const simb::MCTruth &nu,
geo::TPCGeo const* firstTPC = geo.PositionToTPCptr(start_loc);
if (firstTPC) {
geo::TPCID firstTPCID = firstTPC->ID();
fFSPPitch0.push_back(sbn::GetPitch(wireReadout, sce, start_loc, start_dir, wireReadout.Plane(geo::PlaneID(firstTPCID, 0)).View(), *firstTPC, true, true));
fFSPPitch1.push_back(sbn::GetPitch(wireReadout, sce, start_loc, start_dir, wireReadout.Plane(geo::PlaneID(firstTPCID, 1)).View(), *firstTPC, true, true));
fFSPPitch2.push_back(sbn::GetPitch(wireReadout, sce, start_loc, start_dir, wireReadout.Plane(geo::PlaneID(firstTPCID, 2)).View(), *firstTPC, true, true));
fFSPPitch0.push_back(sbn::GetPitch(wireReadout, sce, start_loc, start_dir, wireReadout.Plane(geo::PlaneID(firstTPCID, 0)).View(), geo::PlaneID(firstTPCID, 0), true, true));
fFSPPitch1.push_back(sbn::GetPitch(wireReadout, sce, start_loc, start_dir, wireReadout.Plane(geo::PlaneID(firstTPCID, 1)).View(), geo::PlaneID(firstTPCID, 1), true, true));
fFSPPitch2.push_back(sbn::GetPitch(wireReadout, sce, start_loc, start_dir, wireReadout.Plane(geo::PlaneID(firstTPCID, 2)).View(), geo::PlaneID(firstTPCID, 2), true, true));
Comment on lines +952 to +954
Copy link
Member

@PetrilloAtWork PetrilloAtWork Feb 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Was this an existing bug, or a new one?

Copy link
Contributor Author

@gputnam gputnam Feb 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I needed to change the GetPitch interface in order to fix introduced by the v10 transition.

}
else {
fFSPPitch0.push_back(-1);
Expand Down
2 changes: 1 addition & 1 deletion sbncode/TPCReco/VertexStub/StubBuilder.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -270,7 +270,7 @@ sbn::Stub sbn::StubBuilder::FromVertexHit(const art::Ptr<recob::Slice> &slice,

// See if we can compute a track pitch
if (pfp_ind >= 0 && trks[pfp_ind]) {
stub.trkpitch.push_back(sbn::GetPitch(wireReadout, sce, trks[pfp_ind]->Start(), trks[pfp_ind]->StartDirection(), wireReadout.Plane(vhit_hit.WireID()).View(), geometry.TPC(vhit_hit.WireID()), true, fPositionsAreSCECorrected));
stub.trkpitch.push_back(sbn::GetPitch(wireReadout, sce, trks[pfp_ind]->Start(), trks[pfp_ind]->StartDirection(), wireReadout.Plane(vhit_hit.WireID()).View(), vhit_hit.WireID(), true, fPositionsAreSCECorrected));
}
else {
stub.trkpitch.push_back(-1.);
Expand Down
20 changes: 13 additions & 7 deletions sbncode/TPCReco/VertexStub/StubMergeAlgorithms.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,9 @@
#include "StubMergeAlgorithms.h"

#include "larcorealg/Geometry/WireReadoutGeom.h"
#include "larcore/Geometry/Geometry.h"
#include "larcore/CoreUtils/ServiceUtil.h"
#include "larcorealg/Geometry/GeometryCore.h"

geo::Point_t sbn::GetLocation(const spacecharge::SpaceCharge *sce, geo::Point_t loc_w, geo::TPCID TPC, float xsign) {
if (sce && sce->EnableCalSpatialSCE()) {
Expand Down Expand Up @@ -49,14 +52,17 @@ geo::Point_t sbn::GetLocationAtWires(const spacecharge::SpaceCharge *sce, geo::P
double sbn::GetPitch(
const geo::WireReadoutGeom &wireReadout, const spacecharge::SpaceCharge *sce,
geo::Point_t loc, geo::Vector_t dir,
geo::View_t view, geo::TPCGeo const& tpc,
geo::View_t view, geo::PlaneID const& plane_id,
bool correct_sce, bool track_is_sce_corrected, float xsign) {

double angleToVert = wireReadout.WireAngleToVertical(view, tpc.ID()) - 0.5*::util::pi<>();
const geo::GeometryCore *geo = lar::providerFrom<geo::Geometry>();

double angleToVert = wireReadout.WireAngleToVertical(view, plane_id) - 0.5*::util::pi<>();
Copy link
Member

@PetrilloAtWork PetrilloAtWork Feb 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

By the way, this is the direction from geo->Plane(plane_id).GetIncreasingWireDirection(), which would make cosgamma a simple std::abs(geo::vect::Dot(dirToVert, dir_w)).

Copy link
Contributor Author

@gputnam gputnam Feb 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That sounds right. I prefer this method, because other pieces of code (like the Calorimetry) also use it to compute the pitch. So more people would notice if something breaks.


geo::Vector_t dir_w;

geo::PlaneGeo const& plane = wireReadout.Plane(geo::PlaneID{0, 0, view});
geo::PlaneGeo const& plane = wireReadout.Plane(plane_id);
geo::TPCGeo const& tpc = geo->TPC(plane_id);
// "dir_w" should be the direction that the wires see. If the track already has the field
// distortion corrections applied, then we need to de-apply them to get the direction as
// seen by the wire planes
Expand Down Expand Up @@ -92,8 +98,8 @@ double sbn::GetPitch(
loc_w = sbn::GetLocationAtWires(sce, loc, tpc.DriftDir(), xsign);
}

geo::Point_t locw_traj = (correct_sce) ? sbn::GetLocation(sce, loc_w, tpc.ID(), xsign) : loc_w;
geo::Point_t locw_pdx_traj = (correct_sce) ? sbn::GetLocation(sce, loc_w + pitch * dir_w, tpc.ID(), xsign) : (loc_w + pitch * dir_w);
geo::Point_t locw_traj = (correct_sce) ? sbn::GetLocation(sce, loc_w, plane_id, xsign) : loc_w;
geo::Point_t locw_pdx_traj = (correct_sce) ? sbn::GetLocation(sce, loc_w + pitch * dir_w, plane_id, xsign) : (loc_w + pitch * dir_w);

pitch = (locw_traj - locw_pdx_traj).R();

Expand Down Expand Up @@ -222,8 +228,8 @@ float sbn::StubPeakdQdxOffset(const sbn::StubInfo &A, const sbn::StubInfo &B,
geo::Vector_t dir((end-vtx).Unit());

// Input point, dir are always space charge corrected here
float Apitch = GetPitch(wireReadout, sce, end, dir, A.vhit_hit->View(), geometry.TPC(A.vhit_hit->WireID()), true, true);
float Bpitch = GetPitch(wireReadout, sce, end, dir, B.vhit_hit->View(), geometry.TPC(B.vhit_hit->WireID()), true, true);
float Apitch = GetPitch(wireReadout, sce, end, dir, A.vhit_hit->View(), A.vhit_hit->WireID(), true, true);
float Bpitch = GetPitch(wireReadout, sce, end, dir, B.vhit_hit->View(), B.vhit_hit->WireID(), true, true);

return abs(A.vhit->charge / Apitch - B.vhit->charge / Bpitch);
}
2 changes: 1 addition & 1 deletion sbncode/TPCReco/VertexStub/StubMergeAlgorithms.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ struct StubInfo {
double GetPitch(
const geo::WireReadoutGeom &wireReadout, const spacecharge::SpaceCharge *sce,
geo::Point_t loc, geo::Vector_t dir,
geo::View_t view, geo::TPCGeo const& tpc,
geo::View_t view, geo::PlaneID const& plane_id,
bool correct_sce, bool track_is_sce_corrected, float xsign=1.);

/// Get the location in the presence of space charge
Expand Down
2 changes: 1 addition & 1 deletion sbncode/TPCReco/VertexStub/VertexChargeVacuum_module.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -430,7 +430,7 @@ void sbn::VertexChargeVacuum::produce(art::Event& evt)

// Compute the pitch. Since we have used the corrected vertex and Space-Point position, the
// pt and dir here are space-charge corrected regardless of the input configuration
vhit.pitch = sbn::GetPitch(wireReadout, sce, spXYZ, dir, hit.View(), geo->TPC(hit.WireID()), fCorrectSCE, true);
vhit.pitch = sbn::GetPitch(wireReadout, sce, spXYZ, dir, hit.View(), hit.WireID(), fCorrectSCE, true);

vhit.dqdx = vhit.charge / vhit.pitch;

Expand Down