set_injection_parameters.py 10.2 KB
Newer Older
Marc Arene's avatar
Marc Arene committed
1 2 3
import sys
# cf https://docs.python.org/3/tutorial/modules.html
sys.path.append('../')
4
import numpy as np
5 6
from configparser import ConfigParser

7
import bilby.gw.conversion as conv
8

9
from Headers.PN_Coefficients import * # Headers.Constants already imported in PN_Coeff
10
import Library.python_utils as pu
11
import Library.param_utils as paru
Marc Arene's avatar
Marc Arene committed
12

Marc Arene's avatar
Marc Arene committed
13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
import Library.config as conf
import Library.roq_utils as roqu


LAL_PI = 3.141592653589793238462643383279502884
LAL_MSUN_SI = 1.988546954961461467461011951140572744e30
LAL_MTSUN_SI = 4.925491025543575903411922162094833998e-6

def Compute_fISCO(Rmin, m1, m2):
    """
    Compute the frequency at the Innermost Stable Circular Orbit
    """

    vISCO = np.sqrt(1.0 / Rmin)
    Mt = m1 + m2
    fmax = vISCO**3 / (PI * GEOM * Mt)      # fmax signal
    return fmax

def Compute_LAL_fISCO(m1_SI, m2_SI):
    m1 = m1_SI / LAL_MSUN_SI
    m2 = m2_SI / LAL_MSUN_SI
    m_sec = (m1 + m2) * LAL_MTSUN_SI
    piM = LAL_PI * m_sec
    vISCO = 1 / np.sqrt(6)
    fISCO = vISCO**3 / piM

    return fISCO
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77


def ComputeChirpTime3p5PN(f_low, m1, m2):
    """
    Computes chirp times at 3.5 PN.

    Inputs:
    -------
        f_low: [float] low frequency cutoff of the detector
        sigpar: [array_like] parameters of the system
    Outputs:
    --------

    """
    Mt = m1 + m2
    Mt_sec = Mt * GEOM           # Total mass in s
    etaM = m1 * m2 / Mt**2
    etaM2 = etaM*etaM

    Tau_0 = tf1 * Mt_sec / etaM
    Tau_2 = tf2 + tf3 * etaM
    Tau_3 = tf4 * PI
    Tau_4 = tf5 + tf6 * etaM + tf7 * etaM2
    Tau_5 = PI*(tf9 + tf8*etaM)
    Tau_6 = tf11 + tf12*PI2 + tf10*(EULER + np.log(4.0)) + ( tf13 + tf14*PI2 )*etaM + tf15*etaM2 + tf16*etaM2*etaM
    Tau_6_log = tf10
    Tau_7 = PI*(tf19 + tf18*etaM + tf17*etaM2 )
    v_0 = pow(PI*Mt_sec*f_low,1.0/3.0)
    v_0_2 = v_0*v_0
    v_0_3 = v_0_2*v_0
    v_0_4 = v_0_3*v_0
    v_0_5 = v_0_4*v_0
    v_0_6 = v_0_5*v_0
    v_0_7 = v_0_6*v_0

    Tau_chirp = Tau_0*pow(v_0,-8.0)*(1.0 + Tau_2*v_0_2 + Tau_3*v_0_3 + Tau_4 *v_0_4 +  Tau_5 *v_0_5 +  (Tau_6 + Tau_6_log*np.log(v_0))*v_0_6 +  Tau_7 *v_0_7 )
    return Tau_chirp

Marc Arene's avatar
Marc Arene committed
78 79 80 81 82 83 84 85 86 87 88 89 90 91
def compute_time_from_flow_to_fhigh_3p5PN(f_low, f_high, mass_1, mass_2):
    """
    Computes the time it takes for the binary to go from f_low to f_high, using 3.5 PN approximation.
    """

    # Compute the time it takes from flow to get to merger
    tc_flow = ComputeChirpTime3p5PN(f_low, mass_1, mass_2)

    # Compute the time it takes from fhigh to get to merger
    tc_fhigh = ComputeChirpTime3p5PN(f_high, mass_1, mass_2)

    # Substracting the two gives the time it take for the binary to go from f_low to f_high

    return tc_flow - tc_fhigh
92

93 94 95 96 97 98 99 100
def ini_file_to_dict(file_path='../examples/GW170817.ini'):

    config = ConfigParser()
    config.read(file_path)
    injection_parameters = pu.config_parser_to_dict(config)['parameters']


    injection_parameters = conv.generate_mass_parameters(injection_parameters)
101 102 103
    # injection_parameters['reduced_mass'] = (injection_parameters['mass_1'] * injection_parameters['mass_2']) / (injection_parameters['mass_1'] + injection_parameters['mass_2'])
    injection_parameters['reduced_mass'] = paru.component_masses_to_reduced_mass(injection_parameters['mass_1'], injection_parameters['mass_2'])

104 105 106 107 108 109 110


    # Set the duration and sampling frequency of the data segment that we're going
    # to inject the signal into. For the
    # TaylorF2 waveform, we cut the signal close to the isco frequency
    m1 = injection_parameters['mass_1']
    m2 = injection_parameters['mass_2']
