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

fixed lpy and reformatted various plots

parent 7ce338eb
No related branches found
No related tags found
No related merge requests found
Pipeline #693419 passed
Stage: lint
Hide whitespace changes
Inline Side-by-side
locklost/plugins/ham6_power.py
+ 184
− 184
View file @ 763a7b6a
......@@ -32,192 +32,192 @@ def power_in_ham6(event):
'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)
# 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()
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)
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
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)
# 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()
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)
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,
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
p_after_FS,
label=f'Power after fast shutter (Projection, Energy: {e_after_FS:.1f}J)',
linewidth=2.5,
color='cyan'
)
if ifo == 'H1':
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.axhline(
asa_sat_thresh,
label='AS_A Saturation Threshold',
color='k',
linestyle='--'
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
)
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
)
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()
# 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()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment