Feature/sungbino caf calo updates 2025 a #559

sungbinoh · 2025-08-27T15:20:59Z

Description

Adding phi and |E| for each hit for caf-level calorimetry
Related PR in sbnanaobj: PR151

Have you added a label? (bug/enhancement/physics etc.)
Have you assigned at least 1 reviewer?
Is this PR related to an open issue / project?
Does this PR affect CAF data format? If so, please assign a CAF maintainer as additional reviewer.
Does this PR require merging another PR in a different repository (such as sbnanobj/sbnobj etc.)? If so, please link it in the description.
Are you submitting this PR on behalf of someone else who made the code changes? If so, please mention them in the description.

MC BNB + cosmic flat.caf

phi

efield

Data BNB + light flat.caf

phi

efield

Data Offbeam + light flat.caf

phi

efield

…the CAF SRCalo object per-hit.

…e'r. Use actualy true deposition time to define true calo point time.

…_caf-calo-updates_tested

…s_2025A

PetrilloAtWork

Points to discuss:

I oppose having TrackCaloSkimmer depend on CAFMaker. CAF making has no role in TrackCaloSkimmer. Having a shared type definition (caf::ParticleIDE) is good, but it should live either in Calibration library, or in its own place (and be in sbn namespace rather than caf).
The definition of the electric field correction may imply a bug in the code. The result of an investigation on the matter should be written down as comment in the code.
I ask to reconsider the unit of the angle phi. Having it in degrees means that every computation will need to convert it to radians. Personally I would consider to store its cosine instead, but if an angle is required, I believe it should be stored as radians.

I have left some additional points here and there, mostly C++ style and recommendations.

PetrilloAtWork · 2025-08-28T19:23:55Z

sbncode/Calibration/TrackCaloSkimmer_module.cc

    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;


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

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.

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?

(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)

sbncode/CAFMaker/FillReco.cxx

PetrilloAtWork · 2025-08-28T20:30:44Z

sbncode/CAFMaker/FillReco.cxx

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


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.

Thank you, added more words in the comments!

sbncode/CAFMaker/FillTrue.h

PetrilloAtWork · 2025-08-28T20:59:45Z

sbncode/CAFMaker/FillTrue.h

    float totE;
  };

+  struct ParticleIDE {


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

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

PetrilloAtWork · 2025-08-28T21:04:56Z

sbncode/Calibration/TrackCaloSkimmer.h

    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,


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,

PetrilloAtWork · 2025-08-28T21:23:11Z

sbncode/Calibration/TrackCaloSkimmer_module.cc

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


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.

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

sbncode/CAFMaker/FillReco.cxx

…o this from CAFMaker's side

sungbinoh · 2025-09-02T21:55:35Z

Points to discuss:

I oppose having TrackCaloSkimmer depend on CAFMaker. CAF making has no role in TrackCaloSkimmer. Having a shared type definition (caf::ParticleIDE) is good, but it should live either in Calibration library, or in its own place (and be in sbn namespace rather than caf).

The definition of the electric field correction may imply a bug in the code. The result of an investigation on the matter should be written down as comment in the code.

I ask to reconsider the unit of the angle phi. Having it in degrees means that every computation will need to convert it to radians. Personally I would consider to store its cosine instead, but if an angle is required, I believe it should be stored as radians.

I have left some additional points here and there, mostly C++ style and recommendations.

Hi @PetrilloAtWork , Thank you! I have made several updates following the suggestions.
For points to discuss,

I have moved the ReadoutIDE's definition into the TrackCaloSkimmer.h and in the sbn namespace. It is referred in CAFMAker's side.
Have added additional texts to resolve the confusion
I agree with your point but keeping it as degree to have the exactly same value provided by the GnocchiCalorimetry_module.cc. Converting it to radian should be simple.

Please also have a look to our responses to your additional points.
Thank you!!

PetrilloAtWork

Thank you for the answers. There still a clarification that I would like on the previous comments, but more important is that I noticed something I missed in the first review about the service management.
In short:

link list should only include service interfaces, not implementation (like the Standard ones)
loading of the service should be limited to once per produce() call.

Inline comments explain more.

PetrilloAtWork · 2025-09-03T01:44:12Z

sbncode/CAFMaker/CMakeLists.txt

                  larcorealg::Geometry
                  larcore::Geometry_Geometry_service
                  larcorealg::GeoAlgo
+		  lardata::DetectorInfoServices_DetectorPropertiesServiceStandard_service


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.

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.

PetrilloAtWork · 2025-09-03T02:21:06Z

sbncode/CAFMaker/FillReco.cxx

+    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|


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.

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

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.

PetrilloAtWork · 2025-09-03T02:54:03Z

Points to discuss:
[...]
3. I ask to reconsider the unit of the angle phi. Having it in degrees means that every computation will need to convert it to radians. Personally I would consider to store its cosine instead, but if an angle is required, I believe it should be stored as radians.
[...]

I agree with your point but keeping it as degree to have the exactly same value provided by the GnocchiCalorimetry_module.cc. Converting it to radian should be simple.

It is not about how it was written in the original code (in the same way in which it does not matter whether an energy was stored in MeV or GeV or erg), but rather how people expect it, and how it is expected to be used. CAF promise is that angles would be in radians; there needs to be a real need for that promise to be broken. I am not even sure there is any exception now (maybe some times bound to timestamps stored in nanoseconds).

…ts.txt and FillReco.h headers

sungbinoh · 2025-09-03T15:36:13Z

Points to discuss:
[...]
3. I ask to reconsider the unit of the angle phi. Having it in degrees means that every computation will need to convert it to radians. Personally I would consider to store its cosine instead, but if an angle is required, I believe it should be stored as radians.
[...]

I agree with your point but keeping it as degree to have the exactly same value provided by the GnocchiCalorimetry_module.cc. Converting it to radian should be simple.

It is not about how it was written in the original code (in the same way in which it does not matter whether an energy was stored in MeV or GeV or erg), but rather how people expect it, and how it is expected to be used. CAF promise is that angles would be in radians; there needs to be a real need for that promise to be broken. I am not even sure there is any exception now (maybe some times bound to timestamps stored in nanoseconds).

That makes sense. (sorry, I've not known that there has been CAF promise). Updated phi to rad from degree.

gputnam · 2025-09-03T16:37:41Z

@PetrilloAtWork I believe all of your comments have now been addressed. Please take another look.

PetrilloAtWork

There are some indentation fixes pending.
Worse, GitHub ate a request from yesterday's review.
The rest should be ok.

PetrilloAtWork · 2025-09-03T19:45:30Z

sbncode/CAFMaker/FillReco.cxx

+            //
+            // 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)) {


Please fix this indentation.

Suggested change

if (track.HasValidPoint(traj_point_index)) {

if (track.HasValidPoint(traj_point_index)) {

PetrilloAtWork · 2025-09-03T19:49:55Z

sbncode/Calibration/CMakeLists.txt

                            sbncode_CAFMaker
                            sbnobj::Common_Calibration_dict
                            larevt::SpaceCharge
+			    caf_RecoUtils


Please fix the indentation:

Suggested change

caf_RecoUtils

caf_RecoUtils

PetrilloAtWork · 2025-09-03T19:51:14Z

sbncode/Calibration/TrackCaloSkimmer.h

+// Useful functions
+#include "sbncode/CAFMaker/RecoUtils/RecoUtils.h"


I suggest a comment detailing why they are useful 😉

Suggested change

// Useful functions

#include "sbncode/CAFMaker/RecoUtils/RecoUtils.h"

#include "sbncode/CAFMaker/RecoUtils/RecoUtils.h" // sbn::ReadoutIDE, sbn::PrepTrueHits()...

PetrilloAtWork · 2025-09-03T20:04:48Z

sbncode/CAFMaker/FillReco.cxx

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


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.

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

PetrilloAtWork · 2025-09-03T20:05:25Z

sbncode/CAFMaker/FillReco.cxx

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


These will not be needed after the change below.

Suggested change

#include "larevt/SpaceChargeServices/SpaceChargeService.h"

#include "larcore/CoreUtils/ServiceUtil.h"

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

FNALbuild · 2025-09-04T18:11:12Z

✔️ CI build for LArSoft Succeeded on slf7 for c14:prof -- details available through the CI dashboard

FNALbuild · 2025-09-04T18:15:25Z

✔️ CI build for LArSoft Succeeded on slf7 for e26:prof -- details available through the CI dashboard

FNALbuild · 2025-09-04T18:20:40Z

❌ CI build for SBND Failed at phase build SBND on slf7 for c14:prof -- details available through the CI dashboard

🚨 For more details about the failed phase, check the build SBND phase logs

parent CI build details are available through the CI dashboard

FNALbuild · 2025-09-04T18:22:52Z

❌ CI build for ICARUS Failed at phase build ICARUS on slf7 for c14:prof -- details available through the CI dashboard

🚨 For more details about the failed phase, check the build ICARUS phase logs

parent CI build details are available through the CI dashboard

FNALbuild · 2025-09-04T18:37:05Z

⚠️ CI build for SBND Warning at phase ci_tests SBND on slf7 for e26:prof - ignored warnings for build -- details available through the CI dashboard

🚨 For more details about the warning phase, check the ci_tests SBND phase logs

parent CI build details are available through the CI dashboard

FNALbuild · 2025-09-04T18:39:40Z

❌ CI build for ICARUS Failed at phase ci_tests ICARUS on slf7 for e26:prof - ignored warnings for build -- details available through the CI dashboard

🚨 For more details about the failed phase, check the ci_tests ICARUS phase logs

parent CI build details are available through the CI dashboard

gputnam and others added 6 commits August 6, 2025 19:49

Cleanup code for getting IDEs associated with SimChannels

ba72c99

Draft of code to get costh_drift (cos_phi, e.g.) and the efield into …

90cdf65

…the CAF SRCalo object per-hit.

Incorporate updates to CAF calo-truth matching into calibration ntupl…

98b624d

…e'r. Use actualy true deposition time to define true calo point time.

Lookup TrackHitMeta object correctly.

aa5484c

Merge branch 'feature/gputnam-caf-calo-updates' into feature/sungbino…

5b70b40

…_caf-calo-updates_tested

adding phi

455cb8d

sungbinoh requested a review from PetrilloAtWork August 27, 2025 15:21

sungbinoh self-assigned this Aug 27, 2025

sungbinoh added the enhancement New feature or request label Aug 27, 2025

sungbinoh requested a review from jzennamo August 27, 2025 15:22

sungbinoh mentioned this pull request Aug 27, 2025

adding efield and phi to SRCaloPoint for 2025A SBNSoftware/sbnanaobj#151

Merged

sungbinoh added this to SBND March 2025 production Aug 27, 2025

sungbinoh moved this to Todo in SBND March 2025 production Aug 27, 2025

kjplows added this to SBN software development Aug 27, 2025

kjplows moved this to Testing in SBN software development Aug 27, 2025

Merge branch 'release/SBN2025A' into feature/sungbino_caf-calo-update…

b15f1f8

…s_2025A

PetrilloAtWork requested changes Aug 28, 2025

View reviewed changes

sungbinoh added 5 commits August 29, 2025 12:07

following Gianluca's comments part1

e70f014

following Gianluca's suggestion part 2

2feb51f

ParticleIDE -> ReadoutIDE

c6a3c29

minor bug fixes

035418f

moving ReadoutIDE to TrackCaloSkimmer.h and in sbn namespace. refer t…

a6b090a

…o this from CAFMaker's side

PetrilloAtWork requested changes Sep 3, 2025

View reviewed changes

updating phi to rad from degree, and modifying some lines in CMakeLis…

1e0a36f

…ts.txt and FillReco.h headers

Refactor CAF helper functions utilized in CaloSkimmer into RecoUtils.

46e144f

PetrilloAtWork requested changes Sep 3, 2025

View reviewed changes

SBNSoftware deleted a comment from FNALbuild Sep 4, 2025

kjplows merged commit 8c5781a into SBNSoftware:release/SBN2025A Sep 5, 2025
3 of 6 checks passed

github-project-automation bot moved this from Testing to Done in SBN software development Sep 5, 2025

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

		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<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());

	lardata::DetectorInfoServices_DetectorPropertiesServiceStandard_service
	lardataalg::headers

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

		// Useful functions
		#include "sbncode/CAFMaker/RecoUtils/RecoUtils.h"

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

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

Feature/sungbino caf calo updates 2025 a #559

Feature/sungbino caf calo updates 2025 a #559

Uh oh!

Conversation

sungbinoh commented Aug 27, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Description

MC BNB + cosmic flat.caf

Data BNB + light flat.caf

Data Offbeam + light flat.caf

Uh oh!

PetrilloAtWork left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

PetrilloAtWork Sep 3, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

sungbinoh commented Sep 2, 2025

Uh oh!

PetrilloAtWork left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

PetrilloAtWork commented Sep 3, 2025

Uh oh!

sungbinoh commented Sep 3, 2025

Uh oh!

gputnam commented Sep 3, 2025

Uh oh!

PetrilloAtWork left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

sungbinoh commented Aug 27, 2025 •

edited

Loading

PetrilloAtWork Sep 3, 2025 •

edited

Loading