Atmospheric ray tomography for low-Z materials: implementing new methods on a proof-of-concept tomograph
Gholamreza Anbarjafari, Aivo Anier, Egils Avots, Anzori Georgadze, Andi Hektor, Madis Kiisk, Marius Kutateladze, Tõnu Lepp, Märt Mägi, Vitali Pastsuk, Hannes Plinte, Sander Suurpere
AAtmospheric ray tomography for low-Z materials:implementing new methods on a proof-of-concept tomograph
Gholamreza Anbarjafari b,f , Aivo Anier a , Egils Avots a,b , Anzori Georgadze e,a,c , Andi Hektor a,d, ∗ , Madis Kiisk a,c ,Marius Kutateladze a , T˜onu Lepp a,c , M¨art M¨agi a , Vitali Pastsuk a,c , Hannes Plinte a , Sander Suurpere c a GScan OU, Maealuse 2 /
1, 12618 Tallinn, Estonia b iCV Lab, University of Tartu, Narva mnt 18, Tartu 51009, Estonia c Institute of Physics, University of Tartu, W. Ostwaldi tn 1, 50411 Tartu, Estonia d KBFI, Akadeemia tee 23, 12618 Tallinn, Estonia e Institute for Nuclear Research NANU, 03028 Kyiv, Ukraine f PwC Advisory, Helsinki, Finland
Abstract
Cosmic rays interacting with the atmosphere result in a flux of secondary particles including muons and electrons.Atmospheric ray tomography (ART) uses the muons and electrons for detecting objects and their composition. Thispaper presents new methods and a proof-of-concept tomography system developed for the ART of low-Z materials.We introduce the Particle Track Filtering (PTF) and Multi-Modality Tomographic Reconstruction (MMTR) meth-ods. Based on Geant4 models we optimized the tomography system, the parameters of PTF and MMTR. Based onplastic scintillating fiber arrays we achieved the spatial resolution 120 µ m and 1 mrad angular resolution in the trackreconstruction. We developed a novel edge detection method to separate the logical volumes of scanned object. Weshow its e ff ectiveness on single (e.g. water, aluminum) and double material (e.g. explosive RDX in flesh) objects.The tabletop tomograph we built showed excellent agreement between simulations and measurements. We are ableto increase the discriminating power of ART on low-Z materials significantly. This work opens up new routes for thecommercialization of ART tomography. Keywords: muon tomography, atmospheric ray tomography, plastic fiber scintillators, detection of low-Z materials
1. Introduction
Coulomb scattering based muon tomography wasproposed in 2003 [1]. Since then many research groupshave demonstrated that the method can be successfullyimplemented in di ff erent fields to detect high-Z ma-terials, e.g., nuclear safety [2–4] or security applica-tions [5–7]. In recent years several reviews have ap-peared providing excellent insights into last develop-ments [8–12].Coulomb scattering depends strongly on atomic num-ber Z . It makes the recognition of high-Z materials (e.g.nuclear materials) in muon tomography or more gener-ally, in atmospheric ray tomography (ART), easier thanthe recognition of low-Z materials. Thus ART applica-tions on low-Z materials have merged slowly. Techno-logical challenges are particularly complex in security ∗ Corresponding author.
Email address: [email protected] (Andi Hektor) applications, where low-Z materials should be identi-fied within minute-scale time.However, some dedicated research papers have ap-peared on the feasibility of the low-Z ART for securityapplications. Klimenko et al have studied whether elec-trons can be used for flux attenuation as an additionalsignature in muon tomography [13]. Other accompany-ing physical e ff ects such as muonic X-rays and electro-magnetic showers have also been investigated. The au-thors summarize that the inclusion of the latter requiresadditional detection systems. Cuellar et al reported re-sults based on Geant4 simulations that the stopping ef-fect of the low momentum part of the leptons spectrumcould be used as a source of additional information toscattering tomography, especially for medium and low-Z materials [14]. Later the e ff ect has been demonstratedwith simulated and experimental results [15]. The au-thors have proposed to exploit the ratio of stoppingpower to scattering, which enables to eliminate the sam- Preprint submitted to Nucl. Instrum. Methods Phys. Res., A March 2, 2021 a r X i v : . [ phy s i c s . i n s - d e t ] M a r le thickness as an unknown variable. First, one identi-fies the material from the ratio, then the mean scatteringangle and the known radiation length can be computedto calculate the thickness. The simulated results withmaterials of 1 m for heavy metal plates or cubes of1 m for paper and nylon, as well as experimental re-sults on the pallet-size units with low-Z materials, werepresented with 30, 10 and 5 minute long measurements.The detection of explosives and narcotics has been in-vestigated at the TUMUTY facility based on the Geant4simulations, where the material classification was doneusing the machine learning (ML), Support Vector Ma-chine (SVM) [16]. The authors conclude that ARTallows to discern narcotics and explosives from back-ground and metals, but not be separated from eachother (the various 20 × ×
20 cm size objects with10 to 30 minute measuring time were studied). Theexperimental validations of the numeric results on anRPC-based tomography system concluded that the ac-tual spatial resolution (0.6 mm) was lower than previ-ously were simulated [17]. Four materials were tested:flour (substitute of narcotics), aluminium, steel, andlead. The overall discrimination rate with the SVM clas-sifier could reach 70, 95, and 99% with 1, 5, and 10minutes measurement, accordingly.We have found only a few examples of low-Z ARTin non-security applications. Bikit et al demonstratedthat low-Z materials can be imaged by the simultane-ous detection of cosmic muons and photons created viabremsstrahlung inside the investigated bone and soft tis-sue sample [18, 19]. Another example is inspection ofnuclear waste encapsulated in concrete matrix. The in-vestigations on the containers of nuclear waste showedthat litre-size gas bubbles in concrete due to uraniumdioxide oxidation can be discovered within days of ex-posure [20]. Two recently papers reported that ART candetect reinforcement rebars in concrete. The first studyis based on Monte Carlo simulations [21] and the sec-ond on experimental tests [22].Three main factors influence the ART performance:(i) the accuracy of particle trajectory estimation, (ii)ability to determine particle type and energy, (iii) thetomographic reconstruction methods deployed. Thefourth, a system independent factor, is the measurementtime available. In the early days of muon tomography, ageometry-based reconstruction approach, Point of Clos-est Approach (PoCa) [23], was introduced. PoCa is ap-proximating multiple scattering processes with a singleinteraction point. It shows fairly good results in detec-tion of high-Z materials submerged to low-Z materials.In conventional X-ray Computed Tomography (CT) an-alytical methods such as filtered back projection (FPB) have been historically widely used for its computationallightness. In CT the positions of source and detector arewell known, but non-isotropic exposure and high Pois-son noise in ART reduce its performance considerably.Iterative methods seem to be most suited, either alge-braic [24] or statistical approaches [25, 26], where thebest trajectory estimates are obtained using scatteringand displacement data.Due to the properties of the Coulomb multiple scat-tering, the performance of tomographic reconstructionis greatly improved if the energy of passing particle isestimated. For example, the measured angular distri-bution of passing particles allows to classify the par-ticles according their estimated energy [27]. The pro-posed momentum multi-group model method was ex-perimentally tested and the reconstruction performancewas quantitatively compared against the squared aver-age, median angle squared and the average of squaredangle using the receiver operating characteristic (ROC)curves [28]. The results demonstrate clear superiorityfor multi-group approach. The CRIPT-project presentedan another approach of momentum estimation employ-ing massive iron sheets between their tracking detec-tors [29].The accurate reconstruction of particle tracks is espe-cially important in low-Z ART and in applications hav-ing sub-cubic-meter volumes scanned. The resolutionof scattering angle at mrad is required. The scatteringdensity of muons is directly related to the atomic num-ber Z of material. The density is an order of magnitudelower for materials with low-Z (Ca or lighter elements, Z (cid:46)
20) in comparison to high-Z (Pb or heavier ele-ments, Z (cid:38) µ m (FWHM) and angles of 2 mrad(FWHM) [5]. For large scale systems, a sub-mm posi-tion resolution perpendicular to the drift tube wire hasbeen reported [15]. For an another type, resistive platechambers (RPC), spatial resolution of 0.6 mm on lab-scale has been achieved [17]. The authors show no dataon the spacing between the detector plates, so we cannotestimate the angular resolution.On the sea level, the electrons and positrons consti-tute about 30-40% of the total lepton flux and abouta half of that is in practical usable energy range [13].Due to low energy, this fraction of the lepton spectrumis strongly influenced from the environment around thedetector system, for example, the building constructionstructures above the detector system. However, one canemploy it as an additional source of information for thedetection of small and medium sized objects, such as2assenger luggage or similar scale objects.In this paper we present new methods and a proof-of-concept tomography system developed specifically forthe ART of low-Z materials. We introduce the ParticleTrack Filtering (PTF) and Multi-Modality TomographicReconstruction (MMTR) methods. PTF allows to clas-sify the passing particle events (of charged leptons) intodi ff erent groups so that the absorption and the scatter-ing e ff ect are more isolated. Some aspects of the PTFmethod are described in a patent application by someauthors of the paper [30].In tomographic reconstruction we have developed anovel edge detection method to separate the logical vol-umes of scanned object. We demonstrated the e ff ective-ness of the method in case of single composition mate-rial (e.g. water, aluminum, steel) and double composi-tion material (e.g. explosive RDX in flesh-like material)objects. The MMTR approach enables to combine mul-tiple reconstructions obtained from the di ff erent PTF-classified particle event groups and the scattering andtransmission parameters from the track reconstruction.We have composed the numerical models using theGeant4 software package [31] to optimize the tomog-raphy system and the parameters of PTF and MMTRprocedures. For the tracker system (the hodoscope) wechose plastic scintillating fibers arrays. With this trackerwe were able to achieve a spatial resolution of 120 µ mfor muons and the 1 mrad angular resolution in the par-ticle track reconstruction.To validate our theoretical and numerical tests webuilt a physical proof-of-concept tomography system.The measurements on the prototype showed excellentagreement with the numerical simulations. We are ableto increases significantly the discriminating power ofART for low-Z materials.The structure of the paper is the following. We startwith the description of PTF realized on plastic scintil-lator fiber arrays. In Section 3 we present the MMTRmethod. In both sections we validate the methods withnumerical Geant4 models. In Section 4 we describethe general design principles of the physical proof-of-concept ART system. In Section 5 we give the technicaldescription of the prototype and the first results of themeasurements. In the last section we summarize theresults and give some hints on the next steps we haveconsidered.
2. The particle track classification: the PTF method
In this section we describe the main elements of theParticle Track Filtering (PTF) method. We have devel-oped a virtual hodoscope system in Geant4 addressing the critical factors of the low-Z ART described above.Mainly, we address the accuracy of particle track recon-struction and the ability to classify the type and energyof particle event candidates. The hodoscope consistsof the three position sensitive detector plates composedby plastic scintillator fibers. The detector plate has 4-layered structure: the two orthogonally placed double-layered fiber-mats, which give the x and y position of theparticle hit (see Fig. 1). The fiber diameter is 1.0 mm,the pitch in a single layer is 1.1 mm. The top layer isshifted by a half pitch, its position has been aligned ac-cording to the positions of the lower layer fibers. Thisensures close to 100% geometrical detection e ffi ciencyat the every angle of an incidence and high spatial reso-lution (for details, see Sec. 4). Figure 1: A snapshot from the Geant4 model of the position-sensitivedetector plate composed by the two orthogonally placed double-fiber-layers (gray cylinders). The red line denotes the trajectory of a passingmuon. The greened fibers illustrate the propagation of scintillationlight generated in the fibers by the muon passage.
The hodoscope has three detector plates in order toestimate the energy and the type of the passing particleusing the intrinsic scattering angle θ (in the hodoscope).Each detector plate has a well-defined thickness and ma-terial thus the probability distribution of the intrinsicscattering can be calculated. Though the typical θ ofpassing muons is in the order of mrad it can be mea-sured in the hodoscope. As a result we can classifythe reconstructed particle events (with some statisticalprobability) into separated categories according to theirscattering angle θ in the hodoscope. Fig. 2 illustrates theintrinsic scattering θ in the hodoscope.Fig. 3 shows the simulation results of the filteringspectrum for the hodoscope having the distance 100 mmbetween the two adjacent plates. We used the CRY cos-mic ray event generator to model the atmospheric rayflux consisting of muons and electrons at sea level [32].We fixed the spatial resolution of detector plates at3 igure 2: The intrinsic scattering angle θ we use to classify the particleevents (passing muons / electrons) in the hodoscope. The hodoscopehas the three detector plates (black bold lines numbered as 1, 2, 3).The arrows denote the reconstructed particle trajectory. The angle θ denotes the intrinsic scattering angle of the particle in the plate 2. ff erent number of groups; we callthose PTF groups below. For example, a possible ro-bust PTF classification schema is to classify the particleevents to muons and electrons or the muons with low,medium and high momentum. In Fig. 3 we have sepa-rated the spectrum into the three PTF groups: F1 (dom-inated by muons), F2 (mixed muons and electrons) andF3 (dominated by electrons). This classification schemahas been used in the figures of the reconstruction resultsin the next sections.
3. Tomographic reconstruction: the MMTR method
With the design of the hodoscopes as describedabove, we developed a model of the tomography systemin Geant4 with the interrogation volume 1 × × .The motivation of the size comes from practical applica-tions, e.g., it would be suitable for small luggage, postalparcels or human body. Fig. 4 shows the dimensions ofthe simulated tomography system. The system consistsof the two horizontal and the two vertical hodoscopesforming an enclosed interrogation volume aka the vol-ume of interest (VOI). When a muon / electron passes theentrance hodoscope it provides the entrance incidenceangle, the intrinsic scattering angle for PTF and the en-trance trajectory to the VOI with its entrance location.In the same way one can characterize an exiting particlefrom VOI. Figure 3: The distribution of atmospheric ray muons and electrons asa function of the intrinsic scattering angle θ in the hodoscope (fromthe Geant4 model with the CRY event generator). The distributionshows we can apply the intrinsic scattering angle θ as a discriminat-ing parameter classifying the type and energy range of the hodoscopepassing particle. The latter improves the tomographic reconstructionof scanned samples very significantly. The colored areas denote themuon and electron dominated values of θ (blue, yellow) and the mixedregion (gray). The geometry of the tomography system maximisesthe usage of side flux: on the y -axis the side hodoscopesmeasure the low incidence angle flux and on the x -axisthe elongated hodoscopes with the 2 m length increasethe field of view in the central part of VOI.In our Geant4 simulations the exposure time was es-timated by the internal code of CRY, which calculatesthe elapsed time on the generation of selected numberof particles (all types) per the defined generation area,10 ×
10 m in our case. The tomography system wascentred below the CRY source. The CRY generationarea was placed 30 cm from the upper hodoscope.In order to test PTF, we implement a tomographic re-construction procedure based on the back projection.Thus we fill VOI with a voxel grid. In the back-projection, we use the z -planes as a reference point andcalculate the corresponding x - and y -coordinates of theparticle on a particular z -plane. This process is repeatedfor all z -planes, up to filling the VOI, see Fig. 5. The z -plane is formed from the voxel grid point at the z -coordinate. The amount of the z -planes depends on thevoxel size set by the user. As we focus on low- Z ma-terials and moderate volumes we approximate the tra-jectories of the particles with a straight line betweenthe entering and the exiting locations of a particle (atthe boundary of VOI). Then we allocate so called scorevalue to each voxel when a trajectory passes through thevoxel (see Fig. 6).We can apply the back-projection procedure to the4 Figure 4: A schema of the Geant4 model of the tomography systemwith the four hodoscopes (upper and lower horizontal and left andright vertical bold dark lines), the radiation source (CRY, gray areaabove) and an example of studied objects in the center of the system(dark purple). The numbers 1-6 denote the horizontal detector plates. transmission and the scattering parameters. In the caseof transmission based reconstruction, the count of par-ticle trajectories in the every voxel can be used as thescore value. In the case of scattering based reconstruc-tion then one can choose the score value as the totalscattering angle or the scattering angle per unit lengthof the trajectory. The total scattering angle is defined asthe angle between vectors BC and DE, see Fig. 6. Alter-natively, other back-projection parameters are possible(see Sec. 5).Summing up all particle trajectories one forms thevolume density map (VDM) of VOI. We build manyVDMs, a dedicated VDM per PTF group. In the caseof scattering parameter, we calculate the median or themean scattering angle for each voxel forming the VDMin the scattering based reconstruction.It is relevant that one has di ff erent options how tocombine the detector hit data in the di ff erent plate de-tectors. Not all the reconstructed trajectory candidatesare complete, i.e. with hits in every plate detector ofthe entrance and the exiting hodoscope. One can havedi ff erent strategies how to deal with incomplete trajec-tories. For simplicity, we remove all the candidates thathave been registered only in the entering hodoscope. Ifa candidate has hits in every plate of the entering ho-doscope and at least a hit in a plate of the exiting ho-doscope (detector plates 1, 2, 3 and 4 or 5 or 6, seeFig. 6) it forms so-called incomplete trajectory candi-date. To apply the scattering based backprojection, atleast the two hits in the plates are needed in the entrance Figure 5: The z -planes (blue) filling the volume of interest (VOI) andthe voxel grid (black). The number 3 refers to the closest upper detec-tor plate to VOI and 4 refers to the closest lower detector plate to VOI(see the numbering of the plates in Fig. 4). and the exiting hodoscopes. We consider the completeand incomplete candidates in the reconstruction.After the considering trajectory candidates we applyPTF as described in Sec. 2. We remind that the PTFprocedure classifies the candidates into three or morePTF groups. After that we feed the groups to the opti-mised tomographic reconstruction procedure per group.For example, we can apply di ff erent back-projection pa-rameters on di ff erent groups.The third step of the tomographic reconstruction isthe composing of VDM. To reduce the amount of possi-ble false trajectories we define the upper and lower limitvalues for the total scattering angle. Thus so far our re-construction results the number of VDMs depending onthe set of the PTF classification groups, back projectionparameters etc.Fig. 7 and 8 shows the e ff ect of filtering parameterson the composing of VDM. The figures compare the tworeconstructions of VOI (VDMs) with and without thehodoscope and VOI filtering using scattering and trans-mission parameter, accordingly. The figures show theVDMs at z =
100 cm, where the plane crosses the inter-rogated object at the half height. The object is a cubeof water with the side length of 10 cm centered in themiddle of VOI. The estimated real measurement time isapproximately 18 minutes, the voxel size 1 × × .The upper panels present the reconstruction procedure5 igure 6: The real and reconstructed particle trajectories through VOIin the Geant4 simulated tomography system. The designating scoresto the voxels that the reconstructed trajectory passes through shown. with no filtering. The lower panel of Fig. 7 presentsthe reconstruction, where the PTF extracts the high-est momentum fraction (F1, see Fig. 3) dominated bymuons. The lower panel of Fig. 8 shows the reconstruc-tion results, where the PTF extracts the lowest momen-tum fraction (F3, see Fig. 3) of muons and the most ofthe electrons. We obtained both the figures from thesame Geant4 simulation.The figures demonstrate that despite the reduced fluxin the lower panels the filtering can significantly in-crease the contrast of VDMs. In those examples wehave separated the two physical e ff ects: the absorptionand scattering. Thus the di ff erent settings of the filteringparameters and ranges can potentially increase the dis-crimination capability for object detection and materialclassification. In this subsection we describe the algorithms we de-veloped to extract object shapes from VDMs. A cru-cial property of the algorithms is that they have to op-erate well in high Poisson noise conditions. We alsonotice that due do proprietary information the method-ology of the edge detection procedure is disclosed withsome generality.We start the object reconstruction on the 2D data. Analgorithm enhances the edges of objects from VDM it-eratively. After several iterations, the image sharpnessis enhanced to a level, where we apply an image thresh-old to extract the binary image, the final reconstructed
Figure 7: The e ff ect of hodoscope and VOI filtering. A cross-sectionalGeant4 snapshots on horizontal plane ( z = − × ×
10 cm ) without (upper) and with (lower) PTF. The ho-doscope filtering range 0 < θ < shape of the object. Fig. 9 shows a double-object ex-ample, a smaller inner cube of RDX-explosive with theside length of 10 cm surrounded with a larger cube offlesh-like material with the side length 30 cm. As de-scribed above we apply PTF. In this case the VDM rep-resents the distribution of mean scattering angles. Thespatial resolution of the plates in the Geant4 simulationwas 0.1 mm. The estimated measurement time in theGeant simulation is 18 minutes.After the 2D reconstruction, the results are fed to the3D reconstruction procedure. To di ff erentiate betweenthe object candidates, we introduce a labelling step thatgives a unique index to the connected volumes. This canbe achieved having so-called Connected-ComponentLabelling (CCL), where the subsets of connected com-ponents are uniquely labelled using, e.g., Two Pass Al-gorithm [33]. Fig. 10 presents a half-cut of recon-structed and detected objects (for the same simulationas in Fig. 9). The result indicates clearly that the VDMsof the two materials are distinguishable enough for theedge detection and the two di ff erent logical volumeshave been established. The two di ff erent logical vol-umes are denoted: yellow for RDX, orange and blackedges for flesh.We tested the same object configuration witha shorter, 2-minutes estimated, measurement time.Fig. 11 shows the results. One can clearly see the loss of6 igure 8: The e ff ect of hodoscope and VOI filtering. A cross-sectionalGeant4 snapshots on horizontal plane ( z = − × ×
10 cm ) without (upper) and with (lower) PTF. In the lowerpanel, the hodoscope filtering range is 10 < θ <
30 mrad. contrast. Due to lower statistics, we increased the voxelsize from 1 cm to 9 cm . However, the edge detec-tion algorithm is still able to detect and separate the twomaterials as di ff erent logical volumes.Once we establish and separate the logical volumesof di ff erent object we can apply special algorithms forthe material or other classifications on the reconstructedvolumes. Those will be the subject of future work.
4. Design of the physical prototype
As described above we developed the PTF andMMTR methods based on the virtual Geant4 models.To prove (and improve) the methods we built a sim-plified physical prototype system. Fig. 12 shows theprincipal design of the physical prototype. The proto-type has an upper hodoscope with three detector platesbased on the above described fiber-mat technology. Thelower hodoscope has been replaced with a single detec-tor plate. This design means the system cannot measurethe total scattering angle. However, one can estimate thetotal scattering angle by calculating the angle betweenthe vectors BC and CD, see Fig. 6. This angle dependson the position of the object in a vertical axis, thus onecan do a comparison material classification experiment
Figure 9: Example of the double-object detection algorithm. From topto down, the original reconstructed image down to the detected objectsthrough the iterative object detection procedure in case of double-object system: a cube (10 × ×
10 cm ) of RDX explosive materialsurrounded by a larger cube of flesh (30 × ×
30 cm ). The esti-mated measurement time in the Geant simulation is 18 minutes. Thevoxel size is 1 × × , the spatial resolution of the detector plates0.1 mm. fixing objects into the same position in every measure-ment. An another shortcoming of the minimal designis that the detector plates cover only a small fractionof the total solid angle around VOI and only a limitedfraction of the flux is available. The detector plate sizeis 25 ×
25 cm , the distance between the two adjacentplated in the upper hodoscope is 7.5 cm. The volumeof VOI is 25 × ×
25 cm . Thus the field of view islimited to 64 degrees in x - and y -axis. The limited fieldof view results in low vertical reconstruction accuracy ifno correction algorithm implemented.First, we studied the detection and reconstruction per-formance of the prototype in the Geant4 simulations.We describe the main detection performance parame-ters, detection e ffi ciency and spatial resolution, in thefollowing subsections. In this section we explain how we build up a particletrack candidate form hit candidates (lit fibers). First, weexplain our notations and the coordinate system. For7 igure 10: The object reconstruction and detection of a double ob-ject with spectroscopic filtering relying on the scattering density pa-rameter. The double object: a cube of RDX explosive material(10 × ×
10 cm ) centered in a cube of flesh (30 × ×
30 cm ),the voxel size is 1 × × . The object detection and separation tothe two logical volumes are denoted by yellowish colors. The spatialresolution of the detector plate is 0.1 mm. simplicity, we deal with the x - and y -directions sepa-rately and so below we show only a direction if su ffi -cient. Fig. 13 shows the fiber schematics and the no-tations through the four plates of the prototype in a se-lected direction ( x or y ). As described in Sec. 2 andshown in Fig. 1 a plate has the four fiber layers in totalarranged to the two orthogonal bi-layers. We addressthe fiber layers as h AB where A denotes the plate and B denotes the layer in the plate.For example, h11 denotes the plate 1, the y-directional bi-layer and the first layer. h14 denotes theplate 1, the x-directional bi-layer and the second layer.For the z-coordinate of the hit candidate we define a vir-tual plane between the bi-layers. So the layers h11 andh12 are above the virtual z -plane and the layers h13 andh14 are below the z -plane. We define the component of x or y layers as the fiber-centre o ff sets from the z -plane.The following parameters describe the fiber: the corediameter, the width and the pitch of the cladding. Thepitch is the distance between the adjacent fiber centres.As we know all the positions and other parameters ofthe fibers we can compute the centres of the fibers inthe well-defined ’absolute’ Cartesian coordinate space.We enumerate the fibers starting from 0 to N in a layer,where N is the max number of fibers in the layer. Thusdetecting a fiber lit we can then convert it the ’absolute’coordinates.To reconstruct a plate hit we handle the x - and y -directions separately. For example, Fig. 14 shows the Figure 11: The object reconstruction and detection of a double ob-ject with spectroscopic filtering relying on the scattering density pa-rameter. The double object: a cube of RDX explosive material(10 × ×
10 cm ) centered in a cube of flesh (30 × ×
30 cm ), thevoxel size is 3 × × . The object detection and separation to thetwo logical volumes are denoted by yellowish colors. The estimatedmeasurement time 2 minutes and the spatial resolution of the detectorplate 0.1 mm. reconstructed hit for the y -direction from the three litfibers in the y by-layer. We compute the hit coordinatefinding the weighed arithmetic mean of fiber centres lit.In the figure, red denotes the fibers lit and the blackcross is the hit coordinate reconstructed. The same cal-culation is performed for the x by-layer. The red lineis an example of particle track candidate reconstructedhaving some other hit candidates from the other plates.The x and y coordinates of a hit reconstructed (atleast) in the two planes we can project the results tothe z -plane. We use the linear extrapolation based onthe reconstructed track candidate. Fig. 15 shows the ex-trapolation results (black circled crosses). Combiningthe component values from the both axes projected tosame z plane, we have reconstructed the 3-dimensionalcoordinate of a hit candidate.For the track reconstruction we need a clustering al-gorithm to cluster the lit fibers from the weighted meanof fiber centres computed. Fig. 16 shows some casesof clustered lit fibers to illustrate the approach. Cluster-ing algorithm works in the two stages. First, we clusterthe lit fibers from the single layers. We form clusters ofconsecutively lit fibers. In the layer h11 the fibers are1, 3, 4 and in the layer h12 the fibers 0, 2, 5. Second,we merge the clusters from the di ff erent layers if theyintersect component-wise: either the x or y coordinates8 igure 12: The general spatial design with the detector plates and thestructural elements. An example object (gray cubic) is in the middle ofthe detector between the upper hodoscope (the three horizontal plates)and the lower single detector plate (a plate below the cube). The fiber-mat couplings with SiPMs nor the DAQ-system are not shown. depending on the axis are handled. Fig. 16 shows thisas the dotted intersecting region. We compute the min-imum extreme of a cluster by taking the leftmost fiber,computing its intersections in the fiber layer z -plane andtaking the minimum value. The same logic is followedto calculate the maximum extreme value. We tested the hit and track reconstruction algorithmsdescribed in the last sub-section in the Geant4 simula-tions of the prototype. To estimate the quality and ef-fectiveness of the algorithms we collected the follow-ing information of the muon / electron hits in the virtualdetector plates: (i) the hit coordinates and (ii) the num-ber of fiber(s) providing the signal from the same eventfor the every detector plate and layer. We generated thesimulation data with a statistically su ffi cient number ofevents having the CRY-based source described above.Then we calculated the distance di ff erences between theGeant4 provided coordinate values (i) and the coordi-nates provided by our algorithms calculated for every x and y layer and all the simulated events. We calculatedthe distributions and the standard deviations, some re-sults are presented in Table 1. The numbers in the tableare based on the 13 699 simulated hit events: 10 907muon events and 2792 electron events. The distributionhas its mean at zero and it corresponds to the Laplacedistribution. Figure 13: The labeling schema of the bi-layered scintillating fibersof the detector (showing the y -directional by-layer only).Figure 14: The reconstruction of the hit candidate for the y -directionalby-layer. The red denotes the lit fibers. The cross depicts the recon-structed hit candidate. The red line is an example of particle track can-didate reconstructed having some hit candidates from the other plates. Stdev ( µ m) Detector plate Muons Electrons All
Plate 1 113 1234 567Plate 2 110 726 343Plate 3 153 840 1112Plate 4 109 2769 1261
Average 121 1392 821
Table 1: Standard deviations between the hit locations and recon-structed hit locations from the Geant4 simulations of the prototypetomography system.
In the case of muons, the standard deviation variesslightly above 100 µ m. The reason of this rather largevalue is in a limited set of extreme cases, which have1-2 mm deviations from the actual hit coordinates. Inthe case of electrons, we expect larger deviations dueto more scattering events of electrons in the detectorplates, decreasing the performance of the track recon-struction algorithm. However, the number of strongly9 igure 15: The hit extrapolation to the plate-z. The red fibers representfibers generating signal, the cross depicts the calculated particle hitcoordinate, the circled cross indicates the extrapolated coordinate.Figure 16: An clustering example, the red fibers represent the fibersgenerating signal. deflecting extreme event is considerably small. If weremove some tens of most deflecting events, the stan-dard deviation falls down to a few hundreds of µ m.To test the performance of the tomographic recon-struction in the Geant4 simulations we had di ff erentcube shaped materials in VOI. We fixed the locationof the object in the same position size and had the6 × × cube shaped objects. We centered theobject in the horizontal plane and located the objectsvertically 20 mm from the lower detector plate of theupper hodoscope. Fig. 17 and 18 show the reconstruc-tion results of the three di ff erent materials. To createthe VDMs we applied the mean scattering angle. Theestimated real measurement time was 30 minutes.The horizontal cross-sections have been recon-structed quite accurately, for lead object it is overesti-mated about 10 mm, for Plexiglas 10 mm less than itsactual size. Due to the limited angle tomography, recon-structed objects are vertically elongated about 2-3 timestheir actual size.
5. The physical prototype and the first measure-ments
Based on the design of the prototype tested withGeant4, we constructed the physical prototype to havea physical validation of the proposed detector systems,PTF and MMTR. As in the case of virtual model,the physical prototype (see Fig. 19) has the upper ho-doscope with three detector plates and the lower ho-doscope with single detector plate. We have the detector
Figure 17: The reconstructed logical volumes by the prototype tomog-raphy system. The scanned object (plexiglass, aluminium, iron) hasthe size 6 × × . The voxel size 0 . × . × . and theestimated measurement time 30 minutes. The scale of the coordinatesis denoted in voxels (e.g. 10 corresponds to 5 cm). layers with the 1 mm round single cladding scintilla-tion fibers (Saint-Gobain BCF-12) assembled as shownin Fig. 1. The active area of each detector plate is247 ×
247 mm . We spaced the detector plates in the up-per hodoscope with 75 mm pitch, the distance betweenlower plate of the upper hodoscope and the lower ho-doscope is 250 mm and thus the dimensions of the VOIare 247 × ×
250 mm .We installed the detector plates into a rigid innerframe that fixed the aluminium support frames of thedetector plates. The external frame is used to fix the po-sitions of the SiPM arrays and the DAQ boards. Ontop of the upper hodoscope and below the lower ho-doscope, we assembled the two trigger detectors (scin-tillator plates).We assembled the detector plates using groovedalignment tables (see Fig. 20). This guarantees the accu-rate and controlled positioning of fibers in the fiber mat.The alignment tables have been produced from 1 mmthick sheets of phenolic paper laminate (PFCP201) bymilling the grooves with a pitch of 1.100 mm. We gluedthe fibers on the alignment table using UV-curing glue(Loctite AA 3311). We added TiO2 powder to the glueto reduce optical leakage between the fibers.We assembled the two separately produced double-layered fiber mats orthogonally to each other in orderto form a 2D-sensitive four-layered detector plate. Wesupported the assembly with an aluminium frame thatprovides su ffi cient mechanical rigidity, flatness and fix-ation with the frame. At the corners of each PFCP20110 igure 18: The horizontal cross-sections of reconstructed logical vol-umes at the half-height of the object applying (the same reconstructionschema as in Fig. 17). The scanned object (plexiglass, aluminium,iron) has the size 6 × × . The voxel size 0 . × . × . and the estimated measurement time 30 minutes. The scale of thecoordinates is expressed in voxels (e.g. 10 corresponds to 5 cm). plate we have the alignment holes allow to align thetwo orthogonal 2D fiber mats against each other and toalign the four detector plates with respect to each. Forthe latter we used the alignment rods once the detectorplates were mounted in the frame. The structure en-sures a mechanical tolerance for relative plate positionsbelow ≤ . We composed the data acquisition system of the pro-totype implementing eight CAEN DT5550W boards.We use Ketek (PA3325-WB-0808) and Hamamatsu(S13361-3050AE-08) SiPMs arrays to collect scintilla-tion light from the fibres.Two plastic scintillator plates trigger the data read-out. The active area of a trigger is 30 ×
30 cm , theupper one placed on top of the upper hodoscope and thesecond one placed under the lower hodoscope. The trig-gers have two 6 × SiPMs (Hamamatsu S13360-6050CS) using the SiPM Readout Front-End BoardCAEN A1702. When a muon / electron passes the trig-ger system (both plastic scintillator plates) a trigger sig-nal is created that initiates the readout of all the CAENDT5550W boards. Figure 19: The prototype tomography system is encased into a light-tight box. The upper casing plate removed for the photograph.Figure 20: The assembly of a detector plate. The inset shows themagnified front image of the grooved surface of the support plate andthe fiber positions in the fiber mat.
Having the physical prototype our first objective wasto test the detection e ffi ciency of the completed particletracks. The e ffi ciency is defined as a completed particletrack candidate reconstructed if the triggered signals arecollected from all the 8 bi-layer fiber mats. The averagecount rate from the triggering system is about 3 Hz in-door (the lab room has concrete walls and ceilings andone additional floor above). Considering the triggeringdetector area (0.09 m ) and the field of view (66 de-grees) of the double trigger plate telescope, the averagecount rate is in the expected range. We performed a30-minute measurement with the empty VOI volume toestimate the e ffi ciency of the completed tracks. In thetest we recorded the 6594 triggered events, out of whichthe 3562 were completed tracks. It means the total e ffi -ciency for the completed particle tracks is around 50%.Fig. 21 shows the comparisons of the reconstruc-11ion having the di ff erent filtering conditions on the 2500completed trajectories. The voxel size was 0 . × . × . , the reconstruction parameter was the meanscattering angle. Due to the limited field of view, thehorizontal cross-sectional images are a little unclear,di ff erent from the simulation results of the full model ofthe tomography system. After we applied PTF, the con-trast of the reconstructed images improved, especiallyin the case of light materials. Figure 21: The horizontal cross-sections with the di ff erent filteringconditions for a cube of water and aluminium (6 × × ): nohodoscope nor VOI filtering applied (left); no hodoscope filtering ap-plied (middle); hodoscope and VOI filtering applied (right). Having the optimised PTF and MMTR parameters,we reconstruct the 3D images for the three di ff erentcases: a cubic of aluminium, Plexiglas and ammo-nium nitrate (in granular form) with the dimensions6 × × (see Fig. 22 and 23. The reconstructionparameter for VDMs was the mean scattering angle.The reconstruction results are vertically elongated aswe predicted from the Geant4 simulations and vary be-tween 15-25 cm. Horizontally, the side length vary be-tween 7-8 cm. The dimensions of the objects are overes-timated, which is expected due to limited angle tomog-raphy. Theoretically, one can apply a special correctionif needed.
6. Discussions and conclusions
In this paper, we discussed on the two new methodsin muon tomography: PTF (the classification of particlecandidates in the hodoscope) and MMTR (the creatingof multi-VDMs following the classification schema ofPTF, etc). Our main motivation is to improve object re-construction and material classification for low-Z mate-rials in the case of ART. We applied the methods on thehodoscope system based on plastic scintillating fibers.
Figure 22: The detected logical volumes as a result of the reconstruc-tion of a cube samples (6 × × ) of plexiglass (upper left), alu-minium (upper right) and steel (lower left). The samples were centred2 cm below the lower detector plate of the upper hodoscope. The redregion represents the logical volume for the detected object. The scaleof the coordinates is in cm. We proposed the detector technology to ensure the re-quired detector performance, mainly the spatial accu-racy, the sensitivity and the lightness of the detector sys-tem. First, we presented the results based on the Geant4simulations.Second, we validated the methods on the minimaltable-top prototype tomography system. We tested theperformance of PTF and MMTR on the prototype hav-ing di ff erent low-Z materials. We showed that apply-ing the proposed methods we can successfully detect thelow-Z objects in a 10–30-minute time scale. For exam-ple, we tested the objects composed by water, Plexiglas,ammonium nitrate (mimicking explosive material). OurGeant4 simulations indicate that one can reduce the ex-posure time to a few minutes range having the full-scaletomography system (see Sec. 3) based on the same tech-nology like the physical minimal prototype.The particle tracking algorithm determining the par-ticle hit location in a detector plate and its trajectory inthe hodoscope provide the accuracy 120 µ m for muons.We demonstrated that VDM can be analysed as a set of2D images applying an edge detection algorithm on theVDM. It allows to detect the objects as separated logi-cal volumes. The method works e ffi ciently in the caseof high Poisson noise, hence applicable for short mea-surement times.When logical volumes reconstructed one can be inter-ested in the material classification of volumes. The pos-12 igure 23: The horizontal cross-sectional cuts of the same samples asin Figure 22. The cut is at the half-height of the object and the recon-struction schema same as in Figure 22. The scale of the coordinates isin cm. sible classification procedure can follow as a single pro-cedure or having the multiple procedures, e.g., one perPTF group. For example, if one combines the classifica-tion based the filtered fractions F1 and F3 (see Fig. 1),the discrimination ability of materials increases includ-ing both the scattering and transmission data simultane-ously. Furthermore, establishing the reconstructed ob-jects as discrete entities enables to apply material clas-sification techniques based on machine learning.One of our next steps is to develop the material clas-sification algorithms on VDMs. Some potential meth-ods to start the classification study would be machinelearning based such as the Kullback-Leibler Divergence(KLD) or the Support vector Machine (SVM). Anotherfuture step will be the reconstruction of complex objectsrepresenting more realistic user cases. For example, wecan implement some iterative reconstruction methodsin order to increase the ability of generating VDMs forcomplex scanning objects.The PTF and MMTR methods combined with theproposed detector system improve the reconstruction ofVOI significantly. The combination makes the detectionand potentially classification of low-Z composed mate-rials possible in small-scale detector system. It opens upnew routes for the commercialization of ART tomogra-phy. The developed methods and the detector systemis interesting in the context of industrial security appli-cations such as baggage scanning in airports and cargoand shipping container scanners for illicit materials andcontraband. Acknowledgements
We would like to thank Andrea Giammanco for veryhelpful discussions and comments. This work was con-tinuously and strongly supported by the startup com-pany GScan O ¨U. This work was supported by the En-terprise Estonia, the grant EU48693. AH thanks the Es-tonian Research Council for the grant PRG434 and theEC for the ERDF CoE program project TK133.
References [1] K. N. Borozdin, G. E. Hogan, C. Morris, W. C. Priedhorsky,A. Saunders, L. J. Schultz, M. E. Teasdale, Radiographicimaging with cosmic-ray muons, Nature 422 (2003) 277–277.doi: .[2] A. Clarkson, et al., Characterising encapsulated nu-clear waste using cosmic-ray muon tomography, JINST10 (2015) P03020. doi: . arXiv:1410.7192 .[3] D. Poulson, J. Durham, E. Guardincerri, C. Morris, J. Bacon,K. Plaud-Ramos, D. Morley, A. Hecht, Cosmic ray muon com-puted tomography of spent nuclear fuel in dry storage casks,Nucl. Instrum. Meth. A 842 (2017) 48–53. doi: . arXiv:1604.08938 .[4] D. Mahon, A. Clarkson, S. Gardner, D. Ireland, R. Jebali,R. Kaiser, M. Ryan, C. Shearer, G. Yang, First-of-a-kind muog-raphy for nuclear waste characterization, Phil. Trans. R. Soc. A.377 (2019). doi: .[5] C. L. Morris, et al., Tomographic imaging with cosmic raymuons, Sci. Glob. Secur. 16 (2008) 37–53. doi: .[6] F. Riggi, et al., The muon portal project: Commissioning of thefull detector and first results, Nuclear Instruments and Meth-ods in Physics Research Section A: Accelerators, Spectrome-ters, Detectors and Associated Equipment 912 (2018) 16–19.doi: .[7] A. Harel, D. Yaish, Lingacom muography, Philosophical Trans-actions of the Royal Society A: Mathematical, Physical and En-gineering Sciences 377 (2019) 20180133. doi: .[8] L. Bonechi, R. D’Alessandro, A. Giammanco, Atmosphericmuons as an imaging tool, Rev. Phys. 5 (2020) 100038.doi: . arXiv:1906.03934 .[9] R. Kaiser, Muography: overview and future directions, Phil.Trans. R. Soc. A. 377 (2019). doi: .[10] S. Vanini, P. Calvini, P. Checchia, A. Rigoni Garola, J. Klinger,G. Zumerle, G. Bonomi, A. Donzella, A. Zenoni, Muographyof di ff erent structures using muon scattering and absorption al-gorithms, Philosophical Transactions of the Royal Society A:Mathematical, Physical and Engineering Sciences 377 (2019)20180051. doi: .[11] S. Procureur, Muon imaging: Principles, technologies and ap-plications, Nuclear Instruments and Methods in Physics Re-search Section A: Accelerators, Spectrometers, Detectors andAssociated Equipment 878 (2018) 169–179. doi: , radiation Imaging Techniques and Appli-cations.[12] G. Bonomi, Progress in Muon Tomography, PoS EPS-HEP2017(2017) 609. doi: .[13] A. V. Klimenko, C. L. Morris, W. C. Priedhorsky, K. N.Borozdin, L. J. Schultz, Fusing signatures of di ff erent physical rocesses in muon tomography, in: IEEE Nuclear Science Sym-posium Conference Record, 2005, volume 1, 2005, pp. 307–311. doi: .[14] L. Cu´ellar, K. N. Borozdin, J. A. Green, N. W. Hengartner,C. Morris, L. J. Schultz, Soft cosmic ray tomography for de-tection of explosives, in: 2009 IEEE Nuclear Science Sym-posium Conference Record (NSS / MIC), 2009, pp. 968–970.doi: .[15] G. Blanpied, S. Kumar, D. Dorroh, C. Morgan, I. Blanpied,M. Sossong, S. McKenney, B. Nelson, Material discrimina-tion using scattering and stopping of cosmic ray muons andelectrons: Di ff erentiating heavier from lighter metals as well aslow-atomic weight materials, Nuclear Instruments and Meth-ods in Physics Research Section A: Accelerators, Spectrom-eters, Detectors and Associated Equipment 784 (2015) 352–358. doi: , symposium on Ra-diation Measurements and Applications 2014 (SORMA XV).[16] Z. Yifan, Z. Zhi, Z. Ming, W. Xuewu, Z. Ziran, Discrimi-nation of drugs and explosives in cargo inspections by apply-ing machine learning in muon tomography, High Pow. LaserPart. Beams 30 (2018) 086002. doi: .[17] X.-Y. Pan, Y.-F. Zheng, Z. Zeng, X.-W. Wang, J.-P. Cheng,Experimental validation of material discrimination ability ofmuon scattering tomography at the tumuty facility, Nu-clear Science and Techniques 30 (2019) 120. doi: .[18] I. Bikit, D. Mrdja, K. Bikit, J. Slivka, N. Jovancevic, L. Ol´ah,G. Hamar, D. Varga, Novel approach to imaging by cosmic-raymuons, EPL (Europhysics Letters) 113 (2016) 58001.[19] D. Mrdja, I. Bikit, K. Bikit, J. Slivka, J. Hansman, L. Ol´ah,D. Varga, First cosmic-ray images of bone and soft tissue,EPL (Europhysics Letters) 116 (2016) 48003. doi: .[20] M. Dobrowolska, J. Velthuis, L. Fraz˜ao, D. Kikoła, A noveltechnique for finding gas bubbles in the nuclear waste contain-ers using muon scattering tomography, J. Instrum. 13 (2018)P05015–P05015. doi: .[21] M. Dobrowolska, J. Velthuis, A. Kopp, M. Perry, P. Pearson,Towards an application of muon scattering tomography as atechnique for detecting rebars in concrete, Smart Materialsand Structures 29 (2020) 055015. doi: .[22] E. Niederleithinger, S. Gardner, T. Kind, R. Kaiser, M. Grun-wald, G. Yang, B. Redmer, A. Waske, F. Mielentz, U. E ff ner,et al., Muon tomography of a reinforced concrete block–first ex-perimental proof of concept, arXiv preprint arXiv:2008.07251(2020).[23] L. Schultz, K. Borozdin, J. Gomez, G. Hogan, J. McGill,C. Morris, W. Priedhorsky, A. Saunders, M. Teasdale, Image re-construction and material z discrimination via cosmic ray muonradiography, Nuclear Instruments and Methods in Physics Re-search Section A: Accelerators, Spectrometers, Detectors andAssociated Equipment 519 (2004) 687 – 694. doi: https://doi.org/10.1016/j.nima.2003.11.035 .[24] Z. Liu, S. Chatzidakis, J. M. Scaglione, C. Liao, H. Yang, J. P.Hayward, Muon tracing and image reconstruction algorithmsfor cosmic ray muon computed tomography, IEEE Transactionson Image Processing 28 (2019) 426–435. doi: .[25] L. J. Schultz, G. S. Blanpied, K. N. Borozdin, A. M. Fraser,N. W. Hengartner, A. V. Klimenko, C. L. Morris, C. Orum, M. J.Sossong, Statistical reconstruction for cosmic ray muon tomog-raphy, IEEE Transactions on Image Processing 16 (2007) 1985–1993. doi: . [26] S. Pesente, S. Vanini, M. Benettoni, G. Bonomi, P. Calvini,P. Checchia, E. Conti, F. Gonella, G. Nebbia, S. Squarcia, G. Vi-esti, A. Zenoni, G. Zumerle, First results on material identifica-tion and imaging with a large-volume muon tomography pro-totype, Nuclear Instruments and Methods in Physics ResearchSection A: Accelerators, Spectrometers, Detectors and Associ-ated Equipment 604 (2009) 738 – 746. doi: https://doi.org/10.1016/j.nima.2009.03.017 .[27] C. L. Morris, K. Borozdin, J. Bacon, E. Chen, Z. Luki´c,E. Milner, H. Miyadera, J. Perry, D. Schwellenbach, D. Aberle,W. Dreesen, J. A. Green, G. G. McDu ff , K. Nagamine, M. Sos-song, C. Spore, N. Toleman, Obtaining material identificationwith cosmic ray radiography, AIP Advances 2 (2012) 042128.doi: .[28] J. O. Perry, J. D. Bacon, K. N. Borozdin, J. M. Fabritius, C. L.Morris, Analysis of the multigroup model for muon tomographybased threat detection, Journal of Applied Physics 115 (2014)064904. doi: .[29] V. Anghel, J. Armitage, F. Baig, K. Boniface, K. Boudjem-line, J. Bueno, E. Charles, P.-L. Drouin, A. Erlandson, G. Gal-lant, R. Gazit, D. Godin, V. Golovko, C. Howard, R. Hydo-mako, C. Jewett, G. Jonkmans, Z. Liu, A. Robichaud, T. Stocki,M. Thompson, D. Waller, A plastic scintillator-based muon to-mography system with an integrated muon spectrometer, Nu-clear Instruments and Methods in Physics Research Section A:Accelerators, Spectrometers, Detectors and Associated Equip-ment 798 (2015) 12 – 23. doi: .[30] A. Georgadze, M. Kiisk, M. M¨agi, E. Avots, G. Anbarjafari,Method and apparatus for detection and / or identification of ma-terials and of articles using charged particles, 2020.[31] J. Allison, et al., Geant4 developments and applications, IEEETrans. Nucl. Sci. 53 (2006) 270. doi: .[32] C. Hagmann, D. Lange, D. Wright, Cosmic-ray shower gener-ator (cry) for monte carlo transport codes, in: 2007 IEEE Nu-clear Science Symposium Conference Record, volume 2, 2007,pp. 1143–1146. doi: ..