Analysis of Secondary Particles as a Complement to Muon Scattering Measurements

Abstract

Cosmic ray tomography is an emerging imaging technique utilizing an ambient source of radiation. One common tomography method is based on the measurement of muons scattered by the examined objects, which allows the reconstruction and discrimination of materials with different properties. From the interaction of air shower particles induced through cosmic rays with the material to be scanned, secondary particles, predominantly photons, neutrons and electrons, can be produced, which carry complementary information about the objects and their materials. However, this information is currently not fully exploited or only studied in coincidence with the incoming air shower particles. Therefore, this work presents a novel approach utilizing only the information from secondary particles to reconstruct and discriminate objects made out of a variety of materials. It also includes a detailed analysis of the kinematics of secondary particles and their dependency on material characteristics. In addition, a reconstruction algorithm to produce 3D maps of the examined volume from the measurement of secondary particles is introduced. This results in a successful reconstruction and differentiation of objects in various geometrical compositions.

1. Introduction

1.1. Muon Scattering Tomography and Its Use Cases

Non-destructive imaging techniques have become an essential tool to increase the level of safety and security within our society. Many of these techniques use radiation from artificial sources, such as X-rays [1,2], UV-light [3] or microwaves [4,5]. By contrast, the emerging method of cosmic ray tomography utilizes a natural source of energetic particles induced by cosmic rays [6,7]. Stable, charged particles from astrophysical sources, mainly protons, interact with atoms and molecules at high altitude in earth’s atmosphere and produce a cascade of high-energy ions and particles towards the surface, so-called air showers. At sea level, these air showers contain predominantly muons, electrons, photons, and neutrons [8,9].

Muons from air showers have an average energy of 4 GeV and a flux of around 1 cm⁻²min⁻¹ at sea level [10]. They have, unlike electrons, photons and neutrons, a high material penetration power, as they have a nearly 200 times higher mass than an electron and mostly interact with atoms through Coulomb scattering. This makes muons induced by cosmic rays a perfect natural source for tomography systems [11]. Muon Scattering Tomography takes advantage of the deflection of the muon trajectory due to the Coulomb scattering, with the deflection angle being proportional to the atomic number Z and the density $\rho$ of the material. By measuring the scattering angles for many muons, the properties of an examined object can be reconstructed, allowing to differentiate between materials with distinct atomic numbers and densities [12,13]. Another technique of utilizing cosmic ray muons for imaging of large-scale structures is the so-called Muon Radiography. For this method, the differential muon absorption is measured, which allows for the discrimination between materials, as the absorption is proportional to the atomic number and density of the examined materials [14,15].

Both techniques of Muon Tomography are used in a variety of fields, ranging from the reconstruction of archaeological sites [16,17] over nuclear waste inspection [18,19] and geological studies of volcanoes [20,21] to shipping container scans [22,23].

The development of Muon Scattering Tomography systems has become an important consideration for border security mainly due to its benefits over conventional X-ray technology. The utilization of naturally occurring cosmic rays prevents the need for artificial radiation resulting in a safe and passive scanning technology with relatively short scan times [24,25]. Furthermore, heavily shielded contraband, especially nuclear materials or narcotics, can be reconstructed due to the high penetration power of muons [26,27].

1.2. Complementary Information from Secondary Particles

While passing through an object, muons produce additional particles due to energy loss as a result of the interaction with the surrounding matter. These so-called secondary particles are predominately photons, neutrons and electrons. Similarly produced protons have a low production rate and low penetration power compared to the other particles and are therefore neglected. In addition, other air shower particles, such as photons, neutrons and electrons, also produce secondary particles while interacting with matter [28,29,30,31]. The production processes of secondary particles are dependent on the properties of the object’s material [32,33,34,35,36]. Hence, a kinematic analysis of such particles creates complementary information to the reconstruction through solely Muon Scattering Tomography or Muon Radiography.

The main interactions producing secondary photons are Bremsstrahlung and annihilation. Minor contributions come from inelastic and elastic hadron scattering, as well as the photo-nuclear reaction, muon or hadron capture, and radioactive decays. Secondary neutrons are predominantly a product of inelastic hadron scattering, but also originate from photo-nuclear reaction and muon capture. The origin of secondary electrons is either ionization, the conversion of a photon to an e⁺e⁻ pair or Compton scattering.

Current efforts to utilize secondary particles in cosmic ray tomography are mainly focusing on reconstructing secondary photons, electrons or neutrons in coincidence with incoming air shower muons. While some methods rely only on the measurement of secondary particles to reconstruct the examined objects [37,38,39], other techniques follow a hybrid approach combining the muon scattering and secondary particle measurements for reconstruction [40,41].

1.3. Proposed Work

The aim of this work is to analyze the kinematics of secondary particles detected in a simulated cosmic ray tomography scanning system for shipping containers. As the interaction of secondary particles with the environment limits the range these particles can travel before their detection, the small-scale application of a shipping container scanner is chosen. This work also addresses the question as to whether a reconstruction of the container content using only secondary particle information is possible presenting, a novel and unique approach for a safe and passive scanning technology, which is so far not studied in literature. In addition, the proposed method includes all possible secondary particles resulting from the interaction of all kinds of air shower particles, not just air shower muons. Finally, this work is the first of its kind with a focus on the reconstruction of shipping container content or similar sized objects based on secondary particle measurements. Therefore, the emphasis of this work lies on the following three items.

First, the simplified simulation setup and the analysis of secondary particles within the application of the shipping container scanner is presented. The analysis focuses on the dependency of the reconstructed kinematics of secondary particles on the material properties, particularly the atomic number and the density. A simple method for background reduction from air shower particles and secondary particles from the container walls is shown.

Second, a reconstruction algorithm creating a 3D map of the container content is introduced. This includes the separate detection of secondary photons, neutrons and electrons in the simulation setup. Furthermore, the back-tracing of the detected particles through the container volume, and the final generation of the 3D map based on the different measurements including a background removal method is presented.

Third, the results of various test scenarios are shown. These scenarios are chosen to examine the sensitivity of the reconstruction algorithm to differentiate between objects with distinct material properties and in challenging geometrical compositions.

The paper is organized as follows. Section 2 covers the presentation of the simulation setup, as well as the analysis of the secondary particle kinematics and their dependence on the material properties. Following this part, the different steps of the reconstruction algorithm are shown. The results for the different test scenarios are presented in Section 3. The work is concluded with the discussion of the results and future improvements in Section 4.

2. Materials and Methods

2.1. Simulation Setup

The generation of air showers induced by cosmic rays is implemented using the Cosmic-ray Shower Library (CRY, v1.7) [9]. The altitude is set to sea level and air showers contain muons, electrons, photons, neutrons, protons and pions. The energy of the primary cosmic ray proton ranges from 1 GeV to 100 TeV. The maximum number of particles within an air shower is set to 30 to allow for realistic shower sizes within reasonable computation time. The interaction of particles with matter is simulated with the Geant4 Monte Carlo simulation toolkit (v10.07.p02) [42]. All relevant interaction processes are included, as defined by the internal “Shielding” physics list [43].

The simplified geometry implemented within Geant4 is visualized in Figure 1 and consists of a cubical volume with a side length of 15 m filled with air. The origin of the coordinate system is set to be at the center of the main volume, with the $z-\mathrm{axis}$ pointing in the positive vertical direction. The air showers originate from a 100 m² square area in the $xy-\mathrm{plane}$ located at $z=7.5\mathrm{m}$. The shipping container is modeled by a hollow box with 2 mm thick stainless-steel walls centered at the origin. Its dimensions follow the size of a 20 ft ISO container ($2.44\mathrm{m}\times 6.06\mathrm{m}\times 2.59\mathrm{m}$). Four detector planes surround the container, an upper and lower detector plane in the $xy-\mathrm{plane}$ at $z=\pm 1.75\mathrm{m}$ with an area of $3.5\mathrm{m}\times 7.0\mathrm{m}$, and two detector planes on the sides in the $yz-\mathrm{plane}$ at $x=\pm 1.75\mathrm{m}$ with an area of $7.0\mathrm{m}\times 3.5\mathrm{m}$. The detectors are modeled as 1 nm thick vacuum and show perfect acceptance and reconstruction efficiency to all particles passing through them with perfect spatial and energy resolution. If a particle is detected in one of the four planes, the following information is assigned to the so-called hit:

• Particle type: electron, photon, neutron
• Hit location: $x_h,y_h,z_h$
• Particle momentum: $p_x^h,p_y^h,p_z^h$
• Particle kinetic energy: $E_{kin.}^h$

Within the container, a various number of blocks made out of distinct materials and with different sizes are placed. The exact positions of the volumetric centers of the objects and their compositions used in the following studies are presented in the corresponding sections.

2.2. Analysis of Secondary Particle Kinematics

The production mechanisms for secondary particles are diverse. However, most of the underlying physical interactions are material dependent and therefore, measurements of secondary particles yield information about the properties of the matter, in which they were produced. Hence, it is possible to analyze the kinematics of these particles in order to aid in material discrimination applications. The two main factors affecting the production rate and energy distribution of secondaries are the atomic number Z and the material density $\rho$.

These dependencies are studied by placing two different blocks inside the container, positioned at $x=0\mathrm{m}$, $y=\pm 1.75\mathrm{m}$ and $z=0\mathrm{m}$. A set of 50 million air showers is generated, which is equivalent to an exposure time of around 17 min. After simulating their interactions in the whole volume, all hits generated by neutrons, photons and electrons are reconstructed by registering the kinematic information of the particles, creating three 2D maps of $E_{kin.}^h$ as a function of the hit location, one for the upper detector plane, one for the lower plane and a combined map for both sidewise planes.

To reduce the background contributions arising from primary air shower particles, hits are taken into consideration only if the direction of motion of the particle indicates that the origin of the particle is within the volume enclosed by the detector planes. This technique rejects all hits from air shower particles detected with the upper detector plane, while no hits are rejected in the lower detector plane and 31% of hits from shower particles are eliminated in the sidewise detectors. The container walls, and to a minor degree the air, are additional sources for secondary particles. To decrease this background process, the particle kinematics are analyzed for an empty container scenario with one billion generated air showers. The resulting 2D maps from this setup are then down-scaled to a 50 million generated showers equivalent and subtracted from the 2D maps created with the two blocks setup. The resulting maps show the increased or reduced production of secondary particles from the additional material inside the examined volume.

2.2.1. Kinematic Dependence on Density

To study the impact of different densities on the production of secondary particles, three materials with similar atomic number are chosen: Cesium ($\rho =1.87\mathrm{g}/{\mathrm{cm}}^3,Z=55$), Tin ($\rho =7.31\mathrm{g}/{\mathrm{cm}}^3,Z=50$) and Palladium ($\rho =12.02\mathrm{g}/{\mathrm{cm}}^3,Z=46$). Each element is analyzed separately, with a cube of size 1 m³ located at $y=-1.75\mathrm{m}$ and a thin plate ($1.0\mathrm{m}\times 1.0\mathrm{m}\times 0.1\mathrm{m}$) positioned at $y=1.75\mathrm{m}$. Two blocks with different thickness made out of the same material are studied in order to investigate the effect of self-shielding. If secondary particles are produced within an object, they can also interact with the surrounding matter until leaving the object or getting stopped within. This will affect the number of detectable secondary particles and therefore needs to be taken into account.

Figure 2 shows the distribution of the kinetic energy of photon and neutron hits in the upper detector as a function of the $y-\mathrm{position}$. The cube at $y=-1.75\mathrm{m}$ produces overall significantly more secondary particles than the plate at $y=1.75\mathrm{m}$ as the probability of interaction for air shower particles increases with a higher material volume. For the block and the plate, the number of reconstructed photons increases slightly with increasing density as more atoms become available for potential production of secondary particles. The energy distribution shows a peak at 511 keV coming from annihilation of electrons and positrons, as well as a broad spectrum around this peak originating from Bremsstrahlung. The number of secondary neutrons produced in the block is the highest for the medium density material (Tin), followed by the high density material (Palladium) and the lowest number appearing for the low density element (Cesium). When neutrons are produced in the thin plate, the production rate increases clearly with increasing density. This can be explained with the effect of self-shielding. While a higher density increases the chance for secondary neutron production, it also increases the chance that the produced secondaries interact with the material and get stopped before reaching the detector. If the neutrons are not stopped, they still have a higher chance to lose more energy within the material due to the increased probability of interaction. This can also be seen, as the mean neutron energy decreases with higher density. The self-shielding shows different implications for secondary photons and neutrons, as the underlying processes for their production and their interactions with the surrounding matter are different.

The distributions from the sidewise detector planes in Figure 3 show similar shapes and features as Figure 2. The only difference is that the number of photons from the cube at $y=-1.75\mathrm{m}$ is now higher for the low density material (Cesium) than for the medium (Tin) and high (Palladium) one. This is due to the higher self-shielding in Tin and Palladium and the fact that the path of secondary particles through the material is longer if they reach the sidewise layers of the detector than the upper layer.

In Figure 4 two main regimes are present in the energy distributions measured in the lower detector plane. These two domains can be treated as two independent measurements.

At high energies, around 10 MeV for photons and 100 MeV for neutrons, a deficit of detected particles proportional to the density and the thickness of the object is visible. This deficiency arises because of the absorption of the incoming air shower particles by the examined materials.

At low energies, around $0.5$ MeV for photons and 1 MeV for neutrons, the production of secondary photons appears to decrease with increasing density, as the shielding effect is lower in materials with lower density. The self-shielding is the highest for particles reaching the lower detector as their path through the target material is the longest compared to the upper or sidewise planes. Furthermore, air shower particles lose energy while passing through the objects and therefore, secondary particles created at a later point in time have lower energy and are more likely to get stopped inside the material. The photons reconstructed around $y=0\mathrm{m}$ are the result of the interaction of secondary particles, which were produced within one block, with the other block. Secondary neutrons show an increase in production for materials with higher density, however, the also increasing shielding effect can counteract this behavior, as seen for the bigger cube at $y=-1.75\mathrm{m}$.

Electrons only contribute with the absorption of incoming air shower electrons (see Figure 5), since the container walls are stopping any secondary electron coming from the objects inside. This is also the reason, why no electron distributions are shown for the other detector planes.

2.2.2. Kinematic Dependence on Atomic Number

The dependence of the production of secondary particles on the atomic number is studied by choosing three pairs of materials, with each pair having a similar density:

• Magnesium ($\rho =1.74\mathrm{g}/{\mathrm{cm}}^3,Z=12$) and Cesium ($\rho =1.87\mathrm{g}/{\mathrm{cm}}^3,Z=55$)
• Chromium ($\rho =7.18\mathrm{g}/{\mathrm{cm}}^3,Z=24$) and Ytterbium ($\rho =6.73\mathrm{g}/{\mathrm{cm}}^3,Z=70$)
• Molybdenum ($\rho =10.22\mathrm{g}/{\mathrm{cm}}^3,Z=42$) and Lead ($\rho =11.35\mathrm{g}/{\mathrm{cm}}^3,Z=82$)

For each material pair, 1 m³ cubes are analyzed independently, with the lower atomic number element positioned at $y=-1.75\mathrm{m}$ and the higher atomic number element positioned at $y=1.75\mathrm{m}$.

Figure 6 shows the distribution of the kinetic energy of photon and neutron hits in the upper detector as a function of the $y-\mathrm{position}$ for all three material pairs. The relations with regards to the different densities are similar to the results shown in Figure 2. Secondary photons coming from materials with lower atomic number Z at $y=-1.75\mathrm{m}$ show a lower mean energy and a lower production rate than from materials with high Z at $y=1.75\mathrm{m}$. This is related to the additional positive charge in the nucleus for higher values of Z creating a higher Coulomb force resulting in more Bremsstrahlung photons with higher energies. This effect gets smaller for higher density, as the self-shielding counteracts this effect. In addition, the increased amount of Bremsstrahlung leads also to a higher chance of pair production followed by annihilation increasing the intensity of the peak at 511 keV for higher Z objects. Secondary neutrons show a similar increasing dependence on the atomic number. An atom with higher atomic number also has a bigger nucleus, increasing the chance for secondary neutron production through inelastic scattering processes. An interesting feature visible for Chromium are distinct peaks in the neutron energy, which are hints of dedicated excitation energy levels of the nucleus for specific elements. The distributions for the sidewise layers in Figure 7 are all supporting the observations seen in Figure 6.

The kinetic energy of photon and neutron hits in the lower detector plane are shown in Figure 8. As also seen in the plots from the upper and sidewise planes, the dependencies on the material density are the same as discussed in Section 2.2.1. While the absorption of air shower photons is increasing with higher values of Z, the absorption of air shower neutrons is independent of the atomic number. Since the electromagnetic interaction of photons increases with the higher positive charge of higher Z nuclei, the energy loss for those photons will also be higher resulting in a higher stopping power. However, neutrons can only interact through electroweak and hadronic processes, which have lower cross-sections than electromagnetic interactions, resulting in a lower stopping power for neutrons. Secondary photons show a higher production for higher values of Z only for the low density material, as the shielding effect becomes dominant for higher densities. The rate of secondary neutrons is also increasing proportional with the atomic number, as the nuclei of elements become more unstable with higher Z increasing the chance of neutron emission. The absorption of air shower electrons follows the same principles as for photons and Figure 9 shows therefore also a higher absorption rate for materials with higher atomic number.

2.3. 3D Map Reconstruction

A key component of a tomography system is the reconstruction algorithm. As this work is designed to utilize solely secondary particle data and is independent from muon scattering or absorption based measurements, no relevant reconstruction algorithm is currently available. Therefore, a first, simple procedure has been designed utilizing the measured kinematics of the secondary particles to allow a discrimination between different objects and materials.

2.3.1. Particle Detection

The detection of secondary particles is performed separately for the different detector planes and particle types. Similar to Section 2.2, hits are taken into consideration only if the direction of motion of the particle indicates that the origin of the particle is within the volume enclosed by the detector planes. As also shown in the same section, the material properties impacts the number of reconstructed particles differently for photons, neutrons and electrons, as well as for the different detector planes. While the upper and sidewise detectors can be utilized to detect secondary production from the container content, the lower detector shows the production of low-energy secondary particles as well as the absorption of high-energy incoming air shower particles. The energy boundaries between production and absorption measurements in the lower detector are extracted from Figure 4. Photons with an energy above 1 MeV are used for the absorption measurement, while photons below 1 MeV are part of the production detection. Similarly, neutrons with an energy above 2 MeV are considered for the absorption pattern, neutrons below 2 MeV are declared as part of the production processes. Because the detection of secondary electrons is negligible, this particle type is only taken into account for the absorption measurement in the lower detector plane. Overall, this results in nine separate measurements, labeled as $M1$–$M9$ as shown in

Table 1. The nine separate measurements for the different particle types and detector planes.

	Photons	Neutrons	Electrons
Upper detector	$M1$	$M2$	–
Sidewise detectors	$M3$	$M4$	–
Lower detector—production	$M5$	$M6$	–
Lower detector—absorption	$M7$	$M8$	$M9$

.

2.3.2. Particle Back-Tracing

Using the hit information for each measurement $M1$–$M9$, the particle can be traced back from the position of the hit through the volume enclosed by the detector planes. This is done by starting at the position of the hit $(x_h,y_h,z_h)$ and using the momentum vector $(p_x^h,p_y^h,p_z^h)$ to project the possible trajectory of the particle assuming a straight line path through the volume. As the exact origin of the secondary particle is uncertain, every position on the projected path is equally likely to be the point of production.

The examined volume is discretized into cubical voxels of arbitrary size. To assess, through which voxel the particle trajectory runs, a ray-tracing algorithm [44] is used. Every traversed voxel is marked and assigned a score $s_h=1$. The score is summed up over all number of hits N within each measurement to generate a total voxel score $s_{tot}$:

$$s_{tot}=\sum _{h=1}^N{s}_h$$

(1)

An example of the back-tracing process reduced to 2D for the sake of clarity is visualized in Figure 10.

2.3.3. 3D Map Generation

In addition to secondary particles from the container content, air shower particles and secondaries produced by the container and the air are included in the calculation of the total voxel scores $s_{tot}$. To reduce those background processes, all measurements $M1$–$M9$ are repeated for an empty container. Afterwards, the difference between the voxel scores for the empty and loaded container for each measurement is defined as the final voxel score $s_{fin}$:

$$s_{fin}=s_{tot}^{loaded}-s_{tot}^{empty}$$

(2)

If this score is positive, an increased secondary particle production is measured ($M1$–$M6$), while a negative score represents the absorption of air shower or secondary particles, which is the case for measurements $M7$–$M9$. To allow a combination of all measurements, the final voxel scores of $M7$–$M9$ are transformed to a positive definite quantity $s_{fin}^{\prime }$. This is done by finding the highest value $s_{high}$ separately in each measurement and subtracting each voxel score $s_{fin}$ from it:

$$s_{fin}^{\prime }=s_{high}-s_{fin}$$

(3)

Since the exact origin voxel of a secondary particle is unknown, every voxel on the back-traced particle trajectory is marked. This results in noise in the reconstruction, visible as multiple low-scoring voxels. To remove this noise, a minimum global score threshold $t_{min}$ is applied to the final voxel scores using the maximum final voxel score $s_{max}^{\left({}^{\prime }\right)}$ for each measurement with the following logic resulting in the cleaned voxel score $s_{clean}$:

$$s_{clean}=\left\{\begin{array}{cc}{s}_{fin}^{\left({}^{\prime }\right)}& \mathrm{if}{s}_{fin}^{\left({}^{\prime }\right)}>s_{max}^{\left({}^{\prime }\right)}*t_{min}\hfill \\ 0& \mathrm{otherwise}\hfill \end{array}\right.$$

(4)

The values of $t_{min}$ are manually defined to distinguish between different materials and objects within the container and can also differ between $M1$–$M9$.

The final 3D voxel map is created by taking any combination of the nine measurement planes. The combination is performed by summing up the cleaned voxel scores in each voxel over all utilized measurements to a combined score $s_{comb}$. However, this summation is only done, if all voxel scores from all utilized measurements are non-zero. If any voxel scores from any measurements are zero, the combined voxel score is also set to zero:

$$s_{comb}=\left\{\begin{array}{cc}\sum _{meas.}{s}_{clean}& \mathrm{if}\mathrm{all}{s}_{clean}>0\hfill \\ 0& \mathrm{otherwise}\hfill \end{array}\right.$$

(5)

The set of measurements used for the combination is material dependent as each one contains distinct information about the examined object. Together with the settings of the noise thresholds $t_{min}$, this parameter set allows a discrimination between materials and objects with different densities and atomic numbers.

3. Results

3.1. Test Scenarios

To test the findings from Section 2.2 and the reconstruction algorithm proposed in Section 2.3, three scenarios are presented:

• Separate water and lead block
  Each cube has a size of 1 m³. Both blocks are positioned at $x=0\mathrm{m}$ and $z=0\mathrm{m}$, with the water cube at $y=1.75\mathrm{m}$ and the lead cube at $y=-1.75\mathrm{m}$.
• Lead block in water basin
  A 1 m³ lead cube is located at the center of a water basin ($2.0\mathrm{m}\times 4.0\mathrm{m}\times 2.0\mathrm{m}$) positioned at the origin.
• Multiple lead blocks
  Four lead cubes with a side length of 0.6 m are located at four different positions in the container as listed in

  Table 2. The position of the four lead blocks.

  	$\mathit{x}-\mathbf{Position}$	$\mathit{y}-\mathbf{Position}$	$\mathit{z}-\mathbf{Position}$
  Block 1	$-0.5\mathrm{m}$	$-1.0\mathrm{m}$	$-0.5\mathrm{m}$
  Block 2	$-0.5\mathrm{m}$	$-1.0\mathrm{m}$	$0.5\mathrm{m}$
  Block 3	$-0.5\mathrm{m}$	$1.0\mathrm{m}$	$-0.5\mathrm{m}$
  Block 4	$0.5\mathrm{m}$	$-1.0\mathrm{m}$	$-0.5\mathrm{m}$

  . The location of the blocks is visualized in Figure 11.

Water and lead are chosen as the main materials, as they show a significant difference in their density and atomic number. Therefore, these two materials should be clearly distinguishable from another and are suitable for first test scenarios of the proposed reconstruction method. For each scenario 50 million air showers are generated, which is equivalent to an exposure time of around 17 min, comparable with other studies [45,46,47]. For the subtraction of the empty container, one billion air showers are simulated and afterwards scaled down to an equivalent of 50 million showers. The voxel size is set to 1 dm³ to allow for a sufficient spatial resolution within reasonable computation time and low level of background noise. The parameter sets are chosen manually for each scenario to create a clear visual discrimination between the different materials and objects while maintaining the shape and size of the items. If a measurement out of $M1$–$M9$ is by itself not capable of discriminating between the target objects and background noise, it is defined as negligible and therefore not utilized.

3.2. Separate Water and Lead Block

For scenario 1, two different reconstruction parameter sets are used. The set optimized for lead utilizes all measurements, except the photon production in the lower detector, since this is negligible (see Figure 8). The corresponding noise thresholds $t_{min}$ are listed in

Table 3. The parameter set with the noise thresholds $t_{min}$ optimized for lead.

	Photons	Neutrons	Electrons
Upper detector	20%	10%	–
Sidewise detectors	20%	10%	–
Lower detector—production	–	10%	–
Lower detector—absorption	40%	40%	40%

. All measurements except the neutron production in all detector layers are used for the optimized setup for water, with its noise thresholds shown in

Table 4. The parameter set with the noise thresholds $t_{min}$ optimized for water.

	Photons	Neutrons	Electrons
Upper detector	15%	–	–
Sidewise detectors	15%	–	–
Lower detector—production	15%	–	–
Lower detector—absorption	30%	30%	30%

.

The 3D voxel maps with the lead and water parameter sets are shown in Figure 12. The lead optimized map shows a clear picture of the lead block and no sign of the water block. When applying the water optimized set, the lead block starts to fade out and the water block in front becomes the prominent object in the map.

3.3. Lead Block in Water Basin

The set optimized for lead in Table 3 is not suitable to the more complex setup of the lead block in the water basin, since the secondary particles originating from the lead block also interact with the surrounding water resulting in more noise. Therefore, a higher set of thresholds (see

Table 5. The parameter set with the noise thresholds $t_{min}$ optimized for a lead block in water basin.

	Photons	Neutrons	Electrons
Upper detector	40%	30%	–
Sidewise detectors	40%	30%	–
Lower detector—production	–	–	–
Lower detector—absorption	60%	60%	60%

) is applied. The neutron production measurement from the lower detector is negligible in this scenario and therefore not utilized. The parameter set used for the reconstruction of the water basin is the same as in Table 4.

The 3D voxel map and its 2D projection in the $xy-\mathrm{plane}$ with the lead parameter set are shown in Figure 13. Both show the reconstructed lead block with additional noise around it. This is due to the additional interaction of the secondary particles originating the lead block with the surrounding water. It follows that the reconstruction of the realistic shape of a scanned object is more complicated with material in close proximity.

Figure 14 shows the 3D voxel map and its 2D projection in the $xy-\mathrm{plane}$ with the parameter set optimized for water. Both show the water basin being reconstructed, however a hole emerges at the position of the lead block. This relates to the usage of the photon production measurement with the lower detector. While materials with low density show a significant production of secondary photons in the lower detector, materials with high density lack this kind of signature because of the increased self-shielding effect. Therefore, the water parameter set includes this measurement (see Table 4), unlike the lead optimized set (see Table 5). However, the lead block within the water basin shields the lower detector against secondary photons produced within itself or its vicinity and therefore creates the gap seen in Figure 14.

3.4. Multiple Lead Blocks

Since, in this scenario, multiple lead blocks are shadowing each other in various directions, the lead set in Table 3 is adapted to create a higher $t_{min}$ set optimized for this scenario, which is shown in

Table 6. The parameter set with the noise thresholds $t_{min}$ optimized for multiple lead blocks.

	Photons	Neutrons	Electrons
Upper detector	30%	20%	–
Sidewise detectors	30%	20%	–
Lower detector—production	–	20%	–
Lower detector—absorption	50%	50%	50%

. This is necessary, since the amount of noise increases because of the additional interaction of the secondary particles with the different blocks.

The 3D voxel maps with the regular lead parameter set (see Table 3) and the parameter set optimized for multiple lead blocks (see Table 6) are shown in Figure 15. The map reconstructed with the regular lead parameter set contains all four blocks, though block 1 is depicted significantly smaller than the other objects. The reason for this feature is that block 2 shields the cube underneath from the incoming air shower particles and therefore reduces the probability of secondary particle production in block 1. In addition, artifacts related to noise are visible between the cubes. Those disturbances arise because of the interaction and scattering of secondary particles produced in block 2 with the blocks underneath and vice versa. If the higher noise thresholds according to Table 6 are applied, these artifacts are eliminated. However, block 1 also vanishes for this parameter set due to the lower number of secondary particles produced by it. A more complex combination of the different measurements can improve the 3D map, e.g., if block 1 is mainly reconstructed from the measurements in the lower detector plane.

4. Discussion

In this work, a novel method to analyze and reconstruct the content of a shipping container utilizing solely secondary particles in the context of cosmic ray tomography has been presented. It has been shown that secondary particles induced by air showers from cosmic rays carry information about the properties of materials, particularly the density and atomic number, in which they were produced. This information can be measured with a simplified geometry of four detector planes surrounding the shipping container and used for the reconstruction of the examined objects. Within this simplistic simulation model, a 3D reconstruction algorithm is introduced and validated by three different test scenarios.

As this is the first study of such kind, the aim of this work was to test the general hypothesis if secondary particles can be used independently of the muon measurements within a cosmic ray tomography application. Therefore, a basic and efficient design was chosen, which also allows an easy reproducibility of the results. Since the results successfully confirm the general principle that secondary particles may be used to derive information complementary to more typical approaches to cosmic ray tomography, further iterations of this work with various improvements will follow.

One emphasis will be the incorporation of a more realistic geometry setup. This includes a detailed model of a 20 ft ISO container and complex container content arrangements as seen by custom agencies around the world. Further, a realistic set of detectors will be simulated to measure photons, neutrons and electrons with viable spatial and energy resolutions, as well as reasonable acceptance and efficiencies according to the particle type and its energy. These points will undoubtedly reduce the sensitivity of the method presented in this work.

Additionally, a future focus will lie in the improvement of the reconstruction algorithm. The kinematic analysis shown in Section 2.2 reveals a variety of relations between the properties of secondary particles and the characteristics of the materials. The algorithm proposed in this work only utilizes a small subset of the available information. Therefore, a machine learning approach will be part of the future work, which includes all the kinematic variables to enhance the reconstruction and classification of the examined objects. With the output of this AI-classifier, an automatized material scan can be enabled, which will use the available measurements and additional parameters, like the noise thresholds, to asses the material characteristics and the shapes of the objects in the best possible way. The presented method using secondary particles is designed as a complementary measurement predominantly for the muon scattering approach to cosmic ray tomography. Hence, a combination of both results is planed to achieve unprecedented sensitivity for a cosmic ray tomography system for a shipping container.