Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALInference 4.1.9.1-00ddc7f
•All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
imrtgr_imr_consistency_test.py
Go to the documentation of this file.
1##python
2"""
3Perform the consistency check between the inspiral and ringdown estimates of the mass and spin of the final
4black hole in a binary black hole merger.
5
6P. Ajith, Abhirup Ghosh, Archisman Ghosh, 2015-09-18
7
8$Id:$
9"""
10
11import matplotlib as mpl
12mpl.use('Agg')
13import matplotlib.pyplot as plt
14import numpy as np
15import scipy
16from scipy import interpolate
17from optparse import OptionParser
18import time, os
19import lalinference.imrtgr.imrtgrutils as tgr
20import pickle, gzip
21import sys
22from lalinference import git_version
23
24from scipy.stats import gaussian_kde #rahul: for KDE implementation
25
26from matplotlib import rc
27import matplotlib
28matplotlib.rc('text.latex', preamble=r'\usepackage{txfonts}')
29
30rc('text', usetex=True)
31rc('font', family='serif')
32rc('font', serif='times')
33rc('mathtext', default='sf')
34rc("lines", markeredgewidth=1)
35rc("lines", linewidth=2)
36rc('axes', labelsize=10)
37rc("axes", linewidth=0.5)
38rc('xtick', labelsize=8)
39rc('ytick', labelsize=8)
40rc('legend', fontsize=10)
41rc('xtick.major', pad=6)
42rc('ytick.major', pad=6)
43rc('xtick.minor', size=5)
44rc('ytick.minor', size=5)
45
46def set_tick_sizes(ax, major, minor):
47 for l in ax.get_xticklines() + ax.get_yticklines():
48 l.set_markersize(major)
49 for tick in ax.xaxis.get_minor_ticks() + ax.yaxis.get_minor_ticks():
50 tick.tick1line.set_markersize(minor)
51 tick.tick2line.set_markersize(minor)
52 ax.xaxis.LABELPAD=10.
53 ax.xaxis.OFFSETTEXTPAD=10.
54
55# Module for confidence calculations
56class confidence(object):
57 def __init__(self, counts):
58 # Sort in descending order in frequency
59 self.counts_sorted = np.sort(counts.flatten())[::-1]
60 # Get a normalized cumulative distribution from the mode
61 self.norm_cumsum_counts_sorted = np.cumsum(self.counts_sorted) / np.sum(counts)
62 # Set interpolations between heights, bins and levels
63 self._set_interp()
64 def _set_interp(self):
65 self._length = len(self.counts_sorted)
66 # height from index
67 self._height_from_idx = interpolate.interp1d(np.arange(self._length), self.counts_sorted, bounds_error=False, fill_value=0.)
68 # index from height
69 self._idx_from_height = interpolate.interp1d(self.counts_sorted[::-1], np.arange(self._length)[::-1], bounds_error=False, fill_value=self._length)
70 # level from index
71 self._level_from_idx = interpolate.interp1d(np.arange(self._length), self.norm_cumsum_counts_sorted, bounds_error=False, fill_value=1.)
72 # index from level
73 self._idx_from_level = interpolate.interp1d(self.norm_cumsum_counts_sorted, np.arange(self._length), bounds_error=False, fill_value=self._length)
74 def level_from_height(self, height):
75 return self._level_from_idx(self._idx_from_height(height))
76 def height_from_level(self, level):
77 return self._height_from_idx(self._idx_from_level(level))
78
79######################################################################################
80################################# MAIN PROGRAM #######################################
81######################################################################################
82
83if __name__ == '__main__':
84
85 start_time = time.time()
86
87 # read inputs from command line
88 parser = OptionParser()
89 parser.add_option("-i", "--insp-post", dest="insp_post", help="file containing the posterior samples from the PE inspiral run")
90 parser.add_option("-r", "--ring-post", dest="ring_post", help="file containing the posterior samples from the PE ringdown run")
91 parser.add_option("-m", "--imr-post", dest="imr_post", help="file containing the posterior samples from the full PE IMR run")
92 parser.add_option("-f", "--fit-formula", dest="fit_formula", help="fitting formula to be used for the calculation of final mass/spin [options: 'nospin_Pan2011', 'nonprecspin_Healy2014', 'bbh_average_fits_precessing'", default="nonprecspin_Healy2014")
93 parser.add_option("-p", "--mf-chif-prior", dest="prior_Mfchif_file", help="pickle file containing the interpolation object of the prior in (Mf, chif) used in the PE runs", default=None)
94 parser.add_option("-o", "--out-dir", dest="out_dir", help="output directory")
95 parser.add_option("--insp-fhigh", dest="insp_fhigh", help="Upper cutoff freq for the inspiral analysis")
96 parser.add_option("--ring-flow", dest="ring_flow", help="Lower cutoff freq for the ringdown analysis")
97 parser.add_option("--m1-inj", dest="m1_inj", help="injected value of component mass m1 (if this is an injection)")
98 parser.add_option("--m2-inj", dest="m2_inj", help="injected value of component mass m2 (if this is an injection)")
99 parser.add_option("--chi1-inj", dest="chi1_inj", help="injected value of spin magnitude of mass m1 (if this is an injection)")
100 parser.add_option("--chi2-inj", dest="chi2_inj", help="injected value of spin magnitude of mass m2 (if this is an injection)")
101 parser.add_option("--chi1z-inj", dest="chi1z_inj", help="injected value of z-component of spin of mass m1 (if this is an injection)")
102 parser.add_option("--chi2z-inj", dest="chi2z_inj", help="injected value of z-component of spin of mass m2 (if this is an injection)")
103 parser.add_option("--phi12-inj", dest="phi12_inj", help="injected value of the azimuth angle of (hats2 - hats1) from vecL (if this is an injection)", default=0.)
104 parser.add_option("-w", "--waveform", dest="waveform", help="waveform used for recovery")
105 parser.add_option("-d", "--debug-plots", dest="debug_plots", help="debug plots")
106 parser.add_option("--N_bins", type="int", dest="N_bins", default=201, help="number of bins (default=201)")
107 parser.add_option("--dMfbyMf_lim", type="float", dest="dMfbyMf_lim", default=1., help="absolute value of limit for range of dMfbyMf_vec, defined as [-dMfbyMf_lim, +dMfbyMf_lim]")
108 parser.add_option("--dchifbychif_lim", type="float", dest="dchifbychif_lim", default=1., help="absolute value of limit for range of dchifbychif_vec, defined as [-dchifbychif_lim, +dchifbychif_lim]")
109 parser.add_option("--use_KDE", type="int", dest="MfafKDE", help="use KDE or not after getting samples of Mf, af")
110 parser.add_option("-s", "--pepostproc", dest="pepostproc", default= "cbcbayespostproc", help="package used to generate summary webpage to visualize the output [options: 'cbcbayespostproc', 'pesummary']")
111
112 (options, args) = parser.parse_args()
113 MfafKDE = options.MfafKDE
114 insp_post = options.insp_post
115 ring_post = options.ring_post
116 imr_post = options.imr_post
117 prior_Mfchif_file = options.prior_Mfchif_file
118 out_dir = options.out_dir
119 fit_formula = options.fit_formula
120 debug_plots = options.debug_plots
121 pepostproc = options.pepostproc
122 if options.insp_fhigh is not None:
123 insp_fhigh = float(options.insp_fhigh)
124 else:
125 print('Inspiral cutoff freq not provided. To have it displayed on the results page, please pass command line option --insp-fhigh.')
126 insp_fhigh = np.nan
127 if options.ring_flow is not None:
128 ring_flow = float(options.ring_flow)
129 else:
130 print('Ringdown cutoff freq not provided. To have it displayed on the results page, please pass command line option --ring-flow.')
131 ring_flow = np.nan
132 waveform = options.waveform
133 if waveform is None:
134 print('Recovery approximant not provided. To have it displayed on the results page, please pass command line option --waveform.')
135
136 N_bins = int(options.N_bins) # Number of grid points along either axis (dMfbyMf, dchifbychif) for computation of the posteriors
137 dMfbyMf_lim = float(options.dMfbyMf_lim)
138 dchifbychif_lim = float(options.dchifbychif_lim)
139 lalinference_datadir = os.getenv('LALINFERENCE_DATADIR')
140
141 # create output directory and copy the script
142 os.system('mkdir -p %s' %out_dir)
143 os.system('mkdir -p %s/data' %out_dir)
144 os.system('mkdir -p %s/img' %out_dir)
145 os.system('cp %s %s' %(__file__, out_dir))
146 os.system('cp %s %s/'%(os.path.join(lalinference_datadir, 'imrtgr_webpage_templates/*.*'), out_dir))
147
148 # creating file to save the run command
149 run_command = open('%s/command.txt'%(out_dir),'w')
150 for arg in sys.argv:
151 run_command.write('%s\n' %arg)
152 run_command.write("\n")
153 run_command.write("\n")
154 run_command.write("%s"%git_version.verbose_msg)
155 run_command.close()
156
157 # creating soft links for PE results
158 if pepostproc == 'cbcbayespostproc':
159 insp_post_link = insp_post
160 ring_post_link = ring_post
161 imr_post_link = imr_post
162 # posterior samples data file generated by pesummary is stored in
163 # <path to summary directory>/samples/<posterior_samples>.dat
164 # need to go to the parent directory before creating soft link
165 elif pepostproc == 'pesummary':
166 insp_post_link = os.path.realpath(os.path.dirname(insp_post))
167 ring_post_link = os.path.realpath(os.path.dirname(ring_post))
168 imr_post_link = os.path.realpath(os.path.dirname(imr_post))
169 os.system("sed -i 's/posplots/home/g' %s/result.html"%out_dir)
170 else:
171 raise ValueError('Package to generate summary webpage unknown !')
172
173 insp_posplots = os.path.realpath(os.path.dirname(insp_post_link))
174 ring_posplots = os.path.realpath(os.path.dirname(ring_post_link))
175 imr_posplots = os.path.realpath(os.path.dirname(imr_post_link))
176
177 insp_target = os.path.join(out_dir, 'pe_insp')
178 ring_target = os.path.join(out_dir, 'pe_ring')
179 imr_target = os.path.join(out_dir, 'pe_imr')
180
181 if insp_posplots != insp_target:
182 if os.path.islink(insp_target):
183 print('... removing existing link %s'%(insp_target))
184 os.system('rm %s'%(insp_target))
185 print('... linking %s to %s' %(insp_posplots, insp_target))
186 os.system('ln -s %s %s' %(insp_posplots, insp_target))
187 if ring_posplots != ring_target:
188 if os.path.islink(ring_target):
189 print('... removing existing link %s'%(ring_target))
190 os.system('rm %s'%(ring_target))
191 print('... linking %s to %s' %(ring_posplots, ring_target))
192 os.system('ln -s %s %s' %(ring_posplots, ring_target))
193 if imr_posplots != imr_target:
194 if os.path.islink(imr_target):
195 print('... removing existing link %s'%(imr_target))
196 os.system('rm %s'%(imr_target))
197 print('... linking %s to %s' %(imr_posplots, imr_target))
198 os.system('ln -s %s %s' %(imr_posplots, imr_target))
199
200 # read the injection mass parameters if this is an injection
201 m1_inj = options.m1_inj
202 m2_inj = options.m2_inj
203 chi1_inj = options.chi1_inj
204 chi2_inj = options.chi2_inj
205 chi1z_inj = options.chi1z_inj
206 chi2z_inj = options.chi2z_inj
207 phi12_inj = options.phi12_inj
208
209 if m1_inj == None or m2_inj == None or chi1_inj == None or chi2_inj == None or chi1z_inj == None or chi2z_inj == None or phi12_inj == None:
210 plot_injection_lines = False
211 else:
212 m1_inj = float(m1_inj)
213 m2_inj = float(m2_inj)
214 chi1_inj = float(chi1_inj)
215 chi2_inj = float(chi2_inj)
216 chi1z_inj = float(chi1z_inj)
217 chi2z_inj = float(chi2z_inj)
218 phi12_inj = float(phi12_inj)
219 plot_injection_lines = True
220 q_inj = m1_inj/m2_inj
221 eta_inj = q_inj/(1.+q_inj)**2.
222 Mf_inj, chif_inj = tgr.calc_final_mass_spin(m1_inj, m2_inj, chi1_inj, chi2_inj, chi1z_inj, chi2z_inj, phi12_inj, fit_formula)
223
224 ###############################################################################################
225 # Read the posteriors from the inspiral, ringdown and imr PE runs (after post-processing)
226 ###############################################################################################
227
228 # parameter names assigned according to the package used to generate summary webpage
229 if pepostproc == 'cbcbayespostproc':
230 m1_id, m2_id, a1_id, a2_id, a1z_id, a2z_id, phi12_id = 'm1', 'm2', 'a1', 'a2', 'a1z', 'a2z', 'phi12'
231 elif pepostproc == 'pesummary':
232 m1_id, m2_id, a1_id, a2_id, a1z_id, a2z_id, phi12_id = 'mass_1', 'mass_2', 'a_1', 'a_2', 'spin_1z', 'spin_2z', 'phi_12'
233 else:
234 raise ValueError('Package to generate summary webpage unknown !')
235
236 # read data from the inspiral posterior file
237 insp_data = np.genfromtxt(insp_post, dtype=None, names=True)
238
239 m1_i, m2_i, chi1_i, chi2_i, chi1z_i, chi2z_i = insp_data[m1_id], insp_data[m2_id], insp_data[a1_id], insp_data[a2_id], insp_data[a1z_id], insp_data[a2z_id]
240 # if there is phi12 in the posterior, read the values.
241 if phi12_id in insp_data.dtype.names:
242 phi12_i = insp_data[phi12_id]
243 else:
244 phi12_i = np.zeros(len(m1_i))
245 # compute the final mass and spin
246 Mf_i, chif_i = tgr.calc_final_mass_spin(m1_i, m2_i, chi1_i, chi2_i, chi1z_i, chi2z_i, phi12_i, fit_formula)
247
248 # read data from the ringdown posterior file
249 ring_data = np.genfromtxt(ring_post, dtype=None, names=True)
250
251 m1_r, m2_r, chi1_r, chi2_r, chi1z_r, chi2z_r = ring_data[m1_id], ring_data[m2_id], ring_data[a1_id], ring_data[a2_id], ring_data[a1z_id], ring_data[a2z_id]
252 # if there is phi12 in the posterior, read the values.
253 if phi12_id in ring_data.dtype.names:
254 phi12_r = ring_data[phi12_id]
255 else:
256 phi12_r = np.zeros(len(m1_r))
257 # compute the final mass and spin
258 Mf_r, chif_r = tgr.calc_final_mass_spin(m1_r, m2_r, chi1_r, chi2_r, chi1z_r, chi2z_r, phi12_r, fit_formula)
259
260 # read data from the IMR posterior file
261 imr_data = np.genfromtxt(imr_post, dtype=None, names=True)
262
263 m1_imr, m2_imr, chi1_imr, chi2_imr, chi1z_imr, chi2z_imr = imr_data[m1_id], imr_data[m2_id], imr_data[a1_id], imr_data[a2_id], imr_data[a1z_id], imr_data[a2z_id]
264 # if there is phi12 in the posterior, read the values.
265 if phi12_id in imr_data.dtype.names:
266 phi12_imr = imr_data[phi12_id]
267 else:
268 phi12_imr = np.zeros(len(m1_imr))
269 # compute the final mass and spin
270 Mf_imr, chif_imr = tgr.calc_final_mass_spin(m1_imr, m2_imr, chi1_imr, chi2_imr, chi1z_imr, chi2z_imr, phi12_imr, fit_formula)
271
272
273 print('... read posteriors')
274 ###############################################################################################
275
276 ###############################################################################################
277 # compute the limits of integration for computing delta_Mf and delta_chif
278 ###############################################################################################
279 Mf_lim = max(abs(np.append(np.append(Mf_i, Mf_r), Mf_imr)))
280 chif_lim = max(abs(np.append(np.append(chif_i, chif_r), chif_imr)))
281
282 # the integral used to compute (Delta Mf, Delta af) has limits from -infinity to +infinity. We
283 # are approximating this by setting the limits to (-Mf_lim to Mf_lim) and (-chif_lim to chif_lim)
284 # where Mf_lim and chif_lim are the max values of Mf and chif where the posteriors have nonzero
285 # support. The scipy.signal.correlate2d function requires arguments x_bins and y_bins that need
286 # to be symmetric around zero
287 Mf_bins = np.linspace(-Mf_lim, Mf_lim, N_bins)
288 chif_bins = np.linspace(-chif_lim, chif_lim, N_bins)
289
290 dMf = np.mean(np.diff(Mf_bins))
291 dchif = np.mean(np.diff(chif_bins))
292
293 Mf_intp = (Mf_bins[:-1] + Mf_bins[1:])/2.
294 chif_intp = (chif_bins[:-1] + chif_bins[1:])/2.
295
296
297 print('useKDE=',MfafKDE)
298 if MfafKDE==1:
299 print('replacing lal P(Mfaf) with its KDE pdf')
300 M_i,C_i=np.meshgrid(Mf_intp,chif_intp)
301
302 joint_data=np.vstack([Mf_i,chif_i]);kernel=gaussian_kde(joint_data)
303 f_i = lambda x,y:kernel.evaluate([x,y])
304 print("for inspiral kernel",kernel.integrate_box([-Mf_lim,-chif_lim],[Mf_lim,chif_lim]))
305 P_Mfchif_i = np.vectorize(f_i)(M_i,C_i)/kernel.integrate_box([-Mf_lim,-chif_lim],[Mf_lim,chif_lim])
306
307 joint_data=np.vstack([Mf_r,chif_r]);kernel=gaussian_kde(joint_data)#;M_i,C_i=np.meshgrid(Mf_bins,chif_bins)
308 f_r = lambda x,y:kernel.evaluate([x,y])
309 print("for post-inspiral kernel",kernel.integrate_box([-Mf_lim,-chif_lim],[Mf_lim,chif_lim]))
310 P_Mfchif_r = np.vectorize(f_r)(M_i,C_i)/kernel.integrate_box([-Mf_lim,-chif_lim],[Mf_lim,chif_lim])
311
312 joint_data=np.vstack([Mf_imr,chif_imr]);kernel=gaussian_kde(joint_data)#;M_i,C_i=np.meshgrid(Mf_bins,chif_bins)
313 f_imr = lambda x,y:kernel.evaluate([x,y])
314 print("for imr kernel",kernel.integrate_box([-Mf_lim,-chif_lim],[Mf_lim,chif_lim]))
315 P_Mfchif_imr = np.vectorize(f_imr)(M_i,C_i)/kernel.integrate_box([-Mf_lim,-chif_lim],[Mf_lim,chif_lim])
316
317
318 #rahul: end KDE of Mf,af
319
320 elif MfafKDE==0:
321 print('using default samples, NOKDE')
322 # compute the 2D posterior distributions for the inspiral, ringodwn and IMR analyses
323 P_Mfchif_i, Mf_bins, chif_bins = np.histogram2d(Mf_i, chif_i, bins=(Mf_bins, chif_bins), density=True)
324 P_Mfchif_r, Mf_bins, chif_bins = np.histogram2d(Mf_r, chif_r, bins=(Mf_bins, chif_bins), density=True)
325 P_Mfchif_imr, Mf_bins, chif_bins = np.histogram2d(Mf_imr, chif_imr, bins=(Mf_bins, chif_bins), density=True)
326 # transpose to go from (X,Y) indexing returned by np.histogram2d() to array (i,j) indexing for further
327 # computations. From now onwards, different rows (i) correspond to different values of Mf and different
328 # columns (j) correspond to different values of chif
329 #rahul:Transpose is forbidden if using KDE
330 P_Mfchif_i = P_Mfchif_i.T
331 P_Mfchif_r = P_Mfchif_r.T
332 P_Mfchif_imr = P_Mfchif_imr.T
333 print('computed P_Mfchif without using KDE')
334
335
336 ###############################################################################################
337
338
339 ###############################################################################################
340 # Undo the effect of the prior from the lalinference posterior. Lalinference assumes a #
341 # uniform prior in component masses. We need to assume a uniform prior in Mf, chif #
342 ###############################################################################################
343
344 if prior_Mfchif_file is not None:
345
346 os.system('cp %s %s/data'%(prior_Mfchif_file, out_dir))
347
348 # read the interpolation object, reconstruct the data from the interpolation object
349 f = gzip.open(prior_Mfchif_file,'rb')
350 P_Mfchif_pr_interp_obj = pickle.load(f)
351 P_Mfchif_pr = P_Mfchif_pr_interp_obj(Mf_intp, chif_intp)
352
353 # compute the corrected 2D posteriors in Mf and chif by dividing by the prior distribution
354 P_Mfchif_i = P_Mfchif_i/P_Mfchif_pr
355 P_Mfchif_r = P_Mfchif_r/P_Mfchif_pr
356 P_Mfchif_imr = P_Mfchif_imr/P_Mfchif_pr
357
358 # removing nan's
359 P_Mfchif_i[np.isnan(P_Mfchif_i)] = 0.
360 P_Mfchif_r[np.isnan(P_Mfchif_r)] = 0.
361 P_Mfchif_imr[np.isnan(P_Mfchif_imr)] = 0.
362
363 # removing infinities
364 P_Mfchif_i[np.isinf(P_Mfchif_i)] = 0.
365 P_Mfchif_r[np.isinf(P_Mfchif_r)] = 0.
366 P_Mfchif_imr[np.isinf(P_Mfchif_imr)] = 0.
367
368 print('... computed (prior) corrected posteriors')
369
370 ################################################################################################
371 # compute the posterior of (delta_Mf/Mf, delta_chif/chif)
372 ################################################################################################
373
374 # compute interpolation objects for the Mf,chif posterior and delta_Mf and delta_chif posterior
375 P_Mfchif_i_interp_object = scipy.interpolate.interp2d(Mf_intp, chif_intp, P_Mfchif_i, fill_value=0., bounds_error=False)
376 P_Mfchif_r_interp_object = scipy.interpolate.interp2d(Mf_intp, chif_intp, P_Mfchif_r, fill_value=0., bounds_error=False)
377
378 # defining limits of delta_Mf/Mf and delta_chif/chif.
379 dMfbyMf_vec = np.linspace(-dMfbyMf_lim, dMfbyMf_lim, N_bins)
380 dchifbychif_vec = np.linspace(-dchifbychif_lim, dchifbychif_lim, N_bins)
381
382 # compute the P(dMf/Mf, dchif/chif) by evaluating the integral
383 diff_dMfbyMf = np.mean(np.diff(dMfbyMf_vec))
384 diff_dchifbychif = np.mean(np.diff(dchifbychif_vec))
385 P_dMfbyMf_dchifbychif = np.zeros(shape=(N_bins,N_bins))
386
387 # compute the posterior on the fractional deviation parameters (delta_Mf/Mf, delta_chif/chif).
388 # Approximate the integral in Eq.(6) of the document LIGO-P1500185-v5 by a discrete sum
389 for i, v2 in enumerate(dchifbychif_vec):
390 for j, v1 in enumerate(dMfbyMf_vec):
391 P_dMfbyMf_dchifbychif[i,j] = tgr.calc_sum(Mf_intp, chif_intp, v1, v2, P_Mfchif_i_interp_object, P_Mfchif_r_interp_object)
392
393 # normalization
394 P_dMfbyMf_dchifbychif /= np.sum(P_dMfbyMf_dchifbychif) * diff_dMfbyMf * diff_dchifbychif
395
396 # Marginalization to one-dimensional joint_posteriors
397 P_dMfbyMf = np.sum(P_dMfbyMf_dchifbychif, axis=0) * diff_dchifbychif
398 P_dchifbychif = np.sum(P_dMfbyMf_dchifbychif, axis=1) * diff_dMfbyMf
399
400 # compute the confidence region corresponding to the GR value (delta_Mf/Mf = 0, delta_chif/chif = 0).
401 # the 'confidence' class is defined on top of this script
402 conf_v1v2 = confidence(P_dMfbyMf_dchifbychif)
403 gr_height = P_dMfbyMf_dchifbychif[np.argmin(abs(dMfbyMf_vec)), np.argmin(abs(dchifbychif_vec))] # taking value closest to (0,0)
404 gr_conf_level = conf_v1v2.level_from_height(gr_height)
405 print('... no deviation from GR above %.1f%% confidence level'%(100.*gr_conf_level))
406
407 # creating the parameter table
408 param_table = [['Upper cutoff freq for the inspiral analysis: %s Hz'%insp_fhigh],
409 ['Lower cutoff freq for the ringdown analysis: %s Hz'%ring_flow],
410 ['Waveform approximant: %s'%(waveform)],
411 ['Final mass/spin fitting formula: %s'%(fit_formula)],
412 ['No deviation from GR above %.1f%% confidence level'%(100.*gr_conf_level)]]
413 np.savetxt('%s/summary_table.txt'%(out_dir), np.array(param_table), delimiter='\t', fmt='%s')
414
415 # save results
416 np.savetxt(out_dir+'/data/Mfchif.dat.gz', (Mf_bins,chif_bins))
417 np.savetxt(out_dir+'/data/P_Mfchif_i.dat.gz', P_Mfchif_i)
418 np.savetxt(out_dir+'/data/P_Mfchif_r.dat.gz', P_Mfchif_r)
419 np.savetxt(out_dir+'/data/P_Mfchif_imr.dat.gz', P_Mfchif_imr)
420 np.savetxt(out_dir+'/data/dMfbyMf_vec.dat.gz', dMfbyMf_vec)
421 np.savetxt(out_dir+'/data/dchifbychif_vec.dat.gz', dchifbychif_vec)
422 np.savetxt(out_dir+'/data/P_dMfbyMf_dchifbychif.dat.gz', P_dMfbyMf_dchifbychif)
423 np.savetxt(out_dir+'/data/P_dMfbyMf.dat.gz', P_dMfbyMf)
424 np.savetxt(out_dir+'/data/P_dchifbychif.dat.gz', P_dchifbychif)
425 np.savetxt(out_dir+'/data/GR_confidence.txt', [gr_conf_level])
426
427 #########################################################################################
428
429 #########################################################################################
430 # plotting
431 #########################################################################################
432 #inspiral
433 P_m1m2_i, m1_bins_i, m2_bins_i = np.histogram2d(m1_i, m2_i, bins=50, density=True)
434 P_chi1chi2_i, chi1_bins_i, chi2_bins_i = np.histogram2d(chi1_i, chi2_i, bins=50, density=True)
435
436 P_m1m2_i = P_m1m2_i.T
437 P_chi1chi2_i = P_chi1chi2_i.T
438
439 conf_m1m2_i = confidence(P_m1m2_i)
440 s1_m1m2_i = conf_m1m2_i.height_from_level(0.68)
441 s2_m1m2_i = conf_m1m2_i.height_from_level(0.95)
442
443 conf_chi1chi2_i = confidence(P_chi1chi2_i)
444 s1_chi1chi2_i = conf_chi1chi2_i.height_from_level(0.68)
445 s2_chi1chi2_i = conf_chi1chi2_i.height_from_level(0.95)
446
447 conf_Mfchif_i = confidence(P_Mfchif_i)
448 s1_Mfchif_i = conf_Mfchif_i.height_from_level(0.68)
449 s2_Mfchif_i = conf_Mfchif_i.height_from_level(0.95)
450
451 plt.figure(figsize=(5,5))
452 plt.pcolormesh(m1_bins_i, m2_bins_i, tgr.gf(P_m1m2_i), cmap='YlOrBr')
453 plt.contour(m1_bins_i[:-1], m2_bins_i[:-1], tgr.gf(P_m1m2_i), levels=(s2_m1m2_i,s1_m1m2_i), linewidths=(1,1.5))
454 if plot_injection_lines == True:
455 plt.axvline(x=m1_inj, ls='--', color='k')
456 plt.axhline(y=m2_inj, ls='--', color='k')
457 plt.xlabel(r'$m_1 [M_{\odot}]$')
458 plt.ylabel(r'$m_2 [M_{\odot}]$')
459 plt.xlim([min(m1_i), max(m1_i)])
460 plt.ylim([min(m2_i), max(m2_i)])
461 plt.grid()
462 plt.savefig('%s/img/inspiral_m1m2_thumb.png'%(out_dir), dpi=72)
463 plt.savefig('%s/img/inspiral_m1m2.png'%(out_dir), dpi=300)
464
465 plt.figure(figsize=(5,5))
466 plt.plot(m1_i, m2_i, 'k.', ms=0.1)
467 if plot_injection_lines == True:
468 plt.axvline(x=m1_inj, ls='--', color='k')
469 plt.axhline(y=m2_inj, ls='--', color='k')
470 plt.xlabel(r'$m_1 [M_{\odot}]$')
471 plt.ylabel(r'$m_2 [M_{\odot}]$')
472 plt.xlim([min(m1_i), max(m1_i)])
473 plt.ylim([min(m2_i), max(m2_i)])
474 plt.grid()
475 plt.savefig('%s/img/inspiral_m1m2_scatter_thumb.png'%(out_dir), dpi=72)
476 plt.savefig('%s/img/inspiral_m1m2_scatter.png'%(out_dir), dpi=300)
477
478 plt.figure(figsize=(5,5))
479 plt.pcolormesh(chi1_bins_i, chi2_bins_i, tgr.gf(P_chi1chi2_i), cmap='YlOrBr')
480 plt.contour(chi1_bins_i[:-1], chi2_bins_i[:-1], tgr.gf(P_chi1chi2_i), levels=(s2_chi1chi2_i,s1_chi1chi2_i), linewidths=(1,1.5))
481 if plot_injection_lines == True:
482 plt.axvline(x=chi1_inj, ls='--', color='k')
483 plt.axhline(y=chi2_inj, ls='--', color='k')
484 plt.xlabel(r'$\chi _1$')
485 plt.ylabel(r'$\chi _2$')
486 plt.xlim([min(chi1_i), max(chi1_i)])
487 plt.ylim([min(chi2_i), max(chi2_i)])
488 plt.grid()
489 plt.savefig('%s/img/inspiral_chi1chi2_thumb.png'%(out_dir), dpi=72)
490 plt.savefig('%s/img/inspiral_chi1chi2.png'%(out_dir), dpi=300)
491
492 plt.figure(figsize=(5,5))
493 plt.plot(chi1_i, chi2_i, 'k.', ms=0.1)
494 if plot_injection_lines == True:
495 plt.axvline(x=chi1_inj, ls='--', color='k')
496 plt.axhline(y=chi2_inj, ls='--', color='k')
497 plt.xlabel(r'$\chi _1$')
498 plt.ylabel(r'$\chi _2$')
499 plt.xlim([min(chi1_i), max(chi1_i)])
500 plt.ylim([min(chi2_i), max(chi2_i)])
501 plt.grid()
502 plt.savefig('%s/img/inspiral_chi1chi2_scatter_thumb.png'%(out_dir), dpi=72)
503 plt.savefig('%s/img/inspiral_chi1chi2_scatter.png'%(out_dir), dpi=300)
504
505 plt.figure(figsize=(5,5))
506 plt.pcolormesh(Mf_bins, chif_bins, tgr.gf(P_Mfchif_i), cmap='YlOrBr')
507 plt.contour(Mf_bins[:-1], chif_bins[:-1], tgr.gf(P_Mfchif_i), levels=(s2_Mfchif_i,s1_Mfchif_i), linewidths=(1,1.5))
508 if plot_injection_lines == True:
509 plt.axvline(x=Mf_inj, ls='--', color='k')
510 plt.axhline(y=chif_inj, ls='--', color='k')
511 plt.xlabel(r'$M_f [M_{\odot}]$')
512 plt.ylabel(r'$\chi_f$')
513 plt.xlim([min(Mf_i), max(Mf_i)])
514 plt.ylim([min(chif_i), max(chif_i)])
515 plt.grid()
516 plt.savefig('%s/img/inspiral_Mfchif_thumb.png'%(out_dir), dpi=72)
517 plt.savefig('%s/img/inspiral_Mfchif.png'%(out_dir), dpi=300)
518
519
520 #ringdown
521 P_m1m2_r, m1_bins_r, m2_bins_r = np.histogram2d(m1_r, m2_r, bins=50, density=True)
522 P_chi1chi2_r, chi1_bins_r, chi2_bins_r = np.histogram2d(chi1_r, chi2_r, bins=50, density=True)
523
524 P_m1m2_r = P_m1m2_r.T
525 P_chi1chi2_r = P_chi1chi2_r.T
526
527 conf_m1m2_r = confidence(P_m1m2_r)
528 s1_m1m2_r = conf_m1m2_r.height_from_level(0.68)
529 s2_m1m2_r = conf_m1m2_r.height_from_level(0.95)
530
531 conf_chi1chi2_r = confidence(P_chi1chi2_r)
532 s1_chi1chi2_r = conf_chi1chi2_r.height_from_level(0.68)
533 s2_chi1chi2_r = conf_chi1chi2_r.height_from_level(0.95)
534
535 conf_Mfchif_r = confidence(P_Mfchif_r)
536 s1_Mfchif_r = conf_Mfchif_r.height_from_level(0.68)
537 s2_Mfchif_r = conf_Mfchif_r.height_from_level(0.95)
538
539 plt.figure(figsize=(5,5))
540 plt.pcolormesh(m1_bins_r, m2_bins_r, tgr.gf(P_m1m2_r), cmap='YlOrBr')
541 plt.contour(m1_bins_r[:-1], m2_bins_r[:-1], tgr.gf(P_m1m2_r), levels=(s2_m1m2_r,s1_m1m2_r), linewidths=(1,1.5))
542 if plot_injection_lines == True:
543 plt.axvline(x=m1_inj, ls='--', color='k')
544 plt.axhline(y=m2_inj, ls='--', color='k')
545 plt.xlabel(r'$m_1 [M_{\odot}]$')
546 plt.ylabel(r'$m_2 [M_{\odot}]$')
547 plt.xlim([min(m1_r), max(m1_r)])
548 plt.ylim([min(m2_r), max(m2_r)])
549 plt.grid()
550 plt.savefig('%s/img/ringdown_m1m2.png'%(out_dir), dpi=300)
551 plt.savefig('%s/img/ringdown_m1m2_thumb.png'%(out_dir), dpi=72)
552
553 plt.figure(figsize=(5,5))
554 plt.plot(m1_r, m2_r, 'k.', ms=0.1)
555 if plot_injection_lines == True:
556 plt.axvline(x=m1_inj, ls='--', color='k')
557 plt.axhline(y=m2_inj, ls='--', color='k')
558 plt.xlabel(r'$m_1 [M_{\odot}]$')
559 plt.ylabel(r'$m_2 [M_{\odot}]$')
560 plt.xlim([min(m1_r), max(m1_r)])
561 plt.ylim([min(m2_r), max(m2_r)])
562 plt.grid()
563 plt.savefig('%s/img/ringdown_m1m2_scatter_thumb.png'%(out_dir), dpi=72)
564 plt.savefig('%s/img/ringdown_m1m2_scatter.png'%(out_dir), dpi=300)
565
566 plt.figure(figsize=(5,5))
567 plt.pcolormesh(chi1_bins_r, chi2_bins_r, tgr.gf(P_chi1chi2_r), cmap='YlOrBr')
568 plt.contour(chi1_bins_r[:-1], chi2_bins_r[:-1], tgr.gf(P_chi1chi2_r), levels=(s2_chi1chi2_r,s1_chi1chi2_r), linewidths=(1,1.5))
569 if plot_injection_lines == True:
570 plt.axvline(x=chi1_inj, ls='--', color='k')
571 plt.axhline(y=chi2_inj, ls='--', color='k')
572 plt.xlabel(r'$\chi _1$')
573 plt.ylabel(r'$\chi _2$')
574 plt.xlim([min(chi1_r), max(chi1_r)])
575 plt.ylim([min(chi2_r), max(chi2_r)])
576 plt.grid()
577 plt.savefig('%s/img/ringdown_chi1chi2_thumb.png'%(out_dir), dpi=72)
578 plt.savefig('%s/img/ringdown_chi1chi2.png'%(out_dir), dpi=300)
579
580 plt.figure(figsize=(5,5))
581 plt.plot(chi1_r, chi2_r, 'k.', ms=0.1)
582 if plot_injection_lines == True:
583 plt.axvline(x=chi1_inj, ls='--', color='k')
584 plt.axhline(y=chi2_inj, ls='--', color='k')
585 plt.xlabel(r'$\chi _1$')
586 plt.ylabel(r'$\chi _2$')
587 plt.xlim([min(chi1_r), max(chi1_r)])
588 plt.ylim([min(chi2_r), max(chi2_r)])
589 plt.grid()
590 plt.savefig('%s/img/ringdown_chi1chi2_scatter_thumb.png'%(out_dir), dpi=72)
591 plt.savefig('%s/img/ringdown_chi1chi2_scatter.png'%(out_dir), dpi=300)
592
593 plt.figure(figsize=(5,5))
594 plt.pcolormesh(Mf_bins, chif_bins, tgr.gf(P_Mfchif_r), cmap='YlOrBr')
595 plt.contour(Mf_bins[:-1], chif_bins[:-1], tgr.gf(P_Mfchif_r), levels=(s2_Mfchif_r,s1_Mfchif_r), linewidths=(1,1.5))
596 if plot_injection_lines == True:
597 plt.axvline(x=Mf_inj, ls='--', color='k')
598 plt.axhline(y=chif_inj, ls='--', color='k')
599 plt.xlabel(r'$M_f [M_{\odot}]$')
600 plt.ylabel(r'$\chi_f$')
601 plt.xlim([min(Mf_r), max(Mf_r)])
602 plt.ylim([min(chif_r), max(chif_r)])
603 plt.grid()
604 plt.savefig('%s/img/ringdown_Mfchif.png'%(out_dir), dpi=300)
605 plt.savefig('%s/img/ringdown_Mfchif_thumb.png'%(out_dir), dpi=72)
606
607 #IMR
608 P_m1m2_imr, m1_bins_imr, m2_bins_imr = np.histogram2d(m1_imr, m2_imr, bins=50, density=True)
609 P_chi1chi2_imr, chi1_bins_imr, chi2_bins_imr = np.histogram2d(chi1_imr, chi2_imr, bins=50, density=True)
610
611 P_m1m2_imr = P_m1m2_imr.T
612 P_chi1chi2_imr = P_chi1chi2_imr.T
613
614 conf_m1m2_imr = confidence(P_m1m2_imr)
615 s1_m1m2_imr = conf_m1m2_imr.height_from_level(0.68)
616 s2_m1m2_imr = conf_m1m2_imr.height_from_level(0.95)
617
618 conf_chi1chi2_imr = confidence(P_chi1chi2_imr)
619 s1_chi1chi2_imr = conf_chi1chi2_imr.height_from_level(0.68)
620 s2_chi1chi2_imr = conf_chi1chi2_imr.height_from_level(0.95)
621
622 conf_Mfchif_imr = confidence(P_Mfchif_imr)
623 s1_Mfchif_imr = conf_Mfchif_imr.height_from_level(0.68)
624 s2_Mfchif_imr = conf_Mfchif_imr.height_from_level(0.95)
625
626 plt.figure(figsize=(5,5))
627 plt.pcolormesh(m1_bins_imr, m2_bins_imr, tgr.gf(P_m1m2_imr), cmap='YlOrBr')
628 plt.contour(m1_bins_imr[:-1], m2_bins_imr[:-1], tgr.gf(P_m1m2_imr), levels=(s2_m1m2_imr,s1_m1m2_imr), linewidths=(1,1.5))
629 if plot_injection_lines == True:
630 plt.axvline(x=m1_inj, ls='--', color='k')
631 plt.axhline(y=m2_inj, ls='--', color='k')
632 plt.xlabel(r'$m_1 [M_{\odot}]$')
633 plt.ylabel(r'$m_2 [M_{\odot}]$')
634 plt.xlim([min(m1_imr), max(m1_imr)])
635 plt.ylim([min(m2_imr), max(m2_imr)])
636 plt.grid()
637 plt.savefig('%s/img/imr_m1m2.png'%(out_dir), dpi=300)
638 plt.savefig('%s/img/imr_m1m2_thumb.png'%(out_dir), dpi=72)
639
640 plt.figure(figsize=(5,5))
641 plt.plot(m1_imr, m2_imr, 'k.', ms=0.1)
642 if plot_injection_lines == True:
643 plt.axvline(x=m1_inj, ls='--', color='k')
644 plt.axhline(y=m2_inj, ls='--', color='k')
645 plt.xlabel(r'$m_1 [M_{\odot}]$')
646 plt.ylabel(r'$m_2 [M_{\odot}]$')
647 plt.xlim([min(m1_imr), max(m1_imr)])
648 plt.ylim([min(m2_imr), max(m2_imr)])
649 plt.grid()
650 plt.savefig('%s/img/imr_m1m2_scatter_thumb.png'%(out_dir), dpi=72)
651 plt.savefig('%s/img/imr_m1m2_scatter.png'%(out_dir), dpi=300)
652
653 plt.figure(figsize=(5,5))
654 plt.pcolormesh(chi1_bins_imr, chi2_bins_imr, tgr.gf(P_chi1chi2_imr), cmap='YlOrBr')
655 plt.contour(chi1_bins_imr[:-1], chi2_bins_imr[:-1], tgr.gf(P_chi1chi2_imr), levels=(s2_chi1chi2_imr,s1_chi1chi2_imr), linewidths=(1,1.5))
656 if plot_injection_lines == True:
657 plt.axvline(x=chi1_inj, ls='--', color='k')
658 plt.axhline(y=chi2_inj, ls='--', color='k')
659 plt.xlabel(r'$\chi _1$')
660 plt.ylabel(r'$\chi _2$')
661 plt.xlim([min(chi1_imr), max(chi1_imr)])
662 plt.ylim([min(chi2_imr), max(chi2_imr)])
663 plt.grid()
664 plt.savefig('%s/img/imr_chi1chi2_thumb.png'%(out_dir), dpi=72)
665 plt.savefig('%s/img/imr_chi1chi2.png'%(out_dir), dpi=300)
666
667 plt.figure(figsize=(5,5))
668 plt.plot(chi1_imr, chi2_imr, 'k.', ms=0.1)
669 if plot_injection_lines == True:
670 plt.axvline(x=chi1_inj, ls='--', color='k')
671 plt.axhline(y=chi2_inj, ls='--', color='k')
672 plt.xlabel(r'$\chi _1$')
673 plt.ylabel(r'$\chi _2$')
674 plt.xlim([min(chi1_imr), max(chi1_imr)])
675 plt.ylim([min(chi2_imr), max(chi2_imr)])
676 plt.grid()
677 plt.savefig('%s/img/imr_chi1chi2_scatter_thumb.png'%(out_dir), dpi=72)
678 plt.savefig('%s/img/imr_chi1chi2_scatter.png'%(out_dir), dpi=300)
679
680 plt.figure(figsize=(5,5))
681 plt.pcolormesh(Mf_bins, chif_bins, tgr.gf(P_Mfchif_imr), cmap='YlOrBr')
682 plt.contour(Mf_bins[:-1], chif_bins[:-1], tgr.gf(P_Mfchif_imr), levels=(s2_Mfchif_imr,s1_Mfchif_imr), linewidths=(1,1.5))
683 if plot_injection_lines == True:
684 plt.axvline(x=Mf_inj, ls='--', color='k')
685 plt.axhline(y=chif_inj, ls='--', color='k')
686 plt.xlabel(r'$M_f [M_{\odot}]$')
687 plt.ylabel(r'$\chi_f$')
688 plt.xlim([min(Mf_imr), max(Mf_imr)])
689 plt.ylim([min(chif_imr), max(chif_imr)])
690 plt.grid()
691 plt.savefig('%s/img/imr_Mfchif.png'%(out_dir), dpi=300)
692 plt.savefig('%s/img/imr_Mfchif_thumb.png'%(out_dir), dpi=72)
693
694 # IR overlap
695 plt.figure(figsize=(5,5))
696 CSi = plt.contour(Mf_bins[:-1], chif_bins[:-1], tgr.gf(P_Mfchif_i), levels=(s2_Mfchif_i,s1_Mfchif_i), linewidths=(1,1.5), colors='orange')
697 CSr = plt.contour(Mf_bins[:-1], chif_bins[:-1], tgr.gf(P_Mfchif_r), levels=(s2_Mfchif_r,s1_Mfchif_r), linewidths=(1,1.5), colors='red')
698 CSimr = plt.contour(Mf_bins[:-1], chif_bins[:-1], tgr.gf(P_Mfchif_imr), levels=(s2_Mfchif_imr,s1_Mfchif_imr), linewidths=(1,1.5), colors='k')
699 if plot_injection_lines == True:
700 plt.axvline(x=Mf_inj, ls='--', color='k')
701 plt.axhline(y=chif_inj, ls='--', color='k')
702 plt.xlim([min(np.append(Mf_i, Mf_r)), max(np.append(Mf_i, Mf_r))])
703 plt.ylim([min(np.append(chif_i, chif_r)), max(np.append(chif_i, chif_r))])
704 plt.xlabel(r'$M_f~[M_\odot]$')
705 plt.ylabel(r'$\chi_f$')
706 plt.grid()
707
708 CSi.levels = np.asarray(CSi.levels)
709 CSr.levels = np.asarray(CSr.levels)
710 CSimr.levels = np.asarray(CSimr.levels)
711
712 strs_i = [ 'inspiral', 'inspiral' ]
713 strs_r = [ 'ringdown', 'ringdown' ]
714 strs_imr = [ 'IMR', 'IMR' ]
715 fmt_i = {}
716 fmt_r = {}
717 fmt_imr = {}
718 for l,s in zip(CSi.levels, strs_i):
719 fmt_i[l] = s
720 for l,s in zip(CSr.levels, strs_r):
721 fmt_r[l] = s
722 for l,s in zip(CSimr.levels, strs_imr):
723 fmt_imr[l] = s
724
725 ## Label every other level using strings
726 plt.clabel(CSi,CSi.levels[::2],inline=True,fmt=fmt_i,fontsize=14, use_clabeltext=True)
727 plt.clabel(CSr,CSr.levels[::2],inline=True,fmt=fmt_r,fontsize=14, use_clabeltext=True)
728 plt.clabel(CSimr,CSimr.levels[::2],inline=True,fmt=fmt_imr,fontsize=10)
729
730 plt.savefig('%s/img/IMR_overlap.png'%(out_dir), dpi=300)
731 plt.savefig('%s/img/IMR_overlap_thumb.png'%(out_dir), dpi=72)
732
733 #(dMf/Mf, dchif/chif)
734 conf_v1v2 = confidence(P_dMfbyMf_dchifbychif)
735 s1_v1v2 = conf_v1v2.height_from_level(0.68)
736 s2_v1v2 = conf_v1v2.height_from_level(0.95)
737
738 conf_v1 = confidence(P_dMfbyMf)
739 s1_v1 = conf_v1.height_from_level(0.68)
740 s2_v1 = conf_v1.height_from_level(0.95)
741
742 conf_v2 = confidence(P_dchifbychif)
743 s1_v2 = conf_v2.height_from_level(0.68)
744 s2_v2 = conf_v2.height_from_level(0.95)
745
746 # Calculation of condifence edges
747 left1_v1 = min(dMfbyMf_vec[np.where(P_dMfbyMf>=s1_v1)[0]])
748 right1_v1 = max(dMfbyMf_vec[np.where(P_dMfbyMf>=s1_v1)[0]])
749
750 left2_v1 = min(dMfbyMf_vec[np.where(P_dMfbyMf>=s2_v1)[0]])
751 right2_v1 = max(dMfbyMf_vec[np.where(P_dMfbyMf>=s2_v1)[0]])
752
753 left1_v2 = min(dchifbychif_vec[np.where(P_dchifbychif>s1_v2)[0]])
754 right1_v2 = max(dchifbychif_vec[np.where(P_dchifbychif>s1_v2)[0]])
755
756 left2_v2 = min(dchifbychif_vec[np.where(P_dchifbychif>s2_v2)[0]])
757 right2_v2 = max(dchifbychif_vec[np.where(P_dchifbychif>s2_v2)[0]])
758
759 plt.figure(figsize=(5,5))
760 plt.subplot2grid((3,3), (0,0), colspan=2)
761 plt.plot(dMfbyMf_vec, P_dMfbyMf, color='k', lw=1)
762 plt.axvline(x=left1_v1, color='c', lw=0.5, ls='-.')
763 plt.axvline(x=right1_v1, color='c', lw=0.5, ls='-.')
764 plt.axvline(x=left2_v1, color='b', lw=0.5, ls='-.')
765 plt.axvline(x=right2_v1, color='b', lw=0.5, ls='-.')
766 #plt.xlabel(r'$\Delta M_f/M_f$')
767 plt.ylabel(r'$P(\Delta M_f/M_f)$')
768 #plt.grid()
769
770 plt.subplot2grid((3,3), (1,0), colspan=2, rowspan=2)
771 plt.pcolormesh(dMfbyMf_vec,dchifbychif_vec,P_dMfbyMf_dchifbychif, cmap='YlOrBr')
772 plt.contour(dMfbyMf_vec,dchifbychif_vec,tgr.gf(P_dMfbyMf_dchifbychif), levels=(s2_v1v2,s1_v1v2), linewidths=(1,1.5))
773 plt.plot(0, 0, 'k+', ms=12, mew=2)
774 plt.xlabel(r'$\Delta M_f/M_f$')
775 plt.ylabel(r'$\Delta \chi_f/\chi_f$')
776 plt.xlim([-dMfbyMf_lim,dMfbyMf_lim])
777 plt.ylim([-dchifbychif_lim,dchifbychif_lim])
778 plt.grid()
779
780 plt.subplot2grid((3,3), (1,2), rowspan=2)
781 plt.plot(P_dchifbychif, dchifbychif_vec,'k', lw=1)
782 plt.axhline(y=left1_v2, color='c', lw=0.5, ls='-.')
783 plt.axhline(y=right1_v2, color='c', lw=0.5, ls='-.')
784 plt.axhline(y=left2_v2, color='b', lw=0.5, ls='-.')
785 plt.axhline(y=right2_v2, color='b', lw=0.5, ls='-.')
786 #plt.ylabel(r'$\Delta \chi_f/\chi_f$')
787 plt.xlabel(r'$P(\Delta \chi_f/\chi_f)$')
788 #plt.grid()
789
790 plt.savefig('%s/img/dMfbyMfdchifbychif.png' %(out_dir), dpi=300)
791 plt.savefig('%s/img/dMfbyMfdchifbychif_thumb.png' %(out_dir), dpi=72)
792
793 print('... made summary plots')
794
795 print('... completed in %f seconds' %(time.time()-start_time))
796 #########################################################################################
max
#define max(a, b)
Definition: LALInferenceRemoveLines.c:80
imrtgr_imr_consistency_test.confidence
Definition: imrtgr_imr_consistency_test.py:56
imrtgr_imr_consistency_test.confidence.__init__
def __init__(self, counts)
Definition: imrtgr_imr_consistency_test.py:57
imrtgr_imr_consistency_test.confidence._height_from_idx
_height_from_idx
Definition: imrtgr_imr_consistency_test.py:67
imrtgr_imr_consistency_test.confidence._idx_from_height
_idx_from_height
Definition: imrtgr_imr_consistency_test.py:69
imrtgr_imr_consistency_test.confidence._set_interp
def _set_interp(self)
Definition: imrtgr_imr_consistency_test.py:64
imrtgr_imr_consistency_test.confidence.counts_sorted
counts_sorted
Definition: imrtgr_imr_consistency_test.py:59
imrtgr_imr_consistency_test.confidence.level_from_height
def level_from_height(self, height)
Definition: imrtgr_imr_consistency_test.py:74
imrtgr_imr_consistency_test.confidence.height_from_level
def height_from_level(self, level)
Definition: imrtgr_imr_consistency_test.py:76
imrtgr_imr_consistency_test.confidence._length
_length
Definition: imrtgr_imr_consistency_test.py:65
imrtgr_imr_consistency_test.confidence._idx_from_level
_idx_from_level
Definition: imrtgr_imr_consistency_test.py:73
imrtgr_imr_consistency_test.confidence._level_from_idx
_level_from_idx
Definition: imrtgr_imr_consistency_test.py:71
imrtgr_imr_consistency_test.confidence.norm_cumsum_counts_sorted
norm_cumsum_counts_sorted
Definition: imrtgr_imr_consistency_test.py:61
imrtgr_imr_consistency_test.set_tick_sizes
def set_tick_sizes(ax, major, minor)
Definition: imrtgr_imr_consistency_test.py:46
imrtgr_imr_consistency_test.P_Mfchif_pr_interp_obj
P_Mfchif_pr_interp_obj
Definition: imrtgr_imr_consistency_test.py:350
lalinference_cpnest.int
int
Definition: lalinference_cpnest.py:38