diff --git a/gstlal-burst/lib/gstlal_sngltrigger.c b/gstlal-burst/lib/gstlal_sngltrigger.c index 9cb1dc32f68dc93ff6ab11fd145bfa7cad6f0f0d..30547d30471b956f706e2cff7ca9917a4bc1a283 100644 --- a/gstlal-burst/lib/gstlal_sngltrigger.c +++ b/gstlal-burst/lib/gstlal_sngltrigger.c @@ -1,5 +1,6 @@ /* * Copyright (C) 2012,2013,2015,2016 Chad Hanna + * 2017-2018 Duncan Meacher * * 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 @@ -57,120 +58,129 @@ GstBuffer *gstlal_sngltrigger_new_buffer_from_peak(struct gstlal_peak_state *inp GST_BUFFER_DURATION(srcbuf) = (GstClockTime) gst_util_uint64_scale_int_round(GST_SECOND, length, rate); if (input->num_events) { - if(max_snr) - { - max_snr_index = gstlal_peak_max_over_channels(input); - /* Type casting unsigned int (guint) to int */ - g_assert(max_snr_index >= 0 && max_snr_index < (int)input->channels); - } - guint channel; - for(channel = 0; channel < input->channels; channel++) { - struct GSTLALSnglTrigger *event; - SnglTriggerTable *parent; - double complex maxdata_channel = 0; + /* create a single buffer from max snr channel */ + if(max_snr) { + max_snr_index = gstlal_peak_max_over_channels(input); /* Type casting unsigned int (guint) to int */ - if(max_snr && (int)channel != max_snr_index) - continue; + g_assert(max_snr_index >= 0 && max_snr_index < (int)input->channels); - switch (input->type) - { - case GSTLAL_PEAK_COMPLEX: - maxdata_channel = (double complex) input->interpvalues.as_float_complex[channel]; - break; + /* create a single buffer from max snr channel */ + add_buffer_from_channel(input, channel_name, pad, length, time, rate, chi2, snr_matrix_view, max_snr_index, srcbuf); - case GSTLAL_PEAK_DOUBLE_COMPLEX: - maxdata_channel = (double complex) input->interpvalues.as_double_complex[channel]; - break; + /* add a buffer from each channel */ + } else { - default: - g_assert(input->type == GSTLAL_PEAK_COMPLEX || input->type == GSTLAL_PEAK_DOUBLE_COMPLEX); + guint channel; + for(channel = 0; channel < input->channels; channel++) { + add_buffer_from_channel(input, channel_name, pad, length, time, rate, chi2, snr_matrix_view, channel, srcbuf); } + } + } - if (!maxdata_channel) - continue; - - /* - * allocate new event structure - */ - - - /* - * Populate the SNR snippet if available - * FIXME: only supported for single precision at the moment - */ - if (snr_matrix_view) - { - /* Get the column of SNR we are interested in */ - gsl_vector_complex_float_view snr_vector_view = gsl_matrix_complex_float_column(&(snr_matrix_view->matrix), channel); - /* Allocate an empty time series to hold it. The event takes ownership, so no need to free it*/ - event = gstlal_sngltrigger_new(snr_vector_view.vector.size); - /* Make a GSL view of the time series array data */ - gsl_vector_complex_float_view snr_series_view = gsl_vector_complex_float_view_array((float *) event->snr, event->length); - /* Use BLAS to do the copy */ - gsl_blas_ccopy (&(snr_vector_view.vector), &(snr_series_view.vector)); - } else - event = gstlal_sngltrigger_new(0); - - parent = (SnglTriggerTable *) event; - if (!event) { - /* FIXME handle error */ - } + return srcbuf; +} - /* - * populate - */ - - //*parent = bankarray[channel]; - strcpy(parent->channel, channel_name); - parent->snr = cabs(maxdata_channel); - parent->phase = carg(maxdata_channel); - parent->channel_index = channel; - - XLALINT8NSToGPS(&event->epoch, time); - { - LIGOTimeGPS end_time = event->epoch; - XLALGPSAdd(&end_time, (double) input->interpsamples[channel] / rate); - XLALGPSAddGPS(&parent->end, &end_time); - } - XLALGPSAdd(&event->epoch, (double) (input->samples[channel] - input->pad) / rate); - event->deltaT = 1. / rate; - - /* populate chi squared if we have it */ - parent->chisq = 0.0; - switch (input->type) - { - case GSTLAL_PEAK_COMPLEX: - if (chi2) parent->chisq = (double) *(((float *) chi2 ) + channel); - break; - - case GSTLAL_PEAK_DOUBLE_COMPLEX: - if (chi2) parent->chisq = (double) *(((double *) chi2 ) + channel); - break; - - default: - g_assert(input->type == GSTLAL_PEAK_COMPLEX || input->type == GSTLAL_PEAK_DOUBLE_COMPLEX); - } +void add_buffer_from_channel(struct gstlal_peak_state *input, char *channel_name, GstPad *pad, guint64 length, GstClockTime time, guint rate, void *chi2, gsl_matrix_complex_float_view *snr_matrix_view, int channel, GstBuffer *srcbuf) +{ + struct GSTLALSnglTrigger *event; + SnglTriggerTable *parent; + double complex maxdata_channel = 0; + + switch (input->type) + { + case GSTLAL_PEAK_COMPLEX: + maxdata_channel = (double complex) input->interpvalues.as_float_complex[channel]; + break; + + case GSTLAL_PEAK_DOUBLE_COMPLEX: + maxdata_channel = (double complex) input->interpvalues.as_double_complex[channel]; + break; + + default: + g_assert(input->type == GSTLAL_PEAK_COMPLEX || input->type == GSTLAL_PEAK_DOUBLE_COMPLEX); + } - /* - * add to buffer - */ - - gst_buffer_append_memory( - srcbuf, - gst_memory_new_wrapped( - GST_MEMORY_FLAG_READONLY | GST_MEMORY_FLAG_PHYSICALLY_CONTIGUOUS, - event, - sizeof(*event) + sizeof(event->snr[0]) * event->length, - 0, - sizeof(*event) + sizeof(event->snr[0]) * event->length, - event, - (GDestroyNotify) gstlal_sngltrigger_free - ) - ); - } + if (!maxdata_channel) + return; + + /* + * allocate new event structure + */ + + + /* + * Populate the SNR snippet if available + * FIXME: only supported for single precision at the moment + */ + if (snr_matrix_view) + { + /* Get the column of SNR we are interested in */ + gsl_vector_complex_float_view snr_vector_view = gsl_matrix_complex_float_column(&(snr_matrix_view->matrix), channel); + /* Allocate an empty time series to hold it. The event takes ownership, so no need to free it*/ + event = gstlal_sngltrigger_new(snr_vector_view.vector.size); + /* Make a GSL view of the time series array data */ + gsl_vector_complex_float_view snr_series_view = gsl_vector_complex_float_view_array((float *) event->snr, event->length); + /* Use BLAS to do the copy */ + gsl_blas_ccopy (&(snr_vector_view.vector), &(snr_series_view.vector)); + } else + event = gstlal_sngltrigger_new(0); + + parent = (SnglTriggerTable *) event; + if (!event) { + /* FIXME handle error */ } - return srcbuf; + /* + * populate + */ + + //*parent = bankarray[channel]; + strcpy(parent->channel, channel_name); + parent->snr = cabs(maxdata_channel); + parent->phase = carg(maxdata_channel); + parent->channel_index = channel; + + XLALINT8NSToGPS(&event->epoch, time); + { + LIGOTimeGPS end_time = event->epoch; + XLALGPSAdd(&end_time, (double) input->interpsamples[channel] / rate); + XLALGPSAddGPS(&parent->end, &end_time); + } + XLALGPSAdd(&event->epoch, (double) (input->samples[channel] - input->pad) / rate); + event->deltaT = 1. / rate; + + /* populate chi squared if we have it */ + parent->chisq = 0.0; + switch (input->type) + { + case GSTLAL_PEAK_COMPLEX: + if (chi2) parent->chisq = (double) *(((float *) chi2 ) + channel); + break; + + case GSTLAL_PEAK_DOUBLE_COMPLEX: + if (chi2) parent->chisq = (double) *(((double *) chi2 ) + channel); + break; + + default: + g_assert(input->type == GSTLAL_PEAK_COMPLEX || input->type == GSTLAL_PEAK_DOUBLE_COMPLEX); + } + + /* + * add to buffer + */ + + gst_buffer_append_memory( + srcbuf, + gst_memory_new_wrapped( + GST_MEMORY_FLAG_READONLY | GST_MEMORY_FLAG_PHYSICALLY_CONTIGUOUS, + event, + sizeof(*event) + sizeof(event->snr[0]) * event->length, + 0, + sizeof(*event) + sizeof(event->snr[0]) * event->length, + event, + (GDestroyNotify) gstlal_sngltrigger_free + ) + ); } diff --git a/gstlal-burst/lib/gstlal_sngltrigger.h b/gstlal-burst/lib/gstlal_sngltrigger.h index f3ba91e256699c39f4af93aa2a50a1ba75c1d527..04a4cac5b18458a86e0c8a2b483b1cdeea3250c1 100644 --- a/gstlal-burst/lib/gstlal_sngltrigger.h +++ b/gstlal-burst/lib/gstlal_sngltrigger.h @@ -65,6 +65,7 @@ SnglTriggerTable; */ GstBuffer *gstlal_sngltrigger_new_buffer_from_peak(struct gstlal_peak_state *input, char *channel_name, GstPad *pad, guint64 offset, guint64 length, GstClockTime time, guint rate, void *chi2, gsl_matrix_complex_float_view *snr_matrix_view, GstClockTimeDiff, gboolean max_snr); +void add_buffer_from_channel(struct gstlal_peak_state *input, char *channel_name, GstPad *pad, guint64 length, GstClockTime time, guint rate, void *chi2, gsl_matrix_complex_float_view *snr_matrix_view, int channel, GstBuffer *srcbuf); G_END_DECLS #endif /* __GSTLAL_SNGLTRIGGER_H__ */