Review: LAT Data Analysis

In general, it can be said that there are two types of analyses of LAT Data: a "point" source analysis (e.g., supernova remnant) and an "extended", diffuse source analysis (i.e., source is 2~3^o across); for example, the Large Magellanic Cloud.

Extended Sources. Extended sources are, at a minimum, many light years across and therefore cannot vary over the course of a year (i.e., cannot get brighter or dimmer). As a result, a standard analysis consists of performing a spectral analysis; no light curve is generated.

Before performing a pyLikelihood analysis, you may wish to review some or all of the following topics:

• LAT Data Products
• FITS File Formats (FT1 and FT2)
• FTOOLS
• Event Classes and Instrument Response Functions (IRFs)
• Data Selection
• Lightcurves
• Count Maps
• Likelihood
  • Overview: Performing the actual fit with gtlikelihood
  • Functional Form of the Likelihood
  • Choosing the Data to Analyze – Regions of Interest and Source Regions
  • Model Selection
  • Precomputation of Likelihood Quantities (Livetime Cube)
  • Model Fitting with gtlikelihood
  • Source Detection and Localization
• GRB Spectral Analysis
  • Binned GRB Spectral Analysis
  • LAT Analysis
  • GBM Analysis
  • Joint Fits

  Note: Many of the topics in this section have been taken from David Band's Cicerone and adapted for use in this introduction to the analysis of LAT data..

LAT Data Products

A LAT data analysis requires information about where the LAT was pointing and what was the observing efficiency. This data is provided by the event and spacecraft FITS files, respectively referred to as FT1 and FT2 files.

Origin of the Photon Data

These two filetypes result from the processing of the data downlinked from the Fermi spacecraft (considered to be 'Level 0' data), and are therefore regarded as 'Level 1' data.

LAT Level 1 processing involves reconstructing the interaction of the event in the LAT from the 'hits' in the various detectors, identifying the type of event (e.g., astrophysical photon), and characterizing the event's relevant physical parameters (e.g., direction, energy).

The characterization of an event results in a set of ~200 parameter values that make up the merit ntuple. Since most events are not astrophysical photons, and most of the parameters describing an event are not relevant for data analysis, only a small set of parameters for the counts have been extracted from the event data to form the event file you will normally use.

Note: The full event dataset will always be available for users who wish to examine the full set of parameters of the astrophysical counts, as well as many events that are not included in the event files.

Photon Classification

During reconstruction, the LAT team makes various cuts that will classify the events based on the probability that they result from photons, and on the quality of the reconstruction. The events will be separated into various event classes, and each class will be characterized by its own set of instrument response functions. Note that this reconstruction methodology and event class cuts have evolved, and they are likely to continue to do so.

Event Classes. There are currently three event classes:

• The diffuse class has the smallest point spread function (PSF) and includes the smallest fraction of background counts; by being more restrictive, the diffuse class has the smallest effective area. This class is intended for the study of diffuse emission.

• The source class is a superset of the diffuse class; by loosening the event selection cuts, the effective area is increased at the expense of including a higher fraction of background counts and photons with a larger PSF.

The trade-off between effective area, background and PSF is favorable for analyzing point sources.

• The transient class loosens the selection cuts even further, and is thus a superset of the source class. The effective area is increased, especially below 250 keV, at the expense of a larger PSF and additional background events.

This class is ideal for studying transients on timescales of less than an hour, such as gamma-ray bursts, because almost all the events originate in the source.

Contents of the Event (FT1) Files. For each astrophysical photon the event file contains the following information (the list does not include all possible quantities):

Contents of the Spacecraft (FT2) Files. Spacecraft files contain the following information for 30 second intervals (some intervals may be shorter):

Other Files. In analyzing the LAT data you will use other files, some of which the tools access without your intervention, some of which are intermediate products of the analysis:

FT1 & FT2 Files/Formats

Starting with the signals from different components of the LAT resulting from interactions of charged particles, the LAT team has reconstructed the paths of the electron-positron pair produced when a gamma ray interacts with a tungsten atom in the LAT, and then calculated the gamma ray's arrival time, incident energy, and origin. The resulting data are stored as FITS files.

FITS File Format

The Flexible Image Transport System (FITS) file format was originally developed to provide a standard image file format, but has been expanded to provide standards for many different file types used in astronomy. The files consist of a series of one or more 'Header and Data Units' (HDUs), each of which contains an ASCII header followed by a binary table.

The ASCII header describes the contents of the binary table (e.g., the column names and units); thus, FITS files are largely self-defining. Headers have rows of text consisting of 8 character keywords followed first by the value of the keyword and then by a comment describing the keyword.

The first ('primary') HDU is reserved for images, and often has no binary table and only a simple header (called the primary header) identifying the file (e.g., name, date of creation, mission). Subsequent HDUs are called 'extensions', which tend to contain copious information pertaining to that extension.

There are standard FITS file formats; for example, PHA for a binned spectrum. These standard file formats have required extensions, and required keywords in the headers of these extensions. Similarly, there are standard extensions (with standard keywords)—such as EBOUNDS for storing an energy grid—that can be used in mission-specific file formats.

Note: Fermi FITS files are defined in the Science Data Products File Format Document (GLAST-GS-DOC-0001 in the Fermi ground system document system) and can be found in the project's Science Data Products File Format Document.

FTOOLS

FTOOLS (fcopy, fdump, fmodhead, fplot, and fverify) allow you to examine, copy and manipulate FITS files. The tool 'fv' is a powerful GUI-based utility that should satisfy most needs to examine and modify FITS files. 'fv' will display images, but 'ds9' was designed specifically to display and manipulate images.

Event Classes and Instrument Response Functions (IRFs)

IRFs describe LAT performance in terms of transformation probability from a true physical quantity (i.e., energy and direction of photons) to the corresponding measured quantity. It is important to note that the IRFs depend not only on the instrument itself, but also on the reconstruction algorithms, the background rejection algorithm, and on the selection of events.

At the time this tutorial was written, the latest Instrument Response Function set was P6_v3 (ScienceTools v9r8p2). Pass_6 provides the default cuts for three event classes:

Data Selection

Data selection precedes data exploration. The fundamental LAT data are simple event lists, and data selection involves making cuts on the event lists. Data selections can be made when extracting the LAT event list from the database, and further selections are made using the gtselect tool. Thus, you can extract data for a source spanning a large time range, and use gtselect to break this event list into a series of shorter time ranges. You can extract LAT events from a large spatial area that includes a gamma-ray burst, and after localizing the burst, you can select the counts from a smaller area centered on the burst.

gtselect – When selecting data, cuts are made using the gtselect tool when extracting data from event lists. This tool enables you to extract data for a long span of time and divide the event list into a series of shorter time ranges. LAT events can also be extracted from a large spatial area that includes a GRB. After localizing the burst, counts can then be selected from a smaller area centered on the burst. (See gtselect.)

Note: After selection, use 'fv' and 'ds9' to display count plots as a function of position, time, or energy. When the number of counts becomes excessive, use gtbin to bin the data in time energy, or space; then use fv and ds9 to plot the resulting FITS files.

gtbin – Bins GBM or LAT events list in time, energy, and/or space to produce light curves, spectra, count cubes, or count maps, respectively. (See gtbin.)

Lightcurves

Lightcurves depict source intensity as the detected count rate, not as photon or energy flux. The gtbin tool lightcurve (LC) option requires the name of the input and output files, and then presents three time binning options:

• Linear (LIN) – requires start and stop times in MET.

• Constant signal-to-noise (SNR) ratio – requires start and stop times in MET, plus SNR value.

• Specified time bins (FILE) – requires name of the input FITS file containing the time binning. (Note: Use gtbindef to convert an ASCII file with the times into a properly formatted FITS file.)

The gtbin FITS file output contains the lightcurve, which can now be plotted.

For examples, see

• Explore LAT Data
• Explore the LAT Data (Bursts)
• GBM Gamma-Ray Burst Analysis

Count Maps

Binned count maps for LAT data depict the number of counts in spatial pixels:

• gtbin's 'CMAP' option provides a count map in one energy band, i.e., the energy band is selected when creating the event file input to gtbin, and all events are binned in a single event file.

• gtbin's 'CCUBE' option is used to bin multiple energy bands.

Both the CMAP and CCUBE options require you to specify the:

• range of time for which the event file is to be created, and

• rectilinear spatial grid, including the:
  • number of pixels in each direction,
  • pixel size,
  • coordinate system,
  • center of the grid, and the
  • type of projection.

  Note: Ten projections are provided (see Calabretta & Greisen 2002, A&A, 395, 1077); AIT for 'Aitoff' is suggested.

