Measurement of the average shape of longitudinal profiles of cosmic-ray air showers at the Pierre Auger Observatory

The profile of the longitudinal development of showers produced by ultra-high energy cosmic rays carries information related to the interaction properties of the primary particles with atmospheric nuclei. In this work, we present the first measurement of the average shower profile in traversed atmospheric depth at the Pierre Auger Observatory. The shapes of profiles are well reproduced by the Gaisser-Hillas parametrization within the range studied, for E &gt; 10^{17.8} eV. A detailed analysis of the systematic uncertainties is performed using 10 years of data and a full detector simulation. The average shape is quantified using two variables related to the width and asymmetry of the profile, and the results are compared with predictions of hadronic interaction models for different primary particles.

section_content: Ultra-High Energy Cosmic Rays (UHECRs) are the most energetic particles known in the Universe. The study of the cascades resulting from their interactions with atmospheric nuclei can provide a unique glimpse into hadronic interaction properties at center-of-mass energies more than one order of magnitude above those attained in human-made colliders. The collision of a UHECR with an atmospheric nucleus initiates an extensive air shower of secondary particles developing in the traversed air mass, usually referred to as slant depth, X. 

In a high-energy hadronic interaction, most of the secondary particles are pions, of which around one third are π 0 mesons. These immediately decay into two photons and initiate an electromagnetic cascade that is dominated by e ± and γ. The charged pions, along with the other secondary hadrons produced in smaller numbers such as kaons and protons/neutrons, constitute the hadronic cascade. Since the interactions of e ± and γ yield virtually no hadrons, an air shower can be considered as the sum of two independent cascades: a hadronic cascade, waning as it penetrates further in the atmosphere losing energy via the decay of neutral pions, and an electromagnetic cascade that is constantly being fed by the hadronic counterpart. After only a few generations the vast majority of the total energy of the primary particle has been transferred to the electromagnetic part of the cascade. 

The cascade progresses until the average energy of single e ± and γ particles falls below the critical energy at which energy is lost predominantly by collisions instead of radiative processes. Atmospheric nitrogen molecules are excited by the passage of charged particles, and the subsequent nitrogen de-excitation results in the emission of a quantity of fluorescence light which is proportional to the energy lost by the shower electrons. The energy deposited by an air shower as a function of the traversed depth is known as its longitudinal profile, which can be measured by the detection of fluorescent light at the ground. 

The integral of the longitudinal profile gives a calorimetric measurement of the shower energy. To measure this energy, we need to estimate the full profile. A functional form must be used to extrapolate outside the observed region. The atmospheric depth at which the profile has a maximum, X max, is the variable most sensitive to the cross-section of the first interaction and to the mass of the primary particle. The precise shape of the energy deposit profile, however, has remained largely untested. This paper describes the first measurement of the shape of the longitudinal profile as a function of traversed atmospheric depth for UHECRs with energies above 10 17. 8 eV. While measured profiles for individual events often have large uncertainties, particularly at values far from the maximum, in this work a high precision is achieved by averaging, in energy bins, showers detected by the fluorescence detectors of the Pierre Auger Observatory. 

The motivation for this measurement is three-fold. The first is to cross-check the assumption that shower profiles are well described by the currently used parametrization. The second is to provide a new way to control the quality of the shower reconstruction for the fluorescence detector. The third is to use the shape parameters to make new independent tests on hadronic interaction models and analysis of primary cosmic-ray composition. 

This paper is organized as follows. The Pierre Auger Observatory and the event reconstruction procedure are described in section 2. In section 3, the functional form used to fit the average profiles is defined, with its two variables and their interpretation. In section 4, we show how the average longitudinal profiles are reconstructed. The analysis is validated by comparing the average profiles obtained after full detector simulation and data reconstruction to the ones calculated directly from the simulated energy deposits in the atmosphere. The systematic uncertainties associated with the measurement are estimated in section 5. Finally, the results of the fit are presented in section 6, and the shape variables measured at each energy are compared to the expectations for proton and iron initiated showers obtained from different hadronic interaction models.
section_name: Introduction

section_content: The Pierre Auger Observatory [1] is a hybrid detector, consisting of a 3000 km 2 Surface Detector (SD) overlooked by the Fluorescence Detector (FD). The SD is composed of 1600 water-Cherenkov detectors separated by 1. 5 km. The FD consists of four sites with six telescopes each. The field of view of each telescope spans 30 • in azimuth and ranges from 1. 5 • to 30 • in elevation. Three additional telescopes called HEAT (High Elevation Auger Telescope) cover the elevation range from 30 • to 60 •. This range is important for showers with energies lower than the ones studied in this paper and thus HEAT data were not used here. 

The measurement of atmospheric properties is essential for the reconstruction of the air showers measured by the detectors mentioned above. The molecular properties (temperature, humidity and pressure height profiles) are provided by the Global Data Assimilation System in three-hour intervals [2]. The aerosol content is monitored hourly by calibrated laser shots from two laser facilities located near the center of the SD array, and cross-checked by LIDAR stations at each FD site. Cloud coverage in the shower path is measured by the LIDAR and cloud cameras at each site (every 15 and every 5 minutes, respectively), and complemented with data from the Geostationary Operational Environmental Satellite (GOES). 

The reconstruction of shower profiles from FD data proceeds in the following steps (see e. g. [1] ). Firstly, the shower-detector plane, spanned by the pointing directions of pixels in the shower image, is calculated. The shower axis within this plane is obtained using the timing information of each pixel, as well as the timing of the closest SD station with signal (hybrid reconstruction). 

For each time bin, a vector pointing from the telescope to the shower is defined, and the signals of all photomultipliers (PMTs) pointing to the same direction within a given opening angle are summed. This angle is determined event-by-event by maximizing the ratio of the signal to the accumulated noise from the night sky background. 

Given the reconstructed geometry, the signal at each time bin can be converted into an energy deposited by the shower at a given slant depth. Every time bin i is projected to a path of length l i along the shower track. The slant depth, X i, is inferred by integrating the atmospheric density through l i. During its path from the shower axis to the FD, light is attenuated due to Rayleigh scattering on air and Mie scattering on aerosols. The light emitted on the shower track at time bin i can be calculated from the measured light at the aperture corrected by this attenuation factor. 

The detected photons correspond to different light emission mechanisms and can reach the telescope directly or by scattering in the atmosphere. Fluorescence light is emitted isotropically along the shower track. High-energy charged particles emit Cherenkov light in a forward-concentrated beam. Even if the shower does not point directly to the detector, a fraction of this beam will be scattered into the field of view. This fraction is calculated taking into account the characteristics of both molecular and aerosol scattering in the atmosphere. 

The Cherenkov and fluorescence light produced by an air shower are connected to the energy deposit by a set of linear equations [3]. The profile of energy deposit as a function of slant depth is functionally described by the Gaisser-Hillas parametrization [4],

which has four parameters: the maximum energy deposit, (dE/dX) max, the depth at which this maximum is reached, X max, and shape parameters X 0 and λ. This function is used in the calculation of the Cherenkov beam accumulated up to X i, which determines the number of Cherenkov photons seen at the aperture. The proportionality between the number of fluorescence photons and the energy deposit is given by the fluorescence yield [5], which depends on the molecular properties of the atmosphere. The statistical uncertainty in dE/dX i is calculated from the Poisson uncertainty of photoelectrons detected by the photomultipliers of the fluorescence telescopes. 

The four parameters that describe the shower profile and their uncertainties are obtained from a log-likelihood fit to the number of photoelectrons detected at the PMTs. The calorimetric energy is the integral of f GH, and the total energy of the shower is estimated by correcting for the invisible energy carried away by neutrinos and high-energy muons.
section_name: Event reconstruction at the Pierre Auger Observatory

section_content: The maximum energy deposit of the longitudinal profile, (dE/dX) max, is proportional to the energy of the primary particle and varies by three orders of magnitude in the energy range studied in this work. The X max value of each shower is, on average, a characteristic of the primary particle mass, but it varies greatly also for showers with the same primary, due to the depth at which the first interaction occurs, the large phase space for high-energy interactions and the general stochastic nature of shower development. 

The focus of this work is to separate these two parameters, (dE/dX) max and X max, to isolate the information contained in the profile shape itself. First, each measured shower profile is normalized to its fitted (dE/dX) max. With this rescaling, all showers have a maximum value at 1. Then, the showers are transformed to the same development stage by translating the slant depth by X max, i. e., X ≡ X -X max, thus all profiles are centered at zero. The normalized energy deposit profile, with (dE/dX) ≡ (dE/dX)/ (dE/dX) max, can be described by the Gaisser-Hillas function, written as a function of parameters R and L [6],

where R = λ/|X 0 |, L = |X 0 |λ and X 0 ≡ X 0 -X max. In this notation, the Gaisser-Hillas function is a Gaussian with standard deviation L, multiplied by a term that distorts it, with the asymmetry governed by R, and thus these parameters are less correlated than λ and X 0. 

An equivalent parametrization as a function of the Full Width at Half Maximum, f FWHM, and asymmetry, f, has been reached independently [7]. 

Note that in previous analyses by the HiRes/MIA [8] and HiRes Collaborations [9], the energy deposit profile was constructed as a function of shower age (s = 3X/(X + 2X max )), and the resulting shapes were found to be compatible with a Gaussian distribution with standard deviation σ age. This width, however, is convolved with (and dominated by) the X max value 1. In this work we chose to represent the profiles in atmospheric depth because it preserves the measured event-by-event shape. Also, R and L taken from the fit to the profiles in atmospheric depth have been shown to be sensitive to the mass of the primary cosmic rays and to the properties of the high-energy interactions that occur at the start of the shower. The sensitivity is kept by an average profile [10], i. e., the average shape of a sample of profiles initiated by different primary particles is related to the average mass composition in the sample.
section_name: Average longitudinal shower profile

section_content: The event selection used here is based on the selection criteria developed by the Pierre Auger Collaboration to measure the X max distributions [11], and is applied to data covering the period from January 2004 to March 2015. Good weather conditions are required, with no clouds in the sky and a measured vertical aerosol optical depth lower than 0. 1 at a height of 3 km. 

A good shower profile reconstruction requires an accurate knowledge of the shower geometry, which is obtained by adding information from SD single-station triggers in coincidence with the FD shower trigger. An hybrid geometry reconstruction is performed using the start time of the pulse from the SD station with the largest signal in conjunction with the times from the triggered FD pixels. The SD station used should be within 1500 m of the core, and have an expected trigger probability for the specific geometry above 95% for both proton and iron showers. For showers landing within the SD array, the closest station distance is usually within 750 m of the core, and the expected one-SD station trigger probability is 100% for energies above 10 17. 6 eV. So, the hybrid geometry requirement does not introduce any selection bias. 

On the measured longitudinal profile, strict cuts are made: at least 300 g/cm 2 must be observed, including the X max depth, for which the expected resolution must be less than 1 Shower age is defined as s(X ) = 1 + 2X 3Xmax-X. As L is the profile width in depth, a relation between σage and L can be found by making X = ±L in the previous formula. Doing the derivative with respect to L and Xmax, we find that

L∆(Xmax) ≈ 1/6, where ∆ stands for the proton-iron difference in each variable. So the majority of the composition separation in σage comes from Xmax. Figure 2. Average longitudinal shower profile for simulated events with energies between 10 18. 8 and 10 19. 2 eV. The profiles for proton and iron are shown in blue and red, respectively. In the bottom plot, the ratio of reconstructed profiles to the generated ones is shown. 

