{
"cells": [
{
"cell_type": "markdown",
"id": "63072b91-0d65-4c06-bb17-6ebbc07818a5",
"metadata": {
"tags": []
},
"source": [
"# Energy Calibration\n",
"_Author: Danaé Valdenaire_
\n",
"_Created: Jun. 26, 2026_
\n",
"_Last updated: Jul. 1st, 2026 by Philipp Schreiner_\n",
"\n",
"---\n",
"\n",
"This tutorial walks you through the **energy calibration** using `cait`. But first, a little bit of context is required.\n",
"\n",
"The detector response varies over time due to temperature fluctuations and other instabilities. \n",
"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.\n",
"\n",
"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:\n",
"\n",
"```{admonition} Available objects for energy calibration\n",
"**Testpulse Response:** Describes the TPH over time\n",
"\n",
"- [**`vai.TPRPoly`**](cait.versatile.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).\n",
"\n",
"- [**`vai.TPRCubicSpline`**](cait.versatile.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.\n",
"\n",
"**Transfer Function:** Describes the relation between TPA and TPH\n",
"\n",
"- [**`vai.TFPchip`**](cait.versatile.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.\n",
"\n",
"- [**`vai.TFPoly`**](cait.versatile.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.\n",
"```\n",
"\n",
"```{tip}\n",
"`vai.TPRCubicSpline` + `vai.TFPchip` is the most robust combination when the detector response is not perfectly stable. \n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "582ed15f-df99-48b4-870c-d341ab6a6927",
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"import numpy as np\n",
"import scipy as sp\n",
"\n",
"import cait as ai\n",
"import cait.versatile as vai"
]
},
{
"cell_type": "markdown",
"id": "7727ca08-a6d7-4ed8-a996-5653208b579e",
"metadata": {
"tags": []
},
"source": [
"## First, you need (mock) data\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "025f7690-b9c0-4483-85dd-e96f398764da",
"metadata": {},
"outputs": [],
"source": [
"# Generate mockdata\n",
"fdirh5 = \"tutorial_output\"\n",
"hdf5_name = \"my_first_trigger\"\n",
"os.makedirs(fdirh5, exist_ok=True)\n",
"\n",
"record_length = 2**14\n",
"\n",
"trigger_config = {\n",
" \"trigger_channels\": [\"Ch0\"],\n",
" \"passive_channels\": [\"Ch1\"],\n",
" \"testpulse_channels\": [\"TP0\", \"TP1\"],\n",
" \"controlpulses_above\": [9., 9.],\n",
" \"f_noise\": 1000,\n",
" \"copy_events\": True,\n",
"}\n",
"n_channels = len(trigger_config[\"trigger_channels\"]) + len(trigger_config[\"passive_channels\"])\n",
"\n",
"stream = vai.MockStream(seed=137, rate_Hz=2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bc81cc0a-a9fb-49ef-82b7-3878d4730aad",
"metadata": {},
"outputs": [],
"source": [
"# Load DataHandler, trigger, and compute main parameters\n",
"dh = ai.DataHandler(record_length=record_length, \n",
" nmbr_channels=n_channels, \n",
" sample_frequency=stream.sample_frequency)\n",
"dh.set_filepath(path_h5=\"tutorial_output\", fname=\"my_first_trigger\", appendix=False)\n",
"dh.init_empty()\n",
"\n",
"# Trigger\n",
"dh.trigger_zscore(stream, **trigger_config)\n",
"\n",
"for group in [\"events\", \"testpulses\", \"noise\", \"controlpulses\"]:\n",
" dh.cmp(group)\n",
" \n",
"# Stability\n",
"for g in [\"events\", \"testpulses\", \"noise\"]:\n",
" dh.calc_controlpulse_stability(channel=0, group=g, lb=0.185, ub=0.190)\n",
" dh.calc_controlpulse_stability(channel=1, group=g, lb=0.0184, ub=0.0190)\n",
" \n",
"# Create SEV (particle pulses)\n",
"quality_cuts = ai.cuts.LogicalCut()\n",
"quality_cuts.add_condition(dh[\"events/controlpulse_stability\", 0])\n",
"quality_cuts.add_condition((dh[\"events/onset\", 0]>-0.645)*(dh[\"events/onset\", 0]<-0.625))\n",
"quality_cuts.add_condition((dh[\"events/decay_time\", 0]>8)*(dh[\"events/decay_time\", 0]<12))\n",
"quality_cuts.add_condition((dh[\"events/pulse_height\", 0]>0.035) * (dh[\"events/pulse_height\", 0]<0.04))\n",
"\n",
"event_iterator = dh.get_event_iterator(group=\"events\", channel=0, flag=quality_cuts.get_flag())\n",
"sev = vai.SEV(event_iterator)\n",
"*_, rms = vai.apply(vai.TemplateFit(sev, bl_poly_order=1), event_iterator)\n",
"sev = vai.SEV(event_iterator[:, rms<0.005])\n",
"sev.to_file(\"tutorial_output/SEV\")\n",
"pars, _ = sp.optimize.curve_fit(ai.fit.pulse_template, \n",
" sev.t, sev, p0=[-1, -0.5, 1.5, 1, 0.1, 10], \n",
" bounds=((-10, -10, -10, 0, 0, 0), (10, 10, 10, np.inf, np.inf, np.inf)),\n",
")\n",
"sev_fit = vai.SEV(ai.fit.pulse_template(sev.t, *pars), sev.dt_us)\n",
"sev_fit = sev_fit/np.max(sev_fit)\n",
"sev_fit.to_file(\"tutorial_output/SEV_fit\")\n",
"\n",
"# Create SEV (testpulses)\n",
"quality_cuts_TP = ai.cuts.LogicalCut()\n",
"quality_cuts_TP.add_condition(dh[\"testpulses/controlpulse_stability\", 0])\n",
"quality_cuts_TP.add_condition(dh[\"testpulses/testpulseamplitude\", 0]==1)\n",
"\n",
"event_iterator = dh.get_event_iterator(group=\"testpulses\", channel=0, flag=quality_cuts_TP.get_flag())\n",
"sev_tp = vai.SEV(event_iterator)\n",
"*_, rms = vai.apply(vai.TemplateFit(sev_tp, bl_poly_order=1), event_iterator)\n",
"sev_tp = vai.SEV(event_iterator[:, rms<0.00505])\n",
"sev_tp.to_file(\"tutorial_output/SEV_TP\")\n",
"\n",
"# Create NPS\n",
"noise_cuts = ai.cuts.LogicalCut()\n",
"noise_cuts.add_condition(dh[\"noise/controlpulse_stability\", 0])\n",
"noise_cuts.add_condition(abs(dh[\"noise/pulse_height\", 0])<0.004)\n",
"noise_cuts.add_condition(dh[\"noise/variance\", 0]<1e-5)\n",
"\n",
"_, fit_rms = vai.apply(vai.FitBaseline(model=3, where=1.0), dh.get_event_iterator(\"noise\", 0))\n",
"lb, ub = np.quantile(fit_rms.flatten(), [0.1, 0.9])\n",
"noise_cuts.add_condition((fit_rms.flatten() > lb) * (fit_rms.flatten() < ub))\n",
"noise_cleaned = dh.get_event_iterator(\"noise\")[0, noise_cuts.get_flag()]\n",
"nps = vai.NPS(noise_cleaned.with_processing([vai.RemoveBaseline(), vai.TukeyWindow()]))\n",
"nps.to_file(\"tutorial_output/NPS\")\n",
"\n",
"# OF\n",
"of = vai.OF(sev_fit, nps)\n",
"of.to_file(\"tutorial_output/OF\")\n",
"of_tp = vai.OF(sev_tp, nps)\n",
"of_tp.to_file(\"tutorial_output/OF_TP\")\n",
"\n",
"# Template fit and OF pulse height\n",
"truncation_limit = 0.14\n",
"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}\")\n",
"dh.apply_ofilter(\"events\", of=of, sev=sev_fit, only_channels=0)\n",
"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}\")\n",
"dh.apply_ofilter(\"testpulses\", of=of_tp, sev=sev_tp, only_channels=0)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ddafe0f7-6f5b-4316-81b2-302bfe2b34ac",
"metadata": {},
"outputs": [],
"source": [
"# If you have it already, load the existing DataHandler\n",
"dh = ai.DataHandler(\n",
" record_length=record_length, \n",
" nmbr_channels=n_channels, \n",
" sample_frequency=stream.sample_frequency,\n",
")\n",
"dh.set_filepath(path_h5=\"tutorial_output\", fname=\"my_first_trigger\", appendix=False)"
]
},
{
"cell_type": "markdown",
"id": "843c1eb7-2172-4ca1-8b70-84c9e040b9f8",
"metadata": {
"tags": []
},
"source": [
"## Cleaning Testpulses\n",
"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.\n",
"\n",
"```{tip}\n",
"You can perform a template fit on your testpulses beforehand using the function [**`vai.TemplateFit`**](`cait.versatile.TemplateFit`). This could be helpful for the next step. More details about this in the tutorial **Reconstructing the pulse amplitude**.\n",
"```\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "52d1a243-1981-4479-8e70-888d221b19bc",
"metadata": {},
"outputs": [],
"source": [
"dh.content(\"test*\")"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "ffe5a0e0-930f-4663-9bc2-5b680c06e095",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Unique test pulse heights: [1. 2. 3. 4.]\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "3937f5778cb44d5a99ef32bb8c83b8d3",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"VBox(children=(FigureWidget({\n",
" 'data': [{'name': 'TPA=1',\n",
" 'opacity': 0.8,\n",
" 'show…"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# The testpulse amplitudes (the 'energy' injected by the testpulse)\n",
"tpas = dh[\"testpulses/testpulseamplitude\", 0]\n",
"# The unique testpulsamplitudes (4 in this example)\n",
"unique_tpas = np.unique(tpas)\n",
"print(\"Unique test pulse heights: \", unique_tpas)\n",
"\n",
"# The RECONSTRUCTED amplitudes, TPH\n",
"tp_phs = dh[\"testpulses/templatefit_pars-TL_0.14\", 0, :, 0]\n",
"\n",
"# Our goal is now to remove everything that has not been nicely\n",
"# reconstructed. This can very simply be done by plotting a histogram\n",
"# of the reconstructed amplitudes for each TPA individually and to cut\n",
"# around them by eye (or with some automated cut based on quantiles)\n",
"tph_bounds = [\n",
" (0.050, 0.063),\n",
" (0.110, 0.120),\n",
" (0.156, 0.167),\n",
" (0.208, 0.220),\n",
"]\n",
"\n",
"# Alternative that cuts at the 10% and 90% quantile\n",
"# tph_bounds = [\n",
"# np.quantile(tp_phs[tpas==tpa], [0.1, 0.9])\n",
"# for tpa in unique_tpas\n",
"# ]\n",
"\n",
"vai.Histogram(\n",
" {\n",
" f\"TPA={tpa:.0f}\": dh[\"testpulses/templatefit_pars-TL_0.14\", 0, tpas==tpa, 0]\n",
" for tpa in unique_tpas\n",
" },\n",
" bins=(0, 0.4, 300),\n",
" xlabel=\"Testpulse height, TPH (V)\",\n",
" backend=\"plotly\",\n",
").add_line(\n",
" x=sum([[b[0], b[0], None, b[1], b[1], None] for b in tph_bounds], start=[]),\n",
" y=[0, 100, None, 0, 100, None]*len(tph_bounds),\n",
" name=\"Keep inside\",\n",
");"
]
},
{
"cell_type": "markdown",
"id": "520bd508-e8a2-481a-ab40-df5fdec94797",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "776d0601-3eeb-4c14-9ced-dfb67197dbec",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Survived: 362/576, 62.85 %\n"
]
}
],
"source": [
"tp_quality_cuts = ai.cuts.LogicalCut()\n",
"tp_quality_cuts.add_condition(dh[\"testpulses/controlpulse_stability\", 0])\n",
"tp_quality_cuts.add_condition((dh[\"testpulses/templatefit_rms-TL_0.14\", 0] > 0)*(dh[\"testpulses/templatefit_rms-TL_0.14\", 0] < 0.00585))\n",
"tp_quality_cuts.add_condition(np.abs(dh[\"testpulses/baseline_difference\", 0]) < 0.01)\n",
"tp_quality_cuts.add_condition((dh[\"testpulses/rise_time\", 0] > 0.230)*(dh[\"testpulses/rise_time\", 0] < 0.4))\n",
"\n",
"print(f\"Survived: {tp_quality_cuts.counts()}/{tp_quality_cuts.total()}, {100*tp_quality_cuts.counts()/tp_quality_cuts.total():.2f} %\")"
]
},
{
"cell_type": "markdown",
"id": "8fdc01a2-a735-44de-859b-22844e8095e0",
"metadata": {},
"source": [
"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`**](`cait.versatile.ScatterPreview`) you can click on each point of the scatter plot to visualise each event. "
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "5d0de30d-9c60-4e7c-bdd3-31833a2c22db",
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "325e064ca046430bb5796d9ad549df5b",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"VBox(children=(HBox(children=(Button(description='cut selected', style=ButtonStyle(), tooltip='cut selected'),…"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"flag = tp_quality_cuts.get_flag()\n",
"\n",
"vai.ScatterPreview(\n",
" x=dh[\"testpulses/hours\"][flag], \n",
" y=dh[\"testpulses/templatefit_pars-TL_0.14\", 0, :, 0][flag],\n",
" ev_it=dh.get_event_iterator(\"testpulses\", channel=0, flag=flag).with_processing(vai.RemoveBaseline()),\n",
" xlabel=\"Time (h)\",\n",
" ylabel=\"Testpulse amplitude from template fit (V)\",\n",
" width=500,\n",
" backend=\"plotly\",\n",
");"
]
},
{
"cell_type": "markdown",
"id": "823b7aaa-ae82-482a-8aa1-fd636a8628f7",
"metadata": {},
"source": [
"If you still find a small number of outliers, there is a quick and dirty way to remove them automatically:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "3e4f0714-0878-4c35-bb01-2e6e31d9e44c",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Survived: 361/576, 62.67 %\n"
]
}
],
"source": [
"cond = np.ones(len(tpas), dtype=bool)\n",
"for tpa in unique_tpas:\n",
" q0, q1, q2 = np.quantile(tp_phs[tpas==tpa], sp.stats.norm.cdf([-0.5, 0, 0.5]))\n",
" cond[tpas==tpa] = np.abs((tp_phs[tpas==tpa] - q1)/(q2 - q0)) < 3\n",
"\n",
"tp_quality_cuts.add_condition(cond)\n",
"print(f\"Survived: {tp_quality_cuts.counts()}/{tp_quality_cuts.total()}, {100*tp_quality_cuts.counts()/tp_quality_cuts.total():.2f} %\")"
]
},
{
"cell_type": "markdown",
"id": "7f4102ca-d754-4d9d-b3c4-e975b0051c36",
"metadata": {},
"source": [
"When you are happy with all your cleaning, you can save the cut flag and move on to the calibration part $\\downarrow$"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "20493417-eb66-43ff-9910-3153501bdc8f",
"metadata": {},
"outputs": [],
"source": [
"dh.apply_logical_cut(\n",
" cut_flag=tp_quality_cuts.get_flag(), \n",
" naming=\"cuts_for_cleaning\",\n",
" channel=0,\n",
" type=\"testpulses\",\n",
" delete_old=True,\n",
")"
]
},
{
"cell_type": "markdown",
"id": "ff4b1c3b-49d6-4e01-8b3b-9341a5dcea5f",
"metadata": {
"tags": []
},
"source": [
"## Energy calibration\n",
"We will now use the cleaned **TPAs**, **TPHs** and the **TP timestamps** to construct the 3-dimensional response function as discussed in the introduction."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "560b0dc5-51d9-4cf7-9246-d0fb7541387e",
"metadata": {},
"outputs": [],
"source": [
"# Load cleaning flag from above\n",
"cuts_for_cleaning = dh[\"testpulses/cuts_for_cleaning\", 0]\n",
"\n",
"# Load the testpulse timestamps, TPA and TPH values from before and apply the flag\n",
"tp_ts = dh.get_event_iterator(\"testpulses\", channel=0).timestamps[cuts_for_cleaning] \n",
"tpas = dh[\"testpulses/testpulseamplitude\", 0][cuts_for_cleaning]\n",
"tp_phs = dh[\"testpulses/templatefit_pars-TL_0.14\", 0, :, 0][cuts_for_cleaning]"
]
},
{
"cell_type": "markdown",
"id": "48a1fc57-59df-4d96-b026-4cfd28f99cd3",
"metadata": {},
"source": [
"Now, we are all set to perform the energy calibration with the [**`vai.EnergyCalibration`**](cait.versatile.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."
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "0a998b9e-cd9f-4e92-9b00-dcf96faa40f9",
"metadata": {},
"outputs": [],
"source": [
"my_ecal = vai.EnergyCalibration(\n",
" tp_x=tp_ts,\n",
" tp_phs=tp_phs,\n",
" tpas=tpas,\n",
" testpulse_response=vai.TPRCubicSpline(kernel_length=0.02, remove_outliers=True),\n",
" transfer_function=vai.TFPchip(),\n",
" max_x_gap=0.5,\n",
")"
]
},
{
"cell_type": "markdown",
"id": "3d87185e-78a3-4449-9810-f5d4433689d4",
"metadata": {},
"source": [
"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.\n",
"\n",
"```{warning}\n",
"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.\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "c00bdf35-c698-46dc-bcab-b33306da564d",
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "5ad1a896fd6e4ea9bc2dd683f11786d5",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"HBox(children=(VBox(children=(HBox(children=(FloatSlider(value=0.02, description='kernel_length', max=3.0, ste…"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"my_ecal.preview(width=450, backend=\"plotly\")"
]
},
{
"cell_type": "markdown",
"id": "0d24e8ff-b740-46c4-ba5a-04a5db16d8a0",
"metadata": {},
"source": [
"```{important}\n",
"`my_ecal` is the object that you will use to tranform pulse heights in testpulse equivalent amplitudes and vice-versa. \n",
"- `my_ecal()`: pulse height (PH) $\\rightarrow$ testpulse equivalent amplitude (TPE).\n",
"- `my_ecal.inverse()`: pulse height (PH) $\\leftarrow$ testpulse equivalent amplitude (TPE).\n",
"```\n",
"\n",
"```{tip} \n",
"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.\n",
"```\n",
"\n",
"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$"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "200d71f2-d50b-4c96-9210-7698c2922dc2",
"metadata": {},
"outputs": [],
"source": [
"# Note that here, we load the information about the EVENTS!\n",
"# Above, we used the testpulses to BUILD the energy calibration\n",
"# object, and now we use it on events.\n",
"event_ts = dh.get_event_iterator(\"events\", channel=0).timestamps\n",
"event_phs = dh[\"events/templatefit_pars-TL_0.14\", 0, :, 0]\n",
"\n",
"TPE = my_ecal.inverse(event_ts, event_phs)\n",
"\n",
"# Don't forget to save your calculation in your DataHandler.\n",
"dh.set(\"events\", testpulse_equivalent=TPE, n_channels=2, channel=0, overwrite_existing=True)"
]
},
{
"cell_type": "markdown",
"id": "f69cac3a-fa5d-4206-8d63-d7b97038aca8",
"metadata": {},
"source": [
"In the other direction, it would be:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "6427817f-8250-4605-99b7-dd52780328fd",
"metadata": {},
"outputs": [],
"source": [
"phs = my_ecal(event_ts, TPE)"
]
},
{
"cell_type": "markdown",
"id": "f10863a4-6b10-4962-bd6b-9f7278563304",
"metadata": {},
"source": [
"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$"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "e9a73fa9-7c1d-41b4-85f0-a20c542cb2f1",
"metadata": {},
"outputs": [],
"source": [
"# Save calibration in a .json file\n",
"my_ecal.to_file(\"tutorial_output/my_ecal_function\")\n",
"\n",
"# Load it later like this:\n",
"my_ecal = vai.EnergyCalibration.from_file(\"tutorial_output/my_ecal_function\")"
]
},
{
"cell_type": "markdown",
"id": "b9a1076a-dc1b-4bf9-b68e-09b9304ff78f",
"metadata": {
"tags": []
},
"source": [
"## CPE factor\n",
"**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). \n",
"\n",
"$\\text{CPE} = \\frac{\\mathrm{Peak~energy}}{\\mathrm{Peak~position}} \\left[\\frac{\\mathrm{keV}}{\\mathrm{V}}\\right]$\n",
"\n",
"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)."
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "cecb9567-1946-49f2-bfcf-88c8727a45c6",
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "762d9f4dc85549f5b012527a21c4ad92",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"VBox(children=(FigureWidget({\n",
" 'data': [{'showlegend': False,\n",
" 'type': 'histogram',\n",
" …"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Get data\n",
"tpe = dh[\"events/testpulse_equivalent\", 0]\n",
"\n",
"# ... normally, one would also clean the data to get a nice spectrum\n",
"\n",
"# Fitting the peak with a gaussian to extract the mean\n",
"min_fit, max_fit = 0.565, 0.595\n",
"fit_x = np.linspace(min_fit, max_fit, 100)\n",
"bins = np.linspace(0.5, 1.0, 400)\n",
"\n",
"hist = vai.Histogram(\n",
" tpe,\n",
" bins=bins,\n",
" xlabel=\"Testpulse equivalent amplitude (V_TP)\", \n",
" ylabel=\"Counts\",\n",
" backend=\"plotly\",\n",
")\n",
"\n",
"gauss_fitpar = sp.stats.norm.fit(tpe[(tpe>min_fit)*(tpemin_fit)*(tpe