{ "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.15?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`" ] }, { "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": [ "irfs = load_cta_irfs(\n", " \"$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits\"\n", ")\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],\n", " energy_reco[1:],\n", " data=(bkg_data * solid_angles).to_value(\"s-1\"),\n", " unit=\"s-1\",\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.77719e-11142.5531826.68significance
0.02887531.0647e-11220.6074420.96significance
0.04307357.74208e-12180.2282938.04significance
0.06425324.22411e-12111.9851118.13significance
0.09584711.73695e-1275.5441499.549significance
0.1429761.23576e-1250.3005215.216significance
0.2132798.16227e-1331.625280.8063significance
0.318155.57441e-1320.714932.1657significance
0.4745874.1604e-1313.996313.1734significance
0.7079463.25327e-139.720055.43152significance
1.056052.64287e-137.197562.4464significance
1.575322.29178e-136.002561.44036significance
2.349921.90235e-135.14260.878003significance
3.505392.09706e-1350.522997gamma
5.229032.62358e-1350.337574gamma
7.800193.51438e-1350.161236gamma
11.63565.19109e-1350.0682205gamma
17.3576.78149e-1350.0365443gamma
25.89151.05731e-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.77719e-11 142.553 1826.68 significance\n", "0.0288753 1.0647e-11 220.607 4420.96 significance\n", "0.0430735 7.74208e-12 180.228 2938.04 significance\n", "0.0642532 4.22411e-12 111.985 1118.13 significance\n", "0.0958471 1.73695e-12 75.5441 499.549 significance\n", " 0.142976 1.23576e-12 50.3005 215.216 significance\n", " 0.213279 8.16227e-13 31.6252 80.8063 significance\n", " 0.31815 5.57441e-13 20.7149 32.1657 significance\n", " 0.474587 4.1604e-13 13.9963 13.1734 significance\n", " 0.707946 3.25327e-13 9.72005 5.43152 significance\n", " 1.05605 2.64287e-13 7.19756 2.4464 significance\n", " 1.57532 2.29178e-13 6.00256 1.44036 significance\n", " 2.34992 1.90235e-13 5.1426 0.878003 significance\n", " 3.50539 2.09706e-13 5 0.522997 gamma\n", " 5.22903 2.62358e-13 5 0.337574 gamma\n", " 7.80019 3.51438e-13 5 0.161236 gamma\n", " 11.6356 5.19109e-13 5 0.0682205 gamma\n", " 17.357 6.78149e-13 5 0.0365443 gamma\n", " 25.8915 1.05731e-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": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEKCAYAAADq59mMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZhU1bX38e9qpkZAVBQEQRtEAwjSzBFQwSiOOIFBhVxRlIu+OMR5NhG9ajQxapSEhGg0iUHFGBGcBwYxMijiQNSEoCKDiIqIzKz3j12NTdtdXdVdVaeG3+d5zkPVqVOnFn26e/U+e++1zd0RERFJRlHUAYiISO5R8hARkaQpeYiISNKUPEREJGlKHiIikjQlDxERSVrdqAPIhN13391LSkqiDkNEJKfMnz//c3ffo7LXCiJ5lJSUMG/evKjDEBHJKWb2UVWv6baViIgkTclDRESSpuQhIiJJK4g+DxHJfps3b2bp0qVs2LAh6lAKTnFxMa1bt6ZevXoJv0fJozJ77gkrV35/f4sWsGJF5uMRKQBLly6lSZMmlJSUYGZRh1Mw3J3Vq1ezdOlS2rZtm/D78vq2lZkNNrMJa9asSe6NlSWOePtFpNY2bNhAs2bNlDgyzMxo1qxZ0i2+vE4e7j7F3Uc3bdo06lBEJAFKHNGoydc9r5OHiEhtnX322bz33ns1eu+yZcsYOnTo9uennXYaBx54IHfeeSfXX389L7zwQqrCzDj1eYhI7slgv+Qf/vCHGr+3VatWPPbYYwCsWLGC2bNn89FHVc67yylqeYhI7klTv+S6des49thj6dq1K507d2bSpEkMGDBge4WKiRMnsv/++zNgwADOOeccxo4dC8DIkSO54IIL6Nu3L+3atdueMJYsWULnzp0BGDRoEJ999hmlpaXMnDmTkSNHbj9u7ty59O3bl65du9K7d2/Wrl3LkiVLOPjgg+nevTvdu3dn9uzZALzyyisMGDCAoUOH0qFDB4YPH07ZirCVnWfr1q1cdtll9OrViwMPPJDf/e53tfoalVHLozItWlT9V42IpN9FF8GCBTV774ABle8vLYVf/zruW5955hlatWrF1KlTAVizZg3jx48Hwi2ocePG8cYbb9CkSRMOO+wwunbtuv29y5cvZ9asWfzrX//i+OOP3+F2FcCTTz7Jcccdx4LY/2vixIkAbNq0iWHDhjFp0iR69erF119/TcOGDWnevDnPP/88xcXFfPjhh5x22mnbk9ibb77Ju+++S6tWrejXrx+vvvoqvXv3rvQ8EydOpGnTpsydO5eNGzfSr18/Bg0alNTIqsqo5VGZFSvAPWzr1kHjxjBqlIbpiuS5Ll268MILL3DFFVcwc+ZMyg+2mTNnDoceeii77bYb9erV45RTTtnhvSeeeCJFRUV06tSJlUm0gN5//31atmxJr169ANh5552pW7cumzdv5pxzzqFLly6ccsopO/S79O7dm9atW1NUVERpaSlLliyp8jzPPfccDz74IKWlpfTp04fVq1fz4Ycf1ubLBKjlUb2ddoIhQ+Cxx+A3v4Hi4qgjEsl/1bQQiDc66JVXavyx+++/P/Pnz2fatGlcddVVDBo0aPtrZbeGqtKgQYOEjy3P3Ssd7XTnnXfSokUL3nrrLbZt20Zxud895T+rTp06bNmypcrzuDv33HMPRx55ZMIxJUItj0QMHw5r1kCsKSsi+WnZsmXstNNOjBgxgksvvZQ33nhj+2u9e/dm+vTpfPnll2zZsoXJkyen5DM7dOjAsmXLmDt3LgBr165ly5YtrFmzhpYtW1JUVMRDDz3E1q1ba3SeI488kvHjx7N582YAPvjgA9atW1fruNXySMRhh4XRHX/5S2iFiEi00tQv+fbbb3PZZZdRVFREvXr1GD9+PJdeeikAe+21F1dffTV9+vShVatWdOrUiVTMIatfvz6TJk3i/PPPZ/369TRs2JAXXniB8847jyFDhvDoo48ycOBAGjVqVKPznH322SxZsoTu3bvj7uyxxx488cQTtY7bkmle5aqePXt6rdfzuPhiuPfe0O+x666pCUxEtlu0aBEdO3aMOoy4vvnmGxo3bsyWLVs46aSTOOusszjppJOiDislKvv6m9l8d+9Z2fF5fduqxuVJKjN8OGzaFPo+RKQg/exnP6O0tJTOnTvTtm1bTjzxxKhDioxaHolyh06doHlzmD49NYGJyHa50PLIZ2p5pItZaH3MmAEffxx1NCIikVLySMbpp4d///rXaOMQEYmYkkcy2rWDvn3hz38Ot7FERAqUkkeyRoyAd9+FhQujjkREJDJKHsk65RSoWzfM+RARKVBKHsnafXc4+ujQ71HNjE8RSa/ly+HQQ1V2LgpKHjUxfDh8+mkYeSUikRk3DmbNghtvTOU5x9GhQweOOOIITjvtNO644w5+//vf06tXL7p27cqQIUP49ttvgVCK/dxzz2XgwIG0a9eO6dOnc9ZZZ9GxY0dGjhy5/ZyNGzfmiiuuoEePHhx++OHMmTOHAQMG0K5dO5588kmAKkuwZy13z/utR48enlLr1rk3aeJ+1lmpPa9IAXvvvfe2P77wQvdDD616KyoqK3u941ZUVPV7Lryw+hjmzp3rXbt29W+//da//vprb9++vd9+++3++eefbz/mmmuu8bvvvtvd3c844wwfNmyYb9u2zZ944glv0qSJL1y40Ldu3erdu3f3N998093dAZ82bZq7u5944ol+xBFH+KZNm3zBggXetWtXd3dft26dr1+/3t3dP/jgA0/5761qlP/6lwHmeRW/V1XbqiZ22glOPjnMNr/3XlXaFcmw3r1h8WL4/HPYtg2KisId5X33rd15Z82axQknnEDDhg0BGDx4MADvvPMO1157LV999RXffPPNDhVqBw8ejJnRpUsXWrRoQZcuXQA44IADWLJkCaWlpdSvX5+jjjoKCGXfGzRoQL169ejSpQtLliwBYPPmzYwdO5YFCxZQp04dPvjgg9r9Z9JMyaOmRoyAP/0JnnoKKiz6IiK1U11FdoBzz4UJE8Lfbps2hZql991Xu8/1Kobgjxw5kieeeIKuXbvywAMP8Eq5su9l5dGLiop2KJVeVFTEli1bAKhXr972cunljyt/TLwS7NlIfR41NXAgtGypUVciEVm5EsaMgX/+M/ybik7z/v37M2XKFDZs2MA333yzfUXBtWvX0rJlSzZv3sxf0vQzn2wJ9qip5VFTderAaafBPffAF1/AbrtFHZFIQXn88e8e33tvas7Zq1cvjj/+eLp27co+++xDz549adq0KePGjaNPnz7ss88+dOnShbVr16bmA8tJtgR71FQYsTbeeAN69IDf/Q5Gj079+UUKSLYURiwru/7tt99yyCGHMGHCBLp37x51WGmnwojlpLQke2W6dYOOHUO5EhHJC6NHj6a0tJTu3bszZMiQgkgcNZHXt63cfQowpWfPnuek5QPKKu1eey189BHss09aPkZEMuevKnyakLxueWSEKu2KSAFS8qittm2hXz9V2hVJgULog81GNfm6K3mkwogR8N578NZbUUcikrOKi4tZvXq1EkiGuTurV69Oel5JXvd5ZMwpp8D554c5H6WlUUcjkpNat27N0qVLWbVqVdShFJzi4mJat26d1HuUPFKhWTM45pjQ73HrrWEOiIgkpV69erRt2zbqMCRBum2VKsOHw7JlMH161JGIiKSdkkeqDB4MTZpozoeIFAQlj1Rp2DBUZps8GdavjzoaEZG0UvJIpREj4OuvQ6VdEZE8puSRSgMGqNKuiBQEJY9UqlMnzDifNi1U2hURyVMJJw8za2RmGoNaneHDYfNmePTRqCMREUmbKpOHmRWZ2elmNtXMPgP+BSw3s3fN7HYz2y9zYeaQ0lLo1EmjrkQkr8VrebwM7AtcBezp7m3cvTlwMPBP4FYzG5GBGHNLWaXdWbMgtjaxiEi+iZc8Dnf3ce6+0N23le109y/cfbK7DwEmpT/EHKRKuyKS56pMHu6+ueyxme1qZgeYWTszK6rsGCmnpAT691elXRHJW/H6PJqa2dVm9jbhNtXvgEeAj8zsUTMbmKkgc9KIEbBoESxYEHUkIiIpF++21WPAJ8DB7v4Dd+/v7j3dvQ1wK3CCmY3KSJS56JRToF49zfkQkbxkhVA7v2fPnj5v3rzMf/CJJ8KcOfDJJ6q0KyI5x8zmu3vPyl6rtiS7mfUDFrj7utjoqu7AXe7+UYrjzC977gkrV4bHdct9mVu0gBUroolJRCRFEpkkOB741sy6ApcDHwEPpjWqfFCWOBLdLyKSQxJJHls83Ns6gdDiuAtokt6wvhMb4TXRzB6Lt09ERDInkeSx1syuAkYAU2MlSuolcnIz+6OZfWZm71TYf5SZvW9m/zazK+Odw90Xu/uo6vaJiEjmJJI8hgEbgVHuvgLYC7g9wfM/ABxVfkcs+dwLHA10Ak4zs05m1sXMnqqwNU/0P5JT5syJOgIRkVqpNnm4+wp3/5W7z4w9/9jdE+rzcPcZQMXysr2Bf8daD5uAvwEnuPvb7n5che2zJP8/25nZaDObZ2bzVq1aVdPTpEffvjBuHGzZEnUkIiI1EkVJ9r0I80fKLI3tq5SZNTOz3wLdYrfPKt1XkbtPiM1L6bnHHnukMPwEtWhR+f7mzWHYMLj+ejj0UFi8OLNxiYikQLVDddPAKtlX5WQTd18NjKluX9apbjjuscfCeedB165wzz1wxhmhqKKISA6IouWxFGhT7nlrYFkEcUTr9NNh4ULo0QPOPBN+/GNYvTrqqEREEhKvttXOZnaLmT1kZqdXeO2+WnzmXGA/M2trZvWBU4Ena3G+KpnZYDObsGbNmnScvvb23htefBFuuw3+8Q848EB4/vmooxIRqVa8lsf9hFtMk4FTzWyymTWIvfbDRE5uZg8DrwE/MLOlZjbK3bcAY4FngUXAI+7+bo3/B3G4+xR3H920adN0nD416tSByy+H11+Hpk1h0CD46U9hw4aoIxMRqVKVta3MbIG7l5Z7fg1wDHA88Ly7d89MiLUXWW2rZK1fD1dcEfpAOncORRUPPDDqqESkQMWrbRWv5dGgwtodNwMTgBlAs9SGKAA0bAh33w3TpsHnn4fOdLPvb3vuGXWkIlLg4iWPKcBh5Xe4+5+AS4BN6QwqVbK+z6MqRx8dOtOrovpYIhIxlWTPZvGG7hbAdRORaNW2JPsuwP8AJeWPd/cLUhWgiIjklkQmCU4jLEP7NrAtveGIiEguSCR5FLv7xWmPREREckYiM8wfMrNzzKylme1WtqU9shTI2Q7zMlXVxyoqUqe5iEQqkeSxiVCC/TVgfmzLid7nnJgkGM+KFaFjvPy2cCHUrw8jRsDWrVFHKCIFKpHkcTHQ3t1L3L1tbGuX7sCkCl26hEmEL7wAt9wSdTQiUqASSR7vAt+mOxBJwqhRobDiDTfA9OlRRyMiBSiRDvOtwAIze5mwoiCgobqRMoPf/hbmzQtJZMECiGLNEhEpWIkkjydim2STJk3gkUegTx/4yU9CSZOiKCrsi0ghSiR5PAZscPetsH0N8gbx35IdzGwwMLh9+/ZRh5IeXbvCXXfBmDGhrPtVlS6qKCKScon8qfoi0LDc84bAC+kJJ7VyfrRVIkaPDsvaXncdzJoVdTQiUiASSR7F7v5N2ZPY453SF5IkxQwmTIC2beHUU0M1XhGRNEskeawzs+1rd5hZD2B9+kKSpO28c+j/WLUqrIW+TVVkRCS9EkkeFwGPmtlMM5sJTCKsBCjZpFs3uPPO0HF+xx1RRyMiea7aDnN3n2tmHYAfEJal/Ze7b057ZJK8c8+Fl1+Gq6+G/v2hb9+oIxKRPFVly8PM+pc9dvfN7v6Ou79dljjMbGcz65yJICVBZvCHP8A++4T+j9Wro45IRPJUvNtWQ8xstpldb2bHmllvMzvEzM4ys4eAp9hxFFbWyfnCiDXRtGno/1i5EkaO1KJRIpIWcVcSNLNdgaFAP6AloaN8ETDV3XNmXGjOriRYG/fcAxdcEPo/Lrkk6mhEJAfFW0lQy9DmK3cYOhSefBJmzoQf/jDqiEQkx8RLHqpnka/MYOLEMGz3oIPC8/LbnntGHaGI5DAlj3y2yy5Vz/nQYlIiUgvxRlu1zGQgIiKSO+LN8/hjrMP8FeAZYJa7b8lIVCIiktWqTB7ufrSZFQMDgJOAO8zsY0IiecbdP85MiCIikm3izjB39w3EkgWAmbUFjgZ+Y2Z7unvv9IdYc3lfkl1EJCJJdZi7+3/d/T53Px7oX+0bIlYQJdmr06JFcvtFRBJQ49FW7r4plYFImqxYEeZ8lG2jRkGDBmEJWxGRGtJQ3UJz7bVh+O4tt0QdiYjkMCWPQlNSElofv/89fPRR1NGISI6qNnmY2dtmtrDCNtPM7jSzZpkIUlLs6qvDLPObb446EhHJUYm0PJ4GpgLDY9sUYAawAnggbZFJ+rRpA//7v3D//bB4cdTRiEgOSiR59HP3q2Jrebzt7tcAA9z9NqAkveFJ2lx1FdStC+PGRR2JiOSgRJJHYzPrU/bEzHoDjWNPNeM8V7VsCeedBw8+CB9+GHU0IpJjEkkeo4A/mNl/zey/wB+As82sEaAhO7ns8suhuBh+/vOoIxGRHBM3eZhZEdDO3bsApUA3dz/Q3ee6+zp3fyQjUUp6tGgBY8fCX/8KixZFHY2I5JC4ycPdtwFjY4/XuPtXGYkqRQpyGdpkXXYZNGoEP/tZ1JGISA5J5LbV82Z2qZm1MbPdyra0R5YCKk+SgN13hwsvDOuev/121NGISI6odhnaWD9HRe7u7dITUuoV5DK0yfjiC2jbFg4/HCZPjjoaEckStVqG1t3bVrLlTOKQBOy2G1x8MTz+OLz5ZtTRiEgOSGSG+U5mdq2ZTYg938/Mjkt/aJJRF10Eu+4KN9wQdSQikgMS6fO4H9gE9I09XwrclLaIJBpNm8Kll8KUKTBnTtTRiEiWSyR57OvuvwA2A7j7esDSGpVE4/zzoVkztT5EpFqJJI9NZtYQcAAz2xfYmNaoJBpNmoSJg888A7NnRx2NiGSxRJLHDYRlaNuY2V+AF4HL0xqVROf//T9o3hyuvz7qSEQkiyUy2up54GRgJPAw0NPdX0lvWBKZRo3gyivhxRdh+vSooxGRLFVl8jCzkrLH7r7a3ae6+1Pu/nnsdTOz1ukPUTJuzJhQOPH668PStSIiFcRredxuZpPN7H/M7AAza25me5vZYWY2DngV6JihOCWTGjYMC0bNmAEvvRR1NCKSheLOMDezToQFoPoBLYFvgUXANOAxd9+QiSBrSzPMa2DDBthvP9h7b5g1K6w8KCIFJd4M87rx3uju7wHXpCUqyW7FxXDNNXDuufDss3DUUVFHJCJZJJHRVlKozjoL9tlHfR8i8j15nTxUkr2W6teH666DuXNh6tSooxGRLFJtVd18oD6PWti8GTp2hJ13hvnz1fchUkBqVVU3NuLq2NiqglJo6tWDVatCtd2iopA8yrY994w6OhGJSCIJYTxwOvChmd1qZh3SHJNkm6+/rnz/ypWZjUNEskYiM8xfcPfhQHdgCWFlwdlmdqaZ1Ut3gCIikn0SuhVlZs0I5UnOBt4E7iIkk+fTFpmIiGStuPM8AMzscaAD8BAw2N2Xx16aZGbqhRYRKUDVJg/gD+4+rfwOM2vg7hur6oWXAvLf/4b1z0WkoCRy26qyVQNfS3UgksVatKh8vxl06wZ//3tm4xGRyMWrqrunmfUAGppZNzPrHtsGADtlLEKJ3ooVYYZ5xe0//wn1r04+GX76U9i0KepIRSRD4t22OpLQSd4a+FW5/WuBq9MYk+SKtm1D0cQrroBf/zqsPjhpEpSURB2ZiKRZtTPMzWyIu0/OUDxpoRnmGfD446EWlhk88ACccELUEYlILdVohrmZjYg9LDGziytuaYlUctfJJ8Mbb8C++8KJJ8Ill+g2lkgei9dh3ij2b2OgSSWbyI7atYNXX4WxY+FXv4JDDoGPPoo6KhFJg0RuW+3h7qsyFE9a6LZVBB57DEaNgjp14MEH4bjjoo5IRJJUq8KIwGwze87MRpnZrimOTfLV0KGhCm9JCQweDI0a7VhUUcUVRXJaIrWt9gOuBQ4A5pvZU+X6Q0Sq1r59GIF13nnw7beVH6PiiiI5KaHaVu4+x90vBnoDXwB/SmtUkj+Ki+Hee6OOQkRSLJH1PHY2szPM7GlgNrCckERERKRAJVLb6i3gCeBGd1dZEhERSSh5tPNCWKtWouOu5W1Fcky8SYK/jj180sy+t2UoPskXVRVXBBgzBrZuzVwsIgVi+XI49NBQni7V4rU8Hor9e0fqP1YKTmXfve5w3XVw882wejX85S/QoEHmYxPJUzfeGMrP3Xgj3Hdfas9dZfJw9/mxh6Xuflf518zsQmB6akOpnJm1A64Bmrr70Ni+E4FjgebAve7+XCZikRQzg5tugj32gIsugi+/DOXdd9456shEclrDhrBhw3fPx48PW3ExrF+fms9IZKjuGZXsG5nIyc3sj2b2mZm9U2H/UWb2vpn928yujHcOd1/s7qMq7HvC3c+JxTEskVgki114ITz0EMyYAQMHwmefRR2RSE576y3Yfffvnu+0EwwfHtZuS5UqWx5mdhpwOtC2Qh9HE2B1gud/APgN8GC589YB7gWOAJYCc2PnrwPcUuH9Z7l7vN8k18bOJbluxAjYbbcwM71/f3juOZV2F6mBNWvgjDPCnWCzcCd4w4bQoE9lQYd4fR5lczp2B35Zbv9aYGEiJ3f3GWZWUmF3b+Df7r4YwMz+Bpzg7rcACRVAMjMDbgWedvc3qjhmNDAaYO+9907ktBK1Y46BF16AY4+Fvn1DAuncOeqoRHLGl1/CoEGh5dGnD3TvDqNHw4QJofM8leL1eXwEfAQclNqPZC/gk3LPlwJ9qjrYzJoBNwPdzOyqWJI5HzgcaGpm7d39txXf5+4TgAkQCiOmMH5Jp759YeZMOPJIOPhgmDo17BORuFavhiOOgHffDcvrlK9Fmo4iD/FuW81y9/5mthYo/8vXAHf3mvZqVjagv8pf7u6+GhhTYd/dwN01/HzJdp07h9LugwbB4YeHCr3HHBN1VCJZa9Wq8KPy/vvwj3/AUUel/zOr7DB39/6xf5u4+87ltia1SBwQWhptyj1vDSyrxfkkH5WUhDGGHTvC8cfDn/8cdUQiWWnlyjDO5IMPYMqUzCQOSKy21b5m1iD2eICZXWBmu9TiM+cC+5lZWzOrD5wKpGXSoZkNNrMJa9asScfpJd2aN4eXXw6znH7yk7BOuohst3w5DBgQRlFNmxZuW2VKIkN1JwNbzaw9MBFoC/w1kZOb2cPAa8APzGypmY1y9y3AWOBZYBHwiLu/W6Poq+HuU9x9dNOmTdNxesmEnXcO/R5DhsBPf6o1QURiPv00JI5PPoGnnw6tj0xKpLbVNnffYmYnAb9293vM7M1ETu7up1WxfxowLYk4pZAVF8OkSVC3im9XrQkiBebjj+Gww8KUqGefhX79Mh9DIsljc2zOxxnA4Ni+eukLSaQSdepEHYFIVliyJLQyvvgijGb/4Q+jiSOR21ZnEobr3uzu/zWztkBO9F6qz0NE8sl//hO6AL/6Cl58MbrEAWCFUG29Z8+ePm/evKjDkNqKV7a9AL6PpbB9+GFocaxfH+bSduuW/s80s/nu3rOy1xIZbdXPzJ43sw/MbLGZ/dfMFqc+TJFaeOSRqCMQSZt//Su0ODZuhJdeykziqE4ifR4TgZ8C8wEtuiDRadGi8s7xevVg2LAwQ+raa7WwlOSVd9+FH/0oNK5ffjl7KvYkkjzWuPvTaY9EpDpVrWizcWMo4HP99eFPtIkTwwgtkRy2fDkMHhw6yOvVC4mjY8eoo/pOIsnjZTO7HXgc2Fi2s6qChNnEzAYDg9u3bx91KJJODRrAAw9Ahw5w9dVhxtTf/x5/9UKRLHfhhTB/PjRqBLNnw/77Rx3RjqrtMDezlyvZ7e5+WHpCSj11mBeQyZPDbPQ99oCnnoIuXaKOSCQpFRdyKpPKhZwSVasOc3cfWMmWM4lDCsyQIaEq7+bNoRrv1KlRRySSsGnToGXL8LhsalM6FnJKhURGW7Uws4lm9nTseSczG1Xd+0Qi06MHzJkD++0XiiredZeG8kpW+89/Qv/GsceG/o1jjw3fssXF6VnIKRUSmST4AKEOVavY8w+Ai9IVkEhKtG4dWiDHHx/WRz/vvNAaEcki69aFAYKdOsErr8AvfgFvvw3168OYMfDPf4Z/qxorEqVE+jzmunsvM3vT3bvF9i1w99KMRJgC6vMoYNu2wTXXwK23hvGOjz4Ku+4adVRS4NzDt+Kll4bChiNGwG23QatW1b83k2rV5wGsi63m57GT/RDIiXofKk8iFBXBLbfA/ffDjBmhI12VeSVC77wT/o4ZNgyaNQvL1jz0UPYljuokkjwuJqy3sa+ZvQo8SFgGNuupJLtsN3JkqOmwtYp5rqrMK2n21VfhDmppaVhjfPx4mDcvmoq4qVDtPA93f8PMDgV+QFhC9n13181jyT2HHBJ1BFKAtm0L05CuvBI+/xz+93/hpptCqyOXVdnyMLNeZrYnQGwBpx7AzcAvzWy3DMUnIpJzli8PtaiefhoOOghGjQqT/ObPDy2OXE8cEP+21e+ATQBmdghwK+GW1RpgQvpDE8mwjRurP0YkAVdfHbrYjjkmLNz00ENh8F82FDRMlXjJo467fxF7PAyY4O6T3f06QPU+JP906wavvRZ1FJLDGjYM4y8eeOC7fStWwDnn5F+9zrjJw8zK+kR+BLxU7rVEamKJZJ+q6l3tsgt8803ovbzwwvBYJEm/+U2YGV6WKLJ1dngqxEseDwPTzewfwHpgJoCZtUdDdSVXrVgRBtlX3L78MtS+Pu88uPvuUPf6ueeijlZyhDv86lehhdGsWUge2Tw7PBWqTB7ufjNwCWGGeX//bjZhERqqK/moSZPwp+PMmeEn/8gjwxDfL76o9q1SuLZsgbFj4ZJL4OSTw9Kw2T47PBW0DK1IZTZsCOMpb7sNdtstJJWhQ/PvxrXUytq1cOqpoaDhZZeFQgZFicyeyxG1nWEuUniKi0PymGAo0M8AAAuWSURBVDcP2rSBH/8YTjoJPv006sgkS3z6aZg69Oyz8NvfhrpU+ZQ4qlNA/1WRGujaNdx/+MUvwm+JTp2gaVOVOClwb70FffqEarhPPRUm/hUaJQ+R6tStG+5JvP02dO8OX39d+XEqcVIQpk2D/v3D3wuzZsFRR0UdUTSUPEQS1b49vPRS9cdJ3ho/Pqy7sd9+8PrrcOCBUUcUHSUPkWSow7wgbdsWyqefd16YNT5jRu5VwU21vE4emuchGffoo1q1MM98+20YaPfLX8L558MTT0DjxlFHFb28Th6a5yEZ9+Mfh4H+06dHHYmkwMqVMHBgSBh33RXmj5atLV7o8jp5iKRFVSVOWrSAP/4xjOEcMACOOy6s/CM5pawi7vTpYUTVO++E5HHBBVFHll2UPESSVVWJkxUr4Mwz4cMPw2yxWbPCUN+zzoKlS6OOWhI0blwoMjBoUCi0PGMGHH981FFlH80wF0mX1avh//4vzE4vKgoFF6+8MhRhlKzTsGEoLFBRcTGsX5/5eLKBZpiLRKFZs9DL+v77MGRIKHWy775w551hQqEmGmaVl17a8Y5kw4b5WxE3FZQ8RNKtpAT+/Gd44w3o0QMuvrjqCYWaaJhxW7aEHP+jH4UamGUVcTduzN+KuKmgdTlEMqVbt1Dm/fnnww11idzChWGJ2HnzQr/Gpk3Qrh2MHg0TJoTOc6mckodIph1xRNQRFLyNG0Pdy1tvhV13hUmT4JRTdpwDeu+90cWXC3TbSiTb3HST1hBJo1dfhdLS8GU+/XRYtChMz1HxgOQoeYhkm+uug733DqOzliyJOpq8sXZtmCF+8MFh1vgzz8Cf/hTGNUjy8jp5qDyJZK14Ew0XLgyjs+67LxRjHD4cFizIbHx55umn4YADwq2o888PKw4feWTUUeW2vE4eKk8iWSveRMMuXcKfxIsXw0UXwZNPhs72QYNCZ3sBzM1Klc8/h5/8JBQzbNw43LK66y7VpkqFvE4eIjmtTRu44w745BO45ZawnsigQWFNkV120TyRSpSVFlm+HB5+GDp2hL/9Da6/Ht58Ew46KOoI84eSh0i222WXMDN9yRKYODFMg67qVmyBzxMpKy1y0EGhM7xduzC95uc/hwYNoo4uv6g8iUiu2bYtfmnXAviZrkilRdJD5UlE8klRNT+2J58MDz5YMMN9160LqwSXb1kUF6u0SLopeYjkm7lz4YwzoHlzOPzwUJgxD6v6fvNNKBdWUhJuV+2++3elRTZtUmmRdFPyEMk3H38Mc+bA5ZfDsmVhbGqbNtCrV6jyu2hRuLWVo8UZ164N4wdKSkJXUI8eMHs29O4N554L//wnjBkTBq5J+qjPQyQX7bln5Z3jLVp8/7fm+++H1Yz+/nd4/fWwb//94YMPqj5/Fv5eWLMG7rknFCX+4osw/Pb668OCTZIe6vMQyTfx5olU9IMfwBVXhD/Jly4Nkw/32SfzMdfQV1/BjTeGlsZ110HfvqFhNXWqEkeUlDxECslee4V7O889F/+4yoYupVnZHI2y/Pfll3DDDSFp3HBDeG3+fJgyJdyBk2ipqq6IfF/r1mFJ3TFjwgJWGTBuXFi59+qrQ467+274+ms46aRwe6q0NCNhSILU5yFSqOKVkR06NPSRbN0aikCdey4ceyzUTf3fm1XN0SgqCrPCDzww5R8pCVKfh4h8X7zijI8+GkZt/fzn8M47cOKJYbr2TTelfBjT4sVhNnjZvMc6dUKe+vRTJY5spuQhUqiq63Rv1SrcL1qyBB5/HDp0CD3WbdrAsGHwyispGe7bsmWYk7FtW5jo5x4q0mf5iOGCl9fJQyXZRVKgbt3Q8fDcc2F47wUXhOq+AwembC32lSvDnbHXX9ccjVyhPg8RSd769WHt1jPPrPqYAvjdku/U5yEiqdWwIYwcGXUUEiElDxERSZqSh4iIJE3JQ0RqLt5wX8lrmmEuIjWnYVEFSy0PERFJmpKHiIgkTclDRESSpuQhIiJJU/IQEZGkFUR5EjNbBXxUy9M0BWpbJCvZcyRyfHXHxHu9stcS3bc78Hk1saVatl6DRI6r6vVk9lfcl6vXoCbnSefPQm2uAeTudUjkHPu4+x6VvuLu2hLYgAmZPkcix1d3TLzXK3stiX3zdA1qfx2S2V9xX65eg3RdhyiuQS5fh9qeQ7etEjclgnMkcnx1x8R7vbLXEt0XhWy9BokcV9XryezPhuuQqhiy6Wch164BRPOzsIOCuG0lqWdm87yKapuSGboG2aFQr4NaHlJTE6IOQHQNskRBXge1PEREJGlqeYiISNKUPEREJGlKHiIikjQlD0k5M2tnZhPN7LGoYykkZtbIzP5kZr83s+FRx1OICul7X8lDdmBmfzSzz8zsnQr7jzKz983s32Z2ZbxzuPtidx+V3kgLQ5LX42TgMXc/Bzg+48HmqWSuQSF97yt5SEUPAEeV32FmdYB7gaOBTsBpZtbJzLqY2VMVtuaZDzmvPUCC1wNoDXwSO2xrBmPMdw+Q+DUoGFpJUHbg7jPMrKTC7t7Av919MYCZ/Q04wd1vAY7LbISFJZnrASwlJJAF6A/DlEnyGryX2eiio28wScRefPcXLYRfUntVdbCZNTOz3wLdzOyqdAdXgKq6Ho8DQ8xsPNlTRiNfVXoNCul7Xy0PSYRVsq/K2aXuvhoYk75wCl6l18Pd1wFnZjqYAlXVNSiY7321PCQRS4E25Z63BpZFFIvoemSDgr8GSh6SiLnAfmbW1szqA6cCT0YcUyHT9YhewV8DJQ/ZgZk9DLwG/MDMlprZKHffAowFngUWAY+4+7tRxlkodD2ip2tQORVGFBGRpKnlISIiSVPyEBGRpCl5iIhI0pQ8REQkaUoeIiKSNCUPERFJmpKHFDwz22pmC8ptcUvOZ5KZPRZbI+L1WGwfm9mqcrGWVPG+m8xsXIV9Pc1sYezxi2bWNP3/A8lXmuchBc/MvnH3xik+Z93YRLLanOMA4CZ3P6ncvpFAT3cfm8B7/+7u+5fbdwew2t1vMbNRwO7uflttYpTCpZaHSBXMbImZ/dzM3jCzt82sQ2x/o9gCQXPN7E0zOyG2f6SZPWpmU4DnzKzIzO4zs3dja51MM7OhZvYjM/t7uc85wswerySE4cA/EojzaDN7LRbnJDNrFJvtvMHMesSOMeAU4G+xt/0DOL02Xx8pbEoeItCwwm2rYeVe+9zduwPjgUtj+64BXnL3XsBA4HYzaxR77SDgDHc/jLCyXwnQBTg79hrAS0BHM9sj9vxM4P5K4uoHzI8XeGzxrSuBH8XiXAhcGHv5YULNpbJzLXP3/wK4++dAEzPbJd75RaqikuwisN7dS6t4raxFMJ+QDAAGAcebWVkyKQb2jj1+3t2/iD3uDzzq7tuAFWb2MoS63Wb2EDDCzO4nJJX/qeSzWwKrqom9L2Elu9mhcUF9YFbstYeB6WZ2OSGJPFzhvatin/FVNZ8h8j1KHiLxbYz9u5Xvfl4MGOLu75c/0Mz6AOvK74pz3vsJCzZtICSYyvpH1hMSUzwGPOPuP6n4grsvMbNlwMHASUCPCocUxz5DJGm6bSWSvGeB82P9CJhZtyqOm0VY2a/IzFoAA8pecPdlhPUfriWskV2ZRUD7amKZDRxqZu1isTQys/3Kvf4wcDewyN1XlO00syJgd3ZcDU8kYUoeIt/v87i1muPHAfWAhWb2Tux5ZSYTFg16B/gd8DqwptzrfwE+cfeq1r2eSrmEUxl3XwmMAiaZ2VuEZLJ/uUMeATrzXUd5md7ALHffGu/8IlXRUF2RNDKzxu7+jZk1A+YA/cpaAGb2G+BNd59YxXsbAi/H3pPSX/Jmdi9hDYrpqTyvFA71eYik11OxEU31gXHlEsd8Qv/IJVW90d3Xm9kNwF7AxymO600lDqkNtTxERCRp6vMQEZGkKXmIiEjSlDxERCRpSh4iIpI0JQ8REUmakoeIiCTt/wPEStW2vNVTOwAAAABJRU5ErkJggg==\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(f\"Energy ({t['energy'].unit})\")\n", "plt.ylabel(f\"Sensitivity ({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": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEKCAYAAABgyEDNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZzN9f7A8dd7BsMobmixzkhaVCJkTbYaKrmFVMoty1S3RGnhqrRwLbd0paKR5caUrRRRKolU1iRbfmWXUFTGbsz798f3DGOcmfnOzFm+c+b9fDy+D3O+53y/3zenvH3Wt6gqxhhjTKSJCncAxhhjTDBYgjPGGBORLMEZY4yJSJbgjDHGRCRLcMYYYyKSJThjjDERqUi4AwgkEWkLtAV6xMbGhjscY4wpMA4dOqTAW8AsVZ0V7ngCQSJxHVzJkiX14MGD4Q7DGGMKDBE5pKolwx1HIEVkCy4mJibcoRhjjAkza8EZY4yJyBacTTIxxhgTkayL0hhjTESKqBacqs5S1cTo6OiwPD85OZn4+HiioqKIj48nOTk5LHEYY4yJsBZcOCUnJ5OYmMihQ4cA2Lp1K4mJiQB07tw5nKEZY0yhFFGTTDJ0UfY4cuRISJ8dHx/P1q1bzzgfFxfHli1bQhqLMcbkViROMomoBJcuHLMoo6Ki8PdnKSKkpaWFNBZjjMmtSExwETUGFy5Hjx4lu51THnroIZYvX+43ARpjjAmOiEpwItJWRJJOnDgRsmfu2rWL5s2bc/DgQYoWLXraezExMTRo0IBx48ZRr149rrrqKl555RV+++23kMVnjDGFVUQluFDPovzuu++oV68eq1atYtq0aYwfP564uDhEhLi4OMaOHcs333zDr7/+yqhRoyhRogSPPfYYFSpU4LbbbmPWrFmkpqaGJFZjjMlBtIgk+eYyRAZVjbgjNjZWg23y5MlaokQJrVy5sq5cudL1dWvWrNE+ffroeeedp4Cef/75+sQTT+i6detUVXXSpEkaFxenIqJxcXE6adKkYP0WjDHmJOCgeuDv70AeNskkl9LS0hgwYAADBw6kUaNGvP/++5x//vm5vs/x48eZM2cO48ePZ/bs2aSmplKtWjW2b9/OsWPHTn4uNjaWpKQkW2pgjAmqSJxkYgkuFw4cOECXLl2YMWMGXbt25Y033iAQu6bs3r2bSZMm0a9fP44fP37G+7bUwBgTbJbgPC6Y6+C2bNnCLbfcwtq1a3n55Zfp1asXIhLQZ9hSA2NMuERigrNJJi4sWLCAevXqsX37dj7++GN69+4d8OQGUKVKFb/nK1euHPBnGWNMpIuoBBcMSUlJtGrVirJly7JkyRJuuOGGoD1r0KBBftfTlSpVipSUlKA91xhjIpEluCwcP36cnj17cv/999OqVSsWL17MxRdfHNRndu7cmaSkpNOWGvTo0YP169fTvHlzdu/eHdTnG2NMJImoMbh0eZlkkpycTP/+/dm2bRsVK1akVKlSrFu3jscee4xhw4YRrgoFALNnz6Zjx45UqFCBuXPnUq1atbDFYoyJTJE4BmcJjjMrAaRLTEzkzTffDHR4ebJkyRJuuukmoqOjmTNnDnXq1Al3SMaYCGIJzuPyOouyoFQC2LBhAwkJCezdu5f333+f66+/PtwhGWMihCW4AiK3LbiCND1/586d3Hjjjaxdu5YJEybYAnBjTEBEYoKzSSacmp5fNIvzXlKhQgUWLFhAkyZNuPvuu3n55ZfDHZIxxniSJThOTc//HPgUuBMoU6IEgwYNCnNk/pUuXZpPPvmEjh078vjjj9OnTx/PtTSNMSbcioQ7AC/o3LkzpKWxvGdPbv3rL94BjolQ7KuvoFo1qF8fgrCwOz9iYmKYPHkyF1xwAcOHD+fXX39lwoQJFCtWLNyhGWOMJ9gYXGZpabBgAYwfD9Onw+HDcNllcO+9cM89UL58QGPNL1Vl6NCh9OvXj1atWvHee+9RqlSpcIdljClgbAyuMIiKgubN4e23YdcuGDMGzjkHnnoKKleGm2+G996DDDv+h5OI0LdvX8aPH8/8+fNp1qwZb7zxBvHx8URFRREfH09ycnK4wzTGeF/E1YPzfAtORC4DegHlgHmqOiqna4JSTWDDBpgwwUl8O3dC2bLQuTPcdx/UqhXYZ+XRxx9/TLt27UhNTT1tVqiV3DHG5CQSW3A5JjgRaQx8r6oHReRu4GpghKqeuXDM7UNFxgE3A3tU9YoM51sDI4Bo4C1VHZLhvShgjKp2y+n+wawHR2oqfPaZ04X54YdOS65WLSfR3XUXlCsXnOe6VL58eXbt2nXGea+t6TPGeEskJjg3XZSjgEMichXwJLAVeDufz50AtM54QkSigdeBNkAN4E4RqeF77xZgETAvn8/NvyJFoE0bmDrVacmNHOl0a/bqBRUqQPv28PXXYQsvq/0qt23bFuJIjDEmvNwkuFRfOfN2OC23EcDZ+Xmoqi4E9mU6fQ3ws6puUtVjwGTfM1HVmaraCPBWH1vZsvDww7BiBaxaBQ89BAsXQpMm8Pe/w/r1IQ8pq7V7FStWDHEkxhgTXm4SXIqI9APuBmb7WlqZ10QHQkVge4bXO4CKItJMRF4VkTeBOVldLCKJIrJcRJanpqYGIbwc1KwJr7wCW7bAwIHwxRdwxRVw//3w668hCyOrkjsHDhxg0aJFIYvDGGPCzU2C6wQcBbqp6i6cRPSfIMTib6GZquqXqvqIqt6vqq9ndbGqJgHPA98FoxipayVLQv/+sHGj06IbPx4uugieeQb27w/64/2V3Bk4cCBly5alWbNmDB061BaFG2MKBTeTTIaq6lM5ncv1g0XigY/SJ5mISEPgOVVN8L3uB6Cqg3N776BOMsmtjRudhDdlCpx7Ljz7LCQmQogXZO/fv58ePXowdepU2rRpw9tvv025ME+IMcZ4R2GdZOJvy/o2gQ4EWAZUF5GqIlIMuAOYmZsbiEhbEUk6ceJEEMLLo2rVYPJkWLoULr8cevaEGjWcSSohXKJRqlQpJk+ezBtvvMG8efOoVauWdVkaYyJalglORB4UkdXAJSLyQ4ZjM7A6Pw8VkXeBb3333iEi3VQ1FXgYmAusB6aq6trc3FdVZ6lqYjiLk2apXj1nXG72bChRAjp1ggYNnF1TQkREePDBB1m8eDHFixenWbNmDBkyxLosjTERKcsuShEpDZwDDAb6ZngrRVUzz4D0hLzWgwu5EyecBePPPgs7dji7owwZ4rTwQiRjl2Xr1q2ZOHGidVkaU4hFYhelq51MfDMnzyfD5syq6tmFVZ4ag8vO4cPw6qsweDCkpDj7XT7/PFSqFJLHqyqjR4+md+/enHvuuUyePJkmTZqE5NnGGG+JxASX4xiciDwM7AY+A2b7jo+CHFfhUKKEs8flxo3OQvFJk6B6dfjXv+DAgaA/3rosjTGRzM0kk97AJap6uape6TtqBjuwvPDkJBM3ypaF4cOd/S7bt3dadC1awL7Q9ATXrl2b7777jvbt29OvXz9uuukmRo8ebRs2G2MKNDfLBOYD1/smgRQIBaaLMiszZ0LHjk6Zns8+c5YXhEB6l2XPnj1JS0uzDZuNKUQisYvSTYIbC1yC0zV5NP28qg4Pbmh5V+ATHMCnnzrbfVWtCp9/HtI6dBUqVOBXP7uv2IbNxkSuSExwbroot+GMvxXD2YMy/fCcAttF6c8NN8DHH8PWrdC0KWzfnvM1AeKvGgHYhs3GmILF8/Xg8iIiWnDpvv0WWreGMmVg3jy48MKgPzI+Pp6tW8+shlShQgV++eWXoD/fGBN6hbIFJyLzReSLzEcogjNAw4bOAvH9+52W3IYNQX9kVhs27927l/fffz/ozzfGmEBw00X5OPCE73gG+B5YHsyg8iqiuigzqlMH5s+H48fhuutgzZqgPs7fhs0jRozgqquuon379jz11FOEpWKDMcbkQp66KEVkgapeF4R4AiKiuigz+vFHaNkSjh51JqFcfXVIH3/06FF69erFm2++SYsWLZg8eTLnhmiGpzEmuAprF2WZDEc5EUkALghBbCazSy91CqqWLOmsk1uyJKSPj4mJYfTo0YwbN46vv/6aq6++mqVLl4Y0BmOMcctNF+UKnC7JFTgbJPcBugUzKJONatXgq6+gXDlo1cpJeCF233338c033xAdHc21115LUlISkThZyZhCJlpEknx7+kYEm0VZUO3c6XRXbt3qLAxv1SrkIezdu5fOnTszd+5c7rvvPl5//XVKlCgR8jiMMflXWLsoi4rIIyIy3Xc8LCJFQxFcbkXsJBN/KlRwSu1Ur+5UI5g9O+QhlC1bltmzZ/PMM88wfvx4mjRpYgvBjTGe4WYnk7eAosD/fKfuAU6oavcgx5ZnhaIFl27fPkhIgFWrnMKqt90WljBmzZrFPffcQ3R0NO+88w4JCQlhicMYkzeebMGJnAc0BioAh4E1wHJUXe0I72YMrp6q/kNVv/Ad9wH18hywCawyZZytvOrVg9tvh3feCUsYbdu2ZdmyZVSoUIE2bdrQsWNH4uLibLNmY0zuiTRHZC7OFpFtgPJADeBpYDUizyNSKsfbuGjBfQd0VNWNvtcXAtNVNbRz1HOhULXg0h04AG3bOt2Wb70FXbuGJYyDBw9y/fXX8+2335523jZrNsbbPNWCE/kPMBJ/dUdFigA3A9GovpftbVwkuJbAeGATIEAccJ+qzs9b5MFXKBMcOAVUb70V5s6F116Dhx4KSxhxcXF+9620zZqN8S5PJbgAcVvROwanooAAP6rq0RwuCatCm+DAWQTeqRN8+CH07+9UCI+ODmkIUVFRfpcNiIgVUzXGozyZ4EQe83P2L2AFqt/ndLmbWZQPASVU9QdVXQXEisg/cx9p8BWqWZRZiYmBadOgWzcYNAhuuilkhVPTValSxe952/XEGJNLdYEHgIq+IxFoBoxB5MmcLnYzyaSHqv6Z/kJV/wB65CnUIFPVWaqaGB3iFovnFC0KY8bAm286GzXXrQvf5/iPnYDxt1mziLB3716mTZsWsjiMMQVeWeBqVPug2gcn4Z0LNAXuzeliNwkuSkQk/YWIROPUhjNeJgKJic6uJ8eOOVUJJk4MyaP9bdY8evRoGjRoQKdOnRg1alRI4jDGFHhVgGMZXh8H4lA9TIYC3FlxM8nkP0A8MBpQnObidnWyqScV6jE4f3bvhjvugC+/hIcfhpdfhmKh/zfKoUOH6NSpEx999BEDBgxgwIABZPi3kzEmjDw6BvcMcCvwoe9MW2Am8DKQhGq207LdJLgonH7PVjiTTD4F3lJVzw50WYLzIzUV+vZ1klujRs44XYUKYQgjlR49ejBhwgQefPBBRo4cSaHvUjbGAzyZ4ABE6gBNcPLPIlRdl2uzvSgLmylTnAkoZ5/tJLkmTUIegqrSt29fhg0bRseOHZk4cSIxMTEhj8MYc4qHE1wToDqq4xE5FzgL1c1uLnUzBmciSadOsHixk+CaN4eRIyHE/8gREYYOHcpLL73EtGnTuOmmm0hJSQlpDMaYAkBkAPAU0M93pigwye3lluAKoyuugGXL4MYb4ZFHoEsXOHQo5GH06dOH//3vf3z55Zc0a9aMPXv2hDwGY4yn3QrcAjhdcqo7gbPdXuz5BCcifxeRMSLyoYjcEO54Ikbp0jBjBrz4IiQnO+NymzaFPIwuXbowc+ZM1q9fT+PGjdm82VXPgzGmcDjm2zXC6WYSyVUXqptJJrNO3vyUv3CKoL6pqkdy80DfPcfh7CW2R1WvyHC+NTACiMaZyDIkw3vnAC+pao7FVm0MLpc+/hg6d3a6Kt95B9q0CXkI3377LTfddBPFixfnk08+oWbNmiGPwZjCzJNjcCKPA9WB64HBQFfgHVRHurncTQtuE3AAGOM79gO7gYt9r/NiAtA64wnf+rrXcXaOrgHcKSI1Mnzkad/7JtDatIHlyyEuztn55MUXIcRbajVs2JCvvvqKqKgomjZtyldffRXS5xtjPEj1JWA68B7OdpHPuk1uvus12wNYmNU5YG1O12dz33hgTYbXDYG5GV738x0CDAVa5XC/RJxW5fJixYqpyYODB1XvvlsVVNu2Vf3jj5CHsHXrVr3kkku0ePHi+thjj2lcXJyKiMbFxemkSZNCHo8xhQVwUPP497lXDzddlOuBBPWVLRCRKr5EdJmIrFTV2q6z6en3jQc+Ul8XpYh0AFqrr5CqiNwD1Af+D/gHsAz4XlVH53Rv66LMB1V4/XV49FGIj4cJE6Bx45CG8Pvvv1O/fn02ZRoTtJI7xgSPp7ooRVI4c2jsFNUca8GBuy7KPsAiEZkvIl8CXwGPizPY979sr8wdf1taqKq+qqp1VPWBnJKbbbYcACLObifz5zuVCZo0gX/8w9kNJUTKlSvH8ePHzzh/6NAh+vfvH7I4jDFhonq2L4n9F+iLs9FyJZwlAwPd3iY35XIuxVcuByfx5Ktkjp8WXEPgOVVN8L3uh/Ogwbm9t7XgAuTgQaciwUsvQYkS8MILTo25IkWC/mgruWNMaHmqBZdOZAmq9XM8lwU35XLGqepRVV2lTv2daGBOnoLN3jKguohUFZFiwB04e465Zi24ACtZEv79b1izxtmsuXdvqF0bFi4M+qOzKrlTqVKloD/bGOMZJxDpjEg0IlGIdAZc/wXvpovyFxEZBSen6n9GLlaS+yMi7wLfApeIyA4R6aaqqcDDwFxgPTBVVdfm5r5q5XKC4+KLnaUEM2ZASgpcd52zrGDnzqA90l/JHXAmRfmrFm6MiUh3AbfjzNzfDXT0nXPFbRflUKA0UAcYoqrv5SnUIBORtkDbmJiYHkeO5Hp5nnHj0CEYMgSGDXPqzj33nLMbStGiAX9UcnIy/fv3Z9u2bVSpUoU777yTUaNGUaJECWbNmkXdunUD/kxjCitPdlHmU5YJTkRuy/gSeAZYCnwCoKrvBz26PLIxuBDYuBF69YLZs+Gyy+C116BFi6A/dt26ddx0003s3r2bd955h7///e9Bf6YxhYGnEpzI08AbqO7L4v0WQCyqH2V7m2wS3PhsrlNV7eoy1JCxFlwYfPSRk+g2bYLbb3fK8QR5nGz37t20a9eOpUuX8tJLL/Hoo49aXTlj8sljCa4d8CRwBPgO+A0ojrOrSS3gc+DfqP6W7W3cdFEWNNaCC7EjR5wuy8GDIToannnGWUcXxKKqhw8fpkuXLkyfPp0HHniAkSNHUiQEszuNiVSeSnDpRKoDjYHywGGc+RkLfRW9c77cxULvC3H2h2yAs/DuW6C3uqzHEw6W4MJk82YnsX34oTMxZeRIuCF4+2OnpaXRv39/hgwZQkJCAlOnTqVUKVfrP40xmXgyweWTm1mU7wBTcTJoBWAaMDmYQeWVLRMIs6pV4YMPYM4cZy/LhARo3x62bw/K46Kiohg8eDBjxoxh3rx5NGnSxGZYGmNOctOCW6KZFtWJyGJVbRDUyPLBWnAecPSoMx43cKDTbTlwoLNDSpCWcHz++ed06NDBZlgak0eFqgUnImVEpAwwX0T6iki8iMSJyJPA7NCFaAqkmBj4179g7Vpnu6/evaF+fVixIiiPa9WqFd988w3FixenadOmfPDBB0F5jjGm4Miui3IFzu78nYD7gfnAl8CDwH1BjywPrIvSg6pWdbosp0yBX36Ba65xxulSUgL+qBo1arB48WJq1qzJbbfdxvDhw/1u92WMKSBEhiFSCpGiiMxD5HdE7nZ9eST+BWBdlB71559Oq270aKhY0Vk7165dwB+TcYZly5Yt+emnn9i+fTtVqlRh0KBBVo3AGD882UUp8j2qtRC5Ffg78CgwH9Wr3FzuZpKJMYHxt7/BG2/A11/DOefA3/8Ot94a8EkoJUqUYMqUKdx8883MmzePbdu2oaps3bqVxMREkpOTA/o8Y0zQpG+RdCPwbpYLv7NgCc6EXsOGzljc0KEwdy7UqAEjRkAAu5ajoqJYvXr1Geet5I4xYSZyISJjEZnu4tOzEPkRqAvMQ+RcnMXf7h6VzU4mjVX1axGJyW9pnFCzLsoCZPNm+Oc/4ZNPoE4dePNN59cAsJI7xrjnqotSZBxwM7AHX6kz3/nWOOulo4G3UB3i4oHTUe3g4nPnAPtRPYFILFAK1V05Xkf2LbhXfb9+6+ZGXmCTTAqgIE5CyarkTpkyZfJ9b2MKqQlA69POiEQDrwNtgBrAnYjUQORKRD7KdJyXq6eJdAHaAp19P3cAXO8ekV2CO+7bj7KiiLya+chVkCFi5XIKKBFnH8v16+H++53uyho1nB1R8sFfyZ2oqCj27t3L448/jv1DyJjTFBGR5RmOxDM+oboQyDwOdg3wM6qbUD2GsxFIO1RXo3pzpmNPLmOql+G4FngOuMXtxdkluJtxarMdwVkykPkwJrACPAmlc+fOJCUlERcXh4gQFxfH+PHjeeihh3j55Zdp164d+/fvD/BvwpgCK1VV62Y4klxeVxHI+D/pDt85/0TKIjIaqI1Iv2zvrNozw9EDqA243uTWzU4mV6nqKrc39AIbg4sAx4/DK6849eaio+HFF52dUAK0ofKoUaPo2bMnl156KTNnzuTCCy8MyH2NKahcLxMQiQc+OjkGJ9IRSEC1u+/1PcA1qPYMQpBFgR9QvczNx93MotwrIjNEZI+I7BaR90QkuPVQjClaFJ580tkJ5dprnXG5AO6E8uCDD/Lpp5+yc+dOrrnmGhYsWBCQ+xpTCO0AKmd4XQnYGZA7i8xCZKbv+AjYALgeu3CT4MYDM3E2Wq4IzPKdMyb4qlZ1iqpOnQo7dzqTUHr1ggB0LbZo0YIlS5ZQrlw5WrVqxVtvvRWAgI0pdJYB1RGpikgx4A6cnBEILwEv+47BQFNU+7q92E2CO09Vx6tqqu+YAJybp1CDzGZRRigR6NgRfvwRHnjAKcNTowbMmAH53ImnevXqLF68mJYtW9KjRw8effRRUlNTAxS4MRFG5F2cmfWXILIDkW6opgIP48zZWA9MRXVtQJ6nuiDD8TWqO3IVrosxuM9xpoa+6zt1J3CfqrbMS7yhYGNwEW7JEkhMhB9+gLZtnS2/slgS4FZqaiqPP/44I0aMICEhgSlTplC6dOkABWyM94nIUeBtYJaqzgpzMItQbYJICk4d0pPvAIqqq8KPbhJcFeA1oKHvQd8AvVR1a54CDwFLcIXA8ePOcoIBA5wW3gsvwCOP5HsSypgxY/jnP//JRRddxKxZs7jooosCFLAx3ubJvSjzyTZbNgXb1q3w0EPOOF2tWs5OKNdck69bLliwgPbt25OWlsb06dNp0aJFgII1xrs8leCcUm1Zc7knpe1FaQq2uDiYNQumT4c9e6BBA+jZE/76K8+3vO6661i6dCnly5cnISGBrl27Eh8fT1RUFPHx8bZZszHBl16ubQXwG/B/wE++n11PpbYWnIkc+/fD0087Y3IXXACvvgrt2ztdmHm63X6aNm3KqlWnLwONjY0lKSnJyu6YiOKpFlw6Z0H4TFTn+F63AVqh2sfN5dm24EQkSkRuz3eQxoRCqVJOUluyxElwHTvCzTfDli15vF0p/vjjjzPOW0UCY0Km3snkBqD6MXCd24uzTXCqmoYz/dOYgqNePVi6FIYPhwULnCUFQ4fC0dwXxdiexTZh27Zty2+UxnhNtIgkiUjbcAeSwe+IPI1IPCJxiPQH9rq92M0Y3Gci8riIVBaRMulHnsPNJRG5UETGirvaQcY4ihRxdj9Zvx5uuAH69oXLL8/12rmsKhKce64nl4Iakx8nVDUx7EsETncnzrrrGcAHwHm+c664SXBdgYeAhZzaaHl5rsPMQETG+bb+WpPpfGsR2SAiP4tIXwBV3aSq3fLzPFOIVa4MH3zg1JuLiYHbboPmzWHlSleX+6tIICLs2bOH4cOH+603Z4wJENV9qPZCtbbv6JWbqt45JjhVrernyO/OtBPIVFNI/NQUEpEa+XyOMY6EBFi1CkaNcva3rFMHunaFX3/N9jJ/FQneeustbrvtNvr06UPXrl05moeuT2OMCyLnIvIfROYg8sXJw6UcE5yIxIrI0yKS5HtdXURuzk/Mmk1NIV+L7VRNIZdEJDG9jpFttWT8KlLE2err55/h8cchORmqV3cqFRw6lOVlnTt3ZsuWLaSlpbFlyxa6du3KtGnTePbZZ5kwYQItWrRg9+7dIfyNGFNoJAM/AlWB54EtOHtfuuJ2s+VjQCPf6x3AwFyF6I7fmkIiUlZ8tYMkm9pBqpqUXseoSIBKqpgIVbo0DBsG69ZB69bw7LNwySVOwktLc3WLqKgonn/+eaZOncrKlSupV68eK112expjXCuL6ljguG8/yq5AA7cXu0lw1VR1GHAcQFUP4+wHFmj+7qmquldVH1DVaqo6ONsb2GbLJjeqVXMWiC9YAOedB3ffDQ0bwjffuL5Fx44dWbRoEapKkyZNmD7d5kIZE0DHfb/+ishNiNTGKcfjipsEd0xESuDb8FJEqgHBGHTId00hVZ2lqonR0dEBDcxEuKZNYdkymDABduyAxo3hjjucbcBcuPrqq1m2bBk1a9akY8eOPPfcc6S5bAkaY7I1EJHSQB/gceAt4FG3F7tJcAOAT4DKIpIMzAOezEOgOVkGVBeRqpLHmkLWgjN5FhUF//gH/N//OV2WM2c63Zb/+hekpOR4+QUXXMD8+fPp0qULzz//PLfffju2m44pYLy1Ds6ZeFgd1b9QXYNqc1TroOo6L7jaqktEyuL0ewqwWFV/z3PQzv3eBZoB5YDdwABVHSsiNwL/BaKBcao6KC/3t626TL5t3+4kt0mT4PzzYeBAuO8+yKF3QFUZPnw4TzzxBFdddRUffvhhlmvpjPESj27VNR/V5nm+3GWCuw1ogtNNuUhVZ+T1gcHk+5dH25iYmB5HjhwJdzgmEixdCo89Bl9/DVdcAf37O1uA5ZDo5syZw5133knx4sWZMWMGjRo1yvbzxoSbRxPcIKA0MAU41WpR/c7V5S7qwb0BXMSpgqedgI2q+lAewg0Ja8GZgFKFadPgueecnVEuusjZGeWee6BYsSwvW79+PW3btmX79u3ce++9zJ07l23btlGlShUGDRpkmzUbT/Fogpvv56yi6qqGlZsEtxa4Qn0fFJEoYOiT0uQAAB+8SURBVLWqXp7bWIPNWnAmqNLSnF1R/v1vWLECKlWCJ56A7t0h024n6fbt28e1117LunXrTjtvFQmM13gyweWTm0kmG4CMgwiVgR+CE07+2CxKE1RRUc5WX8uWwdy5cOGF0KsXxMfD4MF+a9CVKVOGAwcOnHHeKhIYE3xZtuBEZBbOmFtpoB6w1Pe6PvCNqrYKVZC5ZV2UJmQWLXKS25w5Trmehx+G3r0hw2bMUVFRfvesFBFbTmA8IxJbcNkluGxr7qjqgqBElA/WRWnCZuVKJ9FNnw7Fi0NiorMdWKVKxMfHs9XPmrrSpUuzd+9erMfBeEGhSnBnfFCkFHByDyzNxY7OoWYtOBM2GzbAkCHO8gIR+Mc/+PDSS7nr2Wc5lGG/y+joaE6cOEHz5s2ZNGkSFSpUCGPQxoCIHAXeBmZ5qmSOSCMgngz5B9W3XV3qYpJJIvAicBhIw1kLpwGoKBA0luBM2G3dCv/5D7z1Fhw/zpb69bl/yxY+27Xr5CzKY8eO8dBDD3HWWWcxadIkbrjhhnBHbQoxT7bgRCYC1YDvgfQdPBTVR1xd7iLB/QQ0zO/i7lCyBGc8Y9cu+O9/4Y03nB1RbrgBOnWCdu2gbFnWrVvH7bffztq1a+nXrx8vvPACtlm4CQePJrj1QI28Fl50M4tyI5B1LREPsa26jOdccIHTZbl1K7zwAvz0E3Tr5pxPSKDGN9+wdM4cunfvzuDBg2nWrBnbt2/P+b7GFA5rgAvyerGbFlxtnJI5S8iwybK6bCKGg7XgjGepOhNSpk1zjo0bnV1RmjdnaZUqdJo8mf3FizNhwgTatvXGloCmcPBoC24+UAtnFv+pTf5Vb3F1uYsEtxRYBKzGGYPz3V//l/toQ8MSnCkQVJ0q4+nJ7qef0KgolsXGMv7AAc5NTOTpkSMpls1uKcYEikcTnP/Z/C5n8btJcN+oaoHaSM8SnClwVGH1apg+nbSpU4nasIE0YOVZZxH3+OOUS0yE8uXDHaWJYJ5McAAi5+OsxQZYiuoe15e6SHCDgK3ALE7vovTcMgFbB2cigiqsXcv6F1+E6dO5LC0NFeG3iy/m9T17+PiPP/hb+fI81KcP7Tp1ctbdpR82QcXkkScTnMjtwH+AL3Fm8F8LPIGqq8rCbhLcZj+nbZmAMSGwZcsW+t5yCxevXs3twBU5XRAdfXrCK14cSpQ481x27/v7vL/PVKwIZ50Vgj8FEwoeTXCrgOtPttpEzgU+R/UqV5fncfalp1mCM5Hk+PHjlCtXjv3793MxUB0o7jsqlinD0OefhyNHnOPw4VM/Z3X4+0z6udyKi4PLL4caNU4/zj47sH8IJug8muBWo3plhtdRwKrTzmV3uYsWXBd/59XlSvJwsARnIk1I9rNUhWPH/CfBzK8PHYLNm2HdOuf48Uc4emqSG5Urn0p2GRNg6dKBidUEnEcT3H+Ampxeru0HVJ9ydbmLBDcyw8viQEvgO1XtkPtoQ8MSnIk0We1nGRsby9KlS7n88jBXr0pNPT3hrV3r/Lp+/ektw4oVT2/p1a8PV17pVGowYeXhrbraA41xxuAWkouC27nuohSR0sBEdbkOIRwswZlIk5ycTGJi4mn7WRYpUoSoqCiOHz9Ohw4dePrpp6lZs2YYo/TjxAnYsuVU4ktPfuvXO61AgHPOgWuvheuuc45atXKsmG4Cz5MtuHzKS4IrCvygqpcFJ6T8swRnIlFycjL9+/c/rSp4QkICr7zyCiNHjiQlJYVbb72VZ599llq1aoU73OylpTktvq+/hgULnGPjRue9UqWgSRMn2TVrBldfbbNDQ8BTCU5kEapNEEnBKdN28h2cvShLubqNiy7K9Lpw4GztVQOYqqp9cx91cNkyAVNY7du3jxEjRjBixAj++usvbrnlFp555hnq1q0b7tDc27EDFi48lfA2bHDOn3UWNG58qoVXty7Y4veA81SCCxA3CS7jSvJUYKuq7ghqVPlkLThTWP3555+MHDmSV155hT/++IMbb7yRZ599lvr164c7tNzbtev0hLd2rXO+RAlo1OhUwqtfH2JiwhtrBPBkghMp4+dsCqrHXV1uywSMiTz79+/ntdde4+WXX2bfvn3ccMMNDBgwgM2bN5/Rzdm5c+dwh+vOb7/BV1+dSng//ODM/IyNhaZNoVUr57BJK3ni0QS3BagM/IHTPfk34FdgD9AD1RXZXu6iBXcbMBQ4z/eA9HpwrvpAw8ESnDGOlJQURo0axUsvvcRvv/1GVFTUacsKYmNjSUpKKjhJLqN9+5yEN28efP65M3EF4NxzoWVLJ9ldfz1UqRLeOAsIjya40cAMVOf6Xt8AtAamAiNQzbZrwk2C+xloq6rrAxJwCFiCM+Z0Bw8epFKlSvz5559nvBcXF8eWLVtCH1Sg/fKLk+jSj127nPPVq59q3TVv7szaNGfwaIJbjmpdv+dEvkc129lUbhLc16raOP+Rho4lOGPOlNVicYCjR49GVtUCVWdJQnqy+/JLOHDA6bqsU+dU665hQ2frMePVBPcpMA+Y7DvTCbgepxW3DNWrs73cRYIbgVNw7gNO32z5/bxHHVyW4Iw5U1aLxQEqVKjAo48+SmJiIqVKeXb0Ie+OH4elS+Gzz5yEt3ixs0avRAlnSUKDBlCvnjNDs5BWbfBogisHDACa+M4sAl4A/gKqoPpztpe7SHDj/ZxWVe2a+2hzT0RKAm8Ax4AvVTU5p2sswRlzJn+LxWNjY3n44YdZvnw5X3zxBaVKleLBBx/kkUceoUKFCmGMNsj273dmaH7+OXzxhTNDM31sskKFU8ku/ShXLrzxhoAnE1w6kbNQPZDr61Q15AcwDmcWzJpM51sDG4Cfgb6+c/fgjAECTHFz/9jYWDXGnGnSpEkaFxenIqJxcXE6adKkk+8tW7ZMb7/9do2KitJixYppt27ddP369WGMNoQOHFBdtEj1v/9V7dxZ9ZJLVJ2OTueIj1ft0EF1yBDVefNU//gj3BEHHHBQw5APsj2gkcI6hW2+11cpvOH2+rAsExCRpsAB4G1VvcJ3Lhr4P5z+1R3AMuBOoB3wsap+LyLvqOpdOd3fWnDG5N3GjRsZPnw448aN48iRI7Rr144nnniCxo0L1FB8/v31F6xcCcuWwfLlzrFp06n3q1c/1dK7+mq46CKne7OALlHw5F6UIkuADsBMVGv7zq3BlzdyvDwcCQ5AROKBjzIkuIbAc6qa4Hvdz/fRHcAfqvqRiExW1TtyurclOGPy77fffuO1117jtddeY9++fTRq1Ignn3ySlJQUnn766YK5li6/9u6FFSucZJee+HZk2PciJsYpIVS1qv+jTBkQCV/82fBkF6XIElTrI7IyQ4Jb5fl6cH4SXAegtap2972+B6gPPAW8BhwBFmkWY3AikggkAhQrVqzO0YylO4wxeXbw4EHGjRvHyy+/zNatWxGR02ZjFui1dIGwaxd8/73Tutu8+fTjjz9O/2ypUhAfn3UCLBm+/OLRBDcdGI6TAxoAjwB1cdHQgWwSnIg8lt2Fqjo8d5Gecf94Tk9wHYGETAnuGlXtmYt72l6UxgRJamoq5cuX5/fffz/jvYhZSxdof/11ZtLLeBw+fPrnzz8fLrwQqlVzfs348wUXBLX706MJrhwwAmiFs8nIp0AvVPe6uTy7LbrTS/JeAtQDZvpetwUW5inY7O3A2ZIlXSVgZ25u4Os3nlWyZMkegQzMGOOU59m71//fK1u3bmX//v2RucQgP0qXdsr/+KvuoAp79pye8DZtco6FCyE52flMuuLFTyW9zEmwalVnyUMkceZl3INqnrsG3CwT+BRor6opvtdnA9NUtXVeH+q7Tzynt+CK4EwyaQn8gjPJ5C5VXZuLe1oLzpggym4t3d/+9jd69uxJr169KFu2bIgji0BHj8K2bU4ZofTEl/7zxo2QeZ5BhQpOspszB84+2/89s+HRFtyXqDbL8+UuEtyPwFWqetT3OgZYpaqX5vmhIu8CzYBywG5ggKqOFZEbgf8C0cA4VR2Ul/vbJBNjgiOrtXR9+/Zl5cqVzJgxg5IlS/Lggw/Sp08fLrjggjBGG8FUnc2nMye+HTvg00/zNJHFowluEFAamAKc+ktd9TtXl7tIcP2B24EZOHXhbsWpB/fvvEUcPNaCMyb4/BVeTZ9gsmbNGgYPHszkyZMpWrQo3bt358knn6SKbXjseR5NcPP9nFVUW7i63M0sShG5GrjW93Khqq50H2HoWQvOmPD6+eefGTJkCG+//TaqSpcuXejbty/Vq1cPd2gmC55McPnkdkpOLLBfVUcAO0SkahBjMsYUcBdddBFvvfUWP//8Mw888ADvvPMOl156KXfeeSerV68mOTmZ+Ph4oqKiiI+PJzk5xx34jMk1N12UA4C6wCWqerGIVMCZZOK5bQ2si9IYb9q1axfDhw9n1KhRHDhwgOjoaE6cOHHy/UK/ls4DIrEF5ybBfQ/UBr5T30pyEflBVWuGIL48sS5KY7xp7969VKtWjb/++uuM92wtXXhFYoLLbh1cumOqqiKicHJ3f0/K0IILdyjGGD/Kli3L/v37/b63detWNm3axIUXXhjiqIznOPsVZ03V1VpsNy24x4HqOJsgDwa6Au+q6quuAg0Da8EZ413ZraUDSEhI4IEHHuDmm2+mSBE3/wY3geCpFpyIv82eFbgKqIRqtKvbuJxFeT1wA85WKXNV9bNchBpyluCM8a6s1tINGTKEffv2MWbMGH755RcqVqxIjx496N69OxUrVgxjxIWDpxJcZiJNgP7AOcAg3FY7yKmeDjDUzTkvHVYPzhhvy64u3fHjx3XGjBmakJCggEZHR+utt96qc+fO1RMnToQx6siGN+vBtVT4UmG+wvW5vd5NF+V3qnp1pnOenGRisyiNiSwbN24kKSmJcePG8fvvv1OtWjXuv/9+7r33Xj799NMsF5yb3PNUPTiRm3BabH8BA1H9Ok+3ySrBiciDwD+BajgVttOdDXyj+dgAM9isi9KYyHL06FHee+89Ro8ezVdffUV0tDMEY0sNAsdTXZQiaTgb8K/CGXs7neotrm6TTYIrjdPfORjom+GtFFXdl8twQ8oSnDGRa+3atTRs2JCUlJQz3rOlBnnnsQR3Xbbvqy5wdRsXXZQNgLV6ejWBGqq6xF2koWcJzpjIFhUVhb+/u0SEtLS0MERU8HkqwaUTKQ5chNOK24hqrsae3GzVNQo4kOH1Qd85zxGRtiKSlLHbwhgTebLavDk2NtbvInJTwIgUQWQYTjfl/4BJwHZEhiFS1O1t3CQ40Qz/VFLVNNwtEA85VZ2lqonp/fPGmMg0aNAgYmNjTztXtGhRDh48SM2aNZk/398m9KYA+Q9QBqiKah2cXbSqAX8DXnJ7EzcJbpOIPCIiRX1HL2BTnkI2xpgA6Ny5M0lJScTFxSEixMXFMX78eL799ltiYmJo0aIFjz76KIcPHw53qCZvbgZ64BsaA0B1P/AgcKPbm7gZgzsPeBVogdMPOg/orap7ch9zaNgYnDGF18GDB3nqqad4/fXXufTSS5k4cSJ169YNd1ie56kxOJH/Q/XiXL+XSY4tOFXdo6p3qOp5qnq+qt7l5eRmjCncSpYsyWuvvcann35KSkoKDRo04Pnnn+f48ePhDs24tw6RLmecFbkb+NHtTdy04C7GmVRyvqpeISI1gVtUdWDu4g0da8EZYwD++OMPHnnkESZNmkTdunWZOHEil156abjD8iSPteAqAu8Dh4EVOL2H9YASwK2o/uLmNm7G4MYA/YDjAKr6A3BHHkI2xpiQOuecc5g4cSLTpk1j8+bN1K5dmxEjRthSAq9T/QXV+sALwBZgG/ACqte4TW7gLsHFqurSTOdSXQcaQrZMwBjjT4cOHVizZg2tWrWid+/etGrVim3btlllca9T/QLVkai+iuq83F7upovyY+BhnCreV4tIB6CbqrbJW8TBZ12Uxhh/VJVx48bRu3dvUlNTSUtL49ixYyffL8zbfXmqizJA3CS4C4EkoBHwB7AZ6KyqWRd0CjNLcMaY7GzatIkaNWpw9OjRM94rrNt9FcoEd/KDTiXvKM24LsGjLMEZY3Ji232dLhITXI5jcCJSVkReBb4CvhSRESJSNvihGWNM8GS13VdW503B42aSyWTgN6A90MH385RgBmWMMcHmb7svgCuuuAKbqBYZ3CS4Mqr6oqpu9h0DcfYDM8aYAivzdl9VqlShZcuWzJ49mzZt2rBvn6erghkX3EwyeQlYDkz1neoAXK6qA4IcW/rzL8Sp7FpaVTu4ucbG4IwxeTV27Fj++c9/UqlSJWbMmEHNmjXDHVJIFMoxOOB+4B3gqO+YDDwmIikisj+7C0VknIjsEZE1mc63FpENIvKziPTN6noAVd2kqt1cxGmMMfnWrVs3FixYwJEjR2jYsCFTp07N+SLjSW72ojxbVaNUtajviPKdO1tVS+Vw+QSgdcYTIhINvA60AWoAd4pIDRG5UkQ+ynScl8fflzHG5FmDBg1Yvnw5tWrVolOnTjz11FM2LlcAuZlF2S3T62gRcdU9qaoLgcwd2dcAP/taZsdwWoTtVHW1qt6c6bBNnY0xYVG+fHnmz5/PAw88wLBhw7jxxhttXK6AcdNF2VJE5ohIeRG5ElgMnJ2PZ1YEtmd4vcN3zi/fMoXRQG0R6ZfN5xJFZLmILE9N9eROYsaYAqZYsWKMGjWKMWPG8OWXX1KvXj1Wr14d7rCCJVpEkkSkbbgDCRQ3XZR34ZQMXw3MxqkF93g+nin+HpPN8/eq6gOqWk1VB2fzuSTgeeA7EX+PMMaYvOnevTsLFizg8OHDNGjQIFLH5U6oaqKqzgp3IIHipouyOtALeA9nV+d7ROTMxSPu7QAqZ3hdCdiZj/udpKqzVDUxOjo6ELczxpiTGjRowIoVK06Oy/Xt29fG5TzOTRflLOAZVb0fuA74CViWj2cuA6qLSFURKYZTemdmPu53klUTMMYEU8ZxuaFDh1K7dm0qV65s1Qg8ys06uFKquj/Tueqq+lOONxd5F2gGlAN2AwNUdayI3Aj8F4gGxqnqoDzG75etgzPGBFv37t0ZO3bsaecKcjWCSFwHl2WCE5EnVXWY7+eOqjotw3v/VtV/hShG13yDo21jYmJ6HDlyJNzhGGMiWHx8PFu3nllUpaBWIyhsCe47Vb0688/+XnuNteCMMcGWXTWCEydOUNAmu0VigstuDE6y+Nnfa2OMKVSyqjqgqrRu3Zr169eHOCKTWXYJTrP42d9rT7BJJsaYUPFXjSA2Npa7776bJUuWcOWVV9K7d2/+/PPPMEVosktwV4nIfhFJAWr6fk5/fWWI4ssVWyZgjAmVzNUI4uLiSEpKYuLEifz00090796dV199lerVq5OUlGRLCsLAdUXvgsTG4IwxXvD999/zyCOP8NVXX1GrVi1GjBhB06ZNwx2WX4VtDK7AsS5KY4yX1KpViwULFjBlyhT27t3LddddR6dOndi2bVu4QysUrAVnjDEhcOjQIYYNG8bQoUMREZ566imeeOIJv1XFw8FacMYYY/IkNjaW5557jg0bNtC2bVuee+45LrvsMnr27ElcXJzthhIEEdWCs4XexpiCYuHChdx9991s3779tPPh2g0lEltwEZXg0lkXpTGmIIiLi/M7HheO3VAswRUQluCMMQVBdruhpKWlhTSWSExwNgZnjDFhktVuKEWKFPG7z6XJHUtwxhgTJv52Q4mJiaFIkSLUqVOHL774IkyRRYaISnC2Ds4YU5D42w1l7NixrFq1ivPPP5/rr7+e4cOH++3GNDmzMThjjPGglJQU7r33Xt5//33uuusuxowZE9Q1czYGZ4wxJiTOPvtspk+fzqBBg3j33Xdp3LgxmzdvDndYBYolOGOM8SgR4V//+hezZ89my5Yt1K1bl88//zzcYRUYluCMMcbj2rRpw7JlyyhfvjwJCQm89NJLNi7ngiU4Y4wpAC666CIWL17MbbfdxhNPPMFdd92FzTXIXkQlOJtFaYyJZGeddRZTp05l8ODBTJkyhUaNGrFp06Zwh+VZEZXgrOCpMSbSiQh9+/Zlzpw5bNu2jbp169K3b1/i4+Ntw+ZMbJmAMcYUUBs3bqRZs2bs2LHjtPN52bA5EpcJWIIzxpgCrEqVKmdUJIDcb9hsCa6AsARnjCksArVhcyQmuIgagzPGmMImqw2bszpfmFiCM8aYAszfhs2xsbEMGjQoTBF5hyU4Y4wpwPxt2ByOiuBe5PkxOBH5O3ATcB7wuqp+mtM1NgZnjDG5Y2NwuSQi40Rkj4isyXS+tYhsEJGfRaRvdvdQ1Q9UtQdwL9ApiOEaY4yJIEFtwYlIU+AA8LaqXuE7Fw38H3A9sANYBtwJRAODM92iq6ru8V33MpCsqt/l9FxrwRljTO5EYguuSDBvrqoLRSQ+0+lrgJ9VdROAiEwG2qnqYODmzPcQEQGGAB9nl9xEJBFIBChWrFhA4jfGGBNgmYadcDHslFfhmGRSEci4KnGH71xWegKtgA4i8kBWH1LVJFWtq6p1ixQJat42xpjCSWQcInvINOyESGtENiDyMzkMO6H6ASEadgpHJhA/57LsJ1XVV4FXXd1YpC3QNiYmJo+hGWOMycYE4DXg7ZNnnGGn18k47CQykyyGnfANOwFP+64LmnAkuB1A5QyvKwE7A3FjVZ0FzBKR7iJyOB+3igbyW5IgL/dwc01On8nq/dycz3yuCJCaQ1zBEI7vIZjfQVbvuT0Xju/Bq99BTp/Ly3te/X8hEN+Bm/uUEJHlGV4nqWrSaZ9QXUgWw074hp3wDTuRxbATGYadcDGnIl9UNagHEA+syfC6CLAJqAoUA1YBlwc7jlzGnBSOe7i5JqfPZPV+bs5nPgcsLyzfQzC/A7d/3tmcC/n34NXvIC9/zjm959X/FwLxHQTyPgrxmuHvdIUOCm9leH2PwmvZXP+IwgqF0QoPBPPPLqgtOBF5F2gGlBORHcAAVR0rIg8Dc3H+RTFOVdcGM448mBWme7i5JqfPZPV+bs4H4vcfCOH4HoL5HWT1nn0Heft8bv+cc3rPq99DoGII1u8lV8NO5GLYKb88v9DbhJ+ILFfVuuGOo7Cz7yH87DsAXxflR/iWfiHSEHgO1QTf634AOF2UYWVbdRk3knL+iAkB+x7Cz76DMy0DqiNSFZFiwB3AzDDHBFgLzhhjjFsZhp2A3cAAVMciciPwX3zDTqh6YqdnS3DGGGMiknVRGmOMiUiW4IwxxkQkS3AmX0TkQhEZKyLTwx1LYSIiJUXkfyIyRkSs8FeY2H//3mYJrhALUDmjTaraLbiRFg65/D5uA6ars6ffLSEPNoLl5nuw//69zRJc4TYBaJ3xhJzaV64NUAO4U0RqiMiVIvJRpuO80Icc0Sbg8vvA2eIufdPyQGzjZE6ZgPvvwXiYbbtfiGkAyhmZwMnN94Gzp2sl4HvsH6oBlcvvYV1oozO5Yf9jmMxyVc5IRMqKyGigtqTvYGACKavv432gvYiMwhvbSUU6v9+D/ffvbdaCM5nltpzRXiDLOn0m3/x+H6p6ELgv1MEUYll9D/bfv4dZC85kFrRyRiZP7PvwBvseCiBLcCazZUB1EakqHttXrpCy78Mb7HsogCzBFWK+ckbfApeIyA4R6aaqqUB6OaP1wFQPljOKSPZ9eIN9D5HD9qI0xhgTkawFZ4wxJiJZgjPGGBORLMEZY4yJSJbgjDHGRCRLcMYYYyKSJThjjDERyRKcKXRE5ISIfJ/hyLYkUCiJyHRfjbElvti2ichvGWKNz+K6gSLyYqZzdUXkB9/P80SkdPB/B8Z4h62DM4WOiBxQ1bMCfM8ivsXA+bnH5cBAVb01w7l7gbqq+rCLa2eo6sUZzr0E7FXVwSLSDSinqkPzE6MxBYm14IzxEZEtIvK8iHwnIqtF5FLf+ZK+IpjLRGSliLTznb9XRKaJyCzgUxGJEpE3RGStr17eHBHpICItRWRGhudcLyLv+wmhM/ChizjbiMi3vjiniEhJ364aR0Skju8zAnQEJvsu+xC4Kz9/PsYUNJbgTGFUIlMXZacM7/2uqlcDo4DHfef6A1+oaj2gOfAfESnpe68h8A9VbYFTZTseuBLo7nsP4AvgMhE51/f6PmC8n7gaAyuyC9xXZLYv0NIX5w9AL9/b7+LskZh+r52quhlAVX8HzhaRv2V3f2MiiZXLMYXRYVWtlcV76S2rFTgJC+AG4BYRSU94xYEqvp8/U9V9vp+bANNUNQ3YJSLzwampIiITgbtFZDxO4uvi59nlgd9yiL0RTkXpb5xGGsWARb733gUWiMiTOInu3UzX/uZ7xp85PMOYiGAJzpjTHfX9eoJT/38I0F5VN2T8oIjUBw5mPJXNfcfjFCY9gpME/Y3XHcZJntkR4BNVvSfzG6q6RUR2AtcCtwJ1Mn2kuO8ZxhQK1kVpTM7mAj1941qISO0sPrcIp8p2lIicDzRLf0NVd+LUD3samJDF9euBi3KI5RvgOhG50BdLSRGpnuH9d4FXgfWquiv9pIhEAeU4vSq1MRHNEpwpjDKPwQ3J4fMvAkWBH0Rkje+1P+/hFMZcA7wJLAH+yvB+MrBdVddlcf1sMiRFf1R1N9ANmCIiq3AS3sUZPjIVuIJTk0vSXQMsUtUT2d3fmEhiywSMCSAROUtVD4hIWWAp0Di9JSUirwErVXVsFteWAOb7rgloIhKR13FqmC0I5H2N8TIbgzMmsD7yzVQsBryYIbmtwBmv65PVhap6WEQGABWBbQGOa6UlN1PYWAvOGGNMRLIxOGOMMRHJEpwxxpiIZAnOGGNMRLIEZ4wxJiJZgjPGGBORLMEZY4yJSP8PyoFfZMI7O8kAAAAASUVORK5CYII=\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(f\"Energy ({t['energy'].unit})\")\n", "ax1.set_ylabel(\"Expected number of bkg counts\")\n", "\n", "ax2 = ax1.twinx()\n", "ax2.set_ylabel(f\"ON region radius ({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": null, "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.7.0" }, "nbsphinx": { "orphan": true } }, "nbformat": 4, "nbformat_minor": 2 }