ML 500 GeV MDM 1 GeV.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "import uproot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "s500_dm1_file =  uproot.open(\"mc_392924.SlepSlep_direct_500p5_1p0_2L8.exactly2lep.root\")\n",
    "s100_dm50_file =  uproot.open(\"mc_392925.SlepSlep_direct_100p0_50p0_2L8.exactly2lep.root\")\n",
    "s200_dm100_file = uproot.open(\"mc_392936.SlepSlep_direct_200p0_100p0_2L8.exactly2lep.root\")\n",
    "sets = [s500_dm1_file, s100_dm50_file, s200_dm100_file]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "SMALL_SIZE = 14\n",
    "BIGGER_SIZE = 16\n",
    "\n",
    "plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=BIGGER_SIZE)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=BIGGER_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize\n",
    "plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title\n",
    "plt.rc('lines', linewidth = 3)\n",
    "plt.rc('figure', figsize= (7, 5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_fig(xTitle=\"xTitle\",yTitle=\"yTitle\"):\n",
    "    fig = plt.figure()\n",
    "\n",
    "    ax = plt.gca()\n",
    "    ax.set_prop_cycle(color=['#002BFF','#A15FFF','#FF0061','#FF8913'])\n",
    "    ax.set_xlabel(xTitle,labelpad=6)\n",
    "    ax.set_ylabel(yTitle)\n",
    "    ax.tick_params(direction='out', length=8, width=2)\n",
    "    ax.yaxis.offsetText.set_fontsize(10)\n",
    "    #ax.xaxis.labelpad = 15\n",
    "    #ax.yaxis.labelpad = 15\n",
    "    ax.tick_params(axis='both', which='major', labelsize=15)\n",
    "    \n",
    "    return fig, ax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mt2(file, mdm, indices = None):\n",
    "    mets = file[b'mini;1'].array(b'met_et')\n",
    "    mets_phi = file[b'mini;1'].array(b'met_phi')\n",
    "    lep_pt = file[b'mini;1'].array(b'lep_pt')\n",
    "    lep_eta = file[b'mini;1'].array(b'lep_eta')\n",
    "    lep_phi = file[b'mini;1'].array(b'lep_phi')\n",
    "    if indices:\n",
    "        mets = np.take(mets, indices)\n",
    "        mets_phi = np.take(mets_phi, indices)\n",
    "        lep_pt = np.take(lep_pt, indices)\n",
    "        lep_eta = np.take(lep_eta, indices)\n",
    "        lep_phi = np.take(lep_phi, indices)\n",
    "\n",
    "    lep_p = np.multiply(lep_pt, np.cosh(lep_eta))\n",
    "    mpts = np.sqrt(mets**2 - mdm**2)\n",
    "    lep_et = lep_pt/lep_p\n",
    "    \n",
    "    # split q in equal increments of magnitude and phi\n",
    "    splittings_mag = []\n",
    "    for mpt in mpts:\n",
    "        splittings_mag.append(np.linspace(0, mpt, 10))\n",
    "    s_phi = np.linspace(0, 2*pi, 10)\n",
    "    \n",
    "    mt2 = []\n",
    "    for i, qs in enumerate(splittings_mag):\n",
    "        #print(i)\n",
    "        max_mtq = []\n",
    "        for q in qs:\n",
    "            eq = np.sqrt(q**2 + mdm**2)\n",
    "            vec_mpt = np.array([[mpts[i]*np.cos(mets_phi[i]), mpts[i]*np.sin(mets_phi[i])]]*len(s_phi))\n",
    "            vec_q = np.column_stack([q*np.cos(s_phi), q*np.sin(s_phi)])\n",
    "            vec_minq = np.array(vec_mpt - vec_q)\n",
    "            vec_lep1 = np.array([[lep_pt[i][0]*np.cos(lep_phi[i][0]), lep_pt[i][0]*np.sin(lep_phi[i][0])]]*len(s_phi))\n",
    "            vec_lep2 = np.array([[lep_pt[i][1]*np.cos(lep_phi[i][1]), lep_pt[i][1]*np.sin(lep_phi[i][1])]]*len(s_phi))\n",
    "            eqmin = np.sqrt(np.array([sum(x) for x in np.multiply(vec_minq, vec_minq)]) + mdm**2)\n",
    "            \n",
    "            p1qt = np.array([sum(x) for x in np.multiply(vec_q, vec_lep1)])\n",
    "            p2mqt = np.array([sum(x) for x in np.multiply(vec_minq, vec_lep2)])\n",
    "            \n",
    "            mtalpha = np.sqrt(2*(lep_et[i][0]*eq - p1qt))\n",
    "            mtbeta =  np.sqrt(2*(lep_et[i][1]*eqmin - p2mqt))\n",
    "            #mtalpha = np.where(np.isnan(mtalpha), -1, mtalpha)\n",
    "            #mtbeta = np.where(np.isnan(mtbeta), -1, mtbeta)\n",
    "            \n",
    "            max_mtq.append(np.where(mtalpha >= mtbeta, mtalpha, mtbeta))\n",
    "        max_mtq = np.array(max_mtq)\n",
    "        mask = (np.isnan(max_mtq))\n",
    "        max_mtq[mask] = -1\n",
    "        max_mtq[mask] = np.max(max_mtq)+1\n",
    "        mt2.append(np.min(max_mtq))\n",
    "    return np.array(mt2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MET distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x19334e8e0f0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdEAAAFfCAYAAAARefYzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABOj0lEQVR4nO3deXwURf7/8dfHSLgvYSGCCCIqKD9AQFlckahcGoFVWTlkOb2+grqC99eDRVkFXFzFA1cXFEREEL+KR4B4sohK0CgQEOQQ5dANi4RbDfX7ozvDZDKTTCY3eT8fj36E6a6uri4gn6nq6ipzziEiIiIFd1xpF0BERKS8UhAVERGJkYKoiIhIjBRERUREYqQgKiIiEiMFURERkRgpiIpEwcyGmZnzt9PDHE8MOt4taP8LZraliMtS5HmKSGwUREUKZi/w5zD7h/jHQj0IXF7EZSiOPEUkBgqiIgWzABhsZpa9w8yqAlcCr4Umds5tdM59WZQFKI48yxszq1zaZRABBVGRgpoFNAXOD9p3ORBHmCAa2vVqZseb2YNmttHMDplZhpn928zOD0ozyMy+NLN9ZrbHzFaZ2fV55NnM70a+3szGm9kOM/vZzBaa2Ukh5almZs+Y2S4z22tmr5vZef75w/K7eTNra2ZvmtluMztoZsvMrEuYe/7BzM42s6VmdsDMNpjZDWHyO8XMZpvZf8zssJmlmdnlIWnG+eVrbWaLzGwf8Kp/7HdmNsfMMv0yzTCzPn76RD/Nk2b2o5lVCsm3hl8HD+d33yKRKIiKFMx3wMfk7NIdArwO7Ivi/DuBW4EngJ7AcOA94AQAP5i+BHwE/BH4E/AcUCeKvO8GWgAjgFuAzsDskDT/9I8/ClwBfBMmTVhm1h74xC/rtXit711Aipl1CEleC3jZv5e+wArgGTO7MCi/JsBnQFu8OukDfAG8ZmZ9whThDbx66QM85u9bAFzi3/sA4Fdgash5TwMNyN0FfjVQHa9+RWLjnNOmTVs+GzAMcBwNUruBKsCJwG9AdyDRT9Mt6LwXgC1Bn98CFuRxnduA/+ZTltA8m/nX/ShMXg5o5H8+AzgC3BGS7gk/3bB8rvsesBaID9oX5+/7v5DyOeDCoH2VgQzgn0H7/gX8B6gXcp0lQFrQ53F+freEpOvh778qZP+b/v7EoH0fAu+FpPsCSC7tf1vayvemlqhIwc3DCwq98VozO/ECTDRWAJea2QQzO9/M4sMcr2tmL5nZZWZWpwDlejvk8yr/58n+z06A+eUPNj+/jP3nvl39c4/43dLH+/mlABeEnHLAOfdB9gfn3GFgQ1BZAHoB7wB7svPz81wEtDWzWiF5vh7y+fdAVpj94e7naeBCMzvNv59zgLOBZ/O4bZF8KYiKFJBzbi/wf3hdukOA2c65I1Ge/jfgAbwuyaXALv85Xn0/74/wunCb4AWH/5hZipm1iSLv/4Z8Puz/rOL/PNH/+VNIuh+jyPsEvFbnfXhdpsHbaLzAH/z7ZHeYPA4HlQW8LtYhYfKb7B+vF3L+jpDPJwK7nXO/RnE/r+N92cl+tnwDsB1YGCatSNSOL+0CiJRTM/FafscBA6M9yf+FPxGYaGYJwGXAFKAa0N9PMx+Yb2Y18LqIJwLJZnZSAYJ1ONlBqAGwOWh/wyjO/RmvK/gpvHvPJYay7cL7IjExwvHtoZcI+bwDL3hXCgmkue7HOfermT0P3Ghmk/Cen/7dOfdbAcsskoOCqEhsluCNEP3ZObcmlgycczuB583sUqB1mOP7gLfMrDnwOF7L7D+xF5nP8ALRn4BJQfv/FEVZ95vZUrxBQF8UMphnS8Yb/LTGOXcwhvM/xWsdX44/WtcX6X6exRuAlN0drwFFUmgKoiIxcM5lUYAWaDYzewP4Cm9Qy26853K98J/Nmdl4vJbUB3gtsZOAm/EG2hQmgOKc+8bMXgYe9LteVwIX4T3bBa+lmZcxeCOTF5nZv/BagvWB9kCcc+6uAhbpfuBz4GMzexLYAtTF+0LR3Dk3Ip/7WWxm/wb+6XeHfwv0wwv0ue7HObfNzBbiBd2FzrnvC1hekVwUREVK1sd4LaVReF24W/FahRP845/hBc3H8J5D/gQsxnsWWRSuw5tZ6Q4gHnjfL8tbwJ68TnTOfeEPyHkAb0RvbbyW8RfAtIIWxDm31cw64o2+/RvwO7wu3tXAi1FmcwXeKy0T8QYZvYlXVy8Q/n7m4QVRDSiSImHOhT5mEJGKxMxuxwtCzZxzW0u7PIVlZk/hvZJ0gj8qOPjYbOAPeC3douiSlgpOLVGRCsTMLsPrLk3D6+7sgvc+6avlMYD6syzVBtbgtax74Y28nRwcQM3s90A7vMFbYxRApagoiIpULHvxZkK6C2+2nm14XbMPlGKZCmM/8BfgVLzBQpuBezj6mky25XgzSr2I986oSJFQd66IiEiMNNmCiIhIjBRERUREYqRnor769eu7Zs2alXYxRESkDFm5cmWGc+53kY4riPqaNWtGampqaRdDRETKEDP7Lq/j6s4VERGJkYKoiIhIjBRERUREYqQgKiIiEiMNLBIpAZmZmfz000/8+mvo+tEiUloqVapEgwYNqFWrVsx5KIiKFLPMzEx+/PFHGjduTNWqVTGz0i6SSIXnnOPgwYNs27YNIOZAqu5ckWL2008/0bhxY6pVq6YAKlJGmBnVqlWjcePG/PTTTzHnoyAqUsx+/fVXqlatWtrFEJEwqlatWqjHLAqiIiVALVCRsqmw/zcVREVERGJU4kHUzI43s7vMbIOZHTazH8zssZA0Zmb3mNn3ZnbQzD42s3Zh8jrTzN4zswNmtt3MxptZXIndTDF6/7Gjm4iIlE2lMTp3BnAx8FdgHdAEODMkzV3AfcDtfpoxQIqZtXbO7QQws7pACpAO9MVblPfveF8M7i3+2xApnN6jSvZ6C58q3vzHjRvH/PnzWb16dfFeSKQMKdGWqJn1AgYA3ZxzzzrnPnLOveScuycoTRW8IPqwc+5J51wK8CfAAaODsrsBqApc4Zxb4pybhheYx5hZ7C/9iAgAw4YNw8y45pprch274447MDMuu+yyAuXZrFkzzCzHdtddd+VIs3XrVnr37k316tWpX78+N998M7/88kuONKtWraJr165UrVqVxo0bM378eJxz+V4/LS2NgQMH0qhRIypXrszJJ5/MpZdeyuuvv86RI0cKdC8FlZiYiJnx0EMP5Tp21VVXYWaMHn30V1x2/ZtZ4H3GCy+8kKeeeirXQJiC5h3OLbfcQseOHalSpQoFWdFq06ZNXHPNNTRt2pTKlSvTqFEjLrzwQl588cVcf2/HopLuzh0BvO+cS88jzXlALeDV7B3Ouf3AQuCSoHSXAIucc5lB+17BC6xdi6zEIhVYkyZNmDt3Lvv37w/s++2335g1axYnn3xyTHnef//97NixI7Dde+/RjqOsrCySkpLYu3cvS5cuZc6cOcyfP5+xY8cG0mRmZtK9e3caNmzIihUreOKJJ5g8eTJTpkzJ87pvvfUWnTp1Ys+ePcyYMYP09HRee+01+vXrx4QJE9i+fXtM91MQTZo0YcaMGTkC/q5du3jzzTdp0qRJrvTdunVjx44dbNmyhcWLF9O7d28eeOABunTpkuPvJJa8Qx05coShQ4cyZMiQqO8nNTWVs88+m9WrVzN16lRWrVrFO++8w3XXXceLL77IihUros6rvCrpINoJWG9mT5pZpv8sc4GZNQpK0xLIAjaEnLvWPxacbl1wAufcVuBASDoAzMzltRXFzYkca9q0acNpp53Gq68GvtPy9ttvU6VKFRITE2PKs2bNmiQkJAS2GjVqBI4tXryYNWvWMGvWLNq3b0/37t2ZNGkSzz33HJmZ3vfl2bNnc+DAAV588UVat27NlVdeyZ133smUKVMitkb379/P8OHDSUpK4p133qFnz56ceuqpnHPOOYwYMYLU1FQaN24cSJ+enk5SUhI1a9akQYMGDBw4kJ07dwKwaNEi4uPj2bVrV45r3HPPPbRt2zbPe7/kkkvYt28fH374YWDfSy+9RKdOnWjevHmu9JUrVyYhIYHGjRvTrl07xowZw4cffsgXX3zBpEmTCpV3qKlTp3LTTTdx+umn55sWvMkKhg4dymmnncYnn3xCnz59OP3002nXrh0DBw7k/fff57zzzguk37ZtGwMGDKBu3brUrVuXpKQkNmzwfs2vX78eM2PVqlU5rvHPf/6T+vXrl+mZvko6iCYAw4B2eN26w4EOwOt2dJxxXWCfcy4r5NzdQDUziw9K93OYa+z2j4lIERg5ciTTp08PfJ4+fTrDhw+P+dWARx99lHr16tGuXTsmTJiQo8tv+fLltGrVKkfLqWfPnhw+fJiVK1cG0nTp0iXHu7c9e/Zk+/btbNmyJew1Fy9eTEZGBnfccUfEcmXfz44dO7jgggto3bo1n3/+OSkpKezbt48+ffpw5MgRunXrRr169Zg3b17gXOccc+bMYfDgwXnee6VKlRgyZEiu+hw5cmSe5wVr3bo1vXr14rXXXivyvAsiLS2N9PR0brvtNo47Lnwoya7TAwcOcOGFF1KlShU++ugjli9fzoknnki3bt04cOAAp59+Oh07dmT27Nk5zp89ezb9+/enUqVKxXIPRaGkg6j5W1/n3DvOubnAn4FzgYuC0oX7OmlhjkVKl2u/c87y2mK6G5EKYNCgQaSmprJhwwZ27txJcnIyw4YNiymvm2++mTlz5vDBBx8wevRoHnvsMW688cbA8Z07d9KwYcMc59SvX5+4uLhASzBcmuzP2WlCrV+/HoAzzjgjsG/VqlXUqFEjsGX/An/mmWdo27YtEydOpFWrVrRp04aZM2eyYsUKUlNTiYuLY8CAATl+4S9btoytW7cyaNCgfOtgxIgRLFiwgMzMTFJTU9myZQv9+vXL97xgZ555Jps2bSqWvKMVrk737NmTo07/9re/AfDKK6/gnGPGjBm0adOGli1b8uyzz7Jv3z7eeustAAYPHsycOXMCvQnff/89S5cuzfeLSWkr6dG5u4FNzrngfpB/A7/gjdB9z09T08ziQlqjdYADzrlfg/KqE+YatQnfQhWRGNStW5fLL7+c6dOnU6dOHRITE2N+HjpmzJjAn9u0aUOtWrXo378/EydOpF69ekDkl9+D94emyf7FW5DW8RlnnEFaWhoAbdu2DXQZrly5ko8//jhHN3O2jRs3cu655zJ48GAef/xxvvvuO5o2bcrs2bNJTEzM0SUcSatWrWjbti1z5swhLS2NAQMGUK1atajLDd79hrvXosi7MGrWrBmo00svvTTQy7By5Uo2b95MzZo1c6Q/cOAAGzduBGDgwIHcdtttLF26lAsuuICXX36Z5s2b07lz5xIrfyxKOoiuBSqH2W9A9tC4dUAc0AL4JihN6DPQdYQ8+zSzJkD1kHQiUkgjRoxg6NCh1KhRg/HjxxdZvp06dQLg22+/pV69eiQkJLBs2bIcaTIyMsjKygq0NhMSEnK1OLPnPg1toWbLfs63bt26wC/l+Ph4WrRoAeQMvkeOHCEpKYlHH300Vz7Z+Xfo0IGWLVvy8ssvc9tttzFv3jwmT54c9X2PGDGCp59+mk2bNrFo0aKoz8uWnp4e8TlnYfOOVnCdnn322QAcd9xxgTqNj48PpD1y5Ajt2rXjlVdeyZXPCSecAECDBg3o1q0bs2fP5oILLmD27NlcffXVxVb+olLS3blvAW3MrH7QvguASsBX/udPgEy811oAMLNqQG/g3aDz3gV6mlnwV5v+wEHgo6IvukjFdfHFFxMfH09GRgZ//OMfiyzf7FbLiSeeCEDnzp1Zu3YtP/zwQyDNkiVLqFy5Mh06dAikWbp0KYcOHcqRplGjRhFfzejRowf16tXj4YcfzrdM7du3Z82aNTRt2pQWLVrk2IJbUldffTWzZ88mOTmZ/fv3c+WVV0Z93/3792f9+vWcdNJJgS8S0Vq9ejXJyckRu2kLk3dBtGvXjlatWjFp0iSyskKHsOTUvn17vv32W+rXr5+rTrODKHhduvPmzWPlypWsWrWqzHflQskH0X8Cu4CFZtbbzAYBs4AU59y/AZxzh4BHgHvMbJSZXQzM88s6NSivacBhYIGZdTOz64BxwJSQ115EpJDMjK+//prNmzdTuXK4zqT8LV++nMcee4y0tDQ2b97Mq6++yo033kifPn0C3cM9evTgrLPOYsiQIXz55ZekpKRw++23c+211waWqho0aBDVqlVj2LBhrF69mgULFvDII48wZsyYiN251atX51//+hfJycn06tWL5ORkNm7cyKpVq5gyZQqHDh0iLs6b7GzUqFHs2bOH/v3789lnn7Fp0yZSUlK47rrr2Lt3byDPwYMHk56ezn333UefPn0KtJRWzZo12bZtG59++mme6Q4fPszOnTvZvn07X331FVOmTCExMZEOHTpw2223FSrvUN9++y1paWls376dX375hbS0NNLS0iK+62lmvPDCC2zcuJHOnTvzxhtvsH79etauXcvzzz/PDz/8EKjTq6++moYNG9K3b18++ugjNm/ezMcff8zYsWMDI3QBLr/8cn799VdGjhzJueeey2mnnVageygVzrkS3fC6ad8B9uM913wBqBuSxoD/BX7Aa1kuBc4Ok9eZwPt+mh3Ag0BcLOXq0KGDK0vem3J0k/ItPT29tIsQk6FDh7qkpKSojz/wwAPurLPOiph+5cqVrlOnTq527dquSpUq7owzznAPPPCA279/f4503333nUtKSnJVq1Z1J5xwghs9erQ7dOhQjjRff/2169Kli6tcubJLSEhw48aNc0eOHMn3nlauXOmuuuoql5CQ4I4//nh3wgknuG7durmZM2e6rKysQLr169e7K6+80tWpU8dVqVLFnX766W706NHu8OHDOfLr0qWLA9ybb76Z77W7du3qRo0aFfXxoUOHOrxBki4uLs7Vq1fPde3a1T3xxBO5ylHQvCOlyb5e8LZ58+Y8z9uwYYMbMWKEa9KkiatUqZKrVauWO//8892TTz6Z4+9t586dbtiwYe53v/udi4+Pd82aNXPDhw93//nPf3Lk9+c//9kB7oknnsjzukUpr/+jQKrLI3aYi2KWj4qgY8eOLjU1tdSun9ccuRfdWnLlkKK3du1aWrVqVdrFEJEI8vo/amYrnXMdI52rVVxERERipCAqIiISIwVRERGRGCmIioiIxEhBVEREJEYKoiIiIjFSEBUREYmRgqiIiEiMFERFRERiVNKruIiIL69ZqopDcc98NW7cOObPn8/q1auL90IiZYhaoiIS1rBhwzAzrrnmmlzH7rjjDsyMyy67rEB5TpgwgT/84Q9Ur1494mTxW7dupXfv3lSvXp369etz880355oEfdWqVXTt2pWqVavSuHFjxo8fTzRTmKalpTFw4EAaNWpE5cqVOfnkk7n00kt5/fXXOXLkSL7nF0ZiYiJmxkMPPZTr2FVXXYWZMXr06MC+7Po3MypVqkSDBg248MILeeqppwJrn8aadzjZ1wrepk2bliNNeaz34qYgKiIRNWnShLlz57J///7Avt9++41Zs2bFtDD34cOHueKKK/jLX/4S9nhWVhZJSUns3buXpUuXMmfOHObPn8/YsWMDaTIzM+nevTsNGzZkxYoVPPHEE0yePJkpU6bkee233nqLTp06sWfPHmbMmEF6ejqvvfYa/fr1Y8KECWzfvr3A91NQTZo0YcaMGTkCz65du3jzzTdp0qRJrvTdunVjx44dbNmyhcWLF9O7d28eeOABunTpkuPvJJa8w3nuuefYsWNHYBs6dGjgWHmu9+KkICoiEbVp04bTTjuNV199NbDv7bffpkqVKiQmJhY4v/HjxzN27NjAIs6hFi9ezJo1a5g1axbt27ene/fuTJo0ieeee47MTG+Fw9mzZ3PgwAFefPFFWrduzZVXXsmdd97JlClTIraK9u/fz/Dhw0lKSuKdd96hZ8+enHrqqZxzzjmMGDGC1NRUGjduHEifnp5OUlISNWvWpEGDBgwcODCwEPiiRYuIj49n165dOa5xzz330LZt2zzv/5JLLmHfvn18+OGHgX0vvfQSnTp1CrvIduXKlUlISKBx48a0a9eOMWPG8OGHH/LFF18wadKkQuUdTp06dUhISAhsVatWDRwrz/VenBRERSRPI0eOZPr06YHP06dPZ/jw4RG7Ywtj+fLltGrVKkfLqWfPnhw+fJiVK1cG0nTp0iXHL/iePXuyfft2tmzZEjbfxYsXk5GRwR133BHx2tn3s2PHDi644AJat27N559/TkpKCvv27aNPnz4cOXKEbt26Ua9ePebNmxc41znHnDlz8l1EulKlSgwZMiRXfY4cOTLP84K1bt2aXr168dprrxV53rfccgv169fnnHPOYdq0aTm6WstzvRcnBVERydOgQYNITU1lw4YN7Ny5k+TkZIYNG1Ys19q5cycNGzbMsa9+/frExcUFWiTh0mR/zk4Tav369QCcccYZgX2rVq2iRo0agW327NkAPPPMM7Rt25aJEyfSqlUr2rRpw8yZM1mxYgWpqanExcUxYMCAQHqAZcuWsXXrVgYNGpTvPY4YMYIFCxaQmZlJamoqW7ZsoV+/fvmeF+zMM89k06ZNRZr3+PHjmTt3LikpKQwYMICxY8fyt7/9LXC8vNd7cdHoXBHJU926dbn88suZPn06derUITExMabnodGK1MIN3h+aJrs7sSCt4zPOOIO0tDQA2rZtGxiss3LlSj7++GNq1KiR65yNGzdy7rnnMnjwYB5//HG+++47mjZtyuzZs0lMTMzRNRlJq1ataNu2LXPmzCEtLY0BAwZQrVq1qMsN3v2Gu9fC5H3fffcF/tyuXTuysrKYMGEC9957b2B/ea734qIgKiL5GjFiBEOHDqVGjRqMHz++2K6TkJDAsmXLcuzLyMggKysr0OpJSEjI1fL56aefAHK1lLKdfvrpAKxbt47OnTsDEB8fT4sWLYCcQeDIkSMkJSXx6KOP5sonO/8OHTrQsmVLXn75ZW677TbmzZvH5MmTo77PESNG8PTTT7Np0yYWLVoU9XnZ0tPTIz7nLGze2Tp16kRmZiY//vgjDRs2PCbqvTioO7cceP+xnJtISbv44ouJj48nIyODP/7xj8V2nc6dO7N27Vp++OGHwL4lS5ZQuXJlOnToEEizdOlSDh06lCNNo0aNaNasWdh8e/ToQb169Xj44YfzLUP79u1Zs2YNTZs2pUWLFjm2mjVrBtJdffXVzJ49m+TkZPbv38+VV14Z9X3279+f9evXc9JJJ9GpU6eozwNYvXo1ycnJEbtpC5N3sLS0NKpUqUKdOnWAY6Pei4OCqIjky8z4+uuv2bx5M5UrV445n61bt5KWlhYYiJKWlkZaWhr79u0DvF+6Z511FkOGDOHLL78kJSWF22+/nWuvvZZatWoB3jPaatWqMWzYMFavXs2CBQt45JFHGDNmTMRuxerVq/Ovf/2L5ORkevXqRXJyMhs3bmTVqlVMmTKFQ4cOERcXB8CoUaPYs2cP/fv357PPPmPTpk2kpKRw3XXXsXfv3kCegwcPJj09nfvuu48+ffoEyheNmjVrsm3bNj799NM80x0+fJidO3eyfft2vvrqK6ZMmUJiYiIdOnTgtttuK1TewRYuXMhzzz3H6tWr2bhxI88//zz3338/1113XeDv+1io9+Kg7lyRUlLcMwgVteDWQKzuv/9+XnzxxcDn7FddPvjgAxITE4mLi+Ptt9/mxhtv5A9/+ANVq1Zl0KBBObr4ateuzZIlSxg1ahQdO3akbt26jB07ljFjxuR57b59+/Lpp58yceJEhg8fTkZGBrVq1aJ9+/bMmDGDq6++GoBGjRqxbNky7r77bnr16sWhQ4c4+eST6dGjR44vEE2bNuX8889n6dKlPPjggwWui9q1a+ebJiUlhRNPPJG4uDjq1KlD69ateeCBB7j++uuJj48vVN7BKlWqxNNPP82YMWM4cuQIzZs3Z/z48YwaNSpHnsdCvRc1i2a2iYqgY8eOLjU1tdSuX5Bu2vL2y7eiW7t2La1atSrtYohIBHn9HzWzlc65jpHOVXeuiIhIjBRERUREYqQgKiIiEiMFURERkRgpiIqIiMRIr7iUIk2cICJSvqklKiIiEiMFURERkRiVaBA1s2Fm5sJsNwSlMTO7x8y+N7ODZvaxmbULk9eZZvaemR0ws+1mNt7M4kryfkREpGIrrWeiFwEHgz4HL4x3F3AfcDuwDhgDpJhZa+fcTgAzqwukAOlAX+BU4O94XwruRaQ86P14yV5v4S3Fmv24ceOYP38+q1evLtbriJQlpdWdu8I592nQ9hOAmVXBC6IPO+eedM6lAH8CHDA66PwbgKrAFc65Jc65acBfgTFmVrqzEYscI4YNG4aZcc011+Q6dscdd2BmXHbZZVHnt2XLFkaOHEnz5s2pWrUqzZs35+677+bgwYM50m3dupXevXtTvXp16tevz80338wvv/ySI82qVavo2rUrVatWpXHjxowfP55opjBNS0tj4MCBNGrUiMqVK3PyySdz6aWX8vrrr3PkyJGo7yUWiYmJmBkPPfRQrmNXXXUVZsbo0Ud/zWXXv5lRqVIlGjRowIUXXshTTz0VWIMz1rzDueWWW+jYsSNVqlSJuCpLNPX+0Ucf0aFDB6pUqULz5s2ZNm1antcNPq9Pnz787ne/o3LlyjRv3px+/frx3nvvRXV+aSlrz0TPA2oBr2bvcM7tBxYClwSluwRY5JzLDNr3Cl5g7VoC5RSpEJo0acLcuXPZv39/YN9vv/3GrFmzCrww97p168jKyuKZZ55hzZo1TJ06lZkzZ3LLLUdbyFlZWSQlJbF3716WLl3KnDlzmD9/PmPHjg2kyczMpHv37jRs2JAVK1bwxBNPMHnyZKZMmZLn9d966y06derEnj17mDFjBunp6bz22mv069ePCRMmsH379gLdTyyaNGnCjBkzcgSeXbt28eabb9KkSZNc6bt168aOHTvYsmULixcvpnfv3jzwwAN06dIlx99JLHmHOnLkCEOHDmXIkCFhj0dT75s3b+bSSy/lvPPO48svv+Tuu+/mpptu4rXXXsvz2s8++ywXXXQRderUYe7cuaxbty6w4PbNN9+cb9lLU2kF0Y1m9puZfWNm1wftbwlkARtC0q/1jwWnWxecwDm3FTgQki4gwrPYwFboOxI5BrVp04bTTjuNV18NfK/l7bffpkqVKiQmJhYor169evHCCy/Qs2dPmjdvTlJSEv/7v/+b4xfs4sWLWbNmDbNmzaJ9+/Z0796dSZMm8dxzz5GZ6X1nnj17NgcOHODFF1+kdevWXHnlldx5551MmTIlYmt0//79DB8+nKSkJN555x169uzJqaeeyjnnnMOIESNITU2lcePGgfTp6ekkJSVRs2ZNGjRowMCBAwMLUi9atIj4+Hh27dqV4xr33HMPbdu2zbMOLrnkEvbt28eHH34Y2PfSSy/RqVOnsItsV65cmYSEBBo3bky7du0YM2YMH374IV988QWTJk0qVN6hpk6dyk033RRYSDtUNPU+bdo0GjVqxNSpU2nVqhXXXnstQ4cODbvQdrbvv/+em2++mVtuuYWZM2dy0UUXccopp9C5c2dGjx6d6/HAJ598QteuXalWrRqNGzfmf/7nfwL/Np599lkaNmzIb7/9luOcQYMG0bdv33zrIBYlHUR34D3v/DPQG/gMmGZm2euS1AX2OeeyQs7bDVQzs/igdD+HyX+3f0xEisjIkSOZPn164PP06dMZPnx4xDUkCyIzM5O6dY/+l12+fDmtWrXK0XLq2bMnhw8fZuXKlYE0Xbp0oWrVqjnSbN++PbBOaajFixeTkZHBHXfcEbEs2fezY8cOLrjgAlq3bs3nn39OSkoK+/bto0+fPhw5coRu3bpRr1495s2bFzjXOcecOXMYPHhwnvdbqVIlhgwZkqs+R44cmed5wVq3bk2vXr1yte6KIu+8RFPvy5cvp0ePHjnO69mzJ6mpqbm6oLPNnz+fX375JeLfTfC/s1WrVtGjRw/69OnDV199xYIFC0hLS2PEiBGA13X9888/k5KSEjhn//79vPHGG/n+3cSqRIOoc26Rc+4h59xi59y7zrkheF2395pZdlnCfZW0MMcipQv7VdQ5Z3ltsd6TyLFu0KBBpKamsmHDBnbu3ElycjLDhg0rdL5bt27l0Ucf5cYbbwzs27lzJw0bNsyRrn79+sTFxQVaguHSZH/OThNq/fr1AJxxxhmBfatWraJGjRqBbfbs2QA888wztG3blokTJ9KqVSvatGnDzJkzWbFiBampqcTFxTFgwIBAeoBly5axdetWBg0alO99jxgxggULFpCZmUlqaipbtmyhX79++Z4X7Mwzz2TTpk259hdF3pFEU++R0vz2229kZGSEzXf9+vXUqlWLhISEwL633norx9/N0qVLAZg8eTL9+/dn7NixnHbaaXTq1IlnnnmG1157jZ9++om6dety6aWX5vi7ef311zn++OPp3bt34SshjLIwY9F84CqgGV5LsqaZxYW0RusAB5xz2V9ldvv7QtUmfAtVRGJUt25dLr/8cqZPn06dOnVITEws8PPQUD/++CM9e/ake/fu3HprzgVyI7Vwg/eHpsnuTixI6/iMM84gLS0NgLZt2wZaSitXruTjjz+mRo0auc7ZuHEj5557LoMHD+bxxx/nu+++o2nTpoHnd8FdwpG0atWKtm3bMmfOHNLS0hgwYADVqlWLutzg3W+4ey2KvPMSTb3H8ncTeuzCCy8kLS2N3bt3c+6555KV5YWDlStX8u233zJ37txc+W/cuJEGDRowePBghg0bxoEDB6hWrRqzZ8+mX79+VKlSpaC3G5WyEESzObznnHFAC+CboGOhz0DXEfLs08yaANVD0olIERgxYgRDhw6lRo0ajB8/vlB57dy5k4suuojWrVsza9asHL9AExISWLZsWY70GRkZZGVlBVo4CQkJuVqcP/30E0CuVlC27Od869ato3PnzgDEx8fTokULIOcv8SNHjpCUlBT2OV52/h06dKBly5a8/PLL3HbbbcybN4/JkydHXQcjRozg6aefZtOmTSxatCjq87Klp6dHfM5Z2LwjiabeI6U5/vjjqVevXth8Tz/9dPbs2cOOHTs48cQTAahevTotWrTI1Xo9cuQI11xzTa4vXkDgC8xll13G8ccfzxtvvMHFF19MSkoKixcvjuGOo1MWRudeCWQA3wGfAJl4r7UAYGbV8J6fvht0zrtATzOrGbSvP967px8Vd4FFKpqLL76Y+Ph4MjIy+OMf/xhzPjt27CAxMZFWrVoxZ84cjj8+5/f4zp07s3btWn744YfAviVLllC5cmU6dOgQSLN06VIOHTqUI02jRo0ivprRo0cP6tWrx8MPP5xvGdu3b8+aNWto2rQpLVq0yLHVrHn0V87VV1/N7NmzSU5OZv/+/Vx55ZVR10P//v1Zv349J510Ep06dYr6PIDVq1eTnJwcsZu2MHnnJZp679y5c47nkdlpOnbsSKVKlcLm269fPypVqlSgv5vQv5cWLVoEntVWrlyZfv36MXv2bObOnUtCQgJduxbfSxslPWPRa2Z2p5ldYmaXmdksvOA33jl3xDl3CHgEuMfMRpnZxcA8v5xTg7KaBhwGFphZNzO7DhgHTAl57UVEioCZ8fXXX7N582YqV64cUx7bt2+na9euJCQk8I9//IOMjAx27tzJzp07A911PXr04KyzzmLIkCF8+eWXpKSkcPvtt3PttddSq5b3CvigQYOoVq0aw4YNY/Xq1SxYsIBHHnmEMWPGROwyrF69Ov/6179ITk6mV69eJCcns3HjRlatWsWUKVM4dOgQcXHehGejRo1iz5499O/fn88++4xNmzaRkpLCddddx969ewN5Dh48mPT0dO677z769OkTKF80atasybZt2/j000/zTHf48GF27tzJ9u3b+eqrr5gyZQqJiYl06NCB2267rVB5h/r2229JS0tj+/bt/PLLL6SlpZGWlhZ4Rzeaer/hhhv44Ycf+Mtf/sLatWt5/vnneeGFFyKWFbxXc/7xj3/w5JNPMnjwYN5//322bNnCl19+GXh9Jvvv5s477+Tzzz/nhhtu4Msvv+Tbb7/lrbfe4vrrr8+R5+DBg1m0aBHTpk1j0KBBHHdc8YW6qLpzzex0oI5z7nP/c1XgfqA13vuaT0Z5vW+AEUATvEFA6cAQ59ysoDSP4AXNu4F6QCrQ3Tn3Y3YC59xuP8A+ifcO6c/AY3iBVKR8KOYZhIpacCssFosXL2bDhg1s2LAh1zPVzZs306xZM+Li4nj77be58cYb+cMf/kDVqlUZNGhQjq7V2rVrs2TJEkaNGkXHjh2pW7cuY8eOZcyYMXlev2/fvnz66adMnDiR4cOHk5GRQa1atWjfvj0zZszg6quvBqBRo0YsW7aMu+++m169enHo0CFOPvlkevTokeMLRNOmTTn//PNZunQpDz74YIHro3bt2vmmSUlJ4cQTTyQuLo46derQunVrHnjgAa6//nri4+MjnhdN3qGuueYaPvroaEfe2WefDRz9u4mm3k855RTeeecdbr31Vp555hkaNWrEE088kW8r/cYbb6RVq1b8/e9/56qrrmLPnj2ccMIJ/P73v+fNN9+kS5cugPfK1ccff8y9995L165dycrKonnz5lx++eU58rvgggto3Lgx6enpvPLKKwWui4KwaGb5MLPFQJpz7g7/8xS8GYRWAW2AvzjnnirOgha3jh07utTU1BK9ZqxLoV2U+3GAlGFr166lVatWpV0MEYkgr/+jZrbSOdcx0rnRtnHbAMv8DI8DhgB3Ouc6AA8B1xWoxCIiIseAaEfn1gGyp+c4G29Cg/n+5w+ByB3eEqBFuEVEji3RtkR/xHvtBKAHsNE5973/uQbwW9izREREjmHRtkTfBB42s9bAMODZoGP/j5xLmYmIiFQI0QbRu4AqQE+8gDoh6FgfYEkRl0vkmBJphhkRKV3RDK7NS1RB1F+O7NoIx84rVAlEjnGVKlXi4MGDRTr9mogUjYMHD0acCCIaUT0TNbNNZhZ2jR8za21m6s4ViaBBgwZs27aNAwcOFPpbr4gUDeccBw4cYNu2bTRo0CDmfKLtzm0GRJqmpArQNOYSiBzjsmey2b59e8TloESk5FWqVImGDRsWaLapUAWZgD7SV+iOaOUUkTzVqlWrUP9RRaRsihhE/YWys+fGccBCM/slJFlV4ASgeOdVEhERKYPyaoluAt7z/zwUbw7b/4SkOYw3/+3zRV80ERGRsi1iEHXOvQG8AYG19sY75zaXULlERETKvGhfcRle3AUREREpb6IeWGRmzYGrgJPxRuQGc865kUVZMBERkbIu2vVE+3J0ceyf8J6FBtPLbyIiUuFE2xJ9CG+1lqudc6GDi0RERCqkaINoc2CsAqiIiMhR0QbRdUC94iyIRC90XdKLbg2fTkREile064neAdzjDy4SERERom+JjsNria41sw3Af0OOO+dc16IsmIiISFkXbRDNAr4pzoKIiIiUN9FOtpBYzOUQEREpd6J9JioiIiIhog6iZtbYzKaYWaqZbTaz1v7+v5hZp+IrooiISNkUVRA1s7OAVcCfge14U//F+4ebArcUS+lERETKsGhbon8H1gKnAFcAFnTsE+D3RVwuERGRMi/a0bnnAwOdc/vMLC7k2I9AQtEWS4pF78dzfl6oDgQRkcKItiV6JI9j9YGDRVAWERGRciXaIPo5EGlN0auAZUVTHBERkfIj2u7cB4EUM1sMvIy39Fk3M7sFuBy4oJjKJyIiUmZFO9nCR2b2R+AfwHR/9yPAFuCPzrnPYrm4mTXGmwmpOlDTObfP32/A3cD/4HUXrwBuds6lhZx/JjAV6Az8DDwP/NU5lxVLecqLNs+HPNu8Vc82RURKQ7QtUZxzbwNvm1kLoAGwyzlX2KkAJwP78IJosLuA+4Db8VaQGYPXEm7tnNsJYGZ1gRQgHegLnIo3ivg44N5ClktERCRf0b4n2ib7z865b51znxQ2gJpZF6AX8GjI/ip4QfRh59yTzrkU4E94Xcijg5LeAFQFrnDOLXHOTQP+Cowxs1qFKZuIiEg0oh1YlGZmX5nZWDM7sbAX9V+TmQqMBzJCDp8H1AJezd7hnNsPLAQuCUp3CbDIOZcZtO8VvMCqFWVERKTYRRtEBwLfAQ8DW81skZldbWbVYrzuDUAV4Kkwx1rirRqzIWT/Wv9YcLp1wQmcc1uBAyHpADAzl9cW432IiEgFFu3AornAXDOrjxdQBwOzgH1m9jowy+92zZeZ1cMb7TvYOferN4Yoh7rAvjCDg3YD1cws3jn3i5/u5zCX2O0fE8g9wYKIiBSZAq3i4pzLcM5Ndc51As7AG63bHUguQDYTgM+cc+/kdakw+yzMsUjpcu13zlleW9SlFxER8UU9OjeY3417rr81AH6L8ryzgBHABWZWx9+d3SVc28yy8FqSNc0sLqQ1Wgc44Jz71f+8298XqjbhW6giIiJFqiBLoZmZdTezmcBOvO7cmsAoop879zSgErAcLwju5uhz0R/wBhutA+KAFiHnhj4DXUfIs08za4L3ukyOZ6UiIiLFIaqWqJlNBgYBJwIb8d7HnOWc21TA6/0buDBkXy/gTuBSYBPeAKZMvNdaHvKvXw3oDfwz6Lx3gdvNrKZzbq+/rz/ePL4fFbBcIiIiBRZtd+4IvFdOZjnnPon1Ys65DODD4H1m1sz/49KgGYseAe4zs90cnWzhOLyWarZpwM3AAjObCDQHxgFTQl57ERERKRbRBtET/RGxJeURvKB5N1APSAW6O+d+zE7gnNttZhcDT+K9Q/oz8BheIBURESl20b7i8gsEZi66AC+wPeuc2+lPA/hjUJdqgTjnXgBeCNnn8EbxTsjn3HTgoliuKyIiUljRPhOtDLwEXMHRV0gW4g0wmgSsx5uqT0REpMKItjt3AtAN+DOwBPgx6Ni7wI0oiBabXKu2iIhImRBtEB0I3Ouce9mf9zbYZqBZkZZKRESkHIj2PdF6eHPXRsqjctEUR0REpPyItiW6GW/h6/fDHDsXb2FtKS3B8+Mu1ALdIiIlJdqW6EzgLjO7Goj39zkzuxC4FZheHIUTEREpy6INopOAt/Gm+vuvv+/fQAqQ7JybGulEERGRY1W074lmAQPM7CmgJ96k87vwAqim2CtLtPSZiEiJKdAqLs65pcDSYiqLiIhIuVKg9URFRETkKAVRERGRGCmIioiIxEhBVEREJEYRg6iZLfBXaMHMhphZvZIrloiISNmXV0u0L3CC/+cZwKnFXxwREZHyI68g+iPeVH9wdPkzERER8eX1nuirwGNmNgUvgH5qZpHSOudcgd45lTJIc/CKiBRIXoHvVmAZcCbwAPACsK0EyiRoDVERkfIgYhB1zjlgHoCZDQMed859VULlkpKgKQJFRAol2rlzTynugoiIiJQ3Ub8namYnmtmjZrbCzDaa2edmNsnMEoqzgCIiImVVVEHUzE4HvgJuBvYBnwP7gVuANDM7rdhKKCIiUkZFO6J2IrAHONc5tyV7p5k1BRb7x68o8tJJVDI2Hf1z/ealVw4RkYom2u7cC4H7ggMogHPuO2Ccf1xERKRCiTaIxgN7Ixzb6x8XERGpUKINomnATWaWI715sy/c6B8XERGpUKJ9JjoeeAtYa2ZzgR1AAvAn4DQgqXiKJyIiUnZF+55ospldBjwE/C9H59JdCVzmnFtcfEUUEREpm6Ke79Y5lwwkm1k1oC6w2zl3oNhKVsFomj8RkfKnwItyO+cOOOe2xRJAzayfmX1iZrvM7JCZfWNm95pZfFAaM7N7zOx7MztoZh+bWbsweZ1pZu+Z2QEz225m480srqBlEhERiVVJr7xSD/gAmAz8DJyL94pMAjDaT3MXcB9wO7AOGAOkmFlr59xOADOrC6QA6Xjrnp4K/B3vS8G9JXMrIiJS0ZVoEHXOPRuy6wMzqwWMMrObgMp4QfRh59yTAGa2HNiCF2SzA+QNQFXgCudcJrDEz2ecmU3y95UJ7z9W2iUQEZHiUuDu3GKwi6PvmZ4H1MJbyxQA59x+YCFwSdA5lwCLQoLlK3iBtWuxllZERMRXKkHUzOLMrJqZnY83H+8z/tJrLYEsYEPIKWv9Y9la4nX1BjjntgIHQtIFX9PltRXRrYmISAWSbxA1s3gz+8LMehThdff721LgI7znn+CN+t3nnMsKSb8bqBY0AKku3jPVULv9YyIiIsUu3yDqnPsFOAX4rQivex7QBRiLNzDoyeBLhklvYY5FShe2Vemcs7y2gt+CiIhUdNEOLFoC9ADeL4qLOue+8P/4bzPLAF40s7/jtSRrmllcSGu0DnDAOfer/3m3vy9UbcK3UKWgeoe8t7rwltIph4hIGRZtEJ0KvGRmxwP/hzftX44Wn3NuU5jzopEdUE/Be84ZB7QAvglKE/oMdB0hzz7NrAlQPSSdiIhIsYl2YNFHQBO8dzY/AtbjDf4J3mL1B//nZuATIBNvTl4A/BmSegPvBp3zLtDTzGoG7esPHPTLJyIiUuyibYkOL4qLmVky3iQJa/BG4f4B77noXOfcRj/NI8B9Zrabo5MtHIfXGs42DW9U7wIzmwg0x5u0YUpZekc0P5rqT0SkfIt2AvoXi+h6K4BhQDO8gUqbgLvxgmK2R/CC5t14MxylAt2dcz8GlWe3mV2MNyBpId5z0MfwAqmIiEiJKNCMRf56omfiBzd/IoSoOefuw5vSL680Dpjgb3mlSwcuKsj1RUREilLUky2Y2ShgJ/A13ijdM/z9/2dmNxdP8URERMquqIKomV0LPI43Mvcqjr63Cd6ECVcWeclERETKuGhbomOAvzvnrgNeDzm2Dr9VKiIiUpFEG0RPARZFOLaf8BMfiIiIHNOiDaIZeCNqwzkD2FYkpRERESlHog2iC4H7zax50D5nZvWBW/GelYqIiFQo0QbRe4HDwGq8yRIc8ATeEmVZwPhiKZ2IiEgZFu1kC7vMrCPwF6AnsNE/90ngsfI0S5DEKHhCek1GLyICFGCyBefcXuBBfxMREanwCjpjUS2gNdAY+AFYo1aoiIhUVFEHUTO7H2+y+BocnWxhr5lNds49VByFExERKcuiCqJm9le8OW+fB14BfgQaAgOBv5rZ8c65ccVVSBERkbIo2pbotXgzFt0etG8N8L6Z7QGuQyuoiIhIBRPtKy61iTxjUbJ/XEREpEKJNoh+BpwT4dg5/nEREZEKJWJ3rr92aLabgdfN7DdgHkefiV4FjAD6FmchRUREyqK8non+hjczUTYDHvE3QvZ/nU9eUkIyNuX8XL95+HQiIlJ4eQW+8eQMoiIiIhIkYhDVKysiIiJ5i3ZgkYiIiIQoyIxFrYB+QBOgSshh55wbWpQFExERKeuinbFoCDAd7xnpT8AvIUn07FRERCqcaFui9wFvACOdcz8XX3FERETKj2iDaAJwgwKoiIjIUdEOLFoGtCrOgoiIiJQ30bZERwMLzGwXsBjYHZrAOXekKAsmIiJS1kUbRH8AvgReinDcFSAvERGRY0K0ge85oD/wf8A6co/OFRERqXCiDaJ9gdudc48XZ2FERETKk2gHFu0H0gt7MTP7k5m9aWbbzGyfma00s4EhaczM7jGz783soJl9bGbtwuR1ppm9Z2YHzGy7mY03s7jCllFERCRa0QbRGcCgIrjeGGAfcCvQB/gAeNnMbgpKcxfee6kTgd5++hQzS8hOYGZ1gRS8Z7F98SbLHwv8tQjKKCIiEpVou3O/Awaa2RIgmfCjc6dHkU9v51xG0Of3zawRXnCdamZV8ILow865JwHMbDmwBW+E8L3+eTcAVYErnHOZwBIzqwWMM7NJ/j4REZFiFW0Qfcb/2RS4OMxxhzctYJ5CAmi2Lzm6qPd5QC3g1aBz9pvZQuASjgbRS4BFIcHyFbzWa1dgYX5lERERKaxog+gpxViG8zj6vLUlkAVsCEmzFm90MEHp3g9O4JzbamYH/GMKoiIiUuyiCqLOue+K4+JmdjFeK3SEv6susM85lxWSdDdQzczinXO/+Ol+DpPlbv9YuGvlOUl+hw4dClByERGRUlxP1MyaAS8DbzjnXgg6FC7YWZhjkdJpRRkRESkR0S6Ftpl8gpNzrnm0FzWzE4B3ga3A4KBDu4GaZhYX0hqtAxxwzv0alK5OmKxrE76FinPOwu3P1rFjx2My+GZsOvrn+lH/DYmISDSifSb6EbmDaD2855n7CHk+mRczqwa8BcQDSc65/UGH1wFxQAvgm6D9Lf1jwelahuTbBKgekk5ERKTYRPtMdFi4/WZWB++Vl5Ro8jGz44F5wGnAH5xzP4Uk+QTIBP4EPOSfUw3vfdF/BqV7F7jdzGo65/b6+/oDB/ECvhSn3iETVy28pXTKISJSygo1abxz7mczmwxMwHu+mZ+ngUuBW4ATzOz3Qce+dM4dMrNHgPvMbDdeq3IM3rPbqUFppwE3460sMxFoDowDpugdURERKSlFsfLKIeCkKNP28H+Gm4P3FLxJFR7BC5p343UZpwLdnXM/Zid0zu32R/Y+ifc6y8/AY3iBtMxq87ymHhYROZbEHET9rtnWeIFrTTTnOOeaRZHG4bVsJ+STLh24KJrrioiIFIdoR+ceIfLo3EwgqchKJCIiUk5E2xIdT+4geghvTt13nXN7irRUIiIi5UC0o3PHFXM5REREyp1Sm7FIRESkvIvYEjWz+wuSkXNufOGLIyIiUn7k1Z07Lorzg5+TKoiKiEiFklcQrZTPue3wZhXqSe6ly6Qi0QxGIlJBRXwm6pzLCrfhzQ70EvAZcCZwnf9TRESkQol6sgV/gvcHgCF4q6jcBjztr+8pIiJS4eQbRM2sAfC/eC3OQ3jPPh8LWX1FRESkwslrdG5t4E7gJrzFrh8HJjrndpdQ2URERMq0vFqim/EWuV6MN4BoB1DXzOqGS+yc2xRuv4iIyLEqryBax//Zk6Orr+QlrtClERERKUfyCqLDS6wUIiIi5VDEIOqce7EkCyIiIlLeFMWi3CJlVu9ROT8vfKp0yiEixyZNQC8iIhIjBVEREZEYqTtXygV1y4pIWaQgWoFkhLzJW7956ZRDRORYoSAqhfb5qpyfzy2dYhSaWrsiUlAKolLulXbwK+3ri0jpURCVCkUBT0SKkoKolBnBAa4sBLe8yhMajIvjGiJS9imISqkpqkAkIlJaFERFoqCALyLhKIgWsfcfK+0SFI1cI27/XwFO7v340T8vvCXnoSiDkYKWiJQHCqISleCgml9AzZG2eIpTpmnwkkjFoSAq5ZJaqiJSFiiIFrM2zz+efyIpNWUpGBekLGrdipQNJT4BvZm1MLNnzewrM8sysw/DpDEzu8fMvjezg2b2sZm1C5PuTDN7z8wOmNl2MxtvZnElcR8iIiKlsYrLWcClwHp/C+cu4D5gItAb2AekmFlCdgIzqwukAA7oC4wHxgJ/LbaSi4iIBCmN7tyFzrk3AMxsPlA/+KCZVcELog8755709y0HtgCjgXv9pDcAVYErnHOZwBIzqwWMM7NJ/j4pZWWpu1REpKiVeBB1zh3JJ8l5QC3g1aBz9pvZQuASjgbRS4BFIcHyFbzWa1dgYZEVuoIIfa2lKNz3ds5nwg8m3RIhZcVQEl8qNDpYpOSUxYFFLYEsYEPI/rVA/5B07wcncM5tNbMD/rEcQdTMXF4X7dChQ6zlLTd+3pvzc/3wyaQcUotfpHSUxSBaF9jnnMsK2b8bqGZm8c65X/x0P4c5f7d/TPLx7VchO6J8Ql4cLdaiVNYCSlkrj4gUndIYWBSNcK1GC3MsUrpc+51zltdWBGUWEZEKpiwG0d1AzTCvqtQBDjjnfg1KVyfM+bUJ30IVEREpUmWxO3cdEAe0AL4J2t/SPxacrmXwiWbWBKgekk6kQtNAI5HiUxaD6CdAJvAn4CEAM6uG977oP4PSvQvcbmY1nXPZQ2b6AweBj0quuCIlT89ZRcqGEg+ifkC81P/YGKhlZv38z+845w6Y2SPAfWa2G69VOQav63lqUFbTgJuBBWY2EWgOjAOmVLR3RINH3dapWXrlEBGpaEqjJdoAmBeyL/vzKXiTKjyCFzTvBuoBqUB359yP2Sc453ab2cXAk3ivs/wMPIYXSCus0NdYylpQDX5vtKK/M1oWqKtXpHBKY7KFLRwdaRspjQMm+Fte6dKBi4qscCISoAArkr+y+ExU8hHa2iyqtFIx6HmqSNEpi6+4iIiIlAtqiYpIgFqpIgWjlqiIiEiM1BKVgBOC1tf5byl8vdLIXREpbxREpdSELpMm5YdG7op41J0rIiISIwVRERGRGKk7t4wq7an8gp+PQuk8I5XyQ927UlEpiBaxNs/rOZ8cm/T6i0huCqLlgGYdEhEpmxREywgFypxCR+7qlRcRKYv0pEtERCRGaolKVEp7IgYRkbJIQVTKBXXvikhZpCBaivQcNHaaIrD8KMjrL3pVRsobdcyJiIjESC1RESlyeb1TGnxMLU0p7xRERaTUaAIHKe/UnSsiIhIjtUSlwMravLoauXvsKqquX3UhS3FREC1BGo0rUnQUGKUsUHeuiIhIjNQSlULLq3u3NLp+83qHVF2/IlKU1BIVERGJkVqickwLbXlK+RXtu6egZ6RSchREi5kGE5Uf6uqtGGKdhlCBWcJREJUiF/octCxTS1WiDZT5BV8F3IpJQVRKVLRLqpW1d1GlfNFMSFJSynUQNbMzgalAZ+Bn4Hngr865rNIsl8SmrLVg82qlBnf1qhtYpOIqt0HUzOoCKUA60Bc4Ffg73ojje0uxaBKlogqapf0aTUHSKsBWDMU10EldxmVPuQ2iwA1AVeAK51wmsMTMagHjzGySv0+kxOUVYNVqrZhKOvgVJogrUBdMeQ6ilwCLQoLlK8BEoCuwsFRKJWVaeXrWqoB7bMrveW1w4Cqq13r0ClDxMedcaZchJmb2E/C0c25cyP79wDjn3OSQ/XneaIcOHUhNTS10ub49WaM9jxXRBtiyEJjzekYb7Xmh5+Y321Ne+eQlOJ/Quru1t74oHCui/TJQkHxClUSr2cxWOuc6RjxejoPor8Dtzrl/hOz/AZjpnLsnZH9+N5oBfFfIYnXwf64sZD7HItVNZKqbyFQ3kaluIivKumnqnPtdpIPluTsXIFxgtHD7nXNW3IXJDtR5fWupqFQ3kaluIlPdRKa6iawk66YMPxHK126gTpj9tfFedxERESlW5TmIrgNaBu8wsyZAdf+YiIhIsSrPQfRdoKeZ1Qza1x84CHxUOkUSEZGKpDwH0WnAYWCBmXUzs+uAccAUvSMqIiIlodyOzoXAtH9PknPav3GlNe1f0MPsYh/EVN6obiJT3USmuolMdRNZSdZNuQ6iZY3+UUemuolMdROZ6iYy1U1kCqIiIiLlQHl+JioiIlKqFERFRERipCAqIiISIwXRImBmZ5rZe2Z2wMy2m9l4M4sr7XIVJzP7k5m9aWbbzGyfma00s4EhaczM7jGz783soJl9bGbtwuR1zNafmTX268eZWY2g/RW2bszseDO7y8w2mNlhM/vBzB4LSVMh68fMBpjZF/6/mW1mNtPMGoWkOebrxsxamNmzZvaVmWWZ2Ydh0hRZPUSbV1jOOW2F2IC6wHa8BcK7461zuh94qLTLVsz3vRx4GbgKuAh4FG/O4puC0tyNN/nFaKAb8A7eRP8JFaX+/Dra6ddNDdWNA5jl39f1eMsWDgb+FpKmwtUP0Mf/d/IkcLFfL1uAL4DjKlLdAH2B74F5wFrgwzBpiqweoskrYllLu7LK++ZX/m6gVtC+O4ADwfuOtQ2oH2bfy8Bm/89VgD3A/UHHqwP/Cf4HfCzXH9AF+C9wG0FBtCLXDdAL+BU4M480FbJ+8NZDXhmyLzuwtqpIdUPOLw3zQ4NoUdZDtHlF2tSdW3iRFgevivct+5jknMsIs/tLoIH/5/OAWsCrQefsx1ss/ZKgc47J+vO7i6YC4/G+0QaryHUzAnjfOZeeR5qKWj+V8H6ZB/vZ/5n9vmOFqBvn3JF8khRlPUSbV1gKooXXkpAJ751zW/G+6bQMe8ax6zwg+5djSyAL2BCSZi056+VYrb8b8L7hhlsquCLXTSdgvZk9aWaZ/nOqBSHP/Spq/UwHupjZEDOrZWanAw8BHwR96aiodROqKOsh2rzCUhAtvLqEX3ptt3+sQjCzi/GeY2QHjbrAPpd7CsbdQDUziw9K93OYLMtt/ZlZPeBBYIxz7tcwSSps3QAJwDCgHTAAGI63gPLrZpbd2qqQ9eOcexuvbv6J1yL9BogDrghKViHrJoyirIdo8wqrvC/KXVZEvTj4scjMmuE9D33DOfdC0KFI9RJ67FirvwnAZ865d/JIU1Hrxvytr3NuF4CZ7cBbeeki4D0/XYWrHzO7EG9hjcfxVqlqiLeoxutm1i3ol3yFq5sIirIeos0rFwXRwqvQi4Ob2Ql4/+G34o0mzLYbqGlmcSHf8OoAB4JaaMdU/ZnZWXjP/S4wszr+7mr+z9pmlkUFrRvfbmBTdgD1/Rv4BTgTL4hW1Pr5O/Cmc+7O7B1mlobXHdkXWEDFrZtQRVkP0eYVlrpzC6/CLg5uZtWAt4B4IMl/GJ9tHV5XVIuQ00KfURxr9Xca3gCR5Xj/OXdztIv7B7zBRhW1bsB7zhSOAdmDSSpq/bQE0oJ3OOe+wXv14lR/V0Wtm1BFWQ/R5hWWgmjhVcjFwc3seLx3uE4DLnHO/RSS5BMgE/hT0DnVgN54dZbtWKu/fwMXhmwT/WOXApOpuHUD3peuNmZWP2jfBXhfPL7yP1fU+vkOaB+8w8xa4Y0k3eLvqqh1E6oo6yHavMIr7feByvuG91B6B7AE7yXd64B9lJOXmgtx3//Ee1ZwM/D7kK2yn+ZuvFFwo/BeHn8b73WPhhWp/vAGiwTeE63IdYP3KsFWvJZ6b2AQ3kv1S0LSVbj6AW7Ba43/3b+fq/EGF20GqlekusF7BNLP35YDa4I+Vyvqeogmr4hlLe3KOhY2vGc57+N9u9mBNzIzrrTLVcz3vMUPDOG2Zn4aA/4XrxvzILAUOLui1R/hg2iFrRu8brN38GaO2Q28ANQNSVPh6se/5/8BvvbrZhswF2he0eoGaFaSv1+izSvcpvVERUREYqRnoiIiIjFSEBUREYmRgqiIiEiMFERFRERipCAqIiISIwVRERGRGCmIioiIxEhBVEREJEYKoiIiIjFSEBUpBWY2zMycv50e5nhi0PFuMeT/FzO7Iv+UOc4Za2ZfBy2Onb2/vplNMLNVZrbPzA6Z2UYzm2lmiQW8xu/M7FczezqPNCP9+040z5dmdntBriNSUhRERUrXXuDPYfYP8Y/F6i9A1EHUX/v0HmC8C5oL1Mxa462uMhyY4+d5Cd5qNM2BD8ysYbTXcc79B29ljP5mFh8h2RC8FU0+8ssyHrjHX7tWpExREBUpXQuAwcGtPzOrClwJvFaC5RgJ/Aq8HlSOSn759gNtnXN/c84tds594Jyb5pw7H2+lkTwXLQ7jReAEICn0gJk1A7oAM4OC+ZvAIeCaAl5HpNgpiIqUrllAU+D8oH2X4y0SHDaImllXM3vPzPaa2X4zW+S3GLOPb/HzvDqoS/iFfMpxDTDXOZcVtO9KvPVi7/RbkLk45152zv03pHxXmNmnZnbAzH42s3lmdnJQkoXAfwnfAv8z3ooaM4OukYW3dq2CqJQ5CqIipes74GNyBpQheC3CfaGJzSwJeM8/NhhvPc6awFIza+InuxzYCSwCOvvbg5EK4Ae4lnjLPwW7GMgCkqO9GTO7AS/4p+Ot/Xg90Br4KHthZOfcL8ArQFKYLtrBwCfOuW9D9n8MnGZmzaMti0hJUBAVKX0zgT+ZWRUzOxFv8eCZEdI+jvessK9z7g3n3BtAL7xgNxbAOfclcBjIcM596m8b87h+Z//nVyH7TwL+45w7GLzTzI4zs+ODNvP31wAmAjOccyOcc+845+biPUNtjNdlnO1FIB7oH5Tv74HT/WOhvvR//j6P+xApcQqiIqVvHlAZ6I33jHEnXmszBzM7DTgVmB0cxIADwHLgghiv38j/Gdpla6EJfe/gPQfN3rKDY2egVpjy/QCsCy6fc+5zYC25W+CHgVfDXDO7bI3CHBMpNQqiIqXMObcX+D+8gDIEmO2cOxImaQP/57/IGcR+BS4D6sVYhCr+z8Mh+78HfucPdAp2E3AO0CdC+VLClO//hSnfTKCzmbXwR+r2B95wzv0cpozZreHQsoiUquNLuwAiAngB5W28L7YDI6TZ5f+8Gy9Qhfolxmtn51uXo8EK4H28wTy9CBq165zbAIGRtOHyGQasCXOd0Fd2ZgET8L48pOGN2A3XlYt/DCAjwnGRUqEgKlI2LMHrxvzZORcuAAF8A2wBznLOPZJPfoeJvtW2zv/ZHNgetP81YCMw0cz+HWmEbpBP8AJlC+dcpGAY4JzbZmbv4Q0magP8CCyOkPwU/+c3+eUrUpIUREXKAP81jkgt0Ow0zsxGAW/43Z+v4rXMGgLnAVudc1P85OlAFzO7DO8Za4ZzbkuErD/HC7rnAv8Out4v/qxHi4A0M3sKWIHX4k3AewUG/Bamcy7Tn1noKTP7Hd6kCnvwBhV1BT50zr0ccu0XgZfwguRjzrnfIpSxE1638Kd5VJFIidMzUZFyxDn3Dt4AnerA83gBbhJeUFselPRuvFbbq3iBb1weeR4C3sAb2BR67Gu8VuJMvNdp3vCv+Te8rt+u/gjc7PTP4j0rPQOvu/Zd4K94X9jTwlx+AZBJyLuhYVwGvOmcO5BHGpESZ0EzfIlIBeXPgfs+0Mw5t7V0S5OTmTUCtgI9nXO5Ri2LlCYFUREBwMwWA+udc6NLuyzBzOwxvGkHLyrtsoiEUneuiGS7GfghdBWXMmAHMKq0CyESjlqiIiIiMVJLVEREJEYKoiIiIjFSEBUREYmRgqiIiEiMFERFRERi9P8BrLEeRGd9wsAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "make_fig(\"Met (GeV)\", \"Number of events\")\n",
    "title(\"Missing energy\")\n",
    "nbins   = 100 # number of bins in histogram\n",
    "xlow    = 0\n",
    "xhigh   = 1000\n",
    "bins    = np.linspace(xlow, xhigh, nbins)\n",
    "for i in range(len(sets)):\n",
    "    n, bins, patches = plt.hist(sets[i][b'mini;1'].array('met_et')/1e3, bins, alpha = 0.7)\n",
    "legend([\"Ml 500 Gev MDM 1 Gev\", \"Ml 100 Gev MDM 50 Gev\", \"Ml 200 Gev MDM 100 Gev\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MT2 Distribution for (500 GeV, 1GeV)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We build something called \"stransverse mass\". Recall regular transverse mass for a 1-> 2 decay:\n",
    "\\begin{align}\n",
    "m_T = \\sqrt{2 (E_{T1}E_{T2} - \\mathbf{p_{T1}\\cdot p_{T2}})}\n",
    "\\end{align}\n",
    "or \n",
    "\\begin{align}\n",
    "m_T = \\sqrt{2 (E_{T1}E_{T2} - |p_{T1}||p_{T2}|\\cos(\\Delta \\phi)}\n",
    "\\end{align}\n",
    "Where \n",
    "\\begin{align}\n",
    "E_T = E\\frac{|p_T|}{|p|}\n",
    "\\end{align}\n",
    "and\n",
    "\\begin{align}\n",
    "m_T \\leq m\n",
    "\\end{align}\n",
    "\n",
    "Here, we have two $\\tilde l-> l + DM$, where the $DM$ goes into $MET$. However, we don't know how the $MET$ is split up between the two unknown $DM$ particles, which means we have to split $\\mathbf{p_{T}}^{miss}$ into all possible combinations $\\mathbf{q_T}$ and $\\mathbf{p_{T}}^{miss} - \\mathbf{q_T}$. Then, we define the the $m_{T2}$ variable (\"stransverse mass\") as:\n",
    "\\begin{align}\n",
    "m_{T2} = \\text{min}_{\\mathbf{q_T}}\\bigg \\{\\text{max} \\big[m_T(\\mathbf{p_T^{\\alpha}}, \\mathbf{q_T}, MDM), m_T(\\mathbf{p_T^{\\beta}}, \\mathbf{p_{T}}^{miss} - \\mathbf{q_T}, MDM)]\\big]\\bigg\\}\n",
    "\\end{align}\n",
    "\n",
    "Where $\\alpha$ and $\\beta$ indicate the two respective $\\tilde{l} -> l + DM$ decays we are looking at.\n",
    "\n",
    "Note that we assume $|p_T| \\approx |p|$.\n",
    "\n",
    "$m_{T2}$ satisfies the inequality:\n",
    "\\begin{align}\n",
    "m_{T2}^2 \\leq m_{\\tilde l}^2 - MDM^2\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Wenjie\\AppData\\Roaming\\Python\\Python37\\site-packages\\ipykernel_launcher.py:40: RuntimeWarning: invalid value encountered in sqrt\n",
      "C:\\Users\\Wenjie\\AppData\\Roaming\\Python\\Python37\\site-packages\\ipykernel_launcher.py:39: RuntimeWarning: invalid value encountered in sqrt\n"
     ]
    }
   ],
   "source": [
    "mt2_0 = mt2(sets[0], 1e3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x19330f4a5c0>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAc8AAAFfCAYAAAAoHkdoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1MUlEQVR4nO3debwcVZ338c+XsAazQQZQiYbNYSI6OkRQXFACImpEUUAQkMGNRxAfQEZAcMImAQUeTHBwQwiCCMhiUIQkCIggCCIwJEFWAYEoeEPIwhZ+zx+n+qaodN9bfW+vud/369Wv233qdNXv3O3Xp+rUOYoIzMzMrLzV2h2AmZlZt3HyNDMzq5OTp5mZWZ2cPM3MzOrk5GlmZlYnJ08zM7M6OXlaTZL2lxTZ401Vtr8/t31HSeNzr/t6XC9pmKSvSbpO0gJJz0n6k6TPSSr1eynpAEn3S3pR0sKGfwNWHCffri9W2b5uFn9IOjFXXvn+je9n/1NqfJ+uqFL3zZKulbRY0jOSfiJpvSr1xkm6VNKzkhZJukzSG+po8xskfVfSfZKWSVoqaZ6ksyX9R9n9NELh9+yDVbaPl/RKtv3zufL8729IWiLpEUmXS9qj2u9Zru63qmyTpIey7T8tEfdhkmZKejJ7z5Q6272GpP8j6XeSeiS9lO3rKkn7Slq9zv3dJenePrZvNpA4hyonTyvjOWDfKuX7ZdsqngTeVXgAnFso+zKwDnAM8L/AF4GPA78Ffgic0l9Akl4H/AC4GdgB2LGuFg1Mre/DJ4FG3DD9Hl79ffqv/MaszdeTvnefAg4itfuqfCKQNBy4DtgS+GwW8xbAbyWt218Qkt4P3AN8CPgf4CPAZNL3+13A7wfcwsHp6/dwcR/v250U94eBY4EXgJ8B10pap8ZxPiNJhfL3AuOBJSXj/QKwAXBFyfq9JI0g/T2cDtxJ+jlOAg4ntfUnpN+7epwHTJC0dY3t+2VfZ9Qb75AUEX74UfUB7E9KCucCDwPKbVsHeJb0RxzAjjX2EcCJVcqHAetVKT8HeB5Yp5/Yts/2vUOD2rpGvn2FbeNz34dXgE0K22fnvg8n5sor37/x/Rx7SlZv9X7qnQEsBEbnyt6XvXe3XNlXgeXA5rmyTYCXgcP6Ocb6wNOkBDm8ynYBX27x7+H7c9//xcC6he33577/n6/y/d+8yj4/mf0sp1X5fT0v+/69v7DtR6SE9gjw0xJxr5Z9XT3b75Q62nxu9newbY3tbwe2q/P7uCHwEnBmje0PAje28mfbzQ/3PK2M84E3knpGFZ8gJcBfDGSHEbE8Iv5ZZdMfgbWAsbXeK+lcUg8MYE52quncbNsakk7MTs+9mH09UdIaufdXTsN+WdKpkp4g9UZG9xP2TcBDwD65fW0MfIDWfFr/GPCriFhYKYiIG4FHgV0L9f4QEQ/k6j1MSoj5etV8gZRAD46IpcWNkXyvWC7pC9lpweclPS3px/nTyZLulbTS74qkbbOfxcf7iQvgMrIPCrn3bwdsRvodLS0ifgFcCXwh66nnPQbcQK6XK2ltUm+/9M85Il6pJ6bcsV5P+h37fkTcWmPfd0bEzYX3/bukX2aneJdJ+r2k9+beswC4BtireMo3q7cp7nWW5uRpZfwVuJFXnzLbD7icvk+XDcT2pN7Vk33UOQE4JHt+EOmU3AnZ6/OAI0n/BD5K6pF8PSsv+gbwJtJp40+QPun356e8+vuwD/A4K5L5YDwmabmkv0o6JX9KMXu+Cek0d9G9wITc6zeXrFfNJODJiLizbNCSpgLfI/XAPwYcQTrle7WkYVm184GPShpTePs+wD+BX5c41FLSh7Xi7+HvSR9q6vVr0ge1iVW2zQA+lfsZfJx0duLSARynXu8nfTC9quwbsuvQNwPrkT4AfRJ4BphdOE17HvAvwM6FXewLLAMuGXDUQ4yTp5U1A9hd0tqSXku61tbQT6mSdgb2AL4TES/XqhcRDwLzspdzI+IPEfGgpK2AvUinTr8ZEddGxHHAcaRP228t7GoB8ImIuCoiroyIZSXCnAFsIemd2et9SafwBnPN8wFSwv8sKelcDBwK/DJXZwzplGlPlff/k/RPs2K9PuoVk1fRxqQPS6+iNMBr9dxDWfl4UrI8ISIOz77nPyH9HLchXSsFuIB0+nKP3D7XAD4N/DwiXuwnrooZwCRJr5e0Vra/gf4ePpp9fW2VbZdm8VZ66vsBV0TEc1XqNtrG2ddH84XZgKX8zyD///vbWf0dIuLSiPg16QPhQ6TrvBW/JP1u7Jfb71qk68JXRMSzjW/OqsnJ08q6hPQpfTLwGeApYE6jdi5pAmkQx/WUGDBUw/uyr8WRkJXX2xfKr6g36UXEQ6Sezr6SJpJ6coP6EBERP42IU7LEMysijiAlpB0lVQZCVQavVIu3OLClnnpl69xLul5WeUzKynci/R+5IP+PHbgVWET2M4mIlU6Fkj4ojKW+799vST39vUm/i+uQPmwMRM3vaUQsJp1Z2VfSRsAH64xzMGr9DL7Oq38GM6D3rMT2pL/RV3I/A5HOBlT+LoiIF4CfAx+TNCor3pV0ycKnbOvg5GmlZJ+4ryD989sPuGCg13SKJG0KzCINSvp4X73OflR6X8VTvk8VtlOjXlkzgD2BzwO3RcR9A9xPX36WfX1H9rWH9E9+pdtSSL3J/PXjnj7qVeuR5j0GVLul5ZNZLAcWyjfIvj7Aq/+xvwSMJF0/rZgBvFvSJtnrfYEHIuIP/cTUK/uwc0H23s8CvxxEb2lc9rXW78EMUtI8FPg7KRG1wmPZ1+LP4VzSz+AdvDrm9UineY9l5Z/BwcCYQi/1PGBtUm8T0t/zk6S/QSuprvuEbMibAfyK9KFrr0bsMBtwM4fUS/lQRCwaxO4qCWQj0shBcq8hXQPKG+ip1ouBM0nXlg7pp+5gpSGgEUslPUK6nlk0gdSrq7i3j3pz+znedaQe79vz1z0j4l4ASa8p1K98Tz9I9cSc/57/AjgL2EfSmaSe48n9xFPNDOAoUhs/NoD3V3yEdJ37jhrbZ5OS5teA0yNi+SCOVY/rSSOBP0ouoUXEU2QfBCXlT3MvzOqfRY3eY/6DbkT8QdJ9pF71laTrn2e0sH2rBCdPq8csUuJYWPlnOhiS/oUVn+Z3ioh/DHKXlQTyaeCkXPlnsq83DnL/AETEQkknk24XuKgR+6yiEnN+tOUvgc9KGlXpbUl6D2kk9C8L9b4jadPsNHPl2uS7SddW+/JD0injaZI+WG3EbcEs0j/uN0REnz2XiHgu+2e9L/AEqfdT1yjZbD/zJZ1FGvhyTb3vB5C0GynxnlmrjRHxiqQTSKeXzxnIcQYiIv4m6QLgS5IurDXiNld/iaTfAf8O/KnkGaEZwInA0aQ84FO2dXLytNKyT6aN6nGuQ/rHNx44ANg464VWzK23FxoR90r6GTAlu+ZzM2kk7rHAzyLi7kbEnh3r+EbtS9KdpH9e95F6mjsBXwF+ExG/zVX9Nml06i+z5D0KOBW4jXR9ruKHpNN1V0o6JtvnCaTTgd/vK5aIeFrS7tn+/izpe8DdpAQ5jnSKL8gmCsgGap0CTJf0r6QPMM9ndXcCflRowwzS79BxwE3ZLTR1i4iD66j+NkljgTVJp0I/SjplOYvUg+3rOGcDZw8kxuya+HhWXB6bIOlT2fNf9/PB5GBWTGzxQ9KHzMrp+PeRzqbkBy8dRvpweI2kH5NOw44F/gMYFhHFD03nk34nvkpKuNVGZ1sfnDytXTYk9dwgXcMq+gADu/3js6QRhgeQZjB6gjQA6bgB7KtV7iP9s3wt6drVg8DxpMTYK+uRfIA068wvgBdJ9yoeXjgtt0TSDqRJFc4nDRyZA/zfbCBMnyJiTjYy+WukW4HGkRLmI6SfyRER8edc/aMlzcvqHpTVfSw75v2F3c8inXp8fdbGVqjcfvE86TTsn0hnJy4d5Cjp/hxM+n2s2J0V1xk3IX0/q4qIRZK2J10a2Dvbz7qkCSzuAD5H7qxHRPxJ0juA/wa+S/pg9Q9SW1dK/hHxmKTfkgZ+udc5AGru746Zmdmqx6NtzczM6uTkaWZmVicnTzMzszo5eZqZmdXJydPMzKxOvlUlM3bs2Bg/fny7wzCzLrJg0YqFeDYcuXYbI7FmueOOO56OiH8pljt5ZsaPH8/tt9/e7jDMrIvc8/iKaXXfsvGoPmpat5K00ipD4ORpZjZgTphDl695mpmZ1anlyTNba+5ISfdLekHS45LOKNSRpKMlPSZpmaQbJb2tyr4mSJojaamkJyQdn1u53szMrCnacdr2J6T5FI8D5pPmzZxQqHMkaTLvI7I6hwGzJW2VLcuDpDGkyZLnkhZz3Qw4jfSB4JjmN8PMzIaqliZPSR8iTcj87xFRdV1BSWuTkufJETE9K7uFNInywaxIjAeSVpHfLVt9Y5akkaQVNU4d5LqQZmb9Gn/kr3qfPzL1I22MxFqt1adtDwCuq5U4M9uRVqC/uFIQEUuAmcAuuXq7ANcUkuRFpIS6fcMiNjMzK2h18twW+Iuk6ZIWZdcqL5P0ulydLYHlrLyU0bxsW77e/HyFiHgUWFqoB4Ck6OvRiMaZmdnQ0OrkuRGwP/A20unb/wS2Bi6XpKzOGGBxtvByXg8wXNKauXoLqxyjJ9tmZmbWFK0eMKTssWtEPAMg6UnS6vM7kBbPhbSYbrX3FrfVqrdSeUSoSt1eEydOdO/TzMxKaXXPswe4p5I4MzcBL7JixG0PMKLKLSejgaUR8VKu3ugqxxhF9R6pmZlZQ7S65zkPWKtKuYBXsufzgWHA5sB9uTrFa5zzKVzblDQOWLdQrytMPqjv7TPPak0cZmbWv1b3PK8C3ippbK7sfcAawF3Z65uBRcDulQqShgOTgatz77sa2FnSiFzZnsAy0mlgMzOzpmh18vwB8AwwU9JkSXsD5wOzI+ImgIh4HpgKHC3pIEmTgEuyWKfl9nU28AJwmaQdJX0RmAKc7ns8zcysmVp62jYiFknaAfgu6Z7MF4ErgUMLVaeSkuVRwPrA7cBOEbEgt6+eLLFOJ90DuhA4g5RAzczMmqbl0/NFxAPAh/upE8BJ2aOvenNJo3TNzMxaxkuSmZkN0K1HT2p3CNYmTp5mZgO04ci12x2CtYnX8zQzM6uTk6eZmVmdfNrWzGyAFix6vve5T+EOLU6eZmYDtO235vQ+93qeQ4tP25qZmdXJydPMzKxOTp5mZmZ1cvI0MzOrk5OnmZlZnZw8zczM6uTkaWZmVicnTzMzszo5eZqZmdXJydPMzKxOnp7PzGyAPCXf0OXk2UKTD2p3BGZm1gg+bWtmZlYn9zy7RF+91plntS4OMzNz8jQzG7B7Hn+29/lbNh7Vxkis1Zw8zcwGaPL0m3qfe/DQ0OJrnmZmZnVy8jQzM6uTk6eZmVmdnDzNzMzq5ORpZmZWJydPMzOzOjl5mpmZ1cnJ08zMrE5OnmZmZnVqafKUtL+kqPI4MFdHko6W9JikZZJulPS2KvuaIGmOpKWSnpB0vKRhrWyPmQ1tG4xYq/dhQ0u7pufbAViWe/1Q7vmRwLHAEcB84DBgtqStIuIpAEljgNnAXGBXYDPgNNKHgWOaHr2ZGXDbN3ZsdwjWJu1Knn+MiMXFQklrk5LnyRExPSu7BXgEOJgVifFAYB1gt4hYBMySNBKYIunUrMzMzKwpOu2a53bASODiSkFELAFmArvk6u0CXFNIkheREur2LYjTzMyGsHYlzwclvSzpPklfypVvCSwH7i/Un5dty9ebn68QEY8CSwv1etW41tr7GHSLzMxsyGj1adsnSdczbwOGAXsBZ0saHhFnAGOAxRGxvPC+HmC4pDUj4sWs3sIq++/JtpmZNd3suQt6n+84YcM2RmKt1tLkGRHXANfkiq6WtBZwjKQzK9WqvFVVttWqV7UXGRGqVl4xceJE9z7NrC6fn3F773Ov5zm0dMI1z0uB9YDxpJ7jiCq3nIwGlkbES9nrnqysaBTVe6RmZmYN0wnJsyJI1zGHAZsXthWvcc6ncG1T0jhg3UI9MzOzhuuE5PlJ4Gngr8DNwCJg98pGScOBycDVufdcDewsaUSubE/SvaM3NDtgMzMb2lp6zVPSL0iDhe4m9TD3zB6HRMQrwPOSpgLHSuphxSQJqwHTcrs6GzgEuEzSKcCmwBTgdN/jaWZmzdbq0bb3AQcA40iDe+YC+0XE+bk6U0nJ8ihgfeB2YKeI6B3WFhE9kiYB00n3gC4EziAlUDMzs6Zq9Wjbo4Gj+6kTwEnZo696c0nT/JmZmbVUJ1zzNDMz6ypOnmZmZnVy8jQzM6tTu1ZVMTPrelu9fmS7Q7A2KZU8Jb0JGB0Rt2Wv1wG+CWxFWt1kevNCNDPrTFd95b3tDsHapOxp2+nAp3KvTwIOB14HnCHpoEYHZmZm1qnKJs+3Ar8HkLQasB/w9YjYGjgR+GJzwjMzM+s8ZZPnaOCZ7PnbSct+XZq9vp40w4+ZmdmQUHbA0ALSZO03AR8EHoyIx7JtrwFebkJsZmYd7cJbH+19vve2b2hjJNZqZZPnL4GTJW0F7A98P7ftLcBDDY7LzKzjHX35Pb3PnTyHlrLJ80hgbWBnUiLNT533MWBWg+MyMzPrWKWSZ0QsAb5QY9t2DY3I6ja5j7HOM89qXRxmZkNFqQFDkh6S9O81tm0lyadtzcxsyCg72nY8sFaNbWsDb2xINGZmZl2gnrlto0b5RNJ6mmZmZkNCzWuekg4FDs1eBjBT0ouFausA6wEXNSc8MzOzztPXgKGHgDnZ888CtwP/KNR5AZgL/KjxoZmZmXWmmskzIq4ErgSQBHB8RDzcorjMzMw6VtlbVf6z2YGYmZl1i9LreUraFNgDeANphG1eRMTnGhmYmVmnm7TlBu0Owdqk7HqeuwKXkEbn/p10rTOv1khcM7NV1o/3f0e7Q7A2KdvzPJG0espnIqI4aMjMzGxIKZs8NwUOd+I0MzMrP0nCfGD9ZgZiZmbWLcr2PP8L+H+Sbo0Iz2NrZgacMesvvc8P3elNbYzEWq1s8pxC6nnOk3Q/8M/C9oiI7RsZmJlZpztzzv29z508h5ayyXM5cF8zAzEzM+sWZSdJeH+T4zAzM+sa9ayqYmZmZtSRPCW9XtLpkm6X9LCkrbLy/ytp2+aFaGZm1llKJU9JbwbuAfYFniBN0bdmtvmNwFcHcvAsIS+WFJJekyuXpKMlPSZpmaQbJb2tyvsnSJojaamkJyQdL2nYQGIxMzMrq2zP8zRgHrAJsBug3LabgXcO8PjfBhZXKT8SOBY4BZic1ZktaaNKBUljgNmkqQF3BY4HDgeOG2AsZmZmpZRNnu8BpkbEYlaex3YBsNHKb+mbpPcCHwK+Uyhfm5Q8T46I6RExG9g9O+7BuaoHkhbj3i0iZkXE2aTEeZikkfXGY2ZmVlbZ5PlKH9vGAsvqOWh2anUaqbf4dGHzdsBI4OJKQUQsAWYCu+Tq7QJcExGLcmUXkRKq7zk1M7OmKZs8bwNqrem5B/D7Oo97IGlZs7OqbNuSdF/p/YXyedm2fL35+QoR8SiwtFDPzMysocpOknAC6ZrjtcCFpFOoO0r6KvAJ4H1lDyhp/Wx/+0TES5KKVcYAiyNieaG8Bxguac2IeDGrt7DKIXqybcXj9rls2tZbb12uAWZmmb22GdfuEKxNyk6ScIOkjwP/DzgnK54KPAJ8PCJureOYJwG3RsSv+zpklTJV2VarntcXNbOmO3m3t7Y7BGuTsj1PIuJXwK8kbQ5sADwTEXVN2Zfd8nIA8D5Jo7Pi4dnXUZKWk3qOIyQNK/Q+RwNLI+Kl7HVPVlY0iio90ohYqYubN3HiRCdcMzMrpVTylPTWiLgbICIeAB4Y4PG2ANYAbqmy7XHgx6TTwsOAzXn1fLrFa5zzKVzblDQOWLdQz8zMrKHKDhj6s6S7JB0u6bWDON5NwAcKj1OybR8m3fd5M7CIdHsKAJKGk+73vDq3r6uBnSWNyJXtSRr5e8MgYjQzM+tT2dO2ewGfAU4Gpkq6DpgBXB4RS8seLCKeBq7Pl0kanz39XXYfKZKmAsdK6iH1Ig8jJfppubeeDRwCXCbpFGBT0tJppxduXzEza4qjLru797mvfw4tZQcM/Rz4uaSxpES6D3A+sFjS5cD52WQGjTKVlCyPIq0jejuwU0QsyMXUI2kSMJ10D+hC4AxSAjUza7qf3fZY73Mnz6Gl9IAh6O05TgOmSdqCNNft50m90rr2ldvnucC5hbIgjco9qZ/3zgV2GMhxzczMBmpAS5Jl1yC3yR4bAC83MigzM7NOVs+SZJK0k6QZwFOk07YjgIMYwNy2ZmZm3arsrSrfBvYGXgs8SFpl5fyIeKiJsVkDTD6o9raZ1SZHNDOzfpW9TnkAaaL28yPi5ibGY2Zm1vHKJs/XZvPJmpmZDXllb1V5EdJMQ6RJ4NcHvh8RT2XT9S2IiOeaF6aZmVnnKHvNcy3gp8BurJh4fSZp4NCpwF9IC1ibmZmt8sqOtj0J2JF0X+eGrFjhBLJp8hocl5mZWceqZ3q+YyLiQknDCtseBsY3NCozsy7w1UlbtDsEa5OyyXN9YF6NbasBazUmHDOz7nHoTm9qdwjWJmVP2z4MvKvGtm149dJhZmZmq7SyyXMGcKSkzwBrZmUh6QPAocA5zQjOzMysE5VNnqcCvyJNyffPrOwmYDbwm4iYVuuNZmZmq5qy93kuBz4t6SzSyNoNgGdIidMLT5vZkPS5c//Y+/zH+7+jjZFYq9W7JNnvgN81KRYzs64yZ/7f2x2CtcmAliQzMzMbypw8zczM6uTkaWZmVicnTzMzszrVTJ6SLstWTEHSfpLWb11YZmZmnauvnueuwHrZ858AmzU/HDMzs87XV/JcwIop+SrLkJmZmQ15fSXPi4EzJC0nJc4/SFpe4/Fya8I1MzNrv74mSTgU+D0wAfhv4Fzgby2IyczMrKPVTJ4REcAlAJL2B86MiLtaFJeZWcf71ife0u4QrE3Kzm27SbMDMTPrNntv+4Z2h2BtUvo+T0mvlfQdSX+U9KCk2ySdKmmjZgZoZmbWaUolT0lvAu4CDgEWA7cBS4CvAn+WtEXTIjQzM+swZVdVOQV4FtgmIh6pFEp6I3Bttn23hkdnZmbWgcomzw8AB+YTJ0BE/FXSFOB7DY7LzKzjfXTaihUar/rKe9sYibVa2eS5JvBcjW3PZdvNzIaU//3bonaHYG1SdsDQn4GvSHpVfUkCvpxt75ekT0m6WdIzkp6XdJ+kYyStmasjSUdLekzSMkk3SnpblX1NkDRH0lJJT0g6XtKwku0xMzMbsLI9z+OBq4B5kn4OPAlsBOwObAF8pOR+1gd+C3wbWAhsA0zJ9nVwVudI4FjgCGA+cBgwW9JWEfEUgKQxwGxgLmkO3s2A00gfBo4pGYuZmdmAlL3P8zeSPgqcCHyDFXPd3gF8NCKuLbmf7xeKfitpJHCQpK8Aa5GS58kRMR1A0i3AI6TkWkmMBwLrALtFxCJgVrafKZJOzcrMzMyaovR9nhHxm4iYCIwAxgEjImKbiLhmkDE8w4prptsBI0nz6laOuwSYCeySe88uwDWFJHkRKaFuP8h4zMzM+lT3YtgRsTQi/hYRSwd6UEnDJA2X9B7SvaP/k00HuCWwHLi/8JZ52baKLUmndPNxPQosLdTLHzP6egy0LWZmNvSUvebZaEtIp2gBZpCubwKMARZHxPJC/R5guKQ1I+LFrN7CKvvtybaZmZk1Td09zwbZDngvcDhpwM/03LZqvUBV2VarXtVeZESor0f9TTAzs6GqLT3PiPhT9vQmSU8D50k6jdRzHCFpWKH3ORpYGhEvZa97srKiUVTvkZqZmTVMu3qeeZVEugnpOuYwYPNCneI1zvkUrm1KGgesW6hnZmbWcP32PLMJDP4AHFn2lpQ6vTv7+jBpse1FpPtHT8yOPxyYDPwg956rgSMkjYiIysxHewLLgBuaEKOZ2Up+tN/EdodgbdJv8oyIFyVtArw82INJ+g1pcoN7SaNq30267vnziHgwqzMVOFZSDysmSVgNmJbb1dmkUbqXSToF2JQ02cLpvsfTzFplxwkbtjsEa5Oy1zxnAR8Erhvk8f4I7A+MJyXjh4CjSMmwYiopWR5FmpHodmCniFhQqRARPZImkQYazSRd5zyDlEDNzMyaqmzynAb8VNLqwBWk6fleNao1Ih7qbycRcSxp6r2+6gRwUvboq95cYIf+jmlmZtZoZZNn5TriYcChNep4UnYzMxsSyibP/2xqFGZmXWibk2b3Pr/tGzu2MRJrtbITw5/X7EBWFZMPancEZtYqf3/uhXaHYG1S132eklaTtJWk7SWt26ygzMzMOlnpGYYkHQT8NzCWNFjoHcCfJF0BXBcR321KhNY0/fWSZ57VmjjMzLpNqZ6npC8AZ5JG2u7BirlmAX4HfLLhkZmZmXWosqdtDwNOi4gvApcXts0H/rWhUZmZmXWwsslzE6DWotdLqD5Ju5mZ2SqpbPJ8mjQrUDX/SpqT1szMbEgomzxnAt+UtGmuLCSNJU2acEWjAzMzM+tUZZPnMcALwP+SJnYP4LvAPNIE78c3JTozM7MOVCp5RsQzwETgZGAN4EHSbS7TgXdFxLNNi9DMzKzDlL7PM1s384TsYWY25M08+D3tDsHapHTyBJA0EtgKeD3wOHCv1880s6HqLRuPancI1ib1zDD0TdLC1a9hxSQJz0n6dkSc2IzgzMzMOlGp5CnpONI6nD8CLgIWABsCewHHSVo9IqY0K0gzM7NOUrbn+QXSDENH5MruBa6T9CzwRWBKg2MzMzPrSGWT5yhqzzD0G+D/NCYcM7PuMf7IX/U+f2TqR9oYibVa2fs8byWtolLNO7LtZmZmQ0LNnqekfGI9BLhc0svAJay45rkHcACwazODNDMz6yR9nbZ9mTSTUIWAqdmDQvnd/ezLzMxsldFXwjueVydPMzMzo4/k6VtPzMzMqis7YMjMzMwy9cww9G/Ap4BxwNqFzRERn21kYGZmZp2q7AxD+wHnkK6B/h14sVDF10bNzGzIKNvzPBa4EvhcRCxsXjhmZmadr2zy3Ag40InTzMysfPL8PfBvwJwmxmJm1lVuPXpSu0OwNimbPA8GLpP0DHAt0FOsEBGvNDIwM7NOt+HI4thJGyrKJs/HgTuBn9bYHnXsy7rE5INqb5t5VuviMDPrNGUT3g+BPYErgPmsPNq2FEm7A/sCW5NWarkP+E5E/CxXR8BRpJVaxgJ/BA6JiD8X9jUBmAa8C1hIWmv0uIhYPpDYzMzMyiqbPHcFjoiIMwd5vMOAh4FDgaeBDwMXShobEdOyOkeSRvceQUrUhwGzJW0VEU8BSBoDzAbmZrFtBpxGmvThmEHGaGZWyoJFz/c+9yncoaVs8lxCSlSDNTkins69vk7S60gJcpqktUnJ8+SImA4g6RbgEdJ110piPBBYB9gtIhYBsySNBKZIOjUrMzNrqm2/tWIMpdfzHFrKTs/3E2DvwR6skDgr7gQ2yJ5vB4wELs69ZwkwE9gl955dgGsKSfIiUkLdfrBxmpmZ9aVsz/OvwF6SZgG/ofpo23MGGMN2rOjVbgksB+4v1JlHuuZKrt51heM/Kmlptm1m8SCS+pwFaeutt64vajMzG7LKJs//yb6+Eah2Y1OQpu+ri6RJpGuWB2RFY4DFVQb99ADDJa0ZES9m9RZW2WVPts3MzKxpyibPTRp9YEnjgQuBKyPi3Nymaj1EVdlWq17VHmZEqFp5xcSJEz0/r5mZlVIqeUbEXxt5UEnrAVcDjwL75Db1ACMkDSv0PkcDSyPipVy90VV2PYrqPVIzM7OGafl6npKGA1cBawIfyQYEVcwHhgGbF962ZbYtX2/Lwn7HAesW6pmZmTVc2SXJHqafZcciYtMS+1kduATYAnh3RPy9UOVmYBGwO3Bi9p7hwGTgB7l6VwNHSBoREc9lZXsCy4Ab+m2QmZnZIJS95nkDKyfP9UkjZRdTGPnah++RJkb4KrCepHfmtt0ZEc9LmgocK6mHFZMkrEaaTajibOAQ0ny7pwCbAlOA032Pp5mZNVvZa577VyuXNJp068rsksf7YPa12kxFm5AmQ5hKSpZHkRL07cBOEbEgF09PNlJ3Oum2lIXAGaQEamZm1lSDmsw9IhZK+jZwEmnkbH/1x5eoE9n+Tuqn3lxgh3KRmpmZNU4jVkJ5Hti4AfsxM+sqnpJv6Bpw8swG/2xFOlV6b6MCMjMz63RlR9u+Qu3RtosAf/wyM7Mho2zP83hWTp7Pk+a8vToinm1oVGZmZh2s7GjbKU2Ow8ys69zz+Ip+w1s2HtXGSKzVGjFgyMxsSJo8/abe5x48NLTUTJ6SvlnPjiLi+MGHY91i8kG1t808q3VxmJm1Q189zykl3p+/DurkaWZmQ0JfE8Ov0c/jHcC1pGXAHmhumGZmZp2jZvKMiOXVHqR5ZH8K3ApMAL6YfTUzMxsSSg8Yypb8+m9gP9J6ml8DvhcRLzYpNjMzs47Ub/KUtAHwDVIP83nStc0zCutwmpmZDRl9jbYdBXwd+ArpuuaZwCkR0dOi2MzMzDpSXz3Ph4FRpEFBJwJPAmMkjalWOSIeanx41o18G4uZrer6Sp6js687s2Idzr4MG3Q0ZmZmXaCv5PmfLYvChoy+eqXgnql1lw1GrNXuEKxNaibPiDivlYGYgU/5Wne57Rs7tjsEa5O+JkkwMzOzKpw8zczM6uTkaWZmVicvSWZmNkCz5y7ofb7jhA3bGIm1mpOnmdkAfX7G7b3PvZ7n0OLTtmZmZnVy8jQzM6uTk6eZmVmdnDzNzMzq5AFDtkrwzERm1kpOntY1+psX18ysVZw8bZXnyejNrNGcPG3Ia0aP1gnZbNXW8gFDkjaX9H1Jd0laLun6KnUk6WhJj0laJulGSW+rUm+CpDmSlkp6QtLxkryuqJmZNVU7ep5vBj4M/AFYs0adI4FjgSOA+cBhwGxJW0XEUwCSxgCzgbnArsBmwGmkDwTHNLMBZmYAW71+ZLtDsDZpR/KcGRFXAki6FBib3yhpbVLyPDkipmdltwCPAAezIjEeCKwD7BYRi4BZkkYCUySdmpWZmTXNVV95b7tDsDZp+WnbiHilnyrbASOBi3PvWQLMBHbJ1dsFuKaQJC8iJdTtGxOtmZnZyjpxkoQtgeXA/YXyedm2fL35+QoR8SiwtFDPzMysoToxeY4BFkfE8kJ5DzBc0pq5egurvL8n2/YqkqKvRyMbYGZmq7ZOvVWlWjJTlW216jkZWlv53tKh4cJbH+19vve2b2hjJNZqnZg8e4ARkoYVep+jgaUR8VKu3ugq7x9FlR5pRGjlqitMnDjRCdfM6nL05ff0PnfyHFo68bTtfGAYsHmhvHiNcz6Fa5uSxgHrFuqZmZk1VCcmz5uBRcDulQJJw4HJwNW5elcDO0sakSvbE1gG3NCCOM3MbIhq+WnbLBF+OHv5emCkpE9lr38dEUslTQWOldTDikkSVgOm5XZ1NnAIcJmkU4BNgSnA6b7H0zqdV4Ex627tuOa5AXBJoazyehPSZAhTScnyKGB94HZgp4hYUHlDRPRImgRMJ90DuhA4g5RAzczMmqblyTMiHmHFyNladQI4KXv0VW8usEPDgjPrcO6xmnWGThxtazaked1Ss87XiQOGzMzMOpqTp5mZWZ2cPM3MzOrka55mZgM0acsN2h2CtYmTp5nZAP14/3e0OwRrE5+2NTMzq5OTp5mZWZ2cPM3MzOrka55mZgN0xqy/9D4/dKc3tTESazUnTzOzATpzzv29z508hxaftjUzM6uTk6eZmVmdnDzNzMzq5ORpZmZWJw8YMltF9LeUmdf7NGsc9zzNzMzq5J6nmQ2Ye7s2VDl5mg0R/SW6WtqRAPuK1QnZOoFP25qZmdXJPU8z69NAe6xDwV7bjGt3CNYmTp5m1jSr+unXk3d7a7tDsDbxaVszM7M6OXmamZnVyadtzayrdOLtMav66WlbmZOnmQ0ZjU5yR112NwCPrwUbv+Drn0OJk6eZWT9qJd17RjyWnqzp5DnUOHma2SploLfW+JYcq4eTp5m1hZNV5+nE68mdysnTzGwV4wFMzefkaWbWRM3qYXdaEhxqCburk6ekCcA04F3AQuBHwHERsbydcZmZdarBJPNmfBDo1g8XXZs8JY0BZgNzgV2BzYDTSBM/HNPG0MzMmq6brhl3U6xldW3yBA4E1gF2i4hFwCxJI4Epkk7NyszMzBqum6fn2wW4ppAkLyIl1O3bE5KZmQ0F3Zw8twTm5wsi4lFgabbtVSRFX48WxWxmZquAbj5tO4Y0SKioJ9tWlzvuuONpSX8dRDxbV3Y1iH10CrelM7ktnae3HYP559EhVpWfCcDW+h7QmLa8sVphNydPgGo9RlUrjwg1M5BK7zUiJjbzOK3gtnQmt6XzrCrtALelXt182rYHGF2lfBTVe6RmZmYN0c3Jcz6Fa5uSxgHrUrgWamZm1kjdnDyvBnaWNCJXtiewDLihPSGZmdlQ0M3J82zgBeAySTtK+iIwBTjd93iamVkzKaJ779LIpuebzqun55vSjun5cheomzowqRXcls7ktnSeVaUd4LbUfYxuTp6dxL94nclt6UyrSltWlXaA21L3MZw8zczM6tPN1zzNzMzawsnTzMysTk6eZmZmdXLybABJEyTNkbRU0hOSjpc0rN1x5UnaXNL3Jd0labmk66vUkaSjJT0maZmkGyW9rUq9trVX0u6Sfinpb5IWS7pD0l7d1o7s+J+SdLOkZyQ9L+k+ScdIWrPb2lKI5fXZzyYkvSZX3vFtkbR/jcUjDuymduRiWF3SkZLul/SCpMclnVGo0/HtkXR9jZ9LSHpXW9oREX4M4kGahP4J0sLcO5HWGV0CnNju2Apx7go8BlwCzAOur1LnKNIkEwcDOwK/Bp4GNuqU9gK3ABcCewA7AN8hzWX8lW5qRxbDl4CTgE8AHwC+nsU9vdvaUmjXhcBT2c/lNd3UFmD/LO4PAO/MPTbopnbk4jg/i+NLpKUa9wG+VajT8e0BJhR+Hu8ErgX+Aazejna0/A9rVXtkP7AeYGSu7L9IS6ONbFdcVeJcLff8UgrJE1gbeBb4Zq5s3eyX88RcWVvbC4ytUnYh8HA3taOP9p1EumdZ3dgW4L3AP4GvkUue3dIWViTP19TY3hXtyI73IeAlYEIfdbqmPYW418x+z/6nXe3wadvB64pFuSPilX6qbAeMBC7OvWcJMJPUxoq2tjcinq5SfCewQfa8K9rRh2dI/xigy9qSnfqaBhxP+sSf11Vt6UM3teMA4LqImNtHnW5qT96HSL3In2WvW94OJ8/Bq2tR7g62JbAcuL9QPo9Xt6MT27sdUPkH0XXtkDRM0nBJ7wEOIX2aDrqvLQeSegBnVdnWbW15UNLL2XXoLxXi65Z2bAv8RdJ0SYuya3yXSXpdIc5uaU/ep4G/Ab/LXre8Hd2+nmcnaOii3G00BlgcK09t2AMMl7RmRLxIh7VX0iTS9dwDsqJubMcSYK3s+QzgiOx517RF0vrACcA+EfGStNLELt3SlieBY4HbgGHAXsDZkoZHxBl0TzsANiKdhr6LlGxGAKcCl0t6Z/YBrZvaA4Ck4cBk4AdZG6AN7XDybIxq0zRVXZS7w9VqR3FbR7RX0njS9c4rI+Lc3Kauagep5zwc2Ab4Jmm+5i9n27qlLScBt0bEr/uo0/FtiYhrgGtyRVdLWgs4RtKZ/cRX3Nbun4myx64R8QyApCdJq07tAMzJ6nVLeyomA69hxSnbipa2w8lz8HpYNRbl7gFGSBpW+PQ2GlgaES/l6o2u8v6WtlfSeqRl6R4ljSCs6Kp2AETEn7KnN0l6GjhP0ml0SVskvZnU83+fpEocwysxSFpOl7SlhktJo7vH013t6AEeqiTOzE3Ai6TRq3PorvZUfBp4ICJuz5W1vB2+5jl4q8qi3PNJp6k2L5QXrxG0vb3ZaZurSANrPpINDMjH1xXtqKGSSDehe9qyBbAG6TainuxRue75OGkQUbe0pS9Bd7VjXo1yAZUBhN3UHiSNIg36KfY6W94OJ8/BW1UW5b4ZWATsXinIXVu4Olevre2VtDrpXtUtgF0i4u+FKl3Rjj68O/v6MN3TlptI90XmH6dk2z4MfJvuaUs1nySNHv4r3dWOq4C3ShqbK3sf6YPOXdnrbmoPpHui12Ll5Nn6drT6/pxV7UG6yPwkMIt0Y+4XgcV03iQJw4FPZY9bgHtzr4dndY4ijTo7CJgE/Ir0T2PDTmkv8ANSD+AQVr5peq1uaUcWw29I90PuAnwQOC6L4aJcna5oS5W27U/1SRI6ui3AL0iTVewCfJQ0yUCw8iQcHd2OLIaRpMsat5CSyN6kiVJmFep1RXtyfzN/rrGtpe1oyx/WqvYgXT+4jvTp5UnSqMNh7Y6rEOP47J9Atcf4rI6Ab5BOtS0jDQN/eye1F3hkVWhHdvwTgP/N/ngXkk7ZfgVYI1enK9pSJZ79WTl5dnxbgG8B92X/hJcBdwD7Fup0fDtyMWxOmmlnCel0+rnAmG5sDzCWNOnDkTW2t7QdXs/TzMysTr7maWZmVicnTzMzszo5eZqZmdXJydPMzKxOTp5mZmZ1cvI0MzOrk5OnmZlZnZw8zczM6uTkaWZmVicnT7MWkrS/pMgeb6qy/f257TtKGp973dfjeknDJH1N0nWSFkh6TtKfJH1OUum/dUnTJM2sUv4GSd+VdJ+kZZKWSpon6WxJ/1Hn9+HtWdz/1UedEyS9ImkTSetIelLS7rXqm7WSp+czayFJ+wM/AZ4DzoyIYwvbzyFN1j8C2Ilsfs7Cbm4hzVH6/VzZItIk4I8DM0hrNS4mrWpyKHBaRBxRIr7NSEtZbRe59RIlvR+4ElgAfA+4mzSX6FtJ89i+KSLW6W//hWPdTfof9JYq2wQ8BDwaEdtnZYeSJv3+t1ixPqNZe7RykmI//BjqD1ZMmH4uadkx5batAzxLSq4B7FhjH0GVVSBI6xmuV6X8HOB5YJ0S8U0D/lgoW5+0OsXvyVbgKWwX8OUBfC8Oz9pSbfLu7bNtB+TKxgAvAHu0++fohx8+bWvWHucDbwTekyv7BCkB/mIgO4yI5RHxzyqb/khaA3FslW29JK0F7ANcWNj0BVICPTgillY5bkTE96rs7wuS7pL0vKSnJf1Y0nq5KhcAy4F9q4SzH2nVi0tzx+kBrgE+31c7zFrBydOsPf4K3MirE8d+wOWk062NtD1pybMn+6n3TmA06VRx3iTgyYi4s+wBJU0lnd6dDXwMOAL4EHC1pGEAEfEUKRnuXSnL3rs26dT15RGxqLDrG4HtszpmbePkadY+M4DdJa0t6bWkxXlnNPIAknYG9gC+ExEv91P9naRTpXcXyjcmJfvivodJWj33UFY+npQsT4iIwyPi2oj4SRbHNqSFmSvOAzYkLQZe8XHSQs7nVYnxTmBNoK4BSmaN5uRp1j6XkE6nTgY+AzxFGujTEJImAD8DrgdOKfGW1wGLIuLF4q5q1L+XtDhx5TEpK9+J9L/lgnxyBW4lDWx6X24fV5J6xcUe+BOkXmvRP3KxmrXN6u0OwGyoiojnJF1BShzjgQsi4pWsAzcokjYFZpEGJX28RK8TYG3SgJyix4AJVco/SRrktDVwdq58g+zrAzWOs37lSUS8IOnnwH6SRgDDScn39Ih4pcp7l2Vf6xrZa9ZoTp5m7TUD+BWpp7ZXI3YoaWNSD3YR8KEq1w1reYY0orXoOmBHSW/PX/eMiHuz472myn4gnYrtqXGcvPOAL5Guc44m/V+qdfq6MuDo6RrbzVrCydOsvWYBFwMLK8loMCT9CytOd+4UEf/oq37BfGANSRtHxOO58h+SrmFOk/TBaiNuC2YBrwBviIhZ/R00Im6R9BdSD3wMcEcf34tNsq/39bdfs2Zy8jRro4hYTuN6nOuQRq+OBw4ANs56oRVz++mF3ph93YY02UIlxqezmX0uB/4sqTJJwivAONI1ygCWZPUflHQKMF3SvwI3kO4zHUc6JfujiPht4dgzgBNI11cP6SPGbYG/RcRDfdQxazonT7NVx4asmI3ogirbP0AaPFRVRDwi6TbSAKbLCtvmSHor8DXSLD/jSAnzkWyfR0TEn3P1j5Y0L6t7UFb3MdLp5PurHP584HjgZdIgp1o+AlzUx3azlvD0fGbWK5s+8EzgtSVOz7aUpG2Bm0nT8/2l3fHY0ObkaWa9sskK7gHOiYjvtDuePEmXAz0RcUC7YzHzfZ5m1iu7BnsA0Gm9zrVJEyR8o92xmIF7nmZmZnVzz9PMzKxOTp5mZmZ1cvI0MzOrk5OnmZlZnZw8zczM6vT/AT7cJyoZLq0oAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "make_fig(\"MT2 (GeV)\", \"Number of events\")\n",
    "title(\"MT2 for Ml 500 Gev MDM 1 GeV\")\n",
    "nbins   = 50 # number of bins in histogram\n",
    "xlow    = 1\n",
    "xhigh   = 700\n",
    "bins    = np.linspace(xlow, xhigh, nbins)\n",
    "n, bins, patches = plt.hist(np.array(mt2_0)/1e3, bins, alpha = 0.7)\n",
    "axvline(np.sqrt(500**2-1**2), linestyle = \"--\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}