#!/usr/bin/python
#
# Project Librarian: Brandon Piotrzkowski
# Graduate Student
# UW-Milwaukee Department of Physics
# Center for Gravitation & Cosmology
# <brandon.piotrzkowski@ligo.org>
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
"""
Script to calculate the joint false alarm rate of the provided candidates,
either the spatiotemporal variant if sky information is provided, or
just temporal if not.
"""
__author__ = "Brandon Piotrzkowski <brandon.piotrzkowski@ligo.org>"
# Global imports.
import numpy as np
from ligo.raven import search
from ligo.gracedb.rest import GraceDb
import argparse
# Command line options.
parser = argparse.ArgumentParser(description='Perform query of GraceDB')
parser.add_argument("-s", "--superevent", metavar="(S,MS)XXXXXXX",
default=None,
help="GraceDB ID of superevent (required)"),
parser.add_argument("-e", "--ext_event", metavar="(E,M)XXXXXXX",
default=None,
help="GraceDB ID of external event (required)"),
parser.add_argument("-w", "--window", nargs=2, metavar="t", type=int,
action="append", default=None,
help="Time window [tl, th] seconds to search around event time (required with th > th)"),
parser.add_argument("-g", "--ext_search", metavar="GRB SubGRB SubGRBTargeted HEN MDC",
default=None,
help="Search of external event."),
parser.add_argument("-S", "--superevent_skymap", metavar="bayestar.fits.gz",
default=None,
help="fits(.gz) filename for superevent sky map"),
parser.add_argument("-E", "--ext_event_skymap", metavar="'glg_healpix_all_bn_v00.fit'",
default=None,
help="fits(.gz) filename for external event sky map"),
parser.add_argument("-u", "--gracedb_url", metavar="https://gracedb.ligo.org/api/",
default='',
help="Gracedb url"),
parser.add_argument("-f", "--far_ext", metavar="1e-6", default=None,
help="False alarm rate of GRB"),
parser.add_argument("-n", "--em_rate", metavar="1e-3", default=None,
help="False alarm rate of GRB"),
parser.add_argument("-t", "--far_gw_thresh", metavar="1.16e-5", default=None,
help="Maximum cutoff for GW FAR considered in the search"),
parser.add_argument("-T", "--far_ext_thresh", metavar="1e-4", default=None,
help="Maximum cutoff for GRB FAR considered in the search"),
parser.add_argument("-m", "--se_moc", dest="se_moc", action="store_true",
help="Assume the superevent sky map is multi-ordered (MOC)."),
parser.add_argument("-M", "--ext_moc", dest="ext_moc", action="store_true",
help="Assume the external sky map is multi-ordered (MOC)."),
parser.add_argument("-r", "--se_ring", dest="se_ring", action="store_true",
help="Assume the superevent map uses RING ordering rather than nested."),
parser.add_argument("-R", "--ext_ring", dest="ext_ring", action="store_true",
help="Assume the external sky map uses RING ordering rather than nested."),
parser.add_argument("-c", "--use_radec", dest="use_radec", action="store_true",
help="Choose to use RA and dec of external sky map if a single pixel sky map."),
parser.add_argument("-p", "--use_preferred_event_skymap", dest="use_preferred_event_skymap", action="store_true",
help="Choose to use the GW sky map in the preferred event rather than the superevent.")
args = parser.parse_args()
print(args)
# Check required options are there
if not args.superevent:
raise AssertionError('superevent graceid not given')
else:
se_id = str(args.superevent)
if not args.ext_event:
raise AssertionError('external event graceid not given')
else:
exttrig_id = str(args.ext_event)
if not args.window:
raise AssertionError('search window not given')
else:
tl, th = int(args.window[0][0]), int(args.window[0][1])
# Load other options
if args.gracedb_url:
gracedb = GraceDb(str(args.gracedb_url))
else:
gracedb = None
ext_search = args.ext_search
se_fitsfile = args.superevent_skymap
ext_fitsfile = args.ext_event_skymap
if args.far_ext:
far_ext = float(args.far_ext)
else:
far_ext = None
if args.em_rate:
em_rate = float(args.em_rate)
else:
em_rate = None
if args.far_gw_thresh:
far_gw_thresh = float(args.far_gw_thresh)
else:
far_gw_thresh = None
if args.far_ext_thresh:
far_ext_thresh = float(args.far_ext_thresh)
else:
far_ext_thresh = None
if se_fitsfile and (ext_fitsfile or args.use_radec):
incl_sky = True
else:
incl_sky = False
results = search.coinc_far(
se_id, exttrig_id, tl, th, ext_search=ext_search,
se_fitsfile=se_fitsfile, ext_fitsfile=ext_fitsfile,
incl_sky=incl_sky, gracedb=gracedb,
far_ext=far_ext, em_rate=em_rate,
far_gw_thresh=far_gw_thresh,
far_ext_thresh=far_ext_thresh,
se_moc=args.se_moc, ext_moc=args.ext_moc,
se_nested=not args.se_ring, ext_nested=not args.ext_ring,
use_radec=args.use_radec,
use_preferred_event_skymap=args.use_preferred_event_skymap)
print(results)
