Commit 6ee2014b authored by Aaron Viets's avatar Aaron Viets

lal_sensingtdcfs: New element in gstlal-calibration to solve for the...

lal_sensingtdcfs:  New element in gstlal-calibration to solve for the time-dependent correction factors of the sensing function
parent 2386bdfb
Pipeline #74387 failed with stages
in 1 minute and 6 seconds
......@@ -20,7 +20,8 @@ lib@GSTPLUGINPREFIX@gstlalcalibration_la_SOURCES = \
gstlal_dqtukey.c gstlal_dqtukey.h \
gstlal_property.c gstlal_property.h \
gstlal_typecast.c gstlal_typecast.h \
gstlal_matrixsolver.c gstlal_matrixsolver.h
gstlal_matrixsolver.c gstlal_matrixsolver.h \
gstlal_sensingtdcfs.c gstlal_sensingtdcfs.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)
......@@ -61,7 +61,7 @@ struct _GSTLALMatrixSolver {
GSTLAL_MATRIXSOLVER_F32 = 0,
GSTLAL_MATRIXSOLVER_F64,
GSTLAL_MATRIXSOLVER_Z64,
GSTLAL_MATRIXSOLVER_Z128,
GSTLAL_MATRIXSOLVER_Z128
} data_type;
/* timestamp book-keeping */
......
This diff is collapsed.
/*
* Copyright (C) 2019 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.,
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*/
#ifndef __GSTLAL_SENSINGTDCFS_H__
#define __GSTLAL_SENSINGTDCFS_H__
#include <glib.h>
#include <gst/gst.h>
#include <gst/base/gstbasetransform.h>
G_BEGIN_DECLS
#define GSTLAL_SENSINGTDCFS_TYPE \
(gstlal_sensingtdcfs_get_type())
#define GSTLAL_SENSINGTDCFS(obj) \
(G_TYPE_CHECK_INSTANCE_CAST((obj), GSTLAL_SENSINGTDCFS_TYPE, GSTLALSensingTDCFs))
#define GSTLAL_SENSINGTDCFS_CLASS(klass) \
(G_TYPE_CHECK_CLASS_CAST((klass), GSTLAL_SENSINGTDCFS_TYPE, GSTLALSensingTDCFsClass))
#define GST_IS_GSTLAL_SENSINGTDCFS(obj) \
(G_TYPE_CHECK_INSTANCE_TYPE((obj), GSTLAL_SENSINGTDCFS_TYPE))
#define GST_IS_GSTLAL_SENSINGTDCFS_CLASS(klass) \
(G_TYPE_CHECK_CLASS_TYPE((klass), GSTLAL_SENSINGTDCFS_TYPE))
typedef struct _GSTLALSensingTDCFs GSTLALSensingTDCFs;
typedef struct _GSTLALSensingTDCFsClass GSTLALSensingTDCFsClass;
/**
* GSTLALSensingTDCFs:
*/
struct _GSTLALSensingTDCFs {
GstBaseTransform element;
/* stream info */
gint rate;
gint channels_in;
gint channels_out;
gint unit_size_in;
gint unit_size_out;
enum gstlal_sensingtdcfs_data_type {
GSTLAL_SENSINGTDCFS_FLOAT = 0,
GSTLAL_SENSINGTDCFS_DOUBLE
} data_type;
/* timestamp book-keeping */
GstClockTime t0;
guint64 offset0;
guint64 next_in_offset;
guint64 next_out_offset;
gboolean need_discont;
/* properties */
gint sensing_model;
double freq1;
double freq2;
double freq4;
};
/**
* GSTLALSensingTDCFsClass:
* @parent_class: the parent class
*/
struct _GSTLALSensingTDCFsClass {
GstBaseTransformClass parent_class;
};
GType gstlal_sensingtdcfs_get_type(void);
G_END_DECLS
#endif /* __GSTLAL_SENSINGTDCFS_H__ */
......@@ -72,6 +72,7 @@
#include <gstlal_property.h>
#include <gstlal_typecast.h>
#include <gstlal_matrixsolver.h>
#include <gstlal_sensingtdcfs.h>
/*
......@@ -108,6 +109,7 @@ static gboolean plugin_init(GstPlugin *plugin)
{"lal_property", GSTLAL_PROPERTY_TYPE},
{"lal_typecast", GSTLAL_TYPECAST_TYPE},
{"lal_matrixsolver", GSTLAL_MATRIXSOLVER_TYPE},
{"lal_sensingtdcfs", GSTLAL_SENSINGTDCFS_TYPE},
{NULL, 0},
};
......
......@@ -1062,7 +1062,7 @@ def compute_exact_kappas_from_filters_file(pipeline, X, freqs, EPICS, rate):
MV_matrix = mkinterleave(pipeline, MV_matrix)
MV_matrix = pipeparts.mkcapsfilter(pipeline, MV_matrix, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=%d" % (rate, (2 * num_stages) * (2 * num_stages + 1)))
kappas = pipeparts.mkgeneric(pipeline, MV_matrix, "lal_matrixsolver")
kappas = list(mkdeinterleave(pipeline, kappas, 2 * num_stages))
kappas = list(mkdeinterleave(pipeline, pipeparts.mkcapsfilter(pipeline, kappas, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=%d" % (rate, 2 * num_stages)), 2 * num_stages))
for i in range(len(kappas)):
kappas[i] = pipeparts.mkcapsfilter(pipeline, kappas[i], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)
if i >= len(kappas) / 2:
......
......@@ -8,11 +8,11 @@ IFO = H
# determines where to look for filters files (e.g., O1, O2, O3, ER10, ER13, ER14, PreER10, PreER13, PreER14)
OBSRUN = O3
START = $(shell echo 1240617600 - 256 | bc)
START = $(shell echo 1246226240 - 256 | bc)
#1238997618
END = $(shell echo 1240617600 + 1024 + 256 | bc)
END = $(shell echo 1246226240 + 120 + 256 | bc)
#1240876818
SHMRUNTIME = 36000
SHMRUNTIME = 600
# How much time does the calibration need to settle at the start and end?
PLOT_WARMUP_TIME = 256
PLOT_COOLDOWN_TIME = 256
......@@ -25,13 +25,13 @@ DCSLINESCONFIGS = ../../config_files/O2/H1/tests/H1DCS_AllCorrections_CleaningLi
#../../config_files/O2/H1/tests/H1DCS_AllCorrections_Cleaning.ini
DCSFCCCONFIGS = ../../config_files/O2/H1/tests/H1DCS_FreqIndepAndFccCorrections_Cleaning.ini
GDSTESTCONFIGS = Filters/O3/GDSFilters/H1GDS_1239476409_testAllCorrections.ini
DCSTESTCONFIGS = Filters/O3/GDSFilters/L1DCS_C01_1239579018_test.ini
GDSSHMCONFIGS = Filters/ER14/GDSFilters/H1GDS_1234630818_latency_test.ini
DCSTESTCONFIGS = Filters/O3/GDSFilters/H1DCS_C01_1239472998_AllCorrectionsTest.ini
GDSSHMCONFIGS = Filters/O3/GDSFilters/H1GDS_lowLatency_1239476409.ini
GDSOLDCONFIGS = Filters/ER14/GDSFilters/L1GDS_1235491416_old.ini
GDSBETTERCONFIGS = Filters/ER14/GDSFilters/L1GDS_1235491416_better.ini
GDSBESTCONFIGS = Filters/ER14/GDSFilters/L1GDS_1235491416_best.ini
all: highpass_filter_ASD_GDS filters_tf_GDS
all: pcal_DCS_transfer_functions
###############################################
### These commands should change less often ###
......@@ -123,8 +123,9 @@ DCS_pcal2darm_plots:
#$(IFO)1_easy_raw_frames.cache $(IFO)1_hoft_DCS_frames.cache $(IFO)1_C00_frames.cache
python pcal2darm_timeseries.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --raw-frame-cache $(IFO)1_easy_raw_frames.cache --gstlal-frame-cache-list '$(IFO)1_C00_frames.cache,$(IFO)1_hoft_DCS_frames.cache' --config-file '$(DCSCONFIGS)' --pcal-channel-name CAL-PCALY_RX_PD_OUT_DQ --gstlal-channel-list 'GDS-CALIB_STRAIN,DCS-CALIB_STRAIN' --labels 'C00,TEST' --pcal-time-advance 0.00006103515625
pcal_DCS_transfer_functions: $(IFO)1_easy_raw_frames.cache $(IFO)1_hoft_DCS_frames.cache $(IFO)1_hoft_DCS_FCC_frames.cache $(IFO)1_hoft_DCS_TEST_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_easy_raw_frames.cache --denominator-channel-name CAL-PCALY_TX_PD_OUT_DQ --denominator-correction y_arm_pcal_corr --numerator-correction arm_length --frequency-max 400 --numerator-frame-cache-list $(IFO)1_hoft_DCS_frames.cache,$(IFO)1_hoft_DCS_FCC_frames.cache,$(IFO)1_hoft_DCS_TEST_frames.cache --numerator-channel-list DCS-CALIB_STRAIN_CLEAN,DCS-CALIB_STRAIN_CLEAN,DCS-CALIB_STRAIN_CLEAN --config-file $(DCSCONFIGS) --use-median --labels 'Only scalar corrections,Cavity pole corrections,All corrections'
pcal_DCS_transfer_functions:
#$(IFO)1_hoft_DCS_frames.cache $(IFO)1_hoft_DCS_TEST_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_easy_raw_frames.cache --denominator-channel-name CAL-PCALY_RX_PD_OUT_DQ --denominator-correction y_arm_pcal_corr --numerator-correction arm_length --frequency-max 400 --numerator-frame-cache-list $(IFO)1_hoft_DCS_frames.cache,$(IFO)1_hoft_DCS_TEST_frames.cache --numerator-channel-list DCS-CALIB_STRAIN,DCS-CALIB_STRAIN --config-file $(DCSCONFIGS) --use-median --labels 'Without SRC Corrections,With SRC Corrections'
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'
......
......@@ -69,7 +69,7 @@ parser.add_option("--gps-end-time", metavar = "seconds", type = int, help = "GPS
parser.add_option("--ifo", metavar = "name", type = str, help = "Name of the interferometer (IFO), e.g., H1, L1")
parser.add_option("--denominator-frame-cache", metavar = "name", type = str, help = "Name of frame cache file that contains denominator of transfer functions")
parser.add_option("--denominator-channel-name", metavar = "name", type = str, default = None, help = "Channel-name of denominator")
parser.add_option("--denominator-name", metavar = "name", type = str, default = 'Pcal(f)', help = "Name of denominator in plot title, in latex math mode")
parser.add_option("--denominator-name", metavar = "name", type = str, default = '{\\rm Pcal}(f)', help = "Name of denominator in plot title, in latex math mode")
parser.add_option("--denominator-correction", metavar = "name", type = str, default = None, help = "Name of filters-file parameter needed to apply a correction to the denominator")
parser.add_option("--numerator-frame-cache-list", metavar = "name", type = str, help = "Comma-separated list of frame cache files that contain numerators of transfer functions")
parser.add_option("--numerator-channel-list", metavar = "list", type = str, default = None, help = "Comma-separated list of channel-names of numerators")
......@@ -177,7 +177,7 @@ if options.config_file is not None:
num_corr = []
denom_corr = []
if options.numerator_correction is not None:
if len(filters[options.numerator_correction]) > 1:
if numpy.size(filters[options.numerator_correction]) > 1:
corr = filters[options.numerator_correction]
# Check the frequency spacing of the correction
corr_df = corr[0][1] - corr[0][0]
......@@ -199,7 +199,7 @@ if options.numerator_correction is not None:
num_corr.append(corr[1][before_idx] + 1j * corr[2][before_idx])
if options.denominator_correction is not None:
if len(filters[options.denominator_correction]) > 1:
if numpy.size(filters[options.denominator_correction]) > 1:
corr = filters[options.denominator_correction]
# Check the frequency spacing of the correction
corr_df = corr[0][1] - corr[0][0]
......@@ -243,7 +243,7 @@ def plot_transfer_function(pipeline, name):
denominator = calibration_parts.caps_and_progress(pipeline, denominator, "audio/x-raw,format=F64LE", "denominator")
denominator = calibration_parts.mkresample(pipeline, denominator, 5, False, int(sample_rate))
if options.denominator_correction is not None:
if len(filters[options.denominator_correction]) == 1:
if numpy.size(filters[options.denominator_correction]) == 1:
denominator = pipeparts.mkaudioamplify(pipeline, denominator, float(filters[options.denominator_correction]))
denominator = pipeparts.mktee(pipeline, denominator)
......@@ -251,11 +251,11 @@ def plot_transfer_function(pipeline, name):
for i in range(0, len(labels)):
numerator = pipeparts.mklalcachesrc(pipeline, location = numerator_frame_cache_list[i], cache_dsc_regex = ifo)
numerator = pipeparts.mkframecppchanneldemux(pipeline, numerator, do_file_checksum = False, skip_bad_files = True, channel_list = map("%s:%s".__mod__, channel_list))
numerator = calibration_parts.hook_up(pipeline, numerator, numerator_channel_list[i], ifo, 1.0)
numerator = calibration_parts.hook_up(pipeline, numerator, numerator_channel_list[i], ifo, 1.0, element_name_suffix = "%d" % i)
numerator = calibration_parts.caps_and_progress(pipeline, numerator, "audio/x-raw,format=F64LE", labels[i])
numerator = calibration_parts.mkresample(pipeline, numerator, 5, False, int(sample_rate))
if options.numerator_correction is not None:
if len(filters[options.numerator_correction] == 1):
if numpy.size(filters[options.numerator_correction] == 1):
numerator = pipeparts.mkaudioamplify(pipeline, numerator, float(filters[options.numerator_correction]))
# Interleave the channels to make one stream
channels = calibration_parts.mkinterleave(pipeline, [numerator, denominator])
......@@ -294,7 +294,7 @@ if options.poles is not None:
for i in range(0, len(real_poles) / 2):
poles.append(float(real_poles[2 * i]) + 1j * float(real_poles[2 * i + 1]))
colors = ['limegreen', 'g', 'y', 'c', 'm', 'b'] # Hopefully the user will not want to plot more than six datasets on one plot.
colors = ['blue', 'limegreen', 'y', 'c', 'm', 'b'] # Hopefully the user will not want to plot more than six datasets on one plot.
for i in range(0, len(labels)):
# Remove unwanted lines from file, and re-format wanted lines
f = open('%s_%s_over_%s_%d-%d.txt' % (ifo, labels[i].replace(' ', '_').replace('/', 'over'), options.denominator_channel_name, options.gps_start_time, data_duration),"r")
......@@ -330,13 +330,13 @@ for i in range(0, len(labels)):
if i == 0:
plt.figure(figsize = (10, 10))
plt.subplot(211)
plt.plot(frequency, magnitude, colors[i % 6], linewidth = 0.75, label = labels[i].replace('_', '\_'))
#leg = plt.legend(fancybox = True)
#leg.get_frame().set_alpha(0.8)
plt.plot(frequency, magnitude, colors[i % 6], linewidth = 0.75, label = r'${\rm %s}$' % labels[i].replace('_', '\_').replace(' ', '\ '))
leg = plt.legend(fancybox = True)
leg.get_frame().set_alpha(0.8)
plt.gca().set_xscale(freq_scale)
plt.gca().set_yscale(mag_scale)
if i == 0:
#plt.title(r'%s $%s$ / $%s$' % (ifo, options.numerator_name.replace('_', '\_'), options.denominator_name.replace('_', '\_')))
plt.title(r'${\rm %s} \ %s \ / \ %s$' % (ifo, options.numerator_name, options.denominator_name))
plt.ylabel(r'${\rm Magnitude}$')
plt.xlim(options.frequency_min, options.frequency_max)
plt.ylim(options.magnitude_min, options.magnitude_max)
......@@ -350,6 +350,6 @@ for i in range(0, len(labels)):
plt.xlim(options.frequency_min, options.frequency_max)
plt.ylim(options.phase_min, options.phase_max)
plt.grid(True, which = "both", linestyle = ':', linewidth = 0.3, color = 'black')
plt.savefig('%s_transfer_functions_%s_%s%s_%d-%d.png' % (ifo, options.numerator_name, options.denominator_name, options.filename_suffix, options.gps_start_time, data_duration))
plt.savefig('%s_transfer_functions_%s_%s%s_%d-%d.pdf' % (ifo, options.numerator_name, options.denominator_name, options.filename_suffix, options.gps_start_time, data_duration))
plt.savefig('%s_%s_over_%s_%d-%d.pdf' % (ifo, numerator_channel_list[-1], options.denominator_channel_name, options.gps_start_time, data_duration))
plt.savefig('%s_%s_over_%s_%d-%d.png' % (ifo, numerator_channel_list[-1], options.denominator_channel_name, options.gps_start_time, data_duration))
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment