{'fraction_non_alphanumeric': 0.031755753242179725, 'fraction_numerical': 0.012034726156390195, 'mean_word_length': 4.607142857142857, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Reconstruction of multiple Compton scattering events in MeV gamma-ray Compton telescopes towards GRAMS: the physics-based probabilistic model

Aimed at progress in mega-electron volt (MeV) gamma-ray astronomy, which has not yet been well-explored, Compton telescope missions with a variety of detector concepts have been proposed so far. One of the key techniques for these future missions is an event reconstruction algorithm that is able to determine the scattering orders of multiple Compton scattering events and to identify events in which gamma rays escape from the detectors before they deposit all of their energies. We revisit previous event reconstruction methods and propose a modified algorithm based on a probabilistic method. First, we present a general formalism of the probabilistic model of Compton scattering describing physical interactions inside the detector and measurement processes. Then, we also introduce several approximations in the calculation of the probability functions for efficient computation. For validation, the developed algorithm has been applied to simulation data of a Compton telescope using a liquid argon time projection chamber, which is a new type of Compton telescope proposed for the GRAMS project. We have confirmed that it works successfully for up to 8-hit events, including correction of incoming gamma-ray energies for escape events. The proposed algorithm can be used for next-generation MeV gamma-ray missions featured by large-volume detectors, e.g., GRAMS.


Introduction
Astronomical observations of gamma rays from a few 100 keV to a few 10 MeV remain to be explored in modern astronomy. COMPTEL, aboard the Compton Gamma-Ray Observatory, observed the MeV gamma-ray universe in 1991-2000 [1] and found more than 30 sources in the energy band from 0.75 to 30 MeV [2]. After this mission, a few space satellite missions, e.g., INTEGRAL/SPI [3] have succeeded, but the number of the detected sources in the MeV band is still limited. These days, stimulated by the dawn of multi-messenger astronomy, including gravitational wave [4] and high-energy neutrino observations [5], this curtained window of electromagnetic waves is drawing increasing attention. Towards high-sensitivity observations in the 2020s and 2030s, several MeV gamma-ray missions have been proposed at this moment: COSI 1 [6], SMILE [7], e-ASTROGAM [8], AMEGO [9], GRAMS [10], etc.
All of the missions above utilize the Compton telescope, one of the most promising techniques for imaging gamma-ray sources in the sub-MeV/MeV bands [11,12,13]. A Compton telescope measures the position and deposited energy at each interaction, and calculates the scattering angle by the kinematics of Compton scattering. Then, the incoming gamma-ray direction is constrained to a ring in the sky, traditionally called an "event circle", and in the following referred to as a "Compton circle". If recoiled electron trajectories can be measured additionally, the gammaray direction is constrained on an arc-shaped region [14,15,16]. After an accumulation of many events, the incoming gamma-ray direction can be identified as intersections of the constrained circles or arcs. This basic principle of Compton imaging [17], called a back-projection method, has been superseded by statistical methods for image reconstruction [18,19,20,21,22,23,24,25].
To calculate the scattering angle of a gamma-ray event, it is required to accurately identify the scattering order of the detected signals and estimate the incident gamma-ray energy. When a gamma ray deposits its total energy, it is sufficient to identify the first and second interactions. If the third interaction is identified additionally, the incident energy can be estimated even if a gamma ray escapes outside the detector before they are absorbed, as discussed later. The scattering order determination becomes complicated as the gamma-ray energy increases since gamma rays can easily be scattered multiple times in a detector and can escape from a detector. These multiple scattering events are con-sidered to be dominant in multi-layered or largevolume detectors. GRAMS is an example since it utilizes a thick and wide-area liquid argon time projection chamber as a high-efficiency Compton telescope [10]. In this project, like many future missions, it is difficult to measure the time-of-flight between interactions of scattered gamma rays considering its timing resolution because the length between interactions gets shorter than a detector with two separated layers like COMPTEL. In many future missions, the scattering order is difficult to determine directly, and thus, one has to determine it based on the detected energies and positions. This paper focuses on this order determination problem of the multiple scattering events.
Several approaches have been studied in a variety of fields, e.g., astronomy [26,27,28,29,30,31,32,33], homeland security [34,35,36,37,38,39,40,41], and nuclear physics [42,43,44]. A widely used method is the mean squared difference (MSD) method [27,40]. For events with n (≥ 3) interactions, the scattering angles of all the interaction sites except for the first and last interaction points, i.e., i-th interactions with 2 ≤ i ≤ n − 1, can be calculated in two ways, namely-one is based on kinematics and the other on geometrical information. Then, by comparing the calculated angles at each site, the most plausible scattering order is estimated using a figure-of-merit based on the MSD of the calculated scattering angles. Another similar and well-formulated method is to calculate the probability of Compton scattering, assuming the scattering order, and choose the one with the highest probability (deterministic method [40]). A Bayesian method [31] and a neural network method [32] are also proposed, and sometimes they outperform the above approaches. However, these approaches require large simulation data sets and a lot of computer resources. One should be careful that when the detector properties change, it may be required to re-compute the training process in a neutral network method.
Usually, events such that gamma rays escape from the detector (escape events) are rejected in the analysis. However, as gamma-ray energy increases, photoabsorption becomes less likely, and the escape events become non-negligible to achieve high detection efficiency. Then, reconstruction methods should also estimate the escape energy and distinguish the events from those that deposit all of their gamma-ray energies in the detector (fully-absorbed events). If the number of interactions is three or more, the gamma-ray energy of escape events can be estimated by combining the position and energy information (see, e.g., [28,30]). For example, using the MSD method, [42] proposed an algorithm for germanium spectrometers in nuclear physics experiments [44]. They demonstrated that the gammaray detection efficiency can be increased by identifying escape events and estimating their escape energies correctly. This consideration is also essential for astronomical observations.
Towards the future missions, especially GRAMS, in this work, we revisit the classical approaches and formulate an event reconstruction algorithm that performs full treatment of both the fully absorbed events and the escape events based on a probabilistic method. To derive the quantities for the event reconstruction deductively, first, we formulate the probability functions of the event occurrence by considering physical and measurement processes in Compton telescopes. In Section 2, we briefly review the problem and describe the basic concept of the present algorithm. The mathematical formulas of the probability functions are shown in Section 3. In Section 4, we implement the algorithm by introducing several approximations to calculate the derived functions' values efficiently. Here we consider a Compton telescope that measures the interaction positions and deposited energies, in which information on the recoil electron trajectories is not contained. A summary of the developed algorithm is given in Section 5. As a numerical test, we apply it to simulation data sets of a liquid argon detector, regarded as a ideal concept model of the GRAMS project, in Section 6. Finally, we summarize our results and discuss further improvements in Section 7.
We have developed the present algorithm to conduct a proof-of-concept study of the GRAMS project, which newly proposes a large-volume liquid argon Compton telescope. The algorithm will be considered as a standard event reconstruction algorithm in the project's data analysis pipeline. Noted that we are also developing another type of a reconstruction algorithm using a multi-task deep neural network as an alternative approach. A detailed description of the algorithm will be presented by a forthcoming publication, which will also show the performance comparison between the neural network method and the physics-based probabilistic method in this paper.

Basic Concept
In a Compton telescope, ideally an incident gamma ray deposits its energy at interaction sites via Compton scattering or photoabsorption. In general, the measured values at each site are expressed by a tuple of several physical quantities, which is denoted by D I . The index I is the label of each tuple in an event that corresponds to a gamma-ray incidence into the telescope, and the scattering order is unknown at this moment. We define a hit as a measured interaction with D I . When a Compton telescope measures deposited energies and interaction positions, D I is a pair of them: where r I and ε I are the measured position and energy, respectively. Since a gamma ray is scattered multiple times or absorbed by a detector, we obtain a list of D I : When the number of hits in an event is n, we refer to the event as an n-hit event. Note that we do not consider events including pair creation in this work, though they are not negligible at gamma-ray energies above ∼ 5 MeV, as discussed in §7.
In the reconstruction of Compton scattering events, the task is divided into the following: 1. to determine the event type, i.e., a fullyabsorbed event or an escape event, for a given event. 2. to determine the scattering order of the detected hits: where (D I ) ordered is a re-ordered list of (D I ) by determining or assuming the scattering order. Here we define the function τ (·) which maps the scattering order to the data label (I), i.e., i-th interaction corresponds to D τ (i) . Figure 1 shows the scattering order candidates and their corresponding map functions τ (·) for a 3hit event. 3. to estimate the incoming gamma-ray energy. 4. to estimate the incoming gamma-ray direction. Related to the second task, the MSD method [26,27,29,40,42] determines the scattering order by calculating the scattering angles redundantly Scattering order candidates of a fully-absorbed event Scattering order candidates of an escape event from kinematics and geometrical information. For example, [27] defines the following quantity: where ϑ geo i and ϑ kin i are i-th scattering angles calculated by kinematics and geometrically, respectively; ∆ cos ϑ i is the measurement uncertainty, and N is the number of hits. This quantity is interpreted as a modified chi-squared value. For all N ! scattering order candidates, χ 2 c are calculated, and the best scattering order is determined as one that yields the smallest χ 2 c . Related to the energy reconstruction, the socalled three-Compton method is usually used to estimate escape gamma-ray energy [28,30]. If the energy loss at the first and second sites (ε 1 , ε 2 ) and the scattering angle at the second interaction (ϑ 2 ) are known, then the incident gamma-ray energy (Ê 0 ) can be calculated aŝ where m e is the electron rest mass and c is the speed of light. For example, [42] calculates the incident gamma-ray energy for escape events using the three-Compton method and defines figure-of-merits for both fully-absorbed and escape events when the gamma-ray generation position is known. In this work, aiming at MeV gamma-ray telescopes in astronomy, we formulate a probabilistic model of Compton scattering in a detector instead of the MSD methods mentioned above. The most plausible event type and scattering order can be determined as those that yield the maximum probability. While the basic concept is similar to the deterministic methods (e.g., [40]), we derive a probabilistic model that describes both fully-absorbed and escape events in a unified framework. Then, our approach can also identify and reconstruct escape events which have been usually rejected. In our probabilistic model, the physical processes of gamma rays and the measurement processes of the detector are explicitly considered.
The main advantage of this approach is that the probabilistic model derives the figure-of-merit for the event reconstruction deductively. The relative ratio between the figure-of-merits for the fullyabsorbed and escape events can be determined based on the physics-based modeling. It is expected that this approach makes the event type determination more accurate compared to heuristic approaches in the MSD methods. We will examine this point in §6.5. Additionally, when more information is obtained, the proposed algorithm can be naturally extensible by including probability functions corresponding to the new information (see Section 7.2).

Probabilistic model of an event type and scattering order
In this section, we formulate a probabilistic model of a sequence of Compton scattering and photoabsorption in a detector, given an event type and scattering order. Here we assume that an incoming gamma ray is Compton-scattered in the detector n − 1 times and photo-absorbed at last, or is scattered n times and escapes from the detector, i.e., a n-hit event. In addition, we assume that all interactions are measured and a list of the measurement values (D I ) = (D 1 , D 2 , ..., D n ) are obtained. Note that in reality, there is a case that some interactions are not detected due to interactions in passive materials, the detection threshold, and multiple scattering within the same spatial resolution element of the detector. Although the probability of such events strongly depends on the actual detector configuration, we ignore these possibilities, and as a result, our algorithm treats those events as fully-absorbed or escape events. Especially, the effect of passive materials on the algorithm performance is important for realistic situations, and thus we will discuss it later in the discussion section (see §7.3.1).
As shown in Figure 2, the gamma ray changes its energy and its direction of travel every time it interacts with the detector. Then, the gamma ray after the i-th interaction in the detector can be described with the quantitiesq i defined aŝ wherer i represents the i-th interaction position;Ê i andp i represent the energy and momentum vector, respectively, of the gamma ray after the i-th interaction (see Figure 2);θ i andφ i describe the direction of travel of the gamma ray. We refer to the quantityq i as gamma-ray state in this work. To distinguish explicitly the parameters in the gammaray state and those measured by experiments, we put hats on the former ones. Note thatq 0 represents the initial gamma-ray state, and the polarization state of incoming and scattered gamma-rays is neglected in this study.
In astronomical observations,q 0 is described by just three parameters: since incoming gamma rays originate from distant sources. In astronomical cases, the incoming gamma-ray direction and energy are usually unknown, and rather an important goal is the reconstruction of the scattering order. In this section, we formulate the probabilistic model assuming the initial gamma-ray state, and in §4.2 we explain the approximations implemented in the current algorithm considering that the initial gamma-ray state is unknown in reality. The gamma-ray stateq i (i ≥ 1) can be determined when we assume the initial gamma ray statê q 0 and the position of each interaction (r i(≥1) ). Except for i = n,θ i andφ i are determined because the momentum of a scattered gamma ray is parallel tô The gamma-ray energyÊ i is determined by the kinematics of Compton scattering: whereθ scat i is the i-th scattering angle which is calculated as Note that Eq. 9 does not take account of the uncertainty ofÊ i due to the finite momentum fluctuation of the target electrons in the detector material, which is also known as the Doppler broadening effect [45]. We will discuss a possible treatment to include this effect in Section 7.
To describe the probabilistic model, we show the graphical representation of a Compton scattering event in Figure 3. The change of the statê q i is a probabilistic process determined by the physics of Compton scattering. Then, the parameters of the state are related to the measured values (D I ) through measurement processes. For example,Ê i−1 −Ê i is measured as the deposited energy at i-th interaction.
As a general expression, given a fully-absorbed event, its initial stateq 0 and the scattering order map τ (·), the conditional probability to produce its interaction positionsr i (1 ≤ i ≤ n − 1) and the data (D I ) is described as The gamma-ray energyÊ i and the direction of travel (θ i ,φ i ) after the i(≥ 1)-th interaction are determined from the model parameters (see Eqs. 8 and 9).
Here P scat q i , D τ (i) |q i−1 is the probability that a gamma ray with a given stateq i−1 is scattered with changing its state toq i , and the interaction is measured as D τ (i) . It is described as the product of three functions: where P path,scat (·) is the probability that a gamma ray ofq i−1 is scattered atr i ; P KN (·) is the probability that the scattering direction is (θ i ,φ i ) and it is determined by the Klein-Nishina formula [46]; P det (·) is a probability function related to the detector response. We will define it in §3.2. The first and second functions correspond to the probability to produceq i and the last one calculates the probability to obtain D τ (i) givenq i−1 andq i . Also, P abs q n , D τ (n) |q n−1 is the probability that a gamma ray ofq n−1 is absorbed atr n , and the in-teraction is measured as D τ (n) : P abs q n , D τ (n) |q n−1 = P path,abs r n |q n−1 P det D τ (n) |q n−1 ,q n , where P path,abs (·) is the probability that a gamma ray withq i−1 is absorbed atr i . In the following subsections, we explain these functions in detail. Besides, the conditional probability corresponding to the escape event is described as The only difference between Eq. 11 and Eq. 14 is the treatment of the last interaction. Here we introduce P esc q n , D τ (n) |q n−1 , which is the probability that a gamma ray with a given stateq n−1 is scattered with changing its state toq n and escape from the detector, and the interaction is measured as D τ (n) . It is described as P esc q n , D τ (n) |q n−1 = P path,scat r n |q n−1 × d(cosθ n )dφ n P KN θ n ,φ n |q n−1 ×P det D τ (n) |q n−1 ,q n P path,esc (q n ) , where P path,esc (q n ) is the probability that a gamma ray with a state ofq n escapes from the detector. Since we cannot know the momentum direction of the escape gamma ray, we integrate the functions overφ n and cosθ n .
3.1. Probability functions related to physical processes In the above equations, P path,scat (r i |q i−1 ) and P path,abs (r i |q i−1 ) describe the probability that a gamma ray ofq i−1 is scattered or absorbed atr i . They are defined as where σ abs and σ scat are the cross sections of photoabsorption and Compton scattering at an energy ofÊ i−1 , respectively; ρ is the number density of the detector material; σ all is the sum of the cross sections of the photon interactions. The term 1/|r i −r i−1 | 2 corresponds to the solid angle of the volume element atr i seen fromr 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. 16 is described as wherel first is the length of the gamma-ray path inside the detector before it arrives atr 1 , and depends on the unknown incoming direction (θ 0 ,φ 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 σ abs/scat/pair and ρ also depend on the posi-tion. In this case, Eq. 16 is modified as where C is the straight line fromr i−1 tor i . The same applies to Eqs. 17 and 18. The function P KN θ i ,φ i |q i−1 corresponds to the probability that a gamma ray ofq i−1 is scattered to the direction described by (θ i ,φ 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Ê i is the energy of the scattered gamma ray calculated by Eq. 9 andθ scat i is the i-th scattering angle defined in Figure 3: A graphical representation of Compton scattering in a detector and measurement in a Compton telescope. The red arrows correspond to probabilistic processes that take place only in escape events while the blue arrow corresponds to that for fully-absorbed events. Note again thatq 0 is unknown and |r 0 | → ∞ in astronomical situations. The approximations implemented for astronomy are explained in § 4.2.
Eq. 10 and dΩ i is equal to d(cosθ i )dφ 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 path,esc (·) as the probability that a gamma ray escapes from the detector without any interaction. It is described as wherel n is the length betweenr n and the position at which the gamma ray escaped from the detector.

Probability functions related to measurements
The function P det (·) 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 det (·) is usually expressed as the product of two detector response functions defined as whereε i is the true deposited energy at i-th interaction, which is equal toÊ i−1 −Ê i . Here P ene (ε τ (i) | ε i ) is the probability that the true deposited energy ofε i is detected as ε τ (i) , and P pos (r τ (i) |r i ) is one that the interaction position ofr i is detected as r τ (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.

Implementation
The task of this algorithm is to determine the event type, the scattering order and the incident energy and to constrain the incoming direction (see §2). In our approach, these parameters are estimated as those that yield the maximum probability defined as Eq. 11 or Eq. 14. Hereafter we notate the estimated incident energy and incoming direction asÊ * 0 and (θ * 0 ,φ * 0 ), respectively. Ideally, this task is achieved by sweeping the parameter space 8 of the incoming gamma-ray energy and direction, and marginalizing out the true interaction positions (r τ (i) ). However, it is not practical in terms of computational time because the parameter space becomes of high dimensions, i.e., (3n + 3) dimensions when the number of interactions is n.
In this section, we implement the order determination algorithm by focusing on a Compton telescope that measures the interaction positions and deposited energies and does not obtain the trajectories of recoiled electrons. First, we introduce specific forms for the detector response terms. Then, we introduce several approximations to eliminate integral calculations and parameter space search.

Detector response
To describe the detector response terms P ene (ε | ε) and P pos (r |r), we adopt the Gaussian function as a simple model. Namely, these functions are formulated as where σε, σx, σŷ, and σẑ 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 σx /ŷ/ẑ 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.

Approximations in the calculation of the probability functions 4.2.1. Fully-absorbed events
For the fully-absorbed events,Ê * 0 is estimated by a good approximation as the sum of the measured energies: The true interaction positions are approximated as the measured positions: i.e., the positional response function is approximated as a delta function: However, with this approximation, the positional errors in Eq. 24 become ignored. To compensate for this, we modify the energy resolution in Eq. 23. At i-th interaction (2 ≤ i ≤ n−1), the scattering angle can be calculated geometrically, and the positional errors (σx, σŷ, and σẑ) produce the uncertainty on it, to which we refer as ∆(cos ϑ scat i ) pos . Considering the propagation of the positional errors as discussed in the classical approach [27], ∆(cos ϑ scat i ) pos can be calculated. Then, the partial derivative of Eq. 9 with respect to cos ϑ scat which means that the positional errors produce uncertainty on the estimation of the gamma-ray energy through the scattering angle calculation. In our algorithm, for 2 ≤ i ≤ n − 1, we modify the energy resolution as Thus, though the approximation by Eq. 26 ignores the positional errors described in Eq. 24, we include them in the energy response function effectively. The incoming direction (θ * 0 ,φ * 0 ) should be also determined. In Compton telescopes, the incoming gamma-ray direction is constrained on a circle in the sky. We calculate the first scattering angle under assumption of the incident energy (Ê * 0 ), the 9 deposited energy (ε τ (1) ), and the first and second interaction positions (r τ (1) , r τ (2) ), and then the circle can be obtained as the solutions of (θ * 0 ,φ * 0 ) that satisfy the Compton kinematics: HereÊ 1 is calculated by Eqs. 9 and 10. The function P path,scat (r 1 |q 0 ) also depends on θ * 0 andφ * 0 . Here we assume a constant value L 0 for the path length in Eq. 18: This should be comparable to the size scale of the detector. Then, the probabilities are the same for the parameters (θ * 0 ,φ * 0 ) that satisfy Eq. 30, i.e., on a Compton circle, and they are considered to be the approximation of the maximum probability for the assumed scattering order. It should be noted that this approximation ignores the dependence of l first on the first interaction position, i.e., detector geometry effect. We examined how the choice of L 0 affects the algorithm performance and found that it does not significantly depend on L 0 (see Appendix A).

Escape events
For escape events, the incident gamma-ray energy and escape energy can be estimated using the scattering angle measured from geometrical information [26]. Here we calculate the incident gammaray energy at each i-th interaction site (2 ≤ i ≤ n − 1), and use the averaged value as the estimation. Based on the three-Compton method [28,30], E * 0 is estimated aŝ where, In the calculation of the probability function for escape events, the term P esc q n , D τ (n) |q n−1 needs the integration over the direction of escape (see Eq. 15). To calculate it approximately without integral computation, here we assume that We also approximatel n in Eq. 21 as a constant value:l where L esc should be also comparable to the detector size scale like L 0 . Again this implementation ignores the detector geometry effect. We also found that the algorithm performance is not significantly affected by the choice of L esc in Appendix A. Then these approximations yield In the last equation, dσ scat dΩ is a fixed value because after the integration the scattering angle ϑ scat n is determined fromÊ n−1 and ε τ (n) with Compton kinematics. Note that to derive the last equation we used change of variables cosθ scat n →ε n using Eq. 9. The term mec 2 (Ên−1−ε τ (n) ) 2 is due to this change of variables. Here the direction of escape is described with the scattering angleθ scat n at the last interaction and the azimuth angleφ scat n along r τ (n) − r τ (n−1) , not θ n andφ n , because it makes the calculation simpler, e.g., d cos ϑ scat n /dε n .

Algorithm Summary
By calculating the probability functions as described in the previous section, the scattering order and the event type (fully-absorbed or escape) can be determined as follows.
1. One selects a scattering order from all candidates. Here the selected order is labeled with k: When the number of the hits is n, the number of the candidates is n!, i.e., 1 ≤ k ≤ n!. 2. One calculates the conditional probabilities given the event type and scattering order described in Eqs. 11 and 14, and uses the approximations introduced in Section 4.2. Then, for each scattering candidate, one obtains the two probabilities P k fullabs and P k escape , which correspond to the fully-absorbed and escape events, respectively. The calculation procedure of P k fullabs and P k escape is described in the following subsections. 3. One determines the scattering order and the event type as a set of them that yield the maximum value in all of the calculated P k fullabs and P k escape .
5.1. Calculation of P k fullabs 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. 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. The probability is calculated using Eq. 11 by applying the approximations described in Eqs. 29, 31 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. The first scattering angle is calculated by Eq. 30 and incoming gamma-ray direction is constrained on a Compton circle in the sky. 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.

Numerical Experiments
In order to test the reconstruction algorithm developed in this work, we apply it to a simulation data set generated by ComptonSoft, a software package of Geant4-based simulation and data analysis for Compton telescopes [48,49,50]. The GRAMS experiment utilizes a large volume detector using a liquid argon time projection chamber, and MeV gamma-ray events are considered to be dominated by multiple Compton scattering events [10]. We choose this type of Compton telescope for a numerical demonstration of the reconstruction algorithm, and evaluate the performance of the algorithm, i.e., the accuracy of the event classification (fully-absorbed or escape events), the prediction accuracy of the incident gamma-ray energy, the scattering order, and the angular resolution.

Simulation setup
We assume a detector with a size of 140×140×20 cm 3 filled with liquid argon as shown in Figure 4. It is the same as proposed in [10]. The energy resolution σε 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 ∼ 1% (1σ) 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 σx and σŷ to be 2/ √ 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 σẑ = 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 σẑ can be different from σx /ŷ . 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 µ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 first in Eq. 31 and L esc in Eq. 37 to be 10 and 20 cm respectively (see Appendix A for this optimization). The computational time for different numbers of hits is described in Table 1. For this table, we used 1 MeV gammaray events in the following subsection. 6.2. Results of event classification and energy reconstruction We simulated 10 8 events of 1 MeV gamma-ray beam incoming from the top of the detector, i.e., θ 0 = π. 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".  Note that 0 hit represents the event that gamma ray escapes from the detector without any interaction or all the produced signals are lower than the threshold.
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, con-firming 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.

Angular resolution and its 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  Figure 7. The orange lines represent the spectra of the total detected energy of events classified as escape events by the algorithm. The red ones represent the spectra of the incident energy of those events, estimated by the algorithm using Eq. 32. Note that the events for the orange and red lines are the same, but the energies shown here are different (one is the total energy deposit, and the other is the reconstructed energy).  where θ K is the estimated first scattering angle: and θ 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 ∼ 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 forwardscattering events. For example, when using events with scattering angles less than 60 degrees at the first interaction, the ARM FWHM of 3-hit fullyabsorbed events is improved from 7.1 degrees to 4.1 degrees. Figure 11 also shows the FWTM of the ARM distributions.

Accuracy of the event reconstruction
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 gammaray 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 "×" 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 fullyabsorbed 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 Appendix B.

Comparison with other algorithms
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 θ 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 figureof-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 fullyabsorbed 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 fullyabsorbed 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 smallnumber-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 • .

Performance in different gamma-ray energies
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  type + order (all interactions) Figure 13: The same as Figure 12, but using truly fully-absorbed events. type + order (all interactions) Figure 14: The same as Figure 12, but using escape events with truly deposited energies below 900 keV.  Figure 15: Comparison of the performance with those of other algorithms. The black lines are the same as Figure 12.
The green and blue lines represent the modified TANGO algorithm [42] and the classical method [29]. The middle panel shows the accuracy of the scatter order of the first two interactions (solid), and all interactions (dotted). The bottom one shows the accuracy for both event type and scattering order when focusing on the scattering order of the first two interactions for fully-absorbed events and the first three for escape events. 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 ∼ 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.

Discussion and Conclusion
In this work, we have formulated probability functions of physical processes and measurements in Compton telescopes and developed the reconstruction algorithm for the multiple Compton scattering events based on the probabilistic method. This algorithm can treat both fully-absorbed and escape events in a unified framework. The developed algorithm can be used in compact or largevolume Compton telescopes aiming at a few MeV.
A promising application is the GRAMS project which uses a large single-material detector with liquid argon. We verified its performance using simulated data sets of a 140 × 140 × 20 cm 3 liquid argon detector and confirmed that the algorithm works very well for up to 8-hit events for 1 MeV gamma rays. The accuracy for the event type is more than 73%, with a maximum of 95% at 8-hit events. While the scattering order of the entire sequence is reconstructed with an accuracy of about 55% at the maximum, the accuracy gets much higher, e.g., ∼80% for 8-hit events, when focusing on the essential part for the reconstruction, namely, the scattering order of the first two or three interactions. For the application to the GRAMS project, we consider using the events from 3 to 8 hits to achieve a large effective area which is one of its advantages. However, towards practical applications, we need further studies when considering the background, event selection, and some processes ignored in the proposed algorithm. In the following, we discuss several aspects for extending and improving the algorithm.

The background reduction and measurement of gamma rays with known energy
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 gammarays, 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 (P fullabs or P 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 positionrandomized 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 ∼ 30%.  Figure 19: The accuracy for both type and order, when the incident gamma-ray energy is known a priori (red). The scattering order of first two and three interactions are considered for fully-absorbed and escape events, respectively. The black lines are the same as Figure 12. The incoming gamma rays are 1 MeV.

Possible extension of the algorithm for electron-tracking Compton telescopes
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, D I is needed to be redefined as where p e,I is the measured momentum direction of a recoiled electron. Then, we modify the detector response function P det D I |q i−1 ,q i by including the electron trajectory information. Here we introduce a function P track (·) that compares the expected direction of a recoiled electron p i −p i−1 and measured one p e,I . The definition of P track (·) would depend on the configuration of Compton telescopes because the obtained electron images vary by the detectors, and in general P det (·) is reformulated as After modifying P det (·), 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.

For further improvement
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.

Effect of passive materials on the algorithm performance
The probabilistic model proposed in §3 assumes that incoming gamma rays always interact with a sensitive part of the detector, and interaction with passive materials is not considered. However, the sensitive volume is surrounded by passive materials inevitably in a realistic detector configuration. For example, the GRAMS project's liquid argon time projection chamber needs a cryostat. Insensitive regions of liquid argon, photomultiplier tubes or Si-PMs (scintillation light detectors), pixels or wires with readout electronics for ionized electron detection can also act as passive materials. These components may deteriorate the performance of the detector and hence that of the reconstruction algorithm since the energy deposits in them cannot be detected. Part of this deterioration may be recovered by considering the interactions with passive materials in the algorithm. Since the design of the GRAMS detector is still under discussion, here we summarize what kinds of scattering patterns occur with passive materials and discuss possible ways to consider them. Figure 20 shows examples of the scattering patterns ignored in the probabilistic model when passive materials are considered. In the pattern (A), an incoming gamma ray is scattered in a passive region as the first interaction. In this case, the detected hits are indistinguishable from an event with a gamma-ray incoming from the dotted arrow in Figure 20. Thus, it is difficult to identify and reconstruct an incident gamma ray by the algorithm, and this pattern should be treated as a background. In order to reduce these events, thin and low-density material is preferable for the top plate of the cryostat in the GRAMS detector. Hereafter once the first interaction occurs in the passive region, we label the event as (A) no matter whether it is scattered in the passive region after that.
The pattern (B) in Figure 20 shows another example that a gamma ray is scattered in a sensitive region at first but suffers from scattering by passive materials as an intermediate interaction. There are two ways to handle such an event. One possibility is to use the fact that these events may yield a conditional probability smaller than events with interactions occurring only in a sensitive volume. Although the current implementation considers these events as fully-absorbed or escape events assuming no scattering in passive materials, the scatter-ing angles calculated redundantly from kinematics and position information are hard to be matched due to the loss of information about the scattering in the passive region. It would result in that the conditional probability calculated by the algorithm becomes smaller. If so, these events could be distinguished and rejected by applying a threshold for the conditional probability. The other way is to include the scattering process with passive materials in the probabilistic model. While it may require complicated mathematical calculations, in principle, it is possible to include the scattering process in the model by marginalizing the scattering position and angle, similar to the marginalization introduced for the escape events in Eq. 38.
Finally, the pattern (C) in Figure 20 is a gammaray that is scattered in a sensitive region except for the last interaction. No matter whether the last interaction in the passive region is Compton scattering or photoabsorption, this pattern can be effectively regarded as an escape event considering only interactions in a sensitive volume. Thus, when the proposed algorithm works correctly, such an event is identified as an escape event, and the undetected gamma-ray energy can be estimated by Eq. 32.
To consider these events more, we performed a simulation with 1 MeV gamma rays from the zenith the same as §6.2, but the sensitive volume is surrounded by a cryostat with a thickness of 4 mm. It consists of stainless steal (SUS304). Table 2 shows the ratio of the scattering patterns obtained from the simulation for different numbers of hits. Only hits in the sensitive region are considered for calculating the number of hits, and the interaction in passive materials is not counted in it. As the number of hits gets small, the ratio of events suffered from scattering in passive materials becomes large. The dominant scattering patterns in this setup are (A) and (C). The former can be reduced by using a thinner cryostat with other material, e.g., aluminum. Also, we found that more than 80% of the pattern (C) events are classified as escape events by the algorithm, and the energy deposit in the passive region is corrected. Figure 21 shows the reconstructed energy spectrum for 3-hit events comparing with the ideal simulation shown in §6.2. We can see that the low energy tail appears mainly due to the scattering pattern (A). Note that here we consider only the cryostat, but in reality, other passive materials exist as listed before. Thus, further considerations are essential, especially after the design of the GRAMS detector is finalized.

Doppler broadening effect and bremsstrahlung
The proposed algorithm also ignores two important factor for physical processes. One is the momentum distribution of electrons in the detector materials. When the initial electron momentum is not zero, the scattered gamma-ray energy differs from Eq. 9. This effect is known as the Doppler broadening effect [45], and limits the angular resolution. Effectively, this effect could be considered by adding a function σ Doppler in Eq. 29, which is an uncertainty in the energy determination at each interaction site due to the Doppler broadening effect. Since it depends on the gamma-ray energy, the scattering angle, and detector materials, the function σ Doppler should be carefully modeled. When the Doppler broadening effect is comparable to or dominates over the position and energy resolutions, then such a modification would improve the reconstruction accuracy. Note that this effect is not significant in the case of 1 MeV gamma-ray simulation shown in § 6 [45].
The other factor is bremsstrahlung from scattered electrons in the detectors. Bremsstrahlung becomes the dominant energy loss at high energies. For example, in the argon detector, it dominates over the ionization losses at a few MeV. Gamma rays emitted by this process make additional signals, which makes the event reconstruction more complex. Decreasing the accuracy above a few MeV is considered to be partially due to bremsstrahlung (see Figure 17). This process is not straightforward to treat because, unlike Compton scattering, it is challenging to identify the site where a bremsstrahlung photon is produced. A possible way is to calculate the probability that a given signal is produced by a bremsstrahlung photon, considering energies and lengths to other interaction sites, and then merge it to corresponding signals if the calculated probability is large.
In a high energy band, typically above 5 MeV, it is also needed to distinguish pair creation events. To consider these factors (the Doppler broadening, bremsstrahlung, and pair creation), other statistical methods using large data sets might be effective, e.g., the deep neural network technique. We expect that combining an analytical method like this work and simulation-based statistical methods can be a promising way to achieve even better reconstruction performance. As a first step, we are also developing an event reconstruction algorithm using a multi-    task neural network. In a subsequent paper [58], we show that a neural network model improves the performance of the event reconstruction, especially for events with a small number of hits around 1 MeV.

Code availability
The code for the reconstruction algorithm is available at https://github.com/ odakahirokazu/ComptonSoft. This approximation is very simple, and we checked how it affects the algorithm's performance. Figure A.22 shows the accuracy for both type and order, with different values of L 0 from 0 cm to 100 cm. The reconstruction accuracy varies by only 3% at most. We also checked the effect of L esc on the performance. Figure A.23 shows the accuracy for both type and order, with L esc from 0 cm to 100 cm. In this case, the performance gets the best with L esc of 10 cm for 3-hit events and 20 cm for 4-hits events, but as long as an extreme value like 0 or 100 cm is not used, the performance does not depend on L esc so much. Thus we conclude that L 0 and L esc do not affect the performance so much, and this simple approximation is valid. In §6, we set L 0 to 10 cm, a half of the detector thickness, and L esc to 20 cm since it maximizes the performance at 4 hits. Here we examine how the reconstruction accuracy is affected by the assumed energy and positional resolutions. In Figure B.24, we show the accuracy when the energy resolution is improved by a factor of 2 or 4. In the latter case, the accuracy is increased by about 12% at most. Moreover mm. In the best pixel resolution, the accuracy is improved by about 5%.