Commit bb9bbc44 authored by Daniel Williams's avatar Daniel Williams
Browse files

Merge branch 'data_release' into 'datarelease' 

Data release - patch1

See merge request !3
parents 25be03e2 b585ee28

JSD/BBH/BBH_data_creation.py

0 → 100755
+134 −0
Original line number Diff line number Diff line 
#!/usr/bin/env python

# coding: utf-8



import numpy as np

import bilby

import matplotlib.pyplot as plt



duration = 1.

sampling_frequency = 256.

Nt = duration*sampling_frequency

ifos_list = ['H1','L1','V1']



# Set up a random seed for result reproducibility.  This is optional!

np.random.seed(10000)

ref_geocent_time = 1126259642.5

start_time = ref_geocent_time - 0.5



        

# The injections of reference BBH.

injection_parameters = dict(

        mass_1=78, mass_2=72, a_1=0., a_2=0., tilt_1=0., tilt_2=0.,

        phi_12=0., phi_jl=0., luminosity_distance=1400, theta_jn=1.51, psi=1.54,

        phase=0., geocent_time=1126259642.5+0.22, ra=0.89, dec=-0.94) 

        



# Fixed arguments passed into the source model.

waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',

                          reference_frequency=20, minimum_frequency=20)



# Create the waveform_generator using a LAL BinaryBlackHole source function.

waveform_generator = bilby.gw.WaveformGenerator(

    duration=duration, sampling_frequency=sampling_frequency,

    frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,

    parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,

    waveform_arguments=waveform_arguments,

    start_time=ref_geocent_time - 0.5)



# Set up interferometers.  In this case we'll use two interferometers (LIGO-Hanford (H1), LIGO-Livingston (L1). These default to their design sensitivity. They are used to produce BBH signal.

ifos = bilby.gw.detector.InterferometerList(['H1', 'L1','V1'])

ifos.set_strain_data_from_power_spectral_densities(

    sampling_frequency=sampling_frequency, duration=duration,

    start_time=ref_geocent_time - 0.5)

ifos.inject_signal(waveform_generator=waveform_generator,

                   parameters=injection_parameters)



# Transfer BBH signal in time domain to frequency domain.

freq_signal = waveform_generator.frequency_domain_strain()

signal_fd = np.zeros((3,int(Nt/2)+1),dtype=complex)

signal_td = np.zeros((3,int(Nt)))

for i in range(3):

    signal_fd[i] = ifos[i].get_detector_response(freq_signal, injection_parameters)

    signal_td[i] = sampling_frequency*np.fft.irfft(signal_fd[i])



# Calculate the optimal SNR in each detector.



np.sqrt(ifos[0].optimal_snr_squared(signal=signal_fd[0])).real

np.sqrt(ifos[1].optimal_snr_squared(signal=signal_fd[1])).real

np.sqrt(ifos[2].optimal_snr_squared(signal=signal_fd[2])).real





# Signal whitening.

def whiten_signal(signal_fd,ifo):

    

    whitened_signal_fd = signal_fd/ifo.amplitude_spectral_density_array

    

    # The following signal whitening will contain 100 sets of white noise since we don't fix the random seed here.

    whitened_signal_td = np.sqrt(2.0*Nt)*np.fft.irfft(whitened_signal_fd) + 1.0*np.random.randn(256)

    

    return whitened_signal_td



# In this study we generate 100 data files, which contain one signal combining 100 noise realisations.

num_samples=100

# Create the inteteferometers to extract their ASDs.

np.random.seed(10000)

ifos_n = bilby.gw.detector.InterferometerList(ifos_list)

ifos_n.set_strain_data_from_power_spectral_densities(

            sampling_frequency=256, duration=duration,

            start_time=ref_geocent_time-duration/2.0)



whitened_signal_td_array=np.zeros((num_samples,3,256))





np.random.seed(10000)

for i in range(num_samples):

    k=0

    for k in range(len(ifos_list)):

        whitened_signal_td_array[i][k,:] = whiten_signal(signal_fd[k],ifos_n[k]) 

        k+=1

    



# Plot the noisy signal.

times = np.linspace(start_time, start_time+duration, int(Nt))

plt.plot(times,whitened_signal_td_array[0][0])

plt.xlabel('GPS time[s]', fontdict={'family':'Times New Roman', 'size':18})

plt.ylabel('strain', fontdict={'family':'Times New Roman', 'size':18})

plt.grid(color="k", linestyle=":")

plt.savefig('noisy_BBH.pdf')





# import the h5py module for creating data file.

import h5py



f1 = h5py.File('data_0.h5py','r')

name_list=[]

data_set_list=[]



for name in f1:

    name_list.append(name)

    data_set_list.append(f1[f'{name}'])



x_data_list=[injection_parameters['mass_1'],injection_parameters['mass_2'],injection_parameters['luminosity_distance'],

        injection_parameters['geocent_time'],injection_parameters['ra'],injection_parameters['dec'],

        injection_parameters['psi']]



# Record the injections in the data file.

x_data_sum = np.zeros((num_samples,7))

for i in range(num_samples):

    x_data_sum[i] = x_data_list

x_data_sum



# Produce the data file for VItamin analysis.

for i in range(num_samples):

    

    with h5py.File(f'data_20211025{i}.h5py','w') as f1:

        

        for k in range(len(name_list)):

            d = f1.create_dataset(name_list[k],data=data_set_list[k])

                

        del f1['x_data']

        del f1['y_data_noisy']

                

        d1 = f1.create_dataset('x_data',data = x_data_sum[i])

        d2 = f1.create_dataset('y_data_noisy',data = whitened_signal_td_array[i])

JSD/BH_capture/BH_capture_data_creation.py

0 → 100755
+261 −0
Original line number Diff line number Diff line 
#!/usr/bin/env python

# coding: utf-8



import minke.sources

import lalsimulation

import lal

import numpy as np

import minke.distribution

import matplotlib.pyplot as plt

import bilby

get_ipython().run_line_magic('matplotlib', 'inline')



# In this study we generate one signal that combines with 100 noise realisations.

num_samples = 100



# The BH capture waveform sets we use in the data production.

nr_path = './h_psi4/{}'

waveforms_merge = {

    1: "h_m1_L0.9_l2m2_r300.dat",

    2: "h_m2_L0.87_l2m2_r300.dat",

    4: "h_m4_L0.5_l2m2_r300.dat",

    8: "h_m8_L0.35_l2m2_r280.dat",

    16: "h_m16_L0.18_l2m2_r300.dat"

}

# The extraction is corresponding to the waveform.

extraction = {

    1: 300,

    2: 300,

    4: 300,

    8: 280,

    16: 300

}





# Specify the waveform by mass ratio.

mass_ratio= 1

nr_waveform = waveforms_merge[mass_ratio]



# Set the injections.

duration = 1

sampling_frequency = 256

Nt=duration*sampling_frequency



ref_geocent_time = 1126259642.5

start_time = ref_geocent_time-duration/2.0

epoch = str(start_time)



# These are the distance injections we obtain from VItamin analysis. 

distance_set = {

    1: 5000,

    4: 2000,

    8: 1500,

    16: 500

}





# Only one set of sky location injections is used in this script. 

total_mass = 150

ra = np.random.uniform(0.89,0.89,num_samples)

dec = np.random.uniform(-0.94,-0.94,num_samples)

psi = np.random.uniform(1.54,1.54,num_samples)

locations = (ra,dec,psi)

distances = np.random.uniform(distance_set[mass_ratio], distance_set[mass_ratio], num_samples)





# the function to generate an injection for a given waveform.

def generate_for_detector(source, ifos, sampling_frequency, epoch, distance, total_mass, ra, dec, psi):

    """

    Generate an injection for a given waveform.

    

    Parameters

    ----------

    source : ``minke.Source`` object

       The source which should generate the waveform.

    ifos : list

       A list of interferometer initialisms, e.g. ['L1', 'H1']

    sampling_frequency : int

       The sampling_frequency in hertz for the generated signal.

    epoch : str

       The epoch (start time) for the injection.

       Note that this should be provided as a string to prevent overflows.

    distance : float

       The distance for the injection, in megaparsecs.

    total_mass : float

       The total mass for the injected signal.

    ra : float

        The right ascension of the source, in radians.

    dec : float

        The declination of the source, in radians.

    psi : float

        The polarisation angle of the source, in radians.

        

    Notes

    -----

    This is very much a draft implementation, and there are a number of 

    things which should be tidied-up before this is moved into Minke,

    including bundling total mass with the source.

    """

    nr_waveform = source.datafile

    

    waveform = {}

    waveform['signal'] = {}

    waveform['times'] = {}

    

    for ifo in ifos:

            

        det = lalsimulation.DetectorPrefixToLALDetector(ifo)

        hp, hx = source._generate(half=True, epoch=epoch, rate=sampling_frequency)[:2]

        h_tot = lalsimulation.SimDetectorStrainREAL8TimeSeries(hp, hx, ra, dec, psi, det)

        waveform['signal'][ifo] = h_tot.data.data.tolist()

        waveform['times'][ifo] = np.linspace(0, len(h_tot.data.data)*h_tot.deltaT, len(h_tot.data.data)).tolist()

                 

    return waveform





waveforms = []



# Set the detector network we use.

ifos_list = ['H1','L1','V1']



# Generate the mock signal.

for distance, location in zip(distances, np.vstack(locations).T):

    

    hyper_object = minke.sources.Hyperbolic(

    datafile = nr_path.format(nr_waveform),

    total_mass = total_mass,

    distance = distance,

    extraction=extraction[mass_ratio]

    )

    hyper_object.datafile = nr_waveform

    

    waveform = generate_for_detector(hyper_object, ifos_list, sampling_frequency, epoch, distance, total_mass, *location)

    waveforms.append(waveform)





# The start time was specified when generating mock signal, while the merger time t0 and signal length could not be. We do the following process to make the signal have a start time of 1126259642.5s in GPS time, and a t0 of 0.22s. 



peak_index = waveforms[0]['signal']['H1'].index(max(waveforms[0]['signal']['H1']))

supposed_peak_index = int((0.5+0.22)/duration*sampling_frequency)

delta_peak = supposed_peak_index - peak_index

delta_start_time = delta_peak/Nt

new_start_time = start_time + delta_start_time

new_epoch = str(new_start_time)



new_waveforms = []



# By reproducing mock signal with a new start time, we make the signal have a supposed t0=0.22s, while maintaining a correct start time.

for distance, location in zip(distances, np.vstack(locations).T):

    

    hyper_object = minke.sources.Hyperbolic(

    datafile = nr_path.format(nr_waveform),

    total_mass = total_mass,

    distance = distance,

    extraction=extraction[mass_ratio]

    )

    hyper_object.datafile = nr_waveform

    

    new_waveform = generate_for_detector(hyper_object, ifos_list, sampling_frequency, new_epoch, distance, total_mass, *location)

    new_waveforms.append(new_waveform)

    

   



# Pad the waveform to fit the restrict of VItamin.

def fix_waveform(ifos,num_samples,new_waveforms,delta_peak):

    

    fixed_waveforms = []

    for i in range(num_samples):

        fixed_waveform = {}

        fixed_waveform['signal'] = {}

        for ifo in ifos_list:            

            peak_index = waveforms[i]['signal'][ifo].index(max(waveforms[i]['signal'][ifo]))

            new_peak_index = new_waveforms[i]['signal'][ifo].index(max(new_waveforms[i]['signal'][ifo]))

            if new_peak_index > peak_index:

                new_waveforms[i]['signal'][ifo] = new_waveforms[i]['signal'][ifo][new_peak_index-peak_index:]

            if new_peak_index < peak_index:

                for i in range(peak_index-new_peak_index):

                    new_waveforms[i]['signal'][ifo] = [0] + new_waveforms[i]['signal'][ifo]           

            

            

            fixed_waveform['signal'][ifo] = new_waveforms[i]['signal'][ifo]   

            for j in range(delta_peak):

                fixed_waveform['signal'][ifo] = [0] + fixed_waveform['signal'][ifo]           

            for k in range(Nt-len(fixed_waveform['signal'][ifo])):

                fixed_waveform['signal'][ifo] = fixed_waveform['signal'][ifo] + [0]        

            fixed_waveform['signal'][ifo] = fixed_waveform['signal'][ifo][0:Nt]            

                

        fixed_waveforms.append(fixed_waveform)

        

    return fixed_waveforms    





