injection.py 3.06 KB
Newer Older
1
2
3
"""
Tutorial to show signal injection using new features of detector.py
"""
4
5
from __future__ import division, print_function
import numpy as np
moritz's avatar
moritz committed
6
import tupak
7
8
9
10
11
import logging


def main():

moritz's avatar
moritz committed
12
    tupak.utils.setup_logger()
13
14
15
16
17

    outdir = 'outdir'
    label = 'injection'

    # Create the waveform generator
18
19
20
21
    waveform_generator = tupak.waveform_generator.WaveformGenerator(time_duration=4, sampling_frequency=2048,
                                                                    frequency_domain_source_model=tupak.source.lal_binary_black_hole,
                                                                    parameters={'reference_frequency': 50.0,
                                                                                'waveform_approximant': 'IMRPhenomPv2'})
22
23
24
25
26

    # Define the prior
    # Merger time is some time in 2018, shame LIGO will never see it...
    time_of_event = np.random.uniform(1198800018, 1230336018)
    prior = dict()
moritz's avatar
moritz committed
27
28
29
30
31
    prior['luminosity_distance'] = tupak.prior.PowerLaw(alpha=2, minimum=100, maximum=5000, name='luminosity_distance')
    prior['geocent_time'] = tupak.prior.Uniform(time_of_event - 0.01, time_of_event + 0.01, name='geocent_time')
    prior['mass_1'] = tupak.prior.Gaussian(mu=40, sigma=5, name='mass_1')
    prior['mass_2'] = tupak.prior.Gaussian(mu=40, sigma=5, name='mass_2')
    tupak.prior.fill_priors(prior, waveform_generator)
32
33
34
35
36
37
38
39
40
41
42
43
44

    # Create signal injection
    injection_parameters = {name: prior[name].sample() for name in prior}
    if injection_parameters['mass_1'] < injection_parameters['mass_2']:
        injection_parameters['mass_1'], injection_parameters['mass_2'] =\
            injection_parameters['mass_2'], injection_parameters['mass_1']
    logging.info("Injection parameters:\n{}".format("\n".join(["{}: {}".format(key, injection_parameters[key])
                                                               for key in injection_parameters])))
    for parameter in injection_parameters:
        waveform_generator.parameters[parameter] = injection_parameters[parameter]
    injection_polarizations = waveform_generator.frequency_domain_strain()

    # Create interferometers and inject signal
moritz's avatar
moritz committed
45
    interferometers = [tupak.detector.get_inteferometer_with_fake_noise_and_injection(
46
47
48
49
            name, injection_polarizations=injection_polarizations, injection_parameters=injection_parameters,
            sampling_frequency=2048, time_duration=4, outdir=outdir) for name in ['H1', 'L1', 'V1']]

    # Define a likelihood
moritz's avatar
moritz committed
50
51
    likelihood = tupak.likelihood.MarginalizedLikelihood(interferometers, waveform_generator, prior=prior,
                                                         distance_marginalization=True, phase_marginalization=True)
52
53

    # Run the sampler
moritz's avatar
moritz committed
54
    result = tupak.sampler.run_sampler(
55
56
57
58
59
60
61
62
63
        likelihood, prior, label=label, sampler='dynesty', npoints=500, resume=False, outdir=outdir, use_ratio=True,
        injection_parameters=injection_parameters)
    truth = [injection_parameters[parameter] for parameter in result.search_parameter_keys]
    result.plot_corner(truth=truth)
    print(result)


if __name__ == "__main__":
    main()