Skip to content

Commit

Permalink
Merge 51115f6 into fd8de27
Browse files Browse the repository at this point in the history
  • Loading branch information
Alex L. Urban committed Apr 11, 2019
2 parents fd8de27 + 51115f6 commit 03b6dba
Show file tree
Hide file tree
Showing 10 changed files with 502 additions and 23 deletions.
9 changes: 8 additions & 1 deletion bin/gwdetchar-scattering
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,10 @@
# along with GW DetChar. If not, see <http://www.gnu.org/licenses/>.

"""Search for evidence of beam scattering based on optic velocity
This tool identifies time segments when evidence for scattering is strong,
to compare predicted fringes against spectrogram measurements for a specific
time please use the command-line module: `python gwdetchar.scattering --help`
"""

from __future__ import division
Expand Down Expand Up @@ -559,7 +563,10 @@ if nscans > 0:
# render HTML
page.h2('Omega Scans')
page.p('The following event times correspond to significant Omicron '
'triggers that occurr during the scattering segments found above:')
'triggers that occurr during the scattering segments found above. '
'To compare these against fringe frequency projections, please '
'use the scattering module:')
page.pre('$ python -m gwdetchar.scattering --help')
page.add(htmlio.scaffold_omega_scans(
omegatimes, args.main_channel, scandir=scandir))
elif args.omega_scans:
Expand Down
2 changes: 1 addition & 1 deletion gwdetchar/lasso/tests/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
# You should have received a copy of the GNU General Public License
# along with GW DetChar. If not, see <http://www.gnu.org/licenses/>.

