In [3]:
% matplotlib inline
from __future__ import (division,
print_function)
import os
import sys
import copy
import glob
import fnmatch
import warnings
# Numpy & Scipy
import scipy
import numpy as np
# Astropy related
from astropy.io import fits
from astropy import wcs
from astropy import units as u
from astropy.table import Table, Column, vstack, join
from astropy.stats import sigma_clip
from astropy.cosmology import FlatLambdaCDM as cosmo
cosmo = cosmo(H0=70, Om0=0.3)
# Matplotlib
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
import matplotlib.ticker as ticker
from matplotlib.ticker import NullFormatter
from matplotlib.ticker import MaxNLocator
# AstroML
from astroML.plotting import hist
# Color map
from palettable.colorbrewer.sequential import Oranges_4, Greens_4, Blues_5
ORG4 = Oranges_4.mpl_colormap
BLU5 = Blues_5.mpl_colormap
GRN4 = Greens_4.mpl_colormap
# Matplotlib default settings
rcdef = plt.rcParams.copy()
rcdef['figure.figsize'] = 12, 10
rcdef['xtick.major.size'] = 8.0
rcdef['xtick.major.width'] = 1.5
rcdef['xtick.minor.size'] = 4.0
rcdef['xtick.minor.width'] = 1.5
rcdef['ytick.major.size'] = 8.0
rcdef['ytick.major.width'] = 1.5
rcdef['ytick.minor.size'] = 4.0
rcdef['ytick.minor.width'] = 1.5
rcdef['legend.numpoints'] = 1
#rc('axes', linewidth=2)
# Shapely related imports
from shapely.geometry import Polygon, LineString, Point
from shapely import wkb
from shapely.ops import cascaded_union
from shapely.prepared import prep
from descartes import PolygonPatch
# Personal tools
from hscUtils import songPlotSetup, removeIsNullCol
from hscUtils import confidence_interval, ma_confidence_interval_1d, confidence_interval_1d
# SDSS pivot wavelength
sdss_u_pivot = 3551.0
sdss_g_pivot = 4686.0
sdss_r_pivot = 6165.0
sdss_i_pivot = 7481.0
sdss_z_pivot = 8931.0
# GALEX pivot wavelength
galex_fuv_pivot = 1535.0
galex_nuv_pivot = 2301.0
# WISE pivot wavelength
wise_w1_pivot = 34000.0
wise_w2_pivot = 46000.0
# HSC pivot wavelength
hsc_g_pivot = 4782.2
hsc_r_pivot = 6101.7
hsc_i_pivot = 7648.0
hsc_z_pivot = 8883.0
hsc_y_pivot = 9750.8
hscFiltWave = np.asarray([hsc_g_pivot, hsc_r_pivot, hsc_i_pivot, hsc_z_pivot, hsc_y_pivot])
"""
Absolute magnitude of the Sun in HSC filters
Right now, just use the DES filters
"""
SUN_G = 5.08
SUN_R = 4.62
SUN_I = 4.52
SUN_Z = 4.52
SUN_Y = 4.51
# Solar stellar metallicity
Z_SUN = 0.02
# definitions for the axes
left, width = 0.155, 0.66
right = left + width
bottom, height = 0.13, 0.86
bottom_h = left_h = left + width + 0.02
recScat = [left, bottom, width, height]
recHist = [right, bottom, 0.18, height]
SBP1 = [0.13, 0.12, 0.865, 0.30]
SBP2 = [0.13, 0.42, 0.865, 0.54]
# Color
BLUE0 = "#92c5de"
BLUE1 = "#0571b0"
RED0 = "#f4a582"
RED1 = "#ca0020"
PURPLE0 = '#af8dc3'
PURPLE1 = '#762a83'
BROWN0 = '#bf812d'
BROWN1 = '#543005'
GREEN0 = '#7fbf7b'
GREEN1 = '#1b7837'
In [56]:
def columnExplode(inputCat, column):
"""Explode the array column into column of each bands."""
table = copy.deepcopy(inputCat)
# Separate the columns
col_g = table[column][:,0]
col_r = table[column][:,1]
col_i = table[column][:,2]
col_z = table[column][:,3]
# Rename them
col_g.name = column + '_G'
col_r.name = column + '_R'
col_i.name = column + '_I'
col_z.name = column + '_Z'
# Add back
table.add_column(col_g)
table.add_column(col_r)
table.add_column(col_i)
table.add_column(col_z)
#
if table[column].shape[1] == 5:
col_y = table[column][:,4]
col_y.name = column + '_Y'
table.add_column(col_y)
# Remove old column
table.remove_column(column)
return table
def getObsLogM2L(inputCat, filtUse='I', amag_sun=SUN_I):
"""Use MSTAR, Z, and MAGGIES_FILT to estimate the log(M2L)."""
table = copy.deepcopy(inputCat)
# Get the distance module
distMod = cosmo.distmod(table['Z']).value
# Get the absolute magnitudes
maggieCol = 'MAGGIES_' + filtUse
absmag = (np.log10(table[maggieCol].data) / -0.4) - distMod
# Estimate the luminosity
logLum = (amag_sun - absmag) / 2.5
logLumCol = 'LOGLUM_' + filtUse + '_OBS'
# Estimate the log(M/L)
logM2L = table['MSTAR'].data - logLum
logM2LCol = 'LOGM2L_' + filtUse + '_OBS'
# Add new columns
table.add_column(Column(logLum, name=logLumCol))
table.add_column(Column(logM2L, name=logM2LCol))
return table
def getKcorLogM2L(inputCat, filtUse='I', amag_sun=SUN_I):
"""Use MSTAR, Z, and MAGGIES_FILT to estimate the log(M2L)."""
table = copy.deepcopy(inputCat)
# Get the absolute magnitudes
absmagCol = 'ABSMAG_' + filtUse
absmag = table[absmagCol].data
# Estimate the luminosity
logLum = (amag_sun - absmag) / 2.5
logLumCol = 'LOGLUM_' + filtUse + '_KCOR'
# Estimate the log(M/L)
logM2L = table['MSTAR'].data - logLum
logM2LCol = 'LOGM2L_' + filtUse + '_KCOR'
# Add new columns
table.add_column(Column(logLum, name=logLumCol))
table.add_column(Column(logM2L, name=logM2LCol))
return table
def combineSedResults(sed_run, model_str, sedDir='./', kcorZ='0.0',
outDir=None):
"""Combine the ISEDFIT and KCORRECT results."""
# Name of the result file
file_result = os.path.join(sedDir, sed_run + '_' + model_str + '.fits.gz')
file_kcorr = os.path.join(sedDir, sed_run + '_' + model_str + '_kcorr.z' +
str(kcorZ) + '.fits.gz')
file_post = os.path.join(sedDir, sed_run + '_' + model_str + '_post.fits.gz')
# Read in the SED fitting results
if not os.path.isfile(file_result):
raise Exception("Can not find the ISEDFIT RESULT: %s" % file_result)
else:
sed_result = Table.read(file_result, format='fits')
if outDir is None:
outDir = sedDir
# Read in the Kcorrect results
if not os.path.isfile(file_kcorr):
raise Exception("Can not find the KCORRECT RESULT: %s" % file_kcorr)
else:
sed_kcorr = Table.read(file_kcorr, format='fits')
# New output catalog
file_output = sed_run + '.fits'
# Select the useful columns
useful_result = ['ISEDFIT_ID', 'RA', 'DEC', 'Z', 'MAGGIES', 'IVARMAGGIES', 'BESTMAGGIES', 'CHI2',
'MSTAR', 'AGE', 'SFRAGE', 'TAU', 'ZMETAL', 'AV', 'MU',
'MSTAR_ERR', 'AGE_ERR', 'SFRAGE_ERR', 'TAU_ERR', 'ZMETAL_ERR', 'AV_ERR', 'MU_ERR']
useful_kcorr = ['ISEDFIT_ID', 'CHI2', 'KCORRECT', 'ABSMAG', 'IVARABSMAG', 'SYNTH_ABSMAG']
# Isolate the useful columns into new tables
new_result = sed_result[useful_result]
new_kcorr = sed_kcorr[useful_kcorr]
# Rename the chi^2 columns
new_result.rename_column('CHI2', 'CHI2_SED')
new_kcorr.rename_column('CHI2', 'CHI2_KCORR')
# Join the two tables
new_combine = join(new_result, new_kcorr, keys='ISEDFIT_ID')
# Get the relative metallicity
zmetalRel = (new_combine['ZMETAL'] / Z_SUN)
new_combine.add_column(Column(zmetalRel, name='ZMET_REL'))
# Explode the columns in SED catalogs
new_combine = columnExplode(new_combine, 'MAGGIES')
new_combine = columnExplode(new_combine, 'IVARMAGGIES')
new_combine = columnExplode(new_combine, 'BESTMAGGIES')
# Magnitude differences
new_combine.add_column(Column(
(2.5 * np.log10(new_combine['BESTMAGGIES_G'] / new_combine['MAGGIES_G'])),
name='DMAG_G'))
new_combine.add_column(Column(
(2.5 * np.log10(new_combine['BESTMAGGIES_R'] / new_combine['MAGGIES_R'])),
name='DMAG_R'))
new_combine.add_column(Column(
(2.5 * np.log10(new_combine['BESTMAGGIES_I'] / new_combine['MAGGIES_I'])),
name='DMAG_I'))
new_combine.add_column(Column(
(2.5 * np.log10(new_combine['BESTMAGGIES_Z'] / new_combine['MAGGIES_Z'])),
name='DMAG_Z'))
try:
new_combine.add_column(Column(
(2.5 * np.log10(new_combine['BESTMAGGIES_Y'] / new_combine['MAGGIES_Y'])),
name='DMAG_Y'))
except Exception:
pass
# Explode the columns in Kcorrect catalogs
new_combine = columnExplode(new_combine, 'KCORRECT')
new_combine = columnExplode(new_combine, 'ABSMAG')
new_combine = columnExplode(new_combine, 'IVARABSMAG')
new_combine = columnExplode(new_combine, 'SYNTH_ABSMAG')
# logM2L in observed frame
new_combine = getObsLogM2L(new_combine, filtUse='I', amag_sun=SUN_I)
new_combine = getObsLogM2L(new_combine, filtUse='G', amag_sun=SUN_G)
new_combine = getObsLogM2L(new_combine, filtUse='R', amag_sun=SUN_R)
new_combine = getObsLogM2L(new_combine, filtUse='Z', amag_sun=SUN_Z)
try:
new_combine = getObsLogM2L(new_combine, filtUse='Y', amag_sun=SUN_Y)
except Exception:
pass
# logM2L in rest frame
new_combine = getKcorLogM2L(new_combine, filtUse='I', amag_sun=SUN_I)
new_combine = getKcorLogM2L(new_combine, filtUse='G', amag_sun=SUN_G)
new_combine = getKcorLogM2L(new_combine, filtUse='R', amag_sun=SUN_R)
new_combine = getKcorLogM2L(new_combine, filtUse='Z', amag_sun=SUN_Z)
try:
new_combine = getKcorLogM2L(new_combine, filtUse='Y', amag_sun=SUN_Y)
except Exception:
pass
# Save the results to new table
new_combine.write(os.path.join(outDir, file_output), format='fits',
overwrite=True)
return new_combine
def toColorArr(data, bottom=None, top=None):
"""
Convert a data array to "color array" (between 0 and 1).
Parameters:
bottom, top :
"""
if top is not None:
data[data >= top] = top
if bottom is not None:
data[data <= bottom] = bottom
return ((data - np.nanmin(data)) /
(np.nanmax(data) - np.nanmin(data))) * 255.0
def toSizeArr(data, bottom=None, top=None, maxSize=40):
"""
Convert a data array to "size array".
Parameters:
bottom, top :
"""
if top is not None:
data[data >= top] = top
if bottom is not None:
data[data <= bottom] = bottom
return ((data - np.nanmin(data)) /
(np.nanmax(data) - np.nanmin(data))) * maxSize
def isedfitPlot(tableUse, colName1, colName2, titleUse=None,
colColor='CHI2_SED', colSize='Z',
cmap=ORG4, maxSize=350, alpha=0.8,
ax=None, xylim=True, nolabel=False,
titleSize=20, labelSize=36,
titleX=0.04, titleY=0.92,
xmargin=7.5, ymargin=7.5, noEdge=True):
#---------------------------------------------------------#
if ax is None:
fig = plt.figure(figsize=(10, 10))
fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
ax1 = fig.add_subplot(111)
else:
ax1 = ax
# Formatting
ax1 = songPlotSetup(ax1)
# Scatter plot
if noEdge:
ax1.scatter(tableUse[colName1],
tableUse[colName2],
c=toColorArr(tableUse[colColor]),
s=toSizeArr(tableUse[colSize], maxSize=maxSize),
alpha=alpha, cmap=cmap, edgecolor='none')
else:
ax1.scatter(tableUse[colName1],
tableUse[colName2],
c=toColorArr(tableUse[colColor]),
s=toSizeArr(tableUse[colSize], maxSize=maxSize),
alpha=alpha, cmap=cmap)
# Label
if not nolabel:
colnameUse1 = colName1.replace('_', '\_')
colnameUse2 = colName2.replace('_', '\_')
ax1.set_xlabel('$\mathrm{%s}$' % colnameUse1, size=labelSize)
ax1.set_ylabel('$\mathrm{%s}$' % colnameUse2, size=labelSize)
# Axis limits
if xylim:
xmin, xmax = np.nanmin(tableUse[colName1]), np.nanmax(tableUse[colName1])
ymin, ymax = np.nanmin(tableUse[colName2]), np.nanmax(tableUse[colName2])
xmargin, ymargin = ((xmax - xmin) / xmargin), ((ymax - ymin) / ymargin)
ax1.set_xlim(xmin-xmargin, xmax+xmargin)
ax1.set_ylim(ymin-ymargin, ymax+ymargin)
# Title
if titleUse is not None:
titleUse = titleUse.replace('_', '\_')
ax1.text(titleX, titleY, '$\mathrm{%s}$' % titleUse, size=titleSize,
transform = ax1.transAxes)
#---------------------------------------------------------#
if ax is None:
return fig
else:
return ax1
#---------------------------------------------------------#
def simpleScatter(xarr, yarr,
xtick=True, ytick=True,
xstr=None, ystr=None, titleUse=None,
carr=None, sarr=None,
cmap=ORG4, maxSize=350, alpha=0.8,
ax=None, xylim=True, nolabel=False,
titleSize=20, labelSize=36,
titleX=0.04, titleY=0.92,
xmargin=7.5, ymargin=7.5,
noEdge=True, colorUse='k', sizeUse=40):
#---------------------------------------------------------#
if ax is None:
fig = plt.figure(figsize=(10, 10))
fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
ax1 = fig.add_subplot(111)
# Formatting
ax1 = songPlotSetup(ax1)
else:
ax1 = ax
# Color Array
if carr is None:
carr = colorUse
else:
carr = toColorArr(carr)
# Size Array
if sarr is None:
sarr = sizeUse
else:
sarr = toSizeArr(sarr, maxSize=maxSize)
# Scatter plot
if noEdge:
ax1.scatter(xarr, yarr, c=carr, s=sarr,
alpha=alpha, cmap=cmap, edgecolor='none',
rasterized=True)
else:
ax1.scatter(xarr, yarr, c=carr, s=sarr,
alpha=alpha, cmap=cmap,
rasterized=True)
# Label
if xstr is not None:
ax1.set_xlabel(xstr, size=labelSize)
if ystr is not None:
ax1.set_ylabel(ystr, size=labelSize)
# Axis limits
if xylim:
xmin, xmax = np.nanmin(xarr), np.nanmax(xarr)
ymin, ymax = np.nanmin(yarr), np.nanmax(yarr)
xmargin, ymargin = ((xmax - xmin) / xmargin), ((ymax - ymin) / ymargin)
ax1.set_xlim(xmin-xmargin, xmax+xmargin)
ax1.set_ylim(ymin-ymargin, ymax+ymargin)
# Title
if titleUse is not None:
titleUse = titleUse.replace('_', '\_')
ax1.text(titleX, titleY, '$\mathrm{%s}$' % titleUse, size=titleSize,
transform = ax1.transAxes)
#---------------------------------------------------------#
if ax is None:
return fig
else:
return ax1
def simpleHist(tableUse, colName, sample1=None, bins='knuth',
alpha=0.5, table2=None, sample2=None,
showLegend=True, useLog=False):
#---------------------------------------------------------#
fig = plt.figure(figsize=(10, 10))
fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
ax1 = fig.add_subplot(111)
# Formatting
ax1 = songPlotSetup(ax1)
# RedMapper
if sample1 is None:
sample1 = 'Sample\ 1'
if useLog:
hist(np.log10(tableUse[colName]), bins=bins, histtype='stepfilled',
color='b', alpha=alpha, normed=True,
label=sample1, ax=ax1)
else:
hist(tableUse[colName], bins=bins, histtype='stepfilled',
color='b', alpha=alpha, normed=True,
label=sample1, ax=ax1)
# Label
colnameUse = colName.replace('_', '\_')
if useLog:
ax1.set_xlabel('$\log\ \mathrm{%s}$' % colnameUse, size=36)
else:
ax1.set_xlabel('$\mathrm{%s}$' % colnameUse, size=36)
ax1.set_ylabel('$\mathrm{\#}$', size=36)
# Secondary sample ?
if table2 is not None:
if sample2 is None:
sample2 = 'Sample\ 2'
if useLog:
hist(np.log10(table2[colName]), bins=bins, histtype='step',
color='orange', normed=True, label=sample2, ax=ax1,
linewidth=4.5)
else:
hist(table2[colName], bins=bins, histtype='step',
color='orange', normed=True, label=sample2, ax=ax1,
linewidth=4.5)
# Legend
if showLegend:
l_handles, l_labels = ax1.get_legend_handles_labels()
ax1.legend(l_handles, l_labels, loc=(0.05, 0.84),
shadow=True, fancybox=True,
numpoints=1, fontsize=26, scatterpoints=1,
markerscale=1.8, borderpad=0.2, handletextpad=0.1)
#---------------------------------------------------------#
return fig
def isedMassPlot(sedRes, sedStr,
xSize=20, ySize=16,
massThresh=10.2, labelSize=26,
cmap=ORG4, alpha=0.75, maxSize=300,
tickMajor=0.4, tickMinor=0.1,
dpi=100, xmargin=7.5,
ymargin1=7.5, ymargin2=7.5,
ymargin3=6.0, ymargin4=6.0,
titleX=0.10, titleY=0.91, titleSize=35,
xLabelSize=40, yLabelSize=40,
saveFig=True, outDir='./'):
"""Summary plots of iSEDFit results (Stellar mass related)"""
#-----------------------------------------------------------------#
# Exclude problematic results
sedRes = sedRes[sedRes['MSTAR'] >= massThresh]
#-----------------------------------------------------------------#
fig = plt.figure(figsize=(xSize, ySize))
fig.subplots_adjust(left=0.10, right=0.90,
bottom=0.08, top=0.99,
wspace=0.02, hspace=0.00)
# logM - Age plot
ax1 = fig.add_subplot(221)
fig1 = isedfitPlot(sedRes, 'MSTAR', 'AGE',
titleUse=sedStr, titleSize=titleSize,
colColor='Z', colSize='DMAG_I',
ax=ax1, nolabel=True, alpha=alpha,
titleX=titleX, titleY=titleY,
cmap=cmap, maxSize=maxSize,
labelSize=labelSize,
xmargin=xmargin, ymargin=ymargin1)
##
ax1.tick_params(axis='x', which='both', labelbottom='off')
ax1.set_ylabel('$\mathrm{Age\ (Gyr)}$', size=yLabelSize)
ax1.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
ax1.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
ax1.text(0.74, 0.14, '$\mathrm{Color:\ }z$', size=26,
transform = ax1.transAxes)
ax1.text(0.74, 0.05, '$\mathrm{Size:\ }{\Delta\ i}$', size=26,
transform = ax1.transAxes)
# logM - SFR Age plot
ax2 = fig.add_subplot(223)
fig2 = isedfitPlot(sedRes, 'MSTAR', 'SFRAGE',
titleUse=None, cmap=cmap, alpha=alpha,
colColor='Z', colSize='DMAG_I',
ax=ax2, nolabel=True, maxSize=maxSize,
labelSize=labelSize,
xmargin=xmargin, ymargin=ymargin2)
##
ax2.set_xlabel('$\log\ (M_{\star}/M_{\odot})$', size=xLabelSize)
ax2.set_ylabel('$\mathrm{Age_{\ SFR}\ (Gyr)}$', size=yLabelSize)
ax2.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
ax2.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
# logM - ZMetal
ax3 = fig.add_subplot(222)
fig3 = isedfitPlot(sedRes, 'MSTAR', 'ZMET_REL',
titleUse=None, cmap=cmap, alpha=alpha,
colColor='Z', colSize='DMAG_I',
ax=ax3, nolabel=True, maxSize=maxSize,
labelSize=labelSize,
xmargin=xmargin, ymargin=ymargin3)
##
ax3.yaxis.set_label_position("right")
ax3.tick_params(axis='y', which='both', labelleft='off',
labelright='on', labelsize=labelSize)
ax3.tick_params(axis='x', which='both', labelbottom='off')
ax3.set_ylabel('$\mathrm{Z}_{\star}/\mathrm{Z}_{\odot}$', size=yLabelSize)
ax3.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
ax3.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
ax3.yaxis.set_major_locator(ticker.MultipleLocator(0.4))
ax3.yaxis.set_minor_locator(ticker.MultipleLocator(0.1))
# logM - Av
ax4 = fig.add_subplot(224)
fig4 = isedfitPlot(sedRes, 'MSTAR', 'AV',
titleUse=None, cmap=cmap, alpha=alpha,
colColor='Z', colSize='DMAG_I',
ax=ax4, nolabel=True, maxSize=maxSize,
labelSize=labelSize,
xmargin=xmargin, ymargin=ymargin4)
##
ax4.tick_params(axis='y', which='both', labelsize=labelSize,
labelleft='off', labelright='on')
ax4.yaxis.set_label_position("right")
ax4.set_xlabel('$\log\ (M_{\star}/M_{\odot})$', size=xLabelSize)
ax4.set_ylabel('$\mathrm{A}_{V}\ \mathrm{(mag)}$', size=yLabelSize)
ax4.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
ax4.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
ax4.yaxis.set_major_locator(ticker.MultipleLocator(0.04))
ax4.yaxis.set_minor_locator(ticker.MultipleLocator(0.01))
#-----------------------------------------------------------------#
massPng = sedStr + '_logm_plots.pdf'
if saveFig:
fig.savefig(os.path.join(outDir, massPng), dpi=dpi)
plt.close(fig)
#-----------------------------------------------------------------#
def isedChi2Plot(sedRes, sedStr,
xSize=20, ySize=16,
massThresh=10.2, labelSize=26,
cmap=ORG4, alpha=0.75, maxSize=300,
tickMajor=None, tickMinor=None,
dpi=100, xmargin=7.5,
ymargin1=7.5, ymargin2=7.5,
ymargin3=6.0, ymargin4=6.0,
titleX=0.10, titleY=0.91, titleSize=35,
xLabelSize=40, yLabelSize=40,
saveFig=True, outDir='./'):
"""Summary plots of iSEDFit results (Chi^2 related)"""
#-----------------------------------------------------------------#
# Exclude problematic results
sedRes = sedRes[sedRes['MSTAR'] >= massThresh]
#-----------------------------------------------------------------#
fig = plt.figure(figsize=(xSize, ySize))
fig.subplots_adjust(left=0.10, right=0.90,
bottom=0.08, top=0.99,
wspace=0.02, hspace=0.00)
# logM - Age plot
ax1 = fig.add_subplot(221)
fig1 = isedfitPlot(sedRes, 'CHI2_SED', 'MSTAR',
titleUse=sedStr, titleSize=titleSize,
colColor='MSTAR', colSize='Z',
ax=ax1, nolabel=True, alpha=alpha,
titleX=titleX, titleY=titleY,
cmap=cmap, maxSize=maxSize,
labelSize=labelSize,
xmargin=xmargin, ymargin=ymargin1)
##
ax1.tick_params(axis='x', which='both', labelbottom='off')
ax1.set_ylabel('$\log\ (M_{\star}/M_{\odot})$', size=yLabelSize)
if tickMajor is not None:
ax1.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
if tickMinor is not None:
ax1.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
ax1.text(0.75, 0.14, '$\mathrm{Size:\ }z$', size=26,
transform = ax1.transAxes)
ax1.text(0.72, 0.06, '$\mathrm{Color:\ Mass}$', size=26,
transform = ax1.transAxes)
# logM - SFR Age plot
ax2 = fig.add_subplot(223)
fig2 = isedfitPlot(sedRes, 'CHI2_SED', 'Z',
titleUse=None, cmap=cmap, alpha=alpha,
colColor='MSTAR', colSize='Z',
ax=ax2, nolabel=True, maxSize=maxSize,
labelSize=labelSize,
xmargin=xmargin, ymargin=ymargin2)
##
ax2.set_xlabel('${\chi}^2$', size=xLabelSize)
ax2.set_ylabel('$\mathrm{Redshift}$', size=yLabelSize)
if tickMajor is not None:
ax2.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
if tickMinor is not None:
ax2.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
# logM - ZMetal
ax3 = fig.add_subplot(222)
fig3 = isedfitPlot(sedRes, 'CHI2_SED', 'AGE',
titleUse=None, cmap=cmap, alpha=alpha,
colColor='MSTAR', colSize='Z',
ax=ax3, nolabel=True, maxSize=maxSize,
labelSize=labelSize,
xmargin=xmargin, ymargin=ymargin3)
##
ax3.yaxis.set_label_position("right")
ax3.tick_params(axis='y', which='both', labelleft='off',
labelright='on', labelsize=labelSize)
ax3.tick_params(axis='x', which='both', labelbottom='off')
ax3.set_ylabel('$\mathrm{Age\ (Gyr)}$', size=yLabelSize)
if tickMajor is not None:
ax3.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
if tickMinor is not None:
ax3.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
# logM - Av
ax4 = fig.add_subplot(224)
fig4 = isedfitPlot(sedRes, 'CHI2_SED', 'AV',
titleUse=None, cmap=cmap, alpha=alpha,
colColor='MSTAR', colSize='Z',
ax=ax4, nolabel=True, maxSize=maxSize,
labelSize=labelSize,
xmargin=xmargin, ymargin=ymargin4)
##
ax4.tick_params(axis='y', which='both', labelsize=labelSize,
labelleft='off', labelright='on')
ax4.yaxis.set_label_position("right")
ax4.set_xlabel('${\chi}^2$', size=xLabelSize)
ax4.set_ylabel('$\mathrm{A}_{V}\ \mathrm{(mag)}$', size=yLabelSize)
if tickMajor is not None:
ax4.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
if tickMinor is not None:
ax4.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
ax4.yaxis.set_major_locator(ticker.MultipleLocator(0.04))
ax4.yaxis.set_minor_locator(ticker.MultipleLocator(0.01))
#-----------------------------------------------------------------#
chi2Png = sedStr + '_chi2_plots.pdf'
if saveFig:
fig.savefig(os.path.join(outDir, chi2Png), dpi=dpi)
plt.close(fig)
#-----------------------------------------------------------------#
def isedRedshiftPlot(sedRes, sedStr,
xSize=20, ySize=16,
massThresh=10.2, labelSize=26,
cmap=ORG4, alpha=0.75, maxSize=300,
tickMajor=None, tickMinor=None,
dpi=100, xmargin=7.5,
ymargin1=7.5, ymargin2=7.5,
ymargin3=6.0, ymargin4=6.0,
titleX=0.10, titleY=0.91, titleSize=35,
xLabelSize=40, yLabelSize=40,
saveFig=True, outDir='./'):
"""Summary plots of iSEDFit results (Redshift related)"""
#-----------------------------------------------------------------#
# Exclude problematic results
sedRes = sedRes[sedRes['MSTAR'] >= massThresh]
#-----------------------------------------------------------------#
fig = plt.figure(figsize=(xSize, ySize))
fig.subplots_adjust(left=0.10, right=0.90,
bottom=0.08, top=0.99,
wspace=0.02, hspace=0.00)
# logM - Age plot
ax1 = fig.add_subplot(221)
fig1 = isedfitPlot(sedRes, 'Z', 'MSTAR',
titleUse=sedStr, titleSize=titleSize,
colColor='MSTAR', colSize='DMAG_I',
ax=ax1, nolabel=True, alpha=alpha,
titleX=titleX, titleY=titleY,
cmap=cmap, maxSize=maxSize,
labelSize=labelSize,
xmargin=xmargin, ymargin=ymargin1)
##
ax1.tick_params(axis='x', which='both', labelbottom='off')
ax1.set_ylabel('$\log\ (M_{\star}/M_{\odot})$', size=yLabelSize)
if tickMajor is not None:
ax1.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
if tickMinor is not None:
ax1.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
ax1.text(0.72, 0.14, '$\mathrm{Color:\ Mass}$', size=26,
transform = ax1.transAxes)
ax1.text(0.72, 0.06, '$\mathrm{Size:\ }{\Delta\ i}$', size=26,
transform = ax1.transAxes)
# logM - SFR Age plot
ax2 = fig.add_subplot(223)
fig2 = isedfitPlot(sedRes, 'Z', 'AGE',
titleUse=None, cmap=cmap, alpha=alpha,
colColor='MSTAR', colSize='DMAG_I',
ax=ax2, nolabel=True, maxSize=maxSize,
labelSize=labelSize,
xmargin=xmargin, ymargin=ymargin2)
##
ax2.set_xlabel('$\mathrm{Redshift}$', size=xLabelSize)
ax2.set_ylabel('$\mathrm{Age\ (Gyr)}$', size=yLabelSize)
if tickMajor is not None:
ax2.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
if tickMinor is not None:
ax2.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
# logM - ZMetal
ax3 = fig.add_subplot(222)
fig3 = isedfitPlot(sedRes, 'Z', 'ZMET_REL',
titleUse=None, cmap=cmap, alpha=alpha,
colColor='MSTAR', colSize='DMAG_I',
ax=ax3, nolabel=True, maxSize=maxSize,
labelSize=labelSize,
xmargin=xmargin, ymargin=ymargin3)
##
ax3.yaxis.set_label_position("right")
ax3.tick_params(axis='y', which='both', labelleft='off',
labelright='on', labelsize=labelSize)
ax3.tick_params(axis='x', which='both', labelbottom='off')
ax3.set_ylabel('$\mathrm{Z}_{\star}/\mathrm{Z}_{\odot}$', size=yLabelSize)
if tickMajor is not None:
ax3.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
if tickMinor is not None:
ax3.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
# logM - Av
ax4 = fig.add_subplot(224)
fig4 = isedfitPlot(sedRes, 'Z', 'AV',
titleUse=None, cmap=cmap, alpha=alpha,
colColor='MSTAR', colSize='DMAG_I',
ax=ax4, nolabel=True, maxSize=maxSize,
labelSize=labelSize,
xmargin=xmargin, ymargin=ymargin4)
##
ax4.tick_params(axis='y', which='both', labelsize=labelSize,
labelleft='off', labelright='on')
ax4.yaxis.set_label_position("right")
ax4.set_xlabel('$\mathrm{Redshift}$', size=xLabelSize)
ax4.set_ylabel('$\mathrm{A}_{V}\ \mathrm{(mag)}$', size=yLabelSize)
if tickMajor is not None:
ax4.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
if tickMinor is not None:
ax4.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
ax4.yaxis.set_major_locator(ticker.MultipleLocator(0.04))
ax4.yaxis.set_minor_locator(ticker.MultipleLocator(0.01))
#-----------------------------------------------------------------#
zredPng = sedStr + '_zred_plots.pdf'
if saveFig:
fig.savefig(os.path.join(outDir, zredPng), dpi=dpi)
plt.close(fig)
#-----------------------------------------------------------------#
def isedMagDiffPlot(sedRes, filterUse,
xSize=11, ySize=32, sedStr=None,
massThresh=10.2, labelSize=26,
cmap=ORG4, alpha=0.75, maxSize=300,
tickMajor=None, tickMinor=None,
dpi=100, xmargin=7.5,
ymargin1=6.0, ymargin2=7.5,
ymargin3=6.0, ymargin4=6.0,
titleX=0.10, titleY=0.91, titleSize=35,
xLabelSize=40, yLabelSize=40,
saveFig=True, outDir='./'):
"""Summary plots of iSEDFit results (Redshift related)"""
#-----------------------------------------------------------------#
# Exclude problematic results
sedRes = sedRes[sedRes['MSTAR'] >= massThresh]
#-----------------------------------------------------------------#
fig = plt.figure(figsize=(xSize, ySize))
fig.subplots_adjust(left=0.17, right=0.99,
top=0.99, bottom=0.04,
wspace=0.02, hspace=0.00)
# X-array
xarr = 2.5 * np.log10(sedRes['BESTMAGGIES_' + filterUse] /
sedRes['MAGGIES_' + filterUse])
xstr = '$\mathrm{Observed}-\mathrm{BestFit}\ \mathrm{(%s\ band)}$' % filterUse
carr = sedRes['Z']
sarr = sedRes['MSTAR']
# 1. diff_Mag v.s. Chi^2 plot
ax1 = fig.add_subplot(411)
yarr1 = sedRes['CHI2_SED']
fig_a = simpleScatter(xarr, yarr1,
xstr=None, ystr=None,
titleUse=sedStr, titleSize=titleSize,
carr=carr, sarr=sarr, ax=ax1,
xmargin=xmargin, ymargin=ymargin1,
cmap=cmap, maxSize=maxSize, alpha=alpha)
##
ax1.tick_params(axis='x', which='both', labelbottom='off')
ax1.set_ylabel('${\chi}^2$', size=yLabelSize)
if tickMajor is not None:
ax1.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
if tickMinor is not None:
ax1.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
ax1.text(0.08, 0.86, '$\mathrm{Size:\ Mass;\ \mathrm{Color:\ }z}$',
size=26, transform = ax1.transAxes)
# 2. diff_Mag v.s. Mstar
ax2 = fig.add_subplot(412)
yarr2 = sedRes['MSTAR']
fig_b = simpleScatter(xarr, yarr2,
xstr=None, ystr=None,
titleUse=None,
carr=carr, sarr=sarr, ax=ax2,
xmargin=xmargin, ymargin=ymargin2,
cmap=cmap, maxSize=maxSize, alpha=alpha)
##
ax2.tick_params(axis='x', which='both', labelbottom='off')
ax2.set_ylabel('$\log\ (M_{\star}/M_{\odot})$', size=yLabelSize)
if tickMajor is not None:
ax2.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
if tickMinor is not None:
ax2.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
# 3. diff_Mag v.s. Age
ax3 = fig.add_subplot(413)
yarr3 = sedRes['AGE']
fig_c = simpleScatter(xarr, yarr3,
xstr=None, ystr=None,
titleUse=None,
carr=carr, sarr=sarr, ax=ax3,
xmargin=xmargin, ymargin=ymargin3,
cmap=cmap, maxSize=maxSize, alpha=alpha)
##
ax3.tick_params(axis='x', which='both', labelbottom='off')
ax3.set_ylabel('$\mathrm{Age\ (Gyr)}$', size=yLabelSize)
if tickMajor is not None:
ax3.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
if tickMinor is not None:
ax3.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
# 4. diff_Mag v.s. A_V
ax4 = fig.add_subplot(414)
yarr4 = sedRes['AV']
fig_d = simpleScatter(xarr, yarr4,
xstr=None, ystr=None,
titleUse=None,
carr=carr, sarr=sarr, ax=ax4,
xmargin=xmargin, ymargin=ymargin4,
cmap=cmap, maxSize=maxSize, alpha=alpha)
##
ax4.set_ylabel('$\mathrm{A}_{V}\ \mathrm{(mag)}$', size=yLabelSize)
ax4.set_xlabel(xstr, size=xLabelSize)
if tickMajor is not None:
ax4.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
if tickMinor is not None:
ax4.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
#-----------------------------------------------------------------#
dmagPng = sedStr + '_dmag_' + filterUse.strip() + '_plots.pdf'
if saveFig:
fig.savefig(os.path.join(outDir, dmagPng), dpi=dpi)
plt.close(fig)
#-----------------------------------------------------------------#
def sedSummaryPlot(sedRes, sedStr, massThresh=10.2, outDir='./',
cmapUse=ORG4, alphaUse=0.7, maxSizeUse=200):
## Mass Plots
isedMassPlot(sedRes, sedStr, cmap=cmapUse, alpha=alphaUse, maxSize=maxSizeUse,
tickMajor=0.4, tickMinor=0.1, xmargin=7.5,
ymargin1=7.5, ymargin2=7.5, ymargin3=6.0, ymargin4=6.0,
saveFig=True, outDir=outDir)
## Chi2 Plots
isedChi2Plot(sedRes, sedStr, cmap=cmapUse, alpha=alphaUse, maxSize=maxSizeUse,
tickMajor=None, tickMinor=None, xmargin=7.5,
ymargin1=7.5, ymargin2=7.5, ymargin3=6.0, ymargin4=6.0,
saveFig=True, outDir=outDir)
## Redshift Plots
isedRedshiftPlot(sedRes, sedStr, cmap=cmapUse, alpha=alphaUse, maxSize=maxSizeUse,
tickMajor=None, tickMinor=None, xmargin=7.5,
ymargin1=7.5, ymargin2=7.5, ymargin3=6.0, ymargin4=6.0,
saveFig=True, outDir=outDir)
## Magnitude differences
### Magnitude difference in Y-band
try:
isedMagDiffPlot(sedRes, 'Y',
xSize=11, ySize=32, sedStr=sedStr,
cmap=cmapUse, alpha=alphaUse, maxSize=maxSizeUse,
tickMajor=None, tickMinor=None,
saveFig=True, outDir=outDir)
except Exception:
pass
### Magnitude difference in I-band
isedMagDiffPlot(sedRes, 'I',
xSize=11, ySize=32, sedStr=sedStr,
cmap=cmapUse, alpha=alphaUse, maxSize=maxSizeUse,
tickMajor=None, tickMinor=None,
saveFig=True, outDir=outDir)
### Magnitude difference in G-band
isedMagDiffPlot(sedRes, 'G',
xSize=11, ySize=32, sedStr=sedStr,
cmap=cmapUse, alpha=alphaUse, maxSize=maxSizeUse,
tickMajor=None, tickMinor=None,
saveFig=True, outDir=outDir)
### Magnitude difference in R-band
isedMagDiffPlot(sedRes, 'R',
xSize=11, ySize=32, sedStr=sedStr,
cmap=cmapUse, alpha=alphaUse, maxSize=maxSizeUse,
tickMajor=None, tickMinor=None,
saveFig=True, outDir=outDir)
### Magnitude difference in Z-band
isedMagDiffPlot(sedRes, 'Z',
xSize=11, ySize=32, sedStr=sedStr,
cmap=cmapUse, alpha=alphaUse, maxSize=maxSizeUse,
tickMajor=None, tickMinor=None,
saveFig=True, outDir=outDir)
def sedMassCompare(sed1, sed2, sample, label1, label2,
location='./', cmapUse=ORG4, colorUse=RED0,
alphaUse=0.75, maxSizeUse=300,
massCol='MSTAR', colorCol='Z', sizeCol='DMAG_I',
errCol='MSTAR_ERR', chi2Col='CHI2_SED',
chi2Thr=10.0, outDir='./',
xtickFormat='$\mathbf{%4.1f}$',
ytickFormat='$\mathbf{%4.2f}$',
xmin=None, ymin=None, xmax=None, ymax=None):
"""Compare the mass estimated from two iSEDFit run"""
# ------------------------------------------------------------------------------#
strRef = '$\mathrm{%s}$' % label1.upper()
strComp = '$\mathrm{%s}$' % label2.upper()
# ------------------------------------------------------------------------------#
cat1 = Table.read(os.path.join(location, sed1 + '.fits'), format='fits')
cat2 = Table.read(os.path.join(location, sed2 + '.fits'), format='fits')
# Exclude problematic results
catRef = cat1[(cat1[massCol] >= 10.0) & (cat2[massCol] >= 10.0)]
catComp = cat2[(cat1[massCol] >= 10.0) & (cat2[massCol] >= 10.0)]
# ------------------------------------------------------------------------------#
xarr = catRef[massCol]
yarr = catComp[massCol] - catRef[massCol]
carr = catRef[colorCol]
sarr = catRef[sizeCol]
# ------------------------------------------------------------------------------#
xerr = catRef[errCol]
sigErr = np.nanstd(xerr)
medErr = np.nanmedian(xerr)
# ------------------------------------------------------------------------------#
chi2Ref = catRef[chi2Col]
fracBadRef = (len(chi2Ref[chi2Ref >= chi2Thr]) * 100.0) / len(chi2Ref)
chi2Comp = catComp[chi2Col]
fracBadComp = (len(chi2Comp[chi2Comp >= chi2Thr]) * 100.0) / len(chi2Comp)
# ------------------------------------------------------------------------------#
xmargin, ymargin = 10.0, 5.5
# ------------------------------------------------------------------------------#
# Mass v.s. Mass Difference
fig = plt.figure(figsize=(13, 9))
ax1 = plt.axes(recScat)
ax2 = plt.axes(recHist)
ax1 = songPlotSetup(ax1, xtickFormat=xtickFormat, ytickFormat=ytickFormat)
ax2 = songPlotSetup(ax2, xtickFormat=xtickFormat, ytickFormat=ytickFormat)
# ------------------------------------------------------------------------------#
# Scatter plot
# Shaded region to show the typical error
minM, maxM = np.nanmin(xarr), np.nanmax(xarr)
sepM = (maxM - minM) / (xmargin - 1.5)
tempArr = np.asarray([minM - sepM, maxM + sepM])
ax1.fill_between(tempArr,
(tempArr * 0.0 - 3.0 * sigErr),
(tempArr * 0.0 + 3.0 * sigErr), alpha=0.05,
facecolor='k', edgecolor='none', interpolate=True)
ax1.fill_between(tempArr,
(tempArr * 0.0 - 1.0 * medErr),
(tempArr * 0.0 + 1.0 * medErr), alpha=0.05,
facecolor='k', edgecolor='none', interpolate=True)
## Horionzital 0.0 Line
ax1.axhline(0.0, linewidth=4.0, linestyle='--', c='k', alpha=0.4)
# Matched ones
fig_a = simpleScatter(xarr, yarr, xstr=None, ystr=None,
titleUse=None, carr=carr, sarr=sarr, ax=ax1,
xmargin=xmargin, ymargin=ymargin,
cmap=cmapUse, maxSize=maxSizeUse, alpha=alphaUse)
# X, Y limits
if (xmin is not None) and (xmax is not None):
ax1.set_xlim(xmin, xmax)
if (ymin is not None) and (ymax is not None):
ax1.set_ylim(ymin, ymax)
# Text information about the comparison
ax1.text(0.05, 0.92,
strRef + '$: <{\chi}^2>=%4.1f; %4.1f$%%' % (np.nanmedian(chi2Ref), fracBadRef),
verticalalignment='bottom', horizontalalignment='left',
fontsize=26.0, transform=ax1.transAxes, color='k')
ax1.text(0.05, 0.86,
strComp + '$: <{\chi}^2>=%4.1f; %4.1f%%$%%' % (np.nanmedian(chi2Comp), fracBadComp),
verticalalignment='bottom', horizontalalignment='left',
fontsize=26.0, transform=ax1.transAxes, color='k')
# Label
ax1.set_xlabel('$\log\ (M_{\star}/M_{\odot})\ ($' + strRef + '$)$', size=40)
ax1.set_ylabel('$\Delta(\log M{\star})\ $' + '(' + strComp + '-' + strRef + ')',
size=40)
# ---------------------------------------------------------------------------
# Histogram
ax2.set_ylim(ax1.get_ylim())
# Horionzital 0.0 Line
ax2.axhline(0.0, linewidth=4.0, linestyle='--', c='k', alpha=0.4)
n, bins, patches=ax2.hist(yarr, bins=30, edgecolor='none',
orientation='horizontal', histtype='stepfilled',
color=colorUse, alpha=0.80, normed=1)
# Axes setup
ax2.tick_params(axis='x', which='minor', bottom='off', top='off')
ax2.xaxis.set_major_formatter(NullFormatter())
ax2.yaxis.set_major_formatter(NullFormatter())
# Save the figure
fig.savefig(os.path.join(outDir, '%s_sed_logmdiff_%s_%s.png' % (sample, label1, label2)),
dpi=200)
In [5]:
results_col = ['ISEDFIT_ID', 'RA', 'DEC', 'Z', 'MAGGIES', 'IVARMAGGIES', 'BESTMAGGIES', 'CHUNKINDX', 'MODELINDX', 'DELAYED', 'BURSTTYPE', 'CHI2', 'TOTALMASS', 'TOTALMASS_ERR', 'MSTAR', 'AGE', 'SFRAGE', 'TAU', 'ZMETAL', 'AV', 'MU', 'OIIIHB', 'NLYC', 'SFR', 'SFR100', 'B100', 'B1000', 'EWOII', 'EWOIIIHB', 'EWNIIHA', 'NBURST', 'TRUNCTAU', 'TBURST', 'DTBURST', 'FBURST', 'MSTAR_50', 'AGE_50', 'SFRAGE_50', 'TAU_50', 'ZMETAL_50', 'AV_50', 'MU_50', 'OIIIHB_50', 'SFR_50', 'SFR100_50', 'B100_50', 'B1000_50', 'EWOII_50', 'EWOIIIHB_50', 'EWNIIHA_50', 'MSTAR_AVG', 'AGE_AVG', 'SFRAGE_AVG', 'TAU_AVG', 'ZMETAL_AVG', 'AV_AVG', 'MU_AVG', 'OIIIHB_AVG', 'SFR_AVG', 'SFR100_AVG', 'B100_AVG', 'B1000_AVG', 'EWOII_AVG', 'EWOIIIHB_AVG', 'EWNIIHA_AVG', 'MSTAR_ERR', 'AGE_ERR', 'SFRAGE_ERR', 'TAU_ERR', 'ZMETAL_ERR', 'AV_ERR', 'MU_ERR', 'OIIIHB_ERR', 'SFR_ERR', 'SFR100_ERR', 'B100_ERR', 'B1000_ERR', 'EWOII_ERR', 'EWOIIIHB_ERR', 'EWNIIHA_ERR']
kcor_col = ['ISEDFIT_ID', 'Z', 'MAGGIES', 'IVARMAGGIES', 'CHI2', 'FLAM_1500', 'CFLUX_3727', 'KCORRECT', 'ABSMAG', 'IVARABSMAG', 'SYNTH_ABSMAG', 'ABSMAG_FILTERLIST']
new_col = ['ISEDFIT_ID', 'RA', 'DEC', 'Z', 'CHI2_SED', 'MSTAR', 'AGE', 'SFRAGE', 'TAU', 'ZMETAL', 'AV', 'MU', 'MSTAR_ERR', 'AGE_ERR', 'SFRAGE_ERR', 'TAU_ERR', 'ZMETAL_ERR', 'AV_ERR', 'MU_ERR', 'CHI2_KCORR', 'MAGGIES_G', 'MAGGIES_R', 'MAGGIES_I', 'MAGGIES_Z', 'MAGGIES_Y', 'IVARMAGGIES_G', 'IVARMAGGIES_R', 'IVARMAGGIES_I', 'IVARMAGGIES_Z', 'IVARMAGGIES_Y', 'BESTMAGGIES_G', 'BESTMAGGIES_R', 'BESTMAGGIES_I', 'BESTMAGGIES_Z', 'BESTMAGGIES_Y', 'KCORRECT_G', 'KCORRECT_R', 'KCORRECT_I', 'KCORRECT_Z', 'KCORRECT_Y', 'ABSMAG_G', 'ABSMAG_R', 'ABSMAG_I', 'ABSMAG_Z', 'ABSMAG_Y', 'IVARABSMAG_G', 'IVARABSMAG_R', 'IVARABSMAG_I', 'IVARABSMAG_Z', 'IVARABSMAG_Y', 'SYNTH_ABSMAG_G', 'SYNTH_ABSMAG_R', 'SYNTH_ABSMAG_I', 'SYNTH_ABSMAG_Z', 'SYNTH_ABSMAG_Y', 'LOGLUM_I_OBS', 'LOGM2L_I_OBS', 'LOGLUM_G_OBS', 'LOGM2L_G_OBS', 'LOGLUM_R_OBS', 'LOGM2L_R_OBS', 'LOGLUM_Z_OBS', 'LOGM2L_Z_OBS', 'LOGLUM_I_KCOR', 'LOGM2L_I_KCOR', 'LOGLUM_G_KCOR', 'LOGM2L_G_KCOR', 'LOGLUM_R_KCOR', 'LOGM2L_R_KCOR', 'LOGLUM_Z_KCOR', 'LOGM2L_Z_KCOR']
print(results_col)
print("---------------------------------------------------------------------------------------------------------")
print(kcor_col)
print("---------------------------------------------------------------------------------------------------------")
print(new_col)
In [6]:
# Location of the SED results
#sedDir = os.path.join(os.getenv("HOME"), 'astro4/massive/dr15b/isedfit')
sedDir = os.path.join(os.getenv("HOME"), 'work/isedfit')
resDir = os.path.join(sedDir, 'result')
figDir = os.path.join(sedDir, 'figure')
# Available models
modNow = [os.path.basename(i).replace('_post.fits.gz', '') for i in
glob.glob(resDir + "/*_post.fits.gz")]
for ii, mod in enumerate(modNow):
print("# %d: %s" % (ii+1, mod))
In [7]:
# 2016-08-12
for ii, mod in enumerate(modNow):
print("## %d - Dealing with SED run : %s" % (ii+1, mod))
# Get the sedStr and sedMod
modSeg = mod.split('_')
if len(modSeg) is 11:
if 'bc03' in modSeg[4]:
sedStr = '_'.join(modSeg[0:6])
sedMod = '_'.join(modSeg[6:11])
elif 'fsps' in modSeg[4]:
sedStr = '_'.join(modSeg[0:5])
sedMod = '_'.join(modSeg[5:11])
else:
warnings.warn("XX Wrong Model? Please check! %s" % mod)
elif len(modSeg) is 12:
if 'bc03' in modSeg[4]:
sedStr = '_'.join(modSeg[0:7])
sedMod = '_'.join(modSeg[7:12])
elif 'fsps' in modSeg[4]:
sedStr = '_'.join(modSeg[0:6])
sedMod = '_'.join(modSeg[6:12])
else:
warnings.warn("XX Wrong Model? Please check! %s" % mod)
else:
warnings.warn("XX Wrong Model? Please check! %s" % mod)
#print(sedStr, sedMod)
# Combine the results
sedRes = combineSedResults(sedStr, sedMod, outDir=sedDir,
sedDir=resDir, kcorZ='0.1')
# Make summary plots
if 'redbcg' in modSeg[1]:
sedSummaryPlot(sedRes, sedStr, cmapUse=ORG4, alphaUse=0.75,
maxSizeUse=300, outDir=figDir)
elif 'redmem' in modSeg[1]:
sedSummaryPlot(sedRes, sedStr, cmapUse=GRN4, alphaUse=0.40,
maxSizeUse=140, outDir=figDir)
elif 'nonbcg' in modSeg[1]:
sedSummaryPlot(sedRes, sedStr, cmapUse=BLU5, alphaUse=0.25,
maxSizeUse=100, outDir=figDir)
In [8]:
# For reference model
## redBCG
redbcgRef = combineSedResults('dr1_redbcg_use_sed5b_fsps2',
'fsps_v2.4_miles_chab_charlot_sfhgrid01',
outDir=sedDir,
sedDir=resDir, kcorZ='0.1')
sedSummaryPlot(redbcgRef, 'redBCG', cmapUse=ORG4, alphaUse=0.75,
maxSizeUse=300, outDir=figDir)
## nonBCG
nonbcgRef = combineSedResults('dr1_nonbcg_use_sed5b_fsps2',
'fsps_v2.4_miles_chab_charlot_sfhgrid01',
outDir=sedDir,
sedDir=resDir, kcorZ='0.1')
sedSummaryPlot(nonbcgRef, 'nonBCG', cmapUse=BLU5, alphaUse=0.25,
maxSizeUse=100, outDir=figDir)
In [41]:
def isedfitPlot(tableUse, colName1, colName2, titleUse=None,
colColor='CHI2_SED', colSize='Z',
cmap=ORG4, maxSize=350, alpha=0.8,
ax=None, xylim=True, nolabel=False,
titleSize=20, labelSize=36,
titleX=0.04, titleY=0.92, vline=None,
xmargin=7.5, ymargin=7.5, noEdge=True,
xmin=None, ymin=None, xmax=None, ymax=None):
#---------------------------------------------------------#
if ax is None:
fig = plt.figure(figsize=(10, 10))
fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
ax1 = fig.add_subplot(111)
# Formatting
ax1 = songPlotSetup(ax1)
else:
ax1 = ax
# Scatter plot
if noEdge:
ax1.scatter(tableUse[colName1],
tableUse[colName2],
c=toColorArr(tableUse[colColor]),
s=toSizeArr(tableUse[colSize], maxSize=maxSize),
alpha=alpha, cmap=cmap, edgecolor='none', rasterized=True)
else:
ax1.scatter(tableUse[colName1],
tableUse[colName2],
c=toColorArr(tableUse[colColor]),
s=toSizeArr(tableUse[colSize], maxSize=maxSize),
alpha=alpha, cmap=cmap, rasterized=True)
# Label
if not nolabel:
colnameUse1 = colName1.replace('_', '\_')
colnameUse2 = colName2.replace('_', '\_')
ax1.set_xlabel('$\mathrm{%s}$' % colnameUse1, size=labelSize)
ax1.set_ylabel('$\mathrm{%s}$' % colnameUse2, size=labelSize)
# Axis limits
if xylim:
if xmin is None:
xmin = np.nanmin(tableUse[colName1])
if xmax is None:
xmax = np.nanmax(tableUse[colName1])
if ymin is None:
ymin = np.nanmin(tableUse[colName2])
if ymax is None:
ymax = np.nanmax(tableUse[colName2])
xmargin, ymargin = ((xmax - xmin) / xmargin), ((ymax - ymin) / ymargin)
ax1.set_xlim(xmin-xmargin, xmax+xmargin)
ax1.set_ylim(ymin-ymargin, ymax+ymargin)
# vline
if vline is not None:
ax1.axvline(vline, linestyle='--', linewidth=10, c='k', zorder=0,
dashes=(40,8), alpha=0.95)
# Title
if titleUse is not None:
titleUse = titleUse.replace('_', '\_')
ax1.text(titleX, titleY, '$\mathrm{%s}$' % titleUse, size=titleSize,
transform = ax1.transAxes)
#---------------------------------------------------------#
if ax is None:
return fig
else:
return ax1
#---------------------------------------------------------#
def isedMassPlot(sedRes, sedStr,
xSize=18, ySize=16,
massThresh=10.2, labelSize=26,
cmap=ORG4, alpha=0.75, maxSize=400,
tickMajor=0.4, tickMinor=0.1,
dpi=100, xmargin=7.5,
ymargin1=7.5, ymargin2=7.5,
ymargin3=6.0, ymargin4=6.0,
titleX=0.42, titleY=0.06, titleSize=50,
xLabelSize=50, yLabelSize=50,
saveFig=True, outDir='./',
titleUse=None, xmin=None, xmax=None,
ymin=None, ymax=None):
"""Summary plots of iSEDFit results (Stellar mass related)"""
#-----------------------------------------------------------------#
# Exclude problematic results
sedRes = sedRes[sedRes['MSTAR'] >= massThresh]
#-----------------------------------------------------------------#
fig = plt.figure(figsize=(xSize, ySize))
fig.subplots_adjust(left=0.12, right=0.90,
bottom=0.09, top=0.99,
wspace=0.02, hspace=0.00)
# logM - Age plot
ax1 = fig.add_subplot(221)
ax1 = songPlotSetup(ax1, xlabel=46, ylabel=46, border=7.0,
xtickFormat='$\mathbf{%4.2f}$',
ytickFormat='$\mathbf{%4.1f}$')
fig1 = isedfitPlot(sedRes, 'MSTAR', 'AGE',
titleUse=None,
colColor='Z', colSize='DMAG_I',
ax=ax1, nolabel=True, alpha=alpha,
cmap=cmap, maxSize=maxSize,
labelSize=labelSize,
xmargin=xmargin, ymargin=ymargin1,
xmin=xmin, ymin=ymin, xmax=xmax, ymax=ymax)
##
ax1.tick_params(axis='x', which='both', labelbottom='off')
ax1.set_ylabel('$\mathrm{Age\ (Gyr)}$', size=yLabelSize)
ax1.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
ax1.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
ax1.text(0.75, 0.14, '$\mathrm{Color:\ }z$', size=45,
transform=ax1.transAxes,
verticalalignment='bottom', horizontalalignment='center')
ax1.text(0.75, 0.04, '$\mathrm{Size:\ }{\Delta\ i}$', size=45,
transform=ax1.transAxes,
verticalalignment='bottom', horizontalalignment='center')
# logM - SFR Age plot
ax2 = fig.add_subplot(223)
ax2 = songPlotSetup(ax2, xlabel=46, ylabel=46, border=7.0,
xtickFormat='$\mathbf{%4.2f}$',
ytickFormat='$\mathbf{%4.1f}$')
fig2 = isedfitPlot(sedRes, 'MSTAR', 'LOGM2L_I_KCOR',
titleUse=None, cmap=cmap, alpha=alpha,
colColor='Z', colSize='DMAG_I',
ax=ax2, nolabel=True, maxSize=maxSize,
labelSize=labelSize,
xmargin=xmargin, ymargin=ymargin2,
xmin=xmin, ymin=0.15, xmax=xmax, ymax=0.59)
##
ax2.set_xlabel('$\mathrm{Redshift}$', size=xLabelSize)
ax2.set_ylabel('$\log\ (M_{\star}/L_{\star})_{i,\ \mathrm{Kcor}}$', size=yLabelSize)
ax2.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
ax2.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
# logM - ZMetal
ax3 = fig.add_subplot(222)
ax3 = songPlotSetup(ax3, xlabel=46, ylabel=46, border=7.0,
xtickFormat='$\mathbf{%4.2f}$',
ytickFormat='$\mathbf{%4.1f}$')
fig3 = isedfitPlot(sedRes, 'MSTAR', 'ZMET_REL',
titleUse=titleUse, cmap=cmap, alpha=alpha,
titleX=titleX, titleY=titleY, titleSize=titleSize,
colColor='Z', colSize='DMAG_I',
ax=ax3, nolabel=True, maxSize=maxSize,
labelSize=labelSize,
xmargin=xmargin, ymargin=ymargin3,
xmin=xmin, ymin=ymin, xmax=xmax, ymax=ymax)
##
ax3.yaxis.set_label_position("right")
ax3.tick_params(axis='y', which='both', labelleft='off',
labelright='on', labelsize=46)
ax3.tick_params(axis='x', which='both', labelbottom='off')
ax3.set_ylabel('$\mathrm{Z}_{\star}/\mathrm{Z}_{\odot}$', size=yLabelSize)
ax3.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
ax3.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
ax3.yaxis.set_major_locator(ticker.MultipleLocator(0.4))
ax3.yaxis.set_minor_locator(ticker.MultipleLocator(0.1))
# logM - Av
ax4 = fig.add_subplot(224)
ax4 = songPlotSetup(ax4, xlabel=46, ylabel=46, border=7.0,
xtickFormat='$\mathbf{%4.2f}$',
ytickFormat='$\mathbf{%4.1f}$')
fig4 = isedfitPlot(sedRes, 'MSTAR', 'AV',
titleUse=None, cmap=cmap, alpha=alpha,
colColor='Z', colSize='DMAG_I',
ax=ax4, nolabel=True, maxSize=maxSize,
labelSize=labelSize,
xmargin=xmargin, ymargin=ymargin4,
xmin=xmin, ymin=ymin, xmax=xmax, ymax=ymax)
##
ax4.tick_params(axis='y', which='both', labelsize=46,
labelleft='off', labelright='on')
ax4.yaxis.set_label_position("right")
ax4.set_xlabel('$\log\ (M_{\star}/M_{\odot})$', size=xLabelSize)
ax4.set_ylabel('$\mathrm{A}_{V}\ \mathrm{(mag)}$', size=yLabelSize)
ax4.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
ax4.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
ax4.yaxis.set_major_locator(ticker.MultipleLocator(0.04))
ax4.yaxis.set_minor_locator(ticker.MultipleLocator(0.01))
#-----------------------------------------------------------------#
massPng = sedStr + '_logm_plots.pdf'
if saveFig:
fig.savefig(os.path.join(outDir, massPng), dpi=dpi)
plt.close(fig)
#-----------------------------------------------------------------#
def isedChi2Plot(sedRes, sedStr,
xSize=20, ySize=16,
massThresh=10.2, labelSize=26,
cmap=ORG4, alpha=0.75, maxSize=300,
tickMajor=None, tickMinor=None,
dpi=100, xmargin=7.5,
ymargin1=7.5, ymargin2=7.5,
ymargin3=6.0, ymargin4=6.0,
titleX=0.10, titleY=0.91, titleSize=35,
xLabelSize=40, yLabelSize=40,
saveFig=True, outDir='./'):
"""Summary plots of iSEDFit results (Chi^2 related)"""
#-----------------------------------------------------------------#
# Exclude problematic results
sedRes = sedRes[sedRes['MSTAR'] >= massThresh]
#-----------------------------------------------------------------#
fig = plt.figure(figsize=(xSize, ySize))
fig.subplots_adjust(left=0.10, right=0.90,
bottom=0.08, top=0.99,
wspace=0.02, hspace=0.00)
# logM - Age plot
ax1 = fig.add_subplot(221)
fig1 = isedfitPlot(sedRes, 'CHI2_SED', 'MSTAR',
titleUse=sedStr, titleSize=titleSize,
colColor='MSTAR', colSize='Z',
ax=ax1, nolabel=True, alpha=alpha,
titleX=titleX, titleY=titleY,
cmap=cmap, maxSize=maxSize,
labelSize=labelSize,
xmargin=xmargin, ymargin=ymargin1)
##
ax1.tick_params(axis='x', which='both', labelbottom='off')
ax1.set_ylabel('$\log\ (M_{\star}/M_{\odot})$', size=yLabelSize)
if tickMajor is not None:
ax1.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
if tickMinor is not None:
ax1.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
ax1.text(0.75, 0.14, '$\mathrm{Size:\ }z$', size=26,
transform = ax1.transAxes)
ax1.text(0.72, 0.06, '$\mathrm{Color:\ Mass}$', size=26,
transform = ax1.transAxes)
# logM - SFR Age plot
ax2 = fig.add_subplot(223)
fig2 = isedfitPlot(sedRes, 'CHI2_SED', 'Z',
titleUse=None, cmap=cmap, alpha=alpha,
colColor='MSTAR', colSize='Z',
ax=ax2, nolabel=True, maxSize=maxSize,
labelSize=labelSize,
xmargin=xmargin, ymargin=ymargin2)
##
ax2.set_xlabel('${\chi}^2$', size=xLabelSize)
ax2.set_ylabel('$\mathrm{Redshift}$', size=yLabelSize)
if tickMajor is not None:
ax2.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
if tickMinor is not None:
ax2.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
# logM - ZMetal
ax3 = fig.add_subplot(222)
fig3 = isedfitPlot(sedRes, 'CHI2_SED', 'AGE',
titleUse=None, cmap=cmap, alpha=alpha,
colColor='MSTAR', colSize='Z',
ax=ax3, nolabel=True, maxSize=maxSize,
labelSize=labelSize,
xmargin=xmargin, ymargin=ymargin3)
##
ax3.yaxis.set_label_position("right")
ax3.tick_params(axis='y', which='both', labelleft='off',
labelright='on', labelsize=labelSize)
ax3.tick_params(axis='x', which='both', labelbottom='off')
ax3.set_ylabel('$\mathrm{Age\ (Gyr)}$', size=yLabelSize)
if tickMajor is not None:
ax3.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
if tickMinor is not None:
ax3.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
# logM - Av
ax4 = fig.add_subplot(224)
fig4 = isedfitPlot(sedRes, 'CHI2_SED', 'AV',
titleUse=None, cmap=cmap, alpha=alpha,
colColor='MSTAR', colSize='Z',
ax=ax4, nolabel=True, maxSize=maxSize,
labelSize=labelSize,
xmargin=xmargin, ymargin=ymargin4)
##
ax4.tick_params(axis='y', which='both', labelsize=labelSize,
labelleft='off', labelright='on')
ax4.yaxis.set_label_position("right")
ax4.set_xlabel('${\chi}^2$', size=xLabelSize)
ax4.set_ylabel('$\mathrm{A}_{V}\ \mathrm{(mag)}$', size=yLabelSize)
if tickMajor is not None:
ax4.xaxis.set_major_locator(ticker.MultipleLocator(tickMajor))
if tickMinor is not None:
ax4.xaxis.set_minor_locator(ticker.MultipleLocator(tickMinor))
ax4.yaxis.set_major_locator(ticker.MultipleLocator(0.04))
ax4.yaxis.set_minor_locator(ticker.MultipleLocator(0.01))
#-----------------------------------------------------------------#
chi2Png = sedStr + '_chi2_plots.pdf'
if saveFig:
fig.savefig(os.path.join(outDir, chi2Png), dpi=dpi)
plt.close(fig)
#-----------------------------------------------------------------#
In [42]:
isedMassPlot(redbcgRef, 'redBCG', cmap=ORG4,
alpha=0.75, maxSize=400, titleUse='cenHighMh',
tickMajor=0.4, tickMinor=0.1, xmargin=7.5,
ymargin1=7.5, ymargin2=7.5,
ymargin3=6.0, ymargin4=6.0,
saveFig=True, outDir=figDir, xmin=11.10)
isedMassPlot(nonbcgRef, 'nonBCG', cmap=BLU5,
alpha=0.30, maxSize=90, titleUse='cenLowMh',
tickMajor=0.4, tickMinor=0.1, xmargin=7.5,
ymargin1=7.5, ymargin2=7.5,
ymargin3=6.0, ymargin4=6.0,
saveFig=True, outDir=figDir, xmin=11.10)
In [43]:
ii = 0
sedArr = []
for mod in modNow:
# Get the sedStr and sedMod
modSeg = mod.split('_')
if len(modSeg) is 11:
if 'bc03' in modSeg[4]:
sedStr = '_'.join(modSeg[0:6])
sedMod = '_'.join(modSeg[6:11])
elif 'fsps' in modSeg[4]:
sedStr = '_'.join(modSeg[0:5])
sedMod = '_'.join(modSeg[5:11])
else:
warnings.warn("XX Wrong Model? Please check! %s" % mod)
elif len(modSeg) is 12:
if 'bc03' in modSeg[4]:
sedStr = '_'.join(modSeg[0:7])
sedMod = '_'.join(modSeg[7:12])
elif 'fsps' in modSeg[4]:
sedStr = '_'.join(modSeg[0:6])
sedMod = '_'.join(modSeg[6:12])
else:
warnings.warn("XX Wrong Model? Please check! %s" % mod)
else:
warnings.warn("XX Wrong Model? Please check! %s" % mod)
sedArr.append(sedStr)
print("# %2d : %s" % (ii+1, sedStr))
ii+=1
In [14]:
redbcg = 'redbcg'
sedRef = 'dr1_%s_use_sed5b_fsps2' % redbcg
for sedRun in sedArr:
if (redbcg in sedRun) and (sedRun != sedRef):
label1 = '_'.join(sedRef.split('_')[4:])
label2 = '_'.join(sedRun.split('_')[4:])
# ------------------------------------------------------------------------------#
#try:
print("## Compare %s with %s" % (sedRun, sedRef))
sedMassCompare(sedRef, sedRun, redbcg, label1, label2,
location=sedDir, cmapUse=ORG4, colorUse=RED0,
alphaUse=0.75, maxSizeUse=300,
chi2Thr=10.0, outDir=figDir, sizeCol='DMAG_Y')
#except Exception:
# print("## Cannot compare %s and %s" % (sedRun, sedRef))