Energy Calibration
Author: Danaé Valdenaire
Created: Jun. 26, 2026
Last updated: Jul. 1st, 2026 by Philipp Schreiner
This tutorial walks you through the energy calibration using cait. But first, a little bit of context is required.
The detector response varies over time due to temperature fluctuations and other instabilities. As a result, the pulse height (PH) recorded for a given event does not directly reflect the deposited energy. To account for these variations, testpulses of known amplitude (testpulse amplitude (TPA)) are injected periodically. By tracking their reconstructed amplitude (testpulse pulse height (TPH)) over time, we build a model of how the detector response changes. This model is then used to map each event’s pulse height to a testpulse-equivalent amplitude (TPE) — a corrected quantity that can finally be converted into a physical energy using calibration sources.
Altogether, the detector response can be understood as a (TPA, TPH)-mapping as a function of time. The (TPA, TPH)-mapping is called transfer function and the evolution of a TPH for a fixed TPA over time is called testpulse response. Depending on the situation, one can use different models as presented below:
Available objects for energy calibration
Testpulse Response: Describes the TPH over time
vai.TPRPoly: models the testpulse height over time with a polynomial. Use this when the detector drift is gradual and regular throughout the run (e.g. essentially flat or linearily rising).vai.TPRCubicSpline: models the testpulse height over time with a cubic spline. Use this when the detector drift shows sudden or irregular changes, as the spline can adapt to local variations better than a polynomial. You can choose the ‘smoothness’ of the splines to account for more or less irregularity.
Transfer Function: Describes the relation between TPA and TPH
vai.TFPchip: builds the transfer function using a piecewise cubic interpolation. Use this when the amplitude correction is not perfectly linear across all testpulse amplitudes, as it connects the control points smoothly without introducing unwanted oscillations.vai.TFPoly: builds the transfer function using a polynomial fit. Use this when the amplitude correction is close to linear and a simple, smooth function is sufficient.
Tip
vai.TPRCubicSpline + vai.TFPchip is the most robust combination when the detector response is not perfectly stable.
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, the SEV, NPS, and OF from the Creating SEV, NPS, OF tutorial, and the outputs of Reconstructing the pulse amplitude. 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)
dh.set_filepath(path_h5="tutorial_output", fname="my_first_trigger", appendix=False)
dh.init_empty()
# Trigger
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)
*_, 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_tp = vai.SEV(event_iterator)
*_, rms = vai.apply(vai.TemplateFit(sev_tp, 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)
_, fit_rms = vai.apply(vai.FitBaseline(model=3, where=1.0), dh.get_event_iterator("noise", 0))
lb, ub = np.quantile(fit_rms.flatten(), [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")
# Template fit and OF pulse height
truncation_limit = 0.14
dh.apply_template_fit("events", sev=sev_fit, bl_poly_order=3, only_channels=0, truncation_limit=truncation_limit, tag=f"TL_{truncation_limit:.2f}")
dh.apply_ofilter("events", of=of, sev=sev_fit, only_channels=0)
dh.apply_template_fit("testpulses", sev=sev_tp, bl_poly_order=1, only_channels=0, truncation_limit=truncation_limit, tag=f"TL_{truncation_limit:.2f}")
dh.apply_ofilter("testpulses", of=of_tp, sev=sev_tp, only_channels=0)
# 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)
Cleaning Testpulses
Since the calibration depends crucially on the quality of our testpulse reconstruction, we need to remove corrupted or unstable testpulses before building the testpulse response. We apply several quality cuts on the testpulse parameters. The goal is to keep only testpulses that are well-reconstructed and representative of the detector’s response.
Tip
You can perform a template fit on your testpulses beforehand using the function vai.TemplateFit. This could be helpful for the next step. More details about this in the tutorial Reconstructing the pulse amplitude.
We will apply a stability cut and some quality cuts on our testpulses. The goal at this step is to obtain a clear behavior of each testpulse amplitude distribution over time. Each testpulse amplitude should look like a line (or a band) for the fit to account for the detector response as closely as possible.
dh.content("test*")
# The testpulse amplitudes (the 'energy' injected by the testpulse)
tpas = dh["testpulses/testpulseamplitude", 0]
# The unique testpulsamplitudes (4 in this example)
unique_tpas = np.unique(tpas)
print("Unique test pulse heights: ", unique_tpas)
# The RECONSTRUCTED amplitudes, TPH
tp_phs = dh["testpulses/templatefit_pars-TL_0.14", 0, :, 0]
# Our goal is now to remove everything that has not been nicely
# reconstructed. This can very simply be done by plotting a histogram
# of the reconstructed amplitudes for each TPA individually and to cut
# around them by eye (or with some automated cut based on quantiles)
tph_bounds = [
(0.050, 0.063),
(0.110, 0.120),
(0.156, 0.167),
(0.208, 0.220),
]
# Alternative that cuts at the 10% and 90% quantile
# tph_bounds = [
# np.quantile(tp_phs[tpas==tpa], [0.1, 0.9])
# for tpa in unique_tpas
# ]
vai.Histogram(
{
f"TPA={tpa:.0f}": dh["testpulses/templatefit_pars-TL_0.14", 0, tpas==tpa, 0]
for tpa in unique_tpas
},
bins=(0, 0.4, 300),
xlabel="Testpulse height, TPH (V)",
backend="plotly",
).add_line(
x=sum([[b[0], b[0], None, b[1], b[1], None] for b in tph_bounds], start=[]),
y=[0, 100, None, 0, 100, None]*len(tph_bounds),
name="Keep inside",
);
Unique test pulse heights: [1. 2. 3. 4.]
In addition to the rough cuts above, we will add some more quality cuts. For the quality cuts, as you’ve probably already seen in the Creating SEV, NPS, OF tutorial, there are no general rules. The best way to do it is to look at your data and try what works well. You can apply cuts on all the parameters in your testpulses in your DataHandler. If you performed a template fit on your testpulses, you can also apply an upper limit on the RMS of the fit dh["testpulses/templatefit_rms-TL_0.14", 0]. This is a good way to remove artefacts. Likewise, the RMS that you obtain from the optimum filter pulse height reconstruction works well.
tp_quality_cuts = ai.cuts.LogicalCut()
tp_quality_cuts.add_condition(dh["testpulses/controlpulse_stability", 0])
tp_quality_cuts.add_condition((dh["testpulses/templatefit_rms-TL_0.14", 0] > 0)*(dh["testpulses/templatefit_rms-TL_0.14", 0] < 0.00585))
tp_quality_cuts.add_condition(np.abs(dh["testpulses/baseline_difference", 0]) < 0.01)
tp_quality_cuts.add_condition((dh["testpulses/rise_time", 0] > 0.230)*(dh["testpulses/rise_time", 0] < 0.4))
print(f"Survived: {tp_quality_cuts.counts()}/{tp_quality_cuts.total()}, {100*tp_quality_cuts.counts()/tp_quality_cuts.total():.2f} %")
Survived: 362/576, 62.85 %
To evaluate the quality of your cuts, you could plot the test pulses over time before and after your cuts. With the function vai.ScatterPreview you can click on each point of the scatter plot to visualise each event.
flag = tp_quality_cuts.get_flag()
vai.ScatterPreview(
x=dh["testpulses/hours"][flag],
y=dh["testpulses/templatefit_pars-TL_0.14", 0, :, 0][flag],
ev_it=dh.get_event_iterator("testpulses", channel=0, flag=flag).with_processing(vai.RemoveBaseline()),
xlabel="Time (h)",
ylabel="Testpulse amplitude from template fit (V)",
width=500,
backend="plotly",
);
If you still find a small number of outliers, there is a quick and dirty way to remove them automatically:
cond = np.ones(len(tpas), dtype=bool)
for tpa in unique_tpas:
q0, q1, q2 = np.quantile(tp_phs[tpas==tpa], sp.stats.norm.cdf([-0.5, 0, 0.5]))
cond[tpas==tpa] = np.abs((tp_phs[tpas==tpa] - q1)/(q2 - q0)) < 3
tp_quality_cuts.add_condition(cond)
print(f"Survived: {tp_quality_cuts.counts()}/{tp_quality_cuts.total()}, {100*tp_quality_cuts.counts()/tp_quality_cuts.total():.2f} %")
Survived: 361/576, 62.67 %
When you are happy with all your cleaning, you can save the cut flag and move on to the calibration part \(\downarrow\)
dh.apply_logical_cut(
cut_flag=tp_quality_cuts.get_flag(),
naming="cuts_for_cleaning",
channel=0,
type="testpulses",
delete_old=True,
)
Energy calibration
We will now use the cleaned TPAs, TPHs and the TP timestamps to construct the 3-dimensional response function as discussed in the introduction.
# Load cleaning flag from above
cuts_for_cleaning = dh["testpulses/cuts_for_cleaning", 0]
# Load the testpulse timestamps, TPA and TPH values from before and apply the flag
tp_ts = dh.get_event_iterator("testpulses", channel=0).timestamps[cuts_for_cleaning]
tpas = dh["testpulses/testpulseamplitude", 0][cuts_for_cleaning]
tp_phs = dh["testpulses/templatefit_pars-TL_0.14", 0, :, 0][cuts_for_cleaning]
Now, we are all set to perform the energy calibration with the vai.EnergyCalibration object. At this stage, we will give as inputs the testpulse amplitude over time: tp_ts and tp_phs. We also have to choose the testpulse_response and the transfer_function described in the introduction.
my_ecal = vai.EnergyCalibration(
tp_x=tp_ts,
tp_phs=tp_phs,
tpas=tpas,
testpulse_response=vai.TPRCubicSpline(kernel_length=0.02, remove_outliers=True),
transfer_function=vai.TFPchip(),
max_x_gap=0.5,
)
The first thing you should always check is its preview: It shows the transfer function over time for each TPA on the left and the transfer function for a fixed time on the right. You can click at any time in the left plot to see what the transfer function looks like at that time. Notice that there are some sliders and checkboxes above the plot: They correspond to the transfer function/testpulse response settings that you can specify. Changing them lets you build them in an interactive way to find the combination that work the best to describe your data.
Warning
Changing parameters using the preview sliders/checkboxes does not alter the EnergyCalibration object. For the changes to take effect, you have to input the chosen values in the EnergyCalibration construction in the cell above.
my_ecal.preview(width=450, backend="plotly")
Important
my_ecal is the object that you will use to tranform pulse heights in testpulse equivalent amplitudes and vice-versa.
my_ecal(): pulse height (PH) \(\rightarrow\) testpulse equivalent amplitude (TPE).my_ecal.inverse(): pulse height (PH) \(\leftarrow\) testpulse equivalent amplitude (TPE).
Tip
If you have a doubt on whether you have to use the normal or the inverse, you can print the argument of my_ecal() or my_ecal.inverse() using TAB. One function takes the pulse height as an argument, the other the testpulse amplitude.
The usual route is that you have pulse heights (PHs) and you want the testpulse equivalent amplitudes (TPEs). For that, you need to do this \(\downarrow\)
# Note that here, we load the information about the EVENTS!
# Above, we used the testpulses to BUILD the energy calibration
# object, and now we use it on events.
event_ts = dh.get_event_iterator("events", channel=0).timestamps
event_phs = dh["events/templatefit_pars-TL_0.14", 0, :, 0]
TPE = my_ecal.inverse(event_ts, event_phs)
# Don't forget to save your calculation in your DataHandler.
dh.set("events", testpulse_equivalent=TPE, n_channels=2, channel=0, overwrite_existing=True)
In the other direction, it would be:
phs = my_ecal(event_ts, TPE)
Just like the SEV, NPS, and OF, the EnergyCalibration can also be saved to a file and loaded later (or shared with other analysts) \(\downarrow\)
# Save calibration in a .json file
my_ecal.to_file("tutorial_output/my_ecal_function")
# Load it later like this:
my_ecal = vai.EnergyCalibration.from_file("tutorial_output/my_ecal_function")
CPE factor
CPE factor stands for Converting Pulse Height to Energy. It relates the testpulse equivalent energies that we just obtained from the procedure above to physical energies. The idea is the following: Now that we have removed detector non-linearities using the testpulse calibration, we can plot the TPE spectrum and spot calibration lines of known energies. The CPE factor is the ratio of the energy of the calibration peak (in keV/or eV) over the position of the peak (in V).
\(\text{CPE} = \frac{\mathrm{Peak~energy}}{\mathrm{Peak~position}} \left[\frac{\mathrm{keV}}{\mathrm{V}}\right]\)
Our simulated data includes an \(^{55}\text{Fe}\) source which produces peaks at 5.89 keV and 6.49 keV respectively. We calibrate with the dominant one at 5.89 keV (if you want to be diligent, you’d fit a double peak).
# Get data
tpe = dh["events/testpulse_equivalent", 0]
# ... normally, one would also clean the data to get a nice spectrum
# Fitting the peak with a gaussian to extract the mean
min_fit, max_fit = 0.565, 0.595
fit_x = np.linspace(min_fit, max_fit, 100)
bins = np.linspace(0.5, 1.0, 400)
hist = vai.Histogram(
tpe,
bins=bins,
xlabel="Testpulse equivalent amplitude (V_TP)",
ylabel="Counts",
backend="plotly",
)
gauss_fitpar = sp.stats.norm.fit(tpe[(tpe>min_fit)*(tpe<max_fit)])
hist.add_line(
x=fit_x,
y=len(tpe[(tpe>min_fit)*(tpe<max_fit)])*np.diff(bins)[0]*sp.stats.norm.pdf(fit_x, *gauss_fitpar),
name=f"μ={gauss_fitpar[0]:.3f} V_TP\nσ={gauss_fitpar[1]:.3f} V_TP",
)
Finally, the CPE factor equals to \(\text{CPE} = \frac{\mathrm{Peak~energy}}{\mathrm{Peak~position}} = \frac{5.89}{0.578} = 10.2 \; \text{keV/V}\).
The calibrated spectrum is shown below:
CPE = 5.89 / gauss_fitpar[0]
vai.Histogram(
CPE * tpe,
bins=CPE * bins,
xlabel="Energy (keV)",
ylabel="Counts",
backend="plotly",
);
For more information on the EnergyCalibration and the TransferFunction and TestpulseResponse objects, please refer to their documentations.