Reconstruction of air shower muon lateral distribution functions using integrator and binary modes of underground muon detectors

The investigation of cosmic rays holds signiﬁcant importance in the realm of particle physics, enabling us to expand our understanding beyond atomic conﬁnes. However, the origin and characteristics of ultra-high-energy cosmic rays remain elusive, making them a crucial topic of exploration in the ﬁeld of astroparticle physics. Currently, our examination of these cosmic rays relies on studying the extensive air showers (EAS) generated as they interact with atmospheric nuclei during their passage through Earth’s atmosphere. Accurate comprehension of cosmic ray composition is vital in determining their source. Notably, the muon content of EAS and the atmospheric depth of the shower maximum serve as the most signiﬁcant indicators of primary mass composition. In this study, we present two novel methods for reconstructing particle densities based on muon counts obtained from underground muon detectors (UMDs) at varying distances to the shower axis. Our methods were analyzed using Monte Carlo air shower simulations. To demonstrate these techniques, we utilized the muon content measurements from the UMD of the Pierre Auger cosmic ray Observatory, an array of detectors dedicated to measuring extensive air showers. Our newly developed reconstruction methods, employed with two distinct UMD data acquisition modes, showcased minimal bias and standard deviation.


Introduction
Cosmic rays constitute a population of highly energetic elementary particles and nuclei with an unknown origin that descend upon Earth from outer space.Their spectrum a e-mail: varada.varma@iteda.cnea.gov.ar(correspondingauthor) follows a nearly power law distribution, spanning from approximately 10 9 eV to 10 20 eV [1].Direct measurement of primary cosmic rays with sufficient flux, which occurs at low energies, is feasible through experiments conducted in space.Nevertheless, for energies surpassing approximately 10 15 eV, the flux weakens, necessitating reliance on the interactions between primary particles and atmospheric molecules to generate secondary particles called extensive air showers (EAS) [2].These showers can be observed during their progression in the atmosphere, either on the Earth's surface or underground.The Pierre Auger Observatory, positioned in the southern hemisphere [3], encompasses detectors capable of investigating cosmic ray showers at all three levels: during their development in the atmosphere, on the surface, and underground.Consequently, these showers are reconstructed to examine the primary particles' three principal observables: energy spectrum, arrival direction, and chemical composition.Notably, the Pierre Auger Observatory employs Underground Muon Detectors (UMD) to directly measure the muon content of the showers.High-energy muons exhibit superior penetration capabilities compared to other secondary particles.Subterranean experiments have demonstrated that the density of muons serve as indications of the primary cosmic ray nuclei's chemical composition and energy spectrum.However, this sensitivity is constrained by a threshold imposed by the thickness of the soil covering.Muons possess a unique sensitivity to composition due to a phenomenon where lighter particles (e.g., protons) exhibit lower efficiency in producing multiple muons when compared to heavier nuclei.Although underground detectors cannot determine the energy and specific type of primary particles on an event-to-event basis, it is possible to derive information about the mass composition by comparing the measured distributions of muon multiplicities with those calculated through precise Monte Carlo simulations employing trial models of the primary spectrum and composition.
The UMD at the Pierre Auger Observatory measures the fall in muon density with increasing distance to the shower axis, known as the muon lateral distribution function (MLDF) [6].Once completed, the UMD will be equipped with a total of 73 muon detectors.Among these detectors, a triangular array with a spacing of 750 m will encompass approximately 20 km 2 of the area, housing 61 muon detectors.The remaining 12 detectors will be arranged in a smaller triangular array with a spacing of 433 m, covering approximately 1 km 2 .For the purpose of this study, we will solely focus on the 750 m array.These detectors within the 750 m array possess the capability to measure showers with energies ranging from around 10 16.5 eV to 10 19 eV and can detect events up to a zenith angle of 45°.Each grid location is equipped with three 10 m 2 modules, segmented into 64 plastic scintillation strips containing embedded wavelength-shifting optical fibers that are coupled to an array of 64 silicon photomultipliers (SiPMs) [7].Collectively, these three modules form a 30 m 2 detector with 192 segments.Data acquisition within the UMD is triggered by the associated Surface Detector (SD) station.To minimize contamination from energetic electrons and gamma particles accompanying muons near the shower core, the detectors are buried at a depth of 2.3 m underground.The soil density at the UMD site measures 540 g cm −2 , corresponding to a shielding equivalent to 22 radiation lengths.This ensures absorption of the electromagnetic component of the Extensive Air Shower (EAS), enabling only muons with energies above approximately ∼ 1 GeV/ cos θ µ (where θ µ denotes the zenith angle of muon motion) to reach the buried detectors.
The UMD encompasses two distinct acquisition modes, specifically the binary mode and the integrator mode, also referred to as the Analog to Digital Converter (ADC) mode.These modes operate simultaneously to provide a broad range of detection capabilities [7].Each mode employs distinct methodologies to transform raw signals into the number of muons.
The counting of muons is accomplished by the binary mode segment-wise, which tallies signals exceeding a certain amplitude threshold [8].Each of the 64 SiPM signals is processed independently by a pre-amplifier, a fast shaper, and a discriminator located in every channel of two 32-channel ASICs [9].Subsequently, the signal from the discriminator undergoes sampling into 64 2048-bit traces with a 3.125 ns sampling time, accomplished via a Field-Programmable Gate Array (FPGA).The binary mode has a total inhibition window of 37.5 ns (12 samples × 3.125 ns) [11].This method of muon counting is independent of the gain or fluctuation of the optical device, is independent of the muon impact position, and does not require a thick scintillator.However, it tends to count two muons that arrive simultaneously on the same strip as one because of the pile-up effect.This statistical error can be corrected as long as the number of strips with simultaneous signals is lower than the total segmentation [8].The binary mode saturates if the majority of strips in a given module simultaneously receive a signal, and it can only detect a limited number of muons at the same time, which limits its ability to probe at distances close to the shower axis.The accuracy of the detector resolution in binary mode is also limited by geometric biases such as the corner clipping bias.Overcounting occurs when a muon passes through the lateral edge of a scintillator bar and is detected in both bars, producing the corner clipping bias [12].As more particles hit the detector, the resolution of the binary mode decreases, and the densities closer to the shower cannot be accurately determined.
The ADC mode determines the charge of all segments and converts them to the number of muons.This is done by dividing the total signal charge by the average charge of a single vertical muon signal [13].During ADC acquisition mode, the 64 SiPM signals are added up analogically, and the resulting sum is amplified using low-gain and highgain amplifiers with an amplification factor ratio of about 4. These two channels are called Low Gain (LG) and High Gain (HG) channels, respectively.The signals are sampled every 6.25 ns (twice the binary sampling period).The fluctuations in the ADC mode reduce with the number of detected muons.However, as it depends on the estimation of particles from an integrated signal, signal fluctuations propagate to uncertainties in the estimation of the number of muons.The ADC mode has a resolution of 60% for single-muon signals, which decreases inversely with the square root of the number of particles.It reaches 10% after several tens of muons, matching the resolution of the binary mode.With increasing muons, the ADC mode provides much better resolution than the binary mode, allowing for better precision in measuring particle densities closer to the core [7].The operating range of the LG and HG channels in the ADC mode overlaps with the binary mode and depends on the muon density [13].
Another experiment that measures the muon content of the air shower is the Akeno Air Shower Observatory.It consisted of a 100 km 2 air shower array called the Akeno Giant Air Shower Array (AGASA) [4].The experiment completed operations in 2004 and later merged with the High Resolution Fly's Eye (HiRes) group, forming the Telescope Array Collaboration [5].AGASA has played a vital role in providing critical information about cosmic rays.Surface and muon detectors with concrete shielding were installed in AGASA.AGASA utilized proportional counters for measuring the muon density through two acquisition modes [14].The first method, known as the on-off density method, determines the density by counting the number of hit counters based on the assumption that the number of particles incident on each counter follows a Poisson distribution around the mean value.This approach is comparable to the binary mode of the UMDs located at the Pierre Auger Observatory.The second method, called the analogue density method, calculates the density by dividing the total energy losses in all counters by the average energy loss of a vertically traversing muon.This technique is similar to the ADC method mentioned earlier in the UMDs located at the Pierre Auger Observatory.We illustrate the reconstruction methods proposed based on the UMDs of the Pierre Auger Observatory, but note that these methods can be extended to any muon detector with similar configurations.
The size of an air shower in surface arrays is characterized by fitting the MLDF.This can be achieved by maximizing a likelihood function and evaluating it at a reference distance [15].The likelihood used is specific to each detector and is based on the response of the detector to incoming particles.For the UMDs at the Auger Observatory, a global estimator of the shower muon density is obtained by reconstructing the MLDF and evaluating it at 450 m.This reference distance of 450 m was proved to be the optimal distance for an array of detectors spaced 750 m apart.This means that muon density fluctuates the least at 450 m, compared to other reference distances, for different samples of the same simulated muon LDF [17].In the past, reconstruction methods depended only on information from the binary mode.The time-independent likelihood method is the most straightforward method to reconstruct the MLDF with binary mode data [15].Stations are divided into three different classes based on the number of activated segments, and an approximate likelihood is used for each class.This single-window likelihood method can reconstruct the MLDF with an array of segmented detectors without considering detector timing, using the exact likelihood of the number of muons in a detector given the number of segments with a signal [16].The time-dependent reconstruction method, on the other hand, uses detector timing and a profiling technique to reconstruct the MLDF from binary mode data [17].This paper presents a hybrid reconstruction method suitable for a detector with both a binary mode and an ADC mode.Additionally, it proposes a new reconstruction method that uses only the ADC information.
Section 2 describes the characterization of the output of the ADC.In Section 3, the simulation of the Monte Carlo shower library and the ADC mode is discussed.Section 4 outlines the likelihoods that were utilized for the reconstruction process.The performance of the two newly introduced methods is discussed in Section 5, along with a comparison with the previous reconstruction methods.The paper concludes in Section 6.

Characterization of the ADC output
The ADC output is linear and can be expressed as the arithmetic sum of individual signals [13].The overall charge associated with n muons can be represented as, where q i represent the charge of a single muon.The number of muons that hit the detector (n) follows a Poisson distribution (P(n|µ)) with parameter µ.
The average number of muons in equation ( 2) is given by µ = ρ µ A cos θ , where ρ µ is the mean muon density of an air shower on the shower plane, A the active area of the UMD module, and θ the zenith angle of the air shower.The probability density function that characterizes the total output of the detector can be expressed as, where f (Q|n) corresponds to the total charge distribution of n muons.The mean and standard deviation of Q are, Here q is the mean charge corresponding to a single incident muon and ε[q] corresponds to the relative error of the charge deposited by a muon, which is the ratio of the standard deviation to the mean (ε[q] = σ [q]/ q ).The estimator of the number of muons µ = Q/ q is an unbiased variable with variance σ As stated in Ref. [6], the uncertainties in estimating the muon density of air showers are dominated by Poisson fluctuations in the muon content.Notably, equation (5) highlights that the ADC contributes to the uncertainty of an ideal Poissonian detector, as the fluctuations in µ are greater than those expected from the Poisson distribution.As mentioned in the introduction, the resolution of the binary mode plateaus as more particles hit the detector, while the resolution of the ADC is better at higher particle densities.
To provide a clearer explanation of the detector's resolution, we compare the resolution of the ADC mode (σ [ µ]/µ) with the binary mode and an ideal Poisson detector in Figure 1.In order to calculate the resolution of the binary mode, we utilize the single window (also known as the infinite window strategy), which was developed in Ref. [16].As depicted in the plot, the binary mode demonstrates similar performance to that of an ideal Poisson detector at low numbers of muons.However, as the number of muons increases, the ADC mode begins to outperform the binary mode.For the ADC mode of the UMDs at the Pierre Auger Observatory we consider the charge of a single muon to be a random variable whose logarithm is normally distributed and the probability density function of charge can be represented using a 2-parameter log-normal distribution of the form, where m is the scale parameter and θ is the shape parameter, and can be calculated in terms of q and ε[q] as, If the shape parameter remains small (θ 2 1), which is the case for the UMDs at the Pierre Auger Observatory, the sum of n log-normal terms is expected to behave like a log-normal variable, regardless of the value of n as cited in Ref. [18].When at least one muon is detected (n ≥ 1), the distribution f (Q) for a total charge Q corresponding to n muons is obtained as, The log-normal character of Q is verified using simulations (see Appendix C).
The parameters m n and θ n are estimated from the param- eters corresponding to the distribution function of one incident muon (m and θ ).
In the linear region of the ADC channel, the probability density function of Q can be derived from equations ( 3) and ( 9) and can be expressed as a compound distribution of the form, where δ (Q) represents the Dirac delta function, accounting for the possibility of zero signal due to no muons and Θ (Q) is the Heaviside function.In the worst-case scenario where all muons arrive in the same time bin, a single 10 m 2 module causes the LG channel to saturate at about 362 muons, while the HG channel saturates at a much lower value of around 85 muons [9].As the UMD consists of three such modules, the detection range is even wider.
In regions of higher particle density, i.e., for large values of µ (or µ), it is expected that the probability density function of total charge for a given value of µ can be approximated by a Gaussian distribution (refer to Appendix A).
where Q and σ [Q] can be obtained from equations ( 4) and ( 5) respectively.

Simulations
To generate the Monte Carlo shower library for developing the reconstruction method, we used CORSIKA v7.7100 [19] with the high-energy hadronic interaction model EPOS-LHC [20] and the low-energy hadronic interaction model FLUKA [21,22].We simulated showers of iron and proton primaries with energies in the range of log 10 (E/eV) = [17.5, 19] in steps of log 10 (E/eV) = 0.25 for zenith angles of 30°and 45°.The azimuth angles were randomly generated within the range of −180°to 180°.To reduce the number of tracked particles, a statistical thinning of 10 −6 was applied during simulation.For each energy and zenith angle, we produced 35 proton and 30 iron showers.Although there is experimental evidence of a muon deficit in simulated air showers [23], this does not affect the comparison of different strategies to reconstruct the MLDF.This deficit in simulations seems to be related to the lack of knowledge of hadronic interactions at higher energies, as the high-energy hadronic interaction models used in simulations extrapolate the accelerator data to cosmic ray energies.
From simulated showers, we obtained the average Muon Lateral Distribution Function (MLDF) and the average arrival time distribution of muons as a function of their distance to the shower axis, for each primary energy, zenith angle, and primary.To account for soil shielding effects, we selected muons with energies greater than 1 GeV/cosθ µ .We fitted the average MLDFs with a KASCADE-Grande MLDF [24] to obtain an input MLDF.In Figure 2, we present the arrival time histogram for iron showers at a zenith angle of 30°and energy 10 18.5 eV at three different distances.These histograms show the fraction of particles arriving in 10 ns time bins with respect to the total number of muons.We observed that muons arrive more evenly spaced in time farther away from the shower axis.
The MLDF fit and muon time distribution averages obtained for each primary particle, energy level, and zenith angle serve as input for the detector simulations.The binary mode was simulated as described in Ref. [16] and includes the pile-up effect simulation.The muons are counted within a 40 ns time window.
The simulation program, which was originally developed in Ref. [16], has now been expanded to include the simulation of the ADC.The normalized signal as a function of time generated by a single muon can be characterized as a function resembling a log-normal distribution, denoted as s(t,t 0 ).
where Θ (t) is the Heaviside function, ∆t is the sampling time (6.25 ns for this work), t 0 is the signal start time, m T = ln(τ/ns) and θ T are the parameters of a log-normal- like function, and N is the number of samples of the signal.Note that s(t,t 0 ) has the dimension of [time] −1 .
Once the ADC reaches saturation, the linear proportionality between the signal charge and the number of injected muons is lost.To determine the saturation level of the ADC signal, one can calculate the maximum of the function s(t,t 0 ) for a specific number of muons required to saturate the ADC channel (N µ ).
In this study, we focus solely on the LG channel of the UMD at the Pierre Auger Observatory.The saturation level is determined by the value of N µ = 1086, which represents three times the number of muons required by a single module to reach saturation [9].The parameters of the log-normal-like function are selected as ln(τ/ns) = 4 and θ T = 0.4.These values correspond to a particular set of lab- oratory measurements performed using four 4-m-long scintillator strips inside an aluminum container.Optical fibers were glued to each strip and coupled through an optical connector to a calibrated SiPM array, which consisted of 64 SiPMs connected to the front end of the electronics kit and sealed into a PVC box.A muon telescope was used to trigger the electronics when the muon crossed the four bars.Muon pulses measured with the ADC were recorded at different positions of the muon telescope along the bars, and the average signal as a function of time of one muon was obtained (see reference [10] for more details).Note that any other suitable parameter values for the of the log-normallike function would produce similar outcomes.Detector saturation establishes the upper limit of the ADC detection range.
The ADC signal in units of ADC counts/time corresponding to n muons is obtained by considering the unsaturated signal, where the n, t i , and q i values are obtained by sampling the Poisson distribution, the corresponding average muon arrival time distribution, and the charge distribution of one muon (see equation ( 6)), respectively.Finally, the signal is obtained from the following expression, The values of the parameters corresponding to the single muon charge distribution of equation ( 6) considered are m = 5 and θ = 0.5 [10].In this case as well, these values correspond to the same particular set of laboratory measurements representing the given muon detector.Figure 3 shows an unsaturated and saturated pulse for the aforementioned scenario.
The total signal or charge in a given detector (Q) in equation (1) is obtained by adding the signal evaluated in each discrete time value, obtained according to the sampling time of the ADC, multiplied by ∆t.

Likelihoods and the reconstruction method
This paper examines two reconstruction methods: the ADC reconstruction method, which relies on information from the ADC only, and the ADCProfile reconstruction method, which combines ADC and binary outputs.The ADCProfile reconstruction method is an extension of the Profile reconstruction method that was developed in [17].The MLDF, denoted by µ(r; p), is dependent on two free parameters in p = (µ 0 , β ) and the distance between the detector and the shower axis (r).Here, µ 0 represents the nor- malization parameter, while β is the slope.The MLDF is factorized by µ 0 and a form function g(r; β ) that rapidly decreases with the distance to the shower axis, µ(r; p) = µ 0 g(r; β ) g(r 0 ; β ) , where g(r; β ) is, Fixed values are taken for α = 0.75, r 1 = 320 m, and r 0 is a reference distance, which is taken as 450 m in this work.The γ parameter is also a fixed value and is obtained from the input MLDF.It is almost constant with energy, zenith angle, and primary type.It describes the behavior of the MLDF at large distances to the shower axis, where the statistics are low and the detector has few muons or is non-triggered.In this work, we fix γ at −4.18.The free-fit parameters in p are obtained by minimizing the function, The sum runs over the detectors in this context.To choose L for ADC reconstruction, the estimated number of muons in the LG channel of the ADC is considered.On the other hand, for ADCProfile reconstruction, L selection depends on the estimated number of muons in the various time windows of the binary traces.
To begin with, let us examine the ADC reconstruction.If the estimated number of muons in the ADC channel, denoted as µ, is less than or equal to 3, the station is labeled as non-triggered.This is done to eliminate any accidental triggers that may be caused by random coincidences.In this case the probability is calculated by integrating equation (12) Here Q min = 3 q and erfc(x) = 1 − erf(x), where erf(x) is the error function.
If a station triggers and the estimated number of muons in the LG ADC channel is less than 200, a compound likelihood function is taken into account, When the estimated number of muons in the ADC channel is equal to or greater than 200 and the station is not saturated, a Gaussian likelihood is taken into consideration (refer to Appendix A).This is done to improve the performance of the reconstruction.
where the mean and the variance are given by equations ( 4) and ( 5) respectively.
If a particular station is saturated, we obtain the likelihood by integrating equation ( 23) from Q to +∞.Here Q refers to the integral of a pulse, where the upper portion has been clipped (refer to the bottom panel of figure 3), where the mean and the variance are given by equations ( 4) and ( 5) respectively.Now let us consider the ADCProfile reconstruction method.In this case, the ADC likelihood is utilized when the binary mode of a given muon detector produces a number of bars with a signal (k) greater than or equal to 124 in at least one inhibition window.If k < 124 in all time windows of the trace, the Profile likelihood of Ref. [17] is used.This choice of k is motivated by the fact that the relative uncertainty in estimating the number of incident muons that fall in the same inhibition windows is less than ∼ 16% for k ≤ 124 [27].The non-triggered stations are defined and treated as described in Ref. [17].
To compare the various likelihoods considered in this study, we plotted the logarithm of the likelihoods, as shown in Figure 4. We compared the compound likelihood in equation (22) to a profile likelihood from Ref. [17], which was utilized in the ADCProfile reconstruction approach.The plots corresponds to a value of µ ML = 423, which is the maximum likelihood estimator of µ obtained from Profile likelihood.In this section, we will evaluate the performance of the new reconstructions.As already explained in section 3, to simulate the number of muons arriving at the detector for each shower, we sampled the Poisson distribution with parameter µ given by the average MLDF.The time distributions were sampled to simulate the arrival times of the muons.Each of these realizations is referred to as an event.We simulated 10,000 events for each primary type, zenith angle, and primary energy, including the simulation of the binary and ADC acquisition modes.The number of muons as a function of the distance to the shower axis for simulated events was adjusted using the Profile, ADC, and ADCProfile likelihoods, and then compared against an ideal counter [16].An ideal counter simply counts the number of muons crossing it, which sets the best-case scenario for the resolution achievable with a muon detector.To minimize the loglikelihoods, we used the numerical minimization software library MINUIT [25], which is implemented in the ROOT data analysis framework [26].Figure 5 displays an example of an MLDF fit achieved using the ADC reconstruction method (top panel) and AD-CProfile reconstruction method (bottom panel).The presented data depict the number of muons in the UMD for iron primaries at 10 19 eV and 30°zenith angle, with each station labeled according to the likelihood utilized in the reconstruction process.The top panel plot indicates that Gaussian likelihood (see equation ( 23)) is used for stations near the shower axis, while Compound likelihood (see equation (22)) is used for other stations.The plot also shows a few non-triggered stations and a saturated station.On the other hand, the bottom panel plot demonstrates that the ADC reconstruction method is used for stations near the shower axis, while the Profile reconstruction method is used for other stations.Fig. 5: The muon lateral distribution function was fitted to the simulated data using the ADC reconstruction method (top panel) and the ADCProfile reconstruction method (bottom panel).The data correspond to the estimated number of muons in the UMD calculated with a simulation of an iron primary at 10 19 eV and 30°zenith angle.The solid lines correspond to the MLDF used as the input for the simulation.Note that the plots corresponds to two different events.Figure 6 displays the percentage of saturated events, as a function of the logarithmic energy.An event is categorized as saturated when it involves at least one saturated muon detector.In both the ADCProfile and ADC approaches, detectors are viewed as saturated when the estimated count of muons in the ADC channel surpasses the saturation limit of the LG channel, as described by equation (15).In the Profile reconstruction technique, a detector is regarded as saturated if, in a single time window, a station contains 192 bars with signal.The graph shows a noticeable increase in the percentage of saturated events in relation to primary energy.Moreover, the fraction of saturated events appears to be similar for all the considered reconstruction methods, primarily due to the use of 40 ns wide inhibition windows in binary mode.However, the fraction of saturated events for Profile reconstruction is smaller when the inhibition window is shorter (refer to Appendix B).Fig. 6: Fraction of saturated events for iron primaries at 30°z enith angle as a function of the logarithm of primary energy for ADC, ADCProfile, and Profile reconstruction methods, including the error bars.The error bars are smaller than the marker size.
The reconstructed µ ML (450) (the maximum likelihood es- timator of the average number of muons at 450 m from the shower axis) varied in reconstructions of the same shower due to fluctuations in detector signals.Figure 7 shows the distributions of the reconstructed µ ML (450) for iron primaries at 10 18.5 eV and 30°zenith angle, considering various reconstruction methods and includes the saturated events.The data corresponding to the ideal counter were fitted with a Gaussian distribution, and the dotted line represents the input value of µ(450).As anticipated, the narrowest distribution corresponds to the ideal counter.
To assess the reconstruction performance, it is necessary to compare the input value µ(450) to the corresponding

Counts
Gaussian fit Ideal counter method Profile method ADC method ADCProfile method Fig. 7: A histogram of the reconstructed µ ML (450) for iron primaries at 10 18.5 eV and 30°zenith angle for the different reconstruction methods considered.The Gaussian distribution is fitted to the ideal counter data.
value fitted to the simulated data.The relative bias is defined as the difference between the input µ(450) and the average of µ ML (450) calculated over all reconstructions, i.e., B = µ(450) − µ ML (450) , and represents the systematic uncertainty in the estimation of µ(450).Figure 8 displays the relative bias and relative standard deviation, ε(450) = σ [ µ ML ](450)/µ(450), of µ ML (450) as functions of the log- arithmic energy for different reconstructions, including the saturated events.Here, µ ML and σ [ µ ML ](450) are esti- mated from the Gaussian distribution fit of µ ML (450) his- tograms.The figure shows that the relative bias for the ADC and ADCProfile reconstruction methods is less than 2% at lower simulated energy and lies within 1% for energies larger than 10 17.5 eV.The value of ε(450) decreases with primary energy due to the increase in the number of muons triggering more detectors as well as the smaller relative fluctuations at different detectors.The ε(450) values for the ADC reconstruction are similar to those for the Profile reconstruction.Specifically, ε(450) for the ADC reconstruction is larger than that for the Profile reconstruction at low energies, but at high energies, it is smaller.The ε(450) for the ADCProfile method takes smaller values or is equal to those for the ADC and Profile reconstructions.Note that the ε(450) for the ADCProfile differs by no more than ∼ 2% from that of the ideal counter in the whole energy range considered.
Figure 9 displays the relative bias and relative standard deviation as a function of the distance to the shower axis for iron primaries at an energy of 10 18.25 eV and a zenith angle of 30°.The UMD reference distance is indicated by a dotted line.The ADC and ADCProfile reconstruction methods exhibit a bias similar to that of an ideal counter at distances nearer to the shower axis.This is because the ADC, Fig. 8: Relative bias and relative standard deviation of the reconstructed µ ML (450) as a function of the logarithm of the primary energy for ADC, ADCProfile, and Profile reconstructions, compared against an ideal counter The plot corresponds to the iron primaries at 30°zenith angle.The error bars are smaller than the marker size for all reconstruction methods.
which dominates in this region, is more suitable for measuring high values of muon density than the binary mode.The ADC and ADCProfile methods yield a smaller standard deviation compared to the Profile method at distances closer to the shower axis.This is because the ADC dominates in the ADCProfile method at distances less than 450 m from the shower axis, whereas the binary mode dominates at larger distances.The ADCProfile reconstruction method has a relative standard deviation comparable to that of the Profile reconstruction method at distances away from the shower axis.
Based on the relative bias and relative standard deviation, the ADCProfile method has the best reconstruction performance.
The coverage probability quantifies the quality of the parameter errors.This is the fraction of events in which the confidence interval includes the true value from the simulated MLDF.In the reconstruction, the 1σ errors of the MLDF normalization µ 0 and the slope parameter β are calculated for each event.The parabolic parameter error is obtained directly from the fit procedure [25].Figure 10 displays the coverage of µ ML (450) for the reconstructions of iron primaries at 30°zenith angle and different energies, using different reconstruction methods.The dotted line at 68% corresponds to the coverage of the 1σ Gaussian confidence interval.As depicted in the figure, the coverage probabilities of all reconstructions are close to each other and close to the Gaussian reference.
To thoroughly examine the impact of saturated events on the reconstruction process, we plot in the top panel of Fig- ADCProfile method Profile method Ideal counter method ADC method Fig. 10: The coverage of the 1σ confidence interval of µ ML (450) as a function of the logarithm of energy for iron primaries at 30°zenith angle.The dotted line at 68% corresponds to the coverage of the 1σ confidence interval of Gaussian likelihood.The ADC, ADCProfile, and Profile reconstruction methods are compared against an ideal counter.ure 11, the reconstructed value of µ ML (450) for iron pri- maries with an energy of 10 18.5 eV and a zenith angle of 30°.The ADC method was used for the reconstruction, and the blue and red histograms indicate the reconstructed candidate and saturated events, respectively.In the bottom panel of the same figure, we have plotted the relative bias and relative standard deviation against energy, after excluding the saturated events for all reconstruction methods.As shown in the figure, the ADC, ADCProfile, and Profile reconstruc-tion methods exhibit better performance and have a relative standard deviation that is closer to that of an ideal counter, particularly at higher energies.Fig. 11: A comparison of the effects of saturated events in the reconstruction procedure.The top panel shows the histograms of the reconstructed µ ML (450) for iron primaries at 10 18.5 eV energy and 30°zenith angle for the ADC reconstruction method.The blue and red histograms corresponds to the reconstructed candidate and saturated events respectively.The bottom panel shows the relative bias and the relative standard deviation of the reconstructed µ ML (450) against the logarithmic energy for iron primaries at 30°z enith angle for different reconstruction methods, excluding the saturated events.
In this study, the saturation level in the LG channel of the ADC was set to 1086 muons.To assess the performance of the ADC reconstruction method at various levels of saturation, 10,000 events were reconstructed at lower (∼ 300 muons) and higher (∼ 2500 muons) saturation levels.Figure 12 compares the relative bias and relative standard deviation of µ ML (450) at different levels of ADC saturation.The plot indicates that although a higher ADC saturation levels result in better relative standard deviation, especially at larger energies, this improvement is minimal and not significant compared to the unrealistic saturation level chosen.It is worth noting that the results achieved for iron nuclei at 45 • , as well for protons at 30 • and 45 • , are similar to those for iron nuclei at 30 .In all cases, the reconstruction outperforms both the ADC and Profile reconstructions across the entire energy range lyzed.The relative standard deviation of µ ML (450) for the ADCProfile reconstruction is the closest to that of ideal counter.

Conclusions
The electronics of segmented detectors can be signed to operate using two modes acquisition: the binary mode and the ADC mode.Currently, both acquisition modes are available for the underground muon detectors at the Pierre Auger Observatory.In the past, muon detectors at AGASA also had both acquisition modes, which could be used independently or in combination to estimate the density of muons hitting a detector.In this study, new method was developed to reconstruct the lateral distribution function based solely on the ADC mode.The underground muon detectors at the Pierre Auger Observatory were used for this case study.It was shown that the performance of the ADC method is similar to the Profile method developed for the binary mode.Specifically, at lower energies, the Profile method outperforms the ADC-based method slightly, whereas at high energies, the ADC-based method performs slightly better than the Profile method.A second method was also developed in this paper based on both acquisition modes, where the Profile likelihood is considered for muon detectors with relatively low values of the muon density, and the ADC likelihood is considered for high values of the muon density.The performance of the combined method is similar or better than the Profile and ADC methods, depending on the primary energy.The detector resolution corresponding to the standard deviation of the estimated number of muons at 450 m from the shower axis obtained with the combined method is quite close to that corresponding to an ideal counter simulated in the range of energies.The performance of the different reconstruction methods was established from the simulations.The combined method can be used to reconstruct the experimental data, and its consistency with the other methods can be verified.saturated events for the Profile reconstruction method.The fraction of saturated events is higher for the Profile method with a 40 ns window and increases rapidly with energy compared to that of the 25 ns window.As shown in equation ( 9), when θ 2 1, the sum of n lognormal terms is expected to exhibit log-normal behavior.Simulations confirmed that Q possesses log-normal characteristics.Figure 15 compares histograms of the sum of n random log-normal variables with corresponding log-normal functions (dotted lines) with parameter values m = 5 and θ = 0.5 [10].
The figure shows that equation ( 9) is an excellent approximation for the LG channel.Therefore, when the variance of the distribution is small (θ 2 1), the sum of n log-normal terms behaves similarly to the sum of narrowly distributed random variables, regardless of the value of n.The figures reveal that, for smaller n values, the distribution has a prominent tail toward positive total charge values, indicating a larger shape parameter.However, for n = 100, the tail of the distribution is almost negligible, and the shape parameter is insignificant.

