gstlal_inspiral_check_livetimes

#!/usr/bin/env python
#
# Copyright (C) 2015 Chad Hanna
#
# 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 2 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.

import sys
import sqlite3
from optparse import OptionParser
from ligo.lw import dbtables
from ligo.lw import ligolw
from ligo.lw import lsctables
from ligo.lw import utils as ligolw_utils
from ligo.lw.utils import segments as ligolw_segments
from ligo import segments
from ligo.segments import utils as segmentsUtils

from gstlal import far

def parse_command_line():
	parser = OptionParser(usage="%prog [options] database.sqlite")
	parser.add_option("--post-marginalized-likelihood-file", metavar = "filename", help = "Path of post marginalized likelihood file")
	options, database = parser.parse_args()

	return options, database


options, database = parse_command_line()

_, rankingstatpdf = far.parse_likelihood_control_doc(ligolw_utils.load_filename(options.post_marginalized_likelihood_file, contenthandler=far.RankingStat.LIGOLWContentHandler))

fapfar = far.FAPFAR(rankingstatpdf.new_with_extinction())

livetime = fapfar.livetime
print "%s: livetime used in making IFAR plots\n" % (livetime,)

connection = sqlite3.connect(database[0])
xmldoc = dbtables.get_xml(connection)

for name, in set(connection.cursor().execute("SELECT name from segment_definer")):
	segs = ligolw_segments.LigolwSegments(xmldoc).get_by_name(str(name)).coalesce()
	print "%s: livetime %s" % (name, abs(segs))
	print "%s: two or more detectors on: %s\n" % (name, abs(segmentsUtils.vote(segs.values(), 2)))

datasegs = ligolw_segments.LigolwSegments(xmldoc).get_by_name(str("datasegments")).coalesce()
vetosegs = ligolw_segments.LigolwSegments(xmldoc).get_by_name(str("vetoes")).coalesce()

coinclivetime=abs(segmentsUtils.vote((datasegs-vetosegs).values(),2))

print "analyzable time (datasegments-vetosegments): %s" % (coinclivetime,)

for end_time, end_time_ns, combined_far, false_alarm_rate in connection.cursor().execute("SELECT end_time, end_time_ns, combined_far, false_alarm_rate from coinc_inspiral JOIN coinc_event on coinc_event.coinc_event_id == coinc_inspiral.coinc_event_id WHERE NOT EXISTS(SELECT * FROM time_slide WHERE time_slide.time_slide_id == coinc_event.time_slide_id AND time_slide.offset != 0) ORDER BY false_alarm_rate ASC LIMIT 5"):
	print "gps=%s: gstlalFAP=%s, gstlalFAR(Hz)=%.7E, gstlalFAR(/yr)=%.7E, coincltFAR(Hz)=%.7E, coincltFAR(/yr)=%.7E" % (end_time + end_time_ns * 1e-9, false_alarm_rate, float(combined_far), float(combined_far) * 86400 * 365.25, float(combined_far) * float(livetime) / float(coinclivetime), float(combined_far) * float(livetime) / float(coinclivetime) * 86400. * 365.25,)