ScintSim1: A new Monte Carlo simulation code for transport of optical photons in 2D arrays of scintillation detectors.

Two-dimensional (2D) arrays of thick segmented scintillators are of interest as X-ray detectors for both 2D and 3D image-guided radiotherapy (IGRT). Their detection process involves ionizing radiation energy deposition followed by production and transport of optical photons. Only a very limited number of optical Monte Carlo simulation models exist, which has limited the number of modeling studies that have considered both stages of the detection process. We present ScintSim1, an in-house optical Monte Carlo simulation code for 2D arrays of scintillation crystals, developed in the MATLAB programming environment. The code was rewritten and revised based on an existing program for single-element detectors, with the additional capability to model 2D arrays of elements with configurable dimensions, material, etc., The code generates and follows each optical photon history through the detector element (and, in case of cross-talk, the surrounding ones) until it reaches a configurable receptor, or is attenuated. The new model was verified by testing against relevant theoretically known behaviors or quantities and the results of a validated single-element model. For both sets of comparisons, the discrepancies in the calculated quantities were all <1%. The results validate the accuracy of the new code, which is a useful tool in scintillation detector optimization.

Introduction

Treatment verification is an important part of modern radiotherapy (RT) in which new technologies allows us to have tight target margins, and hence often escalate the prescribed dose to the target. Because of different sources of uncertainties,[1] anatomical information (just before or during the treatment delivery) is usually needed for these technologies to be truly beneficial.

Megavoltage cone-beam computed tomography (MVCBCT), as an imaging modality that provides three-dimensional (3D) images using the linear accelerator's MV photon beams, has been considered advantageous in producing treatment-time 2D and 3D images; 2D images for comparison with 2D reference images and 3D for visualizing 3D anatomy with soft-tissue contrast at the time of treatment, as well as straight dose calculation in modern image-guided RT.[234]

However, the very low quantum efficiency (QE) at MV energies imposes significant restrictions on the detectors used for imaging. These detectors, referred to as electronic portal imaging devices (EPIDs), have been developed based on liquid-filled matrix ionization chambers, phosphor screen with a lens-camera combination, and currently mostly amorphous silicon flat panel imagers.[5] With the exception of the matrix ionization chambers EPID, these imagers use a convertor layer to convert the radiation to light which is then detected by the underlying photo detector.

Using thick, segmented scintillation crystals with optically-isolated detector elements, the QE of these kind of detectors for 6 MV X-rays was demonstrated to increase from 2 to 3% in the commercial devices to about 18% by using 1 cm thick CsI (Tl) scintillation crystal arrays.[6] In addition, relatively low-dose MVCT projections could be taken using this device at frequent intervals and with a high speed.[3] Since then, empirical studies as well as Monte Carlo simulations have further shown that CsI (Tl) crystals up to 40 mm thick can improve detective quantum efficiency (DQE) values at very low dose (˜20% at zero spatial frequency).[78]

The detection process in scintillation materials consists of two stages: i) Energy deposition following interactions of high energy photons with the scintillation material and ii) production of optical photons and their transport inside the scintillator. In most publications regarding analytical or simulation-based scintillation detector studies, optical photon transport has not been included, for example in.[91011] Inclusion of optical transport will surely affect some performance parameters of these detector models. For example, reflectivity of the walls, absorption and scattering in the crystal, and other optical effects have been shown to affect some aspects of the detector performance.[121314]

One of the main reasons for ignoring optical photon transport in simulations is that only a very limited number of Monte Carlo simulation packages exist for modelling light transport inside scintillation crystals. Some of these models are described in.[15161718] Also, in some models, detailed information (for example, the various quantities of interest describing the history of the scintillation light rays) is not always accessible or configurable.

Here, we present the development and validation work for ScintSim1, an in-house code for simulation of light transport inside multi-element scintillation crystals. The basis for the new code was a single-element model which has been previously reported.[12] The new model will help in optimization of the structure of this kind of detector.

Materials and Methods

Description of the model

An optical Monte Carlo simulation model for one single element of a scintillation detector, written in XL FORTRAN programming language, was described in reference.[12] We revised various aspects of the model and developed it in the MATLAB programming environment to model 2D arrays of scintillation crystals. The code was rewritten and revised based on the existing code for single-element detectors, with the additional capability to model 2D arrays of detector elements.

The MCNP4c Monte Carlo code[19] was first used to simulate X-ray and electron transport in the individual elements of the detector. The depth-dose distribution of the energy deposited in each element in 1 mm layers was then given to the model as an input file.

Each detector element is in the shape of a rectangular parallelepiped with configurable dimensions, material, etc., [Figure 1]. The optical model accepts the deposited energy distribution within the array from existing ionizing radiation transport packages, generates the corresponding number of optical photons, and follows each optical photon history through the detector element (and, in case of cross-talk, the surrounding ones) until it reaches a configurable receptor and is scored (according to its position of incidence) or is attenuated. Beside different parameters that the program keeps recording during a simulation, some of which will be described in the Results section, the code computes the distribution of the incident photons on a 2D pixelated screen with configurable pixel size.

Figure 1

Central element of the array and the simulation coordinate system. The source is located above this detector element and the transmitted optical photons exit from the bottom face

All the detector elements in ScinSim1 were modeled to have a Lambertian-type coating (a class of diffuse surface) on five faces with the exit face modeled as polished. Other parameters such as element dimensions, crystal material, coating reflectivity, etc., are input parameters to the program.

Septal walls between detector elements were modeled to be extremely thin (as is normally the case in practice) and their transparency for the non-reflected rays was modeled to be either 0 or 100%. Thus, depending on these values, an optical photon that does not reflect from a wall is absorbed in it (and ignored) or enters the adjacent element (and the simulation continues in that element, with a finite probability of reaching further elements too).

Testing and validation of the model

The correctness of the individual components of the model was tested by comparing the results with expected values. Some of these tests are listed below. For those tests that are performed during a full simulation, a 3 × 3 × 10 mm3 detector element made of CsI (Tl) crystal was modeled, 10 mm being the thickness dimension.

Randomness of the random number generator (RNG)

The distribution of the generated pseudorandom numbers was tested. Here, we tested and used the ‘rand’ command in the MATLAB R2010a software with its default settings. This command returns pseudorandom numbers with uniform distribution within the interval (0,1). We performed the test by calling 107 random numbers and binning them into 8000 intervals of width 1.25 × 10-4. Then, we compared the mean and standard deviation of the distribution with the expected theoretical values.

Fresnel coefficients of reflection as a function of angle of incidence

These coefficients are used when the optical photons are incident on the polished exit face (scintillator/air interface). For those optical photons that do not undergo internal reflection, the probability of reflection (Pref) is given by Fresnel equations. More details are available in reference.[12]

For an interface that separates two materials with indices of refraction n1 = 2.2 and n2 = 1, we obtained Pref and Fresnel reflection coefficients as a function of angle of incidence and compared them with expected theoretical values.

The spatial distribution of the points of incidence (on the exit face) for the transmitted rays

This was tested by binning the distance of incidence points from each axis into 15 intervals of 0.1 mm width and obtaining the number of transmitted rays on each strip to find the distribution.

Presence of any unexpected bias in the direction cosines (DCs) of the transmitted rays

The x and y DCs are expected to have a mean value of zero [Figure 1] because no preferable direction exists. The z DC is expected to be negative (for a transmitted ray the angle is always larger than 90°). For photons whose distance of incidence points were binned in the previous test, the averaged DCs in each bin was plotted.

Polar angle probability density of rays reflected from a Lambertian surface

Reflections from the five faces of each element which were modeled as Lambertian reflectors were dealt with as follows. The polar angle probability density of the reflected ray for such a reflector can be shown to be pθ(θ) = sin2θ,[12] where θ is the polar angle of the reflected ray relative to the surface. For calculating such a probability density from the uniformly distributed random numbers, we used θ=1/2cos(1−2RN) to compute the required angular distribution. RN represents a random number in this equation. Then, we compared the obtained probability density during a simulation with what was expected (pθ(θ) = sin2θ).

The results of this model were also compared with those of the experimentally-validated single-element code. At least 7.3 million optical photon histories were followed in each case.

Results

We theoretically tested various modules of the new code for an array of parallel crystals. Here are the results of the tests we discussed in the previous section with their corresponding heading numbers:

The results of the test are displayed in Figure 2. The observed mean counts per interval was 1,250.000 with a standard deviation value of 35.353, which are very close to the values expected from a uniform distribution of truly random numbers (107/8,000 = 1250; = 35.355. Also, it can be seen from Figure 2 that the distribution is symmetric around the mean.

Figure 2

Distribution of the random numbers produced by the ‘rand’ command in MATLAB R2010a (107 calls binned into 8,000 intervals)

Figure 3 shows Pref and Fresnel reflection coefficients (rperp and rpara) as functions of the angle of incidence. The subscripts denote rays with E-field perpendicular and parallel to the plane of incidence, respectively. The curves are of the expected shape and the critical angle (θc = 27.0°) and Brewster's angle (θp = 24.5°) both have the expected values. Note that Pref = 1 for θi > θc. For normal incidence, rpara, rperp, and Pref have the expected values of −0.37, 0.37, and 0.14, respectively, which are also seen in Figure 3.

Figure 3

Fresnel coefficients of reflection and the probability of reflection as functions of the angle of incidence on the exit face of a scintillator with n1 = 2.2. The values of the theoretically expected θc and θp are also shown

For a 10 mm thick crystal, the points of incidence were seen to be fairly uniformly distributed on the surface with a small bias towards the center, as shown in Figure 4. This behavior is expected because photons are more probable to fall near the center of the surface rather than on the peripheral regions.

Figure 4

Spatial distribution of the points of incidence (on the exit face) for the transmitted rays

The averaged DCs in each bin are displayed in Figure 5. The mean values of the x and y DCs (cosα and cosβ in Figure 5) ranged from −0.0100 to 0.0067. For z DC, the corresponding value was −0.698 ± 0.002, giving an average transmission angle of 45.71° relative to the normal of the exit face, which confirms no unexpected bias for a 10 mm thick crystal.

Figure 5

Distribution of the mean direction cosine of the transmitted rays against the position of the point of incidence on the exit face

Figure 6 shows the value of pθ(θ) obtained during a simulation, plotted against θ. Also shown are some expected values calculated from pθ(θ) = sin2θ. As can be seen, the agreement with the theoretical expected trend is good, where the equation of the best fit curve to the data is pθ(θ) =1.004sin (1.9996θ + 0.001).

Figure 6

Comparison of the simulated and expected polar angle probability densities of a ray reflected from a Lambertian surface

Some of the scored quantities that were obtained during a simulation are compared with the single-element code in Table 1. As can be seen, there is close agreement between the two sets of data.

Table 1

Comparison of the different parameters calculated during a simulation for the original model (experimentally-validated single-element code) and the new model (an array of the same elements)

Discussion and Conclusion

A number of groups have studied the use of thick-segmented crystalline scintillators as the convertor in MV imaging devices.[789121320] The detector elements in these 2D arrays of scintillation crystals are separated by optically opaque/reflective septal walls. Such arrays could be built thick enough for the detection of high energy photons while maintaining the spatial resolution within an acceptable level. With the exception of some theoretically-demonstrated anomalous behavior, frame averaging can also be utilized to improve signal-to-noise ratio.[21]

A new Monte Carlo simulation program to perform optical photon transport in an array of parallel scintillation crystal detector elements was described. Such a Monte Carlo simulation tool is very useful in optimization of the design of segmented crystalline scintillators. One of the advantages of such an in-house model is the possibility of adding and scoring different quantities according to the special application to achieve insight into the involved processes and their relative importance. Some of these quantities are shown in Table 1.

The distribution of the generated pseudorandom numbers is an important aspect of every Monte Carlo simulation code. Pseudorandom numbers are obtained either by calls to RNGs provided by the computer operating system or computing language, or alternatively using inline coded RNG algorithms. Here, we tested the ‘rand’ command available in the MATLAB R2010a software and found it to have a sufficiently long period and produce acceptable uniformity (lack of bias) in the generated number. This approach is different from that used in the original single-element code where the SLAC-RAN6 RNG was coded inline.[12]

To keep the model simple, general, and independent of the type of the optical detection system; the screen (optical receiver) was modeled to be 100% efficient in detecting optical photons and it could also be placed at a configurable distance from the exit face of the scintillator.

The results of the new model show good agreement with two separate sets of reference data: (i) Various individual modules of the code were tested against relevant theoretically known behaviors or quantities, and (ii) the results of the new code were compared with the previously reported single-element model when identical parameters were set for both models. Some of these results were shown in Table 1 and Figure 5. As the single-element code had been experimentally-validated,[12] this constitutes an indirect method of validating the new model. Direct experimental validation of the new code was not practicable due to the lack of access to the required materials and equipment.

All the data presented in Table 1 depend only on events inside the elements. Otherwise, we could not compare the results of a single element with those of an array. It should be emphasized that the results depend on the detector configuration (e.g., detector element wall coating type and reflectivity or crystal thickness) and the tests were carried out for identical configurations; for other configurations, we would have obtained different values.

The model described here is a special-purpose one for scintillation detectors. To reduce simulation run time, only the main effects were modeled. Moreover, for the same reason and to reduce complexity, unpolarized light was assumed by allocating equal probability to Fresnel intensity reflection coefficients for polarized light perpendicular and parallel to the plane of incidence. This is a reasonable assumption given the very large number of reflections involved. It offers simplification by obviating the need to assign a polarization state to each ray.

Some other optical simulation models such as DETECT2000[15] or PHOTON[16] allow simulation of optical properties as a function of light wavelength. The wavelength is randomly selected from a user-specified spectrum after each interaction. In ScintSim1, the wavelength of the scintillation light is accounted for by using single effective values when the user enters the magnitudes of quantities such as refractive index, surface reflectivity, and optical attenuation length in an input file.

The ability to model surfaces with different levels of roughness would be useful too. Introducing a new parameter to the program describing roughness of the surfaces, as in Geant4 Monte Carlo code,[17] would provide the possibility of studying the effects of surface roughness. In ScintSim1, a diffusely reflecting (Lambertian) surface is modeled but not a rough (non-smooth) one.

Although adding the above-mentioned effects would enhance the flexibility of the model, in any case, the relevance and importance of each additionally modeled effect should be weighed against the extra burden of simulation run time and model complexity.

The code described here adds more capabilities compared to the single-element model. Even further useful features can be added to it, a few of which were mentioned above. As the next step, development of the model to allow a divergent-element geometry is near completion. Moreover, it is interesting and useful to compare the results of our models with those of other available optical Monte Carlo codes. We are currently working on performing such comparisons.