Run 1B: Calo cluster-based time clustering

[michaelmackenzie]

Cluster tracker hits using the stopping target location and a calorimeter cluster. This uses a simple cone from the cluster to select hits consistent in time/space with the calo cluster, using the configured trajectory direction.

[FNALbuild]

Hi @michaelmackenzie,
You have proposed changes to files in these packages:

• CalPatRec

which require these tests: build.

@Mu2e/write have access to CI actions on Run1B.

⌛ The following tests have been triggered for bd8c81a: build (Build queue - API unavailable)

About FNALbuild. Code review on Mu2e/Offline.

[oksuzian]

I'll review this pull request by first loading the PR details and its changes.Let me get the actual diff/files changed in this PR.Here's my review of PR #1733. This PR adds a new CalLineTimePeakFinder module that clusters tracker hits using a calorimeter cluster seed and a straight-line cone from the stopping target. Below are the issues I found:

---

🔴 Bugs / Correctness Issues

1. Missing stopping_target_radius_ initialization in constructor

The constructor initializer list does not initialize stopping_target_radius_ from the config, even though it is declared as a member and used in findTimePeakInCluster. This means it will hold an indeterminate value.

  CalLineTimePeakFinder::CalLineTimePeakFinder(const art::EDProducer::Table<Config>& config) :
    art::EDProducer{config}
    , hit_tag_(config().hit_coll_tag())
    , calo_cluster_tag_(config().calo_cluster_coll_tag())
    , min_tc_hits_(config().min_tc_hits())
    , min_calo_cluster_energy_(config().min_calo_cluster_energy())
    , hit_time_sigma_thresh_(config().hit_time_sigma_thresh())
    , hit_xy_sigma_thresh_(config().hit_xy_sigma_thresh())
    , fit_dir_(config().fit_direction())
    , diag_level_(config().diag_level())

Fix: Add , stopping_target_radius_(config().stopping_target_radius()) to the initializer list (between hit_xy_sigma_thresh_ and fit_dir_).

2. Potential division by zero in finalizeTimeCluster

If weight is zero (which shouldn't normally happen given the isGoodTimeCluster check, but could if all nStrawHits() return 0), the code will divide by zero when computing avg_t0 /= weight and tc._pos = tc_pos / weight.

    avg_t0 /= weight;
    t0sq /= weight;
    tc._pos = tc_pos / weight;
    tc._t0._t0 = avg_t0;
    tc._t0._t0err = std::sqrt(t0sq  - avg_t0*avg_t0);

Suggestion: Add a guard if(weight > 0.) before the divisions.

3. Potential NaN from std::sqrt of negative argument in finalizeTimeCluster

Due to floating-point rounding, t0sq - avg_t0*avg_t0 could be slightly negative (when all hit times are nearly identical), producing NaN from std::sqrt.

Suggestion: Use std::sqrt(std::max(0., t0sq - avg_t0*avg_t0)) to clamp.

4. Spatial cone projection is incorrect when seed_dir.z() is small or negative unexpectedly

The spatial consistency check uses dz to project along the seed direction, but pos_at_hit is computed as cl_pos + dz * seed_dir, which does not project to the same z-plane as the hit when seed_dir is not purely along z. This is because scaling the full 3D direction vector by the z-distance does not land at the hit's z-coordinate — it only works if seed_dir.z() ≈ 1.

      const double dz = (seed_dir.z() > 0.) ? hit_pos.z() - cl_pos.z() : cl_pos.z() - hit_pos.z();
      const CLHEP::Hep3Vector pos_at_hit = cl_pos + dz * seed_dir;
      const double x_y_dist = (hit_pos - pos_at_hit).perp();

Fix: To project along the line to a specific z-plane, the parameter should be t = dz / seed_dir.z(), and then pos_at_hit = cl_pos + t * seed_dir. This correctly parameterizes the line by z. The current code multiplies seed_dir (which has x/y components too) by dz, which overshoots in z for off-axis trajectories.

---

🟡 Design / Robustness Concerns

5. hit_indices vector is declared but never used

A local std::vector<size_t> hit_indices is declared and reserved but never populated or used. Hits are instead pushed directly into tc._strawHitIdxs.

    std::vector<size_t> hit_indices;
    const size_t n_hits = combo_hit_col_->size();
    hit_indices.reserve(100);

Fix: Remove the unused hit_indices variable.

6. Time propagation model assumes speed of light

The time-of-flight calculation uses CLHEP::c_light (the speed of light), but tracker particles (conversion electrons ~105 MeV) travel at β ≈ 0.99995c. While very close to c, if precision matters, this could be made configurable or use a more accurate velocity.

      const double time_at_hit = cl_time + (hit_pos - cl_pos).dot(seed_dir) / CLHEP::c_light;

7. Hits are exclusively assigned to the first matching cluster

The hits_used_in_tcs_ set prevents a hit from being assigned to more than one time cluster. This greedy approach is order-dependent — clusters processed earlier get priority. If two nearby calo clusters legitimately share overlapping hits, the second cluster may be starved of hits and fail the min_tc_hits_ cut.

8. event_ member pointer is stored but only used for diagnostics (if at all)

The raw pointer event_ is set in produce() but is never referenced anywhere else in the code. This is dead state and could be removed.

9. Forward declarations shadow actual includes

At line 39-41, class Calorimeter, class Tracker, and class ModuleHistToolBase are forward-declared inside the mu2e namespace, but the full headers are already included above. The ModuleHistToolBase forward declaration is entirely unnecessary as it's never used.

10. Typo in comment

  // module constructor, parameter defaults are defiend in CalPatRec/fcl/prolog.fcl

"defiend" → "defined"

---

Summary

Severity	Issue	Location
🔴 Bug	stopping_target_radius_ not initialized from config	Constructor (L119-128)
🔴 Bug	Spatial cone projection math is incorrect for off-axis seeds	findTimePeakInCluster (L237-239)
🔴 Bug	Possible division by zero in finalizeTimeCluster	finalizeTimeCluster (L308)
🔴 Bug	std::sqrt of potentially negative value → NaN	finalizeTimeCluster (L321)
🟡 Cleanup	Unused hit_indices variable	findTimePeakInCluster (L216)
🟡 Cleanup	Unused event_ member pointer	Class member
🟡 Cleanup	Unnecessary forward declarations	L39-41
🟡 Design	Greedy hit assignment is order-dependent	findTimePeakInCluster / produce
🟡 Physics	TOF uses exact c instead of βc	findTimePeakInCluster (L228)
⚪ Typo	"defiend" in comment	L117

The most critical issue is #1 — stopping_target_radius_ is never initialized, meaning the cone-based spatial selection will use garbage values and produce incorrect results.

[FNALbuild]

☔ The build tests failed for bd8c81a.

Test	Result	Details
test with	✅	Command did not list any other PRs to include
merge	✅	Merged bd8c81a at b55e0e5
build (prof)	✅	Log file. Build time: 08 min 49 sec
ceSimReco	✅	Log file.
g4test_03MT	✅	Log file.
transportOnly	✅	Log file.
POT	✅	Log file.
g4study	✅	Log file.
cosmicSimReco	✅	Log file.
cosmicOffSpill	✅	Log file.
ceSteps	✅	Log file.
ceDigi	✅	Log file.
muDauSteps	✅	Log file.
ceMix	✅	Log file.
rootOverlaps	✅	Log file.
g4surfaceCheck	✅	Log file.
trigger	❌	Log file. Return Code 2.
check_cmake	❌	Log file. Return Code 1.
FIXME, TODO	✅	TODO (0) FIXME (0) in 0 files
clang-tidy	✅	0 errors 0 warnings
whitespace check	✅	no whitespace errors found

N.B. These results were obtained from a build of this Pull Request at bd8c81a after being merged into the base branch at b55e0e5.

For more information, please check the job page here.
Build artifacts are deleted after 5 days. If this is not desired, select Keep this build forever on the job page.

[michaelmackenzie]

@FNALbuild run build test

[FNALbuild]

⌛ The following tests have been triggered for 5fbd4cb: build (Build queue - API unavailable)

[FNALbuild]

☀️ The build tests passed at 5fbd4cb.

Test	Result	Details
test with	✅	Command did not list any other PRs to include
merge	✅	Merged 5fbd4cb at b55e0e5
build (prof)	✅	Log file. Build time: 04 min 09 sec
ceSimReco	✅	Log file.
g4test_03MT	✅	Log file.
transportOnly	✅	Log file.
POT	✅	Log file.
g4study	✅	Log file.
cosmicSimReco	✅	Log file.
cosmicOffSpill	✅	Log file.
ceSteps	✅	Log file.
ceDigi	✅	Log file.
muDauSteps	✅	Log file.
ceMix	✅	Log file.
rootOverlaps	✅	Log file.
g4surfaceCheck	✅	Log file.
trigger	❌	Log file. Return Code 2.
check_cmake	✅	Log file.
FIXME, TODO	✅	TODO (0) FIXME (0) in 0 files
clang-tidy	✅	0 errors 0 warnings
whitespace check	✅	no whitespace errors found

N.B. These results were obtained from a build of this Pull Request at 5fbd4cb after being merged into the base branch at b55e0e5.

For more information, please check the job page here.
Build artifacts are deleted after 5 days. If this is not desired, select Keep this build forever on the job page.

[oksuzian]

I'll review this pull request by first loading the PR details and its changes.Now let me get the actual file diffs to do a thorough independent review:Now I have the full diff. Let me analyze the actual code carefully. I note from the comments that you previously posted a review, and that the author has since pushed updates (the second build at 5fbd4cb passed check_cmake). Let me review the current version of the code from the latest diff.

---

Code Review: PR #1733 — Calo cluster-based time clustering

This PR adds a new CalLineTimePeakFinder art module (336 lines) that clusters tracker hits by seeding from calorimeter clusters and using a straight-line cone from the stopping target. It also adds the CMake build target and FCL prolog defaults. Here are the issues I found in the current code (head commit 5fbd4cb):

---

🔴 Bugs / Correctness Issues

1. Division by zero when seed_dir.z() is near zero

In findTimePeakInCluster, the trajectory parameter dr is computed as:

      const double dz = hit_pos.z() - cl_pos.z(); // distance along z to the hit
      const double dr = dz/seed_dir.z(); // distance along the trajectory to the hit location in z

If seed_dir.z() is zero or very small (e.g., if the calo cluster happens to be at the same z as the stopping target), this produces a division by zero or inf. There is no guard against this. While seed_dir_mag <= 0 is checked earlier, that does not protect against a zero z-component.

Suggestion: Add a guard like if(std::abs(seed_dir.z()) < 1.e-6) return; after normalizing seed_dir, or check std::abs(dz_target) since a zero dz_target implies the same problem.

2. Cone radius formula allows negative dz hits to pass

The cone radius at a hit's z-position is:

      const double cone = std::abs(dz/dz_target) * stopping_target_radius_;

The use of std::abs(dz/dz_target) means that hits behind the calo cluster (on the wrong side of the cone) get a positive cone radius instead of being rejected. For a downstream track, hits with z > cl_pos.z() are beyond the calorimeter and should not be included in the cone selection. The std::abs treats them equivalently to hits between the target and calorimeter.

Suggestion: Either reject hits where dz has the wrong sign relative to the expected track direction, or compute the cone only for the correct direction. For example, for downstream, reject hits with dz > 0 (beyond the calo cluster).

3. Hits are not deduplicated across calo clusters

The produce() method loops over all calorimeter clusters and builds independent TimeCluster objects, but there is no mechanism to prevent the same tracker hit from being assigned to multiple time clusters. The previous version of this code (from the earlier review comment) referenced a hits_used_in_tcs_ set, but the current code has removed it. This means the same hit can appear in multiple output TimeClusters.

Whether this is a bug or intentional depends on downstream consumer expectations. If downstream modules assume exclusive hit assignment, this is a problem.

---

🟡 Design / Robustness Concerns

4. tracker_ geometry pointer is fetched but never used

In beginRun, tracker_ is populated:

    mu2e::GeomHandle<mu2e::Tracker> th;
    tracker_ = th.get();

However, tracker_ is never referenced anywhere else in the module. This is dead code and a dangling pointer that adds confusion.

Suggestion: Remove tracker_ and its initialization unless it is needed for a planned feature.

5. Unused using CLHEP::HepVector declaration

using CLHEP::HepVector;

HepVector is never used in the module. This should be removed.

6. Unused includes

• Offline/DataProducts/inc/Helicity.hh — Helicity is never used.
• <iostream> — no std::cout/std::cerr calls; only printf is used.
• <set> — no std::set is used in the current code.
• Offline/GeometryService/inc/DetectorSystem.hh — DetectorSystem is never used.

Suggestion: Remove unused #include directives to reduce compile dependencies.

7. Empty beginJob() and endJob() methods

Both beginJob() and endJob() have empty bodies. While harmless, they are unnecessary boilerplate since the base class provides defaults.

  void CalLineTimePeakFinder::beginJob() {
  }

8. printf used instead of art's message facility

The module uses printf throughout for diagnostic output. The art framework's mf::LogDebug, mf::LogInfo, mf::LogWarning, and mf::LogError facilities are the standard way to produce output in art modules, and they integrate with the framework's message configuration and filtering. printf output cannot be suppressed or redirected by FCL.

9. No bounds check on cl.diskID()

    cl_pos.setZ(cl_pos.z() + (cl.diskID() == 0 ? calo_d0_offset_ : calo_d1_offset_));

This assumes diskID() can only be 0 or 1. If a cluster ever reports an unexpected disk ID, the code silently uses calo_d1_offset_ as a default, which could produce incorrect results without any warning.

10. Comment typo in seed direction normalization

    if(seed_dir_mag <= 0.) return; // can't define a seed direction, so return false{

The comment says return false{ — the trailing false{ is leftover from an edit. Should just be return; with a corrected comment.

---

⚠️ CI Status

The latest build at 5fbd4cb passed most tests but the trigger test still fails (Return Code 2). This should be investigated — it may indicate a missing trigger path configuration for the new module, or a pre-existing issue on the Run1B branch.

---

Summary

Severity	Issue	Location
🔴 Bug	Division by zero when seed_dir.z() ≈ 0	findTimePeakInCluster (L233)
🔴 Bug	Cone accepts hits behind calo cluster (std::abs on dz)	findTimePeakInCluster (L243)
🔴 Bug	Same hits can be assigned to multiple time clusters (no dedup)	produce / findTimePeakInCluster
🟡 Cleanup	tracker_ fetched but never used	beginRun (L152)
🟡 Cleanup	Unused using HepVector, includes (Helicity, iostream, set, DetectorSystem)	Top of file
🟡 Cleanup	Empty beginJob()/endJob()	L148, L333
🟡 Style	printf instead of art message facility	Throughout
🟡 Robustness	No validation of cl.diskID()	L211
⚪ Typo	Comment has leftover false{	L218
⚠️ CI	Trigger test failing	Build #2940

The most critical issue is #1 — the potential division by zero in the trajectory parameterization when seed_dir.z() is small. Issue #2 (cone geometry accepting wrong-side hits) can silently add incorrect hits to time clusters and degrade reconstruction quality.

[michaelmackenzie]

I addressed a few of the last minor suggestions. Most of the major ones are not real concerns (e.g. the dz to the target will never be very small, the calo cluster disk ID won't be > 1, we're not worried about tracker hits behind the calo, etc.). I'm allowing each calo cluster to form an independent time cluster, so hits can be duplicated for in-time and in-space calo clusters (as is done with TZClusterFinder).