HPGe DSP Chains

A DSP chain takes the raw waveform table of one detector channel and returns a table with one row per event and one column per extracted parameter. The HPGe chain dsp_icpc and its variant for compressed raw data dsp_icpc_compressed compute the baseline, the pole-zero (PZ) corrected tail, timing points, energy estimators from several filters, the maximum current (A), the charge-drift parameter (Q-drift), late-charge parameters (LQ), pile-up and saturation flags, and an optional machine-learning quality-cut label. The pulser and auxiliary channels use the shorter chains dsp_puls, dsp_puls_compressed, dsp_aux and dsp_aux_compressed.

The layout of the input tables and the DSPConfig are described in Waveforms and Configuration. The individual filters are described in Filters, and the functions that compute single parameters from waveforms in Waveform Parameters.

Calling a Chain

out = dsp_icpc(data, config, τ, pars_filter; f_evaluate_qc = missing)
out = dsp_icpc_compressed(data, config, τ, pars_filter; f_evaluate_qc = missing)

The input columns are

• dsp_icpc: waveform, baseline, timestamp, eventnumber, daqenergy.
• dsp_icpc_compressed: waveform_presummed, waveform_windowed, presum_rate, baseline, timestamp, eventnumber, daqenergy, t_sat_lo, t_sat_hi, deadtime.

The columns baseline, timestamp, eventnumber and daqenergy are the values reported by the FADC. They are copied to the output as blfc, timestamp, eventID_fadc and e_fc. Encoded columns are decoded first.

Processing Steps

The two variants run the same sequence. In the compressed variant the steps that need the full trace use waveform_presummed (the "presummed waveform") and the steps that need fine time resolution use waveform_windowed (the "windowed waveform"); see Uncompressed and Compressed Input in Waveforms and Configuration.

1. Saturation. The samples equal to the lowest ADC value (0) and to the highest value (2^fc_bit_depth - fc_bit_depth, multiplied by presum_rate for presummed waveforms) are counted, together with the longest run of consecutive saturated samples at each end (saturation).
2. Baseline. A linear fit in bl_window gives mean, standard deviation, slope and offset of the baseline (signalstats). The mean is subtracted from the waveforms.
3. Quality-cut label. If f_evaluate_qc is given, it is evaluated on the baseline-subtracted waveform (see below). Otherwise the label is -1.
4. Extrema. The maximum and minimum of the baseline-subtracted waveform, before the PZ correction.
5. Tail. tailstats fits a straight line to the logarithm of the baseline-subtracted waveform in tail_window. The result is the tail decay time (tail_τ) and the mean and standard deviation of the logarithm. See Decay Time.
6. Pole-zero correction. The waveform is deconvolved with InvCRFilter(τ), which turns the exponential tail into a flat top. Steps 7 to 17 work on the PZ-corrected waveform. The tail window is evaluated again (signalstats) to check the correction.
7. t0. get_t0 applies a fast asymmetric trapezoid with the lengths kwargs_pars.t0_flt_pars and finds the first time the output stays above t0_threshold for at least kwargs_pars.t0_mintot.
8. Fixed-length energy filters. Trapezoidal filters with (rise, flat-top) of 10 µs and 4 µs (e_10410), 5 µs and 3 µs (e_535), 3 µs and 1 µs (e_313) and 10 µs and 10 µs (e_101010). The energy is the maximum of the filter output.
9. Rise-time points. t10, t50, t80, t90, t99 are the times at which the PZ-corrected waveform first stays for kwargs_pars.tx_mintot above 10 %, 50 %, 80 %, 90 % and 99 % of the maximum of the 10 µs/10 µs trapezoid (get_threshold). The drift time is t90 - t0.
10. Maximum of the simple derivative. DerivativeFilter gives the current without smoothing. The time of its maximum is t_amax.
11. Q-drift. get_qdrift integrates the waveform and compares two consecutive integration windows that start at t0, delimited by the first and last element of qdrift_int_length.
12. LQ. For every [early, gap, late] triple of lq_windows and for each of the start points t80, t90 and t_amax, get_lq subtracts the charge in the early window from the charge in the late window.
13. Optimized energy filters. Trapezoidal, CUSP and ZAC filters with the optimized (or default) parameters are applied. The energy is not the filter maximum but the filter output at a fixed pick-off time after t50, read with a polynomial signal estimator (kwargs_pars.sig_interpolation_order and sig_interpolation_length). The pick-off time is t50 + rt + ft/2 for the trapezoid and t50 + flt_length/2 for CUSP and ZAC. The maximum of each filter output and its time are stored as well.
14. Current (A). The maximum in current_window of three versions of the current: the Savitzky-Golay derivative with window sg.wl and degree sg_flt_degree (a_sg), the simple derivative (a_raw), and the simple derivative smoothed by a_mwa_passes alternating passes of a moving average of window mwa.wl (a_mwa). Additional columns are written for the fixed windows a_sg_windows and a_mwa_windows. The maxima are interpolated quadratically around the maximum sample.
15. Pile-up. get_triggers counts the intersections of a smoothed signal with a threshold of inTraceCut_std_threshold times its standard deviation in bl_window, with a minimum time over threshold of kwargs_pars.intrace_mintot. dsp_icpc applies it to the Savitzky-Golay derivative; dsp_icpc_compressed applies it to the presummed waveform after a 2 µs/1 µs trapezoidal filter with a minimum time over threshold of 2 µs.
16. Current rise. t50_current is the time at which the SG derivative of the whole waveform first stays above half of its maximum.
17. Inverted waveform. The PZ-corrected waveform is multiplied by -1. The maxima of the 10 µs/4 µs and 3 µs/1 µs trapezoids (e_10410_inv, e_313_inv) and the t0 of the inverted waveform (t0_inv) are stored; they flag events with a negative-going component after the pulse.

In the compressed variant the pole-zero corrected presummed waveform is additionally evaluated in the two windows auxpz1_window and auxpz2_window.

A Synthetic Example

The following example runs both variants on synthetic waveforms: a linear rise of 1 µs followed by an exponential tail, white noise, and baseline 1000. The configuration is the generic example from Waveforms and Configuration.

using LegendDSP
using PropDicts, Unitful, StableRNGs, TypedTables, RadiationDetectorSignals, ArraysOfArrays

grid(start, stop, step) = PropDict(:start => start, :stop => stop, :step => step)
rt_ft_grid = PropDict(:rt => grid(1.0u"µs", 16.0u"µs", 0.5u"µs"), :ft => grid(1.0u"µs", 4.0u"µs", 0.2u"µs"))

config = DSPConfig(PropDict(
    :enc_pickoff_trap => 40.0u"µs", :enc_pickoff_zac => 41.0u"µs", :enc_pickoff_cusp => 41.0u"µs",
    :flt_length_cusp => 38.0u"µs", :flt_length_zac => 38.0u"µs",
    :bl_window => PropDict(:min => 0.0u"µs", :max => 39.0u"µs"),
    :tail_window => PropDict(:min => 70.0u"µs", :max => 110.0u"µs"),
    :current_window => PropDict(:min => 43.0u"µs", :max => 62.0u"µs"),
    :auxpz1_window => PropDict(:min => 70.0u"µs", :max => 90.0u"µs"),
    :auxpz2_window => PropDict(:min => 90.0u"µs", :max => 110.0u"µs"),
    :t0_threshold => 4.0, :inTraceCut_std_threshold => 5, :sg_flt_degree => 3,
    :qdrift_int_length => [2.5u"µs", 5.0u"µs"],
    :lq_windows => [[2.5, 0.0, 2.5]u"µs"],
    :e_grid_trap => rt_ft_grid, :e_grid_zac => rt_ft_grid, :e_grid_cusp => rt_ft_grid,
    :a_grid_wl_sg => grid(30.0u"ns", 350.0u"ns", 32.0u"ns"),
    :a_mwa_windows => [48.0u"ns", 288.0u"ns"], :a_mwa_passes => 3,
    :a_sg_windows => [60.0u"ns", 100.0u"ns"],
    :flt_defaults => PropDict(
        :sg => 100.0u"ns", :mwa => 96.0u"ns",
        :trap => PropDict(:rt => 5.0u"µs", :ft => 2.5u"µs"),
        :zac => PropDict(:rt => 5.0u"µs", :ft => 2.5u"µs"),
        :cusp => PropDict(:rt => 5.0u"µs", :ft => 2.5u"µs"),
        :dplms => PropDict(:rt => 5.0u"µs", :ft => 2.5u"µs"),
    ),
    :kwargs_pars => PropDict(
        :fc_bit_depth => 16,
        :t0_flt_pars => [40.0u"ns", 100.0u"ns", 2000.0u"ns"],
        :t0_mintot => 1500.0u"ns", :tx_mintot => 32.0u"ns", :intrace_mintot => 100.0u"ns",
        :int_interpolation_order => 3, :int_interpolation_length => 100.0u"ns",
        :sig_interpolation_order => 3, :sig_interpolation_length => 700.0u"ns",
    ),
))

The waveform generator and the raw table of the uncompressed layout:

rng = StableRNG(1)
τ_true = 450.0u"µs"

function synthetic_waveform(amplitude; n = 8192, dt = 16.0u"ns", t_rise = 50.0u"µs", rise = 1.0u"µs")
    t = range(0.0u"ns", step = dt, length = n)
    x = ustrip.(u"µs", t) .- ustrip(u"µs", t_rise)
    shape = map(x) do xi
        xi < 0 ? 0.0 : xi < ustrip(u"µs", rise) ? xi / ustrip(u"µs", rise) :
            exp(-(xi - ustrip(u"µs", rise)) / ustrip(u"µs", τ_true))
    end
    RDWaveform(t, 1000 .+ amplitude .* shape .+ 3 .* randn(rng, n))
end

amplitudes = [4500.0, 5000.0, 5500.0, 6000.0]
N = length(amplitudes)
wvfs = ArrayOfRDWaveforms([synthetic_waveform(a) for a in amplitudes])

raw = Table(
    waveform = wvfs, baseline = fill(1000.0f0, N), timestamp = zeros(UInt64, N),
    eventnumber = UInt32.(1:N), daqenergy = zeros(UInt16, N),
)

Running the chain with an empty pars_filter uses the default filter parameters of the configuration:

out = dsp_icpc(raw, config, τ_true, PropDict())
length(columnnames(out))

62

A selection of the output columns. The three energy estimators are close to the amplitude of the baseline-subtracted signal:

Table(amplitude = amplitudes, e_trap = out.e_trap, e_cusp = out.e_cusp, e_zac = out.e_zac)

Table with 4 columns and 4 rows:
     amplitude  e_trap   e_cusp   e_zac
   ┌─────────────────────────────────────
 1 │ 4500.0     4505.05  4504.99  4505.2
 2 │ 5000.0     5005.83  5005.86  5006.0
 3 │ 5500.0     5505.95  5506.31  5506.14
 4 │ 6000.0     6006.73  6006.72  6006.73

The timing points are in µs, and the drift time is t90 - t0. The tail decay time before the pole-zero correction recovers the decay time that was used to generate the waveforms:

Table(t0 = out.t0, t50 = out.t50, t90 = out.t90, drift_time = out.drift_time,
      tail_τ = uconvert.(u"µs", out.tail_τ))

Table with 5 columns and 4 rows:
     t0          t50         t90         drift_time  tail_τ
   ┌───────────────────────────────────────────────────────────
 1 │ 50.0013 μs  50.501 μs   50.8998 μs  898.513 ns  450.03 μs
 2 │ 50.0011 μs  50.5002 μs  50.9001 μs  899.037 ns  450.23 μs
 3 │ 50.0009 μs  50.5001 μs  50.8995 μs  898.557 ns  450.125 μs
 4 │ 50.0009 μs  50.5005 μs  50.9006 μs  899.729 ns  449.736 μs

Optimized filter parameters override the defaults entry by entry:

pars_filter = PropDict(:trap => PropDict(:rt => 8.0u"µs", :ft => 2.0u"µs"))
out_opt = dsp_icpc(raw, config, τ_true, pars_filter)
(default = out.e_trap, optimized = out_opt.e_trap)

(default = [4505.047960587772, 5005.825048591914, 5505.947687137607, 6006.732146839369], optimized = [4505.19301616319, 5005.780116777349, 5506.237880116372, 6006.82619118547])

The compressed variant needs a presummed and a windowed waveform per event. The following lines build them from the full-rate waveforms by summing groups of four samples and by cutting a window of about 30 µs around the rise:

presum_rate = 4

function presummed(w, r)
    s = reshape(w.signal, r, :)
    RDWaveform(range(first(w.time), step = r * step(w.time), length = size(s, 2)), vec(sum(s, dims = 1)))
end
windowed(w, idx) = RDWaveform(w.time[idx], w.signal[idx])

raw_compressed = Table(
    waveform_presummed = ArrayOfRDWaveforms([presummed(w, presum_rate) for w in wvfs]),
    waveform_windowed = ArrayOfRDWaveforms([windowed(w, 2500:4400) for w in wvfs]),
    presum_rate = fill(UInt16(presum_rate), N),
    baseline = fill(1000.0f0, N), timestamp = zeros(UInt64, N),
    eventnumber = UInt32.(1:N), daqenergy = zeros(UInt16, N),
    t_sat_lo = zeros(UInt16, N), t_sat_hi = zeros(UInt16, N), deadtime = zeros(UInt16, N),
)

