Multiparametric approach to the assessment of muon tomographic results for the inspection of a full-scale container

Experimental results from a dataset collected with a full-scale muon tomograph for the inspection of cargo containers were studied in a single scattering scenario with a multiparametric analysis based on the method of the Point Of Closest Approach (Poca). To search for high-Z materials, a 4 dm3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {dm}^3$$\end{document} Pb block was positioned inside the volume to be inspected, in order to quantitatively investigate the appearance of the Poca signal. Signal-to-noise ratio and significance of the Poca signal were investigated by means of mono-dimensional spectra of the Poca components, for different values of the scattering angle between the incoming and outgoing muon tracks and with different angle-dependent weights. A systematic scan of two-dimensional maps was also carried out, as a strategy to search for possible enhancements to the Poca signal. A comparison was also done between the results obtained from the two half-volumes, one containing the Pb block and one left empty, to take into account the response of the detector and some aspects of the Poca strategy.


Introduction
The detection of illicit nuclear materials (e.g. Uranium, Plutonium) which could be transported inside cargo containers, escaping the standard control procedures, has received a large attention during the last years, due to the potential hazard arising from the huge number of containers in transit across the world. Inspection techniques based on the scattering of cosmic muons have been suggested [1] as a promising strategy to provide reliable monitoring systems to detect the presence of hidden quantities of such materials inside large volumes, a e-mail: Francesco.Riggi@ct.infn.it (corresponding author) even when these are embedded in a bulk of objects mainly composed by low-density, low-Z elements. Muon scattering is indeed particularly sensitive to the presence of high-Z elements, thus allowing for a discrimination between low-Z and high-Z materials, by looking at the scattering angle between the incoming and outgoing trajectories of these particles, once they are properly reconstructed inside a tracking detector above and below the volume to be inspected. An important advantage of such technique is the worldwide availability of such natural secondary cosmic ray radiation, whose flux, momentum and angular distribution are reasonably well known at the sea level or moderate altitudes. As compared to traditional radiography techniques based, for instance, on the use of artificially produced X-or γ -rays, or even to neutrons, the use of cosmic ray muons is somewhat limited by the unmodifiable and relatively low impinging flux, thus requiring a longer scan time, which is usually a point of concern for any practical large-scale application.
Notwithstanding, several projects worldwide have designed and built prototypes of muon detectors for such an application, from small size devices to full-scale systems. Various detection technologies have been employed, from scintillator strips [2] and scintillating fibres [3] to gas-based detectors, such as drift tubes [4], drift chambers [5], Resistive Plate Chambers (RPC) [6,7] and Gas Electron Multipliers (GEMs) [8]. Even semiconductor sensors have been employed for the construction of a small size prototype [9]. Some examples also exist of industrial attempts of providing muon tomographs specifically designed for harbour sites [10,11].
In parallel with the development of the hardware needed for such detectors, a large effort in the implementation of existing and new methods of analysis and imaging algorithms has also been pursued, to provide an efficient way of determining the presence of hidden materials from the measured raw data, and to improve the sensitivity and the resolution of the method, with also the goal of ultimately reducing the amount of time required to scan each individual container. Important benchmarks of the algorithms and software methodology are also the amount of false positives and the failure probability in the detection of dangerous materials, since both contribute to the overall efficacy of the technique.
As part of a local project we have built a full-scale prototype (18 m 2 sensitive area) of a muon tomograph with all features to be potentially employed for the inspection of a standard Twenty-Foot Equivalent (TEU) container, a volume of about 6 m × 2.4 m × 2.6 m. A description of this facility and several aspects of its construction and test have been reported in previous papers [12][13][14][15][16][17][18][19][20][21][22]. The main aim of this project was to explore all the main aspects of a large-scale installation for muon tomography, especially the ones related with building a detector with a large sensitive area and a modular structure, operating a large number of photosensors and electronic channels (O(10 4 )), controlling the data flow and the operational conditions even by remote, testing simulation procedures of the entire setup and comparing different reconstruction and imaging strategies. Some of these aspects have been already addressed in our previous papers, already cited.
In this work we analyzed a dataset collected with the full-scale detector, with the aim of evaluating the usefulness of a multiparametric approach used in connection with a simplified reconstruction algorithm based on the Point Of Closest Approach (Poca) method, which is known to be suited to single scattering scenarios. The experimental setup is briefly described in Sect. 2, while the basics of the Poca method are recalled in Sect. 3. Section 4 reports the various analysis strategies, used in order to provide different mono-and bi-dimensional spectra of the Poca components, and an assessment of their performance. Possible improvements to the technique, aiming at the reduction of the required scan time, are discussed in Sect. 5. Some concluding remarks are finally presented in Sect. 6. Fig. 1 The muon scattering process, particularly sensitive to the presence of a scattering center, may be used to provide a muographic image of a hidden volume. The scattering angle is defined by the incoming (a) and outgoing (b) track directions, reconstructed by an upper and lower muon tracker. Eight 3 m × 6 m detection planes (4X and 4Y) are employed in the setup presented here. The red and blue surfaces denote the X and Y detection planes of the upper and lower trackers, segmented into equal 1 m × 3 m modules, as sketched on the right part of the figure

Experimental setup
A description of the Muon Portal project is included in Ref. [22] and references therein. Here we only briefly recall its main components.

Geometrical structure of the Muon Portal
The detection setup is based on four position-sensitive detection planes (henceforth called "logical" planes), two placed above and two below the volume to be inspected. Each logical plane is made up of two physical planes (X and Y). Each physical plane in turn is made up of six detection modules (1 m × 3 m each) installed in such a way to cover both the X-and the Y-coordinates by the same type of modules, as shown in Fig. 1, leaving a negligible dead area between adjacent modules. A total of 48 identical detection modules are used in the Muon Portal. A motorized mechanical structure accommodates all the detection modules on four metal grids (one per logical plane), also allowing to modify the vertical position of the various planes. Figure 2 shows a picture of the muon tomograph during the installation, with some of the detection modules already in position. For the measurements discussed here, the logical planes were positioned at Z = 0, +135 cm, +405 cm and +540 cm respectively. The entire Muon Portal was installed in a big hall, with an approximate size of 7 m × 6 m × 7 m, equipped with an air conditioning system in order to maintain the room temperature within an acceptable range, especially during summer. Special care was also taken to reduce the amount of light entering from outside, which is of paramount importance for the photosensors being used.

The detection modules
Each of the 48 detection modules is segmented into 100 strips of Amcrys extruded plastic scintillator (1 × 1 × 300 cm 3 ), painted with a reflective TiO 2 coating, with two Wave Length Shifter (WLS) Kuraray Y11 1 mm fibres embedded inside two grooves in each strip, to transport the photons to the photosensors (one per fiber) placed at one of the fiber edge, with an optical cookie for a stable coupling between each photosensor and the corresponding fiber. The other end of the fiber was polished and covered with a reflective coating. The upper surface (3 m 2 ) of the strips was also covered by an aluminium adhesive mirror tape. All the strips in a module were finally stretched between two aluminum layers (1 mm thickness). The overall number of channels (fibers and photosensors) is 9600. Extensive tests of the individual strips and WLS fibres have been reported in previous papers [16,17].

Photosensors
The photosensors employed are silicon photomultipliers (SiPM), specifically designed by STMicroelectronics to provide a good photon detection efficiency and cell fill factor, as well as to ensure a low cross-talk and dark count rate. Different SiPM prototypes both with the p-on-n and n-on-p technologies have been produced by STMicroelectronics for this project before the mass production. The layout of the final chip was based on n-on-p technology and embedded 4 independent round-shaped SiPMs (about 2 mm 2 sensitive area), integrating two types of SiPMs in pairs, with a slightly different number of cells and pitch. The model MUON60, with a cell pitch 60 µm × 60 µm, 548 cells and cell fill-factor of 67 % was finally employed for the project. A custom procedure was implemented for the mass characterization of all the 9600 devices, using a custom software implemented with LabView to measure the current-voltage (I-V) curve and find the Breakdown Voltage value for each individual SiPM. Detailed results on this procedure and the results obtained are reported in Refs. [18,20].

Electronic readout and data acquisition
Considering the number of modules and the required granularity necessary to achieve reasonable space and angular resolutions, a large number of channels was required for this detector. A compression technique within each detection module (100 strips) was employed in this prototype, to reduce such number by a factor 10. This is achieved by the use of the two WLS fibers embedded in each strip, optically coupled to an equal number of SiPMs. A proper combination of the signals originating from ten channels grouped together (giving the Tens and the Units) was used to identify the interested strip [15]. Although this largely reduced the amount of electronics required, a net effect is observed on the light collection efficiency, since the scintillation photons produced in each strip are shared between the two fibers.

The point of closest approach (POCA) strategy
The simplest algorithm in Muon Tomography is the Point Of Closest Approach (Poca), which does the simplified assumption that muon scattering occurs at a single point. It therefore searches for the geometrical point of closest approach between the incoming and outcoming reconstructed track directions (see Fig. 1), P POCA = (a + b)/2, where a and b are the vectors defining the direction of the incoming and outcoming tracks, respectively, and P POCA are the coordinates of the Point Of Closest Approach to be found. Such a method is easy to implement and provides results in reduced computation time, useful as a first-order approximation of the solution or as a complementary strategy to be combined with more sophisticated algorithms. It is well known that the Poca approach, since it neglects multiple scattering processes through the volume material, has the drawback of providing low-resolution images, it is quite sensitive to the presence of shielding materials located above or below the potential threat and cannot localize very well materials close to the volume borders [13,23]. Moreover, this method may discard part of the tracks, when the intersection of the incoming and outgoing tracks is located in a voxel which lies outside the volume of interest. This has motivated the implementation of several alternative imaging methods [13,19,[23][24][25][26][27][28][29]. Among others, log-likelihood algorithms, two-point autocorrelation analyses, density-based clustering algorithms and friends-of-friends strategies have been sometimes used, even by ourselves [13,19], since they are based on more realistic physical and statistical assumptions, and allow to tackle the problems encountered with the Poca algorithm. These algorithms have been frequently employed as test tools in the analysis of simulated data, in various scenarios for the location of heavy materials inside a container. However, in most cases only a qualitative evaluation of the imaging results has been done, without a quantitative assessment of the possibilities offered by the various reconstruction strategies.
The aim of this paper is not to compare the Poca approach with more sophisticated strategies, but rather to show that even within the Poca strategy, a multiparametric analysis and the use of multidimensional views of the Poca components may help to largely improve the search of high-Z objects. This approach has been applied to a set of experimental data available from the Muon Portal project, rather than to simulated data, since the limitations and response of a real prototype detector need to be taken into account for any performance evaluation of the algorithm used.

Experimental conditions and event reconstruction
During the Muon Portal commissioning phase several sets of measurements were performed, with a trigger given by the coincidence between the eight detection physical planes (4 for the X-and 4 for the Y-coordinate). Thanks to the architecture of the electronics employed, the two half-planes, left (X = 0-300 cm) and right (X = 300-600 cm), along the 6 m long X-side, could be operated as two independent 3×3 m 2 detectors. For the present investigation we employed such a configuration, also in order to compare the results obtained from the right half-volume, where a Pb block of 4 dm 3 was positioned, to that obtained from the left half-volume, which was left empty. Bias voltage to the SiPMs and individual thresholds were adjusted in such a way to have an average dark count rate (DCR) of about 1 kHz from each SiPM, relatively low with respect to the intrinsically high DCR for these devices (10 6 Hz) when operated at a threshold corresponding to 1 photoelectron (p.e.), but still too high to reach a detection efficiency close to 100% for the single channel. From previous measurements in the lab and from the measured value of the cosmic ray flux, this was estimated to be close to 80%. The global effect of this efficiency on the performance of the entire detector is discussed in Sect. 5.
For the present measurements we inserted a set of Pb blocks with an overall volume of 4 dm 3 (20 (X) × 20 (Y) × 10 (Z) cm 3 ) in between the two inner tracking planes, at Z = 215 cm with respect to the bottom plane, almost at the center of the right half-volume (X = 445 cm, Y = 160 cm), positioned on a light metal stool as a mechanical support, shown in Fig. 3. All the remaining volume delimited by the two left half-planes was left empty.
Raw data consist in a list of hit strips in each detection plane, either along the X-or Y-coordinate. Contiguous hits are grouped into clusters to perform tracking reconstruction.
In a muon tomograph, each track is usually reconstructed correlating the hits measured over several detection planes, at least two above and two below the volume to be inspected. This is also the case of our experimental setup, where the hits are independently measured over 4 physical detection planes along the X-coordinate and 4 detection planes along the Y-coordinate, as described above. Track quality in a particle tracking detector is generally evaluated by the number of clusters contributing to that track and to some χ 2 -like parameter, which quantifies how close are the clusters to a straight line.
For muon tomographs which have several (more than just two) detection planes in the upper and lower parts of the setup, one could impose a 'good track' quality requirement on each of the two independent reconstructed segments. However, in most of the cases (including our setup), only two logical planes are available above and below the hidden volume, so the muon track is reconstructed by the use of 4 clusters only, one per plane.
Therefore a very preliminary selection of cluster combinations (i.e. possible track candidates) was first performed requesting two conditions: (i) an average residual, with respect to the straight line fit (performed considering all four points), smaller than 10 cm; (ii) the Poca coordinates confined inside the scanned volume. We checked, either modifying the value of the cut on the average residual or from GEANT simulations, that this is a reasonable compromise to proceed with the correct association of clusters to a muon track. A larger value of the average residual very unlikely may originate from a genuine muon track being deflected by a high-Z scattering center and would introduce a large number of wrong cluster associations, hence a large number of fake tracks.
In the operating conditions of the measurement used for this study the average cluster multiplicity, i.e. the average number of reconstructed clusters per event in each half-plane, was about 1.5, with a fraction of events having more than two clusters per event smaller than 10%. Moreover, clusters characterized by a cluster size larger than 2 (about 1% of the total) were not used in the tracking reconstruction procedure, since they are unlikely to originate from muon tracks. This resulted in an average number of cluster possible combinations of about 6 and the reconstruction of tracks was simply achieved by selecting the track with the minimum χ 2 among all the possible combinations. It was verified by Monte Carlo simulations that approximately 1 in 10 6 combinations of clusters may be randomly aligned and pass the selection criteria cited above.
An off-line alignment procedure was also undertaken, keeping fixed three of the four (X or Y) coordinates in each logical plane, and shifting the last coordinate by small amounts until the χ 2 was minimized. This procedure was repeated for all the detection planes, and misalignment corrections were incorporated in the tracking and image reconstruction procedure. The results showed that misalignment corrections were of the order of a few mm, compatible with the typical space resolution (about 3 mm) of each plane.
A total of about 690k muon tracks were reconstructed, in about 2 days data taking. Among these, 397 k events gave a set of Poca coordinates within the volume defined by the two inner tracking planes, between Z = 135 and Z = 405 cm, almost equally distributed between the two half-volumes.

Scattering angle distributions
Once the track was reconstructed by means of its four clusters, the scattering angle, defined as the angle between the two track segments in the upper and lower detectors, is a quantity somewhat correlated with the average residual, as seen in Fig. 4, which shows a scatter plot of these two quantities for all the tracks crossing the detectors and producing a Poca information located inside the overall volume to be explored. For most of the tracks, which pass through the volume without suffering large scattering effects, the average residual is between 1 and The overall scattering angle distribution is determined by the response of the detector, which includes its tracking capability, space and angular resolution, and any disuniformity in detection efficiency, since most of the muon tracks do not cross heavy scattering centers, but simply traverse the inner volume. In our measurement the average value of the scattering angle was (3.001 ± 0.006) • for the inclusive data sample (Fig. 5, top). An increase of the average value of the scattering angle is expected within the small region where a high-Z material was located, as shown in Fig. 5 (bottom), where a sensibly higher value, (5.2 ± 0.2) • , was obtained for the subset of tracks with a Poca located inside a volume of 20 (X) × 20 (Y) × 40 (Z) cm 3 surrounding the place where a 4 dm 3 Pb block had been positioned.
The scattering angle is then a useful control variable, and one can think to locate the suspect material by searching for voxels which exhibit an anomalously high value of this quantity. With a voxel size of 20 cm along the X-, Y-and Z-directions, we have 15×15×14=3150 voxels in each semi-volume, which can be individually scanned, evaluating the average value of the scattering angle and its uncertainty for each voxel. The voxel ordering was chosen such as to have voxels 1-15 (16-30,…) identifying the cells along the X-direction for the first (second, …) Y row, and so on, with 225 voxels in each Z section. With the minimal condition that each voxel must contain at least two entries, the distribution of the average scattering angle over all the voxels of the right half-volume (X = 300-600 cm) is shown in Fig. 6 (top).
One may first observe that the distribution of the average scattering angle inside each voxel has a pattern with large fluctuations. Some periodicity, corresponding to the number of voxels (225) in each Z-section, and a general trend with a minimum near intermediate Z-values may also be noticed. Both are consequences of the geometry and segmentation of the detector. A large fraction of the voxels exhibit a scattering angle larger than the average, although with a large uncertainty, since the number of entries which contribute to determine the mean scattering angle inside each voxel is in most cases very small.
For this reason, a filtering of such distribution was performed, by considering both the number of entries in each voxel and the significance of the difference between the average scattering angle evaluated in each voxel and the overall average scattering angle computed The value of the threshold on the number of entries per voxel clearly scales with the total collected statistics, and for our dataset was fixed to 160. After applying a 5σ cut on this difference and selecting only those voxels which have a number of entries larger than the assumed threshold value, the distribution shown in Fig. 6 (bottom, blue symbols) is obtained. Comparing the top and bottom distributions in Fig. 6 we see that only a few voxels survive these cuts. The two spatially contiguous voxels (numbered 802 and 1027, differing by exactly 225 voxels) correspond to the location of the Pb block in two adjacent Z-sections. A few voxels with a false positive are however present in the lowermost Z section. Considering the median of the scattering angle distribution instead of the mean is another option, whose results are also shown in Fig. 6 (bottom, red symbols). The use of the median, which was estimated to be 1.87 • , was seen, at least in this case, to reduce the number of false positives, after applying the same criteria concerning the threshold on the number of entries per voxel and the 5σ difference.
The number of entries in each voxel, even without large scattering effects, is expected to be dependent, due to purely geometrical acceptance considerations, upon the position of the voxel in the XY plane and, to a minor extent, also on the Z-section along the vertical axis. More sophisticated strategies could make use of a threshold dependent on the location of the voxel (for instance its Z-coordinate, or the distance from the center of the detector in the XY plane). We also did some attempt to further improve the analysis strategy, considering a local average of the scattering angle (instead of the global average), i.e. a reference value dependent on the location (X, Y and Z) of the voxel, but we realized that these strategies could be too much adapted to the particular dataset and not of general use.
Several tests, carried out with different thresholds on the number of entries per voxel and with different cuts on the average scattering angle, demonstrated that the presence of the Pb block is always correctly identified, while the amount of false positives may slightly depend on the cuts chosen. By applying the same strategy (with strictly the same cuts) to the left half-volume (X=0-300 cm), only some false positives in the low Z-region are found, without any evidence of signal in the intermediate region. Since the presence of false positives is also related to the detailed response of the detector, it is advisable that a set of measurements is carried out periodically with an empty volume, for calibration purposes. We verified that a similar measurement, carried out a few days later, gave the same response concerning the location of false positives, still evidencing the presence of the Pb block in the right half-volume. The presence of false positives, which in this context also originates from disuniformities in the detector efficiency, may be then considered an intrinsic component of the detector response and could be incorporated in the scanning procedure.

Mono-dimensional POCA spectra
Mono-dimensional distributions of the individual X-, Y-or Z-components of the Poca information, computed selecting different slices in the other two components, are a complementary way to search for meaningful signal enhancements over the background.
To have a realistic estimate of such possible enhancements, Fig. 7 shows the distribution of the Poca X-component, after selecting a YZ region of 40 cm × 40 cm (140-180 cm in Y, 195-235 cm in Z) surrounding the position where the Pb block (20×20×10 cm 3 ) was located (top plot), and the corresponding spectrum of the same component, selecting a different, but with an equal size, YZ region (100-140 cm in Y, 235-275 cm in Z) around an empty volume (bottom plot). A lower cut on the scattering angle, ϑ > 5 • , was imposed, to enhance those regions characterized by a large scattering angle. The effect of choosing different values of such a cut is discussed further below. The number of events in the two spectra cannot be Fig. 6 Distribution of the average scattering angle versus voxel number, as defined in the text, without any constraint apart from requiring at least two entries per voxel (top). Introducing a threshold on the minimum number of entries per voxel (160) and a 5σ cut with respect to the overall average scattering angle, the bottom plot (blue symbols) is obtained, where the location of the Pb block in two spatially contiguous voxels at 802 and 1027, is correctly identified, together with a few false positives distributed in the lowermost Z section. In the bottom plot also the median, instead of the mean, of the scattering angle distribution in each voxel, was considered (red symbols). This did not alter the evidence of the voxels corresponding to the location of the Pb block, but somewhat reduced the number of false positives in the lowermost zone directly compared, since they refer to different regions, although with equal volumes. The presence of a peak around Poca X = 445 cm in the top spectrum, corresponding to the location of the Pb block, is clearly visible, peak which is not present in the lower spectrum. A similar result was observed in the Poca Y component, while only a very broad peak is observed along Z, as expected, due to the fact that the upper and lower segments of the muon track are nearly parallel to each other, thus resulting in an elongation of the reconstructed image along the vertical axis, a well known effect when using the Poca method.
The visibility of a peak in the spectra of the Poca components may be quantified by the signal-to-noise ratio and by its significance. To search for a meaningful enhancement in these variables across the overall volume to be inspected, an array of 210 histograms was considered for the X-component of Poca, considering Y and Z bin sizes of 20 cm (15 sections over Y and 14 over Z). For each histogram, the signal-to-noise ratio S/N and a simple estimator of its significance, S/ √ N , were extracted, considering the signal contained in the bin with the largest number of events in the histogram. The background was estimated as the average of all remaining voxels in the histogram. Also in this case a cut ϑ > 5 • was introduced. Signal-to-noise ratio and significance are plotted in Fig. 8 for the 210 histograms. The top Fig. 7 Top: Distribution of the Poca X-component when selecting an YZ region around the volume where the Pb block was located. Bottom: Background distribution, obtained selecting an empty, equal volume, YZ region. Only events with a scattering angle ϑ > 5 • are considered here plot shows that for several of them (i.e. several bins in the Y-and Z-components), at least 4 or 5, there is some evidence of a large signal-to-noise ratio. However, the significance of the peaks found in each of these histograms is larger than 10 only for two of them, which are actually close in space and correspond to the location of the Pb block in two contiguous bins along the same vertical direction, in two adjacent Z-sections.
An additional, more compact, strategy to look for strong enhancements both in the signalto-noise ratio and in its significance, is to use a scatter plot of these two variables, such as those in Fig. 9, where different values of the lower cut on the scattering angle were considered. Figure 9 clearly shows how the two close points lay in the upper right part of the plot, thus indicative of a meaningful evidence of a well identified peak in the Poca spectrum. Realistic threshold values in the range 4 • − 5 • were seen to give the best results, resulting in a well isolated pair of points with large signal-to-noise ratio and significance values. Similar results Fig. 8 Signal-to-noise ratio (top) and significance (bottom) along a series of 210 Poca X histograms spanning the overall volume to be inspected, for different slices of the Y-and Z-Poca components were also obtained for the Y-component of the Poca. A careful and systematic inspection of all these spectra may then be considered as part of a reasonable strategy to search for evidence of scattering centres from heavy materials.
Even if the optimal value of the cut on the scattering angle is not known a priori, however this value lies in a narrow range, say from 3 • to 6 • . This is also confirmed by GEANT simulations. Smaller values do not evidence high values of both signal-to-noise ratio and significance, while very high values of this cut risk to lower too much the statistics, cutting most of the good events (as can be seen in Fig. 9 when requesting a cut ϑ > 5 • ). The idea behind this analysis is that data processing of the events collected by a muon tomograph could be employed in real time to produce many plots of interest, providing several views of the various distributions, in order to help the user to judge whether there is a potential hazard. Plots like those shown in Fig. 9 could be typical control plots to be shown during the scan procedure. Of course a prompt and automatic detection of abnormal values of some parameter is also advisable, with easy-to-see symbols or alarms. We made some tests concerning the use of such tools, also considering the possibility to process the data in parallel  Signal-to-noise ratio as a function of the cut on the minimum scattering angle. Blue dots and red squares report a comparison between the results obtained weighting the entries in each Poca X histogram with a weight proportional to ϑ 2 and without any weight respectively. Also shown in the plot is the value of the significance (black triangles) on different processors and visualizing the plots of interest on multiple screens. The detailed implementation of this approach needs to be customized for any specific application. As for Fig. 9, one may find some metrics combining the signal-to-noise ratio and significance values, like for instance their normalized quadratic sum, and produce an alarm for high values of this quantity, evaluated after a fine scan across the range of possible cuts on the scattering angle. To investigate the effect of the cut on the scattering angle and make a comparison of the results obtained assuming no weight or a weight proportional to ϑ 2 , we repeated the same analysis on the X-component of the Poca, already shown in Fig. 7, with the same windows on the Y-and Z-components, to extract the signal-to-noise ratio and the significance of the peak corresponding to the location of the Pb blocks. The results are shown in Fig. 10. Blue dots and red squares show the signal-to-noise ratio as a function of the lower cut on the ϑ scattering angle, weighting the entries in each Poca X histogram with a weight proportional to ϑ 2 and without any weight, respectively. While a weight proportional to ϑ 2 is effective when no ϑ cut is introduced, its effect is almost useless for realistic ϑ cuts (of a few degrees), so that the signal-to-noise ratio obtained in the two cases is comparable. The significance is reported as black triangles in the upper part of the plot, and shows that cuts larger than 4-5 • start to worsen it. Therefore, as already outlined, a realistic choice in such a situation could be a cut of the order of 4-5 • , which maximizes both the signal-to-noise ratio and the significance of the observed peak. Optimization procedures should be investigated for each specific experimental setup, in order to find the best conditions to process the collected data.
The previous results were obtained using the total collected statistics (about 200 k reconstructed Poca events inside the right half-volume). A meaningful evidence of a peak in the signal-to-noise ratio may be obtained even with a smaller number of events, hence in a shorter scan time. The evolution of the signal-to-noise ratio is shown in Fig. 11 as a function of the number of Poca reconstructed events. The plot shows that starting from a number of events of the order of 20 k, i.e. a factor 10 smaller than the overall statistics, a signal-to-noise ratio sensibly larger than one is observed. In the present condition, this corresponds to a few hours scan time, which is still too high for most practical applications. A more detailed discussion of this aspect is reported in Sect. 5.

Two-dimensional POCA maps
Examining two-dimensional XY, XZ and YZ maps may provide additional views of the scattering center distribution, allowing for a quantification of possible signal enhancements over the average background, and pointing out the presence of potential hazardous materials. 2D sectional views of 3D distributions are more easily to visualize than full 3D representations, which require sophisticated graphical tools (such as for instance fast rotation 3D techniques) and usually suffer from the difficulty to proper display hidden parts of the distribution. On the other hand, the analysis of 2D sections requires the comparison of different plots, possibly on multiple screens, in order to allow for a direct comparison between different views. With our voxel binning of 20 cm in all coordinates and an overall volume of 300 × 300 × 280cm 3 we should explore a total of 44 two-dimensional maps, which may be organized over several visualization windows either in a static mode (looking at different plots in different parts of the monitor or in different monitors) or even in dynamic mode, flushing in sequence the various maps over the same portion of the screen. Figures 12 and 13 show some of these XY and XZ maps, for a few intervals in Z, by using the surf3 graphical ROOT option [30]. Similar results are also obtained for the YZ maps. The vertical scale was normalized in each case to the map with the highest number of events. Clear evidence of a peak is seen in two of the maps (Fig. 12), corresponding to two sections in Z, and in one map (Fig. 13) corresponding to the location of the Pb block in Y. As expected, the image is more elongated in Z than along X or Y, due to the intrinsic nature of the Poca algorithm. In all these maps a cut on the scattering angle (ϑ > 5 • ) was used. We verified that slightly different values of this cut (4 • − 6 • ) gave always evidence of the observed spots.
Finally, we compared the results of two-dimensional maps of the Poca events measured in the two half-volumes (right and left) with the aim of evaluating the difference between the two XY distributions. The same cut on the scattering angle (ϑ > 5 • ) was employed, as for the previous plots. As an example, Fig. 14 shows some of these maps in different sections of the vertical coordinate. While background is smoothed and reduced by using the

Scan times and possible improvements to the detection efficiency
After an R&D phase to test and choose the individual components, the construction of a full-scale Muon Portal detector prototype was completed. The goals of this prototype were to understand the main aspects of building a large size (about 120 m 3 ) detector, able to accommodate a full-scale 20-foot container, to test all the detector components, from scintillator strips to WLS fibers and silicon photomultipliers, to design a compact electronics readout for several thousands channels, to set up the simulation and reconstruction algorithms and finally to carry out some commissioning measurements to test the overall performance of the setup.
The global performance of a muon tomograph for fast container inspection depends both on hardware aspects (efficiency, spatial and angular resolution of the detector, capability to operate in a stable mode over long time periods,…) and on software aspects (speed and performance of the reconstruction and analysis algorithms, graphics capabilities, easiness of calibration and modification of data collection conditions,..).
The space resolution in each detection plane (either X or Y) amounts to about 3 mm, due to the strip segmentation. The typical angular resolution in the reconstruction of the muon track above and below the inner volume was evaluated to be around 0.3 • for the geometry described in Sect. 2, also taking into account multiple scattering effects, as modelled by Critical aspects of the project were also evidenced, mainly related to the requirement of reaching optimal working conditions for the photosensors collecting the scintillation light. This in turn requires silicon photomultipliers with a high photon detection efficiency and small dark current, as well as to operate in a controlled light environment, with low temperature fluctuations over a large volume, which is now possible thanks to the availability of a variety of new SiPMs on the market.
The light collection efficiency is a key factor, especially in the present conditions, since a channel compression strategy was adopted, which required the sharing of the scintillation light from each strip between two WLS fibers, in order to reconstruct the hit strip by the proper combination of Tens and Units information, reducing the number of electronic readout channels. This solution however strongly reduces the overall detection efficiency, which is given to a first approximation by N p ×N l ×N f , where N p (= 2) is the number of physical planes (X or Y) in each logical plane, N l (= 4) is the number of logical planes, N f (= 2) the number of WLS fibers per strip and is the average individual efficiency for the light collection from each WLS fiber. While for the present set of measurements the detection efficiency of the individual channel was about 80%, the above consideration shows that even with an efficiency = 90%, the overall efficiency of the entire Portal would only be 18%. On the other hand, a small further increase in the individual light collection efficiency (say for instance from 90% to 95%) would increase by a factor 3 the overall efficiency. The expected improvement in case of individual readout of each strip, by summing up the light collected by the two fibers into a single channel, would be even larger. As a matter of fact we have to consider that in the present situation the number of light photons collected from each of the two WLS fibers in the same strip are not really independent quantities, due to the sharing of the photons between the two channels, with an additional reduction of the efficiency if each signal has to overcome the channel threshold.
Of course an improved efficiency would reduce by a corresponding factor the scan time required to have a reliable alarm information, which, in the above-mentioned conditions, was still too high (several hours) for most practical applications. After the completion of the Phase 1 of this Project, work is still in progress, to convey and sum the scintillation light from the two WLS fibers into a single, larger area (3 mm diameter) and small dark current SiPM. From preliminary tests, it is expected that this strategy, combined with the use of a correspondingly higher number of readout channels, would increase the efficiency by at least a factor 20, reducing the scan time to values of the order of 10-15 min, which are compatible with the requirements of large-scale applications.

Concluding remarks
In the present paper we addressed several aspects of the analysis of muon tomographic data which can be performed with the simplest strategy, based on the Point Of Closest Approach method. This basic algorithm, despite its intrinsic simplicity and limited validity, offers the possibility to generate a wide set of different mono-and two-dimensional distributions of the variables of interest, which may be used as multiparametric views of the 3D distribution of the scattering centers, to complement the information provided by each single projection or by a 3D view. The results are promising, since a body of compelling evidence for the presence of a peak in the various distributions, originated from the presence of a few dm 3 Pb block, could be easily obtained in this way. Various sets of cuts and weighting factors for the collected events have been investigated, comparing the signal-to-noise ratio and the significance obtained in different conditions. Due to its fast computing speed, the Poca method allows real-time analysis of the collected events. As discussed in Sect. 5, this also gives the opportunity to display at the same time, even on multiple screens or by the use of different processors, many views of the distributions, designing a customized presentation of the results, in order to facilitate the control of potentially hazardous materials inside the container being scanned. The need to have in real time a multitude of views and plots may be further satisfied by a parallel use of computing resources, since the different projections extracted in the present analysis may be worked out in parallel by different processors, accessing the data as they are measured, and updating the required projections in real time. This is also an important aspect when considering other advanced reconstruction and imaging algorithms, such as the Expectation Maximization-Maximum Likelihood (EM-EL), some variant of clustering algorithms and the Friends-Of-Friends method, on which we have been working [13,19]. These are known to be more CPU time consuming and usually require more powerful graphical tools, but provide better performance with respect to the basic Poca method, especially in case of multiple scattering centers.
Funding Open access funding provided by Università degli Studi di Catania within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.