{ "cells": [ { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "# Likelihood Analysis of a bright point source with fermipy" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "In this tutorial, we will demonstrate how to perform a standard point-source analysis of Fermi-LAT data using fermipy. \n", "\n", "This sample analysis is based on the PG 1553+113 analysis performed by the LAT team and described in [Abdo, A. A. et al. 2010, ApJ, 708, 1310](http://adsabs.harvard.edu/abs/2010ApJ...708.1310A) and closely follows the [Likelihood Analysis with Python](http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/python_tutorial.html) thread. This tutorial assumes you have the most recent ScienceTools installed and [fermipy](http://fermipy.readthedocs.org) installed on top of it. For instructions on installing fermipy and the Fermi ScienceTools you should consult the [fermipy Installation Instructions](http://fermipy.readthedocs.org/en/latest/install.html). We will also make significant use of python, so you might want to familiarize yourself with python including matplotlib and other libraries (there's a beginner's guide at http://wiki.python.org/moin/BeginnersGuide). " ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "## Get the Data" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "For this thread the original data were extracted from the [LAT data server](http://fermi.gsfc.nasa.gov/cgi-bin/ssc/LAT/LATDataQuery.cgi) with the following selections (these selections are similar to those in the paper):\n", "\n", "* Search Center (RA,Dec) = (238.929,11.1901)\n", "* Radius = 30 degrees\n", "* Start Time (MET) = 239557417 seconds (2008-08-04T15:43:37)\n", "* Stop Time (MET) = 256970880 seconds (2009-02-22T04:48:00)\n", "* Minimum Energy = 100 MeV\n", "* Maximum Energy = 300000 MeV\n", "\n", "Once you exectute the query you can download the data and put it in your working directory. You'll need to change the names of the files to match the filenames that the data server gave you. Alternatively you can run with the following tarball which already contains the downloaded files as well as all of the ancillary files that will be generated for the analysis. In this example the output of the analysis will go into a subdirectory called *pg1553*." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "import os\n", "if os.path.isfile('../data/pg1553.tar.gz'):\n", " !tar xzf ../data/pg1553.tar.gz\n", "else:\n", " !curl -OL https://raw.githubusercontent.com/fermiPy/fermipy-extras/master/data/pg1553.tar.gz\n", " !tar xzf pg1553.tar.gz" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "If you would like to download the raw data files used in this analysis (FT1 and FT2) you can do so by executing the following cell." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "if os.path.isdir('../data'):\n", " os.environ['DATADIR'] = '../data'\n", "else:\n", " !mkdir pg1553/data\n", " !mkdir pg1553/data/ft1\n", " !mkdir pg1553/data/ft2\n", " os.environ['DATADIR'] = 'pg1553/data'\n", " !curl -o pg1553/data/ft1/L1504241622054B65347F25_PH00.fits -OL https://raw.githubusercontent.com/fermiPy/fermipy-extras/master/data/ft1/L1504241622054B65347F25_PH00.fits\n", " !curl -o pg1553/data/ft1/L1504241622054B65347F25_PH01.fits -OL https://raw.githubusercontent.com/fermiPy/fermipy-extras/master/data/ft1/L1504241622054B65347F25_PH01.fits\n", "\n", " \n", "#!curl -o pg1553/L1504241622054B65347F25_PH00.fits http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/pyLikelihood/L1504241622054B65347F25_PH00.fits\n", "#!curl -o pg1553/L1504241622054B65347F25_PH01.fits http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/pyLikelihood/L1504241622054B65347F25_PH01.fits\n", "#!curl -o pg1553/L1504241622054B65347F25_SC00.fits http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/pyLikelihood/L1504241622054B65347F25_SC00.fits" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "### Make a file list\n", "\n", "You'll then need to make a file list with the names of your input event files. You can either just make one with a text editor or do the following from the command line." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [], "source": [ "!ls -1 $DATADIR/ft1/*PH*.fits > pg1553/ft1.txt" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "../data/ft1/L1504241622054B65347F25_PH00.fits\n", "../data/ft1/L1504241622054B65347F25_PH01.fits\n" ] } ], "source": [ "!cat pg1553/ft1.txt" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "## Make a config file\n", "\n", "fermipy bases its analysis on a configuration file (in [yaml](http://yaml.org) format). We're just going to use a really simple config file for a standard analysis. There are many many more options which you can use or you can modify these options after the fact within the analysis chain.\n", "\n", "\n", "Make a config file named 'config.yaml' like the following. For more details on the config file see [config.html](http://fermipy.readthedocs.org/en/latest/config.html). You will probably need to customize this a bit since your files might not be in the same place or named the same. The galactic and isotropic diffuse will need to be located on your system (they are included in the science tools or can be downloaded from the FSSC). In the following example we set the path to these files with the environment variable FERMI_DIFFUSE_DIR. If FERMI_DIFFUSE_DIR is not defined fermipy will look for the location of these files within the FSSC STs distribution. \n", "\n", "```\n", "data:\n", " evfile : ft1.txt \n", " scfile : ft2.tx \n", "\n", "binning:\n", " roiwidth : 10.0 \n", " binsz : 0.1\n", " binsperdec : 8\n", "\n", "selection :\n", " emin : 100\n", " emax : 300000\n", " tmin: 239557417 \n", " tmax: 256970880\n", " zmax : 90\n", " evclass : 128\n", " evtype : 3\n", " target : '4FGL J1555.7+1111'\n", "\n", "gtlike:\n", " edisp : True\n", " irfs : 'P8R3_SOURCE_V3'\n", "\n", "model:\n", " src_roiwidth : 15.0\n", " galdiff : '$FERMI_DIFFUSE_DIR/gll_iem_v07.fits'\n", " isodiff : '$FERMI_DIFFUSE_DIR/iso_P8R3_SOURCE_V3_v1.txt'\n", " catalogs :\n", " - '4FGL'\n", "```" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "## Start the analysis\n", "\n", "Next, you create an analysis script and run the setup steps which include running the selections and generating exposure maps etc. This will take a bit.\n", "\n", "This is where the magic happens. fermipy will load the point source model, create an xml file for you which contains the models for all your sources in the region of interest (ROI), decide on all the appropriate cuts and binnings and just go. All of this is configurable from python or from the config file. And, if you need to rerun things, it's smart enough to not overwrite files if it doesn't need to." ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "### Load up some useful modules" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true, "deletable": true, "editable": true }, "outputs": [], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "import matplotlib\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's ignore some deprecation warnings" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "import warnings\n", "from matplotlib import MatplotlibDeprecationWarning\n", "warnings.filterwarnings(\"ignore\", category=MatplotlibDeprecationWarning) " ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "### Import the GTAnalysis module from fermipy\n", "\n", "You start by importing the module and then creating an instance of the analysis object from our config file. When instantiating the analysis object we can override any options defined in the configuration file by passing keyword arguments to the object constructor. Here we explicitly set the verbosity parameter to 3 (INFO) which supresses DEBUG output. When we create the object, it spits out a bunch of information about all of the parameters that were used. You can see there are many more options than the ones we chose." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "WARNING: AstropyDeprecationWarning: astropy.extern.six will be removed in 4.0, use the six module directly if it is still needed [astropy.extern.six]\n", "2021-06-23 09:58:36 INFO GTAnalysis.__init__(): \n", "--------------------------------------------------------------------------------\n", "fermipy version v1.0.1 \n", "ScienceTools version 2.0.18\n" ] } ], "source": [ "from fermipy.gtanalysis import GTAnalysis\n", "gta = GTAnalysis('pg1553/config.yaml',logging={'verbosity': 3})\n", "matplotlib.interactive(True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's take a look at the initial input event files" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['../data/ft1/L1504241622054B65347F25_PH00.fits', '../data/ft1/L1504241622054B65347F25_PH01.fits']\n" ] } ], "source": [ "with open(\"pg1553/ft1.txt\") as f:\n", " input_files = f.readlines()\n", "\n", "input_files = [f.strip(\"\\n\") for f in input_files]\n", "print(input_files)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Read one file in as a table" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "from astropy.table import Table" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "t = Table.read(input_files[0], hdu=1)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Table length=113339\n", "
ENERGY | RA | DEC | L | B | THETA | PHI | ZENITH_ANGLE | EARTH_AZIMUTH_ANGLE | TIME | EVENT_ID | RUN_ID | RECON_VERSION | CALIB_VERSION [3] | EVENT_CLASS [32] | EVENT_TYPE [32] | CONVERSION_TYPE | LIVETIME | DIFRSP0 | DIFRSP1 | DIFRSP2 | DIFRSP3 | DIFRSP4 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
MeV | deg | deg | deg | deg | deg | deg | deg | deg | s | s | ||||||||||||
float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float64 | int32 | int32 | int16 | int16 | bool | bool | int16 | float64 | float32 | float32 | float32 | float32 | float32 |
538.79675 | 258.76797 | 11.601158 | 32.717495 | 26.587275 | 21.099522 | 69.525536 | 39.699387 | 302.11157 | 239646594.1632708 | 3154345 | 239645594 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 130.803790807724 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
765.5338 | 258.56604 | 11.574936 | 32.59832 | 26.756044 | 57.46366 | 200.40448 | 70.16773 | 66.16222 | 239840161.56974736 | 8362449 | 239835659 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 76.85447326302528 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
391.90915 | 258.77057 | 11.551644 | 32.668377 | 26.564293 | 66.689735 | 252.94148 | 38.697117 | 33.092915 | 239898124.18749014 | 2213352 | 239897571 | 0 | 0 .. 0 | False .. True | False .. False | 1 | 56.99404388666153 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
120.571365 | 258.83325 | 11.796055 | 32.94541 | 26.610472 | 44.83362 | 96.635826 | 54.28294 | 276.6544 | 239951043.89974436 | 1003061 | 239950625 | 0 | 0 .. 0 | False .. True | False .. False | 1 | 123.81579771637917 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
500.81796 | 258.63654 | 11.483924 | 32.538124 | 26.655235 | 48.03507 | 98.21642 | 57.537155 | 275.82956 | 239951099.17692286 | 1148299 | 239950625 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 179.09297621250153 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
393.2541 | 259.0739 | 11.939454 | 33.20065 | 26.455925 | 32.59749 | 194.16124 | 39.280994 | 38.334652 | 240030010.17818636 | 10147729 | 240025172 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 24.26987051963806 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
260.86182 | 258.9468 | 12.058491 | 33.264046 | 26.618439 | 61.366016 | 281.5017 | 31.956478 | 287.96527 | 240036675.37049332 | 336220 | 240036563 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 109.3906192779541 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
155.16846 | 258.91968 | 11.72253 | 32.91009 | 26.502935 | 55.22312 | 197.95798 | 61.02289 | 60.150562 | 240138580.18634808 | 13479863 | 240133960 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 78.8684054017067 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
159.39615 | 259.2734 | 11.618165 | 32.96556 | 26.144552 | 22.752687 | 165.9323 | 19.66616 | 13.963443 | 240202427.6246177 | 12010846 | 240196811 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 222.24398529529572 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
133.3081 | 224.78586 | -1.6613442 | 354.95584 | 47.886524 | 58.574303 | 341.0126 | 61.304764 | 282.9876 | 248159241.57002193 | 10403360 | 248154894 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 30.230947762727737 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
314.87418 | 224.90024 | -1.6593533 | 355.0792 | 47.80763 | 50.744946 | 18.809458 | 21.62766 | 77.092995 | 248163738.1249223 | 7501457 | 248160623 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 301.02908006310463 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
145.44543 | 224.17506 | -1.7963057 | 354.16092 | 48.216053 | 48.30415 | 11.835739 | 15.309311 | 345.5459 | 248164144.1878574 | 8388840 | 248160623 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 121.36332374811172 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
2333.3445 | 225.02441 | -1.1435628 | 355.75418 | 48.08519 | 51.67315 | 5.2675304 | 22.999802 | 323.16785 | 248164312.956939 | 8864099 | 248160623 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 87.07305860519409 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
447.18857 | 225.62411 | -2.9184122 | 354.5411 | 46.40254 | 42.133 | 6.514918 | 40.86052 | 96.61587 | 248169145.10990927 | 7693886 | 248166353 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 90.15713900327682 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
459.95984 | 226.6249 | -4.038852 | 354.44446 | 44.902985 | 34.743168 | 5.922209 | 32.011417 | 90.60807 | 248215208.9780431 | 6094041 | 248211903 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 178.89627787470818 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
291.9761 | 225.77574 | -1.9272188 | 355.71527 | 47.00006 | 23.413742 | 358.9999 | 19.086782 | 63.520535 | 248215454.9093352 | 6544969 | 248211903 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 196.8822024166584 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
722.9538 | 224.60498 | -0.7475925 | 355.73172 | 48.661766 | 34.350864 | 332.4631 | 34.535465 | 307.61978 | 248216130.13690475 | 7908530 | 248211903 | 0 | 0 .. 0 | False .. True | False .. False | 1 | 39.67203438282013 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
274.50766 | 225.04887 | -1.8633881 | 355.02316 | 47.55807 | 45.89875 | 337.548 | 23.772915 | 323.13287 | 248227417.63520315 | 6824508 | 248223836 | 0 | 0 .. 0 | False .. True | False .. True | 0 | 88.35028231143951 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
144.3101 | 224.8169 | -0.6086211 | 356.10675 | 48.608936 | 51.036015 | 337.5925 | 54.70884 | 288.92017 | 248227951.11360538 | 8201273 | 248223836 | 0 | 0 .. 0 | False .. True | False .. False | 1 | 53.631305277347565 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |