Skymap tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# This is a test\n",
    "\n",
    "In this tutorial we will learn how:\n",
    "\n",
    "1. to receive the LIGO/Virgo open and public alerts,\n",
    "\n",
    "3. to visualize and compare the sky localizations,\n",
    "\n",
    "4. to cross-match a gravitational-wave sky localization with a galaxy catalog,\n",
    "    \n",
    "5. to search for spatial and temparal coincidents."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LIGO/Virgo Open Public Alerts\n",
    "The third observing run (\"O3\") from LIGO and Virgo Collaborations marked the beginning of the era of public alerts for the field of gravitational-wave astronomy. \n",
    "Automated machine-readable alerts were distributed through the NASA’s Gamma-ray Coordinates Network ([GCN](https://gcn.gsfc.nasa.gov/)) without any human intervention. Shortly afterward, they were vetted by a rapid response team and either confirmed or retracted. Additional information in Section [Alert Timeline](https://emfollow.docs.ligo.org/userguide/analysis/index.html#alert-timeline) in the [Public Alerts User Guide](https://emfollow.docs.ligo.org/userguide/index.html).\n",
    "\n",
    "In the LIGO/Virgo O3 public alerts, fifty-six detection candidates are send to the community and they are collected in the Gravitational-Wave Candidate Event Database ([GraceDB](https://gracedb.ligo.org/)) in the [Public Alert page](https://gracedb.ligo.org/superevents/public/O3/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Receiving GCNs\n",
    "The Gamma-ray Coordinates Network/Transient Astronomy Network (GCN/TAN) is a system for distributing astronomical alerts. There are two main types of alerts: \n",
    "\n",
    "1. GCN Circulars: human-readable astronomical bulletin.\n",
    "2. GCN Notices: A machine-readable alert. \n",
    "\n",
    "Here using the [`pygcn`](https://pypi.org/project/pygcn/) module for receiving GCN/TAN notices in [VOEvent](https://www.ivoa.net/documents/VOEvent/) XML format.\n",
    "Pygcn package implements a simple client that listens for VOEvent XML format notices. By default, it connects to one of the anonymous GCN/TAN server; no sign-up or configuration is necessary to begin receiving alerts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following basic handler function will parse out the URL of the FITS file. The sky map will be read by the [`read_map`](https://healpy.readthedocs.io/en/latest/generated/healpy.fitsfunc.read_map.html?highlight=read_map#healpy.fitsfunc.read_map) method implemented in the [`healpy`](https://pypi.org/project/healpy/) package. Finally we will extract the posterior mean distance and its posterior standard deviation from the image's header. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gcn\n",
    "import healpy as hp\n",
    "\n",
    "@gcn.handlers.include_notice_types(\n",
    "    gcn.notice_types.LVC_EARLY_WARNING,\n",
    "    gcn.notice_types.LVC_PRELIMINARY,\n",
    "    gcn.notice_types.LVC_INITIAL,\n",
    "    gcn.notice_types.LVC_UPDATE,\n",
    "    gcn.notice_types.LVC_RETRACTION)\n",
    "\n",
    "def process_gcn(payload, root):\n",
    "    # Respond only to 'test' events.\n",
    "    # For real observation replace\n",
    "    #  'test' in observation'\n",
    "    if root.attrib['role'] != 'test':\n",
    "        return\n",
    "\n",
    "    # Read all of the VOEvent parameters from the \"What\" section.\n",
    "    params = {elem.attrib['name']:\n",
    "              elem.attrib['value']\n",
    "              for elem in root.iterfind('.//Param')}\n",
    "\n",
    "    if params['AlertType'] == 'Retraction':\n",
    "        print(params['GraceID'], 'was retracted')\n",
    "        return\n",
    "\n",
    "    # Respond only to 'CBC' events. Change 'CBC' to 'Burst'\n",
    "    # to respond to only unmodeled burst events.\n",
    "    if params['Group'] != 'CBC':\n",
    "        return\n",
    "\n",
    "    # Print all parameters.\n",
    "    for key, value in params.items():\n",
    "        print(key, '=', value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### An Offline Example\n",
    "For convenience, we run an offline example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load a VOEvent event.\n",
    "!curl -O https://emfollow.docs.ligo.org/userguide/_static/MS181101ab-2-Preliminary.xml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import lxml.etree\n",
    "\n",
    "# Read the VOevent.\n",
    "payload = open('MS181101ab-2-Preliminary.xml', 'rb').read()\n",
    "root = lxml.etree.fromstring(payload)\n",
    "\n",
    "# Run our handler function.\n",
    "process_gcn(payload, root)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sky localization and interforometer network.\n",
    "Here we show the importance in construction a network of interferometers for pinpointing the gravitational-wave source localization over the sky. The case of the GW190814 will be examined.\n",
    "\n",
    "### Real time improving of the sky localizations of GW190814\n",
    "GW190814 was first identified on 2019 August 14, 21:11:00 UTC as a loud two-detector event in LIGO Livingston and Virgo data (S/N 21.4 and 4.3). A [Notice](https://gcn.gsfc.nasa.gov/notices_l/S190814bv.lvc) was issued through NASA's Gamma-ray Coordinates Network 20 minutes later with a two-detector source localization computed using the rapid Bayesian algorithm BAYESTAR ([Singer & Price 2016](https://journals.aps.org/prd/abstract/10.1103/PhysRevD.93.024013)).\n",
    "\n",
    "Shortly thereafter, reanalyses including LIGO Hanford data were performed. A coincident gravitational-wave signal was identified in all three detectors. Results of these three-detector analyses were reported in the [GCN Circular 25324](https://gcn.gsfc.nasa.gov/gcn3/25324.gcn3) within 2.3 hr of the time of the event providing a three-detector localization [LIGO Scientific Collaboration & Virgo Collaboration 2019](https://iopscience.iop.org/article/10.3847/2041-8213/ab960f).\n",
    "\n",
    "We will apply the method [`ligo_skymap_contour_moc`](https://lscsoft.docs.ligo.org/ligo.skymap/tool/ligo_skymap_contour_moc.html) from [`ligo.skymap`](https://lscsoft.docs.ligo.org/ligo.skymap/) to create the 90% credible regions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download the sky maps of GW190814 from GraceDB.\n",
    "from astropy.utils.data import download_file\n",
    "\n",
    "# Two-detector source localization.\n",
    "gw190814_2ifo = ('https://gracedb.ligo.org/api/superevents/'\n",
    "                 'S190814bv/files/bayestar.multiorder.fits,0')\n",
    "gw190814_2ifo = download_file(gw190814_2ifo)\n",
    "\n",
    "# Three-detector source localization.\n",
    "gw190814_3ifo = ('https://gracedb.ligo.org/api/superevents/'\n",
    "                                      'S190814bv/files/bayestar.multiorder.fits')\n",
    "gw190814_3ifo = download_file(gw190814_3ifo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ligo.skymap.tool.ligo_skymap_contour_moc import main\n",
    "\n",
    "# Define the 90% credible regions of two-detector source localization.\n",
    "main([gw190814_2ifo, '--contour', '90', '--output', 'GW190814 with 2 IFO'])\n",
    "\n",
    "# Define the 90% credible regions of three-detector source localization.\n",
    "main([gw190814_3ifo, '--contour', '90', '--output', 'GW190814 with 3 IFO'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize an Aladin widget at a defined sky's position and with the DSS color survey.\n",
    "import ipyaladin as ipyal\n",
    "\n",
    "aladin = ipyal.Aladin(target='16.91000 -28.1100', fov=180, survey='P/DSS2/color')\n",
    "aladin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show the credible regions in the Aladin widget with 2 different colors.\n",
    "colors = [\"magenta\", \"yellow\"]\n",
    "\n",
    "credible_regions = [\"GW190814 with 2 IFO\", \"GW190814 with 3 IFO\"]\n",
    "\n",
    "for credible_region, color in zip(credible_regions, colors):\n",
    "    aladin.add_moc_from_URL(credible_region, {'color': color, 'opacity': 0.7,\n",
    "                            'adaptativeDisplay': False, 'name':credible_region})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scroll up and visualize the credible region contours in the interactive Aladin widget. From the `Manage Layers`, <img src=\"images/ipyaladin_layer.png\" alt=\"the Layer Button\" style=\"width:30px; display: inline-block;\"/> the 2 credible regions can be independently selected. \n",
    "\n",
    "If you hover the mouse pointer over the credible region plan a message will be displayed above it to show the coverage in sky percentage. How many square degrees? How refined the localization is?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross-Matching the sky localization with a Galaxy Catalog\n",
    "\n",
    "**TO BE FINISHED 3D sky map description**\n",
    "\n",
    "A gravitational-wave sky localization contains additional information related to the distance... \n",
    "\n",
    "Using the [`crossmatch`](https://lscsoft.docs.ligo.org/ligo.skymap/postprocess/crossmatch.html) method in [`ligo.skymap`](https://lscsoft.docs.ligo.org/ligo.skymap/) package, we can extract the galaxies within the 90% credible volume. The galaxy table will be showed in the previous Aladin widget as an new overlay layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astroquery.vizier import VizierClass\n",
    "from astropy.coordinates import SkyCoord\n",
    "from ligo.skymap.io import read_sky_map\n",
    "from ligo.skymap.postprocess import crossmatch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we retrieve the GLADE catalog.\n",
    "\n",
    "**NOTE**. To keep the MyBinder memory, we load a verified subsample from the GLADE catalog.\n",
    "For real applications, download the entire [GLADE](https://vizier.u-strasbg.fr/viz-bin/VizieR?-source=VII/281) catalog."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For real application, uncomment the lines below to load the entire GLADE catalog.\n",
    "\n",
    "#vizier = VizierClass(\n",
    "#    row_limit=-1, columns=['GWGC', '_RAJ2000', '_DEJ2000', 'Dist'])\n",
    "#cat, = vizier.get_catalogs('VII/281/glade2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Comment this line in real application. We use this to save MyBinder memory.\n",
    "from astropy.table import Table\n",
    "cat = Table.read('https://github.com/ggreco77/Tutotest/blob/main/resources/table_tuto?raw=true')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get coordinates and distance columns.\n",
    "coordinates = SkyCoord(cat['_RAJ2000'], cat['_DEJ2000'], cat['Dist'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read the initial GW190814 skymap with 3 IFO.\n",
    "skymap = read_sky_map(gw190814_3ifo, moc=True)\n",
    "\n",
    "# List the galaxies within the 90% credible volume.\n",
    "result = crossmatch(skymap, coordinates)\n",
    "table = cat[result.searched_prob_vol < 0.9]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add the galaxy list in the previous Aladin widget.\n",
    "aladin.add_table(table)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Multi-Messenger Observation of GW170817\n",
    "\n",
    "On August 17, 2017 astronomers around the world were alerted to gravitational waves observed by the Advanced LIGO and Advanced Virgo detectors. This gravitational wave event, now known as GW170817, appeared to be the result of the merger of two neutron stars (BNS). Less than two seconds after the GW170817 signal, NASA's [Fermi](https://fermi.gsfc.nasa.gov/) satellite observed a gamma-ray burst, now known as GRB170817A, and within minutes of these initial detections telescopes around the world began an extensive observing campaign. The [Swope telescope](http://obs.carnegiescience.edu/swope) in Chile was the first to report a bright optical source (SSS17a/AT2017gfo) in the galaxy NGC 4993 and several other teams independently detected the same transient over the next minutes and hours.  \n",
    "GW170817 marks a new era of multi-messenger astronomy, where the same event is observed by both gravitational waves and electromagnetic waves. [The Dawn of Multi-Messenger Astrophysics: Observations of a Binary Neutron Star Merger](http://ligo.org/science/Publication-GW170817MMA/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipyaladin import Aladin\n",
    "\n",
    "from ipywidgets import Layout, Box, widgets\n",
    "\n",
    "a = Aladin(layout=Layout(width='33.33%'), target='NGC4993', fov=0.1)\n",
    "b = Aladin(layout=Layout(width='33.33%'), survey='P/DSS2/red')\n",
    "c = Aladin(layout=Layout(width='33.33%'), survey='P/2MASS/color')\n",
    "\n",
    "# synchronize target between 3 widgets.\n",
    "widgets.jslink((a, 'target'), (b, 'target'))\n",
    "widgets.jslink((b, 'target'), (c, 'target'))\n",
    "\n",
    "# synchronize FoV (zoom level) between 3 widgets.\n",
    "widgets.jslink((a, 'fov'), (b, 'fov'))\n",
    "widgets.jslink((b, 'fov'), (c, 'fov'))\n",
    "\n",
    "items = [a, b, c]\n",
    "\n",
    "box_layout = Layout(display='flex',\n",
    "                    flex_flow='row',\n",
    "                    align_items='stretch',\n",
    "                    border='solid',\n",
    "                    width='100%')\n",
    "box = Box(children=items, layout=box_layout)\n",
    "box"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The three widgtes show the NGC 4993 in 3 different color bands. Press the `Search` button <img src=\"images/ipyaladin_search.png\" alt=\"the Search Button\" style=\"width:30px; display: inline-block;\"/> and type **GW170817**. The reticle will be located to the position of the electromagnetic transient AT2017gfo identified as a kilonova phenomenon."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aladin Sky Atlas\n",
    "**INTERNAL NOTE.** **We are working with the Aladin Beta version but in a few weeks we will use the official release!**\n",
    "\n",
    "Aladin is an interactive sky atlas allowing the user to visualize digitized astronomical images or full surveys, superimpose entries from astronomical catalogues or databases, and interactively access related data and information from the *Simbad* database, the *VizieR* service and other archives for all known astronomical objects in the field.\n",
    "\n",
    "**You now must install the Aladin GUI following those [instructions](<https://aladin.u-strasbg.fr/java/nph-aladin.pl?frame=downloading>).**\n",
    "\n",
    "It would be ideal to increase the maximum memory size used by your Java runtime environment by following the instructions below.\n",
    "\n",
    "Download the [AladinBeta.jar](https://aladin.u-strasbg.fr/java/AladinBeta.jar) from the Aladin download page. Execute it from a terminal by typing:\n",
    "\n",
    "$ java -Xmx2g -jar AladinBeta.jar\n",
    "\n",
    "The flag -Xmx<ammount of memory> specifies the maximum memory allocation pool for a JVM. Here 2GB of memory is allocated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Online Material for working with Aladin\n",
    "\n",
    "**Sky map information.**\n",
    "\n",
    "1. Sky localization of GW170817 binary neutron star merger provided by the LIGO and Virgo Collaborations: \n",
    "  https://dcc.ligo.org/public/0157/P1800381/006/GW170817_skymap.fits.gz\n",
    "\n",
    "\n",
    "2. Error Box of GRB170817 provided by GBM instrument on-board of the Fermi Gamma-ray Space Telescope:              https://gammaray.nsstc.nasa.gov/gbm/science/grbs/grb170817a/gbuts_healpix_systematic.fit\n",
    "\n",
    "**Temporal information.**\n",
    "\n",
    "3. Trigger time of GW170817. DATE-OBS = 2017-08-17T12:41:04.4 / UTC date of the observation.\n",
    "\n",
    "\n",
    "4. Fermi/GBM trigger time of GRB170817: 2017-08-17T12:41:06.  \n",
    "T90 duration = 2.0 ± 0.5 s, starting at T0−0.192 s."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Through the Aladin graphical interface, we will demonstrate the spatial and temporal coincidence between the GW170817 and the short GRB170817. Next, we will retrieve the galaxies inside the intersection region and we search for the host galaxy of the event NGC 4993.**\n",
    "\n",
    "\n",
    "<b>1.</b> Load the sky localization of GW170817.\n",
    "\n",
    "↳ <em>Copy and Paste the [link](https://dcc.ligo.org/public/0157/P1800381/006/GW170817_skymap.fits.gz) in the Command Bar </em> - [Video tutorial]().\n",
    "\n",
    "<b>2.</b> Identify the 90% credible region.\n",
    "   \n",
    "↳ <em>Coverage → Generate a MOC based on... → The current probability skymap → MOC generation window </em> - [Video tutorial]().\n",
    "\n",
    "The plot shows the 90% credible regions. These contours were constructed ranking the pixels from most probable to least, and finally counted how many pixels summed to a given total probability. Place your mouse arrow on the MOC plane (bottom right). You can then read in the description on the top right the solid angle ('Space') associated with the skymap. How many square degrees to you read?\n",
    "   \n",
    "<b>3.</b> Add a time window for coincident with external trigger search [[-1,5](https://emfollow.docs.ligo.org/userguide/analysis/searches.html#coincident-with-external-trigger-search)].\n",
    "   \n",
    "↳ <em>Coverage → Generate a Space-Time MOC based on... → The selected Space MOC... → Properties window </em> - [Video tutorial]()\n",
    "   \n",
    "   \n",
    "<b>4.</b> Load the GRB170817 error box provided by Fermi/GBM.\n",
    "\n",
    "↳ <em>Copy and Paste the [link](https://gammaray.nsstc.nasa.gov/gbm/science/grbs/grb170817a/gbuts_healpix_systematic.fit) in the Command Bar </em>.\n",
    "\n",
    "<b>5.</b> Identify the 90% probability region.\n",
    "    \n",
    "↳ <em>Coverage → Generate a MOC based on... → The current probability skymap → MOC generation window </em> - [Video tutorial]().\n",
    "\n",
    "<b>6.</b> Add the Fermi/GBM trigger time.\n",
    "    \n",
    "↳ <em>Coverage → Generate a Space-Time MOC based on... → The selected Space MOC... → Properties window </em> - [Video tutorial]()\n",
    "    \n",
    "\n",
    "<b>7.</b> Search for temporal and spatial intersection.\n",
    "\n",
    "↳ [Video tutorial]().\n",
    "\n",
    "<b>8.</b> Querying the GLADE galaxy catalog (v2.3) and select the galaxies compatible with the direction of the source.\n",
    "\n",
    "↳ Aladin data collections tree → Select → click on the catalog item → in the popup menu check  → by region&MOC   </em> - [Video tutorial]().\n",
    "\n",
    "<b>9.</b> Filtering the galaxies according with the distance information.\n",
    "\n",
    "↳ Catalog →  Create a filter  → Properties window → Advanced mode → Or enter your filter definition: `${Dist}>31 && ${Dist}<45 {draw}` - [Video tutorial]().\n",
    "\n",
    "<b>10.</b> Search for the host galaxy NGC 4993.\n",
    "\n",
    "↳ <em>Type in the Command bar: NGC 4993</em> - [Video tutorial]()."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "p38workshop",
   "language": "python",
   "name": "p38workshop"
  },
  "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.8.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}