import h5py
import numpy
import pandas
import scipy
from gwpy.timeseries import TimeSeries
from gwpy.spectrogram import Spectrogram
from gwpy.table import EventTable
from gwpy.segments import DataQualityFlag
# figure generated using mermaid
range_file = h5py.File('fig2-O2_O3_range.hdf','r')
print(range_file.keys())
print(range_file['O2'].keys())
print(range_file['O2/LIGO-Hanford'].keys())
o2_range_mpc_hanford = range_file['O2/LIGO-Hanford/range_Mpc'][:]
o2_range_times_hanford = range_file['O2/LIGO-Hanford/time_weeks'][:]
o3_range_mpc_hanford = range_file['O3/LIGO-Hanford/range_Mpc'][:]
o3_range_times_hanford = range_file['O3/LIGO-Hanford/time_weeks'][:]
range_file.close()
asd_file = h5py.File('fig2-O2_O3_asd.hdf','r')
print(asd_file.keys())
print(asd_file['O2'].keys())
print(asd_file['O2/LIGO-Hanford'].keys())
o2_asd_strain_hanford = asd_file["O2/LIGO-Hanford/strain_per_sqrt_Hz"][:]
o2_asd_freq_hanford = asd_file['O2/LIGO-Hanford/freq_Hz'][:]
o3_asd_strain_hanford = asd_file["O3/LIGO-Hanford/strain_per_sqrt_Hz"][:]
o3_asd_freq_hanford = asd_file['O3/LIGO-Hanford/freq_Hz'][:]
asd_file.close()
fig3_strain_file = h5py.File('fig3-H1-HOFT_SECTRO_NORM-1240279218-21600.hdf5','r')
fig3_sqz_file = h5py.File('fig3-H1-SQZ_SECTRO_NORM-1240279218-21600.hdf5','r')
fig3_spectrogram_strain = Spectrogram.read(fig3_strain_file)
fig3_spectrogram_sqz = Spectrogram.read(fig3_sqz_file)
fig3_strain_file.close()
fig3_sqz_file.close()
print(fig3_spectrogram_strain)
print(fig3_spectrogram_sqz)
fig4_blip_data = TimeSeries.read('fig4-blip-STRAIN_DATA-1258745539-60.hdf')
print(fig4_blip_data)
fig4_ex_loud_data = TimeSeries.read('fig4-extremely_loud-STRAIN_DATA-1266658476-60.hdf')
print(fig4_ex_loud_data)
fig4_slow_scat_data = TimeSeries.read('fig4-slow_scattering-STRAIN_DATA-1260189587-60.hdf')
print(fig4_slow_scat_data)
fig4_fast_scat_data = TimeSeries.read('fig4-fast_scattering-STRAIN_DATA-1268452970-60.hdf')
print(fig4_fast_scat_data)
fig5_glitch_rate_file = h5py.File('fig5-scattering_glitch_rate.hdf5','r')
print(fig5_glitch_rate_file.keys())
fig5_l1_pre = fig5_glitch_rate_file['L1_pre'][:]
fig5_l1_post = fig5_glitch_rate_file['L1_post'][:]
fig5_h1_pre = fig5_glitch_rate_file['H1_pre'][:]
fig5_h1_post = fig5_glitch_rate_file['H1_post'][:]
fig5_labels_text = fig5_glitch_rate_file['labels'][:]
fig5_glitch_rate_file.close()
fig5_strain_data = TimeSeries.read('fig5-STRAIN_DATA-1262389396-300.hdf')
print(fig5_strain_data)
fig5_fringe_file = h5py.File('fig5-fring_data.hdf','r')
print(fig5_fringe_file.keys())
print(fig5_fringe_file['fringe_data'].keys())
fig5_fringe_data = fig5_fringe_file['fringe_data/harmonic_1'][:]
fig5_fringe_times = fig5_fringe_file['fringe_data/times'][:]
fig5_fringe_file.close()
print(fig5_fringe_data)
print(fig5_fringe_times)
fig6_blrms_file = h5py.File('fig6-blrms_data.hdf','r')
print(fig6_blrms_file.keys())
print(fig6_blrms_file['blrms'].keys())
fig6_blrms_data_CS = fig6_blrms_file['blrms/CS'][:]
fig6_blrms_times = fig6_blrms_file['blrms/times'][:]
fig6_blrms_file.close()
print(fig6_blrms_data_CS)
print(fig6_blrms_times)
fig6_strain_data = TimeSeries.read('fig6-strain_data.hdf')
print(fig6_strain_data)
fig7_coh_freqs,fig7_coh = numpy.loadtxt('fig7-mag-mag-HxLx.txt',unpack=True)
Nsegs=3433565
fig7_mag_coupling_files = [
'fig7-H1_magnetic_coupling_20200114.txt',
'fig7-H1_magnetic_coupling_20200128.txt',
'fig7-H1_magnetic_coupling_20200204.txt',
'fig7-H1_magnetic_coupling_20200218.txt',
'fig7-H1_magnetic_coupling_20200303.txt',
'fig7-H1_magnetic_coupling_20200324.txt'
]
fig7_dat = pandas.read_csv(fig7_mag_coupling_files[0])
fig7_measured_filter = fig7_dat['flag'] == 'Measured'
fig7_upper_limit_filter = fig7_dat['flag'] == 'Upper Limit'
fig7_measured_freqs_0 = fig7_dat['frequency'][fig7_measured_filter]
fig7_measured_vals_0 = numpy.abs(fig7_dat['factor'][fig7_measured_filter])
fig7_upper_limit_freqs_0 = fig7_dat['frequency'][fig7_upper_limit_filter]
fig7_upper_limit_vals_0 = numpy.abs(fig7_dat['factor'][fig7_upper_limit_filter])
fig8_trigger_file = 'fig8-H1-GDS-CALIB_STRAIN-ETMY-PCAL-TRIGS.xml.gz'
fig8_triggers = EventTable.read(fig8_trigger_file, tablename='sngl_burst',format='ligolw')
print(fig8_triggers)
fig8_triggers['time'] = fig8_triggers['peak_time'] + fig8_triggers['peak_time_ns']*1e-9
fig8_oplev_data = TimeSeries.read('fig8-H1-SUS-ETMY_L3_OPLEV_SUM_OUT_DQ-1176946821-720.gwf',
channel='H1:SUS-ETMY_L3_OPLEV_SUM_OUT_DQ')
print(fig8_oplev_data)
# BLRMS parameters
fig8_flower = 10
fig8_fupper = 50
fig8_step = 1
fig8_oplev_data_ = fig8_oplev_data.bandpass(fig8_flower, fig8_fupper,
fstop=[fig8_flower/2., fig8_fupper*1.5],
filtfilt=False, ftype='butter').rms(fig8_step)
fig8_dq_segment_file = 'fig8-H1-DCH-ETMY_L3_OPLEV_BLRMS_GT65-1176946821-600.hdf5'
fig8_dq_segments = DataQualityFlag.read(fig8_dq_segment_file)
print(fig8_dq_segments)
fig9_histo_file = h5py.File('fig9-L1-HISTOGRAM-PYCBC-1239641067-693023.hdf','r')
print(fig9_histo_file.keys())
fig9_bins = fig9_histo_file["bins"][:]
fig9_raw_tot = fig9_histo_file["raw_hist_no_dq"][:]
fig9_raw_rem = fig9_histo_file["raw_hist_cat2"][:]
fig9_frac_tot = fig9_histo_file["norm_hist_no_dq"][:]
fig9_frac_rem = fig9_histo_file["norm_hist_cat2"][:]
fig9_histo_file.close()
print(fig9_bins)
print(fig9_raw_tot)
print(fig9_raw_rem)
print(fig9_frac_tot)
print(fig9_frac_rem)
fig10_vt_file = h5py.File('fig10-vt_ratio_data.hdf','r')
print(fig10_vt_file.keys())
print(fig10_vt_file['ifar100yr'].keys())
print(fig10_vt_file['time_loss_error'].keys())
fig10_ifar100yr = fig10_vt_file['ifar100yr/vt_ratio'][:]
fig10_ifar100yr_error = fig10_vt_file['ifar100yr/vt_ratio_error'][:]
fig10_ifar1yr = fig10_vt_file['ifar1yr/vt_ratio'][:]
fig10_ifar1yr_error = fig10_vt_file['ifar1yr/vt_ratio_error'][:]
fig10_time_error_lower = fig10_vt_file['time_loss_error/error_lower_bound'][:]
fig10_time_error_upper = fig10_vt_file['time_loss_error/error_upper_bound'][:]
fig10_vt_file.close()
print(fig10_ifar100yr)
print(fig10_ifar100yr_error)
print(fig10_time_error_lower)
print(fig10_time_error_upper)
fig11_glitch_time = 1253878751.9375
fig11_template_dur = 88.73
fig11_opt_snr_val = 29.53
fig11_snr_ts = TimeSeries.read('fig11-SNR_TS.txt')
print(fig11_snr_ts)
fig11_snr_gated_ts = TimeSeries.read('fig11-SNR_TS.txt')
print(fig11_snr_gated_ts)
fig12_strain_data = TimeSeries.read('fig12-S191110af_data.gwf',
channel='H1:DCS-CALIB_STRAIN_CLEAN_C01')
print(fig12_strain_data)
fig12_omc_data = TimeSeries.read('fig12-S191110af_data.gwf',
channel='H1:OMC-PZT1_MON_AC_OUT_DQ')
print(fig12_omc_data)
# these data files are very large
fig13_data_h1 = numpy.load('fig13-spec_10.00_2000.00_H1_7200sSFT_O3_C01_C01gated_and_C01cleangated.npz')
# Only first 20 data point are loaded in this
# example notebook due to the large size.
# Figure in paper uses the full dataset.
fig13_freqs_h1 = fig13_data_h1['freq'][:20]
fig13_asd_c01_h1 = fig13_data_h1['h1_ampwt'][:20]
fig13_asd_c01_gate_h1 = fig13_data_h1['h1_ampwt_gated'][:20]
fig13_asd_c01_clean_gate_h1 = fig13_data_h1['h1_ampwt_clean_gated'][:20]
print(fig13_freqs_h1)
print(fig13_asd_c01_h1)
print(fig13_asd_c01_gate_h1)
print(fig13_asd_c01_clean_gate_h1)
# these data files are very large
fig13_data_l1 = numpy.load('fig13-spec_10.00_2000.00_L1_7200sSFT_O3_C01_C01gated_and_C01cleangated.npz')
# Only first 20 data point are loaded in this
# example notebook due to the large size.
# Figure in paper uses the full dataset.
fig13_freqs_l1 = fig13_data_l1['freq'][:20]
fig13_asd_c01_l1 = fig13_data_l1['l1_ampwt'][:20]
fig13_asd_c01_gate_l1 = fig13_data_l1['l1_ampwt_gated'][:20]
fig13_asd_c01_clean_gate_l1 = fig13_data_l1['l1_ampwt_clean_gated'][:20]
print(fig13_freqs_l1)
print(fig13_asd_c01_l1)
print(fig13_asd_c01_gate_l1)
print(fig13_asd_c01_clean_gate_l1)
fig14_freqs,fig14_coherence = numpy.loadtxt('fig14-HL_32Hz_Coherence.dat',unpack=True)
#notches
fig14_data_gate=dict(f=fig14_freqs,Coh=fig14_coherence,NSeg=1064670)
# stochastic notch list
# https://git.ligo.org/stochastic/stochasticdetchar/-/blob/master/O3/notchlists/HL/notchlist_HL_0.03125_O3B_v10.txt
fig14_freqsToRemove =[60.0,120.0,180.0,240.0,300.0,360.0,420.0,480.0,540.0,
600.0,660.0,720.0,780.0,840.0,900.0,960.0,1020.0,1080.0,
1140.0,1200.0,1260.0,1320.0,1380.0,1440.0,1500.0,1560.0,
1620.0,1680.0,500.0,1000.0,1497.5,17.09375,17.59375,410.3125,
1083.6875,331.90625,434.90625,1083.09375,265.5625,848.9375,
575.15625,108.84375,1390.40625,52.8125,145.34375,1220.40625,
189.96875,763.84375,26.34375,31.4375,37.71875,12.4375,
1991.09375,2991.09375,234.5625,890.125,33.1875,32.0,48.0,
50.0,20.125,20.78125,21.90625,22.34375,27.46875,27.71875,
35.71875,40.90625,40.9375,60.5,61.0,61.5,62.5,63.0,63.5,
25.0,394.6875,1153.09375,31.5,32.5,33.5,34.0,30.0,40.0,100.0,
53.40625,33.5,40.0,31.5,32.5,30.0,100.0,436.53125,20.21875,
20.25,20.34375,174.5625,258.46875,276.6875,409.9375]
fig14_nBinsToRemove =[7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,
7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,1281.0,2561.0,
2721.0,33.0,33.0,33.0,33.0,33.0,33.0,33.0,3.0,7.0,5.0,1.0,35.0,
1.0,9.0,9.0,11.0,5.0,1.0,1.0,7.0,1.0,13.0,21.0,3.0,7.0,1.0,7.0,65.0,5.0,
5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0,17.0,
5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0,1.0,1.0,1.0,
1.0,1.0,1.0,1.0]
fig14_deltaF=0.03125
fig14_notches=[]
for f,n in zip(fig14_freqsToRemove,fig14_nBinsToRemove):
fig14_notches.extend([x for x in numpy.arange(f-(n-1)*fig14_deltaF/2,f+(n+1)*fig14_deltaF/2,fig14_deltaF)])
fig14_notches=numpy.sort(numpy.unique(fig14_notches))
fig14_notch_filter=[x not in fig14_notches for x in fig14_data_gate['f']]
# set up bins for histogram
fig14_bins = numpy.linspace(1e-7,1e-4,1000)
# only look at frequency band used in stochastic search
fig14_f_filter = numpy.logical_and(fig14_data_gate['f']>=20,fig14_data_gate['f']<=1726)
fig14_counts,fig14_edges = numpy.histogram(fig14_data_gate['Coh'][fig14_f_filter],bins=fig14_bins)
# apply notch list
fig14_f_filter = numpy.logical_and(fig14_f_filter,fig14_notch_filter)
fig14_counts_notched,edges=numpy.histogram(fig14_data_gate['Coh'][fig14_f_filter],bins=fig14_bins)
fig15_coherence_file = 'fig15-L1-SQZ-OMC_TRANS_RF3_Q_NORM_DQ_data.mat'
fig15_data = scipy.io.loadmat(fig15_coherence_file)
fig15_freqs = numpy.array(fig15_data['freqs'])[0]
fig15_coh = numpy.array(fig15_data['coh'])[0]
print(fig15_freqs)
print(fig15_coh)