40 g/cm 2. The minimum angle between the shower axis and the pointing vector of any of the pixels with signal has to be larger than 20 •, to minimize direct Cherenkov light that induces larger geometrical uncertainties in the determination of X max. Finally, a fiducial field of view is defined to guarantee a uniform acceptance over a range in X max that covers the vast majority of the observed distributions in data, to reduce possible composition bias from specific geometries. In this work, two additional selection cuts with respect to the criteria developed in [11] were applied. When showers cross two telescopes, differences in alignment cause time residuals in the geometry fit, as the distance to the shower estimated by both telescopes is different. The estimation of the atmospheric attenuation depends on the distance and, therefore, also the energy deposit is affected. Large residuals were found in one of the telescopes (Coihueco 6), so events in which the profiles were measured using this telescope are excluded. Also, events for which the pixel time-fit used to determine the shower axis yielded very large reduced χ 2 values (above 5) are excluded. The first cut affects approximately 3% of the events (649), while in the second only two events are excluded. 

In total, 21532 events are selected for analysis. These are divided into six energy bins. The shower profiles are constructed in 10 g/cm 2 bins in X, in which each energy deposit is accumulated with a weight corresponding to the inverse of its squared error. The average profiles for all energies are shown in figure 1. 

The reconstruction and analysis of the average longitudinal shower profiles was validated with a full detector simulation for energies between 10 17. 8 and 10 20 eV with proton and iron as primary particles. Comparing simulated and reconstructed showers an excellent agreement is found for the central part, but reconstruction increasingly deviates for negative X, and a lower bound should be set for the analysis. At depths beyond X max, showers lose most of the primary information, but the later parts of the profile are necessary to measure the profile width. An upper bound should maximize the relative importance of the profile start while keeping the statistical uncertainty on L and R below the proton/iron separation. While no detailed optimization was attempted for this first measurement, common bounds for all energies were set at [-300,+200] g/cm 2, fulfilling the above requirements. 

Figure 2 shows the generated and reconstructed profiles for proton and iron simulations for energies around 10 19 eV; figure 3 shows the data reconstruction for the same energy bin, together with the different light components along the profile. Note that the chosen fit region is dominated by fluorescence light, directly proportional to the shower energy deposit along X. 

The shower profiles are fitted with equation 3. 1, leaving all parameters unconstrained. In addition to R and L, the normalization and the maximum position are allowed to vary around (dE/dX) = 1 and X = 0 g/cm 2 to account for possible effects of smearing due to X max resolution and corresponding bias in the average energy deposit, and for possible systematic effects in the extraction of X max and (dE/dX) max from the fit to single event profiles. The results for data are within the predictions from full-detector simulations of showers, i. e., at all energies, the normalization is unitary within 0. 5%, and the maximum position is within ±1 g/cm 2. 

To minimize the residual reconstruction biases, the differences between the reconstructed and generated profile parameters are calculated separately in proton and iron simulations, and the average across both primaries are subtracted from the values of L and R obtained in data. This bias correction is compatible with zero for energies above 1 EeV and is lower than 2 g/cm 2 in L and 0. 005 in R (with different signs for each primary), for the lowest energy bin analysed. Half of the applied correction is included as a systematic uncertainty.
section_name: Data selection and Monte Carlo validation

section_content: The reconstruction of longitudinal profiles from the measurements at the FD requires several steps with which systematic uncertainties are associated. This section describes the effects we have considered, which are summarized in table 1. 

Atmospheric conditions play a crucial role in the propagation of the light from the shower to the fluorescence detector [12]. Several systematic uncertainties related to this process were studied. They included the impact of possible patches of clouds in the sky, differences found when separating data by the seasons of the year, and uncertainties in the overall aerosol content and its height dependence. The atmospheric aerosol attenuation, τ A (h), is measured hourly with the central laser facility [13, 14]. The measurement compares the number of photons detected on the FD (as a function of height) in a given laser run with the one detected on a clear reference night. The aerosol height profile on the shower path can be calculated from τ A (h), but has two main sources of uncertainty: the determination of the aerosol profile on the reference night (which is fixed for all showers) and the propagation of the uncertainty on the laser measurement at the FD to the path to each air shower (which varies from bin to bin and depends on geometry and atmospheric conditions). The uncertainty related to the estimation of the aerosol content was found to be the largest in this work, yielding approximately a ±0. 02 and ±5 g/cm 2 uncertainty in R and L, respectively. These values are obtained by recalculating the average profiles after again reconstructing the full data sample, first with τ A (h) changed coherently for all events according to the reference night uncertainties, and then with τ A (h) varied within the non-correlated uncertainties of each event. 

The uncertainties in the determination of the different light components were also considered. First, the reconstruction was repeated changing the fluorescence and Cherenkov yield values within their uncertainties (4% and 5%, respectively), and accounting or not for the multiple scattering corrections. Then, the data was separated according to the fluorescence fraction of the event, which gave a larger systematic difference. When only showers with fluorescence fractions lower than 90% (the average value for the analyzed sample) are used, the resulting average profile is 2 g/cm 2 larger in width, L, while the change in R is negligible. 

The fit of the longitudinal profile of individual events is constrained in the integral of the normalized longitudinal profile, by an energy dependent value taken from shower simulations. Furthermore, the Gaisser-Hillas parameters are constrained by X 0 = -121 ± 172 g/cm 2 and λ = 61 ± 13 g/cm 2, values found previously in a small ensemble of very high-energy events with long tracks in the FD for which the unconstrained fit was possible [11]. The uncertainty in the estimation of the profile constraints was propagated in the reconstruction by shifting the central values by their standard deviation. This contribution is relatively more important in L (around 3 g/cm 2 ) than in R (0. 01). 

The reconstructed shower geometry is given first by the zenith and azimuth angles of the shower-detector plane, and then the shower axis is defined by distance, angle and time references. All these five parameters are varied within ±1σ to obtain the systematic uncertainty associated with the reconstructed shower geometry. To ensure the quality of the reconstruction, the dependence of R and L on the zenith angle or distance from X max to the telescope was also studied, but was found to be small and contained in the geometric systematic uncertainty of around 2 g/cm 2 in L and 0. 02 in R. 

The effect of the uncertainty in telescope alignment was also considered. One of the telescopes had been previously excluded from the analysis by looking directly at residuals in the timing at the crossing between telescopes. Smaller effects were tested by studying the telescope-to-telescope difference of the reconstructed shape. These are smaller than other systematic uncertainties, but at the level of 0. 012 in R and 1. 8 g/cm 2 in L. 

Finally, the uncertainty of the energy scale of 14% [15] and the previously described proton-iron discrepancy from the bias correction are also added, but they are small (below 1 g/cm 2 in L and 0. 01 in R) in comparison with the previously described ones. 

The uncertainties vary with energy and are asymmetric. Table 1 shows the largest value for each category, and table 2 shows the total error per energy bin.
section_name: Systematic uncertainties

section_content: The fit of data profiles with the Gaisser-Hillas parametrization is shown in figure 4. The fitted function follows the data points through the depth range used in this work, [-300, +200] g/cm 2, with residuals within the statistical uncertainties. This is the first experimental demonstration that the Gaisser-Hillas parametrization is an accurate functional form to describe the central part of UHECR longitudinal profiles above 10 17. 8 eV at 1% accuracy. 

The values of R and L obtained from the fit for the six energy intervals studied are shown in table 2, along with their energy dependent statistical and systematic uncertainties. The average energy and number of events, N, in each bin are also listed for reference. 

The existing hadronic interaction models give different predictions for the shape variables. CORSIKA [16] was used to simulate proton, helium, nitrogen and iron showers with the EPOS-LHC [17], QGSJetII-04 [18] and Sibyll2. 3c [19] models. The evolution of R and L with energy, along with their respective systematic and statistical uncertainties, is shown in figure 5. Both the asymmetry, R, and the width, L, in data agree well with the predicted values for all models. For the asymmetry, all models give similar predictions, and the results seem to point to the composition becoming heavier with energy, although current system- atic uncertainties still hinder any composition claim. For L, data is consistent with a linear increase with log 10 (E/eV). L is compatible with the predictions of Sibyll2. 3c for all compositions, but points to a lighter composition if compared with the other two models, that predict smaller values of L. 

To understand better the interplay of the two measured variables, it is interesting to see the results in the (R, L) plane for two fixed energies (figure 6 ). In these plots, all composition scenarios are represented (as a combination of proton, He, N and Fe) for a given energy. For all models, proton has a lower R and larger L than iron, so moving from the top left to the bottom right points within a given model, makes the transition from lighter to heavier primaries. It is interesting, however, to note that the areas have different shapes: some models predict, at a given energy, a smaller L for proton than for Helium, while others do the opposite. 

Since R and L are experimentally correlated, these plots provide a smaller phase space Figure 6. L vs R for the energy bin 10 18 to 10 18. 2 eV (left) and from 10 18. 8 to 10 19. 2 eV (right). The inner dark grey ellipse shows the fitted value for data and its statistical uncertainty, and the outer light grey line the systematic uncertainty. For each hadronic model proton, helium, nitrogen and iron showers were simulated and average profiles were built making all possible combinations. Each of the points represents the value of R and L for a given model and composition combination, so the phase space spanned by each model is contained in its respective colored area. Pure proton is, for each model, on the upper left side (low R and high L) and the transition to iron goes gradually to the lower right side. 

within which to constrain the predictions of the hadronic interaction models. In figure 6, for the 10 18 eV energy bin (left), it can be seen that the average value in data is in the area occupied by most models for a light composition, while at 10 19 eV (right) it is within the predictions for intermediate mass primaries. Within the current experimental resolution, the data are fully compatible with most composition scenarios at the 2σ level for all models. It is, however, interesting to note that Sibyll2. 3c occupies a different phase space than the other two post-LHC models, indicating that a higher precision measurement of the profile shape could provide a test on the hadronic interaction models.
section_name: Results

section_content: In this work, the Pierre Auger Observatory has been used to measure the average shape of the longitudinal profile of air showers. The method was first validated in a full detector simulation of proton and iron primaries, which showed that reconstructed and simulated profiles are in very good agreement for all energies above 10 17. 8 eV. The average longitudinal profiles as a function of atmospheric depth of ultra-high energy cosmic rays have been presented for the first time. They are well described by a Gaisser-Hillas function throughout the fitting range chosen, -300 to +200 g/cm 2 around the shower maximum, to a 1% level precision. 

The systematic uncertainties contributing to the measurement were estimated, and we found that the uncertainty in the atmospheric aerosol content is the main factor affecting the shape of the reconstructed longitudinal profile. In fact, the variations of L and R values obtained when dividing the sample according to different detection conditions are all within the range allowed by the atmospheric uncertainty. The validation of the profile shape around the shower maximum increases our confidence in its use for the determination of the shower calorimetric energy and of the X max parameter. 

The two shape parameters measured in this work were compared with model predictions and found to be compatible with most mass composition scenarios at the two sigma level. However, different hadronic models predict shape parameters which span different regions in the (R, L) plane, even when considering all possible mass compositions. Therefore, a future measurement of the average longitudinal profile shape with smaller systematic uncertainties could provide constraints on high-energy interaction models.
section_name: Summary

section_name: Acknowledgments

section_content: The successful installation, commissioning, and operation of the Pierre Auger Observatory would not have been possible without the strong commitment and effort from the technical and administrative staff in Malargüe. We are very grateful to the following agencies and organizations for financial support:

Argentina -Comisión Nacional de Energía Atómica ; Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT) ; Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET) ; Gobierno de la Provincia de Mendoza ; Municipalidad de Malargüe ; NDM Holdings and Valle Las Leñas ; in gratitude for their continuing cooperation over land access; Australia -the Australian Research Council; Brazil -Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) ; Financiadora de Estudos e Projetos (FINEP) ; Fundação de Amparo à Pesquisa do Estado de Rio de Janeiro (FAPERJ) ; São Paulo Research Foundation (FAPESP) Grants No. 2010/07359-6 and No. 1999/05404-3 ; Ministério da Ciência, Tecnologia, Inovações e Comunicações (MCTIC) ; Czech Republic -Grant No. MSMT CR LTT18004, LO1305, LM2015038 and CZ. 02. 1. 01/0. 0/0. 0/16_013/0001402 ; France -Centre de Calcul IN2P3/CNRS ; Centre National de la Recherche Scientifique (CNRS) ; Conseil Régional Ile-de-France ; Département Physique Nucléaire et Corpusculaire ( PNC-IN2P3/CNRS ); Département Sciences de l'Univers (SDU-INSU/CNRS) ; Institut Lagrange de Paris (ILP) Grant No. LABEX ANR-10-LABX-63 within the Investissements d'Avenir Programme Grant No. ANR-11-IDEX-0004-02 ; Germany -Bundesministerium für Bildung und Forschung (BMBF) ; Deutsche Forschungsgemeinschaft (DFG) ; Finanzministerium Baden-Württemberg ; Helmholtz Alliance for Astroparticle Physics (HAP) ; Helmholtz-Gemeinschaft Deutscher Forschungszentren (HGF) ; Ministerium für Innovation, Wissenschaft und Forschung des Landes Nordrhein-Westfalen ; Ministerium für Wissenschaft, Forschung und Kunst des Landes Baden-Württemberg; Italy -Istituto Nazionale di Fisica Nucleare (INFN) ; Istituto Nazionale di Astrofisica (INAF) ; Ministero dell'Istruzione, dell'Universitá e della Ricerca (MIUR) ; CETEMPS Center of Excellence; Ministero degli Affari Esteri (MAE) ; México -Consejo Nacional de Ciencia y Tecnología (CONACYT) No. 167733 ; Universidad Nacional Autónoma de México (UNAM) ; PAPIIT DGAPA-UNAM; The Netherlands -Ministry of Education, Culture and Science; Netherlands Organisation for Scientific Research (NWO) ; Dutch national e-infrastructure with the support of SURF Cooperative; Poland -National Centre for Research and Development, Grant No. ERA-NET-ASPERA/02/11 ; National Science Centre, Grants No. 2013/08/M/ST9/00322, No. 2016/23/B/ST9/01635 and No. HARMONIA 5-2013/10/M/ST9/00062, UMO-2016/22/M/ST9/00198 ; Portugal -Portuguese national funds and FEDER funds within Programa Operacional Factores de Competitividade through Fundação para a Ciência e a Tecnologia (COMPETE) ; Romania -Romanian Ministry of Research and Innovation CNCS/CCCDI-UESFISCDI, projects PN-III-P1-1. 2-PCCDI-2017-0839/19PCCDI/2018, PN-III-P2-2. 1-PED-2016-1922, PN-III-P2-2. 1-PED-2016-1659and PN18090102 within PNCDI III; Slovenia - Slovenian Research Agency; Spain -Comunidad de Madrid ; Fondo Europeo de Desarrollo Regional (FEDER) funds; Ministerio de Economía y Competitividad ; Xunta de Galicia ; European Community 7th Framework Program Grant No. FP7-PEOPLE-2012-IEF-328826 ; USA -Department of Energy, Contracts No. DE-AC02-07CH11359, No. DE-FR02-04ER41300, No. DE-FG02-99ER41107 and No. DE-SC0011689 ; National Science Foundation, Grant No. 0450696 ; The Grainger Foundation ; Marie Curie-IRSES/EPLANET ; European Particle Physics Latin American Network ; European Union 7th Framework Program, Grant No. PIRSES-2009-GA-246806 ; and UNESCO.
section_name: Acknowledgments