For the CCUBE option, you must also specify the energy bands to be used:

• Linear binning (LIN) – minimum and maximum energies, and the number of bins.

• Logarithmic (LOG) – minimum and maximum energies, and the number of bins; the software will calculate the bin edges. (probable default)

• User specified energy bins (FILE) – name of the FITS file with the energy binning.

Note: Use gtbindef to convert an ASCII file with the energy bands into the properly formatted FITS file.

Output map files can be viewed by fv or ds9.

For examples, see:

• Explore LAT Data
• Explore the LAT Data (Bursts)

Likelihood

The likelihood that is applicable to the LAT data is constructed and then used to find the best fit model parameters, which includes the description of a source's spectrum, its position, and even whether it exists.

The likelihood (L) is the probability of obtaining the data given an input model. For the purposes of this discussion, the input model is the distribution of gamma-ray sources on the sky, and includes their intensity and spectra.

Overview: Performing the actual fit with gtlikelihood

When performing the actual fit with gtlikelihood, the parameter space can be quite large – the spectral parameters from a number of sources must be fit simultaneously; therefore, the SAE provides a choice of three 'optimizers' to maximize the likelihood efficiently. Fitting requires repeatedly calculating the likelihood for different trial parameter sets until a value sufficiently near the maximum is found; the optimizers guide the choice of new trial parameter sets to converge efficiently on the best set. The variation of the likelihood in the vicinity of the maximum can be related to the uncertainties on the parameters, and therefore these optimizers estimate the parameter uncertainties.

Thus likelihood spectral fitting provides the best fit parameter values and their uncertainties. But is this a good fit? When χ2 is a valid statistic, then we know that the value of χ2 is drawn from a known distribution, and we can use the probability of obtaining the observed value as a goodness-of-fit measure. When there are many degrees of freedom (i.e., the number of energy channels minus the number of fitted parameters) then we expect the χ2 per degree of freedom to be ~1 for a good fit. However, when χ2 is not a valid statistic, we usually do not know the distribution from which the maximum likelihood value is drawn, and therefore we do not have a goodness-of-fit measure.

Model Fitting

When it is known that a source is present, the objective is to determine the best value of the spectral model parameters. To determine the best model (i.e., the one having highest probability of resulting in the data), vary the spectral parameters until the likelihood is maximized.

Note: χ2 is -2 times the logarithm of the likelihood in the limit of a large number of counts in each bin, and therefore where χ2 is a valid statistic, minimizing χ2 is equivalent to maximizing the likelihood.

A number of steps are necessary to fit a source's spectra:

1. Select the data. Data must be used from a substantial spatial region around the source(s) being analyzed due to the overlapping of the point spread functions of nearby sources.

2. Select the model. This model selected must include:
   • the position of the source(s) being analyzed;
   • the position of nearby sources;
   • a model of the diffuse emission;
   • the functional form of the source spectra; and
   • values of the spectral parameters.

   In fitting the source(s) of interest, let the parameters for these sources vary, but because the region around these sources includes counts from nearby sources that are not of interest, you might also let the parameters from these nearby sources vary.

3. Precompute a number of quantities that are part of the likelihood computation. As parameter values are varied in searching for the best fit, the likelihood is calculated many times. Therefore, precomputing a number of computation-intensive quantities will greatly speed up the fitting process.

4. Perform the actual fit. The parameter space can be quite large – the spectral parameters from a number of sources must be fit simultaneously; therefore, the SAE provides a choice of three optimizers to maximize the likelihood efficiently.

Fitting requires repeatedly calculating the likelihood for different trial parameter sets until a value sufficiently near the maximum is found, and the optimizers guide the choice of new trial parameter sets to converge efficiently on the best set.

The variation of the likelihood in the vicinity of the maximum can be related to the uncertainties on the parameters; therefore, these optimizers estimate the parameter uncertainties.

While likelihood spectral fitting provides the best fit parameter values and their uncertainties, the question remains, "Is this a good fit?":

When χ2 is a valid statistic, we know that the value of χ2 is drawn from a known distribution, and we can use the probability of obtaining the observed value as a goodness-of-fit measure. When there are many degrees of freedom (i.e., the number of energy channels minus the number of fitted parameters), we expect the χ2 per degree of freedom to be ~1 for a good fit.

