KAHIBARO
Creating an Energy Spectrum
Table of Contents
Understanding the Goal
In this step of the gamma ray detector example you already have a simulation that records energy deposition in the scintillation crystal, and you have applied an energy resolution model in the digitizer chain. The goal now is to turn all these individual recorded events into an energy spectrum, that is, a histogram of counts as a function of measured energy.
An energy spectrum is one of the most common outputs in detector simulations. It allows you to see photopeaks, Compton continua, and to compare simulated detector performance with experimental measurements.
From Simulation Output to Spectrum Input
Before you can create a spectrum you must know exactly which quantity you will histogram and in which file it is stored.
In this project the spectrum is usually built from one of two possible outputs.
First, if you are still at the level of raw energy depositions, you can use an energy deposition actor attached to the detector crystal that writes, for each event or step, the deposited energy in the crystal. This is useful for educational purposes but does not yet represent the real detector response.
Second, and more realistic, you can use the singles output from the digitizer. Each single represents a processed detector event that has already gone through energy summation, blurring and thresholds. This is what you typically want for an energy spectrum that should resemble what a real detector would measure.
In both cases the file you read will usually be a ROOT file produced by GATE, containing a tree with branches such as totalEnergy or energy for each recorded interaction or single. The units are typically in MeV, so you must keep your units consistent when you define the spectrum axis.
Always verify:
- Which output object you use for the spectrum, hits, energy deposition actor, or singles.
- The branch name that contains the measured energy, for example
energyortotalEnergy. - The energy unit in the output file. In GATE and Geant4 this is usually MeV, not keV.
Choosing Energy Binning
An energy spectrum is essentially a histogram. You must choose the energy range and the number of bins. These choices affect the appearance of the spectrum and how easy it is to interpret physical structures such as peaks and Compton edges.
Assume the simulated gamma source has a main line at energy $E_{\gamma}$, for example $662 \,\text{keV}$ for Cs-137 or $511 \,\text{keV}$ for PET-like tests. If the simulated energy values are stored in MeV, then $E_{\gamma}$ in MeV is $E_{\gamma,\text{MeV}} = E_{\gamma} / 1000$.
A practical choice is to set the histogram limits slightly below zero and somewhat above the main line, for instance from zero to 1.5 times the photopeak energy. For a 662 keV source that would be from 0 to about $1.0 \,\text{MeV}$.
You also need to decide how many bins you want. More bins give finer detail but require more events to keep the statistical fluctuations acceptable. For basic spectra 512 or 1024 bins are common choices.
You seldom need to compute a bin width directly, because Python tools like NumPy and Matplotlib can take the minimum and maximum energy and the number of bins and compute the bin width internally. However, conceptually the bin width $\Delta E$ is
$$
\Delta E = \frac{E_{\text{max}} - E_{\text{min}}}{N_{\text{bins}}}.
$$
For a clear spectrum, choose:
• An upper energy bound a bit higher than the expected photopeak.
• A number of bins large enough to see the peak shape but not so large that statistical noise hides physical features.
Building the Spectrum with Python
The usual workflow is to read the simulation output with Python and build the spectrum using NumPy. This is where you connect the Monte Carlo data to a simple analysis script.
A typical structure for a spectrum script is as follows. First, load the energy values from the ROOT file using a library such as uproot. Second, convert energies to your preferred unit if needed. Third, compute a histogram. Finally, optionally save the spectrum to a text file or NumPy array for later use.
A simplified example for singles data might look like this, written in a typical Python style for GATE analysis:
import uproot
import numpy as np
# Open the ROOT file and tree
file = uproot.open("output/singles.root")
tree = file["Singles"] # tree name may differ
# Read the energy array (in MeV)
energies_mev = tree["energy"].array(library="np")
# Define histogram parameters
emin_mev = 0.0
emax_mev = 1.0
nbins = 1024
# Create the histogram
hist, bin_edges = np.histogram(
energies_mev,
bins=nbins,
range=(emin_mev, emax_mev)
)
# Compute bin centers for plotting or saving
bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:])After this, the array hist contains the number of counts in each energy bin, and bin_centers gives the corresponding energies.
Even though this looks like a simple step, some subtle choices matter. If your detector model includes an upper energy threshold that cuts events above a certain energy, you should set $E_{\text{max}}$ to that threshold or slightly above it. If there is a lower threshold, you can start the histogram slightly below that value or directly at the threshold.
The histogram must include only events that correspond to valid detector signals. Always apply the same energy thresholds in your analysis as in your digitizer configuration, or you will misinterpret the resulting spectrum.
Plotting and Saving the Spectrum
Once you have the histogram and bin centers, you can visualize the spectrum. For this you typically use Matplotlib. Because the spectrum may span orders of magnitude in count rate, it is often helpful to plot the vertical axis on a logarithmic scale, but for basic introduction a linear scale is usually sufficient.
Here is a simple plotting example:
import matplotlib.pyplot as plt
plt.figure()
plt.step(bin_centers * 1000, hist, where="mid") # convert MeV to keV
plt.xlabel("Energy [keV]")
plt.ylabel("Counts")
plt.title("Simulated Gamma-Ray Energy Spectrum")
plt.grid(True)
plt.tight_layout()
plt.show()This creates a typical step plot where peaks and continua can be clearly seen. For a single photopeak source you expect a prominent peak near the known gamma energy, with a lower energy tail related to partial energy depositions, Compton events that escape, and your blurring model.
To reuse the spectrum later, you can save the data in a simple text or NumPy format. For example:
np.savetxt(
"output/energy_spectrum.txt",
np.column_stack((bin_centers, hist)),
header="Energy[MeV] Counts"
)This allows you to analyze or plot the spectrum again without re-reading the large ROOT file.
Relating the Spectrum to Detector Performance
The final purpose of creating the energy spectrum is usually not just to see a pretty plot, but to evaluate detector characteristics such as energy resolution and detection efficiency.
From the spectrum you can determine the position of the photopeak and its full width at half maximum (FWHM). The relative energy resolution $R$ at the photopeak energy $E_{\text{peak}}$ is defined as
$$
R = \frac{\text{FWHM}}{E_{\text{peak}}}.
$$
This is often expressed as a percentage:
$$
R_{\%} = \frac{\text{FWHM}}{E_{\text{peak}}} \times 100\%.
$$
In the context of this project you already applied an energy blurring model in the digitizer that mimics a given FWHM at a reference energy. When you create the spectrum you can check whether the simulated spectrum indeed shows a peak width compatible with that configuration.
Key use of the energy spectrum:
• Locate the photopeak position and compare it to the known gamma energy.
• Measure the FWHM to verify the implemented energy resolution.
• Inspect the continuum to understand scattering and escape processes in your detector model.
Even without detailed quantitative analysis, simply inspecting the shape and position of the spectrum is a powerful way to confirm that your gamma ray detector simulation behaves in a physically reasonable way and that the previous steps in the project, such as geometry, physics configuration, and digitizer settings, are correctly implemented.
Views: 10