gstlal-calibration: New plotting script for BNS range.

gstlal-calibration/bin/gstlal_compute_strain

@@ -977,7 +977,7 @@ if ChannelNames["witnesschannellist"] != "None":
 if compute_calib_statevector:
 	noisesub_gate_channel = ChannelNames["noisesubgatechannel"] if "noisesubgatechannel" in ChannelNames else lownoise_channel_name 
 	noisesub_gate_bitmask = int(Bitmasks["noisesubgatebitmask"]) if "noisesubgatebitmask" in Bitmasks else lownoise_bitmask
-if compute_calib_statevector and (any(line_witness_channel_list) or any(witness_channel_list)) and noisesub_gate_channel != obsintent_channel_name and noisesub_gate_channel != lownoise_channel_name and noisesub_gate_channel != hwinj_channel_name and noisesub_gate_bitmask >= 0:
+if compute_calib_statevector and (any(line_witness_channel_list) or any(witness_channel_list)) and noisesub_gate_channel != obsintent_channel_name and noisesub_gate_channel != lownoise_channel_name and noisesub_gate_channel != hwinj_channel_name and noisesub_gate_bitmask > 0:
 	channel_list.append((instrument, noisesub_gate_channel))
 	headkeys.append("noisesubgatechannel")
 

gstlal-calibration/tests/check_calibration/ASD_plots

@@ -41,6 +41,25 @@ if options.analyze_additional_hoft_channel:
 # make asds
 calcs_asd = calcs_data.asd(4,2,method='lal_median')
 hoft_asd = hoft_data.asd(4,2,method='lal_median')
+
+numpy.savetxt("%s_CALIB_STRAIN_ASD_%d_%d%s.txt" % (options.ifo, start_time, end_time, options.filename_suffix), calcs_asd)
+numpy.savetxt("%s_CALIB_STRAIN_CLEAN_ASD_%d_%d%s.txt" % (options.ifo, start_time, end_time, options.filename_suffix), hoft_asd)
+
+calcs_asd_tmp = numpy.loadtxt("%s_CALIB_STRAIN_ASD_%d_%d%s.txt" % (options.ifo, start_time, end_time, options.filename_suffix))
+hoft_asd_tmp = numpy.loadtxt("%s_CALIB_STRAIN_CLEAN_ASD_%d_%d%s.txt" % (options.ifo, start_time, end_time, options.filename_suffix))
+
+calcs_asd_txt = []
+hoft_asd_txt = []
+freq = 0.0
+for i in range(0, len(calcs_asd)):
+	calcs_asd_txt.append([freq, calcs_asd_tmp[i]])
+	hoft_asd_txt.append([freq, hoft_asd_tmp[i]])
+	freq = freq + 0.25
+calcs_asd_txt = numpy.array(calcs_asd_txt)
+hoft_asd_txt = numpy.array(hoft_asd_txt)
+numpy.savetxt("%s_CALIB_STRAIN_ASD_%d_%d%s.txt" % (options.ifo, start_time, end_time, options.filename_suffix), calcs_asd_txt)
+numpy.savetxt("%s_CALIB_STRAIN_CLEAN_ASD_%d_%d%s.txt" % (options.ifo, start_time, end_time, options.filename_suffix), hoft_asd_txt)
+
 if options.analyze_additional_hoft_channel:
 	additional_hoft_asd = additional_hoft_data.asd(4,2)
 

gstlal-calibration/tests/check_calibration/Makefile

@@ -6,22 +6,24 @@
 # which interferometer (H or L)
 IFO = H
 # determines where to look for filters files (e.g., O1, O2, O3, ER10, ER13, ER14, PreER10, PreER13, PreER14)
-OBSRUN = ER13
+OBSRUN = ER14
 
-START = $(shell echo 1232874910 - 1715 | bc)
+START = $(shell echo 1234694144 - 715 | bc)
 #1229094912
 #1225967424
 #1185763328
-END = $(shell echo 1232875986 + 1715 | bc)
+END = $(shell echo 1234694144 + 4096 + 715 | bc)
 #1229099008
 #1225968448
 #1185771520
 SHMRUNTIME = 400
 # How much time does the calibration need to settle at the start and end?
-PLOT_WARMUP_TIME = 1725
-PLOT_COOLDOWN_TIME = 1777
+PLOT_WARMUP_TIME = 715
+#885
+PLOT_COOLDOWN_TIME = 715
+#4521
 
-GDSCONFIGS = Filters/ER13/GDSFilters/H1GDS_noisesub_test_1232874910.ini
+GDSCONFIGS = Filters/ER14/GDSFilters/H1GDS_noisesub_test_1232874910.ini
 DCSCONFIGS = Filters/ER13/GDSFilters/L1DCS_TEST_1231620000.ini
 #../../config_files/O2/H1/tests/H1DCS_AllCorrections_Cleaning.ini
 DCSFCCCONFIGS = ../../config_files/O2/H1/tests/H1DCS_FreqIndepAndFccCorrections_Cleaning.ini
@@ -29,7 +31,8 @@ GDSTESTCONFIGS = ../../config_files/PreER13/H1/H1GDS_TEST_1225558818.ini
 DCSTESTCONFIGS = ../../config_files/O2/H1/tests/H1DCS_AllCorrections_Cleaning_TEST.ini
 GDSSHMCONFIGS = Filters/ER14/GDSFilters/L1GDS_1234630818_latency_test.ini
 
-all: noise_subtraction_ASD_GDS noise_subtraction_tf_GDS filters_tf_GDS
+all: noise_subtraction_range_plots_GDS
+#noise_subtraction_ASD_GDS noise_subtraction_range_plots_GDS
 
 ###############################################
 ### These commands should change less often ###
@@ -143,6 +146,9 @@ noise_subtraction_ASD_DCH_DCS: $(IFO)1_hoft_DCS_frames.cache $(IFO)1_clean_C02_f
 noise_subtraction_ASD_DCS_TEST: $(IFO)1_hoft_DCS_TEST_frames.cache
 	./ASD_comparison_plots --ifo $(IFO)1 --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --raw-frame-cache $(IFO)1_hoft_DCS_TEST_frames.cache --calcs-channel-name DCS-CALIB_STRAIN --hoft-frame-cache $(IFO)1_hoft_DCS_TEST_frames.cache --hoft-channel-name DCS-CALIB_STRAIN_CLEAN
 
+noise_subtraction_ASD_C00: $(IFO)1_C00_frames.cache
+	./ASD_plots --ifo $(IFO)1 --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --raw-frame-cache $(IFO)1_C00_frames.cache --calcs-channel-name GDS-CALIB_STRAIN --hoft-frame-cache $(IFO)1_C00_frames.cache --hoft-channel-name GDS-CALIB_STRAIN_CLEAN
+
 noise_subtraction_tf_GDS: $(IFO)1_hoft_GDS_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_hoft_GDS_frames.cache --numerator-channel-list GDS-CALIB_STRAIN_CLEAN --denominator-channel-name GDS-CALIB_STRAIN --magnitude-min 0.0 --magnitude-max 1.5 --phase-min -20.0 --phase-max 20.0 --numerator-name 'CLEAN' --denominator-name 'STRAIN' --labels 'CALIB_STRAIN_CLEAN / CALIB_STRAIN' --use-median
 
@@ -155,6 +161,9 @@ noise_subtraction_tf_DCH_DCS: $(IFO)1_hoft_DCS_frames.cache $(IFO)1_clean_C02_fr
 noise_subtraction_tf_DCS_DCH: $(IFO)1_hoft_DCS_frames.cache $(IFO)1_clean_C02_frames.cache
 	python plot_transfer_function.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --numerator-frame-cache $(IFO)1_clean_C02_frames.cache --denominator-frame-cache $(IFO)1_hoft_DCS_frames.cache --numerator-channel-list DCH-CLEAN_STRAIN_C02 --denominator-channel-name DCS-CALIB_STRAIN_CLEAN --magnitude-min 0.7 --magnitude-max 1.3 --phase-min -20.0 --phase-max 20.0 --numerator-name 'DCH' --denominator-name 'DCS' --labels 'DCH_CLEAN / DCS_CLEAN' --use-median
 
+noise_subtraction_range_plots_GDS: $(IFO)1_hoft_GDS_frames.cache
+	python plot_BNS_range.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --frame-cache-list '$(IFO)1_hoft_GDS_frames.cache,$(IFO)1_hoft_GDS_frames.cache' --channel-list 'GDS-CALIB_STRAIN,GDS-CALIB_STRAIN_CLEAN' --range-min 0.0 --range-max 140
+
 filters:
 	if [ -d Filters/$(OBSRUN)/GDSFilters ]; then \
 		svn up Filters/$(OBSRUN)/GDSFilters; \

gstlal-calibration/tests/check_calibration/frame_manipulator.py

@@ -299,6 +299,8 @@ channel_list.append("GRD-IFO_OK")
 #channel_list.append("ISC_LOCK_STATUS")
 channel_list.append("GRD-ISC_LOCK_OK")
 channel_list.append("GRD-ISC_LOCK_ERROR")
+channel_list.append("GRD-IFO_INTENT")
+channel_list.append("GRD-IFO_READY")
 temp_list = channel_list
 channel_list = []
 for chan in temp_list:

gstlal-calibration/tests/check_calibration/plot_BNS_range.py

+#!/usr/bin/env python
+# Copyright (C) 2019  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
+from math import pi
+import resource
+import datetime
+import time
+
+from gwpy.timeseries import TimeSeries
+from gwpy.astro import inspiral_range
+
+import matplotlib
+matplotlib.rcParams['font.family'] = 'Times New Roman'
+matplotlib.rcParams['font.size'] = 22
+matplotlib.rcParams['legend.fontsize'] = 18
+matplotlib.rcParams['mathtext.default'] = 'regular'
+matplotlib.use('Agg')
+import glob
+import matplotlib.pyplot as plt
+
+from optparse import OptionParser, Option
+import ConfigParser
+
+
+parser = OptionParser()
+parser.add_option("--gps-start-time", metavar = "seconds", type = int, help = "GPS time at which to start processing data")
+parser.add_option("--gps-end-time", metavar = "seconds", type = int, help = "GPS time at which to stop processing data")
+parser.add_option("--ifo", metavar = "name", type = str, help = "Name of the interferometer (IFO), e.g., H1, L1")
+parser.add_option("--frame-cache-list", metavar = "list", type = str, help = "Comma-separated list of frame cache files containing strain data")
+parser.add_option("--channel-list", metavar = "list", type = str, help = "Comma-separated list of channel names for strain data")
+parser.add_option("--integration-time", metavar = "seconds", type = int, default = 64, help = "Amount of data in seconds to use for each calculation of range")
+parser.add_option("--stride", metavar = "seconds", type = int, default = 32, help = "Time-separation in seconds between consecutive values of range")
+parser.add_option("--range-min", metavar = "Mpc", type = float, default = 0.0, help = "Minimum value for range on plot")
+parser.add_option("--range-max", metavar = "Mpc", type = float, default = 160.0, help = "Maximum value for range on plot")
+
+options, filenames = parser.parse_args()
+
+# parameters
+gps_start_time = options.gps_start_time
+gps_end_time = options.gps_end_time
+ifo = options.ifo
+integration_time = options.integration_time
+stride = options.stride
+range_min = options.range_min
+range_max = options.range_max
+
+num_points = int((gps_end_time - gps_start_time - integration_time + stride) / stride)
+
+# Re-format list of frame caches and channels
+frame_cache_list = options.frame_cache_list.split(',')
+channel_list = options.channel_list.split(',')
+if len(frame_cache_list) != len(channel_list):
+	raise ValueError("--frame-cache-list and --channel-list must be the same length.")
+
+# Make a time vector
+times = []
+for i in range(0, num_points):
+	times.append(i * stride + 0.5 * integration_time)
+
+# Decide what unit ot time to use in the plot
+dur = times[-1]
+if dur > 60 * 60 * 100:
+	t_unit = 'days'
+	sec_per_t_unit = 60.0 * 60.0 * 24.0
+elif dur > 60 * 100:
+	t_unit = 'hours'
+	sec_per_t_unit = 60.0 * 60.0
+elif dur > 100:
+	t_unit = 'minutes'
+	sec_per_t_unit = 60.0
+else:
+	t_unit = 'seconds'
+	sec_per_t_unit = 1.0
+
+for i in range(0, num_points):
+	times[i] = times[i] / sec_per_t_unit
+
+# Collect range data in arrays
+ranges = []
+medians = []
+stds = []
+for i in range(0, len(channel_list)):
+	ranges.append([])
+	for j in range(0, num_points):
+		data = TimeSeries.read(frame_cache_list[i], "%s:%s" % (ifo, channel_list[i]), start = gps_start_time + j * stride, end = gps_start_time + j * stride + integration_time)
+		PSD = data.psd(8, 4, method = 'lal_median')
+		BNS_range = float(inspiral_range(PSD, fmin=10).value)
+		ranges[i].append(BNS_range)
+	medians.append(numpy.median(ranges[i]))
+	stds.append(numpy.std(ranges[i]))
+
+# Make plots
+colors = ['b', 'y', 'r', 'c', 'm', 'g'] # Hopefully the user will not want to plot more than six datasets on one plot.
+plt.figure(figsize = (15, 10))
+for i in range(0, len(channel_list)):
+	plt.plot(times, ranges[i], colors[i % 6], linewidth = 0.5, label = r'%s [median = %0.1f Mpc, $\sigma$ = %0.1f Mpc]' % (channel_list[i].replace('_', '\_'), medians[i], stds[i]))
+	plt.title("%s binary neutron star inspiral range" % ifo)
+	plt.ylabel('Angle-averaged range [Mpc]')
+	plt.xlabel('Time [%s] from %s UTC (%d)' % (t_unit, time.strftime("%b %d %Y %H:%M:%S", time.gmtime(gps_start_time + 315964782)), gps_start_time))
+	plt.ylim(range_min, range_max)
+	plt.grid(True)
+	leg = plt.legend(fancybox = True)
+	leg.get_frame().set_alpha(0.5)
+plt.savefig('%s_BNS_range_%d-%d.png' % (ifo, int(gps_start_time), int(dur)))
+plt.savefig('%s_BNS_range_%d-%d.pdf' % (ifo, int(gps_start_time), int(dur)))
+