LIGO/Virgo Public Alerts User Guide

GW170817/GRB170817A discovery images

Welcome to the LIGO/Virgo Public Alerts User Guide! This document is intended for both professional astronomers and science enthusiasts who are interested in receiving alerts and real-time data products related to gravitational-wave (GW) events.

Three sites (LHO, LLO, Virgo) together form a global network of ground-based GW detectors. The LIGO Scientific Collaboration and the Virgo Collaboration jointly analyze the data in real time to detect and localize transients from compact binary mergers and other sources. When a signal candidate is found, an alert is sent to astronomers in order to search for counterparts (electromagnetic waves or neutrinos).

Advanced LIGO and Advanced Virgo began their third observing run (O3) on April 1, 2019. For the first time, LIGO/Virgo alerts are public. Alerts are distributed through NASA’s Gamma-ray Coordinates Network (GCN). There are two types of alerts: human-readable GCN Circulars and machine-readable GCN Notices. This document provides a brief overview of the procedures for vetting and sending GW alerts, describes their contents and format, and includes instructions and sample code for receiving GCN Notices and decoding GW sky maps.

Contents

Getting Started Checklist

In addition to reading this user guide, we encourage you to sign up for certain e-mail lists and community resources. Make sure that you have completed the items on this checklist.

1. Read This User Guide

Pay particular attention to the Alert Contents section of this user guide to familiarize yourself with the contents of machine-readable LIGO/Virgo GCN Notices. Play with the Sample Code to receive example GCN Notices and practice working with sky localization maps.

2. Subscribe to GCN Circulars

Subscribe to GCN Circulars and review the instructions for posting GCN Circulars. A GCN Circular is a short, public bulletin to rapidly report an astronomical observation. GCN Circulars are distributed by email and archived online. [1] LIGO/Virgo uses GCN Circulars to announce detections, and the astronomy community expects participants to promptly disseminate preliminary reports of follow-up observations of LIGO/Virgo counterparts using GCN Circulars as well.

Important

GCN Circulars can only be posted from registered email addresses. You must sign up for GCN Circulars in advance in order to post to the list.

3. Join the OpenLVEM Community

Sign up to the OpenLVEM mailing list by following the OpenLVEM instructions. LIGO/Virgo will use this list to make announcements and solicit input. It is also a great place to ask questions or discuss issues related to LIGO/Virgo public alerts. Documents relating to teleconferences and in-person meetings are available at OpenLVEM wiki.

4. Visit GraceDB

Familiarize yourself with GraceDB, LIGO/Virgo’s online portal for alerts and real-time results.

[1]GCN Circulars are similar to Astronomer’s Telegrams (ATels). By longstanding convention, the gamma-ray burst and gravitational-wave astronomy communities use GCNs, whereas the supernova community uses ATels.

Observing Capabilities

This section summarizes the observing capabilities of the global gravitational-wave detector network as of early 2019. This as a quick reference to the anticipated commissioning and observing schedule, sensitivity to gravitational-wave transients, and sky localization accuracy, as described in the following external documents:

  • White Paper [1] on gravitational-wave data analysis and astrophysics
  • Living Review [2] on prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo, and KAGRA
  • Working Plan Towards O3 [3]

Timeline

The gravitational-wave observing schedules is divided into Observing Runs, epochs of months to years of operation at fixed sensitivity; down time for construction and commissioning; and transitional Engineering Runs between commissioning and observing runs. The long-term observing schedule is shown below. Since BNS mergers are a well-studied class of gravitational-wave signals, this figure gives the BNS range for each observing run.

Long-term observing schedule

Engineering Run 14 (ER14) started on 2019-03-04. The transition into Observing Run 3 (O3) occurred on 2019-04-01.

During O3, we expect that three facilities (LHO, LLO, and Virgo) will observe for one year. It is possible that the Japanese KAGRA detector may come online and become part of the international gravitational-wave network at some point during O3. The near-term observing schedule is shown below, reproduced from [3].

Current observing schedule

Live Status

There are a handful of public web pages that report live status of the LIGO/Virgo detectors and alert infrastructure.

Sensitivity and Sky Localization Accuracy

The following O3 projections are adapted from the Living Review [2] on prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo, and KAGRA. The BNS range and burst range (luminosity distance of detectable sources, averaged over sky position and source orientation) given the detectors’ anticipated sensitivities are listed in the table below.

Detector Range (Mpc)
BNS Burst
LIGO 120–170 75–90
Virgo 65–85 40–50
KAGRA 8–25  

We expect 1–10 BNS events over the course of O3. For BNS events, the median sky localization accuracy of in terms of the 90% credible area will be 120–180 deg². 12–21% of BNS mergers will be localized to less than 20 deg².

Event Rates

See the rates presentation from recent Town Hall meeting for LIGO and Virgo’s most current estimates of astrophysical rates of compact binary mergers.

[1]LIGO Scientific Collaboration & Virgo Collaboration 2018, The LSC-Virgo White Paper on Gravitational Wave Data Analysis and Astrophysics, LIGO Document T1800058-v1. https://dcc.ligo.org/LIGO-T1800058/public
[2](1, 2) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, Living Rev. Rel., 21, 3. doi:10.1007/s41114-018-0012-9
[3](1, 2) LIGO Scientific Collaboration & Virgo Collaboration 2018, LIGO-Virgo Working Plan Towards O3, LIGO Document G1801056-v4. https://dcc.ligo.org/LIGO-G1801056/public

Data Analysis

In this section we describe the different online searches looking for GW signals, the selection and vetting of candidates, and parameter estimation analysis.

When multiple candidates from different pipelines are close enough together in time, they will be considered as originating from the same physical event and will be grouped into a single superevent. See the following pages for technical details.

Online Pipelines

A number of search pipelines run in a low latency, online mode. These can be divided into two groups, modeled and unmodeled. The modeled (CBC) searches specifically look for signals from compact binary mergers of neutron stars and black holes (BNS, NSBH, and BBH systems). The unmodeled (Burst) searches on the other hand, are capable of detecting signals from a wide variety of astrophysical sources in addition to compact binary mergers: core-collapse of massive stars, magnetar star-quakes, and more speculative sources such as intersecting cosmic strings or as-yet unknown GW sources.

Superevents

Superevents are an abstraction to unify gravitational-wave candidates from multiple search pipelines. Each superevent is intended to represent a single astrophysical event.

A superevent consists of one or more event candidates, possibly from different pipelines, that are neighbors in time. At any given time, one event belonging to the superevent is identified as the preferred event. The superevent inherits properties from the preferred event such as time, significance, sky localization, and classification.

The preferred event may change after a preliminary alert has been sent, but the name of the superevent will stay the same. The naming scheme is described in the alert contents section.

Selection of the Preferred Event

When multiple online searches report events at the same time, the preferred event is decided by applying the following rules, in order:

  1. An event detected in multiple interferometers is preferred over an event from a single interferometer.
  2. An event from modeled CBC searches are preferred over an event from unmodeled Burst searches (see Searches for details on search pipelines).
  3. In the case of multiple CBC events, the event with the highest SNR is preferred. In the case of multiple Burst events, the event with the lowest false alarm rate (FAR) is preferred.

Note

  • A Preliminary GCN is automatically issued for superevents when the false alarm rate is lower than a threshold value.
  • In case of an event created by a pipeline due to an offline analysis, no preliminary GCN will be sent.
  • The SNR is used to select the preferred event among CBC candidates because higher SNR implies better sky location and parameter estimates from low-latency searches.

Candidate Vetting

GWCelery orchestrates and supervises performing basic data quality and detector state checks, grouping of events from individual pipelines into superevents, initiating automated sky localization and parameter estimation, inferring classification and source properties, and sending alerts.

A Data Quality Report is prepared that consists of a semi-automated detector characterization and data quality investigation for each event. It provides a variety of metrics based on auxiliary instrumental and environmental sensors that are used by a rapid response team in order to make a decision of whether to confirm or retract a candidate.

The rapid response team consists of commissioning, computing, and calibration experts from each of the detector sites, pipeline experts, detector characterization experts, and follow-up advocates.

Sky Localization and Parameter Estimation

Immediately after one of the search pipelines reports an event, sky localization and parameter estimation analyses begin. These analyses all use Bayesian inference to calculate the posterior probability distribution over the parameters (sky location, distance, and/or intrinsic properties of the source) given the observed gravitational-wave signal.

There are different parameter estimation methods for modeled (CBC) and unmodeled (Burst) events. However, in both cases there is a rapid analysis that estimates only the sky localization, and is ready in seconds, and a refined analysis that explores a larger parameter space and completes up to hours or a day later.

Modeled Events

BAYESTAR [1] is the rapid CBC sky localization algorithm. It reads in the matched-filter time series from the search pipeline and calculates the posterior probability distribution over the sky location and distance of the source by coherently modeling the response of the gravitational-wave detector network. It explores the parameter space using Gaussian quadrature, lookup tables, and sampling on an adaptively refined HEALPix grid. The sky localization takes tens of seconds and is included in the preliminary alert.

LALInference [2] is the full CBC parameter estimation algorithm. It explores a greatly expanded parameter space including sky location, distance, masses, and spins, and performs full forward modeling of the gravitational-wave signal and the strain calibration of the gravitational-wave detectors. It explores the parameter space using Markov chain Monte Carlo and nested sampling. For all events, there is an automated LALInference analysis that uses the least expensive CBC waveform models and completes within hours and may be included in a subsequent alert. More time-consuming analyses with more sophisticated waveform models are started at the discretion of human analysts, and will complete days or weeks later.

Unmodeled Events

cWB, the burst search pipeline, also performs a rapid sky localization based on its coherent reconstruction of the gravitational-wave signal using a wavelet basis and the response of the gravitational-wave detector network [3]. The cWB sky localization is included in the preliminary alert.

Refined sky localizations for unmodeled bursts are provided by two algorithms that use the same MCMC and nested sampling methodology as LALInference. LALInference Burst (LIB) [4] models the signal as a single sinusoidally modulated Gaussian. BayesWave [5] models the signal as a superposition of wavelets and jointly models the background with both a stationary noise component and glitches composed of wavelets that are present in individual detectors.

[1]Singer, L. P., & Price, L. R. 2016, Phys. Rev. D, 93, 024013. doi:10.1103/PhysRevD.93.024013
[2]Veitch, J., Raymond, V., Farr, B., et al. 2015, Phys. Rev. D, 91, 042003. doi:10.1103/PhysRevD.91.042003
[3]Klimenko, S., Vedovato, G., Drago, M., et al. 2011, Phys. Rev. D, 83, 102001. doi:10.1103/PhysRevD.83.102001
[4]Lynch, R., Vitale, S., Essick, R., Katsavounidis, E., & Robinet, F. 2017, Phys. Rev. D, 95, 104046. doi:10.1103/PhysRevD.95.104046
[5]Cornish, N. J., & Littenberg, T. B. 2015, Class. Quantum Grav., 32, 135012. doi:10.1088/0264-9381/32/13/135012

Inference

For CBC events, we calculate a number of quantities that are inferred from the signal. In preliminary alerts, these quantities are based on the candidate significance and the matched-filter estimates of the source parameters. Once parameter estimation has been completed, updated values will be provided based on samples drawn from the posterior probability distribution.

Classification

The classification consists of five numbers, summing to unity, that give the probability that the source is a BNS, NSBH, BBH merger, or contains at least one MassGap component, or is of terrestrial (i.e. a background fluctuation or a glitch) origin. See the figure in the alert contents section for the boundaries of the source classification categories in the \((m_1, m_2)\) plane.

This assumes that terrestrial and astrophysical events occur as independent Poisson processes. A source-dependent weighting of matched-filter templates is used to compute the mean values of expected counts associated with each of these five categories. The mean values are updated weekly based on observed matched-filter count rates. They are then used to predict the category for new triggers uploaded by search pipelines.

For details, see [1], especially Section II E. and Equation 18 therein.

Properties

The source properties consist of a set of numbers, each between zero and unity, that give the probabilities that the source satisfies certain conditions. These conditions are:

HasNS: At least one of the compact objects in the binary (that is, the less massive or secondary compact object) has a mass that is consistent with a neutron star. Specifically, we define this as the probability that the secondary mass satisfies \(m_2 \leq 3 M_{\odot}\).

HasRemnant: The probability that the source ejected a nonzero mass outside the final remnant compact object. This is calculated using the disk mass fitting formula from [2] (Equation 4).

For preliminary estimates based on the matched-filter pipeline results, the source properties are calculated using a supervised learning technique; mass-dependent rates are the same as those used for the classification. For parameter estimation, updated properties are calculated from posterior samples.

[1]Kapadia, S. J., Caudill, S., Creighton, J. D. E., et al., 2019. arXiv:1903.06881
[2]Foucart, F., Hinderer, T. & Nissanke, S. 2018, Phys. Rev. D, 98, 081501. doi:10.1103/PhysRevD.98.081501

The timeline for distribution of alerts is described below.

Alert Timeline

Here, we describe the sequence of LIGO/Virgo alerts for a single event that will be distributed through the Gamma-ray Coordinates Network (GCN) via notices and circulars (see the Alert Contents and Sample Code sections for details).

Timeline for sending gravitational-wave alerts

Within 1–10 minutes after GW trigger time, the first preliminary notice will be sent fully autonomously. Preliminary notices may or may not include a sky localization. The trigger will be immediately and publicly visible in the GraceDB database. Since the procedure is fully automatic, some preliminary alerts may be retracted after human inspection for data quality, instrumental conditions, and pipeline behavior.

Within 24 hours after the GW trigger time (possibly within 4 hours for BNS or NSBH sources, to be decided), the Initial or Retraction notice and circular will be distributed. It will include an updated sky localization and source classification. At this stage, the event will have been vetted by human instrument scientists and analysts. The candidate will either be confirmed by an Initial notice and circular or withdrawn by a Retraction notice and circular if the data quality is unsuitable.

Within a day, black hole mergers will be fully vetted by experts and retraction or confirmation status will be reported.

Update notice and circulars are sent whenever the sky localization area or significance accuracy improves (e.g. as a result of improved calibration, glitch removal, or computationally deeper parameter estimation). Updates will be sent up until the position is determined more accurately by public announcement of an unambiguous counterpart. At that point, there will be no further sky localization updates until the publication of the event in a peer-reviewed journal.

At any time, we can promote an extraordinary candidate that does not pass our public alert thresholds if it is compellingly associated with a multimessenger signal (e.g. GRB, core-collapse SN). In this case, Initial notices and circulars will be distributed.

Alert Contents

Public LIGO/Virgo alerts are distributing using NASA’s Gamma-ray Coordinates Network (GCN). There are two types of alerts:

GCN Notices are machine-readable packets. They are available as VOEvent XML and several other formats. We strongly recommend receiving LIGO/Virgo alerts in the VOEvent XML format.

GCN Circulars are short human-readable astronomical bulletins. They are written in a certain well-established format and style. You can subscribe to GCN Circulars to receive and post them by email, or you can view them in the public GCN Circulars archive.

Notice Types

For each event, there are up to four kinds of GCN Notices:

A Preliminary GCN Notice is issued automatically within minutes after a gravitational-wave candidate is detected. The candidate must have passed some automated data quality checks, but it may later be retracted after human vetting. There is no accompanying GCN Circular at this stage.

An Initial GCN Notice is issued after human vetting (see Candidate Vetting). If the signal does not pass human vetting (e.g., it is a glitch), then instead of an initial alert there will be a retraction. The initial alert is also accompanied by a GCN Circular which should be considered as the first formal publication of the candidate and can be cited as such.

An Update GCN Notice is issued whenever further analysis leads to improved estimates of the sky localization, significance, or classification. There may be multiple updates for a given event, and updates may be issued hours, days, or even weeks after the event.

Lastly, a Retraction GCN Notice is issued if the candidate is rejected as a result of vetting by human instrument scientists and data analysts. A retraction indicates that the candidate has been withdrawn because it is probably not astrophysical.

All types of GCN Notices except for Retraction notices contain the following information, which are described in further detail below:

  • Name: a unique identifier for the candidate
  • Significance: estimated false alarm rate
  • Sky localization: inferred sky position and (CBC candidates only) distance
  • Inference: inferred source classification and properties (CBC candidates only)

Of the above fields, Retraction notices provide only the name.

Initial and Update notices are accompanied by human-readable GCN Circulars, which restate all of the above information and also may include a data quality assessment.

Notice Contents

The table below is a representation of the contents of a LIGO/Virgo GCN Notice.

Root
IVORN ivo://gwnet/LVC#[{T,M}]SYYMMDDabc-{1,2,3}-{Preliminary,Initial,Update,Retraction}
Role {observation,test}
Who
Date Time sent (UTC, ISO-8601), e.g. 2018-11-01T22:34:49
Author LIGO Scientific Collaboration and Virgo Collaboration
WhereWhen Time of signal (UTC, ISO-8601), e.g. 2018-11-01T22:22:46.654437
What
GraceID GraceDB ID: [{T,M}]SYYMMDDabc. Example: MS181101abc
Packet Type GCN Notice type: {Preliminary,Initial,Update,Retraction}
Notice Type Numerical equivalent of GCN Notice type: {150,151,152,164}
FAR Estimated false alarm rate in Hz
Sky Map URL of HEALPix FITS sky localization file
Group CBC Burst
Pipeline {Gstlal,MBTAOnline,PyCBC,SPIIR} {cWB,oLIB}
CentralFreq N/A Central frequency in Hz
Duration Duration of burst in s
Fluence Gravitational-wave fluence in erg cm\(^{-2}\)
BNS, NSBH, BBH, MassGap, Noise Probability that the source is a BNS, NSBH, BBH, or MassGap merger, or terrestrial (i.e, noise) respectively N/A
HasNS, HasRemnant Probability, under the assumption that the source is not noise, that at least one of the compact objects was a neutron star, and that the system ejected a non-zero amount of neutron star matter, respectively.
Name

The name of an event is its GraceDB ID, a uniquely assigned identifier such as MS181101abc. A GraceDB ID has three parts:

  • Prefix: S for normal candidates and MS or TS for mock or test events respectively. The S stands for superevent.
  • Date: The six-digit UTC date of the event consisting of a two-digit year, month, and day of month.
  • Suffix: A lowercase alphabetic string that is incremented automatically (a, b, …, z, aa, ab, … az, aaa, etc.) whenever a candidate on a given date is added to GraceDB.
Significance

The significance of the event is quantified by its false alarm rate (FAR): the expected rate of events from the pipeline that produced the preferred event with equal or greater significance in the absence of any astrophysical signals.

Sky Localization

The sky localization consists of the posterior probability distribution of the source’s sky position and (for CBC events only) luminosity distance. The GCN Notice and Circular will provide a URL for the sky localization file stored in GraceDB. The sky localization is saved in a FITS file as a HEALPix all-sky image. See our sample code for instructions on working with sky localization files.

Important

We generally provide localizations in two HEALPix formats, distinguished by file extension:

*.fits.gz

A subset of the standard HEALPix-in-FITS format (see semi-official specifications from the HEALPix team and from the gamma-ray community) that is recognized by a wide variety of astronomical imaging programs including DS9 and Aladin. It uses HEALPix implicit indexing and the NESTED numbering scheme. This is the primary and preferred format, and the only format that is explicitly listed in the GCN Notices and Circulars. See the section Working with Sky Maps for details.

*.multiorder.fits

A new variant of the HEALPix format that is designed to overcome limitations of the *.fits.gz format for well-localized events from three-detector operations and future gravitational-wave facilities (see rationale). It uses HEALPix explicit indexing and the NUNIQ numbering scheme, which is closely related to multi-order coverage (MOC) maps in Aladin. This is the internal format that is used by the LIGO/Virgo low-latency alert pipeline. This is an experimental format, and it is currently recommended only for advanced users. See the section Multi-Order Sky Maps (For Advanced Users) for details.

Both formats always use celestial (equatorial, J2000) coordinates.

Inference

The inference section is present for CBC events only. It has two parts:

Classification: Five numbers, summing to unity, giving probability that the source belongs to the following five categories:

The figure below shows the extent of the three astrophysical categories (BNS, NSBH, BBH, and MassGap) in terms of the component masses \(m_1\) and \(m_2\).

Note

By convention, the component masses are defined such that \(m_1 \geq m_2\), so that the primary compact object in the binary (i.e., component 1), is always more massive than the secondary compact object (i.e., component 2).

In the mass diagram below, the upper diagonal region \(m_1 < m_2\) is lightly shaded in order to indicate that the definitions of four mass classes (BNS, NSBH, BBH, MassGap) are symmetric in \(m_1\) and \(m_2\).

Mass parameter space

Properties: Probabilities that the source has each of the following properties, assuming that it is not noise (e.g., assuming that it is a BNS, NSBH, BBH, or MassGap merger):

  • HasNS: The mass of one or more of the binary’s two companion compact objects is consistent with a neutron star.
  • HasRemnant: A non-zero amount of neutron star material remained outside the final remnant compact object (a necessary but not sufficient condition to produce certain kinds of electromagnetic emission such as a short GRB or a kilonova).

All of the quantities in the Classification and Properties sections are model dependent to some extent: the Classification section takes into consideration prior knowledge of astrophysical compact binary merger rates from previous LIGO/Virgo observations, and both the Classification and Properties sections depend on details of neutron star physics (e.g. maximum NS mass, equation of state). See the earlier subsection of the Data Analysis section for implementation details.

Circular Contents

The following information will be present in the human-readable GCN Circulars.

Data Quality Assessment

Circulars may contain concise descriptions of any instrument or data quality issues that may affect the significance estimates or the GW parameter inferences. Unresolved data quality issues could mean that sky localization estimates may shift after they have been mitigated, but does not mean that they will. This is to be considered as advisory information.

Sky Localization Ellipse

Generally, GW sky localizations are irregularly shaped. However, for particularly accurately localized events, the sky localization region can be well described by an ellipse. When the area of the 90% ellipse is less than 1.35 times the area of the smallest possible 90% credible region, the GCN Circular will provide a 90% containment ellipse. For details of the ellipse fitting algorithm, see ligo.skymap.postprocess.ellipse.

The ellipse is described in the format of a DS9 region string. Many tools can read DS9 region strings, including DS9, Aladin, astropy-regions, and pyregion. The region string contains the right ascension, declination, semi-major axis, semi-minor axis, position angle of the semi-minor axis). Here is an example:

icrs; ellipse(03h08m25s, -45d08m14s, 9d, 3d, 112d)

Not Included in Alerts

The alerts will not contain quantitative estimates of intrinsic properties such as masses and spins, nor contain information on the GW strain or reconstructed waveforms. After final analysis, those data products are released through the Gravitational Wave Open Science Center.

Examples

Below are some sample VOEvents to illustrate the formatting of the GCN Notices.

Preliminary
Initial
Update
Retraction
<?xml version="1.0" ?>
<voe:VOEvent xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xmlns:voe="http://www.ivoa.net/xml/VOEvent/v2.0"
xsi:schemaLocation="http://www.ivoa.net/xml/VOEvent/v2.0 http://www.ivoa.net/xml/VOEvent/VOEvent-v2.0.xsd"
 version="2.0" role="test" ivorn="ivo://gwnet/LVC#MS181101ab-1-Preliminary">
    <Who>
        <Date>2018-11-01T22:34:49</Date>
        <Author>
            <contactName>LIGO Scientific Collaboration and Virgo Collaboration</contactName>
        </Author>
    </Who>
    <What>
        <Param name="internal" dataType="int" value="0">
            <Description>Indicates whether this event should be distributed to LSC/Virgo members only</Description>
        </Param>
        <Param name="Packet_Type" dataType="int" value="150">
            <Description>The Notice Type number is assigned/used within GCN, eg type=150 is an LVC_PRELIMINARY notice</Description>
        </Param>
        <Param name="Pkt_Ser_Num" dataType="string" value="1">
            <Description>A number that increments by 1 each time a new revision is issued for this event</Description>
        </Param>
        <Param name="GraceID" dataType="string" value="MS181101ab" ucd="meta.id">
            <Description>Identifier in GraceDB</Description>
        </Param>
        <Param name="AlertType" dataType="string" value="Preliminary" ucd="meta.version">
            <Description>VOEvent alert type</Description>
        </Param>
        <Param name="HardwareInj" dataType="int" value="0" ucd="meta.number">
            <Description>Indicates that this event is a hardware injection if 1, no if 0</Description>
        </Param>
        <Param name="OpenAlert" dataType="int" value="1" ucd="meta.number">
            <Description>Indicates that this event is an open alert if 1, no if 0</Description>
        </Param>
        <Param name="EventPage" dataType="string" value="https://example.org/superevents/MS181101ab/view/" ucd="meta.ref.url">
            <Description>Web page for evolving status of this GW candidate</Description>
        </Param>
        <Param name="Instruments" dataType="string" value="H1,L1,V1" ucd="meta.code">
            <Description>List of instruments used in analysis to identify this event</Description>
        </Param>
        <Param name="FAR" dataType="float" value="9.11069936486e-14" ucd="arith.rate;stat.falsealarm" unit="Hz">
            <Description>False alarm rate for GW candidates with this strength or greater</Description>
        </Param>
        <Param name="Group" dataType="string" value="CBC" ucd="meta.code">
            <Description>Data analysis working group</Description>
        </Param>
        <Param name="Pipeline" dataType="string" value="gstlal" ucd="meta.code">
            <Description>Low-latency data analysis pipeline</Description>
        </Param>
        <Param name="Search" dataType="string" value="MDC" ucd="meta.code">
            <Description>Specific low-latency search</Description>
        </Param>
        <Group type="GW_SKYMAP" name="bayestar">
            <Param name="skymap_fits" dataType="string" value="https://emfollow.docs.ligo.org/userguide/_static/bayestar.fits.gz" ucd="meta.ref.url">
                <Description>Sky Map FITS</Description>
            </Param>
        </Group>
        <Group type="Classification">
            <Description>Source classification: binary neutron star (BNS), neutron star-black hole (NSBH), binary black hole (BBH), MassGap, or terrestrial (noise)</Description>
            <Param name="BNS" dataType="float" value="0.95" ucd="stat.probability">
                <Description>Probability that the source is a binary neutron star merger (both objects lighter than 3 solar masses)</Description>
            </Param>
            <Param name="NSBH" dataType="float" value="0.01" ucd="stat.probability">
                <Description>Probability that the source is a neutron star-black hole merger (primary heavier than 5 solar masses, secondary lighter than 3 solar masses)</Description>
            </Param>
            <Param name="BBH" dataType="float" value="0.03" ucd="stat.probability">
                <Description>Probability that the source is a binary black hole merger (both objects heavier than 5 solar masses)</Description>
            </Param>
            <Param name="MassGap" dataType="float" value="0.0" ucd="stat.probability">
                <Description>Probability that the source has at least one object between 3 and 5 solar masses</Description>
            </Param>
            <Param name="Terrestrial" dataType="float" value="0.01" ucd="stat.probability">
                <Description>Probability that the source is terrestrial (i.e., a background noise fluctuation or a glitch)</Description>
            </Param>
        </Group>
        <Group type="Properties">
            <Description>Qualitative properties of the source, conditioned on the assumption that the signal is an astrophysical compact binary merger</Description>
            <Param name="HasNS" dataType="float" value="0.95" ucd="stat.probability">
                <Description>Probability that at least one object in the binary has a mass that is less than 3 solar masses</Description>
            </Param>
            <Param name="HasRemnant" dataType="float" value="0.91" ucd="stat.probability">
                <Description>Probability that a nonzero mass was ejected outside the central remnant object</Description>
            </Param>
        </Group>
    </What>
    <WhereWhen>
        <ObsDataLocation>
            <ObservatoryLocation id="LIGO Virgo"/>
            <ObservationLocation>
                <AstroCoordSystem id="UTC-FK5-GEO"/>
                <AstroCoords coord_system_id="UTC-FK5-GEO">
                    <Time>
                        <TimeInstant>
                            <ISOTime>2018-11-01T22:22:46.654437</ISOTime>
                        </TimeInstant>
                    </Time>
                </AstroCoords>
            </ObservationLocation>
        </ObsDataLocation>
    </WhereWhen>
    <How>
        <Description>Candidate gravitational wave event identified by low-latency analysis</Description>
        <Description>H1: LIGO Hanford 4 km gravitational wave detector</Description>
        <Description>L1: LIGO Livingston 4 km gravitational wave detector</Description>
        <Description>V1: Virgo 3 km gravitational wave detector</Description>
    </How>
    <Description>Report of a candidate gravitational wave event</Description>
</voe:VOEvent>

<?xml version="1.0" ?>
<voe:VOEvent xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xmlns:voe="http://www.ivoa.net/xml/VOEvent/v2.0"
xsi:schemaLocation="http://www.ivoa.net/xml/VOEvent/v2.0 http://www.ivoa.net/xml/VOEvent/VOEvent-v2.0.xsd"
 version="2.0" role="test" ivorn="ivo://gwnet/LVC#MS181101ab-2-Initial">
    <Who>
        <Date>2018-11-01T22:36:25</Date>
        <Author>
            <contactName>LIGO Scientific Collaboration and Virgo Collaboration</contactName>
        </Author>
    </Who>
    <What>
        <Param name="internal" dataType="int" value="0">
            <Description>Indicates whether this event should be distributed to LSC/Virgo members only</Description>
        </Param>
        <Param name="Packet_Type" dataType="int" value="151">
            <Description>The Notice Type number is assigned/used within GCN, eg type=151 is an LVC_INITIAL notice</Description>
        </Param>
        <Param name="Pkt_Ser_Num" dataType="string" value="2">
            <Description>A number that increments by 1 each time a new revision is issued for this event</Description>
        </Param>
        <Param name="GraceID" dataType="string" value="MS181101ab" ucd="meta.id">
            <Description>Identifier in GraceDB</Description>
        </Param>
        <Param name="AlertType" dataType="string" value="Initial" ucd="meta.version">
            <Description>VOEvent alert type</Description>
        </Param>
        <Param name="HardwareInj" dataType="int" value="0" ucd="meta.number">
            <Description>Indicates that this event is a hardware injection if 1, no if 0</Description>
        </Param>
        <Param name="OpenAlert" dataType="int" value="1" ucd="meta.number">
            <Description>Indicates that this event is an open alert if 1, no if 0</Description>
        </Param>
        <Param name="EventPage" dataType="string" value="https://example.org/superevents/MS181101ab/view/" ucd="meta.ref.url">
            <Description>Web page for evolving status of this GW candidate</Description>
        </Param>
        <Param name="Instruments" dataType="string" value="H1,L1,V1" ucd="meta.code">
            <Description>List of instruments used in analysis to identify this event</Description>
        </Param>
        <Param name="FAR" dataType="float" value="9.11069936486e-14" ucd="arith.rate;stat.falsealarm" unit="Hz">
            <Description>False alarm rate for GW candidates with this strength or greater</Description>
        </Param>
        <Param name="Group" dataType="string" value="CBC" ucd="meta.code">
            <Description>Data analysis working group</Description>
        </Param>
        <Param name="Pipeline" dataType="string" value="gstlal" ucd="meta.code">
            <Description>Low-latency data analysis pipeline</Description>
        </Param>
        <Param name="Search" dataType="string" value="MDC" ucd="meta.code">
            <Description>Specific low-latency search</Description>
        </Param>
        <Group type="GW_SKYMAP" name="bayestar">
            <Param name="skymap_fits" dataType="string" value="https://emfollow.docs.ligo.org/userguide/_static/bayestar.fits.gz" ucd="meta.ref.url">
                <Description>Sky Map FITS</Description>
            </Param>
        </Group>
        <Group type="Classification">
            <Description>Source classification: binary neutron star (BNS), neutron star-black hole (NSBH), binary black hole (BBH), MassGap, or terrestrial (noise)</Description>
            <Param name="BNS" dataType="float" value="0.95" ucd="stat.probability">
                <Description>Probability that the source is a binary neutron star merger (both objects lighter than 3 solar masses)</Description>
            </Param>
            <Param name="NSBH" dataType="float" value="0.01" ucd="stat.probability">
                <Description>Probability that the source is a neutron star-black hole merger (primary heavier than 5 solar masses, secondary lighter than 3 solar masses)</Description>
            </Param>
            <Param name="BBH" dataType="float" value="0.03" ucd="stat.probability">
                <Description>Probability that the source is a binary black hole merger (both objects heavier than 5 solar masses)</Description>
            </Param>
            <Param name="MassGap" dataType="float" value="0.0" ucd="stat.probability">
                <Description>Probability that the source has at least one object between 3 and 5 solar masses</Description>
            </Param>
            <Param name="Terrestrial" dataType="float" value="0.01" ucd="stat.probability">
                <Description>Probability that the source is terrestrial (i.e., a background noise fluctuation or a glitch)</Description>
            </Param>
        </Group>
        <Group type="Properties">
            <Description>Qualitative properties of the source, conditioned on the assumption that the signal is an astrophysical compact binary merger</Description>
            <Param name="HasNS" dataType="float" value="0.95" ucd="stat.probability">
                <Description>Probability that at least one object in the binary has a mass that is less than 3 solar masses</Description>
            </Param>
            <Param name="HasRemnant" dataType="float" value="0.91" ucd="stat.probability">
                <Description>Probability that a nonzero mass was ejected outside the central remnant object</Description>
            </Param>
        </Group>
    </What>
    <WhereWhen>
        <ObsDataLocation>
            <ObservatoryLocation id="LIGO Virgo"/>
            <ObservationLocation>
                <AstroCoordSystem id="UTC-FK5-GEO"/>
                <AstroCoords coord_system_id="UTC-FK5-GEO">
                    <Time>
                        <TimeInstant>
                            <ISOTime>2018-11-01T22:22:46.654437</ISOTime>
                        </TimeInstant>
                    </Time>
                </AstroCoords>
            </ObservationLocation>
        </ObsDataLocation>
    </WhereWhen>
    <How>
        <Description>Candidate gravitational wave event identified by low-latency analysis</Description>
        <Description>H1: LIGO Hanford 4 km gravitational wave detector</Description>
        <Description>L1: LIGO Livingston 4 km gravitational wave detector</Description>
        <Description>V1: Virgo 3 km gravitational wave detector</Description>
    </How>
    <Citations>
        <EventIVORN cite="supersedes">ivo://gwnet/LVC#MS181101ab-1-Preliminary</EventIVORN>
        <Description>Updated localization is now available</Description>
    </Citations>
    <Description>Report of a candidate gravitational wave event</Description>
</voe:VOEvent>

<?xml version="1.0" ?>
<voe:VOEvent xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xmlns:voe="http://www.ivoa.net/xml/VOEvent/v2.0"
xsi:schemaLocation="http://www.ivoa.net/xml/VOEvent/v2.0 http://www.ivoa.net/xml/VOEvent/VOEvent-v2.0.xsd"
 version="2.0" role="test" ivorn="ivo://gwnet/LVC#MS181101ab-3-Update">
    <Who>
        <Date>2018-11-01T22:36:25</Date>
        <Author>
            <contactName>LIGO Scientific Collaboration and Virgo Collaboration</contactName>
        </Author>
    </Who>
    <What>
        <Param name="internal" dataType="int" value="0">
            <Description>Indicates whether this event should be distributed to LSC/Virgo members only</Description>
        </Param>
        <Param name="Packet_Type" dataType="int" value="152">
            <Description>The Notice Type number is assigned/used within GCN, eg type=152 is an LVC_UPDATE notice</Description>
        </Param>
        <Param name="Pkt_Ser_Num" dataType="string" value="3">
            <Description>A number that increments by 1 each time a new revision is issued for this event</Description>
        </Param>
        <Param name="GraceID" dataType="string" value="MS181101ab" ucd="meta.id">
            <Description>Identifier in GraceDB</Description>
        </Param>
        <Param name="AlertType" dataType="string" value="Update" ucd="meta.version">
            <Description>VOEvent alert type</Description>
        </Param>
        <Param name="HardwareInj" dataType="int" value="0" ucd="meta.number">
            <Description>Indicates that this event is a hardware injection if 1, no if 0</Description>
        </Param>
        <Param name="OpenAlert" dataType="int" value="1" ucd="meta.number">
            <Description>Indicates that this event is an open alert if 1, no if 0</Description>
        </Param>
        <Param name="EventPage" dataType="string" value="https://example.org/superevents/MS181101ab/view/" ucd="meta.ref.url">
            <Description>Web page for evolving status of this GW candidate</Description>
        </Param>
        <Param name="Instruments" dataType="string" value="H1,L1,V1" ucd="meta.code">
            <Description>List of instruments used in analysis to identify this event</Description>
        </Param>
        <Param name="FAR" dataType="float" value="9.11069936486e-14" ucd="arith.rate;stat.falsealarm" unit="Hz">
            <Description>False alarm rate for GW candidates with this strength or greater</Description>
        </Param>
        <Param name="Group" dataType="string" value="CBC" ucd="meta.code">
            <Description>Data analysis working group</Description>
        </Param>
        <Param name="Pipeline" dataType="string" value="gstlal" ucd="meta.code">
            <Description>Low-latency data analysis pipeline</Description>
        </Param>
        <Param name="Search" dataType="string" value="MDC" ucd="meta.code">
            <Description>Specific low-latency search</Description>
        </Param>
        <Group type="GW_SKYMAP" name="bayestar">
            <Param name="skymap_fits" dataType="string" value="https://emfollow.docs.ligo.org/userguide/_static/bayestar.fits.gz" ucd="meta.ref.url">
                <Description>Sky Map FITS</Description>
            </Param>
        </Group>
        <Group type="Classification">
            <Description>Source classification: binary neutron star (BNS), neutron star-black hole (NSBH), binary black hole (BBH), MassGap, or terrestrial (noise)</Description>
            <Param name="BNS" dataType="float" value="0.95" ucd="stat.probability">
                <Description>Probability that the source is a binary neutron star merger (both objects lighter than 3 solar masses)</Description>
            </Param>
            <Param name="NSBH" dataType="float" value="0.01" ucd="stat.probability">
                <Description>Probability that the source is a neutron star-black hole merger (primary heavier than 5 solar masses, secondary lighter than 3 solar masses)</Description>
            </Param>
            <Param name="BBH" dataType="float" value="0.03" ucd="stat.probability">
                <Description>Probability that the source is a binary black hole merger (both objects heavier than 5 solar masses)</Description>
            </Param>
            <Param name="MassGap" dataType="float" value="0.0" ucd="stat.probability">
                <Description>Probability that the source has at least one object between 3 and 5 solar masses</Description>
            </Param>
            <Param name="Terrestrial" dataType="float" value="0.01" ucd="stat.probability">
                <Description>Probability that the source is terrestrial (i.e., a background noise fluctuation or a glitch)</Description>
            </Param>
        </Group>
        <Group type="Properties">
            <Description>Qualitative properties of the source, conditioned on the assumption that the signal is an astrophysical compact binary merger</Description>
            <Param name="HasNS" dataType="float" value="0.95" ucd="stat.probability">
                <Description>Probability that at least one object in the binary has a mass that is less than 3 solar masses</Description>
            </Param>
            <Param name="HasRemnant" dataType="float" value="0.91" ucd="stat.probability">
                <Description>Probability that a nonzero mass was ejected outside the central remnant object</Description>
            </Param>
        </Group>
    </What>
    <WhereWhen>
        <ObsDataLocation>
            <ObservatoryLocation id="LIGO Virgo"/>
            <ObservationLocation>
                <AstroCoordSystem id="UTC-FK5-GEO"/>
                <AstroCoords coord_system_id="UTC-FK5-GEO">
                    <Time>
                        <TimeInstant>
                            <ISOTime>2018-11-01T22:22:46.654437</ISOTime>
                        </TimeInstant>
                    </Time>
                </AstroCoords>
            </ObservationLocation>
        </ObsDataLocation>
    </WhereWhen>
    <How>
        <Description>Candidate gravitational wave event identified by low-latency analysis</Description>
        <Description>H1: LIGO Hanford 4 km gravitational wave detector</Description>
        <Description>L1: LIGO Livingston 4 km gravitational wave detector</Description>
        <Description>V1: Virgo 3 km gravitational wave detector</Description>
    </How>
    <Citations>
        <EventIVORN cite="supersedes">ivo://gwnet/LVC#MS181101ab-1-Preliminary</EventIVORN>
        <EventIVORN cite="supersedes">ivo://gwnet/LVC#MS181101ab-2-Initial</EventIVORN>
        <Description>Updated localization is now available</Description>
    </Citations>
    <Description>Report of a candidate gravitational wave event</Description>
</voe:VOEvent>

<?xml version="1.0" ?>
<voe:VOEvent xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xmlns:voe="http://www.ivoa.net/xml/VOEvent/v2.0"
xsi:schemaLocation="http://www.ivoa.net/xml/VOEvent/v2.0 http://www.ivoa.net/xml/VOEvent/VOEvent-v2.0.xsd"
 version="2.0" role="test" ivorn="ivo://gwnet/LVC#MS181101ab-4-Retraction">
    <Who>
        <Date>2018-11-01T23:36:23</Date>
        <Author>
            <contactName>LIGO Scientific Collaboration and Virgo Collaboration</contactName>
        </Author>
    </Who>
    <What>
        <Param name="internal" dataType="int" value="0">
            <Description>Indicates whether this event should be distributed to LSC/Virgo members only</Description>
        </Param>
        <Param name="Packet_Type" dataType="int" value="164">
            <Description>The Notice Type number is assigned/used within GCN, eg type=164 is an LVC_RETRACTION notice</Description>
        </Param>
        <Param name="Pkt_Ser_Num" dataType="string" value="4">
            <Description>A number that increments by 1 each time a new revision is issued for this event</Description>
        </Param>
        <Param name="GraceID" dataType="string" value="MS181101ab" ucd="meta.id">
            <Description>Identifier in GraceDB</Description>
        </Param>
        <Param name="AlertType" dataType="string" value="Retraction" ucd="meta.version">
            <Description>VOEvent alert type</Description>
        </Param>
        <Param name="HardwareInj" dataType="int" value="0" ucd="meta.number">
            <Description>Indicates that this event is a hardware injection if 1, no if 0</Description>
        </Param>
        <Param name="OpenAlert" dataType="int" value="1" ucd="meta.number">
            <Description>Indicates that this event is an open alert if 1, no if 0</Description>
        </Param>
        <Param name="EventPage" dataType="string" value="https://example.org/superevents/MS181101ab/view/" ucd="meta.ref.url">
            <Description>Web page for evolving status of this GW candidate</Description>
        </Param>
    </What>
    <WhereWhen>
        <ObsDataLocation>
            <ObservatoryLocation id="LIGO Virgo"/>
            <ObservationLocation>
                <AstroCoordSystem id="UTC-FK5-GEO"/>
                <AstroCoords coord_system_id="UTC-FK5-GEO">
                    <Time>
                        <TimeInstant>
                            <ISOTime>2018-11-01T22:22:46.654437</ISOTime>
                        </TimeInstant>
                    </Time>
                </AstroCoords>
            </ObservationLocation>
        </ObsDataLocation>
    </WhereWhen>
    <Citations>
        <EventIVORN cite="retraction">ivo://gwnet/LVC#MS181101ab-1-Preliminary</EventIVORN>
        <EventIVORN cite="retraction">ivo://gwnet/LVC#MS181101ab-2-Initial</EventIVORN>
        <EventIVORN cite="retraction">ivo://gwnet/LVC#MS181101ab-3-Update</EventIVORN>
        <Description>Determined to not be a viable GW event candidate</Description>
    </Citations>
</voe:VOEvent>

Sample Code

This section provides Python sample code for receiving and interacting with GCN Notices. GCN Notices are available over several different protocols and in several different formats. LIGO/Virgo strongly recommends using the VOEvent Transport Protocol (VTP) to receive notices in VOEvent XML format because it is anonymous, configuration-free, and easy to parse.

This tutorial will walk you through writing a Python script to receive and process the example LIGO/Virgo GCN notices that are sent every hour. The tutorial is broken into the following subsections:

Prerequisites

In order to run this sample code, you will need Python >= 3.5 on a Unix-like operating system (Linux or macOS) and a few third-party Python packages:

If you are on a Mac and use the MacPorts package manager, you can install all of the above with the following command:

$ sudo port install py37-gcn py37-healpy

Otherwise, the fastest way to install the dependencies is with pip, a package manager that comes with most Python distributions. To install these packages with pip, run the following command:

$ pip install pygcn healpy

Another option is the Anaconda Python distribution. To install these packages using Anaconda, first install conda and then run the following commands:

$ conda config --add channels conda-forge
$ conda install pygcn healpy

Receiving GCNs

Next, we’ll write a GCN handler function that we want PyGCN to call every time it receives a GCN notice. We decorate the handler with @gcn.handlers.include_notice_types to specify that we only want to process certain GCN notice types (LVC_PRELIMINARY, LVC_INITIAL, and LVC_UDPATE).

Events come in two very general flavors: ‘CBC’ for compact binary coalescence candidates detected by matched filtering, and ‘Burst’ for candidates detected by model-independent methods. Your handler can take different actions based on this. The example below will handle only ‘CBC’ events.

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.

Note

Observe in the example below that we do not have to explicitly download the FITS file because the hp.read_map() function works with either URLs or filenames. However, you could download and save the FITS file in order to save it locally using astropy.utils.data.download_file(), requests.get(), urllib.request.urlopen(), or even curl.

The following basic handler function will parse out the URL of the FITS file, download it, and extract the probability sky map:

import gcn
import healpy as hp

# Function to call every time a GCN is received.
# Run only for notices of type
# LVC_PRELIMINARY, LVC_INITIAL, or LVC_UPDATE.
@gcn.handlers.include_notice_types(
    gcn.notice_types.LVC_PRELIMINARY,
    gcn.notice_types.LVC_INITIAL,
    gcn.notice_types.LVC_UPDATE,
    gcn.notice_types.LVC_RETRACTION)
def process_gcn(payload, root):
    # Respond only to 'test' events.
    # VERY IMPORTANT! Replace with the following code
    # to respond to only real 'observation' events.
    # if root.attrib['role'] != 'observation':
    #    return
    if root.attrib['role'] != 'test':
        return

    # Read all of the VOEvent parameters from the "What" section.
    params = {elem.attrib['name']:
              elem.attrib['value']
              for elem in root.iterfind('.//Param')}

    # Respond only to 'CBC' events. Change 'CBC' to "Burst'
    # to respond to only unmodeled burst events.
    if params['Group'] != 'CBC':
        return

    # Print all parameters.
    for key, value in params.items():
        print(key, '=', value)

    if 'skymap_fits' in params:
        # Read the HEALPix sky map and the FITS header.
        skymap, header = hp.read_map(params['skymap_fits'],
                                     h=True, verbose=False)
        header = dict(header)

        # Print some values from the FITS header.
        print('Distance =', header['DISTMEAN'], '+/-', header['DISTSTD'])
Listen for GCNs

Now, we will start the VOEvent client to listen for GCNs using the gcn.listen function. By default, this will connect to the anonymous, public GCN server. You just need to tell gcn.listen what function to call whenever it receives an GCN; in this example, that is the process_gcn handler that we defined above.

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

When you run this script you should receive a sample LIGO/Virgo GCN Notice every hour. For each event received it will print output that looks like what is shown in the Offline Testing example below.

Note

gcn.listen will try to automatically reconnect if the network connection is ever broken.

Offline Testing

Sometimes it is convenient to be able to explicitly call the GCN handler with a sample input, rather than waiting for the next broadcast of a sample alert. You can download the example GCN notices from this documentation and pass it into your GCN handler at any time. First, download the sample GCN notice:

$ curl -O https://emfollow.docs.ligo.org/userguide/_static/MS181101ab-1-Preliminary.xml

Then you can manually invoke your GCN handler using this Python code:

import lxml.etree
payload = open('MS181101ab-1-Preliminary.xml', 'rb').read()
root = lxml.etree.fromstring(payload)
process_gcn(payload, root)

Upon running this, you should see:

internal = 0
Packet_Type = 150
Pkt_Ser_Num = 1
GraceID = MS181101ab
AlertType = Preliminary
HardwareInj = 0
OpenAlert = 1
EventPage = https://example.org/superevents/MS181101ab/view/
Instruments = H1,L1,V1
FAR = 9.11069936486e-14
Group = CBC
Pipeline = gstlal
Search = MDC
skymap_fits = https://emfollow.docs.ligo.org/userguide/_static/bayestar.fits.gz
BNS = 0.95
NSBH = 0.01
BBH = 0.03
MassGap = 0.0
Terrestrial = 0.01
HasNS = 0.95
HasRemnant = 0.91
Distance = 39.76999609489013 +/- 8.308435058808886

Working with Sky Maps

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 hierarchical, progressively refined, all-sky survey images (HiPS). Second, the value stored at each pixel is the probability that the gravitational-wave source is within that pixel.

HEALPix projection

Let’s download an example FITS file with curl:

$ curl -O https://emfollow.docs.ligo.org/userguide/_static/bayestar.fits.gz

We can look at the metadata inside the FITS file by printing its header with tools like funhead from Funtools, imhead from WCSTools, or fitsheader from Astropy:

$ fitsheader bayestar.fits.gz
# HDU 0 in bayestar.fits.gz:
SIMPLE  =                    T / conforms to FITS standard
BITPIX  =                    8 / array data type
NAXIS   =                    0 / number of array dimensions
EXTEND  =                    T

# HDU 1 in bayestar.fits.gz:
XTENSION= 'BINTABLE'           / binary table extension
BITPIX  =                    8 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                   32 / length of dimension 1
NAXIS2  =             50331648 / length of dimension 2
PCOUNT  =                    0 / number of group parameters
GCOUNT  =                    1 / number of groups
TFIELDS =                    4 / number of table fields
TTYPE1  = 'PROB    '
TFORM1  = 'D       '
TUNIT1  = 'pix-1   '
TTYPE2  = 'DISTMU  '
TFORM2  = 'D       '
TUNIT2  = 'Mpc     '
TTYPE3  = 'DISTSIGMA'
TFORM3  = 'D       '
TUNIT3  = 'Mpc     '
TTYPE4  = 'DISTNORM'
TFORM4  = 'D       '
TUNIT4  = 'Mpc-2   '
PIXTYPE = 'HEALPIX '           / HEALPIX pixelisation
ORDERING= 'NESTED  '           / Pixel ordering scheme: RING, NESTED, or NUNIQ
COORDSYS= 'C       '           / Ecliptic, Galactic or Celestial (equatorial)
NSIDE   =                 2048 / Resolution parameter of HEALPIX
INDXSCHM= 'IMPLICIT'           / Indexing: IMPLICIT or EXPLICIT
OBJECT  = 'MS181101ab'         / Unique identifier for this event
REFERENC= 'https://example.org/superevents/MS181101ab/view/' / URL of this event
INSTRUME= 'H1,L1,V1'           / Instruments that triggered this event
DATE-OBS= '2018-11-01T22:22:46.654437' / UTC date of the observation
MJD-OBS =    58423.93248442613 / modified Julian date of the observation
DATE    = '2018-11-01T22:34:49.000000' / UTC date of file creation
CREATOR = 'BAYESTAR'           / Program that created this file
ORIGIN  = 'LIGO/Virgo'         / Organization responsible for this FITS file
RUNTIME =     3.24746292643249 / Runtime in seconds of the CREATOR program
DISTMEAN=    39.76999609489013 / Posterior mean distance (Mpc)
DISTSTD =    8.308435058808886 / Posterior standard deviation of distance (Mpc)
LOGBCI  =    13.64819688928804 / Log Bayes factor: coherent vs. incoherent
LOGBSN  =    261.0250944470225 / Log Bayes factor: signal vs. noise
VCSVERS = 'ligo.skymap 0.1.8'  / Software version
VCSREV  = 'becb07110491d799b753858845b5c24c82705404' / Software revision (Git)
DATE-BLD= '2019-07-25T22:36:58' / Software build date
HISTORY
HISTORY Generated by calling the following Python function:
HISTORY ligo.skymap.bayestar.localize(event=..., waveform='o2-uberbank', f_low=3
HISTORY 0, min_inclination=0.0, max_inclination=1.5707963267948966, min_distance
HISTORY =None, max_distance=None, prior_distance_power=2, cosmology=False, mcmc=
HISTORY False, chain_dump=None, enable_snr_series=True, f_high_truncate=0.95)
HISTORY
HISTORY This was the command line that started the program:
HISTORY bayestar-localize-lvalert -N G298107 -o bayestar.multiorder.fits

There are several useful pieces of information here:

  • COORDSYS=C, telling you that the HEALPix projection is in the Celestial (equatorial, ICRS) frame, as all LIGO/Virgo probability sky maps will be.
  • OBJECT, the unique LIGO/Virgo identifier for the event.
  • REFERENC, a link to the candidate page in GraceDB.
  • INSTRUME, a list of gravitational-wave sites that triggered on the event: H1 for LIGO Hanford, L1 for LIGO Livingston, and V1 for Virgo.
  • DATE-OBS, the UTC time of the event. In the case of a compact binary coalescence candidate, this is the time that the signal from the merger passed through the geocenter.
  • MJD-OBS, same as DATE-OBS, but given as a modified Julian day.

You can view the sky map in many common FITS image viewers such as Aladin:

Aladin screenshot

or DS9 (although DS9 shows HEALPix sky maps in an unusual orientation; see Figure 4 of [1] for more information.

DS9 screenshot

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, and Java. However, since this is a Python tutorial, we are going to demonstrate how to manipulate HEALPix maps with the official Python library, Healpy.

Reading Sky Maps

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:

$ curl -O https://emfollow.docs.ligo.org/userguide/_static/bayestar.fits.gz

Next, we need to read in the file in Python with Healpy:

>>> import healpy as hp
>>> import numpy as np
>>> hpx = hp.read_map('bayestar.fits.gz')
NSIDE = 2048
ORDERING = NESTED in fits file
INDXSCHM = IMPLICIT
Ordering converted to RING

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:

>>> hpx, header = hp.read_map('bayestar.fits.gz', h=True, verbose=False)
Manipulating HEALPix Coordinates

The image data is a 1D array of values:

>>> hpx
array([2.70726059e-66, 1.27374324e-66, 2.62611513e-67, ...,
       2.04700874e-40, 1.05781210e-35, 4.44174764e-31])

Healpy has several useful plotting routines including hp.mollview for plotting a Mollweide-projection all-sky map:

>>> hp.mollview(hpx)
_images/skymaps-2.svg

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:

>>> npix = len(hpx)
>>> sky_area = 4 * 180**2 / np.pi
>>> sky_area / npix
0.0008196227004015301

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

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

>>> nside = hp.npix2nside(npix)
>>> nside
2048

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

>>> ipix = 123
>>> theta, phi = hp.pix2ang(nside, ipix)
>>> ra = np.rad2deg(phi)
>>> dec = np.rad2deg(0.5 * np.pi - theta)
>>> ra, dec
(129.375, 89.81725848475484)

Let’s find which pixel contains the point RA=194.95, Dec=27.98.

>>> ra = 194.95
>>> dec = 27.98
>>> theta = 0.5 * np.pi - np.deg2rad(dec)
>>> phi = np.deg2rad(ra)
>>> ipix = hp.ang2pix(nside, theta, phi)
>>> ipix
13361492
Test if a Sky Location is in the 90% Credible Region

You can easily test if a given sky position is in the 90% credible region. Let’s continue using the sky position from the previous example, for which we have already determined the pixel index.

Use the following simple algorithm to construct a map that gives the credible level of each pixel:

  1. Sort the pixels by descending probability density.
  2. Cumulatively sum the pixels’ probability.
  3. Return the pixels to their original order.

In Python, you can use this simple recipe:

>>> i = np.flipud(np.argsort(hpx))
>>> sorted_credible_levels = np.cumsum(hpx[i])
>>> credible_levels = np.empty_like(sorted_credible_levels)
>>> credible_levels[i] = sorted_credible_levels
>>> credible_levels
array([1., 1., 1., ..., 1., 1., 1.])

Note

Observe that the values in the resulting credible level map vary inversely with probability density: the most probable pixel is assigned to the credible level 0.0, and the least likely pixel is assigned the credible level 1.0.

Tip

This recipe is implemented in the package ligo.skymap as the function find_greedy_credible_levels():

>>> from ligo.skymap.postprocess import find_greedy_credible_levels
>>> credible_levels = find_greedy_credible_levels(hpx)
>>> credible_levels
array([1., 1., 1., ..., 1., 1., 1.])

To check if the pixel that we identified in the previous section is within the 90% credible level, simply test if the value of the credible level map is less than or equal to 0.9 at that pixel:

>>> credible_levels[ipix]
0.9999999999947833
>>> credible_levels[ipix] <= 0.9
False

The credible level map has a value greater than 0.9 at that sky location, therefore the sky location is outside the 90% credible region.

Find the Area of the 90% Credible Region

Since we just found the credible level map, it’s easy to compute the 90% credible area by counting the number of pixels inside the 90% credible region and multiplying by the area per pixel.

In the Python expression below, note that (credible_levels <= 0.9) evaluates to a binary array; when it is summed over, true values are treated as 1 and false values are treated as 0.

>>> np.sum(credible_levels <= 0.9) * hp.nside2pixarea(nside, degrees=True)
30.979279207076633
Most Probable Sky Location

Let’s find the highest probability pixel.

>>> ipix_max = np.argmax(hpx)

What is the probability density per square degree at that position?

>>> hpx[ipix_max] / hp.nside2pixarea(nside, degrees=True)
0.0782516470191411

Where is the highest probability pixel on the sky? Use hp.pix2ang.

>>> theta, phi = hp.pix2ang(nside, ipix_max)
>>> ra = np.rad2deg(phi)
>>> dec = np.rad2deg(0.5 * np.pi - theta)
>>> ra, dec
(194.30419921875, -17.856895095545475)
Integrated Probability in a Circle

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 hp.query_disc. Note that this function 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.

First, we define the RA, Dec, and radius of circle in degrees:

>>> ra = 213.22
>>> dec = -37.45
>>> radius = 3.1

Then we convert to spherical polar coordinates and radius of circle in radians:

>>> theta = 0.5 * np.pi - np.deg2rad(dec)
>>> phi = np.deg2rad(ra)
>>> radius = np.deg2rad(radius)

Then we calculate the Cartesian coordinates of the center of circle:

>>> xyz = hp.ang2vec(theta, phi)

We call hp.query_disc, which returns an array of the indices of the pixels that are inside the circle:

>>> ipix_disc = hp.query_disc(nside, xyz, radius)

Finally, we sum the probability in all of the matching pixels:

>>> hpx[ipix_disc].sum()
3.655661941088471e-10
Integrated Probability in a Polygon

Similarly, we can use the hp.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.

>>> xyz = [[-0.69601758, -0.41315628, -0.58724902],
...        [-0.68590811, -0.40679797, -0.60336181],
...        [-0.69106913, -0.39820114, -0.60320752],
...        [-0.7011786 , -0.40455945, -0.58709473]]
>>> ipix_poly = hp.query_polygon(nside, xyz)
>>> hpx[ipix_poly].sum()
1.128695302404769e-12

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 hp.ud_grade for upsampling or downsampling a sky map and hp.get_interp_val for performing bilinear interpolation between pixels. See the Healpy tutorial for other useful operations.

[1]Calabretta, M. R., & Roukema, B. F. 2007, Mon. Notices Royal Astron. Soc., 381, 865. doi:10.1111/j.1365-2966.2007.12297.x

Multi-Order Sky Maps (For Advanced Users)

For most events, LIGO/Virgo distributes both the standard HEALPix format with the file extension .fits.gz, as well as an experimental multi-resolution HEALPix format, distinguished by the file extension .multiorder.fits.

What Problem Do Multi-Resolution Sky Maps Solve?

The multi-resolution format has been introduced as a forward-looking solution to deal with computational challenges related to highly accurate localizations. We are getting better at pinpointing gravitational-wave sources as more detectors come online and existing detectors become more sensitive. Unfortunately, as position accuracy improves, the size of the standard HEALPix sky maps will blow up. This started being a minor inconvenience in O2 with GW170817. It will get slowly worse as we approach design sensitivity. It’s already a major pain if you are studying future detector networks with simulations.

It is worth reviewing why LIGO/Virgo has adopted HEALPix rather than a more commonplace image format for sky maps in the first place. Gravitational-wave localizations are distinguished from many other kinds of astronomical image data by the following features:

  • The probability regions can subtend large angles.
  • They can wrap around the whole sky.
  • They can have multiple widely separated modes.
  • They can have irregular shapes or interference-like fringes.

As a consequence of these features, it is difficult to pick a good partial-sky projection (e.g. gnomonic, orthographic) in the general case. Traditional all-sky projections have wild variations in pixel size (e.g. plate carée) or shape (e.g. Mollweide, Aitoff) and as well as seams at the projection boundaries. HEALPix was already well-established for specialized uses in astronomy that cannot tolerate such projection artifacts (e.g. cosmic microwave background data sets and full-sky mosaics from optical surveys).

The natural sky resolution varies from one gravitational-wave event to another depending on its SNR and the number of detectors. During early Advanced LIGO and Virgo, HEALPix resolutions of nside=512 (\(6.9'\) per pixel) to nside=2048 (\(1.7''\) per pixel) were adequate. However, as existing gravitational-wave detectors improve in sensitivity and additional detectors come online, finer resolutions will be required.

Fortunately, the increased resolution will come at little to no computational cost for actually producing localizations because most LIGO/Virgo parameter estimation analyses use a simple multi-resolution adaptive mesh refinement scheme that limits them to sampling the sky at only about 20k points.

HEALPix mesh refinement scheme

An illustration of the adaptive mesh refinement scheme, reproduced from [1].

When these multi-resolution meshes are flattened to a single HEALPix resolution, all but the finest nodes in the mesh become long sequences of repeated pixel values. High resolution also does not cost much in terms of disk space because gzip compression can store the long runs of repeated pixel values efficiently. The diagram below illustrates a multi-resolution structure that is fairly typical of gravitational-wave localizations from O1 and O2.

An example HEALPix multi-resolution mesh

An example multi-resolution mesh from a typical two-detector (LHO and LLO) localization produced with BAYESTAR. Reproduced from [1]. The scale bar at bottom right has a length of 10°.

However, the resolution does come at a significant cost in the time it takes to decompress and read the FITS files (already up to tens of seconds for GW170817) and in terms of memory (up to several gigabytes). The time and memory will worsen as localization accuracy improves.

The multi-resolution format is immune to these issues because it is a direct representation of the adaptive mesh produced by the LIGO/Virgo localization algorithms.

The UNIQ Indexing Scheme

Recall from before that three pieces of information are required to specify a HEALPix tile: nside to specify the resolution, ipix to identify a sky position at that resolution, and the indexing scheme.

HEALPix has a couple different indexing schemes. In the RING scheme, indices advance west to east and then north to south. In the NESTED scheme, indices encode the hierarchy of parent pixels in successively lower resolutions. The image below illustrates these two indexing schemes.

HEALPix RING and NESTED indexing schemes

The RING and NESTED indexing schemes of HEALPix. Reproduced from The image below, reproduced from [2].

There is a third HEALPix indexing scheme called UNIQ. The UNIQ indexing scheme is special because it encodes both the resolution and the sky position in a single integer. It assigns a single unique integer to every HEALPix tile at every resolution. If ipix is the pixel index in the NESTED ordering, then the unique pixel index uniq is:

\[\mathit{uniq} = \mathit{ipix} + 4 \, \mathit{nside}^2.\]

The inverse is:

\[\begin{split}\mathit{nside} = 2^{\lfloor \log_2(\mathit{uniq}/4)/2 \rfloor} \\ \mathit{ipix} = \mathit{uniq} - 4 \, \mathit{nside}^2.\end{split}\]
FITS Format for Multi-Order Sky Maps

The FITS format for LIGO/Virgo multi-resolution sky maps uses the UNIQ indexing scheme and is a superset of the FITS serialization for Multi-Order Coverage (MOC) maps specified by IVOA [4] as part of the Hierarchical Progressive Survey (HiPS) capability [3], notably used by Aladin for storing and display all-sky image mosaics.

Let’s download an example multi-order FITS file with curl:

$ curl -O https://emfollow.docs.ligo.org/userguide/_static/bayestar.multiorder.fits

Let’s look at the FITS header:

$ fitsheader bayestar.multiorder.fits
# HDU 0 in bayestar.multiorder.fits:
SIMPLE  =                    T / conforms to FITS standard
BITPIX  =                    8 / array data type
NAXIS   =                    0 / number of array dimensions
EXTEND  =                    T

# HDU 1 in bayestar.multiorder.fits:
XTENSION= 'BINTABLE'           / binary table extension
BITPIX  =                    8 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                   40 / length of dimension 1
NAXIS2  =                19200 / length of dimension 2
PCOUNT  =                    0 / number of group parameters
GCOUNT  =                    1 / number of groups
TFIELDS =                    5 / number of table fields
TTYPE1  = 'UNIQ    '
TFORM1  = 'K       '
TTYPE2  = 'PROBDENSITY'
TFORM2  = 'D       '
TUNIT2  = 'sr-1    '
TTYPE3  = 'DISTMU  '
TFORM3  = 'D       '
TUNIT3  = 'Mpc     '
TTYPE4  = 'DISTSIGMA'
TFORM4  = 'D       '
TUNIT4  = 'Mpc     '
TTYPE5  = 'DISTNORM'
TFORM5  = 'D       '
TUNIT5  = 'Mpc-2   '
PIXTYPE = 'HEALPIX '           / HEALPIX pixelisation
ORDERING= 'NUNIQ   '           / Pixel ordering scheme: RING, NESTED, or NUNIQ
COORDSYS= 'C       '           / Ecliptic, Galactic or Celestial (equatorial)
MOCORDER=                   11 / MOC resolution (best order)
INDXSCHM= 'EXPLICIT'           / Indexing: IMPLICIT or EXPLICIT
OBJECT  = 'MS181101ab'         / Unique identifier for this event
REFERENC= 'https://example.org/superevents/MS181101ab/view/' / URL of this event
INSTRUME= 'H1,L1,V1'           / Instruments that triggered this event
DATE-OBS= '2018-11-01T22:22:46.654437' / UTC date of the observation
MJD-OBS =    58423.93248442635 / modified Julian date of the observation
DATE    = '2018-11-01T22:34:49.000000' / UTC date of file creation
CREATOR = 'BAYESTAR'           / Program that created this file
ORIGIN  = 'LIGO/Virgo'         / Organization responsible for this FITS file
RUNTIME =     3.24746292643249 / Runtime in seconds of the CREATOR program
DISTMEAN=    39.76999609489013 / Posterior mean distance (Mpc)
DISTSTD =    8.308435058808886 / Posterior standard deviation of distance (Mpc)
LOGBCI  =    13.64819688928804 / Log Bayes factor: coherent vs. incoherent
LOGBSN  =    261.0250944470225 / Log Bayes factor: signal vs. noise
VCSVERS = 'ligo.skymap 0.1.8'  / Software version
VCSREV  = 'becb07110491d799b753858845b5c24c82705404' / Software revision (Git)
DATE-BLD= '2019-07-25T22:36:58' / Software build date
HISTORY
HISTORY Generated by calling the following Python function:
HISTORY ligo.skymap.bayestar.localize(event=..., waveform='o2-uberbank', f_low=3
HISTORY 0, min_inclination=0.0, max_inclination=1.5707963267948966, min_distance
HISTORY =None, max_distance=None, prior_distance_power=2, cosmology=False, mcmc=
HISTORY False, chain_dump=None, enable_snr_series=True, f_high_truncate=0.95)
HISTORY
HISTORY This was the command line that started the program:
HISTORY bayestar-localize-lvalert -N G298107 -o bayestar.multiorder.fits

This should look very similar to the FITS header for the standard HEALPix file from the previous section. The key differences are:

  1. The ORDERING key has changed from NESTED to NUNIQ.
  2. The INDXSCHM key has changed from IMPLICIT to EXPLICIT.
  3. There is an extra column, UNIQ, that explicitly identifies each pixel in the UNIQ indexing scheme.
  4. The PROB column has been renamed to PROBDENSITY, and the units have change from probability to probability per steradian.
Reading Multi-Resolution Sky Maps

Now let’s go through some of the same common HEALPix operations from the previous section, but using the multi-resolution format. Instead of Healpy, we will use astropy-healpix because it has basic support for the UNIQ indexing scheme.

First, we need the following imports:

>>> from astropy.table import Table
>>> from astropy import units as u
>>> import astropy_healpix as ah
>>> import numpy as np

Next, let’s read the sky map. Instead of a special-purpose HEALPix method, we just read the FITS file into an Astropy table using Astropy’s unified file read/write interface:

>>> skymap = Table.read('bayestar.multiorder.fits')
Most Probable Sky Location

Next, let’s find the highest probability density sky position. This is a three-step process.

  1. Find the UNIQ pixel index of the highest probability density tile:

    >>> i = np.argmax(skymap['PROBDENSITY'])
>>> uniq = skymap[i]['UNIQ']

    What is the probability density per square degree in that tile?

    >>> skymap[i]['PROBDENSITY'] * (np.pi / 180)**2
0.0782516470191411

  2. Unpack the UNIQ pixel index into the resolution, nside, and the NESTED pixel index, ipix, using the method astropy_healpix.uniq_to_level_ipix(). (Note that this method returns level, which is the logarithm base 2 of nside, so we must also convert from level to nside using astropy_healpix.level_to_nside().)

    >>> level, ipix = ah.uniq_to_level_ipix(uniq)
>>> nside = ah.level_to_nside(level)

  3. Convert from nside and ipix to right ascension and declination using astropy_healpix.healpix_to_lonlat() (which is equivalent to hp.pix2ang):

    >>> ra, dec = ah.healpix_to_lonlat(ipix, nside, order='nested')
>>> ra.deg
194.30419921874997
>>> dec.deg
-17.856895095545468
Probability Density at a Known Position

Now let’s look up the probability density at a known sky position. In this case, let’s use the position of NGC 4993:

>>> ra = 197.4133 * u.deg
>>> dec = -23.3996 * u.deg

Brute Force Linear Search

The following brute force method of looking up a pixel by sky position has a complexity of \(O(N)\), where \(N\) is the number of multi-resolution pixels.

  1. Unpack the UNIQ pixel indices into their resolution and their NESTED pixel index.

    >>> level, ipix = ah.uniq_to_level_ipix(skymap['UNIQ'])
>>> nside = ah.level_to_nside(level)

  2. Determine the NESTED pixel index of the target sky position at the resolution of each multi-resolution tile.

    >>> match_ipix = ah.lonlat_to_healpix(ra, dec, nside, order='nested')

  3. Find the multi-resolution tile whose NESTED pixel index equals the target pixel index.

    >>> i = np.flatnonzero(ipix == match_ipix)[0]
>>> i
13484

    That pixel contains the target sky position.

    >>> skymap[i]['PROBDENSITY'] * (np.pi / 180)**2
0.03467919098907807

Fast Binary Search

The following binary search method of looking up a pixel by sky position exploits the algebraic properties of HEALPix. It has a complexity of \(O(\log N)\) where \(N\) is the number of multi-resolution pixels. It assumes that every sky position is mapped on to exactly one multi-resolution tile, which is true for LIGO/Virgo multi-resolution sky maps.

  1. First, find the NESTED pixel index of every multi-resolution tile, at an arbitrarily high resolution. (nside = 2**29 works nicely because it is the highest possible HEALPix resolution that can be represented in a 64-bit signed integer.)

    >>> max_level = 29
>>> max_nside = ah.level_to_nside(max_level)
>>> level, ipix = ah.uniq_to_level_ipix(skymap['UNIQ'])
>>> index = ipix * (2**(max_level - level))**2

    Sort the pixels by this value.

    >>> sorter = np.argsort(index)

  2. Determine the NESTED pixel index of the target sky location at that resolution.

    >>> match_ipix = ah.lonlat_to_healpix(ra, dec, max_nside, order='nested')

    Do a binary search for that value.

    >>> i = sorter[np.searchsorted(index, match_ipix, side='right', sorter=sorter) - 1]
>>> i
13484

    That pixel contains the target sky position.

    >>> skymap[i]['PROBDENSITY'] * (np.pi / 180)**2
0.03467919098907807
[1](1, 2) Singer, L. P., & Price, L. R. 2016, Phys. Rev. D, 93, 024013. doi:10.1103/PhysRevD.93.024013
[2]Górski, K. M., Wandelt, B. D., et al. 1999. arXiv:astro-ph/9905275
[3]Fernique, P., Allen, et al. 2015, Astron. Astrophys., 578, A114. doi:10.1051/0004-6361/201526075
[4]Fernique, P., Boch, T., et al. 2014, IVOA Recommendation. arXiv:1505.02937

Basic Observability Calculations

Now we are going to teach our GCN handler how to determine whether a gravitational-wave event is observable. We are going to use the astropy.coordinates module. (See also the Astropy example on observation planning in Python.) First, we will need to import a few extra Python modules:

import astropy.coordinates
import astropy.time
import astropy.units as u

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 Astropy to transform those coordinates to a horizontal (altitude–azimuth) 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

You 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. For more advanced functionality, we recommend the astroplan package.

def prob_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 map
    npix = len(m)
    nside = hp.npix2nside(npix)

    # Get time now
    time = 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, now
    frame = 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, now
    altaz = 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!
    return prob

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

@gcn.handlers.include_notice_types(
    gcn.notice_types.LVC_PRELIMINARY,
    gcn.notice_types.LVC_INITIAL,
    gcn.notice_types.LVC_UPDATE)
def process_gcn(payload, root):
    # 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':
    #    return
    if root.attrib['role'] != 'test':
        return

    # Respond only to 'CBC' events. Change 'CBC' to "Burst'
    # to respond to only unmodeled burst events.
    if root.find(".//Param[@name='Group']").attrib['value'] != 'CBC':
        return

    skymap_url = root.find(".//Param[@name='skymap_fits']").attrib['value']

    skymap, header = hp.read_map(skymap_url, h=True, verbose=False)
    prob = prob_observable(skymap, header)
    print('Source has a {:d}% chance of being observable now'.format(
        int(round(100 * prob))))
    if prob > 0.5:
        pass # FIXME: perform some action

Let’s run the new GCN handler now…

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

When you run this script, each time you receive a sample LIGO/Virgo GCN Notice, it will print something like the following (note that probability will change as a function of time):

Source has a 76% chance of being observable now

Distance and 3D Sky Maps

All sky localization FITS files for CBC events are three dimensional: they include both the sky probability map and a directionally dependent distance estimate. This can be useful for identifying possible host galaxies using a galaxy redshift catalog. For further details see [1] and for sample code see [2].

https://dcc.ligo.org/public/0119/P1500071/007/demonstration.png https://dcc.ligo.org/public/0119/P1500071/007/proj3d.png
[1]Singer, L. P., Chen, H.-Y., Holz, D. E., et al. 2016, Astropys. J. Lett., 829, L15. doi:10.3847/2041-8205/829/1/L15
[2]Singer, L. P., Chen, H.-Y., Holz, D. E., et al. 2016, Astropys. J. Supp., 226, 10. doi:10.3847/0067-0049/226/1/10

Additional Tools

The ligo.skymap package includes a number of advanced tools for working with GW probability sky maps.

Appendix

Change Log

Version 10 (2019-07-31)

Data Analysis

  • Add a more detailed description of the RAVEN pipeline. Previously just mentioned types of searches but now has information on external experiments, coincident searches, and coincident false alarm rates.
  • Fixed PyCBC Live reference.

Alert Contents

  • Changed the data type of the UNIQ column of the multi-order sky map format from an unsigned integer to a signed integer as specified by the MOC-in-FITS standard.

    This will improve interoperability with the mrdfits function from the IDL Astronomy User’s Library and the fv FITS Viewer from FTOOLS, both of which were reported to have problems with the old unsigned integer column. It will also make it simpler to work with Numpy indexing operations, since Numpy uses a signed integer type for indexing.

    This change will go into effect in the LIGO/Virgo low-latency alert system on 2019-08-06.

    Users of ligo.skymap should update to version 0.1.8 or newer because older versions will be unable to read old files with unsigned UNIQ columns. The new version of ligo.skymap can read files with either signed or unsigned UNIQ columns.

Version 9 (2019-06-13)

General

Data Analysis

  • Renamed this section from “Procedures” to “Data Analysis” and reordered its subsections to better reflect the chronological order of the steps of the analysis.

Sample Code

  • Add tutorial on working with multi-resolution sky maps.
  • Add sample code to test whether a sky position is in the 90% credible region.
  • Add sample code to find the area of the 90% credible region.

Version 8 (2019-05-22)

Alert Contents

  • Describe the two localization formats, the official *.fits.gz HEALPix-in-FITS format and the experimental multi-resolution HEALPix *.multiorder.fits format.

    Effective 2019-05-28, the multi-resolution file suffix will be renamed from *.fits to *.multiorder.fits. The old *.fits suffix had caused confusion because the multi-resolution format is not the same as the *.fits.gz files without gzip compression.

    The multi-resolution format is currently recommended only for advanced users. Tutorials and sample code will soon be included in an upcoming version of the User Guide.

  • Update the description of the HasNS property in the sample GCN Notices. Previously, it was defined as the probability that at least one object in the binary has a mass that is less than 2.83 solar masses. Now, for consistency with the source classification definitions, it is defined as the probability that at least one object in the binary has a mass that is less than 3 solar masses.

  • Add some shading to the source classification diagram to make it clear that the definitions of the source classes are symmetric under exchange of the component masses, but that by convention the component masses are defined such that \(m_1 \geq m_2\).

Version 7.1 (2019-03-02)

  • Remove the warning on the front page about significant changes to this document occurring before the start of O3.

Version 7 (2019-03-02)

Observing Capabilities

  • Record the official start of O3.

Procedures

  • Add Gravitational Wave/High Energy Neutrino search to the list of multi-messenger search pipelines.

Sample Code

  • Add instructions for installing required packages using the Anaconda Python distribution.

Version 6 (2019-03-08)

Alert Contents

  • Switch to the GW170817 Hanford-Livingston-Virgo localization for the example sky map.

Version 5 (2019-03-01)

Alert Contents

  • Add a human-readable description to the Pkt_Ser_Num parameter.
  • Add <EventIVORN cite="supersedes"> elements to the sample Initial and Update notices in order to cite all prior VOEvents. GraceDB already includes this metadata, but it was missing from the examples.
  • Add MassGap classification for compact binary mergers.

Version 4 (2019-02-15)

General

Getting Started Checklist

  • Update links for OpenLVEM enrollment instructions.

Observing Capabilities

  • Update planned dates for Engineering Run 14 (ER14; starts 2019-03-04) and Observing Run 3 (O3; starts 2019-04-01).
  • Add Live Status section, listing some public web pages that provide real-time detector status.

Sample Code

  • Update the example GCN notice handler so that it does not fail if the notice is missing a sky map, because as we have specified them, LVC_RETRACTION notices never contain sky maps and LVC_PRELIMINARY notices may or may not contain sky maps.
  • When building the documentation, test all of the sample code automatically.

Version 3 (2019-02-13)

Alert Contents

  • Remove the skymap_png parameter from the VOEVents. The sky map plots take longer to generate than the FITS files themselves, so they would have needlessly delayed the preliminary alerts.
  • Change the IVORN prefix from ivo://gwnet/gcn_sender to ivo://gwnet/LVC, because GCN traditionally uses the text after the / to indicate the name of the mission, which most closely corresponds to “LVC,” short for “LIGO/Virgo Collaboration.” Note that this IVORN is used for historical purposes and continuity with the GCN notice format used in O1 and O2, and is likely to change in the future with the commissioning of additional gravitational-wave facilities.
  • Retraction notices now get a separate GCN notice packet type, LVC_RETRACTION=164. The Retraction parameter has been removed from the <What> section.

Version 2 (2018-12-13)

Alert Contents

  • Removed the Vetted parameter from GCN Notices. It was intended to indicate whether the event had passed human vetting. However, it was redundant because by definition Preliminary events have not been vetted and all Initial and Update alerts have been vetted.
  • The type of the Retraction parameter in the GCN Notices was changed from a string (false or true) to an integer (0 or 1) for consistency with other flag-like parameters.
  • Remove the units attribute from parameters that are not numbers.

Sample Code

  • GCN has now begun publicly broadcasting sample LIGO/Virgo GCN Notices. Updated the sample code accordingly with instructions for receiving live sample notices.

Version 1 (2018-11-27)

Getting Started Checklist

  • Updated instructions for joining the OpenLVEM Community.

Observing Capabilities

  • Changed the expected number of BNS events in O3 from 1-50, as stated in the latest version of the Living Review, to 1-10 events, as stated in the more recent rates presentation.

Alert Contents

  • In the example VOEvents, moved the Classification and Inference quantities from the <Why> section to the <What> section so that they validate against the VOEvent 2.0 schema.

Glossary

BBH
Binary black hole system, a binary system composed of two black holes.
BNS
Binary neutron star, a binary system composed of two neutron stars.
BNS range
A figure of merit to describe the sensitivity of a gravitational-wave detector to BNS mergers, defined as the average luminosity distance at which the merger of two \(1.4~M_\odot\) objects would be detectable with a signal to noise ratio of 8. See also Burst range.
Burst range
A figure of merit to describe the sensitivity of a gravitational-wave detector to unmodeled bursts, defined with reference to optimistic models of gravitational wave emissions from stellar collapse as the average luminosity distance at which a monochromatic burst at 150 Hz with a fluence of \(E_\mathrm{GW} = 10^{-2} M_\odot c^2\) would be detectable with a signal to noise ratio of 8. See also BNS range.
CBC
Compact binary coalescence.
FITS
Flexible Image Transport System, a format for astronomical tables, images, and multidimensional data sets.
GCN
The Gamma-ray Coordinates Network (https://gcn.gsfc.nasa.gov), a portal for discoveries and observations of astronomical transients. Historically, GCN has served high-energy satellites but now also other electromagnetic wavelengths and also gravitational-wave, cosmic ray, and neutrino facilities.
GCN Circular
A human-readable astronomical bulletin distributed through GCN.
GCN Notice
A machine-readable alert distributed through GCN.
GraceDB
Gravitational Wave Candidate Event Database (https://gracedb.ligo.org), the official public marshal portal for LIGO/Virgo candidates.
GRB
Gamma-ray burst.
HEALPix
Hierarchical Equal Area isoLatitude Pixelation, a scheme for indexing positions on the unit sphere.
HEN
High Energy Neutrino, particularly in the context of multi-messenger GW+HEN follow-up.
KAGRA
Kamioka Gravitational Wave Detector, an underground gravitational-wave detector in the Kamioka mine in Japan.
LHO
LIGO Hanford Observatory (see LHO observatory home page), site of a 4 km gravitational-wave detector in Hanford, Washington, USA.
LLO
LIGO Livingston Observatory (see LLO observatory home page), site of a 4 km gravitational-wave detector in Livingston, Louisiana, USA.
MassGap
Compact binary systems with at least one compact object whose mass is in the hypothetical “mass gap” between neutron stars and black holes, defined here as 3-5 solar masses.
NSBH
Neutron star black hole, a binary system composed of one neutron star and one black hole.
primary
When referring to the two component compact objects or the masses of the two component compact objects in a binary, the primary is the more massive one, i.e., \(m_1 \geq m_2\). See secondary.
secondary
When referring to the two component compact objects or the masses of the two component compact objects in a binary, the secondary is the less massive one, i.e., \(m_2 \leq m_1\). See primary.
SN
Supernova.
SNR
Signal-to-noise ratio, here applied to gravitational-wave signals. It is defined the square root of the integral over frequency of the power spectral density of the gravitational-wave signal over the integral over frequency of the average power spectral density of the noise.
Terrestrial
Classification for signals in gravitational-wave detectors that are of instrumental or environmental origin. Terrestrial signals are not astrophysical and not due to gravitational waves. Some examples of sources of terrestrial signals are statistical noise fluctuations, detector glitches, and ground motion.
Virgo
Virgo Observatory (see Virgo observatory home page), site of a 3 km gravitational-wave detector in Cascina, Italy.
VOEvent
An XML format for describing astronomical transients. For the specification, see the official VOEvent IVOA Recommendation.
VTP
VOEvent Transport Protocol, a simple TCP-based protocol for sending and receiving VOEvents, used by GCN. For the specification, see the official VTP IVOA recommendation.