In this hands-on activity, you will learn how to produce get the parameters of a spectral model—let’s say, flux and spectral index—for the blazar PG1553+113 observed with Fermi LAT. You will perform a full likelihood analysis for PG1553+113.
The results from the analysis of this dataset were published by the Fermi LAT Collaboration. So you will be able to compare the results of your analysis with those published in that paper.
A number of steps are necessary to fit a source's spectra; these are described in detail below.
- Select the data. The data from a reasonably large spatial region around the source(s) being analyzed must be used because of the overlapping of the point spread functions (PSF) of nearby sources.
- Select the model. This 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.
- Precompute a number of quantities that are part of the likelihood computation. As the parameter values are varied in searching for the best fit, the likelihood is calculated many times. While not strictly necessary, precomputing a number of computation-intensive quantities will greatly speed up the fitting process.
- Finally, perform the fit. The parameter space can be quite large—the spectral parameters from a number of sources must be fit simultaneously—and therefore the likelihood tools provide a choice of three numerical “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 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.
We will use the Enrico python tools instead of using directly the ScienceTools. Enrico is very convenient because it hides away the complexity of the analysis that is going on under the hood. Enrico takes care of calling all the appropriate ScienceTools in the correct order.
First of all, let’s test Enrico’s installation. Open a terminal and type:
enrico_setupcheck
Everything should go well because we are running the activity in the VM, where Enrico is pre-installed. Now that you’ve checked that everything is in place, let’s get the data you will analyze.
Let’s extract the events for a region of radius 15˚ around PG 1553 in the time range 2008-08-05 to 2009-02-21. These are the parameters that you will use for the exercise:
Normally you would go to go to the NASA FSSC server website, input the selection criteria and download the files. However, in order to save time and bandwidth, all the necessary files have already been downloaded and are available in the VM in the folder LAT_day02/data. The data files are the following:
pg1553_PH0#.fits: events filespg1553_SC00.fits: spacecraft file
In the terminal, now type
cd ~/LAT_day02
mkdir flux
cd fluxYou are now in the directory which will hold the results of the likelihood analysis.
Since there is more than one events file, you will need to provide an input file list in your analysis below. This text file can be generated by typing:
ls ~/LAT_day02/data/pg1553_PH* > events.txtEnrico uses configuration files to run the analysis. You can use the enrico_config tool to quickly make a config file called pg1553.conf. It will ask you for the required options and copy the rest from a default config file located at enrico/data/config/default.conf:
Create configuration file
[fermi@localhost flux]$ enrico_config pg1553.conf
[Message]: Please provide the following required options [default] :
Output directory [/home/fermi/LAT_day02/flux] :
Target Name : PG1553
Right Ascension: 238.93
Declination: 11.19
redshift, no effect if null [0] : 0.36
ebl model to used
0=Kneiske, 1=Primack05, 2=Kneiske_HighUV, 3=Stecker05, 4=Franceschini, 5=Finke, 6=Gilmore : 0
Options are : PowerLaw, PowerLaw2, LogParabola, PLExpCutoff, Generic
Generic is design to allow the user to fit with non-supported models
EBL absorption can be added for PowerLaw2, LogParabola, PLExpCutoff
Spectral Model [PowerLaw] : PowerLaw2
ROI Size [15] :
FT2 file [/home/fermi/enrico/Data/download/lat_spacecraft_merged.fits] : /home/fermi/LAT_day02/data/pg1553_SC00.fits
FT1 list of files [/home/fermi/enrico/Data/download/events.lis] : events.txt
tag [LAT_Analysis] : spectrum
Start time [239557418] :
End time [334165418] :
Emin [100] :
Emax [300000] :
IRFs [CALDB] :
evclass [128] :
evtype [3] :
Corresponding IRFs = ('P8R2_SOURCE_V6', ['BACK', 'FRONT'])
Is this ok? [y] :
Corresponding zmax = 90
Some notes:
- Always give the full path for the files in enrico
- We are using the
PowerLaw2model--a simple power-law spectrum with the form
dN/dE \propto E^gamma
where dN/dE is the monochromatic photon flux, gamma is the spectral index.
- Time is give in MET
- Energy is given in MeV
- ROI size is given in degrees
- The redshift of this blazar is z=0.36. This turns out to be important, as we will see in the next activity.
If you want, you can edit this config file by hand to make further adjustments.
The ScienceTools need a source model written in the xml format. You can run enrico_xml to make such model of the sky and store it into a xml file which will be used for the analysis. The options for this step are directly provided in the config file. For the enrico_xml tool, the relevant options are in the [space] and [target] sections. The out file is given by [file]/xml field.
This tool will automatically add the following sources to the xml source model file:
- Your target source
- The galactic (
GalDiffModel) and isotropic (IsoDiffModel) diffuse components that are the dominant background sources in most LAT analysis. - All the LAT sources from the 4-year catalog (3FGL) that are inside the ROI. The spectral parameters of the sources within 3 degrees of our source are left free so they can be fit simultaneously with our source, whereas those further away are fixed to their catalog values.
Create sources model
[fermi@localhost flux]$ enrico_xml pg1553.conf
use the default location of the catalog
use the default catalog
Use the catalog : /home/fermi/enrico/Data/catalog/gll_psc_v16.fit
[Message]: Summary of the XML model generation
Add 41 sources in the ROI of 17.0 ( 15.0 + 2 ) degrees
4 sources have free parameters inside 3.0 degrees
0 source(s) is (are) extended
Iso model file /home/fermi/enrico/Data/diffuse/iso_P8R2_SOURCE_V6_v06.txt
Galactic model file /home/fermi/enrico/Data/diffuse/gll_iem_v06.fits
[Message]: write the Xml file in /home/fermi/LAT_day02/flux/PG1553_PowerLaw2_model.xml
You can explore the PG1553_PowerLaw2_model.xml output file with a text editor.
This is where the magic happens. The enrico_sed tool will take care of executing all the steps necessary for you to get a global likelihood fit and find the best-fit parameters by minimizing the likelihood function. enrico_sed will look into the options you have selected in pg1553.conf and call the following ScienceTool tasks (here is an overview of the LAT analysis tools):
gtselect: Perform event selection (ROI, energy cuts)gtmktime: Perform time selection based on spacecraft file, getting GTIsgtbin: Bin the data, computes counts mapgtltcube: Perform the calculation of the livetime cube. This is the most computationally intensive step, taking.gtexpcube2: Use the previously generated livetime cube and apply it to your ROI to obtain an exposure cube.gtsrcmaps: Create model counts maps for each of the sources in the source model catalog. This is used to speed up the likelihood calculation ofgtlike.gtlike: Maximize the likelihood function numerically and get initial estimate of parameters- Optimize fit: improve parameter estimate
To run the likelihood maximization just issue in the terminal:
enrico_sed pg1553.conf
You will need to be patient because this will take quite a while. On my 2015 MacBook running the analysis in the VM, this took about 50 minutes to complete. Some possibilities of things to do while you wait:
- go get some water or coffee
- google for radiative processes in blazars
- come see our computer cluster at the second floor of IAG
- discuss gamma-ray astrophysics with your colleagues
After the analysis finishes, you should get something similar to this output from enrico_sed.
Inspect the log generated on the terminal by enrico_sed and identify the different tasks that are being performed.
Open the file PG1553_PowerLaw2_239557418_334165418_100_300000_spectrum.results—which holds the results of the likelihood fit—in a text editor and answer the following questions:
- What is the resulting spectral index and uncertainty?
- What is the flux? What do you think are the units?
- What is the approximate statistical significance of the presence of a point source in the location of this blazar, in standard deviations? Is the source statistically significant against background?
The observed counts map is available in the file PG1553_spectrum_CountMap.fits and the model counts map is in PG1553_spectrum_ModelMap.fits. The residual counts map resulting from the subtraction of the model map and the observed map is in PG1553_Residual_Model_cmap.fits.
Open and inspect these three FITS files with ds9 or IPython. Do you think your model accounts well for the observed photons? Why?
You can use enrico_testmodel to compute the log(likelihood) of the PowerLaw, LogParabola and PLExpCutoff models. An ascii file is then produced in the Spectrum folder with the value of the log(likelihood) for each model. Use the Wilk’s theorem to decide which model best describes the data.