{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# How to generate likelihood profiles for a fitted spectrum?" ] }, { "cell_type": "raw", "metadata": {}, "source": [ "
\"Download
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this tutorial we will learn how to generate likelihood profiles for a spectrum using `csspec`. This is done by running `csspec` with non-zero values for the `dll_sigstep` and `dll_sigmax` parameters. As we'll see, the results can then be plotted using the `show_spectrum.py` example script." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The preliminary data/fit\n", "Before we can demonstrate how to generate the likelihood profiles, we need to first simulate some data to work with. For this example we will simulate four hours on a simple Crab source. Once we have the simulated Crab data we will fit the model back to the data. This fitted model will be used to generate the SED and likelihood profile in the next section." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Load the modules\n", "import gammalib\n", "import ctools\n", "import cscripts\n", "\n", "# Define the ctools install directory\n", "import os\n", "ct_dir = os.environ['CTOOLS']\n", "os.environ['CALDB'] = ct_dir + '/share/caldb/'" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Configure some preliminary variables\n", "inmodel = ct_dir + '/share/models/crab.xml'\n", "caldb = 'prod2'\n", "irf = 'North_5h'\n", "emin = 0.1\n", "emax = 100" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Simulate some data\n", "sim = ctools.ctobssim()\n", "sim['inmodel'] = inmodel\n", "sim['caldb'] = caldb\n", "sim['irf'] = irf\n", "sim['edisp'] = False\n", "sim['outevents'] = 'outfile.fits'\n", "sim['prefix'] = 'sim_events_'\n", "sim['ra'] = 83.63\n", "sim['dec'] = 22.151\n", "sim['rad'] = 5\n", "sim['tmin'] = 0\n", "sim['tmax'] = 7200\n", "sim['emin'] = emin\n", "sim['emax'] = emax\n", "sim['deadc'] = 0.98\n", "\n", "# Run the simulation (prevents writing to disk)\n", "sim.run()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Fit the data\n", "fitter = ctools.ctlike(sim.obs())\n", "fitter['edisp'] = False\n", "fitter.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generate the likelhood profiles\n", "Now that we have a fitted source model and observations we can generate the sample of likelihood profiles." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Configure csspec\n", "sed = cscripts.csspec(fitter.obs())\n", "sed['outfile'] = 'sed_likelihood_profile.fits'\n", "sed['caldb'] = caldb\n", "sed['irf'] = irf\n", "sed['srcname'] = 'Crab'\n", "sed['emin'] = emin\n", "sed['emax'] = emax\n", "sed['enumbins'] = 10\n", "sed['debug'] = True\n", "\n", "# Parameters that control the likelihood profile generation\n", "sed['dll_sigstep'] = 1\n", "sed['dll_sigmax'] = 5\n", "\n", "# Execute csspec & save results to file\n", "sed.execute()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have a fitted SED, we can generate a plot from the results. A simple way to accomplish this is by using the `show_spectrum.py` example script. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEWCAYAAABBvWFzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXm8HGWVsJ9T3X2X3OwJJJgNZFVkX5UtgAiMMm7IIiqg\n4iDq4OACAjMEGEBAwW/ADUVQFBBUBAZR4khIInvYwhIIS1YgJJD1rr2c74+qm3Qu93ZXdVdVd3Wf\nJ7/6pavq3e57+9aps7znFVXFMAzDMKrFqfUADMMwjMbABIphGIYRCiZQDMMwjFAwgWIYhmGEggkU\nwzAMIxRMoBiGYRihYALFMAzDCAUTKIZhGEYoJF6giMj7ROT3IvJjEfl0rcdjGIYRBiJylIgsEJGX\nROTsWo/HD5L0lfIichbwiKr+U0TuVNWP13pMhmEY1SAiDvAScDjwOvAYcIKqLqjpwMpQNxqKiFwv\nIitE5JkB18tJ6ZuAE0TkCmBsLIM1DMOIln2Bhaq6WFWzwK1A3b8s141AAW4Ajiy+4Enpa73rOwMn\nishO3r3Pi8hVQFpVvwGcA6yKd8iGYRiRMAlYWnS+zLtW16RrPYB+VHWuiEwbcHmjlAYQkX4pvUBV\nbwJuEpFpIvJzYBhwZayDNgzDiAYZ5Frd+yfqRqAMwWBSet/iAp6w+bdyDYlI3f8yDMOoH1R1sIe6\nL0alM7ounwtSZYWqTiw6XwZMLTqfjOtLqWvqyeQ1GKFKae1cE8txwblnx1a3XPlS94e65/d6ufN6\nneOg9f2UDWuebY7jn+PBrlXLunyO051Rvg9gwoAmHgO28ywwLcAJwF1VDyxi6l1DSaSUnn7QgbHV\nLVe+1P2h7vm9PvB80eIlJccSJtXMcdD6fsqGNc82x5WVqWaO/fYfFCkUKq6rqnkR+TpwH+6L//Wq\n+kJYY4uKugobFpGtgbtVdRfvPAW8iBs69wbwKHBiJRMrIhrGm4cxNKd85avceN1Paz2MhsbmOB6k\nY3RVJi8R0X9PjfZd/n/ya6rqr16oG5OXiNwMPAjsICJLRORUVc0D38CV0s8Bt1YjpWdcchmzZs8J\nZ8DGuzjlc5+t9RAaHpvjaJk1ew4zLrkslLYc8X80CnWloUSJaSiGYfglDA3lrLR/DeWqnGkohrEZ\npv1Fj81xcnBEfB+NQr075UNlxiWXMf2gA5l+8EG1HophGHXIrNlzmDVnbihtNePbupm8DMMwBhCG\nyevs9Cjf5S/PrW0Ik1dTaSiGYRhxkU6l/BcOtAayfmlGrcyICLPvR4/NcXJwAhyNgmkohmEYEdBI\n4cB+aSqBYk75aLF5jR6b42gxp3x1mFPeMAxjAGE45S9uH+e7/H92v90QTvlmFKJGRJh9P3psjpOD\n+VAMwzCMUBCtPDlkUjGTl2EYxgDCMHldOWy87/Lf6VrVECavptJQzClvGEYpQnXKJ148BMc0FCM0\nZs2eY8I6YmyO4yEMDeXqDv8ayn90BtNQRORYYAbwPmAfVX0i8CAjoJH8QYZhGHVDxOnr5wOfBB4I\nddBV0lQmLyNa7M05emyOk4Mz6A7m4aCqLwKI1FeqYhMohmEYEdCMPhQzeRmhYWskosfmODlIoeD7\nGLS+yEwReabomO/9f0zMP4pvmkpDsSgvwzBKEWaUV6ZEtuEXc328mOsrWV9VjwhlIDFiUV6GYRgD\nCCPK6/pRW/ou/6W1b1XUn4jcD3xbVecFrRsFZvIyDMOIAAlwBG5b5BMishTYH/hfEbk3lEFXiQkU\nIzTMvh89NsfJIcqwYVX9s6pOUdV2Vd1KVY8O/ycITlP5UAzDMOIiyrDhesV8KIZhGAMIw4fy2zET\nfJf/3OoVlsvLMAzDGJyhwoEbmabyocy45DKzQUeIzW302BxHy6zZc5hxyWWhtJVKpXwfjYKZvIzQ\nsMSF0WNzHA9hmLx+P3ai7/LHv/NmQ5i8TKAYhmEMIAyBclsAgXJcgwgU86EYhmFEQOKlQwU0lQ/F\niBaz70ePzXFyEBHfR6NgGophGEYENOPbuvlQDMMwBhCGD+WOAOtQPmnrUAzDMIyhcBooHNgvzaiV\nGRFh9v3osTlODlEmh6xXmkpDsf1QDMMoRZj7oTSSoPCL+VAMwzAGEIYP5Z4t3uO7/EdXvh6oPxG5\nAjgG6AVeAU5V1XWBBxoyZvIyDMOIAAnwrwLuA3ZW1d2BhcD3Qh18hZhAMULD7PvRY3OcHKrdU74U\nqvp3Ve2v+DAwOdTBV0hT+VAMwzDiwknF9r7+ReDWuDorhflQDMMwBhCGD2Xmlv6VhiPeWvau/kRk\nJlC8mEUABc5T1bu9MucBe6rqpysda5iYhmIYhhEBpaTRU309PN3XW7K+qh5Rsn2Rk4F/AQ4LPrpo\nMB+KERpm348em+PkIDL0sUdrG6eMGLXxCN62HAV8F/hXVS0tmWLENBTDMIwIiHgdyjVACzDTSy75\nsKqeEW2X5anIhyIiHUCPqubDH1I0mA/FMAy/hOFDuX/CFN/lD12xtHlyeYmIA5wAnATsg7uYplVE\nVgJ/Aa5T1YWRjdIwDCNhOGp7yg/F/cC2uItnJqrqFFXdEjgINwb6+yLyuYjGaCQEs+9Hj81xcnAc\nx/fRKPj1oXxYVbMDL6rqO8AfgT+KSCbUkRmGYSSYClfAJxpfAmUwYVJJmVpjySGjxeY1emyOoyXU\n5JDNJ0/KO+VFZCs2D1g4TFV/G+moIsCc8oZh+CUMp/xDW031Xf6DbyxpCKe8H+PdPsC1wJeB04Cj\nIx2RkVjMvh89NsfJwfZDGQRVvUtEHlHVFQAismX0wzIMw0g20oQ2L8vlZRiGMYAwTF6PTvS/DmXf\nN5toHUo/IrK/qj4c1WAMwzAaBdtTvjwjIxmF0RCYfT96bI6TgyP+j0bBcnkZhmFEgDSSpPBJUIHS\nfDNk+MbWSESPzXFyaEKffGCBMj+SURiGYTQYzShQgvpQDhCREQAicr6I/ElE9oxgXEYCMft+9Ngc\nJwcR8X1U0PZFIvK0iDwpIn8VkYkR/AiBCSpQ/lNV14vIgcCHgeuBn4Y/LMMwjGQjWvB9VMAVqrqb\nqu4B3ANcEO7oKyOoQOnf/+SjuCnr78Hd5MUwzL4fAzbHycFJp3wfQVHVDUWnHUBd5MoP6kNZLiI/\nx9VOLheRVhK0jXBcizibcYWsYRib40T8HBCR/wa+AKwBDo20M58EWikvIsOAo4D5qrrQSxy5i6re\nF9UAw0JENL/A35rMxUuX8+v/dx2Ft1bibLkFJ5/5FaZNmeSvn/dsh/zzIWjCN8lZs+fYG3TE2BzH\nQxgr5RfssJ3v8ju99PK7+hORmcCE4kuAAuep6t1F5c4G2lV1RqXjDYtAGoqqdgF/Kjp/A3gj7EFF\nRee3v1m2zJLOHq5/5AUu7uqlA+gE/vMfs/jSfu9jakdb2fodP78FmTO3KQWKYRibKGWpeKSri0e7\nukvWV9UjfHZ1C64fZYbfsUVFU+Xyum74uLLlHutay9WFHB1F1zqB/3DS7DNsVNn6J81/mGG/+CWc\n973KB2sYRk0JQ0N56X3b+y6/wwsLA/UnItup6sve528AB6nqccFHGi6JWykvItsA5wEj+yfQM8X9\nBHev+wdU9ebB6p624e2y7S+DzYQJ3vlWhZyv+myzPRx4gDtCwzCaloh9qd8XkR1wnfGLgdOj7Mwv\niRMoqvoa8GURua3o8qeA21X1HhG5FRhUoPy0o7yGsrx7LZ2DaCjLnTQ/bS+voZz8rKehNCFm348e\nm+PkUGE4sC9U9djIGq8CXwJFRDpUtVNE0kBBtfqZEpHrgY8BK1R116LrRwE/wo0eu15VL/fR3GTg\nGe9zfqhCR21TXqDs0juCcxct49JsbqMP5dxMmlO2nsyk1vIR0i2pxAS9GYYRIZWEAyedsgJFRL4L\njBcRB7jMO74SQt83ANcAvynqy8HdHfJw4HXgMRG5U1UXiMjngT2AK71ggGJ9cimbhMqQeubUn11a\ndlBTgfe8voIrfnEr+vZqZNwYzjztBKa9Z0LZugDOyOFw0IG+yjYa9uYcPTbHySHqsOF6xM+e8ocA\nDwNZ4FjgI6r65VA6F5kG3N2voYjI/sAFqnq0d34OoMVaioiMBS7BXQvzS1W93POhXAt0A3NV9ZZB\n+rINtgzD8EUYTvnFu+3ku/y0pxc0xAZbfuwzncApqlpQ1duAf0Q4nkm42kY/y7xrG1HVd1T1q6q6\nfb+gUdUuVf2iqn5tMGFixIPlmYoem+PkEGUur3rFz57yjwOPw8YdGwd1eIfEYDMbWlzzKV/5KltP\nmwrA6FGj2H3XXTaaEPr/UO288vOnnplfV+NpxPN+6mU8jXL+o2t/wlPPzN/4fAiDpMgJETmr1H1V\nvcp3WwFXyh+pqn/zXaF8e4OZvGao6lHe+btMXlX0FYvJK+51PY30dmMY9UIYJq9le77Pd/nJT7xQ\nM5OXiPQnltwR2Ae4yzs/BnhUVT/nt62gYcNhPy2FzbWSx4DtPEHzBnACcGLIfUZL1zo01xdLVzJi\nXHJegwyjyXBSyfjbVNULAURkNrCnqq73zmfgrsD3Tc12bBSRm4HpwDgRWYLrjL/BW/V5H5vChl8I\nq88Zl1zG9IMOjDRSpufcr8NLr1HYovyaFYAlnd3c/NwStLsPaW/hsztPZWpHu6+6w078Ohx9dDXD\nDRVbIxE9NsfRMmv2HGbNmRtKWwmM8poAFL8N97F5LrGy1GzHRlX97BDX7wXuDaufYmbEkA7l4bmv\nMXnxa7w8bETZsityOR5Zv5ofFPIb17x8e/lq9hsxhgnp8r+aIyY/SLqOBIphJJ3pBx/E9IMP4sJL\nq7ayJ9F48BvgURG5wzv/BHBjkAaCJod8PUj5ZuSf67rZryfLrL6usmWX5zv5OfmNq/I7gB8U8vzb\n2tVMSg1MAPNuDtP6SnVgb87RY3OcHJLm31TVS0TkXuAgXPfGqar6ZJA2Aj2PRGRv3CxV07y64o5j\n00r3ZidbULbRHJfn15UtewGD5w3bljwX+qife/iRisZoGEb0JEye9JPHzQ+mVLBpV9AX3N8B38E1\nfdXFDmFBiMOHss+INrJr2vlT+/CyZZd1rqWzr+ddecOWtbTxp47yPph//eB+lQ80Asy+Hz02x9ES\npg8laRqKiJwJnAb8EVdZ+K2IXKeq1/htI6hAWamqd5UvVp/E4UM5YPtxtPWsZMLULcqW3bZ7JOc/\nv4j/7slu9KGc35bh398/hSntrWXrp5xkfWENo94J1YcSYXLIiPgSsJ+qdgKIyOXAQ7gpsnwRdB3K\n4bhhvP+HmyoeAFX905CV6oS41qHkn/4HPP4M7PkBX+UXv76CXxflDTs5SN6w1Q5y2GHVDNcwjEEI\nYx3K6oP9ewLGzH6m5qlXRGQ+sI+q9njnbcBjqrqL7zYCCpTfAjsBz7HJ5KWq+kXfjdSIRszlZYso\nDSMawhAoaw7ZzXf50Q88XVF/IvJt4ApgvKq+E7T+gLbOAk4G7sA1eX0cuFFVf+S3jaAmr31UdceA\ndeqGOHwosbL6zfgWUY6fDFI6HbfZ96PH5jhawvShELFJWkQm4ybJXRxGe6p6lYjMAvrTpUcb5QU8\nKCLvV9XnA9arC+LwocRJ4aJzKIwvH17cz+I1G7jpgacpbOjGGd7O5w/ZjWmjywcPAKT3+TgceWSl\nQzWMRBCmDyWGMK+rcYOkwvRr53EjvGKJ8tofeEpEXsP1oVjYcA1Z+eeHWeJkfJV9I5dl5ppVXJ7f\ntIjy7BeXccTo8WyVLt/GXl3jSJURKPbmHD02x8khShOxiBwDLFXV+WH1U4sor6MCljci5K1sjiey\nOV9l53Wv5UeFzRdRXp7P883Vb7OXj62N9wCab/85w6iCKh/0IjKTzVOfCK7mcD5wLnDEgHvVUnWU\nV9CV8qHY6oxwGN65gdOzveULMvQiyvcUcpze+XbZ+rmHyi+iNPt+9NgcJ4mhg2bmvLOeuavXl66t\nesRg10XkA8DWwNPiqieTgXkisq+qvlXxcF2hVLyFep6AgiroSvlfA2eq6hrvfAzwwyREeUHjOeWX\ntg7jZ4U2X2Vfz2+gsyjNC7hmr1dI8d1UeT/KxfvvV1dpXgwjCkJd2JgZ+i/m4AljOHjCmI3n33/t\nTd/tquqzwMSN/bguiD1VdXVFA93EDcAjXi6v/iivXwVpIGjY8JOquke5a/VII4YNz5+0LY/n/b1A\nrMrnWNi9jqu1sNGH8h/isH37SManyouKz33xRDKXXlzdgA0jIYQRNrz+X/xnshjxl0cq7k9EXgX2\nrjZs2GtrTzZFec2JOsrLEZEx/ZLQ29/dXlxrxBaTt2S6+vOhACzvHcWMle8g2RyaSfOFLcYyqbXF\nV13nQ/vHuu7F1rwYiSemTBaq+t4w2hGRVtx1hsNxn+vHiMgxqnqR3zaCCoMf4oYO/wHXQHgccEnA\nNoyQGH/8QdC9wXf5KbhhepUgO09F33i5ZJlZj8xj+n57VdhDUV8TtgEfWlMzYj6UBJG8l6I7gbXA\nPIoyoQQhqFP+NyLyOHAYro3tU0ldk9IIpE4/L7a+8tfOQNeXzoBcWLqS/BMzq+4rffZVMNzfBmWG\nUa8kUMue3L/9eqUEfg30BIgJkTpAhsX30F1z5yN05Uub13YA3ljgOheX9fbx5xVvI7kcmk7ziQnj\nmOzTvDZpv7mkPvbRaofckJh2kiCSl7z1QRHZRVUr3kixqewKjRblFSdvL36Lpwv+/kBW5nPM71zL\nD4sCAL61rpNdOkaxhQ9T1qcefBhMoBjFxOS/mzV7bnipV0qEDdcTXlJIxZUHp3pO/ooWrgeK8koy\njRjlFSf/nLgNf+spnYlhUSHH1k6aRbkN/FTfHaL8VUmxdbp8iPL5X/0CLZdZRNlgNKsPRfO52IQK\ngDNqi6qjvDqP9f976vjDnJplGxaRaaXuB1l/GHQdylmDXF4LzFPVp4K0ZSSLid1dXJQr7aebBUzP\nD72IchvNc2F2bdm+bCdK4130bIBCsvYXkYSYvMJcsB7U5LW3d9ztnX8MeAY4XURuV9UrwhqYUV8s\naxvGz7vLL6L8C7aIMkqaUTsBoLcbLeTLl6snEuKUF5G5qnqgiKxncztdv8lrpN+2gv7dTsZdkbnB\nG8gFwD3AwbihZiZQGpTxmRS7OP6+LlvlR3JW11quKvKhnCUOew3zt4gyIS92RozohjVQJiik7kjI\nF1lVD/T+H1FtW0EFypZA8QYcWWCCqnaLSEVxy0YySI0axXbdpX/FT/T1sGdLG9vRyrS2Fs7uXEc6\nnyeXSnFkx0hfWY0B2N//CuNmo1l9KPR0xbb3T1gkMGy4aoIKlN8BD4vInd75McDNItKBhRI3NNMO\nfx/TNpT2f6x7aw27bzkagN2BoyvsK33IARXWNBoV7VoLPhOh1g0J0VCKTF3FA+4/j8bk5WW1vBHX\nTH6g19npqvq4V+Qkv23VCgsbrpyWIw6FbE/JMh8Oq7PWlvgiehL2Ftm0393uTugr/f0Lg1nPvMgD\nz7wYUmvJiKANw9TVT9DkkPODbFhfT1jYcHLQvp54Inpa25vSLJFE8jNvht6u2PpLH//tqsOGe07x\n/4rVduPfaxY23I+nNJwEbKOqF4vIFGArVX3UbxtBTV5PiMg+qvpYwHpGExCWfb/wzGy0rzuEEZUm\ntfeR0OIv/X+90Mw+FHriEyhhkMCXlZ/gbvt7GHAxsAH4MbCP3waCCpT9gJNEZDFu8I5tAWyEjv7t\nLnRdDNrkLtMTJ1Cals7OQIlQ64IIfShehO1pQP+GWueq6l+rbHY/Vd1TRJ4EUNXVIuIvX5JHUIFS\nelNxo6kJ6815w9OLyb9d9dYOZRmTzYWyb2qcNKV2AtDbDT3Ra62hEr2GcpWqXhVie1kRSeE5f0Rk\nC1yNxTe2BbBRdzz36lr6Vla7+Vx5DsgrgV6/jNrR2YV2JktDiWGlfNgd/A9wB7CliFwCHIu7f71v\ngqZe6XfavFdVLxKRqcDEIE4bo3EJy76fX7uG1/v8rTlYlc/xZG83rZqnV1Ls0drua/EkQF4VZs+B\nBL31N6sPRft6oTdhYcPRayhfE5HPA48D31LV8nmNSqCqvxORecDhuMLqE6r6QpA2gpq8ip02FwHr\ngT8SwGljGOVIrV/Pi2USUQIMz/WySHv4MXgr8nN8J9dLr7SxId1atn5egTlzEyVQmpbubrQrWU75\navUHEZkJTBjQogLn4T6LL1JVFZH/Bq4CvlRlf3up6jxgQdG1Y1T17hLVNiOwU75ap43RuIT15txX\nUNbnyguUJzXHrWxKRNkBXAmcoDm2y5VflV/AW+4S05qXMKJ+mlE7AdCePrQnWRqKpFND3ntg2Soe\nWP52yfqqeoTPrn7BpvyK1fALETm5fz8UETkR+GaQtoMKlKqdNoZRjml93Vyt5R8eQ2U13pMcF+r6\n8h29d3t0/33Qb3yxkmEGQkaMBxn6AWOUJt/Zg26IfmFjqJR4gThkyhYcMmWLjecXP7YwYNMyUVXf\n9E4/BTxbwQgHcizwBxE5CXfx+heAjwRpIKhAqdppU0tspXy0hGXfX9Y6jOu0fDjvgnwXnZp9V1bj\npyXDOalhZeuf84cfMPKOv6IvPVn5YH0iux8KTvUCpVl9KIWeLNodfS6v2SvXMmdlVa6ITUTrQ7lC\nRHbHfaFfBPxbtQ2q6qsicgLwZ2Ap8BFVDRRa50ugiIioy5BOm/4ywX6EeJlx3vdqPQTDByNSDlNa\nnbLltsm2cHZflsthY1bjs4GjMy3kMuW/2ql5/4Rli9CH7q92yOV5/wGQKe/XSRyqsaSv0e4sha7o\nBcqBHe0c2NHOZQuWVd9YhPOiql8Iq62iHRvBfa6PBVLAo95j3fc6Q78ayv0i8kfgTlVdgOe0EZEW\nETkMOBm4HzfXl9GkhPXmnB0xnMldPswbrcNxcm2c2dNJSyFPn5Niv7YOtkz7/Fq/uphcBvSlV6sb\nsA9S+XD28vA1x3G/18WQJqfQm6XQk428n1Bxyr8U1Qkf8/4f+MWRQa6VxK9AOQr4InCLiGwDrAHa\ncKXYfcDVtmOjERYto8cwKVfaYdnPpJYW9hhW3rw1GH2LVpJPCbzyZvnCVZLJF2JbRKm53viESj4f\nS1+5nr7kCZTk8LtSG2wB4WYbVtUe3DC1n4hIBhgPdKuqZVs0NhKWfX/bHceiE6N//C5+eTV5H9Fk\nYbBzvkAY76u+5njNW/HtbpiOJ8izry9Hvi9hOzamkhGEUcsNtlDVLPBGtR0bxlC07/pe6PIRpVUl\n7yx4mL7ueN56CzFaoXTV6/FtRjVibCxJELO9efJ9CQsoTY7JKzRs624jNMLyocgBh0AMD8RVtzxK\nTzYugRKORPEzx7piqZv7Kg62yKIxPDj7erPkkqahJCTbcE022DKMuHB2OTCWvYneLFzG2O4eljrl\nF0FOKWRpyfVxPzk6tECnOBxKmr50i6/6BYgvzcvqd6CnM/p+ANraQwmHLkdfb55c0jSUhAiUMDfY\nMoFihEZYPhQZPiaE0ZRnXQF2ymWZ78PWvaQAnZrlWtQNUdYCXydLR76FkVr+QVdQ0Nlz4aADqxrz\nrNlzmX5wmTbefgvtjN5kCED7MCQGX0GuL0c2awIlLgYsnPSN33Uo31XVK7zPn1HV24vuXaqq5wbt\n2DBqTXehQFaV7nx5dWix9nCrJ0zAXfdyLcoJhR620/JRZoUN70BfF7reX/TaUGj32rJt6OrVsQkU\nGbs2FudzXy5PNqYAitBItg/lL8CeQSv51VBOAK7wPn8PuL3o3lGACRQjcSu4uwvKe8lxta4rW7Z0\nmpfy9dl1P3SPD6DHV7el0CFTRqNLF5QsU1ixCu2Kx+SVGvMOlMhZFRbZXJ6++l433WhUpF75FSgy\nxOeKOzaMWtPmCK+S5idO+TQvrxS66CT3rjQvT5Dmm055DeWC689h1L2zYdkrlQ/YJ9l31qPd8Tjl\nU52daBwCpVAgG2eoXBgkJGx4CH5RSSW/AkWH+DzYudGkJC3PVLsjZEQY5mMjpGm08fVC5yYfCvB1\nhGlOm6/6svhVWP02+mqwJIADeeDlZRyy3eSSZfpWrKbQE0/YcMuatUgmHoGSNA1FEmbyEpGzBjlf\nC8zzu3Ddr0DZTUTW4Woj7QPCzBKzKbclhzSKaXMclqcztPv4w99J8rRqhu+QY7gW2CAOR5GmNwWL\n/Tw4Vq2iMK4Dfeut8mVLoGvWoG+VXky4YV0f+d54wqE7NnThxCFQ8gWyPoIfqmV+to/52ZCEcfKc\n8nt7R3+6+o8BzwCni8jt/X70UvhdKZ9o3a0fSw4ZLUkT1JNb0/S0tzPNV+k0vbSzf9GV/gT7vuq/\nuYpsaxqWVydQDkhBrkwbG9b3xLZmY/yGXlKZGNahFAr0xaCg7JhuYcd0C7eGEXadPIEyGdhTVTcA\niMgFwD3AwcA8NvnRh8RvlNdZpe6r6lV+2jGMemJUyqEtBvs/QG7leiQdjwmkuytHLqaIqHxXH8Qg\nUAoFDW1xaGxELFBE5BvA14AscI+qnlNlk1sCxepZFpigqt0i4mt3M78mr/6FLzvibvd7l3d+DGD7\nyRtA8nwo4zNpsj7S3IfB+pWdOD58LeV4aEMXHxxeOgigM5cn5yMUOgxyXb0Qg6AsqCZvJ78IfSgi\nMh33+fsBVc2JyPgQmv0d8LCI3InrzvgYcLOIdADP+2nAr8nrQgARmY2rEq33zmfgqkSGkTjGj2kj\nr+2x9LV+Qx+FEKKUNvRkWZ0r/bLYnSuQj+ltPt+bxclFb9rJq8b2MyWErwLfV9UcgKquqrZBVb1Y\nRP6Cu1ujAKer6uPe7ZP8tBH09WwCm6tEfd41w0iUdgIw7j0j0Y54fA0r5q8gm62+rx0kzbreXMky\nnYVCbOah7u4cuVRcAiXybsIl2rDhHYCDReRSoBv4TtHDvxpyuJmCFNfkFYigAuU3uLt43eF1+Els\nUy0jobRNGw8byufhCoPu+SvozcdjtOmJUaBkswXIx5BtuKDkkqahVGnyEpGZbP7C3r8/yfm4z+7R\nqrq/iOwD3Aa8t8r+zgROA/7o9fVbEblOVa/x20YggaKql4jIvUD/q+ipqhr9htxGIkiaD0V23AF6\numLpa8PT1PhSAAAbh0lEQVTdz9ITwq6Nz+X62LnMHiTr8zFqKD158iH4hsqRV5KnoZRwys96ZTkP\nvLK8ZHVVPWLopuV04E9eucdEpCAi41S1mtw+XwL2U9VOr4/LgYeAcAWKiHwQeNjbV/4J4IkKBmsY\n9cW0bSGsNQdlWI9DdwgaSldeWSel21mXUwoxrTfuy+YhBoFSKBDbzxQaJQTK9O0mM71ogerFMx8L\n2vqfgcOB2SKyA5CpUpiAq5UUv/XkCZgJxa+GcjLwYxF5Cfgr8NdKMlEajU2StBMAZ8e9IYbFcgDv\nqENnCNlyx5FiVZl28igaU0aknt48hRhWhGe1QC5h8iTi5JA3AL8Skfm4S6K+EFKbj3guDYBPANcH\nacBvlNfpACKyE3A0cKOIjALuxxUw/1TVhO1+YzQ7Mu49sfXVi0OPzyivrQtZtnaDd9i6kGOR4/6Z\nLpI0i3zsveIST199qq76EDE5JXk+lAjXoXg7534+5DavEpEHgANwNZPALo2gPpQFwALgahFpBw4F\nPgNchbtk32hikuZDiZM+Vd6Tz/KaU/5P7kVJ85gKy7WX3cjztDpMklZGSorF+SzTyrTx3kLOVz8A\nL0mal8Qte2lhLdc5wzfd9PEA7y0UUIlpHUrC5EkSUdV5uKviK8KvD+W/StxeqaomTAyjBDl1tYGF\nlA8lXVfI00sXv9u4mVeWr2mOXGEYeVFyZZ6sW2uWl7WyBZtBHd9ZVSQG30ZO3XjWRJFOxv6FRbkZ\n33WLiLYAHiyxTQduVMA44CK/HRqNi2knQ9NbUHJAj49n73J6+f2Azbx+jHI8vWzDsLJtjNJgmXnX\nF/KsoJcLgOfznUyglRE+t/Ud2dvL2tZW331VSl4byylfT8S+BbCq/rD/s4iMAM4ETgVuBX44VD3D\nMFxyqmxHjmuq2MxrL5+beb0Nvh3Y6zWP0rVxN8pOcpxBntWFYYyQ8kJlXK6PVZnSYcxhkEN9ayhT\n81mmFdzSUwtZlni+oMVOmiWpeNYdAYkRKGHiWycTkbHAWbhL8H+Nm4JldVQDM5KH+VBKs5A0P6L8\n2/xiugfdzOtx0nycDJPK/Nl+k07f61BW0sttA7Shn6AcRy8dPtLS5NU1e0VNQfGdy2uRk9kYUHB+\nbw83pYssNnEqOTH4luoNvz6UK4FPAdcBu/SnN64FIrINcB4wUlWPG+qaYdQTLY6Q8v4vx8RCK2eQ\n5ydFm3mdgTCRVtaK0lLmzXdcocCP8ben/FDa0N7kuNBHG6/n0jxRKL9jZbXkoeLFmjXLUhzD+px6\nw6+G8i3cWOfzgfPE/UL3z1Ygp021qOprwJdF5LZS14z4Me1kaFIivOZkSPswg4xJpVlX6OBEetlD\nszwpGSbSykgnxRgffS0kzQXi7yH/mg6+tfHjpPm6jza+4PTRG0P4VU6VStcl1MyZn3ANRUQOAr6p\nqp/2W8evDyX0mRGR63HTI69Q1V2Lrh8F/AhwgOtV9fKw+zaMuGkRWOqk8ettGJ9KMZ5hXJhby7mp\nYBrAGhzSPl+OJ+rQ2pCfNgr499dUQxCT12B1DX+IiAMch+vecIDtg9Sv5QZbN+DmiPlNUT8OcC1u\nSoHXgcdE5E5VXSAinwf2AK5U1TcYPCVA8+mYdYT5UIYmLUKmQhNIcb1FhSxbl1lw6FcTgtLakB+W\nOOlYFhwm0uSVkLBhABEZCfwbcAZuaq1vqeocEVkapJ2abbClqnNFZODuqfsCC1V1MYCI3Ap8HFig\nqjcBN4nIWBH5KbC7iJytqpd7AQOXFF+rZEyGERVtjvhOojg1n2WqF6W0WFJ8ON8DuA/v10VoK9PO\nEicdaMXy2FSKsZ42dE5AbWiRk8GJ4XmtqhX702umoCTE5CUiVwHHA7cDh6vqq0W3A01fvW2wNQko\nlojLcIVM8Vjewd1cpuQ1I35MOxmakSmHTMrfA2ZdqpVnvWiwZwfc29NH/a585bsbpgKGuqpCQWKI\n8qJy01XNTF7JCRu+C/dZ++gAYRKYettga7DfQGhfh1O+8lW2njYVgNGjRrH7rrtsfAjOmj0HwM7t\nPJLzJX3d9JHj/V7q+edz7p9RFOcrs3leKbh7I/Wbxxb5OJ8FG/0mfsoD7EgGVHjVO3+vdz/8c1dj\nm+KliVmq/s7BNZeVKz+v0MtbmmdUmFpFDEkzw0BVZwGzRGQHETkNWAncXUl+RtEA9kUROQ/XYVO8\nwdbvVfWyoB177U3DHfiu3vn+wAxVPco7Pwc3iqxqE5aIqHauqbYZowTmQxmae3b/EL2rqs0uDs9m\n+/hAmYWEL3T2VhTZdEHfGi5sGR2oTpB38GmFLFt7gmFaIcfi/kSUTprFZfxCa3OVrZO/qrCOs5zg\nQahXF9ahqhWrGCKiud/8t+/y6S+cX1V/YSIiY4B/xX3GX6mqvpWGWm+wJWz+nXwM2M4TNG8AJwAn\nVtG+YdQF7Y5DKoQ31jZHGFamnTZHfAuUKUX+miWS4vAif81SH6vKewr+fRuvORle8wTHjL413Jj2\nn/Gjf0/aSohng4LBqAv5EBhvwfqvRSRNwKjrwGEIYW2wJSI3A9OBcSKyBLhAVW8QkW8A97EpbPiF\navvqZ8YllzH9oAPtLToibF6HptVxcEIQKPu2ll+9PsxxfD8F3nZaedvz1wx8M/Tjnu/M5yt+YAcJ\nN64mbDhIwstlmmN5WCtXIozy8gKWdvBOxwCrVdWPi803qpoDbg5Sp2Zxbar62SGu3wvcG0WfM877\nXhTNGkZZRremybXF8+fW0SuxbUa1MhtPOG8B2E5zLBR/c7he86zyEl6+ql2Mp9VXbrJJkmYSaR7T\nMHbyjHQ/lBM29iLyA6Au7PnJCZQ26h7zoQxNR0eafG/1iQkf7uxi/47SukP7BqdsivuwKKhWvNd7\nkHoFL7nmi1peKGwgj9C9MUdZf8LLddrOcB/bB4RGfE7543D3pqo5JlAMIwZGjGhBc9Vn5R1OllEj\nS7fTttIhG9PqC6WKFewBy/rtaxV93D5IwsvP0McwypsMQyOGsGEvPcqbqvpK5J35ILBAEZEJqrqi\n6HxiUvaXNx9KtNi8Dk3L8FY011Z1O4eOKN9GxnFi2ZYX+teHxLOCfQdy/IzypqihEl7u4yPh5Szv\nCIUqBYqIzGTzZRmCK1fPU9W7vWsnArdU1VGIVKKhtIrIl1T1ehGZCnwY+FXI44oE86EYtSI9tgPS\n2Vj6ak07ODGZvPLqf5+SgQSpJ4ST/v8M3xqKv2zNJSlh8pr17EIeeO7lktVV9YhS90UkhZsFPlRn\nfDVUEuW1RET+KSJnAI6qXhvBuIwEYj6UocmMHQEt1WsNDyxbxSGTx5fuy3FA4tFQtinkeT6AX6LY\nWf5Kwb+z3BFhnBZwfLz1jx8i4eV4Wn3VB8JZTl2iq+m7bM/0XTblXbz49r9V0sMRwAuq+nollaOg\nUh/KG8DWVJjHyzCaDRk/HoaFsFVuZx7ZcsuSRdrTaTLZeATKtprjWfw5n9d7zvLfD3CWr9V2RpQR\nSg7u89lPT6MkxXodxvH0shc55pFmC5+CK1ScyF3Ux1NH5i6ozIfSDpwOnAt8SESOVNWKxKvRWJh2\nUoLRo6C1+gfM9DHld0Rpa02Ty1e6e0gw1uA/b1gpZ3lHGVNUWmC1Or73rBolKUYxjAsL6/imU8EG\nYKFoKNFGeanqqZF2UAGVfMN3AH7kLXqZLSIfDHlMkWFOeaNmjJsAfd2xdNU+vJVCTCavbcnz0yp3\nh/TjLEfhZRyfutDmBKmzVHMs05AWNtqOjeVR1acHnD8U3nCixZzy0WI+lBKMGgPZ6rfKnfX0Aqbv\ntlPJMm1tGQq5MBbmled5SXOV+vu5qnGWtznCUdrr3wdSRJA60yTDNDI8nA9h/hKSvj5MbB2KYcSA\nbDkZ8tW/+crS1cjEgdsIbU5bRwua7626L1/jEXB8rgivxlkexIcyWN2akJz09aFhAsUIDdNOhkZG\nbwlavRnq0COPKlsm096KZnuq7ssPi5wMjk93TTXOckfgNdIVWZFqZnlKSPr6MDGBYhgxIMPLO9PD\nwhneDoV4NJQl6Qzpgv8AgDGSYoznLD8rgLM8hZvmPlkaSq06rh1NJVDMKR8t5kMpQSaEkGH8zbHT\n0QL56tO8+MFxHJwKAwAq8YdEzeJCjiWhOeWb6vEK+BQoItKhqp1efvyCagi6ew0wp7zRDKSHt4HG\nk7MqJVKxBhCknghIhQIoSL2tUxm2JsPcbAganpm83o2IfBcYLyIOcJl3fCXqgRnJw7ST6PEzxzJy\nBKTiSb3ipKRiH0WQev1O+Yr6qbBe1dShBhY1fjSUR4CHgSxwLDX8/RiGUR4ZNgx3J/UY+opJQ3Gk\n8udzzZzyFjY8KJ3AKar6c+A2z+xlGO/CfCjR42uORwx3l5bHgDhOLKYo01CSQVnhoKqPi8hhRZc2\nS5kqIpeq6rmhjywCzClvNAXDOmJ7LXdi0lDi4rVClkWFkJzyJlCG5ATgCu/z94Dbi+4dhZvXq+4x\np3y0mKCOHl9zPHwkZKrfHdIPbpRXhXWD+lAqfEAHiSbbNtXCtqkWHugNwylvAmUoZIjPg50bhlFL\n2jsgFU9m3bg0lESavCxseEh0iM+DnRtNivlQosefD2UkZKvfHdIProYSvebgaigVdVNDp3x0HYvI\nbsDPgDZcN8QZqvp4ZB36xK9A2U1E1uG+JLR7n/HO4/nmGobhC+kYDfl4dod0pLIMwFCfPpQwqdRE\n55MrgAtU9T4RORq4Ejg0yg794EugqGrMO9MYScS0k+jxtQ5l1HgIkA6lGiQVjw9FSKDJK9qw4QIw\nyvs8GlgeZWd+aT4jn2E0Ou0doSSi9EOjRXmFSrQC5T+Av4nID3Fl7Yei7MwvJlCM0DAfSvT4mWNp\n6SAu16aD4zt9/ZR8likF1xS3RFIcmHU3HFvqZFiaii4qze/4wu+4un5FZCYwofgS7i/2PODDwJmq\n+mcRORb4Fe4e8zWlqQSKrUMxmoLW+NyaLSkh5/NFfIWTYQWu4BjoPS6XyjJXhXwM8lx/OZ/llUJI\n/qcSGsqsec/wwLz5Jaur6pACQkRuUtUzvXJ/EJHrKx1mmIhqcwRpiYhq55paD8MwGoobd96P7pWr\nIu+nq1DZc+pb3e/ww/axget9u/sdVLViFUNENP/4X32XT+19VKD+ROQ53MiuB0TkcOD7qrpPBUMN\nlabSUAzDCJcWR8jHEJfbW9CKjXgNmnrlNOB/RCQF9FAnCXsb3i9mxMes2XNqPYSGp97mOC1CJoYj\nkYjj/wiIqj6oqnur6h6q+kFVfTKCnyAwpqEYhlExGYF8DM97RyoPM6jdFsAJFYRVYALFCA0Ldoie\nepvjFK6WEn0/gt9A6En5LJO8hZ3LJM3+fW402fJUhuURRpO9C0tfbxiG4R9HJJatfFPi3z6/Ip1h\nRdoVHE8MuBfrAy+pproqMIFihIatQ4meepvjdEwmrxZHkpc00Gm+BCMmUAzDqBgHIRXDwsFUFT6U\nmmEmr8bGFjZGi81r9NTbHKdESMVg2klLPBrKC7k+XsiFtLDRMYHS0NgGW4YRLo7EE8yUIh4N5QPp\nFj6QbuGO3q6q24o423Bd0nwi1IiMelsj0YjU2xynkXgOifcIhQjXodQrTaWhGIYRLhKThpLId/0m\n1FBMoBihUW/2/Uak3uY4JUIhhgenE5MPJVQaSPPwiwkUwzAqxiEeu3kio7xSzRc23Hwi1IiMerPv\nNyL1Nsf9CxsjP4j3CAUR/0eDYBqKYRgVkxJBY1opnzgNxUxehlE59Wbfb0TqbY4zjiAxeOWD5PKq\nGyIUtCKyK/AzoANYBJykqhsi69AnJlAMw6iYdNqBTPRv4k5f5F1EQKSC9pfAWao6V0ROAb4L/FeU\nHfqh+XQyIzLqzb7fiNTbHKfTQibtRH7E7UcJhWh9KDuo6lzv89+BT4cz6OowDcUwjIpJpx0kDg0l\niX7raH1Lz4rIMap6N3AcMDnKzvxiAsUIjXqz7zci9TbH6bSDk45DoAiqCXPLV+mUF5GZwITiS7ix\nCecBXwSuEZH/Au4C6sIo2FQCxZJDGka4pNNCIR1HlJfEslz+yb5enurrDaexEuOd9eAjzHro0ZLV\nVfWIMj0cCSAi2wMfDTi6SJDESf0KERHVzjW1HkZDU297dTQi9TbHKz/xUQqr34m8n1eWrSdfiO9Z\ndfBby1HVikWYiGhh+Yu+yzuTdgzUn4hsoaorRcQBbgDuV9Ubg480XJpKQzEMI1xScZq8Iu8lZKL1\noZwoIl/DNYH9qR6ECZiGYhhGFaw9/hh0zerI+3nx1bXk8/E9qw5Ysax6DeWNhb7LO1ttX1V/9YJp\nKIZhVIyTSaOZ6B8jycyKlXj5EBhbh2KERr2tkWhE6m2OJZ1CMjEcjiAS3xHO5FguL8MwDN9IOgXp\n6PUHQUicF8VyeRlG5dRT9FGjUndznHGQTBwCJYEGpAbSPPxiAsUwjIqRlAMxRHmlUgKFpD2gkzbe\n6mk+ncyIjHqz7zci9TbHsfhP+n0oMR6hzE3cPps6wDQUwzAqxtVQYjB5iSTPgpS4AVePCRQjNOrO\nvt+A1NscO3E55Z0kPp8TN+CqMYFiGEbFSGsGspnI+0mnHSTG1CuhkDwJWDUmUIzQqLc8U41Ivc2x\nZFIQx8JGR5KXw97Chg3DMALQ2ga56DOnpzIpJGlpokxDMYzKqac350al7uY4k4ZM9CYvJ5VEDaXW\nA4gfEyiGYVSMpDOxCJSU44QWzhsfSRtv9SROoIjINrg7lo1U1eO8ax/H3WBmBPArVZ1ZwyE2LfVm\n329E6m6OY9JQUilBUgl7QDehyStxXiNVfU1Vvzzg2p2q+hXgq7j7Kxs14Kln5td6CA1P3c1xv4YS\n8ZFKO6RjPEIhwuSQInKsiDwrInkR2XPAve+JyEIReUFEPhLOD+OPmmkoInI98DFgharuWnT9KOBH\nuMLuelW9PECz5wM/DnWghm/WrF1b6yE0PHU3xy0ZyLdE3o2TctCkmbyijfKaD3wS+PlmXYq8D/el\n+n3AZODvIrK9xrTxVS01lBvw9kTux9vO8lrv+s64u5Lt5N37vIhcJSJb9RcfUPf7wF9U9anIR16G\natJjBK1brnyp+0Pd83u9lmlAqu07SH0/ZcOa58TNcSoN6cGPWa+9MeS9oGUdR3hwfReOI4MeQ90b\n7Lqfa6EQoYaiqi+q6kLe7aj5OHCrquZUdRGwENi36p/FJzUTKKo6Fxi41du+wEJVXayqWeBW3AlC\nVW9S1bOAXhH5KbC7iJwNICLfAA4HjhWRr8T2QwzBrDlzY6tbrnyp+0Pd83t94PmixUtKjiVMqpnj\noPX9lA1rnhM3x5kMZFoGPR545fUh7wUtm0oLD23oJJWWQY+h7g123c+1cJAAR2hMApYWnS/3rsVC\nTbcAFpFpwN39Ji8R+TRwpOcPQUQ+B+yrqv8eQl8JC2I3DKOWVLkF8CJgWoAqK1R14oA2ZgITii/h\n7iF/nqre7ZW5H/iWqj7hnV8LPKiqN3vnvwTuUdU7Kv1ZglBvUV6D/QJDEQSNsF+zYRjJQFW3DqGN\nIyqotgyYUnQ+GXi92rH4pd6ivJYBU4vOY50MwzCMBFL8snwXcIKItHhLLLYDHo1rILUWKAMNiI8B\n24nINBFpAU7AnSDDMAzDQ0Q+ISJLgf2B/xWRewFU9XngNuB54C/AGXFFeEENfSgicjMwHRgHrAAu\nUNUbRORoNg8b/n5NBmgYhmEEoqZOecMwDKNxqLXJq6aIyDYi8ksRua3WY2lURGSYiNwoIj8Xkc/W\nejyNiH2Po0dEPi4i14nILSJSibO8KTANBRCR2/rzghnh4oV+r1bVe0TkVlU9odZjalTsexw9IjIa\nuFJVT6v1WOqRhtBQROR6EVkhIs8MuH6UiCwQkZf6F0Ea1VHBXE9m00KrfGwDTTD2fY6eKubY0juV\noCEECiGncTFKEmiucYXJ5P6icQ0y4QSd443F4hleQxB4juspvVO90hACJcw0LkZpgs41cAduSpwf\nA3fHN9LkEnSORWSsfY+DUcEc11V6p3ql3lbKh8nAnDbLGJAkTVXfwU15b1THkHOtql3AF2sxqAaj\n1Bzb9zgcSs3xNcA1tRhUkmgIDWUIIkvjYrwLm+vosTmOHpvjKmlkgWJpXOLD5jp6bI6jx+a4ShpJ\noFgal/iwuY4em+PosTkOmYYQKF4alweBHURkiYicqqp54BvAfcBzuJvOvFDLcTYCNtfRY3McPTbH\n0WALGw3DMIxQaAgNxTAMw6g9JlAMwzCMUDCBYhiGYYSCCRTDMAwjFEygGIZhGKFgAsUwDMMIBRMo\nhmEYRiiYQDFqiojkReQJEXnS+/+7tR5TPyJyu4hsLSIPe2NbLCJvFY116hD1LhaRiwZc20tEnvY+\n/11ERsTxMxhGnNjCRqOmiMg6VR0Zcpspb9VzNW28H7hYVT9ddO1kYC9V/Xcfde9Q1R2Lrl0JrFLV\ny0XkVGALVb2imjEaRr1hGopRawbdFEpEXhORGSIyT0SeFpEdvOvDvN32HvHuHeNdP1lE7hSR/wP+\nLi4/EZHnReQ+EblHRD4lIoeJyJ+K+vmwiPxxkCGcBNxZdvDuDn8Pisjj4u433q6qzwPdIrJHUdHP\n4O6vAW5+qM/6mRzDSBImUIxa0z7A5PWZontvqepewM+Ab3vXzgP+T1X3Aw4DfiAi7d69PYBPqeqh\nwKeAqar6fuDzwAcBVPUfwE4iMs6rcyrwq0HGdQAwr9TARWQL4BzgMFXdG5gPfNO7fStwolfuAGC5\nqi72xvA2MFxEQtXMDKPWNPIGW0Yy6FLVPYe4d4f3/zzgk97njwDHiMh3vPMWNqUcn6mqa73PBwK3\nA6jqChG5v6jdm4DPiciNwP64AmcgWwEry4z9Q8D7gQdFRIAMMNe7dwswC/gubtbaWwbUXeX1sa5M\nH4aRGEygGPVMr/d/nk3fVQE+raoLiwuKyP5AZ/GlEu3eiLsdcS9wu6oWBinTBbSVGZ8A96rqyQNv\nqOpiEXldRA7CFYYDhWYb0F2mfcNIFGbyMmpNqQf/YPwN2OgUF5Hdhyg3F/i050uZAEzvv6Gqb+Bu\nnHQernAZjBeA7cqM5UHgEBHZxhvLMBEprnMr8D/A86r6VtGYBRjH5tvNGkbiMYFi1Jq2AT6US73r\nQ4UfXgxkROQZEZkPXDREuT/i7sD3HPAbXLPZ2qL7vwOWquqCIer/BTi01MA9IfEl4Pci8hTwT2D7\noiK3ATvzbnPXvsBctRBLo8GwsGGjYRGRDlXtFJGxwCPAAf2agohcAzyhqjcMUbcN+IdXJ9Q/EhG5\nFvi9qs4Js13DqDXmQzEamf8VkdG4zvKLioTJ48AG4KyhKqpqj4hcAEzC1XTC5AkTJkYjYhqKYRiG\nEQrmQzEMwzBCwQSKYRiGEQomUAzDMIxQMIFiGIZhhIIJFMMwDCMUTKAYhmEYofD/AWrMSpO1P8/H\nAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Add the examples script directory to the path\n", "import sys\n", "sys.path.append(ct_dir + '/share/examples/python/')\n", "\n", "# Plot the SED and likelihood profile\n", "import show_spectrum as show_spectrum\n", "show_spectrum.plot_spectrum(sed['outfile'].filename(), '')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the above image we can see the spectral points that are produced by `csspec`. Additionally, behind the spectral points we can see the profiled likelihood values. These are plotted as the change in log-likelihood as a function of differential flux." ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.11" } }, "nbformat": 4, "nbformat_minor": 2 }