Citation |

- Permanent Link:
- http://ufdc.ufl.edu/UFE0022617/00001
## Material Information- Title:
- 3D Deterministic Radiation Transport for Dose Computations in Clinical Procedures
- Creator:
- Al-Basheer, Ahmad
- Place of Publication:
- [Gainesville, Fla.]
- Publisher:
- University of Florida
- Publication Date:
- 2008
- Language:
- english
- Physical Description:
- 1 online resource (140 p.)
## Thesis/Dissertation Information- Degree:
- Doctorate ( Ph.D.)
- Degree Grantor:
- University of Florida
- Degree Disciplines:
- Nuclear Engineering Sciences
Nuclear and Radiological Engineering - Committee Chair:
- Sjoden, Glenn E.
- Committee Members:
- Bolch, Wesley E.
Haghighat, Alireza Dempsey, James F. Ditto, William L. Hawari, Ayman - Graduation Date:
- 8/9/2008
## Subjects- Subjects / Keywords:
- Charged particles ( jstor )
Dosage ( jstor ) Electrons ( jstor ) Energy ( jstor ) Libraries ( jstor ) Numerical quadratures ( jstor ) Photons ( jstor ) Physics ( jstor ) Radiation transport ( jstor ) Simulations ( jstor ) Nuclear and Radiological Engineering -- Dissertations, Academic -- UF - Genre:
- Electronic Thesis or Dissertation
bibliography ( marcgt ) theses ( marcgt ) Nuclear Engineering Sciences thesis, Ph.D.
## Notes- Abstract:
- The main goal of this dissertation was to establish the feasibility of basing megavoltage external photon beam absorbed dose calculations in voxelized phantoms on SN deterministic calculations and pre-calculated electron absorbed dose kernels derived from full-physics Monte Carlo. The SN derived electron absorbed dose kernel method EDK-SN, developed as part of this research, achieves total execution times that are on the order of several times to orders of magnitude faster than conventional full-physics Monte Carlo electron transport methods considering equivalently detailed models and data fidelity. With the rapid movement toward intensity modulated radiation therapy (IMRT), radiation beam intensities have increased dramatically over the past decade, thus heightening the need for further characterization of out-of-field organ absorbed doses, along with their associated biological risks. Assessment of these tissue absorbed doses is complicated by two fundamental limitations. First, anatomic information on the patient is generally restricted to a partial body CT image acquired for treatment planning; consequently, whole-body computational phantoms must be employed to provide the out-of-field anatomy model structure for absorbed dose evaluation. Second, existing methods based on Monte Carlo radiation transport, even with the application significant variance reduction, are quite computationally inefficient at large distances from the primary beam, and point-kernel methods do not properly handle tissue inhomogeneities. Moreover, since absorbed dose are generally tracked in all major organs in the body, variance reduction schemes for Monte Carlo are not all effective in this regard. The outcome of this dissertation is to demonstrate that absorbed dose from high-energy external beams radiation can be accurately computed for whole body and organ-specific absorbed doses. The EDK- SN method implements voxelized phantoms with discrete ordinates (SN) transport computations coupled with directional influences and (pre-computed) full-physics Monte Carlo based electron absorbed dose kernels to yield total body absorbed dose information. This research shows that the deterministic techniques coupled with Monte Carlo based electron absorbed dose-kernels has significant potential for organ absorbed dose evaluation in the clinical management of radiation therapy patients. ( en )
- General Note:
- In the series University of Florida Digital Collections.
- General Note:
- Includes vita.
- Bibliography:
- Includes bibliographical references.
- Source of Description:
- Description based on online resource; title from PDF title page.
- Source of Description:
- This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
- Thesis:
- Thesis (Ph.D.)--University of Florida, 2008.
- Local:
- Adviser: Sjoden, Glenn E.
- Electronic Access:
- RESTRICTED TO UF STUDENTS, STAFF, FACULTY, AND ON-CAMPUS USE UNTIL 2010-08-31
- Statement of Responsibility:
- by Ahmad Al-Basheer.
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- Copyright Ahmad Al-Basheer. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
- Embargo Date:
- 8/31/2010
- Classification:
- LD1780 2008 ( lcc )
## UFDC Membership |

Downloads |

## This item has the following downloads:
albasheer_a.pdf
albasheer_a_Page_094.txt albasheer_a_Page_006.txt albasheer_a_Page_071.txt albasheer_a_Page_105.txt albasheer_a_Page_114.txt albasheer_a_Page_033.txt albasheer_a_Page_035.txt albasheer_a_Page_078.txt albasheer_a_Page_140.txt albasheer_a_Page_054.txt albasheer_a_Page_088.txt albasheer_a_Page_098.txt albasheer_a_Page_095.txt albasheer_a_Page_118.txt albasheer_a_Page_135.txt albasheer_a_Page_134.txt albasheer_a_Page_121.txt albasheer_a_Page_128.txt albasheer_a_Page_109.txt albasheer_a_Page_110.txt albasheer_a_Page_100.txt albasheer_a_Page_016.txt albasheer_a_Page_115.txt albasheer_a_Page_139.txt albasheer_a_Page_108.txt albasheer_a_Page_082.txt albasheer_a_Page_043.txt albasheer_a_Page_013.txt albasheer_a_Page_099.txt albasheer_a_Page_089.txt albasheer_a_Page_065.txt albasheer_a_Page_024.txt albasheer_a_Page_138.txt albasheer_a_Page_047.txt albasheer_a_Page_067.txt albasheer_a_Page_132.txt albasheer_a_Page_034.txt albasheer_a_Page_036.txt albasheer_a_Page_106.txt albasheer_a_Page_066.txt albasheer_a_Page_027.txt albasheer_a_Page_049.txt albasheer_a_Page_112.txt albasheer_a_Page_130.txt albasheer_a_Page_005.txt albasheer_a_Page_028.txt albasheer_a_Page_050.txt albasheer_a_Page_057.txt albasheer_a_Page_096.txt albasheer_a_Page_092.txt albasheer_a_Page_125.txt albasheer_a_Page_040.txt albasheer_a_Page_022.txt albasheer_a_Page_002.txt albasheer_a_Page_084.txt albasheer_a_Page_081.txt albasheer_a_Page_058.txt albasheer_a_Page_103.txt albasheer_a_pdf.txt albasheer_a_Page_012.txt albasheer_a_Page_127.txt albasheer_a_Page_123.txt albasheer_a_Page_042.txt albasheer_a_Page_053.txt albasheer_a_Page_003.txt albasheer_a_Page_018.txt albasheer_a_Page_020.txt albasheer_a_Page_137.txt albasheer_a_Page_102.txt albasheer_a_Page_097.txt albasheer_a_Page_079.txt albasheer_a_Page_080.txt albasheer_a_Page_052.txt albasheer_a_Page_019.txt albasheer_a_Page_031.txt albasheer_a_Page_090.txt albasheer_a_Page_061.txt albasheer_a_Page_129.txt albasheer_a_Page_087.txt albasheer_a_Page_032.txt albasheer_a_Page_004.txt albasheer_a_Page_062.txt albasheer_a_Page_122.txt albasheer_a_Page_124.txt albasheer_a_Page_025.txt albasheer_a_Page_133.txt albasheer_a_Page_011.txt albasheer_a_Page_039.txt albasheer_a_Page_046.txt albasheer_a_Page_073.txt albasheer_a_Page_117.txt albasheer_a_Page_001.txt albasheer_a_Page_008.txt albasheer_a_Page_060.txt albasheer_a_Page_064.txt albasheer_a_Page_059.txt albasheer_a_Page_030.txt albasheer_a_Page_023.txt albasheer_a_Page_051.txt albasheer_a_Page_045.txt albasheer_a_Page_116.txt albasheer_a_Page_029.txt albasheer_a_Page_017.txt albasheer_a_Page_131.txt albasheer_a_Page_041.txt albasheer_a_Page_093.txt albasheer_a_Page_091.txt albasheer_a_Page_037.txt albasheer_a_Page_069.txt albasheer_a_Page_074.txt albasheer_a_Page_136.txt albasheer_a_Page_009.txt albasheer_a_Page_126.txt albasheer_a_Page_086.txt albasheer_a_Page_113.txt albasheer_a_Page_083.txt albasheer_a_Page_044.txt albasheer_a_Page_010.txt albasheer_a_Page_072.txt albasheer_a_Page_014.txt albasheer_a_Page_070.txt albasheer_a_Page_068.txt albasheer_a_Page_038.txt albasheer_a_Page_048.txt albasheer_a_Page_119.txt albasheer_a_Page_063.txt albasheer_a_Page_056.txt albasheer_a_Page_026.txt albasheer_a_Page_077.txt albasheer_a_Page_120.txt albasheer_a_Page_055.txt albasheer_a_Page_085.txt albasheer_a_Page_101.txt albasheer_a_Page_111.txt albasheer_a_Page_015.txt albasheer_a_Page_075.txt albasheer_a_Page_076.txt albasheer_a_Page_007.txt albasheer_a_Page_107.txt albasheer_a_Page_021.txt albasheer_a_Page_104.txt |

Full Text |

-IK~~ rO22 (+ 2 -a 2)/e+(1/a2 2)+(1 +2a*)e2/) (A-23) E = (A-24) 1 + a(1- pu) pu= cos9 ->~dpu= sin Bd6 (A-25) dpu = dE (A-26) hv E = (A-27) hv From this example one can note that; the produced electron in a certain direction has a unique energy for a giving initial photon. So for an angle aGZ aroundiscretized range of angles we can calculate the average energy of the electrons produced in th Pair production Figure illustrates schematically a pair production event in the nucleus electric field. The incident photon gives up all of its energy in the creation of an electron positron pair with kinetic energies Tand T The energy conservation equation, ignoring the vanishing small kinetic energy giving to the nucleus is simply hv = 2moc2 + + T (A-28) = 1.022MIeVy + T- + T+ The average electron/positron average energy is -v h-1.022M2eVy T = (A-29) transport paths relative to a spatially discretized mesh grid interval, 3-D SN methods can be implemented to directly yield very accurate dose distributions. Therefore, the deterministically derived dose using a kerma approximation will be accurate up to the range of electron path lengths comparable to ranges smaller than the photon spatial SN mesh grid interval in anatomical models (typically representative of the anatomical data voxel size). As mentioned earlier, at higher photon energies, where electron ranges exceed the SN mesh grid interval, this is not the case, and interactions yielding energetic electrons must account for the mechanism of electron transport through the problem phase space to deliver an absorbed dose after representative electron transport distances. This work will set the foundation for state of the art 3-D absorbed dose computations for high energy external beams to account for interactions yielding energetic electrons. The EDK-SN methodology will account for electron transport in the phantom by applying a pre-computed electron absorbed dose kernel calculated to a very low variance It is worth comparing the EDK-SN methodology, proposed here, with the traditional convolution-superposition technique, which applies a database of "energy deposition kernels," calculated using Monte Carlo techniques to determine a absorbed dose29,30. We note that these kernels can be used to describe the spread of energy about the primary photon's point of interaction in the phantom, however, these have limited use since the convolution-superposition technique is performed without regard to the incident photon angular component. While, in contrast, the EDK-SN procedure is based on a detailed, full physics Boltzmann transport solution with direction and flux information directly available for coupling with the EDKs along a propagated photon direction of travel. 106 -1.374 106 -1.374 106 -1.374 106 -1.374 106 -1.374 106 -1.374 106 -1.374 105 -1.415 105 -3.00 1 -0.001205 28 -0.2995 -999 -999 -999 -999 -999 -999 -999 -999 -999 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 u=95 u=96 u=97 u=98 u=99 u=100 u=101 u=102 u=103 vol=1.60E+01 $Radii-upper shaft vol=1.70E+01 $Radii-lower shaft vol=2.70E+01 $Radii-distal vol=4.90E+01 $Ulnae-proximal vol=2.70E+01 $Ulnae-upper shaft vol=1.80E+01 $Ulnae-lower shaft vol=7.00E+00 $Ulnae-distal vol=8.60E+01 $Hand vol=1.30E+01 $Teeth vol= 78406.867 $ Air vol= 10210.44 $ Lung -999 imp:p,e=1 u=301 -999 imp:p,e=1 u=302 303 2 -1.00 -999 imp:p,e=1 u=303 vol= 27482.027 $ Tissue 304 101 -1.32 -999 imp:p,e=1 u=304 vol= 2304.37 $ Bone C ******************* Lattice definition ************************************ 1001 0 -100 fill=10000 imp:p,e=1 $ surrounding box 555 0 -200 lat=1 u=10000 imp:p,e=1 fill = 0:59 0:26 0:166 C *************** Image data start from here ******************************* Phantom Input from Ghoast-3-D C ************************************************************* c Surface Cards C ********************************************************* 100 rpp 0 60 0 27 0 167 $(-29:30 -13:13 -83:83) 200 rpp 0 1.0 01.0 0 1.0 999 rpp -25 85 -25 52 -25 192 c *************************************************************** Material Cards C*************************************************************** mode pe sdef x=d1 y=d2 z=d3 erg=d4 AXS= 0 1 0 dir= 0.99763 vec= 0 1 0 sil H 30 40 spl D 0 1 si2 H 0.0 1.0 sp2 D 0 1 si3 H 119 136 sp3 D 0 1 si4 H 0.01 8.0 sp4 D 0.0 1.0 c tallies cards phys:p 4j 1 ctme 500 totnu based on assuming charged particle equilibrium conditions, not suitable for high energy beams. Moreover, the statistical uncertainty associated with each mesh tally site for such computations would require extended periods of running time, especially off-axis from the source and outside the region covered by the primary beam to converge to acceptable levels. These factors limit the Monte Carlo calculations from being fuuly adequate in practical applications where whole body and organ-specific absorbed dose information is needed. Notably, an equivalent Monte Carlo calculation to obtain the flux distribution with statistically reliable results over the entire problem space may require orders of magnitude longer running times compared to parallel deterministic calculations. When comparing *F8 (p mode) transport and *F8 (p, e mode) transport, (Figure 4-5) shows significant differences in the outcomes for the same simple model. This directly results from accessing different calculation modes in MCNP5 with regard to electron physics. If electron transport is turned on (mode p, e), then all photon collisions except coherent scatter lead to electrons that are banked for later transport. Alternatively, if electron transport is turned off (achieved by omitting the 'e' on the MCNP5 Mode card), then a "thick-target bremsstrahlung" model (TTB) is used. This model generates electrons, but assumes that they are locally "slowed to rest". The TTB production model contains many approximations compared to models used in actual electron transport. The choice of which tally is to be used and the detailed physics treatment by which the Monte Carlo code is executed is essential to the success of the EDK-SN method. In conducting computations for generation of the EDKs, we used the full detailed physics treatment, implementing a *F8 (p, e mode) tally, as opposed to the alternative methods that lead to less robust electron physics, essentially employing a collisional kerma. A detail discussion on photons interactions in matter is presented in Appendix A. 4.4 EDK -SN Code System Methodology 4.4.1 Net Current For any SN calculation, the PENTRAN discrete ordinates code preserves angular information explicitly in parallel data storage arrays. This allows one to calculate Jser and to apply the pre- computed electron absorbed dose kernel along the net current direction, Figure 4-6, where: J =J i^+J 7+J. k net nx ny =(4-3) J, = J J-, (4-4.a) LIST OF TABLES Table page 1-1 Comparison between Monte Carlo methods and deterministic techniques.............._.. ....28 1-2 Approximate attenuation of photons with a layer of water equal to the maximum range of secondary charged particles28 .........._.... ........_....._. ............3 1-3 Approximate thickness of water needed to achieve various attenuation amounts compared to maximum energy electron ranges generated by the same beam ...................32 2-1 Total number of directions, cumulative problem time** required for a single energy group, P3 scattering, 126,000 fine mesh cells using 12 processors on a parallel cluster ................ ...............54..__.._....... 3-1 Computed tomography image sources for the development of the UF pediatric phantom series .............. ...............60.... 3-2 Total number of directions, Cumulative Problem Time required for different energy groups, P3 Scattering, and 189k fine mesh cells............... ...............76. 5-1 Comparison of selected organ absorbed dose rate (MeV/g. s) calculated using MCNP5 pulse height tally with (photon, electron mode) and EDK-SN for the UF hybrid 15-year-old male phantoms for a chest 8 MV X-ray flat weighted source ..........115 5-2 Prostate absorbed dose rate (MeV/g.s) calculated for the test problem using EDK-SN and MCNP5 using pulse height tally (*F8) with (photon, electron mode) at different MCNP5 Monte Carlo running times .....__._.. ... ..._.... ...............116. 5-2 Absorbed dose rate (MeV/g. s) calculated for the test problem using EDK-SN and MCNP5 using pulse height tallies (*F8) with (photon, electron mode) .......................117 5-3 Absorbed dose rate (MeV/g. s) calculated for the test problem using MCNP5 using pulse height tally *F8 with (photon, electron mode) and MCNP5 energy deposited tallies (F6) ........._ .......... ...............118.... A-1 Ratios of photon interactions in water for range of energies ..............._ ............ .....124 2.2.2 Monte Carlo Model The MCNP5 Monte Carlo model geometry was defined to be equivalent to that used in the PENTRAN model. The 60C0 source was distributed isotropically over multiple equi-probable sources cells in one energy group. Source spectra were defined in a consistent manner with the various multi-group libraries considered. Volumetric cell flux F4 mesh-tallies and energy deposited F6 tallies were used in the Monte Carlo geometry description that was equivalent to the discretized SN volumes. The MCNP5 F6 tallies were used along the central axis of the model across the lead box, air region, and water phantom. All simulations were performed using a PC cluster composed of 16 Dual Intel Xeon processors at 2.4 GHz with 2 Gb/processor, for a total memory of 32 Gb. MCNP5 tallies along the central axis were obtained using variance reduction techniques. The cell fluxes converged, on average, to within 3 % for the reference case. 2.3 Comparison of Deterministic and Monte Carlo Results Calculations were performed to compare the SN and Monte Carlo results for the scalar fluxes along the model central axis. Figure 2-2 depicts a 3-D scalar flux distribution computed by PENTRAN for the reference case for different energy groups. The highest energy group (group 1) is forward peaked with some ray-effects along the central axis. However, as the radiation down scattered to lower energy groups, it demonstrates more isotropic behavior evident in energy groups 2 and 3. This simulation was based on the use of an S42 angular Legendre- Chebychev quadrature (1848 directions/mesh) with P3 Scattering anisotropy using the CEPXS library. This highlights the fidelity of the information available in a converged deterministic SN computation. 2.3.1 Investigation of Angular Quadrature Medical physics applications represented in SN models often require small mesh sizes, relatively large models, and low density materials. These problems can lead to severe ray-effects, 1 -(x2 2) f (x, y, z) = 2exp( 2 )(1-2) where o, is the angular spread about the axis including the penumbra. The on-axis term may be expressed in terms of its relationship to the depth dose in matter, go(z4f) g~z)= gozz (1-3) where SSD is the source-to-surface distance, and ze, is the effective depth in water, assuming the linear stopping power of the material at depth : is relatively independent of electron energy for normal body tissue. The quantity go(z4f) can be related directly to a depth dose curve measured in a water phantom, assuming a uniform incident beam, by the Eq. 1-4 go (z)I~xX (x x )2 + y_ y 2 L d 14 Do (x, y, z) = 2 x x'y 14 where the limits of the integral depend on the dispersion of the field at depth z. An example of pencil method used for clinical photon radiotherapy is an application of the Ahnesjo pencil beam model6. In this method, combination of measured data and Monte Carlo calculated convolution kernels are used to characterize clinical photon beams and eventually calculate the absorbed dose. The strength of the Ahnesjo pencil beam method is that it is fast enough for clinical use, on the order of minutes. It also uses measured data for input that may be specific to the treatment unit. However, because this method applies simple corrections using pencil beams in homogenous matter, it is not accurate as a complete coupled charged-particle- photon Monte Carlo simulation in heterogeneous media. Moreover, it does not fully represent the directional component of the scattered radiation. will be scattered into the solid angle daZ in the direction 6, B and being related by Eq. (A- 8), we have that (A-9) daZ daZ daZ daZ d~ sin0 d8 (A-10) d~Ze de sin 8e dBe cot 8e = (1 +a) tan( ) (A-1 1) -1 1 = (1+ a) dB (A-12) sin2 8e COS2 8/2 sin 8e = 1l+(1+a)2 tan208/2 1/2 (A-13) e = 27r (A-14) de = -dj > = -1 (A-15) daZ 1 (1 +cos8) sin 8 (A-16) d~Ze (1+ a) sin3 8e for 0.511 MeV photons, a = 1 S= 2 sin(8/ 2) 1+3 sin2 (6 /32) (A-17) dae, d~Ze 2 hv hv h v-8 hv hv -= -cs+ p cos 8e (A-2) c c hv sin B = pc sin 8e (A-3) pc can be written in terms of T by invoking the "law of invariance": p~c= T(T+moc2 (A-4) In which mo is the electron's rest mass. As a result of solving these equations algebraically, one can derive the following Compton's relations hv hv = (A-5) 1+(hvlm0C2)(1-COS 8) co~ 1hymC)a (A-6) 21 T, = hv hv = hv 1-+2 i2 Bi2 (A-7) To obtain the fraction of the scattered photons in a given direction, Klein and Nishina (reference) have carried out a quantum-mechanical treatment of the problem using the Dirac equation for the electron and have obtained the equation derKN 02 hy'hv hv . ~ -v +-si2 (A-8) da 2 hv hv h Equations similar to those for the Compton photon distribution can be obtained for the Compton electron distribution. Since the probability that an electron will be scattered into the solid angle daZ situated in the direction 8e is the same as the probability that a primary quantum CHAPTER 4 ELECTRON ABSORBED DOSE KERNELS TO ACCOUNT FOR SECONDARY PARTICLE TRANSPORT INT DETERMINISTIC SIMULATIONS In previous chapters, Charged Particle Equilibrium (CPE) and its limitations were presented, and that in cases with high energy photons, corrections for secondary electron transport may be significant. As mentioned, this is readily treated in Monte Carlo codes, yet quite difficult to treat explicitly in deterministic codes due to the large electron ranges and added numerical complexities in reaching convergence in photon-electron transport problems The goal here is to produce a methodology to accurately estimate organ absorbed doses and attribute whole body absorbed doses incurred from irradiations characteristic of external beam therapy using the CT-based anatomical patient geometries in 3-D deterministic discrete ordinates (SN) radiation transport models. Accurate absorbed doses will be rendered using pre-computed, high resolution Monte Carlo based electron kernels. In doing so, the flux and the current provided by the SN method will be used to both scale and project the absorbed dose using the Monte Carlo based kernels, and map it to the surrounding voxels. Hence, this approach preserves high fidelity and accuracy, accounting for complete photon transport, with the final absorbed dose delivered from electron transport propagated by coupling with 3-D deterministic photon transport; this is what makes this approach novel. In this chapter, water phantom results using a 0 to 8-MeV step uniform beam will demonstrate that absorbed doses can be accurately obtained within the uncertainties of a full physics Monte Carlo simulation. 4.1 Introduction to The EDK-SN Method Various modalities are applied in absorbed dose calculations for high energy photon beams, including empirical models that are based on measurements made in the clinic, Of theoretically based models, such as the application of Monte Carlo techniques3'48. Alternatively, as mentioned earlier, deterministic methods have been implemented to assess absorbed dose where 'F is the eonegy fluece~, and (ie": is the mass energy absorption coefficient for 7 ray energy E in medium Z. This result is very important, since it allows us to calculate the measurable dose D based on the calculable quantity Kc - 1.5.1 Deterministic Dose Calculations for Diagnostic Application For low energy beams, such as the one used in diagnostic applications, the concept of charged particle equilibrium is of a great importance as a mean of equating absorbed dose D to collisional kerma Kc - For a volume V, containing a smaller volume v, irradiated with an indirectly ionizing radiation field (i.e., Photons and Neutron fields), there are four basic causes for CPE failure in volume v: In-homogeneity of the atomic composition within volume V. In-homogeneity of density in V. Non-uniformity of the field on indirectly ionizing radiation in V. Presence of a non- homogeneous electric or magnetic field in V. These four situations will disrupt CPE by preventing charged particles of a given type and energy leaving the volume v to be replaced by an identical particle of the same energy entering27 In practice, CPE conditions are met if the attenuation of the primary beam is negligible in the maximum range of the secondary electrons generated by the beam. As indicated by Table (1- 2) for a layer of water of 5 cm thickness, this scenario is plausible in the limit of diagnostic medical physics applications, where the energies fall within 10 to 200 keV, and CPE conditions are met for water/tissue volumes less than 5x5x5 cm3, and consequently D is equal to Kc under such conditions. 27F. H. Attix, "Energy imparted, energy transferred and net energy transferred," Phys Med Biol 28, 1385-1390 (1983).44 28N. I. o. S. a. T. NIST, ESTAR(Stopping Power and Range Tables for Electrons ), http://physics.nist.gov/PhysRefData/Star/Tx/SARhm.1 29T. R. Mackle, A. F. Bielajew, D. W. Rogers, and J. J. Battista, "Generation of photon energy deposition kernels using the EGS Monte Carlo code," Phys Med Biol 33, 1-20 (1988).54 30A. Ahnesjo, P. Andreo, and A. Brahme, "Calculation and application of point spread functions for treatment planning with high energy photon beams," Acta Oncol 26, 49-56 (1987).53 31L. A. N. L. X-5 Monte Carlo Team, M~CNP5- A General Monte Calrlo N-Particle transpotr~r~r~t~rt~tr~rt~ code, Version 5, 10 32J. E. White, J. D. Intgersoll, C. Slater, and R. Roussin, BUGLE-96 : A Revised2\~ulti-group Cross Section Library for L WR Applications Based on ENDF B- y7, 46 33L. J. Lorence, J. E. Morel, and G. D. Valdez, Physics Guide to CEPXS: A naultigroup coupled electron-photon cross section generating code, 11 34G. SJODEN, and A. HAGHIGHAT, PENTRAN: Parallel Environment Neutral-Particle 7RA szort Sn in 3-D Cartesian Geonsetry 35C. Yi, and A. Haghighat, PENM~SHExpress, 62 36G. Longoni, "Advanced quadrature sets, acceleration and preconditioning techniques for the discrete ordinates method in parallel computing environments," (University of Florida, Gainesville, FL, 2004), pp. 47 37B. Petrovic, and A. Haghighat, "New directional theta-weighted (DTW) differencing scheme and reduction of estimated pressure vessel fluence uncertainty," in the 9th hIternational Synaposiunt on Reactor Dosintetry, edited by Prague, Czech Republic, 1996), pp. 746-753.48 38G. Sjoden, and A. Haghighat, "The exponential directional weighted (EDW) differencing scheme in 3-D cartesian geometry," in The Joint hIternational Conference on Ma thema'lltliL and Superconsputing for Nuclear Applications, edited by Saratoga Springs, New York, 1997 ), pp. 21267-21276.49 39G. Sj oden, "An Efficient Exponential Directional Iterative Differencing Scheme for Three- Dimensional Sn Computations in XYZ Geometry," Nuc. Science and Eng 155, 179-189 (2007). 12 40G. Sj oden, and A. Haghighat, "Taylor Proj section Mesh Coupling Between 3-D Discontinuous Grids for Sn," Transactions of the American Nuclear Society 74, 178-179 (1996).77 this work focused on the use of 16 energy groups for a flat weighted 8 MeV high energy source term, the parallel EDK-SN computation could be parsed to a maximum of 16 parallel processors. Near linear speedups were achieved. The parallel code follows the following steps in executing the EDK-SN methodology for total dose calculation: Initiate the code on all processors for sets of energy groups, up to one group per processor Sort the locally parallel-allocated fine mesh distributions into a global monotone increasing mesh distributions Create Energy group dose matrix aliased to the location of the DDVs Translate material tags into densities for later calculation Project the Group-dependent EDK based on the group current vector and material at the DDV Scale group dose using group scalar flux, correcting for densities at material interfaces Use MPI-REDUCE to sum all doses for the phase space across the parallel processing grid 4.8 Discussion and Findings We examined a new EDK-SN methodology of coupling pre-computed Monte Carlo based electron absorbed dose kernels (EDKs) with the discrete ordinates (SN) solutions to achieve accurate absorbed doses for megavolt energy photon scenarios, since charged particle equilibrium fails and full photon-electron physics must be accounted for to fully attribute absorbed dose. The EDK-SN methodology was compared with the Monte Carlo method and shown to give results consistent within the statistical uncertainty of the absorbed dose computed using Monte Carlo calculations. By following the EDK-SN approach, the absorbed dose applied to tissue from any incident external high energy photon beam (pencil, fan, areal, etc) can therefore be readily determined based on a coupling of the deterministic SN photon transport solution from the PENTRAN code using the EDK-SN methodology. As a result, this research renders the first computation times more tractable, including voxel down-sampling and parallel computing resources. 1.4.2 Deterministic Methods Several different methods have been developed for numerically solving the Boltzmann equation" ~l. In general these methods depend on discretizing the energy, angular, and spatial domains of the Boltzmann transport equation. 1.4.2.1 Discrete ordinates method (SN method) The most common deterministic approach currently used for photon transport is the discrete ordinates method (SN method) developed in detail by Carlson over 50 years agol6. The discrete ordinates equations are derived by discretizing the BTE over phase space elements (i.e., spatial-energy-direction variables) to construct discrete-variable cell particle balances. Equation (1-7) is the standard Legendre expanded Cartesian multigroup integro-differential formulation of the BTE, and x, y, z are the spatial variables along each axis in the Cartesian coordinates. g= =1 =0 k=1 -I+k [ 1~ r(x,y,z)cos(kp)+0b' (x,y,z)sin(kp) )+q (x,y,z,pi,p) (1-7) Other variables in Equation (1-7) are: pu, r, correspond to the direction cosines along each axis for each ordinate Y, is the energy group angular flux o-, is the total group macroscopic cross section 0-, a,,i~ is the lth moment macroscopic scattering cross section for group g to g Figure 3-4. UF-Series B phantom (11 years old male) with (398 x 242 x 252 voxels); the corresponding PENTRAN "back-thinned" models with (132 x 80 x 84 voxels) and 72 materials (upper right), and another with (79 x 48 x 50 voxels) and 66 materials (lower right) into a 3-D spatial grid distribution and sub-divided into coarse mesh grids, each containing a corresponding number of Eine mesh cells. 3.3 PENTRAN-MP Code System Currently, the PENTRAN-MP code system is composed of 3 levels of calculation, including pre-processing algorithms that include the GHOST-3-D phantom voxel collapsing (back-thinning) code, the DXS x-ray source generation code 45, the PENMSH-XP 3-D mesh distribution and input generator code, and the CEPXS code (a product of Sandia National Laboratories) for generating multigroup cross sections 33. The radiation transport calculations are carried out in PENTRAN, and post-processed using PENDATA to gather parallel distributed data, and 3-D-DOSE to render user specified absorbed dose data. 3.3.1 Pre-processing Codes 3.3.1.1 The collapsing code (GHOST-3-D) The GHOST-3-D, developed by Ghita and Al-Basheer, is a code used to input CT-scan binary data obtained from a high resolution voxelized human phantom. There are several constraints in simulating mega-voxel models, including balancing limited computer memory (even with the benefit of parallel memory storage, as in PENTRAN) with proper discretization, etc. To optimize accuracy and computational work in representing the voxelized phantom, we developed a code, GHOST-3-D, that provides an equivalent model with a reduced number of meshes (subj ect to user-derived accuracy targets), overcoming the mentioned difficulties. 3.3.1.2 The mesh generator code (PENMSH-XP) PENMSH-XP35 by Yi and Haghighat is a 3-D Cartesian-based mesh generator that prepares discretized material and source distributions for the PENTRAN particle transport code. PENMSH-XP is capable of performing two different methodologies in generating PENTRAN input files. In the first methodology, a physical model is partitioned into slices along the z-axis, storage and time limitations is also a must. This issue has been investigated, particularly in the low density portions of the model (mostly in air). Finally, we compare the performance achieved with three cross section libraries available from the nuclear engineering and physics communities. These multi-group gamma-ray libraries were derived from the BUGLE- 96321ibrary, the CEPXS package33, and the SCALES package. 2.1 Code Systems PENTRAN (Parallel Environment Neutral particle Transport) is a 3-D discrete ordinates code that solves the Boltzmann transport equation on distributed parallel computers. It is capable of complete, hybrid phase space decomposition with up to a 97% parallel fractionl9,34 I executing a parallel transport calculation on each distributed node, PENTRAN performs discrete ordinates sweeps only on the local portion of the phase space assigned. As problem sizes increase, more machines may be needed to scale up to accommodate the total parallel memory demand. The inherent parallel scalability of PENTRAN with advanced numerical treatments (including block adaptive mesh structures and adaptive differencing schemes) makes it ideal to apply to challenging problems in medical physics. The PENTRAN code is part of a code system that includes pre- and post-processors, and other utilities to manage large parallel datasets. For the purpose of comparison to the Monte Carlo methods, the well known MCNP5 Monte Carlo code developed by Los Alamos National Laboratory is used to simulate a 60C0 beam for an equivalent unit. MCNP5 is used to calculate the scalar flux along the central axis of the model. For these studies, cell flux mesh-tally (F4) and energy deposition volume tallies (F6) are defined in the Monte Carlo model in a way that is analogous to the Cartesian fine mesh cells used in the PENTRAN. Figure 3-10. Simulation methodology for 3-D absorbed dose computations using PENTRAN- MP code system. Left: UF series B voxelized phantom, middle: correspondent PENTRAN input, right: 3-D absorbed dose distribution after irradiation with a x-ray source Carlo calculations to the right lung. This behavior can be explained by the ray-effects in the low density lung tissue (see Chapter 2). Since the X-ray source is highly directional, and the organ (right lung) is of a low density material, and little primary radiation and very low scattered radiation is streaming through the right lung, making conditions ideal for inducing ray effects in SN computations for the photon component. This is mitigated by using high order quadratures, and can reduce absolute differences; these were reduced to 4% (Table 5-2). However, the increase of the quadrature level costs additional running time, with an increase from 2 hours to 9 hours. This increase in quadrature level was necessary for large low density volumes that experience low levels of scattered radiation. Another method, ordinate splitting, allows for a biased quadrature set to focus on the "important" directions, although the technique was not used here. Table 5-2. Absorbed dose rate (MeV/g. s) calculated for the test problem using EDK-SN and MCNP5 using pulse height tallies (*F8) with (photon, electron mode) Organ MC(*F8) 2o--MC EDK-SH (MC-EDK- (p, e mode) Uncertainty SN)/EDK-SN Right Lung-(Sn32) 5.62E-06 2.0% 6.96E-06 19.33% Right Lung-(Sn62) 5.62E-06 2.0%/ 5.39E-06 4.03% Left Lung 2.70E-04 0.40% 2.79E-04 3.0% Liver 1.00E-05 3.0% 1.02E-05 1.54% Pancreas 1.11E-06 4.0% 1.17E-06 5.1% 5.2.3 MCNP5 Energy Deposited Tally (F6) and Pulse Height Tally (*F8) Results When comparing MCNP5 organ absorbed dose results, for the 15-year old phantom, calculated using energy deposited tally (F6) and pulse height tally (*F8), we noted inconsistencies in the results of the two tallies (Table 5-3). This is in keeping with our earlier investigation of the two tallies earlier in Chapter 4. This inconsistency in the F6 tally suggests that full physics Monte Carlo calculations are essential for accurate absorbed dose estimation in high energy external beam absorbed dose calculations. To achieve maximum speedup with the EDK-SN methodology, parallel programming methods were implemented in the EDK-SN code using energy group decomposition via the Message Passing Interface (MPI)49 library. The MPI library is meant to provide an essential virtual topology, process synchronization, and communication functionality between a set of processes (that have been mapped to nodes/servers/computers) in a language independent way, with library specific syntax (bindings), plus a few features that are library specific. The EDK-SN code has been designed in ANSI FORTRAN-90 (to take advantage of dynamic array allocation) and parallelized for operation on distributed memory, multiple instruction/ multiple data (MIMD) machine architectures using initial standard Message Passing Interface protocols (MPI-1). All Input/ Output (I/O) is performed by each processor in parallel (as required) to include input processing, initialization, processing, and Eile output. The EDK-SN code carries out the computations with electron dose kernels based on the scalar flux scaling for each energy group independently Consequently, the decomposition used is energy domain decomposition. The basic philosophy of this approach is to decompose energy groups to each processor, and carry the total angular (net current vector) and spatial domain information on each processor. The main advantages of the EDK-SN parallel algorithm is that appreciable speed up can be obtained to determine the dose in any voxel throughout the spatial domain. Since little communication is required between the processors, n processors would solve the problem nearly n-times faster than a single unit. Energy group absorbed doses are partitioned for the entire problem, and partitioned energy group assignments are scaled to n processors for coupling photon transport information with EDKs to produce group absorbed doses over the spatial domain. An MPI REDUCE call is made to collect all absorbed doses across processors. Since 5.1.3 EDK-SN Absorbed dose Comparison The SN calculation was performed based on a full photon transport calculation using deterministic techniques in the heterogeneous phantom. In principle, this eliminates any need for any correction for tissue heterogeneities which affect the different scattering components in differing ways. Moreover, the EDK-SN calculation was performed for material specific absorbed doses. This allowed the use of the pre-computed electron absorbed dose kernels for the different densities present in the phantom model. Pulse height tail using photon-electron mode (SBLLBSS) 4.00E-04 3.50E-04 -+- MC-Dose 3 3.00E-04 n-> Myaterial -EDK-SN2. 2.50E-04 -r 1.50E-04 1.5 S1.00E-04 1. 5.00E-05 -I I a9 0.00E+00 0 10 20 30 40 50 X(cm) Figure 5-6. EDK-SN heterogeneous SBLLBSS phantom total absorbed dose rate distributions along the X-axis compared to MCNP *F8 tally. MCNP uncertainty (20) average (8.0%) on-axis. Average absolute relative difference (4.1%) For the phantoms used in this study, calculations were performed to compare the EDK-SN and Monte Carlo pulse height tally results for the absorbed dose deposited along the model central axis. As shown in the Figures 5-6 and 5-7, the EDK-SN results are in very good agreement with the MCNP5 results, and differences between the two are within the statistical uncertainty of the Monte Carlo calculation. The Monte Carlo uncertainties, on the central axis of Table A-1. Ratios of photon interactions in water for range of energies Photon Photon Energy (MeV) Inercton 2.0 3.0 6.0 10.0 14.0 18.0 20.0 % of Incoherent 99.21 % 97.15 % 88.59 % 77.06 % 67.32 % 59.44 % 56.04 % scattering % of Nuclear 0.79 % 2.82 % 10.78 % 21.18 % 29.76 % 36.69 % 39.64 % P.P % of electron 0.0 % 0.03 % 0.61 % 1.76 % 2.9 % 3.90 % 4.35 % P.P APPENDIX A HIGH ENERGY PHOTON TREATMENT (IN CONDENSED ELECTRON TRANSPORT) The photons energy used in radiotherapy treatment planning are 20 MeV or less. Photon interactions that take place at these high energies predominantly Compton Effect (C.E.) or/and Pair Production interactions by different percentages, as shown in Figure A-1 and Table A-1. Att. Coefficient WATER (cm2/g) + Scatt. Coh.(cm2/g) +e Scatt. Incoh.(cm2/g) P.E.A(cm2/g) P.P Nuc. Field (cm2/g) +r P.P Elec. Field (cm2/g) + Total w/ coh +e Total w/o Coh. 1.00E+04 1.00E+02 S1.00E+00 E o 1.00E-04 1.00E-06 1.00E-08 1.00E-03 1.00E-02 1.00E-01 1.00E+00 Energy(MeV) 1.00E+01 1.00E+02 Figure A-1. Total mass attenuation coefficient for photons in water, indicating the contributions associated with phoelectric absorption, Compton scattering and electron-positron pair production G7(20 30)key G6(30-40)key 7.51E+06 6.01E+06 4.51E+06*PETA S3.01E+06 ' 1.51E+06 1.00E+04 0 10 20 30 Y(cm) 2.00E+05 G 8(10-20)key 1.50E+05 .MCNP5 -PENTRAN 9.00E+05 i5.00E+04 ; 0.00E+0 6 .. 0 5 10 Y~cm)15 20 25 Y=166cm X=13Decm 1.41E+06 1,21EI06 4.10E*05 is 5n 1 : IleG: 0i 00E*0 2.04140E0 .0E0 .0*110E Z, .cm)G Z=43.5cm Xr13.0cm 1.01E407 0.10E+06 , 8.10E*06 7.10E+06*Mc8s 3.10E+0 a9 G, 000000t 2.0E0140E0160E 1 0E+1 .0E+2 Z(cm) G5(40-50)key 6.01E+06 5 01E+06 3.01E+06 201E+06 1 01E+06 * MONPS, * NPENTAN O 5 10 Y~cm) 15 20 25 -M,;Npf * PENTRAN i1.76E+06 - 1.00E+04 ; "**** O 5 10 Y(cm) 15 20 25 Figure 3-8. Angular flux distributions using S32 with P3 scattering anisotropy in PENTRAN using the CEPXS library vs. MCNP5 Monte Carlo results for energy groups 5, 6, 7 and 8 along the source central axis. Monte Carlo uncertainties were less than 2%, 3%, 6%, and 30% respectively Figure 3-9. Angular flux distributions using S32 with P3 scattering anisotropy in PENTRAN using the CEPXS library vs. MCNP5 Monte Carlo results for energy groups 1, 3, 5 and 7 along the z axis at different depths in the phantom. Monte Carlo uncertainties were less than 5% at y=4.3 cm, 6% at y=13.5 cm 2.6 Discussions and Findings In this chapter, we highlighted recent studies examining the discretization requirements of the deterministic SN method compared with the Monte Carlo method with application to dosimetry problems. In doing so, we chose to simulate an idealized 3-D 60C0 model using the PENTRAN code. Since radiation transporting in such applications encounters little interaction in low density materials, such as air, attention was given to the quadrature level and mesh size distribution needed for solution accuracy through the air and in a water phantom. As expected, radiation transport in this gamma ray application is highly directional, and that a traditional SN quadrature, such as an S12 leVel-Symmetric set (168 total directions), considered robust for many neutron reactor physics applications, is insufficient for modeling gamma rays in the 60C0 model. This limitation was overcome by the use of the higher levels of quadrature, up to S42 (1848 total directions per cell) with Legendre-Chebychev quadrature sets and P3 anisotropy. We also noted that at least an S32 angular quadrature leads to the elimination of most ray effects and scalar fluxes within the statistical uncertainty of Monte Carlo. Moreover, model mesh sizes were selected to preserve system volumes while taking into account the optical thickness of each material. Smaller Cartesian meshes do better represent the curvature of the model, but they require higher quadrature orders, and therefore larger memory and longer running times. Even in optically thin air regions, a mesh cell interval of on the order of 1 cm was necessary to yield an accurate solution, and more cells were necessary where the particle flux encountered strong gradients; cell intervals may be relaxed near boundaries where detail is not required. Also of importance was the flexibility of the ADS in PENTRAN, permitting more appropriate spatial differencing schemes. Scalar flux distributions obtained from PENTRAN agreed with the Monte Carlo prediction when simulating equivalent geometry models of the 60C0 model with three different multigroup cross-section libraries. Still, realistic representations of the energy spectra SN differencing schemes; these issues were investigated and presented in Chapter 2 in assessing what critical discretizations are required to yield accurate photon transport soluationS42 4.4.3 Euler Angles Euler angles, described in (Figure 4-8), provide a means of representing the spatial orientation of any frame of angular space as a composition of rotations from a given reference frame. The translations afforded by Euler angles were implemented to enable projections of the EDKs for a collection of cells aliased to a single reference pencil beam vector, as modeled in the Monte Carlo simulations (as in Figure 4-2) to generate EDKs to an arbitrary net current vector specifying the direction of the photons derived from SN computations. In the following discussion, we denote a fixed system in lower case (x,y,z) and a corresponding rotated system in upper case letters (X Y,Z). We use the term "Euler Angle" for any representation of three dimensional rotations of the electron kernel absorbed dose matrix that aligns with the direction of the net current for each DDV. y(2) Y casel sh op sn cs 0oc 0ia 1 csin 0 cosip 0os 0.o 1 Figure 4-8. Euler angles, zyz convention. Similar conventions can be obtained by following of the following order in rotation (zyz, xyx, xzx, yzy, yxy) Where : j: dN (t )dtd No 0 where dT -is the electron stopping power (MeV cm2/g p dxC To is the starting energy of the electron (MeV) dN (t ) tf (t) = is the differential distribution of farthest depth of penetration N (t ) is the number of particles penetrating a slab of thickness t Through inspection of these relationships described in Figures 1-1 and 1-2, one can decide what kind of tolerance is acceptable in applying a collisional kerma approximation for a dose calculation based on a particular voxel size, i.e. to determine CPE approximation limits. This also has been demonstrated in Table 1-3; for example at 2 MeV the CPE computed dose will be accurate to approximately 5% using a 1 cm voxel, assuming the CSDA range applies. Table 1-3. Approximate thickness of water needed to achieve various attenuation amounts compared to maximum energy electron ranges generated by the same beam P hotona Range (cm) at which CPE is Average" RewA Energy conserved by e- range (Mle 1 2%0 3% 5% (cm) (Cm') 0.3 0.085 0.17 0.256 0.431 0.056 0.085 0.6 0.112 0.225 0.34 0.572 0.15 0.225 1.0 0.142 0.286 0.43 0.726 0.286 0.43 2.0 0.203 0.41 0.617 1.04 0.617 1.04 c-~~' d3J\ S (+YH: Ir// / / ** LIST OF FIGURES FiMr page 1-1 Comparison between thicknesses required to attenuate different percentages of primary photons and Rcs range for maximum energy secondary electrons produced by the same beam ................. ...............33.____..... 1-2 Comparison between thicknesses required to attenuate different percentage of primary photons and average range for maximum energy secondary electrons produced by the same beam ................. ...............33.____..... 2-1 60C0 unit (left: Mesh distribution generated by PENTRAN package; Right: MCNP5 Cell/Surface schematic input) ................ ...............38........... .... 2-2 60C0 y-rays 3-D Groups 1,2 3 and 4 scalar flux computed by PENTRAN with the CEPXS gamma cross section library ................. ......... ...............40..... 2-3 S12 (Level Symmetric) and S42 (Legendre-Chebychev) quadrature first octants (with 21 and 231 directions/ octant respectively) .............. ...............41.... 2-4 60Co y-ray 3-D Group 1 scalar flux distribution for x-level depths at 0 cm, 5 cm, 10 cm and 15 cm, computed by PENTRAN with the CEPXS gamma cross section library ................. ...............42................. 2-5 60Co y-ray 3-D Group 1 scalar flux distribution for x-level depths at 0 cm, 5 cm, 10 cm and 15 cm, computed by PENTRAN with the CEPXS gamma cross section library; an S42 angular quadrature (1848 directions) with P3 Scattering anisotropy ..........43 2-6 Scalar flux distributions using S12, S22, S32, and S42 with P3 Scattering anisotropy in PENTRAN using the BUGLE-96 gamma library .............. ...............44.... 2-7 60C0 model spatial meshing, axial level............... ...............47. 2-8 Scalar flux comparison of variation in the interval of 3-D mesh cells containing air; fine meshes of 1 cm size (along each axis) in air proved to be best in providing high detail flux distribution............... ..............4 2-9 Scalar flux distributions generated by PENTRAN for a single energy group of each library using the MCNP5 solution as the reference case. ................. .................5 2-10 Relative deviation of energy group 1 discrete ordinates results with Monte Carlo results ............... ...............52........_ ...... 2-11 Absorbed dose rate comparison between SN method using kerma approximation and MC techniques using F6 tally, MC statistical uncertainty were less than 3% on average ................. ...............53........ ...... I ~X(cm) Figure 5-3. Monte Carlo derived total absorbed dose rate distributions along the X-axis computed using *F8 MCNP5 tally for three different phantoms. MCNP uncertainty (20) average (10.0%), 500 Minutes run Pulse height tally using photon-electron mode (SBLLBSS) 4.00E-04 3.5 3.50E-04 -*-MC-Dose a 3.00E-04 - -6- Material 2.5 E 2.50E-04 n 52.00E-04 Z 2~ 1.50E-04 1.5 5.00E-05 -g 0.00E+00 0 10 20 30 40 50 X~rcm) Figure 5-4. Monte Carlo derived total absorbed dose rate distributions along the X-axis computed using *F8 MCNP5 tally for the _SBLLB SS_ phantom. MCNP uncertainty (20) average (7.0%), 1000 Minutes run 500 Minutes (Pulse height tally using photon- electron mode) 4.00E-04 - 3.50E-04 - 3.00E-04 - 2.50E-04 - 2.00E-04 - ~.50Ea4- 1.50E-04 - 5.010E-05 0.00E+00 - -*-- SBLLBSSS -( SSBBSSS -8- SSSSSSS Dose rate using MCNP5 different tallies 3.00E-01 2.50E-01 2.00E-01 -*- iI Energy deposited tally using (p,e mode)" o -a-Pulse height tally using (p,e 1.50E-01 .-Y mode) S1.00E-01 g 5.00E-02 0.00E+00 g 0 5 10 15 20 25 Depth (cm) Figure 4-4. Comparison between F6 (p, e mode) and *F8 (p mode) tallies in MCNP5, statistical uncertainty was less than 1% I ODose rate (p~e mode Vs. p mode) 1 Figure 4-5. Comparison between *F8 (p, e mode) and *F8 (p mode) tallies in MCNP5, statistical uncertainty was less than 1% 0.30 0.25 0.20 0.15 0.10 0.05 0.00 -*Pulse hight tally with (p,e mode) -.Pulse hight tally with (p mode) ~--3Ccrl 101 Depth (cm) calculations20,21, although experience has shown that 3-D deterministic electron transport is quite challenging in order to yield a converged solution. To properly treat 3-D electron transport physics deterministically, yet still achieve reasonably fast and accurate whole body computation times for simulations with high energy photons, we have developed angular-energy dependent transport "electron dose kernels," or EDKs that can be coupled to SN solutions. These kernels were derived via full physics MCNP5 Monte Carlo electron transport simulations, and are applied based on deterministic photon solutions over the problem phase space, thereby accurately accounting for the absorbed dose that is ultimately rendered by charged particle electron transport. This is possible because the PENTRAN discrete ordinates code preserves detailed photon flux and current data on a voxel-by-voxel basis, so that accurate whole body absorbed doses may be rapidly achieved for high energy photon sources by performing a single deterministic SN multigroup photon calculation on a parallel cluster. This is followed by linking the SN-derived photon fluxes and net currents to our Monte-Carlo-based Electron Dose Kernels to account for a full-physics absorbed dose; we call this the 'EDK-SN' methodology for determining 'total body dosimetry'. This 'total body dosimetry' enables dose to be computed across the entire phantom phase space, wherever it is required, presented in the following sections. 4.2 Monte Carlo Based Dose Kernels Using Monte Carlo calculations, electron absorbed dose kernels for pre-determined photon energy groups were generated; energy was deposited in a collective "cloud" of voxels (i, j, k) composed of tissue, where absorbed doses attributed to these voxels are directly caused from incident primary photons in a given energy group that propagate electrons from voxel, (i',jk') the 'dose-driving voxel' or DDV. When producing these 'electron absorbed dose kernels,' the photons having energies lower than the limits of the energy group are cut off as they are R V~(rR, E)+ o(r, E )~(rR, E) dE dn' o- (r,0'- OE'4 E) YlI r, O',E') 0 4xr (1-6) 8[S (r, i',E)Y(r, 0',E ')] +q(r, 1, E) dE However, the full Boltzmann-FP equation is extremely challenging to solve deterministically . 1.4.1 Monte Carlo Method Monte Carlo techniques are a widely used class of computational algorithms for the purpose of statistically modeling various physical and mathematical phenomena. In general, individual probabilistic events are simulated sequentially, where probability distributions governing particle interactions are statistically sampled to describe the total phenomenon. In principle, for each case sampled, there are four required inputs; namely, a random number generator, source characteristics, a simulation phase space, and probability distributions governing the possible outcomes. Results of individual particle histories are generated using pseudo-random numbers to sample the probability distributions. This process is repeated many times, until the statistical precision required for the solution is achieved, based on the law of Large Numbers and the Central Limit Theorem"l. Therefore, the Monte Carlo technique is an effective method for indirectly solving the Boltzmann radiation transport equation via statistical methods, and can readily treat the full physics of coupled photon-electron transport. Simulations for particle histories are repeated many times for large number of histories, where the computer code tracks the desired physical quantities, means, and uncertainties (variances) for all of the simulated histories. Soft tissue (ICRU-44) 1000 +1CRU-44 -rFIT 100 1.00 E-03 1.00E-0n 1.00E-01 1.00 E+00 1.00E+01 1.0+2 0.1 Energy(Mev) Cortical Bone (ICRU-44) 1 00E+03 lCRU-44 'uI 1 00E+02 -I FIT 0 1 00E+01 S10 -03 1.00E-02 ~1.00E-01 1.00 E+00 1.00E+01 1.0+2 1 00E-01 Energy(Mev) DOSE calculates absorbed dose in each voxel yielding a whole-body absorbed dose distribution and dose-volume histograms for critical organs of interest. Energy(Mev) Figure 3-3. Values of the mass energy-absorption coefficient, PdenlP, and the fitting function as a function of photon energy, for Lung (ICRU-44), soft tissue (ICRU-44) and cortical bone (ICRU-44) 3.4 Computational Example For this work, a pediatric patient (an 11i-year-old male) was voxelized with high resolution using 0.94 mm x 0.94 mm x 6.00 mm voxel size (398 voxel x 242 voxel x 252 voxel) as a case 3-1 PENTRAN-MP Code System Structure Flowchart ................. .............................58 3-2 The UF Series B pediatric phantoms series ................ ...............60.............. 3-3 Values of the mass energy-absorption coefficient, IdenlP, and the fitting function as a function of photon energy, for Lung (ICRU-44), soft tissue (ICRU-44) and cortical bone ................ ...............64....... ...... 3-4 UF-Series B phantom (11 years old male) with (398 x 242 x 252 voxels); the corresponding PENTRAN "back-thinned" models with (132 x 80 x 84 voxels) and 72 materials (upper right), and another with (79 x 48 x 50 voxels) and 66 materials (low er ri ght) .............. ...............66.... 3-5 Transversal view of the 11i-year-old male Series B original phantomcollapsed into 189k voxels and 4 materials, collapsed into 900k voxels and 4 materials (lower left) and collapsed into 900k voxels and 72 materials............... ...............6 3-6 x-ray 3-D Group 1, 3, 5 and 7 scalar flux computed by PENTRAN with the CEPXS cross section library, S32 angular quadrature (1088 directions) with P3 Scattering ani sotropy ................. ...............68.._._._....... 3-7 Flux distributions using S32 with P3 scattering anisotropy in PENTRAN using the CEPXS library vs. MCNP5 Monte Carlo results ................. ...............68........... .. 3-8 Angular flux distributions using S32 with P3 scattering anisotropy in PENTRAN using the CEPX S library vs .............. ...............69.... 3-9 Angular flux distributions using S32 with P3 scattering anisotropy in PENTRAN using the CEPXS library vs. MCNP5 Monte Carlo results for energy groups 1, 3, 5 and 7 along the z axis at different depths in the phantom............... ...............69 3-10 Simulation methodology for 3-D absorbed dose computations using PENTRAN-MP code system ................. ...............71................. 3-11 CT view of the 11i-year-old male Series B original phantom collapsed into 189 k voxels and the correspondent absorbed dose distributions 900 k voxels and the correspondent absorbed dose distributions .............. ...............72.... 3-12 CT view of the 11i-year-old male Series B original phantom collapsed into 189 k voxels and the correspondent absorbed dose distributions 900 k voxels and the correspondent absorbed dose distributions .............. ...............73.... 3-13 CT view of the 11i-year-old male Series B original phantom collapsed into 189 k voxels and the correspondent absorbed dose distributions 900 k voxels and the correspondent absorbed dose distributions .............. ...............74.... 3-14 Simulation methodology for 3-D organ absorbed dose computations using PENTRAN-MP code system .............. ...............74.... 5.1.2 PENTRAN Model As shown in Figure 5-5, we simulated equivalent three phantoms to those simulated with Monte Carlo. The heterogeneous models consisted of three key components; an X-ray source (15 x1 x15 cm3) with flat weighted 8 MV source, divided into 16 energy groups, a soft tissue phantom (45x45x45 cm3), eqUIValent layers of bone and lung tissue identical to that used in the Monte Carlo simulations. Models were divided into 3 z-levels; each z-level was divided equally into 9 coarse meshes containing (15 x15 x15) fine mesh cells for each coarse mesh; a total of 91,125 fine mesh cells were used in the SN problems. Although lower quadrature set (S32) aef sufficient to produce accurate SN results for similar calculations (Chapter 1), to assure minimum ray-effects on the SN calculations, these simulations were based on the use of an S62 angUarT Legendre-Chebychev quadrature with P3 Scattering anisotropy using the CEPXS library. SSSSSSS SSSBSSS SBLLBSS Figure 5-5. Three SN phantoms equivalent to the one used in the Monte Carlo simulation. Phantom SSSSSSS is a homogeneous soft-tissue phantom; phantoms SSBBSSS and SBLLBSS are heterogeneous phantoms consisting of combinations of soft and bone tissue, and soft, bone, and lung tissue-equivalent materials, respectively requirements for deterministic calculations discussed in previous sections using different quadrature sets. Table 3-2. Total number of directions, Cumulative Problem Time required for different energy groups, P3 Scattering, and 189k fine mesh cells Number SN Number Parallel Cumulative of (relative of Decomposition Problem Time Energy workload) Processors (hours) groups A E S 4 42 (2) 12 4 1 3 4.1 6 24 (1) 12 8 1 2 2.1 8 32 (2.32) 12 4 1 3 3.69 8 32 (2.32) 16 8 1 2 3.06 Parallel decomposition results for solving the phantom are presented in Table 3-2 based on a constant problem phase space of 189,000 cells with P3 aniSotropy and different numbers of energy groups numbers and several angular quadratures. "A, E, S" refers to the angular, energy group, or spatial decomposition applied to solve the problem using parallel computing. From the table data, it can be observed that increasing the number of directions requires longer running times that essentially scale with processor workload. 3.6 Discussion and Findings We have demonstrated a methodology incorporating a discrete ordinates (SN) approach that is an effective alternative to the Monte Carlo method for accurately solving the transport equation over an entire human phantom. Models here were based on the UF series B pediatric phantoms developed by the University of Florida ALRADS group. We simulated the complete phantom discretized using 900k/186k voxels, and demonstrated that the SN method in the PENTRAN-MP code system, developed as a team effort, readily provides accurate whole body fluxes. Radiation transport throughout the phantom encounters few interactions in low density 1.00E+10 1.00E08 --*-5m Air mesh -*2cm Air mesh -*-1cm Air mesh 1.00E+06 -I ---MCNP5 L.. 1.00E+04 1.00E+02 1.00E+00 0.00 30.00 60.00 90.00 Distance (cmr) Figure 2-8. Scalar flux comparison of variation in the interval of 3-D mesh cells containing air; fine meshes of 1 cm size (along each axis) in air proved to be best in providing high detail flux distribution Using the ADS described earlier with different mesh sizes, we were able to achieve a converged solution. Figure 2-8 depicts the solution behavior versus mesh size in the air region (even though there is very little optical thickness or attenuation by the air); this can be explained by the fact that larger mesh sizes in the air do not adequately propagate the radiation from the source, and therefore result in solutions with significant truncation error, particularly near the source region, while using smaller mesh sizes in air can provide good solutions with fine detail over the problem domain. 2.3.3 Investigation of Cross-Section Libraries: After satisfying the requirements for properly representing angular quadratures and spatial discretization, the use of an appropriate cross section library is also important in solving problems in medical physics. The scalar fluxes resulting from the deterministic calculations using different cross-section libraries have been examined (Figure 2-9) and were compared with V1 A V1 B Figure 4-9. EDK-matrix after rotation to direction cosines A: (-0.5, 0.35, 0.8) and B: (0.25,- 0.40,-0.88) relative to the reference coordinate system (0, 0, 1). Dose Driving Voxel is located at (0,0,0) which appear as unphysical flux oscillations in the scalar flux distribution unless treated properly. In general, this behavior occurs for problems with highly angular dependent sources and/or fluxes, or when the source/sink is localized in a small region of space, in low density or highly absorbent media. Grou 1 rou 2 rup3Gru Fiur 2-2 600yry rus12,3ad4saa xcmue yPNRNwt h CEPXS gamm crs eto irr,(S2aglrqartr (88drcin)wt P3i Scteigaiorp hr sd Tomtgt this rbehg nua udaueoresaeotnnee.Fgr - deit a sceai ofteodnt est itiuin soitdwt 1 0e-ymti (rotatonall lee-ymtic) and S4 LeenreChbyhe quaraur sets on 1/to ui thgue origi ando p-rojec toD e rus1 nd at thea squar "pont"plo ted inteFiue BeyndRA S20, Legede Chbyhe quadratue sets polmus beg nused due t neative wreighs tate occur ine standrd lee- symmetric sets36. Particular attention must be made in deriving quadrature sets to enable preservation of Spherical harmonics (Legendre) scattering moments and numerical quadrature orthogonality relationships. O 2008 Ahmad Al-Basheer Table 5-3. Absorbed dose rate (MeV/g. s) calculated for the test problem using MCNP5 using pulse height tally *F8 with (photon, electron mode) and MCNP5 energy deposited tallies (F6) Organ MC(*F8) 2o--MC MC: (F6) (*F8-F6t)/F6t (p, e mode) Uncertainty Right Lung 5.62E-06 2.0% 4.28E-06 31.3% Left Lung 2.70E-04 0.40% 2.74E-04 1.54%~ Liver 1.00E-05 3.0% 9.14E-06 9.41% Pancreas 1.11E-06 4.0% 1.03E-06 7.84% Wang et al. conducted a similar study using their VIP Man organ MCNP for absorbed dose comparison for the three different tally methods including cell flux tally (F4), energy deposited tally (F6) and pulse height tally (*F8). In their comparison they found that the F4 tallies agree within 10% with the F6 tallies, while the differences between the tallies F8 and F6 are relatively large (28%) 14. These results agree with our findings, and provide a strong proof for the need of applying full physics transport when calculating whole body absorbed dose. 5.3 Discussion and Findings To properly treat 3-D electron transport physics deterministically, yet still to achieve reasonably fast and accurate whole body computation times using high energy photons, we used angular and energy dependent electron transport "absorbed dose kernels." In a methodology called 'EDK-SN' we have demonstrated that by using different electron kernels for different materials (cortical bone, soft tissue and lung), the EDK-SN methodology incorporating a discrete ordinates (SN) approach is an effective alternative to the Monte Carlo method for accurately estimating the absorbed dose for high energy photon beams over an entire human phantom. In doing so, we simulated three phantoms using MCNP5 and PENTRAN-MP. Each phantom had a flat weighted 8 MV x-ray source impinging on 45x45x45 cc soft tissue phantom. Two of these phantoms were slabs of 1-cm-thick discs representative of bone, or lung tissue. Comparisons of MCNP5 Monte Carlo based dose results and those computed using PENTRAN-MP implementing the EDK-SN methodology showed good agreement between the two approaches. 60 N 40 20 ~ "1 ,o'~-~ ~;51~~ Y~cm) Ylcm) 0 cm 5cm 60 A 40 20 0 60 N 40 20 0 ~~ " Y (cm) Y (cm) 10 cm 15cm Figure 2-4. 60Co y-ray 3-D Group 1 scalar flux distribution for x-level depths at 0 cm, 5 cm, 10 cm and 15 cm, computed by PENTRAN with the CEPXS gamma cross section library; an S12 angular quadrature (168 directions) with P3 Scattering anisotropy was used as a reference case S12 S42 Figure 2-3. S12 (Level Symmetric) and S42 (Legendre-Chebychev) quadrature first octants (with 21 and 231 directions/ octant respectively) Using a uniform 1 cm mesh (for reasons discussed later) in the air region, PENTRAN was used to calculate the gamma flux distributions in the system for several different quadrature sets. PENTRAN provides the option of readily using any quadrature that could be prescribed by the user. The flux distributions produced by an initial computation, using S12 leVel-Symmetric quadrature with P3 Scattering anisotropy and 126,000 fine mesh cells~ is shown in Figure 2-4. The ray-effects in the S12 COmputation are evident when compared to scalar fluxes in Figure 2-5 using S42 quadrature, with a similar mesh structure and P3 Scattering (the reference case). While there is evidence of small ray effects in the reference case far from the source, the differences between Figures 2-4 and 2-5 is dramatic. We also found that there was a small difference (on the order of several percent, depending upon location) in flux solutions computed with S42 between P3 and Pi computations, but little difference noted in the problems tested between P3 and P5 anisotropy, and so P3 Was implemented for all other computations. To further investigate the effects of raising the quadrature order, unbiased quadrature sets were tested using S22, S32, etc, up to S42 with P3 aniSotropy. A comparison of the PENTRAN computed fluxes in Group 1 using selected quadratures and equivalent MCNP5 calculations are given in Figure 2-6. stylized phantoms, while maintaining the anatomic realism of segmented voxel phantoms with respect to organ shape, depth, and inter-organ positioning to great accuracy. One study showed that these phantoms match (to within 1%) to that of a newborn of the ICRP Publication 89 reference newborn child. In this study, hybrid computational phantoms representing reference 15-year male body anatomy and anthropometry is used, shown in Figure 5-10. This phantom was matched, by the ALRADS group, to anthropometric data and reference organ mass data provided by the ICRP, where it matched to within 1% overall" Figure 5-10. Frontal views of 3-D rendering of the 50th percentile UFH-NURBS 15-year male and female The availability of hybrid phantoms as a replacement for voxel-based phantoms provides the needed platform for constructing phantoms representing the unique internal and external -1.03 -1.03 -1.03 -0.41 -1.09 -1.04 -1.06 -1.03 -1.02 -1.04 -1.03 -1.05 -1.05 -1.02 -1.03 -1.04 -1.00 -0.001 -1.03 -0.34 -1.03 -1.03 -1.10 -1.10 -1.10 -1.10 -1.519 -1.519 -1.387 -1.351 -1.219 -1.351 -1.196 -1.196 -1.196 -1.219 -1.351 -1.319 -1.319 -1.387 -1.319 -1.352 -1.352 -1.352 -1.352 -1.273 -1.273 -1.273 -1.273 -1.35 -1.415 -1.319 -1.319 -1.378 -1.378 -1.374 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 imp:p,e=1 u=39 u=40 u=41 u=42 u=43 u=44 u=45 u=46 u=47 u=48 u=49 u=50 u=51 u=52 u=53 u=54 u=55 u=56 u=57 u=58 u=59 u=60 u=61 u=62 u=63 u=64 u=65 u=66 u=67 u=68 u=69 u=70 u=71 u=72 u=73 u=74 u=75 u=76 u=77 u=78 u=79 u=80 u=81 u=82 u=83 u=84 u=85 u=86 u=87 u=88 u=89 u=90 u=91 u=92 u=93 u=94 vol=3.40E+01 $Salivary Glands parotidd) vol=1.00E+01 $Scrotum vol=1.40E+02 $SI W vol=7.70E+02 $SI C vol=1.72E+02 $Skin vol=4.10E+01 $Spinal Cord vol=1.21E+02 $Spleen vol=3.30E+01 $Stomach W vol=2.08E+02 $Stomach C vol=1.80E+01 $Testes vol=3.20E+01 $Thymus vol=8.00E+00 $Thyroid vol=5.50E+01 $Tongue vol=2.00E+00 $Tonsil vol=6.00E+00 $Trachea vol=3.00E+00 $Urinary bladder W vol=1.61E+02 $Urinary bladder C vol=4.23E+02 $Air (in body) vol=7.00E+01 $Left Colon W vol=1.95E+02 $Left Colon C vol=2.00E+01 $Salivary Glands submaxillaryy) vol=7.00E+00 $Salivary Glands (sublingual) vol=2.20E+01 $Coastal cartilage of ribs vol=0.00E+00 $Cervical Discs vol=2.80E+01 $Thoracic Discs vol=1.20E+01 $Lumbar Discs vol=7.48E+02 $Cranium vol=3.30E+01 $Mandible vol=1.06E+02 $Scapulae vol=3.00E+01 $Clavicles vol=4.40E+01 $Sternum vol=2.90E+01 $Ribs vol=5.10E+01 $Vertebrae-C vol=2.28E+02 $Vertebrae-T vol=1.95E+02 $Vertebrae-L vol=1.29E+02 $Sacrum vol=4.67E+02 $Os Coxae vol=2.20E+02 $Femur-proximal vol=1.63E+02 $Femur-upper shaft vol=1.33E+02 $Femur-lower shaft vol=2.59E+02 $Femur-distal vol=2.07E+02 $Tibiae-proximal vol=1.01E+02 $Tibiae-upper shaft vol=1.01E+02 $Tibiae-lower shaft vol=7.30E+01 $Tibiae-distal vol=1.40E+01 $Fibulae-proximal vol=4.00E+00 $Fibulae-upper shaft vol=1.40E+01 $Fibulae-lower shaft vol=1.70E+01 $Fibulae-distal vol=2.40E+01 $Patellae vol=4.30E+02 $Ankle+Feet vol=1.50E+02 $Humerus-proximal vol=7.40E+01 $Humerus-upper shaft vol=6.70E+01 $Humerus-lower shaft vol=9.30E+01 $Humerus-distal vol=1.60E+01 $Radii-proximal BUGLE-96 PE NTRAN vs ENDFIB-VI MCNP5 1.00E+10- -o- MCNP5 1DEO1.00E+06 -PNR4 1.00E+04- 1.00E+02- 1.00E 100 0.0 30.0 60.0 Distance (cm) A C EPX(S PENIT RAN i ENDFW-VIIMONP5 SCALE PENITRAN vs ENDF B-UI MONIP5 1.00E+10 1.00E+10 1.00E+08 cr 1.00E+08- PEIITRAl 1.00E+06 \I 1.00E+06 1.00E+04 -II1.00E+04- . 1.00E+02 -I 1.00E+02 1.000+ 00 1 00E+00 000 30.00 60 0 000 3000[ 60100 Distance~cm] Distance~cm) B C Figures 2-9. Scalar flux distributions generated by PENTRAN for a single energy group of each library using the MCNP5 solution as the reference case. These fluxes span the central axis in the three layers of the model. In general, the deterministic results yield a very good match to the Monte Carlo results CHAPTER 6 SUMMARY, CONCLUSIONS AND FUTURE WORK This work began with the primary goal of developing a new parallel 3-D Dose calculation tool based on deterministically rendered high energy photon transport solutions in radiotherapy settings. 6.1 Conclusions and Findings In this dissertation, the EDK-SN method is presented for generating detailed absorbed dose information using high energy external beam source terms. Because detailed whole body and organ-specific absorbed doses can be accurately obtained throughout a voxelized phantom with proper electron interaction physics projected onto the phantom voxels, based on SN-computed photon fluxes and net current vectors, this highlights the strength of the EDK-SN transport based methodology. Overall, based on the problems examined, with an adequate discretization implemented using the EDK-SN method, it is possible to produce absorbed dose results consistent within the statistical uncertainty of MCNP5 Monte Carlo calculations. Overall, our studies examined the capabilities of the discrete ordinates method compared with the Monte Carlo method with application to medical physics problems. Since radiation transporting in such applications encounters few interactions in low density material, such as air, attention was given to minimum quadrature levels and mesh size distributions necessary for solution accuracy in air and in a water phantom. Also, it is evident that multi-group cross-section libraries must also be carefully tested for problem applicability. Of particular importance was the flexibility of the adaptive SN differencing scheme in PENTRAN, permitting more accurate spatial differencing schemes in cells with a range of optical thicknesses. In addition, deterministic models based on the UF series B pediatric phantoms and with the PENTRAN-MP 3 WHOLE BODY DOSIMETRY SIMULATIONS WITH THE UF-B SERIES PHANTOM USING THE PENTRAN-MP CODE SYSTEM FOR LOW ENERGY PHOTONS ............ ..... ._ ...............57.... 3.1 Methodology .................. .. ...............58. 3.2 UF Series B Pediatric Phantoms .............. ...............59.... 3.3 PENTRAN-MP Code System ............ ..... ._ ...............61... 3.3.1 Pre-processing Codes.................. ..........___......__ ............6 3.3.1.1 The collapsing code (GHOST-3-D)............... ..............6 3.3.1.2 The mesh generator code (PENMSH-XP) ................. ............... ....61 3.3.1.3 S ource spectra and cros s- secti ons generate on ................ ..................62 3.3.2 Transport Calculation............... ..............6 3.3.3 Post-Processing Codes ................. ..................... ........6 3.3.3.1 Data management utility (PENDATA) .............. .....................6 3.3.3.2 Dose calculation code (3-D-DOSE) .............. ....................6 3.4 Computational Example............... ...............64 3.5 Computational Results .............. .. .. ........... ......... .. .......6 3.5.1 Flux Comparison of Deterministic and Monte Carlo Results ................... ........67 3.5.2 Dose Calculations and Dose Volume Histograms (DVH). ............... ... ............70 3.5.3 Parallel Computation............... ..............7 3.6 Discussion and Findings .............. ...............76.... 4 ELECTRON ABSORBED DOSE KERNELS TO ACCOUNT FOR SECONDARY PARTICLE TRANSPORT IN DETERMINISTIC SIMULATIONS .............. ...............79 4.1 Introduction to The EDK-SN Method .............. ...............79.... 4.2 Monte Carlo Based Dose Kernels................. ......... ....................8 4.3 EDK Determinations Based on The MCNP5 Monte Carlo Code ............... ... ........._...82 4.3.1 Computational Model............... ...............83. 4.3.2 MCNP5 Tallies Comparison ................. ...............83................ 4.4 EDK -SN Code System Methodology .............. ...............86.... 4.4.1 N et Current............... ..... ............8 4.4.2 MCNP Based Electron Kernels............... ...............87 4.4.3 Euler Angles............... ...............89. 4.4.4 Ordinate Proj section ................. ...............92.____ .... 4.4.5 Final Dose Calculations .............. ...............93.... 4.5 External Source Benchmark............... ........ ...............9 4.6 Comparison of Deterministic and Monte Carlo Results ..............._ ............_.........96 4.7 Parallelization of EDK-SN Code ................. ......... ...............99..... 4.8 Discussion and Findings .............. ...............101.... 5 APPLICATION OF EDK-SN TO MODELS WITH HETEROGENITIES .........................103 5.1 Heterogeneous Tissue-Equivalent Phantoms ......____ ..... ... .__ ........_......103 5.1.1 Monte Carlo Simulations .............. ...............103.... 5.1.2 PENTRAN Model ................. .......... ...............107.... 5.1.3 EDK-SN Absorbed dose Comparison............... ..............10 On the other hand, Monte Carlo calculations for organs outside the x-ray field (such as the prostate) showed significant statistical uncertainty (44%). This poor absorbed dose tally result is due to the fact that the Monte Carlo did not have enough histories for these results to converge to low statistical uncertainty in an acceptable time. Alternatively, as shown in Table 5-2,when running the Monte Carlo for extended running times (~40 hours on our parallel cluster), the Monte Carlo calculated absorbed dose delivered to the prostate began to show convergence to the absorbed dose calculated using EDK-SN. This illustrates that the EDK-SN methodology, executed in 2 hours (SN execution time = 1.5 hr, EDK- SN execution time = 0.5 hr), is capable of providing an accurate solution and a significant speed up for a full physics absorbed dose calculation in human phantoms for high energy external beams. Table 5-2. Prostate absorbed dose rate (MeV/g.s) calculated for the test problem using EDK-SN and MCNP5 using pulse height tally (*F8) with (photon, electron mode) at different MCNP5 Monte Carlo running times Dose 2o--MC (EDK-SN- Running (M~ev/gm.Sec) Uncertainty MC)/EDK-SN Time ED)K-SN 2.29E-08 ----- ----- 2hr MC1t 2.21E-08 44.00% 3.62% 3.9hr MC2 2.34E-08 39.20% 2.14% 4.7hr MC3 2.18E-08 34.40% 5.05% 6;.1hr MC4 2.23E-08 22.20% 2.69%/ 14.5hr MC5 2.46E-08 16.86;% 6.91% 28.15hr MC6 2.406E-08 14.14% 4.58% 38.85hr 5.2.2 Ray-Effects When considering a forward peaked volumetric (10 cm x 1 cm x 17 cm) X-ray source placed on top of the left lung area of the phantom, we noted significant differences between the absorbed dose calculated by EDK-SN and the absorbed dose calculated by full physics Monte (4-4.b) (4-4.c) J = J (4-5) Where ~ju is the photon net current unit vector. Figure 4-6. Deterministic calculation will provide the angular information explicitly 4.4.2 MCNP Based Electron Kernels When generating these electron kernels, a single beamlet of mono-energetic photons, for a certain energy group, interacts at the center of a water phantom, as discussed in the previous section. For the purpose of this study, we generated the EDKs for 16 photon energy groups covering the energy spectrum from 0 to 8 MeV. Based on our previous study, as discussed in chapter 1, these energies will violate CPE for the model considered. Moreover, the width of the energy group is set to be 0.5 MeV, minimizing the effect of energy discretization in an acceptable SN running time (we note that energy group discretizations can be optimized for better Jnet I~-~et scattering medium such as air, or in a strongly absorbing medium such lead, where the photon flux usually is typically strongly angular dependent. In such situations, either a higher order quadrature set (many directions) is required or some ray-effects remedies must be appliednl~l (such as Taylor Projection techniques and adaptive differencing schemeS19). More discussion on discretization and ray-effects are presented in Chapter 2. 1.4.2.2 Monte carlo versus deterministic SN methods The most fundamental difference between the Monte Carlo methods and the deterministic SN transport method is that the SN method solves the transport equation for the average particle behavior throughout the problem phase space, while the Monte Carlo method, the terms of the transport equation are solved statistically by simulating the individual particles and recording some aspects (tallies) of their average behavior, but only at specific tally locations requested by the user. In the discrete ordinates method, the phase space is divided into many small meshes, in which the cross section is assumed to be constant. In the medical physics applications, the model geometry, in this case, the patient is represented using regular voxel-based grid generated from a segmented CT scan data or other digital imaging technique. The types of grids can then be easily represented in an input deck for the SN calculations, making both Monte Carlo methods and deterministic techniques capable of representing the model of study accurately. Another advantage of the MC methods is that they can easily utilize continuous energy interaction data, whereas deterministic calculations use multigroup or energy group averaged data. However, photon attenuation coefficients in the energy range of interest for radiotherapy (RT) are quite smooth with little structure (in contrast to neutron cross sections); and are accurately represented by multigroup data. Recent research we conducted has shown that multi-group deterministic methods produce highly accurate absorbed dose distributions in photon the source, were 4% on average while the average absolute difference at the same cells was 3.7%. ED]K Va. Pulse height tally using photon-electron 4 00E04 -mode (SSSSSSS) S3.50 E-04 n 3 00 E-04- E -e EDK-SN Ea 2.50 E-04 ~5;-*-MC-Dose 2 00 E-04- QL 1.50 E-04 2u 1.00 E-04- a 5 00 E-05 0.00 E+00- ,, O 10 20 30 40 50 X(cm ) Figure 5-7. EDK-SN Homogenous phantom total absorbed dose rate distributions along the X- axis compared to MCNP *F8 tally. MCNP uncertainty (20) average (6%) on-axis. Average absolute relative difference (3.7%) We note that at the interface lung region, Figure 5-6, the EDK-SN derived absorbed dose indicated a lower absorbed dose than Monte Carlo. This can be related to the application of the lung electron kernel over the bone and soft tissue regions. This phenomenon was mitigated by applying a correction factor for voxels located at the interface of different material densities. 5.1.4 Correction Factor Approach The range of a charged particle of a given type and energy in a given medium is the expectation value of the path-length that it follows until it comes to rest. On the other hand, the RCSDA TepreSents the range in the continuous slowing down approximation (CSDA). Figure 5-8 gives the CSDA range RCSDA for electron in cortical bone, soft tissue and lung tissue for electron energies between 0.0 MeV to 17.5 MeV. Figure 5-9 depicts a comparison between the ratio of application of integrated electron transport kernels coupled to the SN photon transport methodology to properly account for the total absorbed dose in high energy photon irradiation scenarios. A fundamental difference between the deterministically guided EDK- SN technique and the existing convolution-superposition technique is the fact that our EDK-SN approach is based on a full physics Monte Carlo generated kernel projected using an accurate SN transport solution that provides a net current direction and flux scaling to scale the absorbed dose using the spectrally dependent EDKs. Overall, based on the water phantom problem examined here, with an adequate discretization using the SN method coupled to properly computed energy dependent electron kernels, it is possible to produce absorbed dose results with the EDK-SN method that are consistent with Monte Carlo calculations. A major advantage of the EDK-SN method is that it provides a detailed, accurate absorbed dose distribution wherever required throughout the system, while incorporating the full physics electron transport, and is generally much faster than Monte Carlo. Since PENTRAN is based on a Cartesian geometry, CT voxel data can readily be mapped onto the problem geometry to define the tissue and organ geometry structure. Subsequent discussion will account for material density variations (encountered at tissue-bone interfaces) where we will apply the EDK-SN methodology to a whole body voxelized phantom to rapidly achieve a full physics whole body, 'total body' absorbed dose. 4.3.1 Computational Model To highlight the differences between the different tallies, used for absorbed dose estimation, in the MCNP5 code, we considered a box of 11x1Ix1 1 cm3 COmposed of a water cube containing a lxlxl cm3 8 MeV forward peaked X-ray source placed at the center of the cube Figure (4-2). The size of this model allows complete deposition of secondary electrons produced from the X- ray source. A flux mesh-tally, using "DE" absorbed dose energy and "DF" absorbed dose function cards, was generated to cover the entire model. Moreover, since no mesh tally is available for *F8 and F6 tallies in the MCNP5 code, a 3-D tally distribution of cells, equivalent in purpose to the MCNP5 F4 mesh tally, were implemented in the model. This model, with the implemented tally sites, was used for highlighting the differences between the different tallies, and consequently was used for reporting a scaled full physics absorbed dose. 4.3.2 MCNP5 Tallies Comparison The source in the center of the model was applied to generate a pencil beam impinging perpendicularly, along the z-axis, on the box of water. As shown in Figures 4-3, 4-4, and 4-5, the absorbed dose is highly forward peaked, with little dose resulting from back-scattered radiation. Absorbed doses resulting from the use of an F4 volume flux tally method, using a "DE" dose energy and "DF" dose function cards, compared to an F6 tally using both the photon-electron mode (p, e mode) and photon only (p mode) were all in agreement, since the two tallies depend on similar physics in predicting absorbed dose. These results are depicted in Figure 4-3. On the other hand, when comparing absorbed doses computed using a *F8 using a fully coupled photon electron transport mode (p, e mode) and the F6 derived dose tallies, the differences were significant, as depicted in Figure 4-4. This is expected, since F6 tallies are analogous to F4 tallies in MCNP, and can be directly aliased to assumptions equivalent to CPE, while *F8 tallies are determined based on full charged particle transport physics, as discussed below31 shows a comparison between Eine meshes side size of 1 cm, 2 cm and 5 cm in the air region with the DTW scheme fixed in air. These results demonstrate that a 1-cm mesh in air is necessary to maintain accuracy. Based on these observations, we designed our PENTRAN model as follows; the model is 30x30x70 cm3, and consists of 45 total coarse mesh blocks arranged among 5 axial levels, containing a total of 126,000 Eine meshes. In recent PENTRAN versions, selection of the Exponential Directional Iterative (EDI) method is made as an improvement over EDW39. The geometry of the model has been represented using the block adaptive meshing capability of the PENTRAN code that allows the user to use a higher mesh density in the regions of interest. Axial level thicknesses were determined based on major components of our model, and the X-Y spatial dimensions of coarse meshes (containing various densities of Eine meshes) were uniformly 10 cm on each side. The DTW differencing scheme was locked for the regions with low density material (air regions), and the adaptive differencing scheme (with adaptive upgrade potentially up to the EDW differencing scheme) was prescribed for higher material density regions. A mesh coupling scheme, Taylor Proj section Mesh Coupling (TPMC)40, was used to link discontinuous meshes using first order projections as radiation sweeps through the problem. The Cartesian mesh sizes for the full size model were selected to preserve component volumes, also taking into account the optical thickness of each material. The Eine-mesh intervals in different problem regions are: Source- 0.5 cm; lead- 0.5 cm; water- 1.0 cm; Air- 1.0 cm. Smaller mesh sizes better represent curved boundaries and irregular structure in the model, but require larger memory and may lead to longer running time (depending upon the spectral radius of the equations). Smaller mesh sizes were selected for the source region of the problem, where the particle flux has a strong gradient. Mesh sizes in the air can be made larger near the edges of the problem where a more detailed solution is not required. 5-1 Three Monte Carlo phantoms used in the ab sorbed dose study ................. ................. .1 04 5-2 Monte Carlo derived total absorbed dose rate distributions along the X-axis computed using *F8 MCNP5 tally for three different phantoms .............. ..................105 5-3 Monte Carlo derived total absorbed dose rate distributions along the X-axis computed using *F8 MCNP5 tally for three different phantoms .............. ..................106 5-4 Monte Carlo derived total absorbed dose rate distributions along the X-axis computed using *F8 MCNP5 tally for the _SBLLBSS_ phantom ................. ...............106 5-5 Three SN phantoms equivalent to the one used in the Monte Carlo simulation. .............107 5-6 EDK-SN heterogeneous SBLLB SS phantom total absorbed dose rate distributions along the X-axis compared to MCNP *F8 tally............... ...............108. 5-7 EDK-SN Homogenous phantom total absorbed dose rate distributions along the X- axis compared to MCNP *F8 tally............... ...............109. 5-8 Continuous slowing down approximation range (Rcsda) for electrons in Cortical Bone, Soft Tissue and Lung as a function of energy ................. ......... ................110O 5-9 Continuous slowing down approximation range (RCSDA) ratio comparison for electrons in Cortical Bone, Soft Tissue and Lung as a function of energy ................... ...11 1 5-10 Frontal views of 3-D rendering of the 50th percentile UFH-NURBS 15-year male and female ................. ...............112................ 5-11 Simulation methodology for EDK-SN computations using PENTRAN-MP code system ........... ..... ._ ...............113... 5-12 Absorbed dose contribution of each energy group is calculated via the EDK-SN procedure, then folded to calculate the total absorbed dose delivered to the 15 year old male phantom from 8 MeV flat weighted source ...........__.......__ ..............115 A-1 Total mass attenuation coefficient for photons in water, indicating the contributions associated with phoelectric absorption, Compton scattering and electron-positron pair production ................. ...............124................ A-2 Sketch of Compton scattering, photon of energy hv collide with an electron. The photon' s incident forward momentum is hvl/c, where c is the speed of light in vacuum ................. ...............125._._._....... A-3 Compton effects differential cross section per unite angle for the number of electrons scattered in the direction 8e............... ...............128.. A-4 Pair Production in the Coulomb force field of the nucleus. ..........._... ........_.._.........130 Resda Co mparison (Cortical Bone, Soft Tissue, Lung) 1.20E+00 1.00E+00 8.00E-01 -*-- Tissue/Bone 6.00E-01 ---Tissue/Lung u --6- Bone lung 4.00E-01 2.00E-01 0 .00E+ 00 0 5 10 15 20 Electron Energy(MeV) Figure 5-9. Continuous slowing down approximation range (RCSDA) ratio comparison for electrons in Cortical Bone, Soft Tissue and Lung as a function of energy The key to this approach is that photon transport is carried with the 3-D- SN code PENTRAN. The use of such code provided us with the directional component of the photon transport in the homogenous medium. By having this information we were able to project the 3- D electron kernel along the net current direction. The assumption of a homogenous medium in photon transport in a heterogeneous problem, often used in other methods, may lead to an underestimation of the scattered component of the photon beam, and misrepresentation of its directional component, especially between the different organs of the human body (Figure 5-9), leading to misrepresentation of the absorbed dose. 5.2 Computational Phantom of 15-Year Old Male A new set of stylized phantoms have been developed by the UF Advanced Laboratory for Radiation Dosimetry Studies (ALRADS). These phantoms represent a totally new class of whole-body patient hybrid phantoms so named as they retain the non-uniform scalability of TABLE OF CONTENTS page ACKNOWLEDGMENTS .............. ...............4..... LIST OF TABLES ........._.._ ..... ._._ ...............8.... LIST OF FIGURES .............. ...............9..... AB S TRAC T ............._. .......... ..............._ 14... CHAPTER 1 INTRODUCTION AND LITERATURE REVIEW ................. .............. ......... .....16 1.1 Overview ............... .... .. ......... ...............16...... 1.2 Initial Clinical Dose Evaluations .............. ...............17.... 1.3 Pencil Beam Transport ................. ...............19................ 1.4 Solving the Boltzmann Equation .............. ...............21.... 1.4.1 Monte Carlo Method ................. ...............23.......... ..... 1.4.2 Deterministic M ethods ................... ...... ........... ...............25...... 1.4.2.1 Discrete ordinates method (SN method) ................. ......................25 1.4.2.2 Monte carlo versus deterministic SN methods .................. ...............27 1.5 Dose and Charged Particle Equilibrium............... ......................2 1.5.1 Deterministic Dose Calculations for Diagnostic Application ................... ........30 1.5.2 High Energy Beam Dose Calculations ................. ...............31........... 1.5.3 Limits of Charged Particle Equilibrium ................. ............... ......... ...31 1.6 Outline of This Research............... ...............35 2 CRITICAL DISCRETIZATION ISSUES IN 3-D DISCRETE ORDINATES MEDICAL PHYSIC S SIMULATION ................. ...............36................ 2.1 Code Systems ................. ...............37........... .... 2.2 Simulation Geometry .............. ...............38.... 2.2.1 PENTRAN Model ................. ...............3.. 8......... .... 2.2.2 M onte Carlo M odel .............. .... ...... ....... .. ........3 2.3 Comparison of Deterministic and Monte Carlo Results ................ ............ .........39 2.3.1 Investigation of Angular Quadrature .........._.._.._ .....__ ......._.. .........39 2.3.2 Investigation of Spatial Discretization (Meshing) .............. .....................4 2.3.3 Investigation of Cross-Section Libraries:.. ...49............ 2.4 Absorbed Dose Calculations ................. ...............52......... ..... 2.5 Parallel Decomposition ................. ...............54............... 2.6 Discussions and Findings ........_................. ............_........5 referred to 'coarse z-levels', and then each z-level is partitioned into single coarse levels of 3-D x-y-z fine mesh grids. Predefined shapes can also be used to generate the desired shapes in a particular model. A second methodology is based on two binary input files, one for the phantom geometry, which can be established from real CT data, and the other containing the source density for each voxel/fine mesh, which can be used to specify a fixed-source distribution in the model. By combining the two files, PENMSH-XP directly generates an input deck for the PENTRAN transport code. 3.3.1.3 Source spectra and cross-sections generation Based on the bremsstrahlung and characteristic parametric functions of the model proposed in 1991 by Tucker, Barnes, and Chakraborty 46 the DXS code, principally developed by Monica Ghita, generates x-ray spectra in the diagnostic range (50-140 keV), according to the input parameters (tube potential, anode angle, type and amount of filtration), and accommodates them to any energy group structure. Optionally, as part of the PENTRAN-MP package, DXS also generates the input files for CEPXS (Sandia National Lab), the x-ray cross-section generator code, corresponding to a desired energy group structure and PN order. Group-dependent cross- sections rendered by CEPXS are then stripped from the output file and adapted to a format suitable for the PENTRAN deterministic calculations using another code developed by Ghita, GREPXS47 3.3.2 Transport Calculation Transport calculations were performed using the PENTRAN SN and MCNP5 Monte Carlo as discussed in Chapter 2. PENTRAN SN MCNP Electron High Re~solution Calculations Dose Kernel (ED=K) i n Di re ctio n One Direction for Projection Set Electro n Dose Generation Code! Using Euler Angles Correspondent Delta Module to Define 1 ~Discreet Angles Parallel EDK- SN Code for 3D- Dose Calculations EDK Sets for Each EG and Each Material 3D- Dose Organ Di stribtutio~n Do3se -PrCoue Model Specific Figure 4-11. A flow-chart for the EDK-SN Code System methodology. The function of this module is to define a AO on which the net current vector will be discretized into. For each specific model, after executing the deterministic calculation using PENTRAN, we generate 3-D fluxes, current vectors, and material distributions for each energy group for a particular model. These 3-D data distributions, along with pre-computed EDKs, will then be LIST OF REFERENCES 1R. Mohan, in XIlth International Conference on the Use of Computers in RadiRRR~~~~~RRRRR~~~~atio Therapy, edited by (Medical Physics Publishing, Madison, WI, 1995, Salt Lack City, 1997), pp. 16-18.35 2T. R. Mackle, P. Reckwerdt, and N. Papanikolaou, "3-D Radiation Treatment Planning and Conformal Therapy," M~edical Physics Publishing Corporation Madison WI, 1995), pp. 34 3B. A. Fraass, J. Smathers, and J. Deye, Summary and recommendations of a National Cancer Institute workshop on issues limiting the clinical use of Monte Carlo dose calculation algorithms for megavoltage external beam radiation therapy," Med Phys 30, 3206-3216 (2003).29 4F. M. Khan, The Physics ofRRRRRRRRRRRRRRRRRRRadito Therapy, 3 (Lippincott Williams & Wilkins, 2003).36 5P. R. Almond, A. E. Wright, and M. L. Boone, "High-energy electron dose perturbations in regions of tissue heterogeneity. II. Physical models of tissue heterogeneities," Radiology 88, 1146-1153 (1967).37 6A. Ahnesjo, M. Saxner, and A. Trepp, "A pencil beam model for photon dose calculation," Med Phys 19, 263-273 (1992).73 7A. L. EDWARDS, J. A. RATHKOPF, and R. K. SMIDT, "Extending the Alias Monte Carlo sampling method to general distributions,," in the International Conference on] Aubeinal'ltliL and Computations and Reactor Physics, edited by (American Nuclear Society, Inc., La Grange Park, IL, 1991), pp. 38 sA. Ahnesj o, "Collapsed cone convolution of radiant energy for photon dose calculation in heterogeneous media," Med Phys 16, 577-592 (1989).52 9T. Knoos, A. Ahnesj o, P. Nilsson, and L. Weber, "Limitations of a pencil beam approach to photon dose calculations in lung tissue," Phys Med Biol 40, 1411-1420 (1995).50 10B. Dionne, "Automated variance reduction technique for 3-D Monte Carlo coupled electron- photon-positron simulation using deterministic importance functions," Nuclear Engineering Sciences (University of Florida, Gainesville, Fl, 2007), pp. 75 11E. E. Lewis, and W. F. Miller, Computational M~ethods of Neutron Transport ( American Nuclear Society, La Grange Park, Ill., USA, 1993).18 12L. J. Shukovsky, "Dose, time, volume relationships in squamous cell carcinoma of the supraglottic larynx," Am J Roentgenol Radium Ther Nucl Med 108, 27-29 (1970).31 13J. J. Fischer, and J. E. Moulder, "The steepness of the dose-response curve in radiation therapy. Theoretical considerations and experimental results," Radiology 117, 179-184 (1975).30 Maximum electron energy CSDA range in liquid water 1.00E+00 1 .00E-01 -*RCSDA range of max. energy e- 1.00E-02 Thickness where 1% att. Y ~ IThickness where 2% att. S1.00E-03 _c: Thickness where 3% att. -*-Thickness where 5% att 1.00E-04 1.00E-02 1.00E-01 1.00E+00 1.00E+01 Energy(MeV) -1. Comparison between thicknesses required to attenuate different percentages of primary photons and Rcsn range for maximum energy secondary electrons produced by the same beam i~- -*-Ave.range of rnax. energy e- -.Range where% atten. Range where 2% atten. Range where 3% atten. - Range where 5% atten. Figure 1-2. Comparison between thicknesses required to attenuate different percentage of primary photons and average range for maximum energy secondary electrons produced by the same beam With proper model discretization of the angle-energy-spatial transport phase space, and at low photon energies yielding primary and secondary electron interactions with small electron Figure 1 Average range of maximum energy electron in liquid water EE E 1.00E-03 a 2 1.00E-04 o 1.0E0 cn0 0OE-02 1.00E-01 1.00E+00 Energy (MeV) 1.00E+01 Dose-rate(along Y-axis) 4.00E-05 3.50E-05 S3.00E-05 ~2.50E-05 S2.00E-05 1.50E-05 o 1.00E-05 5.00E-06 0.00E+00 --+- D(7 5-8 OMev) --* D(7 .0-7 5Mev) --b-- D(6 5-7 OMev) D(6 0-6.5Mev) -- D(5 5-6.0Mev) -*--D(5.0-5.5Mev) -- D(4 5-5.0Mev) ----(4 0-4 sMov) ---- D(3 5-4 OMev) -o-- D(3.0-3 5Mev) - D(2 5-3 OMev) ~-D(2 0-2 5Mev) SD(1 5-2 OMev) D(1 0-1 5Mev) D(0 5-1 OMev) tD(O O-o SMev) 0 10 20 30 40 50 Y(cm) Figure 4-15. Energy group absorbed dose rate distributions (Mev/ g.sec) along the y-axis using the EDK-SN method based on PENTRAN-MP Dose rate along the X-axis 4.0E-04 3.5E-04 --Dose using MnCNP5 S3.0E-04 pulse height tally(*F8) [" 1 I I-*EDK<-SN 6 2.5E-04 r 2.0E-04 1 .5E-04 ij 1.0E-04 5.0E-05 0.0E+00 0 10 20 30 40 50 X(cm) Figure 4-16. EDK-SN Total absorbed dose rate distributions along the X-axis compared to MCNP *F8 tally derived absorbed doses. MCNP uncertainties (20) average (6.0%) in the Source region. The average absolute relative difference between the two methods was 3.7% e" Positron (e+) Figure A-4. Pair Production in the Coulomb force field of the nucleus For hv values well above the threshold energy 2moc2, the electron and positron are strongly forward directed. Their average angle of departure relative to the original photon direction is roughly: moc2 B (radians) (A-30) Pair production in the electron field: In the kinematics of Pair Production in the electron filed (triple production), the photon divides its energy between the positron-electron pair produced and the host electron. The energy conservation equation becomes: hv = 2moc2 + T T- T,- (A-31) = 1.022M~eV + T + T + T The average kinetic energy of the three particles is -v h-1.022M2eVy T = (A-32) Compton Scattering: When passing through matter, photons can interact with the electrons of the matter. The Compton scattering describes an inelastic scattering of a photon with a free electron. If the electron' s initial binding energy is small compared to the photon energy (En << hV ) then the stuck electron can be considered initially free and unbound. Figure 10 shows this inelastic scattering process. The kinematics of this scattering (conservation of momentum and energy) lead to the energy of the scattered photon hv and to the energy of the scattered electron Es: E = hv Y p~~ hv=c, v PY hy'l',Y Figure A-2. Sketch of Compton scattering, photon of energy hv collide with an electron. The photon's incident forward momentum is hvl/c, where c is the speed of light in vacuum Based on the unbound and stationary condition of the electron, the Compton relation can be obtained by solving energy conservation equations T = hv hv (A-1) And conservation of momentum dose calculated via a Monte Carlo full physics absorbed dose in an equivalent model. This P factor is a kernel-specific factor, and once determined, it can be used for any execution implementing the corresponding electron kernels, provided grid spacing is consistently used. Photons Figure 4-1. Procedure to produce EDK based on Monte Carlo Calculations 4.3 EDK Determinations Based on The MCNP5 Monte Carlo Code In the MCNP531 Monte Carlo code, there are three different methods by which one can tally to yield absorbed doses to be used for integrated organ absorbed dose calculations; these are achieved by using either a cell flux tally (F4), an energy deposited tally (F6), or a pulse height tally (*F8), and each has a different basis for computing the absorbed radiation dosel4,31. To illustrate this, we compare absorbed doses computed using the F4, F6 and *F8 tallies, respectively to highlight the differences in these computations, and establish a procedure that ensures that our electron kernels are calculated accurately, and the Monte Carlo transport and subsequent absorbed doses determined as a benchmark are completely based on a detailed physics treatment. Figure 3-11i. CT view of the 11i-year-old male Series B original phantom collapsed into (a) 189 k voxels and the correspondent absorbed dose distributions (b) 900 k voxels and the correspondent absorbed dose distributions 4.4.4 Ordinate Projection In the discrete angle approximations, it is assumed that the EDK-SN is applied on distinct equal weight angles 04~(8, 9)a(u,r,,(). To generate these angles, the following relations should hold for any given octant: AB = (4-10) order The Norde in any 3-D EDK-SN quadrature corresponds to the number of levels from each direction cosine on the unit sphere, and there are Ndir= Norde (Norde + 2) ordinates on the unit sphere, with Ndir, oct oder rder+ 2)/ 8 in each octant, with Norder/2 distinct direction cosine values. Figure 4-10 shows Norde 42 used rendering 231 direction per octant. 8= (2i -1) (4-11) S=( )(4-12) 2 2i where i = N orde For any unit angular direction vector 0Z the ordinate must lie on the unit sphere: pU2 +2 +2 = 1 (4-13) pu = Sin (0)Cos (p) (4-14.a.) r = Sin (0)Sin (9) (4-14.b.) 5 = Cos (0) (4-14.c.) In triple production, the threshold for this process is 4moc2 =2.044 MeV. The higher threshold is required by conservation of momentum as shown in Table A-1 and also derived by Perrin (1933). Triple production becomes significant in radiotherapy treatment planning for energies greater than 14 MeV were it represents ~3% of the total photon interaction. Straggling and Multiple Scattering There is typically a distribution of angle and energy for charged particle passing the same medium. This is a result of multiple scattering and angular straggling which the spread in angular distribution observed in a population of initial identical charged particle after they have traverse equal path lengths. It will be somewhat exaggerated if the particle have passed through a layer of material, since multiple scattering then cause individual differences in scattering angle as well causing them to propagate in a conical angular distribution in accordance with Moliere's theory. However, for a discretized solid angle, one can consider an average angle for these scattered electrons on which approximate transport calculations can be made. The spread of the solid angle can be giving by Fermi angle for straggling especially when considering restricted stopping power including soft collisions, that yield small spread in the angle, and hard collisions resulting in 6 rays with energies less than the cutoff value A . Restricted stopping power (dT/ pdx~) Is the fraction of the collisions stopping power that include the soft collisions plus those hard collisions resulting in 6 rays with energies less than the cutoff value A . d T =kr I 2 (r + 2)2C =~~ k In+ rI- (A-33) One of the most important issues limiting the use of Monte Carlo techniques for clinical absorbed dose calculations is the long execution time required to reduce the statistical variance associated with the calculated absorbed dose at each differential volume. At the same time, the outcome of any radiotherapy procedure depends solely on the accuracy of the radiotherapy treatment planning system, so that a small change in the absorbed dose delivered (A 5%) can result in a dramatic change in local response of the tissue (f 20%)12,13 In order to reduce the statistical variance a, a ac associated with Monte Carlo, large numbers of histories N are required, and consequently this leads to higher running time to yield acceptable results. Improved variance reduction techniques may reduce the running time. However, such techniques are to be used with caution, as improper application of statistical variance reduction, method applied un-conservatively may cause wrong answers. Also, variance reduction is often tied to one objective; this can be problematic when faced with determining a dose at all locations in the phase space. The increasing use of voxelized human body geometries has created challenges in many aspects of applying Monte Carlo simulationS14. With increasing use of modern technologies such as CT and MRI in vivo imaging, anatomical definition has shifted from using simplified surface equations to segmented tomographic images in order to truly reflect the anatomy of a specific human body. The segmented images are represented in the form of a large uniform 3-D matrix of hexahedral voxels, where each voxel represents a constant material with a tag number (for example, CT number). One digital image may contain billions of voxels, making it very challenging be fully represented explicitly using Monte Carlo, particularly considering large models size and number of tally sites, adding significant overhead to the computation. Some special data management treatments are needed for Monte Carlo simulations to make 1 r15 O V(c,! O 1 I-15 o Y(cm) O Ocm 5cm 60 I. 60 40 40 20~ --: I 20 0 34 0 30 30 15 0 L 30 1 Y1e 0 0e 10 cm 15 cm Figure 2-5. 60Co y-ray 3-D Group 1 scalar flux distribution for x-level depths at 0 cm, 5 cm, 10 cm and 15 cm, computed by PENTRAN with the CEPXS gamma cross section library; an S42 angular quadrature (1848 directions) with P3 Scattering anisotropy for dosimetry applications may require fine energy group structures; it is clear the applicability of the cross section library should be investigated for each problem considered. Development of a standardized multi-group cross-section library for these kinds of applications would be of benefit to minimize potential inconsistencies in the application of currently available libraries. The results of parallel deterministic PENTRAN code were obtained within comparable running times to parallel MCNP5 Monte Carlo calculations for tally sites adj acent to the source. Overall, we demonstrated that with an adequate quadrature set, cell mesh interval, and cross-section library, it is possible to produce results with the PENTRAN discrete ordinates code that agree with the Monte Carlo prediction. A major advantage of the SN method is that it provides a detailed, accurate flux distribution at thousands of cells throughout the system, while the Monte Carlo method is best for providing accurate values at a few selected points near the source. Both methods can be advantageous in these types of applications, provided the level of SN discretization is adequate. This chapter laid out the foundations necessary for accurate SN discretizations in general, although only a homogenous phantom was considered. The need to implement heterogeneous media in the SN simulation is essential for any medical physics applications. In the next chapter, we set the foundation for the new PENTRAN-MP code system. In this system, we implement a state of the art series B phantom into our methodology by down-sampling these phantoms into accurate phantoms to generate an equivalent PENTRAN input. Moreover, we calculate the absorbed dose by applying a kerma approximation for diagnostic procedures. Then, the EDK-SN method is developed, presented, and applied in subsequent chapters. primary beam, and point-kernel methods do not properly handle tissue inhomogeneities. Moreover, since absorbed dose are generally tracked in all major organs in the body, variance reduction schemes for Monte Carlo are not all effective in this regard. The outcome of this dissertation is to demonstrate that absorbed dose from high-energy external beams radiation can be accurately computed for whole body and organ-specific absorbed doses. The EDK- SN method implements voxelized phantoms with discrete ordinates (SN) transport computations coupled with directional influences_and (pre-computed) full-physics Monte Carlo based electron absorbed dose kernels to yield total body absorbed dose information. This research shows that the deterministic techniques coupled with Monte Carlo based electron absorbed dose-kernels has significant potential for organ absorbed dose evaluation in the clinical management of radiation therapy patients. factors that were derived from the relative photon/electron density in the path through tissue as compared to that in water4'5. These methods, at their best, were useful for obtaining depth dose information in the center of a broad beam, for a medium consisting of homogeneous layers. Even then, these approaches do not account for the fact that photons interact not only with electrons, but also with the nucleus. For example, the elastic scattering cross section depends on the atomic number of the material. The clinical implication of this is that not all materials behave like scaled-density water; real bone, for example, scatters photons to greater extent than water of equivalent bone density due to the presence of Ca and other high Z elements not found in liquid water. The effect of this was better addressed in the "Pencil Beam" method. 1.3 Pencil Beam Transport. The Pencil Beam method sums the dose distribution from individual small diameter rays called pencil beams. A pencil beam is made up of particles which pass through a differential cross sectional area, 3x~y The off-axis dependence of the dose distribution for each pencil beam is described by thick-target multiple scattering'. The on-axis dose information is obtained from measured depth-dose data. The calculation of absorbed dose assumes that the patient or phantom can be represented by a stack of slabs of different material types. Each individual slab is homogeneous and infinite in lateral extent. Assuming a single pencil beam has normal incidence on such geometry, the absorbed dose at a position (x, y, z) can be expressed as D(x, y, z) = f (x, y, z)g(z) (1-1) where the first factor is the off-axis term, and the second factor is the central-axis term. The off- axis term may be represented by a Gaussian of the form 14B. Wang, G. Xu, J. T. Goorley, and A. Bozkurt, "ISSUES RELATED TO THE USE OF MCNP CODE FOR AN EXTREMELY LARGE VOXEL MODEL VIP-MAN," in The Monte Calrlo M~ethod: Versatility Unbounded In A Dynamic Computing World, edited by (American Nuclear Society, Chattanooga, Tennessee, 2005), pp. 32 1G. I. Bell, and S. Glasstone, Nuclear Reactor Theory, (Robert E. Krieger Publishing CO. Inc, Malabar, FL, USA, 1985).39 16B. G. Carl son, and K. D. Lathrop, Discrete Ordinates Angular Quadrature of the Neutron Transport Equation, (Los Alamos Scientific Laboratory Report, LA-3186, 1965).40 17D. Shedlock, and A. Haghighat, Neutron Analysis of Spent Fuel Storage Installation Using Parallel Computing and Advance Discrete Ordinates and Monte Carlo Techniques," Radiat Prot Dosimetry 116 662-666 (2005).45 IsK. D. Lathrop, "Remedies for Ray Effects," Nuclear Science and Engineering 45, 255-268 (1971).41 19G. Sjoden, and A. Haghighat, "PENTRAN A 3-D Cartesian Parallel SN Code with Angular, Energy, and Spatial Decomposition," in Int 'l. Conf: on Math. M~eth & Supercomp. for Nuclear App., edited by Saratoga Springs, New York, 1997), pp. 553.14 20A. Al-Basheer, M. Ghita, G. Sjoden, W. Bolch, and C. Lee, "Whole Body Dosimetry Simulations using the PENTRAN-MP Sn Code System," in American Nuclear Socity 2007 Annual2~eeting, edited by Bosten, MA, 2007), pp. 66 21G. M. Daskalov, R. S. Baker, R. C. Little, D. W. O. Rogers, and J. F. Williamson, "Two- Dimensional Discrete Ordinates Photon Transport Calculations for Brachytherapy Dosimetry Applications," Nucl. Sci. Eng. 134, 121-134 (2000).19 22G. M. Daskalov, R. S. Baker, D. W. Rogers, and J. F. Williamson, "Multigroup discrete ordinates modeling of 125I 6702 seed dose distributions using a broad energy-group cross section representation," Med Phys 29, 113-124 (2002).20 23M. L. Williams, D. Ilas, E. Sajo, D. B. Jones, and K. E. Watkins, "Deterministic photon transport calculations in general geometry for external beam radiation therapy," Med Phys 30, 3183-3195 (2003).5 24A. Haghighat, International Training Course/Workshop on M~ethodologies for Particle Transport Simulation ofNuclear Systems, 76 25(ICRU) 1980, RadiRRR~~~~~RRRRR~~~~atio quantities and units, (International Commission on Radiation Units and Measurements, Bethesda, MD, 1980).42 26F. H. Attix, "The partition of kerma to account for bremsstrahlung," Health Phys 36, 347-3 54 (1979).43 3D DETERMINISTIC RADIATION TRANSPORT FOR DOSE COMPUTATIONS IN CLINICAL PROCEDURES By AHMAD AL-BASHEER A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2008 cm3 VOxels, equivalent to the size of the EDK-SN voxel and the voxel size to be used in the SN calculation. X-ray Beam SSSSSSS SSBBSSS SBLLBSS Soft Tissue Lung Tissue Cortical Bone Figure 5-1. Three Monte Carlo phantoms used in the absorbed dose study. Phantom SSSSSSS is a homogeneous soft-tissue (' S') phantom; phantoms SSBBSSS and SBLLBSS are heterogeneous phantoms consisting of combinations of soft and bone tissue ('B'), and soft (' S'), bone ('B'), and lung ('L') tissue-equivalent materials, respectively After running the three phantoms for 100 minutes, Figure 5-2, using MCNP5 photon, electron mode and pulse height tallies (*F8) the results uncertainties (20) were 20% on average. A 'continue run' was executed for the three phantoms for 400 minutes more, for 500 minutes total, reducing the Monte Carlo uncertainty (20) down to 10% on average (Figure 5-3). However, the Monte Carlo computed difference in the deposited dose between the two tissues is within 2%, making it impossible to observe any difference between the different tissues in the phantoms, especially for therapeutic photon beams; this is shown in (Figure 5-4). The estimated time required to reduce Monte Carlo uncertainty to (20) 2% on-axis was 12500 minutes; hence, the time required to reduce the Monte Carlo uncertainty off-axis to 2% will be much greater. Radiative losses: (dT /pdx~)rad ZT (A-34) (dT / odx),,,, 820 Angular straggling (Fermi 1950) Electrons suffer multiple scattering events preclude straight transport in the material. Further complication are caused by the statistical fluctuation Figure A-5. Pair Production in the Coulomb force field of the nucleus 2= [N] Z,,,r24/ 0 (A-35) ao = 2 (A-36) moc As shown in the Figure 2-6, the PENTRAN results with S42 Legendre-Chebychev quadratures are in good agreement with the MCNP5 results, and differences between the two are within the statistical uncertainty of the Monte Carlo in air and the water phantom distal from the source. Overall, this figure shows that increasing the angular quadrature order results in a progressive improvement in solution accuracy, and that a substantial number of directions are required for proper representation in air region between the source and the water phantom. We further note that while the S12 leVel-Symmetric quadrature is quite adequate to represent neutron scattering in nuclear engineering neutronics and reactor physics applications, it does not produce accurate results here. For a general problem with a high degree of angular dependence, the ray- effects must be largely mitigated by employing a large number of directions, i.e., a high-order angular quadrature, From our analysis, a minimum of S32 Legendre-Chebychev quadratures should be considered. Investigation of alternative, less computationally intensive methods for mitigating the ray-effects in problems of this type has shown that ordinate biasing methodologies can also be ofbenefitl9 1.OOE+10- 1.OOE+08- -HSn22 Sn32 -MSn42 1.OOE+06 MMCNP5 1.OOE+04- 1.OOE+02- 1.OOE I O O.OOE+OO 3.OOE+01 6.OOE+01 9.OOE+01 Distance (cm) Figure 2-6. Scalar flux distributions using S12, S22, S32, and S42 with P3 Scattering anisotropy in PENTRAN using the BUGLE-96 gamma library study. This geometry data was then collapsed into a 4.7 mm x 4.7 mm x 30.0 mm voxel size (79 x 48 x 50 voxels) model using GHOST-3-D. Subsequently, using PENMSH-XP, the collapsed model was cast into a 3-D spatial distribution and sub-divided into six coarse mesh z-levels containing five coarse meshes in x-y domain, and these contained 6320 fine mesh cells (Figure 3-4, 3-5). An isotropic volumetric (17.8 cm x 3.65 cm x 30 cm) x-ray source was placed above the upper left side of the phantom. Although the model has relatively small fine mesh sizes, using S32 Legendre-Chebychev quadrature, the adaptive differencing, and CEPXS cross-sections with P3 aniSotropy, we were able to eliminate any unphysical oscillations that may result in the deterministic solution. Using 4, 6, and 8 energy groups, the SN model was executed on 16 Opteron processors and was completed in 1.8, 2.1 and 3.8 hours, respectively. With the general purpose radiation transport code MCNP5, we simulated an equivalent model to the one simulated by PENTRAN. The collapsed 3-D binary files composing the UF voxel phantom was converted to ASCII formatted matrix files and then incorporated into a MCNP5 lattice defined input. An input deck example for the MCNP5 model is presented in Appendix B. Volumetric flux (F4) mesh tallies were used in the Monte Carlo geometry and were equivalent to the discretized SN volumes all over the model. All simulations were performed using the parallel cluster Einstein composed of 16 processors at 2.4 GHz with 4 Gb/processor, for a total of 64 Gb of parallel memory used for general purposes. Parallel MCNP5 results in the source region converged with a statistical error less than 4% in the source region on average. Dose volume histograms (DVHs), for a selected number of organs giving the percentage of tissue volume (organ) that receives an absorbed dose greater than the value indicated on the abscissa, were calculated and compared to those results from the equivalent Monte Carlo model. As anticipated, DVH results are in good agreement for the SN and Monte Carlo radiation transport simulations (Figure 3-15). However, Monte Carlo calculations can require substantial run times for low precision at "out of field" distances in the phantom, particularly as photons down-scatter to lower energies. The Monte Carlo results computed values will eventually converge with increased run time as more particles reach sites distal to the source. DVH (Left Lung) DVH-(RIGHT LUNG) In-field Out-of-field 100.00% k I100.00% 75.00% i 1 75.00% S50.00% DVH(L Lung)-SN 50.00% -t DVH-(R Lung)-SN SDVH(L Lung)-MCNP DVH-(R Lung)-MCNP 25.00% i L 125.00% 0.00% II ... .00 ,**** 0.00 0.40 0.80 1.20 1.60 2.00 0.00 0.10 0.20 0.30 0.40 0.50 Dose (mGy) Dose (mGy) a. b. Figure 3-15. Dose Volume Histograms (DVH) computed based on PENTRAN (S32,P3) transport results using the CEPXS library vs. DVHs computed using MCNP5 Monte Carlo results. MCNP5 Monte Carlo results required 16 hours 16 processors compared with the PENTRAN-MP SN solution, requiring <1.9 hours on 12 processors for the whole-body computation 3.5.3 Parallel Computation For the UF 11 year old male phantom problem, a decomposition strategy was used to assist the parallel performance. In Table 3-2, we present the cumulative problem time and memory G2(70-80)key 1 81E+06 1 51E+06 MCNP5 1 21+06- *PENTRANI 9.10E+05 610E+05 1.00E+04 -* 0 5 10 15 20 25 Y(cm) G4(50-60)key 6.01E+06 5101E+06 ., 4.01E+06 i"IMCNP5 x -PENTRAN q3.01E+r06 2.01E+06 1.01E+06 1.00E+04 0 5 10 15 20 25 Y(cm) G1(80-90))key 6.01E+05 4.51E+05 -MCNP5 30E0 ~PENTRAN 1 51E+05 1 00E+03 " 0 5 10 15 20 25 Y(cm) G3 (60-70) key 3.51E+06 2 81E5+06 2 11E+0 i MN k -PENTRAN 1.41E+06 7.10E+05 1.00E+04__ 0 10 20 30 bicm 50 20'e" ,-.,2 Figure 3-6. x-ray 3-D Group 1, 3, 5 and 7 scalar flux computed by PENTRAN with the CEPXS cross section library, S32 angular quadrature (1088 directions) with P3 Scattering amisotropy Figure 3-7. Flux distributions using S32 with P3 scattering anisotropy in PENTRAN using the CEPXS library vs. MCNP5 Monte Carlo results for energy groups 1, 2, 3 and 4 along the source central axis. Monte Carlo uncertainties were less than 5%, 4%, 4%, and 3% respectively Therefore, this is an example that highlights that tissue heterogeneities can be difficult to discern in these problems using the Monte Carlo method. As explained previously, one of the most important issues limiting the use of the Monte Carlo techniques for clinical absorbed dose calculations is achieving a statistically reliable result in a reasonable time. This problem is even greater when a photon-electron mode is required, when CPE is not applicable, causing computer running times needed to reduce the statistical error to increase rapidly. Improved variance reduction techniques may reduce the computation times. However, such techniques are to be used with caution, as improper application of statistical variance reduction may cause wrong answers. Moreover, variance reduction techniques may not reduce the time required to reduce the variance to acceptable levels, especially when accurate absorbed dose estimations are needed over the entire problem. 100 Minutes (Pulse height tally using photon- 4.00E-04 -electrn mode) 3.50E-04 --*-- SBLLBSSS Ilu -o- SSBBSSS (0 3.00E-04 as 2.50E-04- 1.50E-04- ur1.00E-04 - 5.00E-05- 0.00E+00 O 10 20 30 40 50 X(cmR) Figure 5-2. Monte Carlo derived total absorbed dose rate distributions along the X-axis computed using *F8 MCNP5 tally for three different phantoms. MCNP5 uncertainty (20) average (20.0%), 100 Minute run Although the model has relatively small fine mesh sizes, using Legendre-Chebychev S32 quadrature, the DTW differencing scheme, and CEPXS cross-sections with anisotropy level P3, we were able to eliminate any unphysical oscillations that may result in the deterministic solution. Using 16 energy groups, the SN model was executed on 16 processors and completed in 2.0 hours. With the general purpose radiation transport code MCNP5, we simulated an equivalent model to the one simulated using PENTRAN. The collapsed 3-D binary files were converted to ASCII formatted matrix files and then incorporated into an MCNP5 lattice defined input; this matrix was also generated by our down-sampling code GHOST-3-D. 5.2.1 EDK-SN Absorbed Dose Calculations Calculations were performed to calculate the SN results for the scalar fluxes for all energy groups. Using EDK-SN, 3-D absorbed dose distributions were calculated to each voxel in the entire 15 year old male phantom model for a broad-beam external photon field. Figure 5-11 shows the stages of calculation in the PENTRANP-MP code system that resulted in 3-D Absorbed dose distributions for the entire phantom. Figure 5-12 shows absorbed dose contributions from selected energy groups, and the total absorbed dose resulting from adding all energy group-absorbed dose contributions, as prescribed by the EDK-SN procedure. Organ absorbed doses using EDK-SN, for a selected number of organs, were calculated and compared to Monte Carlo predictions, and are presented in Table 5-1. The total execution time for the Monte Carlo was ~4 hrs (wall clock time) on 16 processors, while the total wall clock execution time for the EDK-SN was ~ 2hrs on the same number of processors. As anticipated, the absorbed doses for the organs located directly in the x-ray field (right+ left lung) or partially in the x-ray field (liver) were in good agreement for both the EDK-SN and MCNP5 pulse height tally data. This balance of neutral particles is cast as the linear Boltzmann Transport Equation (BTE). The BTE is an integro-differential equation which describes the behavior of neutral particles in terms of spatial, angular, and energy domains as they interact in a system; the steady-state transport equation is given in Equation (1-5): dE 'WrR dn' o-, (,',E 'a) Yr, ', ')qE,0E (15 0 4xr where y(r7, OE)d'rdEdai is the angular flux with energy between E an~d E dE, at position r within the volume element dF~, and in direction aZ within the solid angle daZ. The scattering term, defined with the double integral term, represents the sum of the particles scattered into dEdaZ from all the dE 'dnZ' after scattering collisions represented by the double differential cross-secti on Gs (F,ER 4 E, iZ' 4 6) q (r-, 6,E ) represent the external source. The left side of Equation (1-5) represents streaming and collision terms (loss), and the right side represents scattering and independent sources (gain). Since it describes the flow of radiation in a 3-D geometry with angular and energy dependence, this is one of the most challenging equations to solve in terms of complexity and model size. For rendering a deterministic computational solution for a large problem in a reasonable time, it is necessary to utilize a robust parallel transport algorithm and a high performance computing system. All the methods described below are different approaches to solving this equation. For cases with charged particle transport, the Boltzmann and Fokker-Plank formalisms (Boltzmann-FP) 10are generally used. A simplified form of this equation is usually used, which is: Collisional kerma, Kc can be defined in terms of the related net energy transferred, Sn ,26,27 and R25 as. dE" /dm (1-10) The energy-transferred in a volume Vis Et;' = (R~n), (Rozer)"" on R + c Er, R', (1-11) in which (R~n), = the radiant energy of uncharged particles entering the volume y, (Rozer) onr= the radiant energy of uncharged particles leaving V, except that which originated from radiative losses of kinetic energy by charged particles while in V. Rr = radiant energy emitted, where u indicates uncharged particles. Charge particle equilibrium (CPE) can be established in a volume when the energy carried into the volume V by charged particles is balanced by the energy carried out of the volume V by the charged particles. Note that, CPE is not established for thin materials irradiated by photons. The range of secondary electrons produced is larger than the material dimensions, so that much of the electron energy is lost simply by electrons leaking out of the material. Under the conditions of CPE, Attix demonstrated that when the energy imparted equals energy transferred, r = 8,,, collisional kerma Kc is equal to the total dose Dar,2 Eq(11). CPE~iP(-2 methodology, we were able to accurately yield a solution analogous to the MCNP5 Monte Carlo calculations using the CPE principle for low energy beams. . We have demonstrated; that by using electron kernels for the different materials (cortical bone, soft tissue and lung) the EDK-SN methodology incorporating a discrete ordinates (SN) approach is an effective alternative to the Monte Carlo method for accurately estimating the absorbed dose for high energy photon beams over an entire human phantom. In doing so, we simulated three phantoms using MCNP5 and PENTRAN. Each phantom had a flat weighted 8 MV x-ray source impinging on a 45x45x45 cc soft tissue phantom. Two of these phantoms have slabs of 1-cm-thick discs of materials representative of bone, or lung tissue. The comparison of the MCNP5 results and the EDK-SN showed good agreement between the two methodologies. We also simulated the complete human phantom using 270,540 voxels, and demonstrated that the EDK-SN method readily provides accurate whole body 3-D absorbed dose distributions. The angular data provided by the SN method was used to guide the absorbed dose mapping process for secondary electrons to the surrounding voxels; this makes our approach truly novel. We were able to accurately yield a solution analogous to the MCNP5 Monte Carlo calculations in the voxelized phantom. We note that that tally results from the 16 energy groups, used to represent the radiation spectrum, from the CEPXS multi-group cross-section library, were compared to continuous energy MCNP5 Monte Carlo results using the ENDF/B-VI data libraries. Parallel deterministic PENTRAN-MP results were obtained within comparable running times to parallel Monte Carlo calculations for tally sites adjacent to the source. An equivalent Monte Carlo calculation to obtain the absorbed dose distributions using photon, electron mode with pulse height tallies converged to reasonable statistical uncertainty over the entire problem Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy 3D DETERMINISTIC RADIATION TRANSPORT FOR DOSE COMPUTATIONS INT CLINTICAL PROCEDURES By Ahmad Al-Basheer August 2008 Chair: Glenn Sjoden Maj or: Nuclear Engineering Sciences The main goal of this dissertation was to establish the feasibility of basing megavoltage external photon beam absorbed dose calculations in voxelized phantoms on SN deterministic calculations and pre-calculated electron absorbed dose kernels derived from full-physics Monte Carlo. The SN derived electron absorbed dose kernel method EDK-SN, developed as part of this research, achieves total execution times that are on the order of several times to orders of magnitude faster than conventional full-physics Monte Carlo electron transport methods considering equivalently detailed models and data Eidelity. With the rapid movement toward intensity modulated radiation therapy (IMRT), radiation beam intensities have increased dramatically over the past decade, thus heightening the need for further characterization of out-of-Hield organ absorbed doses, along with their associated biological risks. Assessment of these tissue absorbed doses is complicated by two fundamental limitations. First, anatomic information on the patient is generally restricted to a partial body CT image acquired for treatment planning; consequently, whole-body computational phantoms must be employed to provide the out-of-field anatomy model structure for absorbed dose evaluation. Second, existing methods based on Monte Carlo radiation transport, even with the application significant variance reduction, are quite computationally inefficient at large distances from the a is the angle between the y-axis and the line of nodes . p is the angle between the z-axis and the Z-axis . y is the angle between the line of nodes and the Y-axis. The above description works for the z-y-z form. Similar conventions are obtained just renaming the axis (zxz, xyx, xzx, yzy, yxy). There are six possible combinations of this kind, and all of them behave in an identical way to the one described before. Just imagine three rotations over the set of axis forming the rotation matrix R. [R ]= [Rz (a)][R, (/7)][Rz (7)] (4-6) Where R is the rotational element matrix Depending on the convention used for the calculating these angles, they can be related to the angles for the net current obtained from the PENTRAN run. Assuming that rotated matrix indices to be known, we generated the vector A for each point from the center of the matrix, where the dose values of this vector form the matrix dose point to be projected onto the SN matnx : [A '] =[R][A ] (4-7) [R -] [A '] =[R1] [R ][A ] (4-8) [A ] =[R ][A '] (4-9) When calculating the original un-rotated A vector for each point, we interpolate the vector location, in 3-D, for the corresponding dose value based on the original EDK matrix calculated with the Monte Carlo model. This generates a 3-D translated (by rotation) electron absorbed dose kernel matrix corresponding to a particular direction in space. For the current code system, this process was performed for 1848 directions, generating 1848 files for each energy group. An example of a rotation of the EDK from a reference (Monte Carlo) EDK is depicted in Figures 4- 9.a and 4-9.b. $Ears $Esophagus $External nose $Eye balls $Gall Bladder W $Gall Bladder C $Heart W $Heart C $Kidney-cortex (L) $Kidney-cortex (R) $Kidney-medulla (L) $Kidney-medulla (R) $Kidney-pelvis (L) $Kidney-pelvis (R) $Larynx $Lens $Liver $Lung (L) $Lung (R) $Nasal layer (anterior) $Nasal layer (posterior) $Oral cavity layer $Pancreas $Penis $Pharynx $Pituitary Gland $Prostate $Rectosigmoid W $Rectosigmoid C UF_15yr voxel MCNP model whole body Matrix size [60,27,167] Voxel resolution=1.0* 1.0 1.0 2008-7 Ahmad Al-Basheer The University of Florida WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW APPENDIX B MCNP5 PHNTOM INPUT DECK cells Card Body composition and density 1 -0.001205 -999 imp:p,e=1 u=1 vol=78406.867 $ Air Residual Soft Tissue Adrenal (L) Adrenal (R) Brain Breast Bronchi Right Colon W Right Colon C -1.03 -1.03 -1.03 -1.04 -0.94 -1.07 -1.03 -0.83 -1.10 -1.03 -1.05 -1.07 -1.03 -1.03 -1.04 -1.06 -1.05 -1.05 -1.04 -1.04 -1.04 -1.04 -1.07 -1.07 -1.06 -0.24C -0.24C -1.02 -1.02 -1.02 -1.02 -1.05 -1.03 -1.03 -1.03 -1.03 -0.49 i i i -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 0 -99! 0 -99! -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 mp:p,e=1 u=2 mp:p,e=1 u=3 mp:p,e=1 u=4 I I I I vol=3.94E+04 $ vol=2.00E+00 $ vol=3.00E+00 $ vol=1.39E+03 $ vol=1.70E+01 $ vol=4.00E+00 $ vol=9.70E+01 $ vol=1.45E+02 $ Svol=3.00E+00 vol=2.40E+01 vol=2.00E+00 Svol=1.10E+01 Svol=0.00E+00 Svol=4.50E+01 Svol=1.86E+02 vol=4.23E+02 Svol=8.00E+01 Svol=8.40E+01 Svol=3.10E+01 vol=3.20E+01 vol=5.00E+00 Svol=6.00E+00 Svol=2.20E+01 Svol=0.00E+00 Svol=1.23E+03 7 vol=1.72E+03 8 vol=2.04E+03 Svol=0.00E+00 Svol=0.00E+00 vol=0.00E+00 Svol=1.08E+02 Svol=2.30E+01 Svol=2.00E+00 Svol=0.00E+00 Svol=4.00E+00 Svol=0.00E+00 Svol=1.39E+02 imp:p,e=1 u imp:p,e=1 u imp:p,e=1 u imp:p,e=1 u imp:p,e=1 u Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 9 imp:p,e=1 9 imp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 Simp:p,e=1 1=5 1=6 1=7 1=8 1=9 u=10 u=11 u=12 u=13 u=14 u=15 u=16 u=17 u=18 u=19 u=20 u=21 u=22 u=23 u=24 u=25 u=26 u=2 u=2 u=29 u=30 u=31 u=32 u=33 u=34 u=35 u=36 u=37 u=38 Table 1-2. Approximate attenuation of photons with a layer of water equal to the maximum range of secondary charged particles28 Primary Radiation Energy Maximum Secondary Photon Attenuation (%i) inl (MeV) Electrons Range (cm) Maximum Electron Range 0.1 0.01387 0.23% 1.0 0.43 3.00% 4.933 10.04% 9.305 15.% 1.5.2 High Energy Beam Dose Calculations Based on the previous discussion, as the energy of the ionizing radiation increases, the penetration power of the secondary charged particles increases more rapidly than the penetration power of the primary radiation in a given volume. This leads to CPE failure, caused by the small number of charged particles produced at deeper depths in volume V, relative to the number of charged particle produced at initial points along the primary radiation field direction. This is indicative of significant attenuation of the primary indirectly ionizing radiation. As a result of this phenomenon, higher numbers of charged particles will be generated close to the entrance compared to the rest of the volume, leading to CPE failure considering dose in that volume. The degree of this failure becomes progressively large for higher energies, so that D no longer equals to Kc A more detailed discussion on this topic is presented in chapter 4. 1.5.3 Limits of Charged Particle Equilibrium Assuming that a maximum energy is transferred from a photon to an electron, we can use this assumption to investigate a CPE requirement for a water phantom region, as defined earlier. Considering the continuous slowing down approximation range of the electron (REsn), and average range of the maximum energy electrons Ross =IT dTL (1-13) 0.5 Of 00,,0 UU r S0O 000 00 0 U CO 0 00 1 1 Figure 4-10. First octant" of S42 equal weight angles used for proj ecting the EDK-SN 4.4.5 Final Dose Calculations The EDK-SN model specific portion, described in Figure 4-11, serves as the critical link in the EDK-SN system to accumulate the absorbed dose in each fine mesh. The algorithm determines the accumulated energy deposited in each voxel for each photon energy group, based on the SN net current vector and is scaled to the appropriate dose by the photon flux. Finally, the code accumulates all the absorbed dose contributions from all the energy groups in to a 3-D total absorbed dose matrix. Figure 4-11 shows a flow chart representing the EDK-SN methodology. An energy dependent EDK was generated using high resolution direction proj section sets and a MC pre-computed EDK for a single direction. These kernels were generated over the unit sphere of ordinates via the Euler Angles described in Section 4.4.3. Concurrently, when generating the EDK, a FORTRAN90 module will be generated for later use in the dose calculation portion of the 5.1.4 Correction Factor Approach ................. ...............109........... .. 5.2 Computational Phantom of 15-Year Old Male ................ .............................111 5.2.1 EDK-SN Absorbed Dose Calculations ................. ............... ......_. ...114 5.2.2 Ray-Effects ................ ...... ......_.._ ....... .. ....... ... ....... ...............1 6 5.2.3 MCNP5 Energy Deposited Tally (F6) and Pulse Height Tally (*F8) Results ................. ...............117____....... 5.3 Discussion and Findings ............ _...... ._ ...............118.. 6 SUMMARY, CONCLUSIONS AND FUTURE WORK .............. ..... ............... 12 6.1 Conclusions and Findings .............. ...............121.... 6.2 Future W ork ............ _...... ._ ...............123... APPENDIX A HIGH ENERGY PHOTON TREATMENT (IN CONDENSED ELECTRON TRANSPORT) ................. ...............124................ B MCNP5 PHNTOM INPUT DECK .............. ...............133.... LI ST OF REFERENCE S ................. ...............136................ BIOGRAPHICAL SKETCH ................. ...............140......... ...... dose values were measured in water, a table of the percentage of the maximum dose as a function of depth (on-axis) in the water could be deduced. A similar table could be generated for off-axis ratios (OARS) by measuring lateral profiles at a variety of depths. The absorbed dose at a given point P in the patient could then be estimated by knowledge of the maximum dose for the treatment, and applying a percent depth dose (PDD) correction for the depth of P in the patient, and an off-axis correction for the distance to P from the axis at that depth4 Tables of PDDs were made available for various sources (incident energies), field sizes, depths in the phantom, and SSDs. To modify the dose for an SSD that was different from the reference SSD, the ratio of inverse distances squared would be applied. For a change in field size, tables of field-size-factors would be subsequently applied. To compensate for the curved surfaces of a patient, as opposed to the flat surface of the phantom, isodose lines could be shifted along rays proj ected from the source. It is important to point out that these techniques are often quite successful for predicting absorbed dose in a homogeneous volume. In fact, in some cases they may be superior to alternative simplified methods (for homogeneous volumes) because they are based on measurements from a particular machine; thus, they circumvent the problem of how to model the photon distributions emerging from the source and incident on the patient/phantom, which can be a significant obstacle in accomplishing computational radiation transport based models. The real challenge for these empirical-based methods, then, lies in correcting the homogeneous water absorbed dose distribution to account for the many inhomogeneities present in the human body, including air cavities, bone, dense muscle mass, and sparse lung tissue. One way of addressing this problem is to apply "equivalent-thickness corrections" to the water measurements. These corrections were based on compressing or stretching isodose lines by Aside from the geometric assumptions of infinite homogeneous slabs, the pencil beam model does not handle changes in scatter from lateral heterogeneities; deviation of up to 5% can be expected for low energy photon beams, and 14% for high energy photon beams8,9 1.4 Solving the Boltzmann Equation At the most fundamental level, radiation transport simulation involves determining the distribution of the particles (in physical space, direction, time, and energy) after they have traversed some distance into a medium. Formally, a way to begin this discussion is by writing the general equation for particle transport from kinetic theory of gases. This theory does not explicitly treat individual particles, but rather describes the aggregate behavior of a continuous particle distribution in a statistical manner. Let n be the expected number of particles within a differential volume element d3r at position r with a velocity (which is directly related to its kinetic energy) that is within a differential velocity d3v, at time t. Then, if a general solution is found for the "particle density function" n(r, v, t) in phase space, then the transport problem is solved. When transporting the radiation through an arbitrary volume V, bounded with a surface area S, the rate of change of number of particles is due to flow out through the surface of V, collision events within V, or sources (scattering and independent sources) of particles located inside of V. (all the terms of rate of change; commonly the first two terms are losses and the last term is production) chan ge in )= to flow + to collisionsI + to sources nat time t houhSin V in V BIOGRAPHICAL SKETCH Ahmad Al-Basheer was born in 1978 in Ramtha City, Jordan, in the heart of the Middle East. Generations of his family have made Ramtha their hometown for centuries. This city, which borders Syria and Jordan, is the gateway to the desert region of the Middle East. Ahmad is the oldest of eight children, and his parents, Khaled and Fatima, are both involved in educational careers; both are school principals. Ahmad was involved not only in traditional educational venues, but also in chess, soccer, reading and poetry. Ahmad competed in several reading and poetry contests throughout his schooling; he won the National Award for high school students for Best Poet. Ahmad graduated from Jordan University for Science and Technology with a Bachelor of Science degree in applied physics in 2000. Upon graduation, he spent one year as a physics teacher at the Yarmouk University Model School. In January, 2002, Ahmad was granted admission to the graduate program in the Department of Nuclear and Radiological Engineering at the University of Florida. He performed his research work under the supervision of Dr. Samim Anghaie, where he investigated a new approach for reducing scattered photons and electron contamination in Cobalt-60 therapy beam. In fall 2004 he j oined the University of Florida Transport Theory Group (UFTTG) to pursue his Ph.D. with Dr. Glenn Sj oden. During his Ph.D., He won several honors, including best papers presented to number of conferences, a number of invited speeches at national conferences, and several recognition form number of organization for his achievements. Ahmad married Ghadeer Al-Basheer, on June 16th 2007 and has one daughter, Jude Al-Basheer who was born in June 14th 2008. Solving the Boltzmann equation directly using deterministic techniques can be significantly more efficient than using the Monte Carlo methods, and therefore it can be used to improve the quality of radiotherapy planning. Relative to the amount of information generated, and with the use of HPC resources, such techniques usually require significantly lower running times (relatively to Monte Carlo) to provide detailed and precise solutions. Moreover, the rapid proliferation of Computed Tomography (CT) and Magnetic Resonance Imaging (MRI) enable human anatomy to be represented using regular images consisting of rectangularly (hexahedral) shaped voxels in a regular matrix. These new developments make Cartesian geometry deterministic techniques ideal for medical physics applications, especially given the demand for high resolution absorbed dose distributions. This research is focused on the development of EDK-SN, a first-of-a-kind capability to assess organ absorbed doses incurred in radiation therapy (RT) applications using anatomical patient models coupled with 3-D deterministic discrete ordinates (SN) radiation transport methods. Absorbed doses delivered as a result of secondary charged particles will be linked using electron kernels to account for secondary electron effects. This methodology is proposed for regions in the phase space where Charged Particle Equilibrium (CPE) is not satisfied, and subsequent absorbed dose estimation using a collisional kerma is not valid. 1.2 Initial Clinical Dose Evaluations In the field of radiation therapy, absorbed dose distributions were historically first approximated without doing any explicit radiation transport at all, but rather by extrapolating from measured results taken for a particular energy, field size, and source-to-surface distance (SSD) on a given machine. An extrapolation was performed by modifying the absorbed dose that was experimentally determined under the reference parameters by ratios representing the change between specific parameters and the parameters of interest. For example, if absorbed Figure 4-2. Monte Carlo model used for comparisons between MCNP5 pulse height tally, energy deposited tally and cell flux tally employing dose conversion factors (DE, DF). The source is located at the dose driving voxel DDV 0 202 v G ui l U l a~ illll IIvIFIl v utig~ S0.15--E o. 10 Cell fact 8 0.05 0.00 0 5 10 15 Depth (cm) Figure 4-3. Comparison between F6 (p, e mode), F6 (p mode), and F4 (DE, DF) tallies in MCNP5, statistical uncertainty was less than 1% |

Full Text |

PAGE 1 1 3D DETERMINISTIC RADIAT ION TRANSPOR T FOR DOSE COMPUTATIONS IN CL INICAL PROCEDURES By AHMAD AL-BASHEER A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2008 PAGE 2 2 2008 Ahmad Al-Basheer PAGE 3 3 PAGE 4 4 ACKNOWLEDGMENTS Any project, no matter how individualized, will almost certainly require input, assistance or encouragement from others; my project is no exception. I extend my sincerest gratitude to Dr. Glenn Sjoden. He has been a constant source of thoughtful guidance in pursuing this project. Because of his input, advice and challenge, I have matured as a researcher and as a student. I would also like to acknowledge the contri bution of my supervisory committee: Dr. Ali Haghighat, Dr. James Dempsey, Dr. Wesley Bolch, Dr. Ayman Hawari a nd Dr. William Ditto. I would also like to take this opportunity to express my sincer e gratitude to my classmates and the faculty members in the Nuclear and Radiological Engineering and Radiation Oncology Departments for their constructiv e suggestions and helpful comment s. During my research, I had the pleasure of working with the University of Florida Transport Theo ry Group (UFTTG). With the help of coworkers and friends, even the tedium became bearable. Special thanks go to my family, and mostly my wife Ghadeer, Mom and Dad; they have definitely been encouraging me for the longest! They will always be dear to my heart. I thank God for all He has done for me. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS...............................................................................................................4LIST OF TABLES................................................................................................................. ..........8LIST OF FIGURES.........................................................................................................................9ABSTRACT...................................................................................................................................14CHAPTER 1 INTRODUCTION AND LITERATURE REVIEW..............................................................161.1Overview.......................................................................................................................161.2Initial Clinical Dose Evaluations.................................................................................. 171.3Pencil Beam Transport.................................................................................................. 191.4Solving the Boltzmann Equation.................................................................................. 211.4.1Monte Carlo Method.........................................................................................231.4.2Deterministic Methods...................................................................................... 251.4.2.1Discrete ordinates method (SN method)..............................................251.4.2.2Monte carlo versus deterministic SN methods.................................... 271.5Dose and Charged Particle Equilibrium........................................................................ 281.5.1Deterministic Dose Calculations for Diagnostic Application........................... 301.5.2High Energy Beam Dose Calculations.............................................................. 311.5.3Limits of Charged Particle Equilibrium............................................................ 311.6Outline of This Research............................................................................................... 352 CRITICAL DISCRETIZATION ISSUES IN 3-D DISCRET E ORDINATES MEDICAL PHYSICS SIMULATION................................................................................... 362.1Code Systems................................................................................................................372.2Simulation Geometry.................................................................................................... 382.2.1PENTRAN Model............................................................................................. 382.2.2Monte Carlo Model...........................................................................................392.3Comparison of Deterministic and Monte Carlo Results............................................... 392.3.1Investigation of Angular Quadrature................................................................392.3.2Investigation of Spatial Discretization (Meshing)............................................ 452.3.3Investigation of Cross-Section Libraries:..........................................................492.4Absorbed Dose Calculations.........................................................................................522.5Parallel Decomposition................................................................................................. 542.6Discussions and Findings.............................................................................................. 55 PAGE 6 6 3 WHOLE BODY DOSIMETRY SIMULATIONS WITH THE UF-B SERIES PHANTOM USING THE PENTRAN-MP CODE SYSTEM FOR LOW ENERGY PHOTONS........................................................................................................................ ......573.1Methodology.................................................................................................................583.2UF Series B Pediatric Phantoms...................................................................................593.3PENTRAN-MP Code System....................................................................................... 613.3.1Pre-processing Codes........................................................................................613.3.1.1The collapsing code (GHOST-3-D)....................................................613.3.1.2The mesh generator code (PENMSH-XP)..........................................613.3.1.3Source spectra and cross-sections generation..................................... 623.3.2Transport Calculation........................................................................................ 623.3.3Post-Processing Codes......................................................................................633.3.3.1Data management utility (PENDATA).............................................. 633.3.3.2Dose calculation code (3-D-DOSE)...................................................633.4Computational Example................................................................................................ 643.5Computational Results.................................................................................................. 673.5.1Flux Comparison of Deterministic and Monte Carlo Results........................... 673.5.2Dose Calculations and Dose Volume Histograms (DVH)................................ 703.5.3Parallel Computation.........................................................................................753.6Discussion and Findings...............................................................................................764 ELECTRON ABSORBED DOSE KERNELS TO ACCOUNT FOR SECONDARY PARTICLE TRANSPORT IN DE TERMINISTIC SIMULATIONS .................................... 794.1Introduction to The EDK-SN Method........................................................................... 794.2Monte Carlo Based Dose Kernels.................................................................................804.3EDK Determinations Based on The MCNP5 Monte Carlo Code................................. 824.3.1Computational Model........................................................................................ 834.3.2MCNP5 Tallies Comparison............................................................................. 834.4EDK -SN Code System Methodology........................................................................... 864.4.1Net Current........................................................................................................ 864.4.2MCNP Based Electron Kernels.........................................................................874.4.3Euler Angles......................................................................................................894.4.4Ordinate Projection........................................................................................... 924.4.5Final Dose Calculations.................................................................................... 934.5External Source Benchmark..........................................................................................954.6Comparison of Deterministic and Monte Carlo Results............................................... 964.7Parallelization of EDK-SN Code...................................................................................994.8Discussion and Findings.............................................................................................1015 APPLICATION OF EDK-SN TO MODELS WITH HETEROGENITIES......................... 1035.1Heterogeneous Tissue-Equivalent Phantoms.............................................................. 1035.1.1Monte Carlo Simulations................................................................................1035.1.2PENTRAN Model........................................................................................... 1075.1.3EDK-SN Absorbed dose Comparison.............................................................. 108 PAGE 7 7 5.1.4Correction Factor Approach............................................................................ 1095.2Computational Phantom of 15-Year Old Male........................................................... 1115.2.1EDK-SN Absorbed Dose Calculations............................................................ 1145.2.2Ray-Effects......................................................................................................1165.2.3MCNP5 Energy Deposited Tally (F6) and Pulse Height Tally (*F8) Results............................................................................................................. 1175.3Discussion and Findings.............................................................................................1186 SUMMARY, CONCLUSIONS AND FUTURE W ORK.................................................... 1216.1Conclusions and Findings...........................................................................................1216.2Future Work................................................................................................................123APPENDIX A HIGH ENERGY PHOTON TREATMEN T (IN CONDENSED ELECTRON TRANSPORT) ......................................................................................................................124B MCNP5 PHNTOM INPUT DECK......................................................................................133LIST OF REFERENCES.............................................................................................................136BIOGRAPHICAL SKETCH.......................................................................................................140 PAGE 8 8 LIST OF TABLES Table page 1-1 Comparison between Monte Carlo me thods and determ inistic techniques....................... 28 1-2 Approximate attenuation of photons with a layer of water equal to the maxim um range of secondary charged particles28..............................................................................31 1-3 Approximate thickness of water need ed to achieve various attenuation amounts com pared to maximum energy electron ranges generated by the same beam................... 32 2-1 Total number of directions, cumulative problem time** required for a single energy group, P3 scattering, 126,000 fine m esh cel ls using 12 processors on a parallel cluster.................................................................................................................................54 3-1 Computed tomography image sources fo r the developm ent of the UF pediatric phantom series................................................................................................................. ..60 3-2 Total number of directions, Cumulative Problem Time required for different energy groups, P3 scattering, and 189k fine mesh cells................................................................. 76 5-1 Comparison of selected organ absorbed dose rate (MeV/g. s) calculated using MCNP5 pulse height tally with (photon, electron mode) and EDK-SN for the UF hybrid 15-year-old male phantoms for a ch est 8 MV X-ray flat weighted source.......... 115 5-2 Prostate absorbed dose rate (MeV/g.s) calculated for the test problem using EDK-SN and MCNP5 using pulse height tally (*F8) with (photon, electron mode) at different MCNP5 Monte Carlo running times................................................................................ 116 5-2 Absorbed dose rate (MeV/g. s) calcu lated for the test problem using EDK-SN and MCNP5 using pulse height tallies (* F8) with (photon, electron mode).......................... 117 5-3 Absorbed dose rate (MeV/g. s) calculate d for the test problem using MCNP5 using pulse height tally *F8 wi th (photon, electron mode) and MCNP5 energy deposited tallies (F6)................................................................................................................... .....118 A-1 Ratios of photon interactions in water for range of energies........................................... 124 PAGE 9 9 LIST OF FIGURES Figure page 1-1 Comparison between thicknesses required to attenuate different percentages of prim ary photons and CSDAR range for maximum ener gy secondary electrons produced by the same beam............................................................................................... 331-2 Comparison between thicknesses required to attenuate different percentage of primary photons and average range for maximum energy secondary electrons produced by the same beam............................................................................................... 332-1 60Co unit (left: Mesh distribution genera ted by PENTRAN package; Right: MCNP5 Cell/Surface schematic input)............................................................................................ 382-2 60Co -rays 3-D Groups 1,2 3 and 4 scalar flux computed by PENTRAN with the CEPXS gamma cross section library................................................................................. 402-3 S12 (Level Symmetric) and S42 (Legendre-Chebychev) quadr ature first octants (with 21 and 231 directions/ oc tant respectively)....................................................................... 412-4 60Co -ray 3-D Group 1 scalar flux distribution for x-level depths at 0 cm, 5 cm, 10 cm and 15 cm, computed by PENTRAN with the CEPXS gamma cross section library........................................................................................................................ .........422-5 60Co -ray 3-D Group 1 scalar flux distribution for x-level depths at 0 cm, 5 cm, 10 cm and 15 cm, computed by PENTRAN with the CEPXS gamma cross section library; an S42 angular quadrature (1 848 directions) with P3 scattering anisotropy.......... 432-6 Scalar flux dist ributions using S12, S22, S32, and S42 with P3 scattering anisotropy in PENTRAN using the BUGLE-96 gamma library............................................................. 442-7 60Co model spatial meshing, axial level............................................................................. 472-8 Scalar flux comparison of variation in the interval of 3-D mesh cells containing air; fine meshes of 1 cm size (along each axis) in air proved to be best in providing high detail flux distribution....................................................................................................... .492-9 Scalar flux distributions generated by PENTRAN for a single energy group of each library using the MCNP5 solution as the reference case................................................... 512-10 Relative deviation of energy group 1 disc rete ordinates results with Monte Carlo results.................................................................................................................................522-11 Absorbed dose rate comparison between SN method using kerma approximation and MC techniques using F6 tally, MC statis tical uncertainty were less than 3% on average...............................................................................................................................53 PAGE 10 10 3-1 PENTRAN-MP Code System Structure Flowchart...........................................................583-2 The UF Series B pediatric phantoms series....................................................................... 603-3 Values of the mass energy-absorption coefficient, and the fitting function as a function of photon energy, for Lung (ICRU-44) soft tissue (ICRU-44) and cortical bone....................................................................................................................................643-4 UF-Series B phantom (11 years old male) with (398 242 252 voxels); the corresponding PENTRAN back-thinned mo dels with (132 80 84 voxels) and 72 materials (upper right), and another with (79 48 50 voxels) and 66 materials (lower right).......................................................................................................................663-5 Transversal view of the 11-year-old male Series B original phantomcollapsed into 189k voxels and 4 materials, collapsed into 900k voxels and 4 materials (lower left) and collapsed into 900k voxels and 72 materials............................................................... 673-6 x-ray 3-D Group 1, 3, 5 and 7 scalar flux computed by PENTRAN with the CEPXS cross section library, S32 angular quadrature (1 088 directions) with P3 scattering anisotropy...........................................................................................................................683-7 Flux distributions using S32 with P3 scattering anisotropy in PENTRAN using the CEPXS library vs. MCNP5 Monte Carlo results............................................................... 683-8 Angular flux distributions using S32 with P3 scattering anisotropy in PENTRAN using the CEPXS library vs............................................................................................... 693-9 Angular flux distributions using S32 with P3 scattering anisotropy in PENTRAN using the CEPXS library vs. MCNP5 Monte Carlo results for energy groups 1, 3, 5 and 7 along the z axis at differe nt depths in the phantom..................................................693-10 Simulation methodology for 3-D absorbed dose computations using PENTRAN-MP code system.................................................................................................................... ....713-11 CT view of the 11-year-old male Seri es B original phantom collapsed into 189 k voxels and the correspondent absorbed dose distributions 900 k voxels and the correspondent absorbed dose distributions........................................................................ 723-12 CT view of the 11-year-old male Seri es B original phantom collapsed into 189 k voxels and the correspondent absorbed dose distributions 900 k voxels and the correspondent absorbed dose distributions........................................................................ 733-13 CT view of the 11-year-old male Seri es B original phantom collapsed into 189 k voxels and the correspondent absorbed dose distributions 900 k voxels and the correspondent absorbed dose distributions........................................................................ 743-14 Simulation methodology for 3-D organ absorbed dose computations using PENTRAN-MP code system.............................................................................................74 PAGE 11 11 3-15 Dose Volume Histograms (DVH) comput ed based on PENT RAN (S32,P3) transport results using the CEPXS library vs. D VHs computed using MCNP5 Monte Carlo results.................................................................................................................................754-1 Procedure to produce EDK base d on Monte Carlo Calculations....................................... 824-2 Monte Carlo model used for comparisons between MCNP 5 pulse height tally, energy deposited tally and cell flux tally empl oying dose conversion factors (DE, DF).............. 844-3 Comparison between F6 (p, e mode), F6 (p mode), and F4 ( DE, DF) tallies in MCNP5, statistical uncertainty was less than 1%..............................................................844-4 Comparison between F6 (p, e mode) and *F 8 (p mode) tallies in MCNP5, statistical uncertainty was less than 1%............................................................................................. 854-5 Comparison between *F8 (p, e mode) and *F 8 (p mode) tallies in MCNP5, statistical uncertainty was less than 1%............................................................................................. 854-6 Deterministic calculation will provid e the angular information explicitly........................ 874-7 For each energy group, using Monte Carlo ra diation transport, a pencil beam of mono-energetic photons interacts at the center of a water phantom.................................. 884-8 Euler angles, zyz convention. Similar conventions can be obtained by following of the following order in rotation ( zyz, xyx, xzx, yzy, yxy )......................................................894-9 EDK-matrix after rotation to direct ion cosines A: (-0.5, 0.35, 0.8) and B: (0.25,0.40,-0.88) relative to the referenc e coordinate system (0, 0, 1)....................................... 914-10 First octant of S42 equal weight angles used for projecting the EDK-SN........................934-11 A flow-chart for the EDK-SN Code System...................................................................... 944-12 External Source impinging on block of water................................................................... 954-13 External Source 3-D absorbed dose di stribution for 2 Energy groups computed by EDK-SN based on PENTRAN-MP A. EG(7.5-8.0Mev) and B. EG(0.5-1.0Mev)............ 974-14 Energy group absorbed dose rate distri butions (in Mev/ g.sec) along the X-axis using the EDK-SN method based on PENTRAN-MP....................................................... 974-15 Energy group absorbed dose rate distri butions (Mev/ g.sec) along the y-axis using the EDK-SN method based on PENTRAN-MP................................................................ 984-16 EDK-SN Total absorbed dose rate distri butions along the X-axis compared to MCNP *F8 tally derived absorbed doses.......................................................................................984-17 Total absorbed dose rate distributions of EDK-SN along the Y-axis compared to MCNP *F8 tally................................................................................................................. 99 PAGE 12 12 5-1 Three Monte Carlo phantoms used in the absorbed dose study....................................... 1045-2 Monte Carlo derived total absorbed do se rate distributions along the X-axis computed using *F8 MCNP5 tally for three different phantoms.................................... 1055-3 Monte Carlo derived total absorbed do se rate distributions along the X-axis computed using *F8 MCNP5 tally for three different phantoms.................................... 1065-4 Monte Carlo derived total absorbed do se rate distributions along the X-axis computed using *F8 MCNP5 tally for the _SBLLBSS_ phantom.................................. 1065-5 Three SN phantoms equivalent to the one us ed in the Monte Carlo simulation..............1075-6 EDK-SN heterogeneous SBLLBSS phantom to tal absorbed dose rate distributions along the X-axis compared to MCNP *F8 tally............................................................... 1085-7 EDK-SN Homogenous phantom total abso rbed dose rate dist ributions along the Xaxis compared to MCNP *F8 tally................................................................................... 1095-8 Continuous slowing down approximation ra nge (Rcsda) for electrons in Cortical Bone, Soft Tissue and Lung as a function of energy....................................................... 1105-9 Continuous slowing down approximation range (RCSDA) ratio comparison for electrons in Cortical B one, Soft Tissue and Lung as a function of energy...................... 1115-10 Frontal views of 3-D rendering of th e 50th percentile UFH-NURBS 15-year male and female........................................................................................................................1125-11 Simulation methodology for EDK-SN computations using PENTRAN-MP code system..............................................................................................................................1135-12 Absorbed dose contribution of each en ergy group is calculated via the EDK-SN procedure, then folded to calculate the to tal absorbed dose deliv ered to the 15 year old male phantom from 8 Me V flat weighted source......................................................115A-1 Total mass attenuation coefficient for photons in water, indicating the contributions associated with phoelectric absorption, Compton scattering and electron-positron pair production.................................................................................................................124A-2 Sketch of Compton scattering, photon of energy h collide with an electron. The photons incident forward momentum is hc where c is the speed of light in vacuum.............................................................................................................................125A-3 Compton effects differential cross section per unite angle for the number of electrons scattered in the direction e..............................................................................................128A-4 Pair Production in the Coulom b force field of the nucleus.............................................. 130 PAGE 13 13 A-5 Pair Production in the Coulom b force field of the nucleus .............................................. 132 PAGE 14 14 Abstract of Dissertation Pres ented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy 3D DETERMINISTIC RADIATION TRANSPORT FOR DOSE COMPUTATIONS IN CLINICAL PROCEDURES By Ahmad Al-Basheer August 2008 Chair: Glenn Sjoden Major: Nuclear Engineering Sciences The main goal of this dissertation was to es tablish the feasibility of basing megavoltage external photon beam abso rbed dose calculations in voxelized phantoms on SN deterministic calculations and pre-calculated electron absorb ed dose kernels derived from full-physics Monte Carlo. The SN derived electron absorbed dose kernel method EDK-SN, developed as part of this research, achieves total execution times that ar e on the order of several times to orders of magnitude faster than conventional full-phys ics Monte Carlo electron transport methods considering equivalently detailed models and data fidelity. With the rapid movement toward intensity mo dulated radiation therapy (IMRT), radiation beam intensities have increased dramatically over the past decade, thus heightening the need for further characterization of out-o f-field organ absorbed doses, along with their associated biological risks. Assessment of these tissue absorbed doses is complicated by two fundamental limitations. First, anatomic information on the patie nt is generally restrict ed to a partial body CT image acquired for treatment planning; conseque ntly, whole-body computational phantoms must be employed to provide the out-of-field anatomy model structure for absorbed dose evaluation. Second, existing methods based on Monte Carlo radiation transport, even with the application significant variance reduction, are quite computati onally inefficient at large distances from the PAGE 15 15 primary beam, and point-kernel methods do not properly handle tissue inhomogeneities. Moreover, since absorbed dose ar e generally tracked in all majo r organs in th e body, variance reduction schemes for Monte Carlo are not all effective in this regard. The outcome of this dissertation is to dem onstrate that absorbed dose from high-energy external beams radiation can be accurately computed for whole body and organ-specific absorbed doses. The EDKSN method implements voxelized phantoms with discrete ordinates (SN) transport computations coupled with directional influences and (pre-computed) full-physics Monte Carlo based electron absorbed dose kernel s to yield total body abso rbed dose information. This research shows that the deterministic t echniques coupled with Monte Carlo based electron absorbed dose-kernels has significant potential fo r organ absorbed dose eval uation in the clinical management of radiation therapy patients. PAGE 16 16 CHAPTER 1 INTRODUCTION AND LITERATURE REVIEW 1.1 Overview W ith the growing implementation of radiat ion applied in differe nt therapeutic and diagnostic medical procedures, comprehensive and accurate ab sorbed dose computations are fundamental to precisely understand the impact and outcome of such procedures. In practice, the clinical outcome in radiothera py can be significantly improved w ith enhancements provided by a treatment planning systems (TPS)1. Current TPS have been based on augmenting clinically measured absorbed dose distributions. Such me thods yield only approximate solutions, and are limited in rendering a precise solution. Detaile d solutions are achiev able by solving the Boltzmann transport equation, which describe the beha vior of neutral partic le radiation transport throughout the phase space of a problem. Often absorbed dose calculation systems employ relatively simplistic yet fast algorithms to pred ict the patient absorbed dose distribution from first principles, using simp lified radiation transport models2. Overall, these techniques are not capable of accurately describing high fidelity whole body detailed radiation transport through the human body, particularly as new procedur es in medical physics are developed. In recent years, Monte Carlo techniques have proven to be ca pable of representing very complex geometries in a rigor ous manner. This gave the Monte Carlo technique a great advantage over previous methods in rendering accurate estimates of tissue absorbed doses, especially when simulating human anatomy and yielded accurate prediction of absorbed dose delivered in various medical physics procedures. However, this accuracy is achieved at a high cost: extremely long running times are typically necessary to achieve a statistically precise solution. This limits the practical use of Monte Carlo simulations for clinical procedures even with High Performance/Parallel Computing (HPC) 3. PAGE 17 17 Solving the Boltzmann equation directly using deterministic techniques can be significantly more efficient than using the Monte Carlo methods, a nd therefore it can be used to improve the quality of radiotherapy planning. Re lative to the amount of information generated, and with the use of HPC resources, such tech niques usually require significantly lower running times (relatively to Monte Carlo) to provide deta iled and precise solutions. Moreover, the rapid proliferation of Computed Tomography (CT) and Magnetic Resonance Imaging (MRI) enable human anatomy to be represented using regular images consisting of rectangularly (hexahedral) shaped voxels in a regular matrix. These new developments make Cartesian geometry deterministic techniques ideal for medical physics applications, especial ly given the demand for high resolution absorbed dose distributions. This research is focused on the development of EDK-SN, a first-of-a-kind capability to assess organ absorbed doses incurred in radiati on therapy (RT) applica tions using anatomical patient models coupled with 3-D deterministic discrete ordinates (SN) radiation transport methods. Absorbed doses delivered as a result of secondary charged particles will be linked using electron kernels to account for secondary electron effects. This methodology is proposed for regions in the phase space where Charged Particle Equilibrium (CPE) is not satisfied, and subsequent absorbed dose estimation us ing a collisional kerma is not valid. 1.2 Initial Clinical Dose Evaluations In the field of radiation therapy, absorbed dose distributions were historically first approxim ated without doing any explicit radiation tr ansport at all, but rather by extrapolating from measured results taken for a particular en ergy, field size, and sour ce-to-surface distance (SSD) on a given machine4. An extrapolation was performed by modifying the absorbed dose that was experimentally determined under the re ference parameters by ratios representing the change between specific paramete rs and the parameters of interest. For example, if absorbed PAGE 18 18 dose values were measured in water, a table of the percentage of the maximum dose as a function of depth (on-axis) in the water could be deduced. A similar table could be generated for off-axis ratios (OARS) by measuring latera l profiles at a variet y of depths. The absorbed dose at a given point P in the patient could then be estimated by knowledge of the maximum dose for the treatment, and applying a percent depth dose (PDD) correction for the depth of P in the patient, and an off-axis correction for the distance to P from the axis at that depth4. Tables of PDDs were made available for vari ous sources (incident en ergies), field sizes, depths in the phantom, and SSDs. To modify th e dose for an SSD that was different from the reference SSD, the ratio of inverse distances s quared would be applied. For a change in field size, tables of field-size-factor s would be subsequently applied. To compensate for the curved surfaces of a patient, as opposed to the flat surf ace of the phantom, isodose lines could be shifted along rays projected from the source4. It is important to point out that these techniques are often quite successful for predicting absorbed dose in a homogeneous volume. In f act, in some cases they may be superior to alternative simplified methods (for homogene ous volumes) because they are based on measurements from a particular machine; thus, th ey circumvent the proble m of how to model the photon distributions emerging from the source and in cident on the patient/phantom, which can be a significant obstacle in accomplishing comput ational radiation transport based models. The real challenge for these empirical-bas ed methods, then, lies in correcting the homogeneous water absorbed dos e distribution to account for th e many inhomogeneities present in the human body, including air cavities, bone, de nse muscle mass, and sparse lung tissue. One way of addressing this problem is to a pply equivalent-thickness corrections to the water measurements. These corrections were ba sed on compressing or stre tching isodose lines by PAGE 19 19 factors that were derived from the relative phot on/electron density in the path through tissue as compared to that in water4,5. These methods, at their best, were useful for obtaining depth dose information in the center of a broad beam, for a medium consisting of homogeneous layers. Even then, these approaches do not a ccount for the fact that photons in teract not only wi th electrons, but also with the nucleus. For example, the elasti c scattering cross section depends on the atomic number of the material. The clinical implication of this is that not all materials behave like scaled-density water; real bone, for example, s catters photons to greater extent than water of equivalent bone density due to the presence of Ca and other high Z elements not found in liquid water. The effect of this was bette r addressed in the Pencil Beam method6. 1.3 Pencil Beam Transport. The Pencil Beam method sums the dose distri bution from individual small diameter rays called pencil beams. A pencil beam is made up of particles which pass through a differential cross sectional area, x y The off-axis dependence of the dose distribution for each pencil beam is described by thick-target multiple scattering7. The on-axis dose information is obtained from measured depth-dose data. The calculation of absorbed dose assumes that the patient or phantom can be represented by a stack of slabs of different material types. E ach individual slab is homogeneous and infinite in lateral extent. Assuming a single pencil be am has normal incidence on such geometry, the absorbed dose at a position (x, y, z) can be expressed as (,,)(,,)() Dxyzfxyzgz (1-1) where the first factor is the off-axis term, and th e second factor is the ce ntral-axis term. The offaxis term may be represented by a Gaussian of the form PAGE 20 20 22 221() (,,)exp()zz x y fxyz (1-2) where z is the angular spread about the axis incl uding the penumbra. The on-axis term may be expressed in terms of its relations hip to the depth dose in matter, 0()effgz 2 0()()eff effSSDz gzgz SSDz (1-3) where SSD is the source-to-surface distance, and effz is the effective depth in water, assuming the linear stopping power of the material at depth z is relatively independe nt of electron energy for normal body tissue. The quantity 0()effgz can be related directly to a depth dose curve measured in a water phantom, assuming a uniform incident beam, by the Eq. 1-4 '2'2 '' 0 0 22() ()() (,,)expzz xygz xxyy Dxyz dxdy (1-4) where the limits of the integral depend on the dispersion of the field at depth z. An example of pencil method used for clinical photon radiotherapy is an application of the Ahnesjo pencil beam model6. In this method, combination of measured data and Monte Carlo calculated convolution kernels are used to char acterize clinical photon beams and eventually calculate the absorbed dose. The strength of the Ahnesjo pencil be am method is that it is fast enough for clinical use, on the order of minutes. It also uses measured data for input that may be specific to the treatment unit. However, because this method applies simple corrections using pencil beams in homogenous matter, it is not accurate as a comple te coupled charged-particlephoton Monte Carlo simulation in he terogeneous media. Moreover, it does not fully represent the directional component of the scattered radiation. PAGE 21 21 Aside from the geometric assumptions of in finite homogeneous slabs, the pencil beam model does not handle changes in scatter from late ral heterogeneities; devi ation of up to 5% can be expected for low energy photon beam s, and 14% for high energy photon beams8,9. 1.4 Solving the Boltzmann Equation At the m ost fundamental level, radiation transport simulation involves determining the distribution of the particles (in physical space, direction, time, and energy) after they have traversed some distance into a medium. Formal ly, a way to begin this discussion is by writing the general equation for particle transport from kinetic theory of gases. This theory does not explicitly treat individu al particles, but rather describes th e aggregate behavior of a continuous particle distribution in a statistical manner. Let n be the expected number of particles within a differential volume element 3drat position r with a velocity (which is directly related to its kinetic energy) that is within a differential velocity3dv, at time t. Then, if a general solution is found for the particle density function n(r, v, t) in phase space, then the transport problem is solved. When transporting the radiati on through an arbitrary volume V, bounded with a surface area S, the rate of change of number of particle s is due to flow out through the surface of V, collision events within V, or sources (scattering and independe nt sources) of particles located inside of V. (all the terms of rate of change; commonly the first two terms are losses and the last term is production) R ateofChangedueChangedueChangedue changeintoflowtocollisionstosources nattimetthroughSinVinV PAGE 22 22 This balance of neutral particles is cast as the linear Boltzmann Transport Equation (BTE). The BTE is an integro-differential equation which describes the behavior of neutral particles in terms of spatial, angular, and energy domains as they interact in a system; the steady-state transport equation is given in Equation (1-5): (,,)(,)(,,) rErErE 04 (,',') (,',')(,,)sdEdrEErEqrE (1-5) where (,,) rEdrdEd is the angular flux with energy between E and E+dE at position r within the volume element dr, and in direction within the solid angle d. The scattering term, defined with the double integr al term, represents the sum of the partic les scattered into dEd from all the '' dEd after scattering collisions repr esented by the double differential cross-section (,',')srEE (,,) qrE represent the external source. The left side of Equation (1-5) represents streaming and collision terms (loss), and the right side represents scattering and independent sources (gain). Since it describes the flow of radiation in a 3-D geometry with angular and energy depe ndence, this is one of the most challenging equations to solve in terms of complexity and model size. For rende ring a deterministic computational solution for a large problem in a reasonable time, it is necessary to utilize a robust parallel transport algorithm and a high perf ormance computing system. All the methods described below are different approaches to solvin g this equation. For cases with charged particle transport, the Boltzmann and Fokker-P lank formalisms (Boltzmann-FP) 10are generally used. A simplified form of this equation is usually used, which is: PAGE 23 23 (,,)(,)(,,) rErErE 04 (,',') (,',') [(,',)(,',')] (,,)sdEdrEErE SrErE qrE E (1-6) However, the full Boltzmann-FP equation is extremely challenging to solve deterministically. 1.4.1 Monte Carlo Method Monte Carlo techniques are a widely used class of computational algorithms for the purpose of statistically modeling various physical and mathem atical phenomena. In general, individual probabilistic events are simulated sequentially, where probability distributions governing particle interactions ar e statistically sampled to desc ribe the total phenomenon. In principle, for each case sampled, there are four required inputs; namely, a random number generator, source characteristics, a simulatio n phase space, and probability distributions governing the possible outcomes. Results of indi vidual particle historie s are generated using pseudo-random numbers to sample the probability distributions. This process is repeated many times, until the statistical precision required for the solution is achieved, based on the law of Large Numbers and the Central Limit Theorem11. Therefore, the Monte Carlo technique is an effective method for indirectly solving the Boltzm ann radiation transport equation via statistical methods, and can readily treat th e full physics of coupled photon-el ectron transport. Simulations for particle histories are repeat ed many times for large number of histories, where the computer code tracks the desired physical qu antities, means, and uncertainti es (variances) for all of the simulated histories. PAGE 24 24 One of the most important issues limiting th e use of Monte Carlo techniques for clinical absorbed dose calculations is the long execution time required to reduce th e statistical variance associated with the calculated absorbed dose at each differential volume. At the same time, the outcome of any radiotherapy procedure depends solely on the accuracy of the radiotherapy treatment planning system, so that a small cha nge in the absorbed dose delivered ( 5%) can result in a dramatic change in local response of the tissue ( 20%)12,13. In order to reduce the statistical variance 1 N associated with Monte Carlo, large numbers of histories N are required, and consequently this l eads to higher running time to yield acceptable results. Improved variance reduction techniques may reduce the running time. However, such techniques are to be used with caution, as improper application of statistical variance reduction, method applied un-conserva tively may cause wrong answers. Also, variance reduction is often tied to one objective; this can be problematic when faced with determining a dose at all locations in the phase space. The increasing use of voxelized human body ge ometries has created challenges in many aspects of applying Monte Carlo simulations14. With increasing use of modern technologies such as CT and MRI in vivo imaging, anatomical de finition has shifted from using simplified surface equations to segmented tomographic images in order to truly reflect the anatomy of a specific human body. The segmented images are represented in the form of a large uniform 3-D matrix of hexahedral voxels, where each vox el represents a constant mate rial with a tag number (for example, CT number). One digital image ma y contain billions of voxels, making it very challenging be fully represented explicitly using Monte Carlo, particularly considering large models size and number of tally sites, addi ng significant overhead to the computation. Some special data management treatments are ne eded for Monte Carlo simulations to make PAGE 25 25 computation times more tractable, includi ng voxel down-sampling and parallel computing resources. 1.4.2 Deterministic Methods Several different methods have been deve loped for numerically solving the Boltzmann equation11,15. In general these methods depend on disc retizing the energy, angular, and spatial domains of the Boltzmann transport equation. 1.4.2.1 Discrete ordinates method (SN method) The most common deterministic approach cu rrently used for photon transport is the discrete ordinates method (SN method) developed in detail by Carlson over 50 years ago16. The discrete ordinates equations are derived by discretizing the BTE over phase space elements (i.e., spatial-energy-direction variables) to construct discrete-variable cell par ticle balances. Equation (1-7) is the standard Legendre expanded Cartes ian multigroup integro-differential formulation of the BTE, and z y x ,, are the spatial variables along each ax is in the Cartesia n coordinates. G g L l l k k l lgl lggs g g gP kl kl zyxPzyx l zyxzyx zyx zyx1 '0 1 ,' ,',! 2,, ,,12 ,,,,,,,,,, ', ',,,cos ,,sin ,,,,kk CglSglgxyzkxyzkqxyz (1-7) Other variables in Equation (1-7) are: correspond to the direction cosines along each axis for each ordinate g is the energy group angular flux g is the total group macroscopic cross section lggs ,', is the l th moment macroscopic scattering cross section for group g to g PAGE 26 26 lg ,' is the l th flux moment for group g k lCg,' and k lSg,' are l th, kth associated cosine, sine flux moments for group g ) ( lP and )(k lP are the Legendre and Associ ated Legendre polynomials. ,,,,gqxyz is independent source for group g The SN methodology is a widely accepted determin istic scheme in yielding numerical solution to neutral particle radiation transport problems, and can be applied to medical physics problems, provided that discretizat ion and truncation errors can be reduced to an acceptable level over the problem space. Depending on the amount of data needed per floating point operation, the SN method is more efficient than the Mont e Carlo method for many problems, and with proper discretization produ ces an accurate solution17. The SN method is widely used in nuclear engineering to obtain a numerical solution of the Boltzmann transport equation. The method solves the BTE for a discrete number of ordinates (directi ons) with associated weights16. The combination of discrete directions and weights is referred to as quadrature set11. A drawback of the SN method is that with the limited number of directions involved, in certain situations, this may lead to so called ray-effects, which appear as unphysical flux oscillations. In general, this behavior occurs for problems with highly angular dependent/localized sources and/or fluxes, or when the source is localized in a small region of space, in low density or highly absorbent media. This can make medical physics applications, whic h often require small mesh sizes, relatively large models, and low density materials, suffer severe ray-effects. The ray-effects can be alleviated naturally by the presence of a distri buted isotropic source, which tends to flatten the flux in the angular domain. However, in medical phys ics, the photon transport problems may exhibit significant ray e ffects. This effect becomes worsen in a low PAGE 27 27 scattering medium such as air, or in a strongl y absorbing medium such lead, where the photon flux usually is typically strongly angular dependent. In such situations, either a higher order quadrature set (many directions) is required or some ray-effects remedies must be applied11,18 (such as Taylor Projection techniques and adaptive differencing schemes19). More discussion on discretization and ra y-effects are presented in Chapter 2. 1.4.2.2 Monte carlo versus deterministic SN methods The most fundamental difference between the Monte Carlo methods and the deterministic SN transport method is that the SN method solves the transport equa tion for the average particle behavior throughout the problem phase space, while the Monte Carlo method, the terms of the transport equation are solved statistically by si mulating the individual pa rticles and recording some aspects (tallies) of their average behavior but only at specific tall y locations requested by the user. In the discrete ordinates method, the phase space is divided into many small meshes, in which the cross section is assumed to be consta nt. In the medical physics applications, the model geometry, in this case, the pati ent is represented using regular voxel-based grid generated from a segmented CT scan data or other digital imaging technique. The types of grids can then be easily represented in an input deck for the SN calculations, making both Monte Carlo methods and deterministic techniques capable of representing the model of study accurately. Another advantage of the MC methods is that they can ea sily utilize continuous energy interaction data, whereas deterministic calcul ations use multigroup or energy group averaged data. However, photon attenuation coefficients in the energy range of interest for radiotherapy (RT) are quite smooth with little structure (in contrast to neutron cross sections); and are accurately represented by multigroup da ta. Recent research we conducted20 has shown that multi-group deterministic methods produce highly accurate absorbed dose distributions in photon PAGE 28 28 transport calculations. Similar conclusions were drawn by Daskalov et. al and Williams et. al2123. Table (1-1) shows a comparison between Monte Carlo and Deterministic techniques24. Table 1-1. Comparison between Monte Ca rlo methods and deterministic techniques 1.5 Dose and Charged Particle Equilibrium Absorbed dose, D, can be defined as a non-st atistical quantity in te rms of a statistical quantity, energy imparted The absorbed dose to matter at a point P is given by / Dddm (1-8) where is the expectation value of the energy imparted by ioni zing radiation to matter of mass m in a finite volume V and inout R RQ (1-9) in which in R = the radiant energy incident on the volume V out R = the radiant energy emerging from V and Q= the net energy derived from rest mass within V (negative for the reverse process). where the radiant energy R can be defined as the energy of particles (excluding rest energy) emitted, transferred, or received25. PAGE 29 29 Collisional kerma,cK can be defined in terms of th e related net energy transferred, n tr 26,27 and R25 as. /n trddm (1-10) The energy-transferred in a volume V is ()()n nonrr trinuoutuu R RRQ = r tru R (1-11) in which ()inuR = the radiant energy of unchar ged particles entering the volume V ()nonr outuR= the radiant energy of uncharged particles leaving V except that which originated from radiative losses of kinetic energy by charged particles while in V r u R = radiant energy emitted, where u indicates uncharged particles. Charge particle equilibrium (CPE) can be es tablished in a volume when the energy carried into the volume V by charged particles is balanced by the energy carried out of the volume V by the charged particles. Note that, CPE is not established for thin materials irradiated by photons. The range of secondary electrons produced is larger than the material dimensions, so that much of the electron energy is lost simply by electrons leaking out of the material. Under the conditions of CPE, Attix demonstrat ed that when the energy imparted equals energy transferred, tr collisional kerma cK is equal to the total dose D26,27 Eq(11). ,.CPE en c EZDK (1-12) PAGE 30 30 where is the energy fluence, and ,en EZ is the mass energy absorption coefficient for ray energy E in medium Z. This result is very im portant, since it allows us to calculate the measurable dose D based on the calc ulable quantity cK 1.5.1 Deterministic Dose Calculatio ns for Diagnostic Application For low energy beams, such as the one used in diagnostic applications, the concept of charged particle equilibrium is of a great im portance as a mean of equating absorbed dose D to collisional kerma cK For a volume V, containing a smaller volume v, irradiated with an indirectly ionizing radiation field (i.e., Photons and Neutron fields), there are four basic causes for CPE failure in volume v: In-homogeneity of the atomic composition within volume V In-homogeneity of density in V Non-uniformity of the field on indi rectly ionizing radiation in V Presence of a nonhomogeneous el ectric or magnetic field in V These four situations will disrupt CPE by prev enting charged particles of a given type and energy leaving the volume v to be replaced by an identical particle of the same energy entering27. In practice, CPE conditions are met if the attenu ation of the primary beam is negligible in the maximum range of the secondary electrons ge nerated by the beam. As indicated by Table (12) for a layer of water of 5 cm thickness, this scenario is plausible in the limit of diagnostic medical physics applications, where the energies fall within 10 to 200 keV, and CPE conditions are met for water/tissue volumes less than 55 cm3, and consequently D is equal to cK under such conditions. PAGE 31 31 Table 1-2. Approximate attenua tion of photons with a layer of water equal to the maximum range of secondary charged particles28 1.5.2 High Energy Beam Dose Calculations Based on the previous discussion, as the ener gy of the ionizing radiation increases, the penetration power of the secondary charged particles increases more rapidly than the penetration power of the primary radiation in a given volume. This leads to CPE failure, caused by the small number of charged particles produced at deeper depths in volume V relative to the number of charged particle produced at initial points along the primary radiation field direction27. This is indicative of significant atte nuation of the primary indi rectly ionizing radiation. As a result of this phenomenon, higher numbers of charged particles will be generated close to the entrance compared to the rest of the volume, leadi ng to CPE failure considering dose in that volume. The degree of this failure becomes progressively large for hi gher energies, so that D no longer equals tocK A more detailed discussion on this topic is presented in chapter 4. 1.5.3 Limits of Charged Particle Equilibrium Assuming that a maximum energy is transferred from a photon to an electron, we can use this assumption to investigate a CPE requirement for a water phantom region, as defined earlier. Considering the continuous slowing dow n approximation range of the electron (CSDAR), and average range of the maximum energy electrons PAGE 32 32 0 0 0 01 .() ()ftdt ttttdt N dNt dt dt (1-14) where dT dx is the electron stopping power (MeV cm2/g) T0 is the starting energy of the electron (MeV) () ()fdNt tt dt is the differential distribution of farthest depth of penetration ()Ntis the number of particles penetrating a slab of thickness t Through inspection of these relationships de scribed in Figures 1-1 and 1-2, one can decide what kind of tolerance is acceptable in applying a collisional kerma approximation for a dose calculation based on a particular voxel size i.e. to determine CPE approximation limits. This also has been demonstrated in Table 13; for example at 2 MeV the CPE computed dose will be accurate to approximately 5% using a 1 cm voxel, assuming the CSDA range applies. Table 1-3. Approximate thickness of water needed to achieve various attenuation amounts compared to maximum energy electron ranges generated by the same beam PAGE 33 33 Figure 1-1. Comparison between thicknesses required to attenua te different percentages of primary photons and CSDAR range for maximum ener gy secondary electrons produced by the same beam Figure 1-2. Comparison between thicknesses required to attenua te different percentage of primary photons and average range for maximum energy secondary electrons produced by the same beam With proper model discretization of the angleenergy-spatial transpor t phase space, and at low photon energies yielding primary and secondary electron interactions with small electron PAGE 34 34 transport paths relative to a spatially discretized mesh grid interval, 3-D SN methods can be implemented to directly yield ve ry accurate dose distributions. Th erefore, the deterministically derived dose using a kerma approximation will be accurate up to the range of electron path lengths comparable to ranges smaller than the photon spatial SN mesh grid interval in anatomical models (typically representative of the anatomical data voxel size). As mentioned earlier, at higher photon ener gies, where electron ranges exceed the SN mesh grid interval, this is not the case, and interacti ons yielding energetic el ectrons must account for the mechanism of electron tr ansport through the problem phase space to deliver an absorbed dose after representative electron transport dist ances. This work will set the foundation for state of the art 3-D absorbed dose computations for high energy external beams to account for interactions yielding ener getic electrons. The EDK-SN methodology will account for electron transport in the phantom by applying a pre-comput ed electron absorbed dose kernel calculated to a very low variance It is worth comparing the EDK-SN methodology, proposed here with the traditional convolutionsuperposition techni que, which applies a database of energy deposition kernels, calculated using Monte Carlo techni ques to determine a absorbed dose29,30. We note that these kernels can be used to describe the spread of energy about the primary photons point of interaction in the phantom, however, these have limited use since the co nvolutionsuperposition technique is performed without regard to the incident photon angular component. While, in contrast, the EDK-SN procedure is based on a detailed, full physics Boltzmann transport solution with direction and flux information directly available for coupling with the EDKs along a propagated photon dire ction of travel. PAGE 35 35 1.6 Outline of This Research The method developed in this dissertation, EDK-SN, includes implementing voxelized phantoms with discrete ordinates (SN) transport computations coupled with multi-directional and (pre-computed) electron absorbed dose kernels (EDK). Particular care was taken when generating the EDKs using full-physics Monte Carlo calculations (with the MCNP5 code) and tallying the results using pulse he ight tally (*F8) tallies in MCNP5 to extremely low statistical uncertainties. The application of the pre-com puted EDKs was based on the photon flux and net current vector for the dose driving voxel (DDV) calculated via the disc rete ordinates parallel code PENTRAN with a voxelized phantom. Chapter 2 introduces code systems, and presents is a study of critical discretization issues in 3-D discrete ordina tes medical physics simulations benchmarked with the Monte Carlo. The discre tization of the phase space variables in the transport equation will be disc ussed, along with the absorbed dose calculation comparisons. Chapter 3 discusses whole body dosimetry simulations with the UF-B series phantoms. It also presents the creation of the PENTRAN-MP code system. Chapter 4 elaborates further on CPE, and compares absorbed dose calculations using different absorbed dose tallies in the MCNP5 code. It also introduces the EDK-SN methodology and the parallel st rategies. In Chapter 5, we apply the new EDK-SN methodology to heterogeneous phantom s, including the whole body 15 year old phantom and compare the results to the full physics (Monte Carlo MCNP5 *F8) tally. Chapter 6 summarizes the research and draws conclusions on the objectives accomplished; aspects for future development are also highlighted. PAGE 36 36 CHAPTER 2 CRITICAL DISCRETIZATION ISSUES IN 3-D DISCRETE ORDINATES MEDICAL PHYSICS SIMULATION Discretization of the BTE m ust be well unders tood in solving this equation for medical physics problems. As discussed, the fundamental difference between the SN and Monte Carlo methods is that the SN method directly solves the BTE an d requires phase space discretization, while the Monte Carlo method applies a statistical approach to represent the radiation transport. It is important to explore discreti zation issues associated with the SN method. To account for discretization effects, 3-D flux distributions were computed in a water phantom exposed to 60Co -rays using the PENTRAN parallel SN code and compared to analogous models implemented using the continuous energy MCNP531 Monte Carlo results with E NDF/B-VI data libraries. A major issue affecting the accuracy of the SN method is selecting an ap propriate discrete set of directions, or the angular quadrature set for which the transport equation is to be solved. Choosing a quadrature set that does not include enough discrete directions to represent the physics can lead to ray-effects. The ray-effects are numerical artif acts sometimes exhibit in the SN computations when there is high radiation streaming and little scattering. With a finite number of directions considered in each fine mesh cell for the radiation emanating at points distal from the source; it often re sults in the spines at points distal from the source. In this work, a number of studies have been conduc ted through comparison with the Monte Carlo simulations to determine suitable quadrature sets required to reasonably li mit the ray-effects in a proposed radiotherapy model that includes radiation transport in ai r to still yield high accuracy. Another set of simulations were accomplished to address anothe r difficulty in performing the SN simulations with regard to spatial discretization. Selection of an adequate 3-D spatial mesh size to represent the different geometric features of the model is also criti cal; at the same time, selection of a differencing scheme capable of representing the physics w ithin tractable memory PAGE 37 37 storage and time limitations is also a must. This i ssue has been investigated particularly in the low density portions of the model (mostly in air). Finally, we compare the performance achieved with three cross section libraries availabl e from the nuclear engineering and physics communities. These multi-group gamma-ray libraries were derived from the BUGLE9632library, the CEPXS package33, and the SCALE5 package. 2.1 Code Systems PENTRAN (Parallel Environment Neutral particle Transport) is a 3-D discrete ordinates code that solves the Boltzmann transport equation on distributed parallel co mputers. It is capable of complete, hybrid phase space decomposition with up to a 97% parallel fraction19,34. In executing a parallel transport calculation on each distributed node, PENTRAN performs discrete ordinates sweeps only on the local portion of th e phase space assigned. As problem sizes increase, more machines may be needed to s cale up to accommodate the total parallel memory demand. The inherent parallel scalability of PENTRAN with advanced numerical treatments (including block adaptive mesh structures and ad aptive differencing schemes) makes it ideal to apply to challenging problems in medical physics. The PENTRAN c ode is part of a code system that includes preand post-proce ssors, and other utilities to manage large parallel datasets. For the purpose of comparison to the M onte Carlo methods, the well known MCNP5 Monte Carlo code developed by Los Alamos National Laboratory is used to simulate a 60Co beam for an equivalent unit. MC NP5 is used to calcula te the scalar flux along the central axis of the model. For these studies, cell flux mesh-ta lly (F4) and energy deposit ion volume tallies (F6) are defined in the Monte Carlo model in a way that is analogous to the Cartesian fine mesh cells used in the PENTRAN. PAGE 38 38 2.2 Simulation Geometry Because the purpose of our investigations was to identify critical discretization issues, the models we considered were simplified to enable us to focus on proper numerical discretization. 2.2.1 PENTRAN Model We considered a box of lead containing a 60Co source and a water tank for use in multiple benchmarking applications (Figure 2-1). As we mentioned earlier, this model consisted of four key components; a 60Co source (2 cm3), a lead source housing 4 cm thick, an air layer, and a homogeneous water phantom (30 cm3) at a source-to-surface distance (SSD) of 34 cm. A 3-D spatial distribution was generated by the PENMSHXP35 mesh generating code in the PENTRAN system. The model was divided into 5 z-levels at 0 cm, 30 cm, 60 cm, 64 cm, and 66 cm; each z-level was divided equally into 9 coa rse meshes containing various numbers of fine spatial mesh cells; a total of 126,000 fine mesh cel ls were used in the SN problem. Figure 2-1. 60Co unit (left: Mesh distribution genera ted by PENTRAN package; Right: MCNP5 Cell/Surface schematic input) PAGE 39 39 2.2.2 Monte Carlo Model The MCNP5 Monte Carlo model geometry was define d to be equivalent to that used in the PENTRAN model. The 60Co source was distributed isotropically over multiple equi-probable sources cells in one energy group. Source spectra were defined in a consistent manner with the various multi-group libraries cons idered. Volumetric cell flux F4 mesh-tallies and energy deposited F6 tallies were used in the Monte Carl o geometry description that was equivalent to the discretized SN volumes. The MCNP5 F6 tallies were used along the central axis of the model across the lead box, air region, and water phantom. All simulations were performed using a PC cluster composed of 16 Dual Intel Xeon processors at 2.4 GHz with 2 Gb/processor, for a total memory of 32 Gb. MCNP5 tallies along the centra l axis were obtained using variance reduction techniques. The cell fluxes c onverged, on average, to within 3 % for the reference case. 2.3 Comparison of Deterministic and Monte Carlo Results Calculations were performed to compare the SN and Monte Carlo results for the scalar fluxes along the model central axis. Figure 2-2 de picts a 3-D scalar flux distribution computed by PENTRAN for the reference case for different energy groups. The highest energy group (group 1) is forward peaked with some ray-effects al ong the central axis. However, as the radiation down scattered to lower energy groups, it demonstrates more isotropic behavior evident in energy groups 2 and 3. This simulation was based on the use of an S42 angular LegendreChebychev quadrature (1848 directions/mesh) with P3 scattering anisotropy using the CEPXS library. This highlights the fidelity of the information available in a converged deterministic SN computation. 2.3.1 Investigation of Angular Quadrature Medical physics applica tions represented in SN models often require small mesh sizes, relatively large models, and low density materials. These problems can lead to severe ray-effects, PAGE 40 40 which appear as unphysical flux oscillations in the scalar flux distribution unless treated properly. In general, this be havior occurs for problems with highly angular dependent sources and/or fluxes, or when the source/sink is localiz ed in a small region of space, in low density or highly absorbent media. Group 1 Group 2 Group 3 Group 4 Figure 2-2. 60Co -rays 3-D Groups 1,2 3 and 4 scalar flux computed by PENTRAN with the CEPXS gamma cross section library, (S42 angular quadrature (1848 directions) with P3 scattering anisotropy where used) To mitigate this problem, high angular quadrature orders are often needed. Figure 2-3 depicts a schematic of the ordinate de nsity distribution associated with S12 level-symmetric (rotationally level-symmetric) and S42 Legendre-Chebychev quadrature sets on 1/8th of a unit sphere (other octants are symmetri c). Note that ordinates (vectors ) in Figure 2-3 would begin at the origin and project to e nd at the square point plot ted in the Figure. Beyond S20, LegendreChebychev quadrature sets must be used due to negative weights that occur in standard levelsymmetric sets36. Particular attention must be made in deriving quadrature sets to enable preservation of Spherical harmonics (Legendre) scattering moments a nd numerical quadrature orthogonality relationships. PAGE 41 41 S12 S42 Figure 2-3. S12 (Level Symmetric) and S42 (Legendre-Chebychev) quadr ature first octants (with 21 and 231 directions/ octant respectively) Using a uniform 1 cm mesh (for reasons di scussed later) in the air region, PENTRAN was used to calculate the gamma flux di stributions in the system for several different quadrature sets. PENTRAN provides the option of r eadily using any quadrature that could be prescribed by the user. The flux distributions produced by an initia l computation, using S12 level-symmetric quadrature with P3 scattering anisotropy and 126,000 fine mesh cells, is shown in Figure 2-4. The ray-effects in the S12 computation are evident when compar ed to scalar fluxes in Figure 2-5 using S42 quadrature, with a similar mesh structure and P3 scattering (the reference case). While there is evidence of small ray e ffects in the reference case far from the source, the differences between Figures 2-4 and 2-5 is dramatic. We al so found that there was a small difference (on the order of several percent, depending upon location) in flux solutions computed with S42 between P3 and P1 computations, but little difference noted in the problems tested between P3 and P5 anisotropy, and so P3 was implemented for all other computat ions. To further investigate the effects of raising the quadrat ure order, unbiased quadrature sets were tested using S22, S32, etc, up to S42 with P3 anisotropy. A comparison of the PENT RAN computed fluxes in Group 1 using selected quadratures and e quivalent MCNP5 calculations are given in Figure 2-6. PAGE 42 42 0 cm 5cm 10 cm 15cm Figure 2-4. 60Co -ray 3-D Group 1 scalar flux distribution for x-level de pths at 0 cm, 5 cm, 10 cm and 15 cm, computed by PENTRAN with the CEPXS gamma cross section library; an S12 angular quadrature (1 68 directions) with P3 scattering anisotropy was used as a reference case PAGE 43 43 0 cm 5 cm 10 cm 15 cm Figure 2-5. 60Co -ray 3-D Group 1 scalar flux distribution for x-level de pths at 0 cm, 5 cm, 10 cm and 15 cm, computed by PENTRAN with the CEPXS gamma cross section library; an S42 angular quadrature (1 848 directions) with P3 scattering anisotropy PAGE 44 44 1.00E+00 1.00E+02 1.00E+04 1.00E+06 1.00E+08 1.00E+10 0.00E+003.00E+016.00E+019.00E+01 Distance (cm)Flux Sn12 Sn22 Sn32 Sn42 MCNP5 As shown in the Figure 2-6, the PENTRAN results with S42 Legendre-Chebychev quadratures are in good agreemen t with the MCNP5 resu lts, and differences between the two are within the statistical uncertainty of the Monte Carlo in air and the water phantom distal from the source. Overall, this figure shows that incr easing the angular quadrature order results in a progressive improvement in solution accuracy, and that a substantial number of directions are required for proper representation in air region between the source and the water phantom. We further note that while the S12 level-symmetric quadrature is quite adequate to represent neutron scattering in nuclear engineering neutronics a nd reactor physics applications, it does not produce accurate results here. For a general problem w ith a high degree of angular dependence, the rayeffects must be largely mitigated by employing a large number of directions, i.e., a high-order angular quadrature, From our analysis, a minimum of S32 Legendre-Chebychev quadratures should be considered. Investigat ion of alternative, less computationally in tensive methods for mitigating the ray-effects in problems of this t ype has shown that ordinate biasing methodologies can also be of benefit19. Figure 2-6. Scalar fl ux distributions using S12, S22, S32, and S42 with P3 scattering anisotropy in PENTRAN using the BUGLE-96 gamma library PAGE 45 45 2.3.2 Investigation of Spatial Discretization (Meshing) In addition to proper resoluti on of the angular variable, the spatial variable must be appropriately discretized to represent the materi als and geometry of the system to yield the solution to the BTE. In PENTRAN, a block ad aptive Cartesian three-dimensional hexahedral mesh arrangement is used to define the spatial-ma terial distribution of the system, i.e. the system is partitioned into coarse mesh blocks containing fine mesh cells having a uniform cross-section in each fine mesh cell. The solutions to the SN equations are obtaine d by sweeping along each cell and discrete direction on the unit sphere; this process is usua lly referred to as a transport sweep 11. In discretizing a spatial mesh, it is of ten necessary to examine, using successive computations, the effect of the coupling between angle, mesh size, cross section, and the spatial differencing scheme used to solve SN equations. Here, we will briefly highlight how spatial differencing schemes are applied in the SN method. Equation (1-6), is integrated over the volume ,,ijkijkVxyz to yield volume and surface averaged flux terms. Rearranging, we can obtain the zeroth moment balance equation (omitting group subscripts for brevity): AA z in z out m y in y out m x in x out mq z y x (2-1) Equation (2-1) contains seven unknowns. Three of which, entrant fluxes (in) are assumed to be known. Also, the collectiv e cell averaged angular volumetric source A q (which includes scatter and fixed source s) is assumed to be known fr om previous source iteration. Therefore, the average angular flux A and the exiting ( out) surface angular fluxes are unknowns, where these values are derived using a set of auxiliary equations which enable PAGE 46 46 coupling of the average and boundary angular fluxes. Linear weighted sc hemes use the auxiliary formulations in Equation (2-2): 1 1 a ax inA x out (2-2a) 1 1 b by inA y out (2-2b) 1 1 c cz inA z out (2-2c) The standard Diamond Differenci ng (DD) scheme results when a=, b=, and c= in equations (2-2a, 2-2b, and 2-2c ); the DD scheme is second or der accurate, but may lead to negative solutions11. Petrovic and Haghighat developed the Directional Theta-Weighted scheme DTW37 where a can be expressed as in in in in () 1 22 mmm A xmyz mm xq xyz a yz (2-3) Similar weights can be derived along the y and z -axes for b and c respectively. DTW can remedy the unphysical oscillations inherent in solutions rendered by the DD scheme, and it is always positive. In cases of high angular stre aming, exponential aux iliary equations may be needed. These are implemented into PENTRAN to represent streaming behavior of the angular flux across the spatial cell. The expo nential Directional Weighted (EDW)38,39 can be employed in this situation with the following auxiliary formulation: (2-4) 111(,,)exp(()/||)exp(()/||)exp(()/||)moimjmkmxyzaPx Py Pz PAGE 47 47 Multiple differencing schemes in PENTRAN have an adaptive structure between successive schemes of increasing accuracy within a coarse mesh containing many fine meshes called adaptive differencing strategy (ADS): if one begins with the Diamond scheme (Diamond Zero (DZ)), the detection of a negative flux-fix-up causes select ion of the Directional Theta-Weighted (DTW) scheme, and if weights (as in Equation 2-3) are beyond certain values, selection of the Exponential Dire ctional Weighted (EDW) scheme is made. In recent PENTRAN versions selection of the Exponential Directi onal Iterative (EDI) method is made as an improvement over EDW. This ADS enables the code to automatically select the most appropriate differencing scheme within a particular coarse mesh We examined the effect of different sizes of mesh cells with various differencing schemes in our model. Figure 2-7 depi cts the fourth z-level mesh structure (plotted at hexahedral cell cente rs) relative to the full-size model; note that sparser meshing was used in the surrounding air region. Figure 2-7. 60Co model spatial meshing, axial level 4 The size of the fine mesh cells in the ai r region of the problem was of particular importance due to low density; smaller meshes yield high accuracy, but can lead to more pronounced ray effects; larger meshes mitigate th e ray-effects, but may be inaccurate. Figure 2-8 PAGE 48 48 shows a comparison between fine meshes side size of 1 cm, 2 cm and 5 cm in the air region with the DTW scheme fixed in air. These results demons trate that a 1-cm mesh in air is necessary to maintain accuracy. Based on these observations, we designed our PENTRAN model as follows; the model is 30x30x70 cm3, and consists of 45 total coarse mesh blocks arranged among 5 axial levels, containing a total of 126, 000 fine meshes. In recent PENTR AN versions, selection of the Exponential Directional Iterat ive (EDI) method is made as an improvement over EDW39. The geometry of the model has been represented using the block adap tive meshing capability of the PENTRAN code that allows the user to use a hi gher mesh density in the regions of interest. Axial level thicknesses were determined based on major components of our model, and the X-Y spatial dimensions of coarse meshes (contai ning various densities of fine meshes) were uniformly 10 cm on each side. The DTW differencing scheme was locked for the regions with low density material (air regi ons), and the adaptive differencing scheme (with adaptive upgrade potentially up to the EDW differencing scheme) was prescribed for higher material density regions. A mesh coupling scheme, Tayl or Projection Mesh Coupling (TPMC)40, was used to link discontinuous meshes using first order projectio ns as radiation sweeps through the problem. The Cartesian mesh sizes for the full size model were selected to preserve component volumes, also taking into account the optical thickness of each material. The fine-mesh intervals in different problem regions are: Source0.5 cm; lead0.5 cm; water1.0 cm; Air1.0 cm. Smaller mesh sizes better represen t curved boundaries and irre gular structure in the model, but require larger memory and may lead to longer running time (depending upon the spectral radius of the equations). Smaller mesh sizes were selected for the source region of the problem, where the particle flux has a strong gradient. Mesh sizes in the air can be made larger near the edges of the problem where a more detailed solution is not required. PAGE 49 49 Figure 2-8. Scalar flux comparison of variation in the interval of 3-D mesh cells containing air; fine meshes of 1 cm size (along each axis) in air proved to be best in providing high detail flux distribution Using the ADS described earlier with different mesh sizes, we were able to achieve a converged solution. Figure 2-8 depicts the solution behavior versus mesh size in the air region (even though there is very little optical th ickness or attenuation by the ai r); this can be explained by the fact that larger mesh sizes in the air do not adequately propagate the radiation from the source, and therefore result in so lutions with significant truncation error, particularly near the source region, while using smaller mesh sizes in air can provide good solutions with fine detail over the problem domain. 2.3.3 Investigation of Cross-Section Libraries: After satisfying the requiremen ts for properly representing a ngular quadratures and spatial discretization, the use of an a ppropriate cross section library is also important in solving problems in medical physics. Th e scalar fluxes resulting from the deterministic calculations using different cross-section libra ries have been examined (Figur e 2-9) and were compared with PAGE 50 50 similar values generated using MCNP5. To best present comparisons with the Monte Carlo results (since Monte Carlo statistical errors were lowest (<3% relative) along the central axis), we consider a slice through the central axis of the model. It is im portant to note that even with the use of volumetric mesh tallies that yield global fluxes across a grid spanning the Monte Carlo geometry, the statistical uncer tainty associated with each mesh tally site, particularly those distal or off-axis from the source region, can routinely exceed 40% or more, limiting their utility in practical applications. Figure 2-10 shows that for this particular prob lem, the relative differences in scalar fluxes between all deterministic solutions and the M onte Carlo reference increased adjacent to the gamma source. Also, for this problem, deterministic results using BUGLE-96 and CEPXS demonstrated good results in this water phantom and results were w ithin the statistical uncertainty of the Monte Carlo method. Subsequent investigation of these solutions could not attribute this behavior to grid density, quadrature, or anisot ropic scattering order. We did determine that the relative fluxes of the Monte Carlo results (converged to <3% relative error), particularly near the source, were strongly dependent upon the precise method by which the geometric volume tallies were defined, and this could easily account for the differences computed near the source. On the other hand, cross section libraries such as CEPXS provide flexibility in defining energy group width not available in BUGLE-96. In general, differences can be attributed to the energy group structure of the libraries. Result s here demonstrate the importance of library benchmarking exercises for specific problem applications. This may also suggest that in some cases the use of an equivalent surface source fluence prescribed on the entering phantom surface may be more practical to implement for dosimetry in the phantom, discussed in the next section. PAGE 51 51 A B C Figures 2-9. Scalar flux dist ributions generated by PENTRAN for a single energy group of each library using the MCNP5 solution as the refe rence case. These fluxes span the central axis in the three layers of the model. In ge neral, the deterministic results yield a very good match to the Mon te Carlo results PAGE 52 52 Figure 2-10. Relative deviation of energy group 1 discrete ordi nates results with Monte Carlo results. Results using PENTRAN with BUGLE-96 were overall best in the water phantom, and differences were within the Monte Carlo statistical uncertainty 2.4 Absorbed Dose Calculations For the energy beam used in this study, the maximum electron range produced by the primary photons is approximately 0.6 cm (as com puted by the Tabata-ItoOkabe (TIO) formula, estimate for electron range, 41), which is less than the size of th e smallest fine mesh used in the water phantom. This allows electrons produced to deposit their energy within the fine mesh, achieving CPE within the phantom volume. By applying this collisional kerma approximation for absorbed dose computation based on the fluxes calculated using the SN method with PENTRAN, we compared the SN based absorbed dose estimates with absorbed doses calculated based on MCNP5 track length (F6) tallies, which estimate the energy deposited averaged over a given cell volume and are equivalent to the co llisional kerma, this is shown in Figure 2-11. PAGE 53 53 Dose Rate Along The X-axis (5 cm Depth)1.00E+00 1.00E+01 1.00E+02 1.00E+03 1.00E+04 0102030X-axis(cm)Dose rate Sn MC Dose Rate Along the X-axis (15cm Depth)1.00E+00 1.00E+01 1.00E+02 1.00E+03 1.00E+04 0102030X-axis(cm)Dose rate Sn MC Water Phantom Dose Rate Along The Z axis1.00E+00 1.00E+01 1.00E+02 1.00E+03 1.00E+04 051015202530Depth(cm)Dose rate MC Sn A B C Figure 2-11. Absorbed dose rate comparison between SN method using kerma approximation and MC techniques using F6 tally, MC statistical uncerta inty were less than 3% on average For the energy range used in this problem, we demonstrate that the absorbed dose rate calculated using the kerma approximation for the SN method using the reference case (P3 anisotropy, CEPXS gamma cross section library, and S42 angular quadrature) is comparable to that calculated using F6 tall y in the MCNP5 code. However, for high energy photons, typically used in radiotherapy, the secondary electrons ge nerated can have large ranges exceeding the fine mesh interval, consequently vi olating the requirements for CPE, and subsequently invalidating PAGE 54 54 the collisional kerma approximation as a basis fo r absorbed dose computation. In order to use deterministic techniques in high energy photon dosim etry applications, some special treatment for secondary electron transport is therefore re quired; this will be solved using the EDK-SN method. 2.5 Parallel Decomposition For solving the 60Co model problem in a parallel environment, different domain decomposition strategies were cons idered. In Table 2-1, we pr esent the CPU times and memory requirements for the PENTRAN calculations us ing different quadrat ure sets. Parallel decomposition results for solution of the 60Co problem model are also presented in Table 2-1 based on a constant number of meshes of 126,000 with P3 anisotropy for a single energy group and different angular quadrature orders. From the table data, it can be observed that increasing the number of directions scales up running times to the processor workload. More work on each machine lowers the relative impact of message passing synchronization and latency, relative to additional computational work, as evidenced wh ere eleven times added workload only required 9.3 times as much computati onal effort (referring to the S42 calculation). This is a useful illustration that emphasizes that w ith parallel computation, added model fidelity may not actually cost additional wall clock time to yield a more robust solution. Table 2-1. Total number of directions, cumula tive problem time** required for a single energy group, P3 scattering, 126,000 fine mesh cells using 12 processors on a parallel cluster Directions and time ratio are referred to the S12 quadrature set. ** Non dedicated running time PAGE 55 55 2.6 Discussions and Findings In this chapter, we highlight ed recent studies examining th e discretization requirements of the deterministic SN method compared with the Monte Ca rlo method with application to dosimetry problems. In doing so, we chose to simulate an idealized 3-D 60Co model using the PENTRAN code. Since radiation transporting in such applications encounters little interaction in low density materials, such as air, attention was given to the quadratur e level and mesh size distribution needed for solution accuracy through the air and in a water phantom. As expected, radiation transport in this gamma ray application is highly directional, and that a traditional SN quadrature, such as an S12 level-symmetric set (168 total directions), considered robust for many neutron reactor physics applications, is in sufficient for modeling gamma rays in the 60Co model. This limitation was overcome by the use of the higher levels of quadrature, up to S42 (1848 total directions per cell) with Legendre-Chebychev quadrature sets and P3 anisotropy. We also noted that at least an S32 angular quadrature leads to the elimin ation of most ray effects and scalar fluxes within the statistical uncertainty of Mon te Carlo. Moreover, model mesh sizes were selected to preserve system volumes while taking into account the optical thickness of each material. Smaller Cartesian meshes do better re present the curvature of the model, but they require higher quadrature orders, and therefore larger memory and longer running times. Even in optically thin air regions, a mesh cell interval of on the order of 1 cm wa s necessary to yield an accurate solution, and more cells were necessary where the particle flux encountered strong gradients; cell intervals may be relaxed near boundaries where deta il is not required. Also of importance was the flexibility of the ADS in PENTRAN, permitting more appropriate spatial differencing schemes. Scalar fl ux distributions obtained from PENTRAN agreed with the Monte Carlo prediction when simulating equi valent geometry models of the 60Co model with three different multigroup cross-section libraries. Still, realistic representations of the energy spectra PAGE 56 56 for dosimetry applications may require fine ener gy group structures; it is clear the applicability of the cross section library shoul d be investigated for each problem considered. Development of a standardized multi-group cross-section library fo r these kinds of applications would be of benefit to minimize potential inconsistencies in th e application of curren tly available libraries. The results of parallel deterministic PENTRAN c ode were obtained within comparable running times to parallel MCNP5 Monte Carlo calculations for tally sites adjacent to th e source. Overall, we demonstrated that with an adequate quadrat ure set, cell mesh interval, and cross-section library, it is possible to produce results with th e PENTRAN discrete ordinates code that agree with the Monte Carlo prediction. A major advantage of the SN method is that it provides a detailed, accurate flux distribu tion at thousands of cells throughout the system, while the Monte Carlo method is best for providing accurate values at a few selected points near the source. Both methods can be advantageous in these types of applications, prov ided the level of SN discretization is adequate. This chapter laid out the founda tions necessary for accurate SN discretizations in general, although only a homogenous phantom was consider ed. The need to implement heterogeneous media in the SN simulation is essential for any medical phys ics applications. In the next chapter, we set the foundation for the new PENTRAN-MP c ode system. In this system, we implement a state of the art series B phantom into our methodology by do wn-sampling these phantoms into accurate phantoms to generate an equivalent PENTRAN input. Moreover, we calculate the absorbed dose by applying a kerma approximation for diagnostic procedures. Then, the EDK-SN method is developed, presented, and applie d in subsequent chapters. PAGE 57 57 CHAPTER 3 WHOLE BODY DOSIMETRY SIMULATIONS WITH THE UF-B SERIES PHANTOM USI NG THE PENTRAN-MP CODE SYSTEM FOR LOW ENERGY PHOTONS This work describes our 3-D SN Medical Physics code syst em PENTRAN-MP in support of ongoing development research in larg e scale 3-D absorbed dose computations42. In the PENTRAN-MP code system, Figure 3-1, an initi al high-resolution data phantom is obtained from state of the art 3-D whole-body voxeli zed human phantoms based on Computed Tomography (CT) images43. A pediatric patient (an 11-year-old male) 3-D CT image is then down sampled (collapsed) into a coarser one. Th e collapsed model was cast into a 3-D spatial distribution and sub-divided into coarse meshes c ontaining fine meshes with constant densities. SN calculations were performed using PENTRAN a nd Monte Carlo results were generated with MCNP5. To accurately compare our deterministic results with the Monte Carlo, we also generated an equivalent MCNP5 31 model implementing conti nuous energy ENDF/B-VI crosssection data libraries. Post pr ocessing in the PENTRAN-MP code system includes seamless parallel data extraction using the PENDATA code, followed by th e application of the 3-D-DOSE code to determine absorbed dose in each voxel of the 3-D phantom based on fitted mass-energy absorption coefficients 44, hence yielding a whole-body absorbed dose distribution. This chapter is presented as follows: Sect ions 3.1, 3.2 and 3.3 highlight the methodology for accomplishing our deterministic calculations. Section 3.4 pres ents a computational human phantom example implementing the methodology of PENTRAN-MP a nd an equivalent MCNP 5 model. Section 3.5 compares the Monte Carlo and deterministic soluti ons and analyzes the parallel performance of the PENTRAN. Finally, Section 3.6 discusses the findings of this Chapter. PAGE 58 58 Figure 3-1. PENTRAN-MP Code System Structure Flowchart 3.1 Methodology In this work, we must determine the radi ation transport flux and absorbed dose in a computational tomographic (or voxelized) phantom using deterministic techniques. Central to this is the application of deta iled anatomical models in whic h we can render accurate absorbed doses. Therefore, coupling of these anatomical models with robust 3-D deterministic transport simulations can yield critical or gan and whole body absorbed dose pr ofiles from patient-specific radiotherapy treatments, includi ng absorbed dose delivered to colla teral tissue located away from the exposure site. To properly attribute absorb ed doses in our models, we must set up a discretized phase-space based on a detailed phantom for the material spatial distribution. Due to PAGE 59 59 interest in fully characteriz ing diagnostic absorbed doses to younger patients, for this presentation, we chose a pediatri c phantom. Also, to ultimately render absorbed doses in the phantom wherever we desire, there are a number of codes that constitute the PENTRAN-MP code system that include preprocessing to establish an x-ray source, appropriate interaction cross sections, phantom spatial mesh, followed by the radiation transport calculation, and postprocessing computation data to yi eld absorbed doses and absorbed dose volume histograms, etc. The details of this methodology are pr esented in the paragraphs below. 3.2 UF Series B Pediatric Phantoms The UF Advanced Laboratory for Radiatio n Dosimetry Studies (ALRADS) has a well established background in work w ith detailed human phantoms and phantom data. Following the development of their CT-Contours segmentation so ftware, and the construction of their first tomographic dosimetry model (UF Newborn) the ALRADS research team focused on constructing an expanded series of pediatric tomographic phantoms (Series B phantoms), this time using live patient CT images from the Sha nds image archive. A total of five additional phantoms have been constructed from CAP (chest abdomen pelvis) and head CT scans of pediatric patients examined at the University of Florida (UF) Shands Childrens Hospital. Patient inclusion criteria included those patients displaying normal or near-normal anatomy at their respective age and gender (e.g ., avoidance of regions of tumo r, edema, organ enlargement) Table 3-1. While the previous phantoms preserved the body dimensions and organ masses as seen in the original patients who were scanne d, comprehensive adjustments were made for the Series B phantoms to better ma tch International Commission on Radiological Protection (ICRP) age-interpolated reference body masses, body heights, sitting heights, and internal organ masses43. PAGE 60 60 Table 3-1. Computed tomogr aphy image sources for the development of the UF pediatric phantom series Figure 3-2. The UF Series B pe diatric phantoms series (Newborn, 9-months, 4-years, 8-years, 11-years and 14-years old patients)43 For this work, a pediatric patient (an 11-year-old male, depicted in Figure 3-2) is voxelized with high resolution using 0.94 mm 0.94 mm 6. 00 mm voxels as a case study. This geometry data was then collapsed into a 4.7 mm 4.7 mm 30.0 mm voxelized model using GHOST-3D, described in section 3.3. Subsequently, using PENMSH-XP, the co llapsed model was cast PAGE 61 61 into a 3-D spatial grid distribution and sub-divi ded into coarse mesh grids, each containing a corresponding number of fine mesh cells. 3.3 PENTRAN-MP Code System Currently, the PENTRAN-MP code system is composed of 3 levels of calculation, including pre-processing algor ithms that include the GHOST -3-D phantom voxel collapsing (back-thinning) code, the DXS x-ray source generation code 45, the PENMSH-XP 3-D mesh distribution and input generato r code, and the CEPXS code (a product of Sandia National Laboratories) for generating multigroup cross sections 33. The radiation transport calculations are carried out in PENTRAN, and post-processed using PENDATA to gather parallel distributed data, and 3-D-DOSE to render user specified absorbed dose data. 3.3.1 Pre-processing Codes 3.3.1.1 The collapsing code (GHOST-3-D) The GHOST-3-D, developed by Ghita and Al-Bas heer, is a code used to input CT-scan binary data obtained from a high resolution voxelized human phantom. There are several constraints in simulating mega-voxel models, including balancing limited computer memory (even with the benefit of parallel memory storag e, as in PENTRAN) with proper discretization, etc. To optimize accuracy and computational wo rk in representing the voxelized phantom, we developed a code, GHOST-3-D, that provides an equivalent model with a reduced number of meshes (subject to user-derived accuracy targ ets), overcoming the mentioned difficulties. 3.3.1.2 The mesh generator code (PENMSH-XP) PENMSH-XP35 by Yi and Haghighat is a 3-D Cartesian-based mesh generator that prepares discretized material and source distribu tions for the PENTRAN particle transport code. PENMSH-XP is capable of performing two di fferent methodologies in generating PENTRAN input files. In the first methodol ogy, a physical model is partitione d into slices al ong the z-axis, PAGE 62 62 referred to coarse z-levels, and then each z-level is partitioned into single coarse levels of 3-D x-y-z fine mesh grids. Predefined shapes can also be used to ge nerate the desired shapes in a particular model. A second met hodology is based on two binary input files, one for the phantom geometry, which can be established from real CT data, and the other containing the source density for each voxel/fine mesh, which can be used to specify a fixed-source distribution in the model. By combining the two files, PENMSH-X P directly generates an input deck for the PENTRAN transport code. 3.3.1.3 Source spectra and cross-sections generation Based on the bremsstrahlung and ch aracteristic parametric functi ons of the model proposed in 1991 by Tucker, Barnes, and Chakraborty 46 the DXS code, principally developed by Monica Ghita, generates x-ray spectra in the diagnos tic range (50-140 keV), according to the input parameters (tube potential, anode angle, type and amount of filtration), and accommodates them to any energy group structure. Optionally, as part of the PENTRAN-MP package, DXS also generates the input files for CEPXS (Sandia Na tional Lab), the x-ray cross-section generator code, corresponding to a desire d energy group structure and PN order. Group-dependent crosssections rendered by CEPXS are then stripped fr om the output file and adapted to a format suitable for the PENTRAN deterministic calculati ons using another code developed by Ghita, GREPXS47. 3.3.2 Transport Calculation Transport calculations were performed using the PENTRAN SN and MCNP5 Monte Carlo as discussed in Chapter 2. PAGE 63 63 3.3.3 Post-Processing Codes 3.3.3.1 Data management utility (PENDATA) Immediately following a transport calculation, a data management utility, PENDATA, seamlessly gathers parallel stored data automatically following a parallel PENTRAN run, and provides several options for the user in stripping results from parall el output files, including data extractions from binary file storage which is an essential utility for a c ode with parallel I/O. 3.3.3.2 Dose calculation code (3-D-DOSE) The 3-D-DOSE code, developed by Al-Basheer is used to re nder absorbed doses based on photon fluence results derived fr om the PENTRAN transport calculations. In doing so, a methodology for determining absorbed dose in the phantom based on photon fluence was developed for situations where electron transport effects do not require specific treatment. CPE is assumed to be established in a volume when the energy carried into the volume by charged particles is balanced by the energy carried ou t of the volume by the charged particles. In the example simulation in Section 3.4, the maximum energy used is 90 kV, for which the range of the secondary electrons produced is smaller than smal lest voxel dimensions used in the model. Moreover, the atomic number of the material Z is relatively low, excluding major production of bremsstrahlung radia tion, thus depositing all th e photon energy in the voxel, achieving CPE. Consequently kerma is equal to absorbed dose in this case. We note this is not true for regions in close proximity to the su rface of the model using higher energies due to absorbed dose buildup, where corrections fo r secondary electron tr ansport may be quite significant, especially for high energy beams. After extracting 3-D scalar fluxes from the phantom model we apply the 3-D-DOSE c ode. Based on fitted mass-energy absorption coefficients 44, shown in Figure 3-3, and the scalar fluxes from the PENTRAN calculations, 3-D- PAGE 64 64 Lung(ICRU-44)1.00E-02 1.00E-01 1.00E+00 1.00E+01 1.00E+02 1.00E+03 1.00E+04 1.00E-031.00E-021.00E-011.00E+001.00E+011.00E+02Energy(Mev)Mass-absor p( cm2/ g) ICRU-44 Fit Cortical Bone (ICRU-44) 1.00E-02 1.00E-01 1.00E+00 1.00E+01 1.00E+02 1.00E+03 1.00E+041.00E-031.00E-021.00E-011.00E+001.00E+011.00E+02Energy(Mev)Mass-absor p( cm2/ g) ICRU-44 FIT Soft tissue (ICRU-44)0.01 0.1 1 10 100 1000 10000 1.00E-031.00E-021.00E-011.00E+001.00E+011.00E+02Energy(Mev)Men/roh ( cm2/ g) ICRU-44 FIT DOSE calculates absorbed dose in each voxel yi elding a whole-body absorbed dose distribution and dose-volume histograms for cr itical organs of interest. Figure 3-3. Values of the mass energy-absorption coefficient, and the fitting function as a function of photon energy, for Lung (ICRU44), soft tissue (ICRU -44) and cortical bone (ICRU-44) 3.4 Computational Example For this work, a pediatric patient (an 11-year -old male) was voxelized with high resolution using 0.94 mm 0.94 mm 6.00 mm voxel size (398 voxel 242 voxel 252 voxel) as a case PAGE 65 65 study. This geometry data was then collapsed into a 4.7 mm 4.7 mm 30.0 mm voxel size (79 48 50 voxels) model using GHOST-3-D. Subs equently, using PENMSH-XP, the collapsed model was cast into a 3-D spatial distribution a nd sub-divided into six coarse mesh z-levels containing five coarse meshes in x-y domain, and these containe d 6320 fine mesh cells (Figure 3-4, 3-5). An isotropic volumetric (17.8 cm x 3.65 cm x 30 cm) x-ray source was placed above the upper left side of the phantom. Although the model has relatively sm all fine mesh sizes, using S32 Legendre-Chebychev quadrature, the adaptive differencing, and CEPXS cross-sections with P3 anisotropy, we were able to eliminate any unphysical oscillations that may result in the deterministic solution. Using 4, 6, and 8 energy groups, the SN model was executed on 16 Opteron processors and was completed in 1.8, 2.1 and 3.8 hours, respectively. With the general purpose radiati on transport code MCNP5, we simu lated an equivalent model to the one simulated by PENTRAN. The collapsed 3-D binary files composing the UF voxel phantom was converted to ASCII formatted matrix files and then incorp orated into a MCNP5 lattice defined input. An input d eck example for the MCNP5 model is presented in Appendix B. Volumetric flux (F4) mesh ta llies were used in the Monte Carlo geometry and were equivalent to the discretized SN volumes all over the model. All simulations were performed using the parallel cluster Einstein composed of 16 processors at 2.4 GH z with 4 Gb/processor, for a total of 64 Gb of parallel memory used fo r general purposes. Parall el MCNP5 results in the source region converged with a st atistical error less than 4% in the source region on average. PAGE 66 66 Figure 3-4. UF-Series B phantom (11 year s old male) with (398 242 252 voxels); the corresponding PENTRAN back-thinned mo dels with (132 80 84 voxels) and 72 materials (upper right), and another with (79 48 50 voxels) and 66 materials (lower right) PAGE 67 67 Figure 3-5. Transversal view of the 11-year-old male Series B original phantom (upper left), collapsed into 189k voxels and 4 materials ( upper right), collapsed into 900k voxels and 4 materials (lower left) and collapse d into 900k voxels and 72 materials (lower right) 3.5 Computational Results 3.5.1 Flux Comparison of Deterministic and Monte Carlo Results Comparison of the SN and Monte Carlo results for the scalar fluxes along the source central axis and in the source region along the z axis at different depths for different energy groups were made. Figure 3-6 depicts a 3-D scal ar flux distribution computed by PENTRAN for the reference case. This simulation was based on the use of an S32 angular Legendre-Chebychev quadrature (1088 directions/mesh) with P3 scattering anisotropy using th e CEPXS with 8 energy groups. PAGE 68 68 Figure 3-6. x-ray 3-D Group 1, 3, 5 and 7 scalar flux computed by PENTRAN with the CEPXS cross section library, S32 angular quadrature (1 088 directions) with P3 scattering anisotropy Figure 3-7. Flux distributions using S32 with P3 scattering an isotropy in PENTRAN using the CEPXS library vs. MCNP5 Monte Carlo result s for energy groups 1, 2, 3 and 4 along the source central axis. Mont e Carlo uncertainties were less than 5%, 4%, 4%, and 3% respectively PAGE 69 69 Figure 3-8. Angular fl ux distributions using S32 with P3 scattering anisotropy in PENTRAN using the CEPXS library vs. MCNP5 Monte Ca rlo results for energy groups 5, 6, 7 and 8 along the source central axis. Monte Carlo uncertainties were less than 2%, 3%, 6%, and 30% respectively Figure 3-9. Angular fl ux distributions using S32 with P3 scattering anisotropy in PENTRAN using the CEPXS library vs. MCNP5 Monte Ca rlo results for energy groups 1, 3, 5 and 7 along the z axis at diffe rent depths in the phantom. Monte Carlo uncertainties were less than 5% at y=4.3 cm, 6% at y=13.5 cm PAGE 70 70 As shown in (Figures 3-7, 3-8, and 3-9), the PENTRAN results with S32 LegendreChebychev quadrature are in very good agreemen t with the MCNP5 predictions in the source region; MCNP5 mesh tallies were converged to <6 % rela tive error, and typi cally < 2% for most of the results for the reference case (using 8 en ergy groups). Overall, this figures show that deterministic techniques are capable of producing accu rate results within th e statistical error of Monte Carlo methods. Since the uncertainty is pr oportional to 1/ (number of histories), to achieve an equivalent-uncer tainty Monte Carlo simulation with converged mesh tally information 70 cm away from the source in energy group 1, we esti mate over 200 hours are needed for MCNP5 executed on 16 processors to achieve fluxes to within 5% for this group; longer times would be require d for lower energy x-rays. 3.5.2 Dose Calculations and Dose Volume Histograms (DVH). Using 3-D-Dose, 3-D absorbed dose distributions were calc ulated in each voxel in the whole model of the 11-year old male phantom model for broad-beam ex ternal photon field. Figure 3-10 shows the stages of calculation in th e PENTRANP-MP code system that results in 3D Dose distributions for the entire phantom. Figures 3-11, 3-12, a nd 3-13 show selected slices of the phantom and the corresponding absorbed dose distribution, while Figure 3-14 depicts selected organs with correspondent 3-D Absorbed dose distributions. PAGE 71 71 Figure 3-10. Simulation met hodology for 3-D absorbed dose computations using PENTRANMP code system. Left: UF series B v oxelized phantom, middle: correspondant PENTRAN input, right: 3-D abso rbed dose distribution after irradiation with a x-ray source a. PAGE 72 72 b. Figure 3-11. CT view of the 11-year-old male Se ries B original phantom collapsed into (a) 189 k voxels and the correspondent absorbed dos e distributions (b) 900 k voxels and the correspondent absorbed dose distributions a. PAGE 73 73 b. Figure 3-12. CT view of the 11-year-old male Se ries B original phantom collapsed into (a) 189 k voxels and the correspondent absorbed dose distributions (b) 900 k voxels and the correspondent absorbed dose distributions a. PAGE 74 74 b. Figure 3-13. CT view of the 11-year-old male Se ries B original phantom collapsed into (a) 189 k voxels and the correspondent absorbed dose distributions (b) 900 k voxels and the correspondent absorbed dose distributions Figure 3-14. Simulation met hodology for 3-D organ absorbed dose computations using PENTRAN-MP code system PAGE 75 75 Dose volume histograms (DVHs), for a selected nu mber of organs giving the percentage of tissue volume (organ) that receives an absorbed dose greater than the value indicated on the abscissa, were calculated and compared to thos e results from the equivalent Monte Carlo model. As anticipated, DVH results are in good agreement for the SN and Monte Carlo radiation transport simulations (Figure 3-15). However, Monte Carlo calculations can require substantial run times for low precision at o ut of field distan ces in the phantom, particularly as photons down-scatter to lower energies The Monte Carlo results comput ed values will eventually converge with increased run time as more pa rticles reach sites distal to the source. a. b. Figure 3-15. Dose Volume Histograms (D VH) computed based on PENTRAN (S32,P3) transport results using the CEPXS library vs. DVHs computed using MCNP5 Monte Carlo results. MCNP5 Monte Carlo results required 16 hour s 16 processors compared with the PENTRAN-MP SN solution, requiring <1.9 hours on 12 processors for the whole-body computation 3.5.3 Parallel Computation For the UF 11 year old male phantom problem, a decomposition strategy was used to assist the parallel performance. In Table 3-2, we present the cumulative problem time and memory PAGE 76 76 requirements for deterministic calculations disc ussed in previous sections using different quadrature sets. Table 3-2. Total number of directions, Cumula tive Problem Time required for different energy groups, P3 scattering, and 189k fine mesh cells Parallel Decomposition Number of Energy groups SN (relative workload) Number of Processors A E S Cumulative Problem Time (hours) 4 42 (2) 12 4 1 3 4.1 6 24 (1) 12 8 1 2 2.1 8 32 (2.32) 12 4 1 3 3.69 8 32 ( 2.32) 16 8 1 2 3.06 Parallel decomposition results for solving the phantom are presented in Table 3-2 based on a constant problem phase space of 189,000 cells with P3 anisotropy and diffe rent numbers of energy groups numbers and several angular quadratures. A, E, S refers to the angular, energy group, or spatial decomposition applied to solve th e problem using parallel computing. From the table data, it can be observed th at increasing the number of dire ctions requires longer running times that essentially scale with processor workload. 3.6 Discussion and Findings We have demonstrated a methodology in corporating a discrete ordinates (SN) approach that is an effective alternative to the Monte Carlo method for accurately solving the transport equation over an entire human phantom. Models he re were based on the UF series B pediatric phantoms developed by the University of Florid a ALRADS group. We simulated the complete phantom discretized using 900k/186k voxe ls, and demonstrated that the SN method in the PENTRAN-MP code system, developed as a te am effort, readily provides accurate whole body fluxes. Radiation transport throughout the phantom encounters few interactions in low density PAGE 77 77 materials, such as air. Attention was given to the quadrature level, using at least an S32 LegendreChebychev angular quadrature, yielding 1088 total di rections per mesh-cell; we were able to accurately obtain a solution analogous to the MCNP5 Monte Carlo calcula tions. Also, the tally results from the eight energy gr oupsrepresenting the spectrum a nd using the CEPXS multi-group cross-section library, were co mpared to the continuous en ergy MCNP5 Monte Carlo results using the ENDF/B-VI data libraries. Overall, with an adequate quadrature, mesh si ze, and cross-section libr ary, it is possible to produce results with the PENTRAN discrete ordinates code that agree within the statistical uncertainty of the on-axis MCNP5 Monte Carlo ca lculations for the proposed model. Parallel deterministic PENTRAN results were obtained within comparable running times to parallel MCNP5 Monte Carlo calculations for ta lly sites adjacent to the source. A major advantage of the SN method is that it provides a detailed, accurate flux distribution at thousands of cells throughout the sy stem, while the Mont e Carlo method only provides highly accurate values for selected points near the source. Though the use of volumetric mesh tallies can provide global fluxes across a grid spanning the Monte Carlo geometry, the statistical error asso ciated with each mesh tally site that would require, according to estimates, more than 200 hr for tallies 70 cm away from the source to converge to less than 5%., limits their utility in prac tical applications. An equivalent Monte Carlo calcula tion to obtain the flux distribution with reasonable statistical sampling error over the entire problem space may require orders of magnitude longer running times compared to parallel deterministic SN calculations. While more tests are needed, the recently developed Exponential Directional Iterative (EDI) SN differencing scheme was excellent in d ealing with fine inhomogeneous material PAGE 78 78 densities inherent in medical physics problems, properly handling numerical radiation transport physics for diagnostic energy beams. Overall, for the energy ranges used in th e problem, the secondary electrons need not specifically be modeled using charged particle transport methods to accu rately determine global absorbed dose distributions. This is due to the fact that for an SN calculation with an appropriate numerical fidelity (evidenced by the excellent agreement with Monte Carlo), secondary electrons that result from photon interac tions in the phantom have equi valent ranges smaller than the minimum SN spatial grid interval used to represent the anatomical structure. Alternatively, at higher photon energies (not consider ed here), this is not the case, and charged particle transport and corresponding interactions yielding energetic electr ons must be considered over the problem phase space to yield an accurate absorbed dose distribution. In Chapter 4, a new EDK-SN approach is presented to account for higher energy beams and sec ondary electron absorbed dose. PAGE 79 79 CHAPTER 4 ELECTRON ABSORBED DOSE KERNELS TO ACCOUNT FOR SECONDARY PARTICLE TRANSPORT IN DETERMINISTIC SIMULATIONS In previou s chapters, Charged Particle Equilibrium (CPE) and its limitations were presented, and that in cases with high energy photons, corre ctions for secondary electron transport may be significant. As mentioned, this is readily treated in Monte Carlo codes, yet quite difficult to treat explicitly in deterministic codes due to th e large electron ranges and added numerical complexities in reaching conve rgence in photon-electron transport problems10 The goal here is to produce a methodology to accurately estimate organ absorbed doses and attribute whole body absorbed doses incurred from irradiations characteristi c of external beam therapy using the CT-based anatomical patient geometri es in 3-D deterministic discrete ordinates (SN) radiation transport models. Accurate absorbed doses will be rendered using pre-computed, high resolution Monte Carlo based electron kernels. In doing so, the flux and the curren t provided by the SN method will be used to both scale and pr oject the absorbed dos e using the Monte Carlo based kernels, and map it to the surrounding voxels. Hence, this a pproach preserves high fidelity and accuracy, accounting for complete photon transpor t, with the final absorbed dose delivered from electron transport propagate d by coupling with 3-D determinis tic photon transport; this is what makes this approach novel. In this chapte r, water phantom results using a 0 to 8-MeV step uniform beam will demonstrate that absorbed doses can be accurately obtained within the uncertainties of a full physics Monte Carlo simulation. 4.1 Introduction to The EDK-SN Method Various modalities are applied in absorb ed dose calculations for high energy photon beams, including empirical models that are based on measurements made in the clinic2, or theoretically based models, such as the application of Monte Carlo techniques3,48. Alternatively, as mentioned earlier, deterministic methods have been implemented to assess absorbed dose PAGE 80 80 calculations20,21, although experience has shown that 3-D deterministic electron transport is quite challenging in order to yield a converged solution. To properl y treat 3-D electron transport physics deterministically, yet still achieve reasonably fast and accurate whole body computation times for simulations with high energy photons, we have developed angular-energy dependent transport electron dose kernels, or EDKs that can be coupled to SN solutions. These kernels were derived via full physic s MCNP5 Monte Carlo el ectron transport simulations, and are applied based on deterministic photon soluti ons over the problem phase space, thereby accurately accounting for the abso rbed dose that is ultimately rendered by charged particle electron transport. This is possible because th e PENTRAN discrete ordi nates code preserves detailed photon flux and current data on a voxel -by-voxel basis, so that accurate whole body absorbed doses may be rapidly achieved for hi gh energy photon sources by performing a single deterministic SN multigroup photon calculation on a parallel cluster. This is followed by linking the SN-derived photon fluxes and net currents to our Monte-Carlo-based El ectron Dose Kernels to account for a full-p hysics absorbed dose; we call this the EDK-SN methodology for determining total body dosimetry. This total body dosimetry enables dose to be computed across the entire phantom phase space, wherever it is required, presented in the following sections. 4.2 Monte Carlo Based Dose Kernels Using Monte Carlo calculations, electron abso rbed dose kernels for pre-determined photon energy groups were generated; energy was deposited in a collective cloud of voxels ( i, j, k ) composed of tissue, where absorbed doses attributed to these voxels are directly caused from incident primary photons in a given energy gr oup that propagate electrons from voxel, ( i',j,k ) the dose-driving voxel or DDV When producing these electron absorbed dose kernels, the photons having energies lower than the limits of the energy group are cut off as they are PAGE 81 81 generated to prevent a cascadi ng contribution to the EDK for lower energy groups. These kernels were then applied along the direction of the ne t current and scaled by the magnitude of the photon flux on a fine mesh basis pr ovided by the PENTRAN discrete ordinates code, as well be explained in later sections. Figure 4-1 illustrates the DDV, with an indicated direction of photon current. The Electron absorbed Dose Fractional kernel in voxel ( i, j, k ), given as the EDFg, due to a primary photon group, can be determined in term s of the initial photon energy for a particular energy group g By partitioning the energy deposited in voxel ( i,j,k ) into multiple energy bins that are aliased to the SN energy group structure, we can pred ict the fractional electron absorbed dose kernel contribution per unit photon flux per source particle: (,,)(,,)(',',')gEDFijkEDKijkijk gg MC (4-1) Therefore, the SN transport-guided absorbed dose rate can then be obtained for any voxel location ( i, j, k) by summing the product of the EDFg( i,j,k ), the particle normalization factor ( ), and the SN scalar flux (,,)NgSijk divided by the voxel mass M (i, j, k): (,,)[(,,)](,,)/(,,)Ng ggsS s D ijkEDFijkijkMijk (4-2) This procedure will enable the detailed info rmation to be obtained on the absorbed dose deposited in the model by scaling the EDKs by the photon flux and pr ojecting them by along photon net current vectors computed for each voxe l. This will directly account for secondary electron physics in each voxel to achieve an accurate total body attributed absorbed dose. The particle normalization factor in equation 4-2 correctly scales the number of voxels linked using the EDK, and is set by the deterministic grid spacing. We note that this factor is specific to the grid used to estimate the EDK-SN. Consequently, to determine this factor, the absorbed dose is calcu lated through the EDK-SN model and should be scaled to the absorbed PAGE 82 82 dose calculated via a Monte Carlo full physics absorbed dose in an equivalent model. This factor is a kernel-specific f actor, and once determined, it can be used for any execution implementing the corresponding electron kernels, provided grid spacing is consistently used. Figure 4-1. Procedure to produce E DK based on Monte Carlo Calculations 4.3 EDK Determinations Based on The MCNP5 Monte Carlo Code In the MCNP531 Monte Carlo code, there are three different methods by which one can tally to yield absorbed doses to be used for inte grated organ absorbed dose calculations; these are achieved by using either a cell flux tally (F4), an energy deposited tally (F 6), or a pulse height tally (*F8), and each has a different basis for computing the absorbed radiation dose14,31. To illustrate this, we compare absorbed doses computed using the F4, F6 and *F8 tallies, respectively to highlight the differences in thes e computations, and establish a procedure that ensures that our electron kernel s are calculated accura tely, and the Monte Carlo transport and subsequent absorbed doses determined as a be nchmark are completely based on a detailed physics treatment. PAGE 83 83 4.3.1 Computational Model To highlight the differences between the different tallies, used for absorbed dose estimation, in the MCNP5 code, we c onsidered a box of 11x11x11 cm3 composed of a water cube containing a 1x1x1 cm3 8 MeV forward peaked X-ray source placed at the center of the cube Figure (4-2). The size of this model allows complete deposition of secondary electrons produced from the Xray source. A flux mesh-tally, using DE absorbed dose energy and DF absorbed dose function cards, was generated to cover the entir e model. Moreover, since no mesh tally is available for *F8 and F6 tallies in the MCNP5 code, a 3-D tally distribution of cells, equivalent in purpose to the MCNP5 F4 mesh tally, were implemented in the model. This model, with the implemented tally sites, was used for highlighti ng the differences between the different tallies, and consequently was used for reporti ng a scaled full physics absorbed dose. 4.3.2 MCNP5 Tallies Comparison The source in the center of the model was app lied to generate a pencil beam impinging perpendicularly, along the z-axis, on the box of wate r. As shown in Figures 4-3, 4-4, and 4-5, the absorbed dose is highly forward peaked, with little dose resulting from back-scattered radiation. Absorbed doses resulting from the use of an F4 volume flux tally method, using a DE dose energy and DF dose function cards, compared to an F6 tally using bo th the photon-electron mode (p, e mode) and photon only (p mode) were a ll in agreement, since the two tallies depend on similar physics in predicting absorbed dose. Thes e results are depicted in Figure 4-3. On the other hand, when comparing absorbed doses com puted using a *F8 using a fully coupled photon electron transport mode (p, e mo de) and the F6 derived dose ta llies, the differences were significant, as depicted in Figure 4-4. This is expected, since F6 tallies are analogous to F4 tallies in MCNP, and can be directly aliased to assumptions equivalent to CPE, while *F8 tallies are determined based on full charged particle transport physics, as discussed below31. PAGE 84 84 Figure 4-2. Monte Carlo model used for comparisons between MCNP5 pulse height tally, energy deposited tally and cell flux tally employing dose conversion factors (DE, DF). The source is located at the dose driv ing voxel DDV Figure 4-3. Comparison between F6 (p, e mode), F6 (p mode), and F4 ( DE, DF) tallies in MCNP5, statistical uncertainty was less than 1% PAGE 85 85 Figure 4-4. Comparison between F6 (p, e mode) a nd *F8 (p mode) tallies in MCNP5, statistical uncertainty was less than 1% Figure 4-5. Comparison between *F8 (p, e mode) and *F8 (p mode) tallies in MCNP5, statistical uncertainty was less than 1% PAGE 86 86 When comparing *F8 (p mode) transport and *F 8 (p, e mode) transport, (Figure 4-5) shows significant differences in the outcomes for the sa me simple model. This directly results from accessing different calculation modes in MCNP5 with regard to electron physics. If electron transport is turned on (mode p, e), then all ph oton collisions except coherent scatter lead to electrons that are banked for later transport. A lternatively, if electron transport is turned off (achieved by omitting the e on the MCNP5 Mode card), then a thick-target bremsstrahlung model (TTB) is used. This model generates electrons, but assumes that they are locally slowed to rest. The TTB production mode l contains many approximations co mpared to models used in actual electron transport. The choice of which tally is to be used and the detailed physics treatment by which the Monte Carlo code is execu ted is essential to th e success of the EDK-SN method. In conducting computations for generati on of the EDKs, we used the full detailed physics treatment, implementing a *F8 (p, e mode) tally, as opposed to the alternative methods that lead to less robust electron physics, es sentially employing a collisional kerma. A detail discussion on photons interactions in matter is presented in Appendix A. 4.4 EDK -SN Code System Methodology 4.4.1 Net Current For any SN calculation, the PENTRAN discrete ordina tes code preserves angular information explicitly in parallel data storage arrays. This allows one to calculate netJ and to apply the precomputed electron absorbed dose kernel along th e net current direction, Figure 4-6, where: (4-3) (4-4.a) nxnynzijk netJJJJ nxnxnx J JJ PAGE 87 87 (4-4.b) (4-4.c) (4-5) Where u j is the photon net current unit vector. Figure 4-6. Deterministic calculation will provide the angular information explicitly 4.4.2 MCNP Based Electron Kernels When generating these electron kernels, a si ngle beamlet of mono-en ergetic photons, for a certain energy group, interacts at the center of a water phantom, as disc ussed in the previous section. For the purpose of this study, we generated the EDKs fo r 16 photon energy groups covering the energy spectrum from 0 to 8 MeV. Based on our previous study, as discussed in chapter 1, these energies will violate CPE for the model considered. Moreover, the width of the energy group is set to be 0.5 MeV, minimizing the effect of energy discretization in an acceptable SN running time (we note that energy group disc retizations can be optimized for better nynyny J JJnznznz J JJ net u net J j J PAGE 88 88 time performance, although this is not considered here). Kernels are computed and then scaled by photon flux in the DDV in dose kernel voxels (pro jected from the DDV) along the net current direction of the deterministica lly computed beams of same en ergy group (Figure 4-7). It is important to note that for each energy group calcu lation, the Monte Carlo model was set to cutoff photons with energies lower than the lower lim it of the energy group being modeled. This is important to prevent any cas cading contribution to the EDK to lower energy photons. Figure 4-7. For each energy group, using Monte Carlo radiation transport, a pencil beam of mono-energetic photons interacts at the center of a water phantom These Electron Dose Kernels (EDKs) can th en be implemented to enable direct attribution of the final equivalent absorbed dose initiated by the photon beamlets. The EDKs can also be placed in various tissue types and ener gies relative to a spec ific direction of photon beam determined from solution of the Boltzma nn transport equation using PENTRAN. This will permit complete attribution of the absorbed dose ultimately due to electron straggling, where the EDK coupled with SN (EDK-SN) provides an equivalent elect ron absorbed dose look-up for each photon beamlet in a discrete ordinates computation. However, to achieve precise solutions using the EDK-SN, the SN calculation should be carried out with great attention to phantom model representation and meshing, proper quadrature sets, cross s ection libraries, and PAGE 89 89 SN differencing schemes; these issu es were investigated and presen ted in Chapter 2 in assessing what critical discretizations are required to yield accura te photon transport soluations42. 4.4.3 Euler Angles Euler angles, described in (Figure 4-8), pr ovide a means of representing the spatial orientation of any frame of angular space as a composition of rotations from a given reference frame. The translations afforded by Euler angles were implemented to enable projections of the EDKs for a collection of cells aliased to a single reference pencil beam v ector, as modeled in the Monte Carlo simulations (as in Figure 4-2) to generate EDKs to an arbitrary net current vector specifying the direction of the photons derived from SN computations. In the following discussion, we denote a fixed system in lower case ( x,y,z ) and a corresponding ro tated system in upper case letters ( X,Y,Z ). We use the term "Euler Angle" for any representation of three dimensional rotations of the electron kernel absorb ed dose matrix that aligns with the direction of the net current for each DDV. Figure 4-8. Euler angles, zyz convention. Similar conventions can be obtained by following of the following order in rotation ( zyz, xyx, xzx, yzy, yxy ) Where : PAGE 90 90 is the angle between the y-axis and the line of nodes is the angle between the z-axis and the Z -axis is the angle between the line of nodes and the Y -axis. The above description works for the zyz form. Similar conventions are obtained just renaming the axis ( zxz, xyx, xzx, yzy, yxy ). There are six possible combinations of this kind, and all of them behave in an identical way to the one described before. Just imagine three rotations over the set of axis form ing the rotation matrix R [][()][()][()]zyzRRRR (4-6) Where R is the rotational element matrix Depending on the convention used for the calculat ing these angles, they can be related to the angles for the net current obtained from th e PENTRAN run. Assuming that rotated matrix indices to be known, we generated the vector A for each point from the center of the matrix, where the dose values of this vector form th e matrix dose point to be projected onto the SN matrix: ['][][] ARA (4-7) 11[]['][][][] R ARRA (4-8) 1[][]['] A RA (4-9) When calculating the original un-rotated A vector for each point, we interpolate the vector location, in 3-D, for the corresponding dose valu e based on the original EDK matrix calculated with the Monte Carlo model. This generates a 3-D translated (by rotation) electron absorbed dose kernel matrix corresponding to a particular direc tion in space. For the current code system, this process was performed for 1848 directions, gene rating 1848 files for each energy group. An example of a rotation of the EDK from a referen ce (Monte Carlo) EDK is depicted in Figures 49.a and 4-9.b. PAGE 91 91 A B Figure 4-9. EDK-matrix after rotation to direct ion cosines A: (-0.5, 0.35, 0.8) and B: (0.25,0.40,-0.88) relative to the reference coordi nate system (0, 0, 1). Dose Driving Voxel is located at (0,0,0) PAGE 92 92 4.4.4 Ordinate Projection In the discrete angle approximations, it is assumed that the EDK-SN is applied on distinct equal weight angles (,)(,,) To generate these angles, the following relations should hold for any given octant: orderN (4-10) The orderN in any 3-D EDK-SN quadrature corresponds to the number of levels from each direction cosine on the unit sphere, and there are (2)dirorderorderNNN ordinates on the unit sphere, with _(2)/8diroctorderorderNNN in each octant, with orderN/2 distinct direction cosine values. Figure 4-10 shows orderN42 used rendering 231 direction per octant. (21) 2i (4-11) (21) () 22j i (4-12) where i = orderN For any unit angular direction vector the ordinate must lie on the unit sphere: 2221 (4-13) ()()SinCos (4-14.a.) ()()SinSin (4-14.b.) ()Cos (4-14.c.) PAGE 93 93 Figure 4-10. First octant of S42 equal we ight angles used for projecting the EDK-SN 4.4.5 Final Dose Calculations The EDK-SN model specific portion, described in Figure 4-11, serves as the critical link in the EDK-SN system to accumulate the absorbed dose in each fine mesh. The algorithm determines the accumulated energy deposited in each voxel for each photon energy group, based on the SN net current vector and is scaled to the approp riate dose by the photon flux. Finally, the code accumulates all the absorbed dose contributions from all the energy groups in to a 3-D total absorbed dose matrix. Figure 4-11 shows a flow ch art representing the EDK-SN methodology. An energy dependent EDK was generated using high reso lution direction projection sets and a MC pre-computed EDK for a single direction. These ke rnels were generated over the uni t sphere of ordinates via the Euler Angles described in Section 4.4.3. Concurrently, when generating the EDK, a FORTRAN90 module will be generated for later use in the dose calcula tion portion of the PAGE 94 94 Figure 4-11. A flow -chart for the EDK-SN Code System methodology. The function of th is module is to define a on which the net current vector will be discretized into. For each specific model, after executing the de terministic calculati on using PENTRAN, we generate 3-D fluxes, current vectors, and mate rial distributions for each energy group for a particular model. These 3-D da ta distributions, along with precomputed EDKs, will then be PAGE 95 95 implemented to calculate the 3-D absorbed dose distribution, for each energy group, scaled by the flux along the current vector direction. Finally the EDK-SN will accumulate these individual energy group dose contributions into a 3D total absorbed dose distribution. In Section 4.7, we describe the parallel code algorithm. 4.5 External Source Benchmark An external source problem for benchmarking purposes is presented to demonstrate the EDK-SN methodology and provide a dire ct comparison with a full -physics Monte Carlo solution. The model considered was a homogenous block of water (45 45 453cm) with a 15 15 1 3cm source (defined for 16 energy groups spanning from 0.0 to 8.0 MeV cast as a flat-weighted energy distribution source term). A 3-D spat ial mesh distribution was generated by the PENMSH-XP mesh generating code in the PENTRAN system19. The model was then subdivided into 3 z-levels; each z-level was then divided equally into 9 coarse meshes, each coarse mesh contained 3375 fine mesh cells; therefore, a to tal of 91,125 fine mesh cells were used in the complete SN problem as shown in Figure 4-12.A. A. B Figure 4-12. External Source impinging on block of water (left: Mesh di stribution generated by PENTRAN package; Right: MCNP5 Cell/Surface schematic input) PAGE 96 96 For the purposes of direct comparison to the Monte Carlo method, the MCNP5 Monte Carlo code is used to simulate an equivalent external sour ce problem as shown in Figure 4-12.B. For validation of the EDK-SN approach, the equivalent MCNP5 model is used to calculate the absorbed dose along the central axis and across the model for these studies using MCNP5 pulse height tallies (*F8), with full physics transport (p, e mode) to al low direct comparison of fluxes and absorbed dose values for the Cartesian fine mesh cells used in PENTRAN-MP and the EDKSN approach. 4.6 Comparison of Deterministic and Monte Carlo Results Using a uniform 1 3cm mesh, PENTRAN was used to calcu late the flux distributions and the net current for each fine mesh cell in th e system. Based on the PENTRAN results and the pre-calculated electron kernels, the EDK-SN code was used to produce 3-D absorbed dose rate distributions for all energy groups. The absorbed dose rates for energy groups (7.5-8.0 MeV) and (0.5-1.0 MeV) are depicted in Figures 4-13A and 4-13B, respectively. These figures provide a clear insight of the 3-D absorbed dose behavi or for these two different energy groups. For example, the high energy photons indicate more fo rward peaked absorbed dose distributions (as expected), and, while still forward peaked, the lower energy photons indicate more isotropic behavior, also as expected. Figures 4-14 and 4-15 depict the absorbed dose rate for each energy group (over a total of 16 energy groups) calculated by EDK-SN along the central axis, and laterally across the problem respectively. Alternatively, Figures 4-16 and 417 show the MCNP Monte Carlo total absorbed dose along the central axis and laterally ac ross the problem after 1000 minutes using 16 processors of same PC cluster. The statistical uncertainty 2 on average, of the *F8 MCNP5 Monte Carlo tally was 6.0% for the study along the central axis (x-axis), and 4.0% along the y- PAGE 97 97 axis, 2 cm away from the source, in the source region. On the other hand, the average absolute relative differences between the EDK-SN absorbed dose calculation and the full physics MCNP derived absorbed doses were 3.7% along the cent ral axis (X-axis), and 2.5% along the Y-axis, 2cm away from the source. A B Figure 4-13. External Source 3-D absorbed dos e distribution for 2 Energy groups computed by EDK-SN based on PENTRAN-MP A. EG(7.5-8.0Mev) and B. EG(0.5-1.0Mev) Figure 4-14. Energy group absorbed dose rate di stributions (in Mev/ g. sec) along the X-axis using the EDK-SN method based on PENTRAN-MP PAGE 98 98 Figure 4-15. Energy group absorbed dose rate distributions (Mev/ g.sec) along the y-axis using the EDK-SN method based on PENTRAN-MP Figure 4-16. EDK-SN Total absorbed dose rate distributions along the X-axis compared to MCNP *F8 tally derived absorbed doses. MCNP uncertainties (2 ) average (6.0%) in the Source region. The average absolute re lative difference betw een the two methods was 3.7% PAGE 99 99 Figure 4-17. Total absorbed dose rate distribu tions of EDK-SN along the Y-axis compared to MCNP *F8 tally. MCNP uncertainty (2 ) average (4.0%) in the Source region. Average absolute relative difference (2.5%) Although small, less than 1%, the ripple apparent in lateral absorbed dose, shown in Figure 4-15 and Figure 4-17, can be explained by the mesh structure used in projecting the absorbed doses along the collective region of elec trons, and low level ray effect, in the SN calculation, causing absorbed dose to oscillate 4.7 Parallelization of EDK-SN Code EDK-SN is designed to run on distri buted memory architectures, where each processor is an independent unit. This system is a distribute d memory computer system, where the processing elements are connected by a network and tare highly scalable. The interconnection scheme among the processors is fundamental for distribute d memory architectures because it affects, in part, the performance of the system. For clustertype architectures, the cluster using identical processors (matched for speed) connected with an Ethernet local area network to enable data transfer among the processor units. PAGE 100 100 To achieve maximum speedup with the EDK-SN methodology, para llel programming methods were implemented in the EDK-SN code using energy group decomposition via the Message Passing Interface (MPI)49 library. The MPI library is meant to provide an essential virtual topology, process sync hronization, and communication f unctionality between a set of processes (that have been mapped to nodes/server s/computers) in a language independent way, with library specific syntax (bindings), plus a few features that ar e library specific. The EDK-SN code has been designed in ANSI FORTRAN-90 (to take advantage of dynamic array allocation) and parallelized for operation on distributed memory, multiple instruction/ multiple data (MIMD) mach ine architectures using initial standard Message Passing Interface protocols (MPI-1). All Input/ Output (I /O) is performed by each processor in parallel (as required) to include input processing, initialization, processing, and file output. The EDK-SN code carries out the computations with electron dose kernels based on the scalar flux scaling for each energy group independe ntly Consequently, the decomposition used is energy domain decomposition. The basic philosophy of this approach is to decompose energy groups to each processor, and carry the total angular (net current vector) and spatial domain information on each processor. The main advantages of the EDK-SN parallel algorithm is that appreciable speed up can be obtained to determine the dose in any voxel throughout the spatial domain. Since little communication is required between the processors, n processors would solve the problem nearly n-times faster than a single unit. Energy group absorbed doses are partitioned for the entire problem, and partitioned energy group assignments are scaled to n processors for coupling photon transport information with EDKs to produce group absorbed doses over the spatial domain. An MPI_REDUCE call is made to collect all absorbed doses across processors. Since PAGE 101 101 this work focused on the use of 16 energy groups for a flat weighted 8 MeV high energy source term, the parallel EDK-SN computation could be pa rsed to a maximum of 16 parallel processors. Near linear speedups were achieved. The parallel code follows the following steps in executing the EDK-SN methodology for total dose calculation: Initiate the code on all pro cessors for sets of energy groups, up to one group per processor Sort the locally parallel-a llocated fine mesh distribu tions into a global monotone increasing mesh distributions Create Energy group dose matrix aliased to the location of the DDVs Translate material tags into densities for later calculation Project the Group-dependent E DK based on the group current v ector and material at the DDV Scale group dose using group scalar flux, correc ting for densities at material interfaces Use MPI-REDUCE to sum all doses for the phase space across the parallel processing grid 4.8 Discussion and Findings We examined a new EDK-SN methodology of coupling precomputed Monte Carlo based electron absorbed dose kernels (EDKs ) with the discrete ordinates (SN) solutions to achieve accurate absorbed doses for megavolt energy photon scenarios, since charged particle equilibrium fails and full phot on-electron physics must be acc ounted for to fully attribute absorbed dose. The EDK-SN methodology was compared with the Monte Carlo method and shown to give results consistent within the statistical uncertainty of the absorbed dose computed using Monte Carlo calculations. By following the EDK-SN approach, the absorbed dose applied to tissue from any incident external high energy pho ton beam (pencil, fan, areal, etc) can therefore be readily determined based on a coupling of the deterministic SN photon transport solution from the PENTRAN code using the EDK-SN methodology. As a result, this research rend ers the first PAGE 102 102 application of integrated electron transport kernels coupled to the SN photon transport methodology to properly account for the total abso rbed dose in high energy photon irradiation scenarios. A fundamental difference between the deterministica lly guided EDKSN technique and the existing convolutionsuperposition te chnique is the fact that our EDK-SN approach is based on a full physics Monte Carlo generated kernel projected using an accurate SN transport solution that provides a net curren t direction and flux scaling to sc ale the absorbed dose using the spectrally dependent EDKs. Overa ll, based on the water phantom problem examined here, with an adequate discreti zation using the SN method coupled to properly computed energy dependent electron kernels, it is po ssible to produce absorbed dose results with the EDK-SN method that are consistent with Monte Carlo calculations A major advantage of the EDK-SN method is that it provides a detailed, accurate absorbed dose distribution wherever required throughout the system, while incorporating the full physics electron transport, and is generally much faster than Monte Carlo. Since PENTRAN is based on a Cartesian geometry, CT voxel data can readily be mapped onto the problem geometry to define the tissue and organ geometry structure. Subsequent discussion will account for material density variations (encountered at tissue-bone interfaces) where we will apply the EDK-SN methodology to a whole body voxelized phantom to rapidly achieve a full physics whol e body, total body absorbed dose. PAGE 103 103 CHAPTER 5 APPLICATION OF EDK-SN TO MODELS WITH HETEROGENITIES This chapter focuses the many aspects that make up the EDK-SN methodology in order to apply it to a voxelized phant om, validated side-by-side using full physics Monte Carlo calculations. The SN calculations are performed based on a full photon transport calculation using deterministic techniques in a heterogene ous phantom. This eliminates any need for any photon corrections for tissue he terogeneities, since the SN solution yields the complete solution to the Boltzmann transport equation. Moreover, because the EDK-SN calculation is performed over the full spectrum of the radiation for a discrete number of ener gy groups, we can thereby generate separate sets of elect ron absorbed dose kernels to account for the different densities present in the phantom model. In this chapter, we will address applying the EDK-SN method to a heterogeneous phantom, including the application to a new state of the art human phantom. 5.1 Heterogeneous Tissue-Equivalent Phantoms Three slab based phantoms were simulated, using the EDK-SN method, as a prelude to more complex calculations for full human phantoms. Each phantom had a flat weighted 0 to 8 MV x-ray source impinging on 45x45x45 cc soft ti ssue (S) phantom. Two of these phantoms have slabs of 1-cm-thick discs of materials representative of bone (B), or lung tissue (L) (Figure 5-1). 5.1.1 Monte Carlo Simulations In order to compare the accuracy of the EDK-SN methodology, a series of Monte Carlo simulations were also performed using the radi ation transport code MC NP5. An X-ray flat weighted spectrum 8 MV used as input to MCNP5 to guarantee no specific preference was giving to any particular energy gr oup. Cell estimates of tissue absorbed dose were calculated in MCNP5 by tallying the energy depos ited by the photons and their seco ndary electron within 1 PAGE 104 104 cm3 voxels, equivalent to the size of the EDK-SN voxel and the voxel size to be used in the SN calculation. Figure 5-1. Three Monte Carlo phantoms used in the absorbed dose study. Phantom SSSSSSS is a homogeneous soft-tissue (S) phantom; phantoms SSBBSSS and SBLLBSS are heterogeneous phantoms consisting of combin ations of soft and bone tissue (B), and soft (S), bone (B), and lung (L) tissue-equivalent materials, respectively50 After running the three phantoms for 100 mi nutes, Figure 5-2, using MCNP5 photon, electron mode and pulse height tallie s (*F8) the results uncertainties (2 ) were 20% on average. A continue run was executed for the three phantoms for 400 minutes more, for 500 minutes total, reducing the Mont e Carlo uncertainty (2 ) down to 10% on average (Figure 5-3). However, the Monte Carlo computed difference in the depos ited dose between the two tissues is within 2%, making it impossible to observe a ny difference between the differe nt tissues in the phantoms, especially for therapeutic photon beams; this is shown in (Figure 5-4). The estimated time required to reduce Monte Carlo uncertainty to (2 ) 2% on-axis was 12500 minutes; hence, the time required to reduce the Monte Carlo uncertainty off-axis to 2% will be much greater. PAGE 105 105 Therefore, this is an example that highlights that tissue heterogeneities can be difficult to discern in these problems using the Monte Carlo method. As explained previously, one of the most im portant issues limiti ng the use of the Monte Carlo techniques for clinical absorbed dose calculat ions is achieving a statistically reliable result in a reasonable time. This problem is even gr eater when a photon-electron mode is required, when CPE is not applicable, causing computer running times needed to reduce the statistical error to increase rapidly. Improved variance re duction techniques may reduce the computation times. However, such techniques are to be us ed with caution, as improper application of statistical variance reduction may cause wr ong answers. Moreover, variance reduction techniques may not reduce the time required to reduce the variance to acceptable levels, especially when accurate absorbed dose esti mations are needed over the entire problem. Figure 5-2. Monte Carlo derived total absorbed dose rate distributions al ong the X-axis computed using *F8 MCNP5 tally for three different phantoms. MCNP5 uncertainty (2 ) average (20.0%), 100 Minute run PAGE 106 106 Figure 5-3. Monte Carlo derived total absorbed dose rate distributions al ong the X-axis computed using *F8 MCNP5 tally for thre e different phantoms. MCNP uncertainty (2 ) average (10.0%), 500 Minutes run Figure 5-4. Monte Carlo derived total absorbed dose rate distributions al ong the X-axis computed using *F8 MCNP5 tally for th e _SBLLBSS_ phantom. MCNP uncertainty (2 ) average (7.0%), 1000 Minutes run PAGE 107 107 5.1.2 PENTRAN Model As shown in Figure 5-5, we simulated equiva lent three phantoms to those simulated with Monte Carlo. The heterogeneous models consiste d of three key components; an X-ray source (15 cm3) with flat weighted 8 MV source, divided into 16 energy groups, a soft tissue phantom (45 cm3), equivalent layers of bone and lung tissue identical to that used in the Monte Carlo simulations. Models were divided into 3 z-levels; each z-level was divided equally into 9 coarse meshes containing (15) fine mesh cells for each coarse mesh; a total of 91,125 fine mesh cells were used in the SN problems. Although lower quadrature set (S32) are sufficient to produce accurate SN results for similar calculations (Chapter 1), to assure minimum ray-effects on the SN calculations, these simulations were based on the use of an S62 angular Legendre-Chebychev quadrature with P3 scattering anisotropy using the CEPXS library. Figure 5-5. Three SN phantoms equivalent to the one us ed in the Monte Carlo simulation. Phantom SSSSSSS is a homogeneous soft-tissue phantom; phantoms SSBBSSS and SBLLBSS are heterogeneous phantoms consisting of combinations of soft and bone tissue, and soft, bone, and lung tissueequivalent materials, respectively PAGE 108 108 5.1.3 EDK-SN Absorbed dose Comparison The SN calculation was performed based on a full photon transport calculation using deterministic techniques in the heterogeneous phantom. In principle, this eliminates any need for any correction for tissue heterogeneities which affect the different sc attering components in differing ways. Moreover, the EDK-SN calculation was performed for material specific absorbed doses. This allowed the use of the pre-computed electron absorbed dose kernels for the different densities present in the phantom model. Figure 5-6. EDK-SN heterogene ous SBLLBSS phantom total absorb ed dose rate distributions along the X-axis compared to MCNP *F8 tally. MCNP uncertainty (2 ) average (8.0%) on-axis. Average absolute relative difference (4.1%) For the phantoms used in this study, calculations were performed to compare the EDK-SN and Monte Carlo pulse heig ht tally results for the absorbed dose deposited along the model central axis. As shown in the Figures 5-6 and 5-7, the EDK-SN results are in very good agreement with the MCNP5 results, and differen ces between the two are within the statistical uncertainty of the Monte Carlo ca lculation. The Monte Carlo uncerta inties, on the central axis of PAGE 109 109 the source, were 4% on average while the aver age absolute difference at the same cells was 3.7%. Figure 5-7. EDK-SN Ho mogenous phantom total absorbed dos e rate distributions along the Xaxis compared to MCNP *F8 tally. MCNP uncertainty (2 ) average (6%) on-axis. Average absolute relative difference (3.7%) We note that at the interface l ung region, Figure 5-6, the EDK-SN derived absorbed dose indicated a lower absorbed dose th an Monte Carlo. This can be related to the application of the lung electron kernel over the bone and soft tissue regions. Th is phenomenon was mitigated by applying a correction factor for voxe ls located at the interface of different material densities. 5.1.4 Correction Factor Approach The range of a charged partic le of a given type and ener gy in a given medium is the expectation value of the path-length that it follow s until it comes to rest On the other hand, the RCSDA represents the range in the continuous slowing down approximation (CSDA). Figure 5-8 gives the CSDA range RCSDA for electron in cortical bone, soft tissue and lung tissue for electron energies between 0.0 MeV to 17.5 MeV. Figure 5-9 depicts a comparison between the ratio of PAGE 110 110 the CSDA ranges for cortical bone to soft tissue, lung to soft tissue, and cortical bone to tissue. One can notice that the behavior of electrons in the different tissues is quite similar for the same electron energy. This implies that there is a possibility of applying a simple correction factor within different tissues to account for the different densities in the phantom for each energy group. Unlike other methodologies, this is possibl e due to the dependence of the EDK-SN on a full photon physics transport calculation accomplished by the SN computation to provide the scattered photon information and a directional (current) component. Figure 5-8. Continuous slowing down approximation range (Rcsda ) for electrons in Cortical Bone, Soft Tissue and Lung as a function of energy PAGE 111 111 Figure 5-9. Continuous slowi ng down approximation range (RCSDA) ratio comparison for electrons in Cortical B one, Soft Tissue and Lung as a function of energy The key to this approach is that photon transport is carried with the 3-DSN code PENTRAN. The use of such code provided us with the directional component of the photon transport in the homogenous medium. By having this information we were able to project the 3D electron kernel along the net cu rrent direction. The assumption of a homogenous medium in photon transport in a heterogene ous problem, often used in other methods, may lead to an underestimation of the scattered component of the photon beam, and misrepresentation of its directional component, especially between the different organs of the human body (Figure 5-9), leading to misrepresentation of the absorbed dose. 5.2 Computational Phantom of 15-Year Old Male A new set of stylized phantoms have been developed by the UF Advanced Laboratory for Radiation Dosimetry Studies (A LRADS). These phantoms represent a totally new class of whole-body patient hybrid phantom s so named as they retain th e non-uniform scalability of PAGE 112 112 stylized phantoms, while mainta ining the anatomic realism of segmented voxel phantoms with respect to organ shape, depth, and inter-organ positioning to great accuracy. One study showed that these phantoms match (to within 1%) to that of a newborn of the ICRP Publication 89 reference newborn child. In this study, hybrid computational phantoms represen ting reference 15-year male body anatomy and anthropometry is used, shown in Figure 5-10. This phantom was matched, by the ALRADS group, to anthropometric data and re ference organ mass data provided by the ICRP, where it matched to within 1% overall51. Figure 5-10. Frontal views of 3-D rendering of the 50th percen tile UFH-NURBS 15-year male and female The availability of hybrid phantoms as a re placement for voxel-based phantoms provides the needed platform for constructing phantoms representing the unique internal and external PAGE 113 113 anatomy of individual pediatri c patients. However, this process would be complicated tremendously if each individual organ of the phantom had to be adjusted to match patientspecific organ volumes as seen under CT imagin g, even when those images are available. The 15-year-old male was voxelized with hi gh resolution using 2 mm 2 mm 2 mm voxel size (302 voxel 139 voxel 836 voxel) as a case study. This geometry data was then collapsed into a 1 cm 1 cm 1 cm voxe l size (60 27 167 voxels), for total of 270540 voxel-model using GHOST-3-D. Subsequently, using PENMSH-XP, the collapsed model was cast into a 3-D spatial distributi on, and sub-divided into 10 coarse mesh z-levels containing 18 coarse meshes in x-y domain, and these contained a number of fine mesh cells (Figure 5-11). An 8 MV flat weighted, forward peaked volumetric (20 cm x 1 cm x 17 cm) X-ray source was placed above the chest area of the phantom. Figure 5-11. Simulati on methodology for EDK-SN computations using PENTRAN-MP code system. Left: UF Anthropmorphic phantom, middle: corresponda nt PENTRAN input, right: EDK-SN 3-D absorbed dose distribution after irradiation with a x-ray source PAGE 114 114 Although the model has relatively small fine mesh sizes, using Legendre-Chebychev S32 quadrature, the DTW differencing scheme, and CEPXS cross-sections with anisotropy level P3, we were able to eliminate any unphysical oscillations that may result in the deterministic solution. Using 16 energy groups, the SN model was executed on 16 processors and completed in 2.0 hours. With the general purpose radiation transport code MCNP5, we simulated an equivalent model to the one simulated using PENTRAN. The collapsed 3-D binary fi les were converted to ASCII formatted matrix files and then incorporat ed into an MCNP5 lat tice defined input; this matrix was also generated by our down-sampling code GHOST-3-D. 5.2.1 EDK-SN Absorbed Dose Calculations Calculations were performed to calculate the SN results for the scalar fluxes for all energy groups. Using EDK-SN, 3-D absorbed dose distributions we re calculated to each voxel in the entire 15 year old male phantom model for a broad-beam external photon field. Figure 5-11 shows the stages of calculation in the PENT RANP-MP code system that resulted in 3-D Absorbed dose distributions for the entire phantom. Figure 5-12 shows absorbed dose contributions from selected en ergy groups, and the total absorbed dose resulting from adding all energy group-absorbed dose contribu tions, as prescribed by the EDK-SN procedure. Organ absorbed doses using EDK-SN, for a selected number of or gans, were calculated and compared to Monte Carlo predictions, and are pr esented in Table 5-1. The total execution time for the Monte Carlo was ~4 hrs (wall clock time) on 16 processors, while the total wall clock execution time for the EDK-SN was ~ 2hrs on the same number of processors. As anticipated, the absorbed doses for the organs loca ted directly in the x-ray field (r ight+ left lung) or partially in the x-ray field (liver) were in good agreement for both the EDK-SN and MCNP5 pulse height tally data. PAGE 115 115 Total-Dose EG1 EG13 EG10 EG7 EG4 EG16 Figure 5-12. Absorbed dose contribution of eac h energy group is calculated via the EDK-SN procedure, then folded to calculate the to tal absorbed dose deli vered to the 15 year old male phantom from 8 MeV flat wei ghted source. The original phantom was collapsed into 270K voxels Table 5-1. Comparison of sel ected organ absorbed dose rate (MeV/g. s) calculated using MCNP5 pulse height tally with (photon, electron mode) and EDK-SN for the UF hybrid 15-year-old male phantoms for a ch est 8 MV X-ray flat weighted source PAGE 116 116 On the other hand, Monte Carlo calculations for organs outside the x-ra y field (such as the prostate) showed significant statistical uncertainty (44%). This poor absorbed dose tally result is due to the fact that the Monte Carlo did not have enough histories for these results to converge to low statistical uncertainty in an acceptable time. Alternatively, as shown in Table 5-2,when running the Monte Carl o for extended running times (~40 hours on our parallel cl uster), the Monte Carlo calculated absorbed dose delivered to the prostate began to show convergence to the absorbed dose calculated using EDK-SN. This illustrates that the EDK-SN methodology, executed in 2 hours (SN execution time = 1.5 hr, EDKSN execution time = 0.5 hr), is capable of providi ng an accurate solution and a significant speed up for a full physics absorbed dose calculation in human phantoms for high energy external beams. Table 5-2. Prostate absorbed dose rate (MeV/g .s) calculated for the test problem using EDK-SN and MCNP5 using pulse height tally (*F8) with (photon, electron mode) at different MCNP5 Monte Carlo running times 5.2.2 Ray-Effects When considering a forward peaked volumet ric (10 cm x 1 cm x 17 cm) X-ray source placed on top of the left lung area of the phant om, we noted significant differences between the absorbed dose calculated by EDK-SN and the absorbed dose ca lculated by full physics Monte PAGE 117 117 Carlo calculations to the right lung. This behavior can be explained by the ray-effects in the low density lung tissue (see Chapter 2). Since the X -ray source is highly directional, and the organ (right lung) is of a low density material, and little primary radiation and very low scattered radiation is streaming through the right lung, making conditions ideal for inducing ray effects in SN computations for the photon component. This is mitigated by using high order quadratures, and can reduce absolute differences; these were reduced to 4% (Table 5-2). However, the increase of the quadrature level costs additional running time, with an increase from 2 hours to 9 hours. This increase in quadrature level was necessary for large low density volu mes that experience low levels of scattered radiation. Another method, ordinate splitting, a llows for a biased quadrature set to focus on the important directions, although the technique was not used here. Table 5-2. Absorbed dose rate (MeV/g. s) calculated for the test problem using EDK-SN and MCNP5 using pulse height tallies (* F8) with (photon, electron mode) 5.2.3 MCNP5 Energy Deposited Tally (F6) and Pulse Height Tally (*F8) Results When comparing MCNP5 organ absorbed dos e results, for the 15 -year old phantom, calculated using energy deposited tally (F6) and pulse heig ht tally (*F8), we noted inconsistencies in the results of the two tallies (T able 5-3). This is in keeping with our earlier investigation of the two tallies earlier in Chapter 4. This in consistency in the F6 tally suggests that full physics Monte Carlo ca lculations are essential for accu rate absorbed dose estimation in high energy external beam absorbed dose calculations. PAGE 118 118 Table 5-3. Absorbed dose rate (MeV/g. s) cal culated for the test pr oblem using MCNP5 using pulse height tally *F8 wi th (photon, electron mode) and MCNP5 energy deposited tallies (F6) Wang et al. conducted a similar study using th eir VIP Man organ MCNP for absorbed dose comparison for the three different tally methods including cell flux tally (F4), energy deposited tally (F6) and pulse heig ht tally (*F8). In thei r comparison they found that the F4 tallies agree within 10% with the F6 tallies, while the differenc es between the tallies F8 and F6 are relatively large (28%) 14. These results agree with our findings, a nd provide a strong proof for the need of applying full physics transport when cal culating whole body absorbed dose. 5.3 Discussion and Findings To properly treat 3-D electron transport physic s deterministically, yet still to achieve reasonably fast and accurate whole body computa tion times using high energy photons, we used angular and energy dependent electron transport absorbed dos e kernels. In a methodology called EDK-SN we have demonstrated that by using different electron kernels for different materials (cortical bone, soft tissue and lung), the EDK-SN methodology incorporating a discrete ordinates (SN) approach is an effective alternativ e to the Monte Carlo method for accurately estimating the absorbed dose for high energy photon beams over an entire human phantom. In doing so, we simulated three phantoms using MCNP5 and PENTRAN-MP. Each phantom had a flat weighted 8 MV x-ray s ource impinging on 45x45x45 cc soft tissue phantom. Two of these phantoms were slabs of 1-cm-thick discs representative of bone, or lung tissue. Comparisons of MCNP5 Monte Carlo based dose results a nd those computed using PENTRAN-MP implementing the EDK-SN methodology showed good agreement between the two approaches. PAGE 119 119 However, the Monte Carlo simulation took 100 minutes for the MC (2 ) uncertainty to be 20% (on average) for tally points along the central axis of the phantom using an MCNP5 *F8 tally (we note here that 20% uncertain ty cant be viewed as converged). Additional runs totaling 500 minutes were made, reducing the Monte Carlo (2 ) uncertainty to 10%, on average. We estimated the approximate time required to reduce Monte Carol (2 ) uncertainty to 2% for the on-axis locations (relative to the source X-ray beam) to be 12,500 minutes. Undoubtedly, the time required to reduce the Monte Carlo uncertainties at locations off-axis from the source beam to 2% will be much greater. Based on the UF hybrid phantoms developed by the University of Florida ALRADS group, we simulated a complete phantom discretized using 270,540 voxels, and demonstrated that; in this case, the EDK-SN method readily provides accurate whole body 3-D absorbed dose distributions. We were able to accurately yield doses analogous to the MCNP5 Monte Carlo calculations. We further noted that ray-effects were evident in regions of the tissue-lung interface, which indicate higher angular quadratures or ordinat e splitting treatments are necessary in lung tissue and low density tissue inte rfaces. Also, tally results from the 16 energy groups, used to represent the radiation spectru m from the CEPXS multi -group cross-section library were compared to continuous energy MCNP5 Monte Carlo results using the ENDF/B-VI data libraries. A major advantage of the EDK-SN method is that it provid es a detailed, accurate 3-D absorbed dose distribution at thousands of ce lls throughout the system, while the Monte Carlo method only provides highly accurate values for sele cted points near the source. Though the use of volumetric mesh tallies, Monte Carlo calcula tions can provide global absorbed doses across a grid spanning the Monte Carlo geometry. However, at present, absorbed doses for this data are PAGE 120 120 based on assuming charged particle equilibrium conditions, not suitable for high energy beams. Moreover, the statistical uncertainty associated w ith each mesh tally site for such computations would require extended periods of running time, es pecially off-axis from the source and outside the region covered by the primary beam to converg e to acceptable levels. These factors limit the Monte Carlo calculations from bei ng fuuly adequate in practical applications where whole body and organ-specific absorbed dose information is needed. Notably, an equivalent Monte Carlo calculation to obtain the flux distri bution with statistically reliable results over the entire problem space may require orders of magnitude longer running times compared to parallel deterministic calculations. PAGE 121 121 CHAPTER 6 SUMMARY, CONCLUSIONS AND FUTURE WORK This work began with the prim ary goal of developing a new parallel 3-D Dose calculation tool based on deterministically rendered high energy photon transpor t solutions in radiotherapy settings. 6.1 Conclusions and Findings In this dissertation, the EDK-SN method is presented for generating detailed absorbed dose information using high energy external beam so urce terms. Because detailed whole body and organ-specific absorbed doses can be accurate ly obtained throughout a voxelized phantom with proper electron interaction physics projected onto the phantom voxels, based on SN-computed photon fluxes and net current vectors, this highlights the strength of the EDK-SN transport based methodology. Overall, based on the problems ex amined, with an adequate discretization implemented using the EDK-SN method, it is possible to produce absorbed dose results consistent within the stat istical uncertainty of MCNP 5 Monte Carlo calculations. Overall, our studies examined the capabilitie s of the discrete ordi nates method compared with the Monte Carlo method with application to medical physic s problems. Since radiation transporting in such applications encounters few interacti ons in low density mate rial, such as air, attention was given to minimum quadrature levels and mesh size distributions necessary for solution accuracy in air and in a water phantom. Also, it is eviden t that multi-gr oup cross-section libraries must also be carefully tested for probl em applicability. Of particular importance was the flexibility of the adaptive SN differencing scheme in PENTRAN, permitting more accurate spatial differencing schemes in cells with a range of optical thicknesses. In addition, deterministic models based on the UF series B pediatric phantoms and with the PENTRAN-MP PAGE 122 122 methodology, we were able to accurately yiel d a solution analogous to the MCNP5 Monte Carlo calculations using the CPE principle for low energy beams. We have demonstrated; that by using electron kernels for th e different materials (cortical bone, soft tissue and lung) the EDK-SN methodology incorporating a discrete ordinates (SN) approach is an effective alte rnative to the Monte Carlo met hod for accurately estimating the absorbed dose for high energy photon beams over an entire human phantom. In doing so, we simulated three phantoms using MCNP5 and PENT RAN. Each phantom ha d a flat weighted 8 MV x-ray source impinging on a 45x45x45 cc soft ti ssue phantom. Two of these phantoms have slabs of 1-cm-thick discs of materials represen tative of bone, or lung ti ssue. The comparison of the MCNP5 results and the EDK-SN showed good agreement between the two methodologies. We also simulated the complete human pha ntom using 270,540 voxels, and demonstrated that the EDK-SN method readily provides accurate whol e body 3-D absorbed dose distributions. The angular data provided by the SN method was used to guide the absorbed dose mapping process for secondary electrons to the surrounding voxels ; this makes our approach truly novel. We were able to accurately yield a solution analogous to the MCNP5 Monte Carlo calculations in the voxelized phantom. We note that that tally results from the 16 energy groups, used to represent the radiati on spectrum, from the CEPXS multi-g roup cross-section library, were compared to continuous energy MCNP5 Monte Carlo results using the ENDF/B-VI data libraries. Parallel deterministic PENTRAN-MP results were obtained within comparable running times to parallel Monte Carlo calculations for ta lly sites adjacent to the source. An equivalent Monte Carlo calculation to obtain the absorbed dose distributions using photon, electron mode with pulse height tallies converg ed to reasonable sta tistical uncertainty ove r the entire problem PAGE 123 123 space may require orders of magnitude longer running times compared to the time needed for complete, whole body parallel deterministic calculat ions. The work presente d in this dissertation has extended the capabilities of the discrete ordinates method to accurately calculate radiation absorbed dose contributions in medical phys ics procedures using the PENTRAN-MP code system with the EDK-SN methodology. 6.2 Future Work The new EDK-SN, integrated into the PENTRAN-MP Code System, represents a leap forward in computational medical physics; large 3-D radiation transport calculations for high energy photons and absorbed dose contributions can now be performed accurately in a fraction of the computation time requi red in the Monte Carlo methodol ogy. The methods described in this dissertation can be further enhanced and developed by studying the following issues: Additional parallel performance assessments. Examine performance for solving additional real-world problems. Optimization study on the EDK-SN. This should be undertaken, in order to identify the necessary level of accuracy wh ich yields the best performance in terms of speed-up. Test alternative method for applying corre ction factors in tissu e inhomogeneities. Conduct timing and performance studies with the new hybrid PENTRAN code implementing Discrete Ordinates and Me thod of Characteristics. Apply for an NIH R01 grant. Long-term plans are to implement PENTRA N-MP as a secondary treatment planning software for documentation of both in-field an d out-of-field organ absorbed doses to normal tissues. This information will be extremely valuab le to radiation epidemiology studies which are currently being conducted using far more simplisti c dosimetry schemes, particularly to out of fields organs. This data will also be of use to future optimization studies where present goals of tumor absorbed dose conformality and near-field ab sorbed dose avoidance to normal tissues will be supplemented by the additiona l goal of secondary cancer risk minimization in out-of-field tissues and organs. PAGE 124 124 Att. Coefficient WATER (cm2/g) 1.00E-08 1.00E-06 1.00E-04 1.00E-02 1.00E+00 1.00E+02 1.00E+04 1.00E-031.00E-021.00E-011.00E+001.00E+011.00E+02 Energy(MeV)Att. Coefficient(cm2/g) Scatt. Coh.(cm2/g) Scatt. Incoh.(cm2/g) P.E.A(cm2/g) P.P Nuc. Field (cm2/g) P.P Elec. Field (cm2/g) Total w/ coh Total w/o Coh. APPENDIX A HIGH ENERGY PHOTON TREATMENT (IN CONDENSED ELECTRON TRANSPORT) The photons energy used in radiotherapy tr eatm ent planning are 20 MeV or less. Photon interactions that take place at these high energies predomin antly Compton Effect (C.E.) or/and Pair Production interactions by different percenta ges, as shown in Figure A-1 and Table A-1. Figure A-1. Total mass attenuation coefficient for photons in wate r, indicating the contributions associated with phoelectric absorption, Compton scattering and electron-positron pair production Table A-1. Ratios of photon interacti ons in water for range of energies Photon Energy (MeV) Photon Interaction 2.0 3.0 6.0 10.0 14.0 18.0 20.0 % of Incoherent scattering 99.21 % 97.15 % 88.59 % 77.06 % 67.32 % 59.44 % 56.04 % % of Nuclear P.P 0.79 % 2.82 % 10.78 % 21.18 % 29.76 % 36.69 % 39.64 % % of electron P.P 0.0 % 0.03 % 0.61 % 1.76 % 2.9 % 3.90 % 4.35 % PAGE 125 125 Compton Scattering: When passing through matter, photons can interact with the electrons of the matter. The Compton scattering describes an inelastic scattering of a photon w ith a free electron. If the electrons initial bindi ng energy is small compared to the photon energy (<< h ) then the stuck electron can be considered initially free and unbound. Figure 10 shows this inelastic scattering process. The kinematics of this sc attering (conservation of momentum and energy) lead to the energy of the scattered photon 'h and to the energy of the scattered electron Figure A-2. Sketch of Compt on scattering, photon of energy h collide with an electron. The photons incident forward momentum is hc where c is the speed of light in vacuum Based on the unbound and stationary condition of the electron, the Compton relation can be obtained by solving energy conservation equations 'Thh (A-1) And conservation of momentum PAGE 126 126 'coscosehh p cc (A-2) 'sinsinehpc (A-3) p c can be written in terms of T by invoking the law of invariance: 2 0(2) p cTTmc (A-4) In which 0m is the electrons rest mass. As a result of solving these equations algebraically, one can derive the following Comptons relations 2 01()(1cos) h h hmc (A-5) 2 0cot(1)tan 2ehmc (A-6) 21 1 12sin/2eThhh (A-7) To obtain the fraction of the scattered photons in a given direction, Klein and Nishina (reference) have carried out a quantum-mechanical treatment of the problem using the Dirac equation for the electron and have obtained the equation 2 2 '' 2 0 'sin 2KNdr hhh dhh h (A-8) Equations similar to those for the Compt on photon distribution can be obtained for the Compton electron distribution. Sinc e the probability that an electron will be scattered into the solid angle 'd situated in the direction e is the same as the probability that a primary quantum PAGE 127 127 will be scattered into the solid angle d in the direction and being related by Eq. (A8), we have that KNKN eedd d ddd (A-9) sin sineeeeddd ddd (A-10) cot(1)tan() 2e (A-11) 2211 (1) sincos2ed (A-12) 1/2 22sin1(1)tan2e (A-13) 2e (A-14) 1e ed dd d (A-15) 31(1cos)sin (1) sine ed d (A-16) for 0.511 MeV photons,1 3/2 22sin(/2)13sin(/2)ed d (A-17) 2 2 '' 3/2 22 0 '2sin(/2)13sin(/2) sin 2KN edr hhh dhh h (A-18) PAGE 128 128 Figure A-3. Compton effects differential cross section per unite angle for the number of electrons scattered in the direction e Example Considering a Compton scattering of 1.022 MeV (2 ), the ratio of recoil electrons with energies above 0.511 MeV to recoil elect ron with energies below 0.511 will be: '(1) 0.5 1(1) Thh hh (A-19) 0.560o (A-20) cot(1)tan 2e (A-21) ()30o e 180 60 60 02sin 2sinKN KNd d d R d d d (A-22) PAGE 129 129 2 22 222 0(22)/(1)(12) 2KNdr d (A-23) 1 1(1) (A-24) cossin dd (A-25) dd (A-26) 'h h (A-2 7) From this example one can note that; the prod uced electron in a cer tain direction has a unique energy for a giving initial photon. So for an angle aroundiscretized range of angles we can calculate the average energy of the electrons produced in th Pair production Figure illustrates schematically a pair production event in the nucleus electric field. The incident photon gives up all of it s energy in the creation of an el ectron positron pair with kinetic energies TandT. The energy conservation equation, ignoring the vanishing small kinetic energy giving to the nucleus is simply 2 02 1.022 hmcTT M eVTT (A-28) The average electron/positron average energy is 1.022 2 hMeV T (A-29) PAGE 130 130 Figure A-4. Pair Production in the Coulomb force field of the nucleus For h values well above th e threshold energy 22 0mc, the electron and positron are strongly forward directed. Their average angle of departure relative to the original photon direction is roughly: 2 0mc T (radians) (A-30) Pair production in the electron field: In the kinematics of Pair Production in the electron filed (triple production), the photon divides its energy between the positron-electron pair produced and the host electron. The energy conservation equation becomes: 2 012 122 1.022 hmcTTT M eVTTT (A-31) The average kinetic energy of the three particles is 1.022 3 hMeV T (A-32) PAGE 131 131 In triple production, the threshold for this process is 2 04 mc=2.044 MeV. The higher threshold is required by conservation of momentum as shown in Table A-1 and also derived by Perrin (1933). Triple production becomes signifi cant in radiotherapy treatment planning for energies greater than 14 MeV were it repres ents ~3% of the total photon interaction. Straggling and Multiple Scattering There is typically a distribution of angle and energy for charged particle passing the same medium. This is a result of multip le scattering and angular straggling which the spread in angular distribution observed in a population of initial identical ch arged particle after they have traverse equal path lengths. It will be some what exaggerated if the particle have passed through a layer of material, since multiple scattering then cause indi vidual differences in scattering angle as well causing them to propagate in a conical angular distribution in accordance with Molieres theory. However, for a discretized solid angle, one can consider an average angle for these scattered electrons on which approximate transport calculations can be made. The spread of the solid angle can be giving by Fermi angle for stra ggling especially when considering restricted stopping power including soft collisions, that yield small spread in the angle, and hard collisions resulting in rays with energies less than the cutoff value Restricted stopping power (/) dTdx Is the fraction of the collisions stopping power that include the soft collisions plus those hard collisions resulting in rays with energies less than the cutoff value 2 2 0(2) 2 ln (,) 2 dT C kG dx Z Imc (A-33) PAGE 132 132 Radiative losses: (/) (/)820rad ionizingdTdx Z T dTdx (A-34) Angular straggling (Fermi 1950) Electrons suffer multiple scattering events preclude straight transport in the material. Further complication are caused by the statistical fluctuation Figure A-5. Pair Production in the Coulomb force field of the nucleus 22 22 0 24322() ln 2 Ta zee NZ Z e T (A-35) 2 0 2 0h a mc (A-36) PAGE 133 133 APPENDIX B MCNP5 PHNTOM INPUT DECK c UF_15yr voxel MCNP model whole body c Matrix size [60,27,167] C Voxel resolution=1.0 1.0 1.0 c 2008-7 c Ahmad Al-Basheer c The University of Florida C ************************************************************************** c c cells Card c C ************************************************************************** c ---------------------------c Body composition and density c ---------------------------1 1 -0.001205 -999 imp:p,e=1 u=1 vol=78406.867 $ Air 2 2 -1.03 -999 imp:p,e=1 u=2 vol=3.94E+04 $Residual Soft Tissue 3 3 -1.03 -999 imp:p,e=1 u=3 vol=2.00E+00 $Adrenal (L) 4 4 -1.03 -999 imp:p,e=1 u=4 vol=3.00E+00 $Adrenal (R) 5 5 -1.04 -999 imp:p,e=1 u=5 vol=1.39E+03 $Brain 6 6 -0.94 -999 imp:p,e=1 u=6 vol=1.70E+01 $Breast 7 7 -1.07 -999 imp:p,e=1 u=7 vol=4.00E+00 $Bronchi 8 8 -1.03 -999 imp:p,e=1 u=8 vol=9.70E+01 $Right Colon W 9 9 -0.83 -999 imp:p,e=1 u=9 vol=1.45E+02 $Right Colon C 10 10 -1.10 -999 imp:p,e=1 u=10 vol=3.00E+00 $Ears 11 11 -1.03 -999 imp:p,e=1 u=11 vol=2.40E+01 $Esophagus 12 12 -1.05 -999 imp:p,e=1 u=12 vol=2.00E+00 $External nose 13 13 -1.07 -999 imp:p,e=1 u=13 vol=1.10E+01 $Eye balls 14 14 -1.03 -999 imp:p,e=1 u=14 vol=0.00E+00 $Gall Bladder W 15 15 -1.03 -999 imp:p,e=1 u=15 vol=4.50E+01 $Gall Bladder C 16 16 -1.04 -999 imp:p,e=1 u=16 vol=1.86E+02 $Heart W 17 17 -1.06 -999 imp:p,e=1 u=17 vol=4.23E+02 $Heart C 18 18 -1.05 -999 imp:p,e=1 u=18 vol=8.00E+01 $Kidney-cortex (L) 19 19 -1.05 -999 imp:p,e=1 u=19 vol=8.40E+01 $Kidney-cortex (R) 20 20 -1.04 -999 imp:p,e=1 u=20 vol=3.10E+01 $Kidney-medulla (L) 21 21 -1.04 -999 imp:p,e=1 u=21 vol=3.20E+01 $Kidney-medulla (R) 22 22 -1.04 -999 imp:p,e=1 u=22 vol=5.00E+00 $Kidney-pelvis (L) 23 23 -1.04 -999 imp:p,e=1 u=23 vol=6.00E+00 $Kidney-pelvis (R) 24 24 -1.07 -999 imp:p,e=1 u=24 vol=2.20E+01 $Larynx 25 25 -1.07 -999 imp:p,e=1 u=25 vol=0.00E+00 $Lens 26 26 -1.06 -999 imp:p,e=1 u=26 vol=1.23E+03 $Liver 27 27 -0.240 -999 imp:p,e=1 u=27 vol=1.72E+03 $Lung (L) 28 28 -0.240 -999 imp:p,e=1 u=28 vol=2.04E+03 $Lung (R) 29 29 -1.02 -999 imp:p,e=1 u=29 vol=0.00E+00 $Nasal layer (anterior) 30 30 -1.02 -999 imp:p,e=1 u=30 vol=0.00E+00 $Nasal layer (posterior) 31 31 -1.02 -999 imp:p,e=1 u=31 vol=0.00E+00 $Oral cavity layer 32 32 -1.02 -999 imp:p,e=1 u=32 vol=1.08E+02 $Pancreas 33 33 -1.05 -999 imp:p,e=1 u=33 vol=2.30E+01 $Penis 34 34 -1.03 -999 imp:p,e=1 u=34 vol=2.00E+00 $Pharynx 35 35 -1.03 -999 imp:p,e=1 u=35 vol=0.00E+00 $Pituitary Gland 36 36 -1.03 -999 imp:p,e=1 u=36 vol=4.00E+00 $Prostate 37 37 -1.03 -999 imp:p,e=1 u=37 vol=0.00E+00 $Rectosigmoid W 38 38 -0.49 -999 imp:p,e=1 u=38 vol=1.39E+02 $Rectosigmoid C PAGE 134 134 39 39 -1.03 -999 imp:p,e=1 u=39 vol=3.40E+01 $Salivary Glands (parotid) 40 40 -1.03 -999 imp:p,e=1 u=40 vol=1.00E+01 $Scrotum 41 41 -1.03 -999 imp:p,e=1 u=41 vol=1.40E+02 $SI W 42 42 -0.41 -999 imp:p,e=1 u=42 vol=7.70E+02 $SI C 43 43 -1.09 -999 imp:p,e=1 u=43 vol=1.72E+02 $Skin 44 44 -1.04 -999 imp:p,e=1 u=44 vol=4.10E+01 $Spinal Cord 45 45 -1.06 -999 imp:p,e=1 u=45 vol=1.21E+02 $Spleen 46 46 -1.03 -999 imp:p,e=1 u=46 vol=3.30E+01 $Stomach W 47 47 -1.02 -999 imp:p,e=1 u=47 vol=2.08E+02 $Stomach C 48 48 -1.04 -999 imp:p,e=1 u=48 vol=1.80E+01 $Testes 49 49 -1.03 -999 imp:p,e=1 u=49 vol=3.20E+01 $Thymus 50 50 -1.05 -999 imp:p,e=1 u=50 vol=8.00E+00 $Thyroid 51 51 -1.05 -999 imp:p,e=1 u=51 vol=5.50E+01 $Tongue 52 52 -1.02 -999 imp:p,e=1 u=52 vol=2.00E+00 $Tonsil 53 53 -1.03 -999 imp:p,e=1 u=53 vol=6.00E+00 $Trachea 54 54 -1.04 -999 imp:p,e=1 u=54 vol=3.00E+00 $Urinary bladder W 55 55 -1.00 -999 imp:p,e=1 u=55 vol=1.61E+02 $Urinary bladder C 56 56 -0.001 -999 imp:p,e=1 u=56 vol=4.23E+02 $Air (in body) 57 57 -1.03 -999 imp:p,e=1 u=57 vol=7.00E+01 $Left Colon W 58 58 -0.34 -999 imp:p,e=1 u=58 vol=1.95E+02 $Left Colon C 59 59 -1.03 -999 imp:p,e=1 u=59 vol=2.00E+01 $Salivary Glands (submaxillary) 60 60 -1.03 -999 imp:p,e=1 u=60 vol=7.00E+00 $Salivary Glands (sublingual) 61 100 -1.10 -999 imp:p,e=1 u=61 vol=2.20E+01 $Coastal cartilage of ribs 62 100 -1.10 -999 imp:p,e=1 u=62 vol=0.00E+00 $Cervical Discs 63 100 -1.10 -999 imp:p,e=1 u=63 vol=2.80E+01 $Thoracic Discs 64 100 -1.10 -999 imp:p,e=1 u=64 vol=1.20E+01 $Lumbar Discs 65 101 -1.519 -999 imp:p,e=1 u=65 vol=7.48E+02 $Cranium 66 101 -1.519 -999 imp:p,e=1 u=66 vol=3.30E+01 $Mandible 67 104 -1.387 -999 imp:p,e=1 u=67 vol=1.06E+02 $Scapulae 68 103 -1.351 -999 imp:p,e=1 u=68 vol=3.00E+01 $Clavicles 69 108 -1.219 -999 imp:p,e=1 u=69 vol=4.40E+01 $Sternum 70 103 -1.351 -999 imp:p,e=1 u=70 vol=2.90E+01 $Ribs 71 102 -1.196 -999 imp:p,e=1 u=71 vol=5.10E+01 $Vertebrae-C 72 102 -1.196 -999 imp:p,e=1 u=72 vol=2.28E+02 $Vertebrae-T 73 102 -1.196 -999 imp:p,e=1 u=73 vol=1.95E+02 $Vertebrae-L 74 108 -1.219 -999 imp:p,e=1 u=74 vol=1.29E+02 $Sacrum 75 103 -1.351 -999 imp:p,e=1 u=75 vol=4.67E+02 $Os Coxae 76 103 -1.319 -999 imp:p,e=1 u=76 vol=2.20E+02 $Femur-proximal 77 103 -1.319 -999 imp:p,e=1 u=77 vol=1.63E+02 $Femur-upper shaft 78 104 -1.387 -999 imp:p,e=1 u=78 vol=1.33E+02 $Femur-lower shaft 79 103 -1.319 -999 imp:p,e=1 u=79 vol=2.59E+02 $Femur-distal 80 106 -1.352 -999 imp:p,e=1 u=80 vol=2.07E+02 $Tibiae-proximal 81 106 -1.352 -999 imp:p,e=1 u=81 vol=1.01E+02 $Tibiae-upper shaft 82 106 -1.352 -999 imp:p,e=1 u=82 vol=1.01E+02 $Tibiae-lower shaft 83 106 -1.352 -999 imp:p,e=1 u=83 vol=7.30E+01 $Tibiae-distal 84 107 -1.273 -999 imp:p,e=1 u=84 vol=1.40E+01 $Fibulae-proximal 85 107 -1.273 -999 imp:p,e=1 u=85 vol=4.00E+00 $Fibulae-upper shaft 86 107 -1.273 -999 imp:p,e=1 u=86 vol=1.40E+01 $Fibulae-lower shaft 87 107 -1.273 -999 imp:p,e=1 u=87 vol=1.70E+01 $Fibulae-distal 88 106 -1.35 -999 imp:p,e=1 u=88 vol=2.40E+01 $Patellae 89 105 -1.415 -999 imp:p,e=1 u=89 vol=4.30E+02 $Ankle+Feet 90 103 -1.319 -999 imp:p,e=1 u=90 vol=1.50E+02 $Humerus-proximal 91 103 -1.319 -999 imp:p,e=1 u=91 vol=7.40E+01 $Humerus-upper shaft 92 104 -1.378 -999 imp:p,e=1 u=92 vol=6.70E+01 $Humerus-lower shaft 93 104 -1.378 -999 imp:p,e=1 u=93 vol=9.30E+01 $Humerus-distal 94 106 -1.374 -999 imp:p,e=1 u=94 vol=1.60E+01 $Radii-proximal PAGE 135 135 95 106 -1.374 -999 imp:p,e=1 u=95 vol=1.60E+01 $Radii-upper shaft 96 106 -1.374 -999 imp:p,e=1 u=96 vol=1.70E+01 $Radii-lower shaft 97 106 -1.374 -999 imp:p,e=1 u=97 vol=2.70E+01 $Radii-distal 98 106 -1.374 -999 imp:p,e=1 u=98 vol=4.90E+01 $Ulnae-proximal 99 106 -1.374 -999 imp:p,e=1 u=99 vol=2.70E+01 $Ulnae-upper shaft 100 106 -1.374 -999 imp:p,e=1 u=100 vol=1.80E+01 $Ulnae-lower shaft 101 106 -1.374 -999 imp:p,e=1 u=101 vol=7.00E+00 $Ulnae-distal 102 105 -1.415 -999 imp:p,e=1 u=102 vol=8.60E+01 $Hand 103 105 -3.00 -999 imp:p,e=1 u=103 vol=1.30E+01 $Teeth 301 1 -0.001205 -999 imp:p,e=1 u=301 vol= 78406.867 $ Air 302 28 -0.2995 -999 imp:p,e=1 u=302 vol= 10210.44 $ Lung 303 2 -1.00 -999 imp:p,e=1 u=303 vol= 27482.027 $ Tissue 304 101 -1.32 -999 imp:p,e=1 u=304 vol= 2304.37 $ Bone C ******************* Lattice definition ************************************ 1001 0 -100 fill=10000 imp:p,e=1 $ surrounding box 555 0 -200 lat=1 u=10000 imp:p,e=1 fill = 0:59 0:26 0:166 C **************** Image data start from here ******************************* Phantom Input from Ghoast-3-D C ********************************************************************************* c c Surface Cards c C ************************************************************************** 100 rpp 0 60 0 27 0 167 $(-29:30 -13:13 -83:83) 200 rpp 0 1.0 0 1.0 0 1.0 999 rpp -25 85 -25 52 -25 192 c ************************************************************************************** Material Cards C*************************************************************************************** mode p e sdef x=d1 y=d2 z=d3 erg=d4 AXS= 0 1 0 dir= 0.99763 vec= 0 1 0 si1 H 30 40 sp1 D 0 1 si2 H 0.0 1.0 sp2 D 0 1 si3 H 119 136 sp3 D 0 1 si4 H 0.01 8.0 sp4 D 0.0 1.0 c -------------------------------------------------------------------c tallies cards c --------------------------------------------------------------------c ---------------------------------------------------phys:p 4j 1 ctme 500 totnu print PAGE 136 136 LIST OF REFERENCES 1R. Mohan, in XIIth International Confer ence on the Use of Computers in Radiation Therapy edited by (Medical Physics Publishing, Madiso n, WI, 1995, Salt Lack City, 1997), pp. 16-18.35 2T. R. Mackie, P. Reckwerdt, and N. Papani kolaou, "3-D Radiation Treatment Planning and Conformal Therapy," Medical Physics Publishing Corporation Madison WI, 1995), pp. 34 3B. A. Fraass, J. Smathers, and J. Deye, "Su mmary and recommendations of a National Cancer Institute workshop on issues limiting the clinical use of Monte Carlo dose calculation algorithms for megavoltage external beam radiation therapy," Med Phys 30, 3206-3216 (2003).29 4F. M. Khan, The Physics of Radiation Therapy, 3 (Lippincott Williams & Wilkins, 2003).36 5P. R. Almond, A. E. Wright, and M. L. Boone "High-energy electron dose perturbations in regions of tissue heterogeneity. II. Physical models of tissue heterogeneities," Radiology 88, 1146-1153 (1967).37 6A. Ahnesjo, M. Saxner, and A. Trepp, "A pencil beam model for photon dose calculation," Med Phys 19, 263-273 (1992).73 7A. L. EDWARDS, J. A. RATHKOPF, and R. K. SMIDT, "Extending the Alias Monte Carlo sampling method to general distributions,," in the International Conference on Mathematics and Computations and Reactor Physics edited by (American Nuclear Society, Inc., La Grange Park, IL, 1991), pp. 38 8A. Ahnesjo, "Collapsed cone convolution of radiant energy for phot on dose calculation in heterogeneous media," Med Phys 16, 577-592 (1989).52 9T. Knoos, A. Ahnesjo, P. Nilsson, and L. Weber, "Limitations of a pencil beam approach to photon dose calculations in l ung tissue," Phys Med Biol 40, 1411-1420 (1995).50 10B. Dionne, "Automated variance reduction tech nique for 3-D Monte Carlo coupled electronphoton-positron simulation using determ inistic importance functions," Nuclear Engineering Sciences (University of Florida, Gainesville, Fl, 2007), pp. 75 11E. E. Lewis, and W. F. Miller, Computational Methods of Neutron Transport ( American Nuclear Society, La Grange Park, Ill., USA, 1993).18 12L. J. Shukovsky, "Dose, time, volume relations hips in squamous cell carcinoma of the supraglottic larynx," Am J Roen tgenol Radium Ther Nucl Med 108, 27-29 (1970).31 13J. J. Fischer, and J. E. Moulde r, "The steepness of the dose-res ponse curve in radiation therapy. Theoretical considerations and e xperimental results," Radiology 117, 179-184 (1975).30 PAGE 137 137 14B. Wang, G. Xu, J. T. Goorley, and A. Bozkurt, "ISSUES RELA TED TO THE USE OF MCNP CODE FOR AN EXTREMELY LARGE VOXEL MODEL VIP-MAN," in The Monte Carlo Method: Versatility Unbounded In A Dynamic Computing World, edited by (American Nuclear Society, Chattanooga, Tennessee, 2005), pp. 32 15G. I. Bell, and S. Glasstone, Nuclear Reactor Theory, (Robert E. Krieger Publishing CO. Inc, Malabar, FL, USA, 1985).39 16B. G. Carlson, and K. D. Lathrop, Discrete Ordinates Angular Quadrature of the Neutron Transport Equation (Los Alamos Scientific Laboratory Report, LA-3186, 1965).40 17D. Shedlock, and A. Haghighat, Neutron Analysis of Spent Fuel Storage Installation Using Parallel Computing and Advance Discrete Ordinates and Monte Ca rlo Techniques," Radiat Prot Dosimetry 116 662-666 (2005).45 18K. D. Lathrop, "Remedies for Ray Effects," Nuclear Science and Engineering 45, 255-268 (1971).41 19G. Sjoden, and A. Haghighat, "PENTRAN A 3D Cartesian Parallel SN Code with Angular, Energy, and Spatial Decomposition," in Intl. Conf. on Math. Meth & Supercomp. for Nuclear App., edited by Saratoga Springs, New York, 1997), pp. 553.14 20A. Al-Basheer, M. Ghita, G. Sjoden, W. Bolch, and C. Lee, "Whole Body Dosimetry Simulations using the PENTRAN-MP Sn Code System," in American Nuclear Socity 2007 Annual Meeting edited by Bosten, MA, 2007), pp. 66 21G. M. Daskalov, R. S. Baker, R. C. Little, D. W. O. Rogers, and J. F. Williamson, "TwoDimensional Discrete Ordinates Photon Trans port Calculations for Brachytherapy Dosimetry Applications," Nucl. Sci. Eng. 134, 121-134 (2000).19 22G. M. Daskalov, R. S. Baker, D. W. Rogers and J. F. Williamson, "Multigroup discrete ordinates modeling of 125I 6702 seed dose distributions using a broa d energy-group cross section representation," Med Phys 29, 113-124 (2002).20 23M. L. Williams, D. Ilas, E. Sajo, D. B. Jones, and K. E. Watkins, "Deterministic photon transport calculations in gene ral geometry for external beam radiation therapy," Med Phys 30, 3183-3195 (2003).5 24A. Haghighat, International Training Course/Workshop on Methodologies for Particle Transport Simulation of Nuclear Systems 76 25(ICRU) 1980, Radiation quantities and units (International Commission on Radiation Units and Measurements, Bethesda, MD, 1980).42 26F. H. Attix, "The partition of kerma to account for bremsstrahlung," Health Phys 36, 347-354 (1979).43 PAGE 138 138 27F. H. Attix, "Energy imparted, energy transferred and net energy transferred," Phys Med Biol 28, 1385-1390 (1983).44 28N. I. o. S. a. T. NIST, ESTAR(Stopping Power and Range Tables for Electrons ) http://physics.nist.gov/PhysRe fData/Star/Text/E STAR.html.71 29T. R. Mackie, A. F. Bielajew, D. W. Rogers, and J. J. Battista, "Generation of photon energy deposition kernels using the EGS Monte Carlo code," Phys Med Biol 33, 1-20 (1988).54 30A. Ahnesjo, P. Andreo, and A. Brahme, "Calculation and application of point spread functions for treatment planning with high energy photon beams," Acta Oncol 26, 49-56 (1987).53 31L. A. N. L. X-5 Monte Carlo Team, MCNP5A General Monte Carlo N-Particle transport code, Version 5 10 32J. E. White, J. D. Intgersoll, C. Slater, and R. Roussin, BUGLE-96 : A Revised Multi-group Cross Section Library for LWR Applications Based on ENDF/B-VI 46 33L. J. Lorence, J. E. Morel, and G. D. Valdez, Physics Guide to CEPXS: A multigroup coupled electron-photon cross section generating code 11 34G. SJODEN, and A. HAGHIGHAT, PENTRAN: P arallel E nvironment N eutral-Particle TRAN sport Sn in 3-D Cartesian Geometry 33 35C. Yi, and A. Haghighat, PENMSH Express 62 36G. Longoni, "Advanced quadrature sets, acceleration and preconditioning techniques for the discrete ordinates method in parallel computing environments ," (University of Florida, Gainesville, FL, 2004), pp. 47 37B. Petrovic, and A. Haghighat, "New directi onal theta-weighted (DTW) differencing scheme and reduction of estimated pressure vessel fluence uncertainty," in the 9th International Symposium on Reactor Dosimetry edited by Prague, Czech Republic, 1996), pp. 746-753.48 38G. Sjoden, and A. Haghighat, "The exponentia l directional weighted (EDW) differencing scheme in 3-D cartesian geometry," in The Joint International Conference on Mathematics and Supercomputing for Nuclear Applications edited by Saratoga Springs, New York, 1997 ), pp. 21267-21276.49 39G. Sjoden, "An Efficient Exponential Directio nal Iterative Differencing Scheme for ThreeDimensional Sn Computations in XYZ Geometry," Nuc. Science and Eng 155, 179-189 (2007).12 40G. Sjoden, and A. Haghighat, "Taylor Project ion Mesh Coupling Between 3-D Discontinuous Grids for Sn," Transactions of the American Nuclear Society 74, 178-179 (1996).77 PAGE 139 139 41T. Tabata, R. Ito, and S. Okabe, "Generalized Semiempirical Equations For The Extrapolated Range of Electrons," Nucl. Instrum. Methods ;103: No. 1, 85-91(1972). (1972).65 42A. Al-Basheer, M. Ghita, and G. Sjoden, "Inves tigation of deterministic simulation parameters for a 3-D radiotherapy Co-60 device model," in ANS Winter Meeting and Nuclear Technology Expo, edited by Albuquerque, NM, 2006), pp. 511-513.13 43C. Lee, C. Lee, J. L. Williams, and W. E. Bolch, "Whole-body voxel phantoms of paediatric patients--UF Series B," Phys Med Biol 51, 4649-4661 (2006).7 44ICRU 1992 Photon, electron, proton and ne utron interaction data for body tissue (International Commission on Radiation Units and Measurements, Bethesda, MD, 1992).15 45M. GHITA, A. Al-BASHEER, and G. SJODEN "Validation of DXS: A Diagnostic X-Ray Source Generator For X-Ray Transport Simulations," in The American Nuclear Society's 14th Biennial Topical Meeting of the Radi ation Protection and Shielding Division edited by Carlsbad, New Mexico, 2006), pp. 9 46D. M. Tucker, G. T. Barnes, and D. P. Chakraborty, "Semiempirical model for generating tungsten target x-ray spectra," Med Phys 18, 211-218 (1991).6 47M. GHITA, DXS: A Diagnostic X-ray Souce Generator 74 48C. M. Ma, E. Mok, A. Kapur, T. Pawlicki, D. Fi ndley, S. Brain, K. Forster, and A. L. Boyer, "Clinical implementation of a Monte Carl o treatment planning system," Med Phys 26, 21332143 (1999).3 49W. Gropp, E. Lusk, and A. Skjellum, Using MPI: Portable Parallel Programming with the Message-Passing Interface, (The MIT Press, Cambridge, Massachussetts, 1994).72 50A. K. Jones, F. D. Pazik, D. E. Hintenlang, and W. E. Bolch, "MOSFET dosimeter depth-dose measurements in heterogeneous tissue-equivalent phantoms at diagnostic x-ray energies," Med Phys 32, 3209-3213 (2005).69 51C. Lee, D. Lodwick, D. Hasenauer, J. L. Williams, C. Lee, and W. E. Bolch, "Hybrid computational phantoms of the male and female newborn patient: NURBS-based whole-body models," Phys Med Biol 52, 3309-3333 (2007).70 PAGE 140 140 BIOGRAPHICAL SKETCH Ahm ad Al-Basheer was born in 1978 in Ramtha City, Jordan, in the heart of the Middle East. Generations of his family have made Ra mtha their hometown for centuries. This city, which borders Syria and Jordan, is the gatewa y to the desert region of the Middle East. Ahmad is the oldest of eight children, a nd his parents, Khaled and Fatima, are both involved in educational careers ; both are school principals. Ahmad was involved not only in traditional educational venues, but also in chess, soccer, reading and poetry. Ahmad competed in several reading and poetry contests throughout his schooling; he won the National Award for high school students for Best Poet. Ahmad graduated from Jordan University fo r Science and Technology with a Bachelor of Science degree in applied physics in 2000. Upon gr aduation, he spent one year as a physics teacher at the Yarmouk University Model Sc hool. In January, 2002, Ahmad was granted admission to the graduate program in the Depart ment of Nuclear and Ra diological Engineering at the University of Florida. He performed his research work under the supervision of Dr. Samim Anghaie, where he investigated a new approach for reducing scattered photons and electron contamination in Cobalt-60 therapy beam. In fall 2004 he joined the Un iversity of Florida Transport Theory Group (UFTTG) to pursue his Ph.D. with Dr. Glenn Sjoden. Duri ng his Ph.D., He won several honors, including best papers presented to numb er of conferences, a number of invited speeches at national conferences, and several recogniti ons form number of organizati on for his achievements. Ahmad married Ghadeer Al-Basheer, on June 16th 2007 and has one daughter, Ju de Al-Basheer who was born in June 14th 2008. |