diff options
Diffstat (limited to 'plot_llh/syst_variations.ipynb')
| -rw-r--r-- | plot_llh/syst_variations.ipynb | 515 |
1 files changed, 515 insertions, 0 deletions
diff --git a/plot_llh/syst_variations.ipynb b/plot_llh/syst_variations.ipynb new file mode 100644 index 0000000..a7161f5 --- /dev/null +++ b/plot_llh/syst_variations.ipynb @@ -0,0 +1,515 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/smandalia/.local/lib/python2.7/site-packages/matplotlib2tikz/__init__.py:31: UserWarning: matplotlib2tikz has been renamed to tikzplotlib (which is Python-3-only). matplotlib2tikz will no longer be supported.\n", + " warnings.warn(\"matplotlib2tikz has been renamed to tikzplotlib (which is Python-3-only). matplotlib2tikz will no longer be supported.\")\n" + ] + } + ], + "source": [ + "\"\"\"Systematic variations of HESE 7yr sample\"\"\"\n", + "from copy import deepcopy\n", + "import numpy as np\n", + "import matplotlib as mpl\n", + "from matplotlib import pyplot as plt\n", + "plt.style.use('./paper.mplstyle')\n", + "\n", + "import matplotlib2tikz" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/cvmfs/icecube.opensciencegrid.org/py2-v3/RHEL_7_x86_64/lib/python2.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: to-Python converter for std::vector<golemfit::UnfoldingSegment, std::allocator<golemfit::UnfoldingSegment> > already registered; second conversion method ignored.\n", + " \"\"\"Entry point for launching an IPython kernel.\n", + "/cvmfs/icecube.opensciencegrid.org/py2-v3/RHEL_7_x86_64/lib/python2.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: to-Python converter for std::vector<golemfit::AstrophysicalNeutrinoModel, std::allocator<golemfit::AstrophysicalNeutrinoModel> > already registered; second conversion method ignored.\n", + " \"\"\"Entry point for launching an IPython kernel.\n", + "/cvmfs/icecube.opensciencegrid.org/py2-v3/RHEL_7_x86_64/lib/python2.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: to-Python converter for nusquids::marray<double, 3u, std::allocator<double> > already registered; second conversion method ignored.\n", + " \"\"\"Entry point for launching an IPython kernel.\n", + "/cvmfs/icecube.opensciencegrid.org/py2-v3/RHEL_7_x86_64/lib/python2.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: to-Python converter for nusquids::marray<double, 5u, std::allocator<double> > already registered; second conversion method ignored.\n", + " \"\"\"Entry point for launching an IPython kernel.\n" + ] + } + ], + "source": [ + "import GolemFitPy as gf\n", + "\n", + "# setup fitter\n", + "datapaths = gf.DataPaths()\n", + "npp = gf.NewPhysicsParams()\n", + "\n", + "sparams = gf.SteeringParams(gf.sampleTag.MagicTau)\n", + "sparams.quiet = False\n", + "sparams.fastmode = False\n", + "sparams.minFitEnergy = 6E4 # GeV\n", + "sparams.maxFitEnergy = 1E7 # GeV\n", + "sparams.load_data_from_text_file = False\n", + "sparams.use_legacy_selfveto_calculation = False\n", + "sparams.do_HESE_reshuffle=False\n", + "\n", + "FITTER = gf.GolemFit(datapaths, sparams, npp)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1.0 2.5 2.0 0.444444447756 0.0 8.0 7.0 1.0 1.0 0.0 10.0 9.99999968266e-21 0.990000009537 0.0 1.0 1.0 1.0 1.0 56400.0 0.0\n", + "convNorm = 1.0\n", + "promptNorm = 1.0\n", + "muonNorm = 1.0\n", + "astroNorm = 8.0\n", + "astroDeltaGamma = 2.5\n" + ] + } + ], + "source": [ + "fitparams = gf.FitParameters(gf.sampleTag.MagicTau)\n", + "print fitparams.anisotropyScale,\n", + "print fitparams.astroDeltaGamma,\n", + "print fitparams.astroDeltaGammaSec,\n", + "print fitparams.astroFlavorAngle1,\n", + "print fitparams.astroFlavorAngle2,\n", + "print fitparams.astroNorm,\n", + "print fitparams.astroNormSec,\n", + "print fitparams.astroParticleBalance,\n", + "print fitparams.convNorm,\n", + "print fitparams.CRDeltaGamma,\n", + "print fitparams.cutoffEnergy,\n", + "print fitparams.darkNorm,\n", + "print fitparams.domEfficiency,\n", + "print fitparams.holeiceForward,\n", + "print fitparams.muonNorm,\n", + "print fitparams.NeutrinoAntineutrinoRatio,\n", + "print fitparams.piKRatio,\n", + "print fitparams.promptNorm,\n", + "print fitparams.timeSplit,\n", + "print fitparams.zenithCorrection\n", + "\n", + "nuisance = ('convNorm', 'promptNorm', 'muonNorm', 'astroNorm', 'astroDeltaGamma')\n", + "for n in nuisance:\n", + " print n, '=', fitparams.__getattribute__(n)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0.111 0.111 0.111 0.111 0.111 0.111\n", + " 0.111 0.111 0.111 0.111 0.111 0.111\n", + " 0.111 0.111 0.111 0.111 0.111 0.111\n", + " 0.111 0.111 0.00184875]\n", + "[0.11109244 0.11109244 0.11109244 0.11109244 0.11109244 0.11109244\n", + " 0.11109244 0.11109244 0.11109244 0.11109244 0.11109244 0.11109244\n", + " 0.11109244 0.11109244 0.11109244 0.11109244 0.11109244 0.11109244\n", + " 0.11109244 0.11109244]\n", + "2635.5112\n" + ] + } + ], + "source": [ + "binning = FITTER.GetEnergyBinsMC()\n", + "print np.diff(np.log10(binning))\n", + "print np.diff(np.log10(np.logspace(np.log10(6E4), np.log10(1E7), 21)))\n", + "\n", + "livetime = FITTER.GetLivetime()/(3600.*24)\n", + "print livetime" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "def parse_exp(x):\n", + " h = np.sum(x, axis=(0,1,2,3))\n", + " return np.concatenate([[h[0]], h])" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "null_exp = parse_exp(FITTER.GetExpectation(fitparams))\n", + "\n", + "d = fitparams.convNorm\n", + "fitparams.convNorm = 1.4\n", + "norm_conv_p = parse_exp(FITTER.GetExpectation(fitparams))\n", + "fitparams.convNorm = 0.6\n", + "norm_conv_n = parse_exp(FITTER.GetExpectation(fitparams))\n", + "fitparams.convNorm = d\n", + "\n", + "d = fitparams.promptNorm\n", + "fitparams.promptNorm = 2.4\n", + "prompt_p = parse_exp(FITTER.GetExpectation(fitparams))\n", + "fitparams.promptNorm = 0\n", + "prompt_n = parse_exp(FITTER.GetExpectation(fitparams))\n", + "fitparams.promptNorm = d\n", + "\n", + "d = fitparams.muonNorm\n", + "fitparams.muonNorm = 1.5\n", + "muon_p = parse_exp(FITTER.GetExpectation(fitparams))\n", + "fitparams.muonNorm = 0.5\n", + "muon_n = parse_exp(FITTER.GetExpectation(fitparams))\n", + "fitparams.muonNorm = d\n", + "\n", + "d = fitparams.astroNorm\n", + "fitparams.astroNorm = 12\n", + "astro_p = parse_exp(FITTER.GetExpectation(fitparams))\n", + "fitparams.astroNorm = 4\n", + "astro_n = parse_exp(FITTER.GetExpectation(fitparams))\n", + "fitparams.astroNorm = d\n", + "\n", + "d = fitparams.astroDeltaGamma\n", + "fitparams.astroDeltaGamma = 3\n", + "dastro_index_p = parse_exp(FITTER.GetExpectation(fitparams))\n", + "fitparams.astroDeltaGamma = 2\n", + "dastro_index_n = parse_exp(FITTER.GetExpectation(fitparams))\n", + "fitparams.astroDeltaGamma = d" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAFDCAYAAABhmOTuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3U+MI9l9H/Dvz/IfwIddTm92oeEC+sNeGIjhALtkjy7xSUMKyey0AFvs2QQ5eptcG8kpcveOBKR7A0ijbumQUyxWC8ghOWSHlAJPawdwyNHNFw3JdWDdrOZGBoYjaD09XF0EKLB/ObxXxeJ/VrPIqiK/H4DbZFWx+COH279+r977PVFVEBER0dX9RtQBEBERJR2TKRER0YKYTImIiBbEZEpERLQgJlMiIqIF/WbUAazS7/7u7+qvfvUrAMD169eRTqcjjoiIiKLUarX+QVVfXvQ8G5VMf//3fx/NZjPqMIiIKCZE5GdhnIfdvERERAtiMiUiIloQkykREdGCmEyJiIgWxGRKRES0ICZTIiKiBTGZEhERLYjJlIiIaEFMpkRERAtiMiUiIloQkykREdGCNiqZdrtdiAhEBMfHx1GHQ0QUO7VaDbVaDY1GY2Rfu92G4zgRRBV/G5VM0+k0VBWqymRKRGuvXC6PbDs9PUWtVoPjOCOJsVwuI5vNolgsolKpoNFooNFoeMn15s2byGQyM1+30WigUChAREaScq1Ww7Vr15DL5VCr1RZ7gzEiqhp1DCuzs7OjXDWGiDbB4eEhGo0GWq3WwLZCoYB8Pu89vnHjBorFIgDg2rVreP78OQCTEHd2dpBKpbxje70eKpXKXK/vOA7q9Tp6vR7q9frAvtPTUxwcHCz8HsMgIi1V3Vn0PBvVMiUiSrJ2uz1Xa67T6Yzd7jiOl0gBoFAoeMmx1+sNtDq3trZweXnpna9Wq82dSF3lchmNRgO9Xm9gu5ug1wmTKRHRmnG7Wf3a7fbIcVtbW143bCqVGkl6bnLd29tDtVoNHEc+n0cmkxnoTu71etja2gp8rrhjMiUiWiONRgN37twZ2X55eTmSxNwWoptET05OvFatu+/09BT5fB7ZbHbuGNrtNnZ2TM9puVzGvXv3vH3NZnOgdbwumEyJaO2JTL75x+A4zvRj/XK5yceVSv3jfJcsV6LX643tRu31el63rctNru72YrHobctkMuh0Onj//fdxcnICoD/Sd5Zms+kl31KphF6v57WAO53OWnbz/mbUARAR0XjuoB/X5eXlyICeVCqFu3fvIpVKoVareYOJho1LYG4S9bdY/ceVy2Wve7dcLuPk5ASpVAqO46Dk/4thilQqhWKxiJOTk7lbpEHOHxdMpkS09uadtFAqDbYqp5m3xZnLzXfcOG6L0NVut9HpdMYmzFktvq2trZFrou7jcc9zHAeFQsG7btpoNLzjgl7zvHv3LnK5HDqdzszn9no9VKtVJlMiIlo9N9G6A40eP36MXq+H09NTFItFZLPZkaR5eXk5trXoJjS3BdzpdAZG+rqDlcYl4U6n410vdWWzWWQyGezt7eHRo0fedsdxkMlkUK1WvZHCzWbTGz2czWbR6XRQLpe9/dVqFeVyGZeXlwPPixqvmRIRrYFisYiDgwPvVigUkEqlcHBw4CXCUqk0cM2zXq+PLeywv78/kKTc66euSYkUgJcEh5XL5YHndTod1Ot15PN5bG9ve3Hl83mvaziTyXiPd3Z2kM/nvRZ4Pp/3ur3jgMmUiGjNOI6DarWKTqeD09PTkdG6tVoNp6en2N7eHukyrtVqA927LjeRARgZyASYBFsul3F4eDg2QZdKpYHXymQyKJfLqNVqePbs2dhz+vmTt7+reNbzVoXdvEREa6ZUKk285jhP5aFxz61UKqjVami322P3p1IpVCqVid2uqVRq4BpwrVZDp9PBwcHBQHGITCbjJctGo5GYaTRsmRIRJYRbN3eZpp2/WCyG9vqZTAYXFxdoNBrY2tpCq9XyWtCFQgG1Wg1bW1vetWA3kXc6Hdy/fx+dTsfbHgeszUtERBuLtXmJiIhigsmUiIhoQUymREREC2IyJSIiWtBGJdNutwsRgYjg+Pg46nCIiGhNbNQ803Q6jW63G3UYRES0ZjaqZUpERLQMG9UyJSKiydwCCG7tW792u41ms5m41VxWhS1TIqI1NK4+7unpKWq1GhzHgeNfFd0e71ZYqlQqaDQaaDQaqNVqaDQauHnz5ki93kkajQYKhQJExFsU3FWr1XDt2jXkcrnYVC8KAysgERGtmcPDQzQaDbR8i64eHh6iUCh4Lc7Dw0PcuHHDKw947do1PH/+HIBJhjs7O15xeXeR8iDLnTmOg3q9PrKYOWCS+jw1gleBFZCIiGiEf6k0P8dxBrpuC4WClxx7vd5Aq3Nra2ug+HytVrvSuqHlchmNRmNkmbRpi5gnFZMpEdEacbtY/dwFw/22tra8Llh3sW8/N7nu7e2hWq1eKZZ8Po9MJjPQpdzr9QaWUFsXTKZERGui0Wjgzp07I9svLy9HEpjbOhxe69S/7/T0FPl8fuxi39O0223s7Jie03K5jHv37nn7ms1mYpZVC4LJlIhoTfR6vbFdqL1eb2QRbTe5utuLxaK3LZPJoNPp4P333/fWIK3VanMPGGo2m14CLpVK6PV6Xiu40+mwm5eIKJFEzM1vd9dsOz/vb3Mcs80//aPbNdvS6cHn53Jmu2+QD46PzTZ/hTX//oB6vd7Um1+tVpu41ui45OUmUX+L1X9cuVz2unfL5TLy+TyKxeLIKOBZUqkUisXiwMLgswR9jTjYqHmmrdbo/0+uSqX//4/jAGNGlXs2aAA0EUWkVquNjIIdlkqlvO7Zaa29ra2tkeTrPh73PMdxUCgUvOumjUbDO+4q1zvv3r2LXC6HTqcz8/m9Xg/VajV581lVdWNuQE5NKhy9VSrqqVTGH+PeXNns5GP29/vHNZvTz9dsKhHRlVWrVT05OfFuxWJRM5mMnpyc6MXFhaqqplKpgefU63XN5/Mj53r+/PnA9ouLi4HH9Xpdnz9/PjGWi4sLbbVaI9szmYxms9mB51YqFa3X61oqlQbOn8lktFqt6sXFhfe4Xq97x7ZarZHnXRWApoaRX8I4SVJuuVzuap/2BEymRBRHlUpFs9nswLaDgwOtVqsTH7uKxaKXgF2ZTMa7P+45ficnJxO3+89zcXGhxWLR2+c/73Ds/iScz+e9+IrF4tTEPo+wkimvmV7Bm2++CRFBuy0Axt/OzsRboWZnZ/JxgNn/5ptvRvJeiGi9OI6DarWKTqeD09PTkdG6tVoNp6en2N7eHrnGWqvVBrp3Xfl83hvpOzyQydXr9VAul3F4eDi2+lKpVBp4vUwmg3K5jFqthmfPnk08r8vfHe3vKp71vJUJIyMn5QZA43mrTGy1Dv2BNrWF6++qJiIKalqrs1qtzmyVBn0ttxVbqVS0Uql4LU63W7ler6vqYEs1n897rdFxreigwJZptG7duhVikn816rdDRDRxNLC7b9r+oDKZDC4uLtBoNLC1tYVWq+W1oguFAmq1Gra2ttBut70WtXv//v376HQ63vY4YG3eReRyQLsNNJvmPmCGxL/3HnB01B8e32oBOztANjs4TD6dBp4+RRrA09CC2scf/MG/wN/+7X8I7YxEROuKtXmvYtzcmEXnmoXg5he/GOLZHPzkJ0ykRESrtFnJNGytlrlc6bZKAdMaVR2ctJ3LmW3Dk7e7XUAV//3RI+j+vrmAWqn0u4AfPDDbbt8e7BqGvdjq33b7NhTAbZyDiIhWa6OKNiCXM12yfudjkk+pNNgqBUyLdIO6xImIaH68Zrpm3F7sDfpnJSK6srCumW5Wy3QD7MOtaZmwUlxERAm2VslUREowc4ZGF+/bEA7cydJMpkREq7I2A5BsIi0DWL9VZwNwsA8H/bUH3YUtxt38l4Xdgc7jbizOREQ03VpdMxWRCoCqqjbG7d+Ma6YfALjt29IEkJtwtAN4LdksgMlLRa3R14SIyLO280xFpCgiYxe+E5EDu79kW6I05Nat/4rrMFNnngAAduDWAG5CoBBkvbrAZRzZY4/Q9o7L2uOaEJhEPCkZExEREKNkKiJ5ETmAaSqNLLBnE2xbVWuq6gDYFpHwalutiQ8++ADdJyaNpq9fH5iLmsua7t9Ws+ltOz46AgAcHx1521q29W6Ob9sbERFNErtuXps0U6paHtr+XFWv+R7nARyqasG3beO7ecMmdq5N3L4nRERh2KipMSKSHbP5EkB+1bFsmid4MeoQiIhiLxHJFGaE7vCidT0AEJGUqvZsl++O3Xa5ydNjwpTGJ1GHQEQUe7G5ZjpDCqNTXtzkugUA9lpqTlXLkxJpt9v1Fuz23479dXRpQBpPkLZDmYiIaLyktEx7Y7a5yXXuZdbT6TS63W44EW0Is0AcERFNM3cyFZHX7d2O/VkCAFX9TthBjXGJ0RG+Kfv64xItERHRygTp5n0HQEZVfwngEYCXADwSka8uJTIf2207nDS3AIwdtUvhqaCECksTEhFNFSSZ1lX1ByLyeQA5Vb2rqh8C+GhJsQ1zhuaVFgBUVvTaG6uEM5RwNrBtUtlBEbOuustdY52lCYlo3QVJps/tzzyAmm97KBMQRSRrizYUAdyx1Y68KTGqegggYysgHQC4UNXapPON4x+AxEFH8ynh36OErw0M2JqmXC55x5XLoy3ahw+XFSkRUXSCDEDKicg1AIew10tF5CZCKixvu3LbAE6nHDNx3zw4ACm4J7c6ePjwIW4DOLe3L6OfUN2/pPwp9gGAXQC7OMMPbat2H6YS8Ae4BeCD5QdORLRCc7dMVfXbMImzrKo/EpGvABhXTIHWyAcffABVxfmDBwCA3du3B0oUuvzbdm+bQvvnDx5425yK6ZF/gldX/yaIiJZs7nKCInJPVe8uOZ6lYjnBaLE0IRHFTRSrxhyKyNsi8sKiL0pERLROgiTTgqp+D8ANEfmKb95pYnAAUrSu2xsR0bqZewCSqj7y/xSR10Xkr2BWafnekuILFQcgRYufPBGtq7lbpm5L1CbR+wB+BDPHtGVbqn+8pBhpTXRxHV22TYloDQWZGlMTEXfkyCmAfVV1lxT5EABE5O2ktFJp9V61bVMOPyKidRO00P07bjfvMBH5YgjxEBERJU6g0bzDidR277pJNId+EXwiIqKNEaRow/cBQEQ+594AtADs2f3fVtUfLSPIsHA0b7SayKGJXNRhEBGFLsgApDdE5BImgbZgSv9d2PuJkE6nvYo8TKarl0MbOfTXbW+1phfNb/m+WaXS+GNYOJ+I4iBIN28ZZrWYlwCUVHVLVT+F1a0aQwmXw39BzhbNPxZBbkdwhGNvfxYtKMRrve7s5CAi6IrAORNc902uqaAEheDPHjKbElH0gi7B5ibOa77tL4YYD62xT9/632jjm0Nb34Mpky8A3Ipebfu4PXTsq75jz9DFP8NfInG1Q4hoDQWpzfsVAGrXNL0J4EV7/6uq+p2lRhkS1uZdL6z1S0SLiqI2bwfA10TkBTuq9x0ReQZge9EgiIiIkixIOcEP0e+Hg6p+SUTesNuJVu521AEQEVlBWqYjVPXDJK0iw6kx68VdrJyIKGpBKyCNcwbgrRDOs3QsdL9ezm3bdDfiOIiIxiZTEXkRwHPMLqMqcxxDtBRftu1SfgGJKGpju3ltAfuaqn7KvQH4EoDXhrbtALizwniJiIhiZ9o10/2hxy/65pkCAFS1DdOCJSIi2lgTr5n6llebJRVSLESBKMS7R0QUpSCjeb/gLhDuso+/EG5IRMG4I7RFSlNr/faPE4i0Rva/8srjqN8KESVUkFVj3gVwKiLPROSxLdhwBozUhyNaiVde/jHEXo14AEBxhtu+yTL7cKAQVFDytl2Hacc+wb8aOFcTOfy3j4+XHzQRraVA80xV9UsACgAcAHdU9Yaq/nIpkRHN8Itf3ICqA1XF7m0zTeb8AaBqbk7FHFfad7cpuk+eAADS13/LO07VrGjzJh5G9VaIKOHmrs27DtLptD59+hQAcHR0xMIN5MnJ1wEALf1GxJEQ0SqFVZt3o5IpC93TJCyaT7SZoih0T0RERGOEUU6QKPGO8OmoQyCiBGMyJQJwjJ9HHQIRJRiTKRGAYxzZn0REwU28ZioinxuzbV9Efioi/ygifycif7TM4IhW5T0c4z2mUiK6omkt0z0RuVDVHwCAiPw5gG0AhwB69v6fioi4xxAREW2iacn0DMAFADdRdlT12779jwA4InLPdwxRImXRsvdykcZBRMk0LZleA5ASkc+p6v/F5GrindCjIlqxFtxpZpxnSkTBTZtnWrQ/L+3PLXeHiLzgOy4xq8Z0u12v0DmrH5FfC1m0sO19P7qSAkSQlq5XCN+REiCCknzdO27X7jyX3YGi+RBB/Xf+MOq3RUQrMq1l2gHwDoA7YsvD+FqpPRE5AfAMpk5vIqTTaXS73ajDoBj61y9/Fx9//E6o5yz8+q9DPR8RxdeVygmKyItAoDVPY4HlBGlVRMxKNaqJ+VuTaCNFWk7QJlFeXCKa6MzeiGgTTE2mIvK2iLwvIv/Rt21fRP4JwHMR+Z9Lj5CIiCjmphVt+BbMNVMB8G/tguCvAyjb258C+A07NYaIfB7gM3iAz0QdBhGtyNRygv5+ZBHJAjgZ6lt2bNIlIp9d/H3UIRDRCk1Lpj/1P1DVtojUxhz3D+GGRJR8u3gAADiPOA4iWo1ZRRsgIm+r6vfstvvuThF5QVV/CdMNTEQ+P8Ru1CEQ0QpNG4DkduEeukXv3akwdmpMT0T+AkB92UESERHF2cSWqU2c79rbyD4RKQBoJm2uKdEq7Hu1TEqRxkFEq3Gl9UxtF++jsIMhWhcOyvYekynRJpg2NeYNO8f0Bd+2b/nmmP6jf/4pEfU52IeD/ajDIKIVmXbNdAemG/eXgLee6TOYNapeA3ADwMsi8tWlRxkSFrqnVSnjDGWcIW0r33clPVAIvyU5QAQ5X9H8Y7kOiOBYjr3jctJi0XyiBJjWzftMVf310HpDjz8C0LZJNhFY6J5W5datW3j48GFo59v69a9COxcRhW9aMv0CBhf9vpxwHGv0Eg354IMPBh6nMfw/Ssv3328M7Dm2NyMHUy+F/6MRxdm0bt57ItK09Xk/B3Od9IvuNVQR+ZyIvL2KIIk2W9veiCiuJiZTO+XlJsy10w7MfNI67OAjmD+qU6r6nVUESkREFFezavN+AlPs/h0RyQD4PEx3b09VP1pBfEQb7wlejDoEIpph7nmmqtqBaaES0QqlwbooRHE3az3Tr4rIMzun9K9E5LO+fW/Yeadc05RoidJ4gjSeRB0GEU0xsWVqp7y8BVNO8BLAl2CmwnxRVf+Pqn4oIp8A+DsA/2Yl0RJtoKdIRx0CEc0wrZv3xtDapd8XkUMA90XkW6r6I5giDkRERBttWjJ9PLxBVXsAvmS7dwF3mhwRLU3Fq+/rTD2OiKIz9Zqpy12CzaWq78Ksd3on/JCIyK+EM5RwNvtAIorMtCXYvm3XM30G4FsiklPVv/Ht/76I5MHFwYmWqoSvAWC7lCjOprZMbQvUAfCaP5H69jdgWqhEtCRn+CbO8E3s2ur357I7UDTfvSNS8ormn8tnARHsyrl3XEkcQASOlPDKKyNXcYhoATO7eVX1k2kFGrg4ONFy3bp1K+QztvDxx18I+ZxEm01UN6d89s7OjjabzajDIIqUHTyITfp/n2gSEWkNzVy5krkrIBHRergedQBEa4jJlGjDcEVfovDNNTVmXXS7XW+AxvHxcdThEEWii+vosn1KFKq5W6Z27dLmuFG9SZFOp9Ht8u9y2myv2rYpr5gShSdIy/RL4za6i4UTERFtqiDJ9H2MX4KtNGYbERHRxggyAKkAUwmpA6BntwmAmwC+E3ZgRLQcTeTsPZbWJgpLkGS6A+AUZjk2v1R44RDRsuXQjjoEorUTJJkequqj4Y22pUpECZGztX7ZLiUKz9zJVFUf2cFGdwB0VPVHIvK6qn64vPCIKGxtfNPe+0akcRCtk7kHIInITQBtmFG9Wbv5IxH542UERkTLlB0olj96y3lzskWcKce18Oabb0b9ZogiF2Q0b0FVX1PVOwA+BLwi91yCjShB/sdrr0HRxhGOvW1ZtKAQ3+Ak4wkARRnXfXWTKihBIdi3i8I9fPhwFWETxVqQZPrjCds595soQf7d7/0eAOD4CFA1t5Zd/yGXdbe1oKpIXzeVkrpP+seW9s2xTgUw4xKJaO5VY+xC4X+hqj8TkS/aa6afA/COXfc09rhqDFG4du0KNOdcgYYSKopVY+4BaImIAui5yzjBzDMlog10HnUARDERZDTvJwBeE5EigM/DjOj9/tIiI6LYO8dtAMBuxHEQRS3wEmyqWltGIESUPF+2bVN28tKmC7QEm4j8uYhcisg/isgzEfmTZQVGRESUFEGWYNsHcAPAHkzB+20AJRERVf3ekuIjIiKKvSDdvNt2jqnrIwANEfnzkGMiooRQb5o5O3ppswXp5n08YTurZhMR0UYLkkyvTdj+ontHRN5eLBwiShLBPgT7UYdBFLkgRRvuA3gDgwuEZ3yPBUBOVV8KNcIQsWgDUbjc+ebz/h4hipuwijYEaZlmALwLwPHd/I8r4KpORBvnAT4DiGBXzr0C+CVxABE4UvIK5qftzq6kB4rltyQHiCAnX4eIsHA+JdLC65n6cW1Tos1y69Yt4OFPQj0nC+dTEs3dzRt3IlKC6XLOAGio6khiZzcvUbwdiymsf6xPI46ENkUUtXljS0QyMFN3HPu4CjMflogS5Bg/jzoEoitZi2QKoAjgwvc4O+lAIoqvYxzZn0TJErtkagvp31DVwzH7DmC6crcAwG2JAngJg6OMISIpVe0tOVwiCtF7No0eRxoFUXCxSaYikodpURYwlBjt/hMAdVVtuI9FpOgrvL+1smCJiIh8gtTmfd3edRNdCQBU9TthBGKTZENEXgKQGnNIaai1WgdwCKAG4NnQsVtslRIlT9abXZeLNA6ioILMM30HQEZVfwngEUzX6iMR+epSIvMRkXHXQC8B5O39Guz/fSKSAtBYdkxEFL4WdtDCwgMriVYuSDdvXVV/ICKfh6l0dAPwRtIu2xZM8vTr2ddPqWpHRFq+ruKR661EFH8tO3aQ7VJKmiDJ9Ln9mYdpCbpWMVE1hdFrom5y3QLQ8w1GYquUKKHcNul6zH6nTRKkmzcnIl+BafV9FwBE5CZWM/Bn3PVP93WHW6wTdbtdr7SZ/3Z8fBxKkES0qDa4EBUl0dwtU1X9tl0gvKyqP7KJNYN+i3WZLjE6KCll45p7oFE6nUa32w0zLiIiomBTY1T1zPewo6rfDzmeSa/bFpHhpLkFdukSrZUn/RUdiRJl7m7eMWuV9kTkpoj8ccgxTeLYgg6uAsxKNUS0JtL4BD/EHW9FGUdKgAhKdkUZEcGu3XkuuwOrz7h3xLdSDVegoVUJ0jId6GZV1Y8AfBTWguB2+ksepjTglohcwBSsb9vXOxSRA5tQMwAufAUbiGgN1H/7X+L418ehnOsBPhP6ijZEk0xdNcZeI80BuAYz5WR4ZEAGQFNV/3RpEYYonU7r06dmNYqjoyMOPCJaZ3bhcqzJyli0HCtZNcZeIz0Tke8C+BDA/aFDOqr64aJBrAoHIBFtjl08AACcRxwHbYZ5u3kPAeRXNeCIiGhRP8Ru1CHQBplrAJKqfjIpkYrI58IMiIiIKGkCrxojIi8MbToEkIhrpkS0OfbhFkUrRRoHbYYgU2P2ReSfYIo09Hw/E/NN9VdA4uAjovXmoAwH5ajDoA0RpGW6DeCaqn7i3ygi3wo3pOXhACSizeFgH0CC/tqnRAtSm7c+nEite2EFQ0QUljLOUMYZ0raYQ1fSA0UeWpIDRJDzFYQ4luuACI7l2DsuJy1ABC3J4ZVXHkf9tiimgiRTnTDYaD+cUIiIwnPr1q2Qz/gJPv74nZDPSetiatGGgQNF7gN4wz506+QKgDdU9VNLiC10Ozs72mw2ow6DiBJIbBGIeX9nUjKspGjDkAyAdzG4HJoAOFg0CCIioiQLkkwPVfXR8EYReRZiPEvljuYFWE6QiILhijY0zdzXTFX1kYi8ICJvi8gXAUBEXk9aOUFVhaoykRJRIGl8gjTGjcGkad588825B3iJ5Lxju5ICRJCW7sAqQvXf+cOo39JYQeaZ3oQpdP8lmKL3gFk1ZlVLsBERRSaNJ0jjSdRhJM7PH76OLL4WyrlKOEPh138dyrnCFqSbt6CqrwFeYoWqfiJuvykR0Rp7inTUISRSC980d3wDt47tzcgBUPvfUf7KACUxSdkZc1zUgiTTH0/YzqFtREQ0Vst2ZOZCONeZTcwOvhHC2cIVJJl+wQ4h/hlsArXzTr8A4Afhh0ZEFB8Vr5ZSHNtF8bWDFoD1b3UFSab3ALRERAH0fL27N0OPiogoZko4s/eYTKNyO+oApghSAelde830Lswi4d9S1Ruq+svlhBY+Fronoqsq4WsohTSQhq7mHPFd7D1IBaR/gqkZfT9JCdSPFZCI6KpYAelqumIGbqV18UVGzsUs+L6r4aXUKCog7anq90XkpohkAFyo6o8WDYCIiNZXGk9DO9eXbbs0jn/OBCna8H3785GqngHoiMhfichXlxYdEVFM3Ea8r9mF5fErrwAiOJfdgVV23DsiJe9y2bl8FhDBrpx7x5XEAUTg2OPSwEZMKgpStOF196eIfBdAB8AnMIUciIjWWpyv14XpP3/8sTedJQxPAbwR+go+8RPkmulP0W9dVwCcTVjfNLZ4zZSIrsq9XtfvajTXUAX7gB3p+wCfwS7+Hrt4gB/CHL8PBw7KcLCPsj3uOoC/fPll3PjFL1b7JuYQ62vD7iySEGOL4popALwzrtg9EdG6+5OX/xM+/vhGKOc6xh/hxsf/K5RzUTwEaZl+xb1umlRsmRJRHIQ5wjVsJfk6AMDR+FUZEjGFM1TDm+sbVst07mS6DtLptD59akaWcQk2IorKEnorwxPj4JbRBb30bl4ReRtACmYx8Kaq/o3d/iKAPMxi4QUAP1XVP1s0kFVIp9PoduP3lyARUVw42AcAr3gizWemRHWQAAAStklEQVRiy9QWaci6SXTCMSkAz1T1U0uKL1Ts5iWiOIhx4y/WsZ3LZwEAu/qz0M65igFIzrRECgCq2hORs2nHEBHRoKa3hkor0jiSZhd/H3UIE01Lpj9174jIPswKOvsADlX1O+OOIyKi2XIxnp5/3VtBNH6lFnbxAEA85/tOS6besjC24tGZiGAokQ4cR0REs+Vswfw4tku7eNXei18/rzt3N46mJdNxn+S4Vmj8PnEiohhr20WuEcNFrrt4EUAc26XxNi2ZvuVbs9RVGLPtLQDDrVUiIkqgV2EK28WxlbTvrSUbv7HG05JpDsC1Mdu3hx5/PrxwiIjW3xE+HXUIieSgbO8lK5mequq7s04gIt8KMR4iorV3jJ9HHcIU4RW5D1uc58BOWzWmMuc55j0uct1u11s6iNWPiCgqxzjCMY6iDmOspr3FURln3mIBcbNR5QRZtIGI4iDOhRHiHFycywnOvZ4pERGFJ4sWIIKW5CCS83rNupICRJCWrrfgtiMlQAQl+bp33K7d6S7i/corj0OJK4cmcjFtm163tzgKugQbEREt6OWXHwMfh9uWCWt5uLZXnSl+4lxZnd28REQJF2bPbIx7eZeydF1Ui4MTEVHMqFeIbvEMeIRje+94ylHReNW2TWOY55lMiYio7xjvefdofkymREQJJ3b+ZRgttmNbUOI4hHNtEiZTIqLEc+deOlOPmsd7tqDE8cJnCl+cl67j1BgiojXwAJ8BRLAr596UmpI4gAgcKXlTatJ2Z1fS3nEiQEtygAiydkWbOMqhHdvl69gyJSJKuFu3bgEPfxLu+WIozkvXcWoMERElAisgERERrbGNSqYsdE9ElFxH+HRsl6/bqGum6XQa3W6cC1IREdEkcV66bqOSKRERJZe7bN1xtGGMxWRKRESJ8J5No8eRRjHeRl0zJSIiWga2TImIKBGy3gzT+C0Tx2RKRESJ0II7HTR+9RGYTImIKBFayAKIY7uUyZSIiBIivu1SDkAKR60G7O0B166ZitG5nHnc64X/WqenQKFgbnt7QLkMOA7Q6ZjXDRq3G3OtFl6Mjcb4z6NQMPdzufBer90GtrfN65ye9rc7jnn9RiOc1wn7fER0BW17ix8m0zAUi0C1CuzYv5vKZfM4lQrvNRoN88u8Xjfndn9WKkAmYxJKO+CXrFgEzs5mHxdUPj/4edy924+51QLeessk18PDxV8rmzXnHnZxYf6YCesPmrDPR0Rrhck0CXo9k3wAk5CGk3Q+b5JqUpRK5qe/JRm2kxPg+XPzB0Mcz0dEgT3Bi3iCF6MOYywm0yQ4PDQJ9e7dyceUSuG2hJfJH2fQ1vRVXyeO5yOiQNL4BGl8EnUYYzGZrsrhYb9rc2+v39IETKIsl/vXQAuFwWuK7nW6TGb6a3z0kfnpv47odqW63cTD1xZd9bo59vDQvP5wF+ysGIPwd5Vms4PXbhsNc/7tbXOd8iqv3emMXkf1v/9azZz79LR/PRcwj4e3zXM+xxl/PletZrafnvbj9197nfb+Haf/2rlc/7NzzKLPKBSCf/5ECZXGE6TxJOowxlPVjbnlcjldqnxeFVCtVEa3p1L9xwcH5rjnz83jTEY1m+3vb7XM/osL8xgwt3p9/ljqdfOcg4PRbScn/W3VqtlWKg0+P5UycbtmxTiO+3lUq4PbS6XRz6lSMduKxf7ru/dnvbb72P++xm1z338+3//s3dcdt80fd9Dz+Z+bzZptrmJx8N9/1vt3X9t9rGqeO/xvRrTm3F+F4Z4TTQ0hv7Blumzttmlt3LnT33ZyYga0pFKmVdLpmEE5rqyZS+VdB3VbpEEGv2xtBYtzeCRwPm/i7nTmi3GaSsW0ttzWJWDev3vt1B/vjRvm5/PnZmDRoq/t575GJtPvsnU/23HbOp2rn8//3LMz0/J3ua1Jf+t00vsHzPvN5/ufBQDcuxfOAC4iCgXnmS6b+8tv+Hrb8C/dx48Hu1+LReCll8z9fN506z1+PH0AzN7e+JGtV+H+cu905otxmr29wcQ5zXBX9qKvPc64a5/b21c716Tz+bld2e6172lJelJX/uGhSb4nJ+bW6czu9idaMxW4v0ecSOMYh8l02Wa1ctxW1o0bwMHB+GNOToD7901CPTkZf0ynM7s1GqRle3lpfvp/YU+LcVnm+XzirNczrf5UCnj0yPx0nODzVfN581m415GnDUYjWlMluFP54pdM2c27bG4XXaMxmMxOT00XsPtLcrjLstHot8RSqX6Ls1AYTcydjulC9Sfaca0lf1fjsIuL/v1ez7x+sWiS6TwxLkuUrx0Gt6v8rbf6/yb+zzoIN4E2m/0/Mog2SAlfQwlfizqM8cK48JqU2/Xr1xWmEpUeHR1d6WL1WNWqGRySSpmr49mseewfYHJwYLaXSub+8KCcgwMzkKVUMjf/QBe/kxNzXD5vXmPWsZlM/zXdQS6pVP851ao5V71ujnHjGHfOeWOs12d/HsPHu4N0stlgr91q9Qc6ZbPm/Qxvq9cHt2Uy5rMY3latmmPDPJ+qiTebNe/h5MQcl8n0nzfP+3dlMsEGohGtEff3d8jnDGUAkphzbYadnR1tNptRh0EUjP/6aKEwvYeBaI2JCADTCAzxnC1V3Zl95HTs5iWKs8PD/uAodx4w0Ya6bW9xxAFIRHH21lvm2rqbVPP5qCMiisx51AFMwWRKFGfZLLt1iaxz2y7djTiOcZhMiYgoEb5s26ZxHOnDa6ZEREQLYjIlIiJaELt5iYgoERTi3YsbtkyJiChRRAQignP5LCCCXTmHiFmVsCRmeUJHSt5xabuzK2nvOBGgJbnZLzYnJlMiIkqEV17+MQT7oZ0vh3Zo52IFJCIi2kg5+Tra+CYrIBEREV1VG98M7VxMpkRERAtiMiUioo10hE+Hdi4mUyIi2kjH+Hlo52IyJSKijXSMo9DOxWRKREQb6T0ch3YuJlMiIqIFMZkSEdFGyqIV2rmYTImIaCO1sHCtBg+TKRERbaQWsqGdi8mUiIg2Unjt0jVLpiJSEpHw/tQgIqI1Fl6h+7VJpiJSAlAGsBV1LEREtFnWJpmqqgMgsUvCHB8fRx0C0ZXx+0tJ9AQvAkA6jHOtdAk2ESkCuKGqh2P2HQDowLYsbXIMev4KgKqqNsbtj/MSbCKCTVoOj9YLv7+USCIQAKoqi55qJS1TEcnbZFkGkBqz/wRAW1VrNolu28RLRES0FGk8Ce1cK0mmqtpQ1VNMvtpbGmpN1mESLwBvYNHBmFs+SBzdbjdQ3PN2Xc1z3CZ2g0X5npf12mGc96rnCPI8fncXw+9uuOcI+rxVfX+fhtPDC2D13bwnAFKq6k+UWQCPVPXa0LZW0Kb3rG5eEdEg73ferqt5jpt1zDp2k0X5npb12mGc96rnCPK8VX53g8aWBPzuhnuOoM9b1fdXBAAkOd28M2wBuBza1gMAERnpEp7EdgvvANjj9BgiIpqlglJo5/rN0M50dSmMTmdxk+sWbGKdRVVrAGozDvt/IvJbY7Y/BTCuDzgtIvP0Dc9z3Kxj5n2tJInyPS3rtcM471XPEeR5q/zuBo0tCfjdDfccQZ+3yu/vPw8Q10RxSKbjkqWbXIdbrAtR1d8O83zLZrutDwFkAEBVw5thTLRktmfpDswofUy6/EIUN/a7+xz2uwuTk25O+x0ch2R6idERvikAUNW5WqVrbAfARwCccdOJiGKuqqoF+4vpDACTKSXFjv86qoiUZjVmIk+mqtoWkeGkuQX+jwcA92z3NVGi2JH2bou0B2Av2oiI5ufvRRGR4jx1DyJPppZjA3YTRwFAJcqAwrRAsYot+0spC6Cmqp3h5xMt2xW/v1m7PwsgD6DByxQUhRCKBc1VonYlydT3P1QRJkFcwPc/l6oe2nmjRZjrgxfr0CLzJcIC+n3v/v0nAOruX0EicuL/o8L9hxWRDoAqgNyqYida9PsLYMv2PHUAtABsryh0ojC+v27N97l6SVc6z3RTjZtfa7c/H5pfmwdwaK8zFQFkbLGLkWOJVuWK398SgG23NWDneC88l48oqKt8f33bWqo6VyMmDvNMN9KEubCXMC14wIxyrtljM0hwEX9aP3N8f+/DjkK3A5B4iYJiY47vrysz7znjcs10E00tVqGqDREp2n/0G+AADoqXWd/fnoi8b1uoKfD7S/Ey8/trt839RyCTaXRmFqvw9d0n/voxrZ0g31+iuJmrWNC8XbwAu3mjtLJiFURLwO8vJVno318m0+iwWAUlGb+/lGShf3+ZTCNipwWxWAUlEr+/lGTL+P4ymUbLGVoEfa2KVdDa4/eXkizU7y/nmS6Rr1hFGeavnnsYqgTjq8CRgRm0MbNsFdEq8PtLSbbq7y+TKRER0YLYzUtERLQgJlMiIqIFMZkSEREtiMmUiIhoQUymREREC2IyJSIiWhCTKRER0YKYTIkSxK4NupFEJCMiVbvYc5RxlESkMlQ9hzYcl2CjxLCLpJcBHABoA3jft3sbQFVVY1EbVkQqAKCq5ZDOlwJwBqAIQKYcN+szKgEoJ7hS0WNVPfVvsO/5EMCFb7P7PcjMWgrOJsW7ALIAHACH/mLndk3Wij3noao69jWZTMnDZEqJoaodAIf2l9/7Y36pVkUkE5NEUcXQUk4iUrpqbPaX+56ITC1ZNsdnVIEpsbYW7PssA9gbSoB5mAQ4sxWrqjURaQB4DvMHWW9ovyMi26p6GG70tE6YTGltqOqeiKiINGxSiTKWcS3kuRcaXhZVbYvITtRxhMG2Ds8AfH5MAmyISHv8M0epas8m1DKGVg6xvQKPQwiZ1hivmdK6aWCO1sgyiUhKRPK20Lb7uIL+4sNRxFTyPWxGFUfIKjCFyyetPxl0BZATAMUx16XvzOoqJmLLlNZNHaZ1AcBrVdyFaVncAFC3rZY8zC/Pjn3Opd1/4e+Ktc8v2eMAcw3u1LffvTaZQn8Jp56Nwd2Wh0mkWXu8tzrFpPh85z9Bv1W0yKLbXqvYXTVj6DNwE09hwmcw7TNswnyGb9l97jXFsu85/uuZl+j/wbNnW8t5e47TAN2p7oogY9kYM7Pex9DxPZh/b3/3+MYO+qIAVJU33hJ1g/nFfDBhX8l8rQeOzQw9Ttn7RZjE4X9+1X9uAC33ePs4D6Die62ib18RQNZ3XH1oX3XCe5kUX909n32c8b+3OT6jOsxApMqk57mfgRsDTOLQADGW3McwA3jc99+a9Fn4tg1/9qUZ7ynj/tu4cfo//zk/k7Hvw7ftxB+XjTM7LRbeeFNVdvPS2knBtuDcqQs6eP20DeCO7/HwtdUKTOvFbblBfd2IaloyJdvK6QA4sVMl3FGjc1+rnRaf7SLeUd/aixr8OnBdVU/VjCj2WntD3Zg9/7nd9+oeM8dneAmgo6o9VW374s2i/1l0AAxcp7WfI3xd4XkA9+d9Y75/k5Guc/vvcSIidTsoLTvndwEw//4ZNy6YRDr3tVfaXOzmpXWzjf41wQyAnpsUrceYnvA66HfrZTE0ItfqwSS6hoiUYboaK3bAy80AsU6Lb2dGnIHo4KjePAD/NcBp3cfzfIbjnl/zvc6kRHkC84fLHkwLMWg3dgOmW3pghLT2u9AVZhqQ240887ugqh13IJKIHE54b0QjmExp3dxBP6G5rS3/6MxZ81Az6P8C7WD8oKEUgI6IFG1r1G1lnWD0etsIex3PbdmOjc/+0s+MPjsUQc57lc8QMPNbt2yL8FLHz7e9D9Oy97digygDaInIPIk4yPuowHT3txCgtUybjd28tDZEpArgntstZxOdm7zcY1K+LjxgNLF4XaL2+amh5xcB1Gx3YWaopfM+Juv4Xitju0Qnxmd/6Xf8sQ7FfSW2+7Yw7/Fzfobj3FBVR1VrOmEkrE2A9wGcXKUr1f4b7AN45I/PxpiHr1UZ5H344s1dobVMG4otU0oM3wjRDIC3RLxCQC/BtPQqOjq/8yaAuyLizRMc+uXesQmyB9OtW9fBwgo53/O3YLoj9+y+HkxCdSvhZFT11P6CPgSw4xZqsF2NTTtFxd91PC0+d5/bWnavY1YB7I/7RT/lMwL6FZDcblA3TneUsQN7vRimxXhiE9bYGO3zy/Z9HgBwfDFdiMhz33vtwPz7DCfWCoaupwZh42jDFKpwX/uZfb0cBv9YmvVd8DvF9D+OiAaI6tSCKkRry62eo6pzt9RoNptk34LpJejZ1rBbYGHf3wr1dZXPc94MzOjdqd3oqxCnWCge2M1LRGHLw9TQdUcK92wCbQJmbu5Q9zhR4rGblzbSUBfllWvm0ijb1V3ydWlvwXRR12139xZM9zg/d1ob7OYlokSwXasnmDw6eFVxlGCux9bn7aKm9cdkSkREtCBeMyUiIloQkykREdGCmEyJiIgWxGRKRES0ICZTIiKiBf1/rK4qTbzz/N0AAAAASUVORK5CYII=\n", + "text/plain": [ + "<matplotlib.figure.Figure at 0x7f11881b6a10>" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fontsize = 16\n", + "\n", + "fig = plt.figure(figsize=(7, 5))\n", + "ax = fig.add_subplot(111)\n", + "\n", + "ax.set_xlim(binning[0], binning[-1])\n", + "ax.set_ylim(1E-1, 40)\n", + "\n", + "ax.set_xscale('log')\n", + "ax.set_yscale('log')\n", + "\n", + "ax.set_xlabel(r'Deposited Energy [GeV]', fontsize=fontsize)\n", + "ax.set_ylabel(r'Events per 2635 days', fontsize=fontsize)\n", + "\n", + "ax.tick_params(axis='x', labelsize=fontsize, which='both', direction='in')\n", + "ax.tick_params(axis='y', labelsize=fontsize, which='both', direction='in')\n", + "\n", + "ax.step(\n", + " binning, null_exp, drawstyle='steps-pre',\n", + " linewidth=2, linestyle='-', color='black'\n", + ")\n", + "\n", + "ax.step(\n", + " binning, norm_conv_p, drawstyle='steps-pre', label=r'$+40\\%\\:N_{\\rm{atm}}$',\n", + " linewidth=2, linestyle='--', color='blue'\n", + ")\n", + "ax.step(\n", + " binning, norm_conv_n, drawstyle='steps-pre', label=r'$-40\\%\\:N_{\\rm{atm}}$',\n", + " linewidth=2, linestyle=':', color='red'\n", + ")\n", + "\n", + "ax.legend(loc='upper right', prop={'size': fontsize})#, bbox_to_anchor=(0.96, 0.96), markerfirst=False)\n", + "\n", + "ax.text(\n", + " 8E4, 0.2, r'${\\rm\\bf IceCube\\:Preliminary}$', fontsize=fontsize, color='r'\n", + ")\n", + "\n", + "matplotlib2tikz.save(\"plots/tex/syst_convnorm.tex\")" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAFDCAYAAABhmOTuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3U+MI1d+H/Dvz3YcYA9aTiuSMRQMr9i62MACWpKji+McPGzBmJkWkF32TILYyCHbRTkwkEN2m5pTUzlkxN49G8tqAwESB4imuXuYlhRsmrOX+GC4ScoBHNjAblNeBMMBLKub0sGLOFj/cnivikU2yWY1i2QV+f0AhSaLxapHDoc/vn+/J6oKIiIiur5fWHYBiIiIko7BlIiIaEYMpkRERDNiMCUiIpoRgykREdGMfmnZBVikL33pS/qzn/0MAHDz5k2k0+kll4iIiJap1Wr9raq+NOt51iqY/sZv/Aaazeayi0FERDEhIj+N4jxs5iUiIpoRgykREdGMGEyJiIhmxGBKREQ0IwZTIiKiGTGYEhERzYjBlIiIaEYMpkRERDNiMCUiIpoRgykREdGMGEyJiIhmtFbBtNvtQkQgIqhUKssuDhERrYi1SnSfTqfR7XaXXQwiIloxa1UzJSJapkajga2tLYgIGo3GwGP1eh03btxALpdDvV5fUgkX7+DgAAcHB9jZ2cHBwUGo55ZKpTmVKry1qpkSES1ToVBAp9NBKpVCtVpFoVDwHysWi+h0Otjb21tiCRerXC6jWq3693O5HABM9R6Uy+VYLanJmikR0YKVSiU0Gg30er2B/alUakklWrxer4cXX3xxYF+pVMKjR4+ufG6n05lXsa6NwZSIaMEKhQIymQxc1/X39Xo9bGxsLLFUi3V+fo5yuTwQGDc2Ni79wBjFay6PEwZTIqIFabfbyOfzAC7XwprN5kCz76rLZDJotVrIZDL+vpOTkyvfg0ajgfv378+7eKGxz5SIaEGazSYcxwEAOI6DcrmMRqPh96UmKZhOO/inWq2Obb7OZrP+7V6vh8ePH6PVak08X6/Xi2VzOIMpESWGyPjHajXAxim4LjDpu161fzuXA9rt0cft7ppzAUCrZY6NSiqVQrFYvDQQadlc1/UD/iS1Wi3S6+7s7ODp06cDNdVh9XodxWIx0utGhc28RERL8vDhQzQaDXQ6nVj0l/Z6PRwdHS38uuVyGeVyeaCmOswbBR1XosGfaCsun89rnIZSE9H66HQ66PV6lwLG5uYmUqkUnj59ilQqhUajgVKp5DePnpycoFqtotFo+FNJqtUqTk5OcHBwgGw2i16vh0wmg2w26z/fqzkeHR2hVCrh/PwcR0dHqNVqE6/h7c9msxNriVE08wKmtplKpfzaebvdHhlU6/X6wGCl09NTtNttlEolFIvFiWWdRERaqpq/1pODVHVttlwup0REy1CtVsfuz2QyA/uy2ax/++joSPf29lRVNZPJ6MXFhV5cXGitVtOjoyP/uEKhMPD8i4sLf//Z2ZmqqhaLRX//uGsE98/bycmJ1mo1/zWdnZ355VBVPTs7G3iNQbVaLZKyAmhqBPGFzbxERHPU6/VQKpVQLpdH1uYcx5nYD5jJZNC2nbqpVMrfTk5OLtXGgjW3YG0w2IR8fn4+8RqL0uv1sLW1hVKphBs3buDGjRvY3NwceA2NRmNk36zrujg6OkKn08HBwcFU02nmjQOQiIjmKJVKoVarjR2w42VDGqfT6fjNnsGgeOvWrYHHzs/Pr93UOeoa3ijjeUmlUtAruhkdxxk5GGrc/mVizZSIKIYajQYajQZOT09RrVbRbrfRbDb9vL17e3vodDpoNBpwXdcPyO12G51OB/V63b/9+PFjdDodf/+4awDA1tYW6vV6LAZEJQkHIBERxUwul7tyvmUSrpEEUQ1AYs2UiChGvKky81w5ZhHXWDesmRIR0dpizZSIiCgm1iqYdrtdiAhEBJVKZdnFISKiFbFWU2PS6TS63e6yi0FERCtmrWqmRERE88BgSkRENCMGUyIiohkxmBIREc2IwZSIiGhGDKZERAvSaDSwtbUFEUGj0Rh4rF6v48aNG8jlcmuVmejg4AD1eh2u68J13WUX59rWamoMEdEyFQoFdDodf6WY4KosxWIRnU4He3t7SyzhYpXLZWxtbfnvQ7lcRr1en7gkXVyxZkpEtGClUgmNRuPSOpzBNUjXgeu6Az8otra2xi5VF3cMpkREC1YoFJDJZAaaNXu93lotezZqMfKNjY1Lzd9JwWBKRLQg7XYb+bzJqV4qlfDo0SP/sWazOdfFuOPm/Pz80o8Hr2Y+XGNPAgZTIkoOEbMFbW+bfcfH/X2ua/Y5Tn9ft2v2pdODz8/lzP7g2p6VitkXzOEdwdqfzWYT2WwWAOA4Dnq9nl8T8/pSk6JUKk21jQuMvV4P5+fnA/u84Dq8PwnWagBSq3X5/6GnVuv/v3NdoFQaf57gqnW5HDCitQIAsLtrzuVdOz9hkZ87d4APPxz/OBGtllQqhWKxeGkg0rK5rgsn+CNkjFn7Nkf9cPCCaBKbu1kzvYa7d+/6q8+02+N/rR4euv5x+XxuwhmP8fnn/y76ghKtGtXBX7OAqZGqmhqqx3HMvuBUi3Ta7Bte7KLVMvtzgf+jlYrZF6yZ5ib9H76ehw8f+gt1xyGA9Ho9HB0dLeRaGxsbl2qt3v0k1dA9a7U4uIhE+mKbAHJ28yqn+wAqdnvX7ssCaNktWDl9BiBtt+cA1unfgmjddDod9Ho9v5nXs7m5iVQqhadPnyKVSqHRaKBUKqFarSKVSuHk5ATVahWNRgPlchnVahXVahUnJyc4ODhANptFr9dDJpNBNpv1n+/VHI+OjlAqlXB+fo6joyPUarWJ1/D2Z7NZZDKZsa+nNKn5LsC7xig3btzAxcWFf7/RaPivbVGiWhwcqro2W87+rkVge2L33Qvs27X7aoF9N+2+Z/b+nTt3VLNZ8zu52VTf/r7Zt7/f39dsmn3ZrA64eVPVnhuoaf9n9+A2/LRxxwGqtZoSUQxVq9Wx+zOZzMC+bOA//dHRke7t7amqaiaT0YuLC724uNBaraZHR0f+cYVCYeD5FxcX/v6zszNVVS0Wi/7+cdfIDn/hzNHe3t7Aaxi+vwgAmhpBfFnLZt7gG7B97x4A4PjJE3+fa3/RObu7/r7us2cAgPTNm1BVfPjhh+Gah1QvD2DodgFVPL+ivO12y28ulnGdvlap5ODu3btTvAtEtAi9Xg+lUgnlcnlkbc5xnIlJCjKZjD+NJJVK+dvJycmlmmOn0/FvB2uDwSbkUYN7gtdYpGq1ik6ng3q9joODA2xubiYyYQOwZgOQkMsBzebgvuAIQI/jDI4CBPr9LXPw33/1V/E7/6cEByUc2n33ABzb7a3AsaYEgmBIfQJg224fAPjoo7kUk4iuIZVKoVarjR2w42VDGqfT6fhNw8GgeOvWrYHHzs/PJzbLTjLqGo1GYyEDo1Yl49Na1kzj5ne++lUAgFur+TXh4ydPAADb9+4N1KQ942rXwK7diCjJGo0GGo0GTk9PUa1W0W630Ww2/by9e3t76HQ6aDQacF3XD8jtdtuv7Xm3Hz9+jE6n4+8fdw3AZCGq1+uxGBCVJGs1ACmfz2tzuGa6YhwxoxddvXpoOxHFUy6XQyuCea3LvkYSRDUAiTXTFeOiBBfTjbIjovjxpsrMc+WYRVxj3axXn+kacG0TL+ulRMlUKBQGposk9RrrZqWCqYg4MMOcFz8sLSZKMM28DKZERIuzMs28NpCWALDXHJfTJnrpR0dtwYHLXsrF4MaZNkREk61MMFVVFyYp0Vq7if+MmzDp0oJzU2dJe8ipNkREk8UumIpIUURGTroSkT37uGNrojSki3+NLl4BAjNRmwAUeWQhdr9gHwK1f2EHLGUBKNpoBo4jIqKrxSaYikhBRPZgvtkvJXK0AbatqnVbC90UkWSmypinO3cAm6XJ23J2Mnar2fT3Vfb3AQCV/X1/X8tOG8pls/6+Z0jjGdJjL0dERDGcZ2qDZkpVS0P7L1T1RuB+AUBZVbcC+2oAjlR15FLt6zDPNHJe52vMPidERFFYq3mmIpIdsfscQHwWAVxR3qo2REQ0XlKmxmzABM+gHgCISEpVe7bJN2/3na/z9JgoXZWEn4iIElIzhelDHZ7y4gXXDQCwfak5VS2NC6TdbndghKu3VYIrvNCQmt2IiGicpNRMeyP2ecH18npCY6TTaXS73WhKtCZq/mwjDp4mIhpn6mAqIq/bm96CeQ4AqOp3oy7UCOe4PMI3Za8/KtBSRBx/UTh3qeUgIoqzMM28bwPIqOoXAJ4CeBHAUxH51lxKFmCbbYeD5gaAkaN2KToOanDYzEtENFGYYHqiqj8QkVcB5FT1oap+DOCTOZVtmDs0r3QL7Mybu0M4OIQDEZOSMGhcekIRwA1UZF2XqQmJaLWFCabeEgMFAMF1eyKZgCgiWZu0oQjgvs125E+JUdUygIzNgLQH4ExVQ60fFByAxEFH0/rQv9VutwYGbk1SKjn+caVSv7+VqQmJaBVNnbRBRL4N019aBeCo6o9E5DaAV1X1j+ZYxsgwaUN47+bzuNdqoQvgrcB+71MTDKlPAGzb7QO7bxemt9UFULLPYv4HIoqLhSdtUNXvwPRTlmwg/QZMOldaYfu/8ivIAdi+d28gRaEnuG/73j0AwPGTJ/4+t2Za4p3d3WUUn4hoIcLUTB+p6sM5l2euWDNdrrSYaUldZU4lIoqHqGqmYeaZlkXkDMBjO6KXKBSzmg0QUTc7EVFshBmAtGX7Rm+JyDcC804TgwOQlqtrNyKiVTN1zVRVnwb/isjrIvJDmFVaEjEAiRmQlov1UiJaVVPXTL2aqA2ijwH8CGaOacvWVL8+pzISERHFWpg+07qIeJWKAwC7qvq5vf8xAIjIN5NSS6Vl4OAvIlpNYRPdv+018w4Tkd+OoDy0wpp+svzWUstBRBS1UKN5hwOpnWt6oao/ApADvyVpghy4xCwRraYwSRu+DwAi8hVvgwmeO/bx79igGlsczbtcOTSRY1MvEa2gMAOQviYi5zABtAWgDeAMCaqNptNpPzMPg+nitZFDGzmIAE5gedRWa3LS/FbgE+Y4g48xcT4RxUGYeaYlmNViXoTJzbuhqr+Ixa0aQ4n3IfZRgUKQPtz2Wwmc/K9DIWhicFmaZ0hDIdjOf9U/Nn/4W1AIdu36qkycT0RxEHYJNi9w3gjs/3KE5aEVdufOHwJ41977IPDIX9m/bZjU+d723O7/i8Cxf2L/luZUSiKi8EKN5hWRr6vqDwB8EridmU/RaNV8+GF/ObeK3YJyGJ3QYVSaDRfAczm297ZnLxwR0QzCBNMOgEMRaajqUxH5HyJyCODxnMpGNNGxvygccyoR0XKFSSf4MYB84P6bIvI1u59o4VgvJaK4CNNneomqfiwiL0RVmHnj1JjV8hYGFywnIlqWsBmQRjkE8CCC88wdE90TEdE8jAymIvJlABe4ujNKpjiGaE740SOieBjZzGsT2NdV9Re9DcCbAF4b2pcHcH+B5SXyKQQKWXYxiIgm9pnuDt3/cmCeKQBAVdswNVgiIqK1NbbPNLC82lVSEZWFKBSxzbxs7CWiZQszmvcNb4Fwj73/RrRFIgrH5Ol1/JHa5vakfL8S2Fr+/pdfPl32SyGihAqzasw7AA5E5DMRORWRz2BG8v7HuZWOaIKXXjrFE2xDIbiHn/r7d/GXUAhq6GfTv4kuFIJnSA+co4nfhUKQRQuffnprYWUnotUSap6pqr4JYAsmm9t9Vb2lql/MpWREV/ibv7mF7Xvm9vGTP/BXBHJrvwcAcHYBVbN1n5nj0jfhH6eqyGW/tKTSE9EqEdX16XFKp9P6/LlJnr6/v8/EDeQTOyh4jf47EBEAEWmpav7qIyeLImlDYjBpA42z76fdr0w4iohotLUKpkTjVPyl4SrLLAYRJdRMuXmJVkUFDKNEdH0MpkQwS5a/e+VRRESjsZmXCACQXXYBiCjBxtZMReQrI/btishPROTnIvJjEfnn8ywc0aJk4SILd9nFIKKEmlQz3RGRM1X9AQCIyLcBbAIoA+jZ278vIuIdQ5RULX/de86NIaLwJgXTQwBnALxA2VHV7wQefwrAFZFHgWOIEqllm3lzSy4HESXTpAFINwCkAs29436yd6IsENEy5NFCHv08va1W/zHHGZ/nNzcUfYOP3b272NdARMszKZgW7d9z+3fDe0BEXggcl5hVY7rdrp/gnNmPaNCHAIBnSEMh2M5/1f+s5A9/CwrBbqBP9R6OoRBU2v9sIHF+cI3Vjz5aygshoiWY1MzbAfA2gPsiJtmaiHxFVf8aQE9EqgA+A5IzaoMZkGicO3f+EB99dC+w5y8Ct//E/i3ZLeh/zrVcRJQM18rNKyJfBkKteRoL+Xxem83msotBa4C5fomSYam5eVX186GmXiIKeIJte+t4qeUgosWYmAFJRL4pIu+LyL8P7NsVkX8AcCEi/23uJSRKoG18gG18sOxiENGCTEra8B5Mn6kA+Jd2QfDX0e84+n0Av2CnxhBRwLbdiGg9TGzmDbYji0gWQHWobdm1QZeIAlgnJVovk4LpT4J3VLUtIvURx/1ttEUiWgW7yy4AES3QpGB6AzD9pqr6R3bfY+9BEXlBVb8A7KQ6IvLtYubBgUSUIJOCqdeE+w0RaajqX3tTYezUmAsRqQGoLaKgREni+vNRnaWWg4gWY2wwtYHzHbtdekxEtgA0kzbXlGgRXNvMy1BKtB6uNc/UNvE+jbowRKuiZBODMZgSrYdJU2O+ZueYvhDY915gjunPg/NPiYiI1tWkpA15mGbcLwB/PdPPYFapeg3ALQAvici35l7KiDDRPS3KTXRxE92BJPgiLYxbfUbEDRyXu/T4yy+fLvslEdEEk4LpZ0Prl/ZU9Tuq+rGqfqKqbVV9BwkazZtOp6GqUFUGU5qrLl5BF68M7Gvid6EQZNFf320fFSgE+4G0g1n8HRSCZmB11U8/vTX/QhPRtU0Kpm8M3T8fedT4dU6J1tedO8DNm/6PN1VFLvslAECraRLgqwKVfXN4ZT/nH9dq/jEAIJdlonyipBi7aoyd/vIUwPcANABk7ENNVf3CLhpeAJBS1e/Ov6iz46oxlERdSQMA0srlA4miFtWqMWNrpnbKy22YvtMOgBO7XYjIzwG0kKBASpRUaTxHGs+XXQwimuCq3LyfwyS7f1tEMgBehWnu7anqJwsoH9HaS9u/rJcSxdfU80xVtQNTQyWiBWKdlCj+rlrP9Fsi8pmdU/pDEfm1wGNfs/NOuaYp0VwxaydR3I2tmdp5pQ9g0gmeA3gTQFtEfltV/5eqfiwinwP4MYB/sZDSEq2hGrxBc8ynRBRXk5p5bw2NcPq+iJQBPBaR91T1RzBJHIhojhwc2lvuUstBRONNCqaXUq6oag/Am7Z5F0Bg9jkRzYVjm3gZSonia2KfqcfOKfXZzEc3ANyPvkhEFHQIB4ds4iWKtUnzTL9ja6DfBnAmIq8PPf59AJ8gQekEiYiI5mFizdTWQF0Ar6nqn494vAFTQyWiObmHY9zD8VDS/HEJ8wERJ3CcM/IYJs4nitaVzbyq+vmkBA1cHJxovo7xFo7x1sA+hUCHGoWeYBsKwT381N+3i7+EQlALNBPfRBf/6dPKXMtMtG6m6jMloiW6cwe4d28gab7HS5ivCmzfM/uOn/yBf5xb+z0AgLMbSK6PCu7io2W8EqKVNTbR/SpionsioGtG4iO9Rv/3icaJKtH91OkEiWg1eKusMpQSRWetmnm73a4/MIOLgxMRUVSmrpmKyDdh1jK9NKo3KdLpNLpdrr1B645dHURRC9PM+yZG/C8UkRdU9YvoikRE89T0R/YygRlRVMI0876P0UuwMTULUYLk0EYO7WUXg2ilhKmZbgF4T0Q6AHp2nwC4DeC7UReMiOYjZxuYWC8lik6YYJoHcACzHFtQKrriENG8tZFbdhGIVk6YYFpW1afDO21NlYiIaG1NHUxV9amIvACzUkxHVX8kIq+r6sfzKx4RRW0fFXurMuEoIgpj6gFIInIbQBtmVG/W7v5ERL4+j4IR0XxU8C4qeNfOuc5dkTQ/F0ia7zJpPtEYYUbzbqnqa6p6H8DHgJ/knkuwESXIf33tNb9OmsXfQSFoDvWjPkMaCsFN/L2/r4b/AoVgN7BM+T0cM2k+EcIF0z8bs59ZyYgS5F/9+Meo2ET4reYfAwBy2cGk+emb5tjusx/6SfOd3V8HALi1YNL8PSbNJ0K4YPqGiPyava0AICJfAfBGxGUiokXJ5UxUbA1NlOl2bVRN9/e5rtnn9KeWd/FXOF5QUYniLMxo3kcAWiKiAHoifuvu7chLRUSJ4K2yyuYpWndhRvN+DuA1ESkCeBVmRO/351YyIiKihAi9BJuq1udRECJKItZJiYCQS7CJyLdF5FxEfi4in4nIv5lXwYgo/hQC5YB+olBLsO0CuAVgBybh/SYAR0REVf9oTuUjIiKKvTDNvJt2jqnnEwANEfl2xGUiooQQ28zLxl5ad2GaecelOeFaTkREtNbCBNMbY/Z/2bshIt+crThERETJE6aZ900RKWNwgfAMgI6IlGDSCuYAsP+UaE08wba9xdQNtN7C1EwzAN4B4Aa24P0auN4w0VrZxgdI409tInwHjriACFxx/ET4aekCIuj6CfMFIi20JAeIICct/9iKVHDyj//psl8WUWgzr2caxLVNidbL6Usv4T98+mlk56vgXQRy6xMlhqiuxjg8EXFgmqAzABqqeimw5/N5bTabCy8bEU2nYtOUVlbke4niT0Raqpqf9TyhMyDFkYhkYKbuuPb+Ecx8WCJKkHft38oyC0F0DSsRTAEUAZwF7mfHHUhEccb/upRMsQumNpH+LVUtj3hsD6YpdwMAvJoogBcxOMoYIpJS1d6ci0tEEcoGFh4nSpLYBFMRKcD8LN3CUGC0j1cBnKhqw7svIsVA4v2NhRWWiOaiBa/rin2mlCxhcvO+bm96gc4BAFX9bhQFsUGyISIvAkiNOMQZqq2eACgDqAP4bOjYDdZKiZKnZZt5c0suB1FYYeaZvg0go6pfAHgK07T6VES+NZeSBYjIqI6UcwAFe7sO+/9PRFIAGvMuExFFL48W8pyuTgkUppn3RFV/ICKvAsip6i3AH0k7bxswwTOoZ6+fUtWOiLQCTcWX+luJiIjmJUwwvbB/CzA1Qc8iOjdSuNwn6gXXDQC9wGAk1kqJiGihwjTz5kTkGzC1vu8BgIjcxmIG/ozq//SuO1xjHavb7QbSmfW3SqUSSSGJaDbPkMYzpJddDKLQpq6Zqup37ALhJVX9kQ2sGfRrrPN0jsuDklK2XFMPNEqn0+h2u1GWi4gilMbzZReB6FpCTY1R1cPA3Y6qfj/i8oy7bltEhoPmBtikS7RSvDopf/JS0kzdzDtirdKeiNwWka9HXKZxXJvQwbMFs1INEa2I5wCeowZXHEAEjrj+ijLbcgyI4HigmwbwDvCOEwGOZRsQwenLLy/7JdGaCNNnOtDMqqqf2FVkIukzFZGszXBUBHBfRPaCU2LsHNOMiBTtcWeBhA1EtALu3LkT2bla+CeRrmhDNMnEVWNsH2kOwA2YKSftoUMyAJqq+vtzK2GE0um0Pn9u+mT29/c58IhohYldgWZVVsai+VjIqjG2j/RQRL4H4GMAj4cO6ajqx7MWYlE4AIlonewuuwC0RqYdgFQGUFjUgCMiolntYubKBtHUpgqmqvo5gJGBVES+oqp/HWWhiIhm5aJkbzlLLQeth9CrxojIC0O7ygAS0WdKROvDtc28DKW0CGGmxuyKyD/AJGnoBf4m5rMazIDEwUdEq60EFyWuj0oLEqZmugnghm3y9YnIe9EWaX44AImIiOYhzDzTk+FAaj2KqjBERFG5iS5uYjgfd2sgucPg5gaOy4085uWXT5f9siimwgRTFZGvjNjP8edEFDtdvIIuXhnY18TvQiHIBtZM3UcFCsE+jv19WfwdFILm0DLln356a76FpsQKuzj4iYj8WERO7dYEUJ1T2YiIru/OHeDmTaiqv+WyXwIAtJqAqtkq++bwyn7OP67V/GMAQC7bP44r2tAkEzMgDRxoAucjDC6HJgD2VPXNOZQtcvl8XpvN5rKLQURJZDMqgRmVVspCMiANKdtcvMMF+WzWQiyKN5oXYDpBIgqHK9rQJFM386rqUxF5QUS+KSK/DQAi8nrS0gl6zTgMpEQUxnO7UTgvv3yKilQAEVQCA7xy0gJE0Boa7NWVNCCCtH+sO7CKUFwHgYWZZ3obJtH9mzBJ7wHgkwUuwUZEtEQ1cNXH8KIetBXXQWBh+kzfU9V37O3bXpOviHwjKTl72WdKRNflislP4ygTQYQRZVfzPLqtl9Fn+mdj9rM3nohWnoNDe4vBNIz+NKTcxOOmcc+fvrQ987miFiaYvmEj+E9hA6idd/oGgB9EXzQiovhwbBMvQ2k4LX/1ntnrXcd4K7JzRS1MMH0EoCUiCqDnjYoFcDvyUhERxcyhTUPOYBpOyw6xmb1eihjXS8MF03dU9TURKQJ4FWZh8ET0lXo4NYaIaLHy9udHFHXJ+NZLQ84zFZEzAI9V9Yt5FWiemOieiK4rzv118RZdM2+chQmmO6r6fRG5LSIZAGeq+qN5FYyIKE7i3F+3PuL73odJ2vB9+/epqh4C6IjID0XkW3MrHRFRTBwDgVT4q+vu3buB1XPgZ1MIJlY4lm1ABNv+cQ4ccQERuOL4x6WlCwXwLKKyKQQKufrAJQiTtOF176+IfA9AB8DnMIkciIhW2lvo99mtsrc/+ghPIj7n//7l34z4jPETZtWYuoj8GMARgJ/ALBR+n029RLQ+dOaaGkTQFcHdu3eX/WJG2rabSb0Kf9kcb/UcVWBbjwFVHPsr8rhw1QFU4ajrH9fVNKCKrf/7J5GUTWzdNI7C9JkCwNujkt0TEa26O3fu4KOPojvfR1GeLELbtl66Dk3aUQqTTjAxaQPHYTpBIooDb4retN+/ixTnlebinE4w9ACkJPPmmYoI55gS0RI17UZhPME2nsR0atLYZl4R+SaAFMxi4E1V/XMR1WDrAAAS40lEQVS7/8sACgAyALYA/ERV/+0CyjozzjMlojho2mxK8PPWxseun+PJmXjcMmzjg2UXYayxzbwi8g8Asl4QHXNMCsBnqvqLcypfpNjMS0SxsG5tqRHZtmU7jrBsi1g1xp0USAFAVXsicjjpGCIiGpSzTbzxq5cCLnYBxLFeihjXSycH0594N0RkFyZP8S6Asqp+d9RxRER0tXYkad/no2TT/8UxmMIG+jiaFEz9NBM249GhiGAokA4cR0RESVeyf+MXTncxc2vs3EwKpqMapUfVQuPXsE5EFGP7qNhblQlHLcfNZRdgAjfGgX5SMH0QWLPUszVi3wMAw7VVIiIao4J3/VtxE+f5DnHuz50UTHMAbozYvzl0/9XoikNEtPoqQ3/jpGvrpukll2OUkp22k7RgeqCq71x1AhF5L8LyEBGtvPjWS4FXbN2U/XfhTMqAVJvyHNMet3TMgERE8ZC1G4VxE13cjGlD9Niaqap+Ms0Jpj0uDpgBiYjiIOtnGaIwunjF3opfvTnMEmxERBSBFvJoIT+wjFtX0oAI0v4ybi5ccQAROOL6x23LMSCCY/84sxzcyy+fRlK2JnJoxnQebBfxHSAVdgk2IiKa0ckv/yY2/v5nkZ1PIcCn5tascmjPfI55iW+9NMQSbKuAuXmJaCVFmE83JybJYUvjVzudx9J1i8jNS0RECSC2rhZFiGnHeEWbOC9bx2BKREQB8W3mjfPSdQymRETk2192ASaIc38ugykRUcI9wba9dTzzuSozn2F+4rx0HafGEBEl3DY+QBp/aqfJOHDEBUTgiuNPqUlLFxBBd2BKTQstyQEiyEkLIkAF+6jEtH7aRi62y9cxmBIRJdzpSy+hgr+N5FzvooI/fOluJOdaJ5waQ0REiVCRivmrlcjOyakxRES0VuK8dN1aNfMy0T0RUXJVEMcwaqxVzZSJ7omIkiu+9dI1C6ZERJRk8V22jsGUiIgSIc5L1zGYEhFRIrTgDbqN3ywUBlMiIkqElm3mjWPaBgZTIiJKhLxNJBi/eumaTY0hIiKaBwZTIiKiGTGYRqFeB3Z2gBs3TEbpXM7c7/Wiv9bBAbC1ZbadHaBUAlwX6HTMdcOW2ytzvR5dGRuN0e/H1pa5nctFd712G9jcNNc5OOjvd11z/UYjmutEfT4iCu0Z0niG9LKLMRL7TKNQLJpta8t82ZZKgONc/bwwvACVzwNHR0AqNfjY5mb4cxaL5u/OTjRl9BQKZvPej4cP+9cCTNDb2QH29oBqdbZrZbPm/Rj+IXF2Zn7MRPWDJurzEVFoaTxfdhHGYs00CXq9fsA7ORkMpIAJXLXa4st1Xd4PjWBNMmrVKnBxMRjE43Q+Igotbbc4YjBNgnLZBNSHD8cf4ziXg2xcBcvZbi/mOnE8HxGF8txuccRguijlsqlden+DTau9nmka9vpAt7YG+xS9frpMZvI1PvnE/A32I5bL/XN4fZijaoQnJ+bYctlc33vetGUMI9hUms0O9t16zeSbm6af8jrX7nQu96MGX3+9bs59cNDvzwXM/eF905zPdUefz1Ovm/0HB/3yB/teJ71+1+1fO5frv3euWfwZW1vh33+ixKrZLYZUdW22XC6nc1UoqAKqtdrl/alU//7enjnu4sLcz2RUs9n+462WefzszNwHzHZyMn1ZTk7Mc/b2Lu+rVvv7jo7MPscZfH4qZcrtuaqMo3jvx9HR4H7Hufw+1WpmX7HYv753+6pre/eDr2vUPu/1Fwr999677qh9wXKHPV/wudms2ecpFgf//a96/d61vfuq5rnD/2ZEK66GXa1hN9JzAmhqBPGFNdN5a7dNbeP+/f6+atUMaEmlTK2k0wEePOg/nrXJnL1+UK9GGmbwy8ZGuHIOD+ApFEy5O53pyjhJrWZqW17tEjCvPzhIyyvvrVvm78WFGVg067WDvGtkMv0mW++9HbWv07n++YLPPTw0NX+PV5sM1k7HvX7AvN5Cof9eAMCjR5dbD4hWnINDODhcdjFG4mjeefO+/Ib724a/dE9PB5tfi0XgxRfN7ULBNOudnk4eALOz0/8CnpX35d7pTFfGSXZ2ph/dPNyUPeu1RxnV93md0dCTzhfkNWV7fd+TgvS4pvxy2QTfatVsnc7Vzf5EK8axTbxxTHfPYDpvV9VyvFrWrVtmqsgo1Srw+LEJqOOmknQ6V9dGw9Rsz8/N3+AX9qQyzss070+c9Xqm1p9KAU+fmr+uG36+aqFg3guvH3nSYDSiFXUI86M8jsGUzbzz5jXRNRqDwezgwDQBe1+Sw02WjUa/JpZK9WucW1uXA3OnY5pQg4F2VG0p2NQ47Oysf7vXM9cvFk0wnaaM87LMa0fBayp/8KD/bxJ8r8PwAmiz2f+RQUSxsFY10263CxEBAOzv76NSqURz4nodeP998yUHmC/+kxPTV5ZK9UfK3r5tki6kUqam5X0htlr9UbReTXBzc7AmViiYfjRvRChgzrOxYY4dDpSZjAmutZoJjqlUv1/00SPz1zu/l2DB64Nrt80Xd/D605TR02iY6457P0Yd75Xp0SMTfKa9drvdL/f775vHM5nBfdmseZ+8ffW6eX4+3/8BUqv1+z69fVGcr1g0Tdzvvw989plpmn7woN/se35ujpv0+j3ej5tZE10QJdQ9HNtb20stxyhiBjOth3w+r03vC54oKYL9o1tbk1sYiFaZrQwhwrglIi1VzV995GRs5iWKs3K5PzjKmwdMtKaO7RZHa9XMS5Q4Dx70m7I3N02TPNGaesv+jWN7KoMpUZxls2zWJUoABlMiIkqIONZJDQZTIiJKBIX4t+KGA5CIiIhmxJopERElgtgaafzqpayZEhFRwhzLNiCCbRGICEQcOGKWJXTFgYiZkpqWLiCCrn+cQKSFluQAEeSkFVmZGEyJiCgRXnrpdNlFGIsZkIiIaC2ZhErMgERERHRt+6hEdi4GUyIiWksVvBvZuRhMiYhoLVUiPBeDKRERraXo6qUMpkREtLaykZ2JwZSIiNZSFm5k52IwJSKitdTCzDNifAymRES0llps5iUiIppNHkwnSEREFBsMpkRERDNaqWAqIo6IRNcITkREK+sZ0pGda2WCqYg4AEoANpZdFiIiir80nkd2rpUJpqrqAkjskjCVSmXZRSC6Nn5+KYnSA39ms9Al2ESkCOCWqpZHPLYHoANbs7TBMez5awCOVLUx6vE4L8EmIlin5fBotfDzS0kkZg02qKrMeq6F1ExFpGCDZQlAasTjVQBtVa3bILppAy8REdGc1CI700KCqao2VPUAQHvMIc5QbfIEJvAC8AcW7Y3YCmHK0e12Q5V72qaraY5bx2awZb7meV07ivNe9xxhnsfP7mz42Y32HGGft6jPby3CnsFFN/NWAaRUNRgoswCequqNoX2tsFXvq5p5RUTDvN5pm66mOe6qY1axmWyZr2le147ivNc9R5jnLfKzG7ZsScDPbrTnCPu8hX1+RSBIUDPvFTYAnA/t6wGAiFxqEh7HNgvnAexwegwREV3FibCZ95ciO9P1pXB5OosXXDdgA+tVVLUOoH7FYf9PRP7RiP3PAYxqA06LyDRtw9Mcd9Ux014rSZb5muZ17SjOe91zhHneIj+7YcuWBPzsRnuOsM9b5Of310OUa6w4BNNRwdILrsM11pmo6i9Heb55s83WZQAZAFDVcX3ORLFjW5buw4zSx7juF6K4sZ/dC9jPLkxMuj3pOzgOwfQcl0f4pgBAVaeqla6wPIBPALijphMRxdyRqm7ZL6ZDAAymlBT5YD+qiDhXVWaWHkxVtS0iw0FzA/yPBwCPbPM1UaLYkfZejbQHYGe5JSKaXrAVRUSK0+Q9WHowtVxbYC9wbCHKCUBLNkOyig37pZQFUFfVzvDziebtmp/frH08C6AAoMFuClqGCJIFTZWidiHBNPAfqggTIM4Q+M+lqmU7b7QI0z94tgo1skAg3EK/7T34eBXAifcrSESqwR8V3j+siHQAHAHILarsRLN+fgFs2JanDoAWgM0FFZ0ois+vl/N9qlbShc4zXVej5tfa/RdD82sLAMq2n6kIIGOTXVw6lmhRrvn5dQBserUBO8d75rl8RGFd5/Mb2NdS1akqMXGYZ7qWxsyFPYepwQNmlHPdHptBgpP40+qZ4vP7GHYUuh2AxC4Kio0pPr+ezLTnjEuf6TqamKxCVRsiUrT/6LfAARwUL1d9fnsi8r6toabAzy/Fy5WfX7tv6h+BDKbLc2WyikDbfeL7j2nlhPn8EsXNVMmCpm3iBdjMu0wLS1ZBNAf8/FKSRf75ZTBdHiaroCTj55eSLPLPL4PpkthpQUxWQYnEzy8l2Tw+vwymy+UOLYK+UskqaOXx80tJFunnl/NM5yiQrKIE86vnEYYywQQycGRgBm1cmbaKaBH4+aUkW/Tnl8GUiIhoRmzmJSIimhGDKRER0YwYTImIiGbEYEpERDQjBlMiIqIZMZgSERHNiMGUiIhoRgymRAli1wZdSyKSEZEju9jzMsvhiEhtKHsOrTkuwUaJYRdJLwHYA9AG8H7g4U0AR6oai9ywIlIDAFUtRXS+FIBDAEUAMuG4q94jB0ApwZmKTlX1ILjDvuYygLPAbu9zkLlqKTgbFB8CyAJwAZSDyc7tmqw1e86yqrr2mgym5GMwpcRQ1Q6Asv3ye3/El+qRiGRiEiiOMLSUk4g41y2b/XLfEZGJKcumeI9qMCnWVoJ9nSUAO0MBsAATAK+sxapqXUQaAC5gfpD1hh53RWRTVcvRlp5WCYMprQxV3RERFZGGDSrLLMuoGvLUCw3Pi6q2RSS/7HJEwdYODwG8OiIANkSkPfqZl6lqzwbUEoZWDrGtAqcRFJlWGPtMadU0MEVtZJ5EJCUiBZto27tfQ3/x4WWUyQncbS6rHBGrwSQuH7f+ZNgVQKoAiiP6pe9f1VRMxJoprZoTmNoFAL9W8RCmZnELwImttRRgvjw79jnn9vGzYFOsfb5jjwNMH9xB4HGvbzKF/hJOPVsGb18BJpBm7fH+6hTjyhc4fxX9WtEsi277tWJv1Yyh98ALPFtj3oNJ72ET5j18YB/z+hRLgecE+zPP0f/Bs2NrywV7joMQzaneiiAj2TJmrnodQ8f3YP69g83jazvoi0JQVW7cErXBfDHvjXnMMR/rgWMzQ/dT9nYRJnAEn38UPDeAlne8vV8AUAtcqxh4rAggGzjuZOixozGvZVz5Trzz2fuZ4Gub4j06gRmIVBv3PO898MoAEzg0RBkd7z7MAB7v9bfGvReBfcPvvXPFa8p4/zZeOYPv/5TvycjXEdhXDZbLljM7qSzcuKkqm3lp5aRga3De1AUd7D9tA7gfuD/ct1qDqb14NTdooBlRTU3GsbWcDoCqnSrhjRqduq92UvlsE3FeA2svavh+4BNVPVAzotiv7Q01Y/aC5/Zeq3fMFO/hOYCOqvZUtR0obxb996IDYKCf1r6PCDSFFwA8nvaFBf5NLjWd23+Pqoic2EFp2Sk/C4D598945YIJpFP3vdL6YjMvrZpN9PsEMwB6XlC0TjE54HXQb9bLYmhErtWDCXQNESnBNDXW7ICX2yHKOql8+SvKGYoOjuotAAj2AU5qPp7mPRz1/HrgOuMCZRXmh8sOTA0xbDN2A6ZZemCEtPab0BVmGpDXjHzlZ0FVO95AJBEpj3ltRJcwmNKquY9+QPNqW8HRmVfNQ82g/wXawehBQykAHREp2tqoV8uq4nJ/2yW2H8+r2Y4sn/3Sz1x+diTCnPc67yFg5rdu2BrhuY6eb/sYpmYfrMWGUQLQEpFpAnGY11GDae5vIURtmdYbm3lpZYjIEYBHXrOcDXRe8PKOSQWa8IDLgcVvErXPTw09vwigbpsLM0M1nfcxXidwrYxtEh1bPvul3wmWdajc12Kbb7emPX7K93CUW6rqqmpdx4yEtQHwMYDqdZpS7b/BLoCnwfLZMhYQqFWGeR2B8uauUVumNcWaKSVGYIRoBsADET8R0IswNb2aXp7feRvAQxHx5wkOfbl3bIDswTTrnuhgYoVc4PkbMM2RO/axHkxA9TLhZFT1wH5BlwHkvUQNtqmxaaeoBJuOJ5XPe8yrLXv9mEcAdkd90U94j4B+BiSvGdQrpzfK2IXtL4apMVZtwBpZRvv8kn2dewDcQJnOROQi8Fo7MP8+w4G1hqH+1DBsOdowiSq8a39mr5fD4I+lqz4LQQeY/OOIaICoTkyoQrSyvOw5qjp1TY2uZoPsA5hWgp6tDXsJFnaDtdBAU/k0583AjN6d2Iy+CHEqC8UDm3mJKGoFmBy63kjhng2gTcDMzR1qHidKPDbz0loaaqK8ds5cusw2dTuBJu0NmCbqE9vcvQHTPM73nVYGm3mJKBFs02oV40cHL6ocDkx/7Mm0TdS0+hhMiYiIZsQ+UyIiohkxmBIREc2IwZSIiGhGDKZEREQzYjAlIiKa0f8H7RjuCsThG0kAAAAASUVORK5CYII=\n", + "text/plain": [ + "<matplotlib.figure.Figure at 0x7f1179a73a10>" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fontsize = 16\n", + "\n", + "fig = plt.figure(figsize=(7, 5))\n", + "ax = fig.add_subplot(111)\n", + "\n", + "ax.set_xlim(binning[0], binning[-1])\n", + "ax.set_ylim(1E-1, 40)\n", + "\n", + "ax.set_xscale('log')\n", + "ax.set_yscale('log')\n", + "\n", + "ax.set_xlabel(r'Deposited Energy [GeV]', fontsize=fontsize)\n", + "ax.set_ylabel(r'Events per 2635 days', fontsize=fontsize)\n", + "\n", + "ax.tick_params(axis='x', labelsize=fontsize, which='both', direction='in')\n", + "ax.tick_params(axis='y', labelsize=fontsize, which='both', direction='in')\n", + "\n", + "ax.step(\n", + " binning, null_exp, drawstyle='steps-pre',\n", + " linewidth=2, linestyle='-', color='black'\n", + ")\n", + "\n", + "ax.step(\n", + " binning, prompt_p, drawstyle='steps-pre', label=r'$N_{\\rm{prompt}}=2.4$',\n", + " linewidth=2, linestyle='--', color='blue'\n", + ")\n", + "ax.step(\n", + " binning, prompt_n, drawstyle='steps-pre', label=r'$N_{\\rm{prompt}}=0$',\n", + " linewidth=2, linestyle=':', color='red'\n", + ")\n", + "\n", + "ax.legend(loc='upper right', prop={'size': fontsize})#, bbox_to_anchor=(0.96, 0.96), markerfirst=False)\n", + "\n", + "ax.text(\n", + " 8E4, 0.2, r'${\\rm\\bf IceCube\\:Preliminary}$', fontsize=fontsize, color='r'\n", + ")\n", + "\n", + "matplotlib2tikz.save(\"plots/tex/syst_prompt.tex\")" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAFDCAYAAABhmOTuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3c9zI+d5J/DvE2eVlA8ShtpR1UCHlUAfdp2LBHB0iU8eUAdZVNXGoGb3spc1G8pe9rA24dnLQHsZkfEfYDRz3MPOEPZBVLkqIcan5JIB4Ow5Zis5CFNlZUjIl1S0ZT97eN9uNED8aqKB7ga+nyoMgUaj8QLE4OHz9vs+r6gqiIiI6Ob+IOkGEBERZR2DKRER0YIYTImIiBbEYEpERLQgBlMiIqIF/WHSDVilb37zm/ov//IvAIA7d+4gn88n3CIiIkpSp9P5Z1W9vehxNiqYfvvb30a73U66GURElBIi8k9xHIfdvERERAtiMCUiIloQgykREdGCGEyJiIgWxGBKRES0IAZTIiKiBTGYEhERLYjBlIiIaEEMpkRERAtiMCUiIloQgykREdGCNiqY9no9iAhEBPV6PenmEBHRmtioYJrP56GqUFUGUyJKpVarhd3dXYgIWq3W0H3NZhO3bt1CqVRCs9lMqIXxazabqNVqc+13fHwMz/PQ7/eD62mwUavGEBGlXblchud5yOVyODo6QrlcDu6rVCrwPA+Hh4cJtjA+rVYL3W4X5+fnKBQKM/e/vLxErVZDrVZDLpfDycnJXI9bBQZTIqIUqlar2N3dRb/fRy6XC7aHr2dduVxGuVzGixcv0O/353rM1dUVLi8vUxNEfRvVzUtElBXlchmFQgGu6wbb+v0+tra2EmxV8nK5XOoCKcBgSkSUKt1uFzs7OwBMdvro0aPgvna7PdTtu4lc10Wz2QzOn6YFu3mJaK2ITL6v0QAcx1x3XaBanbyv6uB6qQR0u+P3OzgwxwKATsfsu4h2uw3HNtJxHNRqNbRareBcalqCaXXamxdydHQUW9e0n62H2+C6bvB+JYnBlIgopXK5HCqVyrWBSGnQaDRW/pyj3bu7u7uo1WoMpkREcQtnlNM4ziBLnaXTmW+/RbPScR48eIBSqQTP8zb6fGm/38etW7dwdXUVZLq5XI5TY4iIaJjnecH5Ul+xWEShUMD+/j6ePn0KwEwpqVarQXZ4enqKarWKy8tLnJ6eotFoBPtcXFyg2+1if38fnU4HuVwOx8fHKBaL6Pf7KBQKKBaLY485LftMopv38PBw6Fie56VmMBKDKRFRSjSbzbFzSP0g5weScrmMXC6HnZ2dYD5qLpdDsVhEo9FAv98P9gEGARkwA3gKhULQbby7u4vz8/Oxx5xmFd28nueh2+2iUqkgl8vh1VdfHbr/9PR0ZjtXhcGUiChh/X4ftVoNruvi4uLiWqByHAcvXry49rhwlhbuAr68vJyYDZ6fn+PBgwdD28IZXvhxo3Nc49btdtFqtdBsNnF5eYnt7W2Uy2UUi0UAJgM/PT1FpVIBYN6H4+Nj5HI5XFxcoFqtBvcljcGUiChhuVwOjUZjYrY3T6Y4TjjA+ucW7969C8/zgoCVZAGEYrGIYrE4saKT4zhDg4tyuVxqqz9xnikRUcZ0u114nodmsxlcf/LkCTzPC7YDwP7+frCPXwDi8PAQnueh1WrBdd0gSE86Js1HdN6hb2tgZ2dH2+120s0gIqKUEJGOqu7M3nM6ZqZEREQLYjAlIiJaEIMpERHRghhMiYiIFrRRwbTX60FEICKo1+tJN4eIiNbERs0zzefz6PV6STeDiIjWzEZlpkRERMvAYEpERLQgBlMiohRptVrY3d2FiKDVag3d12w2cevWLZRKpaDKUdYdHx+j2WzCdV24/irrEzSbTRwfH8PzPPT7/eB6GjCYEhGlSLlcxv7+frAoeFilUsGDBw/Q6XRSU+B9EbVaDcViEZVKBY7j4OLiYuofCZeXl6jVatje3sabb76JQqGQmiXYGEyJiFKoWq2i1Wqh3+8PbV/mKi6r5rpusBQcYJaDm7W029XVFS4uLnB1dZWqPygYTImIUqhcLgfF6X39fn9oJZgs63a717ZtbW1d69oelcvlUpONhjGYEhGlSLfbxc6OqbterVbx6NGj4L52uz2UyWXZ5eXltT8M/Kx7NBsPc10XzWYzOH+aFhs1z5SIKO3a7XawhqfjOKjVami1WiiXy/A8LzXBtFqtzrXf0dHR2K7pfr+Py8vLoW1+cJ20uLmfrYfb4Lru0JqniVHVjbmUSiUlojUHmEvY+++bbZ9+OtjWaJhtBweDbV98YbbduTP8+GLRbG+3B9sePjTbHj4cbAvff0ONRmPodqVS0XK5PPa+LDs/P9dcLje07eLiQgHo1dXVXMc4PT3VQqGwUDsAtDWG+MJuXiKiFHvw4AFarRY8z1ub86WAyUJHu3P925MyWREZekwul0vN1Bh28xLRelG9vu3s7Po2xzGXsHx+/OM7nevb6nVzCSuV5m3lWJ7nBedLfcViEYVCAfv7+3j69CkAMxe1Wq0GI19PT09RrVZxeXmJ09NTNBqNYJ+Liwt0u13s7++j0+kgl8vh+PgYxWIR/X4fhUIBxWJx7DGnjaxdtJu3WCxe2355eTm1G/vw8HDoMZ7npWcwUhzpbVYu7OYlojQ7OjqauH20O7NYLAbdoeVyWS8uLlTVdAv724vFYrB/uVzWq6srbTQaenp6OrR90jGX7fDwcKgto7cvLi6Gbo++P+Vyeej+mwC7eaPrdACR8Zdw4Q3XnbyfCPC9730vWH1GpDNlXze0X2nqMV977VlybwwRJarf76NaraJWq43N+BzHGTunMpylhbuARwf2hJ2fn1/L5sJdpeFjThtVG4ejoyN4nheMzN3e3h56na1Wayg7dhwHx8fHcF03eK9SM9c0joiclUvJDkzwxycAqp/CDEx4H+8qAAWgB/iOKqANHAT73YEZmPAF7gT7AdA2/r0qoEW0g30f4qEqoA/xfrBf0e7XRnHo+b/AHVVAP8PLN/qriog207isU9Vkpn6WGs4uC4WCXl1d6dHR0VA2Fz7OpGOuMzAzXYSELp/ZbX8duv9v7M+T0H6v223PAQDvvfeeCdDFbwIAOu1BiKw/NHvWH5aCN7rT/t8AgFIxHEqB/B3ABXCA3y7npRLR2ul2u0FG519/8uQJPM8LtgPA/v5+sI9fAOLw8BCe56HVasF13aBk4aRj0nzEBObNsLOzo+12O+lmXCMiABpooA0HJ3DQwAnMwIj3cYYzfIAzAB8Ej1AoxDwWg9/fp9jDHj7Ds9u3cfc3v1nlSyAiyiQR6ajqzuw9p9vQzDRd3nvvvdiO1cG/xf/68svYjkdERLMxM10zJss158KJiGi6uDJTzjNdOwdJN4CIaOMwmK6ZAyz8BxYREUXEYLpmXPhz1FJQ+JmIaEOsVTAVEQdmztD1hfI2hGu7eRlKiYhWZ22CqQ2kVQDpqHqckCpMKScGUyKi1VmbqTGq6gJY76G6c7qDHiCCXlDK0JQ97EgJEEEpVAKxLnVABPVQ2cOSmLqLHVsCkaUOiYimS11mKiIVAHdVtTbmvkOYzHMLCAIohdy+/Qz48vXZO0bw5Zd3Yz0eEdG6Sc08UxEpAygC2AXgqWp15P4jAOeq2grdfqaqzdA+DQCn/j6jNmGeadzstFWk5GNCRBSrtauApKotVT0GMGnwkDMSJM8BzLegHt3YF8jjC+STbgYRUaqlrpt3HBEpjtl8CWDyKrIUi7wt7E9ERJNlIpjCnCMdXaCvDwAiklPVvj3XumO3XW7y9Jg4+TlpL9FWEBGlW2q6eWfIwQ46CvGDqz8YqamqJVWtTgqkvV4vNLp1cKnX68trecY9B5ibEhHNkJXMdNxy735wnbyk/Ih8Po9ejzlWNI3ZuxARbbi5g6mIvGWv+kURHABQ1Z/E3agxLmGy07Ccff5xgZZi0gim7rIMBBHRJFG6eT8CUFDV3wJ4CuBVAE9F5IdLaVmI7bYdDZpbAMZOgaH4mMXKT5JuBhFRqkXp5j1X1Z+LyJsASqp6FwBEpLCcpl3jikglNK90F+yDXDrHvsWsjkFENFmUzPTK/iwDaIa2xzKdX0SKtsJRBcCHInIYnhJjKyIVRKRi97sIF2yYR3gAEgcdzefE5qZ7cgaI4Gxo8Bbg1yX0yxOKAGeyB4hgL9jPgSMuIILzP/pO0i+JiCh2UTLTkojcAlCDPYEmIvdwfZTtjdiu3C6A4yn7TLxvHhyAFN3t289iLSf4+dffju1YRERpEamcoIgcwJT6eyoi3wdQAHClqn+5rAbGieUEk8XShESUNnGVE4wymveRqj7wb6vqzxZ9ciIionUQpZu3JiIXAJ7YEb1EkdwJ6iix1i8RrZcoA5B2bXfuXRH5fmjeaWZwAFKyengdPcS7PBwRURrMnZmq6tPwTxF5S0T+CmbJs0ycM+UApGQxLyWidTV3ZupnojaIPgHwSwCfA+jYTPXPltRGWhOv2wsR0bqJcs60KSL+OMxjAAeq+pW9/SsAEJEfZCVLJSIiikvUQvcf+d28o0TkuzG0h9YapyUR0XqKNJp3NJDauaZXqvpLACUAnTgbR+ulHRTL58eEiNbL3OdM/XmlIvKGf4H5Vty39/+FDaqpxdG8ySqhixK4ZjsRrZ8oRRvehlktxj9vKgBeAVBdQruWgqN5k1Wy3bzMS4lo3USZZ1qFWS3mVQCOqm6p6jdgRvQSzdS1uWld6oAI6kEh/BJK0gFE0JHSUNH8nuQBEeSDfV244gAicMTFa689S/plERFFCqbnquoHzluh7a/E2B5aY7dvxx/44izCT0R0U3MXureDjdSuaXoPwCv2+g9V9SdLbWVMWOh+vbBwPhEtKq5C91EyUw/A/xSRl+2o3o9E5AWA7UUbQXQT7+MM7+Ms6WYQEUUqJ/grADuh2++KyNt2O9HKneEDe42pKRElK0pmeo2q/kpEXo6rMcvGqTHr5cxeiIiSFrUC0jgnAO7HcJyl49SY9cK8lIjSYmwwFZFXAFxh9veUzLEPERHRWhvbzWsL2DdV9Rv+BcC7AL41sm0HwIcrbC9RiIJ/yxFRGkw7Z3owcvuV0DxTAICqdmEyWKKVUwgUknQziIgmB9PQ8mqz5GJqCxERUSZFGc37jr9AuM/efifeJhHNR/zcVIAz2QNEsBeUHXTgiAuIwBUnKE+Ylx4ggl6wn0Ckg46UcP5H30n6JRFRRkWZZ/pjEflrESnBFHAo2J/3ltU4omlu334WWznBErrA17Eciog20NzlBIMHiBRh1i71Ji0UnlYsJ0iTlMSsZdPRUsItIaJVSqKcIAAz6EhVT7IWSAEWbaDJ/BVtiIhuIo6iDZnBog1ERLQMGxVMiSZ5iLq9Vp+yFxHReAymRADq+Di4RkQU1UKF7onWRR0Mo0R0cwymRAA+thciopuYGExF5I0x2w5E5Nci8jsR+QcR+Y/LbBzR6hTthYgoumnnTPdF5EJVfw4AIvIjANsAagD69vqfi4j4+xBlVRFu0k0gogybFkxPAFwA8AOlp6p/Ebr/KQBXRB6F9iHKpA78OdtchYaIopsWTG8ByInIG6r6j5j8LePF3iqiFevYLl6WbSCim5g2AKlif17an1v+HSLycmi/zKwawwpINMkOOthBByJAT/KACPJBIXwXrjiACBxxg6L5e3IGiOBsqGg+4O/w2mvPkn5ZRLQi0zJTD8BHAD4UEQHMoCSbpfZF5AjACyA7J5tYAYkmibNovi/u4xFRekUudA8AIvIKEGnN01RgoXtaFbFrlt/gvxcRrVBihe6BIIjya4KIiAgzgqmI/EBEHovI/whtOxCR3wO4EpH/s/QWEmXQp9jDp9hLuhlEtCLTijZ8AnPOVAD8ZxF5JiJvAajay58D+AM7NYaIQvbwGfbwWdLNIKIVmVroPtyPbBcFPxrpW3Zt0CWiED8nPUu0FUS0KtOC6a/DN1S1KyLNMfv9c7xNIso+5qREm2XaOdNbgDlvGtr2xL8SmmsqS2gXUcYd2AsRbYJpmanfhft9EWmp6j/6U2Hs1JgrEWkAaKyioURZcoCFR9oTUYbcaJ4pAIjIPQDtLM015TxTWhlONCXKhLjmmU4dgDTlyV9W1aeLPjnRunJtF6+TcDuIaDWmTY15284xfTm07ZPQHNPfheefEtFAFS6q2am0SUQLmjYAaQemG/e3QLCe6QuYhTW+BeAugNsi8sOltzImLHRPq3QHPUAEvaFC+B10pASIoCSdoGh+XeqACOrBfiWUpAOIoCMlFs0nSrlp3bwvVPUkdLs/cvtzAF0bZDOBhe5pVW7ffgZ8+Xpsx2PRfKJ0mxZM38Hwot+XE/bjCAuiEb/5jR/8FHmM/ifphP711QHU7b9hih1OPiNKvWndvI9EpG3r874Bc570u/45VBF5Y2QOKhER0UaaGEztlJd7MOdOPQDn9nIlIr+D+cM6p6o/WUVDiTbVF8jjC+STbgYRTTGrNu9XMMXuPxKRAoA3Ybp7+6r6+QraR7Tx8niedBOIaIa555mqqgeToRLRCvk5KYfOEaXXrPVMfygiL+yc0r8SkX8Xuu9tO++Ua5oSLdFzeyGi9JqYmdopL/cB/Bima/ddmKkw31XV/6uqvxKRrwD8A4D/tJLWEm0klr8mSrtp3bx3R+oV/kxEagCeiMgnqvpLmCIORLREDfj1pFmckCitpgXTayVXVLUP4F3bvQuMTpUjotg58GulsDwhUVpNPWfqs/NMA6r6Y5j1Tj+Mv0lEFOagAYddvUSpNnUJNrue6QsAnwAoqerfj9xfBvDXqjpXUE4al2CjLOJqbkTLE9cSbFODoM1AXQDfGg2k9v4WTIZKREv0Ps4AEZwNFc0H/Er5fsF8EeBM9gAR7AX7OXDEBUTgigMRsHA+UcxmZpSq+tW0Ag1ZWhycKItu344/8LFwPlG8pnbzrht28xIBeTHlH3rKEoVEcXXzzl0BiYjWQw/+0nCb84c00bJlYuAQEcWnB5YmJIrbRgXTXq8XDN6o1+tJN4coEa/bCxHFZ+5uXrt2aXvcqN6syOfz6PX4NzkREcUrSmb67riN/mLhRJQVbXshorhECaaPMX4JNhYMJcqQNhy0+d+WKFZRRvPuAvhERDwAfbtNANwD8JO4G0ZEy1FCN+kmEK2dKMF0B8AxzHJsYbn4mkNEy1ayXbxcpYIoPlGCaU1Vn45utJkqEWVEF6Wkm0C0duY+Z6qqT0XkZRH5gYh8FwBE5C1V/dXymkdERJR+cwdTEbkHoAszqrdoN38uIn+2jIYR0XI8RB0KQT0ohF9CSTqACDpSGiqa35M8IIJ8sK8LVxxABI64LJpPZEUagKSq3wKCwApV/UrEXyCKiLLgT19qAV/HdzwWzSeKFkz/bsJ2FvgkypDdf/0bAEDdXgYUJYz+h+6F/vU5AFy4AJ7Lmd22F3cziTIlSjB9x1bX/yfY/28i8gaAdwD8PP6mEVHaneEDe41/U9NmixJMHwHoiIgC6Id6d+/F3ioiygTmpUTG3MHULgL+LRGpAHgTgKeqP1tay4go9ZiXEhmR1zNV1eYyGkJERJRVkZZgE5EficiliPxORF6IyH9dVsOIKAsUzEuJoi3BdgDgLoB9mIL32wAcERFV/csltY+IUkwhwTWiTRalm3dbVT8M3f4cQEtEfhRzm4iIiDIlSjCdVOaES1AQbSixGSnzUtp0Uc6Z3pqw/RX/ioj8YLHmEBERZU+UzPRdEalheIHwAgBPRKowa5uWAPD8KRERbZQomWkBwI8BuKFL+HYDXCKRaKN8ij0oBHtBIXwHjriACFxxgoL5eekBIugF+wlEOuhICRBBSTrBviycT1m08HqmYVzblGiz/PFLV7EWzX+IOvAlYCYOEGWHqK7H0AERcWC6oAsAWqp6LbDv7Oxou91eeduIaE5+mdI1+V6i9LM153cWPU7kCkhpJCIFmKk7rr19CjMflogypD7ykygrIlVASrEKgIvQ7eKkHYkovT62F6KsSV1magvp31XV2pj7DmG6crcAwM9EAbyK4VHGEJGcqvaX3FwiihX/DqZsSk0wFZEyzP+kXYwERnv/EYBzVW35t0WkEiq8v7WyxhLRUhThzt6JKIWi1OZ9y171A50DAKr6kzgaYoNkS0ReBZAbs4szkq2eA6gBaAJ4MbLvFrNSouzpwB8HwgFIlC1Rzpl+BKCgqr8F8BSma/WpiPxwKS0LEZFxfT+XAMr2ehOmYAREJAegtew2EVH8Oiiiw65eyqAo3bznqvpzEXkTQElV7wLBSNpl24IJnmF9+/w5VfVEpBPqKr52vpWI0m/H1n1hXkpZEyWYXtmfZZhM0LeKz30O18+J+sF1C0A/NBiJWSkREa1UlG7ekoh8Hybr+ykAiMg9rGbgz7jzn/7zjmasE/V6vVAps8GlXq/H0kgiItpMc2emqvoXdoHwqqr+0gbWAgYZ6zJd4vqgpJxt19wDjfL5PHq9XpztIqIYfYG8vcb/p5QtkabGqOpJ6Kanqj+LuT2TnrcrIqNBcwvs0iVaK3k8T7oJRDcydzfvmLVK+yJyT0T+LOY2TeLagg6+XZiVaohoTeTtRcSFKw4gAkfcYEWZPTkDRHA2dKoG8Hfw9xMBzmQP53/0nYRfEW2KKJnpUDerqn4O4PO4FgS301/KMKUBt0TkAqZgfdc+X01EDm1ALQC4CBVsIKI18PZ77+EXv/hFLMfaw2exrmhDNM3UVWPsOdISgFswU066I7sUALRV9c+X1sIY5fN5ff7cdCM9fPiQA4+I1tieXYHmjCvQ0BRxrRoz1xJsIvJTmHOUj0fu8lT1V4s2YlW4BBvR5hAbTNdlmUlajlUvwVYDUF7VgCMiosUdJN0A2iBzBVNV/QrA2EAqIm+o6j/G2SgiokUdYOFkg2hukVeNEZGXRzbVAGTinCkRbQ4XVXvNSbQdtBmiTI05EJHfwxRp6Id+ZuaTGq6AxMFHROvNxQFcdvXSikTJTLcB3LJdvgER+STeJi0PKyARbY6qXRs1M3/tU6ZFqc17PhpIrUdxNYaIKE530ANE0Bsq8tBBR0qACErSCYo81KUOiKAe7FdCSTqACDpSggjw2mvPkn5JlFJRgqmKyBtjtrMfhYhS5/bt+APfl1/ejf2YtB7mmmcKACLyBMDb9qZfJ1cAvK2q31hC22LHeaZEdFN22io4bXW9rHqeKWCqHf0Yw8uhCYDDRRtBRJR2XNGGpokSTGuq+nR0o4i8iLE9S+WP5gVYTpCIouGKNjTN3OdMVfWpiLwsIj8Qke8CgIi8laVygvl8HqoKVWUgJaJI/BVtKJrXXns29wAv/9KTPCCCfLDvYBWhtK4EFGWe6T2YQvfvwhS9B8yqMatago2IKDHP7YWiiXPQVg93cPr1f4nteHGK0s27q6rfAoLAClX9Svx+UyKitcblk2+iiA4AoK51AHWYf8MUJQDD47p6oX99TjAIzI29lYuLEkz/bsJ2jm0jorXXgD8TgGUgougENZLXO1RECabv2CHE/wT7rth5p+8A+Hn8TSMiSg8HJ/ZaGvOi9OrYs4KlGI71Ps7stb0YjhavKMH0EYCOiCiAfqh3917srSIiShnHdvMylEazY9+xOPLSM3xgr6Uvy41StOGRqj4QkQqAN2EWBs/U+qb5fF6fPzdDCDg1hoiiYNGGm4lzkfYze6y9GH8JSRRtqInIBYAnqvrbRZ84CSx0T0SUXenNS6MF031V/ZmI3BORAoALVf3lshpGRJQmaT5fl2ZfJN2AFYlStOFn9udTVT0B4InIX4nID5fWOiKilDjDB6Fzduvr2Wuv4WxolR0E1RTChRXOZA8QwV6wnwNHXEAErjjBfp/hIMZiF4p05qXRija85f8UkZ8C8AB8BVPIgYhorZ3Zy7q7++WXsebeddRx/tKfxnIshUCRztIGUQYg/RqDPwkaAE4mrG+aWlw1hohuajCDQYMvdAllSZ9iD3v4DHsAPgMAHOAAO3BRhYuDYLHyO+ihh9dx/tKfYvdf/2aVL2Eue2L+ZDjTFHZnL2EUWBIDkADgo3HF7omI1t17772HX/ziF7Ed70++/tvYjhWnz1J8Ttj/4yWNHb1RMtPvZ20qzChmpkSUBnFOF4lbmqcALaNtcWWmkQcgZZm/BJuIcI4pEdEYB3BxwNIUkU3s5hWRHwDIwSwG3lbVv7fbXwFQhlksfBfAr1X1v62grQvjPFMiSof09pC5qNpr6atB/GnQBZ2+oWATu3lF5PcAin4QnbBPDsALVf3GktoXK3bzElEadMRUqi1pJ+GWXOeKCaKOpjA7zegAJHdaIAUAVe2LyMm0fYiIaFgpxTMKq3aVl/TlpUhxXjr9nOmv/SsiciAiPxWR340p0vBrEBHR3Epoo5Tart6qvaTPZ/CnHaXPtMw0mBlrKx6diAhU9SeT9iMiotm6sSxIthx3km7AVAdJN2CiacF0XKf0uCw0hQOoiYjoJtI8RPMAC5/aXJppwfR+aM1S3+6YbfcBjGarREQ0wUPU7bX6lL2S0bO5aXz1dOOT5pHG04JpCcCtMdu3R26/GV9ziIjWXx0fB9fS5nWbm6axy9G13bzpC6XTg+mxqv541gFE5JMY20NEtPbqIz9pPn594zQG02mjeRtzHmPe/RLHCkhElAYf2wutj4mZqap+Ps8B5t0vDVgBiYjSoZh0AyZqByON01dQ4k4wPCp9Z3Tnrs1LRETxKMLFp8ijI6WhBbd7kgdEkA8W3HZNRSIROOIG++3JGSAytIj39773vVjaVkI3tUUlengdPbyedDPGiroEGxERLeiTl/47dr/+W3Riy1AVca0O5xeTSF9eihTnpRGWYFsHrM1LRGspxpq1YusGawrrBi9j6bqkFgcnIqK1ls4u3rRjMCUiyjixs0LjyNcexnCM5UlvzyKDKRERBepJN2CKdjDDNH1d0AymREQUqNvctJ5sM8ZK6yhjgFNjiIgy71PsQSHYC6bKOHDEBUTgihNMqclLDxBBLzSlRqRjFisXQUk6+Bh1fJzKUJrupeuYmRIRZdwfv3QFfB3f8W7ffgbgbnwHjEmal67j1BgiIsqEGGcAhY7JqTFERLRB0rzQJYKJAAAQR0lEQVR03UadM2WheyKi7KrbM7pptFGZKQvdExFlV33kZ5psVGZKRETZleal6zYqMyUioixL79J1DKZERJQJRbhJN2EiBlMiIsqEDvwZLOmb0slgSkREmeCv/5rG0g0MpkRElAk7tsB9+vJSjuYlIiJaGINpHJpNYH8fuHXL1Lsqlcztfj/+5zo+BnZ3zWV/H6hWAdcFPM88b9R2+21uNuNrY6s1/v3Y3TXXS6X4nq/bBba3zfMcHw+2u655/lYrnueJ+3hEtFbYzRuHSsVcdnfNl221CjjO7MdF4QeonR3g9BTI5Ybv296OfsxKxfzc34+njb5y2Vz89+PBg8FzASbo7e8Dh4fA0dFiz1Usmvdj9A+Jiwvzx0xcf9DEfTwiiuwL5O219BXfYWaaBf3+IOCdnw8HUsAErkZj9e26Kf8PjXAmGbejI+DqajiIp+l4RBRZHs+Rx/OkmzEWg2kW1GomoD54MHkfx7keZNMq3M7uEhf7jfv9yMr7S7Sm8vaSRgymq1KrmezS/xnuWu33Tdewfw50d3f4nKJ/nq5QmP4cn39ufobPI9Zqg2P45zDHZYTn52bfWs08v/+4edsYRbirtFgcPnfrd5Nvb5vzlDd5bs+7fh41/PqbTXPs4+PB+VzA3B7dNs/xXHf88XzNptl+fDxof/jc67TX77qD5y6VBu+daxZ+xu5u9PefKKOe20sqqerGXEqlki5VuawKqDYa17fncoPbh4dmv6src7tQUC0WB/d3Oub+iwtz2yzfp3p+Pn9bzs/NYw4Pr287OhpsOz012xxn+PG5nGm3b1Ybx/Hfj9PT4e2Oc/19ajTMtkpl8Pz+9VnP7d8Ov65x2/zXXy4P3nv/ecdtC7c76vHCjy0WzTZfpTL8+5/1+v3n9m+rmseO/s6I1hzQUKAxe8dIx0RbY4gvzEyXrds12caHHw62HR2ZAS25nMlKPA+4f39wf9HWn/TPg/oZaZTBL1tb0do5OoCnXDbt9rz52jhNo2GyLT+7BMzrDw/S8tt79675eXVlBhYt+txh/nMUCoMuW/+9HbfN825+vPBjT05M5u/zs8lwdjrp9QPm9ZbLg/cCAB49ut57QLTmGmijgXbSzRiLo3mXzf/yGz3fNvql++zZcPdrpQK8+qq5Xi6bbr1nz6YPgNnfH3wBL8r/cve8+do4zf7+/KObR7uyF33uccad+7zJaOhpxwvzu7L9c9/TgvSkrvxazQTfoyNz8bzZ3f5Ea8bBib2Wvhq9DKbLNivL8bOsu3fNVJFxjo6AJ09MQJ00lcTzZmejUTLby0vzM/yFPa2NyzLP+5Nm/b7J+nM54OlT89N1o89XLZfNe+GfR542GI1oTTkwvVHpC6UcgLR8fhddqzUczI6PTRew/yU52mXZag0ysVxukHHu7l4PzJ5nulDDgXZcthTuahx1cTG43u+b569UTDCdp43LkuRzx8HvKr9/f/A7Cb/XUfgBtN0e/JFBtEFO4OAEMc/hj8lGZaa9Xg8iAgB4+PAh6vV6PAduNoHHj82XHGC++M/PzbmyXG4wUvbePVN0IZczmZb/hdjpDEbR+png9vZwJlYum/No/ohQwBxna8vsOxooCwUTXBsNExxzucF50UePzE//+H6BBf8cXLdrvrjDzz9PG32tlnneSe/HuP39Nj16ZILPvM/d7Q7a/fixub9QGN5WLJr3yd/WbJrH7+wM/gBpNAbnPv1tcRyvUjFd3I8fAy9emK7p+/cH3b6Xl2a/aa/f5/9xs2ihCyKKnZjBTJthZ2dH2+10nrwmmih8fnR3d3oPA9Ea25MzAMCZ7sV2TBHpqOrO7D2nYzcvUZrVaoPBUf48YKINdYYPcIYPkm7GWBvVzUuUOffvD7qyt7dNlzzRhjqzP+PLS+PDYEqUZsUiu3WJLD8nTePJSXbzEhERLYiZKRERZUQac1KDwZSIiDJBIcG1tGE3LxER0YIYTImIKBMEai4CnMkeIII9EYgIRBw4YpYmdMWBiFmlMC89QAS9YD+BSAcdKZkdYsJgSkREmXD79rNYj3eG92M7FisgERHRRjKJKSsgERERpQKDKRERbaSHqMd2LAZTIiLaSHV8HNuxGEyJiGgj1WM8FoMpERFtpPjyUgZTIiLaWMXYjsRgSkREG6kIN7ZjMZgSEdFG6mDh6aUBBlMiItpIHXbzEhERLWYHndiOxWBKRES0oLUKpiLiiEh8eTsREdEc1iaYiogDoApgK+m2EBFR+n2BfGzHWptgqqougMwuCVOv15NuAtGN8fNLWZTHc/MjBitdgk1EKgDuqmptzH2HADzYzNIGx6jHbwA4VdXWuPvTvASbiGCTlsOj9cLPL2VRXgTPAajqwquE/2EM7ZlJRMowpSZ2YQLm6P1HAM79ICgiRyJSUdXmKtpHRESb53mMx1pJN6+qtlT1GEB3wi7OSDZ5DnP+E0AwsOhwzKUcpR29Xi9Su+ftuppnv03sBkvyNS/rueM47k2PEeVx/Owuhp/deI8R9XGr+/w25nqeeay6m/cIQE5Vw4GyCOCpqt4a2daJmnrP6uYVEY3yeuftuppnv1n7rGM3WZKvaVnPHcdxb3qMKI9b5Wc3atuygJ/deI8R9XGr+vy64qCKk1i6edMwAGkLwOXItj4AiEhu3oPY87E7APY5PYaIiGZxcBLbsVZyznSGHK5PZ/GD6xZsYJ3Fnl+ddY71/4nIvxmz/TmAcX3AeRGZp294nv1m7TPvc2VJkq9pWc8dx3Fveowoj1vlZzdq27KAn914jxH1cav8/P6HCO2aKA3BdFyw9IPraMa6EFV9Kc7jLZvttq4BKACAqk4650yUOrZn6UPYQYeTTr8QpY397F5hMGB2C8C9ad/BaQimlzDZaVgOAFR1rqx0je0A+ByAO246EVHKnarqrv1iOgHAYEpZsRM+jyoizqxkJvFgqqpdERkNmlvgfzwAeMTpQZRFdqS9n5H2Aewn2yKi+YV7Uew0zZl1DxIPppY7Mq90F3GOWU7YAsUqtkJzdJuqem2OLtGy3fDzW7T3FwGUAbR4moKSEEOxoLlK1K6qaIP/H6oCEyAuEPrPpao1O2+0AnN+8GIdMrJFi1X4v1gR8QCcAiitqu1EMRRb2bI9Tx6ADoDtFTWdKJZiQbbm+1y9pCudZ7qpxs2vtduvRubXlgHU7HmmCoCCLXZxbV+iVbnh59cBsO1nA3aO98Jz+YiiusnnN7Sto6pzJTFpmGe6kSbMhb2EyeABM8q5afctIMNF/Gn9zPH5fQI7Ct0OQOIpCkqNOT6/vsK8x0zLOdNNNLVYhaq2RKRif+l3wQEclC6zPr99EXlsM9Qc+PmldJn5+bXb5v4jkME0OTOLVYT67jN//pjWTpTPL1HazFUsaN4uXoDdvElaWbEKoiXg55eyLPbPL4NpclisgrKMn1/Kstg/vwymCbHTglisgjKJn1/KsmV8fhlMk+XaKTC+tSpWQWuPn1/Kslg/v5xnukShYhVVmL96HmGkEkyoAkcBZtDGzLJVRKvAzy9l2ao/vwymREREC2I3LxER0YIYTImIiBbEYEpERLQgBlMiIqIFMZgSEREtiMGUiIhoQQymREREC2IwJcoQuzboRhKRgoic2sWek2yHIyKNkeo5tOG4BBtlhl0kvQrgEEAXwOPQ3dsATlU1FbVhRaQBAKpajel4OQAnACoAZMp+s94jB0A1w5WKnqnqcXiDfc01ABehzf7noDBrKTgbFB8AKAJwAdTCxc7tmqwNe8yaqrr2ORlMKcBgSpmhqh6Amv3yezzmS/VURAopCRSnGFnKSUScm7bNfrnvi8jUkmVzvEcNmBJra8G+ziqA/ZEAWIYJgDOzWFVtikgLwBXMH2T9kftdEdlW1Vq8rad1wmBKa0NV90VERaRlg0qSbRmXIc+90PCyqGpXRHaSbkccbHZ4AuDNMQGwJSLd8Y+8TlX7NqBWMbJyiO0VeBZDk2mN8ZwprZsW5shGlklEciJStoW2/dsNDBYfTqJNTuhmO6l2xKwBU7h80vqTUVcAOQJQGXNe+sNZXcVEzExp3ZzDZBcAgqziAUxmcRfAuc1ayjBfnp59zKW9/yLcFWsf79j9AHMO7jh0v39uMofBEk592wZ/WxkmkBbt/sHqFJPaFzr+EQZZ0SKLbgdZsb9qxsh74Aee3QnvwbT3sA3zHt639/nnFKuhx4TPZ15i8AfPvs2Wy/YYxxG6U/0VQcaybSzMeh0j+/dhft/h7vGNHfRFEagqL7xk6gLzxXw44T7HfKyH9i2M3M7Z6xWYwBF+/Gn42AA6/v72dhlAI/RcldB9FQDF0H7nI/edTngtk9p37h/P3i6EX9sc79E5zECkxqTH+e+B3waYwKER2uj4t2EG8PivvzPpvQhtG33vnRmvqeD/bvx2ht//Od+Tsa8jtO0o3C7bzuK0tvDCi6qym5fWTg42g/OnLujw+dMugA9Dt0fPrTZgshc/c4OGuhHVZDKOzXI8AEd2qoQ/anTuc7XT2me7iHc0tPaiRj8PfK6qx2pGFAfZ3kg3Zj98bP+1+vvM8R5eAvBUta+q3VB7ixi8Fx6AofO09n1EqCu8DODJvC8s9Du51nVufx9HInJuB6UV5/wsAOb3X/DbBRNI5z73SpuL3by0brYxOCdYAND3g6L1DNMDnodBt14RIyNyrT5MoGuJSBWmq7FhB7zci9DWae3bmdHOSHR4VG8ZQPgc4LTu43new3GPb4aeZ1KgPIL5w2UfJkOM2o3dgumWHhohrYMudIWZBuR3I8/8LKiq5w9EEpHahNdGdA2DKa2bDzEIaH62FR6dOWseagGDL1AP4wcN5QB4IlKx2aifZR3h+vm2a+x5PD+zHds++6VfuP7oWEQ57k3eQ8DMb92yGeGljp9v+wQmsw9nsVFUAXREZJ5AHOV1NGC6+zuIkC3TZmM3L60NETkF8MjvlrOBzg9e/j65UBcecD2wBF2i9vG5kcdXADRtd2FhJNN5jMm80HMVbJfoxPbZL30v3NaRdt+I7b7dnXf/Od/Dce6qqquqTZ0wEtYGwCcAjm7SlWp/BwcAnobbZ9tYRiirjPI6Qu0t3SBbpg3FzJQyIzRCtADgvkhQCOhVmEyvodfnd94D8EBEgnmCI1/ung2QfZhu3XMdLqxQCj1+C6Y7ct/e14cJqH4lnIKqHtsv6BqAHb9Qg+1qbNspKuGu42nt8+/zs2X/POYpgINxX/RT3iNgUAHJ7wb12+mPMnZhzxfDZIxHNmCNbaN9fNW+zkMAbqhNFyJyFXqtHszvZzSwNjByPjUK244uTKEK/7lf2OcrYfiPpVmfhbBjTP/jiGiIqE4tqEK0tvzqOao6d6ZGs9kgex+ml6Bvs2G/wMJBOAsNdZXPc9wCzOjdqd3oq5CmtlA6sJuXiOJWhqmh648U7tsA2gbM3NyR7nGizGM3L22kkS7KG9fMpetsV7cT6tLegumiPrfd3Vsw3eN832ltsJuXiDLBdq0eYfLo4FW1w4E5H3s+bxc1rT8GUyIiogXxnCkREdGCGEyJiIgWxGBKRES0IAZTIiKiBTGYEhERLej/Aw4+6xDtGmozAAAAAElFTkSuQmCC\n", + "text/plain": [ + "<matplotlib.figure.Figure at 0x7f1179863390>" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fontsize = 16\n", + "\n", + "fig = plt.figure(figsize=(7, 5))\n", + "ax = fig.add_subplot(111)\n", + "\n", + "ax.set_xlim(binning[0], binning[-1])\n", + "ax.set_ylim(1E-1, 40)\n", + "\n", + "ax.set_xscale('log')\n", + "ax.set_yscale('log')\n", + "\n", + "ax.set_xlabel(r'Deposited Energy [GeV]', fontsize=fontsize)\n", + "ax.set_ylabel(r'Events per 2635 days', fontsize=fontsize)\n", + "\n", + "ax.tick_params(axis='x', labelsize=fontsize, which='both', direction='in')\n", + "ax.tick_params(axis='y', labelsize=fontsize, which='both', direction='in')\n", + "\n", + "ax.step(\n", + " binning, null_exp, drawstyle='steps-pre',\n", + " linewidth=2, linestyle='-', color='black'\n", + ")\n", + "\n", + "ax.step(\n", + " binning, muon_p, drawstyle='steps-pre', label=r'$N_{\\rm{muon}}=1.5$',\n", + " linewidth=2, linestyle='--', color='blue'\n", + ")\n", + "ax.step(\n", + " binning, muon_n, drawstyle='steps-pre', label=r'$N_{\\rm{muon}}=0.5$',\n", + " linewidth=2, linestyle=':', color='red'\n", + ")\n", + "\n", + "ax.legend(loc='upper right', prop={'size': fontsize})#, bbox_to_anchor=(0.96, 0.96), markerfirst=False)\n", + "\n", + "ax.text(\n", + " 8E4, 0.2, r'${\\rm\\bf IceCube\\:Preliminary}$', fontsize=fontsize, color='r'\n", + ")\n", + "\n", + "matplotlib2tikz.save(\"plots/tex/syst_muon.tex\")" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAFDCAYAAABhmOTuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3c9zIumZJ/DvM7O7hzl0U+qYQ1GH6Ua+7Fy2Dar+A1zgg1RyxNqoavdugWeva4uunQiXPBG71ah93wF576sS7Yjp6nLEDFT/AVOA574tevZgeg5TKtqXOWyMnz28+UKCkh8pXiCT/H4iMiSSJHlJIR7eX88rqgoiIiK6vT/adAGIiIjijsGUiIhoSQymRERES2IwJSIiWhKDKRER0ZL+zaYLsE5/8id/ov/yL/8CALh79y7S6fSGS0RERJvU6XT+WVX/dNnzJCqY/vmf/zna7fami0FERBEhIv/XxXnYzEtERLQkBlMiIqIlMZgSEREticGUiIhoSQymRERES2IwJSIiWhKDKRER0ZIYTImIiJbEYEpERLQkBlMiIqIlMZgSEREtKVHBtN/vQ0QgIjg9Pd10cYiIQmm1WigUChARtFqtsfsajQbu3LmDXC6HRqOxoRK60Wg0UKlUAu87OzvD2dkZjo6OcHZ2tuaSTZeoRPfpdBr9fn/TxSAiupV8Po9er4dUKoVqtYp8Pj+8r1gsotfr4eTkZIMlXE6r1UK320Wz2UQmk7lxf6VSQbVaHd7O5XIAEInXnKiaKRHRNiiXy2i1WhgMBmP7U6nUhkrkRj6fx8nJCbLZ7I37BoMB3nvvvbF95XIZz549W1fxZmIwJSKKmXw+j0wmg3q9Ptw3GAyws7OzwVKt1vX1NSqVCnq93nDfzs7OjS8Um8JgSkQUE91uF3t7ewBu1sra7fZYs++2yWQy6HQ6Y82/zWYzMq85UX2mRERx1m63USqVAAClUgmVSgWtVmvYl7rpwFIulxc6rlqt3qpJ2t/8OxgM8Pz5c3Q6ndDnWQUGUyJKDJHp99VqgBenUK8Ds+KC6uj3XA7odoOPOz425wKATscc60oqlUKxWLwxECmMer0+DM4u1Go1Z+ea5+joCK9evQocqLQJbOYlIoqpJ0+eoNVqodfrhe4vHQwGuLy8XFHJVqtSqaBSqQQOVNoU1kyJKDH8NcpZSqVRLXWeRVsZl62V9nq9YX+plc1mkclkhrU0q16vI5PJ4PLyclhb7Ha7uL6+xs7ODi4uLlAoFNDr9dBoNJDNZtHr9YZTT6rVKprNJs7OzpDNZjEYDJDJZOYGr1U38wJmDmqhUBjWxrvdbiSCKoMpEVEMNBqNwPmU5XIZtVptGJx6vR6azSYuLy/R7XbRaDRQLBZxcXGB+/fvD4NQNpsdNhUDZoDPYDDA3t4eLi8vhwHZHl8oFNBsNmeWcdXNvK1WC9fX18jn8xgMBri+vsbFxUUkgimbeYmIImwwGKBcLqNSqQTW/Eql0jAgAiYolstlNBoNvHnzBtfX1wBMk/Dr16+Ry+VuZE+yUqnUcAtKnOCflrIK3W4XZ2dnaDQaeP78Oc7OztD1OqQHgwEKhQLK5TLu3LmDO3fuYHd3d+VlWhRrpkREEZZKpVCr1abW+mw2JKvRaAwzIdlA2uv10Gq1hsfZoGz7We2IYH+/6/3799Hr9Ya1vuvr65UP9slms8hms4E18FQqBV20nX4DWDMlItoimUwGV1dXaLVa2NnZQafTwWAwGO5rtVrDNHyFQgGNRgM7Ozvodrtot9vDvL4nJyfDIFyv18cCNt0kUY70ru3t7Wm73d50MYiIKCJEpKOqe/OPnI01UyIioiUxmBIRES2JwZSIiGhJDKZERERLYjBdQi5ncn0Gbf7sKZ3O9ONEFs+gQkRE0ZSoYDorqPmWBUS9Pjv4ERER+SUqmLrW6Zhcn0GbPzjnctOPUzX3l0rTg/dkTs9Zgd7/vEREtB6JCqazgpq/WbZUmh38iIiI/BIVTKOsXp8evCf7VGcFeodLExJRxLRaLRQKBYjIjfy6jUYDd+7cQS6XG2YxirtFV6GJAubm3TK2mZdBlWj75PN59Hq9YT5e/6LgxWJxmJN3G1QqFcQpYx1rplumXDYbEW2vcrmMVquFwWAwtv+2a4RGTVRWggmDwZSIKGby+TwymQzqvhGHg8FgbNWXOLPN2XHCYEpEFBPdbhd7eyYne7lcxrNnz4b3tdvtsWbfuGq1Wnj06NGmixEa+0y3lJ0P6x99nMsB3jq7Nxwfj/pbOx1gz7eGwv4+8PLlaspJRItrt9soeQMiSqUSKpXKcC3SXq+38WC66ICharU6tUl6MBjEsrmawXTL7O8Dv/mN23O6Ph8RLS+VSqFYLN4YiBRGvV4fBmcXpi1gvqhGo4FiseioNOvFZt4t8/Ll9Dmxt0kyQbRVgtKYHR6afS9ejPbZNGj+QNPvm33p9PjjbV5R/xy201Oz7/R0tG8FeUOfPHmCVquFXq8Xur90MBjg8vLSeZluy45SjivWTImIYqDX6w37S61sNotMJoOjoyO8evVquL9eryOTyeDy8nJYW+x2u7i+vsbOzg4uLi5QKBTQ6/XQaDSQzWbR6/VQqVRQrVZRrVbRbDZxdnaGbDaLwWCATCaDbDY7s4zLNPN2u130ej10vb6o169fYzAY4OzsDMViEZlMZqFzb4yqJmbL5XJK4dg6KhFtVrVanbo/k8kMb19dXWmxWBzed3l5qaqqJycnw987nY6qqmaz2bFzZTIZffv2rb59+1ZrtdrweFXVfD7v7sUsoFar3SjfKgBoq4P4wmZemqndNhsRbcZgMEC5XEalUgms+ZVKpbF+xkwmg3K5jEajgTdv3uD6+hqAaRJ+/fo1crncjexJViqVGm7NZvNGbXBd8z/r9TouLy/R6/VwdnZ2Yz5tFDGY3sLBwQFExOl2cHCw6ZcVKJe7mWifiNYnlUqhVqtBVQMH+NhsSFaj0UC320WxWMTu7i4AEwSfP3+OarWKTqeDq6srABj2s9rg6u93vX///ljwvL6+XltTa6lUQrPZxNu3b3FychKLvtREBdNOp+Mk8P1mBcNbV3FOIkqeTCaDq6srtFot7OzsoNPpYDAYDPe1Wi3kvG/IhUIBjUYDOzs76Ha7aLfbw7y+Jycn6PV6aLVaqNfrYwGbbhJN0JBNEXH2Yvf39/HS0eRL8UYXRvFvYQczcmk3ItpGItJR1b35R845TxQ/wFdlb29Po5g42QZTV9wGevMzQW8TIkoQV8E0Uc28UbW/v+/0fGwyJiJaLwbTCHj58qWz6T+rYue6Tw5GsvuDNn/TsJ0DH9FxVkRES2EwpZkcV5qZmpCIthL7TLdMlAczsf+ViKKGfaZEREQRsVXBVERKIjI7eSQREZFjWxNMRaQEoAxgO5aaJyKi2NiaYKqqdQDb3SEaQhJSHRIRRUXkgqmIFEUkMG+ViJx495e8mihNcD1nFXA3b5VrpBLRtorMeqYikgeQBVAAcGNpAi/ANlW1ZW+LSFFVG+stabS5ynxkuc7ORES0jSJTM1XVlqqeAehOOaRkA6mnCdNHSkREtFGRCaazTBmhew0gv+6y0O1xOTci2laRaeadYwcmePoNAEBEUqo6EJEigD1v37WqTqvh0oZ0+Rchoi0Vi5opgBRuTnmxwXUHAFS1oao5VS1PC6T9fj9wtOrp6enqSk43+PP3lnzDyDqd2bl+O53RsaUSc/0SUXTEpWY6CNhng+tkjXWqdDqNfr/vpkQU2v6++9y8zPVLRFGwcDAVkQ+9X+1I2xIAqOovXRcqwDVM7dQv5T1/UKClCJo30DiXW3zqTL0OnJ8vXyYiIhfCNPP+BEBGVX8P4BWA9wC8EpGfrqRkPl6z7WTQ3AHQCjiciIhorcI08zZV9dci8gGAnKreBwARyaymaDfUJ+aVFgDU1vTcREREU4Wpmb71fuYB+BMlOMlpIyJZETkBUATwyMt2NJwSo6oVABkvA9IJgKuwCRv8A5A46CgcpiYkIppu4fVMReRnMP2lVZgECl+KyAMAH6jqr1ZYRmf2RLQNjHfMHR4CX3wBfP65+R0wHXLlMnB8bH4HgH4fuHcPuHvX/J4QBwcHztIJWq7WWrVzVv2jfImIwlj7eqaq+ilMP2XZC6Q/gkn/l1y53M05G6enZp+/5mvnfExmLEinzf4IB+eXL19CVZ1srnU6DKREFA1hRvM+U9Un9raqfraaIq1QLge0JxaWefHi5nGl0vgESMAEPmZpJyKiAGGaef8AMx3muTeiN3b29va0PRlMaW1s0vxV1FKJiG5j7c28AApe3+h9EfmRb95pbER2AFJQTZjmspmRiIg2beGa6Y0HmmBaBXAZmwFIUa2Z2oiw5TU21zXThFw2IlqhtddMbU1URD4UkecAvgTwNYCOV1P94bKFSaxazWyA6cMVGY0stoKqYYeHZp+/37deT1zCWn/+Xjv4GjC/z8r165fLJe6yEZFDYZI2NETE1gHOAByr6rfe7d8CgIj8OC611Ehx2cR7egp8842780WYy1y/XNGGiJYRZgDSVzDTYl5Nuf97MOkGIxtMI9vM65INzP4qWkREeQASm4yJkslVM2+YmmllMpB6c03fquqXAHIAOOtv0yIYRImItl2YpA2fAYCIvG83mOB55N3/qRdUIyuyo3mJiCjWwiRt+C7MajG2IUwAvAugvIJyrUQi1jO1ry+d3mw5iIgSJMw80zLMajHvweTm3VHVP4YZ0UtRce+e2SLMVdJ8Js4noqgIuwSbDZx3fPvfdVgeWtbdu5suwVT7+/vOk+a7Ot/xsZPTEFFChQmmEJEfquqvAXzt+31d65nSIiLcjP3y5Uun55PJyaJL4LgtIlpGmGbeHoD/JiLveKN6fyIibwDsrqZotJQErGhDRBQVC9dMVfW3APZ8t78vIt/19hPFmv3OMfmdgohoEbfOzTs8gampxmIVmXQ6rd942YGePn3K6TEx5zIJBJM2ECXTJlaNmebcwTnWIp1ODxeqZiBdUMJWtJnM3+tvJS+Vpuf5nazRijDXL1GSBAZTEXlXRP4gIv86Z/sDgOKay0zrdH5uti23v+/+nI4HLhNRhAX2marqtyLSUNVHdp+IPADQ802PgYhkAXyw+mLSxtjVbLbcIgON63WO+iWiYLMGIE3OvHvXH0gBQFW7IpJyXyyKjAQ18RIR3dbUPlPf8mrzMJgSEVGihRmA9JFdINzybn/ktkgUKS9ejC8+HkFMTUhEmxZqaoyI/B3MUms9mMxHPQAP4jI1JhHrmboW4TkjBwcHztMTulprNcKXjYh8NrGeqU3UkIUXUKctFE5b5OHDTZdgKpfpCV2mJgQSM26LiDxLJ22IEyZtoGlcJoAgoviIUtKG2GDShiXY7AR+h4dmn79PtV43+/yjgPt9s49rrBLRlkpUMCVaF85JJUqWRDXzcgASTeO6mZcDkIjigc28REREEcFgSuuTyyVujTN/Mnw/u9xs0ObvbrbLzXIKLFG0TQ2mIvJ+wL5jEfnKS3L/f0TkP66ycLRlul2zJYDrxPlMmk8UbbPmmR6JyJWq/hoARORnAHYBVAAMvN//QkTEHkM0U4L6q+dNgfUv7TZLwiryRLE1K5ieA7gCYANlT1U/9d3/CkBdRJ75jiGajpGBiLbUrD7TOwBSvubeaeMSey4LREREFDezgqld9Pva+7lj7xCRd3zHxWbVmH6/P0xqzqQNG3B6arYIc5U0n4nziZJl6jxTEfkRTABVAHYsYlNV/1FE/gCgCuANgPMQy7VtFOeZbliEJ1+uImk+4GbeaoQvG1HsrTzRvap+NuNxd7xjYhFEKSKePt10CaZymTQfcJs4n9//iKIv1Koxlqp+O9HUSzRfxJt4o4rjtoiib2bSBhH5sYhciMh/9e079pp534rI/155CYmIiCJuVtKGTwD8BKa/9D+LyGsR+RBA2dv+AsAfeVNjiObrdMxKM5NVrXTadAz2+6N9pZLZ588W/+KF2Xd4uJ7yRkSpNJ4ViYiiZ2Yzr79T1lsUvDrRUVv3gi7RfD//uUnlk81uuiSxcn5ufnIVGqLomhVMv/LfUNWuiDQCjvtnt0WirTVtkI+/RmoFrWF2eDg+pJXDXIkoIuYlbYCI/Ni377n9xTcAyd2wRSKayibCn2wln5Ywf7KV3K7bzsT5RO7NCqa2CbdisyDZqTAi8i6AgYj8TwDNVReSKJBqImqlrpPmA0ycT+TarRcHF5EHANpxmmvKpA20Lq4XG3eJreNEIytP2jDnyd9R1VfLPjkREdE2mDU15rveHNN3fPs+8c0x/Vf//FOitTs8TNw0GSKKpll9pnswzbi/B4brmb4BkAPwHQD3AfypiPx05aV0hInut8wXX5gtwpg0nygZZjXzvlHVc9/twcTtrwF0vSAbC+l0Gv2gaRgUT59/vukSTLW/v+80cb7Lc7GvlMi9WcH0I4wv+n095Tj+a9JmRLiJ12XifJdJ84loNWY18z4TkbaXn/d9mH7S79k+VBF5f2IOKhERUSJNDabelJcHMH2nPZj5pE14g48AdACkVPWX6ygo0Q31usnr609c2++buR/p9PixuZzZ3+mM9p2emn0J6z/P5bgSDZFr83LzfguT7P4nIpIB8AFMc+9AVb9eQ/mIpvubvwG++WbTpYidbnfTJSDaPrdO2hBHTNpAU9kaawSrbK4TQDBpA9GIq6QN89Yz/amIvPHmlP6tiPyZ777vevNOuaYpxd/entmIiG5hVtKGnwH4TwA+BvAIwD/CTIX5DwCgqr8FUAdwtPpiEq1YNpu4peFExrubO53ZSfP93c12uVkmzScyZtVM76vqnqqeq+pnqloGsAvgUxH5nnfMm9UXkWgNOp3xaLHFXCfOZ9J8ohl9piLyM1X9dMp9nwD4O5gRvdeq+serK6I77DOlOGLSfKLVWUufqe/J3vffVtWPYdY7fbRsAYiIiOJu1jzTT70BRj8DcCUiH07c/xlMSkGmZ6H4S6dvzk0lIlrQvHmmH3sLgTeC5pWqaktE7qysdETrEoP5qi7TCu7v7ztJeZiwMVtEU81t5lXVb2claIjT4uBEU/3ud2aLoH3XI4bgLnF+gsZtEc10q8XBibZOhJt4XSbNB5g4n2gVFhqARERERNMxmBIBoywE9fpo34sXZt/kUm82W4Hf4aHZ9+LF6ssaIUGXgiiJEhVM+/0+RAQigtOErRRCc7juL63XxwMzEW21hRPde2uXtlX1H1ZbpNVh0gZamwhnM3CZBCLCL5NoIa6SNoQZgPR9ADcikYi8o6q/X7YgRFvl+HjTJSCiNQrTzHsBs0j4pFLAPqJkS1gz72R3c70+O2m+n123nYnzKc7C1EwLAD4RkR6AgbdPADwA8EvXBSOi6Nvfd5/ononzKY7C9Jm2AdQAXE/cVVbV77su2Cqwz5TWpt83PyM4f5WJ84lGNtFnWlHVVwEFCWr6JUq2e/fMT0YFokRYOJiq6isReQdmpZieqn4pIh96i4QTkd/du5suQSxx3BbF1cIDkETkAYAuzKhem976axH54SoKRhRr/f6oqTei7JzrZbcDhyOGEjZui7ZImNG8BVX9jqo+AvBbYJjknvlPiGLEdeJ8V0nzieIsTJ/p30/Zz04hohhxmTjfddJ8uwJNLuf0tEQrF6Zm+pGI/Jn3uwKAiLwP4CPHZSKKv1yOEeEW9vbMRhQ3YWqmzwB0REQBDHzfSB84LxVR3HW7my4BEa3RwjVTb5Hw7wB4AuA5gE9U9T5TCRIFaLfNBgCnp2YCpX9xhU7H7JusvabTZr9/8FKpxLRARBEXenFwVW2soiBEW8VlE+/5ubtzEdFKLJwBCQBE5GcwNdN3YVIKnqjq/1pR2ZxjBiSKJTtXpBS9NNiusykxAxKt29ozIInIMYD7AI5gEt7vAiiJiKjqr5YtCBFNEcEgSkTjwozm3VXVR6r6SlW/VtWWN+f0zqoKR0TJ8/Dh6PdSafrKM5Mt6UHHsKuZ1iVMMH09ZT+HLRKt0osXZkuA/f3xcVrLYj4JWpcwA5Cm1UDftb+IyI/Z5Evk2A9+YH4moCNxMp9EmPSCk5fHcT4JopnCBNPvi0gF4wuEZwD0RKQMk1YwB4DBlMglf7snLaxW23QJKEnCBNMMgI/nHFNeoixEFCQGTbwu0wru7+87SXnIcVu0TkuvZ+rHtU2JkmV/f995onsmzqc4CrWe6QLHbGxtUxEpwTRBZwC0VJWBnWjFXCbNB9zWcCM8PZe2UOgMSFEkIhmYqTt17/YlzHxYovhjJoNbKXudTgymtA5hpsZEWRHAle92dtqBRERErkUumIpIUUSqU+478e4vec261nsArieOTa2ynERrozqqlR4empqqf1BSvW72+atg/b7Zl06PnyuXS1w2A38SBz97KYI2/6W0axIk7LJRSJEJpiKSF5ETmBHBNwKhF2C7qtrwmnN3RaToO2RnTUUlireHD4G/+qtNl2Ll9vfdn5Njo2iahRPdi8iH3q92YE8JAFT1l04LZIJmSlXLE/vfquod3+08zAjjgheEB74+07FjLSa6J4o214nzXWLX9XZyleg+TM30JwAy3vqlr2CaVl+JyE+XLcQ8IhLUB3oNIO/93oBJGGGbd1urLhMREZEVJpg2VfXXIvIBgJyqPvGmwny9orL57WCiTxRmCTiISMqbBtPxaqslAJU1lIkonk5P3SbAJZrh4GCx/uigrdMZHWsXPYhqv3WYqTFvvZ95mJqgtY5GjxRu9ona4LoDXxMvWCslmu0XvzA/GVBpDVz3M0e13zpMMM2JyB2YWl8JAETkAdYz8GcQsM8+72SNdap+vx84Kfzp06c45QcLJcXTp5suQSxxuMVygvqac7nF+6DrdeD83G2ZXAqTAelTb4Hwsqp+KSI/gsk29HbOQ124xs0RvimvXEGBNlA6nUa/33dZLqL44RfHW5lcP5XIL1QGJFX1fy/oqepnjssz7Xm7IjIZNHfAJl0iokjLJiSFzsIDkETkxxO7BiLyQER+6LhM09Qn5pUWAHCRJaKwOp3xkR0RJCJOtgOHo1VKJaYmvI0YvN2cCDPP9KdBc0pdLQjuTX/JwyRt2AHwDCZhfdd3zAlGyez9g44WwnmmRIj0hMmDgwPnq8a4mrMa4cuWGLap3WVwdjXPdGYw9fpIcwDuwOS77U4ckgHQVtW/WLYg65BOp/Wbb74BwEFHlGCr+ESKINcJIBhMt9Nagqnvyf4aprZ4MXFXb5PLroXFmilRcjCYRkPUr9u6MyBVAFyo6mcTW2wCKRFNSKfNJ51/hLudGV/39aC8eGH2HR6OPz4oezzF3sHBzZHLsxIr+N8qds2FaYsLbLOFRvOq6rcAAkfuisj7qvqPLgtFRESbsYqkCK4WHYhyLXfhAUjDB4i8M7GrGpc+UzbzEiUHm3lvJ8qvcxVlW3uiexE5FpE/wCRpGPh+xmawuM2AJCIcfEREoWSz43Mml2n6jHKOWbqdMEkbdgHc8Zp8h0TkE7dFWh1mQCJyzPaj+hcr31KuBz9HNccs3U6YeaYPVPVVwP53JwNsVLGZl8ixCLcJBuXhXtb+/j5evny59HkifNkSV7ZNrGeqIvJ+wP7jZQtBRDH1+edmi6B9V6NefFwllFCNZrCi2wvTzPsTAN/1vu3ZPLkC4LsAbmRGIqIEmJwuEyEuapB+q6jpRtExq0e3EiaYZgB8jPHl0ATAidMSERHRxtRDJWldr1qEs7GHCaaVKX2mbxyWZ6X865kynSCRA/aTlxngQ0lIRkfnovw2CzXP1Jtj+ggmjeCXIvKhqv7DykrnGAcgETkW5dEqjrmctxrly2YDvKv1W10vXuBqEJjlagDSwjVTEXkAs+RZF8DfA/gSwNci8kNV/fWyBSGiGGIH29bZ88KKq0DvdhWg48hOKQrTzFtQ1e8Aw8AKVf1WktIrT0Q3RbmDjSLFZY0+isJMjfn7Kfsj2FBBRGvT75tPuXR6fH8uZ/b7OwZPT80+/3iFTsfsc9WuSLQBYYLpRyLyZ97vCpgk9wA+clwmIiKiWAmTAeldAB2YQOqfHvNAVX+/grI5xwFIRHRbSRmA5LpsUb9um8iA9LHXZ/oEwHMAn6jq/bgEUoCJ7okiL52+2Vy8hY6POXZr24Spmf4BZoWY53EKoH6smRJFXISrbK6XdIsq1kxvJ0zN9EhVfwXgvrcc2/eWfXIiojG/+53ZEsCOu5q2+cdtlUrTj3M9bqvdNhuFs3AwVdXPvJ+vVPUcQE9E/lZEfrqy0hFRsiSkmTfKfv7zA+ztybBLbNnNLfG26AmzOPiH9qeI/DWAHoBvYZI4EBFRCLncaPWYoM1f46zXpx/nfp1V91kRVrGCT9SESdrQEBHbUl1DwELhRERLsclXI5wMwlVty2VaPLd9iTXvXBFOhBtBYQYgfQWgHJTsPi44AIko4iI8AMl1jlnA3WAml5ctwn8CiJhquKq7juK15+bFlFVjiIicifAaW46Tqzs7V7JEN0tW6AFIccZ5pkQRVypFe50toimm1kxF5McAUjDZjtp2qTUvE1IeZrHwAoCvVPW/rKGsS0un0+j3+5suBhHRSoiYyr2/67lcnn68vyk3lwO6HE56a7NqpnUALVX9lX/NUlX9VlU/U9VPYdY2nfGnIiIK4cULs1Eo7gfLumvSToqZwXTewt+qOgBw7rZIRJRYP/iB2fxsdgK/w0Ozzx9463WzL4HNxC9fjqbK+F9+qTR7+o1fp2P3CYCHayz9dpgVTL+yv3gZj/5aRP41IEnDVyAicmF/H3jo6IO83zcb0RrMGs07/CroZTw6FxGo6i+nHUdEtJSgEbNBczSCmoInBy/duzf98RRTdv5x9FofZgXToHdgUC2U71Qiip67dzddAnLODtGJVzB9HDAXqhCw7zGAydoqEdFmsYmX1mhWMM0BuBOwf3fi9gfuikNERDRNdtMFmGpWMD1T1Y/nnUBEPnFYHiIioikcZ/V3aNZo3kXzekU3/9cEZkAiSpBczv1in0RTTK2ZqurXi5xg0eOigBmQiBKE6XxojRbOzUtEFCvtttkizNUC3CKCg4ODTb+cRGMwJaLtZJt5T09NZiR/106nY/ZNNgOn02a/vwWrVDIh7U2DAAASMklEQVT7HK6xuorFslexqDctLswSbERE5IDL5dwALukWBQsvDr4NuDg4EW0jG0xdfJ67PJdrq1i43NXi4GzmJSJaBFe0iYAcorpAOIMpEdEigla0ocUt2h/94oXZd3g4/ngRKLoAojlKm32mRESLcLWaDW0l9pkSEcXcKgYgRTE2iJgarKq7RPfsMyUiIgDup9o4PZ/TTFQlRHHFGIDNvEREsed6qo1TDjNRfQ7bjxq9gWAMpkREi1jFvIwkcNi1dogvnJ3LtUQ18zLRPRHRmjls5j3E5zjE507O5VqiaqZMdE9Et8Ya6cZ9gcP5B21IomqmRES0Zqen43mRtxSDKRERjRwcjGd6qtdNf3HJN4q23zf70unxx+ZyZn/Ht4j3L35hNgeO8Zc4xl86OZdriWrmJSK6NZuRZ9tTCv7mN8A//dPNDES39fQp8Pq1k1PV8T+83/67k/O5xKQNRESLSMpo3gh/aah7f4OSw78BkzYQEa3T55+bDVi+6fP01DSnRlGEE/qXvS2K2MxLRLQIV82egLM+xOSJbqsAgykRUVil0nitFDA10qDmR3+N1Hr6dDXl2nJ3Yac2pmcetwls5nWh0QCOjoA7d0xzTi5nbg8G7p/r7AwoFMx2dASUy6bJqdcLPzG60RiVudFwV8ZWK/h6FAqjCdyunq/bBXZ3zfOcnY321+vm+VstN8/j+nyUbFGeLiIy6h+OmD7uoY97my5GIAZTF4pF4PIS2PP6sMtlczuVcvccrZb5MG82zbntz1oNyGRMQAmbA7NYBM7P3ZXRyufHr8eTJ6MydzrA48cmuFYqyz9XNmvOPenqynyZcfWFxvX5iCi0Pu6ij7vDeO9fBtV2Y0/b/Gw3tsvvDAymcTAYmOADmIA0GaTzeRNU48I2j/lrkq5Vq8Dbt+YLQxTPR8nW6QQ3/0aBamRHLN/DOe4hmlnsGEzjoFIxAfXJk+nHlEpua8Kr5C+nwxUlZj5PFM9HybW3N2q5oRAeApBhvPd3W5dKo+8BQZtfp+P++wKD6bpUKqOmzaOjUU0TMIGyXB71gRYK432Ktp8uk5n9HF9/bX76+xFtU6ptJp7sW7SaTXNspWKef7IJdl4Zw/A3lWaz4323rZY5/+7uqA0n7HP3ejf7Uf2vv9Ew5z47G/XnAub25L5FzlevB5/PajTM/rOzUfn9fa+zXn+9PnruXG507WybVqEQ/vrT5mWzZqPtoaqJ2XK5nK5UPm++BNVqN/enUqPbJyfmuLdvze1MRjWbHd3f6Zj7r67MbfvlqtlcvCzNpnnMycnNfdXqaN/lpdlXKo0/PpUy5bbmlTGIvR6Xl+P7S6Wb16lWM/uKxdHz29/nPbe97X9dQfvs68/nR9fePm/QPn+5w57P/9hs1uyzisXxv/+812+f295WNY+d/JsRufDwodkiqA1o2/+/tCSvftpWB/GFNdNV63ZNbePRo9G+atUMaEmlTK2k1zODciz7jdX2g9oaaZjBLzs74co5ORI4nzfl7vUWK+MstZqpbdnaJWBev7+Nxpb3/n3z8+1bM7Bo2ef2s8+RyYyabO21DdrX693+fP7Hnp+bmr9la5P+2um01w+Y15vPj64FADx75mYAF9GkL74wWwTlvM2V42N35+I801WzH36T/W2TH7qvX483vxaLwHvvmd/zedOs9/r17AEwR0fBI1tvw36493qLlXGWo6Obc/KmmWzKXva5gwT1fe7u3u5c087nZ5uybd/3rCA9rSm/UjHBt1o1W683v9mf6DY+j+Z6ocAokLoaulWvu5vQwGC6avNqObaWdf8+cHISfEy1Cjx/bv7y1WrwMb3e/NpomJrt9bX56f/AnlXGVVnk+kTZYGBq/akU8OqV+Vmvh5+vms+ba2H7kWcNRqPosykHo7i+sstMT46tcLji0tjMu2q2ia7VGg9mZ2emCdh+SE42WbZao5pYKjWqcRYKNwNzr2eaUP2BNqi25G9qnHR1Nfp9MDDPXyyaYLpIGVdlk8/tgm0qf/x49DfxX+swbABttzl4Je6++cZstFEuZyclqmba7/ch3izdp0+f4tRVBpJGA7i4MB9ygPngbzZN+0EqNRop++CBGQ6fSpmalv1A7HRGo2htTXB3d7wmls+bfjQ7IhQw59nZMcdOBspMxgTXWs0Ex1Rq1C/67Jn5ac+fz4+P4O12zQe3//kXKaPVapnnnXY9go63ZXr2zASfRZ+72x2V++LC3J/JjO/LZs11svsaDfP4vb3RFxCb/CKVGu1zcb5i0TRxX1wAb96YpunHj0fNvtfX5rhZr9+yX26mtU5QfPzud5suwXS29WPRrpk1cp2E0eXsJC7BRhR1/v7RQmF2CwPRsqK81JzjspnTcQk2ou1XqYwGR9l5wLQdSqWbOfFevDD7Jvstg3LfHR6afa6XSzs+djvM1aFTb4uiRDXzEsXO48ejpuzdXdMkT7RK/uAeMXbhutNNFmIKNvMSEcVdhPs5XbJjXlzFLZfNvKyZEhHFnR2U6CKY2uk66eitGRrlMewMpkREceeyj/Oet15oBFstI7rODgAGUyKi+HPZz3n3rrtzOWaDqauUgu22u+kxDKZERDQSxaxMHhv3XNWZJ1OSL4NTY4iI4q7fN/2lk/2cuZwZZeNP9XN6avb5k9Z0OqNl/mJARJxtrjCYEhHF3fGxu4ztEba/v+/4jCFXnpqBU2OIiCgeHC8QwKkxRESUPBFeHIDBlIiI4iHCCwQwmBIRUTxEMJGExQFIRERES2IwJSKieCiVIpt/mMGUiIji4fzc8RQgdwkK2WdKRETxUHM3L9RwlEsQDKZERBQXEW3iBdjMS0REtDQGUyIiiocXL8zmjLsMgGzmJSKiePjBD8zPCKbB3apgKiIlAG1V7W66LERE5NjDh5suwVRbE0y9QFoG0Nt0WYiIaAWcNvG6tTV9pqpaBxDbJWFO/WsLEsUM378UY05yFK51CTYRKQK4r6qVgPtOYGqVO8AwOIY9fw3Apaq2gu6P8hJsIoIkLYdH24XvX4ojbwk2qOrSq4SvpZlXRPIAsgAKCGiGFZEqgKYNgiJSFZGiqjbWUT4iIooB8WJeBL+4raWZV1VbqnoGYNrAoNJEbbIJ0/8JwPSHishJwJYPU45+yAVlF226WuS4JDaDbfI1r+q5XZz3tucI8zi+d5fD967bc4R93Prev+6SQKy7mbcKIKWq/kCZBfBKVe9M7OuErXrPa+YVEQ3zehdtulrkuHnHbGMz2SZf06qe28V5b3uOMI9b53s3bNnigO9dt+cI+7h1vX/Fq+nGppl3jh0A1xP7BgAgIilVHSxyEq8/ds/7/ZrTY4iIaF2iEExT8AYd+djgugMvsM7j9a/O62P9fyLybwP2fwMgqA04LSKLtA0vcty8YxZ9rjjZ5Gta1XO7OO9tzxHmcet874YtWxzwvev2HGEft873778PUa6pohBMg4KlDa6TNdalqOq/c3m+VfOarSsAMgDA2jbFiYikADyCN+hwWvcLUdR47923GA2Y3QHwYNZncBSC6TVM7dQvBQCLNvFusT0AXwOoB00nIoq4S1UteB9M5wAYTCku9vz9qCJSmleZ2XgwVdWuiEwGzR3wHw8AnnF6EMWRN9Le1kgHAI42WyKixflbUbxpmnPzHmw8mHrqE/NKCwBcrwK7MUskq9jxzdFtqCpTJdLa3fL9m/XuzwLIA2ixm4I2wUGyoMkxPYHWlbTB/kMVYQLEFXz/XKpa8eaNFmH6B6+2oUa2bLIK+4cVkR6ASwC5dZWdyEGylR2v5akHoANgd01FJ3KSLMjL+b5QK+la55kmVdD8Wm//24n5tXkAFa+fqQgg4yW7uHEs0brc8v1bArBrawPeHO+l5/IRhXWb969vX0dVF6rEbE2i+7jxauuTrmFq8IAZ5dzwjs0gxkn8afss8P59Dm8UujcAiV0UFBkLvH+tzKLnjEqfaRLNS1bREpGi90e/Dw7goGiZm2xFRC68GmoKfP9StCyaLGjhL4EMppszN1mFr+0+9v3HtHXCvH+JomahZEGLNvECbObdpLUlqyBaAb5/Kc6cv38ZTDeHySoozvj+pThz/v5lMN0Qb1oQk1VQLPH9S3G2ivcvg+lm1b0pMNZWJaugrcf3L8WZ0/cv55mukC9ZRRnmW88zTGSC8WXgyMAM2pibtopoHfj+pThb9/uXwZSIiGhJbOYlIiJaEoMpERHRkhhMiYiIlsRgSkREtCQGUyIioiUxmBIRES2JwZSIiGhJDKZEMeKtDZpIIpIRkUtvsedNlqMkIrWJ7DmUcFyCjWLDWyS9DOAEQBfAhe/uXQCXqhqJ3LAiUgMAVS07Ol8KwDmAIgCZcdy8a1QCUI5xpqLXqnrm3+G95gqAK99u+z7IzFsKzguKTwBkAdQBVPzJzr01WWveOSuqWveek8GUhhhMKTZUtQeg4n34XQR8qF6KSCYigeISE0s5iUjptmXzPtyPRGRmyrIFrlENJsXaVvBeZxnA0UQAzMMEwLm1WFVtiEgLwFuYL2SDifvrIrKrqhW3padtwmBKW0NVj0RERaTlBZVNliWohrzwQsOroqpdEdnbdDlc8GqH5wA+CAiALRHpBj/yJlUdeAG1jImVQ7xWgdcOikxbjH2mtG1aWKA2skoikhKRvJdo296uYbT48CbKVPLdbG+qHI7VYBKXT1t/MuwKIFUAxYB+6UfzmoqJWDOlbdOEqV0AGNYqnsDULO4DaHq1ljzMh2fPe8y1d/+VvynWe3zJOw4wfXBnvvtt32QKoyWcBl4Z7L48TCDNescPV6eYVj7f+asY1YqWWXR7WCu2q2ZMXAMbeApTrsGsa9iGuYaPvftsn2LZ9xh/f+Y1Rl94jrzact47x1mI5lS7Ikggr4yZea9j4vgBzN/b3zye2EFfFIKqcuMWqw3mg/lkyn0l87YeOzYzcTvl/V6ECRz+x1/6zw2gY4/3bucB1HzPVfTdVwSQ9R3XnLjvcsprmVa+pj2fdzvjf20LXKMmzECk2rTH2WtgywATODREGUv2NswAHvv6O9OuhW/f5LUvzXlNGfu3seX0X/8Fr0ng6/Dtq/rL5ZUzO6ss3LipKpt5aeuk4NXg7NQFHe8/7QJ45Ls92bdag6m92Job1NeMqKYmU/JqOT0AVW+qhB01unBf7azyeU3Ee+pbe1HD9wM3VfVMzYjiYW1vohlz4D+3fa32mAWu4TWAnqoOVLXrK28Wo2vRAzDWT+tdR/iawvMAni/6wnx/kxtN597foyoiTW9QWnbB9wJg/v4ZWy6YQLpw3yslF5t5advsYtQnmAEwsEHR8xqzA14Po2a9LCZG5HoGMIGuJSJlmKbGmjfg5UGIss4q396ccoai46N68wD8fYCzmo8XuYZBj2/4nmdaoKzCfHE5gqkhhm3GbsE0S4+NkNZRE7rCTAOyzchz3wuq2rMDkUSkMuW1Ed3AYErb5hFGAc3WtvyjM+fNQ81g9AHaQ/CgoRSAnogUvdqorWVVcbO/7QavH8/WbAPL533oZ24+2okw573NNQTM/NYdr0Z4rcHzbZ/D1Oz9tdgwygA6IrJIIA7zOmowzf0dhKgtU7KxmZe2hohcAnhmm+W8QGeDlz0m5WvCA24GlmGTqPf41MTjiwAaXnNhZqKmc4Hper7nynhNolPL533o9/xlnSj3rXjNt4VFj1/wGga5r6p1VW3olJGwXgB8DqB6m6ZU729wDOCVv3xeGfPw1SrDvA5feXO3qC1TQrFmSrHhGyGaAfBYZJgI6D2Yml5Nb87vfADgiYgM5wlOfLj3vAA5gGnWbep4YoWc7/E7MM2RR959A5iAajPhZFT1zPuArgDYs4kavKbGtjdFxd90PKt89j5bW7b9mJcAjoM+6GdcI2CUAck2g9py2lHGdXj9xTA1xqoXsALL6D2+7L3OEwB1X5muROSt77X2YP4+k4G1hon+1DC8cnRhElXY537jPV8O41+W5r0X/M4w+8sR0RhRnZlQhWhr2ew5qrpwTY3m84LsY5hWgoFXG7YJFo79tVBfU/ki583AjN6d2Yy+DlEqC0UDm3mJyLU8TA5dO1J44AXQNmDm5k40jxPFHpt5KZEmmihvnTOXbvKauku+Ju0dmCbqptfcvQPTPM7rTluDzbxEFAte02oV00cHr6scJZj+2OaiTdS0/RhMiYiIlsQ+UyIioiUxmBIRES2JwZSIiGhJDKZERERLYjAlIiJa0v8HHwraZJC47EAAAAAASUVORK5CYII=\n", + "text/plain": [ + "<matplotlib.figure.Figure at 0x7f1179a3fa90>" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fontsize = 16\n", + "\n", + "fig = plt.figure(figsize=(7, 5))\n", + "ax = fig.add_subplot(111)\n", + "\n", + "ax.set_xlim(binning[0], binning[-1])\n", + "ax.set_ylim(1E-1, 40)\n", + "\n", + "ax.set_xscale('log')\n", + "ax.set_yscale('log')\n", + "\n", + "ax.set_xlabel(r'Deposited Energy [GeV]', fontsize=fontsize)\n", + "ax.set_ylabel(r'Events per 2635 days', fontsize=fontsize)\n", + "\n", + "ax.tick_params(axis='x', labelsize=fontsize, which='both', direction='in')\n", + "ax.tick_params(axis='y', labelsize=fontsize, which='both', direction='in')\n", + "\n", + "ax.step(\n", + " binning, null_exp, drawstyle='steps-pre',\n", + " linewidth=2, linestyle='-', color='black'\n", + ")\n", + "\n", + "ax.step(\n", + " binning, astro_p, drawstyle='steps-pre', label=r'$N_{\\rm{astro}}=12$',\n", + " linewidth=2, linestyle='--', color='blue'\n", + ")\n", + "ax.step(\n", + " binning, astro_n, drawstyle='steps-pre', label=r'$N_{\\rm{astro}}=4$',\n", + " linewidth=2, linestyle=':', color='red'\n", + ")\n", + "\n", + "ax.legend(loc='upper right', prop={'size': fontsize})#, bbox_to_anchor=(0.96, 0.96), markerfirst=False)\n", + "\n", + "ax.text(\n", + " 8E4, 0.2, r'${\\rm\\bf IceCube\\:Preliminary}$', fontsize=fontsize, color='r'\n", + ")\n", + "\n", + "matplotlib2tikz.save(\"plots/tex/syst_astro.tex\")" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAFDCAYAAABhmOTuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3c9zIumZJ/Dvs7O7hzl0U+qYQ9EROzbqjdidUxsk/wEu8EHV1RFrS1Vz3xaq2ev2iK65lPqkRvYed8MgH3cjRgXtCHe1HDEL5cMepwTtOY9FzR5K5cOURPsyhw372cObLySQIBLehEz4fiIyJJIk8wUhHt5fzyuqCiIiIprdv1p2AYiIiJKOwZSIiGhODKZERERzYjAlIiKaE4MpERHRnP71sguwSH/6p3+q//Iv/wIAuHv3LtLp9JJLREREy9Rqtf5ZVf9s3vOsVTD9i7/4C1xcXCy7GEREFBMi8n9dnIfNvERERHNiMCUiIpoTgykREdGcGEyJiIjmxGBKREQ0JwZTIiKiOTGYEhERzYnBlIiIaE4MpkRERHNiMCUiIpoTgykREdGc1io379XVFUQEAPD06VMcHR0tt0BERGus2+2iWq0ilUrh8vISAFAulyc+5uTkBJlMBtfX1wCAYrEYeTmnsVbBNJ1O4+rqatnFICIiAMfHxwPBM5fLoVqtjg2QpVIJhUIB+Xy+d7ter2N3d3ch5Z2EzbxERLQU9Xod1Wq1dzuTyaDRaIw9vlqt9gIpABQKBVQqlUjLOK21qpkSEVF8NBoNZDKZ3u1Op4NHjx4FHttut0f2bWxsoNlsRla+MFgzJSKKuXa7jVwuBxFBqVTq7e92u4FBJin8gdQ+j8PDw8Bjr6+vsbGxMbAvlUoBMK/DsrFmSkQUY91uF5VKBbVaDYAZoNNsNpHP59FsNpfWX3hwcDDVceVyuRf0gnS7XTx79gy1Wg2np6cTj7ODjiwbXK+vrydeYyFUdW22XC6nRLS+gPFbpdI/rlKZfKxfNjv+uP39/nEXF7OVuVarjeyreIUNui/Jstls77kNazQamkqlBvZdXl4qAL25uZn5mgAu1EF8YTMvEVGMBdU8NzY2ZhrF6h/sEwfDzbMHBwdja7wbGxsjx9vbS6+Vgs28RLRGVKc7rlg02zRaremOy+WmOy4q3W4XtVrN2bzMeZt5m80mCoUCbm5uRu7vdrsj+7LZ7Mi+6+vrgdG9y8RgSkSUMMNBpFqtIpPJoFar9aaKtNvt3qCds7MzFAoFdDod1Ot1ZLNZdDodlEollMtllMtlNBoNnJycIJvNotvtIpPJIJvNji3DvFNStra2UCwWBwJko9HA7u5ub1+n00G73e7VwIvF4kCNvNFoTB3Uo8ZgSkSUQHYkbKfTQaPRQK1WQ7vd7gWbs7MzbG9v94KurdnZQJTJZNDtdrG1tYVardYLyPb4QqEwcc7nvFKpFA4ODnBycgIAePv2LTKZzEASh2aziVqt1itzuVzGyckJ6vU6Op0ONjc3Y5GwAWAwJSJKHH/fYSaTwcHBAer1Ot6+fdur1T158gTHx8c4Pj7Go0ePAmuZqVSqd3yj0cCTJ08G7u90OgPTV1zLZrMTa7/FYnGkWXrc1Jll4wAkIqIEq9frvabQzc1NACYIPnv2DOVyGa1Wq5f31k4lsYkO/PM2t7e30el0erevr68jDaSrhsGUiChBhueWZjIZXF5eotlsYmNjA61WC91ut7ev2Wwi541+KhQKqNfr2NjYQLvdxsXFBer1OgBT4+t0Omg2m6hWq7cmnKdBotMOb1sBW1tbenFxsexiEBFRTIhIS1W35j0Pa6ZERERzYjAlIiKaE4MpERHRnBhMiYiI5rRWwfTq6goiAhHB0dHRsotDREQrYq2SNqTTaVxdXS27GEREtGLWqmZKREQUhbWqmRIRUbzY3LwvX77E9vb2rekCT05OkMlkeguFu1oFZ14MpkREtBR21RrLZmoaF1BLpRIKhUIvGX+pVJppXdcosJmXiIgWrtvt4r333hvYd3BwgOPj47GPqVarA0vPFQqFuZeCc4XBNC6ePwdEgAcPBveLmM3vwQOz7/nz/r5q1eyLSZMHEdEk19fXKJVKA8n1NzY2BlbE8Wu32yP7NjY2ekn7l43BlIgo5trtNnK5HEQEpVKpt7/b7QYGmSTIZDJotVoDK9M0Go2BmqefXejczy4fNy4ALxL7TOPiwQMgaNGBoH3+GqlVLPZrpVdXwPvvA3fvmt+tXA5ot4GLC/M7ABwdAZ9/Djx9an4HgFYL2NoCslnzOxEtTbfbRaVSQa1WA2AWyG42m8jn8yMryCzSwcHBVMeVy+Ve0BvmX8u02+3i2bNnaI35zOl2u71BR5YNrtfX12OvsSgMpvMKE4zSaeDNG+D1a/M7YALg11+bx7KJloiGNJvNgX7BSqWCarW6xBL1y+HS3t4eXrx4MXYN1aBgaYPrcI11GdjMO4tWy22N7c0b4Je/dHe+dNrUaIcTVLRaZr+tlQImiKv2vwgA5n5V1kpp9cwzBuHqyuyzX4StXM7s9/+/HB2Zff7/qxn/n4JqnhsbGzONYo1DEA5SKpVQKpUGaqrDgvpT7e1l10oB1kxns+UtfWeD0HBqQhuMhgVlX6pWzRZX9oODmaOIEq3b7aJWqzmbl+mimRcA6vX6wHSXdrsdGFSz2ezIea6vr8f2sS4ag+ksJnx7Wjlv3iy7BETuzDIGwbItPsOCapzjvmQ7MhxEqtUqMpkMarVar/m13W73Bu2cnZ2hUCig0+mgXq8jm82i0+n05nmWy2U0Gg2cnJwgm82i2+0ik8lMrCm6aOZtNpu952L7RM/OznrX7XQ6aLfbvRp4sVgcqJE3Go2pg3rkVHVtNiCn5r9hdKtUtKdSCT7Gbmvl9WuzEVFsVHwfWJeXl7q7u6uqquVyWWu1mqqqHh4e9n5vtVqqqprNZgfOk8lk9ObmRm9ubrRSqfSOV1XN5/ORPoebmxsFMLLZ56JqnudwOexzLJfLA6/DrABcqIP4slZ9pjm0oBjsL/kKD6AQ/Lt/6H87/Q//pwqFoIL+N9O7uIJC8Bpp3L9/v7f6jKvt/v37C3sdQkmnzVYsmj4gf5P0LHNj4/o8iRLE33eYyWRwcHCAer2Ot2/f9gblPHnyBC9fvkQulxs7FzOVSvW2RqMxMvjHPwfUtVQqFRiU7KhlwNREG43GwOMODw+xu7uLw8PD2KQSBNZ2AJL4tq8BAP/9f3zcC2z/83/ZZoNT33Hve/ve4Fe/+pXzEkVxztjJZoHHj5ddCqKVUq/Xe02hm5ubAEwQfPbsGcrlMlqtFi4vLwH0R73a4OofBbu9vT0QPK+vr8eOrKVRokF9ACtKRBw92Qu88867+PbbD5ycTbxa3Dr9LYhoNs1mE5lMphfo2u02KpUK9vb20O12e/2Itp8UMMG1WCz2ksTbx967dw+np6e9PkjbZ9rpdJDJZGIzuCdKItJS1a25z7NOH+BbW1t6cXEx93mGWzCt/f1+K6idajqOP2+CSBWm5vvR3GUDgJ2dHZyfnzs5FxHRKnMVTNe0mXc+Ozuuz1gE4K4vMbZNxnGfBkRENCPWTGPA1nRd/Cli3WTs8okSETnAmiklz/6+2QA32WQ4OpiIYoJJG2hxomji/d3v3J+TiCgkNvPGwNo08xIRxQybeQOISFFE1ijXHxERxcHKNPOKSBHAAYDoUnZEZJ1S/RIRraKVqZmqahVA/Npwp+Bf0c1m4gva/F2OdoWo4c2kt/x6vVIdEhEtWeyCqYjsikh5zH2H3v1FryZKgdwGv9jOW33zBvjIl+hi3vzBREQzik0zr4jkAWQBFBDQVOsF2IaqNu1tEdlV1fpiSxqtaccNBa0QNXQmF8XpDWiKpZ2d0WWuZmWDbdByXEREt4jdaF4vaKZU9WBo/42q3vHdzgMoqWrBt68CoGYD7rC4juaNs7UZHcyEEkRrydVo3tjUTCcZM0L3GsDqZ2Gmxfjqq2WXgIgSLBHBFMAGTPD06wKAiKRUtSsiuwC2vH3XqtpecBljwSbP9ycQoikM96kSEYUQuwFIY6RgAqqfDa4bAKCqdVXNqerBuEB6dXUVOFr1yFW/Wwy022YjIqLFSUrNtBuwzwbX4RrrWOl0GldXV25KRKvFjgCePKqLiCjQ1DVTEfnQ297xtk9F5NMoC+dzDVM79UsBgKoGBVqicA4OzGa5SMQf1/m5RORcmGbexwAyqvp7AC8AvAfgxSICqtdsOxw0NwAEjtolCm1np7+ijQtHR0Bc5+cSkXNhmnkbqvoLEfkugJyqbgOAiGSiKdqI6tC80gKAyoKuTavu/HzwdjodPE0maGTX0dHofNerK+aJJFojYWqmN97PPAB/ogQnE/NEJCsihwB2ATz0sh31Po1UtQQg42VAOgRwGTZhg38A0ioNOlqElU9N6Jo/RyQRrbwwNdOciNwBUAJQBAARuYfRUbYz8Zpy2wBOJhwz9r5prMMAJH9LZasFbE2Yinxx0Z9KUywCp6dBRymAcwAfBd0ZWmxTExIRzWHqYKqqPxGRfQAHqvprEfkxgAz6NVaKgSjW385m76PVmr8BItapCYmI5jB1OkEROVbVJxGXJ1JMJ7hca5OaEOiPAl7xlhCipFtGOsGSiFwCeOaN6CWicd68WXYJiGiBwgxAKqjqzwFsi8iPReTDqAoVFQ5AooV5/dpswPxLwz14wDmrRDEXps/0hf+nl8Dh72BWafl5ROVzah0GILnGxVRm5HLR8q+/dncuIorE1MFURD5U1d94NdK/gZki8wxAyxuMpKr6i4jKSZRc1eroyLAHD4K/oQTt44o2RLEXps+0LiL2P/0EwL6qfuvd/gYAROSTpNRSiRKDK9oQxV7YRPePbTPvMBH5gYPyEBERJU6YAUil4UDqDUSyQTQHoOOsZERkBDUTE1GsTD3PtPcAke8M7Sqp6l+5KlCU0um0vvGmLDx9+pQjeqfgcgDSWs0zdYmjwIgi42qeaZgl2L4nItcAWt7WBnDp/Z4I6XQaqgpVZSCl5Njf7+eJdLE0nAin2hA5FqbP9ABmtZhXIvJjVf0S6OXnpRVViWBdHpdpBXd2dnA+vOLLqomiifd3v3N/TqI1FiadoD+A9kbtisiPkjIlhukEl+v+/fuRJLpnszERzWrhzbzeRX/k/frK9/ui1jOlhDs/P+81s7vYiIjiIkww7QD4GxF5xxvV+1hE3gLYjKZoFAccSEpEdLvQo3kHHizyPVX9xmF5IsVm3vDiPJCUo4NnxBVtiHqW0sw7TFW/EZF35i3EojDR/exs/vWgPOx2IGnQViz2j2u1OIg0Ft684ao2RI6FzYAU5BTAIwfniRwT3Ye3swO4GjN0dOTuXDQHu5oNETkTWDMVkXdF5I8i8odbtj8C2F1wmWmBzs9NE+/w5tdqBR+jOtjfysVPYiKdNpuLpeE4Z5UIwJhg6iWwr6vqn9gNwA8BfDC0bwvAwwWWl4jiJJsFHj9edimIlm5Sn+n+0O13VfWVf4eqtgHcOC8VEUWvWjXNB/6Obbs03PPng8cGNUk8f26aJbiqDU1y/76brF253CJKO7OxwdS3vNptUo7KQkREq+b0tJ8Oc4WFyYD0BYC/VdXf+PZ9COAvVfWziMrnFKfGLJfraTYu0xICa5Ka0DXb3+qv3RIliKupMaHmmYrI/0Z/qbWM9/Oeqv5+3oIsAoPpcrkOplGkJ+Sc1ZDiPBGZVk8Ec6SXMs9UVX8IoACgCuChqm4nJZDS8l1cmM0Vl+kJaUZc0YZuUyy6a7mI8Rzp0PNMvUFH7QjKEjmbtAHgeqbLEPPxAzQLrmhDtzk9NT9dvFdiPEd6rnSCScNmXhqHqQmJIhLzfnVXzbwuMiARTcX+LzFxPtEaiWkQdW2u3LxEYZye9lt8iIhCc9n/6hiDKRHFh011SKvj+fPRJCCzivE3cjbzElF8xHSkJs3h44/NTxfjESqV+c8RkbE1UxH5TsC+fRH5rZfk/h9F5D9FWTgiWjOvX8d6xOZasOn//HM551kU4aOPzOZCjJt5J9VM90TkUlV/AQAi8tcANgGUAHS93/9KRMQeQ0Q0FzbxLt8337g9n6sm3pib1Gd66m1WR1Ufq+qXqvpCVateEoftaItIRGvF1dJwa/Ih7tzVlWmS9X+xmXdRBFdc9r86NimY3gGQ8jX3jnt1Oi4LRKsrmzWbZT8Hgzb/52i1Ov44JsshWiMff9zvg42ZSc28dtHva+/nhr1DRN7xpRFMzKoxzIC0XP5scq44Ts3rNHk+E+fPqFodnYxsa0HDgvbFtOZCDrjqe43ApJppB8BjAA9FZB+A+GqpXRE5FpFPYfL0JkI6ne7lYmUgXT7bGhS0+VuTisXJx7qws7Pj5kQ+rpPwU0jVqmmq9L+Z5s0fvA5yufjm/oxxM+/YmqmqfjnhcXe8Y6Zd85Qo1lzXIF0vD0cz+OUv3U61sYE2roHGlXYiU68v3cy5eYeaehOBuXlpUZjrdwWty3Jz6/KlwbOQJdhE5BMROROR/+rbty8ifwRwIyJ/O28BiGYV59YoWkHDI+hWVZz/sYJGb8fE2JqpiHwBII/+QuAKYB/AzwFUAIh3/6WqPllIaefEmulqiXNFgTVToghE8E+/kFVj/BcQkSyA8tBFq17QJSKiVWAHXMVx4FWMv5xOaub9rf+Gtyh4PeC4f3ZaIiIiWp7PPzcbhTKpZnoHMP2mqvpzb98ze6dvAFI8G7CJiFyy02n8OWtX0dOnyy5BIk0KprYJ98ci0lTVf7JTYUTkXZgBSBWY/lMiotW2LivaxLF517LpJGM413TSPNNvAXzmbSP3iUgBwAXnmhLRWuBqNsv39dfLLsFYM61n6jXxvnBdGKIw9veXXQJaK+uyok2c55l+9dWySzDWpKkx34Ople7b5Axes+8h+knvD1X1vy2ioC5wagwtCqfGUGLFec5ZBBaRtGELphnXBtK/BvAWQA7ABzBLr/2Zl583EWyiexFhbl5aCPt+m3e7z+Vxlm9dloZbl+QUjk1q5n2rqv71TLtDt18BaHtBNhHS6TSuVn0k3hqJc2vUzs6O00T3TJofA+vSZxrF8k6u2C8y/sULYmJiBiRV/cx3+8dBye9F5FNV/WmEZXSGzbyrZV1ao9hkvIJiHBRiLaEZkI5F5ALAzwA0YabC/ABe06+3HFt+3gIQEa2dgwPzk8E0nBiPOhzbZ+pNebkH03faAdDwthsR+QOAFoBUUmqltLps99TwcpS2iytoG24aFgHYLUkLs78f38CQTsd35HLQwvExcVtu3m9hFgh/LCIZAN8FcA3Tf/pqAeUjGmtnB3DZlchuSVqYmAYEAOuTnMKxmdczTSL2mdI4ce5/ZZ8pLZQdpBnH2mkEZVvUeqafishbEfmDiPydiPy5777vicgXXNOUiCikq6v45viNczPv+++bLYbGNvN6U14ewSRuuAbwQ5ipMD9Q1X9Q1W9E5FsA/wjgLxdSWiKiVWADAlsbwrl7d9klGGtSn+n2UNX3SxEpAXjmTZv5NUwSByIiCiPGQaE3wjiO/bpxrc1jcjB9ObxDVbsAfug17wJmRC8REYUR46CAUy83TxyDaYxN7DO1vDmlPV4yhzsAHrovEtHiVSpmI1qoXG50TtfRkdnnT3naagXP6UqnzX6XwZn/DDOZtATbT7zE9m8BfCEiOVX9je/+L0UkDy4OTiuAc+eJPHH+Z7BfJmKY8nBizdSrgVYBfOAPpL77mzA1VCKKmKuk+UycHyOtlhmE5K9xHh2Zff6aaS5n9g0Hkasrsz+dXo9E/O222WLo1mZeVf12UoIGLg5OqyDGiVWws7Pj/JxMnE+JdHFhthhi0gYixDtpg2tMAkHUt5CkDURERHQ7BlMiIkqGo6PBvuQYWatgenV11RuAcRTTPwgREY3x+edmi6GJq8b4icgnMGuZjozqTYp0Oo2rOE+WJiKi8Z4+XXYJxpo6mMLk5h0ZvSMi76jq790ViYiIKECMWxTDNPOewSwSPizGM3yJiIiiFyaYFgC0vKXYzrztGYAnEZWNaGFU+9NibIa3oM2fHMZmeBu3MS8CkWOtViyzHwHhmnm3AJzALMfml3JXHKLV8NFHsW6RIkqmLW86aAznSIcJpiVVfTG8U0SCmn6JEmvaL742wxsRLUg2u+wSjDV1MFXVFyLyDsxKMR1V/bWIfKiq30RXPCIiIk9Mm3iBEH2mInIPQBtmVK/9evBKRH4URcGIkqxYjPfiG4C7xPlMmk8UcgCSqn6gqg8BfAP0ktxzCTaiIaen/TWW48Z14nwmzScK12f692P2s9eIKEHOz8+dncsmzSdaiHTa/Ixh8p0wNdPvi8ife78rAIjIdwB833GZiIiIRr15Y7YYClMzPYaZZ6oAur5vpPecl4qIiGjY69fLLsFYYUbzfgvgAxHZBfBdmBG9X0ZWMiIiIj/bzBtDYWqmAABVrUdRECIioqQKtQSbiPy1iFyLyB9E5K2I/OeoCkaUZNlsrOeXEyVTjOechVmCbR/ANoA9mIT3mwCKIiKq+vOIykeUSDGeW06UXHa+WbW63HIECFMz3VTVh6r6QlVfqWrTm3N6J6rCERER9VQqZgOA58/NihIPHgweY1ea8HvwwOx7/ry/r1odPW4OYfpMX47Z33ZRECIiooli2sQLhAum42qg79pfROQTNvkS9b/wMhE+UUQePAj+Bwva56+RWrb/1VHtNEww/aGIlDC4QHgGQEdEDmDSCuYAMJgSEdFaCRNMMwA+u+WYgznKQkQJ5TKt4M7OjtOUh0RjOUxLOPd6pn5c25Rovezs7DhPdM/E+bQw77/v7FRTj+a9LZB6xyxtbVMRKYpI3vuZWVY5iPxEBkfx2wGE4za/XK6/P66rnJ2fn0NVnW1EC3X3rrNThUraEFde8Nz0putUAZSXXSZab45XOQMra0QRcNjMK6vwbVBEDgF0vUAKEblU1c3h47a2tvTi4mLh5SOaxzqNDLZ9r6vwuUTJICItVd2a9zyhc/NGzUukv62qpYD7DmFGE28AgA2eAN7D4ChjiEhKVbsRF5eIiCg+wVRE8gCyAAoYCoze/WUADVVt2tsisutLvL+xsMISEVHy5XLOThUmN++H3q820BUBQFV/6qIgXpBsish7AFIBhxSHaqsNACUAdQBvh47dYK2UVsX+/rJLQLSi2u4S+IUZgPQYQEZVfw/gBUzT6gsR+dRZacYQkaD1N64B5L3f6zAJIyAiKQDNqMtEtCjVaizzehMln8MxNGGaeRuq+gsR+S6AnKpuA72RtFHbgAmefl3v+ilV7YhIy9dUPNLfSkRENGAZzbwAbryfeZiaoLWIYXcpjPaJ2uC6Ad9IXrBWSivGLufm8P+eiBwL08ybE5Efw9T6fgYAInIPixn4E9T/aa87XGMd6+rqCiIysh0dHTkpJFEUtrbMRkSOOfzsn7pmqqo/8RYIP1DVX3uBNYN+jTVK1xgdlJTyyjX1QKN0Oo0rh5N0iYgowT7/3NmpQmVAUtVTX1rBjqr+ZBFLrqlqG6O10w2wSZeIiGb19KmzU00dTEXkk6FdXRG5JyI/claayapeQgerAKCyoGsT0QIFdcfMst2Pa1JjigeHzbxhaqYDzayq+sqrpTrpMxWRrJfhaBfAQxE59E+J8eaYZkRk1zvu0pewgYhWwI7jpMZcgYYWZWJuXq+PNAfgDsyUk+EZrhkAF6r6V5GV0KF0Oq1v3rwBADx9+pQDjygR1ik3r0vM80u3arUgW1tOcvNOleheRH4GUwM9G7qrs8xl18JiontKIgbT2TCY0q1EIMBCE92XAORV9ct5L0hE4fD7H1FEsllnKQWnCqaq+i2AwEAqIt9R1X9yUhoiGsFkDUQRabX6TT9zCr04uIi849/A1H1ERLTmwkyN2ReRP8Ikaej6fhYjKptz/gxIHHxESVEsmi/PQdtwrXXccSJMlk8UpTC5eTcB3PGafHtE5Au3RYoOMyBREr1+7fZ8NqgWE/M1mCgi6bSzU001mhcweXh92Y/8+98dDrBxxdG8ROszOlgc9YX57ezs4Pz83Pl5aUkcjuYN02eqIvKdgP1cupiIYsd1AgiASSBWjsNmnzDNvI8BfM/7tmfz5AqA7wH4qbMSERE54LoGGUVNl5bMYTNvmGCaAfAZBhPOC4BDZ6UhIiJKoDDBtDSmz/Stw/JEyo7mBZhOkIho7TkchTd1n6mqvvDmln4iIj8AABH5MEnpBNPpNFQVqspASkS0APfv33e2CpDzlYBOT52dKsw803swie5/CJP0HgBeLXAJNiIiShjXg7acnq/ibhXPMKN5C6r6gao+BPAN0EszyF55ogRRXf1pMRQ/tlVwns25ZTTzAvj7Mfv5b0lERGstTDD9voj8ufe7AibJPYDvOy4TERFR9J4/d3aqMKN5jwG0REQBdH1zru45Kw0RRc7m8221llsOoqX7+GNnpwpTM/1MVT8A8ATAMwBfqOq2qv7eWWkixkT3RGb5RkdLOBIl20cfOTtVqHmmInIJ4FmSAqgfE90TEVHP8+fO1jMNE0z3VPVLEbknIhkAl6r6ayelICIiSrAwSRu+9H6+UNVTAB0R+TsR+TSy0hER0UK5TrKwLsIkbfjQ/hSRnwHoAPgWJpEDERGtgChWxoliBR8nHAb7MM28dW8kLwBUELBQOBHRqnNV24r72qiRJElYYWFG8wLAY1X996r6UwZSomTa3zcbYKbHiIzf/NNnisXgY1ymSo0z17Urro0aAw6/MMy9agwRJUu16vZ86xITXNYi16kvcV2EHoCUZJxnSjQol+vn6g3abIIHwATh4fuzWbMRrbuxNVMR+QRACmYx8AtV/Y23/10AeZjFwgsAfquq/2UBZZ0b55kSucUsSpRoDx44O5WM62QWkT8CyNogOuaYFIC3qvonzkoUoa2tLb24uFh2MYhozdlm3jgO8lmrsolAgJaqbs17qknNvNVJgRQAVLULwN3qqkRERIvy1VfOTjUpmP7W/iIi+yLyMxH5Q0CSht+CiNaSHdFLlEgOm3knjebt/Yt4GY9ORQSq+tNxxxEREa2jSTXToEbpoFpo/BrWiYiIbuNwntikmumjgLlQhYB9jwAM11aJ6T4/AAARlklEQVSJiIji7eDA2akmBdMcgDsB+zeHbn/XWWmIiIgWZX8fOHUzhnZSMD1R1c9uO4GIfOGkJERERItUrToLppP6TCtTnmPa45aOGZCIiCgKY2umqvpqmhNMe1wcMAMSkVuVxHyVJgrgMB6ESXRPRDSgWFx2CZLNZcL7uC/pFkvvv+/sVGGXYCMiojlFsVg2l3Sbwd27zk7FmikRzcxO02MNNRzXNUgu6TajqytnKbwYTIloZnaaHoMprTs28xIREc2JwZSI5mYT3g+3mOVyg/f5N39tttXq779/f7FlpzWWyzk7FYMpEc0sgnE04DgaWph229mp2GdKRDO7bRxNqzXdeXI5QJXLuc2LA5FCurgAtuZeFxwAa6ZERInneqpNFFN3YslhMy9rpkRECcdkDcvHmikREa0nhzna1yqYMtE9UbxdXJiNaCE+/9zZqdaqmZeJ7onizWEXFtHtnj51FlDXqmZKRETUw2ZeIlpFxSJTE1IyMZgSUWycnpqNaCGmnQg9BQZTIiJaT44SNgAMpkREtK6yWWenYjAlIqL1xGZeIiKi+GAwdaFeB/b2gDt3TKbuXM7c7nbdX+vkBCgUzLa3Z1ZnrlaBTif8JL16vV/met1dGZvN4NejUDC/53LurtduA5ub5jonJ/391aq5frPp5jquz0dEK2WtkjZEZnfXbIWC+bA9OHA/vt8GqK0toFYDUqnB+zY3w59zd9f83NtzU0YrnzebfT2ePOlfCzBBb28PODwEyuX5rpXNmtdj+IvE5aX5MuPqC43r81Egh11YRLdLp52dijXTJOh2+wGv0RgMpIAJXJXK4ss1K/tFw1+TdK1cBm5uBoN4nM5HgVotp91YRJO9eePsVAymSVAqmYD65Mn4Y4rF0SAbV/5yOlycd+J14ng+Ilqu16+dnYrBdFFKJVO7tD/9Tavdrmkatn2ghcJgn6Ltp8tkJl/j1Svz09+PWCr1z2H7MINqhI2GObZUMte3j5u2jGH4m0qz2cG+W9tMvrlp+ilnuXanM9qP6n/+9bo598lJvz8XMLeH901zvmo1+HxWvW72n5z0y+/ve530/KvV/rVzuf5rV62afYVC+NefiAyHzbxQ1bXZcrmcRiqfVwVUK5XR/alU//bhoTnu5sbczmRUs9n+/a2Wuf/y0twGzNZoTF+WRsM85vBwdF+53N9Xq5l9xeLg41MpU27rtjIGsa9HrTa4v1gcfZ0qFbNvd7d/ffv7bde2t/3PK2ifff75fP+1t9cN2ucvd9jz+R+bzZp91u7u4N//tudvr21vq5rHDv/NVoB9qwdtQW+XcRutJgAKx39gABfqIL6wZhq1dtvUNh4+7O8rl82AllTK1Eo6HeDRo/79dhSG7Qe1NdIwg182NsKVc3gATz5vyt3pTFfGSSoVU9uytUvAPH//IC1b3u1t8/PmxgwsmvfafvYamUy/yda+tkH7Op3Zz+d/7Ompqflbtjbpr52Oe/6Aeb75fP+1AIDj49HWgxWws7PsEtBacThQlKN5o2Y//Ib724Y/dF++HGx+3d0F3nvP/J7Pm2a9ly8nD4DZ2+t/AM/Lfrh3OtOVcZK9venftMNN2fNeO0hQ3+cso6Ennc/PNmXbvu9JQXpcU36pZIJvuWy2Tuf2Zv8EOj+f7rhpEuLb74cc0ERjOUwEzWAatdtqObaWtb1tpooEKZeBZ89MQB03laTTub02GqZme31tfvo/sCeVMSrTvD5x1u2aT/VUCnjxwvysVsPPV83nzWth+5EnDUYjANGObaMVYVvNHGAzb9RsE12zORjMTk7Mf7v9kBxusmw2+zWxVKpf4ywURgNzp2PeEP5AG1Rb8jc1Dru87P/e7Zrr7+6aYDpNGaOyzGu7YJvKHz3q/038r3UYNoBeXHBCJpELbOadzdXVFUQEAPD06VMcuVoYtl4Hzs7MhxxgPvgbDdOEkEr1R8reu2eSLqRSpqZlPxBbrf4oWlsT3NwcrInl86YfzY4IBcx5NjbMscOBMpMxwbVSMcExleq3ex0fm5/2/DbBgu2Da7fNB7f/+tOU0Wo2zXXHvR5Bx9syHR+b4DPttdvtfrnPzsz9mczgvmzWvE52X71uHr+11f8CUqn0+z7tPhfn2901/7BnZ8Dbt6Zp+tGjfrPv9bU5btLzt+yXm3kTXRCRc2IGM62Hra0tvbAf8ERJ4e8fLRQmtzBQj/e9GWv0EbfybGXIWdx6/hzy8cctVZ17LTY28xLFWanUHxxl5wETkRsff+zsVAymRHH26JFphrdBNZ9fdomIVsdHHzk71Vr1mRIlTjbLZt0Z7e8vuwQUe8+f9/sD5sRgSkQryc4iIloENvMSERHNicGUiFYSl3OjWzlq4gXYzEtEK2rLm+zAqTG0CAymRESUKOKwRukKm3mJaKWJ9Dd/s2+xOHiffxteRMnuv39/sWWnQTsxXlaIwZSIVlIUn7u/+pX7c9L0zs/Pna5x7RLTCRIRTYHpCVfPcxF8DDCdIBER0aweODwXgykREa0lBlMiIqI5fe3wXJwaQ0Q0heH16Yn8GEyJiKZQLC67BOTaPoBTR+diMy8REa0ll2shMJgSEU2hWuVKNKvG5Z+T80yJiKbAeaarx0tLyHmmREREccBgSkREa+muw3OtVDAVkaKIZJddDiIiir8rh+damWAqIkUABwA2ll0WIiKKPwbTAKpaBZDY0UVHR0fLLgLRzPj+pSR63/xIuzjXQkfzisgugG1VLQXcdwigA69m6QXHsOevAKipajPo/jiP5hUR50sCES3KOrx/OZp39dhFxlV17tXGF5IBSUTyALIACjABc/j+MoCGDYIiUhaRXVWtL6J8RES3YRClSRbSzKuqTVU9AdAec0hxqDbZgOn/BNAbWHQYsOXDlOPqKlwL+bRNV9Mct47NYMt8zlFd28V5Zz1HmMfxvTsfvnfdniPs4xb1/nXZTrnoZt4ygJSq+gNlFsALVb0ztK8Vtup9WzOviGiY5ztt09U0x912zCo2ky3zOUV1bRfnnfUcYR63yPdu2LIlAd+7bs8R9nELe/+KQJCgZt5bbAC4HtrXBQARSalqd5qTeP2xW97v16o6rhZMRBRaLmd+tlrLLQe5k3N4rjgE0xRGp7PY4LoBL7Dexutfva2P9f+JyL8J2P8GwaOk0yIyTdvwNMfddsy010qSZT6nqK7t4ryzniPM4xb53g1btiQY+3xk7jrM7NeOwXkX8d4Nc7yL9+Z/DFGuseIQTIOCpQ2uwzXWuajqv3V5vqh5zdYlABkAYG2bkkREUgAewht0OK77hShuvPfuDfoDZjcA3Jv0GRyHYHoNUzv1SwHAtE28K2wLwCsA1aDpREQxV1PVgvfBdAqAwZSSYsvfjyoixdsqM0sPpqraFpHhoLkB/uMBwDGnB1ESeSPtbY20C2BvuSUimp6/FcWbpnlr3oOlB1NPdWheaQFAZZkFcmmOZBUbvjm6dVUdmaNLFLUZ379Z7/4sgDyAJrspaBkcJAuaKkXtopI22H+oXZgAcQnfP5eqlrx5o7sw/YOXq1AjmzdZhf3DikgHQA1uB58RTeQg2cqG1/LUAdACsLmgohM5SRbk5XyfqpV0rRYHX5ag+bXe/puh+bV5ACWvn2kXQMZLdjFyLNGizPj+LQLYtLUBb4539ONgiYbM8v717Wup6lSVmJVJdJ80Y5aKu4apwQNmlHPdOzaDBCfxp9Uzxfv3GbxR6N4AJHZRUGxM8f61MtOeMy59puvotmQVTRHZ9f7o2+AADoqXW5OtiMiZV0NNge9fipdpkwVN/SWQwXR5bk1W4Wu7T3z/Ma2cMO9foriZKlnQtE28AJt5l2lhySqIIsD3LyWZ8/cvg+nyMFkFJRnfv5Rkzt+/DKZL4k0LYrIKSiS+fynJonj/MpguV9WbAmOtVLIKWnl8/1KSOX3/cp5phHzJKg5gvvUcYygTjC8DRwZm0MataauIFoHvX0qyRb9/GUyJiIjmxGZeIiKiOTGYEhERzYnBlIiIaE4MpkRERHNiMCUiIpoTgykREdGcGEyJiIjmxGBKlCDe2qBrSUQyIlLzFnteZjmKIlIZyp5Da45LsFFieIukHwA4BNAGcOa7exNATVVjkRtWRCoAoKoHjs6XAnAKYBeATDjutteoCOAgwZmKXqrqiX+H95xLAC59u+37IHPbUnBeUHwCIAugCqDkT3burcla8c5ZUtWqd00GU+phMKXEUNUOgJL34XcW8KFaE5FMTAJFDUNLOYlIcdayeR/ueyIyMWXZFK9RBSbF2krwnucBgL2hAJiHCYC31mJVtS4iTQA3MF/IukP3V0VkU1VLbktPq4TBlFaGqu6JiIpI0wsqyyxLUA156oWGo6KqbRHZWnY5XPBqh6cAvhsQAJsi0g5+5ChV7XoB9QBDK4d4rQIvHRSZVhj7TGnVNDFFbSRKIpISkbyXaNverqC/+PAyylT03bxYVjkcq8AkLh+3/mTYFUDKAHYD+qUf3tZUTMSaKa2aBkztAkCvVvEEpmaxDaDh1VryMB+eHe8x1979l/6mWO/xRe84wPTBnfjut32TKfSXcOp6ZbD78jCBNOsd31udYlz5fOcvo18rmmfR7V6t2K6aMfQa2MBTGPMaTHoNL2Bew0fefbZP8cD3GH9/5jX6X3j2vNpy3jvHSYjmVLsiSCCvjJnbnsfQ8V2Yv7e/eXxtB31RCKrKjVuiNpgP5sMx9xXN23rg2MzQ7ZT3+y5M4PA/vuY/N4CWPd67nQdQ8V1r13ffLoCs77jG0H21Mc9lXPka9nze7Yz/uU3xGjVgBiJVxj3Ovga2DDCBQ0OUsWhvwwzgsc+/Ne618O0bfu2LtzynjP3b2HL6X/8pX5PA5+HbV/aXyytndlJZuHFTVTbz0spJwavB2akLOth/2gbw0Hd7uG+1AlN7sTU3qK8ZUU1NpujVcjoAyt5UCTtqdOq+2knl85qIt9S39qKG7wduqOqJmhHFvdreUDNm139u+1ztMVO8htcAOqraVdW2r7xZ9F+LDoCBflrvdYSvKTwP4Nm0T8z3NxlpOvf+HmURaXiD0rJTvhcA8/fP2HLBBNKp+15pfbGZl1bNJvp9ghkAXRsUPS8xOeB10G/Wy2JoRK6nCxPomiJyANPUWPEGvNwLUdZJ5du6pZyh6OCo3jwAfx/gpObjaV7DoMfXfdcZFyjLMF9c9mBqiGGbsZswzdIDI6S134SuMNOAbDPyre8FVe3YgUgiUhrz3IhGMJjSqnmIfkCztS3/6Mzb5qFm0P8A7SB40FAKQEdEdr3aqK1llTHa3zbC68ezNdvA8nkf+pnRRzsR5ryzvIaAmd+64dUIrzV4vu0zmJq9vxYbxgGAlohME4jDPI8KTHN/CyFqy7Te2MxLK0NEagCObbOcF+hs8LLHpHxNeMBoYOk1iXqPTw09fhdA3WsuzAzVdM4wXsd3rYzXJDq2fN6Hfsdf1qFyz8Rrvi1Me/yUr2GQbVWtqmpdx4yE9QLgMwDlWZpSvb/BPoAX/vJ5ZczDV6sM8zx85c3NUFumNcWaKSWGb4RoBsAjkV4ioPdganoVHZ3feQ/AExHpzRMc+nDveAGyC9Os29DBxAo53+M3YJoj97z7ujAB1WbCyajqifcBXQKwZRM1eE2NF94UFX/T8aTy2ftsbdn2Y9YA7Ad90E94jYB+BiTbDGrLaUcZV+H1F8PUGMtewAoso/f4A+95HgKo+sp0KSI3vufagfn7DAfWCob6U8PwytGGSVRhr/3Wu14Og1+Wbnsv+J1g8pcjogGiOjGhCtHKstlzVHXqmhrdzguyj2BaCbpebdgmWNj310J9TeXTnDcDM3p3YjP6IsSpLBQPbOYlItfyMDl07UjhrhdALwAzN3eoeZwo8djMS2tpqIly5py5NMpr6i76mrQ3YJqoG15z9wZM8zhfd1oZbOYlokTwmlbLGD86eFHlKML0xzambaKm1cdgSkRENCf2mRIREc2JwZSIiGhODKZERERzYjAlIiKaE4MpERHRnP4/562EKt2Zs/MAAAAASUVORK5CYII=\n", + "text/plain": [ + "<matplotlib.figure.Figure at 0x7f11796f0e50>" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fontsize = 16\n", + "\n", + "fig = plt.figure(figsize=(7, 5))\n", + "ax = fig.add_subplot(111)\n", + "\n", + "ax.set_xlim(binning[0], binning[-1])\n", + "ax.set_ylim(1E-1, 40)\n", + "\n", + "ax.set_xscale('log')\n", + "ax.set_yscale('log')\n", + "\n", + "ax.set_xlabel(r'Deposited Energy [GeV]', fontsize=fontsize)\n", + "ax.set_ylabel(r'Events per 2635 days', fontsize=fontsize)\n", + "\n", + "ax.tick_params(axis='x', labelsize=fontsize, which='both', direction='in')\n", + "ax.tick_params(axis='y', labelsize=fontsize, which='both', direction='in')\n", + "\n", + "ax.step(\n", + " binning, null_exp, drawstyle='steps-pre',\n", + " linewidth=2, linestyle='-', color='black'\n", + ")\n", + "\n", + "ax.step(\n", + " binning, dastro_index_p, drawstyle='steps-pre', label=r'$\\gamma_{\\rm{astro}}=3.0$',\n", + " linewidth=2, linestyle='--', color='blue'\n", + ")\n", + "ax.step(\n", + " binning, dastro_index_n, drawstyle='steps-pre', label=r'$\\gamma_{\\rm{astro}}=2.0$',\n", + " linewidth=2, linestyle=':', color='red'\n", + ")\n", + "\n", + "ax.legend(loc='upper right', prop={'size': fontsize})#, bbox_to_anchor=(0.96, 0.96), markerfirst=False)\n", + "\n", + "ax.text(\n", + " 8E4, 0.2, r'${\\rm\\bf IceCube\\:Preliminary}$', fontsize=fontsize, color='r'\n", + ")\n", + "\n", + "matplotlib2tikz.save(\"plots/tex/syst_astro_index.tex\")" + ] + } + ], + "metadata": { + "anaconda-cloud": {}, + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.15" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} |