"""Tests for :mod:`gwdetchar.omega`
"""Tests for :mod:`gwdetchar.lasso`
"""

__author__ = 'Alex Urban <alexander.urban@ligo.org>'
30 changes: 30 additions & 0 deletions gwdetchar/scattering/__init__.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,30 @@
# coding=utf-8
# Copyright (C) Alex Urban (2019)
#
# This file is part of the GW DetChar python package.
#
# GW DetChar 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 3 of the License, or
# (at your option) any later version.
#
# GW DetChar 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 GW DetChar. If not, see <http://www.gnu.org/licenses/>.

"""Methods and utilties for optical scattering
"""

__author__ = 'Alex Urban <alexander.urban@ligo.org>'
__credits__ = 'Duncan Macleod <duncan.macleod@ligo.org>' \
'Joshua Smith <joshua.smith@ligo.org>' \
'Andrew Lundgren <andrew.lundgren>@ligo.org>'

# -- imports ------------------------------------------------------------------

# import scattering utils
from .core import *
173 changes: 173 additions & 0 deletions gwdetchar/scattering/__main__.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,173 @@
# coding=utf-8
# Copyright (C) Alex Urban (2019)
#
# This file is part of the GW DetChar python package.
#
# GW DetChar 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 3 of the License, or
# (at your option) any later version.
#
# GW DetChar 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 GW DetChar. If not, see <http://www.gnu.org/licenses/>.

"""Simple command-line interface to gwdetchar.scattering
This module scans through records of optic motion and predicts scattering
fringe frequencies. For those channels with fringes above a user-specified
threshold, a plot is created comparing the fringes to a high-resolution Q-scan
spectrogram.
To identify time segments where scattering is likely, please use the
command-line script: `gwdetchar-scattering --help`
"""

import os
import sys

from matplotlib import use
use('agg') # noqa

from gwpy.time import to_gps

from .. import (cli, const)
from ..omega import highpass
from ..io.datafind import get_data

from . import *
from . import plot

__author__ = 'Alex Urban <alexander.urban@ligo.org>'
__credits__ = 'Joshua Smith <joshua.smith@ligo.org>' \
'Andrew Lundgren <andrew.lundgren>@ligo.org>' \
'Duncan Macleod <duncan.macleod@ligo.org>'

# global variables

ASD_KW = {
'method': 'median',
'fftlength': 8,
'overlap': 4,
}

MOTION_CHANNELS = [channel for key in OPTIC_MOTION_CHANNELS.keys()
for channel in OPTIC_MOTION_CHANNELS[key]]

logger = cli.logger('gwdetchar.scattering')


# -- main function ------------------------------------------------------------

def main(args=None):
"""Parse command-line arguments, process optics, and write plots
"""
# define command-line arguments
parser = cli.create_parser(description=__doc__)
parser.add_argument('gpstime', type=to_gps,
help='GPS time or datestring to analyze')
cli.add_ifo_option(parser, ifo=const.IFO)
parser.add_argument('-d', '--duration', type=float, default=60,
help='Duration (seconds) of analysis, default: 60')
parser.add_argument('-t', '--frequency-threshold', type=float, default=15,
help='critical fringe frequency threshold (Hz), '
'default: %(default)s')
parser.add_argument('-m', '--multipliers', default='1,2,4,8',
help='harmonic numbers to plot predicted motion for, '
'should be given as a comma-separated list of '
'numbers, default: %(default)s')
parser.add_argument('-x', '--multiplier-for-threshold', type=int,
default=4,
help='frequency multiplier to use when applying '
'--frequency-threshold, default: %(default)s')
parser.add_argument('-w', '--primary-channel',
default='GDS-CALIB_STRAIN',
help='name of primary channel (without IFO prefix), '
'default: %(default)s')
parser.add_argument('-W', '--primary-frametype', default='{IFO}_HOFT_C00',
help='frametype from which to read primary channel, '
'default: %(default)s')
parser.add_argument('-a', '--aux-frametype', default='{IFO}_R',
help='frametype from which to read aux channels, '
'default: %(default)s')
parser.add_argument('-o', '--output-dir', type=os.path.abspath,
default=os.curdir,
help='Output directory for analysis, '
'default: %(default)s')
parser.add_argument('-c', '--colormap', default='viridis',
help='name of colormap to use, default: %(default)s')

# parse arguments
args = parser.parse_args(args)
ifo = args.ifo
gps = float(args.gpstime)
gpsstart = gps - 0.5 * args.duration
gpsend = gps + 0.5 * args.duration
primary = ':'.join([ifo, args.primary_channel])
thresh = args.frequency_threshold
multipliers = [int(x) for x in args.multipliers.split(',')]
harmonic = args.multiplier_for_threshold
if '{IFO}' in args.primary_frametype:
args.primary_frametype = args.primary_frametype.format(IFO=ifo)
if '{IFO}' in args.aux_frametype:
args.aux_frametype = args.aux_frametype.format(IFO=ifo)
logger.info('{0} Scattering: {1}'.format(ifo, gps))

# set up spectrogram
logger.debug('Setting up a Q-scan spectrogram of {}'.format(primary))
hoft = get_data(primary, start=gps-34, end=gps+34,
frametype=args.primary_frametype,
verbose='Reading primary channel:'.rjust(30))
hoft = highpass(hoft, f_low=thresh).resample(256)
qspecgram = hoft.q_transform(qrange=(4, 150), frange=(0, 60), gps=gps,
fres=0.1, outseg=(gpsstart, gpsend), **ASD_KW)
qspecgram.name = primary

# process channels
channels = [':'.join([ifo, c]) for c in MOTION_CHANNELS]
data = get_data(
channels, start=gpsstart, end=gpsend, frametype=args.aux_frametype,
verbose='Reading auxiliary sensors:'.rjust(30))
count = 0 # running count of plots written
for channel in channels:
logger.info(' -- Processing {} -- '.format(channel))
try:
motion = data[channel].detrend().resample(128)
except KeyError:
logger.warning('Skipping {}'.format(channel))
pass
# predict scattering frequencies
fringe = get_fringe_frequency(motion, multiplier=1)
ind = fringe.argmax()
fmax = fringe.value[ind]
tmax = fringe.times.value[ind]
logger.debug('Maximum scatter frequency {0:.2f} Hz at GPS second '
'{1:.2f}'.format(fmax, tmax))
if harmonic * fmax < thresh:
logger.warning('No significant evidence of scattering '
'found in {}'.format(channel))
continue
# plot spectrogram and fringe frequency
output = os.path.join(
args.output_dir,
'%s-%s-%s-{}.png' % (
channel.replace(':', '-'), gps, args.duration)
)
logger.debug('Plotting spectra and predicted fringe frequencies')
plot.spectral_comparison(
gps, qspecgram, fringe, output.format('comparison'), thresh=thresh,
multipliers=multipliers, colormap=args.colormap)
plot.spectral_overlay(
gps, qspecgram, fringe, output.format('overlay'),
multipliers=multipliers)
logger.info(' -- Channel complete -- ')
count += 1 # increment counter
logger.info('{0:g} chanels plotted in {1}'.format(count, args.output_dir))


if __name__ == "__main__": # pragma: no-cover
sys.exit(main())
71 changes: 53 additions & 18 deletions gwdetchar/scattering.py → gwdetchar/scattering/core.py
Original file line number Diff line number Diff line change
Expand Up @@ -21,54 +21,89 @@

import numpy

from scipy.signal import savgol_filter

__author__ = 'Duncan Macleod <duncan.macleod@ligo.org>'

OPTIC_MOTION_CHANNELS = {
'BS': ['SUS-BS_M1_DAMP_L_IN1_DQ'],
'BS': ['SUS-BS_M1_DAMP_L_IN1_DQ',
'SUS-BS_M2_WIT_L_DQ'],
'ETMX': ['SUS-ETMX_M0_DAMP_L_IN1_DQ',
'SUS-ETMX_R0_DAMP_L_IN1_DQ'],
'SUS-ETMX_R0_DAMP_L_IN1_DQ',
'SUS-ETMX_L2_WIT_L_DQ'],
'ETMY': ['SUS-ETMY_M0_DAMP_L_IN1_DQ',
'SUS-ETMY_R0_DAMP_L_IN1_DQ'],
'SUS-ETMY_R0_DAMP_L_IN1_DQ',
'SUS-ETMY_L2_WIT_L_DQ'],
'IM1': ['SUS-IM1_M1_DAMP_L_IN1_DQ'],
'IM2': ['SUS-IM2_M1_DAMP_L_IN1_DQ'],
'IM3': ['SUS-IM3_M1_DAMP_L_IN1_DQ'],
'IM4': ['SUS-IM4_M1_DAMP_L_IN1_DQ'],
'ITMX': ['SUS-ITMX_M0_DAMP_L_IN1_DQ',
'SUS-ITMX_R0_DAMP_L_IN1_DQ'],
'SUS-ITMX_R0_DAMP_L_IN1_DQ',
'SUS-ITMX_L2_WIT_L_DQ'],
'ITMY': ['SUS-ITMY_M0_DAMP_L_IN1_DQ',
'SUS-ITMY_R0_DAMP_L_IN1_DQ'],
'MC1': ['SUS-MC1_M1_DAMP_L_IN1_DQ'],
'MC2': ['SUS-MC2_M1_DAMP_L_IN1_DQ'],
'MC3': ['SUS-MC3_M1_DAMP_L_IN1_DQ'],
'SUS-ITMY_R0_DAMP_L_IN1_DQ',
'SUS-ITMY_L1_WIT_L_DQ'],
'MC1': ['SUS-MC1_M1_DAMP_L_IN1_DQ',
'SUS-MC1_M3_WIT_L_DQ'],
'MC2': ['SUS-MC2_M1_DAMP_L_IN1_DQ',
'SUS-MC2_M3_WIT_L_DQ'],
'MC3': ['SUS-MC3_M1_DAMP_L_IN1_DQ',
'SUS-MC3_M3_WIT_L_DQ'],
'OM1': ['SUS-OM1_M1_DAMP_L_IN1_DQ'],
'OM2': ['SUS-OM2_M1_DAMP_L_IN1_DQ'],
'OM3': ['SUS-OM3_M1_DAMP_L_IN1_DQ'],
'OMC': ['SUS-OMC_M1_DAMP_L_IN1_DQ'],
'PR2': ['SUS-PR2_M1_DAMP_L_IN1_DQ'],
'PR3': ['SUS-PR3_M1_DAMP_L_IN1_DQ'],
'PRM': ['SUS-PRM_M1_DAMP_L_IN1_DQ'],
'PR2': ['SUS-PR2_M1_DAMP_L_IN1_DQ',
'SUS-PR2_M3_WIT_L_DQ'],
'PR3': ['SUS-PR3_M1_DAMP_L_IN1_DQ',
'SUS-PR3_M3_WIT_L_DQ'],
'PRM': ['SUS-PRM_M1_DAMP_L_IN1_DQ',
'SUS-PRM_M3_WIT_L_DQ'],
'RM1': ['SUS-RM1_M1_DAMP_L_IN1_DQ'],
'RM2': ['SUS-RM2_M1_DAMP_L_IN1_DQ'],
'ZM1': ['SUS-ZM1_M1_DAMP_L_IN1_DQ'],
'ZM2': ['SUS-ZM2_M1_DAMP_L_IN1_DQ'],
'OPO': ['SUS-OPO_M1_DAMP_L_IN1_DQ'],
'OFI': ['SUS-OFI_M1_DAMP_L_IN1_DQ',
'SUS-OFI_M1_DAMP_T_IN1_DQ'],
'SR2': ['SUS-SR2_M1_DAMP_L_IN1_DQ'],
'SR3': ['SUS-SR3_M1_DAMP_L_IN1_DQ'],
'SRM': ['SUS-SRM_M1_DAMP_L_IN1_DQ'],
'SR2': ['SUS-SR2_M1_DAMP_L_IN1_DQ',
'SUS-SR2_M3_WIT_L_DQ'],
'SR3': ['SUS-SR3_M1_DAMP_L_IN1_DQ',
'SUS-SR3_M3_WIT_L_DQ'],
'SRM': ['SUS-SRM_M1_DAMP_L_IN1_DQ',
'SUS-SRM_M3_WIT_L_DQ'],
'TMSX': ['SUS-TMSX_M1_DAMP_L_IN1_DQ'],
'TMSY': ['SUS-TMSY_M1_DAMP_L_IN1_DQ'],
}

FREQUENCY_MULTIPLIERS = range(1, 5)


def get_fringe_frequency(timeseries, multiplier=2.0):
"""Calculate the scattering fringe frequency from a optic motion timeseries
def get_fringe_frequency(series, multiplier=2.0):
"""Predict scattering fringe frequency from the derivative of a timeseries
Parameters
----------
series : `~gwpy.timeseries.TimeSeries`
timeseries record of relative motion
multiplier : `float`
harmonic number of fringe frequency
Returns
-------
fringef : `~gwpy.timeseries.TimeSeries`
timeseries record of fringe frequency
See Also
--------
scipy.signal.savgol_filter
for an implementation of the Savitzky-Golay filter
"""
velocity = timeseries.diff()
velocity.override_unit('m/s') # just so multiplication works
velocity = type(series)(numpy.zeros(series.size))
velocity.__array_finalize__(series)
velocity[:] = savgol_filter(series.value, 5, 2, deriv=1)
fringef = numpy.abs(multiplier * 2. / 1.064 * velocity *
velocity.sample_rate.value)
fringef.override_unit('Hz')
Expand Down
Loading

0 comments on commit 03b6dba

Please sign in to comment.