gstlal_inspiral_plot_sensitivity

#!/usr/bin/python
#
# Copyright (C) 2012-2014  Stephen Privitera, Kipp Cannon, Chad Hanna, Drew Keppel
#
# 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.

## @file
# A program to plot the sensitivity of a gstlal inspiral search, e.g., VT

import numpy

try:
  import sqlite3
except ImportError:
  # pre 2.5.x
  from pysqlite2 import dbapi2 as sqlite3

import matplotlib
matplotlib.use("agg")
from matplotlib import pyplot

import os
import sys
import copy
from optparse import OptionParser

from gstlal import inspiral_pipe

from glue import segments
from glue.ligolw import ligolw
from glue.ligolw import lsctables
from glue.ligolw import dbtables
from glue.ligolw import utils as ligolw_utils
from glue.ligolw import array as ligolw_array
from glue import segmentsUtils
from glue.lal import CacheEntry
from glue import iterutils

import lal

from pylal import rate
from pylal import imr_utils
from pylal import git_version

__author__ = "Stephen Privitera <sprivite@caltech.edu>, Chad Hanna <channa@perimeterinstitute.ca>, Kipp Cannon <kipp.cannon@ligo.org>"
__version__ = "git id %s" % git_version.id
__date__ = git_version.date

pyplot.rc('font',**{'family':'serif','serif':['Computer Modern Roman']})
matplotlib.rcParams.update({"text.usetex": True})


def chirp_mass(m1,m2):
    m1 = numpy.array(m1)
    m2 = numpy.array(m2)
    mu = (m1*m2)/(m1+m2)
    mtotal = m1+m2
    return mu**(3./5) *mtotal**(2./5)

class upper_limit(object):
    """
    The upper_limit class organizes the calculation of the sensitive search volume
    for a search described by the input database.
    """
    def __init__(self, opts):
        ## Instance variables ######################################
        self.opts = opts
        self.instruments = []
        self.segments = segments.segmentlistdict()
        self.zero_lag_segments = segments.segmentlistdict()
        self.start_time = None
        self.end_time = None
        ############################################################

        # read the zero lag databases first
        for f in self.opts.zero_lag_database:
            if opts.verbose:
                print >> sys.stdout, "\nGathering search data from zero-lag database: %s...." % (f,)
            working_filename = dbtables.get_connection_filename(f, tmp_path=opts.tmp_space, verbose = opts.verbose)
            connection = sqlite3.connect(working_filename)

            # find out which instruments were on and when
            self.get_on_instruments(connection)
            self.get_segments(connection)                   #get single ifo segments with vetoes applied

            # done with zl database
            connection.commit()
            dbtables.discard_connection_filename(f, working_filename, verbose = self.opts.verbose)

        # derived quantities
        self.get_zero_lag_segments()          #make coincident ifo segments from single ifo segments
        self.get_livetime()                   #get the livetime for the chosen set of instruments

        # open up injection databases
        self.found_injections_fars = {}
        self.found_injections_snrs = {}
        self.total_injections = {}
        for f in self.opts.injection_database:
            if opts.verbose:
                print "Reading results of injection analysis from %s ..."%f
            working_filename = dbtables.get_connection_filename(f, tmp_path=opts.tmp_space, verbose = opts.verbose)
            connection = sqlite3.connect(working_filename)
            for inst in self.instruments:
                self.found_injections_fars.setdefault(inst,[])
                self.found_injections_snrs.setdefault(inst,[])
                self.total_injections.setdefault(inst,[])

                # get sim-far relationships
                found, total, missed = imr_utils.get_min_far_inspiral_injections(connection, segments = self.zero_lag_segments[inst], table_name="coinc_inspiral:table")
                self.found_injections_fars[inst] += found
                self.total_injections[inst] += total

                # get sim-snr relationships
                found, total, missed = imr_utils.get_max_snr_inspiral_injections(connection, segments = self.zero_lag_segments[inst], table_name="coinc_inspiral:table")
                self.found_injections_snrs[inst] += found

            dbtables.discard_connection_filename(f, working_filename, verbose = self.opts.verbose)

        # remove injections that violate user-imposed constraints and
        # symmetrize between m1-m2
        for instr in self.instruments:
            self.total_injections[instr] = imr_utils.symmetrize_sims(self.total_injections[instr], "mass1", "mass2")
            self.found_injections_fars[instr] = self.filter_injections(self.found_injections_fars[instr])
            self.found_injections_snrs[instr] = self.filter_injections(self.found_injections_snrs[instr])


    def get_on_instruments(self,connection):
        '''
        Retrieve the sets of instruments which were on during the search.
        '''
    for inst in connection.cursor().execute("""SELECT DISTINCT(ifos) FROM search_summary"""):
            inst =  frozenset(lsctables.instrument_set_from_ifos(inst[0]))
            if not inst in self.instruments:
                self.instruments.append(inst)

        return self.instruments


    def get_segments(self,connection):
        '''
        Retrieve raw single IFO segments from the database and
        subtract vetoes.
        '''
        #
        # get start and end times of the analysis, used for output
        # naming. these may extend beyond the actual analyzed
        # segments because of vetoes.
        #
        start_time = int( connection.cursor().execute('SELECT MIN(start_time) FROM segment').fetchone()[0] )
        end_time = int( connection.cursor().execute('SELECT MAX(end_time) FROM segment').fetchone()[0] )
        if not self.start_time or start_time < self.start_time:
            self.start_time = start_time
        if not self.end_time or end_time < self.end_time:
            self.end_time = end_time


        # get single ifos from ifo sets
        ifos = set()
        for inst in self.instruments:
            ifos |= inst

        # retrieve single-ifo segments
        segs = imr_utils.get_segments( connection, dbtables.get_xml(connection), "%s:table" % opts.coinc_table_name, opts.live_time_program, opts.veto_segments_name, opts.data_segments_name)

        for ifo in ifos:
            self.segments.setdefault(ifo,segments.segmentlist())
            self.segments[ifo] |=  segs[ifo]

        return self.segments


    def get_zero_lag_segments(self):
        '''
        Compute multi-ifo (coincident) segment list from single ifo segments.
        '''
        if self.opts.verbose:
            print >>sys.stdout,"\nForming coincident segments from single IFO segments..."

        for inst in self.instruments[:]:
            # intersect single ifo segments
            self.zero_lag_segments.setdefault(inst,segments.segmentlist())
            self.segments.union(set(self.segments.keys()) - inst)
            self.zero_lag_segments[inst] = self.segments.intersection(inst) - self.segments.union(set(self.segments.keys()) - inst)
            if self.opts.verbose:
                print >>sys.stdout,"\t%s were on for %d seconds (after vetoes, including playground)" % (','.join(sorted(list(inst))),abs(self.zero_lag_segments[inst]))

            # subtract playground segments
            if not self.opts.include_play:
                self.zero_lag_segments[inst] -= segmentsUtils.S2playground(self.zero_lag_segments[inst].extent())
                if self.opts.verbose:
                    print >>sys.stdout,"\t%s were on for %d seconds (after vetoes, excluding playground)" % (','.join(sorted(list(inst))),abs(self.zero_lag_segments[inst]))

            # remove instrument sets that were never on
            if abs(self.zero_lag_segments[inst]) == 0:
                print >>sys.stderr, "No livetime for in %s observation time. Skipping..."%("".join(sorted(list(inst))))
                self.instruments.remove(inst)

        return self.zero_lag_segments


    def get_livetime(self):
        '''Compute the instrument livetimes from the search segments.'''
        self.livetime = {}
        for inst in self.zero_lag_segments.keys():
            self.livetime[inst] = float(abs(self.zero_lag_segments[inst]))/lal.YRJUL_SI

        return True


    def set_bins(self, bin_type, inst):
        if bin_type == "Total_Mass":
            self.bins = imr_utils.guess_distance_total_mass_bins_from_sims(self.total_injections[inst], nbins = opts.mass_bins, distbins = opts.dist_bins)
            self.sim_to_bins = imr_utils.sim_to_distance_total_mass_bins_function
        if bin_type == "Chirp_Mass":
            self.bins = imr_utils.guess_distance_chirp_mass_bins_from_sims(self.total_injections[inst], mbins = opts.mass_bins, distbins = opts.dist_bins)
            self.sim_to_bins = imr_utils.sim_to_distance_chirp_mass_bins_function
        if bin_type == "Mass1_Mass2":
            self.bins = imr_utils.guess_distance_mass1_mass2_bins_from_sims(self.total_injections[inst], mass1bins = opts.mass_bins, mass2bins = opts.mass_bins, distbins = opts.dist_bins)
            self.sim_to_bins = imr_utils.sim_to_distance_mass1_mass2_bins_function
        if bin_type == "Aligned_Spin":
            self.bins = imr_utils.guess_distance_effective_spin_parameter_bins_from_sims(self.total_injections[inst], chibins = opts.spin_bins, distbins = opts.dist_bins)
            self.sim_to_bins = imr_utils.sim_to_distance_effective_spin_parameter_bins_function
        if bin_type == "Mass_Ratio":
            self.bins = imr_utils.guess_distance_mass_ratio_bins_from_sims(self.total_injections[inst], qbins = opts.mass_bins, distbins = opts.dist_bins)
            self.sim_to_bins = imr_utils.sim_to_distance_mass_ratio_bins_function

        return self.bins, self.sim_to_bins


    def filter_injections(self, sims):
        '''
        Remove injections from those found in the database based on
        user-imposed restrictions.
        '''
        for inst in self.instruments:
            newinjs = []
            for far, sim in sims:
                # exclude simulations from certain waveform families
                if sim.waveform in self.opts.exclude_sim_type:
                    continue

                if not self.opts.min_mtotal < sim.mass1 + sim.mass2 < self.opts.max_mtotal:
                    continue

                if self.opts.max_mass_ratio and not (1.0/self.opts.max_mass_ratio <= sim.mass1/sim.mass2 <= self.opts.max_mass_ratio):
                    continue

                # symmetrize in m1-m2
                if sim.mass1 > sim.mass2:
                    m1 = sim.mass1
                    m2 = sim.mass2
                    sim.mass1 = m2
                    sim.mass2 = m1

                # passed all constraint tests, we can use it
                newinjs.append((far,sim))

        return newinjs


def parse_command_line():

    description = '''
   description:

The program gstlal_inspiral_plot_sensitivity computes the sensitive volume of a CBC search from input databases containing triggers from simulation experiments. These triggers need to be ranked by false alarm rate, the detection statistic used in S6 searches. Then injections which register a trigger louder than the loudest event, by false alarm rate, are considered found. All others are considered missed. The efficiency of detecting an event depends on the source parameters, such as its component masses, distance, spin, inclination, sky position, etc. However, lalapps_cbc_svim only considers the dependency of the efficiency on distance and mass, marginalizing over the other parameters. Injections are binned in distance and mass and the estimated efficiency is integrated over distance to convert the efficiency into a physical volume.
'''

    parser = OptionParser(version = git_version.verbose_msg, usage = description)

    # FAR range and resolution
    parser.add_option("--xaxis-points", default = 50, type = "int", help = "Specify the number of FARs/SNRs for which to compute the search volume.")
    parser.add_option("--min-far", default = 1.0e-6/lal.YRJUL_SI, type = "float", help = "Specify the minimum FAR (Hz)") # one per million years is probably detection worthy
    parser.add_option("--max-far", default = 12.0/lal.YRJUL_SI, type = "float", help = "Specify the maximum FAR (Hz)") # one per month is possibly EM-followup worthy
    parser.add_option("--min-snr", default = 7, type = "float", help = "Specify the minimum SNR for volume calculation.")
    parser.add_option("--max-snr", default = 15, type = "float", help = "Specify the maximum SNR for volume calculation.")

    # Input data options
    parser.add_option("--include-play", default = False, action = "store_true", help = "Include playground data in computing the livetime and volume.")
    parser.add_option("--zero-lag-database", default = [], action = "append", help = "Name of database containing the zero lag segments and triggers.")
    parser.add_option("--injection-database", default = [], action = "append", help = "Name of database containing injection parameters and triggers.")
    parser.add_option("--live-time-program", default = "gstlal_inspiral", metavar = "name", help = "Set the name of the program whose rings will be extracted from the search_summary table.  Default = \"gstlal_inspiral\".")
    parser.add_option("--veto-segments-name", default = "vetoes", help = "Set the name of the veto segments to use from the XML document.")
    parser.add_option("--data-segments-name", default = "whitehtsegments", help = "Set the name of the data segments. Default = \"whitehtsegments\".")
    parser.add_option("--coinc-table-name", default = "coinc_inspiral", metavar = "name", help = "Set the name of the table containing coincident triggers.  Default = \"coinc_inspiral\".")

    # Output data options
    parser.add_option("--user-tag", default = "SEARCH", metavar = "name", help = "Add a descriptive tag to the names of output files.")
    parser.add_option("--output-dir", default = "./", metavar = "name", help = "Select a directory to place output files.")
    parser.add_option("-t", "--tmp-space", metavar = "path", help = "Path to a directory suitable for use as a work area while manipulating the database file.  The database file will be worked on in this directory, and then moved to the final location when complete.  This option is intended to improve performance when running in a networked environment, where there might be a local disk with higher bandwidth than is available to the filesystem on which the final output will reside.")
    parser.add_option("--output-path", default = "./", action = "store", help="Choose directory to save output files.")
    parser.add_option("--output-cache", default = None, help = "Name of output cache file. If not specified, then no cache file will be written.")
    parser.add_option("--verbose", action = "store_true", help = "Be verbose.")

    #
    # Binning options
    #
    parser.add_option("--mass-bins", default = 4, metavar = "integer", type = "int", help = "Number of mass bins (per dimension), if mass bin boundaries are not explicitly set.") 
    parser.add_option("--spin-bins", default = 4, metavar = "integer", type = "int", help = "Number of spin bins (per dimension), if spin bin boundaries are not explicitly set.")
    parser.add_option("--dist-bins", default = 15, metavar = "integer", type = "int", help = "Space distance bins evenly and specify the number of distance bins to use.")

    # injection filters
    parser.add_option("--min-mtotal", metavar = "m", type = 'float', default = 2.0, help = "Specify the minimum total mass to consider among the injections found in the DB. This filters all injections outside this total mass range.")
    parser.add_option("--max-mtotal", metavar = "m", type = 'float', default = 100.0, help = "Specify the maximum total mass to consider among the injections found in the DB. This filters all injections outside this total mass range.")
    parser.add_option("--min-mass", metavar = "MM", type = 'float', default = 1.0, help = "Specify the minimum component mass to consider among the injections found in the DB. This filters all injectionswith any component lighter than MM.")
    parser.add_option("--max-mass", metavar = "MM", type = 'float', default = 99.0, help = "Specify the maximum mass to consider among the injections found in the DB. This filters all injections with any component heavier than MM.")
    parser.add_option("--max-mass-ratio", metavar = "q", type = 'float', default = None, help = "Specify the maximum allowed mass ratio. Should be >= 1.")
    parser.add_option("--exclude-sim-type", default = [], action = "append", metavar = "SIM", help = "When computing the search volume, exclude injections made using the SIM waveform family. Example: SpinTaylorthreePointFivePN. Use this option multiple times to exclude more than one injection type.")

    # more options
    # Bin injections in mass by m1-m2
    parser.add_option("--bin-by-mass1-mass2", default = False, action = "store_true", help = "Bin injections by component mass in two dimensions when estimating the search efficiency.")

    # Bin injections by total mass
    parser.add_option("--bin-by-total-mass", default = False, action = "store_true", help = "Bin injections by total mass when estimating the search efficiency.")

    # Bin injections by chirp mass
    parser.add_option("--bin-by-chirp-mass", default = False, action = "store_true", help = "Bin injections by chirp mass when estimating the search efficiency.")

    # aligned spin
    parser.add_option("--bin-by-aligned-spin", default = False, action = "store_true", help = "Bin injections by aligned spin parameter when estimating the search efficiency.")

    # aligned spin
    parser.add_option("--bin-by-mass-ratio", default = False, action = "store_true", help = "Bin injections by mass ratio when estimating the search efficiency.")

    opts, filenames = parser.parse_args()
    opts.injection_database.extend(filenames)

    opts.bin_types = []
    if opts.bin_by_total_mass:
        opts.bin_types.append("Total_Mass")
    if opts.bin_by_chirp_mass:
        opts.bin_types.append("Chirp_Mass")
    if opts.bin_by_mass1_mass2:
        opts.bin_types.append("Mass1_Mass2")
    if opts.bin_by_aligned_spin:
        opts.bin_types.append("Aligned_Spin")
    if opts.bin_by_mass_ratio:
        opts.bin_types.append("Mass_Ratio")

    if opts.max_mass_ratio and (opts.max_mass_ratio < 1):
        raise ValueError, "The maximum mass ratio must be >=1!"

    return opts, filenames


