Simulating the Signal of a Gamma Detector¶

This IPython Notebook explains how to calculate the signal of a gamma detector using the pulse-height tally. It also shows how to create a realistic image of a detector simulation that takes into account the detector's intrinsic resolution.

In particular, this Notebook illustrates the following features:

• setting up a simulation of a gamma detector
• defining a gamma source from a specific isotopic mixture
• reprocessing the simulation output to obtain a realistic pulse-height spectrum

Setting up the Detector Model¶

In the following, we will set up a simplified model of a gamma scintillation detector. In front of this gamma detector, we will place two point sources:

• The first source will launch $10^{12}$ gamma particles with an energy of 200 keV per second toward the detector.
• The second source simulates the emissions of 1 gram of Cs137.

Cs137 is a radioactive isotope that primarily emits gamma particles with an energy of 661.7 keV due to beta decay into Ba137_m1. The Ba137_m1 releases gamma radiation to de-excite. In addition, this source is implemented as an isotropic source.

To obtain the emissions from an arbitrary isotopic source, OpenMC provides the openmc.data.decay_photon_energy feature, which returns the gamma radiation emitted during radioactive decay.

The pulse-height score is determined by the NaI crystal in the detector with the corresponding tally.

Defining the Material¶

Defining the Geometry¶

The following shows the structure of the cylindrical detector.

Defining the Sources¶

The definition of the first source is straightforward. Keep in mind that the energy in the sources is given in the unit of eV.

The definition of the second source is more extensive and makes use of the openmc.data.decay_photon_energy functionality. The openmc.data.decay_photon_energy probability distribution for the different gamma energies, per second per atom. Cs137 has a half-life of 30.17 years, while Ba137_m1 has a much shorter half-life of 2.55 minutes with a branching ratio of 0.946. This implies that 94.6% of Cs137 decays into Ba137_m1.

In the following, the decay scheme of Cs137 can be seen again (public domain).

We assume an equilibrium between Cs137 and Ba137_m1. To calculate the equilibrium amounts of these isotopes, we use the relationship

$\dfrac{T_{1/2}^{p}}{T_{1/2}^{c} \cdot b_{p,c}} = \dfrac{N_{A}^{p}}{N_{A}^{c}}$,

where $T_{1/2}^{p}$ is the half-life of the parent, $T_{1/2}^{c}$ is the half-life of the child, $b_{p,c}$ is the branching ratio, and $N_{A}^{p}$ and $N_{A}^{c}$ are the equilibrium amounts of the parent and child, respectively.

Defining the Settings¶

Defining the Tallies¶

We obtain the values of the pulse-height tallies by utilizing an energy filter with 1000 bins of size 1 keV, ranging from 0 MeV to 1 MeV.

Running the Simulation¶

Extracting the results¶

In a next step we want to extract and plot the results of the pulse height tally.

At this moment, a few points can make things easier to understand:

• We have plotted red vertical lines at for the energies of the strongest sources in the spectrum - the 661.7 keV Cs137 line and the artificial 800 keV line. But also other characteristic peaks of a gamma spectrum can be seen.
• The prominent peak at the first pulse-height value, corresponding to 0 keV, indicates that particles have deposited no energy in the sodium-iodide crystal. This peak does not have a counterpart in real-life detector measurements, as no interaction implies no possibility of measurement. Consequently, we will exclude this peak in subsequent analyses.
• The current simulation statistics are not optimal and could benefit from simulating a substantially larger number of particles.
• In the illustrated example (despite the suboptimal statistics), peaks are displayed as narrow lines. However, in real-world experiments, these peaks are broadened due to the inherent detector resolution. In the next section, we present a method to post-process the pulse-height data to simulate realistic detector signals.

Post processing the results¶

The intrinsic detector resolution is the result of several physical effects such as charge collection statistics or electronic noise. Further it depends on the energy of the incoming particle. The resolution of a detector R(E) is given as the quotient of the full width at half maximum (FWHM) and the energy (E) of the incoming particle:

$R(E) = \frac{\mathrm{FWHM}(E)} {E}$

Usually, this is determined for a detector in several measurements with known sources, and then a function is fitted to these measured values. In this example we will use

$R(E) = \frac{a + b\sqrt{E+cE^{2}}}{E}$.

Further, we assume

$a = 1000$

$b = 4$

$c = 0.0002$

In a next step, we can broaden the results of the pulse-height tallies according to this function. The idea is to go through all the energy bins, determine the respective Gaussian function, and sample according to that. Here we make use of the fact that for a Gaussian

$\sigma(E) = \dfrac{\mathrm{FWHM}(E)}{2 \sqrt{2 \ln{2}}}$.

Current Limitations¶

There are currently still some limitations in the simulation. These will be addressed in the future:

• Some isotopes perform a beta+ decay before emitting a gamma. The resulting positron then annhilates with an electron, producing two 511 keV gamma particles that generate an annhilation peak in the detector. These 511 keV gamma particles have to be inserted manually as source.
• No summation peaks can be simulated yet.
• A comprehensive validation of this method will be published soon.