Marc Arene's avatar
Marc Arene committed
111
    minimum_frequency = conf.minimum_frequency
112 113 114 115 116
    tc_3p5PN = ComputeChirpTime3p5PN(minimum_frequency, m1, m2)
    # We integrate the signal up to the frequency of the "Innermost stable circular orbit (ISCO)"
    R_isco = 6.      # Orbital separation at ISCO, in geometric units. 6M for PN ISCO; 2.8M for EOB
    m1_min = 1
    m2_min = 1
Marc Arene's avatar
Marc Arene committed
117
    # approximant = 'TaylorF2'
Marc Arene's avatar
Marc Arene committed
118 119 120 121
    # approximant = 'SpinTaylorF2'
    # approximant = 'IMRPhenomD'
    # approximant = 'IMRPhenomPv2'
    # approximant = 'IMRPhenomPv2_NRTidal'
Marc Arene's avatar
Marc Arene committed
122 123 124 125
    f_ISCO = Compute_fISCO(R_isco, m1, m2)
    f_ISCO_run = Compute_fISCO(R_isco, m1_min, m2_min)
    # f_ISCO = Compute_LAL_fISCO(m1 * LAL_MSUN_SI, m2 * LAL_MSUN_SI)
    # f_ISCO_run = Compute_LAL_fISCO(m1_min * LAL_MSUN_SI, m2_min * LAL_MSUN_SI)
Marc Arene's avatar
Marc Arene committed
126
    if conf.approximant == 'TaylorF2':
Marc Arene's avatar
Marc Arene committed
127 128
        f_high = f_ISCO
        f_high_run = f_ISCO_run
129
        maximum_frequency_injected_waveform = 0
130
        maximum_frequency_search_waveform = 0
131 132
        # maximum_frequency_ifo = min(4096 * 4, f_high_run)
        maximum_frequency_ifo = f_high_run
Marc Arene's avatar
Marc Arene committed
133 134 135
        sampling_frequency = 2 * f_high_run
        duration = tc_3p5PN + 2
        reference_frequency = 0
Marc Arene's avatar
Marc Arene committed
136
    elif conf.approximant == 'IMRPhenomPv2' and conf.roq:
Marc Arene's avatar
Marc Arene committed
137 138 139
        # cf:
        # - https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom__c.html#gad3f98acfe9527259a7f73a8ef69a2f7b
        # - line 103 of https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/_l_a_l_sim_i_m_r_phenom_d_8c_source.html#l00103
Marc Arene's avatar
Marc Arene committed
140
        roq_directory = conf.roq_b_matrix_directory
Marc Arene's avatar
Marc Arene committed
141 142 143 144 145 146 147 148 149
        params_best_basis, scale_factor = roqu.get_best_segment_params_from_chirp_mass(injection_parameters['chirp_mass'], roq_directory=roq_directory)

        injection_parameters['roq'] = {}
        injection_parameters['roq']['directory'] = roq_directory + str(int(params_best_basis['seglen'])) + 's/'
        injection_parameters['roq']['params_path'] = injection_parameters['roq']['directory'] + 'params.dat'

        rescaled_params = roqu.rescale_params(params_best_basis, scale_factor)

        minimum_frequency = rescaled_params['flow']
Marc Arene's avatar
Marc Arene committed
150 151
        maximum_frequency_injected_waveform = rescaled_params['fhigh']
        # maximum_frequency_injected_waveform = min(2048, rescaled_params['fhigh'])
Marc Arene's avatar
Marc Arene committed
152
        maximum_frequency_ifo = min(conf.maximum_frequency, rescaled_params['fhigh'])
Marc Arene's avatar
Marc Arene committed
153 154
        maximum_frequency_search_waveform = maximum_frequency_injected_waveform
        # maximum_frequency_search_waveform = maximum_frequency_ifo
155

156
        # maximum_frequency_ifo = rescaled_params['fhigh']
Marc Arene's avatar
Marc Arene committed
157
        sampling_frequency = 2 * max(conf.maximum_frequency, rescaled_params['fhigh'])
Marc Arene's avatar
Marc Arene committed
158
        duration = rescaled_params['seglen']
Marc Arene's avatar
Marc Arene committed
159
        reference_frequency = conf.reference_frequency * scale_factor
Marc Arene's avatar
Marc Arene committed
160 161 162 163 164 165 166 167 168 169 170 171 172 173 174


        injection_parameters['roq']['params'] = params_best_basis
        injection_parameters['roq']['rescaled_params'] = rescaled_params
        injection_parameters['roq']['scale_factor'] = scale_factor

        tc_3p5PN = ComputeChirpTime3p5PN(minimum_frequency, m1, m2)

        LAL_MTSUN_SI = 4.925491025543575903411922162094833998e-6
        f_CUT = 0.2
        M_sec = (m1 + m2) * LAL_MTSUN_SI
        M_min_sec = (m1_min + m2_min) * LAL_MTSUN_SI
        f_high = f_CUT / M_sec
        f_high_run = f_CUT / M_min_sec
    else:
Marc Arene's avatar
Marc Arene committed
175 176 177 178 179 180 181 182 183
        # cf:
        # - https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom__c.html#gad3f98acfe9527259a7f73a8ef69a2f7b
        # - line 103 of https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/_l_a_l_sim_i_m_r_phenom_d_8c_source.html#l00103
        LAL_MTSUN_SI = 4.925491025543575903411922162094833998e-6
        f_CUT = 0.2
        M_sec = (m1 + m2) * LAL_MTSUN_SI
        M_min_sec = (m1_min + m2_min) * LAL_MTSUN_SI
        f_high = f_CUT / M_sec
        f_high_run = f_CUT / M_min_sec
184

Marc Arene's avatar
Marc Arene committed
185
        maximum_frequency_ifo = min(conf.maximum_frequency, f_high_run)
186
        maximum_frequency_injected_waveform = maximum_frequency_ifo
187 188
        maximum_frequency_search_waveform = maximum_frequency_ifo

189
        sampling_frequency = 2 * maximum_frequency_ifo
Marc Arene's avatar
Marc Arene committed
190
        duration = 64
Marc Arene's avatar
Marc Arene committed
191
        reference_frequency = conf.reference_frequency
192
    # import IPython; IPython.embed();
193 194
    # maximum_frequency_ifo = min(4096 * 4, f_high_run)
    # maximum_frequency_ifo = max(4096, f_high)
195 196 197
    # sampling_frequency = 4096
    # sampling_frequency = 4397
    start_time = injection_parameters['geocent_time'] - tc_3p5PN
Marc Arene's avatar
Marc Arene committed
198
    # breakpoint()
199 200
    meta_params = dict(
        minimum_frequency = minimum_frequency,
201
        maximum_frequency_injected_waveform = maximum_frequency_injected_waveform,
202
        maximum_frequency_search_waveform = maximum_frequency_search_waveform,
203
        maximum_frequency_ifo = maximum_frequency_ifo,
204
        maximum_frequency_generated_waveform = None,
Marc Arene's avatar
Marc Arene committed
205 206
        sampling_frequency = sampling_frequency,
        duration = duration,
207
        tc_3p5PN = tc_3p5PN,
Marc Arene's avatar
Marc Arene committed
208
        start_time = start_time,
209 210
        f_high = f_high,
        f_high_run = f_high_run,
Marc Arene's avatar
Marc Arene committed
211
        approximant = conf.approximant,
Marc Arene's avatar
Marc Arene committed
212
        reference_frequency = reference_frequency
213 214 215 216
    )

    injection_parameters['meta'] = meta_params

Marc Arene's avatar
Marc Arene committed
217 218 219 220 221 222
    conf.waveform_arguments = {}
    conf.waveform_arguments['minimum_frequency'] = minimum_frequency
    conf.waveform_arguments['waveform_approximant'] = conf.approximant
    conf.waveform_arguments['reference_frequency'] = reference_frequency
    conf.waveform_arguments['maximum_frequency'] = None

Marc Arene's avatar
Marc Arene committed
223
    conf.minimum_frequency = minimum_frequency
Marc Arene's avatar
Marc Arene committed
224 225
    conf.maximum_frequency_injected_waveform = maximum_frequency_injected_waveform
    conf.maximum_frequency_search_waveform = maximum_frequency_search_waveform
226
    conf.maximum_frequency_ifo = maximum_frequency_ifo
Marc Arene's avatar
Marc Arene committed
227
    conf.reference_frequency = reference_frequency
Marc Arene's avatar
Marc Arene committed
228 229 230 231 232 233 234 235
    conf.sampling_frequency = sampling_frequency

    conf.duration = duration
    conf.tc_3p5PN = tc_3p5PN
    conf.start_time = start_time
    conf.f_high = f_high
    conf.f_high_run = f_high_run

Marc Arene's avatar
Marc Arene committed
236 237

    # These two lines have been added to be able to plot the injected value on phi_c on the corner plots. Otherwise the `phi_c` and `tc_3p5PN` keys are not used in the code.
238
    injection_parameters['phi_c'] = (injection_parameters['phase'] * 2) % (2*np.pi)
Marc Arene's avatar
Marc Arene committed
239 240 241
    injection_parameters['tc_3p5PN'] = tc_3p5PN


242
    return injection_parameters
243 244 245 246 247




if __name__=='__main__':
Marc Arene's avatar
Marc Arene committed
248 249 250 251 252 253
    from optparse import OptionParser
    usage = """%prog [options]
    Plots the three waveforms in the time domain"""
    parser=OptionParser(usage)

    parser.add_option("--inj_file", default='../examples/GW170817.ini', action="store", type="string", help="""Injection file path leading to the `.ini` file specifying the parameters of the injection.""")
254

Marc Arene's avatar
Marc Arene committed
255 256 257
    (opts,args)=parser.parse_args()

    injection_parameters = ini_file_to_dict(opts.inj_file)
258 259

    pu.print_dict(injection_parameters, indent=3, align_keys=True)
Marc Arene's avatar
Marc Arene committed
260 261

    breakpoint()