Reconstruction of multiple Compton scattering events in MeV gamma-ray Compton telescopes towards GRAMS: the physics-based probabilistic model

In the above equations, $P_{\mathrm{path,scat}}(\hat{\mbox{\boldmath$r$}}_{i}\mid\hat{\mbox{\boldmath$q$}}_{i-1})$ and $P_{\mathrm{path,abs}}(\hat{\mbox{\boldmath$r$}}_{i}\mid\hat{\mbox{\boldmath$q$}}_{i-1})$ describe the probability that a gamma ray of $\hat{\mbox{\boldmath$q$}}_{i-1}$ is scattered or absorbed at $\hat{\mbox{\boldmath$r$}}_{i}$. They are defined as

where $\sigma_{\mathrm{abs}}$ and $\sigma_{\mathrm{scat}}$ are the cross sections of photoabsorption and Compton scattering at an energy of $\hat{E}_{i-1}$, respectively; $\rho$ is the number density of the detector material; $\sigma_{\mathrm{all}}$ is the sum of the cross sections of the photon interactions. The term $1/|\hat{\mbox{\boldmath$r$}}_{i}-\hat{\mbox{\boldmath$r$}}_{i-1}|^{2}$ corresponds to the solid angle of the volume element at $\hat{\mbox{\boldmath$r$}}_{i}$ seen from $\hat{\mbox{\boldmath$r$}}_{i-1}$. When incoming gamma rays are considered to be parallel light as in astronomical observations, the term with $i=1$ must be removed, and Eq. LABEL:eq_P_mov is described as

where $\hat{l}_{\mathrm{first}}$ is the length of the gamma-ray path inside the detector before it arrives at $\hat{\mbox{\boldmath$r$}}_{1}$, and depends on the unknown incoming direction ($\hat{\theta}_{0}$, $\hat{\phi}_{0}$) in astronomical cases. Later, we will introduce an approximation of this value (see Eq. 31. Additional position dependence from edge effects is neglected). Note that here it is assumed that the detector consists of a single material. If a Compton telescope consists of several detectors with different materials, e.g., semiconductor detectors and scintillators, then $\sigma_{\mathrm{abs/scat/pair}}$ and $\rho$ also depend on the position. In this case, Eq. LABEL:eq_P_mov is modified as

where $C$ is the straight line from $\hat{\mbox{\boldmath$r$}}_{i-1}$ to $\hat{\mbox{\boldmath$r$}}_{i}$. The same applies to Eqs. LABEL:eq_P_mov_abs and 18.

The function $P_{\mathrm{KN}}\left(\hat{\theta}_{i},\hat{\phi}_{i}\mid\hat{\mbox{\boldmath$q$}}_{i-1}\right)$ corresponds to the probability that a gamma ray of $\hat{\mbox{\boldmath$q$}}_{i-1}$ is scattered to the direction described by $(\hat{\theta}_{i},\hat{\phi}_{i})$. Namely it is the normalized differential cross section of Compton scattering. Following Klein-Nishina’s formula [46], it is defined as

where $r_{e}$ is the classical electron radius and $\hat{E}_{i}$ is the energy of the scattered gamma ray calculated by Eq. 9 and $\hat{\vartheta}^{\mathrm{scat}}_{i}$ is the $i$-th scattering angle defined in Eq. 10 and $\mathrm{d}\hat{\Omega}_{i}$ is equal to $\mathrm{d}(\cos\hat{\theta}_{i})\mathrm{d}\hat{\phi}_{i}$. Note that the polarization is not considered in this process to reduce the amount of calculation because the target position is unknown. The treatment of the polarization with known target position is discussed in [47, 42].

Finally we define $P_{\mathrm{path,esc}}\left(\cdot\right)$ as the probability that a gamma ray escapes from the detector without any interaction. It is described as

where $\hat{l}_{n}$ is the length between $\mbox{\boldmath$\hat{r}$}_{n}$ and the position at which the gamma ray escaped from the detector.

The function $P_{\mathrm{det}}(\cdot)$ corresponds to the measurement process, and then it is determined by the detector response and what kind of quantities a Compton telescope can measure. Here we assume that a Compton telescope measures the deposited energy and position of each interaction independently. In this case, $P_{\mathrm{det}}(\cdot)$ is usually expressed as the product of two detector response functions defined as

where $\hat{\varepsilon}_{i}$ is the true deposited energy at $i$-th interaction, which is equal to $\hat{E}_{i-1}-\hat{E}_{i}$. Here $P_{\mathrm{ene}}(\varepsilon_{\tau(i)}\mid\hat{\varepsilon}_{i})$ is the probability that the true deposited energy of $\hat{\varepsilon}_{i}$ is detected as $\varepsilon_{\tau(i)}$, and $P_{\mathrm{pos}}(\mbox{\boldmath$r$}_{\tau(i)}\mid\hat{\mbox{\boldmath$r$}}_{i})$ is one that the interaction position of $\hat{\mbox{\boldmath$r$}}_{i}$ is detected as $\mbox{\boldmath$r$}_{\tau(i)}$. In general, these two can have any functional forms and they should be assumed so as to model the detector response adequately using experiments with calibration sources, simulations or simplifications. We will introduce specific formulas for these functions in the following section.

To describe the detector response terms $P_{\mathrm{ene}}(\varepsilon\mid\hat{\varepsilon})$ and $P_{\mathrm{pos}}(\mbox{\boldmath$r$}\mid\hat{\mbox{\boldmath$r$}})$, we adopt the Gaussian function as a simple model. Namely, these functions are formulated as

where $\sigma_{\hat{\varepsilon}}$, $\sigma_{\hat{x}}$, $\sigma_{\hat{y}}$, and $\sigma_{\hat{z}}$ are the energy and positional resolutions of the detector, respectively.

The response functions in Eqs. 23 and 24 are very simple, e.g., Eqs. 23 considers only the energy resolution of the detector. We note that one may adopt more complex models in principle. For example, the positional resolution often depends on the interaction position due to the charge sharing effect, pixelization, e.t.c. These effects may be modeled by including the dependence of $\sigma_{\hat{x}/\hat{y}/\hat{z}}$ on the interaction position. While these effects vary with detector systems, and we ignore them here, it may improve the algorithm performance by using more realistic positional/energy response functions.

1. 1.

   The incoming gamma-ray energy is calculated as the sum of the detected energy (Eq. 25), and the interaction positions are assumed to be the same as the detected ones (Eq. 26).

2. 2.

   The first scattering angle is calculated by Eq. 30 and the incoming gamma-ray direction is constrained on a Compton circle in the sky.

3. 3.

   The probability is calculated using Eq. 11 by applying the approximations described in Eqs. 29, 31

1. 1.

   The incoming gamma-ray energy is calculated by Eq. 32, and the interaction positions are assumed to be the same as the detected ones (Eq. 26).

2. 2.

   The first scattering angle is calculated by Eq. 30 and incoming gamma-ray direction is constrained on a Compton circle in the sky.

3. 3.

   The probability is calculated using Eq. 14 by applying the approximations described in Eqs. 29, 31, 37, 38.

As mentioned in §2, when the number of hits is two or fewer, the escape gamma-ray energy cannot be estimated since Eq. 32 requires the detection of three or more hits to calculate the second and subsequent scattering angles. Thus, this algorithm can work for 3 or more hit events.

The derived algorithm differs from the MSD methods in that it compares the deposited energy rather than the scattering angle. Effectively, this algorithm calculates the deposited energy at each site assuming the incident gamma-ray energy and interaction positions, and compares it with the measured value directly through the energy response function. Moreover, by adopting the probabilistic method, the physical processes are naturally taken into account in the scattering order determination.

We assume a detector with a size of $140\times 140\times 20$ cm³ filled with liquid argon as shown in Figure 4. It is the same as proposed in [10]. The energy resolution $\sigma_{\hat{\varepsilon}}$ of the detector is set to as follows [10]:

Eq. 40 is also adopted in the probability calculation in Eq. 23. We note that while the energy resolution of liquid argon time projection chambers in the MeV range has not yet been measured, the energy resolution of $\sim 1$% ($1\sigma$) at 1 MeV is measured with a liquid argon scintillation detector [51]. Thus, the argon detector can have the assumed energy resolution in principle.

