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.

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

import math
import numpy

import sqlite3

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

import os
import sys
from optparse import OptionParser
import itertools

from gstlal import inspiral_pipe
from gstlal import plotutil

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 iterutils
from ligo import segments
from ligo.segments import utils as segmentsUtils

import lal
from lal import rate
from lal.utils import CacheEntry

from gstlal import imr_utils

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

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


def source_type_distance_chirp_mass_bins_from_sims(sims, distbins = 200):
    """
    Given a list of the injections, guess at the chirp mass and distance
    bins.
    """
    dist_mchirp_vals = map(imr_utils.sim_to_distance_chirp_mass_bins_function, sims)

    distances = [tup[0] for tup in dist_mchirp_vals]

    if min(distances) == max(distances):
        return rate.NDBins([rate.LinearBins(min(distances)-1.0, max(distances)+1.0, 1), rate.IrregularBins([0.5, 2.0, 4.5, 45., 450])])
    else:
        return rate.NDBins([rate.LinearBins(min(distances), max(distances), distbins), rate.IrregularBins([0.5, 2.0, 4.5, 45., 450])])
    #return rate.NDBins([rate.LinearBins(min(distances), max(distances), distbins), rate.IrregularBins([0.8, 1.74, 8.07, 14.92, 21.77, 100.0, 218.00, 450.0])])

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.instrument_combos = []
        self.segments = segments.segmentlistdict()
        self.zero_lag_segments = segments.segmentlistdict()
        self.start_time = None
        self.end_time = None
        self.searched_instruments = 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_segments(connection)                   #get single ifo segments with vetoes applied
            self.get_searched_instruments(connection)

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

        # derived quantities
        self.get_instrument_combos()
        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.instrument_combos:
                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")
                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")
                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.instrument_combos:
            self.total_injections[instr] = imr_utils.symmetrize_sims(self.total_injections[instr], "mass1", "mass2")
            self.total_injections[instr] = self.filter_injections(self.total_injections[instr])
            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])

            if self.opts.verbose:
                print >>sys.stdout,"\nInstruments: %s" % instr
                print >>sys.stdout,"Number of found injections: %d" % len(self.found_injections_fars[instr])
                print >>sys.stdout,"Number of total injections: %d" % len(self.total_injections[instr])


    def get_searched_instruments(self,connection):
        '''
        Retrieve the sets of instruments which were on during the search.
        '''
        for inst in connection.cursor().execute("""SELECT DISTINCT(ifos) FROM process WHERE ifos NOTNULL""").fetchall():
            if self.searched_instruments is None:
                self.searched_instruments = inst[0].replace(",","")
            else:
                if self.searched_instruments != inst[0].replace(",",""):
                    raise ValueError("Detected different searched instruments between databases")


    def get_instrument_combos(self, min_instruments = 1):
        instruments = self.segments.keys()
        combos = set()
        for i in range(min_instruments, len(instruments,) + 1):
            for choice in itertools.combinations(instruments, i):
                # NOTE the reference IFO is always the first in
                # alphabetical order for any given combination,
                # hence the sort here
                combos.add(frozenset(sorted(choice)))
        # NOTE add an "ALL" key which will be used to provide a summary plot of the entire analysis not broken down by livetime
        combos.add(frozenset(("ALL",)))
        self.instrument_combos = tuple(sorted(list(combos)))
        return self.instrument_combos


    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 + 1e-9 * start_time_ns) FROM segment').fetchone()[0] )
        end_time = int( connection.cursor().execute('SELECT MAX(end_time + 1e-9 * end_time_ns) 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

        # retrieve single-ifo segments
        self.segments += imr_utils.get_segments( connection, dbtables.get_xml(connection), opts.coinc_table_name, opts.live_time_program, opts.veto_segments_name, opts.data_segments_name)

        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.instrument_combos[:]:
            # intersect single ifo segments
            self.zero_lag_segments.setdefault(inst,segments.segmentlist())
            if inst == frozenset(("ALL",)):
                self.zero_lag_segments[inst] = self.segments.union(self.segments.keys())
            else:
                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)" % (','.join(sorted(list(inst))),abs(self.zero_lag_segments[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
        if bin_type == "Source_Type":
            self.bins = source_type_distance_chirp_mass_bins_from_sims(self.total_injections[inst], distbins = opts.dist_bins)
            self.sim_to_bins = imr_utils.sim_to_distance_chirp_mass_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.instrument_combos:
            newinjs = []
            for sim_params in sims:

                if type(sim_params) == tuple:
                    far, sim = sim_params
                else:
                    sim = sim_params

                # 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

                # Exclude precessing waveforms if requested
                if self.opts.remove_precession:
                    if any(s != 0.0 for s in (sim.spin1x, sim.spin1y, sim.spin2x, sim.spin2y)):
                        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
                if type(sim_params) == tuple:
                    newinjs.append((far,sim))
                else:
                    newinjs.append(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 = __version__, 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("--fiducial-far", default = 1.0e-2/lal.YRJUL_SI, type = "float", help = "Specify a fiducial FAR (Hz) for plotting efficiency")
    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("--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 = 3, 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 = 3, 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 = 50, 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 = 0.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 = float("Inf"), 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 = 0.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 = float("Inf"), 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.")
    parser.add_option("--remove-precession", action = "store_true", help = "Exclude precessing injections when plotting sensitivity.")

    # 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.")

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

    # bin by source type
    parser.add_option("--bin-by-source-type", default = False, action = "store_true", help = "Bin injections by source type: BNS mchirp between 0.8 and 2, NSBH between 2 and 4.5, BBH between 4.5 and 45.0, IMBH between 45.0 and 450.0")

    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.bin_by_source_type:
        opts.bin_types.append("Source_Type")

    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.instrument_combos:

        if opts.verbose:
            print "\n\nComputing sensitive volume for %s observation time binning by %s...\n" % ("".join(sorted(list(instr))),bin_type)
        # check for empty injection sets. these tend to crash
        # set_bins, so we just catch it here
        if not UL.total_injections[instr]:
            print >> sys.stderr, "No injections performed in %s time. Skipping..." % "".join(sorted(list(instr)))
            continue

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

        # check for sufficient number of injection sets
        nbins = 1.0
        for bin in bins:
            nbins *= len(bin)

        if len(UL.total_injections[instr]) / nbins <= 1:
            print >> sys.stderr, "Insufficient injections (%d) performed in %s time to compute sensitivity. Perform more injections or request fewer bins." % (len(UL.total_injections[instr]), "".join(sorted(list(instr))))
            print >> sys.stderr, "Skipping ..."
            continue

        if opts.verbose:
            print "\tInjections performed: %d\n" % (len(UL.total_injections[instr]))


        #
        # compute efficiency at some fiducial far
        #
        found_by_far = [s for f, s in UL.found_injections_fars[instr] if f < opts.fiducial_far]
        eff_fid_lo, eff_fid, eff_fid_hi = imr_utils.compute_search_efficiency_in_bins(found_by_far, UL.total_injections[instr], bins, s2b)

        #
        # 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.Array.build(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.Array.build(u"array", fars))
        child.appendChild(xml)

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

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

        #
        # compute volume by far and snr for all mass bins
        #
        vols_lo_far, vols_far, vols_hi_far = [], [], []
        vols_lo_snr, vols_snr, vols_hi_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
            #
            eff_lo, eff, eff_hi = imr_utils.compute_search_efficiency_in_bins(found_by_far, UL.total_injections[instr], bins, s2b)
            vol_lo = imr_utils.compute_search_volume(eff_lo)
            vol_lo.array *= UL.livetime[instr] #preferred unit is Mpc^3*yr
            vols_lo_far.append(vol_lo)

            vol = imr_utils.compute_search_volume(eff)
            vol.array *= UL.livetime[instr] #preferred unit is Mpc^3*yr
            vols_far.append(vol)

            vol_hi = imr_utils.compute_search_volume(eff_hi)
            vol_hi.array *= UL.livetime[instr] #preferred unit is Mpc^3*yr
            vols_hi_far.append(vol_hi)

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

            xml = ligolw.LIGO_LW({u"Name": u"volume_by_far_%d:%s" % (k, bin_type,)})
            xml.appendChild(ligolw_array.Array.build(u"array", vol.array))
            child.appendChild(xml)

            xml = ligolw.LIGO_LW({u"Name": u"volume_hi_by_far_%d:%s" % (k, bin_type,)})
            xml.appendChild(ligolw_array.Array.build(u"array", vol_hi.array))
            child.appendChild(xml)

            #
            # compute volume vs snr
            #
            eff_lo, eff, eff_hi = imr_utils.compute_search_efficiency_in_bins(found_by_snr, UL.total_injections[instr], bins, s2b)
            vol_lo = imr_utils.compute_search_volume(eff_lo)
            vol_lo.array *= UL.livetime[instr] #preferred unit is Mpc^3*yr
            vols_lo_snr.append(vol_lo)

            vol = imr_utils.compute_search_volume(eff)
            vol.array *= UL.livetime[instr] #preferred unit is Mpc^3*yr
            vols_snr.append(vol)

            vol_hi = imr_utils.compute_search_volume(eff_hi)
            vol_hi.array *= UL.livetime[instr] #preferred unit is Mpc^3*yr
            vols_hi_snr.append(vol_hi)

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

            xml = ligolw.LIGO_LW({u"Name": u"volume_by_snr_%d:%s" % (k, bin_type,)})
            xml.appendChild(ligolw_array.Array.build(u"array", vol.array))
            child.appendChild(xml)

            xml = ligolw.LIGO_LW({u"Name": u"volume_hi_by_snr_%d:%s" % (k, bin_type,)})
            xml.appendChild(ligolw_array.Array.build(u"array", vol_hi.array))
            child.appendChild(xml)

        #
        # finish xml document
        #
        output_filename = "%s_%s-SEARCH_VOLUME_BINNED_BY_%s_%s-%d-%d.xml" % (UL.searched_instruments, "".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, fig_far_range, fig_snr, fig_eff = pyplot.figure(), pyplot.figure(), pyplot.figure(), pyplot.figure()
        fig_far.set_size_inches((8., 8. / plotutil.golden_ratio))
        fig_far_range.set_size_inches((8., 8. / plotutil.golden_ratio))
        fig_snr.set_size_inches((8., 8. / plotutil.golden_ratio))
        fig_eff.set_size_inches((8., 8. / plotutil.golden_ratio))
        ax_far = fig_far.gca()
        ax_far_range = fig_far_range.gca()
        ax_snr = fig_snr.gca()
        ax_eff = fig_eff.gca()

        # 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_\mathrm{total} \in [%.2f, %.2f] \,\mathrm{M}_\odot$" % (mlo[0], mhi[0])
            if bin_type == "Chirp_Mass" or bin_type == "Source_Type":
                label = "$M_\mathrm{chirp} \in [%.2f, %.2f] \,\mathrm{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] \,\mathrm{M}_\odot$" % (mlo[0], mhi[0], mlo[1], mhi[1])
            labels.append(label)

            # NOTE create regular plots, and define log x,y scales below
            #      since otherwise, fill_between allocates too many blocks and crashes
            # sensitivity vs far
            center = numpy.array([v[mmid] for v in vols_far])
            if (center == 0).all():
                continue
            lo = numpy.array([v[mmid] for v in vols_lo_far])
            hi = numpy.array([v[mmid] for v in vols_hi_far])
            line, = ax_far.plot( fars, center, label=label, linewidth=2 )
            ax_far.fill_between( fars, lo, hi, alpha=0.5, color=line.get_color())

            # range vs far
            center = (center/(4*math.pi*UL.livetime[instr]/3))**(1./3)
            lo = (lo/(4*math.pi*UL.livetime[instr]/3))**(1./3)
            hi = (hi/(4*math.pi*UL.livetime[instr]/3))**(1./3)
            line, = ax_far_range.plot( fars, center, label=label, linewidth=2 )
            ax_far_range.fill_between( fars, lo, hi, alpha=0.5, color=line.get_color())

            # sensitivity vs snr
            center = numpy.array([v[mmid] for v in vols_snr])
            lo = numpy.array([v[mmid] for v in vols_lo_snr])
            hi = numpy.array([v[mmid] for v in vols_hi_snr])
            line, = ax_snr.plot( snrs, center, label=label, linewidth=2 )
            ax_snr.fill_between( snrs, lo, hi, alpha=0.5, color=line.get_color())

            # plot fiducial efficiency
            # # far vs snr
            ds = (bins[0].lower() + bins[0].upper() )/2
            eff = numpy.array([eff_fid[(d,) + mmid] for d in ds])
            lo = numpy.array([eff_fid_lo[(d,) + mmid] for d in ds])
            hi = numpy.array([eff_fid_hi[(d,) + mmid] for d in ds])
            line, = ax_eff.plot(ds, eff, label=label, linewidth=2 )
            ax_eff.fill_between(ds, lo, hi, alpha=0.5, color=line.get_color())

        # sensitivity vs far
        ax_far.set_xlabel("Combined FAR (Hz)")
        ax_far.set_ylabel(r"Volume $\times$ Time ($\mathrm{Mpc}^3 \mathrm{yr}$)")
        ax_far.set_xscale("log")
        ax_far.set_yscale("log")
        ax_far.set_xlim([min(fars), max(fars)])
        ax_far.invert_xaxis()
        ax_far.legend(loc="lower left")
        ax_far.grid()

        vol_tix = ax_far.get_yticks()
        tx = ax_far.twinx() # map volume to distance
        tx.set_yticks(vol_tix)
        tx.set_yscale("log")
        tx.set_ylim(ax_far.get_ylim())
        tx.set_yticklabels( ["%.3g" % ((float(k)/(4*math.pi*UL.livetime[instr]/3))**(1./3)) for k in vol_tix] )
        tx.set_ylabel("Range (Mpc)")

        # range vs far
        ax_far_range.set_xlabel("Combined FAR (Hz)")
        ax_far_range.set_ylabel("Range (Mpc)")
        ax_far_range.set_xscale("log")
        ax_far_range.set_xlim([min(fars), max(fars)])
        ax_far_range.invert_xaxis()
        ax_far_range.legend(loc="lower left")
        ax_far_range.grid()

        # sensitivity vs snr
        ax_snr.set_xlabel("Network SNR")
        ax_snr.set_ylabel(r"Volume $\times$ Time ($\mathrm{Mpc}^3 \mathrm{yr}$)")
        ax_snr.set_yscale("log")
        ax_snr.set_xlim([min(snrs), max(snrs)])
        ax_snr.set_ylim(ymin=0)
        ax_snr.legend(loc="lower left")
        ax_snr.grid()

        vol_tix = ax_snr.get_yticks()
        tx = ax_snr.twinx() # map volume to distance
        tx.set_yticks(vol_tix)
        tx.set_yscale("log")
        tx.set_ylim(ax_snr.get_ylim())
        tx.set_yticklabels( ["%.3g" % ((float(k)/(4*math.pi*UL.livetime[instr]/3))**(1./3)) for k in vol_tix] )
        tx.set_ylabel("Range (Mpc)")

        # efficiency vs distance
        ax_eff.set_xlabel("Distance (Mpc)")
        ax_eff.set_ylabel("Efficiency")
        ax_eff.set_xscale("log")
        ax_eff.legend(loc="upper right")
        ax_eff.grid()

        ax_snr.set_title("%s Observing (%.2f days)" % ("".join(sorted(list(instr))), UL.livetime[instr]*365.25))
        ax_far.set_title("%s Observing (%.2f days)" % ("".join(sorted(list(instr))), UL.livetime[instr]*365.25))
        ax_far_range.set_title("%s Observing (%.2f days)" % ("".join(sorted(list(instr))), UL.livetime[instr]*365.25))
        ax_eff.set_title("%s Observing ($\mathrm{FAR} < %s\,\mathrm{Hz}$)" % ("".join(sorted(list(instr))), plotutil.latexnumber("%.2e" % opts.fiducial_far)))

        fig_far.tight_layout(pad = .8)
        fig_far_range.tight_layout(pad = .8)
        fig_snr.tight_layout(pad = .8)
        fig_eff.tight_layout(pad = .8)

        # save and close figures
        ifostr = "%s_%s" % (UL.searched_instruments, "".join(sorted(instr)))

        tag = inspiral_pipe.T050017_filename(ifostr, "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.savefig(tag)
        pyplot.close(fig_far)

        tag = inspiral_pipe.T050017_filename(ifostr, "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)
        pyplot.close(fig_far_range)

        tag = inspiral_pipe.T050017_filename(ifostr, "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.savefig(tag)
        pyplot.close(fig_snr)

        tag = inspiral_pipe.T050017_filename(ifostr, "GSTLAL_INSPIRAL_PLOT_SENSITIVITY_%s_EFFICIENCY_BINNED_BY_%s" % (opts.user_tag, bin_type.upper()), (UL.start_time, UL.end_time), ".png", path = opts.output_dir)
        fig_eff.savefig(tag)
        pyplot.close(fig_eff)


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