Given the dearth of bright sources in the extraction region we have selected, our source model file will be fairly simple, comprising only the Galactic and Extragalactic diffuse emission, and point sources to represent the blazars 3C 279 and 3C 273. Here is an example XML file:

The XML file used for this example is 3C279input_model.xml (in the cell below). For more details on the available XML models, see:

• Descriptions of available Spectral and Spatial Models
• Examples of XML Model Definitions for Likelihood

XML for Extended Sources

If your region includes any extended sources, you will need to add one or more extended sources to your XML model.

Download the extended source templates from the LAT Catalog page (look for "Extended Source template archive"). Extract the archive in the directory of your choice and note the path to the template files, which have names like W44.fits and VelaX.fits. You will need to modify your XML model so that the path to the template file is correct.

In addition, you will need to add the modifier map_based_integral='true' to the <spatialModel> atttribute. This modifier tells gtlike to use the proper integration method when evaluating the extended source.

Here is an example of the proper format for an extended source XML entry for Unbinned Likelihood analysis:

6. Compute the diffuse source responses

The diffuse source responses are computed by the gtdiffrsp tool.

The source model XML file must contain all of the diffuse sources to be fit. gtdiffrsp will add one column to the event data file for each diffuse source. The diffuse response depends on the instrument response function (IRF), which must be in agreement with the selection of events, i.e. the event class and event type we are using in our analysis. In this case we are using the standard SOURCE class (evclass=128) and event type 3 (FRONT+BACK), so we should use the IRF P8R3_SOURCE_V3 (which is the most commonly used in standard analysis).

Alternatively, you can just use CALDB for the IRF and the appropriate IRF will be automatically applied according to the previous event selection you performed during the previous steps of the analysis.

If the diffuse responses are not precomputed using gtdiffrsp, then the gtlike tool will compute them at runtime (during the next step). However, as this step is very computationally intensive (often taking ~hours to complete), and it is very likely you will need to run gtlike more than once, it is probably wise to precompute these quantities.

Here is an example of running gtdiffrsp:

Note: If trying to use CALDB in the example above does not work, just put the IRF name in explicitly. The appropriate IRF for this data set is P8R3_SOURCE_V3.

There are a number of IRFs distributed with the Fermitools. Some of these IRFs, like P8R3_SOURCE_V3::FRONT, are designed to address only a subset of the events in the dataset, and thus would require you to make additional cuts before they could be used.

Once you have run gtdiffrsp, you can find the current names of the diffuse response column using the FTOOLS with:

For example:

The first part tells you the IRF name and, after the double underscore, the name of the diffuse component associated with the response.

If you are using a different IRF or diffuse model, you MUST compute a custom diffuse response prior to continuing the analysis.

Now we are ready to perform a likelihood analysis. You can get the resulting file, including the diffuse responses, here.

7. Run gtlike

We are now ready to run the gtlike application.

You may need the results of the likelihood fit to be output in XML model form (e.g. to use in generating a test statistic map).

To obtain an XML model output in addition to the standard results files, use the sfile parameter on the command line (as shown below) to designate the output XML model filename.

Most of the entries prompted for are fairly obvious.

In addition to the various XML and FITS files, the user is prompted for a choice of IRFs, the type of statistic to use, the optimizer, and some output file names.

The statistics available are:

• UNBINNED: This is a standard unbinned analysis, described in this tutorial, to be used for short timescale or low source count data. If this option is chosen then parameters for the spacecraft file, event file, and exposure file must be given.

• BINNED: This analysis is used for long timescale or high-density data (such as in the Galactic plane) which can cause memory errors in the unbinned analysis. See an explanation in the Binned Likelihood Tutorial.

There are five optimizers from which to choose: DRMNGB, DRMNFB, NEWMINUIT, MINUIT and LBFGS.

Generally speaking, the faster way to find the parameter estimates is to use DRMNGB (or DRMNFB) to find initial values and then use MINUIT (or NEWMINUIT) to find more accurate results.

If you have trouble achieving convergence at first, you can loosen your tolerance by setting the hidden parameter ftol on the command line. (The default value for ftol is 0.01.)

The application proceeds by reading in the spacecraft and event data, and if necessary, computing event responses for each diffuse source. If you get an error similar to this:

Caught St11range_error at the top level: Requested energy, 911163, 
lies outside the range of the input file, 34.171, 877938

This means that we are trying to do our analysis beyond the maximum energy covered by the data sample. We should have checked that it was not the case during the data selection (as recommended in step 1).

To solve this issue we need to apply a harder cut in our event selection, setting the maximum energy within the energy range of available photons in your data sample.

Fortunately, this new event selection will not affect the time calculation (previous steps do not need to be repeated).

Here is the output from our fit with gtlike:

Minuit did successfully converge.

(MUCH OUTPUT OMITTED.)

Computing TS values for each source (4 total)
....!

Photon fluxes are computed for the energy range 100 to 500000 MeV

3C 273:
Prefactor: 9.53263 +/- 0.262701
Index: -2.63212 +/- 0.0212044
Scale: 100
Npred: 6369.65
ROI distance: 10.4409
TS value: 7810.96
Flux: 5.85022e-07 +/- 1.14156e-08 photons/cm^2/s

3C 279:
Prefactor: 7.3524 +/- 0.189805
Index: -2.20869 +/- 0.0144918
Scale: 100
Npred: 7420.98
ROI distance: 0
TS value: 13068.8
Flux: 6.08821e-07 +/- 1.08338e-08 photons/cm^2/s

gll_iem_v07:
Prefactor: 1.20506 +/- 0.0160704
Index: 0
Scale: 100
Npred: 95322.2
Flux: 0.000627754 +/- 8.37036e-06 photons/cm^2/s

iso_P8R3_SOURCE_V3:
Normalization: 1.05663 +/- 0.0281431
Npred: 47447.5
Flux: 0.000129821 +/- 3.45716e-06 photons/cm^2/s

WARNING: Fit may be bad in range [100, 549.28] (MeV)
WARNING: Fit may be bad in range [213340, 326604] (MeV)

Total number of observed counts: 156560
Total number of model events: 156560

-log(Likelihood): 1670346.65

Writing fitted model to 3C279output_model.xml

Since we selected plot=yes in the command line, a plot of the fitted data appears. (Note: Certain configurations are not producing a plot after gtlike. The FSSC is working on this issue at the moment. Please click here to generate a plot within python.)

In the first plot, the counts/MeV vs MeV are plotted. The points are the data, and the lines are the models. Error bars on the points represent sqrt(Nobs) in that band, where Nobs is the observed number of counts. The black line is the sum of the models for all sources.

The colored lines follow the sources as follows:

• Black - summed model
• Red - first source in results.dat file (see below
• Green - second source
• Blue - third source
• Magenta - fourth source

In our case the 3C273 is the one in red, 3C279 is the one in green, the extragalactic background is blue, and the galactic background is magenta. If you have more sources, the colors are reused in the same order.

The second plot gives the residuals between your model and the data. Error bars here represent (sqrt(Nopbs))/Npred, where Npred is the predicted number of counts in each band based on the fitted model.

To assess the quality of the fit, look first for the words at the top of the output <Optimizer> did successfully converge. Successful convergence is a minimum requirement for a good fit.

Next, look at the energy ranges that are generating warnings of bad fits. If any of these ranges affect your source of interest, you may need to revise the source model and refit. You can also look at the residuals on the plot (bottom panel). If the residuals indicate a poor fit overall, you should consider changing your model file, perhaps by using a different source model definition and/or adding new sources in the ROI, and then refit the data.

If the fits and spectral shapes are good, but could be improved, you may wish to simply update your model file to hold some of the spectral parameters fixed.

For example, by fixing the spectral model for 3C 273, you may get a better quality fit for 3C 279. Close the plot and you will be asked if you wish to refit the data:

Refit? [y]n 
Elapsed CPU time: 111.819315
prompt>

Here, hitting "return" will instruct the application to fit again. We are happy with the result, so we type n and end the fit.

Results

When it completes, gtlike generates a standard output XML file with the results of the fit.

If you re-run the tool in the same directory, these files will be overwritten by default. Use the clobber=no option on the command line to keep from overwriting the output files.

Unfortunately, the fit details and the value for the -log(likelihood) are not recorded in the automatic output files. You should consider logging the output to a text file for your records by using gtlike ... > fit_data.txt (or something similar) with your gtlike command.

Be aware, however, that this will make it impossible to request a refit when the likelihood process completes.

Example:

gtlike plot=no sfile=3C279output_model.xml > fit_data.txt

In this example, we used the sfile parameter to request that the model results be written to an output XML file, 3C279output_model.xml.

This file contains the source model results that were written to results.dat at the completion of the fit.

NOTE: If you have specified an output XML model file and you wish to modify your model while waiting at the Refit? [y] prompt, you will need to copy the results of the output model file to your input model before making those modifications.

Errors reported with each value in the results.dat file are 1σ estimates (based on the inverse-Hessian at the optimum of the log-likelihood surface).

Other Useful Hidden Parameters

If you are scripting and wish to generate multiple output files without overwriting, the results and specfile parameters allow you to specify output filenames for the results.dat and counts_spectra.fits files respectively.

If you do not specify a source model output file with the sfile parameter, then the input model file will be overwritten with the latest fit. This is convenient as it allows the user to edit that file while the application is waiting at the Refit? [y] prompt so that parameters can be adjusted and set free or fixed. This would be similar to the use of the newpar, freeze, and thaw commands in XSPEC.

8. Make test-statistic maps

Ultimately, one would like to find sources near the detection limit of the instrument. To do this, you model the strongest, most obvious sources (with some theoretical prejudice as to the true source positions, e.g., assuming that most variable high Galactic latitude sources are blazars which can be localized by radio, optical, or X-ray observations), and then create "Test-statistic maps" to search for unmodeled point sources. These TS maps are created by moving a putative point source through a grid of locations on the sky and maximizing -log(likelihood) at each grid point, with the other, stronger, and presumably well-identified sources included in each fit. New, fainter sources are then identified at local maxima of the TS map.

The gttsmap tool can be run with or without an input source model. However, for a useful visualization of the results of the fit, it is recommended you use the output model file from gtlike. The file must be edited so all parameters are fixed (by setting the free attribute to 0 for each parameter) otherwise gttsmap will attempt a refit of the entire model at every point on the grid.

To see the field with the fitted sources removed (i.e. a residuals map), fix all point source parameters before running the TS map. See residuals map wget cell below.

To see the field with the fitted sources included, edit the model to remove all but the diffuse components. See sources map wget cell below.

In both cases, leave the Galactic diffuse prefactor (called Value in the model file) and the isotropic diffuse normalization parameters free during the fit.

Running gttsmap is extremely time-consuming, as the tool is performing a likelihood fit for all events at every pixel position. One way to reduce the time required for this step is to use very coarse binning and/or a very small region. In the following example, we run a TS map for the central 2.5x2.5 degree region of our data file, with .25 degree bins. This results in 100 maximum likelihood calculations. The run time for each of the maps discussed below was a few hours. A 20x20 degree region would require 6400 maximum likelihood calculations and will lengthen the run time significantly.

Here is an example of how to run the gttsmap tool to look for additional sources:

Note: If trying to use CALDB in the example above does not work, just put the IRF name in explicitly. The appropriate IRF for this data set is P8R3_SOURCE_V3.

Because generating TS maps takes a long time, you may wish to download both the residuals and source files in the cells below.

The output from the fit is below.

In the left panel, only diffuse sources were included in the analysis.

The right panel shows the same field, but with the point sources (3C 279 and 3C 273) included in the model, and thus not included in the output image.

This gives a pseudo-residuals map.

The location of the cursor in the image indicates the TS value at each position.

For comparison, we show a 20x20 degree of the same region that requires 6400 maximum likelihood calculations and takes a few days to complete. The bright source clearly seen in the lower part of this image is a recently identified millisecond pulsar.

NOTE: There may be a number of black pixels in these images where the likelihood fit did not converge. You should try adjusting the tolerance or using a different minimizer to keep this from happening.

In this example, the data set was small enough that an unbinned likelihood analysis is possible.

With longer data sets, however, you may run into memory allocation errors, at which time it is necessary to move to binned analysis for your region. A typical memory allocation error looks like:

gtlike(10869) malloc: *** mmap(size=2097152) failed (error code=12)*** error: can't allocate region*** set a breakpoint in malloc_error_break to debug

If you encounter this error, then Binned Likelihood analysis is for you!