This document explains how to receive, filter, and process gravitational-wave (GW) detection candidate alerts from Advanced LIGO and Virgo. We provide sample code in Python and document alternatives for users of other programming environments. You can download this document and run the code samples in IPython Notebook.

The GCN/TAN system may already be familiar to some, as it is has been in use since the early 1990s to transmit times and coordinates of gamma-ray bursts (GRBs) detected by gamma-ray space missions to ground-based observers. However, there are three peculiariaties of GW observations that veteran GCN users should be aware of:

Unlike GRBs, GW sources cannot (always) be localized to a unique sky location; GW position reconstructions are (sometimes) multimodal and non-Gaussian. As a consequence, LIGO-Virgo GCN notices do not contain a RA, Dec, or error radius; instead they contain a pointer URL to a FITS file containing a probability sky map in the HEALPix all-sky projection.

GCN notices will produce LIGO-Virgo candidates that span a range of significances. The significance of a candidate is reported as the false alarm rate (FAR). The FAR of a detection candidate measures approximately how frequently an event of similar strenght occurs due to chance noise fluctuations or instrumental glitches. A FAR of 1/century ($\sim 3 \times 10^{-10}$ Hz) is generally regarded as sufficient for a "confident" detection. However, alerts will produced for all events with FAR≥1/month ($\sim 4\times10^{-7}$ Hz), such that there should be at least one alert per month on average.

The first few GW alerts will be distributed privately to LIGO MOU partners, over a purpose-built GCN connection.

The first step is to sign up for the GCN network. Signing up for LIGO/Virgo GCN notices is slightly different from the standard signup process. However, if you are already receiving GCN notices (for example, for GRB follow-up), then you can reuse your existing GCN configuration and add the LIGO/Virgo notices.

There are several distribution methods for GCN notices. For the purpose of this tutorial, we will focus on VOEvent over VOEvent Transport Protocol, which is among the more convenient methods for autonomous operations. However, you can use any other distribution method of your choice.

To receive VOEvents, you will need a computer with a static IP address and you will need to register the IP address from which you will connect to the GCN network. Do the following two steps to submit a new GCN site configuration:

If you are registering as a new GCN recipient (as opposed to adding LIGO/Virgo notices to an existing site configuration), then click the "New Site or New User" button and fill out the resulting form. Select the VOEVENT2.0 option for the Distribution Method field.

You will need a robot password to download files that are linked from the GCN notices. This is an automatically generated password that you can use to download files from GraceDb in a script. It is different from your log in for the GraceDb web site, but you use the GraceDb web site to manage your robot password. Here is how to obtain a robot password.

replacing albert.einstein@LIGO.ORG with the robot user name and ABCDEabcde0123456789 with the robot password. This will make it so that you can call curl or other HTTP downloader programs without explicitly providing the user name and password. To test your user name and password, you can run the following command:

IPython Notebook only: If you are following along and running commands in the IPython Notebook, you should also run the following command so that plots are rendered inline inside the notebook:

In [2]:

%matplotlib inline

Next, we'll write a function that we want to get called every time a GCN notice is received. We will use the @gcn.handlers.include_notice_typesfunction decorator to specify that we only want to process certain notice types. There are three notice types:

LVC_PRELIMINARY: Provides the time, significance, and basic parameters about a GW detection candidate. No localization information. Sent with a latency of a minute or so.

LVC_INITIAL: A rapid sky localization is available. Sent with a latency of a few minutes.

LVC_UDPATE: A refined sky localization is availaable. Sent with a latency of hours or more.

In the following example, we will process only the last two types (LVC_INITIAL and LVC_UPDATE), which contain links to sky map FITS files. The following handler function will parse out the URL of the FITS file, download it, and extract the probability sky map.

Events come in two very general flavors: 'CBC' or compact binary coalescence candidates detected by matched filtering, and generic 'Burst' candidates detected by model-independent methods. Most users will want to receive only 'CBC' or only 'Burst' events. In this example code, we are going to keep only 'CBC' events.

Very important: Note that mock or 'test' observations are denoted by the role="test" VOEvent attribute. Alerts resulting from real LIGO/Virgo science data will always have role="observation". The sample code below will respond only to 'test' events. When preparing for actual observations, you must remember to switch to 'observation' events.

In [ ]:

defget_skymap(root):""" Look up URL of sky map in VOEvent XML document, download sky map, and parse FITS file. """# Read out URL of sky map.# This will be something like# https://gracedb.ligo.org/apibasic/events/M131141/files/bayestar.fits.gzskymap_url=root.find("./What/Param[@name='SKYMAP_URL_FITS_BASIC']").attrib['value']# Send HTTP request for sky mapresponse=requests.get(skymap_url,stream=True)# Uncomment to save VOEvent payload to file# open('example.xml', 'w').write(payload)# Raise an exception unless the download succeeded (HTTP 200 OK)response.raise_for_status()# Create a temporary file to store the downloaded FITS filewithtempfile.NamedTemporaryFile()astmpfile:# Save the FITS file to the temporary fileshutil.copyfileobj(response.raw,tmpfile)tmpfile.flush()# Uncomment to save FITS payload to file# shutil.copyfileobj(reponse.raw, open('example.fits.gz', 'wb'))# Read HEALPix data from the temporary fileskymap,header=hp.read_map(tmpfile.name,h=True,verbose=False)header=dict(header)# Done!returnskymap,header# Function to call every time a GCN is received.# Run only for notices of type LVC_INITIAL or LVC_UPDATE.@gcn.handlers.include_notice_types(gcn.notice_types.LVC_INITIAL,gcn.notice_types.LVC_UPDATE)defprocess_gcn(payload,root):# Print the alertprint('Got VOEvent:')print(payload)# Respond only to 'test' events.# VERY IMPORTANT! Replce with the following line of code# to respond to only real 'observation' events.# if root.attrib['role'] != 'observation': returnifroot.attrib['role']!='test':return# Respond only to 'CBC' events. Change 'CBC' to "Burst' to respond to only# unmodeled burst events.ifroot.find("./What/Param[@name='Group']").attrib['value']!='CBC':return# Read out integer notice type (note: not doing anythin with this right now)notice_type=int(root.find("./What/Param[@name='Packet_Type']").attrib['value'])# Read sky mapskymap,header=get_skymap(root)

Finally, we will listen for GCNs. You need to tell the gcn.listen function on what port you want to connect to the GCN network; this will be port 8096. You also need to tell gcn.listen which function to call whenever it receives an GCN; this will be the process_gcn function that we just defined.

When you run the following code snippet, the gcn.listen function will continue until you interrupt the program (by pressing the stop button in IPython Notebook, typing Control-C in the terminal, or sending a kill signal to your Python script).

Note: gcn.listen will automatically reconnect to the GCN network if the network connection is broken.

In [ ]:

# Listen for GCNs until the program is interrupted# (killed or interrupted with control-C).gcn.listen(port=8096,handler=process_gcn)

That was fun! Now kill the listener: if you are running the example code in an IPython Notebook, press the stop button. If you copied the example code into a standalone Python script, then kill the script by typing control-C.

Let's take a look at what is inside one of the LIGO/Virgo probability sky maps. They are FITS image files and can be manipulated and viewed with many commonplace FITS tools. However, they are a little unusual in two regards. First, since they are all-sky images, they are stored in the HEALPix projection, a format that is used for Planck all-sky CMB maps and by Aladin for archival all-sky survey images. Second, the value stored at each pixel is the probability that the gravitational-wave source is within that pixel.

Now, let's go through some examples of manipulating HEALPix sky maps programmatically. The HEALPix project provides official libraries for many languages, including C, C++, Fortran, IDL, Java, and Python. However, since this is a Python tutorial, we are going to demonstrate how to manipulate HEALPix maps with the official Python library, Healpy.

First, if you have not already downloaded an example sky map, you can do so now by having Python call curl on the command line:

You can suppress printing informational messages while loading the file by passing the keyword argument verbose=False. You can read both the HEALPix image data and the FITS header by passing the h=True keyword argument:

Each entry in the array represents the probability contained within a quadrilateral pixel whose position on the sky is uniquely specified by the index in the array and the array's length. Because HEALPix pixels are equal area, we can find the number of pixels per square degree just from the length of the HEALPix array:

In [8]:

npix=len(hpx)sky_area=4*180**2/np.pisky_area/npix

Out[8]:

0.013113963206424481

The pix2ang function converts from pixel index to spherical polar coordinates; the function ang2pix does the reverse.

Both pix2ang and ang2pix take, as their first argument, nside, the lateral resolution fo the HEALPix map. You can find nside from the length of the image array by calling npix2nside:

In [9]:

nside=hp.npix2nside(npix)nside

Out[9]:

512

Let's look up the right ascension and declination of pixel number 123. We'll call pix2ang to get the spherical polar coordinates $(\theta, \phi)$ in radians, and then use Numpy's rad2deg function to convert these to right ascension and declination in degrees.

How do we find the probability that the source is contained within a circle on the sky? First we find the pixels that are contained within the circle using query_disc. Note that query_disc takes as its arguments the Cartesian coordinates of the center of the circle, and its radius in radians.

Then, we sum the values of the HEALPix image array contained at those pixels.

In [14]:

# RA, Dec, and radius of circle in degreesra=213.22dec=-37.45radius=3.1# Spherical polar coordinates and radius of circle in radianstheta=0.5*np.pi-np.deg2rad(dec)phi=np.deg2rad(ra)radius=np.deg2rad(radius)# Cartesian coordinates of center of circlexyz=hp.ang2vec(theta,phi)# Array of indices of pixels inside circleipix_disc=hp.query_disc(nside,xyz,radius)# Probability that source is within circlehpx[ipix_disc].sum()

Out[14]:

0.056844148994230181

Similarly, we can use the query_polygon function to look up the indices of the pixels within a polygon (defined by the Cartesian coordinates of its vertices), and then compute the probability that the source is inside that polygon by summing the values of the pixels.

In [15]:

# Vertices of polygonxyz=[[-0.69601758,-0.41315628,-0.58724902],[-0.68590811,-0.40679797,-0.60336181],[-0.69106913,-0.39820114,-0.60320752],[-0.7011786,-0.40455945,-0.58709473]]# Array of indices of pixels inside polygonipix_poly=hp.query_polygon(nside,xyz)# Probability that source is within polygonhpx[ipix_poly].sum()

Out[15]:

0.0011913816551896161

These are all of the HEALPix functions from Healpy that we will need for the remainder of the this tutorial.

Other useful Healpy functions include ud_grade for upsampling or downsampling a sky map, and get_interp_val for performing bilinear interpolation between pixels. See the Healpy tutorial for other useful operations.

The LIGO/Virgo probability sky maps are always in equatorial coordinates. Once we have looked up the coordinates of the HEALPix pixels, we will use the positional astronomy features of Astropy to transform those coordinates to an alt/az frame for a particular site on the Earth at a particular time. Then we can quickly determine which pixels are visible from that site at that time, and integrate (sum) the probability contained in those pixels.

Note: users may want to do something more sophisticated like determine how much of the probability is visible for at least a certain length of time. This example will illustrate one key function of HEALPix (looking up coordinates of the grid with hp.pix2ang) and some of the key positional astronomy functions with Astropy.

In [ ]:

defprob_observable(m,header):""" Determine the integrated probability contained in a gravitational-wave sky map that is observable from a particular ground-based site at a particular time. Bonus: make a plot of probability versus UTC time! """# Determine resolution of sky mapnpix=len(m)nside=hp.npix2nside(npix)# Get time nowtime=astropy.time.Time.now()# Or at the time of the gravitational-wave event...# time = astropy.time.Time(header['MJD-OBS'], format='mjd')# Or at a particular time...# time = astropy.time.Time('2015-03-01 13:55:27')# Geodetic coordinates of observatory (example here: Mount Wilson)observatory=astropy.coordinates.EarthLocation(lat=34.2247*u.deg,lon=-118.0572*u.deg,height=1742*u.m)# Alt/az reference frame at observatory, nowframe=astropy.coordinates.AltAz(obstime=time,location=observatory)# Look up (celestial) spherical polar coordinates of HEALPix grid.theta,phi=hp.pix2ang(nside,np.arange(npix))# Convert to RA, Dec.radecs=astropy.coordinates.SkyCoord(ra=phi*u.rad,dec=(0.5*np.pi-theta)*u.rad)# Transform grid to alt/az coordinates at observatory, nowaltaz=radecs.transform_to(frame)# Where is the sun, now?sun_altaz=astropy.coordinates.get_sun(time).transform_to(altaz)# How likely is it that the (true, unknown) location of the source# is within the area that is visible, now? Demand that sun is at# least 18 degrees below the horizon and that the airmass# (secant of zenith angle approximation) is at most 2.5.prob=m[(sun_altaz.alt<=-18*u.deg)&(altaz.secz<=2.5)].sum()# Done!returnprob

Finally, we need to update our GCN handler to call this function.

In [ ]:

# Function to call every time a GCN is received.# Run only for notices of type LVC_INITIAL or LVC_UPDATE.@gcn.handlers.include_notice_types(gcn.notice_types.LVC_INITIAL,gcn.notice_types.LVC_UPDATE)defprocess_gcn(payload,root):# Print the alertprint('Got VOEvent:')print(payload)# Respond only to 'test' events.# VERY IMPORTANT! Replce with the following line of code# to respond to only real 'observation' events.# if root.attrib['role'] != 'observation': returnifroot.attrib['role']!='test':return# Respond only to 'CBC' events. Change 'CBC' to "Burst' to respond to only# unmodeled burst events.ifroot.find("./What/Param[@name='Group']").attrib['value']!='CBC':returnskymap,header=get_skymap(root)prob=prob_observable(skymap,header)print('Source has a %d%% chance of being observable now'%round(100*prob))ifprob>0.5:pass# FIXME: perform some action

Let's run the new GCN handler now...

In [ ]:

# Listen for GCNs until the program is interrupted# (killed or interrupted with control-C).gcn.listen(port=8096,handler=process_gcn)

Suppose you have performed some EM observations to follow up on a candidate GW event, and you now want to supply the coordinates of those observations to the LV-EM group. The first step is to obtain a robotic access password from GraceDB and add it to a netrc file (see Section 2). One can then use the Python GraceDB client in order to submit a list of coordinates to GraceDB. To install the GraceDB client package:

$ pip install ligo-gracedb

Note: It is highly recommended to install the GraceDB client (as well as other packages described used in this tutorial) inside a virtual environment. This isolates the packages required for interacting with GraceDB from other packages on the system.
Now that the GraceDB client is installed, one can use a script such as this to submit observation records to GraceDB:

In the above example, the lists of ra and dec values were chosen arbitrarily, as well as the start times and durations (or exposure times). Since the comma-separated lists have three elements, there will be three footprints associated with this observation. And since only one value was provided for the raWidth and decWidth, these widths will be assumed to apply to all three footprints.

However, the data argument in the above curl command is rather unwieldy--and this is only with three simple footprints. So using Python (as in the above code example) or your own favorite language and HTTP library is much preferred. We are also planning an email submission mechanism which may be eaiser for some users.

Here is a complete, working GCN processing script. Copy it into a .py file and customize as needed.

In [ ]:

# Python standard library importsimporttempfileimportshutilimportsysimportglob# Third-party importsimportgcnimportgcn.handlersimportgcn.notice_typesimportrequestsimporthealpyashpimportnumpyasnpimportastropy.coordinatesimportastropy.timeimportastropy.unitsasudefget_skymap(root):""" Look up URL of sky map in VOEvent XML document, download sky map, and parse FITS file. """# Read out URL of sky map.# This will be something like# https://gracedb.ligo.org/apibasic/events/M131141/files/bayestar.fits.gzskymap_url=root.find("./What/Param[@name='SKYMAP_URL_FITS_BASIC']").attrib['value']# Send HTTP request for sky mapresponse=requests.get(skymap_url,stream=True)# Uncomment to save VOEvent payload to file# open('example.xml', 'w').write(payload)# Raise an exception unless the download succeeded (HTTP 200 OK)response.raise_for_status()# Create a temporary file to store the downloaded FITS filewithtempfile.NamedTemporaryFile()astmpfile:# Save the FITS file to the temporary fileshutil.copyfileobj(response.raw,tmpfile)tmpfile.flush()# Uncomment to save FITS payload to file# shutil.copyfileobj(reponse.raw, open('example.fits.gz', 'wb'))# Read HEALPix data from the temporary fileskymap,header=hp.read_map(tmpfile.name,h=True,verbose=False)header=dict(header)# Done!returnskymap,headerdefprob_observable(m,header):""" Determine the integrated probability contained in a gravitational-wave sky map that is observable from a particular ground-based site at a particular time. Bonus: make a plot of probability versus UTC time! """# Determine resolution of sky mapnpix=len(m)nside=hp.npix2nside(npix)# Get time nowtime=astropy.time.Time.now()# Or at the time of the gravitational-wave event...# time = astropy.time.Time(header['MJD-OBS'], format='mjd')# Or at a particular time...# time = astropy.time.Time('2015-03-01 13:55:27')# Geodetic coordinates of observatory (example here: Mount Wilson)observatory=astropy.coordinates.EarthLocation(lat=34.2247*u.deg,lon=-118.0572*u.deg,height=1742*u.m)# Alt/az reference frame at observatory, nowframe=astropy.coordinates.AltAz(obstime=time,location=observatory)# Look up (celestial) spherical polar coordinates of HEALPix grid.theta,phi=hp.pix2ang(nside,np.arange(npix))# Convert to RA, Dec.radecs=astropy.coordinates.SkyCoord(ra=phi*u.rad,dec=(0.5*np.pi-theta)*u.rad)# Transform grid to alt/az coordinates at observatory, nowaltaz=radecs.transform_to(frame)# Where is the sun, now?sun_altaz=astropy.coordinates.get_sun(time).transform_to(altaz)# How likely is it that the (true, unknown) location of the source# is within the area that is visible, now? Demand that sun is at# least 18 degrees below the horizon and that the airmass# (secant of zenith angle approximation) is at most 2.5.prob=m[(sun_altaz.alt<=-18*u.deg)&(altaz.secz<=2.5)].sum()# Done!returnprob# Function to call every time a GCN is received.# Run only for notices of type LVC_INITIAL or LVC_UPDATE.@gcn.handlers.include_notice_types(gcn.notice_types.LVC_INITIAL,gcn.notice_types.LVC_UPDATE)defprocess_gcn(payload,root):# Print the alertprint('Got VOEvent:')print(payload)# Respond only to 'test' events.# VERY IMPORTANT! Replce with the following line of code# to respond to only real 'observation' events.# if root.attrib['role'] != 'observation': returnifroot.attrib['role']!='test':return# Respond only to 'CBC' events. Change 'CBC' to "Burst' to respond to only# unmodeled burst events.ifroot.find("./What/Param[@name='Group']").attrib['value']!='CBC':returnskymap,header=get_skymap(root)prob=prob_observable(skymap,header)print('Source has a %d%% chance of being observable now'%round(100*prob))ifprob>0.5:pass# FIXME: perform some action# Listen for GCNs until the program is interrupted# (killed or interrupted with control-C).gcn.listen(port=8096,handler=process_gcn)