In [1]:
%pylab inline
from __future__ import division
import os
import copy
import argparse
import fnmatch
import numpy as np
import scipy
from scipy.interpolate import interp1d
# Astropy
from astropy.io import fits
from astropy import units as u
from astropy.stats import sigma_clip
# AstroML
from astroML.plotting import hist
# Matplotlib related
# Matplotlib default settings
rcdef = plt.rcParams.copy()
pylab.rcParams['figure.figsize'] = 12, 10
pylab.rcParams['xtick.major.size'] = 8.0
pylab.rcParams['xtick.major.width'] = 2.5
pylab.rcParams['xtick.minor.size'] = 4.0
pylab.rcParams['xtick.minor.width'] = 2.5
pylab.rcParams['ytick.major.size'] = 8.0
pylab.rcParams['ytick.major.width'] = 2.5
pylab.rcParams['ytick.minor.size'] = 4.0
pylab.rcParams['ytick.minor.width'] = 2.5
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from matplotlib.ticker import NullFormatter
from matplotlib.ticker import MaxNLocator
# Personal
import hscUtils as hUtil
import galSBP
import cPickle as pickle
amag_sun_g = 5.33
amag_sun_r = 4.67
amag_sun_i = 4.48
amag_sun_z = 4.42
In [ ]:
# Location of the cutouts
""" On Astro3 / For MacOS"""
baseDir = '/Users/songhuang/astro3/hscs'
redBCG = os.path.join(baseDir, 'redmapper')
nonBCG = os.path.join(baseDir, 'nonbcg')
In [2]:
def normProfile(rad, sbp, minRad, maxRad):
"""
Normalize the profile
"""
offset = np.nanmedian(sbp[(rad >= minRad) & (rad <= maxRad)])
return (sbp-offset)
In [3]:
def logAdd(para1, para2):
"""
Doc
"""
return np.log10((10.0 ** np.asarray(para1)) + (10.0 ** np.asarray(para2)))
In [4]:
def loadGalfitOutput(pklFile):
"""
Load a .pkl file for GALFIT model
Return the GalfitParser object
"""
if not os.path.isfile(pklFile):
raise Exception("XXX Can not find the .pkl file : %s") % pklFile
else:
return pickle.load(open(pklFile, 'rb'))
In [5]:
def findProfile(pattern, loc, verbose=False):
"""
Find the prefix of the ellipse profiles
"""
result = []
for root, dirs, files in os.walk(loc):
for name in files:
if fnmatch.fnmatch(name, pattern):
result.append(os.path.join(root, name))
if verbose:
print "### %d files found !" % len(result)
return result
In [6]:
def readProfile(ellFile):
"""
Load the pickle format 1-D profile
"""
if os.path.isfile(ellFile):
return pickle.load(open(ellFile, 'rb'))
else:
raise Exception("XXX Can not find the Ellipse file")
In [7]:
def getEllipProfile(objid, loc, stage='3', verbose=False):
pattern = '*' + str(objid) + '*_img_ellip_' + stage + '.pkl'
result = findProfile(pattern, loc)
if len(result) == 0:
#raise Exception('XXX Can not find the Ellipse profile ! %s' % pattern)
if verbose:
print 'XXX Can not find the Ellipse profile ! %s' % pattern
return None
elif len(result) == 1:
ellProf = readProfile(result[0])
return ellProf
else:
if verbose:
print "XXX Warning: More than 1 profiles found! Use the first one : %s" % result[0]
ellProf = readProfile(result[0])
return ellProf
In [1]:
def pixKpc(redshift, pix=0.168, show=True, npix=1.0):
"""
The corresponding Kpc size of a pixel
"""
pixKpc = pix * npix * hUtil.cosmoScale(redshift)
if show:
print "# %d pixel(s) = %6.3f Kpc" % (npix, pixKpc)
return pixKpc
In [9]:
def correctProf(ellProf, redshift, verbose=False, extinction=0.0, zp=27.0,
amag_sun=None, corCurve=False):
"""
DOC
"""
scale = hUtil.cosmoScale(redshift)
distmod = hUtil.cosmoDistMod(redshift)
if verbose:
print "### REDSHIFT : %6.4f" % redshift
print "### SCALE : %7.4f kpc/arcsec" % scale
print "### DISTMOD : %7.4f mag" % distmod
sma_kpc = ellProf['sma_asec'] * scale
if not corCurve:
abs_mag = -2.5 * np.log10(ellProf['growth_ori']) + zp - extinction - distmod
else:
abs_mag = -2.5 * np.log10(ellProf['growth_cor']) + zp - extinction - distmod
if amag_sun is not None:
abs_sbp = (amag_sun + 21.572 - (ellProf['sbp'] - extinction))/2.5 + 6.0
else:
abs_sbp = ellProf['sbp'] - extinction
return sma_kpc, abs_sbp, abs_mag
In [10]:
def photoCompare(sample, zp=27.0, prefix='hsc_xmm_log11.6', amag_sun=amag_sun_i,
pix=0.168, absMag=True, aStr='extinction_i', verbose=False,
magStr='cModelAbsMag_i', stage='3', plot=True, color='b',
radius=None, corCurve=False, idStr='objid_dr12',
massStr='logMass', zStr='z', extCorrect=True, ax=None,
marker='o'):
"""
DOC
"""
""" Read in the sample catalog """
catFile = prefix + '_' + sample + '.fits'
catData = fits.open(catFile)[1].data
""" Useful columns """
gal_id = catData[idStr]
gal_z = catData[zStr]
if extCorrect:
gal_alam = catData[aStr]
else:
gal_alam = gal_z * 0.0
gal_mag = catData[magStr]
gal_logm = catData[massStr]
""" The "total" magnitude from SBP """
mag_sbp = []
mag_sdss = []
for ii in range(len(catData)):
""" Get the Ellipse profile """
ellProf = getEllipProfile(gal_id[ii], sample, stage=stage)
if ellProf is not None:
""" Get the radius in Kpc, luminosity profile, and growth curve for absolute magnitude """
sma_kpc, abs_sbp, abs_cog = correctProf(ellProf, gal_z[ii], extinction=gal_alam[ii],
zp=27.0, amag_sun=amag_sun, corCurve=corCurve)
""" Distance module """
distmod = hUtil.cosmoDistMod(gal_z[ii])
if radius is None:
mag_use = ellProf['mag_tot'][0] - gal_alam[ii]
if absMag:
mag_use = mag_use - distmod
else:
rsma_use = radius ** 0.25
rsma_common = np.arange(0, 4.5, 0.1)
indexUse = np.isfinite(abs_cog)
intrpFunc = interp1d(sma_kpc[indexUse]**0.25, abs_cog[indexUse], kind='linear')
try:
mag_use = intrpFunc(rsma_use) - gal_alam[ii]
if not absMag:
mag_use = mag_use + distmod
except Exception:
print "### Interpolation range is bad : %d" % gal_id[ii]
mag_use = ellProf['mag_tot'][0] - gal_alam[ii]
if absMag:
mag_use = mag_use - distmod
mag_sbp.append(mag_use)
mag_sdss.append(gal_mag[ii])
if verbose:
print "### %22d -- %6.3f -- %6.3f -- %6.3f" % (gal_id[ii], gal_mag[ii],
ellProf['mag_tot'][0], distmod)
mag_sbp = np.asarray(mag_sbp)
mag_sdss = np.asarray(mag_sdss)
if plot:
if ax is not None:
ax.scatter(mag_sbp, mag_sdss, marker=marker, color=color)
else:
scatter(mag_sbp, mag_sdss, marker=marker, color=color)
return mag_sbp
In [33]:
def getStackProfile(sample, zp=27.0, prefix='hsc_xmm_log11.6', amag_sun=amag_sun_i,
pix=0.168, absMag=True, aStr='extinction_i', verbose=False,
magStr='cModelAbsMag_i', stage='3', plot=True, plotAvg=True,
color='b', width=0.9, singleColor='k',
radius=None, corCurve=False, useMedian=True, alpha=0.1,
minRsma=0.3, maxRsma=3.5, extCorrect=True, idStr='objid_dr12',
zStr='z', massStr='logMass', rsmaStep=0.1, ax=None, root=None,
idSuffix=None, dirSuffix=None, catFile=None):
"""
DOC
"""
""" Read in the sample catalog """
if catFile is None:
catFile = prefix + '_' + sample + '.fits'
if root is not None:
catFile = os.path.join(root, catFile)
if not os.path.isfile(catFile):
raise Exception('XXX Can not find the input catalog : %s', catFile)
else:
catData = fits.open(catFile)[1].data
""" Useful columns """
gal_id = catData[idStr]
gal_z = catData[zStr]
print "### Median Redshift : %5.3f" % np.nanmedian(gal_z)
if extCorrect:
try:
gal_alam = catData[aStr]
except Exception:
print "### Can not find the column : %s" % aStr
gal_alam = gal_z * 0.0
else:
gal_alam = gal_z * 0.0
if magStr is not None:
try:
gal_mag = catData[magStr]
print " ## Median Magnitude : %5.2f" % np.nanmedian(gal_mag)
except Exception:
print "### Can not find the column : %s" % magStr
gal_mag = gal_z * 0.0
else:
gal_mag = gal_z * 0.0
if massStr is not None:
try:
gal_mass = catData[massStr]
print " ## Median Log Stellar Mass : %5.2f" % np.nanmedian(gal_mass)
except Exception:
print "### Can not find the column : %s" % massStr
gal_mass = gal_z * 0.0
else:
gal_mass = gal_z * 0.0
""" Plot """
for ii in range(len(catData)):
if root is not None:
profDir = os.path.join(root, sample)
else:
profDir = sample
if idSuffix is None:
pattern = str(gal_id[ii])
else:
pattern = str(gal_id[ii]) + '*' + idSuffix
if dirSuffix is None:
location = profDir
else:
location = profDir + '_' + dirSuffix
ellProf = getEllipProfile(pattern, location, stage=stage)
if ellProf is not None:
sma_kpc, abs_sbp, abs_cog = correctProf(ellProf, gal_z[ii],
extinction=gal_alam[ii],
zp=27.0, verbose=False,
amag_sun=amag_sun,
corCurve=corCurve)
if plot:
if ax is not None:
ax.plot(sma_kpc**0.25, abs_sbp, singleColor, alpha=alpha, linewidth=width)
else:
plot(sma_kpc**0.25, abs_sbp, singleColor, alpha=alpha, linewidth=width)
intrpFunc = interp1d(sma_kpc**0.25, abs_sbp,
kind='linear')
try:
rsma_common = np.arange(minRsma, maxRsma, rsmaStep)
asbp_intrp = intrpFunc(rsma_common)
try:
asbp_stack = np.vstack((asbp_stack, asbp_intrp))
except NameError:
asbp_stack = asbp_intrp
except Exception as errMsg:
print errMsg
print "### Bad interpolation range: %d " % gal_id[ii]
print " %5.2f -- %5.2f" % (np.nanmin(sma_kpc**0.25),
np.nanmax(sma_kpc**0.25))
if useMedian:
asbp_avg = np.nanmedian(asbp_stack, axis=0)
else:
asbp_avg = np.nanmean(asbp_stack, axis=0)
if plotAvg:
if ax is not None:
print len(rsma_common)
print len(asbp_avg)
ax.plot(rsma_common, asbp_avg, color, linewidth=3.5)
else:
plot(rsma_common, asbp_avg, color, linewidth=3.5)
return rsma_common, asbp_avg, asbp_stack
In [38]:
#def photoCompare(sample, zp=27.0, prefix='hsc_xmm_log11.6', amag_sun=amag_sun_i,
# pix=0.168, absMag=True, aStr='extinction_i', verbose=False,
# magStr='cModelAbsMag_i', stage='3', plot=True, color='b',
# radius=None, corCurve=False, idStr='objid_dr12',
# massStr='logMass', zStr='z', extCorrect=True, ax=None,
# marker='o'):
fig = plt.figure(figsize=(10, 10))
fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
ax1 = fig.add_subplot(111)
maxR = 100.0
print "### Z1"
mag_sbp_z1 = photoCompare('z1', color='k', radius=maxR, ax=ax1)
print "### Z2"
mag_sbp_z2 = photoCompare('z2', color='r', radius=maxR, ax=ax1)
#mag_sbp_z2_b = photoCompare('z2', color='b', radius=100.0)
print "### Z3"
mag_sbp_z3 = photoCompare('z3', color='g', radius=maxR, ax=ax1)
print "### Z4"
mag_sbp_z4 = photoCompare('z4', color='b', radius=maxR, ax=ax1)
ax1.set_xlim(-21.5, -24.5)
ax1.set_ylim(-21.5, -24.5)
ax1.set_xlabel('Absolute Magnitude (HSC)', size=16.0)
ax1.set_ylabel('Absolute Magnitude (SDSS)', size=16.0)
magRange = np.arange(-25.0, -20.0, 0.2)
ax1.plot(magRange, magRange, 'k')
Out[38]:
In [43]:
fig = plt.figure(figsize=(10, 6))
fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
ax1 = fig.add_subplot(111)
hist(mag_sbp_z4, bins='knuth', histtype='stepfilled', alpha=0.5,
normed=True, color='b')
hist(mag_sbp_z3, bins='knuth', histtype='stepfilled', alpha=0.3,
normed=True, color='g')
hist(mag_sbp_z2, bins='knuth', histtype='step', alpha=0.8,
normed=True, color='r', linewidth=3.0)
hist(mag_sbp_z1, bins='knuth', histtype='step', alpha=0.8,
normed=True, color='k', linewidth=3.0)
ax1.set_xlabel('Absolute Magnitude', size=19.0)
ax1.set_ylabel('Fraction', size=19.0)
ax1.set_xlim(-21.5, -24.5)
for tick in ax1.xaxis.get_major_ticks():
tick.label.set_fontsize(19)
for tick in ax1.yaxis.get_major_ticks():
tick.label.set_fontsize(19)
Out[43]:
In [544]:
#def photoCompare(sample, zp=27.0, prefix='hsc_xmm_log11.6', amag_sun=amag_sun_i,
# pix=0.168, absMag=True, aStr='extinction_i', verbose=False,
# magStr='cModelAbsMag_i', stage='3', plot=True, color='b',
# radius=None, corCurve=False, idStr='objid_dr12',
# massStr='logMass', zStr='z', extCorrect=True, ax=None,
# marker='o'):
R1 = 20.0
R2 = 50.0
R3 = 100.0
R4 = 120.0
sample = 'z2'
print "### R1"
mag_sbp_r1 = photoCompare(sample, radius=R1, plot=None)
print "### R2"
mag_sbp_r2 = photoCompare(sample, radius=R2, plot=None)
print "### R3"
mag_sbp_r3 = photoCompare(sample, radius=R3, plot=None)
#print "### R4"
#mag_sbp_r4 = photoCompare(sample, radius=R4, plot=None)
In [555]:
fig = plt.figure(figsize=(10, 10))
fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
ax1 = fig.add_subplot(111)
magRange = np.arange(-25.0, -20.0, 0.2)
ax1.plot(magRange, magRange, 'k', linestyle='--', linewidth=3.0, alpha=0.8)
ax1.scatter(mag_sbp_r1, mag_sbp_r2, marker='o', color='g', s=75,
alpha=0.7)
ax1.scatter(mag_sbp_r1, mag_sbp_r3, marker='o', color='r', s=75,
alpha=0.5)
#ax1.scatter(mag_sbp_r1, mag_sbp_r4, marker='o', color='b')
ax1.set_xlim(-22.2, -24.3)
ax1.set_ylim(-22.2, -24.3)
ax1.set_xlabel('Absolute Magnitude (within 20 Kpc)', size=25.0)
ax1.set_ylabel('Absolute Magnitude (within larger aperture)', size=25.0)
for tick in ax1.xaxis.get_major_ticks():
tick.label.set_fontsize(20)
for tick in ax1.yaxis.get_major_ticks():
tick.label.set_fontsize(20)
ax1.spines['top'].set_linewidth(3.5)
ax1.spines['right'].set_linewidth(3.5)
ax1.spines['bottom'].set_linewidth(3.5)
ax1.spines['left'].set_linewidth(3.5)
ax1.text(-23.45, -23.0, '$0.3<z<0.4$', size=32.0)
ax1.text(-23.25, -22.8, 'log $(M_{s}/M_{\odot})>11.6$', size=32.0)
ax1.text(-23.35, -22.55, 'within 50 Kpc', color='g', size=32.0)
ax1.text(-23.30, -22.35, 'within 100 Kpc', color='r', size=32.0)
Out[555]:
In [559]:
fig = plt.figure(figsize=(12, 7.5))
fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
ax1 = fig.add_subplot(111)
hist(mag_sbp_r1, bins='knuth', histtype='stepfilled', alpha=0.5,
normed=True, color='b')
hist(mag_sbp_r2, bins='knuth', histtype='stepfilled', alpha=0.3,
normed=True, color='g')
hist(mag_sbp_r3, bins='knuth', histtype='stepfilled', alpha=0.6,
normed=True, color='r', linewidth=3.0)
#hist(mag_sbp_r4, bins='knuth', histtype='step', alpha=0.8,
# normed=True, color='k', linewidth=3.0)
ax1.set_xlabel('Magnitude', size=30.0)
ax1.set_ylabel('Fraction', size=30.0)
ax1.set_xlim(-22.01, -24.4)
ax1.set_ylim(0.01, 1.45)
for tick in ax1.xaxis.get_major_ticks():
tick.label.set_fontsize(28)
for tick in ax1.yaxis.get_major_ticks():
tick.label.set_fontsize(28)
ax1.spines['top'].set_linewidth(3.5)
ax1.spines['right'].set_linewidth(3.5)
ax1.spines['bottom'].set_linewidth(3.5)
ax1.spines['left'].set_linewidth(3.5)
ax1.text(-22.1, 1.2, '$<20$ Kpc', size=32, color='b')
ax1.text(-22.1, 1.0, '$<50$ Kpc', size=32, color='g')
ax1.text(-22.1, 0.8, '$<100$ Kpc', size=32, color='r')
Out[559]:
In [36]:
fig = plt.figure(figsize=(12, 12))
fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
ax1 = fig.add_subplot(111)
rsma_common_z1, asbp_avg_z1, asbp_stack_z1 = getStackProfile('z1', minRsma=0.1, maxRsma=3.8,
ax=ax1, alpha=0.1)
rsma_common_z1m, asbp_avg_z1m, asbp_stack_z1m = getStackProfile('z1', prefix='hsc_xmm_mem', minRsma=0.1,
maxRsma=3.8, ax=ax1, alpha=0.7, plotAvg=False,
width=3.0, singleColor='g')
rsma_common_z1b, asbp_avg_z1b, asbp_stack_z1b = getStackProfile('z1', prefix='hsc_xmm_bcg', minRsma=0.1,
maxRsma=3.8, ax=ax1, alpha=0.95, plotAvg=False,
width=3.0, singleColor='r')
print asbp_stack_z1.shape
std_stack_z1 = np.nanstd(asbp_stack_z1, axis=0)
#ax1.fill_between(rsma_common_z1, asbp_avg_z1+std_stack_z1, asbp_avg_z1-std_stack_z1,
# facecolor='b', alpha=0.2)
ax1.text(2.4, 8.6, '$0.15<z<0.30$', size=42)
ax1.set_xlabel('$R^{1/4}$ (Kpc)', size=32)
ax1.set_ylabel('log $(\mu/[L_{\odot} Kpc^{-2}])$', size=32)
ax1.set_xlim(0.0, 4.2)
ax1.set_ylim(3.2, 9.5)
for tick in ax1.xaxis.get_major_ticks():
tick.label.set_fontsize(30)
for tick in ax1.yaxis.get_major_ticks():
tick.label.set_fontsize(30)
ax1.spines['top'].set_linewidth(3.5)
ax1.spines['right'].set_linewidth(3.5)
ax1.spines['bottom'].set_linewidth(3.5)
ax1.spines['left'].set_linewidth(3.5)
In [30]:
fig = plt.figure(figsize=(12, 12))
fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
ax1 = fig.add_subplot(111)
rsma_common_z2, asbp_avg_z2, asbp_stack_z2 = getStackProfile('z2', minRsma=0.1, maxRsma=3.8,
ax=ax1, alpha=0.1)
rsma_common_z2m, asbp_avg_z2m, asbp_stack_z2m = getStackProfile('z2', prefix='hsc_xmm_mem', minRsma=0.1,
maxRsma=3.8, ax=ax1, alpha=0.5, plotAvg=False,
width=3.0, singleColor='g')
rsma_common_z2b, asbp_avg_z2b, asbp_stack_z2b = getStackProfile('z2', prefix='hsc_xmm_bcg', minRsma=0.1,
maxRsma=3.8, ax=ax1, alpha=0.9, plotAvg=False,
width=3.0, singleColor='r')
print asbp_stack_z2.shape
std_stack_z2 = np.nanstd(asbp_stack_z2, axis=0)
#ax1.fill_between(rsma_common_z2, asbp_avg_z2+std_stack_z2, asbp_avg_z2-std_stack_z2,
# facecolor='b', alpha=0.2)
ax1.text(2.4, 8.6, '$0.30<z<0.40$', size=42)
ax1.set_xlabel('$R^{1/4}$ (Kpc)', size=32)
ax1.set_ylabel('log $(\mu/[L_{\odot} Kpc^{-2}])$', size=32)
ax1.set_xlim(0.0, 4.2)
ax1.set_ylim(3.2, 9.5)
for tick in ax1.xaxis.get_major_ticks():
tick.label.set_fontsize(30)
for tick in ax1.yaxis.get_major_ticks():
tick.label.set_fontsize(30)
ax1.spines['top'].set_linewidth(3.5)
ax1.spines['right'].set_linewidth(3.5)
ax1.spines['bottom'].set_linewidth(3.5)
ax1.spines['left'].set_linewidth(3.5)
In [34]:
fig = plt.figure(figsize=(12, 12))
fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
ax1 = fig.add_subplot(111)
rsma_common_z3, asbp_avg_z3, asbp_stack_z3 = getStackProfile('z3', minRsma=0.1, maxRsma=3.8,
ax=ax1, alpha=0.03)
rsma_common_z3m, asbp_avg_z3m, asbp_stack_z3m = getStackProfile('z3', prefix='hsc_xmm_mem', minRsma=0.1,
maxRsma=3.6, ax=ax1, alpha=0.5, plotAvg=False,
width=3.0, singleColor='g')
rsma_common_z3b, asbp_avg_z3b, asbp_stack_z3b = getStackProfile('z3', prefix='hsc_xmm_bcg', minRsma=0.1,
maxRsma=3.6, ax=ax1, alpha=0.9, plotAvg=False,
width=3.0, singleColor='r')
print asbp_stack_z3.shape
std_stack_z3 = np.nanstd(asbp_stack_z3, axis=0)
#ax1.fill_between(rsma_common_z3, asbp_avg_z3+std_stack_z3, asbp_avg_z3-std_stack_z3,
# facecolor='b', alpha=0.2)
ax1.text(2.4, 8.6, '$0.40<z<0.50$', size=42)
ax1.set_xlabel('$R^{1/4}$ (Kpc)', size=32)
ax1.set_ylabel('log $(\mu/[L_{\odot} Kpc^{-2}])$', size=32)
ax1.set_xlim(0.0, 4.2)
ax1.set_ylim(3.2, 9.5)
for tick in ax1.xaxis.get_major_ticks():
tick.label.set_fontsize(30)
for tick in ax1.yaxis.get_major_ticks():
tick.label.set_fontsize(30)
ax1.spines['top'].set_linewidth(3.5)
ax1.spines['right'].set_linewidth(3.5)
ax1.spines['bottom'].set_linewidth(3.5)
ax1.spines['left'].set_linewidth(3.5)