set_injection_parameters.py 12 KB
 Marc Arene committed Apr 23, 2019 1 2 3 ``````import sys # cf https://docs.python.org/3/tutorial/modules.html sys.path.append('../') `````` Marc Arene committed Jan 17, 2019 4 ``````import numpy as np `````` Marc Arene committed Apr 25, 2019 5 6 ``````from configparser import ConfigParser `````` Marc Arene committed Jan 17, 2019 7 ``````import bilby.gw.conversion as conv `````` Marc Arene committed Apr 25, 2019 8 `````` `````` Marc Arene committed Jan 17, 2019 9 ``````from Headers.PN_Coefficients import * # Headers.Constants already imported in PN_Coeff `````` Marc Arene committed Apr 23, 2019 10 ``````import Library.python_utils as pu `````` Marc Arene committed Jun 22, 2019 11 ``````import Library.param_utils as paru `````` Marc Arene committed Apr 23, 2019 12 `````` `````` Marc Arene committed Sep 27, 2019 13 ``````import Library.CONST as CONST `````` Marc Arene committed Jul 21, 2019 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.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 `````` Marc Arene committed Jan 17, 2019 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 committed Jul 21, 2019 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 `````` Marc Arene committed Jan 17, 2019 92 `````` `````` Marc Arene committed Dec 20, 2019 93 ``````def ini_file_to_dict(inj_file_path='../examples/injection_files/GW170817.ini'): `````` Marc Arene committed Apr 23, 2019 94 `````` config = ConfigParser() `````` Marc Arene committed Sep 26, 2019 95 `````` config.read(inj_file_path) `````` Marc Arene committed Apr 23, 2019 96 97 98 `````` injection_parameters = pu.config_parser_to_dict(config)['parameters'] injection_parameters = conv.generate_mass_parameters(injection_parameters) `````` Marc Arene committed Jun 22, 2019 99 100 `````` injection_parameters['reduced_mass'] = paru.component_masses_to_reduced_mass(injection_parameters['mass_1'], injection_parameters['mass_2']) `````` Marc Arene committed Sep 25, 2019 101 `````` # 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. `````` Marc Arene committed Sep 26, 2019 102 `````` # injection_parameters['phi_c'] = (injection_parameters['phase'] * 2) % (2*np.pi) `````` Marc Arene committed Sep 25, 2019 103 104 105 106 107 108 109 110 111 112 113 114 115 `````` # injection_parameters['tc_3p5PN'] = tc_3p5PN return injection_parameters def compute_extended_analysis_dict(mass_1, mass_2, geocent_time, chirp_mass, **analysis_kwargs): approximant = analysis_kwargs['approximant'] roq = analysis_kwargs['roq'] roq_directory = analysis_kwargs['roq_b_matrix_directory'] minimum_frequency = analysis_kwargs['minimum_frequency'] reference_frequency = analysis_kwargs['reference_frequency'] maximum_frequency = analysis_kwargs['maximum_frequency'] ext_analysis_dict = {} `````` Marc Arene committed Sep 27, 2019 116 117 118 119 120 `````` 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 committed Sep 25, 2019 121 122 123 124 125 `````` # approximant = 'TaylorF2' # approximant = 'SpinTaylorF2' # approximant = 'IMRPhenomD' # approximant = 'IMRPhenomPv2' # approximant = 'IMRPhenomPv2_NRTidal' `````` Marc Arene committed Apr 23, 2019 126 127 128 129 `````` # 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 `````` Marc Arene committed Sep 25, 2019 130 `````` tc_3p5PN = ComputeChirpTime3p5PN(minimum_frequency, mass_1, mass_2) `````` Marc Arene committed Apr 23, 2019 131 132 133 134 `````` # 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 committed Sep 25, 2019 135 `````` f_ISCO = Compute_fISCO(R_isco, mass_1, mass_2) `````` Marc Arene committed Jul 31, 2019 136 `````` f_ISCO_run = Compute_fISCO(R_isco, m1_min, m2_min) `````` Marc Arene committed Apr 30, 2020 137 `````` duration = None `````` Marc Arene committed Oct 18, 2020 138 `````` t_end_minus_tc = 2 `````` Marc Arene committed Sep 25, 2019 139 `````` # f_ISCO = Compute_LAL_fISCO(mass_1 * LAL_MSUN_SI, mass_2 * LAL_MSUN_SI) `````` Marc Arene committed Jul 31, 2019 140 `````` # f_ISCO_run = Compute_LAL_fISCO(m1_min * LAL_MSUN_SI, m2_min * LAL_MSUN_SI) `````` Marc Arene committed Sep 25, 2019 141 `````` if approximant == 'TaylorF2': `````` Marc Arene committed Jul 21, 2019 142 143 `````` f_high = f_ISCO f_high_run = f_ISCO_run `````` Marc Arene committed May 29, 2020 144 145 146 147 148 `````` # 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 committed Apr 30, 2020 149 `````` # duration = int(tc_3p5PN) + 1 + 2 `````` Marc Arene committed Sep 25, 2019 150 `````` # reference_frequency = 0 `````` Marc Arene committed May 29, 2020 151 152 153 154 155 156 `````` 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 `````` Marc Arene committed Oct 18, 2020 157 `````` duration = int(tc_3p5PN) + 1 + t_end_minus_tc `````` Marc Arene committed Sep 25, 2019 158 `````` elif approximant == 'IMRPhenomPv2' and roq: `````` Marc Arene committed Jul 21, 2019 159 160 161 `````` # 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 committed Sep 25, 2019 162 `````` params_best_basis, scale_factor = roqu.get_best_segment_params_from_chirp_mass(chirp_mass, roq_directory=roq_directory) `````` Marc Arene committed Jul 21, 2019 163 `````` `````` Marc Arene committed Sep 25, 2019 164 165 166 `````` 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 committed Jul 21, 2019 167 168 169 `````` rescaled_params = roqu.rescale_params(params_best_basis, scale_factor) `````` Marc Arene committed Sep 25, 2019 170 171 `````` # minimum_frequency = max(minimum_frequency, rescaled_params['flow']) # maximum_frequency_ifo = min(int(maximum_frequency), int(rescaled_params['fhigh'])) `````` Marc Arene committed Sep 25, 2019 172 `````` # sampling_frequency = 2 * maximum_frequency_ifo `````` Marc Arene committed Sep 25, 2019 173 `````` # tc_3p5PN = ComputeChirpTime3p5PN(minimum_frequency, mass_1, mass_2) `````` Marc Arene committed Sep 25, 2019 174 175 176 177 `````` # # 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 committed Sep 25, 2019 178 `````` # reference_frequency = reference_frequency * scale_factor `````` Marc Arene committed Sep 25, 2019 179 180 `````` `````` Marc Arene committed Sep 25, 2019 181 `````` # Settings which give me consistent results between the classic and the roq compute logL `````` Marc Arene committed Jul 21, 2019 182 `````` minimum_frequency = rescaled_params['flow'] `````` Marc Arene committed Jul 31, 2019 183 184 `````` maximum_frequency_injected_waveform = rescaled_params['fhigh'] # maximum_frequency_injected_waveform = min(2048, rescaled_params['fhigh']) `````` Marc Arene committed Sep 25, 2019 185 `````` maximum_frequency_ifo = min(maximum_frequency, rescaled_params['fhigh']) `````` Marc Arene committed Jul 31, 2019 186 187 `````` maximum_frequency_search_waveform = maximum_frequency_injected_waveform # maximum_frequency_search_waveform = maximum_frequency_ifo `````` Marc Arene committed Jul 22, 2019 188 `````` # maximum_frequency_ifo = rescaled_params['fhigh'] `````` Marc Arene committed Sep 25, 2019 189 `````` sampling_frequency = 2 * max(maximum_frequency, rescaled_params['fhigh']) `````` Marc Arene committed Jul 21, 2019 190 `````` duration = rescaled_params['seglen'] `````` Marc Arene committed Sep 25, 2019 191 `````` reference_frequency = reference_frequency * scale_factor `````` Marc Arene committed Jul 21, 2019 192 193 `````` `````` Marc Arene committed Sep 25, 2019 194 195 196 `````` `````` Marc Arene committed Sep 25, 2019 197 198 199 `````` 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 committed Jul 21, 2019 200 `````` `````` Marc Arene committed Sep 25, 2019 201 `````` tc_3p5PN = ComputeChirpTime3p5PN(minimum_frequency, mass_1, mass_2) `````` Marc Arene committed Jul 21, 2019 202 203 204 `````` LAL_MTSUN_SI = 4.925491025543575903411922162094833998e-6 f_CUT = 0.2 `````` Marc Arene committed Sep 25, 2019 205 `````` M_sec = (mass_1 + mass_2) * LAL_MTSUN_SI `````` Marc Arene committed Jul 21, 2019 206 207 208 209 `````` 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 committed Jul 11, 2019 210 211 212 213 214 `````` # 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 committed Sep 25, 2019 215 `````` M_sec = (mass_1 + mass_2) * LAL_MTSUN_SI `````` Marc Arene committed Jul 11, 2019 216 217 218 `````` M_min_sec = (m1_min + m2_min) * LAL_MTSUN_SI f_high = f_CUT / M_sec f_high_run = f_CUT / M_min_sec `````` Marc Arene committed Jul 22, 2019 219 `````` `````` Marc Arene committed Sep 25, 2019 220 `````` maximum_frequency_ifo = min(maximum_frequency, f_high_run) `````` Marc Arene committed Jul 22, 2019 221 `````` maximum_frequency_injected_waveform = maximum_frequency_ifo `````` Marc Arene committed Jul 22, 2019 222 223 `````` maximum_frequency_search_waveform = maximum_frequency_ifo `````` Marc Arene committed Jul 22, 2019 224 `````` sampling_frequency = 2 * maximum_frequency_ifo `````` Marc Arene committed Dec 20, 2019 225 `````` # duration = 64 `````` Marc Arene committed Jan 17, 2021 226 227 `````` 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 `````` Marc Arene committed Oct 18, 2020 228 `````` `````` Marc Arene committed Sep 25, 2019 229 `````` `````` Marc Arene committed Dec 20, 2019 230 231 `````` # 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 committed Jan 17, 2021 232 233 `````` # 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 committed Dec 20, 2019 234 `````` # start_time = int(geocent_time - tc_3p5PN) `````` Marc Arene committed Apr 30, 2020 235 236 `````` # if duration is None: # duration = int(geocent_time - start_time) + 1 + 2 `````` Marc Arene committed Sep 25, 2019 237 238 `````` ext_analysis_dict.update( `````` Marc Arene committed Apr 23, 2019 239 `````` minimum_frequency = minimum_frequency, `````` Marc Arene committed Jul 22, 2019 240 `````` maximum_frequency_injected_waveform = maximum_frequency_injected_waveform, `````` Marc Arene committed Jul 22, 2019 241 `````` maximum_frequency_search_waveform = maximum_frequency_search_waveform, `````` Marc Arene committed Jul 22, 2019 242 `````` maximum_frequency_ifo = maximum_frequency_ifo, `````` Marc Arene committed Jul 22, 2019 243 `````` maximum_frequency_generated_waveform = None, `````` Marc Arene committed Jul 21, 2019 244 245 `````` sampling_frequency = sampling_frequency, duration = duration, `````` Marc Arene committed Apr 23, 2019 246 `````` tc_3p5PN = tc_3p5PN, `````` Marc Arene committed Jul 21, 2019 247 `````` start_time = start_time, `````` Marc Arene committed Apr 23, 2019 248 249 `````` f_high = f_high, f_high_run = f_high_run, `````` Marc Arene committed Sep 25, 2019 250 `````` approximant = approximant, `````` Marc Arene committed Jul 21, 2019 251 `````` reference_frequency = reference_frequency `````` Marc Arene committed Apr 23, 2019 252 253 `````` ) `````` Marc Arene committed Sep 25, 2019 254 `````` return ext_analysis_dict `````` Marc Arene committed Apr 24, 2019 255 `````` `````` Marc Arene committed Oct 18, 2020 256 257 258 ``````def next_power_of_two(val): return 2**int(np.log(val)/np.log(2) + 1) `````` Marc Arene committed Dec 20, 2019 259 ``````def get_inj_parameters_and_analysis_dict(inj_file_path='../examples/injection_files/GW170817.ini', **analysis_kwargs): `````` Marc Arene committed Sep 26, 2019 260 261 262 263 264 265 266 267 268 269 270 `````` 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 `````` Marc Arene committed Apr 24, 2019 271 272 `````` if __name__=='__main__': `````` Marc Arene committed Jul 21, 2019 273 274 275 276 277 `````` from optparse import OptionParser usage = """%prog [options] Plots the three waveforms in the time domain""" parser=OptionParser(usage) `````` Marc Arene committed Dec 20, 2019 278 `````` 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.""") `````` Marc Arene committed Apr 24, 2019 279 `````` `````` Marc Arene committed Jul 21, 2019 280 281 282 `````` (opts,args)=parser.parse_args() injection_parameters = ini_file_to_dict(opts.inj_file) `````` Marc Arene committed Apr 24, 2019 283 284 `````` pu.print_dict(injection_parameters, indent=3, align_keys=True) `````` Marc Arene committed Jul 21, 2019 285 286 `````` breakpoint()``````