In [1]:
import os
import h5py
import pylab
import numpy as np
In [2]:
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)

Up to now (2022.04.28) I have performed two sets of PyCBC runs¶

sampleing rate: 1024; Using relative 221 as prior; Dynesty as sampler

  • /work/yifan.wang/ringdown/GW150914/srate1024220

  • /work/yifan.wang/ringdown/GW150914/srate1024220_221

sampleing rate: 1024; Using relative 221 as prior; cpnest as sampler; zoom in around the merger time

  • /work/yifan.wang/ringdown/GW150914/cpnest-shortstride/srate1024220

  • /work/yifan.wang/ringdown/GW150914/cpnest-shortstride/srate1024220_221

2022.08.03: I realize 1024 Hz sampling rate is too low for ringdown overtone. Lucky that I have 2048 Hz run as well

  • /work/yifan.wang/ringdown/GW150914/220_221/

  • /work/yifan.wang/ringdown/GW150914/220

In [3]:
from pycbc import detector
def tref2th1(tc,ra,dec):
    '''
    Convert time from geocenter to Hanford given ra and dec
    '''
    Det = detector.Detector('H1')
    delay = Det.time_delay_from_earth_center(ra,dec,tc)
    return tc+delay

Get results¶

In [4]:
time = np.arange(-5,11,1)
In [5]:
time
Out[5]:
array([-5, -4, -3, -2, -1,  0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10])
In [6]:
lnb_220={}
lnb_221={}
th1 = {}

for t in time:
    if t<0:
        name = 'M'+str(-t)+'MS'
    else:
        name = str(t)+'MS'
    try:
        #f220=h5py.File('../srate1024220/4ogcringdown_output/posterior_files/H1L1V1-EXTRACT_POSTERIOR_GW150914_'+name+'-1126259200-400.hdf','r')
        #f221=h5py.File('../srate1024220_221/4ogcringdown_output/posterior_files/H1L1V1-EXTRACT_POSTERIOR_GW150914_'+name+'-1126259200-400.hdf','r')
        f220 = h5py.File('../220/4ogcringdown_output/posterior_files/H1L1V1-EXTRACT_POSTERIOR_GW150914_'+name+'-1126259200-400.hdf','r')
        f221 = h5py.File('../220_221/4ogcringdown_output/posterior_files/H1L1V1-EXTRACT_POSTERIOR_GW150914_'+name+'-1126259200-400.hdf','r')
        th1[t] = tref2th1(f220.attrs['tref']+f220.attrs['toffset'],f220.attrs['ra'],f220.attrs['dec'])
        lnb_220[t]=f220.attrs['log_evidence']
        lnb_221[t]=f221.attrs['log_evidence']
    except FileNotFoundError:
        pass
In [7]:
th1
Out[7]:
{-5: 1126259462.4187012,
 -4: 1126259462.4197013,
 -3: 1126259462.4207013,
 -2: 1126259462.4217012,
 -1: 1126259462.4227014,
 0: 1126259462.4237013,
 1: 1126259462.4247012,
 2: 1126259462.4257014,
 4: 1126259462.4277012,
 5: 1126259462.4287014,
 6: 1126259462.4297013,
 7: 1126259462.4307013,
 8: 1126259462.4317012,
 9: 1126259462.4327013,
 10: 1126259462.4337013}
In [8]:
lnb_221
Out[8]:
{-5: -566041.1670293057,
 -4: -566039.7479745758,
 -3: -566037.010470901,
 -2: -566019.7984421665,
 -1: -566017.4879168574,
 0: -566019.5626423721,
 1: -566018.4890987495,
 2: -566018.3286641752,
 4: -566019.1909850541,
 5: -566018.0905366667,
 6: -566015.3941483663,
 7: -566014.0099375955,
 8: -566011.4142172891,
 9: -566011.2983283559,
 10: -566009.9052049697}
In [9]:
_,th1value = zip(*th1.items())
x221,y221=zip(*lnb_221.items())
x220,y220 = zip(*lnb_220.items())

Unfortunately these runs use too low sampling rate (1024) so is not trustable¶

np.savetxt('./pycbc_relamp221.txt',np.transpose([x,y]), header='rundir:/work/yifan.wang/ringdown/GW150914/cpnest-shortstride/\nt_H1 \t lnB')

Cotesta paper: t_Hanford = 1126259462.42323

In [10]:
filec = np.loadtxt('cotesta-B.txt')

Starting time comparison¶

In [11]:
t_co = 1126259462.42323
In [12]:
plt.plot((np.array(th1value)-t_co)*1000,np.array(y221)-np.array(y220),marker='o',label='PyCBC Dynesty Run')
plt.plot(filec[:,0],np.log(filec[:,1]),marker='o',label='Cotesta et al.')
plt.legend()
#plt.plot(0,0.4,marker='o')
plt.title('GW150914')
plt.ylabel('log$\_$evidence (221-220)')
plt.xlim(time[0],time[-1]+0.5)
plt.xlabel('$t-t_\mathrm{ref}$ ms')
Out[12]:
Text(0.5, 0, '$t-t_\\mathrm{ref}$ ms')

Zooming in¶

In [13]:
plt.plot((np.array(th1value)-t_co)*1000,np.array(y221)-np.array(y220),marker='o',label='PyCBC Dynesty Run')
plt.plot(filec[:,0],np.log(filec[:,1]),marker='o',label='Cotesta et al.')
plt.legend()
#plt.plot(0,0.4,marker='o')
plt.title('GW150914')
plt.ylabel('log$\_$evidence (221-220)')
plt.xlim(time[0],time[-1]+0.5)
plt.xlabel('$t-t_\mathrm{ref}$ ms')
plt.xlim(-2.5,2.5)
Out[13]:
(-2.5, 2.5)

Plot pyring, ringdown and pycbc in the same figure¶

In [14]:
rd = np.loadtxt('/work/yifan.wang/ringdown/GW150914/maxisi-data-release/mi-B.csv',delimiter=",")
In [15]:
import lal
tM = 69 * lal.MTSUN_SI
t0rd = 1126259462.423
t0pr = t0rd + 0.00023
In [16]:
fig = plt.figure(figsize=[15,6])

ax1 = fig.add_subplot(111)
ax2 = ax1.twiny()

#pyring
ax1.plot(filec[:,0],np.log10(filec[:,1]),marker='o',label='Cotesta et al.')

#ringdown. 0 corresponds to 1126259462.42323
ax1.plot(rd[:,0]*tM*1000-0.23,np.log10(rd[:,1]),marker='o',label='Ringdown (Isi et al.)')
#The 0 epoch in Ringdown paper corresponds to -0.23 in PyRing paper

ax1.plot((np.array(th1value)-t_co)*1000,np.log10(np.exp(np.array(y221)-np.array(y220))),marker='o',label='PyCBC Dynesty Run')

ax1.axvline(-0.23,ls='--',color='grey',label='$t^\mathrm{H1}_\mathrm{start}=1126259462.423$')
ax1.set_xlabel('$\Delta_\mathrm{start}^{H1}$ [ms]')
ax1.legend()
ax1.set_ylabel("log$_{10} \mathcal{B}^{221}_{220}$")

ax2.set_xticks(ax1.get_xticks() )
ax2.set_xbound(ax1.get_xbound())
ax2.set_xticklabels(['{:.1f}'.format(x /(tM*1000)) for x in ax1.get_xticks()])
ax2.set_xlabel('$\Delta t_\mathrm{start}^{H1} [t_\mathrm{M}]$')

#plt.savefig('comparison_addpycbc_rel221.png',bbox_inches='tight')
Out[16]:
Text(0.5, 0, '$\\Delta t_\\mathrm{start}^{H1} [t_\\mathrm{M}]$')

Zoom in¶

In [17]:
fig = plt.figure(figsize=[15,6])

ax1 = fig.add_subplot(111)
ax2 = ax1.twiny()

#pyring
ax1.plot(filec[:,0],np.log10(filec[:,1]),marker='o',label='Cotesta et al.')

#ringdown. 0 corresponds to 1126259462.42323
ax1.plot(rd[:,0]*tM*1000-0.23,np.log10(rd[:,1]),marker='o',label='Ringdown (Isi et al.)')
#The 0 epoch in Ringdown paper corresponds to -0.23 in PyRing paper

ax1.plot((np.array(th1value)-t_co)*1000,np.log10(np.exp(np.array(y221)-np.array(y220))),marker='o',label='PyCBC Dynesty Run')

ax1.axvline(-0.23,ls='--',color='grey',label='$t^\mathrm{H1}_\mathrm{start}=1126259462.423$')
ax1.set_xlabel('$\Delta_\mathrm{start}^{H1}$ [ms]')
ax1.legend()
ax1.set_ylabel("log$_{10} \mathcal{B}^{221}_{220}$")
ax1.set_xlim(-2,2)

ax2.set_xticks(ax1.get_xticks() )
ax2.set_xbound(ax1.get_xbound())
ax2.set_xticklabels(['{:.1f}'.format(x /(tM*1000)) for x in ax1.get_xticks()])
ax2.set_xlabel('$\Delta t_\mathrm{start}^{H1} [t_\mathrm{M}]$')
Out[17]:
Text(0.5, 0, '$\\Delta t_\\mathrm{start}^{H1} [t_\\mathrm{M}]$')
In [ ]: