From 0d05ff5881b9bf8a9798634f257c25178a3600d0 Mon Sep 17 00:00:00 2001 From: JohannMartyn <133750974+JohannMartyn@users.noreply.github.com> Date: Fri, 19 Jun 2026 20:13:22 +0200 Subject: [PATCH 1/4] Correct local baseline subtraction Previously we compared the ADC sample value with and index value. This is corrected now and we compare an index value with an index value. --- UserTools/PhaseIIADCCalibrator/PhaseIIADCCalibrator.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/UserTools/PhaseIIADCCalibrator/PhaseIIADCCalibrator.cpp b/UserTools/PhaseIIADCCalibrator/PhaseIIADCCalibrator.cpp index bcca11115..c6e6f240a 100644 --- a/UserTools/PhaseIIADCCalibrator/PhaseIIADCCalibrator.cpp +++ b/UserTools/PhaseIIADCCalibrator/PhaseIIADCCalibrator.cpp @@ -757,13 +757,15 @@ PhaseIIADCCalibrator::make_calibrated_waveforms_ze3ra_multi( } std::vector cal_data; const std::vector& raw_data = raw_waveform.Samples(); - for (const auto& asample: raw_data){ + const int raw_data_size = static_cast(raw_data.size()); + for(int s=0; s(asample) - baselines.at(j)) * ADC_TO_VOLT); break; - } else if (asample >= RepresentationRegion.back()){ + } else if (s >= RepresentationRegion.back()){ cal_data.push_back((static_cast(asample) - baselines.back()) * ADC_TO_VOLT); break; From 7aae8e427e6293d65ee111f9c72db8e4cf391d61 Mon Sep 17 00:00:00 2001 From: JohannMartyn <133750974+JohannMartyn@users.noreply.github.com> Date: Fri, 19 Jun 2026 20:22:03 +0200 Subject: [PATCH 2/4] Corrected hit time calculation Previously the adc_threshold was subtracted from the peak value. Now the baseline is subtracted. --- UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp b/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp index c77ee2543..ca3f69d2f 100755 --- a/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp +++ b/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp @@ -1046,8 +1046,9 @@ std::vector PhaseIIADCHitFinder::find_pulses_bythreshold( // New approach to hit timing to avoid 2ns bins - 50% threshold above baseline // Look for where the ADC value crosses 50% of the maximum, assign hit time // TODO: consider using an approach recommended by Bob: get time at 50%, get time at 20%, draw straight line in between and find time to zero threshold + double pulse_baseline = calibrated_minibuffer_data.GetBaseline(); const double threshold_percentage = 0.5; - unsigned short threshold_value = ((max_ADC - adc_threshold) * threshold_percentage) + adc_threshold; + unsigned short threshold_value = ((max_ADC - pulse_baseline) * threshold_percentage) + pulse_baseline; double hit_time = peak_sample; bool hit_time_found = false; @@ -1078,7 +1079,6 @@ std::vector PhaseIIADCHitFinder::find_pulses_bythreshold( std::vector trace_y; double pulse_start_time = pulse_start_sample * NS_PER_ADC_SAMPLE; - double pulse_baseline = calibrated_minibuffer_data.GetBaseline(); for (size_t p = pulse_start_sample; p <= pulse_end_sample; ++p) { double ns_time = p * NS_PER_ADC_SAMPLE; From 206793fa75b222bafb246ce1fb15afcc1e0d338a Mon Sep 17 00:00:00 2001 From: JohannMartyn <133750974+JohannMartyn@users.noreply.github.com> Date: Fri, 19 Jun 2026 20:40:40 +0200 Subject: [PATCH 3/4] Added the "NoDubleHits" option for the PhaseIIADCHitFinder This method works similar the the "Fixed_2023_Gains" method, but it fully avoids the incorrect double hit counting. --- .../PhaseIIADCHitFinder.cpp | 201 +++++++++++++++++- 1 file changed, 200 insertions(+), 1 deletion(-) diff --git a/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp b/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp index ca3f69d2f..66e5f376a 100755 --- a/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp +++ b/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp @@ -1117,8 +1117,207 @@ std::vector PhaseIIADCHitFinder::find_pulses_bythreshold( trace_x, trace_y); } } + } + else if(pulse_window_type == "NoDoubleHits"){ + + size_t pulse_start_sample = BOGUS_INT; + size_t pulse_end_sample = BOGUS_INT; + size_t prev_pulse_end_sample = 0; + + // PMT Timing offsets + double timing_offset=0.0; + std::map::const_iterator it = ChannelKeyToTimingOffsetMap.find(channel_key); + if(it != ChannelKeyToTimingOffsetMap.end()){ //Timing offset is available + timing_offset = ChannelKeyToTimingOffsetMap.at(channel_key); + } else { + if(verbosity>v_error && !mc_waveforms){ + std::cout << "PhaseIIADCHitFinder: Didn't find Timing offset for channel " << channel_key << "... setting this channel's offset to 0ns" << std::endl; + } + } + + // loop through samples until we find a pulse, then extract pulse parameters + for (size_t s = 0; s < num_samples; ++s) { + + // if any values are above threshold, we have found a pulse + if ( !in_pulse && raw_minibuffer_data.GetSample(s) > adc_threshold ) { + in_pulse = true; + if(verbosity>v_debug) std::cout << "PhaseIIADCHitFinder: FOUND PULSE" << std::endl; + + // Determine the pulse start by walking back from the threshold crossing + // to the point where the signal was last at/below baseline+1sigma, then + // take 5 samples before that as the start (margin). + if (static_cast(s) - 5 < 0) { + if(verbosity>v_debug) std::cout << "PhaseIIADCHitFinder: Pulse found is VERY EARLY in the minibuffer (< 5 samples)... assigning pulse start as 0" << std::endl; + pulse_start_sample = 0; + } else { + size_t pulsewinleft = s; + bool found_baseline = false; + + while (pulsewinleft > 0) { + double raw_sample_height = raw_minibuffer_data.GetSample(pulsewinleft); + if (raw_sample_height <= baseline_plus_one_sigma) { + pulse_start_sample = (pulsewinleft >= 5) ? (pulsewinleft - 5) : 0; // avoid underflow + found_baseline = true; + break; + } + pulsewinleft--; + } + + // if we walked all the way back to 0 without finding a baseline + // crossing (maybe ringing?), assign pulse start as 0 + // TODO: figure out what is wrong with these pulses + if (!found_baseline) { + if(verbosity>v_debug) std::cout << "PhaseIIADCHitFinder: Baseline crossing was not found... assigning pulse start as 0" << std::endl; + pulse_start_sample = 0; + } + } + + // Clamp the start AFTER it has been computed + // so the window can never overlap the previous pulse. + if (prev_pulse_end_sample > 0 && pulse_start_sample <= prev_pulse_end_sample) { + if(verbosity>v_debug) std::cout << "PhaseIIADCHitFinder: Pulse start would overlap the previous pulse. Clamping to last pulse end + 1" << std::endl; + pulse_start_sample = prev_pulse_end_sample + 1; + } + // once we reach the baseline again (right side of the pulse), we determine the pulse stop (5 samples after baseline + sigma crossing) + // or in case we've reached the end of the minibuffer, force pulse to end + } else if ( in_pulse && ((raw_minibuffer_data.GetSample(s) < baseline_plus_one_sigma) || (s == num_samples - 1)) ) { + + if (s == num_samples - 1) { + if(verbosity>v_debug) std::cout << "PhaseIIADCHitFinder: Pulse found is VERY LATE in the minibuffer (we're at the final sample)... forcing pulse to end" << std::endl; + pulse_end_sample = s; + } else { + pulse_end_sample = (s + 5 < (num_samples - 1)) ? (s + 5) : (num_samples - 1); // ensure we don't exceed the buffer + } + + // double check that pulse start and stop were found successfully + if (verbosity > v_debug) { + std::cout << "PhaseIIADCHitFinder: Pulse start and end determined! (" + << pulse_start_sample << ", " << pulse_end_sample << ")" << std::endl; + } + + const size_t guard = 5; + bool above_thresh = false; // any sample back above trigger within guard + for (size_t g = 1; g <= guard && (s + g) < num_samples; ++g) { + unsigned short val = raw_minibuffer_data.GetSample(s + g); + if (val > adc_threshold) { + above_thresh = true; + break; + } + } + + // We keep moving this pulse until we have well seperated pulses. + // This method is more stable and thus easier to calibrate against MC. + // But we are potentially losing out on multi-hits that could otherwise + // be resolved. For the charge this plays no role, but for timing maybe. + if (above_thresh) { + continue; + } + + s = pulse_end_sample; + prev_pulse_end_sample = pulse_end_sample; + in_pulse = false; + + unsigned long raw_area = 0; // ADC * samples + unsigned short max_ADC = std::numeric_limits::lowest(); + size_t peak_sample = BOGUS_INT; + for (size_t p = pulse_start_sample; p <= pulse_end_sample; ++p) { + raw_area += raw_minibuffer_data.GetSample(p); + if (max_ADC < raw_minibuffer_data.GetSample(p)) { + max_ADC = raw_minibuffer_data.GetSample(p); + peak_sample = p; + } + } + + // The amplitude of the pulse (V) + double calibrated_amplitude = calibrated_minibuffer_data.GetSample(peak_sample); + + // Calculated the charge detected in this pulse (nC) + // using the calibrated waveform + double charge = 0.; + // Integrate the calibrated pulse (to get a quantity in V * samples) + for (size_t p = pulse_start_sample; p <= pulse_end_sample; ++p) { + charge += calibrated_minibuffer_data.GetSample(p); + } + + // Convert the pulse integral to nC + // FIXME: We need a static database with each PMT's impedance + charge *= NS_PER_ADC_SAMPLE / ADC_IMPEDANCE; + // TODO: consider adding code to merge pulses if they occur + // very close together (i.e. if the end of one is just a few samples away + // from the start of another) + + // New approach to hit timing to avoid 2ns bins - 50% threshold above baseline + // Look for where the ADC value crosses 50% of the maximum, assign hit time + // TODO: consider using an approach recommended by Bob: get time at 50%, get time at 20%, draw straight line in between and find time to zero threshold + const double threshold_percentage = 0.5; + double pulse_baseline = calibrated_minibuffer_data.GetBaseline(); + unsigned short threshold_value = ((max_ADC - pulse_baseline) * threshold_percentage) + pulse_baseline; + + double hit_time = peak_sample; + bool hit_time_found = false; + + // Find the first sample (walking back from the peak) below the 50% level + for (size_t p = peak_sample; p > pulse_start_sample; --p) { + if (raw_minibuffer_data.GetSample(p) < threshold_value) { + hit_time = p; + hit_time_found = true; + break; + } + } + + // Perform simple linear interpolation to find exact crossing point + if (hit_time_found) { + if(verbosity>v_debug) std::cout << "Interpolating hit time..." << std::endl; + if (hit_time > pulse_start_sample && hit_time < pulse_end_sample) { + double x1 = hit_time; + double x2 = hit_time + 1.0; + unsigned short y1 = raw_minibuffer_data.GetSample(static_cast(x1)); + unsigned short y2 = raw_minibuffer_data.GetSample(static_cast(x2)); + if (y2 != y1) { // guard against divide-by-zero on a flat segment + hit_time = x1 + (threshold_value - y1) * (x2 - x1) / static_cast(y2 - y1); // linear interpolation + } + } + } - + // extract the x and y points of the pulse (subtract off baseline and "zero" the pulse to the pulse start) + std::vector trace_x; + std::vector trace_y; + double pulse_start_time = pulse_start_sample * NS_PER_ADC_SAMPLE; + + for (size_t p = pulse_start_sample; p <= pulse_end_sample; ++p) { + double ns_time = p * NS_PER_ADC_SAMPLE; + double val_adc = raw_minibuffer_data.GetSample(p); + trace_x.push_back(ns_time - pulse_start_time); + trace_y.push_back(val_adc - pulse_baseline); + } + + if(verbosity>v_debug) std::cout << "PhaseIIADCHitFinder: Hit time [ns] " << hit_time * NS_PER_ADC_SAMPLE << std::endl; + + if (hit_time < 0.0) { + if(verbosity>v_debug) std::cout << "PhaseIIADCHitFinder: Hit time is negative! Defaulting to peak time" << std::endl; + hit_time = peak_sample; + } + + if(verbosity>v_debug) { + std::cout << "PhaseIIADCHitFinder: Pulse properties: " << std::endl; + std::cout << " chanID: " << channel_key << std::endl; + std::cout << " charge: " << ( charge ) << std::endl; + std::cout << " start time: " << ( pulse_start_sample ) << std::endl; + std::cout << " hit time: " << ( hit_time ) << std::endl; + std::cout << " stop time: " << ( pulse_end_sample ) << std::endl; + } + + // Store the freshly made pulse in the vector of found pulses + pulses.emplace_back(channel_key, + ( pulse_start_sample * NS_PER_ADC_SAMPLE )-timing_offset, + ( hit_time * NS_PER_ADC_SAMPLE )-timing_offset, // interpolated hit time + calibrated_minibuffer_data.GetBaseline(), + calibrated_minibuffer_data.GetSigmaBaseline(), + raw_area, max_ADC, calibrated_amplitude, charge, + trace_x, trace_y); + + } + } // ****************************************************************** // Peak windows are defined only by crossing and un-crossing of ADC threshold } else if(pulse_window_type == "dynamic"){ From e59264a19f8684990aca69625b4b6e3f542312da Mon Sep 17 00:00:00 2001 From: JohannMartyn Date: Fri, 2 Oct 2026 09:03:15 +0200 Subject: [PATCH 4/4] Made the NoDoubleHits the default selection in PhaseIIADCHitFinder.cpp --- UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp | 2 +- UserTools/PhaseIIADCHitFinder/README.md | 10 ++++++---- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp b/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp index 66e5f376a..1b6f3c71c 100755 --- a/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp +++ b/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp @@ -18,7 +18,7 @@ bool PhaseIIADCHitFinder::Initialise(std::string config_filename, DataModel& dat adc_threshold_db = "none"; default_adc_threshold = 7; threshold_type = "relative"; - pulse_window_type = "Fixed_2023_Gains"; + pulse_window_type = "NoDoubleHits"; pulse_window_start_shift = -3; pulse_window_end_shift = 25; adc_window_db = "none"; //Used when pulse_finding_approach="fixed_windows" diff --git a/UserTools/PhaseIIADCHitFinder/README.md b/UserTools/PhaseIIADCHitFinder/README.md index 0e8de25b2..eabbc3635 100755 --- a/UserTools/PhaseIIADCHitFinder/README.md +++ b/UserTools/PhaseIIADCHitFinder/README.md @@ -83,10 +83,12 @@ is manipulable using DefaultADCThreshold and DefaultThresholdType config variabl relative to the calibrated baseline ("relative"), or absolute ADC counts ("absolute"). - PulseWindowType [string]: If using "threshold" on pulse finding approach, this toggle defines - how the pulse windows in a waveform are found. There are three options: fixed window ("fixed"), + how the pulse windows in a waveform are found. There are multiple options: fixed window ("fixed"), dynamic window where the pulse windows are defined by crossing and un-crossing threshold ("dynamic"), and ("Fixed_2023_Gains") which implements the same integration window used in the 2023 Gains calibration where the pulse windows are defined by crossing and un-crossing the baseline. + The fourth option is "NoDoubleHits" which follows the approach from "Fixed_2023_Gains", but corrects a bug + there to make sure that pulses are not double counted. - PulseWindowStart [int]: Start of pulse window relative to when adc trigger threshold was crossed. Only used when PulseFindingApproach==threshold and @@ -117,7 +119,7 @@ verbosity 0 UseLEDWaveforms 0 PulseFindingApproach threshold -PulseWindowType Fixed_2023_Gains +PulseWindowType NoDoubleHits DefaultADCThreshold 7 DefaultThresholdType relative @@ -132,7 +134,7 @@ verbosity 0 UseLEDWaveforms 0 PulseFindingApproach threshold -PulseWindowType Fixed_2023_Gains +PulseWindowType NoDoubleHits DefaultADCThreshold 7 DefaultThresholdType relative @@ -147,7 +149,7 @@ verbosity 0 UseLEDWaveforms 0 PulseFindingApproach threshold -PulseWindowType Fixed_2023_Gains +PulseWindowType NoDoubleHits DefaultADCThreshold 7 DefaultThresholdType relative