Engineering, 2013, 5, 12-17
http://dx.doi.org/10.4236/eng.2013.510B003 Published Online October 2013 (http://www.scirp.org/journal/eng)
Copyright © 2013 SciRes. ENG
Monte Carlo Simulation and a Review of the Physics of the
Positron Annihilation Process in PET
Hamed Hamid Muhammed, Ziya Zeng in
1Division of Informatics, Logistics and Management, School of Technology and Health STH,
Royal Institute of Technology KTH, Stockholm, Sweden
2Department of Biophysics, Faculty of M edicine, Karadeniz Technical Uni ve r sity, Trabzon, Turkey
Email: Hamed.Muhammed@sth.kth.se, zzengin@ktu.edu.tr, ziyazeng@kth.se
Received September 2012
ABSTRACT
In this paper, we investigate the physics of the positron annihilation process, which o ccurs in a PET imaging syste m. In
particular, the diffusion of beta particles (positrons) w ithin water was addressed. Beta particles are emitted isotropically
from the same source point with random directions and randomly chosen energy levels. After traversing a certain dis-
tance within water, beta particles lose a certain amount of its energy. The inelastic collisions with atomic electrons are
mainly responsible for the energy dissipation of charged particles, s uch as electrons and positrons (that have low mass).
The energy loss rate due to inelastic process is estimated by using the Beta-Bloch formula. These results help in under-
standing how to develop and implement a computationally efficient Monte Carlo Simulation of the positron annihilation
process.
Keywords: Positron Annihilation Physics; Positron Emission Tomography; Monte Carlo Simulation
1. Introduction
Positron Emission Tomography (PET) is a medical im-
aging technique used in clinical and preclinical investi-
gations to produce a three dimensional image that can
help to show how certain organs perform their activities
(Rocha and Harbert, 1978 [25]); e.g. in clinical oncology
to detect different types of cancerous tumor tissues
(Strauss and Conti, 1991 [26]; Contractor et al., 2012
[8]). The radionuclides that are commonly used in posi-
tron PET, for studies and examinations performed on
humans (i.e. patients) are characterized by: 1) Their short
half-life time, such as Oxygen-15 (~2 min), Nitrogen-13
(~10 min), Carbon-11 (~20 min), and Fluorine-18 (~110
min), to be able to minimize the radiation dose to the
patient. 2) Their ideal decay characteristics. 3) Their abi-
lity to be prepared with very high specific activities. 4)
Their chemical stability (Rocha and Harbert, 1978 [25];
Murray and Williams, 1958 [20]; Ott, 1981 [21]).
Monte Carlo simulation is an important and popular
approach which plays a key role in the research field of
PET imaging (Thompson et al., 1992 [27]; Buvat and
Castiglioni, 2002 [3]; Buvat and Lazaro, 2006 [4]; Barba
et al., 1995 [1]; Chow, 2012 [7]; Espana et al., 2009 [10];
Levin et al., 1995 [17]; Thomson and Kawrakow, 2011
[27]; Y a ma mo t o , 2010 [29 ]), b ecause of the high costs of
the sophisticated hardware equipment used in the PET
imaging system (especially the cyclotrons) in addition to
the high costs of performing real-life experiments, in-
cluding the costs related to the needed materials, labora-
tories and patient facilities as well as all the staff who
should be involved in the different phases of the experi-
mental tasks (performed on humans or phantoms).
According to the number of publications obtained
from Google Scholar for the time period beginning from
1990 and ending in 2011 (between 1990 and 201 1), the
interest in using the Monte Carlo method for the simula-
tion of the PET imaging system increased gradually dur-
ing this time period, as illustrated in Table 1.
Monte Carlo Simulation became a very powerful tool
that is used to achieve better understanding of the impact
of each part of the PET imaging system and to gain de-
tailed information about the performance of the current
design of the simulated PET imaging system and how to
improve the currently simulated design further. And this
is the main advantage of using Monte Carlo Simulation,
which makes it easy to optimize the different system
parts and parameters with high flexibility.
However, the latter issue also has some disadvantages,
mainly due to the high complexity of an actual PET im-
aging system. The consequence of this issue is to try to
simplify and include only a subset of the parameters of
the actual system in the simulation to achieve a result
which is as close as possible to the actual system. Usually,
the most important parameters, which have significant
H. H. MUHAMMED, Z. ZENGIN
Copyright © 2013 SciRes. ENG
13
Table 1. Number of publication s ob t ained from Google scholar (Totally 248 publications between 1990 and 2011).
Year 1990
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
2001
2002
2003
2004
2005
2006
2007
2008
2009
2010
2011
Number of
publications
2
4
7
4
2
3
5
3
8
9
6
7
11
13
19
20
35
30
19
14
27
27
effects on the performance and efficiency of the resulting
PET imaging system, are taken into account.
This paper addresses one of the most important para-
meters of the PET imaging system. Currently, the posi-
tron annihilation process is still playing the most signifi-
cant role in limiting the spatial resolution of the PET
imaging system (Cal-Gonzalez et al., 2009 [5]; Levin
and Hoffman, 1999 [18]; Jodal et al., 2012 [13]; Paga-
netti, 2012 [22]; Peng and Levin, 2012 [23]).This process
begins with the decay process of positron energy which
ends with the positron-electron annihilation process,
which in turn results in producing and emitting two
511-keV gamma photons (in opposite directions) that can
easily penetrate the body of the patient. Each such pair of
gamma photons are detected by the sensors of the PET
camera which surrounds certain part of the patient’s body
where the annihilation process occurred (Rocha and
Harbert, 1978 [25]). Therefor e , it is important to under-
stand the physics behind these processes to be able to
find solutions to this open problem within the field of
PET imaging.
In this paper, the physical laws behind the annihilation
process are investigated and analyzed to be able to un-
derstand how to perform a Monte Carlo simulation of
this process in a computationally efficient way.
2. The Physics of the Annihilation Process
2.1. Radiative Energy Dis sipation Effects on
Positron Trajectory Calculation
Two main types of interactions exist. The first one is
inelastic collisions with atomic electrons of the material.
The other one is elastic scattering from nuclei in without
any radiation or excitation. For the range of energies of
beta (β) particles of our interest (between 104 - 107 eV),
almost all the betas are deflected because of elastic scat-
tering from nuclei (Evan, 1955 [11]; Meyerhof, 1967
[19]). Additionally, there exist less important intera ctions
which are called Emission of Cherenkov Radiation,
Nuclear Reactions and Bremsstrahlung. In practice, we
can ignore the Bremsstrahlung because the resulting
energies are out of the range of interest. The Cherenkov
Radiation is emitted only if the speed v of a charged par-
ticle, such as an electron or a positron which traverses a
dielectric medium, is greater than the speed of light c. In
our study the speed of the β-particle is usually not greater
than c, therefore th e contribution of the Cherenkov Radi-
ation is negligible.
There are two types of elastic atomic collisions: soft
and hard collisions. In a soft collision corresponding to a
not-so-close encounter
( )
ba
, which is responsible
for only the excitation of an atom, the incoming par-
ticle’s trajectory almost never deviate from its original
track. The energy transferred by an excitation is dissi-
pated in the emission of low-energy radiation of frequen-
cies of 1018 - 1021 Hz as well as in atomic and molecular
vibrations. While in a hard collision, the incoming par-
ticle’s traj ectory is changed considerably. When this kind
of process takes place, there occurs a large energy trans-
fer to an atomic electron which is called a recoiling ioni-
zation electron. If the transferred energy is enough to
produce its own secondary ionization, this recoiling elec-
tron is called a knock-on ora delta electron (δ-electron).
2.2. Energy Dissipation by Inelastic Collision
The inelastic collisions with atomic electrons are mainly
responsible for the energy dissipation of charged par-
ticles of low mass, such as electrons and positrons. As a
traversing particle moves from one point to another in
matter, its energy will be decreased through different
interactions with the atomic electrons that it encounters.
The energy transferred from the particle to the atomic
electron will cause ionization and excitation. The rate, at
which this energy loss occurs, depends on several para-
meters, su ch as the type and energy of these particles and
the atomic-number density of the medium through which
the particle is traversing. However, at a high primary
energy level, excitation is much more probable to occur
than ionization. Ionization and excitation processes can
occur in all atomic shells, but mostly in the inner ones. If
an inner-shell electron is ejected, then filling this vacancy
will be accompanied by an emission of an electron from
this atom.
Usually, at most 50% of the kinetic energy of the pri-
mary particle can be transferred to a free atomic electron.
The average energy lost per distance unit traversed by a
β-particle is called the stopping power. This energy loss
rate due to inelastic processes, denoted by
ddEx
(in-
keV/µm), is estimated according to the Bethe-Bloch for-
mula for betas (Leo, 1987 [16]; Knoll, 1989 [14]) as fo l-
lows:
( )
2
2
02
d4π
d
E mc
SrNZ AB
x
β
= =+
(1)
with
H. H. MUHAMMED, Z. ZENGIN
Copyright © 2013 SciRes. ENG
14
and
( )
( )
3
2
2
1
11221ln 2
8
2
B
γγγ
γ


=+−+ −


where
0
r
is the classical electron radius (2.82 × 1013
cm),
( )
12
1
γβ
= −
,
vc
β
=
, N is the atomic-number
density for the medium, Z is the atomic number, and I is
the average excitation potantial of the medium in eV
(which is approximately
12 7IZ= +
for
12
eff
Z<
;e.g.
for water:
7.22
eff
Z
and
94 eVI=
according to Leo
(1987) [16]. The average energy loss S is valid for both
electrons and positrons. It can be noticed from the for-
mulae that A and B change slowly while S increases ra-
pidly as the particle slows down.
2.3. Hard Collisions and Knock-On (δ) Electrons
Production
While a
β
-particle with kinetic energy
T
E
passes
through a medium of atomic-number density
e
N
, hard
electron-positron collisions occur leading to a δ-electron
ray of energy
E
δ
. The number of these δ-electrons
created per unit distance is [24]:
( )
( )
2
,T
T
WE
N EEE
δδ
δ
=
(2)
where
( )
T
WE
is the incoming particle energy and is
given by:
( )
2
0
2
2
ee
T
rmN
WE
π
β
=
(3)
where h (Planck’s constant) = c (speed of light) = 1. Th is
formulais acceptable for
T
EE
δ
. δ-electronsare pro-
duced at energies higher than the cut-off energy level of
50 keV.
2.4. Multiple Coulomb Elastic Scattering from
Nuclei and Moliere’s Theory
With respect to the multiple Moliere’s theory, the proba-
bility of having an electron of momentum p and velocity
v, scattered at angle
θ
and the angular interval
d
θ
after traversing a thickness t within a material of atomic
number Z and atomic-number density N, can be written
as:
( )
( )
22
0
0
d
11
d dexpln
44
f
y yJyyby
θθθ
λλ λ


= −+




(4)
where y is a dummy variable,
0
J
is the zeroth-order
Bessel function,
c
λ θχ
=
, and b is defined by:
2
2
1167
bc
a
e
χ
χ
=
(5)
where
c
χ
characterizes the minumum single scattering
angle which is also called the angle parameter and is
given by:
( )
4
2
2
4( 1)
c
Nte ZZ
pv
π
χ
+
=
(6)
The screening angle
a
χ
describes the scattering by a
single parameter:
2 2
1.167
aa
χχ
=
(7)
and
( )
2
22 2
01.13 3.76
a
χχ α
= +
(8)
where
0
χ
is the critical scattering angle, below which,
deviations from the Rutherford scattering law occur.
0
χ
is given by :
01
3
0
0.885aZ
λ
χ
=
(9)
where
λ
is the electron DeBroglie wavelength, 0
a is
the Bohr radius,
hp
λ
=
is Planck’s constant and
2
Ze hv
α
=
is the fine structure constant.
For the derivation of
( )
,f
θ
Bethe (1953) [2] as-
sumed that
0c
χχ
which is acceptable for relatively
large thicknesstof a material, but it fails for
2
0
~~
b
c
ye
χχ
. The number of collisions 0
Ω which
occur within t is around
20
~
bc
e
χχ
. Moliere’s theory
is only valid for
0
Ω20>
where the constant value B
(defined as:
bB lnB= −
) > 4.5. For instance , in a water-
equivalent tissue at 25˚C, a β-particle of 80 keV energy
can traverse through approximately 0 .14 mm within the
tissue where nearly 30 collisions occur and B 5; H en ce ,
Moliere’s theory can be used here, but not in the case of
a 50 keV β-particle with a traversing distance of 0.06 mm
in water, where around 20 collisions might occur (Levin
and Hoffman, 1999 [18]).
For all angles,
( )
f
θ
can be expanded as a power se-
ries in
1
B
by variable-change to
( )
12
c
vB
θχ
=
to
obtain:
( )
( )
( )
( )
( )
( )
( )
012
12
d
d
f
vvfvBf vBfv
θθθ
−−

=+ ++…

(10)
in which
( )
( )
( )
2 22
0
0
1 dexpln
!4 44
n
n
fv
u uu
u uJvu
n

 
= −−

 

 

(11)
where
12
.u By=
If the scattering angle is much bigger than
0
χ
, the
distribution function
( )
f
θ
turns into the Rutherford
single scattering law; i.e.,
( )()
3
d 2d
R
fB vv
θθ θ
=
.
H. H. MUHAMMED, Z. ZENGIN
Copyright © 2013 SciRes. ENG
15
Hence, the ratio of the total prob ability to the Rutherford
scattering probability becomes:
( )()
( )
12
41
1.
2
R
Rf fvfBf
= =++…
(12)
where
(0)
f
,
(1)
f
and
(2)
f
are given as follows:
( )
( )
( )
( )
( )
45
2
01
4
21 5
2,
v
v
f vefvv
= =
, and
( )
( )()
( )
2
6 24
16 ln0.4
1 924
lnv
fvvvv
−−
+
=−−
(13)
In practice, it is enough to include the terms of
( )
1
f
and
( )
2
f
in the calculations to get an accurate result.
2.5. Initial-Energy Spectrum of Emitted Positron
This energy spectrum is required to be able to determine
the β-particles trajectories for a certain isotope. Fermi
revealed the expressions for the β-energy spectra as well
as for the half-life of β-decay (Fermi, 1934 [12]; Evans,
1955 [11]; Konopinski 196 6 [15]; Wu and Moskowski,
1966 [28]; Cengiz and Almaz, 2004 [6]). For instance,
the end point or maximum energy for
18
F
is 635 KeV
while its mean energy is 250 K eV and its half-life is 110
min. All isotopes are allowed or super-allowed β-decays
and their energy spectra distributions are given by:
( )()
2
max
, ()NEdEgFZEp EEEdE= −
(14)
where
( )
,FZE
is the Fermi function of the energy E of
the emitted positron. Z is the atomic number of the beta
decay daughter. N(E) is the number of decays at an
energy between E and E + dE, g is a coupling constant, E
and Emax are the total and the maximum (end point) β-
energies, respectively (in units of mc2), p is the mo-
mentum of β-particle (in units of mc) (Wu and Mos-
kowski, 1966 [28]; Daniel, 1968 [9]).
A non-relativistic approximation of
( )
,FZE
for al-
lowed β-decays of lighter elements is given by:
( )
2
,2 (1)
allowed
F ZEe
πη
πη
= −
(15)
where
Z Ep
ηα
= −
and 1 137
α
= (the fine structure
constant) (Wu and Moskowski, 1966 [28]; Daniel, 1968
[9]).
3. Monte Carlo Simulation
The Monte Carlo simulation begins with simulating the
medium where the positrons are emitted with energies
randomly chosen according to the distribution described
previously for the selected isotope.
To achieve a sufficient simulation accuracy of the po-
sitron annihilation process, the following three types of
processes (which affect the energies and directions of the
emitted positrons) should be taken into account: 1)
Energy loss by inelastic collisions (i.e. by ionization and
excitation processes), as des cribed in sub section (II.B). 2)
Hard collisions and δ-electrons production, as described
in subsection (II.C). 3) Elastic scattering by nuclei, as
discussed in subsection (II.D).
Figure 1 shows a 2D-projectionof the simulation re-
sults of 100 tracksfor 18F positronsgenerated from a point
source in water. These results correspond to those ob-
tained by Levin and Hoffman (1999) [18].
The simulation algorithm was implemented in Mat
lab®. The needed input parameters were selecting the iso-
tope and the medium in addition to defining the 3D-dis-
tribution of this isotope. The algorithm proposed and im-
plemented in this work performs the computations of the
positron energy loos in small steps (taking into account
the processes discussed previously). Furthermore, our
algorithm has the new feature or functionality that it ge-
nerates, at the beginning of each step, a little region or
neighborhoo d around the current pos ition of the positron.
This feature leads to the advantage that it becomes easy
to simulate what happens within that region with high
spatial resolution.
4. Conclusions
A user friendly implementation of the simulation algo-
rithm for the positron annihilation process was carried
out. The obtained results were satisfactory and of suffi-
cient accuracy. From Figure 1, it becomes very clear
how the positron annihilation process degrades the spa-
tial resolution of the PET imaging system, because each
point source or voxel of the imaged target object is con-
verted into a much larger 3D-spot or clouds. This leads to
a considerably blurred image, because the voxels (each
of which have a certain isotope concentration) of the im-
aged target object are positioned beside each other, while
Figure 1. A 2D-projection of the simulation results of 100 18F
positron tracks emitted from a point source centered in a
little water volume.
H. H. MUHAMMED, Z. ZENGIN
Copyright © 2013 SciRes. ENG
16
what the PET camera sensors detect are 3D-spots or
clouds (each of which is centered at one of these voxels)
that overlap each other exte nsively .
Hence, the next future step for this research work is to
use this efficient simulation tool to study the positron
annihilation process in order to understand and find out
how to reduce the positron range to be able to enhance
the spatial resolution of the PET imaging system.
It is of course also possible to study the other pro-
cesses and parameters that affect the spatial resolution of
the PET imaging system, such as the effect of the area of
the sensing head of the detector, and the effect of the
non-collinearity of the annihilation photons pair, which
means that the angle between these photons is not always
180˚.
REFERENCES
[1] J. Barba, J. Sempau, J. M. Ferntidez-Varea and F. Salvat,
PENELOPE: An Algorithm for Monte Carlo Simulation
of the Penetration and Energy Loss of Electrons and Po-
sitrons in Matter,” Nuclear Instruments and Methods in
Physics Resear ch B, Vol. 10, 1995, pp. 100:31-46.
http://dx.doi.org/10.1016/0168-583X(95)00349-5
[2] H. Bethe, “Moliere’s Theory of Multiple Scattering,”
Physical Review, Vol. 89, 1953, pp. 1256-1266.
http://dx.doi.org/10.1103/PhysRev.89.1256
[3] I. Buvat and I. Castiglioni, “Monte Carlo Simulations in
SPET and PET,” Quarterly Journal of Nuclear Medicine,
Vol. 46, No. 1, 2002, pp. 48-61.
[4] I. Buvat and D. Lazaro, “Monte Carlo Simulations in
Emission Tomography and GATE: An Overview,” Nuc-
lear Instrument s and Methods in Physics Research A, Vol.
569, 2006, pp. 323-329.
http://dx.doi.org/10.1016/j.nima.2006.08.039
[5] J. Cal-Gonzalez, J. L. Herraiz, S. Espana, M. Desco, J. J.
Vaquero and J. M. Udias, “Positron Range Effects in
High Resolution 3D PET Imaging,” Nuclear Science Sym-
posium Conference Record, 2009.
[6] A. Cengiz and E. A lmaz, “Internal Bremsstrahlung Spec-
tra of β-Particle Emi tters Using the Monte Carlo Met hod,”
Radiation Physics and Chemistry, Vol. 70, 2004, pp.
661-668. http://dx.doi. org/10.1016/j.radphyschem.2004.0
3.008
[7] J. C. L. Chow, M. K. K. Leung and D. A. Jaffray, “Monte
Carlo Simulation on a Gold Nanoparticle Irradiated by
Electron Beams,” Physics in Medicine and Biology, Vol.
57, 2012, pp. 3323-3331.
http://dx.doi.org/10.1088/0031-9155/57/11/3323
[8] K. Contractor, A. Challapalli, G. Tomasi, L. Rosso, H.
Wasan, J. Stebbing and Others, “Imaging of Cellular Pro-
liferation in Liver Metastasis by [18F]Fluorothymidine
Positron Emission Tomography: Effect of Therapy,” Phy-
sics in Medicine and Biology, Vol. 57, 2012, pp. 3419-
3433.
[9] H. Daniel, “Shapes of Beta -Ray Spectra,” Review s of Mo-
dern Physics, Vol. 40, 1968, pp. 659-672.
http://dx.doi.org/10.1103/RevModPhys.40.659
[10] S. Espana, J. L. Herraiz, E. Vicente, J. J. Vaquero, M.
Desco and J. M. Udias, “PeneloPET, a Monte Carlo PET
Simulation tool Based on PENELOPE: Features and Va-
lidation,” Physics in Medic ine and Biology, Vol. 54, 2009,
pp. 1723-1742.
http://dx.doi.org/10.1088/0031-9155/54/6/021
[11] R. D. Evan, “The Atomic Nucleus,” McGraw-Hill, New
York, 1955, Chaps 18-22.
[12] E. Fer mi, “Towards the Theory of β-Rays,” Zeitschrift für
Physik, Vol. 88, 1934, p. 161.
http://dx.doi.org/10.1007/BF01351864
[13] L. Jodal, C. Le Loi rec and C. Champion, “Positron Range
in PET Imaging: An Alternative Approach for Assessing
and Correcting the Blurring,” Physics in Medicin e and Bi-
ology, Vol. 57, 2012, pp. 3931-3943.
http://dx.doi.org/10.1088/0031-9155/57/12/3931
[14] G. F. Knoll, “Radiation Detection and Measurement,”
Wiley, Ne w York, 1989, Chap 2.
[15] E. J. Konopinski, “The Theory of Beta Radioactivity,”
Clarendon Press, Oxford, 1966.
[16] W. R. Leo, “Techniques for Nuclear and Particle Physics
Experiments,” Springer, New York, 1987, Chap 2.
http://dx.doi.org/10.1007/978-3-642-96997-3
[17] C. S. Levin, M. Dahlbom and E. J. Hoffman, “A Monte
Carlo Correction for the Effect of Compton Scattering in
3-D PET Brain Imaging,” IEEE Transactions on Nuclear
Science, Vol. 42, 1995, pp. 1181-1185.
http://dx.doi.org/10.1109/23.467880
[18] C. S. Levin and E. J. Hoffman, “Calculation of Positron
Range and Its Effect on the Fundamental Limit of Posi-
tron Emission Tomography System Spatial Resolution,”
Physics in Medicine and Biology, Vol. 44, 1999, pp. 781-
799. http://dx.doi.org/10.1088/0031-9155/44/3/019
[19] W. E. Meyerhof, “Elements of Nuclear Physics,” Mc-
Graw-Hill, New York, 1967.
[20] A. Murray and D. L. Williams, “Organic Synthesis with
Isotopes,” Interscience, New York, 1958.
[21] D. G. Ott , “Synthesis with Stable Isotopes,” Wiley-Inters-
cience, New York, 1981.
[22] H. Paganetti, “Range Uncertainties in Proton Therapy and
the Role of Monte Carlo Simulations,” Physics in Medi-
cine and Biology, Vol. 57, 2012, pp. R99-R117.
http://dx.doi.org/10.1088/0031-9155/57/11/R99
[23] H. Peng and C. S. Levin, “Study of PET Intrinsic Spatial
Resolution and Contrast Recovery Improvement for PET/
MRI Systems,” Physics in Medic ine and Biology, Vol. 57,
2012, pp. N101-N115.
[24] D. M. Ritson, “Techniques of High Energy Physics,”
Interscience, New York, 1961, Chap 1.
[25] A. F. G. Rocha and J. C. Harbert, “Textbook of Nuclear
Medicine: Basic Science,” Lea and Fekigur, Philadelphia,
1978.
[26] L. G. Strauss and P. S. Conti, “The Applications of PET
in Clinical Oncology,” Journal of Nuclear Medicine:
Official Publication, Society of Nuclear Medicine, Vol.
32, No. 4, 1991, p. 623.
H. H. MUHAMMED, Z. ZENGIN
Copyright © 2013 SciRes. ENG
17
[27] R. M. Thomson and I. Kawrakow, “On the Monte Carlo
Simulation of Electron Transport in the Sub-1 keV Ener-
gy Range,” Medical Physics, Vol. 38, 2011, pp. 4531-
4534. http://dx.doi.org/10.1118/1.3608904
[28] C. S. Wu and S. A. Moskowski, “Beta Decay,” Inters-
cience, New York, 1966.
[29] S. Ya mamoto, M. Imaizumi, T. Watabe, H. Watabe, Y.
Kanai, S E. himosegawa and J. Hatazawa, “Development
of a Si-PM-Based High-Resolution PET System for Small
Animals,” Physics in Medicine and Biology, Vol. 55,
2010, pp. 5817-5831.
http://dx.doi.org/10.1088/0031-9155/55/19/013