{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "
\n", "**This is a fixed-text formatted version of a Jupyter notebook.**\n", "\n", " You can contribute with your own notebooks in this\n", " [GitHub repository](https://github.com/gammapy/gammapy-extra/tree/master/notebooks).\n", "\n", "**Source files:**\n", "[spectrum_pipe.ipynb](../_static/notebooks/spectrum_pipe.ipynb) |\n", "[spectrum_pipe.py](../_static/notebooks/spectrum_pipe.py)\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Spectrum analysis with Gammapy (run pipeline)\n", "\n", "In this tutorial we will learn how to perform a 1d spectral analysis.\n", "\n", "We will use a \"pipeline\" or \"workflow\" class to run a standard analysis. If you're interested in implementation detail of the analysis in order to create a custom analysis class, you should read ([spectrum_analysis.ipynb](spectrum_analysis.ipynb)) that executes the analysis using lower-level classes and methods in Gammapy. \n", "\n", "In this tutorial we will use the folling Gammapy classes:\n", "\n", "- [gammapy.data.DataStore](http://docs.gammapy.org/dev/api/gammapy.data.DataStore.html) to load the data to \n", "- [gammapy.scripts.SpectrumAnalysisIACT](http://docs.gammapy.org/dev/api/gammapy.scripts.SpectrumAnalysisIACT.html) to run the analysis\n", "\n", "We use 4 Crab observations from H.E.S.S. for testing." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setup\n", "\n", "As usual, we'll start with some setup for the notebook, and import the functionality we need." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "from gammapy.data import DataStore\n", "from gammapy.scripts import SpectrumAnalysisIACT\n", "\n", "# Convenience classes to define analsys inputs\n", "# At some point we'll add a convenience layer to run the analysis starting from a plain text config file.\n", "from gammapy.utils.energy import EnergyBounds\n", "from gammapy.image import SkyImage\n", "from gammapy.spectrum import models\n", "from regions import CircleSkyRegion\n", "from astropy.coordinates import SkyCoord\n", "import astropy.units as u" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Select data\n", "\n", "First, we select and load some H.E.S.S. data (simulated events for now). In real life you would do something fancy here, or just use the list of observations someone send you (and hope they have done something fancy before). We'll just use the standard gammapy 4 crab runs." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Use observations [23523 23526 23559 23592]\n" ] } ], "source": [ "# TODO: Replace with public data release\n", "store_dir = '$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2'\n", "data_store = DataStore.from_dir(store_dir)\n", "obs_id = data_store.obs_table['OBS_ID'].data\n", "print(\"Use observations {}\".format(obs_id))\n", "\n", "obs_list = data_store.obs_list(obs_id)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Configure the analysis\n", "\n", "Now we'll define the input for the spectrum analysis. It will be done the python way, i.e. by creating a config dict containing python objects. We plan to add also the convenience to configure the analysis using a plain text config file." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "crab_pos = SkyCoord.from_name('crab')\n", "on_region = CircleSkyRegion(crab_pos, 0.15 * u.deg)\n", "\n", "model = models.LogParabola(\n", " alpha = 2.3,\n", " beta = 0,\n", " amplitude = 1e-11 * u.Unit('cm-2 s-1 TeV-1'),\n", " reference = 1 * u.TeV,\n", ")\n", "\n", "flux_point_binning = EnergyBounds.equal_log_spacing(0.7, 30, 5, u.TeV)\n", "\n", "exclusion_mask = SkyImage.read('$GAMMAPY_EXTRA/datasets/exclusion_masks/tevcat_exclusion.fits')" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "config = dict(\n", " outdir = None,\n", " background = dict(\n", " on_region=on_region,\n", " exclusion_mask=exclusion_mask,\n", " min_distance = 0.1 * u.rad,\n", " ),\n", " extraction = dict(containment_correction=False),\n", " fit = dict(\n", " model=model,\n", " stat='wstat',\n", " forward_folded=True,\n", " fit_range = flux_point_binning[[0, -1]]\n", " ),\n", " fp_binning=flux_point_binning\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Run the analysis\n", "\n", "TODO: Clean up the log (partly done, get rid of remaining useless warnings)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:gammapy.scripts.spectrum_pipe:Running SpectrumAnalysisIACT\n", "INFO:gammapy.background.reflected:Running ReflectedRegionsBackgroundEstimator\n", "Region: CircleSkyRegion\n", "center: \n", "radius: 0.15 deg\n", "ObservationList\n", "Number of observations: 4\n", "Info for OBS_ID = 23523\n", "- Start time: 53343.92\n", "- Pointing pos: RA 83.63 deg / Dec 21.51 deg\n", "- Observation duration: 1687.0 s\n", "- Dead-time fraction: 6.240 %\n", "Info for OBS_ID = 23526\n", "- Start time: 53343.95\n", "- Pointing pos: RA 83.63 deg / Dec 22.51 deg\n", "- Observation duration: 1683.0 s\n", "- Dead-time fraction: 6.555 %\n", "Info for OBS_ID = 23559\n", "- Start time: 53345.96\n", "- Pointing pos: RA 85.25 deg / Dec 22.01 deg\n", "- Observation duration: 1686.0 s\n", "- Dead-time fraction: 6.398 %\n", "Info for OBS_ID = 23592\n", "- Start time: 53347.91\n", "- Pointing pos: RA 82.01 deg / Dec 22.01 deg\n", "- Observation duration: 1686.0 s\n", "- Dead-time fraction: 6.212 %\n", "\n", "\n", "INFO:gammapy.background.reflected:Processing observation Info for OBS_ID = 23523\n", "- Start time: 53343.92\n", "- Pointing pos: RA 83.63 deg / Dec 21.51 deg\n", "- Observation duration: 1687.0 s\n", "- Dead-time fraction: 6.240 %\n", "\n", "INFO:gammapy.background.reflected:Found 7 reflected regions\n", "INFO:gammapy.background.reflected:Processing observation Info for OBS_ID = 23526\n", "- Start time: 53343.95\n", "- Pointing pos: RA 83.63 deg / Dec 22.51 deg\n", "- Observation duration: 1683.0 s\n", "- Dead-time fraction: 6.555 %\n", "\n", "INFO:gammapy.background.reflected:Found 7 reflected regions\n", "INFO:gammapy.background.reflected:Processing observation Info for OBS_ID = 23559\n", "- Start time: 53345.96\n", "- Pointing pos: RA 85.25 deg / Dec 22.01 deg\n", "- Observation duration: 1686.0 s\n", "- Dead-time fraction: 6.398 %\n", "\n", "INFO:gammapy.background.reflected:Found 19 reflected regions\n", "INFO:gammapy.background.reflected:Processing observation Info for OBS_ID = 23592\n", "- Start time: 53347.91\n", "- Pointing pos: RA 82.01 deg / Dec 22.01 deg\n", "- Observation duration: 1686.0 s\n", "- Dead-time fraction: 6.212 %\n", "\n", "INFO:gammapy.background.reflected:Found 17 reflected regions\n", "INFO:gammapy.spectrum.extract:Running \n", "INFO:gammapy.spectrum.extract:Process observation\n", " Info for OBS_ID = 23523\n", "- Start time: 53343.92\n", "- Pointing pos: RA 83.63 deg / Dec 21.51 deg\n", "- Observation duration: 1687.0 s\n", "- Dead-time fraction: 6.240 %\n", "\n", "INFO:gammapy.spectrum.extract:Update observation meta info\n", "INFO:gammapy.spectrum.extract:Offset : 0.5000156610786867 deg\n", "\n", "INFO:gammapy.spectrum.extract:Fill events\n", "INFO:gammapy.spectrum.extract:Extract IRFs\n", "/Users/deil/Library/Python/3.6/lib/python/site-packages/astropy/units/quantity.py:634: RuntimeWarning: invalid value encountered in true_divide\n", " result = super().__array_ufunc__(function, method, *arrays, **kwargs)\n", "INFO:gammapy.spectrum.extract:Process observation\n", " Info for OBS_ID = 23526\n", "- Start time: 53343.95\n", "- Pointing pos: RA 83.63 deg / Dec 22.51 deg\n", "- Observation duration: 1683.0 s\n", "- Dead-time fraction: 6.555 %\n", "\n", "INFO:gammapy.spectrum.extract:Update observation meta info\n", "INFO:gammapy.spectrum.extract:Offset : 0.4999843633857045 deg\n", "\n", "INFO:gammapy.spectrum.extract:Fill events\n", "INFO:gammapy.spectrum.extract:Extract IRFs\n", "INFO:gammapy.spectrum.extract:Process observation\n", " Info for OBS_ID = 23559\n", "- Start time: 53345.96\n", "- Pointing pos: RA 85.25 deg / Dec 22.01 deg\n", "- Observation duration: 1686.0 s\n", "- Dead-time fraction: 6.398 %\n", "\n", "INFO:gammapy.spectrum.extract:Update observation meta info\n", "INFO:gammapy.spectrum.extract:Offset : 1.501990838622691 deg\n", "\n", "INFO:gammapy.spectrum.extract:Fill events\n", "INFO:gammapy.spectrum.extract:Extract IRFs\n", "INFO:gammapy.spectrum.extract:Process observation\n", " Info for OBS_ID = 23592\n", "- Start time: 53347.91\n", "- Pointing pos: RA 82.01 deg / Dec 22.01 deg\n", "- Observation duration: 1686.0 s\n", "- Dead-time fraction: 6.212 %\n", "\n", "INFO:gammapy.spectrum.extract:Update observation meta info\n", "INFO:gammapy.spectrum.extract:Offset : 1.5017625673554083 deg\n", "\n", "INFO:gammapy.spectrum.extract:Fill events\n", "INFO:gammapy.spectrum.extract:Extract IRFs\n", "INFO:gammapy.spectrum.fit:Running SpectrumFit\n", "Source model LogParabola\n", "\n", "Parameters: \n", "\n", "\t name value error unit min max frozen\n", "\t--------- --------- ----- --------------- --- --- ------\n", "\tamplitude 1.000e-11 nan 1 / (cm2 s TeV) nan nan False\n", "\treference 1.000e+00 nan TeV nan nan True\n", "\t alpha 2.300e+00 nan nan nan False\n", "\t beta 0.000e+00 nan nan nan False\n", "Stat wstat\n", "Forward Folded True\n", "Fit range [ 0.7 30. ] TeV\n", "Backend sherpa\n", "Error Backend sherpa\n", "/Users/deil/code/gammapy/gammapy/stats/fit_statistics.py:161: RuntimeWarning: divide by zero encountered in log\n", " term2_ = - n_on * np.log(mu_sig + alpha * mu_bkg)\n", "/Users/deil/code/gammapy/gammapy/stats/fit_statistics.py:166: RuntimeWarning: divide by zero encountered in log\n", " term3_ = - n_off * np.log(mu_bkg)\n", "/Users/deil/code/gammapy/gammapy/stats/fit_statistics.py:203: RuntimeWarning: divide by zero encountered in log\n", " term1 = - n_on * (1 - np.log(n_on))\n", "/Users/deil/code/gammapy/gammapy/stats/fit_statistics.py:204: RuntimeWarning: divide by zero encountered in log\n", " term2 = - n_off * (1 - np.log(n_off))\n", "/Users/deil/code/gammapy/gammapy/stats/poisson.py:383: RuntimeWarning: divide by zero encountered in double_scalars\n", " temp = (alpha + 1) / (n_on + n_off)\n", "/Users/deil/code/gammapy/gammapy/stats/poisson.py:384: RuntimeWarning: divide by zero encountered in log\n", " l = n_on * log(n_on * temp / alpha)\n", "/Users/deil/code/gammapy/gammapy/stats/poisson.py:385: RuntimeWarning: divide by zero encountered in log\n", " m = n_off * log(n_off * temp)\n" ] } ], "source": [ "ana = SpectrumAnalysisIACT(\n", " observations=obs_list,\n", " config=config,\n", ")\n", "ana.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Check out the results\n", "\n", "TODO: Nice summary page with all results" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Fit result info \n", "--------------- \n", "Model: LogParabola\n", "\n", "Parameters: \n", "\n", "\t name value error unit min max frozen\n", "\t--------- --------- --------- --------------- --- --- ------\n", "\tamplitude 2.374e-11 1.807e-12 1 / (cm2 s TeV) nan nan False\n", "\treference 1.000e+00 0.000e+00 TeV nan nan True\n", "\t alpha 1.883e+00 1.770e-01 nan nan False\n", "\t beta 1.853e-01 7.857e-02 nan nan False\n", "\n", "Covariance: \n", "\n", "\tname/name amplitude alpha beta \n", "\t--------- --------- -------- ---------\n", "\tamplitude 3.26e-24 2.09e-13 -6.22e-14\n", "\t alpha 2.09e-13 0.0313 -0.0128\n", "\t beta -6.22e-14 -0.0128 0.00617 \n", "\n", "Statistic: 26.058 (wstat)\n", "Fit Range: [ 0.77426368 27.82559402] TeV\n", "\n" ] } ], "source": [ "print(ana.fit.result[0])" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(,\n", " )" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgIAAAGOCAYAAADl32vyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3Xl8nGd97/3Pb3Yto9EuWZslr1nI\n7iQGEiAhtPRAAoeHJS3wQArkwCnQAm1ZSstWDi20UJ6Ww1IoW3OgtE9ZAqRsWWghmw1k8W7Lli3b\nsmTt0oykWa7zx0jK2LEtaTSLpPm+X6/7Jc1o7vv6ybEz37nuazHnHCIiIlKaPMUuQERERIpHQUBE\nRKSEKQiIiIiUMAUBERGREqYgICIiUsIUBEREREqYgoCIiEgJUxAQEREpYQoCIiIiJcxX7AIKob6+\n3nV2dha7DBERkYLZuXPnaedcw0KvK4kg0NnZyY4dO4pdhoiISMGYWc9iXqdbAyIiIiVMQUBERKSE\nKQiIiIiUMAUBERGREqYgICIiUsIUBEREREqYgoCIiEgJUxAQEREpYQoCIiIiJawkVhZcSZxzOOdI\npVLzx9zjzJ/NfX+hY+56mdfNbOdc359dy4WY2aKen3t8oa+LOTwez9O+Zh4iIpJ7CgJZmJqaYmpq\nikQiQTKZJJlMkkqlnvb93Bt95vcLvfnK+Xk8Hrxe79O+Zh4+n2/+69nHOd33MbjpvYX9RUREVhAF\ngSwMDAzQ399f7DJKzlyYyoaZzQcCv98/f7Q+8FeMXPUW/H4/gUAAv9+f46pFRFa2NR0EzOxW4NZN\nmzYVuxTJsOWX7yh2CWfwfv02UsDU7OExD+axp756PGc8lxd3/CA/1xURWcCaDgLOubuBu7dt2/am\nYtciK0cg2kcwdmr+cXjwMQCmy5qYKW8m5VKQhCTJp51r2FPh4KzDOPeYChGRlWxNBwE5v2TKMRF3\nTM6kmJhJMRl3ROMpYnFHNOGIxVNMJRzTScf07NeZJMSTjnjKEU9C0jmSKUg4RyoFjtljdhyEYZiB\nWXp6itdj+Dzgtffj84DfawS8ht8DAa8R8j11BH0eyn1Gud8o93so9xuVAQ/hgFER8ODz5OZN95q7\nb2bnrffm5FrBYJCysrIzjmAweN5BlyIiK4GCwBoSTzmGYymGYkmGp1KMTCUZmUrNH2PTKcZm0l+j\n8YUHLQa8EPIaQZ8RnHvT9hp+jxEKgtc8eD3gNcPrASP9pj/3tpcOBemvKQeJVDo4JJ0jkUqHismZ\nFDMpmJkLHAnHVHLh2sp8RjjoIRL0UDX7tTqUPmpCXmrKPNSWeagNefF7C/NGPD09zfT0NCMjI/PP\neTweysrKKC8vp6ysjIqKCsrKyhQORGTFUBBYRWLxFKcmk5yaTDIwmWQgmj5OR5MMRlOMTqc4+y3U\nA1SF0m+U/8P9Kz+o+V0iQQ/hgIfKQPpTdkXAQ4U//em7wu+hzJ/+VO4p0ptVyqV7IGJxRzTumIyn\ng8vEbO/FxIxjfCbF+HT6dx6KJTk8HGd0OsW5MkRV0ENdmYf6ci8Nc0eFl8YKL3UbX5vf3yWVYnJy\nksnJyfnnzGw+FMwdoVAor3WIiJyPgsAKE0ukOD6W5OREgpPjs18nkpyaSDA2c+a7XNBr1Jd7aCj3\n0lntp67MQ12Zl9qy9CfimpCHcNCDd/YN/Zq7/5Wt299SjF9rSTxmlPmMMh/Uli3+vJRzjM84hmd7\nRIZjSQZjKQZj6bB0aiLJk/0zxBKZf46/Q+X+UzRVemmu9LGu0su6sI+WSi+tYR8VgdwPDnTOEY1G\niUajDAwMAOD1eqmoqKCyspLKykoqKiq0doKIFISCQJ6db4T8aCrIoXg9BxP1HIrX052o53CillPJ\nqjNet847SodvmCt8I7SFRmjzjdDmHaXVO0LEM8X8h/aZ2WM0u3oKaf+zPpWX63rMiASNSNBD53le\n45xjMu4YiCbpn0yHg77JBKcmkhwcivPgsSkyJyhWBz20Vvloq/LSVuWjo8pHR8RPVTC3b9LJZJKx\nsTHGxsaAp3oNwuEwlZWVhMNhvF5vTtsUEQEFgbxLOqMnUcu+eCP74g3sizeyP95Ifyo8/5pym2aj\nb5DrAz10+Yfo9A2x3jdEu2+EkCWW1f5CI+RLjZnN3xLpqn76mgHxlKN/MsmJ8QTHx5McH0vQO5bg\nP49OnTGuojrooSPio7PaR2e1n85qH61hX84GMWb2Gpw6lf7vV1ZWRlVV1Xw4UDAQkVywUljpbtu2\nbW7Hjh05u96xY8fOuaBQ0jl6xxIcGorTPZKgezjOkZEE07M3rr0GbVU+1kfSR3vET3uVj4ZyT0EG\nj+VihHzm6n2ZK/zNLQuceZzt7CWUM5dZnluVce5YaZxzDE+lODqa4OhYgmOjCXpG4xwdTRCf7ULw\neWB9xMfGGj9dNX421vjpiPjw5ygcZDIzKioqCIfDVFVVUVFRoQGIInIGM9vpnNu20OvUI7AMp6NJ\n9g3OcGAozsGhON3DT73ph3xGV7WPWzaU0VXtp6vaR2tVft4UcsHn8xEIBM44Mlfgm3vzL9SbTTKZ\nJJFInHHE4/EzjpmZGeLxeEGWbTYzamfHX1zZHHyqzpTj+HiCntF08Ds8kuAXx6b4cXcMAL8Huqr9\nbKpNH1vr/DRVLP/P0TnHxMQEExMTnDx5Eq/XS1VVFVVVVUQiEa2QKCKLpiCwRL85NsKn/qOHJ06M\nMxRLfxSc+5/987vK5v+Hv67SW7RR9+dzcsvr5ue3h0IhgsHg/NeV1s08t3dAMBhc8LVzoWBmZobp\n6en5r3Pf5zMoeD1GR8RPR8TPjR3pkY3OOU5NJjk0nA6IB4fi/OxwjB8ejAIQCXrYUpcOBRfXB9hY\n41/2FMdkMsnw8DDDw8NA+jZCJBIhEolQWVm5vF9SRNY0BYElis0k2T8Q45L6AFvr/GypC7C+euV9\n0vd6vZSXl59xBK/+NOtWWDjJhblei4qKiqf9zDk3HwzmNouaO+LxeF7qMTOaK300V/p4dns6HCRT\njmNjCfYPxtk3OMO+wTiPnphO1++BTbV+LqoPcGlDgIvq/JT5lzcYMRaLEYvF6Ovrw+fzEYlEqK6u\npqqqSrMRROQMGiOwRM45ent7V9SmQ2fPS6+srFzUJ+lSl0wmmZqaIhaLzX+NRqMkEssboLlYo9Mp\n9p2eYc/pGfacjtM9HCfpwGOwscbPpQ0BntEY4OJ6PyFfbt68PR4P4XCYmpoaIpHI+XdlFJFVT2ME\n8mQlDMiaGyg2N62soqJixXXtrwZzc/fP7kmIx+Pzn6jnRu5PT0/n/BZDJOjhutYQ17WmFxOaSqTY\nNxhn18AMu/pn+P7+Sb6zbxKfwZY6P5c1BrmiOcCmGj/eLHugUqkUo6OjjI6OpmdQVFZSU1NDdXW1\nxhWIlCj1CGThfLMG8mluTvnc9DF17xZWKpUiFosxOTk5Hw6mpqbyOv5gKpFiz+k4T/TP8GT/NN3D\nCRxQ7jcuawxwZXOQK5uCNFbkJgTOhYKamhqFApE1QD0Cq5zH45kfAa5R4MXn8Xie1nuQuXxwNBpl\nYmIip+MOQj4PVzUHuao5CIQZn07xRP8Mvzk1zW/6pnn4eHqMQVuVj6ubA1zVHOTihkDW41XmZiH0\n9vZSWVlJbW0t1dXVun0gssapRyAL+eoR8Pv9VFdXU11dTTgcXhG3IWRpZmZmmJycnH9TjcVieek1\ncM5xfDzJr/qm+fXJaXafniGRSm/GdGVzgG3rQly1LkhkmSsgmhmRSITa2loikYh6okRWEfUIrBKB\nQGD+Hq2mea1+c2sw1NTUAE/1GoyPjzMxMcHk5CSpVGqBqyzMzGir8tFW5eO2LRXEEime7J9hx4lp\ndp6c5sHeaQy4qN7PdS0hrm0Nsq5y6f/cnXOMjIwwMjKC1+ulurqauro6wuHwwieLyKqgHoEsLLdH\nwO/3U1NTQ21t7TmnvMna5Zyb7zGYCwe5CAaZUs5xeCTBjhNTPHpimsMj6VkQ7VU+rm8Nsr0tRGfE\nt6wep0AgQG1tLXV1ddo5UWSFWmyPgIJAFrIJAh6PZ/7Nv6qqauETpCTMBYPx8fH5YJDrf5P9kwke\nOTHNI8en2DMQJwU0V3i5vi3Es9pCbKxZXiioqKigrq6O2tpazV4RWUEUBDIUMwiEw2Hq6+uprq7W\n/VVZUCqVYmJiYn4nwlgsltPrj04leeTENA/1TvFE/wxJB43lXp7VHuLZ7SG6qrMPBR6Ph+rqaurr\n63XrQGQFUBDIUOggEAgEqKuro66uTgv7yLLE4/H5UDA2NpbTxY4mZlI8cnyKX/ZO8fipdChorvRy\nY0eIG9rLaKt6+piCdfu+wsmtr1/w2sFgkPr6eurq6jTjRaRIFAQyFCoIRCIRGhoaqKqq0oh/yYvJ\nycn5BYGi0WjOrjs+neLh41P817EpnuyfwQFd1T5u7Cjjxo4QtWXpLv+l7mA5N+ugvr6eSCSSs3pF\nZGGaNVAgPp+P+vp6GhoaCAQCxS5HVqsvv2hRL6uYPVqAlEuRTCRJJBM56Sm4BviffuhvruRHsYv4\n4eTFfO3xFv758VGuD/Zwa/kurgG2/PIdWV1/wmx+XwiP5fE22R0/yN+1RdYgBYEslZeX09DQQG1t\nre79S1F4zIPH78Hv9+NwZ2zdvJyevkbvBK+t3MFrK3dwOF7L6MQEN6YehdnhCuHBxwCYCjURr2he\n9HXnNoCamZnB5/Olt7f26n9BIsWmf4VZWLdunVZbk9xa5qdYI/2Pee5v5eTkJCMjIwwPDzM9Pb2s\na5cDO5xj72Cc1zz4Qi5OfINYwtHo9XJTpIybuspoKM9utkAoFKKhoYG6ujrNOBApEo0REFnjYrHY\nfChY7iyEa+6+mV/+zs945MQU9x6O8UT/DACXNwW4pauca1uC+L1LHx/j8Xioq6ujsbFR6xKI5IgG\nC2ZQEBBJm5qaYnh4OOtQcPasgf7JJPcdiXLvkRinoymqAsbzOst4wYZyWsLZ9ZqFw2EaGxuprq7O\n6nwRSVMQAMzsVuDWTZs2venAgQPFLkdkRZmammJoaIjh4WGmpqaWda2kczx+aoafdkd59MQ0SQfP\naAjwgg1lXN8WymojpGAwSENDA/X19bptIJIFBYEM6hEQubBoNDofCmZmZpZ1reGpJPcdjvGT7hj9\n0SSRoIdbNqR7CbIZS+DxeKivr6exsVHrcogsgYJABgUBkcUbHx9ncHCQkZERkslk1tdJOcdv+mb4\n0aEoO0+mN0G6piXI72wq5/LGQFZrbUQiEZqamrRyocgiaB0BEclKOBwmHA6TSqUYHR3l9OnTjI+P\nL3lKoseMq9cFuXpdkP7JBD/pjs3fOmgLe/mdTRU8tzNEmW/x02/nFlMqLy+nsbGR2tpaLd4lskzq\nERCRBcXjcYaGhjh9+vSyxhPMJB2/ODbFPQcnOTScoNxnPH9DGf9tUzmNFUv/XOL3+2lsbKShoUHj\nCETOolsDGRQERHJncnKS06dPMzQ0lPUWys459g/F+cGBKA/2ToGD61qDvHhLBRfV+Zf8Kd/r9VJX\nV0dTU5NW+BSZpSCQQUFAJPdSqdR8L8Hk5GTW1xmMJrnnUJSfHIoyEXdsrvVz65ZytreG8C5xtoGZ\nUVNTQ3NzM2VlZVnXJLIWKAhkUBAQya9oNDrfS5DtAMOpRIr7e6a4e/8kfRNJGsu9vGhLOc/vKlvS\nOII5VVVVNDc3a2ChlCwFgQwKAiKFkUqlGBwcZGBgIOtVDJPOsePENN/bN8newTiVfuOFm8r5b5sr\niASXHgjKy8tpbm6mpqYmq3pEVisFgQwKAiKFNzExQX9/PyMjI1lvgrR/cIZv753kkRPTBDxwc1c5\nt20tpymLgYWhUIimpibq6uo000BKgoJABgUBkeKJx+MMDAwwMDCQ9XbJvWMJvrtvkp/3xEgBN3aE\neNlFlbRVLT0QBAKB+ZkG2jlU1jIFgQwKAiLF55xjaGiIU6dOZX3bYDCa5Hv7J/lJd4yZpOO61iAv\nv7iSDTX+JV/L5/PR2NhIY2Ojph7KmqQgkEFBQGRlGRsbo7+/n9HR0azOH51O8cMDk/zwYJRo3HHN\nuiAvv7iCLXVLnzro9XppaGigqalJ24vLmqIgkEFBQGRlmpqaoq+vj6GhoazGEUzGU9xzMMrd+yeZ\nmHFc0RTglZdUclH90gPB3J4GWotA1goFgQwKAiIrWzwe59SpU5w+fTqr6YexeIofHYry3f1RxqZT\nXNEU4FWXVrI1ix4CM6Ouro7m5mZtciSrmoJABgUBkdUhmUwyMDDAqVOnshpYOJVI8aNDMb6zb3I+\nENx+aWVWtwzMjNraWtatW6dAIKuSgkAGBQGR1SWVSnH69GlOnTqV1bbIU4kU/3Eoxnf3TjA2kx5D\ncPul2Q0qnFutcN26dYRCoSWfL1IsCgIZFAREVifnHKdPn6avry+rQBBLpLjnQJTv7JtkMu7Y3hrk\n9meEac9i2iEw30OgQCCrgYJABgUBkdXNOcfg4CAnT57MKhBMzqS4+8Ak398fZTrheM76EK+6tDKr\nHQ+B+R4C7WcgK5mCQAYFAZG1Ya6H4OTJk8Tj8SWfPzad4tt7J7jnYBTn4Lc2lvPyiyuIhLJbR0CB\nQFYyBYEMCgIia8vcGIK+vr6sAsFgNMm3dk9w75EYAa/xki0V3Lq1PKvNjSAdCFpaWnTLQFYUBYEM\nCgIia1MqlaK/v5++vr6sph2eGE/wf54c58HeaSJBD6+8pJJbNpThW+L2x3M0hkBWEgWBDAoCImtb\nIpGgr6+PgYEBUqnUks/fPzjD1x8fZ/fpOOsqvbz28jDXtQSz2pxI0w5lpVAQyKAgIFIaZmZmOHHi\nBIODg0s+1znHr/qm+drjE/SOJbiozs/rrghntQYBPLUw0bp167RSoRSFgkAGBQGR0hKLxejt7WVs\nbGzJ5yZTjp8difEvT04wMp3i2e0hXnNZ9jMMzIyGhgaam5vx+5e+joFIthQEMigIiJSmsbExent7\ns9rtMBZP8Z19k3xv/yTOwYu3VPCyiyoo92c3oNDj8cwHAm1uJIWgIJBBQUCkdM1NOTxx4kRWyxYP\nRpPc9eQ4D/RMEQl6+N1nVHJzVxneLMYPQHq3w6amJm1/LHmnIJBBQUBEkskkJ06cYGBgIKudDg8O\nxfnyb8bYOxinq9rHHVdWcWlD9vf+fT4fzc3NNDQ04PFk18sgciEKAhkUBERkztTUFMeOHctq/IBz\njl/2TvG1x8Y5HUvxrLYQr708TGNF9p/s/X4/LS0t1NXVZTVLQeR8FAQyKAiIyNmGh4fp7e3Nasni\n6YTju/sm+fa+CXDw3y+u5CVbKwh6s38jDwaDtLa2UlNTk/U1RDIpCGRQEBCRc0mlUvT19dHX15fV\n7YKBaJKvPjbOg71TNJZ7uePKMNdmuf7AnPLyclpbW6mqqsr6GiKgIHAGBQERuZCpqSmOHj3K+Ph4\nVuc/0T/NF389Tu9YgiubArzhqipawsubGRAOh2lra6O8vHxZ15HSpSCQQUFARBZjaGiI3t7erPYv\nSKQc/3Ewyjd3TRBPOV66tYKXXVRJ0Le8+/41NTW0trZqlUJZMgWBDAoCIrJYyWSS3t5eTp8+ndX5\nw7EkX3t8nJ8fTd8u+P2rwlzbsry9B8yM+vp6WlpatAaBLJqCQAYFARFZqvHxcY4ePcrU1FRW5+8a\nmOEffzXGsbEE17YEecNVVTSUL2/dgLk1CJqamjTlUBakIJBBQUBEsuGc48SJE5w6dSqrwYSJlOP7\nB6J8a9cEAK+8tJIXby7PenfDOZpyKIuhIJBBQUBEliMajdLT00M0Gs3q/P7JJP/0mzEePTFNR8TH\nm6+pYmuWmxllKisro7W1lUgksuxrydqjIJBBQUBElss5R19fHydPnsyqdwDgkeNTfPHXYwzFUvz2\nxnJ+77JKKrLcuyCTZhjIuSgIZFAQEJFcicViHDlyJOvegVg8xTd2TXDPgSiRkIc3XFXF9tblrT0w\np66ujpaWFm17LECOg4CZ1QAtQAw44pxLLb/EwlEQEJFccs5x8uTJrBcigvTeBZ/bOcrhkQTXtQR5\n49VV1JUtfxMij8dDU1MTzc3NGlBY4pYdBMwsAvwB8LtAABgAQkAT8BDwv51z9+Ws4jwws1uBWzdt\n2vSmAwcOFLscEVljotEoR44cyWqbY4Dk7GDCb+4ax2fGay8Pc8uGMjw56B2YG1BYX1+/7GvJ6pSL\nIPAT4GvA3c65kbN+dg3wWuAJ59yXclBvXqlHQETyJZVKcfz4cfr7+7O+Rt9Egs/tHOOJ/hkuqffz\nlm2RZa9MOKe8vJy2tjbC4XBOrierh8YIZFAQEJF8Gx8f58iRI1ltYgTp2w33Honx1cfGiScdr3pG\nmFs3l+Nd5lTDOdXV1bS1tWmFwhKSsyBgZh7gCp4aI7DLOXcqJ1UWiIKAiBRCMpmkp6eH4eHhrK8x\nHEvy+V+lpxpuqvHzB9dW0RHx56Q+M5sfP+D1Ln88gqxsubg1sBF4N3ALcICnxghsAaLA54GvroaB\ngwoCIlJIg4ODHD16lFTq/P97XLfvK5zc+vpz/sw5xy97p/jir8aIxh2vuKSSl15UseyFiOZo/EBp\nyEUQ+AbwWeA/3VkvMrNG4PeAYefcV3NQb14pCIhIoU1PT9Pd3X3eaYbX3H0zO2+994LXGJ1O8aVf\nj/GLY1NsqPbx1usirM9R7wCkxw+0t7dTWVmZs2vKyqExAhkUBERKyJdfVOwK5jkcM9MzzMSfPm4g\nPPgY43VXLOo6P41t4SMjv8VYKsRbqn7BHZUP47fcdcb6fD6CwSAeO2u64R0/yFkbUniLDQILTjI1\ns1eYWXj2+/eb2b+b2dW5KFJEZC0zjGAwSFlZGWZGINpHePAxwoOPAcx/H4j2XfA6t5Tt5zuNX+IF\nZfv4+7Hn8P8OvJrueF3O6kwkEkxOTjIzM4Nj7X84lDMtZrDg4865y83sBuBjwN8A73POXV+IAnNB\nPQIiUmzxeJzu7m4mJtIbEC3m1sC5PNg7xed3jjKVcPzeM8K8aEs53hxuPBQMBmlvb9f+BWtAznoE\ngOTs1xcBn3XOfZf0AkMiIrJIfr+fLVu2sG7dumVd55ltIf7ut+u5sjnIVx8f5wP3D9E3kchRlemx\nDQcPHuTgwYNMT0/n7Lqyci0mCBw3s88DrwR+aGbBRZ4nIiIZzIyWlhY2bdpE30Wvz/o61SEv735W\nNW+7NkLPSIJ3/XiQn3ZHs17u+FxGR0fZtWsXJ06cuODsB1n9FnNroBx4IelVBA+Y2TrgMufcjwtR\nYC7o1oCIrDQzMzN0d3czOTm5rOsMRJP8wyOjPDkww7Z1Qd6yrYrqUG7XCAgEArS3t1NdXZ3T60p+\nadZABgUBEVmJnHP09vYua3ligJRz/PBAlH9+Ypwyn/GWbRGuaw3lqMqnRCIR2tvbtTrhKpHLMQIi\nIpIHZkZ7ezudnZ3L2inQY8aLt1Tw8VvqqCv38te/HOGzO0aZSuS2S1+3C9YmBQERkSKrq6tj69at\ny/6k3RHx87Hn1/HSrRX87HCMP/7JIAeGstv74HzmtmDevXs3Y2NjOb22FIeCgIjIClBeXs5FF11E\nVVXVsq7j96S3M/7Q82qJpxzvu3eIf9szQTLHt4Gnp6c5cOAA3d3dxOPxnF5bCuu8QcDM2s3sm2b2\nn2b2PjPzZ/zsO4UpT0SkdPh8PjZv3kxzc/Oyr3VpQ4BPvqCeZ7WF+MaTE3zw/iFOR5MLn7hEw8PD\n7Nq1i/7+/pzOWliS+z5WnHbXiAv1CPwTcD/wNmAd8ICZzS1ltT7PdYmIlKzW1la6urqWNW4AoCLg\n4Y+uj/C26yJ0jyR4549P82DvVI6qfEoymeTYsWPs3bt32bMgsvLAXxW+zTXEd4GfNTjnPjf7/dvM\n7DXAz83sNtAalCIi+VRbW0soFOLQoUPMzGR/n9/MeN76Mi6q8/N3D4/yNw+O8PyuMn7/yjAh34WD\nxpZfvmPJ7aWAKb+fYDCIkbsVDxe0gvaYOKcVvG/DhYKA38xCzrkpAOfcP5tZH/AjoKIg1YmIlLDy\n8nIuvvhiDh06NL80cbaaK3385U21/MuuCb69d5K9p2d45/ZqOqtzt5vhnHg8TiKRIBgM4vfl/voA\njPTA6LGnHvf8V/prpB2q1Wm9FBfahvgdwK+ccw+c9fxVwMedcy8oQH05oXUERGQ1c87R09PD4OBg\nTq73+KlpPv3IKJMzKV5/RRW/vTG9KVI+RCIROjo6CATyuDL9ByPwwdH8XX+VWvY6As65T50dAmaf\n//VqCgEiIqudmdHZ2UlbW1tOrnd5U5BPvqCOZzQG+Mdfj/GJB0eYmMnPugBzaw8sd9EkyZ/FbEPc\nZWafnN1++HtzRyGKExGRpzQ1NbFp06ZlDyIEiIS8vO+GGl53eZgdJ6b5k58Msn8wt2sOzEmlUvOD\nCWOxWO4beO57cn/NErKYvQYeA74EPEF6HAgA5+otWKl0a0BE1pJoNMrBgwdzNn9//+AMn3xohKFY\nitdcFubFW8rx5OlWgZnR1NRES0tL3m5HSFrO9hows4edc9fnrLIiUBAQkbUmHo9z8OBBotFoTq43\nMZPif+8Y5eHj01yzLsjbro0QDuZvzblQKMT69euprKzMWxulLpdB4PeAzcCPgfnNqZ1zv1pukYWi\nICAia1EqlaK7u5vR0dwMlHPOcc/BKF99bJzqMg/v2l7Nlro8DvIDGhoaaG1txevN7Y6Jktsg8DHg\ntcAhnro14JxzNy+7ygJREBCRtezo0aMMDAzk7HoHh+L87YMjDMaSvPbyMC/eXJ7XbvxAIMD69euX\nvbyynCmXQWAvcLlzLj+jSApAQUBE1rpTp07R29ubs+tNzKT4h0dHefTENNe3BvmDayNU+PO7PU1d\nXR3t7e3qHciRXG5D/BhQvfySREQkX5qamtiwYUPOPrlXBjy8+1nVvO7yMI+emOZPfzrIkZH8bi40\nODjIrl27GBkZyWs7cqbFBIEmYK+Z/UjTB0VEVq6amhq2bNmSs0/UZsZtWyv48PNqmU443vuzQe49\nkpvBiecTj8c5dOgQ3d3dJBKJvLYlaYu5NfDccz2v6YMiIitTLBbj4MGDy9qj4GwjU0n+7uFRnuif\n4ZauMt5wVRUBb36n//l8Pjo2MHZLAAAgAElEQVQ6OqipqclrO2tVLm8NHAUeds49MPvm/wjQs9wC\nRUQkP8rKyti6dSuhUChn16wOefnz59Twsosq+OnhGO+/b5D+ydxva5wpkUjQ3d2t3oE8W0wQ+Fcy\nFhICkrPPiYjIChUIBNi6dWtO5+l7zXj1ZWHe8+xqTo4n+dOfnuaxU9MLn7hMw8PDGjuQR4sJAr7M\nGQOz3+d3YqmIiCybz+dj8+bNRCKRnF732pYQf31LHTUhLx/5+TD//54JFrrNvFyJRIJDhw5x+PBh\n9Q7k2GKCwICZ3Tb3wMxeApzOX0kiIpIrHo+HjRs3Ultbm9PrtoR9fOz5tdzQEeL/PDnBJx4cIRrP\nz8ZFmYaGhti9e3fOFlGSxQWBNwPvM7OjZnYUeDdwZ37LEhGRXDEzurq6aGxszOl1Qz4Pf3hdhDuu\nSE8xfM/PBjk+nv9P63PLK/f09JBM5necQilYcNbA/AvNKmdfP57fknJPswZERNJOnjzJiRMncn7d\nJ/un+duHRoknHW+/LsJ1rbkbqHghgUCAzs5OwuFwQdpbTZY9a8DMXmNm8z93zk1khgAz22hmNyy/\nVBERKZR169bR3t6e8+s+ozHIJ26poyXs5a9/OcK3dk2QyvO4AYCZmRn2799Pb29v3scprFW+C/ys\nDvi1me0EdgIDQAjYBDyX9DgBbQItIrLKNDY24vF46OnJ7Uzw+nIvf3lTHZ/fOcq/7J7g8Eict18X\noSzPSxNDeonlsbExOjs7KS8vz3t7a8kFbw2YmRe4GXg2sA6IAXuAe5xzRwtSYQ7o1oCIyNMNDw9z\n+PDhnH+Sds7xg9ldDFvCPt7z7GrWVV7oc2fumBktLS00NzcXpL2VLGebDq0FCgIiIuc2OjpKd3c3\nqVTuR/w/0T/N3z44QsrBO7dXc2VzMOdtnE9lZSVdXV0EAqU72z2XKwuKiMgaFYlE2LRpEx5P7t8O\nLmsM8te31FFf7uWj/znM3fsnC3Yff2Jigt27dzM4OFiQ9lYzBQERkRIXDofZvHlzXrb/barw8dGb\na7m2NchXHhvnMzvGiCcLEwaSySRHjhyhu7tb0wwv4EKzBp5pudrPUkREVrTKysq8hYEyn4c/fmY1\nr7ykgvuOxPjAA0OMTBXujXl4eJjdu3czPr7qZr8XxIV6BF4H7DSzb5rZ681MIy9ERNawiooKtmzZ\ngs+X+4F9HjNedWmYP35mNUdGErz7p4McHonnvJ3zmZtmePz4cU0zPMt5g4Bz7s3OuauBDwI1wFfM\n7EEz+19m9pzZGQUiIrKGlJeXs3nz5ryEAYBntoX4y5tqccCf3TvEw8en8tLO+fT19bF3716mpgrb\n7kq24BgB59xe59ynnHMvJD2V8L+AVwAP57s4EREpvPLy8rz1DABsqPHz18+voyPi4+O/HOHfC7Bp\nUaZoNMqePXs4fVrb5sASBws652LOuR865962mCkJIiKyOpWVleU1DNSUefnQ82q5sSPEXU9O8PeP\njhZsECFAKpWip6dHAwnRrAERETmPfIeBoNf4w+si3H5pJQ/0TPGhnw8xNp3/HQwzzQ0knJiYKGi7\nK4mCgIiInFe+w4CZ8YpLKnnn9giHhuK8+2eDHBvL/w6GmeYGEp48ebKg7a4UCgIiInJBZWVleZta\nOOfZ7WV86Hm1zCQc7/vZII+dms5bW+finOPEiRPs37+feLxwsxlWggWDgJmNm9nYWccxM/u2mW0o\nRJEiIlJccwMI8xkGttQF+Njzn1qJ8Kfd0by1dT7j4+Ps3r2b0dHRgrddLIvpEfgk8CdAK9AG/DHw\nj8A3gX/KX2kiIrKSzE0tzGcYaKzw8tGba7msMcBnd47x9cfHC7KdcaZEIsHBgwdLZmvjxQSBFzrn\nPu+cG3fOjTnnvgD8N+fcv5BeX0BEREpERUVF3vYmmFPu9/C+G2r4rQ1lfGffJJ98aITpAs4omHPq\n1Cn27dvHzMxMwdsupMX8l0yZ2SvNzDN7vDLjZ2s/KomIyBkqKyvzHga8HuPOq6t43eVhHuqd5kMP\nDDFa4BkFAJOTk+zevZuRkZGCt10oi/mv+GrgtUA/cGr2+9eYWRnw1jzWJiIiK1Q4HGbDhg3kc0sa\nM+O2rRW865nVHB6O876fDXJivLAzCiC9edGhQ4fW7K2CCwaB2WWEX+Kcu9U5V++ca5j9/uDs4kL/\nVaA6RURkhYlEInkPA5BelviDz6slGk/x3nsH2Xu6OF31a/VWwQWDgHMuCbykQLWIiMgqU11dTWdn\nZ97b2To7oyAc8PDBB4b45bHi7BUwOTnJnj171tSsgsXcGviFmf2Dmd1oZlfPHXmvTEREVoXa2lo6\nOjry3k5zpY+P3VzHhho/n3xohO/vn8x7m+cyN6tgrexkuJilop41+/XDGc850hsQiYiI0NDQQDKZ\n5Pjx43ltJxz08IHn1vLph0f48mPjDESTvO6KMJ483544l76+PiYnJ+nq6sLv9xe8/VxZMAg4524q\nRCEiIrK6NTc3k0wm6evry2s7Qa/xrmdW85XfjPP9A1EGYynefl2EgLfwYWB8fJw9e/bQ1dVFOBwu\nePu5sJiVBZvM7Etmds/s40vM7A35L01ERFab1tZWGhoa8t6O14zfvzLM664I82DvFB/5+RCTM4Wf\nXggQj8c5cOBA3gNQvixmjMBXgB8BLbOP9wN/lK+CRERkdevo6KCmJv/rzZkZt22p4I+uj7B/MM77\n7xtiMFqcLYWdcxw/fpxDhw6tum2NFxME6p1z3wJSAM65BLC6fksRESmorq4uqqqqCtLWjR1l/NmN\nNQxEk7z33kGOjhZv06CRkRH27NlDLBYrWg1LtZggMGlmdcyuImhm24G1M29CRERyzszYuHEjFRUV\nBWnv8qYgH7mplqSD9983xJ4irTUAMD09zd69exkaGipaDUuxmCDwTuB7wEYz+wXwNeBtea0qg5lt\nmB2j8G8Xek5ERFYWj8fDpk2bCIVCBWmvq9rPx26uJRL08OEHhnj0RHHWGgBIpVIcPnyYY8eOrfgp\nhgsGAefcr4Dnkp5G+D+AS51zjy/m4mb2T2bWb2ZPnvX8C81sn5kdNLP3LNB+t3PuDQs9JyIiK4/P\n52Pz5s0Fm17XWOHjozfX0RHx8/FfjPDTw4XfyjhTf38/+/fvJx4v3u2KhZw3CJjZDXPfO+cSzrld\nzrknnXPx2Z9XmdkzFrj+V4AXnnVdL/AZ4HeAS4DfnZ2JcJmZff+sozHL30tERFaIQCCQ9+2LM1UF\nPXzweTVc1hTgszvG+Pc9E0X9VD4xMcGePXuYnCzOAkgLudA6Av+PmX0c+A9gJzAAhIBNwE3AeuBd\nF7q4c+7nZtZ51tPXAQedc90AZvZN0vsZfAx4cRa/wzmZ2Z3AnUBBVrwSEZHzKysrY9OmTezfv78g\nb8plPg/vvaGGzzw6yl1PTjAyneL1RVp4CNJTDPft20d7e3tBplcuxXl7BJxz7wBeBJwEXgF8hPR4\ngc3A551zz3HOPZpFm63AsYzHvbPPnZOZ1ZnZ54CrzOy953vuHPV/wTm3zTm3baX9oYuIlKLKykq6\nuroK1p7fY7z9uggv2lzODw5E+ftHRkmkitcz4Jzj6NGj9PT0rKhxAxdcWdA5Nwz84+yRK+eKY+f9\nE3HODQJvXug5ERFZ+Wpqamhvb+fYsWMLvzgHPGbccUWYqqCHbzw5wWTc8a7t1QR9xekZADh9+jSx\nWIyNGzeuiKWJFzNrINd6gfaMx23AiSLUISIiRdDY2EhTU1PB2jMzXn5xJXdeXcWvTk7zkf8s3iqE\nc+Z2MZyYmChqHVCcIPAosNnMuswsANxOenqiiIiUiLa2toKsPpjptzeW847tEQ4MxvmL+4cYmSru\n2njxeJz9+/dz+vTpotaR1yBgZt8AHgS2mlmvmb1hdmXCt5JetngP8C3n3K581iEiIitPV1cXlZWV\nBW3z2e1lvOeGGk5OJHn/fUMMFGlJ4jnOOXp6eoraM7CYTYfKzezPzewfZx9vNrNFje53zv2uc26d\nc87vnGtzzn1p9vkfOue2OOc2Ouc+urxfQUREViMzK+iCQ3Ouag7yF8+pYXQ6xZ/dO0jvWKKg7Z9L\nMQcPLqZH4MvANPDM2ce9wF/mrSIRESkZXq+XzZs34/NdcOx6zl1UH+Ajz6slkYI/v2+Q7uGVu+BP\nvi0mCGx0zn0ciAM452Kce+S/iIjIkgUCATZt2oTHU9hha53Vfv7yploCPuMD9w+xt4j7ExTTYv7U\nZ8ysjKc2HdpIuodAREQkJyoqKgq6xsCclrCPj95UR3XIw4d/Psxjp0rv7W0xQeADpFcXbDezu4Cf\nAX+a16pyxMxuNbMvjI5qs0QRkZWuurqatra2grdbX+7lIzfV0lzp5X/91zAPHy/eZkXFsJhNh34C\nvAx4PfANYJtz7v78lpUbzrm7nXN3RiKRYpciIiKL0NTUVJQleKtDXj78vFq6qv38zYMj/LwnVvAa\niuVCmw5dPXeQ3lfgJOmFfzpmnxMREcm59vZ2qqqqCt5uZcDDB55TwyX1Af6/R0b5SXdxdy4slAsN\n0/zb2a8hYBvwGOlBgpcDDwM3nOc8ERGRrJkZGzZsYO/evUxNFbabvszv4X031vCJXw7zuZ1jTCcd\nL95cUdAaCu1Cmw7d5Jy7CegBrp7dwOca4CrgYKEKFBGR0uP1etm0aVPBpxUCBL3Gu59Vw/WtQb78\nm3H+fU/xlwHOp8UMFrzIOffE3APn3JPAlfkrSUREBILBIBs3bsSKsHWw32u8c3s1N7SHuOvJCb7x\n5PiK2jEwlxYTtfaY2ReBfyY9hfA1pJcGFhERyavKyko6Ojro6ekpeNs+j/H26yMEvMa/7ZkknoLX\nXlZZlGCST4sJAncAbwH+cPbxz4HP5q0iERGRDPX19cRiMfr7+wvetteMt2yrwu+F7+6bJJFy3HFF\neE2FgQWDgHNuCvjU7CEiIlJwbW1tTE1NMTY2VvC2PWa86aoq/B7j+weixJOON11dhWeNhIEFg4CZ\nPRv4IOkphPOvd85tyF9ZIiIiTynmTIK59l9/RRi/1/j23kkSKXjztiq8ayAMLObWwJeAdwA7geLu\n1ygiIiXL6/WyceNG9u7dSzJZ+LcjM+PVz6jE54F/3T1Jyjn+57WRVR8GFhMERp1z9+S9EhERkQWE\nQiE2bNjAwYMHizKK38y4/dIwXjO+uWuCpIO3XRvB61m9YWAxQeA+M/sE8O9kbDbknPtV3qrKETO7\nFbh106ZNxS5FRERypKqqitbWVnp7e4tWwysuqcRrcNeTE6Qc/OF1qzcMLCYIXD/7dVvGcw64Offl\n5JZz7m7g7m3btr2p2LWIiEjuNDU1EY1GGRoaKloNL7u4Eo/H+Prj46Sc44+ur8a3CsPAYmYN3FSI\nQkRERJZi/fr1TE1NEY0Wb0+Al26twGvwlcfGcW6Ed2xffWHgvEHAzN55oROdc5/MfTkiIiKL4/F4\n2LhxI3v27CGRSBStjlu3pPci+Mpj43zqodUXBi60xHB49thGekGh1tnjzcAl+S9NRETkwgKBABs2\nbCj6Aj+3bqngjivCPHR8mk89NEIitXqWI77QpkMfcs59CKgnvenQu5xz7wKuAdoKVaCIiMiFhMNh\nWltbi10GL84IA59cRWFgMZsOdQAzGY9ngM68VCMiIpKFpqYmamtri11GOgxcGebhVdQzsJhZA18H\nHjGzb5OeLfDfga/mtSoREZElWr9+PbFYjFgsVtQ6Xry5AufSYwb+7uFR3nH9yp5auJhZAx81s3uA\nG2efusM59+v8liUiIrI0mYMHi7HyYKZbt6TDwFcfH8djK3udgcX0CMwtHrTiFxASEZHSFgwG6erq\n4uDBg8Uuhdu2VpByjq8/MYHH4G3XrczliBcVBERERFaLSCTCunXrOHnyZLFL4aUXVZJy6RUIvQZ/\ncG1kxe1aqCAgIiJrTktLC5OTk0XZtvhsL7u4koSDf9k1gddjvPmalbWF8WJmDYiIiKw6XV1dBAKB\nYpcBwCsvqeTlF1fws8Mx/vFXY0XZMOl81nQQMLNbzewLo6OjxS5FREQKzOfzrYjFhubcfmklL91a\nwY+7Y/zTb8ZXTBhY00HAOXe3c+7OSCRS7FJERKQIKioqaG9vL3YZQHoL49dcVsmLN5fzw4NRvv7E\nxIoIAxojICIia1pDQwMTExNF3alwjpnx+ivCJFLw3X2TBLxw+6XhotakICAiImve+vXriUajTE1N\nFbsUzIw3XBVmJun4192TBDzGe7YUr541fWtAREQE0osNbdiwAY9nZbztecx487YqbuwIcdeTE3z9\nkd7i1VK0lkVERAqorKxsxYwXAPCa8bZrI2xvDfL39x/h1Fhxeit0a0BEREpGfX09ExMTDA4OFrsU\nALwe44+2V+OraaWpKlSUGtQjICIiJaWjo4NQqDhvuufi9xhbmyqL1r6CgIiIlJSVNl6g2PSnICIi\nJWeljRcoJgUBEREpSfX19dTW1ha7jKJTEBARkZLV0dFBMBgsdhlFpSAgIiIly+v1rqj9CIpBQUBE\nREpaeXk5ra2txS6jaBQERESk5DU1NVGqG9St6SCgbYhFRGSxOjs78fv9xS6j4NZ0ENA2xCIislg+\nn4/Ozs5il1FwazoIiIiILEVVVRVNTU3FLqOgFAREREQytLa2Ul5eXuwyCkZBQEREJIOZ0dXVVTJL\nEJfGbykiIrIEoVCoZJYgVhAQERE5h/r6eqqrq4tdRt4pCIiIiJzH+vXr1/yUQgUBERGR8yiFKYUK\nAiIiIhdQVVVFY2NjscvIGwUBERGRBbS2thIKhYpdRl4oCIiIiCzA4/HQ1dW1JncpVBAQERFZhPLy\nclpaWopdRs4pCIiIiCxSU1MTlZWVxS4jpxQEREREFsnM6OzsXFOrDq6d30RERKQAgsEgbW1txS4j\nZxQERERElqihoYGqqqpil5ETCgIiIiJZ6OzsxOv1FruMZVvTQcDMbjWzL4yOjha7FBERWWP8fj8d\nHR3FLmPZ1nQQcM7d7Zy7MxKJFLsUERFZg2pra1f9xkRrOgiIiIjk2/r16/H5fMUuI2sKAiIiIsvg\n8/lW9S0CBQEREZFlqqmpoaampthlZEVBQEREJAc6OjpW5S0CBQEREZEcWK23CBQEREREcmQ13iJQ\nEBAREcmh1XaLQEFAREQkh3w+H+3t7cUuY9EUBERERHJsNS00pCAgIiKSBx0dHatiLwIFARERkTzw\n+/2rYrtiBQEREZE8qa+vX/HbFSsIiIiI5NH69evxeFbu2+3KrUxERGQNCAQCtLS0FLuM81IQEBER\nybOmpiYqKiqKXcY5KQiIiIgUwPr16zGzYpfxNAoCIiIiBVBWVkZTU1Oxy3gaBQEREZECWbduHcFg\nsNhlnEFBQEREpEA8Hg/r168vdhlnUBAQEREpoHA4TH19fbHLmLemg4CZ3WpmXxgdHS12KSIiIvNa\nW1tXzA6FazoIOOfuds7dGYlEil2KiIjIPJ/Pt2KWH17TQUBERGSlqqurIxwOF7sMBQEREZFi6ejo\nKPraAgoCIiIiRRIKhVi3bl1Ra1gZIxVERERKVHNzM/F4vGjtq0dARESkiMyMQCBQtPYVBEREREqY\ngoCIiEgJUxAQEREpYQoCIiIiJUxBQEREpIQpCIiIiJQwBQEREZESpiAgIiJSwhQERERESpiCgIiI\nSAkz51yxa8g7MxsAejKeigCjy7jkcs7P5tylnFMPnF7i9UvVcv8eFFox681327m8fi6uVch/40t9\nvf6NL47+fcN651zDgq9yzpXcAXyhWOdnc+5SzgF2FPvPd7Ucy/17UEr15rvtXF4/F9cq5L/xLF6v\nf+MF+ntQKvWW6q2Bu4t4fjbnLrdeObfV9udazHrz3XYur5+LaxXy3/hq+3u4Wqy2P9ei1VsStwZK\niZntcM5tK3YdIpIf+jcuuVaqPQJr2ReKXYCI5JX+jUtOqUdARESkhKlHQEREpIQpCIiIiJQwX7EL\nKIT6+nrX2dlZ7DJEREQKZufOnafdItYRKIkg0NnZyY4dO4pdhoiISMGYWc/Cr9KtARERWcV29gzz\nmfsOsrNnuNilrFol0SMgIiJrz86eYV79xYeYSaQI+Dzc9cbtXLO+pthlrTrqERARkVXpoe5BZhIp\nUg7iiRQPdQ8Wu6RVSUFARERWpe0b6gj4PHgN/D4P2zfUFbukVUm3BkREZFW6Zn0Nd71xOw91D7J9\nQ51uC2RJQUBERFata9bXKAAsk24NiIiIlDAFARERkRKmICAiIlLCFARERERKmIKAiIhICVMQEBER\nKWEKAiIiIiVMQUBERKSEKQiIiIiUsBUXBMzshWa2z8wOmtl7LvC6l5uZM7NthaxPRERkLVlRQcDM\nvMBngN8BLgF+18wuOcfrwsDbgYcLW6GIiMjasqKCAHAdcNA51+2cmwG+CbzkHK/7CPBxYKqQxYmI\niKw1eQsCZvZsM6uY/f41ZvZJM1u/wGmtwLGMx72zz2Ve9yqg3Tn3/QXav9PMdpjZjoGBgSx+AxER\nkbUvnz0CnwWiZnYF8KdAD/C1Bc6xczzn5n9o5gE+Bbxrocadc19wzm1zzm1raGhYfNUiIiIlJJ9B\nIOGcc6S79j/tnPs0EF7gnF6gPeNxG3Ai43EYeAZwv5kdAbYD39OAQRERkez48njtcTN7L/Aa4Dmz\nAwH9C5zzKLDZzLqA48DtwO/N/dA5NwrUzz02s/uBP3bO7chx7SIiIiUhnz0CrwKmgTc45/pI3+v/\nxIVOcM4lgLcCPwL2AN9yzu0ysw+b2W15rFVERKQkWbr3fm3btm2b27FDnQYiIlI6zGync27BW+c5\nvzVgZuNkDPDL/BHgnHNVuW5TREREspPzIOCcW2hAoIiIiKwQ+RwsCICZNQKhucfOuaP5blNEREQW\nJ58LCt1mZgeAw8ADwBHgnny1JyIiIkuXz1kDHyE9z3+/c64LeD7wizy2JyIiIkuUzyAQd84NAh4z\n8zjn7gOuzGN7IiIiskT5HCMwYmaVwM+Bu8ysH0jksT0RERFZonz2CLwEiAHvAP4DOATcmsf2RERE\nZIny1iPgnJvMePjVfLUjIiIi2ctbEDhrYaEA6X0GJrWgkIiIyMqRzx6BMxYWMrOXAtflqz0RERFZ\nunyOETiDc+47wM2Fak9EREQWls9bAy/LeOgBtnHuPQhERESkSPI5fTBzhkCC9MqCL8ljeyIiIrJE\n+RwjcEc255nZC4FPA17gi865vzrr5+8E3kg6XAwAv++c61lmuSIiIiUpH9sQ/z0XuAXgnHv7Bc71\nAp8BXgD0Ao+a2fecc7szXvZrYJtzLmpmbwE+DrwqJ8WLiIiUmHwMFtwB7CS94+DVwIHZ40ogucC5\n1wEHnXPdzrkZ4JucdTvBOXefcy46+/AhoC2HtYuIiJSUnPcIOOe+CmBmrwducs7FZx9/DvjxAqe3\nAscyHvcC11/g9W/gPDsamtmdwJ0AHR0diyldRESk5ORz+mALkLmWQOXscxdi53junLcZzOw1pGci\nfOJcP3fOfcE5t805t62hoWER5YqIiJSefM4a+Cvg12Z23+zj5wIfXOCcXqA943EbcOLsF5nZLcCf\nAc91zk0vv1QREZHSlM9ZA182s3t4qmv/Pc65vgVOexTYbGZdwHHgduD3Ml9gZlcBnwde6Jzrz3HZ\nIiIiJSXntwbM7KLZr1eTvhVwbPZomX3uvJxzCeCtwI+APcC3nHO7zOzDZnbb7Ms+Qfo2w7+a2W/M\n7Hu5/h1ERERKRT56BN5JepDe357jZ44Flhl2zv0Q+OFZz/1Fxve35KBGERERIT+zBu6c/XpTrq8t\nIiIiuZW3WQNm9gozC89+/34z+/fZ+/siIiKyQuRz+uCfO+fGzewG4LeBrwKfy2N7IiIiskT5DAJz\nqwi+CPisc+67QCCP7YmIiMgS5TMIHDezzwOvBH5oZsE8tyciIiJLlM835leSngb4QufcCFAL/Eke\n2xMREZElylsQmN0YqB+4YfapBOnNh0RERGSFyOesgQ8A7wbeO/uUH/jnfLUnIiIiS5fPWwP/HbgN\nmARwzp3gzE2IREREpMjyGQRmnHOO2d0Dzawij22JiIhIFvIZBL41O2ug2szeBPwU+GIe2xNZkp09\nw3zmvoPs7BkudikiIkWTz90H/8bMXgCMAVuBv3DO/SRf7Yksxc6eYV79xYeYSaQI+Dzc9cbtXLO+\npthliYgUXN6CAMDsG/9PAMzMa2avds7dlc82RRbjoe5BZhIpUg7iiRQPdQ8qCIhIScrHNsRVZvZe\nM/sHM/stS3sr0E16bQGRotu+oY6Az4PXwO/zsH1DXbFLEhEpinz0CHwdGAYeBN5IehGhAPAS59xv\n8tCeyJJds76Gu964nYe6B9m+oU69ASJSsvIRBDY45y4DMLMvAqeBDufc+GJONrMXAp8GvMAXnXN/\nddbPg8DXgGuAQeBVzrkjuStfSsU162sUAJZgZ8+wgpPIGpSPIBCf+8Y5lzSzw0sIAV7gM8ALgF7g\nUTP7nnNud8bL3gAMO+c2mdntwF8Dr8pd+SJyNg2uFFm7LD3VP4cXNEsyu4gQYEAZEJ393jnnqi5w\n7jOBDzrnfnv28XtJn/SxjNf8aPY1D5qZD+gDGtwFfpFwOOyuueaa5f1iIiVspOV6RtpvAPOAS1J9\n7BdUn3i42GWJyAU88MADO51z2xZ6Xc57BJxz3mWc3gocy3jcC1x/vtc45xJmNgrUkb4FMc/M7gTu\nBAgGg8soSURCY8ewVBLncVgqRWjs2MIniciqkNfpg1mwczx39if9xbwG59wXgC8AbNu2zd1///3L\nLm410/1dWa4z/w69tNjliMgCzM71dvl0Ky0I9ALtGY/bgBPneU3v7K2BCDBUmPJWJ93flVzQ4EqR\ntSmfSwxn41Fgs5l1mVkAuB343lmv+R7wutnvXw7ce6HxAXLuxXNERERghfUIzN7zfyvwI9LTB//J\nObfLzD4M7HDOfQ/4ElsYe2UAAAYbSURBVPB1MztIuifg9uJVvDrMLZ4TT6S0eI6IiJwh57MGVqJt\n27a5HTt2FLuMotIYARGR0mJmxZk1ICuT7u+KiMi5rLQxAiIiIlJACgIiIgW2s2eYz9x3kJ09w8Uu\nRUS3BkRECknTeWWlUY+AiEgBaTqvrDQKAiIiBTQ3nddraDqvrAi6NSAiUkDXrK/hrjdu13ReWTEU\nBET+b3t3GyJVFcdx/PvLhwqFChNfqD1BlIqR5QNlgRWFUWiooFFUolGBvTMw6k34JgiCorIMHyJj\nNU0jpbLMyKIoM83UkCQMF5E0Q1DSSP+9uFe8DK67O3N3xpnz+8Dgzjn3nPPfcc7O/557516zOvPX\nee18ksQFhSQdBP4oFF0CHKmhy1raV9O2O20up+JOjNahWt8H9dbIeHt67DL7L6Oves7x7m7vOd41\nnt9wZUQM7HSriEjuASxsVPtq2nanDdmlmBv+GjfDo9b3QUrx9vTYZfZfRl/1nONVbO85Xqf3QSrx\npnqy4NoGtq+mba3x2tk12+vayHh7euwy+y+jr3rO8WZ7HzaLZntdGxZvEocGUiLpx+jCtaXNrDl5\njlvZUl0RaGULGx2AmfUoz3ErlVcEzMzMEuYVATMzs4Q5ETAzM0uYEwEzM7OEOREwMzNLmBOBFiep\nn6R3JL0t6aFGx2Nm5ZF0jaRFklY1OhZrXk4EmpCkxZL+lLSjonyipN2S9kialxdPAVZFxOPApLoH\na2bd0p35HRG/R8SsxkRqrcKJQHNaCkwsFkjqBbwO3AsMBx6UNBwYAuzLNztZxxjNrDpL6fr8NquZ\nE4EmFBGbgMMVxWOBPfkewr/AcmAy0E6WDID/v83Oe92c32Y18wdD6xjMmT1/yBKAwcBqYKqkBTTf\ntbfNLHPW+S1pgKQ3gVGSnm1MaNbsejc6ACuNzlIWEXEMmFnvYMysVB3N77+AJ+sdjLUWrwi0jnZg\naOH5EGB/g2Ixs3J5fluPcSLQOjYD10q6WlJfYAbwUYNjMrNyeH5bj3Ei0IQktQHfAddJapc0KyL+\nA+YA64FfgfcjYmcj4zSz7vP8tnrz3QfNzMwS5hUBMzOzhDkRMDMzS5gTATMzs4Q5ETAzM0uYEwEz\nM7OEOREwMzNLmBMBsxYn6aSkbYXHvM5b9TxJeyX9Imm0pDV5bHskHSnEemsHbWdLereibFB++94+\nklZIOizpgfr8NmbNy9cRMGtxko5GRP+S++ydX+Smlj72AqMj4lChbAIwNyLu76TtZcBvwJCIOJ6X\nzQFGRsQT+fNlwKqI+LCWOM1anVcEzBKV75G/IOmnfM/8+ry8n6TFkjZL2ippcl7+mKSVktYCn0m6\nQNIbknZKWifpY0nTJN0laU1hnLslra4hzjGSvpK0RdInkgZFxN/At8B9hU1nAG3VjmOWKicCZq3v\n4opDA9MLdYci4iZgATA3L3sO2BgRY4A7gJck9cvrbgEejYg7gSnAVcBIYHZeB7ARGCZpYP58JrCk\nmsAlXQi8AkyNiJuBZcD8vLqN7MMfSUPzWDZVM45ZynwbYrPW909E3NhB3ek99S1kH+wA9wCTJJ1O\nDC4Crsh//jwiDuc/3wasjIhTwAFJX0J2b9z8+P3DkpaQJQiPVBn7MGAEsEESQC+yO/FBdtOdVyX1\nB6aTXX//VJXjmCXLiYBZ2k7k/57kzN8Dke2B7y5uKGkccKxYdI5+lwBrgeNkyUK15xMI2B4Rt1dW\nRMQxSRuAyWQrA09VOYZZ0nxowMwqrQeeVr4LLmlUB9t9A0zNzxUYBEw4XRER+4H9wPPA0hpi2QUM\nljQ2j6WvpBGF+jbgGeDSiNhcwzhmyXIiYNb6Ks8ReLGT7ecDfYDtknZw5ph8pQ/Ilul3AG8B3wNH\nCvXvAfsiYle1gUfECWAa8LKkn4GtwLjCJp+SHbZYXu0YZqnz1wfNrGqS+kfEUUkDgB+A8RFxIK97\nDdgaEYs6aLuXiq8Plhybvz5o1gVeETCzWqyTtA34GphfSAK2ADeQneXfkYPAF5JGlx2UpBXAeLJz\nFMzsHLwiYGZmljCvCJiZmSXMiYCZmVnCnAiYmZklzImAmZlZwpwImJmZJex/p86PimXvt2cAAAAA\nSUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ana.spectrum_result.plot(\n", " energy_range=ana.fit.fit_range,\n", " energy_power=2,\n", " flux_unit='erg-1 cm-2 s-1',\n", " fig_kwargs=dict(figsize = (8,8)),\n", ")" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Exercises\n", "\n", "Rerun the analysis, changing some aspects of the analysis as you like:\n", "\n", "* only use one or two observations\n", "* a different spectral model\n", "* different config options for the spectral analysis\n", "* different energy binning for the spectral point computation\n", "\n", "Observe how the measured spectrum changes." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.3" }, "nbsphinx": { "orphan": true } }, "nbformat": 4, "nbformat_minor": 2 }