Monte Carlo study of a 3D Compton imaging device with GEANT4
In this paper we investigate, with a detailed Monte-Carlo simulation based on Geant4, the novel approach  to 3D imaging with photon scattering. A monochromatic and well collimated gamma beam is used to illuminate the object to be imaged and the photons Compton scattered are detected by means of a surrounding germanium strip detector. The impact position and the energy of the photons are measured with high precision and the scattering position along the beam axis is calculated. We study as an application of this technique the case of brain imaging but the results can be applied as well to situations where a lighter object, with localized variations of density, is embedded in a denser container. We report here the attainable sensitivity in the detection of density variations as a function of the beam energy, the depth inside the object and size and density of the inclusions. Using a 600 keV gamma beam, for an inclusion with a density increase of 30% with respect to the sorrounding tissue and thickness along the beam of 5 mm, we obtain at midbrain position a resolution of about 2 mm and a contrast of 12%. In addition the simulation indicates that for the same gamma beam energy a complete brain scan would result in an effective dose of about 1 mSv.
keywords:Gamma-ray imaging, Compton imaging, Germanium strip detector, Nondestructive testing
Pacs:29.40.Gx, 29.40.Wk, 07.85.Fv, 81.70.-q
For long time Compton scattering has been widely employed
as an imaging tool. A photon in the energy range 0.12.0 MeV
during its passage through matter, will dominantly Compton scatter against
the electrons of the material with an intensity proportional to electron
(and mass) density. If the energy of the incident and scattered photons are
known, it is possible to locate the scatter volume using the well known
energy–angle relationship. Photon scattering can therefore readily provide
3D density maps of extended objects of unknown composition.
This is a clear advantage of scattering compared to other photon
imaging techniques, like the well established x–ray applications, which
make use of photons transmitted through the material to be imaged.
In this case the integrated density of the object along the path traversed
by the photon beam is provided. Several profiles from different directions
need to be collected and complex algorithms are required in order
to reconstruct the 3D density distribution.
In the last fifty years Compton scattering has found a broad range of applications in medicine, as imaging tool – or to measure the human tissue density . Compton scattering is also widely used in the field of non destructive testing (NDT) for quality control to detect defects like cracks, inclusions and voids inside industrial manufactures –. The possibility of using backscattered radiation has made Compton scattering very attractive also in those situations where the access to the object to be examinated is restricted to one side, like inspection of airframes  and soil density determination .
More recently attention has been focused on reconstructing the photon arrival direction (for instance from implanted nuclides) exploiting Compton scattering in a position sensitive detector –: this method allows a conical region of possible photon origin to be defined, and with several measurements the true origin can be inferred by intersection of different cones.
In spite of its attractive features, photon scattering has not become the dominant approach in the imaging field . There are two major limitations of this technique: the background induced by events multiply scattered inside the object and the attenuation of both incident and scattered photons. Many attempts to overcome these limitations have been pursued. We recall the use of the information of the transmitted beam to determine the attenuation coefficient, or the use of collimators to define a small region of interest in the examinated object in order to reduce the multiple scatter component. Also the Monte Carlo approach has proven to be able to provide corrections for this source of background [19, 7].
In this paper we want to further explore the potential of the novel method to image the electron density in the human body proposed in  by means of a detailed Monte Carlo (MC) simulation using GEANT4  with a realistic set–up. The application presented here refers to brain imaging but the results can be applied as well to situations where a lighter object, with localized variations of density, is embedded in a denser container.
This paper is organized as follows: in Section 2 we will briefly recall the main characteristics of Compton scattering, in Section 3 the experimental technique will be described and in Section 4 details about the simulation will be given. In Section 5 we will describe the results found and the expected performance of the proposed device, in Section 6 we will give an evaluation of the expected dose, and finally in Section 7 we will give our conclusions.
2 Compton Scattering
In the Compton effect a photon of energy interacts with an atomic electron and it is scattered, with a lower energy , at an angle with respect to the incoming direction. In case of scattering by free electrons at rest the conservation of energy and momentum implies the well–known Compton formula which relates initial and final photon energies to the scattering angle:
where keV is the electron mass (natural units with speed of light are used). The kinetic energy of the struck electron is and the maximum is:
The probability of interaction is given by the Klein-Nishina formula :
where m is the electron classical radius. Eq. 1 and 3 are however only approximations. In reality electrons are not free and transitions of bound electrons are allowed only if the energy transfer is larger than the ionization energy. Moreover electrons are not at rest but move with a certain momentum distribution. This produces the so called ”Doppler broadening” i.e. a smearing with a long, non gaussian tail, of the scattering angle. Both the binding effects and the Doppler broadening are however significant at photon energies less than 100 keV, well below energies used in the present study.
3 Experimental technique
The proposed technique  makes use of a well collimated and monochromatic gamma-emitting source which illuminates the part of the body/object to be imaged. An x-y-z orthogonal coordinate system is introduced such that the axis is along the beam direction. Photons Compton scattered inside the body/object are recorded by means of a surrounding germanium strip detector in the position and its energy is measured.
Using Eq. 1 the angle is calculated and the scattering position is reconstructed by the formula:
The distribution so obtained provides, directly and without any inversion algorithm, the density distribution of the material crossed by the beam i.e. a small cylinder with transverse dimensions given by the beam size and length given by the object length. A full imaging is obtained by raster scanning the surface of the object. Typical scattering applications make use of collimators in front of the detector to precisely define the scatter volume (voxel), given by the intersection of the incident and detected photon. The method proposed here is collimator free. This aspect provides larger counting rates and therefore lower dose absorbed by the patient.
The study was developed in the GEANT4  framework.
The set–up used for the simulation (see fig. 1) consists
of a monochromatic gamma source placed at 555 mm from the head
center, which is positioned in , and at 5 mm distance from a collimator.
This is a lead tube 400 mm long and of inner and outer radius 1 mm
and 10 mm, respectively.
The obtained gamma beam is directed onto the head volume.
The modelling of the human head is based on the ”human-phantom”
example provided with the GEANT4 package, the anthropomorphic MIRD
 phantom was used.
The shape of the head is a cut ellipsoid and consists of three
Externally we find a skin layer made by soft tissue, then the skull
volume made by bone tissue with density =1.486 g/cm.
The skull contains the brain volume made by soft tissue with density
=1.040 g/cm (the same density was used for the skin).
To test the capability of the apparatus to detect inclusions of
different density, we positioned inside the brain volume and on the
beam axis, four small target volumes of 2x2 mm in the transverse
beam direction and 5 mm in the longitudinal one.
These targets are made by soft tissue but a greater value for the
density of =1.35 g/cm was used. This value was chosen in
accordance with magnetic resonance measurements of the density of
brain pathological tissues .
The targets (numbered 0–3) were placed at the positions
-45; -15, 0 ,15 mm along the beam line.
Photon detection is achieved by means of a segmented HPGe detector.
Two half–rings of such devices surround the head volume.
The sensors are parallelepiped 5 mm thick, 100 mm wide and 100 mm high
placed according to an hexadecagon shape to cover the largest possible
azimuthal region. A 5 mm thick lead plate placed in the rear part of
the head acts as a beam stopper. Each simulation consisted of
events. Ten different gamma energies were simulated from 100 to 1000
keV in steps of 100 keV.
Interactions of photons and electrons were simulated using the GEANT4
implementations of physics models developed for the PENELOPE code
. This code has been specifically developed for and
transport in matter and great care was given to the description
of the low energy processes including binding effects and Doppler broadening.
Comparisons with exeprimental data showed that PENELOPE provides consistent
results in the energy range from few keV up to about 1 MeV .
Photons were tracked through the described set-up. Step by step the type of interaction, its position and the volumes crossed were recorded together with the energy deposited in the head volumes for dose calculations. Hits produced in the HPGe detectors were also recorded. The energy resolution of the sensors was taken to be gaussian and was simulated taking the standard 2.96 eV energy to produce an electron-hole pair, a Fano factor of 0.05 and a constant electronics noise of 500 eV (see for example [17, 26], for an analysis on energy and angular resolution for CdZnTe detectors see  ). Hits within a distance of less than 3 mm were grouped together to reconstruct a cluster. The sum of hits energies, after the smearing to account for the detector resolution, was taken to be the cluster energy. Its position is given by the energy weighted average of the hits position. Fig. 2 shows the distribution of the number of hits per event in the HPGe detectors (a), their energy (b), their distance (c) and the number of reconstructed clusters for each scattered photon (d) for a beam energy of 600 keV.
The best resolution for the proposed setup is achieved at . Applying error propagation on Eq. 4 we obtain:
where means sum in quadrature. In the above equation
is the resolution on the photon impact
position along the axis given by the strip granularity.
distance of the photon impact point on the detector from the beam axis,
for the events undergoing a single Compton scattering (see
Sec. 5.1) the average value in the simulated beam energy
range (1001000 keV) varies between 141 and 148 mm;
is the Doppler broadening
contribution, it goes from to about depending on the
Both quantities show a steep energy dependence reaching a plateau at about 500 keV.
is the energy resolution of the germanium detector,
the term dependent on the scattering angle is 1 at
and less than 2 in the range .
and depend on
the germanium detector thickness and strip pitch, but their contribution
is suppressed by a factor which is an average
over the accepted range. In the simulated energy range this quantity
goes down from 0.39 to 0.36 reaching its final value at about 500 keV.
Finally the distances in and between the beam and the detector position
divided by are less than 1 mm, both having an average value of about 0.6 mm.
The error induced by the beam divergence and width turns out to be:
where the first term is the dominant one.
Both these sources of uncertainty are determined by the collimator geometry:
mrad and mm.
This formula is valid for a non attenuated beam, beam attenuation effects
are properly taken into account by the simulation.
The resolution improves as decreases, i.e. the detector must be as close as possible to the object to be imaged. The resolution also improves if the gamma energy increases, both because the energy resolution improves and because the Doppler broadening contribution decreases. As it will be shown later, the number of scattered and detected photons decreases as increases and also the background changes. The best gamma energy is thus a compromise among different aspects.
5 Results of the simulation
5.1 Event types and energy distributions
For the analysis we select events which have only one reconstructed
cluster in the HPGe detector and an energy compatible with that of
Compton scattering. We classify events in four types
according to the following scheme:
S1: the primary underwent a single Compton scattering inside the head volume and then deposited all its energy producing one cluster inside the HPGe.
Sn: in this case the has undergone multiple Compton interactions inside the head volume before releasing all its energy in the HPGe.
Escapers: the primary photon after one or more interactions in the head volume gives a signal in the germanium detectors but not all of its energy is released. The or one of its daughter particles escape without being detected.
NoGamma: for these events the signal in the HPGe was not directly produced by the primary photon.
From the above description it is clear that only events of type S1 are those which carry the correct relation between energy and angle
and can be used to determine the position of the scattering, the
other events constitute the background to our signal.
Figure 3 shows the energy distribution recorded by the
HPGe detectors for three different beam energies. The various histograms
correspond to the event type described above and the lines indicates
the energy range used for the event selection. The lower limit represents
the minimal energy required for a Compton scattering (i.e.
) and the upper one is the energy where the S1
component becomes less than the Sn one.
In the plot relative to the 300 keV beam energy is visible, at the same
energy of the beam, the peak of the coherently scattered events
The peak is not visible at higher energies since the cross-section for this process decrease rapidly () at growing energies. The peak visible at the escape line of the lead at 75 keV and 85 keV in the distributions of NoGamma events indicate that these events are mostly due to electrons emitted by photoelectric effect in the beam stopper and backscattered onto the HPGe.
The percentage of the four types of events as a function of beam energy is given in Table 1.
|E(keV)||S1 (%)||Sn (%)||Escaper (%)||NoGamma (%)|
We see that for beam energies up to 500 keV the event composition is rather stable with about 70% of events being of type S1. For higher energies the component due to Escapers rises rapidly and at 1 MeV the background amounts to 50%.
5.2 z distributions
Figure 4 show the distribution of the reconstructed position
for three different beam energies, from top to bottom 300, 600 and 900 keV
respectively. The various histograms refer to the different event types (see
figure caption for details).
In this picture the head volume extends from -75 mm to 75 mm, the beam
is directed from to .
Distributions are shown in the 100 mm range but they extend in a wider range (-200500 mm); outside the interval 80 mm only background events contribute. At lower energies effects of attenuation are dominant and clearly determine the shape of the distribution. For all beam energies the distribution of Sn events is slowly increasing with increasing and shows a long plateau.
The Escaper component, on the contrary, increases rapidly with beam energy
producing a remarkable distortion on the distribution at large values.
Since the detected energy in this case is always smaller than the real one
the angle of the Escaper events is biased towards values greater
than , i.e. towards positive values of .
The presence of localized increase of density in the head volume induces two effects . An increase of the number of scatterings, producing a peak in the distribution, and an increase in the attenuation of primary photons impinging on the downstream material. This last effect is responsible for the change in the slope of the distribution. Inclusions with lower density or voids, will manifest themselves as dips. Four peaks due to the targets and two at the edges of the head volume, due to the bone tissue, are visible in Fig. 4. They tend to disappear at large for lower beam energy. Higher energies on the contrary show a better penetration capability. We notice how the use of scattered photons makes possible the detection of changes in the composition of the object to be imaged independently of beam attenuation. With an imaging technique based on photon transmission, in the case of a low density material enclosed in a dense container, the dominant contribution to the signal would be provided by the container itself.
5.3 Raw signal and resolution
The raw signal, i.e. not corrected for the beam attenuation, produced by the four targets in absence of resolution effects would be the target profile along , a rectangle 5 mm wide with sharp edges. The resolution modifies this profile blurring the borders. To extract the raw signal from the data we used a model of a rectangular signal pulse convoluted with a gaussian (the blurring function) over an exponential background (this corresponds to the signal of the surrounding tissue and that induced by the Sn and Escaper events). The fit provides the mean position and the height of the signal together with the of the gaussian modelling resolution effects. We studied the signal as a function of beam energy in the range 3001000 keV, simulating events for every energy value. For all investigated energies the results of the fit show that the positions of the targets are correctly reproduced within 0.1 mm. Figure 5 shows, as a function of beam energy and for the four targets, the number of events obtained integrating the signal peak.
The signal is visible for all the four targets over the full energy range. It decreases with increasing beam energy and target position. For the two downstream targets (2 and 3) and for energies above 600 keV the decrease is less pronounced due to the loss of resolution and background increase, which make the peak larger. The resolution as a function of beam energy is reported in Table 2.
As a general trend, the slightly decreases at increasing energies and the two deepest targets usually have larger values. The smallness of the signal for higher energies and deeper targets together with the non perfect modeling of the background and high correlation between fit parameters are responsible for departures from the general trend. We also studied the variation of the signal as a function of the target thickness. We used a beam energy of 600 keV to ensure adeguate penetration. Results for target thickness in the range 1–5 mm are shown in Fig. 6. As expected, the signal increases as a function of the target thickness. Already for 1 mm thickness a signal is visible for all the four targets although at a reduced level. Even if with this small change of density it would not be possible to provide a correct image of the object, the presence of the signal would already fulfill the requirements needed for most NDT applications. As an exercise we tried to obtain the target thickness using the obtained resolution as input to the fit function. As it can be seen in Table. 3, the capability to determine the target size improves with increasing target thickness and deteriorates with increasing depth in the object.
We also investigated the signal dependence on the target density, Fig. 7 shows the results. The signal for densities less than 1 g/cm appears as a reversed peak with respect to the exponentiall falling signal of the surrounding material. The results of the simulations indicate that, at least for the first two targets, inhomogeneities of about 6% could be located. Finally, in Fig. 8 we give the number of raw signal events as a function of the number of generated events in the simulation.
5.4 Corrections for background and attenuation
So far we have described the capability of the proposed device to detect localized changes in the density of an extended object. The density variations produce bumps (or dips) in the distribution of the reconstructed z position. To produce an image of the interior of the object, however, corrections for the background induced by multiply scattered events and Escapers, as well as for beam attenuation need to be implemented.
First we consider the background. The component induced by multiply
scattered events affects the whole z range in an almost flat way
(see Fig. 9-a) while the Escaper component grows with beam
energy and affects the forward z region. The background is inherently
connected to geometrical aspects (and to the detection capabilities of
the HPGe for the Escaper component) and does not depend on the density
distribution of the object to be imaged.
We tested this feature by comparing the z distributions of multiply
scattered events and Escapers of the various simulations done with different
size and density of the targets and with targets removed.
We also added, inside the brain volume, other volumes of different
sizes and densities in order to simulate non uniformities in the
Finally we replaced all the head materials with water.
All these density variations modify the scale of the background
distribution but not the shape.
This indicates that the MC approach can be effective in removing the unwanted
background component also in a ”real world” application of this device.
We notice also that, alternatively, the background component induced
by the Escaper events could be eliminated by surrounding the HPGe
detectors with some highly sensitive scintillator, like the NaI type,
operated as a veto.
In our analysis we used the tail of the distribution in the range 150300 mm to normalize the distribution of ”DATA”, i.e. our reference sample, to ”MC”, the simulation without targets. We chose this region because it receives contributions only from background (see Fig. 9-a). With proper normalization we can subtract the background component and compare the obtained distribution to the original S1 event distribution in Fig. 9-b.
To provide an indication of how a correction for the beam attenuation could be implemented (a full treatment of this subject is beyond the scope of this work) we consider the passage of our pencil beam through the object: on its path it will impinge on a small cylinder of material with length equal to length of the object. We subdivide this cylinder in n equal voxels of thickness dz; if is the number of photons entering the -th voxel the number of photons reaching the next voxel along z will be:
where is the attenuation coefficient which, given the photon energy, is dominated by Compton scattering. In the above relation it is also assumed that the voxel thickness is small enough to allow the replacement of the exponential by its linear expansion. The signal induced by singly scattered events in the th voxel will be:
where is the detection efficiency. It accounts for the probability that the scattered photon will not undergo to any other interaction on its way out of the head volume and that it will produce a signal in the HPGe detector. We assume that this efficency is the same for all the voxels, i.e. that does not depend on . To correct the -th signal for beam attenuation the number of scattered in all the preceding voxels must be iteratively added. With some algebra we obtain:
where the index runs from the first voxel to the th one and is
the number of event impinging on the first voxel. Given the simmetry
of the set-up the first and the last voxel should give the same signal since they
are made by the same material. This constraint can be used to
determine the value in the above formula.
The distribution, corrected for background and for attenuation is shown in Fig. 9-c. In spite of the simplicity of the outlined approach we see a reasonable flatness of the regions corresponding to the homogeneous brain tissue in which targets are embedded. At large however an overcorrection appears. This is probably due to a residual dependence of the detection efficiency related to the geometrical acceptance of the detector. Also the ratio of the corrected (for attenuation and resolution effects) signal peaks to baseline, which should equal the ratio of the target to brain density, is reproduced in an acceptable way. A value of about 1.15 was obtained for the first three targets, to be compared to the true one of 1.30. We also provide an estimate of the image contrast which could be obtained by this device. We define the contrast C as:
where is the mean of the four maxima of the signal from the inserts and is the average of the regions between the inclusions, which correspond to the uniform soft tissue. A value of 12% is obtained for a beam energy of 600 keV. A deconvolution was performed to cross-check the results on the expected resolution obtained, see Table 2, with data not corrected for attenuation and background. The fit function was a rectangular pulse convoluted with a gaussian, the true position and thickness of the four inserts were used as input to the fit, the background was modeled with a linear polynomial. The resolution for the four inserts, extracted from the corrected data, agree with those reported in Table 2. Finally we remark that, contrary to other imaging approaches involving Compton scattering, this method provides directly and without the need of any backprojection technique, the density distribution of the object. What is shown in Fig. 9-c is the image corresponding to a portion of the object 2 mm wide in and directions (the beam size) along its coordinate for a given beam position. The full 3D image is obtained by raster scanning the surface of the object and merging together all the acquired images.
6 Dose calculations
To calculate the absorbed dose by the soft (i.e. skin and brain) and skeleton tissue the head volume was divided in cubic voxels of 1 mm side. During the simulation all the particles, primary or secondaries, were followed and step by step the energy deposited in each voxel, provided by the GetTotalEnergyDeposit() method of GEANT4, was stored. The energy released in soft and skeleton tissue, for our standard simulation of events is reported as a function of the beam energy in Fig. 10. The energy deposited increases linearly with beam energy.
Given the highly collimated beam which is used in our application the most
of the energy deposition is concentrated in a narrow cylinder around the beam path
as can be seen in Fig. 11 and Fig. 12, which
show, for a beam energy of 600 keV, a longitudinal projection of the deposited
energy in soft and skeleton tissue, respectively.
Although almost all the voxels in which the head volume is partitioned
are interested by the release of energy, the energy deposited in 1 mm
radius around the beam path amounts to about one half of the total.
This highly non homogeneous energy deposition makes unpracticable the usual
way of calculating the absorbed dose i.e. dividing the total energy released
in the voxels by their total mass. Furthermore the region to be imaged needs
to be illuminated by raster moving the beam on its surface.
During a scan therefore every part of the imaged region will receive an energy
deposit like that of Figs. 11 and 12
(which refer to one scan position), when illuminated by the beam, and an
additional contribution, which depends on the distance from the beam, when
To estimate the absorbed dose for a complete scan of the brain volume we firstly obtained an analytical representation of the dose as a function of the distance from the beam in the transverse plane. We considered voxels of 22 mm (i.e. the beam spot size) in the transverse direction and 130/17 mm deep along the beam for soft/skeleton tissue, respectively. Then we simulated the scan procedure by considering an array of such voxels 80 mm wide and 40 mm high and moving the beam, in steps of 2 mm, upon it. For every beam position the dose delivered to the illuminated voxel and to all the other was recorded. This study indicates that the dose is almost uniformly distributed over the array and reaches a maximum value of 110 mGy, the contribution of soft and skeleton tissue being equally distribued. This value drops of about 15% for voxels on the array border. Using the radiation and tissue weighting factors recommended in  we obtain for a complete scan with photons for each point, at an energy of 600 keV, an effective dose of about 1 mSv which is half of the average dose imparted during a computed tomographic examination of the head .
In this paper we have studied the potential of a novel approach of 3D
imaging with photon scattering.
We studied as an application of this technique the case of brain imaging and show how the proposed device could be able to detect localized density variations inside it. Results on the attainable sensitivity were shown as a function of the beam energy, of the number of impinging photons and of various target parameters: depth, thickness and density.
We studied in detail the signal and background characteristics and outlined a method, relying on the MC approach, to remove the background. A simple model to correct the signal for beam attenuation has also been introduced.
The simulation indicates that a complete scan of the head volume at 600 keV photon energy would result in an effective dose of about 1 mSv which is comparable to the average dose imparted during a computed tomographic examination of the head.
We would like to thank Cinzia Talamonti for very helpful and stimulating discussions.
-  M. Lenti, Nucl. Instrum. Methods A 588, (2008) 457.
-  P.G. Lale, Phys. Med. Biol. 4, (1959) 159.
-  F.T. Farmer, M.P. Collins, Phys. Med. Biol. 16, (1971) 577.
-  J.J. Batista et al., Phys. Med. Biol. 22, (1977) 229.
-  J.J. Batista, M.J. Bronskill, Phys. Med. Biol. 26, (1981) 81.
-  El Khettabi et al., Phys. Med. Biol. 48, (2003) 3445.
-  M.S. Mooney et al., Phys. Med. Biol. 41, (1996) 2399.
-  J.S. Al-Bahri, N.M. Spyrou, Appl. Radiat. Isot. 49, (1998) 1677.
-  G. Harding, J. Kosanetzky, Nucl. Instrum. Methods A 280, (1989) 517.
-  J.M. Sharaf, Appl. Radiat. Isot. 65, (2007) 1330.
-  B. Achmad, E.M.A. Hussein, Appl. Radiat. Isot. 60, (2004) 805.
-  A.C. Ho, E.M.A. Hussein, Appl. Radiat. Isot. 53, (2000) 541.
-  B.L. Evans et al., Nucl. Instrum. Methods A 480, (2002) 797.
-  F.A. Balogun, P.E. Cruvinel, Nucl. Instrum. Methods A 505, (2003) 502.
-  E.A. Wulf et al., IEEE Trans. on Nucl. Sci. 50, (2003) 1182.
-  C.J. Hall et al., Nucl. Instrum. Methods A 513, (2003) 47.
-  J. van der Marel, B. Cederwall, Nucl. Instrum. Methods A 471, (2001) 276.
-  E.M.A. Hussein Nucl. Instrum. Methods B 263, (2007) 27.
-  R.D. Speller, J.A. Horrocks, Med. Phys. 15, (1988) 707.
-  S. Agostinelli et al, GEANT4 Collaboration, Nucl. Instrum. Methods A 506, (2003) 250.
-  J.D. Bjorken and S.D. Drell, Relativistic Quantum Mechanics, McGraw-Hill, New York, 1964.
-  W.S. Snyder, et al, J. Nucl. Med. Suppl., 3, (1969) 5.
-  M. Just, M. Thelen, Radiology 169, (1988) 779.
-  PENELOPE 2008 – Workshop Proceedings, Barcelona, Spain 30 june–3 july 2008. NEA n. 6416.
-  J. Sempau et al., Nucl. Instrum. Methods B 207, (2003) 107.
-  G.F. Knoll, Radiation Detection and Measurement, Wiley and Sons, New York, 2nd edition, 1989.
-  Y.F. Du et al, Nucl. Instrum. Methods A 457, (2001) 203.
-  A.D. Wrixon, J. Radiol. Prot. 28, (2008) 161.
-  F.A. Mettler et al., Radiology 248, (2008) 254.