MATLAB Monte Carlo
User’s Manual
Medical Physics 501
December 15, 2008
Purpose/Introduction
The goal of this project was to create a functional Monte Carlo program using the Matrix Laboratory, MATLAB. Using our understanding of the underlying physics from the Medical Physics 501 course, an Introduction to Radiological Physics and Dosimetry, we set forth to complete this task within the semester as a whole class, each of us taking a portion, leaving our individual fingerprints. The project was broken into three subgroups, two to handle the photon interaction physics, one for beta particle interactions. What resulted appears to be a functional program which returns the dose distribution of a volume as defined by an input DICOM file. We hope that future students can use our program as a base to create an even more powerful Monte Carlo tool, taking into account more advanced physics and concepts that one semester of work would not allow. That being said, we present to you the guide for MATLAB Monte Carlo, version 1.0.
Table of Contents
Background Physics...... 2 User Interface Description...... 6 Program Flowcharts...... 8 Function Descriptions...... 12 List of References...... 20
Background Physics
Photon Interactions with Matter
There are three dominant types of photon interactions with matter in the energy range of interest: the Compton effect, the photoelectric effect, and pair production. In all of these interactions, a portion of the photon’s energy is transferred to electrons, which can then deposit their energy in matter. The type of interaction depends on the initial energy of the photon and the atomic number of the interacting medium.
To determine the type of interaction for our simulation we used the contribution of each interaction type to the total interaction cross section. These cross-sections are defined per unit mass of material giving an expression for the total mass attenuation coefficient for a specific material at specific photon energy. It can be written in units of cm2/g or m2/kg as
in which τ/ρ is the contribution of the photoelectric effect, σ/ρ that of the Compton effect, and к/ρ that of pair production. The mass attenuation coefficient tables come from the NIST webpage (http://physics.nist.gov/PhysRefData/Xcom/html/xcom1.html).
Several materials can be used for the homogenous material in the program, including air, muscle, fat, bone, tungsten, lead, water, and titanium. For each material, tables are used to find the total mass attenuation coefficient (), photoelectric mass attenuation coefficient (), Compton mass attenuation coefficient (), and pair production mass attenuation coefficient (). All of the mass attenuation coefficients are in the unit of cm2/g.
Compton Interaction
A Compton interaction describes the process of a photon scattering off of an atomic electron. Although the electron is theoretically bound, it is assumed to be unbound, and the collision elastic. This allows the use of basic mechanics to describe the final particles.
The Klein-Nishina cross section is used to define the polar scattering angle of the photon. This is done using a rejection technique, and can be found in the function GetCosTheta. The rejection technique favors forward scattered photons at high energies and isotropic scattered photons at low energies. The rejection technique uses a random polar angle to compute the angular differential scattering cross section for Compton scattering, and if the result obeys the known statistics, the angle is kept. If it does not, the process is repeated with a new random angle.
Once a polar scattering angle is found, the energies and scattering angles of both the photon and electron can be determined in main Compton routine, ComptonInteraction. These are calculated using conservation laws. The azimuthal direction is isotropic, and randomly generated. With the new photon and electron angles, the new direction vectors for both particles are calculated. If the energy of the scattered electron is less than a user defined cut-off energy, Δ, the electron energy is deposited. If the electron is not deposited, it is pushed to the electron stack to be transported. The routine is complete when final electron and photon energies and directions have been found.
Photoelectric Interaction
A photoelectric interaction occurs when an atomic electron absorbs a photon and is emitted from the atom. The energy of the electron is equal to hν – Eb, where Eb is the electron shell binding energy. Secondary radiation occurs as a result of the vacancy in the electron cloud in the form of fluorescence photons and Auger electrons.
In order to determine which shell the photon interacts with, tables containing participation function values were used (medium.tables.photoelectric_). If a random number is less than Pk, the interaction will occur in the k-shell. Otherwise, it will occur in the l-shell.
In order to determine the energy of the emitted electron, k-edge and l-edge tables are used (medium._edge). If the photon interacted in the k-shell, an electron is created with the energy hν-(k-edge). Otherwise, the photon interacted in the l-shell, and an electron with energy hν-(l-edge) is created.
Photoelectric emission is isotropic, and both polar and azimuthal angles are generated randomly. From these angles, direction vectors are given to the newly created electron. The electron is now completely defined, with an energy, direction, and position.
At this point, there is still a small amount of unaccounted energy. The effects of the secondary radiation as a result of the electron shell vacancy have been neglected. This was decided because the shell edge energies are small, and any Auger electrons emitted would likely not leave the voxel. The fluorescence photons that are not accounted for should even out over the entire volume. This remaining energy is deposited in the voxel in which the interaction occurred.
Pair Production
In a pair production interaction, an electron-positron pair is created as a result of the highly energetic photon interacting with the nucleus. The photon must exceed a threshold energy of 2moc2. Triplet production can also occur when a photon interacts with the field of an electron. This process is much less common and occurs only at high energies.
Although the electron and positron are not necessarily given the same energy, the PairProductionInteraction routine assumes both particles leave the interaction site with equal energies. These energies can be computed using (hν-2moc2)/2. The momentum of the nucleus as a result of the interaction is neglected.
The polar angles (alpha-positron, beta-electron) associated with the created particles are computed using non-relativistic expressions of kinetic energy and momentum. While the azimuthal angles of the particles are isotropic, they must be 180o apart in order to conserve momentum. Once one random angle was computed, 180o was added to get the other angle. With the angles now determined, direction vectors were calculated, and the electron and positron were both completely defined. The electron and positron were pushed to their respective stacks.
Beta Particle Interactions with Matter
Mass Collision Stopping Power for Electrons and Positrons
The mass collision stopping power is a differential that gives the energy lost per density-distance for a particles traveling through a medium. The mass collision stopping power for an electron or a positron is given by:
where
The dependence on medium manifests itself in the Z/A and I (mean excitation potential) dependence. δ is the polarization effect correction. The polarization effect correction is zero for gases and increases with increasing particle speed. In liquids and solids, the density is roughly 104 higher than in a gas at atmospheric pressure. As an electron travels through a liquid or solid media, the polarization of the medium near the track of the electron weakens the effective electric field strength felt by the particle. Therefore, the mass collision stopping power is reduced, and is accounted for with the δ in the equation above. The dependence on Z/A has the greatest effect on the magnitude of the mass collision stopping power. For example, 1H has a Z/A value of 1 which makes it well suited for attenuating charged particles.
Restricted Mass Collision Stopping Power
The restricted mass collision stopping power takes into account the delta rays that are produced that leave a certain voxel of the medium. Whether or not the delta rays leave the voxel depends on the energy of the delta-ray. Taking this into account is important because delta-rays with high enough energy will leave the voxel and will not deposit their energy in that voxel. This is important for dosimetry measurement considerations. The calculation for the restricted mass collisional stopping power is calculated in a similar to the equation above; however, is replaced by , where .
Mass Radiative Stopping Power
The rate of bremsstrahalung production by electrons and positrons is given by the mass radiative stopping power, in units of MeV cm2/g. , where Z2/A refers to the atomic number and mass number of the medium and (moc2)2 refers to the square of the rest energy of the charged particle. Since the mass of the proton is ~2000 times larger than the mass of an electron, no appreciable bremsstrahalung will be produced. The Z2 dependence indicates that bremsstrahalung is much more likely to occur in high Z materials, as compared to low Z materials.
Continuous Slowing Down Approximation (CSDA) Range
The CSDA can be used to approximate the range of a particle. The RCSDA is defined as:
which is an integral over the inverse of the total mass stopping power from zero to the energy of the particles. Therefore, a medium with a high mass stopping power will result in a smaller range than a medium with a low mass stopping power.
Delta Ray Production
Delta rays are produced when an incoming electron undergoes a hard collision with an unbound or valence shell electron. Since electrons are indistinguishable from one another, the outgoing electron with the highest energy is defined to be the primary electron. The lower energy electron is therefore defined as a delta-ray or knock-on electron. Delta rays can be produced by both electrons and positrons, though the differential cross sections are slightly different. By definition, the maximum kinetic energy a delta ray can have is T/2. For a positron, the maximum kinetic energy imparted to a delta-ray can be T, because the delta ray is distinguishable from the positron. The differential cross section for electrons (Møller cross section) and the differential cross section for positrons (Bhabha cross section) can be seen below:
Where σ- refers to the inelastic cross section for an electron and σ+ refers to the inelastic cross section for the positron. The differential cross sections will be discussed in more detail below.
Graphical User Interface
The Graphical User Interface (GUI) provides a user friendly portal through which to run the Monte Carlo simulation. The GUI and its features are shown below.
Monte Carlo Graphical User Interface
A) DICOM path window
This window displays the path of the first slice of the inputted DICOM file after selection. For DICOM selection, see PET/SPECT toggle buttons.
B) PET toggle button
Click this button if you wish to input a PET DICOM file. After selection of this button, a window will pop up which will allow you to select one or more DICOM files. The path of the first DICOM file will be shown in the DICOM path window.
C) SPECT toggle button
Select this button if you wish to input a SPECT DICOM file. After selection of this button, a window will pop up which will allow you to select one or more DICOM files. The path of the first DICOM file will be shown in the DICOM path.
D) PET tracer menu
If the PET toggle button is selected, the PET tracer menu will allow you to select from a pre-determined set of PET tracers; 15O, 14N, 11C or 18F. The selected PET tracer will affect the outcome of the simulation; each PET tracer option emits positrons with a different initial energy. This initial energy will determine how far the positron travels before annihilation.
E) SPECT tracer menu
If the SPECT toggle button is selected, the SPECT tracer menu will allow you to select from a pre-determined set of SPECT tracers. These tracers include: 99mTc, 123I and 111In. The selected SPECT tracer will affect the outcome the simulation by determining the energy of the initial photon emitted from the tracer.
F) Material Type
This menu allows you to select the material in which you would like to run the simulation. The material type selection will affect the mean free path of the photons created in the simulation as well as the mass stopping power of the charged particles. Charged particles are also affected by the material type as dictated by the radiation yield and angular scattering dependence.
G) T. Cutoff
The T. Cutoff menu allows you to select from a predetermined set of options, the energy below which, individual delta rays are not tracked. See Beta section of user manual for more information on how T. Cutoff is used.
H) ELF
The Energy Loss Fraction (ELF) menu allows you to choose from a predetermined list of the fraction of energy lost by the charged particle before checking again for interactions.
I) Plot Particles
The Plot Particles checkbox allows for a certain number of particles in the simulation to be plotted in a 3D plot after the simulation has run. Checking of this box also requires the number of particles you wish to plot.
J) Number of Particles to Simulate
The number of particles you wish to simulate should be entered here.
K) Run
The run button is enabled all necessary fields have been entered. Run creates a structure in MatLab of all the user defined variables, passes it to the main function of the Monte Carlo code, and begins the simulation.
Default Settings
Material Type: Air
Number of Particles to Simulate: 1000
Number of Particles to Plot: 10
Plot Particles: Off
Cutoff Energy for delta rays: .01 keV
Energy Loss Fraction: 5%
Main Program Flowchart
Photon Interactions
Electrons