Skip to content

Feature/sungbino caf calo updates 2025 a #559

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
kjplows merged 15 commits into from
Sep 5, 2025
Merged

Feature/sungbino caf calo updates 2025 a #559

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
ba72c99
Cleanup code for getting IDEs associated with SimChannels
gputnam Aug 7, 2025
90cdf65
Draft of code to get costh_drift (cos_phi, e.g.) and the efield into …
gputnam Aug 7, 2025
98b624d
Incorporate updates to CAF calo-truth matching into calibration ntupl…
gputnam Aug 7, 2025
aa5484c
Lookup TrackHitMeta object correctly.
gputnam Aug 26, 2025
5b70b40
Merge branch 'feature/gputnam-caf-calo-updates' into feature/sungbino…
sungbinoh Aug 27, 2025
455cb8d
adding phi
sungbinoh Aug 27, 2025
b15f1f8
Merge branch 'release/SBN2025A' into feature/sungbino_caf-calo-update…
kjplows Aug 28, 2025
e70f014
following Gianluca's comments part1
sungbinoh Aug 29, 2025
2feb51f
following Gianluca's suggestion part 2
sungbinoh Aug 29, 2025
c6a3c29
ParticleIDE -> ReadoutIDE
sungbinoh Aug 29, 2025
035418f
minor bug fixes
sungbinoh Aug 29, 2025
a6b090a
moving ReadoutIDE to TrackCaloSkimmer.h and in sbn namespace. refer t…
sungbinoh Sep 2, 2025
1e0a36f
updating phi to rad from degree, and modifying some lines in CMakeLis…
sungbinoh Sep 3, 2025
46e144f
Refactor CAF helper functions utilized in CaloSkimmer into RecoUtils.
gputnam Sep 3, 2025
6df0c30
refer to *sce in CAFMaker_module.cc from FillReco, indentation fixes
sungbinoh Sep 3, 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
9 changes: 5 additions & 4 deletions sbncode/CAFMaker/CAFMaker_module.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -1386,7 +1387,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
}

// Prep truth-to-reco-matching info
std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> id_to_ide_map;
std::map<int, std::vector<caf::ParticleIDE>> id_to_ide_map;
std::map<int, std::vector<art::Ptr<recob::Hit>>> id_to_truehit_map;
std::map<int, caf::HitsEnergy> id_to_hit_energy_map;

Expand Down Expand Up @@ -1984,8 +1985,8 @@ void CAFMaker::produce(art::Event& evt) noexcept {
FindManyPStrict<recob::Vertex>(fmPFPart, evt,
fParams.PFParticleLabel() + slice_tag_suff);

art::FindManyP<recob::Hit> fmTrackHit =
FindManyPStrict<recob::Hit>(slcTracks, evt,
art::FindManyP<recob::Hit, recob::TrackHitMeta> fmTrackHit =
FindManyPDStrict<recob::Hit, recob::TrackHitMeta>(slcTracks, evt,
fParams.RecoTrackLabel() + slice_tag_suff);

art::FindManyP<recob::Hit> fmShowerHit =
Expand Down Expand Up @@ -2239,7 +2240,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);
Expand Down
46 changes: 44 additions & 2 deletions sbncode/CAFMaker/FillReco.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"

@PetrilloAtWork PetrilloAtWork Sep 3, 2025

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

These will not be needed after the change below.

Suggested change
#include "larevt/SpaceChargeServices/SpaceChargeService.h"
#include "larcore/CoreUtils/ServiceUtil.h"

@sungbinoh sungbinoh Sep 3, 2025

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Done! And moved #include "larevt/SpaceCharge/SpaceCharge.h" into FillReco.h from FillReco.cxx.


namespace caf
{
const float ng_filter_cut = 0.5;
Expand Down Expand Up @@ -794,8 +798,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<spacecharge::SpaceChargeService>();

@PetrilloAtWork PetrilloAtWork Sep 3, 2025

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

I am sorry, this is frustrating me and I am sure it will frustrate you too... somehow a comment on this function was not received by GitHub.
So the gist of it is that invoking a service from art at every point is not good, and this should happen only once in CAFMaker::produce().
So my request is to pass the space charge service provider as an argument from upstream:

Suggested change
double GetEfield(const detinfo::DetectorPropertiesData& dprop, const geo::Point_t loc) {
auto const* sce = lar::providerFrom<spacecharge::SpaceChargeService>();
double GetEfield(const detinfo::DetectorPropertiesData& dprop, spacecharge::SpaceCharge const& sce, const geo::Point_t& loc) {

Notes:

  • I am suggesting to pass a reference instead of a pointer to be consistent with dprop service data — that will require changing access below from sce-> to sce..
  • I am also suggesting to pass loc by reference rather than by value.
  • In CAFMaker::produce() there is already the variable sce defined (as pointer), which should be passed through FillTrackPlaneCalo() to here.

@sungbinoh sungbinoh Sep 3, 2025

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Done! *sce is passed from CAFMaker::produce() into FillTrackCalo() as reference, and into FillTrackPlaneCalo() from FillTrackCalo().


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};
Comment thread
sungbinoh marked this conversation as resolved.
Outdated

@PetrilloAtWork PetrilloAtWork Aug 28, 2025

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

It is not clear from SpaceCharge::GetEfieldOffset() whether it returns a "relative" offset (that is, assuming that the field vector is pointing on positive axis) or an absolute one (that is, in detector reference). In the latter case, the expression to the full field would need to include the absolute direction of the nominal field at the location loc.
Please check this, and document the assumption adding a comment above this line.

@sungbinoh sungbinoh Aug 29, 2025

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Thank you, added more words in the comments!

// Convert to Absolute E Field from relative
EFieldOffsets = EField * EFieldOffsets;
Comment thread
sungbinoh marked this conversation as resolved.
Outdated
// 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<art::Ptr<recob::Hit>> &hits,
const std::vector<const recob::TrackHitMeta*>& thms,
bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend,
const detinfo::DetectorPropertiesData &dprop,
caf::SRTrackCalo &srcalo) {
Expand Down Expand Up @@ -832,7 +856,8 @@ namespace caf

// lookup the wire -- the Calorimery object makes this
// __way__ harder than it should be
for (const art::Ptr<recob::Hit> &h: hits) {
for (unsigned i_hit = 0; i_hit < hits.size(); i_hit++) {
const art::Ptr<recob::Hit> &h = hits[i_hit];
if (h.key() == tps[i]) {
p.wire = h->WireID().Wire;
p.tpc = h->WireID().TPC;
Expand All @@ -844,6 +869,21 @@ 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(i_hit)->Index();
unsigned int int_max_as_unsigned_int{std::numeric_limits<int>::max()};
if (traj_point_index != int_max_as_unsigned_int && // invalid
track.HasValidPoint(traj_point_index)) {
Comment thread
sungbinoh marked this conversation as resolved.
Outdated
float costh_drift = track.DirectionAtPoint(traj_point_index).X();
float phi = acos(abs(costh_drift)) * 180. / M_PI;
float efield = GetEfield(dprop, track.LocationAtPoint(traj_point_index));
p.efield = efield;
p.phi = phi;
Comment thread
sungbinoh marked this conversation as resolved.
Outdated
}
}
}

Expand Down Expand Up @@ -897,7 +937,9 @@ namespace caf
}

void FillTrackCalo(const std::vector<art::Ptr<anab::Calorimetry>> &calos,
const recob::Track& track,
const std::vector<art::Ptr<recob::Hit>> &hits,
const std::vector<const recob::TrackHitMeta*>& thms,
bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend,
const detinfo::DetectorPropertiesData &dprop,
caf::SRTrack& srtrack,
Expand All @@ -912,7 +954,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]);
}
}

Expand Down
5 changes: 5 additions & 0 deletions sbncode/CAFMaker/FillReco.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -204,7 +205,9 @@ namespace caf
bool allowEmpty = false);

void FillTrackPlaneCalo(const anab::Calorimetry &calo,
const recob::Track& track,
const std::vector<art::Ptr<recob::Hit>> &hits,
const std::vector<const recob::TrackHitMeta*>& thms,
bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend,
const detinfo::DetectorPropertiesData &dprop,
caf::SRTrackCalo &srcalo);
Expand All @@ -222,7 +225,9 @@ namespace caf
bool allowEmpty = false);

void FillTrackCalo(const std::vector<art::Ptr<anab::Calorimetry>> &calos,
const recob::Track& track,
const std::vector<art::Ptr<recob::Hit>> &hits,
const std::vector<const recob::TrackHitMeta*>& thms,
bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend,
const detinfo::DetectorPropertiesData &dprop,
caf::SRTrack& srtrack,
Expand Down
26 changes: 13 additions & 13 deletions sbncode/CAFMaker/FillTrue.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -136,7 +136,7 @@ namespace caf {
}//FillTrackTruth

// Assumes truth matching and calo-points are filled
void FillTrackCaloTruth(const std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> &id_to_ide_map,
void FillTrackCaloTruth(const std::map<int, std::vector<caf::ParticleIDE>> &id_to_ide_map,
const std::vector<simb::MCParticle> &mc_particles,
const geo::GeometryCore& geometry,
const geo::WireReadoutGeom& wireReadout,
Expand All @@ -163,10 +163,10 @@ namespace caf {

// Load the hits
// match on the channel, which is unique
const std::vector<std::pair<geo::WireID, const sim::IDE*>> &match_ides = id_to_ide_map.at(srtrack.truth.p.G4ID);
const std::vector<caf::ParticleIDE> &match_ides = id_to_ide_map.at(srtrack.truth.p.G4ID);
std::map<unsigned, std::vector<const sim::IDE *>> 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
Expand Down Expand Up @@ -634,15 +634,15 @@ namespace caf {
void FillTrueG4Particle(const simb::MCParticle &particle,
const std::vector<geo::BoxBoundedGeo> &active_volumes,
const std::vector<std::vector<geo::BoxBoundedGeo>> &tpc_volumes,
const std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE *>>> &id_to_ide_map,
const std::map<int, std::vector<caf::ParticleIDE>> &id_to_ide_map,
const std::map<int, std::vector<art::Ptr<recob::Hit>>> &id_to_truehit_map,
const cheat::BackTrackerService &backtracker,
const cheat::ParticleInventoryService &inventory_service,
const std::vector<art::Ptr<simb::MCTruth>> &neutrinos,
caf::SRTrueParticle &srparticle) {

std::vector<std::pair<geo::WireID, const sim::IDE *>> empty;
const std::vector<std::pair<geo::WireID, const sim::IDE *>> &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty;
std::vector<caf::ParticleIDE> empty;
const std::vector<caf::ParticleIDE> &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty;

std::vector<art::Ptr<recob::Hit>> emptyHits;
const std::vector<art::Ptr<recob::Hit>> &particle_hits = id_to_truehit_map.count(particle.TrackId()) ? id_to_truehit_map.at(particle.TrackId()) : emptyHits;
Expand All @@ -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*/;
Expand Down Expand Up @@ -882,8 +882,8 @@ namespace caf {
return ret;
}

std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> PrepSimChannels(const std::vector<art::Ptr<sim::SimChannel>> &simchannels, const geo::WireReadoutGeom &wireReadout) {
std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> ret;
std::map<int, std::vector<caf::ParticleIDE>> PrepSimChannels(const std::vector<art::Ptr<sim::SimChannel>> &simchannels, const geo::WireReadoutGeom &wireReadout) {
std::map<int, std::vector<caf::ParticleIDE>> ret;

for (const art::Ptr<sim::SimChannel> sc : simchannels) {
// Lookup the wire of this channel
Expand All @@ -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});
}
}
}
Expand Down
13 changes: 10 additions & 3 deletions sbncode/CAFMaker/FillTrue.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,13 @@ namespace caf
float totE;
};

struct ParticleIDE {

@PetrilloAtWork PetrilloAtWork Aug 28, 2025

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

This class has not much to to with a particle. It has more to do with the TPC readout. Please consider a better name ("ReadoutIDE"? I am bad in finding names).

@sungbinoh sungbinoh Aug 29, 2025

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Agree! Updated the ParticleIDE to ReadoutIDE, I think that it is a good naming!

geo::WireID wire;
unsigned short tick;
const sim::IDE *ide;
Comment thread
sungbinoh marked this conversation as resolved.
Outdated
Comment thread
sungbinoh marked this conversation as resolved.
Outdated
};


// Helpers
caf::Wall_t GetWallCross( const geo::BoxBoundedGeo &volume,
const TVector3 p0,
Expand Down Expand Up @@ -73,7 +80,7 @@ namespace caf
void FillTrueG4Particle(const simb::MCParticle &particle,
const std::vector<geo::BoxBoundedGeo> &active_volumes,
const std::vector<std::vector<geo::BoxBoundedGeo>> &tpc_volumes,
const std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE *>>> &id_to_ide_map,
const std::map<int, std::vector<caf::ParticleIDE>> &id_to_ide_map,
const std::map<int, std::vector<art::Ptr<recob::Hit>>> &id_to_truehit_map,
const cheat::BackTrackerService &backtracker,
const cheat::ParticleInventoryService &inventory_service,
Expand Down Expand Up @@ -103,7 +110,7 @@ namespace caf
caf::SRTrack& srtrack,
bool allowEmpty = false);

void FillTrackCaloTruth(const std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> &id_to_ide_map,
void FillTrackCaloTruth(const std::map<int, std::vector<ParticleIDE>> &id_to_ide_map,
const std::vector<simb::MCParticle> &mc_particles,
const geo::GeometryCore & geometry,
const geo::WireReadoutGeom& wireReadout,
Expand Down Expand Up @@ -133,7 +140,7 @@ namespace caf
CLHEP::HepRandomEngine &rand,
std::vector<caf::SRFakeReco> &srfakereco);

std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> PrepSimChannels(const std::vector<art::Ptr<sim::SimChannel>> &simchannels, const geo::WireReadoutGeom &wireReadout);
std::map<int, std::vector<ParticleIDE>> PrepSimChannels(const std::vector<art::Ptr<sim::SimChannel>> &simchannels, const geo::WireReadoutGeom &wireReadout);
std::map<int, std::vector<art::Ptr<recob::Hit>>> PrepTrueHits(const std::vector<art::Ptr<recob::Hit>> &allHits,
const detinfo::DetectorClocksData &clockData, const cheat::BackTrackerService &backtracker);
std::map<int, caf::HitsEnergy> SetupIDHitEnergyMap(const std::vector<art::Ptr<recob::Hit>> &allHits, const detinfo::DetectorClocksData &clockData,
Expand Down
3 changes: 2 additions & 1 deletion sbncode/Calibration/TrackCaloSkimmer.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,7 @@
#include "sbnobj/Common/CRT/CRTHitT0TaggingInfo.hh"
#include "sbnobj/Common/CRT/CRTHitT0TaggingTruthInfo.hh"

#include "sbncode/CAFMaker/FillTrue.h"

@PetrilloAtWork PetrilloAtWork Aug 28, 2025

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

This functionality should not depend on CAF making (see the general review comment).

#include "ITCSSelectionTool.h"

namespace sbn {
Expand Down Expand Up @@ -161,7 +162,7 @@ class sbn::TrackCaloSkimmer : public art::EDAnalyzer {
const std::vector<art::Ptr<simb::MCParticle>> &mcparticles,
const std::vector<geo::BoxBoundedGeo> &active_volumes,
const std::vector<std::vector<geo::BoxBoundedGeo>> &tpc_volumes,
const std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> id_to_ide_map,
const std::map<int, std::vector<caf::ParticleIDE>> id_to_ide_map,

@PetrilloAtWork PetrilloAtWork Aug 28, 2025

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

Object parameters should be passed by reference unless there are explicit reasons to do otherwise ([CF-112]:

Suggested change
const std::map<int, std::vector<caf::ParticleIDE>> id_to_ide_map,
const std::map<int, std::vector<caf::ParticleIDE>>& id_to_ide_map,

const std::map<int, std::vector<art::Ptr<recob::Hit>>> id_to_truehit_map,
const detinfo::DetectorPropertiesData &dprop,
const geo::GeometryCore *geo,
Expand Down
33 changes: 17 additions & 16 deletions sbncode/Calibration/TrackCaloSkimmer_module.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> id_to_ide_map;
std::map<int, std::vector<caf::ParticleIDE>> id_to_ide_map;
std::map<int, std::vector<art::Ptr<recob::Hit>>> id_to_truehit_map;
const cheat::BackTrackerService *bt = NULL;

Expand Down Expand Up @@ -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<geo::BoxBoundedGeo> &active_volumes,
const std::vector<std::vector<geo::BoxBoundedGeo>> &tpc_volumes,
const std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE *>>> &id_to_ide_map,
const std::map<int, std::vector<caf::ParticleIDE>> &id_to_ide_map,
const std::map<int, std::vector<art::Ptr<recob::Hit>>> &id_to_truehit_map,
const detinfo::DetectorPropertiesData &dprop,
const geo::GeometryCore *geo,
const geo::WireReadoutGeom *wireReadout) {

std::vector<std::pair<geo::WireID, const sim::IDE *>> empty;
const std::vector<std::pair<geo::WireID, const sim::IDE *>> &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty;
std::vector<caf::ParticleIDE> empty;
const std::vector<caf::ParticleIDE> &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty;

@PetrilloAtWork PetrilloAtWork Aug 28, 2025

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

The current implementation performs the same query twice, and that would be better avoided. A way to avoid it, wrapped in a lambda function:

Suggested change
std::vector<caf::ParticleIDE> empty;
const std::vector<caf::ParticleIDE> &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty;
auto trackReadoutIDEs = [&id_to_ide_map](int trackID) -> std::vector<caf::ParticleIDE> const&
{
static std::vector<caf::ParticleIDE> const empty;
auto const itIDEs = id_to_ide_map.find(trackID);
return (itIDEs == id_to_ide_map.end())? empty: itIDEs->second;
};
std::vector<caf::ParticleIDE> const& particle_ides = trackReadoutIDEs(particle.TrackId());

This is a suggestion.

@gputnam gputnam Aug 29, 2025

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

I think it's ok as is -- this would be a small optimization at the cost of readability.


std::vector<art::Ptr<recob::Hit>> emptyHits;
const std::vector<art::Ptr<recob::Hit>> &particle_hits = id_to_truehit_map.count(particle.TrackId()) ? id_to_truehit_map.at(particle.TrackId()) : emptyHits;
Expand All @@ -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*/;
Expand Down Expand Up @@ -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<unsigned, sbn::TrueHit> 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;
Expand All @@ -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;

@PetrilloAtWork PetrilloAtWork Aug 28, 2025

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

This way of processing leaves me puzzled. I think it originates in 12-year old code that was written that way while interfacing with GEANT4 steps (so, at the same time when the true hits were being created) so that the information would be always valid. But here, in a single loop, it's just a waste times the 6 lines it's used in (7 now).
Not saying it's your duty to fix it, but it would be nice to.
Basically the alternative is either to use a sum in the look, and then divide once at the end, with a dedicated loop, or to use averaging objects (like lar::util::StatCollector and geo::vect::MiddlePointAccumulator).

@gputnam gputnam Aug 29, 2025

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

This is written purposefully. The number of electrons can get quite large, so you get numerical stability issues if you do a weighted sum the naive way. (sum p_i * w_i / sum w_i). I'm not sure if you mean something else, but in any case I think this is ok as is.

@PetrilloAtWork PetrilloAtWork Sep 3, 2025 •
edited
Loading

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

I do mean that.
So, you take truehits[c].time*old_elec + tick*ide->numElectrons, which I will call time_sum, you divide by new_elec, and then at the next loop you multiply again for new_elec (which in the meanwhile has become old_elec) to get time_sum again, and add a number.
How is that ensuring stability compared to storing old_time_sum all the way?

@PetrilloAtWork PetrilloAtWork Sep 3, 2025

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

(I actually do have my answer, but it assumes the internal precision of the FPU stack/registers is larger than the output CPU registry, plus some compiler optimisation)


// 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);
Expand All @@ -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);
Expand All @@ -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();
}
Expand Down Expand Up @@ -1010,7 +1011,7 @@ void sbn::TrackCaloSkimmer::FillTrackTruth(const detinfo::DetectorClocksData &cl
const std::vector<art::Ptr<simb::MCParticle>> &mcparticles,
const std::vector<geo::BoxBoundedGeo> &active_volumes,
const std::vector<std::vector<geo::BoxBoundedGeo>> &tpc_volumes,
const std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> id_to_ide_map,
const std::map<int, std::vector<caf::ParticleIDE>> id_to_ide_map,
const std::map<int, std::vector<art::Ptr<recob::Hit>>> id_to_truehit_map,
const detinfo::DetectorPropertiesData &dprop,
const geo::GeometryCore *geo,
Expand Down