{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "bb5ba338-91f5-440b-bf58-6cdf3f52eb63", "metadata": {}, "outputs": [], "source": [ "# Program to find measure of expected significance as a function\n", "# of a cut value x_cut applied to measured variable x\n", "# by Monte Carlo simulation of the likelihood ratio statistic.\n", "# G. Cowan / RHUL Physics / December 2022\n", "\n", "import numpy as np\n", "import scipy.stats as stats\n", "import matplotlib\n", "import matplotlib.pyplot as plt\n", "plt.rcParams[\"font.size\"] = 14" ] }, { "cell_type": "code", "execution_count": 2, "id": "acbec2e6-5774-470d-96b3-c29609a6c787", "metadata": {}, "outputs": [], "source": [ "# Define pdfs and likelihood-ratio statisic\n", "s_tot = 10.\n", "b_tot = 100.\n", "p_s = s_tot/(s_tot+b_tot)\n", "def f_s(x):\n", " return 3.*(1.-x)**2\n", "def f_b(x):\n", " return 3.*x**2\n", "def q(x):\n", " return np.log(1. + (s_tot/b_tot)*f_s(x)/f_b(x))" ] }, { "cell_type": "code", "execution_count": 3, "id": "946f3259-a296-46c1-82f8-767a1977fd19", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "....................................................................................................\n", "\n" ] } ], "source": [ "# Generate data under b and s+b hypotheses\n", "qb = []\n", "qsb = []\n", "numExp = 1000000\n", "np.random.seed(seed=1234567) # fix random seed\n", "for i in range(numExp):\n", " n = np.random.poisson(b_tot) # first b only\n", " r = np.random.uniform(0., 1., n)\n", " xb = r**(1./3.)\n", " qb.append(np.sum(q(xb)))\n", " n = np.random.poisson(s_tot+b_tot) # then s+b\n", " r1 = np.random.uniform(0., 1., n)\n", " r2 = np.random.uniform(0., 1., n)\n", " xsb = np.where(r2 < p_s, 1. - r1**(1./3.), r1**(1./3.))\n", " qsb.append(np.sum(q(xsb)))\n", " if len(qsb)%(numExp/100) == 0:\n", " print(\".\", end=\"\", flush=True)\n", "print(\"\\n\")\n", "qb = np.array(qb)\n", "qsb = np.array(qsb)" ] }, { "cell_type": "code", "execution_count": 4, "id": "718f6476-6d3d-4bc0-901a-08bfbf6756cd", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "median[q|s+b] = 21.332\n" ] } ], "source": [ "# Make and analyse histograms of q for b and s+b hypotheses\n", "nBins = 400\n", "qMin = 0.\n", "qMax = 100.\n", "qbHist, bin_edges = np.histogram(qb, bins=nBins, range=(qMin,qMax),\n", " density=True)\n", "qsbHist, bin_edges = np.histogram(qsb, bins=nBins, range=(qMin,qMax),\n", " density=True)\n", "med_q_sb = np.median(qsb)\n", "print(\"median[q|s+b] = {:.3f}\".format(med_q_sb))" ] }, { "cell_type": "code", "execution_count": 5, "id": "144414b8-f977-4056-9b59-f51f6d370e9d", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAGpCAYAAABClwgJAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAY49JREFUeJzt3XtcVGX+B/DPgYERhruQYgmkYiqlWF5azfuFMi3y7qoblmbJZuparaXiNfu1brUV1a68Esu0VZetNrPEFBXsghcoREtNUCsLVO7XgfP743SGGWYGZmDu83m/Xrwczjkz55kjM/Od5/k+30cQRVEEEREREZnMw94NICIiInI2DKCIiIiIzMQAioiIiMhMDKCIiIiIzMQAioiIiMhMDKCIiIiIzMQAioiIiMhMDKCIiIiIzKSwdwPIsMbGRvz888/w9/eHIAj2bg4REZFbEEUR5eXl6NKlCzw8jPczMYByMMnJyUhOTkZdXR0uXLhg7+YQERG5pcuXL+OWW24xul/gUi6OqbS0FEFBQbh8+TICAgLs3Rxqo8rKSnTp0gUA8PPPP0OlUtm5RURE1JKysjJ07doVJSUlCAwMNHoce6AclDxsFxAQwADKiXl6empuBwQEMIAiInISraXPMImciIiIyEwMoIiIiIjMxACKiIiIyEwMoIiIiIjMxCRyIivy9PTEhAkTNLeJiMg1MIByMHIdqIaGBns3hSygQ4cO2Lt3r72bQUREFsY6UA6qrKwMgYGBKC0tZRkDIiIiGzH185c5UERERERmYgBFZEWVlZVQqVRQqVSorKy0d3OIiMhCmANFZGVVVVX2bgIREVkYe6CIiIiIzMQeKGozUQSq1brbfBRAK8sHEREROT0GUNQmjSJw/04gv0h3+4BwYM80BlFEROTaOIRHZhONBE8AcPwXoKre9m0iIiKyJQZQZLaq+qbg6dYgIP8J4MSCpv3375R6qIiIiFwVAygyiygCU/c0/b53FqDyBjr6AH3CpG0XS6QgiiVaAQ8PD4wYMQIjRoyAhwdfbkREroI5UGQW7d6nPmGAr5d0WxCkYGr0u1IAlV8kHavytltTHYKPjw8yMjLs3QwiIrIwfiV2MMnJyejTpw8GDhxo76boad77tGeqbrK4x+9BlGzqHvZCEVHLrl69iocffhjh4eHw9PSEIAg4ceIE+vXrp1mIuz3WrFkDQRBQUFBgcP/zzz8PQRCQlZXV4uOcP38eCoUCb775ZrvbRK6BPVAOJjExEYmJiZq1eBxJtdpw75M2Xy9pX34Re6GIqGWiKGLixInIycnBjBkz0KNHDwiCgOPHj+Pbb7/Fli1brN6GU6dOQRAE9OvXr8XjevTogdmzZ2PNmjWYM2cO1ygl9kBR2zTvfZIJgrRP5u69UJWVlQgLC0NYWBiXciFqJiMjAydOnMCf//xnvP/++1i7di1WrVqFF154ASNGjMCgQYOs3oZTp04hOjoafn5+rR779NNPo6ioCK+99prV20WOjwEUmUw7EGqpzpPcCwVIvVDNi226m+LiYhQXF9u7GUQO5+DBgwCAhx56SLPt008/xaVLlzB37lyrn//q1au4evUq+vfvb9Lxt99+O/r164ctW7agsbHRyq0jR8cAikzSPP+pJc17ody5B4qI9G3duhWCIGDDhg0AgJEjR0IQBISEhCA1NRWCIGDKlClG719VVYV169YhOjoaSqUS3bt3x+uvv45jx45BEASsWrXKpHacOnUKANC/f3+kp6dj9OjR8Pf3R8eOHbFgwQJUVFTo3Wf69Om4dOkSvvjiizY8c3IlzIEikzTPf/Jp5S9Hu4dq6h7g01msTk7UEkNLIzkiSyzX1K1bNyQlJWHz5s1QqVR44oknAADh4eF47rnn0KtXLwQFBRm8b3l5OUaPHo3jx49jzJgxmDx5Mn744QcsXrwYcXFxAIDY2FiT2iEHUEeOHMHGjRsxefJk9OvXD5988glSUlJQWlqKXbt26dznD3/4AwCp92zcuHFtePbkKhhAkdmM5T9p81HoJpNXqw0nnRORpFoN9HaCCV5nFrX/tTxixAj0798f69atw7hx47BmzRoAQH5+Pq5fv4777rvP6H0TEhJw6tQp7Ny5EzNnztRs37Bhg6bnydwA6vvvv0deXh4iIiI0j9WvXz/s3r0bZ8+eRa9evTT3GTBgAADg2LFjJj9fck0cwiOzmfLtk8N4RNSSnJwciKKok3905coVAECnTp0M3ufgwYNIS0vD3LlzdYInQAqsACAgIADdunUzqQ1yAPXPf/5TEzwBgEqlwp/+9CcAQF5ens59/P390aFDB01byX2xB4pM0pYAiMN4RKbzUUi9O46uteF7U8nBi3Zv0bVr1wAAwcHBBu+TnJwMQRDw/PPP6+0LCQkBAPTr1w+CCW80ZWVl+PHHHzF48GCMGTNGb3/nzp0BADU1NQbPxYkhxACKWmVOArk2DuNJS7nIXf5cyoVaIgju9frIyckBoBtA+fj4AACqq6sN3ufAgQO47bbb0KNHD719P//8s97jtXZ+URTx4IMPGtx/6dIlAEBkZKTevurqavj6+pp0HnJdfEenVpmbQC7jMJ70gZCdnY3s7GzNhwMRST1QISEhOkNnYWFS/ZPr16/rHV9SUoKysjJ07drV4OOlp6cDMD//yVCABAD/+9//4Ofnh8GDB+tsb2xsRGlpqaat5L4YQJFZTEkg19Z8GM8dgygi0lVXV4f8/Hy9YCcmJgYeHh44d+6c3n28vKTuOXmYT1tNTQ1efvllADC5ppMcQBkK1j766CN8++23ePTRR+HtrbuUwrlz59DY2Ig77rjDpPOQ62IAZQVpaWkYN24cQkJCWlyDyRmZm8MkD+MBLKpJRJK8vDzU19frBTtBQUHo27cvjh8/DrHZty2VSoWIiAjk5OToJHbX1NRgzpw5OH/+PLy8vBATE2NSG+QA6v3334da3fTGVFBQgMWLFyM8PNxgrtXXX38NQJpJSO6NAZQVVFZWYtiwYdi4caO9m2J3zYfx3E1VVRWioqIQFRWFqqoqezeHyCEYSiCXxcfHo7S0FNnZ2Xr7li1bhsbGRgwfPhyLFi3CkiVL0Lt3b5SXl8Pb2xu9e/fW6zEypLa2Fvn5+Rg2bBgKCwsxePBgrFixAo8++ij69u2LkpISpKWlGRymS09Ph6enJyZOnGj+EyeXwgDKCubOnYvVq1dj1KhR9m6KQ9DutXK3ITxRFFFYWIjCwkK9b9RE7qqlAGr+/Pnw9PTE9u3b9fY9+eSTWLt2LXx9ffHOO+8gPT0diYmJeOWVV1BXV2dy/lNeXh7UajWGDBmCQ4cOoWPHjnjttdfw8ccfY9KkSThx4gTuvvtuvftVVVXhww8/xKRJk9ClSxeznjO5HpcJoLZv346FCxdiwIABUCqVEAQBqampLd4nOzsbEyZMQHBwMFQqFQYNGoQdO3bYpsFOxJKf+8yDIqI33ngDoiji9ttv19t38803Y/r06dixY4feAtweHh5YvXo1rly5gpqaGpw+fRrLly/XDOmZmv901113QRRFvPjii7jtttuwf/9+VFZWoqioCO+//77BWX4A8MEHH6CiogJLly418xmTK3KZAGrlypX417/+hcLCQoSHh7d6fEZGBu655x4cPXoUU6dOxRNPPIHi4mLMnj0bL7zwgg1a7BzaWsJAG/OgiMgcGzduREVFBZKTk006Pjc3F4DpM/DaQq1W44UXXsADDzyA4cOHW+085DxcJoBKSUlBQUEBioqK8Pjjj7d4rFqtxvz58yEIAo4cOYItW7Zg8+bNyM3NRUxMDJKSknRmgaxcuRKCILT446raWsJAG8sZEJE5br31Vmzbtg0qlcqk4w3VlLK0K1euYM6cOZrZfkQuU0hz7NixJh978OBBXLhwAfPmzdPp8vX398eqVaswc+ZMbN26VdMTtXz5csyfP9/ibXY25pYw0Maq5ERkjhkzZph8bG5uLiIjI40uQGwJUVFRmjX7iAAXCqDMkZGRAQAYP3683j552+HDhzXbgoKCrPrCBKRZIbW1tZrfy8rKrHq+tmhPwMOq5ERkLS2tSzdy5EgAsPp7OLkftwyg5OG56OhovX3BwcEIDQ01WMjNVNevX8elS5c09Z/y8/NRUlKCiIgIzXpNzW3atAlr165t8zkdnTyM1+cte7fEtgRBQJ8+fTS3ici2Ro4cqQmiiCzJZXKgzFFaWgoACAwMNLg/ICBAc0xbfPzxx+jfvz8eeughAMD999+P/v374+OPPzZ6nxUrVqC0tFTzc/ny5Taf31Fpxw9V9e6RC+Xr64vTp0/j9OnTXDuLiMiFuGUPlLUlJCQgISHBrPsolUoolUrrNMgB3bUFGBAO7JnGXCgiInI+btkDJfc8GetlKisrM9o7ZW3Jycno06cPBg4caJfzW5OPQgqaZMd/YUkDIiJyTm4ZQMm5T4bynG7cuIHi4mKD+VG2kJiYiPz8fIPLGNiDJYfZBEHqcTqxwHKP6eiqqqoQExODmJgYLuVCRORC3DKAkheB3L9/v94+eRsXirRMEc3mBEF39p2r50GJooj8/Hzk5+dzKRciIhfilgHUmDFj0K1bN+zYsUNTgA0AysvLsX79eigUCrNzmFyRJYpotoZLuxARkTNymSTylJQUZGZmAgC+++47zTa55lN8fDzi4+MBAAqFAikpKYiLi8OwYcMwa9YsBAQEIC0tDRcvXsSGDRvQs2dPezwNJCcnIzk5GQ0NDXY5vzHtKaLZHGtCERGRsxNEFxlXSEhIwLZt24zuT0pK0qsi+8033yApKQlffvkl6urqEBMTgyVLlmD27NlWbm3r5ET20tJSBAQE2KUNVfVA7zel22cWWTbIqaxrqgll6cd2JJWVlfDz8wMAVFRUmLw0BRER2Yepn78u0wOVmpqK1NRUs+4zaNAg7Nu3zzoNohaxdAERETkzt8yBIsfiGn2gRETkThhAORhXrgNljCsnkguCgMjISERGRnIpFyIiF8IAysE4Wh0oa5ETyYGmRHJX5Ovri4KCAhQUFHApFyIDrl69iocffhjh4eHw9PSEIAg4ceIE+vXrhwkTJrT78desWQNBEDRrkzqq559/HoIgICsrq8Xjzp8/D4VCgTfffNNGLSNjXCYHipyLuy4uTERNRFHExIkTkZOTgxkzZqBHjx4QBAHHjx/Ht99+iy1btti7iTZz6tQpCIKAfv36tXhcjx49MHv2bKxZswZz5syx2yQjYgBFdqQ9ouWqQ3hEZFxGRgZOnDiBp556Cq+++ioAoKGhAd26dcOIESMwaNAg+zbQhk6dOoXo6GjNrN2WPP3003j33Xfx2muvYeXKlTZoHRnCITwH4445UIDr5kFVV1dj4MCBGDhwIKqrq+3dHCKHcvDgQQDAQw89pNn26aef4tKlS5g7d669mtUmI0eORFRUVJvue/XqVVy9ehX9+/c36fjbb78d/fr1w5YtW9DY2Nimc1L7MYByMI6UA2XtgMYd8qAaGxtx/PhxHD9+nG90RL/bunUrBEHAhg0bAEjBhyAICAkJQWpqKgRBwJQpU4zev6qqCuvWrUN0dDSUSiW6d++O119/HceOHYMgCFi1apXJbdm3bx/uvfdeREREQKlUokuXLhg7dix27tzZ7udpqlOnTgEA+vfvj/T0dIwePRr+/v7o2LEjFixYgIqKCr37TJ8+HZcuXcIXX3xhs3aSLgZQZJA11sFrTs6DIiL30q1bNyQlJUGlUuGmm25CUlISkpKSsGnTJmRkZKBXr14ICgoyeN/y8nKMGDECSUlJiIyMxJIlS9C3b18sXrwY69atAwDExsaa1I5ly5ZhwoQJKCgowKRJk7B06VKMGTMGZ8+eRW5uroWebevkAOrIkSOYMmUKIiIiMH/+fISEhCAlJQWPPPKI3n3+8Ic/AGjqxSPbYw4UGWSLdfAAFtQk0hBFoKHK3q1onadvu1+4I0aMQP/+/bFu3TqMGzdOs0pEfn4+rl+/jvvuu8/ofRMSEnDq1Cns3LkTM2fO1GzfsGGDpufJlADqxx9/xKuvvoopU6Zg165d8PBo6k9obGxEaWlp255cG8gB1Pfff4+8vDxEREQAkJ5Tv379sHv3bpw9exa9evXS3GfAgAEAgGPHjtmsnaSLARS1ypLr4LXEIXOgtD/UPH2AhlbymDx/L1XQUNV0m8gUDVXArtYTiO1uegWgaP+SRDk5ORBFUSfv58qVKwCATp06GbzPwYMHkZaWhoSEBJ3gCZACq1WrViEgIADdunVr9fxnzpyBKIq4/fbbdYInAPDw8EBwcLC5T6nN5ADqn//8pyZ4AgCVSoU//elPSEpKQl5enk4A5e/vjw4dOmiuGdkeAygH44iLCduql2jqHuDTWQ7QKyUHTaIIHBgG3Mgx/b5Bv09BLskFgmOBP3zWtE9dCch5Xhb4Fk/kzOSgQbu36Nq1awBgNHhJTk6GIAh4/vnn9faFhIQAAPr162dS0dq+ffvCx8cHGzZswI8//ogZM2Zg9OjR8PHxafW+LT2+oX0XL140mmBeVlaGH3/8EYMHD8aYMWP09nfu3BkAUFNTo7cvJCQExcXFrbaXrIMBlINJTExEYmKiZjFDVycnkucXNSWS22Vh4fYETdpKtPImbuQAaZ2bfv9PJ6DD77eDY4GxRwGFr9SrxYCKPH2l3h1HZ6Ge1ZycHAC6AZQcvBibsXrgwAHcdttt6NGjh96+n3/+We/xWtK1a1ccPnwYq1evxs6dO/Hee+/Bx8cH06ZNw0svvWS0FwyQFqdvLjU1FSUlJViyZInePmP5XEBTT9yDDz5ocP+lS5cAAJGRkXr7qqurWaDXjhhAkV3ZtaCmqUFT816lsUf1g50WHiPU38Bj3sgBdmvtCBtq+HHJfQiCRYbGnMWpU6cQEhKiM2QVFiZNy71+/bre8SUlJSgrK8PgwYMNPl56ejoA0wMoABg4cCD27duHiooKpKenY/PmzXj33XdRXFyMvXv3Gr2fnLOlLSMjAwUFBQb3tUTuiTMUIAHA//73P/j5+ek9bzlPKyYmxqzzkeUwgCK7s2lBTTlo8vQBPrvLeNCkHSg1z2syFuTce1IvX0oFoOhPzc5vKNAqygJqi5o+QNkjRS6srq4O+fn5GDZsmM72mJgYeHh44Ny5c3r38fKSuqblYT5tNTU1ePnllwHA5FpK2vz8/PDQQw9h0qRJCAkJwZkzZ8x+jLaSAyhDQeNHH32Eb7/9Fk899RS8vb119p07dw6NjY244447bNJO0scyBuRQrFpQUxSB9HukRN2dnvpBTHAsMK1cGka59yTg5ScFNILQ1DvQUlAjH6NQAYJH023tHy8/6bGnlUvn05bWSWrbLj/gszuB+gopb8ohs+uJ2i4vLw/19fV6wU5QUBD69u2L48ePQ2z2d69SqRAREYGcnBzk5eVpttfU1GDOnDk4f/48vLy8TOqROXnyJAoLC/W2Z2RkoLy8HHfffXcbn5n55ADq/fffh1rdVAyvoKAAixcvRnh4uMGcr6+//hqANKOR7IM9UGR3Vs+Dknud1JVAcbMpv/7RUkAj9zTZotdHEJoCKXUlcChOv13aQ3yafKlWAjgiJ2EogVwWHx+PNWvWIDs7W28pl2XLlmHJkiUYPnw4Zs6cCW9vb3z00Ufo2bMnvL290atXL72eGkNee+01vPfeexgyZAj69OmDoKAg5OfnY9++fbjlllvw4osvWuR5tqa2tlbTE3f+/HkMHjwY48ePx2+//Ybdu3dDEAR8/vnnmqFNbenp6fD09MTEiRNt0lbSxx4osjurFtQUG6XenF1+Ug+PtuBYYOJZ3Z4mC6uursbIkSMxcuRI/cRYOZAalyn1ek2vMNwzJQdTn90pPR8iJ9dSADV//nx4enpi+/btevuefPJJrF27Fr6+vnjnnXeQnp6OxMREvPLKK6irqzM5/yk+Ph4zZ87E1atX8f777+PVV1/FDz/8gKVLlyInJ0cnL8ua8vLyoFarMWTIEBw6dAgdO3bEa6+9ho8//hiTJk3CiRMnDPaGVVVV4cMPP8SkSZPQpUsXm7SV9Ali835SsivtMgY//PADSktL7bLadlU90PtN6faZRdafGWeV84miFHQ0H6oLHQKM+twmPTqVlZWaxUErKiqgUpmQJNxScrs9esyIbOyPf/wj9u/fj8LCQpNeM7t27cKMGTPwyiuv6M2CW7NmDdauXdtiKQFn88477+DRRx/F4cOHMXz4cHs3x+XIs+Bb+/xlD5SDcaS18JyeurIp+PCPbspvGpcp9fw4avAh51Jp50v5R0v7ys9JvVHaeVL8DkQuZuPGjaioqEBycrJJx8vLrpgzA89ZqdVqvPDCC3jggQcYPNkZAyhyOBaJB+ShO1nzpHBnIQ/zTTzb8tAeAylyIbfeeiu2bdtmWo8tDNeUclVXrlzBnDlzNLMOyX6YRE56RFEaUrOXNlck1x76+uxOqbcGkAIPZ6+vI3g0lUloPrQnB1LBscC9J6RjiZzcjBkzTD42NzcXkZGRLRasdBVRUVFm15oi62AARTpEEZiyGzjxi23P2+6ZeHKJAoOz7E44V6+TMdqFFuUZfM0DqU96SfucraeNqB1aWg9u5MiRAFquBk7UFvyqSjqq1brB04BwKbixtnbNxBNFqQhl8+BJnmXnij0y2qUQDOVIccYeEQApgFqzZg0DKLI49kCRUScWAB19bNeR0aaK5GKjfkXxyb9KPTAOMkvNqmtVaedIaV8HuTfKVQNIIiI74zsrGeXrZb/4w6SK5GKjFCRoB09hQwFlmMMMYalUKlRWVqKystLkhNg2kXOkmvdGfdKLCeZERFbAAMrBJCcno0+fPhg4cKC9m2Jzch4U0JQHZZQoSj0ucqK4XKbAnRfk1e6NMjSkx6VhiIgshgGUg3HnOlBm5UE1VOnWeJIrirtr8KRN8NAveyDP1NvlJyWeM4giImoXBlDkUEyOf7QDgHtPOmyeT01NDe6//37cf//9qKmpsd2JtYf0mtePKsqSku7ZG0VE1GZMIifn07xIpgP3OjU0NODTTz/V3LYp7Zl68mLK8nqA8r9hQ9172JOIqI0c82s7EYx0jsiJ49pFMj2tOMvNFcj1o5RhUsCkrShLCqyIiMgsDKDIYenNxDOUOO4qRTJtQRCk3qbpFVKpBxmXgiEiMhsDKHIoLc7Ea744MGscmU+7N0rOjZJn6jG5nIjIZPz0IYdidCaeocWBGTy1nSBIvXfaCeZycjmDKCKiVvETiByOXkXy5kN3rrA4sCOQZ+ppD+eldeKQHhGRCRhAkUObukeEWFOkO3THvCfLEQT95HK5ZhQDKSIioxhAkcOR86AENGJz/Z0Q/tupaaeTDd2pVCqIoghRFK27lEt7yMnlzWtGyYEUc6PIDRQUFEAQBCQkJJi03ZpSU1MhCILmZ+bMme1+zDVr1kAQBBQUFLS/gVb0/PPPQxAEZGVltXjc+fPnda5RVFSUbRqohXWgyOEIArBnioiCD+5CjFdO046woRy6sxbtmlHqSiloknv95FIHXn52bSKRu3nwwQcRGxuL22+/3d5NsZlTp05BEAT069evxeNCQkKQlJQEAHj11Vdt0DJ9DKAcTHJyMpKTk21fdNHBCI2VmuCp0S8aHveddJgFgl2adiBVW9RUcPOzOznrkdzSzTffjDNnziAwMNDm546Pj7dpz5cjOHXqFKKjo+Hn1/IXtpCQEKxZswaA1GNnD3w3dDD2XgvPIUZqRBEdDg3T/Foz7qTTrnNXU1ODadOmYdq0abZdyqW95Nwo7VIHn/RiThS5HS8vL/Tq1Qvh4eH2bopTGDlyZJuH065evYqrV6+if//+lm2UlTCAIg1RlIpX2l1DFTxKcgAAp+tjIXo477BdQ0MD9uzZgz179jhfr6Jc6sA/Wvpdrhf12Z1SWQmiNsrIyIAgCFizZg2OHTuGUaNGwd/fH2FhYVi0aBGqq6sBAJ999hmGDh0KlUqFTp064dlnnzX4Ojpy5AgmTZqE0NBQKJVKREdHY+XKlaiqqtI7tqGhAf/3f/+HHj16oEOHDujRowc2bdqExkbDf9OGcqDq6urw+uuvIy4uDl27doVSqcRNN92EyZMn49SpUy0+35MnTyIuLg7+/v4IDAzEQw891Oa8pKqqKqxbtw7R0dFQKpXo3r07Xn/9dRw7dgyCIGDVqlUmP9a+fftw7733IiIiAkqlEl26dMHYsWOxc+fONrWtLeRr179/f6Snp2P06NHw9/dHx44dsWDBAlRUVNisLaZgAEUa1WqpeCUgJXH72GuAV6uHY2rJUUz9j8BOD3sRPKShu+bJ5Z/dxZ4oarevv/4aY8aMQWBgIBYuXIiIiAi89dZbWLBgAXbv3o3Jkyeja9euWLhwIYKCgvDSSy/hxRdf1HmMt99+GyNHjsSxY8cwceJELF68GDfffDM2btyIcePGoa6uTuf4xx57DH/961/R2NiIxMRExMXF4eWXX8ZTTz1lcruvX7+OJUuWoLa2FhMmTMDSpUsxcuRIfPrppxgyZIjREYTjx49j2LBhUCgUWLhwIQYMGIAPP/wQY8eONbuHury8HCNGjEBSUhIiIyOxZMkS9O3bF4sXL8a6desAALGxsSY91rJlyzBhwgQUFBRg0qRJWLp0KcaMGYOzZ88iNzfXrHa1hxxAHTlyBFOmTEFERATmz5+PkJAQpKSk4JFHHrFZW0wikkMqLS0VAYilpaU2O2dlnShGvCr9VNTa7LS6GhtF8dNYUXwfovg+xNterRAjXpXa5owqKipEACIAsaKiwt7NabvGRlGsKxfFj6M1/zdi9a/SdrKoiooKoz/V1dUmH1tVVdXmYysrKw0eZymHDh3SvC4+/PBDzfa6ujqxb9++oiAIYmhoqPjNN99o9pWVlYk33XST2LFjR7G+vl4URVE8ffq0qFAoxP79+4vXrl3TOcemTZtEAOLmzZv1ztuvXz+d53PlyhUxNDRUBCA+/PDDOo9z8eJFve01NTXilStX9J5XXl6e6OfnJ44dO9bo8/3ggw909s2dO1cEIO7cuVOzbevWrSIAcevWrUauoChOnjxZ9PT01LmfKIri+vXrNec6f/68ZntSUpIIQLx48aLO8RcuXBAFQRCnTJkiNjQ06OxraGgQr1+/brQNhowYMUKMjIw06z6yqVOnigDE7t27i4WFhZrtFRUVYvfu3UUA4pkzZ/TuFxkZ2eZzGmLq5y97oMggu6UbaS3X0hAUi2pwoWCHoJ1cLpOLbnI4z6L8/PyM/kyZMkXn2Jtuusnosffdd5/OsVFRUUaPHT58uM6xffr0MXicpY0cORIPPvig5ncvLy9MnToVoihi0qRJGDhwoGafv78/Jk6ciGvXruHKlSsAgH/+859Qq9V47bXXEBISovPYzzzzDMLCwnSGoN59910AwOrVq3XKitx8881m9UAplUrcfPPNettjYmIwatQoHDlyBPX19Xr7hw8fjhkzZuhsk3tVzMl7PXjwINLS0jB37ly9EgfyUGNAQAC6devW6mOdOXMGoiji9ttvh4eHbkjg4eGB4OBgk9vVXnIP1D//+U9ERERotqtUKvzpT38CAOTl5dmsPa3hLDxyHKIoTZ//Xe2oo8A550scd2kKlVROouj3Gi03cqTkcs7QozYwlCwsJ2sbGn6S9/3000+IiorCV199BUDKlTpw4IDe8V5eXjh79qzmd3k4atiwYXrHGtrWkpycHLz00kvIzMzE1atX9QKm4uJivcTzO++8E83dcsstAICSkhKTz52cnAxBEPD888/r7ZMDyX79+kEw4Ztw37594ePjgw0bNuDHH3/EjBkzMHr0aPj4+LR635Ye39C+ixcvGk0wLysrw48//ojBgwdjzJgxevs7d+4MAA41GYcBFDkO7cWCg2MBz6ZviFX1Uk6WE07Ecy1y0U11pdT7VH5O+vnsrt+LnPI/qL1aSpT19PTU+f23334zemzz3oSWEpWbH5ufnw/RBjluAQEBetsUCkWr++Rg5fr16wCAjRs3mnS+0tJSeHh4IDQ0VG9fp06dDNzDsGPHjmH06NEAgPHjx2um3QuCgA8//BC5ubmora3Vu5+hUgjyczJnksmBAwdw2223oUePHnr7fv75ZwCm5z917doVhw8fxurVq7Fz506899578PHxwbRp0/DSSy+1eF3kOkzaUlNTUVJSgiVLlujtCwoKMvpYOTk5EEVRp0dS26VLlwAAkZGRLT8hG2IARY6hWe8Txh4F0PRhfNcWYEA4sGcaP6PtTh7Om3hW6n0qPycFvrVFUg+Vpy//k9rBnIr11jrW19c5hs7lIKusrAz+/v6tHh8YGIjGxkYUFxcjLCxMZ9+vv/5q5F76Nm7ciNraWmRmZmLo0KE6+7766iurJl6XlJSgrKwMgwcPNrg/PT0dgOkBFAAMHDgQ+/btQ0VFBdLT07F582a8++67KC4uxt69e43eT67DpC0jIwMFBQUG97VEHr4zFiD973//g5+fn9HnbQ/scyfH0Lz3SaGCj0IKmmTHf5FmCjoTX19fVFRUoKKiwmk+lEwmL0YsS+sE7PLj0i9kM/KHqTyU1xq5uvXRo0f19hnaZsyFCxcQEhKiFzxVVVXh5MmTRu5lGV5eXgCAa9eu6e2rqanByy+/DMDw8Ghr/Pz88NBDD+Hw4cPw9/fHmTNn2tdYM8gBlNyrqO2jjz7Ct99+i0cffRTe3t42a1NrGECR/RnqfRIEaUmXacCJBfZrWnsJggCVSgWVSmVSPoLTkXOitBVlSb1R6koGUmRVixYtgkKhwJNPPonLly/r7S8pKdGpyyQnIq9btw6VlZWa7T/99BP+8Y9/mHzeyMhI3LhxA6dPn9Zsa2howPLly1FUVNSWp2IylUqFiIgI5OTk6CRU19TUYM6cOTh//jy8vLwQExPT6mOdPHkShYWFetszMjJQXl6Ou+++26Jtb4n8//T+++9DrW76plxQUIDFixcjPDzcYM6XPXEIzwo2bdqE//znP/j+++/h6+uLESNG4KWXXrLLYodOoaFKr/dJJgiAr5ddWkWmkHOiGqqkgEle+kX+N2yoJiAmsrTbb78db775Jp544gncdtttmDBhArp3765JSD58+DASEhLw9ttvA5Bm/c2bNw9bt27FHXfcgYceegi1tbX497//jbvvvhuffPKJSed98sknsX//ftxzzz2YPn06OnTogIyMDPz0008YOXIkMjIyrPispbpNS5YswfDhwzFz5kx4e3vjo48+Qs+ePeHt7Y1evXqZ1FPz2muv4b333sOQIUPQp08fBAUFIT8/H/v27cMtt9yiV3PLWmpra5Gfn49hw4bh/PnzGDx4MMaPH4/ffvsNu3fvhiAI+Pzzz/WGXe2NPVBWcPjwYTz55JP4+uuv8dlnn6GkpAT33XefTlRNWrR7KVr5sHW2Do3a2lokJCQgISHBYFKpSxAEKehVhhnujVJXGr4fkQUsWLAAX375JR588EF8+eWXeOWVV7Bnzx4UFxdj6dKlesnMW7ZswaZNmyAIAt544w3s27cPy5YtM2tB2okTJ2LPnj3o1q0btm/fjh07dqBXr1745ptvbJLk/OSTT2Lt2rXw9fXFO++8g/T0dCQmJuKVV15BXV2dyflP8fHxmDlzJq5evYr3338fr776Kn744QcsXboUOTk5OqUErCkvLw9qtRpDhgzBoUOH0LFjR7z22mv4+OOPMWnSJJw4ccKmvWGmEkRbTLVwc5cvX0ZERARyc3PRt29fk+5TVlaGwMBAlJaWGpyNYg1V9UDvN6XbZxbZqOdHFKXZXHIP1PQKnR6o5u3qEwZ8Ost5OjQqKys19XMqKirMSuR1SqKo3xvlH80yB0RmSE1N1fSUmbOY8K5duzBjxgy88soreoHjmjVrsHbt2hZLCTgr+fm0dUmc5kz9/HWZd7Tt27drSuMrlUoIgtDqCs3Z2dmYMGECgoODoVKpMGjQIOzYscPibSstLQUAvUJvBAOlC/QTrX0UUuAESEvNOFsiuVvR7o1qvhAxC24SmWXevHkQBEGvWKYx8uw/c2bgOavz589DEAQIgmAwj8sWXCYHauXKlSgsLERoaCjCw8NbvaAZGRmIi4uDt7c3Zs6cicDAQKSlpWH27NkoKCjAc889Z5F2NTY24i9/+QsmTJigKZhGvzOSPN6cIAB7pgJ93rJh26h95IWI5TIHchDFniiiVsXGxurUWLr99ttNul9OTo7m/q4uJCRE5xq1VGPKWlwmgEpJSUF0dDQiIyPx4osvYsWKFUaPVavVmD9/PgRBwJEjRzTTPZOSkvCHP/wBSUlJmDZtGqKjpVXoV65c2WqhNkMjoaIoYuHChbh48SKysrLa8excVAvJ4805y5AdaZEXImYQRWSW2NjYNgVBubm5iIyMtEswYWshISFm15qyNJcJoMaOHWvysQcPHsSFCxcwb948nVoZ/v7+WLVqFWbOnImtW7fihRdeAAAsX74c8+fPN6s9oihi0aJFOHDgAI4cOeJwswccghnJ4+SkDAVRrFpOZBXyGoGGjBw5EoB9empclcsEUOaQp5iOHz9eb5+87fDhw5ptQUFBZv3RiaKIxMRE7N27F4cPH0bXrl3b1V6X1Hz4zowPU057cDLNgyhWLSeyuZEjR2qCKLIMt+xHP3fuHABohui0BQcHIzQ0VHNMWyxatAg7d+7Ejh074OPjg6tXr+Lq1auoq6szep/a2lqUlZXp/Li05sN3BpLHjZm6h0GU02HVciJyMW4ZQMmz4gwt7AhI6yvJx7TF22+/jZKSEgwbNgzh4eGan2PHjhm9z6ZNmxAYGKj5cfleKzOH75x1Jp6vry9+++03/Pbbb663lIu5jFUtZ50oInJCbhlAWZsoigZ/Wuo+XbFiBUpLSzU/hpYlcBltGL6TZ+I5G0EQEBYWhrCwMNdcysUcctXy6RXAZK2FWz+7kyUOiMjpuGUAJfc8Getlkoto2ZJSqURAQIDOj0sSRSn/pQ3Dd+4ef7gEY3WiPruLQ3lE5FTcMoCSc58M5TnduHEDxcXFBvOjbCE5ORl9+vTBwIEDbX5uq39+iSKQfk9ThWrA5Wff1dbWIjExEYmJia67lEtbyHWi/H9/ncmJ5QyiiMhJuGUANWLECADA/v379fbJ2+RjbC0xMRH5+fnIzs626XlFUUrOtqqGKqBYKw8sbGiLtZ9cgVqtxptvvok333yTayE2ZyixnMN5ROQk3DKAGjNmDLp164YdO3ZoKrcCQHl5OdavXw+FQmHW+kOuoFotJWcDUrK2j7ULXEz+1eV7n8gEzRPLb+RI5Q7qK9gbRUQOzWXqQKWkpCAzMxMA8N1332m2yTWf4uPjER8fDwBQKBRISUlBXFwchg0bhlmzZiEgIABpaWm4ePEiNmzYgJ49e9rjaTiEPVNtENcoVO06CT9bXYScWK6ulHqf5GKbu/2lwIpBNhE5KJcJoDIzM7Ft2zadbVlZWZolVKKiojQBFACMGjUKmZmZSEpKwq5du1BXV4eYmBisX78es2fPtmXTdSQnJyM5ORkNDQ12a4PVPq8sGPVM3QN8OoufrS5BEAAvP6nY5md3NU0wKMqS8qKUYfyPJiKHI4iGFnEju5NnApaWltpkRl5VPdD7Ten2mUWAr5eFTyCKUg+D/OE4vcLs/CdRBCbsbBpqtEo7LayyshJ+fn4AgIqKCqhUrp3z1W7yLE3tiQbsiSIiGzL189ctc6DIDtpReVzmrLWgyAxymQNtLLZJRA6IARTZhoUWDta+G/tO3QiXfCEiB8MAysHYsw6U1bRj4eCWOMOaeD4+Prh48SIuXrwIHx8fezfHOXj66i/5ciNH6sUkInIQDKAcjL3qQFmVBYbvZM62Jp6HhweioqIQFRUFDw++3EyiveTLNK1FtdWVjh8xE5Hb4Ds6WZ+Fhu8A5kG5DTkXStB6i9IutMlgiojsjAEUWZcVhu+caTJWXV0dnn76aTz99NOoq6uzd3OcT/PhvBs5wE5PYJcf86KIyK4YQDkYl8uBsuDwnTOqr6/H5s2bsXnzZtTX19u7Oc5HHs6bVt60bp6sKIt5UURkNwygHIxL5kDJWMuH2kK70GZwrL1bQ0QEgAEUWZv2EIsVgieO4LgRefHhyb82bWMuFBHZCQMosp7m+U9W4AylDMiCmhfaTOvEXCgisgsGUGQ9Vsp/crZSBmRl8pp57I0iIhtiAOVgXCqJ3ILlC7SxlIGbM1RoM60TZ+YRkU0xgHIwLpNEbqXq41Z6OHIm2jPzQofo7uPMPCKyEYW9G0AuyoblCxy5w8HHxwd5eXma22Qh8sy8cZnS35q6UuqFAhz7D4KIXAZ7oMj6rFy+wJETyT08PBATE4OYmBgu5WINclK5dmI5h/GIyAb4jk7WZ4XgiYnkpMPTt6lG1I0cqUeKSeVEZEUMoMgpOUsieV1dHdasWYM1a9ZwKRdrkvOiZJ/dKSWVy2vnERFZGAMoclrOkEheX1+PtWvXYu3atVzKxdq0/yDKz0n/3sgBPruLPVFEZHEMoByMy5Qx4AcW2Zqh8gaAFERxZh4RWRhn4TmYxMREJCYmoqysDIGBgfZuTtvYoAI5kR55GE8OlsRGYHfA77cZ0BORZbEHiizPhiUMZPx8JAC6s/IErbc3zswjIgtjAEXWZeUSBjJHLmVAdtJ8Zl5tEf9IiMhiGECR5Wl/SFkxeGIpA2pR85l5XHiYiCyIARRZlg3zn5yllAHZkUKlm1guLzzMIIqI2olJ5ATAgp8nNs5/cvRSBh06dMA333yjuU02JvdC1RY1LfWS1kkKqmw0vExErokBFEEUpRwii+MHFDw9PZ2/JIWzEwRAGSYFTUVZ0raiLKlSuZeffdtGRE6LQ3gOxh51oKrVUg4RIOUU+VgqrLZx8MRRGTJK7oma/GvTNuZDEVE7MIByMImJicjPz0d2drZdzr9nqvN2GjniTLy6ujr87W9/w9/+9jcu5WJvck8UZ+YRkQUwgCIdzhY8OfpMvPr6ejzzzDN45plnuJSLI+DMPCKyEAZQ5NQ4E4/MZmhmHpd6ISIzMYAiy7LDN3ln6zUjOzOUD6WulH7YE0VEJuIsPLIcroFHzkJe8kUmlzgIjpWCK4WKkTkRtYg9UGQ5dlgDj6jNPH11h/IA6e93tz/zooioVeyBIuuwUw0ofuaRyeShvIaqpt5T+QuAnBel3UtFRKSFPVBkOTZaA68ljljKgByYPJTn5Qfce1I3L4qIqAXsgSLLsGP+k1zKIL+oqZSBr5ddmqKnQ4cOOHTokOY2ObDmeVFERC1gAEWWYcf8J7mUQZ+3bHZKk3l6emLkyJH2bga1hbpS+jtmMjkRGcAhPLI8O+Q/aZ+OQ3hkEZoim40scUBEehhAORh7rIXXbqIofcDI7PyN3ZHyoOrr65GcnIzk5GRWIncGzWfmFWUBn/QCdvlxZh4R6RBEke8IjqisrAyBgYEoLS1FQECAVc9VVQ/0flO6fWaRmflDogik3wMUH2vaNr3C5rkkoghM2Nm0KLLZz8NKKisr4efnBwCoqKiASsUcG4cnitIaeXJtKG2Tf5XW0+OwHpHLMvXzlz1Q1D4NVbrBU9hQu9R/4pIuZDEtJZNz7Twi+h0DKLKcyb/arf4TwE4BsiDtoTz/aN19XDuPiMBZeGRJXP6CXIV2kU1PH+DAcClwIiL6HXugiIgMkYfyBA/9xYeJyO0xgCIiag2LbBJRM+0awvviiy9w8OBBHDt2DFeuXEFxcTF8fX0RFhaGO+64AyNGjMDEiRPRuXNnS7WXHA2TackdyWU7WGiTyG2ZHUBVVFTgtddew5YtW3Dp0iXIVRA6dOiAkJAQVFdXIy8vD99++y3ef/99KBQKPPDAA1i6dCmGDh3ayqOTU7Hj8i3OQqlU4pNPPtHcJhchlzgIjgXuPSEN8xGRWzHrVf/222+jR48eWLlyJYKCgrBhwwYcPHgQZWVlqKqqwpUrV3Dt2jXU19fj7Nmz2LZtG2bMmIH9+/dj+PDhmDx5Mi5evGit50K2ZsflW1pTVe8YnWMKhQL3338/7r//figUnLPh1JoX2QSkv/9PeknVyonIrZhVSNPLywuzZ8/G008/jZiYGJNPUl1djZ07d2LTpk2YO3cuVq9e3abGuhOnKKSprpQqNAPAtHJpRXs70n4eADAgHNgzjSMsZEGiKH1xEEXgszuB8nPS9uBY4N6T/GMjcgFWKaR59uxZpKammhU8AYCPjw8eeeQRnD17Fg8//LBZ93VG//jHPxATEwM/Pz8EBQVhzJgx+Prrr+3dLOtygA8OH4UUNMmO/wJUq+3XHkBayiU1NRWpqalcysUVyMnkXn7AxLNNNaJu5LA2FJGbMSuA6t69e7tO5unpicjIyHY9hjOIiIjAyy+/jNzcXBw7dgzdu3dHXFwcrl27Zu+muTRBkHqcTixo2mbvYby6ujrMmzcP8+bNQ11dnX0bQ5YleEi9TjIuOEzkVtqV+ZiUlIT//ve/zGtq5qGHHkJcXBy6d++OPn36YPPmzSgtLUVeXp69m+byBEF3CNKRFhYmF6Td88plXojcSrsCqPXr12Pq1Kno0aMHgoODMXLkSCxZsgSpqanIzc2FWm278ZPt27dj4cKFGDBgAJRKJQRBQGpqaov3yc7OxoQJExAcHAyVSoVBgwZhx44dFm1XXV0d/vWvfyE4OBh33HGHRR/b7hz0g8JHAfQJk27nF9l/GI/cCJd5IXIb7ZoWlJWVhXXr1uHzzz9HaWkpjhw5giNHjkD4/VuZl5cXevfujdjYWPTv3x+xsbGIjY21SlL0ypUrUVhYiNDQUISHh6OwsLDF4zMyMhAXFwdvb2/MnDkTgYGBSEtLw+zZs1FQUIDnnnuuXe05evQo7rvvPlRXV6Nz585IT09HSEhIux7ToThwCQN5YeE+b9m7JURE5Kra1QN14sQJ7N+/Hxs2bMDly5dRWVmJc+fO4Y033kC3bt1QV1eHM2fOYNu2bViyZAlGjRqF4ODgdudSGZKSkoKCggIUFRXh8ccfb/FYtVqN+fPnQxAEHDlyBFu2bMHmzZuRm5uLmJgYJCUl4dy5c5rjV65cCUEQWvxpbsCAAcjJycGxY8dw3333Yfr06SguLrb487YbBy5hADhETju5A0OlDcRG5kMRuYF2BVCvvvoq4uPj8dxzz+Hmm2+Gj48PunfvjieeeAJ5eXmYPHkyRo0ahdOnT+Pf//43VqxYgXvvvdcqybRjx441OUH94MGDuHDhAv74xz+if//+mu3+/v5YtWoV1Go1tm7dqtm+fPlyXLx4scWf5nx8fNCjRw8MHjwYKSkp8PDw0HlMlzL2qENHLPwcI6uRFx3WXidvd4BU3oP5UEQurV1DeFeuXMG0adMM7lMqlXj//ffRu3dvfPbZZ1i6dKnRY20tIyMDADB+/Hi9ffK2w4cPa7YFBQUhKCioXecURRG1tbXtegyH5cDBEyAlkn86y+GbSc5KEABlmNQTK/fKAk35UFxDj8gltasHKiIiAidOnDC6X6lU4sEHH8Q777zTntNYnDw8Fx0drbcvODgYoaGhOkN45nr22WeRlZWFwsJCnDp1CgsWLMCVK1cwZcoUo/epra1FWVmZzg+1naMkkiuVSuzatQu7du3iUi6uTO6Jak4UpR8O6RG5nHYFUHPmzMEXX3yB3bt3Gz2mrq4OFy5caM9pLK60tBQAEBgYaHB/QECA5pi2+PnnnzFz5kz07NkTEyZMwK+//oqjR4+id+/eRu+zadMmBAYGan66du3a5vNTUyK5vSkUCkybNg3Tpk3jUi6uTqHSz4dKv0f64ZAekctpVwD19NNP484778Ts2bMxf/58/Pjjjzr78/Pz8cEHHyA0NLRdjXQ27733Hi5fvoza2lr88ssv+PjjjzFw4MAW77NixQqUlpZqfi5fvmyj1rouDtmRTcm9UNPKgaB+0raSXKD4mHSbJQ6IXEq7vhL7+Pjg0KFDePjhh/HOO+8gNTUVt912G7p27Ypr164hJycHDQ0NWLp0qaXaaxFyz5OxXiZ5HRxbUiqVHOJxQWq1Gv/9738BSAVW2Qvl4gRBWuZlXCaw29/erSEiK2r3u7mfnx/+85//4ODBg0hOTsYXX3yBM2fOAJBypJYtW4Ynn3yy3Q21JDn36dy5c7jrrrt09t24cQPFxcUYMmSIPZqG5ORkJCcno6GhwS7nJ8uqra3F9OnTAQAVFRUMoNwFuz+JXF67hvC0jR49Gv/5z39QUlKCGzduoKKiAgUFBVi8eLHBOkn2NGLECADA/v379fbJ2+RjbC0xMRH5+fnIzs62y/ldFVNPyCHwD5HIZVgsgNIWGBgIX1/HKqyobcyYMejWrRt27NiBnJwczfby8nKsX78eCoUCCQkJdmufU3CyDwKuiUcOQTuRnLPziJyay4wnpKSkIDMzEwDw3XffabbJNZ/i4+MRHx8PQJoZlZKSgri4OAwbNgyzZs1CQEAA0tLScPHiRWzYsAE9e/a0x9NwDg68jIs2uZRBflFTKQPthYaJrEauUF6Upbv9Ro6USO7pK83OKz4mHefgxWiJSJ9ZAdTEiROxdu1avbwhU1RXVyM5ORkqlQpPPPGE2fdvTWZmJrZt26azLSsrC1lZ0htYVFSUJoACgFGjRiEzMxNJSUnYtWsX6urqEBMTg/Xr12P27NkWb5+pnCIHysGXcZFxTTyyG3lGnjzrTmyUKpQDTb1OzWfnseAmkVMRRNH0/uN+/fohLy8PI0eOxNy5czF58uRWFwY+fvw4tm/fjh07dqCiogLbtm1zmIrkjkyeCVhaWmqVxZe1VdY1BRlnFpnQS6OulOraANKUbS8/q7avParqgd5vSrdNem4WVllZCT8/6fpUVFRApeKHpFvSfs0A+lXLp1cwgCJyEKZ+/prVA7V161YcP34cL7zwAh555BHMnz8fvXr1wp133olOnTohODgY1dXVuH79Os6dO4fjx4+jtLQUHh4emD59OjZu3IioqKj2PjeyIFGU8oPajMMORObTDp6IyCmZFUANHDgQSUlJ+PHHH7F3715s27YNGRkZ2L59u96xHh4e6Nu3L+Lj4zF//nx06dLFYo12ZbYewqtWS/lBgJQv5OMyWXH67JGr6+3trVlE2tvb2/YNIMdgLCeKiJyWWR+XgiBAFEV4eHhg0qRJmDdvHjZt2oR77rkHV65cwbVr1+Dj44OwsDDExMTYvBilK0hMTERiYqJdinnumeraHUr2WFTYy8uLMzqpKSdKXQkcimvKfyIip2VWAHXTTTfh/Pnzmt+vX7+OX375Bb17925xnTdyDq4YPDWfiVdVD6jYEUT2oF2lvKFKCqbSOtm7VUTURmbVgRo3bhz+/e9/4y9/+QuOHjWw8jiRg2m+qLCt60Gp1Wrs3bsXe/fuhVqttt2JyXEJgpQwrp00zlpQRE7HrB6ol156CadPn8Yrr7yCV199FYIg4K233kJubi769++P2NhYxMbG4pZbbrFWe8kRONmbva+X/epB1dbWYuLEiQC4lAu1IP0eqWdKDqrkWlGu2C1M5CLMejfv1KkTjh8/joyMDHzxxRfYuHEjKioq8OGHH+K///2vZsmWkJAQxMbG6gRVffr0scoTcDUOXwfKSYpoamM9KHJInr5N5QxKcqXFh4P6SftKcllgk8jBmVUHqjkPDw+sWbMGS5cuRU5ODk6dOqX5Nz8/H/X19dJJBMFxAwIHZas6UGbXSdKuZxMcC9x70ine4O1VD4p1oKhF9RVS4GSMg9dZI3JFVqkD1dw777yD7t27w9/fH8OGDcOwYU09E/X19cjLy8PJkyd11psjF8Jvx0Tt09rr58Awp/mSQuRu2hVAtTQ928vLC/3790f//v3bcwpyZHxTJ7Ksyb82rZNXktu0dh6rlBM5HLNm4RERkQXJBTYB6V9lWFOpAyJyaJwS5GCcIomciCxDe9Fh7Vl32r27YqOUe8hZeUQOhQGUg7FnJfJWOeEMPHvz9vbGG2+8oblNpEeuC2XM7t+TWDkrj8ihMIAi0zVUNS2CGhwrfSN2QrbsRPPy8kJiYqLtTkiuQbvEgawoi/lQRA6EOVDUNk78TdjW1ciJzCYP7RGRw2IARW3jZMGTvCYe0FSN3BYaGhqQkZGBjIwMx81rI8fkZK8xInfDAIrcQvM18WylpqYGo0aNwqhRo1BTU2P7BpDz0p6hJ1NXsvuUyEEwgCK3wS/05FTkYbzJvzZtS+skTeRgEEVkdwygHExycjL69OmDgQMH2rspRGRvhmboycnkRGRXDKAcTGJiIvLz85GdnW3vphCRo+JQHpHdMYAiInI2HMojsjsGUEREzohDeUR2xQCK3BK/uBMRUXuwEjm5pal7gE9nWX9mnpeXF1566SXNbSKzyeUMirLs3RIi0sIAityGXEwzv6ipmKavlWMab29vPP3009Y9Cbk2uZyBuhI4FAcUH2vap66U/uVCw0Q2xyE8chv2KqZJ1G6CAHj5AeMy9etC7fID0u8B6is4O4/IhhhAORjWgbIuW39Jb2hoQHZ2NrKzs7mUC7WfobpQgNQrtdtfN5hiIEVkVQygHAzrQLmWmpoaDBo0CIMGDeJSLmQbcjDFMgdEVsUAikzHN2Mi+9NeIy841vhxLHNAZFVMIifTiKL0jZaI7EtOKm+okl6Xu/3t3SIit8QeKDJNQxVwI0e6HRwrfQsmIvuQc6E4847IbhhAkfnGHuUbNxERuTUGUGQ+Bk9EjsHUfCgisjjmQJFpmEBO5HiYD0VkNwygqHVMIG8zLy8vJCUlaW4TWZycDyWKXPKFyIYYQFHrmEDeZt7e3lizZo29m0HuQO6Nqi2SKpQTkVUxB4rMwwRyIsfVvFK5KOov72JoGxGZjQEUmYfBk1kaGxtx+vRpnD59Go2NjfZuDrmb9Huk5V3kquSiqL+NiNqEQ3gOJjk5GcnJyVw3zUVUV1fj9ttvBwBUVFRApTKwjhmRtZTkSv9qVyUvPqa7zdDaekTUKvZAORiuhUdE7eLpa7ykQfMeJ/ZAEbUZAygiIlciCMC9J/SDKLFRfzYth/GI2owBFLmtRhGoqufnB7kgwQO49yQw+dembZ/d1TSbVnYjhwsOE7URAyhyWzFvAb3fBKbuZhBFLkgQAGVYU09U+bmmfQ/9YpcmEbkSBlDkVnwUQJ8w3W3HfwGq1fZpD5FVybWh9LbzrZ+ovfgqIrciCMCeqfZuBZENsfQIkVWwjAG5HVt+nnh5eWH58uWa20RE5BoYQBFZkbe3N/72t7/ZuxlERGRhHMIjInJnnEFB1CYMoIisqLGxEQUFBSgoKOBSLmQfnr5A2NCm38OG6i4IzlpQRG3CAMrKnnjiCQiCgDfeeMPeTaHf+SiAAeG626xVD6q6uhq33norbr31VlRXV1v+BEStkWfiTa+QfsYelZZvkcsb3MgB1BVcYJjITAygrOiTTz7Bl19+iS5duti7KaRFEIA904ATC5q23bWF9aDIhQmCFDQpVNLt5uUNdgdwgWEiMzGAspJff/0VTzzxBN577z3OvnJAggB09NHtiWI9KHIr2r1QMu1Fh4moRS4TQG3fvh0LFy7EgAEDoFQqIQgCUlNTW7xPdnY2JkyYgODgYKhUKgwaNAg7duywSHvmzZuHxYsX44477rDI45HlGeqJInIbxtbMIyKTuEwZg5UrV6KwsBChoaEIDw9HYWFhi8dnZGQgLi4O3t7emDlzJgIDA5GWlobZs2ejoKAAzz33XJvb8sYbb6CiogJ/+ctf2vwYDsWFu/QFAfBlByG5K3nNvNoiIK2TvVtD5FRcpgcqJSUFBQUFKCoqwuOPP97isWq1GvPnz4cgCDhy5Ai2bNmCzZs3Izc3FzExMUhKSsK5c03rRq1cuRKCILT4Izt79izWr1+Pd999Fx4eLnB5RVF/BXcich1yfpRMFJlQTmQCF/iEl4wdOxaRkZEmHXvw4EFcuHABf/zjH9G/f3/Ndn9/f6xatQpqtRpbt27VbF++fDkuXrzY4o/sq6++QlFREXr06AGFQgGFQoHCwkI89dRTiI2NtdjztZmGqqYV3INjdac/u6BG0Xoz8oicQvo9TCgnMoHLDOGZIyMjAwAwfvx4vX3ytsOHD2u2BQUFISgoyKTHjo+Px4ABA3S2xcXFISEhAfPmzWtbgx3F2KMuv65WzFvSvwPCpfyo9j5dhUKBRYsWaW4TObySXOlfOaFcu3eKiDTc8h1dHp6Ljo7W2xccHIzQ0FCdITxzGAq2vLy8EB4ejh49ehi9X21tLWprazW/l5WVten85jLrC6aLBk8+CqBPGJBf1LRNnpHX3vwopVKJ5OTk9j0Ikb2oK6VeZxd97RO1h8sM4ZmjtLQUABAYGGhwf0BAgOYYW9m0aRMCAwM1P127drX6OUURmLrH6qdxeIIA7J0lBVFEbsnT1/BsvLROHMojMsItAyhbKygowJ///OcWj1mxYgVKS0s1P5cvX7Z6u6rVTb0ufcKknhh35SEAn86yfEkDURRRVFSEoqIiiPwQIkfVvLCmNtaGIjLILT8y5Z4nY71MZWVlRnunrEWpVEKpVNr0nNr2TGUvvTVKGlRVVeGmm24CAFRUVEClYj4JOajmbwAP/Ah83E13myhKwRSH9YjcswdKzn0ylOd048YNFBcXG8yPsoXk5GT06dMHAwcOtOl5+V5IRDqaJ4+LImfoEWlxywBqxIgRAID9+/fr7ZO3ycfYWmJiIvLz85GdnW2X8xMRGaSuBIqPSbc5rEfkngHUmDFj0K1bN+zYsQM5OTma7eXl5Vi/fj0UCgUSEhLs1j4iIpvz9AXChkq3w4bq1nxjQV0iPS6TA5WSkoLMzEwAwHfffafZJtd8io+PR3x8PACpHk9KSgri4uIwbNgwzJo1CwEBAUhLS8PFixexYcMG9OzZ0x5PA8nJyUhOTkZDQ4Ndzk9EbkpOJJdznLR7mGqLmgrqEhEAFwqgMjMzsW3bNp1tWVlZyMrKAgBERUVpAigAGDVqFDIzM5GUlIRdu3ahrq4OMTExWL9+PWbPnm3LputITExEYmKiXRLZDWKeA5H7aL6si6x5MjkRuU4AlZqaitTUVLPuM2jQIOzbt886DXIF7LYncl/ykF5Rlr1bQuSQXCaAIitws3XwrEGhUODhhx/W3CZyGvKQnrpS+iJ1I0d6H+BQHhEABlAOx2FzoNxgHTxrUCqVZveMEjkMQQC8/IB7T0pfqEQR2O1v71YROQS3nIXnyBy2jIEbBk9M/yL6nZwbpf0+wBcIuTkGUERGTN3T/s8IURRRWVmJyspKLuVCroXFNMnNMYAi0uKjaFpUOL9IWi+wPaqqquDn5wc/Pz9UVbHwIDk57UWHb+RI+VHqSgZS5JYYQBFpEQRpXUAiMqD5osNc2oXcGAMoB2OvtfCoCdM8iFqg/QIpyZX+bW1pF1FkTxW5HAZQDsZhk8jdlCXyoIjcGhchJhfFAIqoGUvnQRG5JbnXiYsQk4tiAEXUDPOgiNqpea8TkQtiIU0iA9yw7BVR+4ii9NNQpdvrxMrl5KIYQDkYh61ETm3i6emJqVOnam4Tuaz0e6Rim3LgROTiGEA5mMTERCQmJqKsrAyBgYH2bg61U4cOHbB79257N4PIcrQXGQ4dIvU2leQ2zcgjchMMoMg4zpYhoubkWlANVVIwpa7k+njklhhAkWGiyORPIjJMXhtPvk3khjgLjwxrqGpK/gyOlb5pktkqKyshCAIEQUBlZaW9m0NERBbCAIpaN/Yov2USERFpYQBFrWPwREREpIMBlIPhWnhERESOjwGUg+FaeI6HkxGJrIALDJOTYwBF1AouKExkYVxgmFwAAygiA7igMJEVNVRxgWFyeqwDRWSAvKBwn7ek39v6BdnT0xMTJkzQ3CZye6LIiSnkEtgDRWSE9nt8W4fxOnTogL1792Lv3r3o0KGD5RpH5Kw4ZEcuggEUkREcxiOyghs5HLIjl8AAisgIeRiPiIioOQZQDoZ1oBxLe1M1KisroVKpoFKpuJQLkUzdymuBJQ7ICTCAcjCsA+V6qqqqUFXFIQsijY+7Gd/HEgfkJBhAERFR23n6AmFDm34Pjm3f47HEATkJljEgIqK2EwRpwXE50BFFYLd/6/d74MeWe6KIHBx7oIiIqH0EAVCopB9TEwcVqqbbzHkiJ8QAioiI7Is5T+SEGEAREZHlaOdENc+Hkn8PGyodJyvJlf5lzhM5EeZAEVmRh4cHRowYoblN5PK0c6Ka50ONPSrt9/RloEROjwEUkRX5+PggIyPD3s0gsi05J6p5vSd5O5EL4FdiIiIiIjMxgCIyEXNbiWzA2AvN1Jl68nEtHctZf2QBDKCITDR1j/nvt5WVlQgLC0NYWBiXciEyxYFh+ttMrU6ufZyxY1npnCyEAZSD4Vp4jsVHAfQJk27nFwHVavMfo7i4GMXFxZZtGJGrKj+nv83U6uTaxxk7lpXOyUIYQDkYroXnWAQB2DPV3q0gIiJHwwCKqBWmFlYmIiL3wQCKiIiIyEwMoIiIyDmZM5uOM+/IwlhIk4iInI88m674mLQ0zNijph3bfHkZojZiAEVkRR4eHhgwYIDmNpFbkdfFK8rSX//OkOBY4EaOaY9tzmw67WNNfXyiVjCAIrIiHx8fzqgk96W9Lp6nb+szMuS18tSVQFon27SRqI0YQBERkfWYs/4d18ojJ8IxBSIiIiIzMYAisqKqqipERUUhKioKVVWseExE5CoYQFnBmjVrIAiCzo+cSEzuRRRFFBYWorCwECKnTxMRuQzmQFlJv3798Nlnn2l+9/LysmNriIiIyJIYQFmJQqFA586d7d0MsrCqeulfHwWXeCGyG3Wl6ceKIl+sZBUuM4S3fft2LFy4EAMGDIBSqYQgCEhNTW3xPtnZ2ZgwYQKCg4OhUqkwaNAg7NixwyLtOXPmDMLDw9GjRw888sgjuHr1qkUel+zrri1A7zeBqbtZ0JjIbtI6mV7m4MAwvljJKlymB2rlypUoLCxEaGgowsPDUVhY2OLxGRkZiIuLg7e3N2bOnInAwECkpaVh9uzZKCgowHPPPdfmtgwePBipqano1asXrly5gqSkJIwePRqnTp2CUqls8+OSffgogAHhwPFfmrYd/wWoVgO+HJklsozg2NYLbZrDPxooPycVzmypyCZRG7lMD1RKSgoKCgpQVFSExx9/vMVj1Wo15s+fD0EQcOTIEWzZsgWbN29Gbm4uYmJikJSUhHPnzmmOX7lypV5SePMfbffddx+mTZuGO+64A/fddx8+/fRTFBYW4pNPPrHKcyfrEgRgzzTgzCLgxAJ7t4bIRclFNC35eERW5DI9UGPHjjX52IMHD+LChQuYN28e+vfvr9nu7++PVatWYebMmdi6dSteeOEFAMDy5csxf/78NrctNDQU3bp1w8WLF9v8GGRfgtC23iZBENCnTx/NbSIywtKvD77eyMpcJoAyR0ZGBgBg/PjxevvkbYcPH9ZsCwoKQlBQUJvPV1paioKCAkRFRbX5Mcg5+fr64vTp0/ZuBhERWZhbBlDy8Fx0dLTevuDgYISGhuoM4Znr6aefxgMPPICuXbviypUrWLVqFTp16oQJEyYYvU9tbS1qa2s1v5eVlbX5/ERERGRdLpMDZY7S0lIAQGBgoMH9AQEBmmPa4vLly5gxYwZ69uyJWbNm4eabb8YXX3wBX1/jCZKbNm1CYGCg5qdr165tPj8RERFZl1v2QFnbBx98YPZ9VqxYgWXLlml+LysrYxDlAqqqqjBw4EAAUtmMloJoIiJyHm4ZQMk9T8Z6mcrKyoz2TlmLUqlkiQMXJIoi8vPzNbeJiMg1uOUQnpz7ZCjP6caNGyguLjaYH2ULycnJ6NOnj6bXgoiIiByPWwZQI0aMAADs379fb5+8TT7G1hITE5Gfn4/s7Gyrn8toh4gomrdUgpsSRWlpF/k6Nv+diCxAXWna+5HYaPy4ll6U8uPL9+cLmEzklkN4Y8aMQbdu3bBjxw4sXrwYsbGxAIDy8nKsX78eCoUCCQkJdm2jtYkiMHWPkR3p9wDFx2zeJmczdTeQXyxVKd89Tfr9xC/S73umsQwNkUWYumTL7gDj+w4MM/3xw4ZavqgnuSSXCaBSUlKQmZkJAPjuu+802+SaT/Hx8YiPjwcgLfSbkpKCuLg4DBs2DLNmzUJAQADS0tJw8eJFbNiwAT179rTH00BycjKSk5PR0NBg1fNUq4H8Iul2nzBpuRIA0pIH2sFT2FDLLq/gQvKLpX+P/wJcr5aCJ/l3LvNC1ApPX+n9pSjL/PeZ4FhpiRZjwoYCytCm48q10jVau29RlvQ+qFCZ3h5ySy4TQGVmZmLbtm0627KyspCVlQUAiIqK0gRQADBq1ChkZmYiKSkJu3btQl1dHWJiYrB+/XrMnj3blk3XkZiYiMTERJsmsu+ZauTL1uRfAWUYv4kRkeUJgtTT01AlBU8tvc9M/lUKaDx9gIZqqad8t7/xY+X3rbFH9Y+Tz2lqzxaRES4TQKWmpiI1NdWs+wwaNAj79u2zToOciNH3LYWKwVM7CYKAyMhIzW0i0iIIpvX0KFRNxylULedEab9vGXrNmXpOola4TABF5Ih8fX1RUFBg72YQEZGFueUsPEfGMgZERESOjwGUg7FlGQMiIiJqGwZQRFZUXV2NgQMHYuDAgaiurrZ3c4iIyEKYA0VkRY2NjTh+/LjmNhERuQb2QBERERGZiQGUg2ESORERkeNjAOVgmERORETk+BhAEREREZmJARQRERGRmTgLj8jKQkND7d0EIiKyMAZQDiY5ORnJycloaGiwd1PIAlQqFYqKiuzdDCIisjAO4TkYJpETERE5PgZQRERERGZiAEVkRdXV1Rg5ciRGjhzJpVyIiFwIc6CoiSgC6kp7t8IpVdXr/+6jkJZvOXz4MACgoYFLuRA5BXUl4OkLCIJlH1cUgYaqlh/blGMsdS5qF/ZAkUQUgfR7gLRO9m6JUxqWqvv7XVuAqbuBRrFp2+w06TITkYNL6wQcGGbZF6z8HrvLz/hjm3KMpc5F7cYAysHYbSmXhiqg+FjT72FDpW8upMNHAQwIN+3Y478A17VG7U5eBarV1mkXkVvx9JXeowAgONa04wDz3teKsqT3RUvRfo819timHGOpc1G7cQjPwSQmJiIxMRFlZWUIDAy0TyMm/woow9jta4AgAHumSYFQVb3U00RENiYIwNijUmAgisBu/9aPAzicRRbFAIr0KVR8k2mBIAC+XvZuBZGbEwTpvaq1vE35OCIL4xAeERERkZnYA0VkZYI3c8mIiFwNAygiK1KpVOj6EktDEBG5Gg7hEREREZmJARQRERGRmRhAEVlRTU0NfvvX/fjtX/dDrK+xd3OIiMhCmAPlYJKTk5GcnIyGhgZ7N4UsoKGhATX5nwIAxEb+nxIRuQr2QDmYxMRE5OfnIzs7295NISIiIiMYQBERERGZiQEUERERkZkYQBERERGZiQEUERERkZk4C89BiaIIACgrK7PK41fVA42/z6ovKwPUQiXw+4LlKCsDFJwx1hrta2hMeXlTFfLGmjKUlTVAzYWIiVqm1no/Alp+TzLnWGP3McaS74WmtLMtz6Wt5yKj5M9d+XPYGEFs7QiyiytXrqBr1672bgYREZFbunz5Mm655Raj+xlAOajGxkb8/PPP8Pf3hyAIFnvcsrIydO3aFZcvX0ZAQIDFHpd4ba2N19d6eG2ti9fXeqxxbUVRRHl5Obp06QIPD+OZThzCc1AeHh4tRr7tFRAQwBeylfDaWhevr/Xw2loXr6/1WPraBgYGtnoMk8iJiIiIzMQAioiIiMhMDKDcjFKpRFJSEpRKpb2b4nJ4ba2L19d6eG2ti9fXeux5bZlETkRERGQm9kARERERmYkBFBEREZGZGEARERERmYkBFBEREZGZGEC5iezsbEyYMAHBwcFQqVQYNGgQduzYYe9mOYWffvoJr776KsaPH4+IiAh4e3ujc+fOmDJlCr7++muD9ykrK8OyZcsQGRkJpVKJyMhILFu2zGprG7qSl156CYIgQBAEfPXVVwaP4fU133//+1+MGzcOHTt2hI+PD2699VbMmjULly9f1jmO19Y8oigiLS0No0aNQnh4OHx9fXHbbbdh4cKF+PHHH/WO5/XVt337dixcuBADBgyAUqmEIAhITU01enxbruGOHTswaNAgqFQqBAcHY8KECTh+/Hj7Gi6Syzt06JDo7e0t+vn5ifPnzxf/8pe/iLfeeqsIQNy4caO9m+fwnn32WRGA2L17d/GRRx4R//rXv4pTpkwRPT09RQ8PD/Hf//63zvEVFRVibGysCEAcN26c+Oyzz4r33nuvCECMjY0VKyoq7PRMHF9+fr6oVCpFlUolAhC//PJLvWN4fc3T2NgoPvbYY5q/4UWLFonPPvusOHfuXDEiIkI8evSo5lheW/MtW7ZMBCCGh4eLjz/+uPjMM8+IcXFxoiAIor+/v/jdd99pjuX1NSwyMlIEIIaGhmpub9261eCxbbmGGzduFAGIERER4rJly8THHntMDAgIEL29vcVDhw61ud0MoFxcfX292L17d1GpVIonT57UbC8rKxNjYmJEhUIh/vDDD3ZsoeP7z3/+Ix45ckRv+5EjR0QvLy8xJCRErKmp0WxfvXq1CEB85plndI6Xt69evdrqbXZGarVaHDhwoDho0CBxzpw5RgMoXl/z/OMf/xABiImJiaJardbbX19fr7nNa2ueX375RfTw8BCjoqLE0tJSnX2vvPKKCECcN2+eZhuvr2Hp6eliQUGBKIqiuGnTphYDKHOv4Q8//CAqFAqxZ8+eYklJiWZ7Xl6e6OvrK3bv3l3nNWAOBlAu7vPPP9d7Ecs++OADEYC4YsUKO7TMNYwfP14EIGZnZ4uiKH3b79Kli+jn56f3Tai6uloMDg4Wb775ZrGxsdEezXVoGzduFL29vcW8vDzx4YcfNhhA8fqap6qqSgwJCRG7devW6ocEr635vvzySxGAOHv2bL19P/zwgwhAvP/++0VR5PU1VUsBVFuu4YoVK0QA4rZt2/Qe7/HHHxcBiJ9//nmb2socKBeXkZEBABg/frzePnnb4cOHbdkkl+Ll5QUAUCikdbnPnTuHn3/+GUOHDoVKpdI5tkOHDhg+fDh++uknnD9/3uZtdWR5eXlYu3YtVq5ciZiYGKPH8fqaJz09HdevX0d8fDwaGhqQlpaGF198EW+//bbeNeK1NV90dDS8vb2RlZWF8vJynX2ffvopAGD06NEAeH0toS3XsKXPwLi4OABt/wxUtOle5DTOnTsHQHqhNxccHIzQ0FDNMWSeS5cu4cCBA+jcuTPuuOMOAC1fb+3t586dM3qMu1Gr1UhISEDv3r3x17/+tcVjeX3NIyfJKhQK9OvXD99//71mn4eHB5YuXYrNmzcD4LVti44dO2Ljxo14+umn0bt3bzzwwAPw9/fHd999hwMHDuCxxx7Dk08+CYDX1xLacg3PnTsHPz8/dO7cucXj24IBlIsrLS0FAAQGBhrcHxAQgCtXrtiySS6hvr4ec+fORW1tLV566SV4enoCMO16ax9HwAsvvIDc3Fx8/fXXmh49Y3h9zfPbb78BAP7+97/jzjvvxDfffIPevXvj1KlTeOyxx/D3v/8d3bt3xxNPPMFr20bLly9Hly5dsHDhQrz11lua7UOGDMGcOXM0f9O8vu3XlmtYWlqKm266yeTjzcEhPCIzNTY24pFHHsGRI0ewYMECzJ07195Nclq5ubnYsGEDli9fjjvvvNPezXE5jY2NAABvb298+OGHGDhwIPz8/DBs2DDs2bMHHh4e+Pvf/27nVjq3DRs2ICEhAStWrMDly5dRUVGBzMxMqNVqjBo1CmlpafZuIlkJAygXJ0fqxiLssrIyo9E86RNFEQsWLMD27dsxZ84cvP322zr7Tbne2se5u4cffhjdu3fHmjVrTDqe19c88nUYMGAAunTporMvJiYG3bp1w4ULF1BSUsJr2wYHDx7EqlWr8Oc//xnPPfccbrnlFqhUKgwdOhSffPIJfHx8sHTpUgD827WEtlzDwMBAq11zBlAurqUx3hs3bqC4uJjj7SZqbGzEo48+infeeQezZs1CamoqPDx0X0Ktjam3NobvbnJzc3H27Fl06NBBUzxTEARs27YNAPCHP/wBgiDgww8/BMDra67bbrsNABAUFGRwv7y9urqa17YN9u7dCwAYNWqU3r6wsDDccccduHTpks77LK9v27XlGkZHR6OiogJXr1416XhzMAfKxY0YMQKbNm3C/v37MXPmTJ19+/fv1xxDLWtsbMT8+fOxdetWzJgxA++9954m70lbdHQ0unTpgqysLFRWVurMFKmpqcGRI0fQpUsX9OjRw5bNd1iPPvqowe1HjhzBuXPn8MADDyAsLAxRUVEAeH3NJX+wnzlzRm9ffX09zp8/D5VKhbCwMHTu3JnX1kx1dXUAgKKiIoP75e1KpZJ/uxbQlms4YsQIfPnll9i/fz/+9Kc/6Tze559/rjmmTdpU/ICcRn19vditWzdRqVSKp06d0mzXLqT5/fff26+BTqChoUFMSEgQAYjTpk1rtZ4Oi+W1n7E6UKLI62suuVbZli1bdLavW7dOBCDOmTNHs43X1jw7d+4UAYgxMTE6RRpFURRTU1NFAOJdd92l2cbr2zpLF9L8/vvvrVZIUxBFUWxb6EXO4tChQ4iLi4NSqcSsWbMQEBCAtLQ0XLx4ERs2bMDzzz9v7yY6tDVr1mDt2rXw8/PDU089pan5pC0+Ph6xsbEAgMrKStxzzz3IycnBuHHjcNdddyE3Nxf79u1DbGwsMjMz9WqYkK6EhARs27YNX375Je6++26dfby+5rlw4QKGDBmC3377Dffffz969eqFU6dO4eDBg4iMjMRXX32lmeLNa2uehoYGjB07FhkZGQgLC8MDDzyA4OBg5ObmIj09HUqlEgcOHMA999wDgNfXmJSUFGRmZgIAvvvuO5w8eRJDhw7V9CTFx8cjPj4eQNuu4caNG7Fy5UpERERg6tSpqKysxM6dO1FdXY3PP//c4BCsSdoUdpHT+frrr8V7771XDAwMFH18fMQBAwaI27dvt3eznILcG9LST/NvSyUlJeLSpUvFrl27il5eXmLXrl3FpUuX6n1LJcNa6oESRV5fc126dElMSEgQO3furLleiYmJ4q+//qp3LK+teWpqasT/+7//E++8807R19dXVCgU4s033yz+8Y9/1FkHT8brq6+199ikpCSd49tyDbdv3y4OGDBA9PHxEQMDA8V7771X/Oabb9rVbvZAEREREZmJs/CIiIiIzMQAioiIiMhMDKCIiIiIzMQAioiIiMhMDKCIiIiIzMQAioiIiMhMDKCIiIiIzMQAioiIiMhMDKCIiIiIzMQAioiIiMhMDKCIiIiIzMQAioiIiMhMDKCIiCykqqoK69atQ3R0NJRKJbp3747XX38dx44dgyAIWLVqlb2bSEQWorB3A4iIXEF5eTlGjx6N48ePY8yYMZg8eTJ++OEHLF68GHFxcQCA2NhY+zaSiCyGARQRkQUkJCTg1KlT2LlzJ2bOnKnZvmHDBk3PEwMoItchiKIo2rsRRETO7ODBgxgzZgwSEhKwdetWnX1XrlxB165dERAQgJKSEgiCYKdWEpElMQeKiKidkpOTIQgCnn/+eb19ISEhAIB+/foxeCJyIQygiIja6cCBA7jtttvQo0cPvX0///wzAA7fEbkaBlBERO1QUlKCsrIydO3a1eD+9PR0AAygiFwNAygionbw8vICAFy7dk1vX01NDV5++WUAQP/+/W3aLiKyLgZQRETtoFKpEBERgZycHOTl5Wm219TUYM6cOTh//jy8vLwQExNjx1YSkaUxgCIiaqdly5ahsbERw4cPx6JFi7BkyRL07t0b5eXl8Pb2Ru/eveHt7W3vZhKRBbEOFBFROz355JMoLS3Fv/71L7zzzjvo3r07EhMTMWHCBMTExDD/icgFMYAiImonDw8PrF69GqtXr9bZvmvXLgDMfyJyRRzCIyKyktzcXACcgUfkihhAERFZSU5ODgAGUESuiEu5EBFZyS233AKFQoGCggJ7N4WILIwBFBEREZGZOIRHREREZCYGUERERERmYgBFREREZCYGUERERERmYgBFREREZCYGUERERERmYgBFREREZCYGUERERERmYgBFREREZCYGUERERERmYgBFREREZKb/B0jhBj8nkZEjAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Plot histograms of q\n", "fig, ax = plt.subplots(1,1)\n", "plt.yscale(\"log\")\n", "fig.subplots_adjust(bottom=0.15, left=0.15) \n", "plt.xlabel(r'$q$', labelpad=5)\n", "plt.ylabel(r'$f(q)$', labelpad=10)\n", "plt.plot(bin_edges[:-1], qbHist, drawstyle='steps-post',\n", " label=r'$f(q|b)$', color='dodgerblue')\n", "plt.plot(bin_edges[:-1], qsbHist, drawstyle='steps-post',\n", " label=r'$f(q|s+b)$', color='orange')\n", "ax.axvline(med_q_sb, color=\"black\", linestyle=\"dashed\",\n", " label = r'median$[q|s+b]$')\n", "ax.legend(loc='upper right', frameon=False)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 6, "id": "4ed95437-715a-4d7f-ae61-e4fc4995f94b", "metadata": {}, "outputs": [], "source": [ "# Add code to calculate the p-value of the b-only hypothesis\n", "# for the median q from data generated according to s+b.\n", "# Find the corresponding significance Z." ] }, { "cell_type": "code", "execution_count": null, "id": "7b8bb4bd-43ae-4343-bf7f-3364ed466d0f", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.12.10" } }, "nbformat": 4, "nbformat_minor": 5 }