notebook for making DES y3 plots

des_y3_plots.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import fitsio\n",
    "import numpy as np\n",
    "from tqdm.notebook import tqdm\n",
    "import pickle\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "%matplotlib inline\n",
    "sns.set()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Slow Preprocessing of the Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def _cut_catalog(cat, *, post='', s2n, size):\n",
    "    msk = (\n",
    "        (cat['mcal_s2n_r' + post] > s2n) &\n",
    "        (cat['mcal_s2n_r' + post] < 100) &\n",
    "        ((cat['mcal_T_r' + post] / cat['mcal_Tpsf']) > 0.5) &\n",
    "        (cat['mcal_flags'] == 0)\n",
    "    )\n",
    "    return msk\n",
    "\n",
    "def _process_catalog(d):\n",
    "    size = 0.5\n",
    "    R11 = []\n",
    "    R22 = []\n",
    "    g1 = []\n",
    "    g2 = []\n",
    "    s2n_cut = []\n",
    "    tnames = []\n",
    "    \n",
    "    tilenames = np.unique(d['tilename'])\n",
    "    tname_msks = {}\n",
    "    for tilename in tqdm(tilenames, desc='tilenames init', leave=False):\n",
    "         tname_msks[tilename] = d['tilename'] == tilename\n",
    "    \n",
    "    for s2n in tqdm([14, 20, 25, 30, 40, 50], desc='s2n'):\n",
    "\n",
    "        gmsk = _cut_catalog(d, post='', s2n=s2n, size=size)\n",
    "        g1pmsk = _cut_catalog(d, post='_1p', s2n=s2n, size=size)\n",
    "        g1mmsk = _cut_catalog(d, post='_1m', s2n=s2n, size=size)\n",
    "        g2pmsk = _cut_catalog(d, post='_2p', s2n=s2n, size=size)\n",
    "        g2mmsk = _cut_catalog(d, post='_2m', s2n=s2n, size=size)\n",
    "\n",
    "        for tilename in tqdm(tilenames, desc='tilenames s2n', leave=False):\n",
    "            msk = tname_msks[tilename]\n",
    "            \n",
    "            tnames.append(tilename)\n",
    "            s2n_cut.append(s2n)\n",
    "\n",
    "            _msk = gmsk & msk\n",
    "            g1.append(np.mean(d['mcal_g'][_msk, 0]))\n",
    "            g2.append(np.mean(d['mcal_g'][_msk, 1]))\n",
    "\n",
    "            _msk = g1pmsk & msk\n",
    "            _g1_1p = np.mean(d['mcal_g_1p'][_msk, 0])\n",
    "\n",
    "            _msk = g1mmsk & msk\n",
    "            _g1_1m = np.mean(d['mcal_g_1m'][_msk, 0])\n",
    "\n",
    "            _msk = g2pmsk & msk\n",
    "            _g2_2p = np.mean(d['mcal_g_2p'][_msk, 1])\n",
    "\n",
    "            _msk = g2mmsk & msk\n",
    "            _g2_2m = np.mean(d['mcal_g_2m'][_msk, 1])\n",
    "\n",
    "            R11.append((_g1_1p - _g1_1m) / 2.0 / 0.01)\n",
    "            R22.append((_g2_2p - _g2_2m) / 2.0 / 0.01)\n",
    "\n",
    "        tqdm.write(\"%s: %s\" % (s2n, np.mean(g1) / np.mean(R11) / 0.02 - 1))\n",
    "    \n",
    "    return g1, g2, R11, R22, s2n_cut, tnames\n",
    "\n",
    "if False:\n",
    "    data = {}\n",
    "\n",
    "    for fname in tqdm(['g1p', 'g1m'], desc='files'):\n",
    "        tqdm.write('reading catalog: %s' % fname)\n",
    "        d = fitsio.read('e2e-007-grid/mcal_alltiles_%s.fits' % fname)\n",
    "        data[fname] = _process_catalog(d)\n",
    "\n",
    "    with open('e2e-007-grid/data_g1.pkl', 'bw') as fp:\n",
    "        pickle.dump(data, fp)\n",
    "        \n",
    "    data = {}\n",
    "\n",
    "    for fname in tqdm(['g2p', 'g2m'], desc='files'):\n",
    "        tqdm.write('reading catalog: %s' % fname)\n",
    "        d = fitsio.read('e2e-007-grid/mcal_alltiles_%s.fits' % fname)\n",
    "        data[fname] = _process_catalog(d)\n",
    "\n",
    "    with open('e2e-007-grid/data_g2.pkl', 'bw') as fp:\n",
    "        pickle.dump(data, fp)        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now compute `m` and `c`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"e2e-007-grid/data_g1.pkl\", \"rb\") as fp:\n",
    "    data1 = pickle.load(fp)\n",
    "    \n",
    "with open(\"e2e-007-grid/data_g2.pkl\", \"rb\") as fp:\n",
    "    data2 = pickle.load(fp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_m_and_c_pairs(data, plus, minus, swap=False):\n",
    "    tnames = set(data[plus][-1]) & set(data[minus][-1])\n",
    "    print(\"found %d tile pairs\" % len(tnames))\n",
    "    s2ns = set(data[plus][-2]) & set(data[minus][-2])\n",
    "    print(\"found %d s2n thresholds\" % len(s2ns))\n",
    "    \n",
    "    mpairs = {}\n",
    "    cpairs = {}\n",
    "    for s2n in s2ns:\n",
    "        mpairs[s2n] = []\n",
    "        cpairs[s2n] = []\n",
    "        \n",
    "        for tname in tnames:\n",
    "            pind = None\n",
    "            for _pind in range(len(data[plus][0])):\n",
    "                if data[plus][-2][_pind] == s2n and data[plus][-1][_pind] == tname:\n",
    "                    pind = _pind\n",
    "            \n",
    "            mind = None\n",
    "            for _mind in range(len(data[minus][0])):\n",
    "                if data[minus][-2][_mind] == s2n and data[minus][-1][_mind] == tname:\n",
    "                    mind = _mind\n",
    "\n",
    "            assert mind is not None\n",
    "            assert pind is not None\n",
    "        \n",
    "            if swap:\n",
    "                _m = 1\n",
    "                _c = 0\n",
    "            else:\n",
    "                _m = 0\n",
    "                _c = 1\n",
    "            mpairs[s2n].append((\n",
    "                (data[plus][_m][pind] - data[minus][_m][mind]) / 2,\n",
    "                (data[plus][_m+2][pind] + data[minus][_m+2][mind]) / 2,\n",
    "            ))\n",
    "    return mpairs, None\n",
    "\n",
    "def compute_m_and_c_bootstrap(mpairs, cpairs, n_boot=500, seed=10):\n",
    "    rng = np.random.RandomState(seed=seed)\n",
    "\n",
    "    mvals = {}\n",
    "    for s2n in tqdm(mpairs, desc='s2n'):\n",
    "        mvals[s2n] = []\n",
    "        for _ in range(n_boot):\n",
    "            inds = rng.choice(len(mpairs[s2n]), size=len(mpairs[s2n]), replace=True).astype(np.int)\n",
    "            mvals[s2n].append(\n",
    "                np.mean([mpairs[s2n][i][0] for i in inds]) \n",
    "                / np.mean([mpairs[s2n][i][1] for i in inds]) / 0.02 - 1\n",
    "            )\n",
    "        mvals[s2n] = (np.mean(mvals[s2n]), np.std(mvals[s2n]))\n",
    "        \n",
    "    return mvals, None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "found 196 tile pairs\n",
      "found 6 s2n thresholds\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "be56c96de32145b5bc32154c8e63d1ff",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, description='s2n', max=6.0, style=ProgressStyle(description_width='ini…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "found 190 tile pairs\n",
      "found 6 s2n thresholds\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f8b92300e5ae4f60b979958f0da49a1b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, description='s2n', max=6.0, style=ProgressStyle(description_width='ini…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "mpairs1, cpairs1 = build_m_and_c_pairs(data1, 'g1p', 'g1m')\n",
    "mvals1, cvals1 = compute_m_and_c_bootstrap(mpairs1, cpairs1)\n",
    "\n",
    "mpairs2, cpairs2 = build_m_and_c_pairs(data2, 'g2p', 'g2m', swap=True)\n",
    "mvals2, cvals2 = compute_m_and_c_bootstrap(mpairs2, cpairs2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.00952179 -0.00530148 -0.00484235 -0.00383195 -0.00747722 -0.00759824]\n",
      "[0.0030038  0.00255513 0.00259381 0.00290784 0.00286139 0.00304335]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEJCAYAAAC3yAEAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3dfXxMd97/8VcymUQ0cdtJEGrVVatXUe0qpdukFEnkBukdsrGq9Kpuq2Wlbnotl9pWqa0ti95c9qbWTelNwvVIgqVqS9qKrpvuj1xde+mWkEwkKgiZZM7vD82sMYnNcDIzeD8fDw853+/3zHzON2fyzpmTOSfIMAwDERERkwT7uwAREbm+KFhERMRUChYRETGVgkVEREylYBEREVMpWERExFQKFhERMVWIvwsIFOXlZ3A6vf9IT+vWEZw4cboRKro6qss7qss7qss712NdwcFBtGx5U519CpbvOZ3GFQVL7bqBSHV5R3V5R3V550aqS2+FiYiIqRQsIiJiKr8Ey4YNGxgyZAiDBw9m5cqVHv0HDhwgLS2N+Ph4XnzxRaqrqwEoKioiPT2dhIQEJkyYwJkzZwD44osv6NOnD0OHDmXo0KFMnz7dp9sjIiL/5PNgKS4uZuHChaxatYqsrCzee+89/va3v7mNyczMZObMmWzcuBHDMFi7di0As2fPZtSoUeTl5dGtWzeWLl0KwFdffcXYsWPJzs4mOzubuXPn+nqzRETkez4Plp07d3LvvffSokULmjZtSnx8PHl5ea7+o0ePcu7cOXr27AlAWloaeXl5OBwOdu3aRXx8vFs7wP79+/n0009JSUnhqaee4tixY77eLBER+Z7Pg6WkpASbzeZajoqKori4uN5+m81GcXEx5eXlREREEBIS4tYOEBkZSUZGBhs2bCAuLo5Jkyb5aGtERORSPv9zY6fTSVBQkGvZMAy35fr6Lx0HuJZfeuklV9vIkSP51a9+RUVFBZGRkQ2uq3XrCK+3pZbN1vDn8SXV5R3V5R3V5Z0bqS6fB0ubNm0oKChwLdvtdqKiotz67Xa7a7m0tJSoqChatWpFRUUFNTU1WCwW13pOp5O33nqLJ598EovF4lrv4q8b4sSJ01f099w2WyR2e4XX6zU21eWdQKxr3sovsYZamPzInf4uxUMgzheoLm9dTV3BwUH1/kLu87fC+vXrR35+PmVlZVRWVrJp0yZiY2Nd/TExMYSFhbF7924AsrOziY2NxWq10qtXL3JycgDIysoiNjaW4OBgNm/ezMaNG13td955J02bNvX1pomICH4IlujoaCZNmsTo0aMZNmwYycnJ9OjRg/Hjx7N//34AFixYwNy5c0lISODs2bOMHj0agFmzZrF27VqGDBlCQUEBzz//PADz5s3j3XffJSkpiQ8++IBf/vKXvt4sERH5XpDueX+B3grzDdXVcHorzHuqyzvXzVthIiJyfVOwiIiIqRQsIiJiKgWLiIiYSsEiIiKmUrCIiIipFCwiImIqBYuIiJhKwSIiIqZSsIiIiKkULCIiYioFi4iImMrn92MRCTSBfLFHkWuRjlhERMRUChYRETGVgkVEREylYBEREVPp5L34jE6Si9wYdMQiIiKmUrCIiIipFCwiImIqBYuIiJhKwSIiIqZSsIiIiKn058YiIjegxvzzfx2xiIiIqRQsIiJiKgWLiIiYSudYROS6oEsGBQ4dsYiIiKkULCIiYioFi4iImErBIiIipvJLsGzYsIEhQ4YwePBgVq5c6dF/4MAB0tLSiI+P58UXX6S6uhqAoqIi0tPTSUhIYMKECZw5cwaAU6dO8eSTT5KYmEh6ejp2u92n2yMiIv/k82ApLi5m4cKFrFq1iqysLN577z3+9re/uY3JzMxk5syZbNy4EcMwWLt2LQCzZ89m1KhR5OXl0a1bN5YuXQrAr3/9a3r16kVubi6PPPIIL7/8sq83S0REvufzYNm5cyf33nsvLVq0oGnTpsTHx5OXl+fqP3r0KOfOnaNnz54ApKWlkZeXh8PhYNeuXcTHx7u1A2zbto2UlBQAkpOT2b59Ow6Hw8dbJiIi4IfPsZSUlGCz2VzLUVFR7Nu3r95+m81GcXEx5eXlREREEBIS4tZ+6TohISFERERQVlZGdHR0g+uaOnUyJSUlruXBgxN57LFRVFZW8swzT3qMT00dztChaZSVlfHEE2M9+h99dCTx8UM4fvwYL774gkf/6NGPExc3gMOH/86cObM8+sePn8C99/bj4MEDvPbaKx79zz47iZ4972bPni9ZvHihR/8rr/yS6OiOfPbZTt55Z5lH/y9+MZsf/OBWPvlkK++++zuP/pdfnk+bNm3ZuDGHtWtXe/QvWLCIli1bkp39IevXf+TR/5vfvE14eDjvvbeKTZtyAfi25DRBQUH8Ne8mli9fAcAf/rCc7du3ua3bpEkTlix5B4C3317K55/nu/W3aNGCX/1qMQCLFv2KvXv3uPVHR7fhlVdeA2D+/FcoLDzg1t+x4w+YOXMOAC+99As+3fWVqy6AH/7wdl54YQYAM2ZkUlx83G39O+/sycSJPwfg5z9/lpMnT7r19+nTlyeffBqAn/1sPOfOnXPrj419gJ/+9AkAnngiw2Puave9asd5tr73mquuWrX7Xnl5OVOmTPRY3xf73qBBcfXue5mZM+ja9Xaf73u1+9eE5JUe+97F/LHvhYaGUFVV7bHvffPNYbf1fbnvbVkzx22/h4b/3CsvL2f+/F/y3//9tscY8EOwOJ1OgoKCXMuGYbgt19d/6TjAY/nidYKDvTsYs1othIb+czoiIsKw2SKprAxxa68VGdkEmy2SsrKyy/ZXVZ2qs79Zs3Bstki+++6mOvubN7/QX1zctM7+Fi2aYrNF0qJF3f0ANlskzZuH19nfsuVN2GyRNGtWd3+rVhf6IyOb1Nl/880RtGpVf7/NFkl4eDgREWGu/trvV2hoCDZbJAAREZ7rX9zftGmoR39YmPWi/jCP/iZN/tkfHm716A8PD72oP9Strtp1avubNPFcv2nTMFd/WFhd/f98/NDQEJxO9/6IiCZu/Zeq3fesoZY6x9TuWxaLw2/7Xu3/de9bTf2y79V+H+va9y7mr30vNDTEY9/z3Dd9t+9dut9f6G/Yzz2LxYHVavHorxVkGIZRb28j+OijjygoKHCdB1myZAmGYfDMM88AF94KGzNmDJs3bwagoKCARYsWsXz5cvr06cOuXbuwWCwcO3aMn/zkJ2zZsoUBAwawatUq2rRpQ3V1Nb179+bzzz/HarU2uK4TJ07jdHo/FTZbJHZ7hdfrNbZArCtQPxmturyn/cs71+N8BQcH0bp1RN19V1PYlejXrx/5+fmUlZVRWVnJpk2biI2NdfXHxMQQFhbG7t27AcjOziY2Nhar1UqvXr3IyckBICsry7VeXFwcWVlZAOTk5NCrVy+vQkVERMzj82CJjo5m0qRJjB49mmHDhpGcnEyPHj0YP348+/fvB2DBggXMnTuXhIQEzp49y+jRowGYNWsWa9euZciQIRQUFPD8888D8Nxzz7Fnzx6SkpJYtWoVM2fO9PVmiYjI9/xyEcqUlBTXX3HVeuedd1xfd+3alffff99jvZiYGFasWOHR3qJFC958803zCw1Q81Z+CcDU9Lv9XImIiCd98l5EREylYBEREVMpWERExFQKFhERMZWCRURETKVgERERUylYRETEVAoWERExlYJFRERMpWARERFTKVhERMRUChYRETGVgkVEREylYBEREVMpWERExFQKFhERMZWCRURETKVgERERUylYRETEVAoWERExlYJFRERMpWARERFTKVhERMRUChYRETGVgkUkAOX/9TiHir7jq0MnyFy6g/y/Hvd3SSINpmARCTD5fz3OH3IPUl1jAHDi1Hn+kHtQ4SLXDAWLSID58JNDVFU73dqqqp18+MkhP1Uk4h0Fi0iAOXHqvFftIoFGwSISYFo3C/OqXSTQKFhEAkxaXGdCQ9xfmqEhwaTFdfZTRSLeCfF3ASLiru8dbQD4Xc4BqmsMWjcLIy2us6tdJNApWEQCUN872rB9TxHWUAuTH7nT3+WIeEVvhYmIiKkULHJD0wcRRczn87fCioqKyMzM5MSJE3Tq1IkFCxZw0003uY2pqqrixRdf5KuvvqJJkyYsWLCAzp07YxgG8+fP5+OPPyY4OJg5c+bwox/9CIAHH3yQiIgI12O8+eabtG3b1qfbJteW+j6ICOh8hshV8PkRy+zZsxk1ahR5eXl069aNpUuXeoxZsWIF4eHh5ObmMmPGDKZPnw7Axo0bOXToEDk5OSxZsoTp06dTXV1NeXk5VquV7Oxs1z+Fivwr+iCiSOPwabA4HA527dpFfHw8AGlpaeTl5XmM27ZtG6mpqQDcc889lJWVUVRUxCeffMKQIUMIDg6mU6dOtG3blr/85S/s378fwzAYMWIEw4cPJzc315ebJdcofRBRpHH49K2w8vJyIiIiCAm58LQ2m43i4mKPcSUlJdhsNteyzWbj+PHjlJSUEBUV5dEeHh7O/fffz5QpUygtLSU9PZ0uXbrQuXPD/+6/deuIfz2oHjZb5BWveyWsoZYGPa+v67qcbbu/5e/HTuGodjL1rXxGJ97OAz/q4NeabC3DsZdX1tkeCHPX0O+zvwRaXZov7zTmfDVasOTm5jJ37ly3to4dOxIUFOTWdukygGEYbu2GYRAcHIzT6ayzfeDAgQwcOBCA9u3bM2jQID799FOvguXEidM4nUaDx9ey2SKx2yu8Xu9qOKpqAC77vP6oqz615zIc37/tZC+vZPHaPZyqOOfXcxnDftyJP+QedHs7LDQkmGE/7hQQc+eoqsEaagmIWi4VSPtXLc2Xd652voKDg+r9hbzRgiUxMZHExES3NofDQZ8+faipqcFisWC3292OQGpFR0dTUlLCLbfcAkBpaSlRUVG0adOGkpIS17ja9o8//pibb76Z7t27u/pqj4rE/y53LsOfwaIPIoo0Dp+eY7FarfTq1YucnBwAsrKyiI2N9RgXFxdHdnY2AAUFBYSFhdGuXTtiY2PZsGEDNTU1fPPNNxw+fJju3btz9OhRlixZgtPppLS0lK1bt/LAAw/4ctPkMgL5XEbfO9rQuV1zunVuzWtP36dQETGBz3+tnzVrFtOmTWPZsmW0bduW119/HYDVq1dTUlLCc889R0ZGBjNnziQpKYnQ0FDmz58PQEJCAvv27XOd2H/55Zdp0qQJI0aMoLCwkOTkZJxOJ1OmTCEmJsbXmyb1aN0srM4Q0UUVRa5PPg+WmJgYVqxY4dE+cuRI19dhYWHMmzfPY0xQUBBTp05l6tSpbu0hISHMmTPH/GLFFGlxnes8l6GLKopcn3QiQhqdzmWI3FgULOITuqiiyI1D1woTERFTKVhERMRUChYRETGVgkVEREzV4JP3drud7777zq3t3/7t30wvSERErm0NCpa5c+eycuVKt/udBAUFkZ+f32iFiYjItalBwbJ582b+/Oc/07Jly8auR0RErnENOsfygx/8gGbNmjV2LSIich1o0BFLRkYGP/nJT+jTp4/bVYOfeeaZRitMRESuTQ0KlrfffpuIiAgqKgLrfgIiIhJ4GhQslZWVrF69urFrERG5Ivl/Pc6hou+orjHIXLpD16LzswadY+nUqRMHDx5s7FpERLxWe4fS6poLd4A9ceo8f8g9SP5fj/u5ssBVG8RfHTpB5tIdps9Vg45Yjh07xsMPP0xMTAyhoaGu9g0bNphajIiItwL1DqWBqr4gBkybrwYFy+TJk015MhERswXyHUoDkS+CuEHB0rt3b1OeTETEbLpDqXd8EcS6VpiIXNPS4joTGuL+o0x3KK1ffYFrZhArWETkmtb3jjb8NLErIZYg4MIPyJ8mdtX5lXr4Ioh1B0kRuebpDqUN54tbhStYRERuMI0dxHorTERETKVgERERUylYRETEVAoWERExlYJFRERMpWARERFTKVhERMRUChYRETGVguUaU3sfhcJvTzbKfRRERK6WguUaohsaici1QMFyDbncfRRERAKFguUaohsaici1wOfBUlRURHp6OgkJCUyYMIEzZ854jKmqqiIzM5PExESGDx/OoUPuv5EXFhaSlJTk1vbb3/6WhIQE4uPj2bRpU6Nug7/44j4KIiJXy+fBMnv2bEaNGkVeXh7dunVj6dKlHmNWrFhBeHg4ubm5zJgxg+nTp7v6srKyGDduHJWVla62ffv2sX79erKzs1m1ahXz58/n5MmTPtkeX9INjUTkWuDTYHE4HOzatYv4+HgA0tLSyMvL8xi3bds2UlNTAbjnnnsoKyujqKiIiooKtmzZwuuvv+42fvv27QwaNIiwsDBat25N79692bZtW6Nvj6/phkYici3w6f1YysvLiYiIICTkwtPabDaKi4s9xpWUlGCz2VzLNpuN48eP065dOxYvXsyRI0c8xnfv3t1jvDdat47wavzFbLbIK17XW6kPRJL//y7M2dynf3zZsb6sqyGsoRZAdTVUoNZVK9Dq0nx5pzHnq9GCJTc3l7lz57q1dezYkaCgILe2S5cBDMNwazcMg+Dg+g+unE6nR9vlxtflxInTOJ2GV+vAhW+K3V7h9XpXw1FVA3DZ5/VHXf+Ko6oGa6hFdTVQoNYF2r+8dT3OV3BwUL2/kDdasCQmJpKYmOjW5nA46NOnDzU1NVgsFux2O1FRUR7rRkdHU1JSwi233AJAaWlpneNqtWnTBrvd7lq22+106tTJpC0RERFv+PQci9VqpVevXuTk5AAXTsTHxsZ6jIuLiyM7OxuAgoICwsLCaNeuXb2PGxsby6ZNm6isrKSsrIzPPvuMvn37Ns5GiIjIZfn8nvezZs1i2rRpLFu2jLZt27pOxK9evZqSkhKee+45MjIymDlzJklJSYSGhjJ//vzLPmaPHj1ITU3l4Ycfprq6mokTJxIdHe2LzRERkUv4PFhiYmJYsWKFR/vIkSNdX4eFhTFv3rx6H6N9+/Zs3brVrW3s2LGMHTvWvEJFROSK6JP3IiJiKgWLiIiYSsEiIiKmUrCIiIipFCwiImIqBYuIiJhKwSIiIqZSsIiIiKkULCIiYioFi4iImErBIiIiplKwiIiIqRQsIiJiKgWLiIiYSsEiIiKmUrCIiIipfH6jLxFpmKnpd2OzRWK3V/i7FBGvKFjkhqcf4CLm0lthIiJiKh2xiM/oyEDkxqAjFhERMZWCRURETKVgERERUylYRETEVAoWERExlYJFRERMpWARERFTKVhERMRUChYRETGVgkVEREylYBEREVMpWERExFQKFhERMZXPg6WoqIj09HQSEhKYMGECZ86c8RhTVVVFZmYmiYmJDB8+nEOHDrn1FxYWkpSU5Nb24IMPMnToUNe/Y8eONep2iIhI3XweLLNnz2bUqFHk5eXRrVs3li5d6jFmxYoVhIeHk5uby4wZM5g+fbqrLysri3HjxlFZWelqKy8vx2q1kp2d7frXtm1bn2yPiIi482mwOBwOdu3aRXx8PABpaWnk5eV5jNu2bRupqakA3HPPPZSVlVFUVERFRQVbtmzh9ddfdxu/f/9+DMNgxIgRDB8+nNzc3MbfGBERqZNPb/RVXl5OREQEISEXntZms1FcXOwxrqSkBJvN5lq22WwcP36cdu3asXjxYo4cOeI2vqqqivvvv58pU6ZQWlpKeno6Xbp0oXPnzo27QSIi4qHRgiU3N5e5c+e6tXXs2JGgoCC3tkuXAQzDcGs3DIPg4PoPrgYOHMjAgQMBaN++PYMGDeLTTz/1Klhat45o8NhL2WyRV7zulbCGWhr0vL6uq6FUl3dUV8M09HXhL4FWV2POV6MFS2JiIomJiW5tDoeDPn36UFNTg8ViwW63ExUV5bFudHQ0JSUl3HLLLQCUlpbWOa7Wxx9/zM0330z37t1dbbVHRQ114sRpnE7Dq3UAv9xq11FVA3DZ5w3UWwCrLu+oroZzVNVgDbUEXF1wfc5XcHBQvb+Q+/Qci9VqpVevXuTk5AAXTsTHxsZ6jIuLiyM7OxuAgoICwsLCaNeuXb2Pe/ToUZYsWYLT6aS0tJStW7fywAMPNMo2iIjI5fn0HAvArFmzmDZtGsuWLaNt27auE/GrV6+mpKSE5557joyMDGbOnElSUhKhoaHMnz//so85YsQICgsLSU5Oxul0MmXKFGJiYnyxOSIicgmfB0tMTAwrVqzwaB85cqTr67CwMObNm1fvY7Rv356tW7e6lkNCQpgzZ465hYqIyBXRJ+9FRMRUChYRETGVgkVEREylYBEREVP5/OS9iEhjmJp+d0B+XuRGpCMWERExlYJFRERMpWARERFTKVhERMRUOnkvInIDasw/dtARi4iImErBIiIiplKwiIiIqRQsIiJiKgWLiIiYSsEiIiKmUrCIiIipFCwiImIqBYuIiJhKwSIiIqZSsIiIiKkULCIiYioFi4iImErBIiIiptJl869BU9Pv9ncJIiL10hGLiIiYSsEiIiKmUrCIiIipFCwiImIqBYuIiJhKwSIiIqZSsIiIiKkULCIiYip9QPJ7wcFBflm3Maku76gu76gu71xvdV1uvSDDMIwrLUhERORSeitMRERMpWARERFTKVhERMRUChYRETGVgkVEREylYBEREVMpWERExFQKFhERMZWCRURETKVg8dLp06dJTk7myJEjAOzcuZOUlBQGDx7MwoULA6au6dOnM3jwYIYOHcrQoUPZvHmzz2v6zW9+Q1JSEklJScyfPx8IjPmqq65AmK833niDIUOGkJSUxO9+9zsgMOarrroCYb5qzZs3j2nTpgGBMV911RUo85WRkUFSUpKrjr179zbOnBnSYHv27DGSk5ONO+64w/j222+NyspKIy4uzvjHP/5hOBwOY+zYsca2bdv8XpdhGEZycrJRXFzs81pq7dixw3jssceM8+fPG1VVVcbo0aONDRs2+H2+6qpr06ZNfp+vzz//3BgxYoThcDiMyspKo3///saBAwf8Pl911XXo0CG/z1etnTt3Gn369DGmTp0aMK/HS+syDP+/Hg3DMJxOp/HjH//YcDgcrrbGmjMdsXhh7dq1zJo1i6ioKAD27dtHx44d6dChAyEhIaSkpJCXl+f3uiorKykqKmLGjBmkpKSwaNEinE6nT2uy2WxMmzaN0NBQrFYrnTt35vDhw36fr7rqKioq8vt89e7dm3fffZeQkBBOnDhBTU0Np06d8vt81VVXkyZN/D5fACdPnmThwoU89dRTQOC8Hi+tKxBejwB///vfARg7diypqan88Y9/bLQ5U7B44eWXX6ZXr16u5ZKSEmw2m2s5KiqK4uJiv9dVWlrKvffeyyuvvMLatWspKCjg/fff92lNt912Gz179gTg8OHD5ObmEhQU5Pf5qquu+++/3+/zBWC1Wlm0aBFJSUn07ds3YPavS+uqrq4OiPmaOXMmkyZNolmzZkDgvB4vrSsQXo8Ap06dom/fvixZsoTf//73rFmzhqKiokaZMwXLVXA6nQQF/fPS0YZhuC37S4cOHViyZAlRUVGEh4eTkZHBJ5984pdavv76a8aOHcsLL7xAhw4dAma+Lq7r1ltvDZj5mjhxIvn5+Rw7dozDhw8HzHxdXFd+fr7f52vdunW0bduWvn37utoC4fVYV12B8nq86667mD9/PpGRkbRq1YqHH36YRYsWNcqc6X4sV6FNmzbY7XbXst1ud70d5U+FhYUcPnyY+Ph44MLOEhLi+2/17t27mThxIjNmzCApKYkvvvgiIObr0roCYb4OHTpEVVUVt99+O+Hh4QwePJi8vDwsFotrjD/mq666cnJyaNGihV/nKycnB7vdztChQ/nuu+84e/YsR48e9ft81VXXz372M1JTU/3+eiwoKMDhcLhCzzAMYmJiGuc1edVnaW5A/fv3N7799lvj3LlzRmxsrHH48GGjurraeOKJJ4ycnBy/13XgwAEjNjbWOHnypFFVVWWMHTvW2LBhg09rKSoqMvr06WPs3LnT1RYI81VXXYEwX9u2bTPS0tKM8+fPG+fPnzcef/xx43/+53/8Pl911bV48WK/z9fFPvjgA2Pq1KkBsX/VVVcg7F+GYRhbt241hg0bZpw7d86oqKgwUlJSjC+//LJR5kxHLFchLCyMV199lWeffZbz588TFxdHQkKCv8uia9euPPnkk4wcOZLq6moGDx5McnKyT2tYvnw558+f59VXX3W1jRgxwu/zVV9d/p6vuLg49u3bx7Bhw7BYLAwePJikpCRatWrl1/mqq65nnnmGli1b+nW+6qLX4+X179+fvXv3MmzYMJxOJ6NGjeKuu+5qlDnTHSRFRMRUOnkvIiKmUrCIiIipFCwiImIqBYuIiJhKwSIiIqZSsMh1bfXq1aSmprquzpuZmUlRUZGrPyMjo85rIx05coTbb7/ddRXYlJQURowYQU5OTr3PtWfPHjIyMkhJSSE5OZlx48bx9ddfu42ZOHEiX3/9NRkZGWRkZLhdM6qsrIwf/vCHdT72unXrWLlyJQCLFy/mpZde8moeGurIkSPcddddXq83YMAA9u/f79Gel5dHRkaGGaXJNUSfY5Hr1rx58zh48CBvvfUWbdu2xel0sn79eh577DHWrVtHmzZtLrt+kyZNyM7Odi0fPXqUMWPGYLFYXJ+irlVVVcV//Md/8Nvf/pY77rgDgOzsbMaPH8+WLVuwWCxUVVXxj3/8g9tuuw24EERvvvkmTz/99L/clt27d7vWEwl0Cha5Lh0/fpw1a9awbds2mjdvDkBwcDDDhg3jq6++4q233mLWrFlePWZMTAwTJ05k+fLlHsFSWVlJRUUFZ8+edbWlpqYSERFBTU0NFouFnTt3ul1D6umnn2b58uX069fPdWHMumzevJmtW7eyY8cOmjRpAly4Um1GRgZ2u52bb76Z119/naioKAYMGECPHj0oLCxk8uTJ9OjRg5deeoljx47hcDhISkriqaeeorq6mjlz5vDll19itVpp3749c+fOBaCmpoaZM2eyf/9+KioqyMzMJD4+HofDwauvvkp+fj4Wi4UePXowffp0IiIi3Op944032LBhAy1atKBjx45ezbFcH/RWmFyX9u7dy6+SmXYAAAQ8SURBVK233uoKlYv169eP3bt3X9Hjdu3alf/93//1aG/evDmZmZmMGzeOBx98kMzMTD744AP69etHaGgoAH/605948MEHXet06tSJF154gSlTpnD69Ol6n3PQoEEMGDCAMWPGkJ6eDsC3337LG2+8QV5eHs2aNWPdunWu8bfddhu5ubkMGjSIzMxMHnroIT788EPef/99du7cSU5ODnv27OGLL75g/fr1fPjhh3To0IHCwkIAzp8/z3333cdHH33E1KlTee211wBYtmwZJSUlZGdnk52djdPpdN0ordaf/vQnNm3aRFZWFmvWrLnsdsn1S8Ei163q6uo626uqqq74Cq5BQUGuo4ZLPf744+zYsYP//M//xGaz8c477zBs2DAqKiowDIO9e/dy9913u63z6KOP8u///u/813/9l1d13HfffbRq1Qq4EHZlZWWuvtpbKJw9e5Zdu3bxxhtvMHToUB599FGOHTvGwYMH6dKlCxaLhUceeYRf//rXxMfHu2qzWq2uI7KuXbty4sQJALZv386IESOwWq0EBweTkZHBn//8Z7e68vPzGTRoEBEREYSEhPDQQw95tV1yfdBbYXJd6tmzJ9988w12u93tfhMAn3/++RWdoAbYv38/Xbp08WjfvXs3f/nLXxg3bhz9+/enf//+TJ48meTkZHbs2EFUVBTdunUjONjzd7k5c+aQmprK+vXrG1zHxVfHDQoK4uIrMzVt2hS4cBl5wzBYs2YN4eHhwIU/EAgLC+Omm24iOzubL7/8ks8++4znn3+eJ554gri4OKxWq9tj17r0svROpxOHw+FR28W1XHy1Yblx6IhFrkvR0dFkZGQwefJktxsXffDBB2zatInx48d7/Zj/93//x9KlSxk7dqxHX6tWrVi2bBkFBQWuNrvdzunTp+nSpQtbtmxh4MCBdT5u8+bNee211y57v3GLxVLvEVh9IiIi6Nmzp+s+9adOnWLkyJFs2bKFjz/+mDFjxnDXXXfx7LPPus49Xc7999/P6tWrcTgcOJ1OVq5cyX333ec2JjY2lry8PE6dOoXT6XT74we5ceiIRa5bP//5z1m3bh0TJkygqqqKqqoqunfvzpo1a4iJiXGNe+GFF5g+fbpredSoUYwcOZJz584xdOhQ4MKJ/7CwMCZPnswDDzzg8VydOnViyZIlLFy4kOPHjxMWFkZkZCSvvPIKt956Kzt37uTZZ5+tt9bevXszZswY3nzzzTr7Y2Nj3a7I3FALFixgzpw5pKSkUFVVRXJyMqmpqdTU1LB9+3aSk5Np2rQpzZs3Z86cOZd9rAkTJjBv3jyGDRtGdXU1PXr04Be/+IXbmLi4OAoLC3nooYdo1qwZXbt2pby83Ou65dqmqxuLiIip9FaYiIiYSsEiIiKmUrCIiIipFCwiImIqBYuIiJhKwSIiIqZSsIiIiKkULCIiYqr/D1h2kTFXVh8pAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = plt.gca()\n",
    "\n",
    "m1 = []\n",
    "m1err = []\n",
    "\n",
    "m2 = []\n",
    "m2err = []\n",
    "\n",
    "m = []\n",
    "merr = []\n",
    "\n",
    "s2ns = []\n",
    "for s2n in sorted(list(mvals1.keys())):\n",
    "    s2ns.append(s2n)\n",
    "\n",
    "    m1.append(mvals1[s2n][0])\n",
    "    m1err.append(mvals1[s2n][1])\n",
    "    m2.append(mvals2[s2n][0])\n",
    "    m2err.append(mvals2[s2n][1])\n",
    "    \n",
    "    m.append((mvals1[s2n][0] + mvals2[s2n][0]) / 2)\n",
    "    merr.append(np.sqrt(mvals1[s2n][1]**2 + mvals2[s2n][1]**2) / 2)\n",
    "\n",
    "m = np.array(m)\n",
    "merr = np.array(merr)\n",
    "\n",
    "print(m) \n",
    "print(merr)\n",
    "\n",
    "ax.errorbar(s2ns, m, yerr=merr*3, fmt='o')\n",
    "# ax.set_ylim(-0.03, 0.005)\n",
    "ax.hlines(0, 8, 52, linestyle='dashed', color='k')\n",
    "ax.set_xlim(8, 52)\n",
    "ax.set_xlabel(\"OLD S/N threshold\")\n",
    "ax.set_ylabel(\"m\")\n",
    "plt.savefig(\"m-e2e-007-grid.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}