Efficiency Simulation

Author: Danaé Valdenaire
Created: Jun. 26, 2026
Last updated: Jul. 3rd, 2026 by Philipp Schreiner

---

This tutorial walks you through the efficiency simulation using the cait. But first, a little bit of context is required.

Not all particle events end up in the final spectrum. Some fall below the trigger threshold, others are removed by quality cuts. To correctly interpret the measured spectrum — for example when computing dark matter exclusion limits — the probability of an event surviving this analysis chain must be estimated as a function of energy. This quantity is called the survival probability, or efficiency.

The efficiency is estimated by simulation: artificial pulses of known energy are superimposed onto the real data stream at random timestamps, then processed through the same analysis chain as real data: triggering, main parameter calculation, OF pulse height calculation, template fit, and, finally, quality cuts. Since their true energies are known, the surviving fraction can be extracted as a function of energy.

Through this procedure, we can quantify exactly how many events did, e.g., not trigger because they accidentally landed on top of a testpulse (as any real event could, too), how many triggered pulses were removed through the stability cut because they happend during instable detector periods, how many events were discarded because they were in pileup with another pulse, etc.

See also

For more informations on the algorithm, see dh.efficiency_sim_trigger_of.

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, as well as the EnergyCalibration object from the Energy Calibration notebook. 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)

# Energy calibration
tpas = dh["testpulses/testpulseamplitude", 0]
unique_tpas = np.unique(tpas)
tp_phs = dh["testpulses/templatefit_pars-TL_0.14", 0, :, 0]
tph_bounds = [(0.050, 0.063), (0.110, 0.120), (0.156, 0.167), (0.208, 0.220)]

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))
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)

tp_ts = dh.get_event_iterator("testpulses", channel=0).timestamps[tp_quality_cuts.get_flag()] 
tpas = dh["testpulses/testpulseamplitude", 0][tp_quality_cuts.get_flag()]
tp_phs = dh["testpulses/templatefit_pars-TL_0.14", 0, :, 0][tp_quality_cuts.get_flag()]

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,
)
my_ecal.to_file("tutorial_output/my_ecal_function")

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)
dh.set("events", testpulse_equivalent=TPE, n_channels=2, channel=0, overwrite_existing=True)

# 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)

Simulating timestamps and pulse heights

Our first goal now is to randomly litter the existing data stream that we have with artifical events. For that, we draw random timestamps anywhere on the stream and assign random pulse heights from a spectrum that covers the energy range that we are interested in.

Note

The events are later simulated one at a time, i.e., even if you simulate a million timestamps on a stream that is only one hour long, you will not cause ‘artificial pile up’. The idea here is to answer the statistical question ‘If I placed a single event anywhere on the stream, what is the probability that this event survives my analysis chain (as a function of its pulse height)?’, which is exactly the question about the survival probability/efficiency.

Let’s start with the random timestamps: Here, the only thing that you have to keep in mind is that you don’t simulate too close to the stream’s edges (because if you place a pulse on the last sample of the stream, you cannot meaningfully reconstruct it). This of course means that your efficiency simulation is blind to the very first and last few thousand of samples, but compared to a several hour data stream, this is even less than negligible.

An easy way to do it is to use uniform random numbers like so:

# The number of events we want to simulate.
# The more, the better! Set it to 10^5 for example.
# Here, we set it to 1000 such that the automated
# test that runs over this notebook doesn't take forever.
# The outputs that you see are from running with 10^5
N_sim = 1000 # 10**5

sim_ts = np.sort(
    sp.stats.randint.rvs(
        stream.time[0] + 6*stream.dt_us*record_length, # start a bit later than stream start
        stream.time[-1] - 4*stream.dt_us*record_length, # end a bit earlier than stream end
        size=N_sim,
        random_state=137, # specify seed to make analysis reproducible
    )
)

Important

Always sort your timestamps in ascending order. Reason: The algorithm goes through all simulated events and if it always has to jump back and forth, performance suffers.

Next, we will draw random numbers for our pulse heights. What we actually want, however, is to draw random energies, e.g. between 0 and 20 keV (because that’s where our region of interest for the analysis is). To convert between pulse heights and energies, we remember what we did in the Energy Calibration tutorial. There, we saved the EnergyCalibration object that relates pulse heights to testpulse equivalent energies (TPEs). The TPEs are related to the energies through the CPE factor. So what we have to do is

• First, we define the range of energy (in keV) in which we want to run the simulation. We have to stay in the linear range of the detector otherwise the pulse shape of the particle events and the simulated pulses will differ.

• Then, we convert these energies to TPE energies using our e_cal function. We start with a uniform distribution of energies but (depending on how linear our transfer function is) we will get a non-uniform TPE distribution.

• The last step is to convert TPEs to pulse heights, by dividing by the CPE factor. This value is the ratio between the energy of the peak (in keV) and the position of the peak (in V) in the TPE spectrum.

# Energy range of interest
energy_range = (0, 20)

# Simulate uniform random energies
sim_energies = energy_range[0] + np.diff(energy_range)*sp.stats.uniform.rvs(size=N_sim, random_state=23)

# Load TPE energies and remind ourselves, what the CPE factor was
tpe = dh["events/testpulse_equivalent", 0]
min_fit, max_fit = 0.565, 0.595
fit_x = np.linspace(min_fit, max_fit, 100)
bins = np.linspace(0.1, 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",
)

CPE = 5.89 / gauss_fitpar[0]
print(f"CPE factor: {CPE:.2f} keV/V")

CPE factor: 10.18 keV/V

# Convert energies to TPEs
sim_tpes = sim_energies/CPE

# Load the calibration function from before
my_ecal = vai.EnergyCalibration.from_file("tutorial_output/my_ecal_function")

# Convert TPEs to pulse heights.
# Since the conversion depends on time, we also specify
# that we want the conversion at the simulation timestamps.
sim_phs = my_ecal(sim_ts, sim_tpes)

# Let's see what the simulated pulse height spectrum looks like
vai.Histogram(
    sim_phs,
    bins=500,
    xlabel="Simulated pulse height (V)",
    backend="plotly",
);

Note

After converting TPEs to PHs using the EnergyCalibration object, some of the outputs might be nan. This is because the conversion relies on the existance of regular testpulses and when there are no testpulses (yet), especially in the very beginning and end of a stream, nan values are returned.

Very likely, your experiment’s live time is defined between the first and last testpulse. Hence, the simulated timestamps that you lost are not of interest anyways. Be careful though and check with your experiment’s conventions.

Tip

Since the efficiency usually changes more rapidly at lower energies but is more or less flat at higher energies, it might also be useful to simulate an exponential pulse height spectrum rather than a uniform one.

Trigger efficiency

We are now almost ready to run our efficiency simulation function. But we need two more inputs: The trigger threshold that was used in the main analysis (we didn’t explicitly show you the optimum filter trigger in the tutorials, so we will pick 0.001 V here), and the fit parameters of the pulse template. We fitted the template before and maybe you have the parameters saved either in the dh_SEV or in the /tutorial_output/SEV_fit.xy file. Here, we just quickly run the fit again to get them.

Tip

In principle, you can also use the SEV directly (instead of its fit parameters). The dh.efficiency_sim_trigger_of function accepts it (see the documentation). However, we recommend to always use the fit parameters because with them, the pulses can be extrapolated and nicely placed on the stream and are not just cut off as they would be with just the SEV.

Caution

It is essential that you use the same trigger threshold for the efficiency simulation as you did for triggering your data!

# Load OF and obtain pulse shape model fit parameters
of = vai.OF.from_file("tutorial_output/OF")
sev = vai.SEV.from_file("tutorial_output/SEV")
sev_fitpars, _  = 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)),
)

Now, we are ready to run our efficiency simulation function dh.efficiency_sim_trigger_of. Again, it is a good idea to use the preview=True argument to check if everything looks good before committing to running the function.

sim_kwargs = {
    "stream": stream,
    "trigger_channels": ["Ch0"],
    "testpulse_channels": ["TP0"],
    "of": of,
    "thresholds": [0.001],
    "sim_ts": sim_ts[~np.isnan(sim_phs)],
    "sim_phs": sim_phs[~np.isnan(sim_phs)],
    "sev_fitpars": sev_fitpars,
}

dh.efficiency_sim_trigger_of(**sim_kwargs, preview=True)

What you see above in the preview is the following: What is called event is the stream chunk plus the simulated event. The filtered event is the same after applying the optimum filter. The target index shows you where the event was simulated. When you run the trigger, the trigger indices that are found will be compared against this target index and when there is a trigger within tolerance_samples samples around target index, the event counts as survived trigger. You can change the number of tolerance samples, the size of the stream chunk, as well as the placement of the simualted event on the chunk using the respective arguments of dh.efficiency_sim_trigger_of, but the defaults should serve you well in most cases.

The threshold line shows you the trigger threshold that you put and the search window tells you where the trigger looks for samples above the threshold. Note that the very beginning and ending of the chunk are not searched because of reasons explained in trigger_base.

Finally, the points called triggers are the trigger samples that were found by the trigger algorithm.

Once you convinced yourself that everything works as expected, you can run the simulation:

dh.efficiency_sim_trigger_of(**sim_kwargs)

Two groups have been created in the DataHandler: trig-eff-sim contains trigger information and the simulated chunks, events-eff-sim contains the particle traces which survived the procedure and they are to be considered identical to the events group, i.e. all that happend to the events group, every analysis step, every cut, will have to happen to the events-eff-sim group as well. Comparing the final results to what was simulated will tell you exactly what your efficiency was.

To keep track of where we started from, we also save the simulated energies and TPE values in the respective groups such that we can later access them conveniently.

# Information in the group that contains all simulated events
dh.set(
    "trig-eff-sim",
    simulated_energies=sim_energies[~np.isnan(sim_phs)],
    simulated_tpes=sim_tpes[~np.isnan(sim_phs)],
)
# Information in the group that only contains the survived events
surv_flag = dh["trig-eff-sim/flag_survived_trigger"] * dh["trig-eff-sim/flag_survived_tp"]
dh.set(
    "events-eff-sim",
    simulated_energies=sim_energies[~np.isnan(sim_phs)][surv_flag],
    simulated_tpes=sim_tpes[~np.isnan(sim_phs)][surv_flag],
)

Check out the new groups using

dh.content("*-eff-sim")

events-eff-sim
  event                       (1, 76497, 16384)  (stored externally)
  hours                       (76497,)           float64
  reconstructed_phs           (1, 76497)         float32
  simulated_energies          (76497,)           float32
  simulated_phs               (1, 76497)         float32
  simulated_tpes              (76497,)           float32
  simulated_ts                (76497,)           int32
  time_mus                    (76497,)           int32
  time_s                      (76497,)           int32
  |timestamps                 (76497,)
trig-eff-sim
  event                       (1, 99778, 131072)  (stored externally)
  event_timestamps            (99778,)            int32
  flag_survived_tp            (99778,)            bool
  flag_survived_trigger       (99778,)            bool
  hours                       (99778,)            float64
  reconstructed_phs           (1, 99778)          float32
  simulated_energies          (99778,)            float32
  simulated_phs               (1, 99778)          float32
  simulated_tpes              (99778,)            float32
  time_mus                    (99778,)            int32
  time_s                      (99778,)            int32
  |timestamps                 (99778,)
  trigger_flag                (1, 99778)          bool
  trigger_index               (1, 99778)          int32

You can look at the stream chunks that were simulated (only contains trigger channels) \(\downarrow\)

vai.Preview(dh.get_event_iterator("trig-eff-sim").with_processing(vai.RemoveBaseline()), backend="plotly");

You can look at the events that survived (contains trigger and passive channels) \(\downarrow\)

vai.Preview(dh.get_event_iterator("events-eff-sim").with_processing(vai.RemoveBaseline()), backend="plotly");

Let’s have a look at the survival distribution in terms of simulated amplitudes:

simulated_phs = dh["trig-eff-sim/simulated_phs", 0] # pulse height array we simulated
survived_trigger = dh["trig-eff-sim/flag_survived_trigger"] # events triggered by the algorithm
survived_tp = dh["trig-eff-sim/flag_survived_tp"] # events that are not shadowed by a testpulse

h = vai.Histogram(
    {
        "Simulated": simulated_phs, 
        "Triggered": simulated_phs[survived_trigger],
        "Survived TP": simulated_phs[survived_tp*survived_trigger],
    }, 
    bins=np.linspace(0, 0.1, 100),
    xlabel="Simulated pulse height (V)",
    backend="plotly",
)
h.add_line(x=[0.001, 0.001], y=[0, 1000], name="Trigger threshold")
h.add_line(x=[0.035, 0.035], y=[0, 1000], name="Start Fe line")

We can nicely see that events start to be successfully triggered at the trigger threshold. Higher amplitudes have a larger probability of being triggered. The step that we see around 0.035 V marks the beginning of events originating from the \(^{55}\)Fe calibration source (try plotting dh["events/trigger_phs", 0] as a histogram and compare the amplitudes). Many of those events are present in the data and if another (smaller) pulse comes close to an \(^{55}\)Fe event, it is shadowed because the larger event is preferred by the trigger algorithm. All events larger than \(^{55}\)Fe events will win against them, hence their trigger efficiency will be larger.

