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 12 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<sbn::ReadoutIDE>> 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
1 change: 1 addition & 0 deletions sbncode/CAFMaker/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ art_make_library( LIBRARY_NAME sbncode_CAFMaker
larcorealg::Geometry
larcore::Geometry_Geometry_service
larcorealg::GeoAlgo
lardata::DetectorInfoServices_DetectorPropertiesServiceStandard_service

@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.

Missed this... This library should not be necessary — if it is, there is a problem.
Try with:

Suggested change
lardata::DetectorInfoServices_DetectorPropertiesServiceStandard_service
lardataalg::headers

(also, make sure it's appropriately indented) for the library, and add

               lardata::DetectorPropertiesService

below (e.g. around line 69) for the module.
Rationale: DetectorPropertiesServiceStandard is an implementation of a service interface. Which implementation to use is decided run-time by reading the FHiCL services configuration, so it does not make sense to link a specific implementation at build time. What we really need is the "library" pulling in the interface. I hope my suggestion is the LArSoft's right name for the one we need here.

@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 have checked FillReco.h... we do have a problem: lardataalg/DetectorInfo/DetectorPropertiesStandard.h is included directly. It should not be. Ever.
Please replace that line with the appropriate header, which I think is:

#include "lardataalg/DetectorInfo/DetectorPropertiesData.h"

(bonus points if you try to remove also the inclusion of ParticleInventoryService.h on the previous line, since it looks like it's not used — and if it is needed, leave it there)
With my thanks for fixing this pre-existing bug.

larsim::MCCheater_BackTrackerService_service
nusimdata::SimulationBase
larsim::MCCheater_ParticleInventoryService_service
Expand Down
44 changes: 42 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 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|

@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.

Just to make double sure: does the offset as returned by sce always assume the nominal electric field to have positive x direction?
That assumption is an essential ingredient.

@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.

It depends on how the SCE map is saved. For SBN SCE studies, it is assumed the nominal electric field has positive x direction.

@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.

Yes, it depends on that... unfortunately I don't see that assumption documented.
However, I see many places where that assumption is made, and if you haven't noticed anything strange in testing, that will have to do.

EFieldOffsets += geo::Vector_t{1, 0, 0};
// Convert to Absolute E Field from relative
EFieldOffsets *= EField;
// 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,19 @@ 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();
if (track.HasValidPoint(traj_point_index)) {

@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.

Please fix this indentation.

Suggested change
if (track.HasValidPoint(traj_point_index)) {
if (track.HasValidPoint(traj_point_index)) {

@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.

Fixed!

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;
}
}
}

Expand Down Expand Up @@ -897,7 +935,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 +952,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<sbn::ReadoutIDE>> &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<sbn::ReadoutIDE> &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<sbn::ReadoutIDE>> &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<sbn::ReadoutIDE> empty;
const std::vector<sbn::ReadoutIDE> &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<sbn::ReadoutIDE>> PrepSimChannels(const std::vector<art::Ptr<sim::SimChannel>> &simchannels, const geo::WireReadoutGeom &wireReadout) {
std::map<int, std::vector<sbn::ReadoutIDE>> 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
8 changes: 5 additions & 3 deletions sbncode/CAFMaker/FillTrue.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,8 @@
#include "sbnanaobj/StandardRecord/StandardRecord.h"
#include "sbnanaobj/StandardRecord/SRMeVPrtl.h"

#include "sbncode/Calibration/TrackCaloSkimmer.h"

namespace caf
{
struct HitsEnergy {
Expand Down Expand Up @@ -73,7 +75,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<sbn::ReadoutIDE>> &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 +105,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<sbn::ReadoutIDE>> &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 +135,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<sbn::ReadoutIDE>> 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
9 changes: 7 additions & 2 deletions sbncode/Calibration/TrackCaloSkimmer.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,12 @@
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 {
Expand Down Expand Up @@ -112,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;
Expand Down Expand Up @@ -161,7 +166,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<sbn::ReadoutIDE>> 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<sbn::ReadoutIDE>> 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<sbn::ReadoutIDE>> &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<sbn::ReadoutIDE> empty;
const std::vector<sbn::ReadoutIDE> &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 @@ -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<sbn::ReadoutIDE>> 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