The X-Y positions of signals are pixelized with a size of 2 mm, and we set $\sigma_{\hat{x}}$ and $\sigma_{\hat{y}}$ to be $2/\sqrt{12}$ mm in Eq. 24. This is the standard deviation of X or Y positions of events distributed uniformly in a single pixel. The direction of electron drift in the time projection chamber is assumed to be parallel to the Z-axis, and the resolution of Z position is assumed to $\sigma_{\hat{z}}=$ 1 mm. Note that X and Y positions are determined by the pixel positions, but Z position is estimated from the measured electron drift time. Thus $\sigma_{\hat{z}}$ can be different from $\sigma_{\hat{x}/\hat{y}}$. In this simulation, the signals from adjacent pixels are merged into a single signal. The signal position after merging is set to the center of gravity (energy-weighted average) of the pixels before merging. Note that the electron diffusion in the liquid argon is small; in this case, it is about 300 $\mu$m at the most, much smaller than the pixel size [52, 53]. Thus we ignored the effect of electron diffusion in this simulation. The energy threshold of each pixel is set to 25 keV.

In this demonstration, we set $L_{\mathrm{first}}$ in Eq. 31 and $L_{\mathrm{esc}}$ in Eq. 37 to be 10 and 20 cm respectively (see A for this optimization). The computational time for different numbers of hits is described in Table 1. For this table, we used 1 MeV gamma-ray events in the following subsection.

We simulated $10^{8}$ events of 1 MeV gamma-ray beam incoming from the top of the detector, i.e., $\hat{\theta}_{0}=\pi$. The beam has a radius of 110 cm that covers the whole volume of the detector, and it is co-aligned at the detector center. Figure 5 shows the count of the detected events with the number of hits. The number of counts of fully-absorbed events reaches a maximum at 4 hits. Here we focus on the events with the number of hits from 3 to 8. The ratio of the number of the events with more than 8 hits to the total detected events is only 0.13% and they are negligible. Figure 6 shows an example of applying the algorithm to a 3-hit event. The algorithm outputs 12 probabilities corresponding to each scattering order and event type. In this case, the largest value is the leftmost one, and this event is identified as a fully-absorbed event with a scattering order of “123”.

First, we examine the algorithm performance against fully-absorbed events qualitatively. Figure 7 shows the energy spectra of all detected events (blue, solid) and of those classified as fully-absorbed events after applying the algorithm (red, solid). We confirmed that the events peaking around 1 MeV is correctly classified as fully-absorbed events and the component below 1 MeV, which corresponds to escape events, is reduced successfully after the reconstruction algorithm is applied. As the number of hits is increased, the events around 1 MeV are classified as fully-absorbed more accurately. This is because events with more hits have more physical information to constrain the Compton scattering sequence.

Next, the algorithm performance against escape events is checked. Figure 8 shows energy spectra of events classified as escape events. The orange lines are spectra of the sum of detected energies of events classified as escape events. Note that the blue ones are the same as Figure 7. By comparing the blue and orange lines, we can see how accurately the algorithm identifies the escape events. Considering that the events with total deposited energies less than 1 MeV are escape events, and most of these events are included in the orange line, we confirmed that the algorithm successfully identifies the escape events. Furthermore, we can also check whether the energy correction by Eq. 32 works for the escape events correctly. The red lines in Figure 8 are the spectra of estimated incident gamma-ray energy by applying Eq. 32 to the events classified as escape events, i.e., the events in the orange lines. The peak at 1 MeV is clearly reconstructed, confirming that the algorithm correctly estimates the escape energy. Note that the reconstructed spectra of the escape events are broader than those of fully-absorbed events. It is because the positional information is also used for estimating gamma-ray escape energy. Then the reconstructed energy is affected by both the errors of position and energy of each detected hit.

The dependence of the spectra on the event type and the number of hits is investigated quantitatively. Figure 9 shows the full width at half maximum (FWHM) of the reconstructed energy spectra with different numbers of hits. Note that the detector has an energy resolution of 16.6 keV at 1 MeV as FWHM in this simulation. For fully-absorbed events, the energy resolution does not depend on the number of hits so much. On the other hand, for escape events, the FWHM for the escape events shows a peak at 5 hits. When fully-absorbed events are misidentified as escape events, their incident energies are incorrectly reconstructed as above 1 MeV by the energy correction. Then, this misidentification makes a high-energy tail or hump in the reconstructed energy spectra, and from 5 hits, the FWHM catches it. As mentioned before, the positional uncertainty can affect the energy spectra of the escape events, which also contributes to their tail-like structures in both low- and high-energy bands. To evaluate it, we also show the full width at tenth maximum (FWTM) of the spectra in Figure 9. We can see that the FWTM of the escape events becomes smaller as the number of hits increases, which corresponds to the spectrum having shorter tails. We want to note that though the absolute value of the energy resolution varies by the assumed detector response, in general, the energy resolutions of the fully-absorbed and escape events are different, and they have different dependence on the number of hits.

We also investigated the angular resolution of the reconstructed events and its dependence on the number of hits. The angular resolution is evaluated by the angular resolution measure (ARM), which is defined as

where $\theta_{K}$ is the estimated first scattering angle:

and $\theta_{G}$ is the first scattering angle calculated from the incident direction and the reconstructed positions of the first and second hits.

Figure 10 shows distributions of ARM for both fully-absorbed and escape events up to 8 hits, using the same data set in the previous subsection. We normalized the histograms by the total count in each. In both cases, as the number of hits increases, the main peak at 0 degrees becomes sharper, and the tail components and a sub-peak at $\sim 90$ degrees are reduced. It suggests that the scattering order is estimated more accurately with more hits. We will examine this point quantitatively in the following subsection.

The FWHM of the obtained ARM distributions are shown in Figure 11. It can be seen that they are from 4.9 to 7.1 degrees, depending on the number of hits and the event type. In this case, the FWHM is larger for fully-absorbed events with a smaller number of hits, because these events contain a higher percentage of events with large scattering angles and large deposited energies at the first interaction, and the Doppler broadening effect limits the angular resolution. When 1 MeV gamma ray is back-scattered, the scattered gamma ray has an energy of 0.204 MeV. Then the gamma ray is difficult to escape from the detector since its interaction length is much smaller than that of 1 MeV gamma ray, and it is usually absorbed by the detector after few scattering. Note that since the scattering angle at the first interaction is estimated, the angular resolution can be improved by selecting forward-scattering events. For example, when using events with scattering angles less than 60 degrees at the first interaction, the ARM FWHM of 3-hit fully-absorbed events is improved from 7.1 degrees to 4.1 degrees. Figure 11 also shows the FWTM of the ARM distributions.

In order to quantify the performance of the developed algorithm, we calculate the fraction of the 1MeV gamma-ray events reconstructed with the correct event types (“accuracy (type)” in the following figures), the fraction of the events reconstructed with the correct scattering orders (“accuracy (order)”), and the fraction of the events reconstructed correctly for both the event type and scattering order (“accuracy (type + order)”).

We show these quantities for the 1 MeV gamma-ray events in Figure 12. The algorithm identifies the event type correctly with an accuracy of more than 73%, and when the number of hits reaches 8, the accuracy is improved to 95% (“$+$” in the figure on the left). On the other hand, the accuracy for the scattering order is about 55%, and it reaches the peak at 4-5 hits (the dotted line with “$\times$” in the figure on the left). This is because as the number of hits increases, the number of candidates for the scattering order increases with factorial. Then the correct order must be selected from a larger number of possibilities. However, if we focus on the scattering orders of the first two or three interactions (the solid and dashed lines, respectively), the accuracy of the scattering order improves with the number of hits, and it reaches 80% for 8-hit events. From a practical point of view, it is important to identify the first two or three interactions. In the case of fully-absorbed events, the scattering angle and scattered gamma-ray direction at the first interaction can be calculated correctly as long as the first two interactions are correctly identified. Also, if the first three interactions are correctly identified for escape events, the escape energy can be calculated in principle [30]. The accuracy for both event type and scattering order also has a similar trend (the figure on the right). When only considering the first two interactions for fully-absorbed events and the first three for escape events, the algorithm reconstructs the gamma-ray events with an accuracy of 81% for 8-hit events (see the solid line in the right).

In Figures 13 and 14, we extract the truly fully-absorbed or escape events respectively, and show the fractions of correctly reconstructed events among them. For the escape events, we use events in which the sum of the truly deposited energies is less than 900 keV. The accuracy of the fully-absorbed events shows a similar trend to Figure 12. The accuracy for the event type is excellent for 8-hit events while this decreases to 56% for 3-hits events. However, for truly escape events with 3 hits, the event type is correctly identified with an accuracy of 81%. Note that the algorithm performance depends on the detector’s energy resolution and position resolution. It is discussed in B.

We compare the obtained performance with those of other two algorithms. We adopt the TANGO algorithm [42] as reference, which can distinguish fully-absorbed and escape events based on the MSD method. We slightly modified the reference algorithm to compare the cosine of the scattering angle directly rather than the scattering angle for the following reason. The assumed energy resolution here is not as good as the germanium detector in the original TANGO algorithm. Then the cosine of the scattering angle calculated from kinematics (see Eq. 9 and $|\cos\theta_{kin}|$ in [42]) is often larger than one even for a correct scattering order, mainly when the events include a Compton scattering interaction with a scattering angle close to 0 or 180 degrees. In the original implementation, the scattering angle, not the cosine of it, is used for the FoM, and the correct scattering order of such an event is rejected because the arc-cosine cannot be applied to them. We found that it degrades the algorithm’s performance significantly, and using the cosine of scattering angle results in much better performance because we do not have to use the arc-cosine and avoid rejecting the correct scattering order due to the mathematical reason. Also, since the original implementation assumes the source position in their nuclear experiments, we tuned the TANGO algorithm as follows. The comparison of the cosine of scattering angle is skipped at the first interaction and the polarization effect in the Klein-Nishina formula is ignored, and the traveled distance in Eq.6 is fixed to 10 cm for the first interaction. In addition, we also compare our algorithm with a classical algorithm [29] that does not identify the escape events.

Figure 15 shows the accuracy of the reconstruction by our algorithm (this work), the modified TANGO algorithm (“MSD” in the figure), and the classical algorithm when applying them to the 1 MeV gamma-ray simulation data. The black points for the proposed algorithm are the same as Figure 12. This algorithm outperforms the others in the accuracy for the event type. A plausible reason for this improvement is that the ratio of the figure-of-merits of the escapes to fully-absorbed events is calculated more accurately by formulating the probability model. Note that here the classical method reconstructs any event as a fully-absorbed event, and the accuracy for the event type of the classical method is equal to the ratio of truly fully-absorbed events to the total.

As for the accuracy for the scattering order, the TANGO algorithm outperforms very slightly above 6 hits because when events are dominated by fully-absorbed events the MSD algorithm tends to reconstruct small-scattering events with a slightly better accuracy [40]. However, below 6 hits, our method has a better accuracy with a maximum improvement of about 11% at 3-hit events. When only considering the first two interactions for fully-absorbed events and the first three for escape events, our algorithm shows the best accuracy for both event type and scattering order (see the bottom in the figure).

This advantage is clear particularly for small-number-hit events, namely, 3-hit events. Figures 16 shows the reconstructed energy spectra and ARM distributions of 3-hit events. The difference to the modified TANGO algorithm as well as the classical one is noticeable in the energy spectrum, where the peak height differs by about 40%. This improved performance is also observed in the ARM distribution as the sharpest peak around $0^{\circ}$.

To further investigate the proposed algorithm’s performance for different incoming gamma-ray energies, we also calculate the reconstruction accuracy of gamma rays from 500 keV to 5 MeV. The results are shown in Figure 17. The accuracy gets the maximum at 1–2 MeV though it depends on the number of hits. A factor limiting the performance at 5 MeV would be recoil electrons with energies above few MeV, which produce gamma rays via bremsstrahlung. These gamma rays are not considered or filtered out in the algorithm. Note that in this analysis we removed events that include pair creation to investigate the performance against truly Compton scattering events. In the argon detector, at $\sim 12$ MeV the cross section of pair creation becomes comparable to Compton scattering. The fraction of events that include pair creation is 2.3% at 2 MeV and 21% at 5 MeV. In reality, the accuracy of removing such events also affects the performance. We will discuss it in Section 7.