Caution

There are two ways how testpulses can influence the trigger efficiency:

• All events that are within half a record window around a testpulse are automatically flagged as testpulse, i.e. the existence of a testpulse (regardless of its height) adds a constant efficiency loss (or dead time).

• By design of the trigger algorithm, a (smaller) pulse can additionally be shadowed by a testpulse up to a full record window before a testpulse. This effect might be especially pronounced for controlpulses.

In this example, the testpulses have little influence on the efficiency curve as their rate is a factor 10 lower than that of the \(^{55}\)Fe source and furthermore, most of them are above 0.1 V anyways (try plotting dh["testpulses/pulse_height", 0]). In your real-life-data, however, they can also cause steps just like the one we see here for \(^{55}\)Fe.

Good news

The calculation of the efficiency will be correct by design: the shadowing of smaller pulses, the dead time due to test and controlpulses, etc. are automatically taken into account since the exact same triggering/event-building algorithm runs on the simulated chunks as it would run on the actual stream. Therefore, you don’t have to worry about any of it.

Nevertheless, it pays off to familiarize yourself with the trigger algorithm and how the efficiency simulation work in detail to prevent confusion. The most common source of confusing is the testpulse shadowing explained in the info box above.

Finally, we are interested in the survival fraction as a function of energy (not amplitude). Conveniently, we have previously saved the energies that we simulated:

sim_energies = dh["trig-eff-sim/simulated_energies"]
surv_energies = dh["trig-eff-sim/simulated_energies"][survived_tp*survived_trigger]

# Bin the energies
bins = np.linspace(0, 20, 200)
bin_centers = 0.5 * (bins[:-1] + bins[1:])

sim_count, bin_edges = np.histogram(sim_energies, bins=bins)
surv_count, _ = np.histogram(surv_energies, bins=bins)
efficiency = surv_count / np.where(sim_count!=0, sim_count, np.nan)

# And plot them (you get extra points if you do this with error bars yourself!)
vai.Scatter(
    x=bin_centers, 
    y=efficiency,
    xlabel="Pulse heights (V)",
    ylabel="Trigger efficiency",
    backend="plotly",
);

What is also interesting is to compare the reconstruced energies to the simulated ones. For that, we have to use our EnergyCalibration object again:

# Downsample factor (otherwise widget lags)
down = 100
vai.ScatterPreview(
    x=dh["events-eff-sim/simulated_energies"][::down],
    y=CPE*my_ecal.inverse(
        dh["events-eff-sim/timestamps"][::down],
        dh["events-eff-sim/reconstructed_phs", 0][::down],
    ),
    ev_it=dh.get_event_iterator("events-eff-sim")[:, ::down],
    xlabel="Simulated energy (keV)",
    ylabel="Reconstructed energy (keV)",
    width=500,
    backend="plotly",
);

As you can see, not all energies are correctly reconstructed and the most likely cause for that is pile up. Even if there’s no pile up, however, small variations due to noise flucutations exist.

Fitting the data with the analytical expression of the efficiency

As a last step, we fit a (modified) error function model to the efficiency curve. Since we have found two ‘plateaus’ in our efficiency, we also fit the two separately. Often, the point where the efficiency reaches 50% of its plateau value is termed the detector threshold. We choose this to be 50% of the first plateau here but note that this notion might differ, depending on context or your experiment’s conventions.

def efficiency_fct(x, p, E_thr, sigma_thr):
    return p * 0.5 * (1 + sp.special.erf((x - E_thr)/np.sqrt(2)/sigma_thr))

popt1, pcov1 = sp.optimize.curve_fit(
    efficiency_fct, 
    bin_centers[bin_centers<5], 
    efficiency[bin_centers<5],
    bounds=([0]*3, [np.inf]*3),
    method="trf",
)
popt2, pcov2 = sp.optimize.curve_fit(
    efficiency_fct, 
    bin_centers[bin_centers>5], 
    efficiency[bin_centers>5],
    bounds=([0]*3, [np.inf]*3),
    method="trf",
)
trigger_eff_fit1 = efficiency_fct(bin_centers, *popt1)
trigger_eff_fit2 = efficiency_fct(bin_centers, *popt2)

