MUYSC: An end-to-end muography simulation toolbox

Muography is an imaging technique based on attenuation of the directional muon flux traversing geological or anthropic structures. Several simulation frameworks help to perform muography studies by combining specialised codes from the muon generation (CORSIKA and CRY) and the muon transport (GEANT4, PUMAS, and MUSIC) to the detector performance (GEANT4). This methodology is very precise but consumes significant computational resources and time. In this work, we present the end-to-end python-based MUographY Simulation Code. MUYSC implements a muography simulation framework capable of rapidly estimating rough muograms of any geological structure worldwide. MUYSC generates the muon flux at the observation place, transports the muons along the geological target, and determines the integrated muon flux detected by the telescope. Additionally, MUYSC computes the muon detector parameters (acceptance, solid angle, and angular resolution) and reconstructs the 3-dimensional density distribution of the target. We evaluated its performance by comparing it with previous results of several simulation frameworks.


INTRODUCTION
Muography is a non-invasive technique used to explore geological or anthropic structures acquiring two-dimensional density distributions by registering the integrated muon flux passing through the target.Muon tomography can reconstruct three-dimensional density distributions combining multi-projection muography data, considering the observation site, the screening of neighbourhood mountains and the muon detector parameterization.
Muon tomography reconstructs three-dimensional density distributions of the target from two-dimensional density projections.There exist two main techniques: Inversion methods and filteredback projections.Inversion methods combine muography with gravimetry data (Tanaka et al.(2010); Nishiyama et al.(2014); Guardincerri et al.(2017); Vesga-Ramírez et al.(2021)).Barnound and collaborators inverted gravimetric and muographic data of the Puy de Dôme volcano.They conclude that muographic data cannot retrieve a 3D distribution of the volcano from a single observation point -It's necessary to combine it with gravimetric data (Barnoud et al.(2021)).Marteau and collaborators combined muography information (3 and 5 telescopes) with gravity data of the La Soufriére de Guadeloupe lava dome.The results showed density anomalies inside the lava dome with good agreement with electrical tomography data (Marteau et al.(2017); Rosas-Carbajal et al.( 2017)).
Filtered Back Projection methods use several muon radiography recordings around the target.Nagahara and collaborators per-formed a muon tomography of the Omuroyama volcano with an internal density distribution similar to a checkerboard.The reconstruction errors depend on the number of observation points (they evaluated 4, 8, 16, 32, and 64 points) and the exposure time (Nagahara & Miyamoto(2018)).
Different software tools form composite simulation frameworks in muography: from cosmic ray showers to detector response.In the Monte Carlo approach, authors use specialised codes such as CORSIKA (COsmic Ray SImulations for KAscade), CRY, ARTI, and EcoMug (Moussawi et al.(2022); Sarmiento-Cano et al.(2022); Samalan et al.(2022); Taboada et al.(2022)).COR-SIKA simulates the evolution of extensive air showers in the atmosphere step-by-step, giving high accuracy in particle spectrum composition.Its primary disadvantage is the high simulation time (hours or even days to compute the differential muon flux at the observation site).CRY contains parametrizations of all particles, including background, with execution times of the order of minutes.EcoMug generates cosmic muons according to a differential flux characterisation provided by a user or with a default model based on data from the ADAMO experiment.Its execution times are in the order of seconds.Semi-empirical models offer a fast methodology to compute muon fluxes.They fit muon flux models considering the zenithal incoming angle, muon energy, and altitude, obtaining comparable results with Monte Carlo-based simulation frameworks (Lechmann et al.(2021); Lesparre et al.(2010)).
Other Monte Carlo methods (GEANT4, PUMAS, and MU-SIC) simulate the muon-target interaction: GEANT4 provides accurate muon tracking and geometry/material target structure, allowing the detailed study of the muon scattering effects.PUMAS library is used for transporting muon and tau leptons in matter.It operates in a fast-deterministic or detailed Monte Carlo mode (Niess(2022)).MUSIC integrates energy the losses of muons due to: ionization, bremsstrahlung, electron-positron pair production and muon-nucleus inelastic scattering (Kudryavtsev(2009)).
Semi-empirical methods perform muon transport employing models that parameterize the stopping power of muons in different materials based on the Groom dataset (Groom et al.(2001); Lesparre et al.(2010)).
The muon detector's performance depends on geometrical and functional features: acceptance, solid angle, and angular resolution, characterising the overall muon detector performance.These geometrical parameters depend on the pixel area, the number of pixels and the distance between the detection layers (Lesparre et al.(2010)).Generally, muography experiments address the detector response through GEANT4 simulations, providing information about the energy and time-of-flight of muons crossing the detector.GEANT4 allows us to create sensitive detector matrices with specific materials and photo-detector quantum efficiencies.We also can evaluate the installation of passive layers (lead or iron) for background filtering (Vásquez-Ramírez et al.(2020); Moussawi et al.(2022); Samalan et al.(2022)).Semi-empirical models deal with detector response simulations based on data collected from an operating detector.These models include sensor features (photoelectron spectrum, gain, temperature dependency) and signal features (attenuation, delay, and noise) (Peña Rodríguez(2021)).
This paper introduces MUYSC (MUographY Simulation Code), a python-based framework for simulating muon radiography (atmospheric muon generation, muon transport, and geometrical detector parameterization) and muon tomography (threedimensional reconstruction) of geological targets.In section 3, we describe the muon radiography method used by MUYSC, while section 5 displays the methodology employed by the MUYSC muon tomography module.Section 4 is about how MUYSC computes the muon detector parameters.Finally, in section 6, we wrap up some concluding remarks.

MUYSC FRAMEWORK
MUYSC is an open-source, python-based code with three main functions as shown in Figure 1: The muon radiography module inputs are the geological coordinates of the target topography, the observation point coordinates, the projection point coordinates, the telescope angular aperture (zenith and azimuth), and the resolution.It uploads the topography data of the geological target, sets the observation point and identifies the projection planes where the telescope is pointing.By ray tracing, it estimates the traversing distances of muons inside the object from all the possible trajectories to the telescope, considering the angular aperture and resolution.Semi-empirical models estimate the differential muon flux depending on the zenith angle and the observation point altitude.Then, an algorithm estimates the minimum energy required by the muon for traversing a given distance inside the target.With this minimum energy and the muon incidence angle, we calculate the integrated muon flux, the telescope exposure time (assuming an expected muon rate per pixel), the target transmission matrix, and the density.
The telescope parameterization module calculates the detector acceptance, solid angle, detection area (from the pixel area), the number of pixels and the separation distance between the detection matrices.This module also estimates the expected muon rate using the integrated particle flux crossing the target and the telescope acceptance.
The muon tomography module reconstructs a threedimensional density distribution of the target from several (simulated or measured) muon radiography projections, having the observation points as inputs.

MUON RADIOGRAPHY
This module performs muon radiography of any geological target by inputting the topography, observation point and projection point, angular aperture and resolution of the detector.In this section, we describe how MUYSC creates the target topography, the muon flux generation by using semi-empirical models, and the muon propagation algorithm.

Topography data
The topography script loads the volcano topography and converts the global coordinates (latitude, longitude, altitude) into a local metric system.The module also locates the observation and projection points.
MUYSC downloads the target topography from the SRTM (Shuttle Radar Topography Mission) NASA dataset at a resolution of 1 arc-second (30 meters) (Earth Resources Observation And Science (EROS) Center( 2017)).The dataset contains the latitude, longitude, and altitude per point.The script inputs the target area coordinates (Longmax, Longmin, Latmax, Latmin).Then, the user establishes the observation point (coordinates of the muon telescope) and the projection point (coordinates where the telescope is pointing).
MUYSC converts global coordinates (Latitude, longitude) to cartesian coordinates (x, y in meters) as follows, where (Long0, Lat0) are the (longitude, latitude) of the reference point, and Ce is the Earth's circumference.

Ray tracing
MUYSC estimates the muon travelling distance along the target for every azimuth and zenith angle depending on the detector aperture and angular resolution.The algorithm traces a line joining the observation point (xop, yop, zop) with the projection point observation point (xpp, ypp, zpp), and obtains the altitude for P points along the projected line.MUYSC evaluates if the altitude of the ith point is above the projected line.If it is, the distance between the (i − 1)th and ith point is added to the total traversing distance.We get a travelling distance matrix for (2N − 1) 2 trajectories, where N is the number of scintillator strips.

Muon flux models
MUYSC estimates the incident muon flux on the target using semiempirical models that reproduce results of Monte Carlo methods such as CRY (Su et al.(2021)).

Gaisser model
The Gaisser parameterization model of the differential muon flux in [cm −2 sr −1 s −1 GeV −1 ] units and is follows (Gaisser(1990); Lesparre et al.(2010)): (3) where AG is a scale factor, γ the power index, BG a balance factor, which depends on the ratio of muons produced by pions and kaons, rc the ratio of the prompt muons produced by the decay of charmed particles, and θ the muon incidence zenith angle.E0 represents the energy of muons at sea level and Ê0 on top of the atmosphere (E0 ≈ Ê0).E0,π and E0,K are interpreted as the critical energies of pions and kaons for vertical incidence (i.e.θ = 0).Most authors consider them adjustable parameters as shown in Table 1.MUYSC uses Gaisser's parameters.

Gaisser-Tang model
In 2006 Tang et al. (Tang et al.(2006)) introduced a modified version of the original Gaisser model taking into account an overestimation of the incident muon flux at low energies (E0 < 100 cos θ GeV): where R Earth = 6370 km is Earth's radius and Hatm = 32 km is the altitude where muons with large angle trajectories are produced (Lesparre et al.(2010)).The muon energy is, where, and,  2021)).The model is expressed in the form of a fitting formula, where y = log 10 p and the muon momentum p in GeV c −1 verifies, where Eµ = 0.10566 GeV c −2 is the muon mass and c = 1.AB , a0 , a1 , a2, and a3 are the fitting parameters, adjusted with different momentum ranges, as shown in Table 2.

Altitude correction
The semi-empirical models presented above estimate the muon flux at sea level.MUYSC corrects the altitude variation of the muon flux using the ratio, h is the altitude of the observation point altitude (in meters), h0 = 4900 + 750p, and p is the muon momentum in GeV (Lesparre et al.(2012); Hebbeker & Timmermans(2002)).

Muon energy loss
The energy loss is modelled as where a and b are functions depending on the material properties through which the muons propagate and (L) is the density integrated along the trajectory.Lesparre et al. (Lesparre et al.(2010)) proposed the following model for the muon energy loss in standard rock based on the Groom dataset (Groom et al.(2001)), Notice that x = log E, with E the muon energy in GeV and l4 = 0.0154, l3 = −0.0461,l2 = 0.0368, l1 = 0.0801, l0 = 0.2549.Fig. 2 shows the performance of the muon energy loss model for energies from 10 −2 GeV to 10 3 GeV.

Minimum energy estimation
The minimum muon energy for traversing a given amount of rock is, where Eµ is the muon mass and the target opacity.The opacity integrates the material density along the muon path.MUYSC estimates the muon minimum energy, Emin = E(a), by solving the following minimization problem where ˆ = ρL is the target opacity with ρ = 2.65 g cm −3 as the rock average density along the muon path L. The ratio is obtained from the fitting model.Figure 3 shows the minimum energy a muon needs to cross 100 m of standard rock.

Traversing muon flux
Once MUYSC gets the minimum muon energy, it estimates the traversing muon flux through the target.MUYSC integrates the differential muon flux in energy from Emin to infinity for the muon zenith angle as follows, and the discrete version, The integrated muon flux I( , θ) depends on the differential

Target density
The integrated target density ρ along the muon path L, leads to the opacity as where κ is the mass attenuation coefficient and T the muon transmission (Peña Rodríguez(2021)).The ratio between the traversing I and the open sky muon flux I0 defines the muon transmission as follows, We get the target density by assuming an average density along the muon path We simulated the Mount Etna muography from the observation point (14.996547• Latitude, 37.742343 • Longitude) with an angular aperture of 50 • azimuth and 20 • zenith.MUYSC estimated a maximum rock thickness ∼ 1.7 km (θ = 82 • , φ = −25 • ), and a muon flux of 2 − 3 × 10 1 [cm −2 sr −1 day −1 ] in the mountain borderland as shown Fig. 5(d).The result coincides with the work of Carbone et al. (Carbone et al.(2013)), but MUYSC resolves muon radiography with a higher resolution.
MUYSC also reproduced results obtained by the Diaphane experiment at the Puy de Dôme volcano (Menedeu(2016)).MUYSC sets the muon detector at the Col de Ceyssat observation point (45.764167• Latitude, 2.955389 • Longitude).We observed a maximum rock thickness of around 1.2 km (θ = 89 • , φ = 3 • ), and a muon flux of 2 − 3 × 10 1 [cm −2 sr −1 day −1 ] at the volcano top as shown Fig. 5(f).MUYSC can also compute the target transmission matrix, the detector exposure time for a given muon rate per pixel threshold, and the target opacity assuming a homogeneous material along the muon path (e.g. standard rock density ρ = 2.65g/cm 3 ).Fig. 6-left shows the exposure time estimated by MUYSC for the muography of the Cerro Machin volcano from the P1 observation point with a uniform detector acceptance of 6 cm 2 sr −1 and a threshold of 10 muon/pixel (Vesga-Ramírez et al.( 2020)).The exposure time in regions with low rock thickness (∼ 10 1 m) spans a couple of days, but in areas where rock thickness increases to 1 km, the exposure time reaches 1000 days.

DETECTOR PARAMETERIZATION
The detector parameterization module computes the angular resolution and acceptance of any muon telescope by inputting the pixel area, the number of pixels per sensitive matrix and the distance between them.
Telescope acceptance, T (rm,n), depends on the detection area, S(rm,n), and solid angle, δΩ(rm,n)) (Lesparre et al.(2010)), as, with A = d 2 the pixel area, d the pixel side, and rm,n the trajectory between the pixels m and n.The trajectory depends on the distance between the detection panels D and the relative distance E between the pixels m and n as shown in Figure 7: ) where i, j are the front panel m-pixel coordinates and k, l the rear panel n-pixel coordinates.
The detection area depends on the number of panel-activated pixels, NP , per trajectory, rm,n, and the pixel area, A, as The muon hodoscope reconstructs 2Ni−1 × 2Nj−1 trajectories, where Ni is the number of i-bars and Nj corresponds to j-bars.
The traversing flux I is defined as where N (rm,n) is the number of detected particles and ∆T is the detection time.2010)).The angular aperture is about ±40 • , the angular resolution has a maximum equal to 19.6×10 −3 sr and the maximum acceptance in r0,0 is 34.57cm 2 sr.This module outputs the telescope's parameters such as the detection area, the angular resolution, the solid angle, and the acceptance.Fig. 9 shows the Muon Telescope, MuTe, parameters for Nx = Ny = 30 scintillator strips, a pixel side d = 4 cm, and a separation D = 250 cm (Peña-Rodríguez et al.( 2020)).The MuTe maximum detection area is ∼14.4 cm 2 , the solid angle for perpendicular trajectories reaches ∼1.02×10 −3 sr, and the maximum acceptance is (∼3.6 cm 2 sr) for r0,0 (Peña-Rodríguez et al.( 2020)).

MUON TOMOGRAPHY
Muon tomography reconstructs three-dimensional density distributions using several two-dimensional muograms (Nagahara & Miyamoto(2018)), estimating 3D density using reconstruction algorithms.A two-dimensional muogram integrates the object density along the muon path but does not distinguish density anomalies along the direction.There exist three families of reconstruction algorithms: analytical (e.g.Filtered Back Projection, FBP), statistical (e.g.Max Likelihood), and algebraic methods (e.g.Algebraic Reconstruction Techniques, ART).Density analysis can be done using techniques such as the inversion of gravimetry data but require a priori geological information or a combination with other geophysical data.Long exposure times, a low muon flux and given topography result in a low number of unequally spaced projections of the object.Such limitations define the appropriate tomography reconstruction approach.
FBP is computationally efficient but requires many equally spaced projections, which are impossible in muon tomography.The low number of projections and long acquisition times limit three-dimensional density reconstruction of geological targets using muon tomography (Barnoud et al.(2021)).Algebraic methods give good reconstruction results with a lower number of projections compared with FBP.
MUYSC handles the muon tomography reconstruction based on the Algebraic Reconstruction Technique.To implement muon tomography reconstruction in MUYSC, we evaluated the performance of various reconstruction algorithms found in the library the TomoPy library (Gursoy et al.(2014)).The test consisted in performing a tomographic reconstruction on the Shepp-Logan object, with 180 observation points located uniformly between [0, π] radians.The results and execution times are shown in A. It is clear from Fig. A1 that the best results are given by ART and Gridrec algorithms.

Algebraic Reconstruction Technique
An ART does not require an initial model as ordinary inversion analysis methods do (Nagahara & Miyamoto(2018)).This method uses a system of linear equations Ax = b where x is the volume vector to reconstruct, given the projection data b, with where A is a sparse m × n matrix whose values represent the relative contribution of each output pixel to different points in the projection data (m being the number of individual values in the projection data, and n being the number of output voxels).As it is shown    The systematic error ε assesses the reconstruction performance as follows (Nagahara & Miyamoto(2018)), εi,j = Ri,j − Ti,j, R is the sum of perpendicular slices across a given volume axis and T is the original projection data on the axis.
For the top view, It is necessary to estimate the reconstruction region.This is done by converting the angular view of the telescope to latitude, Longitude and height.Fig. 11 shows the original and reconstructed volume with normalized units.The error map shows a minimum value (εmin = −0.327) on the top of the secondary dome and a maximum value (εmax = 0.525) around it, while an error on the central dome reaches (εmin = −0.309).Note that the error variation on the main dome is lower than in the secondary.
We calculated the error average µε and the standard deviation σε by taking the absolute value of ε in the volcano region.We got an error average of µ |ε| = 0.105, with a standard deviation of σ |ε| = 0.099.
For the anterior view, we used a projection from the front of the structure (0.017 rad).Fig. 12 shows the original and reconstructed volume in units of rock thickness.MUYSC fully reconstructs the volcano shape.The systematic error increases in the middle of the principal dome and at the top of the secondary one.The error reaches (εmax = 0.37 km), with µ |ε| = 0.158 km and σ |ε| = 0.0768 km.
The three-dimensional reconstruction depends strongly of the number of projections.The lower the number of projections, the lower the reconstruction performance.

CONCLUSIONS
This work describes MUYSC, a software tool to simulate muon radiography and muon tomography of geological objects.MUYSC simulates muography measurements as follows: • the integrated muon flux passing across the geological object  2022)), taking into account the observational conditions (geographical coordinates and detector elevation) reported by the authors.The high computational performance of MUYSC allows it to solve muon radiography (100 × 100 pixels) of a volcanic structure in no more than 2 minutes.
MUYSC also calculates the telescope parameters as follows: • the angular resolution • the acceptance • the distance matrix of particles crossing along the detector The MUYSC tomography module evaluates the viability of carrying out a three-dimensional density reconstruction of a geological object from a given number of muon radiography projections.We presented the reconstruction of the Cerro Machin volcano from 42 muon radiography projections and the systematic error evaluation.
MUYSC runs on computing facilities from personal notebooks to high-performance clusters.The MUYSC code can be easily obtained from a GitHub repository through the Muon Tele-

Figure 1 .
Figure 1.Workflow of the MUYSC framework with three modules: muon radiography, telescope parameterization, and muon tomography.The solid line boxes represent the sub-modules within each module.The red arrows represent intakes, the blue ones indicate the output and the dashed line tomography box estimates each projection.The muon tomography module reconstructs a three-dimensional density distribution of the target from several (simulated or measured) muon radiography projections.This module loads the geological coordinates of the target topography, the observation point coordinates, the projection point coordinates, the telescope angular aperture (zenith and azimuth), and the resolution.The telescope parameterization module calculates the detector acceptance, solid angle, and detection area from the pixel area, the number of pixels and the separation distance between detection matrices.

Figure 3 .
Figure 3. MUYSC estimation of the minimum muon energy for traversing 100 m of standard rock employing the minimization of the opacity error function.The minimum muon energy was ∼ 68.97 GeV.
) 3.9 Validation of the MUYSC muon radiography module We tested the MUYSC's muon radiography module by comparing the target rock thickness and the traversing integrated muon flux estimated by other experiments.MUYSC reproduced the results obtained with MuTe at the Cerro Machín volcano for the observation point P4 (4.494946 Latitude and -75.388110Longitude) with an angular aperture of 30 • azimuth and 25 • zenith (Vesga-Ramírez et al.(2020)).The maximum rock thickness (1.3 km) occurs along the intersection of both volcano domes (θ = 84 • , φ = 16 • to 19 • ), and the muon flux at the top of the dome reaching 10 2 [cm −2 sr −1 day −1 ] as shown Fig. 5(b).

Figure 6 -
(right)  shows the transmission matrix of Mount Vesuvius.Muon transmission is the ratio between the traversing and open sky muon flux.The muon transmission reaches 5×10 −2 at the volcano border and decreases with the rock thickness till 1 × 10 −3 at the volcano crater.Such MUYSC results agree with the work of Alessandro et al.(D'Alessandro et al.(2018)).

Figure 6 .
Figure 6.(left) Expected exposure time of the Cerro Machin volcano from the P1 observation point (4.492298 • Latitude, -75.381092 • Longitude) for an acceptance of 6 cm 2 sr −1 and a threshold of 10 muon/pixel.Expected muon transmission of Mount Vesuvius.

Figure 7 .
Figure 7. Tracking of muons crossing the muon hodoscope.The activated pixels m(i, j) and n(k, l) defines the rm,n trajectory between the hodoscope layers.The trajectory distance depends on the separation between the layers D and the relative distance between pixels E.

Figure 8 .
Figure 8. Muon telescope acceptances.(Top-left) Gilbert et al. propose a telescope with a maximum acceptance in r 0,0 is 64 cm 2 sr.(Top-right) Carbone et al. detector having a full acceptance in r 0,0 gis 5.54 cm 2 sr.(Middle-left) Lo Presti et al. telescope with an angular aperture of ±45 • , and a full acceptance in r 0,0 is 1.04 cm 2 sr.(Middle-right) Lesparre et al. design a telescope with a maximum acceptance in r 0,0 is 12.1 cm 2 sr, and an angular aperture of ±34 • .(Bottom-left) After changing D = 80cm, the maximum acceptance in r 0,0 is 25 cm 2 sr, and an angular aperture of ±45 • .(Bottom-right) Uchida et al. propose a telescope with a maximum acceptance in r 0,0 is 34.57cm 2 sr, and an angular aperture is ± 40 • .

Figure 9 .
Figure 9. Parameterization of the Muon Telescope for Nx = Ny = 30, a pixel side d = 4 cm, and a separation D = 250 cm.The MuTe maximum detection area is ∼ 14.4 cm 2 , the solid angle for perpendicular trajectories reaches ∼1.02×10 −3 sr, and the maximum acceptance is ∼3.6 cm 2 sr for r 0,0 .

Figure 10 .
Figure 10.Muon tomography based on an Algebraic Reconstruction Model.ART reconstructs a voxel (x k ) of the three dimensional structure from the projection data pixels (b i ) recorded by three different muon detectors.

5. 2
Three-dimensional reconstruction with MUYSCWe tested the MUYSC tomography module using simulated data from the MUYSC muon radiography module.It allows us to generate projections of the target from several points, such as to evaluate the 3D reconstruction performance.MUYSC simulated a muon tomography of the Cerro Machín Volcano with 42 observation points at 1059 m ±40 m. from the interesting point with a mean altitude of 2495 m a.s.l.±10 m.All points were located in topography-suitable places around the volcano dome.Each data projection (100 × 100 pixels) contains information on the muon path in a zenith range of 28 • and azimuth range of 100 • .The tomography module of MUYSC gives a threedimensional density distribution of the Cerro Machin volcano (100 × 100 × 100 voxels).
• the minimum exposure time to get given muon per pixel threshold • the geological object transmission matrix • the geological object opacity We reproduced the results of different muography experiments using the MUYSC framework.MUYSC got rock thickness and integrated muon flux traversing the following volcanoes: Cerro Machín volcano (Vesga-Ramírez et al.(2020)), Mount Etna (Carbone et al.(2013)), the Puy de Dôme (Menedeu(2016)), and the Mount Vesuvius (D'Errico et al.(

Figure 11 .
Figure 11.Top view of the original structure (top), the reconstructed volume (middle), and the systematic error map (bottom).The topography cut ranges 4.479 • to 4.494 • Latitude and -75.396 • to -75.381 • Longitude.The secondary volcano dome reconstructed by MUYSC got a higher systematic error variance (0.825) in comparison with the principal reconstructed dome.

Figure A1 .
Figure A1.Reconstruction of Shepp-Logan object with different algorithms.

Table 1 .
Parameters of the Gaisser's model for muon generation determined by several authors.

Table A1 .
Execution time of TomoPy reconstruction algorithms.