############################ MAIN PROGRAM #####################################
###############################################################################
###############################################################################

#create an empty cache which will store the output files/metadata
cache_list = []


def write_cache(cache_list, fileout):
    # write cache file
    if opts.output_cache is not None:
        fd = open( fileout, "w" )
        for l in cache_list:
            fd.write( str(l) + "\n")
        fd.close()
    return

#
# MAIN
#


# read in command line opts
opts, database = parse_command_line()

# read in data from input database, store in upper limit object
if opts.verbose: print "\n\nSetting up the search volume calculations...\n"
UL = upper_limit(opts)

# get the range of xaxis values to use
fars = numpy.logspace(numpy.log10(opts.min_far), numpy.log10(opts.max_far), opts.xaxis_points)
snrs = numpy.logspace(numpy.log10(opts.min_snr), numpy.log10(opts.max_snr), opts.xaxis_points)


#
# loop over the requested instruments and mass bin types, compute the
# search volume as a function of FAR and SNR threshold
#
for bin_type in opts.bin_types:

    for instr in UL.instruments:

        # check for empty injection sets
        if not UL.total_injections[instr]:
            print >> sys.stderr, "No injections performed in %s time. Skipping..." % "".join(sorted(list(instr)))
            continue

        if opts.verbose:
            print "\n\nComputing sensitive volume for %s observation time binning by %s...\n" % ("".join(sorted(list(instr))),bin_type)

        bins, s2b = UL.set_bins(bin_type,instr)

        #
        # prepare output XML file. record mass bins, fars, snrs and
        # livetime
        #
        xmldoc = ligolw.Document()
        child = xmldoc.appendChild(ligolw.LIGO_LW())

        # write out mass bins
        for j in range(len(bins[1:])):
            xml = ligolw.LIGO_LW({u"Name": u"mass%d_bins:%s" % (j+1,bin_type)})
            edges = tuple( numpy.concatenate((l,u[-1:])) for l,u in zip(bins.lower(), bins.upper()) )
            xml.appendChild(ligolw_array.from_array(u"array", edges[j+1]))
            child.appendChild(xml)

        # write out fars/snrs/livetime
        xml = ligolw.LIGO_LW({u"Name": u"far_bins:%s" % (bin_type,)})
        xml.appendChild(ligolw_array.from_array(u"array", fars))
        child.appendChild(xml)

        xml = ligolw.LIGO_LW({u"Name": u"snr_bins:%s" % (bin_type,)})
        xml.appendChild(ligolw_array.from_array(u"array", snrs))
        child.appendChild(xml)

        xml = ligolw.LIGO_LW({u"Name": u"livetime:%s" % (bin_type,)})
        xml.appendChild(ligolw_array.from_array(u"array", numpy.array([UL.livetime[instr]])))
        child.appendChild(xml)

        #
        # compute volume by far and snr for all mass bins
        #
        vols_far, errs_far, vols_snr, errs_snr = [], [], [], []
        for k, far, snr in zip(range(opts.xaxis_points), fars, snrs):

            #
            # get found/missed injections
            #
            found_by_far = [s for f, s in UL.found_injections_fars[instr] if f < far]
            found_by_snr = [s for rho, s in UL.found_injections_snrs[instr] if rho > snr]

            #
            # compute volume vs far
            #
            vol, err = imr_utils.compute_search_volume_in_bins(found_by_far, UL.total_injections[instr], bins, s2b)
            vol.array *= UL.livetime[instr] #preferred unit is Mpc^3*yr
            err.array *= UL.livetime[instr] #preferred unit is Mpc^3*yr
            vols_far.append(vol)
            errs_far.append(err)

            #
            # write volume and volume errors array to xml
            #
            xml = ligolw.LIGO_LW({u"Name": u"volume_by_far_%d:%s" % (k, bin_type,)})
            xml.appendChild(ligolw_array.from_array(u"array", vol.array))
            child.appendChild(xml)

            xml = ligolw.LIGO_LW({u"Name": u"volume_error_by_far_%d:%s" % (k, bin_type,)})
            xml.appendChild(ligolw_array.from_array(u"array", err.array))
            child.appendChild(xml)

            #
            # compute volume vs snr
            #
            vol, err = imr_utils.compute_search_volume_in_bins(found_by_snr, UL.total_injections[instr], bins, s2b)
            vol.array *= UL.livetime[instr] #preferred unit is Mpc^3*yr
            err.array *= UL.livetime[instr] #preferred unit is Mpc^3*yr
            vols_snr.append(vol)
            errs_snr.append(err)

            #
            # write volume and volume errors array to xml
            #
            xml = ligolw.LIGO_LW({u"Name": u"volume_by_snr_%d:%s" % (k, bin_type,)})
            xml.appendChild(ligolw_array.from_array(u"array", vol.array))
            child.appendChild(xml)

            xml = ligolw.LIGO_LW({u"Name": u"volume_error_by_snr_%d:%s" % (k, bin_type,)})
            xml.appendChild(ligolw_array.from_array(u"array", err.array))
            child.appendChild(xml)

        #
        # finish xml document
        #
        output_filename = "%s-SEARCH_VOLUME_BINNED_BY_%s_%s-%d-%d.xml" % ("".join(sorted(list(instr))), bin_type.upper(), opts.user_tag, UL.start_time, UL.end_time-UL.start_time)
        ligolw_utils.write_filename(xmldoc, output_filename)
        cache_list.append( CacheEntry( "".join(sorted(list(instr))), bin_type,segments.segment(UL.start_time, UL.end_time), "file://localhost%s/%s" % (os.getcwd(), output_filename)) )

        #
        # set up figures for plots FIXME: since we are now writing the
        # data to disk, it may make more sense to carry out the
        # plotting in a separate code
        #
        fig_far_vol, fig_far_range, fig_snr_vol, fig_snr_range = pyplot.figure(), pyplot.figure(), pyplot.figure(), pyplot.figure()
        figs = [fig_far_vol, fig_far_range, fig_snr_vol, fig_snr_range]

        # plot the volume/range versus far/snr for each bin
        mbins = rate.NDBins(bins[1:])
        labels = []
        for mlo, mmid, mhi in zip(iterutils.MultiIter(*mbins.lower()),
                                  iterutils.MultiIter(*mbins.centres()),
                                  iterutils.MultiIter(*mbins.upper())):

            # plot labels
            if bin_type == "Aligned_Spin":
                label = "$\chi \in [%.2f, %.2f]$" % (mlo[0], mhi[0])
            if bin_type == "Mass_Ratio":
                label = "$m_1/m_2 \in [%.2f, %.2f]$" % (mlo[0], mhi[0])
            if bin_type == "Total_Mass":
                label = "$M_{total} \in [%.2f, %.2f] M_\odot$" % (mlo[0], mhi[0])
            if bin_type == "Chirp_Mass":
                label = "$M_{chirp} \in [%.2f, %.2f] M_\odot$" % (mlo[0], mhi[0])
            if bin_type == "Mass1_Mass2":
                if mmid[0] > mmid[1]: # symmetrized sims have m1 < m2
                    continue
                label = "$m_1 \in [%.2f, %.2f], m_2 \in [%.2f, %.2f] M_\odot$" % (mlo[0], mhi[0], mlo[1], mhi[1])
            labels.append(label)

            # volume plots
            fig_far_vol.gca().errorbar( fars, [v[mmid] for v in vols_far], yerr = [e[mmid] for e in errs_far], label=label )
            fig_snr_vol.gca().errorbar( snrs, [v[mmid] for v in vols_snr], yerr = [e[mmid] for e in errs_snr], label=label )

            # range plots
            dist = [ rate.BinnedArray(mbins,array = (v.array/UL.livetime[instr]/(4*numpy.pi/3))**(1./3) ) for v in vols_far ]
            derr = [ rate.BinnedArray(mbins,array = (d.array/3)*(e.array/v.array) ) for e,v,d in zip(errs_far, vols_far, dist) ]
            fig_far_range.gca().errorbar( fars, [d[mmid] for d in dist], yerr = [e[mmid] for e in derr], label=label )

            dist = [ rate.BinnedArray(mbins,array = (v.array/UL.livetime[instr]/(4*numpy.pi/3))**(1./3) ) for v in vols_snr ]
            derr = [ rate.BinnedArray(mbins,array = (d.array/3)*(e.array/v.array) ) for e,v,d in zip(errs_snr, vols_snr, dist) ]
            fig_snr_range.gca().errorbar( snrs, [d[mmid] for d in dist], yerr = [e[mmid] for e in derr], label=label )

        # customize axes
        ax = fig_far_vol.gca()
        ax.set_xlabel("Combined FAR (Hz)")
        ax.set_ylabel("VT (Mpc$^3$ yr)")
        ax.set_xscale("log")
        ax.set_xlim([min(fars), max(fars)])
        ax.invert_xaxis()
        ax.set_ylim(ymin=0)

        ax = fig_snr_vol.gca()
        ax.set_xlabel("Network SNR")
        ax.set_ylabel("VT (Mpc$^3$ yr)")
        ax.set_xlim([min(snrs), max(snrs)])
        ax.set_ylim(ymin=0)

        ax = fig_far_range.gca()
        ax.set_xlabel("Combined FAR (Hz)")
        ax.set_ylabel("Range (Mpc)")
        ax.set_xscale("log")
        ax.set_xlim([min(fars), max(fars)])
        ax.invert_xaxis()
        ax.set_ylim(ymin=0)

        ax = fig_snr_range.gca()
        ax.set_xlabel("Network SNR")
        ax.set_ylabel("Range (Mpc)")
        ax.set_xlim([min(snrs), max(snrs)])
        ax.set_ylim(ymin=0)

        # common to all figures
        for fig in figs:
            ax = fig.gca()
            ax.grid()
            ax.legend(loc="upper right")
            ax.set_title("Search Sensitivity in %s Time" % "".join(sorted(list(instr))))

        # save and close figures
        tag = inspiral_pipe.T050017_filename(instr, "GSTLAL_INSPIRAL_PLOT_SENSITIVITY_%s_VOLUME_VS_FAR_BINNED_BY_%s" % (opts.user_tag, bin_type.upper()), UL.start_time, UL.end_time, ".png", path = opts.output_dir)
        fig_far_vol.savefig(tag)

        tag = inspiral_pipe.T050017_filename(instr, "GSTLAL_INSPIRAL_PLOT_SENSITIVITY_%s_VOLUME_VS_SNR_BINNED_BY_%s" % (opts.user_tag, bin_type.upper()), UL.start_time, UL.end_time, ".png", path = opts.output_dir)
        fig_snr_vol.savefig(tag)

        tag = inspiral_pipe.T050017_filename(instr, "GSTLAL_INSPIRAL_PLOT_SENSITIVITY_%s_RANGE_VS_FAR_BINNED_BY_%s" % (opts.user_tag, bin_type.upper()), UL.start_time, UL.end_time, ".png", path = opts.output_dir)
        fig_far_range.savefig(tag)

        tag = inspiral_pipe.T050017_filename(instr, "GSTLAL_INSPIRAL_PLOT_SENSITIVITY_%s_RANGE_VS_SNR_BINNED_BY_%s" % (opts.user_tag, bin_type.upper()), UL.start_time, UL.end_time, ".png", path = opts.output_dir)
        fig_snr_range.savefig(tag)

        for fig in figs:
            pyplot.close(fig)

# write a cache file describing the files generated during by this program
if opts.output_cache:
    write_cache(cache_list, opts.output_cache)