Jesús Aranda

Jorge A Pérez

Camilo Rueda

Frank D Valencia

Stochastic Behavior and Explicit Discrete Time in Concurrent Constraint Programming ⋆

We address the inclusion of stochastic information into an explicitly timed concurrent constraint process language. An operational semantics is proposed as a preliminary result. Our approach finds applications in biology, among other areas.

Motivation. The study of quantitative information within languages for concurrency has recently gained a lot of momentum. In many applications, quantitative information becomes crucial when refining models with empirical data, and is of the essence for verification purposes. Two main models of quantitative information can be singled out from the vast literature on the subject. Given a computation that can perform different, competing actions, a probabilistic model provides a probability distribution over such actions. In contrast, a stochastic model relates each action to a random variable which determines its duration: given a set of competing actions, the fastest action (i.e. the one with the shortest duration) is executed. Consequently, notions not considered in a probabilistic model (e.g. speed) are fundamental in a stochastic setting. Not surprisingly, areas in which time is essential (e.g. systems biology, performance modeling) have found in languages featuring stochastic information adequate frameworks for analysis.

Concurrent constraint programming (CCP) [START_REF] Saraswat | Concurrent Constraint Programming[END_REF] is a declarative model for concurrency with strong ties to logic. In CCP, systems are described by pieces of partial information called constraints. Processes interact in a shared store; they either add new constraints or synchronize on the already available information. Timed concurrent constraint programming (tcc) [START_REF] Saraswat | Foundations of timed concurrent constraint programming[END_REF] is a declarative framework for reactive systems. In tcc, time is explicitly represented as discrete time units in which computation takes place; tcc provides constructs to control process execution along such units. In the light of stochastic models for quantitative information, the explicit time in tcc poses a legitimate question, that of determining to what extent the notions of stochastic duration and of discrete time unit can be harmoniously conciliated within a CCP-based framework. The question is relevant because it can give clues on clean semantic foundations for quantitative information in CCP, which in turn, should contribute to the development of more effective reasoning techniques over reactive systems in many emerging applications. In this paper, we outline preliminary results on an operational semantics for a tcc language with explicit stochastic durations.

More into details, the proposed semantics aims at an explicit account of stochastically derived events using the description power of timed CCP calculi. This is a feature that in other CCP calculi (e.g. [START_REF] Bortolussi | Constraint-based approaches to stochastic dynamics of biological systems[END_REF]) is handled at best implicitly. We define stochastic events in terms of the time units provided by the calculus: this provides great flexibility for modeling and, as mentioned before, it allows for a clean semantics. Most importantly, by considering stochastic information and adhering to explicit discrete time, it is possible to reason about processes using quantitative logics (both discrete and continuous), while retaining the simplicity of calculi such as ntcc [START_REF] Nielsen | Temporal concurrent constraint programming: Denotation, logic and applications[END_REF] for deriving qualitative reasoning techniques (such as denotational semantics and proof systems). We consider existing qualitative reasoning techniques have a great potential for guiding/complementing the use of (usually costly) quantitative ones. Such an approach for applying qualitative techniques has shown to be useful in the biological context [START_REF] Fages | Formal cell biology in biocham[END_REF].

This work is part of a larger research programme aimed at developing robust CCPbased techniques for analyzing complex applications and systems in computer music, security and biology. As such, it is our objective to formalize stochastic information in tcc in such a way that resulting languages and techniques (i) remain generic enough so to fit well in the target applications, and (ii) be amenable to efficient implementations, in the form of e.g. simulators and model-checkers.

Description. We consider a variant of tcc in which certain processes are annotated with a function λ, which represents the stochastic information in the language (see below). Annotated processes are tell, when and unless. With a slight abuse of notation, in tell and unless processes λ also stands for the constant value 1. We annotate unless as we see it as a counterpart of when processes. A careful definition of unless in the stochastic context, however, is yet to be completely determined. We do not discard that different applications (e.g. biological systems and computer music) need different unless definitions.

P, Q ::= tell λ (c) | when c do (P, λ) | P Q | local x in P | !P | next (P ) | unless λ c next (P )
Operational Semantics. We use the same notion of discrete time as in ntcc and tcc. We assume that there are discrete time units of uniform size, each of them having its own constraint store. At each time unit, some stimuli are received from the environment; the process then executes with such stimuli as input. At the end of the time unit, some output is produced in the form of responses to the environment, and a residual process to be executed in the next time unit is scheduled. Information does not automatically transfer from one time unit to the following.

The operational semantics, given in Table 1, is defined over process-store configurations. We use γ, γ ′ to range over configurations, and assume a structural congruence relation ≡ to identify processes with minor syntactic differences. The rules of the semantics carry both a probability value (denoted p) and a global rate value (denoted r). They decree two kinds of process execution, immediate (probability value equal to 1 and rate value max), and stochastic. In this sense, processes can be either immediate or stochastic. The idea of the semantics is to schedule immediate processes first, and then move to stochastic processes, whose execution involves a certain duration.

Rules for immediate execution resemble analogous rules in tcc and ntcc. The rule IMMTELL adds a constraint to the store as soon as possible. The rule IMMREP specifies that process !P produces a copy P at the current time unit and then persists in the next time unit. There is no risk of infinite behavior within a time unit. In the Rule IMMUNLESS, process P is precluded if c is entailed by the current store d. The rule IMMINT allows for compositional extension.

Rules for stochastic executions consider the aforementioned function λ. Using the current store as parameter, λ describes how the global rate of the whole process varies. We use δ m (P ) to denote a delay process P with duration m: P will be executed at the m-th time unit from the current one. Given probability and rate values for a process, function ∆ determines its duration. The duration can be thus seen as an exponentially distributed random variable that depends on a probability and a rate.

The rule STOTELL defines stochastic tell actions. The rule STOCHOICE defines a choice over a number of guarded processes. Only those enabled processes, i.e., those whose guards entail from the current store, are considered. The rule STOINT defines the simultaneous occurrence of stochastic actions. As usual, the probability value is calculated assuming independence of the actions. Notice that the current store is not affected by stochastic actions; their influence is only noticeable in the following time units. The rules STOUNLESS and STOREP define unless and stochastic replicated actions, resp. The rule NEXT extends stochastic actions to next processes. In the rule LOCAL, local in P behaves like P , except that all the information on x produced by P can only be seen by P and the information on x produced by other processes cannot be seen by P . Notation (local x, c) P expresses that c is the local information produced by process local x in P . The rule STRCONG is self-explanatory.

These rules define behavior within a time unit; internal behavior takes place until reaching a configuration where no further computation is possible (quiescence). We need to define the residual process to be executed in the following time unit. We start by conjecturing that each quiescent configuration γ has a "standard" form:

γ ≡ j∈J next (P j ) k∈K unless c k next (Q k ) i∈I δ mi (P i ), d .
In the following definition we use A to denote the set of delayed processes in a quiescent configuration.

Definition 1 (Future function) Given a quiescent configuration γ, its residual process is given by function F :

F (γ) = j∈J P j k∈K Q k F ′ (A) IMMTELL tell1(d), c -→1,max skip, c ∧ d IMMREP P, c -→1,max P ′ , c ′ ! P, c -→1,max P next (! P ), c ′ IMMUNLESS unless1 c next (P ), d -→1,max skip, d if d |= c IMMINT P, c -→1,max P ′ , c ′ P Q, c -→1,max P ′ Q, c ′ STOTELL tell λ (d), c -→ 1,λ(c) δ m (tell(d)), c with m = ∆(1, λ(c)) STOCHOICE P i∈I when ci do (Pi, λi), c -→p,r δ m (Pj ), c if c |= cj with r = P i∈{j | c|=c j } λi(c); p = λj (c)/r; m = ∆(p, r). STOINT P, c -→p 1 ,r 1 P ′ , c Q, c -→p 2 ,r 2 Q ′ , c P Q, c -→ p ′ ,r ′ P ′ Q ′ , c with p ′ = p1 × p2; r ′ = r1 + r2. STOUNLESS unless λ c next (P ), d -→p,r δ m (unless c next (P )), d with m = ∆(p, r). STOREP P, c -→p,r δ m (P ), c ′ ! P, c -→p,r δ m (P ) next (! P ), c ′ NEXT P, c -→p,r P ′ , c P next (Q), c -→p,r P ′ next (Q), c LOCAL P, c ∧ ∃xd -→p,r P ′ , c ′ (local x, c)P, d -→p,r (local x, c)P ′ , d ∧ ∃xc ′ STRCONG γ1 -→p, r γ2 γ ′ 1 -→p, r γ ′ 2 ifγi ≡ γ ′ i (i ∈ {1, 2})
Table 1. Operational semantics: internal transition rules.

where function F ′ is defined as

F ′ (δ m1 (P 1 ) . . . δ mn (P n )) = G(δ m1 (P 1 )) . . . G(δ mn (P n ))
and where G is defined as

G(δ m (P )) = δ m-1 (P ) if m > 1 P if m = 1.
Unlike other languages like the stochastic π-calculus [START_REF] Priami | Stochastic pi-calculus[END_REF] or sCCP [START_REF] Bortolussi | Constraint-based approaches to stochastic dynamics of biological systems[END_REF], it is worth noticing that in our semantics stochastic actions can evolve simultaneously; there is no a predefined order for execution. This way, for instance, tell λ1 (c 1 ) tell λ2 (c 2 ) evolves into δ m1 (tell(c 1 )) δ m2 (tell(c 2 )) and in the next unit time, the configuration is δ m1-1 (tell(c 1 )) δ m2-1 (tell(c 2 )) (assuming m 1 ,m 2 > 0). This allows to naturally represent the evolution of different components in parallel.

Discussion. Since variables in tcc are logic (i.e. they can be defined at most once in each time unit), a potential source of inconsistencies is the simultaneous execution of several stochastic actions involving the same variables. This could represent a limitation in modeling. Consider for instance the kind of systems in which it is required to deal with quantities of elements of a certain type (as in biological reactions). In such systems, variables could be part of several actions, which would represent the changes over the elements in consideration. An inconsistency caused by two actions simultaneously altering the value of the same variable is clearly an undesirable feature. Therefore, there is the need for enhancing the semantics with a mechanism that imposes some kind of order over those actions related with potential inconsistencies. This would also presuppose modifications over rules calculating duration of stochastic actions, as concurrent actions would be simulated in a specific order. The formal definition of such a consistency mechanism is part of ongoing work.

Applications in Biology. We think that our language and semantics have applications in the biological domain. This is supported by the fact that CCP-based calculi have shown to be convenient for modelling, simulating and verifying several kinds of biological systems [START_REF] Gutiérrez | Timed concurrent constraint programming for analysing biological systems[END_REF][START_REF] Arbeláez | Timed Concurrent Constraint Programming in Systems Biology[END_REF][START_REF] Bortolussi | Constraint-based approaches to stochastic dynamics of biological systems[END_REF]. In [START_REF] Bortolussi | Constraint-based approaches to stochastic dynamics of biological systems[END_REF], stochastic concurrent constraint programming (sCCP) is used to model biochemical reactions and gene regulatory networks. Functional rates in sCCP give considerable flexibility to formulate reactions. However, sCCP does not include an explicit notion of time and does not exploit the logic nature of CCP for verification. Also, sCCP lacks a means of expressing absence of information, which has proven most useful in the biological context [START_REF] Arbeláez | Timed Concurrent Constraint Programming in Systems Biology[END_REF]. The explicitly timed ccp language ntcc [START_REF] Nielsen | Temporal concurrent constraint programming: Denotation, logic and applications[END_REF] provides both a proof system and a means of representing absence of information. In fact, ntcc was used in [START_REF] Gutiérrez | Timed concurrent constraint programming for analysing biological systems[END_REF][START_REF] Arbeláez | Timed Concurrent Constraint Programming in Systems Biology[END_REF] to model different biological systems using two kinds of partial information: behavioral (e.g. the unknown relative speeds on which a system evolves) and quantitative (e.g. the set of possible values that a variable can take). It must be noticed that ntcc does not allow for stochastic or probabilistic information.

Based on the above, we think that the extension to tcc here proposed could serve several purposes in the biological context. The most immediate use is the definition of enhanced models of systems already modeled in ntcc (the Sodium-Potassium pump, regulation and mutation processes in genetic regulatory networks). Also, although it is not evident that every sCCP process can be translated into our language (the tell operator in sCCP has continuation), we are confident we can model most of the biological systems described in [START_REF] Bortolussi | Constraint-based approaches to stochastic dynamics of biological systems[END_REF]. We also plan to analyse the model in [START_REF] Cardelli | A process model of rho gtp-binding proteins in the context of phagocytosis[END_REF], which describes the cycle of Rho GTP-binding proteins in the context of phagocytosis.

⋆ Research partially supported by the COLCIENCIAS project REACT (No. 1251-330-18902) and the INRIA Équipe Associée FORCES. The work of Jesús Aranda has been supported by COLCIENCIAS (Instituto Colombiano para el Desarrollo de la Ciencia y la Tecnología "Francisco José de Caldas"), INRIA Futurs and ÉGIDE (Centre franc ¸ais pour l'accueil et les échanges internationaux).


R Bajou

M Rosas-Carbajal

A Tonazzo

J Marteau

High-resolution structural imaging of volcanoes using improved muon tracking

Keywords:

Muography is increasingly used to image the density distribution of volcanic edifices, complementing traditional geophysical tomographies. Here we present a new muon data processing algorithm, and apply it to a new generation of scintillator-based muon detectors, to image the relative density distribution in La Soufrière de Guadeloupe volcano (Lesser Antilles, France). Our processing method iteratively searches for the best fit of each muon trajectory, accounting for all the hits registered by the detector related to the particular muon event. We test the performance of our algorithm numerically, simulating the interaction of muons with our detector and accounting for its exact assemblage including the scintillator bars and lead shielding. We find that our new data processing mitigates the impact of spurious signals coming from secondary particles, and improves the amount of successfully reconstructed events. The resulting two-dimensional muon images at La Soufrière have higher angular resolution than previous ones and capture the heterogeneous structure of the dome. They show density anomalies located on the summit southern region, which includes a boiling acid lake and degassing fractures, where the rock is the most porous and fumarolic activity is ongoing. This work shows the importance of combining numerical simulations of muon propagation with precise raw data processing to obtain high-quality results. It is also a first step towards fully as-sessing the noise contamination sources when performing muon tomography, and their correction, prior to geophysical interpretations.

INTRODUCTION

Volcanoes are complex systems that vary significantly in terms of size, internal structure, and the dynamics of heat and mass transport from the magma chamber to the surface. These properties thus require dedicated studies for each particular volcano in order to understand their functioning and mitigate their potential hazards. Both active and dormant volcanoes carry potential edifice collapse hazards associated with their geometry and past activity, which usually leads to mechanically weak rocks and collapse-prone flanks [START_REF] Heap | The mechanical behaviour and failure modes of volcanic rocks: a review[END_REF].

Geophysical imaging is used to study the volcano's internal structure and is performed with different methods such as seismic, gravity, magnetotellurics, electrical resistivity tomography, and muography. Each of these methods is sensitive to the spatial distribution of a particular physical property, but they all provide complementary information useful to 1) determine regions of different rock types, in particular places where the rock has been altered or is naturally mechanically weak;

2) understand the magmatic and meteoric fluid pathways, the potential chemical reactions, and the heat transport dynamics inside the volcano; 3) relate processes detected in the surface to the main geological structures; and 4) calibrate numerical models meant to explain the volcano functioning.

When information about possible unstable and/or altered rock is available, failure surfaces can be suggested to test numerically the velocity and final emplacement of debris avalanches [START_REF] Peruzzetto | Modeling of partial dome collapse of la soufrière of guadeloupe volcano: implications for hazard assessment and monitoring[END_REF].

Muography is a relatively novel method used to image the density distribution inside volcanic edifices [START_REF] Tanaka | Japanese volcanoes visualized with muography[END_REF][START_REF] Marteau | Muon tomography applied to activevolcanoes[END_REF][START_REF] D'alessandro | Volcanoes in italy and the role of muon radiography[END_REF]. Its imaging principle relies on particle absorption and deviation in a dense medium. The natural muon flux resulting from the interaction of cosmic rays with the atmosphere arrives at the Earth's surface and decreases as the muons traverse a dense structure. The decrease in the muon flux is a function of the mean density of the crossed medium along the muons' trajectory. A muon detector measures the flux of muons exiting the volcano along different directions or lines of sight. Thus, two-dimensional (2-D) muonbased images (hereafter called muographies) of the mean density in the different lines of sight can be obtained.

Originally used to characterize the natural rock shielding overlying underground particle physics laboratories in the 1950s, muon absorption imaging was then transported towards archaeological studies and turned into an in situ experiment, by [START_REF] Alvarez | Search for hidden chambers in the pyramids[END_REF], during an imaging campaign performed on the Chephren pyramid of the Giza ancient Egyptian complex. This pioneer field measurement and feasibility demonstration lead the path towards the applications that are used these days in volcanic structure imaging, with the first measurements led by Japanese teams in the 1990s [START_REF] Nagamine | Method of probing inner-structure of geophysical substance with the horizontal cosmic-ray muons and possible application to volcanic eruption prediction[END_REF]. Repeating muon measurements at several positions around a volcano, the final goal is to combine the 2-D mean density images to compute a three-dimensional (3-D) density model, as reported by [START_REF] Rosas-Carbajal | Three-dimensional density structure of La Soufrière de Guadeloupe lava dome from simultaneous muon radiographies and gravity data[END_REF][START_REF] Nishiyama | 3d density modeling with gravity and muon-radiographic observations in showa-shinzan lava dome, usu, japan[END_REF][START_REF] Barnoud | Robust Bayesian Joint Inversion of Gravimetric and Muographic Data for the Density Imaging of the Puy de Dôme Volcano (France)[END_REF]. The analogy with medical X-ray CT-scan is thus often made in the literature. The reader is invited to refer to [START_REF] Bonechi | Atmospheric muons as an imaging tool[END_REF] for an exhaustive review of the muon imaging physics and detection technology, and to [START_REF] Lechmann | Muon tomography in geoscientific research -a guide to best practice[END_REF] for a specific state of muography applications in geosciences.

La Soufrière de Guadeloupe is a volcano located in the Lesser Antilles belonging to France's overseas territory and classed as one of the most dangerous in terms of its active hydrothermal system [START_REF] Siebert | Volcanoes of the World[END_REF]. Not only responsible for sudden violent eruptions without magmatic juvenile material involved (the last major crisis occurred in [1976][1977], hydrothermal activity [START_REF] Komorowski | Guadeloupe[END_REF][START_REF] Tamburello | Spatio-temporal relationships between fumarolic activity, hydrothermal fluid circulation and geophysical signals at an arc volcano in degassing unrest: La soufrière of guadeloupe (french west indies)[END_REF][START_REF] Moretti | The 2018 unrest phase at la soufrière of guadeloupe (french west indies) andesitic volcano: Scrutiny of a failed but prodromal phreatic eruption[END_REF][START_REF] Moune | Gas Monitoring of Volcanic-Hydrothermal Plumes in a Tropical Environment: The Case of La Soufrière de Guadeloupe Unrest Volcano (Lesser Antilles)[END_REF] and the resulting rock alteration process plays a major role in the destabilization of the structure, which worsens the risk of slope failure and partial flank collapse (Heap et al., 2021;[START_REF] Ball | The hydrothermal alteration of cooling lava domes[END_REF]. In addition to established geophysical methods such as seismic tomography, gravimetry surveys, muon absorption imaging, referred to in the literature as muography, or simply muography, has now become part of the imaging tools available to map out density profiles of geological bodies such as volcanic dome, and identify density variations relative to acid hydrothermal fluid circulation.

Since the early 2010s, scintillator-based particle trackers, referred to in the following as muon telescopes, were developed at IP2I, Lyon, and deployed on the field at different locations around La Soufrière de Guadeloupe flanks. The goal was to perform dynamic muon absorption imaging of the andesitic lava dome. These campaigns have set the ground for inserting muon imaging into the scope of conventional monitoring at La Soufrière, as a non-invasive tool. Conjointly with other geophysical and geochemical methods, they have indeed enlightened the hydrothermal system dynamics inside the volcanic edifice and its effects on the rock density (Lesparre et al., 2012a;[START_REF] Jourde | Muon dynamic radiography of density changes induced by hydrothermal activity at the La Soufrière of Guadeloupe volcano[END_REF]. Three-dimensional density models have been inferred, notably by joint inversion of muon and gravimetry datasets by [START_REF] Rosas-Carbajal | Three-dimensional density structure of La Soufrière de Guadeloupe lava dome from simultaneous muon radiographies and gravity data[END_REF], putting forward several low-density anomalies on the southern volcano flank, presumably in the areas where the rock is highly porous and mechanically weaken by the acid fluid circulation. Furthermore, as scintillator detectors offer time information on the muon arrival, continuous monitoring of sudden destabilization in the hydrothermal activity has been explored on time scales of hours and days [START_REF] Gonidec | Abrupt changes of hydrothermal activity in a lava dome detected by combined seismic and muon monitoring[END_REF] through the combination of muon data with seismic signal backed up with vent temperature measurements, highlighting strong To this day, muon telescope data at La Soufrière were reconstructed using a fast and simple tracking principle: a straight-trajectory check of the most-energetic hits on the 3 scintillator panels. But muons are not the only particles generating a signal in scintillator panels, and often many hits per panel are observed in what is registered as a single "event". Muons are indeed expected to form a straight track in the detector, and they are minimum ionizing particles, which implies that their associated energy deposit in the scintillator medium -of the order of 2 MeV/(g.cm -2 ) -is not the highest possible on each panel. Therefore this kind of selection is already biased by the presence of other types of particles generated in atmospheric showers. Furthermore, depending on the relative angle between the muon trajectory and the detector orientation, a single muon can traverse more than one scintillator "pixel" per panel. A solid fit procedure of the hits is thus necessary to retrieve the correct muon trajectory and cope with contamination from other atmospheric particles and spurious hits formed by instrumental noise (e.g. optical or electrical cross-talks, dark counts arising in time coincidence).

In this paper, we present a new method to process the raw muography data, which improves the reconstruction of detected muons. Our tracking algorithm is based on a random sampling consensus procedure RANSAC [START_REF] Fischler | Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography[END_REF], which provides the best fit to the particle trajectory con-sidering all the panel hits involved in a single event. To evaluate the performance of our algorithm, we use the GEANT4 Monte-Carlo simulation framework for particle propagation [START_REF] Agostinelli | Geant4-a simulation toolkit[END_REF], and assess the reconstruction efficiency. Finally, we apply our algorithm to novel data acquired at La Soufrière de Guadeloupe by two muon detectors. We provide the relative density muographies and interpret them in terms of the main geological features observed at the volcano.

METHODOLOGY

Muon detection

Scintillator-based particle detectors are well-established tools for muography [START_REF] Tanaka | Development of a two-fold segmented detection system for near horizontally cosmic-ray muons to probe the internal structure of a volcano[END_REF][START_REF] Marteau | Muons tomography applied to geosciences and volcanology[END_REF][START_REF] Saracino | The muraves muon telescope: technology and expected performances[END_REF][START_REF] Lo Presti | Muographic monitoring of the volcano-tectonic evolution of Mount Etna[END_REF][START_REF] Peña-Rodríguez | Design and construction of mute: a hybrid muon telescope to study colombian volcanoes[END_REF]. With their robust, modular, and relatively simplified design compared to high-energy physics detectors, they are appropriate to deal with the main challenge of rough field operations (Lesparre et al., 2012a).

Since 2017 two muon telescopes, respectively equipped with four and three scintillator panels, are installed at two different locations around the La Soufrière lava dome: at the South-West, and East of the volcano whose positions are shown in Fig. 1.

The novel 4-panel configuration (4p), in addition to decreasing the low-energy particles background, allows increasing spatial resolution on the 2D density muography of the central part of the dome scanned region and offers two more 3-panel sub-configurations (3p1 and 3p2), which lead finally to three different muon images.

As represented in Fig. 1, the 4-panel telescope named SNJ allows scanning at the same time a broad angular region using signal recorded in three panels with respective azimuth and zenith angles (ϕ 3p , θ 3p ) ∈ ([-17; 58] • , [43; 107] • ), as well as a narrower angular region centered on the dome bulk and its two main craters, using all 4 panels: (ϕ 4p , θ 4p ) ∈ ([-5; 46] • , [52; 98] • ). The two geometrical configurations have thus different angular acceptances. SNJ includes 4 detection matrices (front, middle1, middle2, and rear) interspaced by 60 cm; each matrix is an array of 16x16 bars, which is 80 cm long, 1 cm thick, and 5 cm wide.

The 3-panel telescope, named BR, is located on the Eastern flank. It has 3 matrices (front, middle and rear) also interspaced by 60 cm: the front and rear matrices are composed of 32x32 bars, 80 cm long, 1 cm thick, and 2.5 cm wide, while its middle matrix ensures signal coincidence without affecting the total number of lines-of-sight and thus presents a lower pixel resolution with a 16x16 bar configuration, similar to SNJ matrices. All the matrices have the same detection area of 80x80 cm 2 .

Each scintillator bar has its own wavelength-shifter (WLS) fiber embedded within the scintillator medium to conduct the light towards a multi-anode photomultiplier (MAPMT), with 64 channels.

One MAPMT per detection matrix takes care of the light signal conversion to electrons and ampli-fication before processing by the readout electronics described in [START_REF] Marteau | Muons tomography applied to geosciences and volcanology[END_REF]. A simple straight-line interpolation of the signals in the first and last plane would provide an angular resolution of 12 mrad (3p) or 8 mrad (4p) for SNJ and of 6 mrad for BR.

In addition to the detection matrices, a 100 mm thick lead shielding panel has been installed on both telescopes behind the middle matrix, compatible with mechanical constraints on the structure. This passive shielding allows for stopping low-energy particles. We estimated its cutoff energy to be around 170 MeV, using the numerical simulations explained later on in this work. The shielding also triggers electromagnetic showers, composed of secondary electrons and photons, characterized by a high hit multiplicity on the rear panel. Detailed descriptions of the muon telescope assembly used in this study can be found in Lesparre et al. (2012a) and [START_REF] Jourde | Experimental detection of upward-going cosmic particles and consequences for correction of density radiography of volcanoes[END_REF]; and description of the light and electronic readout as well as the timing system in [START_REF] Marteau | The OPERA global readout and GPS distribution system[END_REF], [START_REF] Marteau | Muons tomography applied to geosciences and volcanology[END_REF][START_REF] Jourde | Experimental detection of upward-going cosmic particles and consequences for correction of density radiography of volcanoes[END_REF].

Rock density estimation

In muon absorption tomography, the observable is the number of reconstructed particle tracks N detected by the muon telescope during a period of time ∆T , for each line-of-sight r i,j where (i, j) represents a given front-rear pixels combination. Assuming a point-like telescope in comparison to the size of the scanned object, each (i, j) corresponds to a unique (θ, ϕ) angular direction. The detector is characterized by a direction-dependent experimental acceptance T exp , which accounts for the telescope geometry, as well as instrumental and reconstruction effects. To estimate T exp , the telescope is placed in open-sky mode during a calibration run, in a vertical position, that is, the central axis at zenith θ = 0 • , and the following quantity is derived:

T exp (r i,j ) = N sky (r i,j ) ∆T calib × I sky (r i,j ) [cm 2 .sr] (1)
where N sky (r i,j ) is the number of open-sky muons detected in r i,j , ∆T calib is the calibration run duration (around two days). The open-sky muon integral flux I sky (r i,j ) is estimated with CORSIKA [START_REF] Heck | CORSIKA: A Monte Carlo Code to Simulate Extensive Air Showers[END_REF], a Monte-Carlo simulation tool that allows accurate cosmic muons flux estimation.

CORSIKA propagates extensive atmospheric showers, from their generation by primary cosmic-ray, e.g. an extra-galactic proton, interactions on upper atmosphere nuclei to their ground-level detection. Furthermore, this detailed MC framework allows us to consider geomagnetic effects alongside atmospheric parameters (e.g. density variations through a layered atmosphere model) that affect the muon flux, as studied in [START_REF] Cohu | Atmospheric and Geodesic Controls of Muon Rates: A Numerical Study for Muography Applications[END_REF].

Once the acceptance has been calculated, the telescope is oriented towards the volcano and the transmitted muon flux noted I data can be estimated for each line-of-sight r i,j with the following expression:

I data (r i,j ) = N (r i,j ) ∆T × T exp (r i,j ) [cm -2 .s -1 .sr -1 ] (2)
where N (r i,j ) is the number of transmitted muon candidates along r i,j , ∆T is the tomography run duration, and T exp (r i,j ) is the experimental acceptance in r i,j , defined in Eq. 1.

In order to infer density values from I data (r i,j ), a comparison with reference muon flux after attenuation in rock is necessary. The reference transmitted flux noted I calc ( , θ) depends on the zenith angle θ and on the crossed amount of matter, quantified in muography with a parameter called opacity , defined as the integral of density ρ along rock thickness L:

(L) = ˆL ρ(ξ) dξ [g.cm -2 ] (3)
As for the experimental acceptance computation (cf. Eq. 1), I calc ( , θ) is estimated from Monte-Carlo simulation e.g. CORSIKA, computationally intensive but more accurate at the observation site for the reasons mentioned before than semi-analytical models (e.g. [START_REF] Gaisser | Cosmic Rays and Particle Physics[END_REF][START_REF] Gonidec | Abrupt changes of hydrothermal activity in a lava dome detected by combined seismic and muon monitoring[END_REF]) based on the muon production from mesons in the high atmosphere. A computation for each opacity of the muon minimal crossing energy E min is required to obtain the transmitted muon flux for a known density distribution. For this purpose, the work of the Particle Data Group [START_REF] Workman | Review of Particle Physics[END_REF]) is widely used. It provides an exhaustive review of the different interaction processes at stake, and their associated energy deposit contributions, during muon propagation in a dense medium. Finally, the muon energy spectrum Φ(E, θ) can be integrated over [E min ( ); +∞[:

I calc (θ, ) = ˆ+∞ E min ( ) Φ(E, θ) dE [cm -2 .s -1 .sr -1 ] (4)
Once the opacity is estimated from data-model flux comparison, a mean density value can be obtained for each line-of-sight r i,j : ρ i,j = i,j /L i,j .

Data processing

Along with the energy deposits (so-called hits) generated by muons crossing the detector and used to reconstruct their trajectories, other sources induce signals in the scintillator bars: other atmospheric particles from electromagnetic and hadronic components of atmospheric showers, secondary background particles (electrons and photons) produced inside the detection volume, or instrumental noise caused notably by random emission of electrons from the MAPMT photo-cathode and cross-talks. To mitigate their impact on the track reconstruction process, we present a new track reconstruction algorithm that considers the presence of these "outlier" hits. The tracking algorithm was written in Python, with the use of the ransac module available in the scikit-image library (van der [START_REF] Van Der Walt | scikit-image: image processing in Python[END_REF]. The fit procedure is also useful in the context of the 4-panel telescope configuration, to account for detection inefficiencies: for muons forming hits in only 3 panels, but contained in the 4-panel angular region, the trajectory model allows to extrapolate the muon position in the untouched panel.

Our event-by-event track fitting procedure lies on an hypothesize-and-test iterative process, called random sample consensus RANSAC [START_REF] Fischler | Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography[END_REF]. The different steps of the reconstruction workflow are represented in Fig. 2. First, we filter out events presenting a number of hits per panel larger than 10, which are likely to be due to an electromagnetic shower, or impacting less than 3 detection panels. The process starts by randomly sampling within the hit coordinates, a subset of points of a given size, and fitting a hypothesized line model to this subset. It then evaluates which hits of the whole dataset are consistent with this hypothesized model, computes a fit error, and counts the number of hits whose orthogonal distance to the model is below a given distance threshold t (set to one scintillator width). Those hits form the consensus set. The size of the initial sample s and the maximal number of iterations N max were set to maximize the number of muon hits tagged as inliers, using the simulation of an open-sky muon run with the framework described in Section 2.4. After N max iterations, the algorithm outputs the model that minimizes the fit error and maximizes the size of the consensus set. The hits belonging to the consensus are tagged as 'inliers' and the rest of the hits as 'outliers'. This method manages to increase the number of reconstructed muons and to improve the angular precision on the incident particle direction, compared to the former reconstruction algorithm used, relying on the check of the alignment of the most-energetic hits, i.e a check of the "straightness" of the track. A comparison of the two methods has been studied in [START_REF] Bajou | A New Versatile Method for the Reconstruction of Scintillator-Based Muon Telescope Events[END_REF].

Numerical simulations

Developing simulations for particle detectors in muography, as in all particle physics experiments, represents a crucial step to a better understanding of the detector and its response to particles of interest. More precisely, it allows estimating detection efficiency, a key parameter for quantifying the performance of our tracking algorithm. The most widespread tool for this purpose, the GEANT4 toolkit [START_REF] Agostinelli | Geant4-a simulation toolkit[END_REF] offers a framework with an exhaustive set of physics models to simulate particle interactions and tracking through a detailed numerical replica of the detector. The simulated signals can then be formatted as real data to test the reconstruction algorithm. Examples of simulated 10 GeV/c electron and muon events are shown in Fig. 3. The electron event in Fig. 3a is filtered out thanks to the 7-cm thick shielding panel between the middle 1 and middle 2 panels. The muon track in Fig. 3b is reconstructed with RANSAC. Along with the physical muon hits, Fig. 3b features outlier hits that were randomly added to the muon event signal. These "synthetic" hits were not assigned any intensity value (i.e. deposited energy) as this parameter is not considered in the current tracking workflow. This reconstructed event display illustrates the ability of the RANSAC procedure to fit the correct muon trajectory from a data set polluted by hits generated by instrumental effects like a deficient optical coupling between optical fibers, dark current arising from the photo-cathode and dynodes in MAPM, electronic cross-talk between MAPM channels.

RESULTS

Reconstruction performances

In the following section, we characterize the performance of our muon tracking algorithm for a 4panel telescope. For this purpose, we use the GEANT4 simulation setup described in Section 2.4. We simulate an open-sky acquisition and write the muon hits recorded in the scintillator panels following the real telescope raw data format. Then, we use our algorithm to reconstruct the particle tracks following the reconstruction outline depicted in Fig. 2. Because we know the exact particle trajectories from the GEANT4 simulations, we can assess the tracking algorithm's performance. The results presented hereafter were obtained on a pure-muon sample of 10 7 events simulated in GEANT4 with energies E µ ∈ [1, 10 6 ] MeV. Within this sample, we estimated that 7% of the generated muon tracks contained in the 4-panel acceptance cone are being absorbed in the detector shielding, before im- pacting the rear detection matrix. The RANSAC tracking algorithm allows extrapolating the muon intersection position on this rear matrix.

Influence of secondary particles

Secondary electrons and photons produced by high-energy muons inside the detection volume, notably after crossing the shielding panel, impact the quality of the reconstruction of the incident muon trajectory. The RANSAC tracking allows for mitigating the impact of those spurious hits. In Fig. 4 we show the comparison between the real particles obtained in GEANT4 and the hit tagging performed by the algorithm. For initial sample size s = 2 and N max ≥ 20 the true muon hit inlier fraction reaches a plateau at 98.4 ± 0.2%. Besides, the muon tracking result also shows relatively low contamination from electrons and photons: these secondary hits represent 14.3 ± 0.7% of the total hits, and nearly half of this fraction (6.7 ± 0.6%) have been tagged as outliers.

Reconstruction efficiency

To assess the reconstruction efficiency ε reco of our tracking algorithm, we computed for each pixel corresponding to a given line of sight r i,j , the ratio between the number of RANSAC-reconstructed primary muons and the number of generated primary muons in GEANT4 :

ε MC reco (r i,j ) = N reco (r i,j ) N 0 (r i,j ) (5)
where N reco is the number of reconstructed muons, and N 0 is the number of primary muons.

The reconstruction efficiency obtained with the GEANT4 simulation ε MC reco (cf. Eq. 5) can be compared with the reconstruction efficiency estimated with open-sky calibration data, which we defined as the following:

ε DATA reco (r i,j ) = T exp (r i,j ) T geo (r i,j ) (6)
that is the ratio between the experimental detector acceptance T exp (cf. Eq. 1), and the geometrical acceptance T geo , whose analytical expression is only a function of the telescope geometrical parameters: distance between the front and rear detection matrices, number, length, and width of the scintillator bars in each matrix [START_REF] Sullivan | Geometric factor and directional response of single and multi-element particle telescopes[END_REF][START_REF] Thomas | Analytical derivation of the geometric factor of a particle detector having circular or rectangular geometry[END_REF]. Examples of efficiency maps obtained with the simulation and with real open-sky calibration data are shown in Fig. 5.

The reconstruction efficiency depends on the detection efficiency of the scintillator bars ε scint , that is, the probability for a given hit in a scintillator bar to give a signal above the threshold. We assume ε scint to be the same for all X and Y bars over each panel, although small differences from one bar to another occur in reality due to instrumental effects such as a deficient optical coupling between fibers and the light readout system (cf. Sec. 2.1). In Fig. 6, we show how the mean efficiency ε MC reco as a function of the ε scint parameter. From this relation, and using the mean data reconstruction efficiency ε DATA reco , it is straightforward to interpolate the real telescope scintillator efficiency. We consider here the acceptance measurements performed on a 4-panel telescope (see Fig. 1b). As a result, we obtained the following scintillator efficiencies for the 3-panel and 4-panel events respectively: ε 3p1 scint = 0.899 ± 0.012, ε 3p2 scint = 0.900 ± 0.010 and ε 4p scint = 0.868 ± 0.015. These are consistent with the instrumental efficiency of the order of 0.9 usually considered for the technology of plastic scintillator bars with optical fiber readout employed in our detectors (Lesparre et al., 2012a).

Tomography data analysis

We present in this section the results obtained by applying our tracking algorithm (cf. Sec. 2.3) to real muon tomography data acquired by two muon telescopes (cf. Sec. 2.1): SNJ and BR, respectively located at South-West and East of the La Soufrière lava dome (cf. Fig. 1a).

The SNJ tomography dataset contains 1.9 × 10 7 reconstructed events acquired over a total run duration of 145 days in steady operation conditions (excluding dead time superior to one hour), between February and August 2019. The BR acquisition period ran from February 2018 to August 2019 over 547 steady days and contains 1.2 × 10 7 events. Both tomography runs were preceded by a 2-day long open-sky calibration run. Following the methodology explained in Sec. 2, the detector experimental acceptance is first estimated from calibration runs, then the incoming transmitted muon flux for each detector can be calculated, including the 3-and 4-panel configurations of detector SNJ.

To quantify the agreement between the data-estimated and modeled transmitted fluxes for different density values we use a Pearson χ 2 -test :

χ 2 (θ, ρ) = N los i,j I data (r ij ) -I calc (θ, ρL ij ) σ I (r ij ) 2 (7)
where N los is the total number of lines of sight at a given θ, and σ I is the statistical uncertainty on the estimation of I data .

We have chosen to present the following analysis results for both telescopes in terms of relative mean density variations rather than absolute density whose estimation suffers from irreducible background pollution (see Sec. 4). We then consider the quantity ∆ρ defined with respect to the average value ρ 0 = ρ from each telescope muography: ∆ρ = (ρρ 0 )/ρ 0 . Figure 7 shows on a given azimuth range [ϕ c -2.5, ϕ c + 2.5] • (where ϕ c is the central telescope azimuth) how the χ 2 varies as a function of ∆ρ and zenith angle θ for each muon telescope. The chosen azimuth range targets for both tele- scopes the lava dome region beneath the Tarissan crater (TAR on Fig. 1a, roughly 100 m deep from the dome surface which represents a vertical angular range of ∆θ 10 • ). The observed variations reflect the sensitivity to spatial variations of density. On the chosen azimuth range, we observe rather constant relative density in Fig. 7a located between -14.3% and -4.6% on the whole θ range. On the contrary, BR in Fig. 7b presents a stiff increase in ∆ρ from -50.0% for smaller θ up to +18.2% at θ = 86.5 • , which directly echoes the increase in density from the shallow region close to the surface to the deeper part of the dome.

We now focus on the mean density spatial variations. Since the SNJ detector consists of 4 panels, we show in Fig. 8 the results for the 3-and 4-panel estimates. The measurements obtained by this telescope give evidence of significant density variations. The lowest values around -60% with respect to the mean ρ 0 , are found for line of sights along the west flank for azimuth ϕ ≤ 10 • close to the surface, and for the highest zenith values, corresponding to horizontal lines of sight in the east side (-30%). Besides these extreme values found on the borders, we see inside the dome a tendency from west to east from lower to higher average densities. The largest density values are indeed observed on the east side of the dome reaching up to +70% for azimuth ϕ ≥ 40 • .

The 4-panel image shown in Fig. 8 b) is a more detailed, "zoomed" version of Fig. 8 a). We see that both muographies are consistent, although the muon events used for the 4-panel acquisition naturally comprise fewer muon events than the 3-panel radiography. Since this radiography provides a better spatial resolution in the zone of interest, we present this case in Figs. 8 c) and d), a satellite picture of the lava dome showing the points where the different lines of sight intersect the surface.

The color coding corresponds to the density radiography and the thickness of rock traversed for each line of sight, respectively. When observing these figures, it should be noted that the relative density variations correspond to an average value through the whole trajectory inside the volcano, and not to the point where the circle is represented. We can observe in Fig. 8 c) that the high-density anomaly thus corresponds to lines of sight that have traversed the dome flanks in a region directly south of the summit, usually referred to as the "bulge" in the literature [START_REF] Boudon | A new scenario for the last magmatic eruption of la soufrière of guadeloupe (lesser antilles) in 1530 a.d. evidence from stratigraphy radiocarbon dating and magmatic evolution of erupted products[END_REF][START_REF] Rosas-Carbajal | Volcano electrical tomography unveils edifice collapse hazard linked to hydrothermal system structure and dynamics[END_REF]Heap et al., 2021). This region is thought to consist of relatively fresh and dense rocks. Focusing on the lines of sight close to the South Crater and the Tarissan acid pond, we can interpret from Figs.

8 b) and c) that Tarissan is located in a region of overall lower bulk density than the South-Crater (CS). Finally, if we compare the relative density distribution to the rock thickness distribution, it is interesting that the color spatial distributions are not the same, suggesting that the density anomalies observed do not correspond to changes in the amount of rock traversed by the muons.

Figure 9 shows the equivalent tomography and bird view representations for the BR detector.

Although this is a 3-panel detector, the spatial resolution is increased with respect to SNJ because the scintillator bars are smaller. The largest density anomalies are observed for θ < 80

• and φ ∈ [260, 290],
that is, for lines of sight that cross the same "bulge" region as in the SNJ radiography (cf. Fig. 9 b)). We also retrieve on this East muography an overlaying low-density region above zenith 70 • . The lowdensity anomaly also seems to rise up to the surface below the South-Crater region, around azimuth 290 ± 4 • . On this muography, though, the shallow low-density anomaly appears broader and spans over the whole azimuth range above zenith 75 • , reaching the lowest values of -70% above azimuth 310 • , and zenith 65 • . This low-density region close to the surface corresponds to a superficial part of the dome on the north-west flank, which crosses highly-altered and clay-rich material, revealed during the 2009 landslide [START_REF] Rosas-Carbajal | Volcano electrical tomography unveils edifice collapse hazard linked to hydrothermal system structure and dynamics[END_REF]Heap et al., 2021). The main pits and craters, whose positions are projected on the BR muography, are part of a wide network of cavities spanning from the north fault to the Tarissan crater and expanding on the East flank. These structural features are expected to have a significant impact on the bulk density estimates.

DISCUSSION

In this section, we discuss the results obtained regarding the performance of our particle tracking algorithm, the implications of the density radiographies provided, and possible improvements in future works.

Our RANSAC-based muon tracking algorithm was already shown to outperform the straightness check tests used previously by our team in terms of the number of events reconstructed and the precision on the inferred trajectory [START_REF] Bajou | A New Versatile Method for the Reconstruction of Scintillator-Based Muon Telescope Events[END_REF]. Here, we did a more exhaustive performance test using GEANT4 simulations to test our algorithm against perfectly known particle trajectories.

Moreover, our numerical simulations consider the secondary particles produced as a result of the interaction with the lead panel located at the center of our detector (cf. Fig. 4). The processing of the GEANT4-generated data is encouraging in the sense that our algorithm successfully tagged true muon hits as inliers and significantly mitigated secondary hits. The comparison of real open-sky acquisitions with simulated data in equivalent conditions provides comparable results (cf. Fig. 6) and allows us to estimate the scintillator bars' efficiency. Our tracking algorithm is flexible in that it can be adapted to different geometries of scintillator panels.

The main drawback is that it takes more time to run that the straightness check routines: O(10 3 s) for 10 6 events on an 8-core CPU. Still, these are not prohibiting computation times considering that the analyses so far are not performed in real-time, and the significant improvement in data quality they provide.

The density spatial variations inferred from the real data acquisitions in La Soufrière de Guadeloupe can be qualitatively compared with a high-resolution electrical resistivity tomography of the dome performed by [START_REF] Rosas-Carbajal | Volcano electrical tomography unveils edifice collapse hazard linked to hydrothermal system structure and dynamics[END_REF], keeping in mind that rock conductivity is also influenced by other parameters such as rock porosity and fluid content. This work highlighted the presence of two distinct electrical conductivity regions: a low-conductivity shallow region named "bulge" on the south flank composed of non-altered dense andesitic rock, overlying a highly conductive region composed of mechanically weakened rock by hydrothermal activity. This bulge region scanned by both telescopes appears as the reddish high-density regions on the three muographies, for ϕ > 30 • in Fig. 8a and 8b, and ϕ < 290 • in Fig. 9a. On the other hand, the observed low-density region on the south flank is thought to correspond to the superficial part of the high-conductivity anomaly rising up to the dome summit.

As of today, this high electrical conductivity anomaly is understood to be the main pathway for hydrothermal acid fluids circulation, connecting southern thermal springs at the base of the lava dome to the dome summit. Surface manifestation of the arrival of these fluids includes continuous emission of sulfur and chlorine-rich gas through the fumarolic vents around Southern and Tarissan craters (cf. Fig. 1a) [START_REF] Moune | Gas Monitoring of Volcanic-Hydrothermal Plumes in a Tropical Environment: The Case of La Soufrière de Guadeloupe Unrest Volcano (Lesser Antilles)[END_REF]. The continuous circulation of these fluids greatly contributes to the alteration and mechanical weakness of the dome rocks (Heap et al., 2021). Mechanical weakness is in part related to an increase in porosity, and thus it is reasonable that our muographies observe low-density values in highly altered regions as compared to fresh rock zones.

Among the structural features, several cavities of large volume inside the lava dome could have a significant impact on the observed densities. The BR muography featuring a broad low-density anomaly in Fig. 9 seems particularly impacted by the presence of such a network of cavities. In addition to the Tarissan and South craters, two additional pits of large size are located on the East-South-East flank: Breislack (BLK) and Gouffre 56 (G56), roughly 40 m deep (10% of the dome height), according to the description made by [START_REF] Kuster | Guadeloupe, Canyons, Gouffres[END_REF], around azimuth 305 • and 312 • . Another important cavity, the Spallanzani cave of volume 55x35x11m 3 [START_REF] Kuster | Guadeloupe, Canyons, Gouffres[END_REF], whose entrance is located in the north fault (FNO), is also likely to have an impact on the mean density measurement. The effect on the density of such volume was already estimated between 10 to 15% by Lesparre et al. (2012b). As these cavities are currently not taken into account when computing the total dome thickness traversed by the muons, shown in Fig. 8d 9c, these void spaces add to an average bulk density significantly lower than the average value of actual rocks. The western superficial low-density anomaly that appears in Fig. 8a and b may also be caused by the existing network of cavities, expanding to the west. Nevertheless, the cavities delimitation is not well-documented in this area of the dome. As this low-density anomaly concerns the most superficial part of the dome and thus the shortest thicknesses L (cf. Fig. 8d) it is also likely to be affected by a systematic error on telescope orientation, whose impact on density estimation is important on this thinnest part of the dome. For the shallowest lines of sight of Fig. 8a and 8b covering a band of 3 pixels wide below the topographic profile (black line), we estimate that a 1-degree shift in both telescope azimuth and elevation gives a mean relative error ∆L/L = 0.21 ± 0.12, while for the rest of deeper lines of sights, this error reaches down to 0.09 ± 0.04. In the future, we intend to further develop the simulation tools to reproduce the geometry of the main cavities and provide a better characterization of their impact on the mean density estimates.

The telescopes are ideally located to scan what is thought to be the most-altered and collapseprone part of the lava dome on the south flank. A downside of their closeness to the volcano flanks is that for lines of sight close to the horizon, that is, for large zenith values, there is a high probability of exposition to a particular source of noise contamination. This contamination comes from irreducible background sources such as the forward scattering effect of low-energy muons.

