{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5acb6c5e",
   "metadata": {},
   "source": [
    "======================== Import Packages ==========================="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6d06bd06",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys, os, pdb, glob\n",
    "import numpy as np\n",
    "from astropy.io import ascii\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as mtick"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebac2bbf",
   "metadata": {},
   "source": [
    "========================== Define Fuctions ========================="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "139b88b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_twosample(file, ax, color, region, xpos, ypos):\n",
    "\n",
    "    \"\"\"\n",
    "    PURPOSE:    Plot the results of the two-sample tests\n",
    "                Ouput from calc_twosample_test.py\n",
    "\n",
    "    INPUT:      file = file name of two-sample test results (str)\n",
    "                ax = plot axis (matplotlib ax)\n",
    "                color = color for plotting this region (str)\n",
    "                region = region name (str)\n",
    "                xpos = x position to plot region name (float)\n",
    "                ypos = y position to plot region name (float)\n",
    "\n",
    "    OUTPUT:     Updated plot axis (matplotlib ax)\n",
    "                Prints median probability to terminal\n",
    "    \n",
    "    \"\"\"\n",
    "\n",
    "    ### PLOT REGIONS NAME\n",
    "    plt.text(xpos, ypos, region, size=15, color=color, alpha=0.8)\n",
    "\n",
    "    ### PLOT TWO-SAMPLE RESULTS\n",
    "    t = ascii.read(file)\n",
    "    xv = np.sort(t['col1'])\n",
    "    yv = np.arange(len(xv)) / float(len(xv))\n",
    "    ax.plot(xv, yv, color=color, linewidth=5, alpha=0.5)\n",
    "\n",
    "    ### PRINT MEDIAN PROBABILITY\n",
    "    ind = np.where(yv == 0.5)\n",
    "\n",
    "    print(\" \\tMedian p_phi for \" + region + \" = \" + str(xv[ind][0]))\n",
    "\n",
    "    return ax"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c055da3",
   "metadata": {},
   "source": [
    "============================== Code ================================"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c60b1966",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " \tMedian p_phi for Upper Sco = 9.466132e-07\n",
      " \tMedian p_phi for Taurus = 0.7125768\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAukAAAI4CAYAAADXpneBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABqCUlEQVR4nO3deXxU5d3///csWSYbIQuQyL6DiBEE3Cqi4oLVVgWXatGvWGxrN2+t/mxvt6q9tVTt4lJoFVxqrWKs+0JlUWRf3dhEQoAEQkhCtpnJLOf3x5iRkJkhy8xkJvN6Ph48COec68yVMyfDm4vPuS6TYRiGAAAAAMQMc1d3AAAAAEBLhHQAAAAgxhDSAQAAgBhDSAcAAABiDCEdAAAAiDGEdAAAACDGENIBAACAGENIBwAAAGIMIR0AAACIMYR0AAAAIMYkREivr6/XPffcowsuuEA5OTkymUxasGBBm9vX1NRo9uzZys/PV3p6uqZMmaINGzZErsMAAABIaAkR0isrK/W73/1OW7Zs0Yknntiutl6vVxdddJFefPFF/exnP9Mf/vAHVVRU6KyzztKOHTsi1GMAAAAkMmtXdyAaCgoKVF5erj59+mjdunWaMGFCm9suXLhQK1as0CuvvKLp06dLkq644goNHz5c99xzj1588cVIdRsAAAAJKiFG0lNSUtSnT58OtV24cKF69+6tyy67zL8tPz9fV1xxhV5//XU5nc5wdRMAAACQlCAhvTM2btyocePGyWxueakmTpyoxsZGbd++vYt6BgAAgO4qIcpdOqO8vFxnnnlmq+0FBQWSpLKyMp1wwgmdOn95eXmr7dXV1dqyZYtOOukk2Wy2Dp8fAAAAkWG321VSUqLzzz9feXl5YT03If0Y7Ha7UlJSWm1PTU317++MuXPn6r777uvUOQAAANB1XnjhBV1zzTVhPSch/RhsNlvAunOHw+Hf3xk33XSTLrnkklbbN23apFmzZumFF17QqFGjOvUaANrOMAxVNFSosrFS1Y5qOdwOVTZU6rDzcFd3LTYcqpS8RvvamE1SbnhHmBLNuIJxmnBc2yc9iFvbtkleb/vamM3SiBGR6Q8QwFcLD8txyC1JKqncoXuKb9bAgQPD/jqE9GNonhnmaM3bCgsLO33+5tKZQEaNGqVx48Z16jWARGYYhhxuh+xuuw41HlKts1Z2t112l10ur0tur9v/q6SmRI2uRl/D1G9+SUrLS1Oa0rrse4gp+00dC1EdfHgfPkP7D9W4wQnwd4HZLHk87WtjsUhFRRHpDhCI5aNKOZNb3qeRKE0mpB9DUVGRPv74Y3m93hYPj65evVppaWkaPnx4F/YOSGxew6sD9Qd02Hm4Rdh2e93ac3iPap21Kq8vl9vr7uquog3MJrMsJotcXpdMMmlY7jD/NrPJ7Pva/O3XkmSSyfe7yRT0z6H2RevPJpn8/TaZjvj6iO1H7ztyf0ZyRhiuMIBwcDva+b+JHURIP0J5ebkOHz6sIUOGKCkpSZI0ffp0LVy4UMXFxf550isrK/XKK6/o4osvDlivDiByvIZXZXVlemPbGzrUeEgeo52jbgiLZHPSEb+sSjb5vk464utks9X3uzVFycPHKNmSrGRLspIsSf6vj/xlMVn8ARcAYpHhNeRxtPN/EzsoYUL6448/rpqaGpWVlUmS3nzzTe3du1eS9POf/1w9evTQnXfeqWeffVa7du3y1xZNnz5dp5xyiv7f//t/+vLLL5WXl6cnn3xSHo+HBz6BMDIMQ4fsh1RWV6ZaZ63qnHWqb6pXXVOd6px1cnqcavI0MSr+jSNHl4/+1Vw/f3z+8SFHbo/cbjFZZEuyKS0pzRekzUlKsiS1/v2zL2T1mtoXpi0WqbAoMhcCAKLI44zOKLqUQCH9j3/8o3bv3u3/c3FxsYqLiyVJ1157rXr06BGwncVi0TvvvKNf//rX+stf/iK73a4JEyZowYIFGsGDKkCnVDZWalnJMpXXl6uysbKruxN2GckZykzOVHZqtqxmq5IsSUqxpCjHlqMkS1KLgByonCPJnKQUa4p/pPnIIN5lzEkS/3sBIEFFq9RFSqCQXlJScsxjFixYoAULFrTa3rNnT/3jH//QP/7xj/B3DEgQTZ4mHWw4qIqGClU0VOjzis9V11TX1d1ql2RLsgozCzWgxwDZkmwBSzaaR6HTk9NlNSfMRywAJASPPTqlLlIChXQA0WMYhkpqSlRSU6Iqe5X21u5VtaO6q7vVZlazVdmp2fIaXk0dPFVpSWlKtaYqLy1PFrOlq7sHAOgijKQDiDsuj0vbD23XgYYD2lq5VRUNFV3dpRZMMmlA9gDl2nKVmZIpm9Xmf4jRarb6ylHMSbIl2ZSXlte1JSUAgJgUrYdGJUI6gHZwup3aU7tHBxsOqsHVoPqmejU0NeiQ/ZCq7FVd2rc+GX3Uv0d/f+C2mCxKtiQrPz1fg7IHMQIOAOg0t52RdAAxYm/tXm05uEUlNSUqry+X14jeKMLRTDL5675tSTb1Su+lwT0H66Q+JzF1HwAg4hhJB9Blahw12lW9S3tr92p9+fqov35+Wr7MJrMG9RykAT0GKDMlU5nJmUpLSpPVbCWMAwC6DCPpAKKm2l6t3Yd3a2/tXu2q3qVD9kNRff0kc5L6ZvXViX1O1JheY5gRBQAQsxhJBxBRLo9LO6t3annpcu2t3RuV1zTJpLy0PPVK76Ve6b2Un56v7NRs5aflK8mSFJU+AADQGczuAiAsDMNQrbNWja5GOdwOVTRUaEfVDpXUlERl5U6b1abT+5+uvll91Sejj1KtqRF/TQAAIoV50gF0So2jRk+ufVJurztqD3qaZFKOLcc/Wj6452ANzB5IDTkAoNtgJB1Am1Tbq1XRUKFD9kM62HBQVfYqVTRUyO62R+w1TTLpxD4nqnd6b6UnpysjOUPpSenqaeupZEtyxF4XAICu5iGkAwjE4/WooqFCmw9s1tbKrapx1ETttU/odYLG9h6rYbnDovaaAADEEjflLgCaOd1OrdizQlsrt+pAw4GovGafjD4amjNUvdN7K8eWo4LMAlbgBAAkNMMwGEkHEp3X8GpD+QZtq9ymHVU7Iv56aUlpOqnPSerfo7/6ZvVVenJ6xF8TAIB44nUZMryEdCDhVNmr9HnF59pfv19fHvwy4q9ns9o0uOdgDe45WGN6jVGKNSXirwkAQLyK5kJGEiEd6HJOt1OPrnxUTo8z4q91Qq8TlJuWqyE9h+i4rOMoYQEAoI2iWeoiEdKBLmMYhkpqSvTs5mfDfm6TTMpNy1X/Hv01uOdgjcobJYvZEvbXAQAgUbijuNqoREgHoq6hqUEfl36sVXtXhe2cA3oM0Kj8UeqV3kvZqdnKSM5gOkQAAMLIQ7kL0D0ddhzWij0rtL58fadX+8xLy9PA7IEalD1Io/NHs2AQAAARxkg60M1UNFToxc9e7PSc5nlpebpg6AUa3HMwteQAAEQZNelAN/HpgU/1wc4PVN9U3+FzDO45WKPzR+v4/ONlS7KFsXcAAKA9ormQkURIB8LC4/Xo6+qvtbd2rw42Huz0FIrThk1TUZ8i6soBAIgRjKQDccQwDC0vXa4Pd30YlvON6TVGFw+/mDnLAQCIMYR0IA443A59tPsjrdizotPnspqt6t+jv84ZdI6OyzouDL1rrayuTJf86xJJ0rrZ64Ied+/Se/XW9rc0e/xszR4/OyJ9iUXrytbp5S9e1mcVn6naXi1bkk09U3tqWM4wjS8cr2nDpikjOaOruwkA6EKUuwAxrtZZq0dXPhqWc108/GIV9SliDvMu9Pf1f9fc9XMlSYN6DtKY/DGymq3afXi3lpQs0Ye7PtSovFE6ofcJXdxTAEBXcjOSDsQmu8uut7a/pS8OftHpc008bqIuHHohUyd2sS0Ht2jehnmymq166NyHdNbAs1rsP9R4SO/seEeZKZld00EAQMzwMJIOxJbyunL/SGtnXV90vQb0GEA4jxFLSpbIMAxNHTK1VUCXpNy0XP3wxB9Gv2MAgJjDSDoQIw7UH9DCLxfqYOPBDp+jZ2pPFWYWanzheA3uOTiMvYuu2W/O1obyDXrj6je0af8mvfjZi9pVs0tpSWk6re9punnizeqV3qtFm3nr52ne+nm6Z/I9GtRzkP627m/6vOJzeQ2vxvQao5+c/JOgJSS7qnfp2c3Pas2+Nap2VCszOVMnF56sG8fd2Oo6vrntTd237D7NHj9bFwy9QH9b9zetK1unake15kydEzB8N6u2V0uSetp6tvua2F12vfzFy/rvrv9qd81ueQ2vemf01sTCibr6hKvVv0f/Fse/s+Mdvfrlq9pRtUMew6N+Wf10/pDzdc3Ya5jFBwDigIfFjICu43Q7tbRkqdaVrZPL6+rQOdKT0nX56MvVN6tvtwtfL3z6gl758hWd1OckTR4wWZ9XfK63d7yttWVrteD7C1oFdck3X/zvl/9e/bL66bR+p2lv7V6t2bdGG/dv1GPnP6ZT+p7S4vilJUv1mw9/oyZPk4bnDtcJvU7QgYYDWvT1In20+yP95cK/aFzBuFavU1JToh++9kP1SOmhkwtPVq2zVlZz6I+43hm9JUkf7vpQ1xddrxxbTpuuQ2VjpX769k/1dfXXykrJ0smFJyvJnKR9dfv06pZX1a9HP/3ghB/4j//9x79X8ZZiJVuSNaFwglKtqVpfvl5PrH1CH5d+rCcvelKp1tQ2vTYAIPq8bkNeFyPpQJdYsWeFPtj5QYfbD8sZpguGXqDctNww9iq2vLrlVf3p/D/p9P6nS5LcXrfuW3qf3v3qXT28/GE9cv4jrdq8tvU13XDSDfrJyT/xl/ks/HKhHlr+kO5deq9ev+p1/5STZXVlumvJXbKarfrTBX/SxOMm+s+zYs8K/c/7/6O7ltyl/1z5HyVZklq8zgc7P9AVx1+h2067rc0rsl449ELN3zRfB+oP6PsvfV9nDzpbRX2KNCpvlIblDgt6nruX3K2vq7/W1MFTddfku5SWlObfV1ZXpoamBv+fF+9arOItxcpPz9fc7871j7DXN9XrV+/9Spv2b9Lf1v1NvzrlV23qMwAg+txRHkWXJNYWR8Jze92a88mcTgX0H479oa4Ze023DuiSNHXwVH9Al3zTR9522m1Ktabqo9KPdKD+QKs2BZkFumn8TS3q8KePnq4xvcaosrGyxRzz//rsX7K77PrZxJ+1COiSdFq/0zR99HQdqD+g5aXLW71OT1tP/WLSL9oc0CXpuKzj9Nj5j6l3Rm81uhr11va39MBHD+ia4mt0znPn6KHlD6mysbJFmy8qvtCafWuUY8tpFdAlqTCzUMNyh/n//NLnL0mSfjTuRy1KYDKSM3TH6XfIZDKpeEuxmjxNbe43ACC6oj1HukRIRwL7YOcH+tOqP+mBjx5Qg6vh2A2CuPOMOzUkZ0gYexa7zhtyXqttPVJ76JS+p8gwDG3av6nV/rMHnh1wisnzh5wvSS3arNq3SpI0ZeCUgK9/Up+TJCngDDsTCyd2qGRk4nET9Z8r/6M5U+fo8lGXa2TeSFnMFtU567Twy4X6was/0O6a3f7jV+9b7e//0QH9aG6vW59VfCbJN2p/tGG5wzQsZ5gaXY3aVrmt3X0HAERHtOdIlyh3QQL6aPdHWrxrcafOMbb3WJ3S9xQVZhaGqVeRZVLbZpMxjNAjBQUZBQG3N1+HQA/ZFmQeo03Dt23K6sokSRf+s3WgPVKNo6bVtj4ZfUK2CSXJkqQpg6ZoyiDfPw7qnHX6YOcHemLtE6qyV+nhTx7Wkxc9KUn+/y3om9X3mOc97Dgsl8el7NRs2ZJsAY8pyCjQ9kPbO/WAMgAgsrpiJJ2QjoRRUlOi5zY/J6/RsX8Np1hS1K9HP1095uq4W3zoyBFmh9sRdMTZ4XZIkmzWwIEy0pr/kfDd4d8NedyYXmNabWuuaw+HzJRMXT76cuWn5+t/3v8frStbF/K6dQbTcQJA7COkAxGyrGSZlpYslaH2/5CNLxiv4bnDNbjn4FYPK8aLHqk9lGJNkdPt1L7afUHLc/bV7ZP07awnRyuvL29Rb+3fXlcuScpPyw+6L9C5JCk//ds2vdJ7aW/tXt1yyi3qkdojxHcUHScXnixJ8hpe1TnrlGpN9V+bvbV7j9m+R2oPJVmSVOOokd1lDzia3vy/B4GuHQAgNnRFuQs16ejWGpoadO/Se32L1nQgoN915l26eMTFGpE3Im4DuiSZTWad2PtESQr40KXkK+PYfmh7i2OP9t+v/9tqW62zVqv2rZLJZNKJfVq3W1KyJOD/XjQ/qFvUp8i/bdJxk/xtouFY5T3NQTzJkqTs1GxJ3/bx/Z3vq9HVGLK91WzVCb18c8EHejB5Z9VO7ajaobSkNI3IG9He7gMAooQHR4Ew2lu7V3NWzOlQ29P6naZ7Jt8Td2UtoVw95mpJ0oLNC/R5xect9tU31eu+ZffJa3g1ZeCUoCPpH+z8QCv3rPT/2eP16JEVj8jusus7/b8TsC68rK5M89bPa7GteEuxPj3wqXJsOTp70Nn+7deOvVYp1hT9adWfAj430ORp0odff6iKhoq2f+MhPLXuKf151Z8DjopXNFTowY8flCSd2f9M/z/Sju91vE4uPFlV9io9+NGDsrvsLdqV1ZXpq6qv/H++8vgrJUlz18/Vvtp9/u2NrkY9/MnDMgxDl426rNvNqQ8A3UlXTMFIuQu6HY/Xo6c3Pu0vI2iPEbkjdP7Q89u8qE08+c6A72jmiTP13ObndMPrN2hMrzEqzCxUfVO9Nu3fpPqmeg3JGaL/74z/L+g5Lht1mX7x3i90Up+TlJeWp88rPldZXZny0/N1++m3B2xz6chLtWDTAi3etVjDcoZpT+0efXnwS1nNVt171r0t6rz79ein35/9e/128W91+6Lb1a9HPw3KHqRUa6oONhzU1kNbZXfZ9eLlLwZcOKm97C67/vX5v/T8p8+rf4/+GtxzsJItyapoqNDnFZ/L7XWrX49+uvW0W1u0+92U3+knb/9E7+98Xyv3rlRRnyIlW5K1t3avth/arl+d8isNzRkqSTpn8Dm6bNRlKt5SrCsWXtFiMaNqe7VO6H2Cfnzyjzv9vQAAIsdtpyYd6DDDMLSubJ3e3vF2u9vePOHmFrXR3dUvJv1CJxeerFe/fFWfH/xcXxz8QqnWVA3MHqhzBp2jGaNnBJ2FRPKNdI/KG6V/ff4vfV7xuWxJNk0bNk0/m/izoKF5bO+xunjExfrbur/p49KPJfmmPfzxyT/W2N5jWx0/eeBkvTT9Jf3z039q9b7VWrV3laxmq/LT8/Wd/t/R2YPO1qDsQWG5HrPGzdKo/FFatXeVth/aro37N6q+qV7pSek6Pv94TR44OeA16ZXeS89f+rxe/OxFfbjrQ63et1oWk0W90ntp+ujp+k7/77Q4/jff+Y2K+hRp4ZcLtaF8g9xet/pm9dXVY67WNSdcE9aHXgEA4ceDo0AHGIah3Yd3a8GmBe1q1y+rn6aPnh4TDyhG02n9TtNp/U7rcPuLR1ysi0dc3K42Y3uP9U9h2BZ9s/rqjjPuiFh/mmWnZmvasGmaNmxau9umJaXpxnE36sZxN7bp+I6+DgCg63kodwHazuP1aPOBzVqxZ0WrVSGP5XsjvqeTCk6KUM8AAEB34m4kpANtcthxWAs2LVC1o7rdba88/kqNyh8VgV4BAIDuyNVIuQsQ0r7affq49GNtrdzaofY/n/hz5ablhrlXAACgO/MEnSc9cgvSEdIRFw7UH9CSkiUdDudmk1m//c5vu9WUitE07+J5xz7oKLPHz9bs8bMj0BsAAKLH4zLkdQceSbekENKRoAzD0AufvqCd1Ts7fI4fn/zjgPN3AwAAHEuoenRCOhKSYRh6at1THV645rvDv+tf1h0AAKAjQoV0ayohHQnova/ea3dAH9NrjAb3HKxReaNCzvcNAADQFiFH0lPNEXtdQjpijsPt0OJdi7Vm35o2t8mx5Wj66OkqzCyMYM8AAECiCbXaqIWRdCSCJk+TXt/6ur44+EWb2xT1KdLIvJEaljOMh0IBAEDYuYPO7EJNOhLAgfoDeuHTF1TXVNfmNreeeqsyUzIj2CsAAJDoPI7gI+lWyl3QnZUeLtUzG59p8/HDcobpmrHXRLBHAAAAPqFG0s0pkXtdQjq61MbyjXp92+ttPv7E3ifq+yO/H7kOAQAAHCFkTXoKI+nohh5b+ZgOOw+3+fhT+p6i84ecL5MpcvVfAAAARwq+2ig16ehGDMPQst3LtLRkabvaXT3mao3IGxGZTgEAAAThamQkHd2YYRj68uCXeuXLV9rVbkLhBJ096GzmPAcAAF0i5Eg6UzAintU567Twy4XafXh3m9ukWlP1y0m/JJwDAIAu5Qq1mBHlLohXdc46PbLykXa3++mEnxLQAQBAlzIMQ54gD46azCaZkwnpiENew9uhgH735LtlNkWuxgsAAKAt3I2GDG/gkG5NN8uI4GQWhHREhGEY+vOqP7erTa4tVz+b+DNmbwEAADEh1BzpSelmNUXwtQnpCBvDMHTIfkib92/Wx6Uft7ldRnKGrh17rfpk9Ilg7wAAANonVEi32kyEdMQ2wzC0vny9lpcuV42jps3t8tPydfnoywnnAAAgJnkcIaZfTI1saS4hHZ1iGIY+2PmBVu5d2a5215xwjYblDotQrwAAADov1GqjVltky3MJ6egwt9etBz56oN3tJh03iYAOAABinsfRNXOkS4R0dFBHp1Y8Pv94XTjswgj0CAAAILxCj6RT7oIYs2TXEi3bvazd7cb0GqPLR10egR4BAACEHyPpiAsdKW+xmCw6ofcJOnvQ2cpKyYpQzwAAAMLvmCPpEZzehZCOY6px1Ghb5Ta9+9W7bW5js9p06ahLNTB7oJItyRHsHQAAQGQccySdkI5os7vsWrxrsbZUblF9U3272p7a91SdP/T8CPUMAAAgOo45kl4budcmpKMFwzC0bPcyLS1Z2qH2Z/Q/Q+cOPje8nQIAAOgCoUbSrdSkI1o8Xo/u/+j+Drf/8ck/ZmEiAADQbYQaSbcwuwuiYcehHfrnZ//sUNvR+aN12ajLZDVzOwEAgO7DzUg6ukpJTYne++o97a/f3+62yZZkzRg9g4WJAABAt+P1GPI2BR5Jt6SYZTIT0hEBLo9LT298ukPhXJL6ZfXT90Z+T3lpeWHuGQAAQNfzOEKUukR4FF0ipCekioYKPbn2yXa3G5ozVGN7j1X/Hv2VnZod/o4BAADECLc9RKlLhOvRJUJ6wvnv1//V8tLl7W535xl3KsWaEoEeAQAAxB5G0hEVXsOr3y37Xbvbfaf/d3TWwLNkMVsi0CsAAIDYFPKhUUbSEQ5Nnib9/uPft6tNWlKabjnlFiVZkiLUKwAAgNgVciQ9hZF0dJLL49KLn73YrjZXjblKI3JHyGSK/A0IAAAQi0LVpFPugk4xDENPrXtKVfaqNh1/8fCLNb5wfIR7BQAAEPs8IRYySkqLfLlL5F8hRjidTt1xxx0qLCyUzWbTpEmTtGjRoja1/e9//6spU6YoLy9P2dnZmjhxop5//vkI97hzSg+X6r5l97U5oP9y0i8J6AAAAN8IVZMejZH0hAnp119/vR599FFdc801+vOf/yyLxaJp06Zp+fLQM5288cYbOu+889TU1KR7771XDz74oGw2m2bOnKnHHnssSr1vO8Mw9MoXr+iZjc+0uc2dZ9ypnraeEewVAABAfAk9uwsPjobFmjVr9NJLL2nOnDm67bbbJEkzZ87UmDFjdPvtt2vFihVB2z7++OMqKCjQ4sWLlZLim4Lwpptu0siRI7VgwQLdcsstUfke2qLR1ag/fPKHNh9/er/TNXXI1Aj2CAAAID65Q5S7WG2MpIfFwoULZbFYNHv2bP+21NRUzZo1SytXrtSePXuCtq2trVXPnj39AV2SrFar8vLyZLPZItrv9nC4HfrL6r+0+fhLR15KQAcAAAiiqxczSoiQvnHjRg0fPlxZWVkttk+cOFGStGnTpqBtzzrrLH3xxRe666679NVXX2nnzp26//77tW7dOt1+++2R7HabVdmr9NDyh+RwO9p0/JkDztSJfU6McK8AAADil7sxVEhndpewKC8vV0FBQavtzdvKysqCtr3rrru0a9cuPfjgg3rggQckSWlpaXr11Vf1ve99Lyx9Ky8vb7V9y5Ytx2xrGIZe/OxF7aja0abXMsmky0dfrjG9xrS7nwAAAInE1RAipEdhdpeECOl2u71FuUqz1NRU//5gUlJSNHz4cE2fPl2XXXaZPB6P5s2bp2uvvVaLFi3SKaec0qm+zZ07V/fdd1+72x1sOKgn1j7Rrja3n367bEmxU6IDAAAQiwyvIXdj4Jp0k8XEPOnhYrPZ5HQ6W213OBz+/cH87Gc/06pVq7RhwwaZzb5/NV1xxRU6/vjj9ctf/lKrV6/uVN9uuukmXXLJJa22b9myRddee23ANv/Y8A/trd3b5tc4a+BZOnPAmTKbEqK6CQAAoFNcjV7JCBzSkzLMUVnwMSFCekFBgfbt29dqe3OZSWFhYcB2TU1Nevrpp3X77bf7A7okJSUl6cILL9Tjjz+upqYmJScnd6pvgUpxAqm2V+vPq//crvP/fOLPlZuW25GuAQAAJKRQCxlFo9RFSpCQXlRUpCVLlqi2trbFw6PNo+BFRUUB2x06dEhut1sej6fVPpfLJa/XG3BfuDV5mvSfrf/Rlwe/bFe7G066gYAOAADQTqFndon8KLqUILO7TJ8+3V9L3szpdGr+/PmaNGmS+vXrJ0kqLS3V1q1b/cf06tVL2dnZeu2119TU1OTfXl9frzfffFMjR46M+DSMTrdT/9jwj3YH9JvG36T+PfpHqFcAAADdV6iFjKIx/aKUICPpkyZN0owZM3TnnXeqoqJCQ4cO1bPPPquSkhI9/fTT/uNmzpypZcuWyfimBslisei2227T//7v/+qUU07RzJkz5fF49PTTT2vv3r164YUXIt73BZsWqGB428phJOn8IefrlL6nRKVWCgAAoDsKNZIejYdGpQQJ6ZL03HPP6a677tLzzz+v6upqjR07Vm+99ZbOPPPMkO1++9vfatCgQfrzn/+s++67T06nU2PHjtXChQt1+eWXR6n3bXPt2Gs1NGdoV3cDAAAgroVebZSR9LBKTU3VnDlzNGfOnKDHLF26NOD2H/zgB/rBD34QoZ51Xq4tVz+d8FNZzJau7goAAEDc8zi6viY9YUJ6d5SZnKnLRl2mQT0HdXVXAAAAuo1QI+mWVEbSEcJp/U7TeUPO6+puAAAAdDvM7oIOObnwZAI6AABAhMTC7C6E9DgzMm+kvjv8u13dDQAAgG6L2V3QLmf0P0PnDj63q7sBAADQrTG7C9rshpNuYHEiAACAKAhd7sJIOiSZTCb9aNyPdFzWcV3dFQAAgG7PMIyg5S4ms0nmZB4chaSzBpxFQAcAAIgSr8uQ4Qk8km5JNUVtVXdCeowbnje8q7sAAACQMNyNXV+PLhHSAQAAAL+Qc6SnRWcUXSKkAwAAAH7uxlALGTGSDgAAAERdyOkX0wjpAAAAQNSFHkmn3AUAAACIOkbSAQAAgBhDTToAAAAQY5jdBQAAAIgxIedJp9wFAAAAiD5XQ6iRdEI6AAAAEHXuECE9KZ2QDgAAAESVYRjBHxw1mZiCEQAAAIg2j9OQ1x24Jt1qM8lkJqQDAAAAURWy1CUjurGZkA4AAAAo9EOj0axHlwjpAAAAgKRjzOxCSAcAAACiz90QfI50RtIBAACALhB0ZhcR0gEAAIAuESsLGUmEdAAAAEDSsWrSozf9okRIBwAAACTFzmqjEiEdAAAAkMQUjAAAAEDMcTcGn92FKRgBAACAKDMMI+jsLiaLSZYUatIBAACAqHLbDRnewCPp1jSzTCZCOgAAABBVsfTQqERIBwAAAI7x0Gh0R9ElQjoAAAAQcrXRaC9kJBHSAQAAgGMsZERIBwAAAKKOmnQAAAAgxrgags+RTkgHAAAAukDImnRCOgAAABB9IWvS05jdBQAAAIg6atIBAACAGBN6nnRCOgAAABBVhteQ2x74wVGz1SRzMuUuAAAAQFS57YZkBA7p1jSzTCZCOgAAABBVsbaQkURIBwAAQIKLtYdGJUI6AAAAElysPTQqEdIBAACQ4EIvZBT9enSJkA4AAIAEF3ohI0bSAQAAgKijJh0AAACIMa6GwNMvSoR0AAAAoEswBSMAAAAQY0I+OEpNOgAAABB9oWvSmd0FAAAAiCqvx5DbHjikm5NMsiQzkg4AAABEVahSl656aFQipAMAACCBxeJDoxIhHQAAAAnM3Rh8+sWuemhUIqQDAAAggcXiQkYSIR0AAAAJLFS5CyEdAAAA6AIha9LTumb6RYmQDgAAgAQWciEjRtIBAACA6KMmHQAAAIgx1KQDAAAAMcbVEGIKRkI6AAAAEH2UuwAAAAAxxOP0yuMMHNItKWaZrczuAgAAAERVU23wUfTkHl0bkwnpAAAASEihHhpNziSkAwAAAFHnqo/NOdIlQjoAAAASVKiFjLryoVGJkA4AAIAEFarchZF0AAAAoAvE6vSLEiEdAAAACSpUTTohHQAAAOgCTXUhQnoWIR0AAACIupDzpGdaotiT1gjpAAAASDiG1wha7mIym2RN67rVRiVCOgAAABKQq8ErGUbAfUmZZplMhPSocDqduuOOO1RYWCibzaZJkyZp0aJFbW7/73//W6eeeqrS09OVnZ2t0047TYsXL45gjwEAABApIR8azej6iNz1PYiS66+/Xo8++qiuueYa/fnPf5bFYtG0adO0fPnyY7a99957dfXVV6tfv3569NFH9cADD2js2LHat29fFHoOAACAcAs1R3pyZtdHZGtXdyAa1qxZo5deeklz5szRbbfdJkmaOXOmxowZo9tvv10rVqwI2nbVqlX63e9+p0ceeUS33HJLtLoMAACACAo1R7o1retDetf3IAoWLlwoi8Wi2bNn+7elpqZq1qxZWrlypfbs2RO07Z/+9Cf16dNHv/zlL2UYhurr66PRZQAAAESQqyFwPbrU9auNSgkS0jdu3Kjhw4crKyurxfaJEydKkjZt2hS07YcffqgJEyboL3/5i/Lz85WZmamCggI9/vjjYelbeXm5NmzY0OrXli1bwnJ+AAAAtBbLq41KCVLuUl5eroKCglbbm7eVlZUFbFddXa3Kykp98sknWrx4se655x71799f8+fP189//nMlJSXppptu6lTf5s6dq/vuu69T5wAAAED7hKpJJ6RHid1uV0pKSqvtqamp/v2BNJe2HDp0SC+99JKuvPJKSdL06dN1wgkn6IEHHuh0SL/pppt0ySWXtNq+ZcsWXXvttZ06NwAAAAILFdJjodwlIUK6zWaT0+lstd3hcPj3B2snSUlJSZo+fbp/u9ls1pVXXql77rlHpaWl6t+/f4f7VlBQEHCUHwAAAJHTdNgTdF8sjKR3fQ+ioKCgQOXl5a22N28rLCwM2C4nJ0epqanKzc2VxdJyadhevXpJ8pXEAAAAIL646oKPpKdkd31E7voeREFRUZG2b9+u2traFttXr17t3x+I2WxWUVGRDh48qKamphb7muvY8/Pzw99hAAAARIzXY8htDxzSzUkmWVK6PiJ3fQ+iYPr06fJ4PJo3b55/m9Pp1Pz58zVp0iT169dPklRaWqqtW7e2aHvllVfK4/Ho2Wef9W9zOBz65z//qdGjRwcdhQcAAEBscjfG9kOjUoLUpE+aNEkzZszQnXfeqYqKCg0dOlTPPvusSkpK9PTTT/uPmzlzppYtWybD+HbezJtuukn/+Mc/dPPNN2v79u3q37+/nn/+ee3evVtvvvlmV3w7AAAA6AR3Y4g50mNgISMpQiG9oaFBDodDOTk5MplMkXiJdnvuued011136fnnn1d1dbXGjh2rt956S2eeeWbIdjabTYsXL9btt9+uZ555Rg0NDSoqKtLbb7+t888/P0q9BwAAQLg01YV4aDSjm4T0rVu3atGiRfrkk0+0Zs0a7d+/3z+TislkUnZ2tkaMGKHTTz9dZ5xxhs477zz/1IfRlJqaqjlz5mjOnDlBj1m6dGnA7b169dKCBQsi0zEAAABEVaiHRpMy4zik19fX64UXXtD8+fO1bt06SWpRItLMMAxVVVVp5cqVWrVqlR555BFlZWXpyiuv1A033OBf8RMAAACIFrc9eLlLXNakO51OPf7443r44Yd16NAhGYah3NxcTZo0SSeffLJOPPFE5eXlqWfPnkpNTVV1dbWqq6tVUlKitWvXau3atfrss880b948/f3vf9cFF1ygBx98MOjsKgAAAEC4ueqDl7vEZU360KFDVVZWph49emjWrFm68sorNWXKFJnNx/5mZs+eLUnat2+fXn75Zf373//Wu+++q/fff19/+9vfdOONN3bsOwAAAADaoak29std2tULt9uthx56SKWlpZo3b57OOeecNgX0Ix133HG65ZZbtGrVKn3yySeaNm2a9u/f365zAAAAAB3lbgge0pNjJKS3ayS9pKREKSkpYXvxU089VW+88UarhYIAAACASImHedLb1YtwBvQjJScnR+S8AAAAwNFCPThqscVhSAcAAADindseZCTdZJI1NTbW+CGkAwAAIGF4XIa8rsAj6dZUk0zm2AjpEVlx9GgOh0Nz587Ve++9p9LSUklS//79NW3aNN14442y2WzR6AYAAAASnKs+eD26NUbq0aUojKSXlZWpqKhIt9xyi7788kv16dNHffr00datW/XLX/5S48aNU1lZWaS7AQAAAIScIz1WHhqVohDSf/GLX6iiokJvvfWWdu/erQ8//FAffvihdu3apffff18HDx7UL3/5y0h3AwAAAJC7IcRqoxmxE9LDXu7SXM7S7L333tPNN9+sMWPGtNo3YsQIzZo1S08++aT27Nkjw/j2ovXv3z/cXQMAtMWaNdKqVdLnn0v790tVVZLXK+XnSwUF0jnnSIMHd3UvEa/Ky6XFi6W1a6Xt2333V3Ky1L+/NGiQdO65UkZGV/cS3Vio6RdjZbVRKQIhfeDAgTKZWhbc//GPf9Qf//jHY7Y7kscT/L8iAAARtGCBL6g3y8iQ7HZp3z5pxw7p44+lGTOkiy/usi4iTpWVSd/7nnTEoJwyMqTGRmnrVmnDBunDD6XbbpMGDOi6fqJbc4VYyCgpPTYeGpUiENKfeeYZf0g3DEM33XSTpk+frvPOOy/g8R988IEWLlyoefPmhbsrAICOOO0032hmUZHUt69vlNPrlXbulO67zxekXn5ZGjZMGjmyq3uLeNI8AHfmmb5/5E2Y4AvpTU3SsmXS3XdLNTXSo49Kf/iDFKH1WZDYEnYk/frrr2/x56eeekqfffaZ/vKXvyg3N7fFvqqqKj366KMaN26crrvuunB3BQDQEdde23qb2ewL5b/6lW+Us6LCN6JOSEd7ZGdL//qXNHRoy+3JydLUqVJlpXT//b4SmNWrfWEeCDN3Y/Ca9G4d0o/26KOPaurUqRo0aJAuv/xyDRkyRJK0c+dOFRcXy+12a9GiRZHuBgAgHKxWX+1wRYVvxDOQw4eld96RPv3UF7ocDiktLfCxc+dK48dHrLuIMZmZvl/BjBol5eX57puSksAhnfsLnRS63CWBQvrpp5+ulStX6u6779bLL78su90uSbLZbJo6dap+97vfaezYsZHuBgAgHFwuX3iSfA+SHm3/funBB78N8FarxFoYaI+MDF/49gYIUtxfCINQIT2W5kmPymJGJ554ol5//XV5vV4dPHhQkpSfny+zOXYuBAAghNpa6auvpDlzfAHKZJLOPrv1cfPm+QJUaqp0ww2+muOUFCkrS3rgAd/op8Ui/frXvu1HTRqABFdf73tAWfI9D3E07i+EgZuR9NbMZrN69+4dzZcEAHTUmjXST3/acltjo69c4cYbfWUvR/rqK9/sL5I0c6Z06qnf7hs8WPrTn6Tvftd3jqYmafr0iHYfcej1133/W5OaKk2c2HIf9xfCwDCM4CPpJpOstm46u0tZWZlee+01bdy4scWI+UknnaRLL71UhYWF4Xw5AEAkJSVJOTm+r2tqfOUH6enSVVdJgcoUP/3U93turm+GmKNlZUmXXSa98IL03nvSD34Qsa4jDq1dK737ru/r733Pd78cifsLYeBxGjI8gR8ctaaZZDJ3s5De2NioW2+9Vf/4xz/k9XpbLEpkMpk0f/583XLLLbrxxhv1xz/+UWnBHvAAAMSOk06SPvjA97XLJW3Z4puC8e9/lz76SLrlFl9ob1ZW5vt92DBfyUEg48b5QtT27b7p+IIdh8RSWir95je+fwiOHStddFHrY7i/EAbxUuoihSGkO51OnXPOOVqzZo0Mw9DAgQN15pln+kfNy8rK9PHHH2vXrl2aO3euNm3apKVLlyo5ObnTnQcARElSki88/eY3vqC+bZu0cKF05PS5DQ2+34+abreF5v9Rdbt9de49e0auz4gPFRXSz34mVVf7ylZ+/nPfMw9H4/5CGDTVhwjpGd0spD/00ENavXq1MjMzNXfuXF111VUBj3v55Zf1ox/9SKtXr9ZDDz2ku+++u7MvDQCINovF98DoV19Jy5e3DOnNkwFYQ/zV4na3PBcSW1WV77mHsrJvA3pqauBjub8QBvE0kt7p3rz44osymUwhA7okXXHFFZo7d64Mw9A///nPzr4sAKCrNI9OOhy+Oaub9ejh+/3QoeBt9+/3/Z6U5JtqD4mrrs43gl5SIh13nPTkk6HnUOf+QhjEy/SLUhhCemlpqVJSUnTFFVcc89gZM2YoNTVVpaWlnX1ZAEBX+WZiAEktRz0HDfL9vnWrZARZ0W/NGt/vw4d/OzKKxGO3S7/8pa92vFcv6amnfIsYhcL9hTBwxVG5S6d7k5ubq+Tk5DbNeW6xWJSUlKTcUPVkAICu4/GE3u9ySf/9r+/rAQN8c1E3mzDBF4wqK6VPPmnd9tAh6Y03fF+fe254+ov409Qk3Xqrb7aWnBxfQG/L7G/cXwiDhCp3Oe+881RXV6e1a9ce89i1a9eqrq5O559/fmdfFgAQCZs2ST/5iS+IV1d/u93lktat86322Py/od//fsu22dnS1Km+r+fPlz7++Nsa4e3bffXGDofUp490+eUR/kYQk7xe6be/9Y14Z2VJTzzh+8deW3B/IQxcDUH+F0axF9I7/eDovffeqzfffFPXXXed3nvvPfU/enGLb+zZs0fXXXed8vLydO+993b2ZQEAkbJ2re+XJKWl+ep76+t9o+yNjb4H9666Sjr55NZtr7zSVxe8ebNvdchnnvG1b56to0cP6Y9/9J0XiWfTJmnJEt/XTqd0880t99fVfVvKcsop0g9/2HI/9xc6qak2+P8Wxlq5S6dDeklJif7v//5Pv/71rzV69GhdeeWVmjx5cospGD/66CP9+9//VlJSkubMmaNdu3Zp165drc515plndrY7AIDOGDVKuvde30jn1q2+EoL6el/o6dtXKiiQJk/2jVYGkpQk/c//+OZR/+gjae9e32jn4MHSd74jXX/9twskIfEcWUvudPp+HamxMfDXzbi/0ElNtcHLXZJ7xNaMQCbDCPb0RduYzWaZvvkXrGEY/q+PFmqf5Fv0yH3k1EkJbsOGDRo/frzWr1+vcePGdXV3AMBn06Zj160fzWKRiooi0Rt0N9xfiCCvx9C6Bw4GfPDYZDHp5P/ND5lVA4lkXuv0SHr//v3b/Q0BAAAA0eSq9wadGSg50xxzeTYs5S4AAABALAtZj54VW6UuUhhmdwEAAABiXch69KzYi8Sx1yMAAAAgzFyEdAAAACC2NNUFL3dJptwFAAAAiD7KXQAAAIAYEzqkM5KOdjp8+LD/6+rqatXX10uSPB6PKisr5XK5JEl2u12HDh3yH1tTU6O6ujpJktfrVWVlpZqamvzHVlZWtniN5mMNw1BlZaWc3yww4XA4VFlZqebp9Gtra1VbW9viWIfDIUlyOp0tjq2rq2vR/8rKStntdklSU1OTKisr5fV6/cfW1NT4jz106JD/WJfLpcrKSnm+mTu3vr5e1UcsV15VVaXGbxa9aD62ec79hoaGFsdWV1eroaFBkuR2u1tcw8bGxhbXMND1PvIaHut6N1/DQNf76GvYfOzR17C2trbVNTz6eh95DY889shr2Hy9j7yGR1/vo69h87FHX8OqqqpW1/DI611VVRXweh99zx59vWtqatp8vY+8Z4++3s33bLiud/M1DHS9g92zga73kdcw0PXu6D179PWOymdEfb2c35zX4XKpsr7+22tot6v2m+vQfKzjm2P5jOAz4pifEXV1cn1z3samJh365vpKUk1jo+q/+d48Xq8q6+vV9M15+YyIsc+IGM0Rzhq3v63d7vAfa7fb5bY6O3S9j9wXboT0GPfJJ5/4v/7www+1efNmSb4f6OLiYh08eFCStH37dr311lv+Y5ctW6YNGzZI8v2AFBcXa//+/ZKkr7/+Wq+//rr/2OXLl2vNmjWSfDdncXGx9u3bJ0navXu3iouL/T8wK1eu1MqVKyX5friKi4u1e/duSdK+fftUXFzs/4Ffs2aNli9f7n+d119/XV9//bUkaf/+/SouLvb/YG7YsEHLli3zH/vWW29p+/btkqSDBw+quLjY/8GyefNmffjhh/5j3333XW3ZskWS70OluLjY/0P9xRdf6P333/cf+8EHH+izzz6T5PuBLi4u9n+Abd26Ve+++67/2CVLlmjTpk2SfB90R17vHTt26I033mhxvdevXy/J9yFz5PUuKSnRa6+95j/2k08+0erVqyX5fsCLi4u1Z88eSVJpaamKi4v9HzqrVq3SihUr/G2Li4v9056WlZW1uN5r167VRx991OJ679y5U5J04MCBFtd748aNWtK8NLekt99+W9u2bZPk+2ArLi72f+B++umn+u9//+s/9v3339eXX37Z4no3fzB++eWXLa73f//7X3366actrnfzXwLbtm3T22+/3eJ6b9y4UdK39+yBAwckSTt37mxxz3700Uda+82y9c33bFlZmf96FxcX+49dsWKFVq1aJcn3F01xcbFKS0slSXv27FFxcbH/w3j16tUtfuZee+01//Vuvmeb/+JZv359i3v2jTfe0I4dOyR9e882/yW7adOmFtf73Xff1datWyX5As2R1/uzzz7TBx980OJ6f/HFF5J8fxEeeb23bNnS4p6NymeEx6PijRu175s+7K6qUvHGjWqeeXjl119r5Tc/54ak4o0btfub95zPCD4jjvkZsX69P5hvO3BAb3/++bfXe9s2bfzmOjhcLhVv3KgD3/SBz4gY+oyI0Rzx6qvFqq+0f9P2QIuV77/66it9tcd3f7f3M+LI9y3cOr3i6NKlS9XU1KTJkycrJSWl1f5du3Zpx44d6tevn0aNGuW/sc4777zOvGy317yC1eLFizVlyhRJvn+9JSUlKSMjQx6PR9XV1erRo4eSkpJkt9vV2Nio3NxcSb6bzGKxKDMzU16vV1VVVcrKylJycrLsdrsaGhqUl5cnyfdDbTablZmZKcMwdOjQIWVmZiolJUUOh0P19fXKzc2VyWTy/+s3KyvLf2xGRoZSU1PldDpVV1fnP7aurk5er1c9evSQ5PsXcHp6umw2m5qamlRbW6ucnByZzWbV1dXJ4/EoOztbku8vgbS0NNlsNrlcLh0+fFg9e/aUxWJRfX29XC6XevbsKcn3AZaamqq0tDT/sdnZ2bJarWpoaFBTU5P/2OrqaiUnJys9PV1ut1s1NTX+a9jY2Ci73e6/hoGu95HX8FjXu/kaBrreJpOpxTVsPvboa1hbWyvDMFpcw6Ov95HX8MjrfeQ1bL7eR15Dt9vd4nrbbLYW17D52KOvYVVVlVJSUlpcwyOvt9PpVM43y3Ifeb2PvmePvt41NTWyWq1tut5H3rNHX+/mezZc17v5ng10vYPds4Gu95H3bKDr3dF71uFwtLjeEf+M2LhRh2prlZmSopSkJDlcLtU7ncpNT/ddw29GrrJsNt/1bmhQRlqaUidM4DOCz4hjf0asWKEeKSlKsljU2NQke1OTcjMyfNewsVFWs1kZqanyeL2qbmxUVmqqklNSZB8xgs+IWPmMiNEcUbmvWruf8cpsNqmpqUkej1c2W6okyWVx6qRb8zp0vZctW6Zzzz03IiuOdjqkn3vuuVqyZIkWLlyoSy+9tNX+6dOn67XXXtMTTzyhH//4xzKbzTKbzf5/lSKwSC4zCwAdxrLtiCTuL0RIQ5lLX8yrCrgv/bgkHf+jnA6dN5J5rdPlLhdeeKEMw2jx39bN3G63/7/ALrroIv/2Tv67AAAAAGizeHtoVApDSG8O30fWPDVbvny5amtrNXr0aPXr16+zLwUAAAC0W1NtqDnSY/MRzU73auTIkRo0aJD279/vf8Cg2dtvvy2TydRiFB0AAACIpnibI10K0+wuF154oSS1Knl55513JImQDgAAgC4TeiS9m5a7SNJ3v/tdGYbRYuqe3bt3a8uWLerRo4dOP/30cLwMAAAA0G6hRtKTuvNI+pQpU2Sz2bR+/Xr/fJvNo+jnnXeezObY/OYBAADQ/YUsd8mMzZwall6lpKRoypQpMgzD/wDpO++8Qz06AAAAupRhGGqqCVLuYjJ173IXyVd33jwVo9Pp1JIlS2Qymfz16gAAAEC0uRsNGd7A038nZZhltpqi3KO2CVtInzZtmiTfcqmLFi1SY2OjTj75ZP9qVAAAAEC0NdWFeGg0RktdpDCG9AEDBmj06NGqra3V//7v/1LqAgAAgC7nqgvx0GgihHTJN5puGIY+/fRTSUy9CAAAgK4VKqQnxEi69G0oN5lM6t27t8aNGxfO0wMAAADt4moIMZKekSAh/fTTT1dWVpYMw+CBUQAAAHQ5tz14SLemxW5It4b1ZFar/vKXv2jXrl26+OKLAx5z9913y2SKzadoAQAA0L24GwPP7CJJVluChHRJmjlzZsj99957b7hfEgAAAAioqTb47C4JU+4CAAAAxBJnTYgHR7NiNwrHbs8AAACATjC8hpoOh1httEdsrjYqEdIBAADQTbkavDI8gWvSkzNjd7VRiZAOAACAbirkHOkxXOoiEdIBAADQTbnqQ602GrulLhIhHQAAAN1UyJAewzO7SIR0AAAAdFPOmvicflFqZ0g/dOhQRDpRVVUVkfMCAAAgcTmqgof0lJ7dqNxl4MCBuvXWW1VWVhaWFy8uLtaECRP0+OOPh+V8AAAAQLOg0y9KSunZjUbSBw8erMcee0yDBw/W97//ff3rX/9SQ0NDu15w06ZNuvPOOzV48GDNmDFDn3/+uUaOHNmucwAAAADH0lQbvCY9JYbnSJcka3sO3rx5s1588UXdfffdeuONN/Tmm28qKSlJJ5xwgsaPH6+xY8cqLy9PPXv2VEpKimpqalRdXa2SkhKtW7dO69ev14EDB2QYhqxWq2644Qbdc8896tu3b6S+PwAAACQgw2sEnYLRZDbFfE16u0K6JP3gBz/QlVdeqTfeeEPPPPOM3n//fa1fv17r16+XyRR8QnjD8E0kP2jQIF1//fW6/vrr1a9fv473HAAAAAjC1eCV4Q28kFFShlkmc+wuZCR1IKRLksVi0aWXXqpLL71UBw8e1JIlS7RixQqtXbtW+/fvV2VlpZxOp3JycpSXl6cRI0bo9NNP1xlnnKGTTz453N8DAAAA0EKoUpdYX8hI6mBIP1J+fr6uuOIKXXHFFeHoDwAAANBpoR4aTcqK7Xp0qZ0Pjp599tmaMWNGi22lpaXat29fWDsFAAAAdEboh0a72Uj60qVL1adPnxbbBg4cqIKCAoI6AAAAYkZTbYiR9MxuNpJutVrV1NTUanvzQ6EAAABALGg6HKImPQ5G0tvVwz59+qi6ulpff/11pPoDAAAAdFqokfRYnyNdame5y3nnnadnnnlGkyZN0pQpU5SRkSFJOnz4sG644YY2n8dkMunpp59uX08BAACANkqo2V0efPBBrVixQlu3btXChQv92+12uxYsWNDm8xDSAQAAECnxvpCR1M6Q3rt3b3322Wd6//339cUXX8hut+vee+9VRkaGbr311kj1EQAAAGgzt90IupCRNT32FzKSOjBPusVi0bRp0zRt2jRJ8of0e+65J+ydAwAAANorVD16cmbsj6JL7XxwNJB77rknLkbRnU6n7rjjDhUWFspms2nSpElatGhRu88zdepUmUwm/exnP4tALwEAANBZoerRk+IkpHd6xdF4GUG//vrrtXDhQv3qV7/SsGHDtGDBAk2bNk1LlizRGWec0aZzFBcXa+XKlRHuKQAAADrDWRNiJD0OZnaRwjCSHg/WrFmjl156Sf/3f/+nOXPmaPbs2Vq8eLEGDBig22+/vU3ncDgcuvXWW3XHHXdEuLcAAADojKYQIT0lm5AeMxYuXCiLxaLZs2f7t6WmpmrWrFlauXKl9uzZc8xz/OEPf5DX69Vtt90Wya4CAACgk5qCzOwiSSnZ8RF/46OXnbRx40YNHz5cWVlZLbZPnDhRkrRp06aQ7UtLS/XQQw/p4Ycfls1mi1Q3AQAAEAbUpMeJ8vJyFRQUtNrevK2srCxk+1tvvVUnnXSSrrrqqoj0rby8vNX2LVu2hP21AAAAEoGrLkRNelZ8lLskREi32+1KSUlptT01NdW/P5glS5bo1Vdf1erVqyPSt7lz5+q+++6LyLkBAAASjWEEX8hIUlwsZCQlSEi32WxyOp2ttjscDv/+QNxut37xi1/ohz/8oSZMmBCRvt1000265JJLWm3fsmWLrr322oi8JgAAQHflbjTkdQdfyMhsif2FjKQECekFBQXat29fq+3NZSaFhYUB2z333HPatm2b5s6dq5KSkhb76urqVFJSol69eiktLa1TfQtUigMAAID2c1S5g+6Ll1IXKUEeHC0qKtL27dtVW1vbYntzCUtRUVHAdqWlpXK5XDr99NM1aNAg/y/JF+AHDRqkDz74IKJ9BwAAQNs1HQ5e6pKaGz8hPSFG0qdPn64//vGPmjdvnn8KRafTqfnz52vSpEnq16+fJF8ob2xs1MiRIyVJV111VcAAf+mll2ratGn60Y9+pEmTJkXt+wAAAEBoTYdDzJEeJwsZSQkS0idNmqQZM2bozjvvVEVFhYYOHapnn31WJSUlevrpp/3HzZw5U8uWLZNh+OqYRo4c6Q/sRxs0aJC+//3vR6P7AAAAaCNniJH05DiZI11KkJAu+cpT7rrrLj3//POqrq7W2LFj9dZbb+nMM8/s6q4BAAAgTEKtNhpPNekJE9JTU1M1Z84czZkzJ+gxS5cubdO5mkfaAQAAEFsch0KUu2THT0iPnzF/AAAAIASvx5CjKkhIN5mU0jN+om/89BQAAAAIwVXnlYJUPKT0MMuSHD/RN356CgAAAITgDDGzS3IclbpIhHQAAAB0E6EeGo2n6RclQjoAAAC6iVALGSX3iK/YG1+9BQAAAIIIVe4STzO7SIR0AAAAdBMh50hnJB0AAACIPmdN8HIXRtIBAACAKDMMQ86QI+mEdAAAACCqXPVeGZ7Ac6QnZVpktpqi3KPOIaQDAAAg7oUaRU/Jjr/IG389BgAAAI7S1I3q0SVCOgAAALqBkPXohHQAAAAg+ppCzZEeZw+NSoR0AAAAdANNdd1ntVGJkA4AAIBuwBUipCdlxF/kjb8eAwAAAEdxhih3ScqMv8gbfz0GAAAAjuBxeuVuCDySbk4yKSk9/iJv/PUYAAAAOIKjOsRDoz0tMpniayEjiZAOAACAOOc8RkiPR4R0AAAAxDVnVfCQnppDSAcAAACiznGIkXQAAAAgpoQK6am5hHQAAAAg6hyhyl1yrVHsSfgQ0gEAABC3PE6vXHWBQ7rZaorL1UYlQjoAAADiWKhR9JSc+Jx+USKkAwAAII6FrkePz1IXiZAOAACAONYdHxqVCOkAAACIY45D7qD7COkAAABAFwi1kFFKnC5kJBHSAQAAEMcodwEAAABiiKvRK7fdG3CfJcWspPT4jbrx23MAAAAktGPVo8fr9IsSIR0AAABxqruWukiEdAAAAMSp7vrQqERIBwAAQJzqrgsZSYR0AAAAxCl7RYiadEbSAQAAgOjyuIzQI+l5hHQAAAAgquwH3DK8RsB9KT0tsqbGd8yN794DAAAgITWUu4LuSy9MimJPIoOQDgAAgLjTuD94PXpaQXw/NCoR0gEAABCH7AdChPQ+hHQAAAAgqrweQ42hQnpvQjoAAAAQVY5DHnldgR8ataaZlZQR/xE3/r8DAAAAJJRQ86OnH5ckk8kUxd5EBiEdAAAAcSXU/OhpveK/1EUipAMAACDOOA6FWGk0zhcxakZIBwAAQFyxVwQfSU/JIaQDAAAAUWV4DdkPhpjZhXIXAAAAILqaar0yPIFndknKtMhq6x7xtnt8FwAAAEgIDeWuoPtSenaPUheJkA4AAIA40t1XGm1GSAcAAEDccFQFf2g0vYCQDgAAAERdY3nwkXTKXQAAAIAoM7yGnNUhFjKi3AUAAACIrqZar7zuEDO7pHafaNt9vhMAAAB0a437Q6w02k0WMWpGSAcAAEBcCLWIka139yl1kQjpAAAAiBOOyhD16N1kpdFmhHQAAADEBXtliHKXfMpdAAAAgKgyDCPkSLotj5F0AAAAIKpcdV55nN6A+6w2s6xppij3KLII6QAAAIh5DWUhSl1yLTKZCOkAAABAVNkrgof0tMKkKPYkOgjpAAAAiHmNB4KH9PRutNJoM0I6AAAAYl6okXRbN5t+USKkAwAAIMZ53YYch0LM7NLNpl+UCOkAAACIcY4qjwyvEXBfSrZFlpTuF2m733cEAACAbiVUqUtqfvcrdZEI6QAAAIhx9oMh6tG7YamLREgHAABAjGvcn1gPjUqEdAAAAMS4xvIQc6R3w+kXJUI6AAAAYpjb4VVTbeCZXUwWk2zUpAMAAADRFbLUJc8is8UUxd5EDyEdAAAAMStUqYutd/ccRZcI6QAAAIhhDeWuoPvSC5Ki2JPoIqQDAAAgZoV8aLSAkXQAAAAgqjwuQ/bKwA+NSt13ZheJkA4AAIAY1Vjukgwj4L6UHKusqd03ynbf7+woTqdTd9xxhwoLC2Wz2TRp0iQtWrTomO2Ki4t15ZVXavDgwUpLS9OIESN06623qqamJvKdBgAASGC1u0LVo3ffUXQpgUL69ddfr0cffVTXXHON/vznP8tisWjatGlavnx5yHazZ8/Wli1bdO211+ovf/mLLrjgAj3++OM69dRTZbfbo9R7AACAxNMY4qHRzAHd96FRSere/wT5xpo1a/TSSy9pzpw5uu222yRJM2fO1JgxY3T77bdrxYoVQdsuXLhQZ511Vott48eP13XXXad//vOfuvHGGyPZdQAAgITVUBb8odGMvt07pCfESPrChQtlsVg0e/Zs/7bU1FTNmjVLK1eu1J49e4K2PTqgS9Kll14qSdqyZUvY+woAAACpqd4TeqXRXt17rLl7f3ff2Lhxo4YPH66srKwW2ydOnChJ2rRpk/r169fm8+3fv1+SlJeX1+m+lZeXq7y8vNV2/gEAAAASWcipF3tbZbZ2z5VGmyVESC8vL1dBQUGr7c3bysrK2nW+hx9+WBaLRdOnT+903+bOnav77ruv0+cBAADoTkKVunTn+dGbdf/vUJLdbldKSkqr7ampqf79bfXiiy/q6aef1u23365hw4Z1um833XSTLrnkklbbmx9WBQAASEShHhpNL+ze9ehSgoR0m80mp9PZarvD4fDvb4uPP/5Ys2bN0vnnn68HH3wwLH0rKCgIOMoPAACQyEKNpHf36RelBHlwtKCgIGDdd/O2wsLCY55j8+bNuuSSSzRmzBgtXLhQVmv3vzkAAAC6QqI/NColSEgvKirS9u3bVVtb22L76tWr/ftD2blzpy644AL16tVL77zzjjIyMiLVVQAAgISX6A+NSgkS0qdPny6Px6N58+b5tzmdTs2fP1+TJk3yz+xSWlqqrVu3tmi7f/9+nXfeeTKbzXr//feVn58f1b4DAAAkmkR/aFRKkJr0SZMmacaMGbrzzjtVUVGhoUOH6tlnn1VJSYmefvpp/3EzZ87UsmXLZBiGf9sFF1ygr7/+WrfffruWL1/eYoXS3r17a+rUqVH9XgAAALq7hrLEfmhUSpCQLknPPfec7rrrLj3//POqrq7W2LFj9dZbb+nMM88M2W7z5s2SpD/84Q+t9k2ePJmQDgAAEGaN+0M8NFqYGPE1Mb5L+aZbnDNnjubMmRP0mKVLl7baduSoOgAAACLL7fCqqdYbcJ/JnBgPjUoJUpMOAACA+FC32yUFGSRNybHIbOn+D41KhHQAAADEkLqSpqD7MvsnRj26REgHAABADKktCf7QaObA5Cj2pGsR0gEAABAT3A5vyIdGMwcwkg4AAABEVX1pqHp0q1J6WKLco65DSAcAAEBMqA1Vj55Ao+gSIR0AAAAxoi5EPXrWQEI6AAAAEFVuh1cN5SHq0RPooVGJkA4AAIAYELIevacloerRJUI6AAAAYkDIevQEG0WXCOkAAACIAdSjt0RIBwAAQJc6Zj36AEbSAQAAgKiq33OMevTsxKpHlwjpAAAA6GKhSl0SsR5dIqQDAACgi9XuYhGjoxHSAQAA0GWa6jxqKAv10Cgj6QAAAEBUVW52BN2Xkp2Y9egSIR0AAABdqPpLZ9B9WYMTcxRdIqQDAACgiziqQ5e65BWlRrE3sYWQDgAAgC5RszX4KHpKjlUZ/RLzoVGJkA4AAIAuUrM9eEjPGZ0ik8kUxd7EFkI6AAAAos7V6FVdafBSl56jUqLYm9hj7eoOAAAAdLWXT/6qXcenFyTpojcHRKg3iaFmm1OGJ/Aqo9Z0s9ILEzumJvZ3DwAAIGngdzNbbavc5FD9Xpeyh6coe3jLWUYSdVrAcKrZFrzUJXt4Ype6SIR0AAAATby3d6tta+49oPq9Lh13VrqOn53TBb3qvjwuQ4d3Bl9lNHdMYpe6SNSkAwAAIMpqdzbJ6wpc6mJJMStzQOLOj96MkXQAAIB2qNvTpNJ367V/ZaMa9rnUVOtVSo5FvU62afSNPZXZv2XAbChz6e1Ldit/nE1T5h3X6nxfzKvSF/OqNOGeXhp0cZZ/+9sX71ZDuUsz1g7RVy/Xatd/alVX6lLmgCSd92I/7XqzVmvvq9Dxs3MCjvQvmb1PBzfYddEbA5Re+O1Uhod3OrV1QY0qNzvkqHTLmmaWrZdV+eNtGnldtmx5kY+HlZ8GX2W0x9Bkma2JXeoiEdIBAADaZdd/6rTtuRplDUlWzvGpMieZVLurSbvfqVPZsgZN+cdxyh4WvnKN9f93UCVv1Cl/nE2Zg5LkDT4hyjFVbXFoyax98jQZyh6Wopzj0+VxGKrf59KOf9XouLPSIx7S3Q6vDoeYerHnSEpdJEI6AABAuxx3VrqGXJ7VYnRakna9Uau1v6vQpkcqddbfWo+Yd9S+xQ2a+mI/9Rjc+RKQHS8dlqfJ0Im/ytOIa7Nb7KstaVJSRuQroau3OOV1By51MSeblD2CUheJkA4AANAuuScEXqp+0CVZ2vV6rQ6ud6ip3qPkjPDMADPyuuywBHRJclZ7JEm9J9pa7csaGJ1wXLHWHnRfzqhUWZJ5ZFIipAMAALSbq9Gr8o8bVL2tSa5aj39k2F7pkWEYatjrVvLI8IT0wsnpYTmPJOWMStH+FY3a8PBBjflprvKKUmW2RK/+237QrYay4PU6OcdT6tKMkA4AANAOB9Y2atVvDvhHpQNxN3rD9nppfcIX10bM7KmDmxw6uN6upTftkzXNrNwTUlVwRpoGXpwZttH/YA5/HXzaxaQMi7KGUOrSjJAOAADQRq5Gr1b+fwfkqvVq9I9y1P+8DKUVWGVJMclkMmnVbw+o9P06GYFLrgMyjpHnO1z+EaAPSelmnfW3QlVudqjso0YdXG9XxVq7Dqxu1Nb51Zryj+NazU4TTjVbgz8wmj8uuqP6sY6QDgAA0EaVmxxqOuxR33MyNOam1tMeNuxrXcphTvIFT7c9cBpvPODuUF/85w0yah/svCaTSflFNuUX+erSHVVubXrkkErfr9PnT1bp1If6dKg/x+Kocqt2V/CR9B5DGUU/EpX5AAAAbdRU6ytxSevdepyzbk+TqgOMFCdnW2SymNSwzy2vp+Xwttdt6OD64A9ShtI8VWJdaet/GNSVNqlxf9vCf2qOVcff1FOSQq4C2lkH1gT/Pq1pZmX0TQq6PxER0gEAANoos78vSO5d3CDHETXpTXUerbv/YMCpBS1JJuWNTVVTrUdfvXzYv93rMbTpscqQD1KG0nN0iqypZpV/0qiqLd8uDuSs8fXF8Lbuy85XDwd8vfLljZIC/+MjHAyvoUMhFjDKOzFVJjOlLkei3AUAAKCNckanqvekNB1Y3aj3LitV/njfdIwH1zuUnG1W4eR0lS1raNVu9I966qOfObTpkUrtWVSv1Fyrqrc45XF4NfC7mSp5q67dfUlKM2v4D7P15d+rtGTWPuWPt0kmqepzpzIHJil3bGqrYLzz1Vqt/7+DyhqcrKyByTJZpboSl2q2O2VJNmn0j3p27MIcw6HPHCEfps0f33pKyETHSDoAAEA7nPFoH426oadSelq0f0Wjqrc41e+8DJ2zoK+SMwNHq94T03T6o32UMzpVNVudOrjertwTUnTuc31bLYrUHsfP7qkTf5krW2+rKtbadfirJg26JFOTnyyU2dp6ZHrMj3M06JIsSb75yss/apTHYWjw97M09V/9lHdi+MOy12P4R+oDyRyYHPFVTuORyTDa8/wxomXDhg0aP3681q9fr3HjxnV1dwDAZ9MmyRN82rmALBapqCgSvUF3w/3VLVV96WhR5nO0/hdkqs8paVHsUfhEMq8xkg4AAICIObgxeC26JcWsvBMDr+Ca6AjpAAAAiAi33avaEAsY9T7FJquNOBoIVwUAAAARUbHOLsMTuLLanGxSwRnpUe5R/CCkAwAAIOwMb+g54HNGp8qSxLSLwRDSAQAAEHaVnzrkrAn+IDC16KER0gEAABBWhmFo35LW88U3Sy9MUuZAVhgNhZAOAACAsKordanpcPBR9PzxNplMlLqEQkgHAABAWB290umRkjIsyh1LqcuxENIBAAAQNm67V1WfO4Pu73NaGg+MtgEhHQAAAGFzcINdHqc38E6TSbljU6LboThFSAcAAEBYGF4j5Aqj2cOSlZxhiWKP4hchHQAAAGGxd3GDHJXuoPvzx9ui2Jv4RkgHAABAp9V85VT58uDTLqbmWZU9PDmKPYpvhHQAAAB0SlOtR1+/VhvymF4nM+1iexDSAQAA0GGG19BXrxyWuyHIw6KSUnpaKHVpJ0I6AAAAOqz8k0bV73GFPGbgdzOZdrGdCOkAAADokLrSJu39sD7kMfnjbOoxhGkX24uQDgAAgHZz1nj01Suh69Az+iZp4MWZUepR92Lt6g4AAAAgvjTVevTZE4fkdRlBj0nJtmj4Ndk8LNpBhHQAAAC0mddtaPOfDsnwBg/oMpk0+PIsWW0UbXQUVw4AAABtYhiGSt+rCx3QJfWeZFNmP+ZE7wxCOgAAANqk7KNGVayzhzymx5Bk9ZuaEaUedV+EdAAAABxT1ZcO7VsSeiYXk9mkITN6yGyhDr2zqEkHAABASAc32bXr9bqQx5iTTBp1Q09ZUxkDDgeuIgAAAIKqK21S6bv1khG6Dn34tdlKL0iKUq+6P0bSAQAA0IrH6dXexQ06sLrxmMcOuypbWQN4UDScCOkAAABooWaHU6Xv18tR6T7msb0mpKnnSFYUDTdCOgAAACT5Rs+/fr1O1V862nS8rZdVA6Yxk0skENIBAACguj1N+urftXLVe9p0vK2XVWN+ksOKohFCSAcAAEhghmGocrNDu/5T2+Y2PUemaOgVPQjoEURIBwAASFBuu1e7363ToU/bVt4iSXkn2TTwokyZzAT0SCKkAwAAJBiv29ChTx0qfb9eHqe3TW1svawacGGmsgYxi0s0ENIBAAAShL3SrYPr7arc7JC7sW3hXJL6Tc1Q71PSWEk0igjpAAAA3dzhnU7tXdyghn2udrc94eZc2fKJjNHGFQcAAOiG6ve6VPpener3tj+YS9/Unn83k9HzLkJIBwAA6Eaqtzm16z+1ctvbXs5ytN6npGnABZlh7BXai5AOAAAQx1wNXtXuatLOV2slw+jUudILkzTgokxlHJcUpt6howjpAAAAccJt98pZ7VFTrVf1+1yq/bqpQ3XmR0vOsqjf1AzlnpAahl4iHAjpAAAAMcIwDLkbDTXV+oJ402GPGsrcclR55Kh0t2tGlrZIzbOq79npyh6RQu15jEmYkO50OnX33Xfr+eefV3V1tcaOHasHHnhAU6dOPWbbffv26ZZbbtEHH3wgr9erKVOm6LHHHtPgwYOj0HMAANAduBq8cjd65XEacjV45ar3ylXnkavRkPOQW87DvlDudXWuZKUtegxNUeHkNGX2Y87zWJUwIf3666/XwoUL9atf/UrDhg3TggULNG3aNC1ZskRnnHFG0Hb19fWaMmWKDh8+rN/85jdKSkrSY489psmTJ2vTpk3Kzc2N4ncBAADCxWiu3za+KeVu/t1ryPDK98tj+H/3er75s0fyegwZbt/Xrnqv3A5D3iavPE2GvE2GvC7J7fDK3eCV227IWeOR4Yl8+D4WW75Vgy/PUnofas5jXUKE9DVr1uill17SnDlzdNttt0mSZs6cqTFjxuj222/XihUrgrZ98skntWPHDq1Zs0YTJkyQJF144YUaM2aMHnnkEf3+97+PyvcAAACkbe+aZa8y+QK1JJlMMhYd9O08MmgbLUN4qyDeyQcs40lSpkWDv5+prMHJMpkoaYkX5q7uQDQsXLhQFotFs2fP9m9LTU3VrFmztHLlSu3Zsydk2wkTJvgDuiSNHDlS55xzjl5++eWI9hsAALTkapCa6qWmBt/XrkaTb7T6m1ISt90rj8Mrj9PrG9FuMuR1GfK6jW9GxY2ECei2Xlad8LNcnXRrnnoMSSGgx5mEGEnfuHGjhg8frqysrBbbJ06cKEnatGmT+vXr16qd1+vVp59+qhtuuKHVvokTJ+qDDz5QXV2dMjOZRxQAgKggZwaU3MOi7GHJyhyQLFtvq2z5FkJ5nEuIkF5eXq6CgoJW25u3lZWVBWxXVVUlp9N5zLYjRozoVN/Ky8tbbd+0aZMkacuWLR0+NwCE3bZtkreds0uYze1vg8TUhvtr5wGT7NVHhE+TWTIORLhjscNsNSk526LkTLOS0s1K6WlReqFVnh4WVZlMqnJJ2vvNL0Rcc06z2+1hP3dChHS73a6UlJRW21NTU/37g7WT1KG2bTV37lzdd999Qfdfe+21nTo/AAAAImvjxo06/fTTw3rOhAjpNptNTqez1XaHw+HfH6ydpA61baubbrpJl1xySavtq1ev1k9/+lM9/fTTKioq6tRrhNPll1+uV199NabO2962bT2+LceFOibYvkDbt2zZomuvvVYvvPCCRo0adcy+RQvvd9uP4f2OzHl5vyOH97vtx/B+R+a83eH93rRpk2bNmhWR9zohQnpBQYH27dvXantzmUlhYWHAdjk5OUpJSQlYjnKstu3pW6BymmZFRUUaN25cp14jnGw2W0T605nztrdtW49vy3Ghjgm2L1SbUaNG8X6HuS3vd9vxfrf9GN7vyJyX9ztyeL/bfkxH3u+ePXses1/tlRCzuxQVFWn79u2qra1tsX316tX+/YGYzWadcMIJWrduXat9q1ev1uDBgxPuodGbb7455s7b3rZtPb4tx4U6Jti+SF3DSOD9bvsxvN+ROS/vd+Twfrf9GN7vyJyX9/sYjASwatUqQ5IxZ84c/zaHw2EMHTrUmDRpkn/b7t27jS1btrRo+9BDDxmSjLVr1/q3bd261bBYLMYdd9wRsT6vX7/ekGSsX78+Yq+B2MH7nVh4vxML73di4f1OLJF8vxOi3GXSpEmaMWOG7rzzTlVUVGjo0KF69tlnVVJSoqefftp/3MyZM7Vs2bJvFz+Q9NOf/lR///vfddFFF+m2225TUlKSHn30UfXu3Vu33nprV3w7AAAA6OYSIqRL0nPPPae77rpLzz//vKqrqzV27Fi99dZbOvPMM0O2y8zM1NKlS3XLLbfogQcekNfr1VlnnaXHHntM+fn5Ueo9AAAAEknChPTU1FTNmTNHc+bMCXrM0qVLA27v27evXnnllQj1LLCCggLdc889IR8qRffB+51YeL8TC+93YuH9TiyRfL9NhpEga+MCAAAAcSIhZncBAAAA4gkhHQAAAIgxhHQAAAAgxhDSAQAAgBhDSAcAAABiDCEdAAAAiDGEdAAAACDGENIBAACAGENIBwAAAGIMIR0AAACIMYR0AAAAIMYQ0gEAAIAYQ0gHAAAAYgwhPQSn06k77rhDhYWFstlsmjRpkhYtWtSmtvv27dMVV1yh7OxsZWVl6Xvf+56+/vrrgMc+/fTTGjVqlFJTUzVs2DD99a9/jco5t23bpltuuUWnnXaaUlNTZTKZVFJS0qbvL1HE6z3w1FNPacaMGerfv79MJpOuv/76Nn/PaCmW7gF+ZrtGNO4BfmZjX0fvA35u4099fb3uueceXXDBBcrJyZHJZNKCBQva3L6mpkazZ89Wfn6+0tPTNWXKFG3YsKH9HTEQ1FVXXWVYrVbjtttuM+bOnWuceuqphtVqNT7++OOQ7erq6oxhw4YZvXr1Mh5++GHj0UcfNfr162f07dvXqKysbHHs3/72N0OScfnllxvz5s0zfvjDHxqSjIceeiji55w/f75hNpuNMWPGGEVFRYYkY9euXR2/YN1QvN4DAwYMMHJycowLLrjAsFqtxnXXXReW65GIYuke4Ge2a0TjHuBnNvZ19D7g5zb+7Nq1y5Bk9O/f3zjrrLMMScb8+fPb1Nbj8RinnXaakZ6ebtx7773G448/bowePdrIzMw0tm/f3q5+ENKDWL16tSHJmDNnjn+b3W43hgwZYpx66qkh2z788MOGJGPNmjX+bVu2bDEsFotx5513+rc1NjYaubm5xkUXXdSi/TXXXGOkp6cbVVVVET3noUOHjNraWsMwDGPOnDl8cBwlXu8BwzCMkpISw+v1GoZhGOnp6fyF30Gxdg/wMxt90bgHDIOf2VjXmfuAn9v443A4jPLycsMwDGPt2rXtCun//ve/DUnGK6+84t9WUVFhZGdnG1dffXW7+kFID+LXv/61YbFYjMOHD7fY/vvf/96QZJSWlgZtO2HCBGPChAmttp933nnGkCFD/H9+++23DUnG22+/3eK4FStWGJKM559/PqLnPBIfHK3F6z1wNP7C77hYuweOxM9sdETjHjgaP7OxpzP3wZH4uY0/7Q3pM2bMMHr37m14PJ4W22fPnm2kpaUZDoejza9NTXoQGzdu1PDhw5WVldVi+8SJEyVJmzZtCtjO6/Xq008/1cknn9xq38SJE7Vz507V1dX5X0NSq2PHjx8vs9ns3x+Jc+LY4vUeQPjE0j2ArhGNewCxr6P3ARLPxo0bNW7cOJnNLSP2xIkT1djYqO3bt7f5XIT0IMrLy1VQUNBqe/O2srKygO2qqqrkdDrb1La8vFwWi0W9evVqcVxycrJyc3P9x0XinDi2eL0HED6xdA+ga0TjHkDs6+h9gMQTznuFkB6E3W5XSkpKq+2pqan+/cHaSWpTW7vdruTk5IDnSU1NbXFcuM+JY4vXewDhE0v3ALpGNO4BxL6O3gdIPOG8VwjpQdhsNjmdzlbbHQ6Hf3+wdpLa1NZms6mpqSngeRwOR4vjwn1OHFu83gMIn1i6B9A1onEPIPZ19D5A4gnnvUJID6KgoEDl5eWttjdvKywsDNguJydHKSkpbWpbUFAgj8ejioqKFsc1NTXp0KFD/uMicU4cW7zeAwifWLoH0DWicQ8g9nX0PkDiCee9QkgPoqioSNu3b1dtbW2L7atXr/bvD8RsNuuEE07QunXrWu1bvXq1Bg8erMzMzBbnOPrYdevWyev1+vdH4pw4tni9BxA+sXQPoGtE4x5A7OvofYDEU1RUpA0bNsjr9bbYvnr1aqWlpWn48OFtPhchPYjp06fL4/Fo3rx5/m1Op1Pz58/XpEmT1K9fP0lSaWmptm7d2qrt2rVrW3w4b9u2TYsXL9aMGTP8284++2zl5OToqaeeatH+qaeeUlpami666KKInhOhxes9gPCJtXsA0ReNewCxrzP3Abqv8vJybd26VS6Xy79t+vTpOnDggIqLi/3bKisr9corr+jiiy8OWK8eVHvmikw0M2bMMKxWq/HrX//amDt3rnHaaacZVqvVWLZsmf+YyZMnG0dfxtraWmPIkCFGr169jD/84Q/GY489ZvTr188oLCw0KioqWhz7xBNPGJKM6dOnG3//+9+NmTNnGpKMBx98MOLnrKmpMe6//37j/vvvNy644AJDknHrrbca999/v/HXv/41HJcw7sXrPfDGG2/439vk5GTjpJNO8v958+bNYb5K3Vss3QP8zHaNaNwD/MzGvo7eB/zcxqe//vWvxv3332/85Cc/MSQZl112mf99rKmpMQzDMK677rpW89673W7jlFNOMTIyMoz77rvPeOKJJ4zjjz/eyMzMNLZu3dquPhDSQ7Db7cZtt91m9OnTx0hJSTEmTJhgvPfeey2OCfQDaRiGsWfPHmP69OlGVlaWkZGRYXz3u981duzYEfB15s2bZ4wYMcJITk42hgwZYjz22GP+leciec7mZW8D/RowYEAbr1L3Fq/3QPMHR6BfbV2QAT6xdA/wM9s1onEP8DMb+zp6H/BzG58GDBgQ9H1rDuWBQrphGEZVVZUxa9YsIzc310hLSzMmT55srF27tt19MBmGYbR93B0AAABApFGTDgAAAMQYQjoAAAAQYwjpAAAAQIwhpAMAAAAxhpAOAAAAxBhCOgAAABBjCOkAAABAjCGkAwAAADGGkA4AAADEGEI6AAAAEGMI6QAAAECMIaQDAAAAMYaQDgAAAMQYQjoAAAAQYwjpAAAAQIwhpAMAAAAxhpAOAAAAxBhCOgAAABBjCOkAgE4ZOHCgTCaTli5dqt27d+uGG27Qcccdp9TUVA0fPlz33nuvHA5HV3cTAOIKIR0AEBY7duzQ+PHjNX/+fNXV1fm33XfffZo8ebLq6+u7uIcAED8I6QCAsPj1r3+tvLw8rVy5UrW1taqvr9c///lPZWRkaM2aNbr11lu7uosAEDdMhmEYXd0JAED8GjhwoHbv3i2bzaYvvvhCgwYNarH/X//6l37wgx/IbDZr9+7d6tu3bxf1FADiByPpAICwuPLKK1sFdEm6+uqrNXDgQHm9Xr322msB2y5YsEAmk0klJSUR7iUAxAdCOgAgLCZPnhx035lnnilJ2rhxY7S6AwBxjZAOAAiLwsLCY+47ePBgtLoDAHGNkA4AAADEGEI6ACAsysrKjrkvPz+/xfZPPvlEkyZN0uzZsyVJxx9/vC6//HLt27cvch0FgDhASAcAhMVHH30UdN/HH38sSTrppJP82zZv3qxzzz1XDQ0NuvzyyyVJP/3pT7V48WJNmTJFjY2Nke0wAMQwQjoAICz+/e9/a/fu3a22v/zyy9q1a5csFosuvfRS//a7775bFotFixcv1vnnny9Juvnmm/XCCy9ox44dmjt3btT6DgCxhpAOAAiLpKQkXXDBBVqzZo0kye1266WXXtKNN94oSZo1a5Z/jnSPx6NFixbpsssuU69evVqc56KLLlK/fv303nvvRfcbAIAYYu3qDgAAuoc5c+bozjvv1KRJk5SZmSmXyyWHwyFJmjhxoh555BH/sQcPHpTdbtfQoUMDnmvo0KEqLS2NSr8BIBYxkg4ACIthw4Zp/fr1uv7665WRkSGv16uhQ4fq7rvv1rJly5SRkdHVXQSAuMFIOgAgbAYMGKD58+cf87j8/HzZbDbt3Lkz4P6vv/5aI0aMCHf3ACBuMJIOAIg6i8WiqVOnqri4WJWVlS32LVq0SLt379aFF17YRb0DgK5HSAcAdInf/e53crvdOvvss/XBBx9Ikp588kldddVVGjZsmH70ox91cQ8BoOsQ0gEAXeLEE0/UokWLZLPZ9Morr0iSnnjiCU2ePFmLFy9Wenp6F/cQALoOIR0A0GXOOOMMrV69Wn//+98lSV988YWKi4v9UzUCQKLiwVEAQKeUlJR0dRcAoNthJB0AAACIMYR0AAAAIMaYDMMwuroTAAAAAL7FSDoAAAAQYwjpAAAAQIwhpAMAAAAxhpAOAAAAxBhCOgAAABBjCOkAAABAjCGkAwAAADGGkA4AAADEGEI6AAAAEGMI6QAAAECMIaQDAAAAMYaQDgAAAMSY/x93I+6n6KoKrwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "### SETUP PLOT\n",
    "mpl.rc('xtick', labelsize=12) \n",
    "mpl.rc('ytick', labelsize=12)\n",
    "mpl.rc('xtick.major', size=5, pad=10, width=1)\n",
    "mpl.rc('ytick.major', size=5, pad=10, width=1)\n",
    "mpl.rc('axes', linewidth=1)\n",
    "mpl.rc('lines', markersize=5)\n",
    "fig = plt.figure(figsize=(8, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_ylabel(r'$\\mathregular{f(\\leq p_{\\phi})}$', fontsize=17, labelpad=10)\n",
    "ax.set_xlabel(r'$\\mathregular{p_{\\phi}}$', fontsize=17, labelpad=10)\n",
    "ax.set_xlim([1e-5, 1])\n",
    "ax.set_ylim([0, 1])\n",
    "ax.set_xscale('log')\n",
    "fig.canvas.draw()\n",
    "ax.xaxis.set_major_locator(mtick.FixedLocator((ax.get_xticks()/10).tolist()[1:]))\n",
    "labels = list(['', '0.00001', '0.0001', '0.001', '0.01', '0.1', '1.0'])\n",
    "ax.set_xticklabels(labels, size=12)\n",
    "\n",
    "### PLOT SIGMA LIMITS\n",
    "plt.axvline(x=0.0455, color='red', linestyle='-', linewidth=10, alpha=0.2)\n",
    "plt.text(0.033, 0.52, r'$\\mathregular{2\\sigma}$', size=17, color='red', alpha=0.8)\n",
    "plt.axvline(x=0.0027, color='red', linestyle='-', linewidth=10, alpha=0.2)\n",
    "plt.text(0.00195, 0.52, r'$\\mathregular{3\\sigma}$', size=17, color='red', alpha=0.8)\n",
    "plt.axhline(y=0.5, color='gray', linestyle=':', linewidth=1, alpha=0.8)\n",
    "\n",
    "### PLOT TWO-SAMPLE TEST RESULTS\n",
    "### OUTPUT FROM calc_twosample.py)\n",
    "ax = plot_twosample('../output/twosample_test_usc.txt', ax, 'green', 'Upper Sco', 0.00015, 0.9)\n",
    "ax = plot_twosample('../output/twosample_test_tau.txt', ax, 'darkviolet', 'Taurus', 0.07, 0.1)\n",
    "\n",
    "### SAVE FIGURE\n",
    "fig.savefig('../output/figure_09.png', bbox_inches='tight', dpi=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d1bf8585",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close('all')"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "main_language": "python",
   "notebook_metadata_filter": "-all"
  },
  "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
