StripEnergyThresholdFinder for per-strip slow and fast threshold extraction with diagnostics

[JarredMRoberts]

Adds a standalone application, StripEnergyThresholdFinder, for computing per-strip slow and fast energy thresholds. Slow thresholds are determined from ADC spectra using a noise peak and trough method, while fast thresholds are determined from dt0/dt1 timing crossover. The tool reads calibrated HDF5 data via MModuleLoaderMeasurementsHDF, applies strip mapping and energy calibration from a YAML configuration, and produces ROOT diagnostic outputs (energy spectra with thresholds, dt0 vs dt1 per strip, and threshold distributions) along with CSV export files. The implementation is self-contained under apps/ and does not modify existing modules. Tested on COSI datasets with consistent threshold behavior and expected diagnostic results. Target branch is develop/em.

[JarredMRoberts]

I still need to fix all of the code style issues and work on some optimizations to speed the code up a bit.

[fhagemann]

Is this different from #143 or making #143 obsolete?

[fhagemann]

Hi Jarred,
I had an extensive look at this app and finally got this to run.
I replied to all comments in this PR, and opened a separate PR onto your fork/branch, addressing some code changes to

1. get this app running
2. remove the helper classes by replacing them with existing nuclearizer/megalib code

Here is the PR onto your branch with detailed code changes: JarredMRoberts#1

[zoglauer]

I second not to have an additional parser linked.
You can use the nuclearizer/megalib XML file format instead; your file then would be:

<strip-threshold-finder>
    <input>
        <data_files>
            <file>/path/to/data/file.hdf5</file>
        </data_files>
        <calibration_file>/path/to/calibration/file.ecal</calibration_file>
        <tac_calibration_file>/path/to/tac/calibration/file.csv</tac_calibration_file>
        <strip_map>/path/to/strip/map/file.map</strip_map>
    </input>

    <analysis>
        <!-- number of skipped strips with minimum statistics -->
        <min_entries>10</min_entries> 
        
        <!-- If a threshold cannot be determined set threshold to default value -->
        <fallback_threshold_keV>20</fallback_threshold_keV> 
        
        <!-- limit for locating the low-energy noise peak -->
        <!-- most noise peaks should be between 100 and 250 -->
        <!-- Worst case the ADC max should be set to ~1000 -->
        <noise_search_max_adc>1800</noise_search_max_adc>
    </analysis>
    
    <output>
        <prefix>output_file_prefix</prefix>
    </output>
</strip-threshold-finder>

[fhagemann]

One more remark: I got this app to run and got reasonable results for the fast and slow thresholds only when using a dataset taken with NN off.
We might want to always filter out NN events, for this app to also run using datasets taken with NN on.

[fhagemann]

This is a quick review with code-style comments, I will follow up with some more general comments

[fhagemann]

This is a quick review with code-style comments, I will follow up with some more general comments

[fhagemann]

Some general comments:

• do you have some figures to walk us through the pedestal/peak finding algorithms, to get an idea in what cases the pedestal, flattop or peak detection loops would do their thing?
• we might want to avoid having hard-coded values distributed throughout the file and define them as constants at the top of the file (or member variables of the threshold app class) instead.
• In the end, we might want to generate ONE threshold file per threshold type (fast, slow, hardware), instead of splitting them into LV and HV subfiles. I believe that the forward pipeline will require a combined file in the end anyhow

[fhagemann]

Is this ever not the case?
Also, should this hard-coded value become some member variable or be defined at the very beginning of the file, instead of very far into the code?

[fhagemann]

Looking at this in more detail: shouldn't this be:

Suggested change

	// --- dt0 vs dt1 separation ---
	bool is_dt1 = (TAC > 8000); // initial threshold
	// --- dt0 vs dt1 separation ---
	bool is_dt1 = SH->HasFastTiming();

??

[fhagemann]

Tested on data from HP52406-1 (Am-241 lifetime data, taken with HDFv2.2: gse_20260217T115220.hdf5):

With the current TAC > 8000 line, I get nothing in HV strip 32 (no dt1 distribution, everything seems to be classified as dt0), resulting in the fallback fast threshold:

WARNING: FAST threshold could not be determined for Det 0 Side HV Strip 32: insufficient timing populations (dt0=88441, dt1=293, dt1 fraction=0.3302%). Using fallback threshold = 34 keV.

Changing it to SH->HasFastTiming() results in reasonable dt0 and dt1 distributions, but the fast threshold seems to be determined incorrectly:

This is my XML config file (using strip map and ecal file from the unit tests):

<?xml version="1.0" encoding="UTF-8"?>
<StripEnergyThresholdFinder>
  <Input>
    <DataFiles>
      <DataFile>/home/hagemann/Downloads/gse_20260217T115220.hdf5</DataFile>
    </DataFiles>
    <CalibrationFile>/home/hagemann/Software/COSItools/nuclearizer/resource/unittestdata/406-1/hp52406-1.full.ecal</CalibrationFile>
    <StripMap>/home/hagemann/Software/COSItools/nuclearizer/resource/unittestdata/406-1/hp52406-1.stripmap.map</StripMap>
  </Input>
  <Output>
    <Prefix>output_file_prefix</Prefix>
  </Output>
  <Analysis>
    <MinEntries>10</MinEntries>
    <FallbackThresholdKeV>20</FallbackThresholdKeV>
    <NoiseSearchMaxKeV>40</NoiseSearchMaxKeV>
    <FastFallbackKeV>34</FastFallbackKeV>
  </Analysis>
</StripEnergyThresholdFinder>

[JarredMRoberts]

Whatever was broken here is fixed now. I've never been able to reproduce this behavior.

[fhagemann]

Unresolving this, because I don't think this is true.

When using your current is_dt1 classification, you sometimes get ONLY dt0 events (see your Notes for HV strip 32 for all detectors tested with the faulty NICEv1 board).

I get reasonable dt0 and dt1 distributions ONLY if I incorporate the change listed above:

        // --- dt0 vs dt1 separation ---
        bool is_dt1 = SH->HasFastTiming();

[fhagemann]

--> 95b5a46

[fhagemann]

I can confirm that the app compiles and works when using multiple files.
You might want to remove the index from the XML file (it also works without it).

[JarredMRoberts]

I've put together some slides that go into some detail about the app, what the inputs and outputs are, and how the various algorithms work.

(https://docs.google.com/presentation/d/1_sveTiA17fUNNXd2h1dcrlFbWuK2XNRssNNchlbn5ZY/edit?usp=sharing)

README_ Nuclearizer Threshold Finder.pdf

[fhagemann]

Thanks for creating hardware threshold files and uploading them here: https://drive.google.com/drive/folders/1FN4fp0_imv-opuJHoKn1N3DPwUUPe38A?usp=sharing

I opened PR #208 based on those, but had to make some minor adjustments to be consistent with the format of dummy_hardware_thresholds.csv:

1. adding the STRIP_ID from the StripMap to the first column,
2. printing all HV strips first and then all LV strips (this app currently alternates, giving HV0, LV0, HV1, LV1, etc.) -- this one is purely cosmetic, the current format would still work with the threshold file parser.

There have also been some changes to the TAC calibration module in PR #194, and I will take care of resolving the merge conflicts here.

[fhagemann]

I still believe this is wrong and should be

Suggested change

	// --- dt0 vs dt1 separation ---
	bool is_dt1 = (TAC > FAST_TIMING_TAC_THRESHOLD); // initial threshold
	// --- dt0 vs dt1 separation ---
	bool is_dt1 = SH->HasFastTiming();

[fhagemann]

In your notes, you show these spectra for some detectors with HV strip 32 (the ones tested with the faulty NICEv1 board), which I can reproduced with the current is_dt1 = (TAC > FAST_TIMING_TAC_THRESHOLD);:

When I update the line to bool is_dt1 = SH->FastTiming();, this is what I get:

So, using bool is_dt1 = SH->FastTiming(); seems to give a more reasonable classification into dt0 and dt1 events, but the fit result for the fast threshold seems to be off.

[fhagemann]

I plotted this histogram from the same datafile using Julia, and also plotted the fast threshold that's determined with this app (1901 ADC units) in ADC space, and this is consistent with what we are getting in this app:

Something causes the fast threshold to be wrong here, I will try to find out what it is.

[fhagemann]

Ok, I found the reason:
The fast threshold finding algorithm finds the first bin where there are more dt1 events (with fast timing) than dt0 events (with slow timing), but also requires there to be at least 50 events in that bin in general.
(constexpr int FAST_MIN_CROSSOVER_COUNTS = 50;).

What happens here, is that this default of 50 is too low, so it sets the threshold at the first bin where the total number exceeded 50:

We might want to add some check that look backwards if the previous bin (regardless if it had less than 50 counts) had more dt0 than dt1 events, or shift the fast threshold value back down if that isn't the case.

[fhagemann]

This seems to work (at least in Julia, would need to translate this into C++):

using StatsBase

dt1 = fit(Histogram, is_dt1.energy, 1600:1:2100)
dt0 = fit(Histogram, is_dt0.energy, 1600:1:2100)

for i in eachindex(dt0.weights)
    
    # Initially ignore all bins that have fewer than 50 counts
    if dt0.weights[i] + dt1.weights[i] < 50
        continue
    end
    
    # Find the first bin with more than 50 counts with more dt1 than dt0 events
    if dt0.weights[i] < dt1.weights[i]
        println("Initial guess:   ", first(dt0.edges)[i])
        
        # If that was also the case for previous bins: reduce i
        while dt0.weights[i-1] < dt1.weights[i-1]
            i -= 1
        end
        
        println("Backward search: ", first(dt0.edges)[i])
        break
    end
end

Initial guess:   1901
Backward search: 1793

[fhagemann]

--> 95b5a46 and 7dce517

[fhagemann]

Suggested change

	if (n1 > n0) {
	crossoverADC = ADC_val;
	break;
	}
	if (n1 > n0) {
	crossoverADC = ADC_val;
	// Make sure that this is an actual crossover and not just the first bin
	// with more than FAST_MIN_CROSSOVER_COUNTS where n1 was above n0
	while (n1 > n0 && ADC_val >= first_nonzero) {
	crossoverADC = ADC_val;
	ADC_val--;
	if (ADCMap.find(ADC_val) != ADCMap.end()) {
	n0 = ADCMap[ADC_val].first;
	n1 = ADCMap[ADC_val].second;
	}
	}

Adding this backward search seems to resolve the issue above:

[fhagemann]

The vimdiff between the fast threshold file created (left) with and (right) without the backward search suggests that this might impact that results for the fast thresholds for most HV strips

[fhagemann]

For this specific measurement, it looks like the HV spectra seem to feature this dip with <50 counts/bin, such that we might want to employ the backward search:

The LV spectra seem fine, because they exceed 50 counts/bin in the region where we have the crossover:

All major changes I requested are implemented.

[fhagemann]

When trying this out for Cs-137 data to explore the possibility to use Cs-137 data in payload calibrations for threshold determination, I found that the current implementation would give wrong results if there are a lot of events with fast timing in the noise peak (I guess these are NN events featuring a dominant peak-dip structure, that might give a fast-shaper signal that exceeds the threshold, but in the slow shaper average out to a small energy value).

In this case, the condition n1 .> n0 is reached early:

A fix for this might be to start the crossover search only for ADC bins right of the noise peak, i.e. the slow threshold (plus some tolerance):

Suggested change

	for (auto& a : ADCMap) {
	int ADC_val = a.first;
	int n0 = a.second.first;
	int n1 = a.second.second;
	// Require minimum statistics to avoid noise triggers
	if (n0 + n1 < FAST_MIN_CROSSOVER_COUNTS) {
	continue;
	}
	if (n1 > n0) {
	crossoverADC = ADC_val;
	// Make sure that this is an actual crossover and not just the first bin
	// with more than FAST_MIN_CROSSOVER_COUNTS where n1 was above n0
	while (n1 > n0 && ADC_val >= first_nonzero) {
	crossoverADC = ADC_val;
	ADC_val--;
	if (ADCMap.find(ADC_val) != ADCMap.end()) {
	n0 = ADCMap[ADC_val].first;
	n1 = ADCMap[ADC_val].second;
	}
	}
	break;
	}
	}
	const int nbins = FAST_CROSSOVER_WINDOW_BINS;
	// Extend search window for stability
	int searchMax = first_nonzero + FAST_SEARCH_RANGE_ADC;
	for (int ADC = first_nonzero; ADC < searchMax; ADC++) {
	int n_dt0 = 0;
	int n_dt1 = 0;
	for (int a = ADC; a < ADC + nbins; a++) {
	if (ADCMap.find(a) == ADCMap.end()) {
	continue;
	}
	n_dt0 += ADCMap[a].first;
	n_dt1 += ADCMap[a].second;
	}
	for (auto& a : ADCMap) {
	int ADC_val = a.first;
	int n0 = a.second.first;
	int n1 = a.second.second;
	// Require minimum statistics to avoid noise triggers
	if (n0 + n1 < FAST_MIN_CROSSOVER_COUNTS || ADC_val <= m_SlowThresholdsADC[R] + FAST_FIRST_ADC_OFFSET) {
	continue;
	}
	if (n1 > n0) {
	crossoverADC = ADC_val;
	// Make sure that this is an actual crossover and not just the first bin
	// with more than FAST_MIN_CROSSOVER_COUNTS where n1 was above n0
	while (n1 > n0 && ADC_val >= m_SlowThresholdsADC[R] + FAST_FIRST_ADC_OFFSET) {
	crossoverADC = ADC_val;
	ADC_val--;
	if (ADCMap.find(ADC_val) != ADCMap.end()) {
	n0 = ADCMap[ADC_val].first;
	n1 = ADCMap[ADC_val].second;
	}
	}
	break;
	}
	}
	const int nbins = FAST_CROSSOVER_WINDOW_BINS;
	// Extend search window for stability
	int searchMax = m_SlowThresholdsADC[R] + FAST_SEARCH_RANGE_ADC;
	for (int ADC = m_SlowThresholdsADC[R]; ADC < searchMax; ADC++) {
	int n_dt0 = 0;
	int n_dt1 = 0;
	for (int a = ADC; a < ADC + nbins; a++) {
	if (ADCMap.find(a) == ADCMap.end()) {
	continue;
	}
	n_dt0 += ADCMap[a].first;
	n_dt1 += ADCMap[a].second;
	}

(and I also had to set FAST_FIRST_ADC_OFFSET to 20 for this specific dataset I was looking at).

This is what I get with the updated code (but note that FAST_FIRST_ADC_OFFSET might need to be optimized):

[fhagemann]

They way the plotting code is written right now does ignore detector ID.
When running this threshold app on payload calibration data, all detectors will be plotted into one plot:

Can this be split into multiple plots (maybe one per "active" detector)?

[fhagemann]

I had to run this app now several times on the same data, and the part that takes the longest is reading the HDF5 files and applying the strip map / energy calibration. Would it be possible to also support reading ROA files as input (and not just HDF5) to avoid this first step in the threshold app having to be re-run every time I want to change settings or fix bugs in the app.
This becomes especially relevant for payload data, where there is a LOT of data (and we might be able to speed this up per detector by creating detector specific ROA files (#231)).

[fhagemann]

A lot of the plots and the analysis are hard-coded to the Am-241 energy range, which makes the fast threshold extraction from the FM calibration Cs-137 partyhat measurement challenging. It would be nice to have the option to set some of the constants (e.g. constexpr int FAST_SEARCH_RANGE_ADC = 2000;) and some of the member class variables (like double m_NoiseStartMinCounts = 8.0;`), (or preferably all?) through the config file or terminal commands.