From 11d64b174095815fe735b035b83c8285f5e45a2b Mon Sep 17 00:00:00 2001 From: Aaron Viets <aaron.viets@ligo.org> Date: Mon, 25 Jun 2018 13:01:56 -0700 Subject: [PATCH] gstlal-calibration: New pipeline gstlal_compute_kappas which only computes the kappas (much cheaper than computing h(t)). It can also apply a filter to DARM_ERR to simulate the effect of time-dependent correction factors (\kappa_tst, \kappa_pu, \kappa_c, f_cc, f_s, Q), so that tests can be done quickly. --- gstlal-calibration/bin/Makefile.am | 3 +- gstlal-calibration/bin/gstlal_compute_kappas | 1313 +++++++++++++++++ gstlal-calibration/bin/gstlal_compute_strain | 57 +- .../tests/lal_demodulate_test.py | 39 +- 4 files changed, 1387 insertions(+), 25 deletions(-) create mode 100755 gstlal-calibration/bin/gstlal_compute_kappas diff --git a/gstlal-calibration/bin/Makefile.am b/gstlal-calibration/bin/Makefile.am index c5fc06e9c9..be839aa0c4 100644 --- a/gstlal-calibration/bin/Makefile.am +++ b/gstlal-calibration/bin/Makefile.am @@ -1,3 +1,4 @@ dist_bin_SCRIPTS = \ gstlal_compute_strain \ - gstlal_clean_strain + gstlal_clean_strain \ + gstlal_compute_kappas diff --git a/gstlal-calibration/bin/gstlal_compute_kappas b/gstlal-calibration/bin/gstlal_compute_kappas new file mode 100755 index 0000000000..e2df70d866 --- /dev/null +++ b/gstlal-calibration/bin/gstlal_compute_kappas @@ -0,0 +1,1313 @@ +#!/usr/bin/env python +# +# Copyright (C) 2010-2015 Jordi Burguet-Castell, Madeline Wade, 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. + + +""" +This pipeline computes the time-dependent correction factors (kappas) using DARM_ERR and the injection channels PCALY, x_tst, and x_ctrl. It can also read in a filter that is applied to DARM_ERR to simulate the effect of the kappas or undo the effect of the kappas. + +Type gstlal_compute_kappas --help to see the full list of command line options. +""" + +import sys +import numpy +import time +import resource + +from optparse import OptionParser, Option + +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 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 + +# +# Function definition for writing pipeline graph +# + +def write_graph(demux): + pipeparts.write_dump_dot(pipeline, "%s.%s" % (options.write_pipeline, "PLAYING"), verbose = True) + +# +# Make sure we have sufficient resources +# We allocate far more memory than we need, so this is okay +# + +def setrlimit(res, lim): + hard_lim = resource.getrlimit(res)[1] + resource.setrlimit(res, (lim if lim is not None else hard_lim, hard_lim)) +# set the number of processes and total set size up to hard limit and +# shrink the per-thread stack size (default is 10 MiB) +setrlimit(resource.RLIMIT_NPROC, None) +setrlimit(resource.RLIMIT_AS, None) +setrlimit(resource.RLIMIT_RSS, None) +setrlimit(resource.RLIMIT_STACK, 1024*1024) + +# +# Function definition to obtain the current GPS time +# + +def now(): + return lal.LIGOTimeGPS(lal.UTCToGPS(time.gmtime()), 0) + + +################################################################################################### +############################## Program Command Line Options ####################################### +################################################################################################### + +parser = OptionParser(description = __doc__) + +# Append program specific options + +# These options should be used whether the pipeline runs in full calibration mode or partial calibration mode +parser.add_option("--data-source", metavar = "source", help = "Set the data source from [frames|lvshm]. Required.") +parser.add_option("--frame-cache", metavar = "filename", help = "Set the name of the LAL cache listing the LIGO .gwf frame files (optional). This is required iff --data-source=frames") +parser.add_option("--gps-start-time", metavar = "seconds", help = "Set the start time of the segment to analyze in GPS seconds. This is required iff --data-source=frames") +parser.add_option("--gps-end-time", metavar = "seconds", help = "Set the end time of the segment to analyze in GPS seconds. This is required iff --data-source=frames") +parser.add_option("--wings", metavar = "seconds", type = "int", help = "Number of seconds to trim off of the beginning and end of the output. Should only be used if --data-source=frames.") +parser.add_option("--do-file-checksum", action = "store_true", help = "Set this option to turn on file checksum in the demuxer.") +parser.add_option("--ifo", metavar = "name", help = "Name of the IFO to be calibrated.") +parser.add_option("--shared-memory-partition", metavar = "name", help = "Set the name of the shared memory partition to read from. This is required iff --data-source=lvshm.") +parser.add_option("--derr-sample-rate", metavar = "Hz", default = 16384, type = "int", help = "Sample rate of the error signal channel. (Default = 16384 Hz)") +parser.add_option("--kappas-state-sample-rate", metavar = "Hz", default = 16, type = "int", help = "Sample rate for the outgoing DQ vector GDS-KAPPAS_STATE_VECTOR. (Default = 16 Hz)") +parser.add_option("--tst-exc-sample-rate", metavar = "Hz", default = 512, type = "int", help = "Sample rate for the control signals being read in. (Default = 512 Hz)") +parser.add_option("--coh-sample-rate", metavar = "Hz", default = 16, type = "int", help = "Sample rate for the coherence uncertainty channels. (Default = 16 Hz).") +parser.add_option("--buffer-length", metavar = "seconds", type = float, default = 1.0, help = "Set the length in seconds of buffers to be used in the pipeline (Default = 1.0)") +parser.add_option("--frame-duration", metavar = "seconds", type = "int", default = 4, help = "Set the number of seconds for each frame. (Default = 4)") +parser.add_option("--frames-per-file", metavar = "count", type = "int", default = 1, help = "Set the number of frames per frame file. (Default = 1)") +parser.add_option("--frame-size", metavar = "bytes", type = "int", default = 405338, help = "Approximate size in bytes of frame file images; used when writing to shared memory. (Default=405338)") +parser.add_option("--compression-scheme", metavar = "scheme", type = "int", default = 256, help = "Set the compression scheme for the framecpp_channelmux element. (Default=256, no compression)") +parser.add_option("--compression-level", metavar = "level", type = "int", default = 0, help = "Set the compression level for the framecpp_channelmux element. (Default=0)") +parser.add_option("--write-to-shm-partition", metavar = "name", help = "Set the name of the shared memory partition to write to. If this is not provided, frames will be written to a file.") +parser.add_option("--buffer-mode", metavar = "number", type = "int", default = 2, help = "Set the buffer mode for the lvshmsink element. (Default=2)") +parser.add_option("--frame-type", metavar = "name", default = "TEST", help = "Set the frame type as input to the frame writing element. (Default=TEST)") +parser.add_option("--output-path", metavar = "name", default = ".", help = "Set the output path for writing frame files. (Default=Current)") +parser.add_option("--no-dq-vector", action = "store_true", help = "Set this if you want to turn off all interpretation and calculation of a data quality vector.") +parser.add_option("--chan-prefix", metavar = "name", default = "GDS-", help = "Prefix for all output channel names. (Default = GDS)") +parser.add_option("--chan-suffix", metavar = "name", help = "Suffix for all output channel names.") + +# These are debugging options +parser.add_option("--write-pipeline", metavar = "filename", help = "Write a DOT graph description of the as-built pipeline to this file (optional). The environment variable GST_DEBUG_DUMP_DOT_DIR must be set for this option to work.") +parser.add_option("-v", "--verbose", action = "store_true", help = "Be verbose (optional).") + +# These are options specific to the calibration procedure +parser.add_option("--filters-file", metavar="filename", help = "Name of file containing calibration filters (in npz format)") +parser.add_option("--kappas-filters-file", metavar="filename", help = "Name of file containing kappas simulation filters (in npz format)") +parser.add_option("--factors-from-filters-file", action = "store_true", help = "Compute the calibration factors from reference values contained in the filters file instead of from EPICS channels.") +parser.add_option("--no-coherence", action = "store_true", help = "Gate the calibration factors with a pre-computed coherence channel.") +parser.add_option("--coherence-uncertainty-threshold", metavar = "float", type = float, default = 0.0025, help = "Threshold for the coherence uncertainty for each calibration line. (Default = 0.0025)") +parser.add_option("--coherence-time", metavar = "seconds", type = "int", default = 130, help = "Amount of time used in front end to compute coherence of calibration lines. (Default = 130)") +parser.add_option("--coh-unc-sus-line1-channel", metavar="name", default="CAL-CS_TDEP_SUS_LINE1_UNCERTAINTY", help = "Channel name for SUS line 1 coherence uncertainty. (Default=CAL-CS_TDEP_SUS_LINE1_UNCERTAINTY)") +parser.add_option("--coh-unc-pcaly-line1-channel", metavar="name", default="CAL-CS_TDEP_PCALY_LINE1_UNCERTAINTY", help = "Channel name for PCALY line 1 coherence uncertainty. (Default=CAL-CS_TDEP_PCALY_LINE1_UNCERTAINTY)") +parser.add_option("--coh-unc-pcaly-line2-channel", metavar="name", default="CAL-CS_TDEP_PCALY_LINE2_UNCERTAINTY", help = "Channel name for PCALY line 2 coherence uncertainty. (Default=CAL-CS_TDEP_PCALY_LINE2_UNCERTAINTY)") +parser.add_option("--coh-unc-darm-line1-channel", metavar="name", default="CAL-CS_TDEP_DARM_LINE1_UNCERTAINTY", help = "Channel name for DARM line 1 coherence uncertainty. (Default=CAL-CS_TDEP_DARM_LINE1_UNCERTAINTY)") +parser.add_option("--no-kappatst", action = "store_true", help = "Set this to turn off the calculation of \kappa_tst.") +parser.add_option("--no-kappapu", action = "store_true", help = "Set this to turn off the calculation of \kappa_pu.") +parser.add_option("--no-kappap", action = "store_true", help = "Set this to turn off the calculation of \kappa_p.") +parser.add_option("--no-kappau", action = "store_true", help = "Set this to turn off the calculation of \kappa_u.") +parser.add_option("--no-kappac", action = "store_true", help = "Set this to turn off the calculation of \kappa_c.") +parser.add_option("--no-fcc", action = "store_true", help = "Set this to turn off the calculation of f_cc.") +parser.add_option("--no-srcQ", action = "store_true", help = "Set this to turn off the calculation of the SRC Q.") +parser.add_option("--no-fs", action = "store_true", help = "Set this to turn off the calculation of the SRC spring frequency.") +parser.add_option("--act-timing-from-kappapu", action = "store_true", help = "Set this to use the calculated value of \kappa_pu to measure any timing error in the actuation. If this is set, the phase of \kappa_tst will be adjusted accordingly.") +parser.add_option("--act-timing-from-kappatst", action = "store_true", help = "Set this to use the calculated value of \kappa_tst to measure any timing error in the actuation. If this is set, the phase of \kappa_pu will be adjusted accordingly.") +parser.add_option("--factors-averaging-time", metavar = "Sec", type = int, default = 10, help = "Time over which to average the smoothed time-varying calibration factors (\kappas), given in seconds. (Default = 10 seconds)") +parser.add_option("--src-averaging-time", metavar = "Sec", type = int, default = 600, help = "Time over which to average the smoothed SRC detuning parameters fs and Q, given in seconds. (Default = 600 seconds)") +parser.add_option("--compute-factors-sr", metavar = "Hz", type = int, default = 16, help = "Sample rate at which time-dependent correction factors are computed. (Default = 16 Hz)") +parser.add_option("--demodulation-filter-time", metavar = "s", type = int, default = 20, help = "Length in seconds of low-pass FIR filter used in demodulation of the calibration lines. (Default = 20 seconds)") +parser.add_option("--median-smoothing-time", metavar = "s", type = int, default = 128, help = "Time (in seconds) to smooth out \kappas with a median-like method. (Default = 128 s)") +parser.add_option("--kappas-default-to-median", action = "store_true", help = "If set, bad computed kappas will be replaced by the previous computed median in the running median array. Otherwise, they are replaced with the default value.") +parser.add_option("--filter-latency", metavar = "float", type = float, default = 0.0, help = "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. (Default = 0.0)") +parser.add_option("--record-factors-sr", metavar = "Hz", type = int, default = 16, help = "Sample rate at which calibration factors are recorded. (Default = 16 Hz)") +parser.add_option("--expected-kappapu-real", metavar = "float", type = float, default = 1.0, help = "Expected value for the real part of \kappa_pu. (Default = 1.0)") +parser.add_option("--expected-kappap-real", metavar = "float", type = float, default = 1.0, help = "Expected value for the real part of \kappa_p. (Default = 1.0)") +parser.add_option("--expected-kappau-real", metavar = "float", type = float, default = 1.0, help = "Expected value for the real part of \kappa_u. (Default = 1.0)") +parser.add_option("--expected-kappatst-real", metavar = "float", type = float, default = 1.0, help = "Expected value for the real part of \kappa_tst. (Default = 1.0)") +parser.add_option("--expected-kappapu-imag", metavar = "float", type = float, default = 0.0, help = "Expected value for the imaginary part of \kappa_pu. (Default = 0.0)") +parser.add_option("--expected-kappap-imag", metavar = "float", type = float, default = 0.0, help = "Expected value for the imaginary part of \kappa_p. (Default = 0.0)") +parser.add_option("--expected-kappau-imag", metavar = "float", type = float, default = 0.0, help = "Expected value for the imaginary part of \kappa_u. (Default = 0.0)") +parser.add_option("--expected-kappatst-imag", metavar = "float", type = float, default = 0.0, help = "Expected value for the imaginary part of \kappa_tst. (Default = 0.0)") +parser.add_option("--expected-kappac", metavar = "float", type = float, default = 1.0, help = "Expected value for \kappa_c. (Default = 1.0)") +parser.add_option("--expected-fcc", metavar = "Hz", type = float, default = 360.0, help = "Expected value for the coupled cavity pole. (Default = 360.0 Hz)") +parser.add_option("--expected-fs", metavar = "Hz", type = float, default = 8.0, help = "Expected value for the SRC optical spring frequency. (Default = 8.0 Hz)") +parser.add_option("--expected-srcQ", metavar = "float", type = float, default = 28.0, help = "Expected value for the SRC Q. (Default = 28.0)") +parser.add_option("--kappapu-real-ok-var", metavar = "float", type = float, default = 0.2, help = "Values of the real part of \kappa_pu +/- this number will be considered OK. (Default = 0.2)") +parser.add_option("--kappap-real-ok-var", metavar = "float", type = float, default = 0.2, help = "Values of the real part of \kappa_p +/- this number will be considered OK. (Default = 0.2)") +parser.add_option("--kappau-real-ok-var", metavar = "float", type = float, default = 0.2, help = "Values of the real part of \kappa_u +/- this number will be considered OK. (Default = 0.2)") +parser.add_option("--kappatst-real-ok-var", metavar = "float", type = float, default = 0.2, help = "Values of the real part of \kappa_tst +/- this number will be considered OK. (Default = 0.2)") +parser.add_option("--kappapu-imag-ok-var", metavar = "float", type = float, default = 0.2, help = "Values of the imaginary part of \kappa_pu +/- this number will be considered OK. (Default = 0.2)") +parser.add_option("--kappap-imag-ok-var", metavar = "float", type = float, default = 0.2, help = "Values of the imaginary part of \kappa_p +/- this number will be considered OK. (Default = 0.2)") +parser.add_option("--kappau-imag-ok-var", metavar = "float", type = float, default = 0.2, help = "Values of the imaginary part of \kappa_u +/- this number will be considered OK. (Default = 0.2)") +parser.add_option("--kappatst-imag-ok-var", metavar = "float", type = float, default = 0.2, help = "Values of the imaginary part of \kappa_tst +/- this number will be considered OK. (Default = 0.2)") +parser.add_option("--kappac-ok-var", metavar = "float", type = float, default = 0.2, help = "Values of \kappa_c +/- this number will be considered OK. (Default = 0.2)") +parser.add_option("--fcc-ok-var", metavar = "Hz", type = float, default = 50, help = "Values of f_cc +/- this number (in Hz) will be considered OK. (Default = 50 Hz)") +parser.add_option("--fs-ok-var", metavar = "Hz", type = float, default = 5, help = "Values of SRC spring frequency +/- this number (in Hz) will be considered OK. (Default = 5 Hz)") +parser.add_option("--srcQinv-min", metavar = "float", type = float, default = 0.0, help = "Minimum value of SRC Q inverse that will be considered OK. (Default = 0.0)") +parser.add_option("--srcQinv-max", metavar = "float", type = float, default = 0.5, help = "Maximum value of SRC Q inverse that will be considered OK. (Default = 0.5)") +parser.add_option("--exc-channel-name", metavar = "name", default = "CAL-CS_LINE_SUM_DQ", help = "Set the name of the excitation channel. This is only necessary when the calibration factors computation is turned on, which is the default behavior. (Default = CAL-CS_LINE_SUM_DQ)") +parser.add_option("--tst-exc-channel-name", metavar = "name", default = "SUS-ETMY_L3_CAL_LINE_OUT_DQ", help = "Set the name of the TST excitation channel. This is only necessary when the \kappa_tst factors computation is turned on, which is the default behavior. (Default = SUS-ETMY_L3_CAL_LINE_OUT_DQ)") +parser.add_option("--pcal-channel-name", metavar = "name", default = "CAL-PCALY_RX_PD_OUT_DQ", help = "Set the name of the PCal channel used for calculating the calibration factors. (Default = CAL-PCALY_RX_PD_OUT_DQ)") +parser.add_option("--darm-err-channel-name", metavar = "name", default = "CAL-DARM_ERR_WHITEN_OUT_DBL_DQ", help = "Set the name of the error signal channel. (Default = CAL-DARM_ERR_WHITEN_OUT_DBL_DQ)") + +# These are all options related to the reference channels used in the calibration factors computation +parser.add_option("--ref-channels-sr", metavar = "Hz", default = 16, help = "Set the sample rate for the reference model channels used in the calibration factors calculation. (Default = 16 Hz)") +parser.add_option("--EP4-real", metavar = "name", default = "CAL-CS_TDEP_DARM_LINE1_REF_A_TST_REAL", help = "Set the name of the channel containing the real part of A_tst at the ESD line used for the \kappa_a and \kappa_pu calculation. (Default = CAL-CS_TDEP_DARM_LINE1_REF_A_TST_REAL)") +parser.add_option("--EP5-real", metavar = "name", default = "CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_REAL", help = "Set the name of the channel containing the real part of A_pu at the ESD line used for the \kappa_a calculation. (Default = CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_REAL)") +parser.add_option("--EP3-real", metavar = "name", default = "CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_INV_REAL", help = "Set the name of the channel containing the real part of 1/A_pu at the ESD line used for the \kappa_pu calculation. (Default = CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_INV_REAL)") +parser.add_option("--EP4-imag", metavar = "name", default = "CAL-CS_TDEP_DARM_LINE1_REF_A_TST_IMAG", help = "Set the name of the channel containing the imaginary part of A_tst at the ESD line used for the \kappa_a and \kappa_pu calculation. (Default = CAL-CS_TDEP_DARM_LINE1_REF_A_TST_IMAG") +parser.add_option("--EP5-imag", metavar = "name", default = "CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_IMAG", help = "Set the name of the channel containing the imaginary part of A_pu at the ESD line used for the \kappa_A calculation. (Default = CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_IMAG") +parser.add_option("--EP3-imag", metavar = "name", default = "CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_INV_IMAG", help = "Set the name of the channel containing the imaginary part of 1/A_pu at the ESD line used for the \kappa_PU calculation. (Default = CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_INV_IMAG") +parser.add_option("--EP2-real", metavar = "name", default = "CAL-CS_TDEP_REF_CLGRATIO_CTRL_REAL", help = "Set the name of the channel containing the real part of the factors used to compute A(f_ctrl). (Default = CAL-CS_TDEP_REF_CLGRATIO_CTRL_REAL)") +parser.add_option("--EP2-imag", metavar = "name", default = "CAL-CS_TDEP_REF_CLGRATIO_CTRL_IMAG", help = "Set the name of the channel containing the imaginary part of the factors used to compute A(f_ctrl). (Default = CAL-CS_TDEP_REF_CLGRATIO_CTRL_IMAG)") +parser.add_option("--EP6-real", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE2_REF_C_NOCAVPOLE_REAL", help = "Set the name of the channel containing the real part of C_res at the PCal line used for the \kappa_c and f_cc calculation. (Default = CAL-CS_TDEP_PCALY_LINE2_REF_C_NOCAVPOLE_REAL") +parser.add_option("--EP6-imag", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE2_REF_C_NOCAVPOLE_IMAG", help = "Set the name of the channel containing the imaginary part of C_res at the PCal line used for the \kappa_c and f_cc calculation. (Default = CAL-CS_TDEP_PCALY_LINE2_REF_C_NOCAVPOLE_IMAG") +parser.add_option("--EP7-real", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE2_REF_D_REAL", help = "Set the name of the channel containing the real part of D at the PCal line used for the \kappa_c and f_cc calculation. (Default = CAL-CS_TDEP_PCALY_LINE2_REF_D_REAL") +parser.add_option("--EP7-imag", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE2_REF_D_IMAG", help = "Set the name of the channel containing the imaginary part of D at the PCal line used for the \kappa_c and f_cc calculation. (Default = CAL-CS_TDEP_PCALY_LINE2_REF_D_IMAG") +parser.add_option("--EP8-real", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE2_REF_A_TST_REAL", help = "Set the name of the channel containing the real part of A_tst at the PCal line used for the \kappa_c and f_cc calculation. (Default = CAL-CS_TDEP_PCALY_LINE2_REF_A_TST_REAL") +parser.add_option("--EP8-imag", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE2_REF_A_TST_IMAG", help = "Set the name of the channel containing the imaginary part of A_tst at the PCal line used for the \kappa_c and f_cc calculation. (Default = CAL-CS_TDEP_PCALY_LINE2_REF_A_TST_IMAG") +parser.add_option("--EP9-real", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE2_REF_A_USUM_REAL", help = "Set the name of the channel containing the real part of A_pu at the PCal line used for the \kappa_c and f_cc calculation. (Default = CAL-CS_TDEP_PCALY_LINE2_REF_A_USUM_REAL") +parser.add_option("--EP9-imag", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE2_REF_A_USUM_IMAG", help = "Set the name of the channel containing the imaginary part of A_pu at the PCal line used for the \kappa_c and f_cc calculation. (Default = CAL-CS_TDEP_PCALY_LINE2_REF_A_USUM_IMAG") +parser.add_option("--EP1-real", metavar = "name", default = "CAL-CS_TDEP_REF_INVA_CLGRATIO_TST_REAL", help = "Set the name of the channel containing the real part of the \kappa_tst reference factors. (Default = CAL-CS_TDEP_REF_INVA_CLGRATIO_TST_REAL)") +parser.add_option("--EP1-imag", metavar = "name", default = "CAL-CS_TDEP_REF_INVA_CLGRATIO_TST_IMAG", help = "Set the name of the channel containing the imaginary part of the \kappa_tst reference factors. (Default = CAL-CS_TDEP_REF_INVA_CLGRATIO_TST_IMAG)") +parser.add_option("--EP10-real", metavar = "name", default = "CAL-CS_TDEP_ESD_LINE1_REF_A_TST_NOLOCK_REAL", help = "Set the name of the channel containing the real part of A_tst at the ESD line used for removal of the ESD line. (Default = CAL-CS_TDEP_ESD_LINE1_REF_A_TST_REAL") +parser.add_option("--EP10-imag", metavar = "name", default = "CAL-CS_TDEP_ESD_LINE1_REF_A_TST_NOLOCK_IMAG", help = "Set the name of the channel containing the imaginary part of A_tst at the ESD line used for removal of the ESD line. (Default = CAL-CS_TDEP_ESD_LINE1_REF_A_TST_IMAG") +parser.add_option("--EP11-real", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE4_REF_C_NOCAVPOLE_REAL", help = "Set the name of the channel containing the real part of C_res at the PCal line used for the src_Q and f_s calculation. (Default = CAL-CS_TDEP_PCALY_LINE4_REF_C_NOCAVPOLE_REAL") +parser.add_option("--EP11-imag", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE4_REF_C_NOCAVPOLE_IMAG", help = "Set the name of the channel containing the imaginary part of C_res at the PCal line used for the src_Q and f_s calculation. (Default = CAL-CS_TDEP_PCALY_LINE4_REF_C_NOCAVPOLE_IMAG") +parser.add_option("--EP12-real", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE4_REF_D_REAL", help = "Set the name of the channel containing the real part of D at the PCal line used for the src_Q and f_s calculation. (Default = CAL-CS_TDEP_PCALY_LINE4_REF_D_REAL") +parser.add_option("--EP12-imag", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE4_REF_D_IMAG", help = "Set the name of the channel containing the imaginary part of D at the PCal line used for the src_Q and f_s calculation. (Default = CAL-CS_TDEP_PCALY_LINE4_REF_D_IMAG") +parser.add_option("--EP13-real", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE4_REF_A_TST_REAL", help = "Set the name of the channel containing the real part of A_tst at the PCal line used for the src_Q and f_s calculation. (Default = CAL-CS_TDEP_PCALY_LINE4_REF_A_TST_REAL") +parser.add_option("--EP13-imag", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE4_REF_A_TST_IMAG", help = "Set the name of the channel containing the imaginary part of A_tst at the PCal line used for the src_Q and f_s calculation. (Default = CAL-CS_TDEP_PCALY_LINE4_REF_A_TST_IMAG") +parser.add_option("--EP14-real", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE4_REF_A_USUM_REAL", help = "Set the name of the channel containing the real part of A_pu at the PCal line used for the src_Q and f_s calculation. (Default = CAL-CS_TDEP_PCALY_LINE4_REF_A_USUM_REAL") +parser.add_option("--EP14-imag", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE4_REF_A_USUM_IMAG", help = "Set the name of the channel containing the imaginary part of A_pu at the PCal line used for the src_Q and f_s calculation. (Default = CAL-CS_TDEP_PCALY_LINE4_REF_A_USUM_IMAG") + + +# Parse options +options, filenames = parser.parse_args() + +# Sanity checks for command line options +data_sources = set(("frames", "lvshm")) + +if options.data_source not in data_sources: + raise ValueError("--data-source must be one of %s" % ",".join(data_sources)) + +if options.data_source == "frames" and options.frame_cache is None: + raise ValueError("--frame-cache must be specified when using --data-source=frames") + +if options.no_kappatst and options.no_kappapu and options.no_kappap and options.no_kappau and options.no_kappac and options.no_fcc and options.no_fs and options.no_srcQ: + raise ValueError("Must compute at least one of the time-dependent correction factors") + +if options.wings is not None and options.data_source != "frames": + raise ValueError("--wings can only be set when --data-source=frames") + +if options.ifo is None: + raise ValueError("must specify --ifo") + +if options.data_source == "frames" and (options.gps_start_time is None or options.gps_end_time is None): + raise ValueError("must specify --gps-start-time and --gps-end-time when --data-source=frames") + +if int(options.record_factors_sr) > int(options.compute_factors_sr): + raise ValueError("--record-factors-sr must be less than or equal to --compute-factors-sr") + +if not options.factors_from_filters_file and (not options.no_fs or not options.no_srcQ) and ((options.data_source == "frames" and int(options.gps_start_time) < 1175954418) or (options.data_source == "lvshm" and now() < 1175954418)): + raise ValueError("Cannot compute SRC detuning parameters as the needed EPICS channels are not in the frames until GPS time 1175954418. Use command line options --no-srcQ and --no-fs.") + +if options.gps_start_time is not None: + if options.gps_end_time is None: + raise ValueError("must provide both --gps-start-time and --gps-end-time") + if options.data_source == "lvshm": + raise ValueError("cannot set --gps-start-time or --gps-end-time with --data-source=lvshm") + try: + start = lal.LIGOTimeGPS(options.gps_start_time) + except ValueError: + raise ValueError("invalid --gps-start-time %s" % options.gps_start_time) + try: + end = lal.LIGOTimeGPS(options.gps_end_time) + except ValueError: + raise ValueError("invalid --gps-end-time %s" % options.gps_end_time) + if start >= end: + raise ValueError("--gps-start-time must be < --gps-end-time: %s < %s" % (options.gps_start_time, options.gps_end_time)) + # segment from gps start and stop time if given + seg = segments.segment(start, end) + gps_start_time = seg[0] + gps_end_time = seg[1] +elif options.gps_end_time is not None: + raise ValueError("must provide both --gps-start-time and --gps-end-time") + +################################################################################################### +######################################## Setup #################################################### +################################################################################################### + +# Set up instrument and channel name info from command line options +instrument = options.ifo + +# Set up short-cut names for each of the sample rates used throughout the pipeline and establish caps string shortcuts +derr_sr = options.derr_sample_rate # Sample rate for the error signal +kappasstatesr = options.kappas_state_sample_rate # Sample rate for the KAPPA_STATE_VECTOR +cohsr = options.coh_sample_rate # Sample rate for the coherence uncertainty channels +derr_caps = "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % derr_sr +calibstate_caps = "audio/x-raw, format=U32LE, rate=%d, channel-mask=(bitmask)0x0" % kappasstatesr +coh_caps = "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % cohsr +# caps strings for the computation kappas +ref_factors_caps = "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % options.ref_channels_sr +compute_calib_factors_caps = "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0X0" % options.compute_factors_sr +compute_calib_factors_complex_caps = "audio/x-raw, format=Z128LE, rate=%d, channel-mask=(bitmask)0x0" % options.compute_factors_sr + +# Set up smoothing, averaging and integration sample sizes for kappa calulations +integration_samples = int(options.demodulation_filter_time) * options.compute_factors_sr +factors_average_samples = int(options.factors_averaging_time) * options.compute_factors_sr +src_average_samples = int(options.src_averaging_time) * options.compute_factors_sr +median_smoothing_samples = int(options.median_smoothing_time) * options.compute_factors_sr +kappa_gate_attack_length = -integration_samples * (1.0 - options.filter_latency) +kappa_gate_hold_length = -integration_samples * options.filter_latency - (options.filter_latency != 0) * options.coherence_time * options.compute_factors_sr + +# Set up string for the channels suffix and prefix as provided by the user +if options.chan_suffix is not None: + chan_suffix = options.chan_suffix +else: + chan_suffix = "" +chan_prefix = options.chan_prefix + +# How many EPICS will we for the KAPPA_STATE_VECTOR calculation? It depends on the IFO and the time we are calibrating +if options.no_dq_vector: + num_dq_epics = 0 +elif options.ifo == "H1" and int(options.gps_start_time) > 1175976256: + num_dq_epics = 14 +elif options.ifo == "L1" and int(options.gps_start_time) > 1179588864: + num_dq_epics = 10 +else: + num_dq_epics = 9 + +# +# Load in the filters files that contains filter coefficients, etc. +# + +filters = numpy.load(options.filters_file) +if options.kappas_filters_file is not None: + kappas_filters = numpy.load(options.kappas_filters_file) + derr_kappas_filter = kappas_filters["kappas_filter_td"] + kappas_filter_sr = kappas_filters["kappas_filter_sr"] + kappas_filter_delay = kappas_filters["kappas_filter_delay"] +else: + derr_kappas_filter = numpy.ones(1) + kappas_filter_delay = 0 + kappas_filter_sr = derr_sr + +# If we're reading the reference model factors from the filters file, load them +if options.factors_from_filters_file or not options.no_dq_vector: + try: + EP1_real = float(filters["EP1_real"]) + EP1_imag = float(filters["EP1_imag"]) + EP2_real = float(filters["EP2_real"]) + EP2_imag = float(filters["EP2_imag"]) + EP3_real = float(filters["EP3_real"]) + EP3_imag = float(filters["EP3_imag"]) + EP4_real = float(filters["EP4_real"]) + EP4_imag = float(filters["EP4_imag"]) + EP5_real = float(filters["EP5_real"]) + EP5_imag = float(filters["EP5_imag"]) + EP6_real = float(filters["EP6_real"]) + EP6_imag = float(filters["EP6_imag"]) + EP7_real = float(filters["EP7_real"]) + EP7_imag = float(filters["EP7_imag"]) + EP8_real = float(filters["EP8_real"]) + EP8_imag = float(filters["EP8_imag"]) + EP9_real = float(filters["EP9_real"]) + EP9_imag = float(filters["EP9_imag"]) + except: + if options.factors_from_filters_file: + raise ValueError("Cannot compute time-dependent correction factors, as the needed EPICS are not contained in the specified filters file.") + try: + EP10_real = float(filters["EP10_real"]) + EP10_imag = float(filters["EP10_imag"]) + except: + if num_dq_epics > 9: + num_dq_epics = 9 + try: + EP11_real = float(filters["EP11_real"]) + EP11_imag = float(filters["EP11_imag"]) + EP12_real = float(filters["EP12_real"]) + EP12_imag = float(filters["EP12_imag"]) + EP13_real = float(filters["EP13_real"]) + EP13_imag = float(filters["EP13_imag"]) + EP14_real = float(filters["EP14_real"]) + EP14_imag = float(filters["EP14_imag"]) + except: + if options.factors_from_filters_file and (not options.no_srcQ or not options.no_fs): + raise ValueError("Cannot compute SRC spring frequency or Q, as the needed EPICS are not contained in the specified filters file.") + if num_dq_epics > 10: + num_dq_epics = 10 + +# Load all of the kappa dewhitening and correction factors +darm_act_line_freq = float(filters["ka_pcal_line_freq"]) +pcal_corr_at_darm_act_freq_real = float(filters["ka_pcal_corr_re"]) +pcal_corr_at_darm_act_freq_imag = float(filters["ka_pcal_corr_im"]) +pu_act_esd_line_freq = float(filters["ka_esd_line_freq"]) +opt_gain_fcc_line_freq = float(filters["kc_pcal_line_freq"]) +pcal_corr_at_opt_gain_fcc_freq_real = float(filters["kc_pcal_corr_re"]) +pcal_corr_at_opt_gain_fcc_freq_imag = float(filters["kc_pcal_corr_im"]) +esd_act_line_freq = float(filters["ktst_esd_line_freq"]) +try: + src_pcal_line_freq = float(filters["src_pcal_line_freq"]) + pcal_corr_at_src_freq_real = float(filters["src_pcal_corr_re"]) + pcal_corr_at_src_freq_imag = float(filters["src_pcal_corr_im"]) +except: + if not options.no_srcQ or not options.no_fs: + raise ValueError("Cannot compute SRC spring frequency or Q, as the calibration line frequency is not contained in the specified filters file.") +try: + fcc_default = float(filters["fcc"]) +except: + fcc_default = options.expected_fcc +try: + fs_default = float(filters["fs"]) + srcQ_default = float(filters["srcQ"]) +except: + fs_default = options.expected_fs + srcQ_default = options.expected_srcQ + +# +# Set up the appropriate channel list. In this section, we also fill a list called headkeys +# that will be the keys for the dictionary holding each pipeline branch name, and we set up +# a dictionary that will be populated with pipeline branch names based on the channel list. +# + +head_dict = {} +channel_list = [] +headkeys = [] + +# We need DARM_ERR and at least the pcal and tst excitation channels for computing kappas + +channel_list.extend(((instrument, options.darm_err_channel_name), (instrument, options.pcal_channel_name), (instrument, options.tst_exc_channel_name))) +headkeys.extend(("darm_err", "pcal", "tstexc")) + +# Most likely, we need the excitation channel as well +if not options.no_kappac or not options.no_fcc or not options.no_kappapu or not options.no_fs or not options.no_srcQ: + channel_list.append((instrument, options.exc_channel_name)) + headkeys.append("exc") + +# If we are not reading the EPICS records from the filters, we need them from the frames +# Needed for kappa_tst +if not options.factors_from_filters_file and (not options.no_kappatst or not options.no_kappapu or not options.no_kappap or not options.no_kappau or not options.no_kappac or not options.no_fcc or not options.no_fs or not options.no_srcQ) or num_dq_epics > 0: + channel_list.extend(((instrument, options.EP1_real), (instrument, options.EP1_imag))) + headkeys.extend(("EP1_real", "EP1_imag")) +# These are needed for kappa_pu +if not options.factors_from_filters_file and (not options.no_kappac or not options.no_fcc or not options.no_kappapu or not options.no_fs or not options.no_srcQ) or num_dq_epics > 3: + channel_list.extend(((instrument, options.EP2_real), (instrument, options.EP2_imag), (instrument, options.EP3_real), (instrument, options.EP3_imag), (instrument, options.EP4_real), (instrument, options.EP4_imag))) + headkeys.extend(("EP2_real", "EP2_imag", "EP3_real", "EP3_imag", "EP4_real", "EP4_imag")) +# If we are computing either kappa_c or f_cc, we need some more EPICS records +if not options.factors_from_filters_file and (not options.no_kappac or not options.no_fcc or not options.no_fs or not options.no_srcQ) or num_dq_epics > 8: + channel_list.extend(((instrument, options.EP6_real), (instrument, options.EP6_imag), (instrument, options.EP7_real), (instrument, options.EP7_imag), (instrument, options.EP8_real), (instrument, options.EP8_imag), (instrument, options.EP9_real), (instrument, options.EP9_imag))) + headkeys.extend(("EP6_real", "EP6_imag", "EP7_real", "EP7_imag", "EP8_real", "EP8_imag", "EP9_real", "EP9_imag")) + +# EP10 is needed to remove the ESD line, so it only gets used to compare to the filters-file EP10 here +if num_dq_epics > 9: + channel_list.extend(((instrument, options.EP10_real), (instrument, options.EP10_imag))) + headkeys.extend(("EP10_real", "EP10_imag")) + +# These are needed if we compute the optical spring frequency and/or Q-factor of the Signal Recycling Cavity (SRC) +if not options.factors_from_filters_file and (not options.no_fs or not options.no_srcQ) or num_dq_epics > 13: + channel_list.extend(((instrument, options.EP11_real), (instrument, options.EP11_imag), (instrument, options.EP12_real), (instrument, options.EP12_imag), (instrument, options.EP13_real), (instrument, options.EP13_imag), (instrument, options.EP14_real), (instrument, options.EP14_imag))) + headkeys.extend(("EP11_real", "EP11_imag", "EP12_real", "EP12_imag", "EP13_real", "EP13_imag", "EP14_real", "EP14_imag")) + +# If we are using pre-computed coherence to gate kappas +if not options.no_coherence: + if not options.no_kappatst or not options.no_kappapu or not options.no_kappac or not options.no_fcc or not options.no_fs or not options.no_srcQ: + channel_list.extend(((instrument, options.coh_unc_sus_line1_channel), (instrument, options.coh_unc_pcaly_line1_channel), (instrument, options.coh_unc_darm_line1_channel))) + headkeys.extend(("pcaly_line1_coh", "sus_coh", "darm_coh")) + if not options.no_kappac or not options.no_fcc or not options.no_fs or not options.no_srcQ: + channel_list.append((instrument, options.coh_unc_pcaly_line2_channel)) + headkeys.append("pcaly_line2_coh") + + +#################################################################################################### +####################################### Main Pipeline ############################################## +#################################################################################################### + +pipeline = Gst.Pipeline(name="gstlal_compute_kappas") +mainloop = GObject.MainLoop() +handler = simplehandler.Handler(mainloop, pipeline) + +# +# Turn off debugging tools or verboseness +# + +pipeparts.mkchecktimestamps = lambda pipeline, src, *args: src # comment this line out to turn on the checktimestamps debugging +if not options.verbose: + pipeparts.mkprogressreport = lambda pipeline, src, *args: src + +# +# Read in data from frames or shared memory +# + +if options.data_source == "lvshm": # Data is to be read from shared memory; "low-latency" mode + src = pipeparts.mklvshmsrc(pipeline, shm_name = options.shared_memory_partition, assumed_duration = 1) +elif options.data_source == "frames": # Data is to be read from frame files; "offline" mode + src = pipeparts.mklalcachesrc(pipeline, location = options.frame_cache, cache_dsc_regex = instrument) + +# +# Hook up the relevant channels to the demuxer +# + +if options.data_source == "lvshm": + demux = pipeparts.mkframecppchanneldemux(pipeline, src, do_file_checksum = options.do_file_checksum, skip_bad_files = True, channel_list = map("%s:%s".__mod__, channel_list)) + +elif options.data_source == "frames": + demux = pipeparts.mkframecppchanneldemux(pipeline, src, do_file_checksum = options.do_file_checksum, skip_bad_files = False, channel_list = map("%s:%s".__mod__, channel_list)) + +# Write the pipeline graph after pads have been hooked up to the demuxer +if options.write_pipeline is not None: + demux.connect("no-more-pads", write_graph) + +# Get everything hooked up and fill in discontinuities +for key, chan in zip(headkeys, channel_list): + head_dict[key] = calibration_parts.hook_up(pipeline, demux, chan[1], instrument, options.buffer_length) + +# +# TIME-VARYING FACTORS COMPUTATIONS +# + +# Get all the channels we need first and make sure they have the desired sample rate and caps +darm_err = calibration_parts.caps_and_progress(pipeline, head_dict["darm_err"], derr_caps, "darm_err") +darm_err = calibration_parts.mkresample(pipeline, darm_err, 5, False, int(kappas_filter_sr)) + +# Apply a filter to DARM_ERR to simulate time-dependence (if the filter is not provided, is is set to [1]) +darm_err = pipeparts.mkfirbank(pipeline, darm_err, latency = int(kappas_filter_delay), fir_matrix = [derr_kappas_filter[::-1]], time_domain = True) +derrtee = pipeparts.mktee(pipeline, darm_err) + +# Get the injection channels +pcal = calibration_parts.caps_and_progress(pipeline, head_dict["pcal"], derr_caps, "pcal") +pcal = calibration_parts.mkresample(pipeline, pcal, 5, False, int(kappas_filter_sr)) +pcaltee = pipeparts.mktee(pipeline, pcal) +tstexccaps = "audio/x-raw, format=F64LE, rate=%d" % options.tst_exc_sample_rate +tstexc = calibration_parts.caps_and_progress(pipeline, head_dict["tstexc"], tstexccaps, "tstexc") + +if not options.no_kappac or not options.no_fcc or not options.no_kappapu or not options.no_kappap or not options.no_kappau or not options.no_fs or not options.no_srcQ: + exc = calibration_parts.caps_and_progress(pipeline, head_dict["exc"], derr_caps, "exc") + exc = calibration_parts.mkresample(pipeline, exc, 5, False, int(kappas_filter_sr)) + +# Get the EPICS channels +for key in headkeys: + if key.startswith("EP"): + head_dict[key] = calibration_parts.caps_and_progress(pipeline, head_dict[key], ref_factors_caps, key) + head_dict[key] = calibration_parts.mkresample(pipeline, head_dict[key], 0, False, compute_calib_factors_caps) + head_dict[key] = pipeparts.mktee(pipeline, head_dict[key]) + +# Get the coherence channels +if not options.no_coherence: + if not options.no_kappatst or not options.no_kappapu or not options.no_kappac or not options.no_fcc or not options.no_fs or not options.no_srcQ: + pcaly_line1_coh = calibration_parts.caps_and_progress(pipeline, head_dict["pcaly_line1_coh"], coh_caps, "pcaly_line1_coh") + pcaly_line1_coh = calibration_parts.mkresample(pipeline, pcaly_line1_coh, 0, False, compute_calib_factors_caps) + sus_coh = calibration_parts.caps_and_progress(pipeline, head_dict["sus_coh"], coh_caps, "sus_coh") + sus_coh = calibration_parts.mkresample(pipeline, sus_coh, 0, False, compute_calib_factors_caps) + darm_coh = calibration_parts.caps_and_progress(pipeline, head_dict["darm_coh"], coh_caps, "darm_coh") + darm_coh = calibration_parts.mkresample(pipeline, darm_coh, 0, False, compute_calib_factors_caps) + pcaly_line1_coh = pipeparts.mktee(pipeline, pcaly_line1_coh) + sus_coh = pipeparts.mktee(pipeline, sus_coh) + darm_coh = pipeparts.mktee(pipeline, darm_coh) + if not options.no_kappac or not options.no_fcc or not options.no_fs or not options.no_srcQ: + pcaly_line2_coh = calibration_parts.caps_and_progress(pipeline, head_dict["pcaly_line2_coh"], coh_caps, "pcaly_line2_coh") + pcaly_line2_coh = calibration_parts.mkresample(pipeline, pcaly_line2_coh, 0, False, compute_calib_factors_caps) + pcaly_line2_coh = pipeparts.mktee(pipeline, pcaly_line2_coh) + +# Set up computations for \kappa_tst,\kappa_c, \kappa_pu, f_cc, if applicable + +# demodulate the PCAL channel and apply the PCAL correction factor at the DARM actuation line frequency +pcal_at_darm_act_freq = calibration_parts.demodulate(pipeline, pcaltee, darm_act_line_freq, True, options.compute_factors_sr, options.demodulation_filter_time, options.filter_latency, prefactor_real = pcal_corr_at_darm_act_freq_real, prefactor_imag = pcal_corr_at_darm_act_freq_imag) +if not options.no_kappapu or not options.no_kappac or not options.no_fcc or not options.no_srcQ or not options.no_fs: + pcal_at_darm_act_freq = pipeparts.mktee(pipeline, pcal_at_darm_act_freq) + +# demodulate DARM_ERR at the DARM actuation line frequency +derr_at_darm_act_freq = calibration_parts.demodulate(pipeline, derrtee, darm_act_line_freq, True, options.compute_factors_sr, options.demodulation_filter_time, options.filter_latency) +if not options.no_kappapu or not options.no_kappac or not options.no_fcc or not options.no_srcQ or not options.no_fs: + 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, True, options.compute_factors_sr, options.demodulation_filter_time, options.filter_latency) + +# demodulate DARM_ERR at the ESD actuation line frequency +derr_at_esd_act_freq = calibration_parts.demodulate(pipeline, derrtee, esd_act_line_freq, True, options.compute_factors_sr, options.demodulation_filter_time, options.filter_latency) + +# compute kappa_tst, either using reference factors from the filters file or reading them from EPICS channels +if not options.factors_from_filters_file: + EP1 = calibration_parts.merge_into_complex(pipeline, head_dict["EP1_real"], head_dict["EP1_imag"]) + ktst = calibration_parts.compute_kappatst(pipeline, derr_at_esd_act_freq, tstexc_at_esd_act_freq, pcal_at_darm_act_freq, derr_at_darm_act_freq, EP1) +elif options.factors_from_filters_file: + ktst = calibration_parts.compute_kappatst_from_filters_file(pipeline, derr_at_esd_act_freq, tstexc_at_esd_act_freq, pcal_at_darm_act_freq, derr_at_darm_act_freq, EP1_real, EP1_imag) + +ktst = pipeparts.mktee(pipeline, ktst) + +# Put off smoothing \kappa_tst until after \kappa_pu is computed in case we are correcting the phase of \kappa_tst using \kappa_pu + +# If we're also computing \kappa_pu, \kappa_c, f_cc, fs, or Q, keep going +if not options.no_kappac or not options.no_fcc or not options.no_kappapu or not options.np_kappap or not options.no_kappau or not options.no_srcQ or not options.no_fs or (not options.no_kappatst and options.act_timing_from_kappapu): + # demodulate excitation channel at PU actuation line frequency + exc_at_pu_act_freq = calibration_parts.demodulate(pipeline, exc, pu_act_esd_line_freq, True, options.compute_factors_sr, options.demodulation_filter_time, options.filter_latency) + + # demodulate DARM_ERR at PU actuation line frequency + derr_at_pu_act_freq = calibration_parts.demodulate(pipeline, derrtee, pu_act_esd_line_freq, True, options.compute_factors_sr, options.demodulation_filter_time, options.filter_latency) + + # compute the factor Afctrl that will be used in the computation of kappa_pu and kappa_a, either using reference factors from the filters file or reading them from EPICS channels + if not options.factors_from_filters_file: + EP2 = calibration_parts.merge_into_complex(pipeline, head_dict["EP2_real"], head_dict["EP2_imag"]) + EP3 = calibration_parts.merge_into_complex(pipeline, head_dict["EP3_real"], head_dict["EP3_imag"]) + EP4 = calibration_parts.merge_into_complex(pipeline, head_dict["EP4_real"], head_dict["EP4_imag"]) + afctrl = calibration_parts.compute_afctrl(pipeline, derr_at_pu_act_freq, exc_at_pu_act_freq, pcal_at_darm_act_freq, derr_at_darm_act_freq, EP2) + elif options.factors_from_filters_file: + afctrl = calibration_parts.compute_afctrl_from_filters_file(pipeline, derr_at_pu_act_freq, exc_at_pu_act_freq, pcal_at_darm_act_freq, derr_at_darm_act_freq, EP2_real, EP2_imag) + + # \kappa_pu calcuation, which needs to happen for any of the other kappas to be computed + if not options.factors_from_filters_file: + kpu = calibration_parts.compute_kappapu(pipeline, EP3, afctrl, ktst, EP4) + elif options.factors_from_filters_file: + kpu = calibration_parts.compute_kappapu_from_filters_file(pipeline, EP3_real, EP3_imag, afctrl, ktst, EP4_real, EP4_imag) + + kpu = pipeparts.mktee(pipeline, kpu) + + # Put off smoothing \kappa_pu in case we are correcting the phase of \kappa_pu using \kappa_tst + + # If desired, correct the phase of \kappa_pu using \kappa_tst (This assumes that all stages of actuation have the same variable time delay, and that \kappa_tst is doing a better job of measuring it) + if options.act_timing_from_kappatst: + # Find the magnitude of \kappa_pu + kpu = pipeparts.mkgeneric(pipeline, kpu, "cabs") + # Find the phase of \kappa_tst + phi_ktst = pipeparts.mkgeneric(pipeline, ktst, "carg") + # Multiply by the line-frequency ratio to get the phase of \kappa_pu + phi_kpu = pipeparts.mkaudioamplify(pipeline, phi_ktst, pu_act_esd_line_freq / esd_act_line_freq) + # Find the phase factor + kpu_phase_factor = pipeparts.mkgeneric(pipeline, pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, phi_kpu, matrix = [[0.0, 1.0]])), "cexp") + # Multiply by the magnitude of \kappa_pu + kpu = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, kpu, matrix = [[1.0, 0.0]])), kpu_phase_factor)) + kpu = pipeparts.mktee(pipeline, kpu) + + # If desired, correct the phase of \kappa_tst using \kappa_pu (This assumes that all stages of actuation have the same variable time delay, and that \kappa_pu is doing a better job of measuring it) + if options.act_timing_from_kappapu: + # Find the magnitude of \kappa_tst + ktst = pipeparts.mkgeneric(pipeline, ktst, "cabs") + # Find the phase of \kappa_tst + phi_kpu = pipeparts.mkgeneric(pipeline, kpu, "carg") + # Multiply by the line-frequency ratio to get the phase of \kappa_tst + phi_ktst = pipeparts.mkaudioamplify(pipeline, phi_kpu, esd_act_line_freq / pu_act_esd_line_freq) + # Find the phase factor + ktst_phase_factor = pipeparts.mkgeneric(pipeline, pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, phi_ktst, matrix = [[0.0, 1.0]])), "cexp") + # Multiply by the magnitude of \kappa_tst + ktst = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, ktst, matrix = [[1.0, 0.0]])), ktst_phase_factor)) + ktst = pipeparts.mktee(pipeline, ktst) + + # Now apply the gating and smoothing to \kappa_tst and \kappa_pu + if not options.no_kappatst: + smooth_ktst_nogate = pipeparts.mkgeneric(pipeline, ktst, "lal_smoothkappas", default_kappa_re = options.expected_kappatst_real, default_kappa_im = options.expected_kappatst_imag, array_size = median_smoothing_samples, avg_array_size = factors_average_samples, default_to_median = options.kappas_default_to_median, filter_latency = options.filter_latency) + smooth_ktstR_nogate, smooth_ktstI_nogate = calibration_parts.split_into_real(pipeline, smooth_ktst_nogate) + + if not options.no_coherence: + # Gate kappa_tst with the coherence of the PCALY_line1 line + ktst_gated = calibration_parts.mkgate(pipeline, ktst, pcaly_line1_coh, options.coherence_uncertainty_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True) + # Gate kappa_tst with the coherence of the suspension line + ktst_gated = calibration_parts.mkgate(pipeline, ktst_gated, sus_coh, options.coherence_uncertainty_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True) + # Gate kappa_tst with the coherence of the DARM line + ktst_gated = calibration_parts.mkgate(pipeline, ktst_gated, darm_coh, options.coherence_uncertainty_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True) + + # Smooth kappa_tst + smooth_ktst = calibration_parts.smooth_complex_kappas(pipeline, ktst_gated, options.expected_kappatst_real, options.expected_kappatst_imag, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median, options.filter_latency) + + else: + # Smooth kappa_tst + smooth_ktst = calibration_parts.smooth_complex_kappas_no_coherence(pipeline, ktst, options.kappatst_real_ok_var, options.kappatst_imag_ok_var, options.expected_kappatst_real, options.expected_kappatst_imag, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median, options.filter_latency) + + smooth_ktsttee = pipeparts.mktee(pipeline, smooth_ktst) + smooth_ktstR, smooth_ktstI = calibration_parts.split_into_real(pipeline, smooth_ktsttee) + + smooth_ktstRtee = pipeparts.mktee(pipeline, smooth_ktstR) + smooth_ktstItee = pipeparts.mktee(pipeline, smooth_ktstI) + + + + if not options.no_kappapu: + smooth_kpu_nogate = pipeparts.mkgeneric(pipeline, kpu, "lal_smoothkappas", default_kappa_re = options.expected_kappapu_real, default_kappa_im = options.expected_kappapu_imag, array_size = median_smoothing_samples, avg_array_size = factors_average_samples, default_to_median = options.kappas_default_to_median, filter_latency = options.filter_latency) + smooth_kpuR_nogate, smooth_kpuI_nogate = calibration_parts.split_into_real(pipeline, smooth_kpu_nogate) + + if not options.no_coherence: + # Gate kappa_pu with the coherence of the DARM line + kpu_gated = calibration_parts.mkgate(pipeline, kpu, darm_coh, options.coherence_uncertainty_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True) + # Gate kappa_pu with the coherence of the PCALY_line1 line + kpu_gated = calibration_parts.mkgate(pipeline, kpu_gated, pcaly_line1_coh, options.coherence_uncertainty_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True) + # Gate kappa_pu with the coherence of the suspension coherence + kpu_gated = calibration_parts.mkgate(pipeline, kpu_gated, sus_coh, options.coherence_uncertainty_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True) + + # Smooth kappa_pu + smooth_kpu = calibration_parts.smooth_complex_kappas(pipeline, kpu_gated, options.expected_kappapu_real, options.expected_kappapu_imag, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median, options.filter_latency) + + else: + # Smooth kappa_pu + smooth_kpu = calibration_parts.smooth_complex_kappas_no_coherence(pipeline, kpu, options.kappapu_real_ok_var, options.kappapu_imag_ok_var, options.expected_kappapu_real, options.expected_kappapu_imag, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median, options.filter_latency) + + smooth_kputee = pipeparts.mktee(pipeline, smooth_kpu) + smooth_kpuR, smooth_kpuI = calibration_parts.split_into_real(pipeline, smooth_kputee) + + smooth_kpuRtee = pipeparts.mktee(pipeline, smooth_kpuR) + smooth_kpuItee = pipeparts.mktee(pipeline, smooth_kpuI) + + # Finally, compute \kappa_c and f_cc + if not options.no_kappac or not options.no_fcc or not options.no_srcQ or not options.no_fs: + # demodulate the PCAL channel and apply the PCAL correction factor at optical gain and f_cc line frequency + pcal_at_opt_gain_freq = calibration_parts.demodulate(pipeline, pcaltee, opt_gain_fcc_line_freq, True, options.compute_factors_sr, options.demodulation_filter_time, options.filter_latency, prefactor_real = pcal_corr_at_opt_gain_fcc_freq_real, prefactor_imag = pcal_corr_at_opt_gain_fcc_freq_imag) + + # demodulate DARM_ERR at optical gain and f_cc line frequency + derr_at_opt_gain_freq = calibration_parts.demodulate(pipeline, derrtee, opt_gain_fcc_line_freq, True, options.compute_factors_sr, options.demodulation_filter_time, options.filter_latency) + + # Compute the factor S which will be used for the kappa_c and f_cc calculations + # \kappa_tst and \kappa_pu need to be evaluated at the higher pcal line frequency + ktst_at_opt_gain_freq = pipeparts.mkgeneric(pipeline, ktst, "lpshiftfreq", frequency_ratio = opt_gain_fcc_line_freq / esd_act_line_freq) + kpu_at_opt_gain_freq = pipeparts.mkgeneric(pipeline, kpu, "lpshiftfreq", frequency_ratio = opt_gain_fcc_line_freq / pu_act_esd_line_freq) + if not options.factors_from_filters_file: + EP6 = calibration_parts.merge_into_complex(pipeline, head_dict["EP6_real"], head_dict["EP6_imag"]) + EP7 = calibration_parts.merge_into_complex(pipeline, head_dict["EP7_real"], head_dict["EP7_imag"]) + EP8 = calibration_parts.merge_into_complex(pipeline, head_dict["EP8_real"], head_dict["EP8_imag"]) + EP9 = calibration_parts.merge_into_complex(pipeline, head_dict["EP9_real"], head_dict["EP9_imag"]) + S = calibration_parts.compute_S(pipeline, EP6, pcal_at_opt_gain_freq, derr_at_opt_gain_freq, EP7, ktst_at_opt_gain_freq, EP8, kpu_at_opt_gain_freq, EP9) + elif options.factors_from_filters_file: + S = calibration_parts.compute_S_from_filters_file(pipeline, EP6_real, EP6_imag, pcal_at_opt_gain_freq, derr_at_opt_gain_freq, EP7_real, EP7_imag, ktst_at_opt_gain_freq, EP8_real, EP8_imag, kpu_at_opt_gain_freq, EP9_real, EP9_imag) + + S = pipeparts.mktee(pipeline, S) + + SR, SI = calibration_parts.split_into_real(pipeline, S) + + if not options.no_kappac and not options.no_fcc: + SR = pipeparts.mktee(pipeline, SR) + SI = pipeparts.mktee(pipeline, SI) + + # compute kappa_c + if not options.no_kappac or not options.no_srcQ or not options.no_fs: + kc = calibration_parts.compute_kappac(pipeline, SR, SI) + if not options.no_kappac: + kc = pipeparts.mktee(pipeline, kc) + smooth_kc_nogate = pipeparts.mkgeneric(pipeline, kc, "lal_smoothkappas", default_kappa_re = options.expected_kappac, array_size = median_smoothing_samples, avg_array_size = factors_average_samples, default_to_median = options.kappas_default_to_median, filter_latency = options.filter_latency) + + if not options.no_coherence: + # Gate kappa_c with the coherence of all four of the calibration lines + kc_gated = calibration_parts.mkgate(pipeline, kc, pcaly_line2_coh, options.coherence_uncertainty_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True) + kc_gated = calibration_parts.mkgate(pipeline, kc_gated, darm_coh, options.coherence_uncertainty_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True) + kc_gated = calibration_parts.mkgate(pipeline, kc_gated, pcaly_line1_coh, options.coherence_uncertainty_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True) + kc_gated = calibration_parts.mkgate(pipeline, kc_gated, sus_coh, options.coherence_uncertainty_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True) + + # Smooth kappa_c + smooth_kc = calibration_parts.smooth_kappas(pipeline, kc_gated, options.expected_kappac, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median, options.filter_latency) + + else: + # Smooth kappa_c + smooth_kc = calibration_parts.smooth_kappas_no_coherence(pipeline, kc, options.kappac_ok_var, options.expected_kappac, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median, options.filter_latency) + + smooth_kctee = pipeparts.mktee(pipeline, smooth_kc) + + # compute f_cc + if not options.no_fcc or not options.no_srcQ or not options.no_fs: + fcc = calibration_parts.compute_fcc(pipeline, SR, SI, opt_gain_fcc_line_freq) + if not options.no_fcc: + fcc = pipeparts.mktee(pipeline, fcc) + smooth_fcc_nogate = pipeparts.mkgeneric(pipeline, fcc, "lal_smoothkappas", default_kappa_re = fcc_default, array_size = median_smoothing_samples, avg_array_size = factors_average_samples, default_to_median = options.kappas_default_to_median, filter_latency = options.filter_latency) + + if not options.no_coherence: + # Gate f_cc with all four of the calibration lines + fcc_gated = calibration_parts.mkgate(pipeline, fcc, pcaly_line2_coh, options.coherence_uncertainty_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True) + fcc_gated = calibration_parts.mkgate(pipeline, fcc_gated, darm_coh, options.coherence_uncertainty_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True) + fcc_gated = calibration_parts.mkgate(pipeline, fcc_gated, pcaly_line1_coh, options.coherence_uncertainty_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True) + fcc_gated = calibration_parts.mkgate(pipeline, fcc_gated, sus_coh, options.coherence_uncertainty_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True) + + # Smooth f_cc + smooth_fcc = calibration_parts.smooth_kappas(pipeline, fcc_gated, fcc_default, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median, options.filter_latency) + else: + # Smooth f_cc + smooth_fcc = calibration_parts.smooth_kappas_no_coherence(pipeline, fcc, options.fcc_ok_var, fcc_default, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median, options.filter_latency) + + smooth_fcctee = pipeparts.mktee(pipeline, smooth_fcc) + + # compute f_s and Q + if not options.no_fs or not options.no_srcQ: + expected_Xi = complex((fs_default * fs_default - 1j * src_pcal_line_freq * fs_default / srcQ_default) / (src_pcal_line_freq * src_pcal_line_freq)) + Xi_real_ok_var = float((pow(fs_default + options.fs_ok_var, 2) - pow(fs_default, 2.0)) / pow(src_pcal_line_freq, 2)) + Xi_imag_ok_var = float(options.fs_ok_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, True, options.compute_factors_sr, options.demodulation_filter_time, options.filter_latency, prefactor_real = pcal_corr_at_src_freq_real, prefactor_imag = pcal_corr_at_src_freq_imag) + + # demodulate DARM_ERR at SRC detuning line frequency + derr_at_src_freq = calibration_parts.demodulate(pipeline, derrtee, src_pcal_line_freq, True, options.compute_factors_sr, options.demodulation_filter_time, options.filter_latency) + + # 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 + ktst_at_src_freq = pipeparts.mkgeneric(pipeline, ktst, "lpshiftfreq", frequency_ratio = src_pcal_line_freq / esd_act_line_freq) + kpu_at_src_freq = pipeparts.mkgeneric(pipeline, kpu, "lpshiftfreq", frequency_ratio = src_pcal_line_freq / pu_act_esd_line_freq) + if not options.factors_from_filters_file: + EP11 = calibration_parts.merge_into_complex(pipeline, head_dict["EP11_real"], head_dict["EP11_imag"]) + EP12 = calibration_parts.merge_into_complex(pipeline, head_dict["EP12_real"], head_dict["EP12_imag"]) + EP13 = calibration_parts.merge_into_complex(pipeline, head_dict["EP13_real"], head_dict["EP13_imag"]) + EP14 = calibration_parts.merge_into_complex(pipeline, head_dict["EP14_real"], head_dict["EP14_imag"]) + Xi = calibration_parts.compute_Xi(pipeline, pcal_at_src_freq, derr_at_src_freq, src_pcal_line_freq, EP11, EP12, EP13, EP14, ktst_at_src_freq, kpu_at_src_freq, kc, fcc) + elif options.factors_from_filters_file: + Xi = calibration_parts.compute_Xi_from_filters_file(pipeline, pcal_at_src_freq, derr_at_src_freq, src_pcal_line_freq, EP11_real, EP11_imag, EP12_real, EP12_imag, EP13_real, EP13_imag, EP14_real, EP14_imag, ktst_at_src_freq, kpu_at_src_freq, kc, fcc) + + Xi = pipeparts.mktee(pipeline, Xi) + smooth_Xi_nogate = pipeparts.mkgeneric(pipeline, Xi, "lal_smoothkappas", default_kappa_re = float(numpy.real(expected_Xi)), default_kappa_im = float(numpy.imag(expected_Xi)), array_size = median_smoothing_samples, avg_array_size = src_average_samples, default_to_median = options.kappas_default_to_median, filter_latency = options.filter_latency) + + if not options.no_coherence: + # Gate Xi with all coherences. We apply the gating and smoothing here since Q depends on the inverse of Im(Xi), which fluctuates about zero. + Xi_gated = calibration_parts.mkgate(pipeline, Xi, pcaly_line2_coh, options.coherence_uncertainty_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True) + Xi_gated = calibration_parts.mkgate(pipeline, Xi_gated, darm_coh, options.coherence_uncertainty_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True) + Xi_gated = calibration_parts.mkgate(pipeline, Xi_gated, pcaly_line2_coh, options.coherence_uncertainty_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True) + Xi_gated = calibration_parts.mkgate(pipeline, Xi_gated, sus_coh, options.coherence_uncertainty_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True) + + smooth_Xi = calibration_parts.smooth_complex_kappas(pipeline, Xi_gated, float(numpy.real(expected_Xi)), float(numpy.imag(expected_Xi)), median_smoothing_samples, src_average_samples, options.kappas_default_to_median, options.filter_latency) + + else: + smooth_Xi = calibration_parts.smooth_complex_kappas_no_coherence(pipeline, Xi, Xi_real_ok_var, Xi_real_ok_var, float(numpy.real(expected_Xi)), float(numpy.imag(expected_Xi)), median_smoothing_samples, src_average_samples, options.kappas_default_to_median, options.filter_latency) + + if options.no_srcQ: + # the imaginary part is only used to compute Q + smooth_XiR = pipeparts.mkgeneric(pipeline, smooth_Xi, "creal") + smooth_XiR_nogate = pipeparts.mkgeneric(pipeline, smooth_Xi_nogate, "creal") + else: + smooth_XiR, smooth_XiI = calibration_parts.split_into_real(pipeline, smooth_Xi) + smooth_XiR_nogate, smooth_XiI_nogate = calibration_parts.split_into_real(pipeline, smooth_Xi_nogate) + + smooth_sqrtXiR = calibration_parts.mkpow(pipeline, smooth_XiR, exponent = 0.5) + smooth_sqrtXiR_nogate = calibration_parts.mkpow(pipeline, smooth_XiR_nogate, exponent = 0.5) + + if not options.no_fs and not options.no_srcQ: + smooth_sqrtXiR = pipeparts.mktee(pipeline, smooth_sqrtXiR) + smooth_sqrtXiR_nogate = pipeparts.mktee(pipeline, smooth_sqrtXiR_nogate) + + # compute f_s + if not options.no_fs: + smooth_fs = pipeparts.mkaudioamplify(pipeline, smooth_sqrtXiR, src_pcal_line_freq) + smooth_fs_nogate = pipeparts.mkaudioamplify(pipeline, smooth_sqrtXiR_nogate, src_pcal_line_freq) + + if not options.no_dq_vector: + smooth_fs = pipeparts.mktee(pipeline, smooth_fs) + + # compute SRC Q_inv + if not options.no_srcQ: + smooth_sqrtXiR_inv = calibration_parts.mkpow(pipeline, smooth_sqrtXiR, exponent = -1.0) + smooth_sqrtXiR_inv_nogate = calibration_parts.mkpow(pipeline, smooth_sqrtXiR_nogate, exponent = -1.0) + smooth_srcQ_inv = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, smooth_sqrtXiR_inv, pipeparts.mkaudioamplify(pipeline, smooth_XiI, -1.0))) + smooth_srcQ_inv_nogate = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, smooth_sqrtXiR_inv_nogate, pipeparts.mkaudioamplify(pipeline, smooth_XiI_nogate, -1.0))) + + if not options.no_dq_vector: + smooth_srcQ_inv = pipeparts.mktee(pipeline, smooth_srcQ_inv) + + +# +# KAPPA_STATE_VECTOR BRANCH +# + +#FIXME: Add more comments! + +if not options.no_dq_vector: + + # + # KAPPATST BITS BRANCH + # + if not options.no_kappatst: + ktstSmoothInRange = calibration_parts.compute_kappa_bits(pipeline, smooth_ktstRtee, smooth_ktstItee, options.expected_kappatst_real, options.expected_kappatst_imag, options.kappatst_real_ok_var, options.kappatst_imag_ok_var, int(median_smoothing_samples / 2) + factors_average_samples, status_out_smooth = pow(2,9), starting_rate = options.compute_factors_sr, ending_rate = kappasstatesr) + + # + # KAPPAP BITS BRANCH + # + if not options.no_kappap: + kpSmoothInRange = calibration_parts.compute_kappa_bits(pipeline, smooth_kpRtee, smooth_kpItee, options.expected_kappap_real, options.expected_kappap_imag, options.kappap_real_ok_var, options.kappap_imag_ok_var, int(median_smoothing_samples / 2) + factors_average_samples, status_out_smooth = pow(2,10), starting_rate = options.compute_factors_sr, ending_rate = kappasstatesr) + + # + # KAPPAPU BITS BRANCH + # + elif not options.no_kappapu: + kpSmoothInRange = calibration_parts.compute_kappa_bits(pipeline, smooth_kpuRtee, smooth_kpuItee, options.expected_kappapu_real, options.expected_kappapu_imag, options.kappapu_real_ok_var, options.kappapu_imag_ok_var, int(median_smoothing_samples / 2) + factors_average_samples, status_out_smooth = pow(2,10), starting_rate = options.compute_factors_sr, ending_rate = kappasstatesr) + + # + # KAPPAU BITS BRANCH + # + if not options.no_kappau: + kuSmoothInRange = calibration_parts.compute_kappa_bits(pipeline, smooth_kuRtee, smooth_kuItee, options.expected_kappau_real, options.expected_kappau_imag, options.kappau_real_ok_var, options.kappau_imag_ok_var, int(median_smoothing_samples / 2) + factors_average_samples, status_out_smooth = pow(2,11), starting_rate = options.compute_factors_sr, ending_rate = kappasstatesr) + + # + # KAPPAC BITS BRANCH + # + if not options.no_kappac: + kcSmoothInRange = calibration_parts.compute_kappa_bits_only_real(pipeline, smooth_kctee, options.expected_kappac, options.kappac_ok_var, int(median_smoothing_samples / 2) + factors_average_samples, status_out_smooth = pow(2,12), starting_rate = options.compute_factors_sr, ending_rate = kappasstatesr) + + # + # FCC BITS BRANCH + # + if not options.no_fcc: + fccSmoothInRange = calibration_parts.compute_kappa_bits_only_real(pipeline, smooth_fcctee, fcc_default, options.fcc_ok_var, int(median_smoothing_samples / 2) + factors_average_samples, status_out_smooth = pow(2,13), starting_rate = options.compute_factors_sr, ending_rate = kappasstatesr) + + # + # FS BITS BRANCH + # + if not options.no_fs: + fsSmoothInRange = calibration_parts.compute_kappa_bits_only_real(pipeline, smooth_fs, fs_default, options.fs_ok_var, int(median_smoothing_samples / 2) + src_average_samples, status_out_smooth = pow(2,14), starting_rate = options.compute_factors_sr, ending_rate = kappasstatesr) + + # + # SRCQ BITS BRANCH + # + if not options.no_srcQ: + srcQSmoothInRange = calibration_parts.compute_kappa_bits_only_real(pipeline, smooth_srcQ_inv, 1.0 / srcQ_default, [options.srcQinv_min, options.srcQinv_max], int(median_smoothing_samples / 2) + src_average_samples, status_out_smooth = pow(2,15), starting_rate = options.compute_factors_sr, ending_rate = kappasstatesr) + + # + # COHERENCE BITS BRANCH + # + if not options.no_coherence: + if not options.no_kappatst or not options.no_kappapu or not options.no_kappac or not options.no_fcc or not options.no_srcQ or not options.no_fs: + pcaly_line1_coh_ok = pipeparts.mkbitvectorgen(pipeline, pcaly_line1_coh, threshold = options.coherence_uncertainty_threshold, bit_vector = pow(2,18), invert_control = True) + pcaly_line1_coh_ok = pipeparts.mkcapsfilter(pipeline, pcaly_line1_coh_ok, "audio/x-raw, format=U32LE, rate=%d" % cohsr) + pcaly_line1_coh_ok = pipeparts.mkgeneric(pipeline, pcaly_line1_coh_ok, "lal_logicalundersample", required_on = pow(2,18), status_out = pow(2,18)) + pcaly_line1_coh_ok = pipeparts.mkcapsfilter(pipeline, pcaly_line1_coh_ok, calibstate_caps) + + sus_coh_ok = pipeparts.mkbitvectorgen(pipeline, sus_coh, threshold = options.coherence_uncertainty_threshold, bit_vector = pow(2,16), invert_control = True) + sus_coh_ok = pipeparts.mkcapsfilter(pipeline, sus_coh_ok, "audio/x-raw, format=U32LE, rate=%d" % cohsr) + sus_coh_ok = pipeparts.mkgeneric(pipeline, sus_coh_ok, "lal_logicalundersample", required_on = pow(2,16), status_out = pow(2,16)) + sus_coh_ok = pipeparts.mkcapsfilter(pipeline, sus_coh_ok, calibstate_caps) + + darm_coh_ok = pipeparts.mkbitvectorgen(pipeline, darm_coh, threshold = options.coherence_uncertainty_threshold, bit_vector = pow(2,17), invert_control = True) + darm_coh_ok = pipeparts.mkcapsfilter(pipeline, darm_coh_ok, "audio/x-raw, format=U32LE, rate=%d" % cohsr) + darm_coh_ok = pipeparts.mkgeneric(pipeline, darm_coh_ok, "lal_logicalundersample", required_on = pow(2,17), status_out = pow(2,17)) + darm_coh_ok = pipeparts.mkcapsfilter(pipeline, darm_coh_ok, calibstate_caps) + coherence_bits = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, pcaly_line1_coh_ok, sus_coh_ok, darm_coh_ok)) + if not options.no_kappac or not options.no_fcc or not options.no_srcQ or not options.no_fs: + pcaly_line2_coh_ok = pipeparts.mkbitvectorgen(pipeline, pcaly_line2_coh, threshold = options.coherence_uncertainty_threshold, bit_vector = pow(2,19), invert_control = True) + pcaly_line2_coh_ok = pipeparts.mkcapsfilter(pipeline, pcaly_line2_coh_ok, "audio/x-raw, format=U32LE, rate=%d" % cohsr) + pcaly_line2_coh_ok = pipeparts.mkgeneric(pipeline, pcaly_line2_coh_ok, "lal_logicalundersample", required_on = pow(2,19), status_out = pow(2,19)) + pcaly_line2_coh_ok = pipeparts.mkcapsfilter(pipeline, pcaly_line2_coh_ok, calibstate_caps) + coherence_bits = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, coherence_bits, pcaly_line2_coh_ok)) + + # + # EPICS BITS + # + + D_epics_threshold = 0.0 + A_epics_threshold = 0.0 + C_epics_threshold = 0.0 + other_epics_threshold = 0.0 + + # First, check the EPICS that involve only the digital filter D, EP7 and EP12 + if num_dq_epics > 6: + D_epics_threshold += 1.5 + EP7_real_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP7_real"], 1.0 / EP7_real) + EP7_imag_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP7_imag"], 1.0 / EP7_imag) + + # The above values should be close to one (within 1 / 10^4) + EP7_real_check = pipeparts.mkgeneric(pipeline, EP7_real_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + EP7_imag_check = pipeparts.mkgeneric(pipeline, EP7_imag_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + D_epics_check = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, EP7_real_check, EP7_imag_check)) + + if num_dq_epics > 11: + D_epics_threshold += 2.0 + EP12_real_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP12_real"], 1.0 / EP12_real) + EP12_imag_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP12_imag"], 1.0 / EP12_imag) + + # The above values should be close to one (within 1 / 10^4) + EP12_real_check = pipeparts.mkgeneric(pipeline, EP12_real_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + EP12_imag_check = pipeparts.mkgeneric(pipeline, EP12_imag_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + D_epics_check = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, D_epics_check, EP12_real_check, EP12_imag_check)) + + D_epics_bit = pipeparts.mkbitvectorgen(pipeline, D_epics_check, bit_vector = pow(2,21), threshold = D_epics_threshold) + D_epics_bit = pipeparts.mkgeneric(pipeline, D_epics_bit, "lal_logicalundersample", required_on = pow(2,21), status_out = pow(2,21)) + D_epics_bit = pipeparts.mkcapsfilter(pipeline, D_epics_bit, calibstate_caps) + + # Next, check the EPICS that involve only the actuation function A, EP3, EP4, EP8, EP9, EP10, EP13, EP14 + if num_dq_epics > 2: + A_epics_threshold += 1.5 + EP3_real_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP3_real"], 1.0 / EP3_real) + EP3_imag_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP3_imag"], 1.0 / EP3_imag) + + # The above values should be close to one (within 1 / 10^4) + EP3_real_check = pipeparts.mkgeneric(pipeline, EP3_real_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + EP3_imag_check = pipeparts.mkgeneric(pipeline, EP3_imag_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + A_epics_check = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, EP3_real_check, EP3_imag_check)) + + if num_dq_epics > 3: + A_epics_threshold += 2.0 + EP4_real_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP4_real"], 1.0 / EP4_real) + EP4_imag_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP4_imag"], 1.0 / EP4_imag) + + # The above values should be close to one (within 1 / 10^4) + EP4_real_check = pipeparts.mkgeneric(pipeline, EP4_real_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + EP4_imag_check = pipeparts.mkgeneric(pipeline, EP4_imag_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + A_epics_check = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, A_epics_check, EP4_real_check, EP4_imag_check)) + + if num_dq_epics > 7: + A_epics_threshold += 2.0 + EP8_real_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP8_real"], 1.0 / EP8_real) + EP8_imag_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP8_imag"], 1.0 / EP8_imag) + + # The above values should be close to one (within 1 / 10^4) + EP8_real_check = pipeparts.mkgeneric(pipeline, EP8_real_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + EP8_imag_check = pipeparts.mkgeneric(pipeline, EP8_imag_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + A_epics_check = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, A_epics_check, EP8_real_check, EP8_imag_check)) + + if num_dq_epics > 8: + A_epics_threshold += 2.0 + EP9_real_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP9_real"], 1.0 / EP9_real) + EP9_imag_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP9_imag"], 1.0 / EP9_imag) + + # The above values should be close to one (within 1 / 10^4) + EP9_real_check = pipeparts.mkgeneric(pipeline, EP9_real_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + EP9_imag_check = pipeparts.mkgeneric(pipeline, EP9_imag_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + A_epics_check = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, A_epics_check, EP9_real_check, EP9_imag_check)) + + if num_dq_epics > 9: + A_epics_threshold += 2.0 + EP10_real_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP10_real"], 1.0 / EP10_real) + EP10_imag_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP10_imag"], 1.0 / EP10_imag) + + # The above values should be close to one (within 1 / 10^4) + EP10_real_check = pipeparts.mkgeneric(pipeline, EP10_real_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + EP10_imag_check = pipeparts.mkgeneric(pipeline, EP10_imag_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + A_epics_check = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, A_epics_check, EP10_real_check, EP10_imag_check)) + + if num_dq_epics > 12: + A_epics_threshold += 2.0 + EP13_real_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP13_real"], 1.0 / EP13_real) + EP13_imag_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP13_imag"], 1.0 / EP13_imag) + + # The above values should be close to one (within 1 / 10^4) + EP13_real_check = pipeparts.mkgeneric(pipeline, EP13_real_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + EP13_imag_check = pipeparts.mkgeneric(pipeline, EP13_imag_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + A_epics_check = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, A_epics_check, EP13_real_check, EP13_imag_check)) + + if num_dq_epics > 13: + A_epics_threshold += 2.0 + EP14_real_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP14_real"], 1.0 / EP14_real) + EP14_imag_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP14_imag"], 1.0 / EP14_imag) + + # The above values should be close to one (within 1 / 10^4) + EP14_real_check = pipeparts.mkgeneric(pipeline, EP14_real_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + EP14_imag_check = pipeparts.mkgeneric(pipeline, EP14_imag_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + A_epics_check = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, A_epics_check, EP14_real_check, EP14_imag_check)) + + A_epics_bit = pipeparts.mkbitvectorgen(pipeline, A_epics_check, bit_vector = pow(2,22), threshold = A_epics_threshold) + A_epics_bit = pipeparts.mkgeneric(pipeline, A_epics_bit, "lal_logicalundersample", required_on = pow(2,22), status_out = pow(2,22)) + A_epics_bit = pipeparts.mkcapsfilter(pipeline, A_epics_bit, calibstate_caps) + + # Next, check the EPICS that involve only the sensing function C, EP6 and EP11 + if num_dq_epics > 5: + C_epics_threshold += 1.5 + EP6_real_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP6_real"], 1.0 / EP6_real) + EP6_imag_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP6_imag"], 1.0 / EP6_imag) + + # The above values should be close to one (within 1 / 10^4) + EP6_real_check = pipeparts.mkgeneric(pipeline, EP6_real_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + EP6_imag_check = pipeparts.mkgeneric(pipeline, EP6_imag_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + C_epics_check = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, EP6_real_check, EP6_imag_check)) + + if num_dq_epics > 10: + C_epics_threshold += 2.0 + EP11_real_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP11_real"], 1.0 / EP11_real) + EP11_imag_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP11_imag"], 1.0 / EP11_imag) + + # The above values should be close to one (within 1 / 10^4) + EP11_real_check = pipeparts.mkgeneric(pipeline, EP11_real_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + EP11_imag_check = pipeparts.mkgeneric(pipeline, EP11_imag_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + C_epics_check = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, C_epics_check, EP11_real_check, EP11_imag_check)) + + C_epics_bit = pipeparts.mkbitvectorgen(pipeline, C_epics_check, bit_vector = pow(2,23), threshold = C_epics_threshold) + C_epics_bit = pipeparts.mkgeneric(pipeline, C_epics_bit, "lal_logicalundersample", required_on = pow(2,23), status_out = pow(2,23)) + C_epics_bit = pipeparts.mkcapsfilter(pipeline, C_epics_bit, calibstate_caps) + + # Next, check the remaining EPICS that are combinations of D, A, and C, EP1 and EP2 + if num_dq_epics > 0: + other_epics_threshold += 1.5 + EP1_real_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP1_real"], 1.0 / EP1_real) + EP1_imag_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP1_imag"], 1.0 / EP1_imag) + + # The above values should be close to one (within 1 / 10^4) + EP1_real_check = pipeparts.mkgeneric(pipeline, EP1_real_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + EP1_imag_check = pipeparts.mkgeneric(pipeline, EP1_imag_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + other_epics_check = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, EP1_real_check, EP1_imag_check)) + + if num_dq_epics > 1: + other_epics_threshold += 2.0 + EP2_real_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP2_real"], 1.0 / EP2_real) + EP2_imag_check = pipeparts.mkaudioamplify(pipeline, head_dict["EP2_imag"], 1.0 / EP2_imag) + + # The above values should be close to one (within 1 / 10^4) + EP2_real_check = pipeparts.mkgeneric(pipeline, EP2_real_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + EP2_imag_check = pipeparts.mkgeneric(pipeline, EP2_imag_check, "lal_insertgap", bad_data_intervals = [0.9999, 1.0001], replace_value = 0.0, insert_gap = False) + other_epics_check = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, other_epics_check, EP2_real_check, EP2_imag_check)) + + other_epics_bit = pipeparts.mkbitvectorgen(pipeline, other_epics_check, bit_vector = pow(2,24), threshold = other_epics_threshold) + other_epics_bit = pipeparts.mkgeneric(pipeline, other_epics_bit, "lal_logicalundersample", required_on = pow(2,24), status_out = pow(2,24)) + epics_bits = pipeparts.mkcapsfilter(pipeline, other_epics_bit, calibstate_caps) + + # Add the EPICS bits together + if num_dq_epics > 6: + # There are EPICS for D, A, and C + epics_bits = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, epics_bits, A_epics_bit, C_epics_bit, D_epics_bit)) + elif num_dq_epics > 5: + # There are EPICS for A and C + epics_bits = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, epics_bits, A_epics_bit, C_epics_bit)) + elif num_dq_epics > 2: + # There are EPICS for A + epics_bits = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, epics_bits, A_epics_bit)) + + # + # COMBINE ALL BITS TO MAKE GDS-KAPPA_STATE_VECTOR + # + + # We don't know what is going to be in the state vector, but we know there will be something, so make a list and send that to an adder + dq_list = [] + if not options.no_kappatst: + dq_list.append(ktstSmoothInRange) + if not options.no_kappap or not options.no_kappapu: + dq_list.append(kpSmoothInRange) + if not options.no_kappau: + dq_list.append(kuSmoothInRange) + if not options.no_kappac: + dq_list.append(kcSmoothInRange) + if not options.no_fcc: + dq_list.append(fccSmoothInRange) + if not options.no_fs: + dq_list.append(fsSmoothInRange) + if not options.no_srcQ: + dq_list.append(srcQSmoothInRange) + if not options.no_coherence: + dq_list.append(coherence_bits) + if num_dq_epics > 0: + dq_list.append(epics_bits) + + kappastatevector = calibration_parts.mkadder(pipeline, tuple(dq_list)) + + kappastatevector = pipeparts.mkprogressreport(pipeline, kappastatevector, "progress_calibstatevec_%s" % instrument) + dqtagstr = "channel-name=%s:GDS-KAPPA_STATE_VECTOR, instrument=%s" % (instrument, instrument) + kappastatevector = pipeparts.mktaginject(pipeline, kappastatevector, dqtagstr) + + +# +# Produce time-dependent correction factors to be recorded in the frames +# + +record_kappa_caps = "audio/x-raw, format=F32LE, rate=%d" % options.record_factors_sr + +# Resample the \kappa_pu channels at the specified recording sample rate and change them to single precision channels +if not options.no_kappapu: + + kpuRout = pipeparts.mkaudioconvert(pipeline, smooth_kpuRtee) + kpuRout = calibration_parts.mkresample(pipeline, kpuRout, 1, False, record_kappa_caps) + kpuRout = pipeparts.mkprogressreport(pipeline, kpuRout, "progress_kappa_pu_real_%s" % instrument) + + kpuIout = pipeparts.mkaudioconvert(pipeline, smooth_kpuItee) + kpuIout = calibration_parts.mkresample(pipeline, kpuIout, 1, False, record_kappa_caps) + kpuIout = pipeparts.mkprogressreport(pipeline, kpuIout, "progress_kappa_pu_imag_%s" % instrument) + + smooth_kpuR_nogate = pipeparts.mkaudioconvert(pipeline, smooth_kpuR_nogate) + smooth_kpuR_nogate = calibration_parts.mkresample(pipeline, smooth_kpuR_nogate, 1, False, record_kappa_caps) + smooth_kpuR_nogate = pipeparts.mkprogressreport(pipeline, smooth_kpuR_nogate, "progress_kappa_pu_real_nogate_%s" % instrument) + + smooth_kpuI_nogate = pipeparts.mkaudioconvert(pipeline, smooth_kpuI_nogate) + smooth_kpuI_nogate = calibration_parts.mkresample(pipeline, smooth_kpuI_nogate, 1, False, record_kappa_caps) + smooth_kpuI_nogate = pipeparts.mkprogressreport(pipeline, smooth_kpuI_nogate, "progress_kappa_pu_imag_nogate_%s" % instrument) + +# Resample the \kappa_tst channels at the specified recording sample rate and change them to single precision channels +if not options.no_kappatst: + + ktstRout = pipeparts.mkaudioconvert(pipeline, smooth_ktstRtee) + ktstRout = calibration_parts.mkresample(pipeline, ktstRout, 1, False, record_kappa_caps) + ktstRout = pipeparts.mkprogressreport(pipeline, ktstRout, "progress_kappa_tst_real_%s" % instrument) + + ktstIout = pipeparts.mkaudioconvert(pipeline, smooth_ktstItee) + ktstIout = calibration_parts.mkresample(pipeline, ktstIout, 1, False, record_kappa_caps) + ktstIout = pipeparts.mkprogressreport(pipeline, ktstIout, "progress_kappa_tst_imag_%s" % instrument) + + smooth_ktstR_nogate = pipeparts.mkaudioconvert(pipeline, smooth_ktstR_nogate) + smooth_ktstR_nogate = calibration_parts.mkresample(pipeline, smooth_ktstR_nogate, 1, False, record_kappa_caps) + smooth_ktstR_nogate = pipeparts.mkprogressreport(pipeline, smooth_ktstR_nogate, "progress_kappa_tst_real_nogate_%s" % instrument) + + smooth_ktstI_nogate = pipeparts.mkaudioconvert(pipeline, smooth_ktstI_nogate) + smooth_ktstI_nogate = calibration_parts.mkresample(pipeline, smooth_ktstI_nogate, 1, False, record_kappa_caps) + smooth_ktstI_nogate = pipeparts.mkprogressreport(pipeline, smooth_ktstI_nogate, "progress_kappa_tst_imag_nogate_%s" % instrument) + +# Resample the \kappa_c channels at the specified recording sample rate and change it to a single precision channel +if not options.no_kappac: + kcout = pipeparts.mkaudioconvert(pipeline, smooth_kctee) + kcout = calibration_parts.mkresample(pipeline, kcout, 1, False, record_kappa_caps) + kcout = pipeparts.mkprogressreport(pipeline, kcout, "progress_kappa_c_%s" % instrument) + + smooth_kc_nogate = pipeparts.mkaudioconvert(pipeline, smooth_kc_nogate) + smooth_kc_nogate = calibration_parts.mkresample(pipeline, smooth_kc_nogate, 1, False, record_kappa_caps) + smooth_kc_nogate = pipeparts.mkprogressreport(pipeline, smooth_kc_nogate, "progress_kappa_c_nogate_%s" % instrument) + +# Resample the f_cc channels at the specified recording sample rate and change it to a single precision channel +if not options.no_fcc: + fccout = pipeparts.mkaudioconvert(pipeline, smooth_fcctee) + fccout = calibration_parts.mkresample(pipeline, fccout, 1, False, record_kappa_caps) + fccout = pipeparts.mkprogressreport(pipeline, fccout, "progress_f_cc_%s" % instrument) + + smooth_fcc_nogate = pipeparts.mkaudioconvert(pipeline, smooth_fcc_nogate) + smooth_fcc_nogate = calibration_parts.mkresample(pipeline, smooth_fcc_nogate, 1, False, record_kappa_caps) + smooth_fcc_nogate = pipeparts.mkprogressreport(pipeline, smooth_fcc_nogate, "progress_f_cc_nogate_%s" % instrument) + +# Resample the f_s channels at the specified recording sample rate and change it to a single precision channel +if not options.no_fs: + fsout = pipeparts.mkaudioconvert(pipeline, smooth_fs) + fsout = calibration_parts.mkresample(pipeline, fsout, 1, False, record_kappa_caps) + fsout = pipeparts.mkprogressreport(pipeline, fsout, "progress_f_s_%s" % instrument) + + smooth_fs_nogate = pipeparts.mkaudioconvert(pipeline, smooth_fs_nogate) + smooth_fs_nogate = calibration_parts.mkresample(pipeline, smooth_fs_nogate, 1, False, record_kappa_caps) + smooth_fs_nogate = pipeparts.mkprogressreport(pipeline, smooth_fs_nogate, "progress_f_s_nogate_%s" % instrument) + +# Resample the f_s channels at the specified recording sample rate and change it to a single precision channel +if not options.no_srcQ: + srcQ_inv_out = pipeparts.mkaudioconvert(pipeline, smooth_srcQ_inv) + srcQ_inv_out = calibration_parts.mkresample(pipeline, srcQ_inv_out, 1, False, record_kappa_caps) + srcQ_inv_out = pipeparts.mkprogressreport(pipeline, srcQ_inv_out, "progress_SRC_Q_%s" % instrument) + + smooth_srcQ_inv_nogate = pipeparts.mkaudioconvert(pipeline, smooth_srcQ_inv_nogate) + smooth_srcQ_inv_nogate = calibration_parts.mkresample(pipeline, smooth_srcQ_inv_nogate, 1, False, record_kappa_caps) + smooth_srcQ_inv_nogate = pipeparts.mkprogressreport(pipeline, smooth_srcQ_inv_nogate, "progress_SRC_Q_nogate_%s" % instrument) + +# +# CREATE MUXER AND HOOK EVERYTHING UP TO IT +# + +mux = pipeparts.mkframecppchannelmux(pipeline, None) + +if options.frame_duration is not None: + mux.set_property("frame-duration", options.frame_duration) +if options.frames_per_file is not None: + mux.set_property("frames-per-file", options.frames_per_file) +mux.set_property("compression-scheme", options.compression_scheme) +mux.set_property("compression-level", options.compression_level) + +# Link the output DQ vectors up to the muxer, if applicable +if not options.no_dq_vector: + calibration_parts.mkqueue(pipeline, kappastatevector).get_static_pad("src").link(mux.get_request_pad("%s:%sKAPPA_STATE_VECTOR%s" % (instrument, chan_prefix, chan_suffix))) + +# Link the real and imaginary parts of \kappa_tst to the muxer +if not options.no_kappatst: + calibration_parts.mkqueue(pipeline, ktstRout).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_KAPPA_TST_REAL%s" % (instrument, chan_prefix, chan_suffix))) + calibration_parts.mkqueue(pipeline, ktstIout).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_KAPPA_TST_IMAGINARY%s" % (instrument, chan_prefix, chan_suffix))) + calibration_parts.mkqueue(pipeline, smooth_ktstR_nogate).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_KAPPA_TST_REAL_NOGATE%s" % (instrument, chan_prefix, chan_suffix))) + calibration_parts.mkqueue(pipeline, smooth_ktstI_nogate).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_KAPPA_TST_IMAGINARY_NOGATE%s" % (instrument, chan_prefix, chan_suffix))) + +# Link the real and imaginary parts of \kappa_pu to the muxer +if not options.no_kappapu: + calibration_parts.mkqueue(pipeline, kpuRout).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_KAPPA_PU_REAL%s" % (instrument, chan_prefix, chan_suffix))) + calibration_parts.mkqueue(pipeline, kpuIout).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_KAPPA_PU_IMAGINARY%s" % (instrument, chan_prefix, chan_suffix))) + calibration_parts.mkqueue(pipeline, smooth_kpuR_nogate).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_KAPPA_PU_REAL_NOGATE%s" % (instrument, chan_prefix, chan_suffix))) + calibration_parts.mkqueue(pipeline, smooth_kpuI_nogate).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_KAPPA_PU_IMAGINARY_NOGATE%s" % (instrument, chan_prefix, chan_suffix))) + +# Link the \kappa_c to the muxer +if not options.no_kappac: + calibration_parts.mkqueue(pipeline, kcout).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_KAPPA_C%s" % (instrument, chan_prefix, chan_suffix))) + calibration_parts.mkqueue(pipeline, smooth_kc_nogate).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_KAPPA_C_NOGATE%s" % (instrument, chan_prefix, chan_suffix))) + +# Link the f_cc to the muxer +if not options.no_fcc: + calibration_parts.mkqueue(pipeline, fccout).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_F_CC%s" % (instrument, chan_prefix, chan_suffix))) + calibration_parts.mkqueue(pipeline, smooth_fcc_nogate).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_F_CC_NOGATE%s" % (instrument, chan_prefix, chan_suffix))) + +# Link the f_s to the muxer +if not options.no_fs: + calibration_parts.mkqueue(pipeline, fsout).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_F_S%s" % (instrument, chan_prefix, chan_suffix))) + calibration_parts.mkqueue(pipeline, smooth_fs_nogate).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_F_S_NOGATE%s" % (instrument, chan_prefix, chan_suffix))) + +# Link the src_Q to the muxer +if not options.no_srcQ: + calibration_parts.mkqueue(pipeline, srcQ_inv_out).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_SRC_Q_INVERSE%s" % (instrument, chan_prefix, chan_suffix))) + calibration_parts.mkqueue(pipeline, smooth_srcQ_inv_nogate).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_SRC_Q_INVERSE_NOGATE%s" % (instrument, chan_prefix, chan_suffix))) + +# Check that all frames are long enough, that they have all of the channels by requring a certain amount of time from start-up, and that frames aren't written for times requested by the wings option +def check_complete_frames(pad, info, (output_start, frame_duration, wings_start, wings_end)): + buf = info.get_buffer() + startts = lal.LIGOTimeGPS(0, buf.pts) + duration = lal.LIGOTimeGPS(0, buf.duration) + if not (startts % frame_duration == 0): + return Gst.PadProbeReturn.DROP + if startts < output_start: + return Gst.PadProbeReturn.DROP + if duration != frame_duration: + return Gst.PadProbeReturn.DROP + if wings_start is not None and wings_end is not None: + if startts < wings_start or (startts+duration) > wings_end: + return Gst.PadProbeReturn.DROP + return Gst.PadProbeReturn.OK +if options.data_source == "frames": + start = int(options.gps_start_time) +elif options.data_source == "lvshm": + tm = time.gmtime() + start = int(lal.UTCToGPS(tm)) +# start time of first frame file is the desired start time + kappa settling +output_start = start + options.demodulation_filter_time + options.median_smoothing_time + options.factors_averaging_time + +# If the wings option is set, need to also check that frames aren't written during the requested wing time +if options.wings is not None: + wings_start = int(options.gps_start_time) + options.wings + wings_end = int(options.gps_end_time) - options.wings + mux.get_static_pad("src").add_probe(Gst.PadProbeType.BUFFER, check_complete_frames, (lal.LIGOTimeGPS(output_start,0), lal.LIGOTimeGPS(options.frame_duration*options.frames_per_file,0), lal.LIGOTimeGPS(wings_start, 0), lal.LIGOTimeGPS(wings_end, 0))) +else: + mux.get_static_pad("src").add_probe(Gst.PadProbeType.BUFFER, check_complete_frames, (lal.LIGOTimeGPS(output_start,0), lal.LIGOTimeGPS(options.frame_duration*options.frames_per_file,0), None, None)) + +mux = pipeparts.mkprogressreport(pipeline, mux, "progress_sink_%s" % instrument) + +if options.write_to_shm_partition is not None: + pipeparts.mkgeneric(pipeline, mux, "gds_lvshmsink", sync=False, async=False, shm_name = options.write_to_shm_partition, num_buffers=10, blocksize=options.frame_size*options.frame_duration*options.frames_per_file, buffer_mode=options.buffer_mode) +else: + pipeparts.mkframecppfilesink(pipeline, mux, frame_type = options.frame_type, path = options.output_path, instrument = instrument) + +# Run pipeline + +if options.write_pipeline is not None: + pipeparts.write_dump_dot(pipeline, "%s.%s" %(options.write_pipeline, "NULL"), verbose = options.verbose) + +# Seek the pipeline when necessary +if options.data_source == "frames": + if options.verbose: + print >>sys.stderr, "seeking GPS start and stop times ..." + if pipeline.set_state(Gst.State.READY) != Gst.StateChangeReturn.SUCCESS: + raise RuntimeError("pipeline failed to enter READY state") + datasource.pipeline_seek_for_gps(pipeline, gps_start_time, gps_end_time) + +if options.verbose: + print >>sys.stderr, "setting pipeline state to playing ..." +if pipeline.set_state(Gst.State.PLAYING) != Gst.StateChangeReturn.SUCCESS: + raise RuntimeError("pipeline failed to enter PLAYING state") +else: + print "set to playing successfully" +if options.write_pipeline is not None: + pipeparts.write_dump_dot(pipeline, "%s.%s" %(options.write_pipeline, "PLAYING"), verbose = options.verbose) + +if options.verbose: + print >>sys.stderr, "running pipeline ..." + +mainloop.run() + +if pipeline.set_state(Gst.State.NULL) != Gst.StateChangeReturn.SUCCESS: + raise RuntimeError("pipeline could not be set to NULL") diff --git a/gstlal-calibration/bin/gstlal_compute_strain b/gstlal-calibration/bin/gstlal_compute_strain index f8b259e360..8ed27a55a8 100755 --- a/gstlal-calibration/bin/gstlal_compute_strain +++ b/gstlal-calibration/bin/gstlal_compute_strain @@ -18,27 +18,33 @@ """ -This pipeline produces h(t) given DARM_ERR and DARM_CTRL or given DELTAL_RESIDUAL and DELTAL_CTRL. It can be run online in real-time or offline on frame files. It can write h(t) frames to frame files or to a shared memory partition. +This pipeline produces h(t) given DARM_ERR and DARM_CTRL or given DELTAL_RESIDUAL and DELTAL_CTRL. It can be run online in real-time or offline on frame files. It can write h(t) frames to frame files or to a shared memory partition. -The differential arm length resulting from external sources is +The differential arm length resulting from external sources is -\Delta L_{ext} = d_{err}/(\kappa_c C) + (A_tst * \kappa_tst + A_pu * \kappa_pu) d_{ctrl} +\Delta L_{ext} = ((f^2 + f_s^2 - i * f * f_s / Q) / f^2) +* ((1 + i * f / f_cc) / (\kappa_c C_res)) * d_{err} ++ (A_tst * \kappa_tst + A_pu * \kappa_pu) * d_{ctrl} -where C is the sensing function, A_tst is the TST acutuation function, A_pu is the PUM+UIM actuation, \kappa_c is the time dependent gain of the sensing function, \kappa_tst is the time-dependent gain of TST actuation, and \kappa_pu is the time-dependent gain of the PUM/UIM actuation. \Delta L_{ext} is divided by the average arm length (4000 km) to obtain h(t), the external strain in the detectors, +where C is the static portion of the sensing function, A_tst is the TST actuation function, A_pu is the PUM+UIM actuation, \kappa_c is the time-dependent gain of the sensing function, \kappa_tst is the time-dependent gain of TST actuation, and \kappa_pu is the time-dependent gain of the PUM/UIM actuation. \Delta L_{ext} is divided by the average arm length (4000 km) to obtain h(t), the external strain in the detectors, -h(t) = \Delta L_{ext} / L . +h(t) = \Delta L_{ext} / L . -The time-dependent gains (\kappa's) as well as the value for the coupled cavity pole (f_cc) and SRC detuning parameters are calcuated in this pipeline as well. +The time-dependent gains (\kappa's) as well as the value for the coupled cavity pole f_cc and SRC detuning parameters f_s and Q are calcuated in this pipeline as well. -This pipeline will most often be run in a format where it picks up after part of the actuation and sensing functions have been applied to the apporiate channels. In this mode, the input channels are \Delta L_{res} and \Delta L_{ctrl, i}. This pipeline then applies further high frequency corrections to each of these channels, applies the appropriate time delay to each channel, adds the channels together, and divides by L. +This pipeline will most often be run in a format where it picks up after part of the actuation and sensing functions have been applied to the appropriate channels. In this mode, the input channels are \Delta L_{res} and \Delta L_{ctrl, i}. This pipeline then applies further high frequency corrections to each of these channels, applies the appropriate time delay to each channel, adds the channels together, and divides by L. -h(t) = (\Delta L_{res} * (1 / \kappa_c) * corrections + (\Delta L_{ctrl, TST} * \kappa_tst + (\Delta L_{ctrl, P} + \Delta L_{ctrl, U})* \kappa_pu) * corrections) / L +h(t) = (((f^2 + f_s^2 - i * f * f_s / Q) / f^2) +* ((1 + i * f / f_cc) / \kappa_c) * corrections * \Delta L_{res} ++ \kappa_tst * \Delta L_{ctrl, TST} ++ \kappa_pu * (\Delta L_{ctrl, P} + \Delta L_{ctrl, U})) / L -Note: The \kappa's are complex numbers. Only the real part of the computed \kappa's are applied as time-dependent gain corrections. +Note: The actuation \kappa's are complex numbers. Only the real part of the computed \kappa's are applied as time-dependent gain corrections. Further documentation explaining the time domain calibration procedure can be found in LIGO DCC #T1400256 and #P1700236. For a full list of example command lines that were used to create the O1 h(t) frames, see https://wiki.ligo.org/Calibration/GDSCalibrationConfigurationsO1. +For a full list of example command lines that were used to create the O2 h(t) frames, see https://wiki.ligo.org/Calibration/GDSCalibrationConfigurationsO2. Type gstlal_compute_strain --help to see the full list of command line options. """ @@ -192,7 +198,6 @@ parser.add_option("--apply-kappatst", action = "store_true", help = "Set this to parser.add_option("--apply-complex-kappatst", action = "store_true", help = "Set this to have the \kappa_tst factors filter the actuation chain with an adaptive filter that corrects for both magnitude and phase errors.") parser.add_option("--act-timing-from-kappatst", action = "store_true", help = "Set this to use the calculated value of \kappa_tst to measure any timing error in the actuation. If this is set, the phase of \kappa_pu will be adjusted accordingly.") parser.add_option("--apply-kappac", action = "store_true", help = "Set this to have the \kappa_c factors multiply the sensing chain.") -parser.add_option("--apply-fcc", action = "store_true", help = "Set this to have the f_cc time-dependent, frequency-dependent corrections applied.") parser.add_option("--compute-factors-sr", metavar = "Hz", type = int, default = 16, help = "Sample rate at which calibration factors are computed. (Default = 16 Hz)") parser.add_option("--demodulation-filter-time", metavar = "s", type = int, default = 20, help = "Length in seconds of low-pass FIR filter used in demodulation of the calibration lines. (Default = 20 seconds)") parser.add_option("--median-smoothing-time", metavar = "s", type = int, default = 128, help = "Time (in seconds) to smooth out \kappas with a median-like method. (Default = 128 s)") @@ -290,7 +295,7 @@ parser.add_option("--darm-ctrl-channel-name", metavar = "name", default = "CAL-D parser.add_option("--darm-err-channel-name", metavar = "name", default = "CAL-DARM_ERR_WHITEN_OUT_DBL_DQ", help = "Set the name of the error signal channel. (Default = CAL-DARM_ERR_WHITEN_OUT_DBL_DQ)") # These options are specific to the partial calibration mode -parser.add_option("--partial-calibration", action = "store_true", help = "Set this to run the pipeline in partial calibraiton mode.") +parser.add_option("--partial-calibration", action = "store_true", help = "Set this to run the pipeline in partial calibration mode.") parser.add_option("--deltal-tst-channel-name", metavar = "name", default = "CAL-DELTAL_CTRL_TST_DBL_DQ", help = "Set the name of the partially calibrated control channel for the TST branch of the actuation. (Default = CAL-DELTAL_CTRL_TST_DBL_DQ)") parser.add_option("--deltal-pum-channel-name", metavar = "name", default = "CAL-DELTAL_CTRL_PUM_DBL_DQ", help = "Set the name of the partially calibrated control channel for the PUM/UIM branch of the actuation. (Default = CAL-DELTAL_CTRL_PUM_DBL_DQ)") parser.add_option("--deltal-uim-channel-name", metavar = "name", default = "CAL-DELTAL_CTRL_UIM_DBL_DQ", help = "Set the name of the partially calibrated control channel for the PUM/UIM branch of the actuation. (Default = CAL-DELTAL_CTRL_UIM_DBL_DQ)") @@ -409,7 +414,7 @@ td = not options.frequency_domain_filtering # If we are using EPICS from frames and removing calibration lines, we need EP10 to remove the ESD line. Otherwise, we just remove the other lines if possible. if (not options.factors_from_filters_file) and options.remove_callines and (options.data_source == "lvshm" or (options.ifo == "H1" and int(options.gps_start_time) > 1175976256) or (options.ifo == "L1" and int(options.gps_start_time) > 1179588864)): remove_esd_act_line = True -elif not options.factors_from_filters_file: +else: remove_esd_act_line = False # How many EPICS will we for the CALIB_STATE_VECTOR calculation? It depends on the IFO and the time we are calibrating @@ -452,14 +457,18 @@ if options.factors_from_filters_file or not options.no_dq_vector: except: if options.factors_from_filters_file: raise ValueError("Cannot compute time-dependent correction factors, as the needed EPICS are not contained in the specified filters file.") + if num_dq_epics > 0: + num_dq_epics = 0 try: EP10_real = float(filters["EP10_real"]) EP10_imag = float(filters["EP10_imag"]) - if options.factors_from_filters_file: + if options.factors_from_filters_file and options.remove_callines: remove_esd_act_line = True except: if options.factors_from_filters_file: remove_esd_act_line = False + if num_dq_epics > 9: + num_dq_epics = 9 try: EP11_real = float(filters["EP11_real"]) EP11_imag = float(filters["EP11_imag"]) @@ -472,6 +481,8 @@ if options.factors_from_filters_file or not options.no_dq_vector: except: if options.factors_from_filters_file and (not options.no_srcQ or not options.no_fs): raise ValueError("Cannot compute SRC spring frequency or Q, as the needed EPICS are not contained in the specified filters file.") + if num_dq_epics > 10: + num_dq_epics = 10 # Load all of the kappa dewhitening and correction factors darm_act_line_freq = float(filters["ka_pcal_line_freq"]) @@ -498,7 +509,7 @@ try: high_pcal_line_freq = float(filters["high_pcal_line_freq"]) pcal_corr_at_high_line_freq_real = float(filters["high_pcal_corr_re"]) pcal_corr_at_high_line_freq_imag = float(filters["high_pcal_corr_im"]) - if high_pcal_line_freq > 0: + if high_pcal_line_freq > 0 and options.remove_callines: remove_high_pcal_line = True else: remove_high_pcal_line = False @@ -508,7 +519,7 @@ try: roaming_pcal_line_freq = float(filters["roaming_pcal_line_freq"]) pcal_corr_at_roaming_line_real = float(filters["roaming_pcal_corr_re"]) pcal_corr_at_roaming_line_imag = float(filters["roaming_pcal_corr_im"]) - if roaming_pcal_line_freq > 0.0: + if roaming_pcal_line_freq > 0.0 and options.remove_callines: remove_roaming_pcal_line = True else: remove_roaming_pcal_line = False @@ -618,25 +629,25 @@ if not options.no_dq_vector: # If we are computing the factors in the pipeline, we need the reference model EPICS records # Needed for kappa_tst -if options.factors_from_filters_file and (not options.no_kappatst or not options.no_kappapu or not options.no_kappap or not options.no_kappau or not options.no_kappac or not options.no_fcc or not options.no_fs or not options.no_srcQ) or num_dq_epics > 0: +if not options.factors_from_filters_file and (not options.no_kappatst or not options.no_kappapu or not options.no_kappap or not options.no_kappau or not options.no_kappac or not options.no_fcc or not options.no_fs or not options.no_srcQ) or num_dq_epics > 0: channel_list.extend(((instrument, options.EP1_real), (instrument, options.EP1_imag))) headkeys.extend(("EP1_real", "EP1_imag")) # These are needed for kappa_pu -if options.factors_from_filters_file and (not options.no_kappac or not options.no_fcc or not options.no_kappapu or not options.no_fs or not options.no_srcQ) or num_dq_epics > 3: +if not options.factors_from_filters_file and (not options.no_kappac or not options.no_fcc or not options.no_kappapu or not options.no_fs or not options.no_srcQ) or num_dq_epics > 3: channel_list.extend(((instrument, options.EP2_real), (instrument, options.EP2_imag), (instrument, options.EP3_real), (instrument, options.EP3_imag), (instrument, options.EP4_real), (instrument, options.EP4_imag))) headkeys.extend(("EP2_real", "EP2_imag", "EP3_real", "EP3_imag", "EP4_real", "EP4_imag")) # If we are computing either kappa_c or f_cc, we need some more EPICS records -if options.factors_from_filters_file and (not options.no_kappac or not options.no_fcc or not options.no_fs or not options.no_srcQ) or num_dq_epics > 8: +if not options.factors_from_filters_file and (not options.no_kappac or not options.no_fcc or not options.no_fs or not options.no_srcQ) or num_dq_epics > 8: channel_list.extend(((instrument, options.EP6_real), (instrument, options.EP6_imag), (instrument, options.EP7_real), (instrument, options.EP7_imag), (instrument, options.EP8_real), (instrument, options.EP8_imag), (instrument, options.EP9_real), (instrument, options.EP9_imag))) headkeys.extend(("EP6_real", "EP6_imag", "EP7_real", "EP7_imag", "EP8_real", "EP8_imag", "EP9_real", "EP9_imag")) # EP10 is needed to remove the ESD line -if options.factors_from_filters_file and remove_esd_act_line or num_epics_channels > 9: +if not options.factors_from_filters_file and remove_esd_act_line or num_dq_epics > 9: channel_list.extend(((instrument, options.EP10_real), (instrument, options.EP10_imag))) headkeys.extend(("EP10_real", "EP10_imag")) # These are needed if we compute the optical spring frequency and/or Q-factor of the Signal Recycling Cavity (SRC) -if options.factors_from_filters_file and (not options.no_fs or not options.no_srcQ) or num_dq_epics > 13: +if not options.factors_from_filters_file and (not options.no_fs or not options.no_srcQ) or num_dq_epics > 13: channel_list.extend(((instrument, options.EP11_real), (instrument, options.EP11_imag), (instrument, options.EP12_real), (instrument, options.EP12_imag), (instrument, options.EP13_real), (instrument, options.EP13_imag), (instrument, options.EP14_real), (instrument, options.EP14_imag))) headkeys.extend(("EP11_real", "EP11_imag", "EP12_real", "EP12_imag", "EP13_real", "EP13_imag", "EP14_real", "EP14_imag")) @@ -1834,7 +1845,7 @@ if options.remove_callines: # Make sure we have demodulated pcal at the ~300 Hz pcal line if options.no_kappac and options.no_fcc and options.no_srcQ and options.no_fs: pcal_at_opt_gain_freq = calibration_parts.demodulate(pipeline, pcaltee, opt_gain_fcc_line_freq, td, options.compute_factors_sr, options.demodulation_filter_time, options.filter_latency, prefactor_real = pcal_corr_at_opt_gain_fcc_freq_real, prefactor_imag = pcal_corr_at_opt_gain_fcc_freq_imag) - # Reconstruct a calibrated pcal at only the ~300 Hz pcal line + # Reconstruct a calibrated pcal at only the ~330 Hz pcal line pcaly_line2 = calibration_parts.mkresample(pipeline, pcal_at_opt_gain_freq, 3, False, "audio/x-raw, format=Z128LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) pcaly_line2 = pipeparts.mkgeneric(pipeline, pcaly_line2, "lal_demodulate", line_frequency = -1.0 * opt_gain_fcc_line_freq, prefactor_real = 2.0) remove_pcaly_line2 = pipeparts.mkgeneric(pipeline, pcaly_line2, "creal") @@ -1854,8 +1865,8 @@ if options.remove_callines: esd_act_line = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, tstexc_at_esd_act_freq, EP10)) # Reconstruct a calibrated ESD injection at the ~30 Hz ESD line if options.apply_kappatst: - # Multiply by the magnitude of kappa_tst - esd_act_line = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, esd_act_line, pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, pipeparts.mkgeneric(pipeline, smooth_ktsttee, "cabs"), matrix=[[1.0, 0.0]])))) + # Multiply by kappa_tst + esd_act_line = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, esd_act_line, smooth_ktsttee)) esd_act_line = calibration_parts.mkresample(pipeline, esd_act_line, 3, False, "audio/x-raw, format=Z128LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) esd_act_line_remove = pipeparts.mkgeneric(pipeline, esd_act_line, "lal_demodulate", line_frequency = -1.0 * esd_act_line_freq, prefactor_real = 2.0) esd_act_line_remove = pipeparts.mkgeneric(pipeline, esd_act_line_remove, "creal") diff --git a/gstlal-calibration/tests/lal_demodulate_test.py b/gstlal-calibration/tests/lal_demodulate_test.py index 982411ddd1..a159cdeb70 100755 --- a/gstlal-calibration/tests/lal_demodulate_test.py +++ b/gstlal-calibration/tests/lal_demodulate_test.py @@ -27,6 +27,7 @@ import numpy import sys from gstlal import pipeparts +from gstlal import calibration_parts import test_common @@ -95,6 +96,40 @@ def lal_demodulate_02(pipeline, name): return pipeline +def lal_demodulate_03(pipeline, name): + # + # This test checks sensitivity of the demodulation process used in the calibration pipeline to small changes in line frequency + # + + rate_in = 16384 # Hz + rate_out = 16 # Hz + buffer_length = 1.0 # seconds + test_duration = 1000 # seconds + + # + # build pipeline + # + + # Make fake data with a signal + src = test_common.test_src(pipeline, buffer_length = buffer_length, rate = rate_in, test_duration = test_duration, wave = 0, volume = 1.0, freq = 37.00, width = 64) + + # Demodulate it + head = calibration_parts.demodulate(pipeline, src, 37.10, True, rate_out, 20, 0.0) + + # Smoothing +# head = pipeparts.mkgeneric(pipeline, head, "lal_smoothkappas", array_size = 128 * rate_out, avg_array_size = 10 * rate_out) + + # Measure the amplitude of the result + head = pipeparts.mkgeneric(pipeline, head, "cabs") + head = pipeparts.mkaudioamplify(pipeline, head, 2.0) + pipeparts.mknxydumpsink(pipeline, head, "%s_out.dump" % name) + + # + # done + # + + return pipeline + # # ============================================================================= # @@ -104,5 +139,7 @@ def lal_demodulate_02(pipeline, name): # -test_common.build_and_run(lal_demodulate_01, "lal_demodulate_01") +#test_common.build_and_run(lal_demodulate_01, "lal_demodulate_01") #test_common.build_and_run(lal_demodulate_02, "lal_demodulate_02") +test_common.build_and_run(lal_demodulate_03, "lal_demodulate_03") + -- GitLab