fixed_waveforms = fix_waveform(ifos_list,num_samples,new_waveforms,delta_peak)



np.random.seed(10000)

ifos={}

# Use the detector simulated ASD to whiten signal.

ifos = bilby.gw.detector.InterferometerList(ifos_list)

ifos.set_strain_data_from_power_spectral_densities(

            sampling_frequency=sampling_frequency, duration=duration,

            start_time=ref_geocent_time-duration/2.0)



def whiten_signal(signal,ifo):

    

    signal_fd = np.fft.rfft(signal)/Nt

    whitened_signal_fd = signal_fd/ifo.amplitude_spectral_density_array



    # The following signal whitening will contain 100 sets of white noise since we don't fix the random seed here.

    whitened_signal_td = np.sqrt(2.0*Nt)*np.fft.irfft(whitened_signal_fd) + 1.0*np.random.randn(int(Nt))

    

    return whitened_signal_td



np.random.seed(10000)

whitened_signal_td_array=np.zeros((num_samples,3,int(Nt)))



for i in range(num_samples):

    k = 0

    for ifo in ifos_list:

        whitened_signal_td_array[i][k,:] = whiten_signal(fixed_waveforms[i]['signal'][ifo],ifos[k]) 

        k += 1



# Plot the noisy capture signal.

plt.plot(times,whitened_signal_td_array[0][0])

plt.xlabel('GPS time[s]', fontdict={'family':'Times New Roman', 'size':18})

plt.ylabel('strain', fontdict={'family':'Times New Roman', 'size':18})

plt.grid(color="k", linestyle=":")

plt.savefig('noisy_capture.pdf')



# Import the h5py module for creating data file.

import h5py



f1 = h5py.File('data_0.h5py','r')

name_list=[]

data_set_list=[]



for name in f1:

    name_list.append(name)

    data_set_list.append(f1[f'{name}'])



# Calculate m1,m2 injection.

m2=total_mass/(mass_ratio+1)

m1=total_mass-m2



# Record the injections in the data file.

x_data_sum = np.zeros((num_samples,7))

for i in range(num_samples):

    x_data_sum[i] = [m1,m2,distance_set[mass_ratio],0.22,locations[2][i],locations[0][i],locations[1][i]]

x_data_sum



# Produce the data file for VItamin analysis.

for i in range(num_samples):

    

    with h5py.File(f'data_20211013{i}.h5py','w') as f1:

        

        for k in range(len(name_list)):

            d = f1.create_dataset(name_list[k],data=data_set_list[k])

        

        del f1['x_data']

        del f1['y_data_noisy']

        

        d1 = f1.create_dataset('x_data',data = x_data_sum[i])

        d2 = f1.create_dataset('y_data_noisy',data = whitened_signal_td_array[i])

JSD_calculation.py→JSD/JSD_calculation.py

100644 → 100755
+126 −234

File changed and moved.File mode changed from 100644 to 100755.

Preview size limit exceeded, changes collapsed.

README.md

+10 −9
Original line number Diff line number Diff line
@@ -41,10 +41,11 @@ The VItamin analysis outputs the posterior corner plot for the input data. 




## BH capture data creation



We generate BH capture data using the 'vitamin/BH_capture_data_creation.ipynb' notebook.

This notebook also includes code for signal whitening and generating noise realisations, in order to produce timeseries which are suitable for analysis with bilby.

In this notebook, the distributions of sky location are specified and contain 100 samples, the distance injection and the noise realisation is fixed, and the waveform is specified. 

We generate BH capture data using the 'vitamin/BH_capture_data_creation.py' script. 

To create a data file in VItamin format, we need to load the example data file 'data_0.h5py'. This helps us avoid manually creating many attributes that are unrelated to the real analysis. The example data file can be downloaded in this page:

https://hagabbar.github.io/vitamin_c/quickstart.html

The script also includes code for signal whitening and generating noise realisations, in order to produce timeseries which are suitable for analysis with VItamin.

In this script, the distributions of sky location are specified and contain 100 samples, the distance injection and the noise realisation is fixed, and the waveform is specified. 

Therefore, the script generates 100 data files per waveform. 

Then we perform VItamin analysis on them, and look for the posterior which is visually similar to that of BBH. If there is not, we adjust distance injection for another turn of data creation and make the following analysis, until obtaining the appropriate posterior.

The process should also be performed for other waveforms with different mass ratios.
@@ -52,10 +53,10 @@ The process should also be performed for other waveforms with different mass rat 
## Bilby analyses



In this section, we introduce the code realted to Bilby part of our work, which are put in the 'bilby' directory.

We employ 'bilby/PE_dynesty.py' to apply Bayesian inference on BH capture data with the dynesty sampler.

Using this script, we first reproduce one BH capture signal; then inject the signal into simulated detector noise.

We employ 'bilby/non-spinning_analysis.py' and 'bilby/spinning_analysis.py' to apply Bayesian inference on BH capture data with non-spinning and spinning BBH model respectively.

Using these scripts, we first reproduce one BH capture signal; then inject the signal into simulated detector noise.

We then specify the BBH model, prior and likelihood, and finally perform parameter estimation on it. 

We can easily switch the spinning model to the non-spinning one, by setting a delta distribution prior of six spins at zero.

Actually, we can easily switch the spinning model to the non-spinning one, by setting a delta distribution prior of six spins at zero.



## JS divergence analysis
@@ -71,10 +72,10 @@ We add the following codes in line 1491: 
Then VItamin can generate both corner plot and sample data file for the input data.



### BH capture data creation and analysis

We use the script `BH_capture_data_creation.ipynb`, from `jsd` directory, to inject one BH capture signal into 100 noise realisations, and finally generate 100 data file. We do this for all waveforms and obtain their samples data files by VItamin analysis.

We use the script `BH_capture_data_creation.py`, from `jsd` directory, to inject one BH capture signal into 100 noise realisations, and finally generate 100 data file. We do this for all waveforms and obtain their samples data files by VItamin analysis.



### BBH data creation and analysis

We use the script `BBH_data_creation.ipynb`, from `jsd` directory, to inject one BBH signal into 100 noise realisations, and finally generate 100 data files. 

We use the script `BBH_data_creation.py`, from `jsd` directory, to inject one BBH signal into 100 noise realisations, and finally generate 100 data files. 

In this analysis, we have two kinds of BBH file to produce. 

First, it's the reference BBH, of which the injection is set in advance, and can be easily generate with 100 noise realisations. 

The other one is the BBH recovered by the average peak values of posterior in Sec 3.2.

bilby/non-spinning_analysis.py

0 → 100755
+236 −0
Original line number Diff line number Diff line 
import minke.sources

import lalsimulation

import lal

import numpy as np

import minke.distribution

import matplotlib.pyplot as plt

import bilby



outdir = 'outdir_BH_encounter'

label = 'BH_encounter_0907'

bilby.core.utils.setup_logger(outdir=outdir, label=label)





# In the Bilby analysis we only use one data.

num_samples = 1



# The BH capture waveform sets we use in the data production..

nr_path = '../h_psi4/{}'

waveforms_merge = {

    1: "h_m1_L0.9_l2m2_r300.dat",

    2: "h_m2_L0.87_l2m2_r300.dat",

    4: "h_m4_L0.5_l2m2_r300.dat",

    8: "h_m8_L0.35_l2m2_r280.dat",

    16: "h_m16_L0.18_l2m2_r300.dat"

}



# The extraction is corresponding to the waveform.

extraction = {

    1: 300,

    2: 300,

    4: 300,

    8: 280,

    16: 300

}



# Specify the waveform by mass ratio.

mass_ratio= 1

nr_waveform = waveforms_merge[mass_ratio]



# Set the injections.

duration = 1

sampling_frequency = 256

Nt=duration*sampling_frequency

t0=0.22



ref_geocent_time = 1126259642.5

start_time = ref_geocent_time-duration/2.0

epoch = str(start_time)



total_mass = 150

locations=(np.array([3.69]),np.array([-0.94]),np.array([1.54]))

distances=np.array([5000.0])



# Generate BH capture signal.

def generate_for_detector(source, ifos, sampling_frequency, epoch, distance, total_mass, ra, dec, psi):

    """

    Generate an injection for a given waveform.

    

    Parameters

    ----------

    source : ``minke.Source`` object

       The source which should generate the waveform.

    ifos : list

       A list of interferometer initialisms, e.g. ['L1', 'H1']

    sampling_frequency : int

       The sampling_frequency in hertz for the generated signal.

    epoch : str

       The epoch (start time) for the injection.

       Note that this should be provided as a string to prevent overflows.

    distance : float

       The distance for the injection, in megaparsecs.

    total_mass : float

       The total mass for the injected signal.

    ra : float

        The right ascension of the source, in radians.

    dec : float

        The declination of the source, in radians.

    psi : float

        The polarisation angle of the source, in radians.

        

    Notes

    -----

    This is very much a draft implementation, and there are a number of 

    things which should be tidied-up before this is moved into Minke,

    including bundling total mass with the source.

    """

    nr_waveform = source.datafile

    

    waveform = {}

    waveform['signal'] = {}

    waveform['times'] = {}

    

    for ifo in ifos:

            

        det = lalsimulation.DetectorPrefixToLALDetector(ifo)

        hp, hx = source._generate(half=True, epoch=epoch, rate=sampling_frequency)[:2]

        h_tot = lalsimulation.SimDetectorStrainREAL8TimeSeries(hp, hx, ra, dec, psi, det)

        waveform['signal'][ifo] = h_tot.data.data.tolist()

        waveform['times'][ifo] = np.linspace(0, len(h_tot.data.data)*h_tot.deltaT, len(h_tot.data.data)).tolist()

            

    return waveform





waveforms = []



# Set the detector network we use.

ifos_list = ['H1','L1','V1']



# Generate the mock signal.

for distance, location in zip(distances, np.vstack(locations).T):

    

    hyper_object = minke.sources.Hyperbolic(

    datafile = nr_path.format(nr_waveform),

    total_mass = total_mass,

    distance = distance,

    extraction=300

    )

    hyper_object.datafile = nr_waveform

    

    waveform = generate_for_detector(hyper_object, ifos_list, sampling_frequency, epoch, distance, total_mass, *location)

    waveforms.append(waveform)





# The start time was specified when generating mock signal, while the merger time t0 and signal length could not be. We do the following process to make the signal have a start time of 1126259642.5s in GPS time, and a t0 of 0.22s. 



peak_index = waveforms[0]['signal']['H1'].index(max(waveforms[0]['signal']['H1']))

supposed_peak_index = int((0.5+0.25)/duration*256)

delta_peak = supposed_peak_index - peak_index

delta_start_time = delta_peak/256

new_start_time = start_time + delta_start_time

new_epoch = str(new_start_time)



new_waveforms = []



# By reproducing mock signal with a new start time, we make the signal have a supposed t0=0.22s, while maintaining a correct start time.

for distance, location in zip(distances, np.vstack(locations).T):

    

    hyper_object = minke.sources.Hyperbolic(

    datafile = nr_path.format(nr_waveform),

    total_mass = total_mass,

    distance = distance,

    extraction=300

    )

    hyper_object.datafile = nr_waveform

    

    new_waveform = generate_for_detector(hyper_object, ifos_list, sampling_frequency, new_epoch, distance, total_mass, *location)

    new_waveforms.append(waveform)





# Pad the waveform to fit the restrict of VItamin.

def fix_waveform(ifos,num_samples,waveforms,delta_peak):

    

    fixed_waveforms = []

    for i in range(num_samples):

        fixed_waveform = {}

        fixed_waveform['signal'] = {}

        for ifo in ifos_list:

            fixed_waveform['signal'][ifo] = waveforms[i]['signal'][ifo]   

            for j in range(delta_peak):

                fixed_waveform['signal'][ifo] = [0] + fixed_waveform['signal'][ifo]           

            fixed_waveform['signal'][ifo] = fixed_waveform['signal'][ifo][0:256]   

                

        fixed_waveforms.append(fixed_waveform)

        

    return fixed_waveforms    

    

fixed_waveforms = fix_waveform(ifos_list,num_samples,waveforms,delta_peak)





# Specify the model that we use to perform parameter estimation

waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',

                          reference_frequency=20., minimum_frequency=20.)

waveform_generator = bilby.gw.WaveformGenerator(

    duration=duration, sampling_frequency=sampling_frequency,

    frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,

    parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,

    waveform_arguments=waveform_arguments,

    start_time=ref_geocent_time - 0.5)



# Set the detector network and the noise.

np.random.seed(10000)

ifos = bilby.gw.detector.InterferometerList(ifos_list)

ifos.set_strain_data_from_power_spectral_densities(

            sampling_frequency=sampling_frequency, duration=duration,

            start_time=ref_geocent_time-duration/2.0)



# Set the priors. For the non-spinning BBH analysis, we set prior of six spins as zero.

priors = bilby.gw.prior.BBHPriorDict()

priors.pop('chirp_mass')

priors.pop('mass_ratio')



priors['mass_1'] = bilby.core.prior.Uniform(name='mass_1', minimum=30, maximum=160, unit='$M_{\\odot}$')

priors['mass_2'] = bilby.core.prior.Uniform(name='mass_2', minimum=30, maximum=160, unit='$M_{\\odot}$')

priors['luminosity_distance'] = bilby.core.prior.Uniform(name='luminosity_distance', minimum=1000, maximum=3000, unit='$Mpc_{\\odot}$')

priors['psi'] = bilby.core.prior.Uniform(name='psi', minimum=0.0, maximum=np.pi, boundary='periodic')

for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl']:

    priors[key] = 0.0



priors['geocent_time'] = bilby.core.prior.Uniform(

    minimum=ref_geocent_time + 0.15,

    maximum=ref_geocent_time + 0.35,

    name='geocent_time', latex_label='$t_c$', unit='$s$')





# Transfter signal in time domain to frequency domain, then inject it into the simulated detector noise.

def fft_signal(signal):

    

    signal_fd = np.fft.rfft(signal)/Nt

    

    return signal_fd





signal_fd_array=np.zeros((3,129))

k = 0

for ifo in ifos_list:

    signal_fd_array[k,:] = fft_signal(fixed_waveforms[0]['signal'][ifo]) 

    k += 1



k = 0

for det in ifos_list:

    ifos[k].set_strain_data_from_frequency_domain_strain(signal_fd_array[k],

                                                    sampling_frequency=sampling_frequency,

                                                    duration=duration,

                                                    start_time=start_time)

    k += 1



# Set the likelihood.

likelihood = bilby.gw.GravitationalWaveTransient(

    interferometers=ifos, waveform_generator=waveform_generator,priors=priors)



# We use dynesty as the sampler and perform parameter estimation.

result = bilby.run_sampler(

    likelihood=likelihood, priors=priors, label=label,outdir=outdir,sampler='dynesty', 

    nlive=1000, sample='rwalk',walks =100,nact=5,check_point_delta_t=1800,check_point_plot=True)  



result.plot_corner()