Merge branch '223_HAM6_pmon' into 'master'

New plot for H1 for power going into HAM6 vs after the fast shutter

Closes #223

See merge request jameson.rollins/locklost!142

locklost/plugins/__init__.py

@@ -30,6 +30,9 @@ register_plugin(discover_data)
 from .refine import refine_time
 register_plugin(refine_time)
 
+from .ham6_power import power_in_ham6
+register_plugin(power_in_ham6)
+
 from .ifo_mode import check_ifo_mode
 register_plugin(check_ifo_mode)
 

locklost/plugins/ham6_power.py

+import numpy as np
+import matplotlib.pyplot as plt
+
+from gwpy.segments import Segment
+
+from .. import logger
+from .. import config
+from .. import data
+from .. import plotutils
+
+
+#################################################
+
+
+def power_in_ham6(event):
+    '''Calculates the power under the curve for the PEM-CS_ADC_5_19_2K_OUT_DQ
+    channel and compares it to the power that gets through the fast shutter'''
+
+    plotutils.set_rcparams()
+
+    gps = event.gps
+    ifo = config.IFO
+    window = [-3, 3]
+
+    if ifo == 'L1':
+        logger.info('L1 not currently set up for plugin, plugin will not run')
+        return
+    elif ifo == 'H1':
+        if gps < 1414874845:
+            logger.info('The HAM6 power monitor was not installed for '
+                        'locklosses occurring before November 05, 2024 '
+                        '20:00UTC. Plugin not running.')
+            return
+        else:
+            if event.transition_index[0] < config.GRD_NOMINAL_STATE[0]:
+                logger.info('IFO was not fully locked, plugin not running')
+                return
+            else:
+                channels = [
+                    f'{ifo}:PEM-CS_ADC_5_19_2K_OUT_DQ',
+                    f'{ifo}:ASC-AS_A_DC_NSUM_OUT_DQ',
+                    f'{ifo}:OMC-DCPD_SUM_OUT_DQ',
+                    f'{ifo}:ISI-HAM6_GS13INF_V3_IN1_DQ'
+                ]
+
+                # Calibration
+                calPEM = 0.177  # [W/ct]
+                calASA = 1.15e-5  # From data before GPS=1415043099 lockloss
+                asa_sat_thresh = 130000*calASA  # ~1.5W; raw cts before calib
+                pmon_sat_thresh = 32000*calPEM  # ~5.66kW; rails at ~32k (1V)
+
+                # Get the data
+                segment = Segment(window).shift(gps)
+                bufs = data.fetch(channels, segment)
+
+                # Filter through each channel and get data times
+                # PEM-CS_ADC_5_19_2K_OUT_DQ
+                pmon = bufs[0].data*calPEM
+                t_pmon = np.arange(0, len(pmon))/bufs[0].sample_rate + window[0]
+                # ASC-AS_A_DC_NSUM_OUT_DQ
+                asa = bufs[1].data
+                pasa = asa*calASA
+                t_asa = np.arange(0, len(asa))/bufs[1].sample_rate + window[0]
+                # OMC-DCPD_SUM_OUT_DQ
+                dcpd = bufs[2].data
+                t_dcpd = np.arange(0, len(dcpd))/bufs[2].sample_rate + window[0]
+                # ISI-HAM6_GS13INF_V3_IN1_DQ
+                gs13 = bufs[3].data
+                t_gs13 = np.arange(0, len(gs13))/bufs[3].sample_rate + window[0]
+
+                # Find fast shutter trigger time
+                trigger = (np.abs(gs13) > 100).argmax()
+                t_trigger = t_gs13[trigger]
+
+                # Integrate power monitor over t_trigger+[-0.1, +0.2]
+                pmon_integ_range = range(
+                    (t_pmon >= t_trigger-0.1).argmax(),
+                    (t_pmon > t_trigger+0.2).argmax()
+                )
+                e_ham6 = sum(pmon[pmon_integ_range])*(t_pmon[1]-t_pmon[0])
+
+                # Assume that the power while the shutter was closed is negligible
+                # Only integrate when the shutter is open
+                # If pmon doesn't saturate, shutter open/close status is approximated by:
+                #   abs(p(ham6)/p(AS_A_DC_NSUM)) < pmon_sat_thresh/asa_sat_thresh = 3.8k
+                #
+                # This is because:
+                #   When the shutter is closed, the ratio would be:
+                #       ~1/trans_mirror ~ 1.4E4 > 3.8k
+                #       (or larger if AS_A saturates but power monitor doesn't)
+                #   When the shutter is open, the ratio would be
+                #       ~1 < 3.8k
+                #       (or larger if AS_A saturates but power monitor doesn't)
+                #   pmon saturates earlier than AS_A when shutter is closed, because
+                #    pmon_sat_thresh is 5.7kW and asa_sat_thresh/trans_mirror~20kW
+                #
+                # Therefore, unless pmon saturates,
+                #  abs(p(ham6)/p(AS_A_DC_NSUM)) < pmon_sat_thresh/asa_sat_thresh = 3.8k
+                #  is a good criteria to tell the shutter open state
+
+                # Get ratio between pmon and pasa - add 1e-20 to avoid /0
+                pmon_pasa_ratio = np.divide(np.abs(pmon), np.abs(pasa) + 1E-20)
+                ratio_thresh = pmon_sat_thresh/asa_sat_thresh
+                # Getting power downstream of fast shutter
+                p_after_FS = np.multiply(pmon, np.logical_and(pmon_pasa_ratio < ratio_thresh, pmon > 0.01))
+                # Integrate p_aft_fs over the same integration range for pmon
+                e_after_FS = sum(p_after_FS[pmon_integ_range])*(t_pmon[1]-t_pmon[0])
+                # Check if pmon is close to saturation
+                pmon_sat_check = np.max(pmon) > 0.9*pmon_sat_thresh
+
+                # Plot figures
+                logger.info('Creating HAM6 power plot')
+                # Update default font sizes
+                plt.rcParams['font.size'] = 20
+                plt.rcParams['axes.titlesize'] = 35
+                plt.rcParams['axes.labelsize'] = 20
+                plt.rcParams['xtick.labelsize'] = 20
+                plt.rcParams['ytick.labelsize'] = 20
+                plt.rcParams['legend.fontsize'] = 23
+
+                fig, (ax0, ax1, ax2) = plt.subplots(3, sharex=True, figsize=(27, 16), dpi=100)
+
+                # Plot PEM pmon, AS_A, and power after the fast shutter
+                p_peak = np.nanmax(pmon)  # peak power
+                ax0.semilogy(
+                    t_pmon,
+                    pmon,
+                    label=f'Lockloss power monitor (Energy: {e_ham6:.1f}J, Peak power: {p_peak:.1f}W)',
+                    linewidth=5.5
+                )
+                ax0.semilogy(
+                    t_asa,
+                    pasa,
+                    label='ASC-AS_A_DC_NSUM',
+                    linewidth=5.5
+                )
+                ax0.semilogy(
+                    t_pmon,
+                    p_after_FS,
+                    label=f'Power after fast shutter (Projection, Energy: {e_after_FS:.1f}J)',
+                    linewidth=2.5,
+                    color='cyan'
+                )
+                ax0.set_xlim(t_trigger-0.15, t_trigger+0.18)
+                ax0.set_ylim(1E-4, p_peak*2)
+                ax0.set_ylabel('Power [W]')
+                # Display warning if pmon saturated
+                if pmon_sat_check:
+                    logger.info('Power monitor might have railed!')
+                    ax0.text(
+                        t_trigger+0.1,
+                        10,
+                        'Power monitor might have railed!!',
+                        color='red',
+                        size=22
+                    )
+
+                # Plot DCPDs
+                ax1.plot(
+                    t_dcpd,
+                    dcpd,
+                    label='OMC-DCPD_SUM',
+                    lw=4
+                )
+                ax1.set_ylabel('Current [mA]')
+                ax1.set_ylim(35, 45)
+
+                # Plot GS13 for HAM6 that shows time of fast shutter close
+                ax2.plot(
+                    t_gs13,
+                    gs13,
+                    label='HAM6_GS13INF_V3_IN2',
+                    lw=4
+                )
+                ax2.set_ylabel('counts')
+                ax2.set_xlabel(f'Time [s] since Fast Shutter triggered at {str(gps)}')
+
+                for ax in [ax0, ax1, ax2]:
+                    ax.grid(True)
+                    ax.legend(loc='lower left')
+
+                ax0.set_title('Power through HAM6')
+                fig.tight_layout(pad=0.05)
+
+                # Save plot to lockloss directory
+                outfile_plot = 'HAM6_power.png'
+                outpath_plot = event.path(outfile_plot)
+                fig.savefig(outpath_plot, bbox_inches='tight')
+                fig.clf()

locklost/web/templates/event.tpl

@@ -125,7 +125,7 @@
 </div>
 </div>
 
-<!-- ETM_GLITCH -->
+<!-- ETM Glitch -->
 % if os.path.exists(event.path('ETM_GLITCH.png')):
     <hr />
     <div class="container">
@@ -133,7 +133,19 @@
     <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')
+%   include('collapsed_plots.tpl', title ='ETM Glitch plots', id='ETM_Glitch', plots=etm_glitch_plot, size=5, expand=has_tag, section='main')
+    </div>
+    </div>
+% end
+
+<!-- HAM6 power -->
+% if os.path.exists(event.path('HAM6_power.png')):
+    <hr />
+    <div class="container">
+    <br />
+    <div class="row">
+%   ham6_plot = [event.url('HAM6_power.png')]
+%   include('collapsed_plots.tpl', title ='HAM6 power plot', id='ham6', plots=ham6_plot, size=18, expand=False, section='main')
     </div>
     </div>
 % end