Skip to content
Snippets Groups Projects
Commit 0c4f99c0 authored by Timothy O'Hanlon's avatar Timothy O'Hanlon
Browse files

fixed lpy and reformatted various plots

parent 46437756
No related branches found
No related tags found
No related merge requests found
Pipeline #693331 failed
......@@ -22,118 +22,143 @@ def power_in_ham6(event):
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:
if ifo == 'H1'and 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:
elif ifo == 'L1'and gps < 1410812958:
logger.info('The HAM6 power monitor was not installed for '
'locklosses occurring before Sep 19, 2024 '
'15:29 CT. Plugin not running.')
return
else:
if ifo == 'H1':
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'
]
# Get the data
segment = Segment(window).shift(gps)
bufs = data.fetch(channels, segment)
# 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)
elif ifo == 'L1':
channels = [
f'{ifo}:PEM-CS_ADC_4_11_OUT_DQ',
f'{ifo}:ASC-AS_A_DC_SUM_OUT_DQ',
f'{ifo}:OMC-DCPD_SUM_OUT_DQ',
f'{ifo}:ISI-HAM6_GS13INF_V3_IN1_DQ',
f'{ifo}:PEM-CS_ADC_4_11_GAIN'
]
# 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
# Calibration
calPEM = 1 # calibrated in front end [W]
calASA = 1.562e-5 # Calibrated to same power as PEM during lockloss
asa_sat_thresh = (2**17)*calASA # ~1.5W; raw cts before calib
pmon_sat_thresh = (2**15)*bufs[4].data[0]*calPEM # ~5.66kW; rails at ~32k (1V)
# 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
if ifo == 'H1':
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
elif ifo == 'L1':
trigger = (np.abs(gs13) > 500).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
if ifo == 'H1':
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
)
elif ifo == 'L1':
pmon_pasa_ratio = 1 # ASA already calibrated to PEM channel
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
)
if ifo == 'H1':
ax0.semilogy(
t_pmon,
p_after_FS,
......@@ -141,49 +166,59 @@ def power_in_ham6(event):
linewidth=2.5,
color='cyan'
)
ax0.axhline(
asa_sat_thresh,
label='AS_A Saturation Threshold',
color='k',
linestyle='--'
)
if trigger != 0:
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
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
)
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()
# Plot DCPDs
ax1.plot(
t_dcpd,
dcpd,
label='OMC-DCPD_SUM',
lw=4
)
ax1.set_ylabel('Current [mA]')
if ifo == 'H1':
ax1.set_ylim(35, 45)
elif ifo == 'L1':
ax1.set_ylim(15, 55)
# 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 Lockloss 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()
......@@ -65,7 +65,7 @@ def check_seismic(event):
params['threshold'],
linestyle='--',
color='indigo',
label='seismic threshold',
label='Seismic BLRMS Threshold',
lw=5,
)
# setting left y-axis paramters
......@@ -73,7 +73,7 @@ def check_seismic(event):
ax.yaxis.label.set_color('indigo')
ax.tick_params(axis='y', colors='indigo')
ax.set_ylabel('Estimated peak velocity, 30-100 mHz band [nm/s]')
ax.set_ylim(0, max(channel_data[channel].data)+1)
# ax.set_ylim(0, max(channel_data[channel].data)+1)
# setting right y-axis parameters
ax2.spines['right'].set_color('seagreen')
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment