...
 
Commits (64)
......@@ -214,9 +214,13 @@ test:gstlal:
script:
- export GSTLAL_FIR_WHITEN=0
## Get the necessary ROM data:
#- git clone https://git.ligo.org/lscsoft/lalsuite-extra.git ${GSTLAL_DIR}/lalsuite-extra
#- export LAL_DATA_PATH=${GSTLAL_DIR}/lalsuite-extra/data/lalsimulation/
# Get the necessary ROM data:
- git clone https://git.ligo.org/lscsoft/lalsuite-extra.git ${GSTLAL_DIR}/lalsuite-extra
- export LAL_DATA_PATH=${GSTLAL_DIR}/lalsuite-extra/data/lalsimulation/
- git clone https://git.ligo.org/alexander-pace/gstlal-testing-data.git ${GSTLAL_DIR}/gstlal-testing-data
- export LAL_DATA_PATH=${GSTLAL_DIR}/gstlal-testing-data/
# Run doctests
- cd gstlal
......@@ -269,8 +273,6 @@ test:offline:
#- yum -y install texlive-base texlive-dvipng
#- yum -y install texlive-latex texlive-latex-base-doc texlive-latex-fonts texlive-latex-bin texlive-latex-bin-bin
#- yum -y install texlive-pictures texlive-pstricks texlive-pstricks-doc
- yum -y install texlive*
# Get the necessary ROM data:
- git clone https://git.ligo.org/alexander-pace/gstlal-testing-data.git ${GSTLAL_DIR}/gstlal-testing-data
......
......@@ -135,9 +135,15 @@ class PSDHandler(simplehandler.Handler):
firbank.set_property("fir_matrix", template_t)
self.firbank = firbank
else:
template, _ = lalsimulation.GenerateStringCusp(1.0,30,1.0/options.sample_rate)
firbank.set_property("fir_matrix",numpy.zeros((template.data.length,len(template_bank_table)), dtype=float))
firbank.set_property("latency",-(template.data.length - 1) // 2)
# use templates with all zeros during burn-in period, that way we won't get any triggers.
template = [None] * len(template_bank_table)
for i, row in enumerate(template_bank_table):
template[i], _ = lalsimulation.GenerateStringCusp(1.0,30,1.0/options.sample_rate)
template[i] = lal.ResizeREAL8TimeSeries(template[i], -int(32*options.sample_rate - template[i].data.length)//2 ,int(32*options.sample_rate))
template[i] = template[i].data.data
template[i] *= 0.0
firbank.set_property("latency",-(len(template[0]) - 1) // 2)
firbank.set_property("fir_matrix", template)
self.firbank = firbank
return True
return False
......@@ -203,7 +209,7 @@ template_bank_table = lsctables.SnglBurstTable.get_table(xmldoc)
# filter bank
#
head = firbank = pipeparts.mkfirbank(pipeline, head, time_domain = False, block_stride = 4 * options.sample_rate)
head = firbank = pipeparts.mkfirbank(pipeline, head, fir_matrix = numpy.zeros((len(template_bank_table),int(32*options.sample_rate)),dtype=numpy.float64), block_stride = 4 * options.sample_rate)
#
......
......@@ -139,7 +139,7 @@ class StreamListener(object):
logger.info('processing features for timestamp %f, latency = %.3f s' % (self.timestamp, latency))
### format row
percent_missed = 100. * ((self.num_channels - len(features.keys())) / self.num_channels)
percent_missed = 100 * (float(self.num_channels - len(features.keys())) / self.num_channels)
if features.has_key(self.target_channel):
target_snr = features[self.target_channel][0]['snr']
else:
......@@ -151,7 +151,10 @@ class StreamListener(object):
c = self.conn.cursor()
sql = ''' INSERT INTO fxmonitor(timestamp,latency,percent_missed,target_snr) VALUES(?,?,?,?) '''
c.execute(sql,data)
self.conn.commit()
try:
self.conn.commit()
except sqlite3.OperationalError:
pass
def start(self):
"""
......
......@@ -44,6 +44,7 @@
#include <gstlal_string_triggergen.h>
#include <gstlal/gstaudioadapter.h>
#include <gstlal/gstlal_debug.h>
/*
......@@ -290,6 +291,8 @@ static gboolean set_caps(GstBaseTransform *trans, GstCaps *incaps, GstCaps *outc
GSTLALStringTriggergen *element = GSTLAL_STRING_TRIGGERGEN(trans);
gboolean success = gst_audio_info_from_caps(&element->audio_info, incaps);
g_object_set(element->adapter, "unit-size", GST_AUDIO_INFO_WIDTH(&element->audio_info) / 8, NULL);
return success;
}
......@@ -443,6 +446,8 @@ static void finalize(GObject *object)
element->instrument = NULL;
g_free(element->channel_name);
element->channel_name = NULL;
gst_audioadapter_clear(element->adapter);
g_object_unref(element->adapter);
G_OBJECT_CLASS(gstlal_string_triggergen_parent_class)->finalize(object);
}
......@@ -543,6 +548,7 @@ static void gstlal_string_triggergen_class_init(GSTLALStringTriggergenClass *kla
static void gstlal_string_triggergen_init(GSTLALStringTriggergen *element)
{
g_mutex_init(&element->bank_lock);
element->adapter = g_object_new(GST_TYPE_AUDIOADAPTER, NULL);
element->bank_filename = NULL;
element->bank = NULL;
element->instrument = NULL;
......
......@@ -6,6 +6,7 @@
#include <gst/gst.h>
#include <gst/audio/audio.h>
#include <gst/base/gstbasetransform.h>
#include <gstlal/gstaudioadapter.h>
#include <lal/LIGOMetadataTables.h>
G_BEGIN_DECLS
......@@ -31,6 +32,8 @@ typedef struct {
typedef struct {
GstBaseTransform element;
GstAudioAdapter *adapter;
/*
* input stream
*/
......
[InputConfigurations]
# Filters file containing calibration FIR filters, relative to the directory gstlal-calibration/tests/check_calibration/, from which the pipeline is expected to be run
FiltersFileName: H1DCS_newsrcline_1173225472.npz
FiltersFileName: H1DCS_newsrcline_1173225472_v2.npz
# Data source should be set to frames or lvshm
DataSource: frames
FileChecksum: No
......@@ -51,6 +51,8 @@ PipelineGraphFilename: gstlal_compute_strain
Verbose: Yes
# Turn this on to write data presentation timestamps and real-time unix timestamps to file at the beginning and end of the pipeline, to measure latency
TestLatency: No
# Turn this on to compute transfer functions for the filters by comparing output data to input data
TestFilters: Yes
[TDCFConfigurations]
#########################################################
......@@ -76,9 +78,6 @@ ApplyKappaUIM: No
# Set this to have the \kappa_u factors the actuation chain with an adaptive filter that corrects for both magnitude and phase errors.
ApplyComplexKappaUIM: No
# Set this to use a calibration line injected using the UIM stage of actuation to compute \kappa_U. Otherwise, the DARM_ctrl line is used.
UseUIMLine: Yes
ComputeKappaC: Yes
ApplyKappaC: Yes
......@@ -101,8 +100,6 @@ CoherenceTime: 130
###################################################################
# Options related to the computation configurations for the TDCFs #
###################################################################
ComputeFactorsSR: 16
RecordFactorsSR: 16
# Length in seconds of low-pass FIR filter used in demodulation of the calibration lines
DemodulationFilterTime: 20
# Time (in seconds) to smooth out \kappas with a median-like method
......@@ -110,6 +107,8 @@ MedianSmoothingTime: 128
TDCFAveragingTime: 10
#If set to yes, bad computed kappas will be replaced by the previous computed median in the running median array. Otherwise, they are replaced with the default value
TDCFDefaultToMedian: Yes
# If using X-end Pcal, we need a minus sign, so set this to -1.0
PcalSign: 1.0
##################################################
# Options related to updating cavity pole filter #
##################################################
......@@ -181,10 +180,12 @@ DeltaLTSTChannel: CAL-DELTAL_CTRL_TST_DBL_DQ
DeltaLPUMChannel: CAL-DELTAL_CTRL_PUM_DBL_DQ
DeltaLUIMChannel: CAL-DELTAL_CTRL_UIM_DBL_DQ
DeltaLResChannel: CAL-DELTAL_RESIDUAL_DBL_DQ
####################################
# Data Quality Vector Channel Name #
####################################
InputDQChannel: ODC-MASTER_CHANNEL_OUT_DQ
#####################################
# Data Quality Vector Channel Names #
#####################################
ObsIntentChannel: ODC-MASTER_CHANNEL_OUT_DQ
LowNoiseStateChannel: GRD-ISC_LOCK_OK
HWInjChannel: ODC-MASTER_CHANNEL_OUT_DQ
##################################
# Calibration Line Channel Names #
##################################
......@@ -193,6 +194,17 @@ TSTExcChannel: SUS-ETMY_L3_CAL_LINE_OUT_DQ
PUMExcChannel: SUS-ETMY_L2_CAL_LINE_OUT_DQ
UIMExcChannel: SUS-ETMY_L1_CAL_LINE_OUT_DQ
PCALChannel: CAL-PCALY_TX_PD_OUT_DQ
############################################
# Calibration Line Frequency Channel Names #
############################################
DARMExcLineFreqChannel: None
TSTExcLineFreqChannel: SUS-ETMY_L3_CAL_LINE_FREQ
PUMExcLineFreqChannel: SUS-ETMY_L2_CAL_LINE_FREQ
UIMExcLineFreqChannel: SUS-ETMY_L1_CAL_LINE_FREQ
PCALLine1FreqChannel: None
PCALLine2FreqChannel: CAL-CS_TDEP_CAVITY_POLE_PCAL_FREQ
PCALLine3FreqChannel: None
PCALLine4FreqChannel: CAL-CS_TDEP_D2N_SPRING_PCAL_FREQ
#######################################
# Coherence Uncertainty Channel Names #
#######################################
......@@ -208,7 +220,9 @@ CohUncDARMLine1Channel: CAL-CS_TDEP_DARM_LINE1_UNCERTAINTY
PowerLinesChannel: PEM-EY_MAINSMON_EBAY_1_DQ
# Comma-separated list of witness channels to use to subtract noise from h(t)
# Set to None if no witness channels are to be used
WitnessChannelList: IMC-WFS_A_DC_PIT_OUT_DQ,IMC-WFS_B_DC_PIT_OUT_DQ,IMC-WFS_A_DC_YAW_OUT_DQ,IMC-WFS_B_DC_YAW_OUT_DQ;PSL-DIAG_BULLSEYE_YAW_OUT_DQ,PSL-DIAG_BULLSEYE_WID_OUT_DQ,PSL-DIAG_BULLSEYE_PIT_OUT_DQ;ASC-DHARD_P_OUT_DQ,ASC-DHARD_Y_OUT_DQ,ASC-CHARD_P_OUT_DQ,ASC-CHARD_Y_OUT_DQ;LSC-SRCL_IN1_DQ,LSC-MICH_IN1_DQ,LSC-PRCL_IN1_DQ
WitnessChannelList: IMC-WFS_A_DC_PIT_OUT_DQ,IMC-WFS_B_DC_PIT_OUT_DQ,IMC-WFS_A_DC_YAW_OUT_DQ,IMC-WFS_B_DC_YAW_OUT_DQ,PSL-DIAG_BULLSEYE_YAW_OUT_DQ,PSL-DIAG_BULLSEYE_WID_OUT_DQ,PSL-DIAG_BULLSEYE_PIT_OUT_DQ;ASC-DHARD_P_OUT_DQ,ASC-DHARD_Y_OUT_DQ,ASC-CHARD_P_OUT_DQ,ASC-CHARD_Y_OUT_DQ,LSC-SRCL_IN1_DQ,LSC-MICH_IN1_DQ,LSC-PRCL_IN1_DQ
# What channel should we use to gate the noise subtraction and 60-Hz line subtraction
NoiseSubGateChannel: GRD-ISC_LOCK_OK
###############################
# EPICS Records Channel Names #
###############################
......@@ -269,8 +283,10 @@ CalibStateSR: 16
# Sample rate of control channel
# Should be 16384 if using DARM_CTRL and 4096 if using DELTAL_CTRL
CtrlSR: 16384
# Sample rate of ODC channel
ODCSR: 16384
# Sample rate of input DQ channels
ObsIntentSR: 16384
LowNoiseSR: 16
HWInjSR: 16384
# Sample rate of TST excitation channel
TSTExcSR: 512
# Sample rate of PUM excitation channel
......@@ -284,19 +300,19 @@ EPICSRefSR: 16
# Sample rate for power lines channel
PowerLinesChannelSR: 1024
# Sample rates at which transfer functions will be computed and witness channels will be filtered, given as a semicolon-separated list, e.g., 2048;2048;512;2048. This must be given if WitnessChannelList is not None, and it must be the same length.
WitnessChannelSR: 2048;2048;512;2048
WitnessChannelSR: 2048;512
# Sample rates at which to compute and record TDCFs
ComputeFactorsSR: 16
RecordFactorsSR: 16
[Bitmasks]
ObsReadyBitmask: 4
LowNoiseBitmask: 1
ObsIntentBitmask: 2
CBCHWInjBitmask: 16777216
BurstHWInjBitmask: 33554432
DetCharHWInjBitmask: 67108864
StochHWInjBitmask: 8388608
NoiseSubGateBitmask: 2
NoiseSubGateBitmask: 1
[PipelineConfigurations]
BufferLength: 1.0
......@@ -306,9 +322,9 @@ Dewhitening: No
FilterLatency: 1.0
[DataCleaningConfigurations]
####################################################
###################################################
# Options for turning on and off line subtraction #
####################################################
###################################################
# Remove the DC component from the residual and control channels before filtering
RemoveDC: No
# Subtract the calibration lines from the h(t) spectrum
......@@ -321,6 +337,7 @@ RemovePowerLines: Yes
# Amount by which frequency of power lines varies with time
PowerLinesFreqVar: 0.02
# Time over which to average the transfer function between the power mains witness channel and h(t) at 60 Hz and harmonics
PowerLinesTFMedianTime: 1
PowerLinesTFAveragingTime: 128
#######################################
# Options for broadband noise removal #
......@@ -332,11 +349,13 @@ NumWitnessFFTs: 509
# Sets the minimum number of FFTs necessary to produce the first transfer functions and clean data after data flow starts.
MinWitnessFFTs: 509
# The length in seconds of the filters applied to the witness channels before subtracting from h(t)
WitnessFIRLength: 0.5
WitnessFIRLength: 0.5;1.6
# The frequency resolution of the filters applied to the witness channels before subtracting from h(t). It can be advantageous to lower the frequency resolution in order to average over excess noise.
WitnessFrequencyResolution: 1.0
WitnessFrequencyResolution: 0.5
# List of minima and maxima of frequency ranges where the Fourier transform of h(t) will be replaced by a straight line in the calculation of transfer functions between witness channels and h(t) for noise subtraction. Semicolons separate lists for different sets of witness channels. If no notches are desired, use zeros, e.g., \'0;0;0\'. Here is an example using the expected format: \'495.0,515.0,985.0,1015.0;59,60,119,121;0\' This can be useful, e.g., if there are loud lines in the signal that are not present in the witness channels.
WitnessNotchFrequencies: 12.0,15.0,495.0,515.0,985.0,1015.0;12.0,15.0,495.0,515.0,985.0,1015.0;12.0,15.0,495.0,515.0,985.0,1015.0;12.0,15.0,495.0,515.0,985.0,1015.0
WitnessNotchFrequencies: 12.0,15.0,495.0,515.0,985.0,1015.0;0
# Cutoff frequencies for high-pass filters for witness channels
WitnessHighPasses: 12;10
# The amount of time after transfer functions between witness channels and h(t) are finished to begin the calculation of the next set of transfer functions
WitnessTFUpdateTime: 4
# If lock-loss lasts at least this many seconds, transfer functions will revert to those computed at the beginning of a lock stretch. Set to zero to disbale.
......
[InputConfigurations]
# Filters file containing calibration FIR filters, relative to the directory gstlal-calibration/tests/check_calibration/, from which the pipeline is expected to be run
FiltersFileName: H1DCS_newsrcline_1173225472.npz
FiltersFileName: H1DCS_newsrcline_1173225472_v2.npz
# Data source should be set to frames or lvshm
DataSource: frames
FileChecksum: No
......@@ -51,6 +51,8 @@ PipelineGraphFilename: gstlal_compute_strain
Verbose: Yes
# Turn this on to write data presentation timestamps and real-time unix timestamps to file at the beginning and end of the pipeline, to measure latency
TestLatency: No
# Turn this on to compute transfer functions for the filters by comparing output data to input data
TestFilters: Yes
[TDCFConfigurations]
#########################################################
......@@ -76,9 +78,6 @@ ApplyKappaUIM: No
# Set this to have the \kappa_u factors the actuation chain with an adaptive filter that corrects for both magnitude and phase errors.
ApplyComplexKappaUIM: No
# Set this to use a calibration line injected using the UIM stage of actuation to compute \kappa_U. Otherwise, the DARM_ctrl line is used.
UseUIMLine: Yes
ComputeKappaC: Yes
ApplyKappaC: Yes
......@@ -101,8 +100,6 @@ CoherenceTime: 130
###################################################################
# Options related to the computation configurations for the TDCFs #
###################################################################
ComputeFactorsSR: 16
RecordFactorsSR: 16
# Length in seconds of low-pass FIR filter used in demodulation of the calibration lines
DemodulationFilterTime: 20
# Time (in seconds) to smooth out \kappas with a median-like method
......@@ -110,6 +107,8 @@ MedianSmoothingTime: 128
TDCFAveragingTime: 10
#If set to yes, bad computed kappas will be replaced by the previous computed median in the running median array. Otherwise, they are replaced with the default value
TDCFDefaultToMedian: Yes
# If using X-end Pcal, we need a minus sign, so set this to -1.0
PcalSign: 1.0
##################################################
# Options related to updating cavity pole filter #
##################################################
......@@ -193,6 +192,17 @@ TSTExcChannel: SUS-ETMY_L3_CAL_LINE_OUT_DQ
PUMExcChannel: SUS-ETMY_L2_CAL_LINE_OUT_DQ
UIMExcChannel: SUS-ETMY_L1_CAL_LINE_OUT_DQ
PCALChannel: CAL-PCALY_TX_PD_OUT_DQ
############################################
# Calibration Line Frequency Channel Names #
############################################
DARMExcLineFreqChannel: None
TSTExcLineFreqChannel: SUS-ETMY_L3_CAL_LINE_FREQ
PUMExcLineFreqChannel: SUS-ETMY_L2_CAL_LINE_FREQ
UIMExcLineFreqChannel: SUS-ETMY_L1_CAL_LINE_FREQ
PCALLine1FreqChannel: None
PCALLine2FreqChannel: CAL-CS_TDEP_CAVITY_POLE_PCAL_FREQ
PCALLine3FreqChannel: None
PCALLine4FreqChannel: CAL-CS_TDEP_D2N_SPRING_PCAL_FREQ
#######################################
# Coherence Uncertainty Channel Names #
#######################################
......@@ -208,7 +218,9 @@ CohUncDARMLine1Channel: CAL-CS_TDEP_DARM_LINE1_UNCERTAINTY
PowerLinesChannel: PEM-EY_MAINSMON_EBAY_1_DQ
# Comma-separated list of witness channels to use to subtract noise from h(t)
# Set to None if no witness channels are to be used
WitnessChannelList: IMC-WFS_A_DC_PIT_OUT_DQ,IMC-WFS_B_DC_PIT_OUT_DQ,IMC-WFS_A_DC_YAW_OUT_DQ,IMC-WFS_B_DC_YAW_OUT_DQ;PSL-DIAG_BULLSEYE_YAW_OUT_DQ,PSL-DIAG_BULLSEYE_WID_OUT_DQ,PSL-DIAG_BULLSEYE_PIT_OUT_DQ;ASC-DHARD_P_OUT_DQ,ASC-DHARD_Y_OUT_DQ,ASC-CHARD_P_OUT_DQ,ASC-CHARD_Y_OUT_DQ;LSC-SRCL_IN1_DQ,LSC-MICH_IN1_DQ,LSC-PRCL_IN1_DQ
WitnessChannelList: IMC-WFS_A_DC_PIT_OUT_DQ,IMC-WFS_B_DC_PIT_OUT_DQ,IMC-WFS_A_DC_YAW_OUT_DQ,IMC-WFS_B_DC_YAW_OUT_DQ,PSL-DIAG_BULLSEYE_YAW_OUT_DQ,PSL-DIAG_BULLSEYE_WID_OUT_DQ,PSL-DIAG_BULLSEYE_PIT_OUT_DQ;ASC-DHARD_P_OUT_DQ,ASC-DHARD_Y_OUT_DQ,ASC-CHARD_P_OUT_DQ,ASC-CHARD_Y_OUT_DQ,LSC-SRCL_IN1_DQ,LSC-MICH_IN1_DQ,LSC-PRCL_IN1_DQ
# What channel should we use to gate the noise subtraction and 60-Hz line subtraction
NoiseSubGateChannel: ODC-MASTER_CHANNEL_OUT_DQ
###############################
# EPICS Records Channel Names #
###############################
......@@ -284,7 +296,7 @@ EPICSRefSR: 16
# Sample rate for power lines channel
PowerLinesChannelSR: 1024
# Sample rates at which transfer functions will be computed and witness channels will be filtered, given as a semicolon-separated list, e.g., 2048;2048;512;2048. This must be given if WitnessChannelList is not None, and it must be the same length.
WitnessChannelSR: 2048;2048;512;2048
WitnessChannelSR: 2048;512
# Sample rates at which to compute and record TDCFs
ComputeFactorsSR: 16
RecordFactorsSR: 16
......@@ -296,7 +308,7 @@ CBCHWInjBitmask: 16777216
BurstHWInjBitmask: 33554432
DetCharHWInjBitmask: 67108864
StochHWInjBitmask: 8388608
NoiseSubGateBitmask: 2
NoiseSubGateBitmask: 4
[PipelineConfigurations]
BufferLength: 1.0
......@@ -306,9 +318,9 @@ Dewhitening: No
FilterLatency: 1.0
[DataCleaningConfigurations]
####################################################
###################################################
# Options for turning on and off line subtraction #
####################################################
###################################################
# Remove the DC component from the residual and control channels before filtering
RemoveDC: No
# Subtract the calibration lines from the h(t) spectrum
......@@ -321,6 +333,7 @@ RemovePowerLines: Yes
# Amount by which frequency of power lines varies with time
PowerLinesFreqVar: 0.02
# Time over which to average the transfer function between the power mains witness channel and h(t) at 60 Hz and harmonics
PowerLinesTFMedianTime: 1
PowerLinesTFAveragingTime: 128
#######################################
# Options for broadband noise removal #
......@@ -332,11 +345,13 @@ NumWitnessFFTs: 509
# Sets the minimum number of FFTs necessary to produce the first transfer functions and clean data after data flow starts.
MinWitnessFFTs: 509
# The length in seconds of the filters applied to the witness channels before subtracting from h(t)
WitnessFIRLength: 0.5
WitnessFIRLength: 0.5;1.6
# The frequency resolution of the filters applied to the witness channels before subtracting from h(t). It can be advantageous to lower the frequency resolution in order to average over excess noise.
WitnessFrequencyResolution: 1.0
WitnessFrequencyResolution: 0.5
# List of minima and maxima of frequency ranges where the Fourier transform of h(t) will be replaced by a straight line in the calculation of transfer functions between witness channels and h(t) for noise subtraction. Semicolons separate lists for different sets of witness channels. If no notches are desired, use zeros, e.g., \'0;0;0\'. Here is an example using the expected format: \'495.0,515.0,985.0,1015.0;59,60,119,121;0\' This can be useful, e.g., if there are loud lines in the signal that are not present in the witness channels.
WitnessNotchFrequencies: 12.0,15.0,495.0,515.0,985.0,1015.0;12.0,15.0,495.0,515.0,985.0,1015.0;12.0,15.0,495.0,515.0,985.0,1015.0;12.0,15.0,495.0,515.0,985.0,1015.0
WitnessNotchFrequencies: 12.0,15.0,495.0,515.0,985.0,1015.0;0
# Cutoff frequencies for high-pass filters for witness channels
WitnessHighPasses: 12;10
# The amount of time after transfer functions between witness channels and h(t) are finished to begin the calculation of the next set of transfer functions
WitnessTFUpdateTime: 4
# If lock-loss lasts at least this many seconds, transfer functions will revert to those computed at the beginning of a lock stretch. Set to zero to disbale.
......
......@@ -17,7 +17,8 @@ lib@GSTPLUGINPREFIX@gstlalcalibration_la_SOURCES = \
gstlal_transferfunction.c gstlal_transferfunction.h \
gstlal_trackfrequency.c gstlal_trackfrequency.h \
gstlal_adaptivefirfilt.c gstlal_adaptivefirfilt.h \
gstlal_dqtukey.c gstlal_dqtukey.h
gstlal_dqtukey.c gstlal_dqtukey.h \
gstlal_property.c gstlal_property.h
lib@GSTPLUGINPREFIX@gstlalcalibration_la_CPPFLAGS = $(AM_CPPFLAGS) $(PYTHON_CPPFLAGS)
lib@GSTPLUGINPREFIX@gstlalcalibration_la_CFLAGS = $(AM_CFLAGS) $(LAL_CFLAGS) $(GSTLAL_CFLAGS) $(gstreamer_CFLAGS) $(gstreamer_audio_CFLAGS)
lib@GSTPLUGINPREFIX@gstlalcalibration_la_LDFLAGS = $(AM_LDFLAGS) $(LAL_LIBS) $(GSTLAL_LIBS) $(PYTHON_LIBS) $(gstreamer_LIBS) $(gstreamer_audio_LIBS) $(GSTLAL_PLUGIN_LDFLAGS)
......@@ -67,59 +67,59 @@
*/
static void demodulate_float(const float *src, gsize src_size, float complex *dst, guint64 t, gint rate, const int frequency, complex float prefactor)
static void demodulate_float(const float *src, gsize src_size, float complex *dst, guint64 t, gint rate, const gint64 frequency, complex float prefactor)
{
const float *src_end;
guint64 i = 0;
__int128_t t_scaled, integer_arg, scale = 32000000000ULL;
__int128_t t_scaled, integer_arg, scale = 128000000000ULL;
for(src_end = src + src_size; src < src_end; src++, dst++, i++) {
t_scaled = 32 * t + i * scale / rate;
integer_arg = (t_scaled * frequency) % (100 * scale);
*dst = prefactor * *src * cexpf(-2. * M_PI * I * integer_arg / (100.0 * scale));
t_scaled = 128 * t + i * scale / rate;
integer_arg = (t_scaled * frequency) % (1000000 * scale);
*dst = prefactor * *src * cexpf(-2. * M_PI * I * integer_arg / (1000000.0 * scale));
}
}
static void demodulate_double(const double *src, gsize src_size, double complex *dst, guint64 t, gint rate, const int frequency, complex double prefactor)
static void demodulate_double(const double *src, gsize src_size, double complex *dst, guint64 t, gint rate, const gint64 frequency, complex double prefactor)
{
const double *src_end;
guint64 i = 0;
__int128_t t_scaled, integer_arg, scale = 32000000000ULL;
__int128_t t_scaled, integer_arg, scale = 128000000000ULL;
for(src_end = src + src_size; src < src_end; src++, dst++, i++) {
t_scaled = 32 * t + i * scale / rate;
integer_arg = (t_scaled * frequency) % (100 * scale);
*dst = prefactor * *src * cexp(-2. * M_PI * I * integer_arg / (100.0 * scale));
t_scaled = 128 * t + i * scale / rate;
integer_arg = (t_scaled * frequency) % (1000000 * scale);
*dst = prefactor * *src * cexp(-2. * M_PI * I * integer_arg / (1000000.0 * scale));
}
}
static void demodulate_complex_float(const complex float *src, gsize src_size, float complex *dst, guint64 t, gint rate, const int frequency, complex float prefactor)
static void demodulate_complex_float(const complex float *src, gsize src_size, float complex *dst, guint64 t, gint rate, const gint64 frequency, complex float prefactor)
{
const complex float *src_end;
guint64 i = 0;
__int128_t t_scaled, integer_arg, scale = 32000000000ULL;
__int128_t t_scaled, integer_arg, scale = 128000000000ULL;
for(src_end = src + src_size; src < src_end; src++, dst++, i++) {
t_scaled = 32 * t + i * scale / rate;
integer_arg = (t_scaled * frequency) % (100 * scale);
*dst = prefactor * *src * cexpf(-2. * M_PI * I * integer_arg / (100.0 * scale));
t_scaled = 128 * t + i * scale / rate;
integer_arg = (t_scaled * frequency) % (1000000 * scale);
*dst = prefactor * *src * cexpf(-2. * M_PI * I * integer_arg / (1000000.0 * scale));
}
}
static void demodulate_complex_double(const complex double *src, gsize src_size, double complex *dst, guint64 t, gint rate, const int frequency, complex double prefactor)
static void demodulate_complex_double(const complex double *src, gsize src_size, double complex *dst, guint64 t, gint rate, const gint64 frequency, complex double prefactor)
{
const complex double *src_end;
guint64 i = 0;
__int128_t t_scaled, integer_arg, scale = 32000000000ULL;
__int128_t t_scaled, integer_arg, scale = 128000000000ULL;
for(src_end = src + src_size; src < src_end; src++, dst++, i++) {
t_scaled = 32 * t + i * scale / rate;
integer_arg = (t_scaled * frequency) % (100 * scale);
*dst = prefactor * *src * cexp(-2. * M_PI * I * integer_arg / (100.0 * scale));
t_scaled = 128 * t + i * scale / rate;
integer_arg = (t_scaled * frequency) % (1000000 * scale);
*dst = prefactor * *src * cexp(-2. * M_PI * I * integer_arg / (1000000.0 * scale));
}
}
static void demodulate(const void *src, gsize src_size, void *dst, guint64 t, gint rate, enum gstlal_demodulate_data_type data_type, const int frequency, complex double prefactor)
static void demodulate(const void *src, gsize src_size, void *dst, guint64 t, gint rate, enum gstlal_demodulate_data_type data_type, const gint64 frequency, complex double prefactor)
{
switch(data_type) {
......@@ -581,8 +581,9 @@ static void set_property(GObject *object, enum property prop_id, const GValue *v
GST_OBJECT_LOCK(element);
switch (prop_id) {
case ARG_LINE_FREQUENCY:
element->line_frequency = (int) (100.000000000001 * g_value_get_double(value)); /* Make sure truncation does not corrupt it */
case ARG_LINE_FREQUENCY: ;
double freq = g_value_get_double(value);
element->line_frequency = (gint64) (1000000.0 * freq + (freq > 0 ? 0.5 : -0.5)); /* Round to the nearest uHz */
break;
case ARG_PREFACTOR_REAL:
element->prefactor_real = g_value_get_double(value);
......@@ -607,7 +608,7 @@ static void get_property(GObject *object, enum property prop_id, GValue *value,
switch (prop_id) {
case ARG_LINE_FREQUENCY:
g_value_set_double(value, element->line_frequency / 100.0);
g_value_set_double(value, element->line_frequency / 1000000.0);
break;
case ARG_PREFACTOR_REAL:
g_value_set_double(value, element->prefactor_real);
......
......@@ -79,7 +79,7 @@ struct _GSTLALDemodulate {
gboolean need_discont;
/* properties */
int line_frequency; /* in centihertz */
gint64 line_frequency; /* in uHz */
double prefactor_real;
double prefactor_imag;
};
......
This diff is collapsed.
/*
* Copyright (C) 2018 Aaron Viets <aaron.viets@ligo.org>
*
* This program is free software; you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by the
* Free Software Foundation; either version 2 of the License, or (at your
* option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License along
* with this program; if not, write to the Free Software Foundation, Inc.,
* 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
#ifndef __GSTLAL_PROPERTY_H__
#define __GSTLAL_PROPERTY_H__
#include <glib.h>
#include <gst/gst.h>
#include <gst/base/gstbasesink.h>
G_BEGIN_DECLS
#define GSTLAL_PROPERTY_TYPE \
(gstlal_property_get_type())
#define GSTLAL_PROPERTY(obj) \
(G_TYPE_CHECK_INSTANCE_CAST((obj), GSTLAL_PROPERTY_TYPE, GSTLALProperty))
#define GSTLAL_PROPERTY_CLASS(klass) \
(G_TYPE_CHECK_CLASS_CAST((klass), GSTLAL_PROPERTY_TYPE, GSTLALPropertyClass))
#define GST_IS_GSTLAL_PROPERTY(obj) \
(G_TYPE_CHECK_INSTANCE_TYPE((obj), GSTLAL_PROPERTY_TYPE))
#define GST_IS_GSTLAL_PROPERTY_CLASS(klass) \
(G_TYPE_CHECK_CLASS_TYPE((klass), GSTLAL_PROPERTY_TYPE))
typedef struct _GSTLALProperty GSTLALProperty;
typedef struct _GSTLALPropertyClass GSTLALPropertyClass;
/**
* GSTLALProperty:
*/
struct _GSTLALProperty {
GstBaseSink basesink;
/* stream info */
gint rate;
gint unit_size;
enum gstlal_property_data_type {
GSTLAL_PROPERTY_SIGNED = 0,
GSTLAL_PROPERTY_UNSIGNED,
GSTLAL_PROPERTY_FLOAT,
} data_type;
/* timestamp bookkeeping */
GstClockTime t0;
guint64 offset0;
guint64 next_in_offset;
/* filter memory */
gint64 num_in_avg;
/* properties */
gint64 update_samples;
gint64 shift_samples;
gint64 average_samples;
gboolean update_when_change;
double current_average;
};
/**
* GSTLALPropertyClass:
* @parent_class: the parent class
*/
struct _GSTLALPropertyClass {
GstBaseSinkClass parent_class;
};
GType gstlal_property_get_type(void);
G_END_DECLS
#endif /* __GSTLAL_PROPERTY_H__ */
......@@ -69,6 +69,7 @@
#include <gstlal_trackfrequency.h>
#include <gstlal_adaptivefirfilt.h>
#include <gstlal_dqtukey.h>
#include <gstlal_property.h>
/*
......@@ -102,6 +103,7 @@ static gboolean plugin_init(GstPlugin *plugin)
{"lal_trackfrequency", GSTLAL_TRACKFREQUENCY_TYPE},
{"lal_adaptivefirfilt", GSTLAL_ADAPTIVEFIRFILT_TYPE},
{"lal_dqtukey", GSTLAL_DQTUKEY_TYPE},
{"lal_property", GSTLAL_PROPERTY_TYPE},
{NULL, 0},
};
......
......@@ -167,10 +167,12 @@ def list_srcs(pipeline, *args):
# Various filtering functions
#
def demodulate(pipeline, head, freq, td, rate, filter_time, filter_latency, prefactor_real = 1.0, prefactor_imag = 0.0):
def demodulate(pipeline, head, freq, td, rate, filter_time, filter_latency, prefactor_real = 1.0, prefactor_imag = 0.0, freq_update = None):
# demodulate input at a given frequency freq
head = pipeparts.mkgeneric(pipeline, head, "lal_demodulate", line_frequency = freq, prefactor_real = prefactor_real, prefactor_imag = prefactor_imag)
if freq_update is not None:
freq_update.connect("notify::current-average", update_property_simple, head, "current_average", "line_frequency")
head = mkresample(pipeline, head, 5, filter_latency == 0.0, rate)
if filter_latency != 0:
# Remove the first several seconds of output, which depend on start time
......@@ -198,7 +200,7 @@ def remove_harmonics(pipeline, signal, f0, num_harmonics, f0_var, filter_latency
return elem
def remove_harmonics_with_witness(pipeline, signal, witness, f0, num_harmonics, f0_var, filter_latency, compute_rate = 16, rate_out = 16384, num_avg = 2048, noisesub_gate_bit = None):
def remove_harmonics_with_witness(pipeline, signal, witness, f0, num_harmonics, f0_var, filter_latency, compute_rate = 16, rate_out = 16384, num_median = 2048, num_avg = 160, noisesub_gate_bit = None):
# remove line(s) from a spectrum. filter length for demodulation (given in seconds) is adjustable
# function argument caps must be complex caps
......@@ -260,7 +262,7 @@ def remove_harmonics_with_witness(pipeline, signal, witness, f0, num_harmonics,
# Remove worthless data from computation of transfer function if we can
if noisesub_gate_bit is not None:
tf_at_f = mkgate(pipeline, tf_at_f, noisesub_gate_bit, 1, attack_length = -((1.0 - filter_latency) * filter_samples))
tf_at_f = pipeparts.mkgeneric(pipeline, tf_at_f, "lal_smoothkappas", default_kappa_re = 0, array_size = 1, avg_array_size = num_avg, default_to_median = True, filter_latency = filter_latency)
tf_at_f = pipeparts.mkgeneric(pipeline, tf_at_f, "lal_smoothkappas", default_kappa_re = 0, array_size = num_median, avg_array_size = num_avg, default_to_median = True, filter_latency = filter_latency)
# Use gated, averaged transfer function to reconstruct the sinusoid as it appears in the signal from the witness channel
if filter_latency == 0.5:
......@@ -815,13 +817,17 @@ def compute_kappac(pipeline, SR, SI):
kc = mkmultiplier(pipeline, list_srcs(pipeline, S2, mkpow(pipeline, SR, exponent=-1.0)))
return kc
def compute_fcc(pipeline, SR, SI, fpcal2):
def compute_fcc(pipeline, SR, SI, fpcal2, freq_update = None):
#
# f_cc = - (Re[S]/Im[S]) * fpcal2
#
fcc = mkmultiplier(pipeline, list_srcs(pipeline, pipeparts.mkaudioamplify(pipeline, SR, -1.0*fpcal2), mkpow(pipeline, SI, exponent=-1.0)))
fcc = mkmultiplier(pipeline, list_srcs(pipeline, pipeparts.mkaudioamplify(pipeline, SR, -1.0), mkpow(pipeline, SI, exponent=-1.0)))
fcc = pipeparts.mkaudioamplify(pipeline, fcc, fpcal2)
if freq_update is not None:
freq_update.connect("notify::current-average", update_property_simple, fcc, "current_average", "amplification")
return fcc
def compute_Xi_from_filters_file(pipeline, pcalfpcal4, darmfpcal4, fpcal4, EP11_real, EP11_imag, EP12_real, EP12_imag, EP13_real, EP13_imag, EP14_real, EP14_imag, ktst, kpu, kc, fcc):
......@@ -926,7 +932,7 @@ def update_filters(filter_maker, arg, filter_taker, maker_prop_name, taker_prop_
firfilter = filter_maker.get_property(maker_prop_name)[filter_number][::-1]
filter_taker.set_property(taker_prop_name, firfilter)
def clean_data(pipeline, signal, signal_rate, witnesses, witness_rate, fft_length, fft_overlap, num_ffts, min_ffts, update_samples, fir_length, frequency_resolution, filter_taper_length, use_median = False, parallel_mode = False, notch_frequencies = [], noisesub_gate_bit = None, delay_time = 0.0, critical_lock_loss_time = 0, filename = None):
def clean_data(pipeline, signal, signal_rate, witnesses, witness_rate, fft_length, fft_overlap, num_ffts, min_ffts, update_samples, fir_length, frequency_resolution, filter_taper_length, use_median = False, parallel_mode = False, notch_frequencies = [], high_pass = 15.0, noisesub_gate_bit = None, delay_time = 0.0, critical_lock_loss_time = 0, filename = None):
#
# Use witness channels that monitor the environment to remove environmental noise
......@@ -939,22 +945,21 @@ def clean_data(pipeline, signal, signal_rate, witnesses, witness_rate, fft_lengt
witness_tees = []
for i in range(0, len(witnesses)):
witnesses[i] = mkresample(pipeline, witnesses[i], 5, False, witness_rate)
witnesses[i] = highpass(pipeline, witnesses[i], witness_rate)
witness_tees.append(pipeparts.mktee(pipeline, witnesses[i]))
resampled_signal = mkresample(pipeline, signal_tee, 5, False, witness_rate)
transfer_functions = mkinterleave(pipeline, numpy.insert(witness_tees, 0, resampled_signal, axis = 0))
if noisesub_gate_bit is not None:
transfer_functions = mkgate(pipeline, transfer_functions, noisesub_gate_bit, 1)
transfer_functions = pipeparts.mkgeneric(pipeline, transfer_functions, "lal_transferfunction", fft_length = fft_length, fft_overlap = fft_overlap, num_ffts = num_ffts, min_ffts = min_ffts, update_samples = update_samples, make_fir_filters = -1, fir_length = fir_length, frequency_resolution = frequency_resolution, high_pass = 15, update_after_gap = True, use_median = use_median, parallel_mode = parallel_mode, notch_frequencies = notch_frequencies, use_first_after_gap = critical_lock_loss_time * witness_rate, update_delay_samples = int(delay_time * witness_rate), fir_timeshift = 0, filename = filename)
transfer_functions = pipeparts.mkgeneric(pipeline, transfer_functions, "lal_transferfunction", fft_length = fft_length, fft_overlap = fft_overlap, num_ffts = num_ffts, min_ffts = min_ffts, update_samples = update_samples, make_fir_filters = -1, fir_length = fir_length, frequency_resolution = frequency_resolution, high_pass = high_pass / 2.0, update_after_gap = True, use_median = use_median, parallel_mode = parallel_mode, notch_frequencies = notch_frequencies, use_first_after_gap = critical_lock_loss_time * witness_rate, update_delay_samples = int(delay_time * witness_rate), fir_timeshift = 0, filename = filename)
signal_minus_noise = [signal_tee]
for i in range(0, len(witnesses)):
if parallel_mode:
minus_noise = pipeparts.mkgeneric(pipeline, mkqueue(pipeline, witness_tees[i]), "lal_tdwhiten", kernel = numpy.zeros(fir_length), latency = fir_length / 2, taper_length = filter_taper_length, kernel_endtime = 0)
minus_noise = pipeparts.mkgeneric(pipeline, mkqueue(pipeline, highpass(pipeline, witness_tees[i], witness_rate, fcut = high_pass)), "lal_tdwhiten", kernel = numpy.zeros(fir_length), latency = fir_length / 2, taper_length = filter_taper_length, kernel_endtime = 0)
transfer_functions.connect("notify::fir-filters", update_filters, minus_noise, "fir_filters", "kernel", i)
transfer_functions.connect("notify::fir-endtime", update_property_simple, minus_noise, "fir_endtime", "kernel_endtime")
else:
minus_noise = pipeparts.mkgeneric(pipeline, witness_tees[i], "lal_tdwhiten", kernel = numpy.zeros(fir_length), latency = fir_length / 2, taper_length = filter_taper_length)
minus_noise = pipeparts.mkgeneric(pipeline, highpass(pipeline, witness_tees[i], witness_rate, fcut = high_pass), "lal_tdwhiten", kernel = numpy.zeros(fir_length), latency = fir_length / 2, taper_length = filter_taper_length)
transfer_functions.connect("notify::fir-filters", update_filters, minus_noise, "fir_filters", "kernel", i)
signal_minus_noise.append(mkresample(pipeline, minus_noise, 5, False, signal_rate))
......
......@@ -52,31 +52,25 @@ plot.gca().plot(hoft_asd,label='GDS h(t) ASD')
if options.analyze_additional_hoft_channel:
plot.gca().plot(additional_hoft_asd,label='DCS h(t) ASD')
ax = plot.gca()
#ax.set_ylabel = 'Strain [Hz$^{-1/2}$]'
#ax.set_xlabel = 'Frequency [Hz]'
plot.set_ylabel('Strain [Hz $^{-1/2}$]', fontname = 'Times', fontsize = 18)
plot.set_xlabel('Frequency [Hz]', fontname = 'Times', fontsize = 18)
plot.add_legend([r'CALCS h(t) ASD', r'GDS h(t) ASD', r'DCS h(t) ASD'], loc='upper right', fontsize='small')
ax.set_ylabel('Strain [Hz $^{-1/2}$]', fontname = 'Times', fontsize = 18)
ax.set_xlabel('Frequency [Hz]', fontname = 'Times', fontsize = 18)
ax.legend([r'CALCS h(t) ASD', r'GDS h(t) ASD', r'DCS h(t) ASD'], loc='upper right', fontsize='small')
ax.set_xlim(0.5,8192)
ax.set_ylim(1e-24,1e-16)
ax.legend()
plot.save('%s_%s_%s_spectrum_comparison.png' % (options.ifo, start_time, end_time))
diff1 = calcs_asd / hoft_asd
if options.analyze_additional_hoft_channel:
diff2 = calcs_asd / additional_hoft_asd
diff3 = hoft_asd / additional_hoft_asd
plot = diff1.plot(label="ASD ratios", logy = False)
plot = diff1.plot(label="ASD ratios")
if options.analyze_additional_hoft_channel:
plot.gca().plot(diff2)
plot.gca().plot(diff3)
ax = plot.gca()
#ax.set_ylabel = 'Strain [Hz$^{-1/2}$]'
#ax.set_xlabel = 'Frequency [Hz]'
plot.set_ylabel('Strain [Hz $^{-1/2}$]', fontname = 'Times', fontsize = 18)
plot.set_xlabel('Frequency [Hz]', fontname = 'Times', fontsize = 18)
plot.add_legend([r'CALCS h(t) ASD / GDS h(t) ASD', r'CALCS h(t) ASD / DCS h(t) ASD', r'GDS h(t) ASD / DCS h(t) ASD'], loc='upper right', fontsize='small')
ax.set_ylabel('Strain [Hz $^{-1/2}$]', fontname = 'Times', fontsize = 18)
ax.set_xlabel('Frequency [Hz]', fontname = 'Times', fontsize = 18)
ax.legend([r'CALCS h(t) ASD / GDS h(t) ASD', r'CALCS h(t) ASD / DCS h(t) ASD', r'GDS h(t) ASD / DCS h(t) ASD'], loc='upper right', fontsize='small')
ax.set_xlim(10,5000)
ax.set_ylim(0.7, 1.3)
ax.legend()
plot.save('%s_%s_%s_ASD_residual.png' % (options.ifo, start_time, end_time))
......@@ -8,10 +8,10 @@ IFO = H
# determines where to look for filters files (e.g., O1, O2, O3, ER10, ER13, ER14, PreER10, PreER13, PreER14)
OBSRUN = O2
START = $(shell echo 1185767424 - 715 | bc)
START = $(shell echo 1185763328 - 715 | bc)
#1225967424
#$(shell echo 1185763328 - 700 | bc)
END = $(shell echo 1185771520 + 715 | bc)
END = $(shell echo 1185767424 + 715 | bc)
#1225968448
#$(shell echo 1185767424 + 700 | bc)
# 1185771520
......@@ -27,7 +27,7 @@ GDSTESTCONFIGS = ../../config_files/PreER13/H1/H1GDS_TEST_1225558818.ini
DCSTESTCONFIGS = ../../config_files/O2/H1/tests/H1DCS_AllCorrections_Cleaning_TEST.ini
GDSSHMCONFIGS = ../../config_files/PreER13/H1/tests/H1GDS_1222058826_shm2frames.ini
all: noise_subtraction_ASD_DCS_TEST
all: noise_subtraction_ASD_DCS noise_subtraction_tf_DCS filters_tf_DCS
###############################################
### These commands should change less often ###
......@@ -65,6 +65,9 @@ $(IFO)1_C01_frames.cache:
$(IFO)1_C02_frames.cache:
gw_data_find -o $(IFO) -t $(IFO)1_HOFT_C02 -s $(START) -e $(END) -l --url-type file > $@
$(IFO)1_clean_C02_frames.cache:
gw_data_find -o $(IFO) -t $(IFO)1_CLEANED_HOFT_C02 -s $(START) -e $(END) -l --url-type file > $@
$(IFO)1_hoft_GDS_frames.cache: $(IFO)1_easy_raw_frames.cache filters framesdir
GST_DEBUG=3 gstlal_compute_strain --gps-start-time $(START) --gps-end-time $(END) --frame-cache $(IFO)1_easy_raw_frames.cache --output-path Frames/$(OBSRUN)/$(IFO)1/GDS/ --frame-duration=64 --frames-per-file=1 --wings=0 --config-file $(GDSCONFIGS)
ls Frames/$(OBSRUN)/$(IFO)1/GDS/$(IFO)-$(IFO)1GDS-*.gwf | lalapps_path2cache > $@
......@@ -102,6 +105,9 @@ pcal_DCS_transfer_functions: $(IFO)1_easy_raw_frames.cache $(IFO)1_hoft_DCS_fram
lines_ratio_DCS: $(IFO)1_hoft_DCS_frames.cache
python demod_ratio_timeseries.py --ifo $(IFO)1 --gps-end-time $(PLOT_END) --gps-start-time $(PLOT_START) --denominator-frame-cache $(IFO)1_hoft_DCS_frames.cache --numerator-frame-cache $(IFO)1_hoft_DCS_frames.cache --denominator-channel-name 'DCS-CALIB_STRAIN' --numerator-channel-name 'DCS-CALIB_STRAIN_CLEAN' --frequencies '35.9,36.7,331.9,1083.7;60,120,180' --magnitude-ranges '0.0,0.1;0.0,1.0' --phase-ranges '-180.0,180.0;-180.0,180.0' --plot-titles '$(IFO)1 Calibration Line Subtraction;$(IFO)1 Power Mains Line Subtraction'
filters_tf_DCS: $(IFO)1_hoft_DCS_frames.cache
python plot_filters_transfer_function.py --tf-frequency-min 0.5 --tf-frequency-max 8192 --ratio-frequency-min 10 --ratio-frequency-max 8192 --ratio-magnitude-min 0.7 --ratio-magnitude-max 1.3 --tf-phase-min -180 --tf-phase-max 180 --ratio-phase-min -20 --ratio-phase-max 20
latency_test: $(IFO)1_hoft_GDS_SHM_frames.cache
python latency_plot.py --intime-file gstlal_compute_strain_timestamps_in.txt --outtime-file gstlal_compute_strain_timestamps_out.txt --plot-filename-prefix $(IFO)1GDS_latency --plot-title '$(IFO)1 Calibration Latency vs Time'
......@@ -114,17 +120,26 @@ GDS_over_CALCS: $(IFO)1_hoft_GDS_frames.cache $(IFO)1_easy_raw_frames.cache
GDS_over_C02: $(IFO)1_hoft_GDS_frames.cache $(IFO)1_C02_frames.cache
python plot_transfer_function.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --denominator-frame-cache $(IFO)1_C02_frames.cache --denominator-channel-name DCS-CALIB_STRAIN_C02 --denominator-name 'C02' --numerator-frame-cache-list $(IFO)1_hoft_GDS_frames.cache --numerator-channel-list GDS-CALIB_STRAIN --numerator-name 'GDS' --use-median --magnitude-min 0.7 --magnitude-max 1.3 --phase-min -20.0 --phase-max 20.0 --labels 'GDS / C02'
DCS_over_C02: $(IFO)1_hoft_DCS_FCC_frames.cache $(IFO)1_C02_frames.cache
python plot_transfer_function.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --denominator-frame-cache $(IFO)1_C02_frames.cache --denominator-channel-name DCS-CALIB_STRAIN_C02 --denominator-name 'C02' --numerator-frame-cache-list $(IFO)1_hoft_DCS_FCC_frames.cache --numerator-channel-list DCS-CALIB_STRAIN --numerator-name 'DCS' --use-median --magnitude-min 0.7 --magnitude-max 1.3 --phase-min -20.0 --phase-max 20.0 --labels 'DCS / C02'
DCS_over_C02: $(IFO)1_hoft_DCS_frames.cache $(IFO)1_C02_frames.cache
python plot_transfer_function.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --denominator-frame-cache $(IFO)1_C02_frames.cache --denominator-channel-name DCS-CALIB_STRAIN_C02 --denominator-name 'C02' --numerator-frame-cache-list $(IFO)1_hoft_DCS_frames.cache --numerator-channel-list DCS-CALIB_STRAIN --numerator-name 'DCS' --use-median --magnitude-min 0.7 --magnitude-max 1.3 --phase-min -20.0 --phase-max 20.0 --labels 'DCS / C02'
noise_subtraction_ASD_DCS: $(IFO)1_hoft_DCS_frames.cache
./ASD_comparison_plots --ifo $(IFO)1 --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --raw-frame-cache $(IFO)1_hoft_DCS_frames.cache --calcs-channel-name DCS-CALIB_STRAIN --hoft-frame-cache $(IFO)1_hoft_DCS_frames.cache --hoft-channel-name DCS-CALIB_STRAIN_CLEAN
noise_subtraction_ASD_DCH_DCS: $(IFO)1_hoft_DCS_frames.cache $(IFO)1_clean_C02_frames.cache
./ASD_comparison_plots --ifo $(IFO)1 --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --raw-frame-cache $(IFO)1_clean_C02_frames.cache --calcs-channel-name DCH-CLEAN_STRAIN_C02 --hoft-frame-cache $(IFO)1_hoft_DCS_frames.cache --hoft-channel-name DCS-CALIB_STRAIN_CLEAN
noise_subtraction_ASD_DCS_TEST: $(IFO)1_hoft_DCS_TEST_frames.cache
./ASD_comparison_plots --ifo $(IFO)1 --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --raw-frame-cache $(IFO)1_hoft_DCS_TEST_frames.cache --calcs-channel-name DCS-CALIB_STRAIN --hoft-frame-cache $(IFO)1_hoft_DCS_TEST_frames.cache --hoft-channel-name DCS-CALIB_STRAIN_CLEAN
noise_subtraction_tf_DCS: $(IFO)1_hoft_DCS_frames.cache
python plot_transfer_function.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --numerator-frame-cache $(IFO)1_hoft_DCS_frames.cache --denominator-frame-cache $(IFO)1_hoft_DCS_frames.cache --numerator-channel-name DCS-CALIB_STRAIN_CLEAN --denominator-channel-name DCS-CALIB_STRAIN --magnitude-min 0.0 --magnitude-max 1.5 --phase-min -20.0 --phase-max 20.0 --plot-title 'Noise Subtraction Transfer Function'
python plot_transfer_function.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --numerator-frame-cache $(IFO)1_hoft_DCS_frames.cache --denominator-frame-cache $(IFO)1_hoft_DCS_frames.cache --numerator-channel-list DCS-CALIB_STRAIN_CLEAN --denominator-channel-name DCS-CALIB_STRAIN --magnitude-min 0.0 --magnitude-max 1.5 --phase-min -20.0 --phase-max 20.0 --numerator-name 'CLEAN' --denominator-name 'STRAIN' --labels 'CALIB_STRAIN_CLEAN / CALIB_STRAIN' --use-median
noise_subtraction_tf_DCH_DCS: $(IFO)1_hoft_DCS_frames.cache $(IFO)1_clean_C02_frames.cache
python plot_transfer_function.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --numerator-frame-cache $(IFO)1_hoft_DCS_frames.cache --denominator-frame-cache $(IFO)1_clean_C02_frames.cache --numerator-channel-list DCS-CALIB_STRAIN_CLEAN --denominator-channel-name DCH-CLEAN_STRAIN_C02 --magnitude-min 0.7 --magnitude-max 1.3 --phase-min -20.0 --phase-max 20.0 --numerator-name 'DCS' --denominator-name 'DCH' --labels 'DCS_CLEAN / DCH_CLEAN' --use-median
noise_subtraction_tf_DCS_DCH: $(IFO)1_hoft_DCS_frames.cache $(IFO)1_clean_C02_frames.cache
python plot_transfer_function.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --numerator-frame-cache $(IFO)1_clean_C02_frames.cache --denominator-frame-cache $(IFO)1_hoft_DCS_frames.cache --numerator-channel-list DCH-CLEAN_STRAIN_C02 --denominator-channel-name DCS-CALIB_STRAIN_CLEAN --magnitude-min 0.7 --magnitude-max 1.3 --phase-min -20.0 --phase-max 20.0 --numerator-name 'DCH' --denominator-name 'DCS' --labels 'DCH_CLEAN / DCS_CLEAN' --use-median
filters:
if [ -d Filters/$(OBSRUN)/GDSFilters ]; then \
......
......@@ -144,9 +144,16 @@ channel_list.append("CAL-CS_TDEP_PCALY_LINE1_UNCERTAINTY")
channel_list.append("CAL-CS_TDEP_PCALY_LINE2_UNCERTAINTY")
channel_list.append("CAL-CS_TDEP_DARM_LINE1_UNCERTAINTY")
channel_list.append("CAL-CS_LINE_SUM_DQ")
channel_list.append("SUS-ETMY_L1_CAL_LINE_OUT_DQ")
channel_list.append("SUS-ETMY_L2_CAL_LINE_OUT_DQ")
channel_list.append("SUS-ETMY_L3_CAL_LINE_OUT_DQ")
channel_list.append("CAL-PCALY_RX_PD_OUT_DQ")
channel_list.append("CAL-PCALY_TX_PD_OUT_DQ")
channel_list.append("SUS-ETMY_L1_CAL_LINE_FREQ")
channel_list.append("SUS-ETMY_L2_CAL_LINE_FREQ")
channel_list.append("SUS-ETMY_L3_CAL_LINE_FREQ")
channel_list.append("CAL-CS_TDEP_CAVITY_POLE_PCAL_FREQ")
channel_list.append("CAL-CS_TDEP_D2N_SPRING_PCAL_FREQ")
channel_list.append("PEM-EY_MAINSMON_EBAY_1_DQ")
channel_list.append("CAL-CS_TDEP_REF_INVA_CLGRATIO_TST_REAL")
channel_list.append("CAL-CS_TDEP_REF_INVA_CLGRATIO_TST_IMAG")
......@@ -235,11 +242,10 @@ channel_list.append("CAL-CS_TDEP_PCAL_LINE4_REF_D_REAL")
channel_list.append("CAL-CS_TDEP_PCAL_LINE4_UNCERTAINTY")
channel_list.append("CAL-DARM_CTRL_DBL_DQ")
channel_list.append("CAL-DARM_ERR_DBL_DQ")
channel_list.append("SUS-ETMY_L1_CAL_LINE_OUT_DQ")
channel_list.append("SUS-ETMY_L2_CAL_LINE_OUT_DQ")
channel_list.append("CAL-CS_TDEP_SUS_LINE2_UNCERTAINTY")
channel_list.append("CAL-CS_TDEP_SUS_LINE3_UNCERTAINTY")
channel_list.append("GRD-ISC_LOCK_OK")
channel_list.append("GRD-ISC_LOCK_ERROR")
temp_list = channel_list
channel_list = []
for chan in temp_list:
......@@ -267,16 +273,13 @@ def frame_manipulator(pipeline, name):
demux = pipeparts.mkframecppchanneldemux(pipeline, src, do_file_checksum = False, skip_bad_files = True, channel_list = map("%s:%s".__mod__, ifo_channel_list))
# Make a muxer to collect the channels we need
mux = pipeparts.mkframecppchannelmux(pipeline, None)
mux.set_property("frame-duration", frame_length)
mux.set_property("frames-per-file", frames_per_file)
mux.set_property("compression-scheme", 6)
mux.set_property("compression-level", 3)
channelmux_input_dict = {}
for key, chan in zip(channel_list, channel_list):
head_dict[key] = calibration_parts.hook_up(pipeline, demux, chan, ifo, 1.0)
head_dict[key] = pipeparts.mkprogressreport(pipeline, head_dict[key], "before muxer %s" % key)
calibration_parts.mkqueue(pipeline, head_dict[key]).get_static_pad("src").link(mux.get_request_pad("%s:%s" % (ifo, chan)))
channelmux_input_dict["%s:%s" % (ifo, chan)] = calibration_parts.mkqueue(pipeline, head_dict[key])
mux = pipeparts.mkframecppchannelmux(pipeline, channelmux_input_dict, frame_duration = frame_length, frames_per_file = frames_per_file, compression_scheme = 6, compression_level = 3)
mux = pipeparts.mkprogressreport(pipeline, mux, "end")
pipeparts.mkframecppfilesink(pipeline, mux, frame_type = frame_type, path = output_path, instrument = ifo)
......
......@@ -111,20 +111,33 @@ InputConfigs = ConfigSectionMap("InputConfigurations")
# Search the directory tree for files with names matching the one we want.
filters_name = InputConfigs["filtersfilename"]
filters_paths = []
print "\nSearching for %s ..." % filters_name
# Check the user's home directory
for dirpath, dirs, files in os.walk(os.environ['HOME']):
if filters_name in files:
# We prefer filters that came directly from a GDSFilters directory of the calibration SVN
if dirpath.count("GDSFilters") > 0:
filters_paths.insert(0, os.path.join(dirpath, filters_name))
else:
filters_paths.append(os.path.join(dirpath, filters_name))
if not len(filters_paths):
raise ValueError("Cannot find filters file %s in home directory %s or in /ligo/svncommon/CalSVN/aligocalibration/trunk/Runs/*/GDSFilters", (filters_name, os.environ['HOME']))
print "Loading calibration filters from %s\n" % filters_paths[0]
filters = numpy.load(filters_paths[0])
if filters_name.count('/') > 0:
# Then the path to the filters file was given
filters = numpy.load(filters_name)
else:
# We need to search for the filters file
filters_paths = []
print "\nSearching for %s ..." % filters_name
# Check the user's home directory
for dirpath, dirs, files in os.walk(os.environ['HOME']):
if filters_name in files:
# We prefer filters that came directly from a GDSFilters directory of the calibration SVN
if dirpath.count("GDSFilters") > 0:
filters_paths.insert(0, os.path.join(dirpath, filters_name))
else:
filters_paths.append(os.path.join(dirpath, filters_name))
# Check if there is a checkout of the entire calibration SVN
for dirpath, dirs, files in os.walk('/ligo/svncommon/CalSVN/aligocalibration/trunk/Runs/'):
if filters_name in files:
# We prefer filters that came directly from a GDSFilters directory of the calibration SVN
if dirpath.count("GDSFilters") > 0:
filters_paths.insert(0, os.path.join(dirpath, filters_name))
else:
filters_paths.append(os.path.join(dirpath, filters_name))
if not len(filters_paths):
raise ValueError("Cannot find filters file %s in home directory %s or in /ligo/svncommon/CalSVN/aligocalibration/trunk/Runs/*/GDSFilters", (filters_name, os.environ['HOME']))
print "Loading calibration filters from %s\n" % filters_paths[0]
filters = numpy.load(filters_paths[0])
ifo = options.ifo
......
......@@ -149,20 +149,33 @@ if options.config_file is not None:
# Search the directory tree for files with names matching the one we want.
filters_name = InputConfigs["filtersfilename"]
filters_paths = []
print "\nSearching for %s ..." % filters_name
# Check the user's home directory
for dirpath, dirs, files in os.walk(os.environ['HOME']):
if filters_name in files:
# We prefer filters that came directly from a GDSFilters directory of the calibration SVN
if dirpath.count("GDSFilters") > 0:
filters_paths.insert(0, os.path.join(dirpath, filters_name))
else:
filters_paths.append(os.path.join(dirpath, filters_name))
if not len(filters_paths):
raise ValueError("Cannot find filters file %s in home directory %s or in /ligo/svncommon/CalSVN/aligocalibration/trunk/Runs/*/GDSFilters", (filters_name, os.environ['HOME']))
print "Loading calibration filters from %s\n" % filters_paths[0]
filters = numpy.load(filters_paths[0])
if filters_name.count('/') > 0:
# Then the path to the filters file was given
filters = numpy.load(filters_name)
else:
# We need to search for the filters file
filters_paths = []
print "\nSearching for %s ..." % filters_name
# Check the user's home directory
for dirpath, dirs, files in os.walk(os.environ['HOME']):
if filters_name in files:
# We prefer filters that came directly from a GDSFilters directory of the calibration SVN
if dirpath.count("GDSFilters") > 0:
filters_paths.insert(0, os.path.join(dirpath, filters_name))
else:
filters_paths.append(os.path.join(dirpath, filters_name))
# Check if there is a checkout of the entire calibration SVN
for dirpath, dirs, files in os.walk('/ligo/svncommon/CalSVN/aligocalibration/trunk/Runs/'):
if filters_name in files:
# We prefer filters that came directly from a GDSFilters directory of the calibration SVN
if dirpath.count("GDSFilters") > 0:
filters_paths.insert(0, os.path.join(dirpath, filters_name))
else:
filters_paths.append(os.path.join(dirpath, filters_name))
if not len(filters_paths):
raise ValueError("Cannot find filters file %s in home directory %s or in /ligo/svncommon/CalSVN/aligocalibration/trunk/Runs/*/GDSFilters", (filters_name, os.environ['HOME']))
print "Loading calibration filters from %s\n" % filters_paths[0]
filters = numpy.load(filters_paths[0])
# Get the corrections for numerators and denominator, and resample if necessary
num_corr = []
......
#!/usr/bin/env python
# Copyright (C) 2016 Aaron Viets
#
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 2 of the License, or (at your
# option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
#
# =============================================================================
#
# Preamble
#
# =============================================================================
#
import random
import sys
import os
import numpy
import time
from math import pi
import resource
import datetime
import time
import matplotlib
matplotlib.rcParams['font.family'] = 'Times New Roman'
matplotlib.rcParams['font.size'] = 22
matplotlib.rcParams['legend.fontsize'] = 18
matplotlib.rcParams['mathtext.default'] = 'regular'
matplotlib.use('Agg')
import glob
import matplotlib.pyplot as plt
from optparse import OptionParser, Option
import ConfigParser
import gi
gi.require_version('Gst', '1.0')
from gi.repository import GObject, Gst
GObject.threads_init()
Gst.init(None)
import lal
from lal import LIGOTimeGPS
from gstlal import pipeparts
from gstlal import calibration_parts
from gstlal import simplehandler
from gstlal import datasource
from glue.ligolw import ligolw
from glue.ligolw import array
from glue.ligolw import param
from glue.ligolw.utils import segments as ligolw_segments
array.use_in(ligolw.LIGOLWContentHandler)
param.use_in(ligolw.LIGOLWContentHandler)
from glue.ligolw import utils
from ligo import segments
from gstlal import test_common
#
# =============================================================================
#
# Pipelines
#
# =============================================================================
#
def fill_silence_test_01(pipeline, name, line_sep = 0.5):
#
# This test is intended to help get rid of error messages associated with adding two streams that have different start times
#
rate = 16 # Hz
buffer_length = 1.0 # seconds
test_duration = 300.0 # seconds
filter_latency = 1.0
#
# build pipeline
#
head = test_common.test_src(pipeline, rate = rate, test_duration = test_duration, wave = 5)
head = pipeparts.mktee(pipeline, head)
smooth = calibration_parts.mkcomplexfirbank(pipeline, head, latency = int((rate * 40 - 1) * filter_latency + 0.5), fir_matrix = [numpy.ones(rate * 40)], time_domain = True)
smooth = calibration_parts.mkcomplexfirbank(pipeline, smooth, latency = 23, fir_matrix = [numpy.ones(45)], time_domain = True)
#smooth = pipeparts.mkgeneric(pipeline, head, "lal_smoothkappas", array_size = rate * 128, avg_array_size = rate * 10, filter_latency = 1)
#smooth = pipeparts.mkgeneric(pipeline, smooth, "splitcounter", name = "smooth")
#head = pipeparts.mkgeneric(pipeline, head, "splitcounter", name = "unsmooth")
#channelmux_input_dict = {}
#channelmux_input_dict["unsmooth"] = calibration_parts.mkqueue(pipeline, head)
#channelmux_input_dict["smooth"] = calibration_parts.mkqueue(pipeline, smooth)
#mux = pipeparts.mkframecppchannelmux(pipeline, channelmux_input_dict, frame_duration = 64, frames_per_file = 1, compression_scheme = 6, compression_level = 3)
head = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, head, smooth))
#mux = pipeparts.mkgeneric(pipeline, mux, "splitcounter", name = "sum")
#head = calibration_parts.mkgate(pipeline, smooth, head, 0)
pipeparts.mknxydumpsink(pipeline, head, "%s_out.txt" % name)
#pipeparts.mkframecppfilesink(pipeline, mux, frame_type = "H1DCS", instrument = "H1")
#
# done
#
return pipeline
#
# =============================================================================
#
# Main
#
# =============================================================================
#
test_common.build_and_run(fill_silence_test_01, "fill_silence_test_01")
#!/usr/bin/env python
# Copyright (C) 2018 Aaron Viets
#
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 2 of the License, or (at your
# option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
#
# =============================================================================
#
# Preamble
#
# =============================================================================
#
import numpy
import sys
from gstlal import pipeparts
from gstlal import calibration_parts
import test_common
#
# =============================================================================
#
# Pipelines
#
# =============================================================================
#
def lal_property_test_01(pipeline, name):
#
# This test makes a stream changing between 2's and 8's.
# It should square the 2's and take the 8's to the 8th power.
#
rate = 512 # Hz
width = 64
buffer_length = 1.0 # seconds
test_duration = 100.0 # seconds
gap_frequency = 0.1 # Hz
gap_threshold = 0.0
control_dump_filename = "control_property_test_01.dump"
bad_data_intervals2 = [0.0, 1e35]
bad_data_intervals = [-1e35, 1e-35]
head = test_common.test_src(pipeline, buffer_length = buffer_length, rate = rate, width = width, channels = 1, test_duration = test_duration, wave = 0, freq = 0.1, volume = 1)
head = pipeparts.mkgeneric(pipeline, head, "lal_insertgap", bad_data_intervals = bad_data_intervals, insert_gap = False, fill_discont = True, replace_value = 2.0)
head = pipeparts.mkgeneric(pipeline, head, "lal_insertgap", bad_data_intervals = bad_data_intervals2, insert_gap = False, fill_discont = True, replace_value = 8.0)
head = pipeparts.mktee(pipeline, head)
pipeparts.mknxydumpsink(pipeline, head, "%s_in.dump" % name)
lal_prop_exponent = pipeparts.mkgeneric(pipeline, head, "lal_property", update_when_change = True)
head = calibration_parts.mkpow(pipeline, head, exponent = 0.0)
lal_prop_exponent.connect("notify::current-average", calibration_parts.update_property_simple, head, "current_average", "exponent")
pipeparts.mknxydumpsink(pipeline, head, "%s_out.dump" % name)
#
# done
#
return pipeline
#
# =============================================================================
#
# Main
#
# =============================================================================
#
test_common.build_and_run(lal_property_test_01, "lal_property_test_01")
......@@ -691,8 +691,10 @@ for output_file_number, (svd_bank_url_dict, output_url, ranking_stat_output_url,
# assume all instruments have the same templates, just extract them
# from one of the instruments at random
sngl_inspiral_table = banks.values()[0][0].sngl_inspiral_table.copy()
horizon_factors = {}
for bank in banks.values()[0]:
sngl_inspiral_table.extend(bank.sngl_inspiral_table)
horizon_factors.update(bank.horizon_factors)
@bottle.route("/template_bank.xml.gz")
def get_template_bank_xml(sngl_inspiral_table = sngl_inspiral_table):
xmldoc = ligolw.Document()
......@@ -704,6 +706,7 @@ for output_file_number, (svd_bank_url_dict, output_url, ranking_stat_output_url,
output.close()
return outstr
template_ids = frozenset(row.template_id for row in sngl_inspiral_table)
assert set(horizon_factors) == template_ids, "horizon factors are not assigned for each template id"
assert len(template_ids) == len(sngl_inspiral_table), "template IDs are not unique within the template bank"
@bottle.route("/template_ids.txt")
def get_template_ids(template_ids = sorted(template_ids)):
......@@ -783,9 +786,10 @@ for output_file_number, (svd_bank_url_dict, output_url, ranking_stat_ou