In astronomical observations, MeV gamma-ray signals are usually dominated by various background events, e.g., charged particles and neutrons of cosmic rays, radioactivation, albedo gamma-rays, and galactic/extra-galactic diffuse background. Eliminating the background is critical for the achievement of high sensitivity. For example, in the GRAMS project, the neutron background is removed by pulse shape discrimination, and surrounding plastic scintillators identify the charged particle events. Towards the high sensitivity, future works are needed to study how significantly the background contaminates signal events and whether the background events can be distinguished in some way based on background-included simulations.

As a possible idea, the conditional probability obtained from the algorithm ($\mathcal{P}_{\mathrm{fullabs}}$ or $\mathcal{P}_{\mathrm{escape}}$) could be used for background reduction. The black line in Figure 18 shows the distribution of the conditional probability of 3-hit fully-absorbed events using the 1 MeV gamma-ray simulation data set in Section 6.2. To generate events that do not originate from Compton scattering of gamma rays by a simple way, we use the same data set, but for each hit in all events we re-sampled a position from a uniform distribution in the detector. Then, we applied the reconstruction algorithm to these position-randomized events. These events can be interpreted as a very simplified model of the background gamma rays, e.g., produced from nuclear interactions with cosmic rays [54]. The red line in the figure corresponds to the distribution of the obtained conditional probability of the position-randomized events. They have much smaller values than the gamma-ray events. This result indicates that some background events can be removed by setting a threshold in the obtained value, though a more realistic simulation is needed to validate this application.

This work mainly focuses on astronomical application, but our algorithm could be used in other fields, e.g., imaging in nuclear medicine therapy or monitoring high-intensity radiation fields [55, 56, 57]. In these cases, incoming gamma-ray energy is often known a priori. The reconstruction algorithm can be modified by fixing the incoming energy instead of estimating it by Eq. 25 or 32. This additional information can improve the reconstruction performance. Figure 19 shows the accuracy of 1 MeV gamma-ray event when fixing the incoming energy. The reconstruction accuracy is improved, especially for events with a small number of hits by $\sim 30$%.

Though this paper focuses on the case that only the deposited energies and positions are measured, our approach could be extensible even when a Compton telescope can also measure the trajectories of recoiled electrons and estimate their initial momentum directions [14, 15, 16]. In this case, $\mbox{\boldmath$D$}_{I}$ is needed to be redefined as

where $\mbox{\boldmath$p$}_{\mathrm{e},I}$ is the measured momentum direction of a recoiled electron. Then, we modify the detector response function $P_{\mathrm{det}}\left(\mbox{\boldmath$D$}_{I}\mid\hat{\mbox{\boldmath$q$}}_{i-1},\hat{\mbox{\boldmath$q$}}_{i}\right)$ by including the electron trajectory information. Here we introduce a function $P_{\mathrm{track}}(\cdot)$ that compares the expected direction of a recoiled electron $\left(\hat{\mbox{\boldmath$p$}}_{i}-\hat{\mbox{\boldmath$p$}}_{i-1}\right)$ and measured one $\left(\mbox{\boldmath$p$}_{\mathrm{e},I}\right)$. The definition of $P_{\mathrm{track}}(\cdot)$ would depend on the configuration of Compton telescopes because the obtained electron images vary by the detectors, and in general $P_{\mathrm{det}}(\cdot)$ is reformulated as

After modifying $P_{\mathrm{det}}(\cdot)$, the scattering order can be determined in the same way as in Section 5 but with one exception. In this case, even if the incoming gamma-ray direction is constrained on a Compton circle, the conditional probability depends on location in the circle because the direction of electron recoil varies depending on it. Thus, with this additional information, one can constrain the incoming gamma-ray direction as a much narrower region, as the merit of electron-tracking Compton telescopes.

The developed algorithm successfully reconstructs multiple Compton scattering events with high accuracy. However, we should note that we ignored several important factors in the probabilistic model, e.g., interaction with passive materials, Doppler broadening effect, and bremsstrahlung. In the rest of this paper, we discuss them and a possibility for further improvement of this algorithm.