Fermi Gamma-ray Space Telescope

Likelihood Analysis with the LATAnalysisScripts

The LATAnalysisScripts are a set of three python libraries to designed to make the analysis of LAT data straightforwared and powerful. They use standard python tools to make their installation and use easy and intuitive. Furthermore, they check for common errors in your analysis, are well documented and expose special feature found in the python tools to the end user (you!).

The three libraries (quickAnalysis, quickLike and quickPlot plus a utilitly library called quickUtils) use a single configuration file and a common naming convention for all of the analysis files. There is detailed, clear error reporting and logging of all results and commands to the screen and to a log file.

To get you going, we are going to walk through the downloading, and installation of the scripts. Then, we'll determine an upper limit on the GeV emission from Swift J164449.3+573451 similar to what was done in Burrows et al. (Nature 476, page 421). Once you've gone through this tutorial you'll be all set to perform your own analysis with the LATAnalysisScripts.

Note: In this tutorial, commands that start with '[user@localhost]$' are to be executed within a shell environment (bash, c-shell etc.) while commands that start with '>>>' are to be executed within an interactive python environment. Commands are in bold and outputs are in normal font. Text in ITALICS within the boxes are comments from the author.

Get and Install the Scripts

Download the LATAnalysisScripts package from the user's contributed area in the FSSC and the pyds9 package from the SAO website. If you haven't already done so, you should download the make2FGLxml script as well from user's area. You'll also need to have ds9 installed which can be found on the same SAO website that pyds9 is on. If you have any issues installing ds9 or pyds9, contact SAO.

Make sure you've set up the Fermi Science Tools. If they aren't set up properly the scripts will not be installed in the correct places. To verify execute:

[user@localhost]$ which gtlike
~/ScienceTools/i386-apple-darwin10.8.0/bin/gtlike

If the which command doesn't return anything, then the tools aren't set up. Verify that they are before continuing.

Unpack and install pyds9:

[user@localhost]$ tar zxvf pyds9-1.3.tar.gz
[user@localhost]$ cd pyds9-1.3
[user@localhost]$ python setup.py install

Install make2FGLxml.py into the ScienceTools python tree:

[user@localhost]$ cp make2FGLxml.py $FERMI_DIR/lib/python2.6/site-packages/

Unpack and install the LATAnalysisScripts:

[user@localhost]$ tar zxvf LATAnalysisScripts-0.1.8.tar.gz
[user@localhost]$ cd LATAnalysisScripts-0.1.8
[user@localhost]$ python setup.py install

Verify that everything is installed and working:

[user@localhost]$ python
Python 2.6.5 (r265:79063, Nov 9 2010, 12:49:33)
[GCC 4.2.1 (Apple Inc. build 5664)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>>import make2FGLxml
This is make2FGLxml version 04r1.
For use with the gll_psc_v02.fit and gll_psc_v05.fit and later LAT catalog files.
>>>import ds9
>>>import quickAnalysis
>>>exit()

If any of these commands print errors, then something is wrong. Verify that the Fermi ScienceTools are correctly installed and set up and verify that the external tools are installed correctly.

What do I have now?

In addition to the pyds9 and make2FGLxml packages, you have the LATAnalysisScripts which are made up of three libraries and their corresponding command line tools: quickAnalysis, quickLike and quickPlot. Most of what is explained here can be found in the README file in the LATAnalysisScripts package you downloaded. You can also find information on the specific tools by executing the command line tool without any options or by using the 'help' command in python. For example, try:

[user@localhost]$ quickAnalysis
...help text suppressed

and you should see the help text for the command line version of quickAnalysis. Details on all of the functions would be available if you started python, imported the module (using the 'import quickAnalysis' command) and then typing 'help(quickAnalysis)'. Now that you know where to find all of the documentation, you should be good to go, but let's do a quick example analysis to get you started.

Get the Data

For this thread the original data were extracted from the LAT data server with the following selections (these selections are similar to those in the paper):

  • Search Center (RA,Dec) = (251.2054,+57.5808)
  • Radius = 30 degrees
  • Start Time (MET) = 322963202 seconds (2011-03-28T00:00:00)
  • Stop Time (MET) = 323568002 seconds (2011-04-04T00:00:00)
  • Minimum Energy = 100 MeV
  • Maximum Energy = 300000 MeV

We've provided direct links to the event files as well as the spacecraft data file if you don't want to take the time to use the download server. For more information on how to download LAT data please see the Extract LAT Data tutorial.

Initialize the Analysis

The LATAnalysisScripts read everything from a config file called basename.cfg where basename is some prefix of your choosing (all of the files created by the LATAnalysisScripts will be prepended with this prefix). In our case, we'll use 'SwiftJ1644' as our basename. So, first, create a default config file and rename it to SwiftJ1644.cfg:

[user@localhost]$ quickAnalysis -i
Creating example configuration file called example.cfg
2012-02-01 11:56:00,717 - quickAnalysis - INFO - Created quickAnalysis object: irfs=P7SOURCE_V6, verbosity=0, binned=False, base=example, eventclass=2, rad=10, tmin=INDEF, ra=0, emin=100, emax=300000, tmax=INDEF, dec=0, zmax=100, binsize=0.1,
2012-02-01 11:56:00,718 - quickAnalysis - INFO - wrote common config to example.cfg.
2012-02-01 11:56:00,718 - quickAnalysis - INFO - wrote quickAnalysis config to example.cfg.
[user@localhost]$ cp example.cfg SwiftJ1644.cfg

Now, open that config file up in your favorite text editor and edit up the values to match the analysis we're doing here. The final file should look like:

[common]
irfs = P7SOURCE_V6
verbosity = 0
binned = False
base = SwiftJ1644
eventclass = 2

[quickAnalysis]
rad = 10
tmin = 322963202
ra = 251.2054
emin = 100
emax = 300000
tmax = 323568002
dec = 57.5808
zmax = 100
binsize = 0.1

Note that there is a 'common' section that is used in all three of the LATAnalysisScripts module and there is a section for quickAnalysis which is what we are about to run. When we want to use quickLike and quickPlot we'll have to execute them with the '-i' option as well and create the relavant sections for them in this config file. Also note that there's a file in our working directory called 'example_quickAnalysis.log which shows what we just did. From here on out, the scripts will log everything we do in files called SwiftJ1644_quickAnalysis.log, SwiftJ1644_quickLike.log and SwiftJ1644_quickPlot.log.

There are some other things we need to do before we start. We need to create a file that lists all of the photon files called SwiftJ1644.list, rename the spacecraft file SwiftJ1644_SC.fits.

[user@localhost]$ ls -1 *PH* > SwiftJ1644.list
[user@localhost]$ ln -s L120404155224B0489E7F72_SC00.fits SwiftJ1644_SC.fits

quickAnalysis

We're ready to go, just execute quickAnalysis with the -a option (to tell it to run an analysis) and the -n option follwed by the basename so that it knows which config file to load:

[user@localhost]$ quickAnalysis -a -n SwiftJ1644
...output suppressed.

This should take about 10 minutes and print out everyting it's doing to the screen and into the log file. Once you're done, you should have a filtered eventfile that's been run through gtmktime, a livetime cube and an exposure map. Everything you need for a likelihood anlaysis except for a model file. First, download the 2FGL point source catalog fits file from the FSSC website (or copy it there), make links to the diffuse models and execute quickAnalysis with the '-x' switch.

[user@localhost]$ ln -s $FERMI_DIR/refdata/fermi/diffuseModels/gal_2yearp7v6_v0.fits
[user@localhost]$ ln -s $FERMI_DIR/refdata/fermi/diffuseModels/iso_p7v6source.txt
[user@localhost]$ quickAnalysis -x -n SwiftJ1644
Creating XML model file from 2FGL
2012-02-01 13:48:49,305 - quickAnalysis - INFO - Reading from config file (SwiftJ1644.cfg)
2012-02-01 13:48:49,306 - quickAnalysis - INFO - Reading common variables...
2012-02-01 13:48:49,306 - quickAnalysis - INFO - Reading quickAnalysis variables...
2012-02-01 13:48:49,307 - quickAnalysis - INFO - Created quickAnalysis object: irfs=P7SOURCE_V6, verbosity=0, base=SwiftJ1644, binned=False, eventclass=2, rad=20, tmin=322963202, emin=100, emax=300000, tmax=323568002, binsize=0.1, ra=251.2054, zmax=100, dec=57.5808,
2012-02-01 13:48:49,307 - quickAnalysis - CRITICAL - SwiftJ1644_model.xml doesn't exist.
2012-02-01 13:48:49,307 - quickAnalysis - INFO - SwiftJ1644_model.xml doesn't exist, will create a new one.
This is make2FGLxml version 04r1.
For use with the gll_psc_v02.fit and gll_psc_v05.fit and later LAT catalog files.
Creating file and adding sources for 2FGL
Added 33 point sources and 0 extended sources
2012-02-01 13:48:50,298 - quickAnalysis - INFO - NOTE: if there are extended sources in your ROI, make sure the correspoinding diffuse template is in the working directory.

This model file doesn't include the source we are interested in so we should add the following lines to the top of the SwiftJ1644_model.xml model file.

<source name="SwiftJ1644" type="PointSource">
 <spectrum type="PowerLaw2">
  <parameter free="true" max="10000.0" min="0.0001" name="Integral" scale="1e-07" value="1.0"/>
  <parameter free="true" max="5.0" min="0.0" name="Index" scale="-1.0" value="2.0"/>
  <parameter free="false" max="500000.0" min="20.0" name="LowerLimit" scale="1.0" value="100.0"/>
  <parameter free="false" max="500000.0" min="20.0" name="UpperLimit" scale="1.0" value="300000.0"/>
 </spectrum>
 <spatialModel type="SkyDirFunction">
  <parameter free="false" max="360.0" min="-360.0" name="RA" scale="1.0" value="251.2054"/>
  <parameter free="false" max="90.0" min="-90.0" name="DEC" scale="1.0" value="57.5808"/>
 </spatialModel>
</source>

quickLike

Before we do a likelihood analysis with quickLike we need to add the quickLike section to the config file. Execute quickLike with the -i option which will create another example.cfg file. Copy the quickLike section into the SwiftJ1644.cfg file and edit it to match our anaylsis here. The final quickLike section should look like:

[quickLike]
sourcename = SwiftJ1644
model = SwiftJ1644_model.xml
drmtol = 0.1
mintol = 0.0001

Start up python, import the quickLike module and create your quickLike object by passing it the basename and 'True'. The 'True' means that you want the module to read from the config file. There are other options you can pass to the module which are detailed in the help documentation (accessable via 'help(quickLike)'). Usually we would initialize the DRM routine by running the initDRM command and then the fitDRM command but we are going to skip this and go directly to using the NewMinuit algorithm. We'll need to tell the initMIN command to use the model file we created instead of the output of the DRM fit (which doesn't exist).

[user@localhost]$ python
Python 2.6.5 (r265:79063, Nov 9 2010, 12:49:33)
[GCC 4.2.1 (Apple Inc. build 5664)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>>from quickLike import *
>>>qL = quickLike('SwiftJ1644', True)
...output supressed
>>>qL.makeObs()
...output supressed
>>>qL.initMIN(modelFile="SwiftJ1644_model.xml")
...output supressed

We're now going to do the fit, but first, we want to remove the SwiftJ1644 source from our model and save it for use later using the unloadSource function. This way, we can get a better convergance since the Swift source doesn't really have gamma-ray emission and would be hard to fit (another way to do this, is to set all of the parameters for SwiftJ1644 to '0' and fix them in the XML file). The fit will take about five minutes to finish.

>>>qL.unLoadSource('SwiftJ1644')
2012-02-02 12:50:51,589 - quickLike - INFO - Removed SwiftJ1644 from the model and saved it.
>>>qL.fitMIN()
...some output suppressed
2012-02-02 12:55:08,062 - quickLike - INFO - NEWMINUIT Fit Finished. Total TS: -11469.3411544
2012-02-02 12:55:08,063 - quickLike - INFO - NEWMINUIT Fit Status: 156
2012-02-02 12:55:08,063 - quickLike - INFO - NEWMINUIT fit Distance: 1.35106292821e-05
2012-02-02 12:55:08,063 - quickLike - ERROR - NEWMINUIT DID NOT CONVERGE!!!
2012-02-02 12:55:08,063 - quickLike - ERROR - The fit failed the following tests: HasMadePosDefCovar HasPosDefCovar HasAccurateCovar
...some output suppressed

Now, if you've done things similar to what has been done here, you should see the lines above about the fit not converging. The fun now begins where we try and figure out why the fit wasn't converging and what we can do about it. If you execute the 'paramsAtLimit' function you can see that some parameters in the model are at their limits.

>>>qL.paramsAtLimit()
2012-02-02 15:51:32,050 - quickLike - ERROR - The Prefactor (0.000100002150128) of _2FGLJ1531.0+5725 is close (2.1501276315e-05) to its lower limit (0.0001)
2012-02-02 15:51:32,062 - quickLike - ERROR - The Index (4.99999994135) of _2FGLJ1604.6+5710 is close (1.17301164337e-08) to its upper limit (5.0)
2012-02-02 15:51:32,063 - quickLike - ERROR - The Index (4.92274288596) of _2FGLJ1656.5+6012 is close (0.0154514228086) to its upper limit (5.0)
2012-02-02 15:51:32,065 - quickLike - ERROR - The Index (4.98785968966) of _2FGLJ1559.0+5627 is close (0.00242806206775) to its upper limit (5.0)

Since our model is based on the 2FGL which was produced using two years of data and we are only dealing with a few days of data here, there are probably a lot of low TS sources in our ROI. Let's find them and then remove some of them from the model

>>>qL.reLoadSource()
2012-02-02 15:52:42,865 - quickLike - INFO - Reloaded saved source.
>>>qL.removeWeak()
...output supressed

The output above should show that there are a lot of really low TS sources. Let's remove all of the Free sources that have a TS lower than 1 and then try the fit again.

>>>qL.removeWeak(tslimit=1.0, RemoveFree=True)
...output supressed
>>>qL.unLoadSource('SwiftJ1644')
2012-02-02 12:53:51,589 - quickLike - INFO - Removed SwiftJ1644 from the model and saved it.
>>>qL.fitMIN()
...some output suppressed
2012-02-02 12:55:08,062 - quickLike - INFO - NEWMINUIT Fit Finished. Total TS: -11470.3047291
2012-02-02 12:55:08,063 - quickLike - INFO - NEWMINUIT Fit Status: 156
2012-02-02 12:55:08,063 - quickLike - INFO - NEWMINUIT fit Distance: 0.000256727730312
2012-02-02 12:55:08,063 - quickLike - ERROR - NEWMINUIT DID NOT CONVERGE!!!
2012-02-02 12:55:08,063 - quickLike - ERROR - The fit failed the following tests: HasMadePosDefCovar HasPosDe\ fCovar HasAccurateCovar

Looks like it didn't fit again. This time, let's remove all the sources that have a negative TS even if they're fixed.

>>>qL.reLoadSource()
>>>qL.removeWeak(tslimit=0.0, RemoveFree=True, RemoveFixed=True)
>>>qL.unLoadSource('SwiftJ1644')
2012-02-02 17:26:51,589 - quickLike - INFO - Removed SwiftJ1644 from the model and saved it.
>>>qL.fitMIN()
2012-02-02 17:27:01,039 - quickLike - INFO - NEWMINUIT Fit Finished. Total TS: -11470.0702022
2012-02-02 17:27:01,039 - quickLike - INFO - NEWMINUIT Fit Status: 0
2012-02-02 17:27:01,039 - quickLike - INFO - NEWMINUIT fit Distance: 3.52763544413e-06

Looks like that did it. However, when we check if any parameters are at their limits:

>>>qL.paramsAtLimit()
2012-02-02 17:29:45,456 - quickLike - ERROR - The Index (4.99999570933) of _2FGLJ1559.0+5627 is close (8.58133958737e-07) to its upper limit (5.0)

You can see that this source is not being fit well. It's best to reset this source back to it's catalog values, fix the index (this is the parameter that is having issues) and try the fit again. Exit out of python, copy the output of the fitting to a new model file, edit this source to it's catalog value and try the fitting again.

>>>exit()
[user@localhost]$ cp SwiftJ1644_likeMinuit.xml SwiftJ1644_minimal.xml

Now edit SwiftJ1644_minimal.xml by inserting the xml for the SwiftJ1644 source back into it (like we did at the beginning) and replace the paramters of _2FGLJ1559.0+5627 with the parameters found in SwiftJ1644_model.xml. Set the index of this source from free to fixed. Then start python back up.

[user@localhost]$ python
Python 2.6.5 (r265:79063, Nov 9 2010, 12:49:33)
[GCC 4.2.1 (Apple Inc. build 5664)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>>from quickLike import *
>>>qL = quickLike('SwiftJ1644', True)
...output supressed
>>>qL.makeObs()
...output supressed
>>>qL.initMIN(modelFile="SwiftJ1644_minimal.xml")
...output supressed
>>>qL.unLoadSource('SwiftJ1644')
2012-02-02 17:38:14,752 - quickLike - INFO - Removed SwiftJ1644 from the model and saved it.
>>>qL.fitMIN()
...output supressed
2012-02-02 17:39:57,667 - quickLike - INFO - NEWMINUIT Fit Finished. Total TS: -11470.7003776
2012-02-02 17:39:57,668 - quickLike - INFO - NEWMINUIT Fit Status: 0
2012-02-02 17:39:57,668 - quickLike - INFO - NEWMINUIT fit Distance: 3.4341373915e-05
>>>qL.paramsAtLimit()

That seems to have done it. Now, the only thing to do is to add the Swift source back in and calculate the upper limit (the paper used an upper energy limit of 10GeV so that's what we're doing here).

>>>qL.reLoadSource()
2012-02-09 17:41:08,647 - quickLike - INFO - Reloaded saved source.
>>>qL.calcUpper('SwiftJ1644', Emax=10000)
...output supressed
2012-02-02 17:45:50,554 - quickLike - INFO - SwiftJ1644 UL: 3.48e-07 ph/cm^2/s for emin=100.0, emax=10000.0, delta(logLike)=1.35

Note that this is in ph/cm^2/s and not ergs/cm^2/s.

Finished

So, this thread has covered the basics of what you can do with the LATAnalysisScripts. The goal of these scripts is to make the analysis of LAT data more user friendly so if you have suggestions or feature requests, please contact the helpdesk. There are a lot of features we haven't covered here but all of the documentation is at your fingertips. You can access the help within python, using pydoc or you can get them here:

You can also help out! The version of the code in the user's contributed area is a stable, release version. The bleeding edge version of the code is always located at GitHub specifically in the developer's space: LATAnalysisScripts on GitHub. You can download the most recent version there or you can fork it yourself, modify it in any way and even push changes back to the master. We would love your help! Note that the version on GitHub might not always work and might include bugs. Use at your own risk.

Last updated by: Jeremy S. Perkins 04/04/2012