{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to do advanced model manipulation and fitting in Python?"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "<div align=\"right\"><a href=\"../../../_downloads/howto_advanced_python_fitting.ipynb\" download=\"howto_advanced_python_fitting.ipynb\"><img src=\"../../../_static/download-notebook.jpg\" alt=\"Download Notebook\" height=\"40\"></a></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial you will learn to perform some more advanced model manipulation and fitting from Python. You can use these methods in an interactive Python session to explore your dataset and improve the known model of the sky. You can also use them in a Python pipeline to test multiple model hypotheses and select the one that is most appropriate.\n",
    "\n",
    "As usual start by importing the gammalib, ctools, and cscripts Python modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gammalib\n",
    "import ctools\n",
    "import cscripts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now import the matplotlib package for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulated dataset\n",
    "\n",
    "For the tutorial you will simulate a small dataset. Start by defining the Instrument Response Function and energy range that will be used all along."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "caldb = 'prod2'\n",
    "irf   = 'South_0.5h'\n",
    "emin  = 0.1   # TeV\n",
    "emax  = 160.0 # TeV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now proceed to simulate the dataset. It consists of an hour of observations of the Crab nebula region, as usual pointed at a slightly offset position from the target. The input model is different from the one you have been using so far and contains some surprises. Don't look at it until you have completed the exercises at the end of the tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "evfile = 'events.fits'\n",
    "obssim = ctools.ctobssim()\n",
    "obssim['ra']        = 83.63\n",
    "obssim['dec']       = 22.51\n",
    "obssim['rad']       = 5.0\n",
    "obssim['tmin']      = 0\n",
    "obssim['tmax']      = 3600\n",
    "obssim['emin']      = emin\n",
    "obssim['emax']      = emax\n",
    "obssim['caldb']     = caldb\n",
    "obssim['irf']       = irf\n",
    "obssim['inmodel']   = '$CTOOLS/share/models/crab_beyond.xml'\n",
    "obssim['outevents'] = evfile\n",
    "obssim.execute()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You have saved the events on disk in the file `events.fits`. In this way you can easily re-use the code below by substituting `events.fits` with your own event list or observation definition XML file.\n",
    "\n",
    "Let's first peek at the simulated dataset by displaying a skymap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "skymap = ctools.ctskymap()\n",
    "skymap['inobs']       = evfile\n",
    "skymap['emin']        = emin\n",
    "skymap['emax']        = emax\n",
    "skymap['nxpix']       = 40\n",
    "skymap['nypix']       = 40\n",
    "skymap['binsz']       = 0.02\n",
    "skymap['proj']        = 'TAN'\n",
    "skymap['coordsys']    = 'CEL'\n",
    "skymap['xref']        = 83.63\n",
    "skymap['yref']        = 22.01\n",
    "skymap['bkgsubtract'] = 'IRF'\n",
    "skymap['caldb']       = caldb\n",
    "skymap['irf']         = irf\n",
    "skymap.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVYAAAEKCAYAAABJ430PAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJztvXmUZEd95/v55lJ7dbe6JZDUaiyxeeAwgMYyxpY9YARG1mCwMXjABvOMHxqPwYNsPGMEcx7H4+d5g+Hosc5gGbHZeixGgsMiI9q2gNEckLUgtLhhEGCgUVut3qqX2nL5vT9uNmRV/KLqVnVmdmX176OTR51RceNG3HszMjK+v0VmRhAEQdA7Kqe7A0EQBJuNmFiDIAh6TEysQRAEPSYm1iAIgh4TE2sQBEGPiYk1CIKgx8TEGgRB0GNiYg2CIOgxfZtYJe2SdIukPZLul/TaTvlbJH1d0j2SPiFp2wptVCV9VdJn+tXPIAiCXqN+eV5JOg84z8zukjQN3An8MnAB8Pdm1pT0ZgAz+6NMG38AXAJsMbPnrXbO+uikjU5uX1qYGZ5a6R+8slwjWsNlM2X+oNwfljdQ/mRW9b8r2/X0XGr77VbnW2lho1m+b5Xc93Va11rttF+ZMVCpZtp1zlRz2qj419uc+2A1v673jFQa6RjyHXOugddXwKpOH3rxcV1LG8u6MD93mMbiiZIPrs9zf37SDh5ynjGHO+9ZuNnMLj+V83Uj6ZeBfwM8Ani3mX2+V213U+tHowBmtg/Y1/n3MUl7gJ3LBvIV4EXe8ZIuoLgAfwr8QZlzjk5u50nPvWpJWXXRf4pGZtKJon5s0W+46Xz425kPk/fBGfEvc3sknSjkncuZfIoG0qLFHWNu1dlz60nZyDG/3alvHkkLH3zIretOjBPjbl2a6TVvHz+RlFWmJt3DNT3lFPqf8db2tG5rPL0GAK2x9D4sbvPvWf1YOiGMPTSbVnTuI4Cc56O5zb9ei9vS/io3HzntZr/8nS/UXN3WyNJJ/6u3viPTaHkOHmrxDzc/qlTd6nnfPHu1OpLeBzwP2G9mT+oqvxx4O1AF3mtm/83MPgl8UtJZwFuBvkysA9ljlXQhcDFw27I/vRL4m8xhbwP+E+70EQTBsGJAu+R/JfkAsGRVK6kKvBv4ReCJwEslPbGryn/u/L0v9H1ilTQF3ABcZWZHu8rfCDSB651jTn773Fmi/Ssl3SHpjsZCuvoJgmBjYRgNa5V6AWef/Hx3Xlcm7Zl9CTi0rPhpwANm9m0zWwQ+ArxABW8G/sbM7urXGPu2FQAgqU4xqV5vZjd2lb+CYul+mfmbvJcCz5d0BTAGbJH0V2b2suUVzexa4FqAqe27IlRXEAwBa1iNHjCzS9Zxip3A97ve7wV+Cvg94NnAVkmPNbP3rKPtVenbxCpJwHXAHjO7pqv8cuCPgGeYmbM5BWZ2NXB1p/4zgT/0JlXnQGrzS29Ybc6/gbXj6X6qt+cJsHhOum+ZE7rqxxppu6N+u01nb6+6kPa3OucLRxVnH68252/Cjc6k5/IELYDWlnS81f3+o6J2Ol6Njbp1aY8kRRVPFMsIZa196T5vZdQ/l7ZMpKfP3IfGlnRszVH/2tSPucUJNuZfL++pyQll3tzj7cEDrjCXE01tdA0/VPuwEWcYrf6HK/VGb2b2DuDUN4pXoZ9bAZcCLweeJenuzusK4F3ANLC7U/YeAEnnS7qpj/0JgmCD0MZKvSixFZBhL7Cr6/0FwIO9HkeOfloF3Ir/reFOnmb2IHCFU/4F4Au97FsQBKcPA1rlbb7WuxVwO/A4SRcBPwBeAvz6OtpZF+F5FQTBwFnDinVVJH0Y+DLw45L2SvptM2sCrwFuBvYAHzOz+/s2oGX0VbwKgiBYjgGNHu6xmtlLM+U3kfmF3G9iYg2CYKAYtpatgLMl3dH1/tqOJdCGZlNNrJVW6lG1FnfD5rivGs/tSC9Tpek/GJ4HS06d9ep6FgC1o/Pu8W3Poytj2VCbTa0FWhl1uF13rBW2+yEd1HTa3ep7Tnl4vbWZo04prsdQziW7MpPaNI9mFPXqXGoFMZqxIKg47r6eR1drwv9oVRadPmQe0frx9FmoOtYsAFZP72Vrwvc0w2ui5Dy3FlfuLAZZ7/GU9e6xnlY21cQaBMHGp/C82tzExBoEwYARLddgaPMQE2sQBAOlEK9KT6yxxxoEQbAahR1r6Yk19lhPNwZYJubmcuSIWm4cUvy4pVl3UEcQysU9rZ1Iz1c7lIouOjHn92vHlqSsOZlxPXXEtpGZ1B0VcEPQNR6ZngugOZGKPK0xXxSrNNJ2x5yYo9VMPNaq5+pay8RonV9Iy44dd6vWZ5ywfZlwhDaaCkLts9IQhZYTymZT5SjnppqL01qWyoL/PHuu0M1JX+haLuiWX2iuTLtXDW1QNtXEGgTBxmeNK9ahJCbWIAgGiiFa5Z0+Y481CIKgDGvYCog91iAIgtUwxKKVz182jMTEGgTBQCkcBDZ3/KdNNbG2RsXMo5cGUx4/6Cuu0/87Vdqrh/0oxpOj5yRlc2f7KqqbDTWTAC6fFXYpNuEnCGxOp0GemxP+A1ubTa9DJWMV4KnUizvSwNEAi1uclUcu+6zSdufOc9o91z+XZ13hueoCjPzASYj48EG3rh1N77u1MhYi4+m9qIykz0LWjdlR6rWQSWI5mZ6rsc1/Fjyq837A8LZjbbCw3X+eZ89ZWreVCQC+VkK8CoIg6CFmomWxYg2CIOgp7fIr1rAKCIIgWI1CvCo99YRVQBAEwWqEeHUKSNoFfAg4lyJK2LVm9nZJbwF+iSIq5LeA3zKzI8uOHQO+BIx2+vhxM3vTaue0CjQml/7EUNu/gZNOLNPKAT/u6ciBNJlsu566MZ7sQ9KuF4MT373Ro73VF3MaW1PBoTHuj9ccQWnE13KoHErFnFws08pi6g6ac+FtjafX/Oij0syti1v8n4nVhbTdiQMZse5EGhO2uphx4fWoZ7KsOhloW46ImMsIW1l0xCsn7ioAc+nzUcnEWPXiseaQ47Jcnc88owtL23X0x3XR2uQurf382mgCrzOzJwBPB14t6YnAbuBJZvZk4H/TSXO9jAXgWWb2FOCpwOWSnt7HvgZBMCBOel6VeQ0r/czSug/Y1/n3MUl7gJ1m9vmual8BXuQca8DJiBn1zqvviciDIBgM7U1uFTCQ0Um6ELgYuG3Zn14J/E3mmKqku4H9wG4zW37syXpXnsw53ppLI0MFQbCxKIKwlF6xnn3y8915XXmau1+KvotXkqaAG4CrzOxoV/kbKbYLrveOM7MW8FRJ24BPSHqSmd3n1LsWuBZg/JG7YlUbBBscQzTKu7SGVcByJNUpJtXrzezGrvJXAM8DLrNcNrgOZnZE0heAy4FkYg2CYLgwIxwE1oskAdcBe8zsmq7yy4E/Ap5hZqncXtQ5B2h0JtVx4NnAm1c7Z23BOOubS5XfZkYlbzvqamWbH8xZjVTJHTmUyZw6ml5StTKK+nEnGPOhmbRe5rtH21M1OqfaekGm3fMD7SNpHzTnB9semXFcLHOK+gXb07I1fL7qztMyethX1DWXWgB4ij5A4xHTSdniNl9991w6vYy9teMZl1hH6W8fPOzXHU0tJqrjaRlAcyotr2RcWr3nedy5XgBjDy+9l9W5XpgFaC0OAkNJP1eslwIvB+7t7JUCvAF4B4UZ1e5i7uUrZvY7ks4H3mtmVwDnAR+UVKXYB/6YmX2mj30NgmBAGLFiXTdmdiu4X0s3Zeo/CFzR+fc9FGJXEASbkGE2pSpDeF4FQTBQDEXOqyAIgl5SpL/e3FPPphpdZaHFxLeWCgELF2x163rZTCsLfqzLypHUPrYym3EddWJdem6MAHJcLG3BE7Sc2KLA2EQqxlSaTsZRoDqbChma9QW4bDxVD09YyxxfcUSTyYfSssaJjJuq43bpjQtATjzV9pQvXs2e67jVTvl9qM+l460fT/tVy/XLc6vNZKXVllRUmz0/ddUFPzvwZCZLqxbSPmjO72/lwFIhU16m3DWjtcRjjehWQRAEq2GsyfMq7FiDIAjKEBkEgiAIeoiZNn2sgJhYgyAYKIV4FVlagyAIekjkvBo+Ksv2bjKRCNrVdI/HU/SLNtNyL1gw+BlOlXEt9BR1TTsBtHMq/WLabu1E+cyrNuq7bWrSCaxdy7ipbklV6vaYX7d6NLVCmNp/NClrbfWV78b21GqjPeavfNpTad1WJvi056ZadVyAAaa+m7r21v85tdqwnFuvV37+I9268+emVgFzO/x2a04Q8NaEX7ddS+9vNePSmlgxrMViJEMhXsUeaxAEQU8Jz6sgCIIeEp5XQRAEfSCSCQZBEPQQM2hkknxuFjbVxGrVCs2zlm7MVxq57JNOtsxWRulyBIecO2jVc/mr+aKJOfE2zRFurO4f78XgbI77db3MqfURv25txumXI/YBLJ6T9rc15n9oxn/gpM7Zn6aKzRnitLak/WplspO2nbHlsseOzqTPSLuWyRR7InU5Ni+e6khGGHTcVBs7/Iy/jel0DGNHfDfVkSOOm2oms64nauWE28oyIdL2n7qZVLEVEBNrEARBT4lYAUEQBD1kjeZWESsgCIJgdWIrIAiCoOds9pxXffvakLRL0i2S9ki6X9JrO+VvkfR1SfdI+kQnvXWpY4MgGH4Kq4Bqqdew0s8VaxN4nZndJWkauFPSbmA3cLWZNSW9GbiaImvrqsea2T+udEKricVtS5XjsYfcRLBUZp1smWN+BkzPRbPqN+taC9iEH0C7PZmWN85Ky5oTGVfMsfRbv1UvvxJoTGWsArY61gaT/nfw3HanPNOF5liqiI854825njacPlQdV07wM6dWZ323zdEDqdLfGs+4g46lan/t7DT7LG1fkbflLtdAJaPe151MryMH/Wy5lZnU4qKdcQ32grw3tvjjbS7LCtv++qmvxc4EB4G+rVjNbJ+Z3dX59zFgD7DTzD5vZidtkr4CXFD22H71NQiCwdLupMBe7TWsDGSPVdKFFFlXb1v2p1cCH13nsSf/fiVwJcDoeLKrEATBBiOCsPQASVPADcBVZna0q/yNFD/5r1/rsd10bNquBZjedkHGwj8Igo1EWAWcApLqFBPj9WZ2Y1f5K4DnAZeZ+fH3cscGQTDcmIlmTKzrQ5KA64A9ZnZNV/nlFGLVM8zMlYByx5Y78bK3uUyVx1MRQCP+5WhsczKf5jKRHnXEhYyQ4cV0NceV0hOpAJpOHFH5p3JpTPgP94nzUvFowU92S3vUGUPmMzPrtKumU+bfMipOeTWTaHb0iOMOesi/v6OHUvGqdjwVN8F3i23tSEW5Vi4mrRf3NBPbtzrvuF3P+f2i6VycTLu1ubRuLveqZZ6RU2WzbwX082vjUuDlwLMk3d15XQG8C5gGdnfK3gMg6XxJN61ybBAEQ87JPdYyr2GlbytWM7sV3/DmJqcMM3sQuGKVY4Mg2AQM86RZhvC8CoJgoJwJdqwxsQZBMHCG2Ua1DJtqYlXLqM8sFQds3I+LaYujpdv1vGVak5l2nbiWlUwyQa9uddZJEJiJDdpyYo62M55XLWe4C2dlYqxucQSazOWyWlq3PZpR0DwnrZYjwDX8fnnl1VR3AqAxldZd2Oo/7tMjad3xfb4qVlmeXA8/Lq4nQgKu6NnOxEJ1F3UtX9mz+fRCVGb8dtVI26hm6o7MLH3OKwtrUEczmEEzAl0HQRD0ls2+FbC5vzaCINhwnNxjPR1WAZIeLek6SR/veeNdxMQaBMHAMVOpVxkkvU/Sfkn3LSu/XNI3JD0g6fXFee3bZvbbfRjSEmJiDYJg4PQ4CMsHgMu7CyRVgXcDvwg8EXippCf2cgwrERNrEAQDxWxNDgJnS7qj63Vl2p59CTi0rPhpwAOdFeoi8BHgBX0fXIfNJV5ZGtty/hw/FuqIo6jXDjlZRIGRA6nnrWXcX71YojnLhMpiagFQdSwIvPiZ4LuONv3hsrA9/faf3+G7PLYmHeU59xVcS1Xi2oRvBVGtZXxVl9Fu+SdrNtJr2zqRuTZuZlx/BTR7dnq+2qwfm9dG0rqL0+n9zWX8rThuql6bAHIyDGvBjynbnk+tGFTNWBtMTSRlnns1gJa5aCvjnr02RKu8VcB6c17tBL7f9X4v8FOSdgB/Clws6Woz+3/W0faqbK6JNQiCoaDs/ukp4J3AzOwg8Dv9PnlMrEEQDJQBxWPdC+zqen8B8GC/T3qS2GMNgmCwWLHPWuZFiT3WDLcDj5N0kaQR4CXAp/o0ooRYsQZBMHDWoPivuscq6cPAMykm4b3Am8zsOkmvAW4GqsD7zOz+U+jymthUE6valiQJrGfcBStOYjkvEWBRnpblErXJEcVyyeI8V9e2I3S1M+JGey2xWyccN1VPpAI06bjVjvqC1NRE6kq5ddy/jqPVtI2K0n7lsnMeW0j9ag/V/fvQMM8HN5P10/lZWmlkYrcedcSrqbQsFxe3upj+wXNjBlArrWtTTmxgoOKIdb6AB1Z36mZ+u1bmep+Uw9YmXq3entlLM+U3kYmm12821cQaBMFwkDFC8Dhb0h1d76/tpGPa0MTEGgTBwFmDVcB6za1OKzGxBkEwUAphanMHYYmJNQiCgRPRrdaJpF2SbpG0R9L9kl7bKX+LpK9LukfSJyRtyxzvBlYIgmD4GYC51WmlnyvWJvA6M7tL0jRwp6TdwG7gajNrSnozcDVF1tblfIAi8eCHSp/RDC0sVVhHvn/Qr+oEBs7tp2vc8RPNxXJ2rA0qGTdEL9B1czp1pVzY6n//LWxLv/UbvkhO0xOTMwGpx6fSa3POtO/ue/7kTFL2iLFjbt2JSibD6DKOZ6Jq719Is6HWKv4YHnJWRM1Mu5VFR+l3ri1AxXG3bTsey5Vs2lOvrj+G5lTacGvCd4/2XKErXkbYDM3pzLUZXTpFWMZNdi0Yot1/l9bTSt9WrGa2z8zu6vz7GLAH2Glmnzezk0/BVyg8IrzjvcAKQRBsAqzka1gZyB6rpAuBi4Hblv3plcBHB9GHIAg2CCFenTqSpoAbgKvM7GhX+RsptguuP8X2rwSuBBirbTmVpoIgGBRnuh2rpArwFOB8YA6438weKtO4pDrFpHq9md3YVf4K4HnAZWZrMBV26FzkawG2jp83zL8eguCM4Yy1Y5X0GApR6dnAN4GHgTHg8ZJmgT8HPmhm7s67JAHXAXvM7Jqu8ss77T7DzBxn0fVjFdGeWroJXzmUiisANFIhRdu2ulVb56Tl7bq/PV07mIo8avhKhucW67VrTmZPgJajqTWn/e+W9oQTN3Xc79eOqfS2PGbLAbfuRRNp+Xn1I27dbdW03ZbjM36k5StwP6iflZTVMr6jXhbQ/Yu+i2drzsmymvncVxfT61s/kZaNHPPdheszqbtvThBqjnuup7ksvGkb9Vw8VscV+vhOP/5sdZn+1f7HXohX0G5v7q2Ala7S/w38FfAYM3uumb3MzF5kZk8Gng9sBV6+wvGXdv7+LEl3d15XUCj908DuTtl7ACSdL+mHfr2dwApfBn5c0l5Jfc9TEwTBADCKb64yryElu2LNBTbo/G0/8LaVGjazW/GDzbpBEczsQeCKMucPgmC4ObUNwI1PmT3WFzrFM8C9nQk2CIJgbZzpEyvw28BPA7d03j+Twv708ZL+i5n9ZZ/6FgTBpqR8ams2q1UAhY/RE05aAkh6JPA/gJ8CvgTExBoEwdoov2LdXFYBXVy4zLxqP/B4MzskqbzP3CCoiNbEUnWzOuGnLbWR1D1y8VzfKmD+nFQxrc36anTtYcedM7Oh1B5LL7+n7uaCrXuulM1x/1zmZE6dmvQDUp87eTQpe8zEw27dHx/bl5TtrB12626tpK6yHofa/j2bdI7PBcU+tJBmIj0y4QeJboykF7LS9C/6yNFU7R+ZST8G1eOZsTrZW23KH0OlmdZtZz6xVk37mwuQ3nLKG1P+eBvLupC53GvDwDa5VUCZifV/SvoM8Ned9y8CviRpEvDtaoIgCFZkuCZWSWcBu8zsnjL1y0ysrwZeCPwsxdX4IHBDx7D/59fb0SAIzmCGQLyS9AUK09IacDfwsKQvmtkfrHbsqhOrmVln83jGzP5W0gQwBfghjIIgCFZjCCZWYKuZHZX0fwLvN7M3SSq1Yl3VjULSq4CPU3haAewEPrnurgZBcGazNgeB0xmPtSbpPODXgM+s6cASdV4NPI1OZCoz+6akR6y5i4NASjbxzREmANpOtktXOMIXqsb+2Y9PytHjSZFt90Wxxe2pKHb8/FQd8OKuAjSmHCFk3HelHJlIBZZcNtWzR9OxnT/iC1K76mm8211VX7jZWknH23YC2062fU/nhqV9eNgRIQG2jcwlZSMjvgvvYrW8SNSccFxHT6T3x4u1C7ifuJybasXJ6FrJyMWV+fS+147697c9kd6H+jn+gJvj/dkLXYODwOm0CvhjivTZt5rZ7ZIeTeHevyplJtYFM1tUx19dUo1hWcgHQbAxGQ6rgH0dF34AzOzbkq5Z6YCTlImo8EVJbwDGJT2Hwjrg0+vrZxAEAcjKvU4z7yxZllBmxfp6Cu+re4F/R+Hr/97SXQuCIOhmg6cHkPTTwM8A50jqtgDYApSy5C1jFdAG/qLzCoIgOEU2fOSqEQrLpxpFJL6THKWw41+VleKx3ssK3yvdew9BEARrYgOvWM3sixRboB8ws++up42VVqzP6/z/1Z3/n4wJ8BtATwNU9xMb9QP4VmZT5Xp8JqP0L6RBsdtHMgG0q05w4tGzM53zi5fTTL0zAWiNp6qxxnyrgImxdAxbRn3VeLqWlk9k3FGnHa/m6Yr/WE04VgEejUz88+lK2i+vDGB8eYRm8hldrZ7eiGYu2+1YutKSp97P+hlp22OOZcSI/+uy7QSkzsT1puJkadVR/3muOM/z2EHfjXhxy9K+VfzHa+1kxrHBGJV0LXAhXXOlmT1rtQNXisf6XQBJl5rZpV1/er2k/wX8l3V3NwiCM5eTdqzlOJ3Rrf4aeA+FprSmr5Qy4tWkpJ/tBK5G0s8Ame/zIAiC1VmD4n867VibZvY/1nNg2Xis75O0leK7ZoYibXUQBMH62MB7rF18WtLvAp8AfrgXZmaHVjuwjFXAncBTJG0BZGaZzcWlSNoFfAg4l2JH5Voze7uktwC/BCwC3wJ+y8ySKFmdpINvpzBveK+Z/bcy5w2CIOgRr+j8/z92lRnw6NUOXMkq4GXA/3cyC6uZHV3298cA553cInBoAq8zs7skTQN3StoN7AauNrOmpDcDV1Nkbe1uuwq8G3gOsBe4XdKnzOwfVxqMmm3qh5e6MlZmfXHDDqXukc2cIOWdq+ZfusrZO5KyVkacqC6kO/hVp7s5wcLbp1LFXwpMjqaCxbYRXyQay/lNOjQcH5OGn7iXllPedpYu85nj5y295g3LCD9OWDrlfn/W0vPl4tpWHK/Y+kPOczOXcSe98JFp1UdkBNblwVCB2qy/1Wde1uCaf220kN7f8QdTV2yA0cNLXcI9N9v1sAGM/1fFzC5a77ErrVh3AF+VdCdwJz9Kf/1Y4BnAAQrngVyn9gH7Ov8+JmkPsNPMPt9V7Sv4dmFPAx4ws28DSPoI8AJgxYk1CIIhwBgKl1ZJv+mVm9mHVjt2JauAt0t6F/AsilTWTwbmgD3Ay83se2vo4IXAxXQCuXTxSuCjziE7ge93vd9LkQomCILNwBCsWIGf7Pr3GHAZcBfFFueKrLjHamYtip/uu9fbM0lTwA3AVd3bCZLeSLFdcL13mNedTPtXAlcCjNX9KFJBEGwshmQr4Pe633cE/FI5/spYBawbSXWKSfV6M7uxq/wVFA4Il3UyESxnL7Cr6/0FwIPeOTo2bdcCbJ04fwhuVxAEQ7JiXc4s8LgyFfs2saqIM3gdsMfMrukqv5xCrHqGWcbFBm4HHifpIuAHwEuAX1/1pM0WlQMlBShHfKpuy6x4nboayyQpdJIXWtUPItZwksgtnOXE9nTihQLgCFXVmi9ujNZS1aWeVcVS5r3MhcARJ/HfpHwhpEEaI7XlfK8+nAmG+nBrS1J2qDnl1j3RTAWhZmsNmfAyW4BykgG6wUUzcYBxYq9amRhzq5B7xvzKaX8rGa/DyqGlz4gafkzbNTMEE6ukT/OjnlaBJwAfK3NsP1eslwIvB+6VdHen7A3AO4BRYHcnxutXzOx3JJ1PYVZ1Rcdi4DUUQWarwPvM7P4+9jUIggGxQUICluGtXf9uAt81s71lDlx1YpX0X4E/O2lr2slW+Doz+88rHdcxw/K+92/K1H8QuKLr/U25ukEQDDlDYBVgZl+U9Eh+JGKVyh4A5QJd/2K3Ab+ZHaZrAgyCIFgrawh0fdpyXkn6NeAfgBdT5L26TdKphQ3soipp1MwWOicbp/gpHwRBsD6GI1bAG4GfNLP9AJLOAf6WIrnqipSZWP8K+DtJ76e4HK8EPrj+vgZBcEYzPHuslZOTaoeDlPuVXypWwJ91cmk/m2LP9E/M7OZ1dbPfWBtbWBY3dMdZbtX2uWm51curxmr6iroaqSqfy9g5tz09X8MRuVu+AQI2lvahXvetAtqO++tcy1euZ1upon6o5avvk00/7qhbt53W9VxS97f8zKv/tJjGtd274N/fg/NpALaFxczj3kjvT86r18uo2jwntVbIxVj1XE9Hj2Qi0jnqfXZCci2/M5VbzjPa8Adss8ssOZo9Csg6HBPr5yTdDHy48/7fUlL3KWsVsIcihNbfSpqQNG1mx9bR0SAIgnz8iw2ApMcCjzSz/yjphcDPUnx1fRnfoSlh1WWtpFdR7Cn8eadoJ/DJdfU4CIJg4/M24BiAmd1oZn9gZr9PsVp9W5kGyuwXvJrCJvVo50TfBB6xru4GQRDAjzK1rvY6PVxoZvcsLzSzOyjStKxKma2ABTNb7BjzI6nGsOyQBEGw8dj44lVG1QBgvEwDZSbWL0p6AzAu6TnA7wKfLtP4wKnWYPu2JUXtLf51aEw7Lo+Tmdie9VQZ8GJlAowe8JPueXihRF33xtxD6OxTNZv+GOab6a2eddw+AY400uyFo14gUqDtdPiY4+YKUFcqfHiusgeavnj13fk01u23j6VlAA8fT8XxWXWnAAAYSElEQVSrxdlM3NM5T7zyDdgbTlKi2Z3p9Vqc8n8M1ufSmzZ2wBeO5IhPrdGMKOa5ymbcar2R5cz1NbbMsnJ+DW7BK7GxJ9bbJb3KzP6iu1DSb1OEUF2VMhPr6ynSs9wL/DuKfYb3rrGjQRAEP2JjT6xXAZ+Q9Bv8aCK9BBgBfqVMA2XMrdqSPgl80sweXm9PgyAIoFgdb2SrADN7CPgZST8PPKlT/Fkz+/uybayUmkXAm4DX0LkWklrAO80sUl8HQbA+Nv4eKwBmdgtwy3qOXckq4CoKa4CfNLMdZradIor/pZJ+fz0nC4IgADa6VcAps9LE+pvAS83sOycLOjmoXtb5WxAEwfrY5BPrSnusdTM7sLzQzB7uZAbYcLRHqsw/aqlVQHXed8FT2wtY7LfbcqwCvDKA+pgn9fvtVh1vUK8sE2MaLTgZUhf8Wzq7kCrix+p+LJ2KswG20PbV4CO1VBF/qJa6eObw3GoPLPjus/tm03YfPpYJdD2TWiZUjvgXsn4svY7VjHGHY9hAayR9FjKxumnXHAuTRf8Zrc2kmV4r0/49aztu0+2t6b0BN7lv1gW3tex5bh/szUd/GLYCToWVJtaVnMDLO4gHQRAs5wyeWJ8i6ahTLlY2oA2CIMhjG9sqoBeslP66R5bAQRAEyzhNK1ZJk8B/p/jV/QUzKxVUZa30II2Zj6Rdkm6RtEfS/ZJe2yl/ced9W1I2gK2k10q6r1P3qn71MwiCwbOGDAKrtyW9T9J+SfctK79c0jckPSDp9Z3iFwIfN7NXAc/v6aC66GcywSZFbqy7JE0Dd0raDdxHMbg/zx0o6UnAq4CnUXyzfE7SZzsBYLLIjGpj6W+M+kE/+6RHdTbjijmRbthbtXzOnsa0f5k9gUOe52juVJ5O1vIrzy+kYzhS9cc720jrjlZ9gWWslrpjjmTqNtvp9/hCK70IJxZ919OZE6l78vzxjJvq4XQM4w/564jx/eknuD6XcVk+nN6g6mL6u7bpiZgZqsd9pUxH0sic1Uw8VNue+tp6btsA5ghobvZZnFjCvUpV1dsV6weAdwEfOlkgqQq8G3gOsJfCTfVTwAUUXqQAPQoum9K3FauZ7TOzuzr/PkYR03Wnme0xs2+scvgTKLK3zppZE/giJV3JgiDY4JQ1tSqZ88rMvgQcWlb8NOABM/u2mS0CHwFeQDHJXtCp07f5r58r1h8i6ULgYuC2kofcB/yppB3AHEXywjv60rkgCAaKWJO51XpzXu0Evt/1fi+Fg9M7gHdJ+jf0MZhU3ydWSVPADcBVZuZZGSSY2R5JbwZ2A8eBr1FsLXjtXwlcCTA6urUnfQ6CoL8MwI7VT1ZjdgL4rX6fvG9LYYCOI8ENwPVmduNajjWz68zsX5nZv6ZY5rv7q2Z2rZldYmaXjNSdmG5BEGw8+u95tRfY1fX+AuDBU2pxDfTTKkDAdcAeM7tmHcc/ovP/R1GIXR9e+YggCIaG/k+stwOPk3SRpBHgJcCnTqnFNdDPrYBLgZcD90q6u1P2BmAUeCdwDvBZSXeb2XMlnQ+818yu6NS9obPH2gBebWaHVzuhWkbt4NKskjqasQpwggBrzFdRaydS5bvtZNsEaE6l7bZGfSnVc4+sOGVu8GvAquWfvOZiqlKfkO8eOVdxsoNmfrtVq+Utvb0g3M1mOrhWIxNwfC59XDXvX5zaifSaTzzkj2H6+6kjYU4lrzTS8VZm0+ej7mRCLSo7QbWPzTkVgbZzbb0yoDrrjcGvW5lP+6s535nSJpY+I97410yPo1tJ+jDwTAqhay/wJjO7TtJrgJsp7GfeZ2b39+6sK9O3idXMbiVvnPEJp/6DFCLVyfc/16euBUFwuunhxGpmL82U30TJdNW9ZiBWAUEQBN2csS6tQRAE/eJMjm4VBEHQe4Y81moZNtfE2m6jhWWb8KMZt75xR7hxYlqCH6syly2z7bi65jK6jhxPfw8tbnHaXYsbYduv3F5I211slTcKcRKGAiBH6Mq20XDOt+iIOU6cWYD6gnNtnTKA2mxa5rmeAtSPOsJN06/bHnMENEeoygpSHsrc4JpzrrmM++tiKkhVMu3aTOoqawt+u5Ud25cW5ES5tRITaxAEQe9Yo+fVUBITaxAEA8fN4LGJiIk1CILBEnusQRAEvSe2AoIgCHpNTKxDhJS4DLYzVgHtqbS8Ne5fjsUtjjqbMXCuzaaqqZd5FaA1kqrfNcdFs5IKvkUfPPU8E+iaphPceA1G2jnDBNfdNrMcqTh9k9Ov3PXy3FTrx/261YW0DxU3PhpowXHxnM2kacXJfOqp7xnXUxw3U9uSCR5USy05dCzjop0z23D74Cj7OcuEPhEr1iAIgl4TE2sQBEEPOZOztAZBEPSDsGMNgiDoB2vZEx5CNtXE2h6pMv9j25aUVTJujF6W1daY70q5OOWJTJksnocc90Yn/iVAZTyN3Tq5z+nXiH+bPOEnF7rci/2ai/ParjnxWDOusp6wlvuZZ04TXt3c8e65Mh6WaxLmnNirOu74xAIVJ56qObF9vbirADhZVi1T1yaddif9zLotr25GkKrNbEmrZlx4F89ZKqzZod5MGbFiDYIg6CXhIBAEQdB7QrwKgiDoMZt9Yu1nMsFdkm6RtEfS/ZJe2yl/ced9W1I2X7ik3+/Uu0/ShyX5m0tBEAwXRiFelXkNKf1Mf90EXmdmTwCeDrxa0hOB+yiyrn4pd6CkncB/AC4xsydRJAN7SR/7GgTBAJGVew0r/UwmuA/Y1/n3MUl7gJ1mthtAq7vQ1YBxSQ0KP8JVc4K36+L4eUvV0Zx6X59Nf4tU53yJefyAoxr34KZ7arQXjHnscO53U/nvxeZ4WtbK/Aaoem6mGbfa+rF0DKNHMi6tzbS8NufUzTwa89u8gON+3dGjabtjB/1BeMGj2yd8qwA1U7/YyvRUUmbVzL1xAqx77tUADS/j77gfYH3+LKc884xO7E/r1k74/r5z5yztWy478ZoZ4kmzDP1csf4QSRcCFwO3lalvZj8A3gp8j2JynjGzz/erf0EQDI6TDgKbecXa94lV0hRwA3CVmR0tecxZwAuAi4DzgUlJL8vUvVLSHZLuaM5nAlQEQbBxMEPtcq9hpa8Tq6Q6xaR6vZnduIZDnw18x8weNrMGcCPwM15FM7vWzC4xs0tqY5koQUEQbCys5GtI6dseq4pN1OuAPWZ2zRoP/x7wdEkTwBxwGXBHj7sYBMFpYph/5pehn3aslwIvB+6VdHen7A3AKPBO4Bzgs5LuNrPnSjofeK+ZXWFmt0n6OHAXhXXBV4FrVz2jQXVx6R1rTGSUEKWL9ZEjvrhROzqflDW3+MrP4pZUcKhlNvw9t1ovS2tj0h+DJ35lqTgZTjPxSevOhk39RHlBavSILwKOHkyvY/VYWmZVX6AZ2ZEqcOaMC6B+OG3Xi7sKgJPhVDnxycNxU8VzcwU3xmqOinN/czlS284nORt/1rmVlYxLa/3Y0jOu6ZnLYcAQ/8wvQz+tAm4lHx/5E079B4Erut6/CXhTf3oXBMFpZXPPq+F5FQTB4ImtgCAIgh4zzIp/GWJiDYJgsAy54l+GTTWxVprG+MGlO/ZzOzKxTL0YnJkNfDWcGKsLvjJQq3kJ88p7TnlJ8FzvJKBVT8+Vi7Fac0x8c+2OznhBUv12m6NOHxxRLovngZc5Xo20X5XMykdzmYyEDm4yv62pNxVAeyx9nqzuCFK5fjljqDriKED1mCM4Lviiqdqp91YuVm3tuCPWLWQEx0NLvdLkiJVrpXAQ2Nwz66aaWIMgGBI2eXSrmFiDIBg4sWINgiDoJbHHGgRB0GuGOw5AGWJiDYJg8MRWwPCgllE/ukwNzn0zui6emYyuThZNtfy61bnUWqBdy7i0On0YOZyq2bUTvhvk3COceJ0jvqI+fjhVfcce9pVzz4phcZsfM7TlKOLNiUy227PSWKR179pkblnOfdWtW08fbc2ncVcB2hNpvxZ3+Op7Yzpt13MnHTmWcRGdSa+5Wr4i77r2Zq7N2P7UsqBywh+v64LrXC+AdjJF9MaldbOnZtlUE2sQBENCD1esko4C0z9s2XJGh4PjtHcgCIIzkB6FDZS0g2JS/RvgKUWRbu9Ln9dArFiDIBg4apfeC/gJaUlkgZaZdc9bXwMwsysAOnV/oiedPAVixRoEwWAxCgeBMi+408zU9Vq+GNyx7H2TfFS9gbHpV6wjTlxOyLghZjDHjbGdOd4cMaY54df13Grr+4+l9RZ8kanSODspa2zxb2l9JnVj9M4FoIbjrls5y63bmE7Hltvhak46SezmUyGlNpMTXcorHp4BujsuoOLGafXFK29snmDYGvM/27XjaXl7IiMMjpf/eFYPOaKYkySxOGF6HXPJDxMX3tWTgK6KsE3vIBAr1iAIBo9ZudfqHFz2vsYGcD+IiTUIgsHTu4n16QCSPi3pyRTbAF/tZ9fLEBNrEASDZW17rCs3ZbYXOA48j0LIMjM77eLVpt9jDYJg47EGq4BVMbPp1WsNlr6tWCXtknSLpD2S7pf02k75izvv25IuyRz745Lu7nodlXRVv/oaBMEgKbkNMMQCVz9XrE3gdWZ2l6Rp4E5Ju4H7gBcCf5470My+ATwVQFIV+AFOAsIy5IJMt0ccpT+j3lccN9XqrK/Ut1yF1283cb8FePhQ2q8FX90dcVxta1udoM3gKuqazVhMzHvZVNMMqQAjM07g51ygascKono8vQaVo7N+vxzl2hx3VIC246JZnctYZ8wcT8pGnWsLUGlOJGWNKcd9NudK7ajqOZdnr41KJiA1XvbU3MTkuWNnXLT7ot4bQz1plqGfWVr3Afs6/z4maQ+w08x2A6i82cZlwLfM7Lt96WgQBIMnYgWcOpIuBC4GblvH4S8BPrxC21cCVwKMjWxdR/NBEAyasGM9RSRNATcAV5nZ0TUeOwI8H/jrXB0zu9bMLjGzS+r1zM/gIAg2FrHHun4k1Skm1evN7MZ1NPGLwF1m9lBvexYEwWnDLLunu1no28SqYhP1OmCPmV2zzmZeygrbAAlmaUbV3Lees8VrmX3fynzq8lg5nAoeRbtbkqLGVt9lUU7XNJaKMRr1j3fJuX2u5beJvPiz/nWsH/XcQXPtlqvWnvLdST3Xz+ZEGpMW/J+aldnMfZhzYplmhL3aCf98SZu5jL+OINUezYhqTt3cfcARwGzcF/bwXHsz8ViT8/VqETnEq9Ey9HMr4FLg5cCzusymrpD0K5L2Aj8NfFbSzQCSzpd008mDJU0AzwHWs9INgmAjE1sB68PMbiW/RklMp8zsQeCKrvezpJFrgiAYdox8Zo9NQnheBUEwYAws9liDIAh6hxHiVRAEQc8Z4v3TMmyuiVUqHcC65QSvzu0Iay51u2wfmXHrVpxMovWMct3y1OBz023lyglfofYCGecMr5vTjktqZtHgBchYfIRvI9x23Fdrs5mA0g0nwLJzv1rj/vVqbEnLc5lba3Op62d71G+3Mp2OzcYzwacnyn1kqnPlrSXaI+WDrrfrmYy/3vO81bcK0KLj3px5bhKLhV7J3TGxBkEQ9JLhVvzLEBNrEASDxXDTw2wmYmINgmDwxIo1CIKgl4RL61BhFSUujpWWLwx4GUOzMTSd2Jy5sIfmuELWZnzxqbE9je2ZZMUEqjN+DE6bT8Ur25G61AIsnJWKMfWRTMxRJ+Znc6y822VlISNeHZ1zGnCu43ZfKGuOO/ch58F7zBPKMplIt6b3oTmZExzTNqqeGJRxPU1croGKExu4qJwWeQIgQNMRwLxnHEDNtG+1Wf8ZWx4rNuf2vSYMLOxYgyAIekx4XgVBEPSY2GMNgiDoIWZhFRAEQdBzYsUaBEHQSwxrZZIibhI21cRqVbG4demQKo4CmqN2InOzHbdJbfXVdy9gcDsXQNtRiJXLwumgydRNtZFxY1zY6ijq5qvGY0dSa4OJI46iD9B0xpDL/jrntOF8wHIOniNbUsuGiqPIg59Z1wssDn6W1FbOCsIJYF2ZX8MksYaVmpeVloxVgNcvNTPZXx2LhZxFjPXACCBtlBCvgiAIek6YWwVBEPQOAyxWrEEQBD3EItB1EARBz9ns4pVsE5k9SHoY+O4KVc4GDgyoO4MkxjV8DOvYfszMzjmVBiR9jmL8ZThgZpefyvlOB5tqYl0NSXeY2SWnux+9JsY1fGzmsQX9TX8dBEFwRhITaxAEQY850ybWa093B/pEjGv42MxjO+M5o/ZYgyAIBsGZtmINgiDoO5tiYpX0+5Lul3SfpA9LGuv62zslHV/h2KslPSDpG5KeO5gel8cbm6TrJH1N0j2SPi5pyjmuLumDku6VtEfS1aej/znWO67OsU+W9OXO8fd23+/TzamMq3P8oyQdl/SHg+x30GPMbKhfwE7gO8B45/3HgP+j8+9LgL8EjmeOfSLwNWAUuAj4FlA93WNabWzAlq461wCvd479deAjnX9PAP8EXHi6x9SDcdWAe4CndN7v2Cj37FTG1fX3G4C/Bv7wdI8nXut/bYoVK8WHbVxSjWISeVBSFXgL8J9WOO4FFJPPgpl9B3gAeFrfe7s2krGZ2VEAFYm3xincr5djwGTnuHFgETg6mC6XYr3j+gXgHjP7GoCZHTSzjeTGs95xIemXgW8D9w+or0GfGPqJ1cx+ALwV+B6wD5gxs88DrwE+ZWb7Vjh8J/D9rvd7O2UbghXGhqT3A/8M/Avgnc7hHwdOdI77HvBWMzs0iH6vximO6/GASbpZ0l2SVvriHCinMi5Jk8AfAX88sA4HfWPoJ1ZJZ1GsPC8CzqdYpf0m8GL8D+aSw52yDWMmkRnbywDM7Lc6ZXuAf+sc/jSg1alzEfA6SY8eRL9X4xTHVQN+FviNzv9/RdJlg+j3apziuP4Y+H/NLKsHBMPD0E+swLOB75jZw2bWAG6keEgfCzwg6Z+ACUkPOMfuBXZ1vb8AeLDP/V0L3th+5uQfOz+BPwr8qnPsrwOfM7OGme0H/hfFnvNG4FTGtRf4opkdMLNZ4CbgXw2gz2U4lXH9FPBnnef1KuANkl7T/y4H/WAzTKzfA54uaaKzh3UZcI2ZnWtmF5rZhcCsmT3WOfZTwEskjUq6CHgc8A8D6/nqeGPbI+mx8MM9u18Cvp459lkqmASenql3OjiVcd0MPLlzbA14BvCPA+r3aqx7XGb2c13P69uA/2pm7xpc14NeMvRhA83sNkkfB+4CmsBXWcGrRdLzgUvM7P8ys/slfYzig9kEXr2RhJAVxvb3krZQbGV8Dfj3sHRswLuB9wP3deq938zuGfwoUk5lXGZ2WNI1wO0U2zY3mdlnT8c4lnOK9yvYRITnVRAEQY/ZDFsBQRAEG4qYWIMgCHpMTKxBEAQ9JibWIAiCHhMTaxAEQY+JiTXIIqkl6e5OpKZPS9q2Qt2LJZnWECFM0tsk/Wun/JmSPnMK/f7bjhdUEJwWYmINVmLOzJ5qZk8CDgGvXqHuS4FbO/9fFUnbgaeb2ZdOvZsJfwn8bh/aDYJSxMQalOXLZALUdDyKXkQRIu8XSsZHfRHwua42Lpf0dUm3Ai/sKp+U9D5Jt0v6qqQXdMonJH2sE+P0o5Juk3TSZfdTlJzgg6AfxMQarEonBONlFBOWx6UUPvLfAr4AXFGi2UuBOzvtjwF/QeHu+XPAuV313gj8vZn9JPDzwFs6Lrq/Cxw2sycDfwL8xMkDzOwwMCppR9kxBkEviYk1WIlxSXcDB4HtwO5MvZcCH+n8+yOUWy2eBzzc+fe/oJiYv2mFK+BfddX7BeD1nX58ARgDHkUR2eojAGZ2H0Xw6272U0STCoKBM/SxAoK+MmdmT5W0FfgMxR7rO7ordFazvwo8X9IbKfzhd0iaNrNjK7VNMUmeJOdbLeBXzewby87rhXzsZqxzjiAYOLFiDVbFzGaA/wD8oaT6sj8/G/iame3qRGf6MYr0Ir+8SrN7KEI7QhHt6SJJj+m8717x3gz83smJVNLFnfJbgV/rlD0R+JcnD+jUPZciHU0QDJyYWINSmNlXKSIzvQSg89MciknwE8uq30ARDxZJN0nyfpJ/Fnhmp+154Ergsx3x6rtd9f4EqAP3SLqv8x7gvwPnSLqHIvL+PcBM528/AXzFzJrrGmwQnCIR3So4bXQm0eeZ2ZF1HFsF6mY231np/h3weDNblPR2irQ8f9fjLgdBKWKPNTidvI5CiFrzxEqRqO+WztaEgH9vZoudv90Xk2pwOokVaxAEQY+JPdYgCIIeExNrEARBj4mJNQiCoMfExBoEQdBjYmINgiDoMTGxBkEQ9Jj/H03zi6SCZqFuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10895c470>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Slightly smooth the map for display to suppress statistical fluctuations\n",
    "skymap.skymap().smooth('GAUSSIAN',0.02)\n",
    "\n",
    "from matplotlib.colors import SymLogNorm\n",
    "# The SymLogNorm scale is a Log scale for both positive and negative values\n",
    "# and is linear within a certain range around 0\n",
    "\n",
    "ax = plt.subplot()\n",
    "plt.imshow(skymap.skymap().array(),origin='lower',\n",
    "           extent=[83.63+0.02*20,83.63-0.02*20,22.01-0.02*20,22.01+0.02*20],\n",
    "           # boundaries of the coord grid\n",
    "           norm=SymLogNorm(1)) # the scale will be linear within +-1 count \n",
    "ax.set_xlabel('R.A. (deg)')\n",
    "ax.set_ylabel('Dec (deg)')\n",
    "cbar = plt.colorbar()\n",
    "cbar.set_label('Counts')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model fitting and residual inspection\n",
    "\n",
    "First, fit the usual Crab model to the data using a 3D unbinned analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "like = ctools.ctlike()\n",
    "like['inobs']   = evfile\n",
    "like['caldb']   = caldb\n",
    "like['irf']     = irf\n",
    "like['inmodel'] = '$CTOOLS/share/models/crab.xml'\n",
    "like.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As usual it is wise to look at the output from the optimizer. Store the best-fit minus log-likelihood value for later usage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== GOptimizerLM ===\n",
      " Optimized function value ..: 314283.848\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 ......: 10\n",
      " Lambda ....................: 1e-09\n"
     ]
    }
   ],
   "source": [
    "print(like.opt())\n",
    "like1 = like.opt().value()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fit has properly converged after 10 iterations. Let's take a look at the best-fit model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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",
      " 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 ................: 3.53420581059778e-16 +/- 6.13986803854443e-18 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)\n",
      "  Index ....................: -2.27926367934528 +/- 0.0145867159374822 [-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",
      " Number of scale par's .....: 0\n",
      "=== GCTAModelIrfBackground ===\n",
      " Name ......................: CTABackgroundModel\n",
      " Instruments ...............: CTA\n",
      " Observation identifiers ...: all\n",
      " Model type ................: \"PowerLaw\" * \"Constant\"\n",
      " Number of parameters ......: 4\n",
      " Number of spectral par's ..: 3\n",
      "  Prefactor ................: 0.996370165375849 +/- 0.0082156655736933 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient)\n",
      "  Index ....................: -0.00225668505770844 +/- 0.00508689620095184 [-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": [
    "The background model reproduces well the data, since its best-fit normalization is 1 and the index/tilt is 0. The Crab nebula has a best-fit spectral index of 2.28.\n",
    "\n",
    "To assess whether the model reproduces well the data spatially generate a residual map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "resmap = cscripts.csresmap(like.obs())\n",
    "resmap['algorithm'] = 'SIGNIFICANCE'\n",
    "resmap['emin']      = emin\n",
    "resmap['emax']      = emax\n",
    "resmap['nxpix']     = 40\n",
    "resmap['nypix']     = 40\n",
    "resmap['binsz']     = 0.02\n",
    "resmap['proj']      = 'TAN'\n",
    "resmap['coordsys']  = 'CEL'\n",
    "resmap['xref']      = 83.63\n",
    "resmap['yref']      = 22.01\n",
    "resmap.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To easily inspecting the residual maps along the tutorial, define a function to plot them, and apply it to the latest run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_residuals(resid):\n",
    "    # Slightly smooth the map for display to suppress statistical fluctuations\n",
    "    resid.smooth('GAUSSIAN',0.04)\n",
    "    # Plotting\n",
    "    fig = plt.figure()\n",
    "    ax  = plt.subplot()\n",
    "    plt.imshow(resid.array(),origin='lower',\n",
    "               cmap='bwr',vmin=-3,vmax=3,\n",
    "               extent=[83.63+0.02*20,83.63-0.02*20,22.01-0.02*20,22.01+0.02*20])\n",
    "               # Boundaries of the coord grid\n",
    "    ax.set_xlabel('R.A. (deg)')\n",
    "    ax.set_ylabel('Dec (deg)')\n",
    "    cbar = plt.colorbar()\n",
    "    cbar.set_label('Significance ($\\sigma$)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVYAAAEKCAYAAABJ430PAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJztnXu0XVV97z/fhEACBBESL5KER31Vh+Wh0apYrYAtpqhF8RarSKW31Epb6KCjitxRaq0dWloqaltNi4pK8cFDrVAhWpRShUpieHmoz2KR3JKIvAoJhHzvH2sd3Zwz5zlrnrXWPvuc/D5j7JHsueeca8699vntuX9P2SYIgiDojgWzvYAgCIL5RgjWIAiCjgnBGgRB0DEhWIMgCDomBGsQBEHHhGANgiDomBCsQRDMSSQtlvTvkm6UdKukt8/2msZR+LEGQTAXkSRgD9sPSFoEXAucZvu6WV4au8z2AoIgCGaCq1PhA/XTRfVjJE6KvQlWSauAjwL7ATuAtbbPk3QO8HLgYeC7wBtt35OZYyFwA/BD28dOd81ly5b5oAMPbLrASU25w/tsH+pz109sIUtJ3x07mrVBem3DfB9z+1qQUHKV9G1LyV5z60q1l9zHLC1uxH/efjtbtmxptYpjJG9p2Hc93ApsHWhaa3vt+JNaRqwHngz8je3r26ytK/o8sW4HzrC9QdJSYL2kdcA64Ezb2yW9GzgTeEtmjtOAMWCvJhc86MADueGrX31sY+6vZpfJW39ke/rzsn375LZHH22yonJSn/mcUEttLbGtbHtu3q1bm7UBPPxwszZIv4+5NaRI7XfXXdN9Fy9u3jfVnnsfU/e95J6lWLgw3b5oUfN1Jb9IdmQ+pCWLm9B39Qte0Hxshi1Up6UmCLbaXp173fajwGGS9gYuk/RM27e0XmRLejNe2d5ke0P9//upBOQK21fZHv8Tuw5YmRovaSXwK8A/9LXGIAhmiQULmj0aUv/q/TJwTF9LLmEoXgGSDgIOByYe008G/jkz7D3AH1GpEaaa+xRJN0i6YfPmzS1XGgRB70jV0bvJY8pptLw+qSJpCXA0cNsQdjAtvQtWSXsClwCn275voP0sKnXBhYkxxwJ32V4/3fy219pebXv18uXLO1x5EAS90c2J9YnA1ZJuAr4OrLP9+d7X3oBevQJqF4hLgAttXzrQfhJwLHCU0/5eRwCvkLQGWAzsJenjtl9fvIiCnxMl+skSciqtVHuJPjelm8ttN9VeorvNzZsyppS8X6n9PvJIu/GQ1vOWvDcllKgs21Jyz7IbK7FO9oHUicXQ9k1Uv4RHjt5OrLWP2fnAmO1zB9qPoTJWvcL2g6mxts+0vdL2QcAJwL/MSKgGQTCadKxjHTX6XPkRwInAkZI21o81wPuBpcC6uu0DAJL2l3RFj+sJgmBUmOeCtTdVgO1rgZT/UlJ42r4TWJNo/zKVtS8IgvlAR6qAUSYir4IgGC7jXgHzmPm9uyAIRpM4sc4xGt4wJ7QUyoUZt7SiLsx+O3cRn/hYuvi8trX0596alLU/NW/uWjkPgKZ9cxFhqfWWRDiVvOepvrkI0xYBUvW10p8vlSy4L2+BEKxBEAQdEjrWIAiCHgjBGgRB0CFhvAqCIOiBOLHOcTI3MK2TTyv7F5bEgza/GDt2ZPLFtaAvW0MX4b5NbYA5I1PK+JUL920aLlxK0/2WRJPmSK23xLCXN8BN/tx1b0adgtCxBkEQ9EAI1iAIgg6JE2sQBEEPhGANgiDokPAKCIIg6IE4sc4hErqbVOgqFEap7jJ5jnyY6mRya9iesH6XWK5TFvESa3QuoXRfBQ2bhoOWJK/OeRBs2za5rcSDINe3j6TnRcmrM6TWlXsfU32Tni8puigTGzrWIAiCHgjBGgRB0DEhWIMgCDokjFdBEAQdEzrWmSNpFfBRYD9gB7DW9nmSzgFeDjwMfBd4o+17JoxdDFwD7Fav8WLbZze68IQbVlIhNWfMKcl12dYYUxLyuGhRu/E5A01qvbl5UwePnH2jae7VkvDZkr/PEsPgsCv2th2fu5cpku9DwkALBUatUua5YO1zd9uBM2w/HXgecKqkZwDrgGfaPgT4FnBmYuw24EjbhwKHAcdIel6Paw2CYJhEMcGZYXsTsKn+//2SxoAVtq8a6HYdcHxirIEH6qeL6kcmz3oQBHOKnUAVMJTdSToIOBy4fsJLJwP/nBmzUNJG4C5gne2JY8f7nSLpBkk3bN68ubtFB0HQH/P8xNr7yiXtCVwCnG77voH2s6jUBRemxtl+1PZhwErguZKemem31vZq26uXL1/e/QaCIOiWca+AJo85Sq8rl7SISqheaPvSgfaTgGOBo+qf/Vls3yPpy8AxwC09LjcIgmExh0+jTejTK0DA+cCY7XMH2o8B3gK82PaDmbHLgUdqoboEOBp493TXtOGR7c1C7lLW1ZyIT1nvu7BGN7W+5764U+vNWYdTfUu8FUo8CEooqXpaEj5b4kHQdg0ti/h20je1rtaeAsCOCd4vUx+DGtKRjjXnedR64g7o88R6BHAicHOtKwV4G/BeKjeqdZXs5Trbb5K0P/APttcATwQukLSQSl3xKduf73GtQRAMk25OrOOeRxskLQXWS1pn+5tdTN6GPr0CriVd8eGKTP87gTX1/2+iMnYFQTDf6OjEmvM8AuavYA2CIMjSsY51Cs+jWSEEaxAEw6UsV8AySTcMPF9re+1jp0t7Hs0m80qw2s2NKSnFfkn4a1bZX2DIaGq0KBm/MFP4tWRdJSGtKWPG1q3N15AaX+Jlk7sPbQ1SbdOOlnw+cpQY60rWm/rsNzV0dWK8gpIT6xbbq3Mv5jyPZpt5JViDIJgDdOcVkPQ8GgXmtzNZEASjSTeRV+OeR0dK2lg/1vS/+OmJE2sQBMOlO6+AnOfRrBOCNQiC4TOHw1WbML93FwTB6LETZLead4K1D+tmW6+AHKm1pqz6JdbsHCV9SyzMqf3mKqem3rNU37ZVXnOU9M19Pprut8STo2S/uXuT+tyUhLR2Zu1vSgjWIAiCDokTaxAEQQ+EYA2CIOiYEKxBEAQdEuWv5x5NlfB9hJ6WkjI49FW1NFWBNjdvWwNcSahsifEq9X6lKtVC+5DWHKk9pNpyhqO2xsnc5zt1vaEbpJoSOtYgCIIeCMEaBEHQMSFYgyAIOiRUAUEQBD0wzwVrb7uTtErS1ZLGJN0q6bS6/RxJt0m6SdJlkvZuOjYIgnlAlL9uRbLQF7AOONP2dknvBs6kqto67djpioTZky3Sw7b0tw2x3HXXZm1Q5hWQoiRsswuaJiHPhcSW/J21tb6XfG5KQkdL7k9flV7bhgZ3QpxYZ4btTbY31P+/HxgDVti+yvb4n9h1wMqmY/taaxAEQ2Rcx9o+H+vIMpSz9hSFvk4GPjnDseOvnwKcArBq1QGt1hkEwZCYw0KzCb3vLlfoS9JZVD/5LywdO4jttbZX2169777Lu118EATdEyfWduQKfUk6CTgWOMpOa/pGtUhYEAQdMIcNU03obXe5Ql+SjqEyVr3Y9oMlY5vQxniV69uXsr+p8amvfKy53J597bfESNSUEsNRSZhpCX3l0O0rD3BJ9deJdBImuxP4sfa5u1yhr/cDS4F1ddsHACTtL+mKacYGQTAfCFXAzJii0NcViTZs3wmsmWZsEARznZ3gxDq/FR1BEIwmIViDIAg6JgTr3MGerNzvK/Jq2AalFG0jc0oir3JGk9Qcfb23KbqIEkvtoaRoX0kEXNu95fabyrebo8RYN3G9nRmvwisgCIKgQ0LHGgRB0AMhWIMgCDomBGsQBEGH7ASqgPm9uyAIRpOOAgQkfUjSXZJuGcKqGzOvTqwl+VhLQiHbkvt8pNbWNjdojqYVUnPtub6peduGV5bQxXtTEuLZtG9fuV/7YpjX6tgr4CNU0Zwf7WrCLphXgjUIgjlAh6oA29fUqUU7RdIewFbbMzqChWANgmD4NBesyyTdMPB8re21XS9H0gLgBOB1wHOAbcBukjZTheGvtf3tpvOFYA2CYPg0F6xbbK/ucyk1VwNfpCoVdYvtHQCS9gFeArxL0mW2P95kshCsQRAMl9H0Cjja9iOSzrZ903ij7bup8kJfUueIbsROK1hLQvOGaehqWnAvR84I0Zfxqg9DVYmxr4t8uyXGp6b7zY1Phcrm+g4z6rOLENzGjGBIq+3xoOCzJe0O7ANsAD5h+8cT+kzLyH1tBEGwE9Cdu9VFwNeAp0m6Q9JvtlyZga3AlcAq4KuSDi2dZLS+NoIg2DnozivgtZ1M9FNus312/f+LJX0E+ABwZMkkcWINgmC4jHYxwS2Snj3+xPa3gOIqpXFiDYJg+Iye8Wqc3wc+IWk9cDNwCPD90kl6252kVZKuljQm6VZJp9Xt50i6TdJNki6TtHdm/EiGqgVB0JIRPrHavhE4DLiobroaKFY39Hli3Q6cYXuDpKXAeknrgHXAmba3S3o3ld/YWxLjP0JhqFrK2JizZqcoqdKas96nrL65xMJNwyPbhq5C+n3YurV53y6s/02rv3bhFVCyh7ZhuW1DWkuSapd4EHRRdXjiHCWJ2KdkxLwCJMmufIVsbwMurx/JPtPR21eC7U22N9T/vx8YA1bYvsr2+Ef5OmBlZvw1wN19rS8IglliNE+sV0v6PUkHPHap2lXSkZIuAE5qOtlQvjbqWN7DgesnvHQy8MlhrCEIghFi9HSsx1DJo4skHQzcAyyhOnxeBfy17Y1NJ+tdsErakypy4XTb9w20n0WlLriw5fynAKcArFx5wDS9gyCYdUYw8sr2VuBvgb+tI6yWAQ/Zvmcm800rWOvkBIcC+wMPAbfa/u8mk9cLvAS40PalA+0nAccCRzXVWeSoEzKsBTjssNVdlDoLgqBvRkywDlJHWG1qM0dWsEp6EpVR6Wjg28BmYDHwVEkPAh8ELhhPVpAYL+B8YMz2uQPtx9Tzvtj2g20WP/mak3XibUNEoX2uytxXRx/6+5L8s13k4CwxhPRhrOsihLfp+Kmu17TfrrtObluUiUAvqabalqEar0YwpLVrpvra+DPg48CTbP+y7dfbPt72IcArgMcBJ04x/oj69SMlbawfa6gs/UuBdXXbBwAk7S/pivHBPYSqBUEwKoye8apTsl8bU4WK2b4LeM9UE9u+Fkh9v12RaMP2ncCaJtcPgmAOM4I61nHqX9qvA37G9p/WXgL72f73knma6FhflWi+F7i5FrBBEARljKhgpTJg7aDKDfCnwP1UdqLnlEzSRNHxm8DzqSIQAH6Ryv/0qZL+1PbHSi4YBEEwwoL1520/S9I3AGz/WFJCMz41TQTrDuDp454Akv4X8HfAzwPXACFYgyBozgirAoBHJC2kSh+IpOVUMrCIJoL1oAnuVXcBT7V9t6TGiV+HQcn9amvd7KuqZYmVPbWHlNU5R85boaSSaNPxU7W3GZ9LQl7iFZDa2yOZT3Yf9z2335KQ1r7WMJGdwCvgvcBlwBMkvRM4Hvi/pZM02d2/Svo88On6+fHANXUVwxk5zwZBsJMzoidW2xfWma2OojK+/6rtsdJ5mgjWU4FXAS+sL3QBcEnt2P+S0gsGQRCMqmAFsH0bcFubOabdXS1AbwAut306lbvUnm0uGgTBTsxoJmGpl6YLBlOZSnq8pA+VzjPtyiX9FnAxVaQVwArgM6UXCoIg+AkjKliBQwbzA9SFBA8vnaSpKuC51JmpbH9b0hNKLzQM7ObGhZThpouKnylyBqWmIZ6LF6fHl+j/U8acbdvSfVOGm1yIZ6q9bVXbnFGtbahrX/lYSwyOKYNUicExR2q9uftQEio7cd52mT1qRtsrYIGkx49XZpW0DzNIVtVkwDbbD6s2B0rahdoVIQiCYEaMrlfAXwFfkzRurH8N8OelkzTZ3VckvQ1YIumlwJuBfyq9UBAEATDSJ1bbH5V0Az+tyvoq298snafJ7t5KldnqZuC3qYxXxX5dQRAEP2FEdaySdqOqebUXsA9wvKQ/Lp1n2hNrnRbw7+tHEARBO0b4xAp8lioXynogY4WYnqnysd7MFLrUOn1gEARBOaMrWFfaPqbtJFOdWI+t/z21/nc8J8DrgE4TVHeFPTlsMRfGmGpvm9wYyip2pqzBKQ+ARWTiK+9L3IbMhlO5lBflDAiJRXjxkmTXVKXXkuqvudDRUSVlFU99bnL7KvFGKfGuKPFsaOO10ZlXwOgar74q6eds39xmkqnysd4OIOkI20cMvPRWSf9GlVIrCIKgnNE9sb4Q+A1J36dSBYgqTqroF3qTr409JL2wTlyNpBcAe5SuNgiCAOhUx1qXejoPWAj8g+13tZzyZe1X1Twf64ckPY5K53ovVZnYIAiCmdGBYK3T+/0N8FLgDuDrkj43E/eocWzfLunxwFOoavyNc3vJPE28AtYDh0raC5Dte4tWGgRBMEh3J9bnAt+x/b1qWn0CeCUwY8Eq6f8ApwErgY3A86hq7x051biJTOUV8HrgH8ersNq+b8LrTwKeOK4iSIxfBXwU2I8qUexa2+dJOgd4OfAw8F3gjana3TM54u/YMdlwUhKK2UXO0ZROPheymKrOmTRU3ZPJznjffZPbHngg3TdlOcptYrfdJjXpcY9Ldl2y996T2hbvvTTZN7W0EuNX6v6U5JTNhXKW/I13ko+0xfVzhqfUe1Py2c8xcd7O8tE23/Sy2mF/nLV1yXuo8pb818Brd1Al4G/DaVRlWK6z/RJJPwu8vXSSqU6s+wLfqHMTruen5a+fDLwY2EIVPJBjO3CG7Q2SlgLrJa0D1gFn2t4u6d3AmVTlsH9CH0f8IAhGh0d3NP6G2mJ7dea11CRt/Ra22t4qCUm72b5N0tNKJ5nKK+A8Se+nOgIfARwCPASMASfa/sFUE9veBGyq/3+/pDFghe2rBrpdR5U4eyKdH/GDIBgNSpIlTcMdwKqB5yuBO9vOWacN/AywTtKPZzLnlDpW249SnTDXzWiJNZIOokq9df2El04GPpkY0viIL+kU4BSAFSsOaLPMIAiGREeC9evAUyQdDPwQOAH49TYT2j6u/u+fSLoaeBzwhdJ5evfSlbQnVfnY0wf1tJLOolIXXJgalmhLHvFrfctagEMOWR1Zt4JgxOnqxFqrE38XuJLKFvMh27e2n/kn839lpmN7FaySFlEJ1QttXzrQfhJVZNdRdYWCiczoiG9PVsznbmBJDs2meVOheTQVwJ57JLZ+T8LCkzNe3X13szZIW45yVozdd5/cljBSAfCEyal5tWxZsuvSffZJtE62KOWWVWLkSfVN2OSA9Gek7Wchd60+jF85cu9jKiqsaQ7cTiKvEvPOFNtXUCWG6gRJFwCnjRvUa9erv7Jd5GLam2BVlcD1fGDM9rkD7cdQGatebDsXGtv5ET8IgtGhryrHHTCpgoCk4goCTUqz/HmiBsyfNZj7COBE4EhJG+vHGuD9wFIqxfBGSR+o591f0hX1ZrYD40f8MeBTXR7xgyCYPcZ/WTZ5zAIL6lMq0G8FgZfZftv4k1qCr2GanKy1f2vqh0/y2G77TmDNwPNOj/hBEIwGHXoF9MHQKggsrP25tgFIWgJkNEhBEATTM6qCdaCCwEuoDobH2R4rnaeJYP048CVJH6ayzJ8MXFB6oSAIAhjNE6uka22/UNL9VHJOA6/Z9l4l8zXJFfAXkm4Cjq4v9g7bVxaue9bIWe9T1t2cxTYVelriFZAysgNpS30qTDXnFbB5c7M2SHsL5JRYqU1kLP3J+NMC5diey5ZPasvl0E1dKnfPSrw+UuTub+p6KQ+AXBhzqm/JurqoYFsyftS9ArrC9gvrf9Px2IU0VcqOAdttf1HS7pKW2r6/iwUEQbBzkXKLnG0kfcz2iZJOs31e2/maeAX8FnAx8MG6aQVVuFcQBMGM2LGj2WOIPFvSgcDJtefTPoOP0smanFhPpYrdvx7A9rclTfYKD4IgaMAo6liBD1CFrv4MVdKpQaWP6/bGNBGs22w/rFq5JGkX2meQCYJgJ2bUBKvt9wLvlfR3tn+n7XxNBOtXJL0NWCLppcCbgX9qe+E+SNUoyxkhUu1d9E0ZLRZuz1TRfTAReJZqy+VYvT+h5s6FtP7oR5PbcolPUxaWnFIsZc3JxfAm2pWw7O26a7ryT0mO1ZKCeak5UgZLaJ5vN2e86iL3alNy18q9Z03m6Cokd9QE6zhdCFVooGOlyrm6GbgZ+G0qp/0pgwOCIAhyjKsCRknHKmm8pt/9ku6r/x1/JFx1pqaJu9UOSZ8BPmM748sTBEHQjFH0Chiau1WdROVsqph91U2PAu+zHaWvgyCYMaOqCpC0G/Bq4CAG5GOpzJtKFXA6VSKV59je1/Y+VMmmj5D0B8UrDoIgqBk1VcAAn6WqVrId+J+BRxFTqQLeALzU9pbxBtvfq4sMXgX8denFgiAIRtTdapyVto9pO8lUgnXRoFAdx/bmOoH1yCFNtsaWVFPNWYJT7Tmrr7YnsgjnrO+p2M2SsqVNvQoAHnpoctu2jLdCKm4xt4bUvLk1pOZIvAe77p72CkhZs3PeGSmLei4cs2Teph4iJdb/YesbS8J9JwrA+e4VAHxV0s/ZvrnNJFMJ1kzE9rSvBUEQZBnxE+sLgTdK+h6wjcq+ZNuHlEwylWA9NONmIKoy2EEQBDNi1LwCBnhZF5NMVf66wI04CIKgGSN+Yn11ou1eSettb2w6Se9VWoMgCCYywoJ1df0Yjy79FaoafG+S9Gnbf9Fkkj6LCa4CPgrsB+wA1to+T9JrgD8Bng481/YNmfGnAb9FpXr4e9vv6Wut05EyTiiXLiH1Gyf3u6dp31yC0lR77loln+SUhaIk8WkJiXXlDEepSNtUxdHMtEUVe3M07ds0v+lUfdvS9tZ0NcdERvzEui/wLNsPAEg6myq734uokrPMrmCl8gM7w/YGSUuB9ZLWAbcAr+KnaQgnIemZVEL1uVSGsi9Iutz2t3tcbxAEQ2KEBesBPNY4/whwoO2HJGXcaCbTm2C1vQnYVP//fkljwArb6wA0td/G04HrxstjS/oKcBwNvy2CIBhdRjGkdYB/BK6T9Nn6+cuBiyTtAXyz6SRD0bFKOgg4nDqnawNuAd4paV/gIarqrTmVwSnAKQArVx7QdqlBEAyBUT2x2n6HpH+mijoV8KYBdeXrms7Tu2CVtCdwCXC67UZZYmyPSXo3sA54ALiRSrWQ6rsWWAtw6KGrI09sEIw4I65jpRakyYNcU3pQTf+UOkLrEuBC25eWjLV9vu1n2X4RcDcQ+tUgmCeMWq6ARNrA8Uc/aQNnSp0d63xgzPa5Mxj/BNt3STqAytj1/CbjJt6MEqtmSQLg3u5623lzJvVU8ulcXG6qrOwe6TDTohKlKb16Yr8LF6R/eCxePHl8TleXuu8554oSUrcnFT5bkqS6bUJraG+9bxqqOpeqtDb1QILu0wb2eWI9AjgROFLSxvqxRtJxku6gEpSXS7oSQNL+kq4YGH+JpG9S+ZOdavvHPa41CIIhMcRE1+MeSNdM11HScyTtN/D8DZI+K+m8vooJzgjb1/LYglyDXJbofyeVkWr8+S/0tLQgCGaRYXkF2B6DaT2QxvkgcHTd/0XAu4DfAw6jsuEcX3LtiLwKgmDoFJxGl0ka/Am/tjZYd81C2+MF436tvs4lVL+cG4eyjhOCNQiCoVMgWLfYXp17UdIXqaI7J3KW7c8m2nMslLSL7e3AUdQunDXFcnLeCdaJSvySfKydhO+VJLtsWvIzZ1VL9c1VSE2R+4mUMlTtvXe679KErj9l/IK0oSv13mR+J+666+TEuCXbzVHys7QPo0vJnE3zpnZFkUG3IV26W9k+upuZuIiqIvUWKt/5fwWQ9GTg3tLJ5p1gDYJg9Bk1P1bb75T0JeCJwFX2T/wfFlDpWosIwRoEwVAZlvFK0nHA+4DlVB5IG23/cn5dvi7R9q2ZXDsEaxAEQ2VYkVe2LyPhgTQMQrAGQTB0Rk0V0DUhWIMgGDohWOcQ0mRDey66smm1TchZYjPm2VTn3CKaWvVTlndIV0jNkYrnLAl/3WuvdN+Ut0DOKyC13wJXjJLKuk1DT3N9SxJVp+gieXaqvYsqqamw1Ny8E9/frq4fgjUIgqBjQrAGQRB0yIgnuu6EEKxBEAydOLEGQRB0SOhY5xgLFky2m+TsM7n2FKmKrM4k7lLJxCkj0Z57Tm7bujU9PvV7KmcoS/UtMV7lDFIpo1bO0JXaW+paBZafXMhlaopc31Sl1y4qujYd39RwBKUG1jQlQm3i9bqq2hqCNQiCoEPixBoEQdADIViDIAg6JLwCgiAIemC+n1h7q3klaZWkqyWNSbpV0ml1+2vq5zskTZXA9g/qfrdIukhSB5k3gyCYbYZY82rW6PPEuh04w/YGSUuB9ZLW8dMCXx/MDZS0Avh94Bm2H5L0KeAE4CNTXXDBgslGZu3IxTFmJkiQ8gDI3/TJfRfsko67VMoinvqNlLOyp0glk4Z0PGfOxJta15Il6b6pcNuU9R/SngWJaz2yPW0mT701uZ+UJSGtbWmb2zxn6U85eOScPsIrYLTos5jgJmBT/f/7JY0BK2yvg0YFvnYBlkh6BNgduLOvtQZBMFxCsHaApIOAw4Hrm/S3/UNJfwn8gKpMwlW2r8rMfQp1fZoDDjigi+UGQdAjO4Pxqjcd6ziS9gQuAU63fV/DMY8HXgkcDOwP7CHp9am+ttfaXm179fLly7tadhAEPbEz6Fh7FaySFlEJ1QttX1ow9Gjg+7Y3234EuBR4QR9rDIJg+Mx3wdqbKkCVEvV8YMz2uYXDfwA8T9LuVKqAo4Abph5SX3eisaqDEpipKUp+yuQU/oua5mPNGYNS5MqWphac03OnDGC5eVMGqUz4qxdPNoCl0sSmQkwhbXzK3d7UdlN5SHOUhI6mQmVzeWJT43PXSs072yGtXeRjLV3DXKTPE+sRwInAkZI21o81ko6TdAfwfKoCX1cCSNpf0hUAtq8HLgY2ADfX61zb41qDIBgSO4MqoE+vgGtJ+R5VTCrwZftOYM3A87OBs/tZXRAEs8lcFppNiMirIAiGys7gFRCCNQiCoRMn1rlGUy1+ot+jO5pH/HTxwUhFGC3KGYl0sx5qAAALnUlEQVRSpCwZKWsQpBc85PCg1PuYMlTlIqRK7kPbHKu53K1Nc6SWFAjsy/BUUhAxt4Y+BGCkDQyCIOiBEKxBEAQdEifWIAiCHgjBGgRB0CHD8gqQdA7wcuBh4LvAG23f0/+Vh5ArIAiCYCJDChBYBzzT9iHAt4AzW8/YkHl1YrUnW9oXLMiYdwvCVFMW5pLwyJIwwO0Jc/Sui/dI9l3YNJ8rlJmCS8zUJWbuhn8oXVi+U+SWlbo/ubS2KYeHlAdBzqugbU7YXLhvitxntI1XQBc/4YelY52QEe864Pj+r1oxrwRrEARzg1nQsZ4MfHJYFwvBGgTB0CkQrMskDSZgWmv7J3lDJH0R2C8x7izbn637nEVV0eTCma22nBCsQRAMlULj1Rbb2dp4to+earCkk4BjgaPsEgVeO0KwBkEwVIalY5V0DPAW4MW2H+z/ij9lXglWG7Zte2xbzoiQImdYSH3P9WVgmbh+gK1b03132WXy5lJtWXLrKol+TbUXhFL2dYZI3ffc/U31ze031Tdn6EqRmjf3+Uitt8Qg1YVL08R5u7pfQ9Kxvh/YDVhX19i7zvabhnHheSVYgyAYfYboFfDk/q+SJgRrEARDJyKvgiAIOiYEaxAEQYfsDImuewtplbRK0tWSxiTdKum0uv019fMdkpJuFJKeNlAna6Ok+ySd3tdagyAYHlHzqh3bgTNsb5C0FFgvaR1wC/Aq4IO5gbb/AzgMQNJC4Ick6mSlxz72eUkIYc7i2TbRdUnC4ZKQxRJrdtso1Uzu6llPmlyyr5I5SrwCms7ZBX1Y+qG/9ZasYT7RZzHBTcCm+v/3SxoDVtheB6DmAfRHAd+1fXsvCw2CYOiEYO0ASQcBhwPXz2D4CcBFU8x9CnAKwKpVB8xg+iAIhsnOkOi69x8AkvYELgFOt31f4dhdgVcAn871sb3W9mrbq5ctW95usUEQDIXQsbZA0iIqoXqh7UtnMMXLgA22/7vblQVBMFvsDF4BvQlWVUrU84Ex2+fOcJrXMoUaIH3dGV6J/DdkifGq5Fu2achiSchjzrhSYuhKUfKHkLsHbUODi1K/JuYtqdKamzf1nvdVZTU1b65Ybur+5Pq2Sc3b5u9rujXMJ/pUBRwBnAgcOeA2tUbScZLuAJ4PXC7pSgBJ+0u6YnywpN2BlwIzOekGQTCihLtVC2xfC+S+3ya5Ttm+E1gz8PxBYN9+VhcEwWwyl4VmEyLyKgiCobIzeAWEYA2CYOiE8SoIgqBD4sQ6x5Bg0aJmfVM3tiT8tW3oaa695Js8ZcnNXStnIW46b4nHRI62f0wl1vum47voW+KBUELJfkvub5s1hFdAM+aVYA2CYPSJE2sQBEEPhGANgiDomBCsQRAEHRIhrXMMKZ83dCKpG1tys0uqv5YYQlLrL/l2z+2/xLhRYowpCfdtavArMcCV5J8t6Zsz0pRUf02RCustMQjl7mNbw16OPoxXoWMNgiDogRCsQRAEHROCNQiCoENCFRAEQdADIViDIAg6JLwC5ho22j7B9JwxjS5YMNm8m7O4piztfVUH7SsRcqq9xJqd+0PYtq35vE3/mErCNru4ZylLd4nXR4qS97bttSD9PohM2eGSD9n2CX1zpYwLiRNrEARBh4SONQiCoAeGIVglvQN4JbADuAv4jTqhfu/0XqU1CIJgkCGWZjnH9iG2DwM+D/xx6xkbEifWIAiGzjCMV7bvG3i6B+SUzt0jd6SMHgUkbQZun6LLMmDLkJYzTGJfc4+5urcDbS9vM4GkL1DtvwmLga0Dz9faXltwrXcCbwDuBV5ie3PjhbZgXgnW6ZB0g+3Vs72Orol9zT3m896GiaQvAvslXjrL9mcH+p0JLLZ99jDWFaqAIAjmLLaPbtj1H4HLgaEI1jBeBUEwL5H0lIGnrwBuG9a1d7YTa2PdzBwj9jX3mM97GxXeJelpVO5WtwNvGtaFdyodaxAEwTAIVUAQBEHHhGANgiDomHkhWCX9gaRbJd0i6SJJiwdee5+kB6YYe6ak70j6D0m/PJwVNye1N0nnS7pR0k2SLpa0Z2LcIkkXSLpZ0ljtbjIyzHRf9dhDJH2tHn/z4P2ebdrsqx5/gKQHJP3hMNcddIztOf0AVgDfB5bUzz9FFRMMsBr4GPBAZuwzgBuB3YCDge8CC2d7T9PtDdhroM+5wFsTY38d+ET9/92B/wQOmu09dbCvXYCbgEPr5/uOyj1rs6+B1y8BPg384WzvJx4zf8yLEyvVH9sSSbtQCZE7JS0EzgH+aIpxr6QSPttsfx/4DvDc3ldbxqS9uQ7VkyRgCelQPQN71OOWAA8D9yX6zRYz3dcvATfZvhHA9o9sFyTp652Z7gtJvwp8D7h1SGsNemLOC1bbPwT+EvgBsAm41/ZVwO8Cn7O9aYrhK4D/Gnh+R902EkyxNyR9GPh/wM8C70sMvxj4n3rcD4C/tH33MNY9HS339VTAkq6UtEHSVF+cQ6XNviTtAbwFePvQFhz0xpwXrJIeT3XyPBjYn+qU9gbgNaT/MB8zPNE2Mv5nmb29HsD2G+u2MeDXEsOfCzxa9zkYOEPSzwxj3dPRcl+7AC8EXlf/e5yko4ax7uloua+3A39tO2sPCOYOc16wAkcD37e92fYjwKVUH9InA9+R9J/A7pK+kxh7B7Bq4PlKYCj5GhuS2tsLxl+sfwJ/Enh1YuyvA1+w/Yjtu4B/o9I5jwJt9nUH8BXbW2w/CFwBPGsIa25Cm339PPAX9ef1dOBtkn63/yUHfTAfBOsPgOdJ2r3WYR0FnGt7P9sH2T4IeND2kxNjPwecIGk3SQcDTwH+fWgrn57U3sYkPRl+orN7OelQvR8AR6piD+B5mX6zQZt9XQkcUo/dBXgx8M0hrXs6Zrwv278w8Hl9D/Dntt8/vKUHXTLnQ1ptXy/pYmADsB34BlOEC0p6BbDa9h/bvlXSp6j+MLcDp46SIWSKvf2LpL2oVBk3Ar8Dj90b8DfAh4Fb6n4ftn3T8HcxmTb7sv1jSecCX6dS21xh+/LZ2MdEWt6vYB4RIa1BEAQdMx9UAUEQBCNFCNYgCIKOCcEaBEHQMSFYgyAIOiYEaxAEQceEYA2ySHpU0sY6U9M/Sdp7ir6HS7IKMoRJeo+kFyXaf1HS51us+4t1FFQQzAohWIOpeMj2YbafCdwNnDpF39cC19b/ToukfYDn2b6m/TIn8THgzT3MGwSNCMEaNOVrZBLU1BFFx1OlyPulhvlRjwe+MDDHMZJuk3Qt8KqB9j0kfUjS1yV9Q9Ir6/bdJX2qznH6SUnXSxoP2f0cDQV8EPRBCNZgWuoUjEdRCawUR1DFyH8X+DKwpsG0RwDr6/kXA39PFe75Czy2TvxZwL/Yfg7wEuCcOkT3zcCPbR8CvAN49vgA2z8GdpO0b9M9BkGXhGANpmKJpI3Aj4B9gHWZfq8FPlH//xM0Oy0+Edhc//9nqQTzt12FAn58oN8vAW+t1/FlYDFwAFVmq08A2L6FKvn1IHdRZZMKgqEz53MFBL3ykO3DJD0O+DyVjvW9gx3q0+yrgVdIOosqHn5fSUtt3z/V3FRCcpxcbLWAV9v+jwnXTaV8HGRxfY0gGDpxYg2mxfa9wO8Dfyhp0YSXjwZutL2qzs50IFV5kV+dZtoxqtSOUGV7OljSk+rngyfeK4HfGxekkg6v268F/nfd9gzg58YH1H33oypHEwRDJwRr0Ajb36DKzHQCQP3THCoheNmE7pdQ5YNF0hWSUj/JLwd+sZ57K3AKcHltvLp9oN87gEXATZJuqZ8D/C2wXNJNVJn3bwLurV97NnCd7e0z2mwQtCSyWwWzRi1Ej7V9zwzGLgQW2d5an3S/BDzV9sOSzqMqy/OljpccBI0IHWswm5xBZYgqFqxUhfqurlUTAn7H9sP1a7eEUA1mkzixBkEQdEzoWIMgCDomBGsQBEHHhGANgiDomBCsQRAEHROCNQiCoGP+P4PxOlE8RTI2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10839b7f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_residuals(resmap._resmap)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding model components\n",
    "\n",
    "The residual map shows a cluster of positive residuals on the South-West of the Crab nebula. Is this another source? To test this hypothesis add a point source at the position that you can guess from the residual map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "newpntsrc = gammalib.GModelSky(gammalib.GModelSpatialPointSource(83.7,21.9),\n",
    "                               gammalib.GModelSpectralPlaw(1.e-17,-2.,gammalib.GEnergy(3.e5,'MeV')),\n",
    "                               gammalib.GModelTemporalConst(1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You have defined a sky model object that has three components:\n",
    "\n",
    "- a spatial model, which is a point source at the position guessed from the residual map;\n",
    "- a spectral model which is a power law with spectral index 2, and a flux which approximately 1/10 of the Crab nebula;\n",
    "- a temporal model which is a constant.\n",
    "\n",
    "Name this source `Src1`, and free its position (fixed by default when a new source is created), so that you can fit it to the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "newpntsrc.name('Src1')\n",
    "newpntsrc['RA'].free()\n",
    "newpntsrc['DEC'].free()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally append the new source to the model container"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<gammalib.model.GModelSky; proxy of <Swig Object of type 'GModelSky *' at 0x1034d4960> >"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "like.obs().models().append(newpntsrc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and fit the model including the new source to the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "like.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Does the addition of the new source provide a better fit to the data? You can quantify this using the test statistic (TS) given by twice the log-likelihood difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "189.90216831723228\n"
     ]
    }
   ],
   "source": [
    "like2 = like.opt().value()\n",
    "ts    = -2.0 * (like2 - like1)\n",
    "print(ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "TS is expected to be distributed as a $\\chi^2_n$ with n degrees of freedom, where n is the additional number of degrees of freedom in the model including the new source, in our case 4 (RA, Dec, Prefactor, and Index). The integral of the $\\chi^2_n$ from TS to $\\infty$ is the chance probability that the likelihood improved by that much due to statistical fluctuations. A large value, like the one we got, means that the chance probability is very low, thus we are likely to have found a new source.\n",
    "\n",
    "To make sure the new source improves the data/model agreement look again at the residual map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVYAAAEKCAYAAABJ430PAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJztnXmUZVV97z/fHqAZVaB9SHczRNToMgzaErSJKGCCHcSg+IJRJJIXYiQJZJEVRd4KMcYsDQkRNYl2gopKQGUQI0RoDUqIQqTbZrKIigaC9AvdYQ70RH/fH+cU3K7ap+rsOkPdqv591rqr6u579j57n3Pqd3f9RtkmCIIgaI850z2BIAiC2UYI1iAIgpYJwRoEQdAyIViDIAhaJgRrEARBy4RgDYIgaJkQrEEQzEgkLZD0b5JulXSnpPdP95xGUfixBkEwE5EkYBfbj0uaD9wInGH7pmmeGvOmewJBEARTwcWu8PHy7fzyNRQ7xc4Eq6QlwGeBvYGtwArbF0g6D3gDsAm4G3in7YcrxpgL3AL81PZxk51zr7328v777Vd3gvWOA1Kb+j43+lXnylhC1rFbt9Zrg7xr0/Sa5axhTkLJlWrLHbcuOWutOn+qvZW5NrgR/3HPPaxfv77RLI6VvL7msavgTmDDQNMK2ytG35QyYhVwIPDXtm9uMre26HLHugU4y/ZqSbsBqyStBFYCZ9veIunDwNnAeyrGOAMYAXavc8L999uPW7797W0bq/6a5o1f+lNb08/Lli3j26oEzVNPjW+bO7f+salnvupcqaUlllXZXjXuhg312gA2barXBnnXMUVqDVXrXbCgXlvuuE3vWYqq52P+/PrzSj0L2pqYLOR9c45h6ateVeu4iVhPsVuqg2CD7aVVn9t+CjhE0rOBKyW91PYdjSfZkM6MV7bX2l5d/v4YhYBcZPs626N/YjcBi1P9JS0Gfhn4+67mGATBNDFnTr1XTcr/er8JHNvVlHPoxStA0v7AocDYbfqpwD9VdPsI8IcUaoSJxj5N0i2Sblm3bl3DmQZB0DlSsfWu85pwGC0sd6pI2gk4BrirhxVMSueCVdKuwOXAmbYfHWg/h0JdcHGiz3HAA7ZXTTa+7RW2l9peunDhwhZnHgRBZ7SzY30ecL2k24DvAittf7XzudegU6+A0gXicuBi21cMtJ8CHAcc7bS/1zLgeEnLgQXA7pI+b/vtXc636j7usEOzcavUV6n2lB4ypdeDtG4uRweXo7vNMfxMstHYhtR6N2+u37+KlJ434z/LLHL0qQHFQ9PCzbB9G8V/wkNHZzvW0sfsQmDE9vkD7cdSGKuOt/1Eqq/ts20vtr0/cBLwz10L1SAIeqRlHeuw0eXMlwEnA0dJWlO+lgMfB3YDVpZtnwCQtI+kazqcTxAEw8IsF6ydqQJs3wik/JeSwtP2/cDyRPs3Kax9QRDMBlpSBQwzEXkVBEG/jHoFzGJm9+qCIBhOYsc6wxh7wzJuoKrCjBtEqgDMrZrDnIqQm2kmx9KfWlrK0g/pS5ay3ledK8f6XvdcOf2h/uOUIzeqIkwbPnbMqXi+suJRu3J5CMEaBEHQIqFjDYIg6IAQrEEQBC0SxqsgCIIOiB3rDKfiBjpDha+ceNBU+xA8RE1tEDkbjJwQ4JShq8rIlAp1rQr3rRsuXEXV9Updh5zbm/Mo5cw3J4VkpTG1L0LHGgRB0AEhWIMgCFokdqxBEAQdEII1CIKgRcIrIAiCoANixzqDSOhuqqz/OYXt5s0bP8bcjG/cyjkkrN85luAUTcM+Ia84XirZdlVxvLqXrOoapNqrihxu3Di+rcqDoE9ywlRzZE/q2lY+S6nnuepkXYS0ho41CIKgA0KwBkEQtEwI1iAIghYJ41UQBEHLhI516khaAnwW2BvYCqywfYGk84A3AJuAu4F32n54TN8FwA3AjuUcL7N9bq0Tj7lhOZGnVdVB01GqaYNUTn7SVHuqreoZTBljcmwQVaGjbYR+pkhtUubPG28pm1/xVG7ekpVJdBw566raUNUNac157towOOYY5tJjpK/tOCNtKlnvVJjlgrXL1W0BzrL9YuBw4HRJLwFWAi+1fRDwA+DsRN+NwFG2DwYOAY6VdHiHcw2CoE+imODUsL0WWFv+/pikEWCR7esGDrsJODHR18Dj5dv55asiz3oQBDOK7UAV0MvqJO0PHArcPOajU4F/qugzV9Ia4AFgpe2xfUePO03SLZJuWbduXXuTDoKgO2b5jrXzmUvaFbgcONP2owPt51CoCy5O9bP9lO1DgMXAYZJeWnHcCttLbS9duHBh+wsIgqBdRr0C6rxmKJ3OXNJ8CqF6se0rBtpPAY4Dji7/7a/E9sOSvgkcC9zR4XSDIOiLGbwbrUOXXgECLgRGbJ8/0H4s8B7gSNtPVPRdCGwuhepOwDHAh+uc96mt21otc6yzVSI+ZT3PeS5yvAJSVBZ5bZg0uer8qfVWWZ3nz0+3p0iFuqa8K+bOSd+InGTOqWTbOdcx59gUXXkF5JDjKVDpTTLm/ky8DapJSzrWKs+jxgO3QJc71mXAycDtpa4U4H3ARyncqFYWspebbL9L0j7A39teDjwPuEjSXAp1xRdtf7XDuQZB0Cft7FhHPY9WS9oNWCVppe3vtzF4E7r0CriRtHPcNRXH3w8sL3+/jcLYFQTBbKOlHWuV5xEwewVrEARBJS3rWCfwPJoWQrAGQdAvebkC9pJ0y8D7FbZXbDtc2vNoOplVgnXr1uowzbGkFPs5Boec8MiqY+saLaqewVT/5qGN6TGqxs3J3Vr3Oqby30J6k5Oz3qrrmBo3J3Kz6bNURY5RLYec+zuWVoxXkLOQ9baXVn1Y5Xk03cwqwRoEwQygPa+ApOfRMDC7ncmCIBhO2om8GvU8OkrSmvK1vPvJT07sWIMg6Jf2vAKqPI+mnRCsQRD0zwwOV63D7F5dEATDx3aQ3WrWCdYurJs5YYhdhCw2Da/MHTcVeppzvaq8IFL3Ju0VkO6fs96cY5t6AOSEC6eommvqOlTNtaoyborUGlqz9tclBGsQBEGLxI41CIKgA0KwBkEQtEwI1iAIghaJ8tczj7pK+K4MUl2ELDatGFo1rzbWkGMIqRtKmWPMyTF05WyScq5N3TDqqjl0FaZaRe+GqrGEjjUIgqADQrAGQRC0TAjWIAiCFglVQBAEQQfMcsHa2eokLZF0vaQRSXdKOqNsP0/SXZJuk3SlpGfX7RsEwSwgyl83IlnoC1gJnG17i6QPA2dTVG2dtO9kRcLs8VbbHAt3G9Uym4ZdpqqLptqguVdAlSU5p8JpKnRz8+b0sXXnVUUbSa3rzqENj4mmdPXsduWZkEXsWKeG7bW2V5e/PwaMAItsX2d79M/xJmBx3b5dzTUIgh4Z1bE2z8c6tPSy156g0NepwBem2Hf089OA0wCWLNm30TyDIOiJGSw069D56qoKfUk6h+Jf/otz+w5ie4XtpbaX7rnnwnYnHwRB+8SOtRlVhb4knQIcBxxtp+NAhrVIWBAELTCDDVN16Gx1VYW+JB1LYaw60vYTOX3rUNd41VWO1BxSz1Yqr2ZOiGcO8+en21NfdTk5Yatyg9YNpWzj2naRF7dqjK7yxOZUB85ZW47Rc2x7K+Gw24Efa5erqyr09XFgN2Bl2fYJAEn7SLpmkr5BEMwGQhUwNSYo9HVNog3b9wPLJ+kbBMFMZzvYsc5uRUcQBMNJCNYgCIKWCcE6c7DHK/f7zJs6UXvTY1M0jczpM4oI0kaTVDRVTmG8KrpaW918uX0b4JpGu9XNa9ua8Sq8AoIgCFokdKxBEAQdEII1CIKgZUKwBkEQtMh2oAqY3asLgmA4aSlAQNKnJD0g6Y4eZl2bWbVjzcnHmpPbsyk5OVJTx1atoSq8se65qvqn2tsIpazrXdFGpdmm97fKM6ELr4CZlvu1Me16BXyGIprzs20N2AazSrAGQTADaFEVYPuGMrVoq0jaBdhge0pf0SFYgyDon/qCdS9Jtwy8X2F7RdvTkTQHOAl4G/AKYCOwo6R1FGH4K2z/sO54IViDIOif+oJ1ve2lXU6l5Hrg6xSlou6wvRVA0h7Aa4EPSbrS9ufrDBaCNQiCfhlOr4BjbG+WdK7t20YbbT9IkRf68jJHdC1mlWBN3a8qpXxOaF5TQ0jVHOo+W037Q3ODVE7IZI7hJ8d4lZpvzv1tw0BT1+aijNxsbYTw5hg9p72Y4BCGtNoefcLPlbQzsAewGrjU9kNjjpmUofvaCIJgO6A9d6tLgO8AL5J0n6TfaDgzAxuAa4ElwLclHZw7yHB9bQRBsH3QnlfAW1sZ6Bnusn1u+ftlkj4DfAI4KmeQ2LEGQdAvw11McL2kl4++sf0DILtKaexYgyDon+EzXo3ye8ClklYBtwMHAT/JHaSz1UlaIul6SSOS7pR0Rtl+nqS7JN0m6UpJz67oP5ShakEQNGSId6y2bwUOAS4pm64HstUNXe5YtwBn2V4taTdglaSVwErgbNtbJH2Ywm/sPYn+n2EKoWpN7kWOFTUnnLTKQlzXGpyTVDsn/LXK0p/jMZFaW05V2abG4TbCclPkzCvHyt7UWyAn3DeHuuPmzH9ChswrQJLs4sm3vRG4unwlj5mMzr4SbK+1vbr8/TFgBFhk+zrbo4/4TcDiiv43AA92Nb8gCKaJ4dyxXi/pdyXtu+1UtYOkoyRdBJxSd7BevjbKWN5DgZvHfHQq8IU+5hAEwRAxfDrWYynk0SWSDgAeBnai2HxeB/yV7TV1B+tcsEralSJy4Uzbjw60n0OhLri44finAacBLF687yRHB0Ew7Qxh5JXtDcDfAH9TRljtBTxp++GpjDepYC2TExwM7AM8Cdxp+7/qDF5O8HLgYttXDLSfAhwHHF1XZ1FFmZBhBcAhhyxto9RZEARdM2SCdZAywmptkzEqBauk51MYlY4BfgisAxYAL5T0BPBJ4KLRZAWJ/gIuBEZsnz/Qfmw57pG2n2gy+fHnHK8TzzFYVNG0GmoVdb9S2jBebdpU/9icKp4pA0vVsU2NV03vw8aN6WNT96HKSLPTTlM/P8AOO4xvqzJizq8dmT7DjFdDGNLaNhN9bfwp8Hng+bZ/yfbbbZ9o+yDgeOBZwMkT9F9Wfn6UpDXlazmFpX83YGXZ9gkASftIuma0cwehakEQDAvDZ7xqlcqvjYlCxWw/AHxkooFt3wikvt+uSbRh+35geZ3zB0EwgxlCHeso5X/abwN+xvaflF4Ce9v+t5xx6uhY35RofgS4vRSwQRAEeQypYKUwYG2lyA3wJ8BjFHaiV+QMUkfR8RvAKykiEABeQ+F/+kJJf2L7czknDIIgGGLB+vO2XybpewC2H5KU0IxPTB3BuhV48agngKT/Bfwt8PPADUAI1iAI6jPEqgBgs6S5FOkDkbSQQgZmUUew7j/GveoB4IW2H5SUkf64e1L3q437l1MdtKl1NidENOdcOePWrURa1V51bM4c6s4r59pUJSxPeUxUjVs3gXaOE+GCBfWPraKr6zjWq2Y78Ar4KHAl8FxJHwROBP5v7iB1Vvcvkr4KfKl8fyJwQ1nFcErOs0EQbOcM6Y7V9sVlZqujKYzvv2J7JHecOoL1dOBNwBHliS4CLi8d+1+be8IgCIJhFawAtu8C7moyxqSrKwXoLcDVts+kcJfatclJgyDYjhnOJCzl1HTRYCpTSc+R9KnccSaduaTfBC6jiLQCWAR8OfdEQRAETzOkghU4aDA/QFlI8NDcQeqqAg6jzExl+4eSnpt7oj6wxxsScoxMOVVLc0JlcwwsOUaPHCNEylbQRhXP1LhN8892RU6u2twx6pJzH1LtbRjrmuR5bZbZo2S4vQLmSHrOaGVWSXswhWRVdTpstL1J5V+LpHmUrghBEARTYni9Av4S+I6kUWP9W4A/yx2kzuq+Jel9wE6SXge8G/jH3BMFQRAAQ71jtf1ZSbfwTFXWN9n+fu44dVb3XorMVrcDv0VhvMr26wqCIHiaIdWxStqRoubV7sAewImS/ih3nEl3rGVawL8rX0EQBM0Y4h0rcBVFLpRVQEWiycmZKB/r7UygSy3TBwZBEOQzvIJ1se1jmw4y0Y71uPLn6eXP0ZwAbwNaTVDdFvb48MRUuCLkVfFMVTPNsY421dNXPYNNkybnhLRW0UUS8BwreVMrPeQlDK/rIZLjgVB1bJVVv89xx9KaV8DwGq++LennbN/eZJCJ8rHeAyBpme1lAx+9V9K/UqTUCoIgyGd4d6xHAL8u6ScUqgBRxEll/Yde52tjF0lHlImrkfQqYJfc2QZBEACt6ljLUk8XAHOBv7f9oYZDvr75rOrnY/2UpGdR6FwfoSgTGwRBMDVaEKxler+/Bl4H3Ad8V9JXpuIeNYrteyQ9B3gBRY2/Ue7JGaeOV8Aq4GBJuwOy/UjWTIMgCAZpb8d6GPAj2z8uhtWlwBuBKQtWSf8HOANYDKwBDqeovXfURP3GMpFXwNuBfxitwmr70TGfPx943qiKINF/CfBZYG+KRLErbF8g6TzgDcAm4G7gnana3VPZ4m/dChs2bNtWZZBKtaeMVFA/Byc0f15ycpamDFW7VChptOHJ8Y0bKix7KSomMTdlhEhZ1QDPGz/hnPvQRThp1bhtGMVSNM1nWmV4ygm7bhLC29p1qf+HslfpsD/KirLkPRR5S/5z4LP7KBLwN+EMijIsN9l+raSfBd6fO8hEO9Y9ge+VuQlX8Uz56wOBI4H1FMEDVWwBzrK9WtJuwCpJK4GVwNm2t0j6MHA2RTnsp+liix8EwfDw1Nba3zDrbS+t+Cw1SFO/hQ22N0hC0o6275L0otxBJvIKuEDSxym2wMuAg4AngRHgZNv3TjSw7bXA2vL3xySNAItsXzdw2E0UibPH0voWPwiC4SCVLGmK3AcsGXi/GLi/6Zhl2sAvAyslPTSVMSfUsdp+imKHuXJKUyyRtD9F6q2bx3x0KvCFRJfaW3xJpwGnASxatG+TaQZB0BMtCdbvAi+QdADwU+Ak4NeaDGj7hPLXP5Z0PfAs4Gu543TupStpV4rysWcO6mklnUOhLrg41S3Rltzil/qWFQAHHbQ0sm4FwZDT1o61VCf+DnAthS3mU7bvbD7y0+N/a6p9OxWskuZTCNWLbV8x0H4KRWTX0WWFgrFMaYtvj1fM5xiZUtFJVXSVyzQnV2ZqvtpUEd78+OPj256oCKDLSQSaMlRVVMfTzjuPa5u3w47j2jZWLCEn52hOQcS656oidc+qzpVzf5tStYacSMJO8rEmxp0qtq+hSAzVCpIuAs4YNaiXrld/aTvLxbQzwaoigeuFwIjt8wfaj6UwVh1puyo0tvUtfhAEw0NXXhctMK6CgKTsCgJ1SrP8WaIGzJ/WGHsZcDJwlKQ15Ws58HFgNwrF8BpJnyjH3UfSNeVitgCjW/wR4IttbvGDIJg+Rv+zrPOaBuaUu1Sg2woCr7f9vtE3pQRfziQ5WUv/1pSuNLltt30/sHzgfatb/CAIhoMWvQK6oLcKAnNLf66NAJJ2AsYrxoIgCGoyrIJ1oILAayk2hifYHskdp45g/TzwDUmfprDMnwpclHuiIAgCGM4dq6QbbR8h6TEKOaeBz2x795zx6uQK+HNJtwHHlCf7gO1rM+c9beRY5HPGyHkwquaQCm9MzSsrb+qmjDjGqmS1VWb5FDsm/nnJUI4psYi5c9PuGV1UtYXKCNzaY6TaUpcFhjcNad+CbtgEq+0jyp+7tTFe3ds8Amyx/XVJO0vazfZjbUwgCILti5Rb5HQj6XO2T5Z0hu0Lmo5XxyvgN4HLgE+WTYsowr2CIAimxNat9V498nJJ+wGnlp5Pewy+cgers2M9nSJ2/2YA2z+U9NzcEwVBEMBw6liBT1CErv4MRdKpQUWdy/ba1BGsG21vUqkQlDSP5hlkgiDYjhk2wWr7o8BHJf2t7d9uOl4dwfotSe8DdpL0OuDdwD82PXEXpGqUVeW/bFpcL6fYXA45IbG1B2hj4CorUU5VxpyKdwly8qbmROWmqDIypcZIham2YbzKyQPcNAQ3Zw5tMGyCdZQ2hCrU0LFS5FxdB9wO/BaF0/6EwQFBEARVjKoChknHKmm0pt9jkh4tf46+Hp2s/1jquFttlfRl4Mu2101hzkEQBE8zjF4BvblblUlUzqWI2VfZ9BTwMdtR+joIgikzrKoASTsCbwb2Z0A+5sq8iVQBZ1IkUnmF7T1t70GRbHqZpN/PnnEQBEHJsKkCBriKolrJFuB/Bl5ZTKQKeAfwOtvrRxts/7gsMngd8Fe5JwuCIBhSd6tRFts+tukgEwnW+YNCdRTb68oE1kOHVD88MWWdzbEE51hn2wh/rXuuykVkJKTOWnDq2KqbUPOiP5UoKFtFG3+gOQ4TKYt6jodJir71jTnrHXt9m1aZrRp3iPi2pJ+zfXuTQSYSrBPVRs6omxwEQfAMQ75jPQJ4p6QfAxsp7Eu2fVDOIBMJ1oMr3AxEUQY7CIJgSgybV8AAr29jkInKX2e4EQdBENRjyHesb060PSJple01dQcZ0iRmQRDMZoZYsC4tX6PRpb9MUYPvXZK+ZPvP6wzSZTHBJcBngb2BrcAK2xdIegvwx8CLgcNs31LR/wzgNylUD39n+yNTmUcbFTBTCvuqsMCmuVvrnr/qXJXGq5ShquriNAw9zZmD54y/kG3kWO3iPuTQhnGzi/Do6RhjLEO+Y90TeJntxwEknUuR3e/VFMlZplewUviBnWV7taTdgFWSVgJ3AG/imTSE45D0UgqhehiFoexrkq62/cMO5xsEQU8MsWDdl22N85uB/Ww/Kal2BvjOBKvttcDa8vfHJI0Ai2yvBNDEfhsvBm4aLY8t6VvACdT8tgiCYHgZxpDWAf4BuEnSVeX7NwCXSNoF+H7dQXrRsUraHziUMqdrDe4APihpT+BJiuqtVSqD04DTABYv3rfpVIMg6IFh3bHa/oCkf6KIOhXwrgF15dvqjtO5YJW0K3A5cKbtWllibI9I+jCwEngcuJVCtZA6dgWwAuDgg5dGntggGHKGXMdKKUiTG7m6dKCafoYyQuty4GLbV+T0tX2h7ZfZfjXwIBD61SCYJQxbroBE2sDRVzdpA6dKmR3rQmDE9vlT6P9c2w9I2pfC2PXKOv3G3oyu8j5XeQU89VSzcVOq51TIZO64WTG8HeF54xdSVSg2Rd0k01Xk6PVyIng3bx7flvN8VHlBdBUK3SQsta3E130IzboeSNB+2sAud6zLgJOBoyStKV/LJZ0g6T4KQXm1pGsBJO0j6ZqB/pdL+j6FP9npth/qcK5BEPREj4muRz2QbpjsQEmvkLT3wPt3SLpK0gVdFROcErZvZNuCXINcmTj+fgoj1ej7X+hoakEQTCN9eQXYHoFJPZBG+SRwTHn8q4EPAb8LHEJhwzkx59wReRUEQe9k7Eb3kjT4L/yK0mDdNnNtP1j+/qvleS6n+M+5dijrKCFYgyDonQzBut720qoPJX2dIrpzLOfYvirRXsVcSfNsbwGOpnThLMmWk7NOsI5V4ucYr6r+Y8gxkDQlda7K3KBzMsp4JgxVrtDUdPVv2taEoaqNSNkUOWtoatRK3Z+qMbvK19uVMaiLZ79Ndyvbx7QzEpdQVKReT+E7/y8Akg4EHskdbNYJ1iAIhp9h82O1/UFJ3wCeB1xnP+3/MIdC15pFCNYgCHqlL+OVpBOAjwELKTyQ1tj+pep5+aZE2w+mcu4QrEEQ9EpfkVe2ryThgdQHIViDIOidYVMFtE0I1iAIeicE6wxCGm85zikY2kbS5FTIYlOyQlcrDk55AFTpuVLtVetqK8SxDqml5XgKVK2hqaW+KV15rrTxLI4Np26jSuuwJ2Fpg1klWIMgmBmEYA2CIGiRIU903QohWIMg6J3YsQZBELRI6FhnGHPmwM47b9tWZdzIMQg1fQhycnvm9J/bcA05Brgcw09Tcu5NF1VEJ6LuenMMUlXPaFcpdJuEEbd1vUOwBkEQtEjsWIMgCDogBGsQBEGLhFdAEARBB8z2HWtnqn9JSyRdL2lE0p2Szijb31K+3yppogS2v18ed4ekSyQt6GquQRD0R481r6aNLnesW4CzbK+WtBuwStJKninw9cmqjpIWAb8HvMT2k5K+CJwEfGaiE86ZAwvGiF9trR/X5znprL45FvWcEM/UGCnre5UlOBWmKtITaFqVNoem/+a1EVqcas+phprjyZFzbXMs/alw7KoQ7a68XMIrYGp0WUxwLbC2/P0xSSPAItsroVaBr3nATpI2AzsD93c11yAI+iUEawtI2h84FLi5zvG2fyrpL4B7KcokXGf7uoqxT6OsT7Pvvvu2Md0gCDpkezBede5eLWlX4HLgTNuP1uzzHOCNwAHAPsAukt6eOtb2CttLbS9duHBhW9MOgqAjtgcda6eCVdJ8CqF6se0rMroeA/zE9jrbm4ErgFd1MccgCPpntgvWzlQBKpSoFwIjts/P7H4vcLiknSlUAUcDt0zcpTzvWGNVVgnM+iUp28h1mTIE9FkRtorUHHIMP1XUPTbHcJRjvGpqkKpqT12vsXlMJ+pfZbxKjdt3iPbY87WRjzV3DjORLnesy4CTgaMkrSlfyyWdIOk+4JUUBb6uBZC0j6RrAGzfDFwGrAZuL+e5osO5BkHQE9uDKqBLr4AboaJwfaLAl+37geUD788Fzu1mdkEQTCczWWjWISKvgiDole3BKyAEaxAEvRM71plGXS1+4riqm93Vt2vdyKuqJaXmNWdOfetCGwURc4wxqfk2NX61cc9yjIgpo1TTwpRdGZ7aEF5dCMBIGxgEQdABIViDIAhaJHasQRAEHRCCNQiCoEX68gqQdB7wBmATcDfwTtsPd3/mHnIFBEEQjKWnAIGVwEttHwT8ADi78Yg1mVU7Vhs2b9nWKj6nKkw1cdOqvkVTVvKuQjzregpA8/DXnDX0aeXOCVPNoWpeqTDNHXdMH5vKh5q6Dzn3JscLY/Pm+se2cX/HHtvGv/B96VjHZMS7CTix+7MWzCrBGgTBzGAadKynAl/o62QhWIMg6J0MwbqXpMEETCtsP503RNLXgb0T/c6xfVV5zDkUFU0untps8wnBGgRBr2Qar9bbrqyCXrEaAAAKWElEQVSNZ/uYiTpLOgU4Djjazimc1IwQrEEQ9EpfOlZJxwLvAY60/UT3Z3yGWSVYbdi4cdu2NowIOcaFpsaYpqGYXRUIrApTTbVX5exMHdvU7abq2ubkDc3Je5o6tsrQlSJHoDQ1mjZ9llLjtrXn60nH+nFgR2BlWWPvJtvv6uPEs0qwBkEw/PToFXBg92dJE4I1CILeicirIAiClgnBGgRB0CLbQ6LrzkJaJS2RdL2kEUl3SjqjbH9L+X6rpKQbhaQXDdTJWiPpUUlndjXXIAj6I2peNWMLcJbt1ZJ2A1ZJWgncAbwJ+GRVR9v/DhwCIGku8FMSdbLSfbd9n2PRz7G4thF2mfrWbpoMOidss4qU5buu1biqfxU5lVdzkmo3HbeKumvryjujjerAbTxPfc5hJtJlMcG1wNry98ckjQCLbK8EUP2/9KOBu23f08lEgyDonRCsLSBpf+BQ4OYpdD8JuGSCsU8DTgNYsmTfKQwfBEGfbA+JrjtPGyhpV+By4Ezbj2b23QE4HvhS1TG2V9heanvpXnstbDbZIAh6IXSsDZA0n0KoXmz7iikM8Xpgte3/andmQRBMF9uDV0BnglWFEvVCYMT2+VMc5q1MoAZIn3eKZyKv4mcbxquUISLHUJZDygiRc61y5tDUCJhjMMkxquUYr7qqnJpDzhxS17zKsNck326Tv6/J5jCb6FIVsAw4GThqwG1quaQTJN0HvBK4WtK1AJL2kXTNaGdJOwOvA6ay0w2CYEgJd6sG2L4RqPp+G+c6Zft+YPnA+yeAPbuZXRAE08lMFpp1iMirIAh6ZXvwCgjBGgRB74TxKgiCoEVixzrDkGD+/HrH5lRDTZFTLbMqDDLV3jSRcNW56l4XyLOSp3YeOZb6HJpa77s6tk9Pjqp5parHdjGH8Aqox6wSrEEQDD+xYw2CIOiAEKxBEAQtE4I1CIKgRSKkdYYhpZX4KermQoW0waBpztEqqsIQm/ZPtbcRDpoT7lvXOFh1bXPWkHMdm4b7pmgjb2qO8app3tS6/dswXoWONQiCoANCsAZBELRMCNYgCIIWCVVAEARBB4RgDYIgaJHwCphp2GjLmFjTCnPnvHnjTc85iYFz6Cppck7V0lR7juW66g8hFdpbNW7dP6YcS3+O9b/q2JSlO8frI7XenNDkNiztSc8GcjKOVww89timMdcTTGE2MbsEaxAEQ0/oWIMgCDqgD8Eq6QPAGyn24w8Av14m1O+czqu0BkEQDNJjaZbzbB9k+xDgq8AfNR6xJrFjDYKgd/owXtl+dODtLlCldG4fuSVl9DAgaR1wzwSH7AWs72k6fRLrmnnM1LXtZ3thkwEkfY1i/XVYAGwYeL/C9oqMc30QeAfwCPBa2+tqT7QBs0qwToakW2wvne55tE2sa+Yxm9fWJ5K+Duyd+Ogc21cNHHc2sMD2uX3MK1QBQRDMWGwfU/PQfwCuBnoRrGG8CoJgViLpBQNvjwfu6uvc29uOtbZuZoYR65p5zOa1DQsfkvQiCnere4B39XXi7UrHGgRB0AehCgiCIGiZEKxBEAQtMysEq6Tfl3SnpDskXSJpwcBnH5P0+AR9z5b0I0n/LumX+plxfVJrk3ShpFsl3SbpMkm7JvrNl3SRpNsljZTuJkPDVNdV9j1I0nfK/rcP3u/ppsm6yv77Snpc0h/0Oe+gZWzP6BewCPgJsFP5/osUMcEAS4HPAY9X9H0JcCuwI3AAcDcwd7rXNNnagN0HjjkfeG+i768Bl5a/7wz8B7D/dK+phXXNA24DDi7f7zks96zJugY+vxz4EvAH072eeE39NSt2rBR/bDtJmkchRO6XNBc4D/jDCfq9kUL4bLT9E+BHwGGdzzaPcWtzGaonScBOpEP1DOxS9tsJ2AQ8mjhuupjqun4RuM32rQC2/9t2C6X7WmOq60LSrwA/Bu7saa5BR8x4wWr7p8BfAPcCa4FHbF8H/A7wFdtrJ+i+CPjPgff3lW1DwQRrQ9Kngf8H/CzwsUT3y4D/KfvdC/yF7Qf7mPdkNFzXCwFLulbSakkTfXH2SpN1SdoFeA/w/t4mHHTGjBeskp5DsfM8ANiHYpf2DuAtpP8wt+meaBsa/7OKtb0dwPY7y7YR4FcT3Q8DniqPOQA4S9LP9DHvyWi4rnnAEcDbyp8nSDq6j3lPRsN1vR/4K9uV9oBg5jDjBStwDPAT2+tsbwauoHhIDwR+JOk/gJ0l/SjR9z5gycD7xUAv+Rprklrbq0Y/LP8F/gLw5kTfXwO+Znuz7QeAf6XQOQ8DTdZ1H/At2+ttPwFcA7yshznXocm6fh748/J5PRN4n6Tf6X7KQRfMBsF6L3C4pJ1LHdbRwPm297a9v+39gSdsH5jo+xXgJEk7SjoAeAHwb73NfHJSaxuRdCA8rbN7A+lQvXuBo1SwC3B4xXHTQZN1XQscVPadBxwJfL+neU/GlNdl+xcGntePAH9m++P9TT1okxkf0mr7ZkmXAauBLcD3mCBcUNLxwFLbf2T7TklfpPjD3AKcPkyGkAnW9s+SdqdQZdwK/DZsuzbgr4FPA3eUx33a9m39r2I8TdZl+yFJ5wPfpVDbXGP76ulYx1ga3q9gFhEhrUEQBC0zG1QBQRAEQ0UI1iAIgpYJwRoEQdAyIViDIAhaJgRrEARBy4RgDSqR9JSkNWWmpn+U9OwJjj1UkpWRIUzSRyS9OtH+GklfbTDvr5dRUEEwLYRgDSbiSduH2H4p8CBw+gTHvhW4sfw5KZL2AA63fUPzaY7jc8C7Oxg3CGoRgjWoy3eoSFBTRhSdSJEi7xdr5kc9EfjawBjHSrpL0o3Amwbad5H0KUnflfQ9SW8s23eW9MUyx+kXJN0saTRk9yvUFPBB0AUhWINJKVMwHk0hsFIso4iRvxv4JrC8xrDLgFXl+AuAv6MI9/wFtq0Tfw7wz7ZfAbwWOK8M0X038JDtg4APAC8f7WD7IWBHSXvWXWMQtEkI1mAidpK0BvhvYA9gZcVxbwUuLX+/lHq7xecB68rff5ZCMP/QRSjg5weO+0XgveU8vgksAPalyGx1KYDtOyiSXw/yAEU2qSDonRmfKyDolCdtHyLpWcBXKXSsHx08oNzNvhk4XtI5FPHwe0razfZjE41NISRHqYqtFvBm2/8+5ryplI+DLCjPEQS9EzvWYFJsPwL8HvAHkuaP+fgY4FbbS8rsTPtRlBf5lUmGHaFI7QhFtqcDJD2/fD+4470W+N1RQSrp0LL9RuB/l20vAX5utEN57N4U5WiCoHdCsAa1sP09isxMJwGU/5pDIQSvHHP45RT5YJF0jaTUv+RXA68px94AnAZcXRqv7hk47gPAfOA2SXeU7wH+Blgo6TaKzPu3AY+Un70cuMn2liktNggaEtmtgmmjFKLH2X54Cn3nAvNtbyh3ut8AXmh7k6QLKMryfKPlKQdBLULHGkwnZ1EYorIFK0WhvutL1YSA37a9qfzsjhCqwXQSO9YgCIKWCR1rEARBy4RgDYIgaJkQrEEQBC0TgjUIgqBlQrAGQRC0zP8H6TAUXD3xj8wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1094d9978>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "resmap.models(like.obs().models())\n",
    "resmap.run()\n",
    "plot_residuals(resmap._resmap)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The addition of the new point source has flattened the spatial residuals, even though the fit is obviously not perfect yet. From now on fix the position of `Src1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "like.obs().models()['Src1']['RA'].fix()\n",
    "like.obs().models()['Src1']['DEC'].fix()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Modifying the spectral models\n",
    "\n",
    "Next check if the model spectrum satisfactorily reproduces the data in the region around the Crab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "resspec = cscripts.csresspec(like.obs())\n",
    "resspec['algorithm'] = 'SIGNIFICANCE'\n",
    "resspec['mask'] = True\n",
    "resspec['ra'] = 83.63\n",
    "resspec['dec'] = 22.01\n",
    "resspec['rad'] = 0.2\n",
    "resspec['components'] = True\n",
    "resspec['outfile'] = 'resspec_1.fits'\n",
    "resspec.execute()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use an example script to display the residuals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEOCAYAAACetPCkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl4VOXZ+PHvTUAhoKgsLkQSNrEiIbK+KigIrbgAghuQCFRLymulr61Kpa0F6xuLv0ppXarGFvCFCG7ggq0LFHAprYBQRBFRCYqoQJQIJEgI9++PmYlZZiYzc87s9+e65iJz5pznPJMc5p5nu4+oKsYYY0yomsS7AsYYY5KLBQ5jjDFhscBhjDEmLBY4jDHGhMUChzHGmLBY4DDGGBMWCxzGGGPCYoHDGGNMWCxwGGOMCYsFDmOMMWFpGu8KREPbtm01Jycn3tUwxpiksX79+r2q2i6UfVMycOTk5LBu3bp4V8MYY5KGiOwIdV/rqgpBSUkJOTk5NGnShJycHEpKSuJdJWOMiRsLHI0oKSmhsLCQHTt2oKrs2LGDwsLCsIKHBR5jTCqRVEyr3rdvX420q0pEIj6vv9+lL/BUVFTUbMvMzKS4uJj8/PyIz2WMMW4SkfWq2jeUfa3F4SIRafAoKCioEzQAKioqKCgoqLNfINZaMcYkGgsc9ahqnUd2drbf/bKzsxvs6zan3WQWdIwx0WBdVbXceeedDbZt2rSJF154gaqqqpptzZo1Y8SIEeTm5tbZt1u3bowfP77OtpycHHbsaDhZITs7m9LS0prn1kVmjImncLqqUnI6rpt8wWHFihWUl5fTunVrhg4d2iBoAGzbtq1B8OnTpw+7du1qEHj69OnjN1BFItSg4+siKygoqLM9UOD51a9+xSeffELHjh0pKiqygGOMAazFUcc9b93D+1+9H9E5229tT+a+TL+vbdq0KaTAU9+cOXMoLy9vsL1169b87Gc/q3k+c+bMiOrsU/8asNaKMeknnBaHBY4EMPvR2RzYdaDB9nC6yeoLNeiAs8AT6PqxFosxycW6qpLMLZNvCfhaYx/AgYLO0KFD/QadoUOHulr32Y/OblD/+i0W36A+YMHDmBSQ8C0OEbkCuAxoDzyoqq80dkyytTiiJZRv/f4CT6xaK/cW3xs0aBpjYifhu6pEZC5wObBbVc+utX048CcgA/iLqs6q9dqJwL2qekNj5VvgcCbUMY7Zj87m1sJbIz6Pk6DT6rRWFnSMcVEyLACcDwyvvUFEMoAHgUuAs4BxInJWrV1+7X3dRFl+fj7FxcVkZ2cjImRnZ/sdGL9l8i0N1rKEuvbl3uJ7A55/06ZNzJkzh5kzZzJnzhw2bdrUYB9/3XPGmNiIW1eViOQAy3wtDhE5F5ipqhd7n0/37jrL+3hVVZcHKa8QKATo2LFjH39rJ0xsOJmVFcqxwaYxhzqDzVosxtSV8F1V4DdwXAUMV9UfeZ9fBwwAPgAmAmuBjar6cGNlW1dV/IUyvuJk0aO/bi4ns9DCYUHHpKJkDRxXAxfXCxz9VXVquGVb4EgOTgKHU07XvsyYMcOdihiTIJJ1Ou5O4PRaz7OAXXGqi4kBf19aYpGiBSL/4Pd1k0W66t9aKyYVJFKSw7VANxHpJCLHAGOB5+NcJxNjRUVFZGbWXYGfmZlJUVFRnW2RDso7TUjZ6rRWER8LNqhvUkO8puMuAgYDbYEvgRmq+lcRuRT4I57puHNVtShwKYFZV1Vyi3TVudNUKU5Wu4dyrNPcZNZaMdGU8F1VqjouwPa/AX+LcXVMgsnPz49ohbnvGDeCTjir3UM9ttVprRy1OKy1YhJFwq8cj4S1OExj4jUwH+n/N2utmGhLhgWAxpgw2NiKSSTW4jCG0GdzuX1sLLhx3xdrsaS+hB/jiIeqqip27tzJoUOH4l0Vk4CefvppysrK6nQliQht2rRhy5YtNG/enKysLJo1a9bg2KKiIr+D8vVngsWL07EVsBaLqSttWhzbt2/nuOOOo02bNnFdeGYSV1lZGZ999hmHDx/mmGOOoUOHDrRp0wZVpaysjP3799OpUye/x6by/UdsfCU92BiHH4cOHbKgYYJq06YNubm59O3bl9zcXNq0aQN81/II1lrNz8+ntLSUo0ePUlpaGnbQKCkpIScnhyZNmpCTk0NJSYmj9+ImG18x9aVNVxXEN8WFSW7RvHYS/cZXTloLboyvmMSTVoEj3kSEgoICFixYAMCRI0c49dRTGTBgAMuWLQu5nJycHNatW0fbtm0d7WNiL9QAVFFRQUFBAQUFBTXbkrlb2VK0pBYLHLWV74SqysiPb3EitAz8Qd2yZUs2b95MZWUlLVq04NVXX6VDhw6Rn8+YBGeLHlOTBQ63+AJOkMABcMkll/Diiy9y1VVXsWjRIsaNG8frr78OwFdffcX111/Pxx9/XJMqIzc3l7KyMsaNG8eePXvo379/nW+eCxcu5L777uPw4cMMGDCAP//5z2RkZETtbRpnnCR2TEZudHNZayXxpM3geEhaZ0HbbpE9mrXwBI+92wI/VBl78bksfuwvHNq5mU1vr2XAmR3g8EE4uJcZM2ZwzjnnsGnTJu6++24mTJgAeP7jDBw4kA0bNjBy5Eg++eQTALZs2cITTzzBm2++ycaNG8nIyEioQVUTmlATO6YbG5RPXGnZ4ojWgN2Mn4xvdJ/cHmdS+ulnLFqyjEuHXejZqNVQ+TVvvPEGzzzzDAAXXXQRZWVllJeX89prr7FkyRIALrvsMk488UQAVqxYwfr16+nXrx8AlZWVtG/fPgrvzESTkxxbqSyerRWwFksw1uJwUVnzjsFbJSLQthsjx1zNrXf+nnE3/MTTyhFP15K/bgzfYKq/QVVVZeLEiWzcuJGNGzeydetWxzcoMvHhZDpvSUkJIoKIJNxU3kCc1rmx4522VsBaLMGkZYsj0E18du2rpLKqOqIyD1VVs6+yijatjm103+uvv57WrVvTs2dPVq1a5dlYVckF/XIpefQ+7rjlJ6x689+0PeE4jj/8pWd78X38+pYb+fvy1Xz99ddQ9hFDe3dl1JzZ/KxwAu07fY+vvvqK/fv3B7wvhUk9vqm8Pok2ldcfp3UO5XinLQUbXwkuLQNHIKed0CLiYz/ac4BDVdV8tCfwtxRVz34cewKXj7+Bj/Yc4LN9lRzWDKqaHMvMaVP54U+nk3vhCDJbNOexB+4BYMZtNzHuxz+n90VXcOF5/emYdRoAZ3Xvyv9Ou5EfjBjD0SbNaNasGQ8++KAFjhQVjam8JSUlNftlZ2dHpYsslHr7q3M43J6+bLPBgkublCNbtmzhe9/7XtTOWXbgW/ZVVkV07KGqapo3y6BLuwia13u3eQblm0Ue9BqbRmw8on0NNcbJIkR//8+d3vgqVMmWwt4pXysl2e5LH9UkhyJyInC6qm4Ku2YprE2rY0PqpvInWCulUS1OjPxYgMMHPI/KryM/vwWdmIjFPdoDffN38iFc/1in04+TZfpyKndzhRQ4RGQVMNK7/0Zgj4isVtWfR7FuaaWxbq7AmnNCi44RBy0O7o08aIS4dsVET6Jn5vXHaZ0T/T2nQzdXqC2O1qr6jYj8CJinqjNExFocLjmhRTP2RXjsoapq9kHkgaNl28g/+H3dZHu3RXa8tVYc83UpNTZOkUgLD0Otc7SOj7ZgrQVfyy9QC86NacSx6CILNXA0FZFTgWuAX0WxPmkpbt1cTjnpJrPWimsivUd7PL+5R1pnt46Ph9pThnNycvwGOzfunRILoQaOO4GXgTdUda2IdAYi/Jpp3BZ5N5fHCS2aRRa4rLWS1BL9m3sqCXUKcqKPbfiENKtKRM5X1Tcb25Yo4jGrKl6czOYChzO6nHA6ttKshWdRZQyl6jXUWPeJCZ/bM+BiIRqzqu4HeoewzQRRVFTE448/TkZGBk2aNOGRRx5hwIABIR07fPhw/vWvfzFw4MA6KdiddHNBaOtPgknK1gpYi6UWJx9U6Rh04vmeE+X3HTTliIicKyK3AO1E5Oe1HjOBmKRgFZGWIvKYiDwqIknbhl6zZg3Lli3j7bffZtOmTSxfvpzTTz+9zj7V1YFXrd9222019/Fw0wktmtG8WWR/St9q+ZhrcaKzdStVlZG3dkyN+n32yZDqxKlQ37Oq1nkEWpSbnZ3dYF+n546FxlocxwCtvPsdV2v7N8BVkZ5UROYClwO7VfXsWtuHA3/CE5T+oqqzgDHA06r6gog8ASTl1fn555/Ttm1bjj3W8+3cd4OlnJwcrr/+el555RVuuukm+vbty5QpU9izZw8ZGRk89dRTdOnShaFDh36XnsRFSTkw76S1ApG3WA7shnm3Qs+roO8PIz9/CkjGVCdOOXnPTiciJNzvu3608/cAskPZL9QHcAGebq7NtbZlAB8BnfEErP8AZwHTgTzvPo+HUn6fPn20vvfee6/Btljav3+/9urVS7t166b//d//ratWrVJV1ezsbL3nnntq9uvfv78uWbJEVVUrKyv14MGDNa+tXLlSL7vssthWPIgPd+/XD3fvj3c1wndgj+qeD8J+vPfWKtW7T1ede2m830HMARE/kpXb73nhwoU1r2dnZ+vChQtjct4w3u86DfEzPNQxjmNFpBjIoVYrRVUvCvH4OlT1NRHJqbe5P/Chqn4MICKLgVHATiALz8LDgF1rIlIIFAJ07Ngx6PnvfOFd3tv1TSRVD+is045nxogeAV9v1aoV69ev5/XXX2flypVce+21zJo1C4Brr70WgP379/PZZ58xevRoAJo3b+5qHaMhLuMjTkXaYtlzBE7pCV+8A/Mui+zc1lpJW/n5+TUz2BJphXskQg0cTwEPA38BIksf27gOwKe1nu8EBgD3AQ+IyGXAC4EOVtVioBg8s6qiVEdHMjIyGDx4MIMHD6Znz5489thjgOeWshD/Aa9wxXXhYrz0jLiH1hNwICkDR/1rM1nSfjgRr/fs73Mg0X7foQaOI6r6UFRrAv7mr6mqHgRc/Z8WrGUQLVu3bqVJkyZ06+aZQrpx40ays7N55513avY5/vjjycrK4tlnn+WKK67g22+/pbq6usHd4RJFUo6PONX3h5F/8EfaSklAiZ72Ixri+Z4T7fcdauB4QURuBJYC3/o2qupXLtZlJ1B7mlEWsMvF8uPqwIEDTJ06lX379tG0aVO6du1KcXFxnam1AAsWLODHP/4xv/nNb2jWrBlPPfUUnTt3ZtCgQbz//vscOHCArKws/vrXv3LxxRfH6d24Iym7uZxKkW6udFw8GM/3nGi/71AXAG73s1lVtXPEJ/aMcSxT76wqEWkKfAAMBT4D1gLjVfXdcMtOpwWAySpuaegdcHwNrZsH7zwd2bFfvOMZX/nhi5GfPwoSZV1BLKXqOg7XFwCqaidnVapLRBYBg4G2IrITmKGqfxWRm/CkNskA5kYSNExysG6uMKVQN5dJfqGmVZ/gb7uq/l8kJ1XVcQG2/w34WyRlmvQSt/xc8ZQi3Vwm+YU6xtGv1s/N8XQnvQ1EFDiMccLJbC5I0hldaTqbyySmULuqptZ+LiKtAffzXxgTAjfycyUd6+YyCSRorqogKoDYpiY1xhiTEEId43gBz5J28Axcfw94MlqVMibaIhkj2bP/W2Y+soZReR0YPyB4doKE42R8BGyMxNQRaovjXmC293E3cIGq3h61WqWoL774grFjx9KlSxfOOusshgwZQmZmJnl5eZx00kl06tSJvLw8hg0bVnPMnDlzaN68OeXl5TXbVq1aRevWrcnLyyM3N5dhw4axe/fuiOqUk5PD3r17Hb83t61atYrLL78cgPnz5yMirFixoub1pUuXIiI8/XTo01t9ZQbLCPyvN19ncn7g8YT3Pv+G5zZ+FvI5E0LPqzxTeSP1xTuRTyM2KSnUMY7VInIy3w2S293/wqSqjB49mokTJ7J48WLAs3p8//79DBo0iEmTJnH55Zdz1VV1P7QWLVpEv379WLp0KZMmTarZPmjQoJrFg9OnT+fBBx90dJ/iSBw5coSmTUOdX+FMz549WbRoEUOHDgVg8eLF9OrVK6Kygo2RfHpCCzKPaep3jcjhvcdy1qnHR3TOuHIyPgI2RmIaCLWr6hrg98AqPKlB7heR21TVvoaEaOXKlTRr1owpU6bUbMvLywt6zEcffcSBAwf4/e9/z913310ncPioKvv376dr164AvPXWW9x8881UVlbSokUL5s2bR/fu3amuruYXv/gFL7/8MiLC5MmTmTr1uzkPlZWVjB49miuvvJLJkydz1113UVJSwumnn07btm3p06cPt956K4MHD+a8887jzTffZOTIkVx11VVcf/317Nmzh3bt2jFv3jw6duzYIBC2atWKAwcOsGrVKmbOnEnbtm3ZvHkzffr0YeHChYgIL730EjfffDNt27ald++69wgbNGgQr7/+OlVVVXz77bd8+OGHdX5/K1as4NZbb+XIkSP069ePhx56iGOPPTZgmQcPHmTq1Km88847HDlyhJkzZzJq1KiQ/pbvff4N1z6yJqR960vKbi5j6gn16+KvgH6quhtARNoBywELHCHyfUiGY9GiRYwbN45BgwaxdetWdu/eTfv27QF4/fXXycvLo6ysjJYtW3L33XcDcOaZZ/Laa6/RtGlTli9fzi9/+UueeeYZiouL2b59Oxs2bKBp06Z89dV32WIOHDjA2LFjmTBhAhMmTGDdunU888wzbNiwgSNHjtC7d+86dd+3bx+rV68GYMSIEUyYMIGJEycyd+5cfvrTn/Lss88GfV8bNmzg3Xff5bTTTuP888/nzTffpG/fvkyePJl//OMfdO3atSZjsI+IMGzYMF5++WXKy8sZOXIk27d7EhocOnSISZMmsWLFCs444wwmTJjAQw89xJQpUwKWWVRUxEUXXcTcuXPZt28f/fv3r9NFGMiovA6N7hPIe597MjJb4DDJLtTA0cQXNLzKiHxGVvz9/fbv5ra75ZSecMksV4tcvHgxS5cupUmTJowZM4annnqKn/zkJ0Ddrqp77rmHadOm8fDDD1NeXs7EiRPZtm0bIkJVlSetx/Lly5kyZUpN19JJJ51Uc55Ro0Yxbdq0mrw3b7zxBqNGjaJFC89d9kaMGFGnXrU/gNesWcOSJUsAuO6665g2bVqj76t///5kZWUBnlZXaWkprVq1olOnTjVJIAsKCiguLq5z3NixY7nvvvsoLy9n9uzZNcFy69atdOrUiTPOOAOAiRMn8uCDDzJ48OCAZb7yyis8//zz3HvvvYAn+HzyySeN1n38gI4Rf/BH2koxJtGEGjheEpGXgUXe59diK7zD0qNHj7AGcjdt2sS2bdv4/ve/D8Dhw4fp3LlzTeCobeTIkVx55ZUA3HHHHQwZMoSlS5dSWlrK4MGDAU+Xli/PTX3nn38+f//73xk/fjwi0mgeHF8aeH9852jatClHjx6tOffhw4dr9vHdBRE8qeaPHDlS59hA+vfvz+bNm2nRokVNkPCV31h96lNVnnnmGbp3715n+5dffhm0DsaYRgKHiHQFTlbV20RkDDAQzxjHGpL0Fq6A6y2DUFx00UX88pe/5NFHH2Xy5MkArF27loqKCi688MIG+y9atIiZM2cyffr0mm2dOnXym5P/jTfeoEuXLgCUl5fToYOnO2X+/Pk1+/zgBz/g4YcfZvDgwTVdVb5Wx29/+1vuuusubrzxRh566CEGDhzIj3/8Y6ZPn86RI0d48cUXa+pc33nnncfixYu57rrrKCkpYeDAgYBnttb69eu55ppreO6552paPoGceeaZbN++nY8++oguXbqwaNEiv/v97ne/a3CDqzPPPJPS0lI+/PBDunbtyoIFC7jwwguDlnnxxRdz//33c//99yMibNiwgXPOOSdoHd1g4yMmFTTW3fRHYD+Aqi5R1Z+r6s/wtDb+GO3KpRIRYenSpbz66qt06dKFHj16MHPmTE477TS/+y9evLjmToA+o0ePrpmR5Rvj6NWrFwsWLGD27NkATJs2jenTp3P++edTXf3dPbd+9KMf0bFjR3Jzc+nVqxePP/54nbL/+Mc/cujQIaZNm0a/fv0YOXIkvXr1YsyYMfTt25fWrVv7red9993HvHnzyM3NZcGCBfzpT38CYPLkyaxevZr+/fvz73//O2grBTx3OywuLuayyy5j4MCBZGdn+93vkksuYciQIQ2OnTdvHldffTU9e/akSZMmTJkyJWiZd9xxB1VVVeTm5nL22Wdzxx13BK2fG0bldYh4VlZSTgM2KStoWnUR2exLe+7ntXdU1cHk8OixtOrOHThwgFatWlFRUcEFF1xAcXFxg5lO6Sae15CvlfLEj8+N/cl903H9pHS3tOqpc24306oHu+l1i9CrZJJNYWEh7733HocOHWLixIlpHzTSXmMrz4O9ZqvOU05jgWOtiExW1UdrbxSRG4D10auWibf6XVkm/uI2PmKZeU09jQWOm4GlIpLPd4GiL3AMMDrgUcYYV8V1/UiwlefXe2etBbozoa06T0lBA4eqfgmcJyJDAN9Yx4uq+o+o18wYU8PWj5hEEmquqpXAyijXxRhjTBKITYY6Y0xcORkfAVtDYupK3rQhSah+WvVLL72UDz74IOTjW7VqmLHVmMY4WT8CtobENGQtjhgJlFb9yy+/rEmfUV1dTUaG//tEGBMpJ+MjYGMkpiFrccRIoLTq1dXVDBkyhPHjx9Ozp2c95RVXXEGfPn3o0aNHg0R/t9xyC71792bo0KHs2bMnpu/BGGMgCVocInIFcBnQHnhQVV+Jc5UiEiyt+ltvvcXmzZvp1KkTAHPnzuWkk06isrKSfv36ceWVV9KmTRsOHjxI7969mT17Nr/97W+58847eeCBB2L5NkyaamyMJNBrvykrp22rYzk5WhUzcRHVwCEic4HLgd21U5eIyHDgT3juX/4XVQ2YdVBVnwWeFZET8dzC1nHguOete3j/q/edFlPHmSedyS/6/yKiY/v3718TNMCT/2np0qUAfPrpp2zbto02bdrQpEmTmpTmBQUFjBkzxnnFjWmEkzUkFYer2XvgWwscKSbaLY75wAPA//k2iEgG8CDwfWAnntXpz+MJIr+rd/z1te4D8mvvcUkpWFr12gkAV61axfLly1mzZg2ZmZkMHjyYQ4cO+T2usTTkxrgh2BjJk96e10A5tN6928bsUlFUA4eqviYiOfU29wc+VNWPAURkMTBKVX+Hp3VSh3g+HWcBf1fVt92oV6QtAycCpVX33UnPp7y8nBNPPJHMzEzef/99/vWvf9W8dvToUZ5++mnGjh3L448/XpPC3JhEllP1ceQryC3PVUKKxxhHB+DTWs93AgOC7D8VGAa0FpGuqvqwv51EpBAoBOjYMfHmm/vSqt98883MmjWL5s2bk5OTwxVXXFFnv+HDh/Pwww+Tm5tL9+7d+a//+q+a11q2bMm7775Lnz59aN26NU888USs34YxdTSWpfXNFp4U+D0iKTxB81zFMxNwomQhDppW3ZUTeFocy3xjHCJyNXCxqv7I+/w6oL+qTnXrnJZW3USDXUPhu/aRNbz3+TcRrSP5TdltnoH1n66IQs1MfW6mVY+GncDptZ5nAbviUA9jTJQ5HVhv9fUW6+ZKQPEIHGuBbiLSCfgMGAuMj0M9jDFR5mTxYfGcYWRWrkypbq5UEe3puIuAwUBbEdkJzFDVv4rITcDLeGZSzVXVd6NZD2NM8lmReSkrMi/liR9GcNfDeZc1fvOpxliLJaBoz6oaF2D73/Dct9wYYwKKNDnj0IrejGp5KPL1I9ZiCSrhV44bY9KTk/GR+8sHsuLUCFsrYDegaoQFDmNMQrKbVyUuS3IYhIi4ujo7IyODvLw8evToQa9evfjDH/7A0aNHgx5TWlpq9/82xiQUa3HEUIsWLdi4cSMAu3fvZvz48ZSXl3PnnXcGPMYXOMaPt4lnxoTDyc2rLDljcBY4AigpKan5OScnh6KiIvLz810rv3379hQXF9OvXz9mzpzJjh07uO666zh48CAADzzwAOeddx633347W7ZsIS8vj4kTJzJ69Gi/+xljvuNkfARsDUljLHD4UVJSQmFhYc3zHTt21Dx3M3h07tyZo0ePsnv3btq3b8+rr75K8+bN2bZtG+PGjWPdunXMmjWLe++9l2XLlgFQUVHhdz9jzHec3rzK1pAEZ4GjlmDjGRUVFRQUFFBQUOBqvhhfWVVVVdx0001s3LiRjIyMgLeUDXU/Y0zkHK8hSXEWOOLo448/JiMjg/bt23PnnXdy8skn85///IejR4/SvHlzv8fMmTMnpP2MMc5EOkaSDuMjFjhq8X37z8nJYceOHQ1ez87OprS01JVz7dmzhylTpnDTTTchIpSXl5OVlUWTJk147LHHqK6uBuC4445j//79NccF2s8Y4x7LsRWcBQ4/ioqKKCwspKKiomZbZmYmRUVFjsqtrKwkLy+PqqoqmjZtynXXXcfPf/5zAG688UauvPJKnnrqKYYMGVJzc6fc3FyaNm1Kr169mDRpUsD9jDHusRxbwUU9rXo8uJFWvaSkhIKCAsDT0nB7VpVJPpZW3YTC170V6K6IQflybJ3SM7KTn9ITLgl4J+6gEj2telLIz8+vCRxudU8ZY9JD3HJsxYgFjiBSsTVmjImuuObYihELHMYY46J0yLGVVrmqrAVhImXXjjHfSZsWR/PmzSkrK6NNmzauJi40qU9VKSsrszUzJiac5Ng667TjmTEiovlcYUmbwJGVlcXOnTvZs2dPvKtiklDz5s3JysqKdzVMinOaYytW0mY6rjHGmMDCmY6bVmMcxhhjnLPAYYwxJiwWOIwxxoQlJQfHS0tL6ds3pK46Y4wxHr1D3TElA0dOTo7d3MgYY8IgIm+Hum/CdFWJyFwR2S0im2ttO0lEXhWRbd5/T4xnHY0xxiRQ4ADmA8PrbbsdWKGq3YAV3ufGGGPiKGECh6q+BnxVb/Mo4DHvz48BV8S0UsYYYxpImMARwMmq+jmA99/28aqIiFiqEmOMIfEDR8hEpFBE1onIOksrYowx0ZPogeNLETkVwPvv7kA7qmqxqvZV1b7t2rWLWQWNMSbdJHrgeB6Y6P15IvBcPCpRUlJS83NOTk7XUBNcAAASuUlEQVSd58YYk24SJnCIyCJgDdBdRHaKyA3ALOD7IrIN+L73eUyVlJRQWFhY83zHjh0UFhZa8DDGpC3LjluPkwHwVPxdGmPSg2XHNcYYEzUWOOpR1TqP7Oxsv/tlZ2c32NcYY9KBBY5GFBUVkZmZWWdbZmYmRUVFcaqRMcbElwWORuTn51NcXFzzPDs7m+LiYvLz8+NYK2OMiR8bHA+Rb9A8FX9fxhhjg+MJyFKWGGNSRUrejyMarKVhjDEe1uIwxhgTFgscMWApS4wxqcQCR5RZyhJjTKqxWVUusnQlxphkZbOqjDHGRI2rgUNEWopIhptlJpP6KUjCSVlijDHJwlHgEJEmIjJeRF4Ukd3A+8DnIvKuiPxeRLq5U83kZSlLjDGpxmmLYyXQBZgOnKKqp6tqe2AQ8C9glogUODxHUrOUJcaYVONocFxEmqlqldN93BavwfFgLGWJMSaRxWxw3BcQRGRGY/ukOxvLMMakCrdSjswQkUzgJOBtYLGqfu1S2cYYYxKIW7OqFDgEvAycDvxTRHq5VLYxxpgE4laL431V9XVXPS0i84GHgYvcKFxESoH9QDVwJNR+uFRh4yPGmETiVuDYKyJ9VHU9gKp+ICLtXCrbZ4iq7nW5TGOMMWFyq6vqp8BCEVkoIr8QkRJgu0tlpzVLkGiMSTROFwAKgKr+B8gDFnlfWgmMq72PQwq8IiLrRaSw0b1ThCVINMYkIqfrOFYBzwDPqeontbYfAwwEJgIrVXW+o0qKnKaqu0SkPfAqMFVVX6u3TyFQCNCxY8c+O3bscHLKmLMEicaYeIplksPheAasF4nILhF5T0Q+BrbhaXHMcRo0AFR1l/ff3cBSoL+ffYpVta+q9m3Xzu3hFWOMMT5OFwAeUtU/q+r5QDYwFOitqtmqOllVNzqtoDdx4nG+n4EfAJudlptoLEGiMSZZuJYdV1WrVPVzVd3nVpleJwNviMh/gLeAF1X1JZfPkZAsQaIxJhG5NR03alT1YyAtFxP6EiEWFHjyRGZnZ1NUVGQJEo0xcWV3AEwCtgDQGBNtdgdAY4wxUeNK4BCRq2sNYP9aRJaISG83yjaWWdcYk1jcanHcoar7RWQgcDHwGPCQS2UbB0TE0RoRY4ypz63AUe399zLgIVV9DjjGpbKNMcYkELcCx2ci8ghwDfA3ETnWxbJNhCzPlTEmGtz6cL8Gz704hnvXcZwE3OZS2SYCbuS5sm4uY4w/Nh03BTj9cA90Ddg0YGPSR8ym44rIfhH5xs9jv4h846RsE1/WzWWMCcRprqrjVPV4P4/jVPV4typpgnOS58pfa8LSuRtjgnFtAFtEThSR/iJyge/hVtkmfOHkufKNZfgeBQUFVFRU1NmnoqKCgoKCBvsaY9KPK7mqRORHwP8AWcBG4L+ANbh0z3ETPstzZYyJFrdaHP8D9AN2qOoQ4Bxgj0tlmwjVDhKlpaUBg4ab3VzGmNTnVuA4pKqHAETkWFV9H+juUtkmxiyduzEmGLfSqu8UkROAZ4FXReRrYJdLZRsHImkVWDeXMSYY19dxiMiFQGvg76pa5WrhIUq3dRzRYus4jEkf4azjcGtw/Dd+NucBv3WjfGOMMYnDra6qg7V+bg5cDmxxqWwTJ9bSMMb440rgUNXZtZ+LyL3A826UbYwxJrFEK4NtJtDZrcJEZLiIbBWRD0XkdrfKNcYYEz63xjjeAXz9GhlAO+Aul8rOAB4Evg/sBNaKyPOq+p4b5RtjjAmPW2Mcl9f6+Qjwpaoecans/sCHqvoxgIgsBkYBFjiMMSYOHAUOEfl5kNdQ1T84Kd+rA/Bprec7gQHBDigrK2P+/Pl1tvXo0YN+/fpRVVXlN1lfXl4eeXl5VFRU8OSTTzZ4vW/fvpx99tmUl5ezdOnSBq+fe+65dO/enb1797Js2bIGr19wwQV07tyZL774gpdeeqnB60OHDuX000/n008/ZcWKFQ1eHz58OKeccgoff/wxr732WoPXL7/8ctq2bcvWrVtZs2ZNg9dHjx5N69at2bx5M/6mKl9zzTVkZmayceNGNm7c2OD1/Px8mjVrxtq1a3n33XcbvD5p0iQA/vnPf/LBBx8AsHr1agCGDRtWswZk9erVbN++vc6xmZmZXHPNNQAsX76cnTt31nn9+OOPZ8yYMQC89NJLfPHFF3Veb9OmDSNGjADghRdeoKysrM7rp5xyCsOHDwdgyZIlfPNN3cTNWVlZDBs2DIAnn3yyQZ6uTp06ceGFFwKeBJBVVXVnmZ9xxhmcd955AA2uO7BrLx7Xnk+zZs3s2gvx2guH0xbHcd5/u+NJOeIbEB8BNLzCIuMvk16D6T4iUggUAnTo0MGlUxtjjKnPlQWAIvIKcKWq7vc+Pw54SlWHu1D2ucBMVb3Y+3w6gKr+LtAxtgAwMdgCQmPqSuT/EzG7kVMtHYHDtZ4fBnJcKnst0E1EOonIMcBYbKpvwrMbQRnjvkS5nYFbg+MLgLdEZCmebqTRwGNuFKyqR0TkJjz3NM8A5qpqw45OkzAC3QgKsHxXxqQA13JViUgfYKD36WuqusGVgiNgXVWx4+TbTyI2142JlpKSEseJQ6PZ1RWPripUdb2q/sn7iFvQMMakh0TptgmFG7djTqTuX0ctDhF5Q1UHish+6s50EkDjdd9xa3HEV05ODjt27GiwPTs7m9LS0thXyKSkBB9ojvhYf+/HF3hqT5vNzMykuLjYte7fmLU4VHWg99/jVPX4Wo/j4hU0TPzZjaBMtCXSt2+3+VpStR8FBQUN1lpUVFRQUFBQZ79YcaWrSkSu9k7BRUR+LSJLROQcN8o2ySc/P5/i4uKa59nZ2a5+MzLpzY1uHydC+ZBO9dsxu7WOY5Oq5orIQOB3wL3AL1U16ArvaLGuqsSQyF0JJjkk4uSLSK5rp11Nsej+jcfgeLX338uAh1T1OeAYl8o2SSqZvkGZ2EumwW2nnLbCE637163A8ZmIPAJcA/xNRI51sWxjTJqq340TTrdPNDgZW6kdJEpLS8Pquk207l+3PtyvwbNAb7iq7gNOAm5zqWxjjKkRr2/fboytOAlqTgKP21wJHKpaoapLVHWb9/nnqvqKG2UbY1KP02/usfj2HenMpnTofnNrcFyAfKCzqv5WRDoCp6jqW44Lj4ANjhuTuNxakxDtyReJODAfTeEMjrsVOB4CjgIXqer3RORE4BVV7ee48AhY4EhvNpsrsUTrAzjWf+dUX9gaj1lVA1T1J8AhAFX9GptVZYyJoljP2ku0mU3x5FbgqPLeG1wBRKQdnhaIMRFJxr7iZKxzLCTiYrhI/laJNrMpntwKHPcBS4H2IlIEvAHc7VLZxoQslVNRpIpk/uaeSDOb4smtWVUlwDQ8q8Z3AVdjXVUmxuKZisICVuji/c3d6d/KFrbif4FNqA/geGA68ADwAzxZcacCO4DnnJTt5NGnTx81yWvhwoWKp9tTs7OzdeHChQ328b0e6cPt+mZmZtYpPzMz02+9zXei8bdojP2tAgPWaYifsU7Tqj8HfA2sAYYCJ+JpafyPqm6MuGCHbFZV8gp1qqbTsQSH131czptqEnk6LaTf3ypm03FF5B1V7en9OQPYC3RU1f0RF+oCCxzJwe0P4FhNl7TAkRwscIQnltNxq3w/qGo1sD3eQcOkr1gNutZvtifCLCHTkL8uFvtbucNp4OglIt94H/uBXN/PIvKNGxU0qcvJf2x/4jXomsyzhNKN/a1cEupgSDwewEzgM2Cj93FpKMfZ4HjycmPwEgeDrpEeG8qAfqJy8vtKRsn8t4omYjU4Hm0iMhM4oKr3hnOcjXEkt5KSEgoKCgBPq6GoqChmUzWdDNg6HeyNV6qUeL7neEnWekdTPFKOGOMaW2RlTGJLhsBxk4hsEpG53uSJxkRFOi4Mc/Kek3nRYzL+rRJKqH1a0XoAy4HNfh6jgJOBDDwBrgiYG6ScQmAdsK5jx46u9PmZ9BHvhWHx6Hd38p7j/fsy7iNVxjhqE5EcYJmqnt3YvjbGYRqTSGsx3Lo/RTDxTL6YLJ8x6S7m9+OIFhE5VVU/9/78Mzzp28c2dpwFDtOYeAWOeN2bwgKHaUwqDY7/PxF5R0Q2AUOAn8W7QiY11G96p/rCMH/dDeG853T7fZngEjpwqOp1qtpTVXNVdaSv9WGM2+K16jycD+FAIh2kdvKebSFdmgt1MCSZHrYA0EQiXgvD4jlI7eQ920K61EIqDo6Hw8Y4TKTitTAslEWPiXjvbltIlzpSZnA8UhY4TDKK5gC3BQ7TmFQaHDfGePnrMnA6PlK7XCd1MunFAocxCSKSD2EbpDbxYIHDmCQW7/t3m/RkYxzGpAAbazBO2RiHMcaYqGka7woYY5yzloaJJWtxGGOMCYsFDmOMMWGxwGGMMSYsFjiMMcaEJSWn44rIHmBHI7u1Bsod7hPp622BvY2cOxGE8jtKhHNEWkY4x0X7egn2ml0v7pYfSTluXyuh7BfJ9eLkWslW1XYh7RlqNsRUewDFTveJ9HXCyEKZ6L+jRDhHpGWEc1y0r5dGXrPrxcXyIynH7WsllP0iuV5ida2kc1fVCy7s4/T1RBeL+rtxjkjLCOe4aF8vyX6tQPTfg1vlR1KO29dKKPsl7PWSkl1ViU5E1mmIKzSNsevFhCpW10o6tzjiqbjxXYypYdeLCVVMrhVrcRhjjAmLtTiMMcaExQKHMcaYsFjgMMYYExYLHAlARFqKyGMi8qiI2B14TEAi0llE/ioiT8e7LibxicgV3s+V50TkB26Va4EjSkRkrojsFpHN9bYPF5GtIvKhiNzu3TwGeFpVJwMjY15ZE1fhXCuq+rGq3hCfmppEEOb18qz3c2UScK1bdbDAET3zgeG1N4hIBvAgcAlwFjBORM4CsoBPvbtVx7COJjHMJ/RrxZj5hH+9/Nr7uisscESJqr4GfFVvc3/gQ++3xsPAYmAUsBNP8AD7m6SdMK8Vk+bCuV7E4x7g76r6tlt1sA+p2OrAdy0L8ASMDsAS4EoReYjUSD1hnPN7rYhIGxF5GDhHRKbHp2omAQX6bJkKDAOuEpEpbp3Mbh0bW+Jnm6rqQeCHsa6MSWiBrpUywLUPAJMyAl0v9wH3uX0ya3HE1k7g9FrPs4BdcaqLSWx2rZhwxPR6scARW2uBbiLSSUSOAcYCz8e5TiYx2bViwhHT68UCR5SIyCJgDdBdRHaKyA2qegS4CXgZ2AI8qarvxrOeJv7sWjHhSITrxZIcGmOMCYu1OIwxxoTFAocxxpiwWOAwxhgTFgscxhhjwmKBwxhjTFgscBhjjAmLBQ6TVkSkWkQ21nrc3vhRsSEiT3vvt/Fvb90+EZE9teqaE+C4/xWRu+pt6ysim7w/rxCR1tF/ByZd2DoOk1ZE5ICqtnK5zKbeBVhOyugB/K+qjq61bRLQV1VvCuHYpap6Rq1t9wJlqvo7EbkBaKuq9zipozE+1uIwBhCRUhG5U0TeFpF3RORM7/aW3hvnrBWRDSIyyrt9kog8JSIvAK+ISBMR+bOIvCsiy0TkbyJylYgMFZGltc7zfRFZ4qcK+cBzIdTzEhFZ463nEyLS0rtC+JCI9PHuI8DVeFJr4y13vJPfjzG1WeAw6aZFva6q2ndF26uqvYGHgFu9234F/ENV+wFDgN+LSEvva+cCE1X1Ijx3ccwBegI/8r4G8A/geyLSzvv8h8A8P/U6H1gfrOIi0h64HRjqrecm4H+8Ly/Ck5/IV9YuVd0OoKp7geNE5IRg5RsTKkurbtJNparmBXjN1xJYjycQAPwAGCkivkDSHOjo/flVVfXdUGcg8JSqHgW+EJGV4MlrLSILgAIRmYcnoEzwc+5TgT2N1P08PHd3+6enUcExwBve1xYBq0VkGp4AsqjesXu859jXyDmMaZQFDmO+863332q++78hwJWqurX2jiIyADhYe1OQcufhuUHXITzBxd94SCWeoBSMAC+p6nX1X1DVUhHZBQwCRgN96u3S3HsOYxyzripjgnsZmOodN0BEzgmw3xt47uLYREROBgb7XlDVXXjujfBrPPeL9mcL0LWRuvwTuFBEOnvr0lJEutV6fRGem/ZsUdUvfBtFpAnQlrp3iDMmYhY4TLqpP8Yxq5H97wKaAZtEZLP3uT/P4LmZzmbgEeDfQHmt10uAT1X1vQDHv0itYOOPqn4J3AA8ISL/wRNIzqi1y5PA2Xw3KO7TH3hDVauDlW9MqGw6rjEuEZFWqnpARNoAbwHn+775i8gDwAZV/WuAY1sAK73HuPoBLyIP4rk/w2o3yzXpy8Y4jHHPMu/MpWOAu2oFjfV4xkNuCXSgqlaKyAygA/CJy/XaYEHDuMlaHMYYY8JiYxzGGGPCYoHDGGNMWCxwGGOMCYsFDmOMMWGxwGGMMSYsFjiMMcaE5f8DMhLDSE0glyoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1094e80b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import sys\n",
    "import os\n",
    "sys.path.append(os.environ['CTOOLS']+'/share/examples/python/')\n",
    "\n",
    "from show_residuals import plot_residuals\n",
    "plot_residuals('resspec_1.fits','',0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the residuals it is clear that the model does not reproduce the data well. There is an excess at low energies and the model overshoots the data at high energies. You may try to change the spectral model for the Crab nebula.Try for example an exponentially cutoff power-law rather than a simple power law. Let's use as starting values for the Prefactor and Index the current best-fit values and tentatively set the cutoff energy at 1 TeV."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "crab    = like.obs().models()['Crab']\n",
    "expplaw = gammalib.GModelSpectralExpPlaw()\n",
    "expplaw['Prefactor'].value(crab['Prefactor'].value())\n",
    "expplaw['Index'].value(crab['Index'].value())\n",
    "expplaw['PivotEnergy'].value(crab['PivotEnergy'].value())\n",
    "expplaw['CutoffEnergy'].value(1.e6)\n",
    "crab.spectral(expplaw)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now fit the modified model to the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "like.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and use TS to quantify the model improvement (in this case there is only one additional degree of freedom, the cutoff energy)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "319.8230714731617\n"
     ]
    }
   ],
   "source": [
    "like3 = like.opt().value()\n",
    "ts    = -2.0 * (like3 - like2)\n",
    "print(ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model improvement was again significant. Let's look at the fit results for the Crab nebula."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== GModelSky ===\n",
      " Name ......................: Crab\n",
      " Instruments ...............: all\n",
      " Observation identifiers ...: all\n",
      " Model type ................: PointSource\n",
      " Model components ..........: \"PointSource\" * \"ExponentialCutoffPowerLaw\" * \"Constant\"\n",
      " Number of parameters ......: 7\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 ..: 4\n",
      "  Prefactor ................: 3.5947627991423e-16 +/- 7.19602031907029e-18 [0,infty[ ph/cm2/s/MeV (free,scale=1,gradient)\n",
      "  Index ....................: -1.82160233651441 +/- 0.0226092875383086 [-10,10]  (free,scale=1,gradient)\n",
      "  CutoffEnergy .............: 3206969.84019295 +/- 117128.52342759 [0.1,infty[ MeV (free,scale=1,gradient)\n",
      "  PivotEnergy ..............: 300000 MeV (fixed,scale=1,gradient)\n",
      " Number of temporal par's ..: 1\n",
      "  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)\n",
      " Number of scale par's .....: 0\n"
     ]
    }
   ],
   "source": [
    "print(crab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The best-fit cutoff energy is 3.2 TeV, which seems consistent with the behavior in the energy-bands residual maps above. Note that this has changed significantly the Index value, that now is 1.82. Check the residual spectrum again to make sure that the residuals are reduced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEOCAYAAACetPCkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xt4VNW5+PHvS7iGWxXEC0gClEtFQrgrgqAgRVQQKwISQWiN/Cw9bU+tRVsraHPUUzxU1KPGFvBABEXBC6IoKCKVKiARAY0oBEVEMNQIhiiE9/fHzMRJMpPMzN5zS97P88xD9p69114TdvY7a6+13yWqijHGGBOqevGugDHGmORigcMYY0xYLHAYY4wJiwUOY4wxYbHAYYwxJiwWOIwxxoTFAocxxpiwWOAwxhgTloQPHCLyExF5RESeFpH/F+/6GGNMXReXwCEi80XkoIhsr7R+pIgUiMjHIjITQFU/UNXpwDVA33jU1xhjzA8kHilHRORC4Cjwf6p6rnddCvARcAmwD9gETFTVnSIyGpgJPKiqT9RUfuvWrTU9PT1a1TfGmFppy5YtX6nqaTVtVz8WlalMVdeLSHql1f2Bj1V1N4CILAXGADtV9XngeRF5EagxcKSnp7N582Z3K22MMbWciOwNZbtE6uNoC3zmt7wPaCsiQ0Vknog8CqwKtrOIZIvIZhHZfOjQIdcrl5eXR3p6OvXq1SM9PZ28vLyY7GuMMQlHVePyAtKB7X7L44C/+y1fBzwQSdl9+vRRNy1evFhTU1MVKH+lpqbq4sWLo7qvb/+0tDQVEU1LSwt5P2OMCRewWUO5foeyUTReAQLH+cBqv+VbgVsjKTvSwDFr1qwKF/hYvyqLd9CxoGVM3RJq4IhLH0cQm4DOItIB+ByYAFwb3yrFlojUuE1JSQlZWVlkZWVVWD9r1qwKy9u2beOFF17g+PHjAOzdu5epU6eyfPlyMjIyKmzb7Kxm/O6G31VYl5eXR3Z2NiUlJeX7Z2dnAzBp0qSwPpcxpnaJ16iqJcBQoDXwJXCHqv5DREYBfwNSgPmqmhNJ+X379tVIOsfvfedePjz8YZX1y6Ys49uD31ZZ37RNU8Y9Pg6ANgVtSP06tco2c+fOpbi4uMr6li1b8tvf/rbCusoX/2QwJ3dOlaADnsDzxz/+kU8//ZT27duTk5NjAceYBCciW1S1xsce4hI4oi3SwBFM5W/fAKmpqeTm5tZ4MXSyb3p6Onv3Vh3kkJaWRmFhYYV1obRWoiFQsKvc2gFo0KABV1xxRYXWTqCWDljQMSZeLHC4PBzXycUs0n1jFXQi3f++x+7j5uybaywrGCdBB4LfYrOgY0xkLHDUkuc44hF0wtk/Xi0dqLlfB8ILOsbUdaEGjriNqormy+3huMkqXqOq0tLSAo4aS0tLq7BdoG1i9Zo1a1bC/L6MSRSEOKrKWhzGdcnQr+N2a8VpC8+YRGAtDhNXkX77dvLsSqgtnTm5cxLmWR1jEgmJ/gBgNF8WOJJbLIJOvALHnNw5rn5mY9wUauBIpAcAjQE8DxhGcnvHt08ogwk8fyM/CPUW2X2P3cfR/UerbBfq8zqB9rWHLU3SCSW6JNvLWhwmXG6kd6lpfzdS2gRrsRjjBkJscSRSdtz4e2kmLLgs8tfmBfH+BCZCkyZNIjc3l7S0NESEtLS0sDq2Q9m/2VnNHNczWIslGbMvW8bp5GWjqvy9NBMOvB/ZQfdu8PybNiiy/QF6XA19p0a+v0lKod4mmz17NgB33HFH+To3nteJxwOT8crGYKpnDwBGEDhmv7CDnfu/ieiYw0pWMSblLU5v3jii/R0HHgs6SSvUC+Hs2bMd5TPLy8vj2muvrbAciwtwPB8SrY3Xt2iywBHjwPH2nsMADOhwakT7Owo8FnSSXijf/P9w/x/479/8d1zq5+Q6YYHjB4meEscCR4wfAHzi7U95Lv/ziPd3EnjiGnScsIAVlmUfLWPV7oqTYIaSuRkg/e10R60Vt68TTnKpOc3DFi/JcIut1gQOEWkK/C/wPbBOVWvsBUvGJ8edBJ6dX3zDOWe24Mkbzw9/580L4P2nIzquI9Yn5IpwbnNBxf6ReF6Aa3sfh5NWVjyvyQkdOERkPnA5cFBVz/VbPxK4H898HH9X1XtE5Drga1V9QUSeVNXxNZWfjIHDifGPbiwPHpEYk9mWawe0d7lWNXAasOz2XLlQbn/4Aoe/cNKsdO7cuUL/iP+xfZOKpaWlxSRrtNN9nQj180YjcDj5XYcq0QPHhcBR4P98gUNEUoCPgEuAfXhmBJwIjAFeUtV8EXlCVWucFbCuBQ4nrRWnfTNxCTrgLPDUwaDzh/v/EHCisW3btrF27VqKi4tp2bIlw4YNqxI0fPxbK5Ac3/zd5PTzOmnhxXAgQ+IGDgARSQdW+gWO84FZqvpT7/Kt3k33Af9W1ZUislRVJ9RUdl0LHE7E7RZZPNXBoBOofyRUidY/EitutxpiMVWB0991qIEjkVKOtAU+81veBwwA5gEPishlwAvBdhaRbCAboH37qt+Ajx8/zr59+ygtLXWzzkmvVwvodeGPItr30JEmHC87yfp38iPaP7VhCk0bxeEUbHoenHde2Ls1btyYdkVv0uD9JyM77t4NnpeTW3QRBp5xXcYxrsu4mjcMYPbbVW9zmfD5gkO0bzfFQiK1OMYBP1XVX3iXrwP6q+qvwi07UItjz549NG/enFatWsV1eGBtUnT0O74+drzmDQP49rsTAI4Cx4+aNKBVs0YR7x8OVaWoqIgjR47QoUOHyApx2q9z4H04owdMfTHyMiIQqGMdknd0U6Ti+XljdexkbHHsA872W24H7Her8NLSUtLT0y1ouKhVs0YRX7idBB2A0uNlfO2tQyyICK1ateLQoUORF9J3qrPbVAsu8wSPBZdFtr/D22SVO9j79OnD/v37q3Su9+nTp8K2wTrWk01OTk7AW005OTm1+tiBJFLg2AR0FpEOwOfABMDVs82CRuJwEnQAPjl0lNLjZXxyqGruplBE0lqJ+/nT4+rI93Vwm6ykSRqpx06vst7Xib58+XKAoJ3ru3btiqDCiSeet5oS7TZXvEZVLQGGAq2BL4E7VPUfIjIK+Bue4bjzVTWicBroVtUHH3zAT37yE0f1dkpEyMrKYtGiRQCcOHGCM888kwEDBrBy5cqQy0lPT2fz5s20bt3a0TbJzEmLpfR4GY0bpNDptPCTDibCeRQRB7fJln3zIatanuK5TRbAwksXAnD9S9dXeS/97XSg6m2uZOb7AhGna2dUj53Qt6pUdWKQ9auAyIZ+JIGmTZuyfft2jh07RpMmTXj11Vdp27ZtvKuVlJy0WJy0Vg4d+Y6tb38anyHITji4TTZuwWWMU2Bk4OzPC1kIwIIA71vHeu1kadVj7NJLL+XFFz2dm0uWLGHixB9i6OHDh7nyyivJyMjgvPPOY9u2bQAUFRUxYsQIevXqxY033ljh28bixYvp378/mZmZ3HjjjZSVlcX2AyWhHzVpQOMGKRHte7zspKPUMsbUBonUx1EnTJgwgTvvvJPLL7+cbdu2MW3aNN58803A05zv1asXzz77LK+99hqTJ08mPz+f2bNnM2jQIP785z/z4osvkpubC3humzz55JP885//pEGDBtx0003k5eUxefLkeH7EhOektfL5nnrs/OIbxj+6MaL94/bApDEuqpOBI1D6BTeEch83IyODwsJClixZwqhRoyq8t2HDBp555hkALr74YoqKiiguLmb9+vXlHZCXXXYZp5xyCgBr165ly5Yt9OvXD4Bjx47Rpk0bNz+SqSS1YUrEqV12fuHJvJyUgSOU0VwB3+8PRP43V1tGZNU2dTJwxNvo0aO5+eabWbduHUVFReXrA3V4+TrDAo3oUVWmTJnC3XffHb3KmgqaNqof8ZPyvpxikbZWIE4tFgejuUoa7Sf1u7Mi3r+2jMiqbepk4Ij3CI9p06bRsmVLevTowbp168rXX3jhheTl5XH77bezbt06WrduTYsWLcrX/+lPf+Kll17i3//+NwDDhg1jzJgx/Pa3v6VNmzYcPnyYI0eOkJaWFqdPZqozJtPZQIi4tVhq6lif5v1SE+DBxL4LhrBKtgQdkVUd34gsk3jqZOCIt3bt2vHrX/+6yvpZs2YxdepUMjIySE1N5fHHHwc8gW7ixIn07t2bIUOGlKdUOeecc/jLX/7CiBEjOHnyJA0aNOChhx6ywJGgrh3Q3tFF32mLJR6tlXE0Y5w2Czoiqzo2IitxWeCIoaNHqw7/HDp0KEOHDgXg1FNP5bnnnquyTatWrXjllVfKl+fOnVv+8/jx4xk/vmqm+dqY8qGuc9JiSer+FZNwwg4cInIKcLaqbotCfYwxQThpsSRja8XHyWAW61yPjpACh4isA0Z7t88HDonIG6r6n1GsmzHGJXFtrUSYXytYqpNwWOd6dISUckREtqpqLxH5BZ7Wxh0isk1VA8/4EmeJmnLEJL+6eB45mWFyWMkqxqS8xenNG4e979TvPgJgQaMuYe8LMPvT/hHt51MXWytupxypLyJnAtcAf3RUM2NMUnHSWnmgeBBrzxzFk1MjGMK87FL49hBEmJaps+5ml3SMbGestVKdUAPHbGA1sEFVN4lIR8B+q8bUAfHqWylsmEpZ/fpMPS3Ch1qbHWDU0ULGEX4yS6etldou1MDxhf9tKVXdLSL/E6U6VSAiPwF+jSeT7lpVfTgWxzXGOOektVJyOIPUUyM/dsHxYjizK+MiGAqMt0PenngPLNTA8QDQO4R1FYjIfOBy4KBvpj/v+pHA/XjSp/9dVe8JVoaqfgBMF5F6wGMh1jdh5eTk8MQTT5CSkkK9evV49NFHGTBgQEj7jhw5kn/9618MGjQorDTsxsSLs9YK7Nx7HiXfR5bipazhHA5+811E+3Zu/DW7SiObUhlq/22uagOHiJwPDAROExH/EVQt8Fz0a7IQeBD4P78yU4CHgEvwzPq3SUSe95ZXOXfGNFU9KCKjgZnespLWxo0bWblyJe+++y6NGjXiq6++4vvvv6+wTVlZGSkpgX+1v//97ykpKeHRRx+NRXWNiSunT9qXfHeCryLc99phGRHPX+K7zVWbWys1tTgaAs282zX3W/8NUGMCG1Vd751b3F9/4GNV3Q0gIkuBMap6N57WSaByngeeF5EXgSdqOm6i+uKLL2jdujWNGnkys/omWUpPT2fatGm88sorzJgxg759+zJ9+nQOHTpESkoKy5Yto1OnTgwbNqxCihJjajOnT9oPWFCfku9ORPjsSjfGZD4c0fE73/tbx60VJ8+uxCKlUrWBQ1XfAN4QkYWqWnWm9Mi0BT7zW94HBL1XIyJDgauARlQzyZOIZAPZQHlKjkQzYsQI7rzzTrp06cLw4cMZP348Q4YMAaBx48Zs2LABgAEDBjBz5kzGjh1LaWkpJ0+ejGe1jUlKrZs28rQ4vq9py6qcPLvipLXyxMEujoJOrITax9FIRHKBdP99VPXiCI4ZaOLmoAPuVHUdsK6mQlU1F8gFz3Mc1W07+4Ud7Nz/TU1FhuWcs1pwxxXdq92mWbNmbNmyhTfffJPXX3+d8ePHc889nu4dX9qQI0eO8PnnnzN27FjAE1CMMeFr06IRbVo0YsHI8IcCO3vSPvLWSmLfoPpBqIFjGfAI8HfA6RRz+4Cz/ZbbAfsdlpk0UlJSyvNT9ejRozyRYdOmTYH4zGNsjKnI8oJVL9TAccLFYbCbgM4i0gH4HJhAjANtTS2DaCkoKKBevXp07twZgPz8fNLS0nj//ffLt2nRogXt2rXj2Wef5corr+S7776jrKyM1NTUuNTZmLqoruYFC1Woc46/ICI3iciZInKq71XTTiKyBNgIdBWRfSLyc1U9AczA80DhB8BTqroj4k+QRI4ePcqUKVM455xzyMjIYOfOncyaNavKdosWLWLevHlkZGQwcOBADhw4AMDgwYMZN24ca9eupV27dqxevTrGn8AYU5MxmW0dzRKZDHPah5qrak+A1aqqkT/PH0WWq8pEi51HyWPqy1MpOFxA11O7RrT/qI6jGNdlnMu1qp6TvGAQWl9rdVzNVaWqHSKuiTHGxMGojqMi3rfgcAFAzAOH02dXYiXUtOqTA61X1f8LtN4YY+JtXJdxEV/4p75czVS5UeT02ZVYCbVzvJ/fz42BYcC7+D0Rbowxpm4I9VbVr/yXRaQlsCgqNTLGGJPQQh1VVVkJ0NnNihhjjEkOofZxvMAPT3enAD8BnopWpYwxxiSuUFscc4D7vK//Ai5U1ZlRq1UtdeDAASZMmECnTp0455xzuOiii0hNTSUzM5NTTz2VDh06kJmZyfDhw8v3mTt3Lo0bN6a4uLh83bp162jZsiWZmZlkZGQwfPhwDh48GFGd0tPT+eqrSHOIRs+6deu4/HJPzsuFCxciIqxdu7b8/RUrViAiPP106DmB/Mt0so0xdV1IgcOb7PBDPBlyTyGitGF1m6oyduxYhg4dyieffMLOnTuZO3cuq1evJj8/n9GjR/PXv/6V/Px81qxZU77fkiVL6NevHytWrKhQ3uDBg8nPz2fbtm3069ePhx56KNYfiRMnTsTsWD169GDJkiXly0uXLqVnz54xO74x5gchBQ4RuQZ4BxiHZ97xt0WkxrTq5gevv/46DRo0YPr06eXrMjMzGTx4cNB9PvnkE44ePcpf/vKXChdNf6rKkSNHOOWUUwB45513GDhwIL169WLgwIEUFHjGo5eVlXHzzTfTo0cPMjIyeOCBByqUc+zYMUaOHMljj3nmyrrrrrvo1q0bl1xyCRMnTmTOnDkADB06lNtuu40hQ4Zw//33s3fvXoYNG0ZGRgbDhg3j008/BeD666+v0Bpo1swzfee6desYOnQoV199Nd26dWPSpEnl+blefvllunXrxqBBg1i+fHmF+g0ePJh33nmH48ePc/ToUT7++GMyMzPL31+7di29evWiR48eTJs2je+++67aMr/99lumTZtGv3796NWrF88991zQ/wdjTEWhDsf9I9BPVQ8CiMhpwBogstzBddD27dvp06dPWPssWbKEiRMnMnjwYAoKCjh48CBt2njmX37zzTfJzMykqKiIpk2b8l//9V8AdOvWjfXr11O/fn3WrFnDbbfdxjPPPENubi579uxh69at1K9fn8OHD5cf5+jRo0yYMIHJkyczefJkNm/ezDPPPMPWrVs5ceIEvXv3rlD3r7/+mjfeeAOAK664gsmTJzNlyhTmz5/Pf/zHf/Dss89W+7m2bt3Kjh07OOuss7jgggv45z//Sd++fbnhhht47bXX+PGPf1yeLdhHRBg+fDirV6+muLiY0aNHs2ePJ6FBaWkp119/PWvXrqVLly5MnjyZhx9+mOnTpwctMycnh4svvpj58+fz9ddf079//wq3CI0xwYUaOOr5goZXEZGPyIq/l2bCgfdr3i4cZ/SAS4POgBuRpUuXsmLFCurVq8dVV13FsmXL+OUvfwl4voH7po+99957ueWWW3jkkUcoLi5mypQp7Nq1CxHh+PHjAKxZs4bp06dTv77nv/zUU39INTZmzBhuueUWJk2aBMCGDRsYM2YMTZo0ATzBwZ//BXjjxo3l3+Svu+46brnllho/V//+/WnXrh3gaXUVFhbSrFkzOnToUJ4AMisri9zc3Ar7TZgwgXnz5lFcXMx9991XHiwLCgro0KEDXbp0AWDKlCk89NBDDB06NGiZr7zyCs8//3x5S6q0tLS8tWSMqV6ogeNlEVkN+O6XjKeaSZVMVd27dw+rI3fbtm3s2rWLSy65BIDvv/+ejh07lgcOf6NHj+ZnP/sZALfffjsXXXQRK1asoLCwkKFDhwKeW1oigaZCgQsuuICXXnqJa6+9FhGpMbW7LwV8IL5j1K9fv3wCKlWtMEWubwZE8KSZ9/WVBKufT//+/dm+fTtNmjQpDxK+8muqT2WqyjPPPEPXrhXzGH355ZfV1sEYU/Oc4z8GTlfV34vIVcAgPBMxbQTyYlC/6HC5ZRCKiy++mNtuu43HHnuMG264AYBNmzZRUlJSPgugvyVLljBr1ixuvfXW8nUdOnRg796qEzFu2LCBTp06AVBcXEzbtp58NwsXLizfZsSIETzyyCMMHTq0/FaVr9Vx5513ctddd3HTTTfx8MMPM2jQIG688UZuvfVWTpw4wYsvvlhe58oGDhzI0qVLue6668jLy2PQoEGAZ7TWli1buOaaa3juuefKWz7BdOvWjT179vDJJ5/QqVOnoH06d999d5XJrbp160ZhYSEff/wxP/7xj1m0aBFDhgyptsyf/vSnPPDAAzzwwAOICFu3bqVXr17V1tHULQWHCyJOPRKPBImxVNPtpr8BRwBUdbmq/qeq/hZPa+Nv0a4cgIjUE5EcEXlARKbE4pjRICKsWLGCV199lU6dOtG9e3dmzZrFWWedFXD7pUuXls8C6DN27FiWLl0K/NDH0bNnTxYtWsR9990HwC233MKtt97KBRdcQFnZD3Nu/eIXv6B9+/ZkZGTQs2dPnnii4tTtf/vb3ygtLeWWW26hX79+jB49mp49e3LVVVfRt29fWrZsGbCe8+bNY8GCBWRkZLBo0SLuv/9+AG644QbeeOMN+vfvz9tvv11tKwU8Mx3m5uZy2WWXMWjQINLS0gJud+mll3LRRRdV2XfBggWMGzeOHj16UK9ePaZPn15tmbfffjvHjx8nIyODc889l9tvv73a+pm6ZVTHURFn1S04XMCq3bX7hky1adVFZLuqnhvkvfdVtUe1hYvMBy4HDvqXIyIjgfvxPEz4d1UN2gQQkbHAGOAw8KKqrg22rY+lVXfu6NGjNGvWjJKSEi688EJyc3Pp3bt3vKsVd3YemZr4WikLRi6Ic03C51Za9eomvG4SQj0WAg/ilwxRRFKAh4BL8Ewju0lEnscTRO6utP80oCuwUVUfFZGngRoDh3EuOzubnTt3UlpaypQpUyxoGGPK1RQ4NonIDar6mP9KEfk5sKWmwlV1vYikV1rdH/hYVXd7y1oKjFHVu/G0TioQkX388MCh0/nOTYgq38oyxhifmgLHb4AVIjKJHwJFX6AhMDboXtVrC3zmt7wPGFDN9suBB0RkMLA+2EYikg1kA7Rvn/j57I0xJllVGzhU9UtgoIhcBPj6KF5U1dccHDPQ+MigHS2qWgL8vKZCVTUXyAVPH0fEtTPGGFOtUOfjeB143aVj7gPO9ltuB+x3qWxjjDFRFo+nvzcBnUWkg4g0BCYAz8ehHsYYYyIQ1cAhIkvwPCzYVUT2icjPVfUEMANYDXwAPKWqO6JZj0RROa36qFGj+Oijj0Le35co0Bhj4inUlCMRUdWJQdavoo6lLPGlVZ8yZUr5Q3z5+fl8+eWX5ekzysrKSElJiWc1jTGmRsmbqDDJBEurXlZWxkUXXcS1115Ljx6e5ymvvPJK+vTpQ/fu3ask+vvd735H7969GTZsGIcOHYrpZzDGGLDAETPVpVV/5513yMnJYefOnQDMnz+fLVu2sHnzZubNm0dRURHgmUOid+/evPvuuwwZMoTZs2fHrP7GGOMT1VtViered+7lw8Mfulpmt1O78Yf+f4ho3/79+9OhQ4fy5Xnz5pXP+PfZZ5+xa9cuWrVqRb169cpTmmdlZXHVVVc5r7gxxoTJWhwx0r17d7ZsCfywvX8CwHXr1rFmzRo2btzIe++9R69evSgtLQ24X01pyI0xJhrqZIsj0paBE8HSqvtm0vMpLi7mlFNOITU1lQ8//JB//etf5e+dPHmSp59+mgkTJvDEE0+UpzA3xphYqpOBIx58adV/85vfcM8999C4cWPS09O58sorK2w3cuRIHnnkETIyMujatSvnnXde+XtNmzZlx44d9OnTh5YtW/Lkk0/G+mMYY0z1adWTlaVVN9Fi55GpydSXp1JwuCDi+TziOQmUW2nVjTHGhGFUx1ER71twuAAg4WcPtMBhjDEuGtdlXMQX/kinqo01G1VljDEmLBY4jDHGhMUChzHGmLBY4DDGGBOWhA8cInKOiDwlIg+LyNUxPrarT2enpKSQmZlJ9+7d6dmzJ//zP//DyZMnq92nsLDQ5v82xiSUaM/HMV9EDorI9krrR4pIgYh8LCIzayjmUuABVf1/wOSoVTYGmjRpQn5+Pjt27ODVV19l1apVNSYqtMBhjEk00W5xLARG+q8QkRTgITwB4RxgordV0UNEVlZ6tQEWARNE5K9AqyjXt1xeXl75z+np6RWW3dCmTRtyc3N58MEHUVUKCwsZPHgwvXv3pnfv3rz11lsAzJw5kzfffJPMzEzmzp0bdDtjjImVaE/ktF5E0iut7g98rKq7AURkKTBGVe8GLg9S1C+9AWd5tOrqLy8vj+zs7PLlvXv3li9PmjTJteN07NiRkydPcvDgQdq0acOrr75K48aN2bVrFxMnTmTz5s3cc889zJkzh5UrVwJQUlIScDtjjImVeDwA2Bb4zG95HzAg2MbewHMb0BT4azXbZQPZAO3bt4+oYtX1Z5SUlJCVlUVWVhZupmnxlXX8+HFmzJhBfn4+KSkpQaeUDXU7Y4yJlngEjkBX56BXYlUtxBsQqqOquUAueHJVRVq5WNq9ezcpKSm0adOG2bNnc/rpp/Pee+9x8uRJGjduHHCfuXPnhrSdMcZESzxGVe0DzvZbbgfsj0M9qlBVVJW0tLSA76elpbnW2jh06BDTp09nxowZiAjFxcWceeaZ1KtXj0WLFlFWVgZA8+bNOXLkSPl+wbYzxphYiUfg2AR0FpEOItIQmAA8H4d6BJWTk0NqamqFdampqeTk5Dgq99ixY+XDcYcPH86IESO44447ALjpppt4/PHHOe+88/joo4/KJ3fKyMigfv369OzZk7lz5wbdzhhjYiWqt6pEZAkwFGgtIvuAO1T1HyIyA1gNpADzVXVHNOsRLl8HeFZWFuBpaeTk5DjuGK+uddC5c2e2bdtWvnz33XcD0KBBA9auXVth20DbGWNMrER7VNXEIOtXAauieWynJk2aVB44CgsL41sZY4xJIJZWvRq1cZIrY4xxKuFTjhhjjEksdSpwWAvCOGHnjzEedeZWVePGjSkqKqJVq1auJi7ar4E9AAAUcUlEQVQ0dYOqUlRUZM/NmKgrOFwQ8UyA3U7txh/6/8HlGlVVZwJHu3bt2LdvH4cOHYp3VUySaty4Me3atYt3NUwt5mS+8liS2tj87tu3r1r+JmOMCY+IbFHVvjVtV6f6OIwxxjhngcMYY0xYLHAYY4wJS63sHC8sLKRv3xpv0xljjKmodygb1crAkZ6ebpMbGWNMmETk3VC2s1tVxhhjwmKBwxhjTFgscBhjjAmLBQ5jTFISEUsfFCcWOIwxxoTFAocxxpiwWOAwxhgTFgscxpikk5eXV/5zenp6hWUTfRY4jDFJJS8vj+zs7PLlvXv3kp2dbcEjhiytujEmYTkdNVUbr2/RZGnVjTHGRIUFDmNM3NT0LIaqVnmlpaUF3DYtLa3KtiY6LHAYY5JKTk4OqampFdalpqaSk5MTpxrVPRY4jDFJZdKkSeTm5pYvp6WlkZuby6RJk+JYq7rFOseNMXHju00VyXXIyb4mMOscN8YkNHsWI3lZ4DDGxJw9i5Hc7FaVMSaqnDyLUd31yW5VuS/UW1W1cupYY0ztl6wBozYEPLtVFQbL/29M+Jw8i2ESkwUOY0zM2bMYyc3VW1Ui0hQoVdUyl8stBI4AZcCJUO7BGWMSl++Zi6ysLMDT0sjJybFnMZKEo85xEakHTAAmAf2A74BGwCFgFZCrqrscV9ITOPqq6lehbB+tzvHacG/SmERSF/+mEvkzx+o5jteBTsCtwBmqeraqtgEGA/8C7hGRLIfHSAg25twY41RtuY44vVU1XFWPV16pqoeBZ4BnRKSBw2MAKPCKiCjwqKrm1rSDm4KNOQesaW2MCUltuo648hyHiNyhqrNdqE+w8s9S1f0i0gZ4FfiVqq6vtE02kA3Qvn37Pnv37o30WBHXMxGbnsYkskS+beNEsl5HYv0cxx0ikgqcCrwLLFXVf7tUNqq63/vvQRFZAfQH1lfaJhfIBU8fh1vHNsZET20LGHWFW8NxFSgFVgNnA2+JSE83ChaRpiLS3PczMALY7kbZgdiYc2OMU7X9OuJW4PhQVe9Q1adV9TZgDDDXpbJPBzaIyHvAO8CLqvqyS2WHxMacG2Ocqk3XEbduVX0lIn1UdQuAqn4kIqe5UbCq7gZcab1EysacG2Ocqk3XEbc6x3sCS4EtwPtABtBcVUc7LjwC9hyHMSZRJfIcJDF5jkO8n0JV3wMygSXet14HJvpvY4wxpnZweqvqdRF5BnhOVT8FXgReFJGGwCARmYIniCx0eBxjjEkIdufBeeAYCUwDlohIB+BroDGQArwCzFXVfIfHMMaYWqM2BBxHt6pUtVRV/1dVLwDSgGFAb1VNU9UbalvQSKbhcsaY2iWR0pW4lh3Xm3rkC7fKM8YY45Fo6UpsPo4YsUmgjEl+sfjW77tW+L+ysrIoKSmpsF1JSQlZWVkVtosVCxzGGBOCYN/6kzXDrROuPMeRaKL1HIcTNhLDmOSSSIkK09PTCZS4NS0tjcLCQteOE6v5OHwHG+eXT+pPIrJcRHq7UbYxxtR1iZauxK1bVber6hERGQT8FHgceNilso0xJuYiTVIYjbsKkyZNIjf3h2mI0tLSyM3NjVu6ErcCh2+O8cuAh1X1OaChS2UnvUQaRmeMiUy8v/X7B4nCwsK45rhyK3B8LiKPAtcAq0SkkYtlJzXrUDOJzkb8hSbRvvXHk1tJDlPxPEX+vqruEpEzgR6q+orjwiMQr85xp3981nFu4sEGboQnnr+vREly6MoDgKpaAiz3W/4CexjQGGNclSjB3Wl23CMi8k2A1xER+catSiYLJ7N+JcoJYYwxNXGaq6q5qrYI8Gquqi3cqmQyi3eHmjHGXfZFz8UObBE5RUT6i8iFvpdbZScz61CrW5Kto9mNEX/J9pmNc251jv8C+DXQDsgHzgM2qurFjguPgD05buIlmf6ffSP+/HMgpaamhv3FJpk+s6leqJ3jbgWO94F+wL9UNVNEugGzVXW848IjYIHDxEsi/z9HK4VGIn9mE56YphwBSlW11HvgRqr6IdDVpbKNMcYkELcCxz4R+RHwLPCqiDwH7Hep7FrBOtRMvNmIP+MWVwKHqo5V1a9VdRZwO/APYIwbZRuTLJx2NMe6k9mNEX+WTqducuUBQBH5c4DVmcCdbpRvTKJLtBnaQuGrV1ZWFuBpaeTk5IRc32T8zMYdbnWO/85vsTFwOfCBqk5zXHgEErFz3NQeidjJHO19E2luChM9sU45cl+lg88BnnejbGOjVoxxm/1NOROtDLapQMcolW1MXDlJLRNMovcVROMzm+Tl1gyA74vINu9rB1AAzHOjbG/5I0WkQEQ+FpGZbpVrjFucdDTHO/V+pBd4S6dThwX6JhHuC0jze7UF6rtRrrfsFOATPC2YhsB7wDnV7dOnTx+tTQD1/FeZRLZ48eLy/6u0tDRdvHhxlW1870f6ivS40RLPYzthf1OBAZs1lOtyKBsF3Rn+s7qXk7L9jnE+sNpv+Vbg1ur2qU2BI1n/MOuqmi5IbgeOxYsXa2pqaoVtUlNTY3qeJNtF2P6mggs1cDjtHG/u/bcrnpQjvg7xK4D1Dsv2aQt85re8DxhQ3Q5FRUUsXLiwwrru3bvTr18/jh8/HvAWQGZmJpmZmZSUlPDUU09Veb9v376ce+65FBcXs2LFiirvn3/++XTt2pWvvvqKlStXVnn/wgsvpGPHjhw4cICXX365yvvDhg3j7LPP5rPPPmPt2rUAbNy4scLn2Lt3LzfccAPr16/n/PPPr7D/5ZdfTuvWrSkoKGDjxo1Vyh87diwtW7Zk+/btBBpxds0115Camkp+fj75+flV3p80aRINGjRg06ZN7Nixo8r7119/PQBvvfUWH330UYX3GjRoUD4884033mDPnj0V3k9NTeWaa64BYM2aNezbt6/C+y1atOCqq64C4OWXX+bAgQMV3m/VqhVXXHEFAC+88AJFRUUV3j/jjDMYOXIkAMuXL+ebbypm/G/Xrh3Dhw8H4KmnnqqQuwmgQ4cODBkyBPDcVjp+/HiF97t06cLAgQMr/B78/9/8z70FCxZQ2W233cYXX1SdvqZVq1bMmTOnwrk3derUCtssW7asSn1LSkq48cYbWbNmDQDr169n9+7dfPHFFyGfe/5GjhzJGWecwe7du1m/vuqfdatWrSgqKkqKc++JJ56o8jflu014+umnJ/W5V/maB86ve8E4ChyqOhtARF4BeqvqEe/yLGCZk7L9BBoHWOWGrIhkA9kAbdu2denQsSUinH322QwbNgzwXBS+//77CtscO3aMvLy8Kusvv/zymNXTuGvGjBnk5ORUuGg0bNiQn/3sZzXu++2334a1vq556aWXyp9TufDCC9myZUuVv52SkhKysrLKA6iP74JtqnLrOY4PgZ6q+p13uRHwnqp2c6Hs84FZqvpT7/KtAKp6d7B9kvU5DhsrX3fl5eVF9CBeeno6e/furbI+LS2NwsJCt6uZdOxvKjyxTnK4CHhHRGaJyB3A28DjLpW9CegsIh1EpCEwgVr6jEjl+4iWS6ju8A8ShYWFIT95bSObqmd/U9HhVq6qHGAa8G/ga2BqdS2CMMs+AcwAVgMfAE+patUbnbWQXRRMTWyisPDY35RLQulBT7aXjapyD0k2YiaZOfld2/9T6OL9N5XICHFUlaM+DhHZoKqDROQIFTusxXsSx2Xe8WTt4wgmnukRLDVDcrD/p/DY7yuwmOSqUtVB3n+b17StiQ/7AzHGuM2tlCPjRKS59+c/ichyEenlRtnGmJr5biEYEwtujaq6XVWPiMgg4Kd4RlQ94lLZdV68LgqJnnjPmEhZoHXGrcBR5v33MuBhVX0OT14pk6TinXjPGJO43Aocn4vIo8A1wCrvA4DRStluQhRqi8E3Zan/KysrK2Aqi6ysrArbGWPqHrcu7tfgec5ipKp+DZwK/N6lsk0E6nKLwUlQs4BoTM1cSTmSaGrbcNxQuJ1aIZlTWcRrClZjkl1MU46IR5aI/Nm73F5E+rtRtokPe8LWGBOMW7eq/hfPvBkTvctHgIdcKtuEoPKTnU5z8lgqC2NMMG4FjgGq+kugFEBV/42NqoorN1oMkSbe84lHf4GTIcQ2/NiYEIWSl6SmF55suCnAu97l04CtbpQdyas25apywo2cPCRR/iQns+Elwkx6xsQbschV5SMik4DxQG88D/9dDfxJVd2azCksdbFzPBinnb2J3NEcr9FP0fo8xsRbTHJV+ahqnohsAYZ5V40DerhRtokvu0gaYypz1MchIi1E5FYReRBoj6eTvB7wAp5nO4yJmsrN53AGBDjZ15i6zmnn+CKgK/A+8AvgFTy3qcao6hiHZRsXxOtiF4+OZicDAmz4sTFhCKUjJNgLeN/v5xQ8MwA2d1KmGy/rHI+veHY0OxkQYBP8mLqOGE3k9K6q9g62HC/WOR47bj+x7oZE7tA3JpHFqnO8p4h84zsm0MS7HNcZAI0xxkSPoz4OVU1R1RbeV3NVre/3swWNOiBQMzbUjmZjTHKy1OfGdfHuaHYSmCyoGVMzCxzGdZbnypjazdKqm6ixTmpjkktM06obY4ypOyxwGGOMCYsruaqMCcRuMxlTO1mLwyQcmxfDmMRmgcMklLy8PLKzs8uX9+7dS3Z2tgUPYxJIQo+qEpFZwA3AIe+q21R1VU372aiq5OB0Po1EPneNSUYxnY8jyuaq6px4V8IYY4yH3aoyceMkXYm1NoyJn2QIHDNEZJuIzBeRU4JtJCLZIrJZRDYfOnQo2GYmwcU7XYkxpmZxDxwiskZEtgd4jQEeBjoBmcAXwH3BylHVXFXtq6p9TzvttBjV3rjN0pUYk/gSunPcn4ikAytV9dyatrXO8eRnKUeMib1akXJERM70WxwLbI9XXYwxxngk+qiq/xaRTDzTeRYCN8a3OsYYYxI6cKjqdfGug4kPu0VlTOJK6FtVxhhjEo8FDmOMMWGxwGGMMSYsFjiMMcaEJWme4wiHiBwC9tawWUugOITiatquuveDvdca+CqEY8dbqL+jeB8j0jLC2S+UbZ2cK9W9b+eLu+VHUk64+0T7fInWuZKmqjU/QR0oX1BdeAG5bmxX3fvB3gM2x/vzu/k7ivcxIi0jnP1C2dbJuVLd+3a+uFt+JOWEu0+0z5d4nyt1+VbVCy5tV937oR4jUcWi/m4cI9IywtkvlG2dnCuhHiORRbv+bpUfSTnh7hPt8yWu50qtvFWV6ERks4bwWL8xYOeLCV2szpW63OKIp9yaNzGmnJ0vJlQxOVesxWGMMSYs1uIwxhgTFgscxhhjwmKBwxhjTFgscCQAEWkqIo+LyGMiYlPdmWqJSEcR+YeIPB3vupjEJiJXeq8rz4nICLfKtcARJd450g+KyPZK60eKSIGIfCwiM72rrwKeVtUbgNExr6yJu3DOF1Xdrao/j09NTbyFea48672uXA+Md6sOFjiiZyEw0n+FiKQADwGXAucAE0XkHKAd8Jl3s7IY1tEkjoWEfr6Yum0h4Z8rf/K+7woLHFGiquuBw5VW9wc+9n5j/B5YCowB9uEJHmD/J3VSmOeLqcPCOVfE417gJVV916062EUqttryQ8sCPAGjLbAc+JmIPEzyp50w7gl4vohIKxF5BOglIrfGp2omwQS7tvwKGA5cLSLT3TpYQk8dWwtJgHWqqt8CU2NdGZPwgp0vRYBrFwFTKwQ7V+YB89w+mLU4YmsfcLbfcjtgf5zqYhKfnS8mVDE9VyxwxNYmoLOIdBCRhsAE4Pk418kkLjtfTKhieq5Y4IgSEVkCbAS6isg+Efm5qp4AZgCrgQ+Ap1R1RzzraRKDnS8mVIlwrliSQ2OMMWGxFocxxpiwWOAwxhgTFgscxhhjwmKBwxhjTFgscBhjjAmLBQ5jjDFhscBh6hwRKRORfL/XzJr3ig0Redo738bb3rp9KiKH/OqaHmS/v4jIXZXW9RWRbd6f14pIy+h/AlMX2HMcps4RkaOq2szlMut7H8JyUkZ34C+qOtZv3fVAX1WdEcK+K1S1i9+6OUCRqt4tIj8HWqvqvU7qaAxYi8OYciJSKCKzReRdEXlfRLp51zf1Tp6zSUS2isgY7/rrRWSZiLwAvCIi9UTkf0Vkh4isFJFVInK1iAwTkRV+x7lERJYHqMIk4LkQ6nmpiGz01vNJEWnqfUq4VET6eLcRYBye9Np4y73Wye/HGB8LHKYualLpVpX/zGhfqWpv4GHgZu+6PwKvqWo/4CLgryLS1Pve+cAUVb0Yz0yO6UAP4Bfe9wBeA34iIqd5l6cCCwLU6wJgS3UVF5E2wExgmLee24Bfe99egidHka+s/aq6B0BVvwKai8iPqivfmFBYWnVTFx1T1cwg7/laAlvwBAKAEcBoEfEFksZAe+/Pr6qqb1KdQcAyVT0JHBCR18GT21pEFgFZIrIAT0CZHODYZwKHaqj7QDwzvL3laVTQENjgfW8J8IaI3IIngCyptO8h7zG+ruEYxlTLAocxFX3n/beMH/4+BPiZqhb4bygiA4Bv/VdVU+4CPJN0leIJLoH6Q47hCUrVEeBlVb2u8huqWigi+4HBwFigT6VNGnuPYYwjdqvKmJqtBn7l7TdARHoF2W4Dnpkc64nI6cBQ3xuquh/P/Ah/wjNndCAfAD+uoS5vAUNEpKO3Lk1FpLPf+0vwTNzzgaoe8K0UkXpAayrOEmdMRCxwmLqoch/HPTVsfxfQANgmItu9y4E8g2dCne3Ao8DbQLHf+3nAZ6q6M8j+L+IXbAJR1S+BnwNPish7eAJJF79NngLO5YdOcZ/+wAZVLauufGNCYcNxjXGRiDRT1aMi0gp4B7jA981fRB4EtqrqP4Ls2wR43buPqxd4EXkIzxwNb7hZrqmbrI/DGHet9I5cagjc5Rc0tuDpD/ldsB1V9ZiI3AG0BT51uV5bLWgYt1iLwxhjTFisj8MYY0xYLHAYY4wJiwUOY4wxYbHAYYwxJiwWOIwxxoTFAocxxpiw/H84ixcgqRVATQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x109089208>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "resspec = cscripts.csresspec(like.obs())\n",
    "resspec['algorithm'] = 'SIGNIFICANCE'\n",
    "resspec['mask'] = True\n",
    "resspec['ra'] = 83.63\n",
    "resspec['dec'] = 22.01\n",
    "resspec['rad'] = 0.2\n",
    "resspec['components'] = True\n",
    "resspec['outfile'] = 'resspec_2.fits'\n",
    "resspec.execute()\n",
    "\n",
    "plot_residuals('resspec_2.fits','',0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you actually have an excess at the highest energies. In spite of the significant likelihood improvement the model is not quite perfect yet. Is there perhaps a new spectral component appearing at high energies? Test if the addition of a power-law spectral component with a hard index (2) improves things. To do this, add the hard spectral component to the exponentially-cutoff power law using a composite spectral model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "expplaw2  = like.obs().models()['Crab'].spectral().clone() \n",
    "newcomp   = gammalib.GModelSpectralPlaw(1.e-18,-2.,gammalib.GEnergy(1,'TeV'))\n",
    "comp_spec = gammalib.GModelSpectralComposite()\n",
    "comp_spec.append(expplaw2)\n",
    "comp_spec.append(newcomp)\n",
    "like.obs().models()['Crab'].spectral(comp_spec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's fit this new model to the data and check the likelihood improvement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "311.12034105171915\n"
     ]
    }
   ],
   "source": [
    "like.run()\n",
    "like4 = like.opt().value()\n",
    "ts    = -2.0 * (like4 - like3)\n",
    "print(ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The likelihood improvement is again significant, but is this the best possible model of the sky?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises \n",
    "\n",
    "- Test other spectral models for the Crab nebula, such as log-parabola and broken power law. Can you say what the best spectral model is?\n",
    "- What can you say from the residual maps about the spectrum of `Src1`? Try looking at the residuals below and above 300 GeV, or calculate the residual spectrum at the position of this source.\n",
    "- Replace the 2 sources with a single extended source (disk). Can you say which of the two hypotheses is describing the data better?"
   ]
  }
 ],
 "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.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}