GWSignalProcessing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PyCBC Tutorial: Generating Waveforms and Matched Filtering\n",
    "### Alex Nitz\n",
    "\n",
    "We will be using the [PyCBC](http://github.com/ligo-cbc/pycbc) library, which is used to study gravitational-wave data, find astrophysical sources due to compact binary mergers, and study their parameters. These are some of the same tools that the LIGO and Virgo collaborations use to find gravitational waves in LIGO/Virgo data \n",
    "\n",
    "In this tutorial we will walk through how to find a specific signal in LIGO data. We present how to generate the waveform of a gravitational-wave merger and matched filtering, which is optimal in the case of Gaussian noise and a known signal model. In reality our noise is not entirely Guassian, and in practice we use a variety of techniques to separate signals from noise in addition to the use of the matched filter. \n",
    "\n",
    "Additional [examples](http://pycbc.org/pycbc/latest/html/#library-examples-and-interactive-tutorials) and module level documentation are [here](http://pycbc.org/pycbc/latest/html/py-modindex.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Getting the software environment setup\n",
    "\n",
    "PyCBC is installable through pip, but also relies on portions of the [LALSuite](https://git.ligo.org/lscsoft/lalsuite) c-library. A bundled version of this suitable for use with PyCBC is also available on Mac / Linux through pip. These can be installed as follows within the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: pycbc in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (1.13.8)\n",
      "Requirement already satisfied: lalsuite in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (6.54)\n",
      "Requirement already satisfied: ligo-common in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (1.0.3)\n",
      "Requirement already satisfied: numpy<1.15.3,>=1.13.0 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pycbc) (1.15.2)\n",
      "Requirement already satisfied: Mako>=1.0.1 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pycbc) (1.0.7)\n",
      "Requirement already satisfied: cython in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pycbc) (0.29)\n",
      "Requirement already satisfied: decorator>=3.4.2 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pycbc) (4.3.0)\n",
      "Requirement already satisfied: scipy>=0.16.0 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pycbc) (1.1.0)\n",
      "Requirement already satisfied: matplotlib>=1.5.1 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pycbc) (2.2.3)\n",
      "Requirement already satisfied: pillow in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pycbc) (5.3.0)\n",
      "Requirement already satisfied: h5py>=2.5 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pycbc) (2.8.0)\n",
      "Requirement already satisfied: jinja2 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pycbc) (2.10)\n",
      "Requirement already satisfied: mpld3>=0.3 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pycbc) (0.3)\n",
      "Requirement already satisfied: lscsoft-glue>=1.59.3 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pycbc) (2.0.0)\n",
      "Requirement already satisfied: emcee==2.2.1 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pycbc) (2.2.1)\n",
      "Requirement already satisfied: requests>=1.2.1 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pycbc) (2.19.1)\n",
      "Requirement already satisfied: beautifulsoup4>=4.6.0 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pycbc) (4.6.3)\n",
      "Requirement already satisfied: six>=1.10.0 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pycbc) (1.11.0)\n",
      "Requirement already satisfied: ligo-segments in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pycbc) (1.2.0)\n",
      "Requirement already satisfied: astropy<3.0.0,>=2.0.3 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pycbc) (2.0.8)\n",
      "Requirement already satisfied: weave>=0.16.0 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pycbc) (0.17.0)\n",
      "Requirement already satisfied: python-dateutil in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from lalsuite) (2.7.3)\n",
      "Requirement already satisfied: MarkupSafe>=0.9.2 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from Mako>=1.0.1->pycbc) (1.0)\n",
      "Requirement already satisfied: cycler>=0.10 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from matplotlib>=1.5.1->pycbc) (0.10.0)\n",
      "Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from matplotlib>=1.5.1->pycbc) (2.2.2)\n",
      "Requirement already satisfied: pytz in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from matplotlib>=1.5.1->pycbc) (2018.5)\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from matplotlib>=1.5.1->pycbc) (1.0.1)\n",
      "Requirement already satisfied: backports.functools_lru_cache in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from matplotlib>=1.5.1->pycbc) (1.5)\n",
      "Requirement already satisfied: subprocess32 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from matplotlib>=1.5.1->pycbc) (3.5.3)\n",
      "Requirement already satisfied: pyOpenSSL in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from lscsoft-glue>=1.59.3->pycbc) (18.0.0)\n",
      "Requirement already satisfied: chardet<3.1.0,>=3.0.2 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from requests>=1.2.1->pycbc) (3.0.4)\n",
      "Requirement already satisfied: idna<2.8,>=2.5 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from requests>=1.2.1->pycbc) (2.7)\n",
      "Requirement already satisfied: urllib3<1.24,>=1.21.1 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from requests>=1.2.1->pycbc) (1.23)\n",
      "Requirement already satisfied: certifi>=2017.4.17 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from requests>=1.2.1->pycbc) (2018.10.15)\n",
      "Requirement already satisfied: pytest>=2.8 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from astropy<3.0.0,>=2.0.3->pycbc) (3.8.2)\n",
      "Requirement already satisfied: setuptools in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from kiwisolver>=1.0.1->matplotlib>=1.5.1->pycbc) (40.4.3)\n",
      "Requirement already satisfied: cryptography>=2.2.1 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pyOpenSSL->lscsoft-glue>=1.59.3->pycbc) (2.3.1)\n",
      "Requirement already satisfied: py>=1.5.0 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pytest>=2.8->astropy<3.0.0,>=2.0.3->pycbc) (1.7.0)\n",
      "Requirement already satisfied: attrs>=17.4.0 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pytest>=2.8->astropy<3.0.0,>=2.0.3->pycbc) (18.2.0)\n",
      "Requirement already satisfied: more-itertools>=4.0.0 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pytest>=2.8->astropy<3.0.0,>=2.0.3->pycbc) (4.3.0)\n",
      "Requirement already satisfied: atomicwrites>=1.0 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pytest>=2.8->astropy<3.0.0,>=2.0.3->pycbc) (1.2.1)\n",
      "Requirement already satisfied: pluggy>=0.7 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pytest>=2.8->astropy<3.0.0,>=2.0.3->pycbc) (0.8.0)\n",
      "Requirement already satisfied: funcsigs in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pytest>=2.8->astropy<3.0.0,>=2.0.3->pycbc) (1.0.2)\n",
      "Requirement already satisfied: pathlib2>=2.2.0 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pytest>=2.8->astropy<3.0.0,>=2.0.3->pycbc) (2.3.2)\n",
      "Requirement already satisfied: asn1crypto>=0.21.0 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from cryptography>=2.2.1->pyOpenSSL->lscsoft-glue>=1.59.3->pycbc) (0.24.0)\n",
      "Requirement already satisfied: cffi!=1.11.3,>=1.7 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from cryptography>=2.2.1->pyOpenSSL->lscsoft-glue>=1.59.3->pycbc) (1.11.5)\n",
      "Requirement already satisfied: enum34 in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from cryptography>=2.2.1->pyOpenSSL->lscsoft-glue>=1.59.3->pycbc) (1.1.6)\n",
      "Requirement already satisfied: ipaddress in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from cryptography>=2.2.1->pyOpenSSL->lscsoft-glue>=1.59.3->pycbc) (1.0.22)\n",
      "Requirement already satisfied: scandir in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from pathlib2>=2.2.0->pytest>=2.8->astropy<3.0.0,>=2.0.3->pycbc) (1.9.0)\n",
      "Requirement already satisfied: pycparser in /home/idies/miniconda3/envs/py27/lib/python2.7/site-packages (from cffi!=1.11.3,>=1.7->cryptography>=2.2.1->pyOpenSSL->lscsoft-glue>=1.59.3->pycbc) (2.19)\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "!{sys.executable} -m pip install pycbc lalsuite ligo-common --no-cache-dir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Catalog of Binary Mergers ###\n",
    "\n",
    "PyCBC provides an [API](http://pycbc.org/pycbc/latest/html/catalog.html) to look at the catalog of binary mergers a few examples below. Some key information, such as the 'chirp' mass of a binary merge can be retrieved."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### What binary mergers are in the catalog? ####"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GW151012\n",
      "GW170608\n",
      "GW170729\n",
      "GW150914\n",
      "GW151226\n",
      "GW170814\n",
      "GW170817\n",
      "GW170104\n",
      "GW170809\n",
      "GW170818\n",
      "GW170823\n"
     ]
    }
   ],
   "source": [
    "from pycbc import catalog\n",
    "\n",
    "### List the mergers in the catalog\n",
    "for merger_name in catalog.Catalog():\n",
    "    print(merger_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How can I get parameters? ####\n",
    "\n",
    "One can also retrieve some of the basic parameters of each source\n",
    "from the catalog directly as follows. Note that all parameters are given\n",
    "in the *source* frame. This means that they include the effect of redshift."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[15.2    7.9   35.7   28.6    8.9   24.2    1.186 21.5   25.    26.7\n",
      " 29.3  ]\n",
      "GW170817: 1.186\n",
      "[u'files', u'distance', u'a_final', u'mass1', u'tc', u'far_gstlal', u'far_pycbc', u'mass2', u'mchirp', u'snr_gstlal', u'redshift', u'far_cwb', u'utctime', u'L_peak', u'sky_size', u'mfinal', u'E_rad', u'chi_eff', u'snr_pycbc', u'snr_cwb']\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAHrlJREFUeJzt3X94XVWd7/H3h1AkQ5CAQGwD9xZHpop2aG3kMo+OJoiWUZRa8Qcjiso81Xn8Pdi5Vp+rKHKrU5Tx6ujcKr9EhohSqiIaERtQZwRaWxuwVEaEuaRIRQk2GKE/vvePvUIP6Umy22afk3P25/U85+nZ6+wf35Wk+3v2WnutrYjAzMzK64B6B2BmZvXlRGBmVnJOBGZmJedEYGZWck4EZmYl50RgZlZyTgRmZiXnRGBmVnJOBGZmJXdgvQPI48gjj4zZs2fvUf7oo49yyCGH1D6gOnO9y8X1LpeprPe6deseioijJluvIRLB7NmzWbt27R7l/f39dHd31z6gOnO9y8X1LpeprLek+/KsV1jTkKSDJd0m6eeS7pT0sVR+uaRfS9qQXvOKisHMzCZX5BXBY8ApETEsaQbwY0nfTZ8tjYhvFHhsMzPLqbBEENm0psNpcUZ6eapTM7NpptC7hiS1SNoAbAVujIhb00cXStoo6WJJTykyBjMzm5hq8TwCSe3AdcC7gd8BvwEOAlYCv4qIj1fZZgmwBKCjo2NBb2/vHvsdHh6mra2twMinJ9e7XFzvcpnKevf09KyLiK7J1qtJIgCQ9FHg0Yi4qKKsG/hARJw+0bZdXV3hu4Z2c73LxfUul/7+foYOO54VfZvZMjTCrPZWli6cw6L5nXu9L0m5EkGRdw0dla4EkNQKnArcJWlmKhOwCLijqBjMzBrN0Mh2lq0aYHBohAAGh0ZYtmqA1esHCztmkX0EM4E1kjYCt5P1EVwPXCVpABgAjgQ+UWAMZmYN5cFH/sTI9p1PKhvZvpMVfZsLO2aRdw1tBOZXKT+lqGOamTW6x3fuotp39C1DI4Ud03MNmZlNIwe1VD8tz2pvLeyYTgRmZtNIx2EH0zqj5UllrTNaWLpwTmHHdCIwM5tG2ltnsHzxXDrbWxHQ2d7K8sVz9+muobwaYtI5M7MyWTS/s9AT/1i+IjAzKzknAjOzknMiMDMrOScCM7OScyIwMys5JwIzs5JzIjAzKzknAjOzknMiMDMrOScCM7OScyIwMys5JwIzs5JzIjAzKzknAjOzkvM01GY2qdXrB1nRt5ktQyPMam9l6cI5NZ0m2YrlRGBmE1q9fpBlqwaeeKD64NAIy1YNADgZNAk3DZnZhFb0bX4iCYwa2b6TFX2b6xSRTbXCEoGkgyXdJunnku6U9LFUfpykWyXdLelrkg4qKgYz239bhkb2qtwaT5FXBI8Bp0TEicA84DRJJwOfAi6OiOOBh4FzC4zBzPbTrPbWvSq3xlNYIojMcFqckV4BnAJ8I5VfASwqKgYz239LF86hdUbLk8paZ7SwdOGcOkVkU63QPgJJLZI2AFuBG4FfAUMRsSOtcj/g3iazaWzR/E6WL55LZ3srAjrbW1m+eK47ipuIIqL4g0jtwHXAR4DLIuKZqfxY4IaImFtlmyXAEoCOjo4Fvb29e+x3eHiYtra2IkOfllzvcnG9y2Uq693T07MuIromW68mt49GxJCkfuBkoF3Sgemq4BhgyzjbrARWAnR1dUV3d/ce6/T391OtvNm53uVShnpXG6fQzt1NX+9q6vH7LvKuoaPSlQCSWoFTgU3AGuDMtNo5wDeLisHMpr/RcQqDQyMEu8cpDI1sr3dopVHkFcFM4ApJLWQJ55qIuF7SL4BeSZ8A1gOXFBiDmU1z441TePARJ4JaKSwRRMRGYH6V8nuAk4o6rpk1lvHGIzy+c1eNIykvjyw2s7oabzzCQS0+PdWKf9JmVlfjjVPoOOzgOkVUPp50zszqanQ8wh53DT1yd50jKw8nAjOru0XzO/cYoNbf70RQK24aMjMrOScCM7OScyIwMys5JwIzs5JzIjAzKzknAjOzknMiMDMrOY8jMLMJVZsi2g+laS5OBGY2rtEpokdnBx2dIhpwMmgibhoys3GNN0X0ir7NdYrIiuBEYGbjGm+K6PHKrTE5EZjZuMabInq8cmtMTgRmNq7xpoheunBOnSKyIriz2MzGNd4U0e4obi5OBGY2oWpTRFtzcdOQmVnJ+YrArEF4YJcVpbArAknHSlojaZOkOyW9N5WfL2lQ0ob0enlRMZg1i9GBXYNDIwS7B3atXj9Y79CsCRTZNLQDOC8ing2cDLxT0gnps4sjYl563VBgDGZNwQO7rEiFNQ1FxAPAA+n9NkmbAF/Hmu0DD+yyItWks1jSbGA+cGsqepekjZIulXR4LWIwa2Qe2GVFUkQUewCpDbgZuDAiVknqAB4CArgAmBkRb6uy3RJgCUBHR8eC3t7ePfY9PDxMW1tbkeFPS653uQwPD7Oj5SkMPjzCror/rwdIdB7eSnvrjDpGV5wy/76nqt49PT3rIqJrsvUKTQSSZgDXA30R8Zkqn88Gro+I5060n66urli7du0e5f39/XR3d09JrI3E9S6X0XqX7a6hsv++p4KkXImgsD4CSQIuATZVJgFJM1P/AcCrgTuKisGsmXhglxWlyHEELwDeBAxI2pDKPgScJWkeWdPQvcDbC4zBzMwmMWkikPQCYENEPCrpbOB5wGcj4r6JtouIHwOq8pFvFzVrAGVriiqzPHcNfRH4o6QTgX8E7gO+UmhUZlZXHsBWLnkSwY7IepTPILsS+CxwaLFhmVk9eQBbueTpI9gmaRlwNvAiSS1Ac96vZmaAB7CVTZ4rgtcDjwHnRsRvyEYHryg0KjOrKw9gK5c8iWAbWZPQjyT9BTAPuLrYsMysnvxksnLJkwhuAZ4iqRO4CXgrcHmRQZlZfS2a38nyxXPpbG9FQGd7K8sXz/VdQ00qTx+BIuKPks4FPhcR/1QxLsDMmpQHsJVHnisCSfor4I3Ad1JZywTrm5lZA8mTCN4LLAOui4g7JT0DWFNsWGZmViuTNg1FxC1k/QSjy/cA7ykyKDMzq508U0wcRTai+DnAwaPlEXFKgXGZmVmN5Gkaugq4CzgO+BjZRHG3FxiTmZnVUJ5E8LSIuATYHhE3p4fInFxwXGZmViN5bh/dnv59QNIrgC3AMcWFZGZmtZQnEXxC0mHAecDngKcC7y80KjMzq5k8dw1dn94+AvQUG46ZmdXauIlA0v+ZaMOI8C2kZmZNYKIrgneQPU/4GrJ+gWpPGzMzswY3USKYCbyWbBrqHcDXgGsj4uFaBGZmZrUx7u2jEfG7iPjXiOgB3gK0A3dKelOtgjMzs+LlGVn8POAs4KXAd4F1RQdlZma1M1Fn8ceA04FNQC+wLCJ25N2xpGPJHnL/dGAXsDIiPivpCLJmptlko5Rf5+Yms/xWrx9kRd9mtgyNMKu9laUL53i6aNsvE40s/l/AYcCJwHLgZ5I2ShqQtDHHvncA50XEs8lGIr9T0gnAB4GbIuJ4sgfdfHC/amBWIqvXD7Js1QCDQyMEMDg0wrJVA6xeP1jv0KyBTdQ0dNz+7DgiHgAeSO+3SdpE9rzjM4DutNoVQD/wP/fnWGZlsaJvMyPbdz6pbGT7Tlb0bfZVge2zcRNBRNw3VQeRNBuYD9wKdKQkQUQ8IOnoqTqOWbPbMjSyV+VmeSgiij2A1AbcDFwYEaskDUVEe8XnD0fE4VW2WwIsAejo6FjQ29u7x76Hh4dpa2srLvhpyvUul8p6b/7NNh7fuWuPdQ5qOYA5Tz+01qEVyr/v/dfT07MuIromWy/PXEP7TNIM4FrgqohYlYoflDQzXQ3MBLZW2zYiVgIrAbq6uqK7u3uPdfr7+6lW3uxc73KprPdQ6iOobB5qndHC8sVz6W6ypiH/vmtn0mmoJZ0uKc901WO3E3AJsCkiPlPx0beAc9L7c4Bv7u2+zcpq0fxOli+eS2d7KwI621tZvniu+wdsv+S5IngD8FlJ1wKXRcSmnPt+AfAmYEDShlT2IeCTwDWSzgX+i2z0spnltGh+p0/8NqXyzD56tqSnkg0qu0xSAJcBV0fEtgm2+zHjz0/0kn0J1szMpl6uJp+I+ANZW38v2RxEryYbV/DuAmMzM7MayNNH8CpJ1wE/BGYAJ0XE35ANNPtAwfGZmVnB8vQRvAa4OCJuqSyMiD9KelsxYZmZWa1MeEUgqQXoHJsERkXETYVEZWZmNTNhIoiIncAf0zOLzcysCeVpGvoT2S2gNwKPjhb6UZVmZs0hTyL4TnqZmVkTyjOO4IpaBGJmZvWR5wllx5M9j+AE4ODR8oh4RoFxmZlZjeQZUHYZ8EWyB830kD117MoigzIzs9rJkwha022iioj7IuJ84JRiwzIzs1rJdddQmn30bknvAgYBP0zGzKxJ5LkieB/wZ8B7gAVkM4qeM+EWZmbWMPLcNXR7ejsMvLXYcMzMrNbGTQSSvg2M+xzLiHhVIRGZmVlNTXRFcFH6dzHwdOCrafks4N4CYzIzsxoaNxFExM0Aki6IiBdVfPRtSVUnoTMzs8aTp7P4KElPDB6TdBxwVHEhmZlZLeW5ffT9QL+ke9LybODthUVkZmY1leeuoe+laSaelYruiojHig3LzMxqJc8VAWTjB2an9U+URER8pbCozMysZvI8s/hKsjuIXgg8P726cmx3qaStku6oKDtf0qCkDen18v2I3czMpkCeK4Iu4ISIGHdMwTguBz5PNkldpYsj4qI9Vzczs3rIc9fQHWTjCPZKes7x7/c6IjMzq6k8VwRHAr+QdBvwRCfxfowsfpekNwNrgfMi4uF93I+ZmU0BTdbiI+nF1cpHB5xNsu1s4PqIeG5a7gAeIpu64gJgZkS8bZxtlwBLADo6Ohb09vbusc7w8DBtbW2ThdF0XO9ycb3LZSrr3dPTsy4iJu3TJSImfQEdwOnpdXSebdJ2s4E79vazsa8FCxZENWvWrKla3uxc73JxvctlKusNrI0c59g8dw29DrgNeC3wOuBWSWfuQ3JC0syKxVeT9T+YmVkd5ekj+DDw/IjYCiDpKOAHwDcm2kjS1UA3cKSk+4GPAt2S5pE1Dd2LRyibmdVdnkRwwGgSSH5HjruNIuKsKsWX5A3MzMxqI08i+J6kPuDqtPx64LvFhWRWrNXrB1nRt5ktQyPMam9l6cI5LJrfWe+wzOomz1xDSyUtJhtZLGBlRFxXeGRmBVi9fpBlqwYY2b4TgMGhEZatGgBwMrDSmjQRpGmnb4iIVWm5VdLsiLi36ODMptqKvs1PJIFRI9t3sqJvsxOBlVaekcVfB3ZVLO9MZWYNZ8vQyF6Vm5VBnkRwYEQ8PrqQ3h9UXEhmxZnV3rpX5WZlkCcR/FbSE9NJSDqDbHSwWcNZunAOrTNanlTWOqOFpQvn1Ckis/rLc9fQO4CrJP0L2f3/9wNvLjQqs4KM9gP4riGz3fLcNfQr4GRJbWRzE20rPiyz4iya3+kTv1mFPFNMdEi6BPh6RGyTdIKkc2sQm5mZ1UCePoLLgT5gVlr+JfC+ogIyM7PaypMIjoyIa0i3kEbEDrJbSM3MrAnkSQSPSnoaWUcxkk4GHik0KjMzq5k8dw39A/At4M8l/QQ4CtinaajNzGz6yXPX0M/SU8rmkM01tDkithcemZmZ1cS4TUOSni/p6fBEv8AC4ELg05KOqFF8ZmZWsIn6CP4v8DiApBcBnwS+QtY/sLL40MzMrBYmahpqiYjfp/evJ5t++lrgWkkbig/NzMxqYaIrghZJo4niJcAPKz7L08lsZmYNYKIT+tXAzZIeAkaAHwFIeia+fdTMrGmMmwgi4kJJNwEzge9HRKSPDgDeXYvgzMyseBM28UTET6uU/bK4cMzMrNbyjCzeJ5IulbRV0h0VZUdIulHS3enfw4s6vpmZ5VNYIiCbrO60MWUfBG6KiOOBm9KymZnVUWGJICJuAX4/pvgM4Ir0/gpgUVHHNzOzfIq8IqimIyIeAEj/Hl3j45uZ2RjafTNQATuXZgPXR8Rz0/JQRLRXfP5wRFTtJ5C0BFgC0NHRsaC3t3ePdYaHh2lraysg8unN9S4X17tcprLePT096yKia7L1aj0w7EFJMyPiAUkzga3jrRgRK0lTWXR1dUV3d/ce6/T391OtvNm53uXiepdLPepd66ahbwHnpPfnAN+s8fHNzGyMIm8fvRr4D2COpPvTc44/CbxU0t3AS9OymZnVUWFNQxFx1jgfvaSoY5qZ2d6rddOQmZlNM04EZmYl50RgZlZyTgRmZiXnRGBmVnJ+0pgBsHr9ICv6NrNlaIRZ7a0sXTiHRfM76x2WmdWAE0ET2N+T+Or1gyxbNcDI9p0ADA6NsGzVAICTgVkJuGmowY2exAeHRgh2n8RXrx/MvY8VfZufSAKjRrbvZEXf5imO1symIyeCBjcVJ/EtQyN7VW5mzcWJoMFNxUl8VnvrXpWbWXNxImhwU3ESX7pwDq0zWp5U1jqjhaUL5+xXbGbWGJwIGtxUnMQXze9k+eK5dLa3IqCzvZXli+e6o9isJHzXUIMbPVnv762fi+Z3+sRvVlJOBE3AJ3Ez2x9uGjIzKzknAjOzknMiMDMrOScCM7OScyIwMys5JwIzs5JzIjAzK7m6jCOQdC+wDdgJ7IiIrnrEYWZm9R1Q1hMRD9Xx+GZmhpuGzMxKr16JIIDvS1onaUmdYjAzM0ARUfuDSrMiYouko4EbgXdHxC1j1lkCLAHo6OhY0Nvbu8d+hoeHaWtrq0XI04rrXS6ud7lMZb17enrW5emDrUsieFIA0vnAcERcNN46XV1dsXbt2j3K+/v76e7uLi64acr1LhfXu1ymst6SciWCmjcNSTpE0qGj74GXAXfUOg4zM8vU466hDuA6SaPH/7eI+F4d4jAzM+qQCCLiHuDEWh/XzMyq8+2jZmYl50RgZlZyTgRmZiXnRGBmVnJOBGZmJedEYGZWck4EZmYl50RgZlZyTgRmZiXnRGBmVnJOBGZmJedEYGZWck4EZmYl50RgZlZy9XgeQU2sXj/Iir7NbBkaYVZ7K0sXzmHR/M56h2VmNu00ZSJYvX6QZasGGNm+E4DBoRGWrRoAcDIwMxujKZuGVvRtfiIJjBrZvpMVfZvrFJGZ2fTVlIlgy9DIXpWbmZVZUyaCWe2te1VuZlZmTZkIli6cQ+uMlieVtc5oYenCOXWKyMxs+mrKzuLRDmHfNWRmNrm6JAJJpwGfBVqAL0fEJ6f6GIvmd/rEb2aWQ82bhiS1AP8C/A1wAnCWpBNqHYeZmWXq0UdwEvCfEXFPRDwO9AJn1CEOMzMDFBG1PaB0JnBaRPxdWn4T8D8i4l1j1lsCLAHo6OhY0Nvbu8e+hoeHaWtrKz7oacb1LhfXu1ymst49PT3rIqJrsvXq0UegKmV7ZKOIWAmsBOjq6oru7u49Nurv76daebNzvcvF9S6XetS7Hk1D9wPHViwfA2ypQxxmZkZ9moYOBH4JvAQYBG4H/jYi7pxgm98C91X56EjgoSLinOZc73JxvctlKuv93yPiqMlWqnnTUETskPQuoI/s9tFLJ0oCaZuqFZG0Nk/7V7NxvcvF9S6XetS7LuMIIuIG4IZ6HNvMzJ6sKaeYMDOz/Bo9EaysdwB14nqXi+tdLjWvd807i83MbHpp9CsCMzPbTw2bCCSdJmmzpP+U9MF6x1MUSZdK2irpjoqyIyTdKOnu9O/h9Yxxqkk6VtIaSZsk3Snpvam82et9sKTbJP081ftjqfw4Sbemen9N0kH1jrUIklokrZd0fVpu+npLulfSgKQNktamspr/nTdkIijZxHWXA6eNKfsgcFNEHA/clJabyQ7gvIh4NnAy8M70+232ej8GnBIRJwLzgNMknQx8Crg41fth4Nw6xlik9wKbKpbLUu+eiJhXcctozf/OGzIRUKKJ6yLiFuD3Y4rPAK5I768AFtU0qIJFxAMR8bP0fhvZyaGT5q93RMRwWpyRXgGcAnwjlTddvQEkHQO8AvhyWhYlqPc4av533qiJoBP4fxXL96eysuiIiAcgO2kCR9c5nsJImg3MB26lBPVOzSMbgK3AjcCvgKGI2JFWada/9X8G/hHYlZafRjnqHcD3Ja1LE21CHf7OG/UJZbkmrrPGJqkNuBZ4X0T8IfuS2NwiYicwT1I7cB3w7Gqr1TaqYkk6HdgaEeskdY8WV1m1qeqdvCAitkg6GrhR0l31CKJRrwjKPnHdg5JmAqR/t9Y5niknaQZZErgqIlal4qav96iIGAL6yfpI2tMcXdCcf+svAF4l6V6yZt5TyK4Qmr3eRMSW9O9WssR/EnX4O2/URHA7cHy6q+Ag4A3At+ocUy19CzgnvT8H+GYdY5lyqX34EmBTRHym4qNmr/dR6UoASa3AqWT9I2uAM9NqTVfviFgWEcdExGyy/8s/jIg30uT1lnSIpENH3wMvA+6gDn/nDTugTNLLyb41jE5cd2GdQyqEpKuBbrIZCR8EPgqsBq4B/hvwX8BrI2Jsh3LDkvRC4EfAALvbjD9E1k/QzPX+S7LOwRayL2nXRMTHJT2D7JvyEcB64OyIeKx+kRYnNQ19ICJOb/Z6p/pdlxYPBP4tIi6U9DRq/HfesInAzMymRqM2DZmZ2RRxIjAzKzknAjOzknMiMDMrOScCM7OScyKwXCTtTDMk3iHp65L+bJz1bhi9F77WJM2WFJIuqCg7UtJ2SZ8v8LjfSLcCIultaTbJjelnNeEcWJIul3TmROtMsv1bUp1fUlH26lS2z/vdH5J6JR1fj2PbvnEisLxG0gyJzwUeB95R+aEyB0TEy9Oo2ClTMbo0j3uA0yuWXwvcOZXxVJL0HKAlIu5JE6d9GHhhRPwl2ajgjVN8vGo/iwHgrIrlNwA/n8rj7qUvks0bZA3CicD2xY+AZ6Zv4JskfQH4GXBsml/9yPTZXZK+nL4ZXyXpVEk/SfOsnwQg6SRJ/57mof93SXNS+VvSlce3ySblurLy23Xa36uqxDYCbJI0OqXv68kG54xu90plc9yvl/QDSR2p/MXpimdD+uxQSTMl3VJxJfTXVY73RnaP/Dwa2AYMA0TEcET8Ou1/nqSfpiuF61RljnlJH5F0ezrWyjTCGkn9kv63pJvJpmqu9vs4SdIMZfMzPRPYkGO/75H0ixRT73g/h1S+NO1jo3Y/J+EQSd9R9vyEOyS9viKeU/cygVs9RYRffk36AobTvweSnfj+HphNNvL35Ir17iUbBT2b7LkCc8m+cKwDLiWbTOwMYHVa/6nAgen9qcC16f1byOaUOiItv7him8OAX49uV3Hs2WRD9F8FXEQ2P81NaV+fT+sczu6BlH8HfDq9/zbZBGAAbame5wEfTmUtwKFVfi43A3Mr1ukjGw16GfDKivU2Ai9O7z8O/HN6fzlwZnp/RMX6V45uTzbn0BfG+b28Bfg88BmyK6E3ko0+z7PfLcBT0vv2CX4OLyN7jq7S7/J64EXAa4AvVez7sIr3NwIL6v1361e+l68ILK9WZdMjryU70V2Syu+LiJ+Os82vI2IgInaRNc/cFNlZYoDspA3ZSf3ryp7AdjHwnIrtb4w0tD4ibia7CjmarBnk2tg9RfFY3wNemtb72pjPjgH6JA0ASyuO9xPgM5LeQ3ZS3EE2p9VbJZ1PdrLfVuVYM4Hfphh3kj1E6Ezgl8DFks6XdFja581pmyvITqRj9aSrlQGyidcqfxZj6zFWL1mT0BuAq3PudyNwlaSzyZL2eD+Hl6XXerIrv2cBx5P9Hk+V9ClJfx0Rj1Qccyswa5KYbZpwIrC8RvsI5kXEuyN7IBDAoxNsUzkvzK6K5V3sngL9AmBNZH0PrwQOrthm7L6vJPvG+1ayb9xVpdjWkX2jv3bMx58juzqYC7x99HgR8UmyK4RW4KeSnhXZQ4FeBAwCV0p6c5XDjVTGHJnbImI52Un5NePFWUnSwcAXyL7FzwW+xMQ/i7F1vg14LnBkRPwy535fQfakvwXAOkkHVvs5kF0JLK/4/T8zIi5Jx1lAlhCWS/pIRUgHp5+NNQC34Vm9HUZ2ooWsmWMilwO3Ab+JiMk6gD8N3BwRv9OTn2NQebzRGR6R9OcRMQAMSPor4FmSRoDBiPiSstkhnwd8ZcxxNpG1yd8raRbw9EhPVyN73OR9EfGIpIfTt+YfAW8ia1KqNHpyfii185/J7qdz5bUM+FOe/Uo6ADg2ItZI+jHwt0CbpKeN/TmQNXddIOmqiBiW1AlsJzt//D4ivippmCf//v6CAjvpbWo5EVi9/RNwhaR/AH440YoR8aCkTWSzr04oJYpqJ6LzyZqiBoGfAsel8vdJ6gF2Ar8Avkv2jX6ppO1kHcDVrgi+QzY77A/IHi15UUoIfyJrMhq9u+oc4F+V3XZ7D9lVTWW8Q5K+RPbt+l6yZqm9EhHfrVI23n5bgK+mZiuRPRt4SNIFY38OEfGYpGcD/5GS6jBwNlkCXCFpF1li+HuA1AE/EukpWzb9efZRaxjpJDoAPG9Me3TdKHtuwBqyDtad9Y5nOpD0fuAPEXHJpCvbtOA+AmsIkk4F7gI+N12SAEBEjJDdpdOMz9PdV0Psfvi6NQBfEZiZlZyvCMzMSs6JwMys5JwIzMxKzonAzKzknAjMzErOicDMrOT+P/aeFCDIWRHHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import pylab\n",
    "\n",
    "# Either from the catalog as a whole\n",
    "c = catalog.Catalog()\n",
    "mchirp = c.median1d('mchirp')\n",
    "print(mchirp)\n",
    "\n",
    "# or from a specific merger\n",
    "m = catalog.Merger(\"GW170817\")\n",
    "mchirp_gw170817 = m.median1d('mchirp')\n",
    "print('GW170817: {}'.format(mchirp_gw170817))\n",
    "\n",
    "# print parameters that can be read\n",
    "print m.data.keys()\n",
    "\n",
    "pylab.scatter(c.median1d('mass1'), c.median1d('mass2'))\n",
    "pylab.xlabel(\"Primary Mass (Solar Masses)\")\n",
    "pylab.ylabel(\"Secondary Mass\")\n",
    "pylab.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Transform Mass Parameters into the Detector Frame\n",
    "\n",
    "By default the above interface returns parameters in the *source* frame. Due to cosmological redshift, gravitational-waves are stretched as they travel. This causes the observed waveform to be different in the detector frame. This corresponds to an observed change in the mass parameters (for example). However, the relationship is fairly straighforward."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chirp Mass of GW150914\n",
      "Source Frame: 28.6 Solar Masses\n",
      "Detector Frame: 31.174 Solar Masses\n"
     ]
    }
   ],
   "source": [
    "m = catalog.Merger('GW150914')\n",
    "source_mchirp = m.median1d('mchirp')\n",
    "redshift = m.median1d('redshift')\n",
    "det_mchirp = source_mchirp * (1 + redshift)\n",
    "\n",
    "print('Chirp Mass of GW150914')\n",
    "print('Source Frame: {} Solar Masses'.format(source_mchirp))\n",
    "print('Detector Frame: {} Solar Masses'.format(det_mchirp))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Accessing LIGO/Virgo data\n",
    "\n",
    "In this section, we will look into how to read detector data from the LIGO and Virgo instruments using the PyCBC API. It is possible to both get data around specific events, and also from the full data sets which have been released which cover the S5/S6/O1 LIGO observing runs. Data will be returned as [pycbc TimeSeries objects.](http://pycbc.org/pycbc/latest/html/pycbc.types.html#pycbc.types.timeseries.TimeSeries)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Getting Data Around  Specific Binary Merger in the Catalog\n",
    "\n",
    "One can directly retrieve data around a specific even. Typically this data is centered on the event, though restrictions may apply which have not allowed this. This method by default gets the smallest version of the dataset. If additional data or specific versions are required, please see the following two additional ways to access data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Duration: 32.0s Start: 1126259447 End: 1126259479\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAERCAYAAAB7FtAjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztnXecFPX5xz/PFTiOIu1AEPTovYgnKigqoKBYYosa/VliYkyMLVGDPXbUWBKNGjRqYu8aRSmKgooIR+9Nej06R7n6/P7Y2bvZ2ZnZmd3vlL193q8XL25nZ2eenZ35Pt/vU4mZIQiCIAhZQQsgCIIghANRCIIgCAIAUQiCIAiChigEQRAEAYAoBEEQBEFDFIIgCIIAIA0VAhG9QkTbiGihouONJ6LdRPS5YftQIppNRAuJ6D9ElKPifIIgCGEl7RQCgNcAjFR4vCcA/J9+AxFlAfgPgEuYuTeAtQCuVHhOQRCE0JF2CoGZpwLYqd9GRJ20mf4sIvqOiLq7ON7XAPYZNrcAUMbMy7XXkwBckIrcgiAIYSftFIIFYwHcwMzHALgVwPMpHm87gFwiKtJeXwigfYrHFARBCDVpbxcnokYABgF4n4iim+tr750P4AGTj21k5hFWx2RmJqJLADxNRPUBTARQqVRwQRCEkJH2CgGRVc5uZu5vfIOZPwLwUTIHZeYfAZwEAER0OoCuqQgpCIIQdtLeZMTMewGsJqKLAIAi9Ev1uETUSvu/PoC/AHgx1WMKgiCEmbRTCET0NoAfAXQjog1EdA2AywBcQ0TzACwCcK6L430H4H0Aw7TjRU1JtxHREgDzAXzGzJOVfhFBEISQQVL+WhAEQQDScIUgCIIgeENaOZVbtmzJhYWFQYshCIKQVsyaNWs7Mxck2i+tFEJhYSGKi4uDFkMQBCGtIKK1TvYTk5EgCIIAIGCFQES3ENEirYDc20SUF6Q8giAImUxgCoGIjgBwI4AirYBcNoBLgpJHEAQh0wnaZJQDoIFWWjofwKaA5REEQchYAlMIzLwRwN8ArAOwGcAeZp5o3I+IriWiYiIqLikp8VtMQRCEjCFIk1EzRDKKOwBoC6AhEV1u3I+ZxzJzETMXFRQkjJoSBEEQkiRIk9FwAKuZuYSZKxApQjcoQHkEQRAymiAVwjoAxxNRPkXqVg8DsCRAeQShTrFw4x58tXhr0GIIaURgiWnM/BMRfQBgNiK9BuYg0uhGEAQFnPXs9wCANWNGBSyJkC4EmqnMzPcBuC9IGQRBEIQIQYedCoIgCCFBFIIgCIIAQBSCIAiCoCEKQRAEQQAgCqFO0f2eL/HQ54uDFkMQhDRFFEId4lBFNV7+fnXQYgiCkKaIQhAEQRAAiEIQBEEQNEQhCIIgCABEIQiCIAgaohAEQRAEAKIQBAEA8JcP5uON6WuDFkMQAkUUQh2hqpqDFiGtebd4Pe7+ZGHQYghCoIhCqCOUVVYFLYIgCGmOKIQ6AssCQRCEFBGFUEeoFo0gCJ7x6BdLcOPbc4IWw3NEIdQRsoiCFkEZt74/D8c/8nXQYgRG/wcm4h7xZ4SKf039Gf+btyloMTxHFEIdoS4phA9mbcCWvYeCFiMwdh+owOsS8SQEQKAKgYiaEtEHRLSUiJYQ0QlBypPOMNLbZHSwvAoHy/13jK/Yug+Fo8f5fl5BHVOWl+CbZduCFqNOEGhPZQB/BzCemS8konoA8gOWJ21JdxdCj3vHIzuLsOqRM30757It+zDimam+nU/whitfmQEAWDNmVMCSpD+BrRCIqAmAIQD+DQDMXM7Mu4OSRwgev3Mp0kEZPP/tSvz1f4uCFkPIEII0GXUEUALgVSKaQ0QvE1FD405EdC0RFRNRcUlJif9SpglpvkAQLHh8/DK8Nm1N0GJkLBVV1RllUgxSIeQAGADgBWY+GsB+AKONOzHzWGYuYuaigoICv2UMHRt2HcCzX68AG2xExtd1gcLR4/DDyu1BiyFkMJl2/wWpEDYA2MDMP2mvP0BEQQg2/OY/xXhy0nKs23kgaFF84ZM5G4MWQchgyiqrgxbBVwJTCMy8BcB6IuqmbRoGQBoCJ2Dpln0AAKO5ve6tDyJMW7UjaBFCQfGanUGLkJFUVIlC8JMbALxJRPMB9AfwSMDypA37yypjXtdBixEAYOPug0GLEAoufPFHrNuRGavCZDn+ka9xoLwy8Y4uqKyqow+WBYGGnTLzXABFQcqQrsTlodXR+7btYXlBixAa9pVVBC1CqNmy9xBWbitF33ZNlR1TVghCWkCI1QjpnphmxXEdWwQtQmioq6tAlai+RpUZVlZeFEKakmkzF0FwgurhuzLDnjNRCGnKok17Y17v3F8ekCTeotomHHb04cPfLFVTjqEuhiRbofq7lmeYD0EUQpqi9yFMW7kdQ5+cEpwwKaJf7fz7+9Ux701YtNVvcQJFP55d/dpM5ccUnLGjtAxllVUZt0IIupaRkCR6D8IbP6VvZcxP527ETe/MrXn94OfeRx4zMyYs2uL5eZJBxu7UUHX9jnnoK+RmEypkhSCkG18sCOfg5oQ3p6/z/ZxfLdmG696Y7ft5neCFeSeThjSV1y/TlAEgCiFtGf3RArz+45qgxUiZGQEkXO3cX+b7OZ3ixRCUWT6EoCVIb0QhpBlHtaitEH7Pp1IFMxlmrtkVtAiW2LVCTXawy6Qx8lBF6jZ/vxVodTWjPCQlMkQhpBk7SutmNJGfzF4bXoVgNxYlm2tSeihzIrXeVOBP83uVcf1bs9H17i/9PakFohDSiBvenoPSsvR7uH9ctQMfz9lg+l7zhvVsP9uovvq4h5+371d+zFRYuW0fOt/5BdZ7VLAwk3JWVHxXu1WaF3y5MDw+QFEINizatAeFo8dhVUlp0KIAAD5z2OT7kS+W+N5sxo5LX5qOW96dZ/regCPtywyM6HW4cnlOCFn285MTl6OymvHfH9d4Mjv93RuzQnU/eMkJnVo63vdQhXnL1gy5VKaIQrDh1vfnAwBe/zG9wjrHTv0Z7xevD1oMRyQaAFWX5Fi5bR8OWgwEQRGdIX4+f7Pt7HRikjkZc9btxi3vzk28Yx2gR5vGjvb7eM4GdL9nPOatj23SWF3NmLg4mBn7HR/ND6SvuB5RCDYs2RzJBp6ZhqWHnUTv3P7BPFz04jQfpLHG78nY8KemYu765Du1MjNWbN2nUCL9se2vx4w1O5N2eE5aXDcT/DbsSs7MFl2xvjhlVcz216evxR/fmpOyXMnw9oz16Hnf+EDOHUUUggOMZSLSgfx62ZbvVVUz7vx4Ad4r3hBIxM1j45dipNbPePOeQ76fPxU+nrMRpz09VVlZCT1b9h7CkxOXWb4/Y/VOvDMzuZVf2FZFKhjz5VKc+Ng3KR3DmGvg1CzrFUGHzYpCSBMWbtzjav+G9aydsUs278VbP/mfEBblhW9X1TT6ia7C0oWo2eabZRGFcO1/i5V2dXv1hzW27xc7VOArt3mzigkT/5q6Kn6jywE1y1BGvjjEEWh+IAohDZi1difOevZ7dx8y9kvQUVYZntlir7ZNghbBFeO1khfvajP1iYu34mYf7fNOZ7DDn5rqsSTBkxXXFMS9CXKiwZSWY9QQGYYoBAveC5FTdv1O913DChrVt3wvTKUuChpbywkgtFlVKnvtdipo6Hjf8gwKIU2E2djtxOSy3MYH1MDG1JoJiEKw4PYP5gctQg17D7nvlNW5VSPL93q0UTsrf2P62jpbftsP6uWoHYQypVSFsUkU4Cwq7fSnrVdP+zIoic+MwBUCEWUT0Rwi+jxoWcLKC9+a2EoTYPdYqGxLuWLrPtz9yULc+HYwkRmq+ZML84+qEMHdB9Qq0wzRB/FtZBHsd7/ghWlx5dvTjcAVAoCbACwJWogwozpChMyepCTZczCyevl+5XbXn62sqg7d4PWRCwexKvON6kirkF1SzzCz97vNMjaaLLsf7iyPwYxZa3f5Ur7dSwJVCETUDsAoAC8HKUciEpVX8JrdB9ybjK5+dSae/3al6XvZCh1n41NIu+9815eYsrzEdp8wD27JmGbKK6vjkqFU43fphaBQ4VQ2/oYndnae6VwXCXqF8AyA2wFYTrWI6FoiKiai4pIS+8HDK45ub19eISgSjeuPjzePac9W+KtnshMumRIHD41bjHP/+YOnYaGZohD2mdX1cvnVjZdK4eI5KdooNOcmQ2AKgYjOArCNmWfZ7cfMY5m5iJmLCgoKfJIulqw6Foqm0mSUl1v3FEK1w5E+mfpA/9XKoIyd+rNnSiFD9IEpbkudqFKeTu+ZRAwOeIUS5AphMIBziGgNgHcADCWiNwKUx5KgY5N7Ko4KUvl9gjaneYHTzPRU6iy9V7whqVyBIV0TT4r8VAiLNu3Bv6a4D3rwCrffvUmD3JjXyU6WSsvVRCcFrcwDUwjMfAczt2PmQgCXAJjMzJcHJU+YGd6zten2ZO8ds3C9ZAl6iesF+yzCfJdtiZ3R11ccLuqENk0SX28/TUaj/vE9Hv1yqW/nS4TbifqFA9rFvHZbESBKlct2mwfKK019UMmeXxVB+xDSgqDrlScadEf0MlcYVqisINosP7JCsMt7SIWV20rx2g+rUTh6nG1CkUqsks5GPGOY0YfUNBNSsVKmcPQ43P3JAtt93Dr6dxpCfqP3s1tyc5wPpdtLy9Dz3gl43iScfJlP97gVoVAIzPwtM58VtBxhJdF83izaworyymp8t8J9iKgV0cevgUe+hAUb9+Cvn0VC+Z6bbB41pRqnl1N1aW5VeL1CKC2rjGv56Fcy3BvT7WtwuZXCaD5t17yByyPYc/WrM/CDISR75bZIf5UnJlgXMgyKUCiEuoCXD4TVkZM55aNfLrG9EaurGSOfmYrxCzc7Op6fjVdWJ9HprLqa8cK3qyzNQGY4dZQH0Uhl057EZUy8Hpt73zcBl740PWZbKtdiz4EK7CgtS1GqCG6/u9uAEavnfOve+FwSZsY3y0pw2cs/xWz/bkUw0ZJOyGiFsHnPQcuuSW4Yv3ALOtzxRY3m9xs3frBEMu4vr8TSLftw3RuzHR0v+oD87ENXubU73CuEz+ZvwmPjl+Lcf/7g+DNO2jB2a904kBIR363YntDO7IdcswxVQVOZGPR7YCKOeeirVEXScCdHtvHhSdSwyeT9cfM3Y9iTU+K2V1pcE7uf55x+be0F8JiMVggnPDoZv/1vccL9EiWrvDE9Eko4f4M3CUdz16k77sZd9jNMtw92dPf95VUoHD0OB3TRFge1barYm0SdmacmLQcA/FziXJlUag7C9TsP4LSnpsR8pyindC+wHDvem7ken85VVxLbiL7yLTNj6ZbYqKggVi5O7htmxugP5yeVmOe0l7hbXeg24s7s8F8tMW8+ZDWxCHOLzoxWCAAc2dMPt3Hqvl+8vqZsg1ed1d5NUHl11Tbng12iBvNuFUKl4ab/ZE5teWbjMtrrmeuhiipMNWQ+J1NvKDr4nPT4N1ixrRQvfxdfn6ZZfj3Lwef2D+fjpnf8KYn95k/rMPKZ72Iyvt1e54Ub9+D7FdtT6iS3cXfizmU795fjnZnrcdWrMxwdk5nx7+9XY9u+Q6Ymv95HxIdju73DjCajZO7Qjy3KnVRUWq0QwqsRMlYhuEkk+Wn1jriQw72HKrD7QDlu01VFDaKReeHocdiuyP4KuCtr8cmcjfiVwT5aVV2rIOoZIi+8vjz3f7YYV7wyA4t1eQTJnPOGt+fgiwW1PhQr30UYHuy7P1kIINI9LIrb73zWs9/j8n//hF+4MKsZ2eWgvErUhLLrQAVOfGxywv2Xby3Fg58vxsCHvzZ9tuqZpNy7/UnsGkmZ4eY3r6i2WiEEf99YkbEKYft+54Po+p0H40IOBzwwCf0fmBSz7WBFNZ6bvML3gUJl5vHCjc47mJk1hqmvc8jmZMfKVVXNSnw2VkQzf/XmhWR/i/d1qzKrGWCYHuvSstoB2Yv7772Z67F6+37LYyf6XQtHj8Nxj3xd83pDAtMlgBhHs9nYenzHFnHb3EZ+tWjkLszUzdGtTEZ2P0/Q91TGKoRD5alVqjRzGH02bxP+NnE5vlqivt+uHSoTqaf/vCOlz9uFn1Yzo/s9yTcRb5igblK0CKD+QdyRZJ+Gb5bZR4KMnfpz4FmlenJ1s2UvxLr9w/k497nvLb+zqlLgeibr+labzbbNVrNufxOvaheNX7gFJzwavwqqrmbPTMsqyFiFUOXh01ys+AdvVN9+WWtn5nHbuCbVUsz6gcH4DKe6VP79KZ1s31+hRVC9m2Qjeitys+Ov78795aHKQ9BLqNokETWv7j1UiUMW7VcPKewgF0VfXlw/0Yh2MzT7mmWV1SmZbhOtrszePu/oI+K2XfeGeYm2wY9NxrwN1lFiQZshM1YhxIWbKWT6arUK4cZhnW3ft0tMG/DgJGzb53yQdxJyacftH86viSwyKt1BYxLbje1IpBijqP5pm+Tlmm43G3e8jC6yY5Uuikr1mKKfnfe8d4LpPmUemAL13cv0CuH2D+ZjgcWgeuv78zDSmFGuELNJgJVJ0QzVvS9Uk7EKwcrhowKV9e7X7zyAz+bZJ4klGgDdrBKMGaipYHTcJ9PXIRlystTe1lZmJ+Ns7uM5G3yLLgKAMovZut0KYdbaXXGZs3rMZqiVDur0qOyxUXNe3f1jlGDH/jJMXGxeUmbFtlIUjh6Hx8bH11gyPgv6r3vJ2B/xkklEWcx5S1Prbme22gwTGasQvIwIOr5jc2XHOu/5H7AgQSJSotIV89fbf75w9Dhs00JEnXYBe2eGfQkBIJioKyC238NdZ/ZI+XiDOsU7L4HYwWTwmMmYozBfxAnd7jb3x9itEC54YVpc5myizzpRCI0tVlGpoA9pNiqqRZv2YvlW+2RIY+vZ0rJKDHgwNhBk/oY9NWbO6T8nXtmnOmGqSHAtgzZCZpRCWLhxDyYtjiSROLnJkyVXYQea7Q5mJIkmZ2UOBvkbtJ7IDR2YZUr2lWH0R/ZFxgD1fhqnR9PPVlXMyLpZtFXUf72Nuw8qrCGbGqlcdrPVhZPVdGlZhXL7t37wNR45mTpAB0yS216btgY97nUe6BDmkFEVuAvCTXOiGZ5rxoxCpcVNXvTQJBzXwXxG6BTV/QsSsWaHfVKQk2zM6G3esWXDhPuucNjYRVXTELfoFYIKU0ZjCyVptCcnE/7rZKXlllQGLeMnfy4pdZS8ecu787B9Xzl+O6Sj43NVVlUjx2byVD9XFzml4FZasiX1SqKpilHYIt/+ec3UfghBY1VnZHtpOcYtcFbYzYpjjmqW0udV42RMnKE5wp2YefYedFZGIJlidHZMW+UsJPbNn9bVZN1mK/AnbNtnnrNiHKT2JVFaw8lKyy0pKQTDR4c+OQX3/W+Ro88+/MUS0+1WtZcS3Wojeh1eK5eCkXLdTuuB+D/T1jg6RqqKKdHk7aCHeTpOyFyF4KHJKGy4aYjjzNnu7Nq1ctDMxQ1Rc18imFGTddsoL/VF8DsWYazGq/Dh7A0pn0sFqdzZqk0ihyqqYmov6Vm7Y7+lYxyItekv35J68cSdNuZXp0rP6ym8PvciCDJXIXgYZWRlOhjx9FR8FMSgoYnjpMSFk85PTi1BQbceBYD6WvkML1ZtYbUnBx3LrsduxXna01Nx8b+mW76/aXdtNvPl/7Z2hNvxJ102vYpnPsyF6VSQuQrB5xXCoYoqLNu6D396b56v5wUiuQVvz1iHIgclhvWmtEWbzJf6HQsS+xkA7wfM1dv3o3D0uJi6Q0aiIjRtoD4Kxvj1zGrrBIFqp3IqJIqAsyuop0KWj3Q5Aiqi3vQiHaqowux1u0z3a9fMvNHOybqe2GZ5NSN6tQ7M9wZksELwMiTy3ZnxjsJdB1KLX06Ft2eswx0ObNX3fLIwZhY16h/mS32nvYRVX+JfHXdkzOsPZ0VWW3940653g5cPV+yxnYbsqmCwIcnv7lG14bWpXHfVOjwV279q5dSuWX7Kx9CXQr/r44U4//lprj6vXzQf1SJengmLtqLjnV8kLV+qBKYQiKg9EX1DREuIaBER3eTn+VPNyLXDrJaRVSlcFXRIEBmUn+vMjv769LWO69o7QfUD3bJhbCGyn7cntisfqkj9d77+VPOSGUFaZjbuji0OpzdTqowySpUVCXIF7FA9aVNhMjrv+Wk1yW1WPqP1Ow+kPMvf77D/g2qCXCFUAvgzM/cAcDyA64mop18n9ztpysvM6ETscFHZNVHiDOB8Bqrall1aVhUTm/7FAvNMVT0Pfh7px5xKrf+CRvVNtxuvQ6JuV06aMalAlcno39/bZ+06IZVcFNWPjKos/L0H7TPuT3r8G2xKsUSFlxNWOwJTCMy8mZlna3/vA7AEQHyVKI+o8FghGBVOUD8wEFvnJhEqVwgz15jbV5PllR9Wu2qFCdSWndjroqeyEatLYjSHFDQ2Vxx7DlbgnRnrHEdJpYqTFcLuA+Wmv6N+U1SZpkIqkwLVK0wnkx0nfLkw8UQkVYLK8g+FD4GICgEcDSAulICIriWiYiIqLilJvjn1nR/H2tCNnb4AtYlUT02KzaS0yl5eumVvyqWDVc7ErfIz9HxvUwtHT58jDktVnDiWbHber0GvhCuqOOnS3lZXxHjZzWLI563fjX73T/Qk38AKJ7dD/wcmYezUn+O297t/otKeFW5vzTnrdtUozisHFSqTAzB/5pPBrHubE4Y++W1MWXU7f3vGKgQiagTgQwA3M3Pc087MY5m5iJmLCgoK4g/gkLd+inX0mg18+0165ybLx7NjKyCaVcwsLavEyGe+c5U6b4bKW8fJQ5OohkyURBEmXtPlri9jXl8y1jrE0Y4HP1+M+z+Lj1P/x9crYl6bKZwPZvkfZuzUkfvol/HF3wBn4cl2TFu1vWZQdzuunff8tBrTmtPqtk5RtUo/srl753Th6HGu+nr7GaCgJ1CFQES5iCiDN5n5Iz/PbdTAuw+Uo89fJyb83PbSMvz3xzUJ9zPaEM1m8dstMmCDxG6FsGXPIRyqqIppk2lHmOLhU+XVH9bEbTOaDtoeFhtq+Kd359pmx3qFm0E4WqpcT6qK/Fcv/VQzqLtx5BonI6pNRuWKTEaqJu9rtlvfGz86zMpXTZBRRgTg3wCWMPNTfp9ff/MdqqjCYoemiKKHvsK9nzrNaow8cP3un2g6Z7NqNuIWlfNwuxXC8Y9+jYvHTndcH8grdVA4ehzmWMR/B8k+Q2TIR3M2Ysry5M2cbnhtWq0D2EwRD+veyvGxkinBYcVkF90DjSY31fMJVSaj94rVNGAqtYkkcmK69YIgVwiDAfwfgKFENFf7d6ZfJ9df8O73jMf3Dgp4JcuegxVxpZG/WboNI5/5Lm7fbXsP4fP5m1wd/4DC9oVWrSOjERrz1u/GOf2c+f69TEw7z2X8tx+o7IPhlvU7a8NQzcaSXm2dF1xc7SCcFzBfXRi3lbmI7DHKrfruUWUyqqpmZTZ+s2sI1CqvbXsPxeQ+eE1g1U6Z+Xuondy6wpipvGaH2kJsRrbujTUhPffNStP9Lh47Hau378fwHq2RZ9OfWI9V8TWV6GdvOQ5KSm/ecxDPf7Mq4X6CWt6esQ7fLouflbsZvlROTpvlu8gQ1533m2XbcP1bdgmH7lEVWVhZzXjUopCfKqKTvIGPfA0gUqHZDzKq/LUe45LMa6/+fEPLvxKTQXzCoi01FUL3l1U6Vgh+ENsrOfG1MmswLniPVUa6m8WayoVdg3rOhxi9M/zqV2eqE0JjuYLy1wBQVV2NlxXkaNihsnOhGwKPMgqCdTsO4OlJy2O2TVjkbYy4Matx/a54h9LvXq9tzL3Lp3aTiYjao/XL1rrjKs4c3JSQcNODW0+0kKAeN2bDTbu97Te8X5Fp1Q/7/qHKqkCCMjJSIQx54pukwrpueVddv9wOLezLTaz12ITllOg9qTcZ+WGiEtTiZmzZ5aIHtx4zf4GbQe3Mf8T71FRSaFI7KBmuOP4oJcexIzc7K5DKqhmpEJLl4zkbE+/kgFd/WJ2wVaWTLlV+EJ3h6U1GN2rtNoX0wc3YoqKzWJSwlIuuqmZlGcb//NZ731jbwxrgsfHmeSJekrE+hCC5/7PEJQG6W/Tx9ZvjH52ccqKSEDwHy6viEumsKF6TuNm8U8LSMyLZLHUzzPx/qrn9w/men8MMRwqBiLoCuA3AUfrPMPNQj+RSil28r5eMscgEdUJpWSV2hGAgThdl4DZUt67xi/72hfVec9giEgCG92idojS1FCYwjfrFV0v8qSOV7jhdIbwP4EUALwEItulnEpQF1Kf0xSnJLy0fGrcED43zNrStLvHHtzLbjNWltboVZfNG9RLv5ID/zdvkOOHTa5z2TM50nCqESmZ+wVNJPCSorD9B8IuyympMXqpoFsxqyo6EydckQ4AznCqEz4joDwA+BlBjQ2BmdcZGDwmy9LQg+IFT/4ATGME2/1FNXaqp5TVOo4yuRMSHMA3ALO2fP90+FJDvIjlGEDKd6mrGhl0HE++YJqiKDgySm96Zg8WbvDe/ORopmbmD14J4iVmdekEQzJm5dpfnmbh+8qf35gUtQsp8OncTzh/QzvPz2CoEIhrKzJOJ6Hyz9/0uWZ0sZrVdwoYsa4WwEGSRPsGabB/6iyRaIZwMYDKAs03eYwBpoRCMhezCiOgDQRDsyPIhjdhWITDzfdr/V3svinekg1O5451fBC2CIAghZsnmfRjUqaWn53DsbSWiUQB6AciLbmPmB7wQSjV1yUEmCEJmstSHnA5HixAiehHAxQBuQKSHwUWIZC2nBW6yNAVBEMJIrkk1WdU4PcMgZr4CwC5mvh/ACQDaeyeWIAiCoKdedngUQrRQ+QEiagugAkDKoahENJKIlhHRSiIanerxBEEQ6iodWnpfF8qpQviMiJoCeALAbABrALydyomJKBvAPwGcAaAngEuJqGcqxxQEQairnNy1wPNzJHQqE1EWgK+ZeTeAD4nocwB5zLwnwUcTMRDASmb+WTvPOwDOBZC4NrQgCEKGsb/c+6rNCVcIzFwN4End6zIFygAAjgCwXvd6g7ZNEARBMHBAUQtQO5yajCYS0QVESlPlzI4Vl55FRNcSUTERFZeUlCg8vSAIQvqQ5X2ismOF8CdEeiINywdcAAAgAElEQVSUEdFeItpHRKkGxW5AbKRSOwBxXU6YeSwzFzFzUUFBcja0gYXNk5NQEAQhJLRvpqYntB2OFAIzN2bmLGaux8xNtNdNUjz3TABdiKgDEdUDcAmA/6V4TFNeufpYLw4rCILgG80aqmlcZIfTxLSvnWxzAzNXAvgjgAkAlgB4j5kXpXJMKxolaGgvCIIQdrKCLm5HRHkA8gG0JKJmqLX7NwFg38TVAcz8BQAp4iMIgpAAP3wIiabOvwNwMyKD/yzd9n2I5BAIgiAIPqA2psecRCajaQAGAbiVmTsCuB/AQgBTALzlsWyCIAiCjyRSCP8CUMbMzxLREACPAvgPgD0AxnotnEr+eGrnoEUQBKGO8OpVdTNQJZFCyGbmndrfFwMYy8wfMvM9ANJqhPVhtSUIQgYwbfRQZPlh0A+AhAqBiKJ+hmGIdE+LklahO3Xz5xMEwU/q5WShbdMGdXY8SaQQ3gYwhYg+BXAQwHcAQESdETEbpQ+yRBA85PyjpepKJlBeGem+6PdwctEx7Xw5T6IWmg9r+QZtAEzk2k7wWYg0y0kbRB0IXrJl76HEOwl1Bj9yAvQ0yvPHIJPwLMw83WTbcm/E8Q7pYS94ybRVO4IWQfARvyeYP/p0f3nfgicscLhVwhUnmHck/e72U32WJD1p7NMMKhN4/ZqBQYsQevzICdCzevt+X86TMQqhOtz6AE3yctGvfdO47et3HghAmvTj5uFdgxahznBSF+8bsaQ7fvsQcn1onwlkkEKo1GmE1Y+eGaAk5pzavQA5JqFsew9VBCBN+sEhWAFOue2UmNedWzXy9fxm948AvHHNcSi+e7jSY/rtQ+jVNtVaos7IGIWwYOPumr/9Xu45oayyGrnZ8XL5feOlK4s2pVqNPXWOahHb89avWV2UyhAug1+5qihoEXBil5Zo2ai+0mMu3eLv/fYHnxJrM8bw+sNKa6fMjcO64B9fr/BRmngOllchJ6t2AJl0yxC0PiwPu/fLCsEJZsrUT87ofXjctiWbg1dSQTO0e+ugRVBK1Fe1cKN/UfcvXDbAl37KQAatEOxo2cj7OuOJqGZg5pqdNa+7tG6MJnm5aFg/O0Cp0occn2fjRr5cuCXQ89c11owZFbQIpkTNcvqV+1WDCj095xl92nh6fD2iEADsPej/LNx4w5/YuSXKtKQXPWExb4W9FlS9gBWCkFnon8u/ntMrQEnUIk8RgL9NDDatYs2YUWhQz3wlEBY/4aBOLWJet26i1iabKiX7yoIWIXCm3HYKrh5cqORYZ/dz3+7krL7+zWST4f+ONw/tTgZV8w8r38adZ3ZXcwKXiEIIEWazXFKYAtOxZUPMuGuY5fvP/epo0+0L7x+BQZ1b1rzufnhjbN0brgG4lccKqll+rqfHV8FRLRri3rN6KjnWExf2xQkdWyTeUcdzvxqg5Nxe8eAveqNJivkq0ZVBdhIrdzNl3dMieujaIZ0C8YtljEIwc/qFjRyzG0DhPZGXm41WjfMs3z9RN+jrMbYgDTry6erBhVj64MiYbfVzYldYLRvVxx9O6aTsnPee7WygDbpdqyoTY15udlKrhExheE/3zvKjj2wWt23nfuuJVfFdpykPl01EIAqBiJ4goqVENJ+IPiai+IwsxTRNMMNLJjehcV4OXr06cV307oc3dnQ8szBFlSajxQmiXpwOJh0KGibeyUPembE+7lrNWB0bRfbZDYMxrIe6CJczerdBnyMOs3xfituFh6cv7odjCyOD74heaqOcovkudhMrK/Jy4p/vhRutn8nD8nOVh8smIqgVwiQAvZm5L4DlAO7w+oRVCWK0k5lZFTSubzlb7ta6Vgncd3as06lDS/MB9ZdF8RUN/XQqOz1Vu6YNvBUkAQcrqpBt0JTGGHwCKS1nkZ1F2FFqPZsLU338NWNG4aQu5qu9TGDX/gqUllUBiPeF7D1UmdqxDyQfgJKXG/6IwUAUAjNPZOboLzMdgOe1Xb3I2cmmWgt/w3rZGH1GrSNowi1DMLBDcwDxA+1x2nYjZjdMomGmnsmsI1mcDmnhS3+Kd85lkdrEsGwiDO3RKuF+YciYDgv3KPJnuGVVSWlNDsg2j3xdnMRTIArBGb8G8KXVm0R0LREVE1FxSUlJ0iepVqQR9OGiOdlZNYP9/vIqXHdyrM364V/0xtDurdDfUKPIasww257IXn9uAHbe6mrGkc3zfT+vHcOMgzUpVpaEmMTBuPcTfL598wa2JifVVCtQTMkMeqO0mPm/nt0T15zYIeY9p6bTVNEHPKi4DonO4ZS83DAMt/Z4JiERfUVEC03+navb5y4AlQDetDoOM49l5iJmLiooSD5br8rBjfHMxf1dHZOZbaOAurRujFeuOjZuZnBqd/PvYTaAJTLjqBz07AY8PWGcA+cZnMoEQn1F12ZEr9YgojgzlRn7y6tq/v74D4Nq/v7u9qG4YIB3foaHz+sd8zqohcpNw7ugUf0cjDAJ4ujo0Pf06fWDE+7z/GXWEU1WkTuAumixrUn0vzAGPoQRzxQCMw9n5t4m/z4FACK6EsBZAC5jH9bZiXwIAPCLBI7Bt35zXMzrpVv2ua56+PEfBmFkb/N47UIT30Ki49sVNHPyYLk5119GRkxiqcS6XzqwfdKfjWIWOmhUjFmkrtjbv/4vUo9niIvyAY9f0Nc0qsQrLjtOXYx9lLKK+ETJRHRt3RgL7x+BNofF+5mydROOFy4bEBcpBkRW4GZVf/U8fF5vS7MrAHTSKR7jyBJkomc6rBACiZEjopEA/gLgZGb2pb6z06Xj2f3a4rN5m+K2f3nTSejRJn7m4fb20h9jzj2n4WBF7YzS7FiJ8hDsnJmtm8RGQow5v4/tsYzPChEw/Y7avIXfn9IJv9dCObfbOFjtUGHXN6sialwNEJFyR+8RTZ1HlhyqrIrblopD0i1mt/sVJxyF//641vExxoxfqlCiWAWdXz8nKZv66kfPBBFh5/5yy330K7lZa3fht67PkphkQq/9LnaYDEFJ+ByAxgAmEdFcInrR6xM6WSEAQG+L5abl7+/yvtAfp1nDemibIGIn0X1nlyBDBAzt3godWzbE6kfPxCUDj7Q9lvEmP7tv2zilEuVAefyAZ8SrSJcLtP6yen9O47xYUwAhueQhO+xml1sNmdLRa3nL8K645NjIqujvPhZQNLP/n9W3LR78RW+Tvc0pNymlkgrZWVQz+ze7kkafgxnR38BO1+tNn5OXbXMlo1OiK5QuFiXOP7/hxLhtzRs6q5kWpH8ukBUCM/teGGfzHmc2v2tO7IBHv4yfGVnN1N1mEtvNLMzCURPNROzs2gTglasS50lYnevoI1NLD3n9muNw/ZuzMW7B5hiZUuGW4V3xy6J4s1OZYUaeRaR8RmYn+9TlsQEP0Ut50/AuSmVIheYNc7HLZmZtpM8Rh2GBwqqeudmELXsOAogvNTKseyvccYbzcg12ylm/EvGqR0T75vl47ldHmyZzXnbckehtEkDgdEUUZBJt+Ncwipi/wdmNnZOdZVpCoqIqfrbUoWVD15EYdrdn19aN0b99U7xz7fGO9k+4g8tnwahcjrcpXRBULaObhneJGegP11YwcYpRcZQR4M7+XKF4du0WM5NRh5aN4hTn6TYZtypr/wARJR2NzvlozgYAwG0jugEA7hzVw1XFWrufIlf3uxsnOSrVw1l926Jpfvys/5ij4n1Hbmb9bnxVqskYheCGZg3jIxGKdaWpo1w6sD0+mbPR1bHtZvT1crLwyfWDYwbihGOQjT5KtQ7Spt0HrU+bQA+e2s2fm3r6ncOwZsyoOCXuhe/QzSGzA7YXm/082VkUE+ny25M64B+XHo1+7SzCYRVfQ/1sPTpQ/+GUTlj64Eh0Kog3vawZMwoXHmOeomS3ctafp7QscSLaVYMKcePQzkpMnO2aNcDQ7rEh0F1aNaoJ8Jh0yxDbz//zVwMw2KKEjB9kTIMcN5jdbGazwywi7D3oLvPRbZRDov1N6x/VfNbVqeKwc4La6YOJtwyxzMZev8tayaRC/Vxj2Kl6zFaJVhS2iJ8R1svJUmqXv+vMHnj4iyWm71kF7ukjXe4aFUkcm2exelZ9DXNMVnZEZGtK+d/c+ACPRNiaUU3eunVENzSqn4P1Ow/gpMe/cX0+PVNuOzXu/JP+dHLN311a2+dijAq4YqysEEwwqx+i/42jN1XvIw5zFJvuJRtsBthUJbOL47crN92+Wb6l/X7yUm+cfMaicvrorWQwc74u3bLP8efNFGJHw7YnLuyLmXclX7zspK7OZ5KvXzMQAFDYIl4uKzv7HsV9QnKyqKaky6XH2Qc4RCnXlPARTRvEBBHY3dtuq4RG9zaa05LB7Xjw90v6h6qUfMYpBDO7svEGuszkZtXP1Jc+OBIfXHcCju/YwhPThBs+tZlBpRpznWxil763QzLZrsmQqKG92w5cXU2O5+ZhN1tlGrf1a980pUQpO4epcYFwUpeICa+FyWTHKvJoVcn+pGUzg4hqnhe3zt6NBvOl3a2tz3coirPnx38wuvJbtqXUlUwqOKpFQ3z4+0GJd/SJjFMIZst+Y1LPtFXx/Zf1N2D9nGwUFUbCzo5wWOjtg+tOwN2jeriQ1Dn6XAE9bh65FQ+fEbct6NVPKjjNurb8vMksM1WFcJTBjJSTRSmVEs9O4jv2aBNvsrBSStGQykTK1il7D1XUFCF0+70PN4Q/2/nH9MrGSWRPVHkG0a62Y0HD0HRFBDJIIUSvefTHv/7U2rpDVvbumM9b3IBn9Y3UEmpo0fEsSlFhc/zmpI4OJHXP4YeZ5wq4uc/MTDxuTCRWVFYFU0NBn81sFvWRCLMyA256HZhd+2cvPbqmLDMQueapjAUtbHqBW1118/vYXIjoavrYQuusYDcQgJXbIrPwQw5NetG8hQuOia0iYBtlpLuX481H8VcmqpzsSl54hb5AZhjIGIWQa5hN3Taie00hLrsHK8qAo8xj8qM3ZlBNY64dYq1kdriIOTdDRZ/iiYu3pnwMp8y997Sav6MOzDVjRuGD605wdZzFD4xA+2bxTuHebVMrTpeTnYX3rxuke00pzQ4b2ygoK6dydOXzKwsbvn4VFD2GqoWiPjfUrH+4Gcdo5T+M197u3tSv7owTne2l8c9EdnY02c3/ZziLyPG18IOMUQhmq+vztNpF/doZqpGafN5plqHf2JmsUr3BEzUVMiOZ2bgqzGLCAfe+lPx6OcjXzAePnFdb7oN099Dse04zfiwGqx7ZelI1a9l9L6sVQm52Fub/9XQ8eG6t3yB6mEGdWmDVI7WNoqIDuKqBUj9bd5rwdtuIbrjv7J4Y0Ss2WcuuLIl+8pfrwA8WXfkFMaUjql01hYGMUQhmZQyG92yNNWNGob0haeS8o+NLSlst7KIPSxgrgLZ0sPIx8uRF/WoiQc7s4z4EzksH2UtXFHl2bCO52VlYM2ZUzExaPzBaTRDm3Xs6/vPrgWiSl1iZqliBWWFXqaVJXm7MSmC/Fqtv9J1V61YIKiJh9B3sSh02qmlQLxtXD+7gqi5VdjbhYi2b/WKTrHYrglghEPlXFtwJGaMQmmuDo5OQtKHdW+N9g5nB6n4MkT8ojmSKh11wTDtMuGUI1owZhYYB9wc2cloSfWxVYjcmnaP1pTgsPxcnO8w0tcshScQLNuWf3WKVXxH1HQzr0TpOwXxy/WDX7R31Hey87DCXm0UYc0EffHnTSXFZvy1sVvp+KYRXrqqd2GQlyMPwm4xRCP/9daR09fib7TMFoxhtyFbL86AVgl3l8KD8GnUVu+uZTLmBZBVC8d3DcUai1ZuLivJWq9/eRxyGnx85E0O6FqDSoDT6t2+Kh34RaQ17XIfmpqWsjehX6V4GsGVnRXwzZtWJba+KT4/L0O6ta8xUWbpQ3DCQMQqhQ8uGWDNmlGmKvBlxpaCt9gtVjEAsXj50F1mUFDAy0mD7VV0fx0+sFMLce0+zLLFgR7LVWFUr+qi/pKCxSUKmdhOZRYud2KUA/ds3xcPn9bad5UbDQPU+qeE9Ul/t3TjMvHCgXU0k+wlUyiI5pvju4Zhy2ykR5eXfaROSMQrBLcYfyeohDJN2N+JlfPMj5/fB1NtOTbifsSlImK9XIowDxt2jeuDUbgWWzuxEJJvnof/U4xf0xScmjZCS8WnpQ2KNVJo4JRrVz8En1w9G51b2NvCLNDt+k7xctG8eCYJobNLkyC1/Oq2rqQnYzixsV47lUBINgZIlLzcbR2lZ42HKQwiXkThMGH8jKx+C54KEk9zsLBzRzFlSXl3BOCn4zUkdU8otSXYg0Mvxy2PNnabJ9CC0m6umUgrk/nN64fcnd0KzhvVw87Cu+PP789CqsfNmQ3ZUmKxckla0ih9mp9nxYRpDRCFY4HQhF7Ryd9j3xxOMz51ZOQLjoDdv/W4vRfKUoH/rGjySw2mZkXvP6unquPVysnCklqV9wTHtahoceUWypriG9YIZDkNzX0EUgiVOf6TogBetYunEuaaSxZv3+no+PcbB3izbecW22Gxnq8qaqrjv7J6YsGiLJ8cOy9LeiRhOW8YC7v1gv3bQ2SxIkl0hOMkdqeuID8GCeB+C+X7rd0ZaQkerMublZvsaRrbFYSc4PzC7Rgs31iqsf1/pfR7B1YM74J1r3WUmu+HuUT0wwWGkmlNamThz7XDiVE7GZFRXCIvidkqY3MqBKgQiupWImIiC6whhgfGmsrrJVJcIdsv3K7enfIxhhoYeyWL3ILZqXD8mMSld+c1JHdFNcSLRDJclsJ0MH270QXlVxD+wzyZZLFr5tokCZ3CYGdTJukugZ4RHHwSnEIioPYDTAKwLSgY7nP5G6TYbMfL4hX3xr/87RsmxBne2fpiM/oWjWuTj7d8eb7F35nHVoELH+1Ypnv6/O3M9AOC7FdaTi2hF4N96VKDRLyb/+WTb998K4J4M0xAS5ArhaQC3I5xVH+KEsjIZpXGFaACRED03vWytuH1kNzx/mbViOaFT7CJwf1lV2l87ldxzVk8svH+Eo30drRBcKI1qB9GW0RpVZs3j04mOBY1qssrDQpgeg0AUAhGdA2AjM89zsO+1RFRMRMUlJSU+SBdh8x5nrR79sP+9evWxnh17o6KWln84pbOtM89Ya357aZmr8gXXKHBk2q1g/OSGoZ3jSq5nZ5Gr8tpKcfAzjOrbBlNuOwWnKjIvquL4ju5Lc5tlMKdKMnXDooTJyuDZHUhEXwE43OStuwDcCeB0J8dh5rEAxgJAUVGRb6sJY5KKVWXKfu29nzGd2s39Q+j0QTGL4/YCM+e308fgz6d1xQ0WWaluGNq9NX5YGd/8yG/+fHo3/Pn0bkl/PtUqqUYWb3IWqXaUSfvNoBl7RRFe/3EtnpiwzPFnuuga/qhqXzn19lOTfpbCow48VAjMbOopI6I+ADoAmKdpxnYAZhPRQGb2Jl4wCYyZlGatNwF3TVO84KQu5v54p5E2lU7sBQooXrsrbpvz0F41Mvx6cCFG9j4cg8dMVnPAgHASHukm7LQ8RPX43dIkL9f1KqGR7tn+6U53Dn2r7nL5KeQwhGiB4L/JiJkXMHMrZi5k5kIAGwAMCJMyAIB2DrNwo8s9u2YlXlKVYmbazv3+REmZ3fMbHJqrVC2pichxy9N0x43fWVWLzOCovT8ePb+PzX4Rovky/dubN71y8lmVSNhpGqD/4R+yaEKuZ1+Zs/ruyTLQoo1hqgphzrr4mbsXGAf17CzCLJNVgxlhKg+cLri5KzoWhM8U5IaoK4oIuHSgeSc4Pb3aNkGvtk1w79nuMq4BbxRCuUX58SAIXCFoK4XUg+kVow+TvNymQueBcm8VQZTXfn0sbj29a9x2s6JjbvCrE1zUtBWtdnrXmT0cV+0c3iNcjsy6xulaRdoz+5i5/MJPdLLR9jBnq7+83GyMu/EkDDjSfXe/m4en7ssy4iYizGsCVwhhxamZwosZgxn59XLQqkl8QbBUFYLXPpBoI5fogxS9rNlZ5LjEgPR1qOVTk8qmZrgZZPK1FVhY28Qmwq+7o12zBjWVW1USpkZU4ZEkTYk6n3/R3/vYZrN69dUpKoQmDdz3TXbDGX3axFR97KmF/HUqaOTYZLRsy764NqeZSj+Hdm83t8WwHq1w39k9cbFF5dSw49d8odQjs7BZUcigkBVCitTPycaMO4fh8Qv7eX6uU7oW4Ky+sZ2yzNofGmPcjejtrMmEtKbCxce2x8RbhuDELi1ryiEI6lm9fb/jfYkIVw/ukFKkTJD45ZTdbdNLIRXClIcgT6QCWjXJswxLVQkR4blfxfbSNTMZ9UmQTfro+X3wxY0noX3zBjixs79lpIgIXVtHagE5rW0UoudFCCFEsf97RdvD1PRwCDPpOSUQajizd7wj0MnStmfbJvju9qFKZDAzZTnB7zwEQUiFFo3UJLGFGVEICbjihPD2AP7zaV3xu5M7xW2fvHSbbzJMv2NYTU9et+Q7rD+v2iTw8yNnouOdXyg9phAcfq0Q3JRaSVdEIdjgtAVeUKgo55Aqh6ewjE7k64iiuqCa3YOdTG0cIVj88iGEyfnrFeJDEAKjXbN8vHHNcQn389Nk5DSWPex0V9yzIcyovj+swm+T7cSWTohCEALlRItaTHrc1OURIvRr574sQ7qiOk9l0i1D8PiFfQEAtwyvTQZNtldzOiEmozTkhI4tPIuJDhvNG9ZD0wb+JExdeEw73HFGd1/OJahD9TjdolF9/LKoPX6pJaE9/dVyAEBOtigEIYS8fW3mdBqbfc9pvp3rgXN7pW0sfiZDNf97O2CLyUgQfOCpX3qf1OeEMFWdNMNNpVYOZyNCT4iuEFIt9JiITDAZiUIQAuf8Ae2U9XV2ilmSUdif929uPSVoEULJ4s37AAAbd6vp/mdFJoSdikIQQoHfBey++8tQTLxliK/nTIV6OVm+ZMOnIztKy3w5TyaEnYrBVAgFPhWN1Z2PUM/vk/pI2M1fKunW2p8Q2z0HvWsmNe++00NR2yt4CQQBwRT4MisMGCZaNqqPG0OQfBh2jtISHA83KQ+vkn2HvIvsO6xBbigaQYlCEEJBYQAN3Bdvjm0uHwYfwpx7TsPfL+kPAOjcqiEuPy5xB7BMx6+fLdXeI+mAKIQ6yHlHHxG0CK4JojlLr7ZNYl6HwczSrGG9mmKBzPBvtEtjoorc68iqsK8oVRCYD4GIbgDwRwCVAMYx8+1ByVKXmDZ6KNqkYZneIGbnRkd2WHyGesWUrJLKqLBT7Rp5ndDuJuw3XQlkhUBEpwI4F0BfZu4F4G9ByFEXadu0QagabjglCIk7tGwYk5mcE0Inc42SypzxPbR0y4D6UEE9Ab8HMIaZywCAmf2r1yyEkiD6JhORafnwMJGscg+D+csvokozlcq7TnDTpzpdCUohdAVwEhH9RERTiOhYqx2J6FoiKiai4pKSEh9FFPxEYuxrOVAeiWb5afXOmmG93KX9OpMKArZqkocx5/fBy1cUeXoeMRmlABF9RUQLTf6di4jvohmA4wHcBuA9spgKMfNYZi5i5qKCggKvxBUCJjc7C0seGBm0GKFAn3Gb7MLp/VkbFEmTHlwy8Ei08ijstF/7SOXYPhlQQdYzpzIzD7d6j4h+D+AjjqzBZhBRNYCWAGQJkME0cNhBra6jT1DKJNNPWDmj9+GYt3635yapMBBUlNEnAIYC+JaIugKoB2B7QLIIGc6PdwzFpt2HghajhhYNI2GnAwubS9hpCPjdkI645Nj2aJrvf2i03wSlEF4B8AoRLQRQDuBKzgSPjRBK2hzWAG1C1Ckt+iA0aZBTYzKqy2U2wg4RZYQyAAJSCMxcDuDyIM4thJtjC5th5ppdQYsRKFGHMBHVxNaHoc6NUPeRu0wIFV45BtOJ6GI5i1CzXEjD1BIhDRGFIISKR8/vg8cu6BO0GIFySrdWOK1na9w9qmfMakEQvEYUghAqmuTl4uJjM7ugW15uNl66ogjtm+fXRF5dMrB9wFIJmYD0QxCEEJOXm41lD40Up7LgC6IQhFAyqk8bnN6rddBihIL6OanlZzx/2QBFkgh1HVEIdYTvbj8VW/aGJ5Y+Vf4pg5gSRvVpgzP7tAlaDCFNkHVoHaF983wcW9g8aDGEkDCka6TMy8ndpNyL4BxRCIJQB2mtNdlpEIK2jEL6ICYjQaiD3DWqB1o0qo8zeh8etChCGiEKQRDqIE3z62G0rvmPIDhBTEaCIAgCAFEIgiAIgoYoBEEQBAGAKARBEARBQxSCIAiCAEAUgiAIgqAhCkEQBEEAIApBEARB0KB0amVMRCUA1gYthwtaAtgetBBJks6yA+ktfzrLDqS3/OksO2At/1HMnLCwVVophHSDiIqZuShoOZIhnWUH0lv+dJYdSG/501l2IHX5xWQkCIIgABCFIAiCIGiIQvCWsUELkALpLDuQ3vKns+xAesufzrIDKcovPgRBEAQBgKwQBEEQBA1RCIIgCAIAUQjKIKI1RLSAiOYSUbG2rTkRTSKiFdr/zYKW0woL+f9KRBu1bXOJ6Myg5TSDiJoS0QdEtJSIlhDRCWl27c3kD/21J6JuOvnmEtFeIro5Xa69jfyhv/YAQES3ENEiIlpIRG8TUR4RdSCin7Rr/y4R1XN1TPEhqIGI1gAoYubtum2PA9jJzGOIaDSAZsz8l6BktMNC/r8CKGXmvwUllxOI6D8AvmPml7UHIB/AnUifa28m/81Ig2sfhYiyAWwEcByA65Em1z6KQf6rEfJrT0RHAPgeQE9mPkhE7wH4AsCZAD5i5neI6EUA85j5BafHlRWCt5wL4D/a3/8B8IsAZamTEFETAEMA/BsAmLmcmXcjTa69jfzpxjAAq5h5LdLk2hvQy58u5ABoQEQ5iEwiNgMYCuAD7X3X114UgjoYwEQimkVE12rbWjPzZgDQ/m8VmHSJMZMfAP5IRPOJ6JWQLv07AigB8CoRzSGil4moIdLn2rSXb54AAAcfSURBVFvJD4T/2uu5BMDb2t/pcu316OUHQn7tmXkjgL8BWIeIItgDYBaA3cxcqe22AcARbo4rCkEdg5l5AIAzAFxPREOCFsglZvK/AKATgP6I3HRPBiifFTkABgB4gZmPBrAfwOhgRXKFlfzpcO0BAJqZ6xwA7wctSzKYyB/6a68pqXMBdADQFkBDRJ5dI658AqIQFMHMm7T/twH4GMBAAFuJqA0AaP9vC05Ce8zkZ+atzFzFzNUAXkLkO4WNDQA2MPNP2usPEBlg0+Xam8qfJtc+yhkAZjPzVu11ulz7KDHyp8m1Hw5gNTOXMHMFgI8ADALQVDMhAUA7AJvcHFQUggKIqCERNY7+DeB0AAsB/A/AldpuVwL4NBgJ7bGSP/pQa5yHyHcKFcy8BcB6IuqmbRoGYDHS5NpbyZ8O117HpYg1t6TFtdcRI3+aXPt1AI4nonwiItTe998AuFDbx/W1lygjBRBRR0Rm1UDEBPAWMz9MRC0AvAfgSER+wIuYeWdAYlpiI//riCybGcAaAL+L2obDBBH1B/AygHoAfkYkSiQLaXDtAUv5/4H0uPb5ANYD6MjMe7RtaXHfA5byp8t9fz+AiwFUApgD4DeI+AzeAdBc23Y5M5c5PqYoBEEQBAEQk5EgCIKgIQpBEARBACAKQRAEQdAQhSAIgiAAEIUgCEIdhYgu0oq/VRORZZ9hLRt5GxEtNGx/Qis4OJ+IPiaiprr3+hLRj9rxF2iF5fKJaJz2mUVENEa3/1VEVKIrmPcb3XuPaQXqFhLRxSbyPUtEpQ6+bzNNzvlENIOIeju5TnpEIQgZCRG10D2cWyi2uuU0j855NBG9bPN+ARGN9+LcdR0iOoWIXjNsXgjgfABTE3z8NQAjTbZPAtCbmfsCWA7gDu1cOQDeAHAdM/cCcAqACu0zf2Pm7gCOBjCYiPTZw+8yc3/t38vasUYhkkjZH5HCerdp9a2i36sIQFM4404AczV5rwDwd4efq0EUgpCRMPOO6MMJ4EUAT+se1kEenfZOAM/ayFQCYDMRDfbo/BkFMy9h5mUO9psKIC5Pgpkn6uoCTUck8xeIJG7OZ+Z52n47tMzmA8z8jbatHMBs3Wes6AlgCjNXMvN+APOgKSetAusTAG7Xf0CbOHxIRDO1f4N1x/paO/9SAIVE1DrR99cjCkEQDESX59qscwoRvUdEy4loDBFdpi3HFxBRJ20/qwdUf8zGAPpGBxEiOlm3IpkTzRQH8AmAy3z6qoJzfg3gS+3vrgCYiCYQ0Wwiut24s2ZeOhvaAK1xgWbO+YCI2mvb5gE4QzM3tQRwKoDoe38E8D+TpLi/IzKBORbABYgkNUaPdb52/oEAjkJihRRDTuJdBCGj6QegByIzyJ8BvMzMA4noJgA3INK3IPqAfk9ERwKYoH1GTxFiSyDcCuB6Zv6BiBoBOKRtLwbwkGffpo5BRD8BqA+gEYDmRDRXe+svzDxB0TnuQiQb+E1tUw6AEwEcC+AAgK+JaBYzf63tn4NIKYx/MPPP2mc+A/A2M5cR0XWIlKYeyswTiehYANMQqXr7I4BKImoL4CJEzFFGhgPoGalYAQBook0oxgD4u3YNFiCSqVxp8nlLRCEIgj0zozM0IloFYKK2fQEisznA4gFl5n2647RB5IGP8gOAp4joTUQammzQtm9DpHql4ABmPg6IrOYAXMXMV6k8PhFdCeAsAMO4tqzDBkTMPNu1fb5AxA8QXQ2MBbCCmZ/RyblDd9iXADyme+9hAA9rx3oLwApEfBCdAazU7qt8IlrJzJ0RseycwMwHTUS+WjsOAVit/XOMmIwEwR59HZhq3etq1E6oog9o1AdxhEEZAMBBAHnRF8w8BpHaMw0ATCei7tpbedq+QsAQ0UgAfwFwDjMf0L01AUBfzcyTA+BkRArLgYgeAnAYIitH/bH0BfPOAbBE255NkdpPIKK+APoCmMjM45j5cGYuZOZCAAc0ZQBEJiV/1B27v/Z/U6ptmfkbAFOZea+b7ywKQRBSx/QBNbAEkRlfdJ9OzLyAmR9DxEwUVQhdEc7qmmkHEZ1HRBsAnABgHBFN0La31Wb10f3eRsRU042INhDRNdpbzwFoDGCS5ut5EQCYeReApwDMBDAXkdLZ44ioHYC7EHHuzjaEl95IkVDUeQBuBHCVtj0XwHdEtBiRlcXlOke2FTcCKNL8EYsBXKdt7wFgEREtRaSk900uL5kUtxMEMvSOJqJSZm6kmSFuZeaztO3faq+L9e9pzsB/IvJA5iAyM7vO5DwLAAxi5n1E9CwiJqcqRGaXV2n25VsBlDGzZTSSIHiFKARB8AkiugXAvmgMusU+UwGcq81CBcFXxGQkCP7xAmJ9EjEQUQGAp0QZCEEhKwRBEAQBgKwQBEEQBA1RCIIgCAIAUQiCIAiChigEQRAEAYAoBEEQBEHj/wERlNoEkdiN4QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "m = catalog.Merger(\"GW150914\")\n",
    "\n",
    "# Get the time series data around GW150914 from Hanford\n",
    "# 'ts_han' is a pycbc.types.TimeSeries object which contains\n",
    "# gravitational-wave strain in this instance and has metadata\n",
    "# such as the start time, and sample rate.\n",
    "ts_han = m.strain('H1')\n",
    "\n",
    "# And now livingston\n",
    "ts_liv = m.strain('L1')\n",
    "\n",
    "# We can see how much data was returned and its boundaries\n",
    "# Note: All times are given in seconds since the GPS time epoch\n",
    "print(\"Duration: {}s Start: {} End: {}\".format(ts_han.duration, \n",
    "                                              int(ts_han.start_time),\n",
    "                                              int(ts_han.end_time)))\n",
    "\n",
    "# We can directly plot the time series as follows\n",
    "pylab.plot(ts_han.sample_times, ts_han)\n",
    "pylab.ylabel('Strain')\n",
    "pylab.xlabel('Time (s)')\n",
    "pylab.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Getting Data from S5 / S6 / O1\n",
    "\n",
    "In this section we show how to read data from the bulk data release by LIGO. This currently covers the periods of the S5, S6, and O1 analyses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pycbc.frame import query_and_read_frame\n",
    "\n",
    "# Retrieve the approximate time of the merger\n",
    "m = catalog.Merger(\"GW150914\")\n",
    "start = m.time - 32\n",
    "end = m.time + 32\n",
    "\n",
    "# Get 64 seconds of data roughly around GW150914\n",
    "# The start / end time may be any in the publicly available data sets.\n",
    "ts = query_and_read_frame('LOSC', 'H1:LOSC-STRAIN', start, end)\n",
    "\n",
    "# If we wanted to retreive data from the Livingston detector\n",
    "# we'd use the following command instead\n",
    "# ts = query_and_read_frame('LOSC', 'L1:LOSC-STRAIN', start, end)\n",
    "\n",
    "print(\"Returned {}s of data at {}Hz\".format(ts.duration, ts.sample_rate))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Directly Reading gravitational-wave Frame Files\n",
    "\n",
    "If you store LIGO data on your own computer then you can directly read in the data as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# We'll first download some data for this demonstration\n",
    "!curl -O -J -L https://losc.ligo.org/s/events/LVT151012/H-H1_LOSC_4_V2-1128678884-32.gwf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pycbc.frame import read_frame\n",
    "\n",
    "# Read the data directly from the Gravitational-Wave Frame (GWF) file.\n",
    "file_name = \"H-H1_LOSC_4_V2-1128678884-32.gwf\"\n",
    "\n",
    "# LOSC bulk data typically uses the same convention for internal channels names\n",
    "# Strain is typically IFO:LOSC-STRAIN, where IFO can be H1/L1/V1.\n",
    "channel_name = \"H1:LOSC-STRAIN\"\n",
    "\n",
    "start = 1128678884\n",
    "end = start + 32\n",
    "\n",
    "ts = read_frame(file_name, channel_name, start, end)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate your first gravitational waveform !\n",
    "\n",
    "Here we'll generate the gravitational waveform using one of the available waveform approximants. These can be generated as a time series using [`get_td_waveform`](http://pycbc.org/pycbc/latest/html/pycbc.waveform.html#pycbc.waveform.waveform.get_td_waveform). There are some additional examples using this interface [here](http://pycbc.org/pycbc/latest/html/waveform.html). The key parameters are the masses of the binary (given in solar masses), the time between samples (in seconds), the starting gravitational-wave frequency (Hz) and the name of the approximant we'd like to generate. A variety of approximants are available that include different physical effects. A full review of the different models is outside of the scope of this tutorial.\n",
    "\n",
    "In this example, we've chosen to use the 'SEOBNRv4_opt' model. There are many others available as well with different methodologies and which include different physical effects. This is an implementation of the model introduced [in this paper](https://arxiv.org/pdf/1611.03703.pdf). It models the gravitational waveform of inspiralling and merging black holes, and includes the ability for each black hole to spin in the same direction as the orbit (aligned spin)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pycbc.waveform import get_td_waveform\n",
    "\n",
    "# The output of this function are the \"plus\" and \"cross\" polarizations of the gravitational-wave signal \n",
    "# as viewed from the line of sight at a given source inclination (assumed face-on if not provided)\n",
    "hp, hc = get_td_waveform(approximant=\"SEOBNRv4_opt\",\n",
    "                         mass1=10,\n",
    "                         mass2=10,\n",
    "                         delta_t=1.0/4096,\n",
    "                         f_lower=30)\n",
    "\n",
    "pylab.plot(hp.sample_times, hp, label='Plus Polarization')\n",
    "pylab.plot(hp.sample_times, hc, label='Cross Polarization')\n",
    "pylab.xlabel('Time (s)')\n",
    "pylab.legend()\n",
    "pylab.grid()\n",
    "pylab.show()\n",
    "\n",
    "# Zoom in near the merger time#\n",
    "pylab.plot(hp.sample_times, hp, label='Plus Polarization')\n",
    "pylab.plot(hp.sample_times, hc, label='Cross Polarization')\n",
    "pylab.xlabel('Time (s)')\n",
    "pylab.xlim(-.01, .01)\n",
    "pylab.legend()\n",
    "pylab.grid()\n",
    "pylab.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that in the this case, the two polariations differ only by the phase of the signal. This holds for systems where the orbital plane of the binary doesn't precess. In the zoom-in plot, we can see the merger itself and the ringdown that follows."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How does the waveform change with the mass of the binary?\n",
    "\n",
    "Below you can see how the length of the waveform increases for lower mass binary mergers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Component mass of each binary component. We'll simplify here and assume that each \n",
    "# component of the binary has the same mass. Again, units are in solar masses.\n",
    "for m in [5, 10, 30, 100]:\n",
    "    hp, hc = get_td_waveform(approximant=\"SEOBNRv4_opt\",\n",
    "                         mass1=m,\n",
    "                         mass2=m,\n",
    "                         delta_t=1.0/4096,\n",
    "                         f_lower=30)\n",
    "    \n",
    "    pylab.plot(hp.sample_times, hp, label='$M_{\\odot 1,2}=%s$' % m)\n",
    "pylab.legend()\n",
    "pylab.grid()\n",
    "pylab.xlabel('Time (s)')\n",
    "pylab.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Changing the distance of the waveform\n",
    "\n",
    "The distance of the waveform is also selectable when you generate a waveform. The units used are Megaparsecs. Keep in mind that no redshift effects are taken into account here, so there is a simple linear relationship between distance and amplitude"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for d in [100, 500, 1000]:\n",
    "    hp, hc = get_td_waveform(approximant=\"SEOBNRv4_opt\",\n",
    "                         mass1=10,\n",
    "                         mass2=10,\n",
    "                         delta_t=1.0/4096,\n",
    "                         f_lower=30,\n",
    "                         distance=d)\n",
    "    \n",
    "    pylab.plot(hp.sample_times, hp, label='$Distance=%sMpc$' % d)\n",
    "pylab.legend()\n",
    "pylab.grid()\n",
    "pylab.xlabel('Time (s)')\n",
    "pylab.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Looking for a specific signal in the data\n",
    "\n",
    "If you know what signal you are looking for in the data, then matched filtering is known to be the optimal method in Gaussian noise to extract the siganl. Even when the parameters of the signal are unkown,  one can test for each set of parameters one is interesting in finding."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### preconditioning the data \n",
    " \n",
    "The purpose of this is to reduce the dynamic range of the data and  supress low freqeuncy behavior which can introduce numerical artefacts. We may also wish to resample the data if high frequency content is not important."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# As an example we use the GW150914 data\n",
    "from pycbc.catalog import Merger\n",
    "from pycbc.filter import resample_to_delta_t, highpass\n",
    "\n",
    "merger = Merger(\"GW150914\")\n",
    "\n",
    "# Get the data from the Hanford detector\n",
    "strain = merger.strain('H1')\n",
    "\n",
    "# Remove the low frequency content and downsample the data to 2048Hz\n",
    "strain = resample_to_delta_t(highpass(strain, 15.0), 1.0/2048)\n",
    "\n",
    "pylab.plot(strain.sample_times, strain)\n",
    "pylab.xlabel('Time (s)')\n",
    "pylab.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### filter wraparound \n",
    "\n",
    "Note the spike in the data at the boundaries. This is caused by the highpass and resampling stages filtering the data. When the filter is applied to the boundaries, it wraps around to the beginning of the data. Since the data itself has a discontinuity (i.e. it is not cyclic) the filter itself will ring off for a time up to the length of the filter. \n",
    "\n",
    "Even if a visible transient is not seen, we want to avoid filters that act on times which are not causally connect. To avoid this we trim the ends of the data sufficiently to ensure that they do not wraparound the input. We will enforce this requirement in all steps of our filtering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove 2 seconds of data from both the beginning and end\n",
    "conditioned = strain.crop(2, 2)\n",
    "\n",
    "pylab.plot(conditioned.sample_times, conditioned)\n",
    "pylab.xlabel('Time (s)')\n",
    "pylab.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### calculate the power spectral density\n",
    "\n",
    "Optimal matched filtering requires weighting the frequency components of the potential signal and data by the noise amplitude. We can view this as filtering the data with the time series equivelant of 1 / PSD. To ensure that we can control how much applying this filter to the data, we window the time domain equivelant of the PSD to a specific length. This has the effect of losing some information about line behavior in the detector. However, since our signal spans a large frequency range, and lines are narrow, this is a negligible effect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pycbc.psd import interpolate, inverse_spectrum_truncation\n",
    "# Estimate the power spectral density\n",
    "\n",
    "# We use 4 second samples of our time series in Welch method.\n",
    "psd = conditioned.psd(4)\n",
    "\n",
    "# Now that we have the psd we need to interpolate it to match our data\n",
    "# and then limit the filter length of 1 / PSD. After this, we can\n",
    "# directly use this PSD to filter the data in a controlled manner\n",
    "\n",
    "psd = interpolate(psd, conditioned.delta_f)\n",
    "\n",
    "# 1/PSD will now act as a filter with an effective length of 4 seconds\n",
    "# Since the data has been highpassed above 15 Hz, and will have low values\n",
    "# below this we need to informat the function to not include frequencies\n",
    "# below this frequency. \n",
    "psd = inverse_spectrum_truncation(psd, 4 * conditioned.sample_rate,\n",
    "                                  low_frequency_cutoff=15)\n",
    "\n",
    "pylab.loglog(psd.sample_frequencies, psd)\n",
    "pylab.ylabel('$Strain^2 / Hz$')\n",
    "pylab.xlabel('Frequency (Hz)')\n",
    "pylab.xlim(30, 1024)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### make your signal model\n",
    "\n",
    "Conceptually, matched filtering involves laying the potential signal over your data and integrating (after weighting frequencies correctly). If there is a signal in the data that aligns with your 'template', you will get a large value when integrated over."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In this case we \"know\" what the signal parameters are. In a search\n",
    "# we would grid over the parameters and calculate the SNR time series\n",
    "# for each one\n",
    "\n",
    "# We'll assume equal masses, which is within the posterior probability\n",
    "# of GW150914. \n",
    "m = 36 # Solar masses\n",
    "hp, hc = get_td_waveform(approximant=\"SEOBNRv4_opt\",\n",
    "                     mass1=m,\n",
    "                     mass2=m,\n",
    "                     delta_t=conditioned.delta_t,\n",
    "                     f_lower=20)\n",
    "\n",
    "# We will resize the vector to match our data\n",
    "hp.resize(len(conditioned))\n",
    "\n",
    "# The waveform begins at the start of the vector, so if we want the\n",
    "# SNR time series to correspond to the approximate merger location\n",
    "# we need to shift the data so that the merger is approximately at the \n",
    "# first bin of the data.\n",
    "\n",
    "# This function rotates the vector by a fixed amount of time.\n",
    "# It treats the data as if it were on a ring. Note that\n",
    "# time stamps are *not* in general affected, but the true\n",
    "# position in the vector is.\n",
    "#\n",
    "# By convention waveforms returned from `get_td_waveform` have their\n",
    "# merger stamped with time zero, so we can use the start time to \n",
    "# shift the merger into position\n",
    "template = hp.cyclic_time_shift(hp.start_time)\n",
    "pylab.plot(template)\n",
    "pylab.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### calculating the signal-to-noise time series\n",
    "\n",
    "In this section we will now calculate the signal-to-noise time series for our template. We'll take care to handle issues of filter corruption / wraparound by truncating the output time series. We need to account for both the length of the template and 1 / PSD."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from pycbc.filter import matched_filter\n",
    "import numpy\n",
    "\n",
    "snr = matched_filter(template, conditioned,\n",
    "                     psd=psd, low_frequency_cutoff=20)\n",
    "\n",
    "# Remove time corrupted by the template filter and the psd filter\n",
    "# We remove 4 seonds at the beginning and end for the PSD filtering\n",
    "# And we remove 4 additional seconds at the beginning to account for\n",
    "# the template length (this is somewhat generous for \n",
    "# so short a template). A longer signal such as from a BNS, would \n",
    "# require much more padding at the beginning of the vector.\n",
    "snr = snr.crop(4 + 4, 4)\n",
    "\n",
    "# Why am I taking an abs() here?\n",
    "# The `matched_filter` function actually returns a 'complex' SNR.\n",
    "# What that means is that the real portion correponds to the SNR\n",
    "# associated with directly filtering the template with the data.\n",
    "# The imaginary portion corresponds to filtering with a template that\n",
    "# is 90 degrees out of phase. Since the phase of a signal may be \n",
    "# anything, we choose to maximize over the phase of the signal.\n",
    "pylab.figure(figsize=[10, 4])\n",
    "pylab.plot(snr.sample_times, abs(snr))\n",
    "pylab.ylabel('Signal-to-noise')\n",
    "pylab.xlabel('Time (s)')\n",
    "pylab.show()\n",
    "\n",
    "peak = abs(snr).numpy().argmax()\n",
    "snrp = snr[peak]\n",
    "time = snr.sample_times[peak]\n",
    "\n",
    "print(\"We found a signal at {}s with SNR {}\".format(time, \n",
    "                                                    abs(snrp)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aligning and Subtracting the Proposed Signal\n",
    "\n",
    "In the previous section we found a peak in the signal-to-noise for a proposed binary black hole merger. We can use this SNR peak to align our proposal to the data, and to also subtract our proposal from the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pycbc.filter import sigma\n",
    "# The time, amplitude, and phase of the SNR peak tell us how to align\n",
    "# our proposed signal with the data.\n",
    "\n",
    "# Shift the template to the peak time\n",
    "dt = time - conditioned.start_time\n",
    "aligned = template.cyclic_time_shift(dt)\n",
    "\n",
    "# scale the template so that it would have SNR 1 in this data\n",
    "aligned /= sigma(aligned, psd=psd, low_frequency_cutoff=20.0)\n",
    "\n",
    "# Scale the template amplitude and phase to the peak value\n",
    "aligned = (aligned.to_frequencyseries() * snrp).to_timeseries()\n",
    "aligned.start_time = conditioned.start_time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualize the overlap between the signal and data\n",
    "\n",
    "To compare the data and signal are on equal footing, and to concentrate on the frequency range that is important, we will whiten both the template and the data, then bandpass both the data and template between 30-300 Hz. In this way, any signal that is in the data is transformed in the same way that the template is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# We do it this way so that we can whiten both the template and the data\n",
    "white_data = (conditioned.to_frequencyseries() / psd**0.5).to_timeseries()\n",
    "\n",
    "# apply a smoothing of the turnon of the template to avoid a transient\n",
    "# from the sharp turn on in the waveform.\n",
    "tapered = aligned.highpass_fir(30, 512, remove_corrupted=False)\n",
    "white_template = (tapered.to_frequencyseries() / psd**0.5).to_timeseries()\n",
    "\n",
    "white_data = white_data.highpass_fir(30., 512).lowpass_fir(300, 512)\n",
    "white_template = white_template.highpass_fir(30, 512).lowpass_fir(300, 512)\n",
    "\n",
    "# Select the time around the merger\n",
    "white_data = white_data.time_slice(merger.time-.2, merger.time+.1)\n",
    "white_template = white_template.time_slice(merger.time-.2, merger.time+.1)\n",
    "\n",
    "pylab.figure(figsize=[15, 3])\n",
    "pylab.plot(white_data.sample_times, white_data, label=\"Data\")\n",
    "pylab.plot(white_template.sample_times, white_template, label=\"Template\")\n",
    "pylab.legend()\n",
    "pylab.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Subtracting the signal from the data\n",
    "\n",
    "Now that we've aligned the template we can simply subtract it. Let's see below how that looks in the time-frequency plots!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subtracted = conditioned - aligned\n",
    "\n",
    "# Plot the original data and the subtracted signal data\n",
    "\n",
    "for data, title in [(conditioned, 'Original H1 Data'),\n",
    "                    (subtracted, 'Signal Subtracted from H1 Data')]:\n",
    "\n",
    "    t, f, p = data.whiten(4, 4).qtransform(.001, logfsteps=100,\n",
    "                                           qrange=(8, 8),\n",
    "                                           frange=(20, 512))\n",
    "    pylab.figure(figsize=[15, 3])\n",
    "    pylab.title(title)\n",
    "    pylab.pcolormesh(t, f, p**0.5, vmin=1, vmax=6)\n",
    "    pylab.yscale('log')\n",
    "    pylab.xlabel('Time (s)')\n",
    "    pylab.ylabel('Frequency (Hz)')\n",
    "    pylab.xlim(merger.time - 2, merger.time + 1)\n",
    "    pylab.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Challenge Problems!\n",
    "\n",
    "Use the methods demonstrated above to see if you can calculate the SNR\n",
    "time series in the following data sets. What is the SNR of each signal?\n",
    "Which template matched which data?\n",
    "\n",
    "Information that may be useful:\n",
    "\n",
    "* Signals are all placed between 100 and 120 seconds into the frame file.\n",
    "* You may assume mass1 = mass2 (equal mass) and that each component mass is one of 22, 36, or 50.\n",
    "* Each file starts at gps time 0, and ends at gps time 128\n",
    "* The channel name in each file is \"H1:TEST-STRAIN\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download the challenge set files\n",
    "\n",
    "import urllib\n",
    "\n",
    "def get_file(fname):\n",
    "    url = \"http://github.com/ahnitz/odw-storage/raw/master/{}\"\n",
    "    url = url.format(fname)\n",
    "    urllib.urlretrieve(url, fname)\n",
    "    print('Getting : {}'.format(url))\n",
    "\n",
    "files = ['PyCBC_T2_0.gwf', 'PyCBC_T2_2.gwf', 'PyCBC_T2_3.gwf']\n",
    "\n",
    "for fname in files:\n",
    "    get_file(fname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}