Merge branch 'HP_ETM_glitch' into 'master'

Improved ETM_GLITCH plugin using High Passed Data

Closes #224

See merge request jameson.rollins/locklost!143

locklost/plugins/etm_glitch.py

 import numpy as np
 import matplotlib.pyplot as plt
+from scipy import signal
 
 from gwpy.segments import Segment
 
@@ -10,9 +11,9 @@ from .. import plotutils
 
 GLITCH_SEARCH_WINDOW = [-2.0, 0]
 GLITCH_CHANNEL = ['{}:SUS-ETMX_L3_MASTER_OUT_LR_DQ'.format(config.IFO)]
-THRESH = 3 * (10 ** 5)  # High threshold looking for glitches
-LOW_THRESH = (10 ** 5)  # Low threshold for number of zeroes
-LOW_TIME = 60 * (10 ** -3)  # Threshold for amount of time spent within LOW_GLITCH threshold
+THRESH = 2 * (10 ** 4)  # High threshold looking for glitches
+LOW_THRESH = 2 * (10 ** 4)  # Low threshold for number of zeroes
+LOW_TIME = 50 * (10 ** -3)  # Threshold for amount of time spent within LOW_GLITCH threshold
 
 
 ##############################################
@@ -20,10 +21,11 @@ LOW_TIME = 60 * (10 ** -3)  # Threshold for amount of time spent within LOW_GLIT
 
 def check_glitch(event):
     """
-    This will check for any glitches prior to lockloss in two steps:
-    1. Is there a glitch that goes above THRESH (3e5) in the GLITCH_SEARCH_WINDOW (2 seconds prior to lockloss)
+    This will use a 200Hz high-passed verion of GLITCH_CHANNEL to check for
+    any glitches prior to lockloss in two steps:
+    1. Is there a glitch that goes above THRESH in the GLITCH_SEARCH_WINDOW
     If this is true:
-    2. Does the channel remain within LOW_THRESH (1e5) for at least LOW_TIME (60ms)
+    2. Does the channel remain within LOW_THRESH for at least LOW_TIME
     """
     if config.IFO == 'L1':
         logger.info("NOT SET UP FOR LLO.")
@@ -41,43 +43,48 @@ def check_glitch(event):
     thresh_crossing = segment[1]
     srate = buf.sample_rate
     t = np.arange(segment[0], segment[1], 1 / srate)
+    highpass = signal.butter(4, 200, btype='highpass', output='sos', fs=srate)
+    hp_data = signal.sosfilt(highpass, buf.data)
 
     if event.refined_gps:
         lockloss_time = event.refined_gps
     else:
         lockloss_time = event.gps
-    thresh_ref = buf.data[np.where(t < lockloss_time - LOW_TIME)]
+
+    thresh_ref = hp_data[np.where(t < lockloss_time - LOW_TIME)]
     glitches = np.where((thresh_ref <= -THRESH) | (thresh_ref >= THRESH))[0]
+    # discard first 100ms of data to avoid filtering edge effects
+    glitches = glitches[glitches > 0.1 * srate]
 
     if not any(glitches):
         logger.info("No glitches found in {}".format(GLITCH_CHANNEL[0]))
-        return
-
-    glitch_time = t[glitches[0]]
-    logger.info("Glitch found at time: {}".format(glitch_time))
-    low_ref = buf.data[np.where((t >= glitch_time) & (t < lockloss_time))]
 
-    count = 0
-    low_list = []
-    for value in low_ref:
-        if -LOW_THRESH < value < LOW_THRESH:
-            count += 1
-        else:
-            # gone outside of LOW_THRESH threshold, add count to list and reset.
-            low_list.append(count)
-            count = 0
-    low_list.append(count)
-
-    thresh_crossing = min(glitch_time, thresh_crossing)
-    low_max_time = max(low_list) / srate
-    # Checking that the channel remained within LOW_THRESH (1e5) for at least LOW_TIME (60ms)
-    logger.info("After glitch, channel remained for {}ms below threshold.".format(int(1000 * low_max_time)))
-    if (low_max_time >= LOW_TIME):
-        event.add_tag('ETM_GLITCH')
     else:
-        logger.info('Channel did not remain within threshold, no ETM glitch found.')
+        glitch_time = t[glitches[0]]
+        logger.info("Threshold crossed at time: {}".format(glitch_time))
+        low_ref = hp_data[np.where((t >= glitch_time) & (t < lockloss_time))]
+
+        count = 0
+        low_list = []
+        for value in low_ref:
+            if -LOW_THRESH < value < LOW_THRESH:
+                count += 1
+            else:
+                # gone outside of LOW_THRESH threshold, add count to list and reset.
+                low_list.append(count)
+                count = 0
+        low_list.append(count)
+
+        thresh_crossing = min(glitch_time, thresh_crossing)
+        low_max_time = max(low_list) / srate
+        # Checking that the channel remained within LOW_THRESH (1e5) for at least LOW_TIME (60ms)
+        logger.info("Channel remained for {}ms below threshold.".format(int(1000 * low_max_time)))
+        if (low_max_time >= LOW_TIME):
+            event.add_tag('ETM_GLITCH')
+        else:
+            logger.info('Channel did not remain within threshold, no ETM glitch found.')
 
-    # Make plot
+    # Make plots
     fig, ax = plt.subplots(1, figsize=(22, 16))
 
     ax.plot(
@@ -87,29 +94,52 @@ def check_glitch(event):
         alpha=0.8,
         lw=2,
     )
-    ax.axhline(
+
+    ax.grid()
+    ax.set_xlabel('Time [s] since lock loss at {}'.format(event.gps), labelpad=10)
+    ax.set_ylabel('Counts')
+    ax.legend(loc='best')
+    ax.set_title('Glitch check', y=1.04)
+    ax.set_xlim(segment[0] - event.gps, segment[1] - event.gps)
+    ax.set_ylim(-3.9 * (10 ** 5), 3.9 * (10 ** 5))
+    if thresh_crossing != segment[1]:
+        plotutils.set_thresh_crossing(ax, thresh_crossing, event.gps, segment)
+
+    fig.tight_layout()
+    fig.savefig(event.path('ETM_GLITCH.png'))
+
+    # Make High Passed data plot
+    fig1, ax1 = plt.subplots(1, figsize=(22, 16))
+
+    ax1.plot(
+        t - event.gps,
+        hp_data,
+        label='200Hz HP of {}'.format(buf.channel),
+        alpha=0.8,
+        lw=2,
+    )
+    ax1.axhline(
         -THRESH,
         linestyle='--',
         color='black',
         label='Glitch threshold',
         lw=5,
     )
-    ax.axhline(
+    ax1.axhline(
         THRESH,
         linestyle='--',
         color='black',
         lw=5,
     )
-    ax.grid()
-    ax.set_xlabel('Time [s] since lock loss at {}'.format(event.gps), labelpad=10)
-    ax.set_ylabel('Counts')
-    ax.legend(loc='best')
-    ax.set_title('Glitch check', y=1.04)
-    ax.set_xlim(segment[0] - event.gps, segment[1] - event.gps)
-
-    ax.set_ylim(1.3 * -THRESH, 1.3 * THRESH)
+    ax1.grid()
+    ax1.set_xlabel('Time [s] since lock loss at {}'.format(event.gps), labelpad=10)
+    ax1.set_ylabel('Counts')
+    ax1.legend(loc='best')
+    ax1.set_title('Glitch check, 200Hz High Passed', y=1.04)
+    ax1.set_xlim(segment[0] - event.gps, segment[1] - event.gps)
+    ax1.set_ylim(1.9 * -THRESH, 1.9 * THRESH)
     if thresh_crossing != segment[1]:
-        plotutils.set_thresh_crossing(ax, thresh_crossing, event.gps, segment)
+        plotutils.set_thresh_crossing(ax1, thresh_crossing, event.gps, segment)
 
-    fig.tight_layout()
-    fig.savefig(event.path('ETM_GLITCH.png'))
+    fig1.tight_layout()
+    fig1.savefig(event.path('ETM_GLITCH_HP.png'))

locklost/web/templates/event.tpl

@@ -140,8 +140,8 @@
     <br />
     <div class="row">
 %   has_tag = event.has_tag('ETM_GLITCH')
-%   etm_glitch_plot = [event.url('ETM_GLITCH.png')]
-%   include('collapsed_plots.tpl', title ='ETM Glitch plots', id='ETM_Glitch', plots=etm_glitch_plot, size=5, expand=has_tag, section='main')
+%   etm_glitch_plots = [event.url('ETM_GLITCH.png'), event.url('ETM_GLITCH_HP.png')]
+%   include('collapsed_plots.tpl', title ='ETM Glitch plots', id='ETM_Glitch', plots=etm_glitch_plots, size=5, expand=has_tag, section='main')
     </div>
     </div>
 % end