THE FERMI LARGE AREA TELESCOPE ON ORBIT: EVENT CLASSIFICATION, INSTRUMENT RESPONSE FUNCTIONS, AND CALIBRATION

The TKR is the section of the LAT where γ rays ideally convert to e⁺e⁻ pairs and their trajectories are measured. A full description of the TKR can be found in Atwood et al. (2007, 2009). A simplified schematic of the TKR is shown in Figure 2. Starting from the top (farthest from the CAL), the first 12 paired layers are arranged to immediately follow converter foils, which are composed of ∼3% of a radiation length of tungsten. Minimizing the separation of the converter foils from the following SSD planes, and hence the lever arm between the conversion point and the first position measurements, is critical to minimize the effects of multiple scattering. This section of the TKR is referred to as the thin or front section. The next four layers are similar except that the tungsten converters are ∼6 times thicker; these layers are referred to as the thick or back section. The last two layers have no converter; this is dictated by the TKR trigger, which requires hits in three x–y paired adjacent layers (see Section 3.1.1) and is therefore insensitive to γ rays that convert in the last two layers.

Zoom In Zoom Out Reset image size

Thus, the TKR effectively divides into two distinct instruments with notable differences in performance, especially with respect to the PSF and background contamination. This choice was suggested by the need to balance two basic (and somewhat conflicting) requirements: simultaneously obtaining good angular resolution and a large conversion probability. The tungsten foils were designed such that there are approximately the same number of γ rays (integrated over the instrument FoV) converted in the thin and thick sections. In addition to these considerations, experience on-orbit has also revealed that the aggregate of the thick layers (∼0.8 radiation lengths) limits the amount of back-scattered particles from the CAL returning into the TKR and ACD in high-energy events (i.e., the CAL backsplash) and reduces tails of showers in the TKR from events entering the back of the CAL. These two effects help to decrease the background contamination in front-converting events.

After three years of on-orbit experience with the TKR we can now assess the validity of our design decisions. The choice of the solid-state TKR technology has resulted in negligible down time and extremely stable operation, minimizing the necessity for calibrations. Furthermore, the very high signal-to-noise ratio of the TKR analog readout electronics has resulted in a single hit efficiency, averaged over the active silicon surface, greater than 99.8%, with a typical noise occupancy smaller than 10⁻⁵ for a single readout channel. (We note for completeness that the fraction of non-active area presented by the TKR is ∼11% at normal incidence.) As discussed below, this has yielded extremely high efficiency for finding tracks and has been key to providing the information necessary to reject backgrounds.

The efficiency and noise occupancy of the TKR over the first three years of operation are shown in Figure 3. The variations in the average single-strip noise occupancy are dominated by one or a few noisy strips, which have been disabled at different times during the mission. The baseline of 4 × 10⁻⁶ is dominated by accidental coincidences between event readouts and charged particle tracks (see below and Section 2.1.4) and corresponds to an upper limit of ∼3 noise hits per event in the full LAT on average. Since these noise hits are distributed across 16 towers and 36 layers per tower, their effect on the event reconstruction is insignificant.

Zoom In Zoom Out Reset image size

The TKR readout is digital, i.e., the readout is binary, with a single threshold discriminator for each channel, and no pulse height information is collected at the strip level. The individual electronic chains connected to each SSD strip consist of a charge-sensitive preamplifier followed by a simple CR–RC shaper with a peaking time of ∼1.5 μs. Due to the implementation details, the baseline restoration tends to be current-limited, and the signal at the output of the shaper is far from the exponential decay characteristic of a linear network, with the discriminated output being high for ∼10 μs for minimum ionizing particles (MIPs) at the nominal ∼1/4 MIP threshold setting. As a consequence, the latched TKR strip signals are typically present for that amount of time after the passage of an MIP. If the LAT is triggered within this time window, these latent signals will be read out and become part of the TKR event data. The rate of occurrence of these ghost signals (which may result in additional tracks when the events are reconstructed) depends directly on the charged particle background rate. Mitigation against this contamination in the data is discussed in Sections 2.1.4, 3, and 5.2.

Each detector subsystem contributes one or more trigger primitive signals that the LAT trigger system uses to adjudicate whether to read out the LAT detectors (see Section 3.1). The TKR trigger is a coincidence of three x–y paired adjacent layers within a single tower (hence a total of six consecutive SSD detector planes). Due primarily to power constraints, the TKR electronic system does not feature a dedicated fast signal for trigger purposes. Rather, the logical OR of the discriminated strip signals from each detector plane is used to initiate a non-retriggerable one-shot pulse of fixed length (typically 32 clock ticks or 1.6 μs) that is the basic building block of the three-in-a-row trigger primitive (see Section 3.1.1). In addition, the length (or time over threshold) of this layer-OR signal is recorded and included in the data stream. Since the time over threshold depends on the magnitude of the ionization in the SSD, which in turn depends on the characteristics of the ionizing particle, it provides useful information for the background rejection stage.

The efficiency of the three-in-a-row trigger is >99%. This is due in part to the redundancy of this trigger for the vast majority of events (i.e., by passing through many layers of Si, most events have multiple opportunities to form a three-in-a-row). The trigger efficiency, measured in flight using MIP tracks crossing multiple towers, is shown to be very stable in Figure 4.

Zoom In Zoom Out Reset image size

Perhaps the most important figure of merit for the TKR is the resulting PSF for reconstructed γ-ray directions. At low energy, the PSF is dominated by multiple scattering, primarily within the tungsten conversion foils (tungsten accounts for ∼67% of the material in the thin section and ∼92% of the material in the thick section). At high energy, the combination of the strip pitch of 228 μm, the spacing of the planes, and the overall height of the TKR result in a limiting precision for the average conversion of ∼01 at normal incidence. MC simulations predict that the transition to this measurement-precision-dominated regime should occur between ∼3 GeV and ∼20 GeV. We will discuss the PSF in significantly more detail in Section 6.

Details of the CAL can be found in Atwood et al. (2009) and Grove & Johnson (2010). Here, we highlight some key aspects for understanding the LAT performance. The CAL is a three-dimensional (3D) imaging calorimeter. This is achieved by arranging the CsI crystals in each tower module in 8 layers, each with 12 crystal logs (with dimensions 326 mm × 26.7 mm × 19.9 mm) spanning the width of the layer. The logs of alternating layers are rotated by 90° about the LAT boresight, and aligned with the x- and y-axes of the LAT coordinate system.

Each log is read out with four photodiodes, two at each end: a large photodiode covering low energies (<1 GeV per crystal), and a small photodiode covering high energies (<70 GeV per crystal). Each photodiode is connected to a charge-sensitive preamplifier whose output drives (1) a slow (∼3.5 μs peaking time) shaping amplifier for spectroscopy and (2) a fast shaping amplifier (∼0.5 μs peaking time) for trigger discrimination. In addition, the output of each slow shaper is connected to two separate track-and-hold stages with different gains (× 1 and × 8).

The outputs of the four track-and-hold stages are multiplexed to a single analog-to-digital converter. The four gain ranges (two photodiodes × two track-and-hold gains) span an enormous dynamic range, from <2 MeV to 70 GeV deposited per crystal, which is necessary to cover the full energy range of the LAT. A zero-suppression discriminator on each crystal reduces the CAL data by eliminating the signals from all crystals with energies <2 MeV. To minimize CAL data volume, each log end reports only a single range, the lowest gain unsaturated range (the one-range, best-range readout) for most events. For likely heavy ions, each log end reports all four ranges (the four-range readout) for calibrating the energy scale across the different ranges (see Section 3.1.1 for details of how the readout mode is selected).

The CAL provides two inputs to the global LAT trigger. At each log end the output of each fast shaper (for both the large and the small diode) is compared to an adjustable threshold by a discriminator to form two separate trigger request signals. In the standard science configuration, the discriminator thresholds are set at 100 MeV and 1 GeV energy deposition. The 1 GeV threshold is >90% efficient for incident γ rays above 20 GeV.

For each log with deposited energy, two position coordinates are derived simply from the geometrical location of the log within the CAL array, while the longitudinal position is derived from the ratio of signals at opposite ends of the log: the crystal surfaces were treated to provide monotonically decreasing scintillation light collection with increasing distance from a photodiode. Thus, the CAL provides a 3D image of the energy deposition for each event.

Since the CAL is only 8.6 radiation lengths thick at normal incidence, for energies greater than a few GeV shower leakage becomes the dominant factor limiting the energy resolution, in particular because event-to-event variations in the early shower development cause fluctuations in the leakage out the back of the CAL. Indeed, by ∼100 GeV about half of the total energy in showers at normal incidence escapes out the back of the LAT on average. The intrinsic 3D imaging capability of the CAL is key to mitigating the degradation of the energy resolution at high energy through an event-by-event 3D fit to the shower profile. This was demonstrated both in beam tests and in simulations (Baldini et al. 2007; Ackermann et al. 2010), achieving better than 10% energy resolution well past 100 GeV. The imaging capability also plays a critical role in the rejection of hadronic showers (see Section 3.3.7). Furthermore, for events depositing more than ∼1 GeV in the CAL, imaging in the CAL can be exploited to aid in the TKR reconstruction and in determining the event direction (see Section 3.2).

We have monitored the performance of the CAL continuously over its three years of operation. We observe a slow (∼1% per year) decrease in the scintillation light yield of the crystals due to radiation damage in low Earth orbit, as we anticipated prior to launch (see also Section 7.3.2). Although we do not yet correct for this decreased light yield, we have derived calibration constants on a channel-by-channel basis for the mission to date. In 2012 January we started reprocessing the full data set with these updated calibrations, and in the future expect to maintain a quarterly cadence of updates to ensure that the calibrated values do not drift by more than 0.5%.

The third LAT subsystem is the ACD, critically important for the identification of LAT-entering charged CRs. Details of its design can be found in Moiseev et al. (2007) and Atwood et al. (2009). From the experience of the LAT predecessor, the Energetic Gamma-Ray Experiment Telescope (EGRET), came the realization that a high degree of segmentation was required in order to minimize self-veto due to hard X-ray back-scattering (often referred to as backsplash) from showers in the CAL (Esposito et al. 1999).

The ACD consists of 25 scintillating plastic tiles covering the top of the instrument and 16 tiles covering each of the four sides (89 in all). The dimensions of the tiles range between 561 and 2650 cm² in geometrical surface and between 10 and 12 mm in thickness. By design, the segmentation of the ACD does not match that of the LAT tower modules, to avoid lining up gaps between tiles with gaps in the TKR and CAL. The design requirements for the ACD specified the capability to reject entering charged particles with an efficiency >99.97%. To meet the efficiency requirement, careful design of light collection using wavelength-shifting fibers, and meticulous care in the fabrication to maintain the maximum light yield were needed. The result was an average light yield of ∼23 photoelectrons (p.e.) for a normally incident MIP for each of the two redundant readouts for each tile.

The required segmentation inevitably led to less than complete hermeticity, with construction and launch survival considerations setting lower limits on the sizes of the gaps between tiles. Tiles overlap in one dimension, leaving gaps between tile rows in the other. The gaps are ∼2.5 mm and coverage of these gaps is provided by bundles of scintillating fibers (called ribbons), read out at each end by a photomultiplier tube (PMT). The light yield for these ribbons, however, is considerably less than for the tiles: it is typically ∼8 p.e./MIP and varies considerably along the length of a bundle. Therefore, along the gaps the efficiency for detecting the passage of charged particles is lower. However, the gaps comprise a small fraction of the total area (<1%) and accommodating them did not require compromising the design requirements. In addition to the gaps between ACD tile rows, the corners on the sides of the ACD have gaps that are not covered by ribbons and must be accounted for in the reconstruction and event classification (see Section 3.4).

As with the TKR and CAL, the ACD provides information used in the hardware trigger as well as in the full reconstruction of the event. The output of each PMT is connected to (1) a fast shaping amplifier (with ∼400 ns shaping time) for trigger purposes and (2) two separate slow electronics chains (with ∼4 μs shaping time and different gains) to measure the signal amplitude. The use of the fast signal in the context of the LAT global trigger will be discussed in more detail in Section 3. Although the main purpose of the ACD is to signal the passage of charged particles, this subsystem also provides pulse height information. The two independent slow signals and the accompanying circuitry for automatic range selection accommodate a large dynamic range, from ∼0.1 MIP to hundreds of MIPs.

As for both the TKR and the CAL, the performance of the ACD has been continuously monitored over the past three years. The stability of the MIP signal is shown in Figure 5 and in summary shows very little degradation. Note that slight changes in the selection criteria and spectra of the MIP calibration event sample cause small (<0.5%) variations in mean deposited energy per event.

Zoom In Zoom Out Reset image size

The LAT hardware trigger collects information from the LAT subsystems and, if certain conditions are fulfilled, initiates event readout. Because each readout cycle produces a minimum of 26.5 μs of dead time (even more if the readout buffers fill dynamically), the trigger was designed to be efficient for γ rays while keeping the total trigger rate, which is dominated by charged CRs, low enough to limit the dead-time fraction to less than about 10% (which corresponds to a readout rate of about 3.8 kHz). The triggering criteria are programmable to allow additional, prescaled event streams for continuous instrument monitoring and calibration during normal operation. We will defer discussion of the actual configuration used in standard science operations and of the corresponding performance to the more general discussion of event processing in Section 3.1.

To limit the data volume to the available telemetry bandwidth, collected data are passed to the on-board filter. The on-board filter (see Section 3.1.2) consists of a few event selection algorithms running in parallel, each independently able to accept a given event for inclusion in the data stream to be downlinked. Buffers on the input to the on-board filter can store on the order of 100 events awaiting processing. Provided that the on-board filter processes at least the average incoming data rate no additional dead time will be accrued because the on-board filter is busy. The processors used for the on-board filter must also build and compress the events for output, and the time required to make a filter decision varies widely between events. Therefore, the event rate that the on-board filter can handle depends on the number of events passing the filter. In broad terms, the processors will saturate for output rates between 1 kHz and 2.5 kHz, depending on the configuration of the on-board filter. In practice, the average output rate is about 350 Hz, and the amount of dead time introduced by the on-board filter is negligible.

Soon after launch, it became apparent that the LAT was recording events that included an unanticipated background: remnants of electronic signals from particles that traversed the LAT a few μs before the particle that triggered the event. We refer to these remnants as ghosts. An example of such an event is shown in Figure 6.

Zoom In Zoom Out Reset image size

It is important to re-emphasize a point made in the previous subsections: many of the signals that are passed to the hardware trigger and the on-board filter are generated by dedicated circuits whose shaping times are significantly shorter than for the circuits that read out the data from the same sensors for transmission to the ground. Although the signals with longer shaping times used for the ground event processing measure the energy deposited in the individual channels far more precisely, they also suffer the adverse consequence of being more susceptible to ghost particles. Table 1 gives the characteristic times for both the fast signals (i.e., those used in the trigger) and the slow signals (i.e., those used in the event-level analysis) for each of the subsystems.

The ghost events have been the largest detrimental effect observed in flight data. They affect all of the subsystems, significantly complicate the event reconstruction and analysis, and can cause serious errors in event classification.

• 1.  
  They can leave spurious tracks in the TKR that do not point in the same direction as the incident γ ray and might not be associated with the correct ACD and CAL information due to the different time constants of the subsystems.
• 2.  
  They can leave sizable signals in the CAL that can confuse the reconstruction of the electromagnetic shower, degrading the measurement of the shape and direction of the shower itself, and that can cause the energy reconstruction to misestimate the energy of the incident γ ray.
• 3.  
  They can leave significant signals in the ACD that can be much larger than ordinary backsplash from the CAL and that can cause an otherwise reconstructable γ ray to be classified as a CR.

We will discuss the effects on the LAT effective area in more detail in Section 5. In brief, above ∼1 GeV, where the average fractional ghost signal in the CAL is small, the loss of effective area is of the order of 10% or less. This loss is largely due to the fact that valid γ rays can be rejected in event reconstruction and classification if ghost energy depositions in the CAL cause a disagreement between the apparent shower directions in CAL and TKR. At lower γ-ray energies ghost signals can represent the dominant contribution to the energy deposition in the instrument, and the corresponding loss of the effective area can exceed 20%. We also emphasize that, while these figures are averaged over orbital conditions, the fraction of events that suffer from ghost signals, as well as the quantity of ghost signals present in the affected events, varies by a factor of 2–3 at different points in the orbit due to variation in the geomagnetic cutoff and resulting CR rates. Recovering the effective area loss due to ghost signals is one of the original and main motivations of ongoing improvements to the reconstruction (Rochester et al. 2010).

Finally, during extremely bright solar flares (SFs) the pileup of several >10 keV X-rays within a time interval comparable to the ACD signal shaping time can also cause γ rays to be classified as CRs (see Appendix A).

The presence of ghost signals in the LAT data reduced the efficiency of the event selections tuned on simulated data lacking this effect, so we developed a strategy to account for the ghosts. One of the triggers implemented in the LAT data acquisition system, the PERIODIC trigger (see Section 3.1.1), provides us with a sample of ghost events: twice per second the LAT subsystems are read out independently of any other trigger condition. For the CAL and the ACD, all the channels, even those below threshold, are recorded (this is not possible for the TKR). Since these triggers occur asynchronously with the particle triggers, the ghost signals collected are a representative sample of the background present in the LAT.

We merge the signals channel-by-channel from a randomly chosen periodic trigger event into each simulated event. In more detail, we analyze the signals in each periodic trigger, converting the instrument signals to physical units; analog-to-digital converter channels in the ACD and CAL, and time over threshold signals in the TKR, are converted into deposited energy. Since the intensity and spectrum of CRs seen by the LAT depend on the position of the Fermi satellite in its orbit, the ghost events are sorted by the value of McIlwain L (Walt 2005) at the point in the orbit where the event was recorded. Then, during the simulation process, for each simulated event we randomly select a periodic trigger event with similar McIlwain L, and add the energy deposits in this "overlay event" to those of the original simulated event, after which the combined event is digitized and reconstructed in the usual way.

We have used the improved simulations made with this overlay technique to more accurately evaluate the effective area of the LAT, and to better optimize the event selection criteria when developing Pass 7. Furthermore, we are studying ways to identify and tag the ghost signals and remove them from consideration in the event analysis (Rochester et al. 2010).

We have developed interfaces between the software libraries that implement the flux generation and coordinate transformation used by the ScienceTools and our detailed Geant4-based detector simulation. This allows us to re-use the same source models that we use with the ScienceTools within our detailed Geant4-based detector simulation, insuring consistent treatment of Fermi pointing history and the source modeling. One application for the present work was simulating atmospheric γ rays above 10 GeV for comparison with our Earth limb calibration sample (see Section 3.6.3).

In addition to simulating individual sources, we also simulate uniform γ-ray fields that we can use to explore the instrument response across the entire FoV and energy range of the LAT. These will henceforth be referred to as allGamma simulations, and they have the following features.

• 1.  
  The γ rays are generated with an E⁻¹ number spectrum (where E is the γ-ray energy), i.e., the same number of γ rays are generated in each logarithmic energy bin.
• 2.  
  The individual γ rays are randomly generated on a sphere with 6 m² cross-sectional area (i.e., large enough to contain the entire LAT and most of the spacecraft) centered at the origin of the LAT reference frame, i.e., the center of the TKR/CAL interface plane.
• 3.  
  The directions of the γ-rays are sampled randomly across 2π of downward-going (in the LAT reference frame) directions, so as to represent a semi-isotropic incident flux.

For defining the IRFs we simulate 200 million γ rays over the range log₁₀(E/1 MeV) ∈ [1.25, 5.75] and a further 200 million γ rays over the range log₁₀(E/1 MeV) ∈ [1.25, 2.75]. The net result is to produce a data set that adequately samples the LAT energy band and FoV and that far exceeds the statistics of any single point source.

In order to develop our event classification algorithms, and to quantify the amount of residual CR background that remains in each γ-ray event class, we require accurate models of both the CR background intensities and the interactions of those particles with the LAT. Accordingly, we use the same Geant4-based detector simulation described in the previous section to simulate fluxes of the CR backgrounds.

There are three components to the CR-induced background in the LAT energy range.

• 1.  
  Primary CRs. Protons, electrons, and heavier nuclei form the bulk of charged CRs. Due to the shielding by the magnetic field of the Earth, only particles with rigidities above a few GV (spacecraft position dependent) can reach the LAT orbit. Therefore, the background from primary CRs is relevant above a few GeV. The LAT event classification provides a very effective rejection of charged particles, up to an overall suppression factor of 10⁶ for CR protons (with a prominent contribution from the ACD). However, due to the intense flux of primary protons, a significant number of misclassified CRs enter even the P7SOURCE, P7CLEAN, and P7ULTRACLEAN γ-ray classes. Minimum-ionizing primary protons above a few GeV produce a background around 100–300 MeV in the LAT. However, protons that inelastically scatter in the passive material surrounding the ACD (e.g., the micro-meteoroid shield) without entering the ACD can produce secondary γ rays at lower energies, which are then indistinguishable from cosmic γ rays, and which we refer to as irreducible background (see also Section 4.4).
• 2.  
  Charged secondaries from CR interactions. CRs entering the atmosphere produce a number of secondaries in interactions with the atmosphere itself. Charged particles produced in these interactions can spiral back out, constrained by the magnetic field of the Earth, and eventually re-enter the atmosphere. These particles are predominantly protons, electrons, and positrons and are an important backgrounds below the geomagnetic cutoff rigidity. Due to the high efficiency of the ACD in rejecting charged particles, the dominant contribution are to this background in the P7CLEAN and P7ULTRACLEAN event classes positron annihilations in the passive material around the ACD. The resulting γ rays are again indistinguishable from cosmic γ rays.
• 3.  
  Neutral secondaries from CR interactions. γ rays and neutrons are also produced in local CR interactions. Unaffected by the magnetic field, they approach the LAT from the direction of the Earth with intensities peaking at the Earth limb. Neutrons are easily suppressed by the event classification scheme and do not significantly contribute to the background. The γ-ray background is suppressed by pointing away from the Earth and excluding events that originate from near or below the Earth limb. Due to the long tails of the LAT PSF (see Section 6), however, a small fraction of the events originating from the Earth limb are reconstructed with directions outside the exclusion regions. Since the Earth limb is extremely bright (Abdo et al. 2009c), even this small fraction is an important residual background for celestial γ rays. Since the PSF tails broaden with decreasing energy, the main contribution is at energies below a few hundred MeV (see Section 4.5). Similarly, a small fraction of γ rays entering the LAT with incidence angles greater than 90° are mischaracterized as downward-going (∼5% to ∼0.01%, depending on the energy, incidence angle, and the event sample). Because the Earth limb is extremely bright and some part of it is almost always behind the LAT, γ rays from the limb are the dominant component of "back-entering" background contamination. Furthermore, since we misestimate the directions of these back-entering γ rays by >90°, they are often reconstructed outside the Earth exclusion region.

Figure 10 shows the average CR-induced particle intensities at the orbit of the LAT in the model that we use. For comparison the intensity of the extragalactic diffuse γ-ray emission measured by the LAT (Abdo et al. 2010e) is overlaid to demonstrate the many orders of background suppression necessary to distinguish it from particle background. The model was developed prior to launch based on data from satellites in similar orbits and balloon experiments (Mizuno et al. 2004).

Zoom In Zoom Out Reset image size

As the particle rates are strongly dependent on location in geomagnetic coordinates, the details of the orbit model are also important. For tuning the event analysis, or for estimating the background rates for typical integration times of months or years, the simulated time interval must be at least equal to the precession period of the Fermi orbit (53.4 days). Simulating these high particle rates for such a long time interval is quite impractical, in terms of both CPU capacity and disk storage requirements. For studies of background rejection we usually simulate an entire precession period to ensure a proper sampling of the geomagnetic history, but to limit the particle counts we generate events for only a few seconds of simulated time every several minutes, e.g., a typical configuration requires event generation for 4 s every 4 minutes of time in orbit. This partial sampling is a compromise between the limited CPU and disk usage, and the requirement of having good statistics. Considering the LAT background rejection power, in order to have sizable statistics after even the first stages of the event analysis are performed, we must start with a simulated background data set of over 10⁹ CRs.

Each subsystem provides one or more trigger primitives (or trigger requests) as detailed in the following list.

• 1.  
  TKR (also known as "three-in-a-row"). Issued when three consecutive x–y silicon layer pairs—corresponding to six consecutive silicon planes—have a signal above threshold (nominally 0.25 MIPs). This signals the potential presence of a track in a tower. Since many tracks cross between towers and/or have more than three x–y layers within a tower, the TKR trigger request is very efficient.
• 2.  
  CAL_LO. Issued when the signal in any of the CAL crystal ends crosses the low-energy trigger threshold (nominally 100 MeV).
• 3.  
  CAL_HI. Issued when the signal in any of the CAL crystal ends crosses the high-energy trigger threshold (nominally 1 GeV).
• 4.  
  VETO. Issued when the signal in any of the ACD tiles is above the veto threshold (nominally 0.45 MIP). It signals a charged particle crossing the tile. The trigger system has a programmable list of ACD tiles associated with each TKR tower; if ACD shadowing is enabled (which is the case for the nominal configuration) an additional ROI primitive is assembled when a TKR primitive in a tower happens in coincidence with a VETO primitive in the ACD tiles associated to the same tower.
• 5.  
  CNO. Issued when the signal in any of the ACD tiles is above the CNO threshold (nominally 25 MIPs). It indicates the passage of a heavily ionizing nucleus (CNO stands for "carbon, nitrogen, oxygen") and it is primarily used for the calibration of the CAL (see Section 7.3.2).

The LAT has the ability to generate three other trigger primitives. We will refer to these as special primitives. Two of these are unused in flight and will not be discussed here. The third, the PERIODIC trigger, runs at a nominal frequency of 2 Hz during all science data taking and is used for diagnostic and calibration purposes.

We emphasize that although the trigger primitives provided by the TKR and the CAL are tower-based, a trigger initiates a readout of the whole LAT.

Trigger primitives are collected by the Central Trigger Unit. All 256 possible combinations of the eight trigger primitives are mapped into so-called trigger engines. In brief, some trigger requests are allowed to open a trigger window of fixed duration (nominally 700 ns) and the set of primitives collected in this time interval is compared to a table of allowed trigger conditions. In the case where a trigger condition is satisfied, a global trigger (or trigger accept) is issued and event acquisition commences. The LAT readout mode (i.e., enabling the CAL and ACD zero suppression and selecting the CAL one-range or four-range readout) can be individually set for each trigger engine. In addition, trigger engines are scalable: for each trigger condition a prescale is specified, corresponding to the number of valid trigger requests necessary to issue a single global trigger (obviously no prescale is applied to engines intended for γ-ray collection).

The standard trigger engine definitions are listed in Table 3. Trigger engine 7 is particularly important for γ-ray events: it requires no special primitives, a three-in-a-row signal in the TKR (TKR), an absence of ROI shadowing and CNO from the ACD, and no CAL_HI. Engine 6 handles CAL_HI events and engine 9 handles events with enough energy deposition in the CAL to potentially cause ROI shadowing.

Trigger engine 4 is typical for calibration data (e.g., heavy ion collection): it requires no special primitives, three-in-a-row in the TKR (TKR), a highly ionizing passage in the ACD (CNO) in close correspondence (ROI), and at least a moderate energy deposition in the CAL (CAL_LO). Engines 5 and 10 process very few γ rays and are used primarily for calibration and monitoring, so they are prescaled to limit the readout rate and/or to match the downlink bandwidth allocation.

As mentioned before, since limited telemetry bandwidth is available for data downlink, some event filtering is required. The on-board filter allows the coexistence of different filtering algorithms and, in fact, in the nominal science data taking configuration, all the events are presented to three different filters.

• 1.  
  The gamma filter, designed to accept γ rays (with an output average rate of approximately 350 Hz).
• 2.  
  The hip filter, designed to select heavy ion event candidates, primarily for CAL calibration (with an output average rate of approximately 10 Hz).
• 3.  
  The diagnostic filter, designed to enrich the downlinked data sample in events useful to monitor sensor performance and selection biases: specifically ∼2 Hz of periodic triggers and an unbiased sample of all trigger types prescaled by a factor 250 (with an output average rate of approximately 20 Hz).

The gamma filter consists of a hierarchical sequence of veto tests, with the least CPU-intensive tests performed first. If an event fails a test, it is marked for rejection and not passed on for further processing (however, a small subset of events that are marked for rejection are downlinked through the diagnostic filter). Some of the tests utilize rudimentary, two-dimensional tracks found by the on-board processing. The gamma filter has several steps (listed in order of processing).

• 1.  
  Reject events that have patterns of ACD tile hits that are consistent with CRs and do not have the CAL_LO trigger primitive asserted, making it unlikely that the ACD hits were caused by backsplash.
• 2.  
  Accept all events for which the total energy deposited in the CAL is greater than a programmable threshold, currently set to 20 GeV.
• 3.  
  Reject events that have ACD hit tile patterns that are spatially associated with the TKR towers that caused the trigger, provided that the energy deposited in the CAL is less than a programmable threshold (currently set to 350 MeV).
• 4.  
  Reject events for which a significant energy deposition in the CAL (typically >100 MeV) is present but the pattern of hits in the TKR is unlikely to produce any track.
• 5.  
  Reject events for which rudimentary tracks match with individual ACD tiles that were hit, provided the energy deposited in the CAL is less than some programmable amount (typically 5 GeV).
• 6.  
  Reject events that do not have at least one rudimentary track.

Although it may seem strange to apply the requirement that there be any tracks after cutting due to matches between tracks and the ACD, recall that the on-board filter software is highly optimized for speed, and terminates processing of each event as soon as it is possible to reach a decision. Thus, the testing for track matches is performed during the track-finding stage, at the time the tracks are constructed.

The starting point for the energy evaluation is the measured energy depositions in the crystals. The centroid of the energy deposition is determined and the principal axes of the shower are evaluated by means of a principal moment analysis. In the Pass 6 and Pass 7 event reconstruction procedure, the energy deposition is treated as a single quantity, with no attempt to identify contamination from ghost signals. Work to develop an algorithm to separate the CAL energy deposition into multiple clusters and to disentangle ghost signals is ongoing (Rochester et al. 2010). The amount of energy deposited in the TKR is evaluated by treating the tungsten–silicon detector as a sampling calorimeter; this contribution is an important correction at low energies.

We apply three algorithms to estimate the actual energy of an event: a parametric correction (PC), a fit of the shower profile (SP), and a maximum likelihood (LH) approach. The energy assigned to any given event is the energy from one or the other of these algorithms. The three methods were designed to provide the best performance in different parts of the energy and incidence angle phase space (in fact, the LH algorithm was only tuned up to 300 GeV, while the SP algorithm does not work well for events below ∼1 GeV). Accordingly, they provide different energy resolutions and their distributions have slightly different biases (i.e., the most probable values are slightly above or below the true energy) for different energies and incidence angles; more details can be found in Atwood et al. (2009).

In fact, the only significant change in the Pass 7 event reconstruction relative to Pass 6 is to apply separate corrections for the biases of each energy estimation algorithms. We used MC simulations to characterize the deviations of the most probable value of the energy dispersion from the true energy across the entire LAT phase space for the three methods separately. Such deviations were found to be typically of the order of a few percent (with a maximum value of ∼10%) and always significantly smaller than the energy resolution—with LH displaying a negative bias and PC and SP displaying a positive bias in most of the phase space.

We generated correction maps (as functions of γ-ray energy and zenith angle) and in Pass 7 the residual average bias for all the inputs of the final energy assignment (discussed in Section 3.3.2) is less than 1% in the entire LAT phase space.

For the TKR we merge and assemble clusters of adjacent hit strips into track candidates by combinatorial analysis. We have developed two methods: for CAL-seeded pattern recognition (CSPR) the trajectory of the original γ ray is assumed to point at the centroid of the energy released in the CAL; the blind search pattern recognition (BSPR) can be used when there is little or no energy deposit in the CAL sensitive volumes.

Both the CSPR and BSPR algorithms start by considering nearby pairs of TKR clusters in adjacent layers as candidate tracks (the CSPR algorithm limits the candidate tracks to those pointing toward the CAL energy centroid). Both algorithms then proceed by using a Kalman filtering technique (see Kalman 1960; Fruhwirth 1987, for details about Kalman filtering and its use in particle tracking) which tests the hypotheses that additional TKR clusters were generated by the same incident particle and should be associated with each of the candidate tracks, and adds any such clusters to the appropriate candidate tracks. Furthermore, as each cluster is added to candidate tracks the Kalman filter updates the estimated direction and associated covariance matrix of those tracks. Both the CSPR and BSPR algorithms are weighted to consider the best candidate track to be the one that is both pointing toward the CAL centroid and is the longest and straightest. In the CSPR, the main axis of the CAL energy deposition is also considered, candidate tracks for which the TKR and CAL estimated directions differ significantly are disfavored starting at ∼1 GeV, and increasingly so at higher energies. At the completion of the CSPR algorithm the best candidate track is selected and confirmed as a track, and the clusters in it are flagged as used. We iterate the CSPR algorithm until no further tracks can be assembled from the unused TKR clusters, then proceed with the BSPR.

If more than one track is found in a given event, we apply a vertexing algorithm that attempts to compute the most likely common origination point of the two highest quality (i.e., longest and straightest) tracks, and, more importantly, to use that point as a constraint in combining the momenta of the two tracks to obtain a better estimate of the direction of the incoming γ ray.

The ACD phase of the event reconstruction starts by estimating the energy deposited in each of the tiles and ribbons. Subsequently, these energy depositions are associated with each of the tracks found in the TKR. More specifically, each track is projected onto the ACD, and we calculate whether that track intersects each ACD tile or ribbon with non-zero energy deposition. Furthermore, if the track projection does not actually cross an ACD element with non-zero energy deposition, we calculate the distance of closest approach between the track projection and the nearest such ACD element. Finally, we use the distance calculation and the amount of energy deposited to sort all the possible ACD–TKR associations by how likely they are to represent a charged particle passing through the ACD and into the TKR. This ranking is used when considering whether the event should be rejected in later analysis stages.

For each event, we reduce the output of the TKR, CAL, and ACD reconstruction to a set of a few hundred figure-of-merit quantities whose analyzing power has been studied and optimized with MC simulations.

It is important to note that the best track (and to a lesser extent the second best track) plays a particularly important role in later analyses. Specifically, although we calculate figure-of-merit variables—such as the agreement between the track direction and the axis of the CAL shower, or the distance between the track extrapolation and the nearest ACD energy deposition—for all the tracks in the event, many of the multivariate analyses described in the rest of this section consider only those figure-of-merit variables associated with the two best tracks in the event. Furthermore, the figure-of-merit variables associated with the best track tend to carry significantly more weight in the multivariate analysis—which is of primary importance, as most ghost tracks come from protons and heavy ions and tend to be longer and straighter than tracks from e⁺ and e⁻ from a γ-ray conversion.

Early iterations of the event analysis split the events by energy and then applied different selections in the various energy bands. We found that this approach led to large spectral features at the bin boundaries. Therefore, we chose instead to scale many of our figure-of-merit variables with energy so as to remove most of the energy dependence. This allowed us to have a single set of event selection criteria spanning the entire LAT energy range.

The second step in the event analysis is selecting one energy estimate among those available for the event (see Section 3.2.1). We apply a classification tree (CT) analysis (Breiman et al. 1984; Breiman 1996) to select which of the energy reconstruction methods is most likely to provide the best estimate of the event energy.

Because of the numerous edges and gaps between the logs in the CAL, and because of the huge energy range and large FoV of the LAT, the quality of the energy reconstruction can vary significantly from event to event. Therefore, for each event we also apply a second CT analysis to estimate the probability that the estimated energy is within the nominal core of the energy dispersion. Specifically, we define a scaled energy deviation, as described in more detail in Section 7.1.1, from which most of the energy and angular dependence is factored (i.e., the energy dispersion in the scaled deviation is largely energy and angle independent). We then train CTs that provide probability estimates this event is less than 2σ (P_2σ) or 3σ (P_3σ) away from the most probable value of the energy dispersion. Finally, we define a reconstruction quality estimator P_E by combining P_2σ and P_3σ:

$$\begin{equation} P_{\rm E}= \sqrt{P_{2\sigma } P_{3\sigma }}. \end{equation} \tag{ 4 }$$

Large values of P_E indicate that the event is likely to be in the core of the energy dispersion, and so have an accurate energy estimate.

In the Pass 7 analysis we did not use the energy estimates from the LH algorithm: by construction it is a binned estimator, and energy assignments tended to pile up at the bin boundaries (see Section 7.2 for a more detailed discussion of the effects of removing the LH algorithm from consideration in the energy assignment). The removal of the LH energies causes a somewhat degraded energy resolution for those events where it would have been selected; we compensate for the loss of this energy estimator by requiring a slightly more stringent cut on the energy reconstruction quality when defining the standard event classes (e.g., see Section 3.4.3).

The third step selects the measured direction of the incoming γ ray from the available options. Those options are the directions as estimated from (1) the best track, (2) the best vertex, and (3) and (4) the same two options using the centroid of the energy deposition in the CAL as an additional constraint on the direction of the incident γ ray. Again, we use a CT analysis to combine the information about the event and determine which of the methods is most likely to provide the best direction measurement.

As with the energy analysis, the quality of the direction reconstruction can vary significantly from one event to the next. In this case, we have the additional complication that γ rays can convert at different heights in the TKR, giving us anywhere between 6 and 36 track position measurements. Therefore, we use an additional CT analysis to estimate the probability P_core that the reconstructed direction falls within the nominal 68% containment angle as defined by the simplified analytical model:

$$\begin{equation} C_{68}(E) = \sqrt{ \left[c_0.\left(\frac{E}{100\,{\rm MeV}} \right)^{-\beta }\right]^2+c_1^2} \end{equation} \tag{ 5 }$$

(where the values of the coefficients for front- and back-converting events are listed in Table 4). The reader will notice the similarity with the functional expression used for the PSF prescaling described in Section 6.1.

There is a particular subtlety to the event reconstruction that merits a brief discussion. Specifically, the different event reconstruction algorithms we use provide a set of choices for the best energy and direction of each event. The stages of the event-level analysis described in the previous two subsections select the algorithms that are more likely to provide the best estimates. Therefore, although the event reconstruction was unchanged except for the energy unbiasing, changes in the event-level analysis can result in individual events being assigned slightly different directions and/or energy estimates in Pass 7 with respect to Pass 6. Figure 13 illustrates such differences for the P7SOURCE events above 100 MeV that were also included in the P6_DIFFUSE γ-ray class.

Zoom In Zoom Out Reset image size

The next step of the event analysis starts the process of classifying events as γ rays or particle backgrounds, specifically identifying events for which evidence clearly indicates that a charged particle entered the LAT from inside the FoV. We refer to this stage as the charged particle in the FoV (CPF) analysis. To accomplish this, we first apply a series of filters that classify as background those events for which the best track has hits all the way to the edge of the TKR and (1) points to an active region of the ACD that has significant deposited energy or (2) points directly to less sensitive areas in the ACD, such as the corners or gaps between ACD tiles.

After applying these filters we attempt to account for cases where the best reconstructed track does not represent the incoming particle well. This can happen for a variety of reasons, for example, when the bias of the track-finding algorithm toward the longest, straightest track causes us to incorrectly assign a backsplash particle from the CAL as the best track, or when the incoming particle passes largely in the gaps between the TKR towers and we incorrectly select some secondary particle as the best track. We classify as background the events for which the total energy deposited in the ACD is too large to be accounted for as backsplash from the CAL. It is important to note that this particular requirement loses the benefit of the ACD segmentation, so we apply a conservative version to all events, and a tighter version to events that have a track that passes within an energy-dependent distance of significant energy deposition in the ACD. During extremely bright SFs this total ACD energy requirement can cause us to misclassify many γ rays as CRs because of energy deposition in the ACD from SF X-rays (see Appendix A).

All the events that are classified as charged particles at this stage are removed from further consideration as γ-ray candidates and are passed to separate analyses to identify various species of charged particles (see, e.g., Abdo et al. 2009b; Ackermann et al. 2010, 2012a, 2012b).

Finally, we perform another CT analysis to combine all the available information into an estimate of the probability that the event is from charged particle backgrounds (P_CPF). Although the performance of this selection depends on energy and angle, roughly speaking more than 95% of the background in the telemetered data is removed by means of these cuts, and about 10% of the γ-ray sample is lost.

We note in passing that electrons and positrons cause electromagnetic showers that look extremely similar to γ-ray interactions, and the remaining stages of the event analysis are based on the topology (e.g., the shape, density, and smoothness) of the energy deposition in the TKR (Section 3.3.6) and CAL (Section 3.3.7), and have very little additional discriminating power against such backgrounds.

In the next stage of the event analysis, we use information from the TKR to identify CR backgrounds that were not identified as such by the CPF analysis. These events are not immediately removed from the γ-ray analysis but only flagged to allow for removal from the higher purity event classes. In this part of the analysis, we flag events saturating the energy deposition in the TKR planes with a pattern expected for heavily ionizing particles as well as events with very high energy deposition in the first hit layers of the TKR, which is a signature of a particle coming from the CAL that ranges out in the middle of the TKR.

After these two flags are applied, we divide the events among five branches for evaluating their topologies in the TKR: the first one collects events with a vertex solution in the TKR, the remaining ones separate events with one track and events with many tracks, treating separately events converting in the thin and thick section of the TKR. This is done because each of these different topologies has significantly different ratios of γ-ray signal-to-background contamination, and presents accordingly different difficulty for selecting γ rays. For each branch we apply a different initial selection to remove many CR background events and then pass the remaining events to a CT analysis, which estimates the probability (P_TKR) the event comes from γ-ray interactions as opposed to charged particle backgrounds (and more specifically hadronic charged particle backgrounds).

The variables used in the CT analysis were designed to emphasize differences between the characteristics of γ-ray conversions and hadronic CRs. This includes distinguishing between electromagnetic and hadronic showers by counting extra hits near the track, quantifying the complexity of the event in the TKR and how collimated the overall track structure is, looking how deep into the fiducial volume the event starts, and requiring that the ionization along the track is consistent with e^± pairs.

Next we apply a CAL topology analysis that is similar in design to the TKR topology analysis. The first part of this analysis consists of few general cuts flagging events coming from the bottom and sides of the LAT. As with the TKR topology analysis, we then split the events into five branches depending on the topology; for each branch, we apply a cut and a CT stage. Again, CTs are trained to select γ rays versus hadronic charged particle backgrounds, and provide an estimate of the probability that the event is a γ ray (P_CAL).

The variables used in the CAL topology analysis are not only CAL-derived quantities, but also involve comparisons between what is observed in the CAL and what is seen in the TKR. Among the important discriminants are how well the track solution points toward the CAL energy centroid, how well the direction derived from the CAL (via the moments analysis, see Section 3.2.1) agrees with the track direction, the root mean square (rms) width of the CAL shower, and the ratio of the highest energy in a single crystal to the total energy in all crystals.

After the analyses related to the three LAT subsystems, the last stage of our event analysis applies a final CT analysis utilizing all available information, notably including the output of the CT analyses from previous phases of the event analysis. This second-order CT analysis is particularly important for defining event classes with high γ-ray purity, as discussed in Section 3.4.

At this point, we have a number of specifiers of event reconstruction and classification quality on an event-by-event basis.

• 1.  
  Energy reconstruction quality P_E (see Section 3.3.2).
• 2.  
  Direction reconstruction quality P_core (see Section 3.3.3).
• 3.  
  γ-ray probabilities from the CPF analysis (P_CPF, see Section 3.3.5) and the TKR (P_TKR, see Section 3.3.6) and CAL (P_CAL, see Section 3.3.7) topology analysis.
• 4.  
  Overall γ-ray probability P_all from the final classification step of the event analysis.

These are the basic ingredients for defining the standard γ-ray classes described in Section 3.4.

For the definition of the standard γ-ray classes we confine ourselves to events that have some chance to be useful for most standard science analysis, namely events for which we have enough information to reconstruct an energy and a direction. Therefore, we require that the event has a reconstructed track to allow for a direction estimation. Furthermore, we require that the track points into the CAL and crosses at least 4 radiation lengths of material in the CAL and that at least 5 MeV of energy is deposited in the CAL. This second requirement reduces the fiducial volume of the detector by rejecting off-axis γ rays that pass near the top of the TKR and miss the CAL. Events with no tracks in the TKR or with less energy deposited in the CAL are not considered further as candidate γ rays. The remaining data set (i.e., the events that are passed along as potential candidates for the standard γ-ray classes) is still composed almost entirely of background events.

We note that these cuts remove from consideration two classes of events that might be useful for specific, non-standard, analyses. The first class consists of events that deposit all their energy in the TKR, either because they range out before they reach the CAL or simply because they miss the CAL entirely. In general, the energy resolution is much poorer for these events, though at low energies (<100 MeV) the CAL does not improve the energy resolution significantly. These events have been used effectively in the analysis of both gamma-ray bursts (GRBs; Pelassa et al. 2010) and SFs (Ackermann et al. 2012c). The second class consists of events that do not have reconstructed tracks, but have enough information in the CAL to derive an estimate of the event direction (though without the TKR information the angular resolution for these events is highly degraded). These events can occur because the γ ray entered the LAT at a large incidence angle and missed most of the TKR, or simply because the γ ray passed through the TKR without converting.

For the analysis of brief transient sources (e.g., GRBs) a high purity is not required as the time selection itself limits the amount of background counts in the region of interest (ROI). Accordingly, we define a P7TRANSIENT event class with only a few cuts, with the aim of achieving a residual background rate of a few Hz while maintaining a large efficiency for γ rays. The list of cuts is short.

• 1.  
  The event must pass minimal cuts on the quality estimators P_E and P_core (P_E > 0 and P_core > 0).
• 2.  
  The event energy after all corrections must be greater than 10 MeV.
• 3.  
  The event energy after all corrections must not amount to more than five times the energy deposited in the CAL.
• 4.  
  The CPF analysis must not flag the event as a charged particle (P_CPF > 0.1).
• 5.  
  The event must pass a relatively loose cut on P_all (P_all > 0.2).

It is worth noting that although the cuts are optimized somewhat differently in Pass 6, a similar series of cuts are applied to define the P6_TRANSIENT event class, which is roughly equivalent to the P7TRANSIENT event class.

We have several additional considerations when defining the P7SOURCE event class, which is intended for the analysis of point sources. While the γ rays from a point source are clustered on the sky near the true source position, residual CR background can be modeled as an isotropic component and be accounted for in the analysis of the spectrum or position of the source. Although the exact numbers depend on the details of the PSF and the spectrum of the source, in general we require a background rate of less than ∼1 Hz in the LAT FoV to ensure that we maintain a high enough signal-to-background ratio so that this has little impact on source detection and characterization. Furthermore, in contrast to studies of transient sources, for sources with an integral energy flux above 100 MeV of ∼10⁻¹¹ erg cm⁻² s⁻¹, the precision of spectral studies is often limited by systematic uncertainties below 1 GeV (see Section 8.3), increasing the importance of the precision and accuracy of the event energy and direction reconstruction.

Additionally, we note that for the P7SOURCE event class and both higher purity event classes (P7CLEAN and P7ULTRACLEAN), we developed and optimized the cuts using on-orbit data samples (particularly the clean and dirty samples described in Section 3.6.4) as well as MC simulations. In particular, we performed a comparison of many event quality parameters and reconstruction variables between the clean and dirty samples to identify characteristics of the CR background, which is greatly enhanced in the dirty sample, and devised selection criteria to remove it. Some examples of this procedure are described in more detail in Section 4.3.

An important limitation of this technique is obviously that one needs a very low ratio of residual CR background to γ rays in the clean sample to obtain large contrast between the clean and dirty samples. It can only be applied to data samples which already have a γ-ray-to-residual-background ratio that is comparable to 1. Therefore, when optimizing these selections we only used events that are part of the P7TRANSIENT class (see Section 3.4.1) and additionally have a high value of P_all estimator introduced in Section 3.3.8.

The P7SOURCE event class is defined starting from the P7TRANSIENT event class and includes tighter cuts on many of the same quantities. Specifically, the P7SOURCE event class selection requires the following.

• 1.  
  The event must not have been flagged as background based on topological analysis of the CAL (Section 3.3.6) and TKR (Section 3.3.6) reconstruction, and must pass a tighter cut on P_E; the cut on P_E does reject a noticeable number of events, though by definition these events have less accurate energy estimates on average (P_E > 0.1 or P_E > 0.3 for events which originally used the LH energy estimate).
• 2.  
  The event must pass cuts to reject MIPs based on the event topology in the CAL and the energy deposited in the TKR.
• 3.  
  The event must pass a set of cuts on the agreement between the TKR and CAL direction reconstruction.
• 4.  
  The event must pass a tighter, energy-dependent, cut on P_core:

  $$P_{\rm core}> \max (0.025,0.025+0.175(3.0-\log _{10}(E/1\,{\rm MeV})).$$
• 5.  
  The event must pass a tighter cut on P_all:

  $$\begin{eqnarray*} P_{\rm all}\,&{>}&\, \max (0.7,0.996-1.4\times 10^{-4}\nonumber\\ &&\times(\max (5.4-\log _{10}(E/1\,{\rm MeV}),0)^{5.3})). \end{eqnarray*}$$

As with the TRANSIENT event classes, it is worth noting that although the cuts are optimized somewhat differently in Pass 6, a similar series of cuts are applied to define the P6_DIFFUSE event class, which is roughly equivalent to the P7SOURCE selection.

For the analysis of diffuse γ-ray emission, we need to reduce the background contamination to a level of about ∼0.1 Hz across the LAT FoV, to keep it below the extragalactic γ-ray background at all energies. For comparison, the total Galactic diffuse contribution is ∼1 Hz, depending on where the LAT is pointing, though most of that is localized along the Galactic plane.

The selection of P7CLEAN class events starts from the P7SOURCE event class and includes the following additional cuts.

• 1.  
  The event must pass a series of cuts designed to reject specific backgrounds such as CRs that passed through the mounting holes in the ACD, or the gaps along the corners of the ACD, with minimal costs to the γ-ray efficiency (see Section 4.2).
• 2.  
  The event must pass cuts on the topology of the event in the CAL designed to remove hadronic CRs.

As with the two previous event classes, a similar series of cuts was applied to define the P6_DATACLEAN event class, which is roughly equivalent to the P7CLEAN event class.

For the analysis of extragalactic diffuse γ-ray emission, we need to reduce the background contamination even further below the extragalactic γ-ray background rate to avoid introducing artificial spectral features. As illustrated in Figure 28, the residual contamination for the P7ULTRACLEAN class is ∼40% lower than that of the P7CLEAN class around 100 MeV (the residual levels becoming more similar to each other as the energy increases and becoming the same above 10 GeV).

The selection of P7ULTRACLEAN class events is relatively simple, consisting of a tighter, energy-dependent cut on P_all:

$$\begin{eqnarray*} P_{\rm all}\,&{>}&\, 0.996-0.0394(\max (3.26-\rm{log}_{10}(E/1 \,{\rm MeV}),0)^{1.78}) \nonumber\\ &&\qquad\qquad\qquad\qquad {\rm (Front)}\\ P_{\rm all}\,& {>}&\, 0.996-0.006(\max (4.0-\rm{log}_{10}(E/1 \,{\rm MeV}),0)^{3.0})\nonumber\\ &&\qquad\qquad\qquad\qquad {\rm (Back)}. \end{eqnarray*}$$

The Vela pulsar (PSR J0835−4510) has the largest integral flux >100 MeV of any γ-ray source, and has been well studied in the LAT energy range (Abdo et al. 2009e, 2010c). Furthermore, the pulsed nature of high-energy γ-ray emission gives us an independent control on the background. In fact, off-pulse γ-ray emission is almost entirely absent. These factors combine to make the Vela pulsar an almost ideal calibration source. Unfortunately, the spectrum of Vela cuts off at ∼3 GeV and Vela is nearly undetectable above 30 GeV.

The selection criteria we use to define our Vela calibration samples are listed in Table 6. Specific calibration samples used for particular studies may include additional requirements. For example, the "P7TRANSIENT Vela calibration sample" includes all events in the P7TRANSIENT event class that pass the Vela calibration sample criteria. We used the TEMPO2 package⁷⁶ (Hobbs et al. 2006) and a pulsar timing model⁷⁷ derived from data taken with the Parkes radio telescope (Abdo et al. 2009e; Weltevrede et al. 2010) to assign a phase to each γ ray.

We can achieve excellent statistical background subtraction for any distribution (i.e., spectrum, spatial distribution, any discriminating variable used in the event classification) by subtracting the distribution of off-pulse events (defined as phases in the range [0.7, 1.0]) from the distribution of on-peak events (defined as phases in [0.125, 0.175]∪[0.5125, 0.6125]). Figure 15 shows the distribution of phases of γ rays in the P7TRANSIENT Vela calibration sample, including our standard on-peak and off-pulse regions. (Note that the PSF analysis in Section 6.2.1 uses slightly different definitions of on-peak ([0.12, 0.17]∪[0.52, 0.57]) and off-pulse [0.8, 1.0] regions; the slight difference in definition causes no significant change in the results.)

Zoom In Zoom Out Reset image size

Admittedly, all bright pulsars are potential candidates; nonetheless, for this analysis the consequent increase in event statistics does not warrant the procedural complication required to deal with a stacked sample of pulsars (see also Section 6.2).

At 1 GeV the 95% containment radius for the P7SOURCE event class is ∼14 for front-converting events and ∼25 for back-converting events (see Section 6). Given the density of bright γ-ray sources in the sky, above this energy the instrument PSF becomes narrow enough that we can use the angular distance α between a γ ray and the nearest celestial γ-ray source as a good discriminator for background subtraction, particularly at high Galactic latitudes where there are fewer sources and the interstellar diffuse emission is less pronounced. Unfortunately, no single source is bright enough to provide adequate statistics to serve as a good calibrator. However, by considering γ rays from a sample of bright and/or hard spectrum active galactic nuclei (AGNs) that are isolated from other hard sources, we can create a good calibration sample for energies from 1 GeV to 100 GeV.

Table 7 lists the AGNs that we use here, Figure 16 shows their positions, and Figure 17 shows a comparison of their spectral index (Γ) and integral γ-ray flux between 1 GeV and 100 GeV (F₃₅) to other AGNs in the second LAT source catalog (2FGL; Nolan et al. 2012). Note that different but overlapping sets of AGNs are used for A_eff- and PSF-related studies. More information about the source properties can be found in Abdo et al. (2010b) and Nolan et al. (2012).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The selection criteria we use to define AGN calibration samples are listed in Table 6. As with the Vela calibration samples, specific calibration samples used for particular studies may include additional requirements (e.g., the P7TRANSIENT AGN calibration sample).

To use this calibration sample to perform background subtraction, we define signal and background regions in terms of the angular distance α between the γ ray and the closest AGN. Specifically, we use α < 05 for the signal region and the annulus 387288 < α < 4° for the background region. These ranges are chosen so that the background region contains four times the solid angle of the signal region. We adopted the separation between the signal and background regions in order to minimize signal γ rays from the tails of the PSF leaking into the background region. Figure 18 shows the squares of the angular separations between γ rays and the nearest bright AGN for all γ rays in the P7TRANSIENT event class in 6° regions around the 25 AGNs listed as the A_eff calibration sample in Table 7, including our definitions of source and background regions.

Zoom In Zoom Out Reset image size

The Earth's atmosphere is a very bright γ-ray source. Furthermore, at energies above a few GeV the γ-ray flux seen by the LAT is dominated by γ rays from the interactions of primary CR protons with the upper atmosphere. This consideration, together with the narrowness of the PSF at energies >10 GeV, causes the Earth limb to appear as a very bright and sharp feature, which provides an excellent calibration source. Furthermore, we have selected data from 200 orbits during which the LAT was pointed near the Earth limb as the basis of the Earth limb calibration sample. The selection criteria we use to define the Earth limb calibration samples are listed in Table 6. When using the Earth limb as a calibration source we generally limit the energy range to energies >10 GeV, primarily because at lower energies orbital variations in the geomagnetic field significantly affect the γ-ray fluxes (however see Section 7.5 for a counterexample).

For this calibration source, we define our signal region as 1111002 < θ_z < 1129545 and background regions just above and below the limb: 1086629 < θ_z < 1095725 and 1145193 < θ_z < 1154675. Note that these ranges are defined to give the same solid angle (0.06π sr) in the signal and background regions.

Figure 19 shows the zenith angle distribution for all γ rays in the P7TRANSIENT event class for the Earth limb calibration sample.

Zoom In Zoom Out Reset image size

At energies above ∼30 GeV, no single source provides enough γ rays for a good comparison between flight data and MC. However, the combination of bright Galactic sources and Galactic diffuse backgrounds means that there is a very large excess of γ rays coming from the Galactic plane relative to high Galactic latitudes.

The intensity of the γ-ray emission at low latitudes in the inner Galaxy is more than an order of magnitude greater than at high latitudes in the outer Galaxy. In contrast, the intensity of the CR background is approximately isotropic for observation periods longer than the 53.4 day orbital precession period.

Unfortunately, since the Galaxy extends over much of the sky, and since the data set consists of several thousand orbits it is not practical to disentangle the variations of exposure from the spatial variations in Galactic diffuse emission without relying on detailed modeling of the Galactic diffuse emission. Accordingly, we use the Galactic ridge primarily when we are developing our event selections, rather than for precise calibration of the LAT performance.

Specifically, we developed the event classes that require a high γ-ray purity, i.e., the event classes used in the analysis of celestial point sources and diffuse emission (see Section 3.4), in part by tuning our selection criteria to maximize the contrast between regions in the bright Galactic plane and at high Galactic latitudes. This helped to mitigate the risk that insufficient statistics in the MC training sample or limited accuracy of the MC description of the geometry of the detector and the particle interactions in the LAT limited the discriminating power and accuracy of the event classification analysis.

The selection criteria we use to define the Galactic ridge calibration samples are listed in Table 6. In particular, we use the region (|b| < 15, −40° < l < 50°, which was selected to maximize the total γ-ray flux) to define a clean data sample and the region (|b| > 50°, 90° < l < 270°) to define a dirty data sample. The ratio of γ rays to CR background in the clean region is more than an order of magnitude higher than the same ratio in the dirty region. Figure 20 shows the count maps for the P7SOURCE samples for both regions.

Zoom In Zoom Out Reset image size

Furthermore, to give a sense of the statistics of these samples at high energies, Figure 21 shows the Galactic latitude distribution for all γ rays in the P7TRANSIENT event class with energies above 17783 MeV (i.e., log₁₀(E/1 MeV) > 4.25).

Zoom In Zoom Out Reset image size

As we will see in the next sections, the IRFs depend heavily on θ (and, to a lesser extent, on ϕ). Therefore, any detailed comparison between flight data and MC simulations must account for the distribution of observing profile, particularly t_obs(θ). How best to account for the observing profile depends on the particulars of the calibration samples.

For any point source, the observing profile is determined by the position of the source, the rocking angle of the LAT and the amount of time spent in survey mode relative to pointed observations. Figure 22 shows the observing profile for Vela for the first two years of the mission. Rather than produce a dedicated large statistics MC sample for Vela, we re-use our allGamma MC sample, re-weighting the events in that sample so as to match the Vela observing profile.

Zoom In Zoom Out Reset image size

Similarly, we re-weight the allGamma MC to match the summed observing profiles of all of the AGNs of our sample, which is shown in Figure 23. Unfortunately, since AGNs are intrinsically variable, and since the AGNs in this sample span a range of fluxes, this re-weighting technique will not work as well with this sample. On the other hand, by taking a large set of AGNs, we reduce the bias due to the variability of any one particular source. In broad terms, our re-weighted MC sample reproduced the θ distribution of the AGN sample to better than 2% (see Section 5.5). Finally, we note that since the PSF is narrower above 1 GeV, and the ROI around each AGN is only 6° we do not apply the ROI-based θ_z cut when building the AGN calibration samples.

Zoom In Zoom Out Reset image size

Since the Earth limb is a spatially extended source, we cannot apply the re-weighting technique we used for the Vela and AGN samples to account for the observing profile. On the other hand, since the data set consists of only 200 orbits, and the Earth limb emission is well understood above 10 GeV we can produce an MC simulation of the Earth limb emission for those orbits and compare it with the flight data (see Section 2.5.2).

Finally, Figure 24 shows the statistics available for each of the samples. This shows that the calibration sources span most of the LAT energy range, certainly from 30 MeV up to at least 100 GeV.

Zoom In Zoom Out Reset image size

As explained in more detail in Section 2.1.4, after the start of LAT operations, it became apparent that ghost signals led to a significant decrease in effective area with respect to the pre-launch estimates, for which this effect was not considered. The overlay procedure used to account for this effect, first introduced in the P6_V3 set of IRFs, is described in detail in Section 2.5.1 and its impact on the effective area is shown in Figure 34.

Zoom In Zoom Out Reset image size

The effect of ghost signals is corrected on average as described in the previous section. A smaller correction is necessary to account for the detailed dependence of A_eff on the CR rates. To account for this we need an estimator for the rate of CRs entering the LAT; the obvious one is the trigger rate, but technical issues make this choice impractical. A variable that can be easily obtained from the pointing history files and which is linearly correlated with trigger rate is the live time fraction F_l, the fraction of the total observing time in which the LAT is triggerable and not busy reading out a previous event. The average value of F_l is ∼90% and varies between 82% and 92% over the Fermi orbit.

We bin events from a sample of periodic triggers according to the corresponding live time fraction and for each bin we produce a dedicated allGamma simulation (Section 2.5.3); for each of the resulting overlay data sets we derive A_eff. We have found the dependence on the incidence angle to be small and so we choose to neglect it when studying this effect. Furthermore, we have found that at a given energy A_eff varies linearly with the live time fraction. We perform a linear interpolation in each energy bin in accord with

$$\begin{equation} A_{\rm eff}(E, F_l) = A_{\rm eff}(E) \cdot (c_0(E) F_l + c_1(E)) \end{equation} \tag{ 14 }$$

separately for front and back events.

In Figure 35, we plot c₀ and c₁ as a function of energy for front-converting events in the P7SOURCE event class. As shown by the solid lines, we use a simple piecewise linear fit to describe the energy dependence and the fit parameters are stored in the A_eff tables. As we mentioned, the effective area derived from allGamma simulations that have overlaid periodic trigger events is effectively corrected for the average effect of this live time dependence, so we treat this additional correction (which can be positive or negative) as a modulation of the tabulated A_eff. Correction parameters are read in from the A_eff files of each set of IRFs, and used to correct the calculated exposures.

Zoom In Zoom Out Reset image size

The resulting corrections to the average A_eff can reach −30% at 100 MeV, decreasing at higher energies to <5% above ∼10 GeV. The uncertainties in the corrections are much smaller; studies using flight data confirm the MC-based predictions to better than 2%. Note that over a 53.4 day orbital precession period this effect will tend to the overall average correction described in Section 5.2.1, with less than 1% variation across the sky.

Since the correction to the effective area is based on live time fraction, which is very strongly correlated with the CR intensity, it avoids any direct biases from long term changes in the CR intensity associated with the influence of solar activity on the geomagnetic field. However, the correction does not address the possibility that the CR population changes during the solar cycle in such a way as to change the effective area dependence on the live time fraction. Given the small change in the daily averaged LAT trigger accept rate observed in the mission to date (<5%), we neglect this effect.

The tabulated values of A_eff are averaged over the azimuthal angle of incidence and shown in Figure 36. Much of the azimuthal dependence of the effective area is geometrical, due to the square shape of the LAT and the alignment of the gaps along the x- and y-axes. The rms variation of the effective area as a function of ϕ is typically of the order of 5% and exceeds 10% only at low energies (<100 MeV) or far off-axis (θ > 60°) where the effective area is small, and at very high energies (>100 GeV) where the event rate is small.

Zoom In Zoom Out Reset image size

In order to parameterize the azimuthal dependence of A_eff, we fold the azimuthal angle into the [0°, 45°] interval and remap it into [0, 1] using the following transformation:

$$\begin{equation} \xi = \frac{4}{\pi }\left| \left(\phi {\rm mod}\,\frac{\pi }{2} \right)- \frac{\pi }{4}\right| \end{equation} \tag{ 15 }$$

(the transformation maps 0° to 1, 45° to 0, 90° to 1, and so on). The allGamma events are binned in energy and θ and in each bin a histogram of ξ is fitted (see Figure 37 for an example). Front- and back-converting events are treated separately. The fitting function is

$$\begin{equation} f(\xi) = 1 + q_0 \xi ^{q_1}. \end{equation} \tag{ 16 }$$

The absolute scale is not important: we normalize the correction to result in an average multiplicative factor of 1, so that the average A_eff is tabulated. The fitted parameters q₀ and q₁ are stored in the A_eff tables.

Zoom In Zoom Out Reset image size

We note in passing that the default for high-level analyses is to disregard the azimuthal variations when calculating exposures because they average out when intervals of months or longer are considered.⁸² Although the combined θ and ϕ dependence of the observing profile averages out only on year-long timescales, the eightfold symmetry of the LAT combined with the rotation of the x-axis to track the Sun results in effective averaging over ϕ on short timescales. In fact, we have found that ignoring the ϕ dependence of the effective area results in only a small variation of the exposure on 12 hr timescales (<3% rms at all energies).

To measure the efficiency of a cut on flight data, we can perform background subtraction before and after the cut and compare the excess in the signal region in the two cases. To the extent that the background subtraction is correct, the cut efficiency is simply the ratio of the number of background-subtracted events before and after the cut we are testing

$$\begin{equation} \eta = \frac{n_{s,1} - r n_{b,1}}{n_{s,0} - r n_{b,0}}, \end{equation} \tag{ 18 }$$

where the subscripts s and b indicate the signal and background regions, the subscripts 0 and 1 indicate the samples before and after the cut under test, and r is the ratio of the size of the background region to the size of the signal region. Note that the sample after the cut is a subset of the sample before the cut. Therefore, for a reasonable sample size (n_s > 10), the statistical uncertainty of the efficiency is

$$\begin{equation} \hbox{\fontsize{7.9}{9.9}\selectfont{$\delta \eta = \left[ \dsty\frac{[(n_{s,0} - n_{s,1}) + r^{2}(n_{b,0} - n_{b,1})]\eta ^{2} + (n_{s,1} + r^2 n_{b,1})(1-\eta)^{2}}{(n_{s,0} - r n_{b,0})^{2}} \right]^{1/2}.$}} \end{equation} \tag{ 19 }$$

An example application of this technique using data from the Vela calibration data set is shown in Figure 38.

Zoom In Zoom Out Reset image size

In many cases our cuts have significant overlap—two cuts may reject many of the same events. In such cases, measuring the efficiency of the second cut after the first cut has been applied would give a quite different result than measuring the efficiency of the second cut without the first cut. Therefore, whenever possible, we use the background subtraction technique described above to evaluate the efficiency of each step of our event selection both independently of the other steps as well as after all the other cuts have been applied.

Although the track-finding and fiducial cuts are applied midway through our event selection process, we choose to discuss them first for two reasons: (1) we require an event direction to be able to perform a background subtraction on flight data and (2) we express A_eff as a function of E and θ. Therefore, for performance studies we need to apply some minimum event quality and fiducial cuts before considering events for analysis. The standard cuts require at least one track found, with at least 5 MeV of energy in the CAL and that the track extrapolates through at least four radiation lengths in the CAL (see Section 3.4.1).

Most events that fail these cuts have either poorly reconstructed directions, poorly reconstructed energies, or both—which makes it difficult to study the performance as a function of energy and direction. Therefore, we choose to study the efficiency of the track-finding and fiducial cuts by selecting events that almost fail these cuts in the P7SOURCE calibration samples. Specifically, we study the fraction of events that are very close to the cut thresholds and verify that the flight data agree with our MC simulations, see Figure 39.

Zoom In Zoom Out Reset image size

The LAT hardware trigger, trigger configuration, and on-board filter are described in Section 3.1. For our purpose here we first consider the fractions of γ-ray events that have one of the five physics trigger conditions (TKR, ROI, CAL_LO, CAL_HI, CNO) asserted. Furthermore, only two of them (TKR and CAL_HI) effectively serve to initiate a trigger request. Of the others, ROI does not exist without TKR, CAL_LO alone does not open the trigger window and from the point of view of selecting γ-ray events CNO is primarily a veto rather than a trigger. Since the CAL_HI requires at least 1 GeV in a CAL channel, and our fiducial cuts and quality cuts require at least one track, we are effectively using the TKR condition as the primary trigger for γ rays up to such energies that the CAL_HI is very likely to be asserted.

Using the diagnostic filter events described in Section 3.1, we can measure the fractions of all trigger requests that have individual primitives asserted. However, because of the high particle background rates this does not really probe the trigger stage of the γ-ray selection process. On the other hand, for each of the five relevant trigger primitives we measure this efficiency as a function of energy for the allGamma sample as well as the P7TRANSIENT selection. For the P7TRANSIENT selection, we also measure the fractions of events with each trigger primitive asserted and compare these to the MC predictions. These are shown in Figure 40. The only notable discrepancy is that the MC overpredicts the fraction of events having ROI asserted at very high energies; this is likely related to imperfect simulation of the backsplash from the CAL, and since the CAL_LO and CAL_HI are typically asserted for these events, it does not affect the trigger readout decision (see Table 3).

• 1.  
  The CAL_LO trigger condition requires 100 MeV of energy to be deposited in any CAL channel. This condition starts to be asserted for γ rays at ∼300 MeV and reaches full efficiency for events contained within the CAL at ∼1 GeV.
• 2.  
  The CAL_HI trigger condition requires that 1 GeV of energy be deposited in any CAL channel. This condition starts to be asserted for γ rays at ∼3 GeV and reaches full efficiency at ∼15 GeV.
• 3.  
  The CNO condition requires any ACD tile to have a very large signal, consistent with the passage of a heavy ion. The CNO condition actually serves more to veto than to select an event for readout, but only becomes active at very high energies (>100 GeV), where CAL_LO and CAL_HI are already active.
• 4.  
  The ROI condition requires any ACD tile in a predefined region of interest associated with a TKR tower to have a signal above 0.45 MIP. The ROI condition actually serves more to veto than to select an event for readout and becomes active at a few GeV.
• 5.  
  Since the TKR condition serves as the primary trigger, it is very difficult to measure the efficiency of the condition. However, we can estimate how well the MC simulates this efficiency by studying how well it models cases where the events almost fail the trigger request conditions. This is possible by calculating the fraction of events with exactly one combination of three hit layers in a row in the TKR. Most high-energy γ rays have many such combinations and would have triggered the LAT even if one hit had been lost. However, for back-converting low-energy events (<100 MeV) the fraction of single combination events becomes significant.

Zoom In Zoom Out Reset image size

Since the fiducial cuts for our standard event selection require that a track is found, and that the track extrapolates into the CAL, the interplay between the trigger and the standard event selections are actually quite simple.

• 1.  
  At very low energies, the LAT trigger starts to become efficient at ∼10 MeV and follows the efficiency of the TKR condition, which becomes fully efficient by ∼100 MeV.
• 2.  
  At around 1 GeV the CAL_LO condition becomes active. By design, this is considerably lower than the 10 GeV where the ROI starts to be asserted because of backsplash.
• 3.  
  Above ∼10 GeV the CAL_HI condition becomes active and we no longer rely on the TKR condition as the primary driver of the trigger.

Furthermore, events that are rejected because ROI is asserted are extremely unlikely to pass standard event class selections. Taken together, this means that the only part of the LAT energy band where the trigger has strong influence on A_eff is below ∼100 MeV.

Although the gamma filter has many different cuts (see Section 3.1.2), most events that are rejected by the gamma filter would be either rejected by the fiducial cuts (see Section 3.4.1) or by the P7TRANSIENT event selection. Accordingly, we choose to study the efficiency of the gamma filter as a whole, and only for those events which pass all the other selection criteria for the P7TRANSIENT class event sample.

Since we downlink a small fraction of the events that fail the gamma filter (see Section 3.1), we can check that this is indeed the case for this diagnostic sample. However, the large prescale factor (250) for the diagnostic filter and high level of background rejection in the P7TRANSIENT selection severely limit our statistics for this study, and we can do little more than confirm that the gamma filter is highly efficient for events in the P7TRANSIENT sample (see Figure 41).

Zoom In Zoom Out Reset image size

Measuring the efficiency of the selection for the P7TRANSIENT event class is the most technically challenging part of the A_eff validation for two reasons: (1) at the output of the gamma filter the background rates are still great enough to overwhelm almost all traces of γ-ray signals and (2) the event analysis directs events that are tagged as likely due to CRs away from the remainder of the γ-ray selection criteria, which means that many of the CT analyses are never applied, and cannot be used to construct a cleaner sample on which we can measure the efficiency of any of these cuts independently of the other cuts.

Figure 42(a) shows the efficiency of each part of the P7TRANSIENT selection on events from the allGamma sample that have passed the gamma filter. As stated above, the high levels of CR background make it infeasible to use flight data to obtain stringent constraints on the efficiency of these selection criteria. Therefore, similar to what we did for the track-finding and fiducial cuts, we also study the events that almost fail these cuts in the P7TRANSIENT calibration samples: these comparisons are shown in Figure 42(b).

Zoom In Zoom Out Reset image size

Because of the difficulty in validating the efficiency of the P7TRANSIENT event selection criteria, we have chosen to use the consistency checks described in Section 5.5 to estimate the systematic uncertainty in our A_eff representation.

Although the P7TRANSIENT event class is dominated by residual background across the entire LAT energy range, the background levels are at least reduced to the point where the γ-ray signals in the calibration samples described in Section 3.6 are clearly detectable. This makes the validation of the effective area from this point on much easier. We can compare the efficiency of each cut, as measured on flight data, to the MC prediction with the method described in Section 5.3.1. These comparisons are shown in Figure 43.

Zoom In Zoom Out Reset image size

From the consistency checks described in the previous section, we arrive at a rough overall estimate of the uncertainty of A_eff, which is shown in Figure 48. Note that this estimate is assigned simply to account for the largest observed inconsistency, namely the mismatch between the front- and back-converting events. Roughly speaking, for the P7SOURCE_V6 and P7CLEAN_V6 event classes these uncertainties may be quoted as 10% at 100 MeV, decreasing to 5% at 560 MeV, and increasing to 10% at 10 GeV and above. It is important to note that these uncertainties are statements about overall uncertainty of A_eff at various energies, and do not include any statement about what types of deviations we might expect within the stated uncertainty bands, nor about the point-to-point correlations in any systematic biases of A_eff. Those questions are addressed in the next sections.

Zoom In Zoom Out Reset image size

Since our selection criteria are generally scaled with energy (see Section 3.3.1), we expect any biases of A_eff to be highly correlated from one energy band to the next. This point is very important when estimating the size of potential instrumental spectral artifacts.

We have studied the point-to-point correlation of the effective area through the consistency checks described in Section 5.5, where it is evident that the deviations from unity are not independent for neighboring energy bins. In order to quantify this correlation we first scale the values r_i of the data-to-MC ratio R in each of the N energy bins (indexed by i), turning them into normalized deviations

$$\begin{equation} d_i = \frac{(r_i - m)}{\sqrt{\frac{1}{(N - 1)}\sum _{i = 1}^{N}(r_i - m)^2}}, \end{equation} \tag{ 20 }$$

and m is the arithmetic mean of the r_i) and then we construct the metrics

$$\begin{equation} \tau _n = \frac{\sum _{i = 1 + n}^{N}(d_i - d_{i - n})^2}{(N - n)}. \end{equation} \tag{ 21 }$$

The quantity in Equation (21) is related to a reduced χ²_{N − n}(2x); it is small for highly positively correlated deviations (in which case the differences between neighboring bins are generally small), while the expectation value is 2, for all values of n, for normal uncorrelated errors.

Figure 49 shows this metric, for different values of n, calculated on the front versus back consistency check shown in Figure 47. (Only the Vela calibration data set is considered here.) The extremely small value of τ₁ indicates that neighboring logarithmic energy bins are highly positively correlated (as can be naively inferred from the plot in Figure 47), while on the scale of half a decade in energy (n = 4) there is little evidence of a correlation. This implies that the systematic uncertainties on the effective area are not likely to introduce significant spectral features over scales much smaller than half a decade in energy (which is much larger than the LAT energy resolution). The results for all the consistency checks are summarized in Table 8.

Zoom In Zoom Out Reset image size

For many analyses, especially those that are limited by statistics, it is enough to consider the overall uncertainty and allow for worst case deviations within the stated uncertainty bands. However, doing so will result in very conservative systematic error estimates. We will discuss this in more detail in Section 5.7 when we describe techniques to propagate the estimates of the uncertainty of A_eff to uncertainties on quantities such as fluxes and spectral indices. Furthermore, we will come back to the issues of point-to-point correlations and the potential induced spectral features when we discuss the uncertainties associated with the energy reconstruction in Section 7.

As a source moves across the LAT FoV and A_eff changes with the viewing angle, any errors in the A_eff parameterization as a function of θ potentially could induce artificial variability. We have searched for such induced variability with Vela. We split the data set into 12 hr periods (indexed by i) and compared the number of γ rays observed (n_i) during each period with the number of γ rays we predict ($\tilde{n}_{i}$) based on the fraction of the total exposure for Vela that we integrated during that 12 hr period.

On average, our Vela calibration sample contains 230 (176) on-peak (off-pulse) P7SOURCE class events in the 100 MeV to 10 GeV energy band every 12 hr. Since the off-pulse region is twice the size of the on-peak region, background subtraction yields an average on-peak excess of n_i = 142γ rays with an average statistical uncertainty of σ_i = 16γ rays in each time interval.

The exposure calculation requires several points of input.

• 1.  
  The spacecraft pointing and live time history, which are binned in 30 s intervals.
• 2.  
  The P7SOURCE_V6 A_eff parameterization, which we use when deriving A_eff(E, t), the effective area for γ rays for Vela for each 30 s interval.
• 3.  
  The P7SOURCE_V6MC PSF (see Section 6.1) parameterization, which allows us to calculate the energy-dependent containment C(E, t, 15°) within the 15° ROI for each 30 s interval.
• 4.  
  A parameterization of the spectrum of Vela F(E) so that we may correctly integrate A_eff(E, t) and C(E, t, 15°) over the energy range.

In order to minimize dependence on the modeled flux, we calculate the exposure independently for 24 energy bins (which we index by j). The exposure in a single time and energy bin is

$$\begin{equation} \mathcal {E}_{i,j} = \int ^{E_i} \int ^{t_j} A_{\rm eff}(E,t) C(E,t,15^{\circ }) \frac{F(E)}{\int ^{E_i} F(E) dE} dt dE. \end{equation} \tag{ 22 }$$

We can then express the expected number of γ rays in each time and energy bin ($\tilde{n}_{i,j}$) as a fraction of the total number of γ rays in that energy bin (n_j),

$$\begin{equation} \tilde{n}_{i,j} = \frac{\mathcal {E}_{i,j}}{\Sigma _{i} \mathcal {E}_{i,j}} n_{j}. \end{equation} \tag{ 23 }$$

Then we sum n_{i, j} and $\tilde{n}_{i,j}$ across energy bins to find n_i and $\tilde{n}_{i}$.

We have performed this analysis, dividing the first 700 days of the Vela data sample into 1400 12 hr time intervals and using the phase-averaged flux model

$$\begin{equation} F(E) \propto E^{-\Gamma }e^{-\frac{E}{E_0}} \end{equation} \tag{ 24 }$$

with Γ = 1.38 and E₀ = 1.360 GeV (Abdo et al. 2009e). Figure 50 shows the fractional variation,

$$\begin{equation} \frac{n_{i} - \tilde{n}_{i}}{\tilde{n}_{i}}, \end{equation} \tag{ 25 }$$

the normalized residuals,

$$\begin{equation} \frac{n_{i} - \tilde{n}_{i}}{\sigma _{i}}, \end{equation} \tag{ 26 }$$

and the distribution of the normalized residuals for each of the time intervals. The normalized residuals are very nearly normally distributed. Furthermore, the Fourier transformation of the time series (Figure 51) shows only a small peak corresponding to the orbital precession period and is otherwise consistent with Poisson noise. Note that unlike more complicated analyses that involve fitting the flux of a point source, this analysis tests only the accuracy of the A_eff representation (and to a much lesser extent, the 15° containment of the PSF). We attribute the peak in the Fourier spectrum to small, incidence angle-dependent errors in the effective area. As the orbit precesses, the range of incidence angles sampled, and hence the bias in the calculated exposure, varies slightly.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although we performed this analysis with 12 hr time intervals, as noted in Section 2.3, the LAT boresight follows very similar paths across the sky during successive two-orbit periods. Therefore, the level of instrument-induced variability observed with 12 hr time intervals is likely to be indicative of the systematic uncertainties for variability analyses down to 3 hr timescales.

Although the estimate used in Section 5.6.1 that the systematic uncertainty of A_eff is less than the disagreement between the extreme cases is quite conservative for long timescale observations, it is somewhat less conservative for shorter observations. For example, in observations less than the Fermi orbital precession period of ∼53.4 days, a particular region of the sky might be preferentially observed at incidence angles where the bias of A_eff is particularly large, or during parts of the orbit in which Fermi is exposed to particularly high CR background rates and the correction described in Section 5.2.2 leaves some residual bias in the calculated exposure. Finally, we have observed that ignoring the ϕ dependence of the effective area (Section 5.2.3) can induce artificial quarterly periodicity in the fluxes from directions near extremely bright sources, in particular the Vela pulsar.

A somewhat brute force approach to this problem is to generate IRFs that represent worst case scenarios for measuring specific quantities like fluxes or spectral parameters and use these bracketing IRFs to repeat the analysis and extract the variation in the measured quantities. Of course, the nature of the variations between the IRFs depends on the quantity in question.

We address this in a generic way by scaling A_eff by the product of the relative systematic uncertainty

$$\begin{equation} \epsilon (E) = \frac{\delta A_{\rm eff}(E)}{A_{\rm eff}(E)} \end{equation} \tag{ 27 }$$

and arbitrary bracketing functions B(E) taking values in the [ − 1, 1] interval. Specifically, we define modified A_eff as

$$\begin{equation} A_{\rm eff}^{\prime }(E, \theta) = A_{\rm eff}(E, \theta) \cdot \left(1 + \epsilon (E)B(E) \right). \end{equation} \tag{ 28 }$$

The simplest bracketing functions, B(E) = ±1 clearly minimize and maximize A_eff within the uncertainty band. On the other hand, to maximize the effect on the spectral index in a power-law fit, we choose a functional form that changes sign at the pivot or decorrelation energy E₀ (i.e., the energy at which the fitted differential flux and spectral index are uncorrelated):

$$\begin{equation} B(E) = \pm \tanh \left(\frac{1}{k}\log (E/E_{0}) \right). \end{equation} \tag{ 29 }$$

The parameter k controls the slope of the transition near E₀; in practice we use k = 0.13, which corresponds to smoothing over twice the LAT energy resolution of ΔE/E ∼ 0.15. The bracketing IRFs used for effective area studies are listed in Table 9. Figure 52 shows the bracketing functions for c_index_soft and c_index_hard and their effects on the on-axis A_eff.

Zoom In Zoom Out Reset image size

Table 9. Bracketing A_eff and Corresponding Energy-dependent Scaling Functions Used to Create Them

Name	B(E)
c_flux_lo	+1
c_flux_hi	−1
c_index_soft	$+\tanh (\frac{1}{k}\log (E/E_{0}))$
c_index_hard	$-\tanh (\frac{1}{k}\log (E/E_{0}))$

Download table as:  ASCIITypeset image

We have studied the following two sources from the AGN sample to obtain estimates of the effects of instrumental uncertainties on the measured fluxes and spectral parameters.

• 1.  
  PG 1553+113: associated with 2FGL J1555.7+1111, which has a spectral index of 1.66 ± 0.02, making it one of the hardest bright AGNs.
• 2.  
  B2 1520+31: associated with 2FGL J1522.1+3144, which has a spectral index of 2.37 ± 0.02 (when fit with a power law), making it one of the softest bright AGNs.

In each case we used ScienceTools (version v9r25p2) to perform a series of binned maximum likelihood fits to a 20° × 20° region centered at the source position over the energy range 100 MeV–100 GeV. For each individual fit we followed the same procedure.

• 1.  
  Used the γ-ray and time interval selection criteria as for the Vela calibration sample (see Section 3.6).
• 2.  
  Included all 2FGL sources within 20° in our likelihood model, with the same spectral parameterizations as were used in the 2FGL catalog.
• 3.  
  Included models of the Galactic diffuse emission (ring_2year_P76_v0.fits) and the isotropic diffuse emission (isotrop_2year_P76_source_v1.txt) rescaled by the inverse of the function used to rescale the effective area (so as to ensure the same distribution of expected counts from these diffuse sources).
• 4.  
  Freed the spectral parameters for all point sources within 8° from the center of the region, as well as the overall normalizations of both diffuse components.

Following what was done for the 2FGL catalog, we used a simple power law

$$\begin{equation} \frac{dN}{dE} = N_{0} \left(\frac{E}{E_{0}}\right)^{-\Gamma } \end{equation} \tag{ 30 }$$

to model the differential flux from PG 1553+113 and a "log-parabola"

$$\begin{equation} \frac{dN}{dE} = N_{0}\left(\frac{E}{E_{0}}\right)^{-(\alpha + \beta \ln (E/E_{0}))} \end{equation} \tag{ 31 }$$

to model that of B2 1520+31. (All of this will also be relevant for the tests with bracketing PSFs described in Section 6.5.1 and with energy dispersion included in the likelihood fit described in Section 7.4.)

Tables 10 and 11 show the fit results for PG 1553+113 and B2 1520+31, respectively, using these A_eff bracketing functions, as well as the integral counts (F₂₅) and energy (S₂₅) fluxes between 100 MeV and 100 GeV.

The ranges of the fit values indicate propagated uncertainties from the uncertainty in A_eff (Table 12). It is important to note that the systematic error estimates resulting from this technique represent conservative estimates within the instrumental uncertainties, rather than random variations. Furthermore, many of the bracketing IRFs are mutually exclusive, so when considering relative variations between γ-ray sources it is more appropriate to compare how the relative values change for each set of bracketing IRFs.

Alternatively, given a family of plausible A_eff curves, we can use a weighted bootstrap approach (see Efron & Tibshirani 1993, for more details) for propagating the systematic uncertainties on A_eff. The weighed bootstrap approach is closely related to the bracketing IRFs method described in the previous section and to the methods discussed in Lee et al. (2011) in the context of the analysis of Chandra X-ray data.

The basic idea is that, for each trial, the event data are bootstrap resampled using a weighting based on an effective area scaling function that is drawn from a family of plausible curves. The simplest of such families of A_eff curves (Figure 53) can be constructed starting from Equation (28) and multiplying the scaling function (E)B(E) by a normally distributed random number ξ with zero mean—which effectively becomes the parameter controlling the family itself:

$$\begin{equation} A_{\rm eff}^{\prime }(E, \theta , \xi) = A_{\rm eff}(E, \theta)\cdot \left(1 + \xi \epsilon (E)B(E) \right). \end{equation} \tag{ 32 }$$

Zoom In Zoom Out Reset image size

We have found the results of the weighted bootstrap with the family of scaling functions defined in Equation (32) to be in good agreement with those of the bracketing IRFs approach described in the previous section. This should not be surprising as the bracketing IRFs use the ±1σ excursions of the bracketing function and the weighted bootstrap draws from a Gaussian distribution of function scalings that are also based on that same bracketing function. The real benefits of the weighted bootstrap arise when one has families of plausible effective area functions that have more complicated dependencies (e.g., on incidence angle as well as energy) such that exposures would need to be recalculated to apply the bracketing method.

In our parameterized description of the PSF most of the energy dependence is factored into a scaling term:

$$\begin{equation} S_{P}(E) = \sqrt{\left[c_0\cdot \left(\frac{E}{100\,{\rm MeV}} \right)^{-\beta }\right]^2+c_1^2}. \end{equation} \tag{ 34 }$$

Despite a careful investigation, we did not find a simple satisfactory description of the θ dependence to be incorporated in the scaling function.

When building our MC-based PSF, we use a set of scaling function parameters based on pre-launch simulations and confirmed with analysis from beam tests with the Calibration Unit (see Section 7.3.1 for more details). The values of these parameters are shown in Table 13. Note that S_P(E) has same functional form as for C₆₈ in Section 3.3.3, however we have updated the parameters slightly based on the scaling observed in our allGamma sample.

We define the scaled angular deviation x as

$$\begin{equation} x = \frac{\delta v}{S_{P}(E)}. \end{equation} \tag{ 35 }$$

An example of scaled deviation is shown in Figure 55. The effect of the scaling is to make the profile almost independent of energy, in that the maximum is always close to x = 1 for all energy bins while the PSF 68% containment varies by almost two orders of magnitude from 100 MeV to 100 GeV.

Zoom In Zoom Out Reset image size

Before the fit is performed, each scaled deviation histogram is converted into a probability density with respect to solid angle. The result is illustrated in Figure 55.

Note that, although the scaling removes most of the energy dependence, the simulation indicates significant variation of the PSF with θ. At larger incidence angle the tracks must cross more material in each TKR plane. At energies below ∼1 GeV, this degrades the PSF owing to the increased multiple scattering, while at higher energies (above ∼1 GeV), the additional complication of hard scattering processes in the TKR and additional hits in the TKR from the nascent electromagnetic shower complicate the track-finding and degrade the PSF. Figure 56 shows how the scaled containment radii evolve with energy and incidence angle. To test the ϕ dependence of the PSF we have also measured the containment radii independently for events with ξ > 0.33 and ξ < 0.33, where ξ is the folded azimuthal angle defined by Equation (15). The 68% and 95% containment radii for the two ξ ranges differ by <5% for all energies and angles except at high energies (>10 GeV) and low incidence angles (cos θ > 0.7) for back-converting events, and at even higher energies (>100 GeV) and large incidence angles (cos θ < 0.7) for front-converting events. Even then the maximum difference in the 68% containment radius for events is only 10% at 10 GeV and 25% at 100 GeV. In summary, the variations of the containment radii with ϕ are many times smaller than the corresponding variations with θ for all but the highest energies (≳ 100 GeV).

Zoom In Zoom Out Reset image size

The base function for the PSF is the same as used by the XMM-Newton mission (Kirsch et al. 2004; Read et al. 2011):

$$\begin{equation} K(x,\sigma ,\gamma) = \frac{1}{2\pi \sigma ^{2}} \left(1-{\frac{1}{\gamma }}\right) \cdot \left[ 1+{\frac{1}{2\gamma }}\cdot {\frac{x^2}{\sigma ^2}}\right]^{-\gamma }, \end{equation} \tag{ 36 }$$

which Kirsch et al. (2004) refer to as a King function (King 1962), and is isomorphic to the well-studied Student's t-distribution (Student 1908). Note that this function is defined so as to satisfy the normalization condition:

$$\begin{equation} \int ^{\infty }_{0} K(x,\sigma ,\gamma) \, 2\pi x \,dx = 1; \end{equation} \tag{ 37 }$$

the extra 2πx comes from the integral over the solid angle dΩ = sin (x) dxdϕ ∼ 2πx dx. However, at very low energies the PSF widens to the point that the small angle approximation fails by more than a few percent and K(x)sin (x) must be normalized numerically.

To allow for more accurate descriptions of the tails of the distributions, we use the sum of two King functions to represent the dependence of the PSF on scaled deviation for a given incidence angle and energy:

$$\begin{equation} P(x) = f_{\rm core}K(x, \sigma _{\rm core}, \gamma _{\rm core}) + (1 - f_{\rm core}) K(x, \sigma _{\rm tail}, \gamma _{\rm tail}).\qquad \end{equation} \tag{ 38 }$$

The σ and γ values are stored in tables of PSF parameters as SCORE, STAIL, GCORE and GTAIL respectively. Because of the arbitrary normalization used in fitting the PSF function, f_core must be extracted from the NTAIL table parameter, in conjunction with SCORE and STAIL:

$$\begin{equation} f_{\rm core} = \frac{1}{1 + \texttt {NTAIL} \cdot \texttt {STAIL}^2/\texttt {SCORE}^2}. \end{equation} \tag{ 39 }$$

The fitting function has been revised several times since the development of the first preliminary response functions. The version described here is the one currently being used and is different from e.g., that used for P6_V3 IRFs. A description of the fit functions used in the past is given in Section 6.1.2.

The first set of publicly released IRFs, P6_V3, used a slightly different PSF parameterization. Specifically, it allowed for only one σ parameter and fixed the relative normalization of the two King functions by constraining the two to contribute equally at $x_b = 2\sqrt{5}\sigma$. So for P6_V3 we used

$$\begin{equation} P(x)= f_{\rm core} K(x, \sigma , \gamma _{\rm core}) + (1 - f_{\rm core}) K(x, \sigma , \gamma _{\rm tail}), \end{equation} \tag{ 40 }$$

with

$$\begin{equation} f_{\rm core} = \frac{1}{1 + K(x_{b},\sigma ,\gamma _{\rm core}) / K(x_{b},\sigma ,\gamma _{\rm tail})}. \end{equation} \tag{ 41 }$$

In the 1–10 GeV energy range bright pulsars are excellent sources for evaluating the in-flight PSF: not only they are among the brightest γ-ray sources, providing abundant statistics, the pulsed emission is very stable and the angular distributions of γ rays around the true positions can be estimated readily by phase selecting γ rays. The Vela pulsar is the brightest pulsar and an analysis with adequate statistics can be based on Vela alone. In the remainder of this section, we describe the procedure for Vela only; the extension to an arbitrary number of pulsars is straightforward (see Ackermann et al. 2012e).

We use the P7SOURCE Vela calibration sample, which we divided into front- and back-converting subsamples, and bin the data in 4 energy bins per decade. We then calculate the pulsar phase for each γ ray. As mentioned in Section 3.6.1, we define [0.12, 0.17]∪[0.52, 0.57] as the "on" interval and [0.8, 1.0] as "off" (see Figure 15 for the phase histogram). Next we calculate the containment angles: the position of Vela is known with a precision that greatly exceeds the LAT angular resolution, which can be assumed as the true source position. We make histograms of the angular deviations from Vela for the on-peak and off-pulse intervals and normalize them for the relative phase ranges. To estimate the PSF from flight data, we measure the containment radii from the difference between the histograms.

Above ∼10 GeV spectral cutoffs of pulsars leave AGNs as the only attractive sources for studying the PSF. In Ackermann et al. (2012e), we address the potential contribution from pair halos around AGN and conclude that we see no indication that this phenomenon is the explanation for the PSF being broader than predicted. Thus, we treat AGNs as point sources. Many γ-ray sources are considered to be only "associated" with AGNs, as opposed to "firmly identified" (see Nolan et al. 2012, for a discussion of the distinction), because of the limited angular resolution of the LAT. In the present analysis we consider only sources with high-confidence associations. As for pulsars, the positions of AGN are known with high precision from other wavelengths, and the angular distances from the true directions can be calculated. To accumulate enough statistics we stacked several sources and performed a joint analysis. We selected AGNs from among the LAT sources with the highest significances above 10 GeV outside the Galactic plane. This energy limit is set by the source density: below a few GeV the LAT PSF is broad enough that nearby sources frequently overlap.

A significant difference with respect to the pulsar analysis is the necessity of modeling the background in evaluating the distribution of angular deviations. We assume that after the stacking of the sky regions far from the Galactic plane the background count distribution can be assumed to be isotropic. At each energy the background is modeled as a flat distribution normalized by the amplitude in an annulus centered on the stacked data set. The inner radius of this annulus was chosen to be significantly larger than the region containing γ rays from the stacked AGN sample. The uncertainty of the containment radius in each energy bin was set to the rms of a large sample of MC realizations for the signal and background distributions.

We have developed a procedure (described in detail in Ackermann et al. 2012e), to fit our PSF model to the measured containment radii for different energy ranges. Given the statistical limitations we use a single King function, Equation (36). For the same reason we do not measure the dependence of the PSF on the incidence angle, i.e., we calculate an acceptance-weighted average over the incidence angle. We first fit the experimental 68% and 95% containment radii (R₆₈ and R₉₅) with Equation (36). Then we extract a new scaling relation. And finally, we use the fitted (rather than the measured) 68% and 95% containment radii to obtain a new set of PSF parameters for each energy bin. By using the fitted containment radii, this procedure smooths out the statistical fluctuations across the energy bins.

The 68% and 95% angular containment radii for the flight-based P7SOURCE_V6 PSF are shown in Figure 57.

Zoom In Zoom Out Reset image size

We can create custom PSFs in a similar manner to that used for the effective area (see Section 5.7), though with the additional complication that we have to explore variations in the shape of the PSF as a function of energy and incidence angle. Because the dependence of the PSF on γ is not intuitive we choose to express the bracketing functions in terms of the observable quantities R₆₈ and r = R₉₅/R₆₈ rather than in terms of σ and γ. Specifically, in each bin of energy and incidence angle, we can define the bracketing values R'₆₈ and r' in terms of R₆₈ and r:

$$\begin{eqnarray} R^{\prime }_{68} &= R_{68} (1 + \epsilon _{68}(E)B_{68}(E)),\nonumber \\ r^{\prime } &= r (1 + \epsilon _{r}(E)B_{r}(E)). \end{eqnarray} \tag{ 44 }$$

We can then solve for the King function parameters σ' and γ' which would correspond to these values.

We re-analyzed both the B2 1520+31 and the PG 1553+113 ROIs with each of the PSF bracketing functions listed in Table 14 using the procedure described in Section 5.7. Based on the quality of the fits described in Section 6.3, in particular the residuals on the 68% and 95% containment radii, we have assigned

$$\begin{eqnarray} \epsilon _{68}(E) &= 10\%, \nonumber \\ \epsilon _{r}(E) &= 50\%. \end{eqnarray} \tag{ 45 }$$

Table 14. Bracketing PSFs and the Energy-dependent Scaling Functions Used to Create Them

Name	B68(E)	Br(E)
c_scalehi_t_nom	+1	0
c_scalelo_t_nom	−1	0
c_pivothi_t_nom	$\tanh (\frac{1}{k}\log (E/E_{0}))$	0
c_pivotlo_t_nom	$-\tanh (\frac{1}{k}\log (E/E_{0}))$	0
c_nom_t_scalehi	0	+1
c_nom_t_scalelo	0	−1
c_nom_t_pivothi	0	$\tanh (\frac{1}{k}\log (E/E_{0}))$
c_nom_t_pivotlo	0	$-\tanh (\frac{1}{k}\log (E/E_{0}))$

Note. As in Section 5.7 we use k = 0.13, which corresponds to smoothing over ΔE/E ∼ 0.30.

Download table as:  ASCIITypeset image

Tables 15 and 16 show the fit results for PG 1553+113 and B2 1520+31 as well as the integral counts and energy fluxes between 100 MeV and 100 GeV obtained using these PSF bracketing functions. The ranges of the fit values indicate propagated uncertainties from the uncertainty in the PSF (Table 17).

The greater influence of the uncertainty of the PSF on the flux and spectral measurements for the softer source (B2 1520+31) comes about because at lower energies the wider PSF makes resolving sources more difficult and results in greater correlation with the Galactic and isotropic diffuse components in the likelihood fit, which varied by up to ±4% and ±5%, respectively.

For sufficiently long exposures, the LAT can spatially resolve a number of γ-ray sources. In the 2FGL catalog, 12 spatially extended LAT sources were identified (Nolan et al. 2012) using 24 months of LAT data; and several additional extended sources have been recently resolved from these data by Lande et al. (2012) using special techniques for modeling the spatial extension. Understanding the possible spatial extension of LAT sources is important for identifying multiwavelength counterparts, and using a source model with the correct spatial extent produces more accurate spectral fits and avoids biases in the model parameters.

In addition to using a correct spatial model, the accuracy of the PSF can also affect the analyses of extended sources; and using an incorrect PSF will result in biases both in the fitted model parameters and in the significance of any spatial extension that is found. For example, flight data indicate that the MC-based PSF in the P7SOURCE_V6MC IRFs is too narrow at energies >3 GeV (Section 6.2). Fitting the extension of SNR IC 443 (Abdo et al. 2010g) with a uniform disk for data from the first two years of observations and using the flight-determined PSF in P7SOURCE_V6, we find a best fit disk radius of IC 443 of σ = 035 ± 001 (Lande et al. 2012). By contrast, fitting these same data, but using the MC-based PSF in P7SOURCE_V6MC, we find a best fit radius of σ = 039. This corresponds to a ∼10% systematic bias in the measurement of the extension of IC 443. Since IC 443 is a fairly hard source, with photon index Γ = 2.2, the bias found for softer sources would be somewhat smaller.

Use of an imperfect model of the PSF can also affect the assessment of the statistical significance of the extension of the source. We determine the significance of a measured source extension using the test statistic

$$\begin{equation} \mathrm{TS}_\mathrm{ext} = 2\log (L_\mathrm{ext}/L_\mathrm{pt}). \end{equation} \tag{ 46 }$$

This is twice the difference in log-likelihood found when fitting the source assuming it is spatially extended versus assuming it is a point source. The use and validity of this formula is described in Lande et al. (2012). For IC 443, using the P7SOURCE_V6 IRFs, we find TS_ext = 640, corresponding to a formal statistical significance of ∼25σ. Performing the same analysis using P7SOURCE_V6MC, we obtain TS_ext = 1300, corresponding to ∼37σ.

Neglecting the dependence of the PSF on incidence angle (Section 6.3) may introduce time-dependent biases of the estimated flux of the source, due to the different distribution of the source position with respect to the LAT boresight in each time interval. To estimate the potential size of such biases we consider how great an effect on source fluxes we might see if we were to naively use the nominal containment radius for aperture photometry analysis (i.e., an analysis where we simply count the γ rays within a given aperture radius). Specifically, we use the MC-based PSF to extract the containment radii averaged over the first two years of the mission $\bar{R}(E)$ in the direction of the Vela pulsar. Then we split the data into much shorter time intervals (indexed by i) and compute the fraction of γ rays C_i(E; vela) that fall within $\bar{R}(E;{\rm vela})$ for each of those shorter intervals:

$$\begin{eqnarray} &&C_i(E,{\rm vela}) \nonumber\\ &&\quad = \frac{\int \int ^{\bar{R}(E;{\rm vela})}_{0} A_{\rm eff}(E,\theta ; {\rm vela}) t_{\rm obs,i}(\theta ;{\rm vela}) \sin \theta d\theta dr}{\int A_{\rm eff}(E,\theta) t_{\rm obs,i}(\theta ;{\rm vela}) \sin \theta d\theta }. \nonumber\\ \end{eqnarray} \tag{ 47 }$$

Here we consider 68% and 95% containment radii. Figure 64 shows how C_{68, i}(E; vela) and C_{95, i}(E; vela) vary for each 12 hr time interval over the first 700 days of routine science operations. If there were no θ dependence to the PSF this figure would show lines at 0.68 and 0.95. However, we clearly see that ignoring the θ dependence will cause slight misestimates in the fraction of γ rays falling within the mission-averaged $\bar{R}_{68}$ and $\bar{R}_{95}$. In fact, the width of the bands (i.e., the spread of that misestimation) indicate the errors in flux we would expect if we were making a flux estimate based on aperture photometry with aperture cuts at $\bar{R}$. It is worth pointing out explicitly that the C_{95, i}(E) band is significantly narrower, simply because most of the γ rays are already contained within the $\bar{R}_{95}(E)$ and the derivative of the PSF is smaller.

Zoom In Zoom Out Reset image size

Understanding the effect of the PSF dependence on incidence angle will have on a full likelihood analysis is much more complicated, as it depends on the other sources in the likelihood model. However, it is reasonable to take the width of the C₆₈(E) bands as indicative of the magnitude of the effect when analyzing sources in complex regions where source proximity is a issue and the width of the C₉₅(E) bands as indicative when analyzing isolated sources.

The success of scaling the PSF distributions (see Section 6.1.1) prompted a similar approach to the histogramming and fitting of the energy deviations. The scaling function S_D currently employed, derived by fitting the 68% containment of the measured fractional energy distribution (E' − E)/E across the entire E–θ plane, is

$$\begin{eqnarray} S_{D}(E,\theta) &=& c_0(\log _{10} E)^2 + c_1(\cos \theta)^2 + c_2\log _{10} E + c_3\cos \theta\nonumber\\ && + c_4\log _{10} E\cos \theta + c_5. \end{eqnarray} \tag{ 48 }$$

Front- and back-converting events are treated separately. The parameters used since the Pass 6 analysis are listed in Table 18. These are included in the IRF files for Pass 6 and Pass 7 event classes.

The scaling function S_D allows us to define a scaled energy deviation x for which much of the energy and angular dependence is already accounted for

$$\begin{equation} x = \frac{(E^{\prime } - E)}{S_{D}(E, \theta) E}. \end{equation} \tag{ 49 }$$

This scaled value is calculated for each simulated event and histogrammed, as shown in Figure 66. Scaling the distribution achieves two main effects: the width of the core of the distribution is almost constant for all energies and incidence angles, and the distribution thus expressed is significantly simpler to parameterize than the unscaled version.

Zoom In Zoom Out Reset image size

The scaled energy deviation, as illustrated by the example in Figure 66, has a well-defined core accompanied by elongated tails. We found we can effectively fit this structure with several piecewise functions of the form

$$\begin{equation} R(x, x_0, \sigma , \gamma) = N\exp \left(-{\frac{1}{2}} \left|{\frac{x-x_0}{\sigma }}\right|^{\gamma }\right). \end{equation} \tag{ 50 }$$

Indeed we model the energy dispersion D(x) within each (E, θ) bin as four piecewise functions:

$$\begin{equation} D(x) = \left\lbrace \begin{array}{@{}ll@{}}N_{L}R(x, x_0, \sigma _{L}, \gamma _{L}) & \quad {\rm if\; } (x - x_0) < -\tilde{x}\\ \ms N_{l}R(x, x_0, \sigma _{l}, \gamma _{l}) & \quad {\rm if\; } (x - x_0) \in [-\tilde{x}, 0]\\ \ms N_{r}R(x, x_0, \sigma _{r}, \gamma _{r}) & \quad {\rm if\; } (x - x_0) \in [0, \tilde{x}]\\ \ms N_{R}R(x, x_0, \sigma _{R}, \gamma _{R}) & \quad {\rm if\; } (x - x_0) > \tilde{x}.\\ \end{array}\right. \end{equation} \tag{ 51 }$$

The values of the split point $\tilde{x}$ and of the four exponents γ of the energy dispersion parameterization in Equation (51) are fixed as specified in Table 19. Moreover, the relative normalizations are set by requiring continuity at x = x₀ and $\vert x - x_0 \vert = \tilde{x}$ and therefore the fit is effectively performed with a total of six free parameters, which are stored in the IRF FITS files: the overall normalization N_r = N_l (NORM), the centroid position x₀ (BIAS), the two core scales σ_r (RS1) and σ_l (LS1), and the two tail scales σ_R (RS2) and σ_L (LS2).

Table 19. Numerical Values of the Split Point $\tilde{x}$ and of the Four Exponents γ in Equation (51) that are Fixed When Fitting the Scaled Energy Deviations

$\tilde{x}$	γL	γl	γr	γR
1.5	0.6	1.6	1.6	0.6

Download table as:  ASCIITypeset image

The energy resolution is a figure of merit which is customarily used to summarize in a single number the information contained in the energy dispersion parameterization. As illustrated in Figure 67, we define the energy resolution as the half-width of the energy window containing 34% + 34% (i.e., 68%) of the energy dispersion on both sides of its most probable value, divided by the most probable value itself. We note that this prescription gives slightly larger values of energy resolution than using the smallest 68% containment window.

Zoom In Zoom Out Reset image size

Figures 68 and 69 show the energy resolution for P7SOURCE events as a function of energy and incidence angle. As mentioned in Section 2.1.2 the energy resolution has a broad minimum between ∼1 and ∼10 GeV, degrading at lower energies due to the energy deposited in the TKR and at higher energies due to the leakage of the electromagnetic shower out the sides and the back of the CAL. Conversely, the energy resolution tends to improve as the incidence angle increases. This is especially true at high energy, where a longer path length in the CAL implies less shower leakage—though front-converting events more than 55° off-axis tend to exit the sides of the CAL and have worse energy resolution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As stated in Section 1, we factorize the IRFs, effectively assuming that the energy and direction measurements are uncorrelated. We use our allGamma MC sample to test this hypothesis. We locate each event within the cumulative PSF and measured energy distributions:

$$\begin{eqnarray} P_{P} &= \int _{0}^{\alpha } P(\alpha ;E,\hat{v}) d\alpha ,\nonumber \\ P_{D} &= \int _{0}^{E^{\prime }} D(E^{\prime };E,\hat{v}) dE^{\prime }. \end{eqnarray} \tag{ 52 }$$

If the IRFs perfectly described the MC sample, both of these distributions should be flat. Furthermore, if there are no correlations between the energy measurement and the direction estimate, the two-dimensional combined distribution should also be flat. In practice this is nearly the case; the correlation coefficient between P_P and P_D is small (|C_{P, D}| < 0.1) across most of the energy and incidence angle range, as shown in Figure 70(a). However, in many bins of (E, θ), we do observe a small excess of counts with P_P ∼ 1 and P_D ∼ 0. Those values correspond to events that are in the tails of the PSF and have a energy estimate that are significantly lower than the true γ-ray energy. To quantify the magnitude of this effect, we considered the fraction of highly anti-correlated events f_a with P_P > 0.98 and P_D < 0.02 as a function of (E, θ). In the absence of correlations, this fraction should be f_a ∼ 4 × 10⁻⁴. In fact, for highly off-axis γ rays in some energy ranges this fraction approaches f_a ∼ 0.02; however, averaged over the FoV, it is in the 0.001–0.0035 range for both front- and back-converting events at all energies, as shown in Figure 70(b).

Zoom In Zoom Out Reset image size

Although we cannot reproduce this analysis with flight data, we have studied the correlations between the energy and direction quality estimators (P_E, see Section 3.3.2 and P_core, see Section 3.3.3) and found similarly small effects.

In summary, averaged over several orbital precession periods any biases caused by the correlation between the energy dispersion and the PSF are negligible compared to other systematic uncertainties we consider in the paper. However, it is certainly a potential contributor to instrument-induced variability (see Section 8.4).

Since a direct calibration of the LAT with a particle beam was impractical for schedule and cost reasons, the LAT Collaboration assembled for this purpose a dedicated Calibration Unit (CU) composed of two complete flight spare towers, one additional CAL module, and five ACD tiles. An intensive beam test campaign was performed on the CU, between 2006 July and November, with the primary goal of validating the MC simulation used to define the event reconstruction and the background rejection. The CU was exposed to beams of bremsstrahlung γ rays, protons, electrons, positrons and pions—with energies ranging between 50 MeV and 280 GeV (Table 20)—produced by the CERN Proton Synchrotron (PS) and Super Proton Synchrotron (SPS). In addition, the CU was irradiated with 1–1.5 GeV nucleon^{−1 12}C and ¹³¹Xe beams in the Gesellschaft für Schwerionenforschung (GSI) facility with the purpose of studying the instrument response to heavy ions.

A complete review of the results of the beam test campaign is beyond the scope of this paper. The reader is referred to Baldini et al. (2007) for a description of the experimental setup and a review of the main results and to Ackermann et al. (2010) for additional related material.

From the standpoint of the energy measurement one of the primary goals of the beam test campaign was to validate the leakage correction algorithms and, ultimately, the ability of our MC simulation to reliably predict the energy resolution at high energy. The tests confirmed that we understand the overall shower development and the energy resolution to better than 10% up the maximum available electron energy, i.e., 280 GeV (Ackermann et al. 2010). As we will see in the next section, this implies that the effect of systematic uncertainties on the energy resolution itself, when propagated to the high-level spectral analysis, is essentially irrelevant in the vast majority of practical cases.

The most significant discrepancy was observed in the raw energy deposited in the calorimeter, as measured values were on average ∼9% higher when compared with simulations (with smaller fluctuations around this value, slightly dependent on the particle energy and incidence angle). The origin of this discrepancy is unknown, and possibly could be attributed to residual environmental effects not properly accounted for, although an imperfect calibration of the CU cannot be excluded. In Section 7.3.3, we will see that we have indications from flight data that we do understand the absolute energy scale to a precision better than 9%.

The calibration of the CAL crystals is the starting point of the energy reconstruction chain and underlies all the subsequent steps (a detailed description of LAT on-orbit calibrations is given in Abdo et al. 2009a). As explained in Section 2.1.2, the large dynamic range of the LAT CAL (2 MeV to 70 GeV per crystal) is achieved by means of four independent chains of electronics per crystal, with different amplifications (or ranges), as summarized in Table 21. In the most common readout mode (the so-called single-range), the highest gain range that is not saturated is selected as the best estimate and converted into a digital signal that is eventually used in the event reconstruction.

What is relevant for the discussion here is the on-orbit calibration of the crystal response, as determined by the crystal light yield and the linearity of the electronics. The lower energy scales are calibrated using primary protons that do not undergo nuclear interactions in the LAT, which are selected through a dedicated event analysis. For each crystal, the most probable value of the energy deposition, corrected for the path length, is compared with the MC prediction and the conversion factor between digital signal and MeV is computed. This procedure was first tested and validated with sea-level CR muons. The high-energy ranges are calibrated in the energy range of overlap between lex1 and hex8 using events collected by a special on-board filter for selecting heavy ions (see Figure 72), read out in full four-range mode (i.e., with all four energy ranges read out at each log end). The nonlinearity of the electronics is characterized across the entire energy range using a dedicated internal charge injection system and is corrected for in the energy measurement.

Zoom In Zoom Out Reset image size

Though we originally intended to use heavy primary nuclei for an independent calibration of the high-energy ranges, the uncertainty in the scintillation efficiency of heavy nuclei in the CsI(Tl) crystals relative to electromagnetic showers makes the heavy ion peaks unsuitable for an independent cross-check of the absolute scale. The magnitude of the effect was measured with 600–1600 MeV nucleon⁻¹ C and Ni beams (Lott et al. 2006), but we are uncertain how to scale the effect to the typical on-orbit energies (∼5 GeV nucleon⁻¹ or greater) with our desired accuracy of a few percent.

We have used protons and Galactic CRs from Be to Fe at incidence angles ranging from on-axis to 60° off-axis to demonstrate that the crystal scintillation signal is proportional to path length for each species—i.e., the crystal response is piecewise linear over factors of two in signal. Because the peaks from Be to O are closely spaced and easily resolved, this method demonstrates that the CAL energy scale is linear over at least the range ∼180 MeV to ∼1500 MeV per crystal, so any residual error in absolute energy scale applies equally over that entire range. Unfortunately, it is difficult to bridge the gap between protons at 60° (∼22 MeV) and He on-axis (∼45 MeV), and between He at 60° (∼90 MeV) and Be on-axis (∼180 MeV), so we cannot demonstrate linearity in this region with this technique. We are exploring alternative CR event selections and geometries to cover these energy regions. Nonetheless, the relative variations of the peak positions (for both non-interacting protons and heavier nuclei) can be measured with high accuracy, which effectively allows us to monitor the time drift of the absolute energy scale and prompt the update of new calibration constants when necessary. The results in Figure 73 are in good agreement with the pre-launch estimates of the crystal light yield attenuation due to radiation damage, which results predominantly from trapped charged particles in the SAA. A fit with the function

$$\begin{equation} s(t) = (p_0 - p_1) + p_1 e^{-\frac{(t - t_0)}{\tau }} \end{equation} \tag{ 53 }$$

yields a time constant τ of the order of ∼2 years, and predicts an overall shift (∼p₀ − p₁) of the energy scale (which can be corrected for) on the order of 4.5% after 10 years.

Zoom In Zoom Out Reset image size

As mentioned in the beginning of this section, there are essentially no astronomical sources with spectral features that are sharp enough and whose absolute energies are known to such a level of accuracy that they can be effectively exploited for an on-orbit validation of the absolute energy scale.

One exception, perhaps unique in the energy range of the LAT, is the narrow 67.5 MeV pion decay line predicted to originate from interactions of primary CRs with the surface of the Moon (Moskalenko & Porter 2007). Unfortunately, this feature is located at the lower end of the LAT energy range. While it is not inconceivable that the LAT possibly could provide the first evidence for the existence of such a line, the limited energy and angular resolution at these energies, together with the brightness of the limb of the Moon in continuum γ rays, make the lunar pion line impractical as an absolute energy calibrator.

A practical alternative is the geomagnetic rigidity cutoff in the CR electron spectrum. The LAT is effectively shielded from CRs by the Earth's magnetic field at a low rigidity. The convolution of this shielding with the primary power-law spectrum results in a peaked spectral feature, whose shape and absolute energy can be predicted with good precision by taking advantage of accurate models of the geomagnetic field that are available (Finlay et al. 2010) and particle tracing computer programs (Smart & Shea 2005). See Section 4.1 for additional details.

Figure 74 shows an example of such predictions in different parts of the Fermi orbit. The count spectra, averaged over the LAT FoV and folded with the LAT energy resolution, are shown in linear, rather than logarithmic, scale in order to emphasize the peaked spectral shape. Comparison of the predicted and measured peak positions is an in-flight measurement of the systematic bias in the energy scale. Since the magnitude of the Earth's magnetic field varies across the Fermi orbit, this approach has the potential to provide a series of calibration points, specifically between ∼6 GeV and ∼13 GeV. The fact that this is the energy range for which the energy resolution of the LAT is the best is beneficial for the measurement itself. Moreover, since both electrons and γ rays generate electromagnetic showers in the detector, they are effectively indistinguishable from the calorimeter standpoint—so that energy measurements for one species directly apply to the other.

Zoom In Zoom Out Reset image size

We have used this approach for an in-flight calibration of the absolute energy scale using one year of data. The details of the analysis are beyond the scope of this paper and are discussed in Ackermann et al. (2012a). The main conclusion is that the measured cutoff energies exceed the predictions by 2.6% ± 0.5% (stat) ± 2.5% (sys) in the range 6–13 GeV.

Here we summarize the results related to the systematic uncertainties of energy measurements with the LAT; in the next section, we will discuss how these impact the high-level spectral analysis.

The results from the CU beam test campaign indicated that the energy resolution, as predicted by our MC simulations, is accurate to within 10% up to the maximum accessible beam energy (i.e., 280 GeV). As we will see in the next section, this implies that this particular systematic effect is negligible in most practical situations.

Though we were not able to identify the reason for the ∼9% discrepancy in the energy deposited in the CU calorimeter measured at the beam tests (Section 7.3.1), we do have compelling indications, from the measurement of the CR electron geomagnetic cutoff (Section 7.3.3), that in fact the systematic error on the absolute energy scale is smaller than that. We stress that this measurement involves the real LAT, with the real flight calibrations, in its on-orbit environment. Although the measurement derived from the CR electron cutoff applies only to a small portion of the LAT energy range, the ∼9% effect measured at the beam test affected essentially the entire range of energy and incidence angle, with a weak dependence on both (see also the remarks in Section 7.3.2 about the additional evidence for the fractional error on the energy scale being constant over a wide energy range). This evidence is further supported by the excellent internal consistency of the analysis measuring the electron and positron spectra separately (Ackermann et al. 2012b). However, it is more difficult to constrain the energy scale at the low and high ends of the LAT energy range, where the energy resolution degrades by a factor of ∼2. Based on the full body of information currently available we conclude that that the energy scale for the LAT is correct to +20%/ − 50% of the energy resolution of the LAT at a given energy. This corresponds to an uncertainty of +2%/ − 5% on energy scale over the range 1–100 GeV, and increases to +4%/ − 10% below 100 MeV and above 300 GeV.

Finally, the measured energies of the Galactic CR peaks are being monitored to gauge the time stability of the absolute energy scale, which we can control at the 1% level by applying calibration constants on a channel-by-channel basis.