From ba72c99d06f0bf5fb56b353e82f68bc9ea636e0e Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Wed, 6 Aug 2025 19:49:43 -0500 Subject: [PATCH 01/13] Cleanup code for getting IDEs associated with SimChannels --- sbncode/CAFMaker/CAFMaker_module.cc | 2 +- sbncode/CAFMaker/FillTrue.cxx | 26 +++++++++++++------------- sbncode/CAFMaker/FillTrue.h | 13 ++++++++++--- 3 files changed, 24 insertions(+), 17 deletions(-) diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index 2728b395f..55d18bacc 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -1374,7 +1374,7 @@ void CAFMaker::produce(art::Event& evt) noexcept { } // Prep truth-to-reco-matching info - std::map>> id_to_ide_map; + std::map> id_to_ide_map; std::map>> id_to_truehit_map; std::map id_to_hit_energy_map; diff --git a/sbncode/CAFMaker/FillTrue.cxx b/sbncode/CAFMaker/FillTrue.cxx index 491c813c0..1cfb99dee 100644 --- a/sbncode/CAFMaker/FillTrue.cxx +++ b/sbncode/CAFMaker/FillTrue.cxx @@ -136,7 +136,7 @@ namespace caf { }//FillTrackTruth // Assumes truth matching and calo-points are filled - void FillTrackCaloTruth(const std::map>> &id_to_ide_map, + void FillTrackCaloTruth(const std::map> &id_to_ide_map, const std::vector &mc_particles, const geo::GeometryCore& geometry, const geo::WireReadoutGeom& wireReadout, @@ -163,10 +163,10 @@ namespace caf { // Load the hits // match on the channel, which is unique - const std::vector> &match_ides = id_to_ide_map.at(srtrack.truth.p.G4ID); + const std::vector &match_ides = id_to_ide_map.at(srtrack.truth.p.G4ID); std::map> chan_2_ides; - for (auto const &ide_pair: match_ides) { - chan_2_ides[wireReadout.PlaneWireToChannel(ide_pair.first)].push_back(ide_pair.second); + for (auto const &ide_p: match_ides) { + chan_2_ides[wireReadout.PlaneWireToChannel(ide_p.wire)].push_back(ide_p.ide); } // pre-compute partial ranges @@ -634,15 +634,15 @@ namespace caf { void FillTrueG4Particle(const simb::MCParticle &particle, const std::vector &active_volumes, const std::vector> &tpc_volumes, - const std::map>> &id_to_ide_map, + const std::map> &id_to_ide_map, const std::map>> &id_to_truehit_map, const cheat::BackTrackerService &backtracker, const cheat::ParticleInventoryService &inventory_service, const std::vector> &neutrinos, caf::SRTrueParticle &srparticle) { - std::vector> empty; - const std::vector> &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty; + std::vector empty; + const std::vector &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty; std::vector> emptyHits; const std::vector> &particle_hits = id_to_truehit_map.count(particle.TrackId()) ? id_to_truehit_map.at(particle.TrackId()) : emptyHits; @@ -662,9 +662,9 @@ namespace caf { } } - for (auto const &ide_pair: particle_ides) { - const geo::WireID &w = ide_pair.first; - const sim::IDE *ide = ide_pair.second; + for (auto const &ide_p: particle_ides) { + const geo::WireID &w = ide_p.wire; + const sim::IDE *ide = ide_p.ide; if(w.Plane >= 0 && w.Plane < 3 && w.Cryostat < 2){ srparticle.plane[w.Cryostat][w.Plane].visE += ide->energy / 1000. /* MeV -> GeV*/; @@ -882,8 +882,8 @@ namespace caf { return ret; } - std::map>> PrepSimChannels(const std::vector> &simchannels, const geo::WireReadoutGeom &wireReadout) { - std::map>> ret; + std::map> PrepSimChannels(const std::vector> &simchannels, const geo::WireReadoutGeom &wireReadout) { + std::map> ret; for (const art::Ptr sc : simchannels) { // Lookup the wire of this channel @@ -895,7 +895,7 @@ namespace caf { for (const auto &item : sc->TDCIDEMap()) { for (const sim::IDE &ide: item.second) { // indexing initializes empty vector - ret[abs(ide.trackID)].push_back({thisWire, &ide}); + ret[abs(ide.trackID)].push_back({thisWire, item.first, &ide}); } } } diff --git a/sbncode/CAFMaker/FillTrue.h b/sbncode/CAFMaker/FillTrue.h index 419ed6ee9..da05edebe 100644 --- a/sbncode/CAFMaker/FillTrue.h +++ b/sbncode/CAFMaker/FillTrue.h @@ -40,6 +40,13 @@ namespace caf float totE; }; + struct ParticleIDE { + geo::WireID wire; + unsigned short tick; + const sim::IDE *ide; + }; + + // Helpers caf::Wall_t GetWallCross( const geo::BoxBoundedGeo &volume, const TVector3 p0, @@ -73,7 +80,7 @@ namespace caf void FillTrueG4Particle(const simb::MCParticle &particle, const std::vector &active_volumes, const std::vector> &tpc_volumes, - const std::map>> &id_to_ide_map, + const std::map> &id_to_ide_map, const std::map>> &id_to_truehit_map, const cheat::BackTrackerService &backtracker, const cheat::ParticleInventoryService &inventory_service, @@ -103,7 +110,7 @@ namespace caf caf::SRTrack& srtrack, bool allowEmpty = false); - void FillTrackCaloTruth(const std::map>> &id_to_ide_map, + void FillTrackCaloTruth(const std::map> &id_to_ide_map, const std::vector &mc_particles, const geo::GeometryCore & geometry, const geo::WireReadoutGeom& wireReadout, @@ -133,7 +140,7 @@ namespace caf CLHEP::HepRandomEngine &rand, std::vector &srfakereco); - std::map>> PrepSimChannels(const std::vector> &simchannels, const geo::WireReadoutGeom &wireReadout); + std::map> PrepSimChannels(const std::vector> &simchannels, const geo::WireReadoutGeom &wireReadout); std::map>> PrepTrueHits(const std::vector> &allHits, const detinfo::DetectorClocksData &clockData, const cheat::BackTrackerService &backtracker); std::map SetupIDHitEnergyMap(const std::vector> &allHits, const detinfo::DetectorClocksData &clockData, From 90cdf6511b78b801e2952613106b0c5c88d8a8b2 Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Wed, 6 Aug 2025 19:50:20 -0500 Subject: [PATCH 02/13] Draft of code to get costh_drift (cos_phi, e.g.) and the efield into the CAF SRCalo object per-hit. --- sbncode/CAFMaker/CAFMaker_module.cc | 7 +++-- sbncode/CAFMaker/FillReco.cxx | 44 ++++++++++++++++++++++++++++- sbncode/CAFMaker/FillReco.h | 5 ++++ 3 files changed, 52 insertions(+), 4 deletions(-) diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index 55d18bacc..72e2e3405 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -92,6 +92,7 @@ #include "lardataobj/RecoBase/PFParticle.h" #include "lardataobj/RecoBase/Slice.h" #include "lardataobj/RecoBase/Track.h" +#include "lardataobj/RecoBase/TrackHitMeta.h" #include "lardataobj/RecoBase/Vertex.h" #include "lardataobj/RecoBase/Shower.h" #include "lardataobj/RecoBase/MCSFitResult.h" @@ -1906,8 +1907,8 @@ void CAFMaker::produce(art::Event& evt) noexcept { FindManyPStrict(fmPFPart, evt, fParams.PFParticleLabel() + slice_tag_suff); - art::FindManyP fmTrackHit = - FindManyPStrict(slcTracks, evt, + art::FindManyP fmTrackHit = + FindManyPDStrict(slcTracks, evt, fParams.RecoTrackLabel() + slice_tag_suff); art::FindManyP fmShowerHit = @@ -2152,7 +2153,7 @@ void CAFMaker::produce(art::Event& evt) noexcept { FillTrackDazzle(fmTrackDazzle.at(iPart).front(), trk); } if (fmCalo.isValid()) { - FillTrackCalo(fmCalo.at(iPart), fmTrackHit.at(iPart), + FillTrackCalo(fmCalo.at(iPart), *thisTrack[0], fmTrackHit.at(iPart), fmTrackHit.data(iPart), (fParams.FillHitsNeutrinoSlices() && NeutrinoSlice) || fParams.FillHitsAllSlices(), fParams.TrackHitFillRRStartCut(), fParams.TrackHitFillRREndCut(), dprop, trk); diff --git a/sbncode/CAFMaker/FillReco.cxx b/sbncode/CAFMaker/FillReco.cxx index 59d8e69f3..2f5f4b9ae 100644 --- a/sbncode/CAFMaker/FillReco.cxx +++ b/sbncode/CAFMaker/FillReco.cxx @@ -10,6 +10,10 @@ #include "FillReco.h" #include "RecoUtils/RecoUtils.h" +#include "larevt/SpaceCharge/SpaceCharge.h" +#include "larevt/SpaceChargeServices/SpaceChargeService.h" +#include "larcore/CoreUtils/ServiceUtil.h" + namespace caf { @@ -720,8 +724,28 @@ namespace caf } } + // Helper function: get the e field + double GetEfield(const detinfo::DetectorPropertiesData& dprop, const geo::Point_t loc) { + auto const* sce = lar::providerFrom(); + + double EField = dprop.Efield(); + if (sce->EnableSimEfieldSCE()) { + // Gets relative E field Distortions + geo::Vector_t EFieldOffsets = sce->GetEfieldOffsets(loc); + // Add 1 in X direction as this is the direction of the drift field + EFieldOffsets = EFieldOffsets + geo::Vector_t{1, 0, 0}; + // Convert to Absolute E Field from relative + EFieldOffsets = EField * EFieldOffsets; + // We only care about the magnitude for recombination + EField = EFieldOffsets.r(); + } + return EField; + } + void FillTrackPlaneCalo(const anab::Calorimetry &calo, + const recob::Track& track, const std::vector> &hits, + const std::vector& thms, bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend, const detinfo::DetectorPropertiesData &dprop, caf::SRTrackCalo &srcalo) { @@ -770,6 +794,22 @@ namespace caf p.mult = h->Multiplicity(); p.start = h->StartTick(); p.end = h->EndTick(); + + + // Get the trajectory point index from this hit. Again -- this is too hard. + // + // Use this to get the (SCE corrected) efield and the angle to the drift direction + unsigned traj_point_index = thms.at(h.key())->Index(); + unsigned int int_max_as_unsigned_int{std::numeric_limits::max()}; + if (traj_point_index != int_max_as_unsigned_int && // invalid + track.HasValidPoint(traj_point_index)) { + float costh_drift = track.DirectionAtPoint(traj_point_index).X(); + // p.costh_drift = costh_drift; + float efield = GetEfield(dprop, track.LocationAtPoint(traj_point_index)); + // p.efield = efield; + (void) costh_drift; + (void) efield; + } } } @@ -823,7 +863,9 @@ namespace caf } void FillTrackCalo(const std::vector> &calos, + const recob::Track& track, const std::vector> &hits, + const std::vector& thms, bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend, const detinfo::DetectorPropertiesData &dprop, caf::SRTrack& srtrack, @@ -838,7 +880,7 @@ namespace caf if (calo.PlaneID()) { unsigned plane_id = calo.PlaneID().Plane; assert(plane_id < 3); - FillTrackPlaneCalo(calo, hits, fill_calo_points, fillhit_rrstart, fillhit_rrend, dprop, srtrack.calo[plane_id]); + FillTrackPlaneCalo(calo, track, hits, thms, fill_calo_points, fillhit_rrstart, fillhit_rrend, dprop, srtrack.calo[plane_id]); } } diff --git a/sbncode/CAFMaker/FillReco.h b/sbncode/CAFMaker/FillReco.h index c6df1ab39..bd332e109 100644 --- a/sbncode/CAFMaker/FillReco.h +++ b/sbncode/CAFMaker/FillReco.h @@ -14,6 +14,7 @@ #include "lardataobj/RecoBase/Shower.h" #include "lardataobj/RecoBase/Slice.h" #include "lardataobj/RecoBase/Track.h" +#include "lardataobj/RecoBase/TrackHitMeta.h" #include "lardataobj/RecoBase/Vertex.h" #include "lardataobj/RecoBase/Hit.h" #include "lardataobj/RecoBase/SpacePoint.h" @@ -168,7 +169,9 @@ namespace caf bool allowEmpty = false); void FillTrackPlaneCalo(const anab::Calorimetry &calo, + const recob::Track& track, const std::vector> &hits, + const std::vector& thms, bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend, const detinfo::DetectorPropertiesData &dprop, caf::SRTrackCalo &srcalo); @@ -186,7 +189,9 @@ namespace caf bool allowEmpty = false); void FillTrackCalo(const std::vector> &calos, + const recob::Track& track, const std::vector> &hits, + const std::vector& thms, bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend, const detinfo::DetectorPropertiesData &dprop, caf::SRTrack& srtrack, From 98b624d052816da45dd075b1bca863e91671fd6a Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Thu, 7 Aug 2025 12:10:46 -0500 Subject: [PATCH 03/13] Incorporate updates to CAF calo-truth matching into calibration ntuple'r. Use actualy true deposition time to define true calo point time. --- sbncode/Calibration/TrackCaloSkimmer.h | 3 +- .../Calibration/TrackCaloSkimmer_module.cc | 33 ++++++++++--------- 2 files changed, 19 insertions(+), 17 deletions(-) diff --git a/sbncode/Calibration/TrackCaloSkimmer.h b/sbncode/Calibration/TrackCaloSkimmer.h index c35f9ecd5..5efd61e64 100644 --- a/sbncode/Calibration/TrackCaloSkimmer.h +++ b/sbncode/Calibration/TrackCaloSkimmer.h @@ -64,6 +64,7 @@ #include "sbnobj/Common/CRT/CRTHitT0TaggingInfo.cc" #include "sbnobj/Common/CRT/CRTHitT0TaggingTruthInfo.cc" +#include "sbncode/CAFMaker/FillTrue.h" #include "ITCSSelectionTool.h" namespace sbn { @@ -161,7 +162,7 @@ class sbn::TrackCaloSkimmer : public art::EDAnalyzer { const std::vector> &mcparticles, const std::vector &active_volumes, const std::vector> &tpc_volumes, - const std::map>> id_to_ide_map, + const std::map> id_to_ide_map, const std::map>> id_to_truehit_map, const detinfo::DetectorPropertiesData &dprop, const geo::GeometryCore *geo, diff --git a/sbncode/Calibration/TrackCaloSkimmer_module.cc b/sbncode/Calibration/TrackCaloSkimmer_module.cc index 9dc0be1ac..7e9e57740 100644 --- a/sbncode/Calibration/TrackCaloSkimmer_module.cc +++ b/sbncode/Calibration/TrackCaloSkimmer_module.cc @@ -289,7 +289,7 @@ void sbn::TrackCaloSkimmer::analyze(art::Event const& e) // Prep truth-to-reco-matching info // // Use helper functions from CAFMaker/FillTrue - std::map>> id_to_ide_map; + std::map> id_to_ide_map; std::map>> id_to_truehit_map; const cheat::BackTrackerService *bt = NULL; @@ -545,14 +545,14 @@ geo::Point_t WireToTrajectoryPosition(const geo::Point_t &loc, const geo::TPCID sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, const std::vector &active_volumes, const std::vector> &tpc_volumes, - const std::map>> &id_to_ide_map, + const std::map> &id_to_ide_map, const std::map>> &id_to_truehit_map, const detinfo::DetectorPropertiesData &dprop, const geo::GeometryCore *geo, const geo::WireReadoutGeom *wireReadout) { - std::vector> empty; - const std::vector> &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty; + std::vector empty; + const std::vector &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty; std::vector> emptyHits; const std::vector> &particle_hits = id_to_truehit_map.count(particle.TrackId()) ? id_to_truehit_map.at(particle.TrackId()) : emptyHits; @@ -569,9 +569,9 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, trueparticle.plane0nhit = 0; trueparticle.plane1nhit = 0; trueparticle.plane2nhit = 0; - for (auto const &ide_pair: particle_ides) { - const geo::WireID &w = ide_pair.first; - const sim::IDE *ide = ide_pair.second; + for (auto const &ide_p: particle_ides) { + const geo::WireID &w = ide_p.wire; + const sim::IDE *ide = ide_p.ide; if (w.Plane == 0) { trueparticle.plane0VisE += ide->energy / 1000. /* MeV -> GeV*/; @@ -720,10 +720,11 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, // Organize deposition info into per-wire true "Hits" -- key is the Channel Number std::map truehits; - for (auto const &ide_pair: particle_ides) { - const geo::WireID &w = ide_pair.first; + for (auto const &ide_part: particle_ides) { + const geo::WireID &w = ide_part.wire; unsigned c = wireReadout->PlaneWireToChannel(w); - const sim::IDE *ide = ide_pair.second; + const sim::IDE *ide = ide_part.ide; + unsigned short tick = ide_part.tick; // Set stuff truehits[c].cryo = w.Cryostat; @@ -739,6 +740,8 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, truehits[c].p.y = (truehits[c].p.y*old_elec + ide->y*ide->numElectrons) / new_elec; truehits[c].p.z = (truehits[c].p.z*old_elec + ide->z*ide->numElectrons) / new_elec; + truehits[c].time = (truehits[c].time*old_elec + tick*ide->numElectrons) / new_elec; + // Also get the position with space charge un-done geo::Point_t ide_p(ide->x, ide->y, ide->z); geo::Point_t ide_p_scecorr = WireToTrajectoryPosition(ide_p, w); @@ -754,10 +757,10 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, } // Compute widths - for (auto const &ide_pair: particle_ides) { - const geo::WireID &w = ide_pair.first; + for (auto const &ide_part: particle_ides) { + const geo::WireID &w = ide_part.wire; unsigned c = wireReadout->PlaneWireToChannel(w); - const sim::IDE *ide = ide_pair.second; + const sim::IDE *ide = ide_part.ide; geo::Point_t ide_p(ide->x, ide->y, ide->z); geo::Point_t ide_p_scecorr = WireToTrajectoryPosition(ide_p, w); @@ -782,8 +785,6 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, // Compute the time of each hit for (sbn::TrueHit &h: truehits_v) { - h.time = dprop.ConvertXToTicks(h.p.x, h.plane, h.tpc, h.cryo); - double xdrift = abs(h.p.x - wireReadout->Plane(geo::PlaneID(h.cryo, h.tpc, 0)).GetCenter().X()); h.tdrift = xdrift / dprop.DriftVelocity(); } @@ -1010,7 +1011,7 @@ void sbn::TrackCaloSkimmer::FillTrackTruth(const detinfo::DetectorClocksData &cl const std::vector> &mcparticles, const std::vector &active_volumes, const std::vector> &tpc_volumes, - const std::map>> id_to_ide_map, + const std::map> id_to_ide_map, const std::map>> id_to_truehit_map, const detinfo::DetectorPropertiesData &dprop, const geo::GeometryCore *geo, From aa5484cc87f6c97876c5e976a0422e9d006fc816 Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Tue, 26 Aug 2025 16:45:43 -0500 Subject: [PATCH 04/13] Lookup TrackHitMeta object correctly. --- sbncode/CAFMaker/FillReco.cxx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/sbncode/CAFMaker/FillReco.cxx b/sbncode/CAFMaker/FillReco.cxx index 2f5f4b9ae..233f2df28 100644 --- a/sbncode/CAFMaker/FillReco.cxx +++ b/sbncode/CAFMaker/FillReco.cxx @@ -782,7 +782,8 @@ namespace caf // lookup the wire -- the Calorimery object makes this // __way__ harder than it should be - for (const art::Ptr &h: hits) { + for (unsigned i_hit = 0; i_hit < hits.size(); i_hit++) { + const art::Ptr &h = hits[i_hit]; if (h.key() == tps[i]) { p.wire = h->WireID().Wire; p.tpc = h->WireID().TPC; @@ -799,7 +800,7 @@ namespace caf // Get the trajectory point index from this hit. Again -- this is too hard. // // Use this to get the (SCE corrected) efield and the angle to the drift direction - unsigned traj_point_index = thms.at(h.key())->Index(); + unsigned traj_point_index = thms.at(i_hit)->Index(); unsigned int int_max_as_unsigned_int{std::numeric_limits::max()}; if (traj_point_index != int_max_as_unsigned_int && // invalid track.HasValidPoint(traj_point_index)) { From 455cb8da92b9913b9484932a4e0a1d515e7a5c7d Mon Sep 17 00:00:00 2001 From: sungbinoh Date: Wed, 27 Aug 2025 10:12:50 -0500 Subject: [PATCH 05/13] adding phi --- sbncode/CAFMaker/FillReco.cxx | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/sbncode/CAFMaker/FillReco.cxx b/sbncode/CAFMaker/FillReco.cxx index 96873e6ec..c91930c69 100644 --- a/sbncode/CAFMaker/FillReco.cxx +++ b/sbncode/CAFMaker/FillReco.cxx @@ -879,11 +879,10 @@ namespace caf if (traj_point_index != int_max_as_unsigned_int && // invalid track.HasValidPoint(traj_point_index)) { float costh_drift = track.DirectionAtPoint(traj_point_index).X(); - // p.costh_drift = costh_drift; + float phi = acos(abs(costh_drift)) * 180. / M_PI; float efield = GetEfield(dprop, track.LocationAtPoint(traj_point_index)); - // p.efield = efield; - (void) costh_drift; - (void) efield; + p.efield = efield; + p.phi = phi; } } } From e70f014aa2f887010be696f74c338e924c834fe2 Mon Sep 17 00:00:00 2001 From: sungbinoh Date: Fri, 29 Aug 2025 12:05:22 -0500 Subject: [PATCH 06/13] following Gianluca's comments part1 --- sbncode/CAFMaker/FillReco.cxx | 16 +++++++--------- sbncode/CAFMaker/FillTrue.h | 7 ++++--- 2 files changed, 11 insertions(+), 12 deletions(-) diff --git a/sbncode/CAFMaker/FillReco.cxx b/sbncode/CAFMaker/FillReco.cxx index c91930c69..78ffae458 100644 --- a/sbncode/CAFMaker/FillReco.cxx +++ b/sbncode/CAFMaker/FillReco.cxx @@ -804,12 +804,12 @@ namespace caf double EField = dprop.Efield(); if (sce->EnableSimEfieldSCE()) { - // Gets relative E field Distortions + // Gets fractional E field distortions w.r.t. the nominal field on x-axis geo::Vector_t EFieldOffsets = sce->GetEfieldOffsets(loc); - // Add 1 in X direction as this is the direction of the drift field + // Add 1 in X direction as this is the direction of the drift field, not caring if it is +x or -x direction, since we only want |E| EFieldOffsets = EFieldOffsets + geo::Vector_t{1, 0, 0}; // Convert to Absolute E Field from relative - EFieldOffsets = EField * EFieldOffsets; + EFieldOffsets *= EField; // We only care about the magnitude for recombination EField = EFieldOffsets.r(); } @@ -875,14 +875,12 @@ namespace caf // // Use this to get the (SCE corrected) efield and the angle to the drift direction unsigned traj_point_index = thms.at(i_hit)->Index(); - unsigned int int_max_as_unsigned_int{std::numeric_limits::max()}; - if (traj_point_index != int_max_as_unsigned_int && // invalid - track.HasValidPoint(traj_point_index)) { + if (track.HasValidPoint(traj_point_index)) { float costh_drift = track.DirectionAtPoint(traj_point_index).X(); - float phi = acos(abs(costh_drift)) * 180. / M_PI; + float phi = acos(abs(costh_drift)) * 180. / M_PI; float efield = GetEfield(dprop, track.LocationAtPoint(traj_point_index)); - p.efield = efield; - p.phi = phi; + p.efield = efield; + p.phi = phi; } } } diff --git a/sbncode/CAFMaker/FillTrue.h b/sbncode/CAFMaker/FillTrue.h index da05edebe..d5d3e3d58 100644 --- a/sbncode/CAFMaker/FillTrue.h +++ b/sbncode/CAFMaker/FillTrue.h @@ -40,10 +40,11 @@ namespace caf float totE; }; + /// Ioniosation charge (`sim::IDE`) associated to its TPC readout information. struct ParticleIDE { - geo::WireID wire; - unsigned short tick; - const sim::IDE *ide; + geo::WireID wire; ///< Wire on a given plane closest to the drift path of the charge. + unsigned short tick = 0; ///< Time tick at which the charge passes closest to the wire. + const sim::IDE *ide = nullptr; ///< Deposited charge information. }; From 2feb51f7a0ea7b64e16064428a49bb788b04427b Mon Sep 17 00:00:00 2001 From: sungbinoh Date: Fri, 29 Aug 2025 12:10:15 -0500 Subject: [PATCH 07/13] following Gianluca's suggestion part 2 --- sbncode/CAFMaker/FillReco.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sbncode/CAFMaker/FillReco.cxx b/sbncode/CAFMaker/FillReco.cxx index 78ffae458..92691d596 100644 --- a/sbncode/CAFMaker/FillReco.cxx +++ b/sbncode/CAFMaker/FillReco.cxx @@ -807,7 +807,7 @@ namespace caf // Gets fractional E field distortions w.r.t. the nominal field on x-axis geo::Vector_t EFieldOffsets = sce->GetEfieldOffsets(loc); // Add 1 in X direction as this is the direction of the drift field, not caring if it is +x or -x direction, since we only want |E| - EFieldOffsets = EFieldOffsets + geo::Vector_t{1, 0, 0}; + EFieldOffsets += geo::Vector_t{1, 0, 0}; // Convert to Absolute E Field from relative EFieldOffsets *= EField; // We only care about the magnitude for recombination From c6a3c290c859b97bbdb052507b2894e1fd793381 Mon Sep 17 00:00:00 2001 From: sungbinoh Date: Fri, 29 Aug 2025 12:36:51 -0500 Subject: [PATCH 08/13] ParticleIDE -> ReadoutIDE --- sbncode/CAFMaker/FillTrue.cxx | 18 +++++++++--------- sbncode/CAFMaker/FillTrue.h | 8 ++++---- sbncode/Calibration/TrackCaloSkimmer.h | 2 +- sbncode/Calibration/TrackCaloSkimmer_module.cc | 12 ++++++------ 4 files changed, 20 insertions(+), 20 deletions(-) diff --git a/sbncode/CAFMaker/FillTrue.cxx b/sbncode/CAFMaker/FillTrue.cxx index 1cfb99dee..7ef546af7 100644 --- a/sbncode/CAFMaker/FillTrue.cxx +++ b/sbncode/CAFMaker/FillTrue.cxx @@ -136,7 +136,7 @@ namespace caf { }//FillTrackTruth // Assumes truth matching and calo-points are filled - void FillTrackCaloTruth(const std::map> &id_to_ide_map, + void FillTrackCaloTruth(const std::map> &id_to_ide_map, const std::vector &mc_particles, const geo::GeometryCore& geometry, const geo::WireReadoutGeom& wireReadout, @@ -163,7 +163,7 @@ namespace caf { // Load the hits // match on the channel, which is unique - const std::vector &match_ides = id_to_ide_map.at(srtrack.truth.p.G4ID); + const std::vector &match_ides = id_to_ide_map.at(srtrack.truth.p.G4ID); std::map> chan_2_ides; for (auto const &ide_p: match_ides) { chan_2_ides[wireReadout.PlaneWireToChannel(ide_p.wire)].push_back(ide_p.ide); @@ -634,15 +634,15 @@ namespace caf { void FillTrueG4Particle(const simb::MCParticle &particle, const std::vector &active_volumes, const std::vector> &tpc_volumes, - const std::map> &id_to_ide_map, + const std::map> &id_to_ide_map, const std::map>> &id_to_truehit_map, const cheat::BackTrackerService &backtracker, const cheat::ParticleInventoryService &inventory_service, const std::vector> &neutrinos, caf::SRTrueParticle &srparticle) { - std::vector empty; - const std::vector &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty; + std::vector empty; + const std::vector &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty; std::vector> emptyHits; const std::vector> &particle_hits = id_to_truehit_map.count(particle.TrackId()) ? id_to_truehit_map.at(particle.TrackId()) : emptyHits; @@ -882,8 +882,8 @@ namespace caf { return ret; } - std::map> PrepSimChannels(const std::vector> &simchannels, const geo::WireReadoutGeom &wireReadout) { - std::map> ret; + std::map> PrepSimChannels(const std::vector> &simchannels, const geo::WireReadoutGeom &wireReadout) { + std::map> ret; for (const art::Ptr sc : simchannels) { // Lookup the wire of this channel @@ -1487,7 +1487,7 @@ caf::SRTrackTruth MatchTrack2Truth(const detinfo::DetectorClocksData &clockData, art::ServiceHandle bt_serv; // this id is the same as the mcparticle ID as long as we got it from geant4 - std::vector> matches = CAFRecoUtils::AllTrueParticleIDEnergyMatches(clockData, hits, true); + std::vector> matches = CAFRecoUtils::AllTrueReadoutIDEnergyMatches(clockData, hits, true); float total_energy = CAFRecoUtils::TotalHitEnergy(clockData, hits); std::map track_hits_map = caf::SetupIDHitEnergyMap(hits, clockData, *bt_serv.get()); @@ -1586,7 +1586,7 @@ caf::SRTruthMatch MatchSlice2Truth(const std::vector> &hits ret.index = -1; return ret; } - std::vector> matches = CAFRecoUtils::AllTrueParticleIDEnergyMatches(clockData, hits, true); + std::vector> matches = CAFRecoUtils::AllTrueReadoutIDEnergyMatches(clockData, hits, true); std::vector matching_energy(truths.size(), 0.); for (auto const &pair: matches) { art::Ptr truth; diff --git a/sbncode/CAFMaker/FillTrue.h b/sbncode/CAFMaker/FillTrue.h index d5d3e3d58..be5f48e69 100644 --- a/sbncode/CAFMaker/FillTrue.h +++ b/sbncode/CAFMaker/FillTrue.h @@ -41,7 +41,7 @@ namespace caf }; /// Ioniosation charge (`sim::IDE`) associated to its TPC readout information. - struct ParticleIDE { + struct ReadoutIDE { geo::WireID wire; ///< Wire on a given plane closest to the drift path of the charge. unsigned short tick = 0; ///< Time tick at which the charge passes closest to the wire. const sim::IDE *ide = nullptr; ///< Deposited charge information. @@ -81,7 +81,7 @@ namespace caf void FillTrueG4Particle(const simb::MCParticle &particle, const std::vector &active_volumes, const std::vector> &tpc_volumes, - const std::map> &id_to_ide_map, + const std::map> &id_to_ide_map, const std::map>> &id_to_truehit_map, const cheat::BackTrackerService &backtracker, const cheat::ParticleInventoryService &inventory_service, @@ -111,7 +111,7 @@ namespace caf caf::SRTrack& srtrack, bool allowEmpty = false); - void FillTrackCaloTruth(const std::map> &id_to_ide_map, + void FillTrackCaloTruth(const std::map> &id_to_ide_map, const std::vector &mc_particles, const geo::GeometryCore & geometry, const geo::WireReadoutGeom& wireReadout, @@ -141,7 +141,7 @@ namespace caf CLHEP::HepRandomEngine &rand, std::vector &srfakereco); - std::map> PrepSimChannels(const std::vector> &simchannels, const geo::WireReadoutGeom &wireReadout); + std::map> PrepSimChannels(const std::vector> &simchannels, const geo::WireReadoutGeom &wireReadout); std::map>> PrepTrueHits(const std::vector> &allHits, const detinfo::DetectorClocksData &clockData, const cheat::BackTrackerService &backtracker); std::map SetupIDHitEnergyMap(const std::vector> &allHits, const detinfo::DetectorClocksData &clockData, diff --git a/sbncode/Calibration/TrackCaloSkimmer.h b/sbncode/Calibration/TrackCaloSkimmer.h index c3941620d..e6e8f3aa3 100644 --- a/sbncode/Calibration/TrackCaloSkimmer.h +++ b/sbncode/Calibration/TrackCaloSkimmer.h @@ -162,7 +162,7 @@ class sbn::TrackCaloSkimmer : public art::EDAnalyzer { const std::vector> &mcparticles, const std::vector &active_volumes, const std::vector> &tpc_volumes, - const std::map> id_to_ide_map, + const std::map> id_to_ide_map, const std::map>> id_to_truehit_map, const detinfo::DetectorPropertiesData &dprop, const geo::GeometryCore *geo, diff --git a/sbncode/Calibration/TrackCaloSkimmer_module.cc b/sbncode/Calibration/TrackCaloSkimmer_module.cc index 711eff1ca..71b8bc046 100644 --- a/sbncode/Calibration/TrackCaloSkimmer_module.cc +++ b/sbncode/Calibration/TrackCaloSkimmer_module.cc @@ -289,7 +289,7 @@ void sbn::TrackCaloSkimmer::analyze(art::Event const& e) // Prep truth-to-reco-matching info // // Use helper functions from CAFMaker/FillTrue - std::map> id_to_ide_map; + std::map> id_to_ide_map; std::map>> id_to_truehit_map; const cheat::BackTrackerService *bt = NULL; @@ -545,14 +545,14 @@ geo::Point_t WireToTrajectoryPosition(const geo::Point_t &loc, const geo::TPCID sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, const std::vector &active_volumes, const std::vector> &tpc_volumes, - const std::map> &id_to_ide_map, + const std::map> &id_to_ide_map, const std::map>> &id_to_truehit_map, const detinfo::DetectorPropertiesData &dprop, const geo::GeometryCore *geo, const geo::WireReadoutGeom *wireReadout) { - std::vector empty; - const std::vector &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty; + std::vector empty; + const std::vector &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty; std::vector> emptyHits; const std::vector> &particle_hits = id_to_truehit_map.count(particle.TrackId()) ? id_to_truehit_map.at(particle.TrackId()) : emptyHits; @@ -1011,14 +1011,14 @@ void sbn::TrackCaloSkimmer::FillTrackTruth(const detinfo::DetectorClocksData &cl const std::vector> &mcparticles, const std::vector &active_volumes, const std::vector> &tpc_volumes, - const std::map> id_to_ide_map, + const std::map> id_to_ide_map, const std::map>> id_to_truehit_map, const detinfo::DetectorPropertiesData &dprop, const geo::GeometryCore *geo, const geo::WireReadoutGeom *wireReadout) { // Lookup the true-particle match -- use utils in CAF - std::vector> matches = CAFRecoUtils::AllTrueParticleIDEnergyMatches(clock_data, trkHits, true); + std::vector> matches = CAFRecoUtils::AllTrueReadoutIDEnergyMatches(clock_data, trkHits, true); float total_energy = CAFRecoUtils::TotalHitEnergy(clock_data, trkHits); fTrack->truth.depE = total_energy / 1000. /* MeV -> GeV */; From 035418fc4921f6f4567263f39ebcc1563145ba99 Mon Sep 17 00:00:00 2001 From: sungbinoh Date: Fri, 29 Aug 2025 14:01:51 -0500 Subject: [PATCH 09/13] minor bug fixes --- sbncode/CAFMaker/CAFMaker_module.cc | 2 +- sbncode/CAFMaker/FillTrue.cxx | 4 ++-- sbncode/Calibration/TrackCaloSkimmer_module.cc | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index eb506f346..740617d9c 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -1387,7 +1387,7 @@ void CAFMaker::produce(art::Event& evt) noexcept { } // Prep truth-to-reco-matching info - std::map> id_to_ide_map; + std::map> id_to_ide_map; std::map>> id_to_truehit_map; std::map id_to_hit_energy_map; diff --git a/sbncode/CAFMaker/FillTrue.cxx b/sbncode/CAFMaker/FillTrue.cxx index 7ef546af7..f7d374ef7 100644 --- a/sbncode/CAFMaker/FillTrue.cxx +++ b/sbncode/CAFMaker/FillTrue.cxx @@ -1487,7 +1487,7 @@ caf::SRTrackTruth MatchTrack2Truth(const detinfo::DetectorClocksData &clockData, art::ServiceHandle bt_serv; // this id is the same as the mcparticle ID as long as we got it from geant4 - std::vector> matches = CAFRecoUtils::AllTrueReadoutIDEnergyMatches(clockData, hits, true); + std::vector> matches = CAFRecoUtils::AllTrueParticleIDEnergyMatches(clockData, hits, true); float total_energy = CAFRecoUtils::TotalHitEnergy(clockData, hits); std::map track_hits_map = caf::SetupIDHitEnergyMap(hits, clockData, *bt_serv.get()); @@ -1586,7 +1586,7 @@ caf::SRTruthMatch MatchSlice2Truth(const std::vector> &hits ret.index = -1; return ret; } - std::vector> matches = CAFRecoUtils::AllTrueReadoutIDEnergyMatches(clockData, hits, true); + std::vector> matches = CAFRecoUtils::AllTrueParticleIDEnergyMatches(clockData, hits, true); std::vector matching_energy(truths.size(), 0.); for (auto const &pair: matches) { art::Ptr truth; diff --git a/sbncode/Calibration/TrackCaloSkimmer_module.cc b/sbncode/Calibration/TrackCaloSkimmer_module.cc index 71b8bc046..43f018ead 100644 --- a/sbncode/Calibration/TrackCaloSkimmer_module.cc +++ b/sbncode/Calibration/TrackCaloSkimmer_module.cc @@ -1018,7 +1018,7 @@ void sbn::TrackCaloSkimmer::FillTrackTruth(const detinfo::DetectorClocksData &cl const geo::WireReadoutGeom *wireReadout) { // Lookup the true-particle match -- use utils in CAF - std::vector> matches = CAFRecoUtils::AllTrueReadoutIDEnergyMatches(clock_data, trkHits, true); + std::vector> matches = CAFRecoUtils::AllTrueParticleIDEnergyMatches(clock_data, trkHits, true); float total_energy = CAFRecoUtils::TotalHitEnergy(clock_data, trkHits); fTrack->truth.depE = total_energy / 1000. /* MeV -> GeV */; From a6b090acd5175b5729cb64732564bd39a1326eb1 Mon Sep 17 00:00:00 2001 From: sungbinoh Date: Tue, 2 Sep 2025 16:42:19 -0500 Subject: [PATCH 10/13] moving ReadoutIDE to TrackCaloSkimmer.h and in sbn namespace. refer to this from CAFMaker's side --- sbncode/CAFMaker/CAFMaker_module.cc | 2 +- sbncode/CAFMaker/CMakeLists.txt | 1 + sbncode/CAFMaker/FillTrue.cxx | 14 +++++++------- sbncode/CAFMaker/FillTrue.h | 16 +++++----------- sbncode/Calibration/TrackCaloSkimmer.h | 10 +++++++--- sbncode/Calibration/TrackCaloSkimmer_module.cc | 10 +++++----- 6 files changed, 26 insertions(+), 27 deletions(-) diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index 740617d9c..0c11332ab 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -1387,7 +1387,7 @@ void CAFMaker::produce(art::Event& evt) noexcept { } // Prep truth-to-reco-matching info - std::map> id_to_ide_map; + std::map> id_to_ide_map; std::map>> id_to_truehit_map; std::map id_to_hit_energy_map; diff --git a/sbncode/CAFMaker/CMakeLists.txt b/sbncode/CAFMaker/CMakeLists.txt index 76c3903d4..4ef86971c 100644 --- a/sbncode/CAFMaker/CMakeLists.txt +++ b/sbncode/CAFMaker/CMakeLists.txt @@ -29,6 +29,7 @@ art_make_library( LIBRARY_NAME sbncode_CAFMaker larcorealg::Geometry larcore::Geometry_Geometry_service larcorealg::GeoAlgo + lardata::DetectorInfoServices_DetectorPropertiesServiceStandard_service larsim::MCCheater_BackTrackerService_service nusimdata::SimulationBase larsim::MCCheater_ParticleInventoryService_service diff --git a/sbncode/CAFMaker/FillTrue.cxx b/sbncode/CAFMaker/FillTrue.cxx index f7d374ef7..67884c9c7 100644 --- a/sbncode/CAFMaker/FillTrue.cxx +++ b/sbncode/CAFMaker/FillTrue.cxx @@ -136,7 +136,7 @@ namespace caf { }//FillTrackTruth // Assumes truth matching and calo-points are filled - void FillTrackCaloTruth(const std::map> &id_to_ide_map, + void FillTrackCaloTruth(const std::map> &id_to_ide_map, const std::vector &mc_particles, const geo::GeometryCore& geometry, const geo::WireReadoutGeom& wireReadout, @@ -163,7 +163,7 @@ namespace caf { // Load the hits // match on the channel, which is unique - const std::vector &match_ides = id_to_ide_map.at(srtrack.truth.p.G4ID); + const std::vector &match_ides = id_to_ide_map.at(srtrack.truth.p.G4ID); std::map> chan_2_ides; for (auto const &ide_p: match_ides) { chan_2_ides[wireReadout.PlaneWireToChannel(ide_p.wire)].push_back(ide_p.ide); @@ -634,15 +634,15 @@ namespace caf { void FillTrueG4Particle(const simb::MCParticle &particle, const std::vector &active_volumes, const std::vector> &tpc_volumes, - const std::map> &id_to_ide_map, + const std::map> &id_to_ide_map, const std::map>> &id_to_truehit_map, const cheat::BackTrackerService &backtracker, const cheat::ParticleInventoryService &inventory_service, const std::vector> &neutrinos, caf::SRTrueParticle &srparticle) { - std::vector empty; - const std::vector &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty; + std::vector empty; + const std::vector &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty; std::vector> emptyHits; const std::vector> &particle_hits = id_to_truehit_map.count(particle.TrackId()) ? id_to_truehit_map.at(particle.TrackId()) : emptyHits; @@ -882,8 +882,8 @@ namespace caf { return ret; } - std::map> PrepSimChannels(const std::vector> &simchannels, const geo::WireReadoutGeom &wireReadout) { - std::map> ret; + std::map> PrepSimChannels(const std::vector> &simchannels, const geo::WireReadoutGeom &wireReadout) { + std::map> ret; for (const art::Ptr sc : simchannels) { // Lookup the wire of this channel diff --git a/sbncode/CAFMaker/FillTrue.h b/sbncode/CAFMaker/FillTrue.h index be5f48e69..6ebd4d69f 100644 --- a/sbncode/CAFMaker/FillTrue.h +++ b/sbncode/CAFMaker/FillTrue.h @@ -33,6 +33,8 @@ #include "sbnanaobj/StandardRecord/StandardRecord.h" #include "sbnanaobj/StandardRecord/SRMeVPrtl.h" +#include "sbncode/Calibration/TrackCaloSkimmer.h" + namespace caf { struct HitsEnergy { @@ -40,14 +42,6 @@ namespace caf float totE; }; - /// Ioniosation charge (`sim::IDE`) associated to its TPC readout information. - struct ReadoutIDE { - geo::WireID wire; ///< Wire on a given plane closest to the drift path of the charge. - unsigned short tick = 0; ///< Time tick at which the charge passes closest to the wire. - const sim::IDE *ide = nullptr; ///< Deposited charge information. - }; - - // Helpers caf::Wall_t GetWallCross( const geo::BoxBoundedGeo &volume, const TVector3 p0, @@ -81,7 +75,7 @@ namespace caf void FillTrueG4Particle(const simb::MCParticle &particle, const std::vector &active_volumes, const std::vector> &tpc_volumes, - const std::map> &id_to_ide_map, + const std::map> &id_to_ide_map, const std::map>> &id_to_truehit_map, const cheat::BackTrackerService &backtracker, const cheat::ParticleInventoryService &inventory_service, @@ -111,7 +105,7 @@ namespace caf caf::SRTrack& srtrack, bool allowEmpty = false); - void FillTrackCaloTruth(const std::map> &id_to_ide_map, + void FillTrackCaloTruth(const std::map> &id_to_ide_map, const std::vector &mc_particles, const geo::GeometryCore & geometry, const geo::WireReadoutGeom& wireReadout, @@ -141,7 +135,7 @@ namespace caf CLHEP::HepRandomEngine &rand, std::vector &srfakereco); - std::map> PrepSimChannels(const std::vector> &simchannels, const geo::WireReadoutGeom &wireReadout); + std::map> PrepSimChannels(const std::vector> &simchannels, const geo::WireReadoutGeom &wireReadout); std::map>> PrepTrueHits(const std::vector> &allHits, const detinfo::DetectorClocksData &clockData, const cheat::BackTrackerService &backtracker); std::map SetupIDHitEnergyMap(const std::vector> &allHits, const detinfo::DetectorClocksData &clockData, diff --git a/sbncode/Calibration/TrackCaloSkimmer.h b/sbncode/Calibration/TrackCaloSkimmer.h index e6e8f3aa3..a67ab9c85 100644 --- a/sbncode/Calibration/TrackCaloSkimmer.h +++ b/sbncode/Calibration/TrackCaloSkimmer.h @@ -64,12 +64,17 @@ #include "sbnobj/Common/CRT/CRTHitT0TaggingInfo.hh" #include "sbnobj/Common/CRT/CRTHitT0TaggingTruthInfo.hh" -#include "sbncode/CAFMaker/FillTrue.h" #include "ITCSSelectionTool.h" namespace sbn { class TrackCaloSkimmer; enum EDet {kNOTDEFINED, kSBND, kICARUS}; + + struct ReadoutIDE { + geo::WireID wire; ///< Wire on a given plane closest to the drift path of the charge. + unsigned short tick = 0; ///< Time tick at which the charge passes closest to the wire. + const sim::IDE *ide = nullptr; ///< Deposited charge information. + }; } class sbn::TrackCaloSkimmer : public art::EDAnalyzer { @@ -113,7 +118,6 @@ class sbn::TrackCaloSkimmer : public art::EDAnalyzer { int ID; }; - // Represents a "Snippet" of ADCs shared by a set of hits on a wire struct Snippet { geo::WireID wire; @@ -162,7 +166,7 @@ class sbn::TrackCaloSkimmer : public art::EDAnalyzer { const std::vector> &mcparticles, const std::vector &active_volumes, const std::vector> &tpc_volumes, - const std::map> id_to_ide_map, + const std::map> id_to_ide_map, const std::map>> id_to_truehit_map, const detinfo::DetectorPropertiesData &dprop, const geo::GeometryCore *geo, diff --git a/sbncode/Calibration/TrackCaloSkimmer_module.cc b/sbncode/Calibration/TrackCaloSkimmer_module.cc index 43f018ead..a3707eecd 100644 --- a/sbncode/Calibration/TrackCaloSkimmer_module.cc +++ b/sbncode/Calibration/TrackCaloSkimmer_module.cc @@ -289,7 +289,7 @@ void sbn::TrackCaloSkimmer::analyze(art::Event const& e) // Prep truth-to-reco-matching info // // Use helper functions from CAFMaker/FillTrue - std::map> id_to_ide_map; + std::map> id_to_ide_map; std::map>> id_to_truehit_map; const cheat::BackTrackerService *bt = NULL; @@ -545,14 +545,14 @@ geo::Point_t WireToTrajectoryPosition(const geo::Point_t &loc, const geo::TPCID sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, const std::vector &active_volumes, const std::vector> &tpc_volumes, - const std::map> &id_to_ide_map, + const std::map> &id_to_ide_map, const std::map>> &id_to_truehit_map, const detinfo::DetectorPropertiesData &dprop, const geo::GeometryCore *geo, const geo::WireReadoutGeom *wireReadout) { - std::vector empty; - const std::vector &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty; + std::vector empty; + const std::vector &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty; std::vector> emptyHits; const std::vector> &particle_hits = id_to_truehit_map.count(particle.TrackId()) ? id_to_truehit_map.at(particle.TrackId()) : emptyHits; @@ -1011,7 +1011,7 @@ void sbn::TrackCaloSkimmer::FillTrackTruth(const detinfo::DetectorClocksData &cl const std::vector> &mcparticles, const std::vector &active_volumes, const std::vector> &tpc_volumes, - const std::map> id_to_ide_map, + const std::map> id_to_ide_map, const std::map>> id_to_truehit_map, const detinfo::DetectorPropertiesData &dprop, const geo::GeometryCore *geo, From 1e0a36f38b2960bd6cd16ad95704ddbb695de467 Mon Sep 17 00:00:00 2001 From: sungbinoh Date: Wed, 3 Sep 2025 10:33:10 -0500 Subject: [PATCH 11/13] updating phi to rad from degree, and modifying some lines in CMakeLists.txt and FillReco.h headers --- sbncode/CAFMaker/CMakeLists.txt | 3 ++- sbncode/CAFMaker/FillReco.cxx | 2 +- sbncode/CAFMaker/FillReco.h | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/sbncode/CAFMaker/CMakeLists.txt b/sbncode/CAFMaker/CMakeLists.txt index 4ef86971c..2326f8dee 100644 --- a/sbncode/CAFMaker/CMakeLists.txt +++ b/sbncode/CAFMaker/CMakeLists.txt @@ -29,7 +29,7 @@ art_make_library( LIBRARY_NAME sbncode_CAFMaker larcorealg::Geometry larcore::Geometry_Geometry_service larcorealg::GeoAlgo - lardata::DetectorInfoServices_DetectorPropertiesServiceStandard_service + lardataalg::headers larsim::MCCheater_BackTrackerService_service nusimdata::SimulationBase larsim::MCCheater_ParticleInventoryService_service @@ -66,6 +66,7 @@ cet_build_plugin ( CAFMaker art::module hep_concurrency::hep_concurrency lardataobj::RecoBase nurandom::RandomUtils_NuRandomService_service + lardata::DetectorPropertiesService lardata::DetectorClocksService BASENAME_ONLY ) diff --git a/sbncode/CAFMaker/FillReco.cxx b/sbncode/CAFMaker/FillReco.cxx index 92691d596..f4eae4b18 100644 --- a/sbncode/CAFMaker/FillReco.cxx +++ b/sbncode/CAFMaker/FillReco.cxx @@ -877,7 +877,7 @@ namespace caf unsigned traj_point_index = thms.at(i_hit)->Index(); if (track.HasValidPoint(traj_point_index)) { float costh_drift = track.DirectionAtPoint(traj_point_index).X(); - float phi = acos(abs(costh_drift)) * 180. / M_PI; + float phi = acos(abs(costh_drift)); float efield = GetEfield(dprop, track.LocationAtPoint(traj_point_index)); p.efield = efield; p.phi = phi; diff --git a/sbncode/CAFMaker/FillReco.h b/sbncode/CAFMaker/FillReco.h index 8551b6a44..c1b2d7f65 100644 --- a/sbncode/CAFMaker/FillReco.h +++ b/sbncode/CAFMaker/FillReco.h @@ -5,7 +5,7 @@ #include "art/Framework/Services/Registry/ServiceHandle.h" #include "larsim/MCCheater/ParticleInventoryService.h" -#include "lardataalg/DetectorInfo/DetectorPropertiesStandard.h" +#include "lardataalg/DetectorInfo/DetectorPropertiesData.h" // LArSoft includes #include "larcorealg/Geometry/fwd.h" From 46e144f5f4d704c87e49c8b661bd984dc5e7ba1e Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Wed, 3 Sep 2025 11:34:24 -0500 Subject: [PATCH 12/13] Refactor CAF helper functions utilized in CaloSkimmer into RecoUtils. --- sbncode/CAFMaker/CAFMaker_module.cc | 6 +- sbncode/CAFMaker/FillTrue.cxx | 160 +----------------- sbncode/CAFMaker/FillTrue.h | 12 +- sbncode/CAFMaker/RecoUtils/RecoUtils.cc | 148 ++++++++++++++++ sbncode/CAFMaker/RecoUtils/RecoUtils.h | 24 +++ sbncode/Calibration/CMakeLists.txt | 1 + sbncode/Calibration/TrackCaloSkimmer.h | 11 +- .../Calibration/TrackCaloSkimmer_module.cc | 19 +-- 8 files changed, 193 insertions(+), 188 deletions(-) diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index 0c11332ab..0984ed944 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -64,7 +64,7 @@ #include "art/Framework/Services/System/TriggerNamesService.h" #include "nurandom/RandomUtils/NuRandomService.h" #include "lardata/DetectorInfoServices/DetectorClocksService.h" -#include "lardataalg/DetectorInfo/DetectorPropertiesStandard.h" +#include "lardataalg/DetectorInfo/DetectorPropertiesData.h" #include "lardata/DetectorInfoServices/DetectorPropertiesService.h" #include "larcore/CoreUtils/ServiceUtil.h" #include "larevt/SpaceCharge/SpaceCharge.h" @@ -1394,8 +1394,8 @@ void CAFMaker::produce(art::Event& evt) noexcept { if ( !isRealData ) { art::ServiceHandle bt_serv; - id_to_ide_map = PrepSimChannels(simchannels, wireReadout); - id_to_truehit_map = PrepTrueHits(hits, clock_data, *bt_serv); + id_to_ide_map = sbn::PrepSimChannels(simchannels, wireReadout); + id_to_truehit_map = sbn::PrepTrueHits(hits, clock_data, *bt_serv); id_to_hit_energy_map = SetupIDHitEnergyMap(hits, clock_data, *bt_serv); } diff --git a/sbncode/CAFMaker/FillTrue.cxx b/sbncode/CAFMaker/FillTrue.cxx index 67884c9c7..aeeb0646d 100644 --- a/sbncode/CAFMaker/FillTrue.cxx +++ b/sbncode/CAFMaker/FillTrue.cxx @@ -5,7 +5,6 @@ #include "larcorealg/Geometry/WireReadoutGeom.h" #include "larevt/SpaceCharge/SpaceCharge.h" #include "larcore/CoreUtils/ServiceUtil.h" -#include "RecoUtils/RecoUtils.h" #include "CLHEP/Random/RandGauss.h" @@ -701,7 +700,7 @@ namespace caf { } // get the wall if (entry_point > 0) { - srparticle.wallin = GetWallCross(active_volumes.at(cryostat_index), particle.Position(entry_point).Vect(), particle.Position(entry_point-1).Vect()); + srparticle.wallin = sbn::GetWallCross(active_volumes.at(cryostat_index), particle.Position(entry_point).Vect(), particle.Position(entry_point-1).Vect()); } int exit_point = -1; @@ -778,7 +777,7 @@ namespace caf { exit_point = particle.NumberTrajectoryPoints() - 1; } if (exit_point >= 0 && ((unsigned)exit_point) < particle.NumberTrajectoryPoints() - 1) { - srparticle.wallout = GetWallCross(active_volumes.at(cryostat_index), particle.Position(exit_point).Vect(), particle.Position(exit_point+1).Vect()); + srparticle.wallout = sbn::GetWallCross(active_volumes.at(cryostat_index), particle.Position(exit_point).Vect(), particle.Position(exit_point+1).Vect()); } // other truth information @@ -799,8 +798,8 @@ namespace caf { srparticle.endp = (exit_point >= 0) ? particle.Momentum(exit_point).Vect() : TVector3(-9999, -9999, -9999); srparticle.endE = (exit_point >= 0) ? particle.Momentum(exit_point).E() : -9999.; - srparticle.start_process = GetG4ProcessID(particle.Process()); - srparticle.end_process = GetG4ProcessID(particle.EndProcess()); + srparticle.start_process = sbn::GetG4ProcessID(particle.Process()); + srparticle.end_process = sbn::GetG4ProcessID(particle.EndProcess()); srparticle.G4ID = particle.TrackId(); srparticle.parent = particle.Mother(); @@ -871,37 +870,6 @@ namespace caf { return ret; } - std::map>> PrepTrueHits(const std::vector> &allHits, - const detinfo::DetectorClocksData &clockData, const cheat::BackTrackerService &backtracker) { - std::map>> ret; - for (const art::Ptr h: allHits) { - for (int ID: backtracker.HitToTrackIds(clockData, *h)) { - ret[abs(ID)].push_back(h); - } - } - return ret; - } - - std::map> PrepSimChannels(const std::vector> &simchannels, const geo::WireReadoutGeom &wireReadout) { - std::map> ret; - - for (const art::Ptr sc : simchannels) { - // Lookup the wire of this channel - raw::ChannelID_t channel = sc->Channel(); - std::vector maybewire = wireReadout.ChannelToWire(channel); - geo::WireID thisWire; // Default constructor makes invalid wire - if (maybewire.size()) thisWire = maybewire[0]; - - for (const auto &item : sc->TDCIDEMap()) { - for (const sim::IDE &ide: item.second) { - // indexing initializes empty vector - ret[abs(ide.trackID)].push_back({thisWire, item.first, &ide}); - } - } - } - return ret; - } - } // end namespace @@ -1253,126 +1221,6 @@ bool FRFillNueCC(const simb::MCTruth &mctruth, return true; } -caf::Wall_t caf::GetWallCross(const geo::BoxBoundedGeo &volume, const TVector3 p0, const TVector3 p1) { - TVector3 direction = (p1 - p0) * ( 1. / (p1 - p0).Mag()); - std::vector intersections = volume.GetIntersections(p0, direction); - - assert(intersections.size() == 2); - - // get the intersection point closer to p0 - int intersection_i = ((intersections[0] - p0).Mag() < (intersections[1] - p0).Mag()) ? 0 : 1; - - double eps = 1e-3; - if (abs(intersections[intersection_i].X() - volume.MinX()) < eps) { - //std::cout << "Left\n"; - return caf::kWallLeft; - } - else if (abs(intersections[intersection_i].X() - volume.MaxX()) < eps) { - //std::cout << "Right\n"; - return caf::kWallRight; - } - else if (abs(intersections[intersection_i].Y() - volume.MinY()) < eps) { - //std::cout << "Bottom\n"; - return caf::kWallBottom; - } - else if (abs(intersections[intersection_i].Y() - volume.MaxY()) < eps) { - //std::cout << "Top\n"; - return caf::kWallTop; - } - else if (abs(intersections[intersection_i].Z() - volume.MinZ()) < eps) { - //std::cout << "Front\n"; - return caf::kWallFront; - } - else if (abs(intersections[intersection_i].Z() - volume.MaxZ()) < eps) { - //std::cout << "Back\n"; - return caf::kWallBack; - } - else assert(false); - //std::cout << "None\n"; - - return caf::kWallNone; -}//GetWallCross - -//------------------------------------------ - -caf::g4_process_ caf::GetG4ProcessID(const std::string &process_name) { -#define MATCH_PROCESS(name) if (process_name == #name) {return caf::kG4 ## name;} -#define MATCH_PROCESS_NAMED(strname, id) if (process_name == #strname) {return caf::kG4 ## id;} - MATCH_PROCESS(primary) - MATCH_PROCESS(CoupledTransportation) - MATCH_PROCESS(FastScintillation) - MATCH_PROCESS(Decay) - MATCH_PROCESS(anti_neutronInelastic) - MATCH_PROCESS(neutronInelastic) - MATCH_PROCESS(anti_protonInelastic) - MATCH_PROCESS(protonInelastic) - MATCH_PROCESS(hadInelastic) - MATCH_PROCESS_NAMED(kaon+Inelastic, kaonpInelastic) - MATCH_PROCESS_NAMED(kaon-Inelastic, kaonmInelastic) - MATCH_PROCESS_NAMED(kaon+Inelastic, kaonpInelastic) - MATCH_PROCESS_NAMED(kaon-Inelastic, kaonmInelastic) - MATCH_PROCESS_NAMED(sigma+Inelastic, sigmapInelastic) - MATCH_PROCESS_NAMED(sigma-Inelastic, sigmamInelastic) - MATCH_PROCESS_NAMED(pi+Inelastic, pipInelastic) - MATCH_PROCESS_NAMED(pi-Inelastic, pimInelastic) - MATCH_PROCESS_NAMED(xi+Inelastic, xipInelastic) - MATCH_PROCESS_NAMED(xi-Inelastic, ximInelastic) - MATCH_PROCESS(kaon0LInelastic) - MATCH_PROCESS(kaon0SInelastic) - MATCH_PROCESS(lambdaInelastic) - MATCH_PROCESS_NAMED(anti-lambdaInelastic, anti_lambdaInelastic) - MATCH_PROCESS(He3Inelastic) - MATCH_PROCESS(ionInelastic) - MATCH_PROCESS(xi0Inelastic) - MATCH_PROCESS(alphaInelastic) - MATCH_PROCESS(tInelastic) - MATCH_PROCESS(dInelastic) - MATCH_PROCESS(anti_neutronElastic) - MATCH_PROCESS(neutronElastic) - MATCH_PROCESS(anti_protonElastic) - MATCH_PROCESS(protonElastic) - MATCH_PROCESS(hadElastic) - MATCH_PROCESS_NAMED(kaon+Elastic, kaonpElastic) - MATCH_PROCESS_NAMED(kaon-Elastic, kaonmElastic) - MATCH_PROCESS_NAMED(pi+Elastic, pipElastic) - MATCH_PROCESS_NAMED(pi-Elastic, pimElastic) - MATCH_PROCESS(conv) - MATCH_PROCESS(phot) - MATCH_PROCESS(annihil) - MATCH_PROCESS(nCapture) - MATCH_PROCESS(nKiller) - MATCH_PROCESS(muMinusCaptureAtRest) - MATCH_PROCESS(muIoni) - MATCH_PROCESS(eBrem) - MATCH_PROCESS(CoulombScat) - MATCH_PROCESS(hBertiniCaptureAtRest) - MATCH_PROCESS(hFritiofCaptureAtRest) - MATCH_PROCESS(photonNuclear) - MATCH_PROCESS(muonNuclear) - MATCH_PROCESS(electronNuclear) - MATCH_PROCESS(positronNuclear) - MATCH_PROCESS(compt) - MATCH_PROCESS(eIoni) - MATCH_PROCESS(muBrems) - MATCH_PROCESS(hIoni) - MATCH_PROCESS(ionIoni) - MATCH_PROCESS(hBrems) - MATCH_PROCESS(muPairProd) - MATCH_PROCESS(hPairProd) - MATCH_PROCESS(LArVoxelReadoutScoringProcess) - MATCH_PROCESS(Transportation) - MATCH_PROCESS(msc) - MATCH_PROCESS(StepLimiter) - MATCH_PROCESS(RadioactiveDecayBase) - std::cerr << "Error: Process name with no match (" << process_name << ")\n"; - assert(false); - return caf::kG4UNKNOWN; // unreachable in debug mode -#undef MATCH_PROCESS -#undef MATCH_PROCESS_NAMED - -}//GetG4ProcessID -//------------------------------------------- - float ContainedLength(const TVector3 &v0, const TVector3 &v1, const std::vector &boxes) { static const geoalgo::GeoAlgo algo; diff --git a/sbncode/CAFMaker/FillTrue.h b/sbncode/CAFMaker/FillTrue.h index 6ebd4d69f..d4134d55c 100644 --- a/sbncode/CAFMaker/FillTrue.h +++ b/sbncode/CAFMaker/FillTrue.h @@ -33,7 +33,7 @@ #include "sbnanaobj/StandardRecord/StandardRecord.h" #include "sbnanaobj/StandardRecord/SRMeVPrtl.h" -#include "sbncode/Calibration/TrackCaloSkimmer.h" +#include "RecoUtils/RecoUtils.h" namespace caf { @@ -42,13 +42,6 @@ namespace caf float totE; }; - // Helpers - caf::Wall_t GetWallCross( const geo::BoxBoundedGeo &volume, - const TVector3 p0, - const TVector3 p1); - - caf::g4_process_ GetG4ProcessID(const std::string &name); - void FillSRGlobal(const sbn::evwgh::EventWeightParameterSet& pset, caf::SRGlobal& srglobal, std::map& weightPSetIndex); @@ -135,9 +128,6 @@ namespace caf CLHEP::HepRandomEngine &rand, std::vector &srfakereco); - std::map> PrepSimChannels(const std::vector> &simchannels, const geo::WireReadoutGeom &wireReadout); - std::map>> PrepTrueHits(const std::vector> &allHits, - const detinfo::DetectorClocksData &clockData, const cheat::BackTrackerService &backtracker); std::map SetupIDHitEnergyMap(const std::vector> &allHits, const detinfo::DetectorClocksData &clockData, const cheat::BackTrackerService &backtracker); diff --git a/sbncode/CAFMaker/RecoUtils/RecoUtils.cc b/sbncode/CAFMaker/RecoUtils/RecoUtils.cc index d8e9f210f..a2cdc8922 100644 --- a/sbncode/CAFMaker/RecoUtils/RecoUtils.cc +++ b/sbncode/CAFMaker/RecoUtils/RecoUtils.cc @@ -107,3 +107,151 @@ int CAFRecoUtils::GetShowerPrimary(const int g4ID) return primary_id; } + +std::map> sbn::PrepSimChannels(const std::vector> &simchannels, const geo::WireReadoutGeom &wireReadout) { + std::map> ret; + + for (const art::Ptr sc : simchannels) { + // Lookup the wire of this channel + raw::ChannelID_t channel = sc->Channel(); + std::vector maybewire = wireReadout.ChannelToWire(channel); + geo::WireID thisWire; // Default constructor makes invalid wire + if (maybewire.size()) thisWire = maybewire[0]; + + for (const auto &item : sc->TDCIDEMap()) { + for (const sim::IDE &ide: item.second) { + // indexing initializes empty vector + ret[abs(ide.trackID)].push_back({thisWire, item.first, &ide}); + } + } + } + return ret; +} + +std::map>> sbn::PrepTrueHits(const std::vector> &allHits, + const detinfo::DetectorClocksData &clockData, const cheat::BackTrackerService &backtracker) { + std::map>> ret; + for (const art::Ptr h: allHits) { + for (int ID: backtracker.HitToTrackIds(clockData, *h)) { + ret[abs(ID)].push_back(h); + } + } + return ret; +} + +caf::Wall_t sbn::GetWallCross(const geo::BoxBoundedGeo &volume, const TVector3 p0, const TVector3 p1) { + TVector3 direction = (p1 - p0) * ( 1. / (p1 - p0).Mag()); + std::vector intersections = volume.GetIntersections(p0, direction); + + assert(intersections.size() == 2); + + // get the intersection point closer to p0 + int intersection_i = ((intersections[0] - p0).Mag() < (intersections[1] - p0).Mag()) ? 0 : 1; + + double eps = 1e-3; + if (abs(intersections[intersection_i].X() - volume.MinX()) < eps) { + //std::cout << "Left\n"; + return caf::kWallLeft; + } + else if (abs(intersections[intersection_i].X() - volume.MaxX()) < eps) { + //std::cout << "Right\n"; + return caf::kWallRight; + } + else if (abs(intersections[intersection_i].Y() - volume.MinY()) < eps) { + //std::cout << "Bottom\n"; + return caf::kWallBottom; + } + else if (abs(intersections[intersection_i].Y() - volume.MaxY()) < eps) { + //std::cout << "Top\n"; + return caf::kWallTop; + } + else if (abs(intersections[intersection_i].Z() - volume.MinZ()) < eps) { + //std::cout << "Front\n"; + return caf::kWallFront; + } + else if (abs(intersections[intersection_i].Z() - volume.MaxZ()) < eps) { + //std::cout << "Back\n"; + return caf::kWallBack; + } + else assert(false); + //std::cout << "None\n"; + + return caf::kWallNone; +}//GetWallCross + +//------------------------------------------ +caf::g4_process_ sbn::GetG4ProcessID(const std::string &process_name) { +#define MATCH_PROCESS(name) if (process_name == #name) {return caf::kG4 ## name;} +#define MATCH_PROCESS_NAMED(strname, id) if (process_name == #strname) {return caf::kG4 ## id;} + MATCH_PROCESS(primary) + MATCH_PROCESS(CoupledTransportation) + MATCH_PROCESS(FastScintillation) + MATCH_PROCESS(Decay) + MATCH_PROCESS(anti_neutronInelastic) + MATCH_PROCESS(neutronInelastic) + MATCH_PROCESS(anti_protonInelastic) + MATCH_PROCESS(protonInelastic) + MATCH_PROCESS(hadInelastic) + MATCH_PROCESS_NAMED(kaon+Inelastic, kaonpInelastic) + MATCH_PROCESS_NAMED(kaon-Inelastic, kaonmInelastic) + MATCH_PROCESS_NAMED(kaon+Inelastic, kaonpInelastic) + MATCH_PROCESS_NAMED(kaon-Inelastic, kaonmInelastic) + MATCH_PROCESS_NAMED(sigma+Inelastic, sigmapInelastic) + MATCH_PROCESS_NAMED(sigma-Inelastic, sigmamInelastic) + MATCH_PROCESS_NAMED(pi+Inelastic, pipInelastic) + MATCH_PROCESS_NAMED(pi-Inelastic, pimInelastic) + MATCH_PROCESS_NAMED(xi+Inelastic, xipInelastic) + MATCH_PROCESS_NAMED(xi-Inelastic, ximInelastic) + MATCH_PROCESS(kaon0LInelastic) + MATCH_PROCESS(kaon0SInelastic) + MATCH_PROCESS(lambdaInelastic) + MATCH_PROCESS_NAMED(anti-lambdaInelastic, anti_lambdaInelastic) + MATCH_PROCESS(He3Inelastic) + MATCH_PROCESS(ionInelastic) + MATCH_PROCESS(xi0Inelastic) + MATCH_PROCESS(alphaInelastic) + MATCH_PROCESS(tInelastic) + MATCH_PROCESS(dInelastic) + MATCH_PROCESS(anti_neutronElastic) + MATCH_PROCESS(neutronElastic) + MATCH_PROCESS(anti_protonElastic) + MATCH_PROCESS(protonElastic) + MATCH_PROCESS(hadElastic) + MATCH_PROCESS_NAMED(kaon+Elastic, kaonpElastic) + MATCH_PROCESS_NAMED(kaon-Elastic, kaonmElastic) + MATCH_PROCESS_NAMED(pi+Elastic, pipElastic) + MATCH_PROCESS_NAMED(pi-Elastic, pimElastic) + MATCH_PROCESS(conv) + MATCH_PROCESS(phot) + MATCH_PROCESS(annihil) + MATCH_PROCESS(nCapture) + MATCH_PROCESS(nKiller) + MATCH_PROCESS(muMinusCaptureAtRest) + MATCH_PROCESS(muIoni) + MATCH_PROCESS(eBrem) + MATCH_PROCESS(CoulombScat) + MATCH_PROCESS(hBertiniCaptureAtRest) + MATCH_PROCESS(hFritiofCaptureAtRest) + MATCH_PROCESS(photonNuclear) + MATCH_PROCESS(muonNuclear) + MATCH_PROCESS(electronNuclear) + MATCH_PROCESS(positronNuclear) + MATCH_PROCESS(compt) + MATCH_PROCESS(eIoni) + MATCH_PROCESS(muBrems) + MATCH_PROCESS(hIoni) + MATCH_PROCESS(ionIoni) + MATCH_PROCESS(hBrems) + MATCH_PROCESS(muPairProd) + MATCH_PROCESS(hPairProd) + MATCH_PROCESS(LArVoxelReadoutScoringProcess) + MATCH_PROCESS(Transportation) + MATCH_PROCESS(msc) + MATCH_PROCESS(StepLimiter) + MATCH_PROCESS(RadioactiveDecayBase) + std::cerr << "Error: Process name with no match (" << process_name << ")\n"; + assert(false); + return caf::kG4UNKNOWN; // unreachable in debug mode +#undef MATCH_PROCESS +#undef MATCH_PROCESS_NAMED +}//GetG4ProcessID diff --git a/sbncode/CAFMaker/RecoUtils/RecoUtils.h b/sbncode/CAFMaker/RecoUtils/RecoUtils.h index a99121d8e..0b43760a4 100644 --- a/sbncode/CAFMaker/RecoUtils/RecoUtils.h +++ b/sbncode/CAFMaker/RecoUtils/RecoUtils.h @@ -24,11 +24,13 @@ #include "lardataobj/RecoBase/Hit.h" #include "lardataobj/RecoBase/Track.h" #include "lardataobj/RecoBase/Shower.h" +#include "sbnanaobj/StandardRecord/SREnums.h" //#include "lardataobj/AnalysisBase/MVAPIDResult.h" //#include "lardataobj/AnalysisBase/ParticleID.h" #include "larsim/MCCheater/BackTrackerService.h" #include "larsim/MCCheater/ParticleInventoryService.h" #include "larcore/Geometry/Geometry.h" +#include "larcorealg/Geometry/WireReadoutGeom.h" // c++ #include @@ -37,6 +39,25 @@ // ROOT #include "TTree.h" +namespace sim { + class IDE; +} + +namespace sbn { + struct ReadoutIDE { + geo::WireID wire; ///< Wire on a given plane closest to the drift path of the charge. + unsigned short tick = 0; ///< Time tick at which the charge passes closest to the wire. + const sim::IDE *ide = nullptr; ///< Deposited charge information. + }; + + // Helpers + std::map> PrepSimChannels(const std::vector> &simchannels, const geo::WireReadoutGeom &wireReadout); + std::map>> PrepTrueHits(const std::vector> &allHits, + const detinfo::DetectorClocksData &clockData, const cheat::BackTrackerService &backtracker); + caf::Wall_t GetWallCross( const geo::BoxBoundedGeo &volume, const TVector3 p0, const TVector3 p1); + caf::g4_process_ GetG4ProcessID(const std::string &name); +} + namespace CAFRecoUtils{ std::vector> AllTrueParticleIDEnergyMatches(const detinfo::DetectorClocksData &clockData, const std::vector >& hits, bool rollup_unsaved_ids=1); @@ -47,4 +68,7 @@ namespace CAFRecoUtils{ int GetShowerPrimary(const int g4ID); } + + + #endif diff --git a/sbncode/Calibration/CMakeLists.txt b/sbncode/Calibration/CMakeLists.txt index 98059edc5..32dfd809b 100644 --- a/sbncode/Calibration/CMakeLists.txt +++ b/sbncode/Calibration/CMakeLists.txt @@ -25,6 +25,7 @@ cet_build_plugin( TrackCaloSkimmer art::module sbncode_CAFMaker sbnobj::Common_Calibration_dict larevt::SpaceCharge + caf_RecoUtils ) cet_build_plugin(TrackCaloSkimmerSelectStoppingTrack art::tool diff --git a/sbncode/Calibration/TrackCaloSkimmer.h b/sbncode/Calibration/TrackCaloSkimmer.h index a67ab9c85..ead9571ae 100644 --- a/sbncode/Calibration/TrackCaloSkimmer.h +++ b/sbncode/Calibration/TrackCaloSkimmer.h @@ -51,8 +51,8 @@ #include "larevt/SpaceCharge/SpaceCharge.h" #include "larevt/SpaceChargeServices/SpaceChargeService.h" -#include "lardataalg/DetectorInfo/DetectorPropertiesStandard.h" #include "lardata/DetectorInfoServices/DetectorPropertiesService.h" +#include "lardataalg/DetectorInfo/DetectorPropertiesData.h" #include "lardata/DetectorInfoServices/DetectorClocksService.h" #include "larcorealg/Geometry/fwd.h" @@ -66,15 +66,12 @@ #include "ITCSSelectionTool.h" +// Useful functions +#include "sbncode/CAFMaker/RecoUtils/RecoUtils.h" + namespace sbn { class TrackCaloSkimmer; enum EDet {kNOTDEFINED, kSBND, kICARUS}; - - struct ReadoutIDE { - geo::WireID wire; ///< Wire on a given plane closest to the drift path of the charge. - unsigned short tick = 0; ///< Time tick at which the charge passes closest to the wire. - const sim::IDE *ide = nullptr; ///< Deposited charge information. - }; } class sbn::TrackCaloSkimmer : public art::EDAnalyzer { diff --git a/sbncode/Calibration/TrackCaloSkimmer_module.cc b/sbncode/Calibration/TrackCaloSkimmer_module.cc index a3707eecd..ff73432be 100644 --- a/sbncode/Calibration/TrackCaloSkimmer_module.cc +++ b/sbncode/Calibration/TrackCaloSkimmer_module.cc @@ -12,13 +12,10 @@ #include "TrackCaloSkimmer.h" #include "sbnobj/SBND/CRT/CRTTrack.hh" +#include "sbnanaobj/StandardRecord/SREnums.h" #include "art/Utilities/make_tool.h" -// Useful functions -#include "sbncode/CAFMaker/FillTrue.h" -#include "sbncode/CAFMaker/RecoUtils/RecoUtils.h" - #include "larcore/Geometry/WireReadout.h" #include "larcore/Geometry/Geometry.h" #include "larcore/CoreUtils/ServiceUtil.h" @@ -288,15 +285,15 @@ void sbn::TrackCaloSkimmer::analyze(art::Event const& e) // Prep truth-to-reco-matching info // - // Use helper functions from CAFMaker/FillTrue + // Use helper functions from CAFMaker/RecoUtils std::map> id_to_ide_map; std::map>> id_to_truehit_map; const cheat::BackTrackerService *bt = NULL; if (simchannels.size()) { art::ServiceHandle bt_serv; - id_to_ide_map = caf::PrepSimChannels(simchannels, *wireReadout); - id_to_truehit_map = caf::PrepTrueHits(allHits, clock_data, *bt_serv.get()); + id_to_ide_map = sbn::PrepSimChannels(simchannels, *wireReadout); + id_to_truehit_map = sbn::PrepTrueHits(allHits, clock_data, *bt_serv.get()); bt = bt_serv.get(); } @@ -622,7 +619,7 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, // get the wall if (entry_point > 0) { - trueparticle.wallin = (int)caf::GetWallCross(active_volumes.at(cryostat_index), particle.Position(entry_point).Vect(), particle.Position(entry_point-1).Vect()); + trueparticle.wallin = (int)sbn::GetWallCross(active_volumes.at(cryostat_index), particle.Position(entry_point).Vect(), particle.Position(entry_point-1).Vect()); } int exit_point = -1; @@ -690,7 +687,7 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, exit_point = particle.NumberTrajectoryPoints() - 1; } if (exit_point >= 0 && ((unsigned)exit_point) < particle.NumberTrajectoryPoints() - 1) { - trueparticle.wallout = (int)caf::GetWallCross(active_volumes.at(cryostat_index), particle.Position(exit_point).Vect(), particle.Position(exit_point+1).Vect()); + trueparticle.wallout = (int)sbn::GetWallCross(active_volumes.at(cryostat_index), particle.Position(exit_point).Vect(), particle.Position(exit_point+1).Vect()); } // other truth information @@ -711,8 +708,8 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle, trueparticle.endp = ConvertTVector((exit_point >= 0) ? particle.Momentum(exit_point).Vect() : TVector3(-9999, -9999, -9999)); trueparticle.endE = (exit_point >= 0) ? particle.Momentum(exit_point).E() : -9999.; - trueparticle.start_process = (int)caf::GetG4ProcessID(particle.Process()); - trueparticle.end_process = (int)caf::GetG4ProcessID(particle.EndProcess()); + trueparticle.start_process = (int)sbn::GetG4ProcessID(particle.Process()); + trueparticle.end_process = (int)sbn::GetG4ProcessID(particle.EndProcess()); trueparticle.G4ID = particle.TrackId(); trueparticle.parent = particle.Mother(); From 6df0c309cb26d156ae9dcd6c83b62592755a395f Mon Sep 17 00:00:00 2001 From: sungbinoh Date: Wed, 3 Sep 2025 17:31:40 -0500 Subject: [PATCH 13/13] refer to *sce in CAFMaker_module.cc from FillReco, indentation fixes --- sbncode/CAFMaker/CAFMaker_module.cc | 2 +- sbncode/CAFMaker/FillReco.cxx | 19 ++++++++----------- sbncode/CAFMaker/FillReco.h | 4 ++++ sbncode/Calibration/CMakeLists.txt | 2 +- sbncode/Calibration/TrackCaloSkimmer.h | 2 +- 5 files changed, 15 insertions(+), 14 deletions(-) diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index 0984ed944..770f546af 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -2243,7 +2243,7 @@ void CAFMaker::produce(art::Event& evt) noexcept { FillTrackCalo(fmCalo.at(iPart), *thisTrack[0], fmTrackHit.at(iPart), fmTrackHit.data(iPart), (fParams.FillHitsNeutrinoSlices() && NeutrinoSlice) || fParams.FillHitsAllSlices(), fParams.TrackHitFillRRStartCut(), fParams.TrackHitFillRREndCut(), - dprop, trk); + dprop, *sce, trk); } if (fmCRTHitMatch.isValid() && fDet == kICARUS) { diff --git a/sbncode/CAFMaker/FillReco.cxx b/sbncode/CAFMaker/FillReco.cxx index f4eae4b18..5089f8631 100644 --- a/sbncode/CAFMaker/FillReco.cxx +++ b/sbncode/CAFMaker/FillReco.cxx @@ -10,10 +10,6 @@ #include "FillReco.h" #include "RecoUtils/RecoUtils.h" -#include "larevt/SpaceCharge/SpaceCharge.h" -#include "larevt/SpaceChargeServices/SpaceChargeService.h" -#include "larcore/CoreUtils/ServiceUtil.h" - namespace caf { const float ng_filter_cut = 0.5; @@ -799,13 +795,12 @@ namespace caf } // Helper function: get the e field - double GetEfield(const detinfo::DetectorPropertiesData& dprop, const geo::Point_t loc) { - auto const* sce = lar::providerFrom(); + double GetEfield(const detinfo::DetectorPropertiesData& dprop, spacecharge::SpaceCharge const& sce, const geo::Point_t& loc) { double EField = dprop.Efield(); - if (sce->EnableSimEfieldSCE()) { + if (sce.EnableSimEfieldSCE()) { // Gets fractional E field distortions w.r.t. the nominal field on x-axis - geo::Vector_t EFieldOffsets = sce->GetEfieldOffsets(loc); + geo::Vector_t EFieldOffsets = sce.GetEfieldOffsets(loc); // Add 1 in X direction as this is the direction of the drift field, not caring if it is +x or -x direction, since we only want |E| EFieldOffsets += geo::Vector_t{1, 0, 0}; // Convert to Absolute E Field from relative @@ -822,6 +817,7 @@ namespace caf const std::vector& thms, bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend, const detinfo::DetectorPropertiesData &dprop, + const spacecharge::SpaceCharge &sce, caf::SRTrackCalo &srcalo) { // Collect info from Calorimetry @@ -875,10 +871,10 @@ namespace caf // // Use this to get the (SCE corrected) efield and the angle to the drift direction unsigned traj_point_index = thms.at(i_hit)->Index(); - if (track.HasValidPoint(traj_point_index)) { + if (track.HasValidPoint(traj_point_index)) { float costh_drift = track.DirectionAtPoint(traj_point_index).X(); float phi = acos(abs(costh_drift)); - float efield = GetEfield(dprop, track.LocationAtPoint(traj_point_index)); + float efield = GetEfield(dprop, sce, track.LocationAtPoint(traj_point_index)); p.efield = efield; p.phi = phi; } @@ -940,6 +936,7 @@ namespace caf const std::vector& thms, bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend, const detinfo::DetectorPropertiesData &dprop, + const spacecharge::SpaceCharge &sce, caf::SRTrack& srtrack, bool allowEmpty) { @@ -952,7 +949,7 @@ namespace caf if (calo.PlaneID()) { unsigned plane_id = calo.PlaneID().Plane; assert(plane_id < 3); - FillTrackPlaneCalo(calo, track, hits, thms, fill_calo_points, fillhit_rrstart, fillhit_rrend, dprop, srtrack.calo[plane_id]); + FillTrackPlaneCalo(calo, track, hits, thms, fill_calo_points, fillhit_rrstart, fillhit_rrend, dprop, sce, srtrack.calo[plane_id]); } } diff --git a/sbncode/CAFMaker/FillReco.h b/sbncode/CAFMaker/FillReco.h index c1b2d7f65..38c5d5f6b 100644 --- a/sbncode/CAFMaker/FillReco.h +++ b/sbncode/CAFMaker/FillReco.h @@ -10,6 +10,8 @@ // LArSoft includes #include "larcorealg/Geometry/fwd.h" +#include "larevt/SpaceCharge/SpaceCharge.h" + #include "lardataobj/RecoBase/PFParticle.h" #include "lardataobj/RecoBase/Shower.h" #include "lardataobj/RecoBase/Slice.h" @@ -210,6 +212,7 @@ namespace caf const std::vector& thms, bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend, const detinfo::DetectorPropertiesData &dprop, + spacecharge::SpaceCharge const& sce, caf::SRTrackCalo &srcalo); void FillTrackScatterClosestApproach(const art::Ptr closestapproach, @@ -230,6 +233,7 @@ namespace caf const std::vector& thms, bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend, const detinfo::DetectorPropertiesData &dprop, + spacecharge::SpaceCharge const& sce, caf::SRTrack& srtrack, bool allowEmpty = false); diff --git a/sbncode/Calibration/CMakeLists.txt b/sbncode/Calibration/CMakeLists.txt index 32dfd809b..759c5c9f7 100644 --- a/sbncode/Calibration/CMakeLists.txt +++ b/sbncode/Calibration/CMakeLists.txt @@ -25,7 +25,7 @@ cet_build_plugin( TrackCaloSkimmer art::module sbncode_CAFMaker sbnobj::Common_Calibration_dict larevt::SpaceCharge - caf_RecoUtils + caf_RecoUtils ) cet_build_plugin(TrackCaloSkimmerSelectStoppingTrack art::tool diff --git a/sbncode/Calibration/TrackCaloSkimmer.h b/sbncode/Calibration/TrackCaloSkimmer.h index ead9571ae..1096c056d 100644 --- a/sbncode/Calibration/TrackCaloSkimmer.h +++ b/sbncode/Calibration/TrackCaloSkimmer.h @@ -67,7 +67,7 @@ #include "ITCSSelectionTool.h" // Useful functions -#include "sbncode/CAFMaker/RecoUtils/RecoUtils.h" +#include "sbncode/CAFMaker/RecoUtils/RecoUtils.h" // sbn::ReadoutIDE, sbn::PrepTrueHits()... namespace sbn { class TrackCaloSkimmer;