Abstract:
Gamma-ray tracking methods for use with granular, position sensitive detectors identify the sequence of the interactions taking place in the detector and, hence, the position of the first interaction. The improved position resolution in finding the first interaction in the detection system determines a better definition of the direction of the gamma-ray photon, and hence, a superior source image resolution. A PET system using such a method will have increased efficiency and position resolution.

Description:
This application claims priority to U.S. Provisional Patent Application Ser. No. 60/755,470, titled “Gamma-Ray Tracking Method for Pet Systems” filed Dec. 29, 2005, and to U.S. Provisional Patent Application Ser. No. 60/755,469, titled: “Methods for Increasing Sensitivity of Gamma-Ray Imagers” filed Dec. 29, 2005, both incorporated herein by reference. 

   The United States Government has rights in this invention pursuant to Contract No. W-7405-ENG-48 between the United States Department of Energy and the University of California for the operation of Lawrence Livermore National Laboratory. 

   BACKGROUND OF THE INVENTION 
   1. Field of the Invention 
   The present invention relates to gamma ray detection, and more specifically, it relates to techniques for tracking gamma rays. 
   2. Description of Related Art 
   Positron Emission Tomography (PET) is a very important tool used for medical diagnostics and research. A PET system of increased spatial resolution and efficiency will be of great impact, particularly on small animal imaging, where the features that need to be identified are in the range of millimeters. 
   The spatial resolution in most of the state-of-the-art PET systems is limited by two factors: 1) the dimension of the scintillating pixels used in the detector arrays, and 2) the Compton scatterings that the gamma-ray undergoes in the detector array before being absorbed in a photoelectric interaction. The use of heavy scintillators helps reduce the fraction of Compton scatterings, and the average length between scatterings, but Compton scatterings are still important. 
   Position sensitive semiconductor detectors, especially the double sided segmented planar detectors, are alternatives to scintillator arrays. Such detectors have the advantage of a 3D position sensitivity, and potentially superior position resolution and granularity. The introduction of the planar Ge detectors with orthogonal strip segments has allowed the development of new, compact and efficient Compton scatter imagers. These detectors are characterized by a superior energy resolution. Position sensitivity is provided by the orthogonal strip segments of around 2 mm pitch. Despite a limited segment pitch size of 2 mm, position resolution in the sub-millimeter range is demonstrated by using new waveform filtering methods. Some of the methods employed to make that possible are described in U.S. patent application Ser. No. 11/607,554, titled: “Methods for Increasing the Sensitivity of Gamma-Ray Imagers,”) filed Nov. 30, 2006, using Express Mail label number EV424773615US, incorporated herein by reference. Position sensitive detectors normally used in PET systems provide good intrinsic position resolution, but do not correct for multiple interactions caused by Compton scatterings of the 511 keV annihilation photon in the detector. 
   SUMMARY OF THE INVENTION 
   It is an object of the present invention to provide a method for determining which photon interaction of said a plurality of photon interactions in a detector or detector segment of a PET system has the highest probability of being the first photon interaction in the detector or detector segment. 
   It is another object to reduce the blurring effects in the image produce by a PET system by reducing or eliminating photons scattered within the sample. 
   These and other objects will be apparent based on the disclosure herein. 
   A method to increase imaging sensitivity for PET scanners improves the imaging resolution by removing the blurring effects created by the photons interacting multiple times. In order for the methods to be effective, the PET scanners must have good position resolution and granularity, i.e., capability to identify several interactions taking place in the gamma-ray detector. Bayesian statistics is used to identify the point of the first interaction by determining the photon scattering sequence with the highest likelihood of occurrence. If at least one photon multiply scatters in the detector, it is possible to identify photons that have scattered in the sample, before interacting in the detection system. In this way, the image contrast is also improved. The presented methods work with fully absorbed photons as well as with scatter photons (from the Compton background). 
   The methods presented here make use of the data available in a high granularity, high position resolution detection system to determine more accurately the direction of the incident annihilation gamma-rays by finding the most likely scattering sequence of the gamma-ray photon in the detection system, when the gamma-ray interacts multiple times in the detector. This will make available for high resolution imaging photons that otherwise would have been rejected, or would have provided poor image quality. The position sensitive detectors can be either scintillators or semiconductor detectors. However, the detectors must be able to provide acceptable granularity of detection, i.e., they must be able to provide the positions and energies of multiple interactions taking place in the detection system. A finer detection granularity will increase the efficiency of the proposed event reconstruction methods. Previous methods for gamma-ray interaction sequencing were developed for nuclear physics experiments, and for Compton cameras. The methods presented here are optimized for PET scanners. 
   The invention includes a method that can be programmed onto a computer readable medium. An embodiment of the method comprises the steps of listing each position within a first volume of at least two photon interactions of a first plurality of photon interactions resulting from a first annihilation gamma-ray and determining which photon interaction of the at least two photon interactions has the highest probability of being the first photon interaction in the first volume. The first volume can comprise at least one detector and specifically may be a DSSD. Generally simultaneously with the above steps, the method further comprises listing the position within a second volume facing the sample of a second at least two photon interactions of a second plurality of photon interactions produced by a second annihilation gamma-ray emitted at about 180 degrees from the first annihilation gamma-ray, where the sample is located between the first volume and the second volume; and determining which photon interaction of the second at least two photon interactions has the highest probability of being the first photon interaction in the second volume. The highest probability of being the first photon interaction in the first volume is determined with Bayesian probability theory and where the highest probability of being the first photon interaction in the second volume is determined with Bayesian probability theory. When the first photon interaction in the first volume has an energy of less that 511 keV, the first photon interaction in the first volume is assigned a low reliability index and where when the first photon interaction in the second volume has an energy of less that 511 keV, the first photon interaction in the second volume is assigned a low reliability index. The step of listing each position of at least two photon interactions of a first plurality of photon interactions within a first volume comprises collecting a first set of gamma-ray photon interaction produced data from at least one electrode segment of a first double-sided segmented detector (DSSD) and calculating, using the first set of gamma-ray photon interaction produced data, a 3 dimensional position of the gamma-ray photon interaction in the first DSSD, and where the step of listing each position of at least two photon interactions of a second plurality of photon interactions within a second volume comprises collecting a second set of gamma-ray photon interaction produced data from at least one electrode segment of a second double-sided segmented detector (DSSD) and calculating, using the second set of gamma-ray photon interaction produced data, a 3 dimensional position of the gamma-ray photon interaction in the second DSSD. 
   An improved PET system using the proposed methods will be of useful for biomedical investigations of interest for DOE and DHS. Spin-off algorithms can be used in gamma-ray imaging devices, such as Compton cameras, coded apertures and SPECT systems. An improved small animal PET system of increased resolution and sensitivity will be of great interest in biomedical research. Research laboratories, universities and drug manufacturing companies are potential users. A large PET system using the proposed methods will be of interest for biomedical research, as well as for clinical investigations. 

   
     BRIEF DESCRIPTION OF THE DRAWINGS 
     The accompanying drawings, which are incorporated into and form a part of the disclosure, illustrate embodiments of the invention and, together with the description, serve to explain the principles of the invention. 
       FIG. 1  shows a system of two double sided segmented planar Ge detectors acting as a PET system. 
       FIG. 2  shows a system of eight double sided segmented planar Ge detectors in a small animal PET system. 
       FIG. 3  shows an example coincident event. 
       FIG. 4  shows an attenuation curve for a photon in the detector. 
       FIG. 5  shows an exemplary flow chart of analysis steps used for PET imaging using Bayesian statistics. 
   

   DETAILED DESCRIPTION OF THE INVENTION 
   A PET detection system is a coincidence system that aims to detect two annihilation gamma-rays emitted at about 180 degrees from each other by a β+ radioactive tracer located inside a sample. The positions of the first interactions are used to reconstruct the direction of the two annihilation gamma-rays. This direction forms a line of projection, sometimes referred to as a line-of-response (LOR). By acquiring multiple such events, the distribution of the tracer inside the sample is reconstructed by analyzing the intersection of the corresponding projection lines. 
   The methods presented here make use of the data available in a high granularity, high position resolution detection system to determine more accurately the position where the photon interacted first in the detection system. This is achieved by finding the most likely scattering sequence of the gamma-ray photon in the detection system, when the gamma-ray interacts multiple times in the detector. This is made possible by using a gamma-ray tracking algorithm as described below. The algorithm takes as input the list of interaction vertices provided by a position sensitive detector. Six values are associated to each vertex: a 3-dimensional position vector (r), the deposited energy (E), the index of the detector hit (D), and, if available, a factor characterizing the reliability of the identified interaction (ρ). The factor ρ indicates whether the identified vertex is a single, clean interaction or two overlapping interactions. Overlapping interactions are herein referred to as the interactions that take place within the same detection element, such as a segment For each event, the gamma-ray tracking algorithm will provide the vertex of the first Compton scattering and the likelihood associated with the identified vertex. 
     FIGS. 1 and 2  show examples of PET implementations that are suitable to use with the methods described here. These are systems of 2 and 8 Ge planar double sided segmented detectors (DSSDs), respectively.  FIG. 1  shows a sample  30  located between two Ge double sided segmented detectors  32 ,  34 . The radioactive tracer generates the two coincident annihilation gamma-rays  36 ,  38 .  FIG. 2  shows a system that includes eight (8) detectors  41 - 48 . DSSDs have the capability to provide good position resolution and granularity without compromising the energy resolution; however, any type of detector can work, as long as it provides good granularity, position and energy resolution. 
   Bayesian statistics are used to estimate the most likely scattering sequence for both coincident photons. According to the Bayes theorem, the conditional probability for a certain model, M i , given a set of data, X, is expressed as: 
                   P   ⁡     (       M   i     |   X     )       =           P   ⁡     (     M   i     )       ⁢     P   ⁡     (     X   |     M   i       )           P   ⁡     (   X   )         =         P   i         ∑   j     ⁢           ⁢     P   j         .               (   1   )               
P(M i ) is the prior probability known about the model M i . P(X|M i ) is the conditional probability for data X given model M i , and P(X) is the probability of occurrence of data X. If all possible models are known for a given set of data, P(X) can be calculated as the sum over all models:
 
                   P   ⁡     (   X   )       =       ∑   j     ⁢           ⁢       P   ⁡     (     M   j     )       ⁢       P   ⁡     (     X   |     M   j       )       .                 (   2   )               
This probability will provide the normalization factor for P(M i /X), which is important in assessing the reliability of the reconstructed event.
 
   In the present case, the model M i , is defined by a scattering sequence associated with each of the two coincident photons and a total energy for each incident photon of 511 keV. This is the energy of the (e+e−) annihilation gamma-ray. If no Compton scatterings take place in the sample, the measured event will be consistent with the model for which the incident gamma-ray is 511 keV. The calculation of a probability for this model will finally either confirm or reject the model. At the end, the model with the highest probability will determine the scattering sequence for the two photons, and whether the two incident photons have Compton scattered in the sample or not. 
   When the 511 keV photon Compton scatters in the sample, the photon impinging on the detector can have an energy E total  smaller than 511 keV. Using the tracking algorithm, it will be possible to estimate E total  by resolving the kinematics in the Compton scatterings at the first assumed interaction: 
   
     
       
         
           
             
               
                 
                   E 
                   
                     total 
                     , 
                     k 
                   
                 
                 = 
                 
                   
                     
                       
                         E 
                         
                           int 
                           , 
                           k 
                         
                       
                       2 
                     
                     ⁢ 
                     
                       ( 
                       
                         1 
                         + 
                         
                           
                             1 
                             + 
                             
                               
                                 4 
                                 ⁢ 
                                 
                                   m 
                                   0 
                                 
                                 ⁢ 
                                 
                                   c 
                                   2 
                                 
                               
                               
                                 
                                   E 
                                   
                                     int 
                                     , 
                                     k 
                                   
                                 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     1 
                                     - 
                                     
                                       cos 
                                       ⁢ 
                                       
                                           
                                       
                                       ⁢ 
                                       
                                         θ 
                                         k 
                                       
                                     
                                   
                                   ) 
                                 
                               
                             
                           
                         
                       
                       ) 
                     
                   
                   .. 
                 
               
             
             
               
                 ( 
                 3 
                 ) 
               
             
           
         
       
     
   
   Here, θ k  is the measured gammed-ray scattering angle for the assumed sequence for photon k (k=1,2), m 0 , is the rest mass of the electron. The scattering angles are determined by the direction made by the first assumed interaction points for the two photons, and the direction of the scattering photon k (direction made by interaction points  11  and  12  for photon k=1, and interaction points  21  and  22  for photon k=2 in  FIG. 3 ). The calculation of the probability of a model given the observed interactions uses the probability density function for the calculated energies E total,k  as determined by the propagation of the uncertainties of the parameters in the right side of Equation 3. Assuming a normal distribution for the spread of the total energy as calculated using the Compton scattering formula, E total,k   (c) , around the modeled E total,k , the P(X|M i ) term becomes: 
                   P   ⁡     (     X   |     M   i       )       =       1       σ     E   1       ⁢     σ     E   2       ⁢   2   ⁢           ⁢   π       ⁢       e     -     (           (         m   e     ⁢     c   2       -     E     total   ,   1         )     2       σ     E   1     2       +         (         m   e     ⁢     c   2       -     E     total   ,   2         )     2       σ     E   2     2         )         .               (   4   )               
