#!/usr/bin/env python # Copyright (C) 2015 Ian Harry # # 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, write to the Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. """ Take a coinc xml file containing multiple events and upload to gracedb. """ import argparse import logging from ligo.gracedb.rest import GraceDb import numpy as np import h5py import pycbc import lal import lal.series from glue.ligolw import ligolw from glue.ligolw import table from glue.ligolw import lsctables from glue.ligolw import utils as ligolw_utils from glue.ligolw.utils import process as ligolw_process from ligo.segments import segment, segmentlist class LIGOLWContentHandler(ligolw.LIGOLWContentHandler): pass logging.basicConfig(format='%(asctime)s:%(levelname)s : %(message)s', level=logging.INFO) parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--psd-files", nargs='+', required=True, help='HDF file(s) containing the PSDs to upload') parser.add_argument('--input-file', dest='input_file', required=True, type=str, help='Input LIGOLW XML file of coincidences.') parser.add_argument('--testing', action="store_true", default=False, help="Upload event to the TEST group of gracedb.") parser.add_argument('--min-ifar', type=float, metavar='YEARS', help='Only upload events more significant than given IFAR') parser.add_argument('--production-server', action="store_true", default=False, help="Upload event to production graceDB. If not given " "events will be uploaded to playground server.") args = parser.parse_args() if args.production_server: gracedb = GraceDb() else: gracedb = GraceDb(service_url='https://gracedb-playground.ligo.org/api/') lsctables.use_in(LIGOLWContentHandler) xmldoc = ligolw_utils.load_filename(args.input_file, contenthandler=LIGOLWContentHandler) class psd_segment(segment): def __new__(cls, psd, *args): return segment.__new__(cls, *args) def __init__(self, psd, *args): self.psd = psd psds = {} for psd_file in args.psd_files: (ifo, group), = h5py.File(psd_file, "r").items() psd = [group["psds"][str(i)] for i in range(len(group["psds"].keys()))] psds[ifo] = segmentlist(psd_segment(*segargs) for segargs in zip( psd, group["start_time"], group["end_time"])) coinc_table = table.get_table(xmldoc, lsctables.CoincTable.tableName) coinc_inspiral_table = table.get_table(xmldoc, lsctables.CoincInspiralTable.tableName) coinc_event_map_table = table.get_table(xmldoc, lsctables.CoincMapTable.tableName) sngl_inspiral_table = table.get_table(xmldoc, lsctables.SnglInspiralTable.tableName) xmldoc.childNodes[-1].removeChild(sngl_inspiral_table) xmldoc.childNodes[-1].removeChild(coinc_event_map_table) xmldoc.childNodes[-1].removeChild(coinc_inspiral_table) xmldoc.childNodes[-1].removeChild(coinc_table) for event in coinc_table: coinc_event_table_curr = lsctables.New(lsctables.CoincTable) coinc_event_table_curr.append(event) coinc_inspiral_table_curr = lsctables.New(lsctables.CoincInspiralTable) coinc_event_map_table_curr = lsctables.New(lsctables.CoincMapTable) sngl_inspiral_table_curr = lsctables.New(lsctables.SnglInspiralTable) coinc_event_id = event.coinc_event_id for coinc_insp in coinc_inspiral_table: if coinc_insp.coinc_event_id == event.coinc_event_id: coinc_inspiral_table_curr.append(coinc_insp) if args.min_ifar is not None and \ coinc_inspiral_table_curr[0].combined_far > 1./args.min_ifar/lal.YRJUL_SI: continue sngl_ids = [] for coinc_map in coinc_event_map_table: if coinc_map.coinc_event_id == event.coinc_event_id: coinc_event_map_table_curr.append(coinc_map) sngl_ids.append(coinc_map.event_id) psddict = {} for sngl in sngl_inspiral_table: if sngl.event_id in sngl_ids: sngl_inspiral_table_curr.append(sngl) try: psd = psds[sngl.ifo] except KeyError: parser.error( "--psd-files {0}: no PSDs found for detector {1}".format( " ".join(args.psd_files), sngl.ifo)) try: psd = psd[psd.find(sngl.get_end())].psd except ValueError: parser.error( "--psd-files {0}: no PSD found for detector {1} " "at GPS time {2}".format( " ".join(args.psd_files), sngl.ifo, sngl.get_end())) flow = psd.file.attrs['low_frequency_cutoff'] df = psd.attrs["delta_f"] kmin = int(flow / df) fseries = lal.CreateREAL8FrequencySeries( "psd", psd.attrs["epoch"], kmin * df, df, lal.StrainUnit**2 / lal.HertzUnit, len(psd.value) - kmin) fseries.data.data = psd.value[kmin:] / np.square(pycbc.DYN_RANGE_FAC) psddict[sngl.ifo] = fseries xmldoc.childNodes[-1].appendChild(coinc_event_table_curr) xmldoc.childNodes[-1].appendChild(coinc_inspiral_table_curr) xmldoc.childNodes[-1].appendChild(coinc_event_map_table_curr) xmldoc.childNodes[-1].appendChild(sngl_inspiral_table_curr) ligolw_utils.write_filename(xmldoc, "tmp_coinc_xml_file.xml") psd_xmldoc = lal.series.make_psd_xmldoc(psddict) ligolw_utils.write_filename(psd_xmldoc, "tmp_psd.xml.gz", gz=True) if args.testing: r = gracedb.createEvent("Test", "pycbc", "tmp_coinc_xml_file.xml", search="AllSky", offline=True).json() else: r = gracedb.createEvent("CBC", "pycbc", "tmp_coinc_xml_file.xml", search="AllSky", offline=True).json() logging.info("Uploaded event %s.", r["graceid"]) gracedb.writeLog( r["graceid"], "PyCBC PSD estimate from the time of event", "psd.xml.gz", open("tmp_psd.xml.gz", "rb").read(), "psd").json() logging.info("Uploaded file psd.xml.gz to event %s.", r["graceid"]) xmldoc.childNodes[-1].removeChild(coinc_event_table_curr) xmldoc.childNodes[-1].removeChild(coinc_inspiral_table_curr) xmldoc.childNodes[-1].removeChild(coinc_event_map_table_curr) xmldoc.childNodes[-1].removeChild(sngl_inspiral_table_curr)