A muon tracker was developed using three polyvinyl toluene scintillator panels instrumented with photomultiplier tubes (PMTs) mounted at the corners. Panels are mounted in parallel on an aluminum frame which allows for simple adjustment of angle, orientation and separation distance between the panels. The responses of all PMTs in the system are digitized simultaneously at sub-nanosecond sample spacing. Software was developed to adjust settings and implement event rejection based on the number of panels that detected a scintillation event within a 400-nanosecond record.  The relative responses of the PMTs are used to calculate the position of scintillation events within each panel. The direction of the muons through the system can be tracked using the panel strike order. Methods for triangulation by both time-of-flight (TOF) and PMT magnitude response are reported. The time triangulation method is derived and experimentally demonstrated using parallel cables of differing length. The PMTs used in this experiment are only optimized for amplitude discrimination, not for time spread jitter as would be required to implement TOF methods into the scintillator panels. A Gaussian process regression machine learning tool was implemented to learn the relationship between PMT response features and positions from a calibration dataset. Resolution is analyzed using different numbers of PMTs and low-versus-high PMT sensitivities.  Muons traveling in forward and reverse directions through the detector system were counted in all six axis orientations. The muon detector was deployed for 28 days in an underground tunnel and vertical muon counts were recorded.


Muons are elementary particles commonly created by cosmic rays interacting with the upper atmosphere and neutrino interactions within the earth [1-4]. High-energy muons can travel through the atmosphere and penetrate several kilometers into the crust of the Earth. While these particles weakly interact with matter, their lifetime and deflection angles are sensitive to path length and material density, making them suitable for some imaging applications. Density variations of material between the muon source and detector can be measured from the muon flux rate at different locations, in many ways like medical x-ray or computed tomography scans[5-11]. Muons were first used for geologic characterization to measure the overburden of a tunnel [12], which provided a faster and more cost-effective method than drilling for depths greater than 15.24 m (50 ft). Alvarez et al. [13] used muons to search for hidden chambers within the Egyptian pyramids. Muon imaging techniques have also been applied to the investigation of volcanoes [14,15]. In recent years, research in using muons for noninvasive geophysical imaging applications has expanded [9,16-22]. Muon imaging techniques have also been applied to material investigation of objects [8,11,23,24], nuclear fuel verification [25], searching for nuclear threats [26], and transportation monitoring [27-29] .

Much of the work in muon tomography has been done using muons traveling at low zenith angles (close to vertical) where contamination from muons traveling from the direction opposite the target are negligible. For this reason, maintaining low zenith angles simplifies instrumentation timing needs. If, however, the detector is mounted horizontally (zenith angle = 90 degrees), muons will pass through the detector from both directions, contaminating any image which assumes unidirectional paths. With a sub-nanosecond sample rate, precise cross-channel synchronization in the digitizer, length-matched coaxial cables, and careful construction/calibration, the system described here can determine the forward/reverse directionality of particles traveling through the panels. This capability eliminates image contamination from the opposite direction and potentially enables simultaneous imaging from both directions.

One of the most common architectures of muon path detectors employs strips of scintillator material [17,26,30-33], arranged in orthogonal rows with a dedicated photomultiplier tube (PMT) monitoring scintillation events, or “strikes”, on each strip. In these “hodoscope” systems, a planar strike coordinate registers when a set of PMTs respond at the “same” time. This method enables position tracking of muons passing through the detector plane, with resolution limited by the size and number of strip/PMT pairs employed. To avoid space congestion, light produced from an event can be guided from scintillator to PMT using fiber optic lines. Similarly, scintillator strips or solid panels have been used with an installed pattern of wavelength shifting (WLS) fibers on the surface [34-36]. The fibers absorb light from the scintillator and re-emit light in a different wavelength which is monitored with PMTs. Each of these systems has been demonstrated successfully but they are often mechanically complex and therefore costly, especially when resolution and path tracking improvements require substantively more scintillators and/or PMTs.

Another proven hodoscope-type muon tracking system uses a grid of drift chamber tube detectors. Drift chambers consist of wires tensioned in the center of a cylindrical chamber that is filled with specialized gas and subject to a uniform high voltage electric field. Gas ionization cascades caused by passing muons can be detected with coupled instrumentation on the wires. The distance from the tangential muon path to the central wire location can be determined if an accurate drift time calibration is available. In drift tube tracking systems, muon tracking resolution can be improved by determining the best straight line fit to the possible tracks from the activation pattern of a given event [23]. Drift chamber muon tomography systems can be very accurate but are also heavy, require many channels, and must be carefully assembled and maintained to prevent failure due to gas leakage.

A recent simulation study [37] suggests the feasibility of muon detection systems based on the concept of scintillation cameras [38] using large sheets of plastic scintillator material with a low number of PMTs. For a given resolution, these systems can be a great deal less expensive and complex than their hodoscope counterparts. The approximate X and Y positions of scintillation events are interpolated from relative activations of the surrounding grid of PMTs. Thus, detected strike positions in scintillation camera architectures more closely resemble probability clouds than precise pixel resolutions. With multiple stacked panels, a best fit straight-line technique like that of the drift chamber method can be employed to improve tracking resolution. The system described in this report is based on these concepts but expanded into practice to serve as a practical development platform to investigate possibilities and limitations of the technique for muon instrumentation. Additionally, the system described here uses PMTs mounted to the corners of the panel, rather than as a grid on the surface, to simplify panel stacking and employ minimal PMTs.

Compared to layer addition in scintillator strip methods, panels can be added to the system to increase resolution with low additional cost. Furthermore, precise calibration of previously installed panels may make self-calibration of new panels possible by installing them between calibrated panels such that muon trajectories through the new panel are known. As another possible feature, a configurable panel system can be arranged into separate stacks, enabling non-invasive imaging of items by angle deflection from one stack to the next like what has been accomplished with drift chambers [8,23].

The system is implemented with three panels stacked on a modular frame to facilitate simple adjustment of acceptance angle and orientation. For position localization, it uses a regressive machine learning algorithm. Techniques and derivations are discussed for event localization and empirical data and practical observations are presented.

Materials and Methods

Mechanical Description and Frame Assembly

The core features of this muon detector are its three parallel, solid panels of plastic scintillator. Each scintillator panel is square (76x76x2 centimeter, LxWxH) - (30” x 30” x 0.8”) with chamfered corners for PMT mounting. ET Enterprises series 9266B 2-inch photomultiplier tubes are used to convert light signals from scintillation events into electrical signals. It has been suggested that summation of signals from two smaller PMTs mounted at each corner may improve performance and/or lower cost but for simplicity in this early work, a single PMT at each corner was used. Figure 1 depicts the completed muon detector in operation.

Figure 1.Scintillation camera muon tracker system with configurable frame.

The plastic scintillator material is EJ-200 by Eljen Technology and consists of polyvinyl toluene with organic fluors and a scintillation efficiency of 10,000 photons/MeV. The optimum scintillator wavelength band is approximately 400 nm to 500 nm with maximum emission at 425 nm. The rise time is 0.9 ns and the decay time is 2.1 ns with a pulse width of 2.5 ns measured at full width at half maximum (FWHM). Other scintillator options may improve performance in future work but EJ-200 was selected in this prototype evaluation for the versatility afforded by its combination of long optical attenuation length and fast timing.

To minimize the loss of photons produced during a scintillation event, the PVT was covered with 0.0254 mm (0.001”) thick commercial grade aluminum foil on all sides of the panel except for the PMT-accepting chamfers. The aluminum foil was then overlaid with a protective sheet of 0.762mm (0.03 inch) black vinyl. The spacing between PMT centers constitutes a square grid with side length of approximately 71.1 cm (28 in). In other words, adjacent PMTs have center-to-center spacing of L while cross-corner PMTs have spacing of L. Originally, three PMTs per panel were attached to the chamfered section of the PVT using a light guide bonded with EJ-500 clear, colorless optical cement, with the fourth corner left blind. A fourth PMT was later added to the blind corner for further investigation. EJ-500 cement was chosen because of its similar refractive index to the PVT scintillator panel (1.57 for the cement vs. 1.58 for the panel).

Each PVT/PMT subassembly is secured within a square frame that attaches to the larger aggregating frame (upper frame). The scintillator panel assemblies can slide lengthwise along the upper frame and be locked into position, creating a range of inter-panel spacing from approximately 7.62 cm (3 in) to 187.2 cm (73.7 in). Calculating from the inscribed circle of the square panels, this panel spacing range allows acceptance angle settings between approximately 44° and 169°. The overall detector frame can rotate relative to its base frame via a pair of bearing assemblies. The panel normal can be rotated in 5-degree increments through a full 360°.

The detector frame is constructed primarily from 6105-T5 aluminum extrusion. The complete assembly reaches a maximum height of 282 cm (111 in) and weighs approximately 340 kg (750 lb). In a vertical orientation, the assembly footprint is 152.4 cm (60 in) by 165.1 cm (65 in). The structure is designed to be modular, with the ability to install additional panels and to reposition them easily. The framing follows a design intended to facilitate rapid modification for deployment in unsheltered environments (i.e. field deployments). The consequence of providing an option for exterior deployment is a highly rigid frame with load ratings suitable to withstand high force winds, provided that suitable ground anchoring has been achieved. For indoor settings, the base frame is assembled with casters to enable easy movement. Because the lower frame is mounted on wheels and the upper frame is supported at its mass center, the detector can be easily repositioned (translation and rotation) by a single person.

The PMTs require high voltage (HV), low current power; depending on the desired sensitivity, the required voltage can range from 500 V to 1500 V direct current (DC). Small footprint, high-voltage HVM2000/12P power supplies are mounted on the frame near each PMT. The voltage can be adjusted from 20 V up to 2000 V to vary PMT sensitivities. To set PMT gains for these experiments, a mock muon source (discussed below) was placed on the panel near each PMT and the gain was adjusted to just below the signal clipping threshold. After each PMT was set near clipping, the source was placed on the center of the panel and the gains of each PMT were lowered to a point where their responses were nearly identical. Exact gain voltage values were not collected. The power supplies are energized by a 12 V line running along the detector frame and are connected to the PMTs using safe high voltage connectors (SHV). The sensitivity of each PMT is controlled by the excitation voltage provided by the HVM2000 power supply. Because each PMT has an independent power supply, each tube can be set individually.

The signal output of each PMT is connected to the data acquisition system (DAQ) through Bayonet Neill-Concelman (BNC) connectorized coaxial cables that are nearly identical in length and construction. The PMTs output electrical current pulses, which are then converted to voltage for recording. The current to voltage conversion is accomplished by termination with a 50Ω BNC feed-through adapter at the digitizer input.

Figure 2.Block schematic of instrumentation system and model of detector frame.

The instrumentation and controls of the detector are located on a cart holding a direct current (DC) power supply, computer, monitor, and a National Instruments data acquisition system Figure 2 for operation, and an oscilloscope and function generator for analysis purposes; all systems are powered using an uninterruptable power supply. Data is recorded using three National Instruments PXIe-5160 digitizer cards mounted in a PXIe-1073 chassis. Three separate cards were used to allow expansion to a configuration with 12 PMTs instead of 9. Configurable collection parameters including trigger threshold, sample rate, record length, muon discrimination, etc. are controlled with LabVIEW. Unless otherwise specified, data presented here were collected at 1.25 GS/s sample rate per channel. With digital upsampling, the time resolution of an individual pulse can be increased in post-processing as a convenient processing feature but is not relied upon for arrival time determination. The channels are synchronized to within 5 picoseconds, meaning that channel synchronization in the digitizer does not limit timing resolution. An event is triggered when the signal from one PMT on the center panel crosses a user-defined threshold. Each triggered event causes all channels in the system to synchronously record for 400 ns, 200 ns prior to trigger and 200 ns after trigger. Subsequent logic is applied in the collection software to only save event waveforms that meet a set of user-defined conditions (i.e. if greater than 6 PMTs responded in the record). A response on 6 PMTs guarantees the activation of at least 2 panels in the system. System configuration parameters including trigger thresholds and save logic are recorded in a configuration file stored with each dataset.

Panel Strike Order

Tracking the strike order of the panels enables rejection of noise from muons traveling in reverse through the tracker system, simultaneous imaging on both sides of the detector, and assessment of muon counts in either direction through the system.

A common method of PMT pulse timing on scintillators is to use a constant fraction discriminator (CFD). There are several variants of CFD systems, generally implemented using a dedicated separate piece of equipment. One example of an analog CFD operates by splitting each PMT pulse signal, attenuating one half and delaying the second half, and then using an analog comparator on the two signal halves to trigger when the amplitudes cross one another. Effectively, CFD systems trigger when the PMT signal reaches a defined fraction of the initial pulse height. The CFD technique reduces timing errors (“time walk”) which arise from different gain settings or pulse shapes. For example, if two PMT pulses are timed when the signals crossed a set voltage threshold, the PMT signal with higher amplitude may trigger earlier than the lower amplitude PMT signal, even if the signals were simultaneous. However, if the same signals are timed when they cross 20% of their respective pulse heights, time walk errors can be reduced or eliminated.

Because full PMT signal waveforms were digitized with a high sample rate data acquisition system, a CFD technique was applied digitally for each record in post-processing using MATLAB, rather than with a separate piece of equipment.

Figure 3.(A) Example data of one muon traveling through the detector with a recorded event on all nine PMTs. PMTs on Panel 1 are red, Panel 2 blue, and Panel 3 green. (B) The same data normalized by the maximum of each record. (C) Same data passed through bidirectional digital FIR filter. (D) Same data filtered and normalized by pulse peak. Sample spacing is 0.2 ns in these records due to upsampling

Overall, pulse shape was consistent in muon data collected with this system Figure 3, which further reduces time walk error. Normalization by the pulse peak clarifies that all nine PMT pulses are of similar shape and rise time.

With a maximum scintillator panel spacing of 187.2 cm (73.7 in) in this frame and a rated time spread of 2 ns in these PMTs, forward/reverse direction assessment through the detector is still not completely reliable. The calculated minimum time difference between the arrival of signals at PMTs on the two most distant panels is close to the combined time spread of two PMTs. If the earliest response on each panel is used to clock arrival at that panel, then the time difference between panels at the farthest ends of the frame can be as low as approximately 2.8 ns (muon strikes blind corner of the first panel, travels between panels, then strikes near a PMT corner on the third panel). In the described scenario, the muon travels 202.4 cm from the blind corner of the first panel to a PMT corner of the third panel, taking approximately 6.8 ns at 0.994c. However, the time for the scintillation light signal to travel through the panel to the PMT must be considered. The speed of light in the EJ-200 PVT material, based on the published refractive index, is approximately 0.633c. As the published refractive index is likely based on the standard 589 nm light measurement, the peak 425 nm scintillation light in the panel is expected to be even slower. Furthermore, as discussed in the timing section, only a small portion of scintillation light is expected to take a direct path from event to PMT, meaning that much of the signal will experience reflections and travel farther within the panel, further slowing signal propagation. Using liberal assumptions of only direct path photons and the 0.633c speed value, the time for light from a scintillation event to travel one edge length of 76 cm through the panel is approximately 4.0 ns. muon path described above, the difference in signal arrival time at the PMT on panel 1 and that on panel 2 is only 6.8–4.0 = 2.8 ns. However, if the median PMT arrival time on each panel is used for timing, the lowest time difference is approximately 4.4 ns. Considering the same muon path as above, signal arrival at the median PMT on the first panel would be the same value of ~4.0 ns after the scintillation. The muon would travel for the same ~6.8 ns to a PMT corner on panel 3, but now the scintillation light on panel 3 must travel a panel side length to reach the second PMT, and thus the arrival occurs approximately 4.0 ns after the scintillation. Effectively, the signal propagation times in panel 1 and panel 3 cancel and the difference in arrival from panel 1 to panel 3 is measured as 6.8 ns. In the worst-case scenario for the median approach, a muon strikes the blind corner of panel 1 and strikes the center point between two PMTs on an edge of panel 3. In this case, the median arrival on panel 1 is 4.0 ns after scintillation, the muon travels ~191 cm or ~6.4 ns, and the median arrival on panel 3 is ~2.0 ns after panel 3 scintillation. In this case, the system would register arrival time difference of 6.4-4.0+2.0 = 4.4 ns. Compared to 2.8 ns, the 4.4 ns time difference decreases errors expected from PMT jitter. For this reason, and for reduction of effects from noise or statistical variations of input signals, the median arrival on each panel was used for panel timing.

Even with the median arrival, the PMT time spread can affect the determination of forward versus backward traveling muons. For all data analysis in this report regarding forward and reverse muon counts, direction was only recorded if the signal arrivals on distant panels were calculated to be greater than 4.8 ns (6 digitizer samples) apart. This value was chosen because it is one sample longer than the 4 ns window afforded by the ~88% within ±2 ns time spread jitter value. Knowledge of the strike location on each panel would enable minimization of these timing effects and will be included in future work.

Single-Panel Interaction Position Determination

Triangulation by Timing

After a scintillation event, many photons take a long path with reflections within the panel before reaching each PMT. Only a small portion of the photons produced will take a direct path to each photo-detector, making direct triangulation by timing dependent on either the detection of these photons or on uniform pulse shapes from the scintillator panel. In addition to the low signal, the time spread jitter of the PMTs can introduce large errors in photons time-of-flight (TOF) triangulation. The PMTs used in this experiment are not optimized to minimize timing jitter. However, with the assumptions that the direct path photons are detectable or that pulse shapes are consistent throughout the panel, and that the error from time spread jitter is low, the position of the muon interaction point on a panel can be triangulated by direct derivation using the response timing of a minimum of three PMTs. On a single panel, fewer than three PMTs will create non-unique solution possibilities.

Because the zero time of the muon strike is not known, the position calculation depends on the differences in arrival time at each PMT. Knowing the refractive index of the scintillator plastic, differences in signal arrival time at each PMT can be converted to differences in distance traveled from the scintillation location. A system of equations can be solved for the X and Y positions of the scintillation event as shown in in Fig. 4 with accompanying equations.

Figure 4.Derivation diagram for triangulation by difference in distance to each PMT. In these equations, d1, d2, and d3 are not variables. Rather, differences in length (d1 – d2) and (d3 – d2) are measured quantities produced by the instrumentation from the difference in arrival time at their respective photo-detectors and the known speed of light in the scintillator. L is the side length of the panel. The closed-form solution to these simultaneous equations for X and Y is too long to reproduce here.

To experimentally test the derivation in Figure 4, the equations were used to triangulate a simulated muon strike from photons time-of-flight (TOF) information in the following way. Photon TOF was simulated by the propagation delay though varying lengths of coaxial cable. Three different lengths of coaxial cable were connected in parallel to the output of a function generator with BNC connectors. The lengths were selected to coordinate with a given position on the panel (Fig. 5). Short electronic pulses (~30 ns) created on a function generator were sent through all three BNC cables simultaneously using a signal splitter. The other end of each cable was connected to separate inputs on an NI PXie-5160 card of the DAQ. The muon data collection equipment was set with similar parameters to those for recording muon scintillations. Pulses from the BNC cables were digitized. The DAQ sample rate is high enough to detect the difference in signal delay through each of the cables. Because the pulses used were repetitive and identical, random interleaved sampling could be employed on the DAQ to increase the sample rate up to 5 GS/s on each channel.

To convert differences in cable propagation delays into differences in cable lengths, the differences in arrival time were multiplied by the published propagation speed of electric fields in RG-58 coaxial cable; (0.659c) where c is the speed of light in vacuum. The obtained propagation length differences (d1-d2 and d3-d2) were entered into the solution equations to the simultaneous set shown in Figure 4 and the “strike” position was computed. An error of less than 1 cm (0.39 in) was calculated from the computed position to the known correct position. A full statistical analysis was not performed in this proof of concept test. This small error is most likely due to error in the physical measurement of the cable lengths but may also be attributed to variance in construction of the cables, the signal splitter, or error in arrival time computation.

Figure 5.Fig. 5. (A) PVT scintillator panel with 3 PMTs attached at the corners. BNC cable intersection point is shown. (B) Superimposed actual vs computed plot over scintillator panel. (C) Pulses received on the 3 channels recorded. Small differences in arrival time can be seen. (D) Actual BNC cable intersection point compared to point computed from time differences.

The speed of electric field propagation in the coaxial cable used (0.659c) is close to the speed of light through the PVT scintillator plastic (0.633c). This experiment demonstrates the efficacy of the TOF position triangulation derivation and of the muon recording instrumentation to resolve timing on the scale necessary. The 9266B PMTs in the present system have a time spread jitter of ~4 ns (given in discussion with the manufacturer), which is too high to resolve panel strike position with this method reliably;however, with the addition of PMTs with lower jitter, direct TOF triangulation may be possible. For clarification, the time spread jitter was given as the full width half maximum (FWHM) of the expected probability curve of arrival times from a repeated input. The FWHM equates to 2.35 standard deviations which means that ~88% of transit times are expected to fall within ±2 ns. Due to this equipment limitation, the remainder of this paper focuses on triangulation by signal magnitude.

Position Calculation by Amplitude

The approximate position of the muon interaction point on the plastic scintillator can be determined using the magnitudes of the surrounding PMT responses. Initially, a similar derivation to that in Figure 3 was calculated based on the hypothesis that signal magnitude would scale in a simple predictable manner (i.e. inverse square) with distance from the PMT. However, empirical data showed that the signal magnitude did not scale in a simple manner with distance in these panels. With a light pulse source (discussed below) moved outward from a PMT, measured magnitudes vary nonuniformly with source position on the panel. This nonuniform relationship is likely attributed to nonuniform reflective panel covering, small inconsistencies in the light introduction through holes in the covering, PMT mounting angle variance, or other assembly effects. In a subsequent test, coverings were removed entirely from one panel and replaced with opaque black mat roofing paper with no reflective layer. While the signal was more predictable in the black mat panel, local maxima and minima in the data were still observable and repeatable (Fig. 9). With these results, direct mathematical derivation of strike position was abandoned in favor of a machine learning approach for evaluation of the present system.

Using a machine learning approach, the ideal way to both calibrate and determine the precision of the detector is to use a controlled muon beam. Not only would the exact strike position be known, but also the exact path the muons take through the detector. A secondary option would be to use a proven muon tracker placed on top of the new tracker with simultaneously triggered instrumentation to enable knowledge of the direction and location of naturally occurring muons passing through both systems. Without a muon beam or other proven tracker, calibration is limited to a panel-by-panel technique. One such calibration technique is to use a small piece of scintillator material with its own PMT connected to the same instrumentation system as the tracker and to only collect the muons that pass through both the small scintillator piece and the muon detector. This technique would use real muons to calibrate the system but was not used in this report due to lack of available equipment and the long collection times required to generate a sufficiently large dataset across the panel for training a machine learning algorithm. This coincident scintillation calibration technique will be revisited in future work. One discrete panel calibration/precision estimation technique which is used in gamma cameras is to use a source of ionizing radiation that will produce scintillation events in known locations. The PMT responses can be recorded for known source positions on the panel and the weighted sums can be adjusted so the computed position matches the known position of the source. Practical considerations for this work did not allow for use of an ionizing radiation source.

In lieu of an ionizing radiation source or other calibration equipment, a programmable flashing light emitting diode (LED) (470nm) was used to emulate scintillation events at known locations across the panel. The LED pulse parameters were set to closely mimic the PMT response of muons previously recorded with the system. The pulse parameters on the function generator powering the LED were set to a pulse width of 50 ns and voltage of 2.2 V. The LED can be pulsed much faster or slower and can be made brighter or dimmer if desired for a specific calibration dataset.

Figure 6.Hole-grid calibration sheet installed on the PVT panel. The holes are covered with individual tabs of black electrical tape. The flashing LED is mounted inside the black rubber puck shown.

To enable collection of LED calibration data, the original covering of the panel was replaced with a covering constructed of approximately 1 mm (0.04 in)-thick black acrylonitrile butadiene styrene (ABS) with a 9x9-point grid of 1.5 mm (1/16 in) diameter thru-holes (Fig. 6). As can be seen in Fig. 6, all grid points in the calibration panel except for the point being interrogated are covered with black electrical tape.

Figure 7.Original industrial aluminum foil wrapping beneath the outer black plastic layer (Left) and Mylar wrapping installed with calibration hole grid sheet (Right).

Installation of the calibration grid required that the previous covering of black vinyl and aluminum foil was removed, exposing the bare PVT and displaying the construction of the opposite aluminum foil covering Figure 7. The non-uniform surface of the aluminum foil prompted its replacement with Polyethylene Terephthalate, commonly known as Mylar®, which also proved difficult to keep flat without the use of adhesives or rigid cladding. In order to keep the reflective properties of the material the same on both sides of the PVT, Mylar sheeting was installed on the reverse side of the panel as well. Mylar has a reflectivity of approximately 0.98 versus 0.88 for the aluminum foil. Its introduction into the system could potentially be used to improve photon collection efficiency.

The collected calibration dataset from LED pulses was used to train a machine learning algorithm to learn the relationship between PMT data features and the known X and Y positions of the source. This learned relationship of inputs to outputs was then applied to another LED dataset to test the trained system. Data features extracted from the PMT data included pulse height, normalized pulse heights, pulse height differences, FWHM of the pulses, and arrival times. Each of these features was assessed for performance improvement of the machine learning tools but some features (i.e., arrival time and full width half max) showed no improvement, as expected from the consistency of pulse width and the large time spread jitter. Machine learning tools explored included support vector machines (SVM), artificial neural networks (ANN), and gaussian process regression (GPR). One SVM tool tested was the Joachims svm-light tool [39]. All other machine learning tools, including another SVM, were implemented from the machine learning toolbox in MATLAB. Although all these techniques were successful, the GPR method proved most accurate with these data.

While the Mylar is flatter and more reflective than the original aluminum foil, it still proved difficult to eliminate the wrinkles created in the Mylar surface in the compression mounting of the frame. The primary obstacle that arose was a bowing outward of the ABS sheet to which the Mylar was partially glued, which acted to create a series of large wrinkles on the Mylar surface. A better method would be to use a reflective coating on the PVT rather than a mechanically wrapped sheet or to use one of the following techniques: 1) use a more rigid sheet of ABS or backer material; 2) use a framing system that does not use compression; or 3) abandon reflective coating altogether and use only an opaque covering to avoid problems arising from inconsistent reflection.

Even with the wrinkled reflective material and local maxima and minima across the panel, the machine learning algorithm proved capable of accurately localizing LED flashes throughout the surface (shown in Results). This LED calibration data demonstrates that scintillation camera methods can be implemented successfully using few PMTs mounted on the corners of intact panels of scintillator material.

Three-Dimensional Muon Path Computation

As the machine-learned panel models have so far only been implemented using brief LED flashes rather than real scintillation events, three-dimensional muon paths in this report were computed for evaluation and conceptual demonstration only. The GPR model from the LED calibration was applied to each panel and all three X-Y strike positions were computed from real muon data. Panel spacing was entered into a MATLAB program along with the computed X-Y positions of muon strike locations on each panel. A three-dimensional line was fit to the data from each muon and plotted. These data are compared to simulated random muon tracks through the same detector configuration. Extensive improvements will be made by calibration with muon scintillations and confirmation of muon paths through the detector for use in imaging applications; this is the subject of future work.

Underground Tunnel Deployment

In a test field deployment, the muon detector system was deployed for 28 days approximately 59 m underground in the TA-41 tunnel at Los Alamos National Laboratory. Muon counts were recorded in the upward and downward direction throughout the deployment.

The system was driven from Albuquerque, NM to Los Alamos, NM and unloaded from the truck with a fork-lift. The upper frame of the detector was detached from its base frame in order to fit in the tunnel. Casters on the upper frame allowed the system to be wheeled down the flat concrete floor of the tunnel. Without the lower frame, data could only be collected in the vertical orientation. All data was collected in a configuration with the top and bottom panels as far apart as possible (~187 cm) and thus minimum acceptance angle (~44°) of the system, as this allows for the most reliable up versus down assessment.

The position of the detector in the tunnel corresponds to a position of previously published muon work from the same tunnel [7,16]. The position was 62 m into the tunnel with an approximate overburden of 59 m to the surface in volcanic tuff geology.


General System Performance

This prototype muon detector system works well as an early prototype and serves as a valuable platform for concept evaluation and informing future development. The system may be notably improved in future works based on what was learned here.

The mechanical frame assembly made for easy positioning of the detector by rotational orientation, tilt angle, and acceptance angle by a single person if necessary. Simple modifications could be made to the system to allow for all orientations to be mechanized and controlled digitally.

The electronics and instrumentation system performed as expected. Individual gain control on all channels was possible using the small high-voltage power supplies. All PMT channels could be digitized simultaneously with full waveform digitization and sample spacing of 0.8 ns after being triggered by a single PMT in the system. User defined settings in the LabVIEW software defined the signal trigger threshold and discrimination of whether the event was stored by the number of PMTs which responded within the length of the digital record.

Forward-Reverse directions through the detector system were successfully tracked for approximately 78% of muons in the laboratory and tunnel. The time spread of the PMT responses limits muon timing in the detector. Many muons have arrival differences below the 4.8 ns discriminated signal arrival difference and were counted as unresolved direction.

The system was successfully deployed in the TA-41 tunnel at Los Alamos National Lab and remained stable collecting data without interference for a period of 28 days when it was turned off. Upward and downward muon tracks were recorded and counted.

Triangulation of an electromagnetic signal source in the panel by direct mathematical derivation using TOF was successfully demonstrated. The instrumentation was set to 0.2 ns sample spacing in this experiment and the source was a short electrically simulated muon pulse sent through slightly different lengths of coaxial cable. The computed triangulation point was within 1 cm of the location measured with a measuring tape. While the method clearly works using cables, this time triangulation technique was not possible using these PMTs owing to the time spread jitter of the sensors. Low jitter PMTs (<0.1 ns) may enable this technique for muon localization on bulk scintillators in future work.

Muon interaction point with a panel was determined by the relative magnitudes of the responses on each PMT on the panel through machine learning techniques. Machine learning algorithms were able to learn the relationships between inputs and outputs of a training dataset and then apply those learned patterns to a test dataset.

Three-dimensional muon path tracking was implemented as a conceptual demonstration by fitting best fit lines to the X-Y position computed on each individual panel. While tracks were computed successfully, their accuracy was not tested, and we think can be greatly improved with very simple modifications such as adding a fourth PMT to the blind corner of each panel, creating training datasets from real muons using small coincident scintillator blocks instead of LED pulses, replacing wrinkled foil with rigid or painted or nonreflective coverings, using smaller panels, and adding more panels for the line fit.

Directional Determination and Underground Test

With the addition of multiple panels, muon paths can be computed, acceptance angle can be adjusted, direction through the detector can be determined, and scintillation events that do not pass through all (or a given subset) of the panels can be rejected. All such features were tested successfully with the present system.

The direction of a muon through the system can be determined from the synchronized dataset of PMT responses. A constant fraction discriminator technique is used to determine the first arrival of the pulse from each PMT. The direction is determined by the strike order of the panels. In the vertical position, most muons will be travelling downward. In the horizontal position, the nominal ratio of muons traveling forward versus backward is 50/50. The muon detector was configured with three PMTs per panel. For directionality experimentation, the data collection and storage settings were such that only events that struck all 3 panels were recorded. Muon datasets were collected at Sandia National Laboratories in Albuquerque, NM in both the upright and the inverted vertical configurations, and forward and reverse orientations with N-S and E-W alignment. The direction of the muons for the datasets is tabulated in Table 1.

Data Set Up/North/East Traveling Muons Down/South/WestTraveling Muons Unresolved Direction Muons
Vertical Up-Down 93 (<1%) 9,890 (>99%) 9,549
Vertical Down-Up 110 (<1%) 17,131 (>99%) 1,824
Horizontal N-S 761 (~54%) 664 (~46 %) 1,136
Horizontal S-N 685 (~53%) 621 (~47%) 1,358
Horizontal E-W 715 (~49%) 742 (~51%) 1,260
Horizontal W-E 604 (~51%) 591 (~49%) 1,325
Table 1.Tabulated data of direction of muons through the detector in both orientations for vertical and horizontal positions. Each dataset is 1 hour of data collection. The unresolved direction muons are those where difference in arrival times at the outermost panels is less than 4.8 ns.

Similar ratios of muons travelling forward versus backward through the detector are shown when the detector is rotated 180 degrees, for both the vertical and the horizontal muon datasets. As expected, the vertical datasets show over 99% of muons traveling downward. Also as expected, the North-South horizontal datasets show approximately 50%-50% North-traveling versus South-traveling and East-traveling versus West-traveling muons.\

Figure 8.Mean upward traveling muon counts (Left) and mean downward traveling muons (Right) over a 28-day data collection in the TA41 tunnel in Los Alamos, NM. Each bar represents average muon counts over a two-hour period at each time of day. The error bars show standard error about the mean.

However, the muon directions for the horizontal data are not exactly 50%-50%. Comparison of the N-S configuration dataset to the S-N configuration dataset suggests that there may have been slightly more North-traveling muons during the data collection than South. Additionally, the data showing more North-traveling is less exaggerated in the S-N dataset than the N-S dataset while comparison of the E-W dataset to the W-E dataset shows a slight reversal. One interpretation of these data is that the detector system has an inherent bias toward collecting muons traveling in one direction through the system over the other. This may be due to the angle, mount efficiency, or construction of the PMT being used to trigger the acquisition. A longer data collection may be required to determine if this suggested bias is real. A more appropriate trigger method may be to trigger an event on the activation of any PMT in the system rather than on only one. The event discrimination by number of PMTs activated could still be applied and likely fewer muons would be missed. The much higher unresolved direction muons in the Up-Down configuration compared to the Down-Up configuration may also be attributed to a difference in detection efficiency in one direction versus the other on the trigger PMT. Other possible causes for directional bias may be slight variance in cable construction, PMT assembly, gains, or time spread in the detectors. Finally, the number of muons that could not be resolved for direction is higher than expected, meaning the directional percentages have a high amount of error. Lower jitter PMTs or farther panel spacing could be used correct this problem.

Muon count data in the upward and downward directions over the 28-day data collection in the TA-41 tunnel are shown in Fig. 8. Each bin represents a two-hour window of data collection centered at the labeled time of day. The means were computed for each daily time window over the 28-day period with error bars showing standard error. These data indicate approximate counts of 1-2 upward moving muons and 550 downward muons per hour. Muons of unresolved direction are not included.

The muon counts per unit area are lower than those previously reported [16] from the same location. This is likely due to the much lower acceptance angle of the present system. Other possible contributors to fewer muon counts may be missed muons from suboptimal trigger threshold/trigger method and exclusion of the muons with unresolved direction.

Single Panel Machine Learning Calibration

A training dataset was created using pulses of a 470 nm wavelength LED at various positions across the panel. The relationship of signal magnitude and position on the panel was determined to be non-uniform (Fig. 9) and thus the initial mathematical derivations of position were inappropriate for this system. A machine learning technique was adopted.

Figure 9.LED Calibration dataset from each of the four corner-mounted PMTs. The X and Y axes show the calibration grid. The Z-axis is the maximum voltage magnitude averaged over 200 LED flashes at each point.

Features of the PMT data in the training set were extracted including magnitude, normalized magnitude, full-width half-maximum, and time of arrival. Machine learning systems were trained to learn the relationship between the extracted signal features and the known X and Y positions on the panel. The learned input/output relationship was used to analyze a separate XY dataset to determine the efficacy of the algorithm. Support Vector Machines (SVMs), artificial neural networks (ANNs), and Gaussian Process Regression (GPR) machine learning techniques were all attempted with success but GPR provided the most consistent results in these experiments. All relevant figures show GPR results. It was determined that full-width half-maximum and time of arrival features did not improve performance, which was expected due the time spread of the PMTs and the similarity of pulse shapes and rise times.

Actual Position Computed Position Standard Deviation Computed Position Standard Deviation
X Y X Y σ1 σ2 X Y σ1 σ2
15.24 15.24 15.247 15.242 0.187 0.386 15.196 15.24 0.142 0.231
22.86 22.86 20.571 22.809 0.523 1.887 22.479 22.763 0.411 1.297
30.48 30.48 30.698 30.640 0.668 1.577 30.490 30.553 0.297 1.010
38.1 38.1 38.115 37.861 0.347 1.094 38.084 37.777 0.307 1.031
45.72 45.72 45.844 45.798 0.452 1.074 45.758 45.811 0.053 0.939
53.34 53.34 53.23 53.261 0.424 1.978 53.332 53.332 0.058 0.160
60.96 60.96 60.751 60.853 0.238 1.173 60.932 60.921 0.010 0.449
60.96 15.24 60.853 15.344 0.043 0.680 60.934 15.290 0.060 0.307
53.34 22.86 54.604 22.852 0.381 1.811 54.058 22.857 0.469 1.529
45.72 30.48 46.769 32.473 0.538 3.995 46.261 31.280 0.441 2.811
30.48 25.72 30.444 45.961 0.238 0.942 30.457 45.857 0.335 0.670
22.86 53.34 22.867 53.309 0.093 0.243 22.890 53.317 0.129 0.294
15.24 60.96 15.354 60.825 0.025 0.629 15.283 60.909 0.005 0.220
Table 2.Tabulated machine learning performance data comparing results from using 3 PMTs against results from using 4 PMTs at the same gain setting. All values are in cm. With the addition of the 4th PMT, the mean area of the ellipses of these select locations was approximately 42% smaller than in the 3 PMT configuration with the center point ellipse being ~17% smaller, the three initially instrumented corners at ~29% smaller, and largest improvements showing in the previously blind corner at ~96% smaller.

The PMT responses to the calibration dataset from the four corners of the panel were recorded with the muon data recorder hardware and software. Normalized magnitudes of the PMT responses and the known X and Y positions were recorded for each flash of the LED. Data from two-hundred LED flashes from each of the 81 grid positions were used to train a GPR machine learning tool in MATLAB. A separate verification dataset of one-hundred more LED flashes was used to test the trained algorithm. The resultant calculated strike positions of the verification dataset were used to generate statistical Gaussian ellipses to estimate precision in various locations on the panel. Figure 10 depicts the computed positions from the verification dataset using a GPR trained on the training dataset. Configurations using different numbers of PMTs and low versus high gain were tested. The highest precision implementation uses all 4 PMTs at high gain, although 3 PMTs at high gain also performed well.

Figure 10.Results of machine learning algorithm trained on calibration data, tested on separate data. Horizontal and vertical axes are in inches. Four configurations of PMTs are shown: With only 2 PMTs and high gain on the left corners (Top Left), with 4 PMTs and low gain (Top Right), with 3 PMTs and high gain on the top left, bottom left, and bottom right corners (Bottom Left), and with 4 PMTs and high gain (Bottom Right). Red ellipses denote the area where 90% of the test data is computed.

Each red ellipse in Figure 10 shows the elliptical radii where 90% of the predicted data falls for that position. The intervals were calculated based on the verification data set which was not used in training the GPR. The probability cloud for each point is shown in white and is normalized by its height for the display purposes of this plot. Table 2 shows the actual known position, predicted center position, and long and short-axis ellipse standard deviations of the computed verification data interaction point positions plotted in Fig. 10. Spatial resolution of the panel is estimated from these data. In the best configuration, the largest ellipse was approximately 3.9 cm2 with radii of 0.44 cm and 2.8 cm. The mean ellipse area was 0.84 ± 0.33 cm2 and mean short and long radii were 0.22 ± 0.05 cm and 0.86 ± 0.21 cm respectively. The range denotes standard error about the mean. This is interpreted as a <2.8 cm spatial resolution with higher resolution in most places.

With PMTs on each of the four corners, probability clouds in symmetric positions about the panel center should be of similar size. However, the system is sensitive to nonuniformities in the reflective covering around the PVT panel. In addition to a tendency to wrinkle and an unavoidable wrapping seam, the reflective aluminum foil develops slack which creates non-uniform space between the foil and the panel. A large change in PMT response was observed during LED calibration when the foil is pressed flat compared to being left loose. Additionally, variances in coupling efficiency of the PMTs and uniformity of the panel may be sources of asymmetries in recorded data.

3-D Muon Path Tracking

As a conceptual demonstration, 3-dimensional lines are fit to the individual X and Y positions calculated from the GPR model on each panel. The GPR model used was from the LED calibration of one panel in the 3-PMT high-gain configuration. Figure 11 shows a set of 100 random simulated muon tracks and a set of 100 computed muon tracks from real data.

Figure 11.MATLAB simulated muon tracks generated by random XY positions on panels 1 and 3 with a line drawn between (Left) and muon tracks computed from real data as a conceptual demonstration (Right). The model generated from the calibration dataset from the center panel was applied to all three scintillator panels.

The computed muon tracks tended to avoid the edges of the panels rather than randomly distribute across the surface as expected. This grouping away from the edges indicates that the tracks are not always accurate with this initial test. Muons that struck panels near the edges were likely given solutions that fell inside the bounds of the LED training dataset. The hole grid for LED calibration does not extend to the edges of the panel which likely makes extrapolation beyond the grid unreliable. A calibration grid that extends to the outer edges of the panel or an acceptance technique that ignores edge effects is recommended.

There are many ways that muon track accuracy can be improved and confirmed over this conceptual demonstration with minimal or even zero changes to the system. First, while the recorded PMT signals from LED pulses are visually similar in pulse width, shape, and amplitude to those of muons, LED flashes are not a perfect representation of scintillation events; the panel can be calibrated using muons as described in Materials and Methods. Second, the same GPR model parameters were applied to all three panels rather than using individual models generated from the calibration data of each panel; individual calibration sets for each panel should be generated. Third, non-uniformities in the reflective coverings can be eliminated by constructing the panels with reflective adhesive coating or even a nonreflective covering. Finally, only 3 points are used to fit the linear muon track; additional panels could improve resolution. Each of these possibilities are subjects of future work for this detector system.


Stacked plastic scintillator panels with side or corner mounted PMTs may be an attractive technology to muon tracking experimentalists due to their modularity, ease of use, and stacking simplicity. The potential for a low number of PMTs and simple configurability could lead to low-cost systems with many features, which could make muon tracking accessible to a wider base of researchers and laboratories. The system described here has three panels but could easily be modified to have many more. Each additional panel improves the precision computed tracks by making a best straight-line fit to the data. After calibration of the initial panels, additional panels could be automatically calibrated. As an example feature of a many-panel system, an object could be placed between panel stacks and the angle of muon deflection from one stack to the next could be used to image the density structure of the object similar to previous work [8,23] performed using drift tube muon trackers. The frame of the present detector was designed to be robust, with outdoor applications and field deployments in mind but panels could be made smaller and less bulky frames can be constructed; trading some muon flux for potentially improved resolution and system maneuverability.

The LED calibration experiments provide insight into the implementation of scintillation camera or “Anger logic” on large sheets of plastic scintillator with widely spaced PMTs mounted at the corners. Where PMTs with low transit time jitter could theoretically be used to directly triangulate muon interaction points in the present system, direct computation by magnitude appears to be more complex. The traditional gamma camera approach computes a weighted sum (or difference) of response energies of a PMT grid attached to a scintillator crystal. The system is calibrated by placing a collimated gamma source in known locations and adjusting the weighting coefficients until the computed position matches the known position. This technique is effective with a closely spaced grid placed a distance from the edge of the scintillator to avoid edge effects. However, with the system described herein, having PMTs mounted on the corners of a large slab of plastic scintillator, some assumptions of the traditional approach are no longer valid.

The theoretical relationship of measured intensity as a function of distance from the PMT is exponential. While this relationship is supported by experimental data from LED flashes over most of the panel, the data shows deviation near the edges and corners as well as important local maxima and minima in many locations even beyond half-way across the panel. The LED evaluations are not a perfect representation of scintillation events in the panel and error may have been introduced by small discrepancies of introduction hole geometries, PMT mount angle, and/or other assembly imperfections. However, observations from this dataset suggest that calibration of the optical intensity response over the full panel is necessary in practice. The complex but repeatable intensity characterization in the LED test suggests that machine learning from a calibration dataset is an appropriate approach to implementing “Anger logic” on systems like the one described.

Additional complexities are introduced by the practical difficulty of applying a uniform reflective covering to the faces of the PVT panel to improve photon collection efficiency. The tendency of applied foil to crease, wrinkle and bow leads to important nonlinearities across the panel. It was observed in this work that PMT responses can be considerably altered by simply flattening a bowed section of the foil against the panel face. Non-uniform reflective coverings make the intensity versus distance relationship less predictable and could cause large errors in data interpolation. A more uniform construction could be accomplished using rigid reflective materials, or an optical coating deposited directly to the surface of the PVT.

Performance of the machine learning technique to accurately localize LED flashes even with wrinkled reflective panel covering demonstrates its robustness in this application. Observance of the larger localization clouds near the center of the panel suggest that higher resolution could be achieved with slightly smaller panels.

Data in this report suggests that a low-cost muon tracker can be developed with few PMTs and intact panels of plastic scintillator material. The data further suggest that the sophisticated DAQ system employed here may not be necessary for muon tracking in similar systems after initial developments are complete because full PMT waveform digitization on all channels will not be necessary for operation. The system developed here serves as a flexible, expedient platform for development of scintillation camera muon tracker techniques.

Conclusions and Future Work

A scintillation camera type implementation of a muon tracker was developed using solid square slabs of plastic scintillator material with photomultiplier tubes mounted at the corners. The relative response amplitudes of each PMT were used to compute the strike position on the panel by prediction from a machine learning algorithm trained from muon-mimicking LED flashes on a grid. The effects of PMT gain and number were assessed. Muon directionality measurement was tested in Up-Down, North-South, and East-West configurations with measured direction count ratios close to the expected values of >99/<1, 50/50, and 50/50 respectively.

A direct triangulation method based on TOF was derived and effectively demonstrated using cables of varying length. However, the TOF method was determined infeasible with the current system hardware due to timing jitter on the PMTs. The timing jitter also increases the necessary total system length in order to confidently determine the direction a muon goes through the system. A direct mathematical derivation like that of the TOF method was not applied to amplitude triangulation method due to the nonuniformities in the amplitude response of the panel.

It was concluded that scintillation camera techniques are viable for building a muon tracking system with minimal PMTs and other parts. A machine learning algorithm was implemented to learn the relationship of X and Y panel positions to the PMT response magnitudes. While 3 PMTs are all that is required, addition of a fourth PMT on the blind corner improves resolution and reliability. Position probability clouds are larger toward the center of the panel, indicating that slightly smaller panels may improve resolution.

Future work will include PMTs optimized to minimize timing jitter and overall detector size, improved uniformity in the panel coatings, triggering data acquisition on any PMT rather than just one, radiation source and/or coincident scintillator calibration methods, calibration to the edges of each panel, computed muon path verification, and image generation applications. The combination of TOF and machine learned amplitude methods will be investigated to enhance the performance of the system.

In future work, the panel stack will also enable computation of spatial resolution of the integrated system using data from muon track linear fits. Over a long data collection, the distribution of fit residuals on each panel will produce a statistical representation of resolution. With precise calibration data for the machine-learned algorithm on all three panels, time resolution is also expected to improve. Knowledge of strike position enables knowledge of light signal propagation time in the panel which can be corrected for.

The stack of 3 panels and the muon rejection and saving logic enables an estimation of detection efficiency of the center panel by comparison with the top and bottom panels. The method works because any muons that pass through the top and bottom panels must also pass through the middle panel. This was a missed opportunity in the present data because the integrated system was always triggered from the center panel, meaning the data will not show missed muons on the center panel. By comparison, coincident strikes on the top and center panels does not enable estimation of detection efficiency of the bottom panel because muons can take paths that miss it. Future data acquisitions will trigger from a top or bottom panel or by OR logic.


The authors would like to thank the Sandia National Laboratories, Laboratory Directed Research and Development office for funding this work. The authors would also like to thank Dennis King and Elton Wright at Sandia National Laboratories for their indispensable assistance with frame assembly and mobilization. The authors would also like to thank David Bonal from National Instruments for his assistance with hardware specification and custom LabVIEW software development. Special thanks are also given to Charlotte Rowe and Elena Guardincerri at Los Alamos National Laboratory for their assistance and monitoring support in deploying the muon detector in the TA-41 tunnel. Finally, we would like to thank Mattie Hensley and Andrew Wright at Sandia National Laboratories, J. Andrew Green at Mission Support and Test Services and others for their valuable reviews and comments.

Disclosure Statement: No competing financial interests exist.


  1. Crouch M. F., Landecker P. B., Lathrop J. F., Reines F., Sandie W. G., Sobel H. W., Coxell H., Sellschop J. P. F.. Cosmic-ray muon fluxes deep underground: Intensity vs depth, and the neutrino-induced component. Physical Review D. 1978; 18(7)DOI
  2. Allkofer O.C., Carstensen K., Dau W.D.. The absolute cosmic ray muon spectrum at sea level. Physics Letters B. 1971; 36(4)DOI
  3. Lee T. D., Robinson H., Schwartz M., Cool R.. Intensity of Upward Muon Flux due to Cosmic-Ray Neutrinos Produced in the Atmosphere. Physical Review. 1963; 132(3)DOI
  4. Neddermeyer Seth H., Anderson Carl D.. Note on the Nature of Cosmic-Ray Particles. Physical Review. 1937; 51(10)DOI
  5. Bonechi Lorenzo, D’Alessandro Raffaello, Giammanco Andrea. Atmospheric muons as an imaging tool. Reviews in Physics. 2020; 5DOI
  7. Guardincerri Elena, Rowe Charlotte, Schultz-Fellenz Emily, Roy Mousumi, George Nicolas, Morris Christopher, Bacon Jeffrey, Durham Matthew, Morley Deborah, Plaud-Ramos Kenie, Poulson Daniel, Baker Diane, Bonneville Alain, Kouzes Richard. 3D Cosmic Ray Muon Tomography from an Underground Tunnel. Pure and Applied Geophysics. 2017; 174(5)DOI
  8. Schultz L.J., Borozdin K.N., Gomez J.J., Hogan G.E., McGill J.A., Morris C.L., Priedhorsky W.C., Saunders A., Teasdale M.E.. Image reconstruction and material Z discrimination via cosmic ray muon radiography. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment. 2004; 519(3)DOI
  9. Lesparre N., Gibert D., Marteau J., Déclais Y., Carbone D., Galichet E.. Geophysical muon imaging: feasibility and limits. Geophysical Journal International. 2010; 183(3)DOI
  10. Schultz Larry J., Blanpied Gary S., Borozdin Konstantin N., Fraser Andrew M., Hengartner Nicolas W., Klimenko Alexei V., Morris Christopher L., Orum Chris, Sossong Michael J.. Statistical Reconstruction for Cosmic Ray Muon Tomography. IEEE Transactions on Image Processing. 2007; 16(8)DOI
  11. Alvarez L. W., Anderson J. A., Bedwei F. E., Burkhard J., Fakhry A., Girgis A., Goneid A., Hassan F., Iverson D., Lynch G., Miligy Z., Moussa A. H., Sharkawi M., Yazolino L.. Search for Hidden Chambers in the Pyramids. Science. 1970; 167(3919)DOI
  12. Tanaka H, Nagamine K, Kawamura N. Development of the Cosmic-Ray Muon Detection System for Probing Internal-Structure of a Volcano. Hyperfine Interactions. 2001; 138:521-526. DOI
  13. Nagamine K, Iwasaki M, Shimomura K, Ishida K. Method of probing inner-structure of geophysical substance with the horizontal cosmic-ray muons and possible application to volcanic eruption prediction. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment. 1995; 356(2-3)DOI
  14. Bonneville Alain, Kouzes Richard, Yamaoka Jared, Lintereur Azaree, Flygare Joshua, Varner Gary S., Mostafanezhad Isar, Guardincerri Elena, Rowe Charlotte, Mellors Robert. Borehole muography of subsurface reservoirs. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences. 2018; 377(2137)DOI
  15. Lesparre N., Marteau J., Déclais Y., Gibert D., Carlus B., Nicollin F., Kergosien B.. Design and operation of a field telescope for cosmic ray geophysical tomography. Geoscientific Instrumentation, Methods and Data Systems. 2012; 1(1)DOI
  16. Saracino G., Amato L., Ambrosino F., Antonucci G., Bonechi L., Cimmino L., Consiglio L., Alessandro R. D.’, Luzio E. De, Minin G., Noli P., Scognamiglio L., Strolin P., Varriale A.. Imaging of underground cavities with cosmic-ray muons from observations at Mt. Echia (Naples). Scientific Reports. 2017; 7(1)DOI
  17. Schouten Doug. Muon geotomography: selected case studies. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences. 2018; 377(2137)DOI
  18. Gibert Dominique, Beauducel François, Déclais Yves, Lesparre Nolwenn, Marteau Jacques, Nicollin Florence, Tarantola Albert. Muon tomography: Plans for observations in the Lesser Antilles. Earth, Planets and Space. 2010; 62(2)DOI
  19. Marteau J., Gibert D., Lesparre N., Nicollin F., Noli P., Giacoppo F.. Muons tomography applied to geosciences and volcanology. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment. 2012; 695DOI
  20. Pesente S., Vanini S., Benettoni M., Bonomi G., Calvini P., Checchia P., Conti E., Gonella F., Nebbia G., Squarcia S., Viesti G., Zenoni A., Zumerle G.. First results on material identification and imaging with a large-volume muon tomography prototype. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment. 2009; 604(3)DOI
  21. Borozdin K.N., Gomez J.J., Hogan G.E., Morris C.L., Priedhorsky W.C., Saunders A., Schirato R.C., Schultz L.J.. Scattering muon radiography and its application to the detection of high-Z materials. 2003 IEEE Nuclear Science Symposium. Conference Record (IEEE Cat. No.03CH37515). 2003. DOI
  22. Jonkmans G., Anghel V.N.P., Jewett C., Thompson M.. Nuclear waste imaging and spent fuel verification by muon tomography. Annals of Nuclear Energy. 2013; 53DOI
  23. Riggi S., Antonuccio-Delogu V., Bandieramonte M., Becciani U., Costa A., La Rocca P., Massimino P., Petta C., Pistagna C., Riggi F., Sciacca E., Vitello F.. Muon tomography imaging algorithms for nuclear threat detection inside large volume containers with the Muon Portal detector. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment. 2013; 728DOI
  24. Bandieramonte Marilena, Delogu Vincenzo Antonuccio, Becciani Ugo, Costa Alessandro, Massimino Piero, Pistagna Costantino, Riggi Simone, Sciacca Eva, Vitello Fabio, Rocca Paola La, Petta Catia, Riggi Francesco. Automated object recognition and visualization techniques for muon tomography data analysis. 2013 IEEE International Conference on Technologies for Homeland Security (HST). 2013. DOI
  25. Rocca P La, Antonuccio V, Bandieramonte M, Becciani U, Belluomo F, Belluso M, Billotta S, Blancato A A, Bonanno D, Bonanno G, Costa A, Fallica G, Garozzo S, Indelicato V, Leonora E, Longhitano F, Longo S, Presti D Lo, Massimino P, Petta C, Pistagna C, Pugliatti C, Puglisi M, Randazzo N, Riggi F, Riggi S, Romeo G, Russo G V, Santagati G, Valvo G, Vitello F, Zaia A, Zappalà G. Search for hidden high-Z materials inside containers with the Muon Portal Project. Journal of Instrumentation. 2014; 9(01)DOI
  27. Clarkson A., Hamilton D.J., Hoek M., Ireland D.G., Johnstone J.R., Kaiser R., Keri T., Lumsden S., Mahon D.F., McKinnon B., Murray M., Nutbeam-Tuffs S., Shearer C., Staines C., Yang G., Zimmerman C.. The design and performance of a scintillating-fibre tracker for the cosmic-ray muon tomography of legacy nuclear waste containers. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment. 2014; 745DOI
  28. Marteau Jacques, d’Ars Jean de Bremond, Gibert Dominique, Jourde Kevin, Gardien Serge, Girerd Claude, Ianigro Jean-Christophe. Implementation of sub-nanosecond time-to-digital convertor in field-programmable gate array: applications to time-of-flight analysis in muon radiography. Measurement Science and Technology. 2014; 25(3)DOI
  29. Uchida Tomohisa, Tanaka Hiroyuki K. M., Tanaka Manobu. Development of a muon radiographic imaging electronic board system for a stable solar power operation. Earth, Planets and Space. 2010; 62(2)DOI
  30. Ampilogov N. V., Astapov I. I., Barbashina N. S., Borog V. V., Chernov D. V., Dmitrieva A. N., Kompaniets K. G., Petrukhin A. A., Shutenko V. V., Teregulov A. I., Yashin I. I.. Large area scintillation muon hodoscope for monitoring of atmospheric and heliospheric processes. Astrophysics and Space Sciences Transactions. 2011; 7(3)DOI
  31. Bugg W., Efremenko Yu., Vasilyev S.. Large plastic scintillator panels with WLS fiber readout: Optimization of components. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment. 2014; 758DOI
  32. Anghel V., Armitage J., Baig F., Boniface K., Boudjemline K., Bueno J., Charles E., Drouin P-L., Erlandson A., Gallant G., Gazit R., Godin D., Golovko V.V., Howard C., Hydomako R., Jewett C., Jonkmans G., Liu Z., Robichaud A., Stocki T.J., Thompson M., Waller D.. A plastic scintillator-based muon tomography system with an integrated muon spectrometer. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment. 2015; 798DOI
  33. Jo Woo Jin, Kim Hyun-Il, An Su Jung, Lee Chae Young, Baek Cheol-Ha, Chung Yong Hyun. Design of a muon tomography system with a plastic scintillator and wavelength-shifting fiber arrays. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment. 2013; 732DOI
  34. Aguiar P., Casarejos E., Silva-Rodriguez J., Vilan J. A., Iglesias A.. Geant4-GATE Simulation of a Large Plastic Scintillator for Muon Radiography. IEEE Transactions on Nuclear Science. 2015; 62(3)DOI
  35. Anger Hal O.. Scintillation Camera. Review of Scientific Instruments. 1958; 29(1)DOI
  36. Joachims T. Advances in Kernel Methods - Support Vector Learning. MIT Press; 1999.
  37. GEORGE, E. P.. Cosmic rays measure overburden of tunnel. Commonw. Eng. 1955; 455
  38. Oláh László, Barnaföldi Gergely Gábor, Hamar Gergő, Melegh Hunor Gergely, Surányi Gergely, Varga Dezső. Cosmic Muon Detection for Geophysical Applications. Advances in High Energy Physics. 2013; 2013DOI
  39. Borozdin Konstantin N., Hogan Gary E., Morris Christopher, Priedhorsky William C., Saunders Alexander, Schultz Larry J., Teasdale Margaret E.. Radiographic imaging with cosmic-ray muons. Nature. 2003; 422(6929)DOI