The cosmic Ray spectrum above 0.1 EEV measured by the telescope array and tale fluorescence telescopes

The Telescope Array (TA), deployed in central Utah, is the largest cosmic ray detector in the Northern hemisphere. It was initially designed to study Ultra-High Energy Cosmic Rays (UHECRs) with energies above 1019 eV. TA has added an extension known as the Telescope Array Low Energy Extension (TALE) to lower the experiment's energy threshold. The extension includes the addition of high elevation angle telescopes to one of the existing TA telescope stations (Middle Drum), giving it a field of view of 3-59◦ in elevation and 114◦ in azimuth. This enables the site to see more of the shower development for nearby, lower energy events and allows us to study cosmic rays in 1015.3-1018.3 eV energy range. The observatory now consists of 48 telescopes at three sites overlooking an array of 610 scintillation counters. I have used the data collected between 1 June 2014 and 31 January 2015 to make a monocular measurement of the cosmic ray energy spectrum. I have combined the fluorescence data from the original high energy TA telescopes with that of the TALE low energy extension telescopes, focusing in the energy range of 1017.2 - 1019.0 eV. This is a joint measurement using two sets of telescopes. I have found the energy spectrum has a spectral slope of -3.3±0.06 in this energy range. This is in good agreement with measurements made using other TA techniques as well as other external experiments in this energy range.

THE COSMIC RAY SPECTRUM ABOVE 0.1 EEV MEASURED BY THE TELESCOPE ARRAY AND TALE FLUORESCENCE TELESCOPES by Jihee Kim A dissertation submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Physics Department of Physics and Astronomy The University of Utah December 2018 Copyright © Jihee Kim 2018 All Rights Reserved ❚❤❡ ❯♥✐✈❡rs✐t② ♦❢ ❯t❛❤ ●r❛❞✉❛t❡ ❙❝❤♦♦❧ STATEMENT OF DISSERTATION APPROVAL This dissertation of Jihee Kim has been approved by the following supervisory committee members: Gordon B. Thomson , Chair September 21, 2018 Date Approved Orest G. Symko , Member September 21, 2018 Date Approved , Member Douglas R. Bergman September 21, 2018 Date Approved , Member Pearl E. Sandick September 21, 2018 Date Approved Cynthia M. Furse , Member September 21, 2018 Date Approved and by of Peter Trapa the Department of Physics and Astronomy and by David B. Kieda of the Graduate School. , Chair ABSTRACT The Telescope Array (TA), deployed in central Utah, is the largest cosmic ray detector in the Northern hemisphere. It was initially designed to study Ultra-High Energy Cosmic Rays (UHECRs) with energies above 1019 eV. TA has added an extension known as the Telescope Array Low Energy Extension (TALE) to lower the experiment’s energy threshold. The extension includes the addition of high elevation angle telescopes to one of the existing TA telescope stations (Middle Drum), giving it a field of view of 3-59◦ in elevation and 114◦ in azimuth. This enables the site to see more of the shower development for nearby, lower energy events and allows us to study cosmic rays in 1015.3 -1018.3 eV energy range. The observatory now consists of 48 telescopes at three sites overlooking an array of 610 scintillation counters. I have used the data collected between 1 June 2014 and 31 January 2015 to make a monocular measurement of the cosmic ray energy spectrum. I have combined the fluorescence data from the original high energy TA telescopes with that of the TALE low energy extension telescopes, focusing in the energy range of 1017.2 - 1019.0 eV. This is a joint measurement using two sets of telescopes. I have found the energy spectrum has a spectral slope of -3.3±0.06 in this energy range. This is in good agreement with measurements made using other TA techniques as well as other external experiments in this energy range. CONTENTS ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix ACKNOWLEDGMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x CHAPTERS 1. 2. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 1.2 1.3 1.4 1 5 6 8 Cosmic Ray Energy Spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Cosmic Ray Composition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Cosmic Ray Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Outline of the Dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . THE TELESCOPE ARRAY EXPERIMENT . . . . . . . . . . . . . . . . . . . . . . . 10 2.1 Overview of TA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Middle Drum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Mirrors and PMTs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Mirror reflectivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.3 Electronics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.4 Trigger . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.5 Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 TALE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Mirrors and PMTs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Mirror reflectivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.3 Electronics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.4 Trigger . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.5 Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3. 10 13 13 19 20 24 24 27 27 33 33 36 36 EXTENSIVE AIR SHOWER PHYSICS . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.1 Electromagnetic Showers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.2 Hadronic Showers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4. MONTE CARLO SIMULATION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.1 Monte Carlo Program . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1 CORSIKA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.2 Input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.3 Weighting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.4 Light production . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.5 Light propagation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.6 Light attenuation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 50 51 53 55 56 57 4.2 Monte Carlo Event . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Thrown distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.2 Resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.3 Missing energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5. 57 58 61 61 EVENT RECONSTRUCTION AND SELECTION . . . . . . . . . . . . . . . . . . 66 5.1 Data Format . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Reconstruction Programs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.1 Preprocessing pass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.2 Pass 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.3 Pass 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.4 Pass 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.5 Pass 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.6 Pass 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.7 Pass 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.8 Pass 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.9 Postprocessing pass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Event Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Example Event . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 68 69 71 71 72 74 76 77 79 80 83 85 6. DATA/MONTE CARLO COMPARISONS . . . . . . . . . . . . . . . . . . . . . . . . . 89 7. ENERGY SPECTRUM MEASUREMENT . . . . . . . . . . . . . . . . . . . . . . . . . 103 7.1 Spectrum Calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1.1 Aperture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1.2 Detector ontime . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1.3 Exposure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Energy Spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8. 103 105 105 107 107 CONCLUSIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 v LIST OF FIGURES 1.1 The flux of cosmic ray particles measured by various experiments. . . . . . . . . . 2 1.2 The flux of cosmic ray particles measured by various experiments multiplied by E 2.6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 hXmax i observed by Auger, HiRes, and HiRes/MIA experiments as a function of primary cosmic ray energy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.4 The significance map for the direction of the UHECRs observed by the TA surface detectors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1 Map depicting the original, high energy, portion of the Telescope Array experiment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2 Photos of fluorescence detector stations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3 Photos of an SD and SD assembling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.4 A pair of MD telescopes covering viewing angle of 3◦ -31◦ in elevation. . . . . . . 16 2.5 The UV filter transmissivity as a function of wavelength. . . . . . . . . . . . . . . . . . 17 2.6 Photos of a cluster of PMTs with the UV filter and the face of a Philips PMT. 18 2.7 The image of the TA MD electronics rack. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.8 The diagram of the MD data acquisition crate. . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.9 The schematic outline of the MD electronics. . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.10 A model of MD trigger patterns. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.11 Photo of the RXF module. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.12 Photo of RXF placed at the center of the mirror. . . . . . . . . . . . . . . . . . . . . . . . 26 2.13 Map of Telescope Array and TALE detectors. . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.14 A pair of TALE telescopes with viewing angle of 30◦ -59◦ in elevation. . . . . . . . 30 2.15 Photo of the MD (left) and TALE (right) telescope buildings with opened shelter doors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.16 The field of view in the sky of the MD and TALE telescopes. . . . . . . . . . . . . . 32 2.17 Photo of front side of TALE electronics rack. . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.18 Photo of back side of TALE electronics rack. . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.19 The diagram of TALE triggering. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.20 An example of TALE triggering. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.21 Photo of UVLED module placed at the center of a mirror. . . . . . . . . . . . . . . . 38 3.1 Schematic view of an electromagnetic cascade induced by a cosmic ray photon. 40 3.2 Schematic view of a hadronic cascade induced by a cosmic ray proton. . . . . . . 43 4.1 An example of a Monte Carlo shower. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.2 The energy spectrum measured by TA in 2015. . . . . . . . . . . . . . . . . . . . . . . . . 53 4.3 The mean Xmax from HiRes/MIA and HiRes measurements shown with inferred chemical composition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.4 Thrown energy (E) distribution and thrown impact parameter (Rp ) distribution. 59 4.5 Thrown zenith angle (θ) distribution and thrown azimuthal angle (φ) distribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.6 The energy resolution histogram and the energy resolution as a function of energy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.7 The Rp resolution histogram and the Rp resolution as a function of energy. . . 63 4.8 The ψ resolution histogram and the ψ resolution as a function of energy. . . . . 64 5.1 Examples of a track-like event and a random noise event. . . . . . . . . . . . . . . . . 73 5.2 Shower detector plane (SDP) with the shower axis and the detector location. . 75 5.3 An event display shown with the fitted shower detector plane indicated by the black line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.4 Geometry of an extensive air shower. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.5 An example of time versus angle fit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.6 An example of shower profile fit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.7 The λ variation as a function of energy and the X0 variation as a function of energy from the CORSIKA simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.8 The variation in X1 as a function of energy from the CORSIKA simulation. . . . . 82 5.9 CORSIKA simulation of the fractional calorimetric energy as a function of the total shower energy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.10 Event display of an example event observed by both the Middle Drum and TALE telescopes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.11 The time versus angle fit for the same event as Figure 5.10. . . . . . . . . . . . . . . . 87 5.12 The profile fit for the same event as shown in Figure 5.10. . . . . . . . . . . . . . . . . 88 6.1 The Data/Monte Carlo comparison for the event track length distribution. . . . 90 6.2 The Data/Monte Carlo comparison for the event crossing time distribution. . . 91 6.3 The Data/Monte Carlo comparison for the event zenith angle (θ) distribution. 92 6.4 The Data/Monte Carlo comparison for the event ψ angle distribution. . . . . . . 93 6.5 The Data/Monte Carlo comparison for the event azimuthal angle (φ) distribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 6.6 The Data/Monte Carlo comparison for the event impact parameter, Rp , distribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 vii 6.7 The Data/Monte Carlo comparison for the number of good tubes per degree of track length. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 6.8 The Data/Monte Carlo comparison for the number of participating rings in an event. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 6.9 The Data/Monte Carlo comparison for total energy distribution. . . . . . . . . . . . 99 6.10 The Data/Monte Carlo comparison for the shower detector plane angle. . . . . . 100 6.11 The Data/Monte Carlo comparison for Xcore position. . . . . . . . . . . . . . . . . . . . 101 6.12 The Data/Monte Carlo comparison for Ycore position. . . . . . . . . . . . . . . . . . . . 102 7.1 The histogram of the number of data events. . . . . . . . . . . . . . . . . . . . . . . . . . . 104 7.2 The calculated aperture as a function of the energy. . . . . . . . . . . . . . . . . . . . . . 106 7.3 The exposure as a function of energy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 7.4 The flux J(E) of cosmic ray measured by the Middle Drum and TALE stations.109 7.5 The energy spectrum measured by the Middle Drum and TALE stations. . . . . 110 7.6 The energy spectrum measured by the Middle Drum and TALE stations with a single slope of the power law. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 7.7 The energy spectrum measured by the Middle Drum and TALE stations comparing to the other recent TA energy spectra. . . . . . . . . . . . . . . . . . . . . . . 111 7.8 The energy spectrum measured by the Middle Drum and TALE stations (black circles) comparing to the other external experiments. . . . . . . . . . . . . . . . . . . . . 112 viii LIST OF TABLES 2.1 The MD telescope serial numbers are listed with corresponding HiRes-1 telescope serial numbers, PMT manufacturer, and ring numbers. . . . . . . . . . . . 15 2.2 The MD mirror serial numbers are listed with corresponding HiRes-2 mirror serial numbers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3 TALE telescope serial numbers are listed with corresponding HiRes-2 telescope serial numbers, PMT manufacturer, and ring number. . . . . . . . . . . . . . . . . . . . 29 2.4 TALE mirror serial numbers are listed with corresponding HiRes-1 or HiRes-2 mirror serial numbers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 5.1 The cut levels are described for three different data sets: TALE, MD, and 4-ring (TALE and MD). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.2 The number of data events passed through each cut level for the three different data sets: TALE, MD, and 4-ring (TALE and MD). . . . . . . . . . . . . . . . . . . . . . 85 7.1 The number of data events per each energy bin. . . . . . . . . . . . . . . . . . . . . . . . . 104 ACKNOWLEDGMENTS I would like to thank my advisor, Professor Gordon Thomson, for his guidance and support. He is a very nice advisor, and without his help, I couldn’t have gotten through the Ph.D. program. At times, he seems like a very tough, strict, and serious boss, but far more often, he was more like a supportive grandfather to me. I would like to thank the Telescope Array Utah group: Dr. John Matthew, for countless help and reviews of my dissertation; Stanton Thomas, for fun and relaxed conversations amidst research discussions; Jeremy Smith, who taught me how to run detectors and made hardware and soldering fun; Dr. Bob Cady, with whom I worked closely during surface detector maintenance and FD shifts; Professor Charlie Jui, for his many helpful comments on my research; Professor Pierre Sokolsky, for his guidance and instruction; Dr. Dmitri Ivanov, for being a good teacher who helped me with analysis, programming, and physics; Dr. Tareq AbuZayyad, who helped me with the UTAFD program and provided lots of helpful comments on my research; Dr. Monica Allen, for her friendship and support; and Jamie Zvirzdin, for her company in Delta. I would like to thank the Telescope Array collaboration members. There are too many people to be able to mention everyone. Nonetheless, I found they are supportive and providing helpful comments on my research. I would like to thank my friends: Paul Bergeron, Flo Doval, Greg Furlich, Lauren Simonsen, Sangita Baniya, and Jessica Frew. Without these friends, it could have been a worse experience. We went through the Ph.D. program together and had lots of fun. I would like to thank Ganymede, for the dog therapy he provided me paid off. He reduced my stress and allowed me to relax while writing my dissertation. And last, but not least, I would like to thank my family in Korea. They have always supported me in whatever I am doing. I could not have done this without my family. CHAPTER 1 INTRODUCTION Cosmic rays are energetic particles that are accelerated by astrophysical sources. They are emitted from the Sun, from sources inside our Galaxy, and even from sources outside our Galaxy in the case of the highest energy cosmic rays. Their energy range is quite broad, with experiments observing cosmic ray energies across 12 orders of magnitude, as depicted in Figure 1.1. Cosmic rays at the high end of the energy spectrum are called Ultra-High Energy Cosmic Rays (UHECRs), and they are the most energetic particles observed in the Universe. In 1991, the High Resolution Fly’s Eye Experiment (HiRes) observed the highest energy cosmic ray ever seen. It induced an extensive air shower with an energy of 3.2 × 1020 eV (Bird et al. 1993). At present, the Conseil Européen pour la Recherche Nucléaire (CERN) particle accelerator, the Large Hadron Collider (LHC), accelerates protons in counter rotating beams to 6.5 TeV (6.5 × 1012 eV) in each beam, making the collisions occur at 13 TeV (13 × 1012 eV), which are the highest energy particles that humans have created. This implies that very violent astrophysical sources exist in the Universe that can accelerate particles far outside the range of human capabilities. However, the sources of UHECRs remain unknown. We speculate that they come from violent extragalactic astrophysical objects. Thus, the aim of UHECR studies is to understand where they come from, what their chemical composition is, how they are accelerated to such high energies, and how they propagate through the intergalactic medium. 1.1 Cosmic Ray Energy Spectrum As Figure 1.1 shows, the flux of cosmic ray particles follows a simple power law, Flux(E) ∝ E −γ , where γ is the power law’s exponent and is called a spectral index. In the energy spectrum of cosmic ray particles, there are spectral features where the spectral indices change. Near E ∼ 1015.5 eV, there is a kink, known as the “knee”, where the spectrum steepens, and another at E ∼ 1018.7 eV, known as the “ankle”, where it becomes less steep. In addition, at E ∼ 1019.8 eV, the spectrum drops off dramatically at what is known as the “GZK cutoff” (Greisen 1966; Zatsepin and Kuz’min 1966). To emphasize these detailed 2 Figure 1.1: The flux of cosmic ray particles measured by various experiments. It nominally follows a simple power law, E −3 , which is shown by a green dashed line. Historically, there are two known spectral features, the “knee” and the “ankle”, which are located at E ≈ 1015.5 eV and E ≈ 1018.7 eV, respectively. Taken from Hanlon (2008). 3 spectral features, the energy spectrum is often multiplied by a power of E, such as E 2.6 as it is in Figure 1.2. In the energy spectrum, the knee feature is marked by a change in spectral index. Below the knee feature, the spectral index is about 2.7 and above the knee feature, the spectral index is about 3.0. Cosmic ray particles in this energy range are considered to be associated with the galactic sources. Supernova Remnants (SNRs) are recognized as candidates for the acceleration of galactic cosmic rays. When a supernova occurs, it ejects many particles during the explosion and it can expand over thousands of years. It provides the energetic shock fronts. As particles repeatedly cross through the shock fronts, they gain energy until they are able to escape. The knee feature can be explained as the maximum energy of galactic cosmic ray protons accelerated by SNRs. This acceleration through the magnetic field depends on the characteristic magnetic rigidity, Rc = E/Ze. Therefore, the maximum energy through the acceleration mechanism in SNRs can be written as p . Emax (Z) = Ze × Rc = Z × Emax (1.1) This is called a rigidity-dependent effect: heavy charged nuclei can achieve Z times higher Figure 1.2: The flux of cosmic ray particles measured by various experiments multiplied by E 2.6 . This emphasizes the spectral features, which include the knee, the ankle, and the GZK cutoff. Taken from Patrignani et al. (2016). 4 energy than a proton when accelerated in a magnetic field. In other words, protons are p accelerated up to their maximum energy (Emax ) and then their acceleration is cut off as they escape the shock front. The next element, helium, continues and is accelerated until its energy reaches twice that of the proton’s maximum energy. This trend continues through galactic iron nuclei. Above 1017 eV, galactic iron dies out and a new feature is introduced, which is called the “2nd knee” or “iron knee”. This region, 1017 -1018 eV, may be the place where the transition from galactic cosmic rays to extragalactic cosmic rays occurs, and the composition changes from heavy nuclei to light nuclei. Cosmic rays with energy greater than 1018 eV are considered as extragalactic particles. The astrophysical sources for these particles should have stronger magnetic fields and larger sizes so that particles can be accelerated over longer periods of time to greater energies. Examples of possible sources include: Supermassive Black Holes (SMBHs), Gamma Ray Bursts (GRBs), star burst galaxies, and BL-Lacs. At these high energies, two additional spectral features are observed. Near E ≈ 1018.7 eV, a dip in the energy spectrum is observed and is called the “ankle”. The reason that there is a dip in the flux of cosmic ray particles is not fully resolved yet, but there are three plausible hypotheses. It may be due to the pair production process of ultra-high energy protons interacting with Cosmic Microwave Background (CMB) radiation, galactic and extragalactic transition, or an additional rigidity-dependent cycle of extragalactic sources. At the highest end of the energy spectrum, the flux of cosmic ray particles drops dramatically due to the GZK cutoff (or “GZK suppression”). When UHECR protons with sufficient energy propagate through the Universe, they interact with the CMB radiation and produce a Delta resonance. In the decay of the Delta baryon, a pion is produced, thus taking energy from the primary particle. These processes can be described by γCMB + p −→ ∆+ −→ p + π 0 (1.2) γCMB + p −→ ∆+ −→ n + π + . (1.3) or As UHECR protons propagate through the interstellar medium and continuously experience the above processes, a theoretical upper limit on the energy of cosmic rays from sources can be estimated as 5-6×1019 eV for a distance range > 50-100 Mpc. The measured flux of cosmic ray particles (Abbasi et al. 2008) supports this theoretical expectation. Moreover, when cosmic rays travel over a distance larger than 50 Mpc and with energies greater than this threshold, they have a small probability of being observed on Earth. This is referred to as the “GZK horizon”. 5 The work in this dissertation is focused on particles in the energy range from 1017 to 1018 eV, where the transition from galactic to extragalactic sources is suspected to occur. Moreover, features in this energy range may imply a difference in chemical composition between galactic and extragalactic particles. The heavy galactic component of cosmic rays appears to be dying off while the light extragalactic component of cosmic rays begins entering the Milky Way. 1.2 Cosmic Ray Composition The chemical composition of low energy cosmic rays can be found by direct measurement with balloon- and satellite-based cosmic ray experiments. Alternatively, for high energy cosmic rays, direct detection becomes impractical due to the low event rate. At these energies, the chemical composition must be inferred from ground-based telescopes, such as fluorescence telescopes that observe the extensive air shower induced by the initial cosmic ray. Ground-based experiments study the depth at which the shower achieves the maximum number of charged particles, Xmax , in an extensive air shower. Xmax provides an indirect evaluation of the primary particle’s chemical composition. This is because the chemical composition determines the average characteristics of the shower development. For a proton cosmic ray, one nucleon enters the Earth’s atmosphere and collides with an atmospheric gas molecule, initiating an air shower. The secondary particles from this interaction will do the same thing and the shower will grow. It will reach the maximum size of shower (Xmax ) when ionization is the dominant mode of energy loss. Afterwards the size of shower will decrease, but the shower will continue to develop. On the other hand, an iron primary cosmic ray can be treated as the superposition of 56 nucleons. Each nucleon has an energy E0 /56 where E0 is the energy of the incident primary. In this scenario, each of the 56 individual showers all start at the same point and are superposed. In other words, when the proton and iron primary cosmic rays with the same energy respectively enter the Earth’s atmosphere, the energy of the incident particle that initiates a shower is 56 times less for the iron primary cosmic ray than that for the proton primary cosmic ray. Thus, a shower created by an iron cosmic ray reaches the maximum size of shower (Xmax ) earlier than a shower created by a proton cosmic ray does. Additionally, it can be explained by the cross-section of the primary cosmic ray. Low mass cosmic rays have smaller cross-sections, allowing them to penetrate further into the atmosphere before having their first interaction. Alternatively, heavier mass cosmic rays have relatively larger cross-sections, and therefore tend to interact 6 with air molecules earlier, causing the shower development to start sooner. This results in a different Xmax , on average, for cosmic rays of different masses. The High Resolution Fly’s Eye Experiment (HiRes) and the Pierre Auger Observatory (PAO) agree on the average Xmax observed for cosmic ray masses between 1017 to 1018 eV, as shown in Figure 1.3. These are both ground-based experiments that use fluorescence detection to observe the full shower development, including the Xmax parameter. However, at 1018 eV or higher energies, the distribution of Xmax begins to diverge. Currently, the Telescope Array (TA) and PAO have a joint working group to examine this divergence in Xmax and to better understand the composition of UHECRs. 1.3 Cosmic Ray Sources As noted in the previous section, we speculate that UHECRs come from violent extragalactic astrophysical objects, though the sources remain unknown. Many experiments (Abraham et al. 2007; Abu-Zayyad et al. 2012b) study the arrival direction of cosmic ray particles and their correlations with possible sources. Possible sources include galaxy clusters, supermassive black holes in active galactic nuclei (AGNs), jets and lobes of active radio galaxies, starburst galaxies, gamma-ray bursts, pulsars, and BL-Lacs. Unfortunately, as charged particles travel through interstellar medium, their pathways are bent due to intergalactic and galactic magnetic fields. In order to estimate UHECRs’ trajectory, many simulations Figure 1.3: hXmax i observed by Auger, HiRes, and HiRes/MIA experiments as a function of primary cosmic ray energy. Taken from Watson and Pierre Auger Collaboration (2008). 7 and models take into account different strengths of magnetic fields at different distance scales for different extragalactic sources. The deflection angle is estimated to be a few degrees in a case of a 60 EeV (60 × 1018 eV) proton from a source 50 Mpc away, assuming models where the intergalactic magnetic fields have a strength of 1 nG (Abbasi et al. 2014). The Telescope Array (TA), using the TA surface detectors, has studied cosmic rays with energies between 1018 to 1018.5 eV in order to search for anisotropy. This study shows that the data are isotropic, with no excess along the galactic plane (Abbasi et al. 2017). This implies that cosmic rays with energies between 1018 to 1018.5 eV are from extragalactic sources. However, using cosmic rays with energy above 5.7×1019 eV, we found a local cluster of high energy events, which we call “the hotspot”, located near Ursa Major and ∼19◦ from the supergalactic plane (Abbasi et al. 2014), as shown in Figure 1.4. Though there are still no specific identified sources behind the hotspot, this evidence for an anisotropy indicates some localized sources exist. Figure 1.4: The significance map for the direction of the UHECRs observed by the TA surface detectors. UHECR sky maps in equatorial coordinates is shown. (a) The points show the direction of the UHECRs with E > 57 EeV (57 × 1018 eV) observed by the TA surface detectors. Note the the presence and location of the galactic and supergalactic planes (GP and SGP) as well as the galactic center (GC) and anti-galactic center (Anti-GC); (b) Color contours show the number of observed cosmic-ray events summed over a 20◦ radius circle; (c) Number of expected background events from the geometrical exposure summed over a 20◦ radius circle; (d) Significance map is calculated from (b) and (c). The “hot spot” that is indicated with red colors is located near Ursa Major just off of the supergalactic plane. Taken from Abbasi et al. (2014). 8 1.4 Outline of the Dissertation The measurement of the cosmic ray energy spectrum may give a hint of where UHECRs come from. This is because spectral features may indicate the energy distributions of particles at astrophysical sources and their distribution depends on the surrounding environment near sources. However, propagation effects may change spectral features. The energy spectrum is a combination of these effects. The cosmic ray energy spectrum can be calculated according to how many particles are observed per unit area, per solid angle, per unit time, and per unit energy. This can be written as N (Ei ) J(Ei ) = , (1.4) AΩ(Ei ) × T × ∆Ei where J(Ei ) is the flux of cosmic rays, N (Ei ) is the number of observed cosmic rays, AΩ(Ei ) is the energy-dependent aperture calculated from Monte Carlo (MC) simulation, T is the detector ontime, and ∆Ei is the width of energy bin for the i th energy bin, Ei . As seen in Figure 1.2, the energy spectrum is quite steep. For each order of magnitude increase in energy, the flux drops by a factor of about 1000. At the energy of the knee, the event rate is about 1 particle per m2 per year, while at the ankle, the event rate is about 1 particle per km2 per year. Because the event rate is so low, we employ an indirect measurement rather than a direct measurement to detect high energy cosmic rays. The scale of the experiment must be huge, necessitating a ground-based design observing cosmic rays via the extensive air showers they induce when they interact with the atmosphere. In Chapter 2, a general description of the Telescope Array (TA) and the Telescope Array Low Energy Extension (TALE) experiments is given. Additionally, how detectors trigger on cosmic ray events and how the detectors are calibrated are introduced. In Chapter 3, Extensive Air Showers (EASs) in the atmosphere using simplified models of particle interactions are introduced, and weaknesses of such models are discussed. Some important aspects to understand are properties of EASs, and how to associate them with observable variables using a fluorescence detecter are described. In Chapter 4, we discuss the Monte Carlo simulations used. Input parameters, assumptions of energy spectrum and composition, production and propagation of fluorescence and cherenkov lights, thrown parameters’ distributions, resolutions, and the missing energy correction are presented. In Chapter 5, the format of data and Monte Carlo events and the reconstruction chain are introduced. We describe how events are selected and show an example event. In Chapter 6, to verify that we understand our detector using Monte Carlo simulation, data and Monte Carlo event distributions are compared for various parameters. We assess 9 the extent to which the data and Monte Carlo agree. In Chapter 7, the calculation of the cosmic ray flux measured by the main TA Middle Drum and TALE extension fluorescence telescopes is presented. We also discuss comparisons to the energy spectrum measured by other cosmic ray experiments. In Chapter 8, we summarize the result of energy spectrum measurement, discuss the scientific importance of this result, and the direction of the future work. CHAPTER 2 THE TELESCOPE ARRAY EXPERIMENT As described in the previous chapter, observing UHECRs requires the scale of the experiment to be huge. Therefore, the ground-based experiments are designed to observe cosmic rays via the extensive air shower they induce in the atmosphere. The Telescope Array (TA) experiment is one such ground-based cosmic ray experiment. The TA experiment is the largest cosmic ray observatory in the Northern hemisphere and is deployed in Millard Country, Utah. It observes UHECRs with energies greater than 1018 eV, using both fluorescence telescopes and scintillator surface detectors. The Telescope Array was built and is operated by an international collaboration of universities and institutions in the United States, Japan, South Korea, Russia, and Belgium. This chapter focuses on the general description of the TA and its low energy extension on the Telescope Array Low Energy Extension (TALE) experiments. We will discuss the fluorescence detectors in detail as the focus of this dissertation is to measure the cosmic ray energy spectrum using the fluorescence telescopes in monocular mode. The scintillation surface detectors will be discussed in less detail. A map of the detectors is shown in Figure 2.1. 2.1 Overview of TA The TA is referred to as a hybrid detector because it uses both fluorescence telescopes and scintillator surface detectors. Using the fluorescence telescopes, we measure the scintillation light produced by air molecules during the development of cosmic ray induced air showers. Three telescope stations overlook the scintillator array from the periphery of the array. The Middle Drum (MD) fluorescence station (Abu-Zayyad et al. 2012a) is the northern most point of the TA while the Black Rock (BR) and Long Ridge (LR) fluorescence stations (Tameda et al. 2009) are located in the southeast and the southwest, respectively. The MD fluorescence station, which is shown in Figure 2.2a, consists of 14 reutilized telescopes from the High Resolution Fly’s Eye (HiRes-1) telescopes. The BR (Figure 2.2b) and LR fluorescence stations each consists of 12 new telescopes. 11 Figure 2.1: Map depicting the original, high energy, portion of the Telescope Array experiment. Locations of the scintillator Surface Detectors (SD) are shown as red diamonds and the locations of the three Fluorescence Detector (FD) stations are indicated by blue hexagons. The orange square indicates the location of the Central Laser Facility (CLF). The location of the Telescope Array is in the remote area west of Delta, Utah and is about a 2 21 hour drive south of Salt Lake City. Reprinted with permission from R. Cady. 12 (a) (b) Figure 2.2: Photos of fluorescence detector stations. (a) Photo of the Middle Drum fluorescence station. There are seven bays, each with two telescopes that view the southeastern sky from 3◦ -31◦ in elevation. Telescopes are numbered from the west end of the building. Reprinted with permission from J. N. Matthews. (b) Photo of the Black Rock Mesa fluorescence station. There are three bays, each with four telescopes that view the northwestern sky from 3◦ -31◦ in elevation. The Long Ridge fluorescence station is similar to the Black Rock Mesa station, however, its telescopes view the northeastern sky. Reprinted with permission from S. Udo. 13 An array of 507 scintillator surface detectors (Abu-Zayyad et al. 2012c), shown in Figure 2.3a, are deployed in a square grid with 1.2 km spacing, covering an area of about 700 km2 . The scintillator detector array is used to sample the density of shower particles arriving at the ground. Each SD consists of 2 layers of scintillators, shown in Figure 2.3b, 1.25 cm thick and 3 m2 area. Optical fibers are placed into grooves in the scintillator plastic and run to PhotoMultiplier Tubes (PMTs). The Central Laser Facility (CLF) is located equidistant from each fluorescence station at 20.85 km. The CLF provides a known light source, allowing for atmospheric calibration. Every half hour, it fires a measured 355 nm laser pulse at 10 Hz vertically into the sky. By measuring the laser light observed at each telescope, we can determine the light loss due to aerosol scattering in the atmosphere. 2.2 Middle Drum The Middle Drum (MD) fluorescence detector was refurbished from the HiRes-1 telescopes, indicated in Table 2.1, and the HiRes-2 mirrors, indicated in Table 2.2. It is located in the northern end of the TA array and is operated mainly by the University of Utah group with help from the Russian groups. 2.2.1 Mirrors and PMTs The MD fluorescence station is instrumented with 14 telescopes each consisting of a spherical mirror and a camera. The mirror is composed of four segments arranged in a cloverleaf shape with a total area of 5.2 m2 . The mirror segments are adjustable and are aligned to focus light onto the camera. A camera is composed of a 16 × 16 array of hexagonal PMTs. Each PMT is 40 mm flat to flat. They are in hexagonal close pack or honeycomb geometry. Each of 256 PMTs views 0.98◦ of sky and its effective area is 1079 mm2 . Each telescope has a field of view of about 14◦ in elevation and 16◦ in azimuth. The MD telescopes view 3◦ -31◦ in elevation and about 115◦ in azimuth. Seven of 14 Middle Drum telescopes view 3◦ -17◦ in elevation, while the remaining seven view 17◦ -31◦ in elevation. We define a term “ring” to represent a set of telescopes that share the same elevation view. The MD detector has two rings: ring 1 views from 3◦ -17◦ in elevation while ring 2 views 17◦ -31◦ in elevation. Figure 2.4 shows a pair of MD telescopes, one in ring 1 and one in ring 2. Together they view 3◦ -31◦ in elevation. Fluorescence light emitted by the air shower is collected by the mirrors and focused onto a cluster of PMTs that are the pixels of the camera. In front of the PMTs is a UV band-pass filter to improve the signal to noise ratio. This UV filter transmits in the 14 (a) (b) Figure 2.3: Photos of an SD and SD assembling. (a) Photo of an SD used by the TA experiment (here, deployed for the TALE in-fill array). SDs are solar powered, and data are replayed to a communication tower by a wireless antenna. Reprinted with permission from J. N. Matthews. (b) Photo of SD assembling at the University of Tokyo’s Akeno Observatory. JiHee Kim (right) and Lito-san (left) installing fibers into a scintillator detector. The fibers collect scintillation light from the scintillator pads and deliver it to photomultiplier tubes so that it can be measured. Reprinted with permission from J. N. Matthews. 15 Table 2.1: The MD telescope serial numbers are listed with corresponding HiRes-1 telescope serial numbers, PMT manufacturer, and ring numbers. These are listed for historical tracking of the optical components. MD telescope # 01 02 03 04 05 06 07 08 09 10 11 12 13 14 HiRes-1 telescope # 01 06 03 04 09 11 07 08 02 10 05 12 13 14 PMT EMI EMI EMI EMI EMI EMI Philips Philips Philips Philips Philips Philips Philips Philips ring # 1 2 1 2 1 2 1 2 1 2 1 2 1 2 Table 2.2: The MD mirror serial numbers are listed with corresponding HiRes-2 mirror serial numbers. These are listed for historical tracking of the optical components. MD mirror # 01 02 03 04 05 06 07 08 09 10 11 12 13 14 HiRes-2 mirror # 22 17 21 18 36 34 35 33 32 29 31 30 19 15 16 Figure 2.4: A pair of MD telescopes covering viewing angle of 3◦ -31◦ in elevation. Reprinted with permission from J. N. Matthews. 17 300-400 nm wavelength range with high efficiency. This corresponds to fluorescence light produced by the air showers while it filters out most visible light. Figure 2.5 shows the UV filter transmissivity as a function of wavelength and Figure 2.6a shows a cluster of PMTs with the UV filter. There are two different kinds of PMTs employed in the MD station and these PMT cameras also are reutilized from the HiRes-1 experiment. PMTs with EMI 9974KAFL model are used in cameras 1-6 and PMTs with Philips/Photonis XP3062/FL model, shown in Figure 2.6b, are used in cameras 7-14. The EMI tubes have a thicker, more spherical front face and lower mean quantum efficiency compared to the Philips tubes. These effects are taken into account in the calibration processes. Each PMT is optimized to collect UV light and is provided with its own high voltage setting to provide uniform gain. When it comes to detecting light, efficiency is an important factor to be considered. There are two kinds of efficiencies; Quantum Efficiency (QE) and Collection Efficiency (CE). The quantum efficiency is expressed as a fraction and is defined as the number of photoelectrons emitted from the photocathode per incident photon. We have measured quantum efficiency curves for a sample of PMTs, which were provided by Figure 2.5: The UV filter transmissivity as a function of wavelength. This UV filter transmits in the 300-400 nm wavelength range with high efficiency while it filters out most visible light. Taken from Zundel (2016). 18 (a) (b) Figure 2.6: Photos of a cluster of PMTs with the UV filter and the face of a Philips PMT. (a) Photo of JiHee Kim standing next to a cluster of PMTs with the UV filter lowered to show the PMTs in the camera. Reprinted with permission from J. N. Matthews. (b) Photo of the face of a Philips PMT. Reprinted with permission from J. N. Matthews. 19 Philips. We use a mean of these curves in the analysis. The value is 0.278 at a wavelength of 355 nm. The collection efficiency is expressed as a fraction and is defined as the probability that photoelectrons produced at the photocathode will be multiplied effectively at the successive dynode stages. Some photoelectrons reflect off the surface of the dynode or are lost or absorbed, and these have no contribution to multiplication at later dynodes. We set CE to be 0.9, which is provided from PMT manufacturers. The number of photons hitting the PMTs determines signal size. Photons are converted into photoelectrons with the QE of photocathode. These photoelectrons are increased via a high voltage dynode chain and the PMT anode collects the photoelectrons from the multiplier and supplies the output signal as in current. The pre-amplifiers attached to each PMT convert from the current into a voltage. This PMT signal is then passed via twisted-pair cable to the rack of data acquisition electronics. It is converted to the digital from the charge collected from the PMT by the Charge to Digital Converter (QDC) and the pulse is integrated if trigger conditions are met. To accurately measure the amount of the light, we should consider two additional factors; obscuration of the mirror and atmospheric conditions including aerosols. In accounting for the obscuration of the mirrors by the PMT cluster and camera stand, the total effective area of the mirror is reduced to 3.72 m2 . By the ray-tracing, we include fluctuations in the number of photoelectrons reaching the detector, PMT tube profiles, and lateral spread of the electrons in the shower. The weather conditions affect the amount of light transmitted via the atmosphere. Hourly weather observations are made and recorded as the weather code. The weather code consists of 7 digits and records horizon clouds, overhead clouds, cloud thickness, and haziness. By dividing the sky by 20◦ from the ground, we define horizon and overhead sections. Determining 20◦ is by using hands and fingers as a measuring tool. When you hold your hand at arm’s length, you can estimate angles like this: Clench your fist at arms length, and hold it with the back of your hand facing you. The width of your fist represents about 10◦ . Then we examine existence of clouds on two divided sections, cloud thickness, and haziness at the MD station. The MD detector is designed for observing high energy cosmic rays, thus horizon and overhead is an important factor to be considered. Since higher energy events are brighter, the telescopes can see these showers from greater distances. 2.2.2 Mirror reflectivity The mirror reflectivity is an important factor determining the fraction of photons from an air shower reaching the PMTs. Dust from the desert accumulates on mirrors and reduces 20 the fraction of photons reflected by the mirrors. We wash the mirrors at all three sites once a year. This usually happens at the beginning of the summer right around the time of a TA collaboration meeting at the University of Utah. The mirror reflectivity is measured before and after washing using a reflectometer. This is carried out at four points per segment and a total of 16 points per mirror measurements are taken between 360-740 nm in 10 nm steps for a total of 39 different wavelengths. The average mirror reflectance is the average of these measurement and is defined as a constant value, 80% in the MD analysis. 2.2.3 Electronics Each MD telescope has its own electronics rack for data acquisition, shown in Figure 2.7. The data acquisition rack consists of low and high voltage power supplies for the camera, a high voltage distribution crate to allow adjusting individual PMT gain, low voltage power supplies for the data-acquisition electronics, and a VME crate for the data-acquisition electronics. There are three fans are located between electronics for cooling. The camera’s low and high voltage power supplies are located in the rack. The low voltage supply delivers ±15 VDC to power the pre-amplifiers. The high voltage supply delivers about 1400 V to the high voltage distribution crate. The high voltage distribution crate consists of one HV zener board and eight HV peg boards. The zener board has a chain of zener diodes that divide the 1400 V input down in 20 V steps and distributes the stepped voltages to the eight HV peg boards. These HV peg boards allow one of the stepped voltages to be selected for each PMT to balance the PMT gains in a camera. Each channel HV is then distributed to the PMTs in the camera via HV cables. The HV peg boards also have HV divider resistors to allow read back of each HV channel. The LV power supply for the data acquisition electronics supplies +5 VDC and ±12 VDC to the VME crate. In the VME crate, there are five different types of boards: 16 ommatidial boards (OMBs), a trigger board, a crate central processing unit (CPU) board, a programmable pulse generator (PPG), and a voltage monitoring board (“garbage” board). The diagram of the boards in the VME crate is shown on Figure 2.8. Each ommatidial board (the name comes historically from the “Fly’s Eye” experiment, the conical substructure that makes up the eyes of invertebrates) receives the signals from a four by four “sub-cluster” of PMTs in the camera. Four by four of these sub-clusters make up the 16 by 16 PMTs in the camera. The OMB splits each of the 16 PMT signals into a trigger circuit and a charge-to-digital (QDC) circuit. The trigger circuit has a 375 ns low pass filter and discriminator with adjustable threshold. The QDC circuit has a 1.6 µs delay line and two integrator channels. When the trigger discriminator’s input crosses its threshold, the 21 Figure 2.7: The image of the TA MD electronics rack. 22 Figure 2.8: The diagram of the MD data acquisition crate. From the left, there are a CPU board, PPG board, trigger board, 16 OMBs, and garbage board. Reprinted with permission from J.D. Smith. two integrators start integrating the delayed signal. One integrator (Channel A) integrates for 1.2 µs and the second (Channel B) integrates for 5.6 µs. Both integrators will reset after 25 µs if no higher level trigger occurs. The discriminator also starts a third integrator with constant input that stops when hold-off (end of event) occurs. This integrator records the relative time of the PMT trigger and is called the time-to-digital converter (TDC). When three QDC signals in a OMB trigger within 25 µs and two are physically adjacent, a sub-cluster trigger is generated and all active integrators are prevented from reseting for an additional 25 µs. We use Channel B for MD analysis. Figure 2.9 shows a schematic outline of MD data telescope acquisition electronics. The trigger board receives the sub-cluster trigger signals from the 16 OMBs, and if two sub-clusters trigger within 25 µs of each other, a telescope level trigger is generated. When a telescope trigger occurs, the OMBs prevent all triggered integrators from resetting until all the integrators have been digitized. The telescope trigger starts a 50 µs delay to allow additional tube level triggers to occur from the cosmic ray air-shower event. After the 50 µs delay, a hold-off signal is generated that stops and saves all OMB integrators (QDC and TDC) and signals the crate CPU to start digitizing the event. The trigger board also captures coarse time information received from the central timing rack and sends the 23 Figure 2.9: The schematic outline of the MD electronics. The tube signal is split into a discriminator (filter) and a delay line in Channel A and B and Channel B integrates for 5.6 µs and checks if a sub-cluster level and telescope level trigger conditions are sequentially are met. Reprinted with permission from J.D. Smith. hold-off signal to the central timing rack to be precisely time-stamped there. The coarse time recorded from the trigger board can be match off-line with the precise time-stamps from central timing to give each event a precise GPS time-stamp. The TDC data read from the ommatidial boards are PMT trigger times relative to this precise GPS timestamp. Detailed triggering information will be discussed in Section 2.2.4. The CPU (Motorola 68030) receives the hold-off signal and commands the OMBs to start digitizing all the QDC and TDC integrators. The CPU then reads the hold-off time information from the trigger boards and the digitized QDC and TDC values from the OMBs to build an event packet to be sent to a host computer and recorded to disk. The PPG board is used to inject test signals into the pre-amplifiers in the camera. This is done as a diagnostic test each night before taking data to determine if any channels are malfunctioning. In the background, the CPU monitors the individual PMT trigger rates and adjusts the trigger discriminator thresholds up or down to maintain a constant 200Hz PMT trigger rate. This maintains telescope stability with changing night-sky conditions as stars and clouds cross the field of view. The CPU also periodically (each minute) records the PMT trigger rates and thresholds, power supply voltages, temperature, and other environmental conditions (from the “garbage” board) to be recorded in the data stream. 24 2.2.4 Trigger As described in the previous section, there are three levels of event trigger: tube level trigger, sub-cluster level trigger, and telescope level trigger. Figure 2.10 shows three triggering conditions. Initially, an individual tube is above the threshold, which means a tube is triggered, then we call it a “tube level trigger”. If at least three tubes are triggered in 4×4 sub-cluster of tubes within a 25 µs window and at least one pair must be adjacent indicating a hint of a track, a “sub-cluster level trigger” occurs. When at least two sub-clusters are triggered within an additional 25 µs and these sub-clusters are adjacent, a “telescope level trigger” occurs. The DAQ then adds an additional 25 µs hold-off time to collect more data by tubes that triggered. It also sends a signal to neighboring telescopes to save the information, which is called “nearest neighbor trigger”. This helps us to not lose any small signals from the shower in adjacent telescopes that may not be triggered. 2.2.5 Calibration As a photometric calibration, MD uses a standard candle called the Roving Xenon Flasher (RXF). The RXF, shown in Figure 2.11, consists of five parts: a Xenon flash lamp, a UV band-pass filter (300-400 nm), a narrow band filter (355 nm), a neutral density filter, and a teflon diffuser. The xenon flash lamp has a wide spectrum, its light passes first through a piece of band-pass filter 300-400 nm, which is the same as is in front of the Figure 2.10: A model of MD trigger patterns. Starting with single PMT triggers, if there are at least three PMTs in the sub-cluster, with one pair adjacent, then the sub-cluster level trigger is met. If at least two sub-cluster triggers with one pair adjacent, then the telescope level trigger occurs. Taken from Allen (2012). 25 Figure 2.11: Photo of the RXF module. The black box contains the Xenon flash lamp itself. The black cylinder at the front of the RXF contains the optical filters (band pass and narrow band) as well as the teflon diffuser, and a brass aperture. at the rear of the black box is a 1/2” brass rod that is used to mount the RXF in the center of post. The grey box below the black box contains the electronics for charging and firing the RXF flash lamp. Reprinted with permission from J. N. Matthews. camera. Next it passes through a narrow band filter for 355 nm light. The neutral density filter is to reduce effective brightness of flasher. Finally, it passes through a teflon diffuser to uniformly illuminate the PMTs. The RXF calibration is typically performed once per month at the beginning of the run, and while the telescope shelter doors are closed. The RXF emits 1µs long pulses of light at about 1.5 Hz. Approximately 500 RXF events are recorded in order to get a good measurement of the mean and width of the distribution. Before starting the calibration, the RXF module is acclimated to the temperature of the MD telescope, building for about 30 minutes to improve temperature stability. For the calibration, as Figure 2.12 shows, the RXF is mounted in the center of post of the telescope mirror so that it can uniformly illuminate all of the camera’s PMTs. The RXF calibration is performed from telescope 1 to telescope 14 and then comes beck to telescope 1 to check variation of RXF output over time. It usually takes about 5 minutes per telescope. The number of photons produced per xenon flasher event varies only slightly from pulse to pulse, about 0.3% and the temperature dependence is about -3%/10◦ C between -15 - 40◦ C. 26 Figure 2.12: Photo of RXF placed at the center of the mirror. At the bottom of the mirror, one can see the 12 V battery (grey and blue colored box) that is providing power to the RXF and the UVLED, which is normally mounted in the center post and provides a nightly relative calibration. Reprinted with permission from J. N. Matthews. 27 2.3 TALE The Telescope Array has added an extension to lower the experiment’s energy threshold, allowing us to study cosmic rays with energies down below 1016 eV. This extension is known as the Telescope Array Low Energy extension (TALE). High-elevation-angle telescopes and a graded infill array of more closely packed surface detectors were added to help make measurements of lower energy cosmic ray showers. The TALE extension telescopes view the sky above the MD telescopes and increase the elevation viewing angle up to 59◦ . A map showing the locations of the TALE SD detectors is shown in Figure 2.13. Detailed information will be discussed in the following sections. The TALE fluorescence detector was refurbished from the HiRes experiment. The telescopes all came from HiRes-2, which has FADC data acquisition electronics as indicated in Table 2.3. The mirrors for the telescopes came from HiRes-1/2 as indicated in Table 2.4. The TALE telescope shelter is located next to the MD shelter and is operated mainly by the University of Utah group with help from the Russian groups. 2.3.1 Mirrors and PMTs The TALE fluorescence station has 10 telescopes and each telescope consists of a spherical mirror that is composed of four mirror segments arranged in a cloverleaf shape with a total area of 5.2 m2 identical to the MD telescopes. The TALE telescopes view 30◦ -59◦ in elevation and have roughly the same azimuthal field of view as the MD. Five of 10 telescopes view 30◦ -45◦ (ring 3) in elevation, while the rest of the telescopes view 45◦ - 59◦ (ring 4) in elevation. Figure 2.14 shows a pair of TALE mirrors covering viewing angle of 30◦ -59◦ in elevation. As we mentioned in the previous section, the TALE is very similar to MD with respect to mirrors and UV band-pass filter. Unlike PMTs in MD, the TALE has slightly different PMTs with an extra focusing dynode to increase effective area and improve uniformity of the PMTs. These PMTs are an improved version of Philips/Photonis XP3062/FL model. Thus, it has about 10% larger effective area than the MD PMTs, which is 1197 mm2 . Using the MD and TALE detectors together, the site operates to observe air shower events from 3◦ -59◦ in elevation continuously. This enables us to observe the full development of lower energy showers that initiate closer to the observatory. With the larger field of view, we can also see longer tracks, which gives us additional information about shower geometry and profile. Figure 2.15 and 2.16 show MD and TALE telescope buildings and their field of view on the sky. This dissertation focuses on the spectral analysis using the combined MD and TALE detectors. The goal is to determine the flux of cosmic rays and to produce 28 Figure 2.13: Map of Telescope Array and TALE detectors. Locations of the main Telescope Array scintillator Surface Detectors (SD) are shown as red circles and the locations of the three Fluorescence Detector (FD) stations are indicated by the blue circles. The TALE is located at the MD site on the northern edge of the array. Locations of the scintillator Surface Detectors for TALE are shown as yellow circles. The purple circle centered between the three FD stations indicates the location of the Central Laser Facility (CLF). Reprinted with permission from R. Cady. 29 Table 2.3: TALE telescope serial numbers are listed with corresponding HiRes-2 telescope serial numbers, PMT manufacturer, and ring number. These are listed for historical tracking of the optical components. TALE telescope # 15 16 17 18 19 20 21 22 23 24 HiRes-2 telescope # 39 02 03 04 05 06 07 08 09 10 PMT Philips Philips Philips Philips Philips Philips Philips Philips Philips Philips ring # 4 3 4 3 3 4 3 4 3 4 Table 2.4: TALE mirror serial numbers are listed with corresponding HiRes-1 or HiRes-2 mirror serial numbers. These are listed for historical tracking of the optical components. TALE mirror # 15 16 17 18 19 20 21 22 23 24 HiRes mirror # HiRes-1 11 HiRes-1 15 HiRes-2 42 HiRes-2 41 HiRes-1 14 HiRes-1 10 HiRes-2 13 HiRes-2 38 HiRes-1 13 HiRes-1 08 30 Figure 2.14: A pair of TALE telescopes with viewing angle of 30◦ -59◦ in elevation. From the photo, it is apparent that these telescopes are tilted up significantly as compared with the MD telescopes. Reprinted with permission from J. N. Matthews. 31 Figure 2.15: Photo of the MD (left) and TALE (right) telescope buildings with opened shelter doors. The shorter main TA building contains 14 telescopes viewing 3-31◦ elevation. The taller building on the right has 10 TALE telescopes viewing 31-59◦ elevation above the field of view of the main TA telescopes. Reprinted with permission from J. N. Matthews. 32 60 Elevation [Degree] 50 40 30 20 10 0 240 260 280 300 320 340 Azimuthal angle [Degree] Figure 2.16: The field of view in the sky of the MD and TALE telescopes. Each point represents the center of the azimuthal and elevation field of view of one PMT. The field of view of a whole camera is box-like near horizon, but expands wedge-like as elevation increases due to the projection of the sky onto a plane. The bottom two rings belong to MD telescopes and top two rings belong to TALE telescopes. 33 a cross-calibrated monocular spectrum covering the 1017 decade in energy. For the MD and TALE analysis, we begin to consider the parallax effect between two sets of telescopes since the MD and TALE shelters are about 50 m apart. At distances of 1 km where the infill array of surface detectors starts, this is about a 3◦ effect. Distant events would not be expected to show a large parallax effect, but close-by events would be affected. This parallax is taken into account in the event reconstruction. To accurately measure the amount of the light, we should consider two factors; obscuration of the mirror and the atmospheric conditions including aerosols. Obscuration of the mirrors by the PMT cluster, stand, and cable trays reduces the total effective area of the mirror. All this obscuration is taken into account in the event reconstruction. The weather conditions affect the amount of light transmitted by the atmosphere. Hourly-based weather monitoring is carried out as MD observation, thus the same weather code is entered. For the TALE detector, checking overhead is important since TALE telescopes are looking higher up in the atmosphere. 2.3.2 Mirror reflectivity The mirror reflectivity is an important factor determining the fraction of photons from air showers reaching the PMTs. TALE mirrors are washed at the same time MD mirrors are washed and mirror reflectance per mirror is measured before and after washing. For TALE analysis, we use wavelength dependent mirror reflectivity curves obtained by data from May and June of 2015. Post-washed mirror reflectivity that is scaled down by about 3% is used because most of the time, mirrors are exposed to dust from the desert. In general, the reflectivity of HiRes-1 mirrors is 4-5% lower than that of HiRes-2 mirrors because of production differences. 2.3.3 Electronics Two telescope cameras share one electronics rack. Each has its own crate, which contains link modules, power control boards, 16 FADC boards, trigger/host boards, power supplies, and cooling fans. The photos of the front/back side of electronics crate are shown in Figures 2.17 and 2.18. The link module handles all communication and clock distributions to the FADC boards. The power control boards enable power to the FADC electronics rack, communicate to the link module, and monitor voltages, temperature, humidity, and shelter doors. The 16 FADC boards sum each row and column of PMT signals and pass discriminated row and column triggers to the trigger/host board. The FADC boards then digitize each PMT signal channel 34 Figure 2.17: Photo of front side of TALE electronics rack. 35 Figure 2.18: Photo of back side of TALE electronics rack. 36 (16 channels/board) and a high- and low-gain version of the row and column sums (one high- and one low-gain sum of one row and one column on each board) with fast analog to digital convertors with a sample period of 100 ns. The output of the high-gain sum channels is used to determine the trigger condition of the camera. The output of the low-gain sums is used to better understand saturated channels within the camera. If the signal sum is in excess of 12 ADC counts above the pedestal, the signal is discriminated and a trigger bit is set to be high. If it does not meet threshold, the bit is set to be low. It then sends the bit to the trigger/host board. The trigger time window is calculated as 2×tcoincidence + 10 µs and can be up to 200 µs. The trigger/host board monitors the horizontal and vertical signal sums to scan for patterns and saves information if the trigger condition is met. 2.3.4 Trigger The primary trigger occurs when the trigger board finds a 3-fold coincidence of PMT rows or columns on the horizontal or vertical signal sums. It analyzes discriminated signals of PMT row and column signal sums. There are 16 signal sums from PMT rows and 16 signal sums from PMT columns. Each signal sum goes through a logical operator “AND” with the result from a logical operator “OR” of the next two signal sum channels. For example, signal sum channel 0 goes through logical operator “AND” with the result from logical operator “OR” of signal sum channels 1 and 2. Thus, at the first step, 16 signal sum channels pass through this logic and 15 signal sum channels of 2-fold coincidences remain. By following the same process, 14 signal sum channels of 3-fold coincidences appear. The diagram of this trigger logic is shown in Figure 2.19 and an example of that is shown in Figure 2.20. More detailed information can be found in Zundel (2016). 2.3.5 Calibration As a photometric calibration, TALE uses the UVLED calibration. This UVLED module is placed at the center of each mirror to illuminate the cluster. The UVLED calibration is taken at the beginning and end of MD and TALE runs and during runs. Unlike RXF calibration, which is done once a month, UVLED calibration is done twice per night and each mirror has its own UVLED module. As Figure 2.21 shows, the UVLED module is mounted in the center of post of the telescope mirror so that it can uniformly illuminate all of the camera’s PMTs. It emits about 500 shots of 355 nm UV light. Their UVLED modules are designed to be temperature stabilized within 1/10◦ C of 45◦ C. It continuously flashes UVLED light on the PMT cluster during observation time to monitor the gain of detectors through the night for the MD telescopes, but not for the TALE telescopes. Compared 37 Figure 2.19: The diagram of TALE triggering. Starting with 16 signal sums, each signal sum goes through logical operators, thus 15 signal sum channels of 2-fold coincidences remain. After the same process, 14 signal sum channels of 3-fold coincidences appear. Taken from Zundel (2016). Figure 2.20: An example of TALE triggering. The solid lines are discriminated signals while dashed lines are not discriminated signals. Taken from Zundel (2016). 38 Figure 2.21: Photo of UVLED module placed at the center of a mirror. Reprinted with permission from J. N. Matthews. to RXF calibration, UVLED calibration provides nightly detector gains. Additionally, the PMT gain is balanced to 1 photoelectron/ADC at the beginning of runs (1st night of month). CHAPTER 3 EXTENSIVE AIR SHOWER PHYSICS Ultra-High Energy Cosmic Rays (UHECRs) induce a cascade of secondary particles when they collide with air molecules high in the atmosphere. The secondary particles continue to travel down through the atmosphere at nearly the speed of light. These secondary particles produce more particles as they interact with more air molecules, causing the extent of the shower to cover many kilometers across the ground. This phenomenon is referred to as an air shower, and, if such air showers are large enough, they are known as Extensive Air Showers (EASs). In EASs, about 1010 secondary particles are created. Within the showers, photons, electrons, and positrons, which collectively are referred to as the electromagnetic components, carry about 85% of the total energy; muons carry about 10% of the total energy; pions carry about 4%; and neutrinos carry the remainder. As we discuss below, EASs largely develop through electromagnetic and hadronic cascades. We begin with simplified models for the electromagnetic and hadronic cascades and describe important features of EAS development at the discrete development level, but the actual development process is continuous. This is followed by a discussion of the limitation of these simplified models comparing to the reality. 3.1 Electromagnetic Showers In 1954, Heitler (Heitler 1954) presented a very simple model for the development of electromagnetic (EM) cascades. The pure EM showers are composed of only electrons, positrons, and photons, as depicted in Figure 3.1. Electrons and positrons may radiate one photon by bremsstrahlung after travelling one interaction step length, X, defined as X = X0em ln 2, (3.1) where X0em is the radiation length in the medium (e.g. in air, X0em ≈ 37 g/cm2 (Patrignani et al. 2016)) and is proportional to the interaction step length. It is defined this way so that, on average, the particle loses half of its energy, E, by radiation, as described by 1 em E = E e−X/X0 . (3.2) 2 40 Figure 3.1: Schematic view of an electromagnetic cascade induced by a cosmic ray photon. At each interaction step length, a photon creates an electron-positron pair, known as the “pair production”, and an electron or positron emits a photon by interacting with another charged particle, known as the “bremsstrahlung”. 41 Photons, meanwhile, produce an electron and positron through pair production. In general, after one interaction step length, the number of particles is doubled with the energy equally divided between the two outgoing particles. After n interaction step lengths in the atmosphere, the distance travelled is Xn = nX0em ln 2, (3.3) and the size of the shower is Xn em Nn = 2n = e X0 . (3.4) The shower development continues until the individual energy of the particles drops below a critical energy, εem (e.g. in air, εem = 85 MeV (Matthews 2005)), where the rate of energy loss of electrons via bremsstrahlung is equal to the rate of energy loss via ionization. Above εem , bremsstrahlung and pair production are more likely to occur than ionization. Below εem , the probability of energy loss by bremsstrahlung or pair production is less than that of ionization, thus ionization is the dominant mode of energy loss. When the average secondary particles are at the critical energy, the shower reaches its maximum size at a depth, Xmax , in the atmosphere. The depth of the electromagnetic shower maximum happens after n c interaction step lengths, and is given by em Xmax = nc X0em ln 2. (3.5) At this depth, all particles have energy, εem , so that total energy, E0 , is em E0 = εem Nmax . (3.6) Hence, the number of particles at this depth is E0 em (3.7) = 2nc , = Nmax εem and it is proportional to the incoming primary cosmic ray energy. Thus, the amount of light emitted is proportional to the number of charged particles and we are able to measure this light via our fluorescence telescopes. Such a measurement is vital as it gives crucial information about how the air shower develops. We can, therefore, determine n c using Equation 3.7, ln E0 εem . (3.8) ln 2 We can then restate the maximum size at depth in atmosphere, Xmax , with the total energy, nc = E0 . It can be expressed as em Xmax = nc X0em ln 2 = X0em ln E 0 . (3.9) εem We can see that it evolves logarithmically with the primary cosmic ray energy. This has a critical role in identifying the composition of incoming primary cosmic ray particles. 42 The elongation rate, Λ, is the rate of increase of the mean Xmax per decade in energy, defined as dXmax . (3.10) d log10 E0 This variable is also sensitive to the mass composition of the primary cosmic ray particle. Λ≡ To obtain the elongation rate of the EM shower in air, we may apply Equation 3.9, finding that Λem = 2.3X0em . We can see that the elongation rate is dependent upon the radiation length of the medium. For a pure EM shower in air, the elongation rate is about 85 g/cm2 . However, Heitler’s model assumed that all of the processes’ cross-sections are treated as independent of the energy of particles, while in reality, the cross-section depends on the energy of particles. In addition, collision energy losses are neglected, especially for particles with energy greater than a critical energy, εem . From the energy distribution of particles, it is assumed to be equally divided over the secondary particles. However, in reality, the distribution of energy is highly inhomogeneous. In addition, in this simplified approach, the ratio of electrons to photons is overestimated by about a factor of 2 to 3 and it predicts a ratio 2 while actual experimental measurements show a ratio of the order of 1/6. This means that in reality, multiple photons can be emitted by bremsstrahlung and low energy electrons can be absorbed by air molecules and lose energy much faster than photons do. Therefore, a correction factor, g (i.e. a constant value, g = 10), is adopted to more accurately estimate the number of charged particles, which is given as N em Ne± = max . (3.11) g This correction factor is determined by comparing the model to actual experimental measurements and it should be considered an order of magnitude estimate (Matthews 2005). With the correction factor, the elongation rate of the EM shower in air becomes smaller than about 85 g/cm2 that we estimated using Heitler’s simplified model. This means that the actual EM shower develops more rapidly than we expected using the simplified model. 3.2 Hadronic Showers In a Heitler-like model, the development of hadronic showers is similar to the development of electromagnetic showers, which is based on a very approximate model. Consider a single cosmic ray proton striking an air molecule with a total energy, E0 , as depicted in Figure 3.2. The interaction mean free path of the incoming primary cosmic ray is described as 43 Figure 3.2: Schematic view of a hadronic cascade induced by a cosmic ray proton. At each interaction step length, charged pions are created and keep interacting until their energy is equal to the critical energy and neutral pions are produced and quickly decayed to photons, producing electromagnetic sub-showers. 44 A , (3.12) NA ρ σ where λ is the mean free path of the cosmic ray coming through the atmosphere, A is the λ= mass number of the nucleus, NA is Avogadro’s number, ρ is the density of the atmosphere, and σ is the cross-sectional area for the collision. When the primary cosmic ray proton traverses an interaction step length, it produces many secondary particles, mostly pions. The interaction step length, X, can be defined as X = X0had ln 2, (3.13) where X0had is the nuclear interaction length in the medium, given by 35 g/cm2 A1/3 (one interaction length in air or X0had ≈ 90 g/cm2 (Patrignani et al. 2016) and one interaction length for pions, X0π ≈ 120 g/cm2 (Gaisser et al. 2016)). The secondary particles are all created with an equal amount of energy, regardless of the number of interaction step lengths traversed. At each interaction, Nch charged pions and 1 2 Nch neutral pions are generated. The charged pions carry an energy of 32 E0 and the neutral pions carry the remainder of the total energy, 13 E0 . The neutral pions decay to two photons immediately, creating electromagnetic subshowers; their decay is described as π 0 → 2γ. However, the charged pions continue interacting with air molecules until they reach a critical energy, επ , where the charged pions often usually decay to muons and neutrinos. The charged pion decay is described as π + → µ+ + ν µ and π − → µ− + ν̄µ . In other words, when the energy of the charged pions is below the critical energy, επ , it is difficult for them to travel another interaction step length. This critical energy, επ , and the number of interactions depend on the initial energy of the incoming primary cosmic ray particle. Therefore, the hadronic showers are a combination of electromagnetic cascades induced by the neutral pions and hadronic cascades induced by the charged pions. The electromagnetic portion continues to grow via the neutral pions feeding it and the reminder is due to the charged component. The total energy of the hadronic shower can be calculated by considering the number of electrons from electromagnetic sub-cascades and the number of muons from the hadronic cascades because the number of the charged pions is the same as the number of the muons: em E0 = Nmax εem + Nµ επ . (3.14) 45 Measuring the numbers of charged particles (i.e. electrons, positrons, and pions) is therefore of prime importance for estimating the initial energy of the primary particle. Using Equation 3.11, the total energy of the hadronic shower can be described as em εem + Nµ επ E0 = Nmax = gNe± εem + Nµ επ (3.15) επ = gεem Ne± + Nµ , gεem and it can be determined if both the number of electrons and muons are measured. As we see in Equation 3.15, the relation between the total energy and the number of electrons or the number of muons is linear. In air, this relation is E0 ≈ 0.85 (Ne± + 24Nµ ) GeV. (3.16) Here, understanding the relation between the number of muons and primary cosmic ray energy is important in order to predict the missing energy that is carried by muons. The number of muons can be estimated using the condition for when pions reach the critical energy, επ . After n c interaction step lengths, pion energy is at the critical energy, described as En c = E0 3 ( 2 Nch )nc = επ . (3.17) To determine n c , we may rearrange Equation 3.17 as n c E0 3 = Nch επ 2 3 E0 = nc ln Nch ln (3.18) επ 2 0 ln E επ . nc = ln 32 Nch We can estimate the number of muons by considering the number of charged pions. At each interaction, Nch charged pions are created and when they fall below the critical energy, επ , the charged pion decays into one muon and one neutrino as π + → µ+ + νµ or π − → µ− + ν̄µ . Thus, the number of muons can be expressed as Nµ = Nπ± = (Nch )nc . (3.19) Using the expression for n c from Equation 3.18, we can rewrite the number of muons as β E 3 E0 ln ε 0 / ln 2 Nch π Nµ = (Nch ) . (3.20) ⇒ ln Nµ = ln επ In the above, β is defined as ln Nch/ln 32 Nch and this tells us that the number of muons is proportional to a power law of the primary energy with a spectral index β: Nµ ∝ (E0 )β . We can use this to estimate the fraction of energy in the hadronic cascade. The fraction of energy that goes into a hadronic cascade is related to the missing energy that muons and neutrinos take from the shower and is not observable. The energy of a hadronic cascade, 46 Eh is Eh = Nµ επ = (Nch )nc επ . (3.21) The total energy can then be written as n c 3 E0 = επ , (3.22) Nch 2 and thus we find that the energy fraction of a hadronic cascade is 2 n c Eh (Nch )nc επ . (3.23) = 3 = E0 3 ( 2 Nch )nc επ This tells us that higher primary shower energies, which allow for greater numbers of interaction step lengths to be traversed, have smaller hadronic components. To determine the depth at which the maximum size of the hadronic cascade occurs, p Xmax , we consider the depth at which the EM sub-showers of a hadronic cascade reach their maximum size. We know that after the first interaction length, 1 2 Nch neutral pions are created and they decay to Nch photons as π 0 → 2γ. Each photon initiates an electromagnetic shower of energy 1 3 E0 E0 . (3.24) Nch 3Nch From this, we can calculate the maximum size at depth of the electromagnetic component Eem = = p of the shower included by cosmic ray protons in the atmosphere, Xmax , which is E 0 p p Xmax = X0 + X0em ln 3Nch εem E0 (3.25) p em − X0em ln(3Nch ) = X0 + X0 ln εem em = X0p + Xmax − X0em ln(3Nch ), p where we have used Equation 3.9 in the second step. We see that Xmax is composed of the first interaction length of the primary proton, the depth at which the maximum EM sub-showers occur, and the number of the charged pions created in the first interaction in the hadronic cascade. This provides us with a hint on the composition of primary cosmic ray particles. By employing Equation 3.10, we may determine the elongation rate of a hadronic shower, Λp , to be Λp = p dXmax d log10 E0 d X0em ln(3Nch ) dX0p em = +Λ − . d log10 E0 d log10 E0 Using the inelastic proton-air cross-section and proton-proton multiplicity data, we find that Λp in air is ∼ 58 g/cm2 per decade in energy (Matthews 2005). However, the proton-air cross-section that affects the full air shower development depends on the energy of the incoming primary cosmic ray particle. It rises with energy and 47 results in different shower maximums, Xmax , that constrain the composition of the incoming cosmic ray particle (Abbasi et al. 2015). At the present, it is not feasible for current particle accelerators to reach the particle energy of the highest energy cosmic rays. Therefore, the proton-air cross-section must be calculated using the extrapolation of the hadronic models to high energies. Moreover, in this simplified approach, the energy of secondary particles is taken equally; however, in reality, there are variations in the number of particles and energy of secondary particles so that energies of particles are not equally distributed. In addition to cosmic ray protons, we can suppose that heavier nuclei are incident on the atmosphere and induce air showers. In a simple model, heavier nuclei can be considered as the superposition of many nucleons. For a nucleus with mass number A and total energy E0 , each nucleon has an energy E0 /A and A individual showers all starting at the same point are superposed. Using Equations 3.20 and 3.25, we can estimate the number of muons, A , as we see NµA , and the depth where it reaches the maximum size in the atmosphere, Xmax below: NµA = A E0 /A επ β = and E0 επ β A1−β = Nµp A1−β (3.26) E0 /A = + 3Nch εem E0 (3.27) = X0p + X0em ln − X0em ln A 3Nch εem p = Xmax − X0em ln A. The iron nuclei showers contain about 40 % more muons than proton showers at the same A Xmax X0p X0em ln total energy and the depth of maximum size becomes about 80-100 g/cm2 shallower, which means they reach Xmax earlier and higher in the atmosphere. Also the fluctuations in Xmax should be smaller for heavy incoming cosmic rays due to many sub-showers superposed. In other words, protons are expected to have the longer interaction mean free path, resulting in wider distributions, while heavier nuclei, e.g. iron nuclei, have a shorter interaction mean free path, resulting in a more reliable and narrow distribution. Thus, the depth when the showers reach the maximum size is a critical estimator of the composition of incoming cosmic rays. CHAPTER 4 MONTE CARLO SIMULATION For Fluorescence Detector (FD) measurements, the aperture depends significantly on energy. Since higher energy events are brighter, the telescopes can see these showers from greater distances. On the other hand, lower energy events are only visible nearby. The aperture and acceptance of the detector must be calculated via the Monte Carlo (MC) simulations. The aperture, AΩ, can be written as NReconstructed (E) × A0 Ω0 , (4.1) AΩ = NThrown (E) where NReconstructed (E) is the number of reconstructed MC events that pass quality cuts; NThrown (E) is the number of thrown MC events; A0 indicates the effective area that is determined by thrown maximum impact parameter, Rp ; and Ω0 indicates the solid angle that is determined by thrown maximum zenith angle, θ. Here, NReconstructed (E) NThrown (E) represents the acceptance. To perform the aperture calculations, we need a complete simulation of the detector to create MC samples that model the data. We input previous measurements such as the energy spectra and compositions. We use a library of showers generated by the COsmic Ray Simulations for KAscade (CORSIKA) program for the theoretical modeling of the cascade. Events are generated and thrown according to the expected distributions and they undergo atmospheric scattering, and light transmission. Once they arrive at the detector, the detector response, trigger, and readout electronics are simulated by the detector Monte Carlo that is specific to the Telescope Array. The next step is to write the MC events into the same format as the data and analyze the MC event with the exact same programs as are used for data analysis. Finally, we compare analyzed data with the analyzed MC. If their distributions agree, this demonstrates that we understand the data and our detector well enough to make measurements. We will discuss the general function and the features of particle physics applied in the CORSIKA in Section 4.1. We will also compare the thrown distributions with theoretical expectations and examine the resolution of the simulations in Section 4.2. 49 4.1 Monte Carlo Program The CORSIKA program is a detailed Monte Carlo program for simulating the development of extensive air showers using high energy strong and electromagnetic particle interactions initiated by cosmic ray primary particles. The CORSIKA program tracks secondary particles in a “stack”, listing particles’ position and momentum vectors, particle type, and other parameters. In the CORSIKA package, it is possible to select and use one of several different hadronic interaction generators (e.g. QGSJet, DPMJET, SIBYLL, VENUS, EPOS, etc.) and different electromagnetic interaction generators (e.g. EGS4). Figure 4.1 shows an example of a Monte Carlo shower generated using CORSIKA program with QGSJET II3 model. All hadronic and electromagnetic models are based upon particle interactions and cross-sections measured in the laboratory. However, they suffer from the fact that these measurements must be extrapolated to energies that are unfeasible in the laboratory. The knowledge of high energy particle interaction is incomplete and is explained using theoretical models with extrapolation of hadronic interactions to higher energy, which is not yet explored by current experimental data. Accelerators are only capable of reaching energies of about 1013 eV. Even then, collider experiments see all transverse momenta within their acceptance, which is at large angle to the beam lines. The parts of a shower seen by a fluorescence telescope go down the beam pipe at a collider experiment. Thus, for the first few interactions in the EAS, particle interactions in the air shower are not well understood. However, this part of the air shower is difficult to detect using fluorescence detectors due to the limited number of particles in the shower. Sufficient light for FD measurements begins as the shower reaches the depth of maximum development, Xmax . At this stage, the particles’ energies are similar to Fermi National Accelerator Laboratory (FNAL)’s energies, making FD measurements of Xmax crucial as the particle interactions are well understood here. For example, in a shower initiated by a 1017.5 eV primary proton, individual pion energies after about five interactions in the Heitler-Matthews model reach the energy of particles at the FNAL. For FNAL target energies, ECM = 27 GeV, which corresponds to a particle with an energy of about 400 GeV in the shower. This means that after several interactions, experimental results from FNAL provide a fairly good approximation for particle interactions. In other words, when the average energy of shower particles meets the energy of FNAL, hadronic models are able to give more accurate particle interactions and cross-sections from FNAL fixed target experiments. A proton primary with an energy of 1017.5 eV will reach a shower maximum of about 50 Figure 4.1: An example of a Monte Carlo shower. This simulated event is a 1015 eV iron shower generated using CORSIKA program with QGSJET II-3 for the hadronic interaction model. Taken from Zundel (2016). 109 particles and its core will have a diameter of about 50-100 m. However, the shower will be viewed from several kilometers away by the FD telescope. The only part of the shower to have enough particles to be bright enough to be observed will be the shower core. The telescopes are insensitive to the outer parts of the shower with large transverse momentum. For the case of the tail of the shower, secondary particles can be distributed several kilometers from the axis of EAS. The tail of the shower has low density of particles with large transverse momenta, and is poorly measured by fluorescence detectors. 4.1.1 CORSIKA For all Monte Carlo simulation work done for this analysis, the QGSJet (Quark Gluon String model with JETs) hadronic and the EGS4 (Electron Gamma Shower system version 4) electromagnetic interaction models are used. The QGSJetII-3 model is based on GribovRegge theory, which estimates what a Quantum ChromoDynamics (QCD) solution might look like and attempts to explain a slowly rising cross-section of hadronic collisions at high 51 energies. The EGS4 model was developed at Stanford Linear Accelerator Center (SLAC) to simulate electromagnetic cascades in various geometries and energies up to a few thousand GeV and it is a well understood and trusted simulation package. Simulating an air shower initiated by cosmic rays with energy > 1016.5 eV and tracking all secondary particles in CORSIKA is impractical with respect to computation time and size of output file. Alternatively, using the “thin sampling” option, thinning level (Eth ) is set to a certain energy threshold and it reduces the number of particles being tracked in the simulation. If secondary particles fall below the thinning energy, which is defined as Eth ×E0 , where E0 is the primary energy, only one of the secondary particles is tracked and saved with weight factor, while all other such secondary particles are discarded. In our Monte Carlo simulation, CORSIKA showers are generated with a 10−5 eV thinning level. While using the “thin sampling” option helps to reduce computation to a reasonable time, we lose information in the process. Since the thinning algorithm tracks the information of an average particle, it is difficult to express the statistical fluctuations on particle density distributions as a function of distance from the shower core. The FD analysis, however, only looks at the core and these fluctuations are not important to my analysis. 4.1.2 Input For Monte Carlo, we use the mc2k12 program to generate Monte Carlo events using input cards. The energy of shower, shower size, Gaisser-Hillas parameters, unit vectors of track and impact parameter, triggered mirrors and tubes, photoelectrons received by tubes per event, and time information are recorded in the MC04 DST bank. The parameters in the input card are shown below: • setNr shows the information of date [YYYYMMDD] and data file number [##]; [YYYYMMDD##]. • use DB selects to use the database of calibration; [YES or NO]. If YES, the RXF calibration is used for MD and UVLED calibration is used for TALE. If NO, it uses the nominal constant values for calibration. • iseed is the random seed number and is usually a negative number. • detector shows detector configuration setting. It can be set to MD detector only, TALE detector only, or MD and TALE detectors together; ta md.conf, ta tale.conf, or ta md tale.conf, respectively. For this analysis, the configuration is to use MD and TALE data. 52 • shift origin means the coordinates of events and it is set normally to use the Central Laser Facility frame; [YES or NO]. • nevt means the number of events requested to generate. • event type is input as SHOWER, which means event type shower. • ipot spectr means one energy power law spectrum. • gamma is input as 3, which means the generated power law’s exponent and spectral index. • minEnrgy indicates the minimum shower energies in eV (e.g. 1.0e+16). • maxEnrgy indicates the maximum shower energies in eV (e.g. 1.0e+19). • primary indicates the hadronic model and chemical composition of the cosmic ray primary: qgsjetii-03, proton. The hadronic interaction model and primary type are needed. • rpmin : minimum impact parameter in meters (e.g. 10.0). • rpmax : maximum impact parameter in meters (e.g. 50,000.0). • thesh1 : minimum zenith angle in radians. • thesh2 : maximum zenith angle in radians. • phish1 : minimum azimuthal angle in radians. • phish2 : maximum azimuthal angle in radians. • phirp1 : limits on the azimuthal angle of the impact parameter in radians. • phirp2 : limits on the azimuthal angle of the impact parameter in radians. • dxlim : limit of depth tracing. • hceil : the atmosphere ceiling in meters. 53 4.1.3 Weighting In order to make the Monte Carlo events more realistic, we applied the previously measured energy spectra and composition. The Monte Carlo is thrown as the power law with exponent, -3; however, in reality, the energy spectrum does not follow a simple power law and we have been observing some spectral features such as the 2nd knee, ankle, and GZK cutoff above 1016 eV. As the spectral reference, we used the Telescope Array energy spectrum, which is shown in Figure 4.2. It was presented in the International Cosmic Ray Conference (ICRC) 2015. During generation of the energy spectrum, the events are weighted according to this spectrum, which is normalized at 1017.3 eV. The HiRes/MIA experiment observed a change in composition in the energy region from 1017 -1018 eV, which is shown in Figure 4.3a. Using mean the Xmax results, we obtained proton fraction as a function of energy, which is shown in Figure 4.3b. This was used to weight the number of events between primary proton and iron. Figure 4.2: The energy spectrum measured by TA in 2015. Black dots show the TA SD measurement. Blue squares show the TA BR-LR monocular measurement. Red dots indicate the TALE spectrum, which plays the role of a bridge between TALE Cherenkov and SD spectrums. Red open squares show TALE spectrum using Cherenkov events. Note that this plot uses flux times E 3 to take out the overall power law and show the fine structure of the spectrum. Taken from Ivanov (2016). 54 (a) (b) Figure 4.3: The mean Xmax from HiRes/MIA and HiRes measurements shown with inferred chemical composition. (a) The mean Xmax results from HiRes/MIA and HiRes experiments. The black dots show mean Xmax observed. Events between 1017 and 1018 eV were measured by the HiRes prototype in conjunction with the MIA muon array. The events above 1018 eV we remeasured by the HiRes stereo experiment. The red and blue lines represent Monte Carlo simulations of pure proton and iron MC events, respectively. This plot is interpreted to indicate that the composition of cosmic rays changes from heavy to light between 1017 and 1018 eV. Above 1018 eV, the composition looks light/protonic and unchanging. (b) A fit to the proton fraction as a function of energy from the measurement above. 55 4.1.4 Light production There are two kinds of light produced in the air showers. The first is fluorescence light that is emitted due to the interaction between charged secondary particles and nitrogen molecules. The fluorescence light is isotropically emitted between 300-420 nm with three strong peaks in the spectrum at 337, 357, and 391 nm. The second is Cherenkov light that is emitted when charged particles pass through the medium (i.e. air) at the velocity greater than the speed of light in the air. Cherenkov light is emitted in a narrow cone around the particle’s direction of travel. The amount of these light depends on the number of charged particles in the air shower. For the case of the fluorescence light, the number of photons produced per unit length per unit solid angle is given by d2 Nγ Y Ne = , (4.2) dldΩ 4π where Nγ is the number of photons, Y is the fluorescent yield in photons/electron/m, and Ne is the number of electrons in the air shower (Baltrusaitis et al. 1985). The fluorescent yield used here is based upon the measurement in Kakimoto et al. (1996) and it is defined as the number of photons produced by a charged particle per meter of travel in air. The measurement is carried out that electrons with energies 1.4 MeV, 300 MeV, 650 MeV, and 1000 MeV passed through air at various temperatures and pressures to obtain the formula of fluorescent yield which is, dE dx A1 A2 √ + √ , ρ Y = dE (4.3) 1 + ρB1 t 1 + ρB2 t dx 1.4 MeV dE where dE dx is the energy loss of electron, dx 1.4 MeV is the energy loss of electron having an energy of 1.4 MeV, ρ is the air density in kg/m3 , t is the temperature in Kelvin, and A1 , A2 , B1 , and B2 are constants. The fluorescent yield at 337 and 357 nm have the same pressure dependency. On the other hand, the fluorescent yield at 391 nm does not show pressure dependency compared to the other two major wavelengths. Thus, this fluorescent yield has two pressure dependencies containing the density and the square root of the temperature. By using this measurement, we are able to understand the fluorescent yield of air, which is proportional to the electron energy loss. This is important to estimate the energy of the cosmic ray primary particle. Thus, it follows that the systematic uncertainty in energy estimation from integration of the shower profile can be improved from 30% to 10%. In other words, it helps us to calculate an accurate energy estimation for the cosmic ray fluorescence detection technique. The fluorescent yield has a small dependence on altitude and that dependence is taken into account in the Monte Carlo and reconstruction programs. For Cherenkov light, the number of photons produced per unit length is given by 56 Z Z ∞ dNγ m2 c 4 2πα dEf (E) 2δ − dν , (4.4) = dl c E2 Et where Nγ is the number of photons, α is the fine structure constant, c is the speed of light, ν is the frequency of the radiation, and Et is the Cherenkov photo-production threshold √ energy in the air. The threshold energy is given by Et = mc2 / 2δ, δ is n − 1, with n being the index of refraction of air, m is the electron mass, and f (E) is the electron distribution in the air shower, which is defined below, Z F (E) = ∞ f (Ẽ)dẼ. (4.5) E This F (E) is the fraction of electrons in the shower with energies greater than E. By using 1.0 for the shower age, we calculate F (E) where the shower reaches the maximum size. F (E) is evaluated below, 34.8 (40.4 + E)(1 + 10−4 E)2 where E is the electron energy in MeV and it depends on the shower age slightly. F (E) = (4.6) Up to this point, we evaluated how many direct Cherenkov photons per meter are produced along the shower track. The angular distribution of scattered Cherenkov light is defined as dNγ e−θ/θ0 d2 Nγ = , (4.7) dldΩ dl 2π sin θ where θ is the angle of the Cherenkov beam cone and θ0 depends on the Cherenkov threshold energy and is defined as θ0 ≈ 0.83Et−0.67 . 4.1.5 Light propagation The generated light signals in an EAS propagate through the atmosphere and interact with the atmosphere as they propagate. These interactions result in scattering and absorption of the light. The first interaction we consider is Rayleigh scattering, which is a scattering from air molecules. It depends on the local air density (ρ) and wavelength of the light (λ). The second interaction we consider is Mie scattering, which is a scattering mechanism from dust, or aerosols. The third interaction we consider is ozone absorption, which depends upon the ozone concentration as a function of altitude. For the case of Rayleigh scattering, the amount of scattered light from a beam of Nγ photons can be defined as dNγ Nγ 400 4 = −ρ , (4.8) dl xR λ where xR is the mean free path. The mean free path is normalized to 2,974 g/cm2 at the wavelength 400 nm. The angular distribution of scattered light is expressed as d2 Nγ dNγ 3 = (1 + cos2 θ), dldΩ dl 16π where θ is the scattering angle. (4.9) 57 In the case of Mie scattering, the amount of light scattered from a beam of Nγ photons can be defined as where LM Nγ −h/HM dNγ , (4.10) e =− dl LM is 25 km and is the typical Mie scattering mean free path at the wavelength 360 nm and HM is 1.0 km, which is the scale height. The angular distribution of light scattered from the Cherenkov beam cone for angles in the range, 5-60◦ , is expressed as dNγ d2 Nγ (4.11) ≃ 0.80e−θ/θM , dldΩ dl where θM ≃ 26.7◦ . For the case of ozone absorption, the amount of light absorbed is defined as, dNγ = Nγ ρO3 (h)AO3 (λ), (4.12) dl where ρO3 (h) is the density of ozone at a given height, h, and AO3 (λ) is the coefficient at a given wavelength. 4.1.6 Light attenuation To calculate the amount of light arriving at the detector, we need to correct for the attenuation losses between the light source and the detector. There are three kinds of transmission factors defined. The Rayleigh transmission factor, TR , is given by |x1 −x2 | 400 nm 4 − ( λ ) xR , TR = e (4.13) where x1 and x2 are the locations of the light source and the detector in slant depth. The Mie transmission factor or Aerosol transmission factor, TA , is given by HM sec θ (e−h1 /HM −e−h2 /HM ) T A = e LM , (4.14) where h1 and h2 are the heights of the light source and the detector in g/cm2 . The ozone transmission factor, TO3 , is given by TO3 = e−∆xO3 AO3 , (4.15) where ∆xO3 is the integrated ozone density between the light source and the detector over the slant depth. Thus, a light beam propagating from the light source to the detector has an overall transmission factor, T , which is given by, T = T R T A T O3 . 4.2 (4.16) Monte Carlo Event We used the mc2k12 program and input cards to generate Monte Carlo events. Before we generate a big set of Monte Carlo events, we examined what the thrown distributions 58 look like and compare them to the predicted distributions of energy, distance, zenith angle, and azimuthal angle. These variables are uniformly thrown according to the input cards. We then compared the actual thrown distribution of the Monte Carlo set in Section 4.2.1 to the theoretical expectations. Next, we generated a large set of Monte Carlo events and processed them through the same reconstruction programs as the data to determine how well Monte Carlo events are reconstructed in Section 4.2.2. 4.2.1 Thrown distribution The thrown distributions are continuous, but when they are plotted in histograms, they can be expected to be dependent upon the total number of events simulated, the minimum and maximum limits of parameters, and bin size of histograms. The expected thrown distributions are as follows: Since the energy is thrown using the power law, E −3.0 , the energy distribution is N −2 ∆E − 10−∆E ), (4.17) −2 Ei (10 − Emax ) where N is the total number of thrown Monte Carlo events, Emin is the minimum energy f (E) = −2 (Emin simulated, Emax is the maximum energy simulated, Ei is the energy in a given histogram bin, and ∆E is the histogram bin size. Figure 4.4a shows the actual thrown energy distribution (blue) and the theoretical expectation (red). Since the impact parameter distance is thrown uniformly and randomly inside the circle of Rp distance distribution, is shown in Figure 4.4b. It follows 2N ∆Rp Rp , f (Rp ) = 2 (Rp max − Rp2 min ) where N is the total number of thrown Monte Carlo events, Rp simulated, Rp max max , the impact parameter (4.18) min is the minimum distance is the maximum distance simulated, and Rp is the impact parameter distance in a given histogram bin. Since the zenith angle is thrown using sin θ distribution, the zenith angle distribution, shown in Figure 4.5a, is as follows 2N ∆θ f (θ) = sin θi sin , (4.19) (1 − cos θmax ) 2 where N is the total number of thrown Monte Carlo events, θmax is the maximum zenith angle simulated, θi is the zenith angle in a given histogram bin, and ∆θ is the histogram bin size. Since the azimuthal angle is thrown using the flat distribution, the azimuthal angle distribution is N , (4.20) i where N is the total number of thrown Monte Carlo events and i is the total number of f (φ) = histogram bins. Figure 4.5b shows the actual thrown distribution (blue) and the theoretical expectation (red). 59 30000 # of events 25000 20000 15000 10000 5000 0 17.5 17.6 17.7 17.8 17.9 18 18.1 18.2 18.3 18.4 18.5 7000 8000 9000 10000 log (E/eV) 10 (a) 8000 7000 # of events 6000 5000 4000 3000 2000 1000 0 0 1000 2000 3000 4000 5000 6000 Impact parameter, Rp [m] (b) Figure 4.4: Thrown energy (E) distribution and thrown impact parameter (Rp ) distribution. (a) Thrown energy (E) distribution. Blue histogram shows the actual thrown distribution and Red line shows the theoretical expectation. (b) Thrown impact parameter (Rp ) distribution. Blue histogram shows the actual thrown distribution and Red line shows the theoretical expectation. 60 3500 3000 # of events 2500 2000 1500 1000 500 0 0 5 10 15 20 25 30 35 40 45 Zenith angle, θ [degree] (a) 2500 # of events 2000 1500 1000 500 0 0 50 100 150 200 250 300 350 Azimuthal angle, φ [degree] (b) Figure 4.5: Thrown zenith angle (θ) distribution and thrown azimuthal angle (φ) distribution. (a) Thrown zenith angle (θ) distribution. Blue histogram shows the actual thrown distribution and Red line shows the theoretical expectation. (b) Thrown azimuthal angle (φ) distribution. Blue histogram shows the actual thrown distribution and Red line shows the theoretical expectation. 61 4.2.2 Resolution To determine the resolution, the Monte Carlo events were processed through exactly the same reconstruction programs as the data. We compared reconstructed values to thrown values to determine how well Monte Carlo events are reconstructed. We examined three reconstructed distributions: energy, distance (impact parameter), and ψ angle (angle in the shower detector plane) using all events that passed quality cuts. Figure 4.6a shows the energy resolution histogram by evaluating ln(Erec /Ethr ). It shows the mean value is centered at 0, which means no energy bias is observed. σ shows the spread of the distribution is about 8.6%. Figure 4.6b shows the energy resolution as a function of energy to determine if the energy resolution depends on the energy. It shows that its distribution is quite flat over energies, which means that there is no energy bias as a function of energy. Figure 4.7a shows the Rp resolution histogram by evaluating ln(Rp rec /Rp thr ). It shows the mean value is centered at 0, which means no Rp bias is observed and σ shows the spread of the distribution is about 4.1%. Figure 4.7b shows the Rp resolution as a function of energy to determine if the Rp resolution depends on the energy. As in energy resolution, Rp shows no bias as a function of energy. Figure 4.8a shows the ψ resolution histogram by evaluating (ψrec − ψthr ). It shows the mean value is centered at 0.7◦ , which means that there is only a small ψ bias observed and σ shows the spread of the distribution is about 3.8◦ . For the usual monocular fluorescence detection, the ψ resolution tends to be 5◦ at best. This means that our MC set has really good resolution in terms of the geometry of the shower. Figure 4.8b shows the ψ resolution as a function of energy to determine if the ψ resolution depends on the energy. It shows no ψ bias as a function of energy. All of these tests show that Monte Carlo simulation reconstructs reasonably some important parameters comparing to thrown parameters and that it can be reliably reconstructed. 4.2.3 Missing energy We examined about one hundred thousand CORSIKA events of pure proton and iron primaries to calculate the missing energy in a given thrown total energy of showers. Since we stored thrown Giasser-Hillas parameters in the shower library, we calculated the calorimetric energy by evaluating the number of charged particles from the Giasser-Hillas parameterization and the mean energy deposition from Nerling et al. (2006). When we corrected the missing energy, we used the mixed composition weighting of events since we do not know 62 Entries 120 2941 Mean 0.006083 ± 0.001917 RMS 100 χ2 0.08591 ± 0.001355 / ndf 30.66 / 40 # of events 80 60 40 20 0 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 ln(Erec/Ethr) (a) 1 Entries 2941 0.8 χ2 / ndf 14.36 / 9 0.6 p0 0.006924 ± 0.002615 ln(Erec/Ethr) 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 16 16.5 17 17.5 18 18.5 19 log (E/eV) 10 (b) Figure 4.6: The energy resolution histogram and the energy resolution as a function of energy. (a) The energy resolution histogram. The red line shows a Gaussian fit to the histogram. (b) The energy resolution as a function of energy. The red line shows constant fit to the profile. 63 Entries 250 2941 Mean 0.006873 ± 0.0009156 RMS 0.04103 ± 0.0006475 χ2 / ndf 57.56 / 23 # of events 200 150 100 50 0 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 ln(Rp rec/R p thr) (a) 1 Entries 2941 0.8 χ2 / ndf 3.489 / 9 0.6 p0 0.006548 ± 0.001205 ln(Rp rec/R p thr) 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 16 16.5 17 17.5 18 18.5 19 log (E/eV) 10 (b) Figure 4.7: The Rp resolution histogram and the Rp resolution as a function of energy. (a) The Rp resolution histogram. The red line shows a Gaussian fit to the histogram. (b) The Rp resolution as a function of energy. The red line shows constant fit to the profile. 64 Entries Mean RMS χ2 / ndf 250 2941 0.794 ± 0.08612 3.859 ± 0.0609 57.61 / 21 # of events 200 150 100 50 0 -80 -60 -40 -20 ψ 0 rec -ψ thr 20 40 60 80 [degree] ψ rec - ψ thr [degree] (a) 10 Entries 2941 8 χ2 / ndf 5.469 / 9 6 p0 0.7488 ± 0.1166 4 2 0 -2 -4 -6 -8 -10 16 16.5 17 17.5 18 18.5 19 log (E/eV) 10 (b) Figure 4.8: The ψ resolution histogram and the ψ resolution as a function of energy. (a) The ψ resolution histogram. The red line shows a Gaussian fit to the histogram. (b) The ψ resolution as a function of energy. The red line shows constant fit to the profile. 65 the primary type in real data. The detailed application of missing energy is described in Chapter 5. CHAPTER 5 EVENT RECONSTRUCTION AND SELECTION Fluorescence detectors are operated on clear, moonless nights with at least 3 hours of no sun, no moon. The length of the dark periods depends upon the season of the year. They last for about 2 weeks in summer and 3 weeks in winter. The shortest data taking time is 3 hours per night while, in winter, the longest data taking time can be about 12 hours per night. The detectors do not observe the primary cosmic ray particle itself, but view the development of an air shower via its isotropic emission of fluorescence light. This light arrives at our detector and is collected by telescope mirrors that redirect it and focus it onto the PMTs. In data analysis, we retrace that trajectory to determine the geometry of the air shower. With this knowledge, we can then estimate the energy of the primary particle by integrating the number of secondary shower particles. We will discuss how data is stored in Section 5.1 and the methods used for event reconstruction in Section 5.2. In order to obtain good energy resolution and to properly compare data to Monte Carlo, well-reconstructed events must be identified. The criteria used for selecting such events (quality cuts) and the set of fully reconstructed cosmic ray events will be discussed in Section 5.3. Examples of successfully reconstructed data events by both Middle Drum and TALE telescopes will be shown in Section 5.4. 5.1 Data Format We employ a Data Summary Tape (or Data Storage Tape, DST) storage system. A DST file can contain multiple events. The information for each event is organized and stored in multiple DST banks, including information for which PMTs and mirrors are triggered by the event. Each event is itself delineated by a START and a STOP DST bank, in order to identify individual events at each step during reconstruction. Each DST bank is similar to a Fortran common block or a C structure, defining the type for each variable and function. In every processing step, new DST banks are created to store the results of that step. These 67 new DST banks are added alongside the existing collection of DST banks from the previous processing step in the DST file. All event information is stored in DST banks and the same structure is used for both data and Monte Carlo events. The different DST banks used in different reconstruction stages are discussed below: • HRAW1 stores the signal and timing information of all PMTs that observe an event by the Middle Drum telescopes. All mirror and tube information from the mirrors or tubes triggered, calibrated photon count from the roving xenon flasher, number of telescopes, and total number of tubes for an event are recorded. • FRAW1 stores the signal and timing information of all tubes that observe an event by the TALE telescopes. All triggered mirrors and their 320 channels are recorded. Channels 1 to 256 represent the signals of each of the PMTs. Channels 257 to 272 are high-gain summations for columns 1 to 16 and channels 273 to 288 are low-gain summations for columns 1 to 16. Channels 289 to 304 are high-gain summations for rows 1 to 16 and channels 305 to 320 are low-gain summations for rows 1 to 16. • FTRG1 stores the state of the mirrors and a trigger code. It also stores row or column pattern coincidences per event. • FPHO1 is similar to FRAW1 DST bank – it stores the tube signal and timing for TALE. It includes tube timing information and tube signal in number of photoelectrons with the pedestal subtracted. • FSCN1 contains the information from the FADC scan. It includes telescope and tube numbers along with tube flags that tell if tubes are good, rejected by the scan, or have no signal. In addition, the pedestal and number of time slices in the pulse are stored. • HBAR stores tube signal and timing with calibration information. It reads in the Roving Xenon Flasher coefficients and adds them to the data stream. The number of photoelectrons is calculated by taking the QDC and gain from calibration resources. Mirror reflectivity, a quantum efficiency factor for each tube at 337 nm, and the expected transmission value from UV filter curve for each tube are stored as well. • STPS2 is called the “Filtering DST bank”. It includes variables used to calculate the Rayleigh vector. It determines if an event is noise or an actual cosmic ray event. The Rayleigh vector will be further discussed in Section 5.2.4. A small magnitude for the Rayleigh vector indicates a random or noise event while a large magnitude indicates a real cosmic ray event. 68 • STPLN is filled after fitting the shower detector plane during the 4th reconstruction step (Pass 3). It includes the plane normal vector, track length, time difference between last and first good tubes in microseconds, and labeling for good and bad tubes. • TL4RGF is filled after the TALE and MD 4-ring monocular geometry fit. It contains the basic event geometry such as date, time, zenith angle, azimuthal angle, core position, angle in the shower detector plane, distance of closest approach, track length, crossing time for each site, and total number of photoelectrons. It also contains tube-based information such as the tube’s pointing direction, azimuth, and elevation; number of photoelectrons; timing; edge indicator; and flags indicating the quality of the time and plane fits. • HCTIM is filled after the geometry fit during the 5th reconstruction step (Pass 4). It contains a description of shower track geometry. It includes tube and telescope information based on the time fit, and it is similar to the TL4RGF DST bank. It will be updated again after the profile fit during the 6th reconstruction step (Pass 5). • PRFC is filled after the shower profile fit and it contains results of the fitted shower profiles. It contains the four shower profile parameters of the Gaisser-Hillas function : Nmax , Xmax , X0 , and λ. From these parameters, the shower energy may be estimated by integrating the Giasser-Hillas function and Nerling mean energy loss (dE/dx) function. The profile of light along the shower axis has 4 components including scinitillation, direct (unscattered) Cherenkov light, and Cherenkov light scattered by molecules and aerosols. An indicator flag for the success or failure of the profile fit is also stored. • HCBIN is also filled after the shower profile fit. It contains the binned light flux as a function of shower depth in units of photoelectrons/m2 /degree. • MC04 contains the output of the Monte Carlo simulation. For example, energy; GiasserHillas parameters; primary particle; geometry parameters such as track direction, zenith angle, R~p , shower core, Rp , and ψ; and number of sites in the detector that had mirrors and tubes triggered. 5.2 Reconstruction Programs In general, 6 Pass reconstruction chains are used (from Pass 0 to Pass 5), plus a preprocessing stage. In the preprocessing stage, we check if the operation log file exists, make 69 a cut on the data for good weather, and check if the calibration data file exists. After the preprocessing stage, the first reconstruction chain is Pass 0, which uses time matching to recognize events and extract signal information. In Pass 1, calibration information is added to signal information. In Pass 2, pattern recognition is used to distinguish between a random event and a cosmic ray event. In Pass 3, using only cosmic ray events, this step applies the shower detector plane fit over PMT clusters. In Pass 4, tube timing and tube viewing angle are used to determine shower geometry. In Pass 5, we perform a profile fit on the flux of light as a function of the depth in order to estimate the energy of the cosmic ray shower. In this dissertation, we have added a 7th reconstruction chain to the data pipeline. This Pass 6 is designed to reduce the energy bias after Pass 5. During Pass 5, we found that the energy estimated after the missing energy correction showed about a 10% bias. This bias appears to be due to poorly estimated Giasser-Hillas parameters, which are fixed by our reconstruction program in Pass 5. We examined about one hundred thousand CORSIKA events with Giasser-Hillas parameters and found the variations of the Giasser-Hillas parameters as a function of energy. In Pass 6, we take these variations into account when we redo the energy estimate with the profile fit from the first energy estimation, thus it results in no energy bias. A detailed method of Pass 6 will be discussed in Section 5.2.8. We now turn to discussing the details of each pass’s reconstruction method. 5.2.1 Preprocessing pass Before considering event reconstruction, we gather all information needed for the reconstruction itself. This is a preparation stage where we check log and calibration files, determine good-weather data, and calculate the detector ontime. The details are discussed below. • Log: We find an operation log file that is created once per night. The log file summarizes important information such as how many data files are created that night, a list of individual data parts with weather codes collected, a list of calibration and pedestal files generated, how many hours of good data files, and the average weather code for an overall night. • Good weather: We create two lists of data files. One lists all data files regardless of weather codes. The second lists only data files with good weather. For one night, 70 multiple data files can be generated. Each file contains about one hour of data and is supposed to have a weather code containing the conditions at that time. Weather codes are composed of 7 digits. The first four digits represent observation of horizon clouds in the north, east, south, and west, respectively. If no clouds are on the horizon, it is set as 0 and if clouds exist on that horizon, it is set to 1. The next two digits indicate overhead cloud and overhead cloud thickness, respectively. The overhead cloud flag ranges from 0 to 4 according to how overcast the sky is. If it is clear, it is set to 0. If clouds cover less than 20%, it is set to 1. If clouds cover less than 50%, it is set to 2. If clouds cover less than 75%, it is set to 3. If clouds cover more than 75%, it is set to 4. Cloud thickness is set to 0 if stars are visible through overhead clouds, otherwise it is set to 1. The final digit indicates haziness and it ranges from 0 to 2. If visibility is good (i.e. the Milky Way and city lights are clearly visible), it is set to 0. If lights from a city appear to be spread out and fainter, the air may be hazy, and it is set to 1. If it is ambiguous, it is set to 2. This weather is observed by the operators each hour and the code is entered by the Middle Drum operators. The good weather selection follows. First of all, we consider overhead clouds. As long as the overhead is clear, we consider it to be good weather. If the overhead is less than 50% clear, then we begin to consider the effect of horizon clouds. If no clouds on the east and south exist, then we consider it to be good weather. If there are clouds in both the east and south, then overhead coverage should be less than 20% to be considered to be good weather. It these conditions are not met, events are not processed for further analysis. • Ontime: We calculate how long telescopes are on during an observation. We summarize how long telescopes are operating with no issues with average trigger rate and weather codes. It is not always the same for each telescope, because telescopes can be offline or not run for some reason during a night’s observation. This is taken into account in the ontime calculation. • Pedestal: We calculate pedestal values for all PMTs of all telescopes. • RXF: We calculate each tube’s gain is calculated with the effective area for all telescopes based on the UV filter curve and quantum efficiency curve applied. Roving Xenon Flasher calibration data is taken at the beginning of every monthly data 71 collection period for the Middle Drum and TALE. For analysis, we use the closest Roving Xenon Flasher calibration results with respect to the data we analyze. • Calibration: We collect gain and pedestal values for all PMTs of all telescopes based on Roving Xenon Flasher calibration. 5.2.2 Pass 0 Pass 0 is a reconstruction stage that does time matching between event packets and builds events. The Central Timing computer (CT) sends the time and the event packets for each triggered telescope. Time packets are recorded by GPS time stamps and event packets are recorded by telescope trigger timing. All time and events packets are recorded in raw data files. In Pass 0 program, we match triggered events from individual telescopes using GPS timing within 100 µs time window as one event. Some events can be single telescope events while some may be multitelescope events. Event information is written in the HRAW1 bank for Middle Drum by the program HMA and in the FRAW1 bank for TALE by the program TLFDP0. The “H” stands for “HiRes I” and “MA” stands for “match”. The “TLFD” stands for “TALE FD” and “P0” stands for “Pass 0”. Hence, the HRAW1, FRAW1, FPHO1, and FSCN1 DST banks are filled. • HMA performs time matching of event and time packets (raw data) for Middle Drum • TLFDP0 performs in the same way as the HMA program, but this performs event formation for TALE. 5.2.3 Pass 1 Pass 1 is a reconstruction stage that applies calibration information to the raw data using the program HPASS1. Raw TDC values are converted into a time in units of microseconds. Time delays in the cable between mirrors and the Central Timing computer are taken into account. Raw QDC values are converted into the number of photoelectrons per mm2 by HCAL and HPED programs. The gain and pedestal values are obtained from Roving Xenon Flasher (RXF) calibration for Middle Drum and from UVLED calibration for TALE. RXF calibration is collected once a month, at the beginning of the Middle Drum and TALE FD run periods. Meanwhile, UVLED calibration is collected twice a day, at the beginning and end of each Middle Drum and TALE FD run. In general, UVLED calibration information from the beginning of the FD run is used. However, if there are inoperative telescopes or some telescope operation issues that occur during UVLED calibration, the 72 calibration from the end of the FD run is used. The HRAW1 DST bank is updated and the HBAR, MDWEAT, and TLWEAT DST banks are filled at this time. • HPASS1 applies RXF calibration to calibrate the Pass 0 data and writes it to the HBAR bank. • HCAL represents an electronics calibration bank and contains the TDC’s and QDC’s gains and pedestals. • HPED is also an electronics calibration bank and contains the QDC’s pedestals for each PMT. • MDWEAT is the DST bank that contains the weather codes and the data file number information of Middle Drum station. • TLWEAT is the DST bank that stores the weather codes and the data file number information of TALE station. 5.2.4 Pass 2 Pass 2 is a reconstruction stage that removes noise triggered events by applying the Rayleigh filtering done by the STPS2 program. In addition to real cosmic ray events, airplane blinker lights, flashers, or fluctuations in the night sky can cause the telescope to trigger and record data. We can remove these false events by applying spatial pattern recognition to the PMT cluster, which is called a pattern recognition method. Examples of a real cosmic ray event and a false event are shown in Figures 5.1a and 5.1b. The pattern recognition method is based on a Rayleigh probability density function, which is given by r2 r (5.1) f (r) = 2 e− 2σ2 , σ where r is the magnitude of the Rayleigh vector ~r and σ is N/2. Here, N is the total number of steps that corresponds to the number of adjacent pairs with tube’s spatial separation less than 1.5◦ and time separation less than 14 ns. This Rayleigh distribution is often observed when the overall magnitude of a vector is related to its directional components. In order to determine if the event is the result of random noise or a real cosmic ray, we introduce probability of Rayleigh distribution. The probability that a magnitude of a vector ~r is greater than a certain magnitude, Z ∞R, is R2 P (f ; r > R) = f (r) dr = e− 2σ2 . R (5.2) 73 (a) (b) Figure 5.1: Examples of a track-like event and a random noise event. (a) An example of a track-like event. Colour scheme indicates tube’s timings and size of filled circle is proportional to tube’s signal sizes. The event has a track-like or line-like structure and we note the time order of the PMT hits. That is to say there is a general correlation of position in the track with time. (b) An example of random noise event. It may be an air plane event. All tubes treated as bad tubes appear black in the event display. In this case, all PMTs have the same time. 74 From that, Plog variable is introduced and it is a negative log10 probability that indicates if an event is noise. It is calculated as, R2 , (5.3) N ln 10 implies a random event, while a real cosmic ray event should have a large Plog = and a small Plog value of Plog . Events with greater than 2 in Plog are kept and saved for further analysis. This cut ensures a false event rate of 1% or less. In addition, the projection of Poynting vector direction of light onto the surface of a PhotoMultiplier Tube (PMT) cluster gives the direction of an event: upward-going, downward-going, or horizontal events. Most cosmic ray events have a downward-going direction, while laser events for atmospheric calibration have an upward-going direction. All results are filled in the STPS2 DST bank. 5.2.5 Pass 3 Pass 3 is a reconstruction stage that performs a shower detector plane (SDP) fit via the program STPLN. The track across the PMT cluster of triggered tubes represents the air shower axis. Triggered tubes are a combination of noise tubes and actual signal tubes that need to be segregated. To determine a shower detector plane, we need to determine one line and a point: a line that corresponds to the triggered tube’s pointing direction and a point corresponding to the detector position. The geometry of the shower detector plane is shown in Figure 5.2. Plane fitting can be done by minimizing a χ2 function under two assumptions. The first is that the air shower is treated as a pencil beam (i.e. line source), which means it has no lateral distribution. The second is that all telescopes are located at one point. The χ2 function to minimize is X (n̂ · n̂i ) × wi , (5.4) χ2 = σi2 i where n̂ is a shower detector unit normal, n̂i represents the pointing unit vector of tube i, wi indicates the number of photoelectrons seen by tube i, and σi indicates tube angular uncertainty related to the field of view of the PMT, which is to be set as a constant, σ = 1◦ , for all tubes. Through plane fitting, we filter noise tubes out by two standards. First, tubes are treated as bad tubes or noise tubes if tubes are spatially uncorrelated with a track. In addition, if the time of a tube’s triggering is considerably different from the majority of tubes in a track, it is also considered a noise tube. “Considerably different” means that the tube’s time residual is at least 2 µs, or that spatial correlation of triggering happens 5◦ from a plane. This procedure is iterative and it continues until no further tubes are removed. Once 75 Figure 5.2: Shower detector plane (SDP) with the shower axis and the detector location. Taken from Baltrusaitis et al. (1985). 76 a shower detector plane is identified, the normal vector of the shower detector plane, track length, crossing time, and inverse angular speed are calculated. An example is shown in Figure 5.3. All results are filled in the STPLN DST bank. 5.2.6 Pass 4 Pass 4 is a reconstruction stage that performs the geometry fit by the TLMNGEOM program. It requires the HRAW1 DST bank from Middle Drum, the FRAW1 DST bank from TALE, and the STPLN DST bank to proceed. Figure 5.4 shows a diagram of the geometry. We observe tube timing along the shower axis and know each tube’s pointing direction. By fitting the tube viewing angles for the events with the arrival time of light at tubes from the shower, accurate geometrical reconstruction can be calculated. We call this method a time versus angle fit. The Middle Drum and TALE data and Monte Carlo simulations are reconstructed in monocular mode with the geometry determined by the equation π − ψ − χ Rp i , (5.5) · tan ti = tR p + c 2 where ti and χi are the trigger time and pointing direction of tube i, respectively; ψ is the in-plane angle; Rp is the impact parameter of the shower; c is the speed of the light; and tRp is the time when the shower is calculated to be at Rp . In order to determine the best geometry from combined tube measurements, we perform Figure 5.3: An event display shown with the fitted shower detector plane indicated by the black line. 77 A θi χi Rp ψ C B Figure 5.4: Geometry of an extensive air shower. Line AB represents the shower axis and point C represents the detector location. Line BC defines the ground. The angle χi is a measurable variable. The angle θi indicates light emission angle and is related to the other angles as π − ψ − χi . From the geometry fit, Rp and ψ are extracted. least squares fitting and minimize π − ψ − χ o2 X 1n Rp i t − t + χ2 = , (5.6) · tan i 0 2 c 2 σ i i where σi indicates tube signal uncertainty. By time versus angle fit, the three variables, Rp , tRp , and ψ, are determined. They are calculated using information from each telescope station, respectively, and using information from both stations. We use the fit variables from the station that has the longest view of the track. An example of a geometry fit is shown in Figure 5.5. All results are filled in the TL4RGF and the HCTIM DST banks. 5.2.7 Pass 5 Pass 5 is the stage where the shower profile is fit by the program STPFL12. Given the track geometry, we then compare the observed tube signal with the simulated signal calculated based on a trial set of shower parameters. The shower profile fit is done by minimizing χ2 : χ2 = X 1 observed − Si predicted ), 2 (Si σ i i (5.7) where Si observed is an observed tube signal, Si predicted is a predicted tube signal, and σi is an expected variance in the number of sky noise or sky background photoelectrons in the 5.6 µs integration window. An example of shower profile is shown in Figure 5.6. 78 Figure 5.5: An example of time versus angle fit. The data measurement from each PMT is shown as points with error bars. The best fit to the data is shown as the line. Tubes with signals that are significantly off of this line are considered noise tubes and are discarded. Only good tubes are used for a time fit. From this fit, Rp , tRp , and ψ are determined. Figure 5.6: An example of shower profile fit. Black points are data points and the blue line represents a profile fit based on total signal. The red line represents the scintillation contribution of the shower, the green line represents the Rayleigh Cherenkov scattering contribution of the shower, the magenta line represents the aerosol Cherenkov scattering contribution of the shower, and the grey line represents the direct Cherenkov contribution of the shower. 79 The profile of the shower is fit using the Gaisser-Hillas parameterization formula (Gaisser and Hillas 1977) Xmax −X0 λ Xmax −x x − X0 ·e λ , (5.8) Ne (x) = Nmax · Xmax − X0 where Ne (x) is the number of charged particles at a given slant depth, x, in g/cm2 ; Xmax is the slant depth where the number of secondary particles reaches the maximum; Nmax is the maximum number of particles at Xmax ; X0 is a fit parameter associated with the depth of the first interaction; and λ is a fit parameter explaining the width of the shower profile. During fitting, we fix two parameters, X0 = −100 g/cm2 and λ = 70 g/cm2 , while we vary Xmax and Nmax parameters. At this point, the impact parameter, Rp , can be varied slightly for a given ψ in order to find a better profile fit, with new geometry variables updated in the HCTIM DST bank. This process iterates until χ2 is minimized. Additionally, a parameterization of the mean energy deposit (Nerling et al. 2006) is applied: c1 + c4 + c5 · s, (5.9) (c2 + s)c3 where ci (i =1,2,3,4, and 5) are the constants from Nerling et al. (2006). The mean ionization αeff (s) = N · loss rate, αeff , at a given shower age, s, is determined in the condition at Ecut = 1 MeV. However, the value of Ecut in our program is 50 keV, so N indicates a scale normalization constant to correct for the difference in Ecut . By integrating a combined of Gaisser-Hillas parameterization formula and mean energy deposition parameterization, we estimate the calorimetric energy of the shower, which is the visible portion of the energy from total energy of the shower. All results are filled in the PRFC and the HCBIN DST banks. 5.2.8 Pass 6 Pass 6 is a new additional reconstruction stage that fits the shower profile performed by the TL4RPF program. It is a similar reconstruction stage to the Pass 5 STPFL12 program except that a X0 , X1 , and λ look-up table is used. Given the first energy estimation from Pass 5, we use a look-up table of X0 , X1 , λ to take proper shower parameters into account. After that, shower parameters are weighted based on mixed composition as measured by the HiRes/MIA and HiRes experiments. From that, it re-calculates a calorimetric energy of the shower. This reconstruction stage is designed to reduce the energy bias after Pass 5. During Pass 5, we found that the energy estimated after the missing energy correction showed about a 10% bias, apparently due to wrongly estimated Giasser-Hillas parameters, which are fixed 80 by our reconstruction program in Pass 5. We examined about one hundred thousand CORSIKA events of pure proton and iron primaries to look at Giasser-Hillas parameters. We found variations of the Giasser-Hillas parameters as a function of energies for each primary. It appears to be different from the values that we used in Pass 5 as fixed values. Therefore, in Pass 6, we take these variations into account when we redo the energy estimate with the profile fit from the first energy estimation. These variations are shown in Figure 5.7a, 5.7b, and 5.8. By employing variations of shower parameters, we were able to correct the energy bias. 5.2.9 Postprocessing pass From Pass 6, we obtained the calorimetric energy of the shower, which is visible energy calculated from the fluorescence lights detected by telescopes. Some portion of energy, is however, carried by muons and neutrinos, and this causes our detector to under-estimate the energy. Here, we will correct for the missing energy through studying Monte Carlo simulations. We examined more than one hundred thousand CORSIKA events consisting of pure proton shower or pure iron shower to calculate the calorimetric energy of the shower for a thrown total shower energy. We used the same method as before to calculate the calorimetric energy of the showers, using the Giasser-Hillas parameters plus mean energy deposition parameterization. Figure 5.9 shows the fractional calorimetric energy of showers induced by proton and iron primary cosmic rays. As we see, proton primary cosmic rays result in less missing energy than an iron primary. This is because air showers initiated by iron primary cosmic rays are likely to have more muons and neutrinos than showers initiated by protons. As the energy of primary cosmic rays increases, the required missing energy correction decreases because the number of interactions from secondary particles increases. As the number of interactions increases, greater proportions of the shower energy are transferred to the electromagnetic component of the shower. This results in more fluorescence light that is visible to our telescopes and is therefore not missing energy. The total missing energy correction is carried out based on the evolving mixed composition measurement of the HiRes/MIA and HiRes experiments. Given the calorimetric energy of a shower, we estimate the proton fraction and weigh the missing energy correction by the proton and iron fractions at that energy. We then calculate the total energy of the shower. The accuracy of reconstructing the shower energy from Monte Carlo simulations was discussed previously in Chapter 4. 81 Entries χ2 / ndf 100 52381 Entries 51536 255.5 / 42 χ2 / ndf 63.16 / 42 p0 −460.8 ± 57.5 p0 90 p1 106.9 ± 9.3 p1 271.8 ± 19.2 80 p2 −6.623 ± 0.505 p2 −15.92 ± 1.05 p3 0.1289 ± 0.0091 p3 0.3069 ± 0.0190 −1464 ± 117.5 70 60 λ 50 40 30 20 10 0 16 16.5 17 17.5 18 18.5 log (E/eV) 19 19.5 20 20.5 10 (a) 200 100 Entries 52479 Entries 51011 χ2 / ndf 159.9 / 41 χ2 / ndf 60.57 / 41 p0 1.143e+05 ± 7.076e+03 p0 1.339e+05 ± 1.909e+04 p1 −2.532e+04 ± 1.553e+03 p1 −3.026e+04 ± 4.195e+03 p2 2099 ± 127.5 p2 2566 ± 345.0 p3 −77.06 ± 4.64 p3 −96.65 ± 12.58 p4 1.056 ± 0.063 p4 1.363 ± 0.172 X0 0 −100 −200 −300 −400 16 16.5 17 17.5 18 18.5 log (E/eV) 19 19.5 20 20.5 10 (b) Figure 5.7: The λ variation as a function of energy and the X0 variation as a function of energy from the CORSIKA simulation. (a) The variation in λ as a function of energy from the CORSIKA program. The red indicates a proton primary and the blue indicates an iron primary. (b) The X0 variation as a function of energy from the CORSIKA simulation. Again, red (blue) represents the proton (iron) simulation. 82 100 90 80 Entries 52479 Entries 51888 χ2 / ndf 34.45 / 44 χ2 / ndf 42.53 / 44 p0 18.21 ± 0.66 p0 110.7 ± 2.8 p1 −0.4274 ± 0.0357 p1 − 3.612 ± 0.149 70 X1 60 50 40 30 20 10 0 16 16.5 17 17.5 18 18.5 log (E/eV) 19 19.5 20 20.5 10 Figure 5.8: The variation in X1 as a function of energy from the CORSIKA simulation. Red (Blue) represents the proton (iron) simulation. 1 0.95 Ecal/Etot 0.9 0.85 0.8 0.75 0.7 15 16 17 Entries 52479 Entries 51891 χ 2 / ndf 338.6 / 42 χ 2 / ndf 440.1 / 42 p0 3.577 ± 1.056 p0 p1 − 1.18 ± 0.23 p1 3.58 ± 0.82 p2 0.133 ± 0.019 p2 − 0.2786 ± 0.0675 p3 − 0.005822 ± 0.000690 p3 0.009775 ± 0.002454 p4 8.947e− 05 ± 9.391e− 06 p4 − 0.0001299 ± 0.0000334 18 log (Ecal/eV) 19 − 16.62 ± 3.76 20 21 10 Figure 5.9: CORSIKA simulation of the fractional calorimetric energy as a function of the total shower energy. Red indicates the case of a proton primary while blue indicates an iron primary cosmic ray. 83 5.3 Event Selection We analyzed about 8 months of data collected between June 2014 and January 2015. The Middle Drum and TALE data and Monte Carlo events are processed using the same programs. Some events collected can be seen in the Middle Drum station data only, some in the TALE station data only, or some in both. To obtain good resolution and Data/Monte Carlo comparisons, quality cuts are performed on the fully reconstructed showers. These are described in Table 5.1 and are optimized for the fluorescence events observed by Middle Drum and TALE stations. They differ for each event set: TALE, MD, and 4-ring (TALE and MD). An event is retained if: • An event has a good weather flag. • An event has the following geometry conditions. The reconstructed energy should be greater than 1016.5 eV because Monte Carlo events are thrown from 1016 eV. The zenith angle, θ, is < 70◦ and the in-plane angle, ψ, is < 120◦ . Error in the in-plane angle, dψ, depends on the event set and picking events with longer tracks helps us to improve resolution. This is why additional track length and number of good tubes cuts are not needed. The duration of the event depends on the event set because TALE telescopes are designed for nearby events and low energy events and Middle Drum telescopes are designed for far-away events and high energy events. • An event has been successfully reconstructed. Additionally, Xmax bracketing cut is introduced and implies that Xmax should be between the first and last profile bins. In other words, Xmax should be in the field of view of detectors. The slant depth of the first triggered tube must fall greater than 200 g/cm2 and the profile reduced χ2 < 10. • An event has a good profile shape. A clear peak in the profile should be observed and the maximum signal of the tube should be greater than the certain flux threshold, introduced by Sfraction and Speak variables. These two cuts improve in energy resolution. The case that if you observe only the rising or falling parts of air showers in the field of view, it means that Xmax is located out of our field of view and it is hard to trust Xmax fitted by the Gaisser-Hillas function. It may give us inaccurate Xmax , leading to the incorrect energy estimation. Obtaining events with well-defined triangular shape of profiles is the critical approach to achieving reliable energy estimation. Table 5.2 shows the number of data events that are retained after each cut level. As quality cuts are optimized for fluorescence events, a number of events seen by the TALE 84 Table 5.1: The cut levels are described for three different data sets: TALE, MD, and 4-ring (TALE and MD). For some variables, the same cuts are satisfied for all three data sets, but for the other variables, cuts are assigned differently for different data sets because of characteristic data set. The TALE telescopes are designed for nearby events and low energy events and Middle Drum telescopes are designed for far-away events and high energy events. The 4-ring (TALE and MD) can cover the full energy range. TALE Cut level 1 Cut Cut Cut Cut Cut Cut Cut level level level level level level level 2 3 4 5 6 7 8 Cut level 9 Cut level 10 Cut level 11 Cut level 12 Cut level 13 MD 4-ring (TALE and MD) Down-going events only Both MD and TALE running Good weather cut θ < 70◦ Erec > 1016.5 eV no track length cut crossing time > 2 µs crossing time > 6 µs crossing time > 4 µs no number of good tubes cut An event must be successfully reconstructed (profile fit) χ2 /ndf of profile fit < 10 (Xmax bracketing) Xfirst < Xmax < Xlast Xfirst > 200 g/cm2 Xlast -Xfirst > 150 g/cm2 Sfraction > 1.05 Sfraction > 1.05 Sfraction > 1.20 Speak > 50.0 Speak > 80.0 Speak > 50.0 Speak × Rp2 > 3.2 · 108 dψ < 20◦ dψ < 12◦ dψ < 15◦ ψ < 120◦ 85 Table 5.2: The number of data events passed through each cut level for three different data sets: TALE, MD, and 4-ring (TALE and MD). Among retained events, events seen by both Middle Drum and TALE detectors are dominant. Cut level 1 Cut level 2 Cut level 3 Cut level 4 Cut level 5 Cut level 6 Cut level 7 Cut level 8 Cut level 9 Cut level 10 Cut level 11 Cut level 12 Cut level 13 TALE 112642 43810 30439 428 375 350 144 83 71 65 65 MD 9608 5099 4754 1291 1270 1223 935 682 139 123 118 4-ring (TALE and MD) 4034 3177 2805 1113 1109 1084 886 492 461 461 455 telescope only are cut out because the TALE-only set should have events with short tracks and Cherenkov-dominated events. At the final cut level, events seen by Middle Drum and TALE stations are both dominant. 5.4 Example Event An example of the data event that triggered both the Middle Drum and TALE telescopes is shown in Figures 5.10, 5.11, and 5.12. Figure 5.10 shows an event display with the fitted shower detector plane. Figure 5.11 shows time versus angle fit for the same event. Figure 5.12 shows the profile fit to estimate the calorimetric energy of the shower for the same event. This event has about 12 µs long-duration with more than 50◦ of track length. From the time versus angle fit, its in-plane angle, ψ, is determined as 101.1◦ , which means that the direction of this event is coming to detectors at 6 km distance. From the profile fit, we can clearly see Xmax in the field of view with the triangular structure. The scintillation light is the dominant signal and energy of this event is estimated as 1018.56 eV. This is a typical 4-ring event used in this analysis. 86 Figure 5.10: Event display of an example event observed by both the Middle Drum and TALE telescopes. The line fit for the shower detector plane (black line) is overlaid. The event is from 2014/06/05 07:38:59 UTC. 87 Figure 5.11: The time versus angle fit for the same event as Figure 5.10. The time versus angle fit is shown as the line. 88 Figure 5.12: The profile fit for the same event as shown in Figure 5.10. One can see the overall fit to the measured light (blue) as well as the various contributions to the signal from fluorescence, direct Cherenkov, scattered Cherenkov, etc. CHAPTER 6 DATA/MONTE CARLO COMPARISONS We generate a Monte Carlo simulation using the previously measured energy spectra and compositions as input to calculate the aperture of the detector for understanding the acceptance of the detector. We analyze data and Monte Carlo events using the same reconstruction programs and check if observable parameters’ distributions agree with each other. Comparing distributions from data and the Monte Carlo simulation is an important step to understand the detector. The following pages include figures with two panels; the top panels show each observable parameter’s data (hdt) and Monte Carlo (hmc) distributions, and the bottom panels show the ratio of the data to Monte Carlo (hrat) events. The number of Monte Carlo events is larger than the number of data; thus, we normalized area under the curve of the Monte Carlo histogram to that of the data histogram’s area to properly compare the two. Figures 6.1-6.8 show comparisons for the geometric variables: track length, crossing time, zenith angle (θ), in-plane angle (ψ), azimuthal angle (φ), impact parameter (Rp ), the number of good tubes per degree, and the number of participating rings. In Figure 6.1, there are three peaks present and they indicate two-telescope, three-telescope, and four-telescope events. The ratio of the data to Monte Carlo for the track length does not show these peaks. It appears relatively flat, which means that two distributions agree reasonably well with each other. In Figure 6.2, a couple of small crossing time bins do not agree; however in most of the bins, two distributions agree with each other. In Figure 6.3, there is a small shift in the Monte Carlo with respect to the data. This results in a ratio plot with a slope. However, in given statistics, two distributions are reasonable matched. In Figures 6.4, 6.5, and 6.6, the data and Monte Carlo distributions for the ψ, φ, and Rp agree well. In Figure 6.7, the number of good tubes per degree means how many good tubes exist relative to the length of the shower track. It is a measure of how wide or how near the event is. The data and Monte Carlo distributions agree. In Figure 6.8, the number of participating rings indicates what kinds of events are triggered in terms of the ring concept. The first bin indicates 1 ring participated in the event, the second bin indicates 2 rings participated in the event, 90 hmc Entries Mean RMS Integral 45 40 hdt 35 Entries Mean RMS Integral 30 # of events 2942 45.19 12.69 637 25 638 43.58 13.45 637 20 15 10 5 0 0 10 20 30 40 50 60 70 80 90 Track length [degree] hrat Entries Mean RMS Integral χ 2 / ndf const slope 2 1.8 1.6 62 28.03 19.2 174.2 42.39 / 60 0.9438 ± 0.1882 -0.002007 ± 0.003988 Data/MC 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0 10 20 30 40 50 60 70 80 90 Track length [degree] Figure 6.1: The Data/Monte Carlo comparison for the event track length distribution. In the top panel, data (hdt) is shown by the black squares and the Monte Carlo distribution (hmc) is shown by the red histogram. There are three peaks corresponding to 2, 3, and 4 telescope events. In the bottom panel, the ratio of the data to Monte Carlo (hrat) is shown along with the linear best fit line. 91 hmc Entries Mean RMS Integral 80 70 hdt Entries Mean RMS Integral 60 # of events 2942 9.436 4.171 548 50 638 9.072 4.103 548 40 30 20 10 0 0 2 4 6 8 10 12 14 16 18 20 Crossing time [µ s] hrat Entries Mean RMS Integral χ2 / ndf const slope 2 1.8 1.6 18 8.675 5.602 23.14 8.047 / 16 1.126 ± 0.131 -0.01784 ± 0.01205 Data/MC 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0 2 4 6 8 10 12 14 16 18 20 Crossing time [µ s] Figure 6.2: The Data/Monte Carlo comparison for the event crossing time distribution. In the top panel, data (hdt) is shown by the black squares and the Monte Carlo distribution (hmc) is shown by the red histogram. In the bottom panel, the ratio of the data to Monte Carlo (hrat) is shown along with the linear best fit line. 92 hmc Entries Mean RMS Integral 45 40 hdt Entries Mean RMS Integral 35 # of events 2942 33.58 12.8 638 30 25 638 31.39 14.09 638 20 15 10 5 0 0 10 20 30 40 50 60 70 80 90 Zenith angle, θ [degree] hrat Entries Mean RMS Integral χ 2 / ndf const slope 2 1.8 1.6 35 32.47 23.71 47.21 21.75 / 33 1.264 ± 0.151 -0.01045 ± 0.00398 Data/MC 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0 10 20 30 40 50 60 70 80 90 Zenith angle, θ [degree] Figure 6.3: The Data/Monte Carlo comparison for the event zenith angle (θ) distribution. In the top panel, data (hdt) is shown by the black squares and the Monte Carlo distribution (hmc) is shown by the red histogram. In the bottom panel, the ratio of the data to Monte Carlo (hrat) is shown along with the linear best fit line. There is a small offset between the data and the Monte Carlo, which results in a slope in the ratio. For this analysis, which has limited statistics, the agreement between data and Monte Carlo simulation is good enough. However, when a larger data sample is available, it would be good to better understand and fix this. 93 hmc Entries Mean RMS Integral 70 60 hdt Entries Mean RMS Integral 50 # of events 2942 82.37 20.02 638 40 638 81.57 20.18 638 30 20 10 0 0 20 40 60 80 100 120 140 160 180 In-plane angle, ψ [degree] hrat 2 64.93 Mean 1.8 RMS 29.71 Integral 22.64 χ2 / ndf 1.6 1.4 Data/MC 19 Entries 15.72 / 17 const 0.9641 ± 0.2022 slope -0.0002733 ± 0.0023752 1.2 1 0.8 0.6 0.4 0.2 0 0 20 40 60 80 100 120 140 160 180 In-plane angle, ψ [degree] Figure 6.4: The Data/Monte Carlo comparison for the event ψ angle distribution. In the top panel, data (hdt) is shown by the black squares and the Monte Carlo distribution (hmc) is shown by the red histogram. In the bottom panel, the ratio of the data to Monte Carlo (hrat) is shown along with the linear best fit line. 94 hmc Entries Mean RMS Integral 35 2942 189.4 107.4 638 hdt # of events 30 Entries Mean RMS Integral 25 638 178.7 107.2 638 20 15 10 5 0 50 100 150 200 250 300 350 Azimuthal angle, φ [degree] hrat Entries Mean RMS Integral χ 2 / ndf const slope 2 1.8 1.6 36 168.7 102.7 36.84 23.93 / 34 1.038 ± 0.102 -0.000662 ± 0.000451 Data/MC 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0 50 100 150 200 250 300 350 Azimuthal angle, φ [degree] Figure 6.5: The Data/Monte Carlo comparison for the event azimuthal angle (φ) distribution. In the top panel, data (hdt) is shown by the black squares and the Monte Carlo distribution (hmc) is shown by the red histogram. In the bottom panel, the ratio of the data to Monte Carlo (hrat) is shown along with the linear best fit line. 95 hmc Entries Mean RMS Integral 120 hdt 100 # of events 2942 5652 2781 620 Entries Mean RMS Integral 80 638 5694 2781 620 60 40 20 0 0 2000 4000 6000 8000 10000 12000 14000 Impact parameter, Rp [m] hrat Entries Mean RMS Integral χ2 / ndf const slope 2 1.8 1.6 14 7977 4019 14.51 7.884 / 12 0.9708 ± 0.1124 5.182e-08 ± 1.788e-05 Data/MC 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0 2000 4000 6000 8000 10000 12000 14000 Impact parameter, Rp [m] Figure 6.6: The Data/Monte Carlo comparison for the event impact parameter, Rp , distribution. In the top panel, data (hdt) is shown by the black squares and the Monte Carlo distribution (hmc) is shown by the red histogram. The mean impact parameter for this data set is about 6 km. In the bottom panel, the ratio of the data to Monte Carlo (hrat) is shown along with the linear best fit line. 96 hmc Entries Mean RMS Integral 80 70 hdt Entries 638 Mean 1.765 RMS 0.4674 Integral 638 60 # of events 2942 1.785 0.479 638 50 40 30 20 10 0 0 1 2 3 4 5 6 Number of good tubes per degree hrat Entries Mean RMS Integral χ2 / ndf const slope 2 1.8 1.6 Data/MC 1.4 31 2.566 1.179 44.49 18.09 / 29 1.17 ± 0.21 -0.1286 ± 0.1125 1.2 1 0.8 0.6 0.4 0.2 0 0 1 2 3 4 5 6 Number of good tubes per degree Figure 6.7: The Data/Monte Carlo comparison for the number of good tubes per degree of track length. In the top panel, data (hdt) is shown by the black squares and the Monte Carlo distribution (hmc) is shown by the red histogram. In the bottom panel, the ratio of the data to Monte Carlo (hrat) is shown along with the linear best fit line. 97 250 # of events 200 hmc Entries Mean RMS Integral 150 100 hdt Entries Mean RMS Integral 50 0 0.5 2942 3.011 0.8179 638 1 1.5 2 2.5 3 3.5 4 638 2.969 0.8541 638 4.5 Number of participating rings hrat Entries Mean RMS Integral χ2 / ndf const slope 2 1.8 1.6 Data/MC 1.4 4 2.276 1.157 4.629 2.04 / 2 1.133 ± 0.203 -0.04674 ± 0.06453 1.2 1 0.8 0.6 0.4 0.2 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Number of participating rings Figure 6.8: The Data/Monte Carlo comparison for the number of participating rings in an event. In the top panel, data (hdt) is shown by the black squares and the Monte Carlo distribution (hmc) is shown by the red histogram. In the bottom panel, the ratio of the data to Monte Carlo (hrat) is shown along with the linear best fit line. 98 the third bin indicates 3 rings participated in the event, and the forth bin indicates a 4-ring event. Figure 6.9 shows a comparison for the total reconstructed event energy with the missing energy correction applied. Two distributions agree well. Figure 6.10 shows a comparison for the shower detector plane angle. We note that there is a dip structure around 0◦ seen in the Monte Carlo simulation due to the triggering requirement of our detector, which requires multiple columns. This makes the telescope less sensitive to vertical (or horizontal) showers. The Monte Carlo simulation does a reasonable job of modeling the data, given the statistics. However, one can note that the data and the MC do not agree particularly well on the insensitivity in the trigger in the central region. A future study will require additional work to resolve this. Figures 6.11 and 6.12 show comparisons for the shower core position in the CLF coordinates. The data and Monte Carlo distributions agree well. 99 hmc Entries Mean RMS Integral 80 70 2942 17.89 0.3862 638 hdt Entries 638 Mean 17.89 RMS 0.4096 Integral 638 # of events 60 50 40 30 20 10 0 16 16.5 17 17.5 18 18.5 19 19.5 20 Total energy [log10E/eV] hrat Entries Mean RMS Integral χ2 / ndf const slope 2 1.8 1.6 Data/MC 1.4 20 18.34 0.6494 26.09 21.28 / 18 7.912 ± 1.832 -0.3908 ± 0.1018 1.2 1 0.8 0.6 0.4 0.2 0 16 16.5 17 17.5 18 18.5 19 19.5 20 Total energy [log10E/eV] Figure 6.9: The Data/Monte Carlo comparison for total energy distribution. In the top panel, data (hdt) is shown by the black squares and the Monte Carlo distribution (hmc) is shown by the red histogram. The mean energy of this data set is 1017.9 eV. Again, there is a small offset between data and Monte Carlo resulting in a slope in the ratio. It has little impact on this analysis. In the bottom panel, the ratio of the data to Monte Carlo (hrat) is shown along with the linear best fit line. 100 hmc Entries Mean RMS Integral 50 2942 1.099 29.52 638 hdt Entries 638 Mean 0.7638 RMS 27.28 Integral 638 # of events 40 30 20 10 0 -80 -60 -40 -20 0 20 40 60 80 SDP angle [degree] hrat Entries Mean RMS Integral χ 2 / ndf const slope 2 1.8 1.6 Data/MC 1.4 25 4.016 32.97 24.65 40.66 / 23 0.8643 ± 0.0467 -0.001124 ± 0.001404 1.2 1 0.8 0.6 0.4 0.2 0 -80 -60 -40 -20 0 20 40 60 80 SDP angle [degree] Figure 6.10: The Data/Monte Carlo comparison for the shower detector plane angle. In the top panel, data (hdt) is shown by the black squares and the Monte Carlo distribution (hmc) is shown by the red histogram. We note that there is a dip structure around 0◦ seen in the Monte Carlo due to the triggering requirement of our detector, which requires multiple columns. This makes the telescope less sensitive to vertical (or horizontal) showers. The Monte Carlo does a reasonable job of modeling the data, given the statistics. However, one can note that the data and the MC do not agree particularly well on the insensitivity in the trigger in the central region. A future study will require addition work to resolve this. In the bottom panel, the ratio of the data to Monte Carlo (hrat) is shown along with the linear best fit line. 101 hmc Entries Mean RMS Integral 120 hdt 100 # of events 2942 -5346 4460 638 Entries Mean RMS Integral 80 638 -5234 4853 638 60 40 20 0 -30000 -20000 -10000 0 10000 20000 30000 Xcore [m] hrat Entries Mean RMS Integral χ2 / ndf const slope 2 1.8 1.6 Data/MC 1.4 29 -9750 1.418e+04 64.88 18.09 / 27 0.9305 ± 0.0808 -1.194e-06 ± 1.195e-05 1.2 1 0.8 0.6 0.4 0.2 0 -30000 -20000 -10000 0 10000 20000 30000 Xcore [m] Figure 6.11: The Data/Monte Carlo comparison for Xcore position. In the top panel, data (hdt) is shown by the black squares and the Monte Carlo distribution (hmc) is shown by the red histogram. Xcore is in CLF centered coordinates. In the bottom panel, the ratio of the data to Monte Carlo (hrat) is shown along with the linear best fit line. 102 hmc Entries 2942 Mean 1.465e+04 RMS 3984 Integral 638 120 hdt # of events 100 Entries 638 Mean 1.421e+04 RMS 4614 Integral 638 80 60 40 20 0 -30000 -20000 -10000 0 10000 20000 30000 Ycore [m] hrat Entries Mean RMS Integral χ2 / ndf const slope 2 1.8 1.6 23 -2480 7695 76.73 13.42 / 21 1.243 ± 0.215 -1.964e-05 ± 1.382e-05 Data/MC 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -30000 -20000 -10000 0 10000 20000 30000 Ycore [m] Figure 6.12: The Data/Monte Carlo comparison for Ycore position. In the top panel, data (hdt) is shown by the black squares and the Monte Carlo distribution (hmc) is shown by the red histogram. Ycore is in CLF centered coordinates. In the bottom panel, the ratio of the data to Monte Carlo (hrat) is shown along with the linear best fit line. CHAPTER 7 ENERGY SPECTRUM MEASUREMENT As seen in the previous chapter, Data and Monte Carlo distributions agree with each other. This means we understand the cosmic rays and our detector response reasonably well. Thus, we are ready to calculate the energy spectrum using events observed by Middle Drum and TALE telescopes. We begin with reviewing how the cosmic ray energy spectrum is calculated in Section 7.1. We move onto the aperture calculation using the Monte Carlo simulation, detector ontime, and exposure in Sections 7.1.1, 7.1.2, and 7.1.3. Finally, we show a measured cosmic ray energy spectrum firstly by itself and then in comparison to the other experiments in Section 7.2. 7.1 Spectrum Calculation The cosmic ray energy spectrum can be calculated according to how many particles are observed per unit area, per solid angle, per unit time, and per unit energy. This can be written as N (Ei ) , (7.1) AΩ(Ei ) × T × ∆Ei where J(Ei ) is the flux of cosmic rays, N (Ei ) is the number of observed cosmic rays, J(Ei ) = AΩ(Ei ) is the energy-dependent aperture calculated from the Monte Carlo simulation, T is the detector ontime, and ∆Ei is the width of the energy bin for the i th energy bin, Ei . After the quality cuts, we calculate the number of data events retained. The histogram of the number of data events is shown in Figure 7.1 and the number of data events per energy bin is recorded in Table 7.1. Each energy bin size is 0.1 in log scale. Error estimates of the number of data events are calculated using the paper by Feldman and Cousins for small event count (Feldman and Cousins 1998). 104 Figure 7.1: The histogram of the number of data events. Each energy bin size is 0.1 in log scale. These data events are passed after quality cuts and are used in the energy spectrum calculation. The uncertainties on the points are by Feldman and Cousins (1998). Table 7.1: The number of data events per each energy bin. These data events are used in calculating the energy spectrum. log10 (E/eV) 17.25 17.35 17.45 17.55 17.65 17.75 17.85 17.95 18.05 18.15 18.25 18.35 18.45 18.55 18.65 18.75 18.85 Number of events 12 33 52 67 72 71 70 53 33 43 28 26 17 20 15 5 8 105 7.1.1 Aperture For Fluorescence Detector (FD) measurements, the aperture depends significantly on energy. Since higher energy events are brighter, they can be seen farther away. On the other hand, lower energy events are only observed relatively nearby. The aperture and acceptance of the detector must be calculated by the Monte Carlo (MC) technique. The aperture can be written as NReconstructed (E) × A0 Ω0 , (7.2) NThrown (E) where NReconstructed (E) is the number of reconstructed MC events that are reconstructed AΩ = and passed by the quality cuts, NThrown (E) is the number of thrown MC events, A0 indicates the effective area that is determined by the thrown maximum impact parameter, Rp max , and Ω0 indicates the solid angle that is determined by the thrown maximum zenith angle, θ. Here, NReconstructed (E) NThrown (E) represents the acceptance. The calculated aperture is shown in Figure 7.2. The effective area is given by A0 = π Rp2 where Rp max max − Rp2 min , is the maximum of impact parameter and is set to be 50,000 m and Rp (7.3) min is the minimum of impact parameter and is defined as 10 m. The solid angle is given by Ω0 = 2π{1 − cos(θmax )}, (7.4) where θmax is the maximum of zenith angle and is defined as 80◦ . 7.1.2 Detector ontime Ontime is calculated from the time each telescope has triggers permitted during the run. Each data file records how long each telescope is active and this is summed for each run period. For this analysis, we use TALE data epoch II, which collects data observed during this period: June 2014 through January 2015. We only consider data collected when both Middle Drum and TALE runs where operated for observations. We examine all data files and sum up the time that mirror triggers are permitted, resulting in 537.189 hours. This detector ontime, however, includes good weather as well as bad weather. In our quality cuts, we only select data events that passed good weather conditions, thus this weather cut effect is taken into account in the detector ontime calculation. After the weather cut, the detector ontime is reduced to 470.102 hours. 106 Figure 7.2: The calculated aperture as a function of the energy. The aperture is energydependent. It tells what the acceptance of the detector is as the energy of cosmic ray particles increase. We can see that, as expected, the aperture starts small at low energies and grows. At the mean energy of this data set, 1017.9 eV, the aperture is about 4×107 m2 sr. 107 7.1.3 Exposure The exposure is defined as the product of the aperture and the detector ontime and is given by NReconstructed (E) × A0 Ω0 T, (7.5) NThrown (E) where NReconstructed (E) is the number of reconstructed MC events that are passed by quality AΩT = cuts, NThrown (E) is the number of thrown MC events, A0 indicates the effective area in m2 , Ω0 indicates the solid angle in sr, and T is the detector ontime in sec. The calculated exposure is shown in Figure 7.3. 7.2 Energy Spectrum The 4-ring Middle Drum and TALE monocular energy spectrum was calculated using about 8 months of data from June 2014 through January 2015. By looking at three different data sets, TALE, MD, and 4-ring (TALE and MD), events seen by Middle Drum and TALE stations are both dominant and the quality cuts are optimized for fluorescence events. Thus, Figure 7.3: The exposure as a function of energy. It is a product of the aperture and detector ontime. 108 events below 1017 eV are mostly cut out. After cuts, 638 events are retained and used for energy spectrum calculation. The Monte Carlo simulation is thrown up to 1019 eV, so that the measurement is able to explore the region between 1017 and 1019 eV. The flux of cosmic rays is calculated using Equation 7.1 and is shown in Figure 7.4. Next, the cosmic ray flux is multiplied by energy to the third power, E 3 J(E) in order to take out the dominant slope and better show the fine structure of the spectrum. This is shown in Figure 7.5. Figure 7.6 shows the cosmic ray energy spectrum fit to a single slope power law, which is overlaid. It shows -3.3±0.06 spectral index between the second knee and ankle features. This spectral index agrees with the other Telescope Array recent results. Figure 7.7 shows the cosmic ray energy spectrum observed using this analysis compared with other TA energy spectra. One spectrum is the monocular energy spectrum using only the TALE detector, but it includes Cherenkov-dominated events as well as fluorescence events. Another one is the energy spectrum calculated by the Telescope Array surface detectors using 9 years of data. This 4-ring Middle Drum and TALE monocular energy spectrum plays the role of a bridge between the TALE-only detector and surface detector spectrums. The 4-ring Middle Drum and TALE monocular energy spectrum agrees with the high-end of the TALE-only detector and the low-end of the surface detector spectra. Figure 7.8 shows the cosmic ray energy spectrum observed using this analysis comparing the other external experiments: HiRes-2, Icetop, KASCADE-Grande, and Auger experiments. The 4-ring monocular energy spectrum and HiRes-2 agree well in the most of the energy range. The low-end of the 4-ring monocular energy spectrum and Icetop & KASCADE-Grande agree, but around 1018 eV, their statistics become poor. 109 Figure 7.4: The flux J(E) of cosmic ray measured by the Middle Drum and TALE stations. One sees the strong E3 dependance. 110 Figure 7.5: The energy spectrum measured by the Middle Drum and TALE stations. Here the flux has been multiplied by E3 to take out the dominant slope and better show the fine structure of the spectrum. Figure 7.6: The energy spectrum measured by the Middle Drum and TALE stations with a single slope of the power law. The energy spectrum is fit to a single slope of the power law. The spectral index is -3.3 between the second knee (E ∼ 1017.05 eV) and ankle (E ∼ 1018.7 eV). 111 Figure 7.7: The energy spectrum measured by the Middle Drum and TALE stations comparing to the other recent TA energy spectra. The red square points represent the TALE monocular energy spectrum and the green square points represent the TA surface energy spectrum. The black points represent this bridge spectrum measurement. 112 Figure 7.8: The energy spectrum measured by the Middle Drum and TALE stations (black circles) comparing to the other external experiments: HiRes-2 (blue open circles), Icetop (green open circles), KASCADE-Grande (magenta open circles), and Auger (red open circles) experiments. CHAPTER 8 CONCLUSIONS Ultra-High Energy Cosmic Rays (UHECRs) are known as the most energetic particles in the Universe. However, the sources of UHECRs are yet unknown. We speculate that they come from violent extragalactic astrophysical objects. Thus, the aim of UHECR studies is to understand where they come from, what their chemical composition is, how they are accelerated to such high energies, and how they propagate through the intergalactic medium. The work in this dissertation is focused on particles in the energy range from 1017 to 1018 eV, where the transition from galactic to extragalactic sources is suspected to occur. The measurement of the cosmic ray energy spectrum may give a hint of where UHECRs come from. This is because spectral features may indicate the energy distributions of particles at astrophysical sources and their distribution depends on the surrounding environment near sources. Moreover, features in this energy range may imply a difference in chemical composition between galactic and extragalactic particles. The heavy galactic component of cosmic rays appears to be dying off in this energy range while the light extragalactic component of cosmic rays begins entering the Milky Way. We use TALE data epoch II, which was collected between June 2014 through January 2015, when both the Middle Drum and TALE FDs were operated. By looking at three different data sets, TALE, MD, and 4-ring (TALE and MD), we found that events seen by both Middle Drum and TALE stations are the dominant data in this transition energy region. The quality cuts were used to select well-reconstructed events and were optimized for fluorescence events. Thus, events below 1017 eV are mostly cut out due to contamination by Cherenkov light. After cuts, there are 638 events retained and used for a monocular energy spectrum calculation. For Fluorescence Detector (FD) measurements, the aperture depends significantly on energy. Since higher energy events are brighter, they can be seen farther away. On the other hand, lower energy events are only seen nearby. The aperture and acceptance of the detector must be calculated by the Monte Carlo (MC) technique. We create a complete simulation of the shower and detector to generate MC samples, putting previously measured 114 energy spectra and composition into the MC simulation. Next, we wrote the MC events out in the same format as actual data and analyzed them with the same exact programs as the data. Finally, we compared data with MC and showed that the distributions agree, which tells us that we understand the detector. We then used this information to calculate an aperture and an exposure, which were necessary to determine the cosmic ray energy spectrum. In conclusion, the combined monocular MD and TALE 4-ring detector has observed a cosmic ray spectrum giving a spectral index, -3.3±0.06, in the energy range 1017 -1019 eV. We have demonstrated that the combined MD and TALE 4-ring energy spectrum plays an important role as a bridge between the TALE detector and surface detector spectra. The energy spectrum using only the TALE detector was determined for E < 1018.3 eV. Meanwhile, the energy spectrum using the surface detector was determined for E > 1018.3 eV. Thus, a single energy bin is overlapped using these two techniques. The 4-ring analysis is the optimal technique for studying this energy region. By combining the data from the MD and TALE telescopes, a monocular energy spectrum can connect the other TA techniques and shows good agreement with the same overall spectral slope. More statistics are required to examine the energy spectrum features. The TALE fluorescence data was initially collected while the detector trigger was still being developed. The data analyzed is called TALE data epoch I. TALE data epoch II began when stable running commenced, and ended when another trigger improvement was made. We have chosen to examine the earliest stable running to measure the cosmic ray spectrum in order to work out the logistics and to complete an analysis in a reasonable amount of time. However, much more data (about a factor of five) now exists. By adding the new data, the energy spectrum measurement by this techniques can be extended to both higher and lower energies in addition to reducing the error bars by a bit more than a factor of two. This might enable observation of the second knee and ankle features. However, there were trigger changes and other changes to the detector that will require a significant amount of analysis and simulation work. Therefore, we leave that for a future project. REFERENCES Abbasi, R. U., et al. 2008, Phys. Rev. Lett., 100, 101101 —. 2014, The Astrophysical Journal Letters, 790, L21 —. 2015, Phys. Rev. D, 92, 032007 Abbasi, R. U., et al. 2017, Astroparticle Physics, 86, 21 Abraham, J., et al. 2007, Science, 318, 938 Abu-Zayyad, T., et al. 2012a, Astroparticle Physics, 39-40, 109 —. 2012b, The Astrophysical Journal, 757, 26 —. 2012c, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 689, 87 Allen, M. G. 2012, PhD thesis, University of Utah Baltrusaitis, R. M., et al. 1985, Nuclear Instruments and Methods in Physics Research A, 240, 410 Bird, D. J., et al. 1993, Phys. Rev. Lett., 71, 3401 Feldman, G. J., and Cousins, R. D. 1998, Phys. Rev. D, 57, 3873 Gaisser, T. K., Engel, R., and Resconi, E. 2016, Cosmic Rays and Particle Physics (UK: Cambridge University Press) Gaisser, T. K., and Hillas, A. M. 1977, International Cosmic Ray Conference, 8, 353 Greisen, K. 1966, Phys. Rev. Lett., 16, 748 Hanlon, W. F. 2008, PhD thesis, University of Utah Heitler, W. 1954, Quantum theory of radiation (Oxford: Clarendon) Ivanov, D. 2016, PoS, ICRC2015, 349 Kakimoto, F., Loh, E., Nagano, M., Okuno, H., Teshima, M., and Ueno, S. 1996, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 372, 527 Matthews, J. 2005, Astroparticle Physics, 22, 387 Nerling, F., Bluemer, J., Engel, R., and Risse, M. 2006, Astropart. Phys., 24, 421 Patrignani, C., et al. 2016, Chin. Phys., C40, 100001 Tameda, Y., et al. 2009, Nuclear Instruments and Methods in Physics Research A, 609, 227 116 Watson, A. A., and Pierre Auger Collaboration. 2008, Nuclear Instruments and Methods in Physics Research A, 588, 221 Zatsepin, G. T., and Kuz’min, V. A. 1966, Soviet Journal of Experimental and Theoretical Physics Letters, 4, 78 Zundel, Z. J. 2016, PhD thesis, University of Utah