However, when χ2 is not a valid statistic, we usually do not know the distribution from which the maximum likelihood value is drawn; therefore we do not have a goodness-of-fit measure.

Source Localization

The optimizers find the best fit spectral parameters, but not the location. In other words, the fitting tool does not fit the source coordinates. However, a tool is provided that performs a grid search – mapping out the maximum likelihood value over a grid of locations. As explained below (see Source Detection), it is convenient to use a quantity called the 'Test Statistic' TS that is maximized when the likelihood is maximized.

Source Detection

The Test Statistic is defined as TS=-2ln(Lmax,0/Lmax,1), where Lmax,0 is the maximum likelihood value for a model without an additional source (the 'null hypothesis') and Lmax,1 is the maximum likelihood value for a model with the additional source at a specified location. As can be seen, TS is a monotonically increasing function of Lmax,1, which is why maximizing TS on a grid is equivalent to maximizing the likelihood on a grid. In the limit of a large number of counts, Wilkes Theorem states that the TS for the null hypothesis is asymptotically distributed as χ²x (here χ2 is the distribution, not a value of the statistic), where x is the number of parameters characterizing the additional source. This means that TS is drawn from this distribution if no source is present, and an apparent source results from a fluctuation. Thus, a larger TS indicates that the null hypothesis is incorrect (i.e., a source really is present), which can be quantified. Whether the LAT data are in the limit where this assumption applies is under investigation.

Functional Form of the Likelihood

LAT data is binned into a great many bins because the counts are characterized by many variables; therefore, even with many counts, each bin will contain a small number of counts. The observed number of counts in each bin is characterized by the Poisson distribution, and with a small number of counts per bin, the Poisson distribution cannot be approximated by a normal distribution.

The likelihood L is the product of the probabilities of observing the detected counts in each bin. Assume that the expected number of counts in the ith bin is m_i. Note that m_i is a function of the source model, and will differ for different models. The probability of detecting n_i counts in this bin is p_i=m_i^n_i exp[-m_i]/n_i!. The likelihood L is the product of pi for all i. But notice that this product factors into the product of the m_i^n_i/n_i!, which depends on the data (i.e., the values of n_i) and the product of the exp[-m_i]. The product of exp[-m_i] for all i is equal to the exponential of minus the sum of m_i. The sum of m_i is just the total number N_exp of counts that the source model predicts should have been detected.

Therefore, the likelihood L can be factored into exp[-N_exp], which is purely a function of the source model, and the product of m_i^n_i/n_i!, which is a function of both the source model and the data:

L = exp[-N_exp]  ∏i m_i^n_i/n_i!

This likelihood, with finite size bins and n_i that may be greater than 1, is the basis for binned likelihood. Since binning destroys information (i.e., the precise values of the quantities describing a count), there is a tradeoff between the number of bins (and thus the bin size) and the accuracy; smaller bins result in a more accurate likelihood.

If we let the bin sizes get infinitesmally small, then n_i=0 or 1. The likelihood is now the product of exp[-N_exp], as before, and the product of m_i where i is now the index over the counts.

L = exp[-N_exp]  ∏i m_i

Since m_i is calculated using the precise values for each count, and not an average over a finite size bin, this unbinned likelihood is the most accurate.

For a small number of counts, the unbinned likelihood can be calculated rapidly; but as the number of counts increases, the time to calculate the likelihood becomes prohibitive and the binned likelihood must be used.

Choosing the Data to Analyze – Regions of Interest and Source Regions

Analyzing the Spectrum of a Point Source

The data from a substantial spatial region around the source(s) being analyzed must be used because of the overlapping of the point spread functions of nearby sources.

For the greatest accuracy possible in modeling a single source, it would be necessary to model the entire sky. In reality, the influence of sources a very great distance away from your source will be greatly attenuated (i.e., at 100 MeV, 68% of the counts will be within 3.5 degrees due to the large point spread function at low energies). Thus, it is recommended that you include sources from a large 'Source Region', and counts from a smaller 'Region of Interest' (ROI).

Obtain the positions and spectra of sources in the Source Region outside of the ROI from a catalog; these sources are included for their contribution to the counts in the ROI. You may elect to:

• fix the parameters of the sources other than the one you are studying at their catalog values; or

• perform a fit on the parameters of all these sources.

All sources in the source region will be used, and the size of the Source Region appropriate for your needs is determined by your experience and prior experimentation. In the meantime, the recommended default values are ROI+10 and ROI+5 degrees for sources dominated by ~100 MeV and ~1 GeV events, respectively.

All counts in the ROI will also be included and you will also determine the appropriate ROI size based on your experience and prior experimentation. In the meantime, the recommended default values are 20 and 15 degrees for sources dominated by ~100 MeV and ~1 GeV events, respectively.

When you extract the LAT data from the online database, you will be able to choose the time range, spatial region and energy range from which the counts are extracted. The result is a photon file (FT1) with these counts, and a spacecraft file (FT2) describing the position of the spacecraft and the orientation of the LAT over the chosen time range. The spatial region is the Region of Interest (ROI). Further selections can be performed using the gtselect tool.

Note: Data selection criteria are recorded in the photon file with 'Data Sub-Space' (DSS) keywords.

Model Selection

Models used by the SAE are stored in XML files. XML is a language to define and store data; HTML used to define webpages is a form of XML. XML files are ASCII files full of 'tags' delimited by the < and > symbols.

Note: A model includes the position of the source(s) being analyzed, the position of nearby sources, a model of the diffuse emission, the functional form of the source spectra, and values of the spectral parameters. In fitting the source(s) of interest, you will let the parameters for these sources vary, but because the region around these sources includes counts from nearby sources in which you are not interested, you might also let the parameters from these nearby sources vary.

For historic reasons the SAE uses two XML formats, one used for parameter fitting (e.g., by the likelihood tool) and the other for source simulation. The likelihood XML format includes parameter uncertainties but does not allow time dependence, whereas the simulation format does. The simulation format also includes source models that are not in the likelihood format.

Creating and Editing Source Models

A GUI-driven tool called ModelEditor is included in the SAE tools, and is invoked at the command line; its use is fairly intuitive.

While you might want to define each source in your Source Region, it is easier to start with a catalog for the region. After opening a catalog, choose, choose File --> Extract.... A GUI will open that allows you to define your:

• source Region (RA, Dec and radius in degrees), and the
• flux limit,
• file with the catalog, and the
• output file.

Note: Remember, the Source Region should be larger than the Region of Interest.

A GUI with a list of selected sources will then be displayed. When you select a source, you will be presented with the source properties (i.e., a position and a spectral model), which can be edited. For example, the catalog may have assumed the spectrum was a simple power law, but you would like to try a power law with an exponential cutoff.

Note: To save the model for a particular source, click on the Set components button.

Using the pull-down menu, you can:

• delete a selected source,
• add a diffuse, or
• add a point source.

Multiple xml files. The entire source model need not be in the same XML file. The tools accept lists of source XML files, both at the command line and by reading in simple ASCII files listing the XML files.

The model parameters (in the likelihood XML format) have a number of attributes:

Therefore, if the 'value' is '4.3' and the 'scale' is '1e^-9,' the actual parameter value is 4.3x10^-9.

Spatial Models

Select from four models:

Spectral Models

The units for the spectral models are cm^-2 s^-1 MeV^-1 for point sources and cm^-2 s^-1 MeV^-1 sr^-1 for diffuse sources. All energies are in MeV.

See: Likelihood Tutorial: 4. Create a Source Model XML File.

Precomputation of Likelihood Quantities (Livetime Cube)

The computation of the likelihood usually occurs many times as parameter values are varied in searching for the best fit. While not strictly necessary, precomputing a number of computation-intensive quantities will greatly speed up the fitting process. Fitting involves varying model parameters until the best values are found (the methodology is described elsewhere). Fits are done with various model parameters fixed, or with different sources present or absent. Certain quantities need be calculated only once, speeding up the repeated computation of the likelihood.

Livetime Cube

Because livetime cube calculation is computationally intensive, the livetime cubes are provided on different timescales. Thus to analyze data from a given time period, download livetime cubes spanning most of the period, run gtltcube for the time not covered by these pre-packaged livetime cubes, and then sum all these livetime cubes with gtltsum.

The number of counts that a source should produce should therefore depend on the amount of time that the source spent at a given inclination angle during an observation. This livetime quantity, the time that the LAT observed a given position on the sky at a given inclination angle, depends only on the history of the LAT's orientation during the observation, and not on the source model.

The array of these livetimes at all points on the sky is called the 'livetime cube' and, as a practical matter, the livetime cannot be provided as a continuous function of inclination angle or position on the sky.

Thus, the SAE livetimecubes are provided on a healpix grid on the sky, and in inclination angle bins.

• Size of the spatial grid (in degrees).
• The inclination angle binning (in cos τ steps).
• The spacecraft file (FT2).

gtltcube calculates the livetime cube for the entire sky for the time range covered by the spacecraft file; therefore, the same output file can be used for analyzing different regions of the sky over the same time range.

Note: The livetime cube is calculated on a spatial healpix grid.

Livetime Cubes Are Additive. The livetime cube for a time range that is the sum of two time subranges is the sum of the livetime cube for each of the time subranges.

Thus, the livetime cube for five calendar days can be calculated by adding the livetime cubes for each of the five calendar days. Livetime cubes can be added by gtltsum.

Exposure Maps

Exposure maps are used for extended sources such as the diffuse Galactic and Extragalactic backgrounds, not for individual sources.

Remember, a likelihood analysis considers the counts in the Region of Interest resulting from all sources in a Source Region that is a superset of the Region of Interest. Therefore the exposure map must extend over the entire Source Region, and is specific to the Region of Interest.

The likelihood consists of two factors:

• The first is dependent on the detected counts and differs between binned and unbinned likelihood calculations.

• The second is equal to the exponential of minus the expected total number of counts Nexp for the source model.

The exposure map is the total exposure – area multiplied by time – for a given position on the sky, producing counts in the Region of Interest. Since the response function is a function of the photon energy, the exposure map is also a function of this energy.

Thus, the counts produced by a source at a given position on the sky is the integral of the source flux and the exposure map (a function of energy) at that position.

gtexpmap

The exposure map is calculated by the gtexpmap tool, which derives the Region of Interest from the observation's photon file (FT1); remember that the region from which counts were selected is recorded by keywords in the photon file's FITS header.

Inputs. The Source Region is assumed to be centered on the Region of Interest, but the size of the Source Region must be input. The spatial and energy gridding for the Source Region must also be provided.

The tool also needs the livetime spent at each inclination angle at every point in the Source Region; this can be provided by the livetime cube (as described in the previous subsection); or, if the livetime cube file is not provided, the tool will calculate these livetimes from the spacecraft file (FT2).

Note: The exposure map is calculated on a longitude-latitude grid, while the livetime cube is calculated on a spatial healpix grid.

Diffuse Emission and Individual Counts

The term in the likelihood that is dependent on the counts is a function of the probability of observing each count. This probability is the sum of the probability from all sources, both point and diffuse.

The gtdiffrsp tool calculates the contribution of diffuse sources to this probability.

Notes:

• The photon file (FT1) lists the position of all counts.

• The spacecraft file (FT2) provides the history of the inclination angle for all positions from which the diffuse emission originates.

• The source model describes the diffuse emission.

Model Fitting with gtlikelihood

Fitting involves finding the set of parameters that maximizes the likelihood. (See Overview: Performing the actual fit with gtlikelihood.

The maximum is found by iteratively calculating the function for different sets of trial parameters, and estimating derivatives of the function with respect to the parameters. The algorithms choose new trial parameters that are progressively closer to the set that maximizes the function.

Notes:

• The function is calculated for new sets of trial parameters, until the change in the function value between iterations is sufficiently small (or the number of iterations reaches a maximum value).

• While iterating, these algorithms map out the dependence of the function on the parameters, particularly near the function's maximum. The uncertainties on the best fit parameters are related to this dependence.

• Different algorithms vary in how rapidly they converge to the function maximum, the amount of computer memory they require, and the accuracy with which they map out the dependence of the function on the parameters near the maximum (and thus estimate the uncertainty).

Inputs necessary for fitting parameters with gtlikelihood:

• An event file with the counts to be fit. (See Choosing the Data to Analyze - Regions of Interest (ROIs) and Source Regions.)

• A spacecraft file that covers the time range over which the counts in the event file were extracted.

• A source model. (See Model Selection.)

• Precomputed quantities, such as the livetime cube and the exposure map. (See Precomputation of Likelihood Quantities (Livetime Cube.)

When running gtlikelihood,you have a choice of algorithms, called optimizers, for maximizing the likelihood. (See Overview: Performing the actual fit with gtlikelihood.)

A reasonable strategy is to run gtlikelhood with the DRMNFB optimizer until convergence, and then to run gtlikelhood with NEWMINUIT with the best fit parameter values from this first run in order to calculate the uncertainties on the parameter values.

DRMNFB is efficient at finding the maximum likelihood but approximates the parameter dependence near this maximum; consequently, the uncertainties provided by this optimizer may not be reliable.

On the other hand, NEWMINUIT is a conservative optimizer that converges more slowly than these other methods; indeed, it may exhaust the number of permitted iterations before convergence. However, NEWMINUIT more accurately maps out the parameter space near the likelihood maximum, thus providing more reliable uncertainty estimates.

Convergence criteria is controlled by the hidden variable, fit_tolerance, whose definition depends on the particular optimizer, but is approximately the fractional change in the logarithm of the likelihood.

Source Detection and Localization

Currently, the likelihood tool, gtlike does not fit the source position parameters. Instead, the gttsmap tool runs gtlike at each gridpoint on a rectangular grid. The user provides a model of all the other sources in the Source Region, and gtssmap calculates the Test Statistic (TS) for adding an additional source at each gridpoint. A livetimecube and exposure map precomputed for running gtlike can also be used for gttsmap.

If you have a good approximate source location, such as the gridpoint with the maximum TS from gttsmap or a candidate source identification, then you can refine the location by running gtfindsrc, which maximizes the TS in continuous space.

The TS is -2 times the logarithm of the ratio of the likelihood for the model without the additional source (the null hypothesis) to the likelihood for the model with the additional source. Thus the TS is maximized when the likelihood for the model with the source is maximized; therefore, the location with the maximum TS is the best fit source position.

But is the addition of this source significant? By Wilks' Theorem, if there is no additional source then the TS should be drawn from an χ²n distribution, where n is the difference in the degree of freedom between the models with and without the additional source. In the case under consideration, the additional source is characterized by a source intensity and spectral index (the spectrum is assumed to be a power law); thus n=2. Wilks' Theorem is valid asymptotically as the number of counts increases. At the time this was written, studies were underway to determine if the number of LAT counts for a typical analysis is sufficiently large. If Wilks' Theorem is valid, then integrating χ²2 from the observed TS value to infinity gives the probability that the apparent source is a fluctuation. The resulting significance is ~(TS)^1/2σ, and thus TS=25 is equivalent to 5σ.

Overview: GRB Spectral Analysis

The duration of the prompt burst emission – the ~100 keV component – is (relatively) short; at most 10's of seconds. Therefore, the LAT's pointing will not change significantly during the burst, and all the counts can be treated as having one response function. Within a PSF radius of the burst position, less than one non-burst count per minute is expected. The count rate over the FOV is ~2 Hz; the LAT's FOV is approximately 2 steradians; and the PSF has a radius of ~3.5 degrees at 100 MeV. Consequently, we expect ~0.01 Hz non-burst photons, or 0.7 cts/minute within a PSF radius. Therefore, we can treat all counts within 1-2 PSF radii as burst photons.

Multi-source spatial analysis is unnecessary for spectral analysis of LAT data, since:

1. all counts within a PSF radius of the burst originated in the burst, and

2. all the counts have the same response function.

However, spatial analysis might be necessary for localizing the burst.

All of the counts within a PSF radius, and within a given time range, can be binned into a count spectrum (apparent energy is the single dimension); and traditional spectral analysis can be applied to the resulting series of LAT count spectra.

Since both the LAT and GBM data consist of lists of counts, the same temporal binning for burst data from both detector types can be selected; joint fits on the binned one dimensional spectra can then be performed.

Binned GRB Spectral Analysis

Method

GRB spectral analysis takes advantage of the unique properties of the phenomenon of a relatively short point source transient. Normally, the LAT there will not experience competing emission from other sources in the field-of-view. Therefore, the analysis is one dimensional: the input burst flux is determined from the apparent energies of the events that triggered the detector which, for the purposes of this discussion, these events will be called 'counts'. It is also assumed that there are sufficient counts per bin.

For the LAT all the counts from a region 1-2 PSF-radii around the burst position are selected from the time range which includes the burst; these counts are then binned into energy channels.

Notes:

• For the LAT, these counts should all originate from the burst because estimates of the non-burst event rate predict about one count within a PSF-radius per minute.

• For the GBM, the origin of the counts is unknown, and therefore all the counts from the burst's time range are selected. The counts consist of photons originating from the burst and background from other sources, both astrophysical and instrumental. The background is usually estimated from the count rates before and after the burst. The GBM will have already binned the counts into predetermined energy channels.

• The selected counts may be further binned in time. The result is then a series of count spectra that will be analyzed.

In general each count spectrum is fitted independently.

Consider a count spectrum c_i, where the index i runs over energy channels. This count spectrum is the sum of the burst flux convolved with the detector response and the background b_i. We sample the photon flux striking the detector in different energy channels f_j, where the index j runs over energy channels (NOT necessarily the same channels as the count spectrum!). The response function can be simplified into a mapping between the photon's true energy and the count's apparent energy. With the counts and the fluxes expressed as vectors, the response function is a matrix D_ij, the 'Detector Response Matrix' (DRM) in the burst community, or the 'RSP' or 'RMF' in the X-ray astrophysics community. The resulting matrix equation is

c_i = D_ijf_j+b_i

where summation over j is assumed. Since D_ij is not a square matrix, and even if it is, it is usually nearly singular; this equation cannot be solved by inverting D_ij but requires 'forward folding.' Note that for the LAT b_i~0, but for the GBM b_i is substantial, and in many channels will dominate the burst counts.

In forward folding a model flux vector f'_j is folded through the response, resulting in a model count spectrum c'_i. The underlying model flux is usually an analytic function (e.g., a power law) with a small number of spectral parameters (e.g., normalization and spectral index for a power law flux model). The model c'_i is compared to the observed c_i, and then a new model flux vector f'_j is calculated, usually by varying the spectral parameters. This iterative process ends when the model c'_i is sufficiently close to the observed c_i, resulting in best-fit spectral parameters.

'Sufficiently close' is usually determined by minimizing chi2. A sufficiently small value of chi², e.g., comparable to the number of degrees-of-freedom, indicates that the fit is satisfactory. If the number of counts per bin is not large enough to assume that they are drawn from a Gaussian distribution, then the Cash statistic should be used instead of chi².

If two or more detectors observe the same burst at the same time, then the counts recorded by each detector resulted from the same input burst spectrum. Thus, we can require that the count spectra for each detector be fit by the same flux model. The result is a joint fit.

LAT Analysis

The LAT data consist of photons. To use the techniques described above, these events must first be binned. The steps in the analysis are as follows:

1. Extract the photons from a region around the burst at the time of the burst.

The selection of the photons occurs both in the extraction of the data from the database, and in using the gtselect tool.

2. Bin the photons with gtbin.

Energy bins are either chosen interactively (bins that are uniform in photon energy or the logarithm of the photon energy), or through a binning file.

Time bins can be also be chosen interactively, or through a binning file, and can be:

• uniform in time.
• have a constant signal-to-noise ratio per bin, or be
• Bayesian Blocks.

Note: The binning files are created either by a previous run of gtbin, or by using the gtbindef utility. The output of the binning is a standard PHA or PHAII .fits file.

3. Create the DRM with gtrspgen.

The output is a standard RSP .fits file.

4. Fit the resulting spectra with an analysis tool such as XSPEC. The inputs are the PHA and RSP files created above. (No background file!)

GBM Analysis

The GBM data also consists of individual events (here called 'counts'). Once again, they must be binned.

• Bin the counts with gtbin.

The GBM detectors will have already binned the counts in energy; time bins can be chosen interactively, or through a binning file.

Three interactive time binning methods are available: uniform in time; constant signal-to-noise ratio per bin; or Bayesian Blocks.

The binning files are created either by a previous run of gtbin or using the gtbindef utility. The output of the binning is a PHA or PHAII file, standard FITS filetypes.

• Fit the resulting spectra with an analysis tool such as XSPEC.

The input are the PHA or PHAII file created above, and the RSP and background files provided with the burst data.

Joint Fits

A major hurdle for joint fitting has always been getting spectra from different detectors with the same time bins; however, because the Fermi data are event lists, data can be binned with the same time bins.

The binning tool gtbin can output a file with the time bins used to bin an event list, and can read a binning file to bin an event list. Therefore:

• Bin the data from one detector; for example, using constant signal-to-noise ratio bins.

These time bins should be output in a binning file.

• Use the resulting binning file to bin the event data from the other detectors.

• Fit the spectra simultaneously with an analysis tool such as XSPEC; most such tools are capable of performing joint fits.

Note: The relative normalization between detectors should normally be added as a fitted parameter to compensate for errors in the effective areas of the different detectors.