Updated CRV coincidence cluster finder

[ehrlich-uva]

Changed algorithm that pairs hits of both readout sides of CRV counters.

[FNALbuild]

Hi @ehrlich-uva,
You have proposed changes to files in these packages:

• CRVReco
• RecoDataProducts

which require these tests: build.

@Mu2e/write, @Mu2e/fnalbuild-users have access to CI actions on main.

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

About FNALbuild. Code review on Mu2e/Offline.

[FNALbuild]

☀️ The build tests passed at 98c3b17.

Test	Result	Details
test with	✅	Command did not list any other PRs to include
merge	✅	Merged 98c3b17 at c570352
build (prof)	✅	Log file. Build time: 04 min 05 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 1.
check_cmake	✅	Log file.
FIXME, TODO	✅	TODO (0) FIXME (0) in 3 files
clang-tidy	➡️	2 errors 76 warnings
whitespace check	✅	no whitespace errors found

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

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 examining its details and the code changes.Here is my detailed review of PR #1772. This is a significant refactoring of the CRV coincidence cluster finder that changes the algorithm for pairing hits from both readout sides of CRV counters. Below are the issues I found, organized by severity.

---

🔴 Bugs / Potential Correctness Issues

1. Operator precedence bug in clusterProperties — division vs. map access

In CRVReco/src/CrvCoincidenceFinder_module.cc, the new code computes PE-weighted average times like this:

double hitTime0=hitTimes[0].begin()->second/hitPEs[0].begin()->second;
double hitTime1=hitTimes[1].rbegin()->second/hitPEs[1].rbegin()->second;

This divides only the ->second values, which is correct in terms of C++ operator precedence. However, the hitPEs and hitTimes maps are keyed by float (longitudinal position), and entries are accumulated via +=. The code picks begin() for side 0 and rbegin() for side 1 to find "the longitudinal positions farthest away from each other." But if multiple counters contribute to the same readout side at different longitudinal positions, only the first/last map entry is used — the PE-weighted time from other entries is silently ignored. The comment even acknowledges this ("the situation is actually overconstraint, we ignore the other positions"). This could yield incorrect avgHitTime and avgHitPos for clusters spanning multiple counters on the same side.

2. filterHits now only compares hits from the same readout side when computing PE sums for adjacent counters

In the old code, filterHits iterated over all hits (both SiPM sides) to accumulate PEs from adjacent counters. In the new code, filterHits is called separately per readout side, and iterHitAdjacent iterates over the same hitsFiltered list (which is side-specific). This means the PE sum for the PE threshold check now only includes hits from one readout side. If the threshold logic previously relied on combining PEs from both sides of the same counter, this could reject hits that should pass.

3. Using float as std::map key for hitPEs and hitTimes

std::array<std::map<float,double>, CRVId::nSidesPerBar> hitPEs;
std::array<std::map<float,double>, CRVId::nSidesPerBar> hitTimes;

Using float as a map key is fragile because floating-point equality comparisons can fail due to precision. Two hits that should logically be at the same longitudinal position (same counter, same side) could end up in different map entries if there are tiny floating-point rounding differences in the computation longitudinalPos = counterPos[lengthDirection] ± halfLength. This would corrupt the PE-weighted time averaging.

4. Potential dangling reference: CrvCoincidence.hh still referenced in CalPatRec/inc/HlPrint.hh

The PR deletes RecoDataProducts/inc/CrvCoincidence.hh and removes it from classes.h and classes_def.xml. However, CalPatRec/inc/HlPrint.hh still has a forward declaration:

  class CrvCoincidence;

While a forward declaration alone won't cause a build error (unless it's used somewhere), this is dead code that references a now-deleted class. If any downstream code tries to #include and use CrvCoincidence, it will fail. A search shows this reference may appear in other files too.

---

🟡 Logic / Design Concerns

5. One-side-only fallback can overwrite with the wrong value

In the single-readout-side fallback:

if(!hitTimes[0].empty())
{
  // ... computes avgHitTime from side 0
}
if(!hitTimes[1].empty())
{
  // ... computes avgHitTime from side 1 — overwrites the side-0 result!
}

Both if blocks are independent (not else if). If somehow both sides have hits but the earlier two-side check failed (because hitPosOrigin was outside the counter), execution falls through to the else branch (one-side fallback), where both if blocks execute. The second if silently overwrites avgHitTime computed from side 0. This should likely be else if to use only one side, or the logic should be restructured — though in practice the else branch should only be reached when one side is empty.

6. clusterMaxTimeDifference recomputed by scanning _sectorMap every iteration

Inside the per-sector-type loop, clusterMaxTimeDifference is computed by scanning all entries of _sectorMap:

float clusterMaxTimeDifference=0;
for(auto const &sectorMapEntry : _sectorMap)
{
  if(sectorMapEntry.second.sectorType==crvSectorType && sectorMapEntry.second.maxTimeDifference>clusterMaxTimeDifference)
    clusterMaxTimeDifference=sectorMapEntry.second.maxTimeDifference;
}

This is repeated for maxHalfLength and readout side detection as well. This is O(N²) in the number of sector types × sector map size. It would be cleaner to precompute these values once per sector type.

---

🟢 Minor / Style Issues

7. Typos in comments

• "diving" should be "dividing" (appears multiple times):

  //find avg PE-weighted time by diving by total PEs
  

• "orign" should be "origin":

  //calculate the time of the hit orign from both sides
  

• "overconstraint" should be "overconstrained":

  //if we have more than one position at a readout side, the situation is actually overconstraint
  

8. Mixed float and double types

The code inconsistently uses float for some variables (clusterMaxTimeDifference, maxHalfLength, timeDifferenceOppositeSides) and double for others in the same computation chain. While not a bug, this could lead to subtle precision losses and makes the code harder to reason about.

9. readoutSide parameter type mismatch

filterHits takes unsigned int readoutSide, but it's compared against SiPM%CRVId::nSidesPerBar which is int % int = int. The implicit signed-to-unsigned comparison is safe here but could be made explicit.

---

Summary

Severity	Issue	Location
🔴 Bug	float map keys can split same-counter hits	clusterProperties
🔴 Bug	filterHits PE threshold now only checks one side	filterHits
🔴 Bug	Deleted CrvCoincidence.hh still forward-declared	HlPrint.hh
🟡 Logic	One-side fallback can silently overwrite avgHitTime	clusterProperties
🟡 Logic	Multi-counter clusters use only first/last map entry	clusterProperties
🟡 Perf	Repeated _sectorMap scans per sector type	produce()
🟢 Style	Typos in comments ("diving", "orign", "overconstraint")	Throughout
🟢 Style	Inconsistent float/double mixing	Throughout

[ehrlich-uva]

Responses:

1. This is the desired result.
2. This is the desired result.
3. Fixed.
4. Will be done at a later time.
5. This situation (hits at both counter sides in the else branch) won't happen.
6. Done.
7. Fixed.
8. Will be addressed at a later time.
9. No change required.