{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b62b5ce8-fbf5-41f3-b109-c05ac94a0b35",
   "metadata": {},
   "source": [
    "# Planck cluster counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2a2825a6-6e82-42e0-97b9-65d34f8af150",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from classy_sz import Class as Class_sz"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "688b5e1c-419d-42d5-a987-c5a862bd24f8",
   "metadata": {},
   "source": [
    "## Binned calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee13b313-29bf-4e7a-b72a-e8139bab9e10",
   "metadata": {},
   "source": [
    "Get accurate binned cluster calculations in less than a second!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "0ee0775b-80a4-446d-ae42-4e2ab4a8deb0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 10.3 s, sys: 3.3 s, total: 13.6 s\n",
      "Wall time: 2.07 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "class_sz = Class_sz()\n",
    "\n",
    "\n",
    "class_sz.set({\n",
    "\n",
    "'output' : 'sz_cluster_counts',\n",
    "'mass_function' :  'T08M500c',\n",
    "'has_selection_function' : 1,\n",
    "'experiment' : 0,\n",
    "'y_m_relation' : 0,\n",
    "'use_skyaveraged_noise': 0,\n",
    "'use_planck_binned_proba' : 0,\n",
    "\n",
    "\n",
    "'m_pivot_ym_[Msun]': 3e14,    \n",
    "\n",
    "'M_min' : 1e13*0.7, \n",
    "'M_max' : 1e16*0.7,\n",
    "\n",
    "'z_min' : 0.01,\n",
    "'z_max' : 1.02,\n",
    "\n",
    "'omega_b': 0.0224178568132,\n",
    "'omega_cdm':  0.11933148326520002,\n",
    "'H0': 70.,\n",
    "'tau_reio': 0.0561,\n",
    "'ln10^{10}A_s': 2.9799585,\n",
    "'n_s': 0.96,\n",
    "\n",
    "\n",
    "\n",
    "'bin_z_min_cluster_counts' : 0.01,\n",
    "'bin_z_max_cluster_counts' : 1.01,\n",
    "'bin_dz_cluster_counts' : 0.1,\n",
    "\n",
    "'bin_dlog10_snr': 0.25, # fiducial 0.25\n",
    "'log10_snr_min' : 0.7, # fiducial 0.7\n",
    "\n",
    "#the paramater dlny is crucial, it controls the speed of the calculation\n",
    "#important to check that lnymin and lnymax are broad enough\n",
    "'dlny' : 0.07, # fiducial 0.05\n",
    "'lnymin' : -11, # fiducial -11\n",
    "'lnymax' : 1., # fiducial 1.\n",
    "'sigmaM_ym' :0.173, # fiducial 0.173\n",
    "\n",
    "\n",
    "'dlnM_cluster_count_completeness_grid' : 0.05, # 0.01 fiducial\n",
    "\n",
    "\n",
    "'cluster_count_completeness_grid_z_cutoff_low' : 0.5,\n",
    "'cluster_count_completeness_grid_z_cutoff_mid' : 1.5,\n",
    "'dz_cluster_count_completeness_grid_low_z' : 5e-3,\n",
    "'dz_cluster_count_completeness_grid_mid_z' : 1e-2,\n",
    "'dz_cluster_count_completeness_grid_high_z' : 1e-1,\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "'signal-to-noise_cut-off_for_survey_cluster_completeness' : 6.,\n",
    "\n",
    "\n",
    "# X ray mass bias (if applicable)\n",
    "'B' : 1.25,\n",
    "\n",
    "'ndim_redshifts' : 100,\n",
    "\n",
    "'szcc_dof': 0.,\n",
    "'cosmo_model': 1,\n",
    "\n",
    "})\n",
    "class_sz.compute_class_szfast()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "b20a83b0-c3ed-4bba-a816-1bb0209347ca",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ntot 5.1269644816e+02\n"
     ]
    }
   ],
   "source": [
    "dNdzdy_theoretical = class_sz.dndzdy_theoretical()['dndzdy']\n",
    "z_center = class_sz.dndzdy_theoretical()['z_center']\n",
    "z_edges = class_sz.dndzdy_theoretical()['z_edges']\n",
    "log10y_center = class_sz.dndzdy_theoretical()['log10y_center']\n",
    "log10y_edges = class_sz.dndzdy_theoretical()['log10y_edges']\n",
    "\n",
    "N_z,N_y = np.shape(dNdzdy_theoretical)\n",
    "N_clusters_z_theory = []\n",
    "\n",
    "for iz in range(N_z):\n",
    "    N_clusters_z_theory.append(np.sum(dNdzdy_theoretical[iz][0:]))\n",
    "N_clusters_y_theory = []\n",
    "\n",
    "for iy in range(N_y):\n",
    "    N_clusters_y_theory.append(np.sum(np.asarray(dNdzdy_theoretical)[:,iy]))\n",
    "print('Ntot %.10e'%(np.sum(N_clusters_z_theory)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "58039759-2291-44ca-890b-fb509af0a4ff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.787079475530407"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "class_sz.sigma8()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "953f84ee-c3d7-45a2-828e-b2a17c0f236f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'counts')"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAGxCAYAAACDV6ltAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0KElEQVR4nO3de3wU9b3/8feGkAQh2Rgg2d0aJCrKVe5gBC1ibKI21UqrWKRAOfAoBT1cROB4IELVIG3Rg+UiHCWeA8ppe8Ry8UQw3JRGgolpG0BUjIIlm2hjdgM09/n94S9b1xANyWYvmdfz8ZhHM9/5zuxnO7b79jvfmbEYhmEIAADAxMICXQAAAECgEYgAAIDpEYgAAIDpEYgAAIDpEYgAAIDpEYgAAIDpEYgAAIDpEYgAAIDphQe6gGDQ0NCgs2fPKjo6WhaLJdDlAACAFjAMQ5WVlXI4HAoLa9sYD4FI0tmzZ5WYmBjoMgAAQCucOXNGV1xxRZuOEdBAdOjQIf3qV79Sfn6+SkpKtH37dt19991efU6cOKFFixbp4MGDqqurU//+/fW///u/6tWrlySpqqpKCxYs0LZt21RdXa3U1FStW7dOCQkJLa4jOjpa0pf/hcbExPjs+wEAgPbjdruVmJjo+R1vi4AGovPnz2vw4MH62c9+pnvuuafJ9lOnTmns2LGaPn26li9frpiYGB07dkxRUVGePvPmzdPu3bv1+9//XlarVXPmzNE999yjw4cPt7iOxstkMTExBCIAAEKML6a7WILl5a4Wi6XJCNHEiRPVuXNn/fd///dF93G5XOrZs6deeukl/ehHP5Ikvffee+rXr59yc3N1ww03tOiz3W63rFarXC4XgQgAgBDhy9/voL3LrKGhQbt379a1116r1NRUxcfHa/To0Xr11Vc9ffLz81VbW6uUlBRPW9++fdWrVy/l5uY2e+zq6mq53W6vBQAAmFfQBqKysjKdO3dOK1euVFpamvbs2aMf/vCHuueee3Tw4EFJktPpVEREhGJjY732TUhIkNPpbPbYmZmZslqtnoUJ1QAAmFvQBqKGhgZJ0l133aV58+ZpyJAhWrx4sb7//e9rw4YNbTr2kiVL5HK5PMuZM2d8UTIAAAhRQXvbfY8ePRQeHq7+/ft7tffr109vvfWWJMlms6mmpkYVFRVeo0SlpaWy2WzNHjsyMlKRkZHtUjcAAAg9QTtCFBERoZEjR+rkyZNe7e+//76uvPJKSdLw4cPVuXNn5eTkeLafPHlSp0+fVnJysl/rBQAAoSugI0Tnzp3Thx9+6FkvLi5WYWGh4uLi1KtXLy1cuFD33Xefbr75Zt1yyy3Kzs7Wzp07deDAAUmS1WrV9OnTNX/+fMXFxSkmJkYPPvigkpOTW3yHGQAAQEBvuz9w4IBuueWWJu1TpkxRVlaWJOmFF15QZmamPv30U1133XVavny57rrrLk/fxgczvvzyy14PZvymS2Zfx233AACEHl/+fgfNc4gCiUAEAEDoMcVziAAAAPyFQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQIRvdKGmTr0X71bvxbt1oaYu0OUAANAuCEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CET4RvUNhufvvOJyr3UAADoKAhGalV1UopTVBz3rUzcf1din9im7qCSAVQEA4HsEIlxUdlGJZm0pUKm72qvd6arSrC0FhCIAQIdCIEIT9Q2Glu88rotdHGtsW77zOJfPAAAdBoEITeQVl6vEVdXsdkNSiatKecXl/isKAIB2RCBCE2WVzYeh1vQDACDYEYjQRHx0lE/7AQAQ7AIaiA4dOqT09HQ5HA5ZLBa9+uqrzfb9+c9/LovFomeeecarvby8XJMmTVJMTIxiY2M1ffp0nTt3rn0L7+BGJcXJbo2SpZntFkl2a5RGJcX5sywAANpNQAPR+fPnNXjwYK1du/Yb+23fvl1vv/22HA5Hk22TJk3SsWPHtHfvXu3atUuHDh3SzJkz26tkU+gUZlFGen9JahKKGtcz0vurU1hzkQkAgNASHsgPv/3223X77bd/Y5+//e1vevDBB/X666/rzjvv9Np24sQJZWdn6+jRoxoxYoQk6dlnn9Udd9yhX//61xcNUGiZtIF2rX9gmDJ2HPO69d5mjVJGen+lDbQHsDoAAHwroIHo2zQ0NGjy5MlauHChBgwY0GR7bm6uYmNjPWFIklJSUhQWFqYjR47ohz/84UWPW11drerqf/7Iu91u3xffAaQNtGvMNT006LE9kqSsaSN1U5+ejAwBADqcoJ5U/dRTTyk8PFwPPfTQRbc7nU7Fx8d7tYWHhysuLk5Op7PZ42ZmZspqtXqWxMREn9YtSRdq6tR78W71XrxbF2rqfH58f/lq+BmVFEcYAgB0SEEbiPLz8/Uf//EfysrKksXi2x/hJUuWyOVyeZYzZ8749PgAACC0BG0gevPNN1VWVqZevXopPDxc4eHh+uSTT7RgwQL17t1bkmSz2VRWVua1X11dncrLy2Wz2Zo9dmRkpGJiYrwWAABgXkE7h2jy5MlKSUnxaktNTdXkyZM1bdo0SVJycrIqKiqUn5+v4cOHS5L27dunhoYGjR492u81AwCA0BTQQHTu3Dl9+OGHnvXi4mIVFhYqLi5OvXr1Uvfu3b36d+7cWTabTdddd50kqV+/fkpLS9OMGTO0YcMG1dbWas6cOZo4cSJ3mAEAgBYL6CWzd955R0OHDtXQoUMlSfPnz9fQoUO1bNmyFh9j69at6tu3r2699VbdcccdGjt2rDZu3NheJQMAgA4ooCNE48aNk2G0/I3pH3/8cZO2uLg4vfTSSz6sCgAAmE3QTqoGAADwFwIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwvYC+7b4jq28wPH/nFZfrpj491SnMEsCKWueyiHB9vPLOQJcBAEC7YoSoHWQXlShl9UHP+tTNRzX2qX3KLioJYFUAAKA5BCIfyy4q0awtBSp1V3u1O11VmrWlgFAEAEAQIhD5UH2DoeU7j8u4yLbGtuU7j3tdTgMAAIFHIPKhvOJylbiqmt1uSCpxVSmvuNx/RQEAgG9FIPKhssrmw1Br+gEAAP8gEPlQfHSUT/sBAAD/IBD50KikONmtUWru5nqLJLs1SqOS4vxZFgAA+BYEIh/qFGZRRnp/SWoSihrXM9L7h+TziAAA6MgIRD6WNtCu9Q8MU3xMpFe7zRql9Q8MU9pAe4AqAwAAzeFJ1e0gbaBdY67poUGP7ZEkZU0bGbJPqgYAwAwYIWonXw0/o5LiCEMAAAQxAhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADA9AhEAADC9gAaiQ4cOKT09XQ6HQxaLRa+++qpnW21trRYtWqRBgwapa9eucjgc+ulPf6qzZ896HaO8vFyTJk1STEyMYmNjNX36dJ07d87P3wQAAISygAai8+fPa/DgwVq7dm2TbRcuXFBBQYGWLl2qgoICvfLKKzp58qR+8IMfePWbNGmSjh07pr1792rXrl06dOiQZs6c6a+vAAAAOoCAvrrj9ttv1+23337RbVarVXv37vVq++1vf6tRo0bp9OnT6tWrl06cOKHs7GwdPXpUI0aMkCQ9++yzuuOOO/TrX/9aDoej3b8DAAAIfSE1h8jlcslisSg2NlaSlJubq9jYWE8YkqSUlBSFhYXpyJEjzR6nurpabrfbawEAAOYVMoGoqqpKixYt0v3336+YmBhJktPpVHx8vFe/8PBwxcXFyel0NnuszMxMWa1Wz5KYmNiutQMAgOAWEoGotrZW9957rwzD0Pr169t8vCVLlsjlcnmWM2fO+KBKAAAQqgI6h6glGsPQJ598on379nlGhyTJZrOprKzMq39dXZ3Ky8tls9maPWZkZKQiIyPbrWYAABBagnqEqDEMffDBB3rjjTfUvXt3r+3JycmqqKhQfn6+p23fvn1qaGjQ6NGj/V0uAAAIUQEdITp37pw+/PBDz3pxcbEKCwsVFxcnu92uH/3oRyooKNCuXbtUX1/vmRcUFxeniIgI9evXT2lpaZoxY4Y2bNig2tpazZkzRxMnTuQOMwAA0GIBDUTvvPOObrnlFs/6/PnzJUlTpkzRY489ph07dkiShgwZ4rXf/v37NW7cOEnS1q1bNWfOHN16660KCwvThAkTtGbNGr/UDwAAOgaLYRhGoIsINLfbLavVKpfL5TVHCQAABC9f/n4H9RwiAAAAfyAQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQwTQu1NSp9+Ld6r14ty7U1AW6HABAECEQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0yMQAQAA0wtoIDp06JDS09PlcDhksVj06quvem03DEPLli2T3W5Xly5dlJKSog8++MCrT3l5uSZNmqSYmBjFxsZq+vTpOnfunB+/BUJFfYPh+TuvuNxrHQBgbgENROfPn9fgwYO1du3ai25ftWqV1qxZow0bNujIkSPq2rWrUlNTVVVV5ekzadIkHTt2THv37tWuXbt06NAhzZw5019fASEiu6hEKasPetanbj6qsU/tU3ZRSQCrAgAEC4thGEHxr8kWi0Xbt2/X3XffLenL0SGHw6EFCxbo4YcfliS5XC4lJCQoKytLEydO1IkTJ9S/f38dPXpUI0aMkCRlZ2frjjvu0KeffiqHw9Giz3a73bJarXK5XIqJiWmX74fAyS4q0awtBfr6P+iW//+f6x8YprSBdn+XBQBoI1/+fgftHKLi4mI5nU6lpKR42qxWq0aPHq3c3FxJUm5urmJjYz1hSJJSUlIUFhamI0eONHvs6upqud1urwUdU32DoeU7jzcJQ5I8bct3HufyGQCYXNAGIqfTKUlKSEjwak9ISPBsczqdio+P99oeHh6uuLg4T5+LyczMlNVq9SyJiYk+rh7BIq+4XCWuqma3G5JKXFXKKy73X1EAgKATtIGoPS1ZskQul8uznDlzJtAloZ2UVTYfhlrTDwDQMQVtILLZbJKk0tJSr/bS0lLPNpvNprKyMq/tdXV1Ki8v9/S5mMjISMXExHgt6Jjio6N82g8A0DEFbSBKSkqSzWZTTk6Op83tduvIkSNKTk6WJCUnJ6uiokL5+fmePvv27VNDQ4NGjx7t95oRfEYlxclujfJMoP46iyS7NUqjkuL8WRYAIMgENBCdO3dOhYWFKiwslPTlROrCwkKdPn1aFotFc+fO1eOPP64dO3bor3/9q37605/K4XB47kTr16+f0tLSNGPGDOXl5enw4cOaM2eOJk6c2OI7zNCxdQqzKCO9vyQ1CUWN6xnp/dUprLnIBAAwg4Dedn/gwAHdcsstTdqnTJmirKwsGYahjIwMbdy4URUVFRo7dqzWrVuna6+91tO3vLxcc+bM0c6dOxUWFqYJEyZozZo16tatW4vr4Lb7ji+7qEQZO46p1F3tabNbo5SR3p9b7gEgRPny9ztonkMUSAQic6isqtWgx/ZIkrKmjdRNfXoyMgQAIcwUzyECfO2r4WdUUhxhCADgQSACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmRyACAACmFx7oAgB/uSwiXB+vvDPQZQAAghAjRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPQIRAAAwPRaFYjOnDmjTz/91LOel5enuXPnauPGjT4rDAAAwF9aFYh+8pOfaP/+/ZIkp9Op2267TXl5eXr00Ue1YsUKnxYIAADQ3loViIqKijRq1ChJ0u9+9zsNHDhQf/rTn7R161ZlZWX5sj4AAIB216pAVFtbq8jISEnSG2+8oR/84AeSpL59+6qkpMR31QEAAPhBqwLRgAEDtGHDBr355pvau3ev0tLSJElnz55V9+7dfVogAABAe2tVIHrqqaf03HPPady4cbr//vs1ePBgSdKOHTs8l9IAAABChcUwDKM1O9bX18vtduvyyy/3tH388cfq2rWrevbs6bMC/cHtdstqtcrlcikmJibQ5QAAgBbw5e93q0aIxo8fr8rKSq8wJElxcXG677772lQQAACAv7UqEB04cEA1NTVN2quqqvTmm2+2uSgAAAB/Cr+Uzn/5y188fx8/flxOp9OzXl9fr+zsbH3nO9/xXXUAAAB+cEmBaMiQIbJYLLJYLBo/fnyT7V26dNGzzz7rs+Lq6+v12GOPacuWLXI6nXI4HJo6dar+/d//XRaLRZJkGIYyMjK0adMmVVRUaMyYMVq/fr369OnjszoAAEDHdkmBqLi4WIZh6KqrrlJeXp7X5OmIiAjFx8erU6dOPivuqaee0vr16/Xiiy9qwIABeueddzRt2jRZrVY99NBDkqRVq1ZpzZo1evHFF5WUlKSlS5cqNTVVx48fV1RUlM9qAQAAHVer7zLzh+9///tKSEjQ888/72mbMGGCunTpoi1btsgwDDkcDi1YsEAPP/ywJMnlcikhIUFZWVmaOHFiiz6Hu8wAAAg9vvz9vqQRoq/64IMPtH//fpWVlamhocFr27Jly9pUVKMbb7xRGzdu1Pvvv69rr71Wf/7zn/XWW29p9erVkr4csXI6nUpJSfHsY7VaNXr0aOXm5jYbiKqrq1VdXe1Zd7vdPqkXAACEplYFok2bNmnWrFnq0aOHbDabZz6PJFksFp8FosWLF8vtdqtv377q1KmT6uvr9cQTT2jSpEmS5JnUnZCQ4LVfQkKC14Tvr8vMzNTy5ct9UiMAAAh9rQpEjz/+uJ544gktWrTI1/V4+d3vfqetW7fqpZde0oABA1RYWKi5c+fK4XBoypQprT7ukiVLNH/+fM+62+1WYmKiL0oGAAAhqFWB6IsvvtCPf/xjX9fSxMKFC7V48WLPpa9Bgwbpk08+UWZmpqZMmSKbzSZJKi0tld1u9+xXWlqqIUOGNHvcyMhIz8tpAQAAWvVgxh//+Mfas2ePr2tp4sKFCwoL8y6xU6dOnjlLSUlJstlsysnJ8Wx3u906cuSIkpOT270+AADQMbRqhOiaa67R0qVL9fbbb2vQoEHq3Lmz1/bGW+LbKj09XU888YR69eqlAQMG6N1339Xq1av1s5/9TNKX85Xmzp2rxx9/XH369PHcdu9wOHT33Xf7pAYAANDxteq2+6SkpOYPaLHoo48+alNRjSorK7V06VJt375dZWVlcjgcuv/++7Vs2TJFRERI+ueDGTdu3KiKigqNHTtW69at07XXXtviz+G2ewAAQo8vf7+D+jlE/kIgAgAg9AT8bfcAAAAdSavmEDXO4WnOCy+80KpiAAAAAqHVt91/VW1trYqKilRRUXHRl74CAAAEs1YFou3btzdpa2ho0KxZs3T11Ve3uSgAAAB/8tkcorCwMM2fP19PP/20rw4JAADgFz6dVH3q1CnV1dX58pAAAADtrlWXzL76HjDpy2cBlZSUaPfu3W16xxgAAEAgtCoQvfvuu17rYWFh6tmzp37zm9986x1oAAAAwaZVgWj//v2+rgMAACBgWhWIGn322Wc6efKkJOm6665Tz549fVIUAACAP7VqUvX58+f1s5/9THa7XTfffLNuvvlmORwOTZ8+XRcuXPB1jQAAAO2qVYFo/vz5OnjwoHbu3KmKigpVVFToj3/8ow4ePKgFCxb4ukYAAIB21aqXu/bo0UN/+MMfNG7cOK/2/fv3695779Vnn33mq/r8gpe7AgAQegL+ctcLFy4oISGhSXt8fDyXzAAAQMhpVSBKTk5WRkaGqqqqPG3/+Mc/tHz5ciUnJ/usOAAAAH9o1V1mzzzzjNLS0nTFFVdo8ODBkqQ///nPioyM1J49e3xaIABvF2rq1H/Z65Kk4ytSdVlEm24WBQColYFo0KBB+uCDD7R161a99957kqT7779fkyZNUpcuXXxaIAAAQHtrVSDKzMxUQkKCZsyY4dX+wgsv6LPPPtOiRYt8UhwAAIA/tGoO0XPPPae+ffs2aR8wYIA2bNjQ5qIAAAD8qVWByOl0ym63N2nv2bOnSkpK2lwUAACAP7UqECUmJurw4cNN2g8fPiyHw9HmogAAAPypVXOIZsyYoblz56q2tlbjx4+XJOXk5OiRRx7hSdUAACDktCoQLVy4UH//+9/1i1/8QjU1NZKkqKgoLVq0SEuWLPFpgQAAAO2tVYHIYrHoqaee0tKlS3XixAl16dJFffr0UWRkpK/rAwAAaHdteqJbt27dNHLkSF/VAgAAEBCtmlQNAADQkRCIAACA6RGIgBBT32B4/s4rLvdaBwC0DoEICCHZRSVKWX3Qsz5181GNfWqfsot4ICoAtAWBCAgR2UUlmrWlQKXuaq92p6tKs7YUEIoAoA0IREAIqG8wtHzncV3s4lhj2/Kdx7l8BgCtRCACQkBecblKXFXNbjcklbiqlFdc7r+iAKADIRABIaCssvkw1Jp+AABvQR+I/va3v+mBBx5Q9+7d1aVLFw0aNEjvvPOOZ7thGFq2bJnsdru6dOmilJQUffDBBwGsGPC9+Ogon/YDAHgL6kD0xRdfaMyYMercubP+7//+T8ePH9dvfvMbXX755Z4+q1at0po1a7RhwwYdOXJEXbt2VWpqqqqq+DdldByjkuJkt0bJ0sx2iyS7NUqjkuL8WRYAdBgWwzCCdhbm4sWLdfjwYb355psX3W4YhhwOhxYsWKCHH35YkuRyuZSQkKCsrCxNnDixRZ/jdrtltVrlcrkUExPjs/oBX2q8y0yS1+TqxpC0/oFhShto93tdABAovvz9DuoRoh07dmjEiBH68Y9/rPj4eA0dOlSbNm3ybC8uLpbT6VRKSoqnzWq1avTo0crNzW32uNXV1XK73V4LEOzSBtq1/oFhio/xfomyzRpFGAKANgrqQPTRRx9p/fr16tOnj15//XXNmjVLDz30kF588UVJktPplCQlJCR47ZeQkODZdjGZmZmyWq2eJTExsf2+BOBDaQPtemP+dz3rWdNG6q1F4wlDANBGQR2IGhoaNGzYMD355JMaOnSoZs6cqRkzZmjDhg1tOu6SJUvkcrk8y5kzZ3xUMdD+OoX9cybRqKQ4r3UAQOsEdSCy2+3q37+/V1u/fv10+vRpSZLNZpMklZaWevUpLS31bLuYyMhIxcTEeC0AAMC8gjoQjRkzRidPnvRqe//993XllVdKkpKSkmSz2ZSTk+PZ7na7deTIESUnJ/u1VgAAELrCA13AN5k3b55uvPFGPfnkk7r33nuVl5enjRs3auPGjZIki8WiuXPn6vHHH1efPn2UlJSkpUuXyuFw6O677w5s8QAAIGQEdSAaOXKktm/friVLlmjFihVKSkrSM888o0mTJnn6PPLIIzp//rxmzpypiooKjR07VtnZ2YqK4gF1AACgZYL6OUT+wnOIEEou1NSp/7LXJUnHV6Tqsoig/vcaAGg3pnkOEQAAgD8QiAAAgOkRiAAAgOkRiAAAgOkRiAAAgOlxewoQYi6LCNfHK+8MdBkA0KEwQgQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQAAEyPQAQgIC7U1Kn34t3qvXi3LtTUBbocACZHIAIAAKYXUoFo5cqVslgsmjt3rqetqqpKs2fPVvfu3dWtWzdNmDBBpaWlgSsSAACEnJAJREePHtVzzz2n66+/3qt93rx52rlzp37/+9/r4MGDOnv2rO65554AVQkAAEJRSASic+fOadKkSdq0aZMuv/xyT7vL5dLzzz+v1atXa/z48Ro+fLg2b96sP/3pT3r77bcDWDEAAAglIRGIZs+erTvvvFMpKSle7fn5+aqtrfVq79u3r3r16qXc3Nxmj1ddXS232+21AAAA8woPdAHfZtu2bSooKNDRo0ebbHM6nYqIiFBsbKxXe0JCgpxOZ7PHzMzM1PLly31dKgAACFFBPUJ05swZ/eu//qu2bt2qqKgonx13yZIlcrlcnuXMmTM+OzYAAAg9QR2I8vPzVVZWpmHDhik8PFzh4eE6ePCg1qxZo/DwcCUkJKimpkYVFRVe+5WWlspmszV73MjISMXExHgtAPyrvsHw/J1XXO61DgD+FtSXzG699Vb99a9/9WqbNm2a+vbtq0WLFikxMVGdO3dWTk6OJkyYIEk6efKkTp8+reTk5ECUDKAFsotKlLHjmGd96uajslujlJHeX2kD7QGsDIBZBXUgio6O1sCBA73aunbtqu7du3vap0+frvnz5ysuLk4xMTF68MEHlZycrBtuuCEQJQP4FtlFJZq1pUBfHw9yuqo0a0uB1j8wjFAEwO+COhC1xNNPP62wsDBNmDBB1dXVSk1N1bp16wJdFoCLqG8wtHzn8SZhSJIMSRZJy3ce1239beoUZvFzdQDMzGIYhukv3LvdblmtVrlcLuYTAe0o99Tfdf+mb39G2MszblDy1d39UBGAUObL3++gnlQNoGMpq6zyaT8A8BUCEQC/iY9u2eMzWtoPAHyFQATAb0YlxclujVJzs4MskuzWKI1KivNnWQBAIALgP53CLMpI7y9JTUJR43pGen8mVAPwOwIRAL9KG2jX+geGKT4m0qvdZo3ilnsAARPyt90DCD1pA+0ac00PDXpsjyQpa9pI3dSnJyNDAAKGESIAAfHV8DMqKY4wBCCgCEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0eA4RgIC4LCJcH6+8M9BlAIAkRogAAAAIRAAAAAQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAABgegQiAGiDCzV16r14t3ov3q0LNXWBLgdAKxGIAACA6RGIAACA6RGIAACA6RGIAACA6RGIAKAN6hsMz995xeVe6wBCB4EIAFopu6hEKasPetanbj6qsU/tU3ZRSQCrAtAaBCIAaIXsohLN2lKgUne1V7vTVaVZWwoIRUCIIRABwCWqbzC0fOdxXeziWGPb8p3HuXwGhJCgD0SZmZkaOXKkoqOjFR8fr7vvvlsnT5706lNVVaXZs2ere/fu6tatmyZMmKDS0tIAVQygo8srLleJq6rZ7YakEleV8orL/VcUgDYJ+kB08OBBzZ49W2+//bb27t2r2tpafe9739P58+c9febNm6edO3fq97//vQ4ePKizZ8/qnnvuCWDVADqyssrmw1Br+gEIvPBAF/BtsrOzvdazsrIUHx+v/Px83XzzzXK5XHr++ef10ksvafz48ZKkzZs3q1+/fnr77bd1ww03BKJsAB1YfHSUT/sBCLygHyH6OpfLJUmKi4uTJOXn56u2tlYpKSmePn379lWvXr2Um5t70WNUV1fL7XZ7LQDQUqOS4mS3RsnSzHaLJLs1SqOS4vxZFoA2CKlA1NDQoLlz52rMmDEaOHCgJMnpdCoiIkKxsbFefRMSEuR0Oi96nMzMTFmtVs+SmJjY3qUD6EA6hVmUkd5fkpqEosb1jPT+6hTWXGQCEGxCKhDNnj1bRUVF2rZtW5uOs2TJErlcLs9y5swZH1UIwCzSBtq1/oFhio+J9Gq3WaO0/oFhShtoD1BlAFoj6OcQNZozZ4527dqlQ4cO6YorrvC022w21dTUqKKiwmuUqLS0VDab7aLHioyMVGRk5EW3AUBLpQ20a8w1PTTosT2SpKxpI3VTn56MDAEhKOhHiAzD0Jw5c7R9+3bt27dPSUlJXtuHDx+uzp07Kycnx9N28uRJnT59WsnJyf4uF4DJfDX8jEqKIwwBISroR4hmz56tl156SX/84x8VHR3tmRdktVrVpUsXWa1WTZ8+XfPnz1dcXJxiYmL04IMPKjk5mTvMAABAiwR9IFq/fr0kady4cV7tmzdv1tSpUyVJTz/9tMLCwjRhwgRVV1crNTVV69at83OlAAAgVAV9IDKMb3/0fVRUlNauXau1a9f6oSIAANDRBH0gAoBgdllEuD5eeWegywDQRkE/qRoAAKC9EYgAAIDpEYgAAIDpEYgAALpQU6fei3er9+LdulBTF+hyAL8jEAEAANMjEAEAVN/wz0ec5BWXe60DZkAgAgCTyy4qUcrqg571qZuPauxT+5RdVBLAqgD/IhABgIllF5Vo1pYClbqrvdqdrirN2lJAKIJpEIgAwKTqGwwt33lcF7s41ti2fOdxLp/BFAhEAGBSecXlKnFVNbvdkFTiqlJecbn/igIChEAEACZVVtl8GGpNPyCUEYgAwKTio6N82g8IZQQiADCpUUlxslujZGlmu0WS3RqlUUlx/iwLCAgCEQCYVKcwizLS+0tSk1DUuJ6R3l+dwpqLTEDHQSACABNLG2jX+geGKT4m0qvdZo3S+geGKW2gPUCVAf4VHugCAACBlTbQrjHX9NCgx/ZIkrKmjdRNfXoyMgRTYYQIAOAVfkYlxRGGYDqMEAEAdFlEuD5eeWegy2izyqpaRrrQKowQAQA6BN7JhrYgEAEAQh7vZENbEYgAACGNd7LBFwhEAICQxjvZ4AsEIgBASOOdbPAFAhEAIKTxTjb4ArfdAwBCWuM72ZyuqovOI7Loyydvh9I72eobDOUVl6usskrx0VE8G8oPCEQAgJDW+E62WVsKZJG8QlEovpMtu6hEGTuOed0xZ7dGKSO9P69SaUdcMgMAhLyO8k42Hh8QOIwQAQA6hLSBdt3W3xayl5q+7fEBFn35+IDb+ttC6juFyvkgEAEAOoxOYRYlX9090GW0yqU8PiAUvmN2UYmW7zzu9Z2C+dIfl8wAAAgCHenxAY2X/r4e8IL50h+BCACAINBRHh8Qqk8O7zCBaO3aterdu7eioqI0evRo5eXlBbokAABarPHxAc3NsLHoy0tOwf74gFB9cniHCET/8z//o/nz5ysjI0MFBQUaPHiwUlNTVVZWFujSAABokcbHB0hqEopC6fEBoXrpr0MEotWrV2vGjBmaNm2a+vfvrw0bNuiyyy7TCy+8EOjSAABoscbHB9is3pfFQunxAaF66S/k7zKrqalRfn6+lixZ4mkLCwtTSkqKcnNzL7pPdXW1qqv/+YwHt9vd7nUCANASof74gFB9cnjIjxB9/vnnqq+vV0JCgld7QkKCnE7nRffJzMyU1Wr1LImJif4oFQCAFml8fMBdQ76j5Ku7h0wYkkL30l/IB6LWWLJkiVwul2c5c+ZMoEsCAKDDCMVLfyF/yaxHjx7q1KmTSktLvdpLS0tls9kuuk9kZKQiIyMvug0AALRdqF36C/kRooiICA0fPlw5OTmetoaGBuXk5Cg5OTmAlQEAYG6hdOkv5EeIJGn+/PmaMmWKRowYoVGjRumZZ57R+fPnNW3atECXBgAAQkCHCET33XefPvvsMy1btkxOp1NDhgxRdnZ2k4nWAAAAF2MxDCO4np0dAG63W1arVS6XSzExMYEuBwAAtIAvf79Dfg4RAABAWxGIAACA6RGIAACA6RGIAACA6RGIAACA6RGIAACA6RGIAACA6RGIAACA6RGIAACA6XWIV3e0VePDut1ud4ArAQAALdX4u+2Ll24QiCRVVlZKkhITEwNcCQAAuFSVlZWyWq1tOgbvMpPU0NCgs2fPKjo6WhaLJdDldHhut1uJiYk6c+YM744LEpyT4ML5CC6cj+Dy1fMRHR2tyspKORwOhYW1bRYQI0SSwsLCdMUVVwS6DNOJiYnh/1yCDOckuHA+ggvnI7g0no+2jgw1YlI1AAAwPQIRAAAwPQIR/C4yMlIZGRmKjIwMdCn4/zgnwYXzEVw4H8Glvc4Hk6oBAIDpMUIEAABMj0AEAABMj0AEAABMj0AEAABMj0CEdrF27Vr17t1bUVFRGj16tPLy8prtu2nTJt100026/PLLdfnllyslJeUb+6N1LuWcfNW2bdtksVh09913t2+BJnOp56OiokKzZ8+W3W5XZGSkrr32Wr322mt+qrbju9Tz8cwzz+i6665Tly5dlJiYqHnz5qmqqspP1XZshw4dUnp6uhwOhywWi1599dVv3efAgQMaNmyYIiMjdc011ygrK+vSP9gAfGzbtm1GRESE8cILLxjHjh0zZsyYYcTGxhqlpaUX7f+Tn/zEWLt2rfHuu+8aJ06cMKZOnWpYrVbj008/9XPlHdelnpNGxcXFxne+8x3jpptuMu666y7/FGsCl3o+qqurjREjRhh33HGH8dZbbxnFxcXGgQMHjMLCQj9X3jFd6vnYunWrERkZaWzdutUoLi42Xn/9dcNutxvz5s3zc+Ud02uvvWY8+uijxiuvvGJIMrZv3/6N/T/66CPjsssuM+bPn28cP37cePbZZ41OnToZ2dnZl/S5BCL43KhRo4zZs2d71uvr6w2Hw2FkZma2aP+6ujojOjraePHFF9urRNNpzTmpq6szbrzxRuM///M/jSlTphCIfOhSz8f69euNq666yqipqfFXiaZyqedj9uzZxvjx473a5s+fb4wZM6Zd6zSjlgSiRx55xBgwYIBX23333WekpqZe0mdxyQw+VVNTo/z8fKWkpHjawsLClJKSotzc3BYd48KFC6qtrVVcXFx7lWkqrT0nK1asUHx8vKZPn+6PMk2jNedjx44dSk5O1uzZs5WQkKCBAwfqySefVH19vb/K7rBacz5uvPFG5efney6rffTRR3rttdd0xx13+KVmeMvNzfU6f5KUmpra4t+cRrzcFT71+eefq76+XgkJCV7tCQkJeu+991p0jEWLFsnhcDT5Bxyt05pz8tZbb+n5559XYWGhHyo0l9acj48++kj79u3TpEmT9Nprr+nDDz/UL37xC9XW1iojI8MfZXdYrTkfP/nJT/T5559r7NixMgxDdXV1+vnPf65/+7d/80fJ+Bqn03nR8+d2u/WPf/xDXbp0adFxGCFCUFm5cqW2bdum7du3KyoqKtDlmFJlZaUmT56sTZs2qUePHoEuB5IaGhoUHx+vjRs3avjw4brvvvv06KOPasOGDYEuzZQOHDigJ598UuvWrVNBQYFeeeUV7d69W7/85S8DXRragBEi+FSPHj3UqVMnlZaWerWXlpbKZrN9476//vWvtXLlSr3xxhu6/vrr27NMU7nUc3Lq1Cl9/PHHSk9P97Q1NDRIksLDw3Xy5EldffXV7Vt0B9aa/43Y7XZ17txZnTp18rT169dPTqdTNTU1ioiIaNeaO7LWnI+lS5dq8uTJ+pd/+RdJ0qBBg3T+/HnNnDlTjz76qMLCGGvwJ5vNdtHzFxMT0+LRIYkRIvhYRESEhg8frpycHE9bQ0ODcnJylJyc3Ox+q1at0i9/+UtlZ2drxIgR/ijVNC71nPTt21d//etfVVhY6Fl+8IMf6JZbblFhYaESExP9WX6H05r/jYwZM0YffvihJ5hK0vvvvy+73U4YaqPWnI8LFy40CT2NYdXg9aB+l5yc7HX+JGnv3r3f+JtzUZc23xv4dtu2bTMiIyONrKws4/jx48bMmTON2NhYw+l0GoZhGJMnTzYWL17s6b9y5UojIiLC+MMf/mCUlJR4lsrKykB9hQ7nUs/J13GXmW9d6vk4ffq0ER0dbcyZM8c4efKksWvXLiM+Pt54/PHHA/UVOpRLPR8ZGRlGdHS08fLLLxsfffSRsWfPHuPqq6827r333kB9hQ6lsrLSePfdd413333XkGSsXr3aePfdd41PPvnEMAzDWLx4sTF58mRP/8bb7hcuXGicOHHCWLt2LbfdI3g8++yzRq9evYyIiAhj1KhRxttvv+3Z9t3vfteYMmWKZ/3KK680JDVZMjIy/F94B3Yp5+TrCES+d6nn409/+pMxevRoIzIy0rjqqquMJ554wqirq/Nz1R3XpZyP2tpa47HHHjOuvvpqIyoqykhMTDR+8YtfGF988YX/C++A9u/ff9HfhMZzMGXKFOO73/1uk32GDBliREREGFdddZWxefPmS/5ci2EwvgcAAMyNOUQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAAMD0CEQAOoRx48Zp7ty5PutrsVj06quvetbfe+893XDDDYqKitKQIUNaXSeA4EQgAoCLKCkp0e233+5Zz8jIUNeuXXXy5Enl5OQoKytLsbGxgSsQgE/xtnsAQSVY3t7+9Tednzp1SnfeeaeuvPLKAFUEoD0xQgQgoMaNG6c5c+Zo7ty56tGjh1JTU1VUVKTbb79d3bp1U0JCgiZPnqzPP//cs8/58+f105/+VN26dZPdbtdvfvObJsddt26d+vTpo6ioKCUkJOhHP/qR1/aGhgY98sgjiouLk81m02OPPea1/auXzCwWi/Lz87VixQpZLBaNGzdO06ZNk8vlksVikcViabI/gNBCIAIQcC+++KIiIiJ0+PBhrVy5UuPHj9fQoUP1zjvvKDs7W6Wlpbr33ns9/RcuXKiDBw/qj3/8o/bs2aMDBw6ooKDAs/2dd97RQw89pBUrVujkyZPKzs7WzTff3OQzu3btqiNHjmjVqlVasWKF9u7de9H6SkpKNGDAAC1YsEAlJSXasWOHnnnmGcXExKikpEQlJSV6+OGH2+e/HAB+wSUzAAHXp08frVq1SpL0+OOPa+jQoXryySc921944QUlJibq/fffl8Ph0PPPP68tW7bo1ltvlfRluLniiis8/U+fPq2uXbvq+9//vqKjo3XllVdq6NChXp95/fXXKyMjw/P5v/3tb5WTk6PbbrutSX02m03h4eHq1q2b51Ka1WqVxWJpcmkNQGgiEAEIuOHDh3v+/vOf/6z9+/erW7duTfqdOnVK//jHP1RTU6PRo0d72uPi4nTdddd51m+77TZdeeWVuuqqq5SWlqa0tDT98Ic/1GWXXebpc/3113sd2263q6yszJdfC0AI4ZIZgIDr2rWr5+9z584pPT1dhYWFXssHH3zQ5LJXc6Kjo1VQUKCXX35Zdrtdy5Yt0+DBg1VRUeHp07lzZ699LBaLGhoafPJ9AIQeAhGAoDJs2DAdO3ZMvXv31jXXXOO1dO3aVVdffbU6d+6sI0eOePb54osv9P7773sdJzw8XCkpKVq1apX+8pe/6OOPP9a+fft8VmdERITq6+t9djwAgUUgAhBUZs+erfLyct1///06evSoTp06pddff13Tpk1TfX29unXrpunTp2vhwoXat2+fioqKNHXqVIWF/fP/znbt2qU1a9aosLBQn3zyif7rv/5LDQ0NXpfV2qp37946d+6ccnJy9Pnnn+vChQs+OzYA/yMQAQgqDodDhw8fVn19vb73ve9p0KBBmjt3rmJjYz2h51e/+pVuuukmpaenKyUlRWPHjvWahxQbG6tXXnlF48ePV79+/bRhwwa9/PLLGjBggM/qvPHGG/Xzn/9c9913n3r27OmZFA4gNFkMwzACXQQAAEAgMUIEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABMj0AEAABM7/8BissDOrgKpfcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.errorbar(z_center,N_clusters_z_theory,yerr = np.sqrt(N_clusters_z_theory),ls='None',marker='o')\n",
    "plt.xlabel('redshift')\n",
    "plt.ylabel('counts')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "645abe64-5110-4585-be99-4643a90109c9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'counts')"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAGwCAYAAABGogSnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAyVUlEQVR4nO3deVxUZf//8fcgsrjMGAoMJBpqLriWilLebqGgZplWpmZu2a1hZWgud+bWgpW32mbUt9Lu353d1V1WWmqGu6KWSWkaqdmtJoOkyYAGCszvj77Ot7nBlGmGAc/r+XicR5zrXHPmc43mvDnnOueYHA6HQwAAAAbh5+sCAAAAKhLhBwAAGArhBwAAGArhBwAAGArhBwAAGArhBwAAGArhBwAAGIq/rwuoDEpKSnT8+HHVrl1bJpPJ1+UAAIDL4HA4lJeXp8jISPn5Xf7xHMKPpOPHjysqKsrXZQAAADccPXpU9evXv+z+hB9JtWvXlvTbh2c2m31cDQAAuBx2u11RUVHO7/HLRfiRnKe6zGYz4QcAgCqmvFNWmPAMAAAMhfADAAAMhfADAAAMhfADAAAMhfADAAAMhfADAAAMhfADAAAMhfADAAAMhfADAAAMhfADAAAMhfADAAAMhfADAAAMhfADAAAMhfADAAAMhfDjJWfPFemaaZ/ommmf6Oy5Il+XAwAA/hfhBwAAGArhBwAAGArhBwAAGArhBwAAGArhBwAAGArhBwAAGArhBwAAGArhBwAAGArhBwAAGArhBwAAGArhBwAAGArhBwAAGArhx0uKSxzOn3cePuWyDgAAfIfw4wWr92YpfsFG5/rIJV+oy9PrtHpvlg+rAgAAEuHH41bvzdL4f36lbHuhS7stt0Dj//kVAQgAAB8j/HhQcYlDc1bsU1knuC60zVmxj1NgAAD4EOHHg3YePqWs3IKLbndIysot0M7DpyquKAAA4ILw40En8i4efNzpBwAAPI/w40FhtYM82g8AAHge4ceDYqNDFGEJkuki202SIixBio0OqciyAADA7xB+PKian0mz+sdIUqkAdGF9Vv8YVfO7WDwCAADeRvjxsMRWEXr57usVZg50abdagvTy3dcrsVWEjyoDAACS5O/rAq5Eia0idGOTemo9+zNJ0tJRHfWXa0M54gMAQCXAkR8v+X3QiY0OIfgAAFBJEH4AAIChEH4AAIChEH4AAIChEH4AAICh+DT8vPzyy2rTpo3MZrPMZrPi4uK0atUq5/aCggIlJSWpbt26qlWrlgYNGqTs7GyXfRw5ckT9+vVTjRo1FBYWpkceeURFRUUVPRQAAFBF+DT81K9fX/PmzdOuXbv05ZdfqmfPnrr11lv17bffSpIefvhhrVixQu+99542btyo48ePa+DAgc7XFxcXq1+/fjp37py2bdumN998U0uXLtXMmTN9NSQAAFDJmRwOh8PXRfxeSEiInn32Wd1+++0KDQ3VsmXLdPvtt0uSvvvuO7Vo0ULp6enq3LmzVq1apZtvvlnHjx9XeHi4JCk1NVVTp05VTk6OAgICynyPwsJCFRYWOtftdruioqKUm5srs9nskXGcPVekmJlrJEn75iaoRgC3VAIAwJPsdrssFku5v78rzZyf4uJi/etf/9KZM2cUFxenXbt26fz584qPj3f2ad68uRo0aKD09HRJUnp6ulq3bu0MPpKUkJAgu93uPHpUlpSUFFksFucSFRXlvYEBAIBKxefhZ8+ePapVq5YCAwM1btw4LV++XDExMbLZbAoICFCdOnVc+oeHh8tms0mSbDabS/C5sP3CtouZPn26cnNzncvRo0c9OygAAFBp+fxcTLNmzZSRkaHc3Fz9+9//1ogRI7Rx40avvmdgYKACAwMv3REAAFxxfB5+AgIC1KRJE0lS+/bt9cUXX+i5557T4MGDde7cOZ0+fdrl6E92drasVqskyWq1aufOnS77u3A12IU+AAAAv+fz017/raSkRIWFhWrfvr2qV6+utLQ057bMzEwdOXJEcXFxkqS4uDjt2bNHJ06ccPZZu3atzGazYmJiKrz236sR4K8f5/XTj/P6MdkZAIBKxKffytOnT1efPn3UoEED5eXladmyZdqwYYPWrFkji8WiMWPGKDk5WSEhITKbzXrggQcUFxenzp07S5J69+6tmJgYDR8+XM8884xsNptmzJihpKQkTmsBAIAy+TT8nDhxQvfcc4+ysrJksVjUpk0brVmzRr169ZIkLVy4UH5+fho0aJAKCwuVkJCgxYsXO19frVo1rVy5UuPHj1dcXJxq1qypESNGaO7cub4aEgAAqOQq3X1+fMHd+wQAAADfqfL3+QEAAKgIhB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAohB8AAGAoPg0/KSkp6tixo2rXrq2wsDANGDBAmZmZLn26d+8uk8nksowbN86lz5EjR9SvXz/VqFFDYWFheuSRR1RUVFSRQwEAAFWEvy/ffOPGjUpKSlLHjh1VVFSkv/3tb+rdu7f27dunmjVrOvuNHTtWc+fOda7XqFHD+XNxcbH69esnq9Wqbdu2KSsrS/fcc4+qV6+up556qkLHAwAAKj+Tw+Fw+LqIC3JychQWFqaNGzeqa9eukn478tOuXTstWrSozNesWrVKN998s44fP67w8HBJUmpqqqZOnaqcnBwFBARc8n3tdrssFotyc3NlNps9Nh4AAOA97n5/V6o5P7m5uZKkkJAQl/a33npL9erVU6tWrTR9+nSdPXvWuS09PV2tW7d2Bh9JSkhIkN1u17ffflvm+xQWFsput7ssAADAGHx62uv3SkpKNHHiRN14441q1aqVs33o0KFq2LChIiMj9c0332jq1KnKzMzUBx98IEmy2WwuwUeSc91ms5X5XikpKZozZ46XRgIAACqzShN+kpKStHfvXm3ZssWl/b777nP+3Lp1a0VEROimm27SoUOH1LhxY7fea/r06UpOTnau2+12RUVFuVc4AACoUirFaa8JEyZo5cqVWr9+verXr/+HfTt16iRJOnjwoCTJarUqOzvbpc+FdavVWuY+AgMDZTabXRYAAGAMPg0/DodDEyZM0PLly7Vu3TpFR0df8jUZGRmSpIiICElSXFyc9uzZoxMnTjj7rF27VmazWTExMV6pGwAAVF0+Pe2VlJSkZcuW6aOPPlLt2rWdc3QsFouCg4N16NAhLVu2TH379lXdunX1zTff6OGHH1bXrl3Vpk0bSVLv3r0VExOj4cOH65lnnpHNZtOMGTOUlJSkwMBAXw4PAABUQj691N1kMpXZvmTJEo0cOVJHjx7V3Xffrb179+rMmTOKiorSbbfdphkzZricqvrPf/6j8ePHa8OGDapZs6ZGjBihefPmyd//8rIdl7oDAFD1uPv9Xanu8+MrhB8AAKqeK+I+PwAAAN5G+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIZC+AEAAIbi0/CTkpKijh07qnbt2goLC9OAAQOUmZnp0qegoEBJSUmqW7euatWqpUGDBik7O9ulz5EjR9SvXz/VqFFDYWFheuSRR1RUVFSRQwEAAFWET8PPxo0blZSUpO3bt2vt2rU6f/68evfurTNnzjj7PPzww1qxYoXee+89bdy4UcePH9fAgQOd24uLi9WvXz+dO3dO27Zt05tvvqmlS5dq5syZvhgSAACo5EwOh8Ph6yIuyMnJUVhYmDZu3KiuXbsqNzdXoaGhWrZsmW6//XZJ0nfffacWLVooPT1dnTt31qpVq3TzzTfr+PHjCg8PlySlpqZq6tSpysnJUUBAQKn3KSwsVGFhoXPdbrcrKipKubm5MpvNFTNYAADwp9jtdlkslnJ/f1eqOT+5ubmSpJCQEEnSrl27dP78ecXHxzv7NG/eXA0aNFB6erokKT09Xa1bt3YGH0lKSEiQ3W7Xt99+W+b7pKSkyGKxOJeoqChvDQkAAFQylSb8lJSUaOLEibrxxhvVqlUrSZLNZlNAQIDq1Knj0jc8PFw2m83Z5/fB58L2C9vKMn36dOXm5jqXo0ePeng0AACgsvL3dQEXJCUlae/evdqyZYvX3yswMFCBgYFefx8AAFD5VIojPxMmTNDKlSu1fv161a9f39lutVp17tw5nT592qV/dna2rFars89/X/11Yf1CHwAAgAt8Gn4cDocmTJig5cuXa926dYqOjnbZ3r59e1WvXl1paWnOtszMTB05ckRxcXGSpLi4OO3Zs0cnTpxw9lm7dq3MZrNiYmIqZiAAAKDK8Olpr6SkJC1btkwfffSRateu7ZyjY7FYFBwcLIvFojFjxig5OVkhISEym8164IEHFBcXp86dO0uSevfurZiYGA0fPlzPPPOMbDabZsyYoaSkJE5tAQCAUnx6qbvJZCqzfcmSJRo5cqSk325yOGnSJL399tsqLCxUQkKCFi9e7HJK6z//+Y/Gjx+vDRs2qGbNmhoxYoTmzZsnf//Ly3buXioHAAB8x93v70p1nx9fIfwAAFD1XBH3+QEAAPA2wg8AADAUwg8AADAUwg8AADAUwg8AADAUt8LP0aNHdezYMef6zp07NXHiRL366qseKwwAAMAb3Ao/Q4cO1fr16yX99vDQXr16aefOnXr00Uc1d+5cjxYIAADgSW6Fn7179yo2NlaS9O6776pVq1batm2b3nrrLS1dutST9QEAAHiUW+Hn/PnzzkdHfP7557rlllskSc2bN1dWVpbnqgMAAPAwt8JPy5YtlZqaqs2bN2vt2rVKTEyUJB0/flx169b1aIEAAACe5Fb4efrpp/XKK6+oe/fuGjJkiNq2bStJ+vjjj52nwwAAACojt5/tVVxcLLvdrquuusrZ9uOPP6pmzZoKDQ31WIEVgWd7AQBQ9VTos7169uypvLw8l+AjSSEhIRo8eLA7uwQAAKgQboWfDRs26Ny5c6XaCwoKtHnz5j9dFAAAgLf4l6fzN9984/x53759stlszvXi4mKtXr1aV199teeqAwAA8LByhZ927drJZDLJZDKpZ8+epbYHBwfrhRde8FhxAAAAnlau8HP48GE5HA41atRIO3fudJnYHBAQoLCwMFWrVs3jRQIAAHhKucJPw4YNJUklJSVeKQYAAMDbyhV+fu/AgQNav369Tpw4USoMzZw5808XBgAA4A1uhZ//+Z//0fjx41WvXj1ZrVaZTCbnNpPJRPgBAACVllvh54knntCTTz6pqVOneroeAAAAr3LrPj+//PKL7rjjDk/XAgAA4HVuhZ877rhDn332madrAQAA8Dq3Tns1adJEjz32mLZv367WrVurevXqLtsffPBBjxQHAADgaW492DQ6OvriOzSZ9MMPP/ypoioaDzYFAKDqcff7260jP4cPH3bnZQAAAD7n1pwfAACAqsqtIz+jR4/+w+1vvPGGW8UAAAB4m1vh55dffnFZP3/+vPbu3avTp0+X+cBTAJ5z9lyRYmaukSTtm5ugGgFu36gdAAzJrX81ly9fXqqtpKRE48ePV+PGjf90UQAAAN7isTk/fn5+Sk5O1sKFCz21SwAAAI/z6ITnQ4cOqaioyJO7BAAA8Ci3TnslJye7rDscDmVlZemTTz7RiBEjPFIYAACAN7gVfnbv3u2y7ufnp9DQUP3973+/5JVgAAAAvuRW+Fm/fr2n6wAAAKgQf+oa2ZycHGVmZkqSmjVrptDQUI8UBQAA4C1uTXg+c+aMRo8erYiICHXt2lVdu3ZVZGSkxowZo7Nnz3q6RgC/U1zyf4/j23n4lMs6AODS3Ao/ycnJ2rhxo1asWKHTp0/r9OnT+uijj7Rx40ZNmjTJ0zUC+F+r92YpfsFG5/rIJV+oy9PrtHpvlg+rAoCqxa2nuterV0///ve/1b17d5f29evX684771ROTo6n6qsQPNUdVcHqvVka/8+v9N//w5r+978v3329EltFVHRZAOAz7n5/u3Xk5+zZswoPDy/VHhYWxmkvwAuKSxyas2JfqeAjydk2Z8U+ToEBwGVwK/zExcVp1qxZKigocLb9+uuvmjNnjuLi4i57P5s2bVL//v0VGRkpk8mkDz/80GX7yJEjZTKZXJbExESXPqdOndKwYcNkNptVp04djRkzRvn5+e4MC6i0dh4+pazcgotud0jKyi3QzsOnKq4oAKii3Lraa9GiRUpMTFT9+vXVtm1bSdLXX3+twMBAffbZZ5e9nzNnzqht27YaPXq0Bg4cWGafxMRELVmyxLkeGBjosn3YsGHKysrS2rVrdf78eY0aNUr33Xefli1b5sbIgMrpRN7Fg487/QDAyNwKP61bt9aBAwf01ltv6bvvvpMkDRkyRMOGDVNwcPBl76dPnz7q06fPH/YJDAyU1Wotc9v+/fu1evVqffHFF+rQoYMk6YUXXlDfvn01f/58RUZGXnYtQGUWVjvIo/0AwMjcCj8pKSkKDw/X2LFjXdrfeOMN5eTkaOrUqR4pTpI2bNigsLAwXXXVVerZs6eeeOIJ1a1bV5KUnp6uOnXqOIOPJMXHx8vPz087duzQbbfdVuY+CwsLVVhY6Fy32+0eqxfwhtjoEEVYgmTLLShz3o9JktUSpNjokIouDQCqHLfm/Lzyyitq3rx5qfaWLVsqNTX1Txd1QWJiov7xj38oLS1NTz/9tDZu3Kg+ffqouLhYkmSz2RQWFubyGn9/f4WEhMhms110vykpKbJYLM4lKirKYzUD3lDNz6RZ/WMk/d/VXRdcWJ/VP0bV/P57KwDgv7kVfmw2myIiSl9SGxoaqqwsz91v5K677tItt9yi1q1ba8CAAVq5cqW++OILbdiw4U/td/r06crNzXUuR48e9UzBgBcltorQy3dfrzCz67w3qyWIy9wBoBzcCj9RUVHaunVrqfatW7d6dZ5No0aNVK9ePR08eFCSZLVadeLECZc+RUVFOnXq1EXnCUm/zSMym80uC1AVJLaK0OfJ3ZzrS0d11JapPQk+AFAObs35GTt2rCZOnKjz58+rZ8+ekqS0tDRNmTLFq3d4PnbsmE6ePOk86hQXF6fTp09r165dat++vSRp3bp1KikpUadOnbxWB+BLvz+1FRsdwqkuACgnt8LPI488opMnT+r+++/XuXPnJElBQUGaOnWqpk+fftn7yc/Pdx7FkaTDhw8rIyNDISEhCgkJ0Zw5czRo0CBZrVYdOnRIU6ZMUZMmTZSQkCBJatGihRITEzV27Filpqbq/PnzmjBhgu666y6u9AIAAGVy6/EWF+Tn52v//v0KDg7WtddeW+oePJeyYcMG9ejRo1T7iBEj9PLLL2vAgAHavXu3Tp8+rcjISPXu3VuPP/64y92lT506pQkTJmjFihXy8/PToEGD9Pzzz6tWrVqXXQePt0BVcvZckWJmrpEk7ZuboBoBbv0OAwBVnrvf338q/FwpCD+oSgg/APCbCn22FwAAQFVF+AEAAIbC8XKgiqkR4K8f5/XzdRkAUGVx5AcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABgK4QcAABiKT8PPpk2b1L9/f0VGRspkMunDDz902e5wODRz5kxFREQoODhY8fHxOnDggEufU6dOadiwYTKbzapTp47GjBmj/Pz8ChwFAACoSnwafs6cOaO2bdvqpZdeKnP7M888o+eff16pqanasWOHatasqYSEBBUUFDj7DBs2TN9++63Wrl2rlStXatOmTbrvvvsqaggAAKCKMTkcDoevi5Akk8mk5cuXa8CAAZJ+O+oTGRmpSZMmafLkyZKk3NxchYeHa+nSpbrrrru0f/9+xcTE6IsvvlCHDh0kSatXr1bfvn117NgxRUZGlvlehYWFKiwsdK7b7XZFRUUpNzdXZrPZuwMFAAAeYbfbZbFYyv39XWnn/Bw+fFg2m03x8fHONovFok6dOik9PV2SlJ6erjp16jiDjyTFx8fLz89PO3bsuOi+U1JSZLFYnEtUVJT3BgIAACqVSht+bDabJCk8PNylPTw83LnNZrMpLCzMZbu/v79CQkKcfcoyffp05ebmOpejR496uHoAAFBZ+fu6AF8IDAxUYGCgr8sAAAA+UGmP/FitVklSdna2S3t2drZzm9Vq1YkTJ1y2FxUV6dSpU84+AAAAv1dpw090dLSsVqvS0tKcbXa7XTt27FBcXJwkKS4uTqdPn9auXbucfdatW6eSkhJ16tSpwmsGAACVn09Pe+Xn5+vgwYPO9cOHDysjI0MhISFq0KCBJk6cqCeeeELXXnutoqOj9dhjjykyMtJ5RViLFi2UmJiosWPHKjU1VefPn9eECRN01113XfRKLwAAYGw+DT9ffvmlevTo4VxPTk6WJI0YMUJLly7VlClTdObMGd133306ffq0unTpotWrVysoKMj5mrfeeksTJkzQTTfdJD8/Pw0aNEjPP/98hY8FAABUDZXmPj++5O59AgAAgO9ccff5AQAA8AbCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMBTCDwAAMJRKHX5mz54tk8nksjRv3ty5vaCgQElJSapbt65q1aqlQYMGKTs724cVAwCAyq5Shx9JatmypbKyspzLli1bnNsefvhhrVixQu+99542btyo48ePa+DAgT6sFgAAVHb+vi7gUvz9/WW1Wku15+bm6vXXX9eyZcvUs2dPSdKSJUvUokULbd++XZ07d67oUgEAQBVQ6Y/8HDhwQJGRkWrUqJGGDRumI0eOSJJ27dql8+fPKz4+3tm3efPmatCggdLT0/9wn4WFhbLb7S4LAAAwhkodfjp16qSlS5dq9erVevnll3X48GH95S9/UV5enmw2mwICAlSnTh2X14SHh8tms/3hflNSUmSxWJxLVFSUF0cBAAAqk0p92qtPnz7On9u0aaNOnTqpYcOGevfddxUcHOz2fqdPn67k5GTnut1uJwABAGAQlfrIz3+rU6eOmjZtqoMHD8pqtercuXM6ffq0S5/s7Owy5wj9XmBgoMxms8sCAACMoUqFn/z8fB06dEgRERFq3769qlevrrS0NOf2zMxMHTlyRHFxcT6sEgAAVGaV+rTX5MmT1b9/fzVs2FDHjx/XrFmzVK1aNQ0ZMkQWi0VjxoxRcnKyQkJCZDab9cADDyguLo4rvQAAwEVV6vBz7NgxDRkyRCdPnlRoaKi6dOmi7du3KzQ0VJK0cOFC+fn5adCgQSosLFRCQoIWL17s46oBAEBlZnI4HA5fF+FrdrtdFotFubm5zP8BYEhnzxUpZuYaSdK+uQmqEVCpfzcGJLn//V2l5vwAALyjuOT/fg/eefiUyzpwpSH8AIDBrd6bpfgFG53rI5d8oS5Pr9PqvVk+rArwHsIPABjY6r1ZGv/Pr5RtL3Rpt+UWaPw/vyIA4YpE+AEAgyoucWjOin0q6wTXhbY5K/ZxCgxXHMIPABjUzsOnlJVbcNHtDklZuQXaefhUxRUFVADCDwAY1Im8iwcfd/oBVQXhBwAMKqx2kEf7AVUF4QcADCo2OkQRliCZLrLdJCnCEqTY6JCKLAvwOsIPABhUNT+TZvWPkaRSAejC+qz+Marmd7F4BFRNhB8AMLDEVhF6+e7rFWYOdGm3WoL08t3XK7FVhI8qA7yH+5cDgMEltorQjU3qqfXszyRJS0d11F+uDeWID65YHPkBALgEndjoEIIPrmgc+QEAqEaAv36c18/XZQAVgiM/AADAUAg/AADAUAg/AADAUAg/AADAUAg/AADAUAg/AADAUAg/AADAUAg/AADAUAg/AADAUAg/AADAUAg/AADAUAg/AADAUAg/AADAUAg/AADAUAg/AADAUAg/AADAUAg/AADAUAg/AADAUPx9XQAAALgyFZc4tPPwKZ3IK1BY7SDFRoeomp/J12URfgAAgOet3pulOSv2KSu3wNkWYQnSrP4xSmwV4cPKOO0FAAA8bPXeLI3/51cuwUeSbLkFGv/Pr7R6b5aPKvsN4QcAAHhMcYlDc1bsk6OMbRfa5qzYp+KSsnpUDMIPAADwmJ2HT5U64vN7DklZuQXaefhUxRX1Xwg/AADAY07kXTz4uNPPGwg/AADAY8JqB3m0nzcQfgAAgMfERocowhKki13QbtJvV33FRodUZFkuCD8AAMBjqvmZNKt/jCSVCkAX1mf1j/Hp/X4IPwAAwKMSW0Xo5buvl9XiemrLagnSy3df7/P7/FwxNzl86aWX9Oyzz8pms6lt27Z64YUXFBsb6+uyAAAwpMRWEeoVY+UOz97yzjvvKDk5WampqerUqZMWLVqkhIQEZWZmKiwszNflAQBgSNX8TIprXNfXZZRicjgcvrvLkId06tRJHTt21IsvvihJKikpUVRUlB544AFNmzatVP/CwkIVFhY61+12u6KiopSbmyuz2VxhdQMAAPfZ7XZZLJZyf39X+Tk/586d065duxQfH+9s8/PzU3x8vNLT08t8TUpKiiwWi3OJioqqqHIBAICPVfnw8/PPP6u4uFjh4eEu7eHh4bLZbGW+Zvr06crNzXUuR48erYhSAQBAJXBFzPkpr8DAQAUGBvq6DAAA4ANV/shPvXr1VK1aNWVnZ7u0Z2dny2q1+qgqAABQWVX58BMQEKD27dsrLS3N2VZSUqK0tDTFxcX5sDIAAFAZXRGnvZKTkzVixAh16NBBsbGxWrRokc6cOaNRo0b5ujQAAFDJXBHhZ/DgwcrJydHMmTNls9nUrl07rV69utQkaAAAgCviPj9/lrv3CQAAAL5j2Pv8AAAAlAfhBwAAGArhBwAAGMoVMeH5z7ow7clut/u4EgAAcLkufG+Xd/oy4UdSXl6eJPGMLwAAqqC8vDxZLJbL7s/VXvrtpojHjx9X7dq1ZTKZXLZdeOL70aNHDXElGOO98hltzIz3ysZ4r2yXGq/D4VBeXp4iIyPl53f5M3k48qPfngJfv379P+xjNpsN8RftAsZ75TPamBnvlY3xXtn+aLzlOeJzAROeAQCAoRB+AACAoRB+LiEwMFCzZs1SYGCgr0upEIz3yme0MTPeKxvjvbJ5a7xMeAYAAIbCkR8AAGAohB8AAGAohB8AAGAohB8AAGAohB9JL730kq655hoFBQWpU6dO2rlz5x/2X7RokZo1a6bg4GBFRUXp4YcfVkFBQQVV++eVZ7znz5/X3Llz1bhxYwUFBalt27ZavXp1BVb752zatEn9+/dXZGSkTCaTPvzww0u+ZsOGDbr++usVGBioJk2aaOnSpV6v01PKO96srCwNHTpUTZs2lZ+fnyZOnFghdXpKecf7wQcfqFevXgoNDZXZbFZcXJzWrFlTMcV6QHnHu2XLFt14442qW7eugoOD1bx5cy1cuLBiivUAd/7/vWDr1q3y9/dXu3btvFafN5R3zBs2bJDJZCq12Gy2iin4T3Lnz7iwsFCPPvqoGjZsqMDAQF1zzTV64403yvW+hg8/77zzjpKTkzVr1ix99dVXatu2rRISEnTixIky+y9btkzTpk3TrFmztH//fr3++ut655139Le//a2CK3dPecc7Y8YMvfLKK3rhhRe0b98+jRs3Trfddpt2795dwZW758yZM2rbtq1eeumly+p/+PBh9evXTz169FBGRoYmTpyoe++9t8p8QZZ3vIWFhQoNDdWMGTPUtm1bL1fneeUd76ZNm9SrVy99+umn2rVrl3r06KH+/ftfsX+fa9asqQkTJmjTpk3av3+/ZsyYoRkzZujVV1/1cqWeUd7xXnD69Gndc889uummm7xUmfe4O+bMzExlZWU5l7CwMC9V6FnujPfOO+9UWlqaXn/9dWVmZurtt99Ws2bNyvfGDoOLjY11JCUlOdeLi4sdkZGRjpSUlDL7JyUlOXr27OnSlpyc7Ljxxhu9WqenlHe8ERERjhdffNGlbeDAgY5hw4Z5tU5vkORYvnz5H/aZMmWKo2XLli5tgwcPdiQkJHixMu+4nPH+Xrdu3RwPPfSQ1+rxtvKO94KYmBjHnDlzPF+Ql7k73ttuu81x9913e74gLyvPeAcPHuyYMWOGY9asWY62bdt6tS5vupwxr1+/3iHJ8csvv1RITd50OeNdtWqVw2KxOE6ePPmn3svQR37OnTunXbt2KT4+3tnm5+en+Ph4paenl/maG264Qbt27XKeKvrhhx/06aefqm/fvhVS85/hzngLCwsVFBTk0hYcHKwtW7Z4tVZfSU9Pd/l8JCkhIeGinw+qtpKSEuXl5SkkJMTXpVSI3bt3a9u2berWrZuvS/GaJUuW6IcfftCsWbN8XUqFateunSIiItSrVy9t3brV1+V4zccff6wOHTromWee0dVXX62mTZtq8uTJ+vXXX8u1H0M/2PTnn39WcXGxwsPDXdrDw8P13XfflfmaoUOH6ueff1aXLl3kcDhUVFSkcePGVYnTXu6MNyEhQQsWLFDXrl3VuHFjpaWl6YMPPlBxcXFFlFzhbDZbmZ+P3W7Xr7/+quDgYB9VBm+YP3++8vPzdeedd/q6FK+qX7++cnJyVFRUpNmzZ+vee+/1dUleceDAAU2bNk2bN2+Wv78xvt4iIiKUmpqqDh06qLCwUK+99pq6d++uHTt26Prrr/d1eR73ww8/aMuWLQoKCtLy5cv1888/6/7779fJkye1ZMmSy96PMf52eNCGDRv01FNPafHixerUqZMOHjyohx56SI8//rgee+wxX5fncc8995zGjh2r5s2by2QyqXHjxho1alS5J5cBlc2yZcs0Z84cffTRR1VmfoS7Nm/erPz8fG3fvl3Tpk1TkyZNNGTIEF+X5VHFxcUaOnSo5syZo6ZNm/q6nArTrFkzl/kuN9xwgw4dOqSFCxfq//2//+fDyryjpKREJpNJb731lvNp7gsWLNDtt9+uxYsXX/YvqIYOP/Xq1VO1atWUnZ3t0p6dnS2r1Vrmax577DENHz7c+ZtT69atdebMGd1333169NFH5edXec8kujPe0NBQffjhhyooKNDJkycVGRmpadOmqVGjRhVRcoWzWq1lfj5ms5mjPleQf/3rX7r33nv13nvvlTrNeSWKjo6W9Nu/V9nZ2Zo9e/YVF37y8vL05Zdfavfu3ZowYYKk374oHQ6H/P399dlnn6lnz54+rrJixMbGXrFTEyIiInT11Vc7g48ktWjRQg6HQ8eOHdO11157WfupvN/UFSAgIEDt27dXWlqas62kpERpaWmKi4sr8zVnz54tFXCqVasmSXJU8sekuTPeC4KCgnT11VerqKhI77//vm699VZvl+sTcXFxLp+PJK1du/aSnw+qjrffflujRo3S22+/rX79+vm6nApXUlKiwsJCX5fhcWazWXv27FFGRoZzGTdunJo1a6aMjAx16tTJ1yVWmIyMDEVERPi6DK+48cYbdfz4ceXn5zvbvv/+e/n5+al+/fqXvR9DH/mRpOTkZI0YMUIdOnRQbGysFi1apDNnzmjUqFGSpHvuuUdXX321UlJSJEn9+/fXggULdN111zlPez322GPq37+/MwRVZuUd744dO/TTTz+pXbt2+umnnzR79myVlJRoypQpvhzGZcvPz9fBgwed64cPH1ZGRoZCQkLUoEEDTZ8+XT/99JP+8Y9/SJLGjRunF198UVOmTNHo0aO1bt06vfvuu/rkk098NYRyKe94pd/+obzw2pycHGVkZCggIEAxMTEVXX65lXe8y5Yt04gRI/Tcc8+pU6dOznuhBAcHu/wmWVmVd7wvvfSSGjRooObNm0v67VL/+fPn68EHH/RJ/eVVnvH6+fmpVatWLq8PCwtTUFBQqfbKrLx/xosWLVJ0dLRatmypgoICvfbaa1q3bp0+++wzXw2hXMo73qFDh+rxxx/XqFGjNGfOHP3888965JFHNHr06PIdnf9T14pdIV544QVHgwYNHAEBAY7Y2FjH9u3bndu6devmGDFihHP9/PnzjtmzZzsaN27sCAoKckRFRTnuv//+KnWZYXnGu2HDBkeLFi0cgYGBjrp16zqGDx/u+Omnn3xQtXsuXAb638uFMY4YMcLRrVu3Uq9p166dIyAgwNGoUSPHkiVLKrxud7kz3rL6N2zYsMJrd0d5x9utW7c/7F/ZlXe8zz//vKNly5aOGjVqOMxms+O6665zLF682FFcXOybAZSTO3+ff68qXupe3jE//fTTzu+jkJAQR/fu3R3r1q3zTfFucOfPeP/+/Y74+HhHcHCwo379+o7k5GTH2bNny/W+Joejkp+rAQAA8CBDz/kBAADGQ/gBAACGQvgBAACGQvgBAACGQvgBAACGQvgBAACGQvgBAACGQvgBAACGQvgBUCGWLl2qOnXqVOh7zp49W+3atavQ9wRQ+RF+AHjEyJEjZTKZZDKZFBAQoCZNmmju3LkqKirydWnl8vXXX+uWW25xPhfqmmuu0eDBg3XixAlJ0o8//iiTyaSwsDDl5eW5vLZdu3aaPXu2c7179+7OzyQoKEhNmzZVSkpKpX8IMnClI/wA8JjExERlZWXpwIEDmjRpkmbPnq1nn33W12VdtpycHN10000KCQnRmjVrtH//fi1ZskSRkZE6c+aMS9+8vDzNnz//kvscO3assrKylJmZqenTp2vmzJlKTU311hAAXAbCDwCPCQwMlNVqVcOGDTV+/HjFx8fr448/LrPvoUOHdOuttyo8PFy1atVSx44d9fnnn7v0ueaaa/TUU09p9OjRql27tho0aKBXX33Vpc+xY8c0ZMgQhYSEqGbNmurQoYN27Nhx0fds1KiRJkyYUObRl61btyo3N1evvfaarrvuOkVHR6tHjx5auHChoqOjXfo+8MADWrBggfOI0MXUqFHD+ZmMGjVKbdq00dq1a//wNQC8i/ADwGuCg4N17ty5Mrfl5+erb9++SktL0+7du5WYmKj+/fvryJEjLv3+/ve/q0OHDtq9e7fuv/9+jR8/XpmZmc59dOvWTT/99JM+/vhjff3115oyZYpKSkpKvd8333yjLl26aOjQoXrxxRdlMplK9bFarSoqKtLy5csveWpqyJAhzlN7l8PhcGjz5s367rvvFBAQcFmvAeAdhB8AHudwOPT5559rzZo16tmzZ5l92rZtq7/+9a9q1aqVrr32Wj3++ONq3LhxqSNFffv21f33368mTZpo6tSpqlevntavXy9JWrZsmXJycvThhx+qS5cuatKkie68807FxcW57GPbtm3q3r27Jk+erCeeeOKidXfu3Fl/+9vfNHToUNWrV099+vTRs88+q+zs7FJ9TSaT5s2bp1dffVWHDh266D4XL16sWrVqKTAwUF27dlVJSYkefPDBi/YH4H2EHwAes3LlStWqVUtBQUHq06ePBg8e7DIB+Pfy8/M1efJktWjRQnXq1FGtWrW0f//+Ukd+2rRp4/zZZDLJarU6TzVlZGTouuuuU0hIyEVrOnLkiHr16qWZM2dq0qRJlxzDk08+KZvNptTUVLVs2VKpqalq3ry59uzZU6pvQkKCunTposcee+yi+xs2bJgyMjK0detW9enTR48++qhuuOGGS9YBwHsIPwA8pkePHsrIyNCBAwf066+/6s0331TNmjXL7Dt58mQtX75cTz31lDZv3qyMjAy1bt261Gmy6tWru6ybTCbnaa3g4OBL1hQaGqrY2Fi9/fbbstvtlzWOunXr6o477tD8+fO1f/9+RUZGXnRy87x58/TOO+9o9+7dZW63WCxq0qSJOnbsqHfffVcvvvhiqblNACoW4QeAx9SsWVNNmjRRgwYN5O/v/4d9t27dqpEjR+q2225T69atZbVa9eOPP5br/dq0aaOMjAydOnXqon2Cg4O1cuVKBQUFKSEhodTl6ZcSEBCgxo0bl7ra64LY2FgNHDhQ06ZNu+S+atWqpYceekiTJ0/mcnfAhwg/AHzi2muv1QcffKCMjAx9/fXXGjp0aJkTlf/IkCFDZLVaNWDAAG3dulU//PCD3n//faWnp7v0q1mzpj755BP5+/urT58+ys/PL3N/K1eu1N13362VK1fq+++/V2ZmpubPn69PP/1Ut95660XrePLJJ7Vu3TrnROw/8te//lXff/+93n///XKNFYDnEH4A+MSCBQt01VVX6YYbblD//v2VkJCg66+/vlz7CAgI0GeffaawsDD17dtXrVu31rx581StWrVSfWvVqqVVq1bJ4XCoX79+ZR7JiYmJUY0aNTRp0iS1a9dOnTt31rvvvqvXXntNw4cPv2gdTZs21ejRo1VQUHDJmkNCQnTPPfdo9uzZ5Q57ADzD5ODYKwAAMBCO/AAAAEMh/AAAAEMh/AAAAEMh/AAAAEMh/AAAAEMh/AAAAEMh/AAAAEMh/AAAAEMh/AAAAEMh/AAAAEMh/AAAAEP5/wlsjrYRW6VPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.errorbar(log10y_center,N_clusters_y_theory,yerr = np.sqrt(N_clusters_y_theory),ls='None',marker='o')\n",
    "plt.xlabel('Planck SNR')\n",
    "plt.ylabel('counts')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a322fae1-5f4a-4dda-b954-3a040c48c90e",
   "metadata": {},
   "source": [
    "## Unbinned calculation\n",
    "\n",
    "See [Zubeldia & Bolliet (2024)](https://inspirehep.net/literature/2768906).\n",
    "\n",
    "The class_sz unbinned implementation loops over $z$ and computes in each noise patch in parallel. This is sub-optimal.\n",
    "Given the fact that `cosmocnc` superseeds `class\\_sz` on these calculations, we advise you switch to `cosmocnc` if your work requires \n",
    "heavy cluster cosmology calculation. `class_sz` unbinned cluster calculation can always be used for exploration. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "f38aa8f7-b759-410c-a76e-ab6f62767bf2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 24min 26s, sys: 27.2 s, total: 24min 53s\n",
      "Wall time: 3min 52s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "class_sz = Class_sz()\n",
    "class_sz.set({\n",
    "\n",
    "'output' : 'sz_cluster_counts_fft',\n",
    "'mass_function' :  'T08M500c',\n",
    "'has_selection_function' : 1,\n",
    "'experiment' : 0,\n",
    "'y_m_relation' : 0,\n",
    "'use_skyaveraged_noise': 0,\n",
    "'use_planck_binned_proba' : 0,\n",
    "\n",
    "\n",
    "'m_pivot_ym_[Msun]': 3e14,    \n",
    "\n",
    "'M_min' : 1e13*0.7, \n",
    "'M_max' : 1e16*0.7,\n",
    "\n",
    "'z_min' : 0.01,\n",
    "'z_max' : 1.02,\n",
    "\n",
    "'omega_b': 0.0224178568132,\n",
    "'omega_cdm':  0.11933148326520002,\n",
    "'H0': 70.,\n",
    "'tau_reio': 0.0561,\n",
    "'ln10^{10}A_s': 2.9799585,\n",
    "'n_s': 0.96,\n",
    "\n",
    "\n",
    "\n",
    "'sigmaM_ym' :0.173, # fiducial 0.173\n",
    "\n",
    "\n",
    "'szcc_dof': 0.0,\n",
    "'szcc_qtrunc': -1.0,\n",
    "'szcounts_fft_nz': 550,\n",
    "'szcounts_fft_z_max': 1.01,\n",
    "'szcounts_fft_z_min': 0.01,\n",
    "'szcounts_qmax_fft_padded': 500.0,\n",
    "'N_samp_fftw': 8192,\n",
    "\n",
    "\n",
    "'signal-to-noise_cut-off_for_survey_cluster_completeness' : 6.,\n",
    "\n",
    "\n",
    "# X ray based hydrostatic mass mass bias (if applicable)\n",
    "'B' : 1.25,\n",
    "'cosmo_model': 1,\n",
    "\n",
    "})\n",
    "class_sz.compute_class_szfast()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "f6099657-d4ef-4d5a-a51a-14d4711ef941",
   "metadata": {},
   "outputs": [],
   "source": [
    "zmin = class_sz.pars['szcounts_fft_z_min']\n",
    "zmax = class_sz.pars['szcounts_fft_z_max']\n",
    "nz = class_sz.pars['szcounts_fft_nz']\n",
    "z_arr = np.linspace(zmin,zmax,nz)\n",
    "\n",
    "q_threshold = 6.\n",
    "q_max = 100.\n",
    "nq = 5000\n",
    "q_arr = np.geomspace(q_threshold, q_max,nq)\n",
    "\n",
    "get_dndzdq = np.vectorize(class_sz.get_szcounts_dndzdq_at_z_q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "2237868d-d7c9-4610-8712-cac4f145c217",
   "metadata": {},
   "outputs": [],
   "source": [
    "Nz = []\n",
    "for zp in z_arr:\n",
    "    Nz.append(np.trapz(get_dndzdq(zp,q_arr)*q_arr,x=np.log(q_arr)))\n",
    "Nz = np.asarray(Nz)\n",
    "\n",
    "Nq = []\n",
    "Nq = np.asarray(Nq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "28adeb6d-1acb-4c10-832c-b050e5810bd0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "514.1981625268976"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ntot = np.trapz(Nz,x=z_arr)\n",
    "Ntot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "dadbf2f5-394d-436a-a07c-1145f15ae478",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.787079475530407"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "class_sz.sigma8()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "80c9870a-8394-4b3d-897b-97584146009b",
   "metadata": {},
   "outputs": [],
   "source": [
    "z_edges_low = z_edges[:-1]\n",
    "z_edges_high = z_edges[1:]\n",
    "Nz_binned = []\n",
    "for (zm,zp) in zip(z_edges_low,z_edges_high):\n",
    "    z_arr_bin = np.linspace(zm,zp,50)\n",
    "    Nz_in = []\n",
    "    for zpbin in z_arr_bin:\n",
    "        Nz_in.append(np.trapz(get_dndzdq(zpbin,q_arr)*q_arr,x=np.log(q_arr)))\n",
    "    Nz_binned.append(np.trapz(Nz_in,x=z_arr_bin))\n",
    "Nz_binned = np.asarray(Nz_binned)\n",
    "\n",
    "q_edges_low = log10y_edges[:-1]\n",
    "q_edges_low[0] = 0.6\n",
    "q_edges_high = log10y_edges[1:]\n",
    "Nq_binned = []\n",
    "for (qm,qp) in zip(q_edges_low,q_edges_high):\n",
    "    q_arr_bin = np.linspace(qm,qp,50)\n",
    "    Nq_in = []\n",
    "    for qpbin in q_arr_bin:\n",
    "        if 10.**qpbin<q_threshold:\n",
    "            Nq_in.append(0.)\n",
    "        else:\n",
    "            Nq_in.append(np.trapz(get_dndzdq(z_arr,10**qpbin),x=z_arr))\n",
    "    Nq_binned.append(np.trapz(Nq_in,x=10**q_arr_bin))\n",
    "Nq_binned = np.asarray(Nq_binned)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "d4b5a6d6-357a-4fe5-92c9-9a1d4d575d33",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxUAAAJNCAYAAABHt1gkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABmGklEQVR4nO3deVxVdf7H8fdFFhEE1xKLQA3UTDNIG4PEFRsttWyzxSUxzJpyGTXNAjTXManpNxWJgmVmi6XVWICKC4yNpdI4ZUW5zoiNlYIgAsr5/WHcREAuHuBy4fV8PHjIPed7z/lc70e8b75nsRiGYQgAAAAALpOTvQsAAAAA4NgIFQAAAABMIVQAAAAAMIVQAQAAAMAUQgUAAAAAUwgVAAAAAEwhVAAAAAAwhVABAAAAwBRnexfg6IqLi3X06FE1bdpUFovF3uUAAAAAFTIMQ6dOnVLbtm3l5FR98wuECpOOHj0qX19fe5cBAAAA2OzIkSO6+uqrq217dTJUrFq1Stu3b9euXbu0d+9eFRYWKiEhQWPGjKnwOQcOHND8+fOVnJysY8eOqVmzZrruuus0ceJE3XPPPWXGv/XWW3rppZf09ddfy9XVVSEhIZozZ46CgoKqVGvTpk0lnX9jvLy8qvRc1A9FRUVKTk5WeHi4XFxc7F0O6ij6BLaiV2AregW2uLhPcnJy5Ovra/0MW13qZKiYPXu2Dh06pFatWsnHx0eHDh265PiUlBQNHz5cknTHHXeoffv2OnHihP71r39p48aNZULFvHnzNHv2bPn5+WnChAk6deqU1qxZo1tuuUWbNm1SSEiIzbWWHPLk5eVFqGigioqK1KRJE3l5efFDHRWiT2AregW2oldgi4r6pLoP26+ToSI+Pl4BAQHy8/PTwoULNXPmzArHHj58WHfffbeuuuoqbdy4Uddcc02p9WfPni31ODMzU9HR0QoMDNTOnTvl7e0tSZo4caL+8Ic/aPz48fr3v/9drceYAQAAAPVZnfzkPGDAAPn5+dk0dv78+crJydFrr71WJlBIkrNz6dyUkJCgs2fP6plnnrEGCknq3r27Ro4cqX379iktLc3cCwAAAAAakDoZKmxlGIbee+89tWzZUv369dOuXbu0dOlSLVmyRBs3blRxcXGZ52zZskWSFB4eXmbdoEGDJElbt26t0boBAACA+qROHv5kqwMHDujXX3/VTTfdpMjISL3++uul1t9444366KOPSp3ZnpmZKU9PT7Vp06bM9gICAqxjKlJQUKCCggLr45ycHEnnj1crKioy9XrgmEred95/XAp9AlvRK7AVvQJbXNwnNdUvDh0q/ve//0mS9uzZo2+//VYJCQkaNmyYsrOzNX/+fC1btkx33323Pv/8c+tzsrOzdcUVV5S7vZITrbOzsyvc54IFCxQTE1NmeXJyspo0aWLm5cDBpaSk2LsEOAD6BLaiV2AregW2KOmT06dP18j2HTpUlBzedO7cOc2dO9d6ydnmzZvr9ddf17/+9S/985//VFpamkJDQ6tlnzNnztSUKVOsj0suyxUeHs7VnxqooqIipaSkaODAgVx9AxWiT2AregW2oldgi4v7pOQom+rm0KHiwhOthw4dWmb9HXfcoX/+85/68ssvraHC29u7wpmIkr/kC7d7MTc3N7m5uZVZ7uLiwj/oBo4egC3oE9iKXoGt6BXYoqRPaqpXHPpE7Q4dOqhRo0aSpGbNmpVZX7IsPz/fuiwgIEC5ubk6duxYmfEl51KUnFsBAAAAoHIOHSoaN26sW265RZL0zTfflFlfsszf39+6LCwsTNL5cyAulpSUVGoMAAAAgMo5dKiQpMcee0ySFB0dXeqqTN9++60SExPVtGlT3XbbbdblY8eOlbOzs+bNm1fqMKiMjAy9/fbb6ty5c7WdfwEAAAA0BHXynIr4+HjrDej27t1rXVZyj4nQ0FBFRERIku6//3598MEHev/993XDDTdo0KBBys7O1tq1a3XmzBm98cYbat68uXXbgYGBio6O1uzZs3XDDTdoxIgROnXqlNasWSNJWrZsGXfTBgAAAKqgToaKtLQ0rVy5stSy9PR0paenWx+XhAqLxaK3335bt9xyi5YvX664uDi5ubnplltu0axZs8o9lOmZZ56Rv7+/XnzxRb366qtydXXVrbfeqrlz5yooKKhmXxwAAABQz9TJUJGYmKjExESbxzs7O2vy5MmaPHmyzc958MEH9eCDD15GdQAAAAAuVCdDBQAADVVWVpZ+/PFH7dmzR87OZf+b9vHxkY+Pjx0qA4CKESoAAKhDli1bpueff77C9VFRUYqOjq69ggDABpyRDABAHTJ+/Hi98MIL1ouTSOfPNdy1a5d27dqlyMhI+xUHu0lMTJTFYinVFxXZsmWLLBZLlQ4lrwv8/f3Vp08fe5eBy0SoAACgDvHx8VGHDh10ww03WJd1795dQUFBCgoKsvuhT6dPn9aLL76oW2+9VS1atJCLi4uuvPJKDR48WImJiTp79qyk3z/YLlmyxOZtnzhxQu7u7rJYLHrzzTcrHLd//349+uij6tSpk5o0aaLmzZurc+fOGj16tFJTU02Pr0xVtmexWGz6siUsAHUZhz8BAACb/PDDDxoyZIi+//57DRgwQDNnzlSrVq30v//9Txs3btTYsWP1zTffaPHixZe1/bfeeksFBQVq166dVqxYoYcffrjMmC+//FJhYWFycXHRqFGj1KVLF+Xn5yszM1PJyclq2rSp+vbte9njK1PV7VUWjqKiotSqVSt17NjR5hoq07t3b+Xn58vFxaXatglUhlABAAAqlZ+fr9tvv1379+/X2rVrddddd5VaP2PGDH3xxRf64osvLnsfy5cvV9++fTVs2DBNmjRJ+/fvV/v27UuNiYmJ0enTp5WRkVFqNqfEsWPHTI2vTFW399BDD5W7ndOnT6tXr15q1KiR3nnnnWqdgXJyclLjxo2rbXuALTj8CQCAOujHH3+0fh8VFaXMzEw7VnP+JrTfffedpk6dWiZQlOjRo4cmTpx4WdvfvXu3MjIyNHr0aD3wwANydnbWihUryozLzMxUy5Yty/1AL0lt2rQxNb4y1bW9Rx55RP/617+0aNEi9evXz+b9nz17VtHR0fLz85Onp6eeeuopvfPOO6XGlHdORck5GZs3b9aSJUvUoUMHubm5KTAwsMy9wS5nfEFBgebPn68uXbqocePGatasme644w7t2bOnzNgjR47o3nvvlbe3t7y8vHTHHXeU6nc4JkIF6oWsrCzt3r27wq+srCx7lwgANtu0aZNuvvlm6+OEhAR17tzZrifevv/++5KkRx99tEa2v3z5cnl6emrEiBFq1aqVbr/9dq1cuVLFxcWlxnXo0EG//PKLPvjgA5u2W9XxtbG9v/zlL3rnnXd03333aerUqVV67owZM7RmzRpNnDhRUVFROnv2rB5++GGbe2PWrFl68803FRkZqcWLF8vJyUljxowpdYPhqo4vKirSbbfdppiYGPXq1UuxsbF6+umn9c033ygkJERffvmldezJkyfVu3dvffDBB3r44Ye1cOFCNWnSRH379lVeXl6V/i5QxxgwJTs725BkZGdn27uUBi0qKsqQVOFXVFRUje27sLDQWLdunVFYWFhj+4Djo09gq6+//tpwcnIyRo8ebf0Zdvz4cSMiIsJo1KiRkZmZaZe6WrRoYXh5edk8PjU11ZBk/OUvf6l0bH5+vtGsWTNj9OjR1mXr1q0zJBkbNmwoNfYf//iH4eLiYkgyAgICjLFjxxqvvPKK8c0335S77aqOr4zZ7aWkpBiNGjUyunbtauTm5tq834SEBEOScc011xgnT540DOP8z5W33nrLuOaaa4zmzZsbp0+fNgzj97/7hISEMs/v3r27UVBQYF3+n//8x3B1dTXuv//+cvdny/ilS5cakozPPvus1Days7MNX19fIywszLps5syZhiRjxYoVpcY+9dRThqRSY1E9Lv7/p6Y+uzJTgXohMjJSu3btUlpamnUZl2AE4IgSExPl4eGhhQsXWpe5u7vrr3/9q7y8vLR8+XK71JWTk6OmTZvWyLY/+OADnTx5UqNHj7YuGzx4sFq3bl3mEKhevXpp165dGj16tLKzs5WQkKCJEyfquuuuU+/evbV//35T4ytjZnsHDx7U/fffr6ZNm+rDDz+Uh4dHlfYtSY899pi8vb2tjz08PDR+/HidOHHCpitITZw4Ua6urtbHV111lQIDAys8vM6W8atWrVKnTp0UHBysn3/+2fpVWFiogQMHKi0tTfn5+ZKkdevW6corr9SoUaNK7WfGjBk2vX7UXZyojXqh5A6zF06ddu/e/bJ+YAOAPR06dEj+/v5lTrR1d3dX9+7ddfDgQbvU5eXlpVOnTtXItpcvX67WrVvr6quv1g8//GBdHh4ervfee08///yzWrVqZV3etWtX6+E+hw4d0tatWxUfH6/t27dr2LBh2rVrV6kPwlUdX5nL2d7p06d155136sSJE/rkk0/UoUOHKv4tnde5c+cKl9kSkC4+8V2SWrZsqUOHDl32+H379ik/P1+tW7eucL8///yzfH19tX//fvXo0UONGjUqtd7Hx0fNmjWrtH7UXcxUAABQh/j5+engwYPW3+yWyM/PV0ZGhvz9/e1S1/XXX6+cnJwq/2a/MgcOHFBqaqqOHz+uwMBABQQEWL/eeustFRYWatWqVRU+38/PT6NGjdLWrVsVEhKif//739q5c2e1ja+MrdsbP368MjIyFBMToz/+8Y+XvT+zLv4wX8IwjMsebxiGunbtqpSUlAq/LhU4UD8wUwEAQB0yZswYvfDCC3r66aety/Lz8zVz5kzl5ORo3LhxdqlrxIgR2rZtm+Lj4zV//vxq225CQoIMw9CyZcvK/U317NmztWLFCk2aNOmS27FYLLr55puVnp6u//73v5Xut6rjzWxv6dKlWr16tYYNG6ZnnnnG1H727dunYcOGlVkmlT+rUBsCAgJ0/Phx9evXT05Ol/59dfv27ZWZmalz586VCixZWVk6efJkDVeKmsRMBQAAdUhAQIAef/xxvfXWW9ZlgYGBSkhIUHx8vK699lq71BUREaGOHTtqyZIlWr9+fbljdu3apVdeecXmbRYXFysxMVFdu3ZVRESE7r777jJfI0eO1N69e633v0hJSbHetftC+fn5Sk5OliRdd9111uVVHV+Zqm4vNTVV06dPV8eOHfXGG2/IYrHYvK/yvPrqq8rOzrY+zsvLswaysLAwU9u+XKNGjdKxY8e0dOnSctf/9NNP1u+HDRumn376SW+88UapMYsWLarRGlHzmKkAAKCO6d+/v8aMGaPg4GBJ5+9pMGHCBLsFCklq0qSJPvnkEw0ZMkTDhw9XeHi4Bg4cqJYtW+r48eNKTU1VUlKSpk+fXup5mzZt0pkzZ8psr1WrVvL399eRI0cuOfsyYsQIRUdHa/ny5erRo4cmT56sX375RUOHDlXXrl3VpEkTHTlyRKtXr9b333+vUaNGqWvXrtbnV3V8ZaqyvaysLN177706d+6cRowYoY8++qjC7Xbr1k3dunWrdP+tWrXSzTffrLFjx+rcuXN65ZVX9N///lfx8fFq0qSJza+jOj311FNKSUnRtGnTtHnzZvXr109eXl46fPiwNm3apMaNGys1NVWSNH36dK1evVrjx4/Xrl271KVLF23ZskU7duwodd4MHA+hAgCAOujCQ1liYmLqxIUnrr32Wu3Zs0dxcXFau3at5s2bp9zcXLVo0UI33XSTVq5cqQceeKDUcz777DN99tlnZbbVsWNH64fvim6mJ50/lyMwMFBr1qxRbGysli5dqvXr1ystLU1r167VyZMn5e3trW7dumnGjBkaM2ZMqedXdXxlqrK97777Tj///LMkVXrIWFRUlE2hYtGiRdq+fbv+9re/6aefflKbNm20cuXKMldTqk0uLi76+9//rldeeUVvvvmmoqKiJElt27ZVz549S13Vq3nz5tq+fbumTJlina0ICwtTamqq+vfvb5f6UT0sRkVn5sAmOTk58vb2VnZ2try8vOxdToOXl5cnT09PSVJubm6t/CdcVFSkDRs2aPDgwXJxcanx/cEx0SewVUmvhIWFqXnz5pJq7+cZHAs/V2CLi/ukpj67ck4FAAAAAFM4/AkAgDokKytLP/74Y6krIWVkZMjd3V3S7/flQfU5d+6cjh8/Xum4Fi1aVOl+FkBDQqgAAKAOWbZsmZ5//vlSy0JDQ63fR0VFKTo6uparqt+OHDmidu3aVTouNTVVffr0qfmCAAdEqAAAoA4ZP368mjdvrtDQUDk7l/1vmlmK6temTRulpKRUOu6GG26ohWoAx0SoAACgDvHx8VGHDh104403cvJtLWncuLEGDBhg7zIAh8aJ2gAAAABMIVQAAAAAMIVQAQAAAMAUQgUAAAAAUwgVAAAAAEwhVAAAAAAwhVABAAAAwBRCBQAAAABTCBUAAAAATCFUAAAAADCFUAEAAADAFEIFAACAAzhw4ICGDx+u1q1by2KxaNy4cfYuqU5ITEyUxWLRli1bKh27ZcsWWSwWJSYm1nhd1cnf3199+vSxdxmXRKgAAAA2O336tF588UXdeuutatGihVxcXHTllVdq8ODBSkxM1NmzZyX9/uFtyZIlNm/7xIkTcnd3l8Vi0ZtvvlnhuP379+vRRx9Vp06d1KRJEzVv3lydO3fW6NGjlZqaanp8ZaqyPYvFYtOXLR+Ix4wZo61bt2rGjBl68803NX78+CrXDtQUZ3sXAAAAHMMPP/ygIUOG6Pvvv9eAAQM0c+ZMtWrVSv/73/+0ceNGjR07Vt98840WL158Wdt/6623VFBQoHbt2mnFihV6+OGHy4z58ssvFRYWJhcXF40aNUpdunRRfn6+MjMzlZycrKZNm6pv376XPb4yVd1eZeEoKipKrVq1UseOHS+534KCAm3fvl1PPPGE/vznP0uSioqKtGHDBptrh9S7d2/l5+fLxcXF3qXUO4QKAADqosxM6c03pYMHJX9/6ZFHpIAAu5WTn5+v22+/Xfv379fatWt11113lVo/Y8YMffHFF/riiy8uex/Lly9X3759NWzYME2aNEn79+9X+/btS42JiYnR6dOnlZGRoRtuuKHMNo4dO2ZqfGWqur2HHnqo3O2cPn1avXr1UqNGjfTOO+/Ix8fnkvv96aefZBiGWrRoUWmN586dU0FBgZo0aVLp2IbGyclJjRs3tncZ9RKHPwEAUMdcs2mTnLt2lf7yF+ndd8//2amTZMfjwOPj4/Xdd99p6tSpZQJFiR49emjixImXtf3du3crIyNDo0eP1gMPPCBnZ2etWLGizLjMzEy1bNmy3A/0ktSmTRtT4ytTXdt75JFH9K9//UuLFi1Sv379Ljl2zJgx8vPzk3Q+1JQcMrV161Zt2rRJrq6u2rhxo+bOnasOHTqocePGevfdd63P//nnn/X444/L19dXrq6u8vX11eOPP65ffvml1H5Kzk3YtGmT5syZIz8/P7m7u+vmm2/W559/LknaunWrQkND5eHhIR8fH82dO9em1xsdHS2LxaKDBw+WWXfx+QIldWzevFlLlixRhw4d5ObmpsDAQK1cubLCfZw9e1bR0dHy8/OTm5ubunXrpjVr1pQaU945FVXdX1XHFxQUaP78+erSpYsaN26sZs2a6Y477tCePXvKjD1y5IjuvfdeeXt7y8vLS3fccYd+/PHHCl9zXcJMBQAAdUlmprr/3//JYhhl140bJ4WGStdeW+tlvf/++5KkRx99tEa2v3z5cnl6emrEiBHy8PDQ7bffrpUrV2rOnDlycvr9d6AdOnTQd999pw8++KDCcHOhqo6vje395S9/0TvvvKP77rtPU6dOrXR8ZGSkunfvrsmTJ+vOO++07rdTp0766KOPJEl//vOfVVRUpPHjx8vLy8t6OFV2drZuueUW/fDDD3rkkUcUFBSkPXv26NVXX9XmzZu1c+dONW3atNT+nn76aZ07d05PPfWUCgsL9cILLyg8PFxvvPGGxo0bp0cffVQPPvig3n33XT333HNq165dhTMyZsyaNUv5+fmKjIyUm5ubXn31VY0ZM0bXXnutQkJCyoyfMWOG8vLyrME2ISFBI0eO1JkzZzRmzJhq358t44uKinTbbbfpH//4hx5++GE98cQTys7O1rJlyxQSEqJt27bppptukiSdPHlSvXv31pEjRzRhwgRdd9112rp1q/r27av8/HwTf5O1xIAp2dnZhiQjOzvb3qXAMIy8jAxjvmSsloyCqVMN4/vva3yfhYWFxrp164zCwsIa3xccF30CW52dNs045+RkGFLZr0aNDOPpp+1SV4sWLQwvLy+bx6emphqSjL/85S+Vjs3PzzeaNWtmjB492rps3bp1hiRjw4YNpcb+4x//MFxcXAxJRkBAgDF27FjjlVdeMb755ptyt13V8ZUxu72UlBSjUaNGRteuXY3c3Fyb93vgwAFDkhEVFWVdVlhYaPzpT38yJBmBgYFGXl5emefNmjXLkGT87W9/K7X8//7v/wxJxuzZs63LEhISDEnGjTfeaBQUFFiXr1+/3pBkODs7G1988YV1eUFBgdGmTRvjD3/4Q6X1R0VFGZKMAwcOlFnn5+dnhIWFlamje/fuper4z3/+Y7i6uhr3339/qeeXjL/mmmuMkydPWpefPHnSuOaaa4zmzZsbp0+fNgzj975MSEgwvT9bxi9dutSQZHz22WeltpGdnW34+vqWet0zZ840JBkrVqwoNfapp54yJJUaWxUX//9TU59dOfwJ9UdCgtyDgjRN0r2SGi1dKsPOhwsAQFVZDh06HyHKYxjnz7Gwg5ycnDK/0a4uH3zwgU6ePKnRo0dblw0ePFitW7cucwhUr169tGvXLo0ePVrZ2dlKSEjQxIkTdd1116l3797av3+/qfGVMbO9gwcP6v7771fTpk314YcfysPDo0r7vpTHHnus3HMoPvzwQ7Vu3brMDFNkZKRat26tDz/8sNxtubq6Wh/feuutkqSbb77Z+lt1SXJ1dVXPnj2VmZlZXS+jlIkTJ5aq46qrrlJgYGCF+3vsscfk7e1tfezt7a0JEyboxIkTNl1dq6r7s2X8qlWr1KlTJwUHB+vnn3+2fhUWFmrgwIFKS0uzzkKsW7dOV155pUaNGlVqPzNmzKi09rqAw59QP2RmqnjcODkZxu9NbRgyDEPFjzwiJzsdLgAAVWX4+clisZQfLCyW8ydt24GXl5dOnTpVI9tevny5Wrdurauvvlo//PCDdXl4eLjee+89/fzzz2rVqpV1edeuXa3HxB86dEhbt25VfHy8tm/frmHDhmnXrl2lPuxVdXxlLmd7p0+f1p133qkTJ07ok08+UYcOHar4t3RpgYGB5S4/cOCAbrrpJjk7l/7I5+zsrMDAQO3evbvMcy4+Ob558+aSpHbt2pUZ27x58zLnZlSXi+uQpJYtW+rQoUPlju/cuXOZZdddd50k2RQeq7o/W8bv27dP+fn5at26dYX7/fnnn+Xr66v9+/erR48eatSoUan1Pj4+atasWaX12xszFagXfl2yRMXl/AdskVRsGDpRheukA4A9FY8ZI4thqNy5CsM4f16FHVx//fXKycmp8m/2K3PgwAGlpqbq+PHjCgwMVEBAgPXrrbfeUmFhoVatWlXh8/38/DRq1Cht3bpVISEh+ve//62dO3dW2/jK2Lq98ePHKyMjQzExMfrjH/942furSHVe6eniD7WVLbeFxWKpcF3JvU1s3Z9R0UyeSVXdny3jDcNQ165dlZKSUuHXpQKHI2GmAvXCke3bdf0l1h/etk3Na60aADAhIEB7nnhCN/7tb+dnJgzj9z+XL7fbrOuIESO0bds2xcfHa/78+dW23YSEBBmGoWXLlpX729jZs2drxYoVmjRp0iW3Y7FYdPPNNys9PV3//e9/K91vVceb2d7SpUu1evVqDRs2TM8884zpfVVF+/bt9d133+ns2bOlZivOnj2r77//vtzftteEkkvh/vrrr/K/YLbtzJkzysrK0rXV0Nf79u3TsGHDSi375ptvJJU/q1AbAgICdPz4cfXr16/UBQfK0759e2VmZurcuXOlAktWVpZOnjxZw5Wax0wF6oWDl1hnqWQ9ANQ1R/r319l//1uaNk26997zf373nWTDFWxqSkREhDp27KglS5Zo/fr15Y7ZtWuXXnnlFZu3WVxcrMTERHXt2lURERG6++67y3yNHDlSe/futd7/IiUlpdzfbOfn5ys5OVnS74e8XM74ylR1e6mpqZo+fbo6duyoN95445K/sa8Jw4cP1/HjxxUfH19q+bJly3T8+HHdeeedtVJHyeFZGzduLLU8NjZWxcXF1bKPV199VdnZ2dbH2dnZeu2119SsWTOFhYVVyz6qatSoUTp27JiWLl1a7vqffvrJ+v2wYcP0008/6Y033ig1ZtGiRTVaY3VhpgL1Quatt0r79snQ+RBRomQCMrN3bztUBQAmXHuttGCBvauwatKkiT755BMNGTJEw4cPV3h4uAYOHKiWLVvq+PHjSk1NVVJSkqZPn17qeZs2bdKZM2fKbK9Vq1by9/fXkSNHNO4Sh3SNGDFC0dHRWr58uXr06KHJkyfrl19+0dChQ9W1a1c1adJER44c0erVq/X9999r1KhR6tq1q/X5VR1fmapsLysrS/fee6/OnTunESNGWC//Wp5u3bqpW7duNtdhq+nTp+u9997T448/rt27d+vGG2/Unj17tHz5cnXs2LHM+1VTBgwYoI4dO+q5557TL7/8onbt2iktLU2ff/55qfNlzGjVqpVuvvlmjR07VtL5WbDDhw8rPj7ebjcCfOqpp5SSkqJp06Zp8+bN6tevn7y8vHT48GFt2rRJjRs3VmpqqqTz79Xq1as1fvx47dq1S126dNGWLVu0Y8eOavs7qkl1MlSsWrVK27dv165du7R3714VFhYqISHBpmsM79+/X926dVNeXp4iIyP12muvlTvurbfe0ksvvaSvv/5arq6uCgkJ0Zw5cxQUFFTNrwa1Ydif/6zxy5YpXtJZw5BF5491NIqLNV7SrD//2c4VAoDju/baa7Vnzx7FxcVp7dq1mjdvnnJzc9WiRQvddNNNWrlypR544IFSz/nss8/02WefldlWx44drR++L3W/h+uvv16BgYFas2aNYmNjtXTpUq1fv15paWlau3atTp48KW9vb3Xr1k0zZswo81mhquMrU5Xtfffdd/r5558lqdJDxqKiomokVHh7eys9PV1RUVH66KOPlJCQoCuvvFITJkxQTExMjV3R62KNGjXSRx99pCeffFIvv/yyXF1dFR4ebj0XpTosWrRI27dv19/+9jf99NNPCgwM1FtvvVWmJ2uTi4uL/v73v+uVV17Rm2++qaioKElS27Zt1bNnz1JXPGvevLm2b9+uKVOmWGcrwsLClJqaqv79+9ul/qqwGDV1tosJ/v7+OnTokFq1aiUPDw8dOnTIplBRXFysPn36aPfu3ZcMFfPmzdPs2bPl5+enESNG6NSpU1qzZo0KCwu1adOmKjV3Tk6OvL29lZ2dLS8vr6q+VFSjxMRELRg3TmOKi+Uv6VjjxnqtsFAzly+v8n8aVVFUVKQNGzZo8ODBcnFxqbH9wLHRJ7AVvQJb0SuwxcV9UlOfXevkORXx8fE6ePCgjh8/rgkTJtj8vNjYWO3YsUPPP/98hWMyMzMVHR2twMBAffXVV3rhhRf0+uuva9u2bZLOX52huo7tQ+0aM2aM3tuzR7MkPSDpv48/rr9/912NBgoAAADU0cOfBgwYUOXnfPvtt5o9e7Zmzpyp7t27VzguISFBZ8+e1TPPPFPqBindu3fXyJEjlZiYqLS0NPXmGHyHdOF1v2NiYqr1xkIAgPrp3LlzOn78eKXjWrRoUaX7WQANSZ2cqaiqc+fOafTo0QoICNDs2bMvObbkjorh4eFl1g0aNEiStHXr1mqvEQAA1E1HjhyRj49PpV//+Mc/7F0qUGfVyZmKqlqwYIF2796tzz//vNLfIGRmZsrT01Nt2rQpsy4gIMA6piIFBQUqKCiwPs7JyZF0/ni1oqKiyykf1ejC96C23pOSffD+41LoE9iKXql9LVu21KefflrpuOuuu65OvS/0CmxxcZ/UVL84fKj46quvNGfOHE2bNk3BwcGVjs/OztYVV1xR7rqSk1UuvMbxxRYsWKCYmJgyy5OTk2vlcmW//vqrTpw4UeH65s2bW28w0xBdeNnCpKQkNW7cuNb2nZKSUmv7guOiT2AreqXu2bFjh71LKBe9AluU9Mnp06drZPsOHSoKCws1evRoXXvttdZLdNW0mTNnasqUKdbHOTk58vX1VXh4eK1c/WnOnDmXPBF99uzZeuihh2q8jroqLy/P+v2gQYNq5ZyKoqIipaSkaODAgVx9AxWiT2AregW2oldgi4v7pOQom+rm0KFiwYIF2rt3r/7xj3/Izc3NpueUXEKrPCV/yReewH0xNze3cvfl4uJSK/+gJ06cqDvvvFP5+fkKDQ2VJKWlpcnd3V2S5OPj06B/sFz42mvrPbHX/uCY6BPYil6BregV2KKkT2qqVxz6RO09e/aouLhYf/jDH2SxWKxfffv2lSTFxcXJYrFo+PDh1ucEBAQoNzdXx44dK7O9knMpSs6tqIt8fHwUFBRU6gpX3bt3V1BQkIKCguTj42O/4gAAANAgOfRMxcCBA8u9bXlWVpY2bNigTp06KSQkRDfeeKN1XVhYmHbs2KHk5GSNGjWq1POSkpKsYwAAAADYxqFDxeOPP17u8i1btmjDhg0KCwsrc0ftsWPHasmSJZo3b56GDRtmPdQpIyNDb7/9tjp37mw9rAgAAABA5epkqIiPj1daWpokae/evdZlJfeYCA0NVURExGVtOzAwUNHR0Zo9e7ZuuOEGjRgxQqdOndKaNWskScuWLZOTk0MfFQYAAADUqjoZKtLS0rRy5cpSy9LT05Wenm59fLmhQpKeeeYZ+fv768UXX9Srr74qV1dX3XrrrZo7d66CgoIue7sAAABAQ1QnQ0ViYqISExMv+/l9+vSRYRiXHPPggw/qwQcfvOx9AAAAADiP43wAAAAAmEKoAAAAAGAKoQIAAACAKYQKAAAAAKYQKgAAAACYQqgAAAAAYAqhAgAAAIAphAoAAAAAphAqAAAAAJhCqAAAAABgCqECAAAAgCmECgAAAACmECoAAAAAmEKoAAAAAGAKoQIAAACAKYQKAAAAAKYQKgAAAACYQqgAAAAAYIqzvQsAqkNWVpaysrKUn59vXZaRkSF3d3dJko+Pj3x8fOxVHgAAQL3GTIWDsvzwg+ZLWi3JJSpKysy0d0l2FRcXp+DgYIWGhlqXhYaGKjg4WMHBwYqLi7NjdQAAAPUbMxWOKCFB7hERmibJIsnpxRel2Fhp+XJpzBj71mYnkZGRGjp0aIXrmaUAAACoOYQKR5OZKUVEyFJc/Pubd+7c+T/HjZNCQ6Vrr7VXdXbD4U0AAAD2w+FPjmbFCsliKX+dxXJ+tgIAAACoRYQKR3PwoIqLi8tdVVxcLB08WLv1AAAAoMEjVDiYX728VGwY5a4rNgyd8Pau5YoAAADQ0BEqHMwKnT85++JYYfy2nIOfAAAAUNsIFQ5mV06OlnTqJDk5qUjSWUlGo0ayODlpSadO2pWdbe8SAQAA0MAQKhyMv7+/Fv/vfzqxY4eWSHpPUtGkSTrz1Vda9NNP8vf3t3OFAAAAaGgIFQ7mkUceUXZ2tqa++qpmSXpAUs7TT+tPL72knJwcjRs3zt4lAgAAoIHhPhUOJiAgQPHx8YqIiLAuCwwMVE5OjuLj43VtA7xHBQAAAOyLmQoHNGbMGO3evdv6+JFHHtG3336rMQ30btoAAACwL2YqHFSHDh2s38fExMjDw8OO1QAAAKAhY6YCAAAAgCmECgAAAACmECoAAAAAmEKoAAAAAGAKoQIAAACAKYQKAAAAAKYQKgAAAACYQqgAAAAAYAqhAgAAAIAphAoAAAAAphAqAAAAAJhCqAAAAABgCqECAAAAgCl1MlSsWrVKkZGRuummm+Tm5iaLxaLExMQy44qKirR27VqNHj1anTt3lqenp5o2baqbb75Zr776qs6dO1fhPt566y317NlTHh4eat68uW6//Xbt3r27Bl8VAAAAUD8527uA8syePVuHDh1Sq1at5OPjo0OHDpU77scff9Tdd98tT09P9e/fX0OHDlV2drY+/vhjTZw4URs2bNBHH30ki8VS6nnz5s3T7Nmz5efnpwkTJujUqVNas2aNbrnlFm3atEkhISG18TIBAACAeqFOzlTEx8fr4MGDOn78uCZMmFDhuKZNm+pvf/ubjh07pnXr1mnRokV67bXX9P333+umm27SJ598ovfff7/UczIzMxUdHa3AwEB99dVXeuGFF/T6669r27ZtkqTx48eruLi4Rl8fAAAAUJ/UyVAxYMAA+fn5VTruqquu0sSJE+Xh4VFquYeHh6ZMmSJJ2rp1a6l1CQkJOnv2rJ555hl5e3tbl3fv3l0jR47Uvn37lJaWVg2vAgAAAGgY6mSoqA4uLi6SJGfn0kd4bdmyRZIUHh5e5jmDBg2SVDaIAAAAAKhYvQ0VK1askFQ2PGRmZsrT01Nt2rQp85yAgADrGAAAAAC2qZMnapv1+uuv69NPP1W/fv00ePDgUuuys7N1xRVXlPs8Ly8v65iKFBQUqKCgwPo4JydH0vkrURUVFZkt3WYX7qu2943SSv7ueQ9wKfQJbEWvwFb0CmxxcZ/UVL/Uu1DxySef6IknnpCfn59WrVpV7dtfsGCBYmJiyixPTk5WkyZNqn1/FTlz5oz1+6SkJDVu3LjW9o3ypaSk2LsEOAD6BLaiV2AregW2KOmT06dP18j261Wo2LBhg+6++25deeWV2rx5s3x8fMqM8fb2rnAmomTW4cITuC82c+ZM60ngJc/x9fVVeHi4daajNuTl5Vm/HzRoUJmT1VF7ioqKlJKSooEDB1rP5QEuRp/AVvQKbEWvwBYX90nJ593qVm9Cxd///neNGDFCrVq1Umpqqtq3b1/uuICAAO3YsUPHjh0rc15FybkUJedWlMfNzU1ubm5llru4uNTqP+gL91Xb+0b5eB9gC/oEtqJXYCt6BbYo6ZOa6pV6caJ2SaBo0aKFUlNTde2111Y4NiwsTNL5w5UulpSUVGoMAAAAgMo5fKj49NNPNWLECDVv3lypqamXnGWQpLFjx8rZ2Vnz5s0rdRhURkaG3n77bXXu3FmhoaE1XTYAAABQb9TJw5/i4+OtN6Dbu3evdVnJPSZCQ0MVERGhb7/9VnfeeacKCgrUp08fvf3222W25e/vrzFjxlgfBwYGKjo6WrNnz9YNN9ygESNG6NSpU1qzZo0kadmyZXJycvisBQAAANSaOhkq0tLStHLlylLL0tPTlZ6ebn0cERGhY8eOWS/vWhIKLhYWFlYqVEjSM888I39/f7344ot69dVX5erqqltvvVVz585VUFBQ9b4YAAAAoJ6rk6EiMTFRiYmJlY7r06ePDMO4rH08+OCDevDBBy/rufaUlZWlrKws5efnW5dlZGTI3d1dkuTj41PuVa8AAACAmsJxPg4mLi5OwcHBpc77CA0NVXBwsIKDgxUXF2fH6gAAANAQ1cmZClQsMjJSQ4cOrXA9sxQAAACobYQKB8PhTQAAAKhrOPwJAAAAgCmECgAAAACmECoAAAAAmEKoAAAAAGAKoQIAAACAKYQKAAAAAKYQKgAAAACYQqgAAAAAYAqhAgAAAIAphAoAAAAAphAqAAAAAJhCqAAAAABgCqECAAAAgCmECgAAAACmECoAAAAAmEKoAAAAAGAKoQIAAACAKYQKAAAAAKYQKgAAAACYQqgAAAAAYAqhAgAAAIApzvYuAEDNyMrKUlZWVoXrfXx85OPjU4sVAQCA+opQAdRTcXFxiomJqXB9VFSUoqOja68gAABQbxEqgHoqMjJSQ4cOVX5+vkJDQyVJaWlpcnd3lyRmKQAAQLUhVAD1VMnhTXl5edZl3bt3l4eHhx2rAgAA9REnagMAAAAwhVABAAAAwBRCBQAAAABTCBUAAAAATCFUAAAAADCFUAEAAADAFEIFAAAAAFMIFQAAAABMIVQAAAAAMIVQAQAAAMAUQgUAAAAAUwgVAAAAAEwhVAAAAAAwhVABAAAAwBRCBQAAAABTCBUAAAAATCFUAAAAADCFUAEAAADAlDoZKlatWqXIyEjddNNNcnNzk8ViUWJiYoXjc3JyNGXKFPn5+cnNzU3+/v6aNm2acnNzyx1fXFysl19+WV27dpW7u7tat26tkSNHav/+/TX0igAAAID6q06GitmzZ+v111/XoUOH5OPjc8mxeXl5CgsLU2xsrDp16qTJkyerY8eOWrJkifr166czZ86UeU5kZKSefPJJGYahJ598Urfddps++OAD9ejRQ5mZmTX1sgAAAIB6qU6Givj4eB08eFDHjx/XhAkTLjl28eLFysjI0IwZM5SUlKSFCxcqKSlJM2bM0BdffKHY2NhS41NTUxUfH6/evXtr9+7dWrRokd58802tW7dOv/76q5544omafGkAAABAvVMnQ8WAAQPk5+dX6TjDMBQfHy9PT089++yzpdY9++yz8vT0VHx8fKnly5YtkyTNnTtXrq6u1uV//OMf1adPHyUnJ+vw4cPV8CoAAACAhqFOhgpbZWZm6ujRowoJCZGHh0epdR4eHgoJCdH+/ft15MgR6/ItW7ZY111s0KBBkqStW7fWbOEAAABAPeJs7wLMKDn/ISAgoNz1AQEBSkpKUmZmpnx9fZWXl6esrCxdf/31atSoUbnjL9xueQoKClRQUGB9nJOTI0kqKipSUVHRZb8WOK6S972uvv8HkpM1X5K/pB19+sh/zhz5DRhg56oanrreJ6g76BXYil6BLS7uk5rqF4cOFdnZ2ZIkb2/vctd7eXmVGlfV8eVZsGCBYmJiyixPTk5WkyZNbKwc9VFKSoq9Sygj9+WXdf+mTeokySLJ+PJLWQYP1tsDBqgp5w/ZRV3sE9RN9ApsRa/AFiV9cvr06RrZvkOHCnuYOXOmpkyZYn2ck5MjX19fhYeHW0MJGpaioiKlpKRo4MCBcnFxsXc5Voc2blS7TZt08ZycIWnkxo06MHWq/Pr3t0dpDVJd7RPUPfQKbEWvwBYX90nJUTbVzaFDRcmMQ0UzCyV/aSXjqjq+PG5ubnJzcyuz3MXFhX/QDVxd64H/zJkj/3KWW3Q+WPwnJkbX3nZb7RaFOtcnqLvoFdiKXoEtSvqkpnrFoU/UruwciIvPufDw8JCPj48OHDigc+fOVToecGSuR4/KUsE6y2/rAQAAqoPDh4q2bdsqPT1deXl5pdbl5eUpPT1d7dq1k6+vr3V5WFiYdd3FkpKSJEm9e/eu2cKBWlDYtq2MCtYZv60HAACoDg4dKiwWiyIiIpSbm6u5c+eWWjd37lzl5uZq/PjxpZY/+uijks7fx6KwsNC6/NNPP9WWLVsUHh5u0z0ygLrOf84c66FOFzJ0fqai3fPP135RAACgXqqT51TEx8crLS1NkrR3717rsi1btkiSQkNDFRERIUmaPn261q9fr0WLFmnPnj0KCgrS7t27lZycrB49emjSpEmltt23b19FREQoPj5eQUFBGjJkiLKysvTOO++oRYsWevnll2vtdQI1yX/gQKWNG6dey5erWLIeCmWRtGPcOIVykjYAAKgmdTJUpKWlaeXKlaWWpaenlzpkqSRUeHh4aOvWrYqOjtbatWuVmpoqHx8fTZ06VVFRUXJ3dy+z/bi4OHXt2lWvv/66XnrpJXl6eurOO+/UvHnz1KFDh5p9cUAtCo2P17fDh+vDO+6Qv6QrevbUtfPnEygAAEC1qpOhIjExUYmJiTaP9/b2VmxsrGJjY20a7+TkpCeffFJPPvnkZVYIOA7fvn0167fvczdvLnP3eQAAALMc+pwKAAAAAPZHqAAAAABgCqECAAAAgCmECgAAAACmECoAAAAAmEKoAAAAAGAKoQIAAACAKYQKAAAAAKYQKgAAAACYQqgAAAAAYAqhAgAAAIAphAoAAAAAphAqAAAAAJhCqAAAAABgCqECAAAAgCmECgAAAACmECoAAAAAmEKoAAAAAGAKoQIAAACAKYQKAAAAAKYQKgAAAACYQqgAAAAAYAqhAgAAAIApzvYuAEDNyMrKUlZWlvLz863LMjIy5O7uLkny8fGRj4+PvcoDAAD1CDMVQD0VFxen4OBghYaGWpeFhoYqODhYwcHBiouLs2N1AACgPmGmAqinIiMjNXTo0ArXM0sBAACqC6ECqKc4vAkAANQWDn8CAAAAYAqhAgAAAIAphAoAAAAAphAqAAAAAJhCqAAAAABgCqECAAAAgCmECgAAAACmECoAAAAAmEKoAAAAAGAKoQIAAACAKYQKAAAAAKYQKgAAAACYQqgAAAAAYAqhAgAAAIAphAoAAAAAphAqAAAAAJhCqAAAAABgCqECAAAAgCmECgAAAACm1ItQYRiGPvjgA/Xt21c+Pj5q0qSJOnbsqMjISO3fv7/M+JycHE2ZMkV+fn5yc3OTv7+/pk2bptzcXDtUDwAAADi2ehEq/vznP2vEiBH67rvvNHz4cP3pT39Su3bttGzZMnXv3l3//ve/rWPz8vIUFham2NhYderUSZMnT1bHjh21ZMkS9evXT2fOnLHjKwEAAAAcj7O9CzDr2LFjevHFF+Xn56evvvpK3t7e1nWxsbGaMmWKli5dqhUrVkiSFi9erIyMDM2YMUMLFy60jn366ae1aNEixcbGaubMmbX+OgAAAABH5fAzFQcPHlRxcbFCQkJKBQpJuv322yVJx48fl3T+MKn4+Hh5enrq2WefLTX22Weflaenp+Lj42uncAAAAKCecPhQERAQIFdXV6WnpysnJ6fUuk8++USS1L9/f0lSZmamjh49qpCQEHl4eJQa6+HhoZCQEO3fv19HjhypneIBAACAesDhD39q2bKlFi5cqKlTp6pTp04aNmyYvLy89NVXX2nz5s2aOHGinnjiCUnnQ4V0PoiUJyAgQElJScrMzJSvr2+5YwoKClRQUGB9XBJkioqKVFRUVJ0vDQ6i5H3n/cel0CewFb0CW9ErsMXFfVJT/eLwoUKSJk+erKuuukoRERF67bXXrMtDQ0P1wAMPyNn5/MvMzs6WpDKHSZXw8vIqNa48CxYsUExMTJnlycnJatKkyWW/Bji+lJQUe5cAB0CfwFb0CmxFr8AWJX1y+vTpGtl+vQgVc+bM0fPPP685c+booYceUrNmzZSRkaHJkyerT58+Wrt2rYYOHVot+5o5c6amTJlifZyTkyNfX1+Fh4dbQwkalqKiIqWkpGjgwIFycXGxdzmoo+gT2Ipega3oFdji4j65+HSB6mJzqPDy8lJQUJBuvPFGBQcHKygoSJ07d5bFYqmRwmy1ceNGRUVFafLkyXr66aety0NDQ/Xxxx+rffv2mjp1qoYOHWqdoahoJqLkL7mimQxJcnNzk5ubW5nlLi4u/INu4OgB2II+ga3oFdiKXoEtSvqkpnrF5lDh7u6ubdu2adu2bdYg4e7urhtuuEHBwcHWoNGlSxc5OdXe+d+ffvqpJKlv375l1rVp00adOnXSnj17lJubaz2XouTciotVds4FAAAAgLJsDhU//fST/vvf/2rXrl3avXu3du3apV27dmnHjh3asWOHNWg0btxYXbt2tYaMcePG1VjxklRYWCjp98vGXuz48eNycnKSi4uLAgIC1LZtW6WnpysvL6/UFaDy8vKUnp6udu3aVXiSNgAAAICyqjSlcNVVV2no0KGKjo7Wxx9/rKNHj+ro0aP6+OOPFRUVpdtvv10tWrTQzp079dprrykyMrKm6rYKCQmRJC1durTMYU2vvfaa/vOf/6hXr15yc3OTxWJRRESEcnNzNXfu3FJj586dq9zcXI0fP77GawYAAADqE9Mnardp00ZDhgzRkCFDVFBQoA0bNuiVV17Rpk2bqqO+St1zzz169dVXtW3bNgUGBmro0KFq1qyZdu/erc2bN8vd3V1Lly61jp8+fbrWr1+vRYsWac+ePQoKCtLu3buVnJysHj16aNKkSbVSNwAAAFBfmA4Vp0+f1t///ne9//77+vTTT5WXlyfDMBQcHKwRI0ZUR42X1KhRIyUnJys2NlbvvvuuVq9ercLCQl155ZV66KGHNGvWLHXu3Nk63sPDQ1u3blV0dLTWrl2r1NRU+fj4aOrUqYqKipK7u3uN1wwAAADUJ5cVKk6dOqWPP/5Y77//vpKSknTmzBlJUq9evTRixAjddddd8vPzq9ZCL8XNzU1PP/10qas/XYq3t7diY2MVGxtbw5UBAAAA9Z/NoeLkyZNav3693n//fW3cuFEFBQVq1KiRwsLCNGLECN15551q06ZNTdYKAAAAoA6yOVRceeWVOnv2rFxdXdW/f3/dddddGjZsmFq2bFmT9QEAAACo42y++lNRUZEkqU+fPurdu7euueaaWr0fBQAAAIC6yeaZihtvvFFff/21kpKSlJycbF1+zTXXKCgoyHpfiqCgIF1xxRU1UiwAAACAusfmULFr1y4VFRVp79691pvf7d69W//617/04Ycf6sMPP7TeAK9t27bWoPHcc8/VWPEAAAAA7K9KV39ycXGxzkZERERIks6dO6evv/7aeoft3bt366uvvtLHH3+sTz75hFABAAAA1HOm71PRqFEjdevWTd26ddPYsWMlScXFxdq3b5927dplukAAAAAAdZvpUFEeJycndenSRV26dKmJzQMAAACoQ7h8EwAAAABTCBUAAAAATCFUAAAAADCFUAEAAADAFEIFAAAAAFMIFQAAAABMIVQAAAAAMIVQAQAAAMAUQgUAAAAAUwgVAAAAAEwhVAAAAAAwhVABAAAAwBRCBQAAAABTnO1dAADUlqysLGVlZVW43sfHRz4+PrVYEQAA9QOhAkCDERcXp5iYmArXR0VFKTo6uvYKAgCgniBUAGgwIiMjNXToUOXn5ys0NFSSlJaWJnd3d0lilgIAgMtEqADQYJQc3pSXl2dd1r17d3l4eNixKgAAHB8nagMAAAAwhVABAAAAwBRCBQAAAABTCBUAAAAATCFUAAAAADCFUAEAAADAFEIFAAAAAFMIFQAAAABMIVQAAAAAMIVQAQAAAMAUQgUAAAAAUwgVAAAAAEwhVAAAAAAwhVABAAAAwBRCBQAAAABTCBUAAAAATCFUAGhwDm/apPmSVkv6vF8/HUxJsXdJAAA4NEIFgAZl+yOPKHDYME2TdK+ksJ075RserrSICHuXBgCAwyJUAGgwDqak6JaEBDWS5CxZ/3SS1Gv5ch3atMmu9QEA4KgIFQAajIPPPSejnOUWSYakA7Nn13JFAADUD4QKAA2G69GjslSwzvLbegAAUHX1KlR8+OGHGjhwoFq2bKnGjRurXbt2GjlypI4cOVJqXE5OjqZMmSI/Pz+5ubnJ399f06ZNU25urp0qB1AbCtu2LXemQjo/U1HYtm1tlgMAQL3hbO8CqoNhGJowYYJef/11dejQQffff7+aNm2qo0ePauvWrTp06JB8fX0lSXl5eQoLC1NGRobCw8M1cuRI7dmzR0uWLNHWrVu1bds2NW7c2M6vCEBN8J8zR5bwcBlSqRmLksftnn/ePoUBAODg6kWo+Otf/6rXX39dEydO1F//+lc1atSo1PqzZ89av1+8eLEyMjI0Y8YMLVy40Lr86aef1qJFixQbG6uZM2fWWu0Aao//wIFKGzdOvZYvV7F+DxYWSTvGjVNo//52rA4AAMfl8Ic/5efnKyYmRu3bt9dLL71UJlBIkrPz+exkGIbi4+Pl6empZ599ttSYZ599Vp6enoqPj6+VugHYR2h8vDI//lhLJL0naWvPnvrPxo0K5d8+AACXzeFnKpKTk3XixAmNHTtW586d00cffaTvv/9ezZo104ABA3Tttddax2ZmZuro0aMaNGiQPDw8Sm3Hw8NDISEhSkpK0pEjR6yHSwGof3z79tWs377P3by5zM8DAABQNQ4fKnbt2iVJatSokbp166bvv//eus7JyUmTJ0/WkiVLJJ0PFZIUEBBQ7rYCAgKUlJSkzMzMCkNFQUGBCgoKrI9zcnIkSUVFRSoqKjL/guBwSt533n/HceF7VVv/dukT2Ipega3oFdji4j6pqX5x+FDxv//9T5K0dOlSBQUFaefOnercubP27NmjRx99VC+88II6dOigxx57TNnZ2ZIkb2/vcrfl5eUlSdZx5VmwYIFiYmLKLE9OTlaTJk3Mvhw4sJSUFHuXABudOXPG+n1SUlKtXpyBPoGt6BXYil6BLUr65PTp0zWyfYcPFcXFxZIkV1dXrVu3Tm1/uyTkrbfeqvfee0833HCDXnjhBT322GPVsr+ZM2dqypQp1sc5OTny9fVVeHi4NZSgYSkqKlJKSooGDhwoFxcXe5cDG+Tl5Vm/L+9wyJpAn8BW9ApsRa/AFhf3SclRNtXN4UNFyazDTTfdZA0UJa6//nq1b99eP/zwg06ePGkdW9FMRMlfckUzGZLk5uYmNze3MstdXFz4B93A0QOO48L3qbbfN/oEtqJXYCt6BbYo6ZOa6hWHv/pTx44dJUnNmjUrd33J8vz8fOu5FCXnVlyssnMuAAAAAJTl8DMVffv2lSTt27evzLqioiL98MMP8vDwUOvWrdWmTRu1bdtW6enpysvLK3XIQ15entLT09WuXTuu/AQAAABUgcPPVHTo0EHh4eH64YcfytxjYuHChTp58qTuvPNOOTs7y2KxKCIiQrm5uZo7d26psXPnzlVubq7Gjx9fm+UDAAAADs/hZyok6ZVXXtEtt9yi8ePHa926derUqZP27NmjzZs3y8/PT3/5y1+sY6dPn67169dr0aJF2rNnj4KCgrR7924lJyerR48emjRpkv1eCAAAAOCAHH6mQjo/W/Hll19qzJgx2rVrl/76178qMzNTjz/+uHbu3Kk2bdpYx3p4eGjr1q2aNGmS9u3bpxdeeEHffvutpk6dqk2bNsnd3d2OrwQAAABwPPVipkKSfH19lZCQYNNYb29vxcbGKjY2toarAgAAAOq/ejFTAQAAAMB+CBUAAAAATCFUAAAAADCFUAEAAADAFEIFAAAAAFMIFQAAAABMIVQAAAAAMKXe3KcCACqTlZWlrKws5efnW5dlZGRYb3rp4+MjHx8fe5UHAIDDYqYCQIMRFxen4OBghYaGWpeFhoYqODhYwcHBiouLs2N1AAA4LmYqADQYkZGRGjp0aIXrmaUAAODyECoANBgc3gQAQM3g8CcAAAAAphAqAAAAAJhCqAAAAABgCqECAAAAgCmECgAAAACmECoAAAAAmEKoAAAAAGAKoQIAAACAKYQKAAAAAKYQKgAAAACYQqgAAAAAYAqhAgAAAIAphAoAAAAAphAqAAAAAJhCqAAAAABgCqECAAAAgCmECgAAAACmECoAAAAAmEKoAAAAAGAKoQIAAACAKYQKAAAAAKYQKgAAAACYQqgAAAAAYAqhAgAAAIAphAoAAAAAphAqAAAAAJhCqAAAAABgCqECAAAAgCmECgAAAACmECoAAAAAmOJs7wIAAPaVlZWlrKysCtf7+PjIx8enFisCADgaQgUANHBxcXGKiYmpcH1UVJSio6NrryAAgMMhVABAAxcZGamhQ4cqPz9foaGhkqS0tDS5u7tLErMUAIBKESoAoIErObwpLy/Puqx79+7y8PCwY1UAAEdSb0/UXrRokSwWiywWiz7//PMy63NycjRlyhT5+fnJzc1N/v7+mjZtmnJzc+1QLQAAAOC46mWo+Pe//62oqKgKf8uWl5ensLAwxcbGqlOnTpo8ebI6duyoJUuWqF+/fjpz5kwtVwwAAAA4rnoXKoqKijR69Gh1795dd955Z7ljFi9erIyMDM2YMUNJSUlauHChkpKSNGPGDH3xxReKjY2t5aoBAAAAx1XvQsW8efP09ddfa8WKFWrUqFGZ9YZhKD4+Xp6ennr22WdLrXv22Wfl6emp+Pj42ioXAAAAcHj1KlTs3r1b8+bNU1RUlK677rpyx2RmZuro0aMKCQkpc3iUh4eHQkJCtH//fh05cqQ2SgYAAAAcXr0JFQUFBRo1apS6d++u6dOnVzguMzNTkhQQEFDu+pLlJeMAAAAAXFq9uaTsc889p8zMTO3atavcw55KZGdnS5K8vb3LXe/l5VVq3MUKCgpUUFBgfZyTkyPp/LkcRUVFl1U7HFvJ+877j0txhD65sDZ+ptmPI/QK6gZ6Bba4uE9qql/qRajYsWOHlixZoujoaF1//fU1uq8FCxaUe+fZ5ORkNWnSpEb3jbotJSXF3iXAAdTlPrnwyndJSUlq3LixHatBXe4V1C30CmxR0ienT5+uke07fKg4e/asRo8erW7duunpp5+udHzJDEVFMxElMw8VzWTMnDlTU6ZMKTXe19dX4eHh1lkONCxFRUVKSUnRwIED5eLiYu9yUEc5Qp9cePO7QYMGcfM7O3GEXkHdQK/AFhf3Scln3erm8KEiNzfXev6Dq6truWN69eolSfrwww+tJ3BXdM5EZedcuLm5yc3NrcxyFxcX/kE3cPQAbFGX++TCuupynQ0F7wFsRa/AFiV9UlO94vChws3NTePGjSt33bZt25SZmamhQ4eqdevW8vf3V0BAgNq2bav09HTl5eWV+k1cXl6e0tPT1a5dO/n6+tbWSwAAAAAcmsOHCnd39wrvKzFmzBhlZmZq5syZ+sMf/mBdHhERoTlz5mju3LlauHChdfncuXOVm5urWbNm1XjdAAAAQH3h8KHickyfPl3r16/XokWLtGfPHgUFBWn37t1KTk5Wjx49NGnSJHuXCAC17vCmTZovyV/S5/36qcPzz8t/4EA7VwUAcAT15j4VVeHh4aGtW7dq0qRJ2rdvn1544QV9++23mjp1qjZt2iR3d3d7lwgAtWr7I48ocNgwTZN0r6SwnTvlGx6utIgIe5cGAHAA9TpUJCYmyjCMUoc+lfD29lZsbKwOHz6swsJCHTp0SEuWLFHTpk3tUCkA2M/BlBTdkpCgRjo/fV3yp5OkXsuX69CmTXatDwBQ99XrUAEAqNzB556TUc5yiyRD0oHZs2u5IgCAoyFUAEAD53r0qCwVrLP8th4AgEshVABAA1fYtm25MxXS+ZmKwrZta7McAIADIlQAQAPnP2eO9VCnCxk6P1PR7vnna78oAIBDIVQAQAPnP3Cgdowbp2JJRZLO/vZVLGnHuHHy69/frvUBAOo+QgUAQKHx8cr8+GMtkfSepK09e+o/GzcqtIKbiwIAcKEGefM7AEBZvn37atZv3+du3iwPDw+71gMAcBzMVAAAAAAwhVABAAAAwBRCBQAAAABTCBUAAAAATCFUAAAAADCFUAEAAADAFEIFAAAAAFMIFQAAAABMIVQAAAAAMIVQAQAAAMAUQgUAAAAAU5ztXQAAwL6ysrKUlZWl/Px867KMjAy5u7tLknx8fOTj42Ov8gAADoCZCgBo4OLi4hQcHKzQ0FDrstDQUAUHBys4OFhxcXF2rA4A4AiYqQCABi4yMlJDhw6tcD2zFACAyhAqAKCB4/AmAIBZHP4EAAAAwBRCBQAAAABTCBUAAAAATCFUAAAAADCFUAEAAADAFEIFAAAAAFMIFQAAAABMIVQAAAAAMIVQAQAAAMAUQgUAAAAAUwgVAAAAAEwhVAAAAAAwxdneBQAAUNdkZWUpKyurwvU+Pj7y8fGpxYoAoG4jVAAAcJG4uDjFxMRUuD4qKkrR0dG1VxAA1HGECgAALhIZGamhQ4cqPz9foaGhkqS0tDS5u7tLErMUAHARQgUAABcpObwpLy/Puqx79+7y8PCwY1UAUHdxojYAAAAAUwgVAAAAAEwhVAAAAAAwhVABAAAAwBRCBQAAAABTCBUAAAAATCFUAAAAADCFUAEAQAUOb9qk+ZJWS/q8Xz8dTEmxd0kAUCc5fKj473//qxdffFHh4eG65ppr5OrqqjZt2mjEiBH65z//We5zcnJyNGXKFPn5+cnNzU3+/v6aNm2acnNza7l6AEBdtf2RRxQ4bJimSbpXUtjOnfIND1daRIS9SwOAOsfhQ8XLL7+syZMna//+/QoPD9fUqVMVGhqq9evX65ZbbtE777xTanxeXp7CwsIUGxurTp06afLkyerYsaOWLFmifv366cyZM3Z6JQCAuuJgSopuSUhQI0nOkvVPJ0m9li/XoU2b7FofANQ1zvYuwKyePXtqy5YtCgsLK7V8+/bt6t+/vx577DENHz5cbm5ukqTFixcrIyNDM2bM0MKFC63jn376aS1atEixsbGaOXNmrb4GAEDdcvC553R1OcstkgxJB2bPll///rVcFQDUXQ4/U3HXXXeVCRSSdOutt6pv3746ceKE9u7dK0kyDEPx8fHy9PTUs88+W2r8s88+K09PT8XHx9dK3QCAusv16FFZKlhn+W09AOB3Dh8qLsXFxUWS5Ox8fkImMzNTR48eVUhIiDw8PEqN9fDwUEhIiPbv368jR47Ueq0AgLqjsG1bGRWsM35bDwD4ncMf/lSRw4cPa+PGjfLx8VHXrl0lnQ8VkhQQEFDucwICApSUlKTMzEz5+vqWO6agoEAFBQXWxzk5OZKkoqIiFRUVVedLgIMoed95/3Ep9Iljufq552QZPFiGVGrGouTx1VFRNfZe0iuwFb0CW1zcJzXVL/UyVBQVFenhhx9WQUGBFi1apEaNGkmSsrOzJUne3t7lPs/Ly6vUuPIsWLBAMTExZZYnJyerSZMmZkuHA0vhUpOwAX3iOD4fMEAjN25UsX4PFhZJbw8YoKYFBfp6w4Ya3T+9AlvRK7BFSZ+cPn26RrZf70JFcXGxxowZo23btmn8+PF6+OGHq3X7M2fO1JQpU6yPc3Jy5Ovrq/DwcGsoQcNSVFSklJQUDRw40HrIHXAx+sQBDR6sbzds0EfDh8tfUuubbpL/3Lm6r4ZP0KZXYCt6Bba4uE9KjrKpbvUqVBQXF+uRRx7R6tWr9dBDD+m1114rtb5khqKimYiSv+SKZjIkyc3NzXolqQu5uLjwD7qBowdgC/rEsfgPGKBZv32fu2VLmfPxahK9AlvRK7BFSZ/UVK/Um1BRXFyssWPH6o033tDIkSOVmJgoJ6fS56GXnEtRcm7FxSo75wIAAABAWfXi6k8XBor77rtPb775pvU8igsFBASobdu2Sk9PV15eXql1eXl5Sk9PV7t27So8SRsAAABAWQ4fKkoOeXrjjTd0zz33aNWqVeUGCkmyWCyKiIhQbm6u5s6dW2rd3LlzlZubq/Hjx9dG2QAAAEC94fCHP82ZM0crV66Up6enAgMD9fzzz5cZM3z4cHXv3l2SNH36dK1fv16LFi3Snj17FBQUpN27dys5OVk9evTQpEmTavcFAAAAAA7O4UPFwYMHJUm5ubmaN29euWP8/f2tocLDw0Nbt25VdHS01q5dq9TUVPn4+Gjq1KmKioqSu7t7LVUOAAAA1A8OHyoSExOVmJhYped4e3srNjZWsbGxNVMUAAAA0IA4/DkVAAAAAOzL4WcqAACobllZWcrKylJ+fr51WUZGhvUQWR8fH/n4+NirPACoc5ipAADgInFxcQoODlZoaKh1WWhoqIKDgxUcHKy4uDg7VgcAdQ8zFQAAXCQyMlJDhw6tcD2zFABQGqECAICLcHgTAFQNoQIAAFSq5DyTihDEgIaNUAEAACoVFxenmJiYCtdHRUUpOjq69goCUKcQKgAAQKVKzjPJz8+3nsCelpZW6opYABouQgUAAKhUyeFNeXl51mXdu3eXh4eHHasCUFdwSVkAAAAAphAqAAAAAJhCqAAAAABgCqECAAAAgCmECgAAYLPDmzZpvqTVkj7v108HU1LsXRKAOoBQAQAAbLL9kUcUOGyYpkm6V1LYzp3yDQ9XWkSEvUsDYGeECgAAUKmDKSm6JSFBjXT+evQlfzpJ6rV8uQ5t2mTX+gDYF6ECAABU6uBzz8koZ7lFkiHpwOzZtVwRgLqEUAEAACrlevSoLBWss/y2HkDDRagAAACVKmzbttyZCun8TEVh27a1WQ6AOoZQAQAAKuU/Z471UKcLGTo/U9Hu+edrvygAdQahAgAAVMp/4EDtGDdOxZKKJJ397atY0o5x4+TXv79d6wNgX4QKAABgk9D4eGV+/LGWSHpP0taePfWfjRsVGh9v79IA2JmzvQsAAACOw7dvX8367fvczZvl4eFh13oA1A3MVAAAAAAwhVABAAAAwBQOfwIAALhMWVlZysrKqnC9j4+PfHx8arEiwD4IFQAAoFIlH57z8/OtyzIyMuTu7i6p4X54jouLU0xMTIXro6KiFB0dXXsFAXZCqAAAAJUq78NzaGio9fuG+uE5MjJSQ4cOVX5+vvXvIy0trVTYAhoCQgUAAKhUyYfnijTUD88lMzR5eXnWZd27d+eqWGhwCBUAAKBSDfXwJlsd3rRJ8yX5S/q8Xz91eP55+Q8caOeqgNrD1Z8AAABM2P7IIwocNkzTJN0rKWznTvmGhystIsLepQG1hlABAABwmQ6mpOiWhAQ10vnDP0r+dJLUa/lyHdq0ya71AbWFUAEAAHCZDj73nIxyllskGZIOzJ5dyxUB9kGoAAAAuEyuR4/KUsE6y2/rgYaAUAEAAHCZCtu2LXemQjo/U1HYtm1tlgPYDaECAADgMvnPmWM91OlChs7PVLR7/vnaLwqwAy4pCwAAcJn8Bw5U2rhx6rV8uYol66FQFkk7xo1TaP/+dqzOPkruvl4RLk9cPxEqAAAATAiNj9e3w4frwzvukL+kK3r21LXz5zfIQCGVf/f1CzXUu6/Xd4QKAAAAk3z79tWs377P3by5Qd9Ru+Tu6/n5+QoNDZUkpaWlyd3dXVLDvft6fcc5FQAAAKg2Pj4+CgoKUotfftF8SaslnZkyRS1++UVBQUGEinqKmQoAAIDLVHL+QH5+vnVZRkZGqd/KN8QP0dsfeUS3JCRomn67Z8fOnbKEhytt3DiFxsfbuzzUAGYqAAAALlNcXJyCg4Oth/lIUmhoqIKDgxUcHKy4uDg7Vmcf3GW8YWKmAgAA4DKVnD9QkYY4S3Hwued0dTnLL7zLuF8DPYm9PiNUAAAAXKaGenjTpXCX8bJKDpNzO3xYLdevl2tWlgp9fPTLsGEquOaaetFHhAoAAABUm8K2bWUcPlzuuoZ6l/G4uDgdiolRvH6/MaIhqXViosZJ8q8Hl9nlnAoAAABUG+4yXtbEgQO13GIp9zyTFRaLJoaH27W+6tCgQ8UXX3yhwYMHq1mzZvLw8NAf/vAHvfvuu/YuCwAAwGH5DxyoHePGqVhSkaSzv30V6/xdxhvi+RTOb7yhYuPimHU+ZBUbhlzeeKP2i6pmDfbwp9TUVA0aNEiNGzfW/fffr6ZNm2rt2rW67777dOTIEU2dOtXeJQIAADicrKwsNZk4UZ917ap/TZokf0meXbrI/Ykn1KJnT2VlZTn8+QNVdWT7dl1/ifWHt21T81qrpmY0yFBx9uxZjR8/Xk5OTtq2bZu6d+8uSXruuefUs2dPzZo1S3fffbf8/PzsWygAAICDiYuLU0xMTOmFX38tPfaYJCmqHpw/UFUHpQpDheW39TfUVjE1pEGGis2bN+vHH3/U2LFjrYFCkry9vTVr1iyNGTNGK1eu1HPPPWe/IgEAABwQl9ktK/PWW6V9+6znlZQoOSAqs3dvO1RVvRpkqNiyZYskKbyck2IGDRokSdq6dWttlgQAAFAv1IfLo1a3XqNGafyyZYo3DJ3V+WDh5OQkFRdrvMWiiFGj7F2iaQ0yVGRmZkqSAgICyqxr06aNPD09rWMuVlBQoIKCAuvjnJwcSVJRUZGKiopqoFrUdSXvO+8/LoU+ga3oFdiKXnEcn332mRIMQ9uk85eQlXSwuFjLJf1oGLrq00/Vo0ePGtn3xX1SU/1iMYxyTkWv58LDw5WSkqLMzExde+21ZdZfddVVys3NVXZ2dpl10dHRZY8TlLR69Wo1adKkRuoFAACA4/r111914sQJ/fzzz/rnP/+pX3/9VS1atNDNN9+sVq1aqXnz5mrRokWt1HL69Gk98MADys7OlpeXV7Vtt0HOVJgxc+ZMTZkyxfo4JydHvr6+Cg8Pr9Y3Bo6jqKhIKSkpGjhwoFxcXOxdDuoo+gS2oldgK3oFtri4T0qOsqluDTJUeHt7S1K5MxHS+aDQvHn5F/Zyc3OTm5tbmeUuLi78g27g6AHYgj6BregV2IpegS1K+qSmeqVB3vyu5FyK8s6bOHbsmHJzc8s93wIAAABAWQ0yVISFhUmSkpOTy6xLSkoqNQYAAADApTXIUNG/f3+1b99eq1evVkZGhnV5dna25s+fL1dXV42qB5f2AgAAAGpDgzynwtnZWfHx8Ro0aJB69+6t+++/X02bNtXatWt16NAhLVmyRP7+/vYuEwAAAHAIDTJUSFLfvn2VlpamqKgovfPOOyoqKlLXrl21aNEi3XffffYuDwAAAHAYDTZUSFLPnj316aef2rsMAAAAwKE1yHMqAAAAAFQfQgUAAAAAUwgVAAAAAEwhVAAAAAAwhVABAAAAwBRCBQAAAABTCBUAAAAATCFUAAAAADCFUAEAAADAFEIFAAAAAFMIFQAAAABMIVQAAAAAMIVQAQAAAMAUZ3sX4OgMw5Ak5eTk2LkS2EtRUZFOnz6tnJwcubi42Lsc1FH0CWxFr8BW9ApscXGflHxmLfkMW10IFSadOnVKkuTr62vnSgAAAADbnDp1St7e3tW2PYtR3TGlgSkuLtbRo0fVtGlTWSwWe5cDO8jJyZGvr6+OHDkiLy8ve5eDOoo+ga3oFdiKXoEtLu4TwzB06tQptW3bVk5O1XcmBDMVJjk5Oenqq6+2dxmoA7y8vPihjkrRJ7AVvQJb0SuwxYV9Up0zFCU4URsAAACAKYQKAAAAAKYQKgCT3NzcFBUVJTc3N3uXgjqMPoGt6BXYil6BLWqrTzhRGwAAAIApzFQAAAAAMIVQAQAAAMAUQgUAAAAAUwgVAAAAAEwhVAAX+eKLLzR48GA1a9ZMHh4e+sMf/qB3333XpucahqFPP/1Ujz32mLp16yZvb281adJEN9xwg+bPn68zZ87UcPWoTWZ6pTwnTpzQVVddJYvFottuu60aK4W9VVev/O9//9PkyZMVEBCgxo0bq2XLlurVq5deffXVGqgata06+uTo0aN66qmndN1118nDw0NXXnmlQkND9eabb+rcuXM1VDlq06pVqxQZGambbrpJbm5uslgsSkxMrPJ2iouL9fLLL6tr165yd3dX69atNXLkSO3fv/+y6uLqT8AFUlNTNWjQIDVu3Fj333+/mjZtqrVr1+rQoUNasmSJpk6desnnnzlzRu7u7nJzc1OfPn3UtWtXnTlzRklJScrMzFSPHj20ZcsWNWnSpJZeEWqK2V4pz4MPPqj169crLy9PgwYN0meffVYDlaO2VVevZGRkKDw8XCdOnNCQIUPUuXNn5ebmat++fXJ1ddWGDRtq+JWgJlVHn+zfv18333yzfvnlFw0aNEjdunVTTk6O1q1bp2PHjmnMmDFKSEiohVeDmuTv769Dhw6pVatW8vDw0KFDh5SQkKAxY8ZUaTvjx49XfHy8unTpoiFDhujo0aN699135enpqc8//1wBAQFVK8wAYBiGYRQVFRkdOnQw3NzcjD179liXnzx50ggMDDRcXV2NgwcPXnIbhYWFxvPPP2/8+uuvZZbfcccdhiRj8eLFNVE+alF19MrF3n//fUOS8X//93+GJGPQoEHVXDXsobp6JTs727jmmmuM1q1bG1999VW5+4Hjqq4+eeyxxwxJxosvvlhq+YkTJ4xrrrnGkFTln02oe1JSUqzv44IFCwxJRkJCQpW2sXnzZkOS0bt3b6OgoMC6fMOGDYYkIzw8vMp1cfgT8JvNmzfrxx9/1AMPPKDu3btbl3t7e2vWrFkqLCzUypUrL7kNFxcXPfPMM2revHmZ5TNnzpQkbd26tdprR+2qjl650PHjx/XYY4/p4Ycf1pAhQ2qgYthLdfXKK6+8osOHD2vhwoXq1q1bmfXOzs7VWTZqWXX1SclhK4MHDy61vFmzZgoNDZUk/fzzz9VXOOxiwIAB8vPzM7WNZcuWSZLmzp0rV1dX6/I//vGP6tOnj5KTk3X48OEqbZNQAfxmy5YtkqTw8PAy6wYNGiTJXCBwcXGRxH/+9UF198qECRPUqFEjvfTSS9VSH+qO6uqVd955RxaLRSNGjNB3332nl19+WYsXL9ZHH32kwsLCaq0Zta+6+uT666+XpDKHwp08eVLp6elq06aNrrvuOpPVoj7YsmWLPDw8FBISUmbd5X7m4dMN8JvMzExJKvcYwjZt2sjT09M65nKsWLFCUvn/acCxVGevrFq1Sh988IHWrVun5s2bKzs7u1prhX1VR68UFhZq7969at26tV5++WVFRUWpuLjYur59+/Zat26dunbtWr3Fo9ZU18+UadOm6eOPP9bkyZP12WeflTqnokmTJvrwww/l7u5e7fXDseTl5SkrK0vXX3+9GjVqVGZ9SR9W9TMPMxXAb0o+zHl7e5e73svL67I/8H366aeKi4tT586dNW7cuMuuEXVDdfXK0aNH9eSTT2rkyJEaNmxYtdaIuqE6euXXX3/VuXPn9Msvv2jOnDlavHixfvrpJ/3nP//Rs88+qwMHDuiOO+7g6nIOrLp+plx55ZXasWOHbrvtNn322WdavHixXnvtNWVnZ2vUqFG64YYbqrVuOCZb+u3CcbYiVAA17IsvvtB9990nb29vvffee3Jzc7N3SagjIiIi5OLior/+9a/2LgV1WMmsxLlz5zRx4kRNnTpVV1xxha666irNmTNH99xzjw4dOqT333/fzpXC3n744QeFhITo+PHj2r59u06dOqUjR47oueee09y5c9W/f38uK4saQ6gAflOS2CtK5jk5ORWm+op8+eWXCg8Pl5OTk5KSktSlSxfTdcL+qqNXVq5cqU8//VR/+9vf1KpVq2qvEXVDdfTKheuHDh1aZn3Jsi+//PJyy4SdVdf/P2PGjNGhQ4f08ccfKzQ0VJ6enrr66qv19NNP609/+pN27NihNWvWVGvtcDy29NuF42xFqAB+c6ljCI8dO6bc3NwqXbP5yy+/1MCBA1VcXKykpCT16NGj2mqFfVVHr+zZs0eSdM8998hisVi/2rVrJ0lKSkqSxWIpdSUYOJ7q6BUPDw9dddVVks5fxediJcvy8/PNFQu7qY4+OXXqlNLT09W5c2e1adOmzPq+fftK+v1nDxouDw8P+fj46MCBA+XOXF3qHJ9LIVQAvwkLC5MkJScnl1mXlJRUakxlSgLFuXPn9Nlnn+nmm2+uvkJhd9XRK7169dK4cePKfN13332SpKuvvlrjxo3TXXfdVc3VozZV18+Vfv36SZK++eabMutKlvn7+19umbCz6uiTkquAVXTJ2OPHj0sSh+BC0vl+ysvLU3p6epl1JT3Xu3fvqm20yne2AOqpoqIio3379pe8+dCBAwesy48ePWrs27fPOHnyZKntfPnll0azZs0MT09PIy0trZaqR22qrl4pz4EDB7j5XT1SXb2Snp5uSDK6dOlinDhxwro8KyvLuOqqqwwnJyfju+++q+FXg5pSXX3SsWNHQ5KxbNmyUstPnDhhdOrUyZBkpKSk1ORLQS2r7OZ3x48fN/bt22ccP3681PKauPkdoQK4wObNmw0XFxejadOmxvjx440pU6YYfn5+hiRjyZIlpcaOHj26zD/kX375xWjevLkhybjtttuMqKioMl+xsbG1+6JQI8z2SkUIFfVPdfXKlClTDEmGr6+vMXHiRGP8+PHGFVdcYUgy5s+fX0uvBjWlOvpkw4YNhrOzsyHJ6N+/v/HnP//ZGDdunNG6dWtDkjFixIhafEWoKcuWLTNGjx5tjB492ggKCjIkGSEhIdZlF4bKqKgoQ5IRFRVVZjsRERHWX1ZMnz7dePjhhw1XV1ejRYsWl/VLCkIFcJF//vOfxm233WZ4eXkZ7u7uRs+ePY01a9aUGVfeD/WSD4SX+vLz86u9F4MaZaZXKkKoqJ+qq1cSEhKMm266yWjSpInh4eFhhIaGGh988EENV4/aUh19snPnTuOee+4xfHx8DGdnZ8PT09Po0aOH8fLLLxtnz56thVeBmlby/lf0NXr0aOvYS4WKc+fOGS+99JLRpUsXw83NzWjZsqVx3333GT/88MNl1WUxDMOo2gFTAAAAAPA7TtQGAAAAYAqhAgAAAIAphAoAAAAAphAqAAAAAJhCqAAAAABgCqECAAAAgCmECgAAAACmECoAAAAAmEKoAAAAAGAKoQIAAACAKYQKAAAAAKYQKgAAAACYQqgAANQZ119/vSwWS4VfMTEx9i4RAFAOZ3sXAABAiZEjR6qoqKjUssLCQr344ovKz89X79697VQZAOBSLIZhGPYuAgCA8hQUFOiuu+7Sp59+qldeeUUTJkywd0kAgHIwUwEAqJPy8/M1fPhwbdy4UcuWLdO4cePsXRIAoAKECgBAnXP69GkNHTpUqampSkhI0KhRo+xdEgDgEggVAIA6JS8vT0OGDFFaWprefPNNPfDAA/YuCQBQCUIFAKDOOHXqlAYPHqzPP/9cb7/9tu655x57lwQAsAGhAgBQJ2RnZ+u2227Trl279N5772n48OH2LgkAYCNCBQDA7k6cOKHw8HDt3btXH3zwgW6//XZ7lwQAqAIuKQsAsLshQ4Zow4YNCgsLU58+fcqsv+KKKzRx4sTaLwwAYBNCBQDAroqLi+Xt7a3c3NwKxwwbNkzr1q2rvaIAAFVCqAAAAABgipO9CwAAAADg2AgVAAAAAEwhVAAAAAAwhVABAAAAwBRCBQAAAABTCBUAAAAATCFUAAAAADCFUAEAAADAFEIFAAAAAFMIFQAAAABMIVQAAAAAMIVQAQAAAMAUQgUAAAAAU/4fioDVzQ24FCsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,6))\n",
    "\n",
    "\n",
    "plt.errorbar(z_center,N_clusters_z_theory,yerr=np.sqrt(N_clusters_z_theory),\n",
    "            marker='o',\n",
    "            markersize=5.,\n",
    "            fillstyle='none',\n",
    "            ls='none',\n",
    "            label=r'$\\mathrm{CLASS\\_SZ\\,\\,binned}$',capsize=3,c='k')\n",
    "plt.errorbar(np.asarray(z_center),Nz_binned,#yerr=np.sqrt(Nz_binned),\n",
    "             marker='o',markersize=5.,ls='none',\n",
    "            label=r'$\\mathrm{CLASS\\_SZ\\,\\,from\\,\\,unbinned}$',capsize=5,c='red')\n",
    "\n",
    "\n",
    "\n",
    "plt.xlabel('$z$', fontsize=14)\n",
    "plt.ylabel('$N$', fontsize=14)\n",
    "\n",
    "plt.xticks(fontsize=14)\n",
    "plt.yticks(fontsize=14)\n",
    "plt.legend(fontsize=13,frameon=False)\n",
    "plt.grid()\n",
    "plt.tight_layout()\n",
    "# plt.loglog()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "0b5e60df-ec3d-4ee3-bb41-fb8bcf94d2de",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyYAAAJQCAYAAAB7O3BaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABXH0lEQVR4nO3deVxWZf7/8fctAiKbGk7iiKAGaoqaW9No4oql5TI2Vk4qmkqL02KLY18LzUmzHG2mlcLdnGzStGksQAVcsslEGn9Zhrmkiaam3LKEt3J+fzjc0x0gIMvF8no+Hjzkvs51X+dz4M7O2+tc59gsy7IEAAAAAAbVM10AAAAAABBMAAAAABhHMAEAAABgHMEEAAAAgHEEEwAAAADGEUwAAAAAGEcwAQAAAGBcfdMF1HT5+fk6fvy4fH19ZbPZTJcDAAAAVBjLsnT+/Hk1b95c9epV7pwGwaScjh8/rqCgINNlAAAAAJXm6NGjatGiRaXug2BSTr6+vpIu/7L8/PwMVwOUj8PhUEJCgiIjI+Xu7m66HAAAYJjdbldQUJDznLcyEUzKqeDyLT8/P4IJajyHw6GGDRvKz8+PYAIAAJyqYskCi98BAAAAGEcwAQAAAGAcwQQAAACAcQQTAAAAAMYRTAAAAAAYRzABAAAAYBzBBAAAAIBxBBMAAAAAxhFMAAAAABhHMAEAAABgHMEEAAAAgHEEEwAAAADG1TddAAAAuHoZGRnKyMgodntgYKACAwOrsCIAuDoEEwAAarDY2FjNnj272O0xMTGaNWtW1RUEAFeJYAIAQA0WHR2tYcOGKTc3V71795Ykbd++XV5eXpLEbAmAGoNgAgBADVZwqVZ2drazrUuXLvL29jZYFQCUHYvfAQAAaplly5bJZrMpOTm5xL7Jycmy2WxatmxZpddVkUJCQtS3b1/TZaACEUwAAEClycnJ0UsvvaSbb75ZTZo0kbu7u6699loNGTJEy5Yt08WLFyX97+R4wYIFpR777Nmz8vLyks1m08qVK4vtd/DgQU2ZMkXt2rVTw4YN1bhxY7Vv317jx49XUlJSufuXpCzj2Wy2Un2VJnAANQ2XcgEAUAscOHDA+X1MTIyio6MVGhpqsKLLNQ0dOlTffPONBg4cqBkzZiggIEA//PCDNm3apAkTJmjfvn164YUXrmr8t99+W3l5eWrVqpWWLFmisWPHFurz+eefKyIiQu7u7ho3bpw6dOig3NxcpaenKyEhQb6+vurXr99V9y9JWccrKWDFxMQoICBAbdu2LXUNJenTp49yc3Pl7u5eYWMCV4NgAgBADbd06VJNnjzZ5fVLL72kuLg4RUVFGakpNzdXt912mw4ePKi1a9fqd7/7ncv26dOna9euXdq1a9dV72Px4sXq16+fhg8frkceeUQHDx5U69atXfrMnj1bOTk5SktLU+fOnQuNceLEiXL1L0lZx7vnnnuKHCcnJ0c33XST3NzctGbNmgq9qUG9evXUoEGDChsPuFpcylUDZWRkKDU1tdivK93PHgBQu6Snp2vy5MkuswX79+/XhAkTNGnSJJeZlKoUFxen/fv367HHHisUSgr06NFDDzzwwFWNn5qaqrS0NI0fP15jxoxR/fr1tWTJkkL90tPTdc011xQZCiSpWbNm5epfkooab+LEifrPf/6j+fPnq3///qXe/8WLFzVr1iwFBwfL09NTnTp10jvvvOPSp6g1JgVrVLZs2aIFCxaoTZs28vT0VFhYmJYvX15oP2Xtn5eXp7lz56pDhw5q0KCBGjVqpNtvv1179uwp1Pfo0aMaPXq0/P395efnp9tvv13ffvttqX8GqDkIJjVQbGysunXrVuxXbGys6RIBAFVkyZIl8vf3d1mb4eXlpb/97W/y8/PT4sWLjdT13nvvSZKmTJlSKeMvXrxYPj4+GjVqlAICAnTbbbdp+fLlys/Pd+nXpk0bnTlzRuvWrSvVuGXtXxXjvfjii1qzZo3uvPNOPfbYY2V67/Tp0/XOO+/ogQce0LPPPqsLFy7o7rvvLvVC96eeekorV65UdHS0XnjhBdWrV09RUVHasWPHVfd3OBy65ZZbNHv2bN10001atGiR/vSnP2nfvn3q1auXPv/8c2ffc+fOqU+fPlq3bp3Gjh2r559/Xg0bNlS/fv1c7kSHWsJCuWRmZlqSrMzMzCrb5/Hjx63du3db27dvtyRZkqzt27dbu3fvtnbv3m0dP368ympB7XLhwgVr/fr11oULF0yXAqCU7rrrLqtfv35WVlaW8/8JWVlZlmVZVr9+/ay77rrLSF1NmjSx/Pz8St0/KSnJkmS9+OKLJfbNzc21GjVqZI0fP97Ztn79ekuStXHjRpe+n3zyieXu7m5JskJDQ60JEyZYr732mrVv374ixy5r/5KUd7zExETLzc3NCg8Pd/5eS2Pp0qWWJKtly5bWuXPnnO3nzp2zWrZsaTVu3NjKycmxLOt/P/ulS5cWen+XLl2svLw8Z/uxY8csDw+PQp+rsvRfuHChJcn6+OOPXcbIzMy0goKCrIiICGfbjBkzLEnWkiVLXPo+/PDDliSXvqgcVXmuy4xJDRQYGKiuXbuqS5cuzrYuXbqoa9eu6tq1Kw/TAoA6JCQkRF988YVyc3Nd2nNzc5WWlqaQkBAjddntdvn6+lbK2OvWrdO5c+c0fvx4Z9uQIUPUtGnTQpdz3XTTTdq9e7fGjx+vzMxMLV26VA888ICuv/569enTRwcPHixX/5KUZ7zDhw/rrrvukq+vr95///2rejbN/fffL39/f+drf39/3XfffTp79myp7uz1wAMPyMPDw/n617/+tcLCwpSenn7V/VetWqV27dqpW7duOn36tPPrwoULGjRokLZv3+78PK9fv17XXnutxo0b57Kf6dOnl+r4UbMQTAAAqMEmTpyozMxMPfHEE8623NxcPfTQQ7Lb7br33nuN1OXn56fz589XytiLFy9W06ZN1aJFCx04cEAHDhzQkSNHFBkZqQ8++ECnT5926R8eHq5ly5bp5MmTOnz4sJYvX66bb75Z27Zt0/Dhw3XhwoVy9S/J1YyXk5OjkSNH6uzZs1q9erXatGlT9h+UpPbt2xdqu/766yWpVCHrlzcTkKRrrrlGZ86cuer+X331lb7++ms1bdq00NeSJUt06dIl5+/w4MGDCg0NlZubm8uYgYGBatSoUYn1o2bhrlwAANRgoaGhiouL06RJk5xtYWFhstvtiouL03XXXWekro4dO2rr1q1F3imrPA4dOqSkpCRZlqWwsLAi+6xatUqPPPJIkduCg4M1btw4jR07VjfffLN27Nihzz77TL17966Q/iUp7XiTJ09WWlqa5syZo1tvvfWq9lURfhkICliWddX9LctSeHi4Fi5cWOx+mzZtWoYqUVsQTAAAqOGioqLUtWtX552fJk6cqPvuu89YKJGkUaNGaevWrYqLi9PcuXMrbNylS5fKsiy99dZbRf6L+cyZM7VkyZJig0kBm82mG2+8UTt27ND3339f4n7L2r884y1cuFCrV6/W8OHD9X//93/l2s9XX32l4cOHu7Tt27dPUtGzG1UhNDRUp06dUv/+/VWv3pUv3mndurXS09N16dIll9CTkZGhc+fOVXKlqGpcygUAQC3w80t9Zs+ebTSUSNKkSZPUtm1bLViwQBs2bCiyz+7du/Xaa6+Vesz8/HwtW7ZM4eHhmjRpku64445CX3fffbf27t3rfD5KYmKi8+nyP5ebm6uEhARJ/7u06Wr6l6Ss4yUlJenJJ59U27ZttWLFCtlstlLvqyivv/66MjMzna8zMzP1xhtvqFGjRoqIiCjX2Fdr3LhxOnHiRLEzJidPnnR+P3z4cJ08eVIrVqxw6TN//vxKrRFmMGMCAAAqXMOGDfXhhx9q6NChGjFihCIjIzVo0CBdc801OnXqlJKSkhQfH68nn3zS5X2bN2/WTz/9VGi8gIAAhYSE6OjRo1dcNzNq1CjNmjVLixcvVo8ePfToo4/qzJkzGjZsmMLDw9WwYUMdPXpUq1ev1jfffKNx48YpPDzc+f6y9i9JWcbLyMjQ6NGjdenSJY0aNUoffPBBseN26tRJnTp1KnH/AQEBuvHGGzVhwgRJl2ecvvvuO8XFxalhw4alPo6K9PDDDysxMVFPPPGEtmzZov79+8vPz0/fffedNm/erAYNGigpKUmS9OSTT2r16tWaPHmydu/erQ4dOig5OVk7d+5UQECAkfpReQgmNZjtwAHNlRQiyT0mRoqOlkJDDVcFAMBl1113nfbs2aPY2FitXbtWzz33nLKystSkSRN1795dy5cv15gxY1ze8/HHH+vjjz8uNFbbtm2dJ/DFPbBRury2JSwsTO+8844WLVqkhQsXasOGDdq+fbvWrl2rc+fOyd/fX506ddL06dMVFRXl8v6y9i9JWcbbv3+/c9F3SZe/xcTElCqYzJ8/X9u2bdOrr76qkydPKiwsTG+//Xahn3tVcnd317/+9S+99tprWrlypWJiYiRJzZs3V8+ePV3utta4cWNt27ZN06ZNc86aREREKCkpSQMGDDBSPyqPzSpu9RJKxW63y9/fX5mZmfLz86u6HS9dKmvSJF3Kz5dNUj03N9ksS1q8WCrjX5pAAYfDoY0bN2rIkCFyd3c3XQ6AUsjIyFBGRoZyc3OdC6i3b98uLy8vSZfvXsRt5AFcrao812WNSU2Uni5NmiRbfr7qS3KTZLt0ScrPl+69VzpwwHSFAIAqEhsbq27durnc1al3797q1q2bunXrptjYWIPVAUDpcSlXTbRkiVTcYjib7fKsybx5VVsTAMCI6OhoDRs2rNjtzJZUvEuXLunUqVMl9mvSpInLwwYBXBnBpCY6fFgq7go8y7q8HQBQJ3CpVtU7evSoWrVqVWK/pKQk9e3bt/ILAmoJgklNFBJy5RmTkJCqrAYAgDqlWbNmSkxMLLFfwXNlAJQOwaQmmjhReuEFWZIKxRPLurzOBAAAVIoGDRpo4MCBpssAah0Wv9dEoaGX15HUqyeHpIuSLDc3qV69y+2GH6oFAAAAlBXBpKaKilLunj1aIOkfkhyPPCLt38+tggEAAFAjcSlXDWa1aaOn/vt91uzZ8vD2NloPAAAAcLWYMQEAAABgHMEEAAAAgHEEEwAAAADGEUxqsAMHDji/j4mJUXp6usFqAAAAgKtHMKmhli5dqm7durm8bt++vZYtW2auKAAAAOAqEUxqoPT0dE2ePFljx451tu3fv18TJkzQpEmTXGZSAAAAgJqAYFIDLVmyRP7+/lqwYIGzzcvLS3/729/k5+enxYsXG6wOAAAAKDuCSQ10+PBhde7cWQ0aNHBp9/LyUpcuXXT48GEzhQEAgGrj0KFDGjFihJo2bSqbzaYoHsIsSVq2bJlsNpuSk5NL7JucnCybzVbjLpUPCQlR3759TZdRZgSTGigkJERffPGFcnNzXdpzc3OVlpamkJAQM4UBAPALOTk5eumll3TzzTerSZMmcnd317XXXqshQ4Zo2bJlunjxoqT/nQD+/GqAkpw9e1ZeXl6y2WxauXJlsf0OHjyoKVOmqF27dmrYsKEaN26s9u3ba/z48UpKSip3/5KUZTybzVaqr9KcVEdFRSklJUXTp0/XypUrFR0dXebagarEk99roIkTJ+rFF1/UE0884WzLzc3VjBkzZLfbde+99xqsDgCAyw4cOKChQ4fqm2++0cCBAzVjxgwFBATohx9+0KZNmzRhwgTt27dPL7zwwlWN//bbbysvL0+tWrXSkiVLXNZeFvj8888VEREhd3d3jRs3Th06dFBubq7S09OVkJAgX19f9evX76r7l6Ss45UUsGJiYhQQEKC2bdtecb95eXnatm2bpk6dqscff7zU9cJVnz59lJubK3d3d9Ol1AkEkxooNDRUcXFxmjRpkrMtLCxMdrtdcXFxuu666wxWBwAwIj1dWrJEOnxYCgmRJk6UQkONlZObm6vbbrtNBw8e1Nq1a/W73/3OZfv06dO1a9cu7dq166r3sXjxYvXr10/Dhw/XI488ooMHD6p169YufWbPnq2cnBylpaWpc+fOhcY4ceJEufqXpKzj3XPPPUWOk5OTo5tuuklubm5as2aNAgMDr7jfkydPyrIsNWnSpMQaL126pLy8PDVs2LDEvnVNvXr1Cl06j8rDpVw1VFRUlFJTU52vJ06cqK+//prrRwGgLlq6VGrXTnrxRenddy//2a6dZPC6+Li4OO3fv1+PPfZYoVBSoEePHnrggQeuavzU1FSlpaVp/PjxGjNmjOrXr68lS5YU6peenq5rrrmmyFAgSc2aNStX/5JU1HgTJ07Uf/7zH82fP1/9+/e/Yt+oqCgFBwdLuhyMfn75V8H6ik2bNmnOnDlq06aNGjRooHfffdf5/tOnT+vBBx9UUFCQPDw8FBQUpAcffFBnzpxx2U/BWJs3b9azzz6r4OBgeXl56cYbb9Snn34qSUpJSVHv3r3l7e2twMBAzZkzp1THO2vWLNlstiLXzf5y/URBHVu2bNGCBQvUpk0beXp6KiwsTMuXLy92HxcvXtSsWbMUHBwsT09PderUSe+8845Ln6LWmJR1f2Xtn5eXp7lz56pDhw5q0KCBGjVqpNtvv1179uwp1Pfo0aMaPXq0/P395efnp9tvv13ffvttscdc3TFjUoO1adPG+f3s2bPl7e1tsBoAgBHp6dKkSVJ+fuFt994r9e4tGZhJf++99yRJU6ZMqZTxFy9eLB8fH40aNUre3t667bbbtHz5cj377LOqV+9//+7apk0b7d+/X+vWrSs2IP1cWftXxXgvvvii1qxZozvvvFOPPfZYif2jo6PVpUsXPfrooxo5cqRzv+3bt3ee6D/++ONyOByaPHmy/Pz8nJeGZWZm6re//a0OHDigiRMnqmvXrtqzZ49ef/11bdmyRZ999pl8fX1d9venP/1Jly5d0sMPP6wLFy7oL3/5iyIjI7VixQrde++9mjJliv7whz/o3Xff1TPPPKNWrVoVOzNUHk899ZRyc3MVHR0tT09Pvf7664qKitJ1112nXr16Feo/ffp0ZWdnO8Px0qVLdffdd+unn34q1T/0lnV/penvcDh0yy236JNPPtHYsWM1depUZWZm6q233lKvXr20detWde/eXZJ07tw59enTR0ePHtV9992n66+/XikpKerXr1+hdcg1hoVyyczMtCRZmZmZVb7vrKwsS5IlycrKyqry/aP2uXDhgrV+/XrrwoULpksBUFp/+pNlublZllT4y83t8nYDmjRpYvn5+ZW6f1JSkiXJevHFF0vsm5ubazVq1MgaP368s239+vWWJGvjxo0ufT/55BPL3d3dkmSFhoZaEyZMsF577TVr3759RY5d1v4lKe94iYmJlpubmxUeHl6m/9cfOnTIkmTFxMS4tC9dutSSZIWFhVnZ2dmF3vfUU09ZkqxXX33Vpf2VV16xJFkzZ84sNNYNN9xg5eXlOds3bNhgSbLq169v7dq1y9mel5dnNWvWzPrNb35TYv0xMTGWJOvQoUOFtgUHB1sRERGF6ujSpYtLHceOHbM8PDysu+66q8ifQcuWLa1z584528+dO2e1bNnSaty4sZWTk2NZ1v8+l0uXLi33/krTf+HChZYk6+OPP3YZIzMz0woKCnI57hkzZliSrCVLlrj0ffjhhy1JLn3LoyrPdbmUCwCAmuzw4csxpCiWdXm7AXa7vdC/rFeUdevW6dy5cxo/fryzbciQIWratGmhy7luuukm7d69W+PHj1dmZqaWLl2qBx54QNdff7369OmjgwcPlqt/Scoz3uHDh3XXXXfJ19dX77//foVeGXH//fcXuabk/fffV9OmTQvNdEVHR6tp06Z6//33ixzLw8PD+frmm2+WJN14443Of92XJA8PD/Xs2VPp6ekVdRguHnjgAZc6fv3rXyssLKzY/d1///3y9/d3vvb399d9992ns2fPluquZ2XdX2n6r1q1Su3atVO3bt10+vRp59eFCxc0aNAgbd++3Tkbsn79el177bUaN26cy36mT59eYu3VFZdyAQBQk4WESDZb0dtstsvbDfDz89P58+crZezFixeradOmatGihQ4cOOBsj4yM1D/+8Q+dPn1aAQEBzvbw8HDnGoEjR44oJSVFcXFx2rZtm4YPH67du3e7nDCWtX9Jrma8nJwcjRw5UmfPntWHH37ocvl2RQgLCyuy/dChQ+revbvq13c9Raxfv77CwsJc1rcW+OUNBxo3bixJatWqVaG+jRs3LrRWpaL8sg5Juuaaa3TkyJEi+7dv375Q2/XXXy9JpQqgZd1fafp/9dVXys3NVdOmTYvd7+nTpxUUFKSDBw+qR48ecnNzc9keGBioRo0alVh/dUQwAQCgJps4USrudruWdXmdiQEdO3bU1q1bi7xTVnkcOnRISUlJsiyr2JPrVatW6ZFHHilyW3BwsMaNG6exY8fq5ptv1o4dO/TZZ5+pd+/eFdK/JKUdb/LkyUpLS9OcOXN06623XtW+rqQi78D1yxPjktpLw1Zc2Jacz74p7f6s4mYUy6ms+ytNf8uyFB4eroULFxa73yuFlpqOYAIAQE0WGiotXnw5gNhsl8NIwZ+LFxtZ+C5Jo0aN0tatWxUXF6e5c+dW2LhLly6VZVl66623ivxX4ZkzZ2rJkiXFBpMCNptNN954o3bs2KHvv/++xP2WtX95xlu4cKFWr16t4cOH6//+7//Kva+yaN26tfbv36+LFy+6zJpcvHhR33zzTYWGzCspuM3xjz/+6PLg6J9++kkZGRkV8miEr776SsOHD3dp27dvn6SiZzeqQmhoqE6dOqX+/fu73MShKK1bt1Z6erouXbrkEnoyMjJ07ty5Sq60crDGBACAmi4qStq/X3riCWn06Mt/7t9/ud2QSZMmqW3btlqwYIE2bNhQZJ/du3frtddeK/WY+fn5WrZsmcLDwzVp0iTdcccdhb7uvvtu7d271/l8lMTExCL/hT03N1cJCQmS/nf5ztX0L0lZx0tKStKTTz6ptm3basWKFVecOagMI0aM0KlTpxQXF+fS/tZbb+nUqVMaOXJkldRRMBu2adMml/ZFixYpv6g70F2F119/XZmZmc7XmZmZeuONN9SoUSNFRERUyD7Katy4cTpx4kSxMyYnT550fj98+HCdPHlSK1ascOkzf/78Sq2xMjFjAgBAbXDdddK8eaarcGrYsKE+/PBDDR06VCNGjFBkZKQGDRqka665RqdOnVJSUpLi4+P15JNPurxv8+bN+umnnwqNFxAQoJCQEB09elT3XuHytFGjRmnWrFlavHixevTooUcffVRnzpzRsGHDFB4eroYNG+ro0aNavXq1vvnmG40bN07h4eHO95e1f0nKMl5GRoZGjx6tS5cuadSoUfrggw+KHbdTp07q1KlTqesorSeffFL/+Mc/9OCDDyo1NVU33HCD9uzZo8WLF6tt27aFfl+VZeDAgWrbtq2eeeYZnTlzRq1atdL27dv16aefuqwfKo+AgADdeOONmjBhgqTLs3Hfffed4uLijD1s8uGHH1ZiYqKeeOIJbdmyRf3795efn5++++47bd68WQ0aNFBSUpKky7+r1atXa/Lkydq9e7c6dOig5ORk7dy5s8J+RlWNYAIAACrFddddpz179ig2NlZr167Vc889p6ysLDVp0kTdu3fX8uXLNWbMGJf3fPzxx/r4448LjdW2bVvnCfyVngfSsWNHhYWF6Z133tGiRYu0cOFCbdiwQdu3b9fatWt17tw5+fv7q1OnTpo+fXqh51WUtX9JyjLe/v37dfr0aUkq8fK3mJiYSgkm/v7+2rFjh2JiYvTBBx9o6dKluvbaa3Xfffdp9uzZlXantV9yc3PTBx98oIceekgvv/yyPDw8FBkZqZSUlCKfEXI15s+fr23btunVV1/VyZMnFRYWprfffrvQZ7Iqubu761//+pdee+01rVy5UjExMZKk5s2bq2fPni53omvcuLG2bdumadOmOWdNIiIilJSUpAEDBhipv7xsVmWtCLpKP/30k5566il9/vnnOnDggH788Uc1atRIbdq00aRJk3TPPffI3d3d5T12u12zZs3S2rVrdeLECQUGBur3v/+9YmJi5OPjU2gf+fn5evXVV/Xmm2/qwIED8vHx0cCBA/Xcc8+V+ZpCu90uf39/ZWZmys/Pr1zHXlbZ2dnO48vKyuIBiyg3h8OhjRs3asiQIYX+OwMAAHVPVZ7rVrs1JllZWXr99ddls9k0dOhQTZs2TSNHjtT333+viRMn6rbbbnO5tjA7O1sRERFatGiR2rVrp0cffdR5TWv//v2LnA6Ojo7WQw89JMuy9NBDD+mWW27RunXr1KNHj0q7tzYAAACA4lW7S7maNGmizMzMQvfzvnjxogYNGqSEhAR99NFHGjp0qCTphRdeUFpamqZPn67nn3/e2f9Pf/qT5s+fr0WLFmnGjBnO9qSkJMXFxalPnz5KTEx07mfMmDEaMmSIpk6dqvj4+Co4UgAAUBNdunRJp06dKrFfkyZNyvS8E6Cuq3YzJvXq1SvyP+L69es77wRR8DAly7IUFxcnHx8fPf300y79n376afn4+BR5VwlJmjNnjst+br31VvXt21cJCQn67rvvKvSYAABA7XH06FEFBgaW+PXJJ5+YLhWoUardjElx8vPznYvhOnbsKElKT0/X8ePHNXjw4ELrK7y9vdWrVy/Fx8fr6NGjCgoKkiQlJyc7t/3S4MGDlZycrJSUFI0dO7aSjwgAANREzZo1U2JiYon9OnfuXAXVALVHtQ0mFy5c0Ny5c2VZls6cOaPNmzfr66+/1oQJE5x3GihYDxIaGlrkGKGhoYqPj1d6erqCgoKUnZ2tjIwMdezYscinbxaMc6V1Jnl5ecrLy3O+ttvtki4vGnY4HFd3sFfp5/szsX/UPgWfIT5LAFA8Nze3Uj/ngr9PUdNV5We4WgeT2bNnO1/bbDY9/vjjmveze7QXPBTH39+/yDEK7hxQ0K+s/Ysyb948l7oKJCQkVPk9r3++sD8+Pl4NGjSo0v2j9irNvwQCAIDaLycnp8r2VW2DiY+PjyzLUn5+vo4fP65//vOfeuqpp7Rz505t3Lixym/NW2DGjBmaNm2a87XdbldQUJAiIyON3C64QFGXswFl5XA4lJiYqEGDBnG7YAAA4Lw6qCpU22BSoF69emrRooXuv/9+BQQEaPTo0Xruuec0f/5858xHcTMcBT/Ign5l7V8UT09PeXp6Fmp3d3ev8hO5n+/PxP5Re/F5AgAAkqr0fKDa3ZXrSiIjIyVdXsAulbwm5JdrULy9vRUYGKhDhw7p0qVLJfYHAAAAUDVqVDA5fvy4pP8lt9DQUDVv3lw7duxwuaxJunyZ044dO9SqVSvnHbkkKSIiwrntlwqeX9KnT5/KOgQAAAAARah2wWTfvn1FLrLJyclxru0YMmSIpMsL4idNmqSsrCzNmTPHpf+cOXOUlZWlyZMnu7RPmTJF0uXnnFy4cMHZ/tFHHyk5OVmRkZEKDg6u0GOqaBkZGUpNTVVaWpqzLS0tTampqUpNTVVGRoa54gAAAICrYLMsyzJdxM/NmjVLCxcuVO/evRUSEiI/Pz99//33+uijj3TmzBndfPPNio+Pl5eXl6TLMyO9evXSF198ocjISHXt2lWpqalKSEhQjx49lJKS4uxbYPLkyYqLi1OHDh00dOhQZWRkaM2aNfLx8dHOnTsVFhZW6nrtdrv8/f2VmZlZZYvfZ82aVeSdwQrExMRo1qxZVVILaheHw6GNGzdqyJAhrDEBAABVeq5b7YLJ559/rjfffFOffPKJvv/+e2VlZcnf31+dOnXSXXfdpYkTJ6p+fdc1+5mZmZo1a5bWrl2rEydOKDAwUL///e8VExMjX1/fQvvIz8/XK6+8ojfffFMHDhyQj4+PBg4cqOeee05t2rQpU70mgklGRsYVZ0UKnjgLlBXBBAAA/FydDiY1jYlgAlQWggkAAPi5qjzXrXZrTAAAAADUPQQTAAAAAMYRTAAAAAAYRzABAAAAYBzBBAAAAIBxBBMAAAAAxhFMAAAAABhHMAEAAABgHMEEAAAAgHEEEwAAAADGEUwAAAAAGEcwAQAAAGAcwQQAAACAcQQTAAAAAMYRTAAAAAAYRzABAAAAYBzBBAAAAIBxBBMAAAAAxhFMAAAAABhHMAEAAABgHMEEAAAAgHEEEwAAAADGEUwAAAAAGEcwAQAAAGAcwQQAAACAcQQTAAAAAMYRTAAAAAAYRzABAAAAYBzBBAAAAIBxBBMAAAAAxhFMAAAAABhHMAEAAABgHMEEAAAAgHEEEwAAAADGEUwAAAAAGEcwAQAAAGAcwQQAAACAcQQTAAAAAMYRTAAAAAAYRzABAAAAYBzBBAAAAIBxBBMAAAAAxhFMAAAAABhHMAEAAABgHMEEAAAAgHEEEwAAAADGEUwAAAAAGEcwAQAAAGAcwQQAAACAcQQTAAAAAMYRTAAAAAAYRzABAAAAYBzBBAAAAIBxBBMAAAAAxhFMAAAAABhHMAEAAABgHMEEAAAAgHHVLph8//33eumllxQZGamWLVvKw8NDzZo106hRo/Tvf/+7UP9Zs2bJZrMV+3X48OEi9xMfH6+IiAj5+vrKz89P/fr10+bNmyv56AAAAAAUpb7pAn7p5Zdf1vz589WmTRtFRkaqadOmSk9P1/r167V+/XqtXr1ad955Z6H3jR8/XiEhIYXaGzVqVKht1apVGjt2rJo2baqoqChJ0po1azRo0CC9++67uuOOOyr4qAAAAABcSbULJj179lRycrIiIiJc2rdt26YBAwbo/vvv14gRI+Tp6emyPSoqSn379i1x/LNnz+qPf/yjAgIClJqaqhYtWkiSpk+frhtuuEH333+/Bg8eLF9f3wo7JgAAAABXVu0u5frd735XKJRI0s0336x+/frp7Nmz2rt371WP/49//EPnzp3TH//4R2cokaQWLVpo6tSpOn36tN5///2rHh8AAABA2VW7YHIl7u7ukqT69QtP9GzdulXz58/Xiy++qPXr1ysrK6vIMZKTkyVJkZGRhbYNHjxYkpSSklJBFQMAAAAojWp3KVdxvvvuO23atEmBgYEKDw8vtD0mJsbldaNGjfTXv/5V48aNc2lPT0+XJIWGhhYao6CtoE9R8vLylJeX53xtt9slSQ6HQw6Ho5RHA1RPBZ9hPssAAECq2nOCGhFMHA6Hxo4dq7y8PM2fP19ubm7ObZ07d9aSJUvUt29fBQYG6sSJE/rwww/1zDPPKCoqSo0aNdKwYcOc/TMzMyVJ/v7+hfbj5+fn0qco8+bN0+zZswu1JyQkqGHDhld9jEB1kpiYaLoEAABQDeTk5FTZvmyWZVlVtrerkJ+fr7Fjx2r16tWaPHmy3nzzzVK9b/PmzRo0aJA6duyo//znP872sLAwpaeny+FwFLokzOFwyMPDQ506ddIXX3xR5LhFzZgEBQXp9OnTzmAD1FQOh0OJiYkaNGiQ89JJAABQd9ntdgUEBCgzM7PSz3Wr9YxJfn6+Jk6cqNWrV+uee+7RG2+8Uer3DhgwQG3atNHevXtlt9udP8iCmZLMzExdc801Lu8puCyrqNmUAp6enoXuCCZdXv/CiRxqCz7PAABAUpWeD1Tbxe/5+fmaMGGCli9frrvvvlvLli1TvXplKzcgIECS6xTUldaRXGn9CQAAAIDKUy2DSUEoWbFihe68806tXLnSZV1JaWRnZ+vLL7+Ut7e3M6BIct6KOCEhodB74uPjXfoAAAAAqBrVLpgUXL61YsUK/f73v9eqVauKDSXnz5/XN998U6g9NzdXkydP1vnz5zV69GiXtSSjR4+Wv7+/Xn75ZR07dszZfuzYMb3yyisKCAjQyJEjK/7AAAAAABSr2q0xefbZZ7V8+XL5+PgoLCxMf/7znwv1GTFihLp06aIzZ86oXbt26tGjh9q3b69mzZrp5MmT2rRpk44dO6bw8HC9+OKLLu9t3LixXnnlFY0dO1Zdu3bVnXfeKUlas2aNzpw5ozVr1vDUdwAAAKCKVbtgcvjwYUlSVlaWnnvuuSL7hISEqEuXLmrSpIkeeOABffbZZ9q4caPOnj0rLy8vtW/fXg899JCmTp0qLy+vQu+/5557FBAQoLlz52rp0qWy2Wzq1q2bZs6cqYEDB1bm4QEAAAAoQrW/XXB1Z7fb5e/vXyW3UAMqm8Ph0MaNGzVkyBDuygUAAKr0XLfarTEBAAAAUPcQTAAAAAAYRzABAAAAYBzBBAAAAIBxBBMAAAAAxhFMAAAAABhHMAEAAABgHMEEAAAAgHEEEwAAAADGEUwAAAAAGEcwAQAAAGAcwQQAAACAcQQTAAAAAMYRTAAAAAAYRzABAAAAYBzBBAAAAIBxBBMAAAAAxhFMAAAAABhHMAEAAABgHMEEAAAAgHEEEwAAAADGEUwAAAAAGEcwAQAAAGAcwQQAAACAcQQTAAAAAMYRTAAAAAAYRzABAAAAYBzBBAAAAIBxBBMAAAAAxhFMAAAAABhHMAEAAABgHMEEAAAAgHEEEwAAAADGEUwAAAAAGEcwAQAAAGAcwQQAAACAcQQTAAAAAMYRTAAAAAAYRzABAAAAYBzBBAAAAIBxBBMAAAAAxhFMAAAAABhHMAEAAABgHMEEAAAAgHEEEwAAAADGEUwAAAAAGEcwAQAAAGAcwQQAAACAcQQTAAAAAMYRTAAAAAAYRzABAAAAYBzBBAAAAIBxBBMAAAAAxhFMAAAAABhHMAEAAABgHMEEAAAAgHEEEwAAAADGVbtg8v333+ull15SZGSkWrZsKQ8PDzVr1kyjRo3Sv//97yLfY7fbNW3aNAUHB8vT01MhISF64oknlJWVVWT//Px8vfzyywoPD5eXl5eaNm2qu+++WwcPHqzMQwMAAABQjGoXTF5++WU9+uijOnjwoCIjI/XYY4+pd+/e2rBhg377299qzZo1Lv2zs7MVERGhRYsWqV27dnr00UfVtm1bLViwQP3799dPP/1UaB/R0dF66KGHZFmWHnroId1yyy1at26devToofT09Ko6VAAAAAD/Vd90Ab/Us2dPJScnKyIiwqV927ZtGjBggO6//36NGDFCnp6ekqQXXnhBaWlpmj59up5//nln/z/96U+aP3++Fi1apBkzZjjbk5KSFBcXpz59+igxMVEeHh6SpDFjxmjIkCGaOnWq4uPjq+BIAQAAABSwWZZlmS6itAYPHqyEhATt2rVL3bt3l2VZatGihex2u06cOCFvb29n3+zsbDVr1ky/+tWv9O233zrbx4wZo7///e9KSUlRnz59XMbv16+fkpOTdeTIEbVs2bJUNdntdvn7+yszM1N+fn4Vc6CAIQ6HQxs3btSQIUPk7u5uuhwAAGBYVZ7rVrtLua6k4ESpfv3LEz3p6ek6fvy4evXq5RJKJMnb21u9evXSwYMHdfToUWd7cnKyc9svDR48WJKUkpJSWYcAAAAAoAjV7lKu4nz33XfatGmTAgMDFR4eLknO9SChoaFFvic0NFTx8fFKT09XUFCQsrOzlZGRoY4dO8rNza3I/j8ftyh5eXnKy8tzvrbb7ZIu/0uzw+G4uoMDqomCzzCfZQAAIFXtOUGNCCYOh0Njx45VXl6e5s+f7wwVmZmZkiR/f/8i31cw3VTQr6z9izJv3jzNnj27UHtCQoIaNmxYmsMBqr3ExETTJQAAgGogJyenyvZV7YNJfn6+oqKitHXrVk2ePFljx441Ws+MGTM0bdo052u73a6goCBFRkayxgQ1nsPhUGJiogYNGsQaEwAA4Lw6qCpU62CSn5+viRMnavXq1brnnnv0xhtvuGwvmPkoboaj4AdZ0K+s/Yvi6enpvCPYz7m7u3Mih1qDzzMAAJBUpecD1Xbxe35+viZMmKDly5fr7rvv1rJly1Svnmu5Ja0J+eUaFG9vbwUGBurQoUO6dOlSif0BAAAAVI1qGUwKQsmKFSt05513auXKlcUuVm/evLl27Nih7Oxsl23Z2dnasWOHWrVqpaCgIGd7RESEc9svFTy/5Je3EQYAAABQuapdMCm4fGvFihX6/e9/r1WrVhUZSiTJZrNp0qRJysrK0pw5c1y2zZkzR1lZWZo8ebJL+5QpUyRJTz/9tC5cuOBs/+ijj5ScnKzIyEgFBwdX8FEBAAAAuJJq94DFWbNmafbs2fLx8dHDDz/sfGbJz40YMUJdunSRdHlmpFevXvriiy8UGRmprl27KjU1VQkJCerRo4dSUlLk5eXl8v7JkycrLi5OHTp00NChQ5WRkaE1a9bIx8dHO3fuVFhYWKnr5QGLqE14wCIAAPi5qjzXrXaL3w8fPixJysrK0nPPPVdkn5CQEGcw8fb2VkpKimbNmqW1a9cqKSlJgYGBeuyxxxQTE1MolEhSbGyswsPD9eabb+qvf/2rfHx8NHLkSD333HNq06ZNZR0aAAAAgGJUuxmTmoYZE9QmzJgAAICfq8pz3Wq3xgQAAABA3UMwAQAAAGAcwQQAAACAcQQTAAAAAMYRTAAAAAAYRzABAAAAYBzBBAAAAIBxBBMAAAAAxhFMAAAAABhHMAEAAABgHMEEAAAAgHEEEwAAAADGEUwAAAAAGEcwAQAAAGAcwQQAAACAcQQTAAAAAMYRTAAAAAAYRzABAAAAYBzBBAAAAIBxBBMAAAAAxhFMAAAAABhHMAEAAABgHMEEAAAAgHEEEwAAAADGEUwAAAAAGEcwAQAAAGAcwQQAAACAcQQTAAAAAMYRTAAAAAAYRzABAAAAYBzBBAAAAIBx9Uvb0c/PT127dtUNN9ygbt26qWvXrmrfvr1sNltl1gcAAACgDih1MPHy8tLWrVu1detWZxjx8vJS586d1a1bN2dY6dChg+rVYyIGAAAAQOmVOpicPHlS33//vXbv3q3U1FTt3r1bu3fv1s6dO7Vz505nWGnQoIHCw8OdQeXee++ttOIBAAAA1A42y7Ks8gxw4sQJZ0gpCC3ff/+9bDabbDabLl68WFG1Vkt2u13+/v7KzMyUn5+f6XKAcnE4HNq4caOGDBkid3d30+UAAADDqvJct9QzJsVp1qyZhg4dqqFDhyovL08bN27Ua6+9ps2bN1dEfQAAAADqgHIHk5ycHP3rX//Se++9p48++kjZ2dmyLEvdunXTqFGjKqJGAAAAALXcVQWT8+fP65///Kfee+89xcfH66effpIk3XTTTRo1apR+97vfKTg4uEILBQAAAFB7lTqYnDt3Ths2bNB7772nTZs2KS8vT25uboqIiNCoUaM0cuRINWvWrDJrBQAAAFBLlTqYXHvttbp48aI8PDw0YMAA/e53v9Pw4cN1zTXXVGZ9AFCpMjIylJGRUez2wMBABQYGVmFFAADUTaUOJg6HQzabTX379lWfPn3UsmVLnlcCoMaLjY3V7Nmzi90eExOjWbNmVV1BAADUUaUOJjfccIO+/PJLxcfHKyEhwdnesmVLde3a1fnckq5du+pXv/pVpRQLABUtOjpaw4YNU25urnr37i1J2r59u7y8vCSJ2RIAAKpImZ5j4nA4tHfvXucDFlNTU/Wf//xHeXl5lwf770MWmzdv7gwrzzzzTOVUXk3wHBPUJnX5OSbZ2dny8fGRJGVlZcnb29twRQAAmFeV57rlfsDipUuX9OWXX7o8YPGLL75Qbm6ubDabLl26VFG1VksEE9QmBBOCCQAAP1ejHrDo5uamTp06qVOnTpowYYIkKT8/X1999ZV2795d7gIBAAAA1H7lDiZFqVevnjp06KAOHTpUxvAAAAAAahluqwUAAADAOIIJAAAAAOMIJgAAAACMI5gAAAAAMI5gAgAAAMA4ggkAAAAA4wgmAAAAAIwjmAAAAAAwjmACAJK+27xZcyWtlvRp//46nJhouiQAAOoUggmAOm/bxIkKGz5cT0gaLSnis88UFBmp7ZMmmS4NAIA6g2ACoE47nJio3y5dKjdJ9SXnn/Uk3bR4sY5s3my0PgAA6opqGUxWrVql6Ohode/eXZ6enrLZbFq2bFmRfWfNmiWbzVbs1+HDh4t8X3x8vCIiIuTr6ys/Pz/169dPmzkBAeqcw888I6uIdpskS9KhmTOruCIAAOqm+qYLKMrMmTN15MgRBQQEKDAwUEeOHCnxPePHj1dISEih9kaNGhVqW7VqlcaOHaumTZsqKipKkrRmzRoNGjRI7777ru64445yHgGAmsLj+HHZitlm++92AABQ+aplMImLi1NoaKiCg4P1/PPPa8aMGSW+JyoqSn379i2x39mzZ/XHP/5RAQEBSk1NVYsWLSRJ06dP1w033KD7779fgwcPlq+vb3kPA0ANcKF5c1nffVfkNuu/2wEAQOWrlpdyDRw4UMHBwZUy9j/+8Q+dO3dOf/zjH52hRJJatGihqVOn6vTp03r//fcrZd8Aqp+QZ591Xrb1c5Yuz5i0+vOfq74oAADqoGoZTK7G1q1bNX/+fL344otav369srKyiuyXnJwsSYqMjCy0bfDgwZKklJSUSqsTQPUSMmiQdt57r/IlOSRd/O9XvqSd996r4AEDjNYHAEBdUS0v5boaMTExLq8bNWqkv/71rxo3bpxLe3p6uiQpNDS00BgFbQV9ipKXl6e8vDzna7vdLklyOBxyOBxXVzxQTRR8huvaZ/nG11/X17ffrg9GjFCIpKbduytkzhzdOGBAnftZAADwc1X5/8EaH0w6d+6sJUuWqG/fvgoMDNSJEyf04Ycf6plnnlFUVJQaNWqkYcOGOftnZmZKkvz9/QuN5efn59KnKPPmzdPs2bMLtSckJKhhw4blPRygWkisgw8X/Omnn/TUf79/5/HHlZ2Xpy83bjRaEwAApuXk5FTZvmp8MBk5cqTL65CQEE2dOlXt27fXoEGDNHPmTJdgUl4zZszQtGnTnK/tdruCgoIUGRnpDDZATeVwOJSYmKhBgwbJ3d3ddDlVKjs72/n94MGD5e3tbbAaAACqh4Krg6pCjQ8mxRkwYIDatGmjvXv3ym63O0NDwUxJZmamrrnmGpf3FPzgi5pNKeDp6SlPT89C7e7u7nXuRA61V138PP/8eOvi8QMAUJSq/P9hrVn8XpSAgABJrlNQV1pHcqX1JwAAAAAqT60NJtnZ2fryyy/l7e3tDCiSFBERIenympBfio+Pd+kDAAAAoGrU6GBy/vx5ffPNN4Xac3NzNXnyZJ0/f16jR49W/fr/u2Jt9OjR8vf318svv6xjx445248dO6ZXXnlFAQEBhdatAAAAAKhc1XKNSVxcnLZv3y5J2rt3r7Ot4BkkvXv31qRJk3TmzBm1a9dOPXr0UPv27dWsWTOdPHlSmzZt0rFjxxQeHq4XX3zRZezGjRvrlVde0dixY9W1a1fdeeedkqQ1a9bozJkzWrNmDU99BwAAAKpYtQwm27dv1/Lly13aduzYoR07djhfT5o0SU2aNNEDDzygzz77TBs3btTZs2fl5eWl9u3b66GHHtLUqVPl5eVVaPx77rlHAQEBmjt3rpYuXSqbzaZu3bpp5syZGjhwYKUfHwAAAABXNsuyLNNF1GR2u13+/v7KzMzkdsGo8RwOhzZu3KghQ4bUubtSZWdny8fHR5KUlZXF7YIBAFDVnuvW6DUmAAAAAGoHggkAAAAA4wgmAAAAAIyrlovfAaCqZGRkKCMjQ7m5uc62tLQ0540zAgMDFRgYaKo8AADqDGZMANRpsbGx6tatm3r37u1s6927t7p166Zu3bopNjbWYHUAANQdzJgAqNOio6M1bNiwYrczWwIAQNUgmACo07hUCwCA6oFLuQAAAAAYRzABAAAAYBzBBAAAAIBxBBMAAAAAxhFMAAAAABhHMAEAAABgHMEEAAAAgHEEEwAAAADGEUwAAAAAGEcwAQAAAGAcwQQAAACAcQQTAAAAAMYRTAAAAAAYRzABAAAAYBzBBAAAAIBxBBMAAAAAxhFMAAAAABhHMAEAAABgHMEEAAAAgHEEEwAAAADGEUwAAAAAGEcwAQAAAGAcwQQAAACAcQQTAAAAAMYRTAAAAAAYRzABAAAAYBzBBAAAAIBxBBMAAAAAxhFMAAAAABhHMAEAAABgHMEEAAAAgHEEEwAAAADGEUwAAAAAGEcwAQAAAGAcwQQAAACAcQQTAAAAAMYRTAAAAAAYRzABAAAAYBzBBAAAAIBxBBMAAAAAxhFMAAAAABhHMAEAAABgHMEEAAAAgHEEEwAAAADGEUwAAAAAGEcwAQAAAGAcwQQAAACAcQQTAAAAAMYRTAAAAAAYRzABAAAAYFy1DCarVq1SdHS0unfvLk9PT9lsNi1btqzY/na7XdOmTVNwcLA8PT0VEhKiJ554QllZWUX2z8/P18svv6zw8HB5eXmpadOmuvvuu3Xw4MFKOiIAAAAAV1Itg8nMmTP15ptv6siRIwoMDLxi3+zsbEVERGjRokVq166dHn30UbVt21YLFixQ//799dNPPxV6T3R0tB566CFZlqWHHnpIt9xyi9atW6cePXooPT29sg4LAAAAQDGqZTCJi4vT4cOHderUKd13331X7PvCCy8oLS1N06dPV3x8vJ5//nnFx8dr+vTp2rVrlxYtWuTSPykpSXFxcerTp49SU1M1f/58rVy5UuvXr9ePP/6oqVOnVuahAQAAAChCtQwmAwcOVHBwcIn9LMtSXFycfHx89PTTT7tse/rpp+Xj46O4uDiX9rfeekuSNGfOHHl4eDjbb731VvXt21cJCQn67rvvKuAoAAAAAJRWtQwmpZWenq7jx4+rV69e8vb2dtnm7e2tXr166eDBgzp69KizPTk52bntlwYPHixJSklJqdzCAQAAALiob7qA8ihYDxIaGlrk9tDQUMXHxys9PV1BQUHKzs5WRkaGOnbsKDc3tyL7/3zcouTl5SkvL8/52m63S5IcDoccDsdVHwtQHRR8hvksAwAAqWrPCWp0MMnMzJQk+fv7F7ndz8/PpV9Z+xdl3rx5mj17dqH2hIQENWzYsJSVA9VbYmKi6RIAAEA1kJOTU2X7qtHBxIQZM2Zo2rRpztd2u11BQUGKjIx0BhugpnI4HEpMTNSgQYPk7u5uuhwAAGBYwdVBVaFGB5OCmY/iZjgKfpAF/cravyienp7y9PQs1O7u7s6JHGoNPs8AAEBSlZ4P1OjF7yWtCfnlGhRvb28FBgbq0KFDunTpUon9AQAAAFSNGh9Mmjdvrh07dig7O9tlW3Z2tnbs2KFWrVopKCjI2R4REeHc9kvx8fGSpD59+lRu4QAAAABc1OhgYrPZNGnSJGVlZWnOnDku2+bMmaOsrCxNnjzZpX3KlCmSLj/n5MKFC872jz76SMnJyYqMjCzVM1QAAAAAVBybZVmW6SJ+KS4uTtu3b5ck7d27V6mpqerVq5euu+46SVLv3r01adIkSZdnRnr16qUvvvhCkZGR6tq1q1JTU5WQkKAePXooJSVFXl5eLuNPnjxZcXFx6tChg4YOHaqMjAytWbNGPj4+2rlzp8LCwkpdq91ul7+/vzIzM1n8jhrP4XBo48aNGjJkCGtMAABAlZ7rVsvF79u3b9fy5ctd2nbs2OFy+VVBMPH29lZKSopmzZqltWvXKikpSYGBgXrssccUExNTKJRIUmxsrMLDw/Xmm2/qr3/9q3x8fDRy5Eg999xzatOmTeUeHAAAAIBCquWMSU3CjAlqE2ZMAADAz1XluW6NXmMCAAAAoHYgmAAAAAAwjmACAAAAwDiCCQAAAADjCCYAAAAAjCOYAAAAADCOYAIAAADAOIIJAAAAAOMIJgAAAACMI5gAAAAAMI5gAgAAAMA4ggkAAAAA4wgmAAAAAIwjmAAAAAAwjmACAAAAwDiCCQAAAADjCCYAAAAAjCOYAAAAADCOYAIAAADAOIIJAAAAAOMIJgAAAACMI5gAAAAAMI5gAgAAAMA4ggkAAAAA4wgmAAAAAIwjmAAAAAAwjmACAAAAwDiCCQAAAADjCCYAAAAAjCOYAAAAADCOYAIAAADAOIIJAAAAAOMIJgAAAACMI5gAAAAAMI5gAgAAAMA4ggkAAAAA4wgmAAAAAIwjmAAAAAAwjmACAAAAwDiCCQAAAADjCCYAAAAAjCOYAAAAADCOYAIAAADAOIIJAAAAAOMIJgAAAACMI5gAAAAAMI5gAgAAAMA4ggkAAAAA4wgmAAAAAIwjmAAAAAAwjmACAAAAwDiCCQAAAADjCCYAAAAAjCOYAAAAADCOYAIAAADAOIIJAAAAAOMIJgAAAACMqxXBJCQkRDabrcivvn37Fuqfl5enZ599VqGhoWrQoIGaN2+uKVOm6Icffqj64gEAAACovukCKoq/v78eeeSRQu0hISEur/Pz8zV8+HDFx8frN7/5jUaNGqX09HTFxcVp8+bN+vTTT9W0adOqKRoAAACApFoUTBo1aqRZs2aV2G/58uWKj4/X3Xffrbfffls2m02S9MYbb+j+++/XzJkzFRsbW8nVAgAAAPi5WnEpV1m89dZbkqR58+Y5Q4kkRUdHq3Xr1nr77beVm5trqjwAAACgTqo1wSQvL0/Lli3T3Llz9corr+jf//53oT4//fST/v3vf6tt27YKDg522Waz2TRo0CBlZ2fr888/r6qyAQAAAKgWXcp14sQJTZgwwaWtR48e+vvf/642bdpIkr799lvl5+crNDS0yDEK2tPT03XzzTcX2ScvL095eXnO13a7XZLkcDjkcDjKfRyASQWfYT7LAABAqtpzgloRTCZMmKCbb75ZHTt2lI+Pj7755hstXLhQK1eu1IABA7R37175+voqMzNT0uWF8kXx8/OTJGe/osybN0+zZ88u1J6QkKCGDRtWwNEA5iUmJpouAQAAVAM5OTlVtq9aEUxiYmJcXnfp0kUrVqyQJK1cuVJvvfWWpk2bViH7mjFjhstYdrtdQUFBioyMdAYboKZyOBxKTEzUoEGD5O7ubrocAABgWMHVQVWhVgST4kRHR2vlypXasWOHpk2b5pwpKW5GpOAHX9yMiiR5enrK09OzULu7uzsncqg1+DwDAABJVXo+UGsWvxclICBAkpSdnS1Jat26terVq6f09PQi+xe0F7cGBQAAAEDlqNXBpODOXAUPWfTy8lLPnj21f/9+HTlyxKWvZVlKTEyUt7e3unfvXtWlAgAAAHVajQ8mX3/9dZGLcr7++mtNnz5dkjRmzBhn+5QpUyRdXitiWZazPTY2VgcPHtQf/vAHeXl5VXLVAAAAAH6uxq8xeeedd7Rw4UL16dNHwcHB8vb21jfffKONGzfK4XBoxowZ6tOnj7P/+PHjtWbNGv3973/XoUOHFBERoQMHDmjdunVq1aqV/vznPxs8GgAAAKBuqvHBpF+/fvrqq6+0Z88ebdu2TTk5OQoICNCQIUP0wAMPKDIy0qV/vXr1tGHDBj3//PNauXKlFi1apCZNmujee+/Vn//8ZzVt2tTQkQAAAAB1l836+fVMKDO73S5/f39lZmZyu2DUeA6HQxs3btSQIUO4KxeAMsvIyFBGRkax2wMDAxUYGFiFFQEor6o8163xMyYAAKB6iI2NLfIhxAViYmI0a9asqisIQI1CMAEAABUiOjpaw4YNU25urnr37i1J2r59u/OmMsyWALgSggkAAKgQBZdqFTw/TJK6dOkib29vg1UBqClq/O2CAQAAANR8BBMAAAAAxhFMAABAhfpu82bNlbRa0qf9++twYqLpkgDUAAQTAABQYbZNnKiw4cP1hKTRkiI++0xBkZHaPmmS6dIAVHMEEwAAUCEOJybqt0uXyk2X765T8Gc9STctXqwjmzcbrQ9A9UYwAQAAFeLwM8+oqKc22yRZkg7NnFnFFQGoSQgmAACgQngcPy5bMdts/90OAMUhmAAAgApxoXnzImdMpMszJheaN6/KcgDUMAQTAABQIUKefdZ52dbPWbo8Y9Lqz3+u+qIA1BgEEwAAUCFCBg3SznvvVb4kh6SL//3Kl7Tz3nsVPGCA0foAVG8EEwAAUGF6x8Up/Z//1AJJ/5CU0rOnjm3apN5xcaZLA1DN1TddAAAAqF2C+vXTU//9PmvLFnl7exutB0DNwIwJAAAAAOOYMQEAABUiIyNDGRkZys3NdbalpaXJy8tLkhQYGKjAwEBT5QGo5pgxAQAAFSI2NlbdunVT7969nW29e/dWt27d1K1bN8XGxhqsDkB1x4wJAACoENHR0Ro2bFix25ktAXAlBBMAAFAhuFQLQHlwKRcAAAAA4wgmAAAAAIwjmAAAAAAwjmACAAAAwDiCCQAAAADjCCYAAAAAjCOYAAAAADCOYAIAAADAOIIJAAAAAOMIJgAAAACMI5gAAAAAMI5gAgAAAMA4ggkAAAAA4wgmAAAAAIwjmAAAAAAwjmACAAAAwDiCCQAAAADjCCYAAAAAjCOYAAAAADCOYAIAAADAOIIJAAAAAOMIJgAAAACMI5gAAAAAMI5gAgAAAMA4ggkAAAAA4wgmAAAAAIyrb7oAAAAAAP+TkZGhjIyMYrcHBgYqMDCwCiuqGgQTAAAAoBqJjY3V7NmzdZ2kiZJCJB2WtETSAUkxMTGaNWuWsfoqC8EEAAAAqEaio6PVZc8e3f7BB7Ik2SRZkp6U9M/hw3VjdLTZAisJa0wAAACAaiTv//0/3f7BB3LT5VmEgj/rSbp9wwZd2LfPaH2VhWACAAAAVCOHn3lGVhHtBTMnh2bOrOKKqgbBBAAAAKhGPI4fl62Ybbb/bq+NCCYAAABANXKhefMiZ0ykyzMmF5o3r8pyqgzBBAAAAKhGQp591nnZ1s8VLIRv9ec/V31RVYBgAgAAAFQjnh076p/DhytfkkPSxf9+5evyXbk8rr/eaH2VhWACAAAAVCOxsbEauWGD2kpaIOkfkl6U1FbSyA0bFBsba7S+ysJzTAAAAIBqJDo6WsOGDXNpaytp8H+/r41PfZfqeDDZtWuXYmJi9Mknn8jhcCg8PFzTpk3T6NGjTZcGAACAOiowMLDWho8rqbPBJCkpSYMHD1aDBg101113ydfXV2vXrtWdd96po0eP6rHHHjNdIgAAAFBn2CzLKu5uZLXWxYsX1a5dOx07dkyffvqpunTpIknKzMxUz549dfjwYX3zzTcKDg4ucSy73S5/f39lZmbKz8+vkisHKpfD4dDGjRs1ZMgQubu7my4HAAAYVpXnunVy8fuWLVv07bffasyYMc5QIkn+/v566qmndOHCBS1fvtxcgQAAAEAdUyeDSXJysiQpMjKy0LbBgy8vK0pJSanKkgAAAIA6rU6uMUlPT5ckhYaGFtrWrFkz+fj4OPv8Ul5envLy8pyv7Xa7pMuXwDgcjkqoFqg6BZ9hPssAAECq2nOCOhlMMjMzJV2+dKsofn5+zj6/NG/ePM2ePbtQe0JCgho2bFhxRQIGJSYmmi4BAABUAzk5OVW2rzoZTMpjxowZmjZtmvO13W5XUFCQIiMjWfyOGs/hcCgxMVGDBg1i8TsAAHBeHVQV6mQwKZgpKW5WxG63q3HjxkVu8/T0lKenZ6F2d3d3TuRQa/B5BgAAkqr0fKBOLn4vWFtS1DqSEydOKCsrq8j1JwAAAAAqR50MJhEREZIurwv5pfj4eJc+AAAAACpfnQwmAwYMUOvWrbV69WqlpaU52zMzMzV37lx5eHho3Lhx5goEAAAA6pg6ucakfv36iouL0+DBg9WnTx/ddddd8vX11dq1a3XkyBEtWLBAISEhpssEAAAA6ow6GUwkqV+/ftq+fbtiYmK0Zs0aORwOhYeHa/78+brzzjtNlwcAAADUKXU2mEhSz5499dFHH5kuAwAAAKjz6uQaEwAAAADVC8EEAAAAgHEEEwAAAADGEUwAAAAAGEcwAQAAAGAcwQQAAACAcXX6dsEVwbIsSZLdbjdcCVB+DodDOTk5stvtcnd3N10OAAAwrOAct+CctzIRTMrp/PnzkqSgoCDDlQAAAACV4/z58/L396/Ufdisqog/tVh+fr6OHz8uX19f2Ww20+VctR49emjXrl2myyiXmnYM1bFeu92uoKAgHT16VH5+fqbLAVCFquPfSajZ+EzVDCX9nizL0vnz59W8eXPVq1e5q0CYMSmnevXqqUWLFqbLKDc3N7cafyJa046hOtfr5+dXbWsDUDmq899JqJn4TNUMpfk9VfZMSQEWv0OS9OCDD5ouodxq2jHUtHoB1G78nYSKxmeqZqhOvycu5QLgZLfb5e/vr8zMTP6VCwAAVClmTAA4eXp6KiYmRp6enqZLAQAAdQwzJgAAAACMY8YEAAAAgHEEEwAAAADGEUwAAAAAGEcwAQAAFW7kyJFq3Lix7rjjDtOloJbgM1UzlOf3RDABcFXWrVunQYMGqUmTJrLZbDp8+LDpkgBUIw8//LBWrFhhugzUInymaoby/J4IJgCuSnZ2tvr06aNnn33WdCkAqqG+ffvK19fXdBmoRfhM1Qzl+T0RTABclbFjx+rpp59W3759TZcC1Bnz5s1Tjx495Ovrq1/96lcaMWKE9u/fX6H72Lp1q26//XY1b95cNptN69evL7Lfq6++qpCQEDVo0EA33nijPvvsswqtA1Xj9ddfV6dOneTn5yc/Pz/ddNNN+uijjyp0H3ymKtbzzz8vm82mRx55pELHrQ6/J4IJUEutWrVK0dHR6t69uzw9PWWz2bRs2bIrvmfXrl0aMmSIGjVqJG9vb/3mN7/Ru+++WzUFAyhRSkqKHnzwQX366adKTEyUw+FQZGSksrOzi+y/Y8cOORyOQu379u3TyZMni3xPdna2OnfurFdffbXYOtasWaNp06YpJiZGqamp6ty5swYPHqwffvjh6g4MxrRo0ULPP/+8du/erc8//1z9+/fX8OHD9eWXXxbZn8+UWbt27VJsbKw6dep0xX419vdkAaiVgoODLUlWQECA8/ulS5cW23/Lli2Wu7u75evra02ePNmaNm2a830LFiwo9n179+61JFmHDh2q+IMAcEU//PCDJclKSUkptO3SpUtW586drTvuuMO6ePGis/3rr7+2rr32Wmv+/Pklji/Jev/99wu19+zZ03rwwQdd9tW8eXNr3rx5Lv2SkpKsUaNGleGIUB00btzYiouLK9TOZ8qs8+fPW6GhoVZiYqIVERFhPfzww0X2q8m/J2ZMgFoqLi5Ohw8f1qlTp3Tfffddse/Fixc1efJk1atXT1u3btWbb76pv/zlL/riiy8UFhamp556SkeOHKmiygGUVmZmpiSpSZMmhbbVq1dPGzdu1J49ezRu3Djl5+fr22+/Vf/+/TVixAg9+eSTV7XPCxcuaPfu3Ro4cKDLvgYOHKidO3de3YGgWrh06ZLeeecdZWdn66abbiq0nc+UWQ8++KCGDh3q8nMqSk3+PRFMgFpq4MCBCg4OLlXfLVu26Ntvv9WYMWPUpUsXZ7u/v7+eeuopXbhwQcuXL6+kSgFcjfz8fD3yyCPq1auXOnbsWGSf5s2ba8uWLdq+fbvGjBmj/v37a+DAgXr99dever+nT5/WpUuXdO2117q0X3vttTpx4oTz9cCBA/X73/9eGzduVIsWLTjBrMb27t0rHx8feXp66r777tP777+v66+/vsi+fKbMeOedd5Samqp58+aVqn9N/T3Vv+rqANQaycnJkqTIyMhC2wYPHizp8rXtAKqPBx98UP/v//0/bd++/Yr9WrZsqZUrVyoiIkKtW7fW4sWLZbPZKr2+TZs2Vfo+UDHatm2rtLQ0ZWZm6r333tP48eOVkpJSbDjhM1W1jh49qocffliJiYlq0KBBqd9XE39PzJgAUHp6uiQpNDS00LZmzZrJx8fH2afAjz/+qLS0NOcdgfbt26e0tDT9+OOPlV8wUMdNnTpVH374oZKSktSiRYsr9j158qSmTJmi22+/XTk5OXr00UfLte+AgAC5ubkVWkB78uRJNWvWrFxjwwwPDw9dd9116tatm+bNm6fOnTvrr3/9a7H9+UxVrd27d+uHH35Q165dVb9+fdWvX18pKSn629/+pvr16+vSpUtFvq8m/p4IJgCc16n7+/sXud3Pz8/Zp8AHH3ygG264wflk16FDh+qGG27QBx98ULnFAnWYZVmaOnWq3n//fW3ZskWtWrW6Yv/Tp09rwIABat++vdatW6fNmzdrzZo1evzxx6+6Bg8PD3Xr1k2bN292tuXn52vz5s1FrktAzZOfn6+8vLwit/GZqnoDBgzQ3r17lZaW5vzq3r27/vCHPygtLU1ubm6F3lNTf09cygXgqkRFRSkqKsp0GUCd8uCDD2r16tXasGGDfH19ndd1+/v7y8vLy6Vvfn6+br31VgUHB2vNmjWqX7++rr/+eiUmJqp///769a9/XeS/oGZlZenAgQPO14cOHVJaWpqaNGmili1bSpKmTZum8ePHq3v37urZs6deeuklZWdna8KECZV49KgMM2bM0K233qqWLVvq/PnzWr16tZKTkxUfH1+oL58pM3x9fQutI/P29tY111xT5PqyGv17KvN9vADUOPPmzbvi7YLvuOMOS5L1+eefF7ndx8fHCgoKqsQKAZSGpCK/ivtvOyEhwcrNzS3Unpqaah09erTI9yQlJRW5j/Hjx7v0e/nll62WLVtaHh4eVs+ePa1PP/20vIcHAyZOnGgFBwdbHh4eVtOmTa0BAwZYCQkJxfbnM1U9XOl2wZZVc39PNsuyrIqJOACqq+eff14zZszQ0qVLi5zleOqppzRv3jz9/e9/11133eWy7cSJEwoMDFT//v1dpm8BAAAqEmtMACgiIkKSlJCQUGhbwXR+QR8AAIDKQDABoAEDBqh169ZavXq10tLSnO2ZmZmaO3euPDw8NG7cOHMFAgCAWo9LuYBaKi4uzvl8g7179yo1NVW9evXSddddJ0nq3bu3Jk2a5OyflJSkwYMHq0GDBrrrrrvk6+urtWvX6siRI1qwYIEee+wxI8cBAADqBoIJUEtFRUVd8Wnt48eP17Jly1zaPvvsM8XExOiTTz6Rw+FQeHi4pk2bpjvvvLOSqwUAAHUdwQQAAACAcawxAQAAAGAcwQQAAACAcQQTAAAAAMYRTAAAAAAYRzABAAAAYBzBBAAAAIBxBBMAAAAAxhFMAAAAABhHMAEAAABgHMEEAAAAgHEEEwAAAADGEUwAADVCTk6O5syZo7CwMDVo0ECtW7fWX/7yF+3Zs0c2m01PPPGE6RIBAOVQ33QBAACU5Pz58+rfv78+//xz9e/fXyNHjtSBAwf0xBNPKDIyUpJ0ww03GK4SAFAeBBMAQLU3ceJEpaamauXKlbrnnnuc7QsWLHDOlHTp0sVQdQCAimCzLMsyXQQAAMXZvHmzBg4cqPvuu0+vv/66y7ZTp07pV7/6lby8vHT+/Hm5ubkZqhIAUF6sMQEAVGuvvvqqJOnJJ58stK1JkyaSpPDwcEIJANRwBBMAQLW2adMmXXfddWrVqlWhbRkZGZJYXwIAtQHBBABQbZ07d07nz59XUFBQkds3bdokifUlAFAbEEwAANWWu7u7JOnHH38stO3ChQt64YUXJDFjAgC1AcEEAFBteXt7Kzg4WHv37tXevXud7Xl5eYqKitJXX32levXqKTw83GCVAICKQDABAFRrjz/+uPLz89WnTx/df//9evTRR3X99dfr1KlTatCggdq2bauGDRuaLhMAUE48xwQAUK09+OCDOnfunN544w0tWbJErVu3VnR0tIYPH6527dqxvgQAagmCCQCgWrPZbJo5c6Zmzpzp0r5u3TpJrC8BgNqCS7kAADXSF198IYk7cgFAbUEwAQDUSGlpaZKYMQGA2sJmWZZluggAAMoqJCREFy9e1LFjx0yXAgCoAAQTAAAAAMZxKRcAAAAA4wgmAAAAAIwjmAAAAAAwjmACAAAAwDiCCQAAAADjCCYAAAAAjCOYAAAAADCOYAIAAADAOIIJAAAAAOMIJgAAAACM+/86QJCJOLkRHgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,6))\n",
    "\n",
    "plt.errorbar(10**np.asarray(log10y_center),N_clusters_y_theory,yerr=np.sqrt(N_clusters_y_theory),\n",
    "            marker='o',\n",
    "            markersize=5.,\n",
    "            fillstyle='none',\n",
    "            ls='none',\n",
    "            label=r'$\\mathrm{CLASS\\_SZ\\,\\,binned}$',\n",
    "             capsize=3,c='k')\n",
    "plt.errorbar(10**(np.asarray(log10y_center)),Nq_binned,#yerr=np.sqrt(Nq_binned),\n",
    "             marker='o',markersize=5.,ls='none',\n",
    "            label=r'$\\mathrm{CLASS\\_SZ\\,\\,from\\,\\,unbinned}$',\n",
    "             capsize=5,c='red')\n",
    "\n",
    "\n",
    "\n",
    "plt.xlabel('$q$', fontsize=14)\n",
    "plt.ylabel('$N$', fontsize=14)\n",
    "\n",
    "plt.xticks(fontsize=14)\n",
    "plt.yticks(fontsize=14)\n",
    "plt.legend(fontsize=13,frameon=False)\n",
    "plt.grid()\n",
    "plt.tight_layout()\n",
    "plt.xscale('log')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}