Add Smoothing Detector to Event Formation

trace-to-events/src/channels.rs

@@ -3,13 +3,16 @@ use crate::{
     parameters::{
         AdvancedMuonDetectorParameters, DetectorSettings,
         DifferentialThresholdDiscriminatorParameters, FixedThresholdDiscriminatorParameters, Mode,
-        PeakHeightBasis, Polarity,
+        PeakHeightBasis, Polarity, SmoothingDetectorParameters,
     },
     pulse_detection::{
         AssembleIterable, EventsIterable, Real, WindowIterable,
         advanced_muon_detector::{AdvancedMuonAssembler, AdvancedMuonDetector},
-        detectors::differential_threshold_detector::{
-            DifferentialThresholdDetector, DifferentialThresholdParameters,
+        detectors::{
+            differential_threshold_detector::{
+                DifferentialThresholdDetector, DifferentialThresholdParameters,
+            },
+            smoothing_detector::sec_deriv_smoothing_for_peaks,
         },
         threshold_detector::{ThresholdDetector, ThresholdDuration},
         window::{Baseline, FiniteDifferences, SmoothingWindow},
@@ -51,6 +54,13 @@ pub(crate) fn find_channel_events(
             detector_settings.baseline as Real,
             parameters,
         ),
+        Mode::SmoothingDetector(parameters) => find_smoothing_events(
+            trace,
+            sample_time,
+            detector_settings.polarity,
+            detector_settings.baseline as Real,
+            parameters,
+        ),
     };
     tracing::Span::current().record("num_pulses", result.0.len());
     result
@@ -218,3 +228,40 @@ fn find_advanced_events(
     }
     (time, voltage)
 }
+
+#[tracing::instrument(skip_all, level = "trace")]
+fn find_smoothing_events(
+    trace: &ChannelTrace,
+    sample_time: Real,
+    polarity: &Polarity,
+    baseline: Real,
+    parameters: &SmoothingDetectorParameters,
+) -> (Vec<Time>, Vec<Intensity>) {
+    let sign = match polarity {
+        Polarity::Positive => 1.0,
+        Polarity::Negative => -1.0,
+    };
+    let raw = trace
+        .voltage()
+        .unwrap()
+        .into_iter()
+        .map(|v| sign * (v as Real - baseline))
+        .collect::<Vec<Real>>();
+    let time = (0..raw.len())
+        .map(|t| t as Real * sample_time)
+        .collect::<Vec<Real>>();
+
+    let (time, intensity) = sec_deriv_smoothing_for_peaks(
+        &time,
+        &raw,
+        parameters.noise_centile,
+        parameters.kernel_sigma,
+        parameters.nsig_noise,
+        parameters.min_size,
+    )
+    .unwrap();
+    (
+        time.into_iter().map(|t| t as Time).collect(),
+        intensity.into_iter().map(|v| v as Intensity).collect(),
+    )
+}

trace-to-events/src/parameters.rs

@@ -126,6 +126,23 @@ pub(crate) struct AdvancedMuonDetectorParameters {
     pub(crate) min_amplitude: Option<Real>,
 }
 
+/// Encapsulates the parameters specific to the Smoothing detector.
+#[derive(Default, Debug, Clone, Parser)]
+pub(crate) struct SmoothingDetectorParameters {
+    /// Centile of x to use for noise estimation.
+    #[clap(long)]
+    pub(crate) noise_centile: Real,
+    /// Sigma of the Gaussian kernel for smoothing.
+    #[clap(long)]
+    pub(crate) kernel_sigma: Real,
+    /// Number of standard deviations above noise to use as threshold.
+    #[clap(long)]
+    pub(crate) nsig_noise: Real,
+    /// Minimum size of region to consider a peak, if absent all regions are considered.
+    #[clap(long)]
+    pub(crate) min_size: Option<usize>,
+}
+
 /// Specifies which detector is to be used, and wraps the detector-specific options in each variant.
 #[derive(Subcommand, Debug)]
 pub(crate) enum Mode {
@@ -135,4 +152,6 @@ pub(crate) enum Mode {
     DifferentialThresholdDiscriminator(DifferentialThresholdDiscriminatorParameters),
     /// Detects events using differential discriminators. Event lists consist of time and voltage values.
     AdvancedMuonDetector(AdvancedMuonDetectorParameters),
+    /// Detects events using a smoothed second derivative. Event lists consist of time and voltage values.
+    SmoothingDetector(SmoothingDetectorParameters),
 }

trace-to-events/src/pulse_detection/detectors/mod.rs

@@ -2,6 +2,7 @@
 //! They register detections in the form of a stream of events.
 pub mod advanced_muon_detector;
 pub mod differential_threshold_detector;
+pub mod smoothing_detector;
 pub mod threshold_detector;
 
 use super::{EventData, EventPoint, Pulse, Real, RealArray, TracePoint, pulse::TimeValue};