{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Plotting\n", "Here you can explore the different possibilities that the hep_spt package offers for plotting." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import hep_spt\n", "hep_spt.set_style()\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from mpl_toolkits.axes_grid1 import make_axes_locatable\n", "from scipy.stats import norm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting a (non)weighted sample\n", "Use the function \"errorbar_hist\" to plot the same sample without and with weights. In the non-weighted case, we will ask for frequentist poissonian errors, so we will get asymmetric error bars for low values of the number of entries." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAFCCAYAAADVI1hLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dv29Tedvn8c+FIiKRle7cibZaiSBnAW038nBLIKRkJOKKCinz0KSOJRoeKIBmiqUCCtihQHKev4BnkKiQRjhIIRLVPXf2LqaBIV6QtpgmudMkkheL7xb+HuM4/hUf2+fX+yVF2OfY51wh8TfXuc73hznnBAAAgMGciDoAAACAJCOZAgAACIFkCgAAIASSKQAAgBBIpgAAAEIgmQIAAAiBZCrGzCxnZiUz227aNu23lcxsOsLYtnvsz/d6zTHPlzOzspnlh3XMELEM9XsDss7MlszsX2b20D8PPu/bZpbz2/Jm9g8zW+5wjEy2SbRH8UAyFWPOuYqkXySVzKzkt+1JeijpoX8cVWzzPfZvSap02m9mq8c8X0XS1nHeMyq9vjcAx+OcW5e0JmnHP69Iutv0OPjc3XXOvehwjEy2SbRH8UAylQDOuUeSLpjZUtSxSI0rsrZXh32+f1pScYC37gx6TgCx91zS9abnM5L2gsqPbzfaJg20SYgayVRy/Cip1G6Hma36Uu9S0KAEpV+/banbbUEzW/Yl9ryZ3fGl9Gm/vWxm/9MfY9WX3HdVr44F73/o37vaHIPf1+78OUnB8fNNr73Tch41HTMv6W8d4s/7Y+WC2wRN5877+JpvFbT+v+SaXpdveu+/2h2jzfmPxA3geHyFZbrlM9ScYF0IqlRtPnO0SV1ixuhNRB0A+uOcq/gPWUmHG41lv3/LP39oZhXn3JaZVSRV/Hvzki5IWm9z7Bdmdl3Snuql9uvOuT3//leSTvoyvMzsF+fcj36f/Id12p/vgn/cXIY/cn7/2r3m15kvsTefx3+v003bOpWyr0v6uz/P86btRR+rVL/qvNvh/2XJObfmX3dd0pZzbt3/Pwbn3lP9luv3zSduF7fqiS+A43shadnMXkj6zX/9H/lbflL7zxxtUueYRXs0FlSmEiS43Sep+XZfQYdL39v+NYHd1uP4hCvoxB5cuZRU/3BfUL2P1pLqZfb/oXqpPedfW26JqSIp6IRZadOf4cj5O/i+zXlav7e2nHN3JeXM7B86XKq/67+PC5Jaq3K7HR5Pt9vuv892V3nt4gYwmKASlXPO7fl+oRX/OQ7agq6fuYy3SbRHESGZSp4f1VSZkvQPHf5Afa/61ZykRof1Q5xzd51zRf8VdO5c17ck7T/1rQEo+/0V/9r/bBeUc24ruBpq2d6pk/yuVC9ddznP33W0sZhpPZCZrTrnHjnnmhuSJdWv+tblq3HNJe+WuHp25LfO/TX6+v8B0Ftwq0+HE4jnqg+4CT5/tEnq2CbRHkWEZCrG/AftrpndCbY1j3Lxz9dUv9cf9Aso+7JxXvUro1X/oStI+tG6T6ewrvqV3J6kXefcenBV5+//L/ljNo7t31cMSuDBtj7O/0vT+9XuPE3bgv4JeX+u1u9huul92/7/KGgY86o3zDOtsbeJ67rqHf0bDZx960OxKl8ubz5Gu7i7/P8C6K2kw90R1lRPqCS1byuy3CbRHsWDOeeijgEJ5huf9aaOocuSZnySl2hmVnbOFaKOA0D/aJMQBSpTCOs3SXmrjy4JRsEcKa0njf9eLliI4dYAIkGbhLGjMgUg1Zr6wBR8x+CgWrEnKe8HdgDAwKhMAUgtfzVf8J1+D1Ur/La9puoFAAyEZApAavkRXcGAjZwfLRbMqSbVR0PFYmUBAMkV2aSdZva/JMk59+9RxQAgG/yI2GC6j2kdnsdntuW1G30c8r9LekH7BUCKdgb07xYXFxcl3YwwBgDjZVGc1Dn3yA+V/633q/vyXxYXF2+K9gvImrZtGMvJAEitpv5RW6rf0ltV/RZfMNHitFoWq3XO/dDHcTckLQ4xVAAJRp8pAGkWLIsk1ROniuoTQAaTGeaUgmHzAKJFMgUgzdZUnx06WBD8RdOi4EuS9oLnADAobvMBSC2/NFIw8/WLpu2Jnw0bQHz0TKaCdX+kw2sViQnvAAAA+rrNV/RJVI4J7wAAAA7rmkz5CtS2VB9azIR3AAAAh/W6zfc3qXGrb8nf0us64Z1//UYf5/6uzxgBAABiq5/bfDtNo19YrRoAAKBJr8rUtr5VoSqqV6q6TngnMekdAADIjl6VqXUdntzu72LCOwAAgIaulSnnXMXM9ponvJMkM7vAhHd4Uv6gn9/8cWT7zStndatwLoKIAKA/tF8YJnPORXNis43FxcXFjY2NSM6P4Tlz75Uk6dODqxFHggSIZKHjYaP9Sg/aLxxT2zaM5WQAAABCIJkCAAAIgWQKAAAgBJIpAACAEEimAAAAQiCZAgAACIFkCgAAIASSKQAAgBBIpgAAAEIgmQIAAAiBZAoAACAEkikAAIAQSKYAAABCIJkCAAAIgWQKAAAgBJIpAACAEEimAAAAQiCZAgAACIFkCgAAIASSKQAAgBBIpgAAAEIgmQIAAAiBZAoAACAEkikAAIAQSKYAAABCIJkCAAAIgWQKQKqZ2ar/eti07WGwL7rIAKQFyRSA1DKzJUnrzrk1STn/XJJWzWxbUiW66ACkxUTUAQDACOX815rqiVPOb//RObfe7g1mttHHcb8bSnSIzH61ptLmduN5/v5rrVyaU3FhXlOT/GnE8VCZwsD2qzU9Lr9vPM/ff63H5ffar9YijAr4xjm35qtSkpSX9Fvw2MyWzOxORKFhxJ6UP3Tct1+t6dqzdyq9/VaY3D34otLbiq49e9ezDet2bGQTyRTa6tVY0BghScwsL6nsnNuSJOfcI1+Zmm269Se/74deX5L+Of7vAsfx85s/Ou4rbW7r886BqrWvh7ZXa1/1eefgUMXquMdGNlHLRFs/v/lDtwrnOu7vpzG6XTg/8PGBIVtyzj2SJDNbliTn3AtJO/p26w8pc+beq2O/p1r7qqdvPurpm48jiAhpRTKFjgZpiCQaI8SLma02JVJLqvedCkqq85JKUcUGIB16JlNm9tA5d9c3SGt+27KkPUn5oJECgLjxydNDM7sraUb1judbfqqEXUnbwa0/pM+nB1fbbs/ff63dgy8d3zczdVJbPxU67h/0QnMYnpQ/tL3NePPKWar9EeqnMrXqk6ei1Oh7IOfcupnlzCxPY5ROnRoiKdmNEbLD94v6a5vta21ejoxYuTSn0tvKkW4KkjQ5cUIrF09HEFV/bhXO6VbhXKMN7dZOY3z66YD+o3NuvmkY8XXVq1JSvVS+1P5tSLOVS3OanGj/6xP3xghA+t28crbjvuLCvOZmTx1pwyYnTmhu9pSKC/MDHxvZ1E9lKm9m0rdbetOSdpv2z7a+gXlakq9XY1FcmNevv/95pBM6jRGAOOh2y2tqckIvb1xWaXO70bdzZuqkVi6e7mueKW6noVXPylS3IcRIr16NRdAYFRe/DYSamTqp4mJOL29cpjECEGtTkxOHRhxv/VTQ7cJ5JuzEQLr+1nQYQrynekdOqV6l2ml9n5+HpStfvVo8VrSIlaAxCq7suvWRAgAgrXpVpiqSgr5S86rPHvxc3+ZlyTXtBwAAyJyulalOQ4jN7IK/5bfHSD4AAJBlPW8OtxtCzLBiAACAOtbmAwAACIFkCgAAIASSKQAAgBBIpgAAAEIgmQIAAAiBZAoAACAEkikAAIAQSKYAAABCIJkCAAAIgWQKAAAgBJIpAACAEEimAAAAQiCZAgAACIFkCgAAIASSKQAAgBBIpgAAAEIgmQIAAAiBZAoAACAEkikAAIAQSKYAAABCIJkCAAAIgWQKA3tS/qAz9141np+590pn7r3Sk/KHCKMCgPTar9b0uPy+8Tx//7Uel99rv1qLMCpMRB0AkutW4ZxuFc5FHQYApMaT8oeO7ep+taZrz97p885BY9vuwReV3lb06+9/6uWNy5qa7PxnvduxEQ6VKQAAYuLnN3903Ffa3NbnnQNVa18Pba/WvurzzoFKm9sDHxvhUJkCACBGmrtP9Kta+6qnbz7q6ZuPI4gIvZBMAQAy50n5w6FKTZDA3LxyllthODaSKQBA5sS5z+enB1fbbs/ff63dgy8d3zczdVJbPxU67h+k4oX+0GcKQKqZ2ar/eti0bdnMlszsTpSxAcexcmlOkxPt/2xPTpzQysXTY44IAZIpAKllZkuS1p1za5JyPoHKS5Jzbl3SXvAciIObV8523FdcmNfc7KkjCdXkxAnNzZ5ScWF+4GMjHG7zIXNa+0oE6CuRSjn/tSap4h8XJJX9/oqkJUlbkUQHtOjWBk1NTujljcsqbW43OprPTJ3UysXTKi7Md50WodexEQ7JFDIn6CsR9B/o1D8ByecrUoG8pOeSvpe027R9tvk9ZrbRx6G/Cx0cMICpyQndLpxvJFPd+khhfLjNByD1/K28snOOChSAoeu7MmVmd5xzj/zjZUl7kvLBNgCIsaWmtmpP0ox/PC1pp/mFzrkfeh3MV68WhxgfgATrqzLlO3H+zT+m8yaAxDCz1aYLwSXVb/Xl/O6cpPWoYgOQDoPc5ruu+pWd9K3zJgDEjk+eHprZtpn9S5KCW31+3x63/gCE1fM2n5nlnXPrZlb0m6bVpfMmAMSFr6D/tc32tTYvB4CB9NNnaqb3Sw5jNAwAAMiKrrf5gqpUy+aunTcBAACypFdlKmdmuabHwTwtF4JtatN5k9EwAAAgK7pWppxzL5xzL1SvRE37bXTeBAAA8PqaZ8p31lxreQ4AAJB5zIAOAAAQAskUAABACCRTAAAAIZBMAQAAhEAyBQAAEALJFAAAQAgkUwAAACGQTAEAAIRAMgUAABACyRQAAEAIJFMAAAAhkEwBAACEQDIFAAAQAskUMme/WtPj8vvG8/z913pcfq/9ai3CqAAASUUyhdR5Uv7Qcd9+taZrz96p9LbS2LZ78EWltxVde/auZ0LV7dgAgGwimULq/Pzmj477Spvb+rxzoGrt66Ht1dpXfd45UGlze+BjAwCyaSLqAIBROHPv1bHfU6191dM3H/X0zccRRAQASCuSKQAAEuJJ+cOhCnlw4XjzylndKpyLKqzMI5lCKn16cLXt9vz919o9+NLxfTNTJ7X1U6Hj/kEqXgAwLLcK50iaYog+U8iUlUtzmpxo/2s/OXFCKxdPjzkiAEDSUZlKudaScCDuJeEwcd+8crbjvuLCvH79/c8jndAnJ05obvaUigvzAx8bAJBN5pyL5sRmG4uLi4sbGxuRnD9rgttTnW5/xdUo4t6v1lTa3G50NJ+ZOqmVi6dVXJjX1CTXFyNmUQcwDLRfQGa1bcP4y4HMmZqc0O3C+UYy1a2PFAAAvdBnCgAAIASSKQAAgBBIpgAAAEIgmQIAAAiBZAoAACAEkikAqWdm+ZbnD/2/q9FEBCBNSKYApJqZLUn6j5bNq2a2LakSQUgAUoZ5pgCkmnNu3cx2Wzb/6JxbjyQgAKlDMgUgi/JmJkl559yjqIMBkGwkUwAyJ0igzKxgZkvNVSoz2+jjEN+NKjYAyUOfKQCZYmbLZrbsn+5IykUZD4DkozIFIGsq+tbxfF5SqXmnc+6HXgfw1avFYQcGIJl6VqbMbMl/PWzatuy33RlteAAQjq9CXQiqUc65LUn/5p9v++cAMLCulSk/N0vBOXfXzO42z9XiR8jkzCxPYwQgrpxzLyS9aNm2FlE4AFKoazLlk6QgUco557Z8harst1UkLTW9RhIdOAEAQHb01QHd384r+qfTkprnbJkddlAAAABJ0VcHdOfcIzP7xcx+6/P1P/R6DR04AQBAGnStTJlZvqmfVEXSqqQ9STN+27TqQ4sBAAAyqVdlqrk/1LSkv0tal3TBb8v55wAADNWT8gf9/OaPI9tvXjmrW4VzEUQEtNcrmVrTtyHEwagYmdkFv3joHiP5kDStDfSZe68k0UADcXOrcE63Cucan9FPD65GHBHQXq/RfHuqJ1RS09BihhUnx361ptLmduN5/v5rrVyaU3FhXlOT8Z2zdZRxBw00AADDwHIyCfak/KHr/v1qTdeevVPpbaWxbffgi0pvK7r27J32q7WBjx1Wt+OHibvXsQEAGDaSqQRr15egWWlzW593DlStfT20vVr7qs87B4cqP8c9dljdjh8m7l7HBgBg2OJ7nwd9CfoSHFe19lVP33zU0zcfhxxR/waJPQ5xAwDQjMoUAABACFSmEq7b6Jb8/dfaPfjScf/M1Elt/VRou2/QitdxdIo9TNzSeGIHACBAMpViK5fmVHpbOdL3SJImJ05o5eLpCKLqLalxA4DE/FhZxG2+BLt55WzX/cWFec3NntLkxOEf8+TECc3NnlJxYX7gY4fV7fhh4u51bAAYtVuFc4cq758eXNWnB1dJpFKMZCrBen0wpyYn9PLGZRUXc41tM1MnVVzM6eWNy13naxr1h77b8cPE3evYAAAMG7f5Um5qckK3C+cbo9+69TWKk6TGDQDIHipTAAAAIZBMAQAAhEAyBQAAEALJFAAAQAgkUwAAACGQTAEAAIRAMgUAiKX9ak2Py+8bz/P3X+tx+b32q7UIowKOIpkCAETiSflDx3371ZquPXun0ttKY9vuwReV3lZ07dm7nglVt2MDw0YyBQCIRLv16wKlzW193jk4skZntfZVn3cOVNrcHvjYwLAxAzoAIDJn7r069nuqta96+uZjY4UEIGpUpgAAAEKgMgUAiMynB1fbbs/ff63dgy8d3zczdbLrmp2DVLyAQVGZAgDEzsqlOU1OtP8TNTlxQisXT485ov4xCjF7SKYAAJG4eeVsx33FhXnNzZ46klBNTpzQ3OwpFRfmBz52WIxCRCuSKQBAJG4VznXcNzU5oZc3Lqu4mGtsm5k6qeJiTi9vXNbUZPdeKt2OHRajENGKPlMAgFiampzQ7cL5xqi9bn2kxo1RiGhGMpVyT8ofDl3pBA3AzStnR3rlFlZS4wYAZA/JVMrdKpxLZPKR1LgBZAOjENGMPlMAUs/M8i3Pl81syczuRBUT0inJoxAxOJIpAKlmZkuS/qPpeV6SnHPrkvZaEy2gl6SOQsTokEwBSDWfNO02bbouac8/rkhaGntQSLSkjkLE6NBnCkDWTOtwcjUbVSBIpziPQsRokEwBQBMz2+jjZd+NOg4AyUEyBQxR65QOAaZ0iJU9STP+8bSknQhjAZACPZMpM1v1D+edc3f9tmXVG6S8c+7RCOMDEiWY0iEY3txp+DQi9VzSBf84J2m9eadz7odeB/DVq8VhBwYgmbp2QPejYNadc2uScn4oMSNhACSGv/i74P+Vc27Lb1+StBc8B4BB9apM5fzXmuqjXnKSCpLKfn8wEobGCEAsOedeSHrRsm0tonAApFDXZKqlwcmrXh7/Xj1GwtCBEwAAZEVf80z5W3llyuEAAACH9Tuab6mpo3nPkTB04AQAAFnRszJlZqtBIuU7bD5Xve+U1GYkDAAAQJb0M5rvoZltm9m/JEbCAAAANOvVAX1d0l/bbGckDABgpFonwQ3mb2MSXMQNM6ADAGIpmAQXiLu+RvMBAACgPZIpAACAEEimAAAAQqDPFAAAQ5TkjvOtsQeSEHuUSKYAABiiJHecD2IPEsBPD65GHFEycJsPAAAgBJIpYIj2qzU9Lr9vPM/ff63H5ffar9YijAoAMEokU8AxPCl/6Lhvv1rTtWfvVHpbaWzbPfii0tuKrj171zOh6nZsAEB8kUwBx9CuY2agtLmtzzsHqta+HtperX3V550DlTa3Bz42ACC+6IAOHFPQMfM4qrWvevrmo56++TiCiAAAUaIyBQAAEAKVKeCYOg0Vzt9/rd2DLx3fNzN1Uls/FTruH6TiBQCIHpUpYEhWLs1pcqL9R2py4oRWLp4ec0QAgHGgMhUDzDibHDevnO24r7gwr19///NIJ/TJiROamz2l4sL8wMcGAMQXyVQMMONscnRLbqcmJ/TyxmWVNrcbHc1npk5q5eJpFRfmNTXZ/eNG4gwAyUQyBQzR1OSEbhfON5Kpbn2kgHGh+g2MFskUAKQc1W/0a79aOzQnXv7+a61cmuurup5ldEAHACADeq2yEGYVh6yv4ECa2SfK5ACAJPv5zR9d/171s4rD7cL5gY6ddiRTfaJMDgBIukHns2MVh+64zQcAABAClSkAADKi212VMKs4ZH0FBypTAACAVRxCIJkCACADeq2yUFyY19zsqSMJVT+rOGR9BQeSKQAAMqDXaLtgFYfiYq6xbWbqpIqLOb28cbnrPFNZHskn0WcKABAC08akC6s4DIZkCgAwMKaNAbjNBwzVk/KHQ6Naztx7pTP3XmV+dmBEa79a0+Py+8bz/P3Xelx+33VGawD9ozIVA6yFlB7BVfoocDsFnTwpf+j4OxAsEfJ556CxLVgi5Nff/+zZF6bbsQHUUZkaMdZCwrDcKpw7dAvl04Or+vTgKn/oBmBmD/2/q1HHMgztkuxAP0uEDHpsAHWUPUaMtZCAWFo1s2VJxagDGZZBJk1kiRBgOPpKpsws75zbanq+LGlPUt4592hUwaUFayEBsfOjc2693Q4z2+jj/d8NNxwg/dLcVaFnMmVmS5IeSvreP89LknNu3cxyrYlWWtGvCUiVvJlJKbog7DSKLswSIRLLhGB40jzys2cW4JOm3aZN1yWV/eOKpCVJiU6menWwDNuBk7WQgHgJEigzK5jZUnOVyjn3Q6/3++rV4sgCHKKVS3Mqva0c6UogsUQIMCyDdECfltScXM0OKZbI9OpgGbYDZzeshQSMl5kt+64KkrQjKdft9UnQbSmPMEuE9Dq2xLQLgDSiDuhJ7HMQVb+m4sK8fv39zyPJGmshASNT8V+SNC+pFGEsQ9Gtsh4sEVLa3G60UzNTJ7Vy8XTobgpMuwDUDVKZ2pM04x9Pq35lhwGxFhIwXr6P57/56tR2Fvp8BkuEBLZ+Kuh24XxfiRTTLgC9DXJJ8lzSBf84J+nIiJgk9jmIsl8TayEB4+WcW4s6hiRh2oXsaB1xF/zs0zDibpT6Gc23LOmCmS07514457bM7IIf5bcXp6u6UQ27pAMnACALRrmKQ5r1M5rvhaQXLdtieVU36LDLXn2P6NcEIMuYdgHDkOYphlhORr37HtGvCXHAqCnEDaOREcj60mkkU30K04ET6Ee3BiNMQ9Tr2IiHJ+UPOnPv1ZGvqH92UU67gOQY5RRDSRioQCYAxES3tRbDrOHY69iIh7jODj3KaRf4nUyXLC+dRjIFxAijppA0jEYGSKYAAMAQZHnpNJKpGGBeDwQYNQUgjdI+xRDJVAwwrwd6SXtDBCDZsj7FEKP5gJhg1BRGJRgpGIjLSEGkR9anGKIyBcQEo6aybZQTGlL9RhykebBCqipTo5zUkCs7RI25zpKNecSA9EpMK/yk/KHrlVXQGH3eOWhsCxqjX3//s2sZsdexJa7sAITDPGJAeiUmmerVWIRpjGiIAIxDGucRYzQykKBkSsr27KoAEEdU7YGEJVMAkGTMIwakU6KSqSzPrgogvZhHDEi21IzmW7k0d2QOngCNEYCoMY8YkF6JSab6mV110MaIhghJwPQcydbPPGKDTGjY69hAXKS5DTPnXDQnNttYXFxc3NjYGNoxg0nvBpnUEMBYWNQBDMMo2q9A8MemW7cGAJFp24alKsNI8+yqwCi1Dm8PMLwdAHpLVTIFYDDB8HaqIgBwfCRTAAAg0aKurpNMAQCARIu6up6Y0XwAAABxlKrKFGtEAUgq2i8guVKVTLFGFICkov0CBhdMjRTI33+tlUtzY5saidt8AAAgtnpN6rlfrenas3cqva00tu0efFHpbUXXnr3TfrU28LH7RTIFQPvVmh6X3zee5++/1uPy+66NEACMQ7tRes1Km9v6vHNwZG3Lau2rPu8cHKpYHffY/UrVbT4A7T0pf+h4Cym4qvu8c9DYFlzV/fr7nz2XM+l2bAAYhuZlaI6jWvuqp28+NibzHhUqU0AGdLv6CnNV1+vYAJAFVKaAjBjkym5cV3UA0E23eaPy919r9+BLx/0zUyc7Li83aMWrFZUpAACQWCuX5jQ50T6dmZw4oZWLp0ceA5UpICM6XdmFuaqThndlBwDt3Lxytuv+4sK8fv39zyPdFSYnTmhu9pSKC/MDH7tfVKaAjIvDVd24mdmymS2Z2Z2oYwHQXa8BLlOTE3p547KKi7nGtpmpkyou5noOoBnW4BmSKSADul19FRfmNTd76khC1c9VXa9jx5GZ5SXJObcuaS94DiC5piYndLtwvvF866eCbhfOj2XCTmnA23xmtixpT1LeOfdouCEBGLZuV1/BVV1pc7vR0Xxm6qRWLp7ua/bgBE6LcF1S2T+uSFqStBXsNLONPo7x3fDDApBUx65McVUHpM/aZuXQiL3d/f+np28+am2z0uVdiTUtabfp+WxUgQAYjiflD4f6b56590pn7r0a2gznvQxSmep6VQcgeVgX7hvn3A+9XuOrV4sjDwZAX6JuwwZJpnpe1VEmBxBje5Jm/ONpSTsRxgIgBeiADiBrnksKhv3kJK1HGAuAFBikMtXzqo4yOYC4cs5tmdkFM1uStOeco5sCgFAGSaaeS7rgH3NVByBxnHNrUccAID2OfZsvuIrjqg4AAGDAeaa4qgMAAKijAzoAAEAIJFMAAAAhkEwBAACEYM65aE5s9n//8pe//LfvvmPuTiAr3r59+7Nz7t+jjiMs2i8gmzq1YVEmU/9b0n+V9LHXa4ckaPX+OabzEUN8Y4j6/FmO4Z8pSaZov4iBGLIZQ9s2LLJkatyCJW76mVCUGNIdQ9TnJwYcVxx+VsRADMTQGX2mAAAAQiCZAgAACIFkCgAAIASSKQAAgBBIpgAAAEIgmQIAAAiBZAoAACAEkikAAIAQMjNpJwAAwChQmQIAAAiBZAoAACAEkikgQmZ2J+oYAGAQtF/fTEQdQNaZ2R3n3KOo4xgXM1uWtCcpH9X3bWar/uG8c+5uFDH4OJYk/S3C8+cl5STJOfciqjiQbFlqw2i/DsVB+9Ukk5UpM8ub2bL/YLB69CkAAAIlSURBVEQZR9S/jKv+6+GYzpeXJOfcuqS94Pk4+f/zdefcmqScf55VRd8I5aL4WWBwtGGN84+tDaP9ip1YtV+ZTKYUsx9CFCL6UF5X/apOkiqSomgIck3nrfjnY2dmed8oR8L/Ed6WJOfcI+fcVlSxYCC0YeNvw2i/PNqvozKXTMXlhxD1L6Oi+VBOS9ptej47hnMe4pxb842vJOUl/TbuGLyZiM4b+JukWV/hoN9DgtCGNYy7DaP9+ob2q0XmkinF54cQ6S9jjD6UkfBX8+Uo/hDF4I9QYCf4/qO+XYRjoQ1Tttsw2i9JMWu/stoBfcc5t2VmS2a2PIrOa02dBJtVnHPrMfplHPeHck/fGuBpSTtjOGcnSxF2ms2ZWa7pcT6CRnFb366yK6r/gY68Eyf6RhvmjbENo/2qo/1qI5XJVLdGQGP6ITRdMbUzll/GHv8PgXF+KJ9LuuAf5yRF0hib2WrwPZvZ0rj/KAR/+PzPZ3qc526yLim4mstJ+ntEcaAN2rC6mLVhtF+i/eokc8vJ+AZg2Tn3yJfIK1ENq/S/jHcl/RhRuXY1aDDH9aH033NFUq5HYz2q8y9J+kX1P0Yzqv/fx+IKe9z8z2JX9Z9FJoa2pwFt2OHzj7MNo/2Kj7i1X5lLpqT4/RCiwIcSSC7aMNowxEsmkykAAIBhyeJoPgAAgKEhmQIAAAiBZAoAACAEkikAAIAQSKYAAABCIJkCAAAI4f8Dk91TgMGILgoAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Create a random sample\n", "size = 200\n", "smp = np.random.normal(0, 3, size)\n", "wgts = np.random.uniform(0, 1, size)\n", "\n", "fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(10, 5))\n", "\n", "# Make the non-weighted plot\n", "values, edges, ex, ey = hep_spt.errorbar_hist(smp, bins=10, range=(-7, 7), uncert='freq')\n", "centers = (edges[1:] + edges[:-1])/2.\n", "\n", "ax0.errorbar(centers, values, ey, ex, ls='none')\n", "ax0.set_title('Non-weighted sample')\n", "\n", "# Make the weighted plot\n", "values, edges, ex, ey = hep_spt.errorbar_hist(smp, bins=10, range=(-7, 7), weights=wgts)\n", "centers = (edges[1:] + edges[:-1])/2.\n", "\n", "ax1.errorbar(centers, values, ey, ex, ls='none')\n", "ax1.set_title('Weighted sample');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculating the pull of a distribution\n", "Sometimes we want to calculate the distance in terms of standard deviations from a curve to our measurements. This example creates a random sample of events following a normal distribution and overlies it with the original curve. The pull plot is shown below." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm4AAAHXCAYAAAAFhdYzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdfXxbZf3/8dfV3QWWjbE1gICwtQKDqUA3sChSvlJuVETADpAb/Yqw2YogsGYDsUmrMEjGnV8VNkDufghjQxFExVUcqF8LYxWRG2V0wBcEGd0oNIzurtfvj9OkaZs26U1yctL38/E4j6Un5ySfpU3OO9c513UZay0iIiIikv+K3C5ARERERDKj4CYiIiLiEQpuIiIiIh6h4CYiIiLiEQpuIiIiIh5RsMHNGHODMeYGt+sQERERGSlj3S4giw6pqKioAC5yuxARERGRDJh0GxRsi5uIiIhIoVFwExEREfEIBTcRERERj1BwExEREfEIBTcRERERjyjkXqUiIiIFZ+vWrbz55pts2bKFzs5Ot8uRNIqKipgwYQJ77rkn48ePH/7jjUBNIiIikgNbt27ltdde48MPP1Ro84jOzk4+/PBDXnvtNbZu3Trsx1OLm4iIiEe8+eabbN++nZ133pm99tqLMWPGYEzaob/EJdZaduzYwb///W82b97Mm2++yfTp04f1mGpxExER8YgtW7YAsNdeezF27FiFtjxnjGHs2LHstddeQPfvbzgU3ERERDwifnp0zJgxLlcigxH/fY3E6W0FNxEREY9RS5u3jOTvS8FNRERExCMU3EREhikcDmOM6bOEw2G3SxMpKJFIhEgkwrHHHut2Ka5Rr1IRkWEKh8OJ8AZOTzIRGVnNzc2sWbOGFStWuF2Kq9TiJiIiInlv/fr1lJSUABAMBl2uxj0KbiIiIpL3Nm3a5HYJeUHBTURERHKmsbGR0tJSIpEIy5YtY/bs2bS1tQFOq9rChQtZuXIlCxcuTKxvbm5m1apVNDY2EolEaGtr63fbwT5+Y2Mjs2fPTqxfuXIl8+fPZ/369Yma29raEvc3NjbS3NycuK+/OrLGWluQC7C6oqLCiojkCmCdj1WR7HjhhRfsCy+80Gc9uLsM1jXXXGPLysqstdauWLHCvvvuu9Zaa0tKShK3165da6uqqhL7LF261M6bNy/x80DbDvbxly5dmtg+/nMwGEz8XFZWZltaWhL7lZSUZFRHsv5+d72kzTdqcRMRcYl6o8poFr9eraqqiilTprBs2TJKSkqYMmUKAGVlZTQ2NqbcN5NtB/v48e0Bpk6d2qMFL/n+srIyVq1alXEdI029SkVEuoTDYerr6/usD4VCWQlT6o0qI8WLfzrJQQmgpaUFoEfwOe2001Lum8m2g338qVOn9tg+fk1dc3Nzn/vijz2YmkeKgpuISBcFKRH3HHbYYaxfv57KysrEuuTbQ912OPuA04q2fPnyEX3M4dCpUhGRYYrFYoRCocTPgUCAUChELBZzsSoRb6mqqmL9+vU9Lu7v77TjYLYdzj7gBLF4Z4i4ZcuWDesxh8OVFjdjTNBaG+m6XQW0AWUDrRMRyRfxljlwQlt5eXnilAlAa2srkUiEBx54gKamJvx+f5/9REarxsZGli9fTltbG6WlpcybNy9x34oVK1i8eDGHHXYY0H1Ksrm5maVLl9LW1kYkEiEYDPa77WAfP/mxV65cSUlJCUuXLuXpp59m5cqVVFVVsXbtWhYuXJiYsSG5Va2/OrLF5PpUgDGmEphvrZ1rjCkDSqy1K40x84Cnuzbrsc5a29zvA/b/PKsrKioqVq9ePXLFi8iokO5UqTEmcV8oFCISidDR0dFnO5/PRzAYTFw3l7zfYJ5PJO7FF18E4MADD3S5EhmsDH93aWejd/tU6ek4LWsA64HKftaJiOSVeA/QhoaGlKENoKOjg4aGhsS2IiLDldPgZowps9Ymn/ydAiQPhTytn3UiIiIio16ur3Gbmn6T9IwxqzPY7JCReC4RGb22bYO//x1eeAHWrYO2NvjgA4Cf8L3v1bDXXnDZZQHee6+138cIBAJs2LABQK1uIjJsOQtuKVrbwDklGg9zU4CNXbdTrRMRyTonmH0NOIVp06C9PdVWNVx5ZfdtiACpr3Grrq7u97lisRjRaDTxcyAQoKamhtra2kSHBhGRZLlscSsxxpQk3S4DlgNz4uuAeLBLtS7BWnt0uifrapWrGEa9IjKKvPEG/OAHMe64Iwr8Brib9vZpTJ1aQ0VFLR//uJ/iYpg4ER544BHKy7/I+vXw9NO1PP/8A0ALyeFt/HgfpaWl1NbWJtZVVHR/JGXaG1U9UUUkWc6ucbPWrrTWrsRpTZvSta4ZEj1N26y1zanW5apGERldwuEwGzfCd74DJSUxli0rZ+vWCNCKM+1oK5s3R3jppXKCwRgXXgjf/CYcfvga6urgjjvguef8vPZaE6ecEgQm43ysBti6NcjMmU28/XZ3y9njjz+euB2NRmlpaenTsaGjo4OWlpZES1yqmRxEZPTK+XAguaLhQERkINZCUdHXKS6+k9ZWgBBFRRE6O4c7rMe+XHzxq/z0p7BlC0yYAHV1UFsL48cP7Rq3Qv2clsHTcCDeVSjDgYiI5Nw778CXvgQQD22PAVemDG0w2GE9XuO66+Cll+BrX3PC2/e+B87YnDNH8H8hIqORgpuIFJz4fKO9l3A4zOrVcPDB8Mgj4Iw8dA5wDLBjRGvYZx+4805YtQpmzHB6p8LT3H23xVpLcXHxgPsHAgG1tElWDfQ+kfyl4CYiBSccDvcIPdY6YekjHwlTWQlvvQVHHQVwMNbePaggNdgwVVkJzz4LZ50FMJFzzoHvfheqq2vw+Xwp90nXG1VkJPT3PlFwy28KbiJS8Do74eKL4Vvfgh07IBiExx4DeCOxTU1N9oKU3w933w0nnvgQ48bBjTfC2rW1zJhR2uc5fb6evVGTJ68XKUSNjY2UlpYyf/58t0vxBAU3ESlwYzj3XLjhBhg3zukJes01MGZMz1BUW1tLaWn6IAVDC1PGwMMPn8SqVbDrrvCb3/jx+5u46KJgYptAIEAwGOwzMb1INsRisR5/y4FAgFAoRCwWy2kdlZWVLFy4cFD7LFu2LEvV5D8FNxEpYOOAn3PnnbDzzvDb38LXv959b3Io8vv9NDU1EQwOHKR67zdYFRXwv//rXAO3Zo2fxx6rB3YBYMOGDdTX12vwXcmK5L/b+DiCkUgksS4+jmB5eXmf8JZPXyDa2tp6jH842ii4iUhB2r4dnDG+T2PyZPj97+GYYwbex+/39xg3LVtBauZMePxxmD4d1qwB+AOw64g+h0hvyX/bmY4jmGpfN7W1tXH++ee7XYarFNxEpOB0dsJ55wGcAmziscfgM5/J3vPFe+fFZdI7b/p0J7yVlgLMBh5h8+bs1SgC3X+bDQ0NfUJbXO/hb7Ixx25bWxsLFy5k5cqVrFy5sk8LWnNzc+K++fPn09bWBjjXw7W1tdHc3EwkEqGxsTHtPoVGwU1ECoq1zmC3d94J8AHwBWbPzu5zxnvn9V7SnV7aZx9wxgh/DTiC00+PtxSKFLZjjjmGyy67jKqqKqqqqvrcP3fuXEpKSqiqqmL27NksXrwYgKqqKo499ljKysoIBoNUVlam3afQKLiJSEG57jpnGTcO4FTgSZcrGtjeewMcD2zk17+G+fOd8CmSDfEvFYMZ/makxxOMt5JNmTIlsa7UaXpOWLVqFWVlZQDMmTOH5ub0s18OZR8vUnATkYLx6187rW0Q46STQsDvAfd6y2XuX8AX2Wkn+NnPoEAbCiSPZHP4m3Sam5uZOnXqgNuUlJSwcuVKli1bRmNjI5s2bUr7uEPZx4sU3ETE88LhMM89B1/9KlgbY7fdynnkkfS95fKppxw8yX33OcOGXHGFE0JFRtJQh7/pve9wlZWVpQ1Vs2fPpqSkhHnz5vU4HdrbypUrB72P1ym4iYjn1df/DyedBLEYzJoV5f33M+st53ZPud6dGr78ZYO138NaOPNM6JqTWmREDHX4m977Dlc8VCV3Hli7dm3idrwDQvy05/r16xPbNzc3U1JS0qfjQbp9ComCm4h4WmcnwF288grAGp5/PvPecr0NpXfocKTq1NDZeSVz50J7O3z5y/D++1l5apGcDX+TyooVK1i8eDGNjY2JVrP777+fZcuWUVlZSVlZWeKUZ0lJCSUlJSxbtoyysjKqqqrYtGkTy5YtS1wnl26fQmIKdRJjY8zqioqKitVOly0RKVCRCDiDrm8EDgVeH9T++fgZ+MEHzvAlf/87nH463HuvcwpV5MWuZtgDDzxwxB4z/mUlH98LhSTD313ad7pa3ETEs/7yF7j8cuf2ww9Pw9r/y9pk8bk0cSLcf78zx+ny5XDrrW5XJCL5QsFNRDyptRXOOMOZNB4inHiis97N3nIjaf/9YelS5/aFF8Kzz7pbjxSeXF8aICNDwU1EPMdaOPdceOMNOOIIuOKK7mvasj1ZfC6deSZ885vQ0eGE1A8/dLsiKSRDHTha3KXgJiKec8cd8PDDsMsucN998IMf1CXuy9Vk8bnyox85c5u++KIzTEi8laT34oX/i4ycfD7VL32N5O9LnRNExFNeew0+8Qmn1+Xdd8PZZ/e/baFcdL1mjdOy2NkJf/wjVFQUzv9NBudf//oXnZ2d7LfffowdO9btciRD27dvZ926dRQVFXHAAQcMtKk6J4hI4ejsdE6RtrfDKafAWWe5XVFuHHaY0wnDWvjGN5z/v4xOEyZMAODf//4327dvV3DPc9Zatm/fzr///W+g+/c3HIrrIuIZP/0pPPYYBAJw882ja4iM+GwKf/sbLFjgdjXilj333JPXXnuNzZs3s27dOrfLkUEYO3Yse+6557AfRy1uIuIJr7wC8cvWli6F3XZzt55cGz8e7rrL+XfZMoDCndJH+jd+/Hj23XdfdtppJ4qKdAj3gqKiInbaaSf23Xdfxo8fP+zHU4ubiOQ9a6G62ulV+dWvOqdJR6OPfxzC4fjYdTcBn3C3IHHF+PHjmT59uttliEsU10Uk7913Hzz6KOy6K9xwg9vVuGvBAqdzBnwM+L7L1YhIrim4iUhe27QJLrrIuR2Njr5TpL1t2RLjsMNCQAC4gl13DRAKhYjFYm6XJiI5oOAmInktGIR33oGjjnJ6lGaikEaET645FotRXl7Oz38eAVoBS1tbK5FIhPLy8h7hzYv/VxFJT+O4iUje+tOfnMA2frwz4frMmW5XlHvGmMSQD6FQiEgkQkdHR5/tfD4fwWCQ+vr6PvuJiGek7Suv4CYieWn7dpg925mjs64OuvLIqGOGMeZJoX6+ixQwDcArIt50yy3xidVfpaFhJ03vJCKChgMRkTy0aZMz4CzAypXTqapyTg2O1hak+P87EAjQ2tra73a77hpg06YNwPBa6kQkf6nFTUSybrATo9fVOeHtc5+DU0/Nba35rKamBp/P18+9PiZNqmaUZluRUUPBTUSyLhwO92gts9ZirU0Z3J59Fm66CcaMgRtvHF3TWqUSCoUSt2trayktLe0T3nw+H2PGlPJ//1fLihV99xORwqHgJiJ5w1q48EJnMvmaGmemgNEuOdz6/X6ampoIxuf+wjl9GgwGuf76JsDPggWwebOGAxEpVOpVKiI5E7/uqr/PnQcegKoqmDYN1q1zZkrIZL/RqPdrsmMHHH44NDc702KpwU3Ek9SrVES8Yds2uOwy53ZDQ3dok8zETy2DM8PEf/7jbj0ikh0KbiKSF265xWll239/OP98t6vxpiOPhJNPhg8+cFrdRKTwKLiJiOva27uDxuLFMG6cczsWi/W4yD4Q0Lyc6Sxe7LS+3XorvPii29WIyEhTcBORrEsXwJYsceYj3Xvv1znllO59ysvLiUQiif1aWzUvZzozZzotljt2wKJFblcjIiNNwU1ERlyqidH7C2AvvxxjyRJn/RtvnJEY/iMajdLS0tJnXs6Ojg5aWlqIRqOJdfWjdT6sfoRCMHEiPPQQPPGE29WIyEhSr1IRGXGDmRh91qwga9fWc/LJ8OCDmpczU+l62tbXO6efDz8cmpo0Hp6IR6hXqYi4Iz47QkNDQ8rQBk7r2dq1NwHbefDBmbkt0KPis1DE9TcLxaWXwu67w1NPwcqVg5+9QkTyk1rcRGTEDW6ezCLgx0ANkPm8nIFAgA0buuflLNTPsuG4+Waornaue3vuOafTgsbEE8lranETEXfEp7UqLi5Os+U03nqruk+QGGheTp/PR3V19QhVWrjOPRdmzIB//hPuucftakRkJCi4iUhWpZsY/dOfrmaPPfreM9C8nKWlpdTW1ibWaV7O1MaP755BIRx2BjkWEW9TcBOREZfJxOjgw5hS7rkndQAbaF7OpqYm/H5/Yr2u0+rf2Wc7p0pfeQVuv93takRkuHSNm4hkXSwWIxqN0tDQAMC4cQG2basmGKzlmmv8afbWdVnDdf/9cPrpsPfe8MYbPmCLXkuR/JT2GjcFNxHJGSeAHQU8zi67OK1AmcxJquA2PJ2dcOih8OyzABcC/6PXUiQ/qXOCiOSbHwDOcBWaSD43iorg8stjQAhweilo+jARb1JwE5EcqgSOYupUuOgit2spbL1nr/jBD8oxJgJsAjR9mIhXKbiJSE44Z+ac1rZgECZPdrWcgpc8DVh8+jBrNX2YiNfpGjcRyYlHHoETTwR4m1hsdyZOzHxfXeM2eIMbBLknvc4irtE1biLiPmu7xxODqwcV2kREpJuCm4hk3e9+B2vXAvwHuDnj/TKdl1NSy3T2ikAgkNhWRPKbgpuIZJW1MG/e610/RYGOjANYOBxOBIrkRcFtcDR9mEjhUHATkaz64x/hjTc+yrRp0N5+rQJYjmQ6e8WUKZo+TMRLchrcjDGVXcs1SeuqutYFB1onIt70wx86/158MfjTT5IgIyQ5FKeaPmzy5AAQBJoYO1bTh4l4Rc6CmzGmDDjWWtsIlBljyrrW0bWurb91uapRREbWX/7itLjtsgtccIHb1Yxufr+/x1AfbW0bOPjgev7zH7/mMBXxkJwFN2tts7V2YdePJdbaZuB0oK1r3Xqc0TlTrevBGLM63QIcktX/kMgoFO8s0Hvpr5Xmyiudfy+80Alvkj+Mge9/37m9eDFs3epuPSKSmZxf49Z1+nN+149TiA/j7ZjWzzoRyQPxzgJxA12r9vTT8NvfwsSJmiUhX51yCsyaBa+/Dnfe6XY1IpKJsbl+QmttxBizwhjz9DAe4+h023S1ulUM9TlEZHjirW01NTBNX7/yUlGR0+p2xhlw1VXwjW/A2JwfFURkMHJ6jVvS9WrrgXk4p0Sndq2bAmzsZ52IeMg//gEPPgg+H1xyidvVyECqqmD//eHVV2H5crerEZF0cnmqtJKegWw9sBwo6VpXAjT2s05EPOSqq5x/zz8f9tjD3VpkYGPGOHPHAlx9NXR2uluPiAwsl8FtGVBijKkCsNau7OqggDGmEmjr6sDQZ10OaxSRYfrXv5yWm3HjugOB5Lezz4a99oLnnoPf/MbtakRkILnsVdpmrV3WFdjmJ61fZq1ttNYuG2idiLgvFov1GKA1EAgQCoWIxWKJddGoM1vCf/837L23C0VKSgNNHzZhQvcp7cWLnd9f7/0y7UksItllCnVuOmPM6oqKiorVq1e7XYqIZ4XD4cQBOhaLUV5eTktLCx0dHYltfD4fpaWlNDU10d7uZ/p02LrV8tJLhv32c6duGbxYDPbdFzZtgieegM9+tuf98dBXqMcMkTxh0m2gKa9EpF/JA7ZGo9E+oQ2go6ODlpYWotEo//M/8fHAfqHQ5jF+P3znO87txYvdrUVE+qcWNxHpV/Kptcy8i9P3qBxrm7JQkWTTxo2wzz6weTM88wwcfHD3fWpxE8kJtbiJSC5NAf4EPOl2ITIE06bBvHnO7auvdrcWEUlNwU1EBhSfHaG4uHjA7YqKAgD86lefHXA7yW+XXOL0CL7/fnj5ZberEZHeFNxEJCM1NTX4fL6U940b56Ozs5qZM+HEE3NcmIyoj37UGR6ksxOWLHG7GhHpTcFNRPqVPPRHbW0tpaWlfcKbz+ejqKgUqGXBAmcapeT9xHuCQWcS+ttvh5dfTj8EjIjkjoKbiPQreawuv99PU1MTwaRRdQOBAF/5SpAtW5rYYw8/Z5/ddz/xlnA4zMyZzgT0W7fGOOKIciKRSOL+1tZWIpEI5eXlPcKbfuciuaFepSIyaMk9DI89FhobnWmuLrvM5cJk2IwxWGtZswYOPzwERICOPtv5fD6CwWBiyJj4fiIyLGl7lSq4icigxYNbc7OlrAwmToTXX4ddd3W5MBm2nkPA7AK8l/G+hXo8EckhDQciItkTv3h93jyFtsKUeWgTkdxQcBORIdqH5cthzBj47nfdrkVGUqZDwAQCgcS2IpIbCm4iMkQXs2MHnHGGM9q+FJ6amhrGjUs9BIzP56O6ujrHFYmIgpuIDMEU4DwAamvdrURGVu8hYPbbrxToOwRMaWkptUm/fA0BI5IbCm4ikrFwONx18Xo14Ad+zyGHGA0FUUB6DwHz5JNNVFYGgQBgCAQCBINBmpqa8Pv9KfcTkexRr1IRGZSODpg+Hd5+G1atgspKtyuSbHvvPZgy5X1gMmvWwJw5blckUrDUq1RERtb/+39OaDvkEDjmGLerkVzYZReAZYCmwRJxm4KbiGQsef7K2lpnWiQZLW4EtrFiBbzyitu1iIxeCm4ikrGHH4Z//cvpRTp3rtvVSG69AdxLZydcf73btYiMXgpuIpKxaNT59+KLYdw4d2sRNzjNrbfdBhs3ulyKyCil4CYiGfnrX+Evf3FmSDjvPLerkVzp7kkM8A/gUTZvhrlz/+BmWSKjloKbiGQk3tpWXQ1Jo0BIgQuHw4nZEay1rFp1PAAvvHAMHX3nnheRLFNwE5G0XnoJHnwQxo+H73zH7WrETccc4/QofvttuPtut6sRGX0U3EQkrWuvBWvha1+DPfZwuxpxkzHds2Vce63T01hEckfBTUQG9PbbcOedzu1LL3W3FskPc+c6PYv/9S+np7GI5I6Cm4gM6Mc/hi1b4Mtfhpkz3a5G8sG4cfDd7zq3NSCvSG4puIlIv2Ix+MlPnNuaTF6SnXeeM6PCn/8MTU1uVyMyeii4iUi/fvYzePddOOII+Mxn3K5G8smkSU4PY+jucSwi2afgJiIpbd/ePUK+WtsklQsvdHoa//KXsG6d29WIjA4KbiKS0sqV8OqrsN9+cNJJblcj+egjH4Gzz3Z6HF93ndvViIwOCm4i0oe13ae/FiyAMWPcrUfyV7yn8R13wDvvuFqKyKig4CYiffzxj9DcDLvt5ozdJtKfgw6CL34ROjq6O7KISPYouIlIH5GI8+93vgM+n7u1SP6LXwP54x/D5s3u1iJS6BTcRKSHZ5+FRx+FnXfu7jUoMpCjjoLDDoONG51TpiKSPQpuItJDfEDVb34Tpk1ztxbxhuRpsK67DnbscLcekUKm4CYiCa+/Dvfe63RGuOQSt6sRLzn1VJgxA1pa4MEHIRwOY4zps4TDYbdLFfG0tMHNGPM5Y8x0Y8whxpgFxpjp2S9LRNxw443O+G1z58L06W5XI16SHPajUQiFwlhrE/dba7HWKriJDFNGLW7W2leBW6y1S4DZWa1IRFzR1gZLlzq3NeCuDMU3vgFTp8KTTzpTYYnIyMskuBljzOeAP3T9bAfaWES8aelSZ27Sz30Oysrcrka8aOJE+Pa3nduaBkskOzIJbpuAY4GrjDFfAUqyW5KI5NqWLc5pUlBrmwzPBRc4Q8g8/DC8+KLb1YgUnkyC2xzAAKcBjcDfslqRiOTcz38Ob70Fn/gEHH+829WIl+22G3z9687ta691txaRQpRJcGux1i4C1lpr30OnSkUKSmdn9xAgtbXO0A4iw+F0Uohx++0hnO/9EAgECIVCxGIxN0sT8bxMgtvsrmvcZhhjDkGdE0QKym9+Ay+8AHvvDWec4XY14lXJvUX33DPGpEnldHZGiH/Xb21tJRKJUF5e3iO8qZepyOCkDW7W2ihwHHAGUNn1s4gUiPhF5N/9Lowb524t4l319fWJ29FolC1bWoCOHtt0dHTQ0tJCNKnnQvJ+IpKeSR5nJ6MdjJneNTxIXjPGrK6oqKhYvXq126WI5K2nnoJPfQomT3YG35082e2KxKvMMM6xD/Y4JFLA0r6Rxqbcy5jlwPlAKXAN8G7SAx4K7DdCBYqIi+KTyVdXK7SJiHhByuAGLLLWvm+MaQPmW2tfid9hjDk0N6WJSDatWwe/+AWMHw8XXeR2NVII4i1ngUCA1tbWfrcLBAJs2LABGF5LncholPIat3hQ6/r3K73u03AgIgVgyRKwFr72NfjIR9yuRgpJTU0NPp8v5X0+n4/q6uocVyRSODLpVbo++YeunqUi4mELFixh2bIOoJNbbz1AE4DLsIVCocTt2tpaSktL+4S38eN9lJaWUps0ynPyfiKSXtrOCcaY+3Gua2um6xo3a23eX+Omzgki/fve9+CqqwB+CZyqi8NlxMViMaLRKA0NDV1rApSWVvPMM7X4/X5XaxPJY2mvHcgkuH3FWvtA0s/HWGv/MNA++UDBTSS19nbYZx9nUnk4AmhScJOsca5hm8rOO29k82Z49llnhg4RSSltcMtkHLfk0HYIsGaYRYmIi2691QltRx4J0OR2OTIqbOLcc51b8Vk6RGRo0gY3Y8yp8dvW2meAyqxWJCJZs20bXHedczsYdLcWGV0uuQSKipx5cd94w+1qRLyr3+BmjPmKMeZm4HJjzHJjzP1d17sdlrvyRGQk3Xefc9CcOTPGmjXdF4VrHknJthkzoKoKtm+HH/3I7WpEvKvf4NZ1inQhsNBae7q19jRr7WnA4qE+mTFmXtdyTdK6KmNMpTEmONA6ERm6cDiMtfEBd2O0t5cTjUYS92seSRlp4XC4xxhtxhjuv9/53r90Kbz/vluViXjbgKdKrbXvAZuMMYu7lquBFUN5ImNMJdBorV0GlHQFs7Ku52kE2owxZanWDeX5RKRbfX09v/0tPPcc+P1RNm5soaND80hK9jhfFmyvZQ0VFU5oW7bM7b+GKi4AACAASURBVApFvCmTcdwqgWVdy1KgcYjPVUL39XHru34+HWhLWlfZzzoRGaYvfnE1ALFYtE9oi+vo6KChoSExrpvISIsP4XbDDbB1q7u1iHhRf1NeJVvba8qrVUN5oq6WtrgyYDkwG9iUtH4aMCXFuh6MMaszeEoNFCyScDhwNPAe8KG7pcio9vnPw0EHwQsvwL33wte/7nZFIt6SSYvbImPMGmPMo12dE4Z0qjSu69TnKmtt83AeR0QGw2nmWLhwF4qLiwfcMhAIJE5tiYy0oqLuVrdrroHOTnfrEfGaTFrcrkkecNcYc8wwn7PSWhu/KroNmNp1ewqwset2qnUJ1tqj0z1JV6tcxXAKFSkE69YBnJqYTH7ChBoikUjK06WaR1Jy4cwzoa4OXnwRHnoITj7Z7YpEvGPAFjdjzPR4aDPGzDDGfA5oGeqTGWPmxUNbV2eF5TjXutH1b2M/60RkiK69FqAoMZl8f/NI+nyaR1JyY/z47la3q64CNe6KZC5lcDPGrOsaeDdxrVnXdW5rGV6v0muMMS3GmHe7HrM56b42a21zqnVDeT4RgbfegjvuAGPg0kuddX6/n6amJoJJI/AGAgGCwSBNTU095pHUcCCSLd/8JgQCsGYNPPaY29WIeEfKuUqNMedba28xxuyC06tzBtBsrX0sfl+uCx0szVUq4rRqLFkCX/kKrFzZ9/54z1FdzyZuuPJKuOIKOOYYaNS5FREYxlylG8EZx61rIN5ia+1jyfeJSH7btAluusm5fdll7tYiksq3vw2TJsEf/gBPPeV2NSLe0F9wO8wYc0h8AVqTbmvKKxEP+NGP4IMP4PjjYfZst6sR6WvKFKipcW4vHvKcPCKjS3+nSl/GuZ4tVZPdodba/bJd2HDpVKmMZu3tsO++8O678PjjcNRRqbfTqVJx23/+A9Onw5Yt8PzzzhhvIqNY2lOl/Q0HMtda+7eUj2jMocMqSUSybulSJ7QdeWT/oU0kH+yxB5x7rnNa/+qr4a673K5IJL+lPFXaX2hLd5+IuK+jIz4ECFx+eeptUk0AboxRL1JxRW0tjBkDP/85vPqq29WI5LdMZk4QEQ+5/Xbn9NOhh8IJJ6TeJvUE4FbBTVwxYwZ89auwY4fTC1pE+qfgJlJAtm2DSNe8JJdf7ozfJuIFixY5/952G7z9tru1iOQzBTeRAnLvvc6ppgMOgFNOcbsakczNmgUnneSc6r/hBrerEclfCm4iBaKzs3tIhUWLnGuGRLwkPt7gT38KbW3u1iKSrxTcRArEgw/CP/8J++wDZ53ldjUig1deDp/7HLz/vjMOoYj0peAmUgCsdSbrBggGYdw4d+sRGaq6Ouff66+H995ztxaRfKTgJuJh8WE9ioq+wNq1AG9zwQU7qXeoeFZFhbO0tcGPf+x2NSL5J+XMCYVAMyfIaGEtFBU9BRzOkiVw6aVuVyQyPI895kw8P3Wq09lm0iS3KxLJmSFPMi8iHvHb3wIcDrzNt77lcjEiI+C//gs+8xnYtAl+8hO3qxHJLwpuIh5mLdTXx3+KMHGim9WIjAxjIBRybl97LcRi7tYjkk8U3EQ87Je/jPHUUyFgGnAdgUCAUChETEc68bjKSqeXaWurM4+piDgU3EQ8JLnTQXt7jLPPLgciwCYAWltbiUQilJeX9whv6qwgXpPc6haNwubN7tYjki8U3EQ8pL77vCjf+laUDz9sATp6bNPR0UFLSwvRaDTlfiJecfzxcNhh8M47cPPNblcjkh/Uq1TEQ0yPyUenAJkPL1+o73UpbI88AieeCLvvDq+8Ajvt5HZFIlmlXqUihek4QKOTSuGJj00YX0480QBP8/bbsGyZ29WJuE/BTSQP9D5YxZdU16Z1dlrKyx/F6ZDQv0AggLVWLW3iKeFwuMffrLWWX/1qDuDMxatr3WS0U3ATyQOpDlbW2pTBbdUqaGqCnXaqwefzpXw8n89HdXV1tsoVyakvfQnmzIG339ZsCiIKbiIeUlcX4vvfd24vWlRLaWlpn/Dm8/koLS2ltrY2sS4U754n4kHGwA9/6Ny+5hpnEnqR0UrBTcRDysrCPPUU7LYbXHqpn6amJoLBYOL+QCBAMBikqakJv9+fWK/hQMTrjjsOPvtZZzaF6693uxoR9yi4ieSBWCzWo1Us1UC6O3bAFVc4t6+4AiZOBL/f32Oojw0bNlBfX98jtIl4SX/vhQ8+iCVa3a67DjZudKlAEZcpuIm4ILkFLBaLUV5eTiQSSaxLNZDufffBc8/BPvvAvHm5rlgkOwbzXigri3Hccc6p0i9/+c99HifTDj4iXqZx3ERcYIxJdEYIhUJEIhE6Ojr6bOfz+QgGg1xxRT0zZ8L69XDbbXDuuX0fDzRWm3jPYN8LX/xiPZ/6FMBm3nprZ/bYo+/jgd4L4llpx3FTcBNxQc+BdDMxD1gK/Itt2w5g7FhnbTgcTjkrQigUUkuDeMLg3wsAvwRO5sIL4cYbUz9eoR7bpOApuCm4ST4a/MHqdWBv4DSsvT8LFYm4Y2jB7RPAM4wfX8TLL8NHP9r38Qr12CYFTzMniOSr+FhtxcXFA243cWIA2JtDDgFYmYvSRHIq0/dC96DSzwLL2boVGhqc+zLp4CNSCBTcRFxWUzPwQLqdnc5AuldeCaBWBClc6d4LPQeVDjFmDPzsZ/D005l18AENjSPep+Am4oLkloHa2v4H0p00qZQPP6zlM5+Bz39eA+lK4cn0vdB7UGlYx/nnQ2cnnHVWlJaWlj6dGjo6OmhpaSEajSbWpbomVMRLdI2bSB6IxWJEo1Eaus77BAIBzjmnmptvrmXzZj+PPw5HHeVykSI5kOq9UF1dTW1tbY/xCZ1r2XYHXgZmAK0ZP0ehHvekIKhzgoKbeEnyhdXV1XDzzc48jQ895HJhIjmWrpNBd6eGENDAYC4jKNTjnhQEdU4Q8aJ//hNuuQWKiuDqq92uRiQ/WWt5//0wxkwbcLvuTg0KbOJ9Cm4ieWjhQmeKq/POg4MOcrsakdyJz4AQl24GhEmT4IQTaoBMOzWIeJuCm0ge6Hmw+mzXqdEYkyYtcbEqkdwLh8OJ1rHkpXdwS+7UcM89tYwfX0rv8JaqU4M6+IjX6Ro3kTxiLZSXw1NPQTgMOsaIZOb//b8Y55wTBX4CbOy3U4NInlPnBAU38ZLly+GMM2CPPWDdOtDxRiQz1sKnPw1NTQD1WKtvPeJJ6pwg4hVbtsBllzm36+sV2kQGwxjoHq6tltdfd7MakexRcBPJEzfeCK+8AgceCOee63Y1It4SDof57GcNsALYmX32uWfATg0iXqVTpSJ54K23YP/9IRaD3/0Ojj/e7YpEvOnVV2HmTKcF+y9/cU6finiITpWKeMHllzuh7UtfUmgTGY7p02HBAuf2RRc5U2KJFBIFNxGXPfUU3HEHjBsH113ndjUi3rdoEey5Jzz9NNx9t9vViIwsBTcRF3V2woUXOrcvvhg+9jF36xEpBH5/94wjixZBe7u79YiMJAU3ERfdcw88+aQz/McVV7hdjUjhOOssOPxw+M9/YPFit6sRGTkKbiIuicWcqa3AObBMmuRuPSKFpKjI6akNcO210NLibj0iI0XBTcQlV17p9CY97DD42tfcrkak8JSXO++trVvhggucQXpFvE7BTcQFL7wAS5Y4g4b+6EdO64CIjLxIBHbZxRlm55e/dLsakeHT4UIkx6yF6mrYvh3OP99pFRCR7Nh9d7jqKuf2RRc5lyiIeJmCm0iO3XUXPPEEBAK6aFokF+bPhzlz4I03nOnkRLxMwU0khzZt6h4cdMkSmDrV3XpERoMxY+Cmm5xLE66/Hp57zu2KRIZOwU0khxYtgtZWqKiAc85xuxqR0WPOHOcShR07nH81o4J4lYKbSI787//CLbc4MyTEv/2LSO5ceSXsthv8+c9w551uVyMyNDkPbsaYsl4/VxljKo0xwYHWiXjZtm3Ot3yA2lo48EB36xEZjaZMccZ0A+eShbffdrcekaHIaXAzxlQCtyT9XAZgrW0E2owxZanW5bJGkWyIRODZZ2HGDPje99yuRmT0OussOP5453rTCy5wuxqRwctpcOsKY5uSVp0OtHXdXg9U9rNOxLNeeAEaGpzbt9wCO+/sbj0io5kxsHSpM5/pypXwi1+4XZHI4Lh9jdsUega5af2sE/GEcDiMMSZpGcOsWX9l61Y47zw45hi3KxSRffftnoS+pgbefdfdekQGw+3gNiTGmNXpFuAQt+uU0SccDmOT5tW57rodwBHstZcz/IeI5IfqajjySOc6t0sucbsakcy5HdzagPhIVlOAjf2sE/GgksT1bDff7Ey7IyL5oagIbr0VJkyAO+6ARx91uyKRzIx1+fmXA3O6bpcAjV23U61LsNYene6Bu1rdKoZdocggxGIxotFo10+v8OGHAT7+8RqOProW8LtZmoj0csABzkwKixbBvHnwj3/A5MluVyUysFz3Kq0C5nT9i7W2uWt9JdBmrW1OtS6XNYoMRjgcTtyOxWKUl5cTiUS61liglXXrIpSXlxNLmiQxeT8Rcc+ll8Ls2fB//wcXX+x2NSLpmeTrcQqJMWZ1RUVFxerVq90uRQqYMSZxTVsoFCISidDR0dFnO5/PRzAYpL5rosTk/UTEXS+8AGVlsGUL/PKXcPLJblcko1jaodndvsZNxPPiPUgbGhpShjaAjo4OGhoaEtuKSH4Ih8PMmmXYsuUiAE455R2M2V2t4pK31OImMgzDCWGF+t4T8SJjioDfA5WceCI89JCmpRNXqMVNJNustVhr2WWX4gG3CwQCiW1FJD/EYjFCoRDONanHYkyAX/86xI9/HEu3q4grFNxERkBbG0AN4Et5v8/nozo+WamIuGbgDkVgbSsQ4aKLymluVociyT8KbiLDEAqFsNaZFeG992rx+Urx+XqGN5/PR2lpKbW1tT32E5Hci3cQAohGo7S0tKS4NrUDa1v4/OejfPhh3/1E3KRr3ESG6ac/hW9/GyZNgj//OcYDD0Rp6JqcNBAIUF1dTW1tLX6/xnETcdvgrksNAGHg24CuS5Wc0DVuIkPRd85RZ+l9uuSZZ7rHfrr1VvjkJ/09vplv2LCB+vp6hTYRT2rFuQSiyu1CRBIU3ERS6D3naLxTQXJwa2+H006DrVth/nzndjzwxfUX+ETEPfH3c3HxwB2K/H7n/smTVwAzclCZSHoKbiIpdPc0cwQCAUKhUGL2A2udsLZuHXzyk3D99c528cDXe1FwE8k/NTU1fa5JjfP5fFxySTWnngrvvw9wP/0M0yiSUwpuIqTvadba2kok0j111Q03wL33wrhxW1m+HHbayYWiRWTQkr+Q1dbWUlo6cIei226D6dMB5lBT43xpE3GTgpsImfU06+jooKWlhZqaKPEOotu2fZWZM3NZqYgMR/KXNL/fT1NTE8FgMLEuEAgQDAZpamrC7/czZYozDdZOO8Htt8PNN7tQtEgS9SoVYbA9zabhXLS8GLhcPc1ECkD8M6C/9/PPfw5nnQXjxsEf/wif+Uwuq5NRRL1KRUbeu8CjwBVuFyIiw5Rph6Izz3R6kG/bBlVV8OabOS5UpIta3ERwPqzj74VAIEBra2u/2xYVBXjnnQ1MndpzPxEpbNu3w7HHwurV8KlPOS1vur5VRpha3EQGa6CeZuDj/POrmTo1pyWJSB4YOxaWL4d99oEnn4T//m/o7HS7KhltFNxEyKynGfjYZ59SlizR1FUio9Vuu8Ejj8CECR3cfz+MGXPlgIN0i4w0nSoVSSEWixGNRmlouArYARRz/PHVrFypqatEBH7/e/jCF2DHDoBvYO3tbpckhUGnSkWGYsmSJTQ03AO8BXQCdTz6aANLlixxuTIRyQef/nSME04I4cxneidTpvQcpFskWxTcRFKYPz9MScnLQDFf+AJs23aBZkAQGcVSDdL9hz9EcIYGsrz3XitXX909SHeq/URGgoKbSC/vvgvHHQfr18OcOXDffc5FySIyemUySPfWrR28/HIL0Wg05X4iI0HXuIkk+eADp7v/X/8KM2fCn/4EaeahFpFRYHCDdBfhXGLhKNTjrGSFrnETydTWrfCVrzih7aMfdS4+VmgTkcGzODOsiIw8BTcpaPFR0Xsvva872boVzjgDHn0UAgFYtcoJbyIicdZarLUUp/1GV8yhh7aycaNa2mTkKbhJQVuwYAF1dXWJn4uLi6mrq2PBggWJdVu3wumnOxNJT5kCv/sdHHCAG9WKiBcMNEj3hAk+dt21mr/9DY45BkCjdcvIUnCTgpKq51ckEkmsa21tJRLp7vm1dSucdho8+KAT2hoboazMhcJFJK9lMki3z+fjYx8rpamplv32g2eegT32eI6NG3NdrRQyBTcpKJn0/Oro6KClpYWrr45y2mnwq18BbOIPf4DZs3Nbr4h4Q/KXQr/fT1NTE8FgMLEuEAgQDAZpampi//39/PGPsN9+8J//fIRjjoEBpj8WGRT1KpWCMrieX7sAbcAmoBJrm7NTlIgUnHA4nHKoj1AolAh5b74J//Vf8NJLcOCBzjW0unZW0kh7EFNwk4Iy+C77bwInAM+oy76IjLg333TGhXz++e7e6jNnul2V5DENByKjT6Y9v4qKprFu3e5Y+7ccVSYio82ee8ITT8CnPw2vvw5HHglPPeV2VeJlCm5SsAbq+WWMj4suquZjH8txUSIy6kyd6gwx9MUvwsaN8KlPxTDm8wMOUSTSHwU3KSiZ9PwqKvJxwAGlNDTUptxPRGSk7byzM+TQV78aA6LA00DqIYpEBqLgJgWld8+vv/61iSOPDAIBoIgJEwIsWhRkzZom/H5/yv1EREZK8mfLli0xnn22nLFj45PT9x2iKNV+IskU3KRgdXTAUUeto7GxHtgALGLLlne46qoGlixZ4nZ5IjIKpBqiaPv2vkMUaXJ6yZR6lUpBWr8e5s6F5mbnFMVddznzkIqI5NLgerqPBbYnfirU47MMSL1KZfR58EFn9oPmZigpgT//WaFNRLygE/iu20VInlNwk4KxZQtceimccgq89x6cfDKsXQuHHup2ZSIymmU+Of004HpOOcUCu+agMvEiBTcpCP/4Bxx+OFx3HYwdC9deC7/4hTP/qIhIPhhoiCKfz8fcudVMnuz0PoV/8Pvf57Q88QgFN/G0HTsgGoU5c+DZZ51To088AZdcAoO6tEREJAsynZy+tLSUn/2sluZmOOIIgL04/ni44ALYvBkuv/zyHuO+xZfLL788t/8hcZ2Cm3jWP/8JRx8NwSBs3Qrz58Pf/x7/0BMRcV/vIYpOOukkOjr69io96aST8Pv9lJY6Xz4/97k/MHYs/OQncPDBMe6996GUge+hhx7SMCKjjHqViidcfvnlLF68OGnNzsDFwCL22MPPbbfBF77gUnEiIiPMGENzs+Wcc+D550NABOjos53P5yMYDCaGDzHGqDeqt6lXqXhT8rfGWCzGQw/1/ra5GbiWXXct58knY4nQpm+bIlIoysoMzz8/AVhCqtAGTmtdQ0ND4tSpFD4FN8lLvQetfPnllj6nF6CDDz9s4bbbNGiliBSqrThfVEUcCm6St5xvkNNoaLiWLVv0bVNERpdMhxGZOjWQ2LY3dWooPApuklOZf4jsCjQAr6BvmyIymg00jAj4ePfdaubNg1df7XlP6stM1KnB6xTcJKvSX6vW80Nk0yb4/vcBXgW+D0xm3LhpAz5HIND/t00RES/qPYzIpEmTUm7n8+0F1HLLLbDffgBLEwEuPjdqql6sLS2aG9WrFNwkq1JNsJzqQ+Tll1uorIwyfTr88IcAkznuOPjLX+CyywYetLK6ujp7/wERERf0HkZk/fr11NXVEQgEKCoqIhAIUFdXxzvvPMM//+nnnHOgsxNgHjNmbMOYu2houDLFtcEOXWbiXRoORLJqcB8IAWAD8HsgjLX/CzgtdSUlJbzzzjt99wgEWL9+PX6/H3A+7NTkLyKj0UsvwQEH3AWcBYzBaZvJ/BhfqHnAYzQciHhJK1AGHA/8NbF2oG+byaENdJ2GiIxe++8P8HXWrx/DxReDM/dp/3SZiTcpuEnWWWvZts0yefLAPaMCgWKsbU75IeL3+6mvr2fDhg3s2LGDDRs2UF9f3yO0iYgIzJjhzNu8cGENY8emvsxk3Dgf55+vy0y8SMFNsmb7doBjmD8fPvIReP/9GkDXqomIZEtyp4Yrrqhl111TdWrwsW1bKTfeWMuZZ8KDD8L3vveD3BUpw6LgJkMSi8UIhUI9Tl2GQiHeeivGr34F8+Y5YQ0aWbYMWluhtLSWnXfeg1Sn8CdNmkRtbW3i5+QPHxERyUy6Tg3TpgU44YQgZWVNfPCBn3vvhVNOgRtvvIIzzoC77oING9yrX9JT5wTJSPJF//13FhgP7Ac0Ac4pzGnTWvnWt4qZOxc++Un44IMY0WiUm266iY0bNzJt2jSqq6upra3VaU8RkRx65RVYsQLuvx/Wru1ebwzMmePM/3zccc7t8eMze8xYzPmM/+lPf5r4jK+pqdFnfObS9+iLX5hYaAuwuqKiwsrA2tvbbV1dnS0uLrbGGFtcXGzr6upse3t7j+2cPxVrt2yx9txz6+zYsT6L012p1+Kze+9dZxsarP3737v3ExGR/NXSYu2NN1p7wgnWTphgLbRbqLNQbMHYceOKbUVFnf3Nb9rt5s3d+4VCocTt9vZ2O2vWLOvz9Tw++Hw+O2vWrB7HleT9pIf0+SaTjby4jLbglmkAG8yb7N132+2zz1p7xx3Wwg0WnrDwQdcbOVVoS72IiEh+Sz42vP12u91nn1l2zJjeX9B9FmbZoqJ2e/DB1n7zm9bCfLtmjbUdHdbW1dXZsWPHpjwOjB071tbV1SWeo/exIdNj2Cig4OaWbP8RDvVbTvKbZaA3GUywY8bUWbApFqPgJiJSQAZzbHBa4nofF7ZYmDKkY4MbLXV5HBQV3JJl8xc1Us3Fgz11aW1m33I2bbK2udlaONnCd7ta0CaneWMFLKyzsNzCIvu731n7zjvWFhcP3OIWCARS1ikiIvlpMIHLWT5t4UILd1l4wcKODL7UGwsLLcy1MMe+8461nZ25aanz0CldbwY3oAqoBILDeIzV++67b+KVyHaY6h2kej9P8vP190c4uJYz0/XtZrqFMWneLMUpvh3FH2Pw346G8yYTEZH8M/jg1vPY0N5u7ZQp6S6jCaQ4DnVYmJpmvzFdx7ExwzhmDv8YHX/OoTQADWI/7wU3nKHzq7puzwPKhvg4q0cibGT6hzG4P/QiCwdZOMzC0fbhh6297z5rv/SlOltUND7lPsZMsHvsUWdnzLB2yhRrnW83mQawIgvvW3jWwkP2ggusXbLE2kmThtZyluffVkREZJCSP+OHelZloOPsmDHj7Akn1NmLL7b25JOthb9ZeHcQx7D48W6TLSmx9rDDrC0pqbPGTEi5T1HROFtVVWd/+1trn3jCWphj4RMW9ut6vME3Wgzm2Jd8xm+QATNsPRjcrgEqu24PudXNCW4zLEQsXGthYppfziQL91i4z8L99tRTrf3yl63db786W1SUOpkb47O7715nDzrIWnjJwqsW3hzkH2HyMrhvK5MnW7vvvtaOGTPwftOmBWxn5+DeZJk2TwcCAVtUVGQDgUC+XB8gIiKDNBLHhqG0gH3wgbW77jrwMayoKGCnTrXWmOEdM7uXTE7prrPwooV/2E9+0tqyMmv33LPOGpM6DxQV+exBB9XZM8+09mtfsxZ+Zs8919rzzrO2rKwuRUePeKD12SOOqLMLF1p72WWJ18VzwW0pXa1sXcHtmhTbrM5gaYOKQfyiRiZMZbJPUVHAHnigtXPmWAurLTxi4f40zxNfSq3TrDymx5ss26dm1XImIlK4Rvo67XRf6nsHxUyOYdu3267j38csfCqD46Wx8DsLf7LwlP34x6392MesLSoaauDLxX6jPrh91EYizmnBiRMHfuEmTQrYe+6x9t57rXUunDzVOhfyZxKmZlnYz776qrVvvmntggVDC1JDbZ7OxZtMRERGj2wfG0ZiZIThnNLt7xg9YYLPXnhhnV23ztoXX7Rdx/eDLczOMA+caeHrFv7bwrkWzstwv4UWLvNscEs+VVrFsE6VDj7RD/UPY6hBaiRqjD+nApiIiHhRNlvqeu+X66A4hP08F9zKgHldt4MMo3PCUHuVDuUPY6hBykNdlEVERFzlhVO6I7Cft4KbUzfzcE6TzhvGY/Q7jls2wtRIDwqoljMREZGB5esp3Wz2KtUk8xmIT5qridFFRERGr6HmgUHsl3aSeQU3ERERkfyQNrgV5aIKERERERk+BTcRERERj1BwExEREfEIBTcRERERj1BwExEREfEIBTcRERERjyjk4UDe2GWXXfY65JBD3C5FREREJK3HH3/8RmvtdwfappCDWzswHvir27XkkXiKfcbVKvKPXpfU9LqkptelL70mqel1SU2vS2qHADFr7d4DbTQ2R8W4YS2AtfZol+vIG8aY1aDXpDe9LqnpdUlNr0tfek1S0+uSml6X1OKvSzq6xk1ERETEIxTcRERERDxCwU1ERETEIxTcRERERDxCwU1ERETEIxTcRERERDxCwU1ERETEIxTcRERERDyiYGdOEBERESk0anETERER8QgFNxERERGPUHATERER8QgFNxERERGPUHATERER8QgFNxERERGPUHATERER8QgFNxERERGPUHATERER8QgFNxERERGPUHATERER8QgFNxERERGPUHATERER8QgFNxERERGPUHATERER8QgFNxERERGPUHATERER8QgFNxERERGPUHATERER8QgFNxERERGPUHATERER8QgFNxERERGPUHATERER8QgFNxERERGPUHATERER8QgFNxERERGPUHATERER8QgFNxERERGPKNjgZoy5wRhzg9t1iIiIiIyUsW4XkEWHVFRU3aTBewAAEvtJREFUVAAXuV2IiIiISAZMug0KtsVNREREpNAouImIiIh4xKgKbuFwGGNMnyUcDrtdmoiIiEhaxlrrdg1ZYYxZXVFRUbF69epU9wFQqP93ERER8SRd4yYiIiJSKBTcRERERDxCwU1ERETEI/JyHDdjTGXXzWOttQtdLUZEREQkT+Rdi5sxpgwnsDUCZV0/i4iIiIx6edfiZq1tBpq7fizp+llERERk1Mu74BZnjAkC8/u5b3UGD3EIjz8Opm/P2sQgICnuExEREXFFBsOU5d2p0jhrbQSYb4yZ4nYtIiIiIvkg74KbMSb5urb1wLze21hrj063AM9QUeGk116LoWuEuxT3jdQSDoUSz5O8hEOhrD6vFi1atGjRomXgJW+P0ZnkpHybPaDrFGmztbbRGLMUWGWtXTmEx8mLmRM0S4OIiEh+ysNjtCdnTlgGlBhjqgCGEtpEREREClHedU6w1rbhhDcAhTYRERGRLvnY4iYy6oTDYYwxfZZwOOx2aSIikkfy7hq3kaJr3MSL9PciIpI7efiZ68lr3EREREQkBQU3EREREY9QcBMRERHxCAU3kRTUWSA1vS4i3qL3bOFR54Ts15Gz55K+wuEw9fX1fdaHQqGMPrhy/fvzyt+LV+oUEYfes6nl4euizgmjhb5VpRYOh3u8Ia21WGtH/esiUmj0Geht+v1lTi1u2a8jZ8/lxvN5xVBfF/3+UvNKnTL66G8zNa+8LvrMVYubiCfEYjFCoVDi50AgQCgUIhaLuViViMjgqOUs+xTcJKf0pnYk/39jsRjl5eVEIpHEutbWViKRCOXl5T3Cm1dfJ/3eRUYHXZ6SfTpVmv06cvZcbjzfUHnldclWncaYxGOGQiEikQgdHR19tvP5fASDwUQHi+T93JRvr6dInP7GUhvtn7n58nwZ0KlSkXwVb3VqaGhIGdoAOjo6aGhoSGwrhcULLZFeqFFkNFGLW/bryNlzufF8Q+WV1yWbLW5D5ebvNhaLEY1GaWhoAKC4uJiamhpqa2vx+/1p9/fK32eueeF18UKN4J06c220f+bmy/NlQC1uIvkqfu1HcXHxgNsFAoHEtrk22q7FExHJdwpuIinkspdnTU0NPp8v5X0+n4/q6uoRf85MJQ9eHI1GaWlp6XNat6Ojg5aWFqLRaMr9RERk5Iyq4KYhF0afTH/nuW5ZSq6ptraW0tLSPuHN5/NRWlpKbW1tyv1yRdfipafrwNxX6J/vQ/0bK/TXZVSKn4IptAVYve+++9q49vZ2O2vWLOvz+SyQWHw+n501a5Ztb29PbBsKhexwtbe327q6usTzFBcX27q6uh7Pkw3x58t32aoz+Xc3mN95ci11dXV99knet66uLuV+w9H77yUQCOTk7yWdVK9BpstAj1eohvr/88LrMtgaQ6FQyr+Lkfh8jT9+XK4/392U7veQL69Lvr8X3DpGZyB9vslkIy8uwOpcHoy9/mbJtWzVOZzf+UgGlOHUn0+/v+RaiouLB3wdAoFAyv2szesPyRGV7wer4ci3/5vbX7bcku71zJfXJd/+XvLlGJ0BBbdcHYy9/mbJtWx+mOd6yUb9+WIk/q7z/ENyROXbwWok5dv/LV/es7mWSXDLh9clH/9e4vI86Cu45TK4efnNkmv5/GE+nJalkao/X3j11LNbBvv781JLZD4eiPPhMzfXFNyGV4/br0smpaZbXA9Y2Vrodao02wfjfPmjyPcPpmwfqEbid64W0/5lei1evrwfci3d/2EoQXgkrwlL9Xqne/zhvmezeSCOc/PLVq5lEtzivPglVMFNwS3xSvz/9u6fN47jjOP4b4QgYCFEhKOrkiIhyzQM4wDqGNjUK7CSvIGEhlylMGWoEU9sYtNVqiR+B/7zBuKooKsYUCKkj6UqlaAYAsxC3aS4XepI7t3uzszOn93vB1joqLvlPtybmee52T83dDIuvbMMJfZ5BSHec584fU/Izu39W6VL8ihkkAyqT1Lt2j5jjxGh+2yMwo1Z3YvP10r8EBqjvaTI0T1yw3CFm6Qd13VjLIp8VWnqzuL6qXjoK79i75dQScf3Ks/cBq3Q+iSPkmZDhi68QxSzMWKsheizQ7Xpgk42D6ptf6beL665KObRmMwLWr/CTdKXkj5tWD6T9J8uG0i1SDrd29u7sDdck3GXgTJUZ+k6KKf+VNw1TtdEFSpZlVaADZXkQoud/GPr+z50TToh+oJrjF3XC9VnfePsK9db6oTm2zZD7xfXXLScn0vKmRGOqngXbm+vee7nXTaQalFD4dZjx3mvN3TBEDoxjj0JjH292NriTP2p39dQf9/y7+w6E1lqn+VelsPIbWwJkYtS5LCMP9T7FW5rVyzgUGnKwm3obZVSELkkqsvrxYgz9XpDH7IOLfYg2VcJhxN9z3HLuc9O4V6Wrm1sjBeJxF5C/n0ZrudXuOn1odLPlpZiD5X22HHB1htycC2h0YdIcDHizGW9UpSyX3Luf75XlebcZzM/hyiL7bWtl/p0mK58+kIJOSzBet6FW5JDpZLuSNqXdM/jd4y+cKvlPJOVy+BTynp9pZqpy32/+G6vS/8LkXR8ZiJz67OhE7GrMbWxWs4XiYTIRTnnsATr+RVua1eUfuK6bsvv3ZV0p3p8IGnX8fdMpnAraSYr4/MKsljPVSnbG0ucOSQd377nerVf7vfuG1MbC7k/h9ovIXJRSTkswnrhCjdJf1xaPpT0uOu6fRZJH0narx47z7pp5IUbM1l5rjf2E7JLu09djMItVdLJbYxYfm5K97KceuGW+qrS2ohyQ9DC7UNJP62WtyW903XdPoukv6qaZasKt48aXnPaYXkp7VnJTmT5zkoPrDSz0rXq3wfV/3dZ76aVTPVvl/XGvhzZ5kHxqOF1y/vyZ1a6nMQ3qv//bsV6LPksXftDyve9a4xaevygIb7lOB+sWK/PEnt7Y1+W98v6Qngx3ueyP31zUd/1cl2OLv1tXccIzYMVbg3F01vlFW5Hl3baqmScchm6kHJtTJf3USkF39DvOcmq7GW5HYToD0MkHdcYYyf+UPvSZz+59HXX9YZeQo0tuf59oZZcc5H7+xdyxu1LSX/T65vwvj9Q4bZ8qPSOBjhU2pfP4a96nVVSHs5I8ZVQXeO8LNcv416fDNcvQ8SRu9hxxugPfbbnEqdPjCHbpm+fzf2GuEO8dz7rjf1eiK7rpT5NaOibbStw4bbyCtOQixYXJxxUj+9pgIsT2rg0DNe7K8dOHCEG8xhxljJohU6OvnF0Fftq1NjbizS4XjFEsgoRY4lX7cXWN87Y91WjEH79fM0nFw39jQue/da/cNNAV5Cu2d6BFodJDzx+h3Ph5tIwXAem2InDZ3ux46zl/JVJoZOjTxxjTKqZD65XDPE+hIgx19nEnLTFGerDZK430g1tqPYSus92iafm2o8ccoN74abFRQjfSvpG0mNJP2j7hbks8izcYjaKmInDpzHFjjPm++AqRKf2QRJ4/Xwt4uDaO06Xvy9EjCVetRdbjDaW8410QxuqvcTODSG259BevAq3D5ceb0r6XdsvzGVRQYVbLUbiSJ3g+hxCjvk+uIp9SDf1+R2pxBjMU8xKdekPoWJ0PdRW2i1gXI21jaXiGmefHBbjKEeI973nUYC59Szc3rn089tLj7OefVOgQ6W+XwTdZ1s5nzsWOs42sTtnKEOfh0ISWP98rME1ZiFcerFOG/Nvm33izIVrnF3eh1pJuajODbPZzF67dm1dbmivb9Y+ubjZ7s7S8r6kt6rlz102kGpRIee4pR6UuxYasWeWYnfO0IYaXEkCq5+vhRhc2/qDa5whlBBjyu25GrKNldhnXWdah56hnUAu8i7cvpH0Z0l/aVgG+eaEUItP4RbiqlJX9+/fb2wU9+/fb415yM5y2dAzS6VcVboKhVtYfZLqVE7ELyHGFNvrq88VyLUSz6N0lXuh6Joz+8RVi5SLvAu3lV8kv+65HBafwu2yHlOcxco9CUzlUvguv7dWYhJw1RZn6pnr2Lc7cYnRdz1Xue1L17bi+hVNuXyocFXKGBFLhHrAr3AreQlZuE3BWJNAKjEKtxKTgKu+cZZ2ONFF7n0214saQvWhVKeZMOaOHoUb1st1cM1le66GijOXe0nFFrtIKaGdjflvs3bYws11CRFnqqt7XZXSXkaEwg3DonBrFivOXC8uCSXVfepKaGdj/tusHW/h5rtebKXEOSIUbhgGn/7Wy62gXX4u56t0czmsVEI76xtjqj7rasjCrZbrTZdzUkqcI0LhhnHJfRDJtaANPcswlNAF5piTagkx+ohRuNHG2pUS54i01jfXBBRgPp/LGHP+szFGxhjN5/N0QTWYz+eNHS23OHNWv7fHx8d69epV42tevXql4+Pj89eGUko7g7ujo6Pzx4eHh9re3tbGxsaF12xsbGh7e1uHh4eN6wFJdanuSlzEjBsmRD0Oleb8TRTr4mpb1v2+MSntkKerWO9dqiuQc2+bpV3ANCLMuAG46L333rsyw1Db2NjQ3bt3I0d0UT043bx5c+3rZrPZ8ge1yRj7rO7Z2dmF2a3ZbKajoyOdnZ0Nsr3r16/r4cOH5z8/f/5cDx8+1PXr19eu5zo7m+us7vL2z87OdOvWLZ2cnJz/34sXL3RycqJbt25deC9Sxz1JXaq7Ehcx44YJUcun91KuKl3+G3zOP5rKrNQY5NI22/rQ2IXqe/DWXt90eVGJC4UbpqRv0sn1myhySeKIJ5eCgcKtjAuYJqC1vjHWjvMwgzHmdG9vb+/09DR1KMBg5vP5hcM8taOjo06HMOpDNrmOA2dnZ/r44491fHwsaXHY7O7duzo8PGw9lIUy+FxcErLd5t4XhpbL+wC1vhEUbsCElZKsSokT/eVSMEy9jRljzv/22WymFy9erHztbDbT8+fPr6yHIFo7BBcnAACSqg8BcUFKHnK/gGnqKNwAAFlIUTDkepVnbNzfrhwcKgUmrJTDQ6XEif7m8/l5kVTfhuLp06cXbr5cFwxff/31+bmNy+shPM4vTYZz3CjcgNVKKYhKiRP+KBjyQt+LjsKNwg1YrZRBuZQ4EQ7veR54H6Lj4gQAV3FeDwCUicINmKBSvjaJAhNIg76XLw6VAgCywyE6TBSHSgEAAMaCwg0AAKAQFG4AgGxwbhWwHue4AQAA5KH1HLfvxYiiL2PMfvXwtrX2g6TBAAAAZCK7Q6XGmF0tCrZHknarnwEAACYv60Olxpin1trthv8/7bD6zp50o8sLAQAAkrO23NuBGGPuSXo3dRwAAAC5yH3G7XNJv7fWvnRYl4sTAABASfK8OMEYc9Dw38+stY/qc9qstU8kPZN0IOkkZnwAAAA5SlK4WWs/WfP0vqQn1eNNSY+HjwgAACB/OZ7j9omkLWPMHUmy1n6ROB4AAIAsZHcft+p8tnpGjqINAACgkuOMGwAAABpQuAEAABSCwg0AAKAQFG4AAACFoHADAAAoRNbfnODDGPPfGzdu/GhnZyd1KAAAAK2++uqrP1lr/7DuNWMu3L6T9H1J/0gdS0bqKvbfSaPID/ulGfulGfvlKvZJM/ZLM/ZLsx1JZ9baH697UXb3cQvoX5Jkrf1V4jiyYYw5ldgnl7FfmrFfmrFfrmKfNGO/NGO/NKv3SxvOcQMAACgEhRsAAEAhKNwAAAAKQeEGAABQCAo3AACAQlC4AQAAFILCDQAAoBAUbgAAAIUY7TcnAAAAjA0zbgAAAIWgcAMAACjEpAo3Y8y91DEApTLG7KaOAfmifQD+utQpY/6S+QuMMfuSfpk6jpxU+0SSbltrP0gaTEaMMQfVw232y0LVVj6S9IvUsaRkjLkj6aWkXWvtSep4ckH7aMZYchV5Z7WudcqkZtzwWvXp+La19pGkXT4tL1Qd55G19hNJW0uDzKRV7eTb1HGkVPeRal+8pM+8Rvu4irHkKvJOGJMo3Iwxu1VDQcVa+2Tp086WtfZJ0oDysSWpHmCfVT8DkvRbLWbbpEXbmHwixlqMJZeQd1brU6dM5VDpG6kDyFV1PP3d1HHkovp0XNuV9GmqWJCdTV2cVfphqkCQP8aS1cg7jTrXKaMo3JbOI1j2zFr7iNm29ay1J8aYz40x/7TWvmxfYxqqKfy/T+UT4bo+FD0YYESmNpZ0Qd65qG+dMorC7dInm8u2jDFbS493p9KB2gpaaTF1rcU0/oGkSZxs3bFI2Z/SyectfQgLL/X6U/GmpP8ljAXlmNRYss6U806LXnXKKAq3day1X0jnyXozcThRtSTjfUl1w9iU9Hj4iPLQVqQYYw7qgdYYs8+s0/nVlG8aY+7UfWqCPpX0ZvV4S9Lk20WN9tGMseSKyeaddfrWKXzl1UQZYzYl/UaLc3ZuW2s530DnV4J9rsV+eUPSrxlsUasG1mdanFjNLCVWYiy5irwTBoUbAABAISZxOxAAAIAxoHADAAAoBIUbAABAISjcAAAACkHhBgAAUAgKNwAAgEJQuAEAABTi/5d9T/l4ido8AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Create the samples\n", "size=5e3\n", "\n", "sample = norm.rvs(size=int(size))\n", "\n", "values, edges, ex, ey = hep_spt.errorbar_hist(sample, 40, range=(-4, 4), uncert='freq')\n", "centers = (edges[1:] + edges[:-1])/2.\n", "\n", "# Extract the PDF values in each center, and make the pull\n", "ref = norm.pdf(centers)\n", "ref *= size/ref.sum()\n", "\n", "pull, perr = hep_spt.pull(values, ey, ref)\n", "\n", "# Make the reference to plot (with more points than just the centers of the bins)\n", "rct, step = np.linspace(-4., 4., 1000, retstep=True)\n", "pref = norm.pdf(rct)\n", "pref = size*pref/pref.sum()*(edges[1] - edges[0])/step\n", "\n", "fig, (ax0, ax1) = plt.subplots(2, 1, sharex=True, gridspec_kw = {'height_ratios':[3, 1]}, figsize=(10, 8))\n", "\n", "# Draw the histogram and the reference\n", "ax0.errorbar(centers, values, ey, ex, color='k', ls='none', label='data')\n", "ax0.plot(rct, pref, color='blue', marker='', label='reference')\n", "ax0.set_xlim(-4., 4.)\n", "ax0.set_ylabel('Entries')\n", "ax0.legend(prop={'size': 15})\n", "\n", "# Draw the pull and lines for -3, 0 and +3 standard deviations\n", "add_pull_line = lambda v, c: ax1.plot([-4., 4.], [v, v], color=c, marker='')\n", "\n", "add_pull_line(0, 'blue')\n", "add_pull_line(-3, 'red')\n", "add_pull_line(+3, 'red')\n", "\n", "ax1.errorbar(centers, pull, perr, ex, color='k', ls='none')\n", "ax1.set_ylim(-4, 4)\n", "ax1.set_yticks([-3, 0, 3])\n", "ax1.set_ylabel('Pull');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting efficiencies\n", "Let's suppose we build two histograms from the same sample, one of them after having applied some requirements. The first histogram will follow a gaussian distribution with center at 0 and standard deviation equal to 2, with 1000 entries. The second, with only 100 entries, will have the same center but the standard deviation will be 0.5. The efficiency plot would be calculated as follows:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD+CAYAAADBCEVaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAASR0lEQVR4nO3dMWgbaZ/H8d8/mBjiYv3aXHUQ79nEaYOS47IsWIFYVaqAX9K4jmGbnFPk3WabVNkt4ndTBJS3fZtwga0WltgG2xDuDu71vcU2yca65LhiG/tVY4PeiPyv8Mir2JJmxh5ppEffDwRr9Ew8f49GPz9+ZvSMubsAAGE5l3cBAIDsEe4AECDCHQACRLgDQIAIdwAIUO7hbmZ/NLM/5l0HAIRkJO8CJF0pFotFSffyLgQABoy1a8i95w4AyB7hDgABItwBIECEOwAEiHAHgAAR7gAQIMIdAAJEuANAgAh3AAhQP3xCFQjKyuobfb/+y4nn7928pOXSbA4VYRgR7kDGlkuzWi7N6vOvf5QkvXt0K+eKMIwYlgGAABHuABAgwh0AAkS4A0CACHcACBDhDgABItwBIECEOwAEiHAHgAAR7gAQIMIdAAJEuANAgAh3AAgQ4Q4AASLcASBAhDsABIhwB4AAxYa7mS2Y2byZPWjT/m309W7WxQGDaL9W1+PV10fLhYcv9Xj1tfZr9RyrwrDpGO5mVpAkd1+TVG0sH3PXzHYkVbpQH5C7ldU3idfdr9V1++krlTd/ezvsHXxQebOi209fpQr4NNsFjovrud+RVI0eVyTNt1jn9+4+E/0CAILT6mbX7ZS3dvR+90C1+sdPnq/VP+r97oHKWztd2S5wXNwNsscl7TUtT7ZYp2BmklRw9++aG8xsI0ENVxKsA+SqcbPrs6jVP+rJ+ls9WX+bQUVAZ2c+oeru30W99kkza9WzBwD0WFzPvSppIno8Lmm3udHMFiTJ3V9EbdPN7e5+I66AqHdfTFQtkJN3j24lWq/w8KX2Dj60bZ8YO6/tb0qJvlcWfy1geMWF+3NJ16LH05LWJMnMxt29qsNx+MaZoxlJ5W4UCZzGyuqbluPW925e0nJptivbXPxiSuXNyokxd0kaHTmnxesXu7Jd4LiO4e7u22Z2LRpuqbr7dtS0Lulq1H7XzPYk7TS1A7lbLs1quTR71ANO2vs+7t7NS4nXXZqb0U8//3ripOroyDlNTV7Q0txMV7YLHBfXc5e7P2vx3NVO7UBI0vTyx0ZH9MNXX6q8tXN04nRi7LwWr1/U0tyMxkZj33Kn2i5wXPIjDUAiY6Mjul+6fBTuScfYgSwx/QAABIhwB4AAEe4AECDCHQACRLgDQIAIdwAIEOEOAAEi3AEgQIQ7AASIcAeAABHuABAgwh0AAkS4A0CACHcACBDhDgABItwRrP1aXY9XXx8tFx6+1OPV19qv1XOsCugNwh0DYWX1Tar192t13X76SuXNytFzewcfVN6s6PbTV6kCPu22gX5AuGMgtLrRdSflrZ0T9zGVpFr9o97vHqi8tdO1bQP9gNvsYWA0bnR9VrX6Rz1Zf3t0G7ysray++eQXQqPuezcvcV9U9AzhDmRsuTRLiCN3hDsGxrtHtxKvW3j4UnsHH9q2T4ydT3zj6qz+YgB6iTF3BGnxiymNjrQ+vEdHzmnx+sUeVwT0FuGOgXDv5qVU6y/NzWhq8sKJgB8dOaepyQtampvp2raBfmDunm8BZhvFYrG4sbGRax0Iz36trvLWztGJ04mx81q8flFLczMaG2VEEkGwdg2xPXczWzCzeTN7ELNex3ag18ZGR3S/dPloefubku6XLhPsGAodw93MCpLk7muSqo3lFuvNS/rn7MsDAJxGXM/9jqRq9Lgiab675QAAshAX7uOS9pqWJ4+vYGaFqGcPAOgTWQw+TrRrMLONBP//SgY1AACaxPXcq/otvMcl7TY30msHgP4U13N/Lula9Hha0pokmdm4u1clTZvZdKM9Cvvtxn929xtxBUS9+2LKugEAHXTsuTeCOroaptoU3OtR+wt3f6HD3v14NwsFACQXO+bu7s9aPHe1xTon1gMA5IPpBwAgQIQ7AASIcAeAABHuABAgwh0AAkS4A0CACHcACBDhDgABItwBIEDckgbBWll9o+/Xfzla/vzrHyUd3hN1uTSbV1lATxDuCNZyaZYQx9BiWAYAAkS4A0CACHcACBDhDgABItwBIECEOwAEiHAHgAAR7gAQIMIdAAJEuANAgAh3AAgQ4Q4AASLcASBAhDsABIhwB4AAEe4AEKDYm3WY2YKkqqSCu3/Xon0+elhy9z9kXB8A4BQ69tzNrCBJ7r4mqdpYPtZeitoLx9sBAPmIG5a5o8NeuyRVJM03N7r7dlNvfdrdtzOuDwBwCnHDMuOS9pqWJ1utZGYPJC21eH4jQQ1XEqwDAEghkxOq0Vj8kpmNZ/H9AABnE9dzr0qaiB6PS9ptbmwak9/W4bDNXUlHJ13d/UZcAVHvvpi0YABAvLie+3NJ09HjaUlrktTUQ5/Xp+FfybpAAEB6HcO9cYI0utyx2nTCdD36+kzSdHS5pNz9RbcKBQAkF3udu7s/a/Hc1ehrVYcBL0kEOwD0CT6hCgABItwBIECEOwAEiHAHgAAR7gAQIMIdAAJEuANAgAh3AAgQ4Q4AASLcASBAhDsABIhwB4AAEe4AECDCHQACRLgDQIAIdwAIUOzNOoCzWFl9o+/Xfznx/L2bl7Rcms2horCxv9Fg7p5vAWYbxWKxuLGxkWsd6K7Pv/5RkvTu0a2cKxkO7O+hYe0aGJYBgAAR7gAQIMIdAAJEuANAgAh3AAgQ4Q4AASLcASBAhDsABCj2E6pmtiCpKqng7t+1aL8bPZxx9z9kXB8A4BQ69tzNrCBJ7r4mqdpYbmqfl7Tm7s8kTUfLAICcxQ3L3NFhr12SKpKOh/d003OVaBkAkLO4YZlxSXtNy5PNjVGPvaEg6Xlzu5ltJKjhSoJ1AAApZHJCNRquWXX37Sy+HwDgbOJ67lVJE9HjcUm7bdabb3Wy1d1vxBUQ9e6LcesBAJKL67k/12/j6NOS1iTJzMYbK5jZ3Uawc0IVAPpDx3BvDLNEoV1tGnZZb3r+WzPbMbO/dbVSDKT9Wl2PV18fLRcevtTj1dfar9VzrCpc7G80xI65u/szd19rPnnq7lejr2vu/jt3n4m+rnWzWORvZfVN4nX3a3XdfvpK5c3K0XN7Bx9U3qzo9tNXqQInzXZDkfZnZn+jGZ9QRSqtbuHWTnlrR+93D1Srf/zk+Vr9o97vHqi8tdOV7YYi7c/M/kYz7qGK1Bq3cDuLWv2jnqy/1ZP1txlUFK4s9rXE/h5G9NwBIED03JFa0psuFx6+1N7Bh7btE2Pntf1NKdH3yqoHO2jS3OCa/Y1m9NzRNYtfTGl0pPUhNjpyTovXL/a4orCxv9GMcEcq925eSrzu0tyMpiYvnAic0ZFzmpq8oKW5ma5sNxRpf2b2N5qZu+dbgNlGsVgsbmxs5FoHumO/Vld5a+foRN7E2HktXr+opbkZjY0yKpg19vfQsXYNvNroqrHREd0vXT4Km6Rjvjgd9jcaGJYBgAAR7gAQIMIdAAJEuANAgAh3AAgQ4Q4AAeJSSABntrL6puVMkvduXtJyaTaHikC4Aziz5dKslkuzR3PSpJkTB93BsAwABIhwB4AAEe4AECDG3NFVx0+0NcZkOdHWHexvNBDu6KrGiTb0BvsbDQzLAECACHcACBDhDgABig13M1sws3kze9BhnUK2ZQEAzqJjuDdC293XJFVbhbiZzUv6U3fKAwCcRlzP/Y6kavS4Imn++ApR8O9lXBcA4Aziwn1cnwb3ZBdrATCg9mt1PV59fbRcePhSj1dfa79Wz7Gq4cYJVQAnrKy+Sbzufq2u209fqbxZOXpu7+CDypsV3X76KlXAp9kuOosL96qkiejxuKTdNN/czDbi/km6kr5sAN3UavredspbO3q/e6Ba/eMnz9fqH/V+90DlrZ2ubBedxX1C9bmka9HjaUlrkmRm4+5ebfu/AAy8xtQFZ1Grf9ST9bd6sv42g4qQRsdwd/dtM7sWXRFTdfftqGld0lXp8FJJSdfMbMHdXxz7/zfiCoh678VT1A4AaCN2bhl3f9biuatNj19IenF8HfQX7pSDtJLecKPw8KX2Dj60bZ8YO6/tb0qJvlcWfy3gEBOHDQnulINuWfxiSuXNyokxd0kaHTmnxesXc6gKXC0D4IR7Ny8lXndpbkZTkxc0OvJpnIyOnNPU5AUtzc10ZbvojHAHcEKaobqx0RH98NWXWipOHz03MXZeS8Vp/fDVlxobTT5AwBBhdhiWAXBmY6Mjul+6fHRVTNIxdnQPPXcACBDhDgABItwBIECEOwAEiHAHgAAR7kOCKVmB4UK4DyimZEU/WVl988nUAZ9//aM+//pHjpcccZ37gPp+/ZfEH/hIMiXr/dLlzLeL4dGY3gL9g3AfYEzJCqAdhmUAIED03AcYU7ICaIdw76G85lRnSlZg+BDuPZTlnOppp2T96edfT5xUZUpWIFyMuQ8opmQF0Ak99yHBlKxAtvr91pWEOwCcQr/funIow73ff+MCwFkNZbj3+29cADgrTqj2UJ6TdzH3BzBchrLnnpWV1TeJh3Eak3e93z04eq4xeddPP/+a+qqVNNuWmPsDYWKItT167mfQ6qBqJ8nkXd3aNhCq5dLsJ8Oq7x7d0rtHt4Y+2CV67mfG5F0A+hE9dwAIUGzP3cwWJFUlFdz9u7Ttoctj8i6JCbwAdNax525mBUly9zVJ1cZy0vZ+lcdVK4tfTGl0pPXuZvIu4HTyvAKt329dGTcsc0eHvXJJqkiaT9neE3ndcu5f/mki8bpLczOamrxwIuBPM3mXxAReCFNe72VJulP+91y23a3LkePCfVzSXtPyZMr2nsjrqpX//J+9+JUiWU7eJTGBF8KU5xVoad7PWW67W1e+dfVqGTPbSLDalSy2NQhXrTzbqnzyvff2/64n629lMsIaiOT5Xh6EHEkqLtyrkhpjD+OSdlO2owkfJALQK3Hh/lzStejxtKQ1STKzcXevtmtvcPcbcQVEvfti4orb4JZzQBjyvAItpBzpOObu7tuSZGbzkqqNZUnrMe19i6tWgDDk+V4ehByJ/RCTuz9z9zV3f9b03NVO7b2W9pZzWV21whUrQLbyei9L+V391q0cMXfvyjdOXIDZRrFYLG5sbPRsm/u1uspbO/rzf/yv/nbwd/3uwnktXr+opbmZ1FetAMhPnu/lPskRa9swjOEOAIFoG+7MLQMAASLcASBA/TAs83+fffbZP165kslnmQBgaGxubn7v7v/aqq0fwv2/Jf2DpNN+nKvxW+Gv2VSUGepKh7rSoa50Qq3rr30b7mfVmOIgyQemeom60qGudKgrnWGsizF3AAgQ4Q4AASLcASBAhDsABIhwB4AAEe4AECDCHQACRLgDQIAG/kNMAICT6LkDQIAIdwBnYmaFDm0LZjZvZg96WVO07U51fRt9vdu7inqLcO+Sfj6w+vXNmLe4nz3noIqrLZdjKrp/8p/atBUkyd3XJFU7HXe9rCty18x2JFV6VJKkw9cn+vdtm/bMjrGBDvde7qiUdfXlgSX19Zsx1/CK+9lz3jdJtp3LMRXVtNem+Y6kavS4Imm+J0Upti5J+r27z0Tr9UT03mvcb3o6Wm5uz/QYG9hw7/WOSqMfD6yGfnwz9kl4xf3suQVVwm3ndkx1MK5Pj7XJvAppoZBDx29av712lWi5WabH2MCGu3q8ozKWx4GVRF5vxn4Ir7ifPc+gSrLtfj2m+pK7fxcdS5PHO4Zd3OazqDMqSQVJ/3VslUyPsYEN917vqCzlcWD1OcLrjPr0mKpKmogej0vazbGWI9EQ4EK0uKuTHcNub78gadXdt7u5nZFufvNe6NWOOrbNVuO+lSS9ysZB5e4v1IUD6yy1qU/fjNJheEmSmZXMbL4LPfi4nz3PfdNx290+ptIys3F3r0p6Lula9PS0pFyHjJrqqui34b0ZSeUelzLfOJ6PyfQY6+twTxhUPdlRzZr+YkisVwfWGWvr2psx5rXsh/Bq+bP3SVDF1ZZbWEWvzTUzW4heH0lal3TV3bfN7Fr0l0S1xx2wuLrumtmepJ1edwybOirz7r7WtWPM3Qf2n6S7TY/no6/j0ddCo13SA0mFHta1IOlvkhaanvtLc93ROg9y2GdJaptv3rc9qKnla3XstWw8LnfrtWz1s+e9b1LWlssxxb/Er+F89N7bib42Mqsrx9jATj8Q9Qb+TYdjtRM6POG2ZmZ/cfer0Tp3FZ1s9VP0aNE7rV6rFq/lXtTe6i81AE0GNtwBAO0N7NUyAID2CHcACBDhDgABItwBIECEOwAEiHAHgAAR7gAQoP8HjTjgyzb6GX8AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Create a random sample\n", "raw = np.random.normal(0, 2, 1000)\n", "cut = np.random.normal(0, 0.5, 100)\n", "\n", "# Create the histograms (we do not care about the errors for the moment). Note that the two\n", "# histograms have the same number of bins and range.\n", "h_raw, edges = np.histogram(raw, bins=10, range=(-2, 2))\n", "h_cut, _ = np.histogram(cut, bins=10, range=(-2, 2))\n", "centers = (edges[1:] + edges[:-1])/2.\n", "\n", "ex = (edges[1:] - edges[:-1])/2.\n", "\n", "# Calculate the efficiency and the errors\n", "eff = h_cut.astype(float)/h_raw\n", "ey = hep_spt.clopper_pearson_unc(h_cut, h_raw)\n", "\n", "plt.errorbar(centers, eff, ey, ex, ls='none');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Displaying the correlation between variables on a sample\n", "The hep_spt package also provides a way to easily plot the correlation among the variables on a given sample. Let's create a sample composed by 5 variables, two being independent and three correlated with them, and plot the results. Note that we must specify the minimum and maximum values for the histogram in order to correctly assign the colors, making them universal across our plots." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEQCAYAAAC0v9O7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3xcZb3v8c8zuU4uM+mFxt2mCRfb9E4pFystBTYF3aXi62gVlVe3CoIIHBHOFj0ewY3sLXuDspUDyPYCbt2AYtUX2KKUIrW0tceW0hYSSRBp0hRIadrMJJPJ/Tl/zGSS1UySSduZNTP5vl+vvjrredbT/pquWb/1XNZaxlqLiIjIAI/bAYiISHpRYhAREQclBhERcVBiEBERh1y3AxiJMWZzAru9F1hnrf1SksMREZkwTLquSkowMSyeu2Ce/5a7/3eyw8kIU20+/flH3A4jPXRNIZzX5nYUacHXW0ROwWG3w0gP3adQXPiu21GkhYtW3mRGqkvbHoO19qKx9jHGbO6j/8InS15KQUTp7/quatpP/6nbYaSFnL9+jj3Tn3U7jLRwcfMFTKq+z+0w0kL367eyYsH9boeRJm4asUZzDCIi4qDEICIiDkoMIiLioMQgIiIOSgwiIuKgxCAiIg5KDCIi4qDEICIiDkoMIiLioMQgIiIOSgwiIuKgxCAiIg5KDCIi4qDEICIiDkoMCWr+4+u0vNTI/id2jbrf0PpE24hI9tm9r3PEunXr29i0pYN7HzwyaplblBgSEKw/BMCUsyvJLSmIbR+r5aVGjrzUOK42mWjn71uo2dbKMz88OKxuf007n5n9J758yW6+fMlufnL7GwBs/nkzm3/ezJP3NqQ63KSqf+4ADTua+fOjr8Wt3/LdvQDsW/fGsDZDy7LBC8+E2Lk1zGMPt45a/9TjwVjZzq1hdm4N89Dd7p8MT6ZNWzq48rp34tYNJIyVK4rw+3LYva8zbpmblBgS0PxCHbklBQB4p/tjJ/+T3SYT7K9pB2D+sjKKSnNj2wNCrb38pP793Pv8Em68fzarrp1BzbZW5p3v56JPlHPoQCc12+KfODJN81+OAlC1tJzC0rzY9lD7fvU3frR6A/6Kklgbf0UJVUvL8VeUxG2Tiepe7QLg3OVeSnye2PbQ+umVuZy73Mv0yjzqXu1i59YwL2wIce5yL3U1XcPaZLKVK4o4rSr+e9CefLqdMl8OAKdX5fH8ix1xy9ykxJCA3vYu8koLY9s9wfCwfYL1h5hyduW42mSiPz/TQlFp5AA+ZWYBtdsDjvr5y8pin/e/EmJaZSHvHuiK7TdtZmQ7G9Q920hBaT4A/ooSGnY0D9vnsjvO4XPrL6dqaXmsbMt39wEQaGqnfO6k1ASbZM//NkSJL3I6mV6Zx66tw4/3h/4t0it4q7GH6gUFnLvcy213T42W9VK9oCB1AbuoNdDH5LLBU2/L0f64ZW5SYjhJetrc7fqlSkewl+KywSuh9tbeuPvVbGvl3FVTALjoE+Vc9InIibGhJsSpC4uTH2gKdLX1UOjPj213BoYnvNamkGOoqXzuJMoqinnggt842ma69mA/viEntsAxJ7bqBQXMmJnHBxc14CvLcdQ99nArt31rakrilMQoMSQgt6QgduLvbe8iz+d11B/bW0ikTbar2R6g2OfsSu+vaadqfjGnzi9xKarUO++zc6haWk5noIuGHc10BrspKM3jvKvnsvGbu2htah/7D8kCbYE+Snwe1t7g59+/epiDjT2xuquuL+Opx4O0BfpcjDB1yvw5HGmNJM7WYB9TJnnilrkp/iCYOJRfXE2wrhnOhvBbASZHk0BPeyd5JYWE3w4QfjtAT7CTnmAnwfpDI7bJdEW+XELRXkJHWx8lZfEPoYaa0LCy2u0BPv7lqqTGl0oFpXl0BroB6GrrptDvHArZt+4NCv35zL50JoX+AgJN7TT/5SjnXT2XQl8+ZRXF1D/XxHmfneNG+CdVic9DMHpiaw/24z/mxPb0E22svcFPqT+H6ZV5bH4mxDnLIxdL1QsKmD4zj6efaOOq68uG/dnZojXQR5k/h49fUcKuvZHe5ZsNPVxyQRFA3DK3qMeQAN/saUBk1VFuSUFse/etvwag/MJZlF84C4DeUNeobTLdeaumxOYI3j3Qybzz/QCEgoNDSocahw+rbf55M6uunQGQNZPP1R+oJBC94m9tCsXmETqDkWRRPn8yle+LlAWa2imfP9nRfvalMykszUthxMlzyYeKeSvaC3irsSd20o/XC7h4VTElPg+7toYdyWR6ZXb8LCCy9PSlvV2sW98WK7v0Y5FVfEsWReYeN23pwO/LYcmiwrhlbkqLHoMxZhZwHvC4tda6HU88FR9aGPlw9mDZ+37wqWH7xPYboU2mO3V+CftfCVGzrZWi0tzYsNA9n67lzt8siu13yszBq+eaba08+e0GNvzwIKFALzd+b3bK406G8rmTaK45QsOOZgpL82ITyb/8/B9Z+8SllM+dFOs1+CtKKJ87ifK5k/jzo69RVlFMZ6CbRWvOcPlfcXJULyjgtX2RlUYlPk9sIvnmq97hkfUzuOr6Mh57uJXplXkEW/v48Kd8tAX6+MOGEE89HkkoF6/KjrkngDWrS1mzutRRtnPj4KjBdWv9w9rEK3OLcfs8bIw5F/ADPqDeWvvqONpunr1gzoWV//eDSYsvk1zfVU376T91O4y0kPPXz7Fn+rNuh5EWLm6+gEnV97kdRlrofv1WViy43+0w0oLnPfVmpLp06DGUAEXW2l8DGGNKrLUTY0ZORCQNpcMcQx9wpjHmNGNMDnC+MWbhWI1ERCQ50iExbCUSx6eAfwF6rbWvuBuSiMjE5epQkjEmx1rbZ4z5JlAETLPWvulmTCIiE11KegzGmLjr0KJJwdiI0EBSMMaMOCkiIiLJlfTEEJ03uMAYs3hI2QeNMQsA4i1PTdclqyIiE0HSh5KivYJa4H3GmDbgDKAM0FpCEZE0lIoeg8da+w7wPLCWyO1ev7TWWg0ZiYikn6QnBmvtwGMWzwaagNeBU6MJQ0NGIiJpJiWrkowxFcBka+2PjDElwBVAHhD/tVciIuKak5YYBpaejlDdbK39DYC1tt0Y8yyQHW9rERHJMic8lGSMORMik8zR7ZJjl6daa3uO2W7RYy9ERNLTCSUGY8wZwB3GmLui26cCHwF+bIyZEy1Lh7urRUQkQSd60m6w1n4U6DbG/BCYbq39KfBH4E5jzNQhk88iIpIBTjQx9ANYa+8CeoCPRrd/DGwBVmtJqohIZjmhxGCt7Y/e2Yy19gYgaIy5J1q9FXhVS1JFRDLLic4xTAfKB7attXcCecaYbwPt1tpdJxifiIik2IkuVy0BqoC3jDG51tpea+0txphl1to3TuQPNsZsTmC3xWPvIiIi43FciWHgiahE5hi6Aay1sbfBW2u3HbNf0uTYHK44tDSZf0XGyPMapu3/hNthpIW2zklU1X3E7TDSgvV66H79VrfDSAuBjnJ+t++LboeRFi5/z8h1404MxphyIo+0OA2YBlQaYyYBxcAma23zwL4nkhSstRclEMvmHtt/4X3NTcf712SVh+bkUTHnnrF3nAA2/+k76LiIeHRugd5zHPW7fV/U+69jbhqx5nh6DEettc3GmL8A04FKIhPN/qFJQUREMlPCicEYUwX0WWubAKy1wWjvodda2wF0JClGERFJoYRWJRlj5gEXAU8YY84ZUtUFHElCXCIi4pJEewyLgeeITDbfGb1Xoc5a2wg0Jis4ERFJvUTvY2gBlgN7gX8FPg+claygRETEPQklBmvts8CfiCSHbuDfrLW/S2ZgIiLijvHc+dxD5OU65dbafRC5TyEpUYmIiGsSXpVkrX3HGPOwtbZrSJmegyQikmXG9aykoUlBRESyk16iIyIiDkoMIiLioMQgIiIOSgwiIuKgxCAiIg5KDAkKvbyXcF09gU0vjKu+64Ae/SwTx+59nSPWrVvfxqYtHdz74JFRy8R9SgwJGDi5e6tn4/F6h53sR6oP19Xz7iM/S22wKbBxQ5gdL3byyPfbRq1f91goVlb7SjcbN4TZuCGcqjBT4nguGNq27aBt2w6OPLU+VWGmxKYtHVx53Ttx6wYSxsoVRfh9Oeze1xm3LFu88EyInVvDPPZw66j1Tz0ejJXt3Bpm59YwD93tfpJUYkhAaPcePEVeAHKnTqaz7vWE6r3Vs8mdOjm1wSZZ7SvdACy9oBCfzxPbHlpfUZnD0gsKqajMidX/+ME2LrvcS1Nj77A2mep4LhjCdfUUVs+idNlSeluOEK6rT3ncybJyRRGnVcW/Z/bJp9sp8+UAcHpVHs+/2BG3LBvUvRq53evc5V5KfJ7Y9tD66ZW5nLvcy/TKPOpe7WLn1jAvbAhx7nIvdTVdw9qkmhJDAvrDnXiKimLbfaHQuOqzybO/DePzRZ6EUlGZw46tww/g794duQpqauxj3sJ8Nm4IM//MfACu/kIp8xbmpy7gJDqeC4bew0di++VOmUzvYfevDlOhNdDH5LLB003L0f64Zdng+d+GKPFF/l3TK/PYtXV4L/mhf4v8v7/V2EP1ggLOXe7ltrunRst6qV5QkLqA41BikHFpC/bjG/JlDhzzZZ63MJ+KyhyWLXgrtl/N3m4CR/upfaV7xOGnTHQ8Fwyly5ZSuizyjvLupoPkV1akJlhJmfYxviPVCwqYMTOPDy5qwFeW46h77OFWbvvW1JTEORolhgR4vIX0d0S6uf0dYXKKi8dVP5EEA/2U+jxcc2Mpd37lKE0NvQD4J3liPYVsm2c4Hl0HmsivmEHBzImRGMr8ORxpjZwgW4N9TJnkiVs2EbQF+ijxeVh7g59//+phDjb2xOquur6Mpx4P0hboczHC43vn84RTvGQx3Y1NUA29LUcorJ4FQF9HmJwi74j12ajU5yEY/TIHgxb/MV/mdY+HuObGUnx+DxWVuWx8Jox/UuQzgM/noWZvN5dd7k157CfbiVwwdNa9zuQPr05dsC5pDfRR5s/h41eUsGtvZNjxzYYeLrkg0pOKV5bpSoZ8R9qD/cO+I08/0cbaG/yU+nOYXpnH5mdCnLM88n2oXlDA9Jl5PP1EG1ddX5by2Ae4nqIz4dHdA1d14bp6PF5vbLv5gYdHrQ+9vJfuxiZCL+91Ierk+MCHvDQ1Rq5mmhp7Wbo8MhYaDAwfH77sci8+n4fLVkUmnQGCwf7YfEOmK16yODZHcOwFw2j1bdt24F95MUBWTT6vW9/GS3u7WLd+cLjw0o8dBGDJokIgsnLJ78thyaLCuGXZ4JIPFfNWtBfwVmNP7KQfrxdw8apiSnwedm0NO5LJ9Mq81AUch+s9hqGP7jbGmHR9lPfAuDDVg2XTb7tl1Pris86k+KwzUxBd6sxbmE/tvh52vNiJzzc4PHTtJw/zi2emcfUXSnnk+21UVOYSbO1nzVWRq2Sfz8PGDWECR/u5+guZ31uAyAVBd2NT3AuG6bfdErc+XFfP0ac3ENj0Av0dHZxy9VqX/xUnz5rVpaxZXeoo27mxMvb5urX+YW3ilWW66gUFvLYvstKoxOeJTSTffNU7PLJ+BlddX8ZjD7cyvTKPYGsfH/6Uj7ZAH3/YEOKpxyMJ5eJV7g5HG7fOw8aYa4A51tovR7fPBOYCz1prjyb4Z2x+77y5F/Zef00SI80cD83Jo2LOPW6HkRY2/+k73NesmwsBHp1bwIoF97sdRlr43b4vMqn6PrfDSAvnV/1txNEaV4aSjDEFwC6gzRjzDWNMNVAITAV+YIyZ50ZcIiLiUmKw1nZZa/daa78JTAa+Za39f9baB4DfAbcbY9wdZBMRmaBcm3w2xngArLU3A7XGmDuj248AzwE5ozQXEZEkcXNVkm9Icrgd8BpjHjDGzAXWW2uz58EpIiIZxK05hmJgOVAwsFzVWnsb0AC0WWsPuRGXiIi4tFzVWhuKJgTH8lRr7b1uxCMiIoNSmhiMMdOAaUCQyIr/VmPMQaDXWtuYylhERCS+lCUGY0wu4APeAxQDpwB/B/QBSgoiImkiJYnBGOO11oaBv0Z/YYyZBPzBWpsdD+cXEckSSZ98NsYsBb4e/X3os5FyAd2rICKSZpKaGIwxs4nczfwscI8xZtGQyebN1trsfaONiEiGSnaPoQdYCBwFfgLca4z5kjGmyFrbnuS/W0REjkOy5xjeBn4MLAaeBJqAOmttdrzcVUQkCyU1MVhrO40xfcDZwDLgfmttSzL/ThEROTFJX5Vkre0xxvxn9PPEePO5iEgGS8lyVSUEEZHM4fob3EZijNmcwG6Lkx2HiMhEk7aJIVF5xsMXyyvcDiMt9HQaml67ze0w0kJ5Xj636rgAoD3sYcurX3Q7jLTQ21nO0bpb3Q4jPVSNXJW2icFae9FY+xhjNhtP94XexQ+kIKL017DnJu7b3+N2GGnh0bkHWaPXWQLw4Ja7uG9/l9thpIW73tvGnvIX3Q4jLVzOTSPWufk+BhERSUNKDCIi4qDEICIiDkoMIiLioMQgIiIOSgwiIuKgxCAiIg5KDCIi4qDEICIiDkoMIiLioMQgIiIOSgwiIuKgxCAiIg5KDCIi4qDEIHKCdu/rHLFu3fo2Nm3p4N4Hj4xaJpJOlBgS9NLvD1O7vZXf/7BpWF1DTTvXVm/layt38bWVu/jZHX8ds00mC728l3BdPYFNLyRc37ZtB23bdnDkqfWpCjMlNm3p4Mrr3olbN5AwVq4owu/LYfe+zrhl2ULHxaD65w7QsKOZPz/6Wtz6Ld/dC8C+dW8MazO0zC1KDAloqGkHYN75ZRT5cmPbA0KBXn5Yt5xvbTqHz39vDh+8tmLMNpmq60AkyXmrZ+PxemPbo9WH6+oprJ5F6bKl9LYcIVxXn/K4k2XliiJOq4r/vqsnn26nzJcDwOlVeTz/Ykfcsmyg42JQ81+OAlC1tJzC0rzY9lD7fvU3frR6A/6Kklgbf0UJVUvL8VeUxG2TSkoMCdj5zLt4fZEv/9SZhfxle6ujft75ZbHPDa+2c8rMwjHbZKrQ7j14irwA5E6dTGfd62PW9x4+Etsvd8pkeg9PjCGU1kAfk8sGv2ItR/vjlmUDHReD6p5tpKA0HwB/RQkNO5qH7XPZHefwufWXU7W0PFa25bv7AAg0tVM+d1Jqgh2BEkMCwm19FPsHrwrbW+O/PrN2eyvn/MPUcbXJNP3hTjxFRbHtvlBozPrSZUspXbYUgO6mg+RX6l3M2UbHxaCuth4K/fmx7c7A8NeqtjaFHENN5XMnUVZRzAMX/MbR1i1KDCdR7bajFPnS9jXarus60ER+xQwKZmbHCWAsZf4cjrRGegStwT6mTPLELZvoJtpxAXDeZ+dQtbSczkAXDTua6Qx2U1Cax3lXz2XjN3fR2uTu0LPOYgnwluYQCvQCEA72UlKWF3e/xtrQuNtkGo+3kP6OyLh4f0eYnOLihOs7615n8odXpy5Yl7QG+ijz5/DxK0rYtTdytfhmQw+XXBC5Yo5Xlul0XAwqKM2jM9ANQFdbN4X+Akf9vnVvUOjPZ/alMyn0FxBoaqf5L0c57+q5FPryKasopv65Js777Bw3wgfUY0jIuatO4fCByOqRdw90Mjc6p9AR7I3t8+6BzoTaZLriJYtjY8G9LUcorJ4FQF9HeNT6tm078K+8GCBrJhkhsvT0pb1drFvfFiu79GMHAViyqBCIrFzy+3JYsqgwblk20HExqPoDlQSiV/ytTaHYPEJnMJIsyudPpvJ9kbJAUzvl8yc72s++dCaFpe5eSKZNYjDGuD+wNoKq+ZGVA7XbWyny5ca2v/OZVxz7nTKzcMw2mW6gux+uq8fj9ca2mx94eMT6cF09R5/eQNOdd9P4ldvdCTxJ1qwupeW1M1izujRWtnNjZezzdWv9rFxRxHVr/aOWZTodF4MGJo4bdjRTWJoX2/7l5/8Yq6/feID65w7gryihfO4kzvvsHPb96m/UP3eAfeveYNGaM1yLH8BYa10NAMAYMx1YC/zCWrt/HO02z1s468Jb1pWPvfMEEN5zE/c1Z9c9E8fr0bkFrFhwv9thpIUHt9yl4yLqrvdOYs/0Z90OIy18+8xfmJHqXJ9jMMbMAuYBvwJGDFRERFLD9cQA5AGlwBvWWmuMmQb0WGvdvcNDRGSCSoc5hiPAOcAFxphK4GwgbecbRESyneuJwVr7DvAD4H3ADUCttXb4rYIiIpIS6TCUhLW2Fqg1xhibDrPhIiITmOs9hqGUFERE3JdWiUFERNynxCAiIg5KDCIi4qDEICIiDkoMIiLioMQgIiIOSgwiIuKgxCAiIg5KDCIi4qDEICIiDkoMIiLikBYP0YvHGLM5gd0WJzsOEZGJJm0TQ6Jsfz6de25yO4y0cJrXw6NzC9wOIy0capvBg1vucjuMtFDlNTw0x92Xy6eLrnAp73/nYrfDSA9njlyVtonBWnvRWPsYYzb32P4L73v7YAoiSn+PzM/Xe46j9J7jQQ/NyaNizj1uh5EW6mu/TvvpP3U7jDRx/Yg1mmMQEREHJQYREXFQYhAREQclBhERcVBiEBERByUGERFxUGIQEREHJQYREXFQYhAREQclBhERcVBiEBERByUGERFxUGIQEREHJQYREXFQYpDjsntf54h169a3sWlLB/c+eGTUMhFJT2n7PoZ0FNqzF4/XS1dTE2WX/H1C9WO1yUSbtnTwhdsO8fqOU4fVDSSMlSuK+FtDjyOBDC1bsqgwVeEmVejlvXiKvHQfOIh/5fAXwMSrb9u2A4Cew4eZ/OHVKY03mTZuCOPzGWpf7eHqL5SOWN/U2Meaq4oBqH2lm6bGPgAuu9yb0niTaefvWygqzaGhNsSqa2c46vbXtPPP/+MVTpkZeanW/PP9fOauM9j882YADh3o5ONfrkp5zEOpx5CgrgORl754q2dHTvQHmsasH6tNplq5oojTquJfUzz5dDtlvhwATq/K4/kXO+KWZYPjOSbCdfUUVs+idNlSeluOEK6rT3ncyVD7SjcASy8oxOfzxLaH1ldU5rD0gkIqKnNi9T9+sI3LLvfS1Ng7rE2m2l/TDsD8ZWUUlebGtgeEWnv5Sf37uff5Jdx4/2xWXTuDmm2tzDvfz0WfKOfQgU5qtrW6EXqMEkOCQi/vweONXNHkTZlCuL5+zPqx2mSj1kAfk8sGD6uWo/1xy7JBaPcePEWR/9/cqZPprHt9zPrew0di++VOmUzv4ewYWnv2t5HeAEBFZQ47tnYN2+e7dwcBaGrsY97CfDZuCDP/zHwArv5CKfMW5qcu4CT68zOR3gLAKTMLqN0ecNTPX1YW+7z/lRDTKgt590BXbL9pMyPbblJiSFB/OIynqGhwO9QxZv1YbSSz9Yc7Hf+/faHQmPWly5ZSumwpAN1NB8mvrEhNsEnWFuzHNyT5B45J/vMW5lNRmcOyBW/F9qvZ203gaD+1r3TzyPfbUhpvMnUEeykuG+xRt7f2xt2vZlsr566aAsBFnyjnok+UA9BQE+LUhcXJD3QUSgxyUpX5czjSGjkptAb7mDLJE7dsous60ER+xQwKZmZHYhhLMNBPqc/DNTeWcudXjtLUEDlZ+id5Yj2FjRvCboaYcjXbAxT7nEOy+2vaqZpfzKnzS1yKKkKTz6M4/Itf0tfRQUHlTDxeL/0dkSv+/nAYT3GRY9+R6kdrk01aA32U+XP4+BUl7Nob6Qa/2dDDJRdE/s3xyjKdx1s4+P/bESanuDjh+s6617Nq4rnU5yEYTf7BoMV/TPJf93iIa24sxef3UFGZy8ZnwvgnRT4D+HweavZ2Z8UEdJEvl1C0l9DR1kdJWfzTbENNaFhZ7faA6xPPkAY9BmOMcTuGkUy98mOUf/bTlF3y9xSftZielhYAelpa8M6eDUBfR+QqJ179SG0y3br1bby0t4t16we7/5d+7CBAbLXRpi0d+H05LFlUGLcsGxQvWRybI+htOUJh9SxgyDExQn3bth2xFUrZMvn8gQ95Y6uLmhp7Wbo8suImGBg+n3TZ5V58Pg+XrYpMOgMEg/2x+YZMd96qKbE5gncPdDLvfD8AoeDgkNKhxuHLvTf/vDm2gmnCTj4bY641xqy21tro9hnGmPe7Fc9YBrr84bp6PF5vbPudhx4esX6kNpluzepSWl47gzWrB5ck7txYGft83Vo/K1cUcd1a/6hlmW6k/9/mB0Y+JsJ19Rx9egNNd95N41dudyfwJBgYDtrxYic+3+Dw0LWfPAxEJpfXPR5i44Yw6x4LseaqYiqqcvH5PGzcECZwtD8regtAbBioZlsrRaW5se17Pl3r2G9guerAvk9+u4EvX7KbG875c+qCHYGJnpdT+5caUwR8BPgA8F9AHbAQ+CTQANxpre1J4M/Z/N55cy/su+5zyQw3YzwyP58VC+53O4y08OCWu7ivOTuWB5+oh+bkUTHnHrfDSAv1tV+n/fSfuh1GWvj0rO0jjta40mOw1nZYa/8buAP4CvCP1tpnrLVrgXLgZjfiEhERl+cYrLVvAv8TWGiM+UC0+BYgYIzRxLiIiAtcP/laa18zxvwf4NvGGB/wMvBba238xb8iIpJUrq9KArDWvgF8FbgECFhr33E5JBGRCcv1HsMAa22dMeaL1trseGCKiEiGSosewwAlBRER96VVYhAREfcpMYiIiIMSg4iIOCgxiIiIgxKDiIg4KDGIiIiDEoOIiDgoMYiIiIMSg4iIOCgxiIiIQ9o8K+lYxpjNCey2ONlxiIhMNGmbGBKVZzzc/Hcz3A4jLbSHPWx59Ytuh5EWyvPyubU8O16leqJ6Og1Nr93mdhhpIadrCiV/+0e3w0gPs0auStvEYK29aKx9jDGb83O6LrxhRfa8O/dEPLTlLu77m55DCHDPnEN4T/tvt8NICw17buK+/WO+KXdC+NppXTxZUOd2GGnho6PUaY5BREQclBhERMRBiUFERByUGERExEGJQUREHJQYRETEQYlBREQclBhERMRBiUFERByUGERExEGJQUREHJQYRETEQYlBREQclBhERMRBiWEcdu/rHLFu3fo2Nm3p4N4Hj4xaJpKtQi/vJVxXT2DTC+Oq7zrQlIrwUqr5j6/T8lIj+5/YNep+Q+sTbZMKSgwJ2hcu+pYAAAqoSURBVLSlgyuveydu3UDCWLmiCL8vh937OuOWZZPQnsiXvPX5PyRcP1abTPTS7w9Tu72V3/9w+Mmtoaada6u38rWVu/jayl387I6/jtkmUw2c3L3Vs/F4vcNO9iPVh+vqefeRn6U22CQL1h8CYMrZleSWFMS2j9XyUiNHXmocV5tUUWJI0MoVRZxWFf+9Rk8+3U6ZLweA06vyeP7Fjrhl2eJ4TgJjtclEDTXtAMw7v4wiX25se0Ao0MsP65bzrU3n8PnvzeGD11aM2SZThXbvwVPkBSB36mQ6615PqN5bPZvcqZNTG2ySNb9QR25JAQDe6f7Yyf9kt0kmJYaToDXQx+SywR9ly9H+uGXZIvTyHjzeyJc8b8oUwvX1Y9aP1SYT7XzmXby+yMXC1JmF/GV7q6N+3vllsc8Nr7ZzyszCMdtkqv5wJ56ioth2Xyg0rvps0tveRV5pYWy7Jxgetk+w/hBTzq4cV5tUUmKQcesPhx1f8v5Qx5j1Y7XJROG2Por9g73I9tb4r8+s3d7KOf8wdVxtJLv1tKX30LISw0lQ5s/hSGukR9Aa7GPKJE/cMpmYarcdpciXtq9XPyk83kL6OyLJvr8jTE5x8bjqs0luSUHsxN/b3kWez+uoP7a3kEibVMvuozXJWgN9lPlz+PgVJeza2wXAmw09XHJB5Mo4XlmmOvyLX9LX0UFB5Uw8Xu/glzwcxlPs/LeNVD9am0zkLc0hFOgFIBzspaQsL+5+jbWhcbfJNMVLFtPd2ATV0NtyhMLqWQD0dYTJKfKOWJ+Nyi+uJljXDGdD+K0Ak6NJoKe9k7ySQsJvBwi/HaAn2ElPsJNg/aER27hFl7EJWre+jZf2drFufVus7NKPHQRgyaLI2OCmLR34fTksWVQYtyyTTb3yY5R/9tOUXfL3FJ+1mJ6WFgB6Wlrwzp4NRE4CQNz6kdpksnNXncLhA5GrvHcPdDI3OqfQEeyN7fPugc6E2mS6gpkVQGSVkcfrjW03P/DwqPWhl/fS3dhE6OW9LkSdHL7Z04DIqqPckoLY9u5bfw1A+YWzKL8wkhh7Q12jtnFL2vQYjDF51tq0HXBds7qUNatLHWU7Nw5m9evW+oe1iVeWDQpmVtB14MCwL/k7Dz3MjH+6ZcT6eGWZrGp+CQ2vtlO7vZUiXy5V80sA+M5nXuH2X58V2++UmYVjtskGpcuWRj5UD5ZNv+2WUeuLzzqT4rPOTEF0qVXxoYWRD2cPlr3vB58atk9svxHauCUtEoMx5lTgWmPMb621O1wORxLgO//9AAwdCZ3xT7eMWh+vLNOtuPI9w8qOTQprv/neMduIpBPXh5KMMdOABcCzwD3GmEUuhyQiMqG5nhiAYmAhcBT4CXCvMeZmY0yBq1GJiExQ6TCU9DbwY2Ax8CTQBNRZa7tcjUpEZIJyPTFYazuNMX1EplyWAfdba1tcDktEZMJyPTEAWGt7jDH/Gf2sR5GKiLgoLRIDKCGIiKSLdJh8FhGRNKLEICIiDkoMIiLioMQgIiIOSgwiIuKgxCAiIg5KDCIi4qDEICIiDkoMIiLioMQgIiIOSgwiIuJgrLVuxxCXMWZzAru9v7i4MP+cRSbZ4WSEg4FTaerqdjuMtHB6kYdeb7PbYaSF/vYZHOjWU+wBqgrzOZTTNvaOE0Dr3oPfs9Z+KV5dpieG5UAf8KfkRpMRFkd/3+NqFOlBP4tB+lkM0s9i0GKg3Vob9+XrafN01WNZay8aa5+B5JHIvtlOP4tB+lkM0s9ikH4Wg8a68NYcg4iIOCgxiIiIgxKDiIg4KDGIiIiDEoOIiDgoMYiIiIMSg4iIOCgxiIiIQ9re+ZyujDF51toet+NwizHGWB00wxhj8q21E/p5JDo2BhljZgHnAY9n4s9EPYZxMMacCvyzMWapy6GknDHmWmPM6oGD3BhzhjHm/W7HlQ6MMdOBW6LHx4Q19ARojJmwDzAzxpwLVAFhYL7L4RwXJYYEGWOmAQuAZ4F7jDGLXA4pZYwxRUQO8iuNMSuNMTOBauAGY8y/GGPy3I3QPdErw3OBXwET8mRojLnGGHPvkO0ziRwrk1wMy00lQJG19tfW2leNMSVuBzReSgyJKwYWAkeBnwD3GmNuNsYUuBpVClhrO6y1/w3cAXwF+Edr7TPW2rVAOXCzqwG6Kw8oBd6w1r5pjJk2kU6I0eN/F9BmjPmGMaYaKASmAj8wxsxzNUB39AFnGmNOM8bkAOcbYxa6HdR4aI4hQcaYQsBH5KmE24HzgTprbYOrgaWYMWYO8M/Ao9baZ6NXQ5+Mbve6GpwLjDHvAb4K/BrYT2ToYLe1dsI989sY8z2gwlr70ej21cClRC4kJsy8nDHGA9wO9BLpPTxnrf2Du1GNjxLDOESHTP6JyBXR/dbaFpdDcoUx5gzg28DjwMtEHt/7jrtRuSd6VXw5MAX4/gS8WPBYa/ujn+8C+q2134huX01kArbTzRhTxRiTY63ti86xFAHTrLVvuh3XeCkxjJMxZjKAtfaI27G4KTpkcAtwu7X2XbfjSQcTdVWOMaYMCA5JDvcQOSk+CLRYaw+5GV8yjLY6Md5xkGnHhhKDHDct0RRjTDFwMfA80Dlk1dqXgSestU1uxpcM0XmDC4Ej1to90bIPAk3W2lddDe4kSdsX9Uj6U1IQa20oOmziuCK21t47SrOMFh0qqgXeZ4xpA84AyoisWMwKSgwiMm7R5dvTgCCRpcutxpiDQK+1ttHV4JIsOqfyjjHmeSJzjl3APdZam2lDRiPRclURGRdjTC6RFXrvAf4OOCX6+3uAfhdDS4mBuRTgbKAJeB04NZowMj4pgHoMIjIOxhivtTYM/DX6i+h9G3+YSEOLxpgKYLK19kfRJdtXELmn5TV3Izs51GMQkYREHwXz9YFHwgx57EUukZNiVolOMo+k2Vr7GwBrbTuR+YWsmWhXYhCRMRljZhO5mzn2SJghwyabrbUh96I7uaKP9MBa2xfdLjn2sS/HLlW11rZEE0RWUGIQkUT0MPyRMF8yxhRl0wkxevPmHdEb9QYenPkR4MfRu/4H7mzOarqPQUTGNFEeCWOMybXW9hpjbgcqiTzqZbsx5hrgMuBGa+1hd6NMPiUGEUnIRHgkzDGP93gICFtr/1d0+0YgBPxXtqw+GokSg4gkbCI8EmbgeUfRz98Aiq21t0XnHvKstbvcjTD5lBhERKKiL13CWvvWkLL/IPIo7e9ba99wK7ZUUmIQEYmKrr6qstY+NzDfEC1fZq3d5nJ4KZP1s+siImMZck9GP9ANMPT9IgNJYaK8slR3PovIhGaMKSfySIvTiDz/qTJ6N3cxsGnoS5eyfdJ5gIaSRGRCG3h8vDHGB0wnskx1K+C31r7tbnTuUI9BRCYkY0wV0DfwzghrbTDae+i11nYAHa4G6CLNMYjIhBN9HetFwBPGmHOGVHUBWbsUN1HqMYjIRLQYeI7IZPOd0deR1kXfJZHV75NIhHoMIjIRtQDLgb3AvwKfB85yNaI0oslnEZmQjDEzgQ8Bfwa6rbX7XA4pbajHICITVQ+R90iUDySFiXKfwljUYxCRCcsYU2Ct7XI7jnSjxCAiIg4aShIREQclBhERcVBiEBERByUGERFxUGIQEREHJQYREXFQYhAREYf/D0U44OLyrak0AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Create a random sample\n", "a = np.random.uniform(0, 1, 1000)\n", "b = np.random.normal(0, 1, 1000)\n", "c = a + np.random.uniform(0, 1, 1000)\n", "ab = a*b\n", "abc = ab + c\n", "smp = np.array([a, b, c, ab, abc])\n", "\n", "# Calculate the correlation\n", "corr = np.corrcoef(smp)\n", "\n", "# Plot the results\n", "fig = plt.figure()\n", "\n", "hep_spt.corr_hist2d(corr, ['a', 'b', 'c', 'a$\\\\times$b', 'a$\\\\times$b + c'], vmin=-1, vmax=+1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting a 2D profile\n", "When making 2D histograms, it is often useful to plot the profile in X or Y of the given distribution. This can be done as follows:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAD+CAYAAAAtUeIJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOy9eXhd1Xkv/NtHR9I5mmfJlmRZg0dhY2QDtsHYAZshAdLmEtIyJWkamrZfk6ZD+O59+nW4be99oLft/fq1aUNuG2ZCoEkDJAZsEoNjGzC4wjO2JVmW5UHzfI6ms78/1t7aw9nrfY/OliwpXr/n8WNpL629195rrXde76vpug4FBQUFBYVEEZjrASgoKCgoLCwoxqGgoKCgMC0oxqGgoKCgMC0oxqGgoKCgMC0oxqGgoKCgMC0E53oAmqb9bwDQdf3353osCgoKClcSBv1b5+MWjXNBO+eccQBYl4eirTsCn//GXA0gkJEhbYuNjMzavRMB93zu/n76z/azOczmu/ntz/UNlhTTzx4aJtv9wu+39wO/a57DXL7bTGBX7GXN9us6AFu3bgpP+z7vHIjM2Jimi/nAOBQUFBSuWtyyKYTdP1w87X63fa4d7x6IzsKIeCjGoaCgoDCH0AFM6rGk+s0VFONQUFBQmFPoiCXFBuaOdSjGoaCgoDCH0AHEoDQOBQUFBYVpYHKB5QxUjAP+ojS4CJJAVib9bJ/RNX6jeyhwY5/o6Ez63oD/6JtgVSX9BxHacRgoLiTbJ1rb5H25iC1mXmc7Kmu216WfPeP73Rd4VJUbMegYT0LjSM68NTNQjENBQUFhjjGXTCAZKMahoKCgMIfQ9eRMVXNp3VKMQ0FBQWGOMX1D1dxCMQ4FBQWFOYQOHZNJmKp05eNYuPDrqJttJ+lsOjnn3AHL+Oa5/v0PbyLbC9+l708+u7M7+c7gAw98p1NhAgNiPgIfuGenrFpOtk+cOJX0sxcixAHA5PrNFRTjUFBQUJhD6NAwDo3/Q49+cwXFOBQUFBTmEDqAmNI4FBQUFBSmg8k51B6SgWIcCgoKCnMIHckxDqVxKCgoKFylEKaqq4xxaJr2qPFjra7rj/m930LDbBcb4qJruAgVrbuXbKfG5/veXFoMJrInEKGjshAOkc1aRgXZnn9sgGz3k3KEA/dtOWgjTBEfLuLMZ9SXn/fXW8/7evYvG3RoGENKUv3mCr5qjmuath3Abl3XnwRQY/yuoKCgoJAgdF1oHNP9N5cnx30xDgA1AExm0Wz8rqCgoKCQMDRMJvEPCzUc19A0TDQAeMnfcBQUFBSuLogDgNOX4Re0jwMANE1rALBL1/VDrut7Eui+bibGoKCgoLAQoUNDLAnjj18fh821sMP0T2uadh+APgANuq4/Ies7U1FV26mHzHf4dXT6ejaTloNNDcE4GjXm/lRNi1g4lew7uraKbE/bf5xs1zPC9P2XlZLtodOX6fsz3yZy2zVke0rBBmlbejvtWOcCB8A4tznntcasCw5+AxOowAe/dVz8BpxwmOv7e+FKn+MwhP0duq4/pmnaY8bvAABd13drmlajaVqDWxkwMSNRVSbT0DRtu67ru20D2JZA/z0Atvodh4KCgsJCRAwaxvXpR1XFfDAbgyGYTKFG1/VDmqY9DmCXca0Zwn/tyThmIqrqcU3TmjRNY0QsBQUFBQUvTCIw7X8zAU3TvgXgt4xf8wD02Jqlaqlf5/huAPl+7qGgoKBwNUOHlqRzXAOAdYn4kmXWH13Xn9A07WVN0z6czrPVyXEFBQWFOYQOJOkcTx6mT8MwWTUDeBTCKV5g/EkeAKmjTTEOBQUFhbmEDkwmkXLE4ByNifiSPWD3X+QBOAhgNwAzIqTG+N0TC4JxcMWMuCgNDlQUhd9CSn7H7rdYEpeWg4psCnTTkUMpGXTU1WQDnVYjtY2OHAqdpiOPokzUVUolHTmU+X4L2R6rKJG2cSk/9ELagjtRSM8r/WUBRKJks98iWGzUFd2dBLcnuDWLwmq6/SQzrz6Lp3F7drpRWcI5Pn1S7Mc5DuBJAPcb4bfQdf0VANA0bYPhu+6TRVQBC4RxKCgoKPzyQkvS2e0rqqoPgnkAwCu2609693BCMQ4FBQWFOYQ4OX6VZcdVUFBQUEgec3Vy3A8U41BQUFCYYyQTjjuXWBCMw6/zmwPlDOOe7dt5TqT8AMA6QVlHItNfI9rHl5fT92aQMjJOtnMpRzgHdFoPU5OCAeX8BgAtIh8/N3Yu3UnqCJMyhHN+M2PXmbnTuhkHLlfvg1t3fsCt+ZOzW0uE27MzDRGOq0xVCgoKCgoJQtc1jCURVaUnE8I7Q1CMQ0FBQWEOoUNLsnSsYhwKCgoKVyV0IKlwXGWqUlBQULhqoSGWlHNcaRwKCgoKVyWExqGc4zMOv2k3uCgJP1EUXF8udYPvqCkGXGTUSJn8/jnHeqRtADBQX0C2d6+h563wSIxsH82jpbDMyxNke1rvGNk+GaKXf0pUnvhjoJqOqgrV5ZHtE+n0u2VcYtYFAy6ijU2Z4iPizW+0HLdmUw6dIttnG7NRCCo5jWPusCAYh4KCgsIvK2J6koWcVFSVgoKCwtWK5OpxKB+HgoKCwlUKdQBQQUFBQWFaEEkOVTjutKEFUhAIyx1OnLPJb359ysHtt3YA5ygcZWpKcAgdbiXbuZoXuafkTtjhG+m6B+l9tHN66Y9pJ2hvfQ7Zzjm/e5fRyzfcOX27sR1DlWnStqrn28i+40wtkAwm5YcepityDNZkk+3ZzbOb7oUK6uBkZ+7Z3JqdYPY757zmkIxz2x+SOwCoTFUKCgoKVyl0HUk5x/U5VDkU41BQUFCYQ6i06goKCgoK04Iq5KSgoKCgMG3M5ZmMZLAgGIff4vEcqNPdnJOTc+RxdRlSCpf7uj/nnOfG37lO7kjMvESf7I7m0XbZrrVy5zIAhDtpmSncTp/Kn0inHcSDS2j1v+gwfbI88yLZTCLIOL9bfq2IbK/c5c9BO5lBO9cnC2gHNUAHLlC1UALdA2TfiULaeZ3KOeYZ+M0kwd7fp/PdDT3JXFXKVKWgoKBwFSOZXFVzCcU4FBQUFOYQuq5hIpZMVJXSOBQUFBSuSizEk+MLKyWjgoKCwi8ZzKiq6f7zyzg0TXvU+Pe47drjZhvVVzEOBQUFhTmFcI5P95+fk+Oapm0HsFvX9ScB1Bi/A8CjmqY1AWim+s8LU5UemySP+XM1LbgoBy7qKlhVKW1LPdVO9mXrZayk03ZwUVP919O1CYKjdOQTV5Mit5VO60FhIkzbZSeZT5N5kU6LwaXVGC6j5Z4A/eosBpbKo8Iu3yBfM4kg9wwtL164hV7T4Q7uCfTHz24eJNu5dDCjeXLSkROh55WLOIsVimefWQ40r4gnjjX7Yqjb30veg7y/z5QiMx21pevJheP6PDleY/x7EoJJ1BjXP6/r+m6u87xgHAoKCgpu1J0C6k55UMeTyTON+YpkfBx+YGgaJhoAvGT+rGkaADTouv6ErP+CYRyn16WjaV28FFXbGEXtrjkYkIKCwhXBgS3i/01753YcswUdSUZVCWazTtO0Pezf6vo2r+uapjUA2KXr+iHj754wru/QNG27TPtYMIxjWeMoljWO4vS69KnfTdDGGgUFhfmIM/UpaKqPJ5g1n+ios1WHHcwzpfG5jCOaPehI0lQ1M4/fbmMW9wGAruuvAOiGZb6Kw4JhHCZMrcPOOBQUrkZcyI7iUk68I6dsIA3VffSp/fmAumOTqDs2iTfvF2O94wfiXdwp5Te++8vJMCxoSZqqNABolGkTbG9Ne9TGNLZD+DpMp3gtgO/I+s4I49A0rcFUdZLqn5qGYAFT98JAbaO8DoAMlPObA+dY5xzzXNgaV5uAc2KOMakjRspoJynnXKfQX0cv9iy6ZAXpfAaA9D56bEFmKfTX0e3BKP383mvkz888588xP5lOf7vSD3jPfg4CWIkQ9twkPsS2fdZcd62lxzcRpp3fOS102g8qqINb0+6UI2mjE47rbud5rmsoftkIu2cZ5/dMYy6c4wajeFzTtMcAFEA4xQ8Z4bk9AJoomu6bcZgDALDe770SgdI0FBScSPMZPTbX2PKLBWf4mHFc6SSHhu8i3+P6kx5/HgffM6br+m6DQyn4QKL23vmIrlgUPfAwmWSnYfEgE5N7lUNmbqqaTEF1G52o0MTmg7/c3/iMsS/qjk06r2/OR/NNcbQPtQeHsOxDf4kMrySEc1wlOZxV9BeKD5zb/cvlEjftvQe2iynZtFuo71wJ0fmAokAIRR7nBrJoK5sCgMWDISweDOFImfhYay6Jsys5bQtcjZhBNEkYR93+Xs/zHDNhajq9IRNN12fFXZ8NprQQU47MKuNIJEwMwLrp3PPAPWJj3flUfxIjkkMW7lu9ZxS177CnrWYMAwXz7zB/66IxnFsUf1CwIJaGooD1zVpjQwCAqkD8hqMgk7orUlNR2ZmeVN+ygTRkMofguP5hXDkHs8kw5hOaqyfRUhNPnqqbNaxg/FfTwQfXi7V1w0FvclTrYhhXAss+HMayD4fRUSXWQEnrLDLyJH0cc8k5FpzGkdM1O4vIDPd1Y6L1yjENANi4K/7UrfQE7Sc6lnT5f6aMMSy5GETVxbSpf3sbhNNyyyHhXLy8wbl8RiWB0RRxNiVuL5NWel8MbcWjOF8S/00qOlJRMh6e6nuiWDCtVZ0W0+ovBXrGo+idiH92fjANuQiR/XsXSKD3h9eKdbvhY5rJThc1LSmoaQHevk2sjdvenh1yMUj76eM0DTcOPCyyK2x61sryQGkMtXsS1xhmlWEYmONw3KQwq4wjkTAxTdP2YHJyK3VMP2aLbLrhn8T/00qUwdR50qoqpn52m4q4IlJ6YbyN1Y4YY2qKlDvVavOvR4wsK+H2YfkJWgCTGfSXmGCe37c8iEwEsWociEBs0DCEaSA1Deivsv62rk+Mtb9KtI9WOjdV2aiQzkbTxfW0fvF7EUIo8mAMw+sj6IU8eifzUBh5CCNv0CNKJwyM2S5H0gSRH7MRoVhFBHkA8uD1DXSk2L6p2T9l1Lq2aJVcaOjtKJO2AUD2OZrpdGygNcuhSlrbmbTxiKGIiKpq22714VKScOleJkMWaVjfGDCu2cbMFAij4I70W9M0aVwX6yqXK15m268AMFgm3ntqH49EsOzYJJYd6/c2bTOmrEQKwzVtLUHLtpK469V7OrDsINvdBS1J5/gC9nEYh0Y2aJp2n3FwZEFDZiqSSjCGb4Jzbsu0hiUXx1B1cXZNIu35o7iQHy85Le4VUrcJk2HIkDFJt6en+zOz9XYLApJfGP+cMxlCI6gb8TaD1Q77s2uvPu8donnuXCrOn49nPEWBKIpjM+eUbouK96sMTc/MBwAV6cm9e0vlOFqXxEvzVeecjvns4eTmVWbqqugYdZggs6L0uuL8ml5augm/vtDzDUIwrDjk9KXUvtOB2nc6vNunWSHwqtQ4DGax4BmGCdkiNG2e/UXik+V2Gc5rQ+MwndtumM5tU2t46x6xQG5/TUx7pNzJNFoXCQLvxUzOGFVmpxtlVd6bjvJebzPG0PTpVNLgiH9fj9jkXowjmkITgDCTsqH9nJjX8iXeGljmmHf/JUvGsWTJOIaGxLxlZYl5631HqDadgSi6Ujz8M8XCPyMzteWPp6Eg1WI8Y7r3+8ki1uz90wPTT1cBANVtqahuS/U8BzITME1dboyUTc+kxvk1c3tnj4SevEeYwdyMI9H2hKADk0mUjlU+jmng3T9YAQC45e8+mZX7c4vQZBjJouYT+v6mr8GLcZgai9tsRTkx7ZWt2/OFLVzGRDoDgoDIJOlLYdFeFvFu7+sTBDIvz5s4c8Q/jwgM8KtRjPk8/mMyDDeKYyHP75XXKd61sjMdlZ3pOFAvmOamY4Jpuk1VMq3BjFjrikWnfgeAyQSC7S6nRdGZHs90ikfTUGQ7mrq+kdZ4WyrFvCYaHuxGc7UQqGpavBlcW7GYHFkghB+/pleKoumg/CP6pAHXngj0JE1VKhx3GhjLnv/hqRQ4bWHJRfmUyJiOTLIDgAFbVnfTXCVjHKbkLGMcHRmiXcY4+vvFBpcxDo74e2kaJjiN4nKaIKylY95jW1xJL3WOqXLgmG5FB71uOa3B1Dq8wp57xsWz7RoMIL6F7HvApsVwpijTnJUs4zCFGtkaNTUyGePY/PoQeX/ZOQ/Af4qiVa9f8NWeKOayDGwymBeMg6vHYU8Zsu2lgbhrE63+YgMnT1jUvGmrcHiZIbgpq5Y7/ta9SEfLnSEh7giX9PYB8tnRPCfBKI2Ejevi97CtHIgX0xmoptM72NOCFMTS4q6FbH7AkrBoTzX81TGXIFpsOL/N60t+5CQ4+iKxnJYcEtftjnUAyHT5UGIS85CJ4QY67QW6LELTadQdySqzBh0bsH42SV7MJoD31Fvf4UJENIQXW4R2oskKjOgbEu15WeKe6Q3OQypdRwRxLFgniGBHjjO8Nh1injqMbxIL0dpX7hnnty3LFs81z8bYgwB6s02BwBo755wfXkQzAXvKkepm8Z1Sopa2PRkKSn0k7kOrNYZxIGVEfJvsZqfprsagQmZ6HS5liTbiXBdN9bkAgGUH4xmMZ4qicIjMtr3sAPl4FhMdnaTz3I2r0sdxpRGKzO7nMidbdnZDdhjJxJDEpGGiP18skNmwy3JSs/3MhRdkmoQJufQq4NfJP2rkzfJyslOOcwDIC/nTRPOD9Nj7hw0zXJb33xXE04hpoTci5i4/7D131Al8k6HLMBQS307mhObWjczEZPpI3EjrcRJ2Tsv2mxmByl8n0zRk4fduRDPEWgyNeDNiWTvpPHdDT1LjUD4OOrxtxVn/95elJ6jZ14ua3Za24yUR2MEdRuLsxe/tEBvNzATqxnCauL+Xs3ZACFbIkZx95ExRcw3OR3LpkhEYUBXfTjnOATnBnerv0hjccJt53MjNpBlTUSm98UdjYl5lJqm+qGBM3Ht4gWPoR2oFITf9K27M93Wz5x7x7be9JglcmcX8dXu+KLTOO//5clLtiTnNr8Jw3JmCyaE9YTNLHd0k1NhrDjBmDBdk6QkAZz0P7pQ4dxiJsxfn9NAmhOMVgold3xx/kvi9W5wRWW4s7qWZVlQXYw9p3sRrJEW0y8JuIwHjnIfE3zAYFu+WHfH+BpyPJC1NvhEox3ki4DQGDsn2M3F+VJxTqpVkpeU0ppFUY27Gpx9BlSmZDxPcuhnIFustZ9B7fgYzjXmXrH1O4OHaR8PJE0gulJdrTx+m9zvXngh0AJPKVOUPXNTU+RVikU+XcZg4druIM6p/awaOXCcB82ChDBmj8k2e3UcvFU5iPKcL4rVc8yZeZ/JE+9pu7/amTNF+jeSob+MqoVGYJ8vdKBmhCdSiRfLxU45zABidMCT6oPffcRoDpxFwiBom1JCEyKVpNPHmNI2TJeLbN7THf3uOoa9tps8VcOvm4A3i/rKT4x+tEwKBLJyXE3i49q2v0qe3KeLPhfJy7Z96hqYTXHvncnH/4lN04jblHPcJLmqqfr+/QvPt14qNJ2McA4vE4s+56G035XwUfkMX69vlkUd+S2emM9VBwhN0e2iSbs8aoRc/50PxgwuDQpCozvc2x3AaA6cRjI4bjCXVmzifOyP+X77G+/7JHOyzIzwm//YcQ/eLbDq+A1lD9LxzAg/XHmLqrlDEnwvlna0URiY+/nURDbH9L45K/2Yu6nH4xbxgHFogBYGwkIpu+XYrAGexFXvU1KJW45qtP5cWxI76PWIX2DNo2n/+4NFSAJbNMtZ63tH/vftFNcXb/0YUykqHM/1B602CQK18TxAyLrttwYd0eoWRujyyfajSWnDulCFueCUftBdDWhV1tme2OrWj66dCQcX1c3c5F3shxJydM0OAQ7SkuKL6Itke+bvyqZ+9zGB3/0+rHPI/7LoGAPC17e9NXftx+1ry/llplm2840NBfGo3tE5da3rfCgtrjTgZS6jCaW8PZ4hvHkoT32blFvqc0fF/X+H43a01RFzuvqXIclxP3WidHwifEH0C661rHccLyOcXHLOojpdfrXeltf+WDZrXrP75Jy0B7sb3zX7i3ScznGv+ho9gXPceS1x7Br1nUl1RVW7ib4/K2vSzybhrum1Pb/ye2H8zmZXMHolZfCEWd80Lc8kEksG8YBx2hGbAZkih8gRt4srppHP4ZF+iHXGc8/ydHeL/rbvIP0sKZ1MFcVs1PjuSp19QUVMA8M5+YezeujlecuTMYF/fIZfoAMSd/HZj0wbalMCZmpat8rdu/WgN3LPPRsS7LQ17Z+Cl/GrzASdWivdbddJbIOLOefiBVwLF6bQ3/CKRA8NakqYq5RxPGDKb4Uzlz9/8Cn0SVLZATHDO89GQOdneBIyKIHn/GsH0bjzqHeceSiZtwRUEFTUFAKOEOYYzg3E4fFh8s82bkzN1+jU1vfqOUBXu3eodfMGZAf1gknGjUn61+YAL5WL8q05e+WcPMulRuPZEoOvAZDKFnK52U5UdnPNaZjM0c0lxcddXJL8+gVveomebiiAZS6P7Vk/QxK0pJphtbcBbsnQXE3KDY1znOwVRrij21gqoqCkAuGVTn7TtupP0oTAOmZlzmyI9Oko73WX5u0ycDoq5WTYxfa2gimF6lF8N4NfF3puFVC0rAeu3feWJ5BkbF8r7zm8vAQBs/edznu0bnznveT3R9kShTFU+wTmviz6hPXUyhmHi0KfFWQ5Z3LVfcM7zEBNyTkWQ3HDEn3OZkzzHU+h2jnFNxuh2KmoKAELpye+ev3qtAQDwJ/cc8my/9lraw0qZyRLB8cOCMaxe661x3nOLv/U2ocm/DffsIGNm48CuC0bo9ttefoEe/8/vFwztUz+INzdyobyjWTQJzLlMC5hc+5v3C0FVdm4LELaHZExVKhw3JWXKQe3lvLanI1n3/XjJgKrlAQBYaSVsKrokPneswvI+2p15bukn2O00bbglGHfKkfemmWk0vd0pCbkjSFIiFiPMmHLPWNey2ujlE7SdtF8fFJpA2oTVf/gRS8pfZaztCcPJGn3ZeWDy2lYxJ9E8YyOHnGMvLzcoQNCIQGqjI5lOD1SS7YGbaILx5Ou3T/08GB2Iu3brrY1k/7c+tEKgRsfEGvrkhBXscKPNwe02Nb3/iS0JGIAJI19Uf5eQ7o+8Rzunx5iT5mku+ahu3Kk1DNuc3xPjA3HXMhlB2J5uxQu5Z6yfr7kYr7GcecCi9uVGxsUzD4j5KjngnLdrW8V6G6gW14Ou7A/rPxHtI2WiPb2P9gukRpybZDTDOAhiXNc7rYCTW74tmKou8Z2awTh2sKVnw/TenjhhPwovgjYmT1DH45WPwzc45zUHztTVwBRZ4aQfToLhQhM/qRXEdkXTzCdrbFwinITrznmbJtKYcNtU5oxbGmOHDwb9LeSuYSMDbGb8xuQcvJw5hkN5EW0K40xNU0wzSXDvl0qEUnPvfj5T7KmKYe935J6dxtjfgyl0O7duuHXZlS+IflGv9xxse1peeIkLtuHaT28QTETmJ+Wy71JhuHYsMEvV/GMcfsGZujjcvJcmENxhJK5858UysVBljOPYevH8+o/iFzTHdMaD83v5cT6WwTEhaRZ5CHycmY0zxzy/W9iyH9zubcvmiB9nauKY5sUUQbwXTXoTb+79yGcz795jaIYyxuHn2VcCR1eL8W/b5703OfO0H5gBNzLG4Tf7LgCVq2omwDmvd/+ZUP1knHz1m3zZRwrpY/QEcoeROCw/Q3/y87VyxsExHdOUJEPTIjH42oveqvb5VkGAKqq8N+LZItF/aZd3/+5uscELC73HxxGowgw50/WrUYyM+lvqYSabLYe+gPg2MsbBvR/HeCiUDzHJKZlnt+aJZ1f1eT+7e0AQzcIc7/nj1g23Lgt7kvfRcBYIrr3WI+Ouo51IsJgodAB67Mr7ODRNe9T4sVbX9ceMa/cB6APQoOv6E7K+845x+HVeVxym4/HdFfiuNBZfpj/56g/l9l2O6XAmgY4Cce9aybm7ni6acXTmCOInIwBDQ4KxyRhHjUYTqJx0ufmOk6o7xgRxK0nzJm4P3BZvy7aDI34cOKZZNkETb+79KMbDvXshkz2Xe3Z3pni2jHEMRcS6kn07bt1w63LNCXr8R7cKDfaad+L3PmeB4Nq5UH5O02j8NaHpevlm7bjSUVWapm0HsFvX9WZN0142fu8RY9F3a5pWo2lag67rntEm84NxBAJTDqfiNsPhanNA2U9dTkUnMCcx7dDO22PnxcnwgO1awPYsdxJFd62PE3cvFo83CrikFK5ztLud6yNlNMG4fEOR5/Vz14r/7c7vXKOqxPAi6+9Gc5moEVu7KXkOl1mEIviq5VBdYiQhDL4qNuqQy3dtEj/ztLrb+V1iPMq8XrTxEjk2+8ltL5ztph3MX1lpFU547OWNAIA/udvKy/Ltd28j++dWWBFUrR+J+apZbUmQ54dyp35uahLvVFsrvlFukVMSbW01+q8Q73RXwwny2T955iay3X1yvNj49mMetHtwUuyZEliNUSaZQgpjWck5a2n8yyPi3ewRROl9FulIyRcMo/SYWAATrjFWuQo0BUedgklda5rjeug0LTQO3+gMTDi/WgSwVEfFXgq/bVkjvCwQWpUVAGEKavZrGKH9rFz9Hzu96jJy65Enx+fGVFVj/HsSQLPx8w4A5tHkZgDbAcxjxmHD+rf95aLqWCQmoOSi91flNBkuiWL7ekHMZJW/OOd6T5ZYqAVDM//pOScoJ3kWMUkI83W6PU/zl0G2o0swxpKi+Jj7y+fFfJZWeG+wz61v9vXsco+CRI7nG5qiyTim259De7qYu/JR77nLIeqFFDMp4Yd08T2zNG9t6HJMPLs04P1sTksu7aUDPUoGmeSW3f72gsl4vMBZICqbaRMklz2Xy2133V46E4WABlzhqCpd15+0/doA4CUA62FoHQYKZf3nHePwi//cQte74MAlUVz5Gn1ynHOuf1IlFpisNsKgZlRJ0+M3Y1+GYDp5I97TxjlB5xpuqd2NxqPim9y+LT79fb9xqbQirgkAcGMNnQ6fM0UVFtPiW00NLaJz/c9cFOaUukXehKw3TcydjHFQoP7R4l4AACAASURBVJgKAFzQBWNYLmEc/RDPLsX8XDdvfElofrIMtn4ZDwUue+4Hj9YBkPtcZQKsGz5MVes0TdvD31/f5nVd07QGALt0XT+kaYkzogXHOA7dLIYsywFjJhVLFpWnaAmBK8zCOdfzB2nGcj4oNvmq8fhNfrpMtMlyCnFO0IFU8W45HvcGgD6D8eRFvdsppgYkINkyUntRoZzZl5RLmxICZ4fnUFbmT6P49/eWAgAe+9Ujnu2Lo/TcDRvmqMyU6YdxZzLbvESjn82Fw3JaNCfwdOcac9M/8+Soo1acXSpp8hYIOQsFlz03+4K/4wNTmLsIqe02J3gfANM+nAdAmoF13jEOTrroXEw78rikYh/dJbLNrt8pT28xm1h5jpbqsmLyKckdppkOZ4o6myMW+dpub+LTXCjaG9q92ymmBvCSLSe1N6yROyLzCmiGfPyCmNfVi73ntSCbKVbUJ+6fk5fcDub615bRGQ8KxunxXTIc4LUe2ZY5plIeoOtxcCZGLhyW06I5gee4IUhsOeS99mW0wATFeBo/VwbAymbtBmeh4BIo3vjdJrK9rUbQK8okputJRlWJpdYo0yY4aJr2qMk0DOf4SwA2GM01AHbL+s47xiHDQx3v4qGuXwB/7rz+XNHNeK7kloTv07mUljg7KsQnKTnvzYC4wixcJk8OlbLc0wCWX6YJAIfsMSa9QoRup5gawEu2fqV2Ck/vEzm/H//8e57t2Uyq7rNNYr7Wrvee954e0V5Q4P0OXP/7NtFRXRwyAvJvSzGVmQAXDstp0ZzAU9DvLyUKxXiKztBRUX4tFByObxBjqmymTedXupCTwSge1zTtMQgt4/OGuWqD0dYni6gC5gvjGB9HzEgTcPvfxOfHj42M4GmU4WntPuzSXwEA7NDuE43dALrpavdBW+nZht3GQrJFUsVsKQoOfclZb8MdDfGxK/fMSL5TWrtQLlTi2suJOYq5lCFc1JQdXqaoYZtPoMSolzFsBQshbHMNLJkQ7WNGFpWROudiLzSWywjEdXc9jaWu8fT+2xJyvKfvoGPgU0MWER4aEN8pK8f6Ht/+YNvUz+HwWNy18lr6TM+WUktaHG4XyTPvWmIR+J3nVk39fPKkeHeTMfzzmucd9/rW5VsBAE+s+RkA4MsffYl8NjY7BQ/3+5VmOwle6dRPQnO5/J9lU1copiLDpE1+8tJYelZa67fcCIrosdXjKDhprY1r+8x+4tt0r3Gu/VKI80X9hsKRs8u5rq770GQcxvUIvS7saXgAi7GZ1+0RUutN0me7Zk//ca0R/GYXB7j6PvZaQV7Qui1zdsXxbOMa7aS/0qYqXdd3A8j3uP6kx5/HYX4wjhkEd0BQpkmY8CuhUBEeAHCgXqi+MrXeDzhT1FyDs3W3nDaqJy6LH/8Fg57LKuyVlNDfndMYOI0gO5ee9ye2/4xs58C9H4VF6TQhOxUTzGZ5wLvWx2xrLH7BmZe5cx5+8PNHRIivrEQsl13X62yJN1Tp2FlBMDiJOxe3oasjjIJIBC+EfoIXSurxxoVKTEwkZxbyQsOP6HBdzofiN8LjYI1YaF724MOFggDIaoJzpqi5BmfrppDps8bQyZNCxE62Hkd1nT+TxqkjQqRcvsabQPh9Pz9IRmO5EvhC9CCq8jqx5qUOFEQi6AmHcaSkBE0tNXi++o4rMobRTKbWPZNdNyHoSE7jUClHLBz6VaGU2wn4p5eexR19Z1Hd3ofwhCDcxZEIHm1vxI7ss3gzbyl+enYpAD6pWNtyIVVx0VOzhdnQNExUD9KSZ1NEMB5ZXe2jRnFpWRW61lZhQpAVYnprj9B8vcJpAd7W7aVpmChf6k8iy8+f3cqSfsG93/79Ym6TYXwyTcMEp7FwAsseJiM0p/G8fZvY07e9bZGjO/IacUOnc88XRSLY2N6OktwIBrJieK34LnLcAJ/WnLNQUAkUAe/sunZw9YEENCAJ57jKjmtDV118vqV1sU5UD1oLyER4YgLVg31Yl9OJn8ZZ2L1xbLPYJJWnkqu7wMFvaCEljcs27kKBX+c+BY6prVpFR3Q9/iNhI5KFy3K46akvAgD2felpz3aZpjET4ASChYhroxdQPeC952v6e3Bd1lm8ZrgiOMblB1wCRS677p4vikFyB49VIackoMcmp2puXPui4OD2GhxrOjriFpCJ8MQE1nR00A4tm7Ot4nhq3DU73vpjp3PcXW/DvUgzzjjtrnunJBzhy+i8yTuliAkmhRGKG2kJ014XwQvpnZaUPyXx2WjopO35Zq1y85rb+b3CmenBUbsCAKpKzeuC+afX0cRy55Z/JNvv//grZPv4QLyQERiwlnRXDp30sSUrfm5ahq1rdy2x0oa4GctD+37TPRrH9dgAbXcP5NBRNkNpzv6mU37ICLn2ShnCpRGxYyKHJnjjOda6MdeFvZT9wFJrfA3tacY18XvuGScVvB6mMORNHe2aBgD031KLNa/zez490wxxFs9Pbxe/253f2/9C/G9/Wy6FUazVf1W/B4cO4qHhD+OjQDM3xP+xMlX5h1eYa0GEPmSTH40CxqLmHGmJO6uSA+c8P1MgGEFdz8xL362XhWO/qpQpRDNP8bVddwIA/mXHG3Fthz9yRjW5UZ3vzwR4/+2SDHsJgjKzJQIqMIADZ4pqj4k1JzvPMR/XTUJ73kCyWSISQbLZdV/KuQ49i1LxQHsjCiIRdIcz8GL5tdg1fA0e8rrRFQ7H9Yt5xzi80BMOo4hYSL0hS2zmzmlwkB0UMsGpw5zzfCBMt7dkGwnbPPwVHxulQa89PHPBAFcSnK17T1vV7D2bcU5zSNaEZaKjVxC6kvyZN6dwGAa95uYjprPnZxPJZNe9I68Rt3eeRnWn3Sc7gq+ePYjtOafx/JoD+x48csSR5ZKoDDwvMe8Yx/kG4WC1p/Y4UlKCje3tnqprJBjEkZISwBBAGn5KpwRJzFk1e6jppk+OD6bJN3kXkw+Jkxg5ybMtRbTLDiEeOiLuLzvh7Zc4fnt7vKZhQqZpLBRExmjTEKdptJwR69Ye3RWKRvHgvt348sH3UDgygu6MDHzv+o14/qbtiNoI62KNXnPcuuHWBadFnyoV7W4f1wNtu1Ctn8WaDmfU1OmeusT2fALg0ppzKYy4+j5e7ZR/pnqgD5ezMzfHdVKMwx9O3iOSEtkZx8ehxSjNGY6bjEgwiJacPHwcWjzFOGQFoEwk6qyaLcjyQKWNRfGFI7vxxUaLCDy9biPeDmzGuJFLaO3H/k7YcpLnUMBol9C4rm7Tru1kHKGRKL6082f4yv73p8b+r5tvxIs37HAQMM6kcuuS5E9XXxoS0mlZljeR5DSNvYeEwLKlgRY8pM+/IL5d2WLvLVWc508THpw6XW2kHo9G8dJ3/hZV/f2WVDsygt/b9y4+c/woXl21Gk9tvxeAPHdYouDWBadF92fGd/xM/y5sHDyO6v7e+KipvAhas7LRnFeAmr6e+D2fm4+Ps5bgoctGNgkXnsEqPKvVAwC6VtBrjkthxGXX9WpPxCfrgI7koqqUjwN4WD+GR3AC+PNXHNefwSq82LcNb6dfgx1Lj+KB9kbkR6PoDoXxfcNmONHHv4YZE/7s407pprWvGC9iNR4aa8TD44fj+j3VdBuert0x9fuRVYIxyQ4dceF/Xnjgo1fx2U+OxBGB3/3gXdyb04g/Wv0V3DbwAR748CAKIyPoDmfghWXX442cLQASJ0ic5FnhLqTgwrpr4vP2PPrKT3DvkeOo7elxjP0be/bino+POQiYDJG+AC78dQZueO08CkYi6MkI44N7KpD2K2MYy0rscFdk3F+47UVJkSETrxwQZjTZQcHIML2LM9L9bbWltc73e3Dfbsd6MRGemEBNbw82trbgKV9PBEJjUfzae7vx5Y9sGs36jfj+xu0ArO/FadHLLsW3r+s/62AajvH39aAjM4x9pWvQkXkGazo6kB+NojcUEns2Wok39C0Ilo6hp1jDA2c/tPwIS9dj55kiPDwupycvwHKEJ5b2fHqYjn9mClebxpFoqUHyHukhvLjoTryIO7Hz7N8DAO5a+s2p9lhnN8YA/GRgGb4eEcV7Hg7cB1wEgFHERiwp0cvUdde1nbi+7Rxq2ns8Y8L7FqXjB2Nb0BvMxgNth6xFWHkddpZsQkrUWtzdBULiM68N1LuLDQ05rpsFcVJjUXy6ew8ePPWB2IThDDy/4ga8dN3tuKG9RUoEKvv78a+H/glZY2MOe+mjx/fhjtyjeO2pNXh5xT0AvE0CXest00Y6hG9k3JbQZRxCen3krbfxlQM2jWHTjXgm5zZEQyFpe+uni7DpTIuDadjHbhKw79Z6b6R7D34NX/ru27j3P0/gU92XrbkZieBT/96C6j3/hr2bl6L5z6s8Cff/2HAYQ30BfPQXGta82jrFdI7cW4X1f6bjf4xd4/lcBwZiWPb3bfjHN46Ld3sijFfvrMfpb1Zi7+VlU3/WdEl8072XawEA91/jTONzptCZNv0HRxv4Z0Noa4/8+1585d0PrG97yw148YFNiGZaxDnNmNKoIY98+eB7rFS79Bbxrc5fEMy3YrElzHzSYlUD6+gQ1+2n77/yxuv4zNmj8RrNgXdx7/GP8eMVa/DCekOjMQqMmRGC+Sed823J/IJIpxw6hTUxXip/QvscoK1F1DjjEooI6hpdW4p7OnfijrYTqOnvcfoRTh/A7bm5+K9lD6IvVhTHVN4K3Aity6INpR75X+00nM2u69GeiH8mLp7vamIcRi73hEoNXil4mbqu6z3rWGAmzJjwa7PO4/bBU6ge6HUuwub3sb37FN4sX43CvgAeOfdzYK/jFvjX1Tvwr/XWKVavA373dO7Ep88dRdVAv+P+v3X0F7iz7QiKh4flm2hyEukjI3Ar1OGJCVT192PDpRa8vEJc8zIJAARjuP02PPTGbnzuyBFU9TkJxDfe2Ytf+fgwmvLyUNvXJ20vHiLG7lLLvcw5G0+0oqbbW/Jc2t2HS4c6kfIXwIE3dhpjtwj7rr8KomZvBzZ29TuYzsZXzuD8nlzUbGpF858LhuPFeGr+vBVb9p/F0m6bE3MkggdebcTZfWdRtKETP/qGMEevXEnHusrqbMgYwzP/ZQseevld3Nt4HLU9vc5vu/sXuOc/j+HV61bhqa96VzEsHKHDtPOjUQQHx3DHdz/Gr+9unHr2i9vX4c2vXuv420gk3t/X0HVWKsxU9ffjhvYWvLCeHAJS9SjuHNiLB047NeW3guUo6E9cKt/zBcF67FlyqT1d1deH70T+D7JdwtZXTx/A9tyTeKuoDm/2OSt3ysBl1/VqT8Q/8yn7RR3JRVUtYFPVF5BgqUF2IIEx7Ag2oissUoo8c/k7eLHyOuyaWIcxIM6U9ObwMwCAZ1PX4hmb6ln+UY/71qzN8Yb2dgSMn91t1f29WJd1Fk8seRi9eTE8ePqDKQnmhWXX40d1W9l3u673rINp2O9f1d+PkGRsJmRWWDdh9jIJ/OZ/vI57j3mYkt7Zi7uPHENnOOxgCo6x9fUjEkgh27mx50ejpFTNzc2Wc+ccm9BO2C+FMlHR5T22iq5+LDvUieaBSiz7+zZPxrPsUKeDadj7L+3uw/qjbfjJ0Bju/d6HeOjtQ1Njf+62Blz+VjGQE5jSWO594xgKRyJT9w/dtZpkDHcfOobOcIajzf78mu5ebDzRiu8PR/Hgi/vx5XcOWuairdfzUm16Ov7mN19Bdbfz2V/96QfY/v5p/Mfa1XjmVz6FL/14D77yC9u83HwDnvrstmnZ6bvThcZiT+v/mf5duOt8vLD06PF9uCM7F72hEAq9TDbm+G2+sXSPQBZyfISwVd3fi2uzLqBoaFycs3DhucwNeNZuuv5zZ/uzqWvxXJrFdLxy23E+2QvZ2fvdfRZaVJWm+ziyqGnadwB8x0jHux3ADl3XH5vmPfbUhnO2vhpOd0j8gPmh8/FmQQ3euCjEavcBPcB5WNALb+qvSIkvIBg3xe+7QiF0ZGbH2WQjwSCa8wqwc+lq/LD6bs++6X0TePkXf0Fu8hjkzIHDpKZh66//FT535m08fNSyRT97zUb8sO42/OW+J0npJwYgkyD+3Ni49q5QCJdysuMIZCQYRHNhPlZd7kxqbhIZe1cohM7crDjmEAkGcbYwDyX9QyTx6gqFcDk3O04jMvu35WWhsm/I8/5NBfnoDGf4+vZdoRC6sjLjGHckGMRQaqowX07Ga5mRYBAnCwqw0sOEaLa/V16O4kgkzswoxl6A1R0d5LxMaho2PvrX0oCOb7U9S757IuN7Qvuc57Ojy0rx9q7/O+k90xUO48t5D2FH5lHPcxZjA5Z2aQqod2Q+4riHzCf6bOpavFBwE4KBian7u32yr114cmpJa5q2J72uZuuir//OtN/j4j98G6Nnmt9Jth6HH8yqczyRkoYA1mXFRlE9MCIJX+vFuuxLeAOCcWx8xvtU55Rz3YVnsIqVzjgURKPInJiQOvLWd7YgJ7ITy4bOxIUWNqXWsM4yQGwWr03EMbXe9HR8563/1xF9Ujwygt8+9C7uaD6KsqEBUnLkxIZEFGjZ2CPBINpycjwJhClVR4JBknjKnp/I2Kl5W9rdx2pLXP+oHsDSHm+NpbanF9WI1yZmavzQdXG6fHw8jvC35uZiyQA975SWXdvTw85Lb3o6nn0lPqrrdz94F3fnfIziEdqEWTkwgNbsXFQNxjPF1pxcfByumIqU9IKfPZ0fieDx8I+k5yzeXFQ7JajK8FzaOjyXts6TsQQATMSC2Dm4bson+8WcLwKSAC0NyWkcc3lk0C/jSLjUIIXM8XGEY97nKqbUYkO8yLkcH60UaKjH86jH8wDe/EjkGLhj/Z9NtVd2PeVL8jPHIR3f5csoHRqRhhZyanlPKISujMw4c9aUZOkiDvZ2kjD39bDEkQPHuHoIqfhcHk/AhoPBpJlmIqCeHUvg/lT/awmpPBHGkAgoc8ywruPbG27Bw4ffQ2Ekgu5wGM+u3YhXVm3Hvu/9N/K+GRMTJFOm5oVbd1UDCZgwR0fxa9v/G+6+8HM8cPL9qfG/sPJG7B5djwk9DZAkV0jriZB+BG7dRINB8pzFuuxL+OlIpaONs2rY2wPFhXiw9wAe6reKiu28/M8AgOdyN8Z31rUkfRwLN8khWWowERVK07Q9qbEY6SjIj0ali8hEEGO4I7pvykfy3PG/wYs1G/Bm6CZ8HK5Aac6I1BR2OZs2KfiRTGv6enCyoAAZHu3m/Y+UluJ/1n0Zn+7egwdOfWBtouU3YDwYwz0tx+IclYlKlhxxHAkGPSVP8xmcSeFIaSl+/y8fwIMv7seX3jk4Nfantl6P5j8qx7dve4F4OhCamMDpsiJUdcWbe1InJxEkTKnc2BNhmhRx5Pr73bZ+x58/OooX1tyNF9bEm0n9atmhiQmcLiyctXXXGwphLBDCDyvuwg8rnFlu3fnfDmwXZGrTbut5jdlVKM0d8TQfc8KWOUbZ2Nd0dPgu4/B8/iY8n7/Js+0hr4sLzMfhi3FMp9QghfFAAJBoHIDTUXZmswi3rdtvRU3dOfpz3HHxhGMRFUdG8NVP9mN77gm8Vbwcj+Xdjx1FjXig7ZBlczSc77dmHJEuwpbcfJQOD5IaA0AvxMqBAbTm5npvwpxc/Gf+UowHQvhx8Z345q9sA2ClNhnNC+InS27HZ8+9jQePW5LZ86tvxI+X3Iaf/fBPyHGZz5Ftog/Ky1EcGfH0QTQVFOBCbjZSAU9beHNhPt5bVYVoZgj/+pu34jtfFLEiwVRBMrbmNPFO3HAY//GjBiz+X524c+cnKIyIkNqdd66APqbh/p2HybGXREc8fRDNhfko66fnrScUQk9eBiq7BuL6txXloLBvhOzPSbYcY/igvByLh4ewxENba83LRfHQcMIOZDe4yB5Oy+4Nh/F/bfqGp0bw+uJP4a3X/1Ta1/4cqcBhO/k9mir2fvq4t/42UBB//Y30T2H30ptwR3QfHmg+KPZ0OIzvV2/AuBbDXRc/ke7nld30afD8SAR/m7+XLOPw1vlKKWNJClcT4wASLzVIYTg1FZHJSXqRGalgmm+KZxzrBlulh4lEFMV5FPVNOlTHksgIvn5qHwpyJ/Fs/h3SRfhm6Cb8UdcLvjSS/NFR3PPrf4r7TuyOMyu8mXkLxgO2eH2Pc4NjKSG8XP0ZPHbHNgDAhlNWigiOMPeEQujKzJRKjgdWLMVzn78Fj/z7XvzGu5a282+33ICn7tmGaEZInAz/8R78xi9s7TffgBcfcp41aD4p/rdXsUskNDGQE8Cl/16Kp/57qaO9pyeMtv1nUdkfT9jP5eXio2sq8ZOvbMC93/sQD/7s0FRU0/O3NuDVL2/An37zR/SzS0tR+9MY9v5lIa57rR0FIxH0ZoRw6J4KrPx/ojj86WxfDt4PystRGh1GtcQ5f7qhGG/+QYOIynrz+NQ5lFfvWI2/uemz+Od/+D8Jpd0YCwjCmxazCOyB5UtJ53dnOMze29QIvnmvMAhsOWSp/X7X3aGipVPXPlgTjbu/HRt3eR/S+7UL+/HwxXenfi8ZGcHXj72LZxdtxR8v/U3pfn565H+TY48Gg2wZhzuy5PWB3gp14M3oDQCA/kIxJ7nddIqjhRZVNS9Ojg+lhNCSlSaVEBpTrfrKNfviU0IkEjr4ROnn8Hz+JrzxJVFw2x4TDgATSMNPQp/C1yNiIT686o+n2ii1OBGNpDcUwlgw5GlWMA8Imth8UC5F2hmGCZYwl5biv936W55M65VV23F5u1ixTz68A39y0zYAwJIacaArNiZU8mhGCP/y63fiT7aI9ooKMcb0TOfZhhSP1XTiujKURYc9I49ai/LQ0uA+QGlDTgD/9h9bPQnr6W9WYm9EHNB75fc24w/WbwfgLHR0YMVSQpvKx8j1YWTmDWL9346gaOW3AABdX/0u1kPc44P6paRGcyEnG2ma7tneUpiHo9cuxnd/ZzVu+pdjuH/X4anx/2DHWrT94SIRzgvg9J9V4W//zJngMbo/hAN11VLi35xXgIOLRJ77o4vEodOGduuo3ZMP3DUVbutm+E99dhu++IM9WDQ0lBBhTxuL16sSWXffuvO38IUju/FIo7Xunlm3EW8UbMOYQ1iijX65vd5U9bnF2/Dc4m146x7R//bXnH8n28/c2AHagkCG7w/24dqc83hTE4zjwD3iYKib3jigA0gmdd4cMhtf4bgzMgBN25OfsWTr5ur7cXvsfTxw9kNLQqgSJz3Hz3aQ99g5/AwbOnh31e8DAE6vEyk6ljVaRC9WmIOHLuzBwxffiev77KKteG7xtikfiluCeQs34A+HXiYX4nvl5fij27/hOTY343BjNI/m7fdcfh2/euyIlAC8umwNXlojT/nRu40IEx6jbbnpOfShuK1LmgAAgYEJXPePZ/C5t45OEc8f3n4Ner+Vh0COfOZ6JixGORgR3yHblhfJPMUNAAcNhnv99RYDv9iZO3WGxK1NPfNftmDvp/9t6m/rn38AAHDsQcsnc+uhLyNtaAxfeO4AHvn5R6J/RhjPbFuPZ79wM6KZIYSMcxZu/07jN+owmSOfu+pMZ7ZV9/t9f79wopp5wH5jv2Uu+rfNN+LHlTswFhTvfKRMhOusuWQVAeveTK+rYEcaQtEoHt6zG19+3yLs37txI57dth0FH9Hr7jM9P8VnT8jXnf1kuRu5rbSG7vZxuKG7aqO/Y2QE2rrLukbt5+6CGG6/dJI0ZfkK3w+H8XDp1wAA++8WB4I3v26l69nZ8neOcNxQTe3W8q9NPxy3/V++jWhz05yE484PxhEo3Xp9WJy+9gxvy6Kzdz478BSpenaFw3hw9bek7YHugamff36/2Hyf+oHYjHqGOFT3UId3QrVnF21FV6Ee52MBDMkwtwBvVq6aKnP54bWC2G74WDAws/iMDD0bCsn2yXQNaRNRqUahDdG5ngaXyLfIcAPtXF1WTjP00+3ODKYTE2KtBYNi32TnOO9/+oRgVMtWibMJZdlW/KJXWdrfqdxDPv9Pj98z9fO4QUdTbZ9j/D1C2wEQLbHEwHN94qDXkjyxFq+//jTZ92DLErI9dNx5WPNEqlgHU8W0Gpi6McfoIuUhmxn/dK4gWsv6rawGBSctxuJVQW94EZ0YsadekzKenfm3YixFrjmn99GMYzLsXJOti8RYqy6Kycs8TNdOiRXSiQ21yDiC2phUUH2q+Z9IesIxjklNw90lX5O277z07XjG8Vu/S47ZC+3f+ac5YxzzwlQF0CfDX4CVur6/SAw5t8tafDOVghmwGIYbz5XcgudKbom7bko/eSHgclZrXEK2ptQavFZspSQZyqIZtVcEiYlj5YJ41bc7GalpBvuvO4Qt2i55pielA88O2lrEO8lSiEdG5NsxPc3fe5w4QheC4mAyjGTRfk7Y6cuXeL97cBZDKyNB+tt5+dUSQTQUwnfvvBvfvdNpfi35cGbX3LlFYs5MxjET+LXL7zkEwZKREXz9xF4UFOn+AwumWytkDisAutNEaZr2uK7rj2ma9ijlv543jMM8UOMVNWWXPw58Xkjg9rTo3BH/xmzLfjwgXBzImeGS498vM7LPuLKXjRU4Jcv1jfTi94ogMTGSTm9IO8PwwuEaYbtf2yypq5AhJNO6Ee9qehzxO/CheP6mDd7MN4WJYqxbKd+MWzf7m7Bg6txq1mNMWddlE/TctZ4W469aNn0GU9dHMz3KrwbIBRYTbVGxbipD3uuGW3duLdyNJReTJ1MHtoj/N7lyzJmCoJegdkfkLZRmD8c5yCPBIFqy83A5JzPh+kCJYi6c40Y07OMA7FnHHjUS1/4W1XfeMA4TXlFTduR0xkdYvNm3jky7Ppa+eOpv37vF25E2U4gaaz8kIRTZw3SiBFkECQCsPu+v3OxwmGY80RS6nSN+g0P0cjKd7jJk+BDq/+An2wAAf/eZPZ7tq9fSaddbgoL4zrA8ugAAIABJREFUVU8kV4J21wEhMezY5F0pbnGlv602ShCi1pgYe1XAe+wZk/4qRnICy5hOt3PrjtPCOU3Dy49gYjDPZLTez/AS1H56dineClbizvI2PNhxbIqevFiyGm9cqMTthRcSrg/080fEuvjUM97rYgpzwDiM5LTu5H6f13V9t2cHG+Yd4/CKmrJj8yvxSQyBxI/4Z/fRM0QtQoA2JQHAu7f7Y0yyCBIAyGSc1RzWNNF1E2qHacrNEb+N62l/zWyiuSffV/+oRhO39gGjemKON/PuHaSZYnrInylqSZ28bXSWzZGcwFKRTq8bbt1xWjiHgSL5vtj4Lr0PZYLaxEQKXj+3FN/QRSLEB6OfBowigm9cXIHdwVrcseQMHrh4xBJUF12DN7vrEOvLnbrPaCa/ZzUdYJaftN8soEHTNIApkzHvGIeXpjHRQR/YCWRkxPlIzCP+z6auxfOFlg1281v086lFCMRLKIM1ThND6viwcd1b+msrFmJ7ZadQTSbDeeTzclpoB3XvSmtTnygWzG5Vp/XsS78SryYM2SsBdtnNA4K4DRuEKPs954Z3G1M+aVgEL1w0hKuMMzRB6K9zfuu+PrGJ8/IEEe7vst7j4kXxHosWWeP900HL+W2auewO8egYTcyLtl6a+nntkHj3rCwrnLciyzKPtbwl3tV0in+015nLqCJ90nH9+i2fkM8+ctzZ363xpBxyfm03aU6xaSC1AdFqdwON2/zDnQHxx8UxyyTVere19b2YYvYZi9EFDTIxalsOlbu9HCPiG6T1ONese9301jud1+aOGjYuByNOijgUEvfNihrrJeJUvza9bKzniPg/YEt8lOeVBMnW36s9kGFjlMPx1wJZmYgB2Dm8Hl+PfADAEFSHAYSA2JCVMfeWb4s0/rFhptBY8kxgXSI5ARN1oJvMQtO0HZqmbZdpH/OOcSQL00filT13OrLeptfoaBbKlAR4n7Ww43yJ6G8yDjfcESR2NFeLxVfT4s3cIj4dyLMNzhzU3y/ez2QcdoyNMRX2GDMX5yPIYswl2zfSwkt6gBY4jp4R73xNnbcmy2k8FMIx+tldKWJN2RmHHWOT83vdHDEKgXnVugGcgTLzDSGOYRiYDwcADd8GdF1/BSLvYI3sb+cd4xgoFQTTK5khALz7B0JSu+XvvCU6ztTFgTvhSZmSEkFFBy0FUxEkLTXi2TUt3n1XdtDUs7dbLOL8Qm9C0xsRUn1+2JupXQoLSa0s4k2A3BqDGxxxzM2VE8CyMn/mDMpHkAgKcvwRp+PNQu6WMY6l4/TceWkNiaJokv52i7NpUxI37y2VYt6r27zn/YxRLqfulPf92/PFuivv9V53mZFkE6h7n9tytG8Q333Zh/F1NRYIGmcwHLfZ+AcAtQC+I/vDecc43nukAoC84tZYNk14ZU51E3vuEf23vTbztYYTgUzTMEFFkFQ307pTxjgtefb1CMItYxx9UfFNZIyjI0MwcxkBoTQGgCeOsn4AkJ5OE49LF0R72WJv5kT5CADg3Dnx7CVLklsXPeOCuBaken+b1TW0JhtG8lrD5TTx7NIx72dzzCY9SD+bm/fWJWLeZYyjeYVYt3WnvIWuC/ni/jLGIYvGMkER/6Z1YswyxtF0fZZn30SPB3A4drtwjte/Nf+c44aGsUHTtPt0XX/FyD34qOEwb6JyD847xpF9iQ7d2fK3J33dfzRME19OQjlTLzZZ3TFvFZQLPeRARZDITFSJIo8I9QWAvBDNlEtGaMmV0hgAnjj6QcdFmnGEmHk/f55mHJypqXdCED8Z45D1SxSU1tBpVOCTMQ6/4Oa96hw9rzWf0FRxca8/bVJG/AGgtpFWNWsPes+Lafr2QgDAg0MHHRUEp9KmGxUETbRfKxw3JONI0jnul9kYJqlXXNcSyj047xjHpmfb464Fq6yMk1MDlhyo9kwq1m1pIdueFu2avRxl2NpwbglltNzpyGuqFwux0oiNz323ydE+XF/quK4XOqN9+vMFATNNXu70CW50NeSS7f02SdpT6rU5v6dGYlvDqf0WMymBYbIwrITRYuez8iDuGzVMzeltzg1fCuN3Q7j+0gN0JMKTr9/u+H00Jpix6S+ovdGqD36mRTy7rtoiBE3vW+dz8oPi+vBxa77s/b1wXb5VFExbKQ6J3lxqnYbfeW7V1M/Hm8XKC+aHPe+tG+OrrRZr7WS36+O5EF3tdCC7zYjXV58j+5980XKum4Q9YLPujtmWTVQX9w5pFoEPRK1591o3dud6oTHv47alaj9ZXjKRalwTvw8sda4LM/yjq0H8n3nZafar7TN3tRgndzLcvl8BG3MwrpvZHgDLPGa/ZhchlpkCoP2eNud21IiKsvsqJjo68TSW4mltKUazxNjTh4x3GgGCJZZmXb9HRBpy2S+uuuy48w1cUrGQR/1iOzgJhZOuNr1M17J6b4fYZHf8wNuHQx1QHEkVi1dmkuKk3rnGrmPCDLmj3ruK4/lRsWFrw/EpI5pbxca3Mw47uHf2Yjx23LyKTp9SsoheN7L7mhgxaJHMic+ZESnITEgmzuni4cs171Qc833duNMAuSGzDswE3v0dIZzITOd7/3AlAGD7Xxz1bK88kVhNlPngHJ8OFhzjOLpJEJBrDnhPSE5XYlEMMnCLUGbHNcFFeOT00ASIOqB4skQQAHsWVDvyg7TKPzphSPQSm7aXZDqT7buP04wjTZOb0mqqki9KBPCMh4PMBJYozpykU55wZsSeAdE/GSd9OlOdm1s3EUMTkJkah9MMTUFyzogTeAaNA4LZEif4aEbyznG3hh/XzqQ9n9IkJEgbnAFf6RymHEkW845xvPPbIjnc1n/2VtXPrxCLXMY4ZAf3TBzdakS3vMMkkZslyA4OmqAOKIbH6A3ESYwXBsU3q873DmvkJFO/7dtXezMME7KUFQBP8N1mLjc4xnOpV3y7svzkGMvAoHhuTra34BLOoHc5p2nsfk+Yvu6/Pd6MM5JiEGbJCXHZiXIT3Lo5myrmddW497werxB+veubvdOmcAJP4yq6Hse2l+iDpRTx5zR8zkIho0MmZNGdJjqqBL0qaZUnBLsaa47POEyboQz1++navxzOrxaLU8Y4OAlkMNOQjiSpQ/yG97lz6thhP9iXDNJSaMbDSaZ+22WaxkyAMnMBPON5eo+o7fHYrx7xbOdMTe99JJ5rz95rh5nxN1nkZ8sl2zN5YnBru+mssMkipNPzmjFKt3MCTxaR3BIAQhGaqlLEn9Pw/VooOBz6tPAs2nPreUKZqvzBPGkpQ+Upf6qh6aySgZNAPlonJAd7Cmo7qAiP2QYndcvSZZjgJFO/7X7ASfSUmSsRlObRGglnasrO8nfOYzQqKIcsNYksBxYAhCf8vTu3brj8XbLkhyY4gee6k/Q5Eg4U8ec0fM5C4RfFZxPwv8xRVJUfzDvG4XXSMtZJO5wDxVaIFedI83JW2SMuTAnFvJa2/7jjb7Nry43rwlmmV1U42muNKA0zmqrzJme6XHfRHa6QU3ofvaIW77Pa9zYIbcyu8vdX0VMcJATxfubsQxqTmurtTmdaDTfxf+YL/5+j/d6XPg8AePULLwMAvvzM7021nYqJhy0PWFJ1Sr01x0unfrKu2euBeBHmoTEr4qx2tWi3F4e6scwyU1zIKTeuiai/n/3MGappPqnpfVHjIxai580e1QQAFyLi/UyNqTG0OK6PHVm2eVsVjSfMo5XWujp/XvyxWbkRAHIPWX6NpsJ4jaW4kdbsJ0PydcUVH0vrpde8fT8CwLH1Yt3UfyTWjdbt1Oo2v+xcxFqE1i6ja6vI9pCt//77xHzKcuQdeFisC3s0qD3lyHUvi59ZvqA0jtlFR60giiVN3gvbjyMN4CUUr3BhO2TnO0yMpzBJFo0Kdl6prj9cLhahLK0Jp/LPNThzTndErhFxZjAOF9roWiAcfnULPe8c3IWg3PCrMVGYZKwxfjWW2cb5WifjuJIYKKbXy2AZfaA3YSjG4Q/cScvGz4n647LwOM6RloizajZxzUVabR8jAlzGmZoSnMp/PF9I46t7vZ2YXiVI7TgbEe1Lw97tp4OiXVZbgjPn/Mf9P5C2+TWDpc3Q/k4Wk0ylTSowAACaToj+taumLxyUl9Mvb68M6IW9N4t52/ILb3JBCTsAL/Bw91/9YfJmQC6tOTd2Lrx+4zMz47dT4bg+wZ20LDpD+w44R1rCzqokwYX/pcVo6W7TQfkmX/+Jv3ocEwH623Da0CQjFk0wq19W4MlEcUbyIbccYZUVnzLhVbN8OuCYamWuv7mbJGgnx/DNUr3JYoxhupSwA/ACD3f/ymba0EMxBy6tOTd2NrxeklPPhFfS1TiocFz/WP0mnYW04Uf+CD7nrOJyWXHhwlz4H4f0MfkmT/NpUljVQ0uWnDZUxUjFdeOz5xznQBHWRDA+Tn/b53eLeX9wu/e8c0w1GPA3dzUr5W0cw/eLm/cyh14JYQfgBR7u/hwo5rDtaZqecGO/UlAah09UHI6XStnj+jZwBwTX7+yLu6aNWH87Gk5zXnM5v81wYc24PlHo3BTZAxOO6xMuDfh8prhvxbAYZ0qUpniDS2iRKNxpMZrGJSJCZN05i4A704YI4mW3FC/eK2dwQ5XOZ6e4/AyTrj0XcN2/bRfthHygzlmdsrtbMOvCQkM7YBy89sd7EdbQcdp0d7HEukNlrpBqL56wvt3laNnUzyOjwgRqOsVvvbXRca8bjVoRmUbtiNZhurDU2Xed36YpJtZ9bUBoDZNnnNqDOcvmt00ZtSjNta1if6RMWtfszm/3mnPDy4Q5UB3/t3aRK9RnraKphCFGl+EyN5MUv5uiWHazUyjLcLkr3Wl+uvLFs4p6jfQfI869PWWeDsebm0LmJ7G9jj3gxWt32fUbz1LWNnokC78/vSFzKsISsDQPWW6sWa7FNeOYd4zDL7gDghw4H8nWV2lN4oaD9CftCYlNI9vEJ1aKTbLqZLwU1Z4u3ql81LvveHB+iy2cD2RoSLz7FOOwgXPwBlOZ5JXMszmNgNO2TIYhwzFD6q5f4R3UwWksFNImmVPnzJrjTJhzjaOrxfi37fPWLDjztB9wpaxl4ffLPhxOOCRfQ3KH+dQBQBu4qKk3flskEZT5KPweEOQWYchnXYfyIeZ0d7l4/iqPJMC9aWIDyRiHKXnKcDkmGE9pwLv/J7Xi/iuavP0BXH+unfOBFBTIlyPn4OXAPZtD0GfUU/tFMX4Z46jRaMbUMSa+bUna9M88cGuOM2GeLRKLfmmX931OLxG6yLJz3nPEaTzH1or/6w97NqOwJ/lvz1kguGAcrr6PVIOYDpSPwz+4qCkO3AFBjvHMNgpHadPTyhPyTbI4ShMATvLsN4wFpfDewBfLhNQsYxxcf66d84FkZ8uXI+fgvXxe7KLSCu+/457dNSyIY1Fmcon+9h4WBGjLWm8CtHo5LX1yjGlwUnzbEo9vyxF2bs2lMifDO3PGyftfKhLrZpkkOwen8bRXiTmrP+xNCdecoMdPMQfOAsEF43D1fWbqoK/ycfiELGrKbTM0GUDtwSErNfIMgJVQ1tMx5VxoIYfyC/JNXDDur25BiUYTxeVn6DFz/bn2VJ9nMSj0G/u7tMK7nXv24JiYtyKJ0sZJ/CfbBAGSMY6Kxf7Cv4uJfFIcYfeLKqb4WF0rvS45jWf1x/6oJsUcOAsEF4xzxaAYhz/IoqZIm6HNKdZRIV6p5LzT6Zwo4+EkFO4wEhdaOJAqNnnOeHIH0ShwkmeexpQQvUwvB64/185hZMRI1pcRb8uOc5y7UFLu69EozKAnjpL4AeDmNf4IEGfmyyEy2HKEnVtznCmpZJCe10Xd9LrhNJ4KOo8gC4o5cBYIr2AcO7hS1v1F4t191T1XKUeSgx6bRGxEPvlsVJUtRcCh7fGmqImOTlT/tBPVP/XuHmuon/rZlH5ihUKCDHQ7neX1+8UC0UaMBemKqnKHFhacdC64wzeJsW7bJ54T7Ha+d8ciobaXXBTti/dYC79DvBpKbLy17U6rYk+nQSCKwhYBSGFS5UyG5ZJ44RF6NV+6mb63e8u6zUHZR50EpTVbfOtrBsX4J238rzXVKKDVa62FkQrrW5tvPGorqxEoscbvZYqyp/3IM2O0bD6siRJr7goGxVaZyBbX3tnpXR2u45gxjmKnYDEyKtZNRrq4T7YrFVq/kcSwsl+8SYBRULrWW1QjAMEQumyUJBC1THZnI4Ix1NqKhlXuth5w+Cbx7LUfWd8zrYcOLomUy/fk4tfayL7m3koUe6b2jDF3rpQilR87f7en/PCCViVRSw3ordahPs9S1jZB9cDnxf6z57bjUiR5P3T6XeYS84JxzCQSSipGgJN+OAmGOocB8I6+/9wiPwfSeIO8VgcAlLljf10Y1MTYs3VvybM7VxC3wn7vZdFjnPwuGPJuHzHMPRlp3u2cOSh7Qr4cuXfjwD2bA+V/SQSdfWJdVpV634fzX/nRVDMC9Ng5E6WXwGIHt25kVoBE7+8HXLCNW1BzgytlPVPZdZWPwyd2/9k1AOQVtTh4ndOYT+AcfcUX5FJ+0SWmpoNO3/t8UEiRqyTE53itYFZbDnkvi0+qhGS36Zi3o/myIdVXp3m3c+agKiJXFfdubol+us8eNkxRmSnJmRCHdNE/S/PuH06jD7lx/quzOWLu1nbH359j2IvS6QN4nImSE1i4dXNou+DWsozTsvu3VI6jdYlFmE3NozYt02G2ptIIccE2lKAG8LnpZiy7rmIcCxuc9MNJT9Q5jETQ8Au5rbThYFK3nEJWjJ7ugn5aG8ofpN8pnEq356TPvF/HBCfRc8++NBZvzrGD8r8AwAVd9F8uYRwl+f40puwx+dxxDNsvOIGFWzfFbbSWLrt/dVuqZ8XN0GHn5qTSCHEpiihB7UpB05Ms5HS1+zgAoGlrCVq2WWmwTc2jek8HVpy4cl+Ik6446Yk6hzHb4ExRlZO05FnfRBO3lefoMwRlWf7qKlDg3o2T6Dlw5pzOTvH8qirv52T63EoDKYYpatL7/aoH5XPHMWwO7pPZbnACC7du1r9NRzb5FYgo8zSXoogS1K4olMaRHGrf6UDtOx2ebdOZWi9TVyCDSTB3smXqx6JVhspxUiy4CZfTvugTwVgmWoUDMMXl6Fv9sfg/xXCeu+sWuB19F29z1utwI/MSLRHlnbLaT9QLqXfTMYsAZDfTUSO99XJH5Wgu7a8p+wU9ttE8pyTaZ8Tz50W9ieOhcuEcN0uMZp62Zv5Eg3i3LYesMV22maUqYTgvbMFN47nW+LzMOYv3WfdfPJV8wrrWX2WZj7KzxXUzlUf/aqdtu8hwro8aSUFKDtBSeNcdTufzudPi21QvE8/h0qXkH7Xun494x429jsuBemFOsZsYe1Za73a0UHz3td32Z9KmM6qOS5JupCmkt9PZG9zOb3fNCy6YJmZzfnuBq+9jd3575a6jAn1kSCqqag4xI4xD07QGXdcPzcS95hqchLLu+7T33G9o4eGpTRxP0L0IgB2cKWmu0VwoiGVD+/RNVpw5hINfcw4l8SeCgzWC8Mjqcocz5y6BBGUGu9rB1ffhSl0nhDk8Oe6m3Zqm3QegD0CDrutPyPr5fmtN07YDeBzAer/38gPK1LVsJtIC+ITM0VcyEkNZZGYObnGmpLfuoc1wHHGjmBrAM7aciPdyu5AdxaUcyzlpah5LAkFUXRSSL2cO4eDXnDPbKFtMb8UTqeKbrBpPPJS1rXgU50ss/4I5PxUdqSiwnUfhmCI3725N0Y03vhQfsmoHty79gMsUwaU953LXcaWuE4XG1GuZDbhpt6ZpDQCg6/puTdNqKIXAN+MwHuJdV/EKgjJ1gTNVzSBkznWZo880GVwKR9GRYYutNzZrhWHzlhGASubw13xCXY/3PCweDGHxYDxjyLycuJHSq7SsHZz/xavs7nTQelmYS6pKnWaS9vxRXMi35tVkznndAeQXzgwzkxHuys506foYLvO8vCCRUM2LJMHmrvModZ0U5sDH4UG7vwBgl/FzM4DtAGaHcVDQNG1PAn/mfZJqjsAtQi5cmHOuy1AWCXlqHpmdwvjph0GcWQ40r7BMIaaEV/OJjkLIiVvJSBrKIiGSqVV2pksl27KBNE+GsBDBSd0ylPemo7w3fu66liWWLbMzEEVXivXtTc2jLHv2v61s3ktG0lARDUk1xdp16VjWOIrT69LRtM4ao6l51Hyio+6UfF3WHkthSzBzoDJF1O4ZxpnN+VOZbwFr39fs68Xy074enRTmyTmOPAB2RlIo+0Nl3JxhcKGLswmZqajuFFB3yntcvfVy4mY6xxNhal6Mze0c9wOZRtAVi6IHFvEyNY8CpCGXcfACQOuiMZxbZGk25nOWXAwiL4Ht0Tc0hv5hi2mamsd4fprnN5WhZco57tRKi2MhFMfiv314UG7iKxtIQ3Vf4ulfZExRNu8mZJpiSaMwSS1rHMWyxviIJ/PkuGxdapHEmIaM+NceHCJTFMUgkhdKExjanON+s+smBH8+jnWJCOi6rm9L4glSsDtD07RHPS4367q+m+ubyGA1TdujBVK2BsKzZ05KpBCUrPBKzT6xwNyL1NQ8ao9NOqSj9e867xs6TdtIi6NSpg4A6FxHfxd3ShMASO+zCCEVNQUAOS3y1BK9K+lnd6+hGUOpq1a0m/j3LqOXX3pf/P0n0q1ro5VjyEYA2fAmbvn7rb/1Mud0rQkgE0Gs8jpmUAT0V0zGMQaTyOZGU5CXl4q8UAryijxMTlVAR1zSFQsr/9pJ1FruF4R+uXF9oJ5mOul9E8hBACsl7z6aR3ZHcaMt8ue2+Gupp+iDb9G18iJdegZtFhwroNvT9reQ7WbUUs3uEdTsjh9nICPDX10kW0qT8yuElnTNz7yZjFd2XTaK0wUNyUVVzUI4RR+AAuPnPADS3Cks49B1/ckZGtS8hkxCMUP/ZBIKl/fmSkDmeK86l4KSidk7dHclQDl4UzD7Du+8rDTkZXlI70zhpkRxpj4FTfXWe7xpMJCKjtFZ9181V0+ipcYSdd++TTD66mYNK07N6qMXDOr30ILfTGXX9WGqapxBbeIlABuMn2sASJWDmYiqug/ABk3T7tN1/RW/91OYPmSOdwAYXnSFB+MBmTmoeDQNpWO0nZ5y8F5I4IQPZc7Jl0jrVxJ1Lo3VxED97Ndvr2lJQQ0t3F/1qDxBJ3vksusmhDkKx3XTbl3XD2matsGItuqjjljMRFTVKwAUwzBgSoyy3DdXI6oupk2F1drBmapmAjI7PABEZu+gu4LCtDAXznEv2p2ohUk5xxUUFBQkoBIoAnz23cSgA0md45i7QJz5wThSUkgHNpdfn0XYh0mCy+0/IlRZd+ihqXnUHhzzVV6ybL/3oSkTnKMx8yKdYM6dEsXRlzlHkTJK+xiieXR76Qf0ZhvLpyOD3PU83OjeTGt96W10/9zj8vGnMyfZue/efhedamZ4Ce0tXf4d2vY+Gaa949y3DRbmk+0pUfnamCikncOhw8yhOS6YZWU13X6Str+xfsluy5dJJVAEvLPvTjvliCrkdPVCFnroLjqjoKCwcMDV9+Gy7yaCeRRVlTAU41BQUFCQgKvvw+W2Sxjz4wBgwlCMQ0FBQWEukWQ9DmWqUlBQULiaMQdJDv1AMQ4FBQUFCfxm100U8yRXVcKYH4xjctJX5BQbxeDjcGci6Ur8INhNj92M2pIh1M0UvSmkU45QKUm4yKDCwxfJ9u5b6OiVyzfQ0TdFh+moKK49+xy9vLvXkM0kuCJXOS10RFrRYfr++afpqK3Om+iorPyTfsJD+cio1DZpNgo+IISJcrQXSvJC4LwkC7YBd/E1N4LMnprosBOM0rhrXEqRYEkx2R4HFVWloKCgsPDB1fchs+syAo0nlKlKQUFBYWGDq+9DZtedppVCQ3KmKhWOq6CgoHA1Y2EpHIpxKCgoKMwpVDhuctBjk9M/pj+DCBTLa2JMtLaRfVlHGOMI5JzfnKORqosAAJNh2ska6pOnB6fSSgC88ztSTCvT+afp+w8volPC99fR9w8xQRFce3q/fGdSdUwSATcvHEp20esyuqyUbE/rocc/mUF/e6rmxshaOiVzxhn6UJ3fgBTOec3taa4/R6uSomWTC0vlmBeMQ0FBQeGqhdI4FBQUFBSmB5UdV0FBQUFhmlAHABUUFBQUpgfFOKYPLZCCQFjukPLrOGedWa3y9mBVJd2XO+VKtvoHW9uAcc4PE45M8nQwgPwR+lR6DuNg5ZD5fjt9/5Zysr13Je3k5E5Xd66T9089RY8tVlFCtmcy35ZyPgPAeKU8oAPg18XwjXRNi8z36ZoWOlGvg+vrPJkdD8457Rd+nd9++7uh6YCWhKlqLrWUecE4FBQUFK5maCqqSkFBQUEhYehIzlSlNA4FBQWFqxgqV5WCgoKCwnSgoqoUFBQUFKYBdY5jTjCbURJc1JTf9AMpq5aT7XrrebIdK+nomABTr4N8NhPZw927t55OSZL/7x+T7eMN9LfhkDLqb2PltspTovTfUkv3fbeJbOeimsLtdH0aLuKNinpK5P5cVJifdcWl6eFq88xleqJZef4CrMcx29Gi/397d7DbRhWFcfw7LQKWVSvWSKZCLEPUPSzMG6TiCciaFTxDu2JL+gQVfQHUgLoHVV0gISEaCakSC2iwVESAJjksMlbc4NzjO9eTm5n5/ySrcZ07nrHHPpk5Z84FAESOPf9WyMzuNP9u544lcABARfPrOPJvxU+9bWZPJWXPezuIU1UA0Gt1qqpuu/tum4EEDgCorU2O48SGmT2KfsndP1zy35tmJkmb7n4350kvReAonY+jtEVAieyJ6c84/PGn9PKDlid6ds70lnNBy5GoPUSJKPkdzbtw/kwhJ67+9TL5+I2v021BogRyqq1INDZ63aPk9JXgfX35brrdytXH6f3qKCg8iMbnTo+6KEoKix58AAAFp0lEQVR+l87HUTt5nq059dRmXNnTngQLM/vIzKY5Rx+XInAAwKgdtz7keHLO0USSmW1Jkrs/kPRc0iRnPIEDAGpytTtVVXbAsafTpPg7kr7MGUzgAIDKWp2qKuDuj81s28z2JT1198c54wkcAFBVnSvH3X2n7VgCBwDURpPD9SutiqpZdRVNWhO1HImqrsJ1DypY7O3z24LY8z+SY6PKooPCthpHb6Z3z6iqKqpsirYv9d69Fiw7YgfpdQ/3m6AlSCSqmiqpbCptGVJaFVXagiiqlIzem2yudoGjz23VFy5Xf8fdPy9dHgCMTd8mcipqOWJmU0m7zbmySXMfALAq9/a3Skp7VU0kzYPFnjJrgQEAqtLksETRqaozWflNSfcXH1/lUnhJGyXrAAC9N8bkuJltSnqYWwu8qq5bCKSWXzp3QJSoOypMfqeS29IKCeDE84fPHc358M0Pycej9/WNoHCgWJDgLpmnJUouW/LReL/zILkefQ2VfqZS46P9psv5c9Zh7cnvyBCT4+f0at8709dkuqxJ1iqXwjdHJR9EvwcAwzTAGQCji0TMbHuhWVZWoywAGD2XdBS19DxnXCXrqKq6Y2ZPzSx9TgQAsETbiqpLfMSR0hxdBP2lAQBJlaukcvXiynEAGKwhJsfHIFWlEVVYlFaIRMLxhVVZqcejyqDjX54lH4+EVVtBRVhU0ebvpVueRJMlXXnrxvkPHvydHHsctASJWo5EVVde+NqHE4QF21c62VJK1y1JLqUxluMCAFpyb5kc72mOAwCwBhxxAABWN8DrOAAAHXK1q6oiOV5Xl8m2rpPnpUnOVPI/TI4H6x7NNRK1WykVJb8jqbYi4ftWWFRR+p0QvndBy5Q+J6C7/sx1sXz3NpOO10PgAICq2na75VQVAIxTD1uOEDgAoCqS4wCAHC75cYscB0cc41U6t0CU5Cx5/uiq+TAxH1z5XZrYj64cDx8PkpipOTHCsR0WLUjxfB0XPqfEgj4mp3N0snyu4wAArMxbJscLg42ZbUmaSdpcNp9SSumc4wCAUn6cfyvQzNo673A+m99fFYEDAGpylx8dZd8Kjzg+1snRhiTtSZrmDL4Mp6puvtBM3/uj2uvRih1crb0KnXFPlwjar+k5uxUl/K4Ef7dE49uUMC4It2//9cTYf9NjC1+bcPmJdVtlfJeiz0T0utdeftfM7At3/3R+/4Vm+u7w2+zl/Hnyvb/RTL+dtGQa72uS9hfuJ1pB/99lCBy/HelQM/3+c6Xn32j+fdJqdL8u+DyrbNv/Wd+KVBBve7rzeVrXr03JupW+75GuPxNly+922/M9OdRLzdS6mOHmOldmVdUDh7u/X/P559F6SUQePLadba+7Jhfvsm374pHHBZtJut78fE1SVnkmOQ4AGJ/7kibNzxNJuzmDCRwAMDLu/liSzGwqaTa/v6rqp6oAABfP3XfajuWIAwCQhcABAMhC4AAAZCFwAACyEDgAAFnMe9bOFwBQF0ccAIAsBA5ghHLbaAOLCBx4xVi+UMxsy8ymZvZZ7XW5aM3Vwvdqr0cNZrbd3O7UXpc+I3CcMeYdayxfKKWT2PRds9374S8OTLN/7zZXTE+a+2iBwLFg7DvWiL5QiiaxQW9NdPpe7+m0yR8y0avqVZPmtiN2rCErmsQG/XSmN9OmTjrEogUCxwJ2LGD4mlOTD3M7wuIUgWOJoe5YZra95L/3mlNUY1I0iQ16b+rud2uvRJ+NLnCs+OU5yB2rpI3ywNyXdKv5OXsSm74zsy1Jt8xsy90f1F6fi2Rm2/PPtplNR/hH01pw5fgZzY610/w8qh2r+UK5J+mToX+hNH9A7EmaEFDHoSl2+Uon+a3rkm6P6fO9TgSOBexYABAjcAAAsnAdBwAgC4EDAJCFwAEAyELgAABkIXAAALIQOAAAWQgcAIAs/wEY3VfLJtO4TAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Create a random sample\n", "s = 10000\n", "x = np.random.normal(0, 1, s)\n", "y = np.random.normal(0, 1, s)\n", "\n", "# Make the figure\n", "fig = plt.figure()\n", "ax = fig.gca()\n", "\n", "divider = make_axes_locatable(ax)\n", "cax = divider.append_axes('right', size='5%', pad=0.05)\n", "\n", "h, xe, ye, im = ax.hist2d(x, y, (40, 40), range=[(-2.5, +2.5), (-2.5, +2.5)])\n", "\n", "# Calculate the profile together with the standard deviation of the sample\n", "prof, _, std = hep_spt.profile(x, y, xe, std_type='sample')\n", "\n", "eb = ax.errorbar(hep_spt.cfe(xe), prof, xerr=(xe[1] - xe[0])/2., yerr=std, color='teal', ls='none')\n", "eb[-1][1].set_linestyle(':')\n", "\n", "# Calculate the profile together with the default standard deviation (that of the mean)\n", "prof, _, std = hep_spt.profile(x, y, xe)\n", "\n", "ax.errorbar(hep_spt.cfe(xe), prof, xerr=(xe[1] - xe[0])/2., yerr=std, color='r', ls='none')\n", "\n", "fig.colorbar(im, cax=cax, orientation='vertical');" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.8" } }, "nbformat": 4, "nbformat_minor": 2 }