{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "
\n", "**This is a fixed-text formatted version of a Jupyter notebook.**\n", "\n", "- Try online [![Binder](https://mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.9?urlpath=lab/tree/cta_1dc_introduction.ipynb)\n", "- You can contribute with your own notebooks in this\n", "[GitHub repository](https://github.com/gammapy/gammapy/tree/master/tutorials).\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" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "In 2017 and 2018, the [CTA observatory](https://www.cta-observatory.org/) did a first data challenge, called CTA DC-1, where CTA high-level data was simulated assuming a sky model and using CTA IRFs.\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." ] }, { "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.15.4\n", "astropy: 3.0.5\n", "gammapy: 0.9\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": [ "## DC-1 data\n", "\n", "In this and other Gammapy tutorials we will only access a few data files from CTA DC-1 from `$GAMMAPY_DATA/cta-1dc` that you will have if you followed the \"Getting started with Gammapy\" instructions and executed `gammapy download tutorials`.\n", "\n", "Information how to download more or all of the DC-1 data (in total 20 GB) is here:\n", "https://forge.in2p3.fr/projects/data-challenge-1-dc-1/wiki#Data-access\n", "\n", "Working with that data with Gammapy is identical to what we show here, except that the recommended way to do it is to point `DataStore.read` at `$CTADATA/index/gps` or whatever dataset or files you'd like to access there." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/Users/deil/work/gammapy-tutorials/datasets/cta-1dc\r\n" ] } ], "source": [ "!echo $GAMMAPY_DATA/cta-1dc" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "README.md \u001b[34mcaldb\u001b[m\u001b[m \u001b[34mdata\u001b[m\u001b[m \u001b[34mindex\u001b[m\u001b[m make.py\r\n" ] } ], "source": [ "!ls $GAMMAPY_DATA/cta-1dc" ] }, { "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` 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": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "gps_baseline_110380.fits\r\n", "gps_baseline_111140.fits\r\n", "gps_baseline_111159.fits\r\n", "gps_baseline_111630.fits\r\n" ] } ], "source": [ "!ls -1 $GAMMAPY_DATA/cta-1dc/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": 6, "metadata": {}, "outputs": [], "source": [ "from gammapy.data import EventList\n", "\n", "path = \"$GAMMAPY_DATA/cta-1dc/data/baseline/gps/gps_baseline_110380.fits\"\n", "events = EventList.read(path)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "gammapy.data.event_list.EventList" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(events)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "astropy.table.table.Table" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type(events.table)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Row index=0\n", "\n", "\n", "\n", "\n", "\n", "
EVENT_IDTIMERADECENERGYDETXDETYMC_ID
sdegdegTeVdegdeg
uint32float64float32float32float32float32float32int32
1664502403.0454683-92.63541-30.5148540.03902182-0.9077294-0.27276932
" ], "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 664502403.0454683 -92.63541 -30.514854 0.03902182 -0.9077294 -0.2727693 2" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# First event (using [] for \"indexing\")\n", "events.table[0]" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Table length=2\n", "\n", "\n", "\n", "\n", "\n", "\n", "
EVENT_IDTIMERADECENERGYDETXDETYMC_ID
sdegdegTeVdegdeg
uint32float64float32float32float32float32float32int32
1664502403.0454683-92.63541-30.5148540.03902182-0.9077294-0.27276932
2664502405.2579999-92.64103-28.2627280.0307963711.3443842-0.28383982
" ], "text/plain": [ "\n", "EVENT_ID TIME RA DEC ... DETX DETY MC_ID\n", " s deg deg ... deg deg \n", " uint32 float64 float32 float32 ... float32 float32 int32\n", "-------- ----------------- --------- ---------- ... ---------- ---------- -----\n", " 1 664502403.0454683 -92.63541 -30.514854 ... -0.9077294 -0.2727693 2\n", " 2 664502405.2579999 -92.64103 -28.262728 ... 1.3443842 -0.2838398 2" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# First few events (using [] for \"slicing\")\n", "events.table[:2]" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "astropy.time.core.Time" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Event times can be accessed as Astropy Time objects\n", "type(events.time)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "
OBS_IDGLON_PNTGLAT_PNTIRF
degdeg
int64float64float64bytes13
110380359.9999912037958-1.299995937905366South_z20_50h
111140358.49998338300741.3000020211954284South_z20_50h
1111591.50000565682677411.299940468335294South_z20_50h
" ], "text/plain": [ "\n", "OBS_ID GLON_PNT GLAT_PNT IRF \n", " deg deg \n", "int64 float64 float64 bytes13 \n", "------ ------------------ ------------------ -------------\n", "110380 359.9999912037958 -1.299995937905366 South_z20_50h\n", "111140 358.4999833830074 1.3000020211954284 South_z20_50h\n", "111159 1.5000056568267741 1.299940468335294 South_z20_50h" ] }, "execution_count": 53, "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": 54, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'South_z20_50h'}" ] }, "execution_count": 54, "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": 55, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'Galactic latitude (deg)')" ] }, "execution_count": 55, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAGd9JREFUeJzt3XuUZWV55/HvjxboJspFaZVb22AIiBnDpaJ4iQGjQRgFDZCAGQXFMDoyKk5wcFiJqFlLRU1meReNiiZRo0EuirYCEnQUpEGwQWxtUULTLGlRUGKLXJ75Y++CQ1F1andXnTqnq76ftc6qfXnP3s97qvs8td937/dNVSFJUhdbDDsASdLmw6QhSerMpCFJ6sykIUnqzKQhSerMpCFJ6sykIUnqzKQhSerMpCFJ6uxhww5gtu244461fPnyYYchSZuVK6+88mdVtXS6cvMuaSxfvpyVK1cOOwxJ2qwkubFLOZunJEmdmTQkSZ2ZNCRJnZk0JEmdmTQkSZ2ZNCRJnc27W25n4pzv3Mw7Vqxm3e0b2Hn7JZxyyF68YL9dhh2WJE1prr+3TBqtc75zM284exUb7r4XgJtv38Abzl4FYOKQNJKG8b1l81TrHStW3//Bj9tw9728Y8XqIUUkSf0N43vLpNFad/uGjdouScM2jO8tk0Zr5+2XbNR2SRq2YXxvmTRapxyyF0u2XPSgbUu2XMQph+w1pIgkqb9hfG/ZEd4a7zTy7ilJm4thfG+lqgZ28GEYGxsrR7mVpI2T5MqqGpuunM1TkqTOTBqSpM5MGpKkzkwakqTOTBqSpM5MGpKkzkwakqTOTBqSpM5MGpKkzkwakqTOTBqSpM5MGpKkzkwakqTOTBqSpM6GmjSSfDTJrUmunWJ/krw7yZok302y/1zHKEl6wLCvND4OPLfP/kOBPdvXicAH5iAmSdIUhpo0qupS4Od9ihwBfKIalwHbJ9lpbqKTJE007CuN6ewC3NSzvrbdJkkaglFPGplk20Pmp01yYpKVSVauX79+DsKSpIVp1JPGWmC3nvVdgXUTC1XVmVU1VlVjS5cunbPgJGmhGfWkcR7wkvYuqgOBO6rqlmEHJUkL1cOGefIknwIOAnZMshZ4I7AlQFV9ELgAOAxYA/waeOlwIpUkwZCTRlUdO83+Al41R+FIkqYx6s1TkqQRYtKQJHVm0pAkdWbSkCR1ZtKQJHVm0pAkdWbSkCR1ZtKQJHVm0pAkdWbSkCR1ZtKQJHVm0pAkdWbSkCR1ZtKQJHVm0pAkddZpPo0kOwA7AxuAn1TVfQONSpI0kqZMGkm2o5kA6VhgK2A9sBh4TJLLgPdX1dfmJEpJ0kjod6XxOeATwB9V1e29O5IcALw4yR5V9Y+DDFCSNDqmTBpV9Zw++64ErhxIRJKkkTVtn0aS/SfZfAdwY1XdM/shSZJGVZeO8PcD+wPfBQL8frv8qCSvqKqvDDA+SdII6XLL7U+A/apqrKoOAPYDrgWeDZwxwNgkSSOmS9LYu6quG1+pqu/RJJEbBheWJGkUdWmeWp3kA8Cn2/W/AH6QZGvg7oFFJkkaOV2uNI4H1gCvBU4Gbmi33Q0cPKjAJEmjZ9orjarakOT9wBeqavWE3XcOJixJ0iia9kojyeHA1cCX2/V9k5w36MAkSaOnS/PUG4EnA7cDVNXVwPIBxiRJGlFdksY9VXXHwCORJI28Lknj2iQvAhYl2TPJe4BvzsbJkzw3yeoka5KcOsn+45OsT3J1+3r5bJxXkrRpuiSN/wk8EbgL+BTwS5o7qWYkySLgfcChwD7AsUn2maToZ6pq3/b1kZmeV5K06brcPfVr4LT2NZueDKwZf0gwyaeBI4DvzfJ5JEmzpN98GucDNdX+qjp8hufeBbipZ30t8JRJyh2Z5JnAD4CTq+qmiQWSnAicCLBs2bIZhiVJmkq/5ql3Au8CfkwzY9+H29edNGNPzVQm2TYxSZ0PLK+qJwEXAmdNdqCqOrMdG2ts6dKlsxCaJGky/ebT+HeAJG+pqmf27Do/yaWzcO61wG4967sC6ybEcFvP6oeBt8/CeSVJm6hLR/jSJHuMryTZHZiNP+evAPZMsnuSrYBjgAc9NJhkp57Vw4HrZ+G8kqRN1GXAwpOBS5KMj2q7nLb/YCaq6p4kJwErgEXAR6vquiRvBlZW1XnAq9sn0u8Bfk4z5pUkaUhSNWVf9wOFmhFt925Xv19Vdw00qhkYGxurlStXDjsMSdqsJLmyqsamKzdl81SSZ4wvV9VdVXVN+7qr3b9tkt+fnXAlSZuDfs1TRyY5g2agwiuB9cBi4HdphkR/HPC/Bh6hJGlk9Lt76uQkOwBHAUcDO9Hcens98KGq+sbchChJGhV9O8Kr6hc88HyGJGmB63LLrSRJgElDkrQRTBqSpM66TPe6TZK/SfLhdn3PJM8bfGiSpFHT5UrjYzRzaTy1XV8L/N3AIpIkjawuSePxVXUGcDdAVW1g8hFqJUnzXJek8dskS2iHLU/yeJorD0nSAtNlwMI30jwVvluSfwaejgMHStKC1GW6168muQo4kKZZ6jVV9bOBRyZJGjn9pnvdf8KmW9qfy5Isq6qrBheWJGkU9bvSeFf7czEwBlxDc6XxJOBy4BlTvE+SNE9N2RFeVQdX1cHAjcD+7RzcBwD7AWvmKkBJ0ujocvfU3lW1anylqq4F9h1cSJKkUdXl7qnrk3wE+Cea227/G87VLUkLUpek8VLglcBr2vVLgQ8MLCJJ0sjqcsvtb4B/aF+SpAVs2qSR5Me0T4P3qqo9BhKRJGlkdWmeGutZXkwz9esjBxOOJGmUTXv3VFXd1vO6uar+L/CsOYhNkjRiujRP9T4ZvgXNlccjBhaRJGlkdWmeelfP8j3Aj4E/H0w4kqRR1iVpnFBVN/RuSLL7gOKRJI2wLk+Ef67jNknSPNdvlNu9gScC2yX5s55d29LcRSVJWmD6XWnsBTwP2B54fs9rf+CvZuPkSZ6bZHWSNUlOnWT/1kk+0+6/PMny2TivJGnTTHmlUVXnAucmeWpVfWu2T5xkEfA+4DnAWuCKJOdV1fd6ip0A/KKqfjfJMcDbgb+Y7VgkSd30a556fVWdAbwoybET91fVq2d47icDa8Y72ZN8GjgC6E0aRwCnt8ufA96bJFX1kCfUJUmD1+/uqfGRbFcO6Ny7ADf1rK8FnjJVmaq6J8kdwKMAp5uVpCHo1zx1frv466r6bO++JEfPwrkz2Wk3oQxJTgROBFi2bNnMI5MkTarLLbdv6LhtY60FdutZ3xVYN1WZJA8DtgN+PvFAVXVmO7Pg2NKlS2chNEnSZPr1aRwKHAbskuTdPbu2pXkyfKauAPZsHxS8GTgGeNGEMucBxwHfAo4CLrY/Q5KGp1+fxjqa/ozDgSt7tv8KOHmmJ277KE4CVgCLgI9W1XVJ3gysrKrzgH8EPplkDc0VxjEzPa8kadNluj/ck2xZVXfPUTwzNjY2VitXDqrvXpLmpyRXVtXYdOW6jD21PMlbgX3oeRLcSZgkaeHp0hH+MZo5we8BDgY+AXxykEFJkkZTl6SxpKouomnKurGqTsdJmCRpQerSPPWbJFsAP2w7rm8GHj3YsCRJo6jLlcZrgW2AVwMHAC+muQ1WkrTATHulUVVXtIt3Ai8dbDiSpFHW7+G+85lkyI5xVXX4QCKSJI2sflca75yzKCRJm4V+Axb++1wGIkkafV06wiVJAkwakqSNYNKQJHU2bdJI8tUk2/es75BkxWDDkiSNoi5XGjtW1e3jK1X1C3wiXJIWpC5J474k98+hmuRx9Hl+Q5I0f3UZe+o04BtJxm/BfSbtfNySpIWlyzAiX06yP3AgEODkqvrZwCOTJI2cKZunkuzd/twfWEYz/evNwLJ2myRpgel3pfE6mmaod02yr3BODUlacPoNIzLeb3FoVf2md1+SxZO8RZI0z3W5e+qbHbdJkua5fkOjPxbYBViSZD+aTnCAbWkmZZIkLTD9+jQOAY4HdqXp1xhPGr8E/s9gw5IkjaJ+fRpnAWclObKq/m0OY5IkjagufRoHTDL21N8NMCZJ0ojqkjQOnWTsqcMGF5IkaVR1SRqLkmw9vpJkCbB1n/KSpHmqy9hT/wRclORjNA/1vQw4a6BRSZJGUpexp85Isgr4E5o7qN5SVc6nIUkLUJcrDarqS8CXZuukSR4JfAZYDvwE+PO2r2RiuXuBVe3qf1TV4bMVgyRp43WZue/AJFckuTPJb5Pcm+SXMzzvqcBFVbUncFG7PpkNVbVv+zJhSNKQdekIfy9wLPBDYAnwcuA9MzzvETzQL3IW8IIZHk+SNAe6JA2qag2wqKruraqPAQfP8LyPqapb2mPfwtTTxy5OsjLJZUlMLJI0ZF36NH6dZCvg6iRnALcAvzPdm5JcCDx2kl2nbUR8y6pqXZI9gIuTrKqqH01yrhNpZxNctmzZxN2SpFnSJWm8GFgEnAScDOwGHDndm6rq2VPtS/LTJDtV1S1JdgJuneIY69qfNyS5BNgPeEjSqKozgTMBxsbGnL9ckgZk2uapqrqxqjZU1S+r6k1V9bq2uWomzgOOa5ePA86dWKAdrmTrdnlH4OnA92Z4XknSDPQbGn0VzcN8k6qqJ83gvG8D/jXJCcB/AEe35xwDXlFVLweeAHwoyX00ye1tVWXSkKQh6tc89bxBnbSqbqN5WHDi9pU0d2dRVd8E/sugYpAkbbx+Q6PfOJeBSJJG37Ae7pMkbYaG9XCfJGkz1HXsqTVJFlXVvcDHknxzwHFJkkbQwB7ukyTNP12ap17cljsJ+E86PtwnSZp/usynMX4X1W+ANw02HEnSKJvySiPJEUle1bN+eZIb2tdRcxOeJGmU9Gueej3NcB/jtgb+EDgIeOUAY5Ikjah+zVNbVdVNPevfaJ/kvi2JHeGStAD1u9LYoXelqk7qWV06mHAkSaOsX9K4PMlfTdyY5L8D3x5cSJKkUdWveepk4JwkLwKuarcdQNO34Sx6krQA9Ruw8FbgaUmeBTyx3fzFqrp4TiKTJI2cLs9pXAyYKCRJnZ4IlyQJMGlIkjaCSUOS1JlJQ5LUmUlDktSZSUOS1JlJQ5LUmUlDktSZSUOS1JlJQ5LUmUlDktSZSUOS1JlJQ5LUmUlDktTZUJJGkqOTXJfkviRjfco9N8nqJGuSnDqXMUqSHmpYVxrXAn8GXDpVgSSLgPcBhwL7AMcm2WduwpMkTWbaSZgGoaquB0jSr9iTgTVVdUNb9tPAEcD3Bh6gJGlSo9ynsQtwU8/62nbbQyQ5McnKJCvXr18/J8FJ0kI0sCuNJBcCj51k12lVdW6XQ0yyrSYrWFVnAmcCjI2NTVpGkjRzA0saVfXsGR5iLbBbz/quwLoZHlOSNAOj3Dx1BbBnkt2TbAUcA5w35JgkaUEb1i23L0yyFngq8MUkK9rtOye5AKCq7gFOAlYA1wP/WlXXDSNeSVJjWHdPfR74/CTb1wGH9axfAFwwh6FJkvoY5eYpSdKIMWlIkjozaUiSOjNpSJI6M2lIkjozaUiSOjNpSJI6M2lIkjozaUiSOjNpSJI6M2lIkjozaUiSOjNpSJI6M2lIkjobytDo0kJzzndu5h0rVrPu9g3svP0STjlkL16w36RT3ksjzaQhDdg537mZN5y9ig133wvAzbdv4A1nrwIwcWizY/OUNGDvWLH6/oQxbsPd9/KOFauHFJG06Uwa0oCtu33DRm2XRplJQxqwnbdfslHbpVFm0pAG7JRD9mLJlosetG3Jlos45ZC9hhSRtOnsCJcGbLyz27unNB+YNKQ58IL9djFJaF6weUqS1JlJQ5LUmUlDktSZSUOS1JlJQ5LUmUlDktRZqmrYMcyqJOuBG2d4mB2Bn81COMNmPUbPfKnLfKkHzJ+6zLQej6uqpdMVmndJYzYkWVlVY8OOY6asx+iZL3WZL/WA+VOXuaqHzVOSpM5MGpKkzkwakztz2AHMEusxeuZLXeZLPWD+1GVO6mGfhiSpM680JEmdLfikkeToJNcluS/JlHceJPlJklVJrk6yci5j7Goj6vLcJKuTrEly6lzG2EWSRyb5apIftj93mKLcve3v4+ok5811nP1M9xkn2TrJZ9r9lydZPvdRTq9DPY5Psr7n9/DyYcQ5nSQfTXJrkmun2J8k727r+d0k+891jF10qMdBSe7o+X387awHUVUL+gU8AdgLuAQY61PuJ8COw453pnUBFgE/AvYAtgKuAfYZduwTYjwDOLVdPhV4+xTl7hx2rJv6GQP/A/hgu3wM8Jlhx72J9TgeeO+wY+1Ql2cC+wPXTrH/MOBLQIADgcuHHfMm1uMg4AuDjGHBX2lU1fVVtXrYccyGjnV5MrCmqm6oqt8CnwaOGHx0G+UI4Kx2+SzgBUOMZVN0+Yx76/g54E+SZA5j7GJz+LfSSVVdCvy8T5EjgE9U4zJg+yQ7zU103XWox8At+KSxEQr4SpIrk5w47GBmYBfgpp71te22UfKYqroFoP356CnKLU6yMsllSUYpsXT5jO8vU1X3AHcAj5qT6Lrr+m/lyLZJ53NJdpub0Gbd5vD/oqunJrkmyZeSPHG2D74gZu5LciHw2El2nVZV53Y8zNOral2SRwNfTfL9NuvPqVmoy2R/zc75LXT96rERh1nW/k72AC5OsqqqfjQ7Ec5Il894JH4P0+gS4/nAp6rqriSvoLl6etbAI5t9m8Pvo4uraIYDuTPJYcA5wJ6zeYIFkTSq6tmzcIx17c9bk3ye5tJ9zpPGLNRlLdD71+CuwLoZHnOj9atHkp8m2amqbmmbCG6d4hjjv5MbklwC7EfTBj9sXT7j8TJrkzwM2I4hNztMYtp6VNVtPasfBt4+B3ENwkj8v5ipqvplz/IFSd6fZMeqmrWxtWye6iDJ7yR5xPgy8KfApHcvbAauAPZMsnuSrWg6YUfqziOaeI5rl48DHnIFlWSHJFu3yzsCTwe+N2cR9tflM+6t41HAxdX2ZI6Qaesxod3/cOD6OYxvNp0HvKS9i+pA4I7xJtLNSZLHjveNJXkyzXf8bf3ftZGGfTfAsF/AC2n+yrgL+Cmwot2+M3BBu7wHzZ0j1wDX0TQFDT32TalLu34Y8AOav8pHri40bfsXAT9sfz6y3T4GfKRdfhqwqv2drAJOGHbcE+rwkM8YeDNweLu8GPgssAb4NrDHsGPexHq8tf0/cQ3wNWDvYcc8RT0+BdwC3N3+HzkBeAXwinZ/gPe19VxFnzspR7weJ/X8Pi4DnjbbMfhEuCSpM5unJEmdmTQkSZ2ZNCRJnZk0JEmdmTQkSZ2ZNDRUSR6T5F+S3NAO0fKtJC+c5j3Lpxrls8P5jk+yc8/6R5Ls0/G9ByX5wqacd5rjvjnJs9vl1ybZZhOOcedGlk+Si5NsO8m+05P89cbG0L53aZIvb8p7tXkwaWho2oeQzgEurao9quoAmgfIdh3gaY+neW4FgKp6eVUN9aHAqvrbqrqwXX0tsNFJYxMcBlxTPU8Qz4aqWg/ckuTps3lcjQ6ThobpWcBvq+qD4xuq6saqeg/cf0Xx9SRXta+nTTxAvzJJXp9mDpRrkrwtyVE0Dwj+czvXwJIkl6SdeyTN3BFXteUv6hd4mjk/zmkH6rssyZPa7ae3cx5c0l49vbrnPX+T5Ptp5gj51Phf80k+nuSotuzOwNeSfK3dd2fP+49K8vF2eff2quyKJG+ZENsp7fbvJnnTFFX4S3qetE9yWpp5My6kGV5/fPvjk3y5vQr8epK9e7Zf1p7nzROudM5pj6/5aNhPOPpauC/g1cA/9Nm/DbC4Xd4TWNkuL6edT6BPmUOBbwLbtOvjT5VfQs/TvuPrwFKaUU537y0/IZ6DaOcqAN4DvLFdfhZwdbt8enverYEdaYZw2LI9x9XAEuARNE+7/3X7no8DR7XLP6Fn3hZ65gyhGW7k4+3yecBL2uVXjZejGeLmTJonnLcAvgA8c5K63Ag8ol0+gOYp6G2AbWmeUh+P7SJgz3b5KTTDndAe99h2+RUT4twFWDXsf1++BvNaEAMWavOQ5H3AM2iuPv6Q5sv2vUn2Be4Ffm+St01V5tnAx6rq1wBVNd1ggAfSNJP9uGP5ZwBHtmUvTvKoJNu1+75YVXcBdyW5FXhMW/7cqtrQ1vX8aY4/naePnx/4JA8MFPin7es77frDaZLpxME1H1lVv2qX/wj4/PhnlXYWxCQPpxmu5bN5YKqPrdufT+WBeU7+BXhnz7FvpacJUPOLSUPDdB0PfPFRVa9qBx8cn073ZJoxtP6A5q/m30xyjKnKhI0b2npTyk80/v67erbdS/P/bFMnWOqNaXGffb1xvbWqPjTNce9JskVV3dfnWFsAt1fVvt1CfVCcGzbyPdpM2KehYbqYZiKlV/Zs6+0E3g64pf1iezHN9KMTTVXmK8DLxu9ESvLIdvuvaJqHJvoW8MdJdp9QfiqX0rbbJzkI+Fn171T+BvD8JIvbv+D/6xTlJsb30yRPSLIFzYCU4/4fzU0D8OD+gxU09X54G9suaeaAmWg1zUCc43V5YdvH8wjg+XD/MNs/TnJ0e6wk+YP2PZfxQMI/hgf7PTbfUaA1DZOGhqaqiqaJ44+T/DjJt2km8fnfbZH3A8cluYzmi+g/JznMpGWq6ss07f4rk1wNjN9C+nHgg+Md4T2xrAdOBM5Ocg3wmWnCPx0YS/Jd4G08MMz5VHW9oo3nGuBsmqupOyYpeibwpfGOcJo50r9Ak2B7h+p+DfCqJFfQJM7x83yFprnoW0lW0UwlO1mS/CJNHw1VdRVNfa8G/g34ek+5vwROaD+T63hgutfXAq9rf2c7TajLwe3xNQ85yq00R5I8vJoZ1bah+ev+xPYLexix7EQzJ/ZzNvH92wAbqqqSHEPTKX5Eu+9S4Iiq+sXsRaxRYZ+GNHfOTPMg4WLgrGElDGjmXk/y4STbTtOsNpUDaG5ACHA78DJoHu4D/t6EMX95pSFJ6sw+DUlSZyYNSVJnJg1JUmcmDUlSZyYNSVJnJg1JUmf/H5GmALl2iX3HAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "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", "\n", "plt.scatter(\n", " Angle(table[\"GLON_PNT\"], unit=\"deg\").wrap_at(\"180 deg\"), table[\"GLAT_PNT\"]\n", ")\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=110380`." ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [], "source": [ "obs = data_store.obs(obs_id=110380)" ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Info for OBS_ID = 110380\n", "- Start time: 59235.50\n", "- Pointing pos: RA 267.68 deg / Dec -29.61 deg\n", "- Observation duration: 1800.0 s\n", "- Dead-time fraction: 2.000 %\n", "\n" ] } ], "source": [ "print(obs)" ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 58, "metadata": {}, "output_type": "execute_result" } ], "source": [ "obs.events" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Table length=3\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
EVENT_IDTIMERADECENERGYDETXDETYMC_ID
sdegdegTeVdegdeg
uint32float64float32float32float32float32float32int32
1664502403.0454683-92.63541-30.5148540.03902182-0.9077294-0.27276932
2664502405.2579999-92.64103-28.2627280.0307963711.3443842-0.28383982
3664502408.8205513-93.20372-28.5996250.040096291.0049409-0.77697752
" ], "text/plain": [ "\n", "EVENT_ID TIME RA DEC ... DETX DETY MC_ID\n", " s deg deg ... deg deg \n", " uint32 float64 float32 float32 ... float32 float32 int32\n", "-------- ----------------- --------- ---------- ... ---------- ---------- -----\n", " 1 664502403.0454683 -92.63541 -30.514854 ... -0.9077294 -0.2727693 2\n", " 2 664502405.2579999 -92.64103 -28.262728 ... 1.3443842 -0.2838398 2\n", " 3 664502408.8205513 -93.20372 -28.599625 ... 1.0049409 -0.7769775 2" ] }, "execution_count": 59, "metadata": {}, "output_type": "execute_result" } ], "source": [ "obs.events.table[:3]" ] }, { "cell_type": "code", "execution_count": 60, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 60, "metadata": {}, "output_type": "execute_result" } ], "source": [ "obs.gti" ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 61, "metadata": {}, "output_type": "execute_result" } ], "source": [ "obs.aeff" ] }, { "cell_type": "code", "execution_count": 62, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 62, "metadata": {}, "output_type": "execute_result" } ], "source": [ "obs.edisp" ] }, { "cell_type": "code", "execution_count": 63, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 63, "metadata": {}, "output_type": "execute_result" } ], "source": [ "obs.psf" ] }, { "cell_type": "code", "execution_count": 64, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 64, "metadata": {}, "output_type": "execute_result" } ], "source": [ "obs.bkg" ] }, { "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": 65, "metadata": {}, "outputs": [], "source": [ "# TODO: copy an example XML file for tutorials\n", "# This one is no good: $CTADATA/models/models_gps.xml\n", "# Too large, private in CTA, loads large diffuse model\n", "# We need to prepare something custom.\n", "# !ls $CTADATA/models/*.xml | xargs -n 1 basename" ] }, { "cell_type": "code", "execution_count": 66, "metadata": {}, "outputs": [], "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](..\/api/gammapy.spectrum.models.PowerLaw.rst) object." ] }, { "cell_type": "code", "execution_count": 67, "metadata": {}, "outputs": [], "source": [ "# Read XML file and access spectrum parameters\n", "# from gammapy.extern import xmltodict\n", "\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": 68, "metadata": {}, "outputs": [], "source": [ "# Create a spectral model the the right units\n", "# from astropy import units as u\n", "# from gammapy.spectrum.models import PowerLaw\n", "\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": 69, "metadata": {}, "outputs": [], "source": [ "# start typing here ..." ] }, { "cell_type": "markdown", "metadata": {}, "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.0" }, "nbsphinx": { "orphan": true } }, "nbformat": 4, "nbformat_minor": 2 }