If the spread of the total energy is assumed to be a Lorentz distribution of width σ Ek , the term P(X|M i,k ) is calculated by integrating the probability density function ƒ(E):
 
                         P   ⁡     (     X   |     M     i   ,   k         )       =         ∫     E     total   ,     k       -     σ   k       /   2               E     total   ,   k       +       σ   k     /   2         ⁢       f   ⁡     (   E   )       ⁢           ⁢     ⅆ   E         =       ∫     E     total   ,     k       -     σ   k       /   2               E     total   ,   k       +       σ   k     /   2         ⁢         σ   k       2   ⁢     π   [         (         m   e     ⁢     c   2       -   E     )     2     +       (       σ   k     2     )     2       ]         ⁢           ⁢     ⅆ   E     ⁢         .                     (   5   )               
Then, P(X|M i ) is the product of the two probabilities P(X|M i,k ), with k=1,2:
 
   
     
       
         
           
             
               
                 
                   P 
                   ⁡ 
                   
                     ( 
                     
                       X 
                       | 
                       
                         M 
                         i 
                       
                     
                     ) 
                   
                 
                 = 
                 
                   
                     1 
                     
                       π 
                       2 
                     
                   
                   ⁢ 
                   
                     atan 
                     ( 
                     
                       1 
                       
                         2 
                         ⁢ 
                         
                           
                             
                               ( 
                               
                                 
                                   
                                     m 
                                     e 
                                   
                                   ⁢ 
                                   
                                     c 
                                     2 
                                   
                                 
                                 - 
                                 
                                   E 
                                   
                                     total 
                                     , 
                                     1 
                                   
                                 
                               
                               ) 
                             
                             2 
                           
                           
                             σ 
                             
                               E 
                               1 
                             
                             2 
                           
                         
                       
                     
                     ) 
                   
                   ⁢ 
                   
                     
                       atan 
                       ( 
                       
                         1 
                         
                           2 
                           ⁢ 
                           
                             
                               
                                 ( 
                                 
                                   
                                     
                                       m 
                                       e 
                                     
                                     ⁢ 
                                     
                                       c 
                                       2 
                                     
                                   
                                   - 
                                   
                                     E 
                                     
                                       total 
                                       , 
                                       2 
                                     
                                   
                                 
                                 ) 
                               
                               2 
                             
                             
                               σ 
                               
                                 E 
                                 2 
                               
                               2 
                             
                           
                         
                       
                       ) 
                     
                     . 
                   
                 
               
             
             
               
                 ( 
                 6 
                 ) 
               
             
           
         
       
     
   
   The advantage of a Lorentz distribution is that its longer tails better approximate the spread determined by the Compton profile of bounded electrons. 
   Other competing models account for sequences that contain vertices of multiple overlapping interactions. In the present approach, the events that are likely to be described by such models are rejected. 
   Once the model is defined, P(M i ) can be calculated by accounting for the individual probabilities of the gamma-ray to go L j-1,j  distance through the detector between two successive interactions without being absorbed, and then, interacting within an arbitrary interval δ either in a Compton scattering, or in a photoelectric effect if the model assumes that the current interaction is the last in the event, i.e., a photoelectric absorption. 
                   ⁢             P   ⁡     (     M   i     )       =       ∏     k   =   1     2     ⁢           ⁢       ∏     j   =   1     N     ⁢           ⁢       P     L   ,   j   ,   k       ⁢       P     C   ,     P   ⁢           ⁢   h     ,   j   ,   k       .                   (   7   )                 
The individual probabilities P L,j,k  become:
 
                         P     L   ,   j   ,   k       =       ⁢       ∫     L       (       j   -   1     ,   j     )     ,   k         L       (       j   -   1     ,   j     )     ,     k     +   δ             ⁢       μ     total   ,   j   ,   k       ⁢     ⅇ     -     μ     total   ,   j   ,     k   x             ⁢           ⁢     ⅆ   x                     =       ⁢       (     1   -     ⅇ       μ     total   ,   j   ,   k       δ         )     ⁢     ⅇ     -     μ     total   ,   j   ,       k   L     ⁡     (       j   -   1     ,   j     )       ,   k                         ≅       ⁢       μ     total   ,   j   ,   k       ⁢     δ   ·     ⅇ     -     μ     total   ,   j   ,       k   L     ⁡     (       j   -   1     ,   j     )       ,   k                           =       ⁢       P     interaction   ,   k       ·     P     propagation   ,   k                       (   8   )               
Whereas the probability for Compton scattering is determined by the factor:
 
                     P     C   ,   j   ,   k       =         μ     KN   ,   j   ,   k         μ     total   ,   j   ,   k         ⁢       μ       KN   ⁡     (   θ   )       ,   j   ,   k         μ     KN   ,   j   ,   k             ,   .           (   9   )               
The probability for a photoelectric interaction is:
 
                   P       P   ⁢           ⁢   h     ,   j   ,   k       =         μ     Ph   ,   j   ,   k         μ     total   ,   j   ,   k         ..             (   10   )               
μ total,j,k  is the total attenuation coefficient (μ total,j,k =μ Ph,j,k +μ KN,j,k ), μ Ph,j,k  is the photoelectric absorption coefficient, μ KN,j,k  is the total Compton absorption coefficient, and μ KN(θ),j,k  is the equivalent absorption coefficient for a Compton scattering at angle θ.  FIG. 4  shows an attenuation curve for a photon in the detector. μ KN(θ),j,k  is calculated using the Klein-Nishina differential cross section for the Compton scattering:
 
                                 μ       K   ⁢           ⁢     N   ⁡     (   θ   )         ,   j   ,   k       ∝         ∂   σ         ∂   cos     ⁢           ⁢   θ       ⁢   δcos   ⁢           ⁢   θ       =                 π   ⁢           ⁢     r   0     ⁢         E     scat   ,   j   ,   k         E     0   ,   j   ,   k         ⁡     [         E     0   ,   j   ,   k         E     scat   ,   j   ,   k         +       E     scat   ,   j   ,   k         E     0   ,   j   ,   k         -       sin   2     ⁢   θ       ]       ⁢   δcosθ   ⁢         .                       (   11   )               
