ARTICLE IN PRESS Nuclear Instruments and Methods in Physics Research A 588 (2008) 171–175 www.elsevier.com/locate/nima Multiparametric topological analysis (MTA) for the study of the primary CR composition: Performances with Auger simulated data D. D’Ursoa,b,, M. Ambrosioa, C. Aramoa, F. Guarinoa,b, L. Valorea,b, for the Pierre Auger Collaborationc a Sezione INFN di Napoli, Italy Università di Napoli ‘‘Federico II’’, Italy c Pierre Auger Observatory: Av. San Martin Norte 304, (5613) Malargüe, Argentina b Available online 12 January 2008 Abstract We describe the application of a multiparametric analysis to estimate the UHE Cosmic Rays composition. The proposed method, MTA (Multiparametric Topological Analysis), is based on the study of the correlations among different shower observables. This technique is designed to fully exploit the complementarity of Auger fluorescence and ground array data. In the present work, we report the results of the application to Conex showers, fully simulated through the Auger detector, using only parameters describing the longitudinal development of air showers as recorded by fluorescence detector for hybrid data. r 2008 Elsevier B.V. All rights reserved. Keywords: Ultra high energy cosmic rays; Auger observatory; Mass compostion; Multiparametric analysis 1. Introduction The understanding of the nature of the ultra high energy cosmic rays (UHECR) is a crucial point towards the determination of their origin, acceleration and propagation mechanism. In the energy range of the knee (about 1015 eV) is widely accepted that cosmic rays have a galactic origin and their most probable source are the supernova remnants (SNR). Above 1018 eV all observed cosmic rays are presumed to come from extragalactic sources, because there are no known galactic sources able to produce particles up to such energies and they cannot be contained in our galaxy long enough to be accelerated. The energy at which the transition between galactic and extragalactic cosmic rays occurs is still unknown and only a detailed knowledge of the particle spectrum and composition could allow to distinguish among different models (see Ref. [1] for further details). Designed to solve the UHECR puzzle, the Pierre Auger Observatory [2] collects data with unprecedented precision Corresponding author at: Sezione INFN di Napoli, Italy. E-mail address: domenico.durso@na.infn.it (D. D’Urso). 0168-9002/$ - see front matter r 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.nima.2008.01.034 and statistics. The Observatory is a ‘‘hybrid detector’’ composed by an array of 1600 water Cherenkov particle detectors (tanks), distributed over an area of 3000 km2 (Surface Detector, SD) overlooked by 24 fluorescence telescopes (Fluorescence Detector, FD). The combined use of these two detection techniques provides cross-calibration and a better event reconstruction accuracy. Moreover, the combined use of FD and SD observables will help to put tighter constraints on hadronic interaction models. A description of Auger hybrid performances is given in Ref. [2]. Nowdays, the Auger Observatory has recorded more events above 1018 eV than all the previous cosmic ray experiments on UHECR. The most popular method developed so far to infer the primary composition from fluorescence experiments makes use of the depth of the maximum of the cascade development ðX Max Þ and derives the observed mean mass composition as a function of the primary energy comparing the measured shower maxima with the Monte Carlo predictions (Elongation Rate analysis). Within the Auger Collaboration such studies are performed on hybrid events, in which at least one SD station is involved, with a consequent improvement of the reconstruction ARTICLE IN PRESS 172 D. D’Urso et al. / Nuclear Instruments and Methods in Physics Research A 588 (2008) 171–175 accuracy with respect to pure FD data (shower maximum can be measured with an accuracy better than 40 g cm2 [3]). In surface array experiments, informations on primary composition are usually extracted from signal rise time, time curvature of the shower front, muon content and azimuthal signal asymmetry. This kind of analysis is in progress on Auger SD data alone [4]. However the sensitivity of these methods is limited, due to the strong fluctuations of the physical processes which smear out the informations and produce an overlapping of the observed features of showers initiated by different primaries. Such unsatisfactory situation can improve if more parameters are taken into account, exploiting the unique future of Auger hybrid data, that contain informations on the longitudinal profile and lateral distribution of a shower at the same time. This type of arguments triggered our efforts to find methods which make use of a larger amount of informations and which are capable to allow the identification (at least in terms of probabilities) of the cascade origin also for individual showers. A comparison of the performances of such methods applied to the study of the longitudinal development, recently published by our group [5], indicates that multiparametric topological analysis (MTA) analysis performances are very promising and nearly comparable with those of neural networks. Because of the simplicity of the basic algorithms we decided to use MTA for this study. 2. The MTA method Our aim is to infer the UHECR composition using the full information available in the hybrid Auger data. Any shower observable cannot be unambiguously connected to the primary mass: distributions are always mixed also at fixed energy and we can never identify the individual primary mass without uncertainty. However, using a multiparametric approach, comparing the correlation of a set of discriminant variables extracted from real events with the corresponding quantities derived from simulated ones of known energy and composition, it is possible to define a probability for each individual shower to originate from a nucleus of a given mass. Once a set of measured shower parameters sensitive to primary composition is chosen, MTA method [5] relies on a topological analysis of their correlations. Therefore, application to Auger data requires a full and detailed simulation of showers in both SD and FD detectors. CORSIKA [6] showers are ideal for this purpose, but too heavy both in terms of CPU and storage requirements, to reproduce real data with reasonable statistics. Two solutions are possible: the use of CORSIKA with non optimal thinning, or the use of parametric or semianalytical approaches to shower development. For the study of the longitudinal profile of atmospheric cascades is now available a fast Monte Carlo program, Conex [7], based on a semianalytical approach. The distribution of particles at ground is not provided, nevertheless, tank trigger simulation is performed using a parametrized ‘‘Lateral Trigger Probability’’ function [8] and the time of the station with the highest signal is calculated (a detailed description of the hybrid detector simulation program is given in Ref. [9]). This information is sufficient to operate a hybrid geometry reconstruction [2]. A more complete fast simulator, SENECA [10], which also provides reliable informations of particles at ground, will be soon available. For the time being, we restrict the analysis to the use of only FD parameters for hybrid events by means of Conex simulated showers. Simulated data have been produced by the Conex MonteCarlo, version 2r1.4b, fully processed with the Auger simulation-reconstruction software, using: (i) QGSJET-II-3 [11] as high energy hadronic interaction model; (ii) energy spectrum from 1017:5 eV to 1020:2 eV, according to dN=dE / E 3 , for a total of 2  106 events; (iii) zenith distribution from 0 to 60 , according to dN=ðd cos yÞ / cos y; (iv) uniform azimuth distribution. Quality cuts to simulated data, as described in Ref. [3], have been applied. 3. Definition of the parameter space Once a set of observables is chosen, it defines a parameter space, which is divided into cells whose dimensions are related to the experimental accuracy. The set of shower reconstructed parameters provides the association of the event to one of the cells of the space. For each primary mass and for each energy bin, a set of simulated events, fully processed by the detector simulation and reconstruction codes, is used to populate the parameter space. Restricting to the case of shower quantities measured by the FD, the space is defined by a choice of the four parameters describing the fit with a Gaisser–Hillas function to the shower profile:  ðX Max X 0 Þ=l dE dE X  X0 ¼ eðX Max X 0 Þ=l (1) dX dX Max X Max  X 0 where dE=dX and dE=dX Max are the energy deposit at the depth X and at the shower maximum, respectively. As an example, in Fig. 1 a parameter space built with two parameters from the fit, X Max and l is shown. The parameter space is divided in cells with dimensions 20 and 5 g cm2 , respectively, and is populated with Conex simulated showers induced by the two limit masses, proton (triangles) and iron (circles). ARTICLE IN PRESS D. D’Urso et al. / Nuclear Instruments and Methods in Physics Research A 588 (2008) 171–175 In the near future, the composition analysis will be extended to observables measured by the ground array detector. Considering only two primary masses, protons and iron nuclei, in each cell, which in the most general n-dimensional case is defined by ðh1 . . . hn Þ, one can define ðh1 ...hn Þ the total number of showers N tot , the total number of showers induced by protons and iron nuclei, N Pðh1 ...hn Þ and ðh1 ...hn Þ N Fe , respectively, and then derive the associated frequencies: ðh1 ...hn Þ pPðh1 ...hn Þ ¼ N Pðh1 ...hn Þ =N tot ðh1 ...hn Þ ðh1 ...hn Þ ðh1 ...hn Þ pFe ¼ N Fe =N tot (2) which can be interpreted as the probability for a real shower falling into the cell ðh1 . . . hn Þ to be initiated by protons or iron primary nuclei. QGSJET-II E = 1019.1 eV Proton 173 4. Average mixing probabilities In the case of a data set of M showers, the sample may be seen as composed by a mixture of M  pP proton showers and M  pFe iron induced showers, where pj ¼ M X ðh1 ...hn Þm pj =M (3) m¼1 and the sum is performed over all the experimental data, with (h1 . . . hn Þm indicating the cell interested in by the mth event. A second set of showers is used to compute the mixing probabilities Pi!j that an event of mass i is identified as primary j. The mean Pi!j are obtained by computing pj for samples of pure primary composition i. In Fig. 2 average mixing probabilities Pi!j computed for the 2D example (Fig. 1) for proton and iron are shown. The mean probabilities that the primary mass is correctly identified are of the order of 80%. 5. Determination of the primary mass composition Iron 1200 The mixing probabilities Pi!j can be used for the reconstruction of the primary mass composition as the coefficients in the system of linear equations: 1100 XMax [g/cm2] 1000 900 N 0P ¼ N P  PP!P þ N Fe  PFe!P 800 N 0Fe ¼ N P  PP!Fe þ N Fe  PFe!Fe 700 600 500 400 20 30 40 50 60 70 80 90 100 λ [g/cm2] Fig. 1. X Max vs. l parameter space populated by Conex proton (triangles) and iron (circles) induced showers. where N P and N Fe are the true values, which are altered to N 0P and N 0Fe , due to misclassification. Observed composition is also altered by different triggerreconstruction–selection efficiencies for different primaries, P and Fe for proton and iron. This effect is compensated by applying a correction given by the efficiency ratio to the composition values obtained solving the system. In Fig. 3 the ratio of the overall efficiencies for protons P and irons Fe is shown as a function of the energy obtained applying the analysis cuts described in Ref. [3]. Defined Di ¼ N i =N and if DP and DFe are the solutions, the corrected Proton Iron 1 1 0.8 0.8 Probability Probability (4) 0.6 0.4 0.2 0.6 0.4 0.2 0 0 Proton Iron Proton Iron Fig. 2. Average mixing probabilities estimated for the data in Fig. 1: the mean probability that a proton is correctly identified as proton or that is misidentified as iron (left pad); that one that an iron nucleus is misidentified as proton or correctly classified as iron (right pad). ARTICLE IN PRESS D. D’Urso et al. / Nuclear Instruments and Methods in Physics Research A 588 (2008) 171–175 174 Proton - Iron Efficiency Ratio Reconstructed Composition 1.6 1 1.4 0.8 Composition ∈P /∈Fe 1.2 1 0.8 0.6 0.4 0.6 0.2 0.4 17.5 18 18.5 19 log10 Energy [eV] 19.5 0 20 17.5 Fig. 3. Ratio of the overall efficiencies for protons and irons as a function of the energy. abundancies DP and DFe are given by DP D  P =Fe ¼ P  DFe DFe 18.5 19 log10 Energy [eV] 19.5 20 Fig. 4. MTA mass composition (full markers protons, empty markers irons) as a function of the energy for two proton–iron mixing 90–10% (circles), 70–30% (triangles) compared with the expected values (solid and dashed lines). (5) Elongation Rate Comparison with the conditions 850 DP þ DFe ¼ 1 DP þ DFe ¼ 1. The technique has been tested on simulated samples of known mass composition, in energy bins of 0:2 log10 ðE½eVÞ. The MTA analysis has been applied to 3parameters X Max  l  X 0 . The parameter space has been divided in cells whose dimensions are 20, 5 and 50 g cm2 for X Max ; l and X 0 , respectively. In Fig. 4 MTA composition results (full markers protons, empty markers irons) are shown as a function of the energy for two proton–iron mixing 90–10% (circles), 70–30% (triangles) compared with the expected values (solid and dashed lines). MTA sensitivity to mass composition has been compared with the usual elongation rate analysis performed on the same data: for different simulated composition mixture, the elongation rate computed in the usual way has been confronted with an MTA-elongation rate obtained weighting the mean X Max of pure data samples by the MTA abundancies. Let hX Max ðE i ÞiP and hX Max ðE i ÞiFe be the mean values of the X Max distributions of pure proton and iron samples, DP ðE i Þ and DFe ðE i Þ the MTA percentage for the two masses in the ith energy bin, then it is possible to compute the mean X Max at the energy E i (6) XMax [g/cm2] 800 6. Method performances hX Max ðE i Þi ¼ hX Max ðE i ÞiP  DP ðE i Þ þ hX Max ðE i ÞiFe  DFe ðE i Þ. 18 750 700 Expected ER MTA ER 90%P MTA ER 80%P MTA ER 70%P MTA ER 60%P 650 600 17.5 18 18.5 19 log10 Energy [eV] 19.5 20 Fig. 5. Comparison between the standard elongation rate and that one obtained from MTA results for samples with different proton percentage: 90% (full squares), 80% (empty squares), 70% (full circles) and 60% (empty circles). In Fig. 5 is shown the elongation rate obtained for samples with 90% (full squares), 80% (empty squares), 70% (full circles) and 60% (empty circles) of protons, compared with the expected values (thin lines). Elongation rate curves reproduce the expected ones with high accuracy over the whole energy range. Curves corresponding to proton compositions different at level of 10% are easily distinguishable. The test has shown that, once an hadronic model is chosen, in this case QGSJET-II, MTA applied to FD only ARTICLE IN PRESS D. D’Urso et al. / Nuclear Instruments and Methods in Physics Research A 588 (2008) 171–175 parameters has a sensitivity to the input primary mass better than 5%, in the whole energy range and at each tested proton–iron mixing. The use of both FD and SD observables in the analysis will give better results. It is clear that the dominant uncertainty is related to the present indetermination of hadronic models. Hopefully, Auger data will help to constrain the models. References [1] T. Stanev, astro-ph/0611633; D. Allard, et al., Astropart. Phys. 27 (2007) 61. [2] B. Dawson, Pierre Auger Collaboration, in: Proceedings of the 30th ICRC, Mérida, 2007. [3] M. Unger, Pierre Auger Collaboration, in: Proceedings of the 30th ICRC, Mérida, 2007. 175 [4] M.D. Healy, Pierre Auger Collaboration, in: Proceedings of the 30th ICRC, Mérida, 2007. [5] M. Ambrosio, et al., Astropart. Phys. 24 (2005) 355; M. Ambrosio, et al., in: Proceedings of the CRIS 2004, Nucl. Phys. B (Proc. Suppl.) 136 (2004), astro-ph/0410179; M. Ambrosio, et al., in: Proceedings of Vulcano Workshop 2004, astro-ph/0410184. [6] D. Heck, J. Knapp, J.N. Capdevielle, G. Schatz, T. Thouw, Rep. FZKA 6019 (1998) Forshungszentrum Karlsruhe. [7] T. Pierog, et al., Nucl. Phys. B (Proc. Suppl.) 151 (2006), in: Proceedings of the 13th International Symposium on Very HighEnergy Cosmic Ray Interactions at the NESTOR Institute, Pylos, Greece, 6–12 September 2004, astro-ph/0411260. [8] E. Parizot, Pierre Auger Collaboration, in: Proceedings of the 29th ICRC, Pune, vol. 7, 2005, pp. 71–74. [9] L. Prado, et al., Nucl. Instr. and Meth. 545 (2005) 632. [10] H.J. Drescher, G. Farrar, Phys. Rev. D 67 (2003). [11] S. Ostapchenko, Nucl. Phys. B (Proc. Suppl.) 151 (2006) 143.