diff --git a/gstlal-calibration/bin/gstlal_compute_strain b/gstlal-calibration/bin/gstlal_compute_strain
index 47301ef9e4493484e2c39bf622d3c57af8d0887b..e4ab5900f0aece48809c73119fe33fdd33e60e60 100755
--- a/gstlal-calibration/bin/gstlal_compute_strain
+++ b/gstlal-calibration/bin/gstlal_compute_strain
@@ -253,6 +253,7 @@ powerlines_tf_averaging_time = float(DataCleaningConfigs["powerlinestfaveragingt
powerlines_freq_var = float(DataCleaningConfigs["powerlinesfreqvar"])
witness_channel_fft_time = float(DataCleaningConfigs["witnesschannelffttime"])
num_witness_ffts = int(DataCleaningConfigs["numwitnessffts"])
+min_witness_ffts = int(DataCleaningConfigs["minwitnessffts"])
witness_fir_length = float(DataCleaningConfigs["witnessfirlength"])
witness_frequency_resolution = float(DataCleaningConfigs["witnessfrequencyresolution"])
witness_tf_update_time = float(DataCleaningConfigs["witnesstfupdatetime"])
@@ -326,10 +327,11 @@ tdcf_default_to_median = Config.getboolean("TDCFConfigurations", "tdcfdefaulttom
compute_calib_statevector = Config.getboolean("CalibrationConfigurations", "computecalibstatevector")
# Boolean for computing factors from filters file
factors_from_filters_file = Config.getboolean("TDCFConfigurations", "factorsfromfiltersfile")
-# Boolean for removing calibration lines, power lines, DC component
+# Boolean for removing calibration lines, power lines, DC component, using median in witness transfer functions
remove_cal_lines = Config.getboolean("DataCleaningConfigurations", "removecallines")
remove_power_lines = Config.getboolean("DataCleaningConfigurations", "removepowerlines")
remove_dc = Config.getboolean("DataCleaningConfigurations", "removedc")
+witness_tf_use_median = Config.getboolean("DataCleaningConfigurations", "witnesstfusemedian")
# If td is true we will perform filtering in the time domain (direct convolution) in all FIR filtering routines below
td = not Config.getboolean("PipelineConfigurations", "frequencydomainfiltering")
# Boolean for dewhitening
@@ -746,8 +748,7 @@ if compute_kappac or compute_fcc or compute_kappatst or compute_kappapu or compu
if dewhitening:
# dewhiten DARM_ERR at the DARM actuation line frequency
derr_at_darm_act_freq = calibration_parts.complex_audioamplify(pipeline, derr_at_darm_act_freq, derr_dewhiten_at_darm_act_freq_real, derr_dewhiten_at_darm_act_freq_imag)
- if compute_kappapu or compute_kappac or compute_fcc or compute_srcq or compute_fs:
- derr_at_darm_act_freq = pipeparts.mktee(pipeline, derr_at_darm_act_freq)
+ derr_at_darm_act_freq = pipeparts.mktee(pipeline, derr_at_darm_act_freq)
# demodulate the TST excitation channel at the ESD actuation line frequency
tstexc_at_esd_act_freq = calibration_parts.demodulate(pipeline, tstexc, esd_act_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor)
@@ -940,14 +941,20 @@ if compute_fs or compute_srcq:
Xi_imag_ok_var = float(fs_var / (srcQ_default * src_pcal_line_freq))
# demodulate PCAL channel and apply the PCAL correction factor at SRC detuning line frequency
- pcal_at_src_freq = calibration_parts.demodulate(pipeline, pcaltee, src_pcal_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor, prefactor_real = pcal_corr_at_src_freq_real, prefactor_imag = pcal_corr_at_src_freq_imag)
+ if src_pcal_line_freq == darm_act_line_freq:
+ pcal_at_src_freq = pcal_at_darm_act_freq
+ else:
+ pcal_at_src_freq = calibration_parts.demodulate(pipeline, pcaltee, src_pcal_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor, prefactor_real = pcal_corr_at_src_freq_real, prefactor_imag = pcal_corr_at_src_freq_imag)
+ pcal_at_src_freq = pipeparts.mktee(pipeline, pcal_at_src_freq)
if "pcal4" in pcal_line_removal_dict:
# This will save having to demodulate it again
- pcal_at_src_freq = pipeparts.mktee(pipeline, pcal_at_src_freq)
pcal_line_removal_dict["pcal4"][3] = pcal_at_src_freq
# demodulate DARM_ERR at SRC detuning line frequency
- derr_at_src_freq = calibration_parts.demodulate(pipeline, derrtee, src_pcal_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor)
+ if src_pcal_line_freq == darm_act_line_freq:
+ derr_at_src_freq = derr_at_darm_act_freq
+ else:
+ derr_at_src_freq = calibration_parts.demodulate(pipeline, derrtee, src_pcal_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor)
# Compute the factor Xi which will be used for the f_s and src_Q calculations
# \kappa_tst and \kappa_pu need to be evaluated at the SRC pcal line frequency
@@ -1854,7 +1861,7 @@ if witness_channel_list is not None:
witnesses.append(calibration_parts.caps_and_progress(pipeline, head_dict[key], "audio/x-raw, format=F64LE, channels=1, channel-mask=(bitmask)0x0", key))
if len(witnesses) != len(witness_channel_list[i]):
print "WARNING: Not all requested witness channels are being used"
- clean_strain = calibration_parts.clean_data(pipeline, clean_strain, hoft_sr, witnesses, witness_rates[i], witness_fft_samples, witness_fft_overlap, num_witness_ffts, witness_tf_update_samples, witness_fir_samples, witness_frequency_resolution, witness_filter_taper_length, notch_frequencies = witness_notch_frequencies[i], obsready = obsreadytee, chop_time = witness_chop_time, wait_time = witness_wait_time, filename = "transfer_functions_%d.txt" % i)
+ clean_strain = calibration_parts.clean_data(pipeline, clean_strain, hoft_sr, witnesses, witness_rates[i], witness_fft_samples, witness_fft_overlap, num_witness_ffts, min_witness_ffts, witness_tf_update_samples, witness_fir_samples, witness_frequency_resolution, witness_filter_taper_length, use_median = witness_tf_use_median, notch_frequencies = witness_notch_frequencies[i], obsready = obsreadytee, chop_time = witness_chop_time, wait_time = witness_wait_time, filename = "transfer_functions_%d.txt" % i)
witness_chop_time += witness_chop_increment
witness_wait_time += witness_wait_increment
diff --git a/gstlal-calibration/config_files/H1DCS_AllCorrections_Cleaning.ini b/gstlal-calibration/config_files/H1DCS_AllCorrections_Cleaning.ini
new file mode 100644
index 0000000000000000000000000000000000000000..d992703a9eb32981123f58835b3154f77eac8225
--- /dev/null
+++ b/gstlal-calibration/config_files/H1DCS_AllCorrections_Cleaning.ini
@@ -0,0 +1,316 @@
+[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: Filters/O2filters/GDSFilters/H1DCS_newsrcline_1173225472.npz
+# Data source should be set to frames or lvshm
+DataSource: frames
+FileChecksum: No
+# Right now, SkipBadFiles needs to be off when reading from frames
+SkipBadFiles: No
+############################################
+# If reading from frames use these options #
+############################################
+FrameCache: easy_raw_frames.cache
+###################################################
+# If reading from shared memory use these options #
+###################################################
+SHMPartition: LHO_Online
+# Assumed duration of input frames in seconds
+InputFrameDuration: 1
+
+[OutputConfigurations]
+CompressionScheme: 6
+CompressionLevel: 3
+ChanPrefix: DCS-
+# Set to "None" if you do not want a channel suffix
+ChanSuffix: None
+# Data sink should be set to frames or lvshm
+DataSink: frames
+#################################################
+# If writing to shared memory use these options #
+#################################################
+OutputSHMPartition: hoft_test
+BufferMode: 2
+# Use this to approximate the frame size (in bytes) when writing to shared memory
+FrameSize: 405338
+NumBuffers: 10
+###############################################
+# If writing to frame files use these options #
+###############################################
+FrameType: H1DCS_TEST
+OutputPath: .
+
+[CalibrationConfigurations]
+IFO: H1
+# Set calibration mode to Full or Partial
+CalibrationMode: Full
+ComputeCalibStateVector: Yes
+
+[DebuggingConfigurations]
+# If you want to write a pipeline graph, provide the graph name. Otherwise, set name equal to None
+PipelineGraphFilename: gstlal_compute_strain
+Verbose: Yes
+
+[TDCFConfigurations]
+#########################################################
+# Options related to time dependent correction factors #
+#########################################################
+ComputeKappaTST: Yes
+ApplyKappaTST: No
+# Set this to have the \kappa_tst factors filter the actuation chain with an adaptive filter that corrects for both magnitude and phase errors.
+ApplyComplexKappaTST: Yes
+ComputeKappaPU: Yes
+ApplyKappaPU: No
+# Set this to have the \kappa_pu factors the actuation chain with an adaptive filter that corrects for both magnitude and phase errors
+ApplyComplexKappaPU: Yes
+ComputeKappaP: No
+ApplyKappaP: No
+# Set this to have the \kappa_p factors the actuation chain with an adaptive filter that corrects for both magnitude and phase errors.
+ApplyComplexKappaP: No
+ComputeKappaU: No
+ApplyKappaU: No
+# Set this to have the \kappa_u factors the actuation chain with an adaptive filter that corrects for both magnitude and phase errors.
+ApplyComplexKappaU: No
+ComputeKappaC: Yes
+ApplyKappaC: Yes
+ComputeFcc: Yes
+ApplyFcc: Yes
+ComputeSRCQ: Yes
+ApplySRCQ: Yes
+ComputeFs: Yes
+ApplyFs: Yes
+###########################################
+# Options related to the coherence gating #
+###########################################
+UseCoherence: Yes
+CoherenceUncThreshold: 0.004
+# Amount of time used in front-end to compute coherence
+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
+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
+##################################################
+# Options related to updating cavity pole filter #
+##################################################
+# Duration of the Fcc filter in the time domain in seconds
+FccFilterDuration: 0.01
+# Number of seconds to average Fcc values before creating a new Fcc filter
+FccAveragingTime: 60
+# Number of samples to be used when tapering old filter and ramping in new filter
+FccFilterTaperLength: 32768
+############################
+# Nominal values for TDCFs #
+############################
+ExpectedKappaTSTReal: 1.0
+ExpectedKappaTSTImag: 0.0
+ExpectedKappaPUReal: 1.0
+ExpectedKappaPUImag: 0.0
+ExpectedKappaPReal: 1.0
+ExpectedKappaPImag: 0.0
+ExpectedKappaUReal: 1.0
+ExpectedKappaUImag: 0.0
+ExpectedKappaC: 1.0
+ExpectedFcc: 360.0
+ExpectedFs: 6.91
+ExpectedSRCQ: 21.739
+################################
+# Acceptable variance in TDCFs #
+################################
+KappaTSTRealVar: 0.2
+KappaTSTImagVar: 0.2
+KappaPURealVar: 0.2
+KappaPUImagVar: 0.2
+KappaPRealVar: 0.2
+KappaPImagVar: 0.2
+KappaURealVar: 0.2
+KappaUImagVar: 0.2
+KappaCVar: 0.2
+FccVar: 50.0
+FsVar: 5.0
+SRCQInvMin: 0.0
+SRCQInvMax: 0.5
+#######################
+# EPICS records input #
+#######################
+# Set to Yes if EPICS records for TDCF computations should be read from filters file. If set to No, they will be read from frames
+FactorsFromFiltersFile: Yes
+##############################################################################################
+# Updating Sensing and Actuation filters with all frequency-dependent corrections parameters #
+##############################################################################################
+# Length of time (in seconds) between when inverse-sensing FIR filter is updated
+SensingFilterUpdateTime: 60
+# Length of time (in seconds) over which the smoothed time-dependent parameters of the sensing function are averaged before updating the filter
+SensingFilterAveragingTime: 1
+# Number of samples to be used when tapering old inverse sensing filter and ramping in new filter
+SensingFilterTaperLength: 32768
+# Length of time (in seconds) between when the actuation FIR filters are updated
+ActuationFilterUpdateTime: 60
+# Length of time (in seconds) over which the smoothed time-dependent parameters of the actuation function are averaged before updating the filter
+ActuationFilterAveragingTime: 1
+# Number of samples to be used when tapering old actuation filters and ramping in new filters
+ActuationFilterTaperLength: 32768
+
+[ChannelNames]
+#############################
+# Calibration Channel Names #
+#############################
+DARMCtrlChannel: CAL-DARM_CTRL_WHITEN_OUT_DBL_DQ
+DARMErrChannel: CAL-DARM_ERR_WHITEN_OUT_DBL_DQ
+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
+##################################
+# Calibration Line Channel Names #
+##################################
+DARMExcChannel: CAL-CS_LINE_SUM_DQ
+TSTExcChannel: SUS-ETMY_L3_CAL_LINE_OUT_DQ
+PCALChannel: CAL-PCALY_RX_PD_OUT_DQ
+#######################################
+# Coherence Uncertainty Channel Names #
+#######################################
+CohUncSusLine1Channel: CAL-CS_TDEP_SUS_LINE1_UNCERTAINTY
+CohUncPcalyLine1Channel: CAL-CS_TDEP_PCALY_LINE1_UNCERTAINTY
+CohUncPcalyLine2Channel: CAL-CS_TDEP_PCALY_LINE2_UNCERTAINTY
+CohUncDARMLine1Channel: CAL-CS_TDEP_DARM_LINE1_UNCERTAINTY
+###################################
+# Noise Subtraction Channel Names #
+###################################
+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
+###############################
+# EPICS Records Channel Names #
+###############################
+EP1RealChannel: CAL-CS_TDEP_REF_INVA_CLGRATIO_TST_REAL
+EP1ImagChannel: CAL-CS_TDEP_REF_INVA_CLGRATIO_TST_IMAG
+EP2RealChannel: CAL-CS_TDEP_REF_CLGRATIO_CTRL_REAL
+EP2ImagChannel: CAL-CS_TDEP_REF_CLGRATIO_CTRL_IMAG
+EP3RealChannel: CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_INV_REAL
+EP3ImagChannel: CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_INV_IMAG
+EP4RealChannel: CAL-CS_TDEP_DARM_LINE1_REF_A_TST_REAL
+EP4ImagChannel: CAL-CS_TDEP_DARM_LINE1_REF_A_TST_IMAG
+EP5RealChannel: CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_REAL
+EP5ImagChannel: CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_IMAG
+EP6RealChannel: CAL-CS_TDEP_PCALY_LINE2_REF_C_NOCAVPOLE_REAL
+EP6ImagChannel: CAL-CS_TDEP_PCALY_LINE2_REF_C_NOCAVPOLE_IMAG
+EP7RealChannel: CAL-CS_TDEP_PCALY_LINE2_REF_D_REAL
+EP7ImagChannel: CAL-CS_TDEP_PCALY_LINE2_REF_D_IMAG
+EP8RealChannel: CAL-CS_TDEP_PCALY_LINE2_REF_A_TST_REAL
+EP8Imagchannel: CAL-CS_TDEP_PCALY_LINE2_REF_A_TST_IMAG
+EP9RealChannel: CAL-CS_TDEP_PCALY_LINE2_REF_A_USUM_REAL
+EP9ImagChannel: CAL-CS_TDEP_PCALY_LINE2_REF_A_USUM_IMAG
+EP10RealChannel: CAL-CS_TDEP_ESD_LINE1_REF_A_TST_NOLOCK_REAL
+EP10ImagChannel: CAL-CS_TDEP_ESD_LINE1_REF_A_TST_NOLOCK_IMAG
+EP11RealChannel: CAL-CS_TDEP_PCALY_LINE4_REF_C_NOCAVPOLE_REAL
+EP11ImagChannel: CAL-CS_TDEP_PCALY_LINE4_REF_C_NOCAVPOLE_IMAG
+EP12RealChannel: CAL-CS_TDEP_PCALY_LINE4_REF_D_REAL
+EP12ImagChannel: CAL-CS_TDEP_PCALY_LINE4_REF_D_IMAG
+EP13RealChannel: CAL-CS_TDEP_PCALY_LINE4_REF_A_TST_REAL
+EP13ImagChannel: CAL-CS_TDEP_PCALY_LINE4_REF_A_TST_IMAG
+EP14RealChannel: CAL-CS_TDEP_PCALY_LINE4_REF_A_USUM_REAL
+EP14ImagChannel: CAL-CS_TDEP_PCALY_LINE4_REF_A_USUM_IMAG
+
+[SampleRates]
+# Sample rate at which to compute h(t)
+HoftSR: 16384
+# Sample rate at which to compute CALIB_STATE_VECTOR
+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 TST excitation channel
+TSTExcSR: 512
+# Sample rate of coherence channels
+CohSR: 16
+# Sample rate for the EPICS reference channels
+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
+# Sample rates at which to compute and record TDCFs
+ComputeFactorsSR: 16
+RecordFactorsSR: 16
+
+[Bitmasks]
+ObsReadyBitmask: 4
+ObsIntentBitmask: 2
+CBCHWInjBitmask: 16777216
+BurstHWInjBitmask: 33554432
+DetCharHWInjBitmask: 67108864
+StochHWInjBitmask: 8388608
+
+[PipelineConfigurations]
+BufferLength: 1.0
+FrequencyDomainFiltering: No
+Dewhitening: No
+# Latency of all filtering/averaging/median processes (other than calibration model filters) as a fraction of filter length. Value should be set between 0.0 and 1.0.
+FilterLatency: 0.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
+RemoveCalLines: Yes
+# Subtract the power lines from the h(t) spectrum
+RemovePowerLines: Yes
+####################################################
+# Options for running power mains line subtraction #
+####################################################
+# 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
+PowerLinesTFAveragingTime: 128
+#######################################
+# Options for broadband noise removal #
+#######################################
+# The length in seconds of the fast Fourier transforms used to compute transfer functions between witness channels and h(t). The fft's are windowed with Hann windows and overlapped.
+WitnessChannelFFTTime: 4.0
+# The number of ffts to take before averaging the witness -> h(t) transfer functions calculation. The average is taken after the ratio h(f) / witness(f).
+NumWitnessFFTs: 1800
+# Sets the minimum number of FFTs necessary to produce the first transfer functions and clean data after data flow starts.
+MinWitnessFFTs: 400
+# The length in seconds of the filters applied to the witness channels before subtracting from h(t)
+WitnessFIRLength: 0.5
+# 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
+# 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: 495.0,515.0,985.0,1015.0;495.0,515.0,985.0,1015.0
+# 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: 3600
+# The amount of time to use to taper in newly computed FIR filters for witness channels being used for noise subtraction.
+WitnessFilterTaperTime: 10
+# Should the transfer function calculation use a median? If not, an average (mean) is used.
+WitnessTFUseMedian: Yes
+###############################
+# Options for HOFT_CLEAN bits #
+###############################
+# The amount of data from h(t) and cleaned h(t) that is used to compute and compare the rms. This comparison between cleaned and uncleaned h(t) determines whether the HOFT_CLEAN bits of the calibration state vector are on or off.
+CleaningCheckRMSTime: 1.0
+# Minimum of a range of frequencies in which we expect line/noise subtraction to be impactful. The HOFT_CLEAN_LOWFREQ_OK bit of the calibration state vector is determined based on whether rms of the cleaned data is less than that of uncleaned h(t) in this range.
+CleaningCheckRangeLowMin: 15
+# Maximum of a range of frequencies in which we expect line/noise subtraction to be impactful. The HOFT_CLEAN_LOWFREQ_OK bit of the calibration state vector is determined based on whether rms of the cleaned data is less than that of uncleaned h(t) in this range.
+CleaningCheckRangeLowMax: 40
+# Minimum of a range of frequencies in which we expect line/noise subtraction to be impactful. The HOFT_CLEAN_MIDFREQ_OK bit of the calibration state vector is determined based on whether rms of the cleaned data is less than that of uncleaned h(t) in this range.
+CleaningCheckRangeMidMin: 100
+# Maximum of a range of frequencies in which we expect line/noise subtraction to be impactful. The HOFT_CLEAN_MIDFREQ_OK bit of the calibration state vector is determined based on whether rms of the cleaned data is less than that of uncleaned h(t) in this range.
+CleaningCheckRangeMidMax: 500
diff --git a/gstlal-calibration/config_files/H1GDS_LowLatency_AllCorrections_Cleaning.ini b/gstlal-calibration/config_files/H1GDS_LowLatency_AllCorrections_Cleaning.ini
index f85b94800a8796eb8c008a4f05f25949f12fd944..7f8ac8c57d3df5de0b6d66ab40cc8cd0e1fe98aa 100644
--- a/gstlal-calibration/config_files/H1GDS_LowLatency_AllCorrections_Cleaning.ini
+++ b/gstlal-calibration/config_files/H1GDS_LowLatency_AllCorrections_Cleaning.ini
@@ -1,6 +1,6 @@
[InputConfigurations]
# Filters file containing calibration FIR filters
-FiltersFileName: ../../O2filters/GDSFilters/H1GDS_minlatency_1175954418.npz
+FiltersFileName: Filters/O2filters/GDSFilters/H1GDS_minlatency_1175954418.npz
# Data source should be set to frames or lvshm
DataSource: frames
FileChecksum: No
@@ -9,7 +9,7 @@ SkipBadFiles: No
############################################
# If reading from frames use these options #
############################################
-FrameCache: H1_raw_frames.cache
+FrameCache: easy_raw_frames.cache
###################################################
# If reading from shared memory use these options #
###################################################
@@ -288,6 +288,8 @@ PowerLinesTFAveragingTime: 128
WitnessChannelFFTTime: 4.0
# The number of ffts to take before averaging the witness -> h(t) transfer functions calculation. The average is taken after the ratio h(f) / witness(f).
NumWitnessFFTs: 1800
+# Sets the minimum number of FFTs necessary to produce the first transfer functions and clean data after data flow starts.
+MinWitnessFFTs: 400
# The length in seconds of the filters applied to the witness channels before subtracting from h(t)
WitnessFIRLength: 0.5
# 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.
diff --git a/gstlal-calibration/gst/lal/gstlal_transferfunction.c b/gstlal-calibration/gst/lal/gstlal_transferfunction.c
index 1a7367eb572e45dc477c9cb2837dd5ff8e3011ca..2b2f0ebde7b546c4488945553594bc1a66f3e3b5 100644
--- a/gstlal-calibration/gst/lal/gstlal_transferfunction.c
+++ b/gstlal-calibration/gst/lal/gstlal_transferfunction.c
@@ -111,6 +111,7 @@ enum property {
ARG_FFT_LENGTH = 1,
ARG_FFT_OVERLAP,
ARG_NUM_FFTS,
+ ARG_MIN_FFTS,
ARG_USE_MEDIAN,
ARG_UPDATE_SAMPLES,
ARG_UPDATE_AFTER_GAP,
@@ -305,10 +306,10 @@ static void update_median_ ## DTYPE(DTYPE *median_array, DTYPE new_element, gint
i--; \
} \
median_array[i + 1] = new_element; \
- /* We want to shift the location of the median up one slot if the new number of samples is even */ \
- if(num_in_median % 2) \
- (*index_median)++; \
} \
+ /* We want to shift the location of the median up one slot if the new number of samples is even */ \
+ if(num_in_median % 2) \
+ (*index_median)++; \
} else { \
/* The array is full and the new value is smaller than the central value, so we start at the beginning */ \
if(new_element > *median_array) { \
@@ -318,10 +319,10 @@ static void update_median_ ## DTYPE(DTYPE *median_array, DTYPE new_element, gint
i++; \
} \
median_array[i - 1] = new_element; \
- /* We want to shift the location of the median down one slot if the new number of samples is odd */ \
- if(!(num_in_median % 2)) \
- (*index_median)--; \
} \
+ /* We want to shift the location of the median down one slot if the new number of samples is odd */ \
+ if(!(num_in_median % 2)) \
+ (*index_median)--; \
} \
\
return; \
@@ -504,7 +505,7 @@ DEFINE_UPDATE_FIR_FILTERS(double, );
#define DEFINE_FIND_TRANSFER_FUNCTION(DTYPE, S_OR_D, F_OR_BLANK) \
-static gboolean find_transfer_functions_ ## DTYPE(GSTLALTransferFunction *element, DTYPE *src, guint64 src_size) { \
+static gboolean find_transfer_functions_ ## DTYPE(GSTLALTransferFunction *element, DTYPE *src, guint64 src_size, guint64 pts) { \
\
gboolean success = TRUE; \
\
@@ -790,7 +791,7 @@ static gboolean find_transfer_functions_ ## DTYPE(GSTLALTransferFunction *elemen
element->workspace.w ## S_OR_D ## pf.num_leftover = maximum64(0, element->sample_count + 1 - sample_count_next_fft); \
\
/* Finally, update transfer functions if ready */ \
- if(element->workspace.w ## S_OR_D ## pf.num_ffts_in_avg == element->num_ffts) { \
+ if(element->workspace.w ## S_OR_D ## pf.num_ffts_in_avg == element->num_ffts || (element->workspace.w ## S_OR_D ## pf.num_ffts_in_avg - element->workspace.w ## S_OR_D ## pf.num_ffts_dropped >= element->min_ffts && *element->transfer_functions == 0.0)) { \
if(element->use_median) { \
/* Then we still need to fill the autocorrelation matrix with median values */ \
gint64 median_length = element->num_ffts / 2 + 1; \
@@ -801,12 +802,17 @@ static gboolean find_transfer_functions_ ## DTYPE(GSTLALTransferFunction *elemen
/* Length of median array is odd, so the median is just the middle value */ \
for(i = 0; i < fd_fft_length; i++) { \
for(j = 0; j < num_tfs; j++) { \
- for(k = 0; k < num_tfs; k++) { \
- /* Off-diagonal elements come from the median array */ \
- median_array_num = i * medians_per_freq + j * num_tfs + k; \
- element->workspace.w ## S_OR_D ## pf.autocorrelation_matrix[i * elements_per_freq * j * (num_tfs + 1) + k] = element->workspace.w ## S_OR_D ## pf.autocorrelation_median_real[median_array_num * median_length + element->workspace.w ## S_OR_D ## pf.index_median_real[median_array_num]] + I * element->workspace.w ## S_OR_D ## pf.autocorrelation_median_imag[median_array_num * median_length + element->workspace.w ## S_OR_D ## pf.index_median_imag[median_array_num]]; \
- /* Elements along the diagonal are one */ \
- element->workspace.w ## S_OR_D ## pf.autocorrelation_matrix[i * elements_per_freq * j * (num_tfs + 1) + num_tfs] = 1.0; \
+ /* First, the ratios fft(first channel) / fft(other channels) */ \
+ median_array_num = i * medians_per_freq + j; \
+ element->workspace.w ## S_OR_D ## pf.autocorrelation_matrix[i * elements_per_freq + j] = element->workspace.w ## S_OR_D ## pf.autocorrelation_median_real[median_array_num * median_length + element->workspace.w ## S_OR_D ## pf.index_median_real[median_array_num]] + I * element->workspace.w ## S_OR_D ## pf.autocorrelation_median_imag[median_array_num * median_length + element->workspace.w ## S_OR_D ## pf.index_median_imag[median_array_num]]; \
+ /* Elements along the diagonal are one */ \
+ element->workspace.w ## S_OR_D ## pf.autocorrelation_matrix[i * elements_per_freq + num_tfs + j * (num_tfs + 1)] = 1.0; \
+ if(j < num_tfs - 1) { \
+ for(k = 0; k < num_tfs; k++) { \
+ /* Off-diagonal elements come from the median array */ \
+ median_array_num = i * medians_per_freq + (j + 1) * num_tfs + k; \
+ element->workspace.w ## S_OR_D ## pf.autocorrelation_matrix[i * elements_per_freq + num_tfs + 1 + j * (num_tfs + 1) + k] = element->workspace.w ## S_OR_D ## pf.autocorrelation_median_real[median_array_num * median_length + element->workspace.w ## S_OR_D ## pf.index_median_real[median_array_num]] + I * element->workspace.w ## S_OR_D ## pf.autocorrelation_median_imag[median_array_num * median_length + element->workspace.w ## S_OR_D ## pf.index_median_imag[median_array_num]]; \
+ } \
} \
} \
} \
@@ -814,12 +820,17 @@ static gboolean find_transfer_functions_ ## DTYPE(GSTLALTransferFunction *elemen
/* Length of median array is even, so the median is the average of the two middle values */ \
for(i = 0; i < fd_fft_length; i++) { \
for(j = 0; j < num_tfs; j++) { \
- for(k = 0; k < num_tfs; k++) { \
- /* Off-diagonal elements come from the median array */ \
- median_array_num = i * medians_per_freq + j * num_tfs + k; \
- element->workspace.w ## S_OR_D ## pf.autocorrelation_matrix[i * elements_per_freq * j * (num_tfs + 1) + k] = (element->workspace.w ## S_OR_D ## pf.autocorrelation_median_real[median_array_num * median_length + element->workspace.w ## S_OR_D ## pf.index_median_real[median_array_num] - 1] + element->workspace.w ## S_OR_D ## pf.autocorrelation_median_real[median_array_num * median_length + element->workspace.w ## S_OR_D ## pf.index_median_real[median_array_num]] + I * (element->workspace.w ## S_OR_D ## pf.autocorrelation_median_imag[median_array_num * median_length + element->workspace.w ## S_OR_D ## pf.index_median_imag[median_array_num] - 1] + element->workspace.w ## S_OR_D ## pf.autocorrelation_median_imag[median_array_num * median_length + element->workspace.w ## S_OR_D ## pf.index_median_imag[median_array_num]])) / 2.0; \
- /* Elements along the diagonal are one */ \
- element->workspace.w ## S_OR_D ## pf.autocorrelation_matrix[i * elements_per_freq * j * (num_tfs + 1) + num_tfs] = 1.0; \
+ /* First, the ratios fft(first channel) / fft(other channels) */ \
+ median_array_num = i * medians_per_freq + j; \
+ element->workspace.w ## S_OR_D ## pf.autocorrelation_matrix[i * elements_per_freq + j] = (element->workspace.w ## S_OR_D ## pf.autocorrelation_median_real[median_array_num * median_length + element->workspace.w ## S_OR_D ## pf.index_median_real[median_array_num] - 1] + element->workspace.w ## S_OR_D ## pf.autocorrelation_median_real[median_array_num * median_length + element->workspace.w ## S_OR_D ## pf.index_median_real[median_array_num]] + I * (element->workspace.w ## S_OR_D ## pf.autocorrelation_median_imag[median_array_num * median_length + element->workspace.w ## S_OR_D ## pf.index_median_imag[median_array_num] - 1] + element->workspace.w ## S_OR_D ## pf.autocorrelation_median_imag[median_array_num * median_length + element->workspace.w ## S_OR_D ## pf.index_median_imag[median_array_num]])) / 2.0; \
+ /* Elements along the diagonal are one */ \
+ element->workspace.w ## S_OR_D ## pf.autocorrelation_matrix[i * elements_per_freq + num_tfs + j * (num_tfs + 1)] = 1.0; \
+ if(j < num_tfs - 1) { \
+ for(k = 0; k < num_tfs; k++) { \
+ /* Off-diagonal elements come from the median array */ \
+ median_array_num = i * medians_per_freq + (j + 1) * num_tfs + k; \
+ element->workspace.w ## S_OR_D ## pf.autocorrelation_matrix[i * elements_per_freq + num_tfs + 1 + j * (num_tfs + 1) + k] = (element->workspace.w ## S_OR_D ## pf.autocorrelation_median_real[median_array_num * median_length + element->workspace.w ## S_OR_D ## pf.index_median_real[median_array_num] - 1] + element->workspace.w ## S_OR_D ## pf.autocorrelation_median_real[median_array_num * median_length + element->workspace.w ## S_OR_D ## pf.index_median_real[median_array_num]] + I * (element->workspace.w ## S_OR_D ## pf.autocorrelation_median_imag[median_array_num * median_length + element->workspace.w ## S_OR_D ## pf.index_median_imag[median_array_num] - 1] + element->workspace.w ## S_OR_D ## pf.autocorrelation_median_imag[median_array_num * median_length + element->workspace.w ## S_OR_D ## pf.index_median_imag[median_array_num]])) / 2.0; \
+ } \
} \
} \
} \
@@ -841,10 +852,14 @@ static gboolean find_transfer_functions_ ## DTYPE(GSTLALTransferFunction *elemen
element->num_tfs_since_gap++; \
} else \
GST_WARNING_OBJECT(element, "Transfer function(s) computation failed. Trying again."); \
- element->sample_count -= element->update_samples + element->num_ffts * stride + element->fft_overlap; \
- element->workspace.w ## S_OR_D ## pf.num_ffts_in_avg = 0; \
- element->workspace.w ## S_OR_D ## pf.num_ffts_dropped = 0; \
\
+ if(element->workspace.w ## S_OR_D ## pf.num_ffts_in_avg == element->num_ffts) { \
+ element->sample_count = (gint64) (gst_util_uint64_scale_int_round(pts, element->rate, GST_SECOND) + src_size) % (element->update_samples + element->num_ffts * stride + element->fft_overlap); \
+ if(element->sample_count > element->update_samples) \
+ element->sample_count -= element->update_samples + element->num_ffts * stride + element->fft_overlap; \
+ element->workspace.w ## S_OR_D ## pf.num_ffts_in_avg = 0; \
+ element->workspace.w ## S_OR_D ## pf.num_ffts_dropped = 0; \
+ } \
/* Update FIR filters if we want */ \
if(success && element->make_fir_filters) { \
success &= update_fir_filters_ ## DTYPE(element->transfer_functions, num_tfs, element->fir_length, element->rate, element->workspace.w ## S_OR_D ## pf.fir_filter, element->workspace.w ## S_OR_D ## pf.fir_plan, element->workspace.w ## S_OR_D ## pf.fir_window, element->workspace.w ## S_OR_D ## pf.tukey, element->fir_filters); \
@@ -928,157 +943,6 @@ static gboolean start(GstBaseSink *sink) {
}
-/*
- * event()
- */
-
-
-static gboolean event(GstBaseSink *sink, GstEvent *event) {
- GSTLALTransferFunction *element = GSTLAL_TRANSFERFUNCTION(sink);
- gboolean success = TRUE;
- GST_DEBUG_OBJECT(element, "Got %s event on sink pad", GST_EVENT_TYPE_NAME(event));
-
- if(GST_EVENT_TYPE(event) == GST_EVENT_EOS) {
- /* If End Of Stream is here and we have not yet computed any transfer functions, we should use whatever data we have to compute them now. */
- if(element->data_type == GSTLAL_TRANSFERFUNCTION_F32) {
- if(element->workspace.wspf.num_ffts_in_avg > element->workspace.wspf.num_ffts_dropped && *element->transfer_functions == 0.0) {
- /* Then we should use the transfer functions we have now, since they won't get any better */
- gint64 fd_fft_length = element->fft_length / 2 + 1;
- gint64 fd_tf_length = element->fir_length / 2 + 1;
- gint64 num_tfs = element->channels - 1;
- if(element->use_median) {
- /* Then we still need to fill the autocorrelation matrix with median values */
- gint64 median_length = element->num_ffts / 2 + 1;
- int elements_per_freq = element->channels * num_tfs;
- int medians_per_freq = num_tfs * num_tfs;
- gint64 median_array_num;
- gint64 i, j, k;
- if((element->workspace.wspf.num_ffts_in_avg - element->workspace.wspf.num_ffts_dropped) % 2) {
- /* Length of median array is odd, so the median is just the middle value */
- for(i = 0; i < fd_fft_length; i++) {
- for(j = 0; j < num_tfs; j++) {
- for(k = 0; k < num_tfs; k++) {
- /* Off-diagonal elements come from the median array */
- median_array_num = i * medians_per_freq + j * num_tfs + k;
- element->workspace.wspf.autocorrelation_matrix[i * elements_per_freq * j * (num_tfs + 1) + k] = element->workspace.wspf.autocorrelation_median_real[median_array_num * median_length + element->workspace.wspf.index_median_real[median_array_num]] + I * element->workspace.wspf.autocorrelation_median_imag[median_array_num * median_length + element->workspace.wspf.index_median_imag[median_array_num]];
- /* Elements along the diagonal are one */
- element->workspace.wspf.autocorrelation_matrix[i * elements_per_freq * j * (num_tfs + 1) + num_tfs] = 1.0;
- }
- }
- }
- } else {
- /* Length of median array is even, so the median is the average of the two middle values */
- for(i = 0; i < fd_fft_length; i++) {
- for(j = 0; j < num_tfs; j++) {
- for(k = 0; k < num_tfs; k++) {
- /* Off-diagonal elements come from the median array */
- median_array_num = i * medians_per_freq + j * num_tfs + k;
- element->workspace.wspf.autocorrelation_matrix[i * elements_per_freq * j * (num_tfs + 1) + k] = (element->workspace.wspf.autocorrelation_median_real[median_array_num * median_length + element->workspace.wspf.index_median_real[median_array_num] - 1] + element->workspace.wspf.autocorrelation_median_real[median_array_num * median_length + element->workspace.wspf.index_median_real[median_array_num]] + I * (element->workspace.wspf.autocorrelation_median_imag[median_array_num * median_length + element->workspace.wspf.index_median_imag[median_array_num] - 1] + element->workspace.wspf.autocorrelation_median_imag[median_array_num * median_length + element->workspace.wspf.index_median_imag[median_array_num]])) / 2.0;
- /* Elements along the diagonal are one */
- element->workspace.wspf.autocorrelation_matrix[i * elements_per_freq * j * (num_tfs + 1) + num_tfs] = 1.0;
- }
- }
- }
- }
- }
- success &= update_transfer_functions_float(element->workspace.wspf.autocorrelation_matrix, num_tfs, fd_fft_length, fd_tf_length, element->workspace.wspf.sinc_table, element->workspace.wspf.sinc_length, element->workspace.wspf.sinc_taps_per_df, element->use_median ? 1 : element->num_ffts - element->workspace.wspf.num_ffts_dropped, element->workspace.wspf.transfer_functions_at_f, element->workspace.wspf.transfer_functions_solved_at_f, element->workspace.wspf.autocorrelation_matrix_at_f, element->workspace.wspf.permutation, element->transfer_functions);
- if(success) {
- GST_INFO_OBJECT(element, "Just computed new transfer functions");
- /* Let other elements know about the update */
- g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_TRANSFER_FUNCTIONS]);
- /* Write transfer functions to the screen or a file if we want */
- if(element->write_to_screen || element->filename)
- write_transfer_functions(element->transfer_functions, gst_element_get_name(element), element->rate / 2.0 / (fd_tf_length - 1.0), fd_tf_length, num_tfs, element->write_to_screen, element->filename, TRUE);
- } else
- GST_WARNING_OBJECT(element, "Transfer function(s) computation failed. No transfer functions will be produced.");
- /* Update FIR filters if we want */
- if(success && element->make_fir_filters) {
- success &= update_fir_filters_float(element->transfer_functions, num_tfs, element->fir_length, element->rate, element->workspace.wspf.fir_filter, element->workspace.wspf.fir_plan, element->workspace.wspf.fir_window, element->workspace.wspf.tukey, element->fir_filters);
- if(success) {
- GST_INFO_OBJECT(element, "Just computed new FIR filters");
- /* Let other elements know about the update */
- g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_FIR_FILTERS]);
- /* Write FIR filters to the screen or a file if we want */
- if(element->write_to_screen || element->filename)
- write_fir_filters(element->fir_filters, gst_element_get_name(element), element->fir_length, num_tfs, element->write_to_screen, element->filename, TRUE);
- } else
- GST_WARNING_OBJECT(element, "FIR filter(s) computation failed. No FIR filters will be produced.");
- }
- }
- } else {
- if(element->workspace.wdpf.num_ffts_in_avg > element->workspace.wdpf.num_ffts_dropped && *element->transfer_functions == 0.0) {
- /* Then we should use the transfer functions we have now, since they won't get any better */
- gint64 fd_fft_length = element->fft_length / 2 + 1;
- gint64 fd_tf_length = element->fir_length / 2 + 1;
- gint64 num_tfs = element->channels - 1;
- if(element->use_median) {
- /* Then we still need to fill the autocorrelation matrix with median values */
- gint64 median_length = element->num_ffts / 2 + 1;
- int elements_per_freq = element->channels * num_tfs;
- int medians_per_freq = num_tfs * num_tfs;
- gint64 median_array_num;
- gint64 i, j, k;
- if((element->workspace.wdpf.num_ffts_in_avg - element->workspace.wdpf.num_ffts_dropped) % 2) {
- /* Length of median array is odd, so the median is just the middle value */
- for(i = 0; i < fd_fft_length; i++) {
- for(j = 0; j < num_tfs; j++) {
- for(k = 0; k < num_tfs; k++) {
- /* Off-diagonal elements come from the median array */
- median_array_num = i * medians_per_freq + j * num_tfs + k;
- element->workspace.wdpf.autocorrelation_matrix[i * elements_per_freq * j * (num_tfs + 1) + k] = element->workspace.wdpf.autocorrelation_median_real[median_array_num * median_length + element->workspace.wdpf.index_median_real[median_array_num]] + I * element->workspace.wdpf.autocorrelation_median_imag[median_array_num * median_length + element->workspace.wdpf.index_median_imag[median_array_num]];
- /* Elements along the diagonal are one */
- element->workspace.wdpf.autocorrelation_matrix[i * elements_per_freq * j * (num_tfs + 1) + num_tfs] = 1.0;
- }
- }
- }
- } else {
- /* Length of median array is even, so the median is the average of the two middle values */
- for(i = 0; i < fd_fft_length; i++) {
- for(j = 0; j < num_tfs; j++) {
- for(k = 0; k < num_tfs; k++) {
- /* Off-diagonal elements come from the median array */
- median_array_num = i * medians_per_freq + j * num_tfs + k;
- element->workspace.wdpf.autocorrelation_matrix[i * elements_per_freq * j * (num_tfs + 1) + k] = (element->workspace.wdpf.autocorrelation_median_real[median_array_num * median_length + element->workspace.wdpf.index_median_real[median_array_num] - 1] + element->workspace.wdpf.autocorrelation_median_real[median_array_num * median_length + element->workspace.wdpf.index_median_real[median_array_num]] + I * (element->workspace.wdpf.autocorrelation_median_imag[median_array_num * median_length + element->workspace.wdpf.index_median_imag[median_array_num] - 1] + element->workspace.wdpf.autocorrelation_median_imag[median_array_num * median_length + element->workspace.wdpf.index_median_imag[median_array_num]])) / 2.0;
- /* Elements along the diagonal are one */
- element->workspace.wdpf.autocorrelation_matrix[i * elements_per_freq * j * (num_tfs + 1) + num_tfs] = 1.0;
- }
- }
- }
- }
- }
- success &= update_transfer_functions_double(element->workspace.wdpf.autocorrelation_matrix, num_tfs, fd_fft_length, fd_tf_length, element->workspace.wdpf.sinc_table, element->workspace.wdpf.sinc_length, element->workspace.wdpf.sinc_taps_per_df, element->use_median ? 1 : element->num_ffts - element->workspace.wdpf.num_ffts_dropped, element->workspace.wdpf.transfer_functions_at_f, element->workspace.wdpf.transfer_functions_solved_at_f, element->workspace.wdpf.autocorrelation_matrix_at_f, element->workspace.wdpf.permutation, element->transfer_functions);
- if(success) {
- GST_INFO_OBJECT(element, "Just computed new transfer functions");
- /* Let other elements know about the update */
- g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_TRANSFER_FUNCTIONS]);
- /* Write transfer functions to the screen or a file if we want */
- if(element->write_to_screen || element->filename)
- write_transfer_functions(element->transfer_functions, gst_element_get_name(element), element->rate / 2.0 / (fd_tf_length - 1.0), fd_tf_length, num_tfs, element->write_to_screen, element->filename, TRUE);
- } else
- GST_WARNING_OBJECT(element, "Transfer function(s) computation failed. No transfer functions will be produced.");
- /* Update FIR filters if we want */
- if(success && element->make_fir_filters) {
- success &= update_fir_filters_double(element->transfer_functions, num_tfs, element->fir_length, element->rate, element->workspace.wdpf.fir_filter, element->workspace.wdpf.fir_plan, element->workspace.wdpf.fir_window, element->workspace.wdpf.tukey, element->fir_filters);
- if(success) {
- GST_INFO_OBJECT(element, "Just computed new FIR filters");
- /* Let other elements know about the update */
- g_object_notify_by_pspec(G_OBJECT(element), properties[ARG_FIR_FILTERS]);
- /* Write FIR filters to the screen or a file if we want */
- if(element->write_to_screen || element->filename)
- write_fir_filters(element->fir_filters, gst_element_get_name(element), element->fir_length, num_tfs, element->write_to_screen, element->filename, TRUE);
- } else
- GST_WARNING_OBJECT(element, "FIR filter(s) computation failed. No FIR filters will be produced.");
- }
- }
- }
- }
-
- success = GST_BASE_SINK_CLASS(gstlal_transferfunction_parent_class)->event(sink, event);
-
- return success;
-}
-
-
/*
* set_caps()
*/
@@ -1685,9 +1549,11 @@ static GstFlowReturn render(GstBaseSink *sink, GstBuffer *buffer) {
element->sample_count = element->update_samples;
if(element->data_type == GSTLAL_TRANSFERFUNCTION_F32) {
element->workspace.wspf.num_ffts_in_avg = 0;
+ element->workspace.wspf.num_ffts_dropped = 0;
element->workspace.wspf.num_leftover = 0;
} else {
element->workspace.wdpf.num_ffts_in_avg = 0;
+ element->workspace.wdpf.num_ffts_dropped = 0;
element->workspace.wdpf.num_leftover = 0;
}
}
@@ -1734,9 +1600,11 @@ static GstFlowReturn render(GstBaseSink *sink, GstBuffer *buffer) {
element->sample_count = minimum64(element->sample_count + mapinfo.size / element->unit_size, element->update_samples);
if(element->data_type == GSTLAL_TRANSFERFUNCTION_F32) {
element->workspace.wspf.num_ffts_in_avg = 0;
+ element->workspace.wspf.num_ffts_dropped = 0;
element->workspace.wspf.num_leftover = 0;
} else {
element->workspace.wdpf.num_ffts_in_avg = 0;
+ element->workspace.wdpf.num_ffts_dropped = 0;
element->workspace.wdpf.num_leftover = 0;
}
} else {
@@ -1755,14 +1623,14 @@ static GstFlowReturn render(GstBaseSink *sink, GstBuffer *buffer) {
memset(element->workspace.wspf.autocorrelation_matrix, 0, element->channels * (element->channels - 1) * (element->fft_length / 2 + 1) * sizeof(*element->workspace.wspf.autocorrelation_matrix));
/* Send the data to a function to compute fft's and transfer functions */
- success = find_transfer_functions_float(element, (float *) mapinfo.data, mapinfo.size);
+ success = find_transfer_functions_float(element, (float *) mapinfo.data, mapinfo.size, GST_BUFFER_PTS(buffer));
} else {
/* If we are just beginning to compute new transfer functions with this data, initialize memory that we will fill to zero */
if(!element->workspace.wdpf.num_ffts_in_avg)
memset(element->workspace.wdpf.autocorrelation_matrix, 0, element->channels * (element->channels - 1) * (element->fft_length / 2 + 1) * sizeof(*element->workspace.wdpf.autocorrelation_matrix));
/* Send the data to a function to compute fft's and transfer functions */
- success = find_transfer_functions_double(element, (double *) mapinfo.data, mapinfo.size);
+ success = find_transfer_functions_double(element, (double *) mapinfo.data, mapinfo.size, GST_BUFFER_PTS(buffer));
}
if(!success) {
@@ -1770,9 +1638,11 @@ static GstFlowReturn render(GstBaseSink *sink, GstBuffer *buffer) {
element->sample_count = element->update_samples;
if(element->data_type == GSTLAL_TRANSFERFUNCTION_F32) {
element->workspace.wspf.num_ffts_in_avg = 0;
+ element->workspace.wspf.num_ffts_dropped = 0;
element->workspace.wspf.num_leftover = 0;
} else {
element->workspace.wdpf.num_ffts_in_avg = 0;
+ element->workspace.wdpf.num_ffts_dropped = 0;
element->workspace.wdpf.num_leftover = 0;
}
}
@@ -1918,6 +1788,10 @@ static void set_property(GObject *object, enum property id, const GValue *value,
element->num_ffts = g_value_get_int64(value);
break;
+ case ARG_MIN_FFTS:
+ element->min_ffts = g_value_get_int64(value);
+ break;
+
case ARG_USE_MEDIAN:
element->use_median = g_value_get_boolean(value);
break;
@@ -2019,6 +1893,10 @@ static void get_property(GObject *object, enum property id, GValue *value, GPara
g_value_set_int64(value, element->num_ffts);
break;
+ case ARG_MIN_FFTS:
+ g_value_set_int64(value, element->min_ffts);
+ break;
+
case ARG_USE_MEDIAN:
g_value_set_boolean(value, element->use_median);
break;
@@ -2160,7 +2038,6 @@ static void gstlal_transferfunction_class_init(GSTLALTransferFunctionClass *klas
GstBaseSinkClass *gstbasesink_class = GST_BASE_SINK_CLASS(klass);
gstbasesink_class->start = GST_DEBUG_FUNCPTR(start);
- gstbasesink_class->event = GST_DEBUG_FUNCPTR(event);
gstbasesink_class->set_caps = GST_DEBUG_FUNCPTR(set_caps);
gstbasesink_class->render = GST_DEBUG_FUNCPTR(render);
gstbasesink_class->stop = GST_DEBUG_FUNCPTR(stop);
@@ -2213,6 +2090,14 @@ static void gstlal_transferfunction_class_init(GSTLALTransferFunctionClass *klas
1, G_MAXINT64, 16,
G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
);
+ properties[ARG_MIN_FFTS] = g_param_spec_int64(
+ "min-ffts",
+ "Minimum number of FFTs",
+ "If EOS is reached before transfer functions are computed, this sets the\n\t\t\t"
+ "minimum number of FFTs necessary to produce transfer function at EOS.",
+ 1, G_MAXINT64, 16,
+ G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT
+ );
properties[ARG_USE_MEDIAN] = g_param_spec_boolean(
"use-median",
"Use Median",
@@ -2381,6 +2266,11 @@ static void gstlal_transferfunction_class_init(GSTLALTransferFunctionClass *klas
ARG_NUM_FFTS,
properties[ARG_NUM_FFTS]
);
+ g_object_class_install_property(
+ gobject_class,
+ ARG_MIN_FFTS,
+ properties[ARG_MIN_FFTS]
+ );
g_object_class_install_property(
gobject_class,
ARG_USE_MEDIAN,
diff --git a/gstlal-calibration/gst/lal/gstlal_transferfunction.h b/gstlal-calibration/gst/lal/gstlal_transferfunction.h
index 49bfd37e387d09e1765be82c077caa06bdfc17e6..4ce6c403f73d51a0fc1d3bba9689363093de6a73 100644
--- a/gstlal-calibration/gst/lal/gstlal_transferfunction.h
+++ b/gstlal-calibration/gst/lal/gstlal_transferfunction.h
@@ -144,6 +144,7 @@ struct _GSTLALTransferFunction {
gint64 fft_length;
gint64 fft_overlap;
gint64 num_ffts;
+ gint64 min_ffts;
gboolean use_median;
gint64 update_samples;
gboolean update_after_gap;
diff --git a/gstlal-calibration/python/calibration_parts.py b/gstlal-calibration/python/calibration_parts.py
index ff2f115e03878b773b536848dd4c2927f5b080a3..ef2fd923348257eed17196d8200f87985ba32825 100644
--- a/gstlal-calibration/python/calibration_parts.py
+++ b/gstlal-calibration/python/calibration_parts.py
@@ -758,7 +758,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, update_samples, fir_length, frequency_resolution, filter_taper_length, notch_frequencies = [], obsready = None, chop_time = 0.0, wait_time = 0.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, notch_frequencies = [], obsready = None, chop_time = 0.0, wait_time = 0.0, filename = None):
#
# Use witness channels that monitor the environment to remove environmental noise
@@ -781,8 +781,8 @@ def clean_data(pipeline, signal, signal_rate, witnesses, witness_rate, fft_lengt
if obsready is not None:
transfer_functions = mkgate(pipeline, transfer_functions, obsready, 1, attack_length = -witness_rate * chop_time)
if(chop_time):
- transfer_functions = pipeparts.mkgeneric(pipeline, transfer_functions, "lal_insertgap", chop_length = int(1000000000 * chop_time))
- transfer_functions = pipeparts.mkgeneric(pipeline, transfer_functions, "lal_transferfunction", fft_length = fft_length, fft_overlap = fft_overlap, num_ffts = num_ffts, update_samples = update_samples, make_fir_filters = -1, fir_length = fir_length, frequency_resolution = frequency_resolution, high_pass = 15, update_after_gap = True, notch_frequencies = notch_frequencies, filename = filename)
+ transfer_functions = pipeparts.mkgeneric(pipeline, transfer_functions, "lal_insertgap", chop_length = int(1000000000 * 10))
+ 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, notch_frequencies = notch_frequencies, filename = filename)
signal_minus_noise = [signal_tee]
for i in range(0, len(witnesses)):
minus_noise = pipeparts.mkgeneric(pipeline, mkqueue(pipeline, witness_tees[i], min_length = wait_time), "lal_tdwhiten", kernel = default_fir_filter, latency = fir_length / 2, taper_length = filter_taper_length)
diff --git a/gstlal-calibration/tests/Makefile.am b/gstlal-calibration/tests/Makefile.am
index 223adf09e530833b7dc5fd498ad4d04542139604..dfc3829d2371cb111a96577061f606956d77082f 100644
--- a/gstlal-calibration/tests/Makefile.am
+++ b/gstlal-calibration/tests/Makefile.am
@@ -1,7 +1,15 @@
+# This is a trick taken from the gst-python automake setup.
+# All of the Python scripts will be installed under the exec dir,
+# which prevents the module from getting spread across lib and lib64
+# on e.g. CentOS.
+pkgpythondir = $(pkgpyexecdir)
+
+pkgpython_PYTHON = \
+ test_common.py
+
EXTRA_DIST = \
L-L1_TEST_FRAME-10-10.gwf \
- lal_logical_undersampler_test_01.py \
- test_common.py
+ lal_logical_undersampler_test_01.py
#TESTS = lal_logical_undersampler_test_01.py
diff --git a/gstlal-calibration/tests/check_calibration/ASD_comparison_plots b/gstlal-calibration/tests/check_calibration/ASD_comparison_plots
new file mode 100755
index 0000000000000000000000000000000000000000..dcb2b4c8f0bb87c741311b3a17816e071a25c0e3
--- /dev/null
+++ b/gstlal-calibration/tests/check_calibration/ASD_comparison_plots
@@ -0,0 +1,66 @@
+#!/usr/bin/env python
+
+import matplotlib
+matplotlib.use('Agg')
+from gwpy.timeseries import TimeSeriesDict
+from gwpy.timeseries import TimeSeries
+import glob
+from math import pi
+from gwpy.plotter import BodePlot
+import numpy
+from optparse import OptionParser, Option
+from glue import datafind
+
+parser = OptionParser()
+
+parser.add_option("--ifo", metavar = "name", help = "Name of the IFO")
+parser.add_option("--gps-start-time", metavar = "seconds", help = "Set the GPS start time.")
+parser.add_option("--gps-end-time", metavar = "seconds", help = "Set the GPS end time.")
+parser.add_option("--raw-frame-cache", metavar = "name", help = "Raw frame cache file")
+parser.add_option("--gds-frame-cache", metavar = "name", help = "GDS frame cache file.")
+parser.add_option("--gds-channel-name", metavar = "name", default = "GDS-CALIB_STRAIN", help = "Channel name for h(t) channel in GDS frames. (Default = GDS-CALIB_STRAIN)")
+parser.add_option("--calcs-channel-name", metavar = "name", default = "CAL-DELTAL_EXTERNAL_DQ", help = "Channel name for h(t) channel in raw frames. (Default = CAL-DELTAL_EXTERNAL_DQ)")
+
+options, filenames = parser.parse_args()
+
+start_time = int(options.gps_start_time)
+end_time = int(options.gps_end_time)
+
+# Grab CALCS data
+calcs_data=TimeSeries.read(options.raw_frame_cache, '%s:%s' % (options.ifo, options.calcs_channel_name), start = start_time, end = end_time)
+
+# grab GDS/DCS data
+gds_data = TimeSeries.read(options.gds_frame_cache, "%s:%s" % (options.ifo, options.gds_channel_name), start = start_time, end = end_time)
+
+# make asds
+calcs_asd = calcs_data.asd(4,2)
+gds_asd = gds_data.asd(4,2)
+
+#dewhiten CALCS
+calcs_asd = calcs_asd.filter([30]*6, [0.3]*6, 1e-12 / 4e3)
+
+#plot spectrum
+plot=calcs_asd.plot(label='CALCS h(t) ASD')
+plot.gca().plot(gds_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)
+ax.set_xlim(10,8192)
+ax.set_ylim(1e-24,1e-16)
+ax.legend()
+plot.save('spectrum_comparison.png')
+
+diff = calcs_asd / gds_asd
+plot = diff.plot(label="ASD ratio CALCS / DCS", logy = False)
+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)
+ax.set_xlim(10,5000)
+ax.set_ylim(0.7, 1.3)
+ax.legend()
+plot.save('CALCS_DCS_residual.png')
+
diff --git a/gstlal-calibration/tests/check_calibration/Makefile b/gstlal-calibration/tests/check_calibration/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..196465400edda982a1bcc03bfa485d28e7410094
--- /dev/null
+++ b/gstlal-calibration/tests/check_calibration/Makefile
@@ -0,0 +1,45 @@
+
+#################################
+### Inputs for user to modify ###
+#################################
+
+START = 1186159035
+END = 1186159557
+IFO = H1
+GDSCONFIGS = ../../config_files/H1GDS_LowLatency_AllCorrections_Cleaning.ini
+DCSCONFIGS = ../../config_files/H1DCS_AllCorrections_Cleaning.ini
+
+all: DCS_pcal2darm_plots
+
+################################################
+### These commands should changes less often ###
+################################################
+
+$(IFO)_raw_frames.cache:
+ gw_data_find -o H -t H1_R -s $(START) -e $(END) -l --url-type file > $@
+
+$(IFO)_C00_frames.cache:
+ gw_data_find -o H -t H1_HOFT_C00 -s $(START) -e $(END) -l --url-type file > $@
+
+$(IFO)_C01_frames.cache:
+ gw_data_find -o H -t H1_HOFT_C01 -s $(START) -e $(END) -l --url-type file > $@
+
+$(IFO)_C02_frames.cache:
+ gw_data_find -o H -t H1_HOFT_C02 -s $(START) -e $(END) -l --url-type file > $@
+
+$(IFO)_hoft_GDS_frames.cache: $(IFO)_raw_frames.cache
+ GST_DEBUG=3 gstlal_compute_strain --gps-start-time $(START) --gps-end-time $(END) --frame-duration=64 --frames-per-file=1 --wings=0 --config-file $(GDSCONFIGS)
+ ls H-H1GDS_TEST*.gwf | lalapps_path2cache > $@
+
+$(IFO)_hoft_DCS_frames.cache: $(IFO)_raw_frames.cache
+ GST_DEBUG=3 gstlal_compute_strain --gps-start-time $(START) --gps-end-time $(END) --frame-duration=64 --frames-per-file=1 --wings=0 --config-file $(DCSCONFIGS)
+ ls H-H1DCS_TEST*.gwf | lalapps_path2cache > $@
+
+GDS_pcal2darm_plots: $(IFO)_raw_frames.cache $(IFO)_hoft_GDS_frames.cache
+ python pcal2darm_timeseries.py --ifo=$(IFO) --raw-frame-cache $(IFO)_raw_frames.cache --calibrated-frame-cache $(IFO)_hoft_GDS_frames.cache --config-file $(GDSCONFIGS) --calibrated-channel-list='GDS-CALIB_STRAIN,GDS-CALIB_STRAIN_CLEAN'
+
+DCS_pcal2darm_plots: $(IFO)_raw_frames.cache $(IFO)_hoft_DCS_frames.cache
+ python pcal2darm_timeseries.py --ifo=$(IFO) --raw-frame-cache $(IFO)_raw_frames.cache --calibrated-frame-cache $(IFO)_hoft_DCS_frames.cache --config-file $(DCSCONFIGS) --calibrated-channel-list='DCS-CALIB_STRAIN,DCS-CALIB_STRAIN_CLEAN'
+
+clean:
+ rm *.gwf *.cache *.png
diff --git a/gstlal-calibration/tests/check_calibration/Makefile.DCS_calibration b/gstlal-calibration/tests/check_calibration/Makefile.DCS_calibration
new file mode 100644
index 0000000000000000000000000000000000000000..c44c0d632944da85bbc9bb7d86300ff1cf94b9eb
--- /dev/null
+++ b/gstlal-calibration/tests/check_calibration/Makefile.DCS_calibration
@@ -0,0 +1,32 @@
+# time of GW170817 1187008882
+#SHELL=/bin/bash
+START = 1187007882 # Thu Aug 17 12:24:24 GMT 2017
+END = 1187008582 # Thu Aug 17 12:36:04 GMT 2017
+PLOT_START = $(shell echo $(START) + 200 | bc)
+PLOT_END = $(shell echo $(END) - 100 | bc)
+#PLOTSTART = 1187008082
+#PLOTEND = 1187008482
+FILTER = 'L1DCS_1175961600.npz'
+IFO = L
+#Note: Make sure there is no space after the L or H
+
+all: Filter_for_Frames DCS_CALCS_plots
+
+Filter_for_Frames:
+ cp Filters/GDSFilters/$(FILTER) Frames/DCS_calibration
+
+raw_frames.cache:
+ cd Frames/DCS_calibration; gw_data_find -o $(IFO) -t $(IFO)1_R -s $(START) -e $(END) -l --url-type file > $@
+
+hoft_frames.cache: raw_frames.cache
+ cd Frames/DCS_calibration; gstlal_compute_strain --data-source frames --frame-cache raw_frames.cache --gps-start-time $(START) --gps-end-time $(END) --frame-duration=4 --frames-per-file=1 --filters-file $(FILTER) --ifo $(IFO)1 --frame-type $(IFO)1_TEST --compression-scheme=6 --compression-level=3 --full-calibration --control-sample-rate 16384 --factors-from-filters-file --expected-fcc 376.0 --no-fs --no-srcQ --coherence-uncertainty-threshold 0.02 --kappas-default-to-median --apply-kappatst --apply-kappapu --apply-kappac --verbose
+ cd Frames/DCS_calibration; ls *.gwf | lalapps_path2cache > $@
+
+DCS_CALCS_plots: raw_frames.cache hoft_frames.cache
+ ./ASD_comparison_plots --ifo $(IFO)1 --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --raw-frame-cache Frames/DCS_calibration/raw_frames.cache --gds-frame-cache Frames/DCS_calibration/hoft_frames.cache
+
+
+clean:
+ cd Frames/DCS_calibration; rm *.cache *.gwf
+ rm *.png
+
diff --git a/gstlal-calibration/tests/check_calibration/Makefile.all_tests b/gstlal-calibration/tests/check_calibration/Makefile.all_tests
new file mode 100644
index 0000000000000000000000000000000000000000..c838873925c0479feaef7383fbe7ce0d96406c6d
--- /dev/null
+++ b/gstlal-calibration/tests/check_calibration/Makefile.all_tests
@@ -0,0 +1,20 @@
+# This makefile can be used to run all of the tests instead of running them one by one.
+
+all: response_function DCS_calibration timeserieskappas
+
+response_function:
+ make -f Makefile.response_function
+
+DCS_calibration:
+ make -f Makefile.DCS_calibration
+
+timeserieskappas:
+ make -f Makefile.timeserieskappas
+
+clean:
+ rm *.pdf
+ cd Frames/response_function/C00_no_kappas; rm *.gwf *.cache
+ cd Frames/response_function/C01_no_kappas; rm *.gwf *.cache
+ rm *.png
+ cd Frames/DCS_calibration; rm *.gwf *.cache
+
diff --git a/gstlal-calibration/tests/check_calibration/Makefile.response_function b/gstlal-calibration/tests/check_calibration/Makefile.response_function
new file mode 100644
index 0000000000000000000000000000000000000000..150b0564bb1821fb87960dab73b4cee466f692b7
--- /dev/null
+++ b/gstlal-calibration/tests/check_calibration/Makefile.response_function
@@ -0,0 +1,48 @@
+#Livingston:
+START = 1183905388
+END = 1183906844
+FILTER_C00 = L1GDS_1175954418.npz
+FILTER_C01 = L1DCS_highpass_1175961600.npz
+IFO = L
+
+#Hanford:
+#START = 1186200018
+#END = 1186201218
+#FILTER_C00 = H1GDS_1175954418.npz
+#FILTER_C01 = H1DCS_highpass_1173225472.npz
+#IFO = H
+
+#Note: Make sure that there is no space after the L or H!
+#Another note: Livingston cluster cannot run this for IFO = H. Make sure that the cluster you're using has the raw frames you want.
+
+PLOT_START = $(shell echo $(START) + 408 | bc)
+PLOT_END = $(shell echo $(END) - 444 | bc)
+
+all: Filters_for_Frames C00_hoft_GDS_frames.cache C01_hoft_GDS_frames.cache response_function_bode_plot
+
+Filters_for_Frames:
+ cp Filters/GDSFilters/$(FILTER_C00) Frames/response_function/C00_no_kappas
+ cp Filters/GDSFilters/$(FILTER_C01) Frames/response_function/C01_no_kappas
+
+C00_raw_frames.cache:
+ cd Frames/response_function/C00_no_kappas; gw_data_find -o $(IFO) -t $(IFO)1_R -s $(START) -e $(END) -l --url-type file > $@
+
+C01_raw_frames.cache:
+ cd Frames/response_function/C01_no_kappas; gw_data_find -o $(IFO) -t $(IFO)1_R -s $(START) -e $(END) -l --url-type file > $@
+
+C00_hoft_GDS_frames.cache: C00_raw_frames.cache
+ cd Frames/response_function/C00_no_kappas; gstlal_compute_strain --data-source frames --frame-cache C00_raw_frames.cache --gps-start-time $(START) --gps-end-time $(END) --frame-duration=4 --frames-per-file=1 --filters-file $(FILTER_C00) --ifo $(IFO)1 --frame-type $(IFO)1_TEST --compression-scheme=6 --compression-level=3 --partial-calibration --control-sample-rate=4096 --no-srcQ --no-fs --wings 400 --no-kappac --no-kappapu --no-kappatst --no-fcc --verbose
+ cd Frames/response_function/C00_no_kappas; ls *.gwf | lalapps_path2cache > $@
+
+C01_hoft_GDS_frames.cache: C01_raw_frames.cache
+ cd Frames/response_function/C01_no_kappas; gstlal_compute_strain --data-source frames --frame-cache C01_raw_frames.cache --gps-start-time $(START) --gps-end-time $(END) --frame-duration=4 --frames-per-file=1 --filters-file $(FILTER_C01) --ifo $(IFO)1 --frame-type $(IFO)1_TEST --compression-scheme=6 --compression-level=3 --full-calibration --no-srcQ --no-fs --wings 400 --no-kappac --no-kappapu --no-kappatst --no-fcc --verbose
+ cd Frames/response_function/C01_no_kappas; ls *.gwf | lalapps_path2cache > $@
+
+response_function_bode_plot:
+ python response_function.py --gps_start_time $(PLOT_START) --gps_end_time $(PLOT_END) --dt 6.103515625e-05 --ifo $(IFO)1 --c00_hoft_frames_cache Frames/response_function/C00_no_kappas/C00_hoft_GDS_frames.cache --c01_hoft_frames_cache Frames/response_function/C01_no_kappas/C01_hoft_GDS_frames.cache --raw_frames_cache Frames/response_function/C01_no_kappas/C01_raw_frames.cache --darm_err_channel_name CAL-DARM_ERR_WHITEN_OUT_DBL_DQ --c00_hoft_channel_name GDS-CALIB_STRAIN --c01_hoft_channel_name GDS-CALIB_STRAIN --calcs_hoft_channel_name CAL-DELTAL_EXTERNAL_DQ --response_file Filters/GDSFilters/$(FILTER_C01)
+
+clean:
+ rm *.pdf
+ cd Frames/response_function/C00_no_kappas; rm *.gwf *.cache
+ cd Frames/response_function/C01_no_kappas; rm *.gwf *.cache
+
diff --git a/gstlal-calibration/tests/check_calibration/Makefile.timeserieskappas b/gstlal-calibration/tests/check_calibration/Makefile.timeserieskappas
new file mode 100644
index 0000000000000000000000000000000000000000..72d74b926695a724721a58f33869eba05363ecb0
--- /dev/null
+++ b/gstlal-calibration/tests/check_calibration/Makefile.timeserieskappas
@@ -0,0 +1,22 @@
+START = 1183161618
+#END = 1183194810
+#END = 1183196880
+END = 1183248018
+CHANNEL1 = L1:GDS-CALIB_KAPPA_C
+CHANNEL2 = L1:GDS-CALIB_KAPPA_PU_REAL
+CHANNEL3 = L1:GDS-CALIB_KAPPA_TST_REAL
+CHANNEL4 = L1:GDS-CALIB_F_CC
+CHANNEL5 = L1:GDS-CALIB_KAPPA_PU_IMAGINARY
+CHANNEL6 = L1:GDS-CALIB_KAPPA_TST_IMAGINARY
+CHANNEL7 = L1:GDS-CALIB_SRC_Q_INVERSE
+CHANNEL8 = L1:GDS-CALIB_FS
+IFO = L
+
+all: Time_Series_of_Kappas
+
+Time_Series_of_Kappas:
+ python timeserieskappas.py --gps_start_time $(START) --gps_end_time $(END) --channel_name1 $(CHANNEL1) --channel_name2 $(CHANNEL2) --channel_name3 $(CHANNEL3) --channel_name4 $(CHANNEL4) --channel_name5 $(CHANNEL5) --channel_name6 $(CHANNEL6)
+
+clean:
+ rm *.png
+
diff --git a/gstlal-calibration/tests/check_calibration/README b/gstlal-calibration/tests/check_calibration/README
new file mode 100644
index 0000000000000000000000000000000000000000..1833f0f996a911b848f58ebf2f940dbedb47dc80
--- /dev/null
+++ b/gstlal-calibration/tests/check_calibration/README
@@ -0,0 +1,12 @@
+As of now, this subdirectory contains three tests of the calibration. Other tests will be added soon.
+
+How to use:
+1. The user will open the Makefile for the test they are interested in running (e.g. vi Makefile.response_function)
+2. The user will enter the desired inputs at the top of the Makefile for that test (e.g. START = 1183905388)
+3. The user will run the makefile (e.g. make -f Makefile.response_function)
+4. The user will look at the plots to check the calibration
+
+NOTE: Filters in the Filters subdirectory might not be up-to-date. Checkout the filters from the SVN if desired.
+Example: svn co --username=albert.einstein https://svn.ligo.caltech.edu/svn/aligocalibration/trunk/Runs/O1/GDSFilters/
+
+ANOTHER NOTE: Make sure that the cluster you are using has the raw frames for whichever interferometer you choose. For example, IFO = H in Makefile.response_function would not work on the Livingston cluster. Also make sure that the detector is in lock for whichever start and end times you choose as the inputs. To see where lock stretches are, look at the summary pages (e.g. https://ldas-jobs.ligo.caltech.edu/~alexander.urban/O2/calibration/C02/L1/day/20170808/).
diff --git a/gstlal-calibration/tests/check_calibration/pcal2darm_timeseries.py b/gstlal-calibration/tests/check_calibration/pcal2darm_timeseries.py
new file mode 100644
index 0000000000000000000000000000000000000000..1cf0572b9fec74a294a3321618d2d01431340ca4
--- /dev/null
+++ b/gstlal-calibration/tests/check_calibration/pcal2darm_timeseries.py
@@ -0,0 +1,199 @@
+#!/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 sys
+import os
+import numpy
+import time
+import resource
+
+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 gstlal import pipeparts
+from gstlal import calibration_parts
+from gstlal import test_common
+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 glue import segments
+
+parser = OptionParser()
+parser.add_option("--ifo", metavar = "name", help = "Name of the interferometer (IFO), e.g., H1, L1")
+parser.add_option("--raw-frame-cache", metavar = "name", help = "Raw frame cache file")
+parser.add_option("--calibrated-frame-cache", metavar = "name", help = "Calibrated frame cache file")
+parser.add_option("--config-file", metavar = "name", help = "Configurations file used to produce GDS/DCS calibrated frames, needed to get pcal line frequencies and correction factors")
+parser.add_option("--pcal-channel-name", metavar = "name", default = "CAL-PCALY_TX_PD_OUT_DQ", help = "Name of the pcal channel you wish to use")
+parser.add_option("--calibrated-channel-list", metavar = "list", default = None, help = "Comma-separated list of calibrated channels to compare to pcal")
+parser.add_option("--calcs-channel-list", metavar = "list", default = None, help = "Comma-separated list of calibrated channels in the raw frames to compare to pcal")
+
+options, filenames = parser.parse_args()
+
+# Read the config file
+def ConfigSectionMap(section):
+ dict1 = {}
+ options = Config.options(section)
+ for option in options:
+ try:
+ dict1[option] = Config.get(section, option)
+ if dict1[option] == -1:
+ DebugPrint("skip: %s" % option)
+ except:
+ print("exception on %s!" % option)
+ dict[option] = None
+ return dict1
+
+Config = ConfigParser.ConfigParser()
+Config.read(options.config_file)
+
+InputConfigs = ConfigSectionMap("InputConfigurations")
+filters = numpy.load(InputConfigs["filtersfilename"])
+
+ifo = options.ifo
+
+# Set up channel lists
+channel_list = []
+calibrated_channel_list = []
+calcs_channel_list = []
+channel_list.append((ifo, options.pcal_channel_name))
+if options.calibrated_channel_list is not None:
+ calibrated_channels = options.calibrated_channel_list.split(',')
+ for channel in calibrated_channels:
+ channel_list.append((ifo, channel))
+else:
+ calibrated_channels = []
+if options.calcs_channel_list is not None:
+ calcs_channels = options.calcs_channel_list.split(',')
+ for channel in calcs_channels:
+ channel_list.append((ifo, channel))
+else:
+ calcs_channels = []
+
+# Read stuff we need from the filters file
+frequencies = [float(filters['ka_pcal_line_freq']), float(filters['kc_pcal_line_freq']), float(filters['high_pcal_line_freq'])]
+pcal_corrections= [float(filters['ka_pcal_corr_re']), float(filters['ka_pcal_corr_im']), float(filters['kc_pcal_corr_re']), float(filters['kc_pcal_corr_im']), float(filters['high_pcal_corr_re']), float(filters['high_pcal_corr_im'])]
+demodulated_pcal_list = []
+try:
+ arm_length = float(filters['arm_length'])
+except:
+ arm_length = 3995.1
+
+# demodulation and averaging parameters
+filter_time = 20
+average_time = 1
+rate_out = 1
+
+#
+# =============================================================================
+#
+# Pipelines
+#
+# =============================================================================
+#
+
+
+def pcal2darm(pipeline, name):
+
+ # Get pcal from the raw frames
+ raw_data = pipeparts.mklalcachesrc(pipeline, location = options.raw_frame_cache, cache_dsc_regex = ifo)
+ raw_data = pipeparts.mkframecppchanneldemux(pipeline, raw_data, do_file_checksum = False, skip_bad_files = True, channel_list = map("%s:%s".__mod__, channel_list))
+ pcal = calibration_parts.hook_up(pipeline, raw_data, options.pcal_channel_name, ifo, 1.0)
+ pcal = calibration_parts.caps_and_progress(pipeline, pcal, "audio/x-raw,format=F64LE,channels=1,channel-mask=(bitmask)0x0", "pcal")
+ pcal = pipeparts.mktee(pipeline, pcal)
+
+ # Demodulate the pcal channel at the lines of interest
+ for i in range(0, len(frequencies)):
+ demodulated_pcal = calibration_parts.demodulate(pipeline, pcal, frequencies[i], True, rate_out, filter_time, 0.5, prefactor_real = pcal_corrections[2 * i], prefactor_imag = pcal_corrections[2 * i + 1])
+ demodulated_pcal_list.append(pipeparts.mktee(pipeline, demodulated_pcal))
+
+ # Check if we are taking pcal-to-darm ratios for CALCS data
+ for channel in calcs_channels:
+ calcs_deltal = calibration_parts.hook_up(pipeline, raw_data, channel, ifo, 1.0)
+ calcs_deltal = calibration_parts.caps_and_progress(pipeline, calcs_deltal, "audio/x-raw,format=F64LE,channels=1,channel-mask=(bitmask)0x0", channel)
+ calcs_deltal = pipeparts.mktee(pipeline, calcs_deltal)
+ for i in range(0, len(frequencies)):
+ # Demodulate DELTAL_EXTERNAL at each line
+ demodulated_calcs_deltal = calibration_parts.demodulate(pipeline, calcs_deltal, frequencies[i], True, rate_out, filter_time, 0.5)
+ # Take ratio \DeltaL(f) / pcal(f)
+ deltaL_over_pcal = calibration_parts.complex_division(pipeline, demodulated_calcs_deltaL, demodulated_pcal_list[i])
+ # Take a running average
+ deltaL_over_pcal = pipeparts.mkgeneric(pipeline, deltaL_over_pcal, "lal_smoothkappas", array_size = 1, avg_array_size = int(rate_out * average_time))
+ # Write to file
+ pipeparts.mknxydumpsink(pipeline, deltaL_over_pcal, "%s_%s_over_%s_at_line%d.txt" % (ifo, channel, options.pcal_channel_name, i + 1))
+
+ # Check if we are taking pcal-to-darm ratios for gstlal calibrated data
+ if options.calibrated_channel_list is not None:
+ # Get calibrated channels from the calibrated frames
+ hoft_data = pipeparts.mklalcachesrc(pipeline, location = options.calibrated_frame_cache, cache_dsc_regex = ifo)
+ hoft_data = pipeparts.mkframecppchanneldemux(pipeline, hoft_data, do_file_checksum = False, skip_bad_files = True, channel_list = map("%s:%s".__mod__, channel_list))
+ for channel in calibrated_channels:
+ hoft = calibration_parts.hook_up(pipeline, hoft_data, channel, ifo, 1.0)
+ hoft = calibration_parts.caps_and_progress(pipeline, hoft, "audio/x-raw,format=F64LE,channels=1,channel-mask=(bitmask)0x0", channel)
+ deltal = pipeparts.mkaudioamplify(pipeline, hoft, arm_length)
+ deltal = pipeparts.mktee(pipeline, deltal)
+ for i in range(0, len(frequencies)):
+ # Demodulate \DeltaL at each line
+ demodulated_deltal = calibration_parts.demodulate(pipeline, deltal, frequencies[i], True, rate_out, filter_time, 0.5)
+ # Take ratio \DeltaL(f) / pcal(f)
+ deltaL_over_pcal = calibration_parts.complex_division(pipeline, demodulated_deltal, demodulated_pcal_list[i])
+ # Take a running average
+ deltaL_over_pcal = pipeparts.mkgeneric(pipeline, deltaL_over_pcal, "lal_smoothkappas", array_size = 1, avg_array_size = int(rate_out * average_time))
+ # Write to file
+ pipeparts.mknxydumpsink(pipeline, deltaL_over_pcal, "%s_%s_over_%s_at_line%d.txt" % (ifo, channel, options.pcal_channel_name, i + 1))
+
+ #
+ # done
+ #
+
+ return pipeline
+
+#
+# =============================================================================
+#
+# Main
+#
+# =============================================================================
+#
+
+
+test_common.build_and_run(pcal2darm, "pcal2darm")
+
+
diff --git a/gstlal-calibration/tests/check_calibration/response_function.py b/gstlal-calibration/tests/check_calibration/response_function.py
new file mode 100644
index 0000000000000000000000000000000000000000..8a4aaddb5e0b31a0696c6922781887cee8304f95
--- /dev/null
+++ b/gstlal-calibration/tests/check_calibration/response_function.py
@@ -0,0 +1,171 @@
+##########################################################
+# Create a Bode plot of the response function as derived #
+# from different h(t) pipelines #
+##########################################################
+
+import matplotlib as mpl; mpl.use('Agg')
+from gwpy.plotter import BodePlot
+from gwpy.timeseries import TimeSeriesDict
+from gwpy.timeseries import TimeSeries
+from gwpy.frequencyseries import FrequencySeries
+import numpy
+import ConfigParser
+import sys
+from optparse import OptionParser, Option
+
+parser = OptionParser()
+
+parser.add_option("--gps_start_time", metavar = "seconds", help = "Set the GPS start time.")
+parser.add_option("--gps_end_time", metavar = "seconds", help = "Set the GPS end time.")
+parser.add_option("--dt", metavar = "seconds", help = "Sampling time interval.")
+parser.add_option("--ifo", metavar = "name", help = "Name of the IFO")
+parser.add_option("--c00_hoft_frames_cache", metavar = "name", help = "_")
+parser.add_option("--c01_hoft_frames_cache", metavar = "name", help = "_")
+parser.add_option("--raw_frames_cache", metavar = "name", help = "_")
+parser.add_option("--darm_err_channel_name", metavar = "name", help = "_")
+parser.add_option("--c00_hoft_channel_name", metavar = "name", help = "_")
+parser.add_option("--c01_hoft_channel_name", metavar = "name", help = "_")
+parser.add_option("--calcs_hoft_channel_name", metavar = "name", help = "_")
+parser.add_option("--response_file", metavar = "name", help = "_")
+
+options, filenames = parser.parse_args()
+
+start = int(options.gps_start_time)
+end = int(options.gps_end_time)
+dt = float(options.dt)
+ifo = options.ifo
+C00_hoft_frames_cache = options.c00_hoft_frames_cache
+C01_hoft_frames_cache = options.c01_hoft_frames_cache
+raw_frames_cache = options.raw_frames_cache
+DARM_ERR_channel_name = options.darm_err_channel_name
+C00_hoft_channel_name = options.c00_hoft_channel_name
+C01_hoft_channel_name = options.c01_hoft_channel_name
+CALCS_hoft_channel_name = options.calcs_hoft_channel_name
+response_file = numpy.load(options.response_file)
+
+#response_real = response_file['response_real_function']
+response_real = response_file['response_function'][1]
+#response_imag = response_file['response_imaginary_function']
+response_imag = response_file['response_function'][2]
+# I had to save the real and imaginary parts separately (there might be a way to save the complex number
+# into an npz file from Matlab, but I couldn't figure that out)
+response = response_real + 1j*response_imag
+# This is the same frequency vector as we'll use below for the h(t) data
+#freqs = response_file['frequency_vector']
+freqs = response_file['response_function'][0]
+f0 = freqs[0]
+df = freqs[1]-freqs[0]
+
+response[0] = 0 # zero out the DC component
+response = numpy.ndarray.flatten(response) # flatten to a 1D array
+
+response_fs = FrequencySeries(response, f0=f0, df=df) # packagae as a FrequecySeries object
+
+# Read in DARM_ERR data and compensate for the model jump delay be reading in data from start+dt to end+dt
+#DARM_ERR_data = TimeSeriesDict.read(raw_frames_cache, channels = ["%s:%s" % (ifo, DARM_ERR_channel_name)], start = start+dt, end = end+dt)
+DARM_ERR_data = TimeSeries.read(raw_frames_cache, "%s:%s" % (ifo, DARM_ERR_channel_name), start = start+dt, end = end+dt)
+#DARM_ERR_data = TimeSeriesDict.read(raw_frames_cache, "%s:%s" % (ifo, DARM_ERR_channel_name), start+dt, end+dt)
+# Read in CALCS data
+#CALCS_data = TimeSeriesDict.read(raw_frames_cache, channels = ["%s:%s" % (ifo, CALCS_hoft_channel_name)], start = start, end = end)
+CALCS_data = TimeSeries.read(raw_frames_cache, "%s:%s" % (ifo, CALCS_hoft_channel_name), start = start, end = end)
+# Read in the C00 data without kappas applied
+#C00_data = TimeSeriesDict.read(C00_hoft_frames_cache, channels = ["%s:%s" % (ifo, C00_hoft_channel_name)], start = start, end = end)
+C00_data = TimeSeries.read(C00_hoft_frames_cache, "%s:%s" % (ifo, C00_hoft_channel_name), start = start, end = end)
+# Read in the C01 data without kappas applied
+#C01_data = TimeSeriesDict.read(C01_hoft_frames_cache, channels = ["%s:%s" % (ifo, C01_hoft_channel_name)], start = start, end = end)
+C01_data = TimeSeries.read(C01_hoft_frames_cache, "%s:%s" % (ifo, C01_hoft_channel_name), start = start, end = end)
+
+# Pick out channel from dictionary
+#DARM_ERR_data = DARM_ERR_data["%s:%s" % (ifo, DARM_ERR_channel_name)]
+#CALCS_data = CALCS_data["%s:%s" % (ifo, CALCS_hoft_channel_name)]
+#DARM_ERR_data = DARM_ERR_data["%s:%s" % (ifo, DARM_ERR_channel_name)]
+
+# We want \Delta L, not strain, to find the response function
+#C00_data = C00_data["%s:%s" % (ifo, C00_hoft_channel_name)] * 3995.1
+#C01_data = C01_data["%s:%s" % (ifo, C01_hoft_channel_name)] * 3995.1
+
+dur = end - start
+averaging_time = 16
+chunk_start = start
+chunk_end = start + averaging_time
+
+CALCS_tf_data = numpy.zeros(len(response)) + 1j*numpy.zeros(len(response))
+C00_tf_data = numpy.zeros(len(response)) + 1j*numpy.zeros(len(response))
+C01_tf_data = numpy.zeros(len(response)) + 1j*numpy.zeros(len(response))
+
+N = 0
+
+while (chunk_end <= end):
+ # Correct for model jump delay in darm_err
+ DARM_ERR_chunk = DARM_ERR_data.crop(chunk_start+dt, chunk_end+dt, True)
+ DARM_ERR_chunk = DARM_ERR_chunk.detrend()
+ DARM_ERR_chunk_fft = DARM_ERR_chunk.average_fft(4, 2, window = 'hann')
+
+ CALCS_chunk = CALCS_data.crop(chunk_start, chunk_end, True)
+ CALCS_chunk = CALCS_chunk.detrend()
+ CALCS_chunk_fft = CALCS_chunk.average_fft(4, 2, window = 'hann')
+ CALCS_chunk_fft = CALCS_chunk_fft.filter([30]*6, [0.3]*6, 1e-12)
+
+ C00_chunk = C00_data.crop(chunk_start, chunk_end, True)
+ C00_chunk = C00_chunk.detrend()
+ C00_chunk_fft = C00_chunk.average_fft(4, 2, window = 'hann')
+
+ C01_chunk = C01_data.crop(chunk_start, chunk_end, True)
+ C01_chunk = C01_chunk.detrend()
+ C01_chunk_fft = C01_chunk.average_fft(4, 2, window = 'hann')
+
+ CALCS_chunk_tf = CALCS_chunk_fft / DARM_ERR_chunk_fft
+ C00_chunk_tf = C00_chunk_fft / DARM_ERR_chunk_fft
+ C01_chunk_tf = C01_chunk_fft / DARM_ERR_chunk_fft
+
+ CALCS_tf_data += CALCS_chunk_tf.value
+ C00_tf_data += C00_chunk_tf.value
+ C01_tf_data += C01_chunk_tf.value
+
+ chunk_start += averaging_time
+ chunk_end += averaging_time
+ N += 1
+
+
+CALCS_tf_data = CALCS_tf_data / N
+CALCS_tf = FrequencySeries(CALCS_tf_data, f0=f0, df=df)
+C00_tf_data = C00_tf_data / N
+C00_tf = FrequencySeries(C00_tf_data, f0=f0, df=df)
+C01_tf_data = C01_tf_data / N
+C01_tf = FrequencySeries(C01_tf_data, f0=f0, df=df)
+
+# Make plot that compares all of the derived response functions and the model response
+plot = BodePlot(response_fs, frequencies=freqs, dB = False, linewidth=2)
+plot.add_frequencyseries(CALCS_tf, dB=False, color='#4ba6ff', linewidth=2)
+plot.add_frequencyseries(C00_tf*3995.1, dB = False, color='#ee0000',linewidth=2)
+plot.add_frequencyseries(C01_tf*3995.1, dB = False, color="#94ded7", linewidth=2)
+plot.add_legend([r'Reference model response function', r'Front-end response function', r'Low-latency \texttt{gstlal} response function', r'High-latency \texttt{gstlal} response function'], loc='upper right', fontsize='x-small')
+plot.maxes.set_yscale('log')
+plot.paxes.set_yscale("linear")
+plot.save('%s_all_tf.pdf' % ifo)
+
+# Make a plot that compares the ratios of each derived response to the model
+ratio_CALCS = CALCS_tf / response_fs
+ratio_C00 = C00_tf / response_fs
+ratio_C01 = C01_tf / response_fs
+plot = BodePlot(ratio_CALCS, frequencies = freqs, dB = False, color='#4ba6ff', linewidth=2)
+plot.add_frequencyseries(ratio_C00*3995.1, dB = False, color='#ee0000',linewidth=2)
+plot.add_frequencyseries(ratio_C01*3995.1, dB = False, color="#94ded7", linewidth=2)
+plot.add_legend([r'Front-end response / Reference model response', r'Low-latency \texttt{gstlal} response / Reference model response', r'High-latency \texttt{gstlal} response / Reference model response'], loc='upper right', fontsize='small')
+plot.maxes.set_yscale('linear')
+plot.paxes.set_yscale('linear')
+plot.save('%s_all_tf_ratio.pdf' % ifo)
+
+plot = BodePlot(ratio_CALCS, frequencies = freqs, dB = False, color='#4ba6ff', linewidth=2)
+plot.add_frequencyseries(ratio_C00*3995.1, dB = False, color='#ee0000',linewidth=2)
+plot.add_frequencyseries(ratio_C01*3995.1, dB = False, color="#94ded7", linewidth=2)
+plot.add_legend([r'Front-end response / Reference model response', r'Low-latency \texttt{gstlal} response / Reference model response', r'High-latency \texttt{gstlal} response / Reference model response'], loc='upper right', fontsize='small')
+plot.maxes.set_yscale('linear')
+plot.paxes.set_yscale('linear')
+plot.maxes.set_ylim(.9, 1.1)
+plot.maxes.set_xlim(10, 5000)
+plot.paxes.set_ylim(-5, 5)
+plot.paxes.set_xlim(10, 5000)
+plot.save('%s_all_tf_ratio_zoomed.pdf' % ifo)
+
+
diff --git a/gstlal-calibration/tests/check_calibration/timeserieskappas.py b/gstlal-calibration/tests/check_calibration/timeserieskappas.py
new file mode 100644
index 0000000000000000000000000000000000000000..1a639c97c6bd160502caf1bf48e05306df16ba36
--- /dev/null
+++ b/gstlal-calibration/tests/check_calibration/timeserieskappas.py
@@ -0,0 +1,85 @@
+import matplotlib as mpl; mpl.use('Agg')
+from matplotlib import pyplot as mpl
+from gwpy.timeseries import TimeSeries
+from gwpy.segments import DataQualityFlag
+from optparse import OptionParser, Option
+
+parser = OptionParser()
+
+parser.add_option("--gps_start_time", metavar = "seconds", help = "Set the GPS start time.")
+parser.add_option("--gps_end_time", metavar = "seconds", help = "Set the GPS end time.")
+parser.add_option("--channel_name1", metavar = "name", help = "_")
+parser.add_option("--channel_name2", metavar = "name", help = "_")
+parser.add_option("--channel_name3", metavar = "name", help = "_")
+parser.add_option("--channel_name4", metavar = "name", help = "_")
+parser.add_option("--channel_name5", metavar = "name", help = "_")
+parser.add_option("--channel_name6", metavar = "name", help = "_")
+
+options, filenames = parser.parse_args()
+
+start = int(options.gps_start_time)
+end = int(options.gps_end_time)
+channel_name1 = str(options.channel_name1)
+channel_name2 = str(options.channel_name2)
+channel_name3 = str(options.channel_name3)
+channel_name4 = str(options.channel_name4)
+channel_name5 = str(options.channel_name5)
+channel_name6 = str(options.channel_name6)
+
+data1 = TimeSeries.get(channel_name1, start, end)
+data2 = TimeSeries.get(channel_name2, start, end)
+data3 = TimeSeries.get(channel_name3, start, end)
+data4 = TimeSeries.get(channel_name4, start, end)
+data5 = TimeSeries.get(channel_name5, start, end)
+data6 = TimeSeries.get(channel_name6, start, end)
+
+segs = DataQualityFlag.query('L1:DMT-CALIBRATED:1', start, end)
+
+plot1 = TimeSeries.plot(data1)
+mpl.ylabel('Correction value')
+title1 = channel_name1
+title1 = title1.replace('_', '\_')
+mpl.title(title1)
+plot1.add_state_segments(segs, plotargs=dict(label='Calibrated'))
+plot1.savefig('plot_'+channel_name1+'.png')
+
+plot2 = TimeSeries.plot(data2)
+mpl.ylabel('Correction value')
+title2 = channel_name2
+title2 = title2.replace('_', '\_')
+mpl.title(title2)
+plot2.add_state_segments(segs, plotargs=dict(label='Calibrated'))
+plot2.savefig('plot_'+channel_name2+'.png')
+
+plot3 = TimeSeries.plot(data3)
+mpl.ylabel('Correction value')
+title3 = channel_name3
+title3 = title3.replace('_', '\_')
+mpl.title(title3)
+plot3.add_state_segments(segs, plotargs=dict(label='Calibrated'))
+plot3.savefig('plot_'+channel_name3+'.png')
+
+plot4 = TimeSeries.plot(data4)
+mpl.ylabel('Correction value')
+title4 = channel_name4
+title4 = title4.replace('_', '\_')
+mpl.title(title4)
+plot4.add_state_segments(segs, plotargs=dict(label='Calibrated'))
+plot4.savefig('plot_'+channel_name4+'.png')
+
+plot5 = TimeSeries.plot(data5)
+mpl.ylabel('Correction value')
+title5 = channel_name5
+title5 = title5.replace('_', '\_')
+mpl.title(title5)
+plot5.add_state_segments(segs, plotargs=dict(label='Calibrated'))
+plot5.savefig('plot_'+channel_name5+'.png')
+
+plot6 = TimeSeries.plot(data6)
+mpl.ylabel('Correction value')
+title6 = channel_name2
+title6 = title6.replace('_', '\_')
+mpl.title(title6)
+plot6.add_state_segments(segs, plotargs=dict(label='Calibrated'))
+plot6.savefig('plot_'+channel_name6+'.png')
+