Commit 4a87cdc2 authored by Naresh Adhikari's avatar Naresh Adhikari
Browse files

Add joint-far section and update docstrings

parent 7e364e4f
Pipeline #413796 passed with stages
in 6 minutes and 30 seconds
......@@ -49,6 +49,7 @@ extensions = [
'sphinx.ext.extlinks',
'sphinx.ext.graphviz',
'sphinx.ext.intersphinx',
'sphinx.ext.mathjax',
'sphinx.ext.napoleon',
'sphinx.ext.viewcode']
......@@ -252,4 +253,3 @@ autosectionlabel_prefix_document = True
# 'arxiv': ('https://arxiv.org/abs/%s', 'arXiv:'),
# 'doi': ('https://doi.org/%s', 'doi:')
#}
Joint FAR calculation
---------------------
=====================
One of the main methods to distinguish between a real or noise triggers of any
candidate is to calculate the False Alarm Rate (FAR). The main motivation behind
using the joint FAR between gravitational wave (GW)
and electromagnetic (EM) triggers are to indicate which candidate should receive
additional follow-up.
Given the information about the gravitational wave candidate and its
associated astrophysical counterpart, the coincident false alarm rate
:math:`FAR_{c}` is calculated depending on the search method.
Two search methods are:
Untargeted search method
------------------------
In the untargeted search method, we look for independent gravitational wave and electromagnetic wave (EM) candidates.
GRBs of higher statistical significance are considered and GWs can be of very low significance.
So, for the untargeted search analysis, the joint coincident FAR is
.. math::
FAR_{c} = FAR_{gw} \frac{R_{grb} \Delta t}{I_{\Omega}}
where :math:`FAR_{gw}` is the false alarm rate of the GW candidate given by one of the
GW pipelines such as GstLAL, MBTAOnline, PyCBC Live, SPIIR, or cWB.
:math:`R_{grb}` are the
rate of the unique GRB candidates given by all participating EM signal
detectors, :math:`\Delta t` is the coincident time window of the GW-GRB signal, and
:math:`I_{\Omega}` is the sky map overlap integral outlined below.
Targeted search method
----------------------
In the targeted search method, we respond to GRB triggers coming from GRB experiments and look for gravitational wave events
happened around that time. Because of this,
we consider sub-threshold GRB candidates, which now should dominate the rate compared to real detections.
Since the GRB candidates were searched for only around GW times making this a targeted search and
the joint coincident FAR for this search method is calculated as
.. math::
FAR_c = \frac{Z}{I_{\Omega}} \left( 1 + \ln(\frac{Z_{max}}{Z}) \right)
where Z is the joint ranking statistic defined as:
.. math::
Z = FAR_{gw} FAR_{grb} \Delta t
and the max Z value is calculated using the maximum GW and GRB FAR thresholds.
Skymap overlap integral
-----------------------
In RAVEN, we have built functionality to use different skymap formats for the
skymap overlap integral as first derived in `Ashton et. al. (2018) <https://iopscience.iop.org/article/10.3847/1538-4357/aabfd2/pdf>`_
:
.. math::
I_{\Omega}(x_{gw}, x_{grb}) = \int \frac{p({\Omega}| x_{gw}, \mathcal{H}_{gw}^s)
p({\Omega}| x_{grb}, \mathcal{H}_{grb}^s)}{p({\Omega}| \mathcal{H}^s)} d{\Omega}
Two methods to incorporate the skymap in calculating the skymap overlap integral
(:math:`I_{\Omega}`) are listed as follows:
Hierarchical Equal Area isoLatitude Pixelization (HEALPix) method
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
In this method, the sky map overlap integral is calculated as
.. math::
I_{\Omega}(x_{gw}, x_{grb}) = N_{pix} \sum_{i=0}^{N_{pix}} P(i| x_{gw}, \mathcal{H}_{gw}^s)
\cdot P(i| x_{grb}, \mathcal{H}_{grb}^s)
where :math:`N_{pix}` is the number of pixels in a skymap.
Multi-Ordered Coverage (MOC) method
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
In this method, we use RA and DEC of given pixel :math:`\Omega_{i}` to match
pixels of MOC skymaps and calculate the overlap integral as
.. math::
I_{\Omega}(x_{1}, x_{2}) = 4 \pi \sum_{i_{1}}^{N_{pix}} p(\Omega_{1}(i_{1})| x_{1}, \mathcal{H}^s)
\cdot p(\Omega_{1}(i_{1})| x_{2}, \mathcal{H}^s) \cdot \Delta A (i_{1})
where :math:`p(\Omega_{1}(i)| x_{a}, \mathcal{H}^s)` is the normalized probability density
such that :math:`\sum_{i}p(\Omega_{1}(i)| x_{a}, \mathcal{H}^s)\Delta A = 1.`
The UNIQ ordering system has an advantage over the standard HEALPix system where equal-area pixels are used.
MOC skymaps can increase resolution around the higher probability regions and vice versa, reducing file size.
Moreover, we find using MOC skymaps gives the same result compared to the HEALPix method with much lower latency.
......@@ -47,7 +47,29 @@ from gracedb_sdk import Client
def query(event_type, gpstime, tl, th, gracedb=None, group=None,
pipelines=[], searches=[], se_searches=[]):
""" Query for coincident events of type event_type occurring within a
window of [tl, th] seconds around gpstime. """
window of [tl, th] seconds around gpstime.
Parameters
----------
event_type: str
"Superevent" or "External"
gpstime: float
Event's gps time
tl: float
Start of coincident time window
th: float
End of coincident time window
gracedb: class
SDK REST API client for HTTP connection
group: str
"CBC", "Burst", or "Test",
pipelines: list
List of external trigger pipeline names
searches: array
List of superevent searches
se_searches: array
List of external trigger searches
"""
# Perform a sanity check on the time window.
if tl >= th:
......@@ -102,8 +124,29 @@ def search(gracedb_id, tl, th, gracedb=None, group=None, pipelines=None,
searches=[], se_searches=[], event_dict=None):
""" Perform a search for neighbors coincident in time within
a window [tl, th] seconds around an event. Uploads the
results to the selected gracedb server. """
results to the selected gracedb server.
Parameters
----------
gracedb_id: str
ID of the trigger used by GraceDB
tl: float
Start of coincident time window
th: float
End of coincident time window
gracedb: class
SDK REST API client for HTTP connection
group: string
"CBC", "Burst", or "Test",
pipelines: list
List of external trigger pipeline names
searches: array
List of superevent searches
se_searches: array
List of external trigger searches
event_dict: dict
Dictionary of the gracedb event
"""
# Identify neighbor types with their graceid strings.
types = {'G': 'GW', 'E': 'External', 'S': 'Superevent',
'T': 'Test'}
......@@ -226,7 +269,7 @@ def skymap_overlap_integral(se_skymap, exttrig_skymap=[],
or a flattened sky map with Nested or Ring ordering can be used.
Either a mutli-ordered (MOC) external sky map with UNIQ ordering,
flattened sky map with Nested or Ring ordering,
or a position indicated by RA/dec can be used.
or a position indicated by RA/DEC can be used.
Parameters
----------
......@@ -393,26 +436,26 @@ def odds_ratio(skymap_overlap, tl, th, bayes_gw=0., bayes_grb=0.,
Sky map overlap integral (also could be a generic
overlap integral)
tl: float
start of coincident time window
Start of coincident time window
th: float
end of coincident time window
bayes_gw: array-like
End of coincident time window
bayes_gw: array
Noise vs signal Bayes factor for GW
bayes_grb: array-like
bayes_grb: array
Noise vs signal Bayes factor for GRB
r_c: array-like
r_c: array
Rate of coincidences (per unit time, using
astropy units)
r_gw: array-like
r_gw: array
Rate of real GW candidates considered for search
(per unit time, using astropy units)
r_grb: array-like
r_grb: array
Rate of real GRB candidates considered for search
(per unit time, using astropy units)
r_n_gw: array-like
r_n_gw: array
Rate of noise GW candidates in search
(per unit time, using astropy units)
r_n_grb: array-like
r_n_grb: array
Rate of noise GRB candidates in search
(per unit time, using astropy units)
......@@ -436,7 +479,50 @@ def coinc_far(se_id, ext_id, tl, th, grb_search='GRB', se_fitsfile=None,
""" Calculate the significance of a gravitational wave candidate with the
addition of an external astrophyical counterpart in terms of a
coincidence false alarm rate. This includes a temporal and a
space-time type. """
space-time type.
Parameters
----------
se_id: str
GraceDB ID of superevent
ext_id: str
GraceDB ID of external event
tl: float
Start of coincident time window
th: float
End of coincident time window
grb_search: str
Determines joint FAR method.
"GRB", "SubGRB", or "SubGRBTargeted"
se_fitsfile: str
Superevent's skymap file name
ext_fitsfile: str
External event's skymap file name
incl_sky: bool
If True, uses skymaps in the joint FAR calculation
gracedb: class
SDK REST API client for HTTP connection
far_grb: float
GRB false alarm rate
em_rate: float
Detection rate of external events
se_dict: float
Dictionary of superevent
ext_dict: float
Dictionary of external event
se_moc: bool
If True, assumes multi-order coverage (MOC) superevent skymap
ext_moc: bool
If True, assumes multi-order coverage (MOC) external event skymap
use_radec: bool
If True, use ra and dec for single pixel external skymap
se_nested: bool
If True, assumes GW skymap uses nested ordering, otherwise
assumes ring ordering
ext_nested: bool
If True, assumes external skymap uses nested ordering, otherwise
assumes ring ordering
"""
# Create the SE and ExtTrig objects based on string inputs.
se = SE(se_id, fitsfile=se_fitsfile, gracedb=gracedb, event_dict=se_dict,
......@@ -531,7 +617,50 @@ def calc_signif_gracedb(se_id, ext_id, tl, th, grb_search='GRB',
se_moc=False, ext_moc=False,
use_radec=False, se_nested=False, ext_nested=False):
""" Calculates and uploads the coincidence false alarm rate
of the given superevent to the selected gracedb server. """
of the given superevent to the selected gracedb server.
Parameters
----------
se_id: str
GraceDB ID of superevent
ext_id: str
GraceDB ID of external event
tl: float
Start of coincident time window
th: float
End of coincident time window
grb_search: str
Determines joint FAR method.
"GRB", "SubGRB", or "SubGRBTargeted"
se_fitsfile: str
Superevent's skymap file name
ext_fitsfile: str
External event's skymap file name
incl_sky: bool
If True, uses skymaps in the joint FAR calculation
gracedb: class
SDK REST API client for HTTP connection
far_grb: float
GRB false alarm rate
em_rate: float
Detection rate of external events
se_dict: float
Dictionary of superevent
ext_dict: float
Dictionary of external event
se_moc: bool
If True, assumes multi-order coverage (MOC) superevent skymap
ext_moc: bool
If True, assumes multi-order coverage (MOC) external event skymap
use_radec: bool
If True, use ra and dec for single pixel external skymap
se_nested: bool
If True, assumes GW skymap uses nested ordering, otherwise
assumes ring ordering
ext_nested: bool
If True, assumes external skymap uses nested ordering, otherwise
assumes ring ordering
"""
# Create the SE and ExtTrig objects based on string inputs.
se = SE(se_id, fitsfile=se_fitsfile, gracedb=gracedb, event_dict=se_dict,
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment