Eirik Hatlen / ctools

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using ctools from Python"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "<div align=\"right\"><a href=\"../../../_downloads/ctools-from-python.ipynb\"><img src=\"../../../_static/download-notebook.jpg\" alt=\"Download Notebook\" height=\"40\"></a></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook you will learn how to use the ctools and cscripts from Python instead of typing the commands in the console.\n",
    "\n",
    "ctools provides two Python modules that allow using all tools and scripts as Python classes. To use ctools from Python you have to import the `ctools` and `cscripts` modules into Python. You should also import the `gammalib` module, as ctools without GammaLib is generally not very useful.\n",
    "\n",
    "__Warning:__ Always import gammalib before you import ctools and cscripts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import gammalib\n",
    "import ctools\n",
    "import cscripts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulating events\n",
    "\n",
    "As an example we will simulate an observation of an hour of the Crab nebula."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sim = ctools.ctobssim()\n",
    "sim['inmodel'] = '${CTOOLS}/share/models/crab.xml'\n",
    "sim['outevents'] = 'events.fits'\n",
    "sim['caldb'] = 'prod2'\n",
    "sim['irf'] = 'South_0.5h'\n",
    "sim['ra'] = 83.63\n",
    "sim['dec'] = 22.01\n",
    "sim['rad'] = 5.0\n",
    "sim['tmin'] = '2020-01-01T00:00:00'\n",
    "sim['tmax'] = '2020-01-01T01:00:00'\n",
    "sim['emin'] = 0.03\n",
    "sim['emax'] = 200.0\n",
    "sim.execute()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first line generates an instance of the ctobssim tool as a Python class. User parameters are then set using the `[ ]` operator. After setting all parameters the `execute()` method is called to execute the `ctobssim` tool. On output the `events.fits` FITS file is created. Until now everything is analogous to running the tool from the command line, but in Python you can easily combine the different blocks into more complex workflows.\n",
    "\n",
    "Remember that you can consult the manual of each tool to find out how it works and to discover all the input parameters that you can set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r\n",
      "ctobssim\r\n",
      "========\r\n",
      "\r\n",
      "Simulate event list(s).\r\n",
      "\r\n",
      "\r\n",
      "Synopsis\r\n",
      "--------\r\n",
      "\r\n",
      "This tool simulates event list(s) using the instrument characteristics \r\n",
      "specified by the instrument response function(s) and an input model. The \r\n",
      "simulation includes photon events from astrophysical sources and background\r\n",
      "events from an instrumental background model.\r\n",
      "\r\n",
      "By default, ctobssim creates a single event list. ctobssim queries a pointing\r\n",
      "direction, the radius of the simulation region, a time interval, an energy\r\n",
      "interval, an instrumental response function, and an input model. ctobssim uses\r\n",
      "a numerical random number generator for the simulations with a seed value\r\n",
      "provided by the hidden seed parameter. Changing this parameter for\r\n",
      "subsequent runs will lead to different event samples.\r\n",
      "\r\n",
      "ctobssim performs a safety check on the maximum photon rate for all model \r\n",
      "components to avoid that the tool locks up and requests huge memory \r\n",
      "resources, which may happen if a mistake was made in setting up the input \r\n",
      "model (for example if an error in the flux units is made). The maximum allowed\r\n",
      "photon rate is controlled by the hidden maxrate parameter, which by default \r\n",
      "is set to 1e6.\r\n",
      "\r\n",
      "ctobssim can also generate multiple event lists if an observation definition \r\n",
      "file is specified on input using the hidden inobs parameter. In that \r\n",
      "case, simulation information will be gathered from the file, and for each \r\n",
      "observation an event list will be created.\r\n",
      "\r\n",
      "For each event file, the simulation parameters will be written as data\r\n",
      "selection keywords to the FITS header. These keywords are mandatory for any\r\n",
      "unbinned maximum likelihood analysis of the event data.\r\n",
      "\r\n",
      "\r\n",
      "General parameters\r\n",
      "------------------\r\n",
      "\r\n",
      "inobs [file]\r\n",
      "    Input event list or observation definition XML file. If provided (i.e. the\r\n",
      "    parameter is not blank or NONE), the pointing definition and eventually the\r\n",
      "    response information will be extracted from the input file for event\r\n",
      "    simulation.\r\n",
      "\r\n",
      "inmodel [file]\r\n",
      "    Input model XML file.\r\n",
      " \t \t \r\n",
      "caldb [string]\r\n",
      "    Calibration database.\r\n",
      " \t \t \r\n",
      "irf [string]\r\n",
      "    Instrumental response function.\r\n",
      " \t \t \r\n",
      "(edisp = no) [boolean]\r\n",
      "    Apply energy dispersion?\r\n",
      " \t \t \r\n",
      "outevents [file]\r\n",
      "    Output event list or observation definition XML file.\r\n",
      " \t \t \r\n",
      "(prefix = sim_events_) [string]\r\n",
      "    Prefix for event list in observation definition XML file.\r\n",
      "\r\n",
      "(startindex = 1) [integer]\r\n",
      "    Start index of event list in observation definition XML file.\r\n",
      "\r\n",
      "(seed = 1) [integer]\r\n",
      "    Integer seed value to be used for Monte Carlo simulations. Keep this \r\n",
      "    parameter at the same value for repeatable simulations, or increment \r\n",
      "    this value for subsequent runs if non-repeatable simulations are\r\n",
      "    required.\r\n",
      " \t \t \r\n",
      "ra [real]\r\n",
      "    Right Ascension of CTA pointing (J2000, in degrees).\r\n",
      " \t \t \r\n",
      "dec [real]\r\n",
      "    Declination of CTA pointing (J2000, in degrees).\r\n",
      " \t \t \r\n",
      "rad [real]\r\n",
      "    Radius of CTA field of view (simulation cone radius) (in degrees).\r\n",
      " \t \t \r\n",
      "tmin [time]\r\n",
      "    Start time (UTC string, JD, MJD or MET in seconds).\r\n",
      " \t \t \r\n",
      "tmax [time]\r\n",
      "    Stop time (UTC string, JD, MJD or MET in seconds).\r\n",
      "\r\n",
      "emin [real]\r\n",
      "    Lower energy limit of simulated events (in TeV).\r\n",
      " \t \t \r\n",
      "emax [real]\r\n",
      "    Upper energy limit of simulated events (in TeV).\r\n",
      " \t \t \r\n",
      "(deadc = 0.95) [real]\r\n",
      "    Average deadtime correction factor.\r\n",
      "\r\n",
      "(maxrate = 1.0e6) [real]\r\n",
      "    Maximum photon rate for source models. Source models that exceed this\r\n",
      "    maximum photon rate will lead to an exception as very likely the\r\n",
      "    specified model normalisation is too large (probably due to the\r\n",
      "    a misinterpretation of units). Note that ctools specifies intensity\r\n",
      "    units per MeV.\r\n",
      "\r\n",
      "\r\n",
      "Standard parameters\r\n",
      "-------------------\r\n",
      "\r\n",
      "(publish = no) [boolean]\r\n",
      "    Specifies whether the event list(s) should be published on VO Hub.\r\n",
      "\r\n",
      "(chatter = 2) [integer]\r\n",
      "    Verbosity of the executable:\r\n",
      "     chatter = 0: no information will be logged\r\n",
      "     \r\n",
      "     chatter = 1: only errors will be logged\r\n",
      "     \r\n",
      "     chatter = 2: errors and actions will be logged\r\n",
      "     \r\n",
      "     chatter = 3: report about the task execution\r\n",
      "     \r\n",
      "     chatter = 4: detailed report about the task execution\r\n",
      " \t \t \r\n",
      "(clobber = yes) [boolean]\r\n",
      "    Specifies whether existing files should be overwritten.\r\n",
      " \t \t \r\n",
      "(debug = no) [boolean]\r\n",
      "    Enables debug mode. In debug mode the executable will dump any log file output to the console.\r\n",
      " \t \t \r\n",
      "(mode = ql) [string]\r\n",
      "    Mode of automatic parameters (default is \"ql\", i.e. \"query and learn\").\r\n",
      "\r\n",
      "(logfile = ctobssim.log) [string]\r\n",
      "    Name of log file.\r\n",
      "\r\n",
      "\r\n",
      "Related tools or scripts\r\n",
      "------------------------\r\n",
      "\r\n",
      "csobsdef\r\n"
     ]
    }
   ],
   "source": [
    "!ctobssim --help"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a Jupyter notebook a code line starting with `!` is executed in the shell, so you can do the operation above just from the command line.\n",
    "\n",
    "One of the advantages of using ctools from Python is that you can run a tool using"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sim.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main difference to the `execute()` method is that the `run()` method will not write the output (i.e., the simulated event list) to disk. Why is this useful? Well, after having typed `sim.run()` the `ctobssim` class still exists as an object in memory, including all the simulated events. You can always save to disk later using the `save()` method.\n",
    "\n",
    "The `ctobssim` class has an `obs()` method that returns an observation container that holds the simulated CTA observation with its associated events. To visualise this container, type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== GObservations ===\n",
      " Number of observations ....: 1\n",
      " Number of models ..........: 2\n",
      " Number of observed events .: 202766\n",
      " Number of predicted events : 0\n"
     ]
    }
   ],
   "source": [
    "print(sim.obs())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is one CTA observation in the container and to visualise that observation type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== GCTAObservation ===\n",
      " Name ......................: \n",
      " Identifier ................: \n",
      " Instrument ................: CTA\n",
      " Event file ................: events.fits\n",
      " Event type ................: EventList\n",
      " Statistic .................: cstat\n",
      " Ontime ....................: 3599.99999976158 s\n",
      " Livetime ..................: 3527.99999976635 s\n",
      " Deadtime correction .......: 0.98\n",
      " User energy range .........: undefined\n",
      "=== GCTAPointing ===\n",
      " Pointing direction ........: (RA,Dec)=(83.63,22.01)\n",
      "=== GCTAResponseIrf ===\n",
      " Caldb mission .............: cta\n",
      " Caldb instrument ..........: prod2\n",
      " Response name .............: South_0.5h\n",
      " Energy dispersion .........: Not used\n",
      " Save energy range .........: undefined\n",
      "=== GCTAEventList ===\n",
      " Number of events ..........: 202766 (disposed in \"events.fits\")\n",
      " Time interval .............: 58849.0008007407 - 58849.0424674074 days\n",
      " Energy interval ...........: 0.03 - 200 TeV\n",
      " Region of interest ........: RA=83.63, DEC=22.01 [0,0] Radius=5 deg\n"
     ]
    }
   ],
   "source": [
    "print(sim.obs()[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The observation contains a CTA event list that is implement by the GammaLib class `GCTAEventList`. You can access the event list using the `events()` method. To visualise the individual events you can iterate over the events using a for loop. This will show the simulated celestial coordinates `(RA, DEC)`, the coordinate in the camera system `[DETX, DETY]`, the energies and the terrestrial times (TT) of all events. Let's peek at the first events of the list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dir=RA=83.6182556152344, DEC=22.2843074798584 [0.00478756989312195,-0.000189682185722814] Energy=63.157357275486 GeV Time=315532804.699509 s (TT)\n"
     ]
    }
   ],
   "source": [
    "events = sim.obs()[0].events()\n",
    "for event in events[:1]:\n",
    "    print(event)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use this feature to inspect some of the event properties, for example look at their energy spectrum. For this we will use the Python packages matplotlib. If you do not have matplotlib you can use another plotting package of your choice or skip this step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEKCAYAAAAcgp5RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEpJJREFUeJzt3X2wXVV5x/HvryBaXwpVqCIQgwbRdGwRr2jVdqC+BTFi\nrS9knCotNaUWR6vOmI5Oa2dsB23HjlaqTYWiHQ2iRZtILFWrg61oCRgUpGjK4BDGmiiK0lot+vSP\ns6On1/ty7j3n3H2z+H5m7mSftffZ+8ma3Cf7PHudtVJVSJLa9VN9ByBJmi4TvSQ1zkQvSY0z0UtS\n40z0ktQ4E70kNc5EL0mNM9FLUuNM9JLUuEP7DgDgyCOPrLVr1/YdhiQdVK655pqvV9VRix23KhL9\n2rVr2bVrV99hSNJBJclXRjmu19JNko1Jtt5xxx19hiFJTes10VfVjqrafPjhh/cZhiQ1zYexktQ4\nE70kNc5EL0mNM9FLUuNM9JLUOBO9JDVuVXxhahxrt1ze27VvOf+M3q4tSaPyjl6SGjfxRJ/k1CSf\nSvKOJKdO+vySpKUZKdEnuSjJviTXz2rfkOSmJHuSbOmaC7gTuBewd7LhSpKWatQ7+ouBDcMNSQ4B\nLgBOB9YDm5KsBz5VVacDrwH+eHKhSpKWY6REX1VXArfPaj4F2FNVN1fV94FLgDOr6ofd/m8C95xY\npJKkZRln1M0xwK1Dr/cCj0vyHODpwBHA2+Z7c5LNwGaANWvWjBGGJGkhEx9eWVWXAZeNcNxWYCvA\nzMxMTToOSdLAOKNubgOOG3p9bNc2Muejl6TpGyfRXw2ckOT4JIcBZwHbl3IC56OXpOkbdXjlNuAq\n4MQke5OcU1V3AecBVwA3ApdW1Q3TC1WStBwj1eiratM87TuBncu9eJKNwMZ169Yt9xSSpEW4lKAk\nNc65biSpcb0mekfdSNL0WbqRpMZZupGkxlm6kaTGWbqRpMZZupGkxpnoJalx1uglqXHW6CWpcZZu\nJKlxJnpJatzEV5haioN99sq1Wy7v5bq3nH9GL9eVdHCyRi9JjbN0I0mNM9FLUuNM9JLUOBO9JDXO\nb8ZKUuMcdSNJjbN0I0mNM9FLUuNM9JLUOBO9JDXORC9JjTPRS1LjTPSS1Di/MCVJjfMLU5LUOEs3\nktQ4E70kNc5EL0mNM9FLUuN6XRxcy9PXouTgwuTSwcg7eklqnIlekhpnopekxpnoJalxU0n0Se6T\nZFeSZ07j/JKk0Y2U6JNclGRfkutntW9IclOSPUm2DO16DXDpJAOVJC3PqHf0FwMbhhuSHAJcAJwO\nrAc2JVmf5KnAF4F9E4xTkrRMI42jr6ork6yd1XwKsKeqbgZIcglwJnBf4D4Mkv93k+ysqh9OLGJJ\n0pKM84WpY4Bbh17vBR5XVecBJDkb+Pp8ST7JZmAzwJo1a8YIQ5K0kKmNuqmqi6vqwwvs31pVM1U1\nc9RRR00rDEm62xsn0d8GHDf0+tiubWQuPCJJ0zdOor8aOCHJ8UkOA84Cti/lBC48IknTN1KNPsk2\n4FTgyCR7gT+qqguTnAdcARwCXFRVNyzl4kk2AhvXrVu3tKjVm74mVHMyNWn5Rh11s2me9p3AzuVe\nvKp2ADtmZmZestxzSJIW5hQIktS4XhO9D2Mlafp6TfQ+jJWk6bN0I0mNs3QjSY2zdCNJjbN0I0mN\nM9FLUuOs0UtS46zRS1LjLN1IUuNM9JLUOBO9JDXOh7GS1DgfxkpS4yzdSFLjRlp4ROqbK1tJy+cd\nvSQ1zkQvSY1z1I0kNc5RN5LUOEs3ktQ4E70kNc5EL0mNM9FLUuNM9JLUOBO9JDXOcfSS1Lhe57qp\nqh3AjpmZmZf0GYc0n77m2AHn2dHkWLqRpMaZ6CWpcSZ6SWqciV6SGmeil6TGmeglqXEmeklqnIle\nkhpnopekxk080Sd5ZJJ3JPlAkt+d9PklSUszUqJPclGSfUmun9W+IclNSfYk2QJQVTdW1bnA84En\nTj5kSdJSjHpHfzGwYbghySHABcDpwHpgU5L13b5nAZcDOycWqSRpWUZK9FV1JXD7rOZTgD1VdXNV\nfR+4BDizO357VZ0OvHCSwUqSlm6c2SuPAW4der0XeFySU4HnAPdkgTv6JJuBzQBr1qwZIwxJ0kIm\nPk1xVX0S+OQIx20FtgLMzMzUpOOQDnZ9TZHs9MjtGWfUzW3AcUOvj+3aRubCI5I0feMk+quBE5Ic\nn+Qw4Cxg+1JOUFU7qmrz4YcfPkYYkqSFjDq8chtwFXBikr1Jzqmqu4DzgCuAG4FLq+qG6YUqSVqO\nkWr0VbVpnvadjDGEMslGYOO6deuWewpJ0iJ6nQLB0o0kTZ9z3UhS43pN9I66kaTps3QjSY2zdCNJ\njbN0I0mNs3QjSY2zdCNJjTPRS1LjrNFLUuOs0UtS4yzdSFLjTPSS1Dhr9JLUOGv0ktQ4SzeS1DgT\nvSQ1zkQvSY0z0UtS4xx1I0mNc9SNJDXO0o0kNc5EL0mNM9FLUuMO7TsASavL2i2X93btW84/o7dr\nt8w7eklqnIlekhpnopekxvmFKUlqnF+YkqTGWbqRpMaZ6CWpcSZ6SWqciV6SGmeil6TGmeglqXEm\neklqnIlekhpnopekxk1lmuIkzwbOAH4GuLCq/mka15EkLW7kO/okFyXZl+T6We0bktyUZE+SLQBV\n9aGqeglwLvCCyYYsSVqKpZRuLgY2DDckOQS4ADgdWA9sSrJ+6JDXdfslST0ZOdFX1ZXA7bOaTwH2\nVNXNVfV94BLgzAy8EfhIVV071/mSbE6yK8mu/fv3Lzd+SdIixn0Yewxw69DrvV3by4CnAM9Ncu5c\nb6yqrVU1U1UzRx111JhhSJLmM5WHsVX1VuCt0zi3JGlpxr2jvw04buj1sV3bSFx4RJKmb9xEfzVw\nQpLjkxwGnAVsH/XNLjwiSdO3lOGV24CrgBOT7E1yTlXdBZwHXAHcCFxaVTcs4Zze0UvSlI1co6+q\nTfO07wR2LufiVbUD2DEzM/OS5bxfkrQ4p0CQpMb1mugt3UjS9PWa6H0YK0nTZ+lGkhpn6UaSGmfp\nRpIaZ+lGkhpnopekxlmjl6TGWaOXpMZZupGkxpnoJalxJnpJapwPYyWpcT6MlaTGWbqRpMZNZXFw\nSVqOtVsu7+W6t5x/Ri/XXSne0UtS40z0ktQ4R91IUuMcdSNJjbN0I0mNM9FLUuNM9JLUOBO9JDXO\nRC9JjTPRS1LjHEcvSY1zHL0kNc7SjSQ1zkQvSY0z0UtS40z0ktQ4E70kNc5EL0mNM9FLUuNM9JLU\nOBO9JDVu4ok+yUOTXJjkA5M+tyRp6UZK9EkuSrIvyfWz2jckuSnJniRbAKrq5qo6ZxrBSpKW7tAR\nj7sYeBvw7gMNSQ4BLgCeCuwFrk6yvaq+OOkgJWma1m65vLdr33L+GVO/xkh39FV1JXD7rOZTgD3d\nHfz3gUuAMyccnyRpTOPU6I8Bbh16vRc4JskDkrwDeHSSP5jvzUk2J9mVZNf+/fvHCEOStJBRSzcj\nq6pvAOeOcNxWYCvAzMxMTToOSdLAOHf0twHHDb0+tmsbmQuPSNL0jZPorwZOSHJ8ksOAs4DtSzmB\nC49I0vSNOrxyG3AVcGKSvUnOqaq7gPOAK4AbgUur6obphSpJWo6RavRVtWme9p3AzuVePMlGYOO6\ndeuWewpJ0iJcM1aSGtdrovdhrCRNX6r6H9mYZD/wlTFPcyTw9QmEM23GOVnGOVnGOVnTjvMhVXXU\nYgetikQ/CUl2VdVM33Esxjgnyzgnyzgna7XE6TTFktQ4E70kNa6lRL+17wBGZJyTZZyTZZyTtSri\nbKZGL0maW0t39JKkORy0iT7J85LckOSHSeZ9qp3kliRfSLI7ya6VjLG7/qhx/sRqXSspyf2TfDTJ\nl7s/f3ae43rpz8X6JwNv7fZ/PsnJKxXbEmI8NckdXd/tTvKHKx1jF8ecK8YN7e+9L7s4Fouz9/5M\nclySTyT5Yvd7/vI5jum/P6vqoPwBHgmcCHwSmFnguFuAI1dznMAhwH8ADwUOA64D1q9wnG8CtnTb\nW4A3rpb+HKV/gGcAHwECPB747CqM8VTgw339WxyK41eAk4Hr59nfa18uIc7e+xM4Gji5274f8KXV\n9m+zqg7eO/qqurGqbuo7jsWMGOdqWK3rTOBd3fa7gGev8PUXMkr/nAm8uwY+AxyR5OhVFuOqUHOv\nGDes774ERoqzd1X11aq6ttv+DoMJHo+ZdVjv/XnQJvolKOBjSa5JsrnvYOYx52pdKxzDA6vqq932\nfwIPnOe4PvpzlP7puw9Hvf4Tuo/vH0ny8ysT2pL13ZdLsWr6M8la4NHAZ2ft6r0/J77C1CQl+Rjw\noDl2vbaq/mHE0zypqm5L8nPAR5P8e3enMDETinPqFopz+EVVVZL5hmNNvT8bdi2wpqruTPIM4EPA\nCT3HdDBbNf2Z5L7A3wOvqKpv9xHDQlZ1oq+qp0zgHLd1f+5L8kEGH7EnmpgmEOfYq3WNYqE4k3wt\nydFV9dXuY+W+ec4x9f6cwyj9syJ9uIBFrz+cAKpqZ5K/SnJkVa22OVv67suRrJb+THIPBkn+PVV1\n2RyH9N6fTZduktwnyf0ObANPA+Z8gt+zsVfrmoDtwIu77RcDP/FJpMf+HKV/tgMv6kY4PB64Y6gU\ntRIWjTHJg5Kk2z6Fwe/fN1YwxlH13ZcjWQ392V3/QuDGqnrzPIf13599PrEe5wf4NQa1ru8BXwOu\n6NofDOzsth/KYPTDdcANDEopqy7O+vGT+S8xGLnRR5wPAD4OfBn4GHD/1dSfc/UPg0Xoz+22A1zQ\n7f8CC4zE6jHG87p+uw74DPCElY6xi2Mb8FXgf7t/m+estr4cMc7e+xN4EoPnVp8Hdnc/z1ht/ek3\nYyWpcU2XbiRJJnpJap6JXpIaZ6KXpMaZ6CWpcSZ6TUWSOyd0nn9M8q0kH57VfnySz3YzAr6vG7t+\nYN+zD8xkmOT1SW4bmuFwd5IjJhHbOJKcnWR/kncmefpQbHdmMAPm7iTvnue9SfKVJA+b1f62JK9K\nclo3m+Lurv1H/aG7J4dXaiqS3FlV953AeZ4M3Bv4nap65lD7pcBlVXVJkncA11XV27t9nwaeVVVf\nT/J64M6q+vNxY5kV1yFV9YMx3n82g/HU581q/yTw6qpacAroJG9i8MWbPzkQD4P5VE6pqr1J1gEf\nqKqTui/1fA54fFX9z3Jj1sHLO3qtmCRrk/xzNwnVx5Os6dofluQzGcxz/4bhTwNV9XHgO7POE+BX\ngQ90TT+abTPJw4Hv1SJfg+/uqC/rPjF8uUucB/Y9LclVSa5N8v4M5jE5MBf/G5NcCzwvyWO7v8vu\nJH+Wbt70JFcmOWnofP+S5BeX2WeHJnlzkn/rrvXb3a5tDL59e8BpwJerau/sc9Tgbu5TDL7Io7sh\nE71W0l8C76qqXwDeA7y1a38L8JaqehSDb0Au5gHAt6rqru718GyAT2Qw2dWw3x8qjXxiqP0k4AXA\no4AXZLCIxJHA64CnVNXJwC7glUPv+UZVnVxVlwB/y+CTxknA8N39hcDZ8KP/eO5VVdeN8Peay2Zg\nX1WdAjwW+L0ka6rqc8A9kqzvjjuLQfKfzy7gl5cZgw5yJnqtpF8C3ttt/x2Dr48faH9/t/3e2W9a\noqOB/bPa/qKqTup+Thtq/3hV3dGVM74IPITBwhDrgX/tatwv7toPeB9AV+e/X1VdNUfc7weemcFk\nV78FXDzG3+dpwG92sXwWOIIfz9C4DdjUXWcjP/6EM5d9DKaz0N3Qqp69UprHNxgs3nBod1c/PBvg\nd4HDRzzP94a2f8Dg9yHAR6tq0zzv+a/FTlpV/53kowwWnHg+8JgR45lLgJd2JazZtgEfZvAfwDWL\nlKvuxaBvdDfkHb1W0qf5cV35hQzqxjCYkOrXu+2zZr9ptq7m/AnguV3T8GybNwLrxojxM8ATu4eZ\nB2bsfPgcMXwL+E6Sx80T9zsZlKaurqpvjhHPFcBLkxzaxXNikp/uYvgSg+cXb2Dhsg3Aw1mdM7dq\nBZjoNS33TrJ36OeVwMsYlCE+D/wGcGAh5VcAr+za1wF3HDhJkk8xKIU8uTvP07tdr+nes4dBzf7C\nrv1K4NHdA9sDhmv0uzNYCWhOVbWfQX19WxfPVcAj5jn8HOBvurLKfYbjrqprgG8zqOOP468ZzCi6\nu3vY+3b+/yfxbQzWJP7QIuc5Dbh8zFh0kHJ4pXqX5N7Ad6uqkpwFbKqqZa+3muQtwI6q+tjEgpz7\nOvetqju77S3A0VX18u71gxksCP+IqvrhHO89mzmGV04wtuHhlQ8GLq6qp03jWlr9vKPXavAYBnes\nnwdeCrxqzPP9KYOx99N2RvcJ4XoGI1reAJDkRQzq5q+dK8l3vgucnuSdkw4qyWnAB4EDNfvjgFdP\n+jo6eHhHL0mN845ekhpnopekxpnoJalxJnpJapyJXpIaZ6KXpMb9H8jQkiomH+TPAAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x114dbb128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "#this will visualize plots inline\n",
    "\n",
    "ax = plt.subplot()\n",
    "ax.set_yscale('log')\n",
    "ax.set_xlabel('Log10(Energy [TeV])')\n",
    "\n",
    "energies = []\n",
    "for event in events:\n",
    "    energies.append(event.energy().log10TeV())\n",
    "    \n",
    "n, bins, patches = plt.hist(energies)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting a model to the observations\n",
    "\n",
    "We can use the observation in memory to directly run a likelihood fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "like = ctools.ctlike(sim.obs())\n",
    "like.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is very compact. Where do we define the model fit to the data? Where are the user parameters? `ctlike` doesn’t in fact need any parameters as all the relevant information is already contained in the observation container produced by the ctobssim class. Indeed, you constructed the ctlike instance by using the ctobssim observation container as constructor argument.\n",
    "\n",
    "An observation container, implemented by the `GObservations` class of GammaLib, is the fundamental brick of any ctools analysis. An observation container can hold more than events, e.g., in this case it also holds the model that was used to generate the events.\n",
    "\n",
    "Many tools and scripts handle observation containers, and accept them upon construction and return them after running the tool via the `obs()` method. Passing observation containers between ctools classes is a very convenient and powerful way of building in-memory analysis pipelines. However, this implies that you need some computing ressources when dealing with large observation containers (for example if you want to analyse a few 100 hours of data at once). Also, if the script crashes the information is lost.\n",
    "\n",
    "To check how the fit went you can inspect the optimiser used by `ctlike` by typing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== GOptimizerLM ===\n",
      " Optimized function value ..: 757618.425\n",
      " Absolute precision ........: 0.005\n",
      " Acceptable value decrease .: 2\n",
      " Optimization status .......: converged\n",
      " Number of parameters ......: 10\n",
      " Number of free parameters .: 4\n",
      " Number of iterations ......: 2\n",
      " Lambda ....................: 1e-05\n"
     ]
    }
   ],
   "source": [
    "print(like.opt())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You see that the fit converged after 2 iterations. Out of 10 parameters in the model 4 have been fitted (the others were kept fixed). To inspect the fit results you can print the model container that can be access using the `models()` method of the observation container:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== GModels ===\n",
      " Number of models ..........: 2\n",
      " Number of parameters ......: 10\n",
      "=== GModelSky ===\n",
      " Name ......................: Crab\n",
      " Instruments ...............: all\n",
      " Instrument scale factors ..: unity\n",
      " Observation identifiers ...: all\n",
      " Model type ................: PointSource\n",
      " Model components ..........: \"PointSource\" * \"PowerLaw\" * \"Constant\"\n",
      " Number of parameters ......: 6\n",
      " Number of spatial par's ...: 2\n",
      "  RA .......................: 83.6331 [-360,360] deg (fixed,scale=1)\n",
      "  DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1)\n",
      " Number of spectral par's ..: 3\n",
      "  Prefactor ................: 5.73901839854211e-16 +/- 5.75150678022232e-18 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)\n",
      "  Index ....................: -2.46442704108736 +/- 0.00738985717701118 [-0,-5]  (free,scale=-1,gradient)\n",
      "  PivotEnergy ..............: 300000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)\n",
      " Number of temporal par's ..: 1\n",
      "  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)\n",
      "=== GCTAModelIrfBackground ===\n",
      " Name ......................: CTABackgroundModel\n",
      " Instruments ...............: CTA\n",
      " Instrument scale factors ..: unity\n",
      " Observation identifiers ...: all\n",
      " Model type ................: \"PowerLaw\" * \"Constant\"\n",
      " Number of parameters ......: 4\n",
      " Number of spectral par's ..: 3\n",
      "  Prefactor ................: 0.995734795709415 +/- 0.00725726198557001 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient)\n",
      "  Index ....................: -0.00198344390952578 +/- 0.00257263793207067 [-5,5]  (free,scale=1,gradient)\n",
      "  PivotEnergy ..............: 1000000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)\n",
      " Number of temporal par's ..: 1\n",
      "  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)\n"
     ]
    }
   ],
   "source": [
    "print(like.obs().models())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, in this way we can fetch the minimum of the optimized function (the opposite of the likelihood logarithm) to compare different model hypotheses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "757618.4249159902\n"
     ]
    }
   ],
   "source": [
    "like1 = like.opt().value()\n",
    "print(like1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose you want to repeat the fit by optimising also the position of the point source. This is easy from Python, as you can modify the model and fit interactively. Type the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== GModels ===\n",
      " Number of models ..........: 2\n",
      " Number of parameters ......: 10\n",
      "=== GModelSky ===\n",
      " Name ......................: Crab\n",
      " Instruments ...............: all\n",
      " Instrument scale factors ..: unity\n",
      " Observation identifiers ...: all\n",
      " Model type ................: PointSource\n",
      " Model components ..........: \"PointSource\" * \"PowerLaw\" * \"Constant\"\n",
      " Number of parameters ......: 6\n",
      " Number of spatial par's ...: 2\n",
      "  RA .......................: 83.6334964430607 +/- 0.000614768905318078 [-360,360] deg (free,scale=1)\n",
      "  DEC ......................: 22.0149467750683 +/- 0.000559548446961086 [-90,90] deg (free,scale=1)\n",
      " Number of spectral par's ..: 3\n",
      "  Prefactor ................: 5.73898570353797e-16 +/- 5.75155055360256e-18 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)\n",
      "  Index ....................: -2.46439956654049 +/- 0.00739170138357582 [-0,-5]  (free,scale=-1,gradient)\n",
      "  PivotEnergy ..............: 300000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)\n",
      " Number of temporal par's ..: 1\n",
      "  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)\n",
      "=== GCTAModelIrfBackground ===\n",
      " Name ......................: CTABackgroundModel\n",
      " Instruments ...............: CTA\n",
      " Instrument scale factors ..: unity\n",
      " Observation identifiers ...: all\n",
      " Model type ................: \"PowerLaw\" * \"Constant\"\n",
      " Number of parameters ......: 4\n",
      " Number of spectral par's ..: 3\n",
      "  Prefactor ................: 0.995729516865969 +/- 0.0072572344873009 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient)\n",
      "  Index ....................: -0.00198584806013523 +/- 0.00257264072051407 [-5,5]  (free,scale=1,gradient)\n",
      "  PivotEnergy ..............: 1000000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)\n",
      " Number of temporal par's ..: 1\n",
      "  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)\n"
     ]
    }
   ],
   "source": [
    "like.obs().models()['Crab']['RA'].free()\n",
    "like.obs().models()['Crab']['DEC'].free()\n",
    "like.run()\n",
    "print(like.obs().models())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can quantify the improvement of the model by comparing the new value of the optimized function with the previous one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0637600922491401\n"
     ]
    }
   ],
   "source": [
    "like2 = like.opt().value()\n",
    "ts = -2 * (like2 - like1)\n",
    "print(ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The test statistic TS is expected to be distributed as a chi square with a number of degrees of freedom equal to the additional number of degrees of freedom of the (second) more complex model with respect to the (first) more parsimonious one, for our case two degrees of freedom. The chance probability that the likelihood improved that much because of pure statistical fluctuations is the integral from TS to infinity of the chi square with two degrees of freedom. In this case the improvement is negligible (i.e., the chance probability is very high), as expected since in the model the source is already at the true position."
   ]
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}