scat = vai.Scatter(
    x=bin_centers, 
    y=efficiency,
    xlabel="Simulated energy (keV)",
    ylabel="Trigger efficiency",
    backend="plotly",
)
scat.add_line(
    x=bin_centers,
    y=trigger_eff_fit1,
    name="Fit plateau < Fe",
)
scat.add_line(
    x=bin_centers[bin_centers>5],
    y=trigger_eff_fit2[bin_centers>5],
    name="Fit plateau > Fe",
)
scat.add_line(
    x=[popt1[1]]*2,
    y=[0, 1],
    name=f"Threshold: {popt1[1]:.2f} keV"
)

Cut efficiency

The cut efficiency is similar to the trigger effienciency. But in this case, you want to evaluate how many particle events are removed by your quality cuts. To do so, we will use the group of simulated events, events-eff-sim, as if it was the events group, and apply the exact same analysis chain and quality cuts as for the ‘real’ events.

This means that we have to calculate main parameters, re-applied optimum filter amplitudes, (truncated) template fit parameters, stable detector periods, and quality cuts, just like before:

g = "events-eff-sim"
truncation_limit = 0.14

# THE EXACT SAME PARAMTERS AS FOR THE 'EVENTS' GROUP
trunc_tf_kwargs = {
    "group": g,
    "sev": sev,
    "bl_poly_order": 3,
    "only_channels": 0,
    "truncation_limit": truncation_limit,
    "tag": f"TL_{truncation_limit:.2f}"
}
# THE EXACT SAME PARAMTERS AS FOR THE 'EVENTS' GROUP
of_kwargs = {
    "group": g,
    "of": of,
    "sev": sev,
    "only_channels": 0,
}

dh.cmp(g)
dh.apply_template_fit(**trunc_tf_kwargs)
dh.apply_ofilter(**of_kwargs)

# THE EXACT SAME BOUNDS AS FOR THE 'EVENTS' GROUP
ph_bounds = (0.185, 0.190)
l_bounds = (0.0184, 0.0190)
dh.calc_controlpulse_stability(channel=0, group=g, lb=ph_bounds[0], ub=ph_bounds[1])
dh.calc_controlpulse_stability(channel=1, group=g, lb=l_bounds[0], ub=l_bounds[1])

# THE EXACT SAME CUTS AS FOR THE 'EVENTS' GROUP
quality_cuts = ai.cuts.LogicalCut()
quality_cuts.add_condition(dh[f"{g}/controlpulse_stability", 0])

Caution

Here, you would add all the quality cuts that you previously performed on the ‘events’ group. For simplicity, we only put the testpulse stability here. The objective is to apply exactly the same cuts to the simulated events as to the real events.

Notice that here, we do not mean the quality cuts done when creating the SEV/NPS! Those can be as strict as you have to be to get a clean SEV/NPS. We refer to the quality cuts that reject artefacts in your final energy spectrum which you use to derive your physics results.

The rest is conceptually the same as for the trigger efficiency: We bin the surviving events and compare them to the simulated ones:

surv_energies_cut = surv_energies[quality_cuts.get_flag()]

surv_count_cut, _ = np.histogram(surv_energies_cut, bins=bins)
efficiency_cut = surv_count_cut / np.where(sim_count!=0, sim_count, np.nan)

vai.Line(
    x=bin_centers, 
    y={
        "Survived trigger (and TP)": efficiency,
        "Survived cuts": efficiency_cut,
    },
    xlabel="Simulated energy (keV)",
    ylabel="Efficiency",
    backend="plotly",
);

Tip

Instead of applying all quality cuts at once, you can also add one after the other to visualize the effect of each on the final cut efficiencies.

Feel free to fit the analytic expression from above as an exercise!

Congratulations! You completed your first raw data analysis!

Note on convention

Be careful concerning which efficiency you report. Here, we calculated the survival probability as a function of simulated energies. Double check with your community/experiment whether that’s the physical quantity of interest or if it is, e.g., the survival probability as a function of reconstructed energies, etc.

Note

If you have several channels that are triggered, or if one of them is read out as ‘passive channel’, the procedure is essentially the same as for the one-channel example that we demonstrated above. Just use dh.efficiency_sim_trigger_of with the same arguments that you also used in dh.trigger_of. The events-eff-sim group will then have multiple channels and look just like your ‘real’ events group, such that the full analysis chain of all channels can be applied as before.