out_c = dsp_icpc_compressed(raw_compressed, config, τ_true, PropDict())
Table(amplitude = amplitudes, e_trap = out_c.e_trap, e_trap_over_rate = out_c.e_trap ./ presum_rate,
      t50 = out_c.t50, t50_pre = out_c.t50_pre, trigger_multiplicity = out_c.trigger_multiplicity)

Table with 6 columns and 4 rows:
     amplitude  e_trap   e_trap_over_rate  t50         t50_pre     ⋯
   ┌────────────────────────────────────────────────────────────────
 1 │ 4500.0     18021.2  4505.29           50.5011 μs  50.4766 μs  ⋯
 2 │ 5000.0     20024.4  5006.09           50.5002 μs  50.4764 μs  ⋯
 3 │ 5500.0     22025.0  5506.24           50.5001 μs  50.4764 μs  ⋯
 4 │ 6000.0     24028.2  6007.05           50.5005 μs  50.4766 μs  ⋯

The energies of the compressed variant are in units of presummed ADC counts, presum_rate times larger than the ADC amplitude. The time points t0 to t99 come from the windowed waveform and t50_pre from the presummed waveform; the latter has a coarser time resolution.

Pulser and auxiliary channels are processed by the shorter chains described in Pulser and Auxiliary Chains.

Output Columns

Times are Unitful quantities; their unit is given below. Columns without a unit are in ADC counts (for the compressed variant see the note on presummed values above), counts, or dimensionless. Entries marked "c" exist only in dsp_icpc_compressed; entries marked "u" exist only in dsp_icpc.

FADC Values

Baseline and Tail

A flat PZ-corrected tail has tailslope close to zero. A negative slope indicates that τ is too large, a positive slope that it is too small (see Decay Time).

Extrema

Timing

Energy

For the CUSP and ZAC filters rt is the shaping time σ and ft the length of the flat top; the total length of the filter is flt_length_cusp or flt_length_zac. The decay-time constant of these filters is set to a very long value so that the filters do not apply a CR correction of their own.

Current, Q-drift and LQ

The ratio of a_sg (or a_mwa) to an energy estimator is the A/E parameter used for pulse-shape discrimination; see the A/E Pulse-Shape Parameter tutorial.

Pile-up and Saturation

For the presummed waveforms of the compressed variant the saturation counts are made on the presummed waveform, with the highest value scaled by presum_rate.

Quality-Cut Label

f_evaluate_qc is missing by default and the chain then sets qc_label = -1 for every event without running any classifier. A function can be passed instead. The chain prepares its input in the following way: the baseline-subtracted waveform (before the PZ correction) is smoothed by repeated HaarAveragingFilter passes (5 passes for dsp_icpc, 2 passes for dsp_icpc_compressed; see get_qc_classifier and get_qc_classifier_compressed), normalized to a maximum absolute value of one, and collected in a matrix with one column per event. The function receives this matrix and returns a tuple (labels, decision_values). The labels are stored in qc_label as integers. Machine-Learning Quality Cuts describes how to train such a function with get_qc_ml_func.

# A stand-in for a trained classifier: label every event with 1.
toy_classifier(x) = (ones(size(x, 2)), zeros(size(x, 2)))

out_qc = dsp_icpc(raw, config, τ_true, PropDict(); f_evaluate_qc = toy_classifier)
out_qc.qc_label |> collect

4-element Vector{Int64}:
 1
 1
 1
 1

Pulser and Auxiliary Chains

Besides the HPGe detectors, the LEGEND-200 data acquisition records channels that are not detectors. A test pulser injects charge pulses of a fixed shape and amplitude at regular intervals; its analog output is digitized like a detector channel and shows a step with an exponential tail. Auxiliary channels carry digital or analog signals of the data acquisition and the veto systems. Their waveforms are short spikes on a constant baseline. Examples are the baseline channel, which is read out by forced triggers and records the electronic baseline without a signal, and the muon veto channel. These channels have no pole-zero correction, no filter optimization, and no energy calibration from a gamma line. Their DSP output is used to select events, for example to identify pulser events and forced-trigger events and to check the trigger thresholds and timing.

The production dataflow selects the chain of each channel in the section additional_detectors of the DSP configuration. The section maps a channel name to the name of a LegendDSP.jl function. In the production configuration, one analog pulser channel uses dsp_puls_compressed; the baseline channel, the muon veto channel and a further pulser channel use dsp_aux_compressed. The section raw_fields lists the raw columns that each function reads. The functions take the same DSPConfig as the HPGe chains (see Use in Juleana.jl).

All four functions take a table with the columns baseline, timestamp, eventnumber, daqenergy and the waveform column of the table; dsp_aux_compressed also needs presum_rate. The FADC values baseline, timestamp, eventnumber and daqenergy are copied to the output columns blfc, timestamp, eventID_fadc and e_fc.

Pulser Chain

dsp_puls works on a waveform that has the shape of a detector pulse:

1. The baseline statistics are determined in bl_window with signalstats and the mean is subtracted (blmean, blsigma, blslope, bloffset).
2. e_max is the maximum of the baseline-subtracted waveform.
3. t50 is the first time at which the waveform reaches half of e_max (get_threshold).
4. e_10410 is the maximum of the waveform filtered with TrapezoidalChargeFilter(10u"µs", 4u"µs"). The waveform is not pole-zero corrected, so the trapezoid output depends on the decay time of the pulse.

dsp_puls_compressed runs the same steps on waveform_presummed. The presummed waveform is not rescaled, so the amplitudes are in the units of the presummed samples.

Auxiliary Chain

dsp_aux does not read the configuration. It determines the maximum sample and its time with extremestats. e_max is the maximum minus the baseline reported by the FADC, and t_max is the time of the maximum sample. For dsp_aux_compressed the FADC baseline is multiplied by presum_rate before it is subtracted, because the samples of the presummed waveform are sums of presum_rate samples. t_max keeps the unit of the time axis of the waveform; the pulser chain returns t50 in µs.

A Synthetic Example

The example runs both chains on one synthetic waveform each. The pulser waveform is the charge pulse of the example above. The auxiliary waveform is a spike of 0.3 µs width on a constant baseline. The table of the auxiliary channel has the FADC baseline 100.

puls_raw = Table(
    waveform = ArrayOfRDWaveforms([synthetic_waveform(5000.0)]), baseline = [1000.0f0], timestamp = [UInt64(0)],
    eventnumber = [UInt32(1)], daqenergy = [UInt16(0)],
)
out_puls = dsp_puls(puls_raw, config)

t_aux = range(0.0u"ns", step = 16.0u"ns", length = 1024)
spike(t) = 800.0 * exp(-0.5 * ((ustrip(u"µs", t) - 6.0) / 0.3)^2)
wvf_aux = RDWaveform(t_aux, 100.0 .+ spike.(t_aux) .+ 2.0 .* randn(rng, length(t_aux)))
aux_raw = Table(
    waveform = ArrayOfRDWaveforms([wvf_aux]), baseline = [100.0f0], timestamp = [UInt64(0)],
    eventnumber = [UInt32(1)], daqenergy = [UInt16(0)],
)
out_aux = dsp_aux(aux_raw, config)

(puls = (e_max = only(out_puls.e_max), t50 = only(out_puls.t50), e_10410 = only(out_puls.e_10410)),
 aux = (e_max = only(out_aux.e_max), t_max = only(out_aux.t_max)))

(puls = (e_max = 5002.734662659296, t50 = 50.50078415817402 μs, e_10410 = 4944.864317641659), aux = (e_max = 803.7111780257216, t_max = 6000.0 ns))

The figure marks e_max and t50 of the pulser waveform and e_max and t_max of the auxiliary waveform.

using CairoMakie

wvf_puls = puls_raw.waveform[1]
fig = Figure(size = (700, 600))
ax1 = Axis(fig[1, 1], ylabel = "Signal (ADC counts)", title = "Pulser: dsp_puls")
lines!(ax1, ustrip.(u"µs", wvf_puls.time), wvf_puls.signal)
hlines!(ax1, only(out_puls.blmean) + only(out_puls.e_max), color = :red, linestyle = :dash, label = "blmean + e_max")
vlines!(ax1, ustrip(u"µs", only(out_puls.t50)), color = :green, linestyle = :dash, label = "t50")
xlims!(ax1, 45, 60)
axislegend(ax1, position = :rb)
ax2 = Axis(fig[2, 1], xlabel = "Time (µs)", ylabel = "Signal (ADC counts)", title = "Auxiliary: dsp_aux")
lines!(ax2, ustrip.(u"µs", wvf_aux.time), wvf_aux.signal)
hlines!(ax2, only(aux_raw.baseline) + only(out_aux.e_max), color = :red, linestyle = :dash, label = "baseline + e_max")
vlines!(ax2, ustrip(u"µs", only(out_aux.t_max)), color = :green, linestyle = :dash, label = "t_max")
axislegend(ax2, position = :rt)
fig