Figures for GW190521 publications

This notebook generates (almost) all the Figures in the publications
GW190521: A Binary Black Hole Merger with a Total Mass of 150 Msun avaliable through PRL, arXiv, and LIGO-P2000020.
Properties and astrophysical implications of the 150 Msun binary black hole merger GW190521 avaliable through ApjL, arXiv, and LIGO-P2000021.

The data used in this notebook must be downloaded from the public LIGO DCC page LIGO-P2000158.

Figures 2,3,4,5 in LIGO-P2000020 display the parameter estimation results using only the preferred waveform model (NRSur PHM) while the corresponding figures 1,2,6,9 in LIGO-P2000021 also display results from two other waveform models (Phenom PHM and SEOBNR PHM). This notebook generates the figures in LIGO-P2000021.

Figure 1 from LIGO-P2000020 is produced with a separate notebook, GW190521_discovery_Fig1.ipynb, also available from LIGO-P2000158.

Figures 11 and 12 from LIGO-P2000021 are produced with a separate notebook, GW190521_Hierarchical_plots.ipynb, also available from LIGO-P2000158.

Data release for GW190521

In public DCC page LIGO-P2000158, you will find several data files:

  • GW190521_posterior_samples.h5 - contains parameter estimation (PE) samples from many different runs:
    • NRSur7dq4 is our "preferred" waveform model (referred to as NRSur PHM in the papers). These samples contain the info for Table 1 and Figures 2, 3, 4, 5 of LIGO-P2000020.
    • IMRPhenomPv3HM and SEOBNRv4PHM are alternative waveform models for comparison. Those samples, along with those from NRSur7dq4, contain the info for Table 1 and Figures 1, 2, 3, 4, 5, 6, 8, 9 of LIGO-P2000021.
  • GW190521_studies_posterior_samples.h5 - contains parameter estimation (PE) samples from runs used for specialized studies; these should not be used in preference to the previous runs for most applications.
    • The three sets of PE samples for NRSur7dq4_Lmax2, Lmax3, and Lmax4 are used for Figure 7 of LIGO-P2000021.
    • The set of PE samples directly using NR is used for Figure 8 of LIGO-P2000021, and nowhere else in either paper. It has luminosity_distance set to 1 Mpc and redshift close to 0 for all samples, and hence all masses labeled _source are not actually corrected for redshift. The radiated energy also does not have the redshift correction. The remnant BH properties and peak luminosity are not computed using spin evolution and are thus slightly less accurate than for other runs.
    • The set of PE samples for NRSur7dq4_nospin and NRSur7dq4_spinaligned are used for computing Bayes Factors for spin and for higher multipoles, and for the hierarchical analysis described in Section 5.2.1 and Figures 11 and 12 of LIGO-P2000021.
  • GW190521_Ringdown_samples.tgz is a tarball containing nine more h5 files with posterior samples for the ringdown analysis described in Section 3.2 and Figure 9 of LIGO-P2000021.
  • GW190521_Implications_figure_data.tgz is a tarball containing additional data needed to make Figures 5, 10, 11, 12, and 13 in LIGO-P2000021, including skymaps fits files.
  • GW190521_discovery_Fig1.tgz is a tarball containing additional data needed to make Figure 1 in LIGO-P2000020.
  • GW190521_discovery_figs_pdf.tgz - tarball containing all the figures from the GW190521 discovery paper LIGO-P2000020, in pdf.
  • GW190521_Implications_Figures_pdf.tgz - tarball containing all the figures in the GW190521 implications paper LIGO-P2000021, in pdf.
  • GW190521_md5sums.txt - containing md5sums for each of the h5 files.
In [1]:
# download the data used in this notebook by un-commenting out these shell commands (keep the !),
# or excecuting the commands in a shell, in the same directory as this notebook.

#!curl https://dcc.ligo.org/public/0168/P2000158/004/GW190521_posterior_samples.h5 -O
#!curl https://dcc.ligo.org/public/0168/P2000158/004/GW190521_studies_posterior_samples.h5 -O
#!curl https://dcc.ligo.org/public/0168/P2000158/004/GW190521_Ringdown_samples.tar -O
#!tar xvzf Ringdown_samples_GW190521.tgz
#!curl https://dcc.ligo.org/public/0168/P2000158/004/GW190521_Implications_figure_data.tar -O
#!tar xvzf GW190521_Implications_figure_data.tgz
#!curl https://dcc.ligo.org/public/0168/P2000158/004/GW190521_md5sums.txt -O
In [2]:
# generic python imports
from __future__ import division
import numpy as np
from IPython.display import HTML, display

import scipy
from scipy.special import erf
from scipy.stats import gaussian_kde
from scipy import integrate

import os,sys,subprocess
import h5py 
import hashlib  
import corner

import matplotlib
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import matplotlib.patches as patches
from matplotlib.lines import Line2D as ld

# figure fonts
matplotlib.rc('text', usetex=True)  
matplotlib.rcParams['text.latex.preamble'] = [r"\usepackage{bm}"]
matplotlib.rcParams['mathtext.fontset'] = 'stix'
matplotlib.rcParams['font.family'] = 'STIXGeneral'

# LIGO-specific software; install with conda or pip
import lal                                          # https://pypi.org/project/lalsuite/
from pesummary.gw.file.read import read as GWread   # https://pypi.org/project/pesummary/
from ligo.skymap.io import fits                     # https://pypi.org/project/ligo.skymap/
In [3]:
# put figures in a sub-directory:
figpath = 'GW190521_Implications_Figures_pdf/' 
if not os.path.exists(figpath):
    os.makedirs(figpath)
In [4]:
#===============
# PE SAMPLES 
#===============

# read in PE sample files, check md5sum against https://dcc.ligo.org/public/0166/P2000158/002/GW190521_md5sums.txt

#------------------------------------
# a. Inspiral-Merger-Ringdown waveform model Posterior PE samples
#------------------------------------

posterior_samples_file='GW190521_posterior_samples.h5'
print("md5sum",posterior_samples_file,"=",hashlib.md5(open(posterior_samples_file,'rb').read()).hexdigest())
samples_file = h5py.File(posterior_samples_file,'r')

nrsur = samples_file['NRSur7dq4']['posterior_samples']
phenom = samples_file['IMRPhenomPv3HM']['posterior_samples']
seob = samples_file['SEOBNRv4PHM']['posterior_samples']

prior = samples_file['NRSur7dq4']['priors']['samples']
 
md5sum GW190521_posterior_samples.h5 = 8af9bce0b55b5ebed7853dbfaa69a2d5
In [5]:
#------------------------------------
# b. Higher-Mode Analysis - NRSur7dq4 runs with RIFT Pipeline
#------------------------------------

###HM_samples_fn = '/home/charlie.hoy/public_html/public_release/S190521g/v2/NRSur7dq4_restrictL/samples/posterior_samples.h5'
###HM_samples_fn = 'NRSur7dq4_restrictL_posterior_samples.h5'
HM_samples_fn = 'GW190521_studies_posterior_samples.h5'
print("md5sum",HM_samples_fn,"=",hashlib.md5(open(os.path.join('./',HM_samples_fn),'rb').read()).hexdigest())
HM_samples_file = h5py.File(HM_samples_fn,'r')

#l=2
nrs_l2 = HM_samples_file['NRSur7dq4_Lmax2']['posterior_samples']

#l=2:3
nrs_l3 = HM_samples_file['NRSur7dq4_Lmax3']['posterior_samples']

#l=2:4
nrs_l4 = HM_samples_file['NRSur7dq4_Lmax4']['posterior_samples']
md5sum GW190521_studies_posterior_samples.h5 = 880c79fbfc4058e9dd8a2bbf5841c691
In [6]:
#------------------------------------
# c. Numerical Relativity Runs from RIFT Pipeline
#------------------------------------

###NR_samples_fn = '/home/charlie.hoy/public_html/public_release/S190521g/v1/NR/samples/posterior_samples.h5'
###NR_samples_fn = 'GW190521_NR_posterior_samples.h5'
NR_samples_fn = 'GW190521_studies_posterior_samples.h5'
print("md5sum",NR_samples_fn,"=",hashlib.md5(open(os.path.join('./',NR_samples_fn),'rb').read()).hexdigest())
NR_samples_file = h5py.File(NR_samples_fn,'r')
nr_sim = NR_samples_file['NR']['posterior_samples']
md5sum GW190521_studies_posterior_samples.h5 = 880c79fbfc4058e9dd8a2bbf5841c691
In [7]:
#------------------------------------
# d. Ringdown Runs from pyRing Pipeline
#------------------------------------
# Copied tarball from https://dcc.ligo.org/public/0166/P2000158/002/GW190521_Ringdown_samples.tar
# check the md5sums against https://dcc.ligo.org/public/0166/P2000158/002/GW190521_md5sums.txt

ringdown_samples_dir = 'GW190521_Ringdown_samples/'

fnh5 = 'GW190521A_PROD0_DS_1mode_5M_pesummary_metafile.h5'
print("md5sum",fnh5,"=",hashlib.md5(open(os.path.join(ringdown_samples_dir,fnh5),'rb').read()).hexdigest())
damped_sinusoid_5M   = GWread(os.path.join(ringdown_samples_dir,fnh5))
Damped_sinusoid_5M   = damped_sinusoid_5M.samples_dict['GW190521A_PROD0_DS_1mode_5M']

fnh5 = 'GW190521A_PROD0_DS_1mode_10M_pesummary_metafile.h5'
print("md5sum",fnh5,"=",hashlib.md5(open(os.path.join(ringdown_samples_dir,fnh5),'rb').read()).hexdigest())
damped_sinusoid_10M   = GWread(os.path.join(ringdown_samples_dir,fnh5))
Damped_sinusoid_10M   = damped_sinusoid_10M.samples_dict['GW190521A_PROD0_DS_1mode_10M']

fnh5 = 'GW190521A_PROD0_DS_1mode_15M_pesummary_metafile.h5'
print("md5sum",fnh5,"=",hashlib.md5(open(os.path.join(ringdown_samples_dir,fnh5),'rb').read()).hexdigest())
damped_sinusoid_15M   = GWread(os.path.join(ringdown_samples_dir,fnh5))
Damped_sinusoid_15M   = damped_sinusoid_15M.samples_dict['GW190521A_PROD0_DS_1mode_15M']

fnh5 = 'GW190521A_PROD0_Kerr_220_10M_pesummary_metafile.h5'
print("md5sum",fnh5,"=",hashlib.md5(open(os.path.join(ringdown_samples_dir,fnh5),'rb').read()).hexdigest())
kerr_220   = GWread(os.path.join(ringdown_samples_dir,fnh5))
Kerr_220   = kerr_220.samples_dict['GW190521A_PROD0_Kerr_220_10M']

fnh5 = 'GW190521A_PROD0_Kerr_222_0M_pesummary_metafile.h5'
print("md5sum",fnh5,"=",hashlib.md5(open(os.path.join(ringdown_samples_dir,fnh5),'rb').read()).hexdigest())
kerr_222   = GWread(os.path.join(ringdown_samples_dir,fnh5))
Kerr_222   = kerr_222.samples_dict['GW190521A_PROD0_Kerr_222_0M']

fnh5 = 'GW190521A_PROD0_MMRDNP_10M_pesummary_metafile.h5'
print("md5sum",fnh5,"=",hashlib.md5(open(os.path.join(ringdown_samples_dir,fnh5),'rb').read()).hexdigest())
kerr_HM   = GWread(os.path.join(ringdown_samples_dir,fnh5))
Kerr_HM   = kerr_HM.samples_dict['GW190521A_PROD0_MMRDNP_10M']

fnh5 = 'GW190521A_NRSur7dq4_ftau_220_posterior_samples.h5'
print("md5sum",fnh5,"=",hashlib.md5(open(os.path.join(ringdown_samples_dir,fnh5),'rb').read()).hexdigest())
nrsur_ftau   = GWread(os.path.join(ringdown_samples_dir,fnh5))
nrsur_ftau   = nrsur_ftau.samples_dict['NRSur7dq4']

fnh5 = 'GW190521A_SEOBNRv4PHM_ftau_220_posterior_samples.h5'
print("md5sum",fnh5,"=",hashlib.md5(open(os.path.join(ringdown_samples_dir,fnh5),'rb').read()).hexdigest())
seob_ftau   = GWread(os.path.join(ringdown_samples_dir,fnh5))
seob_ftau   = seob_ftau.samples_dict['SEOBNRv4PHM']

fnh5 = 'GW190521A_IMRPhenomPv3HM_ftau_220_posterior_samples.h5'
print("md5sum",fnh5,"=",hashlib.md5(open(os.path.join(ringdown_samples_dir,fnh5),'rb').read()).hexdigest())
phenom_ftau   = GWread(os.path.join(ringdown_samples_dir,fnh5))
phenom_ftau   = phenom_ftau.samples_dict['IMRPhenomPv3HM']
md5sum GW190521A_PROD0_DS_1mode_5M_pesummary_metafile.h5 = c971e21f2bc317cba5a6a02e3829f8d1
md5sum GW190521A_PROD0_DS_1mode_10M_pesummary_metafile.h5 = 6ef505c83a1f6304c822defe63485d65
md5sum GW190521A_PROD0_DS_1mode_15M_pesummary_metafile.h5 = 70f55c2d74990f0602cc8aa6a49d961d
md5sum GW190521A_PROD0_Kerr_220_10M_pesummary_metafile.h5 = 28d380a482c824d244c331e66ea8efbf
md5sum GW190521A_PROD0_Kerr_222_0M_pesummary_metafile.h5 = 64dfaabbf133e887e9ce73be7c54f50d
md5sum GW190521A_PROD0_MMRDNP_10M_pesummary_metafile.h5 = b309aa330d1c72e7edcf4145ba9347f8
md5sum GW190521A_NRSur7dq4_ftau_220_posterior_samples.h5 = 3cd72b881b5ba462a3f9d525af8d2897
md5sum GW190521A_SEOBNRv4PHM_ftau_220_posterior_samples.h5 = d46d82dc245135c3b0cfb6698b2386fc
md5sum GW190521A_IMRPhenomPv3HM_ftau_220_posterior_samples.h5 = 20d7006063a869456564ab161e6ba609
In [8]:
#========================
# Bounded KDE function, used to smoothe 1-D posterior distributions.
#========================

class Bounded_kde(gaussian_kde):
    r"""Represents a one-dimensional Gaussian kernel density estimator
    for a probability distribution function that exists on a bounded
    domain."""

    def __init__(self, pts, low=None, high=None, *args, **kwargs):
        """Initialize with the given bounds.  Either ``low`` or
        ``high`` may be ``None`` if the bounds are one-sided.  Extra
        parameters are passed to :class:`gaussian_kde`.

        :param low: The lower domain boundary.

        :param high: The upper domain boundary."""
        pts = np.atleast_1d(pts)

        assert pts.ndim == 1, 'Bounded_kde can only be one-dimensional'
        
        super(Bounded_kde, self).__init__(pts, *args, **kwargs)

        self._low = low
        self._high = high

    @property
    def low(self):
        """The lower bound of the domain."""
        return self._low

    @property
    def high(self):
        """The upper bound of the domain."""
        return self._high

    def evaluate(self, xs):
        """Return an estimate of the density evaluated at the given
        points."""
        xs = np.atleast_1d(xs)
        assert xs.ndim == 1, 'points must be one-dimensional'

        pdf = super(Bounded_kde, self).evaluate(xs)

        if self.low is not None:
            pdf += super(Bounded_kde, self).evaluate(2.0*self.low - xs)

        if self.high is not None:
            pdf += super(Bounded_kde, self).evaluate(2.0*self.high - xs)

        return pdf

    __call__ = evaluate
In [9]:
#========================
# PLOTTING FUNCTIONS 
#========================

#-------------------------------------------
# a. Plot with 2D contour + x,y histogram
#-------------------------------------------
def triangle_plot_2d_axes(
    xbounds, ybounds, figsize=(8, 8),
    width_ratios=[4, 1], height_ratios=[1, 4], wspace=0.0, hspace=0.0,
    grid=False,high1d=1):
    """Initialize the axes for a 2d triangle plot.
    """
    high1d = high1d

    fig = plt.figure(figsize=figsize)
    gs = gridspec.GridSpec(
        2, 2,
        width_ratios=width_ratios, height_ratios=height_ratios,
        wspace=wspace, hspace=hspace)

    ax1 = plt.subplot(gs[0])
    ax3 = plt.subplot(gs[2])
    ax4 = plt.subplot(gs[3])

    ax1.minorticks_on()
    ax3.minorticks_on()
    ax4.minorticks_on()

    if grid:
        ax1.grid(which='major', ls='-')
        ax1.grid(which='minor', ls=':')
        ax3.grid(which='major', ls='-')
        ax3.grid(which='minor', ls=':')
        ax4.grid(which='major', ls='-')
        ax4.grid(which='minor', ls=':')

    # Get rid of tick labels
    ax1.xaxis.set_ticklabels([])
    ax4.yaxis.set_ticklabels([])

    # Use consistent x-axis and y-axis bounds in all 3 plots
    ax1.set_ylim(0, high1d)
    ax1.set_xlim(xbounds[0], xbounds[1])
    ax3.set_xlim(xbounds[0], xbounds[1])
    ax3.set_ylim(ybounds[0], ybounds[1])
    ax4.set_xlim(0, high1d)
    ax4.set_ylim(ybounds[0], ybounds[1])

    return fig, ax1, ax3, ax4


def create_fig_and_axes(xbounds, ybounds, figsize=(9.7, 9.7),high1d=1):
    fig, ax1, ax3, ax4 = triangle_plot_2d_axes(
        xbounds, ybounds, figsize=figsize, width_ratios=[4, 1],
        height_ratios=[1, 4], wspace=0.0, hspace=0.0,high1d=high1d)
    return fig, ax1, ax3, ax4





#------------------------------------------------------------
# b. Plotting the main result: NRSur PHM (filled Contour)
#------------------------------------------------------------

def add_samples_to_fig(name, parameter_1,parameter_2,zorder=10,norm_factor_x=1,norm_factor_y=1,bounded_1=False,low_1=None,high_1=None,bounded_2=False,low_2=None,high_2=None):
    
    #name: name of pe result to plot
    #parameter 1 and parameter 2: parameters to plot in the 2D plot
    #zorder: if plottting multiple results, a larger zorder places the contour on top
    #norm_factor_x,y: parameter to scale the corner 1D histograms
    
    x = samples_dict[name][parameter_1]
    y = samples_dict[name][parameter_2]
    xlow, xhigh = xlims
    xsmooth = np.linspace(xlow, xhigh, 1000)
    ylow, yhigh = ylims
    ysmooth = np.linspace(ylow, yhigh, 1000)

    norm_factor_x=0.95*norm_factor_x
    norm_factor_y=0.95*norm_factor_y
    
    c = color_dict[name]
    label = label_dict[name]
    alpha = alpha_dict[name]
    lw = lw_dict[name]
    bb = bins_dict[name]    


    if bounded_1==True:                     
        if parameter_1 == 'chi_p':
            guess = lambda x: x * (1.-x)
            kde = gaussian_kde(x, weights=1./guess(x))
            guess_norm = integrate.quad(guess, low_1, high_1)[0]

            ax1.plot(xsmooth, norm_factor_x* kde(xsmooth) * guess(xsmooth) / guess_norm, color=c, lw=lw, label=label,
             zorder=zorder)
        
            ax1.fill_between(xsmooth, 0, norm_factor_x* kde(xsmooth) * guess(xsmooth) / guess_norm, color=c, alpha=alpha,
                     zorder=zorder)
            
        elif parameter_1 == 'theta_jn':
            guess = lambda x: x * (np.pi-x)
            kde = gaussian_kde(x, weights=1./guess(x))
            guess_norm = integrate.quad(guess, low_1, high_1)[0]

            ax1.plot(xsmooth, norm_factor_x* kde(xsmooth) * guess(xsmooth), color=c, lw=lw, label=label,
             zorder=zorder)
        
            ax1.fill_between(xsmooth, 0, norm_factor_x* kde(xsmooth) * guess(xsmooth), color=c, alpha=alpha,
                     zorder=zorder)
         
        else:       
            print("Using bounded KDE")   
            kde = Bounded_kde(x,low=low_1,high=high_1)
            ax1.plot(xsmooth, norm_factor_x*kde(xsmooth), color=c, lw=lw, label=label,
                     zorder=zorder)
            ax1.fill_between(xsmooth, 0, norm_factor_x*kde(xsmooth), color=c, alpha=alpha,
                             zorder=zorder)            
        
    else:         
        kde = gaussian_kde(x)

        ax1.plot(xsmooth, norm_factor_x*kde(xsmooth), color=c, lw=lw, label=label,
                 zorder=zorder)
        ax1.fill_between(xsmooth, 0, norm_factor_x*kde(xsmooth), color=c, alpha=alpha,
                         zorder=zorder)
    ax1.axvline(np.quantile(x, 0.05), color=c, ls='dashed')
    ax1.axvline(np.quantile(x, 0.95), color=c, ls='dashed')

    if bounded_2==True:        
        print("Using bounded KDE")    
        kde = Bounded_kde(y,low=low_2,high=high_2)    
    else:         
        kde = gaussian_kde(y)
        
    ax4.plot(norm_factor_y*kde(ysmooth), ysmooth, color=c, lw=lw, label=label,
             zorder=zorder)
    ax4.fill_betweenx(ysmooth, 0, norm_factor_y*kde(ysmooth), color=c, alpha=alpha,
                      zorder=zorder)
    ax4.axhline(np.quantile(y, 0.05), color=c, ls='dashed')
    ax4.axhline(np.quantile(y, 0.95), color=c, ls='dashed')
    

###### For the case of m1, m2 reflect samples at m1=m2 to correctly treat the boundary #######
    
    if parameter_1 == 'mass_1_source' and parameter_2 == 'mass_2_source':       
        print("Symmetrising samples around m1=m2")    
        # symmetrize the samples over 1<->2
        x_temp = np.concatenate((x, y))
        y_temp = np.concatenate((y, x))        
        x=x_temp
        y=y_temp    
    else:        
        print("No reflexion done")
    
    my_range = [[xlow, xhigh], [ylow, yhigh]]
    corner.hist2d(x, y, ax=ax3, range=my_range, color=c,
                  plot_datapoints=False, plot_density=True,smooth=True,
                  levels=[0.9],fill_contours=False,bins=bb, lw=4)

    
#----------------------------------------------------------------------------
# c. Plotting other models: Phenom PHM, SEOBNR PHM (transparent contours)
#-----------------------------------------------------------------------------

def add_samples_to_fig_nofilled(name, parameter_1,parameter_2,zorder=10,norm_factor_x=1,norm_factor_y=1,bounded_1=False,low_1=None,high_1=None,bounded_2=False,low_2=None,high_2=None):
    
    x = samples_dict[name][parameter_1]
    y = samples_dict[name][parameter_2]
    xlow, xhigh = xlims
    xsmooth = np.linspace(xlow, xhigh, 1000)
    ylow, yhigh = ylims
    ysmooth = np.linspace(ylow, yhigh, 1000)

    norm_factor_x=0.95*norm_factor_x
    norm_factor_y=0.95*norm_factor_y
    
    c = color_dict[name]
    label = label_dict[name]
    alpha = alpha_dict[name]
    lw = lw_dict[name]
    bb = bins_dict[name]

    if bounded_1==True:                     
        if parameter_1 == 'chi_p':
            guess = lambda x: x * (1.-x)
            kde = gaussian_kde(x, weights=1./guess(x))
            guess_norm = integrate.quad(guess, low_1, high_1)[0]

            ax1.plot(xsmooth, norm_factor_x* kde(xsmooth) * guess(xsmooth) / guess_norm, color=c, lw=lw, label=label,
             zorder=zorder)
            
        elif parameter_1 == 'theta_jn':
            guess = lambda x: x * (np.pi-x)
            kde = gaussian_kde(x, weights=1./guess(x))
            guess_norm = integrate.quad(guess, low_1, high_1)[0]

            ax1.plot(xsmooth, norm_factor_x* kde(xsmooth) * guess(xsmooth), color=c, lw=lw, label=label,
             zorder=zorder)
                  
        else:       
            print("Using bounded KDE")   
            kde = Bounded_kde(x,low=low_1,high=high_1)          
            ax1.plot(xsmooth, norm_factor_x*kde(xsmooth), color=c, lw=lw, label=label,
                     zorder=zorder)            
    else:         
        kde = gaussian_kde(x)    
    
        ax1.plot(xsmooth, norm_factor_x*kde(xsmooth), color=c, lw=lw, label=label,
                 zorder=zorder)
    ax1.axvline(np.quantile(x, 0.05), color=c, ls='dashed')
    ax1.axvline(np.quantile(x, 0.95), color=c, ls='dashed')

    
    if bounded_2==True:        
        print("Using bounded KDE")    
        kde = Bounded_kde(y,low=low_2,high=high_2)    
    else:         
        kde = gaussian_kde(y)    
    
    ax4.plot(norm_factor_y*kde(ysmooth), ysmooth, color=c, lw=lw, label=label,
             zorder=zorder)
    ax4.axhline(np.quantile(y, 0.05), color=c, ls='dashed')
    ax4.axhline(np.quantile(y, 0.95), color=c, ls='dashed')

###### For the case of m1, m2 reflect samples at m1=m2 to correctly treat the boundary #######
    
    if parameter_1 == 'mass_1_source' and parameter_2 == 'mass_2_source':       
        print("Symmetrising samples around m1=m2")    
        # symmetrize the samples over 1<->2
        x_temp = np.concatenate((x, y))
        y_temp = np.concatenate((y, x))        
        x=x_temp
        y=y_temp    
    else:        
        print("No reflexion done")
        
    my_range = [[xlow, xhigh], [ylow, yhigh]]
    corner.hist2d(x, y, ax=ax3, range=my_range, color=c,
                  plot_datapoints=False, plot_density=False,smooth=True,#levels=(np.exp(-0.5),np.exp(-1)),
                  levels=[0.9],no_fill_contour=False,bins=bb, lw=4,
                 )

    

#----------------------------------------------------------------------------
# e. Plotting prior (1D histograms)
#-----------------------------------------------------------------------------
       
    
def add_prior(name, parameter_1,parameter_2,zorder=10,norm_factor_x=1,norm_factor_y=1,bounded_1=False,low_1=None,high_1=None,bounded_2=False,low_2=None,high_2=None):
        
    x = samples_dict[name][parameter_1]
    y = samples_dict[name][parameter_2]
    
    xlow, xhigh = xlims
    xsmooth = np.linspace(xlow, xhigh, 10000)
    ylow, yhigh = ylims
    ysmooth = np.linspace(ylow, yhigh, 10000)

    norm_factor_x=0.95*norm_factor_x
    norm_factor_y=0.95*norm_factor_y
    
    #c = color_dict[name]
    label = 'Prior'
    alpha = 0.0
    lw = 3
    c='black'

    if bounded_1==True:        
        print("Using bounded KDE")    
        kde = Bounded_kde(x,low=low_1,high=high_1)    
    else:         
        kde = gaussian_kde(x)        
       
    ax1.plot(xsmooth, norm_factor_x*kde(xsmooth), color=c, lw=lw, label=label,
             zorder=zorder)
    ax1.fill_between(xsmooth, 0, norm_factor_x*kde(xsmooth), color=c, alpha=alpha,
                     zorder=zorder)

    if bounded_2==True:        
        print("Using bounded KDE")    
        kde = Bounded_kde(y,low=low_2,high=high_2)    
    else:         
        kde = gaussian_kde(y)    
    
    ax4.plot(norm_factor_y*kde(ysmooth), ysmooth, color='black', lw=lw, label=label,
             zorder=zorder)
    #ax4.axis('off')
    ax4.fill_betweenx(ysmooth, 0, norm_factor_y*kde(ysmooth), color=c, alpha=alpha,
                      zorder=zorder)

    my_range = [[xlow, xhigh], [ylow, yhigh]]
    

    
    
def add_prior_bounded(name, parameter_1,parameter_2,zorder=10,norm_factor_x=1,norm_factor_y=1,bounded_1=False,low_1=None,high_1=None,bounded_2=False,low_2=None,high_2=None):
    
    x = samples_dict[name][parameter_1][()]
    y = samples_dict[name][parameter_2][()]
    
    xlow, xhigh = xlims
    xsmooth = np.linspace(xlow, xhigh, 10000)
    ylow, yhigh = ylims
    ysmooth = np.linspace(ylow, yhigh, 10000)

    norm_factor_x=0.95*norm_factor_x
    norm_factor_y=0.95*norm_factor_y
    
    #c = color_dict[name]
    label = 'Prior'
    alpha = 0.3
    lw = 3
    ls = 'solid'
    kde = gaussian_kde(x)
    
    c='black'
    
    if bounded_1==True:                     
        if parameter_1 == 'chi_p':
            guess = lambda x: x * (1.-x)
            kde = gaussian_kde(x, weights=1./guess(x))
            guess_norm = integrate.quad(guess, low_1, high_1)[0]

            ax1.plot(xsmooth, norm_factor_x* kde(xsmooth) * guess(xsmooth) / guess_norm, color=c, lw=lw, label=label,
             zorder=zorder)
            
        elif parameter_1 == 'theta_jn':
            guess = lambda x: x * (np.pi-x)
            kde = gaussian_kde(x, weights=1./guess(x))
            guess_norm = integrate.quad(guess, low_1, high_1)[0]

            ax1.plot(xsmooth, norm_factor_x* kde(xsmooth) * guess(xsmooth), color=c, lw=lw, label=label,
             zorder=zorder)
                  
        else:       
            print("Using bounded KDE")   
            kde = Bounded_kde(x,low=low_1,high=high_1)          
            ax1.plot(xsmooth, norm_factor_x*kde(xsmooth), color=c, lw=lw, label=label,
                     zorder=zorder)            
    else:         
        kde = gaussian_kde(x)    
    
        ax1.plot(xsmooth, norm_factor_x*kde(xsmooth), color=c, lw=lw, label=label,
                 zorder=zorder)

    
    if bounded_2==True:        
        print("Using bounded KDE")    
        kde = Bounded_kde(y,low=low_2,high=high_2)    
    else:         
        kde = gaussian_kde(y)    
    
    ax4.plot(norm_factor_y*kde(ysmooth), ysmooth, color=c, lw=lw, label=label,
             zorder=zorder)

    
#----------------------------------------------------------------------------
# f. Plotting Ring Down modes: Phenom PHM, SEOBNR PHM (transparent contours)
#-----------------------------------------------------------------------------

def add_samples_to_fig_rd(name, parameter_1,parameter_2,zorder=10,norm_factor_x=1,norm_factor_y=1,bounded=False):
    
    x = samples_dict[name][parameter_1]
    y = samples_dict[name][parameter_2]
    xlow, xhigh = xlims
    xsmooth = np.linspace(xlow, xhigh, 1000)
    ylow, yhigh = ylims
    ysmooth = np.linspace(ylow, yhigh, 1000)

    norm_factor_x=0.95*norm_factor_x
    norm_factor_y=0.95*norm_factor_y
    
    c = color_dict[name]
    label = label_dict[name]
    alpha = alpha_dict[name]
    lw = lw_dict[name]
    #ls = ls_dict[name]
        
    kde = gaussian_kde(x)
    bb = bins_dict[name]   
    
    ax1.plot(xsmooth, norm_factor_x*kde(xsmooth), color=c, lw=lw, label=label,
             zorder=zorder)
    #ax1.axis('off')
    ax1.axvline(np.quantile(x, 0.05), color=c, ls='none')
    ax1.axvline(np.quantile(x, 0.95), color=c, ls='none')

    kde = gaussian_kde(y)
    
    ax4.plot(norm_factor_y*kde(ysmooth), ysmooth, color=c, lw=lw, label=label,
             zorder=zorder)

 #                     zorder=zorder)
    ax4.axhline(np.quantile(y, 0.05), color=c, ls='none')
    ax4.axhline(np.quantile(y, 0.95), color=c, ls='none')
    
    my_range = [[xlow, xhigh], [ylow, yhigh]]
    corner.hist2d(x, y, ax=ax3, range=my_range, color=c, ls='dashed',
                  plot_datapoints=False, plot_density=False,smooth=True,#levels=(np.exp(-0.5),np.exp(-1)),
                  levels=[0.9],no_fill_contour=False,bins=bb, lw=4, alpha=0.2
                 )
    
In [10]:
#========================
# Fig-1a: m1-m2 plot 
#========================

samples_dict = dict(IMR=phenom,NRS=nrsur,SEOB=seob)
for key in samples_dict:
    samples_dict[key] = samples_dict[key]

color_dict = dict(NRS="#7570b3", IMR="#d95f02", SEOB="#1b9e77")
label_dict = dict(NRS="NRSur PHM", IMR="Phenom PHM", SEOB="SEOBNR PHM")
alpha_dict = dict(NRS=0.2,IMR=0.2, SEOB=0.2)
lw_dict = dict(NRS=3.5, IMR=1.5, SEOB=1.5)
bins_dict = dict(NRS=50, IMR=50, SEOB=50)


parameter_1 = 'mass_1_source'
parameter_2 = 'mass_2_source'

xlims = [50, 150]
ylims = [20, 120]

norm_factor = 1.60


fig, ax1, ax3, ax4 = create_fig_and_axes(xlims, ylims)
ax3.set_xlabel(r"$m_1 \ [M_\odot]$",fontsize=24)
ax3.set_ylabel(r"$m_2 \ [M_\odot]$",fontsize=24)
ax3.tick_params(labelsize=10)

my_array_x=np.linspace(xlims[0],xlims[1],1000)
my_array_y=np.linspace(ylims[0],ylims[1],1000)

r1x=1/np.max(scipy.stats.gaussian_kde(nrsur[parameter_1])(my_array_x))
r2x=1/np.max(scipy.stats.gaussian_kde(phenom[parameter_1])(my_array_x))
r3x=1/np.max(scipy.stats.gaussian_kde(seob[parameter_1])(my_array_x))

r1y=1/np.max(scipy.stats.gaussian_kde(nrsur[parameter_2])(my_array_y))
r2y=1/np.max(scipy.stats.gaussian_kde(phenom[parameter_2])(my_array_y))
r3y=1/np.max(scipy.stats.gaussian_kde(seob[parameter_2])(my_array_y))

norm_factor_x = np.min([r1x,r2x,r3x])
norm_factor_y = np.min([r1y,r2y,r3y])

add_samples_to_fig_nofilled("SEOB",parameter_1,parameter_2,zorder=-10,norm_factor_x=norm_factor_x,norm_factor_y=norm_factor_y)
add_samples_to_fig_nofilled("IMR",parameter_1,parameter_2, zorder=10,norm_factor_x=norm_factor_x,norm_factor_y=norm_factor_y)
add_samples_to_fig("NRS",parameter_1,parameter_2,zorder=20,norm_factor_x=norm_factor_x,norm_factor_y=norm_factor_y)


ax3.set_xlim(*xlims)
ax1.set_yticklabels([],fontsize=10)
ax4.set_xticklabels([],fontsize=20)
ax3.fill_between(np.arange(50.3,119.7),np.arange(50.3,119.7),119.7,color='white',alpha=1,zorder=50)
#ax3.fill_between(np.arange(20,60),np.arange(20,60),60,color='white',alpha=1,zorder=100)

ax3.plot(np.arange(0,200),np.arange(0,200), 200, c='silver', lw=1, ls='-', alpha=0.4,zorder=100)

leg = ax3.legend(*ax4.get_legend_handles_labels(), loc=2, frameon=False, prop={'size': 20}, bbox_to_anchor=(-0.01, 1.02))
leg.set_zorder(100)
ax3.tick_params(axis='both', labelsize=20)


print("Saving")
fig.savefig(figpath+'m1m2_vf.pdf',format='pdf', transparent=True, bbox_inches='tight')
Symmetrising samples around m1=m2
Symmetrising samples around m1=m2
Symmetrising samples around m1=m2
Saving
In [11]:
#========================
# Fig-1b: chi_p-chi_eff plot 
#=======================

samples_dict = dict(IMR=phenom,NRS=nrsur,SEOB=seob, PRIOR=prior)
for key in samples_dict:
    samples_dict[key] = samples_dict[key]

color_dict = dict(NRS="#7570b3", IMR="#d95f02", SEOB="#1b9e77")
label_dict = dict(NRS="NRSur PHM", IMR="Phenom PHM", SEOB="SEOBNR PHM", PRIOR="Prior")
alpha_dict = dict(NRS=0.2,IMR=0.2, SEOB=0.2)
lw_dict = dict(NRS=3.5, IMR=1.5, SEOB=1.5)
bins_dict = dict(NRS=50, IMR=50, SEOB=50)


parameter_1 = 'chi_p'
parameter_2 = 'chi_eff'


xlims = [0, 1]
ylims = [-1, 1]

norm_factor = 1.60

fig, ax1, ax3, ax4 = create_fig_and_axes(xlims, ylims)
ax3.set_xlabel(r"$\chi_\mathrm{p}$",fontsize=24)
ax3.set_ylabel(r"$\chi_\mathrm{eff}$",fontsize=24)
ax3.tick_params(labelsize=10)

my_array_x=np.linspace(xlims[0],xlims[1],1000)
my_array_y=np.linspace(ylims[0],ylims[1],1000)

r1x=1/np.max(scipy.stats.gaussian_kde(nrsur[parameter_1])(my_array_x))
r2x=1/np.max(scipy.stats.gaussian_kde(phenom[parameter_1])(my_array_x))
r3x=1/np.max(scipy.stats.gaussian_kde(seob[parameter_1])(my_array_x))
r4x=1/np.max(scipy.stats.gaussian_kde(prior[parameter_1])(my_array_x))

r1y=1/np.max(scipy.stats.gaussian_kde(nrsur[parameter_2])(my_array_y))
r2y=1/np.max(scipy.stats.gaussian_kde(phenom[parameter_2])(my_array_y))
r3y=1/np.max(scipy.stats.gaussian_kde(seob[parameter_2])(my_array_y))
r4y=1/np.max(scipy.stats.gaussian_kde(prior[parameter_2])(my_array_x))

norm_factor_x = np.min([r1x,r2x,r3x, r4x])
norm_factor_y = np.min([r1y,r2y,r3y, r4y])

add_samples_to_fig_nofilled("SEOB",parameter_1,parameter_2,zorder=-10,norm_factor_x=norm_factor_x,norm_factor_y=norm_factor_y,bounded_1=True,low_1=0,high_1=1)
add_samples_to_fig_nofilled("IMR",parameter_1,parameter_2, zorder=10,norm_factor_x=norm_factor_x,norm_factor_y=norm_factor_y,bounded_1=True,low_1=0,high_1=1)
add_samples_to_fig("NRS",parameter_1,parameter_2,zorder=-10,norm_factor_x=norm_factor_x,norm_factor_y=norm_factor_y,bounded_1=True,low_1=0,high_1=1)
add_prior_bounded("PRIOR",parameter_1,parameter_2,zorder=-10,norm_factor_x=norm_factor_x,norm_factor_y=norm_factor_y,bounded_1=True,low_1=0,high_1=1,bounded_2=True,low_2=-1,high_2=1)


ax3.set_xlim(*xlims)
ax1.set_yticklabels([],fontsize=10)
ax4.set_xticklabels([],fontsize=20)
ax3.legend(*ax4.get_legend_handles_labels(), loc=2, frameon=False,prop={'size': 20})
ax3.tick_params(axis='both', labelsize=20)

print("Saving")
fig.savefig(figpath+'chis_vf.pdf',format='pdf', transparent=True, bbox_inches='tight')
No reflexion done
No reflexion done
No reflexion done
Using bounded KDE
Saving
In [12]:
#========================
# Fig-3: mf-af plot 
#========================

samples_dict = dict(IMR=phenom,NRS=nrsur,SEOB=seob)
for key in samples_dict:
    samples_dict[key] = samples_dict[key]

color_dict = dict(NRS="#7570b3", IMR="#d95f02", SEOB="#1b9e77")
label_dict = dict(NRS="NRSur PHM", IMR="Phenom PHM", SEOB="SEOBNR PHM")
alpha_dict = dict(NRS=0.2,IMR=0.2, SEOB=0.2)
lw_dict = dict(NRS=3.5, IMR=1.5, SEOB=1.5)
bins_dict = dict(NRS=50, IMR=50, SEOB=50)


parameter_1 = 'final_mass_source'
parameter_2 = 'final_spin'

xlims = [100, 220]
ylims = [0.25, 1]

norm_factor = 1.60


fig, ax1, ax3, ax4 = create_fig_and_axes(xlims, ylims)
ax3.set_xlabel(r"$M_\mathrm{f} \ [M_\odot]$",fontsize=24)
ax3.set_ylabel(r"$\chi_\mathrm{f} $",fontsize=24)
ax3.tick_params(labelsize=10)

my_array_x=np.linspace(xlims[0],xlims[1],1000)
my_array_y=np.linspace(ylims[0],ylims[1],1000)

r1x=1/np.max(scipy.stats.gaussian_kde(nrsur[parameter_1])(my_array_x))
r2x=1/np.max(scipy.stats.gaussian_kde(phenom[parameter_1])(my_array_x))
r3x=1/np.max(scipy.stats.gaussian_kde(seob[parameter_1])(my_array_x))

r1y=1/np.max(scipy.stats.gaussian_kde(nrsur[parameter_2])(my_array_y))
r2y=1/np.max(scipy.stats.gaussian_kde(phenom[parameter_2])(my_array_y))
r3y=1/np.max(scipy.stats.gaussian_kde(seob[parameter_2])(my_array_y))

norm_factor_x = np.min([r1x,r2x,r3x])
norm_factor_y = np.min([r1y,r2y,r3y])

add_samples_to_fig_nofilled("SEOB",parameter_1,parameter_2,zorder=-10,norm_factor_x=norm_factor_x,norm_factor_y=norm_factor_y,bounded_2=True,low_2=0.,high_2=1)
add_samples_to_fig_nofilled("IMR",parameter_1,parameter_2, zorder=10,norm_factor_x=norm_factor_x,norm_factor_y=norm_factor_y,bounded_2=True,low_2=0.,high_2=1)
add_samples_to_fig("NRS",parameter_1,parameter_2,zorder=20,norm_factor_x=norm_factor_x,norm_factor_y=norm_factor_y,bounded_2=True,low_2=0.,high_2=1)

ax3.set_xlim(*xlims)
ax1.set_yticklabels([],fontsize=10)
ax4.set_xticklabels([],fontsize=20)
ax3.tick_params(axis='both', labelsize=20)
ax3.legend(*ax4.get_legend_handles_labels(), loc=3, frameon=False,prop={'size': 20})


print("Saving")
fig.savefig(figpath+'mfaf_vf.pdf',format='pdf', transparent=True, bbox_inches='tight')
Using bounded KDE
No reflexion done
Using bounded KDE
No reflexion done
Using bounded KDE
No reflexion done
Saving
In [13]:
#========================
# Fig-4: kick_mag_vf plot 
#========================
vf = nrsur['final_kick']
vf_prior = prior['final_kick']

label_fontsize = 14
legend_fontsize = 11

galaxy_colors = [u'#1b9e77', u'#c7e9b4', u'#d95f02', u'#253494']
post_prior_colors = ['#7570b3', 'k']


# ---------------------------------------------------------------------------
def set_y_ticks(ax, yticks_str):
    ax.set_yticklabels(yticks_str)
    ax.set_yticks([float(tmp) for tmp in yticks_str])

# ---------------------------------------------------------------------------
def plotter(ax, data, data_prior, xlabel, color=None):

    x_dense = np.linspace(0, 3500, 1000)
    data_kde = Bounded_kde(data, low=0.0)(x_dense)
    data_prior_kde = Bounded_kde(data_prior, low=0.0)(x_dense)
#    data_kde = gaussian_kde(data, bw_method=0.2)(x_dense)
#    data_prior_kde = gaussian_kde(data_prior, bw_method=0.2)(x_dense)
    #norm_factor = 0.95/np.max(data_kde)

    #ax.hist(data, density=True, histtype='step', bins=25, \
    #    color=post_prior_colors[0], label='Posterior', lw=1.5)
    plt.plot(x_dense, data_kde, color=post_prior_colors[0], \
        label='NRSur PHM', lw=1.5)

    #ax.hist(data_prior, density=True, histtype='step', \
    #    bins=25, color=post_prior_colors[1], label='Prior', lw=1.5)
    plt.plot(x_dense, data_prior_kde, color=post_prior_colors[1], \
        label='Prior', lw=1.5, zorder=-1)

    ax.set_xlabel(xlabel, fontsize=label_fontsize)
    ax.set_ylabel("Probability density", fontsize=label_fontsize)

# ---------------------------------------------------------------------------
def plot_rect(ax, yval, vmin, vmax, color, label, height=0.3):
    rect = patches.Rectangle((vmin, yval), vmax-vmin, height, linewidth=1, \
        edgecolor=None, facecolor=color, alpha=0.8, label=label)
    ax.add_patch(rect)

fig = plt.figure(figsize=(5, 3.5))
gs = gridspec.GridSpec(2, 1, height_ratios=[1, 4])
plt.subplots_adjust(wspace=0, hspace=0)

ax = plt.subplot(gs[1])
plotter(ax, vf, vf_prior, '$v_{\mathrm{f}}$ [km/s]')
ax.set_ylim(0, 1.7e-3)
ax.set_xlim(0, 3500)

# These four are dummies, so we get them into the legend
plot_rect(ax, 0.4, 0, 150, galaxy_colors[0], 'Globular clusters')
plot_rect(ax, 0.05, 10, 1000, galaxy_colors[1], 'Nuclear star clusters')
plot_rect(ax, 0.4, 400, 2500, galaxy_colors[2], 'Elliptical galaxies')
ax.scatter(580, 0.85, label='Milky Way', marker='*', s=120, zorder=10, \
    color=galaxy_colors[3])

ax.legend(loc='upper right', fontsize=legend_fontsize, frameon=0)

ax = plt.subplot(gs[0])
plot_rect(ax, 0.4, 0, 150, galaxy_colors[0], 'Globular clusters')
plot_rect(ax, 0.05, 10, 1000, galaxy_colors[1], 'Nuclear star clusters')
plot_rect(ax, 0.4, 400, 2500, galaxy_colors[2], 'Elliptical galaxies')
ax.scatter(580, 0.85, label='Milky Way', marker='*', s=120, zorder=10, \
    color=galaxy_colors[3])
ax.set_ylim(0, 1)
ax.set_xlim(0, 3500)
ax.axis('off')

print("Saving")
fig.savefig(figpath+'kick_mag_vf_h5.pdf',format='pdf', transparent=True, bbox_inches='tight')
Saving
In [14]:
#========================
# Fig-6: theta_JN - D_L plot 
# In this plot, theta_JN is in radians, but elsewhere, it is in degrees.
#=======================

samples_dict = dict(IMR=phenom,NRS=nrsur,SEOB=seob, PRIOR=prior)
for key in samples_dict:
    samples_dict[key] = samples_dict[key]

color_dict = dict(NRS="#7570b3", IMR="#d95f02", SEOB="#1b9e77")
label_dict = dict(NRS="NRSur PHM", IMR="Phenom PHM", SEOB="SEOBNR PHM", PRIOR="Prior")
alpha_dict = dict(NRS=0.2,IMR=0.2, SEOB=0.2)
lw_dict = dict(NRS=3.5, IMR=1.5, SEOB=1.5)
bins_dict = dict(NRS=50, IMR=50, SEOB=50)


parameter_1 = 'theta_jn'
parameter_2 = 'luminosity_distance'

xlims = [0, 3.14159]
ylims = [1000, 10000]

norm_factor = 1.60


fig, ax1, ax3, ax4 = create_fig_and_axes(xlims, ylims)
ax3.set_xlabel(r"$\theta_\mathrm{JN}\ [\mathrm{radians}]$",fontsize=24) 
ax3.set_ylabel(r"$D_L \ [\mathrm{Mpc}]$",fontsize=24)
ax3.tick_params(labelsize=10)

my_array_x=np.linspace(xlims[0],xlims[1],1000)
my_array_y=np.linspace(ylims[0],ylims[1],1000)

r1x=1/np.max(scipy.stats.gaussian_kde(nrsur[parameter_1])(my_array_x))
r2x=1/np.max(scipy.stats.gaussian_kde(phenom[parameter_1])(my_array_x))
r3x=1/np.max(scipy.stats.gaussian_kde(seob[parameter_1])(my_array_x))
r4x=1/np.max(scipy.stats.gaussian_kde(prior[parameter_1])(my_array_x))

r1y=1/np.max(scipy.stats.gaussian_kde(nrsur[parameter_2])(my_array_y))
r2y=1/np.max(scipy.stats.gaussian_kde(phenom[parameter_2])(my_array_y))
r3y=1/np.max(scipy.stats.gaussian_kde(seob[parameter_2])(my_array_y))
r4y=1/np.max(scipy.stats.gaussian_kde(prior[parameter_2])(my_array_x))

norm_factor_x = np.min([r1x,r2x,r3x, r4x])
norm_factor_y = np.min([r1y,r2y,r3y, r4y])


add_samples_to_fig_nofilled("SEOB",parameter_1,parameter_2, zorder=-10,norm_factor_x=norm_factor_x,norm_factor_y=norm_factor_y,bounded_1=True,low_1=0,high_1=np.pi)
add_samples_to_fig_nofilled("IMR",parameter_1,parameter_2, zorder=10,norm_factor_x=norm_factor_x,norm_factor_y=norm_factor_y,bounded_1=True,low_1=0,high_1=np.pi)
add_samples_to_fig("NRS",parameter_1,parameter_2,zorder=20,norm_factor_x=norm_factor_x,norm_factor_y=norm_factor_y,bounded_1=True,low_1=0,high_1=np.pi)
#add_prior("PRIOR",parameter_1,parameter_2,zorder=-10,norm_factor_x=norm_factor_x,norm_factor_y=norm_factor_y,bounded_1=True,low_1=0,high_1=np.pi)
add_prior_bounded("PRIOR",parameter_1,parameter_2,zorder=-10,norm_factor_x=norm_factor_x,norm_factor_y=norm_factor_y,bounded_1=True,low_1=0,high_1=np.pi)


ax3.set_xlim(*xlims)
ax1.set_yticklabels([],fontsize=10)
ax4.set_xticklabels([],fontsize=20)
ax3.legend(*ax4.get_legend_handles_labels(), loc=9, frameon=False,prop={'size': 20})
ax3.tick_params(axis='both', labelsize=20)

print("Saving")
fig.savefig(figpath+'dtheta_vf.pdf',format='pdf', transparent=True, bbox_inches='tight')
No reflexion done
No reflexion done
No reflexion done
Saving
In [15]:
#========================
# Fig-7: m1 - z plot for HMs 
#=======================

samples_dict = dict(l2=nrs_l2,l3=nrs_l3,l4=nrs_l4)
for key in samples_dict:
    samples_dict[key] = samples_dict[key]

color_dict = dict(l4="#253494", l3="#41b6c4", l2="#c7e9b4")

label_dict = dict(l2=r"$\ell = 2$", l3=r"$\ell \leq 3$", l4=r"$\ell \leq 4$")
alpha_dict = dict(l4=0.2,l3=0.2, l2=0.2)
lw_dict = dict(l4=3.5, l3=1.5, l2=1.5)
bins_dict = dict(l4=50, l3=50, l2=50)


parameter_1 = 'mass_1_source'
parameter_2 = 'redshift'

xlims = [60, 160]
ylims = [0.1, 1.3]

norm_factor = 1.60


fig, ax1, ax3, ax4 = create_fig_and_axes(xlims, ylims)
ax3.set_xlabel(r"$m_1 \ [M_\odot]$",fontsize=24)
ax3.set_ylabel(r"Redshift $z$",fontsize=24)
ax3.tick_params(labelsize=10)

my_array_x=np.linspace(xlims[0],xlims[1],1000)
my_array_y=np.linspace(ylims[0],ylims[1],1000)

r1x=1/np.max(scipy.stats.gaussian_kde(nrs_l2[parameter_1])(my_array_x))
r2x=1/np.max(scipy.stats.gaussian_kde(nrs_l2[parameter_1])(my_array_x))
r3x=1/np.max(scipy.stats.gaussian_kde(nrs_l4[parameter_1])(my_array_x))

r1y=1/np.max(scipy.stats.gaussian_kde(nrs_l2[parameter_2])(my_array_y))
r2y=1/np.max(scipy.stats.gaussian_kde(nrs_l3[parameter_2])(my_array_y))
r3y=1/np.max(scipy.stats.gaussian_kde(nrs_l4[parameter_2])(my_array_y))

norm_factor_x = np.min([r1x,r2x,r3x])
norm_factor_y = np.min([r1y,r2y,r3y])

add_samples_to_fig("l2",parameter_1,parameter_2,zorder=-10,norm_factor_x=norm_factor_x,norm_factor_y=norm_factor_y)
add_samples_to_fig("l3",parameter_1,parameter_2, zorder=10,norm_factor_x=norm_factor_x,norm_factor_y=norm_factor_y)
add_samples_to_fig("l4",parameter_1,parameter_2,zorder=20,norm_factor_x=norm_factor_x,norm_factor_y=norm_factor_y)

ax3.set_xlim(*xlims)
ax1.set_yticklabels([],fontsize=10)
ax4.set_xticklabels([],fontsize=20)
ax3.tick_params(axis='both', labelsize=20)
ax3.legend(*ax4.get_legend_handles_labels(), loc=0, frameon=False,prop={'size': 20})


print("Saving")
fig.savefig(figpath+'m1HM_vf.pdf',format='pdf', transparent=True, bbox_inches='tight')
No reflexion done
No reflexion done
No reflexion done
Saving