Previous studies [START_REF] Nishiyama | Monte carlo simulation for background study of geophysical inspection with cosmic-ray muons[END_REF]Gómez et al., 2017;[START_REF] Chen | Monte Carlo simulation for background noise study of the Wudalianchi volcanic density imaging[END_REF] have indeed enlightened that for such large targets as volcanoes, flux measurements are strongly impacted by lowenergy muons that are forward-scattered on the volcano flank facing the detector. These "soft" muons typically in the sub-GeV energy range can mimic dome-outgoing muon tracks and be reconstructed as such. This results in a measured flux that is over-estimated, which in turn leads to a lower opacity and thus lower average density estimate along the impacted line-of-sight. Indeed, this can be clearly observed in Fig. 8a and b for zenith angles higher than 85 degrees. Rock thicknesses in this area are so large that no muon events should be detected, therefore the low-density estimates arise from a large amount of forward scattered muons. Another source of contamination lies in the muon flux produced in the atmosphere below the telescope altitude. This backward muon flux can enter the telescope by its rear matrix, which especially impacts the most horizontal lines of sight above θ > 80 • and can significantly blur muon radiography close to the horizon, accounting for up to 50% of the detected events [START_REF] Jourde | Experimental detection of upward-going cosmic particles and consequences for correction of density radiography of volcanoes[END_REF]. The ideal rejection for this upward flux requires an ns clocking system for a precise time-of-flight measurement.

Gómez et al. ( 2017) have shown for a simplified cone-shaped dome geometry that the forward scattering of muons increases with the zenith angle θ, exceeding 50% of the total number of reconstructed events for θ superior to 85°. This background noise is also location-dependent and thus requires specific estimation for each muon imaging study. To cope with this irreducible background noise and infer realistic absolute density values, the simulation setup presented in Sec. 2.4 is expected to provide the right framework to accurately estimate the percentages of forward-scattered to outgoing muons. We leave these developments for a future contribution.

CONCLUSION

We presented a new muon tracking algorithm to process data from 3-and 4-panel scintillator-based muon detectors. Our algorithm successfully captures muon hits and rejects outlier hits formed by secondary particles in the telescope data. We have also developed simulation tools that show that our new algorithm outperforms previously used ones based on straight-line analyses. Furthermore, we can use these simulations to reconstruct our scintillator's efficiency and study the experimentally calibrated acceptance against the simulation predictions.

Our algorithm was successfully applied to estimate muon density radiographies of La Soufrière with data acquired by a new generation of muon telescopes. The density variations have an improved spatial resolution with respect to previous works and are overall in agreement with previous studies.

Although the absolute mean density estimates are biased by irreducible background sources, the relative mean density radiographies still reflect the heterogeneous structure of the dome, thanks to long and steady runs of data acquisition obtained in the two different detectors.

Future work will focus on correcting the two-dimensional muographies from systematic effects such as muon forward scattering, and in the following, on performing data inversion on the corrected density maps to obtain a realistic three-dimensional density model of the dome. The forward scattering correction will be a natural improvement using the GEANT4 simulations developed in this study.

The topography of the terrain surrounding each detector has key importance in this matter. It can now be reproduced with the implementation of the digital terrain model in the presented GEANT4 framework. As an alternative to GEANT4, we also mention the PUMAS library [START_REF] Niess | The PUMAS library[END_REF]) that could also be used for scattering correction studies. This code is based on backward Monte-Carlo to perform muon transport in matter and was initially developed for muography applications. Furthermore, implementing the geometry of the main known geological structural features of the lava dome, fractures, faults, and craters is also possible and will greatly improve the quality of our radiographies and therefore of the three-dimensional density models obtained with them.

Figure 1 .

 1 Figure 1. (a) Map of La Soufrière lava dome featuring locations of the two telescopes SNJ and BR and their respective angular acceptance. The SNJ telescope has a 4-panel configuration and covers two regions at the same time: a broader region in light yellow shade with its 3-panel sub-configuration, and a tightened region using the 4-panel configuration. The main structural features, notably the Tarissan (TAR) and South (CS) craters, and active fumarolic vents on the dome summit are encompassed by both telescope angular apertures. (b) Photo of the SNJ telescope during its assembly, located at the South-West of the dome, composed of 4 detection matrices and an additional lead shielding panel at the center.

Figure 2 .

 2 Figure 2. Data processing algorithm flowchart based on the RANSAC method.

Figure 3 .

 3 Figure 3. Three-dimensional representations of the hits formed in the telescope scintillator panels by: (a) an electromagnetic shower, triggered by the interaction of the primary electron with the shielding panel located between the Middle1 and Middle2 panels; (b) a simulated muon reconstructed with the RANSAC tracking algorithm, artificially polluted with synthetic outlier hits. The color scale represents the sum of energy deposits in X and Y bars, noted dE. The size of the XY hit is also proportional to dE. This information is however not used in the tracking algorithm.

Figure 4 .

 4 Figure 4. Results of RANSAC tracking on simulated muon sample: (a) fraction of muon hits tagged as inliers as a function of the maximal number of iterations N max , for different sizes of initial sample s; (b) composition of the XY hits formed by muons µ, and secondary particles: electrons e and photons γ, resulting from the interaction of muons with the detector. The inner circle represents the recorded hit composition from GEANT4 particle identification. The outer circle represents the proportion of hits tagged by the tracking algorithm as inliers and outliers. The RANSAC parameters have been set to t = 50 mm, N max = 20, and s = 2.

Figure 5 .

 5 Figure 5. Reconstruction efficiency maps for a 4-panel detector estimated from (a) Monte-Carlo simulated muons with GEANT4 reconstructed with RANSAC considering a scintillator efficiency of ε scint = 0.87, and (b) a real 4-panel detector data.

Figure 6 .

 6 Figure 6. Mean reconstruction efficiency ε reco as a function of scintillator efficiency parameter ε scint , obtained from the muon MC simulation in GEANT4 (colored bands). The efficiency is derived for each of the three 4panel telescope configurations: 3p1, 3p2, and 4p. The figure features calibration open-sky data points for each of the configurations. The hatch regions correspond to the intersection of the MC efficiency bands with the data uncertainty bands of the same configuration.

Figure 7 .

 7 Figure 7. χ 2 distribution (cf. Eq. 7) between estimated and calculated fluxes as a function of (θ, ∆ρ) on a 5 • -wide azimuth range around telescope central azimuth, for (a) SNJ and (b) BR datasets.

Figure 8 .Figure 9 .

 89 Figure 8. Relative density variation two-dimensional maps estimated with the SNJ telescope (South-West) for (a) the three-panel configuration of SNJ (b) the four-panel configuration of SNJ. The map in (b) is projected onto the terrain model in (c): each line of sight (black line) features its density value, plotted at the entry point in the dome, and corresponds to a given dome thickness represented in (d). The main dome craters and pits positions of Fig. 1 are represented in blue lozenge markers.

ACKNOWLEDGMENTS

This work and the data shown were acquired in the context of two projects supported by the French Agence Nationale de la Recherche, projects DIAPHANE (ANR-14-ce04-0001) and MEGaMu (ANR-19-CE05-0033). We thank the editor C. Sens-Schönfelder, reviewer R. Nishiyama, and an anonymous reviewer, whose comments have contributed to improving the quality of our manuscript. The simulation and data analysis were performed on the S-CAPAD/DANTE platform of IPGP, France.

DATA AVAILABILITY

The data presented in this work are not publicly available.