svd_bank.py: move the make_horizon_distance function from gstlal_inspiral

gstlal-inspiral/python/svd_bank.py

@@ -58,6 +58,7 @@ Attributes = ligolw.sax.xmlreader.AttributesImpl
 from gstlal import cbc_template_fir
 from gstlal import misc as gstlalmisc
 from gstlal import templates
+from gstlal import reference_psd
 
 
 class DefaultContentHandler(ligolw.LIGOLWContentHandler):
@@ -439,3 +440,29 @@ def svdbank_templates_mapping(filenames, contenthandler, verbose = False):
 			print >>sys.stderr, "%d/%d:" % (n, len(filenames)),
 		mapping[filename] = sum((bank.sngl_inspiral_table for bank in read_banks(filename, contenthandler, verbose = verbose)), [])
 	return mapping
+
+
+def make_horizon_distance_func(banks):
+	"""
+	Takes a dictionary of objects returned by read_banks keyed by instrument
+	"""
+	# span is [15 Hz, 0.85 * Nyquist frequency]
+	# find the Nyquist frequency for the PSD to be used for each
+	# instrument.  require them to all match
+	sngl_inspiral_table = banks.values()[0][0].sngl_inspiral_table.copy()
+	for bank in banks.values()[0]:
+		sngl_inspiral_table.extend(bank.sngl_inspiral_table)
+	nyquists = set(max(rate for bank in banklist for rate in bank.get_rates()) // 2 for instrument, banklist in banks.items())
+	assert len(nyquists) == 1, "all banks must have the same Nyquist frequency to define a consistent horizon distance function (got %s)" % ", ".join("%g" % rate for rate in sorted(nyquists))
+	# assume default 32 s PSD.  this is not required to be correct, but
+	# for best accuracy it should not be larger than the true value and
+	# for best performance it should not be smaller than the true
+	# value.
+	deltaF = 1. / 32.
+	# FIXME (from Chad) What is the 5/3 for???
+	# pick (m1, m2) from the median template ranked by Mchirp^(5/3)
+	# to provide the canonical waveform model.  See Maggiore equation
+	# (4.3).
+	assert len(sngl_inspiral_table) > 0, "no templates:  must have templates to define horizon distance function"
+	median_row = sorted(sngl_inspiral_table, key = lambda row: row.mchirp**(5./3.))[len(sngl_inspiral_table) // 2]
+	return reference_psd.HorizonDistance(15.0, 0.85 * max(nyquists), deltaF, median_row.mass1, median_row.mass2)