Reconstruction of inclined air showers detected with the Pierre Auger Observatory
Abstract
We describe the method devised to reconstruct inclined cosmicray air showers with zenith angles greater than detected with the surface array of the Pierre Auger Observatory. The measured signals at the ground level are fitted to muon density distributions predicted with atmospheric cascade models to obtain the relative shower size as an overall normalization parameter. The method is evaluated using simulated showers to test its performance. The energy of the cosmic rays is calibrated using a subsample of events reconstructed with both the fluorescence and surface array techniques. The reconstruction method described here provides the basis of complementary analyses including an independent measurement of the energy spectrum of ultrahigh energy cosmic rays using very inclined events collected by the Pierre Auger Observatory.
42]A. Aab 65]P. Abreu 54]M. Aglietta 95]M. Ahlers 83]E.J. Ahn 29]I. Al Samarai 17]I.F.M. Albuquerque 1]I. Allekotte 87]J. Allen 89]P. Allison 11, 8]A. Almela 58]J. Alvarez Castillo 76]J. AlvarezMuñiz 41]R. Alves Batista 45]M. Ambrosio 59]A. Aminaei 96, 0]L. Anchordoqui 65]S. Andringa 45]C. Aramo 73]F. Arqueros 1]H. Asorey 65]P. Assis 31]J. Aublin 76]M. Ave 32]M. Avenier 10]G. Avila 69]A.M. Badescu 12]K.B. Barber 38]J. Bäuml 38]C. Baus 89]J.J. Beatty 35]K.H. Becker 12]J.A. Bellido 32]C. Berat 1]X. Bertou 39]P.L. Biermann 31]P. Billoir 73]F. Blanco 31]M. Blanco 35]C. Bleve 38, 36]H. Blümer 27]M. Boháčová 53]D. Boncioli 23]C. Bonifazi 54]R. Bonino 63]N. Borodai 81]J. Brack 66]I. Brancus 65]P. Brogueira 82]W.C. Brown 42]P. Buchholz 75]A. Bueno 45]M. Buscemi 56, 76, 90]K.S. CaballeroMora 44]B. Caccianiga 31]L. Caccianiga 46]M. Candusso 39]L. Caramete 47]R. Caruso 54]A. Castellina 49]G. Cataldi 65]L. Cazon 48]R. Cester 57]A.G. Chavez 90]S.H. Cheng 54]A. Chiavassa 18]J.A. Chinellato 27]J. Chudoba 45]M. Cilmo 12]R.W. Clay 49]G. Cocciolo 45]R. Colalillo 44]L. Collica 49]M.R. Coluccia 65]R. Conceição 9]F. Contreras 12]M.J. Cooper 90]S. Coutu 79]C.E. Covault 90]A. Criss 91]J. Cronin 39]A. Curutiu 34, 33]R. Dallier 18]B. Daniel 5, 3]S. Dasso 36]K. Daumiller 12]B.R. Dawson 24]R.M. de Almeida 47]M. De Domenico 59, 61]S.J. de Jong 23]J.R.T. de Mello Neto 49]I. De Mitri 24]J. de Oliveira 16]V. de Souza 74]L. del Peral 29]O. Deligny 36]H. Dembinski 86]N. Dhital 46]C. Di Giulio 50]A. Di Matteo 86]J.C. Diaz 18]M.L. Díaz Castro 97]P.N. Diep 65]F. Diogo 18]C. Dobrigkeit 60]W. Docters 58]J.C. D’Olivo 97, 29]P.N. Dong 81]A. Dorofeev 36]Q. Dorosti Hasankiadeh 4]M.T. Dova 27]J. Ebr 36]R. Engel 40]M. Erdmann 42]M. Erfani 83, 18]C.O. Escobar 65]J. Espadanal 8, 11]A. Etchegoyen 91]P. Facal San Luis 59, 62, 61]H. Falcke 91]K. Fang 87]G. Farrar 18]A.C. Fauth 83]N. Fazzini 79]A.P. Ferguson 23]M. Fernandes 86]B. Fick 8]J.M. Figueira 8]A. Filevich 70, 71]A. Filipčič 92]B.D. Fox 69]O. Fratu 42]U. Fröhlich 38]B. Fuchs 91]T. Fuji 31]R. Gaior 7]B. García 76]S.T. Garcia Roca 30]D. GarciaGamez 73]D. GarciaPinto 47]G. Garilli 75]A. Gascon Bravo 34]F. Gate 37]H. Gemmeke 31]P.L. Ghia 23]U. Giaccari 44]M. Giammarchi 64]M. Giller 40]C. Glaser 83]H. Glass 4]F. Gomez Albarracin 1]M. Gómez Berisso 10]P.F. Gómez Vitale 65]P. Gonçalves 38]J.G. Gonzalez 81]B. Gookin 54]A. Gorgi 92]P. Gorham 17]P. Gouffon 59, 61]S. Grebe 89]N. Griffith 53]A.F. Grillo 12]T.D. Grubb 3]Y. Guardincerri 45]F. Guarino 19]G.P. Guedes 4]P. Hansen 1]D. Harari 12]T.A. Harrison 81]J.L. Harton 36]A. Haungs 40]T. Hebbeker 36]D. Heck 42]P. Heimann 36]A.E. Herve 12]G.C. Hill 83]C. Hojvat 91]N. Hollon 36]E. Holt 42, 63]P. Homola 59, 61]J.R. Hörandel 28]P. Horvath 28, 27]M. Hrabovský 38]D. Huber 36]T. Huege 47]A. Insolia 67]P.G. Isar 96] K. Islo 35]I. Jandt 59, 61]S. Jansen 4]C. Jarne 8]M. Josebachuili 35]A. Kääpä 38]O. Kambeitz 35]K.H. Kampert 83]P. Kasper 38]I. Katkov 30]B. Kégl 36]B. Keilhauer 85]A. Keivani 18]E. Kemp 86]R.M. Kieckhafer 36]H.O. Klages 37]M. Kleifges 9]J. Kleinfeller 40]R. Krause 35]N. Krohm 37]O. Krömer 35]D. KruppkeHansen 40]D. Kuempel 37]N. Kunka 52]G. La Rosa 79]D. LaHurd 54]L. Latronico 94]R. Lauer 40]M. Lauscher 34]P. Lautridou 32]S. Le Coz 14]M.S.A.B. Leão 32]D. Lebrun 83]P. Lebrun 22]M.A. Leigui de Oliveira 31]A. LetessierSelvon 29]I. LhenryYvon 38]K. Link 55]R. López 76]A. Lopez Agüera 32]K. Louedec 75]J. Lozano Bahilo 35, 77]L. Lu 8]A. Lucero 38]M. Ludwig 23]H. Lyberis 52]M.C. Maccarone 12]M. Malacari 54]S. Maldera 34]J. Maller 27]D. Mandat 83]P. Mantsch 4]A.G. Mariazzi 34]V. Marin 75]I.C. Mariş 49]G. Marsella 49]D. Martello 34, 33]L. Martin 56]H. Martinez 55]O. Martínez Bravo 29] D. Martraire 3]J.J. Masías Meza 36]H.J. Mathes 35]S. Mathys 94]A.J. Matthews 85]J. Matthews 46]G. Matthiae 38]D. Maurel 13]D. Maurizio 80]E. Mayotte 83]P.O. Mazur 80]C. Medina 58]G. MedinaTanco 38]M. Melissas 8]D. Melo 48]E. Menichetti 37]A. Menshikov 60]S. Messina 92]R. Meyhandan 25]S. Mićanović 6]M.I. Micheletti 40]L. Middendorf 73]I.A. Minaya 44]L. Miramonti 66]B. Mitrica 75]L. MolinaBueno 1]S. Mollerach 91]M. Monasor 30]D. Monnier Ragaigne 32]F. Montanet 54]C. Morello 4]J.C. Moreno 90]M. Mostafá 22]C.A. Moura 18, 21]M.A. Muller 40]G. Müller 31]M. Münchmeyer 48]R. Mussa 54 ‡]G. Navarra 75]S. Navas 27]P. Necesal 58]L. Nellen 59, 61]A. Nelles 35]J. Neuser 76 d]D. Newton 42]M. Niechciol 35]L. Niemietz 40]T. Niggemann 86]D. Nitz 26]D. Nosek 26]V. Novotny 28]L. Nožka 42]L. Ochilo 91]A. Olinto 65]M. Oliveira 76]V.M. OlmosGilbaja 73]M. Ortiz 74]N. Pacheco 18]D. Pakk SelmiDei 27]M. Palatka 2]J. Pallotta 38]N. Palmieri 35]P. Papenbreer 76]G. Parente 76]A. Parra 72]S. Pastor 96, 88]T. Paul 27]M. Pech 63]J. Pȩkala 55] R. Pelayo 20]I.M. Pepe 49]L. Perrone 43]R. Pesce 93]E. Petermann 40]C. Peters 50, 51]S. Petrera 43]A. Petrolini 81]Y. Petrov 3]R. Piegaia 36]T. Pierog 3] P. Pieroni 65]M. Pimenta 47]V. Pirronello 8]M. Platino 40]M. Plum 36]A. Porcelli 63]C. Porowski 91]P. Privitera 27]M. Prouza 1]V. Purrello 2]E.J. Quel 35]S. Querchfeld 79]S. Quinn 35]J. Rautenberg 34]O. Ravel 8]D. Ravignani 34]B. Revenu 27]J. Ridky 52, 76]S. Riggi 42]M. Risse 2]P. Ristori 50]V. Rizi 87]J. Roberts 76]W. Rodrigues de Carvalho 76]I. Rodriguez Cabo 46, 76]G. Rodriguez Fernandez 9]J. Rodriguez Rojo 74]M.D. RodríguezFrías 74]G. Ros 73]J. Rosado 28]T. Rossler 36]M. Roth 1]E. Roulet 5]A.C. Rovero 37]C. Rühle 12]S.J. Saffi 66]A. Saftoiu 29]F. Salamida 55]H. Salazar 90]F. Salesa Greus 46]G. Salina 8]F. Sánchez 75]P. SanchezLucas 65]C.E. Santo 65]E. Santos 17]E.M. Santos 80]F. Sarazin 35]B. Sarkar 65]R. Sarmento 9]R. Sato 40]N. Scharf 49]V. Scherini 36]H. Schieler 41]P. Schiffer 37]A. Schmidt 60]O. Scholten 92, 59, 61]H. Schoorlemmer 27]P. Schovánek 36]A. Schulz 59]J. Schulz 4]S.J. Sciutto 52]A. Segreto 31]M. Settimo 85]A. Shadkam 13]R.C. Shellard 1]I. Sidelnik 41]G. Sigl 68]O. Sima 64]A. Śmiałkowski 36]R. Šmída 93]G.R. Snow 90]P. Sommers 12]J. Sorokin 9]R. Squartini 88]Y.N. Srivastava 71]S. Stanič 89]J. Stapleton 63]J. Stasielak 40]M. Stephan 32]A. Stutz 8]F. Suarez 29]T. Suomijärvi 5]A.D. Supanitsky 85]M.S. Sutherland 88]J. Swain 64]Z. Szadkowski 36]M. Szuba 1]O.A. Taborda 8]A. Tapia 32]M. Tartare 97]N.T. Thao 18] V.M. Theodoro 3]J. Tiffenberg 61, 59]C. Timmermans 15]C.J. Todero Peixoto 66]G. Toma 36]L. Tomankova 65]B. Tomé 48]A. Tonachini 76]G. Torralba Elipe 34]D. Torres Machado 27]P. Travnicek 47]E. Trovato 76]M. Tueros 36]R. Ulrich 36]M. Unger 40]M. Urban 58]J.F. Valdés Galicia 76]I. Valiño 45]L. Valore 59]G. van Aar 60]A.M. van den Berg 59]S. van Velzen 41]A. van Vliet 55]E. Varela 58]B. Vargas Cárdenas 92]G. Varner 73]J.R. Vázquez 76]R.A. Vázquez 30]D. Veberič 46]V. Verzi 27]J. Vicha 8]M. Videla 57]L. Villaseñor 96]B. Vlcek 71]S. Vorobiov 4]H. Wahlberg 8, 11]O. Wainberg 40] D. Walz 77]A.A. Watson 37]M. Weber 40]K. Weidenhaupt 36]A. Weindl 38]F. Werner 90]B.J. Whelan 88]A. Widom 80]L. Wiencke 63 ‡]B. Wilczyńska 63]H. Wilczyński 36]M. Will 91]C. Williams 40]T. Winchen 35]D. Wittkowski 8]B. Wundheiler 59]S. Wykes 91 a]T. Yamamoto 86]T. Yapici 84]P. Younk 85]G. Yuan 42]A. Yushkov 75]B. Zamorano 76]E. Zas 71, 70]D. Zavrtanik 70, 71]M. Zavrtanik 87 c]I. Zaw 56 b]A. Zepeda 91]J. Zhou 37]Y. Zhu 18]M. Zimbres Silva 42]M. Ziolkowski
Reconstruction of inclined air showers detected with the Pierre Auger Observatory
The Pierre Auger Collaboration

