In [1]:
%matplotlib inline
In [2]:
import numpy as np
import h5py as h5
import hackh5_utils as wave;
from os import system as bash
reload(wave);
pycbclalmatplotlibnumpyCopy required files locally. THis includes all h5 and gwf(i.e. frame) files
A1. navigate to https://github.com/llondon6/nrutils_dev/tree/master/review/hackh5_review
A2. click on each file, and select the download option
- Alternatively -
* clone the entire nrutils_dev repository
B1. git clone https://github.com/llondon6/nrutils_dev.git
B2. cd nrutils_dev/review/hackh5_review/
B3. jupyter notebook
In [3]:
# Define which files to be used (use full path is the fie is not in your current directory)
# This is CFUIB0020
hacked_location = r'hacked_q1.2_dc2dcp2.h5'
unhacked_location = r'q1.2_dc2dcp2.h5'
# # This is CFUIB0029
# hacked_location = r'hacked_q1.2_base_96.h5'
# unhacked_location = r'q1.2_base_96.h5'
# Print the keys and vals in each file
dove = h5.File(unhacked_location,'r')
hack = h5.File(hacked_location,'r')
def printh5(H):
for k in [f for f in H.attrs.keys() if 'hat' in f ]:
print '%16s\t=\t%s' % ( k, H.attrs[k] )
#
print '%s\n## The CORRECT h5 File Contains:\n%s'%(60*'-',60*'-')
printh5(dove)
print '%s\n## The HACKED h5 File Contains:\n%s'%(60*'-',60*'-')
printh5(hack)
bash('cp %s ./' % hacked_location)
bash('cp %s ./' % unhacked_location)
Out[3]:
In [4]:
# These waveform parameters we taken from the XML file
# ligo.arcca.cf.ac.uk:/home/spxll/JOKI/Projects/NR2H5/pe/q1.2_dc2dcp2/07_02_2016/rerun_inc1.570796_pol1.42892_M74.097900798
# These waveform parameters are common to all PE runs performed after the
# skylocation bug fix in PyCBC.
#
# These XML files are also avilable from the DCC at
# https://dcc.ligo.org/DocDB/0122/P1500259/011/CFUIB0020_inclinations_polarizations.tar.gz
# as part of the data packages submitted for review
params1 = type('Parameters', (object,), {})
params1.approximant = 'NR_hdf5'
params1.ra = 1.94972503
params1.dec = -1.26157296
params1.polarization = 1.4398966328953218
params1.inclination = 1.5707963267948966
params1.geocentric_end_time = 1126259462.0
params1.instruments = ['H1','L1']
params1.order = 'pseudoFourPN'
params1.sample_rate = 16384
params1.delta_t = 1.0/params1.sample_rate
params1.taper = 'TAPER_NONE'
params1.f_lower = 30.1657820617
params1.distance = 40
params1.numrel_data = unhacked_location
h = h5.File(params1.numrel_data,'r')
params1.mass1 = 40.4166258924
params1.mass2 = 33.6809044161
params1.spin1x = h.attrs['spin1x_lal']
params1.spin1y = h.attrs['spin1y_lal']
params1.spin1z = h.attrs['spin1z_lal']
params1.spin2x = h.attrs['spin2x_lal']
params1.spin2y = h.attrs['spin2y_lal']
params1.spin2z = h.attrs['spin2z_lal']
In [5]:
# Set up parameters for second waveform object
params2 = type('Parameters', (object,), {})
params2.approximant = params1.approximant
params2.ra = params1.ra
params2.dec = params1.dec
params2.polarization = params1.polarization
params2.inclination = params1.inclination
params2.geocentric_end_time = params1.geocentric_end_time
params2.instruments = params1.instruments
params2.order = params1.order
params2.sample_rate = params1.sample_rate
params2.delta_t = params1.delta_t
params2.taper = params1.taper
params2.f_lower = params1.f_lower
params2.distance = params1.distance
params2.numrel_data = hacked_location
h = h5.File(params2.numrel_data,'r')
params2.mass1 = params1.mass1
params2.mass2 = params1.mass2
params2.spin1x = h.attrs['spin1x']
params2.spin1y = h.attrs['spin1y']
params2.spin1z = h.attrs['spin1z']
params2.spin2x = h.attrs['spin2x']
params2.spin2y = h.attrs['spin2y']
params2.spin2z = h.attrs['spin2z']
In [6]:
# Define configuration
params1.inclination = np.pi*0.5
params1.polarization = np.pi*0.5
params1.coa_phase = 0
params2.inclination = params1.inclination
params2.polarization = params1.polarization
params2.coa_phase = params1.coa_phase
def printipc(par):
for k in [ s for s in par.__dict__.keys() if s in ['inclination','polarization','coa_phase'] ]:
print '%s\t=\t%9.4f (radians)' % (k,par.__dict__[k])
# print '%s\t=\t%9.4f (degrees)' % (k,par.__dict__[k]*180/np.pi)
print '%s\n## The NON-HACKED Parameters are:\n%s'%(60*'-',60*'-')
printipc(params1)
print '%s\n## The HACKED Parameters are:\n%s'%(60*'-',60*'-')
printipc(params2)
# Plot TD detector response
outPath = "./mov/H1L1-h5-comparison-1.png"
wave.plot_td_waveform_resp(params1,params2,outPath=outPath,ant=False,resi=False,norm=False)
In [7]:
# Define configuration
params1.inclination = np.pi*0.5
params1.polarization = np.pi*0.5
params1.coa_phase = 0
h = h5.File(unhacked_location,'r')
hh =h5.File(hacked_location,'r')
inc,orbPhase,pol = wave.SimInspiralNRWaveformGetCorrectedRotationAnglesFromH5File(
hh,h,params1.inclination,params1.coa_phase,params1.polarization)
params2.inclination = inc
params2.polarization = pol
params2.coa_phase = orbPhase # Set to 0 to test effect of not correcting coa_phase. This is refernced in the last cell when frame files are compared.
print '%s\n## The NON-HACKED Parameters are:\n%s'%(60*'-',60*'-')
printipc(params1)
print '%s\n## The HACKED Parameters are:\n%s'%(60*'-',60*'-')
printipc(params2)
# print 'inc1 = %f' % params1.inclination
# print 'inc2 = %f' % params2.inclination
# print 'pol1 = %f' % params1.polarization
# print 'pol2 = %f' % params2.polarization
outPath = "./mov/H1L1-h5-comparison.png"
# Plot TD detector response
wave.plot_td_waveform_resp(params1,params2,outPath=outPath,ant=False,resi=False,norm=False)
In [8]:
#
ang_range = np.linspace(-np.pi,np.pi,2e3)
#
h = h5.File(unhacked_location,'r')
hh =h5.File(hacked_location,'r')
# VARYING INCLINATION
inc = np.zeros( ang_range.shape )
orb = np.zeros( ang_range.shape )
pol = np.zeros( ang_range.shape )
for k,ang in enumerate(ang_range):
params1.inclination = ang # np.pi*0.5
params1.coa_phase = 0
# The value below is chosen to coincide with worst-case scenarios
params1.polarization = np.pi/2
inc[k],orb[k],pol[k] = wave.SimInspiralNRWaveformGetCorrectedRotationAnglesFromH5File(
hh,h,params1.inclination,params1.coa_phase,params1.polarization)
#
from matplotlib.pyplot import *
pol = np.unwrap(pol)
inc = np.unwrap(inc)
#
fig = figure( figsize=6*np.array((2,1)) )
subplot(1,2,1)
plot( 180*ang_range/np.pi, 180*(inc-ang_range)/np.pi, 'm', label='shift in inclination' )
xlabel('non-hacked INCLINATION')
ylabel('degrees')
legend(frameon=False)
subplot(1,2,2)
plot( 180*ang_range/np.pi, 180*(pol-params1.polarization)/np.pi, 'b', label='shift in polarization' )
xlabel('non-hacked INCLINATION')
ylabel('degrees')
legend(frameon=False)
Out[8]:
In [9]:
#
h = h5.File(unhacked_location,'r')
hh =h5.File(hacked_location,'r')
# VARYING INCLINATION
inc = np.zeros( ang_range.shape )
orb = np.zeros( ang_range.shape )
pol = np.zeros( ang_range.shape )
for k,ang in enumerate(ang_range):
params1.polarization = ang
params1.coa_phase = 0
# The value below is chosen to coincide with worst-case scenarios
params1.inclination = np.pi*0.5
inc[k],orb[k],pol[k] = wave.SimInspiralNRWaveformGetCorrectedRotationAnglesFromH5File(
hh,h,params1.inclination,params1.coa_phase,params1.polarization)
#
from matplotlib.pyplot import *
pol = np.unwrap(pol)
inc = np.unwrap(inc)
#
fig = figure( figsize=6*np.array((2,1)) )
subplot(1,2,1)
plot( 180*ang_range/np.pi, 180*(inc-params1.inclination)/np.pi, 'm', label='shift in inclination' )
xlabel('non-hacked POLARIZATION')
ylabel('degrees')
legend(frameon=False)
subplot(1,2,2)
plot( 180*ang_range/np.pi, 180*(pol-ang_range)/np.pi, 'b', label='shift in polarization' )
xlabel('non-hacked POLARIZATION')
ylabel('degrees')
legend(frameon=False)
Out[9]:
In [10]:
# ************************************************** #
# Compare Hacked and Shifted-Unhacked frame files.
# The hacked frame files have been shifted in
# inlincation and polarization to match the unhacked, but
# since there was no shift in coa_phase, which PE optimizes
# over, the responses will differ. However, PE runs should
# be the same -- this is being checked externally.
# ************************************************** #
# Compare Hanford Frame files
frame1_location = 'H-H1HWINJ_NR_hdf5_hacked_q1.2_dc2dcp2_inc1.7570_pol1.4948_M74.097900798.gwf'
frame2_location = 'H-H1HWINJ_NR_hdf5_q1.2_dc2dcp2_inc1.5708_pol1.5708_M74.097900798.gwf'
labels = ['hacked','unhacked']
wave.plot_gwfs(frame1_location, frame2_location, labels)
# Compare Livingston Frame files
frame1_location = 'L-L1HWINJ_NR_hdf5_hacked_q1.2_dc2dcp2_inc1.7570_pol1.4948_M74.097900798.gwf'
frame2_location = 'L-L1HWINJ_NR_hdf5_q1.2_dc2dcp2_inc1.5708_pol1.5708_M74.097900798.gwf'
labels = ['hacked','unhacked']
wave.plot_gwfs(frame1_location, frame2_location, labels)
# ************************************************** #
# Repeat code from cell 3 for convineince
# ************************************************** #
# Define configuration
params1.inclination = np.pi*0.5
params1.polarization = np.pi*0.5
params1.coa_phase = 0
h = h5.File(unhacked_location,'r')
hh =h5.File(hacked_location,'r')
inc,orbPhase,pol = wave.SimInspiralNRWaveformGetCorrectedRotationAnglesFromH5File(
hh,h,params1.inclination,params1.coa_phase,params1.polarization)
params2.inclination = inc
params2.polarization = pol
params2.coa_phase = 0*orbPhase
print '%s\n## The NON-HACKED Parameters are:\n%s'%(60*'-',60*'-')
printipc(params1)
print '%s\n## The HACKED Parameters are:\n%s'%(60*'-',60*'-')
printipc(params2)
# print 'inc1 = %f' % params1.inclination
# print 'inc2 = %f' % params2.inclination
# print 'pol1 = %f' % params1.polarization
# print 'pol2 = %f' % params2.polarization
outPath = "./mov/H1L1-h5-comparison.png"
# Plot TD detector response
wave.plot_td_waveform_resp(params1,params2,outPath=outPath,ant=False,resi=False,norm=False)
In [ ]: