From dfb49e80b2b1e10c45fddc6d0f4fba618720ae30 Mon Sep 17 00:00:00 2001 From: Madeline Wade <madeline.wade@ligo.org> Date: Thu, 14 Sep 2017 09:52:31 -0700 Subject: [PATCH] Splitting off cleaning tasks into a new program called gstlal_clean_strain --- gstlal-calibration/bin/gstlal_clean_strain | 1907 ++++++++++++++++++ gstlal-calibration/bin/gstlal_compute_strain | 345 +--- 2 files changed, 1925 insertions(+), 327 deletions(-) create mode 100644 gstlal-calibration/bin/gstlal_clean_strain diff --git a/gstlal-calibration/bin/gstlal_clean_strain b/gstlal-calibration/bin/gstlal_clean_strain new file mode 100644 index 0000000000..657faccdbf --- /dev/null +++ b/gstlal-calibration/bin/gstlal_clean_strain @@ -0,0 +1,1907 @@ +#!/usr/bin/env python +# +# Copyright (C) 2010-2015 Jordi Burguet-Castell, Madeline Wade +# +# 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 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 + +\Delta L_{ext} = d_{err}/(\kappa_c C) + (A_tst * \kappa_tst + A_usum * \kappa_pu) d_{ctrl} + +where C is the sensing function, A_tst is the TST acutuation function, A_usum is the PUM+UIM+TOP 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 . + +The time-dependent gains (\kappa's) as well as the value for the coupled cavity pole (f_cc), the time-dependent gain of the PUM actuation (\kappa_pu) and the overall time-depenent gain of the actuation (\kappa_a) 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}. 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, USUM} * \kappa_pu) * corrections) / L + +Note: The \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. + +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. + +Type gstlal_compute_strain --help to see the full list of command line options. +""" + +import sys +import os +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 + +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) + +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("--dq-channel-name", metavar = "name", default = "ODC-MASTER_CHANNEL_OUT_DQ", help = "Set the name of the data quality (or state vector) channel. (Default=ODC-MASTER_CHANNEL_OUT_DQ)") +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("--frame-segments-file", metavar = "filename", help = "Set the name of the LIGO light-weight XML file from which to load frame segments. This is required iff --data-source=frames") +parser.add_option("--frame-segments-name", metavar = "name", help = "Set the name of the segments to extract from the segment tables. This is required iff --frame-segments-file is given") +parser.add_option("--hoft-sample-rate", metavar = "Hz", default = 16384, type = "int", help = "Sample rate at which to generate strain data. This should be less than or equal to the sample rate of the error and control signal channels. (Default = 16384 Hz)") +parser.add_option("--control-sample-rate", metavar = "Hz", default = 16384, type = "int", help = "Sample rate of the control signal channels. (Default = 16384 Hz)") +parser.add_option("--odc-sample-rate", metavar = "Hz", default = 16384, type = "int", help = "Sample rate of the ODC state vector channel. (Default = 16384 Hz)") +parser.add_option("--calib-state-sample-rate", metavar = "Hz", default = 16, type = "int", help = "Sample rate for the outgoing DQ vector GDS-CALIB_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("--frequency-domain-filtering", action = "store_true", help = "Set this to perform filtering routines in the frequency domain instead of using direct convolution.") +parser.add_option("--obs-ready-bitmask", metavar = "bitmask", type = "int", default = 4, help = "Bitmask used on ODC state vector in order to determine OBSERVATION_READY bit information. (Default=4)") +parser.add_option("--obs-intent-bitmask", metavar = "bitmask", type = "int", default = 2, help = "Bitmask used on ODC state vector in order to determine OBSERVATION_INTENT bit information. (Default=2)") +parser.add_option("--hw-inj-cbc-bitmask", metavar = "bitmask", type = "int", default = 16777216, help = "Bitmask used on ODC state vector in order presence of CBC hardware injection. (Default=16777216)") +parser.add_option("--hw-inj-burst-bitmask", metavar = "bitmask", type = "int", default = 33554432, help = "Bitmask used on ODC state vector in order presence of burst hardware injection. (Default=33554432)") +parser.add_option("--hw-inj-detchar-bitmask", metavar = "bitmask", type = "int", default = 67108864, help = "Bitmask used on ODC state vector in order presence of DetChar hardware injection. (Default=67108864)") +parser.add_option("--hw-inj-stoch-bitmask", metavar = "bitmask", type = "int", default = 8388608, help = "Bitmask used on ODC state vector in order presence of stochastic hardware injection. (Default=8388608)") +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 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("--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-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("--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("--apply-kappapu", action = "store_true", help = "Set this to have the \kappa_pu factors multiply the actuation chain.") +parser.add_option("--apply-kappatst", action = "store_true", help = "Set this to have the \kappa_tst factors multiply the actuation chain.") +parser.add_option("--apply-kappac", action = "store_true", help = "Set this to have the \kappa_c factors multiply the sensing chain.") +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)") +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("--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-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-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 = 330.0, help = "Expected value for the coupled cavity pole. (Default = 330.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("--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("--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 = 10, help = "Values of SRC spring frequency +/- this number (in Hz) will be considered OK. (Default = 10 Hz)") +parser.add_option("--srcQ-ok-var", metavar = "float", type = float, default = 20.0, help = "Values of SRC Q +/- this number will be considered OK. (Default = 20)") +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("--dewhitening", action = "store_true", help = "Dewhitening should be used on the relevant channels, since the incoming channels are whitened and single precision.") +parser.add_option("--low-latency", action = "store_true", help = "Run the pipeline in low-latency mode. This uses minimal queueing. Otherwise, maximal queueing is used to prevent the pipeline from locking up.") +parser.add_option("--remove-callines", action = "store_true", help = "Remove calibration lines at known freqencies from h(t) using software.") +parser.add_option("--remove-powerlines", action = "store_true", help = "Remove 60 Hz spectral lines and some harmonics caused by power lines using witness channel PEM-EY_MAINSMON_EBAY_1_DQ.") +parser.add_option("--powerlines-channel-name", metavar = "name", default = "PEM-EY_MAINSMON_EBAY_1_DQ", help = "Set the name of the channel used as input for 60 Hz power lines to be removed. (Default = PEM-EY_MAINSMON_EBAY_1_DQ)") +parser.add_option("--remove-jitter-imc", action = "store_true", help = "Remove laser beam jitter using 4 IMC channels. This can significantly reduce noise in the spectrum.") +parser.add_option("--imc-a-pitch-channel-name", metavar = "name", default = "IMC-WFS_A_DC_PIT_OUT_DQ", help = "Set the name of one of the channels used as input from the IMC for removal of beam jitter noise. (Default = IMC-WFS_A_DC_PIT_OUT_DQ)") +parser.add_option("--imc-b-pitch-channel-name", metavar = "name", default = "IMC-WFS_B_DC_PIT_OUT_DQ", help = "Set the name of one of the channels used as input from the IMC for removal of beam jitter noise. (Default = IMC-WFS_B_DC_PIT_OUT_DQ)") +parser.add_option("--imc-a-yaw-channel-name", metavar = "name", default = "IMC-WFS_A_DC_YAW_OUT_DQ", help = "Set the name of one of the channels used as input from the IMC for removal of beam jitter noise. (Default = IMC-WFS_A_DC_YAW_OUT_DQ)") +parser.add_option("--imc-b-yaw-channel-name", metavar = "name", default = "IMC-WFS_B_DC_YAW_OUT_DQ", help = "Set the name of one of the channels used as input from the IMC for removal of beam jitter noise. (Default = IMC-WFS_B_DC_YAW_OUT_DQ)") +parser.add_option("--remove-jitter-psl", action = "store_true", help = "Remove laser beam jitter using the bullseye photodiode with 3 PSL channels. This can significantly reduce noise in the spectrum.") +parser.add_option("--bullseye-width-channel-name", metavar = "name", default = "PSL-DIAG_BULLSEYE_WID_OUT_DQ", help = "Set the name of one of the channels used as input from the bullseye photodiode for removal of beam jitter noise. (Default = PSL-DIAG_BULLSEYE_WID_OUT_DQ)") +parser.add_option("--bullseye-pitch-channel-name", metavar = "name", default = "PSL-DIAG_BULLSEYE_PIT_OUT_DQ", help = "Set the name of one of the channels used as input from the bullseye photodiode for removal of beam jitter noise. (Default = PSL-DIAG_BULLSEYE_PIT_OUT_DQ)") +parser.add_option("--bullseye-yaw-channel-name", metavar = "name", default = "PSL-DIAG_BULLSEYE_YAW_OUT_DQ", help = "Set the name of one of the channels used as input from the bullseye photodiode for removal of beam jitter noise. (Default = PSL-DIAG_BULLSEYE_YAW_OUT_DQ)") +parser.add_option("--remove-angular-control", action = "store_true", help = "Remove noise caused by angular control. Uses 4 ASC channels.") +parser.add_option("--asc-dhard-pitch-channel-name", metavar = "name", default = "ASC-DHARD_P_OUT_DQ", help = "Set the name of one of the channels used as input from the ASC to remove angular control noise. (Default = ASC-DHARD_P_OUT_DQ)") +parser.add_option("--asc-dhard-yaw-channel-name", metavar = "name", default = "ASC-DHARD_Y_OUT_DQ", help = "Set the name of one of the channels used as input from the ASC to remove angular control noise. (Default = ASC-DHARD_Y_OUT_DQ)") +parser.add_option("--asc-chard-pitch-channel-name", metavar = "name", default = "ASC-CHARD_P_OUT_DQ", help = "Set the name of one of the channels used as input from the ASC to remove angular control noise. (Default = ASC-CHARD_P_OUT_DQ)") +parser.add_option("--asc-chard-yaw-channel-name", metavar = "name", default = "ASC-CHARD_Y_OUT_DQ", help = "Set the name of one of the channels used as input from the ASC to remove angular control noise. (Default = ASC-CHARD_Y_OUT_DQ)") +parser.add_option("--remove-length-control", action = "store_true", help = "Remove noise caused by length control. Uses 3 LSC channels.") +parser.add_option("--lsc-srcl-channel-name", metavar = "name", default = "LSC-SRCL_IN1_DQ", help = "Set the name of one of the channels used as input from the LSC to remove length control noise. (Default = LSC-SRCL_IN1_DQ)") +parser.add_option("--lsc-mich-channel-name", metavar = "name", default = "LSC-MICH_IN1_DQ", help = "Set the name of one of the channels used as input from the LSC to remove length control noise. (Default = LSC-MICH_IN1_DQ)") +parser.add_option("--lsc-prcl-channel-name", metavar = "name", default = "LSC-PRCL_IN1_DQ", help = "Set the name of one of the channels used as input from the LSC to remove length control noise. (Default = LSC-PRCL_IN1_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") + +# These options are specific to the full calibration mode +parser.add_option("--full-calibration", action = "store_true", help = "Set this to run the pipeline in full calibration mode.") +parser.add_option("--darm-ctrl-channel-name", metavar = "name", default = "CAL-DARM_CTRL_WHITEN_OUT_DBL_DQ", help = "Set the name for the control signal channel. (Default = CAL-DARM_CTRL_WHTIEN_OUT_DBL_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 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("--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)") +parser.add_option("--deltal-res-channel-name", metavar = "name", default = "CAL-DELTAL_RESIDUAL_DBL_DQ", help = "Set the name of the partially calibrated residual channel. (Default = CAL-DELTAL_RESIDUAL_DBL_DQ).") + +# 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.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.frame_segments_file is not None and options.data_source != "frames": + raise ValueError("can only give --frame-segments-file if --data-source=frames") + +if options.frame_segments_name is not None and options.frame_segments_file is None: + raise ValueError("can only specify --frame-segments-name if --frame-segments-file is given") + +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 options.full_calibration is None and options.partial_calibration is None or (options.full_calibration is not None and options.partial_calibration is not None): + raise ValueError("must specify one (and only one) mode of the pipeline: either --full-calibration or --partial-calibration") + +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" or options.data_source == "white": + raise ValueError("cannot set --gps-start-time or --gps-end-time with --data-source=lvshm or --data-source=white") + 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 + +# Make segment list if a frame segments file is provided, other set frame_segments to None +if options.frame_segments_file is not None: + # Frame segments from a user defined file + frame_segments = ligolw_segments.segmenttable_get_by_name(utils.load_filename(options.frame_segments_file, contenthandler = datasource.ContentHandler), options.frame_segments_name).coalesce() + if seg is not None: + # clip frame segments to seek segment if it exists (not required, just saves some meory and I/O overhead) + frame_segments = segments.segmentlistdict((instrument, seglist & segments.segmentlist([seg])) for instrument, seglist in frame_segments.items()) +else: + frame_segments = None + +# Set up short-cut names for each of the sample rates used throughout the pipeline and establish caps string shortcuts +hoftsr = options.hoft_sample_rate # Sample rate for h(t) +calibstatesr = options.calib_state_sample_rate # Sample rate for the CALIB_STATE_VECTOR +odcsr = options.odc_sample_rate # Sample rate of the ODC channel that is read in +ctrlsr = options.control_sample_rate # Sample rate of the control channel (such as DARM_CTRL or DELTAL_CTRL) +cohsr = options.coh_sample_rate # Sample rate for the coherence uncertainty channels +hoft_caps = "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr +ctrl_caps = "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % ctrlsr +calibstate_caps = "audio/x-raw, format=U32LE, rate=%d, channel-mask=(bitmask)0x0" % calibstatesr +odc_caps = "audio/x-raw, format=U32LE, rate=%d, channel-mask=(bitmask)0x0" % odcsr +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 +median_smoothing_samples = int(options.median_smoothing_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 + +# If td is true we will perform filtering in the time domain (direct convolution) in all FIR filtering routines below +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.ifo == "H1" and options.data_source == "frames" and int(options.gps_start_time) > 1175954418) or (options.ifo == "H1" and options.data_source == "lvshm" and now() > 1175954418) or (options.ifo == "L1" and options.data_source == "frames" and int(options.gps_start_time) > 1180184418) or (options.ifo == "L1" and options.data_source == "lvshm" and now() > 1180184418)): + remove_esd_act_line = True +elif not options.factors_from_filters_file: + remove_esd_act_line = False + +# +# Load in the filters file that contains filter coefficients, etc. +# + +filters = numpy.load(options.filters_file) + +# If we're reading the reference model factors from the filters file, load them +if options.factors_from_filters_file: + 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"]) + try: + EP10_real = float(filters["EP10_real"]) + EP10_imag = float(filters["EP10_imag"]) + remove_esd_act_line = True + except: + remove_esd_act_line = False + 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 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.") + +# 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"]) + if src_pcal_line_freq > 10.0: + remove_src_pcal_line = True + else: + remove_src_pcal_line = False +except: + remove_src_pcal_line = False + 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: + 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: + remove_high_pcal_line = True + else: + remove_high_pcal_line = False +except: + remove_high_pcal_line = False +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: + remove_roaming_pcal_line = True + else: + remove_roaming_pcal_line = False +except: + remove_roaming_pcal_line = False +try: + fcc_default = float(filters["fcc"]) +except: + fcc_default = options.expected_fcc +if options.dewhitening: + try: + derr_dewhiten_at_darm_act_freq_real = float(filters["ka_pcal_whitener_re"]) + derr_dewhiten_at_darm_act_freq_imag = float(filters["ka_pcal_whitener_im"]) + derr_dewhiten_at_pu_act_freq_real = float(filters["ka_esd_whitener_re"]) + derr_dewhiten_at_pu_act_freq_imag = float(filters["ka_esd_whitener_im"]) + derr_dewhiten_at_opt_gain_fcc_freq_real = float(filters["kc_pcal_whitener_re"]) + derr_dewhiten_at_opt_gain_fcc_freq_imag = float(filters["kc_pcal_whitener_im"]) + derr_dewhiten_at_esd_act_freq_real = float(filters["ktst_esd_whitener_re"]) + derr_dewhiten_at_esd_act_freq_imag = float(filters["ktst_esd_whitener_im"]) + except: + derr_dewhiten_at_darm_act_freq_real = 1.0 + derr_dewhiten_at_darm_act_freq_imag = 0.0 + derr_dewhiten_at_pu_act_freq_real = 1.0 + derr_dewhiten_at_pu_act_freq_imag = 0.0 + derr_dewhiten_at_opt_gain_fcc_freq_real = 1.0 + derr_dewhiten_at_opt_gain_fcc_freq_imag = 0.0 + derr_dewhiten_at_esd_act_freq_real = 1.0 + derr_dewhiten_at_esd_act_freq_imag = 0.0 + +# If we're performing partial calibration, load the deltal filters +if options.partial_calibration: + reschaindelay = int(filters["res_corr_delay"]) + reschainfilt = filters["res_corr_filter"] + tstdelay = pumuimdelay = int(filters["ctrl_corr_delay"]) + tstfilt = pumuimfilt = filters["ctrl_corr_filter"] + tstchainsr = pumuimchainsr = int(filters["ctrl_corr_sr"]) + if options.dewhitening: + tstdewhitensr = int(filters["deltal_tst_dewhiten_sr"]) + pumuimdewhitensr = int(filters["deltal_pumuim_dewhiten_sr"]) + tstdewhitendelay = int(filters["deltal_tst_dewhiten_delay"]) + pumuimdewhitendelay = int(filters["deltal_pumuim_dewhiten_delay"]) + tstdewhiten = filters["deltal_tst_dewhiten"] + pumuimdewhiten = filters["deltal_pumuim_dewhiten"] + resdewhitendelay = int(filters["deltal_res_dewhiten_delay"]) + resdewhiten = filters["deltal_res_dewhiten"] + +# If we're performing full calibration, load the actuation, sensing filters +if options.full_calibration: + tstchainsr = int(filters["actuation_tst_sr"]) + pumuimchainsr = int(filters["actuation_pumuim_sr"]) + tstdelay = int(filters["actuation_tst_delay"]) + pumuimdelay = int(filters["actuation_pumuim_delay"]) + tstfilt = filters["actuation_tst"] + pumuimfilt = filters["actuation_pumuim"] + reschaindelay = int(filters["inv_sens_delay"]) + reschainfilt = filters["inv_sensing"] + if options.dewhitening: + ctrldewhitendelay = int(filters["dewhiten_ctrl_delay"]) + ctrldewhiten = filters["dewhiten_ctrl"] + ctrldewhitensr = int(filters["dewhiten_ctrl_sr"]) + resdewhitendelay = int(filters["dewhiten_err_delay"]) + resdewhiten = filters["dewhiten_err"] + +# If we're removing 60 Hz lines from the spectrum, load another filter +if options.remove_powerlines: + try: + powerlinessr = int(filters["powerlines_sr"]) + powerlinesdelay = int(filters["powerlines_delay"]) + powerlinesfilt = filters["powerlines_filt"] + except: + raise ValueError("Cannot remove 60 Hz lines because the filters file does contain the needed information") + +# If we're removing laser beam jitter noise from the spectrum, load several more filters +if options.remove_jitter_imc: + try: + imcapitsr = int(filters["jitter_imc_a_pit_sr"]) + imcayawsr = int(filters["jitter_imc_a_yaw_sr"]) + imcbpitsr = int(filters["jitter_imc_b_pit_sr"]) + imcbyawsr = int(filters["jitter_imc_b_yaw_sr"]) + imcapitdelay = int(filters["jitter_imc_a_pit_delay"]) + imcayawdelay = int(filters["jitter_imc_a_yaw_delay"]) + imcbpitdelay = int(filters["jitter_imc_b_pit_delay"]) + imcbyawdelay = int(filters["jitter_imc_b_yaw_delay"]) + imcapitfilt = filters["jitter_imc_a_pit_filt"] + imcayawfilt = filters["jitter_imc_a_yaw_filt"] + imcbpitfilt = filters["jitter_imc_b_pit_filt"] + imcbyawfilt = filters["jitter_imc_b_yaw_filt"] + except: + raise ValueError("Cannot remove beam jitter using imc inputs because the filters file does contain the needed information") +if options.remove_jitter_psl: + try: + bullseyewidsr = int(filters["jitter_bullseye_wid_sr"]) + bullseyepitsr = int(filters["jitter_bullseye_pit_sr"]) + bullseyeyawsr = int(filters["jitter_bullseye_yaw_sr"]) + bullseyewiddelay = int(filters["jitter_bullseye_wid_delay"]) + bullseyepitdelay = int(filters["jitter_bullseye_pit_delay"]) + bullseyeyawdelay = int(filters["jitter_bullseye_yaw_delay"]) + bullseyewidfilt = filters["jitter_bullseye_wid_filt"] + bullseyepitfilt = filters["jitter_bullseye_pit_filt"] + bullseyeyawfilt = filters["jitter_bullseye_yaw_filt"] + except: + raise ValueError("Cannot remove beam jitter using bullseye inputs because the filters file does contain the needed information") +if options.remove_angular_control: + try: + ascdpitsr = int(filters["asc_d_pit_sr"]) + ascdyawsr = int(filters["asc_d_yaw_sr"]) + asccpitsr = int(filters["asc_c_pit_sr"]) + asccyawsr = int(filters["asc_c_yaw_sr"]) + ascdpitdelay = int(filters["asc_d_pit_delay"]) + asccyawdelay = int(filters["asc_d_yaw_delay"]) + asccpitdelay = int(filters["asc_c_pit_delay"]) + ascdyawdelay = int(filters["asc_c_yaw_delay"]) + ascdpitfilt = filters["asc_d_pit_filt"] + ascdyawfilt = filters["asc_d_yaw_filt"] + asccpitfilt = filters["asc_c_pit_filt"] + asccyawfilt = filters["asc_c_yaw_filt"] + except: + raise ValueError("Cannot remove angular control noise using ASC inputs because the filters file does contain the needed information") +if options.remove_length_control: + try: + lscsrclsr = int(filters["lsc_srcl_sr"]) + lscmichsr = int(filters["lsc_mich_sr"]) + lscprclsr = int(filters["lsc_prcl_sr"]) + lscsrcldelay = int(filters["lsc_srcl_delay"]) + lscmichdelay = int(filters["lsc_mich_delay"]) + lscprcldelay = int(filters["lsc_prcl_delay"]) + lscsrclfilt = filters["lsc_srcl_filt"] + lscmichfilt = filters["lsc_mich_filt"] + lscprclfilt = filters["lsc_prcl_filt"] + except: + raise ValueError("Cannot remove length control noise using LSC inputs because the filters file does contain the needed information") + + +# Set up queue parameters +if options.low_latency: + queue_factor = 1 +else: + queue_factor = 0 +long_queue = queue_factor * max(float(len(reschainfilt) - 1) / hoftsr, float(len(tstfilt) - 1) / tstchainsr, float(len(pumuimfilt) - 1) / pumuimchainsr) +short_queue = -1.0 * queue_factor + +time_kc_before_res = (len(reschainfilt) - 1) / int(hoftsr * options.buffer_length) * options.buffer_length +time_res_before_kc = float(options.buffer_length * hoftsr - len(reschainfilt) + 1) / hoftsr +time_kc_ahead_of_res = float(reschaindelay) / hoftsr +time_res_ahead_of_kc = options.buffer_length - time_kc_ahead_of_res +kc_queue_length = queue_factor * max(time_kc_before_res, time_kc_ahead_of_res) if max(time_kc_before_res, time_kc_ahead_of_res) != 0 else -1 * queue_factor +reskc_queue_length = queue_factor * max(time_res_before_kc, time_res_ahead_of_kc) if max(time_res_before_kc, time_res_ahead_of_kc) != 0 else -1 * queue_factor + +time_ktst_before_tst = (len(tstfilt) - 1) / int(tstchainsr * options.buffer_length) * options.buffer_length +time_tst_before_ktst = float(options.buffer_length * tstchainsr - len(tstfilt) + 1) / tstchainsr +time_ktst_ahead_of_tst = float(tstdelay) / tstchainsr +time_tst_ahead_of_ktst = options.buffer_length - time_ktst_ahead_of_tst +ktst_queue_length = queue_factor * max(time_ktst_before_tst, time_ktst_ahead_of_tst) if max(time_ktst_before_tst, time_ktst_ahead_of_tst) != 0 else -1 * queue_factor +tst_queue_length = queue_factor * max(time_tst_before_ktst, time_tst_ahead_of_ktst) if max(time_tst_before_ktst, time_tst_ahead_of_ktst) != 0 else -1 * queue_factor + +time_kpu_before_pu = (len(pumuimfilt) - 1) / int(pumuimchainsr * options.buffer_length) * options.buffer_length +time_pu_before_kpu = float(options.buffer_length * pumuimchainsr - len(pumuimfilt) + 1) / pumuimchainsr +time_kpu_ahead_of_pu = float(pumuimdelay) / pumuimchainsr +time_pu_ahead_of_kpu = options.buffer_length - time_kpu_ahead_of_pu +kpu_queue_length = queue_factor * max(time_kpu_before_pu, time_kpu_ahead_of_pu) if max(time_kpu_before_pu, time_kpu_ahead_of_pu) != 0 else -1 * queue_factor +pumuim_queue_length = queue_factor * max(time_pu_before_kpu, time_pu_ahead_of_kpu) if max(time_pu_before_kpu, time_pu_ahead_of_kpu) != 0 else -1 * queue_factor + +time_res_before_ctrl = float(((len(tstfilt) - 1) / int(tstchainsr * options.buffer_length)) * hoftsr * options.buffer_length - len(reschainfilt) + 1) / hoftsr +time_ctrl_before_res = float(((len(reschainfilt) - 1) / int(hoftsr * options.buffer_length) + 1) * tstchainsr * options.buffer_length - len(tstfilt) + 1) / tstchainsr +time_res_ahead_of_ctrl = float(tstdelay) / tstchainsr - float(reschaindelay) / hoftsr +time_ctrl_ahead_of_res = options.buffer_length - time_res_ahead_of_ctrl +res_queue_length = queue_factor * max(time_res_before_ctrl, time_res_ahead_of_ctrl) if max(time_res_before_ctrl, time_res_ahead_of_ctrl) != 0 else -1 * queue_factor +ctrl_queue_length = queue_factor * max(time_ctrl_before_res, time_ctrl_ahead_of_res) if max(time_ctrl_before_res, time_ctrl_ahead_of_res) != 0 else -1 * queue_factor + +# +# 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 = [] + +# If we are computing the CALIB_STATE_VECTOR, we need the ODC state vector +if not options.no_dq_vector: + channel_list.append((instrument, options.dq_channel_name)) + headkeys.append("odcstatevector") + +# If we are computing the factors in the pipeline, we need the reference model EPICS records +if not options.factors_from_filters_file: + # Needed for kappa_tst + if not options.no_kappatst or not options.no_kappapu or not options.no_kappac or not options.no_fcc: + 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.no_kappac or not options.no_fcc or not options.no_kappapu: + 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.no_kappac or not options.no_fcc: + 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.remove_callines and remove_esd_act_line: + 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.no_fs or not options.no_srcQ: + 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") + +# We also need excitation channels for computing kappas +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 or (options.remove_callines and remove_esd_act_line): + channel_list.append((instrument, options.tst_exc_channel_name)) + headkeys.append("tstexc") +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") + +# We need to make sure we have DARM_ERR and the PCAL channel for computing \kappas +if not options.no_kappac or not options.no_fcc or not options.no_kappatst or not options.no_kappapu or not options.no_fs or not options.no_srcQ: + channel_list.append((instrument, options.pcal_channel_name)) + headkeys.append("pcal") + if options.partial_calibration: + channel_list.append((instrument, options.darm_err_channel_name)) + headkeys.append("darm_err") + +# For full calibration we need DARM_ERR and DARM_CTRL as our input channels +if options.full_calibration: + channel_list.extend(((instrument, options.darm_err_channel_name), (instrument, options.darm_ctrl_channel_name))) + headkeys.extend(("res", "ctrl")) +# For partial calibration we need DELTAL_TST, DELTAL_PUM, DELTAL_UIM, and DELTAL_RES +elif options.partial_calibration: + channel_list.extend(((instrument, options.deltal_res_channel_name), (instrument, options.deltal_tst_channel_name), (instrument, options.deltal_pum_channel_name), (instrument, options.deltal_uim_channel_name))) + headkeys.extend(("res", "tst", "pum", "uim")) + +# If we are removing additional noise from the spectrum (beam jitter, angular control, 60 Hz lines, etc.), we need more channels +if options.remove_powerlines: + channel_list.append((instrument, options.powerlines_channel_name)) + headkeys.append("powerlines") +if options.remove_jitter_imc: + channel_list.extend(((instrument, options.imc_a_pitch_channel_name), (instrument, options.imc_a_yaw_channel_name), (instrument, options.imc_b_pitch_channel_name), (instrument, options.imc_b_yaw_channel_name))) + headkeys.extend(("imc_a_pitch", "imc_a_yaw", "imc_b_pitch", "imc_b_yaw")) +if options.remove_jitter_psl: + channel_list.extend(((instrument, options.bullseye_width_channel_name), (instrument, options.bullseye_pitch_channel_name), (instrument, options.bullseye_yaw_channel_name))) + headkeys.extend(("bullseye_width", "bullseye_pitch", "bullseye_yaw")) +if options.remove_angular_control: + channel_list.extend(((instrument, options.asc_dhard_pitch_channel_name), (instrument, options.asc_dhard_yaw_channel_name), (instrument, options.asc_chard_pitch_channel_name), (instrument, options.asc_chard_yaw_channel_name))) + headkeys.extend(("asc_dhard_pitch", "asc_dhard_yaw", "asc_chard_pitch", "asc_chard_yaw")) +if options.remove_length_control: + channel_list.extend(((instrument, options.lsc_srcl_channel_name), (instrument, options.lsc_mich_channel_name), (instrument, options.lsc_prcl_channel_name))) + headkeys.extend(("lsc_srcl", "lsc_mich", "lsc_prcl")) + + +#################################################################################################### +####################################### Main Pipeline ############################################## +#################################################################################################### + +pipeline = Gst.Pipeline(name="gstlal_compute_strain") +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) + +# When reading from disk, clip the incoming data stream(s) to segment list if one is provided +if options.data_source == "frames" and frame_segments is not None: + for key in headkeys: + currenthead = head_dict[key] + head_dict[key] = calibration_parts.mkgate(pipeline, currenthead, pipeparts.mksegmentsrc(pipeline, frame_segments[instrument]), 1, long_queue, long_queue) + +# +# TIME-VARYING FACTORS COMPUTATIONS +# + +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) + +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) + 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) + +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 or (options.remove_callines and remove_esd_act_line): + 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_fs or not options.no_srcQ: + exc = calibration_parts.caps_and_progress(pipeline, head_dict["exc"], hoft_caps, "exc") + +# If we use the coherence channels multiple times in the pipeline, we need to tee the channels +if not options.no_coherence: + lowfreq_coh_use = int(not options.no_kappatst) + int(not options.no_kappapu) + int(not options.no_kappac) + int(not options.no_fcc) + int(not options.no_srcQ) + int(not options.no_fs) + int(not options.no_dq_vector and (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)) + highfreq_coh_use = int(not options.no_kappac) + int(not options.no_fcc) + int(not options.no_srcQ) + int(not options.no_fs) + int(not options.no_dq_vector and (not options.no_kappac or not options.no_fcc or not options.no_srcQ or not options.no_fs)) + if lowfreq_coh_use >= 2: + pcaly_line1_coh = pipeparts.mktee(pipeline, pcaly_line1_coh) + sus_coh = pipeparts.mktee(pipeline, sus_coh) + darm_coh = pipeparts.mktee(pipeline, darm_coh) + if highfreq_coh_use >= 2: + pcaly_line2_coh = pipeparts.mktee(pipeline, pcaly_line2_coh) + +# Set up computations for \kappa_a, \kappa_tst,\kappa_c, \kappa_pu, f_cc, if applicable +if not options.no_kappac or not options.no_fcc or not options.no_kappatst or not options.no_kappapu or not options.no_srcQ or not options.no_fs: + + # pcal excitation channel, which will be demodulated + pcal = calibration_parts.caps_and_progress(pipeline, head_dict["pcal"], hoft_caps, "pcal") + pcaltee = pipeparts.mktee(pipeline, pcal) + + # DARM_ERR channel, which will have followed different paths if we're doing full vs. partial calibration + if options.full_calibration: + darm_err = calibration_parts.caps_and_progress(pipeline, head_dict["res"], hoft_caps, "darm_err") + else: + darm_err = calibration_parts.caps_and_progress(pipeline, head_dict["darm_err"], hoft_caps, "darm_err") + derrtee = pipeparts.mktee(pipeline, darm_err) + + # 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, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_darm_act_freq_real, 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 or options.remove_callines: + 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, td, compute_calib_factors_complex_caps, integration_samples) + if options.dewhitening: + # dewhiten DARM_ERR at the DARM actuation line frequency + derr_at_darm_act_freq = calibration_parts.complex_audioamplify(pipeline, derr_at_darm_act_freq, derr_dewhiten_at_darm_act_freq_real, derr_dewhiten_at_darm_act_freq_imag) + if 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, td, compute_calib_factors_complex_caps, integration_samples) + if options.remove_callines and remove_esd_act_line: + tstexc_at_esd_act_freq = pipeparts.mktee(pipeline, tstexc_at_esd_act_freq) + + # demodulate DARM_ERR at the ESD actuation line frequency + derr_at_esd_act_freq = calibration_parts.demodulate(pipeline, derrtee, esd_act_line_freq, td, compute_calib_factors_complex_caps, integration_samples) + if options.dewhitening: + # dewhiten DARM_ERR at the ESD actuation line frequency + derr_at_esd_act_freq = calibration_parts.complex_audioamplify(pipeline, derr_at_esd_act_freq, derr_dewhiten_at_esd_act_freq_real, derr_dewhiten_at_esd_act_freq_imag) + + # 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"], long_queue, short_queue) + 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, long_queue, short_queue) + 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, long_queue, short_queue) + + if not options.no_kappatst or not options.no_kappac or not options.no_fcc or not options.no_srcQ or not options.no_fs: + ktst = pipeparts.mktee(pipeline, ktst) + 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) + 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, short_queue, long_queue, attack_length = -integration_samples, 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, short_queue, long_queue, attack_length = -integration_samples, 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, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) + + if not options.no_dq_vector: + ktst_gated = pipeparts.mktee(pipeline, ktst_gated) + smooth_ktstRdq, smooth_ktstIdq = calibration_parts.track_bad_complex_kappas(pipeline, ktst_gated, options.expected_kappatst_real, options.expected_kappatst_imag, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) + + # Smooth kappa_tst + smooth_ktstR, smooth_ktstI = 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) + + else: + if not options.no_dq_vector: + smooth_ktstRdq, smooth_ktstIdq = calibration_parts.track_bad_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) + + # Smooth kappa_tst + smooth_ktstR, smooth_ktstI = 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) + + smooth_ktstRtee = pipeparts.mktee(pipeline, smooth_ktstR) + smooth_ktstItee = pipeparts.mktee(pipeline, smooth_ktstI) + +# If we're also computing \kappa_c, f_cc, or \kappa_pu, keep going +if not options.no_kappac or not options.no_fcc or not options.no_kappapu or not options.no_srcQ or not options.no_fs: + # demodulate excitation channel at PU actuation line frequency + exc_at_pu_act_freq = calibration_parts.demodulate(pipeline, exc, pu_act_esd_line_freq, td, compute_calib_factors_complex_caps, integration_samples) + + # demodulate DARM_ERR at PU actuation line frequency + derr_at_pu_act_freq = calibration_parts.demodulate(pipeline, derrtee, pu_act_esd_line_freq, td, compute_calib_factors_complex_caps, integration_samples) + if options.dewhitening: + # dewhiten DARM_ERR at the PU actuation line frequency + derr_at_pu_act_freq = calibration_parts.complex_audioamplify(pipeline, derr_at_pu_act_freq, derr_dewhiten_at_pu_act_freq_real, derr_dewhiten_at_pu_act_freq_imag) + + # 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"], long_queue, short_queue) + EP3 = calibration_parts.merge_into_complex(pipeline, head_dict["EP3_real"], head_dict["EP3_imag"], long_queue, short_queue) + EP4 = calibration_parts.merge_into_complex(pipeline, head_dict["EP4_real"], head_dict["EP4_imag"], long_queue, short_queue) + 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, long_queue, short_queue) + 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, long_queue, short_queue) + + # \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, long_queue, short_queue) + 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, long_queue, short_queue) + + if not options.no_kappapu or not options.no_srcQ or not options.no_fs: + kpu = pipeparts.mktee(pipeline, kpu) + 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) + 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, short_queue, long_queue, attack_length = -integration_samples, 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, short_queue, long_queue, attack_length = -integration_samples, 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, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) + + if not options.no_dq_vector: + kpu_gated = pipeparts.mktee(pipeline, kpu_gated) + smooth_kpuRdq, smooth_kpuIdq = calibration_parts.track_bad_complex_kappas(pipeline, kpu_gated, options.expected_kappapu_real, options.expected_kappapu_imag, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) + + # Smooth kappa_pu + smooth_kpuR, smooth_kpuI = 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) + + else: + if not options.no_dq_vector: + smooth_kpuRdq, smooth_kpuIdq = calibration_parts.track_bad_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) + + # Smooth kappa_pu + smooth_kpuR, smooth_kpuI = 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) + + 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 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, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_opt_gain_fcc_freq_real, pcal_corr_at_opt_gain_fcc_freq_imag) + if options.remove_callines: + pcal_at_opt_gain_freq = pipeparts.mktee(pipeline, pcal_at_opt_gain_freq) + + # 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, td, compute_calib_factors_complex_caps, integration_samples) + if options.dewhitening: + # dewhiten DARM_ERR at optical gain and f_cc line frequency + derr_at_opt_gain_freq = calibration_parts.complex_audioamplify(pipeline, derr_at_opt_gain_freq, derr_dewhiten_at_opt_gain_fcc_freq_real, derr_dewhiten_at_opt_gain_fcc_freq_imag) + + # Compute the factor S which will be used for the kappa_c and f_cc calculations + if not options.factors_from_filters_file: + EP6 = calibration_parts.merge_into_complex(pipeline, head_dict["EP6_real"], head_dict["EP6_imag"], long_queue, short_queue) + EP7 = calibration_parts.merge_into_complex(pipeline, head_dict["EP7_real"], head_dict["EP7_imag"], long_queue, short_queue) + EP8 = calibration_parts.merge_into_complex(pipeline, head_dict["EP8_real"], head_dict["EP8_imag"], long_queue, short_queue) + EP9 = calibration_parts.merge_into_complex(pipeline, head_dict["EP9_real"], head_dict["EP9_imag"], long_queue, short_queue) + S = calibration_parts.compute_S(pipeline, EP6, pcal_at_opt_gain_freq, derr_at_opt_gain_freq, EP7, ktst, EP8, kpu, EP9, long_queue, short_queue) + 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, EP8_real, EP8_imag, kpu, EP9_real, EP9_imag, long_queue, short_queue) + + 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, long_queue, short_queue) + 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) + + if not options.no_coherence: + # Gate kappa_c with all four of the calibration lines + kc_gated = calibration_parts.mkgate(pipeline, kc, pcaly_line2_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) + kc_gated = calibration_parts.mkgate(pipeline, kc_gated, darm_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) + kc_gated = calibration_parts.mkgate(pipeline, kc_gated, pcaly_line1_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) + kc_gated = calibration_parts.mkgate(pipeline, kc_gated, sus_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) + + if not options.no_dq_vector: + kc_gated = pipeparts.mktee(pipeline, kc_gated) + smooth_kcdq = calibration_parts.track_bad_kappas(pipeline, kc_gated, options.expected_kappac, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) + + smooth_kc = calibration_parts.smooth_kappas(pipeline, kc_gated, options.expected_kappac, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) + + else: + 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) + + if not options.no_dq_vector: + smooth_kcdq = calibration_parts.track_bad_kappas_no_coherence(pipeline, kc, options.kappac_ok_var, options.expected_kappac, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) + + smooth_kctee = pipeparts.mktee(pipeline, smooth_kc) + + # compute f_cc + if not options.no_fcc or not options.nosrcQ or not options.no_fs: + fcc = calibration_parts.compute_fcc(pipeline, SR, SI, opt_gain_fcc_line_freq, long_queue, short_queue) + 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) + + 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, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) + fcc_gated = calibration_parts.mkgate(pipeline, fcc_gated, darm_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) + fcc_gated = calibration_parts.mkgate(pipeline, fcc_gated, pcaly_line1_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) + fcc_gated = calibration_parts.mkgate(pipeline, fcc_gated, sus_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) + + if not options.no_dq_vector: + fcc_gated = pipeparts.mktee(pipeline, fcc_gated) + smooth_fccdq = calibration_parts.track_bad_kappas(pipeline, fcc_gated, fcc_default, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) + + smooth_fcc = calibration_parts.smooth_kappas(pipeline, fcc_gated, fcc_default, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) + else: + 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) + if not options.no_dq_vector: + smooth_fccdq = calibration_parts.track_bad_kappas_no_coherence(pipeline, fcc, options.fcc_ok_var, fcc_default, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) + + smooth_fcctee = pipeparts.mktee(pipeline, smooth_fcc) + + if not options.no_fs or not options.no_srcQ: + # demodulate PCAL channel and apply the PCAL correction factor at SRC detuning line frequency + pcal_at_src_freq = calibration_parts.demodulate(pipeline, pcaltee, src_pcal_line_freq, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_src_freq_real, pcal_corr_at_src_freq_imag) + if options.remove_callines and remove_src_pcal_line: + pcal_at_src_freq = pipeparts.mktee(pipeline, pcal_at_src_freq) + + # demodulate DARM_ERR at SRC detuning line frequency + derr_at_src_freq = calibration_parts.demodulate(pipeline, derrtee, src_pcal_line_freq, td, compute_calib_factors_complex_caps, integration_samples) + + # Compute the factor Xi which will be used for the f_s and src_Q calculations + if not options.factors_from_filters_file: + EP11 = calibration_parts.merge_into_complex(pipeline, head_dict["EP11_real"], head_dict["EP11_imag"], long_queue, short_queue) + EP12 = calibration_parts.merge_into_complex(pipeline, head_dict["EP12_real"], head_dict["EP12_imag"], long_queue, short_queue) + EP13 = calibration_parts.merge_into_complex(pipeline, head_dict["EP13_real"], head_dict["EP13_imag"], long_queue, short_queue) + EP14 = calibration_parts.merge_into_complex(pipeline, head_dict["EP14_real"], head_dict["EP14_imag"], long_queue, short_queue) + Xi = calibration_parts.compute_Xi(pipeline, pcal_at_src_freq, derr_at_src_freq, src_pcal_line_freq, EP11, EP12, EP13, EP14, ktst, kpu, kc, fcc, long_queue, short_queue) + 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, kpu, kc, fcc, long_queue, short_queue) + + XiR, XiI = calibration_parts.split_into_real(pipeline, Xi) + + if options.no_srcQ: + # the imaginary part is only used to compute Q + pipeparts.mkfakesink(pipeline, XiI) + + sqrtXiR = pipeparts.mkpow(pipeline, XiR, exponent = 0.5) + if not options.no_fs and not options.no_srcQ: + sqrtXiR = pipeparts.mktee(pipeline, sqrtXiR) + + # compute f_s + if not options.no_fs: + fs = pipeparts.mkaudioamplify(pipeline, sqrtXiR, src_pcal_line_freq) + fs = pipeparts.mktee(pipeline, fs) + smooth_fs_nogate = pipeparts.mkgeneric(pipeline, fs, "lal_smoothkappas", default_kappa_re = options.expected_fs, array_size = median_smoothing_samples, avg_array_size = factors_average_samples, default_to_median = options.kappas_default_to_median) + + if not options.no_coherence: + # Gate f_s with all four of the calibration lines + fs_gated = calibration_parts.mkgate(pipeline, fs, pcaly_line2_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) + fs_gated = calibration_parts.mkgate(pipeline, fs_gated, darm_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) + fs_gated = calibration_parts.mkgate(pipeline, fs_gated, pcaly_line1_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) + fs_gated = calibration_parts.mkgate(pipeline, fs_gated, sus_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) + + if not options.no_dq_vector: + fs_gated = pipeparts.mktee(pipeline, fs_gated) + smooth_fsdq = calibration_parts.track_bad_kappas(pipeline, fs_gated, options.expected_fs, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) + + smooth_fs = calibration_parts.smooth_kappas(pipeline, fs_gated, options.expected_fs, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) + + else: + smooth_fs = calibration_parts.smooth_kappas_no_coherence(pipeline, fs, options.fs_ok_var, options.expected_fs, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) + + if not options.no_dq_vector: + smooth_fsdq = calibration_parts.track_bad_kappas_no_coherence(pipeline, fs, options.fs_ok_var, options.expected_fs, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) + + if not options.no_dq_vector: + smooth_fs = pipeparts.mktee(pipeline, smooth_fs) + + # compute SRC Q_inv + if not options.no_srcQ: + sqrtXiR_inv = pipeparts.mkpow(pipeline, sqrtXiR, exponent = -1.0) + sqrtXiR_inv = pipeparts.mkaudioamplify(pipeline, sqrtXiR_inv, -1.0) + srcQ_inv = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, [sqrtXiR_inv, long_queue], [XiI, short_queue])) + srcQ_inv = pipeparts.mktee(pipeline, srcQ_inv) + smooth_srcQ_inv_nogate = pipeparts.mkgeneric(pipeline, srcQ_inv, "lal_smoothkappas", default_kappa_re = 1.0 / options.expected_srcQ, array_size = median_smoothing_samples, avg_array_size = factors_average_samples, default_to_median = options.kappas_default_to_median) + + if not options.no_coherence: + # Gate SRC_Q with all four of the calibration lines + srcQ_inv_gated = calibration_parts.mkgate(pipeline, srcQ_inv, pcaly_line2_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) + srcQ_inv_gated = calibration_parts.mkgate(pipeline, srcQ_inv_gated, darm_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) + srcQ_inv_gated = calibration_parts.mkgate(pipeline, srcQ_inv_gated, pcaly_line1_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) + srcQ_inv_gated = calibration_parts.mkgate(pipeline, srcQ_inv_gated, sus_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) + + if not options.no_dq_vector: + srcQ_inv_gated = pipeparts.mktee(pipeline, srcQ_inv_gated) + smooth_srcQdq = calibration_parts.track_bad_kappas(pipeline, srcQ_inv_gated, 1.0 / options.expected_srcQ, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) + + smooth_srcQ_inv = calibration_parts.smooth_kappas(pipeline, srcQ_inv_gated, 1.0 / options.expected_srcQ, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) + + else: + smooth_srcQ_inv = calibration_parts.smooth_kappas_no_coherence(pipeline, srcQ_inv, options.srcQ_ok_var, 1.0 / options.expected_srcQ, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) + + if not options.no_dq_vector: + smooth_srcQdq = calibration_parts.track_bad_kappas_no_coherence(pipeline, srcQ_inv, options.srcQ_ok_var, 1.0 / options.expected_srcQ, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) + + if not options.no_dq_vector: + smooth_srcQ_inv = pipeparts.mktee(pipeline, smooth_srcQ_inv) + +# +# Calibration Line Removal +# + +if options.remove_callines: + # if we didn't compute the kappas, we still need to get the pcal channel + if options.no_kappatst and options.no_kappapu and options.no_kappac and options.no_fcc and options.no_srcQ and options.no_fs: + pcal = calibration_parts.caps_and_progress(pipeline, head_dict["pcal"], hoft_caps, "pcal") + pcaltee = pipeparts.mktee(pipeline, pcal) + # Demodulate pcal at the ~30 Hz pcal line + pcal_at_darm_act_freq = calibration_parts.demodulate(pipeline, pcal_tee, darm_act_line_freq, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_darm_act_freq_real, pcal_corr_at_darm_act_freq_imag) + + # Reconstruct a calibrated (negative) pcal at only the ~30 Hz pcal line + pcaly_line1 = calibration_parts.mkresample(pipeline, pcal_at_darm_act_freq, 3, False, "audio/x-raw, format=Z128LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) + pcaly_line1 = pipeparts.mkgeneric(pipeline, pcaly_line1, "lal_demodulate", line_frequency = -1.0 * darm_act_line_freq, prefactor_real = 2.0) + remove_pcaly_line1, trash = calibration_parts.split_into_real(pipeline, pcaly_line1) + pipeparts.mkfakesink(pipeline, trash) + + # 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, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_opt_gain_fcc_freq_real, pcal_corr_at_opt_gain_fcc_freq_imag) + # Reconstruct a calibrated (negative) pcal at only the ~300 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, trash = calibration_parts.split_into_real(pipeline, pcaly_line2) + pipeparts.mkfakesink(pipeline, trash) + + # Add the first two components together. We will add this to h(t) to remove these lines + remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_pcaly_line1, long_queue], [remove_pcaly_line2, short_queue])) + + if remove_esd_act_line: + # Make sure we have demodulated the ESD excitation channel at the ~30 Hz ESD line + if options.no_kappac and options.no_fcc and options.no_kappatst and options.no_kappapu and options.no_srcQ and options.no_fs: + tstexc_at_esd_act_freq = calibration_parts.demodulate(pipeline, tstexc, esd_act_line_freq, td, compute_calib_factors_complex_caps, integration_samples) + if options.factors_from_filters_file: + esd_act_line = calibration_parts.complex_audioamplify(pipeline, tstexc_at_esd_act_freq, EP10_real, EP10_imag) + else: + # EP10 was read from the frames + EP10 = calibration_parts.merge_into_complex(pipeline, head_dict["EP10_real"], head_dict["EP10_imag"], long_queue, short_queue) + esd_act_line = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, [tstexc_at_esd_act_freq, long_queue], [EP10, short_queue])) + # Reconstruct a calibrated (negative) ESD injection at the ~30 Hz ESD line + if options.apply_kappatst: + # Multiply by the real part of kappa_tst + esd_act_line = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, [esd_act_line, long_queue], [pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, smooth_ktstRtee, matrix=[[1.0, 0.0]])), short_queue])) + 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, trash = calibration_parts.split_into_real(pipeline, esd_act_line_remove) + pipeparts.mkfakesink(pipeline, trash) + # Add into the total line removal stream + remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_from_strain, short_queue], [esd_act_line_remove, long_queue])) + + if remove_high_pcal_line: + # Demodulate pcal at the ~1kHz pcal line + pcaly_line3 = calibration_parts.demodulate(pipeline, pcaltee, high_pcal_line_freq, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_high_line_freq_real, pcal_corr_at_high_line_freq_imag) + # Reconstruct a calibrated (negative) pcal at only the ~1kHz pcal line + pcaly_line3 = calibration_parts.mkresample(pipeline, pcaly_line3, 3, False, "audio/x-raw, format=Z128LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) + pcaly_line3 = pipeparts.mkgeneric(pipeline, pcaly_line3, "lal_demodulate", line_frequency = -1.0 * high_pcal_line_freq, prefactor_real = 2.0) + remove_pcaly_line3, trash = calibration_parts.split_into_real(pipeline, pcaly_line3) + pipeparts.mkfakesink(pipeline, trash) + # Add into the total line removal stream + remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_from_strain, long_queue], [remove_pcaly_line3, short_queue])) + + if remove_roaming_pcal_line: + # Demodulate pcal at the ~3kHz pcal line + pcaly_line4 = calibration_parts.demodulate(pipeline, pcaltee, roaming_pcal_line_freq, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_roaming_line_real, pcal_corr_at_roaming_line_imag) + # Reconstruct a calibrated (negative) pcal at only the ~3kHz pcal line + pcaly_line4 = calibration_parts.mkresample(pipeline, pcaly_line4, 3, False, "audio/x-raw, format=Z128LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) + pcaly_line4 = pipeparts.mkgeneric(pipeline, pcaly_line4, "lal_demodulate", line_frequency = -1.0 * roaming_pcal_line_freq, prefactor_real = 2.0) + remove_pcaly_line4, trash = calibration_parts.split_into_real(pipeline, pcaly_line4) + pipeparts.mkfakesink(pipeline, trash) + # Add into the total line removal stream + remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_from_strain, long_queue], [remove_pcaly_line4, short_queue])) + + if remove_src_pcal_line: + # Make sure we have demodulated pcal at the ~8 Hz pcal line + if options.no_fs and options.no_srcQ: + pcal_at_src_freq = calibration_parts.demodulate(pipeline, pcaltee, src_pcal_line_freq, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_src_freq_real, pcal_corr_at_src_freq_imag) + # Reconstruct a calibrated (negative) pcal at only the ~3kHz pcal line + pcaly_line0 = calibration_parts.mkresample(pipeline, pcal_at_src_freq, 3, False, "audio/x-raw, format=Z128LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) + pcaly_line0 = pipeparts.mkgeneric(pipeline, pcaly_line0, "lal_demodulate", line_frequency = -1.0 * src_pcal_line_freq, prefactor_real = 2.0) + remove_pcaly_line0, trash = calibration_parts.split_into_real(pipeline, pcaly_line0) + pipeparts.mkfakesink(pipeline, trash) + # Add into the total line removal stream + remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_from_strain, long_queue], [remove_pcaly_line0, short_queue])) + +if options.remove_powerlines: + powerlines = calibration_parts.caps_and_progress(pipeline, head_dict["powerlines"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % powerlinessr, "powerlines") + powerlines = pipeparts.mkfirbank(pipeline, powerlines, latency = int(powerlinesdelay), fir_matrix = [powerlinesfilt[::-1]], time_domain = td) + powerlines = calibration_parts.mkresample(pipeline, powerlines, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) + if options.remove_callines: + remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_from_strain, long_queue], [powerlines, short_queue])) + else: + remove_from_strain = powerlines + +if options.remove_jitter_imc: + imc_a_pitch = calibration_parts.caps_and_progress(pipeline, head_dict["imc_a_pitch"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % imcapitsr, "imc_a_pitch") + imc_a_yaw = calibration_parts.caps_and_progress(pipeline, head_dict["imc_a_yaw"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % imcayawsr, "imc_a_yaw") + imc_b_pitch = calibration_parts.caps_and_progress(pipeline, head_dict["imc_b_pitch"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % imcbpitsr, "imc_b_pitch") + imc_b_yaw = calibration_parts.caps_and_progress(pipeline, head_dict["imc_b_yaw"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % imcbyawsr, "imc_b_yaw") + + imc_a_pitch = pipeparts.mkfirbank(pipeline, imc_a_pitch, latency = int(imcapitdelay), fir_matrix = [imcapitfilt[::-1]], time_domain = td) + imc_a_yaw = pipeparts.mkfirbank(pipeline, imc_a_yaw, latency = int(imcayawdelay), fir_matrix = [imcayawfilt[::-1]], time_domain = td) + imc_b_pitch = pipeparts.mkfirbank(pipeline, imc_b_pitch, latency = int(imcbpitdelay), fir_matrix = [imcbpitfilt[::-1]], time_domain = td) + imc_b_yaw = pipeparts.mkfirbank(pipeline, imc_b_yaw, latency = int(imcbyawdelay), fir_matrix = [imcbyawfilt[::-1]], time_domain = td) + + imc_a_pitch = calibration_parts.mkresample(pipeline, imc_a_pitch, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) + imc_a_yaw = calibration_parts.mkresample(pipeline, imc_a_yaw, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) + imc_b_pitch = calibration_parts.mkresample(pipeline, imc_b_pitch, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) + imc_b_yaw = calibration_parts.mkresample(pipeline, imc_b_yaw, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) + + if options.remove_callines or options.remove_powerlines: + remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_from_strain, long_queue], [imc_a_pitch, long_queue], [imc_a_yaw, long_queue], [imc_b_pitch, long_queue], [imc_b_yaw, long_queue])) + else: + remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [imc_a_pitch, long_queue], [imc_a_yaw, long_queue], [imc_b_pitch, long_queue], [imc_b_yaw, long_queue])) + +if options.remove_jitter_psl: + bullseye_width = calibration_parts.caps_and_progress(pipeline, head_dict["bullseye_width"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % bullseyewidsr, "bullseye_width") + bullseye_pitch = calibration_parts.caps_and_progress(pipeline, head_dict["bullseye_pitch"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % bullseyepitsr, "bullseye_pitch") + bullseye_yaw = calibration_parts.caps_and_progress(pipeline, head_dict["bullseye_yaw"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % bullseyeyawsr, "bullseye_yaw") + + bullseye_width = pipeparts.mkfirbank(pipeline, bullseye_width, latency = int(bullseyewiddelay), fir_matrix = [bullseyewidfilt[::-1]], time_domain = td) + bullseye_pitch = pipeparts.mkfirbank(pipeline, bullseye_pitch, latency = int(bullseyepitdelay), fir_matrix = [bullseyepitfilt[::-1]], time_domain = td) + bullseye_yaw = pipeparts.mkfirbank(pipeline, bullseye_yaw, latency = int(bullseyeyawdelay), fir_matrix = [bullseyeyawfilt[::-1]], time_domain = td) + + bullseye_width = calibration_parts.mkresample(pipeline, bullseye_width, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) + bullseye_pitch = calibration_parts.mkresample(pipeline, bullseye_pitch, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) + bullseye_yaw = calibration_parts.mkresample(pipeline, bullseye_yaw, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) + + if options.remove_callines or options.remove_powerlines or options.remove_jitter_imc: + remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_from_strain, long_queue], [bullseye_width, long_queue], [bullseye_pitch, long_queue], [bullseye_yaw, long_queue])) + else: + remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [bullseye_width, long_queue], [bullseye_pitch, long_queue], [bullseye_yaw, long_queue])) + +if options.remove_angular_control: + asc_dhard_pitch = calibration_parts.caps_and_progress(pipeline, head_dict["asc_dhard_pitch"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % ascdpitsr, "asc_dhard_pitch") + asc_dhard_yaw = calibration_parts.caps_and_progress(pipeline, head_dict["asc_dhard_yaw"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % ascdyawsr, "asc_dhard_yaw") + asc_chard_pitch = calibration_parts.caps_and_progress(pipeline, head_dict["asc_chard_pitch"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % asccpitsr, "asc_chard_pitch") + asc_chard_yaw = calibration_parts.caps_and_progress(pipeline, head_dict["asc_chard_yaw"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % asccyawsr, "asc_chard_yaw") + + asc_dhard_pitch = pipeparts.mkfirbank(pipeline, asc_dhard_pitch, latency = int(ascdpitdelay), fir_matrix = [ascdpitfilt[::-1]], time_domain = td) + asc_dhard_yaw = pipeparts.mkfirbank(pipeline, asc_dhard_yaw, latency = int(ascdyawdelay), fir_matrix = [ascdyawfilt[::-1]], time_domain = td) + asc_chard_pitch = pipeparts.mkfirbank(pipeline, asc_chard_pitch, latency = int(asccpitdelay), fir_matrix = [asccpitfilt[::-1]], time_domain = td) + asc_chard_yaw = pipeparts.mkfirbank(pipeline, asc_chard_yaw, latency = int(asccyawdelay), fir_matrix = [asccyawfilt[::-1]], time_domain = td) + + asc_dhard_pitch = calibration_parts.mkresample(pipeline, asc_dhard_pitch, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) + asc_dhard_yaw = calibration_parts.mkresample(pipeline, asc_dhard_yaw, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) + asc_chard_pitch = calibration_parts.mkresample(pipeline, asc_chard_pitch, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) + asc_chard_yaw = calibration_parts.mkresample(pipeline, asc_chard_yaw, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) + + if options.remove_callines or options.remove_powerlines or options.remove_jitter_imc or options.remove_jitter_psl: + remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_from_strain, long_queue], [asc_dhard_pitch, long_queue], [asc_dhard_yaw, long_queue], [asc_chard_pitch, long_queue], [asc_chard_yaw, long_queue])) + else: + remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [asc_dhard_pitch, long_queue], [asc_dhard_yaw, long_queue], [asc_chard_pitch, long_queue], [asc_chard_yaw, long_queue])) + +if options.remove_length_control: + lsc_srcl = calibration_parts.caps_and_progress(pipeline, head_dict["lsc_srcl"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % lscsrclsr, "lsc_srcl") + lsc_mich = calibration_parts.caps_and_progress(pipeline, head_dict["lsc_mich"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % lscmichsr, "lsc_mich") + lsc_prcl = calibration_parts.caps_and_progress(pipeline, head_dict["lsc_prcl"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % lscprclsr, "lsc_prcl") + + lsc_srcl = pipeparts.mkfirbank(pipeline, lsc_srcl, latency = int(lscsrcldelay), fir_matrix = [lscsrclfilt[::-1]], time_domain = td) + lsc_mich = pipeparts.mkfirbank(pipeline, lsc_mich, latency = int(lscmichdelay), fir_matrix = [lscmichfilt[::-1]], time_domain = td) + lsc_prcl = pipeparts.mkfirbank(pipeline, lsc_prcl, latency = int(lscprcldelay), fir_matrix = [lscprclfilt[::-1]], time_domain = td) + + lsc_srcl = calibration_parts.mkresample(pipeline, lsc_srcl, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) + lsc_mich = calibration_parts.mkresample(pipeline, lsc_mich, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) + lsc_prcl = calibration_parts.mkresample(pipeline, lsc_prcl, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) + + if options.remove_callines or options.remove_powerlines or options.remove_jitter_imc or options.remove_jitter_psl or options.remove_angular_control: + remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_from_strain, long_queue], [lsc_srcl, long_queue], [lsc_mich, long_queue], [lsc_prcl, long_queue])) + else: + remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [lsc_srcl, long_queue], [lsc_mich, long_queue], [lsc_prcl, long_queue])) + +# +# CONTROL BRANCH +# + +# zero out filter settling samples +tst_filter_settle_time = 0.0 +tst_filter_latency = 0.0 +pumuim_filter_settle_time = 0.0 +pumuim_filter_latency = 0.0 + +# The reverse of the filters will be used in all filtering below due to the definition of the filtering procedure employed by lal_firbank +if options.partial_calibration: + # enforce caps on actuation channels and set up progress report if verbose is on + tst = calibration_parts.caps_and_progress(pipeline, head_dict["tst"], ctrl_caps, "tst") + tsttee = pipeparts.mktee(pipeline, tst) + pum = calibration_parts.caps_and_progress(pipeline, head_dict["pum"], ctrl_caps, "pum") + pumtee = pipeparts.mktee(pipeline, pum) + uim = calibration_parts.caps_and_progress(pipeline, head_dict["uim"], ctrl_caps, "uim") + uimtee = pipeparts.mktee(pipeline, uim) + + # add together the PUM and UIM actuation channels; this may change in the future... + pumuim = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [pumtee, long_queue], [uimtee, short_queue])) + + # if you need to, dewhiten the TST and PUM/UIM chains + if options.dewhitening: + pumuim = calibration_parts.mkstockresample(pipeline, pumuim, "audio/x-raw, format=F64LE, rate=%d" % pumuimdewhitensr) + pumuim = pipeparts.mkfirbank(pipeline, pumuim, latency = int(pumuimdewhitendelay), fir_matrix = [pumuimdewhiten[::-1]], time_domain = td) + pumuim_filter_settle_time += float(len(pumuimdewhiten)-pumuimdewhitendelay)/pumuimdewhitensr + pumuim_filter_latency += float(pumuimdewhitendelay)/pumuimdewhitensr + tst = calibration_parts.mkstockresample(pipeline, tsttee, "audio/x-raw, format=F64LE, rate=%d" % tstdewhitensr) + tst = pipeparts.mkfirbank(pipeline, tst, latency = int(tstdewhitendelay), fir_matrix = [tstdewhiten[::-1]], time_domain = td) + tst_filter_settle_time += float(len(tstdewhiten)-tstdewhitendelay)/tstdewhitensr + tst_filter_latency += float(tstdewhitendelay)/tstdewhitensr + else: + tst = tsttee + +if options.full_calibration: + # enforce caps on actuation channels and set up progress report, if verbose is on + ctrl = calibration_parts.caps_and_progress(pipeline, head_dict["ctrl"], hoft_caps, "ctrl") + darmctrltee = pipeparts.mktee(pipeline, ctrl) + + if options.dewhitening: + # dewhiten the DARM_CTRL channel + ctrl = calibration_parts.mkstockresample(pipeline, darmctrltee, "audio/x-raw, format=F64LE, rate=%d" % ctrldewhitensr) + ctrl = pipeparts.mkfirbank(pipeline, ctrl, latency = int(ctrldewhitendelay), fir_matrix = [ctrldewhiten[::-1]], time_domain = td) + tst_filter_settle_time += float(len(ctrldewhiten)-ctrldewhitendelay)/ctrldewhitensr + tst_filter_latency += float(ctrldewhitendelay)/ctrldewhitensr + pumuim_filter_settle_time += float(len(ctrldewhiten)-ctrldewhitendelay)/ctrldewhitensr + pumuim_filter_latency += float(ctrldewhitendelay)/ctrldewhitensr + # tee off DARM_CTRL to be filtered with PUM/UIM and TST filters separately + ctrltee = pipeparts.mktee(pipeline, ctrl) + else: + ctrltee = pipeparts.mktee(pipeline, darmctrltee) + tst = ctrltee + pumuim = ctrltee + +# resample what will become the TST actuation chain to the TST FIR filter sample rate +tst = calibration_parts.mkstockresample(pipeline, tst, "audio/x-raw, format=F64LE, rate=%d" % tstchainsr) +# filter TST chain with the TST acutaiton filter +tst = pipeparts.mkfirbank(pipeline, tst, latency = int(tstdelay), fir_matrix = [tstfilt[::-1]], time_domain = td) +if options.low_latency: + tst = calibration_parts.mkinsertgap(pipeline, tst, bad_data_intervals = [-1e35, 1e35], block_duration = 1000000000 * options.buffer_length, remove_gap = False) +tst_filter_settle_time += float(len(tstfilt)-tstdelay)/tstchainsr +tst_filter_latency += float(tstdelay)/tstchainsr +# resample the TST actuation chain to the full sample rate +if tstchainsr != pumuimchainsr or options.apply_kappatst or options.apply_kappapu: + tst = calibration_parts.mkstockresample(pipeline, tst, hoft_caps) + +# resample what will become the PUM/UIM actuation chain to the PUM/UIM FIR filter sample rate +pumuim = calibration_parts.mkstockresample(pipeline, pumuim, "audio/x-raw, format=F64LE, rate=%d" % pumuimchainsr) +# filter the PUM/UIM chain with the PUM/UIM actuation filter +pumuim = pipeparts.mkfirbank(pipeline, pumuim, latency = int(pumuimdelay), fir_matrix = [pumuimfilt[::-1]], time_domain = td) +if options.low_latency: + pumuim = calibration_parts.mkinsertgap(pipeline, pumuim, bad_data_intervals = [-1e35, 1e35], block_duration = 1000000000 * options.buffer_length, remove_gap = False) +pumuim_filter_settle_time += float(len(pumuimfilt)-pumuimdelay)/pumuimchainsr +pumuim_filter_latency += float(pumuimdelay)/pumuimchainsr +# resample the PUM/UIM actuation chain to the full sample rate +if tstchainsr != pumuimchainsr or options.apply_kappapu or options.apply_kappatst: + pumuim = calibration_parts.mkstockresample(pipeline, pumuim, hoft_caps) + +# apply kappa_tst +if options.apply_kappatst: + # Only apply the real part of \kappa_tst as a correction to A_tst + ktst_for_tst = calibration_parts.mkresample(pipeline, smooth_ktstRtee, 3, False, hoft_caps) + tst = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, [ktst_for_tst, ktst_queue_length], [tst, tst_queue_length])) +# apply kappa_pu +if options.apply_kappapu: + # Only apply the real part of \kappa_pu as a correction to A_pu + kpu_for_pu = calibration_parts.mkresample(pipeline, smooth_kpuRtee, 3, False, hoft_caps) + pumuim = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, [kpu_for_pu, kpu_queue_length], [pumuim, pumuim_queue_length])) + +# Add the TST and PUM/UIM chains together to form the full actuation chain +ctrl = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [tst, long_queue], [pumuim, short_queue])) +if tstchainsr != hoftsr or not options.apply_kappatst or not options.apply_kappapu: + ctrl = calibration_parts.mkstockresample(pipeline, ctrl, hoft_caps) + +# +# RESIDUAL BRANCH +# + +# zero out res filter settle time +res_filter_settle_time = 0.0 +res_filter_latency = 0.0 + +# The reverse of the filters will be used in all filtering below due to the definition of the filtering procedure employed by lal_firbank + +# enforce caps on the residual branch and hook up progress report if verbose is on +if options.full_calibration: + if not options.no_kappatst or not options.no_kappapu or not options.no_kappac or not options.no_fcc: + restee = derrtee + else: + res = calibration_parts.caps_and_progress(pipeline, head_dict["res"], hoft_caps, "res") + restee = pipeparts.mktee(pipeline, res) +if options.partial_calibration: + res = calibration_parts.caps_and_progress(pipeline, head_dict["res"], hoft_caps, "res") + restee = pipeparts.mktee(pipeline, res) + +# apply the residual chain filter +res = pipeparts.mkfirbank(pipeline, restee, latency = int(reschaindelay), fir_matrix = [reschainfilt[::-1]], time_domain = td) +if options.low_latency: + res = calibration_parts.mkinsertgap(pipeline, res, bad_data_intervals = [-1e35, 1e35], block_duration = 1000000000 * options.buffer_length, remove_gap = False) +res_filter_settle_time += float(len(reschainfilt)-reschaindelay)/hoftsr +res_filter_latency += float(reschaindelay)/hoftsr +if options.dewhitening: + res = pipeparts.mkfirbank(pipeline, res, latency = int(resdewhitendelay), fir_matrix = [resdewhiten[::-1]], time_domain = td) + res_filter_settle_time += float(len(resdewhiten)-resdewhitendelay)/hoftsr + res_filter_latency += float(resdewhitendelay)/hoftsr + +# Apply factors to actuation and sensing chains, if applicable +if options.apply_kappac: + kc_modify_res = calibration_parts.mkresample(pipeline, smooth_kctee, 3, False, hoft_caps) + res = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, [res, reskc_queue_length], [pipeparts.mkpow(pipeline, kc_modify_res, exponent = -1.0), kc_queue_length])) + +filter_settle_time = max(res_filter_settle_time, tst_filter_settle_time, pumuim_filter_settle_time) +filter_latency = max(res_filter_latency, tst_filter_latency, pumuim_filter_latency) + +# +# CONTROL + RESIDUAL = H(T) +# + +# Add control and residual chains and divide by L to make h(t) +strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [res, res_queue_length], [ctrl, ctrl_queue_length])) +# Remove the calibration lines, if we want +if options.remove_callines or options.remove_powerlines or options.remove_jitter_imc or options.remove_jitter_psl or options.remove_angular_control or options.remove_length_control: + remove_from_strain = pipeparts.mkaudioamplify(pipeline, remove_from_strain, -1.0) + strain = pipeparts.mktee(pipeline, strain) + cleaned_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [strain, short_queue], [remove_from_strain, long_queue])) +# Divide by L in a way that is compatitble with old and new filters files, since old filter files don't recored "arm length" +try: + strain = pipeparts.mkaudioamplify(pipeline, strain, 1.0/float(filters["arm_length"])) +except KeyError: + strain = pipeparts.mkaudioamplify(pipeline, strain, 1.0/3994.5) + +strain = pipeparts.mkprogressreport(pipeline, strain, "progress_hoft_%s" % instrument) + +if options.remove_callines: + try: + cleaned_strain = pipeparts.mkaudioamplify(pipeline, cleaned_strain, 1.0/float(filters["arm_length"])) + except KeyError: + cleaned_strain = pipeparts.mkaudioamplify(pipeline, cleaned_strain, 1.0/3994.5) + + cleaned_strain = pipeparts.mkprogressreport(pipeline, cleaned_strain, "progress_hoft_cleaned_%s" % instrument) +# Put the units back to strain before writing to frames +straintagstr = "units=strain,channel-name=%sCALIB_STRAIN%s,instrument=%s" % (chan_prefix, chan_suffix, instrument) +cleaned_straintagstr = "units=strain,channel-name=%sCALIB_STRAIN_CLEAN%s,instrument=%s" % (chan_prefix, chan_suffix, instrument) +if not options.no_dq_vector: + straintee = pipeparts.mktee(pipeline, strain) + strain = pipeparts.mktaginject(pipeline, straintee, straintagstr) +else: + strain = pipeparts.mktaginject(pipeline, strain, straintagstr) +if options.remove_callines: + cleaned_strain = pipeparts.mktaginject(pipeline, cleaned_strain, cleaned_straintagstr) + +# +# CALIB_STATE_VECTOR BRANCH +# + +#FIXME: Add more comments! + +if not options.no_dq_vector: + # FIXME: When the ODC is written as unsigned ints, this piece can be removed + odcstatevector = calibration_parts.caps_and_progress(pipeline, head_dict["odcstatevector"], odc_caps, "odc_%s" % instrument) + odctagstr = "channel-name=%s:%s, instrument=%s" % (instrument, options.dq_channel_name, instrument) + odcstatevector = pipeparts.mktaginject(pipeline, odcstatevector, odctagstr) + odcstatevectortee = pipeparts.mktee(pipeline, odcstatevector) + + # + # GAP BIT BRANCH + # + + nogap = pipeparts.mkbitvectorgen(pipeline, odcstatevectortee, threshold=1, bit_vector = 1) + nogap = pipeparts.mkcapsfilter(pipeline, nogap, odc_caps) + nogap = pipeparts.mkgeneric(pipeline, nogap, "lal_logicalundersample", required_on = 1, status_out = 512) + nogap = pipeparts.mkcapsfilter(pipeline, nogap, calibstate_caps) + + # + # OBSERVATION-INTENT BIT BRANCH + # + + obsintent = pipeparts.mkgeneric(pipeline, odcstatevectortee, "lal_logicalundersample", required_on = options.obs_intent_bitmask, status_out = 2) + obsintent = pipeparts.mkcapsfilter(pipeline, obsintent, calibstate_caps) + obsintenttee = pipeparts.mktee(pipeline, obsintent) + + # + # OBSERVATION-READY BIT BRANCH + # + + obsready = pipeparts.mkgeneric(pipeline, odcstatevectortee, "lal_logicalundersample", required_on = options.obs_ready_bitmask, status_out = 4) + obsready = pipeparts.mkcapsfilter(pipeline, obsready, calibstate_caps) + obsreadytee = pipeparts.mktee(pipeline, obsready) + + # + # H(t)-PRODUCED BIT BRANCH + # + + htproduced = pipeparts.mkbitvectorgen(pipeline, straintee, bit_vector = 8, threshold = 0) + htproduced = pipeparts.mkcapsfilter(pipeline, htproduced, "audio/x-raw, format=U32LE, rate=%d" % hoftsr) + htproduced = pipeparts.mkgeneric(pipeline, htproduced, "lal_logicalundersample", required_on = 8, status_out = 8) + htproduced = pipeparts.mkcapsfilter(pipeline, htproduced, calibstate_caps) + + # + # FILTERS-OK BIT BRANCH + # + + # Set the FILTERS-OK bit based on observation-ready transitions + filtersok = pipeparts.mkbitvectorgen(pipeline, obsintenttee, bit_vector=16, threshold=2) + filtersok = pipeparts.mkcapsfilter(pipeline, filtersok, calibstate_caps) + filtersok = calibration_parts.mkgate(pipeline, filtersok, obsreadytee, 4, long_queue, short_queue, attack_length = -int(filter_settle_time * calibstatesr), hold_length = -int(filter_latency * calibstatesr)) + filtersok = pipeparts.mkbitvectorgen(pipeline, filtersok, bit_vector = 16, nongap_is_control = True) + filtersok = pipeparts.mkcapsfilter(pipeline, filtersok, calibstate_caps) + + # + # NO-INVALID-INPUT BRANCH + # + + # Check if any of the input data channels had to be replaced by zeroes because they were < 1e-35 + resok = pipeparts.mkcapsfilter(pipeline, restee, hoft_caps) + resok = pipeparts.mkbitvectorgen(pipeline, resok, threshold=1e-35, bit_vector=1) + resok = pipeparts.mkcapsfilter(pipeline, resok, "audio/x-raw, format=U32LE, rate=%d" % hoftsr) + resok = pipeparts.mkgeneric(pipeline, resok, "lal_logicalundersample", required_on = 1, status_out = 1) + resok = pipeparts.mkcapsfilter(pipeline, resok, calibstate_caps) + if options.partial_calibration: + tstok = pipeparts.mkcapsfilter(pipeline, tsttee, ctrl_caps) + tstok = pipeparts.mkbitvectorgen(pipeline, tstok, threshold=1e-35, bit_vector=1) + tstok = pipeparts.mkcapsfilter(pipeline, tstok, "audio/x-raw, format=U32LE, rate=%d" % ctrlsr) + tstok = pipeparts.mkgeneric(pipeline, tstok, "lal_logicalundersample", required_on = 1, status_out = 1) + tstok = pipeparts.mkcapsfilter(pipeline, tstok, calibstate_caps) + pumok = pipeparts.mkcapsfilter(pipeline, pumtee, ctrl_caps) + pumok = pipeparts.mkbitvectorgen(pipeline, pumok, threshold=1e-35, bit_vector=1) + pumok = pipeparts.mkcapsfilter(pipeline, pumok, "audio/x-raw, format=U32LE, rate=%d" % ctrlsr) + pumok = pipeparts.mkgeneric(pipeline, pumok, "lal_logicalundersample", required_on = 1, status_out = 1) + pumok = pipeparts.mkcapsfilter(pipeline, pumok, calibstate_caps) + uimok = pipeparts.mkcapsfilter(pipeline, uimtee, ctrl_caps) + uimok = pipeparts.mkbitvectorgen(pipeline, uimok, threshold=1e-35, bit_vector=1) + uimok = pipeparts.mkcapsfilter(pipeline, uimok, "audio/x-raw, format=U32LE, rate=%d" % ctrlsr) + uimok = pipeparts.mkgeneric(pipeline, uimok, "lal_logicalundersample", required_on = 1, status_out = 1) + uimok = pipeparts.mkcapsfilter(pipeline, uimok, calibstate_caps) + noinvalidinput = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [resok, short_queue], [tstok, long_queue], [pumok, long_queue], [uimok, long_queue])) + noinvalidinput = pipeparts.mkbitvectorgen(pipeline, noinvalidinput, threshold=4, bit_vector=33554432) + if options.full_calibration: + ctrlok = pipeparts.mkbitvectorgen(pipeline, darmctrltee, threshold=1e-35, bit_vector=1) + ctrlok = pipeparts.mkcapsfilter(pipeline, ctrlok, "audio/x-raw, format=U32LE, rate=%d" % ctrlsr) + ctrlok = pipeparts.mkgeneric(pipeline, ctrlok, "lal_logicalundersample", required_on = 1, status_out = 1) + ctrlok = pipeparts.mkcapsfilter(pipeline, ctrlok, calibstate_caps) + noinvalidinput = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [resok, long_queue], [ctrlok, long_queue])) + noinvalidinput = pipeparts.mkbitvectorgen(pipeline, noinvalidinput, threshold=2, bit_vector=33554432) + noinvalidinput = pipeparts.mkcapsfilter(pipeline, noinvalidinput, calibstate_caps) + noinvalidinput = pipeparts.mktee(pipeline, noinvalidinput) + # inputs that are replaced with zeros affect h(t) for a short time before and after the zeros, so we also must account for this corrupted time. + noinvalidinput = calibration_parts.mkgate(pipeline, noinvalidinput, noinvalidinput, 33554432, long_queue, short_queue, attack_length = -int(filter_settle_time * calibstatesr), hold_length = -int(filter_latency * calibstatesr)) + + # + # KAPPA-SMOOTHING-SETTLED BIT BRANCH + # + if not options.no_kappac or not options.no_kappatst or not options.no_kappapu or not options.no_fcc: + smoothingok = pipeparts.mkbitvectorgen(pipeline, obsreadytee, bit_vector=1024, threshold=4) + smoothingok = pipeparts.mkcapsfilter(pipeline, smoothingok, calibstate_caps) + smoothingok = calibration_parts.mkgate(pipeline, smoothingok, obsreadytee, 4, long_queue, short_queue, attack_length =-(median_smoothing_samples+factors_average_samples + integration_samples)) + smoothingok = pipeparts.mkbitvectorgen(pipeline, smoothingok, bit_vector = 1024, nongap_is_control = True) + smoothingok = pipeparts.mkcapsfilter(pipeline, smoothingok, calibstate_caps) + + # + # KAPPATST BITS BRANCH + # + if not options.no_kappatst: + ktstSmoothInRange, ktstMedianUncorrupt = calibration_parts.compute_kappa_bits(pipeline, smooth_ktstRtee, smooth_ktstItee, smooth_ktstRdq, smooth_ktstIdq, options.expected_kappatst_real, options.expected_kappatst_imag, options.kappatst_real_ok_var, options.kappatst_imag_ok_var, long_queue, short_queue, status_out_smooth = 2048, status_out_median = 4096, starting_rate = options.compute_factors_sr, ending_rate = calibstatesr) + + # + # KAPPAPU BITS BRANCH + # + if not options.no_kappapu: + kpuSmoothInRange, kpuMedianUncorrupt = calibration_parts.compute_kappa_bits(pipeline, smooth_kpuRtee, smooth_kpuItee, smooth_kpuRdq, smooth_kpuIdq, options.expected_kappapu_real, options.expected_kappapu_imag, options.kappapu_real_ok_var, options.kappapu_imag_ok_var, long_queue, short_queue, status_out_smooth = 8192, status_out_median = 16384, starting_rate = options.compute_factors_sr, ending_rate = calibstatesr) + + # + # KAPPAC BITS BRANCH + # + if not options.no_kappac: + kcSmoothInRange, kcMedianUncorrupt = calibration_parts.compute_kappa_bits_only_real(pipeline, smooth_kctee, smooth_kcdq, options.expected_kappac, options.kappac_ok_var, status_out_smooth = 131072, status_out_median = 262144, starting_rate = options.compute_factors_sr, ending_rate = calibstatesr) + + # + # FCC BITS BRANCH + # + if not options.no_fcc: + fccSmoothInRange, fccMedianUncorrupt = calibration_parts.compute_kappa_bits_only_real(pipeline, smooth_fcctee, smooth_fccdq, fcc_default, options.fcc_ok_var, status_out_smooth = 524288, status_out_median = 1048576, starting_rate = options.compute_factors_sr, ending_rate = calibstatesr) + + # + # FS BITS BRANCH + # + if not options.no_fs: + fsSmoothInRange, fsMedianUncorrupt = calibration_parts.compute_kappa_bits_only_real(pipeline, smooth_fs, smooth_fsdq, options.expected_fs, options.fs_ok_var, status_out_smooth = 67108864, status_out_median = 134217728, starting_rate = options.compute_factors_sr, ending_rate = calibstatesr) + + # + # SRCQ BITS BRANCH + # + if not options.no_srcQ: + srcQSmoothInRange, srcQMedianUncorrupt = calibration_parts.compute_kappa_bits_only_real(pipeline, smooth_srcQ_inv, smooth_srcQdq, options.expected_srcQ, options.srcQ_ok_var, status_out_smooth = 268435456, status_out_median = 536870912, starting_rate = options.compute_factors_sr, ending_rate = calibstatesr) + + # + # 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 = 8388608, 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 = 8388608, status_out = 8388608) + 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 = 2097152, 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 = 2097152, status_out = 2097152) + 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 = 4194304, 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 = 4194304, status_out = 4194304) + 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, short_queue], [sus_coh_ok, long_queue], [darm_coh_ok, long_queue])) + 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 = 16777216, 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 = 16777216, status_out = 16777216) + 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, short_queue], [pcaly_line2_coh_ok, long_queue])) + + # + # H(T)-OK BIT BRANCH + # + + # First combine higher order bits to determine h(t)-OK + higherbits = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [filtersok, long_queue], [htproduced, short_queue], [obsreadytee, long_queue], [noinvalidinput, long_queue])) + htok_threshold = 28+33554432 + if options.apply_kappatst: + higherbits = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [higherbits, short_queue], [ktstSmoothInRange, long_queue])) + htok_threshold += 2048 + if options.apply_kappapu: + higherbits = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [higherbits, short_queue], [kpuSmoothInRange, long_queue])) + htok_threshold += 8192 + if options.apply_kappac: + higherbits = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [higherbits, short_queue], [kcSmoothInRange, long_queue])) + htok_threshold += 131072 + higherbitstee = pipeparts.mktee(pipeline, higherbits) + + # Now calculate h(t)-OK bit + htok = pipeparts.mkbitvectorgen(pipeline, higherbitstee, bit_vector = 1, threshold = htok_threshold) + htok = pipeparts.mkcapsfilter(pipeline, htok, calibstate_caps) + + # + # HW INJECTION BITS + # + + hwinjcbc = pipeparts.mkgeneric(pipeline, odcstatevectortee, "lal_logicalundersample", required_on = int(options.hw_inj_cbc_bitmask), status_out = 64) + hwinjcbc = pipeparts.mkcapsfilter(pipeline, hwinjcbc, calibstate_caps) + + hwinjburst = pipeparts.mkgeneric(pipeline, odcstatevectortee, "lal_logicalundersample", required_on = int(options.hw_inj_burst_bitmask), status_out = 128) + hwinjburst = pipeparts.mkcapsfilter(pipeline, hwinjburst, calibstate_caps) + + hwinjdetchar = pipeparts.mkgeneric(pipeline, odcstatevectortee, "lal_logicalundersample", required_on = int(options.hw_inj_detchar_bitmask), status_out = 256) + hwinjdetchar = pipeparts.mkcapsfilter(pipeline, hwinjdetchar, calibstate_caps) + + hwinjstoch = pipeparts.mkgeneric(pipeline, odcstatevectortee, "lal_logicalundersample", required_on = int(options.hw_inj_stoch_bitmask), status_out = 32) + hwinjstoch = pipeparts.mkcapsfilter(pipeline, hwinjstoch, calibstate_caps) + + + # + # COMBINE ALL BITS TO MAKE GDS-CALIB_STATE_VECTOR + # + + calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [nogap, long_queue], [higherbitstee, short_queue], [obsintenttee, long_queue], [htok, long_queue], [hwinjcbc, long_queue], [hwinjburst, long_queue], [hwinjdetchar, long_queue], [hwinjstoch, long_queue])) + if not options.no_kappatst or not options.no_kappapu or not options.no_kappac or not options.no_fcc: + calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [smoothingok, long_queue])) + if not options.no_coherence: + calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [coherence_bits, long_queue])) + if not options.no_kappatst: + if not options.apply_kappatst: + calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [ktstSmoothInRange, long_queue], [ktstMedianUncorrupt, long_queue])) + else: + calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [ktstMedianUncorrupt, long_queue])) + if not options.no_kappapu: + if not options.apply_kappapu: + calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [kpuSmoothInRange, long_queue], [kpuMedianUncorrupt, long_queue])) + else: + calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [kpuMedianUncorrupt, long_queue])) + if not options.no_kappac: + if not options.apply_kappac: + calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [kcSmoothInRange, long_queue], [kcMedianUncorrupt, long_queue])) + else: + calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [kcMedianUncorrupt, long_queue])) + if not options.no_fcc: + calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [fccSmoothInRange, long_queue], [fccMedianUncorrupt, long_queue])) + if not options.no_fs: + calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [fsSmoothInRange, long_queue], [fsMedianUncorrupt, long_queue])) + if not options.no_srcQ: + calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [srcQSmoothInRange, long_queue], [srcQMedianUncorrupt, long_queue])) + + calibstatevector = pipeparts.mkprogressreport(pipeline, calibstatevector, "progress_calibstatevec_%s" % instrument) + dqtagstr = "channel-name=%s:GDS-CALIB_STATE_VECTOR, instrument=%s" % (instrument, instrument) + calibstatevector = pipeparts.mktaginject(pipeline, calibstatevector, dqtagstr) + +# Resample the \kappa_a channels at the specified recording sample rate and change them to single precision channels +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, calibstatevector, short_queue).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_STATE_VECTOR%s" % (instrument, chan_prefix, chan_suffix))) + calibration_parts.mkqueue(pipeline, odcstatevectortee, long_queue).get_static_pad("src").link(mux.get_request_pad("%s:%s" % (instrument, options.dq_channel_name))) + strain_queue_length = long_queue +else: + strain_queue_length = short_queue + +# Link the strain branch to the muxer +calibration_parts.mkqueue(pipeline, strain, strain_queue_length).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_STRAIN%s" % (instrument, chan_prefix, chan_suffix))) +if options.remove_callines: + calibration_parts.mkqueue(pipeline, cleaned_strain, strain_queue_length).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_STRAIN_CLEAN%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, long_queue).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, long_queue).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, long_queue).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, long_queue).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, long_queue).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, long_queue).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, long_queue).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, long_queue).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, long_queue).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, long_queue).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, long_queue).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, long_queue).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, long_queue).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, long_queue).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, long_queue).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, long_queue).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 + either filter latency or kappa settling (if computing kappas), whichever is bigger +if not options.no_kappatst or not options.no_kappapu or not options.no_kappac or not options.no_fcc: + output_start = start + max(int(filter_settle_time), options.demodulation_filter_time + options.median_smoothing_time + options.factors_averaging_time) +else: + output_start = start + int(filter_settle_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 657faccdbf..ba28c742c4 100644 --- a/gstlal-calibration/bin/gstlal_compute_strain +++ b/gstlal-calibration/bin/gstlal_compute_strain @@ -1,6 +1,6 @@ #!/usr/bin/env python # -# Copyright (C) 2010-2015 Jordi Burguet-Castell, Madeline Wade +# 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 @@ -22,21 +22,21 @@ This pipeline produces h(t) given DARM_ERR and DARM_CTRL or given DELTAL_RESIDUA The differential arm length resulting from external sources is -\Delta L_{ext} = d_{err}/(\kappa_c C) + (A_tst * \kappa_tst + A_usum * \kappa_pu) d_{ctrl} +\Delta L_{ext} = d_{err}/(\kappa_c C) + (A_tst * \kappa_tst + A_pu * \kappa_pu) d_{ctrl} -where C is the sensing function, A_tst is the TST acutuation function, A_usum is the PUM+UIM+TOP 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 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, h(t) = \Delta L_{ext} / L . -The time-dependent gains (\kappa's) as well as the value for the coupled cavity pole (f_cc), the time-dependent gain of the PUM actuation (\kappa_pu) and the overall time-depenent gain of the actuation (\kappa_a) 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 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}. 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 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. -h(t) = (\Delta L_{res} * (1 / \kappa_c) * corrections + (\Delta L_{ctrl, TST} * \kappa_tst + \Delta L_{ctrl, USUM} * \kappa_pu) * corrections) / 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 Note: The \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. +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. @@ -44,7 +44,6 @@ Type gstlal_compute_strain --help to see the full list of command line options. """ import sys -import os import numpy import time import resource @@ -73,6 +72,10 @@ 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) @@ -91,6 +94,10 @@ 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) @@ -165,6 +172,7 @@ parser.add_option("--factors-averaging-time", metavar = "Sec", type = int, defau parser.add_option("--apply-kappapu", action = "store_true", help = "Set this to have the \kappa_pu factors multiply the actuation chain.") parser.add_option("--apply-kappatst", action = "store_true", help = "Set this to have the \kappa_tst factors multiply the actuation chain.") 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)") @@ -191,27 +199,6 @@ parser.add_option("--tst-exc-channel-name", metavar = "name", default = "SUS-ETM 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("--dewhitening", action = "store_true", help = "Dewhitening should be used on the relevant channels, since the incoming channels are whitened and single precision.") parser.add_option("--low-latency", action = "store_true", help = "Run the pipeline in low-latency mode. This uses minimal queueing. Otherwise, maximal queueing is used to prevent the pipeline from locking up.") -parser.add_option("--remove-callines", action = "store_true", help = "Remove calibration lines at known freqencies from h(t) using software.") -parser.add_option("--remove-powerlines", action = "store_true", help = "Remove 60 Hz spectral lines and some harmonics caused by power lines using witness channel PEM-EY_MAINSMON_EBAY_1_DQ.") -parser.add_option("--powerlines-channel-name", metavar = "name", default = "PEM-EY_MAINSMON_EBAY_1_DQ", help = "Set the name of the channel used as input for 60 Hz power lines to be removed. (Default = PEM-EY_MAINSMON_EBAY_1_DQ)") -parser.add_option("--remove-jitter-imc", action = "store_true", help = "Remove laser beam jitter using 4 IMC channels. This can significantly reduce noise in the spectrum.") -parser.add_option("--imc-a-pitch-channel-name", metavar = "name", default = "IMC-WFS_A_DC_PIT_OUT_DQ", help = "Set the name of one of the channels used as input from the IMC for removal of beam jitter noise. (Default = IMC-WFS_A_DC_PIT_OUT_DQ)") -parser.add_option("--imc-b-pitch-channel-name", metavar = "name", default = "IMC-WFS_B_DC_PIT_OUT_DQ", help = "Set the name of one of the channels used as input from the IMC for removal of beam jitter noise. (Default = IMC-WFS_B_DC_PIT_OUT_DQ)") -parser.add_option("--imc-a-yaw-channel-name", metavar = "name", default = "IMC-WFS_A_DC_YAW_OUT_DQ", help = "Set the name of one of the channels used as input from the IMC for removal of beam jitter noise. (Default = IMC-WFS_A_DC_YAW_OUT_DQ)") -parser.add_option("--imc-b-yaw-channel-name", metavar = "name", default = "IMC-WFS_B_DC_YAW_OUT_DQ", help = "Set the name of one of the channels used as input from the IMC for removal of beam jitter noise. (Default = IMC-WFS_B_DC_YAW_OUT_DQ)") -parser.add_option("--remove-jitter-psl", action = "store_true", help = "Remove laser beam jitter using the bullseye photodiode with 3 PSL channels. This can significantly reduce noise in the spectrum.") -parser.add_option("--bullseye-width-channel-name", metavar = "name", default = "PSL-DIAG_BULLSEYE_WID_OUT_DQ", help = "Set the name of one of the channels used as input from the bullseye photodiode for removal of beam jitter noise. (Default = PSL-DIAG_BULLSEYE_WID_OUT_DQ)") -parser.add_option("--bullseye-pitch-channel-name", metavar = "name", default = "PSL-DIAG_BULLSEYE_PIT_OUT_DQ", help = "Set the name of one of the channels used as input from the bullseye photodiode for removal of beam jitter noise. (Default = PSL-DIAG_BULLSEYE_PIT_OUT_DQ)") -parser.add_option("--bullseye-yaw-channel-name", metavar = "name", default = "PSL-DIAG_BULLSEYE_YAW_OUT_DQ", help = "Set the name of one of the channels used as input from the bullseye photodiode for removal of beam jitter noise. (Default = PSL-DIAG_BULLSEYE_YAW_OUT_DQ)") -parser.add_option("--remove-angular-control", action = "store_true", help = "Remove noise caused by angular control. Uses 4 ASC channels.") -parser.add_option("--asc-dhard-pitch-channel-name", metavar = "name", default = "ASC-DHARD_P_OUT_DQ", help = "Set the name of one of the channels used as input from the ASC to remove angular control noise. (Default = ASC-DHARD_P_OUT_DQ)") -parser.add_option("--asc-dhard-yaw-channel-name", metavar = "name", default = "ASC-DHARD_Y_OUT_DQ", help = "Set the name of one of the channels used as input from the ASC to remove angular control noise. (Default = ASC-DHARD_Y_OUT_DQ)") -parser.add_option("--asc-chard-pitch-channel-name", metavar = "name", default = "ASC-CHARD_P_OUT_DQ", help = "Set the name of one of the channels used as input from the ASC to remove angular control noise. (Default = ASC-CHARD_P_OUT_DQ)") -parser.add_option("--asc-chard-yaw-channel-name", metavar = "name", default = "ASC-CHARD_Y_OUT_DQ", help = "Set the name of one of the channels used as input from the ASC to remove angular control noise. (Default = ASC-CHARD_Y_OUT_DQ)") -parser.add_option("--remove-length-control", action = "store_true", help = "Remove noise caused by length control. Uses 3 LSC channels.") -parser.add_option("--lsc-srcl-channel-name", metavar = "name", default = "LSC-SRCL_IN1_DQ", help = "Set the name of one of the channels used as input from the LSC to remove length control noise. (Default = LSC-SRCL_IN1_DQ)") -parser.add_option("--lsc-mich-channel-name", metavar = "name", default = "LSC-MICH_IN1_DQ", help = "Set the name of one of the channels used as input from the LSC to remove length control noise. (Default = LSC-MICH_IN1_DQ)") -parser.add_option("--lsc-prcl-channel-name", metavar = "name", default = "LSC-PRCL_IN1_DQ", help = "Set the name of one of the channels used as input from the LSC to remove length control noise. (Default = LSC-PRCL_IN1_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)") @@ -363,12 +350,6 @@ chan_prefix = options.chan_prefix # If td is true we will perform filtering in the time domain (direct convolution) in all FIR filtering routines below 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.ifo == "H1" and options.data_source == "frames" and int(options.gps_start_time) > 1175954418) or (options.ifo == "H1" and options.data_source == "lvshm" and now() > 1175954418) or (options.ifo == "L1" and options.data_source == "frames" and int(options.gps_start_time) > 1180184418) or (options.ifo == "L1" and options.data_source == "lvshm" and now() > 1180184418)): - remove_esd_act_line = True -elif not options.factors_from_filters_file: - remove_esd_act_line = False - # # Load in the filters file that contains filter coefficients, etc. # @@ -395,12 +376,6 @@ if options.factors_from_filters_file: EP8_imag = float(filters["EP8_imag"]) EP9_real = float(filters["EP9_real"]) EP9_imag = float(filters["EP9_imag"]) - try: - EP10_real = float(filters["EP10_real"]) - EP10_imag = float(filters["EP10_imag"]) - remove_esd_act_line = True - except: - remove_esd_act_line = False try: EP11_real = float(filters["EP11_real"]) EP11_imag = float(filters["EP11_imag"]) @@ -513,75 +488,6 @@ if options.full_calibration: resdewhitendelay = int(filters["dewhiten_err_delay"]) resdewhiten = filters["dewhiten_err"] -# If we're removing 60 Hz lines from the spectrum, load another filter -if options.remove_powerlines: - try: - powerlinessr = int(filters["powerlines_sr"]) - powerlinesdelay = int(filters["powerlines_delay"]) - powerlinesfilt = filters["powerlines_filt"] - except: - raise ValueError("Cannot remove 60 Hz lines because the filters file does contain the needed information") - -# If we're removing laser beam jitter noise from the spectrum, load several more filters -if options.remove_jitter_imc: - try: - imcapitsr = int(filters["jitter_imc_a_pit_sr"]) - imcayawsr = int(filters["jitter_imc_a_yaw_sr"]) - imcbpitsr = int(filters["jitter_imc_b_pit_sr"]) - imcbyawsr = int(filters["jitter_imc_b_yaw_sr"]) - imcapitdelay = int(filters["jitter_imc_a_pit_delay"]) - imcayawdelay = int(filters["jitter_imc_a_yaw_delay"]) - imcbpitdelay = int(filters["jitter_imc_b_pit_delay"]) - imcbyawdelay = int(filters["jitter_imc_b_yaw_delay"]) - imcapitfilt = filters["jitter_imc_a_pit_filt"] - imcayawfilt = filters["jitter_imc_a_yaw_filt"] - imcbpitfilt = filters["jitter_imc_b_pit_filt"] - imcbyawfilt = filters["jitter_imc_b_yaw_filt"] - except: - raise ValueError("Cannot remove beam jitter using imc inputs because the filters file does contain the needed information") -if options.remove_jitter_psl: - try: - bullseyewidsr = int(filters["jitter_bullseye_wid_sr"]) - bullseyepitsr = int(filters["jitter_bullseye_pit_sr"]) - bullseyeyawsr = int(filters["jitter_bullseye_yaw_sr"]) - bullseyewiddelay = int(filters["jitter_bullseye_wid_delay"]) - bullseyepitdelay = int(filters["jitter_bullseye_pit_delay"]) - bullseyeyawdelay = int(filters["jitter_bullseye_yaw_delay"]) - bullseyewidfilt = filters["jitter_bullseye_wid_filt"] - bullseyepitfilt = filters["jitter_bullseye_pit_filt"] - bullseyeyawfilt = filters["jitter_bullseye_yaw_filt"] - except: - raise ValueError("Cannot remove beam jitter using bullseye inputs because the filters file does contain the needed information") -if options.remove_angular_control: - try: - ascdpitsr = int(filters["asc_d_pit_sr"]) - ascdyawsr = int(filters["asc_d_yaw_sr"]) - asccpitsr = int(filters["asc_c_pit_sr"]) - asccyawsr = int(filters["asc_c_yaw_sr"]) - ascdpitdelay = int(filters["asc_d_pit_delay"]) - asccyawdelay = int(filters["asc_d_yaw_delay"]) - asccpitdelay = int(filters["asc_c_pit_delay"]) - ascdyawdelay = int(filters["asc_c_yaw_delay"]) - ascdpitfilt = filters["asc_d_pit_filt"] - ascdyawfilt = filters["asc_d_yaw_filt"] - asccpitfilt = filters["asc_c_pit_filt"] - asccyawfilt = filters["asc_c_yaw_filt"] - except: - raise ValueError("Cannot remove angular control noise using ASC inputs because the filters file does contain the needed information") -if options.remove_length_control: - try: - lscsrclsr = int(filters["lsc_srcl_sr"]) - lscmichsr = int(filters["lsc_mich_sr"]) - lscprclsr = int(filters["lsc_prcl_sr"]) - lscsrcldelay = int(filters["lsc_srcl_delay"]) - lscmichdelay = int(filters["lsc_mich_delay"]) - lscprcldelay = int(filters["lsc_prcl_delay"]) - lscsrclfilt = filters["lsc_srcl_filt"] - lscmichfilt = filters["lsc_mich_filt"] - lscprclfilt = filters["lsc_prcl_filt"] - except: - raise ValueError("Cannot remove length control noise using LSC inputs because the filters file does contain the needed information") - # Set up queue parameters if options.low_latency: @@ -649,11 +555,6 @@ if not options.factors_from_filters_file: 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.remove_callines and remove_esd_act_line: - 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.no_fs or not options.no_srcQ: 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))) @@ -693,24 +594,6 @@ elif options.partial_calibration: channel_list.extend(((instrument, options.deltal_res_channel_name), (instrument, options.deltal_tst_channel_name), (instrument, options.deltal_pum_channel_name), (instrument, options.deltal_uim_channel_name))) headkeys.extend(("res", "tst", "pum", "uim")) -# If we are removing additional noise from the spectrum (beam jitter, angular control, 60 Hz lines, etc.), we need more channels -if options.remove_powerlines: - channel_list.append((instrument, options.powerlines_channel_name)) - headkeys.append("powerlines") -if options.remove_jitter_imc: - channel_list.extend(((instrument, options.imc_a_pitch_channel_name), (instrument, options.imc_a_yaw_channel_name), (instrument, options.imc_b_pitch_channel_name), (instrument, options.imc_b_yaw_channel_name))) - headkeys.extend(("imc_a_pitch", "imc_a_yaw", "imc_b_pitch", "imc_b_yaw")) -if options.remove_jitter_psl: - channel_list.extend(((instrument, options.bullseye_width_channel_name), (instrument, options.bullseye_pitch_channel_name), (instrument, options.bullseye_yaw_channel_name))) - headkeys.extend(("bullseye_width", "bullseye_pitch", "bullseye_yaw")) -if options.remove_angular_control: - channel_list.extend(((instrument, options.asc_dhard_pitch_channel_name), (instrument, options.asc_dhard_yaw_channel_name), (instrument, options.asc_chard_pitch_channel_name), (instrument, options.asc_chard_yaw_channel_name))) - headkeys.extend(("asc_dhard_pitch", "asc_dhard_yaw", "asc_chard_pitch", "asc_chard_yaw")) -if options.remove_length_control: - channel_list.extend(((instrument, options.lsc_srcl_channel_name), (instrument, options.lsc_mich_channel_name), (instrument, options.lsc_prcl_channel_name))) - headkeys.extend(("lsc_srcl", "lsc_mich", "lsc_prcl")) - - #################################################################################################### ####################################### Main Pipeline ############################################## #################################################################################################### @@ -799,7 +682,7 @@ if not options.no_coherence: if highfreq_coh_use >= 2: pcaly_line2_coh = pipeparts.mktee(pipeline, pcaly_line2_coh) -# Set up computations for \kappa_a, \kappa_tst,\kappa_c, \kappa_pu, f_cc, if applicable +# Set up computations for \kappa_tst,\kappa_c, \kappa_pu, f_cc, if applicable if not options.no_kappac or not options.no_fcc or not options.no_kappatst or not options.no_kappapu or not options.no_srcQ or not options.no_fs: # pcal excitation channel, which will be demodulated @@ -828,8 +711,6 @@ if not options.no_kappac or not options.no_fcc or not options.no_kappatst or not # demodulate the TST excitation channel at the ESD actuation line frequency tstexc_at_esd_act_freq = calibration_parts.demodulate(pipeline, tstexc, esd_act_line_freq, td, compute_calib_factors_complex_caps, integration_samples) - if options.remove_callines and remove_esd_act_line: - tstexc_at_esd_act_freq = pipeparts.mktee(pipeline, tstexc_at_esd_act_freq) # demodulate DARM_ERR at the ESD actuation line frequency derr_at_esd_act_freq = calibration_parts.demodulate(pipeline, derrtee, esd_act_line_freq, td, compute_calib_factors_complex_caps, integration_samples) @@ -936,8 +817,6 @@ if not options.no_kappac or not options.no_fcc or not options.no_kappapu or not if not options.no_kappac or not options.no_fcc or not options.no_srcQ or not options.no_fs: # demodulate 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, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_opt_gain_fcc_freq_real, pcal_corr_at_opt_gain_fcc_freq_imag) - if options.remove_callines: - pcal_at_opt_gain_freq = pipeparts.mktee(pipeline, pcal_at_opt_gain_freq) # 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, td, compute_calib_factors_complex_caps, integration_samples) @@ -1017,6 +896,7 @@ if not options.no_kappac or not options.no_fcc or not options.no_kappapu or not smooth_fcctee = pipeparts.mktee(pipeline, smooth_fcc) + # compute f_s and Q if not options.no_fs or not options.no_srcQ: # demodulate PCAL channel and apply the PCAL correction factor at SRC detuning line frequency pcal_at_src_freq = calibration_parts.demodulate(pipeline, pcaltee, src_pcal_line_freq, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_src_freq_real, pcal_corr_at_src_freq_imag) @@ -1104,178 +984,6 @@ if not options.no_kappac or not options.no_fcc or not options.no_kappapu or not if not options.no_dq_vector: smooth_srcQ_inv = pipeparts.mktee(pipeline, smooth_srcQ_inv) -# -# Calibration Line Removal -# - -if options.remove_callines: - # if we didn't compute the kappas, we still need to get the pcal channel - if options.no_kappatst and options.no_kappapu and options.no_kappac and options.no_fcc and options.no_srcQ and options.no_fs: - pcal = calibration_parts.caps_and_progress(pipeline, head_dict["pcal"], hoft_caps, "pcal") - pcaltee = pipeparts.mktee(pipeline, pcal) - # Demodulate pcal at the ~30 Hz pcal line - pcal_at_darm_act_freq = calibration_parts.demodulate(pipeline, pcal_tee, darm_act_line_freq, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_darm_act_freq_real, pcal_corr_at_darm_act_freq_imag) - - # Reconstruct a calibrated (negative) pcal at only the ~30 Hz pcal line - pcaly_line1 = calibration_parts.mkresample(pipeline, pcal_at_darm_act_freq, 3, False, "audio/x-raw, format=Z128LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) - pcaly_line1 = pipeparts.mkgeneric(pipeline, pcaly_line1, "lal_demodulate", line_frequency = -1.0 * darm_act_line_freq, prefactor_real = 2.0) - remove_pcaly_line1, trash = calibration_parts.split_into_real(pipeline, pcaly_line1) - pipeparts.mkfakesink(pipeline, trash) - - # 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, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_opt_gain_fcc_freq_real, pcal_corr_at_opt_gain_fcc_freq_imag) - # Reconstruct a calibrated (negative) pcal at only the ~300 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, trash = calibration_parts.split_into_real(pipeline, pcaly_line2) - pipeparts.mkfakesink(pipeline, trash) - - # Add the first two components together. We will add this to h(t) to remove these lines - remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_pcaly_line1, long_queue], [remove_pcaly_line2, short_queue])) - - if remove_esd_act_line: - # Make sure we have demodulated the ESD excitation channel at the ~30 Hz ESD line - if options.no_kappac and options.no_fcc and options.no_kappatst and options.no_kappapu and options.no_srcQ and options.no_fs: - tstexc_at_esd_act_freq = calibration_parts.demodulate(pipeline, tstexc, esd_act_line_freq, td, compute_calib_factors_complex_caps, integration_samples) - if options.factors_from_filters_file: - esd_act_line = calibration_parts.complex_audioamplify(pipeline, tstexc_at_esd_act_freq, EP10_real, EP10_imag) - else: - # EP10 was read from the frames - EP10 = calibration_parts.merge_into_complex(pipeline, head_dict["EP10_real"], head_dict["EP10_imag"], long_queue, short_queue) - esd_act_line = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, [tstexc_at_esd_act_freq, long_queue], [EP10, short_queue])) - # Reconstruct a calibrated (negative) ESD injection at the ~30 Hz ESD line - if options.apply_kappatst: - # Multiply by the real part of kappa_tst - esd_act_line = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, [esd_act_line, long_queue], [pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, smooth_ktstRtee, matrix=[[1.0, 0.0]])), short_queue])) - 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, trash = calibration_parts.split_into_real(pipeline, esd_act_line_remove) - pipeparts.mkfakesink(pipeline, trash) - # Add into the total line removal stream - remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_from_strain, short_queue], [esd_act_line_remove, long_queue])) - - if remove_high_pcal_line: - # Demodulate pcal at the ~1kHz pcal line - pcaly_line3 = calibration_parts.demodulate(pipeline, pcaltee, high_pcal_line_freq, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_high_line_freq_real, pcal_corr_at_high_line_freq_imag) - # Reconstruct a calibrated (negative) pcal at only the ~1kHz pcal line - pcaly_line3 = calibration_parts.mkresample(pipeline, pcaly_line3, 3, False, "audio/x-raw, format=Z128LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) - pcaly_line3 = pipeparts.mkgeneric(pipeline, pcaly_line3, "lal_demodulate", line_frequency = -1.0 * high_pcal_line_freq, prefactor_real = 2.0) - remove_pcaly_line3, trash = calibration_parts.split_into_real(pipeline, pcaly_line3) - pipeparts.mkfakesink(pipeline, trash) - # Add into the total line removal stream - remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_from_strain, long_queue], [remove_pcaly_line3, short_queue])) - - if remove_roaming_pcal_line: - # Demodulate pcal at the ~3kHz pcal line - pcaly_line4 = calibration_parts.demodulate(pipeline, pcaltee, roaming_pcal_line_freq, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_roaming_line_real, pcal_corr_at_roaming_line_imag) - # Reconstruct a calibrated (negative) pcal at only the ~3kHz pcal line - pcaly_line4 = calibration_parts.mkresample(pipeline, pcaly_line4, 3, False, "audio/x-raw, format=Z128LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) - pcaly_line4 = pipeparts.mkgeneric(pipeline, pcaly_line4, "lal_demodulate", line_frequency = -1.0 * roaming_pcal_line_freq, prefactor_real = 2.0) - remove_pcaly_line4, trash = calibration_parts.split_into_real(pipeline, pcaly_line4) - pipeparts.mkfakesink(pipeline, trash) - # Add into the total line removal stream - remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_from_strain, long_queue], [remove_pcaly_line4, short_queue])) - - if remove_src_pcal_line: - # Make sure we have demodulated pcal at the ~8 Hz pcal line - if options.no_fs and options.no_srcQ: - pcal_at_src_freq = calibration_parts.demodulate(pipeline, pcaltee, src_pcal_line_freq, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_src_freq_real, pcal_corr_at_src_freq_imag) - # Reconstruct a calibrated (negative) pcal at only the ~3kHz pcal line - pcaly_line0 = calibration_parts.mkresample(pipeline, pcal_at_src_freq, 3, False, "audio/x-raw, format=Z128LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) - pcaly_line0 = pipeparts.mkgeneric(pipeline, pcaly_line0, "lal_demodulate", line_frequency = -1.0 * src_pcal_line_freq, prefactor_real = 2.0) - remove_pcaly_line0, trash = calibration_parts.split_into_real(pipeline, pcaly_line0) - pipeparts.mkfakesink(pipeline, trash) - # Add into the total line removal stream - remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_from_strain, long_queue], [remove_pcaly_line0, short_queue])) - -if options.remove_powerlines: - powerlines = calibration_parts.caps_and_progress(pipeline, head_dict["powerlines"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % powerlinessr, "powerlines") - powerlines = pipeparts.mkfirbank(pipeline, powerlines, latency = int(powerlinesdelay), fir_matrix = [powerlinesfilt[::-1]], time_domain = td) - powerlines = calibration_parts.mkresample(pipeline, powerlines, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) - if options.remove_callines: - remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_from_strain, long_queue], [powerlines, short_queue])) - else: - remove_from_strain = powerlines - -if options.remove_jitter_imc: - imc_a_pitch = calibration_parts.caps_and_progress(pipeline, head_dict["imc_a_pitch"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % imcapitsr, "imc_a_pitch") - imc_a_yaw = calibration_parts.caps_and_progress(pipeline, head_dict["imc_a_yaw"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % imcayawsr, "imc_a_yaw") - imc_b_pitch = calibration_parts.caps_and_progress(pipeline, head_dict["imc_b_pitch"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % imcbpitsr, "imc_b_pitch") - imc_b_yaw = calibration_parts.caps_and_progress(pipeline, head_dict["imc_b_yaw"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % imcbyawsr, "imc_b_yaw") - - imc_a_pitch = pipeparts.mkfirbank(pipeline, imc_a_pitch, latency = int(imcapitdelay), fir_matrix = [imcapitfilt[::-1]], time_domain = td) - imc_a_yaw = pipeparts.mkfirbank(pipeline, imc_a_yaw, latency = int(imcayawdelay), fir_matrix = [imcayawfilt[::-1]], time_domain = td) - imc_b_pitch = pipeparts.mkfirbank(pipeline, imc_b_pitch, latency = int(imcbpitdelay), fir_matrix = [imcbpitfilt[::-1]], time_domain = td) - imc_b_yaw = pipeparts.mkfirbank(pipeline, imc_b_yaw, latency = int(imcbyawdelay), fir_matrix = [imcbyawfilt[::-1]], time_domain = td) - - imc_a_pitch = calibration_parts.mkresample(pipeline, imc_a_pitch, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) - imc_a_yaw = calibration_parts.mkresample(pipeline, imc_a_yaw, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) - imc_b_pitch = calibration_parts.mkresample(pipeline, imc_b_pitch, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) - imc_b_yaw = calibration_parts.mkresample(pipeline, imc_b_yaw, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) - - if options.remove_callines or options.remove_powerlines: - remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_from_strain, long_queue], [imc_a_pitch, long_queue], [imc_a_yaw, long_queue], [imc_b_pitch, long_queue], [imc_b_yaw, long_queue])) - else: - remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [imc_a_pitch, long_queue], [imc_a_yaw, long_queue], [imc_b_pitch, long_queue], [imc_b_yaw, long_queue])) - -if options.remove_jitter_psl: - bullseye_width = calibration_parts.caps_and_progress(pipeline, head_dict["bullseye_width"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % bullseyewidsr, "bullseye_width") - bullseye_pitch = calibration_parts.caps_and_progress(pipeline, head_dict["bullseye_pitch"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % bullseyepitsr, "bullseye_pitch") - bullseye_yaw = calibration_parts.caps_and_progress(pipeline, head_dict["bullseye_yaw"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % bullseyeyawsr, "bullseye_yaw") - - bullseye_width = pipeparts.mkfirbank(pipeline, bullseye_width, latency = int(bullseyewiddelay), fir_matrix = [bullseyewidfilt[::-1]], time_domain = td) - bullseye_pitch = pipeparts.mkfirbank(pipeline, bullseye_pitch, latency = int(bullseyepitdelay), fir_matrix = [bullseyepitfilt[::-1]], time_domain = td) - bullseye_yaw = pipeparts.mkfirbank(pipeline, bullseye_yaw, latency = int(bullseyeyawdelay), fir_matrix = [bullseyeyawfilt[::-1]], time_domain = td) - - bullseye_width = calibration_parts.mkresample(pipeline, bullseye_width, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) - bullseye_pitch = calibration_parts.mkresample(pipeline, bullseye_pitch, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) - bullseye_yaw = calibration_parts.mkresample(pipeline, bullseye_yaw, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) - - if options.remove_callines or options.remove_powerlines or options.remove_jitter_imc: - remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_from_strain, long_queue], [bullseye_width, long_queue], [bullseye_pitch, long_queue], [bullseye_yaw, long_queue])) - else: - remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [bullseye_width, long_queue], [bullseye_pitch, long_queue], [bullseye_yaw, long_queue])) - -if options.remove_angular_control: - asc_dhard_pitch = calibration_parts.caps_and_progress(pipeline, head_dict["asc_dhard_pitch"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % ascdpitsr, "asc_dhard_pitch") - asc_dhard_yaw = calibration_parts.caps_and_progress(pipeline, head_dict["asc_dhard_yaw"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % ascdyawsr, "asc_dhard_yaw") - asc_chard_pitch = calibration_parts.caps_and_progress(pipeline, head_dict["asc_chard_pitch"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % asccpitsr, "asc_chard_pitch") - asc_chard_yaw = calibration_parts.caps_and_progress(pipeline, head_dict["asc_chard_yaw"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % asccyawsr, "asc_chard_yaw") - - asc_dhard_pitch = pipeparts.mkfirbank(pipeline, asc_dhard_pitch, latency = int(ascdpitdelay), fir_matrix = [ascdpitfilt[::-1]], time_domain = td) - asc_dhard_yaw = pipeparts.mkfirbank(pipeline, asc_dhard_yaw, latency = int(ascdyawdelay), fir_matrix = [ascdyawfilt[::-1]], time_domain = td) - asc_chard_pitch = pipeparts.mkfirbank(pipeline, asc_chard_pitch, latency = int(asccpitdelay), fir_matrix = [asccpitfilt[::-1]], time_domain = td) - asc_chard_yaw = pipeparts.mkfirbank(pipeline, asc_chard_yaw, latency = int(asccyawdelay), fir_matrix = [asccyawfilt[::-1]], time_domain = td) - - asc_dhard_pitch = calibration_parts.mkresample(pipeline, asc_dhard_pitch, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) - asc_dhard_yaw = calibration_parts.mkresample(pipeline, asc_dhard_yaw, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) - asc_chard_pitch = calibration_parts.mkresample(pipeline, asc_chard_pitch, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) - asc_chard_yaw = calibration_parts.mkresample(pipeline, asc_chard_yaw, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) - - if options.remove_callines or options.remove_powerlines or options.remove_jitter_imc or options.remove_jitter_psl: - remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_from_strain, long_queue], [asc_dhard_pitch, long_queue], [asc_dhard_yaw, long_queue], [asc_chard_pitch, long_queue], [asc_chard_yaw, long_queue])) - else: - remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [asc_dhard_pitch, long_queue], [asc_dhard_yaw, long_queue], [asc_chard_pitch, long_queue], [asc_chard_yaw, long_queue])) - -if options.remove_length_control: - lsc_srcl = calibration_parts.caps_and_progress(pipeline, head_dict["lsc_srcl"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % lscsrclsr, "lsc_srcl") - lsc_mich = calibration_parts.caps_and_progress(pipeline, head_dict["lsc_mich"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % lscmichsr, "lsc_mich") - lsc_prcl = calibration_parts.caps_and_progress(pipeline, head_dict["lsc_prcl"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % lscprclsr, "lsc_prcl") - - lsc_srcl = pipeparts.mkfirbank(pipeline, lsc_srcl, latency = int(lscsrcldelay), fir_matrix = [lscsrclfilt[::-1]], time_domain = td) - lsc_mich = pipeparts.mkfirbank(pipeline, lsc_mich, latency = int(lscmichdelay), fir_matrix = [lscmichfilt[::-1]], time_domain = td) - lsc_prcl = pipeparts.mkfirbank(pipeline, lsc_prcl, latency = int(lscprcldelay), fir_matrix = [lscprclfilt[::-1]], time_domain = td) - - lsc_srcl = calibration_parts.mkresample(pipeline, lsc_srcl, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) - lsc_mich = calibration_parts.mkresample(pipeline, lsc_mich, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) - lsc_prcl = calibration_parts.mkresample(pipeline, lsc_prcl, 3, False, "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) - - if options.remove_callines or options.remove_powerlines or options.remove_jitter_imc or options.remove_jitter_psl or options.remove_angular_control: - remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_from_strain, long_queue], [lsc_srcl, long_queue], [lsc_mich, long_queue], [lsc_prcl, long_queue])) - else: - remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [lsc_srcl, long_queue], [lsc_mich, long_queue], [lsc_prcl, long_queue])) - # # CONTROL BRANCH # @@ -1418,11 +1126,6 @@ filter_latency = max(res_filter_latency, tst_filter_latency, pumuim_filter_laten # Add control and residual chains and divide by L to make h(t) strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [res, res_queue_length], [ctrl, ctrl_queue_length])) -# Remove the calibration lines, if we want -if options.remove_callines or options.remove_powerlines or options.remove_jitter_imc or options.remove_jitter_psl or options.remove_angular_control or options.remove_length_control: - remove_from_strain = pipeparts.mkaudioamplify(pipeline, remove_from_strain, -1.0) - strain = pipeparts.mktee(pipeline, strain) - cleaned_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [strain, short_queue], [remove_from_strain, long_queue])) # Divide by L in a way that is compatitble with old and new filters files, since old filter files don't recored "arm length" try: strain = pipeparts.mkaudioamplify(pipeline, strain, 1.0/float(filters["arm_length"])) @@ -1431,13 +1134,6 @@ except KeyError: strain = pipeparts.mkprogressreport(pipeline, strain, "progress_hoft_%s" % instrument) -if options.remove_callines: - try: - cleaned_strain = pipeparts.mkaudioamplify(pipeline, cleaned_strain, 1.0/float(filters["arm_length"])) - except KeyError: - cleaned_strain = pipeparts.mkaudioamplify(pipeline, cleaned_strain, 1.0/3994.5) - - cleaned_strain = pipeparts.mkprogressreport(pipeline, cleaned_strain, "progress_hoft_cleaned_%s" % instrument) # Put the units back to strain before writing to frames straintagstr = "units=strain,channel-name=%sCALIB_STRAIN%s,instrument=%s" % (chan_prefix, chan_suffix, instrument) cleaned_straintagstr = "units=strain,channel-name=%sCALIB_STRAIN_CLEAN%s,instrument=%s" % (chan_prefix, chan_suffix, instrument) @@ -1446,8 +1142,6 @@ if not options.no_dq_vector: strain = pipeparts.mktaginject(pipeline, straintee, straintagstr) else: strain = pipeparts.mktaginject(pipeline, strain, straintagstr) -if options.remove_callines: - cleaned_strain = pipeparts.mktaginject(pipeline, cleaned_strain, cleaned_straintagstr) # # CALIB_STATE_VECTOR BRANCH @@ -1694,7 +1388,6 @@ if not options.no_dq_vector: dqtagstr = "channel-name=%s:GDS-CALIB_STATE_VECTOR, instrument=%s" % (instrument, instrument) calibstatevector = pipeparts.mktaginject(pipeline, calibstatevector, dqtagstr) -# Resample the \kappa_a channels at the specified recording sample rate and change them to single precision channels 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 @@ -1798,8 +1491,6 @@ else: # Link the strain branch to the muxer calibration_parts.mkqueue(pipeline, strain, strain_queue_length).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_STRAIN%s" % (instrument, chan_prefix, chan_suffix))) -if options.remove_callines: - calibration_parts.mkqueue(pipeline, cleaned_strain, strain_queue_length).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_STRAIN_CLEAN%s" % (instrument, chan_prefix, chan_suffix))) # Link the real and imaginary parts of \kappa_tst to the muxer if not options.no_kappatst: -- GitLab