{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "**This is a fixed-text formatted version of a Jupyter notebook**\n", "\n", "- Try online [![Binder](https://mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.11?urlpath=lab/tree/cta_sensitivity.ipynb)\n", "- You can contribute with your own notebooks in this\n", "[GitHub repository](https://github.com/gammapy/gammapy/tree/master/tutorials).\n", "- **Source files:**\n", "[cta_sensitivity.ipynb](../_static/notebooks/cta_sensitivity.ipynb) |\n", "[cta_sensitivity.py](../_static/notebooks/cta_sensitivity.py)\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Estimation of the CTA point source sensitivity" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "This notebook explains how to estimate the CTA sensitivity for a point-like IRF at a fixed zenith angle and fixed offset using the full containement IRFs distributed for the CTA 1DC. The significativity is computed for a 1D analysis (On-OFF regions) and the LiMa formula. \n", "\n", "We use here an approximate approach with an energy dependent integration radius to take into account the variation of the PSF. We will first determine the 1D IRFs including a containment correction. \n", "\n", "We will be using the following Gammapy class:\n", "\n", "* [gammapy.spectrum.SensitivityEstimator](..\/api/gammapy.spectrum.SensitivityEstimator.rst)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setup\n", "As usual, we'll start with some setup ..." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import astropy.units as u\n", "from astropy.coordinates import Angle\n", "from gammapy.irf import load_cta_irfs\n", "from gammapy.spectrum import SensitivityEstimator, CountsSpectrum" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Define analysis region and energy binning\n", "\n", "Here we assume a source at 0.7 degree from pointing position. We perform a simple energy independent extraction for now with a radius of 0.1 degree." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "offset = Angle(\"0.5 deg\")\n", "\n", "energy_reco = np.logspace(-1.8, 1.5, 20) * u.TeV\n", "energy_true = np.logspace(-2, 2, 100) * u.TeV" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Load IRFs\n", "\n", "We extract the 1D IRFs from the full 3D IRFs provided by CTA. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "filename = (\n", " \"$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits\"\n", ")\n", "irfs = load_cta_irfs(filename)\n", "arf = irfs[\"aeff\"].to_effective_area_table(offset, energy=energy_true)\n", "rmf = irfs[\"edisp\"].to_energy_dispersion(\n", " offset, e_true=energy_true, e_reco=energy_reco\n", ")\n", "psf = irfs[\"psf\"].to_energy_dependent_table_psf(theta=offset)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Determine energy dependent integration radius\n", "\n", "Here we will determine an integration radius that varies with the energy to ensure a constant fraction of flux enclosure (e.g. 68%). We then apply the fraction to the effective area table.\n", "\n", "By doing so we implicitly assume that energy dispersion has a neglible effect. This should be valid for large enough energy reco bins as long as the bias in the energy estimation is close to zero." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "containment = 0.68\n", "energy = np.sqrt(energy_reco[1:] * energy_reco[:-1])\n", "on_radii = psf.containment_radius(energy=energy, fraction=containment)\n", "solid_angles = 2 * np.pi * (1 - np.cos(on_radii)) * u.sr" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "arf.data.data *= containment" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Estimate background \n", "\n", "We now provide a workaround to estimate the background from the tabulated IRF in the energy bins we consider. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "bkg_data = irfs[\"bkg\"].evaluate_integrate(\n", " fov_lon=0 * u.deg, fov_lat=offset, energy_reco=energy_reco\n", ")\n", "bkg = CountsSpectrum(\n", " energy_reco[:-1], energy_reco[1:], data=(bkg_data * solid_angles)\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute sensitivity\n", "\n", "We impose a minimal number of expected signal counts of 5 per bin and a minimal significance of 3 per bin. We assume an alpha of 0.2 (ratio between ON and OFF area).\n", "We then run the sensitivity estimator." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "sensitivity_estimator = SensitivityEstimator(\n", " arf=arf, rmf=rmf, bkg=bkg, livetime=\"5h\", gamma_min=5, sigma=3, alpha=0.2\n", ")\n", "sensitivity_table = sensitivity_estimator.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Results\n", "\n", "The results are given as an Astropy table. A column criterion allows to distinguish bins where the significance is limited by the signal statistical significance from bins where the sensitivity is limited by the number of signal counts.\n", "This is visible in the plot below." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Table length=19\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
energye2dndeexcessbackgroundcriterion
TeVerg / (cm2 s)
float64float64float64float64str12
0.01935723.77783e-11142.5531826.68significance
0.02887531.06474e-11220.6074420.97significance
0.04307357.74203e-12180.2282938.04significance
0.06425324.22421e-12111.9851118.13significance
0.09584711.73696e-1275.5441499.549significance
0.1429761.23576e-1250.3005215.216significance
0.2132798.1623e-1331.625280.8063significance
0.318155.57432e-1320.714932.1657significance
0.4745874.1604e-1313.996313.1734significance
0.7079463.25326e-139.720055.43152significance
1.056052.64284e-137.197562.4464significance
1.575322.2918e-136.002561.44036significance
2.349921.90236e-135.14260.878003significance
3.505392.09701e-1350.522997gamma
5.229032.62358e-1350.337574gamma
7.800193.51431e-1350.161236gamma
11.63565.19097e-1350.0682205gamma
17.3576.78151e-1350.0365443gamma
25.89151.05727e-1250.0155002gamma
" ], "text/plain": [ "\n", " energy e2dnde excess background criterion \n", " TeV erg / (cm2 s) \n", " float64 float64 float64 float64 str12 \n", "--------- ------------- ------- ---------- ------------\n", "0.0193572 3.77783e-11 142.553 1826.68 significance\n", "0.0288753 1.06474e-11 220.607 4420.97 significance\n", "0.0430735 7.74203e-12 180.228 2938.04 significance\n", "0.0642532 4.22421e-12 111.985 1118.13 significance\n", "0.0958471 1.73696e-12 75.5441 499.549 significance\n", " 0.142976 1.23576e-12 50.3005 215.216 significance\n", " 0.213279 8.1623e-13 31.6252 80.8063 significance\n", " 0.31815 5.57432e-13 20.7149 32.1657 significance\n", " 0.474587 4.1604e-13 13.9963 13.1734 significance\n", " 0.707946 3.25326e-13 9.72005 5.43152 significance\n", " 1.05605 2.64284e-13 7.19756 2.4464 significance\n", " 1.57532 2.2918e-13 6.00256 1.44036 significance\n", " 2.34992 1.90236e-13 5.1426 0.878003 significance\n", " 3.50539 2.09701e-13 5 0.522997 gamma\n", " 5.22903 2.62358e-13 5 0.337574 gamma\n", " 7.80019 3.51431e-13 5 0.161236 gamma\n", " 11.6356 5.19097e-13 5 0.0682205 gamma\n", " 17.357 6.78151e-13 5 0.0365443 gamma\n", " 25.8915 1.05727e-12 5 0.0155002 gamma" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Show the results table\n", "sensitivity_table" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# Save it to file (could use e.g. format of CSV or ECSV or FITS)\n", "# sensitivity_table.write('sensitivity.ecsv', format='ascii.ecsv')" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEOCAYAAACjJpHCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VNX5x/HPE0ASAVFRIoiyiBYXDIRNxQVcwA03pFahlYpStajVat3tr2KrtbWuSEVpXX8WQUtF0apVWVwqiwgqLv1RVARcUJF9fX5/nAmGmEwmmeXO8n2/XvfFzJ07dx5ykzw595zzHHN3RERE6qso6gBERCS3KZGIiEhSlEhERCQpSiQiIpIUJRIREUmKEomIiCRFiURERJKiRCIiIklRIhERkaQokYiISFIaRh1AJuy0007erl27qMMQEckps2bN+tLdd67tuIJIJO3atWPmzJlRhyEiklPM7KNEjtOtLRERSYoSiYiIJEWJREREklIQfSQikv02bNjAokWLWLt2bdShFJzi4mLatGlDo0aN6vV+JZLq7LILfPbZ9/eXlsLSpZmPR6QALFq0iGbNmtGuXTvMLOpwCoa7s2zZMhYtWkT79u3rdY68vrVlZgPMbMzy5cvr9sbqkki8/SKStLVr19KiRQslkQwzM1q0aJFUSzCvE4m7T3L34c2bN486FBFJgJJINJL9uud1IhERSdbZZ5/Nu+++W6/3Ll68mFNPPXXL89NPP53999+fW2+9leuuu44XXnghVWFGSn0kIpJ7MtiPed9999X7va1bt2bChAkALF26lFdffZWPPkpojl9OUYtERHJPmvoxV61axXHHHUdZWRn77bcf48aNo0+fPlsqY4wdO5a99tqLPn36cM455zBixAgAhg4dyoUXXshBBx1Ehw4dtiSPhQsXst9++wHQr18/Pv/8c7p06cK0adMYOnToluNmzJjBQQcdRFlZGT179mTFihUsXLiQQw45hPLycsrLy3n11VcBePnll+nTpw+nnnoqnTp1YvDgwbh7jefZtGkTl112GT169GD//ffnnnvuSeprVB21SKpTWlrzXzsikn6/+AXMmVO/9/bpU/3+Ll3gttvivvXZZ5+ldevWPP300wAsX76c0aNHA+E21ciRI5k9ezbNmjXj8MMPp6ysbMt7lyxZwvTp03nvvfc44YQTtrqlBfDkk09y/PHHMyf2/xo7diwA69ev57TTTmPcuHH06NGDb7/9lpKSElq2bMnzzz9PcXExH374IaeffvqWhPbmm2/yzjvv0Lp1a3r37s0rr7xCz549qz3P2LFjad68OTNmzGDdunX07t2bfv361XuEVnXUIqnO0qXgHrZVq6BpUxg2TEN/RfJc586deeGFF7j88suZNm0alQfqvPHGGxx22GHsuOOONGrUiEGDBm313pNOOomioiL22WcfPqtDy+j999+nVatW9OjRA4DtttuOhg0bsmHDBs455xw6d+7MoEGDtuqn6dmzJ23atKGoqIguXbqwcOHCGs/z3HPP8eCDD9KlSxd69erFsmXL+PDDD5P5Mn2PWiS12XZbOOUUmDAB7roLioujjkgk/9XSciDeKKOXX673x+61117MmjWLyZMnc+WVV9KvX78tr1XcPqpJ48aNEz62MnevdtTUrbfeSmlpKW+99RabN2+muNLvnsqf1aBBAzZu3FjjedydO++8k/79+yccU12pRZKIIUNg+XKINXdFJD8tXryYbbfdliFDhnDppZcye/bsLa/17NmTKVOm8PXXX7Nx40Yef/zxlHxmp06dWLx4MTNmzABgxYoVbNy4keXLl9OqVSuKiop46KGH2LRpU73O079/f0aPHs2GDRsA+OCDD1i1alVKYq+gFkkiDj88jBJ55BEYODDqaEQkTf2Y8+bN47LLLqOoqIhGjRoxevRoLr30UgB23XVXrrrqKnr16kXr1q3ZZ599SMUctW222YZx48ZxwQUXsGbNGkpKSnjhhRc4//zzGThwIOPHj6dv3740adKkXuc5++yzWbhwIeXl5bg7O++8MxMnTkw67sqsLk2wXNW9e3dPej2SSy6BUaNCP8kOO6QmMBHZYv78+ey9995RhxHXypUradq0KRs3buTkk0/mrLPO4uSTT446rJSo7utvZrPcvXtt783rW1v1LpFSncGDYf360FciIgXpf/7nf+jSpQv77bcf7du356STToo6pKygFkmi3GGffaBlS5gyJTWBicgWudAiyWdqkWSCWWiVTJ0KH38cdTQiIllDiaQuzjgj/Pu//xttHCIiWUSJpC46dICDDoKHHw63ukRERImkzoYMgXfegblzo45ERCQrKJHU1aBB0LBhmFMiIiJKJHW2005wzDGhn6SWmaYikl5LlsBhh6kMXtSUSOpj8GD49NMwgktEIjNyJEyfDtdfn8pzjqRTp04cddRRnH766fzxj3/k3nvvpUePHpSVlTFw4EBWr14NhPLx5513Hn379qVDhw5MmTKFs846i7333puhQ4duOWfTpk25/PLL6datG0ceeSRvvPEGffr0oUOHDjz55JMANZaNzwnunvdbt27dPKVWrXJv1sz9rLNSe16RAvbuu+9ueXzRRe6HHVbzVlRUUZ57662oqOb3XHRR7THMmDHDy8rKfPXq1f7tt996x44d/Q9/+IN/+eWXW465+uqr/Y477nB39zPPPNNPO+0037x5s0+cONGbNWvmc+fO9U2bNnl5ebm/+eab7u4O+OTJk93d/aSTTvKjjjrK169f73PmzPGysjJ3d1+1apWvWbPG3d0/+OADT/nvrVpU/vpXAGZ6Ar9jVWurPipXBB41ShWBRTKsZ09YsAC+/BI2b4aionDXeY89kjvv9OnTOfHEEykpKQFgwIABALz99ttcc801fPPNN6xcuXKrSroDBgzAzOjcuTOlpaV07twZgH333ZeFCxfSpUsXttlmG44++mgglKpv3LgxjRo1onPnzixcuBCADRs2MGLECObMmUODBg344IMPkvvPZJASSX0NGQIPPABPPQVVFrARkeTUVkUe4LzzYMyY8Hfc+vWhnurddyf3uV7DsP6hQ4cyceJEysrKuP/++3m5Uqn6ipLuRUVFW5V3LyoqYuPGjQA0atRoS4n3ysdVPiZe2fhspz6S+urbF1q10ugtkYh89hmcey68/nr4NxUd7gcffDCTJk1i7dq1rFy5cstKiStWrKBVq1Zs2LCBR9L0M1/XsvHZRC2S+mrQAE4/He68E776CnbcMeqIRArKE09893jUqNScs0ePHpxwwgmUlZXRtm1bunfvTvPmzRk5ciS9evWibdu2dO7cmRUrVqTmAyupa9n4bKKijcmYPRu6dYN77oHhw1N/fpECki1FGytKxa9evZpDDz2UMWPGUF5eHnVYaaeijTVIaRn56nTtCnvvHUqmiEheGD58OF26dKG8vJyBAwcWRBJJVl7f2nL3ScCk7t27n5OWD6ioCHzNNfDRR9C2bVo+RkQy539VlLXO8rpFkhGqCCwiBU6JJFnt20Pv3qoILJIChdBnm42S/borkaTCkCHw7rvw1ltRRyKSs4qLi1m2bJmSSYa5O8uWLUtq3kpe95FkzKBBcMEFYU5Jly5RRyOSk9q0acOiRYv44osvog6l4BQXF9OmTZt6v1+JJBVatIBjjw39JDfdFOaYiEidNGrUiPbt20cdhtSDbm2lyuDBsHgxTJkSdSQiIhmlRJIqAwZAs2aaUyIiBUeJJFVKSkLVuMcfhzVroo5GRCRjlEhSafBg+PbbUBFYRKRAKJGkkioCi0gBUiJJpYqKwJMnh4rAIiIFIOFEYmZNzEzjWmszZAhs2ADjx0cdiYhIRtSYSMysyMzOMLOnzexz4D1giZm9Y2Z/MLM9MxdmDunSRRWBRaSgxGuRvATsAVwJ7OLuu7l7S+AQ4HXgJjMbkoEYc4tZaJVMnw6xtZhFRPJZvERypLuPdPe57r65Yqe7f+Xuj7v7QGBc+kPMQaoILCIFpMZE4u4bKh6b2Q5mtq+ZdTCzouqOkUratYODD1ZFYBEpCPH6SJqb2VVmNo9wK+se4DHgIzMbb2Z9MxVkTho8GObPhzlzoo5ERCSt4t3amgB8Ahzi7j9w94Pdvbu77wbcBJxoZsMyEmUuGjQIGjXSnBIRyXtWCLX/u3fv7jNnzsz8B594IsyYAZ98oorAIpJzzGyWu3ev7bhay8ibWW9gjruvio3SKgdud/ePUhBn/tplF/jss/C4YaUvc2kpLF0aTUwiImmQyITE0cBqMysDfgV8BDyY1qgqiXXwjzWzCfH2ZZ2KJJLofhGRHJVIItno4f7XiYSWyO1As0RObmZ/MbPPzeztKvuPNrP3zew/ZnZFvHO4+wJ3H1bbPhERiUYiKySuMLMrgSHAobEyKY0SPP/9wF1UasHE3j8KOApYBMwwsyeBBsCNVd5/lrt/nuBniYhIBBJpkZwGrAOGuftSYFfgD4mc3N2nAlWrF/YE/hNrVawH/gac6O7z3P34Klt+JpEZM6KOQEQkZWpNJO6+1N3/5O7TYs8/dvdk+kh2JQwrrrAotq9aZtbCzP4MdI21jKrdV837hpvZTDOb+cUXXyQRbhocdBDccANs3Bh1JCIiSUvk1laqWTX7ahyD7O7LgHNr21fN+8YAYyAM/617mEkqLa2+Y71lSzjiCLj2Wnj2WXjoIWjfPuPhiYikShTrkSwCdqv0vA2wOII40mvp0lAeper22WehBtcjj8C8eVBWBg88oFIqIpKzokgkM4A9zay9mW0D/Ah4MoI4onXGGTB3LnTtCkOHwg9/CMuWRR2ViEidxau1tZ2Z3WhmD5nZGVVeuzuRk5vZo8BrwA/MbJGZDXP3jcAI4J/AfOAxd3+n/v+FuJ8/wMzGLF++PB2nT17btvDii3DTTfCPf8D++8MLL0QdlYhIndRYIsXMHgc+JBRsPAvYAJzh7uvMbLa7l2cuzOREViKlLmbPDoUe33sPLr4Yfvc7KC6OOioRKWCJlkiJd2trD3e/wt0nuvsJwGzgRTNrkbIo5Tvl5TBrFowYAbfeCj17hj4UEZEsFy+RNK6y9shvCaOgpgJKJumw7bZw550weTJ8/nm41WX2/W2XXaKOVERki3iJZBJweOUd7v4A8EtgfTqDSpWs7yOpyTHHxG+NqF6XiGQRlZHPZlbdlJuYArhuIhKtVJaR3x74CdCu8vHufmEyAYqISH5IZGb7ZMLIrXnA5vSGIyIiuSaRRFLs7pekPRIREclJicxsf8jMzjGzVma2Y8WW9shSIGc72yuUlla/v6hIHe4ikjUSSSTrCWXjXwNmxbac6Ll290nuPrx58+ZRh1I/1dXrmjsXttkGfvxj2Kw7jSISvUQSySVAR3dv5+7tY1uHdAcmNejcGe64A55/Hm6sug6YiEjmJZJI3gFWpzsQqYOzz4bTT4frroOpU6OORkQKXCKd7ZuAOWb2EmGlREDDfyNlBvfcAzNnhoQyZw7svHPUUYlIgUokkUyMbZJNmjWD8eOhV6/QXzJ5cuiEFxHJsEQSyQRgrbtvAjCzBkDjtEaVImY2ABjQsWPHqENJj7IyuO02OO88uPlmuOKKqCMSkQKUyJ+w/wJKKj0vAXJi0YycH7WViJ/9LCyKdc01MH161NGISAFKJJEUu/vKiiexx9umLySpEzO4915o1y70l3z5ZdQRiUiBSSSRrDKzLYtYmVk3YE36QpI622670F/y+edw5pmaXyIiGZVIIvkFMN7MppnZNGAcYalcySZdu8Kf/hQ63W+5JepoRKSA1NrZ7u4zzKwT8APAgPfcfUPaI5O6O/98ePlluPJKOPhgOPDAqCMSkQJQY4vEzA6ueOzuG9z9bXefV5FEzGw7M9svE0FKgszgvvugbVs47TT46quoIxKRAhDv1tZAM3vVzK4zs+PMrKeZHWpmZ5nZQ8BTbD2aK+vkfNHG+mjeHMaNC3W6hg7VAlgiknZxV0g0sx2AU4HeQCtCJ/t84Gl3z5mxpjm7QmIy7rgDLroo9JdcolUARKTuEl0hUUvt5it3GDgQJk0K80t69Yo6IhHJMYkmEtXUyFdmMHZsGAp8wAHheeVtl12ijlBE8oQSST7bYYea55RoYSwRSZF4o7ZaZTIQERHJTfHmkfwl1tn+MvAsMN3dN2YkKhERyRk1JhJ3P8bMioE+wMnAH83sY0JSedbdP85MiCIiks3izmx397XEEgeAmbUHjgHuMrNd3L1n+kOsv7wvIy8ikgXq1Nnu7v9197vd/QTg4FrfELGCKCNfm9LSuu0XEamjeo/acvf1qQxE0mTp0jCnpGIbNgwaNw7L9IqIpICG/xaaa64JQ4JvvDHqSEQkTyiRFJp27UKr5N574WONlxCR5NWaSMxsnpnNrbJNM7NbzaxFJoKUFLvqqjC7/YYboo5ERPJAIi2SZ4CngcGxbRIwFVgK3J+2yCR9dtsNhg+Hv/4VFiyIOhoRyXGJJJLe7n5lbC2See5+NdDH3X8PtEtveJI2V14JDRvCyJFRRyIiOS6RRNLUzLaUjjWznkDT2FPNdM9VrVvDeefBgw/Chx9GHY2I5LBEEskw4D4z+6+Z/Re4DzjbzJoAGvqTyy6/HIqL4frro45ERHJY3ERiZkVAB3fvDHQBurr7/u4+w91XuftjGYlS0qO0FEaMgEcegfnzo45GRHJU3ETi7puBEbHHy939m4xElSIFudRuXV12GTRpAr/5TdSRiEiOSuTW1vNmdqmZ7WZmO1ZsaY8sBVQiJQE77RSW5H3sMZg3L+poRCQH1brUbqxfpCp39w7pCSn1CnKp3br46ito3x6OPBIefzzqaEQkS6RsqV13b1/NljNJRBKw445w8cXwxBPw5ptRRyMiOSaRme3bmtk1ZjYm9nxPMzs+/aFJRl18MWy/Pfz611FHIiI5JpE+kr8C64GDYs8XAaqtkW+aN4dLL4VJk2DGjKijEZEckkgi2cPdbwY2ALj7GsDSGpVE48ILoUULtUpEpE4SSSTrzawEcAAz2wNYl9aoJBrNmsGvfgXPPAOvvRZ1NCKSIxJJJL8mLLW7m5k9AvwL+FVao5Lo/Pzn0LIlXHdd1JGISI5IZNTW88ApwFDgUaC7u7+c3rAkMk2awBVXwAsvwNSpUUcjIjmgxkRiZu0qHrv7Mnd/2t2fcvcvY6+bmbVJf4iSceeeC61awbXXhuV5RUTiiNci+YOZPW5mPzGzfc2spZntbmaHm9lI4BVg7wzFKZlUUhIWv5o6FV56KepoRCTLxZ3Zbmb7EBaz6g20AlYD84HJwAR3X5uJIJOlme31sHYt7Lkn7L47TJ8eVlQUkYKS6Mz2hvFedPd3gatTFpXkjuJiuPrqsGbJc89B//5RRyQiWSqRUVtSqM46C9q2VV+JiMSV14lEZeSTtM02IYnMmAFPPx11NCKSpWqt/psP1EeShA0boFOnUEJl1iz1lYgUkJRV/42N3DoutlqiFJpGjeDLL0NV4KKikEgqtl12iTo6EckCiSSH0cAZwIdmdpOZdUpzTJJtvv22+v2ffZbZOEQkKyUys/0Fdx8MlAMLCSsmvmpmPzWzRukOUEREsltCt6vMrAWhRMrZwJvA7YTE8nzaIhMRkZwQdx4JgJk9AXQCHgIGuPuS2EvjzEw92CIiBa7WRALc5+6TK+8ws8buvi6R3nzJc//9b1jvXUQKViK3tqpbDVGLVRSS0tLq95tBeTlMnJjZeEQkq8Sr/ruLmXUDSsysq5mVx7Y+wLYZi1Cit3RpmNledfu//4OOHeHkk+GSS2D9+qgjFZEIxLu11Z/Qwd4G+FOl/SuAq9IYk+SK9u1DQcfLLoNbb4VXX4Vx40JZFREpGLXObDezge7+eIbiSQvNbM+Axx8PtbmKiuCBB+CEE6KOSESSlPTMdjMbEnvYzswuqbqlLFLJDwMHwuzZ0KEDnHgiXHppKK8iInkvXmd7k9i/TYFm1WwiW9tjj3B76+c/h1tugUMPhY8/jjoqEUmzRG5t7ezuX2QonrTQra0IjB8Pw4ZBw4bw4INw/PFRRyQidZSyoo3Aq2b2nJkNM7MdUhCbFIJBg8KtrrZtYcAAaNJk64KPKvwokjcSqbW1J3ANsC8wy8yeqtR/IlKzjh3htdfCKourV1d/jAo/iuS8hGptufsb7n4J0BP4CnggrVFJ/iguhrvvjjoKEUmjRNYj2c7MzjSzZ4BXgSWEhCIiIpJQra23gInA9e6u0igiIrKVRBJJBy+E9XglOu5awlckh8WbkHhb7OGTZva9LUPxSb6oqfAjhM74TZsyF4tIgViyBA47LJTLS6d4LZKHYv/+Mb0hxGdmHYCrgebufmps30nAcUBLYJS7PxdhiJKI6r6T3eGaa+B3v4Nly+Dhh6Fx48zHJpKnrr8+lMO7/vr0jnmpsUXi7rNiD7u4+5TKG9AlkZOb2V/M7HMze7vK/qPN7H0z+4+ZXRHvHO6+wN2HVdk30d3PIRSVPC2RWCQLmcFvfxsKPk6YAMceCytWRB2VSM4rKQk/Xn/+M2zeDKNHh+clJen5vESG/55Zzb6hCZ7/fuDoyjvMrAEwCjgG2Ac43cz2MbPOsTkqlbeWtZz/mti5JJf94hfw0EMwdSr07Quffx51RCI5be5c2Gmn755vuy0MHhzWoUuHGm9tmdnpwBlA+yp9Is2AZYmc3N2nmlm7Krt7Av9x9wWxz/kbcKK73wgkVEfDzAy4CXjG3Wcn8h7JckOGwI47wqmnwiGHwHPPqRy9SD0sXw4/+Um4W2wW7havXQvbbZe+QhLx+kgq5ozsBNxSaf8KYG4Sn7kr8Eml54uAXjUdbGYtgN8CXc3syljCuQA4EmhuZh3d/c/VvG84MBxg9913TyJcyZhjj4UXXoDjjoODDgrJZN99o45KJGd8/TX07w9z5kDPntCtGwwfDmPGhI73dKm1aGPSHxBaJE+5+36x54OA/u5+duz5j4Ge7n5BumJQ0cYc8/bb4adhzRp46qmQVEQkrmXLoF+/8OMzYUIocZesVKxHMj327woz+7bStsLMvk0itkXAbpWetwEWJ3E+yTf77QevvBJu8h55JDzzTNQRiWS1L76AI46Ad96BiRNTk0TqIt6orYNj/zZz9+0qbc3cfbskPnMGsKeZtTezbYAfAZqXIltr1y6MW9x777Da4iOPRB2RSFb6/HM4/HB4/3148kk45pjMx5BIra09zKxx7HEfM7vQzLZP5ORm9ijwGvADM1tkZsPcfSMwAvgnMB94zN3fqf9/Ie7nDzCzMcuXL0/H6SXdWraEl14KC2QNGQK33x51RCJZZckS6NMHFiyAp58Ot7aikMjCVnOA7kA7wi//J4EfuPuxaY8uRdRHkuPWrg1jF594ovrXS0vTP3VXJMt8+mloiXz6KUyeHP7eSrVULmy1OdaKOBm4zd0vBlolG6BIwoqL4bHHan5da5pIgfnkk1D6ZMkS+Oc/05NE6iKRoo0bYnNKzgQqunAapS8kkWo0aBB1BCJZYeHC0BJZtiyMkD/ggKgjSqxF8lPgQOC37v5fM2sPPJzesFJDfSQikk8WLAgtka+/hn/9KzuSCGRgHkk2UB9JnohXar4Avo+lsH34YWiJrF4d5u127Zr+z0xZH4mZ9Taz583sAzNbYGb/NbMFqQlTJEXi9aGI5Lj33w8tkbVr4cUXM5NE6iKRPpKxwMXALECLRkh0Skur71hv1AhOOw0++ACuvlqLZEleeffd0BJxD6Ph99sv6oi+L5FEstzdNbVYolfTEN916+Ccc+Daa+G99+C++8JIL5EctmRJmKG+cGH4W+mll8L83GyUSCJ5ycz+ADwBrKvYmQtVd81sADCgY8eOUYci6dS4MTzwAHTqFFokCxaEOhEta1uFQCR7XXQRzJoFTZrAq6/CXntFHVHNEpmQ+FI1u93dD09PSKmnzvYCMmFCqKHdsmUo+JiN9wFE4igpCX0hVRUXhzqmmZSyznZ371vNljNJRArMqaeGBbLWrw9VgydPjjoikYQ98wy0ik33rpg6le5FqVIhkVFbpWY21syeiT3fx8yG1fY+kch07w5vvAEdO4abzLffruHBktX+7/9CbdJjjw39IcceG75li4vTvyhVKiQyIfF+Qo2t1rHnHwC/SFdAIinRpk1omZxwQljK9/zzYcOGqKMS2crq1WGMyL77hs703/8e5s0L3X7nnguvvx7+zfZScon0kcxw9x5m9qa7d43tm+PuXTISYQqoj6SAbd4MV10VfkKPOALGj4cddog6Kilw7qE775e/DHWzBg+Gm2+G1q1rf28mpbJo46rYcrceO/EBQE7UHFGJFKGoCG66Cf7619BC2XnnMM+k6pbN9w0kr7zzTvib5oc/hB13hGnT4OGHsy+J1EUiieQSQun4PczsFeBBwprpWc/dJ7n78ObNm0cdikRt6NBQV2JTDXNqVUFY0uybb8Jd1rKysKb6qFFheO/BB0cdWfJqnUfi7rPN7DDgB4AB77u7bjZL7om61rYUpM2b4f774corw5K4w4fDDTeElaTzRbw123uY2S4AsfVIugG/BW4xsx0zFJ+ISM5ZsiTUxnrmGTjwQBg2LAwinDkT/vzn/EoiEP/W1j3AegAzOxS4iXBbazkwJv2hiWTYunW1HyOSgKuuCl1yxx4LH38MDz4I06dDeXnUkaVHvETSwN2/ij0+DRjj7o+7+7WAao5I/unaFV57LeooJIeVlISxG/ff/92+pUvD7ax8riUaN5GYWUUfyhHAi5VeS6RGl0j2KS2tfv/228PKldC7dyhytHJlZuOSvHDXXWFGekXSyIVZ6akQL5E8Ckwxs38Aa4BpAGbWEQ3/lVy1dGkYxF91+/rrMC7z/PPhjjtCja7nnos6WskR7vCnP4Ui1DvGepBzZVZ6KtSYSNz9t8AvCTPbD/bvZi4WoeG/ko+aNQt/Uk6bFn4L9O8fhg1/9VWtb5XCtXEjjBgRJheeckroXD/vvNyZlZ4KWmpXpDpr14Yxmr//ffgT8667QkHIfL7RLXW2YgX86EehNuhll4W5r0WJzM7LEamc2S5SeIqLQyKZORN22y1MQz7lFFi8OOrIJEt8+mmYmvTPf4YhvTffnF+/v5RaAAALdklEQVRJpC4K9L8tkqCysnCP4uab4dlnYZ99oHlzlVkpcG+9Bb16haq9Tz0FP/tZ1BFFS4lEpDYNG4b7FvPmhSHC335b/XEqs1IQJk8OZU3MwtyQo4+OOqLoKZGIJKpjR3jxxdqPk7w1enRY4mbPPeHf/4b99486ouygRCJSF+psL0ibN8Oll4bR4cceG2at53K13lTL60SieSSScePHazXGPLN6dRiwd8stYZjvxInQtGnUUWWXvE4kmkciGffDH4aJBFOnRh2JpMBnn0HfviF53HYb3Hnnd2upy3fyOpGIpEVNZVZKS+Evf4FFi0Lp1wED4O23MxubJK2icu+UKXDAAeES/v3voXKOVE+JRKSuaiqzsnQp/PSn8OGHYWbatGlh+PCwYSG5SE4YOTJcun79wrzUKVPgxBOjjiq7aWa7SLosWwa/+12YFV9UFJbHu/zyUCBSsk5JSUgcVRUXw5o1mY8nG2hmu0jUWrQIPbTvvw8DB4ZWyh57wK23hsmLmtSYVV58ceu7liUlhVG5NxWUSETSrV07ePhhmD0bunWDSy6pefKiJjVm3MaNId8fcUSoz2kWWiHr1hVG5d5U0LoiIpnStWsoTf/88+EGvERu7tzQhTVzJpxwAqxfDx06hIWoxowJHe9SOyUSkUw76qioIyh469aFmpw33QQ77ADjxsGgQVvPNx01Krr4co1ubYlkmxtu0BooafTKK9ClS/gyn3EGzJ8fpv+oaEH9KZGIZJtrr4Xddw8TFxYujDqavLFiBVxwARxySJit/uyz8MADYUyEJCevE4lKpEjWijepce7cMMrr7rtDocjBg2HOnMzGl2eeeQb23TfcrrrggrCqcv/+UUeVP/I6kahEimSteJMaO3cOfyovWBDmnjz5ZOio79cvdNQXwNyvVPnyS/jxj0OhxaZNw22t229XraxUy+tEIpLTdtsN/vhH+OQTuPHGsB5Kv35QXh4mNWoeyvdUlDdZsgQefRT23hv+9je47jp4881QBk1ST4lEJNttvz1ccUXoLxk7Nky/rul2bYHPQ6kob3LggaEjvUOHMH3nN7+Bxo2jji5/qUSKSK7ZvDl+CdoC+JmuSuVN0kMlUkTyVVEtP7annAIPPlgwQ4hXrQorIVducRQXq7xJJimRiOSbN96AM8+Eli3hyCND0cg8rD68ciX8/vehAs3IkbDTTt+VN1m/XuVNMkmJRCTffPxxWFD8V7+CTz8N41132w169AjViOfPD7e/crRw5IoVYexBu3ah66i8PIzG6tkTzjsPXn8dzj03DICTzFAfiUgu2mWX6jvWS0u//xv0vffCEn8TJ4YEA7DXXvDBBzWfPwt/LyxfHlYovPXWcNfumGPCaKwDDog6svylPhKRfBZvHkpVnTqFP91ffz3c4ho1Ctq2zXzM9fTNN3D99aEFcu21cNBB4e7d5MlKItlCiUSkkOy6K5x/fqhCHE91Q6DSrGIOSEUu/Ppr+PWvQwL59a/h0ENDld5Jk8JdOskeqv4rIt/Xpk1YNvjcc8NiXBkwciRMnw5XXRXy3R13wLffwsknh5ZI164ZCUPqQX0kIoUqXrnbU0+Fv/8dNm0KRanOOw+OOw4apv5vz5rmgBQVhcmEZWUp/0hJkPpIRCS+eIUjx48Po79+8xt4+2046aQwTfyGG1I+HGrBgjALvWKOZYMGIWd9+qmSSK5QIhEpVLV12LduHYZFLVwITzwROu2vvTYMJT7tNHj55ZQMIW7VKsz52Lw5TCp0D1X0s3wUslSS14lEZeRFUqBhw9BR8dxzYcjwhReGKsR9+6Zs7fnPPgt3z/79b80ByUXqIxGRuluzJqxP+9Of1nxMAfxuyXfqIxGR9CkpgaFDo45CsoQSiYiIJEWJREREkqJEIiL1F28IsRQMzWwXkfrT8CpBLRIREUmSEomIiCRFiURERJKiRCIiIklRIhERkaQURIkUM/sC+CjJ0zQHki3aVddzJHJ8bcfEe7261xLdtxPwZS2xpVq2XoNEjqvp9brsr7ovV69Bfc6Tzp+FZK4B5O51SOQcbd1951rP5O7aEtiAMZk+RyLH13ZMvNere60O+2bqGiR/Heqyv+q+XL0G6boOUVyDXL4OqbqW7q5bW3UwKYJzJHJ8bcfEe7261xLdF4VsvQaJHFfT63XZnw3XIVUxZNPPQq5dA4jmZ6FGBXFrS1LPzGZ6AlVBJX10DbKDroM626X+xkQdgOgaZImCvw5qkYiISFLUIhERkaQokYiISFKUSEREJClKJJJyZtbBzMaa2YSoYykkZtbEzB4ws3vNbHDU8RSiQv3eVyKRrZjZX8zsczN7u8r+o83sfTP7j5ldEe8c7r7A3YelN9LCUMfrcQowwd3PAU7IeLB5qi7XoFC/95VIpKr7gaMr7zCzBsAo4BhgH+B0M9vHzDqb2VNVtpaZDzmv3U+C1wNoA3wSO2xTBmPMd/eT+DUoSFohUbbi7lPNrF2V3T2B/7j7AgAz+xtworvfCByf2QgLS12uB7CIkEzmoD8SU6aO1+DdzEaXHfTNJonYle/+0oXwC2vXmg42sxZm9megq5ldme7gClBN1+MJYKCZjSZ7Snnkq2qvQaF+76tFIomwavbVOJPV3ZcB56YvnIJX7fVw91XATzMdTIGq6RoU5Pe+WiSSiEXAbpWetwEWRxSL6HpkA12DSpRIJBEzgD3NrL2ZbQP8CHgy4pgKma5H9HQNKlEika2Y2aPAa8APzGyRmQ1z943ACOCfwHzgMXd/J8o4C4WuR/R0DWqnoo0iIpIUtUhERCQpSiQiIpIUJRIREUmKEomIiCRFiURERJKiRCIiIklRIpGCZ2abzGxOpS1umfxMMrMJsTUu/h2L7WMz+6JSrO1qeN8NZjayyr7uZjY39vhfZtY8/f8DKQSaRyIFz8xWunvTFJ+zYWzSWjLn2Be4wd1PrrRvKNDd3Uck8N6/u/telfb9EVjm7jea2TBgJ3f/fTIxioBaJCI1MrOFZvYbM5ttZvPMrFNsf5PYYkczzOxNMzsxtn+omY03s0nAc2ZWZGZ3m9k7sbVaJpvZqWZ2hJn9vdLnHGVmT1QTwmDgHwnEeYyZvRaLc5yZNYnNsl5rZt1ixxgwCPhb7G3/AM5I5usjUkGJRARKqtzaOq3Sa1+6ezkwGrg0tu9q4EV37wH0Bf5gZk1irx0InOnuhxNWLGwHdAbOjr0G8CKwt5ntHHv+U+Cv1cTVG5gVL/DYQmJXAEfE4pwLXBR7+VFCDaiKcy129/8CuPuXQDMz2z7e+UUSoTLyIrDG3bvU8FpFS2EWITEA9ANOMLOKxFIM7B57/Ly7fxV7fDAw3t03A0vN7CUItcbN7CFgiJn9lZBgflLNZ7cCvqgl9oMIK/S9GhodbANMj732KDDFzH5FSCiPVnnvF7HP+KaWzxCJS4lEJL51sX838d3PiwED3f39ygeaWS9gVeVdcc77V8LiU2sJyaa6/pQ1hCQVjwHPuvuPq77g7gvNbDFwCHAy0K3KIcWxzxBJim5tidTdP4ELYv0OmFnXGo6bTlixsMjMSoE+FS+4+2LC+hXXENYEr858oGMtsbwKHGZmHWKxNDGzPSu9/ihwBzDf3ZdW7DSzImAntl7lT6RelEhEvt9HclMtx48EGgFzzezt2PPqPE5YAOlt4B7g38DySq8/Anzi7jWt8/00lZJPddz9M2AYMM7M3iIklr0qHfIYsB/fdbJX6AlMd/dN8c4vkggN/xVJIzNr6u4rzawF8AbQu6JlYGZ3AW+6+9ga3lsCvBR7T0p/4ZvZKMIaGlNSeV4pTOojEUmvp2Ijo7YBRlZKIrMI/Sm/rOmN7r7GzH4N7Ap8nOK43lQSkVRRi0RERJKiPhIREUmKEomIiCRFiURERJKiRCIiIklRIhERkaQokYiISFL+H2UAApYVRf4RAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot the sensitivity curve\n", "t = sensitivity_estimator.results_table\n", "\n", "is_s = t[\"criterion\"] == \"significance\"\n", "plt.plot(\n", " t[\"energy\"][is_s],\n", " t[\"e2dnde\"][is_s],\n", " \"s-\",\n", " color=\"red\",\n", " label=\"significance\",\n", ")\n", "\n", "is_g = t[\"criterion\"] == \"gamma\"\n", "plt.plot(\n", " t[\"energy\"][is_g], t[\"e2dnde\"][is_g], \"*-\", color=\"blue\", label=\"gamma\"\n", ")\n", "\n", "plt.loglog()\n", "plt.xlabel(\"Energy ({})\".format(t[\"energy\"].unit))\n", "plt.ylabel(\"Sensitivity ({})\".format(t[\"e2dnde\"].unit))\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We add some control plots showing the expected number of background counts per bin and the ON region size cut (here the 68% containment radius of the PSF)." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0.01, 0.5)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEOCAYAAAD7WQLbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcTfX/wPHXewbDKAotxjJja1GJIkRCslRSIZXqW5apviX6auGraOFrKUoqGoVvKKFUlr4qiaRsSbb8spNQFGM35v3749xhjDszZ2bucubO+/l4nIe5595zzptb3j7rW1QVY4wxJtJEhTsAY4wxJhgswRljjIlIluCMMcZEJEtwxhhjIpIlOGOMMRHJEpwxxpiIZAnOGGNMRLIEZ4wxJiJZgjPGGBORLMEZY4yJSIXCHUAgiUhroDXQNTY2NtzhGGNMvnHo0CEF3gGmq+r0cMcTCBKJe1EWL15cDx48GO4wjDEm3xCRQ6paPNxxBFJEtuBiYmLCHYoxxpgwsxacMcaYiGzB2SQTY4wxEcm6KI0xxkSkiGrBqep0VU2Mjo4Oy/MnTpxIQkICUVFRJCQkMHHixLDEYYwxJsJacOE0ceJEEhMTOXToEABbtmwhMTERgI4dO4YzNGOMKZAiapJJui7KrkeOHAnpsxMSEtiyZcsZ5+Pj49m8eXNIYzHGmJyKxEkmEZXg0oRjFmVUVBT+/ixFhNTU1JDGYowxORWJCS6ixuDC5ejRo2S1c8qjjz7K0qVL/SZAY4wxwRFRCU5EWotI0okTJ0L2zJ07d9KkSRMOHjxI4cKFT3svJiaGevXqMWbMGOrUqcOVV17Jq6++yh9//BGy+IwxpqCKqAQX6lmUP/74I3Xq1GHFihVMmTKFsWPHEh8fj4gQHx/Pu+++y8KFC/n9998ZOXIkxYoV41//+hdxcXHccccdTJ8+nZSUlJDEaowx2YgWkSTfXIbIoKoRd8TGxmqwTZo0SYsVK6YVKlTQ5cuXu75u1apV2rNnTz3//PMV0AsuuECfeuopXbNmjaqqTpgwQePj41VEND4+XidMmBCs34IxxpwEHFQP/P0dyMMmmeRQamoq/fr1o3///lx77bV8/PHHXHDBBTm+z/Hjx5k1axZjx45l5syZpKSkUKVKFbZt28axY8dOfi42NpakpCRbamCMCapInGRiCS4HDhw4wP3338+0adPo1KkTb731FoHYNWXXrl1MmDCB3r17c/z48TPet6UGxphgswTnccFcB7d582ZuvfVWVq9ezdChQ+nevTsiEtBn2FIDY0y4RGKCs0kmLsybN486deqwbds2Pv/8c3r06BHw5AZQsWJFv+crVKgQ8GcZY0yki6gEFwxJSUk0a9aM0qVLs2jRIpo3bx60Zw0YMMDveroSJUqQnJwctOcaY0wksgSXiePHj9OtWzceeughmjVrxg8//MBFF10U1Gd27NiRpKSk05YadO3albVr19KkSRN27doV1OcbY0wkiagxuDS5mWQyceJE+vTpw9atWylXrhwlSpRgzZo1/Otf/2LIkCGEq0IBwMyZM2nfvj1xcXHMnj2bKlWqhC0WY0xkisQxOEtwnFkJIE1iYiJvv/12oMPLlUWLFnHzzTcTHR3NrFmzuPrqq8MdkjEmgliC87jczqLML5UA1q1bR4sWLdizZw8ff/wxN954Y7hDMsZECEtw+UROW3D5aXr+jh07uOmmm1i9ejXjxo2zBeDGmICIxARnk0w4NT2/cCbnvSQuLo558+bRsGFD7r33XoYOHRrukIwxxpMswXFqev5XwBfA3UCpYsUYMGBAmCPzr2TJkvzvf/+jffv2PPnkk/Ts2dNzLU1jjAm3QuEOwAs6duwIqaks7daN2/ft433gmAhFvv0WqlSBunUhCAu78yImJoZJkyZx4YUXMmzYMH7//XfGjRtHkSJFwh2aMcZ4go3BZZSaCvPmwdixMHUqHD4Ml14KDzwA990HZcsGNNa8UlUGDx5M7969adasGR999BElSpQId1jGmHzGxuAKgqgoaNIE3nsPdu6E0aPh3HPhmWegQgW45Rb46CNIt+N/OIkIvXr1YuzYscydO5fGjRvz1ltvkZCQQFRUFAkJCUycODHcYRpjvC/i6sF5vgUnIpcC3YEywBxVHZndNUGpJrBuHYwb5yS+HTugdGno2BEefBBq1gzss3Lp888/p02bNqSkpJw2K9RK7hhjshOJLbhsE5yINAB+UtWDInIvcBUwXFXPXDjm9qEiY4BbgN2qenm68y2B4UA08I6qDkr3XhQwWlU7Z3f/YNaDIyUFvvzS6cL89FOnJVezppPo7rkHypQJznNdKlu2LDt37jzjvNfW9BljvCUSE5ybLsqRwCERuRJ4GtgCvJfH544DWqY/ISLRwJtAK6A6cLeIVPe9dyuwAJiTx+fmXaFC0KoVTJ7stORGjHC6Nbt3h7g4aNsWvvsubOFltl/l1q1bQxyJMcaEl5sEl+IrZ94Gp+U2HDg7Lw9V1fnA3gynrwHWq+pGVT0GTPI9E1X9TFWvBbzVx1a6NDz2GCxbBitWwKOPwvz50LAh3HYbrF0b8pAyW7tXrly5EEdijDHh5SbBJYtIb+BeYKavpZVxTXQglAO2pXu9HSgnIo1F5HUReRuYldnFIpIoIktFZGlKSkoQwstGjRrw6quweTP07w9ffw2XXw4PPQS//x6yMDIruXPgwAEWLFgQsjiMMSbc3CS4DsBRoLOq7sRJRC8HIRZ/C81UVb9R1cdV9SFVfTOzi1U1CXgB+DEYxUhdK14c+vSBDRucFt3YsVC1Kjz3HOzfH/TH+yu5079/f0qXLk3jxo0ZPHiwLQo3xhQIbiaZDFbVZ7I7l+MHiyQAM9ImmYhIfeB5VW3he90bQFUH5vTeQZ1kklMbNjgJ78MP4bzzoG9fSEyEEC/I3r9/P127dmXy5Mm0atWK9957jzJhnhBjjPGOgjrJxN+W9a0CHQiwBKgmIpVEpAhwF/BZTm4gIq1FJOnEiRNBCC+XqlSBSZNg8WK47DLo1g2qV3cmqYRwiUaJEiWYNGkSb731FnPmzKFmzZrWZWmMiWiZJjgReUREVgIXi8jP6Y5NwMq8PFREPgC+9917u4h0VtUU4DFgNrAWmKyqq3NyX1WdrqqJ4SxOmqk6dZxxuZkzoVgx6NAB6tVzdk0JERHhkUce4YcffqBo0aI0btyYQYMGWZelMSYiZdpFKSIlgXOBgUCvdG8lq2rGGZCekNt6cCF34oSzYLxvX9i+3dkdZdAgp4UXIum7LFu2bMn48eOty9KYAiwSuyhd7WTimzl5Aek2Z1ZVzy6s8tQYXFYOH4bXX4eBAyE52dnv8oUXoHz5kDxeVRk1ahQ9evTgvPPOY9KkSTRs2DAkzzbGeEskJrhsx+BE5DFgF/AlMNN3zAhyXAVDsWLOHpcbNjgLxSdMgGrV4N//hgMHgv5467I0xkQyN5NMegAXq+plqnqF76gR7MByw5OTTNwoXRqGDXP2u2zb1mnRNW0Ke0PTE1yrVi1+/PFH2rZtS+/evbn55psZNWqUbdhsjMnX3CwTmAvc6JsEki/kmy7KzHz2GbRv75Tp+fJLZ3lBCKR1WXbr1o3U1FTbsNmYAiQSuyjdJLh3gYtxuiaPpp1X1WHBDS338n2CA/jiC2e7r0qV4KuvQlqHLi4ujt/97L5iGzYbE7kiMcG56aLcijP+VgRnD8q0w3PybRelP82bw+efw5Yt0KgRbNuW/TUB4q8aAdiGzcaY/MXz9eByIyJacGm+/x5atoRSpWDOHKhcOeiPTEhIYMuWM6shxcXF8dtvvwX9+caY0CuQLTgRmSsiX2c8QhGcAerXdxaI79/vtOTWrQv6IzPbsHnPnj18/PHHQX++McYEgpsuyieBp3zHc8BPwNJgBpVbEdVFmd7VV8PcuXD8OFx/PaxaFdTH+duwefjw4Vx55ZW0bduWZ555hrBUbDDGmBzIVReliMxT1euDEE9ARFQXZXq//AI33ABHjzqTUK66KqSPP3r0KN27d+ftt9+madOmTJo0ifNCNMPTGBNcBbWLslS6o4yItAAuDEFsJqNLLnEKqhYv7qyTW7QopI+PiYlh1KhRjBkzhu+++46rrrqKxYsXhzQGY4xxy00X5TKcLsllOBsk9wQ6BzMok4UqVeDbb6FMGWjWzEl4Ifbggw+ycOFCoqOjue6660hKSiISJysZU8BEi0iSb0/fiGCzKPOrHTuc7sotW5yF4c2ahTyEPXv20LFjR2bPns2DDz7Im2++SbFixUIehzEm7wpqF2VhEXlcRKb6jsdEpHAogsupiJ1k4k9cnFNqp1o1pxrBzJkhD6F06dLMnDmT5557jrFjx9KwYUNbCG6M8Qw3O5m8AxQG/us7dR9wQlW7BDm2XCsQLbg0e/dCixawYoVTWPWOO8ISxvTp07nvvvuIjo7m/fffp0WLFmGJwxiTO55swYmcDzQA4oDDwCpgKaqudoR3MwZXR1X/oapf+44HgTq5DtgEVqlSzlZederAnXfC+++HJYzWrVuzZMkS4uLiaNWqFe3btyc+Pt42azbG5JxIE0Rm42wR2QooC1QHngVWIvICIiWyvY2LFtyPQHtV3eB7XRmYqqqhnaOeAwWqBZfmwAFo3drptnznHejUKSxhHDx4kBtvvJHvv//+tPO2WbMx3uapFpzIy8AI/NUdFSkE3AJEo/pRlrdxkeBuAMYCGwEB4oEHVXVu7iIPvgKZ4MApoHr77TB7NrzxBjz6aFjCiI+P97tvpW3WbIx3eSrBBYjbit4xOBUFBPhFVY9mc0lYFdgEB84i8A4d4NNPoU8fp0J4dHRIQ4iKivK7bEBErJiqMR7lyQQn8i8/Z/cBy1D9KbvL3cyifBQopqo/q+oKIFZE/pnzSIOvQM2izExMDEyZAp07w4ABcPPNISucmqZixYp+z9uuJ8aYHKoNPAyU8x2JQGNgNCJPZ3exm0kmXVX177QXqvoX0DVXoQaZqk5X1cToELdYPKdwYRg9Gt5+29mouXZt+Cnbf+wEjL/NmkWEPXv2MGXKlJDFYYzJ90oDV6HaE9WeOAnvPKAR8EB2F7tJcFEiImkvRCQapzac8TIRSEx0dj05dsypSjB+fEge7W+z5lGjRlGvXj06dOjAyJEjQxKHMSbfqwgcS/f6OBCP6mHSFeDOjJtJJi8DCcAoQHGai9vUyaaeVKDH4PzZtQvuugu++QYeewyGDoUiof83yqFDh+jQoQMzZsygX79+9OvXj3T/djLGhJFHx+CeA24HPvWdaQ18BgwFklDNclq2mwQXhdPv2QxnkskXwDuq6tmBLktwfqSkQK9eTnK79lpnnC4uLgxhpNC1a1fGjRvHI488wogRIyjwXcrGeIAnExyAyNVAQ5z8swBV1+XabC/KgubDD50JKGef7SS5hg1DHoKq0qtXL4YMGUL79u0ZP348MTExIY/DGHOKhxNcQ6AaqmMROQ84C9VNbi51MwZnIkmHDvDDD06Ca9IERoyAEP8jR0QYPHgwr7zyClOmTOHmm28mOTk5pDEYY/IBkX7AM0Bv35nCwAS3l3s+wYnIbSIyWkQ+FZHm4Y4nIlx+OSxZAjfdBI8/DvffD4cOhTyMnj178t///pdvvvmGxo0bs3v37pDHYIzxtNuBWwGnS051B3C224vDkuBEZIyI7BaRVRnOtxSRdSKyXkR6AajqJ6raFWdKaIcwhBuZSpaEadPgpZdg4kRnXG7jxpCHcf/99/PZZ5+xdu1aGjRowKZNrnoejDEFwzHfrhFON5NIjrpQ3UwymX7y5qfswymC+raqHsnJA333bAQcAN5T1ct956KB/wNuBLYDS4C7VXWN7/2hwERV/TG7+9sYXA59/jl07Oh0Vb7/PrRqFfIQvv/+e26++WaKFi3K//73P2rUqBHyGIwpyDw5BifyJFANJy8MBDoB76M6ws3lblpwG3GS0WjfsR/YBVzke51jqjofyLi9xjXAelXdqKrHgElAG3EMBj53k9xMLrRqBUuXQny8s/PJSy9BiLfUql+/Pt9++y1RUVE0atSIb7/9NqTPN8Z4kOorwFTgI5ztIvu6TW7gLsHVUtV7fLuETFfVe4FrVPVRIJAVBcoB29K93u471w1niUI7EXk4s4tFJFFElorI0pSUlACGVUBUrgwLFzotub594bbb4O+/s78ugC677DIWLlzIhRdeSPPmzenZsycJCQlWcseYgkz1S1SfQvVJVL/MyaVuuijXAi3UV7ZARCoCs1X1UhFZrqq1chOziCQAM9J1Ubb3PaeL7/V9OIm0W07vbV2UeaAKb74JTzwBCQkwbhw0aBDSEP7880/q1q3LxgxjglZyx5jg8VQXpUgyZw6NnaKabS04cNeC6wksEJG5IvIN8C3wpDiDff/N8sqc2Q5USPe6PLAjJzewzZYDQMTZ7WTuXKcyQcOG8I9/OLuhhEiZMmU4fvz4GecPHTpEnz59QhaHMSZMVM/2JbHXgF44vXnlcZYM9Hd7m5yUy7kEX7kcQPNaMsdPC64QziSTG4DfcCaZ3KOqq3N6b2vBBcjBg05FgldegWLF4MUXnRpzhQoF/dFWcseY0PJUCy6NyCJU62Z7LhNuyuWMUdWjqrpCnfo70cCsXAV76p4fAN8DF4vIdhHprKopwGPAbGAtMDmnyc1acAFWvDj85z+wapWzWXOPHlCrFsyfH/RHZ1Zyp3z58kF/tjHGM04g0hGRaESiEOkIuP4L3k0X5W8iMhJARM4FviQHK8n9UdW7VbWsqhZW1fKq+q7v/CxVvUhVq6jqgFzc18rlBMNFFzlLCaZNg+RkuP56ZzLKjhz1IOeIv5I74Gzz5a9auDEmIt0D3Ikzc38X0N53zpVsE5yqPgfsF5FROBstD1XVsbmLNbisBRdEIs7MyjVr4Lnn4KOP4OKLnc2b/YyX5ZW/kju9evUiOTmZunXrsnSp6/1WjTH5lepmVNugWgbV81C9DdXNbi/PdAxORO5I/xJ4DlgM/M95rn6c+6iDy8bgQmDDBujeHWbOhEsvhTfegKZNg/7YNWvWcPPNN7Nr1y7ef/99brvttqA/05iCwFNjcCLPAm+hmnG9dNr7TYFYVGdkeZssElxWrTRV1U4uQw0ZEWkNtI6Jiel65EiON1gxuTFjhpPoNm6EO+90WnRBHifbtWsXbdq0YfHixbzyyis88cQTVlfOmDzyWIJrAzwNHAF+BP4AiuLsalIT+Ar4D6p/ZHkbK5dj8uzIERgyBAYOhOhopwvziSeCWlT18OHD3H///UydOpWHH36YESNGUCgEszuNiVSeSnBpRKoBDYCywGGcCYjzfRW9s7/cxULvysBwoB7OwrvvgR7qsh5POFiCC5NNm5zE9umnzsSUESOgefAKQKSmptKnTx8GDRpEixYtmDx5MiVKuFr/aYzJwJMJLo/czKJ8H5iMk0HjgCk4+0R6jk0yCbNKleCTT2DWLGcvyxYtoG1b2LYt+2tzISoqioEDBzJ69GjmzJlDw4YNbYalMeYkNy24RZphUZ2I/KCq9YIaWR5YC84Djh51xuP693e6Lfv3d3ZICdISjq+++op27dpRrFgxpk+fTu3atYPyHGMiVYFqwYlIKREpBcwVkV4ikiAi8SLyNDAzdCGafCkmBv79b1i92tnuq0cPqFsXli0LyuOaNWvGwoULKVq0KI0aNeKTTz4JynOMMflHVl2Uy3BqvnUAHgLmAt8AjwAPBj2yXLAuSg+qVMnpsvzwQ/jtN7jmGmecLjk54I+qXr06P/zwAzVq1OCOO+5g2LBhfrf7MsbkEyJDECmBSGFE5iDyJyL3ur48Ev8CsC5Kj/r7b6dVN2oUlCvnrJ1r0ybgj0k/w/KGG27g119/Zdu2bVSsWJEBAwZYNQJj/PBkF6XIT6jWROR24DbgCWAuqle6udzNJBNjAuOcc+Ctt+C77+Dcc52dUW6/PeCTUIoVK8aHH37ILbfcwpw5c9i6dSuqypYtW0hMTLS6csaEk0hlRN5FZKqLTxf2/XoT8EGmC78zYQnOhF79+s5Y3ODBMHs2VK8Ow4dDALuWo6KiWLly5RnnreSOMXkgMgaR3YisynC+JSLrEFmPSK8s76G6EdXOLp84HZFfgNrAHETOw1n87S7cLHYyaaCq34lITF5L44SadVHmI5s2wT//Cf/7H1x9Nbz9tvNrAFjJHWPcc9VFKdIIOAC8h6/UGSLROKXObsSp67kEuBun8szADHfohOpu33VTUW3nIrBzgf2onkAkFiiB6k43v6esWnCv+3793s2NvMAmmeRDQZyEklnJnVKlSuX53sYUSKrzgYzdhNcA630ts2M466TboLoS1VsyHLtz9DyR+4HWQEffz+0A17tHZJXgjvv2oywnIq9nPHIUZIhYuZx8SsTZx3LtWnjoIae7snp1Z0eUPPBXcicqKoo9e/bw5JNPYv8QMuY0hURkaboj0eV15YD0A+nbfef8EymNU52mFiK9s7l3nXTHdcDzwK0u4yKrzftuAZoBTXGWDBgTXGmTUO67z0l0t93mHK+/DhUq5Ph2abMl+/Tpw9atW6lYsSIvvvgiixcvZujQofzyyy+8//77tr2XMY4UVc3NDgn+djrPfHq+6h7gYVd3Vu12+pOkJDDedWAudjK5UlVXuL2hF9gYXAQ4fhxefRWef97Z/eSll5ydUAK0ofLIkSPp1q0bl1xyCZ999hmVK1cOyH2Nya9cLxMQSQBmpBuDqw88j2oL32unVaaacfwtEEEWBn5G9VI3H3czi3KPiEwTkd0isktEPhKR4NZDMaZwYXj6aWcnlOuuc8blArgTyiOPPMIXX3zBjh07uOaaa5g3b15A7mtMAbQEqIZIJUSKAHcBnwXkziLTEfnMd8wA1gGuxy7cJLixOMHG4fSrTvedMyb4KlVyiqpOngw7djiTULp3h/3783zrpk2bsmjRIsqUKUOzZs145513AhCwMRFM5AOciYcXI7Idkc6opgCPAbNxytlMRnV1gJ74CjDUdwwEGqGa9TKE9OG66KJcoRlWjYvIT6paMxfBBpUVPI1w+/Y5O6GMHAlxcU45nttucyap5MHff//NXXfdxezZs+nRowcvv/yy1ZYzBY4ndzLJIzctuD9E5F4RifYd9wJ7gh1YbtgsyghXsiS8+SZ8/z2ULg133OFs9ZXHEjnnnHMOM2bMoHv37rz22mvccsst7Nu3L0BBG5NvRItIkq+hEF4iC3y/JiOyP93hvHZ7GxctuIrAG0B9nJkxC4Huqrol18EHmU0yKQCOH3eWE/Tr57TgXnwRHn88z5NQRo8ezT//+U+qVq3K9OnTqVq1aoACNsbbIrEFZ5stm/xtyxZ49FFnnK5mTWcnlGuuydMt582bR9u2bUlNTWXq1Kk0bdo0QMEa412eSnBOqbbMudyT0vaiNPlbfDxMnw5Tp8Lu3VCvHnTr5ozX5dL111/P4sWLKVu2LC1atKBTp04kJCQQFRVFQkKCbdZsTPCllWtbBvyBsxXYr76fXU+lthaciRz798OzzzpleC680Fkg3rZtrieh7N+/n0aNGrFixenLQGNjY0lKSrKyOyaieKoFl8bZ8eQzVGf5XrcCmqHa083lWbbgRCRKRO7Mc5DGhEKJEk5SW7TISXDt28Mtt8Dmzbm8XQn++uuvM85bRQJjQqbOyeQGoPo5cL3bi7NMcKqairO+IWxEpLKIvCvuagcZA3XqwOLFMGwYzJvn7Gs5eDAczXlRjG2Z1KrbmseZm8Z4kHdmUZ7yJyLPIpKASDwifcjBLH43Y3BfisiTIlJBREqlHbkOFxCRMb6dUVZlON9SRNaJyHrx1RRS1Y3qvnaQMY5ChZzdT9auhebNoVcvuOwymDYNctAtn1lFgvPOOy9QkRrjFSdUNVFVp4c7kHTuBs4DpgGfAOf7zrniJsF1Ah4F5uMM7qUN/uXFOKBl+hPi1BR6E2gFVAfuFpHqeXyOKegqVIBPPnHqzcXEOGvnmjSB5ctdXe6vIoGIsHv3boYNG+a33pwxJkBU96LaHdVavqN7Tqp6Z5vgVLWSnyNPO9NqFjWFfC22UzWFjAmEFi1gxQpnF5TVq52iqp06we+/Z3lZx44dSUpKIj4+HhEhPj6ed955hzvuuIOePXvSqVMnjuai69MY44LIeYi8jMgsRL4+ebiUbYITkVgReVZEknyvq4nILXmJORN+awqJSGnx1Q6SLGoHiUhiWh2jlJSUIIRn8r1CheDhh2H9enjySZg4EapVcyoVHDqU6WUdO3Zk8+bNpKamsnnzZjp16sSUKVPo27cv48aNo2nTpuzatSuEvxFjCoyJwC9AJeAFYDPO5s6uuN1s+Rhwre/1dqB/jkJ0x29NIVXdo6oPq2oVzaL8gqomqWptVa1t+wiaLJUsCUOGwJo10LIl9O0LF1/sJLzUVFe3iIqK4oUXXmDy5MksX76cOnXqsNxlt6cxxrXSqL4LHEd1HqqdgHpuL3aT4Kqo6hDgOICqHsZ/Msqr7UD6qpblgR05uYGItBaRJKvUbFypUsVZID5vHpx/Ptx7L9SvDwsXur5F+/btWbBgAapKw4YNmTrVJvsaE0DHfb/+jsjNiNTCyQ2uuElwx0SkGL4KrSJSBQjGoMMSoJqIVJJc1hSyzZZNrjRqBEuWwLhxsH07NGgAd93lbAPmwlVXXcWSJUuoUaMG7du35/nnnyfVZUvQGJOl/r4q3j2BJ4F3gCdcX62qWR7AjcA8nC1SJuL0gTbO7rps7vkB8DtOdt4OdPadvwlnS5YNQJ9c3Lc1kBQTE6PG5MqBA6p9+6oWK6YaE6Pau7fq/v2uLj18+LDef//9Cmjbtm31wIEDQQ7WmMABjgBJQGvNw9/vATsgWuGJvNzD1VZdIlIap99TgB9U9U/XGTQMbKsuk2fbtjm15yZMgAsugP794cEHIZveAVVl2LBhPPXUU1x55ZV8+umnma6lM8ZLPLpV11xUm+T2crebLV8P3AA0Aa7L7cOCzcbgTMBUqADjxzvbflWtCl27OtUKJk2CLP77EhF69uzJjBkz2LhxI3Xq1GFhDsb0jDGnWYjIG4hch8hVJw+X3NSDewuoitOtCNAB2KCqj+Y65CAa/W7sAAAgAElEQVSzFpwJKFWYMgWef97ZGaVqVWdnlPvugyJFMr1s7dq1tG7dmm3btvHAAw8we/Zstm7dSsWKFRkwYIBt1mw8xbMtuDMpqq5qWLlJcKuBy319tIhIFLBSVS/LaazB5ttDrXVMTEzXI0eOhDscE2lSU51dUf7zH1i2DMqXh6eegi5dIMNuJ2n27t3Lddddx5o1a047bxUJjNd4MsHlkZsuynVA+kGECsDPwQknb9RmUZpgiopytvpasgRmz4bKlaF7d0hIgIED/dagK1WqFAcOHDjjvFUkMCb4Mm3Bich0nKUBJYE6wGLf67rAQlVtFqogc8q6KE3ILFjgJLdZs5xyPY89Bj16QLrNmKOiovzuWSkitpzAeEYktuCySnBZ1txR1XlBiSgPrIvShM3y5U6imzoVihaFxERnO7Dy5UlISGCLnzV1JUuWZM+ePViPg/GCApXgzvigSAng5B5YmoMdnUPNWnAmbNatg0GDnOUFIvCPf/DpJZdwT9++HEq332V0dDQnTpygSZMmTJgwgbi4uDAGbQyIyFHgPWC6eqlkjsi1QALp8g+q77m61MUkk0TgJeAwkIqzFk41jxUFgskSnAm7LVvg5ZfhnXfg+HE2163LQ5s38+XOnSdnUR47doxHH32Us846iwkTJtC8efNwR20KME+24ETGA1WAn4C09TmK6uOuLneR4H4F6nt9cXd6luCMZ+zcCa+9Bm+9BcnJTvHVDh2gTRsoXZo1a9Zw5513snr1anr37s2LL76IbRZuwsGjCW4tUD23hRfdzKLcAGReS8RDbKG38ZwLL3S6LLdsgRdfhF9/hc6dnfMtWlB94UIWz5pFly5dGDhwII0bN2bbtm3Z39eYgmEVcGFuL3bTgquFUzJnEek2WVaXTcRwsBac8SxVZ0LKlCnOsWGDs/1XkyYsrliRDpMmsb9oUcaNG0fr1q3DHa0pQDzagpsL1MSZxX9qk3/VW11d7iLBLQYWACtxxuB899f/5jza0LAEZ/IFVafKeFqy+/VXNCqKJbGxjD1wgPMSE3l2xAiKZLFbijGB4tEE5382v8tZ/G4S3EJVvTbLD3mMJTiT76jCypUwdSqpkycTtW4dqcDys84i/sknKZOYCGXLhjtKE8E8meAARC7AWYsNsBjV3a4vdZHgBgBbgOmc3kXpuWUCtg7ORARVWL2atS+9BFOncmlqKirCHxddxJu7d/P5X39xTtmyPNqzJ206dHDW3aUdNkHF5JInE5zIncDLwDc4M/ivA55C1VVlYTcJbpOf07ZMwJgQ2Lx5M71uvZWLVq7kTuDy7C6Ijj494RUtCsWKnXkuq/f9fd7fZ8qVg7POCsGfggkFjya4FcCNJ1ttIucBX6F6pavLczn70tMswZlIcvz4ccqUKcP+/fu5CKgGFPUd5UqVYvALL8CRI85x+PCpnzM7/H0m7VxOxcfDZZdB9eqnH2efHdg/BBN0Hk1wK1G9It3rKGDFaeeyutxFC+5+f+fV5UrycLAEZyJNSPazVIVjx/wnwYyvDx2CTZtgzRrn+OUXOHpqkhsVKpxKdukTYMmSgYnVBJxHE9zLQA1OL9f2M6rPuLrcRYIbke5lUZzCpz+qarucRxsaluBMpMlsP8vY2FgWL17MZZeFuXpVSsrpCW/1aufXtWtPbxmWK3d6S69uXbjiCqdSgwkrD2/V1RZogDMGNx/Vaa4vzWkXpYiUBMary3UI4WAJzkSaiRMnkpiYeNp+loUKFSIqKorjx4/Trl07nn32WWrUqBHGKP04cQI2bz6V+NKS39q1TisQ4Nxz4brr4PrrnaNmTWcs0YSUJ1tweZSbBFcY+FlVLw1OSHlnCc5EookTJ9KnT5/TqoK3aNGCV199lREjRpCcnMztt99O3759qVmzZrjDzVpqqtPi++47mDfPOTZscN4rUQIaNnSSXePGcNVVNjs0BDyV4EQWoNoQkWScMm0n38HZi7KEq9u46KJMqwsHztZe1YHJqtor51EHly0TMAXV3r17GT58OMOHD2ffvn3ceuutPPfcc9SuXTvcobm3fTvMn38q4a1b55w/6yxo0OBUC692bbDF7wHnqQQXIG4SXPqV5CnAFlXdHtSo8shacKag+vvvvxkxYgSvvvoqf/31FzfddBN9+/albt264Q4t53buPD3hrV7tnC9WDK699lTCq1sXYmLCG2sE8GSCEynl52wyqsddXW7LBIyJPPv37+eNN95g6NCh7N27l+bNm9OvXz82bdp0Rjdnx44dwx2uO3/8Ad9+eyrh/fyzM/MzNhYaNYJmzZzDJq3kikcT3GagAvAXTvfkOcDvwG6gK6rLsrzcRQvuDmAwcL7vAWn14Fz1gYaDJThjHMnJyYwcOZJXXnmFP/74g6ioqNOWFcTGxpKUlJR/klx6e/c6CW/OHPjqK2fiCsB558ENNzjJ7sYboWLF8MaZT3g0wY0CpqE62/e6OdASmAwMRzXLrgk3CW490FpV1wYk4BCwBGfM6Q4ePEj58uX5+++/z3gvPj6ezZs3hz6oQPvtNyfRpR07dzrnq1U71bpr0sSZtWnO4NEEtxTV2n7PifyEapazqdwkuO9UtUHeIw0dS3DGnCmzxeIAR48ejayqBarOkoS0ZPfNN3DggNN1efXVp1p39es7W48Zrya4L4A5wCTfmQ7AjTituCWoXpXl5S4S3HCcgnOfcPpmyx/nPurgsgRnzJkyWywOEBcXxxNPPEFiYiIlSnh29CH3jh+HxYvhyy+dhPfDD84avWLFnCUJ9epBnTrODM0CWrXBowmuDNAPaOg7swB4EdgHVER1fZaXu0hwY/2cVlXtlPNoc05EigNvAceAb1R1YnbXWIIz5kz+FovHxsby2GOPsXTpUr7++mtKlCjBI488wuOPP05cXFwYow2y/fudGZpffQVff+3M0Ewbm4yLO5Xs0o4yZcIbbwh4MsGlETkL1QM5vk5VQ34AY3BmwazKcL4lsA5YD/TynbsPZwwQ4EM394+NjVVjzJkmTJig8fHxKiIaHx+vEyZMOPnekiVL9M4779SoqCgtUqSIdu7cWdeuXRvGaEPowAHVBQtUX3tNtWNH1YsvVnU6Op0jIUG1XTvVQYNU58xR/euvcEcccMBBDUM+yPKAaxXWKGz1vb5S4S2314dlmYCINAIOAO+p6uW+c9HA/+H0r24HlgB3A22Az1X1JxF5X1Xvye7+1oIzJvc2bNjAsGHDGDNmDEeOHKFNmzY89dRTNGiQr4bi827fPli+HJYsgaVLnWPjxlPvV6t2qqV31VVQtarTvZlPlyh4ci9KkUVAO+AzVGv5zq3ClzeyvTwcCQ5ARBKAGekSXH3geVVt4Xvd2/fR7cBfqjpDRCap6l3Z3dsSnDF598cff/DGG2/wxhtvsHfvXq699lqefvppkpOTefbZZ/PnWrq82rMHli1zkl1a4tuebt+LmBinhFClSv6PUqVAJHzxZ8GTXZQii1Cti8jydAluhefrwflJcO2Alqraxff6PqAu8AzwBnAEWKCZjMGJSCKQCFCkSJGrj6Yv3WGMybWDBw8yZswYhg4dypYtWxCR02Zj5uu1dIGwcyf89JPTutu06fTjr79O/2yJEpCQkHkCLB6+/OLRBDcVGIaTA+oBjwO1cdHQgSwSnIj8K6sLVXVYziI94/4JnJ7g2gMtMiS4a1S1Ww7uaXtRGhMkKSkplC1blj///POM9yJmLV2g7dt3ZtJLfxw+fPrnL7gAKleGKlWcX9P/fOGFQe3+9GiCKwMMB5rhbDLyBdAd1T1uLs9qi+60krwXA3WAz3yvWwPzcxVs1rbjbMmSpjywIyc38PUbTy9evHjXQAZmjHHK8+zZ4//vlS1btrB///7IXGKQFyVLOuV//FV3UIXdu09PeBs3Osf8+TBxovOZNEWLnkp6GZNgpUrOkodI4szLuA/VXHcNuFkm8AXQVlWTfa/PBqaoasvcPtR3nwROb8EVwplkcgPwG84kk3tUdXUO7mktOGOCKKu1dOeccw7dunWje/fulC5dOsSRRaCjR2HrVqeMUFriS/t5wwbIOM8gLs5JdrNmwdln+79nFjzagvsG1ca5vtxFgvsFuFJVj/pexwArVPWSXD9U5AOgMVAG2AX0U9V3ReQm4DUgGhijqgNyc3+bZGJMcGS2lq5Xr14sX76cadOmUbx4cR555BF69uzJhRdeGMZoI5iqs/l0xsS3fTt88UWuJrJ4NMENAEoCHwKn/lJX/dHV5S4SXB/gTmAaTl2423Hqwf0ndxEHj7XgjAk+f4VX0yaYrFq1ioEDBzJp0iQKFy5Mly5dePrpp6loGx57nkcT3Fw/ZxXVpq4udzOLUkSuAq7zvZyvqsvdRxh61oIzJrzWr1/PoEGDeO+991BV7r//fnr16kW1atXCHZrJhCcTXB65nZITC+xX1eHAdhGpFMSYjDH5XNWqVXnnnXdYv349Dz/8MO+//z6XXHIJd999NytXrmTixIkkJCQQFRVFQkICEydmuwOfMTnmpouyH1AbuFhVLxKROJxJJp7b1sC6KI3xpp07dzJs2DBGjhzJgQMHiI6O5sSJEyffL/Br6TwgEltwbhLcT0At4Ef1rSQXkZ9VtUYI4ssV66I0xpv27NlDlSpV2Ldv3xnv2Vq68IrEBJfVOrg0x1RVRUTh5O7+npSuBRfuUIwxfpQuXZr9+/f7fW/Lli1s3LiRypUrhzgq4znOfsWZU3W1FttNC+5JoBrOJsgDgU7AB6r6uqtAw8BacMZ4V1Zr6QBatGjBww8/zC233EKhQm7+DW4CwVMtOBF/mz0rcCVQHtVoV7dxOYvyRqA5zlYps1X1yxyEGnKW4IzxrszW0g0aNIi9e/cyevRofvvtN8qVK0fXrl3p0qUL5cqVC2PEBYOnElxGIg2BPsC5wADcVjvIrp4OMNjNOS8dVg/OGG/Lqi7d8ePHddq0adqiRQsFNDo6Wm+//XadPXu2njhxIoxRRza8WQ/uBoVvFOYq3JjT6910Uf6oqldlOOfJSSY2i9KYyLJhwwaSkpIYM2YMf/75J1WqVOGhhx7igQce4Isvvsh0wbnJOU/VgxO5GafFtg/oj+p3ubpNZglORB4B/glUwamwneZsYKHmYQPMYLMuSmMiy9GjR/noo48YNWoU3377LdHRzhCMLTUIHE91UYqk4mzAvwJn7O10qre6uk0WCa4kTn/nQKBXureSVXVvDsMNKUtwxkSu1atXU79+fZKTk894z5Ya5J7HEtz1Wb6vOs/VbVx0UdYDVuvp1QSqq+oid5GGniU4YyJbVFQU/v7uEhFSU1PDEFH+56kEl0akKFAVpxW3AdUcjT252aprJHAg3euDvnOeIyKtRSQpfbeFMSbyZLZ5c2xsrN9F5CafESmEyBCcbsr/AhOAbYgMQaSw29u4SXCi6f6ppKqpuFsgHnKqOl1VE9P6540xkWnAgAHExsaedq5w4cIcPHiQGjVqMHeuv03oTT7yMlAKqITq1Ti7aFUBzgFecXsTNwluo4g8LiKFfUd3YGOuQjbGmADo2LEjSUlJxMfHIyLEx8czduxYvv/+e2JiYmjatClPPPEEhw8fDneoJnduAbriGxoDQHU/8Ahwk9ubuBmDOx94HWiK0w86B+ihqrtzHnNo2BicMQXXwYMHeeaZZ3jzzTe55JJLGD9+PLVr1w53WJ7nqTE4kf9D9aIcv5dBti04Vd2tqnep6vmqeoGq3uPl5GaMKdiKFy/OG2+8wRdffEFycjL16tXjhRde4Pjx4+EOzbi3BpH7zzgrci/wi9ubuGnBXYQzqeQCVb1cRGoAt6pq/5zFGzrWgjPGAPz11188/vjjTJgwgdq1azN+/HguueSScIflSR5rwZUDPgYOA8tweg/rAMWA21H9zc1t3IzBjQZ6A8cBVPVn4K5chGyMMSF17rnnMn78eKZMmcKmTZuoVasWw4cPt6UEXqf6G6p1gReBzcBW4EVUr3Gb3MBdgotV1cUZzqW4DjSEbJmAMcafdu3asWrVKpo1a0aPHj1o1qwZW7dutcriXqf6NaojUH0d1Tk5vdxNF+XnwGM4VbyvEpF2QGdVbZW7iIPPuiiNMf6oKmPGjKFHjx6kpKSQmprKsWPHTr5fkLf78lQXZYC4SXCVgSTgWuAvYBPQUVUzL+gUZpbgjDFZ2bhxI9WrV+fo0aNnvFdQt/sqkAnu5AedSt5Rmn5dgkdZgjPGZMe2+zpdJCa4bMfgRKS0iLwOfAt8IyLDRaR08EMzxpjgyWy7r8zOm/zHzSSTScAfQFugne/nD4MZlDHGBJu/7b4ALr/8cmyiWmRwk+BKqepLqrrJd/TH2Q/MGGPyrYzbfVWsWJEbbriBmTNn0qpVK/bu9XRVMOOCm0kmrwBLgcm+U+2Ay1S1X5BjS3t+ZZzKriVVtZ2ba2wMzhiTW++++y7//Oc/KV++PNOmTaNGjRrhDikkCuQYHPAQ8D5w1HdMAv4lIskisj+rC0VkjIjsFpFVGc63FJF1IrJeRHpldj2Aqm5U1c4u4jTGmDzr3Lkz8+bN48iRI9SvX5/Jkydnf5HxJDd7UZ6tqlGqWth3RPnOna2qJbK5fBzQMv0JEYkG3gRaAdWBu0WkuohcISIzMhzn5/L3ZYwxuVavXj2WLl1KzZo16dChA88884yNy+VDbmZRds7wOlpEXHVPqup8IGNH9jXAel/L7BhOi7CNqq5U1VsyHLapszEmLMqWLcvcuXN5+OGHGTJkCDfddJONy+UzbroobxCRWSJSVkSuAH4Azs7DM8sB29K93u4755dvmcIooJaI9M7ic4kislRElqakeHInMWNMPlOkSBFGjhzJ6NGj+eabb6hTpw4rV64Md1jBEi0iSSLSOtyBBIqbLsp7cEqGrwRm4tSCezIPzxR/j8ni+XtU9WFVraKqA7P4XBLwAvCjiL9HGGNM7nTp0oV58+Zx+PBh6tWrF6njcidUNVFVp4c7kEBx00VZDegOfISzq/N9InLm4hH3tgMV0r0uD+zIw/1OUtXpqpoYHR0diNsZY8xJ9erVY9myZSfH5Xr16mXjch7npotyOvCcqj4EXA/8CizJwzOXANVEpJKIFMEpvfNZHu53klUTMMYEU/pxucGDB1OrVi0qVKhg1Qg8ys06uBKquj/DuWqq+mu2Nxf5AGgMlAF2Af1U9V0RuQl4DYgGxqjqgFzG75etgzPGBFuXLl149913TzuXn6sRROI6uEwTnIg8rapDfD+3V9Up6d77j6r+O0QxuuYbHG0dExPT9ciRI+EOxxgTwRISEtiy5cyiKvm1GkFBS3A/qupVGX/299prrAVnjAm2rKoRnDhxgvw22S0SE1xWY3CSyc/+XhtjTIGSWdUBVaVly5asXbs2xBGZjLJKcJrJz/5ee4JNMjHGhIq/agSxsbHce++9LFq0iCuuuIIePXrw999/hylCk1WCu1JE9otIMlDD93Pa6ytCFF+O2DIBY0yoZKxGEB8fT1JSEuPHj+fXX3+lS5cuvP7661SrVo2kpCRbUhAGrit65yc2BmeM8YKffvqJxx9/nG+//ZaaNWsyfPhwGjVqFO6w/CpoY3D5jnVRGmO8pGbNmsybN48PP/yQPXv2cP3119OhQwe2bt0a7tAKBGvBGWNMCBw6dIghQ4YwePBgRIRnnnmGp556ym9V8XCwFpwxxphciY2N5fnnn2fdunW0bt2a559/nksvvZRu3boRHx9vu6EEQUS14GyhtzEmv5g/fz733nsv27ZtO+18uHZDicQWXEQluDTWRWmMyQ/i4+P9jseFYzcUS3D5hCU4Y0x+kNVuKKmpqSGNJRITnI3BGWNMmGS2G0qhQoX87nNpcsYSnDHGhIm/3VBiYmIoVKgQV199NV9//XWYIosMEZXgbB2cMSY/8bcbyrvvvsuKFSu44IILuPHGGxk2bJjfbkyTPRuDM8YYD0pOTuaBBx7g448/5p577mH06NFBXTNnY3DGGGNC4uyzz2bq1KkMGDCADz74gAYNGrBp06Zwh5WvWIIzxhiPEhH+/e9/M3PmTDZv3kzt2rX56quvwh1WvmEJzhhjPK5Vq1YsWbKEsmXL0qJFC1555RUbl3PBEpwxxuQDVatW5YcffuCOO+7gqaee4p577sHmGmQtohKczaI0xkSys846i8mTJzNw4EA+/PBDrr32WjZu3BjusDwrohKcFTw1xkQ6EaFXr17MmjWLrVu3Urt2bXr16kVCQoJt2JyBLRMwxph8asOGDTRu3Jjt27efdj43GzZH4jIBS3DGGJOPVaxY8YyKBJDzDZstweUTluCMMQVFoDZsjsQEF1FjcMYYU9BktmFzZucLEktwxhiTj/nbsDk2NpYBAwaEKSLv8HyCE5HbRGS0iHwqIs3DHY8xxniJvw2bw1ER3IuCOgYnImOAW4Ddqnp5uvMtgeFANPCOqg5yca9zgVdUtXN2n7UxOGOMyZlIHIMLdoJrBBwA3ktLcCISDfwfcCOwHVgC3I2T7AZmuEUnVd3tu24oMFFVf8zuuZbgjDEmZyIxwRUK5s1Vdb6IJGQ4fQ2wXlU3AojIJKCNqg7Eae2dRkQEGAR87ia5GWOM8TCR24CbgfOBN1H9IliPCscYXDkg/aKN7b5zmekGNAPaicjDmX1IRBJFZKmILE1JSQlMpMYYY04RGYPIbkRWZTjfEpF1iKxHpFeW91D9BNWuwANAh6DFSpBbcJkQP+cy7SdV1deB17O7qaomAUngdFHmOjpjjDGZGQe8Abx38owz7PQm6YedRD4jk2EnfMNOwLO+64ImHAluO1Ah3evywI5A3FhEWgOtY2JiAnE7Y4wx6anOJ5NhJ3zDTviGnchk2Il0w04EedgpHAluCVBNRCoBvwF3AfcE4saqOh2YLiJdRORwHm4VDeS1JEFu7uHmmuw+k9n7OTmf8VwhIBz9vuH4HoL5HWT2nttz4fgevPodZPe53Lzn1f8XAvEduLlPMRFZmu51kq9nLDv+hp3qZvH5tGGnkohURXWUi2fkjqoG7QA+AH4Hjvt+051952/CmUm5AegTzBhyGXdSOO7h5prsPpPZ+zk5n/EcsLSgfA/B/A7c/nlncS7k34NXv4Pc/Dln955X/18IxHcQyPsoJCisSve6vTrLvdJe36cwItR/Tv6OYM+ivDuT87OAWcF8dh5ND9M93FyT3Wcyez8n5wPx+w+EcHwPwfwOMnvPvoPcfT6nf87ZvefV7yFQMQTr9xK0Yae8isjNlk1gichSVa0d7jgKOvsews++A/CNwc0gbfMOkUI4PXI34Aw7LQHuQXV1mCI8yfNbdRlPcNMPb4LPvofwK9jfgcgHwPfAxYhsR6QzqinAY8BsYC0w2QvJDawFZ4wxJkJZC84YY0xEsgRnjDEmIlmCM8YYE5EswZk8EZHKIvKuiEwNdywFiYgUF5H/+molWuGvMLH//r3NElwBJiJjRGS3ZNg4VURaisg6EVkv2Wycqqob1UWNPpO9HH4fdwBT1dm09taQBxvBcvI92H//3mYJrmAbB7RMf0JObZzaCqgO3C0i1UXkChGZkeE4P/QhR7RxuPw+cBbTpm2PFIhtnMwp43D/PRgPC8delMYjNAD1+kzg5OT7wNk9ojzwE/YP1YDK4fewJrTRmZyw/zFMRjmq1ycipUVkFFBLRHoHO7gCKLPv42OgrYiMxBvbSUU6v9+D/ffvbdaCMxnltF7fHiDTQrQmz/x+H6p6EHgw1MEUYJl9D/bfv4dZC85k5NmNUwso+z68wb6HfMgSnMnoZL0+ESmCU6/vszDHVJDZ9+EN9j3kQ5bgCjBJt3GqiGwXkc7qZ+NU9cjGqZHOvg9vsO8hcthmy8YYYyKSteCMMcZEJEtwxhhjIpIlOGOMMRHJEpwxxpiIZAnOGGNMRLIEZ4wxJiJZgjMFjoicEJGf0h1ZlgQKJRGZ6qsxtsgX21YR+SNdrAmZXNdfRF7KcK62iPzs+3mOiJQM/u/AGO+wdXCmwBGRA6p6VoDvWci3GDgv97gM6K+qt6c79wBQW1Ufc3HtNFW9KN25V4A9qjpQRDoDZVR1cF5iNCY/sRacMT4isllEXhCRH0VkpYhc4jtf3FcEc4mILBeRNr7zD4jIFBGZDnwhIlEi8paIrPbVy5slIu1E5AYRmZbuOTeKyMd+QugIfOoizlYi8r0vzg9FpLhvV40jInK17zMCtAcm+S77FLgnL38+xuQ3luBMQVQsQxdlh3Tv/amqVwEjgSd95/oAX6tqHaAJ8LKIFPe9Vx/4h6o2xamynQBcAXTxvQfwNXCpiJzne/0gMNZPXA2AZVkF7isy2wu4wRfnz0B339sf4OyRmHavHaq6CUBV/wTOFpFzsrq/MZHEyuWYguiwqtbM5L20ltUynIQF0By4VUTSEl5RoKLv5y9Vda/v54bAFFVNBXaKyFxwaqqIyHjgXhEZi5P47vfz7LLAH9nEfi1ORemFTiONIsAC33sfAPNE5GmcRPdBhmv/8D3j72yeYUxEsARnzOmO+n49wan/PwRoq6rr0n9QROoCB9OfyuK+Y3EKkx7BSYL+xusO4yTPrAjwP1W9L+MbqrpZRHYA1wG3A1dn+EhR3zOMKRCsi9KY7M0GuvnGtRCRWpl8bgFOle0oEbkAaJz2hqruwKkf9iwwLpPr1wJVs4llIXC9iFT2xVJcRKqle/8D4HVgraruTDspIlFAGU6vSm1MRLMEZwqijGNwg7L5/EtAYeBnEVnle+3PRziFMVcBbwOLgH3p3p8IbFPVNZlcP5N0SdEfVd0FdAY+FJEVOAnvonQfmQxczqnJJWmuARao6oms7m9MJLFlAsYEkIicpaoHRKQ0sBhokNaSEpE3gOWq+m4m1xYD5vquCWgiEpE3cWqYzQvkfY3xMhuDMyawZvhmKhYBXkqX3JbhjNf1zOxCVT0sIv2AcsDWAMe13JKbKWisBWeMMUpHSygAAAAySURBVCYi2RicMcaYiGQJzhhjTESyBGeMMSYiWYIzxhgTkSzBGWOMiUiW4IwxxkSk/wcFqmyC10Hw4AAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot expected number of counts for signal and background\n", "fig, ax1 = plt.subplots()\n", "# ax1.plot( t[\"energy\"], t[\"excess\"],\"o-\", color=\"red\", label=\"signal\")\n", "ax1.plot(\n", " t[\"energy\"], t[\"background\"], \"o-\", color=\"black\", label=\"blackground\"\n", ")\n", "\n", "ax1.loglog()\n", "ax1.set_xlabel(\"Energy ({})\".format(t[\"energy\"].unit))\n", "ax1.set_ylabel(\"Expected number of bkg counts\")\n", "\n", "ax2 = ax1.twinx()\n", "ax2.set_ylabel(\"ON region radius ({})\".format(on_radii.unit), color=\"red\")\n", "ax2.semilogy(t[\"energy\"], on_radii, color=\"red\", label=\"PSF68\")\n", "ax2.tick_params(axis=\"y\", labelcolor=\"red\")\n", "ax2.set_ylim(0.01, 0.5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercises\n", "\n", "* Also compute the sensitivity for a 20 hour observation\n", "* Compare how the sensitivity differs between 5 and 20 hours by plotting the ratio as a function of energy." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.0" }, "nbsphinx": { "orphan": true } }, "nbformat": 4, "nbformat_minor": 2 }