Monte Carlo simulation of the spatial response function of a SPECT measurement device for nuclear fuel bundles

Stockholm University

Hospital physics program

Bachelor’s Thesis

Monte Carlo simulation of the spatial response function of a SPECT

measurement device for nuclear fuel bundles

Author:

Emilie Dul

emilie.m.dul@gmail.com

Supervisors:

Peter Jansson Staffan Jacobsson Sv¨ ard

October 1, 2017

Abstract

The PGET device is currently being developed for partial-defect verification purposes on nuclear fuel

assemblies. It Comprises CdTe detector elements in a heavy tungsten-alloy collimator, for which

collimator slit openings define the field-of-view. This study aims at calculating the spatial response

function of this device for further deployment in tomographic reconstruction algorithms. In this work,

the detector response for 2 different sources (662 keV from Cesium-137 and 1274 keV from Europium-

154) was simulated using the MCNPX software package. In the simulations, energy windows used in

measurements with the PGET device were deployed. The results show the expected characteristics

with strong response for a source position directly in front of the collimator slit opening and decreasing

response as the source is moved into the penumbra and umbra region. The uncertainty of the simulated

response function was less than 3.5 % for both sources. Separate simulations were made to quantify

contributions from septal penetration and scattering from the collimator material into the detector

for the energy windows covering the full -energy peak. These contributions were found to be around

3% for the source of Cesium-137 and 6% for the source of Europium-154.

Abbreviations

IAEA : International Atomic Energy Agency NPT: Non-Proliferation Treaty

PGET: Passive Gamma-Emission Tomography (an existing device of the IAEA) UGET: Universal Gamma-Emission Tomography (a notional device)

CdTe: Cadmium Telluride. A detector material.

LSF: Line-Spread Functions.

CSDA: Continuous Slowing Down Approximation

Contents

1 Introduction 5

2 Background 8

2.1 Geometry of the detector and collimator slit in the PGET device . . . . 8

2.2 Line-Spread functions and energy dependencies . . . . 10

3 Materials and Methods 11 3.1 Simulation procedure . . . . 12

3.2 Modelling gamma-transport using MCNPX . . . . 13

3.2.1 Cell and surface definition . . . . 13

3.2.2 Material definition . . . . 14

3.2.3 Source particles, energy bins and cut-off . . . . 15

3.3 Execution of MCNPX simulations using Matlab . . . . 16

3.4 Extraction of MCNPX tally data using Matlab . . . . 16

4 Results 16 4.1 Gamma-ray energy spectra . . . . 16

4.2 LSFs for a source at x=160 mm . . . . 18

4.2.1 LSFs for the source of Cesium-137 . . . . 18

4.2.2 LSFs for the source of Europium-154 . . . . 19

4.2.3 Direct exposure vs septal penetration and scattering . . . . 19

4.3 Full response function . . . . 20

4.3.1 LSFs for the source of Cesium-137 . . . . 20

4.3.2 LSFs for the source of Europium-154 . . . . 22

5 Discussion 24 5.1 Gamma-ray energy spectrum . . . . 24

5.2 Lateral cut-off of the spatial response function . . . . 24

5.3 Contribution from septal penetration . . . . 24

6 Conclusion 25

7 Acknowledgement 25

8 Appendix 28

8.1 Appendix 1: MCNPX input file (emilie.i) . . . . 28

8.2 Appendix 2: Matlab execution file . . . . 30

8.3 Appendix 3: Matlab extraction file . . . . 31

1 Introduction

Nuclear weapons are the most powerful weapons ever invented. Although they have rarely been used in warfare, they have demonstrated unprecedented destructive powers. The proliferation of nuclear weapons is therefore a matter of intense concern in geopolitics.

The United Nations has worked for decades to eradicate those weapons by means of bilateral non- proliferation treaties. So far, few countries possess a nuclear arsenal. The proliferation of nuclear weapons is regulated by the International Atomic Energy Agency (IAEA). This organization has the responsibility of conducting safeguards to verify the compliance of signatories to the Non-proliferation treaty (NPT) [1].

Nuclear weapon manufacturing requires fissile material, which can be achieved using one of two main options: uranium-235 or plutonium-239. Uranium-235 is available in the earth’s crust as 0.7 % of the natural uranium. Its extraction and trade is closely supervised by the IAEA. Nations can acquire Uranium-235 for energy production in power plants. The amount of Uranium-235 acquired by the nation is recorded by the IAEA, and safeguards inspections verifying that Uranium-235 is indeed used as nuclear fuel are conducted. Pu-239 is created in nuclear fuel when uranium is irradiated in nuclear reactors. Accordingly, the non-proliferation of Pu-239 can also be controlled by verifying the presence and properties of the fuel.

One of the safeguards activities conducted by the IAEA is partial-defect verification of nuclear fuel

assemblies, which refers to the detection of possible diversion of a fraction of the fuel rods from a used

nuclear fuel assembly (a nuclear fuel assembly typically contains between 100 and 300 fuel rods) .

For this purpose, tomographic measurement devices can be used. The PGET device (Passive Gamma

Emission Tomography device), pictured in figure 1, has been developed for the purpose of partial-

defect verification of used nuclear fuel assemblies and the first prototype was successfully tested on

research reactor fuel in the European Joint Research Center in Ipra, Italy, in 2012 [2]. These tests

were followed by additional successful tests on commercial fuel in the Finnish reactors of Olkiluoto

(2013) and Loviisa (2014).

Figure 1: The PGET device (Passive Gamma Emission Tomography device), developed for partial- defect verification. It consists of two detector heads rotating around a cavity where nuclear fuel bundles can be placed

The processes applied in PGET are similar to those of Single Photon Emission Tomography (SPECT) as illustrated in figure 3. Two detector heads rotate around the measured nuclear fuel assemblies, obtaining projections of the gamma radiation field at many positions relative to the fuel rods. An image describing the radiation source distribution inside the fuel rods is then reconstructed using an adequate tomographic algorithm. Finally, partial-defect verification can be executed by analyzing the resulting image and counting the number of fuel rods in it.

Many fission products are found in a nuclear core following neutron irradiation. Their half-life vary from few seconds to several years as seen on the table below.

Figure 2: Characteristic fission products and associated gamma-ray emissions from spent fuel in 0.4

and 2.5 MeV energy region. [3]

Since the PGET device has tapered collimator slits, the amount of radiations reaching the detector is high. In order to avoid saturation issues, it is necessary to perform the measurement several years after the end of irradiation. Isotopes of interest for this study must therefore have a half-life longer than a few years, in order to still be present in the core when the study is performed. Another isotope selection criteria is the energy of the radiation emitted gamma radiation. Low-energy gamma may be absorbed by the surrounding fuel rods.

The data collection and analyses are typically based on one of the two long-lived gamma-emitting fission products Cs-137 (T

=30.1 years) or Eu-154 (T

=8.6 years). According to table 2, both emit gamma rays at relatively high energies; 661.65 keV and 1274.43 keV, respectively, which enables the escape and detection of gamma-rays also from the assemblies’ innermost sections, although atten- uation may limit the transmission, in particular for the lower Cesium-137 energy. These two isotopes dominate the gamma radiation from fuel assemblies with cooling times longer than 5 to 10 years.

Figure 3: Schematic view of the PGET device. The field of view of each detector element is defined by slit openings in a tungsten-alloy collimator

An important feature of the gamma tomography device is the collimator design. The characteristics of the collimator and its slit openings determine the region of the fuel that contributes to the acquired gamma-ray intensity in each detector-element position. A collimator with a small slit width will require more time to obtain a number of counts that is statistically confident, while it will also enable high spatial resolution of the resulting images. Larger slit openings will shorten the required measurement time, while worsening the resolution. Accordingly, the choice of collimator slit dimensions will be a trade-off between precision and time. In order to reduce the assay time to under 5 minutes, the PGET device uses tapered collimator slits (see section 2).

The choice of detector material is also crucial. The PGET device uses detector elements of Cadmium

Telluride (CdTe), which offers relatively poor full-energy peak efficiency [4]. This implies that this

device will have relatively poor capacity to quantify the nuclide composition inside a fuel rod. However,

the small size of the detector elements also means that many detectors can be stacked up in the detector

head, thus leading to shorter assay time. Accordingly, the PGET device offers good performance in

terms of collecting tomographic images relatively quickly, but worse performance in terms of analyzing

the rod inventory of the specific nuclides. An alternative device, UGET (Universal Gamma Emission Tomography device) has been suggested, for which simulation studies indicate better performance for quantitative nuclide measurements [3]. However, the PGET device currently offers the capabilities requested by the IAEA when it comes to time-efficient imaging, and it is foreseen to be the device used by the IAEA for partial-defect verification.

The purpose of this project is to simulate the particle transport in the PGET device in order to obtain the spatial response of this system for further deployment in so-called model-based, or algebraic, to- mographic reconstruction algorithms [5]. The spatial response has been calculated using the MCNPX software package, in terms of so-called line-spread functions (LSF). The details of the geometry of the PGET device are discussed in section 2.1. LSFs are discussed in section 2.2. The software tools and simulation procedures used to obtain the LSFs are presented in section 3, and results are given in section 4.

2 Background

2.1 Geometry of the detector and collimator slit in the PGET device

As described above, PGET is a device capable of addressing partial-defect detection of nuclear fuel assemblies. Taking advantage of the large number of high-energy gamma-radiations inside used nuclear fuel, it can image the fuel rods in the assembly structure in a non-destructive way. In order to do so, the collimated detector array rotates around the nuclear fuel assemblies at a selected axial level, measures the gamma radiation field and tomographic reconstructions are performed based on the recorded data.

Of importance for the quality of the resulting image is the spatial response of the data collecting system for each detector element and it also constitutes important input in algebraic reconstruction algorithms [5]. The spatial response will be governed by the shape of the collimator slit opening to which it is attached. In PGET, all slit openings are shaped the same way, as illustrated in figure 4 and 5. Also introduced in figure 4 and 5 is the system of coordinates used to define the spatial response function. In this experiment, the x-axis is along the collimator slit axis, the z-axis is parallel to the nuclear fuel assembly symmetry axis, and the y-axis is orthogonal to the two . The slit is 100 mm long along the x-axis, and 1.5 mm wide along the y-axis. The dimensions of the slits in the z-directions are 10 mm at the detector window, and 70 mm at the front opening of the collimator slit.

The dimensions of the detector elements are 5*2*5 mm (not to scale in figures 4 and 5). The collimator

walls are made of a tungsten alloy. There is also additional shielding around the detector, and above

the collimator, not included in figure 4. This shielding consists of two steel plates of dimensions

32*2*20 mm, placed vertically along the detector walls. An additional steel plate of dimensions

132*20*42 mm is located on top of the collimator.

Figure 4: Side view of a PGET collimator slit in the xz plane, the detector element is illustrated in blue, to the right.

Figure 5: Top view of a detector element and a collimator slit. The region labelled as 1 is the full- exposure region. Region 2 is the penumbra region, and region 3 is the umbra region.

Nuclear fuel bundles have a typical height of about 4 m. Since they are axially symmetric, it is relevant

to introduce a two-dimensional response function. Based on this approach, the resulting reconstructed

images will also be two-dimensional.

2.2 Line-Spread functions and energy dependencies

At energies under 1.3 MeV, predominant photon interactions are photoelectric effect, Compton scat- tering, and pair production.

• The photoelectric effect dominates at low-energy and high atomic number Z. In this process, the incoming photon energy is completely absorbed by an atom, and an electron is emitted.

• Compton scattering is the interaction of an incoming radiation with an electron in the absorbing material. The energy of the incoming photon is partly absorbed. The photon is deflected and an electron is ejected. Compton scattering is the dominant interaction process for photon energies from 100 keV to 10 MeV.

• Pair production occurs at energies higher than 1.022 MeV. When a photon passes nearby the nucleus of an atom, it may interact with its magnetic field by disappearing and forming a pair of one electron and one positron.

Since the detector is made of equal amounts of Cadmium (Z=48) and Telluride (Z=52), its effective atomic number is Z=50. According to the following diagram, the dominant interactions for energies under 1.3 MeV are photoelectric effects and Compton scattering in the detector.

Figure 6: Dominant interactions for different energies and atomic numbers [6]

In order to obtain the spatial response of PGET, the response of a detector element to a line source, parallel to the z-axis and thus also parallel to the fuel assemblies’ symmetry axis, has been simulated at various source positions, as described in section 3. Of particular interest were the regions located in the regions labelled 1, 2 and 3 in figure 5. Region 1 is the region of full exposure. Photons emitted by a source in this region may reach the whole detector element undisturbed, without passing any distance through the collimator material. Region 2 is the penumbra region, which may only expose part of the detector without passing any distance through the collimator material. Region 3 is the umbra region, were gamma rays must pass the collimator material to reach the detector.

A perfect collimator material would only allow radiations from the full-exposure and penumbra regions

to hit the detector. However, in authentic materials, gamma-ray transmission through the material

will occur. This effect is called septal penetration, since it primarily occurs through short distances by the corners of the collimator slit. In addition, gamma-ray scattering into the detector from gamma- ray interactions into the collimator material will occur. Both these components will give additional contributions to the detector response in all three regions, and in the umbra region (3 in figure 5), they will be the only contributions to the measured radiation.

Since gamma rays emitted by Europium-154 have a higher energy than those emitted by Cesium-137, they are less attenuated when crossing the collimator material. Accordingly, septal penetration is higher for the Europium photons than for Cesium photons.

While direct exposure and septal penetration imply that photons will hit the detector element with preserved energy, one should note that the detected signal will also depend on the energy response of the detector element; few gamma rays deposit all their energy in the detector material and a large fraction will undergo processes like Compton scattering. Consequently, only part of the incoming radiation energy may be deposited in the detector element. Furthermore, gamma-ray scattering into the detector from surrounding materials will add to the complexity of the energy distribution of the recorded detector pulses. Accordingly, the spatial response will depend on the energy windows used when counting the number of events in the detector elements.

The PGET device typically uses the following energy windows:

• 400-600 keV

• 600-700 keV

• 700-1000 keV

• 1000 keV and upwards

Due to the complexity of the response of the PGET instrument, Monte Carlo simulations offer an adequate means to explore it, as further described below.

3 Materials and Methods

The Monte Carlo N-Particles software package (MCNPX) [7] was used for all simulations of gamma-

ray transport in this work. MCNPX is based on the Monte Carlo method. This algorithmic method

aims at calculating a value using physics processes which are stochastic by nature. In this method,

the rules of propagation of particles in a material are described by probabilities of different types

of interaction with atoms in the material. This method is statistical, therefore a high number of

source particles must be simulated in order to obtain a statistically confident result. In the case

when relatively inefficient measurement geometries are simulated, such as the PGET geometry, the

simulations may require heavy computations, and time consumption may become an issue. In order

to limit the computation time in this work, the simulations have been executed on a computer cluster,

using 64 processors in parallel. Various techniques have also been used to limit time consumption,

such as limiting the number of particles and introducing cut-off energies, as further described below.

3.1 Simulation procedure

The MCNPX software package requires the definition of cells, surfaces and materials in order to define the geometry of the problem [7]. The source type and the particles considered must also be defined.

This information is provided to MCNPX in an input file. In this work, the execution of MCNPX input files and the analysis of the resulting output files is done using Matlab, as explained in figure 7:

Figure 7: The simulation procedure used in this work

In order to limit simulation time without increasing uncertainties, the universe of the simulation (the outer limit of the geometry defined in the MCNPX input file) was limited to a box of dimensions 46*2*23 cm. This implies that particles travelling outside of this box will no longer be followed.

The angle of emission was also limited so that only gamma rays emitted in a direction within 22 degrees from the negative x-axis were simulated. Altogether, these limitations imply that large- angle scattering into the detector will be limited in the simulations, while small-angle scattering and full-energy transport will be covered more adequately. Accordingly, the results can be expected to realistically represent experimental data for relatively high-energy bins, whereas the results for low- energy bins may not be realistic. However, as noted in section 2.2, PGET typically only records data at energies higher than 400 keV.

The MCNPX software package also offers the possibility to use importance sampling as a variance

reduction method. When the importance of a cell is, for example, set at 3, the software divides each

particle entering the cell into 3 particles. Each of these 3 particles is followed but they are each

considered as one third of a particle when tallying. This method improves statistics in such a cell

with only a minor increase in simulation time. Since the focus of this work is on the energy deposited

in the detector, the importance in the detector cell was set at 3, while other cell had an importance

of 1. When exploring the magnitude of various contributions to the detector counts, some executions

were also made using importance 0 in the collimator material (see section 4.2).

Defining cut-off energies as explained in 3.2.3, also contributed to reducing simulation time.

3.2 Modelling gamma-transport using MCNPX

In this simulation, only the photon interactions described in section 2.2 were taken into accounts.

Doppler effects, photonuclear interactions, coherent scattering and Bremsstrahlung radiations were excluded. An example MCNPX input file is attached in appendix 1. A description of the contents of this file is given below.

3.2.1 Cell and surface definition

The geometry of the device was defined using 7 cells as illustrated in figures 8 and 9 below:

• The cell in blue is the collimator material, i.e. a box of tungsten that encloses the collimator slit.

• The collimator slit, in red, is filled with air.

• The green cell is the detector cell. It represents a Cadmium Telluride detector element.

• The yellow cells represents plates of steel which are part of the construction material that supports the structure. 2 thin plates of steel are situated around the detector, parallel to the z-axis. A third plate is placed on top of the collimator.

One additional cell, which surrounds the other cells and the source, defines the ”universe” of the MCNPX simulations.

Figure 8: Cell definitions as seen from above (x-y plane) including collimator material (blue), slit

opening (red), steel (yellow) and detector (green). The nuclear fuel bundles are placed along the z-axis

during measurements.

Figure 9: Cell definitions as seen from the side (x-z plane), including collimator material (blue), slit opening (red), steel (yellow) and detector (green). The red cell above the steel plate represents the air surrounding the device. During measurements, the detector heads rotate along the x-y plane.

3.2.2 Material definition

Four materials are of interest in these simulations of the PGET device. Their composition as well as their density are required for MCNPX to calculate how photons are transported in these materials.

• Air is present in the collimator slit and around the device. Air was described as a mixture of Nitrogen (78 %), oxygen (21%) and argon (1%). The density of the air was set to 0.00129 g per cm

• The collimator walls are made of a tungsten alloy, this material has a density of 18 g per cm

. The composition of this tungsten alloy was 96 % tungsten, 3 % Nickel and 1 % Copper.

• The detector is made of Cadmium Telluride. The density of CdTe was assumed to be 5.85 g per cm

. This material consists of equal proportions of Cadmium atoms and Tellurium atoms

• The steel plates of the device were simulated to be made of Iron (68%), Chromium (18%) and Nickel (14%). The density of steel was defined as 8.05 g per cm

.

In the MCNPX input file, the materials can either be described by their composition in terms of weight

fraction or in terms or atomic fraction. Here, all materials were described by their weight fraction

except for the detector which was described in terms of atomic fraction. (The former compositions

are defined using minus signs in the data cards as seen in the input file in appendix 1).

3.2.3 Source particles, energy bins and cut-off

For both simulations, the number of particles simulated was set at 2 billions. This number was selected by evaluating the trade-off between uncertainty and time consumption. A higher number of particles would give better statistics. However, the simulation would take longer. The source is further defined at the event of executing the MCNPX files using Matlab, as described in section 3.3.

The parameter of interest in this MCNPX simulation is the energy deposited in the CdTe detector for each source position. Since the energy of gamma-rays emitted from Cesium is 662 keV, the energy deposited in the detector for each event will be between 0 and 662 keV. Here, the energy deposition in the detector was tallied using an F8 tally, divided into 14 energy bins between 0 and 700 keV, where each bin was 50 keV wide. Similarly, the energy deposited into the detector for the simulation with Europium-154 will be between 0 and 1274 keV. Energy deposition events were tallied in the same way as for the Cesium simulation, with 26 energy bins of 50 keV width between 0 and 1300 keV. When extracting the data using Matlab, as described in section 3.4, bins were summed corresponding to the energy windows used in PGET (see section 2.2).

In order to reduce simulation time, one can instruct the software to stop following a photon when its energy is lower than a certain value. This cut-off energy depends on the properties of the simulated device. Here, a cut-off value of 100 keV was used for photons. Photons entering the detector with energies below 100 keV are not of interest in PGET. Furthermore, photons released in the detector at such low energy are very unlikely to escape the detector. The latter can be verified by calculating the half-value layer of 100 keV photons in Cadmium-Telluride:

The formula for the half value layer [8] is:

HV L = ln2

µ (1)

According to the NIST database [9], the mass attenuation coefficient for 100 keV photons in Cadmium- Telluride is

=1.671 cm

/g. Cadmium-Telluride has a density of ρ = 5.85 g/cm

. The attenuation coefficient of 100 keV photons in the detector is thus:

µ = 1.671 ∗ 5.85 = 9.775cm

Accordingly, the HVL is:

HV L = ln2

9.775 = 0.0709cm

The HVL is thus smaller than the detector dimensions. This means that photons created in the detector at an energy lower than 100 keV are not likely to escape it. Accordingly, the use of a cut-off energy at 100 keV is justified.

The cut-off energy for electrons was set at 300 keV because electrons are heavier particles that do not travel long distances. Indeed, according to the ESTAR database [10], the CSDA range of 300 keV electrons in Cadmium-telluride is 0.1487 g/cm

. Since the density of the Cadmium Telluride material in this study is 5.85 g/cm

, the distance travelled by 300 keV electrons in the detector will be 0.1487 * 5.85 = 0.025 cm = 0.25 mm. This is smaller than the detector’s dimensions (5*2*5 mm).

It is accordingly unlikely that electrons released inside the detector with such a low energy escape the

detector, justifying also this cut-off energy.

3.3 Execution of MCNPX simulations using Matlab

Since Monte Carlo methods are based on statistical processes, the simulations must be run for a large number of source particles in order to obtain good statistics. The number of particles simulated for each positions in this study was 2 billions and the number of source positions was 552. It was thus necessary to run the simulation on several processors to achieve acceptable computation times.

Accordingly, the simulations were executed on a computer cluster, enabling the use of 64 processors in parallel.

A Matlab function (attached in Appendix 2) was created, which calls the MCNPX software and runs the simulations in one batch, sequentially for all 552 source positions. For each position, an output file is created, containing the number of energy deposition events having occurred in the detector for each energy bin.

The source was positioned in 16 positions along the x-axis (see the geometry definition in section 5) from x= 20 mm to x= 320 mm. Along the y-axis, it was placed in steps of 0.15 mm from y= 0 mm to y= 8.40 mm. When simulating the full response function, the number of y-positions was increased with distance in x to cover the gradually expanding regions of full exposure, penumbra and umbra illustrated in 5.

3.4 Extraction of MCNPX tally data using Matlab

The MCNPX output files come in different formats, depending on the input files and how tallying is executed. In order to extract the wanted information from the 552 output files, a Matlab function was written. This function reads each file and extracts the energy response and uncertainty for each energy bin. The results were placed in a Matlab array of variables, which could be used for analyzing and plotting.

4 Results

4.1 Gamma-ray energy spectra

Figure 10 presents a spectrum of a Cesium-137 source obtained through the simulation of the PGET

device with MCNPX. This spectrum is similar in shape to a spectrum obtained experimentally with

the PGET device, which is presented in figure 11. The similarities of the two spectra give confidence

in the simulations, in particular for energy regions above 400 keV, which constitute the tomographic

data collected by PGET (see section 2.2). At lower energies, the simulations underestimate the count

rates, which may be expected due to the limitations of the simulations according to section 3. The

x-ray peak (labeled 1 on figure 11) visible on the experimental spectrum is not visible on the simulated

spectrum because these types of interactions are not taken into account by MCNPX. The agreement

between simulated and experimental spectra is further discussed in section 5.

Figure 10: Simulated spectrum of a source of Cesium-137. The part of the spectrum labeled 2 is the Compton edge and the peak labeled 3 is the photopeak.

Figure 11: Experimental spectrum of a source of Cesium-137 collected using the PGET device. The

peak labeled 1 is the x-ray peak. The part of the spectrum labeled 2 is the Compton edge and the peak

labeled 3 is the photopeak. By courtesy of T. White, IAEA.

4.2 LSFs for a source at x=160 mm

A simulation on a reduced number of positions (57 in total) was run for evaluation purposes. It covered all grid points at x=160 mm, from y=0 to y=0.84 mm. One goal of this simulation was to estimate the number of source particles needed to obtain significant statistics. In addition, these simulations were used to investigate (1) at which distance from the collimator slit axis (i.e. at which values in y) the detector response is small enough for the response function to be truncated, based on the decrease in number of detector counts from a source positioned at a certain y value as compared to the counts from a source at y=0, and; (2) evaluate the magnitude of count rate contributions from septal penetration and scattering (see discussion in section 2.2) as compared to direct exposure.

4.2.1 LSFs for the source of Cesium-137

Figure 12: Counts based on energy deposition in the detector (in arbitrary units), for different source positions along the y axis, for the energy deposition events between 400 and 600 keV (Blue curve) and and between 600 and 700 keV (red curve). Error bars represent the stochastic uncertainty due to Monte Carlo characteristics at the 1 sigma level.

The spatial response function for a source at x=160 mm is presented in figure 12. First, it was verified that the lateral cut-off of the response in the y-direction used in the tomographic analysis codes was adequate. At x=160 mm, the lateral cut-off had been defined at 0.48 mm, and these simulations show that a source in this position contributes with a detector count rate 300 times smaller than the count rate of the same source in the position y=0 when using the detector window 600-700 keV. For the 400-600 keV energy window, the corresponding number is 140 times.

Second, the statistics of the simulation was evaluated. The simulation time for 2 billion particles was

approximately 10 minutes per source position, and it was found that the relative uncertainty of all

individual data points at y <0.5 mm was smaller than 3.5 %, which was considered acceptable. In

addition, the MCNPX software runs a number of statistical checks to inform the users of the reliability

of the simulation’s results. Here, all 10 statistical checks were passed for the data at y < 0.5 mm. It

was concluded that the number of particle was high enough for the purpose of this simulation.

4.2.2 LSFs for the source of Europium-154

Figure 13: Counts based on energy deposition in the detector (in arbitrary units), for different source positions along the y axis, for the energy deposition events between 700 and 1000 keV (Blue curve) and and between 1000 and 1300 keV (red curve). Error bars represent the stochastic uncertainty due to Monte Carlo characteristics at the 1 sigma level

The spatial response for the source of Europium-154 at x=160 mm and varying y positions is presented in figure 13. For this gamma-ray energy (1274 keV), the count rate for a source at y=0mm is about 32 times larger than the count rate for a source at 0.48 mm in the energy window 700-1000 keV and respectively about 110 times larger for the energy window 1000-1300 keV. In this case, a lateral cut-off at 0.48 mm was also justified. Furthermore, statistics were considered to be generally even better than for Cesium-137.

4.2.3 Direct exposure vs septal penetration and scattering

As explained in section 3.1, it is possible to set an importance to a material in order to reduce variance.

By setting importance 0 to a cell, photons entering this cell are no longer followed by the MCNPX software. Here, such a procedure has been used in order to analyze the influence of septal penetration on the energy deposition in the detector.

In figure 14, results are presented from simulations that were executed with two importance options

for both sources. The energy windows considered were 600-700 keV for the source of Cesium and

1000-1300 keV for the source of Europium. In the first simulation, all cells were defined with an

importance of 1. And in the second simulation, the collimator walls were defined with an importance of 0, so that all photons entering the collimator walls were discarded by the software.

Figure 14: Simulation for Cesium (green and blue curve) and Europium (red and magenta curves).

The green and magenta curves were run with importance 0 on the collimator material cell. This means that septal penetration and scattering in the collimator material were excluded in these simulations, offering opportunity to analyze their relative contributions relative to the direct exposure of the detector.

4.3 Full response function

4.3.1 LSFs for the source of Cesium-137

The spatial response of the detector for sources at various distances x from the collimator slit opening

was obtained by simulating a set of 552 source positions. Some results for Cesium-137 are presented

in figure 15. The energy window selected here was 600-700 keV. The origin of the diagram (x=y=0)

represents the center of the collimator slit and the response is plotted for distances source-detector

between 20 and 180 mm. Since the response is expected to be symmetric on both sides of the y-axis,

only positive values of y were represented. Here, the source positions on the y-axis vary between 0 and

6 mm. As expected, response of the detector decreases as the source moves away from the detector,

as well as away from the collimator slit axis.

Figure 15: An example of full Line-Spread Functions and associated uncertainty for the PGET device for a source at various distances x from the detector. The energy of the source in this simulation is

21

In the lateral direction (y), a sharp decrease can be observed at y=0.75 mm. When the source is moved beyond that plane, it enters the penumbra region, as explained in figure 5. When the source leaves the penumbra region and enters the umbra region, there are still contributions to the detector response from septal penetration and scattering.

As expected, the uncertainty increases with the distance to the detector and with the distance from the collimator slit axis (see right part of figure 15). The diagrams are cut at the limit of the response function defined in the tomographic analysis code, and uncertainty increases rapidly when approaching this limit. Indeed, the limit is in the umbra region so that particles must pass the collimator material to reach the detector, and particles entering the collimator wall material have a low probability of reaching the detector. Therefore, counting statistics for a source near this limit is poor. However, for most source positions within the full exposure region, the uncertainty is low, showing that the number of particles simulated was large enough for statistics to be sufficient in this region.

4.3.2 LSFs for the source of Europium-154

Gamma-rays emitted from the Europium-154 source are more penetrating as compared to those emit- ted from Cesium-137. This means that the contribution from gamma-ray emitted by a source in the umbra or penumbra region should be higher than for photons emitted by the source of Cesium.

However, the detector element will also be less efficient in detecting the gamma-rays hitting it, due to the higher energy. An example of the response function for Europium-154 is presented in figure 16.

Here, the energy window between 1000 and 1300 keV has been selected for the analysis.

Figure 16: An example of a Line-Spread Function and associated relative uncertainty for the PGET

device. The energy of the source in this simulation is 1274 keV and the energy bin plotted is 1000 to 23

One may note that the uncertainties for the source of Europium are generally lower than for the source of Cesium in the umbra and penumbra regions due to the higher transmission through the collimator material and thus the larger number of events in the detector. In the full-exposure region, though, the uncertainties are somewhat larger due to the lower full-energy detection efficiency at higher gamma-ray energy.

5 Discussion

5.1 Gamma-ray energy spectrum

As seen in figures 10 and 11, the gamma-ray energy spectrum of the Cesium source obtained in this simulations is in agreement with the spectrum obtained experimentally by the PGET device for energies used in tomographic assessment (> 400 keV).

However, at lower energies, discrepancies occur. This can be explained by the approximations made when defining the geometry in the simulations. The PGET device is a complex instrument. Its structure comprises many elements that were not included in the simulation for simplification purposes.

Since the geometric approximations used primarily will exclude large-angle scattering into the detector, which would give rise to detector events at relatively low energies, these approximations will imply an underestimation of the count rates at low energies. Still, since tomographic assessment is made at energies > 400 keV (see section 2.2), the discrepancies at lower energies can be considered acceptable.

5.2 Lateral cut-off of the spatial response function

The LSFs obtained in this work will constitute an important piece of information to the algebraic tomographic reconstruction codes used when analyzing data from PGET. However, to obtain higher efficiencies in the reconstructions, the LSFs will be truncated in the y-direction. In order to evaluate the validity of this truncation, the count rate at the selected lateral cut-off was compared to the count-rate at y= 0 mm, as presented in section 4.2. As seen in figure 12, the energy deposited when the source is located at y= 0.48 mm is 300 times inferior that the energy deposited at y=0 mm for the source of Cesium in the energy bin 600-700 keV. For the source of Europium, the energy deposited when the source is located at y= 0.48 mm is 110 times inferior that the energy deposited at y=0 mm in the energy bin 1000-1300 keV. Truncating the LSFs at the selected y position is thus a valid approximation.

5.3 Contribution from septal penetration

Simulating the spatial response of the PGET device also opened up a possibility to study different contributions to the response. In this work, contributions from direct exposure were isolated by setting importance 0 to the collimator material in a separate set of simulations (see section 4.2.3):

Consequently, the magnitude of other contributions such as septal penetrations and scattering in the

collimator material could be analyzed. A previous study has been performed on a device called UGET

[3]. This device is similar to PGET but its collimator slit openings are not tapered. In the study

performed on the UGET device, it was found that 10-15% of the count rate comes from the gamma

rays subject to septal penetration, while in this work, contributions from septal penetration in PGET were found to be around 3% for the source of Cesium and 6% for the source of Europium.

These lower values can be explained by the fact that the collimator in the PGET device is tapered in the (x,z) plane and the detector element is smaller in the z-direction than the back opening of the collimator slit. Accordingly, gamma rays travelling inside the tungsten material below or above the detector are therefore very unlikely to hit the detector, as illustrated on figure 17.

Figure 17: Septal penetration in the (xz) plane

6 Conclusion

Enhanced techniques for partial-defect verification is one of the focus area identified in the IAEA’s long term R&D plan. In 2013, this institution stated the need for ”more sensitive and less intrusive alternatives to existing non-destructive assay instruments” [11]. The PGET device has been developed for the purpose of partial-defect verification. In this work, the spatial response of this device was calculated in order to be incorporated in tomographic image reconstruction algorithms. The values calculated by MCNPX had low uncertainty (< 3, 5 %), which is considered low enough for these data to be used.

The contribution from septal penetration and scattering in the collimator material was found to be about 3 % for the source of Cesium and 6 % for the source of Europium, which is much lower than the corresponding contributions for the UGET device. Since the sensitivity to these effects is smaller for PGET than for UGET, a simplified model that takes only direct contributions into account, may represent the response function with decent accuracy. Future studies may be performed based on the tomographic reconstruction of data using the two models to evaluate the performance of such a simplified model.

7 Acknowledgement

I would like to thank all the people that supported and guided me during this project.

Firstly, I would like to thank my supervisors. Staffan Jacobsson Sv¨ ard’s mentoring and encourage-

ments have been especially valuable during this research project. His meticulous scrutiny of my report and his attention to details have helped me improving my writing skills. In addition, thanks to Peter Jansson, who patiently guided me in using the MCNPX software package and the operating systems used at Uppsala University.

Lastly, I would like to thank Peter Andersson, for providing the template MCNPX input files and

Matlab functions this work was derived from.

References

[1] United Nation Office for Disarmement Affairs. The non-proliferation treaty.

www.un.org/disarmament/wmd/nuclear/npt/. Accessed: 2017-04-30.

[2] Levai F. Turunen A. Berndt R. Vaccaro S. Schwalbach P. Honkamaa, T. A Prototype for passive gamma emission tomography. IAEA-CN–220, 46(26):287, 2015.

[3] Jacobsson Sv¨ ard S. Mozin V. Jansson P. Miller E. Honkamaa T. Deshmukh N. White T. Wittman R. Trellue H. Grape S. Davour A. Vaccaro S. Andersson P. Smith, L.E. A viability study of gamma emission tomography for spent fuel verification: JNT 1955 phase I technical report. Technical report, Pacific Northwest National Laboratory, USA, 2016. PNNL-25995.

[4] G. Knoll. Radiation detection and measurement, chapter 13. John Wiley and sons, 2000.

[5] S. Jacobsson Sv¨ ard S, Holcombe and S. Grape. Applicability of a set of tomographic reconstruc- tion algorithms for quantitative spect on irradiated nuclear fuel assemblies. Nuclear Instruments and Methods in Physics Research Section A, 783:128–141, 2015.

[6] Ritenour ER. Hendee WR. Medical imaging physics, chapter 4, page 63. Mosby-Year Book, 3rd ed edition, 1992.

[7] Radiation Safety Information Computer Center (RSICC). MCNPX User’s Manual. Los Alamos National Laboratory, Los Alamos, New Mexico, version 2.4.0 edition, 9 2002.

[8] E. Podgorsak. Radiation physics for medical physicists, chapter 7. Springer, 2010.

[9] National Institute For Standards and Technology. X-ray mass attenuation coeffient for cadmium telluride. http://physics.nist.gov/PhysRefData/XrayMassCoef/ComTab/telluride.html.

Accessed: 2017-04-24.

[10] National Institute For Standards and Technology. CSDA range of 300 keV electrons in Cadmium Telluride. http://http://physics.nist.gov/cgi-bin/Star/e_table.pl. [Accessed 2017-06- 12].

[11] IAEA. Department of safeguards long-term r&d plan, 2012-2023. Technical report, International

Atomic Energy Agency., Vienna, Austria, 2013. STR-375.

8 Appendix

8.1 Appendix 1: MCNPX input file (emilie.i)

c Response function calculation for PGET c

c =============================

C source threads in Z axis are iteratively changed in separate files C see call to file Rst1.rc below

c =============================

c --- CELL CARDS --- c =============================

c C

100 3082 -18.0 -300 320 310 imp:n,p,e=1 $main shielding box 110 4852 -5.85 -320 imp:n,p,e = 3 $ CdTe in det

130 999 -0.00129 -310 imp:n,p,e = 1 $ air in collimator 125 0 -290 300 320 350 360 370 380 imp:n,p,e = 1 $ air around 150 2624 -8.05 -350 imp:n,p,e = 1

160 2624 -8.05 -360 imp:n,p,e = 1 170 2624 -8.05 -370 imp:n,p,e = 0 180 999 -0.00129 -380 imp:n,p,e = 0 120 0 290 imp:n,p,e=0

C c

c =============================

c --- SURFACE CARDS --- c =============================

C

290 box -13.2 -1 -5 46 0 0 0 2 0 0 0 23 $ air

300 box 0 -1 -5 -10 0 0 0 2 0 0 0 10 $ tungsten ytterbox

310 ARB 0 0.075 3.5 -10 0.075 0.5 -10 0.075 -0.5 0 0.075 -3.5 0 -0.075 3.5 &

-10 -0.075 0.5 -10 -0.075 -0.5 0 -0.075 -3.5 1234 5678 2376 1485 1265 3784 $ Collimator 1 320 box -10 -0.1 -0.25 -0.5 0 0 0 0.2 0 0 0 0.5 $ detector

350 box -10 -0.1 -1 -3.2 0 0 0 -0.2 0 0 0 2 $ steel plate 1 360 box -10 0.1 -1 -3.2 0 0 0 0.2 0 0 0 2 $ steel plate 2 370 box 0 -1 5 -13.2 0 0 0 2 0 0 0 4.2 $ steel plate above

380 box 0 -1 9.2 -13.2 0 0 0 2 0 0 0 8.8 $ air above upper steel plate c

c =============================

c --- DATA CARDS --- c =============================

c

c Tungsten alloy, from Peter’s template m3082 28058 -.02383

28060 -.00918

28061 -.0004 28062 -.00127 28064 -.00033 29063 -.01038 29065 -.00462 74180 -.00114 74182 -.25175 74183 -.13595 74184 -.29108 74186 -.28009 c CdTe (detector)

m4852 48000 0.5 52000 0.5 c Steel

m2624 26000 -0.68 24000 -0.18 28000 -0.14 c Air

m999 8000 -0.232 7000 -0.7547 18000 -.0133 c

MODE P

PHYS:P 100 1 1 0 1 CUT:P 0.1

CUT:e 0.3 NPS 2e9 c c

c Cell tally for energy deposition F8:P 110

E8 0.3 19i 1.3 Read file=Rst1.rc c SDEF erg=1.274436 c SDEF erg=0.661657 sp2 0 1

sp4 0 1 c

c Mesh tally: used for plotting during tests c tmesh

c rmesh91:P

c cora91 -10 40i 33

c corb91 -1 10i 1

c corc91 -5 20i 18

c endmd

8.2 Appendix 2: Matlab execution file

function res = Csrun()

% Response functions for the PGETcollimator

% Defining grid for source locations x= 0:20:320;

y = 0:0.15: 8.40;

% Defining z length as function of x, (same number of elements as x) z =44.4:8.4:170.4;

% Mcnpx uses cm x=x/10;

y=y/10;

z=z/10;

% emission angle with some margin theta = atan(42/100)

% Defining end index of y values for the loop below endy = 12:3:57;

for i = 11:length(x)

for j = 1:endy(i); % length(y)

% Text in Read Card file

RCtext = [’SDEF par=p ext=d2 erg=0.661657 axs= 0 0 1 \n VEC = -1 0 0 DIR = d4\n pos=’ , num2str([x(i),y(j)]),’ 0

\nsi2 H ’,num2str([0, z(i)]),’\nsi4 ’,num2str(cos(theta)),’ 1\n’];

% Name of output files ( + "o" in the end) OUTtext = [’C’,num2str(i),’_’,num2str(j)];

%Create file

fir = fopen(’rst1.rc’, ’w’) fprintf(fir, RCtext)

fclose(fir)

% Run MCNPX - 64 parallel jobs

dos([’mpirun -np 64 mcnpx i=emilie.i n=’, OUTtext]) % end

end

8.3 Appendix 3: Matlab extraction file

function [LSF, unc] = GetLSF()

% This function extracts the simulated spatial response from MCNP results files

% ********************************************

% Set geometric variables identical to the MCNP execution

% ********************************************

x= 20:20:320;

y = 0:0.15: 8.40;

z = 44.4:8.4:170.4;

ymax = 8.40 zmax = 170.4;

endy = 12:3:57;

% x= 20:20:180;

% y=0:0.15: 5.25;

% z= 44.4:8.4:111.6;

theta = atan(42/100);

x=x/10;

y=y/10;

z=z/10;

% ********************************************

% Define fraction of solid angle covered by emission.

sa_fraction = 0.5*(1-cos(theta)) for i = 1:9 %length(x)

for j = 1:endy(i)

% ********************************************

% Set filenames identical to the MCNP execution

filename = [’E’,num2str(i),’_’,num2str(j),’o’]

% ********************************************

% Stop execution if file does not exist if exist( filename, ’file’) == 0

error([’ No such file: ’,filename])

else

% Read file if it exists

% ********************************************

% Call adequate read function for the tally you want.

[energy,resp,fel] = readF8tally(filename);

% ********************************************

try

% ********************************************

% Sum energy bins you want to use.

% (Start by using one bin, e.g. (8:8).

temp_energy= energy(8);

temp_resp = (resp(8))+ (resp(9));

temp_fel = sum(fel(8:9));

% ********************************************

% Weight result with z-length and solid angle LSF(i,j) = temp_resp * z(i) * sa_fraction;

% ********************************************

% Calculate uncertainty.

% (Check: How is uncertainty expressed?)

%temp_fel = sum(unc(2:2)*resp(2:2))/temp_resp temp_fel =

sqrt(sum((fel(8:9).*resp(8:9)).*(fel(8:9).*resp(8:9))))/temp_resp unc(i,j) = temp_fel;

% ********************************************

%catch

% error([’ Interpretation error in ’,filename]) end

end end end

[x1,y1]=meshgrid(0:0.15:5.25,20:20:180);

% Present the data in a figure.

subplot(1,2,1);

surf(x1,y1, LSF);

title(’LSF for x between 20 and 180 mm (energy bin:600-700keV)’,’FontSize’,12,’FontWeight’,’bold’) zlabel(’Energy deposited in the detector in arbitrary units’,’FontSize’,12,’FontWeight’,’bold’) xlabel(’position on the y axis’,’FontSize’,12,’FontWeight’,’bold’);

ylabel(’Position on the x-axis’,’FontSize’,12,’FontWeight’,’bold’);

subplot(1,2,2);

surf(x1,y1,unc);

title(’Uncertainty’,’FontSize’,12,’FontWeight’,’bold’);

xlabel(’position on the y axis’,’FontSize’,12,’FontWeight’,’bold’);

ylabel(’Position on the x-axis’,’FontSize’,12,’FontWeight’,’bold’)