gstlal-burst/bin/gstlal_cs_triggergen

@@ -116,7 +116,7 @@ class PipelineHandler(simplehandler.Handler):
 		self.update_psd = dict.fromkeys(triggergens, 0)
 		self.reference_psd = reference_psds
 		# create a StreamBurca instance, initialized with the XML document and the coincidence parameters
-		self.streamburca = streamburca.StreamBurca(xmldoc, process.process_id, options.delta_t, min_instruments = 2, verbose = options.verbose)
+		self.streamburca = streamburca.StreamBurca(xmldoc, process.process_id, options.delta_t, BBCoincDef = streamburca.burca.StringCuspBBCoincDef, min_instruments = 2, verbose = options.verbose)
 
 
 	def appsink_new_buffer(self, elem):

gstlal-burst/gst/lal/gstlal_string_triggergen.c

@@ -101,13 +101,13 @@ G_DEFINE_TYPE_WITH_CODE(
 
 static unsigned autocorrelation_channels(const gsl_matrix_float *autocorrelation_matrix)
 {
-	return autocorrelation_matrix->size1;
+	return autocorrelation_matrix ? autocorrelation_matrix->size1 : 0;
 }
 
 
 static unsigned autocorrelation_length(const gsl_matrix_float *autocorrelation_matrix)
 {
-	return autocorrelation_matrix->size2;
+	return autocorrelation_matrix ? autocorrelation_matrix->size2 : 0;
 }
 
 
@@ -198,6 +198,19 @@ static gsl_vector_float *gstlal_autocorrelation_chi2_compute_norms_string(const
  */
 
 
+static GstClockTime buffer_pts(GSTLALStringTriggergen *element, guint64 offset)
+{
+	gint channel;
+	/* start time of SNR used to construct buffer */
+	GstClockTime t = element->t0 + gst_util_uint64_scale_int_round(offset + (autocorrelation_length(element->autocorrelation_matrix) - 1) / 2 - element->offset0, GST_SECOND, GST_AUDIO_INFO_RATE(&element->audio_info));
+	/* pin buffer start time at earliest unrealized trigger */
+	for(channel = 0; channel < element->num_templates; channel++)
+		if(element->bank[channel].snr > 0.0 && (GstClockTime) XLALGPSToINT8NS(&element->bank[channel].peak_time) < t)
+			t = XLALGPSToINT8NS(&element->bank[channel].peak_time);
+	return t;
+}
+
+
 static GstFlowReturn trigger_generator(GSTLALStringTriggergen *element, GstBuffer *outbuf)
 {
 	float *snrdata;
@@ -209,37 +222,35 @@ static GstFlowReturn trigger_generator(GSTLALStringTriggergen *element, GstBuffe
 	guint sample;
 	gint channel;
 
+	g_mutex_lock(&element->bank_lock);
+
+	/* check that autocorrelation vector has odd number of samples.
+	 * NOTE:  autocorrelation_length() returns 0 if the
+	 * autocorrelation_matrix is not set, so this g_assert also tests
+	 * for a missing autocorrelation_matrix.  in set_property(), if the
+	 * conversion from GValueArray fails the matrix will be left set to
+	 * NULL, so this is also catching those failures. */
+	g_assert(autocorrelation_length(element->autocorrelation_matrix) & 1);
+
+	/* do we have enough SNR to search for triggers? */
+	offset = gst_audioadapter_offset(element->adapter);
 	length = get_available_samples(element);
 	if(length < autocorrelation_length(element->autocorrelation_matrix)) {
-		/* FIXME:  PTS and duration are not necessarily correct.
-		 * they're correct for now because we know how this element
-		 * is used in the current pipeline, but in general this
-		 * behaviour is not correct.  right now, the adapter can
-		 * only not have enough data at the start of a stream, but
-		 * for general streams the adapter could get flushed in mid
-		 * stream and then we might need to worry about what the
-		 * last reported buffer's end time was.  maybe.  maybe not
-		 */
-		GST_BUFFER_PTS(outbuf) = element->t0;
+		GST_BUFFER_PTS(outbuf) = buffer_pts(element, offset);
 		GST_BUFFER_DURATION(outbuf) = 0;
 		GST_BUFFER_OFFSET_END(outbuf) = GST_BUFFER_OFFSET(outbuf) + ntriggers;
+		g_mutex_unlock(&element->bank_lock);
 		return GST_FLOW_OK;
 	}
 
-	g_mutex_lock(&element->bank_lock);
-	snrsample = snrdata = g_malloc(length * element->num_templates * sizeof(*snrdata));
-
 	/* copy samples */
-	offset = gst_audioadapter_offset(element->adapter);
+	snrsample = snrdata = g_malloc(length * element->num_templates * sizeof(*snrdata));
 	gst_audioadapter_copy_samples(element->adapter, snrdata, length, NULL, NULL);
 
 	/* compute the chisq norm if it doesn't exist */
 	if(!element->autocorrelation_norm)
 		element->autocorrelation_norm = gstlal_autocorrelation_chi2_compute_norms_string(element->autocorrelation_matrix, NULL);
 
-	/* check that autocorrelation vector has odd number of samples */
-	g_assert(autocorrelation_length(element->autocorrelation_matrix) & 1);
-
 	/* find events.  earliest sample that can be a new trigger starts a
 	 * little bit in from the start of the adapter because we are
 	 * re-using data from the last iteration for \chi^2 calculation.
@@ -323,10 +334,7 @@ static GstFlowReturn trigger_generator(GSTLALStringTriggergen *element, GstBuffe
 	 * obtain PTS and DURATION of output buffer.
 	 */
 
-	GST_BUFFER_PTS(outbuf) = element->t0 + gst_util_uint64_scale_int_round(offset + (autocorrelation_length(element->autocorrelation_matrix) - 1) / 2 - element->offset0, GST_SECOND, GST_AUDIO_INFO_RATE(&element->audio_info));
-	for(channel = 0; channel < element->num_templates; channel++)
-		if(element->bank[channel].snr > 0.0 && (GstClockTime) XLALGPSToINT8NS(&element->bank[channel].peak_time) < GST_BUFFER_PTS(outbuf))
-			GST_BUFFER_PTS(outbuf) = XLALGPSToINT8NS(&element->bank[channel].peak_time);
+	GST_BUFFER_PTS(outbuf) = buffer_pts(element, offset);
 	GST_BUFFER_DURATION(outbuf) = element->t0 + gst_util_uint64_scale_int_round(offset + length - (autocorrelation_length(element->autocorrelation_matrix) - 1) / 2 - element->offset0, GST_SECOND, GST_AUDIO_INFO_RATE(&element->audio_info)) - GST_BUFFER_PTS(outbuf);
 
 	GST_BUFFER_OFFSET_END(outbuf) = GST_BUFFER_OFFSET(outbuf) + ntriggers;

gstlal-burst/python/streamburca.py

@@ -71,15 +71,16 @@ class backgroundcollector(object):
 
 
 class StreamBurca(object):
-	def __init__(self, xmldoc, process_id, delta_t, min_instruments = 2, verbose = False):
+	def __init__(self, xmldoc, process_id, delta_t, BBCoincDef, min_instruments = 2, verbose = False):
 		self.delta_t = delta_t
+		self.BBCoincDef = BBCoincDef
 		self.min_instruments = min_instruments
 		self.verbose = verbose
 		self.set_xmldoc(xmldoc, process_id)
 
 
 	def set_xmldoc(self, xmldoc, process_id):
-		self.coinc_tables = burca.StringCuspCoincTables(xmldoc, burca.StringCuspBBCoincDef)
+		self.coinc_tables = burca.StringCuspCoincTables(xmldoc, self.BBCoincDef)
 		self.sngl_burst_table = lsctables.SnglBurstTable.get_table(xmldoc)
 		self.process_id = process_id
 		self.time_slide_graph = snglcoinc.TimeSlideGraph(
@@ -110,16 +111,9 @@ class StreamBurca(object):
 		newly_reported = []
 		flushed = []
 		flushed_unused = []
-		for node, events in self.time_slide_graph.pull(newly_reported = newly_reported, flushed = flushed, flushed_unused = flushed_unused, coinc_sieve = coinc_sieve, event_collector = self.backgroundcollector, flush = flush, verbose = False):
-			# for exact template match
-			if not burca.StringCuspCoincTables.ntuple_comparefunc(events, node.offset_vector):
-				# construct row objects for coinc tables
-
-				coinc, coincmaps, coinc_burst = self.coinc_tables.coinc_rows(self.process_id, node.time_slide_id, events, u"sngl_burst")
-
-				# finally, append coinc to tables
-
-				self.coinc_tables.append_coinc(coinc, coincmaps, coinc_burst)
+		for node, events in self.time_slide_graph.pull(newly_reported = newly_reported, flushed = flushed, flushed_unused = flushed_unused, coinc_sieve = coinc_sieve, event_collector = self.backgroundcollector, flush = flush):
+			# construct row objects for coinc tables
+			self.coinc_tables.append_coinc(*self.coinc_tables.coinc_rows(self.process_id, node.time_slide_id, events, u"sngl_burst"))
 
 		# add singles into the noise model
 		if flushed and noninjections:

gstlal-calibration/gst/lal/gstlal_transferfunction.c

@@ -169,6 +169,8 @@ enum property {
 	ARG_FFT_WINDOW_TYPE,
 	ARG_FIR_WINDOW_TYPE,
 	ARG_USE_FIR_FFT,
+	ARG_DF,
+	ARG_F0,
 	ARG_FAKE
 };
 
@@ -2218,6 +2220,14 @@ static void set_property(GObject *object, enum property id, const GValue *value,
 		element->use_fir_fft = g_value_get_boolean(value);
 		break;
 
+	case ARG_DF:
+		element->df = g_value_get_double(value);
+		break;
+
+	case ARG_F0:
+		element->f0 = g_value_get_double(value);
+		break;
+
 	default:
 		G_OBJECT_WARN_INVALID_PROPERTY_ID(object, id, pspec);
 		break;
@@ -2385,6 +2395,14 @@ static void get_property(GObject *object, enum property id, GValue *value, GPara
 		g_value_set_boolean(value, element->use_fir_fft);
 		break;
 
+	case ARG_DF:
+		g_value_set_double(value, element->df);
+		break;
+
+	case ARG_F0:
+		g_value_set_double(value, element->f0);
+		break;
+
 	default:
 		G_OBJECT_WARN_INVALID_PROPERTY_ID(object, id, pspec);
 		break;
@@ -2835,6 +2853,23 @@ static void gstlal_transferfunction_class_init(GSTLALTransferFunctionClass *klas
 		FALSE,
 		G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
 	);
+	properties[ARG_DF] = g_param_spec_double(
+		"df",
+		"Frequency spacing",
+		"Frequency spacing of FFT in Hz.  This cannot be less than the default,\n\t\t\t"
+		"which is the inverse of the temporal duration of the FFT length.",
+		0, G_MAXDOUBLE, 0,
+		G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
+	);
+	properties[ARG_F0] = g_param_spec_double(
+		"f0",
+		"Starting frequency",
+		"Set this to shift the frequencies to be measured.  If unset, measurement\n\t\t\t"
+		"frequencies are multiples of df.  If set, measurement frequencies are\n\t\t\t"
+		"shifted by this amount.",
+		0, G_MAXDOUBLE, 0,
+		G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
+	);
 
 
 	g_object_class_install_property(
@@ -2962,6 +2997,16 @@ static void gstlal_transferfunction_class_init(GSTLALTransferFunctionClass *klas
 		ARG_USE_FIR_FFT,
 		properties[ARG_USE_FIR_FFT]
 	);
+	g_object_class_install_property(
+		gobject_class,
+		ARG_DF,
+		properties[ARG_DF]
+	);
+	g_object_class_install_property(
+		gobject_class,
+		ARG_F0,
+		properties[ARG_F0]
+	);
 }
 
 

gstlal-calibration/gst/lal/gstlal_transferfunction.h

@@ -200,6 +200,8 @@ struct _GSTLALTransferFunction {
 	enum gstlal_transferfunction_window_type fft_window_type;
 	enum gstlal_transferfunction_window_type fir_window_type;
 	gboolean use_fir_fft;
+	double df;
+	double f0;
 };
 
 

gstlal-calibration/lib/gstlal_firtools.c

@@ -2108,7 +2108,7 @@ LONG DTYPE *dpss_ ## LONG ## DTYPE(guint N, double alpha, double compute_time, L
 	static GMutex mutex; \
 	guint i, j, start; \
 	static gpointer *warehouse = NULL; \
-	static guint * warehouse_N = NULL; \
+	static guint *warehouse_N = NULL; \
 	static double *warehouse_alpha = NULL; \
 	static double *warehouse_compute_time = NULL; \
 	gboolean add_to_warehouse = FALSE; \

gstlal-calibration/tests/line_subtraction_ringing.py

@@ -26,6 +26,7 @@
 
 import matplotlib; matplotlib.use('Agg')
 import numpy as np
+import os
 import matplotlib.patches as mpatches
 from matplotlib import rc
 rc('text', usetex = True)
@@ -58,14 +59,17 @@ from ticks_and_grid import ticks_and_grid
 
 
 # Constants
-freq = 410.3		# Hz
+freq = 10.3		# Hz
 test_duration = 400000	# seconds
-amplitude_mod = 0.001
-phase_mod = 0.001
-avg_time = 128		# seconds
-other_lines = [410.301, 410.295, 410.31, 410.28]	# Hz
-other_mod = 0.001
-
+amplitude_mod = 0.000
+phase_mod = 0.000
+avg_time = 512		# seconds
+other_lines = [10.301, 10.295, 10.31, 10.28]	# Hz
+other_mod = 0.000
+fft_length = 16000 * 64
+fft_overlap = 8000 * 64
+num_ffts = (test_duration - 1000 - 8000) // (16000 - 8000)
+tf_length = fft_length // 2 + 1
 
 #
 # =============================================================================
@@ -84,11 +88,11 @@ def line_subtraction_ringing_01(pipeline, name):
 	#
 
 	# The noise floor
-	noise = test_common.test_src(pipeline, rate = 2048, test_duration = test_duration, wave = 5, src_suffix = '0')
+	noise = test_common.test_src(pipeline, rate = 64, test_duration = test_duration, wave = 5, src_suffix = '0')
 	noise = pipeparts.mktee(pipeline, noise)
 
 	# The witness channel used to subtract the line, in this case a pure sinusoid.
-	witness = test_common.test_src(pipeline, rate = 2048, test_duration = test_duration, wave = 0, freq = freq, src_suffix = '1')
+	witness = test_common.test_src(pipeline, rate = 64, test_duration = test_duration, wave = 0, freq = freq, src_suffix = '1')
 
 	# The signal in the data, which may not be a pure sinusoid, but could have amplitude and/or phase modulation.
 	amplitude = test_common.test_src(pipeline, rate = 16, test_duration = test_duration, wave = 0, src_suffix = '2', volume = 0)
@@ -104,11 +108,11 @@ def line_subtraction_ringing_01(pipeline, name):
 	p_mod = pipeparts.mkgeneric(pipeline, p_mod, "lal_smoothkappas", array_size = 1, avg_array_size = 16 * 2500, default_kappa_re = 0)
 	p_mod = pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, p_mod, matrix = [[0.0, phase_mod * 4 * 50]]))
 	signal = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, amplitude, pipeparts.mkgeneric(pipeline, p_mod, "cexp")))
-	signal = calibration_parts.mkresample(pipeline, signal, 4, False, "audio/x-raw, format=Z128LE, rate=2048")
+	signal = calibration_parts.mkresample(pipeline, signal, 4, False, "audio/x-raw, format=Z128LE, rate=64")
 	signal = pipeparts.mkgeneric(pipeline, signal, "lal_demodulate", line_frequency = -1.0 * freq, prefactor_real = -2.0)
-	signal = pipeparts.mkcapsfilter(pipeline, signal, "audio/x-raw, format=Z128LE, rate=2048")
+	signal = pipeparts.mkcapsfilter(pipeline, signal, "audio/x-raw, format=Z128LE, rate=64")
 	signal = pipeparts.mkgeneric(pipeline, signal, "creal")
-	signal = pipeparts.mkcapsfilter(pipeline, signal, "audio/x-raw, format=F64LE, rate=2048")
+	signal = pipeparts.mkcapsfilter(pipeline, signal, "audio/x-raw, format=F64LE, rate=64")
 
 	# Add noise to signal
 	strain = pipeparts.mktee(pipeline, calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, signal, noise)))
@@ -129,28 +133,32 @@ def line_subtraction_ringing_01(pipeline, name):
 		p_mod = pipeparts.mkgeneric(pipeline, p_mod, "lal_smoothkappas", array_size = 1, avg_array_size = 16 * 2500, default_kappa_re = 0)
 		p_mod = pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, p_mod, matrix = [[0.0, other_mod * 4 * 50]]))
 		signal = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, amplitude, pipeparts.mkgeneric(pipeline, p_mod, "cexp")))
-		signal = calibration_parts.mkresample(pipeline, signal, 4, False, "audio/x-raw, format=Z128LE, rate=2048")
+		signal = calibration_parts.mkresample(pipeline, signal, 4, False, "audio/x-raw, format=Z128LE, rate=64")
 		signal = pipeparts.mkgeneric(pipeline, signal, "lal_demodulate", line_frequency = -1.0 * line_freq, prefactor_real = -2.0)
-		signal = pipeparts.mkcapsfilter(pipeline, signal, "audio/x-raw, format=Z128LE, rate=2048")
+		signal = pipeparts.mkcapsfilter(pipeline, signal, "audio/x-raw, format=Z128LE, rate=64")
 		signal = pipeparts.mkgeneric(pipeline, signal, "creal")
-		signal = pipeparts.mkcapsfilter(pipeline, signal, "audio/x-raw, format=F64LE, rate=2048")
+		signal = pipeparts.mkcapsfilter(pipeline, signal, "audio/x-raw, format=F64LE, rate=64")
 
 		# Add signal to strain
 		strain = pipeparts.mktee(pipeline, calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, strain, signal)))
 
 	# Subtract signal using same function that is used in gstlal_compute_strain
-	clean_strain = calibration_parts.remove_lines_with_witnesses(pipeline, strain, [[witness]], [[freq]], [0], [], filter_latency = 1, compute_rate = 16, rate_out = 2048, num_median = 1 * 16, num_avg = avg_time * 16, noisesub_gate_bit = None)
+	clean_strain = calibration_parts.remove_lines_with_witnesses(pipeline, strain, [[witness]], [[freq]], [0], [], filter_latency = 1, compute_rate = 16, rate_out = 64, num_median = 1 * 16, num_avg = avg_time * 16, noisesub_gate_bit = None)
 
-	#strain = calibration_parts.mkresample(pipeline, strain, 4, False, "audio/x-raw, format=F64LE, rate=2048")
-	#clean_strain = calibration_parts.mkresample(pipeline, clean_strain, 4, False, "audio/x-raw, format=F64LE, rate=2048")
+	#strain = calibration_parts.mkresample(pipeline, strain, 4, False, "audio/x-raw, format=F64LE, rate=64")
+	#clean_strain = calibration_parts.mkresample(pipeline, clean_strain, 4, False, "audio/x-raw, format=F64LE, rate=64")
 	# Remove the initial data
 	strain = calibration_parts.mkinsertgap(pipeline, strain, insert_gap = False, chop_length = 600 * 1000000000)
 	clean_strain = calibration_parts.mkinsertgap(pipeline, clean_strain, insert_gap = False, chop_length = 600 * 1000000000)
-	strain = pipeparts.mkprogressreport(pipeline, strain, "strain")
-	clean_strain = pipeparts.mkprogressreport(pipeline, clean_strain, "clean_strain")
+	strain = pipeparts.mktee(pipeline, pipeparts.mkprogressreport(pipeline, strain, "strain"))
+	clean_strain = pipeparts.mktee(pipeline, pipeparts.mkprogressreport(pipeline, clean_strain, "clean_strain"))
 	# Compute ASDs
-	pipeparts.mkgeneric(pipeline, strain, "lal_asd", fft_samples = 16000 * 2048, overlap_samples = 2000 * 2048, window_type = 3, filename = "strainASD_%davg.txt" % avg_time)
-	pipeparts.mkgeneric(pipeline, clean_strain, "lal_asd", fft_samples = 16000 * 2048, overlap_samples = 2000 * 2048, window_type = 3, filename = "clean_strainASD_%davg.txt" % avg_time)
+	pipeparts.mkgeneric(pipeline, strain, "lal_asd", fft_samples = 16000 * 64, overlap_samples = 8000 * 64, window_type = 3, filename = "strainASD_%davg.txt" % avg_time)
+	pipeparts.mkgeneric(pipeline, clean_strain, "lal_asd", fft_samples = 16000 * 64, overlap_samples = 8000 * 64, window_type = 3, filename = "clean_strainASD_%davg.txt" % avg_time)
+
+	# Compute transfer function
+	interleaved = calibration_parts.mkinterleave(pipeline, [clean_strain, strain])
+	pipeparts.mkgeneric(pipeline, interleaved, "lal_transferfunction", fft_length = fft_length, fft_overlap = fft_overlap, num_ffts = num_ffts, use_median = False, update_samples = 1e15, fft_window_type = 3, filename = "lineSubRingingTF_%davg.txt" % avg_time)
 
 	#
 	# done
@@ -180,7 +188,31 @@ clean = np.loadtxt("clean_strainASD_%davg.txt" % avg_time)
 print("Done loading ASD txt file 2 of 2")
 clean = np.transpose(clean)[1]
 
-# Plot against frequency
+# Get TF data
+print("Loading TF txt file")
+# Remove unwanted lines from TF file, and re-format wanted lines
+f = open("lineSubRingingTF_%davg.txt" % avg_time, "r")
+lines = f.readlines()
+f.close()
+os.system("rm lineSubRingingTF_%davg.txt" % avg_time)
+f = open("lineSubRingingTF_%davg.txt" % avg_time, "w")
+for j in range(3, 3 + tf_length):
+	f.write(lines[j].replace(' + ', '\t').replace(' - ', '\t-').replace('i', ''))
+f.close()
+
+TF = np.loadtxt("lineSubRingingTF_%davg.txt" % avg_time)
+print("Done loading TF txt file")
+
+TFfvec = []
+TFmag = []
+TFphase = []
+for j in range(0, len(TF)):
+	TFfvec.append(TF[j][0])
+	tf_at_f = (TF[j][1] + 1j * TF[j][2])
+	TFmag.append(abs(tf_at_f))
+	TFphase.append(np.angle(tf_at_f) * 180.0 / np.pi)
+
+# Plot ASDs against frequency
 plt.figure(figsize = (10, 6))
 colors = ['red', 'blue']
 labels = ['strain', 'clean']
@@ -188,14 +220,32 @@ plt.plot(fvec, strain, colors[0], linewidth = 0.75)
 plt.plot(fvec, clean, colors[1], linewidth = 0.75)
 patches = [mpatches.Patch(color = colors[j], label = r'$%s$' % labels[j]) for j in range(len(labels))]
 plt.legend(handles = patches, loc = 'upper right', ncol = 1)
-plt.title('%d s Average (Modulated Lines)' % (avg_time))
+plt.title('ASD with %d s Average' % (avg_time))
 plt.ylabel(r'${\rm ASD}\ \left[{\rm strain / }\sqrt{\rm Hz}\right]$')
 plt.xlabel(r'${\rm Frequency \ [Hz]}$')
 ticks_and_grid(plt.gca(), xmin = freq - 0.03, xmax = freq + 0.02, ymin = 0.0001, ymax = 10000, xscale = 'linear', yscale = 'log')
 
 if any(other_lines):
-	plt.savefig("lineSubRingingOtherModLinesASD_%davg.png" % (avg_time))
+	plt.savefig("lineSubRingingOtherLinesASD_%davg.png" % (avg_time))
 else:
 	plt.savefig("lineSubRingingASD_%davg.png" % (avg_time))
 
+# Plot transfer function against frequency
+plt.figure(figsize = (10, 10))
+plt.subplot(211)
+plt.plot(TFfvec, TFmag, 'green', linewidth = 0.75)
+ticks_and_grid(plt.gca(), xmin = freq - 0.03, xmax = freq + 0.02, ymin = 0, ymax = 2, xscale = 'linear', yscale = 'linear')
+plt.ylabel('Magnitude')
+plt.title('TF with %d s Average' % (avg_time))
+plt.subplot(212)
+plt.plot(TFfvec, TFphase, 'green', linewidth = 0.75)
+ticks_and_grid(plt.gca(), xmin = freq - 0.03, xmax = freq + 0.02, ymin = -180, ymax = 180, xscale = 'linear', yscale = 'linear')
+plt.ylabel('Phase [deg]')
+plt.xlabel('Frequency [Hz]')
+
+if any(other_lines):
+	plt.savefig("lineSubRingingOtherLinesTF_%davg.png" % (avg_time))
+else:
+	plt.savefig("lineSubRingingTF_%davg.png" % (avg_time))
+
 

gstlal-inspiral/bin/gstlal_inspiral

@@ -373,8 +373,10 @@ def parse_command_line():
 		raise ValueError("must supply exactly as many --svd-bank options as --output")
 	if options.ranking_stat_output and len(options.ranking_stat_output) != len(options.output):
 		raise ValueError("must supply either none or exactly as many --ranking-stat-output options as --output")
-	if options.likelihood_snapshot_interval and not options.ranking_stat_output:
+	if (options.likelihood_snapshot_interval and not options.ranking_stat_output) and not options.injections:
 		raise ValueError("must set --ranking-stat-output when --likelihood-snapshot-interval is set")
+	if options.ranking_stat_output and options.injections:
+		raise ValueError("must not set --ranking-stat-output when --injections is set")
 	if options.ranking_stat_output is None or len(options.ranking_stat_output) == 0:
 		options.ranking_stat_output = [None] * len(options.output)
 
@@ -403,9 +405,12 @@ def parse_command_line():
 
 	if options.data_source in ("lvshm", "framexmit"):
 		missing_options = []
-		for option in ["job_tag", "ranking_stat_input", "ranking_stat_pdf", "zerolag_rankingstat_pdf"]:
+		for option in ["job_tag", "ranking_stat_input", "ranking_stat_pdf"]:
 			if getattr(options, option) is None:
 				missing_options.append("--%s" %option.replace("_","-"))
+		if not options.injections:
+			if options.zerolag_rankingstat_pdf is None:
+				missing_options.append("--zerolag-rankingstat-pdf")
 
 		if missing_options:
 			raise ValueError("missing required option(s) %s when --data-source is lvshm or framexmit" % ", ".join(missing_options))
@@ -776,7 +781,7 @@ for output_file_number, (svd_bank_url_dict, output_url, ranking_stat_output_url,
 	# the new injection file to disk before the pipeline starts up
 	#
 
-	if options.injections and options.reference_psd:
+	if (options.injections and options.reference_psd) and not options.data_source in ("lvshm", "framexmit"):
 		if options.verbose:
 			print("calculating expected SNR for injections...", file=sys.stderr)
 

gstlal-inspiral/bin/gstlal_ll_inspiral_event_plotter

@@ -168,7 +168,7 @@ class EventPlotter(events.EventProcessor):
 		payload = json.loads(message.value())
 
 		time = LIGOTimeGPS(payload['time'], payload['time_ns'])
-		coinc_fileobj = io.StringIO(payload['coinc'])
+		coinc_fileobj = io.BytesIO(payload['coinc'].encode('utf-8'))
 		xmldoc = ligolw_utils.load_fileobj(coinc_fileobj, contenthandler=ligolwcontenthandler)
 		coinc_fileobj.close()
 		sngl_inspiral_table = lsctables.SnglInspiralTable.get_table(xmldoc)
@@ -193,7 +193,7 @@ class EventPlotter(events.EventProcessor):
 		elif message.topic() == self.upload_topic:
 			self.events[event_key]['gid'] = payload['gid']
 			self.events[event_key]['coinc'] = xmldoc
-			psd_fileobj = io.StringIO(payload['psd'])
+			psd_fileobj = io.BytesIO(payload['psd'].encode('utf-8'))
 			self.events[event_key]['psd'] = ligolw_utils.load_fileobj(psd_fileobj, contenthandler=series.PSDContentHandler)
 			psd_fileobj.close()
 

gstlal-inspiral/bin/gstlal_ll_inspiral_event_uploader

@@ -32,7 +32,7 @@ import json
 import logging
 import math
 from optparse import OptionParser
-from io import StringIO
+import io
 import sys
 import time
 
@@ -363,9 +363,9 @@ class EventUploader(events.EventProcessor):
 			return None
 
 		# gzip and upload psd
-		psd_fobj = StringIO.StringIO()
+		psd_fobj = io.BytesIO()
 		with gzip.GzipFile(fileobj=psd_fobj, mode="wb", compresslevel=3) as f:
-			f.write(event['favored']['psd'])
+			f.write(event['favored']['psd'].encode("utf-8"))
 		self.upload_file("strain PSDs", "psd.xml.gz", "psd", psd_fobj.getvalue(), graceid)
 		del psd_fobj
 

gstlal-inspiral/bin/gstlal_ll_inspiral_workflow

@@ -60,9 +60,9 @@ else:
 	marg_pdf = DataCache.generate(DataType.DIST_STAT_PDFS, config.all_ifos)
 
 	# generate dag layers
-	dag.filter_online(svd_banks, dist_stats, zerolag_pdfs, marg_pdf)
 	if config.filter.injections:
 		dag.filter_injections_online(svd_banks, dist_stats, zerolag_pdfs, marg_pdf)
+	dag.filter_online(svd_banks, dist_stats, zerolag_pdfs, marg_pdf)
 	dag.marginalize_online(marg_pdf)
 	dag.track_noise()
 	if config.services.kafka_server:

gstlal-inspiral/python/dags/layers/inspiral.py

@@ -1323,7 +1323,7 @@ def filter_online_layer(config, dag, svd_bank_cache, dist_stat_cache, zerolag_pd
 	if config.upload.before_merger:
 		common_opts.append(Option("upload-time-before-merger"))
 
-	if config.filter.delay_uploads:
+	if config.upload.delay_uploads:
 		common_opts.append(Option("delay-uploads"))
 
 	if config.filter.cap_singles:
@@ -1430,7 +1430,7 @@ def filter_injections_online_layer(config, dag, svd_bank_cache, dist_stat_cache,
 	if config.upload.before_merger:
 		common_opts.append(Option("upload-time-before-merger"))
 
-	if config.filter.delay_uploads:
+	if config.upload.delay_uploads:
 		common_opts.append(Option("delay-uploads"))
 
 	if config.filter.cap_singles:

gstlal-inspiral/python/inspiral.py

@@ -63,6 +63,7 @@ import time
 import http.client
 import tempfile
 import os
+import bisect
 from urllib.parse import urlparse
 
 from ligo.lw import ligolw
@@ -406,7 +407,13 @@ class CoincsDocument(object):
 		#
 
 		if injection_filename is not None:
-			ligolw_add.ligolw_add(self.xmldoc, [injection_filename], contenthandler = LIGOLWContentHandler, verbose = verbose)
+			# sort and add the sim_insprial
+			xmldoc_sim_inspiral = ligolw_utils.load_filename(injection_filename, verbose = verbose, contenthandler = LIGOLWContentHandler)
+			self.sim_inspiral_table = lsctables.SimInspiralTable.get_table(xmldoc_sim_inspiral)
+			self.sim_inspiral_table.sort(key=lambda s:s.geocent_end_time+1.e-9*s.geocent_end_time_ns)
+			self.xmldoc.childNodes[-1].appendChild(self.sim_inspiral_table)
+		else:
+			self.sim_inspiral_table = None
 
 		#
 		# insert time slide offset vectors.  remove duplicate
@@ -650,12 +657,23 @@ class GracedBWrapper(object):
 		# check for no-op
 		if not gracedb_ids:
 			return
-		fobj = io.StringIO()
+		fobj = io.BytesIO()
 		ligolw_utils.write_fileobj(xmldoc, fobj)
 		self.__upload_aux_data(message, filename, tag, fobj.getvalue(), gracedb_ids)
 		del fobj
 
-	def do_alerts(self, last_coincs, psddict, rankingstat_xmldoc_func, seglistdicts, get_p_astro_func):
+	def nearest_sim_table(self, gps_time, sim_inspiral):
+		# search through the sim inspiral table
+		# find the row with geo-centric ends time nearest to the specific gps time (using by bisect)
+		# return the new sim inspiral table with a single row corresponding to the nearest time 
+		sim_inspiral_table = lsctables.New(lsctables.SimInspiralTable)
+		end_times = [inj.time_geocent for inj in sim_inspiral]
+		idx = bisect.bisect_left(end_times, gps_time)
+		row = sim_inspiral[idx]
+		sim_inspiral_table.append(row)
+		return sim_inspiral_table
+
+	def do_alerts(self, last_coincs, psddict, rankingstat_xmldoc_func, seglistdicts, get_p_astro_func, sim_inspiral_table = None):
 		gracedb_ids = []
 
 		# no-op short circuit
@@ -724,6 +742,10 @@ class GracedBWrapper(object):
 			print("sending %s to gracedb ..." % filename, file=sys.stderr)
 			message = io.BytesIO()
 			xmldoc = last_coincs[coinc_event.coinc_event_id]
+
+			if sim_inspiral_table:
+				xmldoc.childNodes[-1].appendChild(self.nearest_sim_table(end_time, sim_inspiral_table))
+
 			# give the alert all the standard inspiral
 			# columns (attributes should all be
 			# populated).  FIXME:  ugly.
@@ -990,7 +1012,7 @@ class GracedBWrapper(object):
 								print(resp, file=sys.stderr)
 							break
 					print("gracedb upload of %s failed on attempt %d/%d" % (filename, attempt, self.retries), file=sys.stderr)
-					print(resp_json, file=sys.stderr)
+					print(resp, file=sys.stderr)
 					time.sleep(random.lognormal(math.log(self.retry_delay), .5))
 				else:
 					print("gracedb upload of %s failed" % filename, file=sys.stderr)

gstlal-inspiral/python/inspiral_pipe.py

@@ -66,7 +66,7 @@ from gstlal.datafind import DataType
 
 
 warnings.warn(
-	"all functionality within this module has been captured in gstlal.segutil and"
+	"all functionality within this module has been captured in gstlal.segments and"
 	"gstlal.dags.layers.inspiral, and will be removed from gstlal in the future",
 	DeprecationWarning,
 )

gstlal-inspiral/python/lloidhandler.py

@@ -1498,7 +1498,7 @@ class Handler(simplehandler.Handler):
 		assert self.fapfar is not None
 
 		# do alerts
-		self.gracedbwrapper.do_alerts(last_coincs, self.psds, self.__get_rankingstat_xmldoc_for_gracedb, self.segmentstracker.seglistdicts, self.__get_p_astro_json)
+		self.gracedbwrapper.do_alerts(last_coincs, self.psds, self.__get_rankingstat_xmldoc_for_gracedb, self.segmentstracker.seglistdicts, self.__get_p_astro_json, sim_inspiral_table = self.coincs_document.sim_inspiral_table)
 
 
 	def web_get_sngls_snr_threshold(self):

gstlal-inspiral/python/stats/inspiral_intrinsics.py

@@ -112,14 +112,14 @@ class SourcePopulationModel(object):
 		"""
 		if filename is not None:
 			with h5py.File(filename, 'r') as model:
-				coefficients = model['coefficients'].value
-				snr_bp = model['SNR'].value
+				coefficients = model['coefficients'][()]
+				snr_bp = model['SNR'][()]
 				try:
-					model_ids = model['template_id'].value
+					model_ids = model['template_id'][()]
 				except KeyError:
                                         # FIXME: assume sequential order if model['event_id'] doesn't exist
-                                        #model_ids = numpy.arange(numpy.shape(model['coefficients'].value)[-1])
-					model_ids = model['event_id'].value
+                                        #model_ids = numpy.arange(numpy.shape(model['coefficients'][()])[-1])
+					model_ids = model['event_id'][()]
 			# PPoly can construct an array of polynomials by just
 			# feeding it the coefficients array all in one go, but then
 			# it insists on evaluating all of them at once.  we don't

gstlal-inspiral/tests/dtdphitest.py

@@ -81,10 +81,10 @@ i = 0
 # the following cov matrix elements were derived from running `gstlal_inspiral_compute_dtdphideff_cov_matrix --psd-xml share/O3/2019-05-09-H1L1V1psd_new.xml.gz --H-snr 5 --L-snr 7.0 --V-snr 4.0 --K-snr 4.0`
 data = h5py.File(pdf_fname)
 # pdb.set_trace()
-refsnr = dict((ifo, data["SNR"][ifo].value) for ifo in data["SNR"])
-transmat_dtdphi = numpy.array([[data["transtt"][combo].value, data["transtp"][combo].value], [data["transpt"][combo].value, data["transpp"][combo].value]])
+refsnr = dict((ifo, data["SNR"][ifo][()]) for ifo in data["SNR"])
+transmat_dtdphi = numpy.array([[data["transtt"][combo][()], data["transtp"][combo][()]], [data["transpt"][combo][()], data["transpp"][combo][()]]])
 covmat_dtdphi = numpy.linalg.inv(numpy.dot(transmat_dtdphi.T, transmat_dtdphi))
-sigmadd = 1. / (data["transdd"][combo].value)**2
+sigmadd = 1. / (data["transdd"][combo][()])**2
 rd_slice = (refhorizon[ifo1] / refsnr[ifo1]) / (refhorizon[ifo2] / refsnr[ifo2])
 while i < ndraws:
 	t,d,p = random_source()

gstlal-ugly/bin/Makefile.am

@@ -22,7 +22,6 @@ dist_bin_SCRIPTS = \
 	gstlal_inspiral_treebank \
 	gstlal_inspiral_treebank_dag \
 	gstlal_inspiral_bankviz \
-	gstlal_query_gwosc_segments \
 	gstlal_segments_trim \
 	gstlal_glitch_population \
 	gstlal_vetoes_from_burst_inj \

gstlal-ugly/python/events.py

@@ -68,7 +68,11 @@ class EventProcessor(object):
 			'bootstrap.servers': kafka_server,
 			'group.id': '-'.join([self._name, tag])
 		}
-		self.producer = Producer(self.kafka_settings)
+		self.producer_settings = {
+			'message.max.bytes': 5242880, # 5 MB
+			**self.kafka_settings,
+		}
+		self.producer = Producer(self.producer_settings)
 		if not self.is_source:
 			self.consumer = Consumer(self.kafka_settings)
 			self.consumer.subscribe([topic for topic in input_topic])

gstlal/bin/Makefile.am

@@ -12,6 +12,10 @@ dist_bin_SCRIPTS = \
 	gstlal_psd_polyfit \
 	gstlal_psd_workflow \
 	gstlal_psd_xml_from_asd_txt \
+	gstlal_query_dqsegdb_segments \
+	gstlal_query_dqsegdb_veto_segments \
+	gstlal_query_gwosc_segments \
+	gstlal_query_gwosc_veto_segments \
 	gstlal_reference_psd \
 	gstlal_spectrum_movie \
 	gstlal_stage_frames \