{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "
\n", "**This is a fixed-text formatted version of a Jupyter notebook.**\n", "\n", " You can contribute with your own notebooks in this\n", " [GitHub repository](https://github.com/gammapy/gammapy-extra/tree/master/notebooks).\n", "\n", "**Source files:**\n", "[cta_1dc_introduction.ipynb](../_static/notebooks/cta_1dc_introduction.ipynb) |\n", "[cta_1dc_introduction.py](../_static/notebooks/cta_1dc_introduction.py)\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![CTA first data challenge logo](images/cta-1dc.png)\n", "\n", "# CTA first data challenge (1DC) with Gammapy\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "The [CTA observatory](https://www.cta-observatory.org/) has started a first data challenge (\"CTA 1DC\"), focusing on high-level data analysis. Gammapy is a prototype for the CTA science tools (see [Gammapy proceeding from ICRC 2017](https://github.com/gammapy/icrc2017-gammapy-poster/blob/master/proceeding/gammapy-icrc2017.pdf)), and while many things are work in progress (most importantly: source and background modeling and cube analysis), you can use it already to analyse the simulated CTA data.\n", "\n", "The main page for CTA 1DC is here:\n", "https://forge.in2p3.fr/projects/data-challenge-1-dc-1/wiki (CTA internal)\n", "\n", "There you will find information on 1DC sky models, data access, data organisation, simulation and analysis tools, simulated observations, as well as contact information if you have any questions or comments.\n", "\n", "### This tutorial notebook\n", "\n", "This notebook shows you how to get started with CTA 1DC data and what it contains.\n", "\n", "You will learn how to use Astropy and Gammapy to:\n", "\n", "* access event data\n", "* access instrument response functions (IRFs)\n", "* use index files and the ``gammapy.data.DataStore`` to access all data\n", "* use the observation index file to select the observations you're interested in\n", "* read model XML files from Python (not integrated in Gammapy yet)\n", "\n", "This is to familiarise ourselves with the data files and to get an overview.\n", "\n", "### Further information\n", "\n", "How to analyse the CTA 1DC data with Gammapy (make an image and spectrum) is shown in a second notebook [cta_data_analysis.ipynb](cta_data_analysis.ipynb). If you prefer, you can of course just skim or skip this notebook and go straight to second one.\n", "\n", "More tutorial notebooks for Gammapy are [here](../index.ipynb),\n", "the Gammapy Sphinx docs are at at http://docs.gammapy.org\n", "If you have a Gammapy-related question, please send an email to to Gammapy mailing list at http://groups.google.com/group/gammapy (registration required) or if you have an issue or feature request, file an issue here: https://github.com/gammapy/gammapy/issues/new (free Github account required, takes 1 min to set up).\n", "\n", "Please note that the 1DC data isn't public and results should be shared on CTA Redmine, as described here: https://forge.in2p3.fr/projects/data-challenge-1-dc-1/wiki#Sharing-of-analysis-results. Unfortunately there's no good way yet to share Jupyter notebooks in CTA, as there is e.g. on Github and nbviewer which can show any public noteook." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Notebook and Gammapy Setup\n", "\n", "Before we get started, please execcute the following code cells.\n", "\n", "The first one configures the notebooks so that plots are shown inline (if you don't do this, separate windows will pop up). The second cell imports and checks the version of the packages we will use below. If you're missing some packages, install them and then select \"Kernel -> Restart\" above to restart this notebook.\n", "\n", "In case you're new to Jupyter notebooks: to execute a cell, select it, then type \"SHIFT\" + \"ENTER\"." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "numpy: 1.13.3\n", "astropy: 3.0.dev20763\n", "gammapy: 0.7.dev5344\n" ] } ], "source": [ "import numpy as np\n", "import astropy\n", "import gammapy\n", "\n", "print('numpy:', np.__version__)\n", "print('astropy:', astropy.__version__)\n", "print('gammapy:', gammapy.__version__)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Getting the 1DC data\n", "\n", "You can find infos how to download the 1DC data here: \n", "https://forge.in2p3.fr/projects/data-challenge-1-dc-1/wiki#Data-access\n", "\n", "Overall it's quite large (20 GB) and will take a while to download. You can just download a subset of the \"gps\", \"gc\", \"egal\" or \"agn\" datasets, the ones you're interested in.\n", "\n", "**In this tutorial, we will only use the \"gps\" (Galactic center survey) dataset, so maybe you could start by downloading that first.**\n", "\n", "While you wait, we strongly recommend you go over some [CTA basics](https://www.youtube.com/playlist?list=PLq3ItKoU0hy7ED6m1eve6WIFs7ibcv3YR) as well as some [Python basics](https://www.youtube.com/playlist?list=PL-Qryc-SVnnu0tPV623TFnrAEQLykkZF5) to prepare yourself for this tutorial.\n", "\n", "Got the data?\n", "\n", "Assuming you've followed the instructions, you should have the ``CTADATA`` environment variable pointing to the folder where all data is located. (Gammapy doesn't need or use the ``CALDB`` environment variable.)\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/Users/deil/work/cta-dc/data/1dc/1dc\r\n" ] } ], "source": [ "!echo $CTADATA" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[34mcaldb\u001b[m\u001b[m \u001b[34mdata\u001b[m\u001b[m \u001b[34mindex\u001b[m\u001b[m \u001b[34mmodels\u001b[m\u001b[m \u001b[34mobs\u001b[m\u001b[m\r\n" ] } ], "source": [ "!ls $CTADATA" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can find a description of the directories and files here:\n", "https://forge.in2p3.fr/projects/data-challenge-1-dc-1/wiki/Data_organisation\n", "\n", "A very detailed specification of the data formats is here:\n", "http://gamma-astro-data-formats.readthedocs.io/\n", "\n", "But actually, instead of reading those pages, let's just explore the data and see how to load it with Gammapy ...\n", "\n", "Before we start, just in case the commands above show that you don't have `CTADATA` set:\n", "\n", "* you either have to exit the \"jupyter notebook\" command on your terminal, set the environment variable (I'm using bash and added the command `export CTADATA=/Users/deil/work/cta-dc/data/1dc/1dc` to my `~/.profile` file and then did `source ~/.profile`), then re-start Jupyter and this notebook.\n", "* or you can set the environment variable by uncommentting the code in the following cell, setting the correct path, then executing it." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# import os\n", "# os.environ['CTADATA'] = '/Users/deil/work/cta-dc/data/1dc/1dc'\n", "# !echo $CTADATA\n", "# !ls $CTADATA" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Let's have a look around at the directories and files in `$CTADATA`.\n", "\n", "We will look at the `data` folder with events, the `caldb` folder with the IRFs and the `index` folder with the index files. At the end, we will also mention what the `model` and `obs` folder contains, but they aren't used with Gammapy, at least not at the moment.\n", "\n", "## EVENT data\n", "\n", "First, the EVENT data (RA, DEC, ENERGY, TIME of each photon or hadronic background event) is in the `data/baseline` folder, with one observation per file. The \"baseline\" refers to the assumed CTA array that was used to simulate the observations. The number in the filename is the observation identifier `OBS_ID` of the observation. Observations are ~ 30 minutes, pointing at a fixed location on the sky." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[34mcaldb\u001b[m\u001b[m \u001b[34mdata\u001b[m\u001b[m \u001b[34mindex\u001b[m\u001b[m \u001b[34mmodels\u001b[m\u001b[m \u001b[34mobs\u001b[m\u001b[m\r\n" ] } ], "source": [ "!ls $CTADATA" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[34magn\u001b[m\u001b[m \u001b[34megal\u001b[m\u001b[m \u001b[34mgc\u001b[m\u001b[m \u001b[34mgps\u001b[m\u001b[m\r\n" ] } ], "source": [ "!ls $CTADATA/data/baseline" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "gps_baseline_110000.fits\r\n", "gps_baseline_110001.fits\r\n", "gps_baseline_110002.fits\r\n" ] } ], "source": [ "!ls $CTADATA/data/baseline/gps | head -n3" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 3270\n", " 11G\t/Users/deil/work/cta-dc/data/1dc/1dc/data/baseline/gps\n" ] } ], "source": [ "# There's 3270 observations and 11 GB of event data for the gps dataset\n", "!ls $CTADATA/data/baseline/gps | wc -l\n", "!du -hs $CTADATA/data/baseline/gps" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's open up the first event list using the Gammapy `EventList` class, which contains the ``EVENTS`` table data via the ``table`` attribute as an Astropy `Table` object." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "\n" ] } ], "source": [ "from gammapy.data import EventList\n", "events = EventList.read('$CTADATA/data/baseline/gps/gps_baseline_110000.fits')\n", "print(type(events))\n", "print(type(events.table))" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Row index=0\n", "\n", "\n", "\n", "\n", "\n", "
EVENT_IDTIMERADECENERGYDETXDETYMC_ID
sdegdegTeVdegdeg
uint32float64float32float32float32float32float32int32
1662774403.255-173.287-62.40990.04861491.607950.2580162
" ], "text/plain": [ "\n", "EVENT_ID TIME RA DEC ENERGY DETX DETY MC_ID\n", " s deg deg TeV deg deg \n", " uint32 float64 float32 float32 float32 float32 float32 int32\n", "-------- ------------- -------- -------- --------- ------- -------- -----\n", " 1 662774403.255 -173.287 -62.4099 0.0486149 1.60795 0.258016 2" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# First event (using [] for \"indexing\")\n", "events.table[0]" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Table length=2\n", "\n", "\n", "\n", "\n", "\n", "\n", "
EVENT_IDTIMERADECENERGYDETXDETYMC_ID
sdegdegTeVdegdeg
uint32float64float32float32float32float32float32int32
1662774403.255-173.287-62.40990.04861491.607950.2580162
2662774459.184-172.62-63.03460.06110350.9790580.5551662
" ], "text/plain": [ "\n", "EVENT_ID TIME RA DEC ENERGY DETX DETY MC_ID\n", " s deg deg TeV deg deg \n", " uint32 float64 float32 float32 float32 float32 float32 int32\n", "-------- ------------- -------- -------- --------- -------- -------- -----\n", " 1 662774403.255 -173.287 -62.4099 0.0486149 1.60795 0.258016 2\n", " 2 662774459.184 -172.62 -63.0346 0.0611035 0.979058 0.555166 2" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# First few events (using [] for \"slicing\")\n", "events.table[:2]" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "# Event times can be accessed as Astropy Time objects\n", "print(type(events.time))" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
OBS_IDGLON_PNTGLAT_PNTIRF
int64float64float64bytes13
110380359.999991204-1.29999593791South_z20_50h
110381359.999991204-1.29999593791South_z20_50h
110382359.999991204-1.29999593791South_z20_50h
110383359.999991204-1.29999593791South_z20_50h
110384359.999991204-1.29999593791South_z20_50h
" ], "text/plain": [ "\n", "OBS_ID GLON_PNT GLAT_PNT IRF \n", "int64 float64 float64 bytes13 \n", "------ ------------- -------------- -------------\n", "110380 359.999991204 -1.29999593791 South_z20_50h\n", "110381 359.999991204 -1.29999593791 South_z20_50h\n", "110382 359.999991204 -1.29999593791 South_z20_50h\n", "110383 359.999991204 -1.29999593791 South_z20_50h\n", "110384 359.999991204 -1.29999593791 South_z20_50h" ] }, "execution_count": 52, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Look at the first few\n", "table[['OBS_ID', 'GLON_PNT', 'GLAT_PNT', 'IRF']][:5]" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'South_z20_50h'}" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Check which IRFs were used ... it's all south and 20 deg zenith angle\n", "set(table['IRF'])" ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0,0.5,'Galactic latitude (deg)')" ] }, "execution_count": 54, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAGRpJREFUeJzt3X2UZVV55/HvDxQaovhGqwi0DYYE\nMWMEKogvMWA0I6h0jJCAGRXFMDgyKk7GYcaViJq1VJZOZvkuGhVNokYTBXxDBQk6itIwYPMi2vIS\nWljSoqJERIFn/jin5FpU3drdVbfu7arvZ6276rzse86zb3XXc8/e5+ydqkKSpBbbjTsASdK2w6Qh\nSWpm0pAkNTNpSJKamTQkSc1MGpKkZiYNSVIzk4YkqZlJQ5LU7F7jDmCx7brrrrV27dpxhyFJ25SL\nLrroB1W1er5yyy5prF27lvXr1487DEnapiS5rqWczVOSpGYmDUlSM5OGJKmZSUOS1MykIUlqZtKQ\nJDVbdrfcLsTakz99j23XvuHpY4hEktos9d8trzR6s33ww7ZL0riN4++WSUOS1MykIUlqZtKQJDUz\naUiSmpk0enPdbeDdU5Im1Tj+bqWqRnbwcZiamipHuZWkLZPkoqqamq+cVxqSpGYmDUlSM5OGJKmZ\nSUOS1MykIUlqZtKQJDUzaUiSmpk0JEnNTBqSpGYmDUlSM5OGJKmZSUOS1MykIUlqZtKQJDUba9JI\n8r4kNyW5bI79SfKWJBuTfDPJAUsdoyTpbuO+0vgA8LQh+w8D9ulfxwPvXIKYJElzGGvSqKrzgR8O\nKbIO+GB1LgDun2S3pYlOkjTTuK805rM7cP3A+qZ+myRpDCY9aWSWbfeYnzbJ8UnWJ1m/efPmJQhL\nklamSU8am4A9B9b3AG6YWaiqTquqqaqaWr169ZIFJ0krzaQnjTOB5/V3UR0M3FJVN447KElaqe41\nzpMn+TBwCLBrkk3Aq4F7A1TVu4DPAIcDG4GfAS8YT6SSJBhz0qiqY+bZX8BLligcSdI8Jr15SpI0\nQUwakqRmJg1JUjOThiSpmUlDktTMpCFJambSkCQ1M2lIkpqZNCRJzUwakqRmJg1JUjOThiSpmUlD\nktTMpCFJambSkCQ1a5pPI8kDgIcBtwHXVtVdI41KkjSR5kwaSe5HNwHSMcAOwGZgFfCQJBcA76iq\nLy1JlJKkiTDsSuPjwAeB36+qHw/uSHIg8Nwke1fV340yQEnS5JgzaVTVU4fsuwi4aCQRSZIm1rx9\nGkkOmGXzLcB1VXXH4ockSZpULR3h7wAOAL4JBPidfvlBSU6oqs+PMD5J0gRpueX2WmD/qpqqqgOB\n/YHLgKcAp44wNknShGlJGvtW1eXTK1V1BV0SuXp0YUmSJlFL89RVSd4JfKRf/zPg20l2BH45ssgk\nSROn5UrjWGAj8HLgJODqftsvgUNHFZgkafLMe6VRVbcleQfwqaq6asbuW0cTliRpEs17pZHkCOAS\n4HP9+mOSnDnqwCRJk6eleerVwEHAjwGq6hJg7QhjkiRNqJakcUdV3TLySCRJE68laVyW5DnA9kn2\nSfJW4KuLcfIkT0tyVZKNSU6eZf+xSTYnuaR/vWgxzitJ2jotSeO/Ao8Cbgc+DPyE7k6qBUmyPfB2\n4DBgP+CYJPvNUvSjVfWY/vXehZ5XkrT1Wu6e+hnwqv61mA4CNk4/JJjkI8A64IpFPo8kaZEMm0/j\nLKDm2l9VRyzw3LsD1w+sbwIeO0u5Zyd5EvBt4KSqun5mgSTHA8cDrFmzZoFhSZLmMqx56k3Am4Fr\n6Gbse0//upVu7KmFyizbZiaps4C1VfVo4IvA6bMdqKpO68fGmlq9evUihCZJms2w+TT+FSDJ66rq\nSQO7zkpy/iKcexOw58D6HsANM2K4eWD1PcAbF+G8kqSt1NIRvjrJ3tMrSfYCFuPr/IXAPkn2SrID\ncDTwaw8NJtltYPUI4MpFOK8kaSu1DFh4EnBekulRbdfS9x8sRFXdkeRE4Gxge+B9VXV5ktcC66vq\nTOCl/RPpdwA/pBvzSpI0Jqmas6/77kLdiLb79qvfqqrbRxrVAkxNTdX69evHHYYkbVOSXFRVU/OV\nm7N5KskTp5er6vaqurR/3d7v3yXJ7yxOuJKkbcGw5qlnJzmVbqDCi4DNwCrgN+mGRH848N9GHqEk\naWIMu3vqpCQPAI4EjgJ2o7v19krg3VX1laUJUZI0KYZ2hFfVj7j7+QxJ0grXcsutJEmASUOStAVM\nGpKkZi3Tve6c5K+SvKdf3yfJM0YfmiRp0rRcabyfbi6Nx/Xrm4C/GVlEkqSJ1ZI0HlFVpwK/BKiq\n25h9hFpJ0jLXkjR+kWQn+mHLkzyC7spDkrTCtAxY+Gq6p8L3TPIPwBNw4EBJWpFapnv9QpKLgYPp\nmqVeVlU/GHlkkqSJM2y61wNmbLqx/7kmyZqqunh0YUmSJtGwK4039z9XAVPApXRXGo8Gvg48cY73\nSZKWqTk7wqvq0Ko6FLgOOKCfg/tAYH9g41IFKEmaHC13T+1bVRumV6rqMuAxowtJkjSpWu6eujLJ\ne4G/p7vt9j/hXN2StCK1JI0XAC8GXtavnw+8c2QRSZImVssttz8H/rZ/SZJWsHmTRpJr6J8GH1RV\ne48kIknSxGppnpoaWF5FN/XrA0cTjiRpks1791RV3Tzw+l5V/R/gyUsQmyRpwrQ0Tw0+Gb4d3ZXH\nfUcWkSRpYrU0T715YPkO4BrgT0cTjiRpkrUkjeOq6urBDUn2GlE8kqQJ1vJE+Mcbt0mSlrlho9zu\nCzwKuF+SPxnYtQvdXVSSpBVm2JXGbwPPAO4PPHPgdQDwF4tx8iRPS3JVko1JTp5l/45JPtrv/3qS\ntYtxXknS1pnzSqOqzgDOSPK4qvraYp84yfbA24GnApuAC5OcWVVXDBQ7DvhRVf1mkqOBNwJ/ttix\nSJLaDGueemVVnQo8J8kxM/dX1UsXeO6DgI3TnexJPgKsAwaTxjrglH7548DbkqSq7vGEuiRp9Ibd\nPTU9ku36EZ17d+D6gfVNwGPnKlNVdyS5BXgQ4HSzkjQGw5qnzuoXf1ZVHxvcl+SoRTh3ZjvtVpQh\nyfHA8QBr1qxZeGSSpFm13HL7Pxu3balNwJ4D63sAN8xVJsm9gPsBP5x5oKo6rZ9ZcGr16tWLEJok\naTbD+jQOAw4Hdk/yloFdu9A9Gb5QFwL79A8Kfg84GnjOjDJnAs8HvgYcCZxrf4Ykjc+wPo0b6Poz\njgAuGtj+U+CkhZ6476M4ETgb2B54X1VdnuS1wPqqOhP4O+BDSTbSXWEcvdDzSpK2Xub74p7k3lX1\nyyWKZ8GmpqZq/fpR9d1L0vKU5KKqmpqvXMvYU2uTvB7Yj4EnwZ2ESZJWnpaO8PfTzQl+B3Ao8EHg\nQ6MMSpI0mVqSxk5VdQ5dU9Z1VXUKTsIkSStSS/PUz5NsB3yn77j+HvDg0YYlSZpELVcaLwd2Bl4K\nHAg8l+42WEnSCjPvlUZVXdgv3gq8YLThSJIm2bCH+85iliE7plXVESOJSJI0sYZdabxpyaKQJG0T\nhg1Y+K9LGYgkafK1dIRLkgSYNCRJW8CkIUlqNm/SSPKFJPcfWH9AkrNHG5YkaRK1XGnsWlU/nl6p\nqh/hE+GStCK1JI27kvxqDtUkD2fI8xuSpOWrZeypVwFfSTJ9C+6T6OfjliStLC3DiHwuyQHAwUCA\nk6rqByOPTJI0ceZsnkqyb//zAGAN3fSv3wPW9NskSSvMsCuNV9A1Q715ln2Fc2pI0oozbBiR6X6L\nw6rq54P7kqya5S2SpGWu5e6przZukyQtc8OGRn8osDuwU5L96TrBAXahm5RJkrTCDOvT+I/AscAe\ndP0a00njJ8D/Gm1YkqRJNKxP43Tg9CTPrqp/XsKYJEkTqqVP48BZxp76mxHGJEmaUC1J47BZxp46\nfHQhSZImVUvS2D7JjtMrSXYCdhxSXpK0TLWMPfX3wDlJ3k/3UN8LgdNHGpUkaSK1jD11apINwB/S\n3UH1uqpyPg1JWoFarjSoqs8Cn12skyZ5IPBRYC1wLfCnfV/JzHJ3Ahv61X+rqiMWKwZJ0pZrmbnv\n4CQXJrk1yS+S3JnkJws878nAOVW1D3BOvz6b26rqMf3LhCFJY9bSEf424BjgO8BOwIuAty7wvOu4\nu1/kdOCPF3g8SdISaEkaVNVGYPuqurOq3g8cusDzPqSqbuyPfSNzTx+7Ksn6JBckMbFI0pi19Gn8\nLMkOwCVJTgVuBH5jvjcl+SLw0Fl2vWoL4ltTVTck2Rs4N8mGqvruLOc6nn42wTVr1szcLUlaJC1J\n47nA9sCJwEnAnsCz53tTVT1lrn1Jvp9kt6q6McluwE1zHOOG/ufVSc4D9gfukTSq6jTgNICpqSnn\nL5ekEZm3eaqqrquq26rqJ1X1mqp6Rd9ctRBnAs/vl58PnDGzQD9cyY798q7AE4ArFnheSdICDBsa\nfQPdw3yzqqpHL+C8bwD+KclxwL8BR/XnnAJOqKoXAY8E3p3kLrrk9oaqMmlI0hgNa556xqhOWlU3\n0z0sOHP7erq7s6iqrwL/YVQxSJK23LCh0a9bykAkSZNvXA/3SZK2QeN6uE+StA1qHXtqY5Ltq+pO\n4P1JvjriuCRJE2hkD/dJkpafluap5/blTgT+ncaH+yRJy0/LfBrTd1H9HHjNaMORJE2yOa80kqxL\n8pKB9a8nubp/Hbk04UmSJsmw5qlX0g33MW1H4PeAQ4AXjzAmSdKEGtY8tUNVXT+w/pX+Se6bk9gR\nLkkr0LArjQcMrlTViQOrq0cTjiRpkg1LGl9P8hczNyb5z8A3RheSJGlSDWueOgn4ZJLnABf32w6k\n69twFj1JWoGGDVh4E/D4JE8GHtVv/nRVnbskkUmSJk7LcxrnAiYKSVLTE+GSJAEmDUnSFjBpSJKa\nmTQkSc1MGpKkZiYNSVIzk4YkqZlJQ5LUzKQhSWpm0pAkNTNpSJKamTQkSc1MGpKkZiYNSVKzsSSN\nJEcluTzJXUmmhpR7WpKrkmxMcvJSxihJuqdxXWlcBvwJcP5cBZJsD7wdOAzYDzgmyX5LE54kaTbz\nTsI0ClV1JUCSYcUOAjZW1dV92Y8A64ArRh6gJGlWk9ynsTtw/cD6pn7bPSQ5Psn6JOs3b968JMFJ\n0ko0siuNJF8EHjrLrldV1Rkth5hlW81WsKpOA04DmJqamrWMJGnhRpY0quopCzzEJmDPgfU9gBsW\neExJ0gJMcvPUhcA+SfZKsgNwNHDmmGOSpBVtXLfcPivJJuBxwKeTnN1vf1iSzwBU1R3AicDZwJXA\nP1XV5eOIV5LUGdfdU58APjHL9huAwwfWPwN8ZglDkyQNMcnNU5KkCWPSkCQ1M2lIkpqZNCRJzUwa\nkqRmJg1JUjOThiSpmUlDktTMpCFJambSkCQ1M2lIkpqZNCRJzUwakqRmJg1JUrOxDI0urTRrT/70\nPbZd+4anjyESaWG80pBGbLaEMWy7NMlMGpKkZiYNSVIzk4YkqZlJQ5LUzKQhjdhcd0l595S2Rd5y\nKy0BE4SWC680JEnNTBqSpGYmDUlSM5OGJKmZSUOS1MykIUlqlqoadwyLKslm4LoFHmZX4AeLEM64\nWY/Js1zqslzqAcunLgutx8OravV8hZZd0lgMSdZX1dS441go6zF5lktdlks9YPnUZanqYfOUJKmZ\nSUOS1MykMbvTxh3AIrEek2e51GW51AOWT12WpB72aUiSmnmlIUlqtuKTRpKjklye5K4kc955kOTa\nJBuSXJJk/VLG2GoL6vK0JFcl2Zjk5KWMsUWSByb5QpLv9D8fMEe5O/vfxyVJzlzqOIeZ7zNOsmOS\nj/b7v55k7dJHOb+GehybZPPA7+FF44hzPknel+SmJJfNsT9J3tLX85tJDljqGFs01OOQJLcM/D7+\netGDqKoV/QIeCfw2cB4wNaTctcCu4453oXUBtge+C+wN7ABcCuw37thnxHgqcHK/fDLwxjnK3Tru\nWLf2Mwb+C/Cufvlo4KPjjnsr63Es8LZxx9pQlycBBwCXzbH/cOCzQICDga+PO+atrMchwKdGGcOK\nv9Koqiur6qpxx7EYGutyELCxqq6uql8AHwHWjT66LbIOOL1fPh344zHGsjVaPuPBOn4c+MMkWcIY\nW2wL/1aaVNX5wA+HFFkHfLA6FwD3T7Lb0kTXrqEeI7fik8YWKODzSS5Kcvy4g1mA3YHrB9Y39dsm\nyUOq6kaA/ueD5yi3Ksn6JBckmaTE0vIZ/6pMVd0B3AI8aEmia9f6b+XZfZPOx5PsuTShLbpt4f9F\nq8cluTTJZ5M8arEPviJm7kvyReChs+x6VVWd0XiYJ1TVDUkeDHwhybf6rL+kFqEus32bXfJb6IbV\nYwsOs6b/newNnJtkQ1V9d3EiXJCWz3gifg/zaInxLODDVXV7khPorp6ePPLIFt+28PtocTHdcCC3\nJjkc+CSwz2KeYEUkjap6yiIc44b+501JPkF36b7kSWMR6rIJGPw2uAdwwwKPucWG1SPJ95PsVlU3\n9k0EN81xjOnfydVJzgP2p2uDH7eWz3i6zKYk9wLux5ibHWYxbz2q6uaB1fcAb1yCuEZhIv5fLFRV\n/WRg+TNJ3pFk16patLG1bJ5qkOQ3ktx3ehn4I2DWuxe2ARcC+yTZK8kOdJ2wE3XnEV08z++Xnw/c\n4woqyQOS7Ngv7wo8AbhiySIcruUzHqzjkcC51fdkTpB56zGj3f8I4MoljG8xnQk8r7+L6mDglukm\n0m1JkodO940lOYjub/zNw9+1hcZ9N8C4X8Cz6L5l3A58Hzi73/4w4DP98t50d45cClxO1xQ09ti3\npi79+uHAt+m+lU9cXeja9s8BvtP/fGC/fQp4b7/8eGBD/zvZABw37rhn1OEenzHwWuCIfnkV8DFg\nI/ANYO9xx7yV9Xh9/3/iUuBLwL7jjnmOenwYuBH4Zf9/5DjgBOCEfn+At/f13MCQOyknvB4nDvw+\nLgAev9gx+ES4JKmZzVOSpGYmDUlSM5OGJKmZSUOS1MykIUlqZtLQWCV5SJJ/THJ1P0TL15I8a573\nrJ1rlM+G8x2b5GED6+9Nsl/jew9J8qmtOe88x31tkqf0yy9PsvNWHOPWLSyfJOcm2WWWfack+cst\njaF/7+okn9ua92rbYNLQ2PQPIX0SOL+q9q6qA+keINtjhKc9lu65FQCq6kVVNdaHAqvqr6vqi/3q\ny4EtThpb4XDg0hp4gngxVNVm4MYkT1jM42pymDQ0Tk8GflFV75reUFXXVdVb4VdXFF9OcnH/evzM\nAwwrk+SV6eZAuTTJG5IcSfeA4D/0cw3slOS89HOPpJs74uK+/DnDAk8358cn+4H6Lkjy6H77Kf2c\nB+f1V08vHXjPXyX5Vro5Qj48/W0+yQeSHNmXfRjwpSRf6vfdOvD+I5N8oF/eq78quzDJ62bE9t/7\n7d9M8po5qvDnDDxpn+RV6ebN+CLd8PrT2x+R5HP9VeCXk+w7sP2C/jyvnXGl88n++FqOxv2Eo6+V\n+wJeCvztkP07A6v65X2A9f3yWvr5BIaUOQz4KrBzvz79VPl5DDztO70OrKYb5XSvwfIz4jmEfq4C\n4K3Aq/vlJwOX9Mun9OfdEdiVbgiHe/fnuATYCbgv3dPuf9m/5wPAkf3ytQzM28LAnCF0w418oF8+\nE3hev/yS6XJ0Q9ycRveE83bAp4AnzVKX64D79ssH0j0FvTOwC91T6tOxnQPs0y8/lm64E/rjHtMv\nnzAjzt2BDeP+9+VrNK8VMWChtg1J3g48ke7q4/fo/ti+LcljgDuB35rlbXOVeQrw/qr6GUBVzTcY\n4MF0zWTXNJZ/IvDsvuy5SR6U5H79vk9X1e3A7UluAh7Slz+jqm7r63rWPMefzxOmzw98iLsHCvyj\n/vX/+vX70CXTmYNrPrCqftov/z7wienPKv0siEnuQzdcy8dy91QfO/Y/H8fd85z8I/CmgWPfxEAT\noJYXk4bG6XLu/sNHVb2kH3xwejrdk+jG0Ppdum/NP5/lGHOVCVs2tPXWlJ9p+v23D2y7k+7/2dZO\nsDQY06oh+wbjen1VvXue496RZLuqumvIsbYDflxVj2kL9dfivG0L36NthH0aGqdz6SZSevHAtsFO\n4PsBN/Z/2J5LN/3oTHOV+Tzwwuk7kZI8sN/+U7rmoZm+BvxBkr1mlJ/L+fTt9kkOAX5QwzuVvwI8\nM8mq/hv80+coNzO+7yd5ZJLt6AaknPZ/6W4agF/vPzibrt736WPbPd0cMDNdRTcQ53RdntX38dwX\neCb8apjta5Ic1R8rSX63f88F3J3wj+bX/Rbb7ijQmodJQ2NTVUXXxPEHSa5J8g26SXz+R1/kHcDz\nk1xA94fo32c5zKxlqupzdO3+65NcAkzfQvoB4F3THeEDsWwGjgf+JcmlwEfnCf8UYCrJN4E3cPcw\n53PV9cI+nkuBf6G7mrpllqKnAZ+d7ginmyP9U3QJdnCo7pcBL0lyIV3inD7P5+mai76WZAPdVLKz\nJclP0/XRUFUX09X3EuCfgS8PlPtz4Lj+M7mcu6d7fTnwiv53ttuMuhzaH1/LkKPcSkskyX2qm1Ft\nZ7pv98f3f7DHEctudHNiP3Ur378zcFtVVZKj6TrF1/X7zgfWVdWPFi9iTQr7NKSlc1q6BwlXAaeP\nK2FAN/d6kvck2WWeZrW5HEh3A0KAHwMvhO7hPuB/mzCWL680JEnN7NOQJDUzaUiSmpk0JEnNTBqS\npGYmDUlSM5OGJKnZ/wdzeY9oTv4FDQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Check the pointing positions\n", "# The grid pointing positions at GLAT = +- 1.2 deg are visible\n", "from astropy.coordinates import Angle\n", "plt.scatter(Angle(table['GLON_PNT'], unit='deg').wrap_at('180 deg'), table['GLAT_PNT'])\n", "plt.xlabel('Galactic longitude (deg)')\n", "plt.ylabel('Galactic latitude (deg)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Load data\n", "\n", "Once you have selected the observations of interest, use the `DataStore` to load the data and IRF for those observations. Let's say we're interested in `OBS_ID=110000`." ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [], "source": [ "obs = data_store.obs(obs_id=110000)" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Info for OBS_ID = 110000\n", "- Start time: 59215.50\n", "- Pointing pos: RA 186.16 deg / Dec -64.02 deg\n", "- Observation duration: 1800.0 s\n", "- Dead-time fraction: 2.000 %\n", "\n" ] } ], "source": [ "print(obs)" ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 57, "metadata": {}, "output_type": "execute_result" } ], "source": [ "obs.events" ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Table length=5\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
EVENT_IDTIMERADECENERGYDETXDETYMC_ID
sdegdegTeVdegdeg
uint32float64float32float32float32float32float32int32
1662774403.255-173.287-62.40990.04861491.607950.2580162
2662774459.184-172.62-63.03460.06110350.9790580.5551662
3662774479.309-170.593-63.47150.0461420.5105491.451442
4662774490.724-175.637-63.64090.04935670.36688-0.7959082
5662774506.865-172.124-63.40130.05054010.6073640.77022
" ], "text/plain": [ "\n", "EVENT_ID TIME RA DEC ENERGY DETX DETY MC_ID\n", " s deg deg TeV deg deg \n", " uint32 float64 float32 float32 float32 float32 float32 int32\n", "-------- ------------- -------- -------- --------- -------- --------- -----\n", " 1 662774403.255 -173.287 -62.4099 0.0486149 1.60795 0.258016 2\n", " 2 662774459.184 -172.62 -63.0346 0.0611035 0.979058 0.555166 2\n", " 3 662774479.309 -170.593 -63.4715 0.046142 0.510549 1.45144 2\n", " 4 662774490.724 -175.637 -63.6409 0.0493567 0.36688 -0.795908 2\n", " 5 662774506.865 -172.124 -63.4013 0.0505401 0.607364 0.7702 2" ] }, "execution_count": 58, "metadata": {}, "output_type": "execute_result" } ], "source": [ "obs.events.table[:5]" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 59, "metadata": {}, "output_type": "execute_result" } ], "source": [ "obs.aeff" ] }, { "cell_type": "code", "execution_count": 60, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 60, "metadata": {}, "output_type": "execute_result" } ], "source": [ "obs.edisp" ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 61, "metadata": {}, "output_type": "execute_result" } ], "source": [ "obs.psf" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's an example how to loop over many observations and analyse them: [cta_1dc_make_allsky_images.py](https://github.com/gammasky/cta-dc/blob/master/data/cta_1dc_make_allsky_images.py)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Model XML files\n", "\n", "The 1DC sky model is distributed as a set of XML files, which in turn link to a ton of other FITS and text files." ] }, { "cell_type": "code", "execution_count": 62, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "model_Arp220.xml\n", "model_agn.xml\n", "model_bkg.xml\n", "model_dm_dSphs.xml\n", "model_dm_gc.xml\n", "model_fermi_bubbles.xml\n", "model_galactic_binaries.xml\n", "model_galactic_bright.xml\n", "model_galactic_composite.xml\n", "model_galactic_extended.xml\n", "model_galactic_pulsars.xml\n", "model_galactic_pwn.xml\n", "model_galactic_snr.xml\n", "model_iem.xml\n", "models_agn.xml\n", "models_egal.xml\n", "models_gc.xml\n", "models_gps.xml\n" ] } ], "source": [ "!ls $CTADATA/models/*.xml | xargs -n 1 basename" ] }, { "cell_type": "code", "execution_count": 63, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " \r\n", " \r\n", " \r\n", " \r\n", " \r\n", " \r\n", " \r\n", " \r\n", " \r\n", " \r\n", " \r\n", " \r\n", " \r\n", " \r\n", " \r\n", " \r\n", " \r\n", " \r\n", " \r\n", "\r\n" ] } ], "source": [ "# This is what the XML file looks like\n", "!tail -n 20 $CTADATA/models/models_gps.xml" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At the moment, you cannot read and write these sky model XML files with Gammapy.\n", "\n", "There are multiple reasons why this XML serialisation format isn't implemented in Gammapy yet (all variants of \"it sucks\"):\n", "\n", "* XML is tedious to read and write for humans.\n", "* The format is too strict in places: there are many use cases where \"min\", \"max\", \"free\" and \"error\" aren't needed, so it should be possible to omit them (see e.g. dummy values above)\n", "* The parameter \"scale\" is an implementation detail that very few optimisers need. There's no reason to bother all gamma-ray astronomers with separating value and scale in model input and result files.\n", "* The \"unit\" is missing. Pretty important piece of information, no?\n", " All people working on IACTs use \"TeV\", why should CTA be forced to use \"MeV\"?\n", "* Ad-hoc extensions that keep being added and many models can't be put in one file (see extra ASCII and FITS files with tables or other info, e.g. pulsar models added for CTA 1DC). Admittedly complex model serialisation is an intrinsically hard problem, there is not simple and flexible solution.\n", "\n", "Also, to be honest, I also want to say / admit:\n", "\n", "* A model serialisation format is useful, even if it will never cover all use cases (e.g. energy-dependent morphology or an AGN with a variable spectrum, or complex FOV background models).\n", "* The Gammapy model package isn't well-implemented or well-developed yet.\n", "* So far users were happy to write a few lines of Python to define a model instead of XML.\n", "* Clearly with CTA 1DC now using the XML format there's an important reason to support it, there is the legacy of Fermi-LAT using it and ctools using / extending it for CTA.\n", "\n", "So what now?\n", "\n", "* In Gammapy, we have started to implement a modeling package in `gammapy.utils.modeling`, with the goal of supporting this XML format as one of the serialisation formats. It's not very far along, will be a main focus for us, help is welcome.\n", "* In addition we would like to support a second more human-friendly model format that looks something like [this](https://github.com/gammapy/gamma-cat/blob/b651de8d1d793e924764ffb13c8ec189bce9ea7d/input/data/2006/2006A%2526A...457..899A/tev-000025.yaml#L11)\n", "* For now, you could use Gammalib to read the XML files, or you could read them directly with Python. The Python standard library contains [ElementTree](https://docs.python.org/3/library/xml.etree.elementtree.html) and there's [xmltodict](https://github.com/martinblech/xmltodict) which simply hands you back the XML file contents as a Python dictionary (containing a very nested hierarchical structure of Python dict and list objects and strings and numbers.\n", "\n", "As an example, here's how you can read an XML sky model and access the spectral parameters of one source (the last, \"Arp200\" visible above in the XML printout) and create a [gammapy.spectrum.models.PowerLaw](http://docs.gammapy.org/dev/api/gammapy.spectrum.models.PowerLaw.html) object." ] }, { "cell_type": "code", "execution_count": 64, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[OrderedDict([('@name', 'Prefactor'),\n", " ('@value', '6'),\n", " ('@error', '0'),\n", " ('@scale', '1e-20'),\n", " ('@min', '0'),\n", " ('@free', '1')]),\n", " OrderedDict([('@name', 'Index'),\n", " ('@value', '2.2'),\n", " ('@error', '-0'),\n", " ('@scale', '-1'),\n", " ('@min', '-4.54545'),\n", " ('@max', '4.54545'),\n", " ('@free', '1')]),\n", " OrderedDict([('@name', 'Scale'),\n", " ('@value', '1'),\n", " ('@scale', '1000000'),\n", " ('@free', '0')])]" ] }, "execution_count": 64, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Read XML file and access spectrum parameters\n", "from gammapy.extern import xmltodict\n", "filename = os.path.join(os.environ['CTADATA'], 'models/models_gps.xml')\n", "data = xmltodict.parse(open(filename).read())\n", "data = data['source_library']['source'][-1]\n", "data = data['spectrum']['parameter']\n", "data" ] }, { "cell_type": "code", "execution_count": 65, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PowerLaw\n", "\n", "Parameters: \n", "\n", "\t name value error unit min max frozen\n", "\t--------- ---------- ----- --------------- --- --- ------\n", "\t index -2.200e+00 nan nan nan False\n", "\tamplitude 6.000e-20 nan 1 / (cm2 MeV s) nan nan False\n", "\treference 1.000e+06 nan MeV nan nan True\n" ] } ], "source": [ "# Create a spectral model the the right units\n", "from astropy import units as u\n", "from gammapy.spectrum.models import PowerLaw\n", "par_to_val = lambda par: float(par['@value']) * float(par['@scale'])\n", "spec = PowerLaw(\n", " amplitude=par_to_val(data[0]) * u.Unit('cm-2 s-1 MeV-1'),\n", " index=par_to_val(data[1]),\n", " reference=par_to_val(data[2]) * u.Unit('MeV'),\n", ")\n", "print(spec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercises\n", "\n", "* Easy: Go over this notebook again, and change the data / code a little bit:\n", " * Pick another run (any run you like)\n", " * Plot the energy-offset histogram of the events separately for gammas and background\n", "\n", "* Medium difficulty: Find all runs within 1 deg of the Crab nebula.\n", " * Load the \"all\" index file via the `DataStore`, then access the ``.obs_table``, then get an array-valued ``SkyCoord`` with all the observation pointing positions.\n", " * You can get the Crab nebula position with `astropy.coordinates.SkyCoord` via ``SkyCoord.from_name('crab')`` \n", " * Note that to compute the angular separation between two ``SkyCoord`` objects can be computed via ``separation = coord1.separation(coord2)``.\n", "\n", "* Hard: Find the PeVatrons in the 1DC data, i.e. the ~ 10 sources that are brightest at high energies (e.g. above 10 TeV).\n", " * This is difficult, because it's note clear how to best do this, and also because it's time-consuming to crunch through all relevant data for any given method.\n", " * One idea could be to go brute-force through **all** events, select the ones above 10 TeV and stack them all into one table. Then make an all-sky image and run a peak finder, or use an event cluster-finding method e.g. from scikit-learn.\n", " * Another idea could be to first make a list of targets of interest, either from the CTA 1DC sky model or from gamma-cat, compute the model integral flux above 10 TeV and pick candidates that way, then run analyses." ] }, { "cell_type": "code", "execution_count": 66, "metadata": {}, "outputs": [], "source": [ "# start typing here ..." ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## What next?\n", "\n", "* This notebook gave you an overview of the 1DC files and showed you have to access and work with them using Gammapy.\n", "* To see how to do analysis, i.e. make a sky image and spectrum, see [cta_data_analysis.ipynb](cta_data_analysis.ipynb) next." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.3" }, "nbsphinx": { "orphan": true } }, "nbformat": 4, "nbformat_minor": 1 }