Skip to content
Snippets Groups Projects
Commit 15054bfc authored by Wei Changfeng's avatar Wei Changfeng
Browse files

Add new file

parents
No related branches found
No related tags found
No related merge requests found
import numpy as np
import bilby
import matplotlib.pyplot as plt
from gwpy.timeseries import TimeSeries
outdir = 'outdir'
label = 'bbh_test_onewaveform'
sampling_frequency = 16384.
dt = 1.0/sampling_frequency
#Load data.
data = np.loadtxt('./mdc-master/injection-h_m1_L0.9_l2m2_r300.dat')
N = len(data)
duration = N*dt
time = np.arange(0,duration,dt)
#The strain data Daniel gave me only include one column. This is to add time column to data.
f = open('injection-h_m1_L0.9_l2m2_r300.txt','w')
for i in range(0,N):
f.write("%.20e %.20e\n" % (time[i],data[i]))
f.close()
plt.plot(time,data)
#Translate "txt" file to "gwf" file, because BILBY only reads "gwf" file.
#Read new data which has 2 columns as a timeseries.
gwffile = TimeSeries.read('injection-h_m1_L0.9_l2m2_r300.txt')
gwffile.t0 = .0
#Give channel a name.
gwffile.name='H1:GWOSC-4KHZ_R1_STRAIN'
gwffile.write('injection-h_m1_L0.9_l2m2_r300.gwf')
#This is to check if "gwf" file is succeed to produce.
#data_2=TimeSeries.read('injection-h_m1_L0.9_l2m2_r300.gwf','H1:GWOSC-4KHZ_R1_STRAIN')
# inject the signal into H1 interferometer.
#ifos = bilby.gw.detector.get_empty_interferometer('H1')
ifos = bilby.gw.detector.InterferometerList(['H1'])
ifos[0].set_strain_data_from_frame_file(frame_file='injection-h_m1_L0.9_l2m2_r300.gwf',
sampling_frequency=sampling_frequency, duration=duration,
start_time=.0,channel='H1:GWOSC-4KHZ_R1_STRAIN')
#ifos[0].set_strain_data_from_power_spectral_densities(sampling_frequency=sampling_frequency, duration=duration, start_time=-0.5)
print(ifos[0])
#Usually for 1 event, we have strain data from 2 or more interferometers, so the following code is needed:
#for i in range(len(ifos)):
# ifos[i].set_strain_data_from_frame_file(frame_file='new_injection-h_m16_L0.18_l2m2_r300.gwf',sampling_frequency=sampling_frequency, duration=duration, start_time=.0,channel='')
# ifos[i].strain_data_from_gwpy_time_series('injection-h_m16_L0.18_l2m2_r300.gwf')
# ifos[i].set_strain_data_from_power_spectral_densities(
# sampling_frequency=sampling_frequency, duration=duration,
# start_time=-0.5)
# Generate waveform
waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
reference_frequency=50., minimum_frequency=10.)
waveform_generator = bilby.gw.WaveformGenerator(
duration=duration, sampling_frequency=sampling_frequency,
frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
waveform_arguments=waveform_arguments)
#Use customized BBH priors.
priors = bilby.core.prior.PriorDict(filename='bilby_gw_prior_files_binary_black_holes.prior')
priors['geocent_time'] = bilby.prior.Uniform(
minimum=.0-0.1,
maximum=.0+0.1,
name='geocent_time', latex_label='$t_c$', unit='$s$')
#likelihood
likelihood = bilby.gw.GravitationalWaveTransient(
interferometers=ifos, waveform_generator=waveform_generator, priors=priors,
distance_marginalization=True, phase_marginalization=True, time_marginalization=True)
result = bilby.run_sampler(
likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000,
injection_parameters=None, outdir=outdir, label=label)
# Make a corner plot.
result.plot_corner()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment