Codes for making plots for O3 anisotropic stochastic paper¶
In [23]:
import numpy as np
import healpy as hp
import h5py
import scipy.io as sio
import matplotlib # requires matplotlib version >= 3.2.2
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
from scipy.io import loadmat
import matplotlib.patches as mpatches
from matplotlib.lines import Line2D
import warnings
warnings.filterwarnings("ignore", category=UserWarning)
Parameters for plotting¶
In [24]:
fontsize = 24
matplotlib.rcParams.update({
"font.size": fontsize,
"axes.titlesize": fontsize,
"axes.labelsize": fontsize,
"xtick.labelsize": fontsize,
"ytick.labelsize": fontsize,
"xtick.major.size": fontsize * 0.8,
"ytick.major.size": fontsize * 0.8,
"legend.fontsize": fontsize * 0.9,
"font.family": "Times new Roman",
"figure.dpi": 100,
"savefig.dpi": 300,
"text.usetex": True,
"path.simplify": True,
"figure.figsize": (8, 6)
})
# Class to deal in an easy way with number of digits in colorbars (used for ULs and sigma maps)
# It allows to choose the number of digits for tickslabels and keeping 10^N at the right of the colorbar
# instead of having 1eN for each tickslabel.
class FormatScalarFormatter(matplotlib.ticker.ScalarFormatter):
def __init__(self, fformat="%1.1f", offset=True, mathText=True):
self.fformat = fformat
matplotlib.ticker.ScalarFormatter.__init__(self,useOffset=offset,
useMathText=mathText)
def _set_format(self):
self.format = self.fformat
if self._useMathText:
#self.format = '$%s$' % matplotlib.ticker._mathdefault(self.format)
self.format = '$%s$' % ('\\mathdefault{%s}' % self.format)
Variable definitions¶
In [25]:
# variable corresponding to spectral shapes of signal models
alphas = ["0", "2over3", "3"]
Figure 2: Top row (BBR SNR maps)¶
In [26]:
# SNR sky-maps for BBR analysis
# Create figure
fig = plt.figure(figsize=(24, 8))
# Limits for the colorbars
vlims = [2.55, 3.72, 3.63]
# Loop over the *combined_maps.hdf5 files for the three spectral indices 0, 2/3 and 3
for alpha, vlim, index in zip (alphas, vlims, range(0,len(alphas))):
# SNR *.hdf5 file
snr_file = "./datafiles/bbr_a" + alpha + "_o1_o2_o3_combined_maps.hdf5"
# Read the file and create a dictionary-like variable; the file has three variables SNR (map_snr_total), signal estimate (ptEst_total) and its sigma (sig_total)
combined_snr = h5py.File(snr_file, "r")
# Select the map of interest (SNR map)
snr_map = np.array(combined_snr["map_snr_total"])
print("alpha:%s, max snr:%.1f" % (alpha, np.max(snr_map)))
# Close the file
combined_snr.close()
# Plot the sky-map for a given alpha as subplot using healpy Mollweide projection
hp.mollview(np.real(np.squeeze(snr_map)),rot=(180,0,0),flip='astro', title="",nest=False, format = "%.2f", cbar=False, sub=(1, 3, index+1), margins=[0.0125,0.0125, 0.0125, 0.0125], min=-vlim, max=vlim)
# Add graticule
hp.graticule(45,90)
# Add useful annotations on the map for right ascension and declination
plt.annotate('12h',(0.46,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
plt.annotate('6h',(0.725,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
plt.annotate('18h',(0.21,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
plt.annotate(r'$45^{\circ}$',(0.03,0.8),xycoords='axes fraction', fontsize=fontsize)
plt.annotate(r'$0^{\circ}$',(-0.045,0.47),xycoords='axes fraction', fontsize=fontsize)
plt.annotate(r'$-45^{\circ}$',(-0.002,0.14),xycoords='axes fraction', fontsize=fontsize)
# Trick to make colorbars all the same (healpy colorbar does not leave much freedom of customization and it is not used)
fig = plt.gcf()
ax = plt.gca()
image = ax.get_images()[0]
cmap = fig.colorbar(image, ax=ax, shrink=0.7, orientation="horizontal", pad=0.15, aspect = 40, ticks=[-2,0,2])
plt.show()
Figure 2: Bottom row (BBR Upper limit maps)¶
In [27]:
# Upper limit sky-maps for BBR analysis
fig = plt.figure(figsize=(24, 8))
# Loop over the upper limit *.mat files for the three spectral indices 0, 2/3 and 3
for alpha, index in zip (alphas, range(0,len(alphas))):
# upper limit *.mat file
ul_dir = "./datafiles/bbr_a"+ alpha +"_zl_limits.mat"
# Read the file and create a dictionary-like variable; the file has variables correponding to upplimits from O1+O2, just O3 (for three baselines) and O1+O2+O3
combined_ul = sio.loadmat(ul_dir)
ul_map = np.squeeze(np.array(combined_ul["flux_ul_combined"]))
print("alpha:%s, min UL:%.1e, max UL:%.1e" % (alpha, np.min(ul_map), np.max(ul_map)))
# Plot the sky-map for a given alpha as subplot using healpy Mollweide projection
hp.mollview(np.real(np.squeeze(ul_map)),rot=(180,0,0),flip='astro', title="",nest=False, format = "%.2g", cbar=False, sub=(1, 3, index+1), margins=[0.0125,0.0125, 0.0125, 0.0125])
# Add graticule
hp.graticule(45,90)
# Add useful annotations on the map for right ascension and declination
plt.annotate('12h',(0.46,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
plt.annotate('6h',(0.725,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
plt.annotate('18h',(0.21,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
plt.annotate(r'$45^{\circ}$',(0.03,0.8),xycoords='axes fraction', fontsize=fontsize)
plt.annotate(r'$0^{\circ}$',(-0.045,0.47),xycoords='axes fraction', fontsize=fontsize)
plt.annotate(r'$-45^{\circ}$',(-0.002,0.14),xycoords='axes fraction', fontsize=fontsize)
fig = plt.gcf()
ax = plt.gca()
image = ax.get_images()[0]
# Further colorbar customization
nticks = 4
interval = 1./(nticks-1)
diff= np.max(ul_map) -np.min(ul_map)
ticks_no_round = [np.min(ul_map)+(diff)*interval*x for x in np.arange(0,nticks)]
cmap = fig.colorbar(image, ax=ax, ticks=ticks_no_round, format=FormatScalarFormatter("%.1f"), shrink=0.7, orientation="horizontal", pad=0.15, aspect = 30)
# Adding units as colorbar label
cmap.set_label(r'[erg cm$^{-2}$ Hz$^{-1}$ s$^{-1}$]', fontsize = 1.2*fontsize, labelpad = 20)
cmap.ax.tick_params(labelsize=fontsize)
plt.show()
Figure 3: Top row (SHD SNR maps)¶
In [28]:
# Limits for the colorbars
vlims = [2.18, 3.72, 3.57]
fig = plt.figure(figsize=(24,8))
# Loop over the shd *.mat files for the three spectral indices 0, 2/3 and 3
for i, (alpha, vlim) in enumerate(zip(alphas, vlims)):
ax = fig.add_subplot(130+(i+1), projection="mollweide")
pproc_file = "./datafiles/shd_ul_a"+ alpha + "_o1_o2_o3_shd_map_files.mat"
shd_data = loadmat(pproc_file)
# Get the signal estimate (map0) and its sigma (sigmaPix0) for different pixels on the sky
p0_map = shd_data["map0"]
sigma_map = shd_data["sigmaPix0"]
sig_map = sigma_map.reshape(360,181).T
# Calculate SNR maps
snr_map = (p0_map / sigma_map).reshape(360, 181).T
print("alpha:%s, max snr:%.1f" % (alpha, np.max(snr_map)))
# Make the X, Y mesh and plot
x = np.linspace(-np.pi,np.pi,360)
y = np.linspace(-np.pi/2,np.pi/2,181)
X,Y = np.meshgrid(x,y)
image = ax.pcolormesh(X,Y,snr_map[:, ::-1], vmin=-vlim, vmax=vlim)
# formate x and y tickes
raval = [-np.pi/2,0,np.pi/2]
ralab = ['18h','12h','6h']
ax.set_xticks(raval)
xtick_labels = ax.set_xticklabels(ralab, color="w",fontweight='bold',fontsize=20)
declval = [-np.pi/4,0,np.pi/4]
decllab = [r'$-45^{\circ}$',r'$0^{\circ}$',r'$45^{\circ}$']
ax.set_yticks(declval)
ax.set_yticklabels(decllab,fontsize=20)
plt.grid(True, color="k", linestyle=":")
x = np.array([-np.radians(180), np.radians(180)])
y = np.array([0, 0])
plt.plot(x, y, ls='-', color="k", lw=0.8)
cmp = plt.colorbar(image, shrink=0.7, ax=ax, orientation="horizontal", aspect=40)
cmap.set_label('', fontsize = 10)
for t in cmap.ax.get_xticklabels(): # reduce the font size of xlabel
t.set_fontsize(10)
plt.show()
Figure 3: bottom row (SHD UL maps)¶
In [29]:
fig = plt.figure(figsize=(24,8))
# Loop over the upper limit *.mat files for the three spectral indices 0, 2/3 and 3
for i, alpha in enumerate(alphas):
ax = fig.add_subplot(130+(i+1), projection="mollweide")
# upper limit *.mat file
pproc_file = "./datafiles/shd_a" + alpha + "_limits_o1_o2_o3_HLV_ZL.mat"
shd_data = loadmat(pproc_file)
omega_ul_map = shd_data["omega_ul"].reshape(360, 181).T
print("alpha:%s, min UL:%.1e, max UL:%.1e" % (alpha, np.min(omega_ul_map), np.max(omega_ul_map)))
# Make the X, Y mesh and plot
x = np.linspace(-np.pi,np.pi,360)
y = np.linspace(-np.pi/2,np.pi/2,181)
X,Y = np.meshgrid(x,y)
image = ax.pcolormesh(X,Y,omega_ul_map[:, ::-1])
# formate x and y tickes
raval = [-np.pi/2,0,np.pi/2]
ralab = ['18h','12h','6h']
ax.set_xticks(raval)
xtick_labels = ax.set_xticklabels(ralab, color="w",fontweight='bold',fontsize=20)
declval = [-np.pi/4,0,np.pi/4]
decllab = ['$-45^{\circ}$','$0^{\circ}$','$45^{\circ}$']
ax.set_yticks(declval)
ax.set_yticklabels(decllab,fontsize=20)
plt.grid(True, color="k", linestyle=":")
x = np.array([-np.radians(180), np.radians(180)])
y = np.array([0, 0])
plt.plot(x, y, ls='-', color="k", lw=0.8)
# Further colorbar customization
nticks = 4
interval = 1./(nticks-1)
diff= np.max(omega_ul_map) -np.min(omega_ul_map)
ticks_no_round = [np.min(omega_ul_map)+(diff)*interval*x for x in np.arange(0,nticks)]
cmap = fig.colorbar(image, ax=ax, ticks=ticks_no_round, format=FormatScalarFormatter("%.1f", mathText=False), shrink=0.7, orientation="horizontal", pad=0.15, aspect = 40)
# Adding units as colorbar label
cmap.set_label(r'$[\mathrm{sr}^{-1}]$', fontsize = fontsize,labelpad=20.)
for t in cmap.ax.get_xticklabels():
t.set_fontsize(20)
plt.show()
Figure 4: C_l upper limits¶
In [30]:
H0 = 2.2005e-18 # Hubble constant
fref = 25.0 # reference frequency
const = 2.0 * np.pi**2 * fref**3 / (3 * H0**2) # normalization constant
colors = ['xkcd:black', 'xkcd:sky blue', 'xkcd:orange']
fig = plt.figure()
ax = fig.add_subplot(111)
# Loop over the upper limit *.mat files for the three spectral indices 0, 2/3 and 3
for alpha, color in zip(alphas, colors):
pproc_file = "./datafiles/shd_ul_a" + alpha + "_o1_o2_o3_shd_map_files.mat"
shd_data = loadmat(pproc_file)
lmax = shd_data['L']
Cell_ul = shd_data["omega_cl_ul"]
Cell = shd_data['omega_cl']
ax.semilogy(np.arange(lmax + 1), np.sqrt(Cell_ul)[0], color=color, linestyle="None", marker="v", markersize = 12, label=r"$\alpha=%s$" % alpha.replace("over", "/"))
ax.grid(which="major", linestyle="-")
ax.grid(which="minor", linestyle=":")
ax.legend(loc="upper right")
ax.set_xticks(np.arange(0, 17, 2))
ax.set_yticks([1e-8, 1e-9, 1e-10])
ax.set_ylim((0.8e-10, 0.2e-7))
ax.set_xlabel("$\ell$")
ax.set_ylabel("$C_\ell^{1/2}[\mathrm{sr}^{-1}]$")
plt.show()
Figure 5: Fisher Matrix comparison¶
In [31]:
baselines = ["HL", "HLV"]
colors = ['xkcd:black', 'xkcd:sky blue', 'xkcd:orange']
fkeep = 2/3 # keep 2/3 of dominant eigen values
yvals = []
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111)
# Loop over the *.mat files for the three spectral indices 0, 2/3 and 3; for HL and HLV network
for i, (alpha, color) in enumerate(zip(alphas, colors)):
for ii, (baseline, style) in enumerate(zip(baselines, ["dotted", "solid"])):
matfile = "./datafiles/shd_ul_a" + alpha + "_shd_map_files_" + baseline + ".mat"
Fisher = loadmat(matfile)["fisher0"]
eigs, _ = np.linalg.eigh(Fisher) # calculate eigen values
eigs = eigs[::-1]
yvals.append(eigs/eigs[0])
if alpha == '2over3':
lab = r'$\alpha=2/3$' + ', ' + baseline
else:
lab = r'$\alpha=$%s'%alpha + ', ' + baseline
ax.loglog(np.arange(len(eigs)), eigs/eigs[0], color=color, linestyle=style, label=lab)
ax.axvline(len(eigs) * fkeep, linestyle="--", color=color)
_, _ = ax.get_legend_handles_labels()
legend_elements1 = [mpatches.Patch(color=color, label=r"$\alpha=%s$" % alpha.replace("over", "/")) for color, alpha in zip(colors, alphas)]
legend_elements2 = [Line2D([0], [0], color="k", lw=2, label=label, linestyle=ls) for ls, label in zip(["-", ":"], ["HLV", "HL"])]
leg1 = ax.legend(handles=legend_elements1, loc='lower left', frameon=True, fontsize=18)
leg2 = ax.legend(handles=legend_elements2, loc=(0.03, 0.3), frameon=True, fontsize=18)
ax.add_artist(leg1)
for line in leg2.get_lines():
line.set_linewidth(2)
ax.set_xlabel("Mode index",fontsize=20)
ax.set_ylabel("Eigenvalues",fontsize=20)
plt.show()
Figure 6: NBR plots¶
In [32]:
# function to NBR plot upper limits for the three directions of interest (Sco X-1, SN 1987a, Galactic Center )
def plotUL(point):
# point = {scox1,sn1987a,gc}
nbr_file = "./datafiles/NBR_" + point + "_ul.hdf5"
struct = h5py.File(nbr_file, "r") # load the corresponding file
freqs = np.array(struct["f_all"])
ul = np.array(struct["ul_" + point])
onesigma = np.array(struct["onesigma_" + point])
plt.xlim([np.amin(freqs),np.amax(freqs)]);
plt.ylim([5e-26,1e-22]);
plt.grid(which = "both", lw = 0.2, ls = "--")
plt.loglog(freqs, ul, c = "dimgrey", lw = 0.3, label = "HLV: O1+O2+O3: 95 % UL" )
plt.loglog(freqs, onesigma, c = "k", lw = 0.5, label = "$1\sigma$ sensitivity")
plt.ylabel("Strain amplitude ($h_0$)")
plt.xlabel("Frequency [Hz]")
plt.legend()
if point == "scox1":
plt.title("Sco X-1")
plt.savefig('ul_and_onesig-scox1-HLV.png',bbox_inches='tight')
if point == "sn1987a":
plt.title("SN 1987a")
plt.savefig('ul_and_onesig-sn-HLV.png',bbox_inches='tight')
if point == "gc":
plt.title("Galactic Center")
plt.savefig('ul_and_onesig-gc-HLV.png',bbox_inches='tight')
return "Plotting " + point
In [33]:
plotUL("scox1")
Out[33]:
In [34]:
plotUL("sn1987a")
Out[34]:
In [35]:
plotUL("gc")
Out[35]:
Figure 7: Sigma maps for single baselines (BBR)¶
In [36]:
# HL, HV, LV 1-sigma maps
# Normalization constant such that the 1-sigma quantity corresponds to GW flux
G = 6.6726e-8
c = 2.9979e10
fref = 25.0
flux_factor = (c**3) *np.pi*(fref**2)/(4*G)
# Same as for ULs maps, but with 9 subplots of 1-sigma sensitivity for three baselines and three spectral indices
fig = plt.figure(figsize=(24,18))
baselines = ["HL", "HV", "LV"]
for ii, baseline in enumerate(baselines):
print("Baseline:%s " % (baseline))
for jj, alpha in enumerate (alphas):
map_dir = "./datafiles/bbr_a{0}_{1}_Map_dirty32.hdf5".format(alpha, baseline)
maps = h5py.File(map_dir, "r") # load file
print("alpha:%s, max snr:%.1f" % (alpha, np.max(np.array(maps['map_snr_total']))))
map_sigma = flux_factor*np.array(maps['sig_total']) # express 1-sigma in GW flux
hp.mollview(np.real(np.squeeze(map_sigma)),rot=(180,0,0),flip='astro', title="",nest=False, format = "%.2f", cbar=False, sub=(3, 3, 3*ii+ jj+1), margins=[0.0125,0.0125, 0.0125, 0.0125])
hp.graticule(45,90)
plt.annotate('12h',(0.46,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
plt.annotate('6h',(0.725,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
plt.annotate('18h',(0.21,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
plt.annotate(r'$45^{\circ}$',(0.03,0.8),xycoords='axes fraction', fontsize=fontsize)
plt.annotate(r'$0^{\circ}$',(-0.045,0.47),xycoords='axes fraction', fontsize=fontsize)
plt.annotate(r'$-45^{\circ}$',(-0.002,0.14),xycoords='axes fraction', fontsize=fontsize)
if alpha=='0':
if baseline=='HL':
plt.annotate('$\\bf{HL}$',(-0.18,0.45),xycoords='axes fraction',fontsize=30)
if baseline=='HV':
plt.annotate('$\\bf{HV}$',(-0.18,0.45),xycoords='axes fraction',fontsize=30)
if baseline=='LV':
plt.annotate('$\\bf{LV}$',(-0.18,0.45),xycoords='axes fraction',fontsize=30)
fig = plt.gcf()
ax = plt.gca()
image = ax.get_images()[0]
nticks = 4
interval = 1./(nticks-1)
diff= np.max(map_sigma) -np.min(map_sigma)
ticks_no_round = [np.min(map_sigma)+(diff)*interval*x for x in np.arange(0,nticks)]
fmt = FormatScalarFormatter("%.1f")
cmap = fig.colorbar(image, ax=ax, ticks=ticks_no_round, format=fmt, shrink=0.7, orientation="horizontal", pad=0.15, aspect = 40)
cmap.set_label(r'[erg cm$^{-2}$ Hz$^{-1}$ s$^{-1}$]', fontsize = 1.2*fontsize, labelpad = 20)
plt.show()
Figure 8: Sigma maps for single baselines (SHD)¶
In [37]:
baselines = {'HL','HV','LV'}
fig = plt.figure(figsize=(24,18))
# Same as for ULs maps, but with 9 subplots of 1-sigma sensitivity for three baselines and three spectral indices
for ii, baseline in enumerate(baselines):
print("Baseline:%s " % (baseline))
for i, alpha in enumerate(alphas):
ax = fig.add_subplot(330+(i+1)+3*ii, projection="mollweide")
pproc_file = "./datafiles/shd_ul_a" + alpha + "_shd_map_files_" + baseline + ".mat"
shd_data = loadmat(pproc_file) # load data file
p0_map = shd_data["map0"]
sigma_map = shd_data["sigmaPix0"]
sig_map = sigma_map.reshape(360,181).T * const # express 1-sigma as a normalized quantity
snr_map = (p0_map / sigma_map).reshape(360, 181).T
print("alpha:%s, max snr:%.1f" % (alpha, np.max(snr_map)))
x = np.linspace(-np.pi,np.pi,360)
y = np.linspace(-np.pi/2,np.pi/2,181)
X,Y = np.meshgrid(x,y)
# make plot
image = ax.pcolormesh(X,Y,sig_map[:, ::-1], vmin=np.min(sig_map), vmax=np.max(sig_map))
raval = [-np.pi/2,0,np.pi/2]
ralab = ['18h','12h','6h']
ax.set_xticks(raval)
xtick_labels = ax.set_xticklabels(ralab, color="w",fontsize=20)
declval = [-np.pi/4,0,np.pi/4]
decllab = ['$-45^{\circ}$','$0^{\circ}$','$45^{\circ}$']
ax.set_yticks(declval)
ax.set_yticklabels(decllab,fontsize=20)
plt.grid(True, color="k", linestyle=":")
x = np.array([-np.radians(180), np.radians(180)])
y = np.array([0, 0])
plt.plot(x, y, ls='-', color="k", lw=0.8)
nticks = 4
interval = 1./(nticks-1)
diff= np.max(sig_map) -np.min(sig_map)
ticks_no_round = [np.min(sig_map)+(diff)*interval*x for x in np.arange(0,nticks)]
cmap = fig.colorbar(image, ax=ax, ticks=ticks_no_round, format=FormatScalarFormatter("%.1f", mathText=False), shrink=0.7, orientation="horizontal", pad=0.15, aspect = 40)
cmap.set_label(r'$[\mathrm{sr}^{-1}]$', fontsize = fontsize, labelpad=20.)
for t in cmap.ax.get_xticklabels():
t.set_fontsize(20)
if alpha=='0':
if baseline=='HL':
plt.annotate('$\\bf{HL}$',(-0.18,0.45),xycoords='axes fraction',fontsize=24)
if baseline=='HV':
plt.annotate('$\\bf{HV}$',(-0.18,0.45),xycoords='axes fraction',fontsize=24)
if baseline=='LV':
plt.annotate('$\\bf{LV}$',(-0.18,0.45),xycoords='axes fraction',fontsize=24)
plt.show()
Figure 9: Sigma maps for single baselines (NBR)¶
In [38]:
# function to NBR 1-sigma plots for the three directions of interest (Sco X-1, SN 1987a, Galactic Center) and three baselines
def plotSigma(point):
# point = {scox1,sn1987a,gc}
maps_mat_output = h5py.File("./datafiles/NBR_" + point + '_sigma_HLV.hdf5', "r")
tempstruct = h5py.File("./datafiles/NBR_scox1_ul.hdf5", "r")
f_all = np.array(tempstruct["f_all"])
plt.loglog(f_all,np.array(maps_mat_output["sigma_HL"]), label='HL', lw = 2,color='black',ls = '-')
plt.loglog(f_all, np.array(maps_mat_output["sigma_HV"]), label='HV',lw = 2,color='skyblue',ls = '-.')
plt.loglog(f_all, np.array(maps_mat_output["sigma_LV"]), label='LV',lw = 2,color='orange',ls=':')
plt.legend(bbox_to_anchor=(1,1), loc="upper right")
plt.grid(which = "both", lw = 0.2, ls = "--")
plt.xlabel('Frequency [Hz]')
plt.ylabel('Sigma(f)')
if point == "scox1":
plt.title("Sco X-1")
if point == "sn1987a":
plt.title("SN 1987a")
if point == "gc":
plt.title("Galactic Center")
plt.tight_layout()
plt.savefig('NBR_sigma_' + point + '.png',bbox_inches='tight')
plt.show()
In [39]:
plotSigma("scox1")
In [40]:
plotSigma("sn1987a")
In [41]:
plotSigma("gc")
