set_injection_parameters.py 13.6 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
import Library.CONST as CONST
Marc Arene's avatar
Marc Arene committed
14 15
import Library.roq_utils as roqu

16 17
import core.utils as cut

Marc Arene's avatar
Marc Arene committed
18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41

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
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 78 79


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
80 81 82 83 84 85 86 87 88 89 90 91 92 93
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
94

Marc Arene's avatar
Marc Arene committed
95
def ini_file_to_dict(inj_file_path='../examples/injection_files/GW170817.ini'):
96
    config = ConfigParser()
97
    config.read(inj_file_path)
98 99 100
    injection_parameters = pu.config_parser_to_dict(config)['parameters']

    injection_parameters = conv.generate_mass_parameters(injection_parameters)
101 102
    injection_parameters['reduced_mass'] = paru.component_masses_to_reduced_mass(injection_parameters['mass_1'], injection_parameters['mass_2'])

Marc Arene's avatar
Marc Arene committed
103
    # 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.
104
    # injection_parameters['phi_c'] = (injection_parameters['phase'] * 2) % (2*np.pi)
Marc Arene's avatar
Marc Arene committed
105 106 107 108
    # injection_parameters['tc_3p5PN'] = tc_3p5PN
    return injection_parameters

def compute_extended_analysis_dict(mass_1, mass_2, geocent_time, chirp_mass, **analysis_kwargs):
109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
    """
    First, this function is clearly messy and needs some non-negligible re-writing... The idea is essentially to retrieve from the input parameters the GPS start_time of the analysis and the duration of the segment. Then some checks are performed, such as on the sampling frequency; the interferometers names are put in a list etc...

    Parameters
    ----------
    mass_1: float
        Mass, in solar mass units, of the first component of the binary.
    mass_2: float
        Mass, in solar mass units, of the second component of the binary.
    geocent_time: float
        GPS time, as measured at the center of the earth, of the merger.
    chirp_mass: float
        Chirp mass of the binary.
    analysis_kwargs: dictionary
        Additional options set in the '[analysis]' section of the config.ini file. The dictionary must contain the keys ['approximant', 'roq',  'minimum_frequency_ifos', 'ifos', 'reference_frequency', 'maximum_frequency',


    Returns
    -------
    ext_analysis_dict: dictionary
        Dictionary containing most importantly the duration of the segment analyzed and `tc_3p5PN` which is the duration to coalescence, and the GPS start_time.
    """
Marc Arene's avatar
Marc Arene committed
131 132 133 134

    approximant = analysis_kwargs['approximant']
    roq = analysis_kwargs['roq']
    roq_directory = analysis_kwargs['roq_b_matrix_directory']
135 136 137
    if type(analysis_kwargs['minimum_frequency_ifos']) == str:
        minimum_frequency_ifos = cut.convert_string_to_dict(analysis_kwargs['minimum_frequency_ifos'])
    else:
138
        minimum_frequency_ifos = {ifo: analysis_kwargs['minimum_frequency'] for ifo in analysis_kwargs['ifos']}
139
    minimum_frequency_all = min(minimum_frequency_ifos.values())
Marc Arene's avatar
Marc Arene committed
140 141 142 143
    reference_frequency = analysis_kwargs['reference_frequency']
    maximum_frequency = analysis_kwargs['maximum_frequency']

    ext_analysis_dict = {}
144 145 146 147 148
    ifo_chosen = analysis_kwargs['ifos'].split(',')
    if set.intersection(set(CONST.IFOS_POSSIBLE), set(ifo_chosen)) != set(ifo_chosen):
        raise ValueError("IFOs wrongly chosen: you must choose between {}. Example: '--ifos=H1,V1'. ".format(CONST.IFOS_POSSIBLE))
    else:
        ext_analysis_dict['ifos'] = ifo_chosen
Marc Arene's avatar
Marc Arene committed
149 150 151 152 153
    # approximant = 'TaylorF2'
    # approximant = 'SpinTaylorF2'
    # approximant = 'IMRPhenomD'
    # approximant = 'IMRPhenomPv2'
    # approximant = 'IMRPhenomPv2_NRTidal'
154 155 156 157

    # 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
158
    tc_3p5PN = ComputeChirpTime3p5PN(minimum_frequency_all, mass_1, mass_2)
159 160 161 162
    # 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
163
    f_ISCO = Compute_fISCO(R_isco, mass_1, mass_2)
Marc Arene's avatar
Marc Arene committed
164
    f_ISCO_run = Compute_fISCO(R_isco, m1_min, m2_min)
Marc Arene's avatar
Marc Arene committed
165
    duration = None
166
    t_end_minus_tc = 2
Marc Arene's avatar
Marc Arene committed
167
    # f_ISCO = Compute_LAL_fISCO(mass_1 * LAL_MSUN_SI, mass_2 * LAL_MSUN_SI)
Marc Arene's avatar
Marc Arene committed
168
    # f_ISCO_run = Compute_LAL_fISCO(m1_min * LAL_MSUN_SI, m2_min * LAL_MSUN_SI)
Marc Arene's avatar
Marc Arene committed
169
    if approximant == 'TaylorF2':
Marc Arene's avatar
Marc Arene committed
170 171
        f_high = f_ISCO
        f_high_run = f_ISCO_run
Marc Arene's avatar
Marc Arene committed
172 173 174 175 176
        # maximum_frequency_injected_waveform = 0
        # maximum_frequency_search_waveform = 0
        # # maximum_frequency_ifo = min(4096 * 4, f_high_run)
        # maximum_frequency_ifo = f_high_run
        # sampling_frequency = 2 * f_high_run
Marc Arene's avatar
Marc Arene committed
177
        # duration = int(tc_3p5PN) + 1 + 2
Marc Arene's avatar
Marc Arene committed
178
        # reference_frequency = 0
Marc Arene's avatar
Marc Arene committed
179 180 181 182 183 184
        maximum_frequency_ifo = min(maximum_frequency, f_high_run)
        maximum_frequency_injected_waveform = maximum_frequency_ifo
        maximum_frequency_search_waveform = maximum_frequency_ifo

        sampling_frequency = 2 * maximum_frequency_ifo
        # duration = 64
185
        duration = int(tc_3p5PN) + 1 + t_end_minus_tc
Marc Arene's avatar
Marc Arene committed
186
    elif approximant == 'IMRPhenomPv2' and roq:
Marc Arene's avatar
Marc Arene committed
187 188 189
        # 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
190
        params_best_basis, scale_factor = roqu.get_best_segment_params_from_chirp_mass(chirp_mass, roq_directory=roq_directory)
Marc Arene's avatar
Marc Arene committed
191

Marc Arene's avatar
Marc Arene committed
192 193 194
        ext_analysis_dict['roq'] = {}
        ext_analysis_dict['roq']['directory'] = roq_directory + str(int(params_best_basis['seglen'])) + 's/'
        ext_analysis_dict['roq']['params_path'] = ext_analysis_dict['roq']['directory'] + 'params.dat'
Marc Arene's avatar
Marc Arene committed
195 196 197

        rescaled_params = roqu.rescale_params(params_best_basis, scale_factor)

Marc Arene's avatar
Marc Arene committed
198 199
        # minimum_frequency = max(minimum_frequency, rescaled_params['flow'])
        # maximum_frequency_ifo = min(int(maximum_frequency), int(rescaled_params['fhigh']))
Marc Arene's avatar
Marc Arene committed
200
        # sampling_frequency = 2 * maximum_frequency_ifo
Marc Arene's avatar
Marc Arene committed
201
        # tc_3p5PN = ComputeChirpTime3p5PN(minimum_frequency, mass_1, mass_2)
Marc Arene's avatar
Marc Arene committed
202 203 204 205
        # # int() function because we need to make sure that sampling_freq * duration = integer to the 14 decimal precision
        # duration = min(int(tc_3p5PN + 2 + 1), int(rescaled_params['seglen']))
        # maximum_frequency_injected_waveform = rescaled_params['fhigh']
        # maximum_frequency_search_waveform = maximum_frequency_ifo
Marc Arene's avatar
Marc Arene committed
206
        # reference_frequency = reference_frequency * scale_factor
Marc Arene's avatar
Marc Arene committed
207 208


Marc Arene's avatar
Marc Arene committed
209
        # Settings which give me consistent results between the classic and the roq compute logL
210
        minimum_frequency_all = rescaled_params['flow']
Marc Arene's avatar
Marc Arene committed
211 212
        maximum_frequency_injected_waveform = rescaled_params['fhigh']
        # maximum_frequency_injected_waveform = min(2048, rescaled_params['fhigh'])
Marc Arene's avatar
Marc Arene committed
213
        maximum_frequency_ifo = min(maximum_frequency, rescaled_params['fhigh'])
Marc Arene's avatar
Marc Arene committed
214 215
        maximum_frequency_search_waveform = maximum_frequency_injected_waveform
        # maximum_frequency_search_waveform = maximum_frequency_ifo
216
        # maximum_frequency_ifo = rescaled_params['fhigh']
Marc Arene's avatar
Marc Arene committed
217
        sampling_frequency = 2 * max(maximum_frequency, rescaled_params['fhigh'])
Marc Arene's avatar
Marc Arene committed
218
        duration = rescaled_params['seglen']
Marc Arene's avatar
Marc Arene committed
219
        reference_frequency = reference_frequency * scale_factor
Marc Arene's avatar
Marc Arene committed
220 221


Marc Arene's avatar
Marc Arene committed
222 223 224



Marc Arene's avatar
Marc Arene committed
225 226 227
        ext_analysis_dict['roq']['params'] = params_best_basis
        ext_analysis_dict['roq']['rescaled_params'] = rescaled_params
        ext_analysis_dict['roq']['scale_factor'] = scale_factor
Marc Arene's avatar
Marc Arene committed
228

229
        tc_3p5PN = ComputeChirpTime3p5PN(minimum_frequency_all, mass_1, mass_2)
Marc Arene's avatar
Marc Arene committed
230 231 232

        LAL_MTSUN_SI = 4.925491025543575903411922162094833998e-6
        f_CUT = 0.2
Marc Arene's avatar
Marc Arene committed
233
        M_sec = (mass_1 + mass_2) * LAL_MTSUN_SI
Marc Arene's avatar
Marc Arene committed
234 235 236 237
        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
238 239 240 241 242
        # 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
Marc Arene's avatar
Marc Arene committed
243
        M_sec = (mass_1 + mass_2) * LAL_MTSUN_SI
Marc Arene's avatar
Marc Arene committed
244 245 246
        M_min_sec = (m1_min + m2_min) * LAL_MTSUN_SI
        f_high = f_CUT / M_sec
        f_high_run = f_CUT / M_min_sec
247

Marc Arene's avatar
Marc Arene committed
248
        maximum_frequency_ifo = min(maximum_frequency, f_high_run)
249
        maximum_frequency_injected_waveform = maximum_frequency_ifo
250 251
        maximum_frequency_search_waveform = maximum_frequency_ifo

252
        sampling_frequency = 2 * maximum_frequency_ifo
Marc Arene's avatar
Marc Arene committed
253
        # duration = 64
Marc Arene's avatar
Marc Arene committed
254 255
        # duration = int(tc_3p5PN) + 1 + t_end_minus_tc # formula for 59 sec
        duration = next_power_of_two(tc_3p5PN + t_end_minus_tc) # formula for 64 sec
256

Marc Arene's avatar
Marc Arene committed
257

Marc Arene's avatar
Marc Arene committed
258 259
    # Use an integer because frame files have an integer start_time
    # hence when extracting a subpart of it, the index will be exact
Marc Arene's avatar
Marc Arene committed
260 261
    start_time = geocent_time + t_end_minus_tc - duration # new formula which should be used
    # start_time = geocent_time  - tc_3p5PN # formula used for reruns of chap 8
Marc Arene's avatar
Marc Arene committed
262
    # start_time = int(geocent_time - tc_3p5PN)
Marc Arene's avatar
Marc Arene committed
263 264
    # if duration is None:
    #     duration = int(geocent_time - start_time) + 1 + 2
Marc Arene's avatar
Marc Arene committed
265 266

    ext_analysis_dict.update(
267 268
        minimum_frequency_ifos = minimum_frequency_ifos,
        minimum_frequency_all = minimum_frequency_all,
269
        maximum_frequency_injected_waveform = maximum_frequency_injected_waveform,
270
        maximum_frequency_search_waveform = maximum_frequency_search_waveform,
271
        maximum_frequency_ifo = maximum_frequency_ifo,
272
        maximum_frequency_generated_waveform = None,
Marc Arene's avatar
Marc Arene committed
273 274
        sampling_frequency = sampling_frequency,
        duration = duration,
275
        tc_3p5PN = tc_3p5PN,
Marc Arene's avatar
Marc Arene committed
276
        start_time = start_time,
277 278
        f_high = f_high,
        f_high_run = f_high_run,
Marc Arene's avatar
Marc Arene committed
279
        approximant = approximant,
Marc Arene's avatar
Marc Arene committed
280
        reference_frequency = reference_frequency
281 282
    )

Marc Arene's avatar
Marc Arene committed
283
    return ext_analysis_dict
284

285 286 287
def next_power_of_two(val):
    return 2**int(np.log(val)/np.log(2) + 1)

Marc Arene's avatar
Marc Arene committed
288
def get_inj_parameters_and_analysis_dict(inj_file_path='../examples/injection_files/GW170817.ini', **analysis_kwargs):
289 290 291 292 293 294 295 296 297 298 299

    injection_parameters = ini_file_to_dict(inj_file_path)

    ext_analysis_dict = compute_extended_analysis_dict(injection_parameters['mass_1'], injection_parameters['mass_2'], injection_parameters['geocent_time'], injection_parameters['chirp_mass'], **analysis_kwargs)

    # 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.
    injection_parameters['phi_c'] = (injection_parameters['phase'] * 2) % (2*np.pi)
    injection_parameters['tc_3p5PN'] = ext_analysis_dict['tc_3p5PN']

    return injection_parameters, ext_analysis_dict

300 301

if __name__=='__main__':
Marc Arene's avatar
Marc Arene committed
302 303 304 305 306
    from optparse import OptionParser
    usage = """%prog [options]
    Plots the three waveforms in the time domain"""
    parser=OptionParser(usage)

Marc Arene's avatar
Marc Arene committed
307
    parser.add_option("--inj_file", default='../examples/injection_files/GW170817.ini', action="store", type="string", help="""Injection file path leading to the `.ini` file specifying the parameters of the injection.""")
308

Marc Arene's avatar
Marc Arene committed
309 310 311
    (opts,args)=parser.parse_args()

    injection_parameters = ini_file_to_dict(opts.inj_file)
312 313

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

    breakpoint()