Fig. 1 :
Fig. 1: Comparison of the detector resolution for the ADC mode, binary mode and an ideal Poisson detector.The variable x in σ [x] refers to the estimator of the number of muons corresponding to each mode.

Fig. 2 :
Fig. 2: Average muon arrival time histogram showing the fraction of muons in each time bin at three different distances to the shower axis for iron showers at a zenith angle of 30°and energy of 10 18.5 eV.

Fig. 3 :
Fig. 3: Simulated unsaturated and saturated pulses at the LG channel of the ADC for an iron shower of 10 19 eV and θ = 30 • .The distances to the shower axis of the unsaturated and saturated pulses are ∼ 876 m and ∼ 211 m, respectively.

Fig. 4 :
Fig. 4: Comparison of different likelihoods used in this work.The plots correspond to iron primary at 10 18 eV at 30°zenith angle.The maximum likelihood estimator of µ is obtained from Profile likelihood and has the value µ = 423

Fig. 12 :
Fig.12:A comparison of the relative bias and relative standard deviation of the reconstructed µ ML (450) with the ADC reconstruction method at three different saturation levels of the LG channel as a function of the logarithm of energy.
of the z score for LG channel.

Fig. 14 :
Fig.14: Fraction of saturated events for primaries at 30°zenith as a function of the logarithm of primary energy for Profile reconstruction method at a 25 ns and at 40 ns inhibition window with error bars.The error bars are smaller than the marker size.

Fig. 15 :
Fig.15: Comparison of the simulated charge distribution with the analytical expression in equation (9) for the LG channel of the UMDs at the Auger Observatory.
Relative bias and relative standard deviation of the reconstructed µ ML (r) as a function of distance to the shower axis for iron primaries at 30°zenith angle and 10 18.25 eV energy The dotted line represents the reference distance of the UMD.The ADC, ADCProfile, and Profile reconstruction methods are compared against an ideal counter.The error bars are smaller than the marker size for all reconstruction methods.