{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Process PV catalogues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from os.path import join\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from h5py import File\n",
    "from astropy.coordinates import match_coordinates_sky, SkyCoord\n",
    "from astropy import units as u\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "SPEED_OF_LIGHT = 299_792.458"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Supernovae data\n",
    "\n",
    "### LOSS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "a2dir = \"/mnt/extraspace/rstiskalek/catalogs/PV/PV_Supranta/A2\"\n",
    "\n",
    "names = [\"z_CMB\", \"mB\", \"x1\", \"c\", \"e_mB\", \"e_x1\", \"e_c\", \"RA\", \"DEC\"]\n",
    "dtype = [(n, np.float32) for n in names]\n",
    "data = np.genfromtxt(join(a2dir, \"loss.csv\"), delimiter=\",\", skip_header=1,\n",
    "                     usecols=[5 + n for n in range(len(names))])\n",
    "\n",
    "loss_data = np.empty(len(data), dtype=dtype)\n",
    "for i, n in enumerate(names):\n",
    "    loss_data[n] = data[:, i]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Foundation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "names = [\"z_CMB\", \"RA\", \"DEC\", \"x1\", \"mB\", \"c\", \"peak\", \"e_peak\", \"e_x1\", \"e_mB\", \"e_c\"]\n",
    "dtype = [(n, np.float32) for n in names]\n",
    "data = np.genfromtxt(join(a2dir, \"foundation.csv\"), delimiter=\",\", skip_header=1,\n",
    "                     usecols=[3 + n for n in range(len(names))])\n",
    "\n",
    "foundation_data = np.empty(len(data), dtype=dtype)\n",
    "for i, n in enumerate(names):\n",
    "    foundation_data[n] = data[:, i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pantheon+, all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {},
   "outputs": [],
   "source": [
    "fpath_full = \"/mnt/extraspace/rstiskalek/catalogs/PV/Pantheon+SH0ES.dat\"\n",
    "fpath_group = \"/mnt/extraspace/rstiskalek/catalogs/PV/pantheon+_groups.hdf5\"\n",
    "\n",
    "data_full = np.genfromtxt(fpath_full, names=True, dtype=None, encoding=None)\n",
    "data_group = {}\n",
    "\n",
    "with File(fpath_group, \"r\") as f:\n",
    "    for key in f.keys():\n",
    "        data_group[key] = f[key][:]\n",
    "\n",
    "# data = data[data[\"zCMB\"] < 0.1]\n",
    "\n",
    "keys = [\"zCMB\", \"mB\", \"mBERR\", \"x1\", \"x1ERR\", \"c\", \"cERR\", \"RA\", \"DEC\",\n",
    "        \"VPEC\", \"VPECERR\", \"biasCor_m_b\", \"biasCorErr_m_b\"]\n",
    "pantheon_data = {}\n",
    "for key in keys:\n",
    "    pantheon_data[key] = data[key]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum 2D separation is: 0.0005450609539575544 deg\n"
     ]
    }
   ],
   "source": [
    "match_coord = SkyCoord(data_full[\"RA\"], data_full[\"DEC\"], unit=\"deg\")\n",
    "catalogue_coord = SkyCoord(data_group[\"RAJ2000\"], data_group[\"DEJ2000\"],\n",
    "                           unit=\"deg\")\n",
    "\n",
    "idxs, sep2d, sep3d = match_coordinates_sky(match_coord, catalogue_coord, nthneighbor=1)\n",
    "print(f\"Maximum 2D separation is: {sep2d.max()}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {},
   "outputs": [],
   "source": [
    "zcmb_supernovae = data_group[\"zcmb\"][idxs]\n",
    "zcmb_group = data_group[\"zGcmb\"][idxs]\n",
    "\n",
    "pantheon_data[\"zCMB_SN\"] = zcmb_supernovae\n",
    "pantheon_data[\"zCMB_Group\"] = zcmb_group\n",
    "\n",
    "mask = pantheon_data[\"zCMB\"] < 0.1\n",
    "for key in pantheon_data.keys():\n",
    "    pantheon_data[key] = pantheon_data[key][mask]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAoi0lEQVR4nO3dfUxVd4L/8c8FBFynQCsjFEqhpW7drTw4Vh60WZkpKU7YUDbZlZL4MMZ0d8iumjA1lcbVTua3QTd1a60GY9PZacP4MG66rivGBnGYtAq1Ag5rZ2qmZIoWChZnAYdBsHD2D3/cchWQi/fe8733vl/JyamH77l8z4nt/fT76LAsyxIAAIDBQuyuAAAAwL0QWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxguzuwKeMDo6qs7OTj3wwANyOBx2VwcAAEyDZVm6ceOGEhISFBIydRtKQASWzs5OJSUl2V0NAAAwA1evXtUjjzwyZZmACCwPPPCApNsPHBUVZXNtAADAdPT39yspKcn5PT6VgAgsY91AUVFRBBYAAPzMdIZzMOgWAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAwJSqG9u1bMcZVTe221YHAgsAAJhSVX2bOnoHVVXfZlsdCCwAAGBKZXmpSoyZrbK8VNvq4LAsy7Ltt3tIf3+/oqOj1dfXp6ioKLurAwAApsGd729aWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgCA4UxYuM1uBBYAAAxnwsJtdptRYNm3b59SUlIUGRmp7OxsnT9/fsryR48e1YIFCxQZGam0tDSdPHnS5ec/+MEP5HA4XI4VK1bMpGoAAAQcExZus5vbgeXIkSMqLy/X9u3b1dzcrIyMDBUUFOjatWsTlj937pxKS0u1fv16tbS0qLi4WMXFxbp06ZJLuRUrVujLL790HocOHZrZEwEAEGBW5STr7JbvaVVOst1VsY3bK91mZ2dryZIl2rt3ryRpdHRUSUlJ2rBhg7Zs2XJX+ZKSEg0MDOjEiRPOazk5OcrMzNT+/fsl3W5h6e3t1bFjx2b0EKx0CwCA//HaSrfDw8NqampSfn7+Nx8QEqL8/Hw1NDRMeE9DQ4NLeUkqKCi4q3x9fb3mzZunJ598UmVlZbp+/fqk9RgaGlJ/f7/LAQAAApdbgaWnp0cjIyOKi4tzuR4XF6eurq4J7+nq6rpn+RUrVujdd99VXV2ddu7cqV/96lf6/ve/r5GRkQk/s7KyUtHR0c4jKSnJnccAAAD3YNrMJCNmCb3wwgsqKipSWlqaiouLdeLECX388ceqr6+fsHxFRYX6+vqcx9WrV31bYQAAApxpM5PcCiyxsbEKDQ1Vd3e3y/Xu7m7Fx8dPeE98fLxb5SXp8ccfV2xsrD777LMJfx4REaGoqCiXAwAAeI5pM5PcCizh4eFavHix6urqnNdGR0dVV1en3NzcCe/Jzc11KS9JtbW1k5aXpC+++ELXr1/Xww8/7E71AACAh5g2M8ntLqHy8nK99dZbeuedd/Tb3/5WZWVlGhgY0Lp16yRJa9asUUVFhbP8pk2bdOrUKe3atUuffvqpXn31VV24cEH/9E//JEn64x//qM2bN6uxsVGff/656urq9Pzzz+uJJ55QQUGBhx4TAAD4szB3bygpKdFXX32lbdu2qaurS5mZmTp16pRzYO2VK1cUEvJNDlq6dKkOHjyorVu36pVXXtH8+fN17NgxLVy4UJIUGhqq1tZWvfPOO+rt7VVCQoKee+45/eQnP1FERISHHhMAAPgzt9dhMRHrsAAA4H+8tg4LAACAHQgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAAB+obqxXct2nFF1Y7vdVYENCCwAAL9QVd+mjt5BVdW32V0V2IDAAgDwC2V5qUqMma2yvFS7qwIbOCzLsuyuxP3q7+9XdHS0+vr6FBUVZXd1AADANLjz/U0LCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAANtVN7Zr2Y4zqm5sd+u+jYdalFpRo42HWrxUM5iCwAIAsF1VfZs6egdVVd/m1n01rZ0asW6fEdgILAAA2y1OflChjttndxSmJyjUcfuMwBZmdwUAAGhq/1+NWLfP7thTukh7Shd5qVYwCS0sAADbleWlKjFmtsryUu2uCgzlsCzLsrsS96u/v1/R0dHq6+tTVFSU3dUBAADT4M73Ny0sAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABANwXlseHLxBYAAD3heXx4QsEFgDAfWF5fPgCC8cBAABbsHAcAAAIKDMKLPv27VNKSooiIyOVnZ2t8+fPT1n+6NGjWrBggSIjI5WWlqaTJ09OWvaHP/yhHA6Hdu/ePZOqAQCAAOR2YDly5IjKy8u1fft2NTc3KyMjQwUFBbp27dqE5c+dO6fS0lKtX79eLS0tKi4uVnFxsS5dunRX2f/8z/9UY2OjEhLoBwUAAN9wO7D827/9m1588UWtW7dOf/mXf6n9+/frz/7sz/TTn/50wvJvvPGGVqxYoc2bN+sv/uIv9JOf/ETf+c53tHfvXpdyHR0d2rBhg37+859r1qxZM3saAAAQkNwKLMPDw2pqalJ+fv43HxASovz8fDU0NEx4T0NDg0t5SSooKHApPzo6qtWrV2vz5s166qmn3KkSAAAIAmHuFO7p6dHIyIji4uJcrsfFxenTTz+d8J6urq4Jy3d1dTn/vHPnToWFhWnjxo3TqsfQ0JCGhoacf+7v75/uIwAAAD9k+yyhpqYmvfHGG/rZz34mh8MxrXsqKysVHR3tPJKSkrxcSwAAYCe3AktsbKxCQ0PV3d3tcr27u1vx8fET3hMfHz9l+Q8++EDXrl3To48+qrCwMIWFham9vV0/+tGPlJKSMuFnVlRUqK+vz3lcvXrVnccAAAB+xq3AEh4ersWLF6uurs55bXR0VHV1dcrNzZ3wntzcXJfyklRbW+ssv3r1arW2turixYvOIyEhQZs3b9b7778/4WdGREQoKirK5QAAAIHLrTEsklReXq61a9fq6aefVlZWlnbv3q2BgQGtW7dOkrRmzRolJiaqsrJSkrRp0yYtX75cu3btUmFhoQ4fPqwLFy7owIEDkqS5c+dq7ty5Lr9j1qxZio+P15NPPnm/zwcAmKHqxnZV1bepLC9Vq3KS7a4OgpzbgaWkpERfffWVtm3bpq6uLmVmZurUqVPOgbVXrlxRSMg3DTdLly7VwYMHtXXrVr3yyiuaP3++jh07poULF3ruKQAAHldV36aO3kFV1bcRWGA79hICAEyIFhZ4mzvf3wQWAABgCzY/BAAAAYXAAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAD+z8VCLUitqtPFQi91VAXyGwAIAfqamtVMj1u0zECwILADgZwrTExTquH0GggW7NQMAAFuwWzMAAAgoBBYAAGA8AgsAADAegQUAfITpyMDMEVgAwEeYjgzMHIEFAHyE6cjAzDGtGQAA2IJpzQAAIKAQWAAAgPEILAAAwHgEFgAIYtWN7Vq244yqG9vtrgowJQILAASxqvo2dfQOqqq+ze6qAFMisACAAexq6SjLS1VizGyV5aX69PcC7mJaMwAYYNmOM+roHVRizGyd3fI9u6sD+ATTmgHAz9DSAUyNFhYAAGALWlgAAEBAIbAAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAD2L3Y8A7CCwA4EGvvX9ZHb2Deu39y3ZXBQgoBBYAAGA8AgsAeNBLBU8qMWa2Xip40u6qAAGFzQ8BwEeqG9tVVd+msrxUrcpJtrs6gO3Y/BAADFRV36aO3kFV1bfZXRXA74TZXQEA8DcbD7WoprVTD0SGqXfwa6UnRuv4hmfueV9ZXqqzhQWAewgsAOCmmtZOjVhS7+DXkqTWjr5p3bcqJ5muIGCG6BICADcVpico1CHFzL79/3zpidE21wgIfAy6BQAAtmDQLQAACCgEFgAAYDwCCwCIPYAA0xFYAAS86YQR1kgBzEZgARDwpgojY2FmcfKDSoyZzRopgKFYhwVAwJtqwbaxMCNJZ7d8z9dVAzBNtLAA8FvTHXeyKidZZ7d8b8JF28ryUmlZAfwAgQWA35rOuJN7hZqpwgwAcxBYAPiljYda1Nk7qNmzQlSWlzppMGEwLRAYCCwA/FJNa6csScNfj2pVTvKkwYQuHyAwzCiw7Nu3TykpKYqMjFR2drbOnz8/ZfmjR49qwYIFioyMVFpamk6ePOny81dffVULFizQnDlz9OCDDyo/P18fffTRTKoGwE+5uw7K2H4+hekJkiYPJnT5AIHB7b2Ejhw5ojVr1mj//v3Kzs7W7t27dfToUV2+fFnz5s27q/y5c+f0V3/1V6qsrNRf//Vf6+DBg9q5c6eam5u1cOFCSdLBgwc1b948Pf744xocHNTrr7+uo0eP6rPPPtO3v/3te9aJvYQA/7dsxxl19A4qMWY2s3WAIOHO97fbgSU7O1tLlizR3r17JUmjo6NKSkrShg0btGXLlrvKl5SUaGBgQCdOnHBey8nJUWZmpvbv3z/lA5w+fVrPPvvsPetEYAH8X3Vju3PqMa0hQHDw2uaHw8PDampqUn5+/jcfEBKi/Px8NTQ0THhPQ0ODS3lJKigomLT88PCwDhw4oOjoaGVkZExYZmhoSP39/S4HAP9G1w2AqbgVWHp6ejQyMqK4uDiX63Fxcerq6prwnq6urmmVP3HihL71rW8pMjJSr7/+umpraxUbGzvhZ1ZWVio6Otp5JCUlufMYAADAzxgzS+i73/2uLl68qHPnzmnFihVauXKlrl27NmHZiooK9fX1OY+rV6/6uLYAAMCX3AossbGxCg0NVXd3t8v17u5uxcfHT3hPfHz8tMrPmTNHTzzxhHJycvT2228rLCxMb7/99oSfGRERoaioKJcDAAAELrcCS3h4uBYvXqy6ujrntdHRUdXV1Sk3N3fCe3Jzc13KS1Jtbe2k5cd/7tDQkDvVAwAAAcrtzQ/Ly8u1du1aPf3008rKytLu3bs1MDCgdevWSZLWrFmjxMREVVZWSpI2bdqk5cuXa9euXSosLNThw4d14cIFHThwQJI0MDCgf/mXf1FRUZEefvhh9fT0aN++fero6NDf/d3fefBRAQCAv3I7sJSUlOirr77Stm3b1NXVpczMTJ06dco5sPbKlSsKCfmm4Wbp0qU6ePCgtm7dqldeeUXz58/XsWPHnGuwhIaG6tNPP9U777yjnp4ezZ07V0uWLNEHH3ygp556ykOPCQAA/Jnb67CYiHVYgODBei1A4PDaOiwAYDc2MwSCE4EFgF+ZajNDd/cjAuA/CCwAfO5+gsVUK+LS+gIELgILAJ/zVrCYqvUFgH9ze5YQANyvsrxU58DZOxW9+aFaO/qUnhit4xuecetzV+UkMxAXCFAEFgA+Nz5Y3Dnrp7WjT5KcZwCQ6BICYLOx7qHX3r+sZTvOKDEmUpKUnhhtc80AmIQWFgC2GuseGhj6Wh29g0qMma3PdxTaXS0AhiGwALDN+O4gSZOOawEAAgsA24yfLTTZVGUAkBjDAsBGTEMGMF3sJQRgxtjXB8D9YC8hAD5xrwXgit78UClbalT05oc+rhmAQENgATAj1Y3tGhj6WjGzZ03apcOaKgA8hcACYEaq6tvUO3hLcyLCJu0OGltLhTVVANwvZgkBmJGpltcf4+7S+gAwGQbdAgAAWzDoFgAABBQCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHnsJAZixojc/VGtHn9ITo9k3CIBX0cICYMZaO/pczgDgLQQWAHfZeKhFqRU12nioZcpy6YnRLmcA8Ba6hIAgMt0unJrWTo1Yt89Zjz2kqvo2leWlalVOsks5uoEA+AotLEAQmawL584WlcL0BIU6bp+r6tvU0Tuoqvo2n9cXAMYQWIAgcmcXTnVju5btOKP//vU3LSqStKd0kdoqC7WndJHK8lKVGDNbZXmpttUbAByWZVl2V+J+9ff3Kzo6Wn19fYqKirK7OoDxqhvbVVXfpoGhr9U7eEuzZ4Vo+OtRFaYnaE/pIrurByBIuPP9zRgWIAiNdfPEzJ7lbD25c3wKAJiEwAIEobK81EkH0gKAiegSAgLAWBcPAQSAP3Hn+5tBt0AAuNdMnrHBtdWN7T6uGQB4BoEFCAD3msnD1GQA/o4xLEAAWJWTPGlXUHVjuwaGvlbM7FlMTQbgt2hhAfzMdJfNH1NV36bewVuaExHG+BYAfovAAniIu0FipsYvmz8dLPwGIBAQWAAPcTdIzNT4ZfOnY1VOss5u+R6tKwD8GoEF8BB3g4S7xmb6ZD32kHPZ/HvxVasPAHgbgQXwkPH773jDZDN9ppqy7KtWHwDwNgIL4CcmG4sy1ZRlb7f6AICvsNItYLDprGDLKrcA/BUr3QIBYqLWkzu7gBhUCyAYEFgAg03UDfTa+5fV0Tuo196/bGPNAMC3WOkWMNhUK9gCQDChhQXwMy8VPKnEmNl6qeBJu6sCAD7DoFsAAGALBt0CAYYF4AAEOwIL4AdYAA5AsCOwAH6ABeAABDvGsACG2XioRTWtnSpMT/DaMv8AYALGsAAGmmrPn/Ho/gGAuxFYAB+Z7uaFdP8AwN1mFFj27dunlJQURUZGKjs7W+fPn5+y/NGjR7VgwQJFRkYqLS1NJ0+edP7s1q1bevnll5WWlqY5c+YoISFBa9asUWcn/3eJwDLdzQu9veszAPgjtwPLkSNHVF5eru3bt6u5uVkZGRkqKCjQtWvXJix/7tw5lZaWav369WppaVFxcbGKi4t16dIlSdKf/vQnNTc365//+Z/V3Nys9957T5cvX1ZRUdH9PRlgmMn2/JksyAAAvuH2oNvs7GwtWbJEe/fulSSNjo4qKSlJGzZs0JYtW+4qX1JSooGBAZ04ccJ5LScnR5mZmdq/f/+Ev+Pjjz9WVlaW2tvb9eijj96zTgy6BQDA/3ht0O3w8LCampqUn5//zQeEhCg/P18NDQ0T3tPQ0OBSXpIKCgomLS9JfX19cjgciomJmfDnQ0ND6u/vdzkAAEDgciuw9PT0aGRkRHFxcS7X4+Li1NXVNeE9XV1dbpW/efOmXn75ZZWWlk6atiorKxUdHe08kpKS3HkMwAjTnTUEADBsltCtW7e0cuVKWZalqqqqSctVVFSor6/PeVy9etWHtQQ8Y7JZQwCAu4W5Uzg2NlahoaHq7u52ud7d3a34+PgJ74mPj59W+bGw0t7erjNnzkzZlxUREaGIiAh3qg4YpywvVVX1bQy2BYBpcKuFJTw8XIsXL1ZdXZ3z2ujoqOrq6pSbmzvhPbm5uS7lJam2ttal/FhY+d3vfqfTp09r7ty57lQL8EuTzRoCANzNrRYWSSovL9fatWv19NNPKysrS7t379bAwIDWrVsnSVqzZo0SExNVWVkpSdq0aZOWL1+uXbt2qbCwUIcPH9aFCxd04MABSbfDyt/+7d+qublZJ06c0MjIiHN8y0MPPaTw8HBPPSsAAPBTbgeWkpISffXVV9q2bZu6urqUmZmpU6dOOQfWXrlyRSEh3zTcLF26VAcPHtTWrVv1yiuvaP78+Tp27JgWLlwoSero6NDx48clSZmZmS6/65e//KXy8vJm+GgAACBQsPkhAACwBZsfAgCAgEJgAWaIdVQAwHcILMAMsY4KAPgOgQWYITYtBADfYdAtgtLGQy2qae1UYXqC9pQusrs6ABCUGHQL3ENNa6dGrNtnX9l4qEWpFTXaeKjFZ78TAAIFgQVBqTA9QaGO2+fput9BtnaEJAAIFAQWBKU9pYvUVlnoVnfQ/Q6ynUlIAgDc5vZKt0Cwut/NCveULmK8DADMEINuAQCALRh0CwAAAgqBBbgPzPwBAN8gsAD3gZk/AOAbBBbgPjDzBwB8g0G3AADAFgy6BQAAAYV1WBDwxu8blPXYQ861VFblJNtdNQDANNHCgoA3fmDs/a5WCwCwB4EFtrrf/Xmm89lPJUQ7B8aW5aUqMWb2jFerBQDYg8ACW3mzxWPssz/p7NOPn1+oPaWLtConWYuTH9T2/7rE2ikA4EcILLCVN1s8yvJSFeqQRiy5BCLWTgEA/8OgW9hqVU6y1wa/rspJ1i8+vqrWjj7NnRMu6XY3UXhYiG7eGmXtFADwI7SwwO9Ntjx+dWO7Wjv6JEmfdN4+V9W3afDWqBJiZrNzMgD4EQIL/N7xX9/u4jn+a9cunvHdQGOtKQy6BQD/RJcQAlZZXupda654swsKAOA9BBb4zPgF3NztjqlubHcJH+M/qygjwfnP4xFOACBw0CUEn7mf2Tl3Tn8e/1l7SheprbKQMSkAEMAILPCZ+9nZ+M6xJ2OfMWKJ9VQAIAiwWzP8zlh30Mj//5sb6pDaKgvtrRQAwG3s1oyAdmdYYT0VAAh8BBYYbaK9hsa6looyEhi7AgBBgi4hGG3ZjjPq6B1UYsxsnd3yPburAwDwILqEEBCqG9s1MPS1Zs8K0cDQ117Z0RkA4B9YhwXGqqpvU+/gLUnS4K1Rbf+vS5LE2ioAEIRoYYFxxvYGmjsnXKGOb66PWNJr71++a0wLACDwEVhgnLFZQJ909unHzy9UYsxsFWUkKDFmtiS5LCAHAAgOdAnBOIXp3yy1f+fy+uOX6AcABA9mCcFtd+7rAwDATDBLCF51574+AAB4G4EFbrtzXx9PGBtoy75AAICJ0CUEI6RsqXH+8+c72BcIAIIBXUIwHi0qAAB3EFhgi7GpyzWtnZKk9MRolzMAAOMxrRm2GD91WZKOb3jG67+T2U0A4L8YwwKfMCEssJEiAJiFMSwwjglTob0xuwkA4Bt0CcEnyvJSbV+h9s5VcwEA/oMuIdhi46EW5xiWPaWL7K4OAMAGdAnBeHfOEgIAYCoEFkxLdWO7lu04o+rGdo98XmF6gkIdcs4SAgBgKnQJYVqYYQMA8DS6hOBxzLABANiJFha4zYQ1VQAA/o8WFniVCWuqAACCC4EFbqN7CADga3QJYUpFb36o1o4+pSdG+2S/HwBA8KBLCB7T2tHncgYAwA4EFkwpPTHa5QwAgB3YSwgTGr90Pl1BAAC7zaiFZd++fUpJSVFkZKSys7N1/vz5KcsfPXpUCxYsUGRkpNLS0nTy5EmXn7/33nt67rnnNHfuXDkcDl28eHEm1YIHjS2df/zXnR5b3RYAgJlyO7AcOXJE5eXl2r59u5qbm5WRkaGCggJdu3ZtwvLnzp1TaWmp1q9fr5aWFhUXF6u4uFiXLl1ylhkYGNAzzzyjnTt3zvxJ4FHjl8xn+jIAwG5uzxLKzs7WkiVLtHfvXknS6OiokpKStGHDBm3ZsuWu8iUlJRoYGNCJEyec13JycpSZman9+/e7lP3888/12GOPqaWlRZmZmdOuE7OEvIMF4gAA3uTO97dbY1iGh4fV1NSkiooK57WQkBDl5+eroaFhwnsaGhpUXl7ucq2goEDHjh1z51e7GBoa0tDQkPPP/f39M/6sYDZ+nMqe0kXO68t21Kmj96YSYyJ1dsuzNtYQAIDb3OoS6unp0cjIiOLi4lyux8XFqaura8J7urq63Co/HZWVlYqOjnYeSUlJM/6sYDY2TqWmtdPlekfvTZczAAB288tpzRUVFerr63MeV69etbtKfqkwPUGhDtfxKpKUGBPpcgYAwG5udQnFxsYqNDRU3d3dLte7u7sVHx8/4T3x8fFulZ+OiIgIRUREzPh+3LandJFLV9AYuoEAAKZxq4UlPDxcixcvVl1dnfPa6Oio6urqlJubO+E9ubm5LuUlqba2dtLyAAAAd3J74bjy8nKtXbtWTz/9tLKysrR7924NDAxo3bp1kqQ1a9YoMTFRlZWVkqRNmzZp+fLl2rVrlwoLC3X48GFduHBBBw4ccH7mH/7wB125ckWdnbfHUly+fFnS7daZ+2mJwdQmG3QLAIBp3B7DUlJSotdee03btm1TZmamLl68qFOnTjkH1l65ckVffvmls/zSpUt18OBBHThwQBkZGfqP//gPHTt2TAsXLnSWOX78uBYtWqTCwkJJ0gsvvKBFixbdNe0ZnjXZoFsAAEzDbs1BjBYWAICd3Pn+JrAAAABbuPP97ZfTmgEAQHAhsPi56sZ2Ldtxhg0KAQABjcDi56rq29TROzjlBoVFb36olC01KnrzQx/WDAAAzyGw+LmyvFQlxsxWWV7qpGVaO/pczgAA+BsCix+b7m7K6YnRLmcAAPyN2wvHwRzju4OmCizHNzzjw1oBAOB5tLD4qerGdg0Mfa2Y2bOm7A4CACAQ0MLip37835/o1oilWaGOKVtXAAAIBLSw+KlbI5bLGQCAQEZg8VMMpAUABBO6hAxW9OaHzqnIRRmu+/0wkBYAEExoYTHY+HVTjv+aHZUBAMGLwGKwO7t7WH4fABCsCCwGO77hGZfQMtXy+wAABDICi+GOb3hG/6944T2X3wcAIJA5LMvy+3mx/f39io6OVl9fn6KiouyuDgAAmAZ3vr9pYQEAAMYjsAAAAOMRWDyourFdy3acYTYPAAAeRmDxoPG7J7tj2Y46pWyp0bIddV6qGQAA/o3A4kFleakzms3T0XvT5QwAAFyxNL8HrcpJntHOyYkxkerovanEmEgv1AoAAP9HYLFBdWO7qurbVJaXqlU5yTq75Vm7qwQAgNHoErLBTMe6AAAQrAgsNpjpWBcAAIIVK936wMZDLapp7VRheoL2lC6yuzoAABiBlW4NU9PaqRHr9hkAALiPwOIDhekJCnXcPgMAAPcxS8gDxmb9LE5+UE3t/+uc/SO5dgdlPfaQlu04o7K8VP3i46tq7ehTemK0jm94xuYnAADAbAQWDxib9dPVN6gR6/afxwLL+O6gpvb/dc4O6ugdlCS1dvTZWXUAAPwCXUIeMDbrpzA94a7ZP+O7g8bPDkpPjJYk5xkAAEyOWUIAAMAWzBICAAABhcACAACMR2ABAADGI7DcQ3Vju5btOKONh1q0bMcZVTe2T/jzO68DAADPIbDcw9gU5JrWzgk3LGQjQwAAvI/Acg9TTVke/3M2MgQAwHuY1gwAAGzBtGYAABBQCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGC/M7gp4wtiG0/39/TbXBAAATNfY9/bY9/hUAiKw3LhxQ5KUlJRkc00AAIC7bty4oejo6CnLOKzpxBrDjY6OqrOzUw888IAcDofd1fFL/f39SkpK0tWrVxUVFWV3dQIS79j7eMfexzv2vmB6x5Zl6caNG0pISFBIyNSjVAKihSUkJESPPPKI3dUICFFRUQH/L4jdeMfexzv2Pt6x9wXLO75Xy8oYBt0CAADjEVgAAIDxCCyQJEVERGj79u2KiIiwuyoBi3fsfbxj7+Mdex/veGIBMegWAAAENlpYAACA8QgsAADAeAQWAABgPAILAAAwHoElgO3bt08pKSmKjIxUdna2zp8/P2X5o0ePasGCBYqMjFRaWppOnjzp8vP33ntPzz33nObOnSuHw6GLFy96sfb+wZPv+NatW3r55ZeVlpamOXPmKCEhQWvWrFFnZ6e3H8Nonv57/Oqrr2rBggWaM2eOHnzwQeXn5+ujjz7y5iMYzdPvd7wf/vCHcjgc2r17t4dr7V88/Y5/8IMfyOFwuBwrVqzw5iOYwUJAOnz4sBUeHm799Kc/tT755BPrxRdftGJiYqzu7u4Jy589e9YKDQ21/vVf/9X6zW9+Y23dutWaNWuW9T//8z/OMu+++6714x//2HrrrbcsSVZLS4uPnsZMnn7Hvb29Vn5+vnXkyBHr008/tRoaGqysrCxr8eLFvnwso3jj7/HPf/5zq7a21mpra7MuXbpkrV+/3oqKirKuXbvmq8cyhjfe75j33nvPysjIsBISEqzXX3/dy09iLm+847Vr11orVqywvvzyS+fxhz/8wVePZBsCS4DKysqy/vEf/9H555GRESshIcGqrKycsPzKlSutwsJCl2vZ2dnWP/zDP9xV9ve//z2BxfLuOx5z/vx5S5LV3t7umUr7GV+8476+PkuSdfr0ac9U2o946/1+8cUXVmJionXp0iUrOTk5qAOLN97x2rVrreeff94r9TUZXUIBaHh4WE1NTcrPz3deCwkJUX5+vhoaGia8p6GhwaW8JBUUFExaPtj56h339fXJ4XAoJibGI/X2J754x8PDwzpw4ICio6OVkZHhucr7AW+939HRUa1evVqbN2/WU0895Z3K+wlv/h2ur6/XvHnz9OSTT6qsrEzXr1/3/AMYhsASgHp6ejQyMqK4uDiX63Fxcerq6prwnq6uLrfKBztfvOObN2/q5ZdfVmlpaVBsgHYnb77jEydO6Fvf+pYiIyP1+uuvq7a2VrGxsZ59AMN56/3u3LlTYWFh2rhxo+cr7We89Y5XrFihd999V3V1ddq5c6d+9atf6fvf/75GRkY8/xAGCYjdmoFAc+vWLa1cuVKWZamqqsru6gSc7373u7p48aJ6enr01ltvaeXKlfroo480b948u6vm15qamvTGG2+oublZDofD7uoErBdeeMH5z2lpaUpPT1dqaqrq6+v17LPP2lgz76KFJQDFxsYqNDRU3d3dLte7u7sVHx8/4T3x8fFulQ923nzHY2Glvb1dtbW1Qdm6Inn3Hc+ZM0dPPPGEcnJy9PbbbyssLExvv/22Zx/AcN54vx988IGuXbumRx99VGFhYQoLC1N7e7t+9KMfKSUlxSvPYTJf/bf48ccfV2xsrD777LP7r7TBCCwBKDw8XIsXL1ZdXZ3z2ujoqOrq6pSbmzvhPbm5uS7lJam2tnbS8sHOW+94LKz87ne/0+nTpzV37lzvPIAf8OXf49HRUQ0NDd1/pf2IN97v6tWr1draqosXLzqPhIQEbd68We+//773HsZQvvo7/MUXX+j69et6+OGHPVNxU9k96hfecfjwYSsiIsL62c9+Zv3mN7+x/v7v/96KiYmxurq6LMuyrNWrV1tbtmxxlj979qwVFhZmvfbaa9Zvf/tba/v27XdNpbt+/brV0tJi1dTUWJKsw4cPWy0tLdaXX37p8+czgaff8fDwsFVUVGQ98sgj1sWLF12mLA4NDdnyjHbz9Dv+4x//aFVUVFgNDQ3W559/bl24cMFat26dFRERYV26dMmWZ7STN/47cadgnyXk6Xd848YN66WXXrIaGhqs3//+99bp06et73znO9b8+fOtmzdv2vKMvkJgCWBvvvmm9eijj1rh4eFWVlaW1djY6PzZ8uXLrbVr17qU/8UvfmH9+Z//uRUeHm499dRTVk1NjcvP//3f/92SdNexfft2HzyNmTz5jsemi090/PKXv/TRE5nHk+94cHDQ+pu/+RsrISHBCg8Ptx5++GGrqKjIOn/+vK8exzie/u/EnYI9sFiWZ9/xn/70J+u5556zvv3tb1uzZs2ykpOTrRdffNEZgAKZw7Isy562HQAAgOlhDAsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxvs/J07oijOKd2kAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.scatter(pantheon_data[\"zCMB_SN\"], pantheon_data[\"zCMB_Group\"], s=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tully-Fisher galaxies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SFI++ galaxies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf_folder = \"/mnt/extraspace/rstiskalek/catalogs/PV/PV_Supranta/tf\"\n",
    "\n",
    "names = [\"RA\", \"DEC\", \"z_CMB\", \"mag\", \"eta\", \"e_mag\", \"e_eta\"]\n",
    "\n",
    "dtype = [(n, np.float32) for n in names]\n",
    "data = np.genfromtxt(join(tf_folder, \"sfi_gals_tf.csv\"), delimiter=\",\", skip_header=1,\n",
    "                     usecols=[2 + n for n in range(len(names))])\n",
    "\n",
    "sfi_gals = np.empty(len(data), dtype=dtype)\n",
    "for i, n in enumerate(names):\n",
    "    sfi_gals[n] = data[:, i]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SFI++ galaxies masked"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "names = [\"RA\", \"DEC\", \"z_CMB\", \"mag\", \"eta\", \"e_mag\", \"e_eta\"]\n",
    "\n",
    "dtype = [(n, np.float32) for n in names]\n",
    "data = np.genfromtxt(join(tf_folder, \"sfi_gals_tf_masked.csv\"), delimiter=\",\", skip_header=1,\n",
    "                     usecols=[2 + n for n in range(len(names))])\n",
    "\n",
    "sfi_gals_masked = np.empty(len(data), dtype=dtype)\n",
    "for i, n in enumerate(names):\n",
    "    sfi_gals_masked[n] = data[:, i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SFI++ groups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "names = [\"RA\", \"DEC\", \"z_CMB\", \"r_hMpc\", \"e_r_hMpc\"]\n",
    "\n",
    "dtype = [(n, np.float32) for n in names]\n",
    "data = np.genfromtxt(join(tf_folder, \"sfi_grps.csv\"), delimiter=\",\", skip_header=1,\n",
    "                     usecols=[1 + n for n in range(len(names))])\n",
    "\n",
    "sfi_groups = np.empty(len(data), dtype=dtype)\n",
    "for i, n in enumerate(names):\n",
    "    sfi_groups[n] = data[:, i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cross $\\texttt{SFI++ galaxies}$ and $\\texttt{SFI++ groups}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "[$3^\\circ40{}^\\prime57.6210022{}^{\\prime\\prime}$\n",
       " $0^\\circ45{}^\\prime09.98210907{}^{\\prime\\prime}$\n",
       " $1^\\circ47{}^\\prime48.02945137{}^{\\prime\\prime}$ ...\n",
       " $0^\\circ59{}^\\prime48.39719296{}^{\\prime\\prime}$\n",
       " $1^\\circ27{}^\\prime20.76476097{}^{\\prime\\prime}$\n",
       " $2^\\circ01{}^\\prime03.90094757{}^{\\prime\\prime}$]"
      ],
      "text/plain": [
       "<Angle [3.6826725, 0.7527728, 1.7966748, ..., 0.996777 , 1.455768 ,\n",
       "        2.0177503] deg>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "match_coord = SkyCoord(sfi_gals[\"RA\"], sfi_gals[\"DEC\"], unit=\"deg\")\n",
    "catalogue_coord = SkyCoord(sfi_groups[\"RA\"], sfi_groups[\"DEC\"],\n",
    "                           unit=\"deg\")\n",
    "\n",
    "idxs, sep2d, sep3d = match_coordinates_sky(match_coord, catalogue_coord, nthneighbor=1)\n",
    "sep2d.to(u.degree)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAr9UlEQVR4nO3deXTTdb7/8VdKF5a2qWXpAoWWfd+hVh1EqLQFUaT3Kg7XU5haHWxR6EGlKpvjPSguMDooeq4Djld0ZFQQVLxYVrHsmzhQpQMXmC5AmTYUpa3t9/eHl/wmw9raNOmnz8c5OYfk+03yTr8nznO+328Sm2VZlgAAANDg+Xh6AAAAANQNwg4AAMAQhB0AAIAhCDsAAABDEHYAAACGIOwAAAAMQdgBAAAYgrADAAAwhK+nB/AG1dXVys/PV1BQkGw2m6fHAQAAcLIsS+fOnVNkZKR8fK6+T46wk5Sfn6+oqChPjwEAAHBFJ06cULt27a66DmEnKSgoSNLPf7Dg4GAPTwMAAPD/ORwORUVFOXvlagg7yXn4NTg4mLADAABe6XpOF+PDEwAAAIYg7AAAAAxB2AEAABiCsAMAADAEYQcAAGAIwg4AAMAQhB0AAIAhCDsAAABDEHYAAACGIOwAAAAMQdgBAAAYgrADAAAwBGEHAABgCMIOAADAEIQdAACAIQg7AAAAQ/h6egDUneiZn9Zo/WPPjXHTJAAAwBPYYwcAAGAIwg4AAMAQhB0AAIAhOMeuEeOcPAAAzMIeOwAAAEMQdgAAAIYg7AAAAAxB2AEAABiCsAMAADAEYQcAAGAIwg4AAMAQhB0AAIAhCDsAAABDEHYAAACGIOwAAAAMQdgBAAAYgrADAAAwBGEHAABgCMIOAADAEIQdAACAIQg7AAAAQxB2AAAAhvBo2M2fP19DhgxRUFCQ2rRpo3Hjxik3N9dlnQsXLig9PV0tW7ZUYGCgkpOTVVRU5LLO8ePHNWbMGDVv3lxt2rTRY489pp9++qk+XwoAAIDHeTTsNm3apPT0dG3btk3r1q1TZWWlRo0apfPnzzvXmT59ulavXq0VK1Zo06ZNys/P1/jx453Lq6qqNGbMGFVUVOjrr7/W22+/rWXLlmn27NmeeEkAAAAeY7Msy/L0EBedPn1abdq00aZNmzRs2DCVlpaqdevWWr58uf7t3/5NknT48GH16NFDOTk5uvHGG/X555/rjjvuUH5+vsLCwiRJS5Ys0RNPPKHTp0/L39//ms/rcDhkt9tVWlqq4OBgt75Gd4qe+albH//Yc2Pc+vgAAOBSNekUrzrHrrS0VJIUGhoqSdq9e7cqKysVHx/vXKd79+5q3769cnJyJEk5OTnq06ePM+okKSEhQQ6HQ99+++1ln6e8vFwOh8PlAgAA0NB5TdhVV1dr2rRpuvnmm9W7d29JUmFhofz9/RUSEuKyblhYmAoLC53r/HPUXVx+cdnlzJ8/X3a73XmJioqq41cDAABQ/7wm7NLT03Xw4EG9//77bn+urKwslZaWOi8nTpxw+3MCAAC4m6+nB5CkjIwMrVmzRps3b1a7du2ct4eHh6uiokIlJSUue+2KiooUHh7uXGfHjh0uj3fxU7MX1/lXAQEBCggIqONXAQAA4Fke3WNnWZYyMjL08ccfa/369YqJiXFZPmjQIPn5+Sk7O9t5W25uro4fP664uDhJUlxcnL755hudOnXKuc66desUHBysnj171s8LAQAA8AIe3WOXnp6u5cuXa9WqVQoKCnKeE2e329WsWTPZ7XalpqYqMzNToaGhCg4O1tSpUxUXF6cbb7xRkjRq1Cj17NlT999/vxYsWKDCwkI9/fTTSk9PZ68cAABoVDwadq+//rokafjw4S63L126VJMmTZIkLVy4UD4+PkpOTlZ5ebkSEhL02muvOddt0qSJ1qxZoylTpiguLk4tWrRQSkqKnnnmmfp6GQAAAF7Bq77HzlP4Hrvrw/fYAQBQ/xrs99gBAACg9gg7AAAAQxB2AAAAhiDsAAAADEHYAQAAGIKwAwAAMARhBwAAYAjCDgAAwBCEHQAAgCEIOwAAAEMQdgAAAIYg7AAAAAxB2AEAABiCsAMAADAEYQcAAGAIX08PgIYjeuanNVr/2HNj3DQJAAC4HPbYAQAAGIKwAwAAMARhBwAAYAjCDgAAwBCEHQAAgCEIOwAAAEMQdgAAAIYg7AAAAAxB2AEAABiCsAMAADAEYQcAAGAIwg4AAMAQhB0AAIAhCDsAAABDEHYAAACGIOwAAAAMQdgBAAAYgrADAAAwBGEHAABgCMIOAADAEIQdAACAIQg7AAAAQxB2AAAAhiDsAAAADEHYAQAAGIKwAwAAMARhBwAAYAjCDgAAwBCEHQAAgCEIOwAAAEMQdgAAAIYg7AAAAAxB2AEAABiCsAMAADAEYQcAAGAIwg4AAMAQhB0AAIAhCDsAAABDEHYAAACGIOwAAAAM4evpAYCLomd+WqP1jz03xk2TAADQMLHHDgAAwBCEHQAAgCEIOwAAAEMQdgAAAIYg7AAAAAxB2AEAABiCsAMAADAEYQcAAGAIwg4AAMAQhB0AAIAhCDsAAABDEHYAAACGIOwAAAAMQdgBAAAYgrADAAAwBGEHAABgCMIOAADAEIQdAACAIQg7AAAAQ/h6egBcWfTMTz09AgAAaEDYYwcAAGAIwg4AAMAQHg27zZs3a+zYsYqMjJTNZtPKlStdlk+aNEk2m83lkpiY6LLO2bNnNXHiRAUHByskJESpqakqKyurx1cBAADgHTx6jt358+fVr18//eY3v9H48eMvu05iYqKWLl3qvB4QEOCyfOLEiSooKNC6detUWVmpyZMn68EHH9Ty5cvdOjs8r6bnIB57boybJgEAwDt4NOySkpKUlJR01XUCAgIUHh5+2WWHDh3S2rVrtXPnTg0ePFiS9Oqrr2r06NF68cUXFRkZWeczAwAAeCuvP8du48aNatOmjbp166YpU6aouLjYuSwnJ0chISHOqJOk+Ph4+fj4aPv27Vd8zPLycjkcDpcLAABAQ+fVYZeYmKg//elPys7O1vPPP69NmzYpKSlJVVVVkqTCwkK1adPG5T6+vr4KDQ1VYWHhFR93/vz5stvtzktUVJRbXwcAAEB98OrvsZswYYLz33369FHfvn3VqVMnbdy4USNHjqz142ZlZSkzM9N53eFwEHcAAKDB8+o9dv+qY8eOatWqlY4cOSJJCg8P16lTp1zW+emnn3T27Nkrnpcn/XzeXnBwsMsFAACgoWtQYXfy5EkVFxcrIiJCkhQXF6eSkhLt3r3buc769etVXV2t2NhYT40JAADgER49FFtWVubc+yZJR48e1b59+xQaGqrQ0FDNmzdPycnJCg8PV15enh5//HF17txZCQkJkqQePXooMTFRaWlpWrJkiSorK5WRkaEJEybwiVgAANDoeHSP3a5duzRgwAANGDBAkpSZmakBAwZo9uzZatKkiQ4cOKA777xTXbt2VWpqqgYNGqQtW7a4fJfdu+++q+7du2vkyJEaPXq0brnlFr355pueekkAAAAe49E9dsOHD5dlWVdc/sUXX1zzMUJDQ/kyYgAAADWwc+wAAABwZYQdAACAIQg7AAAAQxB2AAAAhvDqX55AwxY981NPjwAAQKPCHjsAAABDEHYAAACGIOwAAAAMQdgBAAAYgrADAAAwBGEHAABgCMIOAADAEIQdAACAIQg7AAAAQxB2AAAAhiDsAAAADEHYAQAAGIKwAwAAMARhBwAAYAjCDgAAwBCEHQAAgCEIOwAAAEMQdgAAAIYg7AAAAAxB2AEAABiCsAMAADAEYQcAAGAIwg4AAMAQhB0AAIAhCDsAAABDEHYAAACGIOwAAAAMQdgBAAAYgrADAAAwBGEHAABgCMIOAADAEIQdAACAIXw9PUBjEj3zU0+PAAAADMYeOwAAAEPUKuz27Nmjb775xnl91apVGjdunJ588klVVFTU2XAAAAC4frU6FPvQQw9p5syZ6tOnj/72t79pwoQJuvvuu7VixQr98MMPWrRoUR2PCfxyNT0Ufuy5MW6aBAAA96jVHrvvvvtO/fv3lyStWLFCw4YN0/Lly7Vs2TJ9+OGHdTkfAAAArlOtws6yLFVXV0uSvvzyS40ePVqSFBUVpTNnztTddAAAALhutQq7wYMH69lnn9U777yjTZs2acyYnw9ZHT16VGFhYXU6IAAAAK5PrcJu4cKF2rNnjzIyMvTUU0+pc+fOkqS//OUvuummm+p0QAAAAFyfWn14ol+/fi6fir3ohRdekK8vX40HAADgCbXaY9exY0cVFxdfcvuFCxfUtWvXXzwUAAAAaq5WYXfs2DFVVVVdcnt5eblOnjz5i4cCAABAzdXouOknn3zi/PcXX3whu93uvF5VVaXs7GzFxMTU3XQAAAC4bjUKu3HjxkmSbDabUlJSXJb5+fkpOjpaL730Up0NBwAAgOtXo7C7+N11MTEx2rlzp1q1auWWoQAAAFBztfoI69GjR+t6DsDr8BNkAICGptbfTZKdna3s7GydOnXKuSfvoj/+8Y+/eDAAAADUTK3Cbt68eXrmmWc0ePBgRUREyGaz1fVcAAAAqKFahd2SJUu0bNky3X///XU9DwAAAGqpVt9jV1FRwU+HAQAAeJlahd0DDzyg5cuX1/UsAAAA+AVqdSj2woULevPNN/Xll1+qb9++8vPzc1n+8ssv18lwAAAAuH61CrsDBw6of//+kqSDBw+6LOODFAAAAJ5Rq7DbsGFDXc8BAACAX6hW59gBAADA+9Rqj91tt9121UOu69evr/VAAAAAqJ1ahd3F8+suqqys1L59+3Tw4EGlpKTUxVwAAACooVqF3cKFCy97+9y5c1VWVvaLBgIAAEDt1Ok5dv/xH//B78QCAAB4SJ2GXU5Ojpo2bVqXDwkAAIDrVKtDsePHj3e5blmWCgoKtGvXLs2aNatOBgMAAEDN1Crs7Ha7y3UfHx9169ZNzzzzjEaNGlUngwEAAKBmahV2S5cures5AAAA8AvVKuwu2r17tw4dOiRJ6tWrlwYMGFAnQwEAAKDmahV2p06d0oQJE7Rx40aFhIRIkkpKSnTbbbfp/fffV+vWretyRgAAAFyHWn0qdurUqTp37py+/fZbnT17VmfPntXBgwflcDj0yCOP1PWMAAAAuA612mO3du1affnll+rRo4fztp49e2rx4sV8eAIAAMBDarXHrrq6Wn5+fpfc7ufnp+rq6l88FAAAAGquVmE3YsQIPfroo8rPz3fe9ve//13Tp0/XyJEj62w4AAAAXL9ahd0f/vAHORwORUdHq1OnTurUqZNiYmLkcDj06quv1vWMAAAAuA61OscuKipKe/bs0ZdffqnDhw9Lknr06KH4+Pg6HQ4AAADXr0Z77NavX6+ePXvK4XDIZrPp9ttv19SpUzV16lQNGTJEvXr10pYtW9w1KwAAAK6iRmG3aNEipaWlKTg4+JJldrtdDz30kF5++eU6Gw4AAADXr0Zht3//fiUmJl5x+ahRo7R79+7rfrzNmzdr7NixioyMlM1m08qVK12WW5al2bNnKyIiQs2aNVN8fLy+//57l3XOnj2riRMnKjg4WCEhIUpNTVVZWVlNXhYAAIARahR2RUVFl/2ak4t8fX11+vTp63688+fPq1+/flq8ePFlly9YsECvvPKKlixZou3bt6tFixZKSEjQhQsXnOtMnDhR3377rdatW6c1a9Zo8+bNevDBB6//RQEAABiiRh+eaNu2rQ4ePKjOnTtfdvmBAwcUERFx3Y+XlJSkpKSkyy6zLEuLFi3S008/rbvuukuS9Kc//UlhYWFauXKlJkyYoEOHDmnt2rXauXOnBg8eLEl69dVXNXr0aL344ouKjIysycsDAABo0Gq0x2706NGaNWuWyx6zi3788UfNmTNHd9xxR50MdvToURUWFrp80tZutys2NlY5OTmSpJycHIWEhDijTpLi4+Pl4+Oj7du318kcAAAADUWN9tg9/fTT+uijj9S1a1dlZGSoW7dukqTDhw9r8eLFqqqq0lNPPVUngxUWFkqSwsLCXG4PCwtzLissLFSbNm1clvv6+io0NNS5zuWUl5ervLzced3hcNTJzAAAAJ5Uo7ALCwvT119/rSlTpigrK0uWZUmSbDabEhIStHjx4ktCzBvNnz9f8+bN8/QYAAAAdarGX1DcoUMHffbZZ/rHP/6hI0eOyLIsdenSRTfccEOdDhYeHi7p5w9s/PN5e0VFRerfv79znVOnTrnc76efftLZs2ed97+crKwsZWZmOq87HA5FRUXV4fQAAAD1r1Y/KSZJN9xwg4YMGaKhQ4fWedRJUkxMjMLDw5Wdne28zeFwaPv27YqLi5MkxcXFqaSkxOUrVtavX6/q6mrFxsZe8bEDAgIUHBzscgEAAGjoavWTYnWlrKxMR44ccV4/evSo9u3bp9DQULVv317Tpk3Ts88+qy5duigmJkazZs1SZGSkxo0bJ+nnnzFLTExUWlqalixZosrKSmVkZGjChAl8IhYAADQ6Hg27Xbt26bbbbnNev3h4NCUlRcuWLdPjjz+u8+fP68EHH1RJSYluueUWrV27Vk2bNnXe591331VGRoZGjhwpHx8fJScn65VXXqn31wIAAOBpNuviJyAaMYfDIbvdrtLSUrcelo2e+anbHhued+y5MZ4eAQBgoJp0Sq3PsQMAAIB3IewAAAAMQdgBAAAYgrADAAAwBGEHAABgCMIOAADAEIQdAACAIQg7AAAAQxB2AAAAhiDsAAAADEHYAQAAGIKwAwAAMARhBwAAYAjCDgAAwBCEHQAAgCF8PT0AYIromZ/WaP1jz41x0yQAgMaKPXYAAACGIOwAAAAMQdgBAAAYgrADAAAwBGEHAABgCMIOAADAEIQdAACAIQg7AAAAQxB2AAAAhiDsAAAADEHYAQAAGIKwAwAAMARhBwAAYAjCDgAAwBCEHQAAgCEIOwAAAEMQdgAAAIYg7AAAAAxB2AEAABiCsAMAADAEYQcAAGAIwg4AAMAQhB0AAIAhCDsAAABD+Hp6AKCxip75aY3WP/bcGDdNAgAwBXvsAAAADEHYAQAAGIKwAwAAMARhBwAAYAjCDgAAwBCEHQAAgCEIOwAAAEMQdgAAAIYg7AAAAAxB2AEAABiCsAMAADAEYQcAAGAIwg4AAMAQhB0AAIAhCDsAAABDEHYAAACGIOwAAAAMQdgBAAAYgrADAAAwBGEHAABgCMIOAADAEIQdAACAIQg7AAAAQxB2AAAAhiDsAAAADEHYAQAAGIKwAwAAMARhBwAAYAjCDgAAwBCEHQAAgCEIOwAAAEMQdgAAAIYg7AAAAAxB2AEAABiCsAMAADCEr6cHAHB9omd+WqP1jz03xk2TAAC8FXvsAAAADEHYAQAAGMKrw27u3Lmy2Wwul+7duzuXX7hwQenp6WrZsqUCAwOVnJysoqIiD04MAADgOV4ddpLUq1cvFRQUOC9fffWVc9n06dO1evVqrVixQps2bVJ+fr7Gjx/vwWkBAAA8x+s/POHr66vw8PBLbi8tLdVbb72l5cuXa8SIEZKkpUuXqkePHtq2bZtuvPHG+h4VAADAo7x+j93333+vyMhIdezYURMnTtTx48clSbt371ZlZaXi4+Od63bv3l3t27dXTk6Op8YFAADwGK/eYxcbG6tly5apW7duKigo0Lx58/SrX/1KBw8eVGFhofz9/RUSEuJyn7CwMBUWFl71ccvLy1VeXu687nA43DE+AABAvfLqsEtKSnL+u2/fvoqNjVWHDh30wQcfqFmzZrV+3Pnz52vevHl1MSIAAIDX8PpDsf8sJCREXbt21ZEjRxQeHq6KigqVlJS4rFNUVHTZc/L+WVZWlkpLS52XEydOuHFqAACA+tGgwq6srEx5eXmKiIjQoEGD5Ofnp+zsbOfy3NxcHT9+XHFxcVd9nICAAAUHB7tcAAAAGjqvPhQ7Y8YMjR07Vh06dFB+fr7mzJmjJk2a6L777pPdbldqaqoyMzMVGhqq4OBgTZ06VXFxcXwiFgAANEpeHXYnT57Ufffdp+LiYrVu3Vq33HKLtm3bptatW0uSFi5cKB8fHyUnJ6u8vFwJCQl67bXXPDw1AACAZ9gsy7I8PYSnORwO2e12lZaWuvWwbE1/xB2oT8eeG+PpEQAAl1GTTmlQ59gBAADgygg7AAAAQxB2AAAAhiDsAAAADEHYAQAAGIKwAwAAMARhBwAAYAjCDgAAwBCEHQAAgCEIOwAAAEN49W/FAqg/Nf3JO36CDAC8D3vsAAAADEHYAQAAGIKwAwAAMARhBwAAYAjCDgAAwBCEHQAAgCEIOwAAAEMQdgAAAIYg7AAAAAxB2AEAABiCsAMAADAEYQcAAGAIwg4AAMAQhB0AAIAhCDsAAABDEHYAAACGIOwAAAAMQdgBAAAYgrADAAAwBGEHAABgCF9PDwCgYYqe+WmN1j/23Bg3TQIAuIg9dgAAAIYg7AAAAAxB2AEAABiCsAMAADAEYQcAAGAIwg4AAMAQhB0AAIAhCDsAAABDEHYAAACGIOwAAAAMQdgBAAAYgrADAAAwBGEHAABgCMIOAADAEIQdAACAIQg7AAAAQxB2AAAAhiDsAAAADEHYAQAAGMLX0wMAaByiZ37q1sc/9twYtz4+ADQE7LEDAAAwBGEHAABgCMIOAADAEJxjB6BRquk5f5zDB6AhYI8dAACAIdhjBwDXgT18ABoC9tgBAAAYgrADAAAwBIdiARjB3V+ADAANAXvsAAAADEHYAQAAGIKwAwAAMARhBwAAYAjCDgAAwBCEHQAAgCEIOwAAAEMQdgAAAIbgC4oBoAHit2sBXA577AAAAAxB2AEAABiCsAMAADAE59gBgBeo6TlzAHA57LEDAAAwBGEHAABgCMIOAADAEJxjBwD4xfhePcA72CzLsjw9hKc5HA7Z7XaVlpYqODjYbc/DydEAUDuEIBqzmnSKMYdiFy9erOjoaDVt2lSxsbHasWOHp0cCAACoV0Yciv3zn/+szMxMLVmyRLGxsVq0aJESEhKUm5urNm3aeHo8AMAv5O4jHuwRhCmM2GP38ssvKy0tTZMnT1bPnj21ZMkSNW/eXH/84x89PRoAAEC9afB77CoqKrR7925lZWU5b/Px8VF8fLxycnI8OBkAAN6JD7tcXUP++zT4sDtz5oyqqqoUFhbmcntYWJgOHz582fuUl5ervLzceb20tFTSzycnulN1+Q9ufXwAQO24+7//3qam/3vE3+fq3P33ufj41/N51wYfdrUxf/58zZs375Lbo6KiPDANAMDT7Is8PYF34+9zdfX19zl37pzsdvtV12nwYdeqVSs1adJERUVFLrcXFRUpPDz8svfJyspSZmam83p1dbXOnj2rli1bymaz1fmMDodDUVFROnHihFu/TgU1w3bxTmwX78M28U5sF+/kju1iWZbOnTunyMjIa67b4MPO399fgwYNUnZ2tsaNGyfp51DLzs5WRkbGZe8TEBCggIAAl9tCQkLcPKkUHBzMm88LsV28E9vF+7BNvBPbxTvV9Xa51p66ixp82ElSZmamUlJSNHjwYA0dOlSLFi3S+fPnNXnyZE+PBgAAUG+MCLt7771Xp0+f1uzZs1VYWKj+/ftr7dq1l3ygAgAAwGRGhJ0kZWRkXPHQq6cFBARozpw5lxz+hWexXbwT28X7sE28E9vFO3l6u/BbsQAAAIYw4pcnAAAAQNgBAAAYg7ADAAAwBGFXDxYvXqzo6Gg1bdpUsbGx2rFjh6dHatTmzp0rm83mcunevbunx2pUNm/erLFjxyoyMlI2m00rV650WW5ZlmbPnq2IiAg1a9ZM8fHx+v777z0zbCNyre0yadKkS947iYmJnhm2kZg/f76GDBmioKAgtWnTRuPGjVNubq7LOhcuXFB6erpatmypwMBAJScnX/Kl/ahb17Ndhg8ffsn75be//a3bZyPs3OzPf/6zMjMzNWfOHO3Zs0f9+vVTQkKCTp065enRGrVevXqpoKDAefnqq688PVKjcv78efXr10+LFy++7PIFCxbolVde0ZIlS7R9+3a1aNFCCQkJunDhQj1P2rhca7tIUmJiost757333qvHCRufTZs2KT09Xdu2bdO6detUWVmpUaNG6fz58851pk+frtWrV2vFihXatGmT8vPzNX78eA9Obb7r2S6SlJaW5vJ+WbBggfuHs+BWQ4cOtdLT053Xq6qqrMjISGv+/PkenKpxmzNnjtWvXz9Pj4H/I8n6+OOPnderq6ut8PBw64UXXnDeVlJSYgUEBFjvvfeeByZsnP51u1iWZaWkpFh33XWXR+bBz06dOmVJsjZt2mRZ1s/vDT8/P2vFihXOdQ4dOmRJsnJycjw1ZqPzr9vFsizr1ltvtR599NF6n4U9dm5UUVGh3bt3Kz4+3nmbj4+P4uPjlZOT48HJ8P333ysyMlIdO3bUxIkTdfz4cU+PhP9z9OhRFRYWurxv7Ha7YmNjed94gY0bN6pNmzbq1q2bpkyZouLiYk+P1KiUlpZKkkJDQyVJu3fvVmVlpcv7pXv37mrfvj3vl3r0r9vlonfffVetWrVS7969lZWVpR9++MHtsxjzBcXe6MyZM6qqqrrkFzDCwsJ0+PBhD02F2NhYLVu2TN26dVNBQYHmzZunX/3qVzp48KCCgoI8PV6jV1hYKEmXfd9cXAbPSExM1Pjx4xUTE6O8vDw9+eSTSkpKUk5Ojpo0aeLp8YxXXV2tadOm6eabb1bv3r0l/fx+8ff3v+T3znm/1J/LbRdJ+vWvf60OHTooMjJSBw4c0BNPPKHc3Fx99NFHbp2HsEOjk5SU5Px33759FRsbqw4dOuiDDz5QamqqBycDvNuECROc/+7Tp4/69u2rTp06aePGjRo5cqQHJ2sc0tPTdfDgQc4J9jJX2i4PPvig8999+vRRRESERo4cqby8PHXq1Mlt83Ao1o1atWqlJk2aXPLppKKiIoWHh3toKvyrkJAQde3aVUeOHPH0KJCc7w3eN96vY8eOatWqFe+depCRkaE1a9Zow4YNateunfP28PBwVVRUqKSkxGV93i/140rb5XJiY2Mlye3vF8LOjfz9/TVo0CBlZ2c7b6uurlZ2drbi4uI8OBn+WVlZmfLy8hQREeHpUSApJiZG4eHhLu8bh8Oh7du3877xMidPnlRxcTHvHTeyLEsZGRn6+OOPtX79esXExLgsHzRokPz8/FzeL7m5uTp+/DjvFze61na5nH379kmS298vHIp1s8zMTKWkpGjw4MEaOnSoFi1apPPnz2vy5MmeHq3RmjFjhsaOHasOHTooPz9fc+bMUZMmTXTfffd5erRGo6yszOX/tR49elT79u1TaGio2rdvr2nTpunZZ59Vly5dFBMTo1mzZikyMlLjxo3z3NCNwNW2S2hoqObNm6fk5GSFh4crLy9Pjz/+uDp37qyEhAQPTm229PR0LV++XKtWrVJQUJDzvDm73a5mzZrJbrcrNTVVmZmZCg0NVXBwsKZOnaq4uDjdeOONHp7eXNfaLnl5eVq+fLlGjx6tli1b6sCBA5o+fbqGDRumvn37une4ev8cbiP06quvWu3bt7f8/f2toUOHWtu2bfP0SI3avffea0VERFj+/v5W27ZtrXvvvdc6cuSIp8dqVDZs2GBJuuSSkpJiWdbPX3kya9YsKywszAoICLBGjhxp5ebmenboRuBq2+WHH36wRo0aZbVu3dry8/OzOnToYKWlpVmFhYWeHttol9sekqylS5c61/nxxx+thx9+2Lrhhhus5s2bW3fffbdVUFDguaEbgWttl+PHj1vDhg2zQkNDrYCAAKtz587WY489ZpWWlrp9Ntv/DQgAAIAGjnPsAAAADEHYAQAAGIKwAwAAMARhBwAAYAjCDgAAwBCEHQAAgCEIOwAAAEMQdgAAAIYg7ACgDhw7dkw2m835e5DuMmnSJNlsNtlsNq1cubLe5xk+fLjz+d39WgHUHGEHoF6dPn1aU6ZMUfv27RUQEKDw8HAlJCRo69atnh7tuk2aNOmS362NiopSQUGBevfu7fbnT0xMVEFBgZKSktz+XP/qo48+0o4dO+r9eQFcH19PDwCgcUlOTlZFRYXefvttdezYUUVFRcrOzlZxcbGnR1NlZaX8/Pxqdd8mTZooPDy8jie6vItB7AmhoaFyOBweeW4A18YeOwD1pqSkRFu2bNHzzz+v2267TR06dNDQoUOVlZWlO++802W9Bx54QK1bt1ZwcLBGjBih/fv3O5fPnTtX/fv31xtvvKGoqCg1b95c99xzj0pLS53r7Ny5U7fffrtatWolu92uW2+9VXv27HGZx2az6fXXX9edd96pFi1a6D//8z9VVVWl1NRUxcTEqFmzZurWrZt+//vfuzz322+/rVWrVjkPSW7cuPGyhz43bdqkoUOHKiAgQBEREZo5c6Z++ukn5/Lhw4frkUce0eOPP67Q0FCFh4dr7ty5tfrb7tixQwMGDFDTpk01ePBg7d2795J1Dh48qKSkJAUGBiosLEz333+/zpw541x+7tw5TZw4US1atFBERIQWLlyo4cOHa9q0abWaCUD9I+wA1JvAwEAFBgZq5cqVKi8vv+J6//7v/65Tp07p888/1+7duzVw4ECNHDlSZ8+eda5z5MgRffDBB1q9erXWrl2rvXv36uGHH3YuP3funFJSUvTVV19p27Zt6tKli0aPHq1z5865PNfcuXN1991365tvvtFvfvMbVVdXq127dlqxYoX++te/avbs2XryySf1wQcfSJJmzJihe+65x3k4tKCgQDfddNMlr+Hvf/+7Ro8erSFDhmj//v16/fXX9dZbb+nZZ591We/tt99WixYttH37di1YsEDPPPOM1q1bV6O/a1lZme644w717NlTu3fv1ty5czVjxgyXdUpKSjRixAgNGDBAu3bt0tq1a1VUVKR77rnHuU5mZqa2bt2qTz75ROvWrdOWLVsuiWEAXs4CgHr0l7/8xbrhhhuspk2bWjfddJOVlZVl7d+/37l8y5YtVnBwsHXhwgWX+3Xq1Ml64403LMuyrDlz5lhNmjSxTp486Vz++eefWz4+PlZBQcFln7eqqsoKCgqyVq9e7bxNkjVt2rRrzpyenm4lJyc7r6ekpFh33XWXyzpHjx61JFl79+61LMuynnzySatbt25WdXW1c53FixdbgYGBVlVVlWVZlnXrrbdat9xyi8vjDBkyxHriiSeuOMvlnvuNN96wWrZsaf3444/O215//XWXeX73u99Zo0aNcrnfiRMnLElWbm6u5XA4LD8/P2vFihXO5SUlJVbz5s2tRx999KqvFYD3YI8dgHqVnJys/Px8ffLJJ0pMTNTGjRs1cOBALVu2TJK0f/9+lZWVqWXLls49fIGBgTp69Kjy8vKcj9O+fXu1bdvWeT0uLk7V1dXKzc2VJBUVFSktLU1dunSR3W5XcHCwysrKdPz4cZd5Bg8efMmMixcv1qBBg9S6dWsFBgbqzTffvOR+13Lo0CHFxcXJZrM5b7v55ptVVlamkydPOm/r27evy/0iIiJ06tSpGj9X37591bRpU+dtcXFxLuvs379fGzZscPmbdu/eXZKUl5env/3tb6qsrNTQoUOd97Hb7erWrVuNZgHgWXx4AkC9a9q0qW6//XbdfvvtmjVrlh544AHNmTNHkyZNUllZmSIiIrRx48ZL7hcSEnLdz5GSkqLi4mL9/ve/V4cOHRQQEKC4uDhVVFS4rNeiRQuX6++//75mzJihl156SXFxcQoKCtILL7yg7du31+alXtO/fljDZrOpurq6zp+nrKxMY8eO1fPPP3/JsoiICB05cqTOnxNA/SPsAHhcz549nd/JNnDgQBUWFsrX11fR0dFXvM/x48eVn5+vyMhISdK2bdvk4+Pj3MO0detWvfbaaxo9erQk6cSJEy4fFLiSrVu36qabbnI5X++f9xRKkr+/v6qqqq76OD169NCHH34oy7Kce+22bt2qoKAgtWvX7ppz1ESPHj30zjvv6MKFC869dtu2bXNZZ+DAgfrwww8VHR0tX99L/9PfsWNH+fn5aefOnWrfvr0kqbS0VN99952GDRtWp/MCcB8OxQKoN8XFxRoxYoT++7//WwcOHNDRo0e1YsUKLViwQHfddZckKT4+XnFxcRo3bpz+53/+R8eOHdPXX3+tp556Srt27XI+VtOmTZWSkqL9+/dry5YteuSRR3TPPfc4vwakS5cueuedd3To0CFt375dEydOVLNmza45Y5cuXbRr1y598cUX+u677zRr1izt3LnTZZ3o6GgdOHBAubm5OnPmjCorKy95nIcfflgnTpzQ1KlTdfjwYa1atUpz5sxRZmamfHzq9j+9v/71r2Wz2ZSWlqa//vWv+uyzz/Tiiy+6rJOenq6zZ8/qvvvu086dO5WXl6cvvvhCkydPVlVVlYKCgpSSkqLHHntMGzZs0LfffqvU1FT5+Pi4HE4G4N0IOwD1JjAwULGxsVq4cKGGDRum3r17a9asWUpLS9Mf/vAHST8fivzss880bNgwTZ48WV27dtWECRP0v//7vwoLC3M+VufOnTV+/HiNHj1ao0aNUt++ffXaa685l7/11lv6xz/+oYEDB+r+++/XI488ojZt2lxzxoceekjjx4/Xvffeq9jYWBUXF7vsvZOktLQ0devWTYMHD1br1q0v++XKbdu21WeffaYdO3aoX79++u1vf6vU1FQ9/fTTtf3zXVFgYKBWr16tb775RgMGDNBTTz11ySHXyMhIbd26VVVVVRo1apT69OmjadOmKSQkxBmaL7/8suLi4nTHHXcoPj5eN998s3r06OFy7h4A72azLMvy9BAAUBNz587VypUrG+VPWk2aNEklJSVX/TmxunL+/Hm1bdtWL730klJTU523Hzt2TDExMdq7d6/69+/v9jkAXD/22AFAA7NmzRoFBgZqzZo1dfq4e/fu1Xvvvae8vDzt2bNHEydOlCTnYXJJSkpKUq9ever0eQHUHT48AQANyIIFC5yHcyMiIur88V988UXl5ubK399fgwYN0pYtW9SqVSvn8v/6r//Sjz/+KEnOD1kA8B4cigUAADAEh2IBAAAMQdgBAAAYgrADAAAwBGEHAABgCMIOAADAEIQdAACAIQg7AAAAQxB2AAAAhiDsAAAADPH/AA+MOhbeQxuBAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m = sep2d.value < 0.5\n",
    "plt.figure()\n",
    "plt.hist(sep2d, bins=\"auto\")\n",
    "plt.xlabel(\"Separation [deg]\")\n",
    "plt.ylabel(\"Counts\")\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"../../plots/sfi_gals_to_sfi_groups.png\", bbox_inches=\"tight\",\n",
    "            dpi=450)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2MTF "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [],
   "source": [
    "names = [\"RA\", \"DEC\", \"mag\", \"e_mag\", \"z_CMB\", \"r_hMpc\", \"e_rhMpc\", \"M\", \"eta\", \"e_eta\"]\n",
    "\n",
    "dtype = [(n, np.float32) for n in names]\n",
    "data = np.genfromtxt(join(tf_folder, \"twomtf_k.csv\"), delimiter=\",\", skip_header=1,\n",
    "                     usecols=[2 + n for n in range(len(names))])\n",
    "\n",
    "twomtf_gals = np.empty(len(data), dtype=dtype)\n",
    "for i, n in enumerate(names):\n",
    "    twomtf_gals[n] = data[:, i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Write to HDF5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {},
   "outputs": [],
   "source": [
    "outdir = \"/mnt/extraspace/rstiskalek/catalogs\"\n",
    "fname = \"PV_compilation.hdf5\"\n",
    "\n",
    "with File(join(outdir, fname), 'w') as f:\n",
    "    # Write LOSS\n",
    "    grp = f.create_group(\"LOSS\")\n",
    "    for name in loss_data.dtype.names:\n",
    "        grp.create_dataset(name, data=loss_data[name])\n",
    "\n",
    "    # Write Foundation\n",
    "    grp = f.create_group(\"Foundation\")\n",
    "    for name in foundation_data.dtype.names:\n",
    "        grp.create_dataset(name, data=foundation_data[name])\n",
    "\n",
    "    # Write SFI gals\n",
    "    grp = f.create_group(\"SFI_gals\")\n",
    "    for name in sfi_gals.dtype.names:\n",
    "        grp.create_dataset(name, data=sfi_gals[name])\n",
    "\n",
    "    # Write SFI gals masked\n",
    "    grp = f.create_group(\"SFI_gals_masked\")\n",
    "    for name in sfi_gals_masked.dtype.names:\n",
    "        grp.create_dataset(name, data=sfi_gals_masked[name])\n",
    "\n",
    "    # Write SFI groups\n",
    "    grp = f.create_group(\"SFI_groups\")\n",
    "    for name in sfi_groups.dtype.names:\n",
    "        grp.create_dataset(name, data=sfi_groups[name])\n",
    "\n",
    "    # Write 2MTF gals\n",
    "    grp = f.create_group(\"2MTF\")\n",
    "    for name in twomtf_gals.dtype.names:\n",
    "        grp.create_dataset(name, data=twomtf_gals[name])\n",
    "\n",
    "    # Write Pantheon\n",
    "    grp = f.create_group(\"Pantheon+\")\n",
    "    for name in pantheon_data.keys():\n",
    "        grp.create_dataset(name, data=pantheon_data[name])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv_csiborg",
   "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}