In [8]:
% matplotlib inline
from __future__ import (division,
print_function)
import os
import sys
import copy
import fnmatch
import warnings
import collections
import numpy as np
import scipy
try:
from scipy.stats import scoreatpercentile
except:
scoreatpercentile = False
from scipy.interpolate import interp1d
import cPickle as pickle
# Astropy
from astropy.io import fits
from astropy import units as u
from astropy.stats import sigma_clip
from astropy.table import Table, Column
from astropy.utils.console import ProgressBar
# AstroML
from astroML.plotting import hist
from astroML.density_estimation import KNeighborsDensity
try:
from sklearn.neighbors import KernelDensity
use_sklearn_KDE = True
except:
import warnings
warnings.warn("KDE will be removed in astroML version 0.3. Please "
"upgrade to scikit-learn 0.14+ and use "
"sklearn.neighbors.KernelDensity.", DeprecationWarning)
from astroML.density_estimation import KDE
use_sklearn_KDE = False
from sklearn.neighbors import KDTree
from sklearn.neighbors import BallTree
# Matplotlib related
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.patches import Ellipse
from matplotlib.ticker import NullFormatter
from matplotlib.ticker import MaxNLocator
from matplotlib.collections import PatchCollection
# Matplotlib default settings
rcdef = plt.rcParams.copy()
rcdef['figure.figsize'] = 12, 10
rcdef['xtick.major.size'] = 8.0
rcdef['xtick.major.width'] = 2.5
rcdef['xtick.minor.size'] = 4.0
rcdef['xtick.minor.width'] = 2.5
rcdef['ytick.major.size'] = 8.0
rcdef['ytick.major.width'] = 2.5
rcdef['ytick.minor.size'] = 4.0
rcdef['ytick.minor.width'] = 2.5
# Personal
import hscUtils as hUtil
import galSBP
import coaddCutoutGalfitSimple as gSimple
# Cosmology
import cosmology
c=cosmology.Cosmo(H0=70.0, omega_m=0.3, omega_l=0.7, flat=1)
# 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
# Personal tools
from hscUtils import songPlotSetup, removeIsNullCol
from hscUtils import confidence_interval, ma_confidence_interval_1d, confidence_interval_1d
## Constants
# 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
# 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]
EC1 = [0.135, 0.066, 0.862, 0.30]
EC2 = [0.135, 0.366, 0.862, 0.30]
EC3 = [0.135, 0.666, 0.862, 0.30]
REC = [0.12, 0.11, 0.87, 0.87]
# Universal RSMA ar
RSMA_COMMON = np.arange(0.4, 4.2, 0.02)
# Color
BLUE0 = "#92c5de"
BLUE1 = "#0571b0"
RED0 = "#f4a582"
RED1 = "#ca0020"
PURPLE0 = '#af8dc3'
PURPLE1 = '#762a83'
BROWN0 = '#bf812d'
BROWN1 = '#543005'
GREEN0 = '#7fbf7b'
GREEN1 = '#1b7837'
# 3-sigma
SIGMA1 = 0.3173
SIGMA2 = 0.0455
SIGMA3 = 0.0027
In [2]:
def normProf(sma, sbp, minSma, maxSma, integrate=True,
divide=False):
"""
Naive method to normalize the profile.
Parameters:
sbp : Array for surface brightness profile
sma : Radius range
minSma : Minimum SMA
maxSma Maximum SMA
"""
if integrate:
indInt = np.where(sma ** 4.0 <= minSma)
isoArea = (np.pi * ((sma[indInt] ** 4.0) ** 2.0))
isoRing = np.append(isoArea[1], [isoArea[1:] - isoArea[:-1]])
intNorm = (np.log10(np.nansum((10.0 ** sbp[indInt]) * isoRing)))
else:
intNorm = np.nanmedian(sbp[(sma >= minSma) &
(sma <= maxSma)])
if divide:
return (sbp / intNorm)
else:
return (sbp - intNorm)
def pixKpc(redshift, pix=0.168, show=True, npix=1.0):
"""
Get the corresponding Kpc size of a pixel.
Parameters:
"""
pixKpc = pix * npix * hUtil.cosmoScale(redshift)
if show:
print("# %d pixel(s) = %6.3f Kpc" % (npix, pixKpc))
return pixKpc
def logAdd(para1, para2):
""" Useful for adding magnitudes. """
return np.log10((10.0 ** np.asarray(para1)) +
(10.0 ** np.asarray(para2)))
def errAdd(err1, err2):
"""Add error quadral..."""
return np.sqrt((err1 ** 2.0) +
(err2 ** 2.0))
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 getLuminosity(mag, redshift, extinction=None,
amag_sun=None):
"""Get the absolute magnitude or luminosity."""
distmod = hUtil.cosmoDistMod(redshift)
absMag = (mag - distmod)
if extinction is not None:
absMag -= extinction
if amag_sun is not None:
absMag = ((amag_sun - absMag) / 2.5)
return absMag
def getStackProfiles(sample, loc, name='GAMA',
idCol='ID_USE', tabCol='sum_tab', save=True):
"""Get the stacks of the profiles."""
print("## Sample %s : Will deal with %d galaxies" % (name, len(sample)))
profiles = []
with ProgressBar(len(sample), ipython_widget=True) as bar:
for g in sample:
try:
gFile = os.path.join(loc, g['sum_tab'].replace('./', '')).strip()
gProf = Table.read(gFile, format='fits')
""" Add extra information """
try:
gProf.meta['KCORRECT_I'] = g['KCORRECT_I']
gProf.meta['KCORRECT_G'] = g['KCORRECT_G']
gProf.meta['KCORRECT_R'] = g['KCORRECT_R']
gProf.meta['KCORRECT_Z'] = g['KCORRECT_Z']
gProf.meta['KCORRECT_Y'] = g['KCORRECT_Y']
gProf.meta['LOGM2LI'] = g['LOGM2L_I_OBS']
gProf.meta['LUM_10'] = g['lum_10']
gProf.meta['LUM_100'] = g['lum_100']
gProf.meta['LUM_120'] = g['lum_120']
except Exception:
print("## WARNING: Some metadata may not be available !")
continue
except Exception:
print("## Missing: %s" % gFile)
continue
profiles.append(gProf)
bar.update()
if save:
outPkl = os.path.join(loc, (name + '_profs.pkl'))
hUtil.saveToPickle(profiles, outPkl)
print("## Save %s to %s" % (name, outPkl))
return profiles
def organizeSbp(profiles, col1='muI1', col2='KCORRECT_I',
kind='sbp', norm=False, r1=9.9, r2=10.1, divide=False,
col3=None, col4=None, justStack=False, integrate=False,
sun1=SUN_G, sun2=SUN_R, normArr=None,
index=None, extCat=None):
""" Get the stack of individual profiels, and their med/avg. """
# Surface brightness profile / luminosity profiles
if kind.strip() == 'sbp':
if col2 is not None:
if norm:
stack = np.vstack(normProf(p['rKpc'],
np.asarray(p[col1] + (p.meta[col2] / 2.5)),
r1, r2, divide=divide, integrate=integrate)
for p in profiles)
else:
stack = np.vstack(np.asarray(p[col1] + (p.meta[col2] / 2.5))
for p in profiles)
else:
print("## NO KCORRECTION APPLIED !!")
if norm:
stack = np.vstack(normProf(p['rKpc'], p[col1],
r1, r2, divide=divide, integrate=integrate)
for p in profiles)
else:
stack = np.vstack(np.asarray(p[col1]) for p in profiles)
# Mass profiles
elif kind.strip() == 'mass':
if norm and (normArr is None):
stack = np.vstack(normProf(p['rKpc'],
np.asarray(p[col1] + p.meta[col2]),
r1, r2, divide=divide, integrate=integrate)
for p in profiles)
elif norm and (normArr is not None):
stack = np.vstack((np.asarray(p[col1] + p.meta[col2]) - normArr[i]) for (i, p)
in enumerate(profiles))
else:
stack = np.vstack(np.asarray(p[col1] + p.meta[col2]) for p in profiles)
# Color profiles
elif kind.strip() == 'color':
cSun = (sun1 - sun2)
if col3 is None or col4 is None:
print("## NO KCORRECTION APPLIED !!")
if norm:
stack = np.vstack(normProf(p['rKpc'],
np.asarray(cSun - 2.5 * (p[col1] - p[col2])),
r1, r2, divide=divide, integrate=integrate)
for p in profiles)
else:
stack = np.vstack(np.asarray(cSun - 2.5 *(p[col1] - p[col2])) for p in profiles)
else:
if norm:
stack = np.vstack(normProf(p['rKpc'],
np.asarray(cSun - 2.5 * (p[col1] - p[col2]) -
(p.meta[col3] - p.meta[col4])),
r1, r2, divide=divide, integrate=integrate)
for p in profiles)
else:
stack = np.vstack(np.asarray(cSun - 2.5 * (p[col1] - p[col2]) -
(p.meta[col3] - p.meta[col4]))
for p in profiles)
# Luminosity profiles
elif kind.strip() == 'lum':
if col2 is None:
stack = np.vstack(np.asarray(p[col1]) for p in profiles)
else:
stack = np.vstack(np.asarray(p[col1] - p.meta[col2]) for p in profiles)
else:
raise Exception("## WRONG KIND !!")
if index is not None:
stack = np.vstack(p[index] for p in stack)
if not justStack:
""" Get the median and 1-sigma confidence range """
medProf = confidence_interval(stack, axis=0, alpha=np.asarray([SIGMA1, 1.0]),
metric=np.nanmedian, numResamples=1000,
interpolate=True)
avgProf = confidence_interval(stack, axis=0, alpha=np.asarray([SIGMA1, 1.0]),
metric=np.nanmean, numResamples=1000,
interpolate=True)
stdProf = confidence_interval(stack, axis=0, alpha=np.asarray([SIGMA1, 1.0]),
metric=np.nanstd, numResamples=1000,
interpolate=True)
return stack, medProf, avgProf, stdProf
else:
return stack
def loadPkl(filename):
try:
import cPickle as pickle
except:
warnings.warn("## cPickle is not available!!")
import pickle
if os.path.isfile(filename):
pklFile = open(filename, 'rb')
data = pickle.load(pklFile)
pklFile.close()
return data
else:
warnings.warn("## Can not find %s, return None" % filename)
return None
def updateKcorrect(cat, zCol='z_use', magType='mag_cmodel',
filters=['g', 'r', 'i', 'z', 'y']):
"""Update the K-correction for each band."""
cat = copy.deepcopy(cat)
for f in filters:
magCol = f + magType
extCol = 'a_' + f
absCol = 'ABSMAG_' + f.upper()
kcorCol = 'KCORRECT_' + f.upper()
newKcor = getLuminosity(cat[magCol], cat[zCol],
extinction=cat[extCol]) - cat[absCol]
try:
cat[kcorCol] = newKcor
except Exception:
cat.add_column(Column(newKcor, name=kcorCol))
return cat
def updateMass(cat, m2l='LOGM2L_I_OBS',
apertures = ['5', '10', '15', '25', '30', '40',
'50', '60', '75', '100', '120',
'150', 'max']):
"""Update the stellar masses at each aperture."""
cat = copy.deepcopy(cat)
for aper in apertures:
try:
cat.add_column(Column(cat['lum_' + aper] + cat[m2l],
name=('logm_' + aper)))
except Exception:
print("## Can not update mass for lum_%s" % aper)
return cat
def kdtreeMatch(sample1, sample2, name='kdmatched',
mlim1=11.50, mlim2=12.0, zlim1=0.20, zlim2=0.50,
massCol='logm_100', zCol='z_use',
k1=5, k2=5, k3=5, k4=5, leaf=7,
lamLimit=None, lamCol='LAMBDA_CLUSTER',
pcenLimit=None, pcenCol='P_CEN_1',
massCut=None, zCut=None,
massMargin1=0.005, massMargin2=0.010,
zMargin1=-0.01, zMargin2=-0.005,
plot=True, save=True, folder=None,
unique=True, ballTree=False, metric='l1',
figX=12, figY=18,
mmin=11.21, mmax=12.3, zmin=0.19, zmax=0.51,
sizeCol='logm_10', minSize=10.4,
colorCol='LAMBDA_CLUSTER', minColor=20, maxColor=65,
xLabel='$\mathrm{Redshift}$',
yLabel='$\log\ (M_{\star}/M_{\odot})\ (100\ \mathrm{Kpc})$',
legend1='$\Lambda \leq 20\ \mathrm{Cen}$',
legend2='$\Lambda > %d\ \mathrm{Cen}$',
legend3='$\Lambda \leq 20\ \mathrm{Matched}$',
massKernel=0.06, zKernel=0.025,
prefix1='redBCG', prefix2='nonBCG',
mErrCol='MSTAR_ERR', mass2=None):
"""Match two samples using K-Dtree."""
# Sample1 used for matching (should be the smaller one, e.g. the redBCG)
sampleUse1 = sample1[(sample1[massCol] >= mlim1) &
(sample1[massCol] <= mlim2) &
(sample1[zCol] >= zlim1) &
(sample1[zCol] <= zlim2)]
# Additional parameter cut: By default is Lambda and P_CEN
if lamLimit is not None:
sampleUse1 = sampleUse1[(sampleUse1[lamCol] >= lamLimit)]
if pcenLimit is not None:
sampleUse1 = sampleUse1[(sampleUse1[pcenCol] >= pcenLimit)]
print("# Sample1 Size: ", len(sampleUse1))
## Sample2 used for matching (the larger sample)
sampleUse2 = sample2[(sample2[massCol] >= (mlim1 + massMargin1)) &
(sample2[massCol] <= (mlim2 + massMargin2)) &
(sample2[zCol] >= (zlim1 + zMargin1)) &
(sample2[zCol] <= (zlim2 + zMargin2))]
print("# sample Sample Size: ", len(sampleUse2))
## Isolate the parameters used for matching
if mass2 is None:
data2 = np.stack((np.asarray(sampleUse2[massCol]),
np.asarray(sampleUse2[zCol])), axis=-1)
else:
data2 = np.stack((np.asarray(sampleUse2[massCol]),
np.asarray(sampleUse2[mass2]),
np.asarray(sampleUse2[zCol])), axis=-1)
if not BallTree:
dataTree = KDTree(data2, leaf_size=leaf, metric=metric)
else:
dataTree = BallTree(data2, leaf_size=leaf, metric=metric)
if massCut is None and zCut is None:
if mass2 is None:
data1 = np.stack((np.asarray(sampleUse1[massCol]),
np.asarray(sampleUse1[zCol])), axis=-1)
else:
data1 = np.stack((np.asarray(sampleUse1[massCol]),
np.asarray(sampleUse2[mass2]),
np.asarray(sampleUse1[zCol])), axis=-1)
dist, indAll = dataTree.query(data1, k=k1)
indAll = indAll.ravel()
elif massCut is not None and zCut is None:
if mass2 is None:
data1a = np.stack((np.asarray(sampleUse1[sampleUse1[massCol] <= massCut][massCol]),
np.asarray(sampleUse1[sampleUse1[massCol] <= massCut][zCol])), axis=-1)
data1b = np.stack((np.asarray(sampleUse1[sampleUse1[massCol] > massCut][massCol]),
np.asarray(sampleUse1[sampleUse1[massCol] > massCut][zCol])), axis=-1)
else:
data1a = np.stack((np.asarray(sampleUse1[sampleUse1[massCol] <= massCut][massCol]),
np.asarray(sampleUse1[sampleUse1[massCol] <= massCut][mass2]),
np.asarray(sampleUse1[sampleUse1[massCol] <= massCut][zCol])), axis=-1)
data1b = np.stack((np.asarray(sampleUse1[sampleUse1[massCol] > massCut][massCol]),
np.asarray(sampleUse1[sampleUse1[massCol] > massCut][mass2]),
np.asarray(sampleUse1[sampleUse1[massCol] > massCut][zCol])), axis=-1)
dist1, ind1 = dataTree.query(data1a, k=k1)
dist2, ind2 = dataTree.query(data1b, k=k2)
indAll = np.hstack([ind1.ravel(), ind2.ravel()])
elif massCut is None and zCut is not None:
if mass2 is None:
data1a = np.stack((np.asarray(sampleUse1[sampleUse1[zCol] <= zCut][massCol]),
np.asarray(sampleUse1[sampleUse1[zCol] <= zCut][zCol])), axis=-1)
data1b = np.stack((np.asarray(sampleUse1[sampleUse1[zCol] > zCut][massCol]),
np.asarray(sampleUse1[sampleUse1[zCol] > zCut][zCol])), axis=-1)
else:
data1a = np.stack((np.asarray(sampleUse1[sampleUse1[zCol] <= zCut][massCol]),
np.asarray(sampleUse1[sampleUse1[zCol] <= zCut][mass2]),
np.asarray(sampleUse1[sampleUse1[zCol] <= zCut][zCol])), axis=-1)
data1b = np.stack((np.asarray(sampleUse1[sampleUse1[zCol] > zCut][massCol]),
np.asarray(sampleUse1[sampleUse1[zCol] > zCut][mass2]),
np.asarray(sampleUse1[sampleUse1[zCol] > zCut][zCol])), axis=-1)
dist1, ind1 = dataTree.query(data1a, k=k1)
dist2, ind2 = dataTree.query(data1b, k=k2)
indAll = np.hstack([ind1.ravel(), ind2.ravel()])
else:
if mass2 is None:
data1a = np.stack((np.asarray(sampleUse1[(sampleUse1[zCol] <= zCut) &
(sampleUse1[massCol] <= massCut)][massCol]),
np.asarray(sampleUse1[(sampleUse1[zCol] <= zCut) &
(sampleUse1[massCol] <= massCut)][zCol])), axis=-1)
data1b = np.stack((np.asarray(sampleUse1[(sampleUse1[zCol] <= zCut) &
(sampleUse1[massCol] > massCut)][massCol]),
np.asarray(sampleUse1[(sampleUse1[zCol] <= zCut) &
(sampleUse1[massCol] > massCut)][zCol])), axis=-1)
data1c = np.stack((np.asarray(sampleUse1[(sampleUse1[zCol] > zCut) &
(sampleUse1[massCol] <= massCut)][massCol]),
np.asarray(sampleUse1[(sampleUse1[zCol] > zCut) &
(sampleUse1[massCol] <= massCut)][zCol])), axis=-1)
data1d = np.stack((np.asarray(sampleUse1[(sampleUse1[zCol] > zCut) &
(sampleUse1[massCol] > massCut)][massCol]),
np.asarray(sampleUse1[(sampleUse1[zCol] > zCut) &
(sampleUse1[massCol] > massCut)][zCol])), axis=-1)
else:
data1a = np.stack((np.asarray(sampleUse1[(sampleUse1[zCol] <= zCut) &
(sampleUse1[massCol] <= massCut)][massCol]),
np.asarray(sampleUse1[(sampleUse1[zCol] <= zCut) &
(sampleUse1[massCol] <= massCut)][mass2]),
np.asarray(sampleUse1[(sampleUse1[zCol] <= zCut) &
(sampleUse1[massCol] <= massCut)][zCol])), axis=-1)
data1b = np.stack((np.asarray(sampleUse1[(sampleUse1[zCol] <= zCut) &
(sampleUse1[massCol] > massCut)][massCol]),
np.asarray(sampleUse1[(sampleUse1[zCol] <= zCut) &
(sampleUse1[massCol] > massCut)][mass2]),
np.asarray(sampleUse1[(sampleUse1[zCol] <= zCut) &
(sampleUse1[massCol] > massCut)][zCol])), axis=-1)
data1c = np.stack((np.asarray(sampleUse1[(sampleUse1[zCol] > zCut) &
(sampleUse1[massCol] <= massCut)][massCol]),
np.asarray(sampleUse1[(sampleUse1[zCol] > zCut) &
(sampleUse1[massCol] <= massCut)][mass2]),
np.asarray(sampleUse1[(sampleUse1[zCol] > zCut) &
(sampleUse1[massCol] <= massCut)][zCol])), axis=-1)
data1d = np.stack((np.asarray(sampleUse1[(sampleUse1[zCol] > zCut) &
(sampleUse1[massCol] > massCut)][massCol]),
np.asarray(sampleUse1[(sampleUse1[zCol] > zCut) &
(sampleUse1[massCol] > massCut)][mass2]),
np.asarray(sampleUse1[(sampleUse1[zCol] > zCut) &
(sampleUse1[massCol] > massCut)][zCol])), axis=-1)
dist1, ind1 = dataTree.query(data1a, k=k1)
dist2, ind2 = dataTree.query(data1b, k=k2)
dist3, ind3 = dataTree.query(data1c, k=k3)
dist4, ind4 = dataTree.query(data1d, k=k4)
indAll = np.hstack([ind1.ravel(), ind2.ravel(), ind3.ravel(), ind4.ravel()])
## Unique elements:
indUnique = np.unique(indAll)
print("# All and Unique Matched Sample", len(indAll), len(indUnique))
## Matched samples:
sampleMatchA = sampleUse2[indAll]
sampleMatchU = sampleUse2[indUnique]
## Save results:
if folder is None:
folder = './'
if save:
sampleUse1.write(os.path.join(folder, prefix1 + '_' + name + '.fits'),
format='fits', overwrite=True)
if unique:
sampleMatchU.write(os.path.join(folder, prefix2 + '_' + name + '.fits'),
format='fits', overwrite=True)
else:
sampleMatchA.write(os.path.join(folder, prefix2 + '_' + name + '.fits'),
format='fits', overwrite=True)
## Plot, now the plots are by default used for redBCG-nonBCG match
if plot:
fig1 = plt.figure(figsize=(figX, figY))
fig1.subplots_adjust(left=0.12, right=0.985, wspace=0.05,
bottom=0.05, top=0.995, hspace=0.24)
mx = np.linspace(mmin, mmax, 100)
zx = np.linspace(zmin, zmax, 100)
# Redshift - Mass plot
ax1 = fig1.add_subplot(311)
ax1 = songPlotSetup(ax1)
## Mass limits
ax1.fill_between([zlim1, zlim2], [mlim1, mlim1], [mlim2, mlim2],
facecolor='k', edgecolor='k', alpha=0.05, zorder=0)
## Sample2
p1 = ax1.scatter(sample2[zCol], sample2[massCol],
alpha=0.20, facecolor=BLUE0, edgecolor='none',
label=legend1)
## Sample1
p2 = ax1.scatter(sample1[zCol], sample1[massCol],
facecolor='none', s=((sample1[sizeCol] - minSize) * 100.0),
cmap=ORG4, alpha=0.60, marker='s', linewidth=2.0,
edgecolor=RED0, label=None)
if lamLimit is None:
lamLimUse = 20
else:
lamLimUse = lamLimit
p3 = ax1.scatter(sampleUse1[zCol], sampleUse1[massCol],
edgecolor=ORG4(0.8),
s=((sampleUse1[sizeCol] - minSize) * 250.0),
cmap=ORG4, alpha=1.00,
c=toColorArr(sampleUse1[colorCol], bottom=minColor, top=maxColor),
label=legend2 % lamLimUse, marker='s')
## Matched GAMA sample
p4 = ax1.scatter(sampleMatchU[zCol], sampleMatchU[massCol],
alpha=0.70, facecolor=BLUE1, edgecolor='none', s=90,
label=legend3)
## Legend
ax1.legend(loc=(0.70, 0.035), shadow=True, fancybox=True,
numpoints=1, fontsize=22, scatterpoints=1,
markerscale=0.9, borderpad=0.25, handletextpad=0.1)
legend = ax1.get_legend()
legend.legendHandles[1].set_color(ORG4(0.8))
legend.legendHandles[0].set_sizes([150])
legend.legendHandles[1].set_sizes([200])
legend.legendHandles[2].set_sizes([150])
## Label
ax1.set_xlabel(xLabel, size=28)
ax1.set_ylabel(yLabel, size=25)
## Axis limits
ax1.set_xlim(zmin, zmax)
ax1.set_ylim((mmin - 0.44), mmax)
# ------------------------------------------------------------------------------------------------------#
# Mass Plot
ax2 = fig1.add_subplot(312)
ax2 = songPlotSetup(ax2)
## KDE for sample 1
sampleMKde = KernelDensity(massKernel, kernel='gaussian')
sampleMKde.fit(sampleUse1[massCol][:, None])
mDens1 = np.exp(sampleMKde.score_samples(mx[:, None]))
## KDE for sample 2
sampleMKde = KernelDensity(massKernel, kernel='gaussian')
sampleMKde.fit(sampleMatchA[massCol][:, None])
mDens2A = np.exp(sampleMKde.score_samples(mx[:, None]))
sampleMKde = KernelDensity(massKernel, kernel='gaussian')
sampleMKde.fit(sampleMatchU[massCol][:, None])
mDens2U = np.exp(sampleMKde.score_samples(mx[:, None]))
## Histogram
aa, _, _ = hist(sampleUse1[massCol], bins='knuth', ax=ax2,
normed=True, zorder=1,
histtype='stepfilled', edgecolor='none', facecolor=RED0,
alpha=0.6, label=legend2 % lamLimUse)
bb, _, _ = hist(sampleMatchU[massCol], bins='knuth', ax=ax2,
normed=True, zorder=1,
histtype='stepfilled', edgecolor='none', facecolor=BLUE0,
alpha=0.4, label=legend3)
## Density plot
ax2.plot(mx, mDens1, '-', color='r', zorder=3, linewidth=6.0, alpha=0.7)
ax2.plot(mx, mDens2A, '--', color='b', zorder=3, linewidth=3.0, alpha=0.5)
ax2.plot(mx, mDens2U, '-', color='b', zorder=3, linewidth=4.0, alpha=0.7)
## X, Y Limits
ax2.set_xlim(mmin, mmax)
ylim = np.nanmax(np.hstack([aa, bb])) + 0.72
ax2.set_ylim(0.02, ylim)
## Legend
ax2.legend(loc=(0.650, 0.732), shadow=True, fancybox=True,
numpoints=1, fontsize=24, scatterpoints=1,
markerscale=1.2, borderpad=0.3, handletextpad=0.2)
## X, Y Lables
ax2.set_xlabel(yLabel, size=26)
ax2.set_ylabel('$\mathrm{Normalized\ \#}$', size=40)
## Highlight the median
### sample 1
ax2.plot([np.nanmedian(sampleUse1[massCol]),
np.nanmedian(sampleUse1[massCol])],
[ylim-0.35, ylim-0.02], linewidth=5.0, c='r')
ax2.plot([np.percentile(sampleUse1[massCol], 25),
np.percentile(sampleUse1[massCol], 25)],
[ylim-0.35, ylim-0.02], linewidth=5.0, c='r', linestyle=':')
ax2.plot([np.percentile(sampleUse1[massCol], 75),
np.percentile(sampleUse1[massCol], 75)],
[ylim-0.35, ylim-0.02], linewidth=5.0, c='r', linestyle=':')
### sample 2
ax2.plot([np.nanmedian(sampleMatchU[massCol]),
np.nanmedian(sampleMatchU[massCol])],
[ylim-0.70, ylim-0.37], linewidth=5.0, c='b')
ax2.plot([np.percentile(sampleMatchU[massCol], 25),
np.percentile(sampleMatchU[massCol], 25)],
[ylim-0.70, ylim-0.37], linewidth=5.0, c='b', linestyle=':')
ax2.plot([np.percentile(sampleMatchU[massCol], 75),
np.percentile(sampleMatchU[massCol], 75)],
[ylim-0.70, ylim-0.37], linewidth=5.0, c='b', linestyle=':')
ax2.yaxis.set_major_formatter(NullFormatter())
# ------------------------------------------------------------------------------------------------------#
# Redshift Plot
ax3 = fig1.add_subplot(313)
ax3 = songPlotSetup(ax3)
## KDE for Sample1
sampleMKde = KernelDensity(zKernel, kernel='gaussian')
sampleMKde.fit(sampleUse1[zCol][:, None])
zDens1 = np.exp(sampleMKde.score_samples(zx[:, None]))
## KDE for Sample2
sampleMKde = KernelDensity(zKernel, kernel='gaussian')
sampleMKde.fit(sampleMatchA[zCol][:, None])
zDens2A = np.exp(sampleMKde.score_samples(zx[:, None]))
sampleMKde = KernelDensity(zKernel, kernel='gaussian')
sampleMKde.fit(sampleMatchU[zCol][:, None])
zDens2U = np.exp(sampleMKde.score_samples(zx[:, None]))
## Histogram
aa, _, _ = hist(sampleUse1[zCol], bins='knuth', ax=ax3, normed=True, zorder=1,
histtype='stepfilled', edgecolor='none', facecolor=RED0,
alpha=0.6)
bb, _, _ = hist(sampleMatchU[zCol], bins='knuth', ax=ax3, normed=True, zorder=1,
histtype='stepfilled', edgecolor='none', facecolor=BLUE0,
alpha=0.4)
## Density plot
ax3.plot(zx, zDens1, '-', color='r', zorder=3, linewidth=6.0, alpha=0.7)
ax3.plot(zx, zDens2A, '--', color='b', zorder=3, linewidth=3.0, alpha=0.5)
ax3.plot(zx, zDens2U, '-', color='b', zorder=3, linewidth=4.0, alpha=0.7)
## X, Y Limits
ax3.set_xlim(zmin, zmax)
ylim = np.nanmax(np.hstack([aa, bb])) + 0.72
ax3.set_ylim(0.02, ylim)
## X, Y Lables
ax3.set_xlabel(xLabel, size=28)
ax3.set_ylabel('$\mathrm{Normalized\ \#}$', size=40)
## Highlight the median
### Sample 1
ax3.plot([np.nanmedian(sampleUse1[zCol]),
np.nanmedian(sampleUse1[zCol])],
[ylim-0.35, ylim-0.02], linewidth=5.0, c='r')
ax3.plot([np.percentile(sampleUse1[zCol], 25),
np.percentile(sampleUse1[zCol], 25)],
[ylim-0.35, ylim-0.02], linewidth=5.0, c='r', linestyle=':')
ax3.plot([np.percentile(sampleUse1[zCol], 75),
np.percentile(sampleUse1[zCol], 75)],
[ylim-0.35, ylim-0.02], linewidth=5.0, c='r', linestyle=':')
### Sample 2
ax3.plot([np.nanmedian(sampleMatchU[zCol]),
np.nanmedian(sampleMatchU[zCol])],
[ylim-0.70, ylim-0.37], linewidth=5.0, c='b')
ax3.plot([np.percentile(sampleMatchU[zCol], 25),
np.percentile(sampleMatchU[zCol], 25)],
[ylim-0.70, ylim-0.37], linewidth=5.0, c='b', linestyle=':')
ax3.plot([np.percentile(sampleMatchU[zCol], 75),
np.percentile(sampleMatchU[zCol], 75)],
[ylim-0.70, ylim-0.37], linewidth=5.0, c='b', linestyle=':')
ax3.yaxis.set_major_formatter(NullFormatter())
# ------------------------------------------------------------------------------------------------------#
plt.show()
fig1.savefig(name + '_a.png', dpi=200)
# plt.close(fig1)
## Return results
if unique:
return sampleUse1, sampleMatchU
else:
return sampleUse1, sampleMatchA
def plotMassGrowth(parent1, parent2,
sample1a, sample1b, sample1c,
sample2a, sample2b, sample2c,
col1='logm_100', col2='logm_100', col3='logm_10',
nResample=500, xSep1=11.55, xSep2=11.75, xSep3=11.95,
lamLimit=30, pcenLimit=0.7, sizeCol1=None, sizeNorm=0.1,
colorCol1='z_use', colorLow=0.2, colorUpp=0.55,
showBin1=True, showBin2=True,
showHist1=False, showHist2=True, yKernel=0.05,
xMin=11.35, xMax=12.39, yMin=0.01, yMax=0.79,
xLabel='$\log\ (M_{\star}/M_{\odot})\ (100\ \mathrm{Kpc})$',
yLabel='$\Delta(\log M{\star})_{\mathrm{100\ kpc}-\mathrm{10\ kpc}}$',
outPng='mass_growth', save=True, outlineBin1=True):
"""Plot logM v.s. mass growth."""
# Statistics of the two subsamples
med_1ax = confidence_interval(sample1a[col1], axis=0, alpha=[SIGMA1, 1.0],
metric=np.nanmedian, numResamples=nResample,
interpolate=True)
med_1ay = confidence_interval((sample1a[col2] - sample1a[col3]),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.nanmedian, numResamples=nResample,
interpolate=True)
med_2ax = confidence_interval(sample2a[col1], axis=0, alpha=[SIGMA1, 1.0],
metric=np.nanmedian, numResamples=nResample,
interpolate=True)
med_2ay = confidence_interval((sample2a[col2] - sample2a[col3]),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.nanmedian, numResamples=nResample,
interpolate=True)
med_1bx = confidence_interval(sample1b[col1], axis=0, alpha=[SIGMA1, 1.0],
metric=np.nanmedian, numResamples=nResample,
interpolate=True)
med_1by = confidence_interval((sample1b[col2] - sample1b[col3]),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.nanmedian, numResamples=nResample,
interpolate=True)
med_2bx = confidence_interval(sample2b[col1], axis=0, alpha=[SIGMA1, 1.0],
metric=np.nanmedian, numResamples=nResample,
interpolate=True)
med_2by = confidence_interval((sample2b[col2] - sample2b[col3]),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.nanmedian, numResamples=nResample,
interpolate=True)
# ------------------------------------------------------------------------------------ #
fig = plt.figure(figsize=(13, 9))
ax1 = plt.axes(recScat)
ax2 = plt.axes(recHist)
ax1 = songPlotSetup(ax1)
ax2 = songPlotSetup(ax2)
# ---------------------------------------------------------------------------
# Scatter plot
## Mass separation of two bins
ax1.axvline(xSep1, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)
ax1.axvline(xSep2, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)
ax1.axvline(xSep3, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)
## Horizontal line for 0.0
ax1.axhline(0.0, linewidth=4.5, linestyle='-', c='k', alpha=0.2)
# Parent samples
p1 = ax1.scatter(parent2[col1],
parent2[col2] - parent2[col3], s=20.0,
alpha=0.10, facecolor=BLUE0, edgecolor='none',
label='$\Lambda \leq 20\ \mathrm{Cen;\ All}$')
if sizeCol1 is None:
size1 = 200.0
else:
size1 = ((parent1[sizeCol1] - sizeNorm) * 600.0)
p3 = ax1.scatter(parent1[col1],
parent1[col2] - parent1[col3], edgecolor=ORG4(0.4),
s=size1, alpha=0.5, facecolor='none',
label='$\Lambda > %d\ \mathrm{Cen;\ All}$' % lamLimit,
marker='s', linewidth=1.5)
# Matched ones
p2 = ax1.scatter(sample2c[col1],
sample2c[col2] - sample2c[col3], s=40.0,
alpha=0.70, facecolor=BLUE0, edgecolor=BLUE0,
label='$\Lambda \leq 20\ \mathrm{Cen;\ Use}$')
if sizeCol1 is None:
size2 = 200.0
else:
size2 = ((parent1[sizeCol1] - sizeNorm) * 600.0)
if colorCol1 is None:
color1 = ORG4(1.0)
else:
color1 = toColorArr(sample1c[colorCol1], bottom=colorLow, top=colorUpp)
p4 = ax1.scatter(sample1c[col1],
sample1c[col2] - sample1c[col3], edgecolor=ORG4(0.6),
s=size2, cmap=ORG4, alpha=0.90,
c=color1,
label='$\Lambda > %d\ \mathrm{Cen;\ Use}$' % lamLimit,
marker='s')
# Median values
if showBin1:
ax1.errorbar(med_2ax[2], med_2ay[2], marker='+', ms=1, mec='k',
yerr=0.01, color=BLUE1, ecolor=BLUE1, capthick=3.5, capsize=8,
alpha=0.8, linewidth=4.0, fmt='h', elinewidth=2.0, label=None,
zorder=100)
ax1.errorbar(med_1ax[2], med_1ay[2], marker='+', ms=1, mec='k', linewidth=4.0,
yerr=0.02, mfc=RED1, ecolor=RED1, capthick=3.5, capsize=8,
alpha=0.8, fmt='h', elinewidth=2.0, label=None, zorder=100)
if outlineBin1:
ax1.scatter(med_2ax[2], med_2ay[2], marker='^', s=400, facecolor=BLUE1,
edgecolor='k', linewidth=3.0, zorder=102, alpha=0.9,
label=None)
ax1.scatter(med_1ax[2], med_1ay[2], marker='p', s=420, facecolor=RED1,
edgecolor='k', linewidth=3.0, zorder=102,
label=None)
else:
ax1.scatter(med_2ax[2], med_2ay[2], marker='^', s=400, facecolor=BLUE1,
edgecolor='none', linewidth=3.0, zorder=102, alpha=0.9,
label=None)
ax1.scatter(med_1ax[2], med_1ay[2], marker='p', s=420, facecolor=RED1,
edgecolor='none', linewidth=3.0, zorder=102,
label=None)
if showBin2:
ax1.errorbar(med_2bx[2], med_2by[2], marker='+', ms=1, mec='k',
yerr=0.01, color=BLUE1, ecolor=BLUE1, capthick=3.5, capsize=8,
alpha=0.8, fmt='h', elinewidth=2.0, label=None, zorder=100)
ax1.errorbar(med_1bx[2], med_1by[2], marker='+', ms=1, mec='k',
yerr=0.03, mfc=RED1, ecolor=RED1, capthick=3.5, capsize=8,
alpha=0.8, fmt='h', elinewidth=2.0, label=None, zorder=100)
ax1.scatter(med_2bx[2], med_2by[2], marker='^', s=420, facecolor=BLUE1,
edgecolor='k', linewidth=3.0, zorder=102, alpha=0.9,
label='$[\Lambda \leq 20]$')
ax1.scatter(med_1bx[2], med_1by[2], marker='p', s=420, facecolor=RED1,
edgecolor='k', linewidth=3.0, zorder=102,
label='$[\Lambda > %d]$' % lamLimit)
# Legend
ax1.legend(loc=(0.68, 0.025), shadow=True, fancybox=True,
numpoints=1, fontsize=18, scatterpoints=1,
markerscale=0.9, borderpad=0.25, handletextpad=0.1)
legend = ax1.get_legend()
legend.legendHandles[0].set_color(BLU5(0.2))
legend.legendHandles[0].set_alpha(0.5)
legend.legendHandles[2].set_color(BLU5(0.5))
legend.legendHandles[2].set_alpha(0.9)
legend.legendHandles[3].set_color(ORG4(0.8))
legend.legendHandles[0].set_sizes([25])
legend.legendHandles[2].set_sizes([80])
legend.legendHandles[3].set_sizes([200])
# Label
ax1.set_xlabel(xLabel, size=41)
ax1.set_ylabel(yLabel, size=41)
# Axis limits
ax1.set_xlim(xMin, xMax)
ax1.set_ylim(yMin, yMax)
# ---------------------------------------------------------------------------
# Histogram
ax2.set_ylim(ax1.get_ylim())
## Horizonatal line for 0
ax2.axhline(0.0, linewidth=4.5, linestyle='-', c='k', alpha=0.2)
# Parameters used for histograms
Y1a = sample1a[col2] - sample1a[col3]
Y2a = sample2a[col2] - sample2a[col3]
Y1b = sample1b[col2] - sample1b[col3]
Y2b = sample2b[col2] - sample2b[col3]
Y1c = sample2c[col2] - sample2c[col3]
Y2c = sample1c[col2] - sample1c[col3]
yy = np.linspace(0.0, 1.0, 200)
MKde = KernelDensity(yKernel, kernel='gaussian')
# Show underlying historgrams of combined sample
n, bins, patches = hist(Y1c, bins='knuth', ax=ax2, edgecolor='none',
orientation='horizontal', histtype='stepfilled',
color=BLUE0, alpha=0.70, normed=1)
n, bins, patches = hist(Y2c, bins='knuth', ax=ax2, edgecolor='none',
orientation='horizontal', histtype='stepfilled',
color=ORG4(0.6), alpha=0.40, normed=1, linewidth=4.0)
# KDE densities for bin1
if showHist1:
MKde.fit(Y1a[:, None])
MDens1 = np.exp(MKde.score_samples(yy[:, None]))
MKde.fit(Y2a[:, None])
MDens2 = np.exp(MKde.score_samples(yy[:, None]))
ax2.plot(MDens1, yy, '--', color=RED1, zorder=3, linewidth=5.0, alpha=0.7)
ax2.plot(MDens2, yy, '--', color=BLUE1, zorder=3, linewidth=5.0, alpha=0.7)
# KDE densities for bin2
if showHist2:
MKde.fit(Y1b[:, None])
MDens3 = np.exp(MKde.score_samples(yy[:, None]))
MKde.fit(Y2b[:, None])
MDens4 = np.exp(MKde.score_samples(yy[:, None]))
ax2.plot(MDens3, yy, '-', color=RED1, zorder=3, linewidth=5.0, alpha=0.7)
ax2.plot(MDens4, yy, '-', color=BLUE1, zorder=3, linewidth=5.0, alpha=0.7)
# Setup axis
ax2.tick_params(axis='x', which='minor', bottom='off', top='off')
ax2.yaxis.set_major_formatter(NullFormatter())
ax2.xaxis.set_major_formatter(NullFormatter())
plt.show()
if save:
fig.savefig(outPng + '.png', dpi=200)
return fig
lamLimit = 30
def plotMassProfile(profSample1, profSample2,
col1='muI1', col2='LOGM2LI', matchR=100.0,
norm=False, integrate=False, divide=True,
normR1=10.0, normR2=12.0, showAvg=False,
diffColor1=RED0, diffColor2=RED1, showLMask=False,
xmin=1.02, xmax=4.25, ymin=4.01, ymax=9.79,
dmin=-0.199, dmax=0.399,
vline1=10.0, vline2=100.0, alpha1=0.07, alpha2=0.30,
lw1=2.5, lw2=2.5, color1='k', color2=RED0,
highlight1=False, highlight2=True, lamLimit=30,
mass1=11.55, mass2=11.95, z1=0.2, z2=0.5,
label1="$\Lambda \leq 20;\ \mathrm{Cen}$",
label2="$\Lambda > %d;\ \mathrm{Cen}$" % lamLimit,
showInfo1=True, showInfo2=True, showLegend=True,
rPsfKpc=5.5, kpcArr=[2.0, 5.0, 10.0, 20.0, 50.0, 100.0, 200.0],
save=True, outPng='mass_prof'):
"""Plot the median mass profiles."""
# Median profiles
bm_sm, bm_mm, bm_amg, bm_stdm = organizeSbp(profSample1, col1=col1,
col2=col2, kind='mass',
norm=norm, r1=normR1, r2=normR2,
divide=divide, integrate=integrate)
gu_sm, gu_mm, gu_amg, gu_stdm = organizeSbp(profSample2, col1=col1,
col2=col2, kind='mass',
norm=norm, r1=normR1, r2=normR2,
divide=divide, integrate=integrate)
if showLMask:
# Larger mask
bm_sm_b, bm_mm_b, bm_amg_b, bm_stdm_b = organizeSbp(profSample1, col1=col1,
col2=col2, kind='mass',
norm=norm, r1=normR1, r2=normR2,
divide=divide, integrate=integrate)
gu_sm_b, gu_mm_b, gu_amg_b, gu_stdm_b = organizeSbp(profSample2, col1=col1,
col2=col2, kind='mass',
norm=norm, r1=normR1, r2=normR2,
divide=divide, integrate=integrate)
# Random mix sample
mixM_sm = np.vstack([gu_sm, bm_sm])
randM_sm = []
ii = 0
while ii < 2000:
mprof = np.nanmedian(mixM_sm[np.random.choice(len(mixM_sm), len(bm_sm),
replace=False)], axis=0)
randM_sm.append(mprof)
ii += 1
# Integration check
indInt = np.where((RSMA_COMMON ** 4.0) <= matchR)
isoAreaB = (np.pi * ((RSMA_COMMON[indInt] ** 4.0) ** 2.0))
isoRingB = np.append(isoAreaB[1], [isoAreaB[1:] - isoAreaB[:-1]])
isoAreaG = (np.pi * ((RSMA_COMMON[indInt] ** 4.0) ** 2.0))
isoRingG = np.append(isoAreaG[1], [isoAreaG[1:] - isoAreaG[:-1]])
print("# Sample1: ", np.log10(np.nansum((10.0 ** bm_amg[2][indInt]) * isoRingB)))
print("# Sample2: ", np.log10(np.nansum((10.0 ** gu_amg[2][indInt]) * isoRingG)))
# --------------------------------------------------------------------------------------- #
## Setup up figure
fig = plt.figure(figsize=(13, 18))
ax1 = plt.axes(SBP2)
ax1 = songPlotSetup(ax1)
ax3 = plt.axes(SBP1)
ax3 = songPlotSetup(ax3, ylabel=20)
# --------------------------------------------------------------------------------------- #
## Mark the two interesting radius
if highlight1:
ax1.axvline(vline1 ** 0.25, linewidth=5.5, c='k', linestyle='--',
zorder=0, alpha=0.5, dashes=(30, 6))
else:
ax1.axvline(vline1 ** 0.25, linewidth=4.0, c='k', linestyle='--',
zorder=0, alpha=0.2)
if highlight2:
ax1.axvline(vline2 ** 0.25, linewidth=5.5, c='k', linestyle='--',
zorder=0, alpha=0.5, dashes=(30, 6))
else:
ax1.axvline(vline2 ** 0.25, linewidth=4.0, c='k', linestyle='--',
zorder=0, alpha=0.2)
# --------------------------------------------------------------------------------------- #
## Individual profiles
for gg in gu_sm:
ax1.plot(RSMA_COMMON, gg, c=color1, alpha=alpha1, linewidth=lw1)
for bb in bm_sm:
ax1.plot(RSMA_COMMON, bb, c=color2, alpha=alpha2, linewidth=lw2)
# --------------------------------------------------------------------------------------- #
## Median profiles
ax1.fill_between(RSMA_COMMON, gu_mm[0], gu_mm[1],
facecolor=BLUE1, edgecolor='none', alpha=0.9,
zorder=1005, label=label1)
ax1.fill_between(RSMA_COMMON, bm_mm[0], bm_mm[1],
facecolor=RED1, edgecolor='none', alpha=0.9,
zorder=1005, label=label2)
ax1.plot(RSMA_COMMON, gu_mm[2], linestyle='-', linewidth=8.0,
c=BLUE1, alpha=1.0, zorder=1010)
ax1.plot(RSMA_COMMON, bm_mm[2], linestyle='-', linewidth=7.0,
c=RED1, alpha=1.0, zorder=1010)
# --------------------------------------------------------------------------------------- #
## Y Lables
if norm:
ax1.set_ylabel('$\mathrm{Normalized}\ \ \log ({\mu}_{\star}/[M_{\odot}\ \mathrm{Kpc}^{-2}])$',
size=41)
else:
ax1.set_ylabel('$\log ({\mu}_{\star}/[M_{\odot}\ \mathrm{Kpc}^{-2}])$', size=45)
## Remove the X-axis label
ax1.xaxis.set_major_formatter(NullFormatter())
# --------------------------------------------------------------------------------------- #
## X, Y limits
ax1.set_xlim(xmin, xmax)
ax1.set_ylim(ymin, ymax)
# --------------------------------------------------------------------------------------- #
## Region affected by PSF
# z = 0.4 : 1"=5.4 Kpc
# z = 0.5 : 1"=6.1 Kpc
ax1.fill_between([0.0, rPsfKpc ** 0.25],
[ymin * 0.9, ymin * 0.9], [ymax * 1.1, ymax * 1.1],
facecolor='k', edgecolor='k', alpha=0.05, zorder=0)
## Label the PSF region
ax1.text(0.05, 0.20, '$\mathrm{PSF}$', rotation='vertical',
verticalalignment='bottom', horizontalalignment='left',
fontsize=50.0, transform=ax1.transAxes, weight='bold',
color='k', alpha=0.4)
# --------------------------------------------------------------------------------------- #
## Legend
if showLegend:
ax1.legend(loc=(0.675, 0.675), shadow=True, fancybox=True,
numpoints=1, fontsize=30, scatterpoints=1,
markerscale=1.2, borderpad=0.3, handletextpad=0.2)
# --------------------------------------------------------------------------------------- #
## Information of the sample
if showInfo1:
ax1.text(0.97, 0.83, '$%5.2f < z < %5.2f$' % (z1, z2),
verticalalignment='bottom', horizontalalignment='right',
fontsize=40.0, transform=ax1.transAxes, weight='bold',
color='k', backgroundcolor='w')
if showInfo2:
matchStr = str(int(matchR)).strip()
ax1.text(0.97, 0.90,
'$%5.2f < \log\ (M_{\star, %s}/M_{\odot}) < %5.2f$' % (mass1, matchStr, mass2),
verticalalignment='bottom', horizontalalignment='right',
fontsize=40.0, transform=ax1.transAxes, weight='bold',
color='k', backgroundcolor='w')
# --------------------------------------------------------------------------------------- #
## Secondary Axis
ax2 = ax1.twiny()
kpcs = np.asarray(kpcArr)
kpcTicks= (kpcs ** 0.25)
ax2.set_xlim(xmin, xmax)
ax2.set_xticks(kpcTicks)
ax2.set_xticklabels(['{:g}'.format(kpc) for kpc in kpcs], fontsize=30)
for tick in ax2.xaxis.get_major_ticks():
tick.label.set_fontsize(30)
for tick in ax2.yaxis.get_major_ticks():
tick.label.set_fontsize(30)
ax2.text(0.92, 1.0035, '$\mathrm{Kpc}$',
verticalalignment='bottom', horizontalalignment='left',
fontsize=32.0, transform=ax2.transAxes)
# --------------------------------------------------------------------------------------- #
## Highlight zero
ax3.axhline(0.0, linewidth=4.0, c='k', linestyle='-', zorder=0, alpha=0.3)
## Mark the two interesting radius
if highlight1:
ax3.axvline(vline1 ** 0.25, linewidth=5.5, c='k', linestyle='--',
zorder=0, alpha=0.5, dashes=(30, 6))
else:
ax3.axvline(vline1 ** 0.25, linewidth=4.0, c='k', linestyle='--',
zorder=0, alpha=0.2)
if highlight2:
ax3.axvline(vline2 ** 0.25, linewidth=5.5, c='k', linestyle='--',
zorder=0, alpha=0.5, dashes=(30, 6))
else:
ax3.axvline(vline2 ** 0.25, linewidth=4.0, c='k', linestyle='--',
zorder=0, alpha=0.2)
# --------------------------------------------------------------------------------------- #
## Random Mixed Sample
ax3.fill_between(RSMA_COMMON,
np.percentile(randM_sm - gu_mm[2], 3, axis=0),
np.percentile(randM_sm - gu_mm[2], 97, axis=0),
facecolor='k', edgecolor='k', alpha=0.05, zorder=0)
ax3.fill_between(RSMA_COMMON,
np.percentile(randM_sm - gu_mm[2], 31, axis=0),
np.percentile(randM_sm - gu_mm[2], 69, axis=0),
facecolor='k', edgecolor='k', alpha=0.30, zorder=0)
ax3.plot(RSMA_COMMON, np.percentile(randM_sm - gu_mm[2], 50, axis=0),
c='k', linewidth=4.5, linestyle='--', alpha=0.9)
# --------------------------------------------------------------------------------------- #
## Difference between sample1 and sample2
diffLabel = '$[$' + label2 + '$]-[$' + label1 + '$]$'
ax3.fill_between(RSMA_COMMON,
(bm_mm[0] - gu_mm[1]), (bm_mm[1] - gu_mm[0]),
facecolor=diffColor1, edgecolor='k', alpha=0.50, zorder=1)
ax3.plot(RSMA_COMMON, bm_mm[2] - gu_mm[2], c=diffColor2, linewidth=5.5,
linestyle='-', label=diffLabel)
if showLMask:
ax3.plot(RSMA_COMMON, bm_mm_b[2] - gu_mm_b[2], c=diffColor2, linewidth=4.5,
linestyle='--')
if showAvg:
ax3.plot(RSMA_COMMON, bm_amg[2] - gu_amg[2], c=diffColor2, linewidth=4.0,
linestyle='-.')
# --------------------------------------------------------------------------------------- #
## X, Y- Label
ax3.set_xlabel('$R^{1/4}\ (\mathrm{Kpc})$', size=40)
ax3.set_ylabel('$\Delta$', size=32)
# --------------------------------------------------------------------------------------- #
## X, Y Limits
ax3.set_xlim(xmin, xmax)
ax3.set_ylim(dmin, dmax)
# --------------------------------------------------------------------------------------- #
## Region affected by PSF
# z = 0.4 : 1"=5.4 Kpc
# z = 0.5 : 1"=6.1 Kpc
dSep = (dmax - dmin) / 8.0
ax3.fill_between([0.0, rPsfKpc ** 0.25],
[dmin - dSep, dmin - dSep], [dmax + dSep, dmax + dSep],
facecolor='k', edgecolor='k', alpha=0.05, zorder=0)
# --------------------------------------------------------------------------------------- #
## Legend
if showLegend:
ax3.legend(loc=(0.255, 0.86), shadow=True, fancybox=True,
numpoints=1, fontsize=20, scatterpoints=1,
markerscale=1.2, borderpad=0.3, handletextpad=0.2)
# --------------------------------------------------------------------------------------- #
## Axes setup
ax3.minorticks_on()
ax3.tick_params(axis='y', which='minor', left='off', right='off')
# --------------------------------------------------------------------------------------- #
plt.show()
if save:
fig.savefig(outPng + '_b.png', dpi=300)
return fig
def plotKcorrection(cat, alphaUse=0.6, sizeUse=60,
ylim0=-0.5, ylim1=3.0, save=True,
outDir='./', catStr='hsc_massive'):
"""Make a summary: redshift-kcorrection figure."""
fig = plt.figure(figsize=(12, 12))
ax1 = plt.axes(REC)
ax1 = songPlotSetup(ax1)
# ---------------------------------------------------------------------------
# Mark a few redshift
ax1.axvline(0.2, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)
ax1.axvline(0.4, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)
ax1.axvline(0.5, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)
# Scatter plots
ax1.scatter(cat['Z'], cat['KCORRECT_G'] + 0.1, c='b',
marker='o', facecolor='none', edgecolor='b',
s=sizeUse, alpha=alphaUse, label='$\mathrm{HSC-}g$')
ax1.scatter(cat['Z'], cat['KCORRECT_R'] + 0.1, c='c',
marker='o', facecolor='none', edgecolor='c',
s=sizeUse, alpha=alphaUse, label='$\mathrm{HSC-}r$')
ax1.scatter(cat['Z'], cat['KCORRECT_I'] + 0.1, c='g',
marker='o', facecolor='none', edgecolor='g',
s=sizeUse, alpha=alphaUse, label='$\mathrm{HSC-}i$')
ax1.scatter(cat['Z'], cat['KCORRECT_Z'] + 0.1, c='orange',
marker='o', facecolor='none', edgecolor='orange',
s=sizeUse, alpha=alphaUse, label='$\mathrm{HSC-}z$')
ax1.scatter(cat['Z'], cat['KCORRECT_Y'] + 0.1, c='r',
marker='o', facecolor='none', edgecolor='r',
s=sizeUse, alpha=alphaUse, label='$\mathrm{HSC-Y}$')
# Label
ax1.set_xlabel('$\mathrm{Redshift}$', size=40)
ax1.set_ylabel('$\mathrm{K-correction}\ (\mathrm{mag})$', size=40)
# Axis limits
xmin, xmax = np.nanmin(cat['Z']), np.nanmax(cat['Z'])
xsep = (xmax - xmin) / 8.0
kcorAll = np.asarray([cat['KCORRECT_G'], cat['KCORRECT_R'], cat['KCORRECT_I'],
cat['KCORRECT_Z'], cat['KCORRECT_Y']])
ymin, ymax = np.nanmin(kcorAll), np.nanmax(kcorAll)
ymin = ymin if ymin > ylim0 else ylim0
ymax = ymax if ymax < ylim1 else ylim1
ysep = (ymax - ymin) / 8.0
ax1.set_xlim(xmin - xsep, xmax + xsep)
ax1.set_ylim(ymin - ysep / 2.0, ymax + ysep)
# Legend
ax1.legend(loc=(0.023, 0.75), shadow=True, fancybox=True,
numpoints=1, fontsize=22, scatterpoints=1,
markerscale=1.2, borderpad=0.3, handletextpad=0.2)
legend = ax1.get_legend()
legend.legendHandles[0].set_alpha(1.0)
legend.legendHandles[1].set_alpha(1.0)
legend.legendHandles[2].set_alpha(1.0)
legend.legendHandles[3].set_alpha(1.0)
legend.legendHandles[4].set_alpha(1.0)
## Running median
XX = cat['z_use']
medBins = np.linspace(XX.min(), XX.max(), 20)
dltBins = (medBins[1] - medBins[0])
indBins = np.digitize(XX, medBins)
YY1 = (cat['KCORRECT_G'] + 0.1)
medRunG = [np.nanmedian(YY1[indBins == kk]) for kk in range(20)]
YY2 = (cat['KCORRECT_R'] + 0.1)
medRunR = [np.nanmedian(YY2[indBins == kk]) for kk in range(20)]
YY3 = (cat['KCORRECT_I'] + 0.1)
medRunI = [np.nanmedian(YY3[indBins == kk]) for kk in range(20)]
YY4 = (cat['KCORRECT_Z'] + 0.1)
medRunZ = [np.nanmedian(YY4[indBins == kk]) for kk in range(20)]
YY5 = (cat['KCORRECT_Y'] + 0.1)
medRunY = [np.nanmedian(YY5[indBins == kk]) for kk in range(20)]
ax1.plot((medBins - 0.5 * dltBins), medRunG, c='b', linestyle='--', linewidth=7.0)
ax1.plot((medBins - 0.5 * dltBins), medRunR, c='c', linestyle='--', linewidth=7.0)
ax1.plot((medBins - 0.5 * dltBins), medRunI, c='g', linestyle='--', linewidth=7.0)
ax1.plot((medBins - 0.5 * dltBins), medRunZ, c='orange', linestyle='--', linewidth=7.0)
ax1.plot((medBins - 0.5 * dltBins), medRunY, c='r', linestyle='--', linewidth=7.0)
plt.show()
if save:
fig.savefig(os.path.join(outDir, catStr + '_z_kcorrect.png'), dpi=200)
In [3]:
# Location of the data
newDir = '/Users/songhuang/work/hscs/gama_massive/sbp/'
# Location for figures
figDir = 'figure'
# Location for subsamples
sampleDir = 'sample'
# The SED models
# 'fsps1', 'fsps2', 'fsps3', 'fsps4', 'bc03a'
sedMod = 'fsps2'
# Catalog files for BCG and NonBCG
redbcgStr = 'redbcg_old_s15b_sed_' + sedMod
nonbcgStr = 'nonbcg_old_s15b_sed_' + sedMod
redbcgFile = redbcgStr + '.fits'
nonbcgFile = nonbcgStr + '.fits'
# Output
redbcgPrep = redbcgFile.replace('.fits', '_prep.fits')
nonbcgPrep = nonbcgFile.replace('.fits', '_prep.fits')
try:
redbcgTab
except NameError:
pass
else:
del redbcgTab
try:
nonbcgTab
except NameError:
pass
else:
del nonbcgTab
# Location for the SBP summary file
redbcgDir = os.path.join(newDir, 'redbcg')
nonbcgDir = os.path.join(newDir, 'gama')
# Two summary catalogs
redbcgCat = os.path.join(newDir, redbcgFile)
nonbcgCat = os.path.join(newDir, nonbcgFile)
prefix1 = 'redbcg'
prefix2 = 'nonbcg'
# Find and read in the catalogs
if not os.path.isfile(redbcgCat):
raise Exception("## Can not find catalog for redBCG galaxies: %s" % redbcgCat)
else:
redbcgTab = Table.read(redbcgCat, format='fits')
if not os.path.isfile(nonbcgCat):
raise Exception("## Can not find catalog for nonBCG galaxies: %s" % nonbcgCat)
else:
nonbcgTab = Table.read(nonbcgCat, format='fits')
print("## Deal with %i galaxies in redBCG sample" % len(redbcgTab))
print("## Deal with %i galaxies in nonBCG sample" % len(nonbcgTab))
# For redBCG sample
redbcgTab = updateMass(redbcgTab, m2l='LOGM2L_I_OBS')
redbcgTab = updateKcorrect(redbcgTab, zCol='z_use', magType='mag_cmodel')
# For nonBCG sample
nonbcgTab = updateMass(nonbcgTab, m2l='LOGM2L_I_OBS')
nonbcgTab = updateKcorrect(nonbcgTab, zCol='z_use', magType='mag_cmodel')
# Save the results
redbcgTab.write(os.path.join(newDir, redbcgPrep), format='fits',
overwrite=True)
nonbcgTab.write(os.path.join(newDir, nonbcgPrep), format='fits',
overwrite=True)
In [4]:
plotKcorrection(redbcgTab, alphaUse=0.40, sizeUse=60,
outDir=figDir, catStr=redbcgStr)
plotKcorrection(nonbcgTab, alphaUse=0.15, sizeUse=30,
outDir=figDir, catStr=nonbcgStr)
In [23]:
#-----------------------------------------------------------------#
fig = plt.figure(figsize=(16, 16))
fig.subplots_adjust(left=0.10, right=0.90,
bottom=0.08, top=0.99,
wspace=0.02, hspace=0.00)
# logM - (g-r) Color
ax1 = fig.add_subplot(221)
ax1 = songPlotSetup(ax1)
ax1.scatter(nonbcgTab['logm_100'], (nonbcgTab['ABSMAG_G'] - nonbcgTab['ABSMAG_R']),
marker='o', edgecolor='none', cmap=BLU5, s=40,
facecolor=toColorArr(nonbcgTab['z_use'], bottom=0.1), alpha=0.2)
ax1.scatter(redbcgTab['logm_100'], (redbcgTab['ABSMAG_G'] - redbcgTab['ABSMAG_R']),
marker='o', edgecolor='none', cmap=ORG4, s=150,
facecolor=toColorArr(redbcgTab['z_use'], bottom=0.1), alpha=0.8)
ax1.tick_params(axis='x', which='both', labelbottom='off')
ax1.set_ylabel('$(g-r)\ \mathrm{K-corrected}$', size=40)
ax1.set_xlim(10.6, 12.3)
ax1.set_ylim(0.15, 1.49)
# logM - (g-z) Color
ax2 = fig.add_subplot(223)
ax2 = songPlotSetup(ax2)
ax2.scatter(nonbcgTab['logm_100'], (nonbcgTab['ABSMAG_G'] - nonbcgTab['ABSMAG_Z']),
marker='o', edgecolor='none', cmap=BLU5, s=40,
facecolor=toColorArr(nonbcgTab['z_use'], bottom=0.1), alpha=0.2)
ax2.scatter(redbcgTab['logm_100'], (redbcgTab['ABSMAG_G'] - redbcgTab['ABSMAG_Z']),
marker='o', edgecolor='none', cmap=ORG4, s=150,
facecolor=toColorArr(redbcgTab['z_use'], bottom=0.1), alpha=0.8)
ax2.set_xlabel('$\log\ (M_{\star}/M_{\odot})$', size=40)
ax2.set_ylabel('$(g-z)\ \mathrm{K-corrected}$', size=40)
ax2.set_xlim(10.6, 12.3)
ax2.set_ylim(0.55, 2.69)
# logM - (r-z) Color
ax3 = fig.add_subplot(222)
ax3 = songPlotSetup(ax3)
ax3.scatter(nonbcgTab['logm_100'], (nonbcgTab['ABSMAG_R'] - nonbcgTab['ABSMAG_Z']),
marker='o', edgecolor='none', cmap=BLU5, s=40,
facecolor=toColorArr(nonbcgTab['z_use'], bottom=0.1), alpha=0.2)
ax3.scatter(redbcgTab['logm_100'], (redbcgTab['ABSMAG_R'] - redbcgTab['ABSMAG_Z']),
marker='o', edgecolor='none', cmap=ORG4, s=150,
facecolor=toColorArr(redbcgTab['z_use'], bottom=0.1), alpha=0.8)
ax3.tick_params(axis='y', which='both', labelsize=30,
labelleft='off', labelright='on')
ax3.yaxis.set_label_position("right")
ax3.set_ylabel('$(r-z)\ \mathrm{K-corrected}$', size=40)
ax3.set_xlim(10.6, 12.3)
ax3.set_ylim(0.25, 1.39)
# logM - (i-y) Color
ax4 = fig.add_subplot(224)
ax4 = songPlotSetup(ax4)
ax4.scatter(nonbcgTab['logm_100'], (nonbcgTab['ABSMAG_I'] - nonbcgTab['ABSMAG_Y']),
marker='o', edgecolor='none', cmap=BLU5, s=40,
facecolor=toColorArr(nonbcgTab['z_use'], bottom=0.1), alpha=0.2)
ax4.scatter(redbcgTab['logm_100'], (redbcgTab['ABSMAG_I'] - redbcgTab['ABSMAG_Y']),
marker='o', edgecolor='none', cmap=ORG4, s=150,
facecolor=toColorArr(redbcgTab['z_use'], bottom=0.1), alpha=0.8)
ax4.tick_params(axis='y', which='both', labelsize=30,
labelleft='off', labelright='on')
ax4.yaxis.set_label_position("right")
ax4.set_xlabel('$\log\ (M_{\star}/M_{\odot})$', size=40)
ax4.set_ylabel('$(i-y)\ \mathrm{K-corrected}$', size=40)
ax4.set_xlim(10.6, 12.3)
ax4.set_ylim(0.05, 0.89)
fig.savefig(os.path.join(figDir, redbcgStr + '_mass_color.png'), dpi=200)