Monte Carlo calculations of Linear Energy Transfer based on the PENELOPE code
Master’s thesis
Leah Dickhoff
Supervised by Anders Lundin Torbj¨orn B¨ack
June 2020
I would like to express my sincere thanks to my supervisor Anders Lundin for his constant guidance, help and feedback throughout this project. Mostly his dedication, but also the openness and expertise of the whole Elekta team has been crucial to this project and to my professional development. I further owe my supervisor Tobj¨orn B¨ack at KTH gratitude, both for bringing me into contact with Elekta, as well as for his support and valuable advice regarding this thesis.
This thesis work has been carried out at Elekta in order to implement a Linear Energy Transfer (LET) scoring module for protons and ions into the in-house Monte Carlo code Pegasos. The module is divided into two parts for separate calculations of the dose-averaged LET and the track-averaged LET.
Evaluations of the results for a specific proton beam were split up into the LET of primary particles, of each kind of generated secondaries, as well as of all generations. Secondary particles manifest a high LET, and particularly secondary protons contribute to a post-Bragg peak LET tail.
The study of the individual interaction types allowed to conclude that the main contribution emanates from electronic interactions, while the inclusion of nuclear reactions necessitates further restrictions on e.g. scored lengths.
The carbon-12 ion LET was found to be on average 20 times larger than the proton LET of similar range beams, and the incorporation of secondaries leads to a substantial well-averaged LET tail.
Dissimilarities between the dose- and track-averaged LET were uncovered to be caused by the larger dependence of the dose-averaged LET on the variance of the deposited energy. The track-averaged LET displayed a more stable averaging, while the dose-averaged LET showed fluctuations and elevated values in case of polyenergetic particles.
As conveyed with the Relative Biological Effectiveness (RBE) calculations on a CT-based geometry, the implemented module opens possibilities for an incorporation of the LET in planning systems. Further potential works in- clude comparisons with future lineal energy modules, as well as studies of the LET’s biological implications through its relationship with the RBE.
1 Introduction 5
1.1 Proton therapy and carbon ion therapy . . . 5
1.2 Elekta . . . 7
1.3 Aim of thesis . . . 8
2 Proton and Ion Transport Physics 10 2.1 Charged particle interactions . . . 10
2.1.1 Inelastic Coulomb scattering . . . 10
2.1.2 Elastic Coulomb scattering . . . 12
2.1.3 Non-elastic nuclear interactions . . . 13
2.1.4 Bremsstrahlung . . . 13
2.2 Basic concepts . . . 14
2.2.1 Stopping power . . . 14
2.2.2 Range . . . 15
2.2.3 Cross section and mean free path . . . 16
3 Physical Quantities in Radiotherapy 18 3.1 Fluence . . . 18
3.2 Dose . . . 19
3.3 Linear Energy Transfer . . . 19
3.3.1 Track- and dose-averaged LET . . . 20
3.4 Relative Biological Effectiveness . . . 21
3.5 Lineal energy . . . 24
4 Monte Carlo Methods for Charged Particle Transport 26 4.1 Monte Carlo techniques . . . 26
4.2 Code systems . . . 27
4.2.1 PENELOPE . . . 27
3
CONTENTS 4
4.2.2 PENH . . . 28
4.2.3 Geant4 . . . 28
4.3 Pegasos . . . 28
4.3.1 Input . . . 29
4.3.2 Simulation . . . 29
4.3.3 Mixed Monte Carlo . . . 31
4.4 Linear energy transfer calculations . . . 32
4.4.1 Track-averaged LET . . . 32
4.4.2 Dose-averaged LET . . . 33
5 Methods 35 5.1 Proton LET . . . 35
5.1.1 Code implementation . . . 35
5.1.2 Particle generations and particle types . . . 35
5.1.3 Types of interaction . . . 37
5.1.4 Choice of scoring method . . . 37
5.1.5 Pegasos input . . . 39
5.2 Carbon ion LET . . . 42
5.3 Relative Biological Effectiveness . . . 42
5.4 Statistics . . . 44
5.4.1 Standard deviations . . . 44
5.4.2 Chi-square statistic . . . 47
6 Results and Discussion 49 6.1 Proton LET . . . 49
6.1.1 Primary particle LET . . . 49
6.1.2 Interaction types . . . 58
6.1.3 Secondary particles . . . 62
6.2 Carbon Ion LET . . . 68
6.2.1 Primary particle LET . . . 68
6.2.2 Post-Bragg peak tail . . . 70
6.3 Proton RBE . . . 71
7 Conclusions 75 7.1 Outlook . . . 77
Introduction
1.1 Proton therapy and carbon ion therapy
Radiation therapy is an indispensable tool in treating cancer and, together with surgery and chemotherapy, is part of most cancer patients’ treatment options. The principle behind targeting cancerous cells with radiation is sim- ple: the particles deposit energy as they pass through the body and thereby inflict damage on the DNA strands which possibly kills the cells. There are two different types of treatment classifications in which this principle is applied. Firstly, in brachytherapy, the radiation source is surgically placed near the targeted tumour, and radiates from within. Secondly, in external beam therapy, an external source is directed towards the target volume from outside the body, so most often penetrates through tissue in order to reach the cancerous cells.
Even though the first experiments of radiation therapy range back to the late 1800’s, the idea of using accelerated ions for cancer therapy was originally only advanced in 1946 by Wilson [1], who then was a student of Ernest Orlando Lawrence at the University of California in Berkeley. He realised the benefits of the dose-depth distribution resulting from ion beams (see Figure 1.1), since the deposition of energy per unit track length increased with depth. Therefore, as compared to e.g. photons, a much lesser dose is given to the normal, healthy tissue for the same dose given to the target area. Hence, it is not as imperative to use multiple beams from different angles (as is the case in photon treatments), and only a few beams can be used to achieve the same target dose. This is due to the fact that ions deposit less energy when they have a high velocity - as they slow down, they become
5
CHAPTER 1. INTRODUCTION 6
Figure 1.1: Relative dose depositions by photons, protons and carbon ions of different energies in water. Curves obtained with Elekta’s in-house Monte Carlo system Pegasos.
increasingly slow so deposit most of their energy at the end of their track.
The resulting peak in dose is called the Bragg peak. By stacking multiple beams of different energy and intensity, the width of the Bragg peak can be increased to adapt to each individual clinical case, leading to a Spread-Out Bragg Peak (SOBP). Lateral differences in the target tumor can also be taken into account through the use of modulators.
Currently, the most common type of charged particle therapy is proton therapy, using the lightest and simplest charged particle: the hydrogen ion.
One of the reasons for its popularity is the approximately non-existent dose beyond the Bragg peak. This type of radiation treatment is closely followed by carbon ion therapy, using carbon-12 ions. Since carbons are heavier ions, they are particularly powerful at inflicting damage on the cells - though both the healthy and the cancerous ones. They are also known to generate a non-zero post-Bragg peak dose, caused by numerous secondary particles generated during nuclear reactions between the carbon ions and the nuclei in the patient’s body.
The first proton therapy treatments were conducted in 1954 and 1957 at the Berkeley Radiation Laboratory, and in Uppsala in Sweden, respectively, where they used nuclear physics research equipment including particle ac- celerators. Over the years, more and more clinical proton therapy centres were built (the first one in 1989 at the Clatterbridge Cancer Centre in the UK), and the number of facilities in operation has since increased to 113 (as of April 2020) [2]. The development of proton therapy has been con- siderably slower than e.g. photon therapy since the need for synchrotrons or cyclotrons as particle accelerators makes the necessary machinery more expensive and the setup more technically difficult. For carbon ion therapy, this problem becomes even more substantial, since larger energies are needed for the carbon ions to reach the same range as protons. Hence, both pro- ton therapy and carbon ion therapy are currently lacking overall evidence of cost-competitiveness [3].
1.2 Elekta
Elekta is an international Stockholm-based company and provider of medical radiation treatment systems, within branches of radiosurgery, neurosurgery, radiotherapy and brachytherapy. Elekta was founded in 1972 by Lars Leksell, a Professor of Neurosurgery at Karolinska Institutet, and now counts more than 3000 employees worldwide. The company’s most famous product is the Leksell Gamma Knife, used for stereotactic radiosurgery of brain tumours. It attains its high precision and power through numerous gamma emitting 60Co samples which are focused on the target. Other Elekta products include linear accelerators (linacs), the Versa HD for whole body radiotherapy, software for treatment delivery and planning called Monaco , and brachytherapy products. Finally, Elekta Unity is a real timeR magnetic resonance radiation therapy MR-linac allowing for treatments of high mobility organs such as the lungs.
Monte Carlo codes for detailed reference simulations for Elekta’s products are run using the in-house code system called Pegasos. It has been adapted by Anders Lundin from the PENELOPE code, which was originally developed by F. Salvat [4]. At that point, Pegasos only handled electron, photon and positron transport which were needed for reference simulations for the gamma knife and the linacs. Later, the PENH library (also developed by Salvat [5]) was added to Pegasos by thesis interns supervised by Anders Lundin. For
CHAPTER 1. INTRODUCTION 8 instance, Ludvig Hult investigated the physics (especially the eikonical model for electromagnetic elastic scattering) used in the added PENH [6]. Martin Bj¨ornmalm implemented nuclear reactions and the secondaries thereof in PENH, and benchmarked it against experimental data from Prague Proton Therapy centre [7]. David Blomqvist added neutron transport [8], and Oscar Pastor Serrano included light ion transport into Pegasos [9].
1.3 Aim of thesis
It has been established that proton therapy is the kind of charged particle therapy most used, followed by carbon ion therapy. Such charged particle therapy works by the particles depositing energy in the cells they pass through. Thereby, DNA strands are broken, which can lead to cell mutations and cell deaths. As such, the concept of dose as deposited energy per unit mass serves as a good quantification of the inflicted damage, and is used in most of today’s radiation therapy planning systems. A scoring module of the dose has therefore also already been introduced in the earliest versions of Pegasos. However, the real biological effect of any kind of radiation cannot accurately be measured by just the dose, since it depends on numerous other factors as well, such as the type of radiation, the quality of radiation, the type of tissue, the cell type, or the specific biological repair mechanisms. In 1931 [10], this lead to the introduction of a quantity which is meant to take into account all these properties, and thereby measures a radiation’s biological effect: the Relative Biological Effectiveness (RBE) (see section 3.4).
Numerous RBE models have been developed over the years, and, because of the wide range of influencing factors which should be included, they all have one thing in common: their complexity. However, most of these models show a close dependency of the RBE on more easily computable quantities.
One example is the Linear Energy Transfer (LET), roughly defined as the energy that a charged particle deposits per unit travel distance. Thus, calculations of the LET values for a specific beam serve as a good first step towards calculations of the RBE, and thereby quantifications of the actual biological effect each beam causes. Therefore, the aim of this thesis is the following:
To introduce a scoring module for LET calculations into Pegasos, and study its behaviour by mainly focusing on proton therapy, and briefly on carbon ion therapy.
The structure of the resulting work is as follows: First, the proton transport physics are described, by separately explaining all kinds of interaction mechanisms, as well as later all related concepts. Then, an overview of some physical quantities commonly used in relation with radiotherapy is given. Next, Monte Carlo methods are introduced by describing different code systems, while Elekta’s Pegasos is described in detail. In the same chapter, a summary of LET calculations within common Monte Carlo codes is given. Subsequently, the methods used for the LET module implementation are put forward, and in the following chapters, the corresponding results are presented and discussed.
Chapter 2
Proton and Ion Transport Physics
2.1 Charged particle interactions
When ions pass through matter, they interact with the atoms they encounter through multiple different mechanisms. These interactions are both due to the ion’s positive charge and the associated electronic interactions with the medium’s charged electrons and ions, as well as the actual collisions with the medium’s nuclei. As such, ions undergo four distinct mechanisms:
inelastic Coulomb scattering, elastic Coulomb scattering, non-elastic nuclear interactions, and bremsstrahlung. Their respective target in the medium, ejectiles of the interactions, and influence on the primary particle and on the dose are given in Table 2.1, while Figure 2.1 gives a schematic overview of the different non-negligible mechanisms.
2.1.1 Inelastic Coulomb scattering
During inelastic scattering, the internal state of projectile and target is changed. Inelastic Coulomb scattering takes place when the projectile interacts with atomic electrons and is the main mechanism for energy loss of charged particles passing through matter. Since the proton rest mass is around 1836 times larger than the electron rest mass, and the carbon-12 ion even 21874 times larger, little energy is lost per single proton/ion-electron interaction, and the projectile charged particle moves in approximately straight line. In order to model multiple interactions of a particle along
10
p
p e
(a)
p
p
(b)
p
p’
n
recoil nucleus γ
(c)
Figure 2.1: Charged particle interaction mechanisms: (a) inelastic Coulomb interaction, (b) elastic Coulomb interaction, (c) inelastic nuclear interaction (p: charged particle, e: electron, n: neutron, γ: gamma rays).
CHAPTER 2. PROTON AND ION TRANSPORT PHYSICS 12
Interaction type
Target Ejectiles Influence
on primary particle
Influences on dose
Inelastic Coulomb scattering
Atomic electrons
Primary particle, ionized electrons
Quasi- continuous energy loss
Range in target
Elastic Coulomb scattering
Atomic nucleus Primary particle, recoil nucleus
Change in direction
Lateral penumbra
Non-elastic nuclear interactions
Atomic nucleus Secondary particle, heavier ions, neutrons, and gamma rays
Removal of primary ion
Fluence of primary particle, stray neutrons, prompt gamma rays
Bremsstrahlung Atomic nucleus Primary particle, Bremsstrahlung photon
Energy loss, change in direction
Negligible
Table 2.1: Charged particle interaction types, targets, ejectiles, and influence on primaries and dose. Adapted from [3].
its track, the continuous slowing down approximation (CSDA) is often used.
It assumes a continuous loss of kinetic energy of a particle passing through a medium, thus neglecting any fluctuations in energy loss. The CSDA assumes that the energy loss rate at every point of a particle track equals the total stopping power, and thus leads to the equation for the stopping power as presented in section 2.2.1.
The main consequence of this type of scattering is the range in the target, since the loss of energy is continuous and not significant per single interaction.
Another important effect is the ionization of the atomic electrons, which can lead to the production of secondary electrons, also known as δ-rays, that also contribute to the dose received by the patient.
2.1.2 Elastic Coulomb scattering
Elastic scattering refers to interactions during which the internal state of the target and projectile is left intact. Elastic Coulomb scattering occurs when a charged particle passes through matter while being in close proximity to one
of the nuclei, as it is then deflected by a repulsive force from the positively charged proton(s) of that nucleus.
In order to model such interactions, the static field approximation is used.
Assumptions made in order to use this approximation include an electronic atomic density that is constant in time, which is always valid for small projectiles and small changes in direction and which does not influence the modelling of the transport. The approximation also assumes that there is no photon exchange after recoil of the target. The static field approximation reflects the Coulomb interaction by the Coulomb potential
V (r) = Z0Ze2
r φ(r) (2.1)
for a projectile at a distance r from a point-like charged nucleus, where φ(r) corresponds to the screening of the nuclear charge by atomic electrons.
Since elastic Coulomb scattering leads to a change in direction of the pro- jectile, while conserving the kinetic energy, its main influence on dosimetric quantities is the lateral penumbra, defined as the sharpness of the fall-off of the dose profile.
2.1.3 Non-elastic nuclear interactions
When an incoming ion hits the atomic nucleus, non-elastic nuclear reactions may occur. During these kind of interactions, the nucleus as well as the projectile ion are irreversibly transformed: the ion is absorbed by the nucleus and one or more other particle(s) is/are ejected from that nucleus. Secondary particles created by such ion-induced nuclear reactions include neutrons, protons, deuterons, tritons, 3He, 4He, and other ions [3] [11].
Within Pegasos, nuclear reaction data for the protons is generated using the code TALYS. The probability of a nuclear reaction is determined by the generated cross sections, and the energies of the secondaries by the emission spectra.
2.1.4 Bremsstrahlung
Bremsstrahlung occurs when the projectile charged particle interacts with the electromagnetic field of the atomic nucleus. This interaction can cause the charged particle to be slowed down and hence lose kinetic energy. Obeying the law of conservation of energy, this lost energy is converted to gamma
CHAPTER 2. PROTON AND ION TRANSPORT PHYSICS 14 radiation, which is then emitted. The deceleration by a nucleus of charge Ze on a particle of charge ze and mass M is proportional to ZzeM2, and the intensity of the x-rays scales with the square of this amplitude [12]. Thus, the intensity of the emitted x-rays varies as:
IB ∝ Z2z2
M2 . (2.2)
This implies that bremsstrahlung is an important factor for electrons be- cause of their low mass. It should therefore probably be considered for the secondary electrons generated by charged particle-induced nuclear interac- tions or Coulomb scattering. For protons at therapeutic energies, however, bremsstrahlung is in fact negligible, which is why it is not included as part of the proton nor carbon-12 ion interaction types in Pegasos.
2.2 Basic concepts
2.2.1 Stopping power
The stopping power of charged particles is defined as the average energy loss per unit path length of a particle in a material, and is comprised of three parts: the electronic, the radiative and the nuclear stopping powers [13]. Its equation is found in the CSDA approximation, and the electronic stopping power for the specific case of protons is given by [5] in accordance with [13]:
Sel= N2πZ02e4 mev2 Z
2 ln 2mec2 I
+ f (γ) − δF
(2.3) with
f (γ) = 2[ln (β2γ2) − β2] + ln R + β2
1 − R − 2γme
MR + β2γ2m2e M2R2
where N = ANAρ
mol, ρ is the mass density of the material, Amol is the molar mass, NA is Avogadro’s number, Z0e (with Z0 = +1) is the particle charge, e is the elementary charge, me is the electron mass, v is the particle velocity in units of c, c is the speed of light, I is the mean excitation energy taken from [14], δF is the Fermi density correction, β = vc, γ = q
1 1−β2, R =h
1 + mMe2
+ 2γmMei−1
, and M is the particle rest mass.
To intuitively understand this stopping power dependence, note simply that
− Sel∝z v
2
ρZ, (2.4)
so a charged particle’s average energy loss per unit path length increases linearly with density ρ and atomic number Z, quadratically with the charge z, and is inversely proportional to the particle energy (E ∝ v2). Another notation for the stopping power Sel as lost energy over length is dEdx, which is sometimes used in the following parts for better understanding.
The stopping power varies over a large range of materials, so in cases in which this is not preferable, the mass stopping power can be used. It is obtained by dividing the stopping power by the density ρ of the material:
− dE
d = −1 ρ
dE
dx = −1
ρSel= z2Z
Af (β, I). (2.5) Hence, it varies little for materials with similar ZA, so, for instance, mass stopping powers for water can be scaled by density and then used for materials containing above all light elements, like hydrocarbons or tissue.
The mass stopping power is given in units of MeVcm2/g or multiples.
Comparing −dEd of different materials reveals that, in general and per traversed distance scaled by the density, heavy atoms do not slow down the traversing charged particles as efficiently, since their inner shell atomic electrons are more tightly bound, so are not as receptive to absorbing energy.
2.2.2 Range
The range of charged particles in a medium is defined by the depth at which half the incident particles have stopped in the medium. It is an average value, since the energy losses of individual particles vary slightly due to range straggling [3]. Occasionally, charged particles that underwent nuclear interactions - and are thus absorbed in the material - are ignored, so the range is then only calculated over the number of particles which do not undergo absorption earlier in their trajectory. The range R of charged particles with initial kinetic energy E can analytically be obtained directly from stopping power calculations as follows:
R = Z E
0
dE0 dx
−1
dE0. (2.6)
CHAPTER 2. PROTON AND ION TRANSPORT PHYSICS 16
2.2.3 Cross section and mean free path
The total cross section σ corresponds to the likelihood per unit distance of travel that any interaction takes place. As such, it can be expressed as:
σ = N
Φ, (2.7)
where N equals the number of interactions per unit time per target and Φ is the flux, i.e. the number of incident particles per unit area and per unit time.
The cross sections of a certain particle in a certain medium are often also specified per type of interaction, e.g. the elastic cross section corresponds to the probability that an elastic collision occurs [15].
For a more detailed calculation, for instance within Monte Carlo codes, the differential cross section provides more information. For instance, the doubly differential cross section gives the probability that a particle with kinetic energy E undergoes an interaction while travelling over a distance l, while its direction after the interaction is within the solid angle dΩ and its energy is within dW [16]:
σ = Z dσ
dΩdΩ
= Z E
0
Z d2σ
dΩdWdΩdW.
(2.8)
If the cosine of the scattering angle θ is taken as the independent variable, the transport cross sections describing the cross sections related to the average momentum transferred in one interaction can be established. For elastic interactions, they are defined by:
σel,l = Z
[1 − Pl(cos θ)]dσel
dΩdΩ, (2.9)
where Pl are Legendre polynomials. Hence, the first transport cross section is obtained with P0(x) = 1 and is therefore given by:
σel,1= Z
[1 − cos θ]dσel dΩdΩ
= 2σel Z 1
0
1 − cos θ
2 pel 1 − cos θ 2
d 1 − cos θ 2
= 2σelD1 − cos θ 2
E ,
(2.10)
in which pel(1−cos θ2 ) is the normalised probability density function (PDF) of
1−cos θ
2 in a single collision [16] [5].
Another important quantity that can be derived from the cross section is the mean free path (MFP), which corresponds to the average distance of travel of a charged particle between two consecutive collisions. It can be derived using the assumption that a medium has a homogeneous distribution with N molecules per unit volume. Then, the interaction probability per unit path length is equal to N σ, with the probability of an interaction occurring in the interval (l, l + dl) being N σF (l), where
F = Z ∞
l
p(l0)dl0 (2.11)
is the probability that a particle travels a distance l without interaction.
From there, the PDF of l is given as p(l) = N σ
Z ∞ l
p(l0)dl0, with p(∞) = 0
=⇒ p(l) = N σe−lN σ.
(2.12)
Hence, the MFP equals the expectation value of path length l:
λ = hli
= Z ∞
0
lp(l)dl
= 1 N σ,
(2.13)
and the first transport mean free path (for e.g. elastic interactions) is [4]
λel,1 = 1
N σel,1. (2.14)
Hence, the inverse mean free path (IMFP) can be used as a quantity which is linearly proportional to the interaction probability:
λ−1 = N σ. (2.15)
Chapter 3
Physical Quantities in Radiotherapy
3.1 Fluence
According to the International Commission on Radiological Units and Mea- surements (ICRU), the fluence corresponds to the number of particles traversing a sampling sphere [17], in units of m−2. A such, the fluence is often defined by the number of particles divided by the cross-sectional area of that sphere. It can also be calculated per unit time, e.g. as in a run.
Within Monte Carlo systems, a slightly altered definition of fluence is used, which has also been published by the ICRU: the length of the particle track steps within the sampling volume (of any shape) are utilized in order to calculate the fluence. For instance, Salvat defines the fluence Φ(r) at point r as
Φ(r) ≡ dNin
dA , (3.1)
where dNin is the total number of particles hitting a small sphere of cross- sectional area dA, centred at r [4]. At the same time, the hypothesis holds that the average fluence in a finite volume V is given by:
Φ = 1 V
X
i
(path length of particle i in V). (3.2)
18
3.2 Dose
The absorbed dose, or commonly (and therefore in this thesis) just referred to as dose, is defined by the energy E deposited by radiation particles in matter per unit mass M , and is given in units of Gy (gray) equivalent to the SI base units of J kg−1:
D = E
M. (3.3)
The absorbed dose is meant to be a deterministic quantity that reflects the severity of direct tissue damage caused by radiation.
The equivalent dose takes into account the effect of different radiation types by incorporating weighting factors WR for the radiation type at hand R:
H =X
R
WR· D. (3.4)
It is is given in Sv (sievert), also equivalent to J kg−1; note that the Gy is used for physical quantities, whereas the Sv is used for biological effects.
Then, the effective dose includes the type of tissue or organ upon which the radiation is impending by taking the tissue-weighted sum of the equivalent doses:
E =X
T
WT · H. (3.5)
Its unit is also Sv and is a stochastic measure of the overall health risk to the whole body due to radiation. The official dose limits, which in Sweden are set by the Swedish Radiation Safety Authority, are effective dose values.
3.3 Linear Energy Transfer
The first definition of the linear energy transfer (LET) was given by the ICRU in 1962 [18] and stated that the quantity reflected the average energy locally imparted to the medium by a charged particle of specified energy in traversing a [certain] distance. Since the term ’locally’ referred to values below either a maximum distance or a maximum energy, it was a rather imprecise definition.
The same report states that the difference between stopping power and LET is that the former disregards information as to where the energy is absorbed, whereas the latter treats the energy lost within a volume. Over the following years up to 1970, the need for a preciser LET definition was fulfilled by
CHAPTER 3. PHYSICAL QUANTITIES IN RADIOTHERAPY 20 specifying it to include a maximum energy (and not a distance) cutoff value, which was denoted by ∆ [19]. Nowadays’ revised definition of the LET includes this idea; it is given by the ICRU in 2011 and states that [20]:
The linear energy transfer or restricted linear electronic stopping power, L∆, of a material, for charged particles of a given type and energy, is the quotient of dE∆ by dl, where dE∆ is the mean energy lost by the charged particles due to electronic interactions in traversing a distance dl, minus the mean sum of the kinetic energies in excess of ∆ of all the electrons released by the charged particles, thus
L∆ = dE∆
dl . Unit: J m−1.
This means that the LET can also be given by:
LET∆= L∆= Sel− dEke,∆
dl , (3.6)
where Sel is the electronic stopping power previously given in equation (2.3), and dEke,∆ denotes the secondary electron (δ-rays) kinetic energies which are higher than ∆.
From this last equation follows a quantity called the unrestricted LET, obtained by taking the limit ∆ → ∞. All δ-rays are then included and we have:
LET∞= Sel. (3.7)
3.3.1 Track- and dose-averaged LET
This LET∞ value thus equals Sel and is easily specified for a monoenergetic ion beam. However, in a realistic clinical beam used for proton or carbon ion therapy, there are often different initial particle energies with an energy spread, as well as energy straggling, which leads to an LET value that differs per point in the matter. As such, there is a need to average these LET values in order to get a distribution of LET per small volume. There are currently two different approaches: the track-averaged LET (LETt) and the dose-averaged LET (LETd). The difference between the LETt and LETd is
the weighting of the averaged Sel of the primary charged particles. The LETt weights the averaged Sel by particle tracks, or fluence, whereas LETdweights it by dose, or local energy transfer contributions by electronic interactions.
Track- and dose-averaged LET are commonly given by [21] [22] [23] [24]:
LETt= R∞
0 Sel(E)Φ(E, z)dE R∞
0 Φ(E, z)dE , (3.8)
LETd = R∞
0 Sel(E)D(E, z)dE R∞
0 D(E, z)dE , (3.9)
where Φ(E, z) is the fluence (or the particle spectrum) of primary charged particles with kinetic energy E at location z, and D(E, z) the corresponding deposited dose. The latter calculated per bin can be approximated using the assumption of energy loss with the continuously slowing down approximation, and therefore the stopping power and fluence [25]:
D(E, z) = SelΦ(E, z)
ρ(z) , (3.10)
where ρ(z) is the density at z. Note that the assumption used neglects the escape of delta-rays and the nuclear stopping power. Thus, from equation 3.9 we have for the dose-averaged LET:
LETd= R∞
0 Sel2(E)Φ(E, z)dE R∞
0 Sel(E)Φ(E, z)dE. (3.11) The dose-averaged LET is expected to be more closely related to the RBE [22]
[26] [27] [28], which is why mostly the LETd will be discussed in this thesis.
Conventional LET calculation methods specifically used within Monte Carlo codes are given later in section 4.4.
3.4 Relative Biological Effectiveness
The Relative Biological Effectiveness (RBE) is defined as the ratio of doses to reach the same endpoint X (i.e. the same level of a specific effect) when comparing two types of radiation [27]. Typical endpoints include cell killing,
CHAPTER 3. PHYSICAL QUANTITIES IN RADIOTHERAPY 22 double strand breaks (DBS), cell mutations, chromosome aberrations, and foci formation. For the RBE of protons:
RBE = Dreference(X)
Dprotons(X) , (3.12)
where the reference radiation is most often chosen to be 60Co photon radiation. In vivo RBE values from the first proton therapy measurements were found to have an average RBE of 1.1 relative to such high energy photons, which is why this value is taken as a generic clinical value. Thus, protons are said to be 10% more effective than 60Co gamma rays. This assumption disregards that the RBE depends on numerous other factors, including
• Macroscopic dose D,
• Proton energy distribution,
• Endpoint X,
• Proton beam properties (such as range or modulation width),
• Tissue and cell type.
It is therefore unsurprising that researchers have proven this generic RBE value of 1.1 to be unrepresentative, or in fact even an underestimation, of the actual clinical RBE. For instance, in vitro RBE values relative to the same 60Co radiation were found to be 1.2 [29]. Regarding carbon ions, their RBE variation within each treatment field is considered too significant to even set a generic value; hence, their RBE has to be calculated separately for each case.
Note that this concept of RBE does not hold on the microscopic scale, since there then are inhomogeneities in the energy distribution of the individual particle tracks. Moreover, D is assumed to be homogeneous in the chosen volume, which does not hold for real dose distribution in e.g.
organs.
Regarding the dependence of the RBE on the cell and tissue type, the linear-quadratic (LQ) model is the most commonly used model, since it repro- duces the experimental data fairly well. It follows from the characterisation of the biological system by parameters α and β [30], and the fact that the cell killing by a dose D is given by the DNA double strand break (DSB) yield Y = αD + βD2 [31]. As such, the LQ model corresponds to finding α and β to fit
N N0
surv
= e−αD−βD2, (3.13)
where
N N0
surv is the cell survival fraction [32]. It corresponds to the probability that a cell survives after being exposed to a dose of radiation and follows the hypothesis that for a same dose, the number of cancerous cells is decreased logarithmically (‘log cell kill’). The model includes two tissue- dependent parameters: α, which increases linearly with the dose and includes all biological repair mechanisms that are independent of the dose rate, and β, which increases with the dose squared, and takes into account dose rate dependent mechanisms, which are so-called sub-lethal lesions that originate from misrepairs through mutual cell interaction. The α/β ratio is hence considered to represent the cells’ radiosensitivity. Some critics of the model question whether it truly represents the underlying biological mechanisms or whether it simply is a useful empirical fit. Other drawbacks include that it is limited to a dose range of 1-10 Gy, so is not applicable to lower doses (such as for organs at risk), nor to higher doses (sometimes needed for ocular tumour treatments). The model also excludes heavier particles like alpha particles. Finally, the RBE calculated by the LQ model is very dependent in small variations of the reference beam properties, which vary substantially depending on the source [23].
The relation between the RBE and the LET is not clearly defined either, but it is known that the RBE increases with the LET (mostly because of increasing α, as β stays constant up to an LET of around 20 keV/µm, if one uses the LQ model) [23]. Biologically this is expected as the repair mechanisms of the damaged cells decreases with increasing density of ionisations. However, the RBE decreases with increasing LET after a LET value of around 100 keV/µm due to saturation effects of the energy deposition events. In fact, this value of 100 keV/µm corresponds to one ionisation event per 2 nm, a length which equals the diameter of a DNA strand, so is therefore considered as the optimal LET for the killing of cells [33]. Above this value, the RBE decreases with increasing LET due to the overkill effect, during which one particle deposits much more energy than is necessary to kill a cell, so less cells are killed per absorbed dose [33].
An important point is also that the RBE differs for different particles of same LET, due to a difference in the number of produced complex lesions and due to different dose patterns [34]. Note that such lesions are mainly caused by primary particles for high LET, and by secondaries for low LET.
Moreover, particles of low LET mostly act through indirect action, during which the radiation hits the water and organic molecules in the cells and
CHAPTER 3. PHYSICAL QUANTITIES IN RADIOTHERAPY 24 thereby releases highly reactive radicals, such as hydroxyl (HO−) and alkoxy (RO−2). On the other hand, high LET particles mainly act through direct action, during which radiation hits the DNA strands directly, causing breaks in the molecular structure, which can lead to cell death or abnormalities formed during cell repair mechanisms [35].
3.5 Lineal energy
The lineal energy is a quantity widely used within microdosimetry and is defined as the energy s imparted to the medium in a given volume by a single energy deposition event, divided by the mean chord length l of the volume [36] [37]:
y = s
l . (3.14)
The mean chord length can be calculated as l = 4V
A , (3.15)
where V is the volume and A is the surface area of the given volume. The latter is most often defined as a sphere, where then l = 23d.
Since y is defined as a stochastic quantity, it can be sampled from a probability density function f (y), which in turn can be derived after multiple particle tracks have passed through the volume of interest. In essence, the mean-frequency lineal energy can be calculated by [38]
yF = Z ∞
0
yf (y)dy. (3.16)
It follows that the dose-averaged lineal energy, which is often used with relation to the biological impact of the particles, and with the dose probability density d(y), can be found with
yD = Z ∞
0
yd(y)dy. (3.17)
As such, the lineal energy y differs from the LET in quite some aspects. First of all, since y is a statistical distribution, it is measurable, for instance using a low pressure Tissue Equivalent Proportional Counter (TEPC), in which
a microscopic tissue is simulated by a macroscopic detector cavity loaded with a tissue equivalent gas at low pressure [39]. The LET, on the other hand, is merely calculable. Generally, charged particles are characterised by an energy range, which can be represented by the LET with y distributions.
Another way of looking at it is that the LET focuses on the energy lost by a particle, while y deals with the deposited energy per specified volume [36]. Furthermore, the energy deposition for the LET are calculated by only including the primary particle, while s includes the total energy deposited by the primary as well as all of its generated secondaries, within a single deposition event. Finally, the scored length for the LET is defined as the true step length of the primary particle, while l in y is a mean length which is constant for any particle type or any event [21]. All in all, ion beams with an identical LET have a much larger yD value for light particles, as compared to heavier ones.
In addition, multiple recent studies suggest that the lineal energy could be a means to unveil the physical mechanisms that govern the RBE dependence on the LET. For instance, it has been found that the rate of yD calculated for a nanometer-scale sphere roughly equals the from cell survival curves derived α coefficients of the LQ model [37] [40].
Chapter 4
Monte Carlo Methods for Charged Particle Transport
4.1 Monte Carlo techniques
Monte Carlo simulations repeatedly sample from random probability distri- butions in order to generate numerical results. The method was first sug- gested by Metropolis around 1940 as an approach to resolve neutron diffusion and multiplication problems [16]. Metropolis and Ulam [41] characterised Monte Carlo methods as
a statistical approach to the study of differential equations, or more generally, of integro-differential equations that occur in various branches of the natural sciences.
Monte Carlo techniques use statistical sampling of algorithm generated distributions, as well as statistical inference as a means to find an estimate (e.g. mean) of a parameter of a distribution, corresponding to the solution to a mathematical problem. This use of statistics, added to any possible algorithm errors, implies that solutions always have a statistical uncertainty, defined by the width of the used distribution. Since particle transport is inherently random, Monte Carlo techniques constitute a self-evident choice of method [16].
Within charged particle simulations, there are in general three categories of computational Monte Carlo techniques used: Detailed, Condensed, and
26
Technique Characteristics (Dis-)advantages Detailed MC Modelling of each collision. 3High precision
7Computationally very expensive Condensed
MC
Energy loss and angular scattering are sampled from a probability dis- tribution.
3High computational speed 7Geometry handling problems
Mixed MC
Use of hard and soft events, as above and below an energy and angle threshold.
3Good computational speed 3Good precision
Table 4.1: Different Monte Carlo (MC) techniques used for charged particle transport simulations.
Mixed Monte Carlo. A summary of their characteristics, as well as advan- tages and disadvantages is given in Table 4.1. The method used at Elekta and therefore for this thesis is a Mixed Monte Carlo method, meaning that for the interactions between the charged particles and matter, a distinction is made between hard events and soft events. Hard events include all inter- actions with an energy deposition larger than a set cutoff energy value, as well as a scattering angle larger than a cutoff angle, whereas soft events cor- respond to all other interactions. Energy depositions and scattering angles of hard events are computed in detail. For soft events, energies and angles are sampled from a probability distribution and deposited at a hinge at a random position between previous and current event.
4.2 Code systems
4.2.1 PENELOPE
PENELOPE (Penetration and ENErgy LOss of Positrons and Electrons) is a general-purpose Mixed Monte Carlo code system. It was developed by Francesc Salvat and collaborators at Universitat de Barcelona, and the code was first published in 1995 [42], while the currently used version was released in 2014 [4]. The code is at present maintained by OECD Nuclear Energy Agency. PENELOPE is written in FORTRAN 77, and meant to be accompanied by a user written master program. The code is written for electron, positron and photon transport for energies of some eV up to 1 GeV.
Electron and positron transport mechanisms are true to the Mixed Monte
CHAPTER 4. CHARGED PARTICLE MONTE CARLO 28 Carlo nature, whereas photon interactions cannot really be divided into hard and soft events and are therefore simulated in detail.
4.2.2 PENH
PENH is the PENELOPE extension written in order to include proton trans- port, and was also written by Salvat [5] and published in 2013. Electromag- netic interactions and generalities are described in [5], whereas [43] explains the integration of some nuclear reactions into the code. PENH can include energies up to tens of GeV, and the proton transport mechanisms mimic the ones for electrons and positrons, with a few distinctions. Elastic scattering is desribed by the eikonal scattering model, whereas inelastic scattering is handled by the Sternheimer-Liljequist model.
4.2.3 Geant4
Geant4 (GEometry ANd Tracking) was initially developed independently at CERN and KEK and is another toolkit that uses Monte Carlo methods in order to simulate the passage of particles through matter. It is written in C++ and includes a wide range of physical processes such as electromagnetic, hadronic and optical ones. Details concerning the physics, particle types, energy ranges, etc. are explained in [44], and tips on how to use and implement the code are given in [45].
4.3 Pegasos
Elekta’s internal Monte Carlo code system is called Pegasos and is a modified version of PENELOPE. It has been extended with PENH (then called Pegasos-PENH ); both PENELOPE and PENH are described above. The Pegasos-PENH code handles the combination of electron, positron, photon and proton transport, while Ionos also includes light ion transport for ions up to carbon-12. Pegasos is used for simulation of Elekta products and prototypes, such as the Gamma knife or the MR linac (MRI guided linear accelerator), so all in all for dose calculation, hardware design and radiation shielding studies. Pegasos includes subroutines written in C++ that allow for parallel computing. These subroutines include scoring modules that allow for calculations of e.g. dose and fluence. A graphical user interface
(GUI) called Hermes has also been developed at Elekta in Java and serves to create input data for Pegasos that make the use of complex CAD geometries possible. The user can choose how many nodes to run the code on. If more nodes are chosen, the Hermes application calls PegasosMPI, which is a multiprocessing program using the Message Passing Interface (MPI). The system automatically combines the results (e.g. from all nodes) into a .tot text file that can readily be handled by e.g. MATLAB for postprocessing.
4.3.1 Input
Hermes takes its input from a simulation file, which contains numerous adaptable parameters, and builds material and geometry data files for the simulation in Pegasos. Input parameters are divided into subcategories, such as beam properties, like type of source (point source, body source, part source, parallel beam, divergent beam, phase space, moving source), particle type, beam shape (circular, elliptic, rectangular), spectrum type (monoenergetic, Gaussian, discrete, Co-60, histogram), beam radius, number of particles, source position, and beam direction. As to the geometry, primitive geometries can be created or more complex ones imported; for each of them, the material can be specified and material files imported from the
’pendbase’ data base. Variables like cutoff energies and absorption energies are specified per material and per particle type. Scoring input parameters correspond to geometry and size of the scoring volume, number of bins, and of course scoring module to be used.
4.3.2 Simulation
Pegasos first takes the files generated by this user input in order to simulate the particle interactions. Each simulation shower (i.e. the simulation of one particle) mimics the real particle transport mechanisms and the physics behind them by considering the particles to roughly consecutively undergo a step and then an interaction event, another step and another interaction event, and so on1. More precisely, one whole shower for each of the primary particles from the source follows the following steps, in which important subroutines are given in capitals:
1The design of the tracking algorithm was implemented such that the effect of the geometry on the transport physics in minimized [4], it thus only requires knowledge of the particle state and the material in which the particle currently is.
CHAPTER 4. CHARGED PARTICLE MONTE CARLO 30 1. Particle generation.
2. Location and scoring of particle.
3. Cleaning of the secondary particle stack.
4. Simulation of particle transport.
(a) Start: Start of the simulation in the current medium.
(b) Jump: Determination of segment length.
(c) Step: Moving of particle to end of step length and scoring of fluences.
(d) Check if a material border was crossed. If yes, go to (a).
(e) Check if particle left the simulation region. If yes, go to (h).
(f) Knock: Simulation of particle interaction and scoring of energy deposition.
(g) Check if particle was absorbed. If no, go to (b).
(h) Check if secondary particles were generated. If yes, popping of secondary stack and start of secondary particle simulation by going to (a).
5. End shower.
In each Knock event, the particle loses energy, changes its direction of movement and possibly produces secondary particles. Moreover, true to the Mixed Monte Carlo nature of the code, the Knock subroutine differentiates between soft and hard events. As such, inelastic and elastic multiple Coulomb collisions (each below an energy deposition and scattering angle cutoff value, as mentioned before) are simulated by sampling from a probability distribution and therefore called artificial soft events. A more detailed explanation follows in section 4.3.3.
Another type of interaction that needs specification is the delta- interaction, during which no energy deposition occurs, but it simply acts as an empty Knock call with no alteration of the physical state of the parti- cle. The delta-interaction works as a means to limit step lengths. Currently, there are in total nine types of interaction: artificial soft events, hard elas- tic collisions, hard inelastic collisions, bremsstrahlung emissions, ionisations of inner shells, absorptions, delta-interactions, elastic nuclear events, and inelastic nuclear events.
As to particle types, of interest in this case above all in terms of sec- ondary particle generation, Proteus distinguishes between electrons, photons, positrons, protons, neutrons, deuterons, tritons, Helium-3 particles, and alpha particles, while Ionos also handles lithium-6, lithium-7, boron-10,
boron-11, carbon-11 and carbon-12. There are of course different data bases and interaction algorithms for each of these.
4.3.3 Mixed Monte Carlo
As introduced in section 4.1, Pegasos-PENH is a Mixed Monte Carlo code, meaning that the distinction is made between hard and soft events. It is a method exclusively used in charged particle simulations, making use of the fact that small scattering angle interactions are more probable and frequent than large angle interactions. For e.g. photons or neutrons, the ’condensed’
interactions are not needed, since the probability of an interaction occurring is not related to the deflection angle.
Thus, in order to find a good trade-off between accuracy and speed, the following simulation parameters are used in Pegasos to distinguish between hard and soft events [4]:
• C1H: the average angular deflection over various soft elastic scattering events along a distance corresponding to the mean free path between consecutive hard elastic interactions;
• C2H: the maximum average fractional energy loss by soft events between consecutive hard elastic interactions;
• WCCH: the cutoff energy loss for hard inelastic interactions.
The smaller C1H and C2H are, the more accurate - but slow - the simulation is. More precisely, the hard elastic mean free path λ(h)el is given in [4] by:
λ(h)el (E) = max
λel(E), min
C1H · λel,1(E), C2H · E S(E)
, (4.1)
where S(E) is the total stopping power at particle energy E, λel the elastic mean free path, and λel,1the first transport mean free path given in equation 2.14. Essentially, when the energy E increases, λel approaches a constant value, while λel,1 increases too; hence, the MFP for hard elastic events in equation 4.1 increases with increasing E. In this way, hard interactions are more frequent when the scattering is stronger. Moreover, for low E, λ(h)el = λel, and the simulation is entirely detailed.
As a cutoff energy, the variable WCCH mostly has an impact on the energy distributions in the simulation; if its value is high, the simulation is fast, but the distributions may be distorted.
CHAPTER 4. CHARGED PARTICLE MONTE CARLO 32
4.4 Linear energy transfer calculations
This section serves as an overview over the LET calculation methods most commonly used within Monte Carlo codes in literature. Note that even though most sources cited below make use of Geant4, the methods are of course applicable to most other Monte Carlo code systems.
Calculations of LET within Monte Carlo systems nearly always use the unrestricted LET (LET∞). Most authors do not give an implicit reason for this, but state that is common practice. Only Guan et al. (2015) claim that smoother LET curves are obtained when LET∞ is used, because single δ-rays cause spikes [21]. They however set a cutoff value for the production of δ-rays, but make sure that it is larger than the range of the δ-rays in order to locally include them in the calculation. Grassberger and Paganetti (2011) indicate that the use of LET∞ is only valid if charged particle equilibrium exists, which is not true in the field penumbra or at interfaces [11].
4.4.1 Track-averaged LET
Since the LETt weights Sel by the fluence (see section 3.3.1) and the fluence within Monte Carlo methods corresponds to the path length li of each particle i scaled by the bin volume (see section 3.1), the LETt can be calculated by [21]:
LETt= P
iωi· LETi· li P
iωi· li
,
where ωi is the statistical weight attributed to the primary particle, and LETi = li
i is calculated as the mean energy loss per unit path length according to the charged particle kinetic energy at each step, divided by the step length li of particle i. As such, for the nominator, we have li· LETi = li· li
i = i. Hence:
LETt = P
iωii P
iωili
, (4.2)
and LETt is analogous to summing all energies over all travelled steps for each particle in each bin [22] [46].
4.4.2 Dose-averaged LET
Regarding the LETd calculation methods, the first, most intuitive one stems from taking into account the general definition of the dose-weighted LET given in equation 3.11, as well as the fact that for the unrestricted LET we have Sel,i= LETi:
LETd = P
iωiLETiDi P
iωiDi , (4.3)
for each primary particle i (in literature also called beamlet or field i), and where ωi is again the statistical weight of the primary particle. This method is directly used by [26], [47] and [48].
As presented in section 3.2, the definition of the dose within radiotherapy physics is given by the sum of all deposited energies s for each step s over the mass M of recipient: Di =
P
ss
M = Mi. Since M is constant for one recipient, it then follows that:
LETd= P
iωiLETii P
iωii , (4.4)
where, again, LETi is the mean energy loss per unit length according to the particle energy, divided by the step length. Note that LETi can also be extracted or interpolated from stopping power tables given the energy of the particle at the start of each step. The above relation is used in [21], [24] and [46].
If one then simplifies the linear energy transfer LETi of each particle by its literal definition of deposited energy per path length, LETi = li
i, then:
LETd= P
iωil2i
i
P
iωii, (4.5)
where i is the energy deposited by the primary particle i plus the kinetic energy of released δ-rays; as in [24], [46] and [49]. There is some criticism that this method is highly dependent on the secondary electron cutoff, but this is mainly due to the specific Monte Carlo code (e.g. Geant4), in which particle step lengths can be shortened greatly due to voxel size limitations.
Finally, the last method that is commonly used in literature is obtained by taking a different definition of dose, including stopping power Si, fluence (or number of particles) φi, and density ρi to get: Di = Sρiφi
i . Then:
LETd= P
iSi2φi(z) P
iSiφi(z). (4.6)
CHAPTER 4. CHARGED PARTICLE MONTE CARLO 34 Most papers, e.g. [22] and [46], calculate φi and get Si from available tables.
Since the LETd scoring corresponds to squaring the energy deposition and then dividing that by the path length, it is rather sensitive to very high energy depositions, as well as very small path lengths. This is why some authors impose a minimum scored length before the LETd is scored, as for instance in [24], all scored lengths which are smaller than 2 nm are discarded.
This will avoid numerical instability and unnatural spikes in the LETdvalues.
Methods
In this chapter, the methods used in order to generate the results presented in chapter 6 are laid out. Since the focus lies on protons, most choices and justifications are given for protons, while the ones for carbon ion are briefly given in section 5.2. Finally, the mathematical derivations for the statistical LET evaluations are put forward.
5.1 Proton LET
5.1.1 Code implementation
Scoring modules for the calculations of the LET are written in C++ and added to the existing Pegasos code. The modules consist of one main class called ’LETBox’ which includes the functions and variables necessary to calculate the LET within rectangular scoring boxes, and two subclasses for calculations of the LETt and LETd specifically.
5.1.2 Particle generations and particle types
This subsection serves as an identification and explanation of the idea and methods behind the scoring of LET, both in terms of particle types and particle generations.
By the ICRU publications concerning LET specified in section 3.3 (e.g.
[19], [18] and [20]), and following most papers published on the computation of LET, the LET is understood as a physical quantity which differs per particle type. The idea is to get a measure of each type of radiation and its
35
CHAPTER 5. METHODS 36 consequences on e.g. the RBE and thus the biological damage it causes in a patient’s targeted body volume. Hence, the concept is often only directly logically applicable to one kind of particles per calculation, i.e. the primary particles at hand. As such, the LET module’s default settings only track and score the primary particles.
Even though it has been established that the LET of the primaries is the main goal, these primaries still generate numerous secondaries when they interact with the medium they pass through. Such secondaries include many kinds of particles, of which electrons, photons, positrons, protons, tritons, deuterons, Helium-3 particles, and alpha particles can be specifically tracked in proton simulations in Pegasos. Note that for the carbon-12 simulations, lithium-6, lithium-7, boron-10, boron-11, carbon-11 and carbon-12 are added to these particle types.
The LET of secondaries of each specific type is calculated separately, so both particle generation and particle type can be specified by the user with input flags. Since, however, a clinical proton (or carbon ion) beam obviously does not only include the primary particles, but also numerous secondaries, the overall LET for all particle generations is computed as well.
Regarding which types of secondaries should be accounted for, the definition of LET given above only specifies how electrons should be included in its calculation: if their kinetic energy exceeds the cutoff, it should be subtracted from the energy deposited locally by the primaries. However, Monte Carlo calculations in literature mainly address the unrestricted LET, in which δ-rays are included locally in the LET. As explained before, Pegasos is set up in a way to calculate the energy deposited locally by each particle in each Knock event, and automatically subtracts any energy that is carried away by possible secondaries, since each of them is simulated, too. Hence, the written LET module adds the kinetic energy of any generated electron to the energy scored in the bin in which it was generated, and cancels its further simulation.
As to the other types of secondary particles, the LET concept does not make sense in relation to photons, since their physical state and path remains unchanged until a possible single event of absorption (e.g. Compton effect, photoelectric effect or pair production). Note that when some authors still refer to a photon LET, the LET calculated from the δ-rays generated by the photons is meant. In the same way, photons are still tracked in Pegasos- PENH, but only the other types of secondaries that originate from their interactions are directly scored.
Particle type Composition Mass [u] Charge
Proton n/a 1,0073 +1
Neutron n/a 1.0087 0
Deuteron 1 p+, 1 n0 2.0141 +1 Triton 1 p+, 2 n0 3.0160 +1 Helium-3 2 p+, 1 n0 3.0160 +2 Alpha 2 p+, 2 n0 4.0015 +2
Table 5.1: Types of secondary particles, their compositions in terms of number of protons p+ and neutrons n0, their masses, and their charges.
As a conclusion, for proton therapy it only makes sense to calculate the secondary particle LET for hadrons and heavier ions: protons, neutrons, deuterons, tritons, Helium-3 particles, and alpha particles, which all have been generated by the primary proton beam of interest interacting with the medium it passes through. An overview of their respective compositions, masses and charges is given in Table 5.1.
5.1.3 Types of interaction
Since, by the latest ICRU definition [20], the LET includes only the energy losses originating from the electronic interactions, all nuclear reactions should be excluded from the standard LET calculations. Therefore, the default settings of the LET module do not include them, but the user can manually change this setting with a boolean input flag.
However, an analysis of the LET including nuclear reactions is included in this thesis, since their energy deposition cannot be disregarded in reality.
Though the nuclear events are expected to be seldom, even their individual energy contribution are expected to give rise to large LET values.
5.1.4 Choice of scoring method
This subsection serves as an explanation of the scoring method that is used for any particle generation, and possibly particle type, of interest. However, in order to stay comprehensible, detailed explanations are only given for the main LET calculation of primary protons, but can of course be generalised to other (secondary) particles.
