|Home | About | Journals | Submit | Contact Us | Français|
We have developed a Monte-Carlo photon-tracking and readout simulator called SCOUT to study the stochastic behavior of signals output from a simplified rectangular scintillation-camera design. SCOUT models the salient processes affecting signal generation, transport, and readout of a scintillation camera. Presently, we compare output signal statistics from SCOUT to experimental results for both a discrete and a monolithic camera. We also benchmark the speed of this simulation tool and compare it to existing simulation tools. We find this modeling tool to be relatively fast and predictive of experimental results. Depending on the modeled camera geometry, we found SCOUT to be 4 to 140 times faster than other modeling tools.
SCOUT is a Monte-Carlo simulation tool to model Scintillation Camera OUTput in response to incident high-energy photons. SCOUT is not meant to be a general-purpose, high-fidelity model, such as Geant4 [1, 2], nor does it model physics of the source, object, or collimator. Instead, it models only salient processes affecting signal generation, transport, and readout of rectangular scintillation-camera geometries (figure 1). Narrowing the scope of application in this way enables us to simplify the photon tracking routine and helps to speed up simulation of scintillation camera response.
SCOUT originated at the University of Arizona, where it was used to study how basic camera design affects the joint distribution of multiple output signals  and to examine methods of calibrating detector response as a function of depth . This tool is now being supported at the University of Washington, where we are using it for surface-treatment optimization of the depth-of-interaction detector design based on light sharing of paired crystals (dMiCE) [5–7]. We are further interested in SCOUT as a complimentary utility to the Simulation System for Emission Tomography (SimSET) , a simulation tool also maintained by the University of Washington for rapid photon tracking in tomographic systems.
In support of the dMiCE project, we have expanded the variety of modeled optical-surface treatments that SCOUT supports and we have added options for modeling arrayed scintillation crystals, dual-ended optical readout, and noise processes for Geiger-Müller avalanche photodiodes. Furthermore, we have performed extensive code verification and model validation of SCOUT.
Presently, we report several validation and benchmark tests, comparing with experimental results [5–7, 9, 10] and other simulation tools (Geant4 [1, 2] and Detect2000 [11, 12]). We find this modeling tool to be relatively fast and predictive of experimental results for simple detector and readout geometries.
SCOUT is written in ANSI-C and adopts a procedural programming paradigm. For given camera and radiation parameters, a simulation of the scintillation camera output entails tracing the generation, transport and readout for an ensemble of gamma rays with fixed input position and/or direction. Successive states of the signal path (gamma-ray interaction, scintillation, optical transport, photodetection, and acquisition) and sub-states of the light paths (each volume and boundary of the detector) are treated in separate subroutines; thus, in addition to being recorded for output, the state of an event at each stage can be traced at runtime. Development and testing of SCOUT was aided using process-flow diagrams (e.g., figures 2 and and6),6), compile-time verbosity settings, and optional photon-track history logging. Version control of SCOUT is managed with SVN and code changes are logged therein.
SCOUT uses a 64-bit Mersenne-Twister algorithm [13,14] as coded by  for uniform pseudorandom number generation. For generating a non-uniform distribution, we either use the Direct method (inverse of the cumulative distribution function, when available) or Composition-Rejection method [16–18]. For example, the composition-rejection algorithm used for the Klein-Nishina distribution of Compton scatter was obtained from the EGS code system report . For several standard distributions, such as exponential and binomial distributions, we used source code derived from the RANDLIB.C random number library  updated with the Mersenne-Twister algorithm for uniform random samples.
SCOUT presently assumes the modeled camera is a rectangular array (one or more) of scintillation crystals with optional photodetectors (PD) and optical spacer (OS) on either or both the front and back surfaces (figure 1). Multiple photodetectors on front or back surfaces are arranged as a regular rectangular array. The active area of each PD is assumed to be uniform within a support region centered at its PD-grid location. This support region can be circular, square, or an intersection of a concentric circle and square.
SCOUT generates signals as a series of random processes, including gamma-ray transport, scintillation yield, optical transport, photodetection, and electronic acquisition (figure 2). Photon energies up to 1.022MeV (pair production is not modeled) are input at a specified position/direction and subsequently interact and/or scatter in a fixed or stochastic manner within the detector volume. The number of scintillation photons is randomly generated from an energy-dependent normal distribution at each interaction site. Optical photons are tracked until they are absorbed or detected. Inputs of the photodetectors are conditioned (accounting for random amplification, saturation, crosstalk, and after-pulsing) and electronically read out (accounting for conversion gain and additive Gaussian noise).
The response to an ensemble of gamma rays of specified input direction and energy between 0. 1keV and 1.022 MeV are simulated by each execution of SCOUT. These photons can be propagated to a random depth of interaction (DOI) within the detector starting from its entrance face, or they can be forced to interact at a specified location in the detector. Using beam-width parameters, the position of each photon can also be randomly offset laterally about the mean input location to simulate a uniform-square or uniform-circular beam profile (figure 3). Beam divergence with depth is not modeled.
The pathlength through the detector of a high-energy photon (>0. 1keV) is sampled from a truncated exponential distribution (figure 4b) described by user-specified energy-dependent attenuation data (figure 4a). High-energy photons that exit the detector volume are assumed to escape detection (i.e. no back scatter is simulated). Photon interactions can be forced to be photoelectric or they may be randomly chosen according to the energy-dependent branching ratios for the Photoelectric Effect, incoherent (Compton) scatter, or coherent (Rayleigh) scatter. Attenuation data for this purpose is available through NIST . For simplicity, we presently assume that photoelectric secondaries (i.e. characteristic X-rays) are reabsorbed local to the photoelectron. We sample the energy of a Compton photon (figure 5) as described by Nelson and Nelson  using the Butcher and Messel form of the Klein-Nishina differential cross-section . The deflection angle of the Compton photon is then constrained by the Compton Equation and the azimuth scatter angle is uniformly sampled over the interval [−π, π]. Scattered photons with energy falling below 1 keV are assumed to fully deposit their energy without further transport. Rayleigh scatter is modeled in a similar fashion as Compton scatter, using the Composition-Rejection method as described by Nelson, Hirayama, and Rogers .
A normally distributed random number of optical photons are generated at each simulated interaction site. The mean and variance of the sampled distribution are a function of the deposited energy and interaction type. These scintillation statistics are linearly interpolated from user-supplied data; such data can account for non-proportional effects [31–32]. We assume the ensemble radiance of the scintillation light from an interaction site to be isotropic. Detection times of these optical photons are computed as the scintillation delay time plus travel time in each material. The travel time in material i is Li/nic, where Li is travel distance, ni is refractive index, and c is the speed of light in vacuum. Creation time, t, is randomly sampled from a double-exponential pulse shape (τr + τf) τf−2 e−t/τf (1 − e−t/τr), where the configurable parameters, τr and τf, are the rise and fall time constants, respectively.
Scintillation light is tracked through the crystal and window layers of the detector and can be redirected in the bulk or at the surfaces of these volumes; this process continues for each generated optical photon until it is absorbed by the bulk, photodetector, or other surface (figure 6). Scatter and attenuation within the bulk of the scintillator or windows are modeled by user-specified linear attenuation coefficients for each material. At the surfaces of these volumes, light can be absorbed with specified probability or be redirected according to a given reflector or interface model. A unique optical-surface model can be specified for each of the scintillator and window surfaces, including those between crystal-array elements.
Current models for reflectors and interfaces supported by SCOUT are illustrated in figure 7 and a measure of their scattered light distributions are shown in figure 8. For each optical model in Figure 8, an ensemble of 107 optical photons were simulated to isotropically radiate from a point 10 mm from this surface. These photons were then grouped by angle-of-incidence (1-degree bins) and a 2D polar histogram of their outgoing direction was generated. The specular reflector model simply reverses the component of the photon-direction vector normal to a surface. The specular transmission model uses Snell's Law. The Lambertian reflector and interface models yield a radiant intensity that obeys Lambert's Cosine Law (see Fig. 7). The multi-faceted diffuse (MFD) model consists of locally specular surfaces whose slopes are randomly sampled from a zero-mean normal distribution of configurable variance; the MFD model assumes facet sizes larger than is typical for coherent effects, and multiple facet interactions are accounted for only when a photon is incoherently redirected back into the surface. The Diffuse-Cone model yields a uniform radiant intensity within an opening angle about the specular directions (both reflected and refracted); this opening angle is a configurable linear function of the incident angle . The dual-specular and dual-MFD models are similar to their standard counterparts except that an optical coupling material (of negligible thickness, but distinct refractive index) is included at the interface. The Ideal Retroreflector simply reverses the incident direction. However, the Cube-corner Retroreflector is an accurate geometric model of a rectangular cube-corner array; cube size and orientation are configuration parameters. The isotropic reflector and interface models isotropically randomize the photon direction in the appropriate hemisphere. The bidirectional lookup is for use with experimentally calibrated scatter distribution data (e.g. [23–25]).
The probability of reflection at interfacing volumes is computed using the Fresnel equations; we assume unpolarized light for this purpose and average reflectivity over both polarizations. Except for the MFD and 3-layer-MFD, we measure the angle of incidence with respect to the average (bulk) surface normal. For the MFD and 3-layer-MFD models, the angle of incidence is measured from the local surface normal. For the bidirectional lookup model, reflectivity is built into the lookup table.
For more flexibility in optical surface modeling, a mixture of any two models described above can be used for each surface. For side-facing surfaces of the scintillator (perpendicular to the gamma-ray entrance surface as shown in figure 1), the mixture fraction of surface treatments can be configured as a polynomial function of depth. Furthermore, the scintillator sides can be subdivided into two regions that are symmetric about a vertical midline. A different depth-dependent mixture model can be specified for each region of each side of each crystal element. Additionally, for each surface we model a probability for absorption and a probability of an “air gap”. Randomly air gaps are a real issue that can affect the refractive properties at interfaces and bounding surface treatments of a detector.
Rectangular arrays of photodetectors (PD) can be specified on the front and/or back surfaces of the detector; these arrays are configured separately. For each PD array, the user specifies its width and the number of identical pixels it is divided into.
To be detected, a photon must refract into the PD material, fall within an active support region, and produce an electric signal. The refractive index of the PD material is assumed to be spatially uniform and spectrally averaged. The active area of each PD is assumed to be uniform within a support region centered at its PD-grid location. This support region can be circular, square, or an intersection of a concentric circle and square. Spectral-averaged quantum efficiency is used to randomly determine if a photon entering this support region will produce an electrical output signal. For single-sided readout configurations and when no optical-photon-track output is required, SCOUT is able to apply this quantum efficiency loss before optical photon tracking in order to reduce simulation time.
If the user flags the use Geiger-Müller APDs (GMAPD), then additional configuration parameters are used to describe GMAPD-specific processes; these include the number of microcells, microcell fill factor, pixel crosstalk probability, after-pulse probability, after-pulse characteristic delay time, bias recovery time, and a choice of GMAPD models. The choices of GMAPD models are either stochastic  or event-driven , which offer a tradeoff in simulation speed and accuracy.
Output signals are an amplified integral of signals resulting from photoelectrons, thermal electrons, crosstalk and after-pulsing. The user-specified integration period begins at the first photon arrival time. Each integrated input charge undergoes random amplification characterized by user-supplied mean gain and gain variance. The gain and variance of after-pulse events are scaled by (1 − e−t/τap), where τap is the after-pulse characteristic delay time (assumed to be exponentially distributed). An additive Gaussian variance term can also be specified to model subsequent independent electronic-noise processes.
Optional outputs of the SCOUT tool include: gamma-ray track history, optical photon track histories, track statistics, list-mode data of photodetector-input and amplifier-output signals, and signal statistics. Due to the large number of optical photons per gamma ray, the option to output optical track information is mutually exclusive with output of gamma-ray track or signal data. However, gamma-ray track and signal data can be output together.
We examined the fidelity of SCOUT simulated signals as compared with measurements made for both discrete and monolithic detector geometries.
We consider a dMiCE crystal pair [5–7], which is optically isolated from other crystal pairs by 3 layers of Teflon tape. The material and optical properties of the interface for this crystal pair is something that our lab hopes to optimize using a tool such as SCOUT. However, here we validate the SCOUT model by comparing its simulated detector response for a configuration that we have also experimentally measured. The experimental interface for this purpose (figure 9) consists of a triangle-shaped mirror film (Enhanced Specular Reflector, VM2002 by 3M, St. Paul, MN), which changes the amount of light shared by the two crystals versus interaction depth (zDOI).At the base of this triangular reflective interface, is a shorter rectangular area with a double layer of the VM2002 material. The entire interface (reflector and open area) is optically coupled using a high-index mounting media (Melt-mount 1.704 ® by Cargille Labs, Cedar Grove, NJ) .
In figure 9, we illustrate the experimental dMiCE detector setup used to calibrate the mean detector response (MDR) of photopeak events as a function of zDOI We also indicate the model parameters used to simulate a comparable response from SCOUT. The standard deviation of electronic noise was measured to be 0.4 Analog-to-Digital Units (ADU) from the root-mean-square deviation about the pedestal at zero bias. By varying power from a fast (50ns) low light-emitting diode (LED), gain is measured from the slope of variance over mean response at low light conditions (when less than 5% of microcells fire). An average conversion gain of 3.9 ADU/photoelectron was used for each of the channels in the model. We used a 200-nsec integration period and manufacturer specification (Zecotek, MAPD-3N) for dark count rate is 8 MHz per pixel. We use a photodetection efficiency of 26% (accounting for fill-factor, spectrally-averaged quantum efficiency for LFS, and avalanche probability). The GM-APD pixels were oversized for the crystals used (3-mm vs. 2-mm width, respectively), but the crystal boundary was positioned at the boundary of two photosensors. Specifications for the PDE, afterpulsing, crosstalk, and recovery time that were used in simulation are also shown in figure 9.
Optical-model parameters that we did not have information for were selected to minimize the square error of the simulated versus experimental MDR. These parameters included the amount of air gaps in the interface; the probabilities of absorption and transmittance of the mirror film; and the mixture of specular and Lambertian models of the side reflectors. The resulting values (also shown in figure 9) were quite plausible given other anecdotal observations we made concerning these parameters. For instance, the prediction of air gaps in the interface was initially not expected, but subsequently visually verified with the aid of a confocal microscope.
In figure 10 we show the experimental procedure for identifying the crystal that was hit (gamma-ray interacted in) and the abutted crystal. Charge-integrated signals (QDC) of both photodetectors (PD1 and PD2) are recorded for 25,000 events at several zDOI Photopeak events within each crystal are selected using a threshold (indicated on the scatter plot of figure 10 by dash lines). For each hit-crystal spectrum, we selected a threshold at the mid-point between the photopeak and the Compton backscatter-peak, which are well separated. Using the data in just the upper left and lower right quadrants formed by these thresholds, we are able to clearly identify the hit and abutted crystal for photopeak events. In this way, we are also able to clearly distinguish photopeak events in the abutted crystal from non-photopeak interactions in that same crystal. The mean responses for the hit and abutted channels at each depth are then computed from this filtered data. We found the computed photopeak means to be relatively insensitive to a few percent variation of this threshold. For the interface treatment used, we find this hit-channel identification method to be unambiguous even for zDOI near the optically coupled region of the interface. If a truly transparent region of the interface existed, then the crystals might have to be separately illuminated to independently calibrate their responses. The photodetector signals that are filtered in this way are then used to calibrate hit and abutted detector response to 511-keV photoelectric gamma-ray interactions as a function zDOI.
Repeating this acquisition process for simulated data, we then compare simulated to experimental signal spectra in figure 11 and depth dependence of the MDR for photopeak events in figure 12. Simulated results are scaled by a measured conversion gain of 0.48 ADU/photoelectron for this comparison. The signal-to-noise ratio computed by repeated trials was better than 140 for all measured values and better than 430 for all simulated values. Therefore, error bars were not included in figure 12. The residual difference in simulated and measured depth dependence seen in figure 12 is thought to be due to minor discrepancies in both the mirror film shape and depth alignment of the coincidence beam.
Sensitivity plots of the objective function (total square error of simulated versus experimental MDR summed over position and photodetectors) on the optimized model parameters are shown in figures 13(a–b) and 13(c–d). To fit these parameters, we first localized the least-square-error parameters using a Nelder-Mead search algorithm with 1000 simulated gamma rays per depth. The best-fit parameter values were then refined using 10,000 gamma rays per depth and computing the parabolic minimum of the square error along each dimension. The global minimum least sum-square error for these parameters is found at 6.9% mirror-film absorption, 98.7% mirror-film reflectivity, mixture of 69.0% specular and 31.0% Lambertian for side reflectors, and 54% probability of an air gap between crystals. The high air-gap probability (visually confirmed afterwards by confocal microscope) is due to poor adherence of the Melt-Mount media (made difficult by multiple film layers) and not degassing the Melt-Mount.
We also compared SCOUT simulated output with experiment for a continuous crystal detector called cMiCE [9–10]. All simulation parameters in this case were selected based on prior knowledge or empirical observation. A diagram of the detector and its experimental calibration procedure are shown in figure 14. A 50-mm-square, 8-mm-thick Lu1.2 Y0.8 SiO5:Ce scintillator (LYSO) from Saint Gobain, Nemours FR, was read out using an 8-by-8 multi-anode photomultiplier tube (H8500 from Hamamatsu). We used spectral-average quantum efficiency for LYSO of 22% and a 91% fill-factor for the active area. All surfaces, except the readout side, were roughened (300 grit), the sides were blackened, and the top coated white . We simulated the sides to be 95% absorptive with 80% air gap and the top to be a MF diffuse reflector. For comparing to the experimental spectral response for an ensemble of collimated 511-keV gamma rays, we used NIST attenuation data  for LYSO with SCOUT to randomly select interaction type (Compton or photoelectric) and interaction depth. The experimental mean detector response (MDR) for photopeak events was measured with a physically collimated and coincidence-collimated22 Na point source. For a simulated comparison of the MDR for photopeak events, an ensemble of forced-photoelectric interactions at fixed depths was used.
We represent this calibration process in simulation with a 0.6-mm-diameter dual-energy (511-keV and 1274-keV) gamma-ray beam that is normally incident as shown in figure 14. Twenty thousand events were simulated at each beam position. The ratio of 511-keV to 1274-keV gamma-ray interactions that we simulated was 18:1, which accounts for the branching ratio, the solid angle subtended by the collimating aperture, and the probability ratio for attenuation of these energies. We do not account for collimator penetration or scatter. We compare the resulting signal spectra in figure 15 and the MDR of photopeak events versus 3D position in figures 16 and and17.17. To better capture the relative spatial dependence of the photodetector input signals, the MDR for each channel was divided by its gain. Reported values are in units of photoelectrons. Determination of the MDR as a function 3D-position from 2D-calibration data is accomplished by the maximum-likelihood clustering method .
We compared the CPU time per interacting gamma ray of SCOUT with that of Geant4 [1, 2] and Detect2000 [11, 12]. Since the optical track history depends on detector geometry, we examine two scintillator sizes. Also, since Detect2000 does not track gamma rays, we used a fixed gamma-interaction position (crystal center) and forced photoelectric interaction for comparison with Detect2000. However, for comparison with Geant4, we used a random interaction depth and allowed scatter. Results are presented in Table 1. The majority of CPU time spent in these simulations is due to optical photon tracking (simply because there are so many optical photons generated per gamma ray). Gamma-ray tracking and readout-noise simulation take relatively little time.
The speed improvement in SCOUT is attributed to four reasons. First, since SCOUT assumes a fixed rectangular geometry, intersections of photons with surfaces can be determined by a few trigonometric tests as opposed to iterating over all bounding surfaces. Second, there is no time stepping in SCOUT; bulk attenuation and scatter are randomly sampled for each track segment from a cumulative attenuation distribution as a function of pathlength. Third, when optical photon track data is not requested for output, the number of photons tracked is reduced by a binomial selection of the maximum photodetection efficiency (PDE) of any one photodetector in the camera (the difference in PDE for lower-PDE photodetectors is tested later). Finally, SCOUT source code has been heavily optimized to minimize redundant calculations and it diligently uses lookup tables where possible.
The SCOUT modeling tool has proven to be easily adaptable for the arrayed-scintillator studies we are conducting. Several options and combinations of surface treatments permit significant flexibility and realism in modeling optical transport in a scintillation camera. Furthermore, a comparison with experimental data has validated SCOUT accuracy in modeling the underlying physical processes in scintillation camera signal generation. Simulated results from this model have reproduced measured signal statistics with fair accuracy. We have conducted sufficient verification of this modeling tool to proceed with our surface-treatment optimization study.
The present simplicity of SCOUT model permits it to run relatively fast. In one second of a single 2.26 GHz Intel Xeon processor, SCOUT tracked 4 gamma rays and ~60,000 optical photons with about 2 million diffuse-surface interactions. This speed will enable us to conduct a wide range of parametric studies on detector design. Current limitations of SCOUT that we aim to address in a computationally efficient manner for future releases include: non-rectangular detector geometries, tracking of fluorescent X-rays, spatial and ensemble variations of photodetector properties, spectral and polarity effects on light production and transport, and incident-beam divergence.
†This work was supported in part by NIH/NIBIB grants EB002117 and EB002035, and by DOE grant DE-FG02-08ER64676.