{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2b23f073",
   "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",
    "# Performance improvements by L. Moureaux / UHH IExp / February 2025\n",
    "\n",
    "import itertools as it\n",
    "from multiprocessing import Pool\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": "3b565a2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#  Define pdfs and likelihood-ratio statisic\n",
    "s_tot = 10.\n",
    "b_tot = 100.\n",
    "ps = 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 -2.*np.log1p(s_tot / b_tot * f_s(x) / f_b(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a3488b82",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate(s_tot, b_tot, seed, num):\n",
    "    # Get a rng\n",
    "    rng = np.random.default_rng(seed)\n",
    "\n",
    "    # How many events do we get per toy?\n",
    "    n_toy = np.random.poisson(s_tot + b_tot, num)\n",
    "\n",
    "    # What's the total number of events in all toys?\n",
    "    n_tot = n_toy.sum()\n",
    "\n",
    "    # Randomly assign signal events\n",
    "    ps = s_tot / (s_tot + b_tot)\n",
    "    signal = np.random.choice(n_tot, int(ps * n_tot))\n",
    "\n",
    "    # Generate x from uniform r in [0, 1]\n",
    "    #  for background:  x = r^{1/3}\n",
    "    #  for signal:      x = 1 - r^{1/3}\n",
    "    # Use float32 for efficiency\n",
    "    x = rng.random(n_tot, dtype=np.float32)\n",
    "    x **= 1 / 3\n",
    "    x[signal] = 1 - x[signal]\n",
    "\n",
    "    # Ignore divide-by-zero errors caused by arithmetic precision.\n",
    "    with np.errstate(divide='ignore'):\n",
    "        # Compute the test statistic for all events\n",
    "        q_ = q(x)\n",
    "\n",
    "    # Sum for each toy\n",
    "    return np.add.reduceat(q_, np.cumsum(n_toy) - 1)\n",
    "\n",
    "def parallel_generate(s_tot, b_tot, seeds, num, batch_size=10_000):\n",
    "    batches = (num + batch_size - 1) // batch_size\n",
    "    \n",
    "    # https://numpy.org/doc/stable/reference/random/parallel.html#seedsequence-spawning\n",
    "    seeds = seeds.spawn(batches)\n",
    "\n",
    "    # Arguments for each batch\n",
    "    args = zip(it.repeat(s_tot), \n",
    "               it.repeat(b_tot),\n",
    "               seeds,\n",
    "               it.repeat(batch_size))\n",
    "\n",
    "    q = Pool().starmap(generate, args)\n",
    "\n",
    "    return np.concatenate(q)[:num]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "fed66562",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate data under b and s+b hypotheses\n",
    "numExp = 10_000_000\n",
    "\n",
    "seeds = np.random.SeedSequence(123456)\n",
    "\n",
    "qb = parallel_generate(0, b_tot, seeds, numExp)\n",
    "qsb = parallel_generate(s_tot, b_tot, seeds, 100_000)  # We never need as many toys for qsb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "bce0e84a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "median[q|s+b]   = -41.182\n"
     ]
    }
   ],
   "source": [
    "med_q_sb = np.median(qsb)\n",
    "print(\"median[q|s+b]   = {:.3f}\".format(med_q_sb))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "68563526",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$f(q)$')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmYAAAG8CAYAAAB0YWleAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABq8UlEQVR4nO3de1yUZfo/8M8gMHI+iJYHxASPmGIJbpaCh6RMC63UUhPN0uSnmWltpuLZttzWNac15SuYpqkbtVtmeQRC28QCi7A8ooLrKiogMHJ8fn9MzzDDzAAzzOGZmc/79eLl8Bxm7sc5Xdz3dV+3TBAEAURERERkcy62bgARERERqTAwIyIiIpIIBmZEREREEsHAjIiIiEgiGJgRERERSQQDMyIiIiKJYGBGREREJBEMzIiIiIgkwtXWDSD96urqcPXqVfj4+EAmk9m6OURERNQMgiDgzp076NChA1xcjO//YmAmUVevXkVwcLCtm0FEREQmuHLlCjp16mT0eQzMJMrHxweA6on19fW1cWuIiEhKysvL0aFDBwCqP+S9vLxs3CISlZaWIjg4WP09biwGZhIlDl/6+voyMCMiIi2tWrVS3/b19WVgJkGmpiEx+V9iFAoFevfujcjISFs3hYiIiKxMJgiCYOtGkK7S0lL4+fmhpKSEPWZERKSlsrISM2fOBAB89NFHkMvlNm4RiVr6/c2hTCIiIjsjl8uRkpJi62aQBXAok4iIiEgi2GNGRERkZwRBQEVFBQDA09OT9S4dCHvMiIiI7ExFRQW8vb3h7e2tDtDIMTAwIyIiIpIIBmZEREREEsHAjIiIiEgiGJhZQGpqKh599FEEBgZCJpMhPz/f1k0iIiIiO8DAzALKy8sxePBgrF692tZNISIiIjvCchkWMGXKFADAb7/9ZuOWEBERkT1xmB6zHTt2YObMmRgwYADkcjlkMlmTVZGzsrIwatQoBAQEwMvLC1FRUdi5c6d1GkxERGSiVq1a4ZlnnsEzzzyjtaA52T+H6TFbvHgxLl26hKCgILRv3x6XLl1q9Pi0tDTExsbC3d0dEydOhJ+fH1JTUzFp0iTk5+dj0aJFVmo5ERGRcVq3bo29e/fauhl2p7BYidvlVVrbArzc0dHfw0Yt0uUwgVlSUhK6deuGkJAQvPPOO3jrrbcMHltTU4MZM2ZAJpMhIyMD/fv3BwAkJibioYceQmJiIp599ll069YNgCroaypfjGvBExERSVduYQme3fQ9lNW1Wts93Frh0OvRkgnOHCYwGzFiRLOPPXLkCM6fP49p06apgzIA8PHxwZIlSzBx4kQkJydjzZo1AIAFCxZgxowZZm8zERERWVZhsRLnrpdh1vYfAQDbpkehjZc7AODc9TLM252Dc9fLGJjZUlpaGgBg5MiROvvEbenp6ept/v7+8Pf3t2ibKisrUVlZqf69tLTULPerr9tWKizRfXzt2jW8+eabOHDgAK5fv466ujqcPHkS06dPR8eOHfH111+36P6XLVuG5cuX4+LFi+jSpYvO/rfffhtr1qxBZmYmHn74YYP3c+7cOfTs2RMbNmzA7NmzW9QmInI+5eXl8Pb2BgCUlZXBy8vLxi2SpsJiJUb8NR3K6lp4uLXC3lkPoU9HP/X+AC93eLi1wqztP0qm18wpA7OzZ88CgHqoUlNAQACCgoLUx5ji1q1buHz5srp+WV5eHoqLi9G5c2cEBgbqPWft2rVYvny5yY+pj+YLUorM3X0sCAJGjx6NnJwcTJgwAWFhYZDJZDh58iR+/vlnbNmyxSyP05js7GzIZDL069ev0ePCwsIwadIkLFu2DJMnT4avr6/F20ZE5GzOXS+DsroW6ydEIPK+QJ3vm47+Htg05UFM3XpCMr1mThmYlZSUAAD8/Pz07vf19UVBQYHJ9//vf/8b06ZNU//+xBNPAACSk5MRHx+v95y33noL8+fPV/9eWlqK4OBgk9sAALfLq9QvyLB23i26L3MTu49vl1eZ7Y2QlpaGH3/8Ea+++irWr18PAKitrUXXrl0RHR2NqKgoszxOY7Kzs9GtWzf1X7KNWbhwIT7++GNs2LABixcvtnjbiIicSWGxErO2/wgPt1Z6gzJRWDtvSfWaOWVgZmnx8fEGAzBD5HI55HK5RdoT1s5bq+vWUR05cgQAMHbsWPW2r7/+GpcvX8bSpUst/vjXrl3DtWvXEB0d3azj+/Tpg379+mHLli1YtGgRXFwcpnoNEZHNib1l26ZHNRpsafaambOzwFRO+U0g9pSJPWcNlZaWGuxNszSFQoHevXsjMjLSJo9vj5KTkyGTybBq1SoAQExMDGQyGQIDA5GSkgKZTIann37a4PkVFRVYsWIFunXrBrlcjtDQUHzwwQc4fvw4ZDIZlixZ0qx2ZGdnAwD69++PgwcPYtiwYfDx8UGbNm3w0ksvoaysTOec8ePH4/Llyzh8+LAJV05ERPpo9pY1Z8RInAwgBU4ZmIm5ZfryyG7fvo2ioiK9+WfWkJCQgLy8PGRlZdnk8e1R165dkZiYCC8vL7Rr1w6JiYlITEzE2rVrkZaWhp49exqcvHHnzh1ER0cjMTERISEhmDdvHvr27Yu5c+dixYoVAICIiIhmtUMMzDIyMvD000+jc+fOmDFjBgIDA5GUlITp06frnPPQQw8BqO/tIyKilhN7yzZNedCoHrCbEpgs55SBmTjUdODAAZ194rbmDkeR7UVHR2P+/PmoqKjAoEGDsGzZMixbtgyDBw/GrVu38MADDxg8Nz4+HtnZ2di1axcOHTqEv/zlL/j888+xcuVKfPvttwCMD8x+//135ObmIiUlBX/729+Qk5OD0NBQ7N27V2eZrgEDBgAAjh8/bsKVExFRQ8b2lgHaszMLi5UWbmHjnDIwGz58OLp27YqdO3ciJydHvf3OnTtYuXIlXF1djc4RMxcOZZomJycHgiBo1aUTJ3Dcc889es85cuQIUlNTMWXKFEycOFFrn/j8+/r6omvXrs1qgxiYffTRR+jcubN6u5eXF1544QUAQG5urtY5Pj4+aN26dYsmmxCR82nVqhVGjRqFUaNGcUmmBsSJb8b0lol5ZsrqWpuXmHKY5P+kpCRkZmYCAH755Rf1NrFmWVxcHOLi4gAArq6uSEpKQmxsLAYPHoznnnsOvr6+SE1NxcWLF7Fq1Sp0797dFpeBhIQEJCQk2DTPzR6JQZFm79bNmzcBqEqg6KNQKCCTyfD222/r7BPLmvTr1w8ymazJxy8tLcWFCxcwcOBADB8+XGf/vffeCwC4e/eu3scqKipq8jGIiEStW7fGvn37bN0MyRGLyQLG541JJc/MYQKzzMxMbNu2TWvbsWPHcOzYMQBAly5d1IEZAAwdOhSZmZlITEzEnj17UFVVhfDwcKxcuRKTJk2yZtPJDMSeT83AzMND9ZeSUqm/W/rQoUPo0aMHwsLCdPZdvXpV5/6aenxBEPDUU0/p3X/58mUAQEhIiM4+pVIJT0/PZj0OERHp17CYbICJgda562U2XT/TYQKzlJQUpKSkGHVOVFQU9u/fb5kGmUihUEChUKC2VppFYaUqOzsbgYGBWkOIbdu2BaAq+NtQcXExSktLMXDgQL33d/DgQQDG55fpC7wA4Msvv4S3t7fO49XV1aGkpATh4eHNehwiItJPs3ZnY3XLDBHzzObtzrHp+plOmWMmZZyVabyqqirk5eXpBFHh4eFwcXHRO/vWzc0NQP1wp6a7d+/i/fffBwCtnLXGiIGZviDwX//6F37++We8+OKLcHfX/gvu7NmzqKurw/3339+sxyEiAlRLMnl5ecHLywvl5eW2bo4kiDMqw9p5mxRQdfT3wKHXo7F+QoRNc80YmJHdy83NRXV1tU4Q5e/vj759++LkyZMQBEFrn5eXFzp37oycnBythPy7d+9i8uTJOHfuHNzc3JrdkyUGZp988glqamrU2/Pz8zF37ly0b99eby7bDz/8AICzgInIeBUVFaioqLB1MyRBcyamqUOYgCo4s/VKOQzMyO7pS/wXxcXFoaSkRG8P5Pz581FXV4chQ4Zg9uzZmDdvHnr16oU7d+7A3d0dvXr10unh0qeyshJ5eXkYPHgwLl26hIEDB+Ktt97Ciy++iL59+6K4uBipqanqoVVNBw8eRKtWrTB69GjjL5yIiACYNhNTqhwmx8xRWCLHTJyhIiXmbFNjgdmMGTOwcuVK7NixQ2etzDlz5qCkpASbN2/G1q1bERoaioSEBIwaNQrh4eHNzi/Lzc1FTU0NBg0ahC1btmDOnDnYsGEDPD09MWbMGCxfvlzvBIOKigp88cUXGDNmDDp06GD0dRMRkYo4jCmVmZUtwcBMYsxZLkMzkVGKWtrlLNq4cSM2btyod1/Hjh0xfvx47Ny5E2vXroWXl5d6n4uLC5YuXaqzjuaePXsAND+/7MEHH9QaKtVXuFifTz/9FGVlZXjttdeadTwREeky1zBmQ7ZaBYCBmQMTExltXSzPEGtNR169ejVSU1OhUCjwxhtvNHn8qVOnADR/RqYpampqsGbNGjz55JMYMmSIxR6HiMjRicOYTS1W3lyaqwDYYmYmAzMH19Hfw+7H21vqvvvuw7Zt25pdxFVfTTRzKygowOTJkzFlyhSLPQYRkTMw9zCmuArA1K0ncLu8ioEZkSVMmDCh2ceeOnUKISEhBhc+N4cuXbpg2bJlFrt/InJsLi4u6tncLi7OO4/PUsOYtsxVY2AmMSwwa3uNrVsZExMDABYN2oiImuLh4aFectCZmXsYUwoYmEkM18qUtpiYGHVwRkREttOSdTGljIEZERER2RVzrYspRQzMiIiI7Ex5eTm6dOkCQLXCiGYpIGfQ0nUxpYyBGRERkR1q7kxzR2bqupjNZYtaZs47lUOiFAoFevfujcjISFs3hYiISJIsHTBp1jIrLFZa9LEaYmAmMQkJCcjLy9O7tiMREZGzs1SJDE1iLTNlda3Vi7RzKJOIiIjshrVKZNhqpid7zIiIiMjuOFKJDE0MzIiIiMhu2GpxcWvhUCYREZGdcXFxwYABA9S3nYU18stsjYEZERGRnfHw8HDKSWKOuARTQ84TZhMREZFDcNT8MoCBmeSwjhkREZF+jp5fBnAoU3LMvoh58RWg4mbL78cSPNsA/sFmvctr167hzTffxIEDB3D9+nXU1dXh5MmTmD59Ojp27Iivv/66Rfe/bNkyLF++HBcvXlQvhyJFb7/9NtasWYPMzEw8/PDDBo87d+4cevbsiQ0bNmD27NlWbCERtURFRQV69+4NAMjLy4Onp6eNW2R5tsovs3YwyMDMkRVfARRRQHWFrVuin5snkHDCbMGZIAgYPXo0cnJyMGHCBISFhUEmk+HkyZP4+eefsWXLFrM8jj3Izs6GTCZDv379Gj0uLCwMkyZNwrJlyzB58mT4+vpaqYVE1BKCIODSpUvq240qvqL6t+FnraHtEmXt/DLN6v+HXo+2Wk4bAzNHVnFTFZSN2wIEdbd1a7QVnQFSX1K10UwfCmlpafjxxx/x6quvYv369QCA2tpadO3aFdHR0YiKijLL49iD7OxsdOvWDd7e3k0eu3DhQnz88cfYsGEDFi9ebIXWEZFZFRcAXj3+uH2lfnvFTaCiCNg9RfW75h/CV3OA5MdVtyenAm4aQYcFRjPMyVr5ZWL1/6lbT+B2eRUDMzKjoO5Ahwhbt8Lijhw5AgAYO3asetvXX3+Ny5cvY+nSpbZqlkliYmKQn5+P/Px8o8+9du0arl27hujo6GYd36dPH/Tr1w9btmzBokWLnGrqPZGkiakomoGSuE3QCBI2DwGmfKK6vXuK4VGSy99rB2viccmPaR/n5glM2y+p743CYiXOXS+z+uPaYpIBP4HJ7iUnJ0Mmk2HVqlUAVEGNTCZDYGAgUlJSIJPJ8PTTTxs8v6KiAitWrEC3bt0gl8sRGhqKDz74AMePH4dMJsOSJUua3Zb9+/fjscceQ+fOnSGXy9GhQweMGDECu3btavF1Nld2djYAoH///jh48CCGDRsGHx8ftGnTBi+99BLKynQ/3MaPH4/Lly/j8OHDVmsnETVCTEXZHK3699L3wLlD9ds2D6k/tloJ7Hha9SMGW26ewOTPgJfTgWnfqH5PfUl17o4/Pg8nfKLarnns5M9U+5IfV7VB8+dqjupHs1fOCgqLlRjx13TM253j0PXLROwxI7vXtWtXJCYmYt26dfDy8sIrr7wCAGjfvj0WLVqEnj17wt/fX++5d+7cwbBhw3Dy5EkMHz4c48aNw5kzZzB37lzExsYCACIiIprVjvnz5+Nvf/sbevTogTFjxsDHxweFhYU4evQoTp06heeee84cl9skMTDLyMjA6tWrMW7cOPTr1w9fffUVkpKSUFJSgj179mid89BDDwFQ9To++uijVmknEWloOFGr6IwqyBqyEMh4r75Xy81TlZ7y2Zz6Y6d8AfgFqm57tqnfrjkcmXBC+/7FXriEE7rHTtiuCt5++wo4vEK3B87M+cFNEXPL1k+IQOR9gQ5bv0zEwIzsXnR0NPr3748VK1bg0UcfxbJlywCoZirdunULjz/+uMFz4+PjkZ2djV27dmHixInq7atWrVL3lDUnMLtw4QLWr1+Pp59+Gnv27NEaDqyrq0NJSYlpF2cCMTD7/fffkZubi86dOwNQXVO/fv2wd+9e/Pbbb+jZs6f6HLGC+PHjx63WTiL6g5jvpS8ACh0OfK9Q/T5hOxDUQxUQtekHrPjjPdx5IODl1fhj+AfrD6T0bQvqoXrsb/6s3ZYJ24GKW6qet6LfVefqG261kLB23g4flAEMzMhB5OTkQBAE9O/fX72toKAAAHDPPffoPefIkSNITU1FfHy8VlAGqAK2JUuWwNfXF127dm3y8U+fPg1BENCnTx+dHC0XFxcEBAQYe0kmEwOzjz76SB2UAYCXlxdeeOEFJCYmIjc3Vysw8/HxQevWrdX/Z0RkAZqzIMWAploJ7Bin2j75M8AzqP74Rnq1ZAHB6nIZMpnMvO0UH1MMuDS3F19RBWm7p6gmDewYpwooLdiL5gy1yzQxMCOHIAYjmr1bN2+quu0NBUUKhQIymQxvv/22zr7AQNWwQL9+/Zr1ode3b194eHhg1apVuHDhAiZMmIBhw4bBw6Ppv+4au399+xqroVZaWooLFy5g4MCBGD58uM7+e++9FwBw9+5dnX2BgYEoKipqsr1EZAIxZwzQDmhEkz8DwkboP1dPsOPp6Ylff/3VAg3VeExDPWziUOfl49rDrWIvmhk5w9qYDTEwkxiFQgGFQoHa2lpbN8Wu5OTkANAOzMSgSKlU6j3n0KFD6NGjB8LCwnT2Xb16Vef+GhMcHIz09HQsXboUu3btwvbt2+Hh4YFnn30W7777rsFeOwBITEzU2ZaSkoLi4mLMmzdPZ5+hfDmgvufwqaee0rv/8uXLAICQkBCdfUql0imKVBJZRcMaYUW/1wdi5w+rbj/2jiqHC1ANH9oLcahTbHv7iPpetGn7Vb1sZhredIa1MRtiYCYxZq/87ySys7MRGBioNXTXtm1bAMCtW7d0ji8uLkZpaSkGDhyo9/4OHjwIoPmBGQBERkZi//79KCsrw8GDB7Fu3Tp8/PHHKCoqwr59+wyeJ+bEaUpLS0N+fr7efY0Rew71BV4A8OWXX8Lb21vnusU8uPDwcKMej4j00KwRJg5DirXEAFXvkpsn0HO06geQdN0wHf7BqgBMvMb2/ep/36xRpseMw5uOvDZmQwzMyO5VVVUhLy8PgwcP1toeHh4OFxcXnD17VuccNzc3APXDnZru3r2L999/HwC0ctaay9vbG2PHjsWYMWMQGBiI06dPG30fphIDM33B6L/+9S/8/PPPePXVV+Hurv0hd/bsWdTV1eH++++3SjuJHI6YM9awRpi4rbpCNVwZ1MMsvUkVFRXqNZWzsrKs39vdIUI7901ziBPQHt4UmXC9zpZfBrCOGTmA3NxcVFdX6wRR/v7+6Nu3L06ePKmzZImXlxc6d+6MnJwc5ObmqrffvXsXkydPxrlz5+Dm5tasHqSffvpJvTSKprS0NNy5cwd/+tOfTLwy44mB2SeffIKamhr19vz8fMydOxft27fXm1P3ww8/AECzi9ISkQbNmmNiYDJmg+rfy9+rAjU3z/oZlR0iWtyLJAgC8vLykJeX1/SSTJbSMA9NHOIUZ5O6ear+P9b3Uf3/GFn/TEr5Zeeul6GwWH9ajLmxx8wZFJ2xdQt0mbFN+hL/RXFxcVi2bBmysrJ0lmSaP38+5s2bhyFDhmDixIlwd3fHv/71L3Tv3h3u7u7o2bOnTs+SPhs2bMD27dsxaNAg9O7dG/7+/sjLy8P+/fvRqVMnvPPOO2a5zqZUVlaqew7PnTuHgQMHYuTIkbh+/Tr27t0LmUyGb7/9Vj3Eq+ngwYNo1aoVRo8ebZW2EjmM4iuq4Etz+TtxJqNYckKspG9Pw5WmaDiDVBzerK5Q/YhL8DVzjU4p5JeJ62WKxW2tsWYmAzNH5tmmvtqzFLl5ak/FNlFjgdmMGTOwcuVK7NixQycwmzNnDkpKSrB582Zs3boVoaGhSEhIwKhRoxAeHt7s/LK4uDhUV1fjxIkTyM7ORnV1Nbp06YLXXnsNf/7zn9GmTcuvsTlyc3NRU1ODQYMGYcuWLZgzZw42bNgAT09PjBkzBsuXL9c70aGiogJffPEFxowZgw4dOlilrUR2qWFAoVl/zM0T6PyQ/qKuEl970qw0r1Mc7vzvKWD3pD+K5mqUB2lm/pkt88s6+nvg0OvRyLp4C/N251hlzUwGZo5MsxaNFJnpw2rjxo3YuHGj3n0dO3bE+PHjsXPnTqxduxZeGkUYXVxcsHTpUp11NMWq+M3NL4uLi0NcXJxpjTcgLS3N6HMefPBBrSGNAwcONOu8Tz/9FGVlZXjttdeMfkwip9Gw3EV1eX1Cv5g71vDzzFDJCWciXr++TgKxB80AqeSXdfT3wO123lZ7PAZmjo4fDFi9ejVSU1OhUCjwxhtvNHn8qVOnABg3I9Ne1dTUYM2aNXjyyScxZMiQpk8gclYVN3UX/ZbgYt+SJHYSXP6+2SM4UsovszYGZuTw7rvvPmzbtq3ZxVP11URzVAUFBZg8eTKmTJnS9MFEzkZz6LKiwefHuC26Q5dkmH+w7uhNw/9TDVLIL7MVBmYWsHbtWnz22Wf4/fff4enpiejoaLz77rsGq7WT5U2YMKHZx546dQohISGNFnJ1FF26dDG6VhqRQzKUPwYAY/4OfPlq/bH68smsTCaTqesVmn1JJksR857FnkexIG0jPY7OVL9MxMDMAtLT0zFnzhxERkaisrISb775Jh5//HH88ssvcHXlf7nUNbZeZExMDIDGq+8TkZ3RzB8TZxVqLiqe+pIqoHg5vX7Cko17yjw9PZGfn2/TNhhNc9ZmxU3V/3Hy4xZbY9NeMUqwgG+++Ubr9//7v/9D586dkZeXh759+9qoVWQOMTEx6uCMiByEZv7Y5e9V/5q5ICz9Qfw/1CxI22ASQGGxEueul9mogbbnMAVmd+zYgZkzZ2LAgAGQy+WQyWRISUlp9JysrCyMGjUKAQEB8PLyQlRUFHbu3Gn2tpWUlACoXxibiIgkRDPXKfWl+h4yMxaEJT08g1T/Fp1RDyUXFisx4q/p6rphzpb4DzhQj9nixYtx6dIlBAUFoX379norsWtKS0tDbGws3N3dMXHiRPj5+SE1NRWTJk1Cfn4+Fi1aZJZ21dXV4fXXX8eoUaPQqVMns9wnERGZydWc+sr8k1MBtz8SzSXeQ6ZUKtUzqTMyMuDhYYcJ8pq1Nv9YV/N2uS+U1bVYPyECkfcFOl3iP+BAPWZJSUnIz8/HjRs3MGvWrEaPrampwYwZMyCTyZCRkYEtW7Zg3bp1OHXqFMLDw5GYmKi1vuLixYshk8ka/dFHEATMnDkTFy9ebLL3joiIrKz4Sn2C/7T9QMhDqt4xO+ghq6urw8mTJ3Hy5EnU1dXZujmmEXPOxm1RDR1f/h5uZYUAgLB23k4ZlAEO1GM2YsSIZh975MgRnD9/HtOmTdMqIurj44MlS5Zg4sSJSE5Oxpo1awAACxYswIwZM4xqjyAImD17Ng4dOoSMjAy9y+AQEZENFf1en0vGWmS20aAAbTdXDwyQLUTrIj/AK0TyAbIlOExgZgyxqvrIkSN19onb0tPT1dv8/f2NmoUnCAISEhKwb98+pKenIzi46RdWZWUlKisr1b+XlpY2+/GIiMgIxVdUQZnm4uJkOxoFaF1SX8I/5SuAz6Ee3nS24MxhhjKNIQ5TduvWTWdfQEAAgoKCtIYyjTV79mzs2rULO3fuhIeHB65du4Zr166hqsrw8hJr166Fn5+f+qc5wRwRETWh+Ioqj+xqjuq2WBpjx9Oq/c6wuLg98A/G9dYh6l9Lo16rX/jcyThlj5k4S9LPz0/vfl9f30ZrWTVl06ZNAIDBgwdrbT969KjBUgtvvfUW5s+fr/69tLSUwRkRUUuIQZhYCsPNExi+VPU7K/dLzm3BB96CHHJXF/jeNwA4YesW2YZTBmaWprmQdHPJ5XLI5XILtIaIyEmJOWTjtgCegaqhy2/+LInK/aSr2rsjRlS+h+TnI9HD548RpqIzkpoha42F1Z0yMBN7ysSes4ZKS0sN9qZZmkKhgEKhQG1trU0en4jIbhVfqR/6qiiqzyETg7CEEw5VLDYoKMjWTTCrm+VVuIogVHt3BDxLdUpp2PI5C/Byh4dbK8za/iMOvR5t0RmjTpljJuaW6csju337NoqKivTmn1lDQkIC8vLykJWVZZPHJyKyS+Kw5eZo1Y++HDIHKhbr5eWFGzdu4MaNG/Dy8rJ1c1qssFiJWdt/rC8q27CUho1zzTr6e2DTlAehrK7FbQv3mjllYBYdHQ0AOHDggM4+cZt4DBERSZiY3P/fU/XDli+nq34STrAMhp24XV4FZXUtNk15sL43yj8YCOpu24ZpsNaC6k45lDl8+HB07doVO3fuxNy5cxEREQEAuHPnDlauXAlXV1fEx8fbpG0cyiQiaqaGyf2A6oucwZjdslbwI2UO02OWlJSE+Ph4xMfHY+/evTrbvvjiC/Wxrq6uSEpKQl1dHQYPHoyXX34ZCxYsQL9+/fDrr79i2bJl6N7dNlE6hzIdQ35+PmQymU6Ab2i7JaWkpGitUjFx4sQW3+eyZcsgk8mQn5/f8gZa0Ntvvw2ZTIZjx441ety5c+e0/o+6dOlinQZSy4iLjw9ZaOuWWJ1SqURMTAxiYmKgVCpt3RzLKymoL3ni4BymxywzMxPbtm3T2nbs2DH1B3KXLl0QFxen3jd06FBkZmYiMTERe/bsQVVVFcLDw7Fy5UpMmjTJmk0nsoqnnnoKERER6NOnj62bYjXZ2dmQyWTo169fo8cFBgYiMTERALB+/XortIxMJn4xa+aJtY9QJYgDqsR+J1BXV6cuhG63SzJpaHK24+4/vpfdPFV5gw7cK+owgVlKSorR61FGRUVh//79lmkQkR4dO3bE6dOnbTLrNy4uzmZD9LaSnZ2Nbt26wdvbu9HjAgMDsWzZMgDgurZSJg5dAsCE7UDFLdVtv06qfDLAIRL7nY1O4r8hj70DHF6hWt/UgVcEcJjAzFEwx8yxubm5oWfPnrZuht2IiYlBfn6+SUOm4oobnMjjQMShS6B+1qWbp8OUv3BWYuL/tulRumUoPNvU94b2HA0EdVM99xU3HfY5d5gcM0fBHDPjpaWlQSaTYdmyZTh+/DiGDh0KHx8ftG3bFrNnz1bnX3zzzTd4+OGH4eXlhXvuuQdvvvmm3gA4IyMDY8aMQVBQEORyObp164bFixejoqJC59ja2lr85S9/QVhYGFq3bo2wsDCsXbvW4NCCvhyzqqoqfPDBB4iNjUVwcDDkcjnatWuHcePGITs7u9Hr/emnnxAbGwsfHx/4+flh7NixJud9VVRUYMWKFejWrRvkcjlCQ0PxwQcf4Pjx45DJZFiyZEmz72v//v147LHH0LlzZ8jlcnTo0AEjRozArl27TGqbKcT/u/79++PgwYMYNmwYfHx80KZNG7z00ksoKyuzWlvITCqKVP9O+ER75qWDfkE7G72J/2LZDPF59nSs2m36sMfMwZWXlxvc16pVK7Ru3bpZx7q4uMDDw8OkYysqKvSuhmDu2js//PAD/vKXvyA2NhYzZ87E0aNH8Y9//AOlpaV46qmnMHXqVDz55JMYOHAg9u3bh3fffRe+vr54++231fexadMmzJ49GwEBARgzZgzatm2LrKwsrF69GkePHsXRo0fh7l7/4fHyyy9j69atuO+++5CQkIC7d+/i/fffx/Hjx5vd7lu3bmHevHkYPHgwRo0ahYCAAFy4cAH//ve/sX//fmRkZCAyMlLnvJMnT+K9995DTEwMZs6ciezsbHzxxRf45ZdfkJubq/XcNuXOnTsYNmwYTp48ieHDh2PcuHE4c+YM5s6di9jYWABQz15uyvz58/G3v/0NPXr0wJgxY+Dj44PCwkIcPXoUp06dwnPPPdfsdrWEGJhlZGRg9erVGDduHPr164evvvoKSUlJKCkpwZ49e6zSFjKD4iv1BWPb92Mw5kyc7bkWSFI2btwo9OrVS+jevbsAQCgpKWnR/QEw+DNq1CitYz09PQ0eGx0drXVsUFCQwWMHDBigdWxISIje48zl6NGj6vv84osv1NurqqqEvn37CjKZTAgKChJOnDih3ldaWiq0a9dOaNOmjVBdXS0IgiD8+uuvgqurq9C/f3/h5s2bWo+xdu1aAYCwbt06ncft16+fUFZWpt5eUFCg/v+ZOnWq1v1cvHhRZ/vdu3eFgoICnevKzc0VvL29hREjRhi83k8//VRr35QpUwQAwq5du9TbkpOTBQBCcnKygf9BQRg3bpzQqlUrrfMEQRBWrlypfqxz586ptycmJgoAhIsXL2odf/78eUEmkwlPP/20UFtbq7WvtrZWuHXrlsE26BMdHS2EhIQYdY7omWeeEQAIoaGhwqVLl9Tby8rKhNDQUAGAcPr0aZ3zQkJCTH5MsqDCbEFI9BWEswdt3RJJKCsrU783NT9/7E3B7Qrh858KhJA3vxJ+KShu+gTxdVCYbemm6filoLhZ7SwpKWnR9zeHMiWGQ5mmi4mJwVNPPaX+3c3NDc888wwEQcCYMWO0ep18fHwwevRo3Lx5U71g/UcffYSamhps2LABgYGBWvf9xhtvoG3btlpDcR9//DEAYOnSpVq9fx07dsSrr77a7HbL5XJ07NhRZ3t4eDiGDh2KjIwMVFdX6+wfMmQIJkyYoLVt+vTpAGDU6+fIkSNITU3FlClTdEppiEOuvr6+6Nq1a5P3dfr0aQiCgD59+sDFRfvjxcXFBQEBAc1uV0uJPWYfffQROnfurN7u5eWFF154AQCQm5trtfaQkcTCsQ1LJDjBUFZzeXp6wtPT09bNMFlhsRIj/pqOebtzmk78dyIcynRwjeXRtGrVSuv369evGzy24ZdsY3lMDY/Ny8szaWF3Y/Xv319nW/v27QHoH4YT9xUWFqJLly74z3/+A0CVi3bo0CGd493c3PDbb7+pfz916hQAYPDgwTrH6tvWmJycHLz77rvIzMzEtWvXdAKxoqIidXtFDzzwgM79dOrUCQBQXFzc7MdWKBSQyWRaQ7oiMUDt168fZDJZk/fVt29feHh4YNWqVbhw4QImTJiAYcOGaQ1tG9LY/evbd/HiRYP1xkpLS3HhwgUMHDgQw4cP19l/7733AgDu3r3bZLvISjTLYDQsHOvmCQxfaru2SZCXl1ejKSX2QEz6Xz8hApH3BRq3/qSYb+iAGJg5OGPyuCx1rLX+ovP19dXZ5urq2uQ+MQi6dUs19X716tXNerySkhK4uLjoXUj4nnvuaV6jARw/fhzDhg0DAIwcOVJd3kEmk+GLL77AqVOnUFlZqXOevpIb4jUZM6v30KFD6NGjB8LCwnT2Xb16FUDz88uCg4ORnp6OpUuXYteuXdi+fTs8PDzw7LPP4t133230/0WsI6YpJSUFxcXFmDdvns4+f39/g/eVk5MDQRC0elA1Xb58GQAQEhLS+AWRdWiWwZicCpRcqV9eyTNQlVv2zZ/rZ2CSQwlr5938oEycpbl7isPWM2NgRvQHMXgrLS2Fj49Pk8f7+fmhrq4ORUVFaNu2rda+//3vf81+3NWrV6OyshKZmZl4+OGHtfb95z//UffMWUJxcTFKS0sxcOBAvfsPHjwIoPmBGQBERkZi//79KCsrw8GDB7Fu3Tp8/PHHKCoqwr59+wyeJ9YR05SWlob8/Hy9+xojDmMaCry+/PJLeHt7G7xusjLNMhjJj6n+dfMEOj9UPyuv4ibLYpDq+Z+2X1XLzEHrmTHHTGIUCgV69+6tdxYeWZb4JS0OaTZFrCb/3Xff6ezTt82Q8+fPIzAwUCcoq6iowE8//dTs+zGFm5sbAODmzZs6+8QZpoD+YeKmeHt7Y+zYsUhPT4ePjw9Onz7dssYaQQzMxF5QTf/617/w888/48UXX9SaYUs2UnxFtdyOyM0TmPyZ9heuf7CqZ8TBvoBb4u7du3jiiSfwxBNP2O2QfJPV/g3pEKEqMFxdoQrYi6841FJNDMwkhsn/tjN79my4urpizpw5uHJF901eXFysVVdMTCBfsWKFVq5HYWEh/v73vzf7cUNCQnD79m38+uuv6m21tbVYsGABbty4YcqlNJuXlxc6d+6MnJwcrUT4u3fvYvLkyTh37hzc3NwQHh7e5H399NNPuHTpks72tLQ03LlzB3/605/M2vbGiM/TJ598gpqaGvX2/Px8zJ07F+3bt9ebU0dWJg5hisvtTPhEFZCFjWAQ1oTa2lp8/fXX+Prrr+2yIHmzq/0bIk4Cufy96jWkiHKY4IxDmUR/6NOnDz788EO88sor6NGjB0aNGoXQ0FB1Inl6ejri4+OxadMmAKpZoNOmTUNycjLuv/9+jB07FpWVldi9ezf+9Kc/4auvvmrW486ZMwcHDhzAI488gvHjx6N169ZIS0tDYWEhYmJikJaWZsGrVtUdmzdvHoYMGYKJEyfC3d0d//rXv9C9e3e4u7ujZ8+ezepZ2rBhA7Zv345Bgwahd+/e8Pf3R15eHvbv349OnTrhnXfeseh1iCorK5GXl4fBgwfj3LlzGDhwIEaOHInr169j7969kMlk+Pbbb3WGn8nKiq+ovlQ1E/xZn8xpNFrtvznEXLNv/ly/zUFWA2BgRqThpZdeQkREBN5//31kZGTg3//+N/z8/NC5c2e89tprmDp1qtbxW7ZsQffu3bFlyxZs3LgRnTp1wvz58zF+/PhmB2ajR4/GP//5T6xZswY7duyAp6cnhg0bhs8//xwrVqywxGVqmTNnDkpKSrB582Zs3boVoaGhSEhIwKhRoxAeHt7s/LK4uDhUV1fjxIkTyM7ORnV1Nbp06YLXXnsNf/7zn9GmjXWStnNzc1FTU4NBgwZhy5YtmDNnDjZs2ABPT0+MGTMGy5cv1zvRgaxIc9alm6cq4d+vk0N8qZJx9Fb7bw7N3MOSgvpeV0dgUvUzsriWFqgjEjWnwKw+u3fvFgAIf/vb33T2GSow6whYYNYKzh5UFQk9tVsQbl+2dWvskr0XmE37/Xrzi8o2xUrFh8UCs2m/X2/0OBaYJaJmmTZtGmQymU4RWUPE2aDGzMi0V+fOnYNMJoNMJtObJ0ctJBaLFZO0xaWVxFmX5FRanF/WkGYJDQvmmQV4ucPDrRVmbf8RhcVKiz0OhzIlRqFQQKFQ2GUyJ0lTRESEVo2wPn36NOu8nJwc9fmOLjAwUOv/qLEaaWSkhsOW47aobk/+jEGZk2pxfllD/sGqWZo7nrZonllHfw9smvIgpm49gdvlVeZpux4MzCQmISEBCQkJKC0t1VtAlMhYERERJgVXp06dQkhIiFMEKYGBgUbXSqNmKL4C/PeUKhAbshDIeA8o+l21j0srOT2T88v0sdLryaxtNoCBGRHpJa4hqk9MTAwA9iyRHsVX/igYqwR2jKufddk+QtVjdtjyE1qcgZeXl1WWuiPrY2BGREaLiYlRB2dEag3XuBSJpTDE4SYiMojJ/0REZB7i0kpDFtZvE4vG+gcDQT1UQRrXvHRqJlf8dxLsMSMiIvOoKFL966eRfK1Zn0ysPSXeJpPdvXsXU6ZMAQBs374drVu3tnGLmsfsMzIbEl+Ddow9ZkREZLyG6xNqlsForDfMP5hBmRnU1tbin//8J/75z3/a1Sx+cUbmpikPmndWo5VKZlgDAzOJ4SLmRCR5Yi6Z5vqE4jDmhO2qXjKiRph9dqNYMqO6QrXUlx0HZwzMJIaLmBOR5IlBmOaXoDiE5BlU33vBXDKyJjGHMfUlu17UnDlmRERkutSX6m+LgRhzycgAiyb+i6+7y9+rXpdFv9vl64+BGRERNU2sTwaoFo0GVDMuqyvqg7MJ27UT/Yk0WDzxH6h/3Yn5ZtP2Ax0iLPNYFsLAjIiIDCu+ogrENIvFihrmkrGaPzXC7EsxGeIfrArIkh9X/YjlWuwEc8yIiEg/Mck/+TFVUDZui2qNS+aPUQtYY1kjdIionwwg9vTaCfaYERGRfmKSvyiou+oLr2H+mJun6l8Galbj6emJsrIy9W3Sw057cBmYERGRcTSHhZjobxMymQxeXl62bkazFRYrce56mW0e3M6KzjIwIyKilmFARo0oLFZixF/ToayutWzif0OaRWftKM+MOWZERER2prKyEvHx8YiPj0dlZaWtm9MoMel//YQIHHo92rKJ/5o0i87aUZ4ZAzOJYeV/IrKphkstkSTV1NRg27Zt2LZtG2pqamzdnGYJa+dtvaBMZId5ZgzMJIaV/4nIZvQttUREVsUcMyIiZyQGXpp5N5qzMBsO/bA8BpnIotX+HRADMyIiZyP2jAGGk6KLztSXwZjwCdC+n90kT5N0WKXav4NhYEZE5Gz09YxV3KxfagnQXgPTrxODMjKJ1ar9OxAGZkREzqykQLVsjWYh2YZrYBK1kFWq/Tem6IxqKN4O/sBg8j8RkTMTe8+GLKzf5tcJ6PwQl14i+yfWMkt9STV8f+l7yU9sYY8ZEZEzK/njS8qvQU8CK/pLmqenJ65fv66+TQaIr+PL36uCs+THVIFaCwvOWnJCA3vMiIgcWfEV4GpOfS9B8RXtXLKM9wz3ivkHMyiTKJlMhrZt26Jt27aQyWS2bo60+Qer1nkVtaDgbICXOzzcWmHW9h9RWKw0UwO1MTCzgL///e8IDw+Ht7c3/P39MXz4cPzwww+2bhYRORtx9uXm6PphHEUUsHuSav+ET4CX01W9B+37ceiSzM7RSmV09PfApikPQlldi9sWujYGZhbQuXNnvP/++zh16hSOHz+O0NBQxMbG4uZN+1kSgogcgGb+WHUFcP6wdpK/XyegQ0R9z1jCCbtaU9CZVVZWIiEhAQkJCZJdkklSpTLEXDMzsPREBgZmFjB27FjExsYiNDQUvXv3xrp161BSUoLc3FxbN42InFH7CNWXUsZ7jR/HoUu7UVNTgw8//BAffvihZJdkEktlbJryoO1LZYh/eEz4xLbtaAaHCcx27NiBmTNnYsCAAZDL5ZDJZEhJSWn0nKysLIwaNQoBAQHw8vJCVFQUdu7cadZ2VVVVYfPmzQgICMD9999v1vsmImoWv06qL6WX0+3ii4kci81LZYj8g1XvBYlzmFmZixcvxqVLlxAUFIT27dvj0qVLjR6flpaG2NhYuLu7Y+LEifDz80NqaiomTZqE/Px8LFq0qEXt+e677/D4449DqVTi3nvvxcGDBxEYGNii+yQiMpnYG6Y5pMNcMrIgR8svsxaH6TFLSkpCfn4+bty4gVmzZjV6bE1NDWbMmAGZTIaMjAxs2bIF69atw6lTpxAeHo7ExEScPXtWffzixYshk8ka/WlowIAByMnJwfHjx/H4449j/PjxKCoqMvt1ExEZhblkZAWSyi+zMw4TmI0YMQIhISHNOvbIkSM4f/48nn/+efTv31+93cfHB0uWLEFNTQ2Sk5PV2xcsWICLFy82+tOQh4cHwsLCMHDgQCQlJcHFxUXrPomIbIa5ZGRhksovszMOM5RpjLS0NADAyJEjdfaJ29LT09Xb/P394e/v36LHFARBsjNniMjBFF9RzcgsOmPrlpCTk0x+mR1xysBMHKbs1q2bzr6AgAAEBQVpDWUa680338STTz6JTp064datW/jwww9RUFCAp59+2uA5lZWVWoFbaWmpyY9PRE7sao722pesS0ZkV5wyMCspKQEA+Pn56d3v6+uLgoICvfua4+rVq5g4cSKuX7+OwMBAREZG4rvvvkOvXr0MnrN27VosX77c5MckIickVvP3D66v6L9jnGrb5M8AzyC7WbiZjOPh4aFOo/HwkNZQYWGxEueul9m6GXbLKQMzS9u+fbvR57z11luYP3+++vfS0lIEB/PDlIgMEKv6A8DkVFVAJvaSTf4MCBthu7aRxbm4uKBLly62boaOwmIlRvw1HcrqWib+m8gpAzOxp0zsOWuotLTUYG+apcjlcsjlcigUCigUCtTW1lr18YnIzohV/QFVLplmRX/PINu0iZyemPS/fkIEIu8LZOK/CRxmVqYxxNwyfXlkt2/fRlFRkd78M2tISEhAXl4esrKybPL4REQkfVVVVVi4cCEWLlyIqirp1QsLa+ct7aCsQrrlq5wyMIuOjgYAHDhwQGefuE08hoiISGqqq6uxbt06rFu3DtXV1bZujv0QCyzvnlKfoykxThmYDR8+HF27dsXOnTuRk5Oj3n7nzh2sXLkSrq6uiI+Pt0nbFAoFevfujcjISJs8PhERkcPyDwYmbFcN/VfctHVr9HKYHLOkpCRkZmYCAH755Rf1NrFmWVxcHOLi4gAArq6uSEpKQmxsLAYPHoznnnsOvr6+SE1NxcWLF7Fq1Sp0797dFpeBhIQEJCQk2CTPjYiIqCXsYhkmiedgOkxglpmZiW3btmltO3bsGI4dOwYA6NKlizowA4ChQ4ciMzMTiYmJ2LNnD6qqqhAeHo6VK1di0qRJ1mw6EVHziIVjWZeMJIjLMJmHwwRmKSkpSElJMeqcqKgo7N+/3zINIiIylWZ9Ms1tiijVEIybJzBui23aRmSAOCNz2/QoaSf+S5xT5phJGXPMiJxY8RXg0veqAEwRpZ2cLJbHGLJQ9e9/c+r3lUgziZmcE5dhahmH6TFzFMwxI3JSmj1iIn3Jye0jVD1mGe/Vb9O8zSWYiOwaAzMiIinQLBgrKilQrXsJAI+9o/rXrxOQcEI716xh3hmXYHJ4Hh4eyM3NVd+WArtI/LcDDMyIiKTqvzn1wdqXc+u3+wdrB18MxJyOi4sLwsPDbd0MNSb+mw9zzCSGOWZETm7CJ8C0b+qHK908VdvcPDlMSZIlJv5vmvKg/ST+S7T6P3vMJIY5ZkROzq8T0CFCe7jSP1j1O8DeMQKgWpJpzZo1AIBFixbB3V0avVR2kfivWf0/4YTk3lPsMSMikiL/YFWAJn5pNBy+JKdWXV2N5cuXY/ny5VySyVgSr/7PwIyIiIici4Sr/zMwkxjmmBE5keIrkl1ImcgYdjsjs+iM5N6DDMwkJiEhAXl5ecjKyrJ1U4jIksS6ZQ0LyRLZGbuckSnmmaW+JLn3IAMzIiJbEOuWSTTPhai57HJGpjihZtwWyb0HOSuTiEgKJDp1n6i57GJGpib/4PqATELvP/aYERFZS/EV4GqO7rBJSYFq6j7rlBFZl2bpDIkMZ7LHjIjIGjTXwnTzVA2hiMQK/5M/Y0kMapbWrVvjxIkT6tu2UlisxLnrZTZ7/BYTS2fseFrVe2bE++/c9TIEeLmbffiWgZnEKBQKKBQK1NbW2ropRGQK8a9u8QNe/F3MKRuyUFXRXzOnRazwH9TDum0lu9WqVSubz94vLFZixF/Toayuta/E/4aMLJ0R4OUOD7dWmLc7Bx5urXDo9WizBmcMzCSGlf+J7JjYKwbUV+oXfxd7yPwa/EU+4RNVtX+xwj+RnRCT/tdPiEDkfYH2k/jfQh39PXDo9WhkXbyFebtzcLu8ioEZEZEkib1i4m2g/vei3/WfIy7BRGSEqqoq/P3vfwcAvPrqqzZdkimsnbfTBGWijv4euN3O2yL3zcCMiMgSKoq0h0gOr7BdW8jhVFdX44033gAAzJ49WzJrZVLLcVYmEZEl7J6imm0JAMOX2rYtRGQ3GJgREZnb8KXaRSslvC4fkansdhkmiWNgRkRkbs0JxCRU0JLIWHa5DFNTJPKeZGAmMVzEnMjBiQUtmXNGdswul2EyRGJFZhmYSQwXMSeyE8VXTPsQ9+ukKmhJ5ADsbhkmfcQisxJZM5OzMomIjNWwXpmx9ceCeqj+Qge4BBORFEgoD7RFgdnhw4dx5MgRHD9+HAUFBSgqKoKnpyfatm2L+++/H9HR0Rg9ejTuvfdec7WXiMiyGlbu16dhvTJjAzP/4PoCtCwqSyZo3bo1jh49qr5NjsPowKysrAwbNmzAli1bcPnyZQiCAED1wggMDIRSqURubi5+/vlnfPLJJ3B1dcWTTz6J1157DQ8//LDZL4CIyGxa2hNmDAZk1AKtWrVCTEyMzR6fMzItx6gcs02bNiEsLAyLFy+Gv78/Vq1ahSNHjqC0tBQVFRUoKCjAzZs3UV1djd9++w3btm3DhAkTcODAAQwZMgTjxo3DxYsXLXUtREQtI/aEVVcYrtSvT/EV4GpOfd0yIgfmkDMyJcSoHrM5c+Zg0qRJWLhwIcLDww0eJ5PJ0L17d3Tv3h1TpkyBUqnErl27sHbtWmzfvh1Ll7LYIhFJ3O4pzes1KykAkh+vH9rUR5z1Jd4maqHq6mps3rwZAPDyyy/Dzc3Nao8tzsjcNj3K/mdkSpBRgdlvv/2G0NBQox/Ew8MD06dPx9SpU1FQwL8oiUjihi9VlbNoTv6Y2Ms2ZCGQ8Z72vpI/8tX8OjGnjMyqqqoK/+///T8AQHx8vFUDM5FDzMiUIKMCM1OCMk2tWrVCSEhIi+6DiMjiTJmh1T6ivlcsqLvqdsZ7qn892zAgI6JmadGszMTERERERCAiIgL33XefudpERGR/GvaKJZxQ9aYxKCMiI7QoMFu5ciVkMhkAwNfXF/369VMHav3790d4eDhcXVkqjYjsRMMlWZqzRIvmMZoBmH8wAzJySA49I9McyzIVtyxlq0VR07Fjx7BixQp8++23KCkpQUZGBjIyMtTBmpubG3r16qUO1MSgzdfXt0WNdmQKhQIKhQK1tbW2bgqRcym+okr4d/OsH4pszgQALq1ETsRhZ2RqLsvUklI5xVdQtzmmRU1p0ZJMP/74Iw4cOIBVq1bhypUrKC8vx9mzZ7Fx40Z07doVVVVVOH36NLZt24Z58+Zh6NChCAgIaHGumiPjkkxENiIm8U/YDoQ8pH+JlobLMA3nDHNyLg61RqYmMy3LdP36VbjU3m1RU1rUY7Z+/XrExcVh0aJF6m2hoaF45ZVXMH36dEyaNAnl5eV4//33kZubi1OnTiE7Oxs///xzixpNRGQxYuJ/wwkAmsVnx21R/cullchJOeSMTDMsy1SqrEZL12FoUWBWUFCAZ599Vu8+uVyOTz75BL169cI333yD1157zeCxRESS13AZJoBlMMhm5HI5vvrqK/VtchwtGsrs3LkzfvzxR4P75XI5nnrqKWzdurUlD0NEZHuGkoKZ5E824OrqiieeeAJPPPGE1SbZFRYrce56mVUey5m1KDCbPHkyDh8+jL179xo8pqqqCufPn2/JwxAR2ZY4MYDISRUWKzHir+mYtzvH8RL/JaZFgdnChQvxwAMPYNKkSZgxYwYuXLigtT8vLw+ffvopgoJaPm5LRGQz4jAmk/1JIqqrq5GSkoKUlBRUV1db/PHEpP/1EyJw6PVox0r8l5gW9X96eHjg6NGjmDp1KrZu3YqUlBT06NEDwcHBuHnzJnJyclBbW4vXXnvNXO0lIrIdMyQHE5lDVVUVpk2bBgB49tlnrbYkU1g7bwZlFtbigWlvb2989tlnOHLkCBQKBQ4fPozTp08DUOWgzZ8/H3PmzGlxQ4mIzKL4Sn1Ffk1FZ/QfX3SmfuYlkZNy6KKyEmO2jMFhw4Zh2LBhAICSkhK4ubnB05MfZq+88go2bdqEDz74QL3gLBHZiFjyQpxd2ZC4riVQX3Ay9SXrtY9Ighy2qKxEtSjHzBA/Pz8GZQC++uorfP/99+jQoYOtm0JEQH2u2JCFuvvGbdGu+C2udynWLCNyUg5bVFaiLBKYEfC///0Pr7zyCrZv3261sX8iaiY/PeUtgrrrlr3wD1ZtJyLHLCorQUYFZqNHj260blljlEol1q1bh3/84x8mnd+UHTt2YObMmRgwYADkcjlkMhlSUlIaPScrKwujRo1CQEAAvLy8EBUVhZ07d5qlPdOmTcPcuXNx//33m+X+iIiIyPEZlWN25coVREVFISYmBlOmTMG4ceOaXJD85MmT2LFjB3bu3ImysjJs27atRQ02ZPHixbh06RKCgoLQvn17XLp0qdHj09LSEBsbC3d3d0ycOBF+fn5ITU3FpEmTkJ+fr7XMlLE2btyIsrIyvP766ybfBxFJhJhrJt4mIrIgowKz5ORknDx5EmvWrMH06dMxY8YM9OzZEw888ADuueceBAQEQKlU4tatWzh79ixOnjyJkpISuLi4YPz48Vi9ejW6dOlikQtJSkpCt27dEBISgnfeeQdvvfWWwWNramowY8YMyGQyZGRkoH///gCAxMREPPTQQ0hMTMSzzz6Lbt26AVAFfatXr2708QVBAAD89ttvWLlyJX744Qe4uHCkmMjuiblmQIsWNyYyJ7lcjj179qhvk+MwKjCLjIxEYmIiLly4gH379mHbtm1IS0vDjh07dI51cXFB3759ERcXhxkzZlg8AX7EiBHNPvbIkSM4f/48pk2bpg7KAMDHxwdLlizBxIkTkZycjDVr1gAAFixYgBkzZjTrvv/zn//gxo0bCAsLU2+rra3Fq6++iqSkJOTk5DS7nUQkEWLuGQMzkghXV1euP+2gjArMZDIZBEGAi4sLxowZg2nTpmHt2rV45JFHUFBQgJs3b8LDwwNt27ZFeHg4/Pz8LNXuFklLSwMAjBw5UmefuC09PV29zd/fH/7+/s2677i4OAwYMEBrW2xsLOLj49XFAPWprKxEZWWl+vfS0tJmPR4R/UGzPpm+tSuLrwAlBeZ5rJIr5rkfIpKeojOGP0eswKjArF27djh37pz691u3buG///0vevXqhV69epm9cZZy9uxZAFAPVWoKCAhAUFCQ+hhj6Qvi3Nzc0L59e61etIbWrl2L5cuXm/SYRE5Psz6Zm6d22YuG+4H6vDFD9cwMEc/LeE+75hmRldXU1ODzzz8HAIwdO9aiC5k7TXFZzdqF+j5HrMSoZ/LRRx/Frl270K5dO8TFxVmoSZZXUlICAAZ79Hx9fVFQYKa/rJvprbfewvz589W/l5aWIjjYNtE6kd3RrE+W8Z7qd80PVHG/yK+T6kP3v6eA3ZOa/zhivlljPXNEVlBZWYnx48cDAMrKyiwWmDlVcVnx/X35e1Vw1vBzxEqMeibfffdd/Prrr/jb3/6G9evXQyaT4R//+AdOnTqF/v37IyIiAhEREejUqZOl2muX8vPzmzxGLpczgZOopfTVJzPEP9i0nDH/YAZk5DTE4rLbpkc5R3FZUz8XzMiowOyee+7ByZMnkZaWhsOHD2P16tUoKyvDF198gc8//xwymQwAEBgYiIiICK1grXfv3ha5AFOIPWViz1lDpaWlNsuPUygUUCgUqK2ttcnjExERNcTistZjUt9nTEwMYmJisHr1aixcuBCvvfYacnJykJ2drf43IyMDhw8fBqCaNCClQEPMLTt79iwefPBBrX23b99GUVERBg0aZIumISEhAQkJCTYNDomIiMg2WjQovXXrVoSGhsLHxweDBw/G4MGD1fuqq6uRm5uLn376SXIlIqKjo7F27VocOHAAEydO1Np34MAB9TFEJDHFf8yGbGooUZxVBaiGJcw1G5OInEtzP3PMqEWBWXx8vMF9bm5u6N+/v1adMKkYPnw4unbtip07d2Lu3LmIiIgAANy5cwcrV66Eq6tro9dmSRzKJDJAnFkJGJ4tpTmrioioJUoKgOTHVbetOEPTcvNrrSwpKQmZmZkAgF9++UW9TaxZFhcXp55J6urqiqSkJMTGxmLw4MF47rnn4Ovri9TUVFy8eBGrVq1C9+62WbiYQ5lEBmjOrDQ0W0qcbSnOqgLqZ2oSkdGcplSGPs35zLEAhwnMMjMzddbhPHbsGI4dOwYA6NKli1aJj6FDhyIzMxOJiYnYs2cPqqqqEB4ejpUrV2LSJCOmzxORtDScVWXMTE0iO+Hu7o7k5GT1bUtwqlIZEuIwgVlKSgpSUlKMOicqKgr79++3TINMxKFMIiJqipubm8VTbpyuVIaJzN2ryFW2JSYhIQF5eXnIysqydVOIiIhYKsOAAC93eLi1wqztP6KwWGm2+2VgRkREZGdqamqwb98+7Nu3DzU1NbZujlPq6O+BTVMehLK6FrfN2GvmMEOZRCQxLZ1m3vD8iiLDx7EcBjmZyspKjB49GoBll2SixlmiN5HPpMQwx4wcQnNKWxhzPgDsnmL4OGMWI9dcfFwsryHeJiKyMQ5lSgxzzMghiNPMqytMW3eu4fni78OX6j+uuSZ8oh0oiosWW7FGEZE9KCxW4tz1Mls3wymxx4yI7IdnUMvO9+ukG4AxICPSUlisxIi/pkNZXctSGTbAwIyIiIjUzl0vg7K6FusnRCDyvkCWyrAyDmUSERERAO2isgzKbIOBmcQoFAr07t0bkZGRtm4KERE5GbGo7KYpDzIosxEOZUoM18okIqKmuLu7Y+PGjerb5saisrbDwIyIiMjOuLm5ISEhwdbNIAvgUCYRERGRRLDHjIiso6mVAIqv1Nc806zkX3TGsu0iskO1tbX47rvvAACDBw9Gq1atbNwiB1RyxSYPy8BMYlj5nxxSUysBNFbBP/Ul1b+aFfubolnRP6g7q/uTw7l79y6GDh0KQLUkk5eXl41b5EDEz4+M92zy8AzMJIbJ/+SQNCv0V9zUDczE/eO2qH4Xg7EJn6iKwgKqD8vmriIgVvTXd5uIqDHiZ0bFTVXv/e5JVn14BmZEJB1B3bV/9+sEdIio/92Y5Z00gzAGZERkDP9gm31uMPmfiIiISCIYmBERERFJBAMzIiIiQmGxEoXFSls3w+kxx4yIiMjJFRYrMeKv6VBWsyKArTEwIyIisjNubm5499131bdbSlwjEwA83FohgEsy2QwDM4lhHTMiImqKu7s7Fi5caPb7/WjKg+jT0Y8LmNsQc8wkJiEhAXl5ecjKyrJ1U4iIyMl09PdgUGZj7DEjIiKyM7W1tfjpp58AAA888ACXZHIgDMyIiIjszN27dxEVpVrmjEsyORYOZRIRERFJBAMzIiIiIolgYEZEROTECouVOHe9zNbNoD8wx4yIiMhJaRaWZf0yaWBgRkRE5KTEwrLrJ0Qg8r5AlsqQAAZmRMYovgJU3AQ82wD+wbZuTb3iK9q/m9q2htcn3q+++9M8FlDdFpUUaN+urtC/T1R0xrh2isfruy8iarab5VUAgLB23gzKGlN0xmqf+wzMJIaV/yWs+AqgiFIFGW6eQMIJaQRnmu0CTG9bw+ubnArsGKfa1/D+Gj5mY3ZPavx3kZtnfZDn5qn6V/xd5NlGtS/1pcbPJ3Jwbm5uSExMVN82RWGxErO2/8ghzMZofuZY6XOfgZnEJCQkICEhAaWlpfDz87N1c0hTxU1VIDJkIZDxnup3KQRmYrtE1RWmta3h9RWdqb/fhvfX8FgAGLcFCOquOq9h4DRui+pDTQzKxGM1af41mnBC9W/Da/APVu2ruKn9OBM+Adr3k8bzQWQF7u7uWLZsWYvuQxzG3DY9ir1lhoifOZe/V33eWOFzn4EZkbH8HPzL35jr0zw2qDvQIUL/cQ2DsMaOBRr/4PMP1t3v14lBGZGJ2rC3rHH+wdqpGhbGwIyIiMjO1NXV4fTp0wCAXr16wcWF1a8cBQMzIiIiO6NUKtGnTx8AXJLJ0TDEJiIiIpIIBmZEREREEsHAjIiIiEgiGJgRERE5IbG4LEkLAzMLWLZsGWQymdbPgAEDbN0sIiIiACwuK2WclWkh/fr1wzfffKP+3dTKzERERObG4rLSxcDMQlxdXXHvvffauhlEROSA3NzcsGDBAvVtU7G4rPQ4zFDmjh07MHPmTAwYMAByuRwymQwpKSmNnpOVlYVRo0YhICAAXl5eiIqKws6dO83SntOnT6N9+/YICwvD9OnTce3aNbPcLxERkbu7O9577z289957cHc3Prhifpl0OUyP2eLFi3Hp0iUEBQWhffv2uHTpUqPHp6WlITY2Fu7u7pg4cSL8/PyQmpqKSZMmIT8/H4sWLTK5LQMHDkRKSgp69uyJgoICJCYmYtiwYcjOzoZcLjf5fomIiFqK+WXS5jA9ZklJScjPz8eNGzcwa9asRo+tqanBjBkzIJPJkJGRgS1btmDdunU4deoUwsPDkZiYiLNnz6qPX7x4sU4yf8MfTY8//jieffZZ3H///Xj88cfx9ddf49KlS/jqq68scu1ERORc6urqkJ+fj/z8fNTV1Rl1rphftmnKg8wvkyCH6TEbMWJEs489cuQIzp8/j2nTpqF///7q7T4+PliyZAkmTpyI5ORkrFmzBgCwYMECzJgxw+S2BQUFoWvXrrh48aLJ90FERCRSKpW47777AJi+JBPzy6TJYQIzY6SlpQEARo4cqbNP3Jaenq7e5u/vD39/f5Mfr6SkBPn5+ejSpYvBYyorK1FZWan+vbS01OTHIyIiIvvkMEOZxhCHKbt166azLyAgAEFBQVpDmcZauHAhvvvuO+Tn5yMzMxNxcXG45557MGrUKIPnrF27Fn5+fuqf4OBgkx+fiIiI7JNTBmYlJSUAAD8/P737fX191ceY4sqVK5gwYQK6d++O5557Dh07dsThw4fh6elp8Jy33noLJSUl6p8rV66Y/PhERERkn5xyKNPSPv30U6PPkcvlnLFJRETk5Jyyx0zsKTPUK1ZaWmqwN83SFAoFevfujcjISJs8PhEREdmOUwZmYm6Zvjyy27dvo6ioSG/+mTUkJCQgLy8PWVlZNnl8IiIish2nDMyio6MBAAcOHNDZJ24TjyEiIpIaV1dXzJ49G7Nnz4arK7OSHIlTBmbDhw9H165dsXPnTuTk5Ki337lzBytXroSrqyvi4+Nt0jYOZRIRUVPkcjkUCgUUCgXzkx2Mw4TZSUlJyMzMBAD88ssv6m1izbK4uDjExcUBUP2lkZSUhNjYWAwePBjPPfccfH19kZqaiosXL2LVqlXo3r27LS4DCQkJSEhIsGmeGxEREdmGwwRmmZmZ2LZtm9a2Y8eO4dixYwCALl26qAMzABg6dCgyMzORmJiIPXv2oKqqCuHh4Vi5ciUmTZpkzaYTEREZRRAEFBUVAVCtLtNwacDGcAFzaXOYwCwlJQUpKSlGnRMVFYX9+/dbpkEmEruma2trbd0UIiKSqIqKCrRr1w6AcUsycQFz6XPKHDMp46xMIiKyFC5gLn0MzIiIiJwMFzCXLgZmRERETqCwWIlz18ts3QxqgsPkmDkK5pgREZG55RaW4NlN30NZXcv8MoljYCYxLJdBRETmVFisxLObvgcAbJsehbB23swvkzAGZkRERA5MTPjfNj0K0d3b2ro5DsmcJUgYmBEREdkZV1dXTJ06VX27OZjwb34BXu7wcGuFWdt/xMcvRqHslhLtWnifDMyIiIjsjFwuN7p2J5lfR38PbJryIKZuPYFnN32PcNlF7GrhfXJWpsRwrUwiIjKXwmIlCouVtm6GQ9PsiVwwskeL7489ZhLD5H8iImqKIAioqKgAAHh6eupdkqmwWIkRf02Hspqz/K0lOLDlkyrYY0ZERGRnKioq4O3tDW9vb3WA1pCY9A+AJTLsCHvMiIiIHNhHUx5En45+LJFhJ9hjRkRE5MA6+nswKLMjDMyIiIiIJIKBmcRwViYREZHzYmAmMQkJCcjLy0NWVpatm0JERERWxsCMiIiISCI4K5OIiMjOtGrVCs8884z6NjkOBmZERER2pnXr1ti7d6+tm0EWwKFMIiIiIolgYEZERORgCouVOHe9zNbNIBNwKFNiFAoFFAoFamu5thkREelXXl4Ob29vAEBZWRm8vLzU+zTXyORSTPaHgZnEcBFzIiJqCXGNzPUTIhB5XyCr/tsZDmUSERE5oLB23gzK7BADMyIiIiKJYGBGREREJBEMzIiIiIgkgoEZERERkURwViYREZGdadWqFUaNGqW+TY6DgRkREZGdad26Nfbt22frZpAFcCiTiIiISCIYmEmMQqFA7969ERkZaeumEBERkZUxMJOYhIQE5OXlISsry9ZNISIiiSovL4eXlxe8vLxQXl5u6+aQGTHHjIiIyA5VVFTo3X6zvMrKLSFzYmBGRETkIAqLlZi1/UcuXm5FAV7u8HBTzYz19XBr8f0xMCMiInIQ4gLm26ZHcZ1MK+no74FDr0cDANpV/I7SFt4fAzMiIiIHIQ5jtmFvmVWpg2D9o8tGYfI/ERGRA+AwpmNgjxkREZED4DCmY2BgRkREZGdcXFwQHR2tvq2Jw5j2jUOZFnL58mWMHz8eAQEB8PLyQmRkJAoLC23dLCIicgAeHh5IS0tDWloaPDzYO+ZI2GNmATdv3sQjjzyCxx57DIcOHYK/vz9+/fVXyOVyWzeNiIiIJIyBmQX85S9/wX333YfNmzert4WGhtqwRURERGQPHGYoc8eOHZg5cyYGDBgAuVwOmUyGlJSURs/JysrCqFGj1MONUVFR2LlzZ4vb8uWXX+KBBx7A008/jXbt2iEyMhKpqaktvl8iIiJAtSRT27Zt0bZtWy7J5GAcJjBbvHgxNm/ejEuXLqF9+/ZNHp+WloZHHnkE3333HZ555hm88sorKCoqwqRJk7BmzZoWteXixYv48MMPER4ejm+//RYTJkzAs88+i4yMjBbdLxERkaioqAhFRUW2bgaZmcMEZklJScjPz8eNGzcwa9asRo+tqanBjBkzIJPJkJGRgS1btmDdunU4deoUwsPDkZiYiLNnz6qPX7x4MWQyWaM/murq6hAZGYkVK1agf//+WLBgAUaPHq01tElERETUkMMEZiNGjEBISEizjj1y5AjOnz+P559/Hv3791dv9/HxwZIlS1BTU4Pk5GT19gULFuDixYuN/mi699570bNnT61tvXr1wuXLl1twhUREROTonDL5Py0tDQAwcuRInX3itvT0dPU2f39/+Pv7N/v+Bw0apNXjBgBnzpxpduBIREREzskpAzMxaOrWrZvOvoCAAAQFBekEVsZ47bXX8PDDD+O9997D2LFjcejQIXz55ZfqgFCfyspKVFZWqn8vLW3pMqhERERkbxxmKNMYJSUlAAA/Pz+9+319fdXHmGLgwIHYu3cvkpOTcf/99+Mf//gH9u7di4cfftjgOWvXroWfn5/6Jzg42OTHJyIiIvvklD1m1jB27FiMHTu22ce/9dZbmD9/vvr30tJSBmdERKSXi4sLBgwYoL5NjsMpAzOxp8xQr1hpaanB3jRLkcvlkMvlUCgUUCgUqK2tterjExGR/fDw8EBWVpatm0EW4JRhtphbpi+P7Pbt2ygqKtKbf2YNCQkJyMvL4xuOiIjICTllYBYdHQ0AOHDggM4+cZt4DBERkZQVFiuRW1iCwmKlrZtCZuCUQ5nDhw9H165dsXPnTsydOxcREREAgDt37mDlypVwdXVFfHy8TdrGoUwiImpKRUUFevfujdo6AZ7P/x2VcLN1k8hMHCYwS0pKQmZmJgDgl19+UW8TS1TExcUhLi4OAODq6oqkpCTExsZi8ODBeO655+Dr64vU1FRcvHgRq1atQvfu3W1xGUhISEBCQoJN8tyIiMg+CIKAS5cuAQCCq2rx6mO98MGRczZuFZmDwwRmmZmZ2LZtm9a2Y8eO4dixYwCALl26qAMzABg6dCgyMzORmJiIPXv2oKqqCuHh4Vi5ciUmTZpkzaYTERG1SEd/D1s3gczEYQKzlJQUpKSkGHVOVFQU9u/fb5kGERERERnJKZP/pUyhUKB3796IjIy0dVOIiIjIyhiYSQzLZRARETkvBmZEREREEuEwOWZERESOrrBYidvlVWgtq0H3nr2grKoFZLZuFZkTAzOJYR0zIiLSp7BYiRF/TYey+o/vh6fegwsAD7dWCPByt2nbyHw4lCkxzDEjIiJ9bpdXQVldiznDwtTb1k+IwKHXo1kuw4EwMCMiIrIjmkFYWDtvBmUOhoEZERGRnamrvourSbMRN+xPqKiosHVzyIyYYyYxzDEjIqImCUD1zcs4f1O1PBM5DvaYSQxzzIiIiJwXAzMiIiIiiWBgRkRERCQRDMyIiIiIJIKBGRERkcQUFitRWKy0dTPIBjgrU2I4K5OIyLmJFf4B6C0eG+DlDg/3VnD1a4d7fVtDJpOptrm1Uu8n+8XATGISEhKQkJCA0tJS+Pn52bo5RERkZWKFf/F2w8Cso78Hjvw5Fvhzvnqfp6cqiBP3k/1iYEZERGRn9AVfDMgcA3PMiIiIiCSCgRkREZGdUSqViIyMRGRkJJRKThJwJBzKJCIisjN1dXU4efKk+jY5DvaYEREREUkEAzOJUSgU6N27NyIjI23dFCIiIrIyBmYSw0XMiYiInBcDMyIiIiKJYGBGREREJBGclUlERGSHgoKCbN0EsgAGZkRERHbGy8sLN27csHUzyAI4lElEREQkEQzMiIiIiCSCgRkREZGdUSqViImJQUxMDJdkcjDMMSMiIrIzdXV1SE9PV98mx8HAjIiISCIKi5UoLK7vATt3vQwBXu7qfeT4GJhJjEKhgEKhQG1tra2bQkREVlRYrMSIv6ZDWV3/+T9vd47tGkQ2wRwzieGSTEREzul2eZVWUPbRlAexfkKE7RpENsEeMyIiIgnq6O9h6yaQDbDHjIiIiEgi2GNGRERkhzw9PW3dBLIABmZERER2xsvLC+Xl5bZuBlkAhzKJiIiIJIKBGREREZFEMDAjIiKyM3fv3sUTTzyBJ554Anfv3rV1c8iMmGNGRERkZ2pra/H111+rb5PjYI8ZERERkUQwMCMiIiKSCAZmRERERBLBwIyIiIhIIhiYEREREUkEZ2VKlCAIAIDS0lIbt4TU7pQBlQJQplT9e6cMkMLzI7ar4TZj29bw+sR/9d1fw2M1jzHUHsDw/ZlC83Gk8lwQtUDZnVLUVVZo/Q5Aa5u4vdSl/uu7tLSUMzMtTfy8aeqz5k4ZSv/4XBK/x40lE0w9kyzqwoULCA0NtXUziIiIyATnz59H165djT6PPWYSFRgYCAC4fPky/Pz8bNwa6yktLUVwcDCuXLkCX19fWzfHanjdvG5nwOvmdTuDkpISdO7cWf09biwGZhLl4qJK//Pz83OqF7TI19eX1+1EeN3OhdftXJz1usXvcaPPM3M7iIiIiMhEDMyIiIiIJIKBmUTJ5XIkJiZCLpfbuilWxevmdTsDXjev2xnwuk27bs7KJCIiIpII9pgRERERSQQDMyIiIiKJYGBGREREJBEMzIiIiIgkgoGZlZWXl2PHjh0YP348unfvDg8PD/j7+yM6Ohq7du0yeF5dXR02btyIvn37wsPDA23btsX48eNx9uxZg+dkZWVh1KhRCAgIgJeXF6KiorBz505LXFazZGRkYMGCBRg6dCj8/Pwgk8kQHx9v8PiYmBjIZLJGf7Zv3651TpcuXQweO2vWLAtfoX7GXndaWlqj1/yf//xH73n2/nxnZmbi9ddfx4MPPog2bdqgdevW6NmzJ958800UFxfrPccRnm/AMd7f+uTn5zf5Hm7VqpXWOaa+/qUmPj7e4DX07NlT7zmmvA6kxJTvN0d5vgHzvSdZ+d/KvvvuO0yZMgVt2rTB8OHD8fTTT+P69etITU3F888/j+PHj+ODDz7QOW/WrFnYsmULevfujTlz5uB///sfdu/ejQMHDuD48ePo3bu31vFpaWmIjY2Fu7s7Jk6cCD8/P6SmpmLSpEnIz8/HokWLrHXJalu3bsW2bdvg6emJzp07N7lAe3x8PGJiYnS2V1dXY+3atXBxccHw4cN19vv5+WHevHk62wcMGGBq01vE2OsWRUdH673+Tp066WxzhOf7mWeeQVFRER555BG88MILkMlkSEtLw7vvvovPPvsMx48fR7t27XTOc4Tn2xHe3/r4+/sjMTFR776TJ09i3759iI2N1bvfmNe/lL366qvw9/fX2hYUFKT3WGNfB1Jj6vcbYP/Pt1nfkwJZVU5OjvDJJ58IVVVVWtuvXbsmhISECACEEydOaO07cuSIAEAYPHiwcPfuXfX2Q4cOCTKZTBgyZIjW8dXV1UJoaKggl8uFn376Sb29tLRUCA8PF1xdXYUzZ85Y4Ooal5WVJeTm5go1NTXC999/LwAQpk6davT9/POf/xQACGPGjNHZFxISIoSEhLS8sWZk7HUfPXpUACAkJiY26/4d5fl+5513hKtXr2ptq6urE1555RUBgDB79mydcxzh+XaU97exRo8eLQAQPvvsM63txr7+pWrq1KkCAOHixYvNOt7Y14EUmfL95gjPt7nfkxzKtLJ+/frh+eefh5ubm9b2e+65BzNnzgQApKena+3bsmULAGDVqlVaBeuGDx+O2NhYZGRk4MyZM+rtR44cwfnz5/H888+jf//+6u0+Pj5YsmQJampqkJycbPZra8qAAQMQHh6uM3RhrKSkJADAiy++aI5mWZy5rtsQR3m+33zzTbRv315rm0wmw5IlSwDovi+kytjrdpT3tzGuXr2K/fv3o127dhgzZoytmyMJxr4OpMiU7zdHYO73JIcyJUR8Mbu6aj8taWlp8PLywsMPP6xzTmxsLL755hukp6eje/fu6uMBYOTIkTrHi9vs9c1RUFCAAwcO4N5778UTTzyh95jKykps27YNhYWFCAgIwKBBg9CvXz8rt7Tlzp49iw0bNqCiogIhISF49NFH9Q6BOPLzDRh+X4js/fl2xvd3SkoKamtr8cILL+h8iYua+/qXun379uHOnTuQy+Xo27cvYmJi9Abtxr4O7E1T72N7fr7N/Z5kYCYRtbW1+PjjjyGTyTBixAj19vLycvz3v/9Fnz599L6Zu3XrBgBayaHibXGfpoCAAAQFBdlNMmlDycnJqKurQ3x8vME3+LVr13SSrR977DFs377dbt7oALBz506txFEPDw8sX74cCxcu1DrOkZ9vQJWzBej/0APs+/l2xve3IAjq57SxXu/mvv6l7v/9v/+n9Xv37t2xa9cuPPDAA+ptprwO7Imh7zdN9vx8m/s9yaFMiViyZAl++eUXTJs2DX369FFvLykpAaBKcNbH19dX67jmnqN5vL0QBEHdHWzoA3369OlIS0vDjRs3UFpaiv/85z94/PHH8c033+DJJ5+EYAcrkLVt2xbvvfceTp8+jfLychQWFmLHjh0IDAzEG2+8gY8++kjreEd9vgEgJycHy5cvR7t27fDGG2/o7Lf359sZ39/p6ek4f/48HnnkEb2zE419/UtVdHQ0PvvsM1y5cgVKpRKnT5/GvHnzcP78eYwcORJXr15VH2vK68CeGPp+Axzj+Tb7e9LcSXDOok2bNgKAZv8cPXrU4H199NFHAgChf//+wp07d7T2FRYWCgCEhx9+WO+5GRkZAgDh5ZdfVm979NFHBQDC2bNn9Z7TtWtXwd3d3fiLFsx33aYk/x86dEgAIERHRxvV5traWuGRRx4RAAhfffWVUeeKbHndol9++UVwd3cX7rnnHqG2tla93VGf7wsXLggdO3YU5HK5cOTIkWafZ0/Pt9Te34aY8/Nu8uTJAgAhOTnZqDYYev1bkjmvW7Ro0SIBgLBgwQL1NlNeB5Zkre+3xtji+TaVud+THMo00XPPPYc7d+40+/h7771X7/bk5GTMmjUL999/Pw4ePAhvb2+t/WIEbijaFqfia0bqzTnHUGTfFHNdtynEpP8ZM2YYdZ6LiwumTZuGzMxMHDt2zGBuWmNsed2iPn36YODAgfjuu+9w7tw5da6JIz7fly5dwtChQ3Hjxg189tlnGDp0aLPPtafnW2rvb0PM9f9RXFyMzz77DL6+vhg/frxRbTD0+rckS7wOXnzxRaxZswbHjh1TbzPldWBJ1vp+a4wtnm9Tmfs9ycDMRIZqsRhj69ateOmll9C7d28cPnwYbdq00TnGy8sL7du3x8WLF1FbW6uTf6BvbFszH+HBBx/UOv727dsoKirCoEGDTGqzOa7bFLdv38bnn38Of39/PP3000afL+YaVVRUmPT4trruhvRdh6M93/n5+Rg6dCiuXr2KvXv3YvTo0Ubfh70831J7fxtirv+PnTt3QqlU4oUXXoCnp6fR57f0eTWWJV4H+q7BlNeBJVnr+60p1n6+TWXu9yRzzGxk69atmDFjBnr27IkjR46gbdu2Bo+Njo5GeXm51l9Yom+//VZ9jObxAHDgwAGd48Vtmsfbgx07dqCyshKTJk2Ch4eH0ef/8MMPAFSV4u1VTU0NfvrpJ8hkMnTu3Fm93ZGe7/z8fMTExKCwsBC7d+/GU089ZdL92NPz7Uzv7//7v/8DYHyvN2D49W9vDL02jX0dSJkx32+G2NPzbfb3pLnGWKn5kpKSBJlMJvTq1Uu4du1ak8drFh6srKxUb2+sAGXXrl0FuVwuZGdnq7drFrv7/fffzXY9pjA256hfv34CAK3ifQ39+uuvwu3bt3W2f/fdd0Lr1q0FuVwuXLp0ycQWm0dzrvv48eNCXV2d1rbq6mph3rx5AgDhscce09nnCM/3xYsXhZCQEMHV1VWn6Kg+jvJ8O+L7W5/s7GwBgNC3b99GjzP29S9F//3vf4Vz587pbC8oKBB69uwpABA+/fRTrX3Gvg6kytjvN0d4vs39npQJgsSnLTmYI0eOYMSIERAEATNnztQ7Nh8REYG4uDitbS+99BKSkpLQu3dvPPHEE+qlOlq3bq13qY6jR48iNjYWcrkczz33HHx9fZGamoqLFy9i1apVePvtty15mXplZmaq88Ru3LiBr7/+GqGhoXjkkUcAAD179sSf//xnnfN+/PFHDBgwAA888AB+/PFHg/e/bNkyvPvuuxg+fDi6dOkCuVyO3NxcHDhwAC4uLti0aZNJf6m3lLHXLa7/OGjQIHTs2BHFxcXIyMjA77//js6dOyMjIwMhISFaj+EIz3eXLl1w6dIl/OlPfzK4TM+yZcu0bjvC8w04xvu7KXPmzMHGjRuxYcMGzJkzx+Bxprz+pSYtLQ3Dhg1TzzwNDAxEfn4+vvrqK5SXl2Pq1KlITk6GTCbTOs/Y14HUmPL95gjPN2Dm96SZA0dqQnJycpMzXPT9dV1bWyts2LBBCA8PF+RyudCmTRvhmWeeaTQK/+GHH4THHntM8PPzEzw8PIQBAwYIO3bssODVNa6pazc021JckufDDz9s9P7T0tKE8ePHC2FhYYKPj4/g5uYmdOrUSZg4caLwww8/WOCKmsfY637nnXeEmJgYoUOHDoK7u7vg6ekp9O3bV3j77beFW7duGXwce3++m3pfNPy4cpTnWxAc4/3dGKVSKQQEBAhyubzR17AgmP76l5LLly8LM2bMEPr27SsEBAQIrq6uQps2bYRHH31Up6dMkymvAykx5fvNEZ5vkbnek+wxIyIiIpIIJv8TERERSQQDMyIiIiKJYGBGREREJBEMzIiIiIgkgoEZERERkUQwMCMiIiKSCAZmRERERBLBwIyIiIhIIhiYEREREUkEAzMiIiIiiWBgRkRERCQRDMyIiIiIJIKBGRGRlVRUVGDFihXo1q0b5HI5QkND8cEHH+D48eOQyWRYsmSJrZtIRDbmausGEBE5gzt37mDYsGE4efIkhg8fjnHjxuHMmTOYO3cuYmNjAQARERG2bSQR2RwDMyIiK4iPj0d2djZ27dqFiRMnqrevWrVK3VPGwIyIZIIgCLZuBBGRIzty5AiGDx+O+Ph4JCcna+0rKChAcHAwfH19UVxcDJlMZqNWEpEUMMeMiMjCFAoFZDIZ3n77bZ19gYGBAIB+/foxKCMiBmZERJZ26NAh9OjRA2FhYTr7rl69CoDDmESkwsCMiMiCiouLUVpaiuDgYL37Dx48CICBGRGpMDAjIrIgNzc3AMDNmzd19t29exfvv/8+AKB///5WbRcRSRMDMyIiC/Ly8kLnzp2Rk5OD3Nxc9fa7d+9i8uTJOHfuHNzc3BAeHm7DVhKRVDAwIyKysPnz56Ourg5DhgzB7NmzMW/ePPTq1Qt37tyBu7s7evXqBXd3d1s3k4gkgHXMiIgsbM6cOSgpKcHmzZuxdetWhIaGIiEhAaNGjUJ4eDjzy4hIjYEZEZGFubi4YOnSpVi6dKnW9j179gBgfhkR1eNQJhGRjZw6dQoAZ2QSUT0GZkRENpKTkwOAgRkR1eOSTERENtKpUye4uroiPz/f1k0hIolgYEZEREQkERzKJCIiIpIIBmZEREREEsHAjIiIiEgiGJgRERERSQQDMyIiIiKJYGBGREREJBEMzIiIiIgkgoEZERERkUQwMCMiIiKSCAZmRERERBLBwIyIiIhIIv4/epRpeezW2SAAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot histograms of q\n",
    "hist_range = -200, 0\n",
    "bins = 400\n",
    "\n",
    "plt.hist(qb, range=hist_range, bins=bins, histtype='step', density=True, label=r'$f(q|b)$')\n",
    "plt.hist(qsb, range=hist_range, bins=bins, histtype='step', density=True, label=r'$f(q|s+b)$')\n",
    "\n",
    "plt.axvline(med_q_sb, color=\"black\", linestyle=\"dashed\", label = r'median$[q|s+b]$')\n",
    "\n",
    "plt.legend(loc='upper left', frameon=False)\n",
    "\n",
    "plt.xlim(*hist_range)\n",
    "plt.xlabel(r'$q$', labelpad=5)\n",
    "\n",
    "plt.yscale(\"log\")\n",
    "plt.ylabel(r'$f(q)$', labelpad=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "fb2a1953",
   "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": 8,
   "id": "edb93068",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of toys passing the median = 80\n",
      "Expected p-value = 8e-06\n",
      "Expected significance = 4.314451021808664\n"
     ]
    }
   ],
   "source": [
    "stat = (qb < med_q_sb).sum()\n",
    "print(f'Number of toys passing the median = {stat}')\n",
    "pval = (qb < med_q_sb).mean()\n",
    "print(f'Expected p-value = {pval}')\n",
    "significance = stats.norm.isf(pval)\n",
    "print(f'Expected significance = {significance}')"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}