%matplotlib inline
import numpy as np
import ringdown as rd
import arviz as az
from pycbc.types import TimeSeries
import h5py
from scipy.interpolate import interp1d

import matplotlib as mpl
import matplotlib.pyplot as plt

# PLOTTING OPTIONS
fig_width_pt = 3*246.0  # Get this from LaTeX using \showthe\columnwidth
inches_per_pt = 1.0/72.27               # Convert pt to inch
golden_mean = (np.sqrt(5)-1.0)/2.0         # Aesthetic ratio
fig_width = fig_width_pt*inches_per_pt  # width in inches
fig_height = fig_width*golden_mean      # height in inches
fig_size =  [fig_width,fig_height]

params = { 'axes.labelsize': 24,
          'font.family': 'serif',
          'font.serif': 'Computer Modern Raman',
          'font.size': 24,
          'legend.fontsize': 20,
          'xtick.labelsize': 24,
          'ytick.labelsize': 24,
          'axes.grid' : True,
          'text.usetex': True,
          'savefig.dpi' : 100,
          'lines.markersize' : 14,
          'figure.figsize': fig_size}

mpl.rcParams.update(params)

Plot waveform in time domain

fit = rd.Fit(model='mchi_aligned')
fit.result = az.from_netcdf('/work/yifan.wang/ringdown/GW150914/maxisi-data-release/my_gw150914_0603_srate2048_pystan.nc')

h = fit.result.posterior.h_det.stack(samples=('chain', 'draw'))
lp = fit.result.sample_stats.lp.stack(samples=('chain','draw'))
imaxlp = np.argmax(lp.values)
hmax = h.values[0,:,imaxlp]

hpycbc = TimeSeries(hmax,delta_t=1/2048)

plt.plot(hpycbc.sample_times,hpycbc)

[<matplotlib.lines.Line2D at 0x7fbc6e112048>]

Plot waveform and PSD in frequency domain

def read_strain(file, dname):
    with h5py.File(file, 'r') as f:
        t0 = f['meta/GPSstart'][()]
        T = f['meta/Duration'][()]
        h = f['strain/Strain'][:]
    
        dt = T/len(h)
    
        raw_strain = rd.Data(h, index=t0 + dt*np.arange(len(h)), ifo=dname)
        
        return raw_strain

Plot posterior of A221 for srate = 2048 and 4096

fit4096 = rd.Fit(model='mchi_aligned')
fit2048 = rd.Fit(model='mchi_aligned')

fit4096.result = az.from_netcdf('/work/yifan.wang/ringdown/GW150914/maxisi-data-release/my_gw150914_0602_srate4096.nc')
fit2048.result = az.from_netcdf('/work/yifan.wang/ringdown/GW150914/maxisi-data-release/my_gw150914_0603_srate2048_pystan.nc')

A4096 = fit4096.result.posterior.A.stack(samples=('chain','draw'))
A2048 = fit2048.result.posterior.A.stack(samples=('chain','draw'))

bins = np.linspace(0,5e-21,100)
plt.hist(A4096[1].values,alpha=0.5,bins=bins,density=True,label='srate 4096')
plt.hist(A2048[1].values,alpha=0.5,bins=bins,density=True,label='srate 2048')
plt.legend()
plt.xlabel('A221')

Text(0.5, 0, 'A221')

Plot the $2\times |h(f)| \times \sqrt{f}$

def plot_snr_eachbin(path='/work/yifan.wang/ringdown/GW150914/maxisi-data-release/my_gw150914_0602_srate4096.nc',
               srate=4096):
    
    fit = rd.Fit(model='mchi_aligned')
    fit.result = az.from_netcdf(path)
    
    #Read in strain data
    ifos = ['H1', 'L1']
    input_path = '/work/yifan.wang/ringdown/GW150914/maxisi-data-release/{i}-{i}1_GWOSC_16KHZ_R1-1126257415-4096.hdf5'
    raw_strain = {i: read_strain(input_path.format(i=i[0]), i) for i in ifos}
    for s in raw_strain.values():
        fit.add_data(s)
    T = 0.2
    fit.set_target(1126259462.4083147, ra=1.95, dec=-1.27, psi=0.82, duration=T)
    # condition data
    fit.condition_data(ds=int(round(raw_strain[ifos[0]].fsamp/srate)), flow=20)
    fit.compute_acfs()
    
    #read in maxL waveform template
    hs = fit.result.posterior.h_det.stack(samples=('chain', 'draw'))
    lp = fit.result.sample_stats.lp.stack(samples=('chain','draw'))
    imaxlp = np.argmax(lp.values)
    for ifo in ['H1','L1']:
        if ifo == 'H1':
            hsmax = hs.values[0,:,imaxlp]
            psd = fit.data['H1'].get_psd()
        elif ifo == 'L1':
            hxmas = hs.values[1,:,imaxlp]
            psd = fit.data['L1'].get_psd()
        #save the waveform into a PyCBC TimeSeries type
        h = TimeSeries(hsmax,delta_t=1/srate)
        hfreq = h.to_frequencyseries()
        #plot figures for H1 and L1 separately
        plt.figure()
        
        #psdinterp = interp1d(psd.index.values,psd.values)
        f = hfreq.sample_frequencies
        #df = f[1]-f[0]
        snr_integrand = 2 * abs(hfreq.data) * f**0.5
        plt.plot(f,snr_integrand,label=str(ifo)+' $2|h(f)|\sqrt{f}$ (maxL waveform)')
        plt.plot(psd.index.values[1:],np.sqrt(psd.values[1:]),label=str(ifo)+' ASD')
        plt.xscale('log')
        plt.yscale('log')
        plt.legend()
        #optimal_snr = np.sqrt(2*np.sum(snr_integrand))
        plt.title('sample rate '+str(srate)+' '+str(ifo))
        plt.savefig(str(ifo)+'_srate'+str(srate)+'.png',bbox_inches='tight')

plot_snr_eachbin(path='/work/yifan.wang/ringdown/GW150914/maxisi-data-release/my_gw150914_0602_srate4096.nc',
               srate=4096)

plot_snr_eachbin('/work/yifan.wang/ringdown/GW150914/maxisi-data-release/my_gw150914_0603_srate2048_pystan.nc',2048)