Department of Physics and Astronomy, Lehman College, City University of New York, New York, USA

Centro Atómico Bariloche and Instituto Balseiro (CNEAUNCuyoCONICET), San Carlos de Bariloche, Argentina

Centro de Investigaciones en Láseres y Aplicaciones, CITEDEF and CONICET, Argentina

Departamento de Física, FCEyN, Universidad de Buenos Aires y CONICET, Argentina

IFLP, Universidad Nacional de La Plata and CONICET, La Plata, Argentina

Instituto de Astronomía y Física del Espacio (CONICETUBA), Buenos Aires, Argentina

Instituto de Física de Rosario (IFIR)  CONICET/U.N.R. and Facultad de Ciencias Bioquímicas y Farmacéuticas U.N.R., Rosario, Argentina

Instituto de Tecnologías en Detección y Astropartículas (CNEA, CONICET, UNSAM), and National Technological University, Faculty Mendoza (CONICET/CNEA), Mendoza, Argentina

Instituto de Tecnologías en Detección y Astropartículas (CNEA, CONICET, UNSAM), Buenos Aires, Argentina

Observatorio Pierre Auger, Malargüe, Argentina

Observatorio Pierre Auger and Comisión Nacional de Energía Atómica, Malargüe, Argentina

Universidad Tecnológica Nacional  Facultad Regional Buenos Aires, Buenos Aires, Argentina

University of Adelaide, Adelaide, S.A., Australia

Centro Brasileiro de Pesquisas Fisicas, Rio de Janeiro, RJ, Brazil

Faculdade Independente do Nordeste, Vitória da Conquista, Brazil

Universidade de São Paulo, Escola de Engenharia de Lorena, Lorena, SP, Brazil

Universidade de São Paulo, Instituto de Física, São Carlos, SP, Brazil

Universidade de São Paulo, Instituto de Física, São Paulo, SP, Brazil

Universidade Estadual de Campinas, IFGW, Campinas, SP, Brazil

Universidade Estadual de Feira de Santana, Brazil

Universidade Federal da Bahia, Salvador, BA, Brazil

Universidade Federal de Pelotas, Pelotas, RS, Brazil

Universidade Federal do ABC, Santo André, SP, Brazil

Universidade Federal do Rio de Janeiro, Instituto de Física, Rio de Janeiro, RJ, Brazil

Universidade Federal Fluminense, EEIMVR, Volta Redonda, RJ, Brazil

Rudjer Bošković Institute, 10000 Zagreb, Croatia

Charles University, Faculty of Mathematics and Physics, Institute of Particle and Nuclear Physics, Prague, Czech Republic

Institute of Physics of the Academy of Sciences of the Czech Republic, Prague, Czech Republic

Palacky University, RCPTM, Olomouc, Czech Republic

Institut de Physique Nucléaire d’Orsay (IPNO), Université Paris 11, CNRSIN2P3, Orsay, France

Laboratoire de l’Accélérateur Linéaire (LAL), Université Paris 11, CNRSIN2P3, France

Laboratoire de Physique Nucléaire et de Hautes Energies (LPNHE), Universités Paris 6 et Paris 7, CNRSIN2P3, Paris, France

Laboratoire de Physique Subatomique et de Cosmologie (LPSC), Université GrenobleAlpes, CNRS/IN2P3, France

Station de Radioastronomie de Nançay, Observatoire de Paris, CNRS/INSU, France

SUBATECH, École des Mines de Nantes, CNRSIN2P3, Université de Nantes, France

Bergische Universität Wuppertal, Wuppertal, Germany

Karlsruhe Institute of Technology  Campus North  Institut für Kernphysik, Karlsruhe, Germany

Karlsruhe Institute of Technology  Campus North  Institut für Prozessdatenverarbeitung und Elektronik, Karlsruhe, Germany

Karlsruhe Institute of Technology  Campus South  Institut für Experimentelle Kernphysik (IEKP), Karlsruhe, Germany

MaxPlanckInstitut für Radioastronomie, Bonn, Germany

RWTH Aachen University, III. Physikalisches Institut A, Aachen, Germany

Universität Hamburg, Hamburg, Germany

Universität Siegen, Siegen, Germany

Dipartimento di Fisica dell’Università and INFN, Genova, Italy

Università di Milano and Sezione INFN, Milan, Italy

Università di Napoli ”Federico II” and Sezione INFN, Napoli, Italy

Università di Roma II ”Tor Vergata” and Sezione INFN, Roma, Italy

Università di Catania and Sezione INFN, Catania, Italy

Università di Torino and Sezione INFN, Torino, Italy

Dipartimento di Matematica e Fisica ”E. De Giorgi” dell’Università del Salento and Sezione INFN, Lecce, Italy

Dipartimento di Scienze Fisiche e Chimiche dell’Università dell’Aquila and INFN, Italy

Gran Sasso Science Institute (INFN), L’Aquila, Italy

Istituto di Astrofisica Spaziale e Fisica Cosmica di Palermo (INAF), Palermo, Italy

INFN, Laboratori Nazionali del Gran Sasso, Assergi (L’Aquila), Italy

Osservatorio Astrofisico di Torino (INAF), Università di Torino and Sezione INFN, Torino, Italy

Benemérita Universidad Autónoma de Puebla, Puebla, Mexico

Centro de Investigación y de Estudios Avanzados del IPN (CINVESTAV), México, Mexico

Universidad Michoacana de San Nicolas de Hidalgo, Morelia, Michoacan, Mexico

Universidad Nacional Autonoma de Mexico, Mexico, D.F., Mexico

IMAPP, Radboud University Nijmegen, Netherlands

KVI  Center for Advanced Radiation Technology, University of Groningen, Netherlands

Nikhef, Science Park, Amsterdam, Netherlands

ASTRON, Dwingeloo, Netherlands

Institute of Nuclear Physics PAN, Krakow, Poland

University of Łódź, Łódź, Poland

Laboratório de Instrumentação e Física Experimental de Partículas  LIP and Instituto Superior Técnico  IST, Universidade de Lisboa  UL, Portugal

’Horia Hulubei’ National Institute for Physics and Nuclear Engineering, Bucharest Magurele, Romania

Institute of Space Sciences, Bucharest, Romania

University of Bucharest, Physics Department, Romania

University Politehnica of Bucharest, Romania

Experimental Particle Physics Department, J. Stefan Institute, Ljubljana, Slovenia

Laboratory for Astroparticle Physics, University of Nova Gorica, Slovenia

Institut de Física Corpuscular, CSICUniversitat de València, Valencia, Spain

Universidad Complutense de Madrid, Madrid, Spain

Universidad de Alcalá, Alcalá de Henares (Madrid), Spain

Universidad de Granada and C.A.F.P.E., Granada, Spain

Universidad de Santiago de Compostela, Spain

School of Physics and Astronomy, University of Leeds, United Kingdom

Case Western Reserve University, Cleveland, OH, USA

Colorado School of Mines, Golden, CO, USA

Colorado State University, Fort Collins, CO, USA

Colorado State University, Pueblo, CO, USA

Fermilab, Batavia, IL, USA

Los Alamos National Laboratory, Los Alamos, NM, USA

Louisiana State University, Baton Rouge, LA, USA

Michigan Technological University, Houghton, MI, USA

New York University, New York, NY, USA

Northeastern University, Boston, MA, USA

Ohio State University, Columbus, OH, USA

Pennsylvania State University, University Park, PA, USA

University of Chicago, Enrico Fermi Institute, Chicago, IL, USA

University of Hawaii, Honolulu, HI, USA

University of Nebraska, Lincoln, NE, USA

University of New Mexico, Albuquerque, NM, USA

University of Wisconsin, Madison, WI, USA

University of Wisconsin, Milwaukee, WI, USA

Institute for Nuclear Science and Technology (INST), Hanoi, Vietnam

Deceased

Now at Konan University

Also at the Universidad Autonoma de Chiapas on leave of absence from Cinvestav

Now at NYU Abu Dhabi

