Skip to content
Snippets Groups Projects

Plot horizon distance from ranking statistics

Merged Plot horizon distance from ranking statistics
plot_psd_horizon into master
1 unresolved thread
Merged ChiWai Chan requested to merge plot_psd_horizon into master
1 unresolved thread
Viewing commit bf15b904
Show latest version
1 file
+ 494
− 494
Compare changes
  • Side-by-side
  • Inline
gstlal-inspiral/bin/gstlal_inspiral_plot_sensitivity
+ 494
− 494
@@ -50,345 +50,345 @@ __date__ = "" # FIXME
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)
"""
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]
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])])
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 source_type_distance_duration_bins_from_sims(sims, distbins = 200):
"""
Given a list of the injections, guess at the template duration and
distance bins.
"""
dist_mchirp_vals = map(imr_utils.sim_to_distance_chirp_mass_bins_function, sims)
"""
Given a list of the injections, guess at the template duration 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]
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., 2., 10., numpy.inf])])
else:
return rate.NDBins([rate.LinearBins(min(distances), max(distances), distbins), rate.IrregularBins([0., 2., 10., numpy.inf])])
if min(distances) == max(distances):
return rate.NDBins([rate.LinearBins(min(distances)-1.0, max(distances)+1.0, 1), rate.IrregularBins([0., 2., 10., numpy.inf])])
else:
return rate.NDBins([rate.LinearBins(min(distances), max(distances), distbins), rate.IrregularBins([0., 2., 10., numpy.inf])])
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)
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(f"Gathering search data from zero-lag database: {f}....", file=sys.stderr)
with dbtables.workingcopy(f, tmp_path=opts.tmp_space, discard = True, verbose = opts.verbose) as working_filename:
connection = sqlite3.connect(str(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()
# 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(f"Reading results of injection analysis from {f} ...", file=sys.stderr)
with dbtables.workingcopy(f, tmp_path=opts.tmp_space, discard = True, verbose = opts.verbose) as working_filename:
connection = sqlite3.connect(str(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
# 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("\nInstruments: %s" % instr, file=sys.stdout)
print("Number of found injections: %d" % len(self.found_injections_fars[instr]), file=sys.stdout)
print("Number of total injections: %d" % len(self.total_injections[instr]), file=sys.stdout)
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("\nForming coincident segments from single IFO segments...", file=sys.stdout)
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("\t%s were on for %d seconds (after vetoes)" % (','.join(sorted(list(inst))),abs(self.zero_lag_segments[inst])), file=sys.stdout)
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
if bin_type == "Duration":
self.bins = source_type_distance_duration_bins_from_sims(self.total_injections[inst], distbins = opts.dist_bins)
self.sim_to_bins = imr_utils.sim_to_distance_duration_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
"""
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(f"Gathering search data from zero-lag database: {f}....", file=sys.stderr)
with dbtables.workingcopy(f, tmp_path=opts.tmp_space, discard = True, verbose = opts.verbose) as working_filename:
connection = sqlite3.connect(str(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()
# 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(f"Reading results of injection analysis from {f} ...", file=sys.stderr)
with dbtables.workingcopy(f, tmp_path=opts.tmp_space, discard = True, verbose = opts.verbose) as working_filename:
connection = sqlite3.connect(str(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
# 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("\nInstruments: %s" % instr, file=sys.stdout)
print("Number of found injections: %d" % len(self.found_injections_fars[instr]), file=sys.stdout)
print("Number of total injections: %d" % len(self.total_injections[instr]), file=sys.stdout)
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("\nForming coincident segments from single IFO segments...", file=sys.stdout)
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("\t%s were on for %d seconds (after vetoes)" % (','.join(sorted(list(inst))),abs(self.zero_lag_segments[inst])), file=sys.stdout)
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
if bin_type == "Duration":
self.bins = source_type_distance_duration_bins_from_sims(self.total_injections[inst], distbins = opts.dist_bins)
self.sim_to_bins = imr_utils.sim_to_distance_duration_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:
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.")
# binning options
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.")
parser.add_option("--bin-by-total-mass", default = False, action = "store_true", help = "Bin injections by total mass when estimating the search efficiency.")
parser.add_option("--bin-by-chirp-mass", default = False, action = "store_true", help = "Bin injections by chirp mass when estimating the search efficiency.")
parser.add_option("--bin-by-aligned-spin", default = False, action = "store_true", help = "Bin injections by aligned spin parameter when estimating the search efficiency.")
parser.add_option("--bin-by-mass-ratio", default = False, action = "store_true", help = "Bin injections by mass ratio when estimating the search efficiency.")
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")
parser.add_option("--bin-by-duration", default = False, action = "store_true", help = "Bin injections by template duration 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.bin_by_source_type:
opts.bin_types.append("Source_Type")
if opts.bin_by_duration:
opts.bin_types.append("Duration")
if opts.max_mass_ratio and (opts.max_mass_ratio < 1):
raise ValueError("The maximum mass ratio must be >=1!")
return opts, filenames
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.")
# binning options
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.")
parser.add_option("--bin-by-total-mass", default = False, action = "store_true", help = "Bin injections by total mass when estimating the search efficiency.")
parser.add_option("--bin-by-chirp-mass", default = False, action = "store_true", help = "Bin injections by chirp mass when estimating the search efficiency.")
parser.add_option("--bin-by-aligned-spin", default = False, action = "store_true", help = "Bin injections by aligned spin parameter when estimating the search efficiency.")
parser.add_option("--bin-by-mass-ratio", default = False, action = "store_true", help = "Bin injections by mass ratio when estimating the search efficiency.")
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")
parser.add_option("--bin-by-duration", default = False, action = "store_true", help = "Bin injections by template duration 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.bin_by_source_type:
opts.bin_types.append("Source_Type")
if opts.bin_by_duration:
opts.bin_types.append("Duration")
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 #####################################
@@ -400,13 +400,13 @@ 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
# 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
@@ -434,173 +434,173 @@ snrs = numpy.logspace(numpy.log10(opts.min_snr), numpy.log10(opts.max_snr), opts
#
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("No injections performed in %s time. Skipping..." % "".join(sorted(list(instr))), file=sys.stderr)
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("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)))), file=sys.stderr)
print("Skipping ...", file=sys.stderr)
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)) )
#
# plot the volume/range versus far/snr for each bin
# FIXME: since we are now writing the data to disk,
# it may make more sense to carry out the plotting
# in a separate code
#
far_vols = (vols_lo_far, vols_far, vols_hi_far)
snr_vols = (vols_lo_snr, vols_snr, vols_hi_snr)
eff_fids = (eff_fid_lo, eff_fid, eff_fid_hi)
fig_far = plotsens.plot_sensitivity_vs_far(far_vols, fars, UL.livetime, instr, bins, bin_type)
fig_far_range = plotsens.plot_range_vs_far(far_vols, fars, UL.livetime, instr, bins, bin_type)
fig_snr = plotsens.plot_sensitivity_vs_snr(snr_vols, snrs, UL.livetime, instr, bins, bin_type)
fig_eff = plotsens.plot_fiducial_efficiency(eff_fids, opts.fiducial_far, instr, bins, bin_type)
# save and close figures
ifostr = "%s_%s" % (UL.searched_instruments, "".join(sorted(instr)))
for fig, plot_type in zip((fig_far, fig_far_range, fig_snr, fig_eff), ("VOLUME_VS_FAR", "RANGE_VS_FAR", "VOLUME_VS_SNR", "EFFICIENCY")):
plot_desc = f"GSTLAL_INSPIRAL_PLOT_SENSITIVITY_{plot_type}_{opts.user_tag}_BINNED_BY_{bin_type.upper()}"
tag = datafind.T050017_filename(ifostr, plot_desc, (UL.start_time, UL.end_time), ".png", path = opts.output_dir)
fig.savefig(tag)
pyplot.close(fig)
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("No injections performed in %s time. Skipping..." % "".join(sorted(list(instr))), file=sys.stderr)
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("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)))), file=sys.stderr)
print("Skipping ...", file=sys.stderr)
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)) )
#
# plot the volume/range versus far/snr for each bin
# FIXME: since we are now writing the data to disk,
# it may make more sense to carry out the plotting
# in a separate code
#
far_vols = (vols_lo_far, vols_far, vols_hi_far)
snr_vols = (vols_lo_snr, vols_snr, vols_hi_snr)
eff_fids = (eff_fid_lo, eff_fid, eff_fid_hi)
fig_far = plotsens.plot_sensitivity_vs_far(far_vols, fars, UL.livetime, instr, bins, bin_type)
fig_far_range = plotsens.plot_range_vs_far(far_vols, fars, UL.livetime, instr, bins, bin_type)
fig_snr = plotsens.plot_sensitivity_vs_snr(snr_vols, snrs, UL.livetime, instr, bins, bin_type)
fig_eff = plotsens.plot_fiducial_efficiency(eff_fids, opts.fiducial_far, instr, bins, bin_type)
# save and close figures
ifostr = "%s_%s" % (UL.searched_instruments, "".join(sorted(instr)))
for fig, plot_type in zip((fig_far, fig_far_range, fig_snr, fig_eff), ("VOLUME_VS_FAR", "RANGE_VS_FAR", "VOLUME_VS_SNR", "EFFICIENCY")):
plot_desc = f"GSTLAL_INSPIRAL_PLOT_SENSITIVITY_{plot_type}_{opts.user_tag}_BINNED_BY_{bin_type.upper()}"
tag = datafind.T050017_filename(ifostr, plot_desc, (UL.start_time, UL.end_time), ".png", path = opts.output_dir)
fig.savefig(tag)
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)
write_cache(cache_list, opts.output_cache)