Cosmic-Ray Tomography for Border Security

: A key task for customs workers is the interception of hazardous, illegal and counterfeit items in order to protect the health and safety of citizens. However, it is estimated that only a small fraction of cargo is inspected and an even smaller fraction of trafﬁcked goods are detected. Today, the most widely used technology for scanning vehicles, ranging from vans and trucks to railcars, is γ ray and X-ray radiography. New technologies are required to overcome current technological shortcomings, such as the inability to detect the target material composition, the usage of harmful ionising radiation sources and the resultant low throughput. Cosmic ray tomography (CRT) is a promising technology for cargo screening. Cosmic ray muons have average energies of around 10,000 times larger than a typical X-ray and therefore can penetrate relatively large and dense materials. By analysing muon scattering, it is possible to identify materials hidden inside shielding that is too thick or deep for other imaging methods. CRT is also completely passive, exploiting naturally occurring secondary cosmic radiation, and is therefore safe for humans and animals. Contrary to conventional X-ray-or γ -ray-based imaging techniques, CRT also allows material differentiation and anomaly localisation within the cargo or vehicle through the provision of 3D images. This article reviews the current state-of-the-art technology in CRT, critically assessing the strengths and weaknesses of the method, and suggesting further directions for development.


Introduction
The movement of illegal goods poses a significant threat to our society, demanding effective border controls. A large majority of the drugs trafficked between countries are synthetic. They are hard to regulate by authorities and have a significant negative impact on the economy. These drugs could also contain other substances that the buyer is unaware of, causing serious and sometimes deadly side effects. Firearms, ammunition and explosives are smuggled across borders, falling into the hands of drug cartels, terrorists and other transnational criminal organisations and individuals that utilize these weapons for criminal activity and acts of violence. There is also a crucial need to inspect cargo contents for special nuclear materials (SNM), radiological and biological substances and other materials that may be used for weapons of mass destruction. Table 1

Transnational Crime Estimated Annual Value (USD)
Drug Trafficking [1] 426-652 billion It is important to stress that the available data corresponds only to successful interceptions. Impressive figures have been communicated by the border control authorities. Around 1.1 million seizures of illicit drugs were reported in 2019 in Europe [5]. A large-scale joint action performed by EU law enforcement officers in the same year, for a duration of only 10 days, led to the identification of 476 potential victims of human trafficking (16 of which were minors), as well as the seizure of EUR 85 million worth of drugs and more than 1.2 tonnes of tobacco [6]. Data concerning drug seizures from these sources show that, although the majority of them involve small quantities, most of the total weight seized is attributed to a relatively small number of multi-kilogram consignments, underlying the fact that the detection of even just few illegal shipments can have a significant impact. However, based on a recent test conducted in the USA [7], it is estimated that only 5-10% of illegal goods are detected with conventional X-ray scanning technology.
Several techniques are employed to screen containers and other cargo with the aim of detecting illicit goods. Ion mobility spectrometry (IMS) and gas chromatography IMS (GC-IMS) are commonly used at border check-points together with X-ray scanners. IMS/GC-IMS devices are mainly used as desktop instruments for checking relatively small items and are not so efficient for large cargo inspection in lorries. Powerful X-ray technology (e.g., [8]) requires intensive high-energy radiation to penetrate through the walls of the cargo container and is subject to regulations on the use of radiation sources. Cargo inspection systems are based on linac accelerators to provide X-ray pulses with two energies (typically 4 and 6 MeV, or 6 and 9 MeV), alternating with a frequency of up to 400 Hz (see Figure 1). This approach may also be harmful for people potentially hidden in the cargo.
Millimetre wave technology can only detect objects hidden in or behind human clothes (standard airport full-body scanners) or other types of thin materials, but cannot penetrate through the metallic walls of containers. Figure 1. Conventional high energy X-ray cargo inspection system with material discrimination (similar to, e.g., [8]).
In recent years, to search for SNM, there has been a growth in the usage of so-called 'radiation portals' for non-intrusive cargo inspection, based on γ or neutron detectors [9]. These are passive devices that aim to detect the radiation emitted by SNM, based on generic assumptions about possible radiation sources [10]. Active devices are also in use, which produce intense pulsed emission of γ rays or neutrons to inspect the volume of interest. Active methods, however, pose complications in terms of radiological safety, including activation of the materials. The intrinsic weakness of any such methods, both passive and active, is that they can fail if the objects being searched are well shielded. This can be the case if the SNM is intentionally smuggled, for example, by enclosing it in a heavy metal casing or by filling the cargo with scrap metal. This led to the idea, first developed in a seminal paper by Borozdin et al. in 2003 [11], to use muons for non-intrusive inspection. This technique is usually referred to as muon tomography. In this paper, we will refer to this as cosmic ray tomography (CRT), to encompass more general variants of this method that also exploit other particles from cosmic ray showers as well as muon-induced secondary particles.
This article reviews the state-of-the-art in CRT for border security applications. Section 2 introduces the physical basis of the method. Section 3 reviews the existing muon portals. Section 4 reviews the relevant theoretical knowledge about the interactions of muons and other cosmogenic particles and the available simulation tools. Section 5 discusses the requirements for instrumentation and assesses the strengths and weaknesses of the main detection technologies. Section 6 is devoted to the data reconstruction and imaging tools. Section 7 elaborates on the problems to be solved for CRT to fully become a viable solution in security applications. The conclusions and prospects are summarised in Section 8.

Physical Principles of Muon Tomography
Muons are fundamental sub-atomic particles produced around 15-30 km above sea level in the atmosphere due to the decay of predominantly pions and kaons originating from the interactions of primary cosmic rays with nuclei in the upper atmosphere. At sea level, most of their energy spectrum is in the order of GeV. Muons belong to the leptons, a family of elementary particles that includes also electrons, τ-leptons and neutrinos, defined by being insensitive to the strong force. This implies that they do not undergo nuclear interactions, unlike protons, neutrons and other hadrons. Due to their high mass (about 200 times larger than the electron), at relevant energies (i.e., up to hundreds of GeV) muons do not suffer significantly from other types of energy losses, and can therefore penetrate through matter without losing much of their initial energy.
Their main energy loss mechanism is the ionisation of the atoms in traversed materials. At the typical momenta of cosmogenic muons (i.e., on the order of GeV), their mean energy loss rate (dE/dx) is close to the minimum of the Bethe function and is approximately proportional to the density of the material ρ [12]. Moreover, like any other charged particle, muons are deflected from their initial trajectory when interacting electromagnetically with the nuclei in materials, which offers a possibility to discriminate the elements by their different nuclear electric charge (or atomic number, Z).
These two physical processes, energy loss and deflection (illustrated in Figures 2 and 3), are the foundations of muon absorption radiography and muon scattering tomography, respectively. Muon absorption intrinsically yields 2D information in the form of a map of ρ projections. Therefore, a 3D map can only be obtained by combining different points of view. Scattering, instead, is intrinsically 3D. Muon tomography/radiography (sometimes also called muography) is increasingly being used for studying and monitoring geological and human-made structures, and for border control and security applications [13,14]. The former application relies on the absorption of cosmic ray muons passing through rock or other material, as the goal is typically to find cavities or volumes with different ρ, whereas the latter uses primarily muon scattering (deflection) to discriminate different materials. More details about the physical processes of interest for muon imaging are elaborated in Section 4.2.1.
This article describes border control and security applications and hence, from this point on, will focus mainly on muon scattering, as illustrated in Figure 3.

State-of-the-Art Cosmic Ray Tomography in Security Applications
According to published information, there are more than ten research groups developing CRT systems. Yet, significantly fewer initiatives have led so far to the development of CRT for the inspection of large cargo and trucks in security applications. In this section, we give just some representative examples, whose technical details are reviewed in other sections of this paper, in particular in Section 5.
Pioneering work in applying CRT to security applications for cargo inspection has been performed by the Los Alamos National Laboratory and its spin-off technology company, Decision Sciences [15][16][17]. They were the first group to demonstrate the effective detection of SNM, and, according to their information, they operate commercial portals in The Bahamas, Singapore and at the USA-Mexico border. Cosmic muons are tracked with detectors based on drift tubes, see Section 5.2.1. Their system is able to detect high-Z materials in a short time period, thus allowing the identification of SNM or the shielding of radioactive materials. Its operational capability for low Z-materials is described in Chapter 7 of ref. [14].
Drift tubes have also been chosen for the same purpose by the Israeli company Lingacom [18], and by an INFN team in Italy [19]. The Lingacom setup combines X-ray data (primary system) and CRT (secondary system) to reconstruct the three-dimensional scattering density map of the scanned volume with a proprietary algorithm [18]. The INFN project [19] reuses drift tubes originally built for the use in Compact Muon Solenoid (CMS) experiments at CERN. Their methodology includes a rough estimate of muon momentum based on the scattering induced by a slab of steel downstream of the volume of interest. We elaborate on this in Section 7.3.
Another seminal security-related project, started more than a decade ago, is the Muon Portal collaboration [20], which plans to build a large-area muon detector for non-invasive inspection of shipping containers in ports, searching for the presence of potential fissile (U and Pu) threats. Their detector is based on plastic scintillating strips, see Section 5.4.
Plastic scintillating strips are also used in the Cosmic Ray Inspection and Passive Tomography (CRIPT) project [21], which developed a large-size demonstrator also including a momentum estimation based on scattering in slabs of steel (Section 7.3). They also demonstrated the possibility to operate their muon scanner simultaneously with γ and neutron detectors, to combine CRT, γ ray spectroscopy and neutron counting in a joint data analysis [22].
A Spanish company, Muon Tomography Systems, declare on their home page (muon.systems/en) three main applications: civil engineering, industrial applications and security. Regarding security applications, unfortunately no publications are available. Though according to their home page, they have initiated the Muon Cargo project that aims to make CRT a true system for innocuous and non-invasive inspection of goods for the detection of threats and fraud at ports and borders. The project will be carried out in a port environment with financial support from the Spanish government agency Puertos del Estado. Regarding the detection systems used, some light could be shed from the two available publications [23,24] by Muon Tomography Systems about CRT. The publication describes a few industrial applications with multi-wire proportional chambers (MWPC), with a size of 1 × 1 m 2 ; the wires are spaced 4 mm apart and embedded in a mixture of argon-CO 2 gas.
One decade ago, the UK Atomic Weapons Establishment and Bristol University built demonstrators based on resistive plate chambers (RPC), see Section 5.2.3, aiming to search for special nuclear materials [25][26][27][28]. RPC detectors were also chosen for the Tsinghua University MUon Tomography facilitY (TUMUTY) [29,30], aiming to extend the application of the scattering method to the identification of illicit drugs and explosives, which are materials characterised by a low Z. We elaborate on this topic in Section 7.2. More recently, RPCs have been chosen as the detection technology for the TECNOMUSE project [31,32].

Muons at the Surface of the Earth
The muon flux at the Earth's surface is about 0.0167 cm −2 s −1 for horizontal detectors, with an average energy of 3-4 GeV for vertical muons and 6-7 GeV for muons coming from all directions. For GeV muons that dominate at the Earth's surface, the flux is the highest at the zenith, with an angular dependence of ∼cos 2 θ, where θ is the angle between the muon direction and the vertical.

Muon Interactions
Muons lose energy via excitation and ionisation (as shown in Figure 4), as well as via 'radiative processes', such as bremsstrahlung (as shown in Figure 5), electron-positron pair production and production of hadrons via inelastic scattering. The muon energy loss in radiative processes is proportional to the muon energy and dominates over the ionisation energy loss for muons above a few hundred GeV (depending on the material where muons propagate). However, for low-energy muons in the GeV range, such as those more abundant at the surface of the Earth, the ionisation process is dominant. More details on muon energy losses can be found in the Review of Particle Physics [12].

Multiple Coulomb Scattering
Multiple Coulomb scattering is the process that contributes the most to the muon deflection when a muon passes through different objects. Single scattering, also called Rutherford scattering, is described by the Rutherford formula. If the scattering amount is large, the combined effect of single scattering events (multiple scattering) can be approximated by Molière theory [33], with the scattering angle following a Gaussian distribution for small angles and the Rutherford formula at large angles. For small scattering angles, the scattering angle distribution in a plane perpendicular to particle motion is given by [12]: with where θ is the scattering angle in the plane perpendicular to the particle motion, p is the particle (muon) momentum in MeV, c is the speed of light (β = v/c), v is the particle velocity, z is the particle charge in units of the electron charge (z = 1 for muons), x is the distance travelled by the particle and X 0 is the radiation length of the material measured in the same units as x (either metres or g/cm 2 ). The formula for the parameter θ 0 (the characteristic scattering angle) is valid 10% for 10 −3 < x/X 0 < 10 2 . These conditions are satisfied for the application of CRT considered here. The characteristic muon scattering angle, θ 0 , is inversely proportional to the muon momentum, which makes low-energy muons (such as those abundant on the surface of the Earth) particularly attractive for measuring scattering angles in CRT applications. The accuracy of material identification may also improve significantly if the muon momentum is known, but this is quite difficult to achieve in practice, given the cost of such equipment. For accurate modelling of particle scattering, a detailed Lewis' theory [34] based on transport equations for charged particles is usually employed.

Muon Absorption and Negative Muon Capture
As muons pass through materials, they lose energy and can stop within the material if they have lost all kinetic energy. As the energy loss for GeV muons is dominated by ionisation, it is roughly proportional to the muon path measured in g/cm 2 (track length times density) and to the Z/A of the material. For ultra-relativistic muons, the dependence on the energy/momentum is relatively small. Observing the deficit of the muon flux from certain directions compared to others may be a sign of a high material density, given that for most materials, Z/A ≈ 0.5. This feature can complement the more conventional technique of measuring the scattering angle [35,36]. The main problem lies in the fact that muons are unstable particles with a mean lifetime of 2.2 µs. Some low-energy muons will decay on the flight-even in the air-contributing to muon disappearance.
For hydrogen, Z/A ≈ 1 and hydrogenous materials can, in principle, be discriminated by higher absorption pattern for a fixed muon path, measured in g/cm 2 (distance times density). However, without knowing a priori material density or muon path in the material in g/cm 2 , this feature is difficult to explore.
Upon stopping, negative muons will be captured into an atomic orbit, resulting in the formation of a muonic atom. As the muon mass is much higher than the electron mass, the muon orbit has a much smaller radius that overlaps with the nucleus, in particular, in high A materials. This may result in a muon being absorbed by a nucleus (a weak process) before it decays with an emission of a neutrino and a neutron, also accompanied by neutrons and γ rays from nuclear de-excitation. This process is sometimes called muon-induced fission or nucleus spallation. The probability of muon absorption by a nucleus increases with A. There is no simple approximation for the dependence of neutron multiplicity on the atomic number Z or mass A, but the general trend is that the number of emitted neutrons per muon absorption event also increases with A. An additional signature for high A materials is that the neutron and γ emission correlate in time with a disappearing muon.

Secondary Particles
In the context of CRT, the study of secondary particles remains relatively unexplored when compared with the vast body of research related to muon transmission and muon scattering-based tomographic measurements. Secondary particle analysis, in contrast to these more conventional techniques, deals with particles (predominantly photons, neutrons and electrons) produced as a result of cosmic ray interactions within the target material. The subsequent production and interaction of secondary particles are highly material specific, depending on both the atomic mass and density of the contents, which influences the ability of the secondary particles to escape, to be detected and therefore to be analysed. A representation of some processes contributing to secondary particle production is given in Figure 6. The study and analysis of secondary particles have the potential to significantly increase the ability of tomographs to discriminate between material composition, especially at low atomic masses. Naturally, it is not possible to identify all secondary particles and measure their properties using only detectors optimised for muon identification and tracking. Therefore, modified or additional detectors will be required to fully take advantage of this additional information, see Section 5.5. Current research relating to material imaging using secondary particles is limited to relatively few publications, some of which will be highlighted in the following paragraphs.
One of the first experimental verifications of the use of muon-induced electrons and photons (measured in coincidence with incoming muon tracks) for low-Z material imaging was in ref. [37], where, building on prior studies in ref. [38], a variety of small material samples (lead, copper and aluminium) were analysed. The results verified the feasibility of the technique, with the best results obtained for low-Z materials due to their improved transmissibility for γ rays. Using a similar approach, bones and soft tissues were imaged in ref. [39]. Following this, ref. [40] showed through simulation that materials can be discriminated by considering the ratio and number of induced γ rays and electrons, which was found to vary by more than two orders of magnitude depending on the material considered. Ref. [41] builds upon this work, proposing an updated detector configuration and a new three-dimensional reconstruction technique.
More recently, ref. [42] presented a method to combine information obtained from muon scattering tomography with coincident muon-induced photons and electrons to allow the analysis of low, medium and high density materials. Results from the simulation indicate an enhancement in material discrimination ability when utilising the proposed hybrid approach. These works were based on earlier investigations which laid the foundation for low atomic mass material differentiation using induced secondary particles. In the field of analytical chemistry, in the 1970s, it was proposed that muonic X-rays could be used to quantitatively analyse low-Z material samples [43,44]. Following this, ref. [45] demonstrated the use of X-ray fluorescence spectroscopy to perform non-destructive elemental analysis of extraterrestrial materials using a muon beam originating from a J-PARC MUSE pulsed muon source.
Although many of the prior works discussed have clear applications in the security domain, one of the first explicit applications of secondary particle analysis for this use case was in ref. [46]. It was proposed to use muon-induced neutrons and photons in coincidence with muon tracking information to analyse warheads for the presence of nuclear material. A similar methodology was used in ref. [47], whereby secondary neutrons produced in coincidence with incoming cosmic ray muon tracks (neutron tagging) were used to detect stand-alone nuclear materials. It was found that when imaging a cube of enriched uranium, neutron tagging successfully detected the object within one hour with a high accuracy. However, a precise determination of the object shape required tens of hours.
The use of secondary particles for the analysis of shipping containers was also addressed by refs. [35,36]. In this work, the use of muonic X-rays stopped muons and their subsequent capture, and material reconstruction was simulated using incident cosmic ray neutrons. High atomic mass shielding and nuclear materials were considered as targets. When considering muonic X-ray detection alone, material identification was found to be possible in ideal conditions. However, in a realistic scenario, this was found to be unfeasible. Several scenarios were then considered, whereby muon disappearance was observed in conjunction with various secondary particles. Stand-alone muon disappearance was shown to be sensitive to material density in realistic conditions, and the sensitivity was reduced when requiring secondary particle coincidence. When considering the detection of excess neutrons and γ rays produced as a result of cosmic ray neutron interactions, it was again found to be feasible only under ideal conditions. Given realistic detector scenarios, the technique struggles to detect shielded nuclear materials. Though the results of the studies performed in refs. [35,36] seem initially pessimistic, recent and potentially further advances in detector technology and analysis algorithms may help to make this approach feasible.
Most recently, in ref. [48], a novel approach was proposed to utilise stand-alone information from secondary particles to reconstruct materials in three dimensions and subsequently discriminate them for the purpose of shipping container analysis. The results from the simulations proved in principle that secondary particles carry information relating to target material properties, particularly density and atomic number. Only ideal detectors and acceptances, however, were considered for this study. Nevertheless, the analysis of this information provides a method for material discrimination which is complementary in nature to muon transmission or scattering tomographic techniques, and which can be used to enhance the sensitivity of these tomographic approaches, particularly when the results are combined in a hybrid approach.
Complementary to the body of work related to passively induced secondary particle production, studies have been performed analysing γ rays and neutrons as a result of active neutron-and γ-ray-based inspection of sea containers and materials, for example, in [49][50][51][52]. A full discussion of these works is not within the scope of this article; however, many techniques and results overlap with the topics discussed within this section.

Secondary Particle Interactions
The following sections will describe the underlying physical processes which give rise to secondary particles, as well as their resulting kinematic properties. This is limited to a discussion of photon, neutron and electron production, since other secondary particles are produced at a rate which is not statistically significant for material identification.

Muon Interactions
Muon interactions are described fully in Sections 4.2.1-4.2.3, where the most interesting process in terms of secondary particle production is negative muon capture and absorption.
Further particle production may occur as a result of interactions between particles emitted during these processes and the target materials. Photo-nuclear or neutron-nuclear processes may result in the production of additional neutrons and γ rays. In addition, nuclear materials may be caused to fission as a result of neutron absorption, again resulting in additional neutrons and γ rays.

Neutron Interactions
Cosmic ray neutrons, also produced during air showers, interact alongside cosmic ray muons within target materials. Neutrons are uncharged particles; therefore, they do not interact electromagnetically, meaning they have a longer path length in a given material. There are two main categories of neutron interactions which are governed by a strong force: elastic or inelastic scattering and absorption.
Elastic scattering occurs when a neutron transfers kinetic energy to a nucleus, with the total energy of the system remaining constant. The energy loss of the neutron is given by 2EA/(A + 1) 2 , where E is the kinetic energy of the neutron and A is the atomic number of the nucleus. From this, it is evident that a larger amount of energy is lost for materials with smaller atomic numbers.
More interesting in terms of secondary particle production is the process of inelastic scattering, which results in the non-conservation of the energy of the system. In this case, some energy is used to elevate the nucleus to an excited state or to produce new particles. Upon relaxation of the nucleus to the ground state, element specific γ rays are emitted. Inelastic scattering occurs above a specific threshold energy, which is given by where A is the atomic number of the nucleus, E t is the inelastic threshold energy and 1 is the energy of the first excited state of the target nucleus. The cross-section for inelastic scattering is relatively low for small and light nuclei, with t generally decreasing with an increasing atomic number. Inelastic scattering therefore plays an important role at high neutron energy and in heavy elements.
Neutron absorption is most important in thermal energy, with the cross-section falling off with increasing neutron energy. Resonant peaks, however, feature in the cross-section at higher energies, occurring when the neutron energy is close to that of the energy of the first excited state of the nucleus. When neutron absorption occurs, there is a probability of neutron-induced fission (for radioactive materials) or spallation of the target material, resulting in the emission of excess neutrons. In addition, radiative capture may occur, resulting in the emission of neutrons and γ rays.

Additional Interactions
In addition to the dominant muon and neutron interactions, other air shower particles, such as electrons, protons and photons, play a minor role in the production of secondary particles. The impact of such interactions has not been well studied in the context of cosmic ray tomography, and is not expected to be a statistically significant source of secondary particle production. Nevertheless, the dominant secondary particle production methods are briefly described in the paragraph below for each case.
Protons interact via multiple scattering, ionisation or elastic/inelastic scattering, producing electrons, γ rays and hadrons. Very rarely, antiprotons may be captured at rest, generating electrons and causing nuclear decay, resulting in the production of fission products and γ rays. Electrons and positrons interact through multiple scattering, ionisation, bremsstrahlung and electron-nuclear processes, again producing secondary electrons, γ rays and hadrons. In addition, electron-positron annihilation results in the production of two 511 keV γ rays. Photons produce electron-positron pairs through conversion processes. Compton scattering processes and the photoelectric effect may also produce recoil electrons or nuclear excitation. Photo-nuclear processes may result in nuclear decay along with hadron or alpha particle emission.

Simulation Codes
To effectively reconstruct the image of a scanned object and identify its composition, measurements should be compared to a reliable model. Simulation codes include muon generators, muon propagation codes and detector response models. Different simulation steps and their interdependence are sketched in Figure 7 and described in the next subsections.

Muon Generators
A number of codes exist to sample cosmic ray muons according to their energy spectra and angular distributions at sea level, at higher altitudes or deep underground. For the goals of CRT, we will mention here only a few codes that sample muons at the sea level. These codes can be split into two main categories: (1) simulations of whole atmospheric showers with the production and propagation of all particles and (2) muon generators based on parameterisations of or fits to measured muon energy spectra and angular distributions.
The most widely used code belonging to the first class is CORSIKA [53]. This uses parameterisations of the measured primary cosmic ray spectra of different species and accurate models for particle interactions at high energies tuned to the measurements to ensure high precision of calculated particle spectra at sea level. Similar simulations of the development of atmospheric showers can also be performed with multi-particle transport codes, such as FLUKA [54,55] and GEANT4 [56,57], but require user input to choose models for physical processes and atmospheric properties, e.g., varying density. The main advantage of this approach is to record all particles that reach sea level, not only muons, so they can be later used to simulate effects in the cargo. The main disadvantage is a large processing time due to propagation of many particles through the whole atmosphere.
A widely known code CRY [58] uses libraries of particles previously generated with the code MCNPX [59] through the simulation of atmospheric showers. The libraries include all particles and the code is fast. Its disadvantages include underestimation of particle fluxes, since only showers induced by primary protons are simulated, and a large bin width in energy and zenith angle which affect the accuracy of the results.
Simple muon generators based on parameterisation of muon fluxes tuned to the measurements are available and widely used for various purposes. The well-known Gaisser's parameterisation of the muon flux [60] was modified by a number of users to include low energy muons by adding muon energy losses in the atmosphere and the probability of muon decay, and extended to large zenith angles by taking into account the curvature of the Earth and it atmosphere. One such parameterisation was used previously for some CRT projects, underground particle physics experiments, such as DUNE and LZ, and for other muon simulations (see, for instance, [35,[61][62][63]). A similar code, MUSIBO, is now included in the SilentBorder [64] simulation toolkit. Figure 8 shows the muon energy spectrum (left) and the muon zenith angle distributions (right) as generated by MUSIBO on a surface of a box surrounding the inspected volume and a set of muon detectors. Other parameterisations of the muon flux include those described in refs. [65,66]. More Monte Carlo generators using a similar approach of parameterising fits to the measurements are available, i.e., EcoMug [67], CMSCGEN [68] and others. Some codes allow muon sampling on different surfaces (e.g., a box or a sphere) whereas others, such as CRY, are limited to sampling on a horizontal surface. Hence, for the case of muon panels positioned horizontally and vertically to cover a larger solid angle for muon detection, the application of CRY may require a larger computation time to generate and transport a specific muon sample through the setup.

Particle Transport and Detection
Muon transport can be carried out using multi-purpose toolkits, such as GEANT4 [56,57] and FLUKA [54,55], or with specially designed, house-built codes such as MUSIC [69,70] or PROPOSAL [71]. For muon scattering tomography, when the thickness of material for muon transport is small and the CPU usage is not an issue, multi-purpose toolkits have a few key advantages over specially designed codes. One reason is that GEANT4 (as an example) can be used in all-particle transport and in the simulation for the detector response. Furthermore, these codes can propagate all particles, not just muons or charged leptons in a more general case, to evaluate their effects in detectors by simulating, for instance, light collection. On short distances, usually assumed for muon-scattering tomography at sea level, the codes, such as GEANT4, have a clear advantage of being sufficiently accurate and fast to fulfil their goals.
Simulation of muons and detection of other particles can be carried out in many applications with the GEANT4 toolkit [56,57]. GEANT4 can handle scintillation and Cherenkov photon transport but cannot transport relatively slow drifting electrons in gas tubes, for which a specially designed tool can be used and added to the GEANT4 framework. Other light transport codes exist but will not be discussed here.

Requirements for Detectors
While X-ray scanners use a controlled intensity beam in the form of a fan, CRT uses natural cosmic radiation with a flux intensity of 10,000 muons/min/m 2 . In order to scan large objects such as shipping containers or trucks, CRT detectors must have a large area to register the maximum number of muons that pass through the object under study. This gives rise to the main requirement for the detection technology-the ability to create tracking detectors with dimensions ≈ 3 × 6 m. Since in a commercial scenario, the required number of such muon scanners is hundreds or even thousands of units, the life cycle cost of the scanner and its maintenance should be comparable to the cost of X-ray scanners. The spatial resolution required to provide a high detail of the contents of a cargo container is determined by the accuracy required to detect a small change in the angle between the trajectories of the incoming and outgoing muons. The difference in scattering angles between medium and high Z materials is on the order of ∼ten milliradians, and between medium and low Z is on the order of one milliradian.
The scanning time required to detect contraband should not exceed a few minutes to ensure that a large number of consignments are scanned. Since the CRT must work almost non-stop, it is necessary to guarantee the stability of the system for a long time. To scan sea containers or trucks at checkpoints, a muon tomograph must work reliably in conditions close to field conditions, resist weather changes and ensure reliable operation in a wide range of ambient temperatures. Due to field operation conditions, the system should have a preferably robust design and a long-term performance stability; gas handling or high voltage systems are not the preferred choice. Due to the requirement of a high tracking accuracy, alignment calibration procedures should be built in.
Various tracking detectors are applied for cosmic muon tomographies, such as scintillators, gaseous detectors or solid-state detectors.

Drift Tubes
Drift tubes are gas-based detectors, usually aluminium, with a thin central anode wire running through the centre. The incoming muon ionises the gas, creating an electron-ion pair, and under the action of the applied voltage, the electron avalanche drifts to the anode wire. This drift time is measured and used to determine the position information. The ability of drift tubes to measure coordinates along the trajectory of charged particles determines their use as muon detectors. One of the methods used by drift tubes to determine positions is to measure the drift time of electrons resulting from ionisation, which is the result of the collisions of muons with gas particles. Drift time is defined as the time required for the ionising charge (electrons) to drift from their origin to the anode wire, as shown in the schematic drawing in Figure 9. The electron drift time can be converted into a drift radius, r, which is the distance between the anode wire and the muon track. The muon track can be reconstructed by using at least three drift tube layers for each spatial coordinate. Drift tube planes have been selected by Los Alamos National Laboratory (LANL) as shown in Figure 10 for use in the Large Muon Tracker [15]. The drift tube detectors are made of aluminium and are filled with a 47.5% argon, 42.5% tetrafluoromethane, 7.5% ethane and 2.5% helium gas mixture. They have obtained spatial and angular resolutions of 0.4 mm and 2 mrad (FWHM) with a muon tracking efficiency of close to 100%.  The research on CRT at LANL was started in 2001 to address the problem of nuclear smuggling in a way that would not expose people in a vehicle to additional doses of radiation and maximise the use of safe cosmic rays given their high muon penetrating power.

Barrel Chambers (Drift Tubes)
The INFN collaboration has constructed a Muon Scattering Tomography (MST) demonstrator using the detector barrel originally produced for the Compact Muon Solenoid (CMS) experiment at the Large Hadron Collider. The detector consists of two 300 cm × 250 cm drift chambers (orientated horizontally), vertically separated by a distance of about 160 cm, as shown in Figure 11. A gas mixture of 85% argon and 15% carbon dioxide is used. The chambers enclose a volume of more than 11 m 3 ; one of the largest volumes for MST tests in the world. The scattering angle projection on the x-y plane (with y-axis vertical) is referred to as the ∆φ angle, and the projection on the z-y plane denotes the ∆θ angle. The obtained track spatial and angular resolutions for high energy muons are about 300 µm for single-point position resolution and approximately 1.2 mrad and 10 mrad for ∆φ and ∆θ, respectively. Experimental validations of MST have proven the capability of the system to discriminate different materials [73]. Data were recorded for six blocks of different materials placed between the chambers. Each block is a cube with sides of a length of 10 cm, except tungsten, which is 8 cm high. The measured scattering densities for tungsten are shown in Figure 11 (right) using simulated and measured data as a function of the nominal scattering density. The reconstructed densities do not follow a linear trend but clearly show a saturation effect. In particular, the measured scattering density of tungsten is not much larger than that of lead.

Resistive Plate Chambers
The University of Bristol, in collaboration with the Atomic Weapons Establishment (AWE), have developed a prototype CRT based on Resistive Plate Chambers (RPCs) and evaluated its performance for nuclear safeguards. In the Bristol demonstrator [28], the RPC detectors were composed of six x-y layers, three above and three below the testing volume, with each layer containing two RPC detectors with orthogonal strip orientations. The gas mixture used was composed of 60% argon, 30% freon gas R1 34A and 10% iso-butane. Each RPC has a position resolution of 0.5 mm and an average detection efficiency of 95%. Experimental results were obtained from scanning 50 mm × 50 mm × 50 mm blocks of aluminium, iron and tungsten. The results show that the detectors work as expected and prove the technique's potential. However, the scanning time required to obtain good images amounts to several hours, due to the relatively slow time response of the RPCs and requires further developments to be practical for homeland security applications.

The Multigap Resistive Plate Chambers
Since 2011, the cosmic muon tomography group at Tsinghua University has been developing a prototype of the Tsinghua University cosmic ray MUon Tomography facilitY (TUMUTY) [29]. There are six layers of 73.6 cm × 73.6 cm large scale 2D position sensitive Multigap Resistive Plate Chambers (MRPC) detectors in TUMUTY, with 224 channels of induced signals on each dimension, and a total of 2688 channels for the whole facility. A gas mixture which contains 90% freon gas Rl 34A, 5% iso-butane and 5% sulfur hexafluoride is used as the working gas. From the inside out, it is composed of 5 slices of inner glass with 0.55 mm thickness, 2 slices of outer glass with 0.7 mm thickness, carbon films, Mylar films, PCB plates and honeycomb plates. The six layers of the gas gaps divided by 5 slices of the inner glass are 0.25 mm thick each. Carbon films adhered to the outer glass are used as highvoltage electrodes, while Mylar films act as insulators between the high-voltage electrodes and the PCB plates. The charge-induced signals are acquired from the Cu readout strips attached to the PCB plates.
An experimental validation of the material discrimination ability of muon scattering tomography at the TUMUTY facility was performed with flour, aluminium, steel and lead, chosen because they represent materials with atomic numbers ranging from low to high. A supervised machine learning algorithm called a support vector machine (SVM) uses a set of training data to build a model that can predict the class label of new unseen data. It is commonly used for classification problems. When using SVM classifiers, the overall classification accuracies of these four materials at 1, 5, and 10 minute measurements are 70, 95, and 99%, respectively.

Gas Electron Multiplier Detectors
The Melbourne Muon Tomography System is based on gas electron multiplier (GEM) detectors, used to measure multiple scattering of atmospheric cosmic ray muons in matter for the detection of heavily shielded high-Z radioactive materials (U and Pu) in cargo or vehicles [74]. The GEM detector is a gas detector with a microstructure for detecting charged particles. It uses a thin Kapton foil coated with copper layers on both sides, into which a lattice of chemically etched holes is punched. A voltage is applied to the GEM foil and the resulting strong electric field at the holes can cause an avalanche of ions and electrons to flow through the holes. The electrons are collected by a suitable device, in this case, a reading plane with x-y stripes. A prototype CRT system using four GEM detectors and the corresponding electronics to read 768 channels was used to test the applicability of GEM detectors in CRT. They were operated using 70% argon and 30% carbon dioxide as the gas mixture. Using several thousand cosmic ray muons recorded by the station, a clear separation of samples with medium and high Z (Fe, Pb and Ta) and restoration of their shape was shown. This demonstrates that GEM-based muon tomography is possible in principle. Unfortunately, this method was not further developed due to the complexity of this technology for the production of large-scale detectors.

Solid State Detectors
Silicon-based tracking detectors are used at most large collider experiments and provide precise charged particle tracking with a spatial resolution better than thirty microns. However, the disadvantage of such a system is that individual sensors usually have an active area of several tens of cm 2 that require arranging. They can be arranged to cover tens of square meters in cylindrical or flat panel geometries; however, the channel count and subsequent cost of readout electronics are significant for large applications. In addition, the high-density electronics also require active liquid cooling in most cases, and transportation of a cooling system along with the detector in the field may not be feasible. Despite these drawbacks, there is a CRT prototype based on silicon strip sensors from the ATLAS semiconductor tracker [75]. Moreover, in experimental tests, objects made of different materials were reconstructed, proving the viability of this technology (see Figure 12).
Ref. [76] argues that silicon microstrip detectors are competitive with drift tubes once the maintenance costs are factored in and that its superior position resolution makes it an attractive choice. Moreover, silicon detectors are much less sensitive than gaseous ones to environmental conditions.

Scintillators
Scintillators are another type of detector that are widely applied to CRT. Scintillation detectors have been commercially produced and used for several decades. This is due to the fact that they are robust, easy to build, can be easily optimised through engineering and machinery, are relatively cheap, have little sensitivity to environmental conditions and provide a detection efficiency close to 100%.
The amount of light emitted is proportional to the energy deposited by the muon. The photodetector measures the number of photons produced inside the scintillator. A photomultiplier tube (PMT) or a silicon photomultiplier (SiPM) are used as photodetectors. Scintillation detectors are made from organic or inorganic materials in a liquid, gaseous or solid state. For example, CRT uses only organic solid-state scintillators made of plastic, which is easy for machinery and allows the creation of detectors of various geometries. Scintillation plates can consist of large areas, up to 2 m long, and are commercially produced. The method of direct light collection from the surface of the plate using a PMT or SiPM, however, does not allow one to obtain a spatial resolution of the muon track better than ∼1 cm. To address this issue, light from large-area bulk scintillation plates is collected not directly by a photodetector but indirectly by means of wavelength-shifting fibres (WLS). WLS fibres absorb primary light from the scintillation plate and re-emit photons at higher wavelengths. WLS is applied in the form of a grid on both sides of the scintillation plate, providing the positional sensitivity of the detector. This method also makes it possible to significantly reduce the number of expensive photodetectors and electronic channels. Another way of achieving a high spatial resolution for muon track is using segmented detectors. For this reason, the high-energy physics community has developed detectors based on strips of plastic scintillators of a triangular or square shape. The use of plastic scintillating fibres up to several meters long assembled in multi-row mats has also become widespread. This technology makes it possible to ensure the accuracy of the determination of the hit position resolution for muons at the level of 150 µm [77]. The UT and LT detection modules are separated by a vertical distance of 1 m and the SPEC module is separated by 0.5 m. This makes it possible to achieve angular resolutions of ∼7 mrad, with a high true positive fraction of >90% and a low false positive fraction <10% [21]. CRIPT is a system designed to detect smuggled nuclear materials such as uranium or plutonium. Each super-layer uses plastic scintillating bars of triangular shape for muon detection. The bars are combined into layers that are positioned orthogonally with respect to one another. The energy deposited by a muon in the top and adjacent bottom bars depends on the length of the muon's path through the bars. This makes it possible to reconstruct the muon hit-position using the weighted average of the strip positions, with the weight given by the energy deposited. The CRIPT spectrometer is designed to provide an estimate of the momenta of atmospheric muons, as elaborated upon in Section 7.3.

Plastic Scintillating Strips: Square Bars
The Muon Portal Project was started in 2012 by a collaboration of research centres and industrial partners [78]. Their aim is the production of a full-scale prototype for the scanning of travelling containers. The experimental setup is based on four XY detector planes, each providing the X and Y position measurements; two placed below and two placed above the volume to be inspected. The overall size of the detector fits that of a real twenty-foot equivalent unit container, which has a size of 6 m × 3 m × 3 m. For a suitable implementation of the detection setup, each plane consists of six modules (1 m × 3 m each) in a proper geometry to cover both the xand y-coordinates by the same type of modules without leaving any dead area. Each detection module is segmented into 100 strips of extruded plastic scintillators (1 × 1 × 300 cm 3 ), with WLS fibres to transport the light produced in the scintillator to the silicon photo multiplier (SiPM) placed at one of the fibre ends. The spatial resolution, on the order of a few mm, will be suitable to provide a good tracking capability for each muon, allowing the reconstruction of the incoming and outgoing tracks and consequently the scattering angle with a geometrical angular resolution of about 3 mrad. After the construction of the Muon Portal detector was completed, the first set of measurements was made with a block of lead a few dm thick. Preliminary results have shown that the angular resolution and muon tracking procedures are adequate. One of the main problems is still related to the overall detection efficiency.
Currently, it takes too much time (4-8 h) to collect sufficient statistics (50-100 thousand tracks) to obtain an acceptable image of a lead block of several dm due to the small global efficiency which is only (0.9) 16 × 100% = 18%. It is planned to increase the global efficiency up to 44%.

Scintillating Fibres
Recently, a novel scintillating fibre tracker for the CRT of low Z materials has been developed in a collaboration of Estonian institutions: GScan OU and University of Tartu, Estonia [77].
A prototype scintillating fibre detector system has been developed using fibres of 1 mm in diameter arranged in a double layer array forming a detector plate (see Figure 13). The active area of each detector plate is 247 × 247 mm 2 . The detector plates in the upper hodoscope are placed with 75 mm pitch, the distance between lower plate of the upper hodoscope and the lower hodoscope is 250 mm and thus the dimensions of the VOI are 247 × 247 × 250 mm 3 . The data acquisition system of the prototype implements eight CAEN DT5550W boards to readout Ketek (PA3325-WB-0808) and Hamamatsu (S13361-3050AE-08) SiPMs arrays to collect scintillation light from the fibres. A particle tracking algorithm determines the particle hit location in a detector plate and its trajectory in the hodoscope provides an accuracy of 120 µm per detector plate. Experimental tests were performed with objects made of Plexiglas, aluminium and iron with a size of 6 × 6 × 6 cm 3 and a measurement time of 30 min. Recorded data were analyzed using particle track filtering (PTF) and multi-modality tomographic reconstruction (MMTR) methods. All objects were successfully detected, which proves the viability of the developed detector concept.

Plastic Scintillator with WLS-Fibres Readout
A feasibility study of a prototype muon tomography system based on a plastic scintillator plates with scintillation light readouts using WLS fibres was developed in a collaboration of Korean institutions [79]. The muon detection part consists of four detector modules with a top and bottom configuration. The incoming and outgoing muon trajectories are measured using two detectors above and below the target object, respectively. Each muon detector consists of a plastic scintillator plate (50 × 50 cm 2 ), 64 WLS fibres attached to the top and bottom sides of the plastic scintillator, and silicon photomultipliers (SiPMs) (0.3 × 0.3 cm 2 ) connected to both ends of each WLS fibre. Figure 14 provides a schematic view of a plastic scintillator detector with a WLS fibre readout. The developed muon detector can detect the X-and Y-positions simultaneously using a position determination algorithm. The X-and Y-positions are determined by analysing the distribution of scintillation photons between WLS fibres on the top and bottom sides of a scintillator plate using Anger logic [80].  Table 2 lists various institutions involved in the development of cosmic-ray muon tomography and detection methods.

Detectors for Secondary Particles
A wide range of instrumentation is suitable for the detection of secondary particles (Section 4.2.4). This section will present the most relevant detector types considered so far and discuss other potential options, taking into account features required for security applications. Since secondary particle analysis is a relatively new concept, no stringent requirements have yet been defined regarding detector design parameters. The energies of secondary particles produced within a material are generally centred around 1 MeV; however, if absorption features of cosmic ray photons and neutrons are considered, then the upper limit on the energy range should be extended to at least around 100 MeV [48]. Studies have not yet been performed to determine the spacial resolution required for indepth secondary particle analysis, but given the difficulty in locating the origin or the secondary particle emission it is not expected that state-of-the-art resolutions are necessary.
In refs. [37,38] secondary photon detection was achieved using high purity germanium detectors (HPGe), similar to earlier works described in refs. [43][44][45]. The advantages of HPGe detectors for γ ray detection in particular are the excellent energy resolution (about 1% [81]) and high radio purity that leads to high sensitivity in detecting γ lines. The main issues include the relatively small detector volume, poor operation at low temperatures (about 100 • K and high cost. In their follow up work in ref. [39] the HPGe was replaced by a more versatile, robust and cost effective plastic scintillation detector.
Four plastic scintillator plates of 0.5 × 0.5 × 0.05 m 3 were arranged in a box geometry around the target object. With this setup, the achieved timing resolution was around 2.5 times shorter than with the HPGe setup. This approach was also considered in refs. [40,41]. In ref. [42], a similar box geometry was also considered, with the EJ200 plastic scintillator proposed as the detector material. No large-scale experiments have been performed which deal with the detection of muon-induced γ rays and electrons in this context.
Ref. [82] studied a variety of options for the detection of muon-induced fission neutrons for high atomic mass materials, concluding that for their application the EJ301 liquid scintillator yielded optimal results. Other detection materials considered were polyethylenemoderated 3 He-based drift tubes and plastic scintillators. In ref. [47], two EJ301 liquid scintillator detectors were used. This type of scintillator, as well some others, provides a powerful method to separate γ-ray-and neutron-induced events using the pulse shape discrimination (PSD) technique. The detectors utilised are small, just 12.5 cm in diameter and 5 cm thick, providing only limited coverage for secondary particle detection. The use of plastic scintillators sandwiched between two drift tube layers to reconstruct muon tracks for neutron and γ ray detection has been proposed in ref. [46].
Large scale monitors of radiation at borders, such as radiation portal monitors (RPMs), generally utilise panels of polyvinyl toluene plastic scintillators for γ ray detection and 3 He-based neutron detector tubes [83]. Generally, such monitors are utilised to detect only elevated levels of photons and neutrons compared with the background flux and therefore do not necessarily provide any positional resolution, which is required for a more detailed secondary particle analysis to be made.
As is evident from existing experimental setups, plastic scintillation detectors are widely used for all manners of secondary particle detection. This is due to their flexibility, robustness, high efficiency and reasonable price. Sufficient spacial and angular resolution can be also achieved with relative ease using a detector composed of bars or slats. In addition, the construction of a large area detector is fairly simple, which is not necessarily the case when considering an alternative such as liquid scintillator detectors. The disadvantages of plastic and most other scintillators include relatively poor energy resolution for γ line identification and the large volumes required for high-efficiency particle detection which affects the cost of the setup.
Thallium-doped sodium iodide (NaI(Tl)) detectors are also widely used for photon detection, providing improved energy resolution compared with plastic scintillators (5-10% depending on γ ray energy for NaI [81] compared with up to 12% for 0.662 MeV γ rays for plastic scintillators [84]), as well as an improved efficiency and lower price when compared with HPGe detectors. Good energy resolution in both HPGe and NaI detectors can allow for spectroscopic identification of certain elements under suitable conditions. The cost of HPGe and NaI(Tl) detectors is significantly higher than for plastic scintillators and large-area coverage may be a challenge given the scale of many security-related tasks [81,85].
Zinc sulphide (ZnS) scintillator detectors are also widely used for neutron detection. ZnS combined with powdered boron nitrite (BN) or lithium flouride (LiF) have been considered as a replacement for 3 He-based neutron detectors in RPMs due to a global shortage of 3 He [86]. Boron-based powders have been shown to provide less scintillation light output, and they can be manufactured with a significantly reduced cost compared with LiF and have a performance similar to 3 He. Boron-(B), lithium-(Li) or gadolinium (Gd)loaded plastic and liquid scintillators have also been considered for neutron detection, for example, in refs. [87][88][89][90], providing an often comparable or improved detection efficiency to 3 He-based neutron detectors but with a lower γ ray rejection level.
When considering stand-alone secondary particle reconstruction, as discussed in ref. [48], in order to utilise the full range of information from secondary particles, two energy regimes must be considered. Widely available plastic scintillators are well suited for detecting secondary particles with energies of around 1 MeV. Higher energy particles, in particular neutrons, require a moderator layer to first slow them down or thermalise them. This layer can also act as a conversion medium for high energy γ rays. It could be possible to benefit from a combination of a moderation/conversion material with secondary particle detectors, which could be placed in alternating layers between muon detectors to provide a crude measurement of muon momentum (as discussed later in Section 7.3), in addition to secondary particle detection and tracking.

Data Reconstruction from Simulated and Real Events
MST is a method that uses information on muon trajectories through the VOI in order to create a 3D scattering density map. The MST working principle relies primarily on measuring the deflection of muons penetrating the object to be imaged. Since the true path of muons traversing the VOI cannot be known, reconstruction algorithms make assumptions and model muon paths depending on the model through the VOI and associate muon deflection or passage within a certain region of space in VOI. Predicting the path is crucial for accurate reconstruction algorithms as they affect the image quality dependent on the reconstruction algorithm. Path estimation algorithms usually require the initial and final position and direction of the particles. Several path models are available and will be discussed below, the most common path models are presented in Figure 15 and described in Section 6.1. A common way of performing the reconstruction is by dividing the VOI into 3D elementary volumes called voxels, usually on the order of a cm. For each muon, the reconstruction algorithm will identify which voxels are most likely compatible with the true muon path. Once the traversed voxels are identified, they receive a score reflecting their scattering or density properties based on the scattering angle, localisation, displacement and momentum measurements, if available. Several of these algorithms have been developed or borrowed from conventional CT tomography with different levels of complexity and realism, each using different assumptions. A 3D scattering density map can then be used as an input for various image processing algorithms such as clustering, filtering, feature detection or material identification. After processing every muon event, one obtains a scattering density map of the VOI, as shown in Figure 16, which can then be used as input for more advanced clustering and material identification algorithms based on machine learning techniques. Such post reconstruction algorithms are described in Section 6.2.

Path Models and Reconstruction Algorithms
The most widely used reconstruction algorithms are briefly summarised in this section with reference to those used in CRT. Reconstruction algorithms are generally implemented by introducing concept of voxels, i.e., the VOI is divided into elementary 3D volumes called voxels. The voxels are typically on the order of cms. However, the size can vary and depends on the particle flux density, required resolution and computation time available and is usually subject to optimisation. Each voxel represents a value consisting of a path length or a total scattering angle. However, there are scenarios where the threat object is comparable with the size of the grid or voxel and is not localised in the voxel; then, it is necessary to reconstruct the tracks and visualise the contents of the container independent of the grid and 3D voxels. The presence of a three-dimensional grid also limits the process of automatic identification of objects. Developed methods based on the analysis point data help to detect the presence of high-Z materials inside the container [91].

Point of Closest Approach
PoCA is a commonly used path model and a reconstruction algorithm in MST. As its name indicates, this approach assumes that muon scattering occurs at a single point in space which is located at the closest point between the incoming and outgoing muon trajectories. In 2D space, this point would be the intersection between two straight lines. However, in 3D space, two nonparallel lines might not intersect and the closest point is located in the middle of the (only) segment perpendicular to both lines. Its position can be computed from the parametric equations of the incoming and outgoing tracks solving the following equation system [92]: where P = p + td is the parametric equation of tracks and p 1 , d 1 and p 2 , d 2 are the known points and directions of incoming and outgoing tracks, respectively. The score s reflects the scattering power and can be associated with the voxel where the PoCA is located. Moreover, if a muon momentum measurement, p, is available, s ∝ p θ and, if not, s ∝ θ, where θ is the measured scattering angle between incoming and outgoing tracks. A first noise reduction step can be implemented within the algorithm by applying cuts on θ and rejecting events with too high or too low scattering angles. Once every muon event has been processed, one can use statistical estimators, such as the variance or the median, to compute the total scattering power of a voxel and build the 3D scattering density map of the volume, as shown in Figure 16. The PoCA algorithm is very popular because of its simplicity; however, it has the fundamental shortcoming of assuming a single hard scattering per event, which is often too unrealistic. Besides, the PoCA point of an event with almost parallel incoming and outgoing tracks can be reconstructed outside the VOI, as shown in Figure 16 on the left.
Such events have to be discarded. In cases of non perfect spatial resolution, up to 40% of muon PoCA points can be reconstructed outside of the VOI.

Straight Line Path
SLP is a classical approach utilising the entry and exit information of each muon track pair penetrating the VOI. This model assumes that the muon trajectory follows a straight line trajectory or that scattering is negligible. As a result, SLP is the easiest and computationally the fastest among the path estimation algorithms [93]. In voxelized space, SLP can be efficiently calculated using the Siddon algorithm [94]. This algorithm is reasonably accurate and commonly used for X-ray CT; however, the divergence of trajectories is not negligible in charged particle-related tomography. One drawback of this method is that it is mostly applicable to low-Z scattering mediums. However, in cases where the voxel size is larger, or in the same range as the displacement of the true path of the muon, it can show good results in practice. It can also be successfully employed in the case of high momentum muons with low scattering angles, though it assumes the ability for muon estimation. Moreover, estimation algorithms taking scattering into account show superior imaging quality when applied in CRT.

Most Likely Path
MLP is a probabilistic method to estimate the particle's lateral displacement and angular deflection at a certain depth [95]. It represents the most probable trajectory of muons that enter and exit at the same points and angles. Moreover, it constructs the MLP as a concatenation of these transverse positions by repeating the estimation in a series of depths across the medium. Thus, it exploits the fact that all MCS events are independent and that the angular dispersion due to MCS is approximately Gaussian.
A recent study [96] on CRT based on Bayesian theory and a Gaussian approximation of MCS explored the estimation of the most likely path of a cosmic ray muon travelling through a medium undergoing MCS when the muon's incoming and outgoing positions and directions are known. The study investigated the likelihood of a muon traversing an object with y as a vector that contains the lateral displacement and angle of a muon. From MCS theory, the scattering angle and lateral displacement follow a joint Gaussian distribution. The maximum a posteriori probability that a muon will have displacement and angle y j given the exit point is expressed according to the Bayesian theorem.
Then, the objective function is to maximise this probability with respect to y, to obtain the path for a muon as shown in equation below: where R A and R B are matrices that translate entry and exit points to the corresponding depth, and the covariance matrices Σ A and Σ B are also calculated at that depth.

Clustering Algorithm
There are several clustering algorithms: partitioning methods, hierarchical methods, grid-based, binned clustering, density-based clustering and others (see [97] for an extended review). High-Z materials are imaged with the PoCA method as regions of higher densities of unspecified shape compared with the background. It is therefore a natural choice to employ density-based clustering methods in tomography reconstruction. A binned clustering algorithm developed in [98] takes into account the degree of spatial clustering of muon scattering vertices. A higher density of vertices corresponds to high-Z materials, as strong muon scatterings take place with a greater frequency in such materials. The VOI is divided into cubic voxels and the scattering vertex for the muon is obtained at a minimal distance between the tracks. Furthermore, the scattering vertices inside each voxel are sorted by the scattering angle of the corresponding muon; the vertices corresponding to the largest scattering angles are retained. High values of n improve the contrast between high and low Z materials. For each of the pairs of vertices i, j in each voxel, a metric value d ij is calculated as: where θ i is the scatter angle of muon i and p i is the normalised momentum. Finally, the algorithm returns the median of the distribution log(d ij ) of all weighted metric distances as the final discriminator value for that voxel.

Angle Statistic Reconstruction
The ASR algorithm has a similar structure as the PoCA algorithm, but uses a different approach for determining which voxels are triggered by an event. One defines D r , the distance between the centre of a voxel and the incoming/outgoing muon track, as: where P 1,2 = p 1,2 + t 1,2 d 1,2 is the parametric equation of tracks, p 1 , d 1 and p 2 , d 2 are the known points and directions of incoming and outgoing tracks, respectively. and c is the centre of a voxel.
The D r will be compared to a threshold value d, and every voxel which satisfies D r < d will be triggered. The main advantage that ASR has over the PoCA approach is that every incoming and outgoing track will trigger at least one voxel, as illustrated in Figure 16. Where the PoCA fails at reconstructing low scattering angle events, ASR is able to extract relevant information by assigning a low scattering score to every triggered voxel. As a consequence, 3D scattering density maps obtained from ASR are smoother than with PoCA and are less sensitive to noise. A detailed analysis of the ASR reconstruction algorithm is reported in [99], showing the improvements in the reconstruction performance over PoCA.

Filtered Back Projection
Muon imaging techniques are naturally three-dimensional imaging methods. However, in some special cases, the three-dimensional imaging problem can be reduced to two dimensions. The filtered back projection (FBP) algorithm is often referred to as the convolution method using a one-dimensional integral equation to reconstruct a two-dimensional image. The FBP image reconstruction algorithm is based on the Radon transform [100]. In FBP, one applies a high pass filter to the projection information stored in the detector bins before back-projecting it into the space domain.
This method is not widely used for muon imaging purposes [101], because the FBP method requires a muon beam to calculate the muon beam attenuation and the scattering angle distribution width. This shortcoming can be accounted for by using a relatively long measuring time.

Algebraic Reconstruction Technique
The algebraic reconstruction technique (ART) [102] is an iterative reconstruction technique used in CT, also known as the Kaczmarz algorithm [103] in numerical linear algebra. It reconstructs an image from a series of angular projections (a sinogram). It allows for solving a system of linear equations with a large number of both equations and unknowns. Following the ART algorithm, we can solve the problem by sequentially projecting the point on each hyperplane at each iteration. An advantage of ART over other reconstruction methods is that it is relatively easy to incorporate prior knowledge into the reconstruction process. For example, current CRT reports an image reconstruction technique which involves the simultaneous algebraic reconstruction (SART) method [104], a superior implementation of the ART algorithm. In the SART algorithm, a single iteration is obtained by averaging the projections for all the hyperplanes.

Probabilistic Algorithms
The standard reconstruction algorithm is an EM method [73]. This technique iteratively maximises the posterior probability of the scattering density in the sample volume given the data. As a result, it can generate high-quality reconstructions from large amounts of data, yet can also perform better for shorter scan times. In particular, since it is a maximum likelihood technique, the accuracy of the reconstruction depends crucially on the accuracy of the statistical model on which it is based. To make the approach tractable, this has been kept rather simplistic. For instance, the algorithm ignores the uncertainty in our knowledge of the path taken by the muon or the path length in each voxel. Moreover, a noise model has to be determined for the specific design of the system, which can frequently change in a developmental system, so EM can prove impractical to implement effectively. EM also suffers because it is an iterative method, requiring significant optimisation and computing power to generate real-time reconstructions. Thus, a practical alternative to EM is desirable, which is fast, simple and has improved performance over simple algorithms such as PoCA. The implementation of EM in CRT is fully exploited [105] and represents a substantial improvement over previous heuristic methods.

Machine Learning Techniques
In very recent years, machine learning (ML)-based image reconstruction is seen as the third wave [106] in medical tomography, bringing enhanced imaging quality in comparison to conventional analytical and iterative algorithms. Two aspects in [107] point out the huge potential of ML for CRT: Firstly, artificial neural networks (NN) and similar supervised algorithms based on statistical learning are superior to the traditional reconstruction methods because there is no need to know the exact forward modelling (particle path); the modelling can be learnt from training if enough data are provided. This is particularly useful for scattering-based reconstruction, where modelling is only a rough estimation of reality. Secondly, for applications with a strict time limit, the resulting images are very noisy. It has been shown that the acquired tomography signal was increased 10-fold in low dose CT using deep convolutional networks. Few publications on ML applications in CRT can be found to date, which we summarise below.
To the best of our knowledge, the earliest published application of ML in CRT is the study by the CRIPT team reported in ref. [108] in 2012. They used only information from scattering, taking one input per voxel. Their simulation study, based on GEANT4, compared several supervised ML algorithms.
By means of a cluster analysis and a Markov Monte Carlo simulation to improve the PoCA algorithm (Section 6.1.1), the time of examining the existence of high-Z materials with the Bristol apparatus could be shortened to 60 (90) s to clear 64% (88%) of containers [109].
The Bristol team also published a study in 2016 [110], similar in concept to the aforementioned CRIPT study [108], where they only exploited scattering information. This study was not successful in demonstrating a significant advantage of ML over a simple method based on angles. They also verified the impact of momentum information, finding (somewhat surprisingly) that it makes no difference. The same team, however, proposed several methodological improvements in ref. [111], in particular through a multistep approach which sequentially applied clustering, edge-detection and finally material classification. The latter study was intended to serve as an example for material identification in the context of nuclear waste assays, but its conclusions can be equally applied to security applications.
In ref. [29], the TUMUTY team investigated to what extent drugs and explosives could be discriminated against using only muon scattering information. To differentiate various materials, an SVM classifier was used, a supervised learning model in machine learning. With 60 s of data, they reported probabilities of successful classification of 75%, 50%, 73% and 76% for flour, aluminum, steel and lead, respectively. In addition to the SVM classifier, the study also involved a polynomial kernel function and whitening transformation. Furthermore, they assessed a cases study where objects were placed side by side and vertically suspended, examining how the misclassification deteriorates, thus requiring longer exposure times.
ML, and in particular NNs, can also be used for image reconstruction. The image quality can be superior than other fast image reconstruction methods such as the FBP method discussed previously. In ref. [107], an artificial neural network classifier was used to improve the capability of imaging nuclear spent fuel casks. Yang showed that the acquired X-ray tomographic signal was increased at least by a factor of 10 in a low dose CT study by using a deep convolutional neural network. Moreover, it is possible to perform material classification by means of clustering methods used in machine learning as discussed in the paper [107].
A study based on determining the muon beam energies using multiple Coulomb scattering datasets with deep neural network (DNN) structures was performed with the aim to improve the muon energy resolution. It has been observed that the DNN significantly improves the resolution compared to the that obtained with the fit method, as described in the paper [112].
The binary classification technique of machine learning can be used to obtain a clear reconstructed image by subtracting the background from the data. On the other hand, multiclass classification can be used to discriminate and identify different materials under test using supervised machine learning algorithms. The image reconstructed using PoCA contains many false positives, which act as background and result in a smeared image. Using the binary classification technique of machine learning, the signal and background can be separated, thus resulting in a clear reconstructed image. Similarly, material identification can be performed using a multiclass classification technique, where each class corresponds to a different material. The total number of classes in multiclass classification is equal to the number of materials under which classification needs to be performed in addition to the background class [113].

Outstanding Issues of CRT in Security Applications
This section is devoted to critically assessing the main shortcomings of CRT in order security applications.

Statistical Limitation
The most fundamental issue for CRT is that it suffers from a much slower response (meaning the minimum data acquisition time needed to form a sufficiently discernible image of the cargo content to decide whether to inspect the cargo) than techniques based on other types of radiation, due to the modest natural muon flux of about 150-200 m −2 s −1 . The slow response may not satisfy typical border control requirements, which include a fast object detection. Radiation portals are able to provide screening times on the order of seconds, whereas the typical scanning time for a cargo with CRT is on the order of 1 to 30 min [14]. In spite of this, CRT portals can find their place in combination with more traditional scanning techniques, or as a second line of defence to further investigate ambiguous signals from other techniques.

Low-Z Identification Challenges
The ability of CRT to discriminate between elements is not limited to SNM. Important low-Z targets include narcotics, cigarettes and explosives. The prevention of human trafficking, and hence the detection of human bodies, can be seen as belonging to the same low-Z identification problem. However, the movement towards low-Z material identification has been slow.
Low-Z materials hidden within other high-Z ones, or even worse within other low-Z ones, pose significant challenges to CRT due to the small average deflection. This makes it a much more challenging problem than high-Z identification, for which CRT was originally designed.
Due to the Z dependence of the scattering density parameter, discussed in Section 2 and illustrated in Figure 17, identification of low-Z materials by muon scattering within a limited measurement time is technologically challenging, due to the small average deflection that demands excellent positional resolution. To date, very few research teams have published dedicated investigations on low-Z detection. In 2009, Cuellar et al. [114] reported results based on simulated data, concluding that the low momentum part of the charged leptons spectrum (including electrons and positrons) could be used based on their stopping effect as additional information in scattering tomography to help resolve medium and low-Z materials. The proposed effect was demonstrated by experimental modelling as well. Simulated results with a material size of 1 m 3 (paper and nylon), and experimental results for pallet-size low-Z materials were presented with 30, 10 and 5 × 2 min durations. They concluded that the proposed 'stopping power' method is very useful but geometry dependent.
The TUMUTY [29,30] team performed simulated and real data tests with four materials: flour, aluminum, steel and lead, chosen as representative of a wide range of Z and density values, with flour in particular being considered fairly similar to drugs regarding its relevant physical properties.
Decision Sciences [16] claims the ability of its DISCOVERY system to identify various low-Z materials in security control contexts. Examples reported in ref.
[14] include the discovery of 1900 kg of marijuana (with a total street value of approximately USD 20 million) hidden inside metallic rolls through CRT images acquired in 4.5 min, while X-ray inspection had not provided any actionable information due to its inability to penetrate the steel. They also demonstrated, through test scenes with data acquisition times of 10 min, the ability to recognise a cocaine surrogate hidden inside banana pallets and different chemicals used as surrogates or precursors of explosives hidden inside pallets of bottled water and charcoal [14].
In conclusion, searching for small-Z materials in shipments poses formidable challenges to CRT due to the small average deflection with respect to the high-Z identification problems that are more customary for this technique. This issue will be overcome by the dedicated work of a few teams. This issue can be potentially overcome by increasing the detector resolution, and by including additional experimental observables in a multivariate discrimination, which may include an estimation of the incoming muon momentum (Section 7.3) and the detection and tagging (i.e., particle identification) of additional cosmic ray particles or muon-induced 'secondaries' to complement the muon signal (Section 5.5).

Muon Momentum Estimation
As described in Section 4.2.2, scattering amplitude scale inversely with momentum, p (or, equivalently, energy, E). Thus, a low-energy muon going through a low-density material can mimic the behaviour of a high-energy muon traversing a high-density material. Including momentum measurements in reconstruction algorithms would allow for reaching new levels of precision while reducing the acquisition time. Even though momentum measurements based on tracking particles in a magnetic field are common practice in high energy physics, maintaining a sufficient magnetic field in a large volume comes with a very high cost, making its usage commercially unrealistic in MST. Efforts have been made to develop muon momentum measurement modules based on multiple Coulomb scattering [21,115,116], Cherenkov detectors [117] or time of flight measurements [118].
The RMS of the muon scattering angle distribution, as shown in Equation (2), depends on the muon momentum, its path length in the material and the material radiation length. Thus, by measuring muon deflections through successive known thicknesses of a known material, as shown in Figure 18, one can estimate its momentum. As an example, the CRIPT spectrometer [21] utilises alternating scintillator detection planes and 10 cm thick steel slabs, which act as a scattering material. The deflection measurements are combined with an energy loss measurement, and both are used as inputs to the momentum inference algorithm. The same momentum estimation technique is also used in the INFN demonstrator [19,116]. In practice, information related to scattering can be extracted in various ways, e.g., by the angle between the initial (line connecting the first two hits) and final (line connecting the last two hits) trajectories of the same muon, the χ 2 of the track [116], a Bayesian approach as in [96], or a Kalman filter fit where p is a parameter [119].
Since muon deflections are reconstructed from position measurements, the spatial resolution of the detection system drives the precision of the momentum measurement. In order to overcome this limitation, Ankowski and others [120] proposed the use of a Kalman filter to disentangle deflections originating from the scattering itself from 'ghost' deflections induced by the detector resolution.
Overall, the scattering-based momentum measurement technique performs well at classifying low momentum muons (1 < p < 5 GeV/c), since they undergo large deflections that can be easily measured, but reaches its limit with high energy muons (p > 5 GeV/c), for which deflections are smaller than the angular resolution of the detector. Sub-GeV muons might also be stopped in the scattering material and not reach the last detection plane, causing the loss of the event.
NEWCUT [115] is a dedicated experiment in Japan for the measurement of differential spectra of cosmic ray muons, focusing on the low-energy region that is most important for any cosmic muon imaging applications, including but not limited to security (in fact, the main motivation of the authors is to improve the resolution of their volcano studies). It consists of 19 MWPC and 10 to 15 lead slabs, the latter acting as absorbers/scatterers. An interesting feature of their methodology is that scattering-related information is combined with the energy lost by the muons at each plane traversed, as measured through the signal amplitude induced in the wires of the MWPCs, which is proportional to dE/dx, which is in turn sensitive to momentum when the latter is not larger than a few GeV. A neural-networkbased (NN) regression algorithm is used for combining the input observables, achieving a good accuracy for the range of p < 6 GeV/c, Figure 19. At higher momenta, a saturation effect is observed, see Figure 19B, meaning that the NN has insufficient information to predict the momentum in that regime. This is due to the fact that dE/dx has reached a plateau, and the energy dependence of the scattering has reached a point where the typical deflection angle is smaller than the angular resolution of the detection system. The momentum threshold above which a muon emits Cherenkov radiation is function of the refractive index of the medium it is traversing. By tuning the gas mixture and pressure of the medium, one can choose the momentum threshold, p th , of a Cherenkov radiator module. A muon traversing a Cherenkov radiator with p > p th will emit Cherenkov radiation and release an optical signal. Assembling radiators with increasing p th , Bae J. and Chatzidakis S. [117] were able to categorise muon momentum into several momentum levels from 0.2 GeV/c to 7 GeV/c, with an interval of 1 GeV and a true classification rate around 85%.
Momentum can also be measured from the time of flight (TOF) of a muon through the detector. This solution demands an excellent timing resolution as well as a long travel distance inside the detector (as in [118]). For example, travelling 1 m takes 3.35 ns for a 1 GeV muon and 3.39 ns for a 0.5 GeV muon. Therefore, it takes a state-of-the-art time resolution of 50 ps (challenging but achievable with some detector technologies, such as scintillators and multi-gap RPCs) to give a fair p measurement for sub-GeV muons. The TOF can be a useful input feature to be combined with other parameters in a multivariate analysis.
By themselves, scattering, Cherenkov signals and TOF lead to large uncertainties in p. Their statistical combination would reduce the uncertainty, as they are statistically uncorrelated, due to their different sources of systematic uncertainties. These alternative ways to access the same physical information may be sufficiently complementary; thus, their combination by machine learning could be significantly more advantageous than using them alone.

Conclusions
With a significant amount of illegal materials trafficked through our borders, border security teams are under increased pressure to track and seize illicit cargo. Not all lorries or containers are currently inspected and only a small fraction of illegal goods are found. Current scanning systems at ports and other country entry points (mainly X-ray scanners) cannot always detect or identify contraband, in particular, if they are hidden behind valid and properly declared objects. The CRT technique for border security has been developing since the early 2000s and has the potential to complement existing technologies. Several projects are ongoing to optimise this method, including the choice of muon detectors and their configuration; development of fast muon reconstruction algorithms, image processing and user-friendly interfaces; and the integration of muon scanners into border control facilities. Analysis of muon scattering has also recently been complemented by studying muon absorption and detection of other particles to improve material identification. One of the CRT projects that has received funding from the Horizon 2020 funding scheme is SilentBorder [64]. Although the potential of CRT as a cargo scanner has been demonstrated with small setups and some models tested at ports; generally, CRT projects still need to reduce their scanning/response time and improve the discrimination between legal and illicit materials with similar compositions and densities. Combining data from CRT and X-ray scanners and applying artificial intelligence for automated analysis of images will allow improved detection efficiency for hidden or non-declared contraband.