Now at University of Liverpool
Keywords: Pierre Auger Observatory, ultrahigh energy cosmic rays, inclined extensive air showers, shower reconstruction
Contents
1 Introduction
The Pierre Auger Observatory is a hybrid instrument combining an array of particle detectors, the Surface Detector (SD) array, to sample the air shower front as it reaches the ground and Fluorescence Detector (FD) telescopes to capture the ultraviolet light emitted by the nitrogen as showers develop in the atmosphere. The FD is used to monitor the atmosphere, on dark clear nights, above the 3000 km area over which the SD is laid out. The site is located near the town of Malargüe, in the Argentinian province of Mendoza, at an altitude of about 1400 m above sea level and at an average latitude of S [1].
The SD stations are waterCherenkov detectors which are sensitive to inclined particles, so that the array is also sensitive to very inclined showers. Since the beginning of its deployment the SD of the Pierre Auger Observatory has been routinely recording events with zenith angles up to . However, the reconstruction of events with zenith angles exceeding requires a different method to the one used for events nearer the vertical, due to an asymmetry induced in the lateral distribution of the shower particles by the geomagnetic field.
The showers produced by cosmicray hadrons at large zenith angles traverse much larger atmospheric depths than vertical showers and, as a result, their shower maximum occurs higher in the atmosphere. By the time the shower front reaches the ground, the generation of electrons and photons in the main cascading process is basically finished, and the bulk of them have been absorbed in the atmosphere. Most of the particles that reach the ground are energetic muons from charged pion decays in the showering process accompanied by a smaller electron and photon component that stems from the muons themselves, primarily from decays in flight. These muons travel long distances deviating in the geomagnetic field, so that when they reach the ground the characteristic cylindrical symmetry of the showers is lost.
Although the interest in inclined showers dates to the early days of extensive air shower measurements [2], it was only in the beginning of the 2000s, when the patterns of the muons at the ground level were sufficiently understood [3], that methods were developed to reconstruct data in a reliable way [4]. The analysis of these showers is of particular interest because it enhances the exposure of the detector by 30% and extends the sky coverage to regions that are otherwise unobservable. It provides the basis of a completely independent measurement of the cosmic ray spectrum above eV [5, 6, 7, 8, 9], to be updated in a forthcoming publication. Since inclined showers are mainly composed of muons, their study provides an almost uncontaminated measurement of the muon content of the shower, while for events with zenith angles less than about the muonic component can not be disentangled from the electromagnetic activity, and thus must be inferred indirectly in the absence of shielded detectors, which is intrinsically more challenging. Therefore, the inclined data also provide complementary information to constrain the nature of the arriving particles [10, 11]. In addition, they also constitute the background against which the search for highenergy neutrinos with inclined showers must be made [12, 13].
This article deals with inclined showers as follows. Section 2 describes details of the Surface Detector of the Pierre Auger Observatory that are of relevance for inclinedshower reconstruction. Section 3 deals with the modeling needed for the reconstruction procedure, namely the muon distribution at the ground level, the treatment of the electromagnetic contribution to the signal and also the signal response of the surface detectors to the passage of muons. In section 4 the details of the reconstruction procedure are described and the uncertainties addressed. The energy calibration procedure is discussed in section 5, and section 6 summarizes the conclusions.
2 Inclined events in the Auger Observatory
The SD consists of more than 1600 waterCherenkov stations arranged on a hexagonal grid, each station at a distance of 1.5 km from its six nearest neighbors. Each station is a cylindrical water tank of 10 m surface area, filled to 1.2 m with purified water contained inside a diffuselyreflective liner. The volume of water is viewed by three 9 inch photomultipliers that detect the Cherenkov light emitted during the passage of charged particles. The signal is digitized in time slots of 25 ns using a Flash AnalogtoDigital Converter (FADC) running at 40 MHz. All the stations are controlled remotely, and data are transmitted from the detectors to a central station by Local Area Network radio links. Synchronization to 8 ns of relative precision is provided by commercial GPS systems. A full description of the SD can be found in [1]. The stations are calibrated online continuously by identifying the maximum in the raw signal histograms of signals produced by atmospheric background muons which are sampled at regular intervals. The integrated charge signal of this maximum is related to that of a vertical muon traversing the detector through its center, Vertical Equivalent Muon or VEM, which provides the signal unit with 3% accuracy [14]. The total signal at each station in VEM units, , is obtained by integrating the signal traces in time. Examples of FADC traces in VEM units are shown in figure 1.
When shower electrons and photons ( MeV) reach the surface detectors they are absorbed in the water and provide a Cherenkov light signal which is approximately proportional to the total deposited energy. Shower muons are more energetic ( MeV) and travel right through the tanks, typically giving signals proportional to their track length. For vertical showers initiated by protons or nuclei, the contribution of electrons and photons to the signal is often comparable or even larger than that due to the muons [15]. However, for inclined showers muons dominate the signal due to the electromagnetic component being largely absorbed by the atmosphere. Moreover the detector increases its relative response to muons with respect to electrons and photons as the zenith angle increases since the muon signal scales with the track length and the diameter of the tank exceeds its height.
The FD is distributed in four buildings on the perimeter of the surface array. Each building contains six telescopes that together cover in azimuth and nearly in elevation. The FD is fully described in [16]. Whereas the SD measurements are performed with a duty cycle of almost 100%, the FD only operates on clear moonless nights and has a duty cycle of 13%. The FD allows a calorimetric measurement of the shower energy deposited in the atmosphere, in contrast to the SD. The positions of the triggered pixels and the arrival time of the light are used to extract the shower direction. The profile function of the energy deposited in the atmosphere is determined from the signals at the triggered pixels [17] after taking into account the separate fluorescence and Cherenkov light contributions, as well as light attenuation and dispersion in the atmosphere. The atmospheric conditions are monitored regularly using several techniques to provide relevant data on attenuation and light dispersion parameters [18]. The electromagnetic energy released by the shower in the atmosphere is obtained by fitting the longitudinal profile to a GaisserHillas function [19] and integrating over the range of atmospheric depths. The total energy of the primary particle is derived from the calorimetric energy by adding the invisible energy which accounts for the energy carried by penetrating particles. Most of the detected FD events also induce triggers of at least one SD station and are called hybrid events. The geometry of the shower axis can be determined more precisely for these events by using timing information from the FD pixels, coupled with the arrival time of the shower at the SD station with highest signal. A fraction of these hybrid events, called “golden hybrid” events, have sufficient SD stations triggered to allow for an independent reconstruction with SD techniques.
2.1 Selection of inclined events with the Surface Detector
The SD trigger (described in detail in [20]) has been designed in a hierarchical pattern to identify a cosmicray event and reject random coincidences. The first two triggers, T1 and T2, apply at station level. T1 and T2 require that the pulse height of the signal exceeds two preset values (larger for T2 than for T1). Alternatively a second T2 trigger can also be obtained with smaller pulses that are spread in time (T2ToT), which is a characteristic of the electromagnetic component in showers of near vertical incidence (see left panel of figure 1). The T1 trigger is stored temporarily while the data of all T2 triggers are sent to the central station. T2 triggered stations that lie in a close and compact configuration are checked for coincidences in time. If sufficient coincident stations are found, the third level trigger (T3) is passed and the acquisition process starts. In this process the T1 triggers stored at the local stations are also downloaded and added to the event.
Starting with the T3 triggers, a higherlevel trigger hierarchy is implemented to satisfy the precision requirements of the consequent analysis. The T4 and T5 filters are respectively applied offline for selecting physical events, and to ensure that they fall on a region of the array where the surface detectors were operational at the time of the event, to guarantee their quality. They are defined differently for inclined events. The “inclined T4” condition ensures that the start time of the signals^{1}^{1}1Start time is defined as the arrival time of the first shower particle into the SD station. of at least four nearby T2 triggered stations are compatible with a shower front moving at the speed of light. Accidental triggers are removed by eliminating stations that have times outside the window corresponding to the passage of the shower front, as determined by the rest of the stations. Stations are eliminated one by one (selecting the ones with the highest timing offsets) starting from the T3 trigger selection, until a satisfactory configuration with four or more stations in a compact arrangement is found. This removes a large number of showers reconstructed with incorrect arrival directions due to random coincidences. Events that have stations which are all aligned (according to the hexagonal pattern of the array) are excluded since the arrival direction cannot be reconstructed accurately from these starttime data.
The inclined T5 condition is defined so as to select events that fall inside a region of the array in which all stations were operational (“active”), requiring that the station that is nearest to the reconstructed core and its six adjacent stations are all operational. This active region defines an “active unit cell” shaped as a hexagon and forms the basis of the array aperture calculation [20]. The T5 condition serves two purposes: (i) it avoids the reconstruction of events that fall near to the edge of the array or in regions where a station is temporarily not fully operational, which can have large uncertainties, and (ii) it ensures that the shower falls inside the active area that has been considered to establish the aperture. Note that the T5 condition is the basic trigger of the selection criteria used for the subsequent analysis of inclined showers.
2.2 General characteristics of inclined SD events
As the zenith angle of the showers increases, and they traverse larger atmospheric depths, the showers that reach the ground are increasingly attenuated. It appears that this would result in a reduced efficiency to detect inclined showers with the SD with increasing zenith angle. On the contrary, there is a compensating geometric effect since the density of stations in the shower plane (perpendicular to the shower axis) increases rapidly in proportion to . The signal reduction due to the muon attenuation over long path lengths is relatively small and is compensated by some of the surface detectors being closer to the shower axis where the signal is higher (compared to more vertical events) and thus more likely to trigger. Inclined events can have a large number of triggered stations as illustrated in figure 2.
This effect manifests itself when considering the number of triggered stations of events with a given energy as a function of zenith angle. For instance, considering an energy around 10 EeV (the energy at which the SD is fully efficient for events out to ), a vertical event has on average 8 triggered stations. This number is around 13 for a event of the same energy, while at it increases to 25. This effect is illustrated in figure 3 which shows the average number of stations as a function of shower energy at a fixed zenith angle (left panel) and the average number of stations as a function of zenith angle at a fixed energy (right panel).
3 Modeling particle distributions and detector responses
Showers with require specific reconstruction methods since they are dominated by muons displaying complex density patterns at the ground level. The method used for the reconstruction of inclined showers is based on a fit of the measured signals to the expected pattern, and requires modeling of a twodimensional distribution of the muon number densities at the ground, the response of the detectors to the passage of muons, and the treatment of the electromagnetic component of the signal in the detectors.
3.1 Number density of muons at the ground
The number of muons per unit area, i.e. the muon density , is a two dimensional function of position coordinates relative to the shower axis position , that is , projected onto the shower plane. Descriptions of at the ground level have been obtained and understood with the aid of comprehensive simulations of extensive air showers [3, 21].
As muons traverse the atmosphere they lose energy by ionization and hard muon interactions, namely bremsstrahlung, pair production and nuclear interactions via photonuclear processes. Below the critical energy for the muons, which is of order 500 GeV, ionization losses dominate. Decay probability in flight becomes important by effectively removing muons below an energy threshold which depends on both the muon energy and the distance traveled [22]. As the zenith angle increases from to , the distance traveled to the ground level increases from km to values over an order of magnitude larger, and the average energy of the muons at the ground level rises accordingly. The average energy of muons at the ground level, produced by a primary hadron at zenith angle (), is about 25 (60) GeV.
In the Earth’s magnetic field , the muons are also deflected by the Lorentz force which laterally separates the positively and negatively charged muons [2]. Typical asymmetries are illustrated in Figs. 4 and 5. The separation proceeds in the direction perpendicular to the plane defined by the shower axis and the magnetic field. The magnitude of the separation depends mainly on the component of perpendicular to the shower direction, the muon energy and the distance traveled [3]. The resulting signal patterns at the ground thus have a quite strong dependency on the arrival directions. As the zenith angle changes, the large variations in the distance traveled by the muons are to a large extent responsible for changes in the patterns at the ground level. As the azimuthal direction of the shower changes, smaller differences in the patterns are also observed, due to the varying angle between the typical muon velocity and the field. Two different approaches have been used to obtain these distributions. One is based on a transformation of cylindrically symmetric patterns, exploiting the anticorrelation between muon energy and angle to the shower axis [3]. The other relies on continuous parameterizations in zenith angle and position in the shower plane that are fitted to results obtained from simulations [21]. Both approaches have been shown to reproduce the average profile of a given set of simulated showers with an accuracy better than 5%.
When the arrival direction and the nature of the primary particle are fixed, the muon number density has been shown to scale nearly linearly with shower energy ( with typically being in the range –) [3, 21]. There are some differences between the distributions depending on the assumed nature of the primary particle, its energy and the hadronic interaction model used in the simulations. It has also been shown that these differences are manifested primarily by an overall normalization of the muon densities, and the shapes of these functions are approximately the same for a given arrival direction [23], i.e., weakly dependent on both shower energy and composition. Both characteristics are illustrated in figure 4.
The universal shape of the muon distribution and the scaling between muon number density and shower energy provide the basis of the fitting procedure. The reconstruction of the shower size is based on the fit of measured signals to the expected muon patterns. Details are fully described below.
According to the scaling property mentioned above, the expected muon number density at the ground can be written as:
(3.1) 
Here is a measure of the shower size and is the relative normalization of a particular event with respect to a reference muon distribution, , conventionally chosen to be the average muon density for primary protons with eV obtained with a chosen shower model, QGSJetII03 [24] reference in our case. The dependence of these functions on the zenith and azimuth angles is indicated explicitly.
Two sets of muon distributions were generated for comparison purposes, one following [3] using the Aires [25] package for shower simulations with the QGSJet01 [26] model for hadronic interactions, and the other following [21] with Corsika [27] and QGSJetII03. Examples of such distributions are shown in figure 5 for different zenith angles.
The actual value obtained for depends upon the particular choice of composition and hadronic model made for the reference distribution. The proton showers obtained with the QGSJetII03 model was chosen as a reference for the data analysis presented in this work. For instance, the muon densities for proton showers derived using different highenergy hadronic interaction models scale with the approximate factors relative to the reference distribution given in table 1. A primary composition different from protons would also enhance the shower size scaling with the muon content [32]. In the extreme case of iron, the corresponding approximate factors are also given in table 1. These factors are indicative of the large expected uncertainties associated with the hadronic models and the unknown composition. Nevertheless, these uncertainties do not have a large impact on the measurement of the energy spectrum where the energy scale had been inferred from a subsample of events measured simultaneously with the FD and SD (see section 5). This analysis mimics the procedure to provide an absolute energy calibration [33] used for the reconstruction of events with zenith angle less than . Most of the uncertainties associated with the unknown primary composition and hadronic model, as well as many uncertainties associated with the reconstruction, are absorbed in this robust and reliable calibration procedure.
Model  Proton  Iron 

QGSJet01  1.10  1.46 
QGSJetII03  1.00  1.32 
QGSJetII04  1.19  1.58 
Epos 1.99  1.25  1.65 
Epos LHC  1.22  1.61 
Sibyll 2.1  0.90  1.20 
3.2 Electromagnetic component
As inclined showers develop in the atmosphere, the EM cascade produced by the neutral pions in the hadronic collisions is attenuated, leaving a front dominated by muons. The remaining EM component (electrons, positrons and photons) originates from two different mechanisms. First is due to the tail of the hadronic cascade and decreases very rapidly as the zenith angle increases. Its magnitude depends on composition and interaction models. This component increases with the primary energy and varies according to the fluctuations in shower maximum. This dependence has to be considered as a source of systematic uncertainty in the electromagnetic correction, as will be discussed in section 4.2.3. The second component is produced by the muons themselves and closely follows the muon density patterns. It is mainly due to the showering of the electrons from muon decay in flight, although highenergy muons very close to the core also contribute through pair production and bremsstrahlung. More details can be found in [15].
In the reconstruction procedure (see section 4) the signals measured with the SD are compared to the reference muon distribution, and for this purpose it is necessary to extract the signal induced by muons from the total signal at each detector. This is achieved by subtracting the EM component from the detector signal using the average ratio, , of the electromagnetic to the muonic contributions:
(3.2) 
This ratio was studied using MonteCarlo simulations. Parameterizations of the average electromagnetic () and muonic () contributions to the signal were obtained separately in terms of distance to the shower axis , and of zenith angle . Proton simulations at eV with QGSJet01 was chosen as a reference. The accuracy of this parameterization is better than 5%.
In figure 6 the ratio , averaged over the polar angle (with respect to the shower axis projected onto the shower plane), is shown as a function of for different . The ratio can be seen to decrease as increases from to . At larger and near the shower axis, can be seen to increase slightly due to hard muon interaction processes. At distances from the shower axis exceeding 1 km, changes only weakly with distance, typically lying between 15% and 30%.
3.3 Signal response of the surface detector stations
To reconstruct the position of the shower core and shower size , the measured signals are fitted to the model predictions. The fit requires the evaluation of the probability density function (PDF) of the signal response of each surface detector to the expected number of muons. The PDF distributions are functions of the signal deposited by muons in the detector, which is in turn estimated from the measured signal , using the average ratio of electromagnetictomuonic signals described in the previous section (eq. (3.2)):
(3.3) 
The PDF for each surface detector is constructed from the expected number of muons, assuming Poisson statistics and using a simpler PDF corresponding to the passage of a fixed number of muons.
The basic prerequisite for these probability densities is the signal distribution for a detector hit by a single muon. These were obtained with high statistics using a module in the official software framework [34] of the Pierre Auger Observatory, which simulates the detector response and interfaces it with the Geant4 package [35]. Multiple histograms for the signal response were generated for discrete values of , and different relative positions of the station with respect to the shower core [36]. They are illustrated in Figs. 7 and 8. The signal is mainly due to the Cherenkov light emission by the muon tracks, and the PDF is closely related to the tracklength distributions of the muons inside the detectors, which have a strong dependence on the zenith angle, as explicitly shown in figure 7.
The simulations include also contributions that do not scale with the muon track length and are dependent on the relative position of the stations. Delta rays (i.e., scattered electrons within the detector) account for an enhancement of the signal of order 20%, which increases with energy [4]. For highenergy muons bremsstrahlung, pair production and nuclear interactions inside the detector [4] give contributions which appear as harder tails of the response functions for stations close to the shower core. Lowenergy muons have a reduced Cherenkov efficiency (enhanced by energy loss), becoming zero at the Cherenkov threshold. The average signal increases as the distance to the shower axis is reduced (as illustrated in figure 8) due to the increase in muon energy. The contribution, to the average signal, of the direct Cherenkov light that hits the photomultipliers without any reflection from the walls of the detector, ranges from 3% at to 10% at .
Response histograms for single muons are used to obtain the response of the detector to the passage of multiple () muons by convolution in an iterative fashion, when is between 1 and 8. When a Gaussian approximation is used with an average and a standard deviation of , in terms of the averages () and standard deviations () of the corresponding single muon histograms. As the number of muons grows, histograms rapidly become Gaussianlike, as anticipated from the Central Limit Theorem, see Fig 9. We label these PDF distributions as:
(3.4) 
explicitly indicating the dependence on , and . The statistical uncertainty of these histograms is below . A total of 12 960 histograms is needed to describe the signal distributions for the passage of muons () in 15 bins of , 12 bins of polar angle and 9 bins of distance .
4 Reconstruction
The reconstruction of inclined showers is performed first by fitting the arrival directions of the events using the start times of the signals, and then fitting the muon density patterns at the ground level to the signals measured by the surface detectors.
The reconstruction was also validated using large samples of simulated events. Extensive libraries of isotropically arriving proton and iron showers were made with a thinning level [37] of . A set of events simulated with Aires and QGSJet01 and a spectral index of in the energy range , and a zenith angle between and (100 000 events each) was used. A second and third set were simulated with Corsika using QGSJet04 or Epos LHC, with a spectral index of between and . Each of these sets covers proton and iron in the same energy interval and is further split into three equal subintervals in (each energy subinterval has 30 000 events). Further libraries were also generated with QGSJetII03 and Epos 1.99. All of these showers have subsequently undergone a full simulation of the detector with random impact points in the SD array, within the framework, to generate databases of simulated events.
Unless otherwise indicated, the criteria used to select the simulated events for comparison match the selection that will be used for the measurement of the spectrum with inclined data. Events with are required to pass the inclined T4 and T5 conditions, with also to ensure the SD array is fully efficient.
4.1 Angular reconstruction
The directions of the incident cosmic rays are determined from the relative arrival times of the shower front in the triggered stations. The start times of the signals are fitted to those expected from a shower front with curvature, using a standard minimization procedure.
The angular accuracy depends on the number of triggered stations, on sizes of their signals, and also on the zenith angle of the shower itself. Typically, as the shower becomes more inclined, the arrival direction is reconstructed better, although this trend is inverted for events above . As the shower becomes dominated by the muons, its front develops a smaller time spread than for vertical showers, and the start time can be established more precisely. In addition, as the zenith angle increases, the curvature of the shower front is reduced since the muons are produced further away from the ground and have higher energies. Finally, the number of triggered stations in a shower tends to increase as the zenith angle rises, due simply to the projection of the station positions into the shower plane.
Equivalent precision on angular reconstruction has been obtained with the different models used to describe the shower front and the variance reported in [22, 38].
The angular reconstruction was tested with simulations. The left panel of figure 10 displays the difference between true and reconstructed angles, as well as the angular resolution as a function of zenith angle for showers simulated with QGSJet01 and for model [22]. The angular resolution, defined as the angle at 68% of the cumulative distribution function of the space angle (angular separation between the true and reconstructed directions), is obtained by fitting this distribution with a Gaussian resolution function, , and corresponds to [39]. The right panel of figure 10 displays the angular resolution as a function of shower size, comparing the performance using different shower simulations and hadronic models. The reconstructed zenith and azimuth angles have a bias less than and , respectively, and the angular resolution better than , improving to better than for the highest energies.
The angular reconstruction of the SD was studied using hybrid data, data from a region of the array with stations on a 750 m grid, events with stations that are duplicated at given positions, and different samples of simulated SD events. These approaches have yielded compatible results.
4.2 Shower size and core position reconstruction
4.2.1 Procedure
Once the arrival direction is established, the expected number of muons at each station can be obtained multiplying the corresponding muon number density (eq. (3.1)) by the detector area projected onto the shower plane:
(4.1) 
Estimates of the position of the shower core and the shower size are obtained by fitting the expected number of muons to the muonic part of the measured signal (eq. (3.3)). The fitting is performed using a maximumlikelihood method including the information from nontriggered and saturated stations. This constrains the shower size and reduces any selection bias due to the threshold trigger. The loglikelihood function is the logarithm of the combined probability to obtain the measured signals in all detectors. For each triggered and nontriggered participating station, this is taken as the product of the probability densities of getting the muonic signal when muons are expected. Therefore, the loglikelihood function is given by:
(4.2) 
and depends on the three free parameters which enter the calculation of (eq. (4.1)), two for the shower core position and one for the shower size .
The measured signal when muons are expected, can be produced by different number of muons , each with different probability density function (see eq. (3.4) and figure 9). To compute each local probability , the sum over all the probabilities for all possible numbers of muons to produce the measured signal has to be computed. Each probability must be weighted by the Poisson probability of muons entering the station when are expected. Then, the overall probability of obtaining in a given detector becomes:
(4.3) 
where the average trigger probability in terms of the signal at each station is also included. , estimated using the signal distributions of the triggered stations in the data, is shown in figure 11.
The sum naturally accounts for Poisson fluctuations in the number of muons but additional signal fluctuations enter through . The infinite sum is truncated for practical reasons.
For stations that have no signal the probability distribution has to be replaced by the probability that the detector does not trigger. The probability that a station does not trigger can be obtained by integrating the detector response functions , weighted by the probability function for not triggering. This is in turn obtained by summing the Poisson probabilities of having muons traversing the station, weighted by the corresponding probabilities for not triggering:
To limit the computing time of the minimization procedure, only detector stations that have no signal and are within four concentric hexagons around all the triggered stations are included in the likelihood maximization procedure. For instance, the number of nontriggering stations entering the fit is on average times the number of stations with signal. The number of stations without signal used for the reconstruction has a minor effect on the best fit value of , and becomes negligible when sufficient stations have been considered. The average change in shower sizes when switching from four to two hexagons are below 5%.
4.2.2 Performance
An example of a reconstructed event from the data is shown in the left panel of figure 12 indicating the position of the shower core and the contour plots of the average muon number density, in relation to the measured signals in the stations. The expected muon densities can be obtained by multiplying the density of the model distribution by the reconstructed value of . The values of the expected and measured signals are shown in the right panel of figure 12 for the same event. For each triggered station, the expected muonic signal is obtained by summing over all the average signal responses for all possible numbers of muons with Poisson weights, . Then, the expected signal is estimated from using the average ratio of electromagnetic and muonic signals: (analogous to eq. (3.3)).
The performance is evaluated by studying the reconstruction of events obtained with full simulation. In the left panel of figure 13 the expected and measured signals are compared for a sample of simulated events. The reconstruction procedure achieves an agreement between expected and measured signals at the 5% level for expected signals above 10 VEM. Below this number, the comparison is not relevant due to the trigger effects and associated upward fluctuations of the signal. This is the reason why the “nontriggered” stations should be considered in the reconstruction, as explained in the previous section 4.2.1.
The reconstruction of the core position depends on the zenith angle, nevertheless this dependency is reduced when converted to the shower plane. The average distance between the reconstructed and true core positions in this plane is 108 m, ranging between 80 to 160 m as the zenith angle increases from to .
The relative uncertainty in as obtained from the fit is shown in the right panel of figure 13 for data and simulated events. It decreases from about 13% at ( EeV) to about at ( EeV).
An approximate estimate of the true value of can be made for each simulated shower and compared to the reconstructed value. The number of muons that have reached the ground is obtained and normalized to the number of muons in the reference distribution (i.e. the integral of in eq. (3.1) over an area),
(4.4) 
The total number of muons in the reference distribution is shown in figure 14 where it is clear that the muon number is attenuated with the zenith angle, whereas, by definition, is independent of the zenith angle. The muon number is essentially independent of azimuth angle. For a given particle species and arrival direction the ratio, , scales with shower energy, on average, and is normalized to eV, so that it can be compared directly with .
The difference between and , and its standard deviation, are illustrated in the left and right panels of figure 15, respectively, for protons in a variety of models. The relative difference is on average less than 6%. This implies that is a good estimator of the number of muons in the shower. The relative difference observed can be interpreted as a bias in the measurement of the number of muons. It can be corrected for if is to be used as a measurement of the number of muons in the shower [40, 41]. However, such a bias has no relevance for the spectrum calculation, since the energy is finally obtained with the calibration procedure using events with their energy determined by the FD. The standard deviation decreases as the shower size increases, since for larger showers there are more stations entering the fit. Relative standard deviation ranges from 16 to 19% for the smaller shower sizes considered () and drops to 5 to 9% for the largest events with .
The presented standard deviation is an upper limit of the uncertainty of the fit to the muon distribution due to the showertoshower fluctuations. The showertoshower fluctuations imply variations in both the total number of muons and the shape of the muon distribution of a given shower with fixed arrival direction, energy and primary mass. If the showers only fluctuated by changing the scale of the muon distribution, the fitted value of would directly reflect the changes of normalization of the showertoshower fluctuations.^{2}^{2}2The total number of muons fluctuates in a broad range between 3 and 28%, depending primarily on composition and to a lesser extent on hadronic model assumptions [42, 6]. The minimum values are obtained with pure iron, and the maximum with a protoniron mixture. But the fluctuations in the shape of the lateral profile are not taken into account in the fit which uses the average profile of the muon distribution. As a result, an unknown part of the showertoshower fluctuations propagates into the standard deviation shown in the right panel of figure 15. This is expected to be relatively small since the standard deviation at the larger shower sizes is well below the value of the fluctuations in for protons, which is about 20%.
In addition, the shown standard deviation includes the angular uncertainty propagated into , nevertheless it has only a small effect.^{3}^{3}3This is since the angular accuracy is quite good. A systematic uncertainty of in the zenith angle reconstruction corresponds to an uncertainty in of 2% (3%) for a shower with the zenith angle of (). This is consistent with the uncertainties in obtained in the fitting procedure and shown in the right panel of figure 13.
The standard deviation is slightly lower (by ) for showers simulated with QGSJetII04 and Epos LHC than for those simulated with QGSJet01. This is possibly due to shower fluctuations and to the fact that the reference muon distributions were simulated with QGSJetII03.
4.2.3 Systematic uncertainties
The EM correction to the detector signal (section 3.2) depends on the primary energy, on the composition and on the hadronic interaction model, and becomes largest in the reconstruction of shower size at zenith angles close to . The actual value obtained for depends upon the particular choice of energy, composition and hadronic model made for the reference EM correction, which correspond to proton simulations at eV with QGSJet01 in this work. The systematic effect in associated with the unknown composition and the hadronic model was evaluated by computing alternative EM corrections (as defined in eq. (3.2)) with simulated showers using different combinations of energy ( and eV), composition (proton and iron) and interaction models (QGSJet01 and Sibyll). A sample of data was reconstructed using these variations of the EM correction, leading to corresponding new values of the shower size, . The relative differences in shower size with respect to the reference provide estimates of the systematic uncertainties.
In figure 16 the relative differences for several scenarios are shown as a function of the shower zenith angle. The largest uncertainty in corresponds to the case of the correction computed with 100 EeV proton showers simulated with Sibyll. It decreases from 12% at to less than 3% in absolute value at . This is a conservative estimate since it corresponds to an extreme situation in which the majority of data having energies of a few EeV are reconstructed with the electromagnetic correction corresponding to a 100 EeV shower.
There can be other possible systematic effects on the shower size which depend on zenith angle and are not absorbed in the calibration procedure. Different hadronic interaction models and primary compositions may have different muon attenuations that is manifested in the zenith angle dependence of the muon distributions. These have an intrinsic dependence on the zenith angle which could in principle differ from that of the data. In addition, there are implicit uncertainties which can have zenith angle, dependence due to the accuracy of the models used for the muon distributions, the detector response and, most importantly, the electromagnetic correction.
The extent of these systematic uncertainties can be tested by exploring the zenith angle distribution of events above given thresholds of shower size. If the arrival directions are isotropically distributed and the detector has an efficiency that is independent of zenith angle , the events should have a flat distribution in . Possible systematic effects associated with the zenith angle will appear as deviations from a uniform distribution. This test is quite sensitive due to the steeply falling spectrum. For a spectral index it can be shown that a relative systematic shift of shower size in a given zenith angle bin should lead to a relative increase in the corresponding bin of .
Figure 17 displays the distribution of events with in bins of equal geometrical exposure. The array is fully efficient for showers of this size. The plot indicates that systematic deviations are within a 15% band with a standard deviation of . Assuming a spectral index as obtained with events below , and for energies above 4 EeV [33] this systematic uncertainty corresponds to a potential shift in of less than , well within the combined estimated uncertainties of the muon distributions, detector response and electromagnetic corrections.
5 Energy calibration and resolution
It is possible to relate to the shower energy using simulations, but due to the lack of knowledge of composition and hadronic models the corresponding systematic uncertainties become quite large (see section 3.1). Alternatively, the correlation between and shower energy reconstructed with the FD can be obtained from a subset of “golden hybrid” events (events for which both FD and SD reconstructions are possible), similarly to what has been used to calibrate events below [33]. The energy scale inferred from this subset is applied to all the inclined showers recorded with the SD.
The golden hybrid events with zenith angle greater than are required to pass the inclined T4 and T5 conditions and, in addition, to satisfy a set of FD quality cuts specifically designed to ensure an accurate reconstruction of the arrival direction and of the longitudinal profile. The cuts are adapted versions of those used in the calibration of events with [43]. The station closest to the shower core which is used for the geometrical reconstruction must be at a distance less than 750 m. For a precise estimate of the energy, we require an adequate monitoring of the atmospheric conditions (vertical aerosol optical depth up to 0.1; cloud coverage less than 25% in the FD field of view, distance of the cloud layer to the measured profile greater than 50 g/cm, and thickness of the cloud layer less than 100 g/cm). Furthermore, we exclude a residual contamination of shower profiles distorted by clouds and aerosols by requiring a GaisserHillas fit with a residual smaller than 3 and a negative value of the parameter of the fitted GaisserHillas profile.^{4}^{4}4In the GaisserHillas function, is a parameter not to be confused with the depth of the first interaction. Showers in both data and simulation are found to be best described by negative values [44, 45]. Moreover, the maximum accepted uncertainty of is 150 g/cm. In addition to the quality selection criteria, a fiducial cut on the FD field of view (FOV) is applied [46], ensuring that it is large enough to observe all plausible values of the shower maximum . This “fiducial FOV cut” includes a restriction on the minimum viewing angle of the light in the FD telescope (). Finally, only events with FD energies greater than eV are accepted to ensure a trigger probability of nearly 100% for the SD and FD detectors. For the period from 1 January 2004 to 31 December 2012 the sample has 223 hybrid events with .
To describe the correlation it is sufficient to perform a powerlaw fit to the shower size as a function of the calorimetric hybrid energy ,
(5.1) 
which is then inverted to give the energy conversion. The slope parameter is related to () as discussed in section 3.1. The fit must be handled with care since it is performed on a subset of the data used to calculate the spectrum. A bias can be expected at low energies due to the threshold trigger effects of the SD stations.
The fit is based on a tailored maximumlikelihood method [47] that takes into account the effect of the cut on energies greater than eV, which would otherwise distort a standard leastsquares fit. The ability to include both the uncertainties of the reconstructions of (see figure 13 in section 4.2.1) and , without relying on approximations, is another advantage of this approach. We achieve this by regarding the data as drawn from a twodimensional probability density function that describes the random fluctuations away from the ideal curve given by eq. (5.1). The PDF is centered around this curve and its width is given by the uncertainties of the reconstructed values of and . Showertoshower fluctuations of also contribute to the spread and are taken into account. Their relative size is assumed to be constant in the energy range of interest and fitted to the data as an extra parameter. Since the data are described by a twodimensional PDF, we are able to model the effect of an energy cut by setting the probability to observe events below the cut to zero. The effect of event migration in and out of the accepted region is also treated correctly. The method was validated with extensive MonteCarlo studies and was found to fit the calibration curve without bias. The details will be presented in a dedicated article.
The result of this fit is illustrated in figure 18 and the best resulting parameters are and . The calibration accuracy at the highest energies is limited by the number of events (the most energetic is at eV). An alternative fitting method based on leastsquares minimization was used to calibrate the selected hybrid set as a crosscheck. In this method the trigger bias is removed with an elliptical cut in the showersizeenergy plane to ensure that trigger effects can be ignored. The parameters obtained from the best fit are compatible with the maximumlikelihood fit.
Eq. (5.1) can be used to convert to shower energy for all the inclined events recorded with the SD. The statistical uncertainties of the calibration constants and were converted to uncertainties in in the left panel of figure 19. The latter range between 1.5% at eV and 5% at eV.
To be detected with the FD the showers must have a depth of maximum well within the field of view, and to be detected with the SD the shower core must fall inside the array. The two simultaneous conditions suppress events with large zenith angles due to the geometrical layout of the Observatory [48]. The systematic uncertainty in , accounting for the different angular distributions of the golden hybrid events and the full inclined data set used to calculate the spectrum, is on average .
In addition, there is an overall systematic uncertainty of 14% from the FD energy measurement which relies on the knowledge of the fluorescence yield used in the energy determination of the FD events (3.6%), atmospheric conditions (3.4 to 6.2%), absolute detector calibration (9.9%), invisible energy (3 to 1.5%), stability of the energy scale (5%), and shower reconstruction (6.5 to 5.6%) [49].
The total systematic uncertainty in ranges between 14% at eV and 17% at eV, and is obtained by adding the uncertainties described above in quadrature.
The right panel of figure 19 illustrates the ratio of the inferred SD energy (after the calibration procedure) and the reconstructed FDHybrid energy . The resolution in the SD energy can be inferred from this ratio distribution [50], by fixing the FD energy resolution to 7.6% [49]. The resulting average SD energy resolution is for the selected golden hybrid set. This SD energy resolution is attributed to the combined effect of showertoshower fluctuations and the reconstruction uncertainty of the fitting procedure (see figure 13 in section 4.2.1).
The fact that the value of the calibration parameter is not exactly equal to one reflects the fact that the signals detected with the SD are larger than expected from the reference distribution using protons. This is similar to the situation for events at zenith angles less than , where the measured signal exceeds that expected for proton showers at the same energy, but the energy of the muons involved in inclined showers is significantly higher [3]. The obtained value of implies that corresponds to an FD energy of 5.75 EeV, about 42% less than the energy of 10 EeV used for the proton reference distributions. Equivalently, it implies that the number of muons in a 10 EeV shower is larger than the average muon content of the reference proton showers simulated with QGSJetII03 [51]. This result has important implications concerning composition and hadronic models, which will be addressed in a separate article.
6 Summary and conclusions
We have developed a procedure to reconstruct inclined showers detected with the SD of the Pierre Auger Observatory. After reconstructing the arrival direction in a standard way using the start times of the recorded signals in the triggered stations, the signals of each event are fitted to the expected shape of the muon distribution at the ground. This requires modeling of the corresponding muon patterns, the response of the detectors to the passage of particles, and the amount of electromagnetic component of inclined showers, all of which were obtained with the aid of comprehensive simulations. The reference patterns for the muon distribution for each arrival direction are based on proton simulations with QGSJetII.03 at a fixed energy of 10 EeV. A size parameter has been introduced to scale the SD measurements to the reference distributions. is used as an energy estimator, calibrating it with the calorimetric energy measured by the FD, using a subsample of quality hybrid events which can be reconstructed independently using both the FD and SD techniques. In this process the simulated proton showers used to obtain the reference distributions of muons at the ground are shown to have a large deficit of muons with respect to the measured data.
The performance of the reconstruction was studied for events in the zenith angle range between and , and with exceeding 0.7, to ensure nearly full efficiency for the trigger, selection and reconstruction of events that fall in active regions of the array. The angular resolution was shown to be better than . The SD energy resolution was obtained, comparing the reconstructed SD and FD energies for the subset of events used in the energy calibration. The average SD energy resolution was shown to be about 19.3% attributed both to the reconstruction procedure and to the showertoshower fluctuations. Systematic uncertainties due to the calibration procedure were estimated to be energy dependent, and to fall in the range from 1.5 to 5%. In addition, there is an overall systematic uncertainty of 14% associated with the FD energy assignment. The procedure to analyze inclined events has a resolution in size and angular accuracy which is comparable to that obtained for events with . The analysis of these events opens the possibility to measure the cosmic ray spectrum in an independent way, to explore primary composition, to study the fidelity of current hadronic models extrapolated to the high energies relevant to these events, and to analyze arrival directions from parts of the sky that are not accessible using events with zenith angles less than . An independent measurement of the energy spectrum of ultrahigh energy cosmic rays using very inclined showers is to be presented in a forthcoming publication.
Acknowledgments
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: Comisión Nacional de Energía Atómica, Fundación Antorchas, 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, Argentina; the Australian Research Council; 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 # 2010/073596, # 1999/054043, Ministério de Ciência e Tecnologia (MCT), Brazil; MSMTCR LG13007, 7AMB14AR005, CZ.1.05/2.1.00/03.0058 and the Czech Science Foundation grant 1417501S, Czech Republic; Centre de Calcul IN2P3/CNRS, Centre National de la Recherche Scientifique (CNRS), Conseil Régional IledeFrance, Département Physique Nucléaire et Corpusculaire (PNCIN2P3/CNRS), Département Sciences de l’Univers (SDUINSU/CNRS), France; Bundesministerium für Bildung und Forschung (BMBF), Deutsche Forschungsgemeinschaft (DFG), Finanzministerium BadenWürttemberg, HelmholtzGemeinschaft Deutscher Forschungszentren (HGF), Ministerium für Wissenschaft und Forschung, Nordrhein Westfalen, Ministerium für Wissenschaft, Forschung und Kunst, BadenWürttemberg, Germany; Istituto Nazionale di Fisica Nucleare (INFN), Ministero dell’Istruzione, dell’Università e della Ricerca (MIUR), Gran Sasso Center for Astroparticle Physics (CFA), CETEMPS Center of Excellence, Italy; Consejo Nacional de Ciencia y Tecnología (CONACYT), Mexico; Ministerie van Onderwijs, Cultuur en Wetenschap, Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO), Stichting voor Fundamenteel Onderzoek der Materie (FOM), Netherlands; National Centre for Research and Development, Grant Nos.ERANETASPERA/01/11 and ERANETASPERA/02/11, National Science Centre, Grant Nos. 2013/08/M/ST9/00322 and 2013/08/M/ST9/00728, Poland; Portuguese national funds and FEDER funds within COMPETE  Programa Operacional Factores de Competitividade through Fundação para a Ciência e a Tecnologia, Portugal; Romanian Authority for Scientific Research ANCS, CNDIUEFISCDI partnership projects nr.20/2012 and nr.194/2012, project nr.1/ASPERA2/2012 ERANET, PNIIRUPD20113014517, and PNIIRUPD201130062, the Minister of National Education, Programme for research  Space Technology and Advanced Research  STAR, project number 83/2013, Romania; Slovenian Research Agency, Slovenia; Comunidad de Madrid, FEDER funds, Ministerio de Educación y Ciencia, Xunta de Galicia, Spain; The Leverhulme Foundation, Science and Technology Facilities Council, United Kingdom; Department of Energy, Contract No. DEAC0207CH11359, DEFR0204ER41300, and DEFG0299ER41107, National Science Foundation, Grant No. 0450696, The Grainger Foundation, USA; NAFOSTED, Vietnam; Marie CurieIRSES/EPLANET, European Particle Physics Latin American Network, European Union 7th Framework Program, Grant No. PIRSES2009GA246806; and UNESCO.
References
 [1] Pierre Auger collaboration, J. Abraham et al., Properties and performance of the prototype instrument for the Pierre Auger Observatory, Nucl. Instrum. Meth. A 523 (2004) 50.
 [2] A. M. Hillas, in Proceedings of the 11 ICRC, Budapest, Hungry, Acta Phys. Acad. Sci. Hung. 29 (1970) 533.
 [3] M. Ave, R. A. Vazquez, and E. Zas, Modelling horizontal air showers induced by cosmic rays, Astropart. Phys. 14 (2000) 91 [astroph/0011490].
 [4] M. Ave, R. A. Vazquez, E. Zas, J. A. Hinton, and A. A. Watson, The rate of cosmic ray showers at large zenith angles: A step towards the detection of ultrahigh energy neutrinos by the Pierre Auger Observatory, Astropart. Phys. 14 (2000) 109 [astroph/0003011].
 [5] V. M. Olmos Gilbaja, A Complete method to obtain the energy spectrum of inclined cosmic rays detected with the Pierre Auger Observatory (PhD thesis), University of Santiago de Compostela, Spain (2009).
 [6] H. Dembinski, Measurement of the ultrahigh energy cosmic ray flux using data from very inclined air showers at the Pierre Auger Observatory (PhD thesis), RWTH Aachen University, Germany (2009).
 [7] R. A. Vazquez, for the Pierre Auger collaboration, The cosmic ray flux observed at zenith angles larger than with the Pierre Auger Observatory, in Proceedings of the 31 ICRC, Łódź, Poland (2009) [arXiv:0906.2189].
 [8] H. Dembinski, for the Pierre Auger collaboration, The cosmic ray spectrum above eV as measured with inclined showers recorded at the Pierre Auger Observatory, in Proceedings of the 32 ICRC, Beijing, China, 2 (2011) 101 [arXiv:1107.4809].
 [9] A. Schulz, for the Pierre Auger collaboration, The measurement of the energy spectrum of cosmic rays above eV with the Pierre Auger Observatory, in Proceedings of the 33 ICRC, Rio de Janeiro, Brasil (2013) [arXiv:1307.5059].
 [10] M. Ave, J. A. Hinton, R. A. Vazquez, A. A. Watson and E. Zas, New constraints from Haverah Park data on the photon and iron fluxes of UHE cosmic rays, Phys. Rev. Lett. 85 (2000) 2244 [astroph/0007386].
 [11] M. Ave, J. A. Hinton, R. A. Vazquez, A. A. Watson, and E. Zas, Sensitivity of the Auger Observatory to ultra high energy photon composition through inclined showers, Phys. Rev. D 67 (2003) 043005 [astroph/0208228].
 [12] V. S. Berezinsky and A. Y. .Smirnov, Cosmic neutrinos of ultrahigh energies and detection possibility, Astrophys. Space Sci. 32 (1975) 461.
 [13] K. S. Capelle, J. W. Cronin, G. Parente and E. Zas, On the detection of ultrahighenergy neutrinos with the Auger Observatory, Astropart. Phys. 8 (1998) 321 [astroph/9801313].
 [14] Pierre Auger collaboration, X. Bertou et al., Calibration of the surface array of the Pierre Auger Observatory, Nucl. Instrum. Meth. A 568 (2006) 839.
 [15] I. Valino, J. AlvarezMuniz, M. Roth, R. A. Vazquez, and E. Zas, Characterisation of the electromagnetic component in ultrahigh energy inclined air showers, Astropart. Phys. 32 (2010) 304 [arXiv:0910.2873].
 [16] Pierre Auger collaboration, J. Abraham et al., The Fluorescence Detector of the Pierre Auger Observatory, Nucl. Instrum. Meth. A 620 (2010) 227 [arXiv:0907.4282].
 [17] M. Unger, B. R. Dawson, R. Engel, F. Schussler, and R. Ulrich, Reconstruction of Longitudinal Profiles of UltraHigh Energy Cosmic Ray Showers from Fluorescence and Cherenkov Light Measurements, Nucl. Instrum. Meth. A 588 (2008) [arXiv:0801.4309].
 [18] Pierre Auger collaboration, J. Abraham et al., A Study of the Effect of Molecular and Aerosol Conditions in the Atmosphere on Air Fluorescence Measurements at the Pierre Auger Observatory, Astropart. Phys. 33 (2010) 108 [arXiv:1002.0366].
 [19] T. K.Gaisser and A. M. Hillas, Reliability of the method of constant intensity cuts for reconstructing the average development of vertical showers in Proceedings of the 15 ICRC, Plotdiv, Bulgaria, 8 (1977) 353.
 [20] Pierre Auger collaboration, J. Abraham et al., Trigger and aperture of the surface detector array of the Pierre Auger Observatory, Nucl. Instrum. Meth. A 613 (2010) 29 [arXiv:1111.6764].
 [21] H. P. Dembinski, P. Billoir, O. Deligny, and T. Hebbeker, A phenomenological model of the muon density profile on the ground of very inclined air showers, Astropart. Phys. 34 (2010) 128 [arXiv:0904.2372].
 [22] L. Cazon, R. A. Vazquez, A. A. Watson, and E. Zas, Time structure of muonic showers, Astropart. Phys. 21 (2004) 71 [astroph/0311223].
 [23] M. Ave, J. A. Hinton, R. A. Vazquez, A. A. Watson and E. Zas, Constraints on the ultrahighenergy photon flux using inclined showers from the Haverah Park array, Phys. Rev. D 65 (2002) 063007 [astroph/0110613].
 [24] S. Ostapchenko, Nonlinear screening effects in high energy hadronic interactions, Phys. Rev. D 74 (2006) 014026 [hepph/0505259].
 [25] S. J. Sciutto, The AIRES system for air shower simulations: An Update, in Proceedings of the 27 ICRC, Hamburg, Germany, 1 (2001) 237, (Copernicus Gesellschaft, Hamburg, 2001) [astroph/0106044].
 [26] N. N. Kalmykov and S. S. Ostapchenko, The Nucleusnucleus interaction, nuclear fragmentation, and fluctuations of extensive air showers, Phys. Atom. Nucl. 56 (1993) 346 [Yad. Fiz. 56N3 (1993) 105]; N. N. Kalmykov, S. S. Ostapchenko and A. I. Pavlov, EAS and a quark  gluon string model with jets, Bull. Russ. Acad. Sci. Phys. 58 (1994) 1966 [Izv. Ross. Akad. Nauk Ser. Fiz. 58N12 (1994) 21].
 [27] D. Heck, G. Schatz, T. Thouw, J. Knapp and J. N. Capdevielle, CORSIKA: A Monte Carlo code to simulate extensive air showers, Report FZKA 6019, Forschungszentrum Karlsruhe, 1998.
 [28] S. Ostapchenko, Monte Carlo treatment of hadronic interactions in enhanced Pomeron scheme: I. QGSJETII model, Phys. Rev. D 83 (2011) 014018 [arXiv:1010.1869].
 [29] T. Pierog and K. Werner, EPOS Model and Ultra High Energy Cosmic Rays, Nucl. Phys. Proc. Suppl. 196 (2009) 102–105 [arXiv:0905.1198].
 [30] T. Pierog, I. .Karpenko, J. M. Katzy, E. Yatsenko and K. Werner, EPOS LHC : test of collective hadronization with LHC data [arXiv:1306.0121].
 [31] R. Engel, T. K. Gaisser, P. Lipari, and T. Stanev, Nucleusnucleus collisions and interpretation of cosmicray cascades Phys. Rev. D 46 (1992) 5013; R. S. Fletcher, T. K. Gaisser, P. Lipari, and T. Stanev, SIBILL: An event generator for simulation of high energy cosmic ray cascades, Phys. Rev. D 50 (1994) 5710.
 [32] M. Ave, J. A. Hinton, R. A. Vazquez, A. A. Watson and E. Zas, A New approach to inferring the mass composition of cosmic rays at energies above 10**18 eV, Astropart. Phys. 18 (2003) 367 [astroph/0112071].
 [33] Pierre Auger collaboration, J. Abraham et al., Observation of the suppression of the flux of cosmic rays above eV, Phys. Rev. Lett. 101 (2008) 061101 [arXiv:0806.4302].
 [34] S. Argiro, S. L. C. Barroso, J. Gonzalez, L. Nellen, T. C. Paul, T. A. Porter, L. Prado, Jr. and M. Roth et al., The Offline Software Framework of the Pierre Auger Observatory, Nucl. Instrum. Meth. A 580 (2007) 1485 [arXiv:0707.1652].
 [35] GEANT4 collaboration, S. Agostinelli et al., GEANT4: A Simulation toolkit, Nucl. Instrum. Meth. A 506 (2003) 250; J. Allison et al., Geant4 developments and applications, IEEE Trans. Nucl. Sci. 53 (2006) 270.
 [36] G. Rodríguez Fernández, Horizontal air showers at the Pierre Auger Observatory (PhD thesis), University of Santiago de Compostela, Spain (2007).
 [37] A. M. Hillas, in Proceedings of the 17 ICRC, Paris, France, 8 (1981) 193.
 [38] Pierre Auger collaboration, C. Bonifazi et al., A model for the time uncertainty measurements in the Auger surface detector array, Astropart. Phys. 28 (2008) 523 [arXiv:0705.1856].
 [39] Pierre Auger collaboration, C. Bonifazi, The angular resolution of the Pierre Auger Observatory, Nucl. Phys. Proc. Suppl. 190 (2009) 20 [arXiv:0901.3138].
 [40] I. Valiño, for the Pierre Auger collaboration, A measurement of the muon number in showers using inclined events recorded at the Pierre Auger Observatory, in Proceedings of the 33 ICRC, Rio de Janeiro, Brasil (2013) [arXiv:1307.5059].
 [41] Pierre Auger collaboration, A measurement of the muon number in showers using inclined events detected at the Pierre Auger Observatory, in preparation.
 [42] J. AlvarezMuniz, R. Engel, T. K. Gaisser, J. A. Ortiz and T. Stanev, Atmospheric shower fluctuations and the constant intensity cut method, Phys. Rev. D 66 (2002) 123004 [astroph/0209117].
 [43] R. Pesce, for the Pierre Auger collaboration, Update of the energy calibration of data recorded with the surface detector of the Pierre Auger Observatory, in Proceedings of the 32 ICRC, Beijing, China, 2 (2011) 214 [arXiv:1107.4809].
 [44] C. Song, Longitudinal profile of extensive air showers, Astropart. Phys. 22 (2004) 151.
 [45] HIRES collaboration, T. AbuZayyad et al., A Measurement of the average longitudinal development profile of cosmic ray air showers between 10**17eV and 10**18eV, Astropart. Phys. 16 (2001) 1 [astroph/0008206].
 [46] Pierre Auger collaboration, J. Abraham et al., Measurement of the Depth of Maximum of Extensive Air Showers above 10 eV, Phys. Rev. Lett. 104 (2010) 091101 [arXiv:1002.0699].
 [47] H. Dembinski, for the Pierre Auger collaboration, The cosmic ray spectrum above 4 EeV as measured with inclined showers recorded at the Pierre Auger Observatory, in Proceedings of the 32 ICRC, Beijing, China, 2 (2011) 101 [arXiv:1107.4804].
 [48] C. K. Guerard, for the Pierre Auger collaboration, Acceptance of the Pierre Auger Southern Observatory fluorescence detector to neutrinolike air showers, in Proceedings of the 27 ICRC, Hamburg, Germany, 1 (2001) 760, (Copernicus Gesellschaft, Hamburg, 2001).
 [49] V. Verzi, for the Pierre Auger collaboration, The Energy Scale of the Pierre Auger Observatory, in Proceedings of the 33 ICRC, Rio de Janeiro, Brasil (2013) [arXiv:1307.5059].
 [50] R. C. Geary, The frequency distribution of the quotient of two normal variates, J. R. Stat. Soc. 93 (1930) 442–446.
 [51] G. Rodriguez, for the Pierre Auger collaboration, Inclined showers at the Pierre Auger Observatory: reconstruction, energy calibration and implications for the muon content, in Proceedings of the 32 ICRC, Beijing, China, 2 (2011) 95 [arXiv:1107.4809].