|Year : 2014 | Volume
| Issue : 1 | Page : 18-23
ScintSim1 : A new Monte Carlo simulation code for transport of optical photons in 2D arrays of scintillation detectors
Mohammad Amin Mosleh-Shirazi1, Zinat Zarrini-Monfared2, Sareh Karbasi3, Ali Zamani4
1 Physics Unit, Department of Radiotherapy and Oncology; Radiology Research Center, Shiraz University of Medical Sciences, Shiraz, Iran
2 Department of Medical Physics and Engineering; Student Research Committee, Shiraz University of Medical Sciences, Shiraz, Iran
3 Physics Unit, Department of Radiotherapy and Oncology, Shiraz University of Medical Sciences, Shiraz, Iran
4 Department of Medical Physics and Engineering, Shiraz University of Medical Sciences, Shiraz, Iran
|Date of Submission||16-Sep-2013|
|Date of Decision||14-Dec-2013|
|Date of Acceptance||15-Dec-2013|
|Date of Web Publication||20-Jan-2014|
Mohammad Amin Mosleh-Shirazi
Department of Radiotherapy and Oncology, Namazi Hospital, Shiraz University of Medical Sciences, Shiraz 71936-13311
Source of Support: None, Conflict of Interest: None
| Abstract|| |
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.
Keywords: Image-guided radiotherapy, Monte Carlo simulation, optical photons, scintillation detectors, segmented scintillators
|How to cite this article:|
Mosleh-Shirazi MA, Zarrini-Monfared Z, Karbasi S, Zamani A. ScintSim1 : A new Monte Carlo simulation code for transport of optical photons in 2D arrays of scintillation detectors. J Med Phys 2014;39:18-23
|How to cite this URL:|
Mosleh-Shirazi MA, Zarrini-Monfared Z, Karbasi S, Zamani A. ScintSim1 : A new Monte Carlo simulation code for transport of optical photons in 2D arrays of scintillation detectors. J Med Phys [serial online] 2014 [cited 2019 Dec 11];39:18-23. Available from: http://www.jmp.org.in/text.asp?2014/39/1/18/125481
| 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,  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. ,,
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.  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.  In addition, relatively low-dose MVCT projections could be taken using this device at frequent intervals and with a high speed.  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). ,
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. ,, 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. ,,
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. ,,, 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.  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.  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  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.
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 mm 3 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 10 7 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 (P ref ) is given by Fresnel equations. More details are available in reference. 
For an interface that separates two materials with indices of refraction n1 = 2.2 and n2 = 1, we obtained P ref 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.
|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|
Click here to view
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θ,  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 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θ).
|Figure 2: Distribution of the random numbers produced by the 'rand' command in MATLAB R2010a (107 calls binned into 8,000 intervals)|
Click here to view
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 (10 7 /8,000 = 1250; . Also, it can be seen from [Figure 2] that the distribution is symmetric around the mean.
- [Figure 3] shows P ref and Fresnel reflection coefficients (r perp and r para ) 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 P ref = 1 for θi > θc . For normal incidence, r para , r perp , and P ref have the expected values of −0.37, 0.37, and 0.14, respectively, which are also seen in [Figure 3].
- 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.
- 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 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 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|
Click here to view
|Figure 4: Spatial distribution of the points of incidence (on the exit face) for the transmitted rays|
Click here to view
|Figure 5: Distribution of the mean direction cosine of the transmitted rays against the position of the point of incidence on the exit face|
Click here to view
|Figure 6: Comparison of the simulated and expected polar angle probability densities of a ray reflected from a Lambertian surface|
Click here to view
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) |
Click here to view
| Discussion and Conclusion|| |
A number of groups have studied the use of thick-segmented crystalline scintillators as the convertor in MV imaging devices. ,,,,, 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. 
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. 
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,  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  or PHOTON  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,  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.
| Acknowledgments|| |
This work is part of a thesis submitted by Ms Zinat Zarrini-Monfared to Shiraz University of Medical Sciences for a degree of Master of Science in Medical Physics. The authors acknowledge the receipt of a research grant for this project from Shiraz University of Medical Sciences. Guidance from Professor W Swindell during the development of the previously-published initial single-element model and active participation and contribution of Professor PM Evans in its experimental validation are gratefully acknowledged.
| References|| |
|1.||Evans PM. Anatomical imaging for radiotherapy. Phys Med Biol 2008;53:R151-91. |
|2.||Boyer AL, Antonuk L, Fenster A, Van Herk M, Meertens H, Munro P, et al. A review of electronic portal imaging devices (EPIDs). Med Phys 1992;19:1-16. |
|3.||Mosleh-Shirazi MA, Evans PM, Swindell W, Webb S, Partridge M. A cone-beam megavoltage CT scanner for treatment verification in conformal radiotherapy. Radiother Oncol 1998;48:319-28. |
|4.||Munro P. Portal imaging technology: Past, present, and future. Semin Radiat Oncol 1995;5:115-33. |
|5.||Antonuk LE. Electronic portal imaging devices: A review and historical perspective of contemporary technologies and research. Phys Med Biol 2002;47:R31-65. |
|6.||Mosleh-Shirazi MA, Evans PM, Swindell W, Symonds-Tayler JR, Webb S, Partridge M. Rapid portal imaging with a high-efficiency, large field-of-view detector. Med Phys 1998;25:2333-46. |
|7.||Sawant A, Antonuk LE, El-Mohri Y, Zhao Q, Wang Y, Li Y, et al. Segmented crystalline scintillators: Empirical and theoretical investigation of a high quantum efficiency EPID based on an initial engineering prototype CsI (TI) detector. Med Phys 2006;33:1053-66. |
|8.||Wang Y, Antonuk LE, Zhao Q, El-Mohri Y, Perna L. High-DQE EPIDs based on thick, segmented BGO and CsI: Tl scintillators: Performance evaluation at extremely low dose. Med Phys 2009;36:5707-18. |
|9.||Sawant A, Antonuk LE, El-Mohri Y, Zhao Q, Li Y, Su Z, et al. Segmented crystalline scintillators: An initial investigation of high quantum efficiency detectors for megavoltage x-ray imaging. Med Phys 2005;32:3067-83. |
|10.||Wang Y, Antonuk LE, El-Mohri Y, Zhao Q, Sawant A, Du H. Monte Carlo investigations of megavoltage cone-beam CT using thick, segmented scintillating detectors for soft tissue visualization. Med Phys 2008;35:145-58. |
|11.||Liu L, Antonuk LE, Zhao Q, El-Mohri Y, Jiang H. Countering beam divergence effects with focused segmented scintillators for high DQE megavoltage active matrix imagers. Phys Med Biol 2012;57:5343-58. |
|12.||Mosleh-Shirazi MA, Swindell W, Evans PM. Optimization of the scintillation detector in a combined 3D megavoltage CT scanner and portal imager. Med Phys 1998;25:1880-90. |
|13.||Monajemi TT, Fallone BG, Rathee S. Thick, segmented CdWO4-photodiode detector for cone beam megavoltage CT: A Monte Carlo study of system design parameters. Med Phys 2006;33:4567-77. |
|14.||Wang Y, Antonuk LE, El-Mohri Y, Zhao Q. A Monte Carlo investigation of Swank noise for thick, segmented, crystalline scintillators for radiotherapy imaging. Med Phys 2009;36:3227-38. |
|15.||Cayouette F, Laurendeau D, Moisan C. DETECT2000: An improved Monte-Carlo simulator for the computer aided design of photon sensing devices. Proc. SPIE 4833, Applications of Photonic Technology 5, 2003:69-76. |
|16.||Tickner J, Roach G. PHOTON - An optical Monte Carlo code for simulating scintillation detector responses. Nucl Instr Methods Phys Res Sec B Beam Inter Mater Atoms 2007;263:149-55. |
|17.||Blake SJ, Vial P, Holloway L, Greer PB, McNamara AL, Kuncic Z. Characterization of optical transport effects on EPID dosimetry using Geant4. Med Phys 2013;40:041708. |
|18.||Yang X, Downie E, Farrell T, Peng H. Study of light transport inside scintillation crystals for PET detectors. Phys Med Biol 2013;58:2143-61. |
|19.||Briesmeister J. MCNP-a general Monte Carlo N-particle transport code, version 4C. LA-13709-M. Los Alamos: Los Alamos National Library; 2000. |
|20.||El-Mohri Y, Antonuk LE, Zhao Q, Choroszucha RB, Jiang H, Liu L. Low-dose megavoltage cone-beam CT imaging using thick, segmented scintillators. Phys Med Biol 2011;56:1509-27. |
|21.||Swindell W, Mosleh-Shirazi MA. Noise reduction by frame averaging: A numerical simulation for portal imaging systems. Med Phys 1995;22:1405-11. |
[Figure 1], [Figure 2], [Figure 3], [Figure 4], [Figure 5], [Figure 6]