A complete method to obtain the energy spectrum of inclined cosmic rays detected with the Pierre Auger Observatory

Universidade de Santiago de Compostela Departamento de Física de Partículas A complete method to obtain the energy spectrum of inclined cosmic rays detected with the Pierre Auger Observatory. Víctor Manuel Olmos Gilbaja Santiago de Compostela, marzo - 2009 Universidade de Santiago de Compostela Departamento de Física de Partículas A complete method to obtain the energy spectrum of inclined cosmic rays detected with the Pierre Auger Observatory. Memoria presentada para optar al Grado de Doctor en Física por D. Víctor Manuel Olmos Gilbaja Santiago de Compostela, marzo 2009 Fdo. Víctor Manuel Olmos Gilbaja D. Enrique Zas Arregui, catedrático de Física Teórica del Departamento de Física de Partículas de la Universidad de Santiago de Compostela, CERTIFICA que la memoria titulada A complete method to obtain the energy spectrum of inclined cosmic rays detected with the Pierre Auger Observatory ha sido realizada bajo su dirección por D. Víctor Manuel Olmos Gilbaja en el Departamento de Física de Partículas de la Universidad de Santiago de Compostela y constituye el trabajo de tesis que presenta para optar al Grado de Doctor en Física. Santiago de Compostela a 3 de marzo de 2009 Fdo. Enrique Zas Arregui A mi niña Agradecimientos Después de una revisión bibliográca entre los agradecimientos de otras tesis doctorales más extensa si cabe que la realizada para llevar a cabo este trabajo que constituye mi tesis doctoral, he encontrado aún más similitudes entre todas ellas en esta sección, que las existentes entre los distintos textos en los primeros capítulos. Se tiende a agrupar a los seres a los que se siente agradecido de un modo muy similar independientemente del país de origen del espécimen o del sexo del doctorando: pareja, familia, amigos, director de tesis, compañeros de trabajo y agencia nanciadora. Así mismo, el motivo del agradecimiento a cada uno de estos colectivos se parece de una manera cuando menos sorprendente. En mi caso, he decido no desentonar demasiado con las referencias bibliográcas consultadas en materia de agradecimientos. Este mismo principio es el que he seguido en el capítulo dedicado al maravilloso mundo del rayo cósmico y al experimento con el cual unos cuantos cientícos han decidido aprender algo más sobre ellos, muy a pesar, de que los directores de tesis se empeñen en que innovemos en este capítulo presente en TODAS las tesis sobre temática similar. Antes de comenzar, me gustaría hacer una aclaración porque la gente tiende a dar mucha importancia al orden en el que aparecen mencionados, por esto quiero dejar claro, que el orden de aparición no tiene nada que ver con la importancia que cada uno de ellos tiene en mi vida. Dicho esto vamos al toro no sin pedir disculpas a la gente de quien me olvide y a aquellos que no sean agradecidos por todo aquello que han hecho por mí. En primer lugar a mi niña. No existe manera de describir todos los agradecimientos que te mereces, ni lo mucho que te quiero. Por aguantarme, por estar a mi lado, por reírte conmigo y de mí, por ser como eres, por salir y por quedarte en casa, por animarme día a día, mes a mes, por apoyarme siempre y en todo momento en mis decisiones, por los cuatro meses y por el n de semana de después, por tu habilidad natural con la electrónica, por todo ello y por mucho más. 16 CONTENTS List of Figures 2.1 Cosmic ray energy spectrum as measured by dierent experiments. . . . 29 2.2 Hillas plot of possible accelerating sources . . . . . . . . . . . . . . . . 30 2.3 Propagation features of Ultra High Energy Cosmic Rays . . . . . . . . . 31 2.4 Heitler toy model for cascade development . . . . . . . . . . . . . . . . 34 2.5 Longitudinal development of muonic and electromagnetic components of a 1019 eV shower ............................. 37 2.6 Layout of the hybrid detector . . . . . . . . . . . . . . . . . . . . . . . 39 2.7 Hybrid-stereo event detection . . . . . . . . . . . . . . . . . . . . . . . 39 2.8 Picture and schema of a uorescence detector building . . . . . . . . . . 40 2.9 Picture and schema of a uorescence telescope . . . . . . . . . . . . . . 41 2.10 Fluorescence camera and mercedes ring . . . . . . . . . . . . . . . . . . 41 2.11 Deep water Cherenkov tank . . . . . . . . . . . . . . . . . . . . . . . . 43 2.12 VEM calibration from atmospheric muons . . . . . . . . . . . . . . . . . 45 2.13 Trigger chain diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.14 Cosmic ray ux as measured by the Pierre Auger Observatory . . . . . . 50 2.15 Correlation of arrival direction of highest energy cosmic rays with AGN . 51 2.16 Elongation rate as measured by the uorescence detector of the Pierre AugerObservatory............................. 53 2.17 Photon limit as measured by the Pierre Auger Observatory . . . . . . . . 54 2.18 Neutrino limit as measured by the Pierre Auger Observatory . . . . . . . 55 3.1 Isolated muon in FADC trace . . . . . . . . . . . . . . . . . . . . . . . 67 3.2 Shower front propagation . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.3 Concentric shower front sketch . . . . . . . . . . . . . . . . . . . . . . 71 3.4 Arrival time distribution . . . . . . . . . . . . . . . . . . . . . . . . . . 77 18 LIST OF FIGURES 3.5 Muon delay from time model . . . . . . . . . . . . . . . . . . . . . . . 78 3.6 Picture of a doublet of the surface detector array . . . . . . . . . . . . . 79 3.7 Zenith angle distribution of doublet events . . . . . . . . . . . . . . . . 80 3.8 Signal, distance to core and start time dierence distribution for multiplet events ................................... 81 3.9 Start time variance for three models . . . . . . . . . . . . . . . . . . . . 83 3.10 Start time variance for LCB model . . . . . . . . . . . . . . . . . . . . 84 3.11 Start time variance for vertical models . . . . . . . . . . . . . . . . . . 85 3.12 Start time variance for vertical models . . . . . . . . . . . . . . . . . . 86 3.13 Picture of the Central Laser Facility . . . . . . . . . . . . . . . . . . . . 87 3.14 Comparison of the hybrid zenith angle with surface only reconstructions . 88 3.15 Space angle between hybrid and surface only reconstructions . . . . . . . 89 3.16 Study of some of the space angle systematics . . . . . . . . . . . . . . . 90 4.1 χ2 map for a high multiplicity event . . . . . . . . . . . . . . . . . . . . 97 4.2 Distance to the barycentre of reconstructed core . . . . . . . . . . . . . 99 4.3 χ2 map with and without zero signal stations . . . . . . . . . . . . . . . 100 4.4 Likelihood function maximization . . . . . . . . . . . . . . . . . . . . . 104 4.5 Likelihood function maximization . . . . . . . . . . . . . . . . . . . . . 105 4.6 Inputs for muon maps generation . . . . . . . . . . . . . . . . . . . . . 107 4.7 Muonmaps ................................108 4.8 Signal probability distribution . . . . . . . . . . . . . . . . . . . . . . . 111 4.9 EMratio..................................114 4.10 Attenuation of the total number of muons from maps . . . . . . . . . . 116 4.11Correctionfactor..............................117 4.12 Attenuation curve for inclined events . . . . . . . . . . . . . . . . . . . 118 4.13 Correlation between log10(E) and log10(N19) for the selected 286 events.121 4.14 N19 uncertainty from the maximum likelihood procedure . . . . . . . . . 122 4.15 N19 uncertainty from zenith angle uncertainty . . . . . . . . . . . . . . 124 4.16Ellipticalcut................................125 4.17 Relative uncertainties on EFD and on N19 ................126 4.18 Correlation between log10(E) and log10(N19) for the selected 117 events used to get the calibration line . . . . . . . . . . . . . . . . . . . . . . . 127 LIST OF FIGURES 19 4.19 ESD resolution...............................128 4.20Orthogonalpull ..............................128 4.21 Calibration by FD eye . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 4.22 Calibration for dierent low energy cuts . . . . . . . . . . . . . . . . . . 131 5.1 T4eciencycurve.............................139 5.2 Unit cell of the surface detector array . . . . . . . . . . . . . . . . . . . 140 5.3 Dierential ux of cosmic rays as a function of energy . . . . . . . . . . 142 5.4 Energy spectrum multiplied by E3 .....................143 5.5 Relative uncertainty due to the calibration process . . . . . . . . . . . . 145 5.6 N19 uncertainty due to muon maps . . . . . . . . . . . . . . . . . . . . 146 5.7 Systematic uncertainty in energy due to muon maps . . . . . . . . . . . 147 5.8 Systematic uncertainty in N19 due to the parameterization of the electromagneticsignal.............................148 20 LIST OF FIGURES Chapter 1 Introduction In 1912, Victor F. Hess found that the rate of discharge of the ionization chambers that he ew on several balloon ights increased with altitude. This observation was interpreted as an evidence that ionizing radiation was continuously reaching the Earth from outer space. From the charge spectrum of this radiation it was shown that nuclei from hydrogen to iron were present in cosmic rays. Moreover, new elementary particles such as positrons, muons or pions were discovered in cosmic rays. Twenty six years after Hess discovery, P. Auger and collaborators showed that particle cascades of high energy caused the correlated detections of particles at counters widely separated at ground. As low uxes demand large collection areas, extensive air showers discovered by Auger (and suspected by Rossi two years earlier) allow the detection of Ultra High Energy Cosmic Rays (UHECRs) using the atmosphere as target medium as well as a calorimeter. To characterize the properties of UHECRs and to enlighten the knowledge about their origin, their chemical composition and their energy spectrum, the Pierre Auger Observatory was born at the beginning of the 1990's. It is a hybrid detector that combines a surface detector array of 1600 Cherenkov water tanks with a uorescence detector of 24 telescopes to measure with unprecedented accuracy and statistics and with a full sky coverage the extensive air showers initiated by cosmic rays with energies above 1018 eV . The southern site of this observatory is already built in the Argentinian province of Mendoza and has been taking high quality data for more than 5 years. The northern site will be built in the American state of Colorado. The Pierre Auger Observatory 22 1 Introduction is well suited for the detection of inclined showers because the water Cherenkov tanks used in the surface detector act like volume detectors. The transverse area of 4.3m2 for completely horizontal showers is enough to detect particle uxes at relatively large distances from the shower axis. Besides, the fact that the number of stations within a given distance to shower axis rapidly increases as the zenith angle rises above 60◦ , improves the capacity of detecting inclined showers. The present work is devoted to the analysis of inclined air showers detected by the surface detector of the Pierre Auger Observatory and to determine the ux of UHECRs. It is presented as a complete guide of the reconstruction procedure of inclined events, so my contributions are sometimes merged with the contributions of other people to give a global view of the whole procedure. For this reason I will briey point out how my main contributions to obtaining the energy spectrum t in the general overview of this project which has taken over ten years. The rst analysis of inclined showers was done for the Haverah Park data in 2000 by M. Ave et al. [1]. A FORTRAN code was developed to analyze inclined showers being called et . Some predictions about the Auger Observatory were already done at that time. A new C++ code was developed within the USC astroparticle physics group coordinated by R. Váquez for the analysis of inclined air showers detected by the surface detector of the Pierre Auger Observatory. This code was based on the older FORTRAN version but it includes all the works developed within the group during those years. For instance, the time model for the propagation of muons of the shower done in L. Cazón PhD thesis [2], a new set of the proles of the muon density at ground and the simulation of the response to crossing muons of the Auger surface detector tanks done by G. Rodríguez in his PhD thesis [3] and the parameterizations of the electromagnetic signal in inclined showers derived by I. Valiño as a part of her doctoral work [4]. All these previous works established the basis of the next step towards the measurement of the energy spectrum of cosmic rays using inclined events, the energy calibration. This calibration between the shower size as measured by the surface detector and the calorimetric energy measurement performed by the uorescence detector has been one of my main contributions to the estimate of the cosmic ray ux using inclined events, being the rst energy spectrum of inclined cosmic rays based on this calibration presented in the collaboration meeting of November 2005 [5]. This has been followed up steadily since then and I have contributed signicanly to this task in which the group of USC 23 has played a prominent role. Besides the energy calibration, the exposure is a crucial ingredient to obtain the energy spectrum. I have detemined the energy saturation of the trigger eciency following one of the approaches already used in the analysis of vertical events using hybrid events with a zenith angle over 60◦ . Above this threshold, the exposure can be calculated from geometrical arguments and related to the work done by another analysis group of the Pierre Auger Observatory [6]. Apart from my contributions to the measurement of energy spectrum using inclined showers, other of my contributions to the analysis of inclined cosmic ray events are the checks of the model for the arrival time variance using twin tanks, the collaboration with the software development and bug xing of the et program or the development of the rst seed of the code ported to the oine which I did in collaboration with J. González and M. Roth during a stay within the Karlsruhe group that ended with a presentation of the rst inclined spectrum within the oine framework in the Chicago analysis meeting in 2006. The thesis is organized as follows: In chapter 2 a brief introduction to the cosmic rays and to the extensive air showers as well as a description of both detectors of the Pierre Auger Observatory are given. The latest published results of the Pierre Auger collaboration are briey discussed. In chapter 3 we describe the general procedure for the reconstruction of the incoming direction of the detected inclined showers. The time model that it is used to describe the shower front and to obtain the variance on the measurement of the start time is presented. Making use of the data collected in twin stations, the time variance is shown to describe data with a high accuracy. In chapter 4, the algorithms for the reconstruction of the impact point and of the energy are presented. Two methods are applied consecutively: a minimum χ2 procedure and a maximum likelihood method. The inputs for both of them, such as the tank response to muons or the electromagnetic contribution to the total signal are also discussed. In chapter 5, the energy spectrum of UHECRs is obtained. Firstly we obtain the acceptance of the surface detector to inclined events assuming that above the energy of trigger saturation the exposure of the array is completely geometrical. From the les containing the information on the trigger rate, its active area each second is computed and adding these areas for the whole period of data taking the exposure to inclined events is obtained. Once the exposure is obtained, the ux is directly calculated from the data. Finally in 24 1 Introduction chapter 6, main conclusions of this thesis are presented. Bibliography [1] M. Ave, R. A. Vázquez & E. Zas, Astropart. Phys. 14 (2000) 91. M. Ave, R. A. Vázquez, E. Zas, J. A. Hinton & A. A. Watson, Astropart. Phys. 14 (2000) 109. [2] L. Cazón PhD Thesis, Modelling the muon time distribution in extensive air showers, Universidad de Santiago de Compostela (2004). [3] G. Rodríguez-Fernández PhD Thesis, Horizontal air showers at the Pierre Auger Observatory, Universidad de Santiago de Compostela (2007). [4] I. Valiño PhD. Thesis, Detection of Horizontal Air Showers and Neutrino induced Showers with the Pierre Auger Observatory, Pierre Auger Collaboration Internal Document, GAP Note 2008-024. [5] N. Busca, V.M. Olmos-Gilbaja, P. Privitera, G. Rodriguez, R.A. Vazquez & E. Zas, An estimate of the cosmic ray spectrum using inclined data of the Pierre Auger Observatory, Pierre Auger Collaboration Internal Document, GAP Note 2006-025. [6] Acceptance working group. ipnweb.in2p3.fr/ auger/AugerProtected/AcceptWork.html 32 2 Cosmic rays features and the Auger observatory microwave background but also with the infrared background. The distance of production points of nuclei with energy above 1020 eV is reduced also to ≈100 Mpc [20]. The interactions of cosmic rays with background radiation during their propagation from their sources to the detection point at the Earth, have a direct eect in the energy spectrum. The observed ux should suer a great reduction beyond EGZK ≈5·1019 eV . This expected suppression is often referred to as the GZK cuto. The evidence of events above the limit imposed by the GZK eect proves that if it exists, the spectrum should not cut o sharply but it should present a smooth suppression. The GZK eect also implies a nearby origin for particles with energy above the GZK cuto. If the distance to their sources is limited to several tens of Mpc, charged particles of energies above EGZK should traverse cosmic magnetic elds with little deviation and in consequence they should point back to their sources. Such cosmic ray astronomy would make possible to identify the sources with known astrophysical objects or to establish the existence of new sources which are not visible at lower energies. 2.1.2 Extensive air showers In 1934, Rossi noticed that the coincidences between several counters exceeded the expected rate from pure chance coincidences [21]. In 1938, Pierre Auger and collaborators, after a systematic investigation, discovered extensive air showers (EAS) [22] which are huge cascades initiated by cosmic rays when they enter the Earth's atmosphere. Due to the very low ux of high energy cosmic rays, EAS are the only way they can be observed because using EAS, the eective area of detection is orders of magnitude larger than that of detectors that can be own on balloons or satellites. When the primary cosmic ray enters the Earth's atmosphere, it interacts with a nucleus of the air typically in the rst ≈100 g cm−2 from the top of the atmosphere. Part of the primary energy is used in the production of secondary particles which acquire transverse momenta traveling with an angle with respect to the direction of the primary particle. If the decay mean free path is longer than the interaction mean free path, the interaction process dominates creating more and more particles of lower energies. Decay can also contribute to the particle creation process. These particles also interact or decay as the cascade propagates in the atmosphere. For a hadron primary, secondary particles are mainly mesons that may decay into 2.1 Cosmic rays and extensive air showers 33 muons, electrons and photons. It is convenient to group the particles of the air shower into three categories often referred to as components : electromagnetic, muonic and hadronic. Neutral pions created during the shower development, almost immediately decay into two photons that initiate an electromagnetic subshower. These cascades from π0 decays continously feed the electromagnetic component of the shower. Charged mesons produce muons and neutrinos when they decay. Muons form the muonic component while neutrinos go through undetected carrying away a fraction of the primary energy. Finally, the hadrons that have not decayed form the hadronic component of the air shower. Electromagnetic and hadronic components increase in number of particles reaching a maximum after which the shower size attenuates. The muonic component does not suer a great attenuation after its maximum as muons lose energy mainly by ionization and not too many are lost by decay. The general features of the air showers development can be easily reproduced with the toy model proposed by Heitler in [23] for pure electromagnetic cascades. In the schematic diagram of g. 2.4 each segment represents a particle that at each vertex shares half of its energy with a new particle. Splitting occurs after an interaction length λ . After n interaction lengths the total number of particles is N(x) = 2x/λ being the energy per particle E(x) = E0/N(x) where E0 is the energy of the rst particle. The creation of new particles ends when the average particle energy is below the critical energy Ec . After this point, particles are ignored as they are more likely to lose energy, be absorbed or decay. The maximum number of particles within this model is N(xmax) = E0/Ec being reached at Xmax =λln(E0/Ec) ln2 . The basic features of this model which hold for electromagnetic cascades and also approximately for hadronic showers are: Nmax ∝E0and Xmax ∝ln(E0) Electromagnetic cascades The theory of electromagnetic cascades was developed at the end of the thirties after the discovery of the positron in 1932 [24]. The processes that dominate electromagnetic showers are pair production and bremsstrahlung and their behavior can be described quite accurately in terms of quantum electrodynamics. For particles with energies greater than 34 2 Cosmic rays features and the Auger observatory Figure 2.4: Heitler model for the development of electromagnetic cascades. After an interaction length λ each particle (straight line) creates a new one sharing half of its energy with it. their critical energy in air, collision losses and Compton scattering can be neglected. High energy photons interact with an air nucleus creating an electron-positron pair. Electrons emit high energy photons by bremsstrahlung which produce again more electronpositron pairs. At each stage of the shower, the number of particles increases while the energy per particle decreases. With several approximations to solve the cascade equations, the total number of charged particles for an electromagnetic shower initiated by a photon of energy E0 as a function of depth t in units of radiation lengths is given by the Greisen formula [25]: Ne(E0, t) = 0.31 β1/2 0 et(1−3 2ln(s)) (2.3) where s≈3t t+2β0 is called the age of the shower, β0= ln(E0/0) and 0 is the critical energy. Shower maximum is reached at s= 1 , the number of particles grows for s < 1 and decreases for s > 1 . This formula is an approximation of the average shower behavior and individual showers behave dierently due to uctuations in the interaction point and in the shower development. The lateral spread of the shower particles depends on the transverse momentum transmitted to the secondary particles created by pair production or bremsstrahlung and 2.1 Cosmic rays and extensive air showers 35 particularly on the Coulomb scattering of electrons (multiple elastic scattering). An approximate solution for the lateral distribution function (LDF) of electrons at a depth t is given by the NKG formula [25,26]: ρe(r, t) = Ne(t)C(s) r2 1r r1s−21 + r r1s−9/2 (2.4) where Ne is the total number of electrons, C(s) is a normalization coecient, s is the age of the shower and r1 is the Molière length in units of radiation lengths from the theory of multiple scattering. The dierences of the electromagnetic component of a shower initiated by a hadron with a pure electromagnetic cascade are quite subtle. Gaisser and Hillas proposed in [27] a parameterization for the number of electrons at a depth X initiated by a hadron as: Ne(X) = Nmax X−X0 Xmax −X0Xmax−X0 λ exp Xmax −X λ (2.5) where Nmax is the number of particles at the shower maximum, Xmax is the depth of the shower maximum in g cm−2 , and X0 and λ are parameters that depend on the point of rst interaction and on the shower development. Regarding the lateral spread of the electromagnetic component of an extensive air shower initiated by a hadron, the NKG function stated in eq. 2.4 is a reasonable approximation to it. However, dierent exponents are given from experimental measurements of the LDF that can not be described by just one age parameter. The muonic component and inclined cascades The muonic component that accounts for ∼10% of the total number of particles in an extensive air shower is due to the decay of charged pions and kaons according to the following decay modes: π±→µ±+ν (2.6) K±→µ±+ν (2.7) →π±+π0 (2.8) 36 2 Cosmic rays features and the Auger observatory As muons travel practically in straight lines from their production points, the lateral spread of the muonic component is mainly due to the transverse momentum of the parent pions and kaons. This lateral spread has a dependence on the energy as the opening angles for mesons are smaller for high energy particles. There is no standard parameterization for the LDF of muons as there is the NKG formula for the electrons. In any case, one of the rst LDF for muons as a function of the depth t was also given by Greisen in [28]: ρµ(t, r) = Nµ(t)r rG−0.75 1 + r rG−2.5 (2.9) where rG= 320 m plays a similar role than the Molière radius. As this value is greater than the Molière unit for electrons ( ≈100 m ), muons spread to larger distances from the shower core than electrons, being the main part of the particle density at large distances to the core. Besides its importance at large distances from the shower core, the muonic component becomes crucial for inclined showers. The grammage that a shower must cross before hitting the ground at the southern site of the Pierre Auger Observatory varies from ≈880 g cm−2 for a completely vertical shower to ≈1760 g cm−2 for a 60◦ zenith shower. After the shower maximum the electromagnetic component is rapidly attenuated in the atmosphere as it is shown in g. 2.5 while the attenuation length for the muonic component is much larger. For horizontal showers, the muonic component dominates at ground level because of the absorption of the electromagnetic cascades from π0 decays [29]. However, there is still an electromagnetic component due to muon decay and hard muon interactions. This component follows the behavior of the muonic one as also shown in g. 2.5. For this reason, the muonic component and the muonic signal are crucial for the reconstruction of inclined showers with a surface detector as the main part of the signal is due to the muons at ground. 2.2 Detection techniques of extensive air showers Array of detectors 2.2 Detection techniques of extensive air showers 37 Figure 2.5: Longitudinal development of the muonic and electromagnetic components of an extensive air shower of 1019 eV . From ≈2500 g cm−2 the electromagnetic component is due to muon decay. The rst technique used to detect extensive air shower and the one they were discovered with, is to sample the secondary particles of the shower when they reach the ground with an array of detectors. When a shower triggers an array of detectors, the particle densities and the arrival times in all detectors are registered. The shower arrival angle is generally obtained by tting the arrival times of the shower particles to a shower front moving at the speed of light. The shower energy is obtained comparing to shower simulations of the signals produced in each detector. One of the biggest problems analyzing data from arrays of particle detectors is the need to use simulations. These require assumptions about the nature of the primary cosmic ray and also about the behavior of interactions at energies well beyond those being explored with current particle accelerator experiments. The latter are usually made by choosing a given model for hadronic interactions. These are large inherent uncertainties associated to the simulations because of these unknowns. On the other hand, their duty cycle of ∼100% is essential to gain the enough statistics needed to study the very low ux of high energy cosmic rays. Fluorescence telescopes The detection of the nitrogen uorescence was proposed to follow the longitudinal development of extensive air showers in [30]. The idea is to detect the isotropic nitrogen 38 2 Cosmic rays features and the Auger observatory uorescence light emitted by the deionization of these air molecules that get ionized at the passage of the air shower. Every electron of the shower produces in average about 4 uorescence photons per meter so only high energy showers with a huge amount of electrons produce enough light to generate a detectable signal. Emitted photons have wavelengths between 300 nm and 400 nm corresponding to energies of the transition levels of the nitrogen molecule. This detection technique also requires a clean atmosphere and a moonless nights for operation. The arrival direction of the showers is determined from the timing measurements and with the information of the pointing direction of the light detection system. The integral of the total light emitted along the longitudinal prole is a measurement of the energy of the shower though several corrections have to be made to take into account the attenuation of the light in the atmosphere due to absorption and scattering processes. 2.3 The Pierre Auger Observatory The Pierre Auger Observatory has been designed as a hybrid detector to take advantage of both techniques and to reduce systematic eects due to any of them. An example of the relative strengths are the high statistics provided by the array of detectors with its large detection area and its almost 100% duty cycle and the nearly calorimetric measurement of energy of the uorescence telescopes. Besides combining the strengths of each detection technique the hybrid approach provides some new features as the intercalibration of energy measurements with each technique, the enhancement of the sensitivity to composition and a better resolution in the determination of the arrival direction of the detected air showers. In the next sections, each of the detectors that constitute this hybrid observatory is briey described. To maximize the hybrid performance, the uorescence detector eyes are placed at the edges of the surface detector array as shown in g. 2.6, allowing the detection of the events by the two techniques simultaneously. All showers above E≈8·1018eV that fall inside the surface detector area when the uorescence detector is active are detected by the two techniques simultaneously. Events which are recorded and can be reconstructed with both techniques deserve a special consideration because they are the key for intercalibration of the energy measurements. They are referred to 2.3 The Pierre Auger Observatory 39 as Golden Hybrid Events . One of these events is shown in g. 2.7. Figure 2.6: Layout of the hybrid detector of the Pierre Auger Observatory. Red dots are the stations of the surface array. The four buildings of the uorescence detector are also shown. Figure 2.7: Golden hybrid event detected simultaneously with the uorescence and surface detector techniques. Seventeen SD tanks are triggered by the shower. The event is also stereo because it is detected with two eyes of the uorescence detector. 40 2 Cosmic rays features and the Auger observatory 2.3.1 The Fluorescence Detector Description of the Fluorescence Detector The Fluorescence Detector of the Pierre Auger Observatory is a set of four eyes each one housing six telescopes which are enclosed in a building. Each telescope has a eld of view (FoV) of 30◦ in azimuth x 28.6◦ in elevation between 2◦ and 30.6◦ above the horizon. One of the buildings is shown in g. 2.8. Figure 2.8: Picture (left) and schematic view (right) of Los Leones, one of the four buildings that constitute the uorescence detector. The six telescopes are also shown in the schematic view. The telescopes, as the one shown in g. 2.9, use adapted Schmidt optics to achieve a high optical quality in the whole FoV reducing the coma aberration [32]. Each telescope consists of an aperture system, a spherical mirror ( 3.4m curvature radius) and a camera with 440 photomultipliers (PMTs) placed at the focal surface [33]. Fluorescence light enters the telescope through a 0.85 m radius diaphragm equipped with a corrector ring to increase the eective radius to 1.1m [34]. After traversing an UV lter to remove much of the sky light background. The light is focused onto the camera by a 3.5m x 3.5m spherical mirror. The camera, described in [35] and shown in the left panel of g. 2.10, is an array of 440 hexagonal PMTs placed at the telescope focal surface. To adapt the camera to the spherical focal surface, a honeycomb conguration has been adopted arranging the PMTs in a 22x20 matrix. Each pixel of the camera has a FoV of about 1.5◦ . Due to the spacing between PMTs and to the ineciency of the region of contact between adjacent PMTs, a system of `Wiston cones [36] complements the PMTs at the camera to optimize the 2.3 The Pierre Auger Observatory 41 Figure 2.9: Picture (left) and schematic view (right) of one of the uorescence telescopes. light collection. The unit element of this system is the so-called mercedes star shown in the right panel of g. 2.10 which has been shown to increase the light collection eciency from 50% up to 93% [37]. Figure 2.10: PMT camera (left) and mercedes ring (right). Trigger levels The trigger system of the uorescence detector is a four level trigger described in [38,39]. The rst level trigger is located on the front-end board and makes decisions at the PMT 48 2 Cosmic rays features and the Auger observatory Besides this pattern search, an external third level trigger condition has been set up to use the hybrid performance of the Auger observatory as described in the FD detector section. Whenever a T3 is detected by the uorescence detector, all the stations are requested by the CDAS to send their data within a time window of 120 µs of the event time. This event time has to be corrected by the time it takes to the light to travel from the impact point to the detection eye. This time has been xed to the time it takes the light to travel from the center of the array (Celeste tank) to the corresponding eye. The SD event is agged as an FD T3 trigger. • Physical trigger (T4) This trigger is intended to select the events produced by extensive air showers from the bulk of detected showers than could have been produced by chance coincidences of local triggers. This trigger has a dierent implementation for vertical and for inclined showers. The physical trigger for vertical showers uses two main characteristics of vertical air showers. The rst one is a certain compactness of triggered tanks, the second one is the fact that most of the detected signals are spread enough in time to fulll the ToT trigger condition. The requirement of a compact 3 ToT (3 neighbor tanks in a triangular conguration) ensures that 99 % of the selected showers are physical events [52]. The requirements of the vertical T4 is too restrictive for inclined events because of their elongated proles and the fact that ToT triggers are typically due to the electromagnetic component of the shower which is practically absent in inclined showers. A new T4 is demanded to select inclined events out of the T3 data. The fourth level trigger is proposed in [53] as it was dened in [54]. It is based on the search of a plane shower front compatible in time with the maximum number of triggered stations. Three criteria are used to accept a conguration. A physical incoming direction has to be found, the residuals of the t have to be below a certain level and the stations are asked to lie within a cylinder around the shower axis taking into account that the number of stations is rather proportional to the area of the ground footprint 2.4 Recent results of the Pierre Auger Observatory 49 of the shower. Values of the residual tolerance and of the cylinder radius can be tuned to improve the selection criteria. • Quality trigger (T5) This last trigger level tries to ensure a good reconstruction of the parameters of the detected air showers as its incoming direction and its energy. This trigger excludes events that hit the ground too close to a hole or the edge of the array. Besides ensuring a high quality of the reconstruction parameters, the fth level trigger constrains the exposed area because it is very hard to compute the acceptance of the array if we take into account events that hit the ground at an impact point outside the array [55]. This trigger level ensures that a portion of the array at the area where the air shower hits the ground is properly working at the time it arrives. For vertical showers it is enough that the rst hexagon around the tank with the highest signal is active. For inclined showers, and due to their elongated footprints, the highest signal tank is not so representative of the impact point so the station closest to the reconstructed core is asked to be surrounded by a hexagon of active tanks. 2.4 Recent results of the Pierre Auger Observatory During the ve years that the south site of the Pierre Auger Observatory has been taking data, several results have been published. In this section a brief comment about ve of them concerning the cosmic ray ux, the study of their arrival directions and some hints about composition is given. 2.4.1 Energy spectrum The measurement of the energy spectrum of cosmic rays above 2.5·1018 eV , derived from 20,000 events recorded at the Pierre Auger Observatory was published in [56] and it is shown in g. 2.14. The spectral index γ of the ux J E−γ , at energies between 4·1018 eV and 4·1019 eV was quoted to be 2.69 ±0.02(stat)±0.06(syst) , steepening to 4.2±0.4(stat)±0.06(syst) at higher energies. 50 2 Cosmic rays features and the Auger observatory This measurement is consistent with the prediction by Greisen and by Zatsepin and Kuzmin [17]. With an exposure twice that of HiRes [57] and 4 times that of AGASA [58], our evidence supports the recent report of the former. Figure 2.14: Upper panel: The dierential ux J as a function of energy. Lower Panel: The fractional dierences between Auger (black dots) and HiRes I (grey squares) data compared with predictions of a spectrum with an index of γ=−2.69 . In both panels vertical lines are only statistical uncertainty. 2.4 Recent results of the Pierre Auger Observatory 51 2.4.2 Arrival directions A great challenge for astroparticle physics today is the identication of the sources of UHECR. If these cosmic rays are nuclei, only sources within ≈200 Mpc from the Earth can contribute to the cosmic ray ux. The ux from more distant sources is attenuated mainly due to interactions with the microwave background photons for protons [17] or to photo-desintegration processes for nuclei [59] as has been commented for the suppression of the cosmic ray ux. As nearby sources are not uniformly distributed and given the small magnetic deections of the trajectories of such high energetic cosmic rays, their arrival directions should be anisotropic pointing back to their origin. A recent result of the Pierre Auger Observatory shows a correlation of the arrival direction of very high energy cosmic rays with energies above 57 EeV with 3.2◦ windows at the directions of nearby ( z < 0.018 ) active galactic nuclei [60]. Active Galatic Nuclei (AGN) which are very active galaxies with supermassive black holes located in their centers have been postulated as sites for cosmic ray acceleration [61]. Fig. 2.15 shows the published arrival direction distribution together with the AGN locations from the Veron-Veron Cetty catalog with z < 0.018 [62]. Figure 2.15: Celestial sphere in galactic coordinates with the supergalactic plane indicated by the dashed line. Circles are 3.2◦ radius centered in the arrival direction of the highest energy cosmic rays detected by the observatory. Red stars indicate the position of 472 AGNs with redshift below 0.018 from the 12th catalog of quasar and active nuclei [62]. 52 2 Cosmic rays features and the Auger observatory This result is an evidence of the anisotropy of UHECR which has been tested against a given catalog. Further tests are also in progress. 2.4.3 Mass composition The cosmic ray composition has been studied by the Auger observatory using the longitudinal proles of the showers [63]. The depth of the shower maximum Xmax is directly measured by the uorescence detector of the observatory. The average Xmax has been compared with predictions from shower simulations with dierent hadronic models. The average of Xmax at a certain energy is related to the mean logarithmic mass as: < Xmax >=Dp[ln(E/E0)−< lnA >] + cp (2.10) where Dp is the elongation rate of a proton and cp is the average depth of proton with energy E0 . Dierent models give slightly dierent values for Dp and cp . With more than 4300 events after all the selection cuts, the resulting mean Xmax as a function of energy is shown in g. 2.16. A simple linear t to all the data gives an elongation rate of 54 ±2g cm−2/decade but the t is not very good. Allowing a break in the elongation rate, the ts shown in grey are much better giving elongation rates of 71 ±5g cm−2/decade below 1018.35 eV and 40 ±4g cm−2/decade above that energy cut. 2.4 Recent results of the Pierre Auger Observatory 53 Figure 2.16: Average shower maximum as a function of energy. The change in the elongation rate at 1018.35 eV may be due to a change in the composition of the cosmic rays. Numbers are the available events to nd the average elongation rate at that energy bin. Other parameters such as the uctuations of Xmax or the time structure of shower particles can be also used to infer mass composition [64]. 2.4.4 Photon limit Extensive air showers originated by photons have a dierent development in the atmosphere than those originated by hadrons. They develop over a larger atmospheric depth and the position of their shower maximum Xmax is deeper. As Xmax can be directly measured by the uorescence detector or extracted from measurements of the surface detector, limits to the photon fraction can be set by the Auger observatory as shown 54 2 Cosmic rays features and the Auger observatory in g. 2.17. Some models of UHECR production have been discarded with these limits due to the high photon fraction they predicted [65,66]. Figure 2.17: Photon limit measured by the Auger observatory. Black arrows are limits placed using measurements of the surface detector [65] while arrows tagged as FD are using measurements of the uorescence detector [66]. Other results from experiments as Haverah Park (HP), AGASA (A), Yakustk (Y) and AGASA-Yakustk (AY) are also shown as well as the photon ux predictions from several top-down models. 2.4.5 Neutrino limit A limit to the neutrino ux has been set by the Pierre Auger Observatory [67] studying upcoming showers from ντ interaction with the Earth's crust. As ντ may interact in the Earth's crust under the surface detector of the observatory, the resulting τ lepton can 2.4 Recent results of the Pierre Auger Observatory 55 enter the atmosphere above the array and its decay may be observed. The resulting shower will travel almost horizontally and with dierential features from usual hadronic showers. No ντ candidates have been observed yet and an upper limit to the neutrino ux has been set as shown in g. 2.18. Figure 2.18: Neutrino limit at the 90% of the condence level measured by the Auger observatory. Other results are also shown as well as the neutrino ux from GZK interactions. Also several studies are being carried to set a limit to the neutrino ux searching for deep horizontal showers induced by neutrinos in the atmosphere at large zenith angles [68]. 56 2 Cosmic rays features and the Auger observatory Bibliography [1] V. F. Hess, Phys. Z. 13 (1912) 1804. [2] R. Ulrich PhD. Thesis, Measurement of the proton-air cross section using hybrid data of the Pierre Auger observatory, Pierre Auger Collaboration Internal Document, GAP Note 2008-004. [3] G. A. Medina-Tanco, astroph/0607543. V. Berezinsky, astroph/07102750 [4] N. L. Grigorov et al. , Yad. Fiz. 11 (1970) 1058. T. Abu-Zayyad et al. , Astropart. Phys. 23 (2005) 157. V. A. Derbina et al. , Astrophys. J. 628 (2005) 41. T. Antoni et al. , Astropart. Phys. 24 (2005) 1. D. J. Bird et al. , Astrophys. J. 424 424 (1994) 491. T. T. Yamamoto, Proc. 30th International Cosmic Ray Conference, Mérida-México, 4 (2007) 335. [5] D. J. Bird et al. , Astrophys. J. 424 (1994) 491. [6] S. C. Corbato et al. , Nucl. Phys. Proc. Suppl. 28B (1992) 36. [7] J. Linsley, Phys. Rev. Lett. 97 (1955) 1292. [8] M.A. Lawrence et al. , J. Phys. G 17 (1991) 733. [9] M.M. Winn et al. , J. Phys. G 12 (1986) 653. 64 3 Arrival direction reconstruction. At the beginning of the fties, it was shown that the arrival direction of the cosmic rays can be estimated from the measurements of the arrival time of the signals at dierent points of an extensive air shower sampled by an array of detectors at ground level [3]. Since then, many other similar ground array experiments have been built, and all of them have used the arrival times of the shower front for the determination of the extensive air shower incoming direction. The incoming shower direction is obtained by tting a time associated to the arrival of the shower front at each triggered tank ( tmeas ) to that predicted according to a given model ( texp ) that depends on the incoming direction of the shower. Usually this tting is done by a minimum χ2 method: χ2= i=N X i=0 (∆ti)2 σ2 i = i=N X i=0 (tmeas i−texp i)2 σ2 i (3.1) where i runs through all the array detectors with a detected signal and σi corresponds to the uncertainty in the times. Within the Collaboration of the Pierre Auger Observatory, some works have been published with several reconstruction methods of the arrival direction [48]. They dier on the choices of the expected arrival time, texp i , and of the uncertainty on the arrival time, σi . In the analysis of inclined events presented in the following chapters, a dierent choice than the one made for vertical events will be used. This one is based on the model proposed in [2] and described in section 3.5 that predicts both the curvature of the shower front, texp i , and the time variance, σi . The three magnitudes involved in the χ2 minimization in the dierent approaches are addressed in the next sections. Section 3.2 deals with the start time and its measurement. Dierent procedures to estimate the expected time are explained in section 3.3. Section 3.4 is devoted to the variance on the start time. In order to justify the use of the time model described in 3.5, its prediction for the variance will be compared with to that used in vertical events in section 3.6 and the angular resolution obtained with our approach will be validated using hybrid events in section 3.7. 3.2 Start time Shower particles follow dierent distributions of arrival time at ground depending on the zenith angle of the shower and on the distance to the shower axis, so several choices can 3.2 Start time 65 be made to assign a measured time to the shower front at any distance from the shower axis. At a rst glance the mean of the distribution of arrival times could be thought as the natural choice, however this is not favored because the actual time distributions are not known. The onset of the time distribution is thought to be a better choice because an estimator of this time can be extracted directly from the data: the start time of the signal. However, it is a biased estimator due to the stochastic nature of the arrival of shower particles. Particles of the shower front will always arrive at the detector later than expected from the onset of the time distributions which ultimately must be related to causality. The measured start time of the signal is the arrival time of the rst particle which will become a better estimator of the onset of the time distribution as the number of detected particle increases. This is often referred as the sampling eect and it is potentially important for inclined showers because a single muon can produce a second level trigger. 3.2.1 Expectations The shower front is behind the shower plane and practically tangent to it at the shower axis so the direction normal to this plane can be determined from timing data at relatively small distances from the shower core. In general, it is convenient to measure the arrival time of the shower front at a given point with respect to the arrival time of this plane, tplane . Due to causality, shower particles can not reach the ground before the shower plane. The more the particles travel in straight lines at speed of light and near the shower axis, the closer they are to this imaginary front. Far from the shower axis, they accumulate more delay due to the angular deections in their interactions with the medium and to their subluminal velocities. The arrival time of the rst particle can be writen as: T1=tplane +T0+ts (3.2) where tplane is the arrival time of the shower plane, T0 is the delay of the shower front with respect to this plane and ts is the delay of the rst particle within the shower front. The delay of the shower front with respect to the shower plane, T0 , is related to the curvature of the shower front, whereas the delay of the rst particle with respect 66 3 Arrival direction reconstruction. to the shower front, ts can be related to the variance of the start time. This work will follow this distinction between curvature and variance, so the expected arrival time of the shower front at a given point will be texp =tplane +T0 with a variance V[t] = V[ts] . 3.2.2 Start time measurement The detection time used in eq. 3.1, tmeas , is determined from the trigger time of each tank. Detected events are built on a basis of local triggers as it has been described in the chapter devoted to the Auger observatory. For each tank that passes a T1 trigger, a trigger time is recorded from the time provided by the GPS receiver which is part of each station of the surface detector. This trigger time is related to the instant when the signal fullls one of the T1 trigger conditions. An algorithm that runs through the recorded signal traces searches the Start Bin associated to the onset of the signal. For each of the three PMTs, sections of the ADC trace with at least two consecutive bins with 3 counts over the baseline ( ≈0.06 V EMpeak ) are searched. Two sections are merged if they are closer than 20 bins and the ratio of their charges is greater than 30% , if not they are considered as dierent signals. Once the sections with signal have been found for each PMT, an inclusive merge of them is made between active PMTs: sections with bins in common are taken as parts of a single signal and the earliest (latest) of the bins is assumed to be the start (end) bin of the merged signal. The section with the greatest area ∗area peak is selected as the tank signal and its start and end bins are considered as those of the signal. This algorithm improves the start time determination for traces having signals previous to the trigger time and smaller than the trigger threshold and reduces uctuations on the start time associated with direct light. In addition to these small signals before the trigger time, it is possible that an isolated atmospheric muon crosses the tank just before the shower front. The start time of the signal will be misrecorded being assigned to the arrival of the isolated muon and not to the arrival of the shower particles. A trace cleaning algorithm has been developed in an attempt to keep the part of the trace that belongs to the detected shower and to remove the signal due to these isolated muons (See g. 3.1). The trace cleaning procedure described in detail in [9] is already implemented in the data acquisition software used for the work presented in this thesis [10]. The signal trace is searched for dierent 3.3 Fitting the start time 67 segments separated by gaps without signal. Segments of consecutive bins are classied as acceptable if their number of bins above a tuning threshold ( Nabove ) is greater than a minimum number ( Nmin ), or if the sum of the bin signals ( S ) is greater than a minimum signal ( Smin ). Both minima, Nmin and Smin , as well as the threshold level, can be tuned depending on the aim of the study: evaluation of the start time, rejection of exotics background, etc. The main segment which will dene the start bin is the one with the highest product of the number of bins above the threshold times their total signal. After this treatment of the recorded signal, a good fraction of the random muons is eliminated and the start time of the signal is a better estimator of the arrival time of the shower front. Figure 3.1: FADC trace showing an isolated muon prior to the main signal which distorts the signal start time measurement. 3.3 Fitting the start time The expected arrival times of the cascade at the locations of the stations are given by an assumption about the behavior of the shower front. The front of the cascade consists of the particles that rst pass through a given point which generate start time of the signal. Several hypotheses can be made to estimate the expected arrival time of the rst particles that reach each detector of the surface array being the simplest one a plane front ( T0= 0 in eq. 3.2). Other possibilities include spherical fronts, expanding spherical fronts, parabolic fronts, etc. Fig. 3.2 shows the propagation of the actual shower front and two approximations to its real behavior, the plane front and the spherical front. In 68 3 Arrival direction reconstruction. the work presented in this thesis we will use a model that predicts the mean delay of muons that reach ground level from [2]. Figure 3.2: The particle swarm of the cascade is shown in green with the shower axis in black. Plane and spherical approximations are shown in red. Two detectors of the array are shown in blue: the former triggered by the early part of the shower and the latter by the late part. Their distances to the shower core measured on the shower plane are also shown. 3.3.1 Shower plane front In the reconstruction of inclined showers performed in this work, the rst determination of the arrival direction is performed assuming that the shower front is a plane moving at the speed of light. It has been shown that for a given distance to the core, the delay of the particles with respect to this shower plane decreases with increasing zenith angle [11], so the assumption of a plane shower front is a better approximation for inclined showers than for vertical showers. The expected time can be written for each station i in terms of its position, ~ri as: texp i=t0−ˆn·(~ri−~r0) c (3.3) 3.3 Fitting the start time 69 where t0 is the arrival time of the shower at ground and ~r0 its impact point, ˆn is an unitary vector pointing towards the incoming direction of the shower and c is the speed of light. The impact point of the shower must be obtained in a more accurate way by other means. At the beginning of the reconstruction procedure, ~r0 , is approximated using the barycentre of the tank positions weighted by the square root of the integrated signals. The time dierence between predicted and measured times, for a plane front assumption can be expressed as: ∆ti(t0,ˆn) = tmeas i−texp i=ti−[t0−ˆn·(~ri−~r0) c] (3.4) This dierence between measured and predicted times is minimized with respect to t0 and ˆn , using a minimum χ2 procedure. For this minimization procedure the additional constraint that Pi=3 i=1 n2 i= 1 will be used. This constraint allows us to perform a three parameter t. Using this constraint, the expected time can be written for each station as: texp i=t0−ˆn·(~ri−~r0) c =t0−nx(xi−x0) + ny(yi−y0) + p1−n2 x−n2 y(zi−z0) c (3.5) 3.3.2 Front curvature Dierent studies of the time structure of extensive air showers [1114] have shown that arrival time at ground of the shower front cannot be modeled by a plane ( T06= 0 in eq. 3.2). A simple time delay using a spherical shower front with constant curvature can be implemented as T0 i=d2 i 2Rc , where di is the distance of the station i to the shower axis on the shower plane, R is the radius of curvature of the shower front and c is the speed of light. The radius of curvature, R , can be either tted through the minimum χ2 procedure or may be given in advance if the number of stations available is not high enough to allow a four parameters t. The produced delay with respect to the shower plane is symmetric around the shower axis. For this shower front model, the core position 70 3 Arrival direction reconstruction. is essential because the axis is set-up in the core position and adjusted in direction to t the station times. The time t and LDF t are entangled, since change in the core position requires an axis readjustment. This fact leads to a combined "global" t. However, we know that the shower front is not spherical and has not a constant curvature. On the contrary, the delay T0 in eq. 3.2 depends not only on the zenith angle of the shower and on the distance to the shower core but also on the azimuth of the observation point. In the analysis of vertical events, a spherical shower front model describing a dierence in early and late parts of the shower is used. This model incorporates the asymmetric radius of curvature by construction, since the front is not modeled as xed-radius sphere moving along the axis but it is described as concentrically growing spheres, i.e. like in an explosion. The radius of curvature of the shower front depends on time, and hence on the position at ground, as shown in g. 3.3. For this concentric shower front, no notion of the shower core is needed. The result of the t of the station timings is the shower "center. The shower axis, that it is not tted, is obtained as a normalized vector between the shower "center" and the core position after the core reconstruction. With this model, the LDF part and timing part are tted separately. Both ts are only entangled by the time variance model that contains a dependence on the zenith angle that couples the two parts of the reconstruction. The radius of curvature when the shower-front goes through the core position corresponds to the one obtained from the constant curvature t. The model employed in this thesis work allows us to calculate the mean delay of the muons with respect to the shower plane, which relates to the curvature of the front being the early-late asymmetry of the curvature taken into account. This delay does not arises from any assumption about the shower front structure (plane, spherical, parabolic,...) but from physical arguments about the propagation of the muonic component of the shower in the atmosphere. Measured times are then corrected for the mean delay T0 as predicted by the model and shown in g. 3.5 for certain measurement conditions. These corrected times are tted to a plane front to obtain the incoming direction of the shower. This time model needs the core position to calculate the delays but, besides being these delays based on physical assumptions about the passage of muons through the atmosphere, the model also provides the start time variance to be used in the t. 3.4 Start time variance 71 Figure 3.3: Concentric shower front arriving at two tanks (blue squares) of the surface array. The red line is the shower axis and the yellow star tags the shower center. Dierent radii of curvature of the shower front will be seen from the early or the late part of the shower. 3.4 Start time variance Drawing meaningful physical results from the χ2 tting procedure requires a proper evaluation of the uncertainty of the incoming direction of the detected air shower. The precision of the arrival direction reconstruction, depends not only on the precision of the clock which measures the arrival times, but also on the uctuations in the arrival time at ground of the particles which are part of the shower front. The procedure to obtain the arrival direction from the surface detector data, based on the tting of the start times of the signals in the tanks to the predicted times from a propagation model of the shower front, has to properly take into account the uncertainties in the start time measurement. The uncertainty in the reconstructed direction, extracted from the tting procedure, is directly related to these uncertainties so to have 72 3 Arrival direction reconstruction. an adequate estimate of that uncertainty, uctuations in the arrival time have to be properly modeled or parameterized. The start time variance for each signal that enters eq. 3.1, σi , takes into account these uctuations and should give a proper weight to each station in the reconstruction procedure. The uncertainty in the star time can be considered as a delay in the arrival time of the rst particle with respect to a theoretical shower front ( ts in eq. 3.2) and it is related to the arrival time distribution of the particles at ground. The variance of the start time has a contribution due to the precision of the time measurement device. This contribution is related to the GPS resolution ( 10 ns ) [15] and to the FADC sampling ( 25 ns ) has been estimated in [16] as: V[ts]timing = 102+25 √122 ≈150 ns2 In addition to this estimate there are uncertainties associated to the position of the impact point of the particle in the tank which are not accounted for. Besides, there are uctuations associated to the arrival of the rst particles of the shower front at the tanks. These uctuations, or sampling eect, increase with the distance to the shower core and decrease with the particle density, thus they dominate for stations far away from the shower core. The contribution to the variance of the start time due to the sampling eect, can be obtained from statistical models with a set of assumptions about the distribution of particles in the shower front as it is done in [16,17] or can be computed from the actual distributions obtained from a model based on the work described in [2]. The parameterizations for the start time uncertainty proposed and tested for the Auger Observatory surface detector tanks in [16,17] extend, in principle, to the whole zenith angle interval so they can be used, a priori, for the analysis of inclined events. Given a normalized distribution of arrival times of the particles in the shower front, f(ts) , the probability that the rst particle out of n arrives before a time τ is the complementary to the probability that all the n particles arrive later than τ . F1(τ) = 1 −(1 −F(τ))n (3.6) where F(τ) is the probability that a single particle arrives before a time τ and it is given by: 3.4 Start time variance 73 F(τ) = Zτ 0 f(ts)dts (3.7) The distribution of the delay of the rst particle, τ1 , is just the derivative of F1(τ) at τ1 g(τ1) = dF1 dτ τ=τ1 =n(1 −F(τ1))n−1f(τ1) (3.8) Assuming dierent time distributions of the particles within the shower front ( f(ts) ), the expected value and the variance of the delay of the rst particle with respect to the assumed front can be evaluated directly from the previous distribution g(τ1) . E[τ1]sampling =Z∞ −∞ τ1g(τ1)dτ1 (3.9) V[τ1]sampling =Z∞ −∞ (τ1−E[τ1])2g(τ1)dτ1 (3.10) For instance, assuming a constant rate of particles within a time interval T (box distribution), the arrival time distribution would be f(ts) = 1 T and hence the expected value and the variance for the delay of the rst particle would be: E[τ1]sampling =T n+ 1 (3.11) V[τ1]sampling =T n+ 12n n+ 2 (3.12) Both models proposed for vertical shower analysis [16,17] are based on this assumption of a constant rate arrival of particles. They dier on the time interval considered ( T in eq. 3.11) and on the dependence on the number of particles. The model proposed in [16] is given by: V[ts]sampling ≡V[τ1]sampling =T50 +C n+ 1 2n n+ 2 (3.13) where n is the number of particles crossing the tank calculated from the total integrated signal and T50 is the risetime of the signal dened as the time it takes the integrated 80 3 Arrival direction reconstruction. Figure 3.7: Distribution of the zenith angle of the reconstructed events with at least one doublet or triplet triggered. Both stations in the same doublet or at least two in the same triplet are requested to be red by the shower. Signals and distances to the reconstructed core of both tanks in the doublet and the dierence between both time measurements are shown in g. 3.8 where solid red and dashed blue lines are used respectively for the rst or the second detector triggered in the pair. From these plots, the following cuts are applied to enhance the quality of the data set getting rid of the outliers of the distributions which are assumed to be uctuations, badly reconstructed events or events for which random muons spoil the start time determination in one or more tanks. 1. The signal in each tank must be below 500 V EM to get rid of saturation eects (top panel g. 3.8). 2. The distance of both stations to the shower core on the shower plane must be over 100 m to get rid of the steepness of the lateral development of the shower close to its core (middle panel g. 3.8). 3. The dierence between start time measurements must be below 100 ns to ensure that both tanks are measuring the same shower and not an isolated muon or a coinciding smaller air shower (bottom panel g. 3.8). 3.6 Testing variance models for inclined events with twin tanks 81 Around 10 % of the events are rejected when applying these cuts. Figure 3.8: Top: distribution of the signals in each tank of the doublet. Center: distance to core distribution. Bottom: start time dierence distribution. 82 3 Arrival direction reconstruction. For each model of the start time uncertainty [2,16,17], the dierence in the residuals of the time t is compared with the variance. The dierence in the start time measurement has to take into account the shower front propagation which is model dependent whereas, as the two tanks in the doublet are close, the dierence in the time t residulas is practically independent of the shower front model assumed for the angular reconstruction and, if the uncertainties in the time measurement are correctly described, the distribution of δ dened by eq. 3.19 should have unit variance. δ=∆T pV[T1] + V[T2] (3.19) where ∆T is the dierence between the two residuals, and V[T1] and V[T2] are the variances for both start times predicted by the model. Models [17] and [16] given by eqs. 3.16 and 3.13 need a rst guess of the incoming direction of the shower to calculate the start time variances. These parameterizations are done on a number-of-particles basis and they use the mean track length, which is zenith angle dependent, to calculate the number of particles from the measured signals 1 . This rst guess of the arrival direction has been obtained for this analysis making a time t to a plane front with a constant start time variance of V[T] = 625 ns2 . Model [2] besides needing the incoming direction of the shower to compute the start time uncertainty also needs the core position. In this case, the shower core position is reconstructed with the algorithm used for the event reconstruction, so an iterative angular-core reconstruction is implemented. The uncertainty introduced in the measurement by the detection device has been introduced as a constant term of 20 ns that has been added in quadrature to the uncertainty associated to the width of the shower front which is directly predicted by the model. This term takes into account the uncertainty introduced in the start time measurement by the light collection in the tank, by the GPS time determination and by the procedure of signal digitalization. To check the validity of the proposed models [2, 16, 17] for the angular range of interest to the analysis of inclined events, the distribution of δ is used. This distribution is shown in g. 3.9 for all the three models. 1 n=S RT L(θ) . RT L(θ) = T L(θ) T L(0) 3.6 Testing variance models for inclined events with twin tanks 83 Figure 3.9: ∆T √V[T] distribution for the three evaluated models. Top: [17], center: [16] and bottom: [2]. From plots in g. 3.9 we see that both parameterizations from twin tanks measure- 84 3 Arrival direction reconstruction. ments (models [17] and [16]) are underestimating the start time variance for inclined events, producing a width of the δ distribution greater than one. The model based in [2] shows RMS distribution close to unity. To validate the proposed model for dierent variables, we plot the width of the δ distribution as a function of the mean signal and of the distance to the core on the shower plane. In g. 3.10 we see the at behavior around the unit value which conrms the validity of the model for inclined events for a great range of values. Figure 3.10: RMS of ∆T √V[T] distribution for the [2] model. Top: signal dependence, bottom: distance to core dependence. The RMS dependence on signal and distance to the core for models [17] and [16] 3.6 Testing variance models for inclined events with twin tanks 85 are shown in gs. 3.11 and 3.12. A worse behavior is seen for large signals and for tanks close to the shower axis. Figure 3.11: RMS of ∆T √V[T] distribution for the [17] model. Top signal dependence, bottom distance to core dependence. 86 3 Arrival direction reconstruction. Figure 3.12: RMS of ∆T √V[T] distribution for the [16] model. Top signal dependence, bottom distance to core dependence. 3.7 Angular resolution studies with hybrid events A procedure to take advantage of the hybrid performance of the Pierre Auger observatory, has been set up to check the reconstruction chain. Besides, the use of high quality hybrid events, will allow us to improve the reconstruction algorithms. To extract the angular resolution of the hybrid events, articial showers were generated by laser shots. The Central Laser Facility, shown in g. 3.13, is located in the middle of the array at about 30 km from each uorescence detector. It contains a remotely controlled laser which produces almost vertical showers (within 0.01◦ ) sending at the same time a pulse of light to a surface station to generate an articial hybrid 3.7 Angular resolution studies with hybrid events 87 event [18]. The reconstruction of laser shots is done with the same algorithm used to reconstruct real hybrid events with only one surface station. The obtained angular resolution for laser shots is ∼0.3◦ [19]. Figure 3.13: The CLF with the Celeste SD station. The laser events have geometrical characteristics that are favorable for a better hybrid reconstruction of the geometry compared to real hybrid events. The angular resolution of actual hybrid events has to be modied to account for other uncertainties arising from the particular geometry of the shower development with respect to the uorescence detector that are ignored for laser shot events. One of them is the uncertainty of the hybrid shower axis which can be estimated as the ratio between the core uncertainty and the shower axis lever arm 2 . The mean accuracy for real hybrid events can be quoted as ∼0.6◦ from a mean shower lever arm of 7.5km [20], so we can use the incoming direction from the hybrid reconstruction as a very good estimate of the true one and compare the SD-only incoming direction with it. Quality cuts imposed to select high quality hybrid events are: 1. Golden events from 1st Jan 2004 to 31st May 2008. 2. FD arrival direction reconstruction from ADST v5r1. 3. Reduced χ2 of the time t lower than 5 . 4. Zenith and azimuth uncertainties below 5◦ . 2 the distance along the shower axis from ground to the point seen in the camera by the pixel with the highest elevation. 88 3 Arrival direction reconstruction. With those cuts, more than 2000 golden hybrid events are selected to test the performance of the angular reconstruction algorithms. The comparison between zenith angles of the SD only reconstruction and the hybrid reconstruction is shown in g. 3.14. A better resolution of the zenith angle is shown by the analysis with the curved shower front than that observed for reconstruction with a simple plane shower front. Figure 3.14: Zenith angle comparison between hybrid and surface only reconstructions. Left: plane shower front. Right: shower front with curvature and time variance from [2]. The space angle between both reconstruction is shown in g. 3.15. The distribution has been tted to a Gaussian resolution function dp ∝e−ψ2/2σ2d(cos(ψ)) dφ where ψ is the space angle between both reconstructions. The angular resolution (AR) is related to the σ parameter as AR = 1.5σ . 3.7 Angular resolution studies with hybrid events 89 Figure 3.15: Space angle between hybrid and surface only reconstructions. Left: plane shower front. Right: shower front with curvature from [2]. Considering a resolution of ∼0.6◦ for the hybrid events, the angular resolution for the surface detector can be quoted for this study as ∼1.1◦ . As it is shown in g. 3.16 no systematic is present in the space angle due to the zenith angle, the distance to the axis of the surface detector tank included in the hybrid t, to the uorescence eye or to the acquisition date. 96 4 Impact point and shower size reconstruction ∂ ∂N19 i=n X i=0 (Smeas i−N19Sexp i(1,ˆxi,ˆyi))2 σ2 i!= 0 i=n X i=0 −2(Smeas i−N19Sexp i(1,ˆxi,ˆyi))Sexp i(1,ˆxi,ˆyi) σ2 i = 0 i=n X i=0 Smeas iSexp i(1,ˆxi,ˆyi) σ2 i−N19 i=n X i=0 Sexp i(1,ˆxi,ˆyi)Sexp i(1,ˆxi,ˆyi) σ2 i = 0 N19 = i=n X i=0 Smeas iSexp i(1,ˆxi,ˆyi) σ2 i i=n X i=0 (Sexp i(1,ˆxi,ˆyi))2 σ2 i!−1 (4.4) The N19 value that minimizes the χ2 as computed in eq. 4.2 can be associated to the core position using the precedent equation. Using the barycentre of the triggered tanks weighted by their signals as a rst estimate of the shower core position, a search of the reconstructed core is performed around it. Sampling X and Y values across a grid around the barycentre, the N19 value that minimizes χ2 according to eq. 4.4 and its associated χ2 value according to eq. 4.2, can be calculated for each point of the grid. Two grids are set up for the search of the core position: the rst one has a spacing of 50 m and the second one of 20 m . This last choice matches the spacing used for the muon maps generation. This second minimization is not built around the barycentre but around the position in the rst grid that minimizes χ2 according to eq. 4.2. The core position (X, Y ) and the associated N19 of this second minimal χ2 search are used as input for the nal reconstruction using the maximum likelihood method. Just as an example, a χ2 map for a real event obtained using the grid search procedure is shown in g. 4.1. For this high multiplicity event the minimum of the χ2 associated with the position of the shower core is clearly visible. 4.1 Core position reconstruction with a minimum χ2 method 97 Figure 4.1: χ2 map for a high multiplicity event resulting from the rst grid search (see text) of minimal chi2 . 4.1.1 Signal uncertainty The uncertainty in the signal measurement to be used in eq. 4.2, σi , can be evaluated from signal measurements of doublets as stated in [5] and parameterized for vertical showers as: σi= (0.32 + 0.42secθ)pSmeas i (4.5) In our case, the uncertainty in the detected signal for this χ2 method is based on the Poisson statistics on the number of muons neglecting the electromagnetic correction to the signal and it is given by: 98 4 Impact point and shower size reconstruction σi=pNµi Sexp i Nµi =pTL(θ)qSexp i (4.6) where Nµi is the expected number of muons, TL(θ) is the mean track length for a given zenith angle and Sexp i=NµiTL(θ) is the expected signal. The use of the expected signals instead of the measured ones is done to allow a similar treatment of triggered and zero signal stations. Predicted signals for stations far away from shower axis are not very reliable because muon density maps are restricted in size to a 8km ×8km square. Outside this area the expected signal is hard to handle and it is set to zero. Typically the stations that are outside this area do not have signal and zero signal stations with a null expected signal are not included in the χ2 calculation. 4.1.2 The barycentre as rst guess of the core position The choice of the barycentre as a rst estimate of the shower core compared to the nal reconstructed core is shown in g. 4.2 where the average distance from the reconstructed core to the barycentre measured in the ground plane is plotted versus the reconstructed zenith angle for a large set of T5 events. The mean distance for all the considered events is below 500 m spanning from 200 m for events with a zenith angle of 60◦ to over 1500 m for events with a zenith angle greater than 80◦ , although the statistics rather poor above 84◦ . From g. 4.2 we see that as the reconstructed zenith angle increases, the relation of the impact point with the barycentre of the signals is fainter. This is one of the reasons to not use the vertical denition of the highest trigger level ( T5 or quality trigger) for inclined events. 4.1.3 Importance of zero signal stations In the calculation of N19 and of the core position using the χ2 method, zero signal stations play a very important role to conne the impact point of the showers, above all for low multiplicity events (events with a small number of tanks with signal). In g. 4.3 χ2 maps for the same event using and not using zero signal stations are shown. As it is clearly seen in the left panel, when no zero signal stations are used, there is a two minima structure which can lead to a bad reconstruction of the core position. The instability is 4.1 Core position reconstruction with a minimum χ2 method 99 Figure 4.2: Average distance on ground of the reconstructed core to the barycentre. Barycentre of signals seems to be a good choice as a rst estimation of the shower impact point at ground. eliminated just using all zero signal stations within 1600 m from all the triggered tanks. This is shown in the χ2 map on the right panel that only displays one minimum very close to the true core position which is surrounded by the four triggered stations. This situation is rather common particularly for low multiplicity events. Fig. 4.3 also gives a hint about the importance of the high level trigger which demands that all the stations around that closest to the shower core the core are active to ensure an accurate reconstruction. If the shower core is not surrounded by a complete hexagon of active stations, as can happen when the shower strikes the ground near the edges of the surface detector or close to a position without a working tank, the core position and hence the shower size can be misreconstructed. This situation is somehow reproduced when zero signal stations are not accounted for in the reconstruction and the undesirable result is clearly seen in the left panel of g. 4.3. 100 4 Impact point and shower size reconstruction Figure 4.3: χ2 map for N19 and core reconstruction without (left) and with (right) zero signal stations. Black dots are the triggered stations while the grey dots on the right panel are the zero signal stations considered in the χ2 minimization. 4.2 Shower size estimate using maximum likelihood The method of maximum likelihood estimation (MLE) is a statistical tool rst introduced in [6] and used to calculate the best way of tting a mathematical model to some experimental data. With this method we try to maximize a probability density which must depend on the parameters we want to obtain. In our case, where the parameters we want to obtain are the shower size ( N19 ) and the position of the shower core ( X , Y ), this probability function is built from the product, for all the considered tanks, of the probability of having a measured muonic signal Sµmeas i 1 when a number of muons Nµi is expected. P(N19, X, Y ) = n Y i=1 P(Sµmeas i, Nµi) (4.7) This probability, P , depends on the searched parameters, N19 and (X, Y ) because the expected number of muons, Nµi , depends on them. Nµi depends, for a given tank, on its position with respect to the shower axis measured on the transverse plane, (ˆx, ˆy) , and on the shower size characterized by the parameter N19 . The maximum likelihood 1 We use the term measured muonic signal to refer to the muon signal extracted from measurement obtained after the electromagnetic correction has been subtracted. 4.2 Shower size estimate using maximum likelihood 101 procedure will search a shower size value and a shower core position that produce an expected number of muons at each station which maximizes the total probability stated in eq. 4.7. The two dimensional distributions of the muons projected onto the transverse plane, also called muon maps, are obtained from simulations. The same assumptions made for the minimum χ2 procedure are used here and, for a xed arrival direction, the shape of the muon maps is assumed to be independent of energy which only aects the global normalization N19 . Also, the shape of the maps is assumed not to depend on composition nor on hadronic models [3,4]. We can factor out the normalization N19 relative to a reference distribution (see section 4.3.1). The signal due to the muonic component, Sµmeas i , is obtained subtracting the average electromagnetic component in the measured signal. In this work, this is parameterized as a factor to apply to the signal, which depends on zenith angle and on the distance to the shower core, EMCorrection(ˆx, ˆy) . It will be discussed in section 4.3.3. Nµi≡Nµi(N19,ˆxi,ˆyi) = N19 Nµi(1,ˆxi,ˆyi) Sµmeas i≡Sµmeas i(ˆxi,ˆyi) = Smeas i 1 + EMCorrection(ˆxi,ˆyi) The measured signal Sµmeas i when Nµi muons are expected, can be produced by dierent numbers of muons, each with dierent probability. We dene Pk(Sµ) as the probability distribution of k muons to produce a muonic signal Sµ when traversing a surface detector tank. To compute each local probability, P(Sµmeas i, Nµi) , the sum over the probabilities for all possible numbers of muons k to produce the measured signal, Pk(Sµ) , has to be computed. Each probability must be multiplied by the Poissonian probability of having k muons when Nµi are expected. The overall probability becomes then: P(N19, X, Y ) = n Y i=1 P(Sµmeas i, Nµi) (4.8) = n Y i=1 k=∞ X k=1 Poisson(k, Nµi)Pk(Sµmeas i) (4.9) 102 4 Impact point and shower size reconstruction The probability functions Pk(Sµmeas) are discussed in section 4.3.2. A special comment needs to be made about the treatment of zero signal stations and stations with a saturated signal. To a rst approximation, only an upper o lower limit to the signal is available by the measured value of these stations. For no signal stations, the only information we have is that the measured signal is below a trigger threshold ( Sth ) given, in principle, by the rst level trigger. The contribution of a zero signal station to the total likelihood probability can be evaluated as the probability of having k muons when Nµ are expected multiplied by the probability that these k muons produce a signal below Sth . The probability demands an integral of the probability distribution over all the possible signals below the threshold (i.e. in S ): Pzero(Nµi)≡ P(S < Sth, Nµi) (4.10) =ZSth 0P(S, Nµi)dS (4.11) =ZSth 0 k=∞ X k=0 Poisson(k, Nµi)Pk(S)dS (4.12) = k=∞ X k=0 Poisson(k, Nµi)ZSth 0 Pk(S)dS (4.13) The T1 threshold is of about 1.75 V EM while the ToT is 1.6V EM so the number of muons that can produce such a faint signal has to be quite low. The sum over the possible number of muons, k , in practice can be reduced to below 3. The number of stations which have no signal and are included in the t is selected ensuring that all stations within a distance of 6700 m to a triggered tank are included in the maximum likelihood t. When a signal is saturated, eectively only a lower limit of the actual signal is available. The probability is similarly evaluated as the probability that k muons produce a signal above Sµmeas i : 4.2 Shower size estimate using maximum likelihood 103 Psat(Sµmeas i, Nµi) = P(S > Sµmeas i, Nµi) (4.14) =Z∞ Sµmeas iP(S, Nµi)dS (4.15) =Z∞ Sµmeas i k=∞ X k=0 Poisson(k, Nµi)Pk(S)dS (4.16) = k=∞ X k=0 Poisson(k, Nµi)Z∞ Sµmeas i Pk(S)dS (4.17) Saturated signals are likely to be produced by large numbers of muons. As will be shown, Pk(Sµ) can be accurately approximated by a Gaussian function for k > 8 . Besides, the Poissonian distribution can be also approximated by Gaussian distribution for k high enough and hence this integral can be parameterized as: Psat(Sµmeas i, Nµi) = 1 21−tgh Sµmeas i−S Nµi 2S2√Nµi (4.18) Saturation eects could arise from saturation of the sampling ADC and/or of the PMT. The rst one occurs at typically ρ≈100 V EM/m2 and the second one at ρ≈500 V EM/m2 . A procedure to take into account these two kinds of saturation eects and to deduce the actual signal from the FADC trace obtained has been developed in [7] but it has not been used in this work. Saturated signals have been treated as lower limits of the measured signals. In the implementation of this method we minimize the opposite of the natural logarithm of the total probability. The upper limits of the involved summations are truncated to reduce computation time. −log(P) = − n X i log(P(Sµmeas i)) (4.19) =− n X i log( k=u−lim X k=l−lim Poisson(k, Nµi)Pk(Sµmeas i)) (4.20) Fig. 4.4 shows the likelihood function behavior around the minimum for two real 104 4 Impact point and shower size reconstruction events. Figure 4.4: Likelihood function maximization for two events detected by the Auger surface detector. The negative of the natural logarithm of the likelihood function is minimized with respect to the logarithm of the size parameter N19 . The impact of zero signal stations in the likelihood maximization is shown in g. 4.5. Besides an overall decrease in absolute P for the event which is irrelevant, we note that adding zero signal stations has a signicant eect on the extracted shower size. For the event shown, zero signal stations represent a reduction of over ∼20% in the N19 estimate. This is characteristic of the lower energy events 2 . 2 GAP note in progress 4.3 Calculation of the expected signal 105 Figure 4.5: Comparison between likelihood maximization with (blue) and without (red) zero signal stations is shown. A reduction of ∼20% in the N19 estimate can be observed for this event. 4.3 Calculation of the expected signal The expected signals are needed for the comparison with measured ones in the χ2 minimization or to evaluate the probabilities of having a given measured signal needed in the maximum likelihood method. They are obtained from the expected number of muons and from the signal that they produce traversing an surface detector tank. The minimization of the χ2 given by eq. 4.2 only needs the mean signal produced by an expected number of muons and its uncertainty while for the maximum likelihood method the complete signal probability distribution for this number of muons has to be given to construct the complete probability function as stated in eq. 4.7. 112 4 Impact point and shower size reconstruction 2. Signal distribution probability for k muons. To calculate the signal probability when k muons enter the tank, Pk(Sµ) ,we have to perform the autoconvolution of the probability distribution function of a single muon. P2(Sµ) = Z∞ 0 dS0P1(S0µ)P1(Sµ −S0µ) ··· ··· Pk(Sµ) = Z∞ 0 dS0P1(S0µ)Pk−1(Sµ −S0µ) Although the shape of the signal probability for a single muon is complicated, the function rapidly becomes Gaussian allowing a simple parameterization of its mean value and its width as functions of the number of muons. We have chosen to use the gaussian approximation when the expected number of muons is greater than 8. Pk(Sµ) = Gauss(ζk(E, θ), σk(E, θ)), if k > 8 (4.22) As it is shown in [11], on the one hand, the mean signal of k muons, ζk(E, θ) , grows linearly with the number of muons, k ,for xed zenith angle and energy. On the other hand, its uncertainty, σk(E, θ) , grows with its square root. The normalization of these dependences are functions of the zenith angle and of the muon energy that can be factorized. The mean and width of the gaussian can be parameterized as follows: ζk(E, θ) = f1(E)TL(θ)k (4.23) θk(E, θ) = f2(E)σ2 TL(θ)√k (4.24) where k is the number of muons, E is their energy, θ is the entering zenith angle, < TL > is the mean track length at that zenith, σ2 TL is the track length variance 4.3 Calculation of the expected signal 113 and f1 and f2 are two energy dependent expressions given by: f1(E)=0.92 + 0.085 log10(E) (4.25) f2(E)=0.23 + 0.052 log10(E) (4.26) 4.3.3 Electromagnetic part of the signal Signals produced at the surface detector tanks are dominated by muons in horizontal showers. The electromagnetic component due to the π0 decays is almost completely absorbed before reaching the ground particularly for angles above 65◦ . In any case, this component remains typically at low distances to the shower axis. In addition, there is a small electromagnetic component in inclined showers due to muon decay and to hard muon interactions. This component, that contributes typically about 20% of the muonic signal [1] has been studied in [16]. Both the size and its two dimensional behavior have to be taken into account. The measured signal is corrected subtracting the electromagnetic correction to obtain the muonic signal in order to compare the muonic signal measured to that expected according to the prediction for the number of muons. The electromagnetic component due to muon decay has been shown to follow the two dimensional distribution of the muons [1]. The signal induced by the electromagnetic component was parameterized taking the ratio of the electromagnetic signal, SEM , to the muonic signal, Sµ as a function of the zenith angle and on the distance to the shower axis. The parameterization of this ratio [16], is obtained using shower simulations with AIRES and a tank simulation code called S1000 USC also described in that work. In g. 4.9 the ratio of electromagnetic to muonic signal is shown as a function of the distance to the shower axis for several zenith angles. On the one hand, close to the shower core, this ratio decreases with zenith angle up to 72◦ increasing again from there. This behavior is due to the π0 shower component which rapidly attenuates as the zenith angle increases from 60◦ to 72◦ to be taken over by the hard muon processes that dominate the electromagnetic signal close to the shower core for very inclined showers. On the other hand, far from the shower core, the contribution of the electromagnetic component is only due to muon decay and becomes practically constant. 114 4 Impact point and shower size reconstruction Figure 4.9: EM ratio. A simple parameterization of this ratio has been made in [16]: SEM Sµ (r, θ) = A(θ)rC(θ)−B(θ) log10(r) (4.27) The assumption that the ratio SEM Sµ is equal at the same distance from the shower axis in the shower plane regardless of the azimuthal angle is only an approximation. There is an azimuthal asymmetry in the signal due to the several eects, the most important are the so-called geometrical eect, the longitudinal development eect and ground screening [17]. At θ > 60◦ there is a small early-late asymmetry in the muonic signal that increases slowly with distance to the core. This is due to the long attenuation length of the muons in the atmosphere. At θ < 70◦ there is an important early-late asymmetry in the EM signal that increases rapidly with the distance to the shower axis. However, 4.4 Size correction due to the binning in the muon maps 115 at θ > 70◦ , the asymmetry increases slowly with the distance to the core in the same way as the asymmetry in the muonic signal since the EM component comes mainly from muon decay. The ratio SEM Sµ has a clear azimuthal asymmetry at θ < 70◦ in agreement with the asymmetries obtained for the muon and electromagnetic components. The ratio with these eects taken into account has been also parameterized in [16] as: SEM Sµ (r, θ, ζ) = SEM Sµ (r, θ) (1 + Aasym(r, θ, ζ)) (4.28) where θ is the zenith angle of the shower, ζ is the azimuth angle in the shower plane, r is the distance to the shower axis, SEM Sµ(r, θ) is given by eq. 4.27 and the parameter Aasym characterizes the azimuthal asymmetry. This correction is has not been used in this work. The eect of the magnetic eld in this parameterization has also been studied in [16] can be introduced a posteriori. It is negligible below 80◦ although its eect can be greater than 20% for showers with a zenith angle over 86◦ . It is neglected in this work. In addition, the electromagnetic correction changes with model, compositions and with energy. A systematic study of these eects is in progress but will not be reported in this thesis. 4.4 Size correction due to the binning in the muon maps Muon maps are binned for practical reasons, so detected showers are reconstructed using the available map closest to the zenith angle of the shower. As the shower size drops with zenith angle in an implicit way for each set of maps, a correction factor, κ , must be included in the analysis to account for this eect due to the dierence between the reconstructed zenith angle and the zenith angle of the used map. As a result, the N19 value for a shower arriving with a zenith angle θ and with a reconstructed shower size of N19rec computed with a map generated with a zenith angle of θmap can be obtained from the expression: N19 = N19rec κ(θ, θmap) (4.29) The parameterization of the κ factor has been obtained for the set of muon maps used 116 4 Impact point and shower size reconstruction in this work directly from full Monte Carlo air shower simulations by comparing the total number of muons as a function of the zenith angle for simulated showers of the same energy. In g. 4.10, a parameterization of the total number of muons with respect to a 60◦ shower is plotted versus the zenith angle. As an example, all the detected showers with a zenith angle between 63◦ and 65◦ (dashed red lines) will be reconstructed with the 64◦ map. As the total number of muons of a 63◦ ( 65◦ ) shower is greater (smaller) than that of the used map, the size parameter obtained will be then smaller (greater) than the true one and hence the correction factor κ will be greater (smaller) than one for reconstructed zenith angles below (above) the nominal value of the used map. Figure 4.10: Attenuation of the total number of muons from maps. The correction factor to be applied in each case shown in g. 4.11 can be parameterized as follows: κ(θ, θmap) = 10d(θ)−d(θmap) (4.30) 4.5 Zenith angle dependence evaluation 117 where d(θ) = a+b[1 −cos(θ)] + c[1 −cos(θ)]2 (4.31) and a= 0.395 , b=−0.04 and c=−1.54 . The correction factor increases for larger zeniths and reaches over 12% for very inclined events as it is shown in g. 4.11. Figure 4.11: Correction factor. As a corollary, the uncertainty in the determination of the zenith angle propagates to the estimated N19 . The attenuation shown in g. 4.10 can be used in a similar fashion to calculate the eect of this propagation of the angular uncertainty to the N19 and hence to the estimated energy. 4.5 Zenith angle dependence evaluation As the development stage of an extensive air shower is a function of the atmospheric slant depth traveled before reaching ground, showers with the same energy but with 118 4 Impact point and shower size reconstruction dierent zenith angle of incidence, will have dierent number of muons at ground and will show dierent shower sizes. Following analyses of more vertical showers we can refer to the reduction of the signals with zenith angle as attenuation. The attenuation curve is just a function that describes the dependence of this reduction on the zenith angle. In our case the size parameter N19 is relative to a given map corresponding to the arrival direction of the reconstructed event. As the maps include an implicit attenuation, N19 only depends on energy assuming that maps have the correct attenuation. This assumption can be tested using the sin2(θ) distributions. Assuming an isotropical arrival of cosmic rays and 100% eciency of the detector above a given energy threshold, the zenith angle distribution of events should proportional to sin2(θ) . Above a certain energy threshold, which is related to the saturation of the detector eciency, this distribution should become at. Equivalently a Constant Intensity Cut method [1820] can be made by chosing the N19 value above which the event rate (or intensity) obtained in equal bins of sin2(θ) is constant. The resulting graph for N19 should be at if the array is 100% ecient. Figure 4.12: Solid lines are ts to a rst order polynomial function. Legend shows the slopes of this ts, all of them compatible with a constant attenuation curve. Reduced χ2 values for each t are also shown. 4.6 Energy calibration 119 The corresponding plot obtained for the ve dierent intensities shown in the previous gure are practically at what conrms that attenuation is adequately taken into account by the maps within the statistical accuravy of the method. 4.6 Energy calibration As the shape of the muon maps is practically independent of energy, the size of the detected showers is measured with respect to a reference simulation corresponding to 1019 eV proton air showers. This size parameter, denoted as N19 , gives the relative number of muons of the detected shower with respect to the reference for the same arrival direction. N19 will be shown to have a good correlation with the uorescence energy measurement for a subset of hybrid events. This correlation allows us to use N19 as an estimator of the true energy of the cosmic ray. 4.6.1 Hybrid events selection for calibration curve In order to assign an energy value to each N19 reconstructed value, we make use of the golden hybrid events reconstructed with both the FD and the SD to correlate the energy measured by the FD detector to the reconstructed N19 obtained with the SD in a similar way to what is done in vertical events [21]. From January 2004 to June 2008, 2665 T4 golden hybrid events have been detected with a zenith angle above 60◦ . A set of high quality golden hybrid events is chosen by imposing a collection of selection criteria. Quality criteria are applied to both the FD and the SD reconstructions. 1. FD cuts. (a) The events must have a Cherenkov fraction below 50% to avoid the uncertainties from the model of Cherenkov light. (b) The tank with largest signal must be closer than 750 m to the shower axis. (c) The reduced χ2 of the Gaisser-Hillas t must be below 4 to get rid of events with proles distorted due to clouds or fog. Proles without a clear maximum are rejected by the requirement that the χ2 of a linear t to the prole exceeds the χ2 of the Gaisser-Hillas t by at least 4 . 120 4 Impact point and shower size reconstruction (d) A precise reconstruction of the shower prole, and hence a reliable estimate of its energy, requires the observation of a signicant fraction of the prole and not only of its rising or falling section. For this reason, the shower maximum is required to be within the eld of view and determined with an uncertainty below 50 g/cm2 . (e) The uncertainty in the reconstructed energy must be below 40% . This is a useful parameter to reject badly reconstructed events because the reconstruction algorithm propagates uncertainties on the light ux and the geometrical uncertainties. The magnitudes of these cuts are similar to those applied in the calibration with vertical showers [21]. Table 1 shows the number of events rejected by each cut and the number of remaining events after applying the cuts consecutively. Cut label a b c d e Only 801 520 1189 1822 367 Remaining 1864 1440 759 439 439 Table 4.2: First line: FD cut label. Second line: number of events rejected by the each of the cuts applied on their own to the initial 2665. Third line: number of remaining events after applying cuts consecutively. 2. et cuts. (a) The relative uncertainty in N19 must be below 40% . (b) The inclined T5 quality trigger condition must be satised to get rid of badly reconstructed events due to a misplacement of the shower core. From the 439 events selected by the uorescence cuts, 421 have a statistical uncertainty on N19 below 40% and 286 out of them fulll the T5 trigger condition. A good linear correlation between the size parameter N19 and the uorescence energy EFD is shown in g. 4.13 for the 286 selected events. The energy of the event can be determined with the surface detector with this procedure that eliminates many of the 4.6 Energy calibration 121 dependences on simulation or primary composition uncertainties from the shower size parameter by the linear relation: log(N19) = A+Blog(E) (4.32) Figure 4.13: Linear correlation between log10(E) and log10(N19) for the 286 selected events. Statistical uncertainties on N19 and on EFD have been assigned to each event. On the one hand, for the statistical uncertainty of the SD reconstruction, σN19 , three contributions have been taken into account: σ2 N19 = (σmle N19)2+ (σθ N19)2+ (σsh N19)2 (4.33) where σmle N19 , is the statistical uncertainty obtained with the maximum likelihood method employed for the N19 reconstruction. σθ N19 accounts for the uncertainty on the zenith angle reconstruction which can be estimated using the correction factor described in 128 4 Impact point and shower size reconstruction Figure 4.19: ESD energy resolution. Figure 4.20: Orthogonal pull for the selected hybrid events. 4.6 Energy calibration 129 As a further check, the calibration procedure has been performed for three eyes separately. Loma Amarilla has not been used because it has not enough statistics (only 9 events). Results are shown in g. 4.21 and given in table 4.6.2 showing a good agreement within the statistical uctuations. Thus, no systematic on the calibration curve is introduced when using events detected with dierent eyes of the uorescence detector. Figure 4.21: Calibration curves by FD eye. Leones (red), Coihueco (blue) and Morados (pink) calibration curves show a good agreement between them within the statistical uncertainties that are represented by shaded areas around the mean values of the calibration. 130 4 Impact point and shower size reconstruction Eye a σab σb N Leones -0.87 0.05 1.01 0.05 39 Coihueco -0.83 0.06 1.00 0.05 35 Los Morados -0.79 0.05 0.96 0.05 34 Table 4.3: Calibration parameters by eye. To study a possible bias introduced by the low energy event rejection, the calibration procedure has been performed for ve dierent low energy thresholds for rejection: 1018.7eV , 1018.9eV , 1019.1eV , 1019.3eV and 1019.5eV . As it is shown in g. 4.22 all of them are compatible within the uncertainties. The resulting parameters as well as the number of events in the ts are shown in table 4.6.2. log10(Ecut FD/eV )a σab σb N 18.7 -0.82 0.03 0.98 0.03 115 18.9 -0.78 0.04 0.95 0.04 80 19.1 -0.75 0.08 0.93 0.06 49 19.3 -0.87 0.11 1.01 0.08 39 19.5 -0.81 0.31 0.98 0.20 19 Table 4.4: Calibration parameters for dierent low energy cuts. No systematic is seen due to the uorescence eye or to the low energy cut selection. All the values of the parameters of the calibration curve are compatible between them within the uncertainties. 4.6 Energy calibration 131 Figure 4.22: Calibration curves for dierent low energy cuts. The reference one with the cut at 1018.7eV (red) is compared with the others. Top left Ecut = 1018.9eV (dark blue), top right Ecut = 1019.1eV (green), bottom left Ecut = 1019.3eV (yellow) and bottom right Ecut = 1019.5eV (light blue). The calibration curves show a good agreement between them within the statistical uncertainties which are represented by shaded areas around the mean values of the calibration. 132 4 Impact point and shower size reconstruction Bibliography [1] M. Ave, R. A. Vázquez & E. Zas, Astropart. Phys. 14 (2000) 91. [2] M. Ave, R. A. Vázquez, E. Zas, J. A. Hinton & A. A. Watson, Astropart. Phys. 14 (2000) 109. [3] M. Ave, J. A. Hinton, R. A. Vázquez & E. Zas, Phys. Rev. D 65 (2002) 063007. [4] M. Ave, J. A. Hinton, R. A. Vázquez, A. A. Watson & E. Zas, Phys. Rev. D 67 (2003) 043005. [5] M. Ave et al. , Nucl. Instrum. Meth. A 578 (2007) 180. [6] R. A. Fisher, Philos. Trans. Roy. Soc. London Ser. A 222 , (1922) 309. [7] M. Aglietta, I. De Mitri, S. Maglio, S. Maldera, I. C. Maris, D. Martello, G. Navarra & M. Roth, Recovery of saturated signals of the surface detector, Pierre Auger Collaboration Internal Document, GAP Note 2008-030. [8] D. Veberic & M. Roth, SD Reconstruction Oine Reference Manual, Pierre Auger Collaboration Internal Document, GAP Note 2005-035. [9] P. Billoir, O. Deligny & A. Letessier-Selvon, A complete prcedure for the reconstruction of inclined air showers, Pierre Auger Collaboration Internal Document, GAP Note 2003-003. [10] H. Dembinski, T. Hebbeker & M. Leuthold, A comparison of Monte-Carlo generated Muon Maps with near horizontal SD showers, Pierre Auger Collaboration Internal Document, GAP Note 2007-124. 134 BIBLIOGRAPHY [11] G. Rodríguez-Fernández PhD Thesis. Horizontal air showers at the Pierre Auger Observatory, (2007). [12] L. Cazón, R. A. Vázquez, A. A. Watson & E. Zas, Astropart. Phys 21 (2004) 71. [13] S. J. Sciutto Proc. 26th International Cosmic Ray Conference, Salt Lake City-USA, (1999). S. J. Sciutto astro-ph 9911331. [14] D. Heck et al. , The CORSIKA Air Shower Simulation Program, Report FZKA 6019 (1998). [15] http://geant4.web.cern.ch/geant4 [16] I. Valiño PhD Thesis. Detection of horizontal air showers and neutrino induced showers with the Pierre Auger Observatory, Pierre Auger Collaboration Internal Document, GAP Note 2008-024. [17] C. Pryke, Asymmetry of Air Showers at Ground Level, Pierre Auger Collaboration Internal Document, GAP Note 1998-034. X. Bertou & Pierre Billoir, On the origin of the asymmetry of ground densities in inclined showers, Pierre Auger Collaboration Internal Document, GAP Note 2000-017. M. T. Dova, Asymmetries Observed in Giant Air Showers Using Water Cherenkov Detectors, Proc. 28th International Cosmic Ray Conference, Tsukuba-Japan, 8 (2003) 34. [18] D. M. Edge et al. , J. Phys. A 6 (1973) 1612. [19] J. Alvarez-Muñiz et al. , Phys. Rev. D, 66 (2002) 123004. [20] M. Ave et al. , Astropart. Phys. 19 (2003) 47. [21] Vertical calibration [22] I. C. Maris PhD Thesis, Measurement of the Ultra High Energy Cosmic Ray Flux using Data of the Pierre Auger Observatory, Pierre Auger Collaboration Internal Document, GAP Note 2008-026. BIBLIOGRAPHY 135 [23] C. di Giulio et al. , Study of systematic uncertainties on the S38 energy calibration 136 BIBLIOGRAPHY Chapter 5 Energy spectrum of UHECR In this section, the energy spectrum of Ultra High Energy Cosmic Rays is derived using the surface detector data with a reconstructed zenith angle between 60◦ and 80◦ . The energy spectrum of the cosmic rays is the dierential ux of cosmic rays at a certain energy E : J(E) . It can be obtained dividing the number of cosmic ray events detected, N(E) , in an energy interval ∆E by the exposure of the detector at that energy, η(E) . J(E)∆E=N(E) η(E) (5.1) As the T5 trigger condition discussed in section 2.3.2 is applied to the data to guarantee a high quality of the reconstruction parameters, the total detection area of the surface detector is constrained. This constraint allows us to calculate the active detection area as a multiple of the area of an elementary array cell, acell . Above a given energy threshold, that will be chosen to select events triggering the surface detector independently of their incoming direction or of their energy, the aperture of this elementary area, Acell , becomes purely geometrical. Thus, counting the number of active cells, Ncell , in a given period, ∆t , we can obtain the exposure of the whole detector as η=Ncell Acell ∆t . 5.1 Saturation energy of the trigger eciency The trigger eciency, dened as the probability that an extensive air shower induces a given trigger level, depends on some parameters of the shower such as energy and 144 5 Energy spectrum of UHECR The uncertainty in the acceptance has several contributions mainly from the trigger eciency dependence, the detector response, the non regularity of the array or the rejection of unstable periods. All these contributions have a global eect that has be quantied below ∼1% in the analysis of events below 60◦ [6] being the contribution of the non regularity of the array ∼0.05% and that due to the rejection of bad periods of data taking ∼0.5% . We shall use this value in our results. Besides, the reconstruction program fails to reconstruct less than 0.5% of the T4 events. These problems may be due to ineciencies of the selection or reconstruction algorithms. In this work a 1% systematic uncertainty will be considered conservatively to account for this eciency loss. Uncertainties in energy The uncertainty in the energy derived from the surface detector measurement ESD , has two contributions. The rst one arises from the calibration t and the second one from the uncertainty in N19 . Besides the uncertainties in N19 due to the maximum likelihood method, to the propagation of the uncertainty in the zenith angle reconstruction and to the shower to shower uctuations that have been accounted for in the calibration procedure on an event by event basis, other sources of uncertainties have been identied. Those include the muon maps or the electromagnetic signal parameterization. Calibration uncertainty The uncertainty due to the calibration is obtained propagating uncertainties directly in eq. 4.32: σE(cal) E= ln(10) qσ2 A+σ2 Blog2(N19) + 2ρσAσBlog(N19) (5.5) where σA and σB are the uncertainties in the tted parameters and ρ is the correlation coecient of the t. As it is shown in g. 5.5 this uncertainty amounts to less than 4% at 1020 eV . 5.4 Flux uncertainties 145 Figure 5.5: Relative uncertainty due to the calibration process. Uncertainty due to muon maps To evaluate the systematic uncertainty due to muon maps we have reconstructed the same set of events with the same reconstruction program but using a dierent set of muon maps. Comparing the N19 value obtained using the USC muon maps [7] with the value obtained using the Aachen muon maps [8], an uncertainty on N19 below 10% as shown in g. 5.6 can be quoted. 146 5 Energy spectrum of UHECR Figure 5.6: N19 uncertainty due to muon maps. Most of this systematic uncertainty is practically reabsorbed during the calibration procedure as it is shown in g. 5.7 where the N19 have been converted to energy using for each set their corresponding calibration curve. A systematic uncertainty below 0.5% is observed in the reconstructed energy although the distribution shows non-gaussian tails that should be studied with more detail. The energy spectrum should be insensitive to the choice of muon map after the calibration with hybrid events as it was already mentioned in [8]. However, the fact that the shower size is about 10% greater when measured with respect to one set of maps than when measured with respect to another, may be important to draw any conclusions about the number of muons in data compared to those predicted by simulations. 5.4 Flux uncertainties 147 Figure 5.7: Systematic uncertainty in energy due to muon maps. Uncertainty due to the electromagnetic correction of the signal Proceeding in a similar way to evaluate the uncertainty introduced by the parameterization of the electromagnetic signal, we reconstruct a sample of events increasing or decreasing the electromagnetic contribution within the uncertainty of the parameterization quoted to be 10% in the parameterization of the Monte Carlo simulations that have been used [] The resulting uncertainty is below 5% in N19 as shown in g. 5.8. This uncertainty should be pretty reabsorbed through a calibration procedure and it is thus expected to be small. Some additional studies are needed to evaluate eects of composition, model assumptions and energy. 148 5 Energy spectrum of UHECR Figure 5.8: Systematic uncertainty in N19 due to the parameterization of the electromagnetic signal. Table 5.4 summarizes the systematic uncertainties in the energy spectrum derived with inclined events from the surface detector. Source Uncertainty Saturation of trigger eciency 2% Acceptance 1% Reconstruction algorithms 1% N19 uncertainty <14% Calibration <5% Maps 1% Electromagnetic correction 5% TOTAL <16% Table 5.2: Systematic uncertainties in the cosmic ray ux. Some of these uncertainties are energy dependent so they do not aect the energy spectrum as a global factor but they have to be properly propagated to it. Bibliography [1] N. G. Busca PhD. Thesis, The ultra high energy cosmic ray ux from the southern Pierre Auger observatory data, Pierre Auger Collaboration Internal Document, GAP Note 2006-108. [2] I. C. Maris PhD. Thesis, Measurement of the Ultra High Energy Cosmic Ray Flux using Data of the Pierre Auger Observatory, Pierre Auger Collaboration Internal Document, GAP Note 2008-026. [3] C. Bonifazi & A. Letessier-Selvon, Event selection using the T5 time distribution, Pierre Auger Collaboration Internal Document, GAP Note 2006-042. C. Bonifazi & P. Ghia, Selection of data periods and calculation of the SD geometrical acceptance, Pierre Auger Collaboration Internal Document, GAP Note 2006-101. [4] G. J. Feldman & R. D. Cousins, Phys. Rev. D 57 (1997) 3873. [5] J.Abraham et al. , Phys. Rev. Lett. 101 (2008) 061101. [6] C. Bonifazi & P. Ghia, Selection of data periods and calculation of the SD geometrical acceptance, Pierre Auger Collaboration Internal Document, GAP Note 2006-101. D. Allard, Determination of the aperture of the PAO surface detector, Proc. 29th International Cosmic Ray Conference, Pune-India, 7 (2005) 71. [7] M. Ave, R. A. Vázquez & E. Zas, Astropart. Phys. 14 (2000) 91. 150 BIBLIOGRAPHY [8] H. Dembinski, T. Hebbeker & M. Leuthold, A comparison of Monte-Carlo generated Muon Maps with near horizontal SD showers, Pierre Auger Collaboration Internal Document, GAP Note 2007-124. Chapter 6 Conclusions In this work, the complete reconstruction chain of events of the surface detector of the Pierre Auger Observatory with a zenith angle over 60◦ has been completely described and the the cosmic ray ux measured using those events has been presented. The high statistics of the surface detector allows a good estimate of the energy spectrum of UHECRs above 1018.7eV . It has been obtained using data collected from 1st January 2004 to 31st August 2008 accounting for an exposure of ∼2700 km2sry . Regarding the reconstruction of the incoming direction of the cosmic ray, the time model for the distribution of the arrival time of muons in the shower front has been validated. This model predicts the distributions of the time delay of the particles in the shower front with respect to a plane moving at the speed of light in time with the primary cosmic ray. The mean time delay can be associated to the front curvature and the variance in the arrival time of the rst particle can be also predicted. The time variance predicted by the model has been compared to data making use of doublets. These are made of two tanks placed at practically the same location. The dierence in the residuals of the time t of the two stations has been compared to the predictions obtained from the model. The performance of the time variance has been shown to be in a very good agreement with data. The variance has been shown to agree with data up to a distance of ∼4000 m to the shower core from and up to signals as large as ∼50 V EM without any cuts. The agreement with data is much better than that presented by the model of time variance employed for the analysis of vertical events. In the procedure to reconstruct the shower core position and the shower size, the shower size is estimated as the relative number of muons with respect to a reference 152 6 Conclusions 1019 eV proton shower. This parameter, called N19 , has been shown to have a great correlation with the calorimetric measurement of the energy by the uorescence detector in an analogous way to the calibration procedure for events below 60◦ , playing N19 a role similar to S38 . N19 is used as an estimator of the energy and it is calibrated using hybrid events. The energy of an event collected by the surface detector needs several assumptions about high energy interactions and about the cosmic ray primary. In this work N19 has been related to the energy measurement of the energy for a good quality sample of golden hybrid events to get rid of all of those assumptions. The calibration procedure requires a number of cuts to avoid biases due to threshold eects and noise due to badly reconstructed events. A set of quality cuts has been designed to have the calibration t. They are: • Cherenkov fraction below 50% . • The tank with largest signal must be closer than 750 m to the shower axis. • The reduced χ2 of the Gaisser-Hillas t must be below 4 and the χ2 of a linear t to the prole exceeds the χ2 of the Gaisser-Hillas t by at least 4 . • The shower maximum is required to be within the eld of view and determined with an uncertainty below 50 g/cm2 . • The uncertainty in the reconstructed energy must be below 40% . for the uorescence detector reconstruction and: • The relative uncertainty in N19 must be below 40% . • The inclined T5 quality trigger condition must be satised. for the surface detector reconstruction. This correlation has been used to obtain the relation between N19 and the calorimetric energy measured with the uorescence detector using a linear t in log scale ( log10(E[EeV ]) = a+blog10(N19) ). This t has been performed on a subsample of 153 almost 300 events assigning to each one the statistical uncertainties on both EFD and N19 . The values of the tted parameters are: a= 0.84 ±0.01(stat) b= 1.01 ±0.03(stat) This correlation has been stable since the rst presentation of this method when it was obtained with only 17 events. A study has been made of systematic eects of the cuts in the calibration curve. Using this parameterization, the so-called calibration curve, the expected calorimetric energy of any event registered by the surface detector can be obtained. The energy resolution of the method has been estimated to be 21% . A global 22% of systematic uncertainty in the absolute scale must be kept in mind in spite of not being included in this analyses. It comes from the systematics in the FD energy measurement. The energy saturation of the trigger eciency that allows a geometrical calculation of the aperture has been evaluated from hybrid events. The threshold energy for full eciency of the surface detector has been estimated in 1018.7eV . The rst energy spectrum using inclined events has been presented using more than 4 years of data with an incoming zenith angle between 60◦ and 80◦ . As a result, an evidence of a ux suppression around 5·1019 eV has been found. The ux of cosmic rays has found to be not well described using an overall power-law. The rst part of the spectrum ( 18.7< log(E[eV ]) <19.7 ) has been tted to a power-law ( J∝E−γ ) obtaining a spectral index of γ= 2.70 ±0.02 . The number of events expected from the extrapolation of the power-law function J∝E−2.7 at higher energies has turned out to be 53 ±2 for log(E[eV ]) >19.7 whereas the observed number has been 34 events. This ux suppression, combined with the correlation of arrival directions of the highest energy events with nearby AGN, supports the existence of the GZK cut-o but this claim will have to be conrmed with mass composition studies, a larger statistics and smaller systematic uncertainties.