Skip to content
Snippets Groups Projects
Commit e0ef8342 authored by Aaron Viets's avatar Aaron Viets
Browse files

lal_calibcorr_test.py: improvements to plotting and time delays

parent 65170262
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
tests/asd_derr 0 → 100644
+ 32769
− 0
View file @ e0ef8342
This diff is collapsed.
gds-dcs-filter-generation @ 018e42e5
Subproject commit c0b4e6d26f69994e0c5019c5fa3fc71a5ee3d984
Subproject commit 018e42e57c9e729058420376ceebefbc4fc33090
tests/lal_calibcorr_test.py
+ 84
− 86
View file @ e0ef8342
......@@ -49,7 +49,7 @@ from lal import LIGOTimeGPS
from ligo import segments
gps_start_time = 1269120000
gps_end_time = 1269130000
gps_end_time = 1269127000
chop_time = 5000
filters_name = "H1GDS_1324762762.npz"
......@@ -144,6 +144,14 @@ def simulate_error(df, fmax, amplitude):
return mag * np.exp(1j * phase)
# Get the filters
filters = np.load("/home/aaron.viets/src/gstlal/gstlal-calibration/tests/check_calibration/gds-dcs-filter-generation/filters/PreER15/GDSFilters/H1GDS_1324762762.npz")
# Errors in the models
C_error = np.ones(32769) #simulate_error(0.25, 8192, 0.02)
tst_error = np.ones(32769) #simulate_error(0.25, 8192, 0.02)
pum_error = np.ones(32769) #simulate_error(0.25, 8192, 0.02)
uim_error = np.ones(32769) #simulate_error(0.25, 8192, 0.02)
#
# Read real data from raw frames and calibrate it with real GDS calibration filters.
......@@ -155,8 +163,7 @@ def simulate_error(df, fmax, amplitude):
def lal_calibcorr_01(pipeline, name):
# Get the filters and stuff we need from them.
filters = np.load("/home/aaron.viets/src/gstlal/gstlal-calibration/tests/check_calibration/gds-dcs-filter-generation/filters/PreER15/GDSFilters/H1GDS_1324762762.npz")
# Get stuff we need from the filters.
if "res_highpass" in filters:
invsens_highpass = filters["res_highpass"]
invsens_highpass_delay = int(filters['res_highpass_delay'])
......@@ -236,11 +243,6 @@ def lal_calibcorr_01(pipeline, name):
uim_indices = np.asarray(uim_indices)
# Errors in the models
C_error = np.ones(32769) #simulate_error(0.25, 8192, 0.02)
tst_error = np.ones(32769) #simulate_error(0.25, 8192, 0.02)
pum_error = np.ones(32769) #simulate_error(0.25, 8192, 0.02)
uim_error = np.ones(32769) #simulate_error(0.25, 8192, 0.02)
C_error_calibcorr = []
tst_error_calibcorr = []
pum_error_calibcorr = []
......@@ -416,7 +418,7 @@ def lal_calibcorr_01(pipeline, name):
#
test_common.build_and_run(lal_calibcorr_01, "lal_calibcorr_01", segment = segments.segment((LIGOTimeGPS(0, 1000000000 * gps_start_time), LIGOTimeGPS(0, 1000000000 * gps_end_time))))
#test_common.build_and_run(lal_calibcorr_01, "lal_calibcorr_01", segment = segments.segment((LIGOTimeGPS(0, 1000000000 * gps_start_time), LIGOTimeGPS(0, 1000000000 * gps_end_time))))
# Compute Models including errors
......@@ -439,71 +441,68 @@ uim_with_error = np.array((invsens_model[0], np.real(uim_with_error), np.imag(ui
R_with_error = np.array((invsens_model[0], np.real(R_with_error), np.imag(R_with_error)))
# Identify any files that have filters transfer functions in them
tf_files = [f for f in os.listdir(options.tf_file_directory) if (os.path.isfile(f) and ('GDS' in f or 'DCS' in f) and 'npz' in f and 'filters_transfer_function' in f and '.txt' in f)]
os.system("rm *_ratio_magnitude.txt *_ratio_phase.txt *_reformatted.txt")
tf_files0 = [f for f in os.listdir('.') if (os.path.isfile(f) and ('GDS' in f or 'DCS' in f) and 'npz' in f and 'filters_transfer_function' in f and '.txt' in f)]
tf_files = []
j = 0
k = 0
indices = []
# Put the TF files in order
for i in range(0, len(tf_files)):
if '_filters_transfer_function_' in tf_files[i] and "_res_" in tf_files[i] and not ("calibcorr" in tf_files[i]):
response_file = tf_files.pop(i)
for i in range(0, len(tf_files0)):
if '_filters_transfer_function_' in tf_files0[i] and "_res_" in tf_files0[i] and not ("calibcorr" in tf_files0[i]):
response_file = tf_files0[i]
tf_files.append(response_file)
j += 1
indices.append(j)
for i in range(0, len(tf_files)):
if '_filters_transfer_function_' in tf_files[i] and "_res_" in tf_files[i] and "calibcorr" in tf_files[i]:
response_file = tf_files.pop(i)
k += 1
for i in range(0, len(tf_files0)):
if '_filters_transfer_function_' in tf_files0[i] and "_res_" in tf_files0[i] and "calibcorr" in tf_files0[i]:
response_file = tf_files0[i]
tf_files.append(response_file)
j += 1
indices.append(j)
for i in range(0, len(tf_files)):
if '_filters_transfer_function_' in tf_files[i] and "_tst_" in tf_files[i] and not ("calibcorr" in tf_files[i]):
response_file = tf_files.pop(i)
k += 1
indices.append(k)
for i in range(0, len(tf_files0)):
if '_filters_transfer_function_' in tf_files0[i] and "_tst_" in tf_files0[i] and not ("calibcorr" in tf_files0[i]):
response_file = tf_files0[i]
tf_files.append(response_file)
j += 1
indices.append(j)
for i in range(0, len(tf_files)):
if '_filters_transfer_function_' in tf_files[i] and "_tst_" in tf_files[i] and "calibcorr" in tf_files[i]:
response_file = tf_files.pop(i)
k += 1
for i in range(0, len(tf_files0)):
if '_filters_transfer_function_' in tf_files0[i] and "_tst_" in tf_files0[i] and "calibcorr" in tf_files0[i]:
response_file = tf_files0[i]
tf_files.append(response_file)
j += 1
indices.append(j)
for i in range(0, len(tf_files)):
if '_filters_transfer_function_' in tf_files[i] and "_pum_" in tf_files[i] and not ("calibcorr" in tf_files[i]):
response_file = tf_files.pop(i)
k += 1
indices.append(k)
for i in range(0, len(tf_files0)):
if '_filters_transfer_function_' in tf_files0[i] and "_pum_" in tf_files0[i] and not ("calibcorr" in tf_files0[i]):
response_file = tf_files0[i]
tf_files.append(response_file)
j += 1
indices.append(j)
for i in range(0, len(tf_files)):
if '_filters_transfer_function_' in tf_files[i] and "_pum_" in tf_files[i] and "calibcorr" in tf_files[i]:
response_file = tf_files.pop(i)
k += 1
for i in range(0, len(tf_files0)):
if '_filters_transfer_function_' in tf_files0[i] and "_pum_" in tf_files0[i] and "calibcorr" in tf_files0[i]:
response_file = tf_files0[i]
tf_files.append(response_file)
j += 1
indices.append(j)
for i in range(0, len(tf_files)):
if '_filters_transfer_function_' in tf_files[i] and "_uim_" in tf_files[i] and not ("calibcorr" in tf_files[i]):
response_file = tf_files.pop(i)
k += 1
indices.append(k)
for i in range(0, len(tf_files0)):
if '_filters_transfer_function_' in tf_files0[i] and "_uim_" in tf_files0[i] and not ("calibcorr" in tf_files0[i]):
response_file = tf_files0[i]
tf_files.append(response_file)
j += 1
indices.append(j)
for i in range(0, len(tf_files)):
if '_filters_transfer_function_' in tf_files[i] and "_uim_" in tf_files[i] and "calibcorr" in tf_files[i]:
response_file = tf_files.pop(i)
k += 1
for i in range(0, len(tf_files0)):
if '_filters_transfer_function_' in tf_files0[i] and "_uim_" in tf_files0[i] and "calibcorr" in tf_files0[i]:
response_file = tf_files0[i]
tf_files.append(response_file)
j += 1
indices.append(j)
for i in range(0, len(tf_files)):
if '_response_filters_transfer_function_' in tf_files[i] and not ("calibcorr" in tf_files[i]):
response_file = tf_files.pop(i)
k += 1
indices.append(k)
for i in range(0, len(tf_files0)):
if '_response_filters_transfer_function_' in tf_files0[i] and not ("calibcorr" in tf_files0[i]):
response_file = tf_files0[i]
tf_files.append(response_file)
j += 1
indices.append(j)
for i in range(0, len(tf_files)):
if '_response_filters_transfer_function_' in tf_files[i] and "calibcorr" in tf_files[i]:
response_file = tf_files.pop(i)
k += 1
for i in range(0, len(tf_files0)):
if '_response_filters_transfer_function_' in tf_files0[i] and "calibcorr" in tf_files0[i]:
response_file = tf_files0[i]
tf_files.append(response_file)
j += 1
indices.append(j)
k += 1
indices.append(k)
# Plot limits
freq_min = 10
......@@ -519,9 +518,9 @@ ratio_mag_max = 1.1
ratio_phase_min = -10
ratio_phase_max = 10
j = 0
k = 0
for tf_file in tf_files:
j += 1
k += 1
model_jump_delay = 0.0
# Determine what transfer function this is
if '_tst_' in tf_file:
......@@ -577,7 +576,7 @@ for tf_file in tf_files:
f.close()
# Read data from re-formatted file and find frequency vector, magnitude, and phase
data = numpy.loadtxt(tf_file.replace('.txt', '_reformatted.txt'))
data = np.loadtxt(tf_file.replace('.txt', '_reformatted.txt'))
os.remove(tf_file.replace('.txt', '_reformatted.txt'))
filters_tf = []
frequency = []
......@@ -586,10 +585,10 @@ for tf_file in tf_files:
df = data[1][0] - data[0][0] # frequency spacing
for j in range(0, len(data)):
frequency.append(data[j][0])
tf_at_f = (data[j][1] + 1j * data[j][2]) * numpy.exp(-2j * numpy.pi * data[j][0] * model_jump_delay)
tf_at_f = (data[j][1] + 1j * data[j][2]) * np.exp(-2j * np.pi * data[j][0] * model_jump_delay)
filters_tf.append(tf_at_f)
magnitude.append(abs(tf_at_f))
phase.append(numpy.angle(tf_at_f) * 180.0 / numpy.pi)
phase.append(np.angle(tf_at_f) * 180.0 / np.pi)
# Find frequency-domain models in filters file if they are present and resample if necessary
model_tf = []
......@@ -604,17 +603,17 @@ for tf_file in tf_files:
index = 0
# This is a linear resampler (it just connects the dots with straight lines)
while index < tf_length:
before_idx = int(numpy.floor(cadence * index))
after_idx = int(numpy.ceil(cadence * index + 1e-10))
before_idx = int(np.floor(cadence * index))
after_idx = int(np.ceil(cadence * index + 1e-10))
# Check if we've hit the end of the model transfer function
if after_idx >= len(model[0]):
if before_idx == cadence * index:
model_tf.append(model[1][before_idx] + 1j * model[2][before_idx])
model_magnitude.append(abs(model_tf[index]))
model_phase.append(numpy.angle(model_tf[index]) * 180.0 / numpy.pi)
model_phase.append(np.angle(model_tf[index]) * 180.0 / np.pi)
ratio.append(filters_tf[index] / model_tf[index])
ratio_magnitude.append(abs(ratio[index]))
ratio_phase.append(numpy.angle(ratio[index]) * 180.0 / numpy.pi)
ratio_phase.append(np.angle(ratio[index]) * 180.0 / np.pi)
index = tf_length
else:
before = model[1][before_idx] + 1j * model[2][before_idx]
......@@ -623,57 +622,56 @@ for tf_file in tf_files:
after_weight = cadence * index - before_idx
model_tf.append(before_weight * before + after_weight * after)
model_magnitude.append(abs(model_tf[index]))
model_phase.append(numpy.angle(model_tf[index]) * 180.0 / numpy.pi)
model_phase.append(np.angle(model_tf[index]) * 180.0 / np.pi)
ratio.append(filters_tf[index] / model_tf[index])
ratio_magnitude.append(abs(ratio[index]))
ratio_phase.append(numpy.angle(ratio[index]) * 180.0 / numpy.pi)
ratio_phase.append(np.angle(ratio[index]) * 180.0 / np.pi)
index += 1
numpy.savetxt(tf_file.replace('.txt', '_ratio_magnitude.txt'), numpy.transpose(numpy.array([frequency, ratio_magnitude])), fmt='%.5e', delimiter='\t')
numpy.savetxt(tf_file.replace('.txt', '_ratio_phase.txt'), numpy.transpose(numpy.array([frequency, ratio_phase])), fmt='%.5f', delimiter='\t')
np.savetxt(tf_file.replace('.txt', '_ratio_magnitude.txt'), np.transpose(np.array([frequency, ratio_magnitude])), fmt='%.5e', delimiter='\t')
np.savetxt(tf_file.replace('.txt', '_ratio_phase.txt'), np.transpose(np.array([frequency, ratio_phase])), fmt='%.5f', delimiter='\t')
# Filter transfer function plots
if j == 1 or j in np.array(indices) + 1:
if k == 1 or k in np.array(indices) + 1:
plt.figure(figsize = (10, 8))
plt.subplot(221)
plt.plot(frequency, magnitude, color, linewidth = 1.0, label = r'${\rm %s \ %s \ %s}$' % (ifo, cal_version, component))
leg = plt.legend(fancybox = True)
leg.get_frame().set_alpha(0.5)
plt.ylabel(r'${\rm Magnitude \ [m/ct]}$')
ticks_and_grid(plt.gca(), xmin = freq_min, xmax = freq_max, ymin = mag_min, ymax = mag_max, xscale = options.tf_frequency_scale, yscale = options.tf_magnitude_scale)
ticks_and_grid(plt.gca(), xmin = freq_min, xmax = freq_max, ymin = mag_min, ymax = mag_max, xscale = 'log', yscale = 'log')
ax = plt.subplot(223)
plt.plot(frequency, phase, color, linewidth = 1.0)
plt.ylabel(r'${\rm Phase [deg]}$')
plt.xlabel(r'${\rm Frequency \ [Hz]}$')
ticks_and_grid(plt.gca(), xmin = freq_min, xmax = freq_max, ymin = phase_min, ymax = phase_max, xscale = options.tf_frequency_scale, yscale = 'linear')
ticks_and_grid(plt.gca(), xmin = freq_min, xmax = freq_max, ymin = phase_min, ymax = phase_max, xscale = 'log', yscale = 'linear')
# Plots of the ratio filters / model
plt.subplot(222)
plt.plot(frequency, ratio_magnitude, color, linewidth = 1.0, label = r'${\rm %s \ %s / Model}$' % (ifo, cal_version))
plt.plot(frequency, ratio_magnitude, color, linewidth = 1.0, label = r'${\rm %s \ %s / Actual}$' % (ifo, cal_version))
leg = plt.legend(fancybox = True)
leg.get_frame().set_alpha(0.5)
ticks_and_grid(plt.gca(), xmin = ratio_freq_min, xmax = ratio_freq_max, ymin = ratio_mag_min, ymax = ratio_mag_max, xscale = options.ratio_frequency_scale, yscale = options.ratio_magnitude_scale)
ticks_and_grid(plt.gca(), xmin = ratio_freq_min, xmax = ratio_freq_max, ymin = ratio_mag_min, ymax = ratio_mag_max, xscale = 'log', yscale = 'linear')
ax = plt.subplot(224)
plt.plot(frequency, ratio_phase, color, linewidth = 1.0)
plt.xlabel(r'${\rm Frequency \ [Hz]}$')
ticks_and_grid(plt.gca(), xmin = ratio_freq_min, xmax = ratio_freq_max, ymin = ratio_phase_min, ymax = ratio_phase_max, xscale = options.ratio_frequency_scale)
ticks_and_grid(plt.gca(), xmin = ratio_freq_min, xmax = ratio_freq_max, ymin = ratio_phase_min, ymax = ratio_phase_max, xscale = 'log', yscale = 'linear')
if j in indices:
# Now add the model
if k in indices:
# Now add the true model
plt.subplot(221)
plt.plot(frequency, model_magnitude, 'orangered', linewidth = 1.0, ls = '--', label = r'${\rm %s \ Model \ %s}$' % (ifo, component))
plt.plot(frequency, model_magnitude, 'orangered', linewidth = 1.0, ls = '--', label = r'${\rm %s \ Actual \ %s}$' % (ifo, component))
leg = plt.legend(fancybox = True)
leg.get_frame().set_alpha(0.5)
#plt.title(plot_title)
plt.ylabel(r'${\rm Magnitude \ [m/ct]}$')
ticks_and_grid(plt.gca(), xmin = freq_min, xmax = freq_max, ymin = mag_min, ymax = mag_max, xscale = options.tf_frequency_scale, yscale = options.tf_magnitude_scale)
ticks_and_grid(plt.gca(), xmin = freq_min, xmax = freq_max, ymin = mag_min, ymax = mag_max, xscale = 'log', yscale = 'log')
ax = plt.subplot(223)
plt.plot(frequency, model_phase, 'orangered', linewidth = 1.0, ls = '--')
plt.ylabel(r'${\rm Phase \ [deg]}$')
plt.xlabel(r'${\rm Frequency \ [Hz]}$')
ticks_and_grid(plt.gca(), xmin = freq_min, xmax = freq_max, ymin = phase_min, ymax = phase_max, xscale = options.tf_frequency_scale)
ticks_and_grid(plt.gca(), xmin = freq_min, xmax = freq_max, ymin = phase_min, ymax = phase_max, xscale = 'log', yscale = 'linear')
plt.savefig(tf_file.replace('.txt', '_ratio.png'))
plt.savefig(tf_file.replace('.txt', '_ratio.pdf'))
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment