Skip to content

MSubModuleStripReadout: saturated strips are lost, added neighbor strips have no opposite strip, sub-threshold guard ring hits are kept, ADC truncated #281

Open
#290
Open
MSubModuleStripReadout: saturated strips are lost, added neighbor strips have no opposite strip, sub-threshold guard ring hits are kept, ADC truncated#281
#290
Assignees
fhagemann
Labels
DEEDevelopment related to the detector effects engine

Description

@zoglauer
zoglauer
opened on Oct 19, 2026

This issue was written by Claude under the direction of @zoglauer.

The unit test UTNSubModuleStripReadout (PR #247) found five problems in MSubModuleStripReadout:

nuclearizer/src/MSubModuleStripReadout.cxx

Lines 126 to 299 in 988005f

bool MSubModuleStripReadout::AnalyzeEvent(MReadOutAssembly* Event)
{
// Main data analysis routine, which updates the event to a new level
// Get low-voltage and high-voltage hits
for (auto* Hits : { &Event->GetDEEStripHitLVListReference(), &Event->GetDEEStripHitHVListReference() }) {
set<MReadOutElementDoubleStrip> TriggeredStripHits;
set<MReadOutElementDoubleStrip> NeighborCandidateStripHits;
set<MReadOutElementDoubleStrip> NeighborStripHits;
set<MReadOutElementDoubleStrip> DeadStrips;
for (auto SH = Hits->begin(); SH != Hits->end(); ) {
// If the user wants it applied, apply the FWHM Guassian energy resolution
if (m_ApplyResolutionCalibration == true) {
// Look up the FWHM fit for this strip
if (m_ResolutionCalibration.count(SH->m_ROE) == 1) {
double Sigma = m_ResolutionCalibration[SH->m_ROE]->Eval(SH->m_Energy);
// Smear the hit energy using a Gaussian distribution
SH->m_Energy = gRandom->Gaus(SH->m_Energy, Sigma);
// If energy is lower than zero now, floor it to zero
if (SH->m_Energy < 0) {
SH->m_Energy = 0;
}
} else {
// The fit wasn't found! Handle the error
if (g_Verbosity >= c_Warning) {
cout << m_Name << ": Warning - No resolution calibration fit found for strip ID " << SH->m_ROE.GetStripID() << endl;
}
// Note, if no resolution calibration is found then the energy remains unsmeared
}
}
// Apply the inverse energy calibration
// Look up the fit using the ecal
TF1* Fit = m_Calibration[SH->m_ROE];
if (Fit != nullptr) {
// Apply the inverse energy calibration using ROOT's poly inverter (keV -> ADC) in the allowed ADC range
double calculatedADC = Fit->GetX(SH->m_Energy, 0., m_MaxADCRange);
// Apply hardware limits
if (calculatedADC > m_MaxADCRange) calculatedADC = m_MaxADCRange;
if (calculatedADC < 0) calculatedADC = 0;
SH->m_ADC = static_cast<unsigned int>(calculatedADC);
// Apply the hardware threshold to determine if a strip hit is a nearest-neighbor strip or not
SH->m_IsNearestNeighbor = SH->m_ADC < m_HardwareThresholdMap[SH->m_ROE];
// Keep track of all triggered strips and their nearest neighbors
if (SH->m_IsNearestNeighbor == false) {
TriggeredStripHits.insert(SH->m_ROE);
// Only keep track of nearest neighbors for non-GR events
if (SH->m_IsGuardRing == false) {
if (SH->m_ROE.GetStripID() > 0) {
MReadOutElementDoubleStrip Left = SH->m_ROE;
Left.SetStripID(Left.GetStripID() - 1);
NeighborCandidateStripHits.insert(Left);
}
if (SH->m_ROE.GetStripID() < 63) {
MReadOutElementDoubleStrip Right = SH->m_ROE;
Right.SetStripID(Right.GetStripID() + 1);
NeighborCandidateStripHits.insert(Right);
}
}
}
++SH;
} else {
// If no calibration exists in the .ecal file for this strip set it to ADC value of 0
if (g_Verbosity >= c_Warning) cout << m_Name << ": No inverse calibration found for element " << SH->m_ROE << endl;
DeadStrips.insert(SH->m_ROE);
// TODO: Currently, strip hits without a valid energy calibration are erased
// For double-wide / shorted strips, their information should be added
// to the active strip it is connected to
SH = Hits->erase(SH);
}
}
// Iterate through the list once more to remove sub-threshold hits and add
for (auto SH = Hits->begin(); SH != Hits->end(); ) {
// Iterate only through nearest neighbor candidates
if (SH->m_IsNearestNeighbor == true) {
bool HasTriggeredNeighbor = false;
// Check that one of its neighbors has triggered ...
NeighborCandidateStripHits.erase(SH->m_ROE);
if (HasTriggeredNeighbor == false && SH->m_ROE.GetStripID() > 0) {
MReadOutElementDoubleStrip Left = SH->m_ROE;
Left.SetStripID(Left.GetStripID() - 1);
HasTriggeredNeighbor = TriggeredStripHits.count(Left) > 0;
}
if (HasTriggeredNeighbor == false && SH->m_ROE.GetStripID() < 63) {
MReadOutElementDoubleStrip Right = SH->m_ROE;
Right.SetStripID(Right.GetStripID() + 1);
HasTriggeredNeighbor = TriggeredStripHits.count(Right) > 0;
}
if (HasTriggeredNeighbor == true) {
NeighborStripHits.insert(SH->m_ROE);
++SH;
} else {
// ... or remove it otherwise
SH = Hits->erase(SH);
}
}
else {
++SH;
}
}
// Add "empty" strip hits for neighbors that did not trigger
// Assume an empty baseline (Energy == 0) and no timing (Timing == 0)
for (const auto& NSH : NeighborCandidateStripHits) {
// Skip dead strips and triggered strips
if (DeadStrips.count(NSH) > 0) continue;
if (TriggeredStripHits.count(NSH) > 0) continue;
if (NeighborStripHits.count(NSH) > 0) continue;
// Skip strips with no energy calibration, they are most probably dead or shorted
TF1* Fit = m_Calibration[NSH];
if (Fit == nullptr) {
DeadStrips.insert(NSH);
continue;
}
if (g_Verbosity >= c_Info) {
cout << "Event ID " << Event->GetID() << ": Adding nearest neighbor candidate: " << (NSH.IsLowVoltageStrip() ? "LV" : "HV") << NSH.GetStripID() << endl;
}
MDEEStripHit SH;
SH.m_ROE = NSH;
SH.m_Energy = 0;
SH.m_FastPeakTime = 3200; // some value that will result in m_Timing of 1000ns, to pass TAC cuts ...
SH.m_Timing = 0;
SH.m_ADC = 0;
SH.m_TAC = 0;
SH.m_HasFastTiming = false;
SH.m_IsNearestNeighbor = true;
SH.m_IsGuardRing = false;
// Note, if no resolution calibration is found then the energy remains unsmeared
if (m_ApplyResolutionCalibration == true) {
if (m_ResolutionCalibration.count(NSH) == 1) {
double Sigma = m_ResolutionCalibration[NSH]->Eval(SH.m_Energy);
SH.m_Energy = max(gRandom->Gaus(SH.m_Energy, Sigma), 0.0);
}
}
// Apply the inverse energy calibration using ROOT's poly inverter (keV -> ADC) in the allowed ADC range
double calculatedADC = Fit->GetX(SH.m_Energy, 0., m_MaxADCRange);
SH.m_ADC = static_cast<unsigned int>(clamp(calculatedADC, 0.0, m_MaxADCRange));
if (NSH.IsLowVoltageStrip() == true) {
Event->AddDEEStripHitLV(SH);
} else {
Event->AddDEEStripHitHV(SH);
}
NeighborStripHits.insert(NSH);
}
}
return true;
}

1. Saturated strips are lost

For an energy beyond the calibration at m_MaxADCRange, Fit->GetX(Energy, 0, m_MaxADCRange) finds no root and returns NaN ("Interval does not contain a root"). The clamping (L173-L174) does not catch NaN, and static_cast<unsigned int>(NaN) is undefined behavior (0 on arm64) (L170-L176). The strip then counts as below threshold and is removed as a nearest neighbor without triggered neighbor. In the test, a 9000 keV hit (ADC 18000) disappears instead of being saturated at 16383. Energies above the range should be set to m_MaxADCRange. The same applies to the added neighbor hits (L286-L287).

2. Added neighbor strips have no opposite strip ID

The "empty" neighbor hits created for non-triggered neighbors (L266-L275) do not set m_OppositeStripID, which stays uninitialized (in the test: 1833997776 instead of 33, the opposite strip of the triggering strip 10). It should be taken from the triggering strip.

3. A sub-threshold guard ring hit next to a triggered strip is kept

Guard ring hits are excluded from adding nearest-neighbor candidates (L185-L186), but not from being kept as nearest neighbor (L213-L240): a sub-threshold guard ring hit (strip 64) next to triggered strip 63 is kept as neighbor. A kept strip-64 hit later raises the guard ring veto, although the guard ring did not trigger. Sub-threshold guard ring hits should be removed.

4. The ADC value is truncated instead of rounded

static_cast<unsigned int>(calculatedADC) truncates (L176, L287); GetX() returns e.g. 199.99999 for an exact ADC of 200, which gives 199. The result depends on the root-finder precision; lround() would make it stable.

5. Minor: m_ApplyResolutionCalibration is not initialized in the constructor

Minor, since MModuleDEESMEX::Initialize() always sets it before initializing the sub-module (MModuleDEESMEX.cxx L115).

Done when: all problems described in this issue are fixed, all ExpectedToFail() instances related to this issue are removed from the unit tests, and the unit tests pass on both macOS and Linux.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

Labels

DEEDevelopment related to the detector effects engine

Type

No type

Projects

No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions