From a82ccc795180b67bf30e8be24a466a1ca28f8617 Mon Sep 17 00:00:00 2001 From: Aaron Viets <aaron.viets@ligo.org> Date: Tue, 30 Oct 2018 20:13:08 -0700 Subject: [PATCH] gstlal_compute_strain: Option to choose any bit of the ODC state vector to gate the noise subtraction. --- gstlal-calibration/bin/gstlal_compute_strain | 23 +++++++++++++------ ...CS_FreqIndepAndFccCorrections_Cleaning.ini | 5 +++- .../python/calibration_parts.py | 12 +++++----- .../tests/check_calibration/Makefile | 10 ++++---- 4 files changed, 31 insertions(+), 19 deletions(-) diff --git a/gstlal-calibration/bin/gstlal_compute_strain b/gstlal-calibration/bin/gstlal_compute_strain index 94a44eccab..0d1b471b0a 100755 --- a/gstlal-calibration/bin/gstlal_compute_strain +++ b/gstlal-calibration/bin/gstlal_compute_strain @@ -1814,7 +1814,7 @@ if compute_calib_statevector: # # Set the FILTERS-OK bit based on observation-ready transitions - filtersok = pipeparts.mkbitvectorgen(pipeline, obsintenttee, bit_vector=pow(2,3), threshold=2) + filtersok = pipeparts.mkbitvectorgen(pipeline, obsintenttee, bit_vector = pow(2,3), threshold=2) filtersok = pipeparts.mkcapsfilter(pipeline, filtersok, calibstate_caps) filtersok = calibration_parts.mkgate(pipeline, filtersok, obsreadytee, 4, attack_length = -int(filter_settle_time * calibstate_sr), hold_length = -int(filter_latency * calibstate_sr)) filtersok = pipeparts.mkbitvectorgen(pipeline, filtersok, bit_vector = pow(2,3), nongap_is_control = True) @@ -2169,12 +2169,21 @@ if remove_cal_lines: clean_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, straintee, calibration_lines)) +# Pick a bit from the ODC state vector to gate the power mains and noise subtraction +noisesub_gate_bitmask = 2 if config_version < 1 else int(Bitmasks["noisesubgatebitmask"]) +if compute_calib_statevector and (remove_power_lines or witness_channel_list is not None) and noisesub_gate_bitmask >= 0: + noisesubgate = pipeparts.mkgeneric(pipeline, odcstatevectortee, "lal_logicalundersample", required_on = noisesub_gate_bitmask, status_out = pow(2,28)) + noisesubgate = pipeparts.mkcapsfilter(pipeline, noisesubgate, calibstate_caps) + noisesubgatetee = pipeparts.mktee(pipeline, noisesubgate) +else: + noisesubgatetee = None + # Next, remove 60 Hz power lines and harmonics if remove_power_lines: if not remove_cal_lines: clean_strain = straintee powerlines = calibration_parts.caps_and_progress(pipeline, head_dict["powerlines"], "audio/x-raw, format=F64LE, channels=1, channel-mask=(bitmask)0x0", "powerlines") - clean_strain = calibration_parts.remove_harmonics_with_witness(pipeline, clean_strain, powerlines, 60, 5, powerlines_freq_var, filter_latency_factor, compute_rate = compute_factors_sr, rate_out = hoft_sr, num_avg = powerlines_tf_averaging_time * compute_factors_sr, obsready = None if not compute_calib_statevector else obsreadytee) + clean_strain = calibration_parts.remove_harmonics_with_witness(pipeline, clean_strain, powerlines, 60, 5, powerlines_freq_var, filter_latency_factor, compute_rate = compute_factors_sr, rate_out = hoft_sr, num_avg = powerlines_tf_averaging_time * compute_factors_sr, noisesub_gate_bit = noisesubgatetee) # Remove excess noise using any provided witness channels if witness_channel_list is not None: @@ -2190,10 +2199,6 @@ if witness_channel_list is not None: if not (remove_cal_lines or remove_power_lines): clean_strain = straintee - # If possible, gate the data being used to compute transfer functions to be sure we are locked - if not compute_calib_statevector: - obsreadytee = None - for i in range(0, len(witness_channel_list)): # Length of ffts used to compute FIR filters witness_fft_samples = int(witness_channel_fft_time * witness_rates[i]) @@ -2212,7 +2217,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, 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, delay_time = witness_delay_time, wait_time = witness_wait_time, critical_lock_loss_time = critical_lock_loss_time, filename = None if witness_tf_filename is None else "%s_%d.txt" % (witness_tf_filename, 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], noisesub_gate_bit = noisesubgatetee, delay_time = witness_delay_time, wait_time = witness_wait_time, critical_lock_loss_time = critical_lock_loss_time, filename = None if witness_tf_filename is None else "%s_%d.txt" % (witness_tf_filename, i)) witness_delay_time += witness_delay_increment witness_wait_time += witness_wait_increment @@ -2254,6 +2259,10 @@ if compute_calib_statevector and (remove_cal_lines or remove_power_lines or witn # Add these into the CALIB_STATE_VECTOR calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, calibstatevector, clean_hoft_ok_lowfreq, clean_hoft_ok_midfreq)) +# Check if we gated the cleaning with any ODC bits, and if so, add that to the CALIB_STATE_VECTOR +if noisesubgatetee is not None: + calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, calibstatevector, noisesubgatetee)) + # # Produce time-dependent correction factors to be recorded in the frames # diff --git a/gstlal-calibration/config_files/O2/H1/tests/H1DCS_FreqIndepAndFccCorrections_Cleaning.ini b/gstlal-calibration/config_files/O2/H1/tests/H1DCS_FreqIndepAndFccCorrections_Cleaning.ini index db82f6aac7..8102826caa 100644 --- a/gstlal-calibration/config_files/O2/H1/tests/H1DCS_FreqIndepAndFccCorrections_Cleaning.ini +++ b/gstlal-calibration/config_files/O2/H1/tests/H1DCS_FreqIndepAndFccCorrections_Cleaning.ini @@ -1,6 +1,6 @@ [InputConfigurations] # Track what "version" of config file this is, so that the pipeline knows which options are present in this file -ConfigVersion: 0 +ConfigVersion: 1 # Filters file containing calibration FIR filters, relative to the directory gstlal-calibration/tests/check_calibration/, from which the pipeline is expected to be run FiltersFileName: H1DCS_newsrcline_1173225472.npz # Data source should be set to frames or lvshm @@ -199,6 +199,8 @@ PCALChannel: CAL-PCALY_TX_PD_OUT_DQ # Coherence Uncertainty Channel Names # ####################################### CohUncSusLine1Channel: CAL-CS_TDEP_SUS_LINE1_UNCERTAINTY +CohUncSusLine2Channel: CAL-CS_TDEP_SUS_LINE2_UNCERTAINTY +CohUncSusLine3Channel: 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 @@ -296,6 +298,7 @@ CBCHWInjBitmask: 16777216 BurstHWInjBitmask: 33554432 DetCharHWInjBitmask: 67108864 StochHWInjBitmask: 8388608 +NoiseSubGateBitmask: 2 [PipelineConfigurations] BufferLength: 1.0 diff --git a/gstlal-calibration/python/calibration_parts.py b/gstlal-calibration/python/calibration_parts.py index e3d1713122..62af24a649 100644 --- a/gstlal-calibration/python/calibration_parts.py +++ b/gstlal-calibration/python/calibration_parts.py @@ -195,7 +195,7 @@ def remove_harmonics(pipeline, signal, f0, num_harmonics, f0_var, filter_latency return elem -def remove_harmonics_with_witness(pipeline, signal, witness, f0, num_harmonics, f0_var, filter_latency, compute_rate = 16, rate_out = 16384, num_avg = 2048, obsready = None): +def remove_harmonics_with_witness(pipeline, signal, witness, f0, num_harmonics, f0_var, filter_latency, compute_rate = 16, rate_out = 16384, num_avg = 2048, noisesub_gate_bit = None): # remove line(s) from a spectrum. filter length for demodulation (given in seconds) is adjustable # function argument caps must be complex caps @@ -256,8 +256,8 @@ def remove_harmonics_with_witness(pipeline, signal, witness, f0, num_harmonics, tf_at_f = pipeparts.mkgeneric(pipeline, tf_at_f, "lal_insertgap", chop_length = 1000000000) # Removing one second # Remove worthless data from computation of transfer function if we can - if obsready is not None: - tf_at_f = mkgate(pipeline, tf_at_f, obsready, 1, attack_length = -((1.0 - filter_latency) * filter_samples)) + if noisesub_gate_bit is not None: + tf_at_f = mkgate(pipeline, tf_at_f, noisesub_gate_bit, 1, attack_length = -((1.0 - filter_latency) * filter_samples)) tf_at_f = pipeparts.mkgeneric(pipeline, tf_at_f, "lal_smoothkappas", default_kappa_re = 0, array_size = 1, avg_array_size = num_avg, default_to_median = True, filter_latency = filter_latency) # Use gated, averaged transfer function to reconstruct the sinusoid as it appears in the signal from the witness channel @@ -920,7 +920,7 @@ def update_filters(filter_maker, arg, filter_taker, maker_prop_name, taker_prop_ firfilter = filter_maker.get_property(maker_prop_name)[filter_number][::-1] filter_taker.set_property(taker_prop_name, firfilter) -def clean_data(pipeline, signal, signal_rate, witnesses, witness_rate, fft_length, fft_overlap, num_ffts, min_ffts, update_samples, fir_length, frequency_resolution, filter_taper_length, use_median = False, notch_frequencies = [], obsready = None, delay_time = 0.0, wait_time = 0.0, critical_lock_loss_time = 0, filename = None): +def clean_data(pipeline, signal, signal_rate, witnesses, witness_rate, fft_length, fft_overlap, num_ffts, min_ffts, update_samples, fir_length, frequency_resolution, filter_taper_length, use_median = False, notch_frequencies = [], noisesub_gate_bit = None, delay_time = 0.0, wait_time = 0.0, critical_lock_loss_time = 0, filename = None): # # Use witness channels that monitor the environment to remove environmental noise @@ -940,8 +940,8 @@ def clean_data(pipeline, signal, signal_rate, witnesses, witness_rate, fft_lengt resampled_signal = mkresample(pipeline, signal_tee, 5, False, witness_rate) transfer_functions = mkinterleave(pipeline, numpy.insert(witness_tees, 0, resampled_signal, axis = 0)) - if obsready is not None: - transfer_functions = mkgate(pipeline, transfer_functions, obsready, 1) + if noisesub_gate_bit is not None: + transfer_functions = mkgate(pipeline, transfer_functions, noisesub_gate_bit, 1) transfer_functions = pipeparts.mkgeneric(pipeline, transfer_functions, "lal_transferfunction", fft_length = fft_length, fft_overlap = fft_overlap, num_ffts = num_ffts, min_ffts = min_ffts, update_samples = update_samples, make_fir_filters = -1, fir_length = fir_length, frequency_resolution = frequency_resolution, high_pass = 15, update_after_gap = True, use_median = use_median, notch_frequencies = notch_frequencies, use_first_after_gap = critical_lock_loss_time * witness_rate, update_delay_samples = delay_time * witness_rate, filename = filename) signal_minus_noise = [signal_tee] for i in range(0, len(witnesses)): diff --git a/gstlal-calibration/tests/check_calibration/Makefile b/gstlal-calibration/tests/check_calibration/Makefile index bb202116ab..99829fcabf 100644 --- a/gstlal-calibration/tests/check_calibration/Makefile +++ b/gstlal-calibration/tests/check_calibration/Makefile @@ -8,8 +8,8 @@ IFO = H # determines where to look for filters files (e.g., O1, O2, O3, ER10, ER13, ER14, PreER10, PreER13, PreER14) OBSRUN = O2 -START = 1186938402 -END = 1186946716 +START = 1186944000 +END = 1186945024 SHMRUNTIME = 400 # How much time does the calibration need to settle at the start and end? PLOT_WARMUP_TIME = 8000 @@ -22,7 +22,7 @@ GDSTESTCONFIGS = ../../config_files/O2/H1/tests/H1GDS_LowLatency_AllCorrections_ DCSTESTCONFIGS = ../../config_files/O2/H1/tests/H1DCS_AllCorrections_Cleaning.ini GDSSHMCONFIGS = ../../config_files/PreER13/H1/tests/H1GDS_1222058826_shm2frames.ini -all: pcal_DCS_transfer_functions +all: DCS_over_C02 ############################################### ### These commands should change less often ### @@ -92,8 +92,8 @@ latency_test: $(IFO)1_hoft_GDS_SHM_frames.cache GDS_over_C02: $(IFO)1_hoft_GDS_frames.cache $(IFO)1_C02_frames.cache python plot_transfer_function.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --numerator-frame-cache $(IFO)1_hoft_GDS_frames.cache --denominator-frame-cache $(IFO)1_C02_frames.cache --numerator-channel-name GDS-CALIB_STRAIN --denominator-channel-name DCS-CALIB_STRAIN_C02 --magnitude-min 0.7 --magnitude-max 1.3 --phase-min -20.0 --phase-max 20.0 --plot-title 'GDS Test / C02 Transfer Function' -DCS_over_C02: $(IFO)1_hoft_DCS_frames.cache $(IFO)1_C02_frames.cache - python plot_transfer_function.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --numerator-frame-cache $(IFO)1_hoft_DCS_frames.cache --denominator-frame-cache $(IFO)1_C02_frames.cache --numerator-channel-name DCS-CALIB_STRAIN --denominator-channel-name DCS-CALIB_STRAIN_C02 --magnitude-min 0.7 --magnitude-max 1.3 --phase-min -20.0 --phase-max 20.0 --plot-title 'DCS Test / C02 Transfer Function' +DCS_over_C02: $(IFO)1_hoft_DCS_FCC_frames.cache $(IFO)1_C02_frames.cache + python plot_transfer_function.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --numerator-frame-cache $(IFO)1_hoft_DCS_FCC_frames.cache --denominator-frame-cache $(IFO)1_C02_frames.cache --numerator-channel-name DCS-CALIB_STRAIN --denominator-channel-name DCS-CALIB_STRAIN_C02 --magnitude-min 0.7 --magnitude-max 1.3 --phase-min -20.0 --phase-max 20.0 --plot-title 'DCS Test / C02 Transfer Function' noise_subtraction_tf_DCS: $(IFO)1_hoft_DCS_frames.cache python plot_transfer_function.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --numerator-frame-cache $(IFO)1_hoft_DCS_frames.cache --denominator-frame-cache $(IFO)1_hoft_DCS_frames.cache --numerator-channel-name DCS-CALIB_STRAIN_CLEAN --denominator-channel-name DCS-CALIB_STRAIN --magnitude-min 0.0 --magnitude-max 1.5 --phase-min -20.0 --phase-max 20.0 --plot-title 'Noise Subtraction Transfer Function' -- GitLab