The total Compton scattering coefficient is:
 
   
     
       
         
           
             
               
                 
                   μ 
                   
                     
                       K 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       N 
                     
                     , 
                     j 
                     , 
                     k 
                   
                 
                 ∝ 
                 
                   
                     ∫ 
                     
                       - 
                       1 
                     
                     1 
                   
                   ⁢ 
                   
                     
                       
                         ∂ 
                         σ 
                       
                       
                         
                           ∂ 
                           cos 
                         
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         θ 
                       
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       ⅆ 
                       cos 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     
                       θ 
                       . 
                     
                   
                 
               
             
             
               
                 ( 
                 12 
                 ) 
               
             
           
         
       
     
   
   Finally, introducing P(M 1 ) and P(X|M i ) in P i , the likelihood of the model becomes: 
   
     
       
         
           
             
               
                 
                   
                     
                       
                         P 
                         i 
                       
                       = 
                       
                         
                           P 
                           interactions_i 
                         
                         · 
                         
                           P 
                           propagation_i 
                         
                         · 
                         
                           P 
                           kinematics_i 
                         
                       
                     
                   
                 
                 
                   
                     
                       
                         
                           P 
                           interactions_i 
                         
                         = 
                         
                           
                             ∏ 
                             
                               k 
                               = 
                               1 
                             
                             2 
                           
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             
                               ∏ 
                               inter_j 
                             
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               
                                 ( 
                                 
                                   1 
                                   - 
                                   
                                     e 
                                     
                                       - 
                                       
                                         k 
                                         
                                           total 
                                           , 
                                           j 
                                           , 
                                           
                                             k 
                                             δ 
                                           
                                         
                                       
                                     
                                   
                                 
                                 ) 
                               
                               ⁢ 
                               
                                 
                                   μ 
                                   
                                     
                                       K 
                                       ⁢ 
                                       
                                           
                                       
                                       ⁢ 
                                       N 
                                     
                                     , 
                                     j 
                                     , 
                                     k 
                                   
                                 
                                 
                                   μ 
                                   
                                     total 
                                     , 
                                     j 
                                     , 
                                     k 
                                   
                                 
                               
                               ⁢ 
                               
                                 
                                   μ 
                                   
                                     
                                       K 
                                       ⁢ 
                                       
                                           
                                       
                                       ⁢ 
                                       
                                         N 
                                         ⁡ 
                                         
                                           ( 
                                           θ 
                                           ) 
                                         
                                       
                                     
                                     , 
                                     j 
                                     , 
                                     k 
                                   
                                 
                                 
                                   μ 
                                   
                                     
                                       K 
                                       ⁢ 
                                       
                                           
                                       
                                       ⁢ 
                                       N 
                                     
                                     , 
                                     j 
                                     , 
                                     k 
                                   
                                 
                               
                             
                           
                         
                       
                       , 
                     
                   
                 
                 
                   
                     
                       
                         
                           P 
                           propagation_i 
                         
                         = 
                         
                           e 
                           
                             
                               - 
                               
                                 
                                   ∑ 
                                   
                                     k 
                                     = 
                                     1 
                                   
                                   2 
                                 
                                 ⁢ 
                                 
                                   
                                     ∑ 
                                     inter_i 
                                   
                                   ⁢ 
                                   
                                       
                                   
                                   ⁢ 
                                   
                                     k 
                                     
                                       total 
                                       , 
                                       j 
                                       , 
                                       
                                         
                                           k 
                                           L 
                                         
                                         ⁡ 
                                         
                                           ( 
                                           
                                             
                                               j 
                                               - 
                                               1 
                                             
                                             , 
                                             j 
                                           
                                           ) 
                                         
                                       
                                       , 
                                       k 
                                     
                                   
                                 
                               
                             
                             ⁢ 
                             
                                 
                             
                           
                         
                       
                       , 
                     
                   
                 
                 
                   
                     
                       
                         P 
                         kinematics_i 
                       
                       = 
                       
                         
                           1 
                           
                             
                               ∏ 
                               
                                 k 
                                 = 
                                 1 
                               
                               2 
                             
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             
                               
                                 σ 
                                 
                                   E 
                                   
                                     j 
                                     , 
                                     k 
                                   
                                 
                               
                               ⁢ 
                               
                                 
                                   2 
                                   ⁢ 
                                   π 
                                 
                               
                             
                           
                         
                         ⁢ 
                         
                           e 
                           
                             
                               - 
                               
                                 
                                   ∑ 
                                   
                                     k 
                                     = 
                                     1 
                                   
                                   2 
                                 
                                 ⁢ 
                                 
                                   
                                     ∑ 
                                     
                                       interm_inter 
                                       ⁢ 
                                       _i 
                                     
                                   
                                   ⁢ 
                                   
                                     
                                       
                                         
                                           ( 
                                           
                                             
                                               
                                                 
                                                   E 
                                                   
                                                     total 
                                                     , 
                                                     j 
                                                     , 
                                                     k 
                                                   
                                                 
                                               
                                             
                                             
                                               
                                                 
                                                   - 
                                                   
                                                     E 
                                                     
                                                       total 
                                                       , 
                                                       j 
                                                       , 
                                                       
                                                         k 
                                                         
                                                           ( 
                                                           c 
                                                           ) 
                                                         
                                                       
                                                     
                                                   
                                                 
                                               
                                             
                                           
                                           ) 
                                         
                                         2 
                                       
                                       
                                         σ 
                                         
                                           
                                             E 
                                             ⁢ 
                                             
                                                 
                                             
                                             ⁢ 
                                             j 
                                           
                                           , 
                                           k 
                                         
                                         2 
                                       
                                     
                                     . 
                                   
                                 
                               
                             
                             , 
                           
                         
                       
                     
                   
                 
               
             
             
               
                 ( 
                 13 
                 ) 
               
             
           
         
       
     
   
   Monte Carlo simulated events were used for testing the described algorithms. A large, double sided segmented Ge detector was assumed as detection system. 
   The results show that the proposed tracking algorithm succeeds to reliably make available for imaging the events that consist of multiple interactions. The simulation results indicate that 60% more events can be made available for imaging. 
   Thus, the algorithm, as generally depicted in  FIG. 5 , generates a list of the interactions ( 50 ) including their energies and positions, checks to see if the total energy deposited is less than or equal to 1,022 keV ( 52 ) and if it is equal to 1,022 keV, then the probabilities are calculated for all possible combinations ( 54 ). The best probability is then calculated ( 56 ) and a determination is made whether the interaction is to be considered a viable event ( 58 ) or a scattered event ( 60 ). The viable events are then usable for imaging. 
   A digital computer system can be programmed to perform the method of this invention. Once programmed to perform particular functions pursuant to instructions from program software that implements the method of this invention, such digital computer system in effect becomes a special-purpose computer particular to the method of this invention. The techniques necessary for this are well-known to those skilled in the art of computer systems. 
   Computer programs implementing the method of this invention will commonly be distributed to users on a distribution medium such as floppy disk or CD-ROM. From there, they will often be copied to a hard disk or a similar intermediate storage medium. When the programs are to be run, they will be loaded either from their distribution medium or their intermediate storage medium into the execution memory of the computer, configuring the computer to act in accordance with the method of this invention. All these operations are well-known to those skilled in the art of computer systems. 
   The term “computer-readable medium” encompasses distribution media, intermediate storage media, execution memory of a computer, and any other medium or device capable of storing for later reading by a computer a computer program implementing the method of this invention. 
   The foregoing description of the invention has been presented for purposes of illustration and description and is not intended to be exhaustive or to limit the invention to the precise form disclosed. Many modifications and variations are possible in light of the above teaching. The embodiments disclosed were meant only to explain the principles of the invention and its practical application to thereby enable others skilled in the art to best use the invention in various embodiments and with various modifications suited to the particular use contemplated. The scope of the invention is to be defined by the following claims.