Reconstructing the pulse amplitude

Author: Àfrica González Pedraza, Danaé Valdenaire
Created: May 28th, 2026
Last updated: June 28th, 2026 by Philipp Schreiner

---

Previously, we have estimated the amplitudes of our pulses by the pulse_height main parameter. It is very easy to calculate and doesn’t need any special input. It is calculated by applying a moving average to your pulses and searching its maximum. Nothing fancy and perfectly useful to get a first idea of the energy spectrum.

Here, we learn two more sophisticated ways: This tutorial teaches you how to reconstruct the pulse amplitudes employing:

• Standard Event Fit (SEV fit), also called Template Fit

• Optimum Filter (OF)

While both the OF and the truncated SEV fit methods provide similar, reliable event amplitude reconstruction up to the truncation limit (TL), they have distinct advantages:

• The SEV fit allows you to reconstruct amplitudes even when an event pushes the detector out of its linear range, i.e. above the so-called truncation limit (TL). This is impossible with the OF.

• The OF enhances the signal-to-noise ratio. Therefore, it outperforms the SEV fit method when precisely reconstructing sub-keV/low-energy events, where the signal-to-noise ratio is low.

Usually, both methods are applied and used in their respective regimes.

import os
import numpy as np
import scipy as sp
import cait as ai
import cait.versatile as vai

First, you need (mock) data

Make sure that you have the DataHandler that you obtained after triggering and the SEV, NPS, and OF from the Creating SEV, NPS, OF tutorial. If you don’t have them or want to start fresh, run the following cells (after deleting the old HDF5 file if you want to restart from the beginning). Otherwise, you can skip this step.

# Generate mockdata
fdirh5 = "tutorial_output"
hdf5_name = "my_first_trigger"
os.makedirs(fdirh5, exist_ok=True)

record_length = 2**14

trigger_config = {
    "trigger_channels": ["Ch0"],
    "passive_channels": ["Ch1"],
    "testpulse_channels": ["TP0", "TP1"],
    "controlpulses_above": [9., 9.],
    "f_noise": 1000,
    "copy_events": True,
}
n_channels = len(trigger_config["trigger_channels"]) + len(trigger_config["passive_channels"])

stream = vai.MockStream(seed=137, rate_Hz=2)

# Load DataHandler, trigger, and compute main parameters
dh = ai.DataHandler(record_length=record_length, 
                    nmbr_channels=n_channels, 
                    sample_frequency=stream.sample_frequency)

# Set the path to the desired HDF5 file
dh.set_filepath(path_h5="tutorial_output", fname="my_first_trigger", appendix=False)
dh.init_empty()

dh.trigger_zscore(stream, **trigger_config)

for group in ["events", "testpulses", "noise", "controlpulses"]:
    dh.cmp(group)
    
# Stability
for g in ["events", "testpulses", "noise"]:
    dh.calc_controlpulse_stability(channel=0, group=g, lb=0.185, ub=0.190)
    dh.calc_controlpulse_stability(channel=1, group=g, lb=0.0184, ub=0.0190)
    
# Create SEV (particle pulses)
quality_cuts = ai.cuts.LogicalCut()
quality_cuts.add_condition(dh["events/controlpulse_stability", 0])
quality_cuts.add_condition((dh["events/onset", 0]>-0.645)*(dh["events/onset", 0]<-0.625))
quality_cuts.add_condition((dh["events/decay_time", 0]>8)*(dh["events/decay_time", 0]<12))
quality_cuts.add_condition((dh["events/pulse_height", 0]>0.035) * (dh["events/pulse_height", 0]<0.04))

event_iterator = dh.get_event_iterator(group="events", channel=0, flag=quality_cuts.get_flag())
sev = vai.SEV(event_iterator)
fitpar, _, rms = vai.apply(vai.TemplateFit(sev, bl_poly_order=1), event_iterator)
sev = vai.SEV(event_iterator[:, rms<0.005])
sev.to_file("tutorial_output/SEV")
pars, _  = sp.optimize.curve_fit(ai.fit.pulse_template, 
    sev.t, sev, p0=[-1, -0.5, 1.5, 1, 0.1, 10], 
    bounds=((-10, -10, -10, 0, 0, 0), (10, 10, 10, np.inf, np.inf, np.inf)),
)
sev_fit = vai.SEV(ai.fit.pulse_template(sev.t, *pars), sev.dt_us)
sev_fit = sev_fit/np.max(sev_fit)
sev_fit.to_file("tutorial_output/SEV_fit")

# Create SEV (testpulses)
quality_cuts_TP = ai.cuts.LogicalCut()
quality_cuts_TP.add_condition(dh["testpulses/controlpulse_stability", 0])
quality_cuts_TP.add_condition(dh["testpulses/testpulseamplitude", 0]==1)

event_iterator = dh.get_event_iterator(group="testpulses", channel=0, flag=quality_cuts_TP.get_flag())
sev = vai.SEV(event_iterator)
fitpar, _, rms = vai.apply(vai.TemplateFit(sev, bl_poly_order=1), event_iterator)
sev_tp = vai.SEV(event_iterator[:, rms<0.00505])
sev_tp.to_file("tutorial_output/SEV_TP")

# Create NPS
noise_cuts = ai.cuts.LogicalCut()
noise_cuts.add_condition(dh["noise/controlpulse_stability", 0])
noise_cuts.add_condition(abs(dh["noise/pulse_height", 0])<0.004)
noise_cuts.add_condition(dh["noise/variance", 0]<1e-5)

noise_events = dh.get_event_iterator("noise")[0, noise_cuts.get_flag()].with_processing(vai.RemoveBaseline())
_, fit_rms = vai.apply(vai.FitBaseline(model=3, where=1.0), dh.get_event_iterator("noise", 0))
lb, ub = np.quantile(fit_rms, [0.1, 0.9])
noise_cuts.add_condition((fit_rms.flatten() > lb) * (fit_rms.flatten() < ub))
noise_cleaned = dh.get_event_iterator("noise")[0, noise_cuts.get_flag()]
nps = vai.NPS(noise_cleaned.with_processing([vai.RemoveBaseline(), vai.TukeyWindow()]))
nps.to_file("tutorial_output/NPS")

# OF
of = vai.OF(sev_fit, nps)
of.to_file("tutorial_output/OF")
of_tp = vai.OF(sev_tp, nps)
of_tp.to_file("tutorial_output/OF_TP")

# If you have it already, load the existing DataHandler
dh = ai.DataHandler(
    record_length=record_length, 
    nmbr_channels=n_channels, 
    sample_frequency=stream.sample_frequency,
)
dh.set_filepath(path_h5="tutorial_output", fname="my_first_trigger", appendix=False)

Standard event fit (template fit)

In the following, the (truncated) SEV fit method is used to reconstruct the energy of the pulses. It consists of fitting the SEV to each event, where the resulting amplitude is proportional to the energy deposition. The fit can also account for a polynomial baseline and a slight time offset between template and the pulse, which are fitted simultaneously with the scaled SEV.

When the detector exits the linear range of the phase transition, the pulse saturates. In this case, only samples below a truncation limit (TL) are included in the fit.

The procedure is as follows:

1. Approximate the TL: Fit the SEV without truncation and plot the fit RMS versus the SEV amplitude. The TL corresponds to the amplitude at which the RMS starts increasing, marking the onset of saturation.

2. Apply the truncated SEV fit to all events surviving the stability cuts.

3. Quality cuts: Remove events with a high RMS or an inconsistent amplitude. Ideally, fit the RMS distribution with a Gaussian or Rayleigh function and cut within a chosen number of standard deviations around the mean.

Events surviving these cuts are used for the energy spectrum reconstruction.

Template fit

We start by fitting the SEV to the entire event spectrum without truncation to spot where the detector enters the non-linear regime. The truncation limit is our benchmark to distinguish between the linear and the non-linear regime.

The easiest and most convenient way to do it is through the dedicated DataHandler function. You supply it with the SEV, the group (events, testpulses, …) that you want to act on, the SEV you’d like to use, the baseline polynomial degree (0 for constant baseline, etc.), as well as many other things for which we refer you to the documentation of dh.apply_template_fit.

One of the most useful features is that you can set preview=True when calling the function. If you do, an interactive plot opens where you can see the effects of the function on the events without calculating and storing anything (yet) in the DataHandler. You can use it if you are not certain on what the input parameters should be.

# Load the SEV we created previously
sev = vai.SEV.from_file("tutorial_output/SEV")

# You can put them into the function directly but
# I think it's convenient to have them in a dictionary
tf_kwargs = {
    "group": "events",
    "sev": sev,
    "bl_poly_order": 3,
    "only_channels": 0,
}

# Don't run the fit yet, just show what it would do
dh.apply_template_fit(**tf_kwargs, preview=True)

If you think it’s doing the right thing, you run it without the preview argument:

dh.apply_template_fit(**tf_kwargs)

The output are new datasets templatefit_pars, templatefit_rms, and templatefit_shift which are detailed in vai.TemplateFit.

Next, visualize the SEV fit amplitude vs. the RMS. While scatter plots are faster to render, they contain less information; therefore, density plots are recommended. It is important to zoom in on the different regions. Finally, select the TL below the amplitude where the RMS starts to increase.

fitted_amplitude = dh["events/templatefit_pars", 0, :, 0]
fit_rms = dh["events/templatefit_rms", 0]

# Always make sure to only consider POSITIVE RMS VALUES
# because when the fit fails for some events, the respective
# RMS values are set to -404.
plot_flag = fit_rms > 0

vai.ScatterPreview(
    x=fitted_amplitude[plot_flag], 
    y=fit_rms[plot_flag],
    ev_it=dh.get_event_iterator("events", 0, flag=plot_flag),
    xlabel="Template fit amplitude (V)",
    ylabel="Template fit rms (V)",
    width=500,
    backend="plotly",
).scatter.add_line(x=[0.14, 0.14], y=[0.004, 0.01], name="Truncation limit");

By clicking on the events to preview what they look like, we learn that the band at 0.0075 V corresponds to pileup events. The band around 0.005 corresponds to ‘good’ events. For those, we can see that the RMS starts to ‘take off’ around 0.15 V. We decide to put our truncation limit to 0.14 V.

Plotting with Viztool

datasets = {
    "Time (h)": ["hours", None, None],
    "Pulse Height Phonon (V)": ["pulse_height", 0, None],
    "Amplitude (V) ": ["templatefit_pars", 0, 0, None],
    "Shift ": ["templatefit_shift", 0, None],
    "RMS ": ["templatefit_rms", 0,  None],
}

ai.VizTool(
    datahandler=dh,
    group="events",
    datasets=datasets,
    bins=100,
).show()

If you don’t have any saturated pulses, you can stop here. Otherwise, you’ll no do the truncated fit.

Template fit with truncation limit

Now, you can run the template fit by adding the truncation limit. The fit will reconstruct all the saturated pulses so you can access the energies above the saturation. Don’t forget to verify if the fit is done properly using the preview functionality. Click through the events: The ones below the TL will not be affected, but for ones above the TL, you’ll see a difference.

truncation_limit = 0.14

trunc_tf_kwargs = {
    **tf_kwargs,
    "truncation_limit": truncation_limit,
    "tag": f"TL_{truncation_limit:.2f}"
}

dh.apply_template_fit(**trunc_tf_kwargs, preview=True)

Again, once you’re happy with what it’s doing, you can run the fit for all events:

dh.apply_template_fit(**trunc_tf_kwargs)

Note

If you wish to re-do the fit, you must delete the previous fit before re-applying the fit with the new parameters. You can do this by running: `dh.drop(“events”, “templatefit_pars_TL_xxx”)

But of course, you can also just keep the old result and choose a new ‘tag’.

For more information on the template fit function, especially to see how it works with two channels simultaneously (correlated fit), check out dh.apply_template_fit, vai.TemplateFit and vai.TemplateFitCorrelated.

Reconstructing the amplitude with the optimum filter (OF)

The energy reconstruction in the linear range of the detector, i.e. below the TL, can be best performed by filtering the recorded traces with the optimum filter (OF). The OF is a frequency filter that considers both the noise power spectrum (NPS) and the standard event (SEV), enhancing the characteristic frequencies present in the signal while suppressing those dominant in the noise. While the template fit is powerful to access the higher energies, the OF can reconstruct efficiently the low energies. In a sense, those two methods are complementary.

of = vai.OF.from_file("tutorial_output/OF")
sev = vai.SEV.from_file("tutorial_output/SEV_fit")
nps = vai.NPS.from_file("tutorial_output/NPS")

of_kwargs = {
    "group": "events",
    "of": of,
    "sev": sev,
    "only_channels": 0,
}

dh.apply_ofilter(**of_kwargs, preview=True)

# Output: of_ph, of_max_val, of_eval_pos, of_max_pos, of_rms, of_peak_rms

<cait.versatile.plot.highlevel.preview.Preview at 0x282763950>

dh.apply_ofilter(**of_kwargs)

We can also apply the OF on testpulses (we will need it later in the Energy calibration tutorial).

Warning

Be careful to use the OF with the testpulse template and not the event template. If the template you used to create the OF doesn’t match the pulse shape that you’re filtering, weird things will happen.

sev_tp = vai.SEV.from_file("tutorial_output/SEV_TP")
of_tp = vai.OF.from_file("tutorial_output/OF_TP")

tf_tp_kwargs = {
    "group": "testpulses",
    "of": of_tp,
    "sev": sev_tp,
    "only_channels": 0,
}

of_tp_kwargs = {
    "group": "testpulses",
    "sev": sev_tp,
    "bl_poly_order": 1,
    "only_channels": 0,
    "truncation_limit": truncation_limit,
    "tag": f"TL_{truncation_limit:.2f}",
}

dh.apply_ofilter(**tf_tp_kwargs)
dh.apply_template_fit(**of_tp_kwargs)

With VizTool

Inspect the RMS vs. OF amplitude plot alongside the other parameters. Here, the RMS is computed between the filtered SEV and the filtered event traces, while the peak RMS is calculated similarly but restricted to data points around the peak position. When you are satisfied, apply cuts on the OF amplitude (using the TL as an upper bound) and on the RMS or peak RMS to keep only well-reconstructed events.

datasets = {
    "Time (h)": ['hours', None, None],
    "Pulse Height Phonon (V)": ["pulse_height", 0, None],
    "OF Amplitude (V) ": ["of_ph", 0, None],
    "OF Max Val ": ["of_max_val", 0, None],
    "OF Eval Pos ": ["of_eval_pos", 0, None],
    "OF Max Pos ": ["of_max_pos", 0, None],
    "OF RMS ": ["of_rms", 0, None],
    "OF Peak RMS ": ["of_peak_rms", 0, None],
}

viz = ai.VizTool(
    datahandler=dh,
    group="events",
    datasets=datasets,
    bins=100,
)
viz.show()

For more information on the OF pulse height reconstruction, especially to see how it works with two channels simultaneously (correlated evaluation), check out dh.apply_ofilter, and vai.OFPulseHeight.

Comparing the two

Finally, we plot the results of the OF reconstruction with the one from the template fit. We can clearly see that the amplitudes agree below the truncation limit but that the OF reconstruction fails above. At the same time, the OF reconstruction has less spread (higher resolution) at low energies.

clean_flag = (dh["events/templatefit_rms", 0]<0.0053) * (dh["events/templatefit_rms", 0]>0) * (dh["events/of_rms", 0]>0)

vai.Scatter(
    x=dh["events/templatefit_pars-TL_0.14", 0, :, 0][clean_flag],
    y=dh["events/of_ph", 0][clean_flag],
    xlabel="Truncated template fit amplitude (V)",
    ylabel="Optimum filter amplitude (V)",
    xrange=(0, 0.25),
    backend="plotly",
).add_line(x=[0, 0.2], y=[0, 0.2])

We see that both methods agree below the trucation limit, but as we approach it, the OF starts to underestimate the pulse height.