In [2]:
%pylab 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
# Matplotlib related
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from matplotlib.ticker import NullFormatter
from matplotlib.ticker import MaxNLocator
# 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
# Personal
import hscUtils as hUtil
import galSBP
import coaddCutoutGalfitSimple as gSimple
import matplotlib.patches as mpatches
from matplotlib.patches import Ellipse
from matplotlib.collections import PatchCollection
import cosmology
c=cosmology.Cosmo(H0=70.0, omega_m=0.3, omega_l=0.7, flat=1)
In [3]:
# Absolute magnitude of sun in HSC filters
# Actuall borrowed from DES filters
# Values from magsun.data in FSPS
amag_sun_des_g = 5.08
amag_sun_des_r = 4.62
amag_sun_des_i = 4.52
amag_sun_des_z = 4.52
amag_sun_des_y = 4.51
amag_sun_sdss_g = 5.45
amag_sun_sdss_r = 4.67
amag_sun_sdss_i = 4.56
# Based on http://www.baryons.org/ezgal/filters.php
amag_sun_ukiss_y = 4.515
# Extinction correction factor for HSC
## A\_lambda = Coeff * E(B-V)
a_hsc_g = 3.233
a_hsc_r = 2.291
a_hsc_i = 1.635
a_hsc_z = 1.261
a_hsc_y = 1.076
#
SIGMA1 = 0.3173
SIGMA2 = 0.0455
SIGMA3 = 0.0027
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'
In [4]:
# Code for Get Bootstrap mean or median
def _confidence_interval_1d(A, alpha=.05, metric=np.mean, numResamples=10000, interpolate=True):
"""Calculates bootstrap confidence interval along one dimensional array"""
if not isinstance(alpha, collections.Iterable):
alpha = np.array([alpha])
N = len(A)
resampleInds = np.random.randint(0, N, (numResamples,N))
metricOfResampled = metric(A[resampleInds], axis=-1)
confidenceInterval = np.zeros(2*len(alpha),dtype='float')
if interpolate:
for thisAlphaInd, thisAlpha in enumerate(alpha):
confidenceInterval[2*thisAlphaInd] = scoreatpercentile(metricOfResampled,
thisAlpha*100/2.0)
confidenceInterval[2*thisAlphaInd+1] = scoreatpercentile(metricOfResampled,
100-thisAlpha*100/2.0)
else:
sortedMetricOfResampled = np.sort(metricOfResampled)
for thisAlphaInd, thisAlpha in enumerate(alpha):
confidenceInterval[2*thisAlphaInd] = sortedMetricOfResampled[int(round(thisAlpha*numResamples/2.0))]
confidenceInterval[2*thisAlphaInd+1] = sortedMetricOfResampled[int(round(numResamples -
thisAlpha*numResamples/2.0))]
return confidenceInterval
def _ma_confidence_interval_1d(A, alpha=.05, metric=np.mean, numResamples=10000, interpolate=True):
A = np.ma.masked_invalid(A, copy=True)
A = A.compressed()
confidenceInterval = _confidence_interval_1d(A, alpha, metric, numResamples, interpolate)
return confidenceInterval
def confidence_interval(A, axis=None, alpha=.05, metric=np.mean, numResamples=10000, interpolate=True):
"""Return the bootstrap confidence interval of an array or along an axis ignoring NaNs and masked elements.
Parameters
----------
A : array_like
Array containing numbers whose confidence interval is desired.
axis : int, optional
Axis along which the confidence interval is computed.
The default is to compute the confidence interval of the flattened array.
alpha: float or array, optional
confidence level of confidence interval. 100.0*(1-alpha) percent confidence
interval will be returned.
If length-n array, n confidence intervals will be computed
The default is .05
metric : numpy function, optional
metric to calculate confidence interval for.
The default is numpy.mean
numResamples : int, optional
number of bootstrap samples. The default is 10000.
interpolate: bool, optional
uses scipy.stats.scoreatpercentile to interpolate between bootstrap samples
if alpha*numResamples/2.0 is not integer.
The default is True
Returns
-------
confidenceInterval : ndarray
An array with the same shape as `A`, with the specified axis replaced by one twice the length of the alpha
If `A` is a 0-d array, or if axis is None, a length-2 ndarray is returned.
"""
if interpolate is True and scoreatpercentile is False:
print("need scipy to interpolate between values")
interpolate = False
A = A.copy()
if axis is None:
A = A.ravel()
outA = _ma_confidence_interval_1d(A, alpha, metric, numResamples, interpolate)
else:
outA = np.apply_along_axis(_ma_confidence_interval_1d, axis, A, alpha,
metric, numResamples, interpolate)
return outA
def normProf(sma, sbp, minSma, maxSma, divide=False):
"""
Naive method to normalize the profile.
Parameters:
sbp : Array for surface brightness profile
sma : Radius range
minSma : Minimum SMA
maxSma Maximum SMA
"""
offset = np.nanmedian(sbp[(sma >= minSma) &
(sma <= maxSma)])
if divide:
return (sbp / offset)
else:
return (sbp-offset)
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_b_I'] = g['KCORRECT_b_I']
gProf.meta['KCORRECT_c_I'] = g['KCORRECT_c_I']
gProf.meta['KCORRECT_G'] = g['KCORRECT_G']
gProf.meta['KCORRECT_b_G'] = g['KCORRECT_b_G']
gProf.meta['KCORRECT_c_G'] = g['KCORRECT_c_G']
gProf.meta['KCORRECT_R'] = g['KCORRECT_R']
gProf.meta['KCORRECT_b_R'] = g['KCORRECT_b_R']
gProf.meta['KCORRECT_c_R'] = g['KCORRECT_c_R']
gProf.meta['KCORRECT_Z'] = g['KCORRECT_Z']
gProf.meta['KCORRECT_b_Z'] = g['KCORRECT_b_Z']
gProf.meta['KCORRECT_c_Z'] = g['KCORRECT_c_Z']
gProf.meta['KCORRECT_Y'] = g['KCORRECT_Y']
gProf.meta['KCORRECT_b_Y'] = g['KCORRECT_b_Y']
gProf.meta['KCORRECT_c_Y'] = g['KCORRECT_c_Y']
gProf.meta['LOGM2LI_A'] = g['logm2lI_A']
gProf.meta['LOGM2LI_B'] = g['logm2lI_B']
gProf.meta['LOGM2LI_C'] = g['logm2lI_C']
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_c_I',
kind='sbp', norm=False, r1=9.9, r2=10.1, divide=False,
col3=None, col4=None, justStack=False,
sun1=amag_sun_des_g, sun2=amag_sun_des_r):
""" Get the stack of individual profiels, and their med/avg. """
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)
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)
for p in profiles)
else:
stack = np.vstack(np.asarray(p[col1]) for p in profiles)
elif kind.strip() == 'mass':
if norm:
stack = np.vstack(normProf(p['rKpc'],
np.asarray(p[col1] + p.meta[col2]),
r1, r2, divide=divide) for p in profiles)
else:
stack = np.vstack(np.asarray(p[col1] + p.meta[col2]) for p in 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) 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) 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)
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 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
In [5]:
newDir = '/Users/songhuang/work/hscs/gama_massive/sbp/'
bcgFile = 'redbcg_1d_160211.fits'
memFile = 'redmem_1d_160211.fits'
gamaFile = 'gama_1d_160211.fits'
try:
bcgTab
except NameError:
pass
else:
del bcgTab
try:
memTab
except NameError:
pass
else:
del memTab
try:
gamaTab
except NameError:
pass
else:
del gamaTab
# Two summary catalogs
bcgCat = os.path.join(newDir, bcgFile)
memCat = os.path.join(newDir, memFile)
gamaCat = os.path.join(newDir, gamaFile)
if not os.path.isfile(bcgCat):
raise Exception("## Can not find catalog for BCGs : %s" % bcgCat)
else:
bcgTab = Table.read(bcgCat, format='fits')
if not os.path.isfile(memCat):
raise Exception("## Can not find catalog for cluster members : %s" % memCat)
else:
memTab = Table.read(memCat, format='fits')
if not os.path.isfile(gamaCat):
raise Exception("## Can not find catalog for GAMA galaxies : %s" % gamaCat)
else:
gamaTab = Table.read(gamaCat, format='fits')
print("## Deal with %i galaxies in redBCH sample" % len(bcgTab))
print("## Deal with %i galaxies in redMEM sample" % len(memTab))
print("## Deal with %i galaxies in GAMA sample" % len(gamaTab))
In [6]:
def doubleSchechter(logm, logm0=10.91,
logphi1=-2.97, logphi2=-2.79,
alpha1=-0.46, alpha2=-1.58):
phi1 = (10.0 ** logphi1)
phi2 = (10.0 ** logphi2)
dlogm = (logm - logm0)
term1 = np.log(10.0) * np.exp(-1.0 * (10.0 ** dlogm))
term2 = phi1 * (10.0 ** ((alpha1 + 1.0) * dlogm))
term3 = phi2 * (10.0 ** ((alpha2 + 1.0) * dlogm))
return term1 * (term2 + term3)
In [7]:
def bernardiLF(logL, model):
L = (10.0 ** logL)
if model is 'cmodel':
phi1 = 0.928E-2
L1 = 0.3077E9
alpha = 1.918
beta = 0.433
phi2 = 0.964E-2
L2 = 1.8763E9
gamma = 0.479
elif model is 'sersic':
phi1 = 1.343E-2
L1 = 0.0187E9
alpha = 1.678
beta = 0.300
phi2 = 0.843E-2
L2 = 0.8722E9
gamma = 1.058
elif model is 'serexp':
phi1 = 1.348E-2
L1 = 0.3223E9
alpha = 1.297
beta = 0.398
phi2 = 0.820E-2
L2 = 0.9081E9
gamma = 1.131
elif model is 'serg2d':
phi1 = 1.902E-2
L1 = 6.2456E9
alpha = 0.497
beta = 0.589
phi2 = 0.530E-2
L2 = 0.8263E9
gamma = 1.260
else:
raise Exception("# Wrong model!")
term1 = phi1 * beta * ((L / L1) ** alpha)
term2 = (np.exp(-1.0 * ((L / L1) ** beta))) / (scipy.special.gamma(alpha / beta))
term3 = phi2 * ((L / L2) ** gamma) * (np.exp(-1.0 * L / L2))
return term1 * term2 + term3
In [10]:
lfB12 = Table.read('iband-lf-B12.txt', format='ascii')
In [8]:
area1 = 65.0
area2 = 90.0
area3 = 92.0
# 0.15 < z < 0.40
gamaZ1 = gamaTab[(gamaTab['Z'] >= 0.15) & (gamaTab['Z'] <= 0.4)]
memZ1 = memTab[(memTab['Z'] >= 0.15) & (memTab['Z'] <= 0.4)]
bcgZ1 = bcgTab[(bcgTab['Z'] >= 0.15) & (bcgTab['Z'] <= 0.4)]
dGama1 = (c.V(np.nanmin(gamaZ1['Z']), np.nanmax(gamaZ1['Z'])) / ((360.0 ** 2.0) / np.pi)) * area1
dMem1 = (c.V(np.nanmin(memZ1['Z']), np.nanmax(memZ1['Z'])) / ((360.0 ** 2.0) / np.pi)) * area3
dBcg1 = (c.V(np.nanmin(bcgZ1['Z']), np.nanmax(bcgZ1['Z'])) / ((360.0 ** 2.0) / np.pi)) * area2
# 0.15 < z < 0.50
gamaZ2 = gamaTab[(gamaTab['Z'] >= 0.15) & (gamaTab['Z'] <= 0.5)]
memZ2 = memTab[(memTab['Z'] >= 0.15) & (memTab['Z'] <= 0.5)]
bcgZ2 = bcgTab[(bcgTab['Z'] >= 0.15) & (bcgTab['Z'] <= 0.5)]
dGama2 = (c.V(np.nanmin(gamaZ2['Z']), np.nanmax(gamaZ2['Z'])) / ((360.0 ** 2.0) / np.pi)) * area1
dMem2 = (c.V(np.nanmin(memZ2['Z']), np.nanmax(memZ2['Z'])) / ((360.0 ** 2.0) / np.pi)) * area3
dBcg2 = (c.V(np.nanmin(bcgZ2['Z']), np.nanmax(bcgZ2['Z'])) / ((360.0 ** 2.0) / np.pi)) * area2
print(dGama1, dBcg1, dMem1)
print(dGama2, dBcg2, dMem2)
In [21]:
fig = plt.figure(figsize=(14, 8))
rec = [0.11, 0.12, 0.87, 0.87]
ax1 = plt.axes(rec)
ax1.axvline(10.9, linewidth=8.0, color='k', linestyle='--',
alpha=0.2, zorder=0)
# ---------------------------------------------------------------------------
# Histogram
## cModel
temp = ax1.hist(gamaZ1['lumI_cmodel'], bins=24, range=[10.1, 11.9],
orientation='vertical', histtype='stepfilled',
color='k', alpha=0.2, normed=1,
label='$\Lambda \leq 20\mathrm{\ cModel},\ z=0.40$')
temp = ax1.hist(bcgZ1['lumI_cmodel'], bins=24, range=[10.1, 11.9],
orientation='vertical', histtype='stepfilled',
color='r', alpha=0.2, normed=1,
label='$\Lambda > 20\mathrm{\ cModel},\ z=0.40$')
## SBP 100
temp = ax1.hist(gamaZ1['lum_100'], bins=24, range=[10.1, 11.9], linewidth=5.0,
orientation='vertical', histtype='step',
color='k', alpha=0.9, normed=1,
label='$\Lambda \leq 20\mathrm{\ 100 Kpc},\ z=0.40$')
temp = ax1.hist(bcgZ1['lum_100'], bins=24, range=[10.1, 11.9], linewidth=5.0,
orientation='vertical', histtype='step',
color='r', alpha=0.9, normed=1,
label='$\Lambda > 20\mathrm{\ 100 Kpc},\ z=0.40$')
## SBP 100
temp = ax1.hist(gamaZ2['lum_100'], bins=24, range=[10.1, 11.9], linewidth=5.0,
orientation='vertical', histtype='step',
color='k', alpha=0.9, normed=1, linestyle=':',
label='$\Lambda \leq 20\mathrm{\ 100 Kpc},\ z=0.50$')
temp = ax1.hist(bcgZ2['lum_100'], bins=24, range=[10.1, 11.9], linewidth=5.0,
orientation='vertical', histtype='step',
color='r', alpha=0.9, normed=1, linestyle=':',
label='$\Lambda > 20\mathrm{\ 100 Kpc},\ z=0.50$')
## SBP 150
temp = ax1.hist(gamaZ1['lum_max'], bins=24, range=[10.1, 11.9], linewidth=5.0,
orientation='vertical', histtype='step',
color='k', alpha=0.9, normed=1, linestyle='--',
label='$\Lambda \leq 20\mathrm{\ 150 Kpc},\ z=0.40$')
temp = ax1.hist(bcgZ1['lum_max'], bins=24, range=[10.1, 11.9], linewidth=5.0,
orientation='vertical', histtype='step',
color='r', alpha=0.9, normed=1, linestyle='--',
label='$\Lambda > 20\mathrm{\ 150 Kpc},\ z=0.40$')
ax1.legend(loc=(0.02, 0.55), shadow=True, fancybox=True,
numpoints=1, fontsize=16, scatterpoints=1,
markerscale=1.2, borderpad=0.25, handletextpad=0.25)
ax1.set_xlim(10.21, 11.9)
ax1.set_ylim(0.02, 3.75)
# Axes setup
# Minor Ticks on
ax1.minorticks_on()
# Axes Thickness
for axis in ['top','bottom','left','right']:
ax1.spines[axis].set_linewidth(3.5)
# Tick Label Size
for tick in ax1.xaxis.get_major_ticks():
tick.label.set_fontsize(24)
for tick in ax1.yaxis.get_major_ticks():
tick.label.set_fontsize(24)
# Tick Length and Width
ax1.tick_params('both', length=10, width=3.0, which='major')
ax1.tick_params('both', length=6, width=2.5, which='minor')
# Label
ax1.set_xlabel('$\log\ (L_{\star}/L_{\odot})\ \mathrm{(HSC)}$', size=35)
ax1.set_ylabel('$\mathrm{Normalized\ \#}$', size=39)
fig.savefig('../figure/hscMassive_lum_distri_norm.png', dpi=300)
plt.show()
In [22]:
fig = plt.figure(figsize=(14, 8))
rec = [0.11, 0.12, 0.87, 0.87]
ax1 = plt.axes(rec)
ax1.axvline(10.9, linewidth=8.0, color='k', linestyle='--',
alpha=0.2, zorder=0)
# ---------------------------------------------------------------------------
# Histogram
## cModel
nGama1, bGama1, pGama1 = ax1.hist(gamaZ1['lumI_cmodel'],
bins=24, range=[10.1, 11.9],
orientation='vertical', histtype='stepfilled',
color='k', alpha=0.2, normed=0,
label='$\Lambda \leq 20\mathrm{\ cModel},\ z=0.4$')
nBcg1, bBcg1, pBcg1 = ax1.hist(bcgZ1['lumI_cmodel'],
bins=24, range=[10.1, 11.9],
orientation='vertical', histtype='stepfilled',
color='r', alpha=0.2, normed=0,
label='$\Lambda > 20\mathrm{\ cModel},\ z=0.4$')
## SBP 100
nGama2, bGama2, pGama2 = ax1.hist(gamaZ1['lum_100'],
bins=24, range=[10.1, 11.9], linewidth=5.0,
orientation='vertical', histtype='step',
color='k', alpha=0.9, normed=0,
label='$\Lambda \leq 20\mathrm{\ 100 Kpc},\ z=0.4$')
nBcg2, bBcg2, pBcg2 = ax1.hist(bcgZ1['lum_100'],
bins=24, range=[10.1, 11.9], linewidth=5.0,
orientation='vertical', histtype='step',
color='r', alpha=0.9, normed=0,
label='$\Lambda > 20\mathrm{\ 100 Kpc},\ z=0.4$')
## SBP 100; up to z=0.5
nGama2b, bGama2b, pGama2b = ax1.hist(gamaZ2['lum_100'],
bins=24, range=[10.1, 11.9], linewidth=5.0,
orientation='vertical', histtype='step',
color='k', alpha=0.9, normed=0, linestyle=':',
label='$\Lambda \leq 20\mathrm{\ 100 Kpc},\ z=0.5$')
nBcg2b, bBcg2b, pBcg2b = ax1.hist(bcgZ2['lum_100'],
bins=24, range=[10.1, 11.9], linewidth=5.0,
orientation='vertical', histtype='step',
color='r', alpha=0.9, normed=0, linestyle=':',
label='$\Lambda > 20\mathrm{\ 100 Kpc},\ z=0.5$')
## SBP 150
nGama3, bGama3, pGama3 = ax1.hist(gamaZ1['lum_150'],
bins=24, range=[10.1, 11.9], linewidth=5.0,
orientation='vertical', histtype='step',
color='k', alpha=0.9, normed=0, linestyle='--',
label='$\Lambda \leq 20\mathrm{\ 150 Kpc},\ z=0.4$')
nBcg3, bBcg3, pBcg3 = ax1.hist(bcgZ1['lum_150'],
bins=24, range=[10.1, 11.9], linewidth=5.0,
orientation='vertical', histtype='step',
color='r', alpha=0.9, normed=0, linestyle='--',
label='$\Lambda > 20\mathrm{\ 150 Kpc},\ z=0.4$')
ax1.legend(loc=(0.64, 0.42), shadow=True, fancybox=True,
numpoints=1, fontsize=22, scatterpoints=1,
markerscale=1.2, borderpad=0.25, handletextpad=0.25)
ax1.set_xlim(10.21, 11.9)
ax1.set_ylim(1, 1660)
# Axes setup
# Minor Ticks on
ax1.minorticks_on()
# Axes Thickness
for axis in ['top','bottom','left','right']:
ax1.spines[axis].set_linewidth(3.5)
# Tick Label Size
for tick in ax1.xaxis.get_major_ticks():
tick.label.set_fontsize(24)
for tick in ax1.yaxis.get_major_ticks():
tick.label.set_fontsize(24)
# Tick Length and Width
ax1.tick_params('both', length=10, width=3.0, which='major')
ax1.tick_params('both', length=6, width=2.5, which='minor')
# Label
ax1.set_xlabel('$\log\ (L_{\star}/L_{\odot})\ \mathrm{(HSC)}$', size=35)
ax1.set_ylabel('$\mathrm{Number}$', size=39)
fig.savefig('../figure/hscMassive_lum_distri.png', dpi=300)
plt.show()
In [23]:
binSize = bGama1[1] - bGama1[0]
print(binSize)
dnGama1 = np.log10(nGama1 / (dGama1 * binSize))
dnBcg1 = np.log10(nBcg1 / (dBcg1 * binSize))
dnGama2 = np.log10(nGama2 / (dGama1 * binSize))
dnBcg2 = np.log10(nBcg2 / (dBcg1 * binSize))
dnGama3 = np.log10(nGama3 / (dGama1 * binSize))
dnBcg3 = np.log10(nBcg3 / (dBcg1 * binSize))
dnGama2b = np.log10(nGama2b / (dGama2 * binSize))
dnBcg2b = np.log10(nBcg2b / (dBcg2 * binSize))
In [34]:
fig = plt.figure(figsize=(13, 11))
rec = [0.12, 0.12, 0.87, 0.87]
ax1 = plt.axes(rec)
#ax1.axvline(10.9, color='k', linewidth=8.0, alpha=0.4)
# Baldry+2012
ax1.plot(lfB12['col1'] - 0.04, np.log10(lfB12['col3']), c='k', linestyle='-',
linewidth=10.0, alpha=0.3, label='$\mathrm{Baldry\ et\ al.\ (2012)}$')
# GAMA
ax1.plot((bGama2[1:] + bGama2[0:-1])/2.0, dnGama1, c=BLUE0,
linewidth=8.0, label='$\Lambda \leq 20\mathrm{\ cModel},\ z=0.40$',
linestyle='-', alpha=0.6)
"""
ax1.plot((bGama2[1:] + bGama2[0:-1])/2.0, dnGama2, c=BLUE1,
linewidth=5.0, label='$\Lambda \leq 20\mathrm{\ 100 Kpc},\ z=0.40$',
linestyle='-', alpha=1.0)
ax1.plot((bGama2b[1:] + bGama2b[0:-1])/2.0, dnGama2b, c=BLUE1,
linewidth=5.0, label='$\Lambda \leq 20\mathrm{\ 100 Kpc},\ z=0.50$',
linestyle='--', alpha=0.8)
ax1.plot((bGama3[1:] + bGama3[0:-1])/2.0, dnGama3, c=BLUE1,
linewidth=4.0, label='$\Lambda \leq 20\mathrm{\ 150 Kpc},\ z=0.40$',
linestyle='-.', alpha=0.8)
"""
# BCG
ax1.plot((bBcg1[1:] + bBcg1[0:-1])/2.0, dnBcg1, c=RED0,
linewidth=8.0, label='$\Lambda > 20\mathrm{\ cModel},\ z=0.40$',
linestyle='-', alpha=0.5)
"""
ax1.plot((bBcg2[1:] + bBcg2[0:-1])/2.0, dnBcg2, c=RED1,
linewidth=5.0, label='$\Lambda > 20\mathrm{\ 100 Kpc},\ z=0.40$',
linestyle='-', alpha=1.0)
ax1.plot((bBcg2b[1:] + bBcg2b[0:-1])/2.0, dnBcg2b, c=RED1,
linewidth=5.0, label='$\Lambda > 20\mathrm{\ 100 Kpc},\ z=0.50$',
linestyle='--', alpha=0.8)
ax1.plot((bBcg3[1:] + bBcg3[0:-1])/2.0, dnBcg3, c=RED1,
linewidth=4.0, label='$\Lambda > 20\mathrm{\ 150 Kpc},\ z=0.40$',
linestyle='-.', alpha=0.8)
"""
# Axes setup
# Minor Ticks on
ax1.minorticks_on()
# Axes Thickness
for axis in ['top','bottom','left','right']:
ax1.spines[axis].set_linewidth(3.5)
# Tick Label Size
for tick in ax1.xaxis.get_major_ticks():
tick.label.set_fontsize(24)
for tick in ax1.yaxis.get_major_ticks():
tick.label.set_fontsize(24)
# Tick Length and Width
ax1.tick_params('both', length=10, width=3.0, which='major')
ax1.tick_params('both', length=6, width=2.5, which='minor')
ax1.legend(loc=(0.65, 0.56), shadow=True, fancybox=True,
numpoints=1, fontsize=19, scatterpoints=1,
markerscale=1.2, borderpad=0.5, handletextpad=0.34)
ax1.set_xlim(10.25, 11.9)
ax1.set_ylim(-6.79, -1.01)
# Label
ax1.set_xlabel('$\log\ (L_{\star}/L_{\odot})\ (\mathrm{HSC}-i\ \mathrm{band})}$', size=40)
ax1.set_ylabel('$\mathrm{d}N/\mathrm{d}\ {\logL_{\star}}\ [{\mathrm{Mpc}^{-3}}{\mathrm{dex}^{-1}}]$',
size=45)
fig.savefig('../figure/hscMassive_lum_function_1.png', dpi=300)
plt.show()
In [9]:
mfM13 = Table.read('smf_Moustakas_13.txt', format='ascii')
mfB13 = Table.read('smf_Bernardi_13_serExp.txt', format='ascii')
mfD15 = Table.read('smf_dSouza_15.txt', format='ascii')
massArr = np.arange(7.0, 13.0, 0.1)
In [10]:
mfM13.colnames, mfB13.colnames, mfD15.colnames
Out[10]:
In [11]:
fig = plt.figure(figsize=(14, 8))
rec = [0.11, 0.12, 0.87, 0.87]
ax1 = plt.axes(rec)
ax1.axvline(11.4, linewidth=4.0, color='g', linestyle='-',
alpha=0.2, zorder=0)
# ---------------------------------------------------------------------------
# Histogram
## cModel
temp = ax1.hist(gamaZ1['lumI_cmodel'] + gamaZ1['logm2lI_C'],
bins=24, range=[10.6, 12.6],
orientation='vertical', histtype='stepfilled',
color='k', alpha=0.2, normed=1,
label='$\Lambda \leq 20\mathrm{\ cModel},\ z=0.40$')
temp = ax1.hist(bcgZ1['lumI_cmodel'] + bcgZ1['logm2lI_C'],
bins=24, range=[10.6, 12.6],
orientation='vertical', histtype='stepfilled',
color='r', alpha=0.2, normed=1,
label='$\Lambda > 20\mathrm{\ cModel},\ z=0.40$')
## SBP 100
temp = ax1.hist(gamaZ1['lum_100'] + gamaZ1['logm2lI_C'],
bins=24, range=[10.6, 12.6], linewidth=5.0,
orientation='vertical', histtype='step',
color='k', alpha=0.9, normed=1,
label='$\Lambda \leq 20\mathrm{\ 100 Kpc},\ z=0.40$')
temp = ax1.hist(bcgZ1['lum_100'] + bcgZ1['logm2lI_C'],
bins=24, range=[10.6, 12.6], linewidth=5.0,
orientation='vertical', histtype='step',
color='r', alpha=0.9, normed=1,
label='$\Lambda > 20\mathrm{\ 100 Kpc},\ z=0.40$')
## SBP 100
temp = ax1.hist(gamaZ2['lum_100'] + gamaZ2['logm2lI_C'],
bins=24, range=[10.6, 12.6], linewidth=5.0,
orientation='vertical', histtype='step',
color='k', alpha=0.9, normed=1, linestyle=':',
label='$\Lambda \leq 20\mathrm{\ 100 Kpc},\ z=0.55$')
temp = ax1.hist(bcgZ2['lum_100'] + bcgZ2['logm2lI_C'],
bins=24, range=[10.6, 12.6], linewidth=5.0,
orientation='vertical', histtype='step',
color='r', alpha=0.9, normed=1, linestyle=':',
label='$\Lambda > 20\mathrm{\ 100 Kpc},\ z=0.55$')
## GAMA
temp = ax1.hist(gamaZ1['logms_gama'] + np.log10(gamaZ1['fluxscale_gama']),
bins=24, range=[10.6, 12.6], linewidth=5.0,
orientation='vertical', histtype='step',
color='k', alpha=0.9, normed=1, linestyle='--',
label='$\Lambda \leq 20\mathrm{\ GAMA},\ z=0.40$')
temp = ax1.hist(bcgZ1['logms_gama'] + np.log10(bcgZ1['fluxscale_gama']),
bins=24, range=[10.6, 12.6], linewidth=5.0,
orientation='vertical', histtype='step',
color='r', alpha=0.9, normed=1, linestyle='--',
label='$\Lambda > 20\mathrm{\ GAMA},\ z=0.40$')
ax1.legend(loc=(0.73, 0.55), shadow=True, fancybox=True,
numpoints=1, fontsize=16, scatterpoints=1,
markerscale=1.2, borderpad=0.25, handletextpad=0.25)
ax1.set_xlim(10.55, 12.59)
ax1.set_ylim(0.01, 2.99)
# Axes setup
# Minor Ticks on
ax1.minorticks_on()
# Axes Thickness
for axis in ['top','bottom','left','right']:
ax1.spines[axis].set_linewidth(3.5)
# Tick Label Size
for tick in ax1.xaxis.get_major_ticks():
tick.label.set_fontsize(24)
for tick in ax1.yaxis.get_major_ticks():
tick.label.set_fontsize(24)
# Tick Length and Width
ax1.tick_params('both', length=10, width=3.0, which='major')
ax1.tick_params('both', length=6, width=2.5, which='minor')
# Label
ax1.set_xlabel('$\log\ (M_{\star}/M_{\odot})\ \mathrm{(HSC)}$', size=35)
ax1.set_ylabel('$\mathrm{Normalized\ \#}$', size=39)
fig.savefig('../figure/hscMassive_mass_distri_norm.png', dpi=300)
plt.show()
In [12]:
fig = plt.figure(figsize=(14, 8))
rec = [0.11, 0.12, 0.87, 0.87]
ax1 = plt.axes(rec)
ax1.axvline(11.4, linewidth=4.0, color='g', linestyle='-',
alpha=0.2, zorder=0)
# ---------------------------------------------------------------------------
# Histogram
## SBP 100
temp = ax1.hist(gamaZ1['lum_100'] + gamaZ1['logm2lI_A'],
bins=24, range=[10.6, 12.6], linewidth=5.0,
orientation='vertical', histtype='step',
color='k', alpha=0.9, normed=1, linestyle='--',
label='$\Lambda \leq 20\mathrm{\ 100 Kpc,\ Model\ A}$')
temp = ax1.hist(bcgZ1['lum_100'] + bcgZ1['logm2lI_A'],
bins=24, range=[10.6, 12.6], linewidth=5.0,
orientation='vertical', histtype='step',
color='r', alpha=0.9, normed=1, linestyle='--',
label='$\Lambda > 20\mathrm{\ 100 Kpc,\ Model\ A}$')
## SBP 100
temp = ax1.hist(gamaZ1['lum_100'] + gamaZ1['logm2lI_B'],
bins=24, range=[10.6, 12.6], linewidth=5.0,
orientation='vertical', histtype='step',
color='k', alpha=0.9, normed=1, linestyle=':',
label='$\Lambda \leq 20\mathrm{\ 100 Kpc,\ Model\ B}$')
temp = ax1.hist(bcgZ1['lum_100'] + bcgZ1['logm2lI_B'],
bins=24, range=[10.6, 12.6], linewidth=5.0,
orientation='vertical', histtype='step',
color='r', alpha=0.9, normed=1, linestyle=':',
label='$\Lambda \leq 20\mathrm{\ 100 Kpc,\ Model\ B}$')
## SBP 100
temp = ax1.hist(gamaZ1['lum_100'] + gamaZ1['logm2lI_C'],
bins=24, range=[10.6, 12.6], linewidth=5.0,
orientation='vertical', histtype='step',
color='k', alpha=0.9, normed=1, linestyle='-',
label='$\Lambda \leq 20\mathrm{\ 100 Kpc,\ Model\ C}$')
temp = ax1.hist(bcgZ1['lum_100'] + bcgZ1['logm2lI_C'],
bins=24, range=[10.6, 12.6], linewidth=5.0,
orientation='vertical', histtype='step',
color='r', alpha=0.9, normed=1, linestyle='-',
label='$\Lambda \leq 20\mathrm{\ 100 Kpc,\ Model\ C}$')
ax1.legend(loc=(0.73, 0.65), shadow=True, fancybox=True,
numpoints=1, fontsize=16, scatterpoints=1,
markerscale=1.2, borderpad=0.25, handletextpad=0.25)
ax1.set_xlim(10.55, 12.59)
ax1.set_ylim(0.01, 2.69)
# Axes setup
# Minor Ticks on
ax1.minorticks_on()
# Axes Thickness
for axis in ['top','bottom','left','right']:
ax1.spines[axis].set_linewidth(3.5)
# Tick Label Size
for tick in ax1.xaxis.get_major_ticks():
tick.label.set_fontsize(24)
for tick in ax1.yaxis.get_major_ticks():
tick.label.set_fontsize(24)
# Tick Length and Width
ax1.tick_params('both', length=10, width=3.0, which='major')
ax1.tick_params('both', length=6, width=2.5, which='minor')
# Label
ax1.set_xlabel('$\log\ (M_{\star}/M_{\odot})\ \mathrm{(HSC)}$', size=35)
ax1.set_ylabel('$\mathrm{Normalized\ \#}$', size=39)
fig.savefig('../figure/hscMassive_mass_distri_norm_mod.png', dpi=300)
plt.show()
In [13]:
fig = plt.figure(figsize=(14, 8))
rec = [0.11, 0.12, 0.87, 0.87]
ax1 = plt.axes(rec)
ax1.axvline(11.4, linewidth=4.0, color='g', linestyle='-',
alpha=0.2, zorder=0)
# ---------------------------------------------------------------------------
# Histogram
## cModel
nGama1, bGama1, pGama1 = ax1.hist(gamaZ1['lumI_cmodel'] + gamaZ1['logm2lI_B'],
bins=20, range=[10.6, 12.6],
orientation='vertical', histtype='stepfilled',
color='k', alpha=0.2, normed=0,
label='$\Lambda \leq 20\mathrm{\ cModel},\ z=0.40$')
nBcg1, bBcg1, pBcg1 = ax1.hist(bcgZ1['lumI_cmodel'] + bcgZ1['logm2lI_B'],
bins=8, range=[10.8, 12.4],
orientation='vertical', histtype='stepfilled',
color='r', alpha=0.2, normed=0,
label='$\Lambda > 20\mathrm{\ cModel},\ z=0.40$')
## SBP 100
nGama2, bGama2, pGama2 = ax1.hist(gamaZ1['lum_100'] + gamaZ1['logm2lI_B'],
bins=20, range=[10.6, 12.6], linewidth=5.0,
orientation='vertical', histtype='step',
color='k', alpha=0.9, normed=0,
label='$\Lambda \leq 20\mathrm{\ 100 Kpc},\ z=0.40$')
nBcg2, bBcg2, pBcg2 = ax1.hist(bcgZ1['lum_100'] + bcgZ1['logm2lI_B'],
bins=8, range=[10.8, 12.4], linewidth=5.0,
orientation='vertical', histtype='step',
color='r', alpha=0.9, normed=0,
label='$\Lambda > 20\mathrm{\ 100 Kpc},\ z=0.40$')
## SBP 100; up to z=0.55
nGama2b, bGama2b, pGama2b = ax1.hist(gamaZ2['lum_100'] + gamaZ2['logm2lI_B'],
bins=20, range=[10.6, 12.6], linewidth=5.0,
orientation='vertical', histtype='step',
color='k', alpha=0.9, normed=0, linestyle=':',
label='$\Lambda \leq 20\mathrm{\ 100 Kpc},\ z=0.50$')
nBcg2b, bBcg2b, pBcg2b = ax1.hist(bcgZ2['lum_100'] + bcgZ2['logm2lI_B'],
bins=8, range=[10.8, 12.4], linewidth=5.0,
orientation='vertical', histtype='step',
color='r', alpha=0.9, normed=0, linestyle=':',
label='$\Lambda > 20\mathrm{\ 100 Kpc},\ z=0.50$')
## SBP 150
nGama3, bGama3, pGama3 = ax1.hist(gamaZ1['logms_gama'] + np.log10(gamaZ1['fluxscale_gama']),
bins=20, range=[10.6, 12.6], linewidth=5.0,
orientation='vertical', histtype='step',
color='k', alpha=0.9, normed=0, linestyle='--',
label='$\Lambda \leq 20\mathrm{\ GAMA},\ z=0.40$')
nBcg3, bBcg3, pBcg3 = ax1.hist(bcgZ1['logms_gama'] + np.log10(bcgZ1['fluxscale_gama']),
bins=8, range=[10.8, 12.4], linewidth=5.0,
orientation='vertical', histtype='step',
color='r', alpha=0.9, normed=0, linestyle='--',
label='$\Lambda > 20\mathrm{\ GAMA},\ z=0.40$')
ax1.legend(loc=(0.63, 0.42), shadow=True, fancybox=True,
numpoints=1, fontsize=22, scatterpoints=1,
markerscale=1.2, borderpad=0.25, handletextpad=0.25)
ax1.set_xlim(10.51, 12.59)
ax1.set_ylim(1, 1690)
# Axes setup
# Minor Ticks on
ax1.minorticks_on()
# Axes Thickness
for axis in ['top','bottom','left','right']:
ax1.spines[axis].set_linewidth(3.5)
# Tick Label Size
for tick in ax1.xaxis.get_major_ticks():
tick.label.set_fontsize(24)
for tick in ax1.yaxis.get_major_ticks():
tick.label.set_fontsize(24)
# Tick Length and Width
ax1.tick_params('both', length=10, width=3.0, which='major')
ax1.tick_params('both', length=6, width=2.5, which='minor')
# Label
ax1.set_xlabel('$\log\ (M_{\star}/M_{\odot})\ \mathrm{(HSC)}$', size=35)
ax1.set_ylabel('$\mathrm{Number}$', size=39)
fig.savefig('../figure/hscMassive_mass_distri_modB.png', dpi=300)
plt.show()
In [14]:
binSize1 = bGama1[1] - bGama1[0]
binSize2 = bBcg1[1] - bBcg1[0]
print(binSize1, binSize2)
dnGama1 = np.log10(nGama1 / (dGama1 * binSize1))
dnBcg1 = np.log10(nBcg1 / (dBcg1 * binSize2))
dnGama2 = np.log10(nGama2 / (dGama1 * binSize1))
dnBcg2 = np.log10(nBcg2 / (dBcg1 * binSize2))
dnGama3 = np.log10(nGama3 / (dGama1 * binSize1))
dnBcg3 = np.log10(nBcg3 / (dBcg1 * binSize2))
dnGama2b = np.log10(nGama2b / (dGama2 * binSize1))
dnBcg2b = np.log10(nBcg2b / (dBcg2 * binSize2))
In [15]:
fig = plt.figure(figsize=(13, 11))
rec = [0.12, 0.12, 0.87, 0.87]
ax1 = plt.axes(rec)
ax1.axvline(11.4, linewidth=4.0, color='g', linestyle='-',
alpha=0.2, zorder=0)
# Leauthaud+2015
ax1.plot(massArr, np.log10(doubleSchechter(massArr)), c='k',
linestyle='-', linewidth=13.0, alpha=0.15,
label='$\mathrm{Leauthaud\ et\ al.\ (2015)}$')
# Moustakas+2013
ax1.plot(mfM13['logM'], mfM13['logPhi'], linestyle='-', linewidth=13.0, alpha=0.15,
c='b', label='$\mathrm{Moustakas\ et\ al.\ (2013)}$')
# Bernardi+2013; SerExp
ax1.plot(mfB13['col1'], mfB13['col2'], linestyle='-', linewidth=13.0, alpha=0.15,
c='g', label='$\mathrm{Bernardi\ et\ al.\ (2013)}$')
# d'Souza+2015;
ax1.plot(mfD15['col1'], mfD15['col2'], linestyle='-', linewidth=13.0, alpha=0.15,
c='orange', label='$\mathrm{d\'Souza\ et\ al.\ (2015)}$')
# GAMA
ax1.plot((bGama1[1:] + bGama1[0:-1])/2.0, dnGama1, c='k',
linewidth=8.0, label='$\Lambda \leq 20\mathrm{\ cModel},\ z=0.40$',
linestyle='-', alpha=0.5)
ax1.plot((bGama2[1:] + bGama2[0:-1])/2.0, dnGama2, c='k',
linewidth=4.0, label='$\Lambda \leq 20\mathrm{\ 100 Kpc},\ z=0.40$',
linestyle='-', alpha=1.0)
ax1.plot((bGama2b[1:] + bGama2b[0:-1])/2.0, dnGama2b, c='k',
linewidth=5.0, label='$\Lambda \leq 20\mathrm{\ 100 Kpc},\ z=0.50$',
linestyle='--')
ax1.plot((bGama3[1:] + bGama3[0:-1])/2.0, dnGama3, c='k',
linewidth=4.0, label='$\Lambda \leq 20\mathrm{\ GAMA},\ z=0.40$',
linestyle=':', alpha=0.8)
# BCG
ax1.plot((bBcg1[1:] + bBcg1[0:-1])/2.0, dnBcg1, c='r',
linewidth=8.0, label='$\Lambda > 20\mathrm{\ cModel},\ z=0.40$',
linestyle='-', alpha=0.5)
ax1.plot((bBcg2[1:] + bBcg2[0:-1])/2.0, dnBcg2, c='r',
linewidth=4.0, label='$\Lambda > 20\mathrm{\ 100 Kpc},\ z=0.40$',
linestyle='-', alpha=1.0)
ax1.plot((bBcg2b[1:] + bBcg2b[0:-1])/2.0, dnBcg2b, c='r',
linewidth=5.0, label='$\Lambda > 20\mathrm{\ 100 Kpc},\ z=0.50$',
linestyle='--')
ax1.plot((bBcg3[1:] + bBcg3[0:-1])/2.0, dnBcg3, c='r',
linewidth=4.0, label='$\Lambda > 20\mathrm{\ GAMA},\ z=0.40$',
linestyle=':', alpha=0.8)
# Axes setup
# Minor Ticks on
ax1.minorticks_on()
# Axes Thickness
for axis in ['top','bottom','left','right']:
ax1.spines[axis].set_linewidth(3.5)
# Tick Label Size
for tick in ax1.xaxis.get_major_ticks():
tick.label.set_fontsize(24)
for tick in ax1.yaxis.get_major_ticks():
tick.label.set_fontsize(24)
# Tick Length and Width
ax1.tick_params('both', length=10, width=3.0, which='major')
ax1.tick_params('both', length=6, width=2.5, which='minor')
ax1.legend(loc=(0.66, 0.50), shadow=True, fancybox=True,
numpoints=1, fontsize=19, scatterpoints=1,
markerscale=1.2, borderpad=0.3, handletextpad=0.34)
ax1.text(0.52, 0.93, '$\mathrm{BC03\ SSP\ +\ Chabrier\ IMF}$',
verticalalignment='bottom', horizontalalignment='left',
fontsize=30.0, transform=ax1.transAxes)
ax1.set_xlim(10.51, 12.59)
ax1.set_ylim(-6.9, -2.41)
# Label
ax1.set_xlabel('$\log\ (M_{\star}/M_{\odot})\ \mathrm{(HSC)}$', size=35)
ax1.set_ylabel('$dN/d\logM_{\star}\ [{\mathrm{Mpc}^{-3}}{\mathrm{dex}^{-1}}]$',
size=39)
fig.savefig('../figure/hscMassive_mass_function_modB.png', dpi=300)
plt.show()
In [16]:
fig = plt.figure(figsize=(14, 8))
rec = [0.11, 0.12, 0.87, 0.87]
ax1 = plt.axes(rec)
ax1.axvline(11.4, linewidth=4.0, color='g', linestyle='-',
alpha=0.2, zorder=0)
# ---------------------------------------------------------------------------
# Histogram
## cModel
nGama1, bGama1, pGama1 = ax1.hist(gamaZ1['lumI_cmodel'] + gamaZ1['logm2lI_C'],
bins=20, range=[10.6, 12.6],
orientation='vertical', histtype='stepfilled',
color='k', alpha=0.2, normed=0,
label='$\Lambda \leq 20\mathrm{\ cModel},\ z=0.40$')
nBcg1, bBcg1, pBcg1 = ax1.hist(bcgZ1['lumI_cmodel'] + bcgZ1['logm2lI_C'],
bins=8, range=[10.8, 12.4],
orientation='vertical', histtype='stepfilled',
color='r', alpha=0.2, normed=0,
label='$\Lambda > 20\mathrm{\ cModel},\ z=0.40$')
## SBP 100
nGama2, bGama2, pGama2 = ax1.hist(gamaZ1['m100_c'],
bins=20, range=[10.6, 12.6], linewidth=5.0,
orientation='vertical', histtype='step',
color='k', alpha=0.9, normed=0,
label='$\Lambda \leq 20\mathrm{\ 100 kpc},\ z=0.40$')
nBcg2, bBcg2, pBcg2 = ax1.hist(bcgZ1['m100_c'],
bins=8, range=[10.8, 12.4], linewidth=5.0,
orientation='vertical', histtype='step',
color='r', alpha=0.9, normed=0,
label='$\Lambda > 20\mathrm{\ 100 kpc},\ z=0.40$')
## SBP 100; up to z=0.55
nGama2b, bGama2b, pGama2b = ax1.hist(gamaZ2['m100_c'],
bins=20, range=[10.6, 12.6], linewidth=5.0,
orientation='vertical', histtype='step',
color='k', alpha=0.9, normed=0, linestyle=':',
label='$\Lambda \leq 20\mathrm{\ 100 kpc},\ z=0.55$')
nBcg2b, bBcg2b, pBcg2b = ax1.hist(bcgZ2['m100_c'],
bins=8, range=[10.8, 12.4], linewidth=5.0,
orientation='vertical', histtype='step',
color='r', alpha=0.9, normed=0, linestyle=':',
label='$\Lambda > 20\mathrm{\ 100 kpc},\ z=0.55$')
## SBP 100; use ModelA
nGama2c, bGama2c, pGama2c = ax1.hist(gamaZ2['m100_a'],
bins=20, range=[10.6, 12.6], linewidth=5.0,
orientation='vertical', histtype='step',
color='k', alpha=0.9, normed=0, linestyle=':',
label='$\Lambda \leq 20\mathrm{\ 100 kpc},\ \mathrm{modA}$')
nBcg2c, bBcg2c, pBcg2c = ax1.hist(bcgZ2['m100_a'],
bins=8, range=[10.8, 12.4], linewidth=5.0,
orientation='vertical', histtype='step',
color='r', alpha=0.9, normed=0, linestyle=':',
label='$\Lambda > 20\mathrm{\ 100 kpc},\ \mathrm{modA}$')
## SBP 150
nGama3, bGama3, pGama3 = ax1.hist(gamaZ1['m150_c'],
bins=20, range=[10.6, 12.6], linewidth=5.0,
orientation='vertical', histtype='step',
color='k', alpha=0.9, normed=0, linestyle='--',
label='$\Lambda \leq 20\mathrm{\ 150kpc},\ z=0.40$')
nBcg3, bBcg3, pBcg3 = ax1.hist(bcgZ1['m150_c'],
bins=8, range=[10.8, 12.4], linewidth=5.0,
orientation='vertical', histtype='step',
color='r', alpha=0.9, normed=0, linestyle='--',
label='$\Lambda > 20\mathrm{\ 150kpc},\ z=0.40$')
ax1.legend(loc=(0.66, 0.27), shadow=True, fancybox=True,
numpoints=1, fontsize=22, scatterpoints=1,
markerscale=1.2, borderpad=0.25, handletextpad=0.25)
ax1.set_xlim(10.51, 12.59)
ax1.set_ylim(1, 1690)
# Axes setup
# Minor Ticks on
ax1.minorticks_on()
# Axes Thickness
for axis in ['top','bottom','left','right']:
ax1.spines[axis].set_linewidth(3.5)
# Tick Label Size
for tick in ax1.xaxis.get_major_ticks():
tick.label.set_fontsize(24)
for tick in ax1.yaxis.get_major_ticks():
tick.label.set_fontsize(24)
# Tick Length and Width
ax1.tick_params('both', length=10, width=3.0, which='major')
ax1.tick_params('both', length=6, width=2.5, which='minor')
# Label
ax1.set_xlabel('$\log\ (M_{\star}/M_{\odot})\ \mathrm{(HSC)}$', size=35)
ax1.set_ylabel('$\mathrm{Number}$', size=39)
fig.savefig('../figure/hscMassive_mass_distri.png', dpi=300)
plt.show()
In [17]:
fig = plt.figure(figsize=(14, 8))
rec = [0.11, 0.12, 0.87, 0.87]
ax1 = plt.axes(rec)
ax1.axvline(11.4, linewidth=4.0, color='g', linestyle='-',
alpha=0.2, zorder=0)
# ---------------------------------------------------------------------------
# Histogram
## cModel
nBcg4, bBcg4, pBcg4 = ax1.hist(bcgZ1['lumI_cmodel'] + bcgZ1['logm2lI_C'],
bins=20, range=[10.6, 12.6],
orientation='vertical', histtype='stepfilled',
color='r', alpha=0.2, normed=0,
label='$\Lambda > 20\mathrm{\ cModel},\ z=0.40$')
nMem4, bMem4, pMem4 = ax1.hist(memZ1['lumI_cmodel'] + memZ1['logm2lI_C'],
bins=20, range=[10.6, 12.6],
orientation='vertical', histtype='stepfilled',
color='g', alpha=0.2, normed=0,
label='$\Lambda > 20\mathrm{\ cModel},\ z=0.40$')
nBcg5, bBcg5, pBcg5 = ax1.hist(bcgZ1['m100_c'],
bins=20, range=[10.6, 12.6],
orientation='vertical', histtype='step',
color='r', alpha=0.2, normed=0,
label='$\Lambda > 20\mathrm{\ cModel},\ z=0.40$')
nMem5, bMem5, pMem5 = ax1.hist(memZ1['m100_c'],
bins=20, range=[10.6, 12.6],
orientation='vertical', histtype='step',
color='g', alpha=0.2, normed=0,
label='$\Lambda > 20\mathrm{\ cModel},\ z=0.40$')
ax1.legend(loc=(0.66, 0.27), shadow=True, fancybox=True,
numpoints=1, fontsize=22, scatterpoints=1,
markerscale=1.2, borderpad=0.25, handletextpad=0.25)
ax1.set_xlim(10.51, 12.59)
ax1.set_ylim(1, 390)
# Axes setup
# Minor Ticks on
ax1.minorticks_on()
# Axes Thickness
for axis in ['top','bottom','left','right']:
ax1.spines[axis].set_linewidth(3.5)
# Tick Label Size
for tick in ax1.xaxis.get_major_ticks():
tick.label.set_fontsize(24)
for tick in ax1.yaxis.get_major_ticks():
tick.label.set_fontsize(24)
# Tick Length and Width
ax1.tick_params('both', length=10, width=3.0, which='major')
ax1.tick_params('both', length=6, width=2.5, which='minor')
# Label
ax1.set_xlabel('$\log\ (M_{\star}/M_{\odot})\ \mathrm{(HSC)}$', size=35)
ax1.set_ylabel('$\mathrm{Number}$', size=39)
#fig.savefig('../figure/hscMassive_mass_distri.png', dpi=300)
plt.show()
In [18]:
binSize1 = bGama1[1] - bGama1[0]
binSize2 = bBcg1[1] - bBcg1[0]
print(binSize1, binSize2)
dnGama1 = np.log10(nGama1 / (dGama1 * binSize1))
dnBcg1 = np.log10(nBcg1 / (dBcg1 * binSize2))
dnGama2 = np.log10(nGama2 / (dGama1 * binSize1))
dnBcg2 = np.log10(nBcg2 / (dBcg1 * binSize2))
dnGama3 = np.log10(nGama3 / (dGama1 * binSize1))
dnBcg3 = np.log10(nBcg3 / (dBcg1 * binSize2))
dnGama2b = np.log10(nGama2b / (dGama2 * binSize1))
dnBcg2b = np.log10(nBcg2b / (dBcg2 * binSize2))
dnGama2c = np.log10(nGama2c / (dGama1 * binSize1))
dnBcg2c = np.log10(nBcg2c / (dBcg1 * binSize2))
In [19]:
dnGama4 = nGama1 / (dGama1 * binSize1)
dnGama5 = nGama2 / (dGama2 * binSize1)
dnBcg4 = nBcg4 / (dBcg1 * binSize1)
dnMem4 = nMem4 / (dMem1 * binSize1)
dnBcg5 = nBcg5 / (dBcg1 * binSize1)
dnMem5 = nMem5 / (dMem1 * binSize1)
dnAll4 = np.log10(dnGama4 + dnBcg4 + dnMem4)
dnAll5 = np.log10(dnGama5 + dnBcg5 + dnMem5)
nGamaError = np.sqrt(nGama2)
nBcgError = np.sqrt(nBcg5)
nMemError = np.sqrt(nMem5)
In [20]:
nBcgError[nBcgError < 1.0] = 1.0
nMemError[nMemError < 1.0] = 1.0
nGamaError[nGamaError < 1.0] = 1.0
In [21]:
dnGamaUpp = (nGama2 + nGamaError) / (dGama2 * binSize1)
dnBcgUpp = (nBcg5 + nBcgError) / (dBcg1 * binSize1)
dnMemUpp = (nMem5 + nMemError) / (dMem1 * binSize1)
dnAllUpp = np.log10(dnGamaUpp + dnBcgUpp + dnMemUpp)
dnAllLow = dnAll5 - (dnAllUpp - dnAll5)
In [23]:
fig = plt.figure(figsize=(13, 11))
rec = [0.12, 0.12, 0.87, 0.87]
ax1 = plt.axes(rec)
ax1.axvline(11.5, linewidth=5.0, color='k', linestyle='--',
alpha=0.1, zorder=0)
hTerm1 = (2.0 * np.log10(0.72))
hTerm2 = (3.0 * np.log10(0.72))
# GAMA
ax1.plot((bGama1[1:] + bGama1[0:-1])/2.0, dnGama1, c=BLUE0,
linewidth=9.0, label='$\Lambda \leq 20\mathrm{\ cModel},\ z=0.40$',
linestyle='-', alpha=0.5)
ax1.plot((bGama2[1:] + bGama2[0:-1])/2.0, dnGama2, c=BLUE1,
linewidth=4.0, label='$\Lambda \leq 20\mathrm{\ 100 kpc},\ z=0.40$',
linestyle='-', alpha=1.0)
ax1.plot((bGama2b[1:] + bGama2b[0:-1])/2.0, dnGama2b, c=BLUE1,
linewidth=4.0, label='$\Lambda \leq 20\mathrm{\ 100 kpc},\ z=0.50$',
linestyle='--', alpha=0.8)
ax1.plot((bGama2c[1:] + bGama2c[0:-1])/2.0, dnGama2b, c=BLUE1,
linewidth=4.0, label='$\Lambda \leq 20\mathrm{\ 100 kpc},\ \mathrm{modA}$',
linestyle='-.', alpha=0.8)
ax1.plot((bGama3[1:] + bGama3[0:-1])/2.0, dnGama3, c='k',
linewidth=4.0, label='$\Lambda \leq 20\mathrm{\ 150 kpc},\ z=0.40$',
linestyle=':', alpha=0.8)
# BCG
ax1.plot((bBcg1[1:] + bBcg1[0:-1])/2.0, dnBcg1, c=RED0,
linewidth=9.0, label='$\Lambda > 20\mathrm{\ cModel},\ z=0.40$',
linestyle='-', alpha=0.5)
ax1.plot((bBcg2[1:] + bBcg2[0:-1])/2.0, dnBcg2, c=RED1,
linewidth=4.0, label='$\Lambda > 20\mathrm{\ 100 kpc},\ z=0.40$',
linestyle='-', alpha=1.0)
ax1.plot((bBcg2b[1:] + bBcg2b[0:-1])/2.0, dnBcg2b, c=RED1,
linewidth=4.0, label='$\Lambda > 20\mathrm{\ 100 kpc},\ z=0.50$',
linestyle='--')
ax1.plot((bBcg2c[1:] + bBcg2c[0:-1])/2.0, dnBcg2c, c=RED1,
linewidth=4.0, label='$\Lambda > 20\mathrm{\ 100 kpc},\ \mathrm{modA}$',
linestyle='-.')
ax1.plot((bBcg3[1:] + bBcg3[0:-1])/2.0, dnBcg3, c='r',
linewidth=4.0, label='$\Lambda > 20\mathrm{\ 150 kpc},\ z=0.40$',
linestyle=':', alpha=0.8)
# Axes setup
# Minor Ticks on
ax1.minorticks_on()
# Axes Thickness
for axis in ['top','bottom','left','right']:
ax1.spines[axis].set_linewidth(3.5)
# Tick Label Size
for tick in ax1.xaxis.get_major_ticks():
tick.label.set_fontsize(24)
for tick in ax1.yaxis.get_major_ticks():
tick.label.set_fontsize(24)
# Tick Length and Width
ax1.tick_params('both', length=10, width=3.0, which='major')
ax1.tick_params('both', length=6, width=2.5, which='minor')
ax1.legend(loc=(0.66, 0.50), shadow=True, fancybox=True,
numpoints=1, fontsize=19, scatterpoints=1,
markerscale=1.2, borderpad=0.3, handletextpad=0.34)
#ax1.text(0.52, 0.93, '$\mathrm{BC03\ SSP\ +\ Chabrier\ IMF}$',
# verticalalignment='bottom', horizontalalignment='left',
# fontsize=30.0, transform=ax1.transAxes)
ax1.set_xlim(10.51, 12.59)
ax1.set_ylim(-6.9, -2.41)
# Label
ax1.set_xlabel('$\log\ (M_{\star}/M_{\odot})\ \mathrm{(HSC)}$', size=35)
ax1.set_ylabel('$dN/d\logM_{\star}\ [{\mathrm{Mpc}^{-3}}{\mathrm{dex}^{-1}}]$',
size=39)
fig.savefig('../figure/hscMassive_mass_function_3.png', dpi=300)
plt.show()
In [122]:
fig = plt.figure(figsize=(13, 11))
rec = [0.12, 0.12, 0.87, 0.87]
ax1 = plt.axes(rec)
ax1.axvline(11.5, linewidth=5.0, color='k', linestyle='--',
alpha=0.1, zorder=0)
hTerm1 = (2.0 * np.log10(0.72))
hTerm2 = (3.0 * np.log10(0.72))
# GAMA
"""
ax1.plot((bGama1[1:] + bGama1[0:-1])/2.0, dnGama1, c=BLUE0,
linewidth=9.0, label='$\Lambda \leq 20\mathrm{\ cModel},\ z=0.40$',
linestyle='-', alpha=0.5)
ax1.plot((bGama2[1:] + bGama2[0:-1])/2.0, dnGama2, c=BLUE1,
linewidth=4.0, label='$\Lambda \leq 20\mathrm{\ 100 kpc},\ z=0.40$',
linestyle='-', alpha=1.0)
ax1.plot((bGama2b[1:] + bGama2b[0:-1])/2.0, dnGama2b, c=BLUE1,
linewidth=4.0, label='$\Lambda \leq 20\mathrm{\ 100 kpc},\ z=0.50$',
linestyle='--', alpha=0.8)
ax1.plot((bGama2c[1:] + bGama2c[0:-1])/2.0, dnGama2b, c=BLUE1,
linewidth=4.0, label='$\Lambda \leq 20\mathrm{\ 100 kpc},\ \mathrm{modA}$',
linestyle='-.', alpha=0.8)
ax1.plot((bGama3[1:] + bGama3[0:-1])/2.0, dnGama3, c='k',
linewidth=4.0, label='$\Lambda \leq 20\mathrm{\ 150 kpc},\ z=0.40$',
linestyle=':', alpha=0.8)
#BCG
ax1.plot((bBcg1[1:] + bBcg1[0:-1])/2.0, dnBcg1, c=RED0,
linewidth=9.0, label='$\Lambda > 20\mathrm{\ cModel},\ z=0.40$',
linestyle='-', alpha=0.5)
ax1.plot((bBcg2[1:] + bBcg2[0:-1])/2.0, dnBcg2, c=RED1,
linewidth=4.0, label='$\Lambda > 20\mathrm{\ 100 kpc},\ z=0.40$',
linestyle='-', alpha=1.0)
ax1.plot((bBcg2b[1:] + bBcg2b[0:-1])/2.0, dnBcg2b, c=RED1,
linewidth=4.0, label='$\Lambda > 20\mathrm{\ 100 kpc},\ z=0.50$',
linestyle='--')
ax1.plot((bBcg2c[1:] + bBcg2c[0:-1])/2.0, dnBcg2c, c=RED1,
linewidth=4.0, label='$\Lambda > 20\mathrm{\ 100 kpc},\ \mathrm{modA}$',
linestyle='-.')
ax1.plot((bBcg3[1:] + bBcg3[0:-1])/2.0, dnBcg3, c='r',
linewidth=4.0, label='$\Lambda > 20\mathrm{\ GAMA},\ z=0.40$',
linestyle=':', alpha=0.8)
ax1.plot((bGama1[1:] + bGama1[0:-1])/2.0, np.log10(dnMem4), c='g',
linewidth=9.0, label='$\Lambda > 20\mathrm{\ cModel},\ \mathrm{Mem}$',
linestyle='--', alpha=0.5)
ax1.plot((bGama1[1:] + bGama1[0:-1])/2.0, np.log10(dnMem5), c='g',
linewidth=5.0, label='$\Lambda > 20\mathrm{\ 100 kpc},\ \mathrm{Mem}$',
linestyle='-', alpha=0.5)
ax1.plot((bGama1[1:] + bGama1[0:-1])/2.0, dnAll4, c='k',
linewidth=9.0, label='$\Lambda > 20\mathrm{\ cModel},\ \mathrm{All}$',
linestyle='--', alpha=0.5)
"""
ax1.fill_between((bGama1[1:] + bGama1[0:-1])/2.0, dnAllLow, dnAllUpp,
alpha=0.15, facecolor='k', edgecolor='none', label=None)
ax1.plot((bGama1[1:] + bGama1[0:-1])/2.0, dnAll5, c='k',
linewidth=4.5, label='$\mathrm{This\ Work;}\ \mathrm{\ 100\ kpc}$',
linestyle='-', alpha=0.9)
# Leauthaud+2015
ax1.plot(massArr, np.log10(doubleSchechter(massArr)), c='c',
linestyle='-', linewidth=5.0, alpha=0.70,
label='$\mathrm{Leauthaud\ et\ al.\ (2015)}$')
# Moustakas+2013
ax1.plot(mfM13['logM'], mfM13['logPhi'], linestyle='-',
linewidth=6.0, alpha=0.8,
c='b', label='$\mathrm{Moustakas\ et\ al.\ (2013)}$')
# Bernardi+2013; SerExp
ax1.plot(mfB13['col1'] - hTerm1, mfB13['col2'] + hTerm2, linestyle='-',
linewidth=6.0, alpha=0.80,
c='r', label='$\mathrm{Bernardi\ et\ al.\ (2013)}$')
# d'Souza+2015;
ax1.plot(mfD15['col1'] - hTerm1, mfD15['col2'] + hTerm2, linestyle='-',
linewidth=6.0, alpha=0.80,
c='g', label='$\mathrm{d\'Souza\ et\ al.\ (2015)}$')
# Axes setup
# Minor Ticks on
ax1.minorticks_on()
# Axes Thickness
for axis in ['top','bottom','left','right']:
ax1.spines[axis].set_linewidth(3.5)
# Tick Label Size
for tick in ax1.xaxis.get_major_ticks():
tick.label.set_fontsize(24)
for tick in ax1.yaxis.get_major_ticks():
tick.label.set_fontsize(24)
# Tick Length and Width
ax1.tick_params('both', length=10, width=3.0, which='major')
ax1.tick_params('both', length=6, width=2.5, which='minor')
ax1.legend(loc=(0.61, 0.70), shadow=True, fancybox=True,
numpoints=1, fontsize=21, scatterpoints=1,
markerscale=1.2, borderpad=0.4, handletextpad=0.5)
#ax1.text(0.52, 0.93, '$\mathrm{BC03\ SSP\ +\ Chabrier\ IMF}$',
# verticalalignment='bottom', horizontalalignment='left',
# fontsize=30.0, transform=ax1.transAxes)
ax1.set_xlim(10.51, 12.59)
ax1.set_ylim(-6.9, -2.41)
# Label
ax1.set_xlabel('$\log\ (M_{\star}/M_{\odot})\ \mathrm{(HSC)}$', size=35)
ax1.set_ylabel('$dN/d\logM_{\star}\ [{\mathrm{Mpc}^{-3}}{\mathrm{dex}^{-1}}]$',
size=39)
fig.savefig('../figure/hscMassive_mass_function_4.png', dpi=300)
plt.show()
In [ ]: