Position Reconstruction in a Dual Phase Xenon Scintillation Detector
We studied the application of statistical reconstruction algorithms, namely maximum likelihood and least squares methods, to the problem of event reconstruction in a dual phase liquid xenon detector. An iterative method was developed for in-situ reconstruction of the PMT light response functions from calibration data taken with an uncollimated -ray source. Using the techniques described, the performance of the ZEPLIN-III dark matter detector was studied for 122 keV -rays. For the inner part of the detector (100 mm), spatial resolutions of 13 mm and 1.6 mm FWHM were measured in the horizontal plane for primary and secondary scintillation, respectively. An energy resolution of 8.1% FWHM was achieved at that energy. The possibility of using this technique for improving performance and reducing cost of scintillation cameras for medical applications is currently under study.
Anumber of applications require measurement of the interaction coordinates within a particle detector. In the low energy region 1 MeV, these include medical radionuclide imaging, gamma-ray astronomy and direct dark matter search experiments. In the latter instance, which motivated the present work, event localization per se is not relevant for detection of dark matter particles, but position sensitivity is important for efficient reduction of the radiation background and correct identification of the candidate events.
ZEPLIN-III is a dual phase (liquid/gas) xenon detector built to identify and measure galactic dark matter in the form of Weakly Interacting Massive Particles (WIMPs). It operated at the Boulby mine (UK) between 2006 and 2011. The detector measures both the scintillation light (S1) and the ionisation charge generated in the liquid by interacting particles and radiation. The ionisation charge drifts upwards to the liquid surface by means of a strong electric field and is extracted into a thin layer of gaseous xenon where it generates UV photons by electroluminescence (S2). Both the scintillation and electroluminescence light are measured by a PMT array and the ratio between S1 and S2 allows to discriminate nuclear recoils (expected to be produced by elastic scatter of WIMPs off xenon nuclei) from the electron recoils from and -ray backgrounds. The details on liquid xenon detector technology as well as on operation of dual phase detectors can be found in recent review papers [1, 2].
The self-shielding property of liquid xenon reduces the rate of background in the interior of the liquid. Using accurate position reconstruction to select only events in an inner "fiducial" volume therefore improves sensitivity to the WIMP signal. While the depth of the interaction can be inferred very accurately (few tens of m FWHM) from the electron drift time in the liquid (the delay between S1 and S2), the position in the horizontal plane has to be reconstructed from the light distribution pattern across the PMT array. Another reason for analysis of the light distribution is the need to eliminate the multiple scatter events that can mimic the WIMP interactions if one of the scatters has occurred in a dead volume of liquid xenon from where no charge can be extracted.
The active volume of ZEPLIN-III is a flat layer of liquid xenon (40 cm in diameter and 3.6 cm thick) above a compact hexagonal array of 31 2-inch vacuum ultraviolet-sensitive PMTs (ETL D730/9829Q) immersed directly in the liquid . Such a flat geometry makes it (from the point of view of position reconstruction) rather similar to the well-studied scintillation camera, which is widely used in areas as diverse as medical research and experimental astrophysics [4, 5]. The position of an event in a scintillation camera is traditionally found by the Anger method which consists in calculating a centroid of the PMT response .
Statistical reconstruction algorithms by maximum likelihood and weighted least squares methods have gained popularity following the pioneering work of Gray and Macovski in 1976 . They offer better precision along with the possibility of checking if the input data correspond to a valid event. These methods require knowledge of the light response functions (LRF) that characterise the response of a given PMT as a function of position of an isotropic light source inside the sensitive volume of the detector. Typically, the LRFs are either measured directly (e.g. by means of a moving collimated radioactive source) or calculated from the detector geometry, either analytically or by means of a Monte Carlo simulation.
In the present work, a method of reconstructing LRFs in situ from the calibration data obtained by irradiating the detector by -rays from an uncollimated radioactive source was developed. Based on the set of reconstructed LRFs, the positions and light yields of scintillation events in the detector can be readily found using either maximum likelihood or weighted least squares methods. This procedure was applied to the WIMP-search data taken with ZEPLIN-III [8, 9, 10, 11].
Ii Experimental Setup
The target region of the detector is shown in Fig. 1. The electric field in the active xenon volume (3.9 kV/cm in the liquid and 7.8 kV/cm in the gas) is defined by a cathode wire grid 36 mm below the liquid surface and an anode plate in the gas phase, 4 mm above the liquid. A second wire grid is located 5 mm below the cathode grid just above the PMT array. This grid defines a reverse field region which suppresses the collection of ionisation charge for events just above the array and helps to isolate the PMT input optics from the external high electric field. The PMTs are powered by a common high voltage supply, with the outputs roughly equalised by means of attenuators (Phillips Scientific 804). The PMT signals are digitised at 2 ns sampling by 8-bit flash ADC (ACQIRIS DC265). To expand the dynamic range of the system, each PMT signal is recorded by two separate ADC channels: one directly and one after amplification by a factor of 10 by fast amplifiers (Phillips Scientific 770). The acquired waveforms were analysed by a dedicated software that searched for pulses above a certain threshold and stored them in a parametrised form . Subsequently, an event filtering tool was used to retain events with a fast S1 signal preceding a wider S2 one. All multiple scatter events containing more than one S2 are filtered out.
A Co radioactive source was used for calibrating the energy response of the detector. This source emits 122 keV and 136 keV -rays which are rapidly absorbed in liquid xenon (with attenuation length 4 mm for these energies) mostly by photoelectric capture . Consequently, most of the interactions can be considered point-like with full energy deposit. The source was positioned at approximately 190 mm above the liquid surface and as close as possible to the detector axis. The calibration was performed daily to monitor the detector stability. There were also several dedicated runs aimed at acquiring sufficient data to train the positioning algorithms. Before the second science run, a specially-designed rectangular copper grid was placed inside the chamber, above the sensitive volume (Fig. 2). The grid structure is 386 mm in diameter, and was manufactured by diamond wire cutting from a 5.1 mm thick copper plate; the void pitch is 30 mm and the straight sections are 5 mm wide. The thickness of the grid was chosen such that it would attenuate the -ray flux from the calibration source by approximately a factor of 2, creating a shadow image that can be used to verify and fine-tune position reconstruction.
Iii Event reconstruction methods
The problem of event reconstruction consists in finding the energy (or, rather, the light signal intensity ) and the position of an event given a set of the corresponding PMT output signals . For an event at position producing photons the probability of the -th PMT detecting photons is well approximated by the Poisson distribution :
where is the expectation for a number of photons detected by the -th PMT out of initial ones with being LRFs – the fraction of the photons emitted by a light source at position that produce a detectable signal in the -th PMT. The corresponding output signal in the general case is a random variable with an expectation value proportional to . The probability distribution for depends on the single photoelectron response of the corresponding PMT and can be quite complex [15, 16]. However, in a few special cases it can be approximated by simple functions. These special cases include:
Photon counting. If is small (say, less than 10) and the PMT has a narrow single photoelectron distribution then can be calculated (almost) unambiguously from by rounding the ratio , where is the average single photoelectron response of the PMT.
Normal distribution. If is large (say, 25 or more) and the single photoelectron distribution of the PMT is reasonably symmetric then, following from the central limit theorem, is approximately normally distributed with the mean equal to .
Iii-a Centroid and corrected centroid
The centroid method of position estimation is the oldest method used by Anger in the first gamma camera in 1957 . It is still widely in use due to its simplicity and robustness. The position estimate is found as the weighted average of PMT coordinates with weights determined by the light distribution across the PMT array:
where are the coordinates of the axis of i-th PMT, is the measured charge and is a flat-fielding coefficient which compensates for variations in gain and quantum efficiency across the PMT array. As one can see from equations (2), no information on LRFs and probability distribution is necessary for application of this method. On the other hand, while the centroid method works reasonably well close to the centre of the detector (up to 100 mm from the centre in ZEPLIN-III), it becomes increasingly biased for events in the periphery. Another disadvantage is that it gives no indication regarding the match of the actual light distribution to the expected one.
If there exists one-to-one mapping between the true position and the one reconstructed by the centroid method then it is possible to invert this mapping to obtain the unbiased "corrected" estimate from the biased centroid one. In practice, this is often done by building a look-up table for a number of known positions on a rectangular grid and interpolating between these points. Another possibility is to use Monte Carlo simulation to calculate the forward mapping and then to use numerical methods to invert it. The latter method was employed in the ZEPLIN-III event filtering routine. It was also used to obtain the first approximation in the iterative LRF reconstruction procedure.
Iii-B Maximum likelihood
The maximum likelihood (ML) technique [4, 5, 17] consists in finding the set of parameters that maximises the likelihood of obtaining the experimentally measured outcome. For the case of photon counting when are known for each PMT, the likelihood function can be easily calculated from the Poisson distribution (1):
Taking into account that , one can write 
where does not depend on neither or . If the LRFs are known, the best estimates and can be found in a straightforward way by maximising function (4). The best estimate of at given , can be found analytically:
By substituting for into (4) one obtains , which is a function of the position only. Then and are found by maximising either analytically or by numerical methods. As a bonus, for the 2D case can be visualised as a colour map, which is very useful for either debugging or checking the validity of a given event.
Iii-C Weighted least squares
If can be considered normally distributed the more flexible weighted least squares (WLS) method can be used instead of ML . In this case the parameter estimates are found by minimising the weighted sum of squared residuals :
where is the expected PMT output charge and is the weighting factor which is reciprocal to the variance of . The best estimates and are obtained by finding the global minimum of
The and minimisations can be separated, as in the likelihood case, reducing by one the dimensionality of the problem.
Under an assumption that is measured exactly and the variance of is only due to statistical fluctuations in the number of detected photoelectrons, the WLS method becomes equivalent to ML . However, in a real detector the measured differs from the true value because of electronic noise. The variance of is also typically higher then expected from Poisson statistics due to non-zero width of the single photoelectron distribution. Compared to the ML method, the WLS makes it much easier to account for these and other factors. Most importantly, it makes it possible to reduce the weights for those PMTs with less well known light response.
Iii-D Method choice
The choice of the WLS method for S2 reconstruction is straightforward: due to its high light output, the S2 signal statistic is quasi-normal except for the PMTs far from low energy events. These PMTs may be either ignored or clustered together so that the photoelectron statistic per cluster is quasi-normal too.
In the case of S1, the total collected charge (from the whole PMT array) is equivalent, depending on the event position, to 1–2 photoelectrons per keV; this means that in the region of interest for WIMP searches (<50 keV) the S1 distribution is too far from normal to use the WLS method with confidence. Consequently, the ML method was used.
Iv Reconstruction of light response functions
The ML and WLS methods described above rely on the knowledge of the LRFs . There are several methods for obtaining the LRFs described in the literature. The most straightforward of these is the direct measurement, by scanning the detector with a moving well-collimated -ray source [14, 20]. Unfortunately, a combination of several factors made this method impractical for the ZEPLIN-III detector. Because of the cryostat, the source could not be placed closer than 190 mm to the liquid surface, and therefore a long collimator was required to achieve good position resolution. However, the available space above the detector was extremely limited and installation of a scanning system would require reduction in the shield thickness that was deemed unacceptable.
Alternatively, the LRFs can be calculated from the detector geometry by means of a Monte Carlo simulation tuned to reproduce the experimental data . Although we attempted a similar method with ZEPLIN-III , the difficulty in reproducing the exact shape of the LRFs by Monte Carlo and time constraints limited its use.
Yet another approach, explored in , is to choose a suitable parametric model for the LRF and adjust the parameters so as to minimize the mean sum of squared residuals in the WLS method for a population of calibration events. This method works well for low parameter count models proposed in that study. However, the LRF shape for S2 in the ZEPLIN-III detector proved to be more complex, requiring a model with at least five parameters in the simplest case:
where is the distance from the PMT axis and , , , and are adjustable parameters. As a result, the parameter adjustment becomes much more difficult due to increased dimensionality of the problem. To overcome this difficulty we developed an iterative method of LRF reconstruction described below.
Iv-a Method description
In this method the detector is irradiated by a non-collimated monoenergetic gamma source and the PMT responses are recorded event by event. Even if the gamma source is not collimated, it is still possible to obtain an estimate for each event position using the centroid or the corrected centroid method, at least for the central part of the detector. After a sufficiently large event sample is acquired, making an additional assumption that the LRF depends smoothly on and assuming that all the events produce the same amount of light, one can obtain the first approximation for the LRF by fitting the PMT response to the events at different by a smooth function of .
This first approximation can now be used to obtain better estimates for the positions of the events in the sample using ML or WLS method. Compared to the centroid estimates, these new estimates are less biased, especially in the case of peripheral events. Fitting again the PMT response as a function of coordinates using the updated event positions gives a second approximation .
The above steps are repeated until some convergence criterion is reached. This can be the fact that the reconstructed dataset has attained some quality that the physical calibration events are known to possess, for example monoenergeticity or some known distribution in the plane. Another option is to iterate until the change in the LRFs on the next step falls below a pre-defined tolerance.
Some additional regularization may be necessary to force the iteration to converge. One is the choice of a smoothing function. Another is the use of some a priori known property of the LRF; for example in the case of a PMT with a circular photocathode it is reasonable to assume that the LRF has axial symmetry , where is the distance from the PMT axis. This type of regularization was used in LRF reconstruction for ZEPLIN-III, the applicability of it will be discussed in section IV-C.
Iv-B ZEPLIN-III example
In order to collect the data necessary for reconstruction of the S2 LRFs, the detector was irradiated with -rays from a Co source. The top plot in Fig. 3(a) shows the - distribution of the estimated Co event positions obtained with a corrected centroid algorithm. Clearly, the events on the periphery tend to be misplaced closer to the centre of the PMT array. The situation deteriorates in the bottom-right corner where one of the PMTs was not functioning. However, for the central part of the array, approximately up to 100 mm from the centre, the centroid performance is good enough to be used for reconstructing the first approximation for the LRFs. This is demonstrated in the bottom plot of Fig. 3(a), where the area of PMT response is plotted versus the distance from its axis, calculated from the event position estimated by centroid. The resulting scatter plot was fitted using linear least squares technique with a cubic spline (the smooth curve on the plot) which was used as a first approximation for the LRF for a given PMT. Then the set of LRFs obtained in this way was used to re-calculate positions of the -ray interactions using the WLS method, producing the position distribution shown on the top plot of Fig. 3(b), and the cycle was repeated. After 5 iterations, the LRFs converged to the final shape shown in Fig. 3(c).
As one can see, the final distribution of the estimated event positions clearly shows the projected image of the copper grid with no significant distortions even in the region close to the non-functioning PMT. Note the ring of events in the periphery of the detector. The analysis indicates that they are well reconstructed as the sum of squared residuals is compatible with that for the events from the main population. Our interpretation is that these events occured near the edge of the field cage where non-uniform electric field pushed the extracted charge even further to the periphery. While the centroid algorithm fails to separate them from the main population (they are actually reconstructed closer to the center than some other peripheral events), the WLS method allows to unambiguously identify them.
The important advantage of the method described above is its ability to handle many more parameters than it was required by the original five-parameter model (8). This means that one can use (together with appropriate regularization) much more flexible non-parametric LRF representations such as look-up table  or cubic spline  previously used only in conjunction with the direct scan method. The cubic spline has an additional advantage of being a smooth function and, as we have found, with appropriate choice of knots it does not require any additional regularization. For this reason, a cubic spline representation for axially symmetric LRFs was adopted for both S1 and S2 LRFs. The knot placement was adjusted experimentally to cover the region of most rapid change in the response function with a denser grid.
Naturally, the assumption about axial symmetry of the PMT response is only an approximation. Several factors, most notably non-uniformities of the PMT photocathode and spatial dependence of the light collection efficiency (especially near the detector edge) can produce considerable deviations from symmetry. Such deviations from the model lead to systematic errors in estimated event position and energy. Fortunately, for such poorly reconstructed events the minimized sum of squared residuals tends to be above average. Plotting against and for a calibration dataset reveals the areas of the detector where the actual light response does not conform to the model (or rather the model is not good enough). Examining such plots for ZEPLIN-III we have found no increase in value near the detector edge which means that the light collection efficiency is indeed axially symmetric for both inner and outer PMTs. This is explained by poor reflectivity of the detector construction materials to xenon scintillation light (175nm) and confirmed by Monte Carlo simulations.
On the other hand, LRF of several PMTs have shown deviation from axial symmetry at small , most probably due to photocathode non-uniformity. This was mitigated by introducing uncertainty for the corresponding LRF that effectively reduced the weight function in the sum (7) (in other words, the contribution of i-th PMT) in the regions where its response was less symmetric. To improve convergence of first iterations, such “bad” PMTs can be temporary ignored by setting .
The method described above assumes that every accepted calibration event produces the same amount of scintillation light, independently of its position in the detector. In fact, the number of scintillation photons per event is affected by systematic and statistic fluctuations. In the case of a well-designed dual phase detector, the systematic fluctuations are negligible for the following reasons:
the liquid scintillator stays uniform due to convection flow and diffusion;
the light yield for S2 depends on the field strength, the gas pressure (both uniform across the detector sensitive volume) and the gas gap width, very well controlled by measuring duration of S2 pulses;
almost all small-angle scatters can be eliminated by applying a cut on the width of S2 pulses.
As for systematic fluctuations, it was empirically found that those with up to ~20% rms do not prevent correct reconstruction of the LRFs.
V-a Spatial resolution
Spatial resolution for both S1 and S2 was measured with Co calibration source. In the central part of the chamber, right below the source, the -rays cross the copper grid at normal incidence creating the sharpest contrast between open and shadow areas. In the reconstructed event distribution, this transition is smeared due to finite spatial resolution and, to some extent, by scattering in the 7-mm anode plate located below the grid. In other words, the sharpness of the edges of the projected image gives an upper limit for the spatial resolution of the detector for S2 signals. In Fig. 4, the distribution of the -positions of the reconstructed events is demonstrated for a narrow patch in the inner part of the detector (100 mm). The distribution is fitted with a convolution of a step function with the Gaussian giving resolution of 1.6 mm FWHM. The resolution worsens towards the edge of the fiducial volume due to combination of lower light collection and edge effects, becoming ~3 mm FWHM at =150 mm.
The spatial resolution for S1 can be estimated by comparing independently reconstructed coordinates for S1 and S2, Fig. 5(a). The difference between the two, shown in Fig. 5(b), is approximately normally distributed with FWHM of 15.0 mm for the whole fiducial volume and 13.0 mm for events with 100 mm. As the contribution of S2 resolution is obviously negligible, these values correspond to the spatial resolution for S1. Note that no energy selection was performed in these measurements so the ~10% admixture of 136 keV present in Co spectrum might marginally improve the results compared to what would be obtained using a pure 122 keV source.
V-B Energy resolution
As demonstrated in , there is strong anti-correlation between scintillation light and extracted charge for electron recoils in liquid xenon under an applied electric field. The reason for this is that part of the scintillation light comes from recombination. For the less dense electron tracks, the electron extraction efficiency is higher while recombination (and scintillation output) is lower. Thus, fluctuation of the electron track density from event to event leads to variations in light and charge outputs, which in a dual phase detector leads in turn to anti-correlated variations of S1 and S2 even for events of the same energy. Consequently, the best energy estimate for a dual phase detector is a linear combination of S1 and S2 light outputs. Fig. 6 shows the relationship between scaled light outputs for S1 and S2 for the events produced by -rays from the Co source. The scaling factors were chosen so that the mean of the distribution is at 125 units for both S1 and S2. One can see that there is indeed anti-correlation with S1 varying approximately by a factor of 3 more than S2.
A more detailed analysis of the plot on Fig. 6 yields the coefficients of the linear combination with the best energy resolution: . Using this formula, an energy resolution of 10.6% FWHM was obtained at 122 keV for the whole fiducial volume – see Fig. 7(a). For the central spot with 50 mm, where the effects from Compton scattering of incoming -rays in copper are minimal, the resolution is 8.1% FWHM and the two lines of the Co source are clearly resolved as shown in Fig. 7(b).
Position sensitivity is crucially important for a modern dark matter detector as it allows one to drastically reduce the background by considering only events inside an inner fiducial volume away from any detector surfaces. A position-sensitive detector also offers better energy resolution as it becomes possible to apply a position-dependent correction to the energy. In the case of a scintillation detector, the optimal performance of the position estimation algorithm depends on how well the set of the PMT LRFs describes the detector response to scintillation events.
In the present work, a novel method for iterative reconstruction of the light response functions from the calibration data acquired with uncollimated -ray source was developed and its suitability has been proven for the real detector. Using the reconstructed LRFs and applying the weighted least squares and maximum likelihood methods to position and energy reconstruction, the performance of the ZEPLIN-III detector was studied for 122 keV -rays. The measured performance for the inner part of the detector (100 mm) is as follows:
spatial resolution of 13 mm FWHM in the horizontal plane for scintillation signal (S1);
spatial resolution of 1.6 mm FWHM for electro-luminescence signal (S2);
energy resolution of 8.1% FWHM for the combined (S1 and S2) signal.
A more detailed description of the implementation of the position reconstruction algorithms and their impact on the WIMP search with ZEPLIN-III will be published as a separate paper. The developed method can also be applied in scintillation cameras for medical imaging for correction of non-uniformities and improving non-linearity associated with both the scintillation crystal and the PMT array as well as those due to the position reconstruction algorithm. The success of the new method in mitigating significant performance irregularities suggests that hardware components may be subject to less stringent requirements, thereby reducing the cost of scintillation cameras. The method can also be of advantage for regular quality control of gamma cameras.
The UK groups acknowledge the support of the Science & Technology Facilities Council (STFC) for the ZEPLIN–III project and for maintenance and operation of the underground Palmer laboratory which is hosted by Cleveland Potash Ltd (CPL) at Boulby Mine, near Whitby on the North-East coast of England. The project would not be possible without the co-operation of the management and staff of CPL. We also acknowledge support from a Joint International Project award, held at ITEP and Imperial College, from the Russian Foundation of Basic Research (08-02-91851 KO a) and the Royal Society. LIP–Coimbra acknowledges financial support from Fundação para a Ciência e Tecnologia (FCT) through the project-grants CERN/FP/109320/2009, CERN/FP/116374/2010 and PTDC/FIS/67002/2006, as well as the postdoctoral grants SFRH/BPD/27054/2006, SFRH/BPD/47320/2008 and SFRH/BPD/63096/2009. This work was supported in part by SC Rosatom, contract #H.4e.126.96.36.1999 from 10.03.2011.
-  E. Aprile and T. Doke, “Liquid xenon detectors for particle physics and astrophysics,” Reviews of Modern Physics, vol. 82, pp. 2053–2097, July 2010.
-  V. Chepel and H. Araújo, “Liquid noble gas detectors for low energy particle physics,” submitted to Journal of Instrumentation, 2012. http://arxiv.org/abs/1207.2292.
-  D. Y. Akimov et al., “The ZEPLIN-III dark matter detector: Instrument design, manufacture and commissioning,” Astroparticle Physics, vol. 27, no. 1, pp. 46–60, 2007.
-  J. Joung, R. Miyaoka, S. Kohlmyer, and T. Lewellen, “Implementation of ML based positioning algorithms for scintillation cameras,” IEEE Trans. on Nuc. Sci., vol. 47, no. 3, Part 3, pp. 1104–1111, 2000.
-  W. Cook, M. Finger, and T. Prince, “A thick anger camera for gamma-ray astronomy,” IEEE Trans. on Nuc. Sci., vol. 32, no. 1, pp. 129–133, 1985.
-  H. O. Anger, “Scintillation camera,” Review of Scientific Instruments, vol. 29, p. 27, 1958.
-  R. Gray and A. Macovski, “Maximum a posteriori estimation of position in scintillation cameras,” IEEE Trans. on Nuc. Sci., vol. 23, no. 1, pp. 849–852, 1976.
-  V. N. Lebedenko, H. M. Araujo, E. J. Barnes, A. Bewick, et al., “Results from the first science run of the ZEPLIN-III dark matter search experiment,” Phys. Rev. D, vol. 80, no. 5, 2009.
-  V. N. Lebedenko, H. M. Araujo, E. J. Barnes, A. Bewick, et al., “Limits on the Spin-Dependent WIMP-Nucleon Cross Sections from the First Science Run of the ZEPLIN-III Experiment,” Phys. Rev. Lett., vol. 103, no. 15, 2009.
-  D. Y. Akimov, H. M. Araujo, E. J. Barnes, V. A. Belov, et al., “Limits on inelastic dark matter from ZEPLIN-III,” Phys. Lett. B, vol. 692, no. 3, pp. 180–183, 2010.
-  D. Akimov, H. Araújo, E. Barnes, V. Belov, et al., “WIMP-nucleon cross-section results from the second science run of ZEPLIN-III,” Phys. Lett. B, vol. 709, pp. 14–20, Mar. 2012.
-  F. Neves et al., “ZE3RA: the ZEPLIN-III Reduction and Analysis package,” Journal of Instrumentation, vol. 6, no. 11, p. P11004, 2011.
-  M. Berger, J. Hubbell, S. Seltzer, J. Chang, et al., “NIST XCOM: Photon Cross Sections Database.” http://physics.nist.gov/xcom.
-  H. H. Barrett, W. C. J. Hunter, B. W. Miller, S. K. Moore, et al., “Maximum-Likelihood Methods for Processing Signals From Gamma-Ray Detectors,” IEEE Trans. on Nuc. Sci., vol. 56, no. 3, Part 2, pp. 725–735, 2009.
-  J. Prescott, “A statistical model for photomultiplier single-electron statistics,” Nuclear Instr. and Meth., vol. 39, no. 1, pp. 173–179, 1966.
-  H. Gale and J. Gibson, “Methods of calculating pulse height distribution at output of a scintillation counter,” J. of Sci. Instr., vol. 43, no. 4, pp. 224–228, 1966.
-  N. Clinthorne, W. Rogers, L. Shao, and K. Koral, “A hybrid maximum-likelihood position computer for scintillation cameras,” IEEE Trans. on Nuc. Sci., vol. 34, no. 1, pp. 97–101, 1987.
-  T. Ling, T. H. Burnett, T. K. Lewellen, and R. S. Miyaoka, “Parametric positioning of a continuous crystal PET detector with depth of interaction decoding,” Physics in Medicine and Biology, vol. 53, no. 7, pp. 1843–1863, 2008.
-  G. Cowan, Statistical Data Analysis. Oxford University Press, USA, June 1998.
-  T. D. Milster, L. A. Selberg, H. H. Barrett, A. L. Landesman, and R. H. Seacat, “Digital position estimation for the modular scintillation camera,” IEEE Trans. on Nuc. Sci., vol. 32, pp. 748 –752, Feb. 1985.
-  F. Neves, V. Solovov, V. Chepel, M. Lopes, et al., “Position reconstruction in a liquid xenon scintillation chamber for low-energy nuclear recoils and -rays,” Nuclear Instr. and Meth. A, vol. 573, pp. 48–52, Apr. 2007.
-  A. Lindote, H. Araujo, J. Pinto da Cunha, D. Akimov, et al., “Preliminary results on position reconstruction for ZEPLIN III,” Nuclear Instr. and Meth. A, vol. 573, pp. 200–203, Apr. 2007.
-  J. Joung, R. S. Miyaoka, and T. K. Lewellen, “cMiCE: a high resolution animal PET using continuous LSO with a statistics based positioning scheme,” Nuclear Instr. and Meth. A, vol. 489, pp. 584–598, Aug. 2002.
-  E. Conti, R. DeVoe, G. Gratta, T. Koffas, et al., “Correlated fluctuations between luminescence and ionization in liquid xenon,” Phys. Rev. B, vol. 68, p. 054201, 2003.