In [1]:
%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
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()
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
# 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, Blues_5
ORG4 = Oranges_4.mpl_colormap
BLU5 = Blues_5.mpl_colormap
In [2]:
# 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
# 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)
In [9]:
# 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,
index=None):
""" 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 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
In [84]:
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
#
bcgDir = os.path.join(newDir, 'redbcg')
memDir = os.path.join(newDir, 'redmem')
gamaDir = os.path.join(newDir, 'gama')
# 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 [93]:
bcgClean = bcgTab[(bcgTab['m100_c'] >= 10.0) &
(bcgTab['c82_120'] <= 14.0) &
(bcgTab['r90_max'] <= 220.0) &
(bcgTab['P_CEN_1'] >= 0.8)]
print(len(bcgClean))
gamaClean = gamaTab[(gamaTab['c82_120'] >= 5.5) &
(gamaTab['c82_120'] <= 14.0) &
(gamaTab['gz_kC'] >= 1.58) &
(gamaTab['m100_c'] >= 10.0) &
(gamaTab['r90_max'] <= 160.0) &
(gamaTab['ur_rest_sed'] >= 2.1) &
((gamaTab['r50_120'] - gamaTab['r20_120']) >=
(32.0 * (gamaTab['m30_c'] - gamaTab['m10_c']) - 1.4))]
print(len(gamaClean))
memClean = memTab[(memTab['c82_120'] >= 5.0) &
(memTab['c82_120'] <= 14.0) &
(memTab['gz_kC'] >= 1.45) &
(memTab['m100_c'] >= 10.0) &
(memTab['r90_max'] <= 160.0) &
((memTab['r50_120'] - memTab['r20_120']) >=
(32.0 * (memTab['m30_c'] - memTab['m10_c']) - 1.4))]
print(len(memClean))
In [94]:
gamaUse = gamaClean[(gamaClean['z_use'] >= 0.20) &
(gamaClean['z_use'] <= 0.45)]
gamaUse2 = gamaClean[(gamaClean['z_use'] >= 0.20) &
(gamaClean['z_use'] <= 0.40)]
print(len(gamaUse))
print(len(gamaUse2))
memUse = memClean[(memClean['z_use'] >= 0.20) &
(memClean['z_use'] <= 0.40)]
print(len(memUse))
bcgUse = bcgClean[(bcgClean['z_use'] >= 0.20) &
(bcgClean['z_use'] <= 0.45) &
(bcgClean['LAMBDA_CLUSTER'] >= 25.0)]
bcgUse2 = bcgClean[(bcgClean['z_use'] >= 0.20) &
(bcgClean['z_use'] <= 0.40) &
(bcgClean['LAMBDA_CLUSTER'] >= 25.0)]
print(len(bcgUse))
print(len(bcgUse2))
In [98]:
print("########################################################\n")
#gamaM1 = gamaClean[(gamaClean['m100_c'] >= 11.56) &
# (gamaClean['m100_c'] < 11.70) &
# (gamaClean['Z'] >= 0.20) &
# (gamaClean['Z'] <= 0.40)]
#gamaM2 = gamaClean[(gamaClean['m100_c'] >= 11.72) &
# (gamaClean['m100_c'] < 11.95) &
# (gamaClean['Z'] >= 0.20) &
# (gamaClean['Z'] <= 0.50)]
gamaM1 = gamaClean[(gamaClean['m100_c'] >= 11.55) &
(gamaClean['m100_c'] < 11.70) &
(gamaClean['z_use'] >= 0.20) &
(gamaClean['z_use'] <= 0.45)]
gamaM2 = gamaClean[(gamaClean['m100_c'] >= 11.71) &
(gamaClean['m100_c'] < 11.90) &
(gamaClean['z_use'] >= 0.20) &
(gamaClean['z_use'] <= 0.50)]
gamaM1b = gamaClean[(gamaClean['m100_c'] >= 11.55) &
(gamaClean['m100_c'] < 11.70) &
(gamaClean['z_use'] >= 0.20) &
(gamaClean['z_use'] <= 0.40)]
gamaM2b = gamaClean[(gamaClean['m100_c'] >= 11.70) &
(gamaClean['m100_c'] < 11.90) &
(gamaClean['z_use'] >= 0.20) &
(gamaClean['z_use'] <= 0.40)]
print(len(gamaM1), np.nanmedian(gamaM1['m100_c']), np.nanmedian(gamaM1['z_use']))
print(len(gamaM2), np.nanmedian(gamaM2['m100_c']), np.nanmedian(gamaM2['z_use']))
print(len(gamaM1b), np.nanmedian(gamaM1b['m100_c']), np.nanmedian(gamaM1b['z_use']))
print(len(gamaM2b), np.nanmedian(gamaM2b['m100_c']), np.nanmedian(gamaM2b['z_use']))
print("########################################################\n")
memM1 = memClean[(memClean['m100_c'] >= 11.55) &
(memClean['m100_c'] < 11.70) &
(memClean['z_use'] >= 0.20) &
(memClean['z_use'] <= 0.40)]
memM2 = memClean[(memClean['m100_c'] >= 11.70) &
(memClean['m100_c'] < 11.90) &
(memClean['z_use'] >= 0.20) &
(memClean['z_use'] <= 0.40)]
print(len(memM1), np.nanmedian(memM1['m100_c']), np.nanmedian(memM1['Z']))
print(len(memM2), np.nanmedian(memM2['m100_c']), np.nanmedian(memM2['Z']))
print("########################################################\n")
#bcgM1 = bcgClean[(bcgClean['m100_c'] >= 11.45) &
# (bcgClean['m100_c'] < 11.68) &
# (bcgClean['Z'] >= 0.20) &
# (bcgClean['Z'] <= 0.40) &
# (bcgClean['P_CEN_1'] >= 0.8)]
#bcgM2 = bcgClean[(bcgClean['m100_c'] >= 11.70) &
# (bcgClean['m100_c'] < 11.95) &
# (bcgClean['Z'] >= 0.20) &
# (bcgClean['Z'] <= 0.50) &
# (bcgClean['P_CEN_1'] >= 0.8) &
# (bcgClean['LAMBDA_CLUSTER'] >= 30.0)]
bcgM1 = bcgClean[(bcgClean['m100_c'] >= 11.55) &
(bcgClean['m100_c'] < 11.70) &
(bcgClean['z_use'] >= 0.20) &
(bcgClean['z_use'] <= 0.40) &
(bcgClean['LAMBDA_CLUSTER'] >= 25.0)]
bcgM2 = bcgClean[(bcgClean['m100_c'] >= 11.70) &
(bcgClean['m100_c'] < 11.90) &
(bcgClean['z_use'] >= 0.20) &
(bcgClean['z_use'] <= 0.50) &
(bcgClean['LAMBDA_CLUSTER'] >= 25.0)]
bcgM1b = bcgClean[(bcgClean['m100_c'] >= 11.50) &
(bcgClean['m100_c'] < 11.70) &
(bcgClean['z_use'] >= 0.20) &
(bcgClean['z_use'] <= 0.40) &
(bcgClean['LAMBDA_CLUSTER'] >= 25.0)]
bcgM2b = bcgClean[(bcgClean['m100_c'] >= 11.70) &
(bcgClean['m100_c'] < 11.90) &
(bcgClean['z_use'] >= 0.20) &
(bcgClean['z_use'] <= 0.40) &
(bcgClean['LAMBDA_CLUSTER'] >= 25.0)]
print(len(bcgM1), np.nanmedian(bcgM1['m100_c']), np.nanmedian(bcgM1['z_use']))
print(len(bcgM2), np.nanmedian(bcgM2['m100_c']), np.nanmedian(bcgM2['z_use']))
print(len(bcgM1b), np.nanmedian(bcgM1b['m100_c']), np.nanmedian(bcgM1b['z_use']))
print(len(bcgM2b), np.nanmedian(bcgM2b['m100_c']), np.nanmedian(bcgM2b['z_use']))
In [51]:
print("########################################################\n")
gamaCM1 = gamaClean[(gamaClean['m10_c'] >= 11.19) &
(gamaClean['m10_c'] < 11.35) &
(gamaClean['Z'] >= 0.20) &
(gamaClean['Z'] <= 0.40)]
gamaCM1b = gamaCM1[(gamaCM1['c82_120'] >= 7.0) &
(gamaCM1['r90_120'] / gamaCM1['r50_120'] >= 2.3) &
(gamaCM1['gz_kC'] >= 1.6) &
(gamaCM1['ur_rest_sed'] >= 2.2)]
gamaCM2 = gamaClean[(gamaClean['m10_c'] >= 11.35) &
(gamaClean['m10_c'] < 11.55) &
(gamaClean['Z'] >= 0.20) &
(gamaClean['Z'] <= 0.40)]
print(len(gamaCM1), np.nanmedian(gamaCM1['m10_c']),
np.nanmedian(gamaCM1['Z']))
print(len(gamaCM2), np.nanmedian(gamaCM2['m10_c']),
np.nanmedian(gamaCM2['Z']))
print(len(gamaCM1b), np.nanmedian(gamaCM1b['m10_c']),
np.nanmedian(gamaCM1b['Z']))
print("########################################################\n")
memCM1 = memClean[(memClean['m10_c'] >= 11.18) &
(memClean['m10_c'] < 11.35) &
(memClean['Z'] >= 0.20) &
(memClean['Z'] <= 0.40)]
memCM2 = memClean[(memClean['m10_c'] >= 11.35) &
(memClean['m10_c'] < 11.55) &
(memClean['Z'] >= 0.20) &
(memClean['Z'] <= 0.50)]
print(len(memCM1), np.nanmedian(memCM1['m10_c']),
np.nanmedian(memCM1['Z']))
print(len(memCM2), np.nanmedian(memCM2['m10_c']),
np.nanmedian(memCM2['Z']))
print("########################################################\n")
bcgCM1 = bcgClean[(bcgClean['m10_c'] >= 11.15) &
(bcgClean['m10_c'] < 11.34) &
(bcgClean['Z'] >= 0.20) &
(bcgClean['Z'] <= 0.40) &
(bcgClean['LAMBDA_CLUSTER'] >= 30.0)]
bcgCM1b = bcgClean[(bcgClean['m10_c'] >= 11.15) &
(bcgClean['m10_c'] < 11.34) &
(bcgClean['m100_c'] <= 11.8) &
(bcgClean['Z'] >= 0.20) &
(bcgClean['Z'] <= 0.39)]
bcgCM2 = bcgClean[(bcgClean['m10_c'] >= 11.35) &
(bcgClean['m10_c'] < 11.55) &
(bcgClean['Z'] >= 0.20) &
(bcgClean['Z'] <= 0.50) &
(bcgClean['LAMBDA_CLUSTER'] >= 30.0)]
print(len(bcgCM1), np.nanmean(bcgCM1['m10_c']),
np.nanmedian(bcgCM1['Z']))
print(len(bcgCM2), np.nanmean(bcgCM2['m10_c']),
np.nanmedian(bcgCM2['Z']))
print(len(bcgCM1b), np.nanmean(bcgCM1b['m10_c']),
np.nanmedian(bcgCM1b['Z']))
In [66]:
gM1_mass = confidence_interval(gamaM1['m100_c'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gCM1_mass = confidence_interval(gamaM1['m10_c'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gM1_mdif1 = confidence_interval((gamaM1['lum_100'] - gamaM1['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gM1_mdif2 = confidence_interval((gamaM1['lum_75'] - gamaM1['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gM1_mdif3 = confidence_interval((gamaM1['lum_50'] - gamaM1['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gM1_mdif4 = confidence_interval((gamaM1['lum_25'] - gamaM1['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gM1_r20 = confidence_interval(gamaM1['r20_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gM1_r50 = confidence_interval(gamaM1['r50_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gM1_r90 = confidence_interval(gamaM1['r90_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gM1_c82 = confidence_interval(gamaM1['c82_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gM1_c52 = confidence_interval((gamaM1['r50_120'] / gamaM1['r20_120']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gM2_mass = confidence_interval(gamaM2['m100_c'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gCM2_mass = confidence_interval(gamaM2['m10_c'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gM2_mdif1 = confidence_interval((gamaM2['lum_100'] - gamaM2['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gM2_mdif2 = confidence_interval((gamaM2['lum_75'] - gamaM2['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gM2_mdif3 = confidence_interval((gamaM2['lum_50'] - gamaM2['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gM2_mdif4 = confidence_interval((gamaM2['lum_25'] - gamaM2['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gM2_r20 = confidence_interval(gamaM2['r20_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gM2_r50 = confidence_interval(gamaM2['r50_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gM2_r90 = confidence_interval(gamaM2['r90_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gM2_c82 = confidence_interval(gamaM2['c82_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
gM2_c95 = confidence_interval((gamaM2['r50_120'] / gamaM2['r20_120']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
In [67]:
bM1_mass = confidence_interval(bcgM1['m100_c'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bCM1_mass = confidence_interval(bcgM1['m10_c'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bM1_mdif1 = confidence_interval((bcgM1['lum_100'] - bcgM1['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bM1_mdif2 = confidence_interval((bcgM1['lum_75'] - bcgM1['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bM1_mdif3 = confidence_interval((bcgM1['lum_50'] - bcgM1['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bM1_mdif4 = confidence_interval((bcgM1['lum_25'] - bcgM1['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bM1_r20 = confidence_interval(bcgM1['r20_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bM1_r50 = confidence_interval(bcgM1['r50_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bM1_r90 = confidence_interval(bcgM1['r90_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bM1_c82 = confidence_interval(bcgM1['c82_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bM1_c52 = confidence_interval((bcgM1['r50_120'] / bcgM1['r20_120']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bM2_mass = confidence_interval(bcgM2['m100_c'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bCM2_mass = confidence_interval(bcgM2['m10_c'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bM2_mdif1 = confidence_interval((bcgM2['lum_100'] - bcgM2['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bM2_mdif2 = confidence_interval((bcgM2['lum_75'] - bcgM2['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bM2_mdif3 = confidence_interval((bcgM2['lum_50'] - bcgM2['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bM2_mdif4 = confidence_interval((bcgM2['lum_25'] - bcgM2['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bM2_r20 = confidence_interval(bcgM2['r20_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bM2_r50 = confidence_interval(bcgM2['r50_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bM2_r90 = confidence_interval(bcgM2['r90_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bM2_c82 = confidence_interval(bcgM2['c82_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
bM2_c95 = confidence_interval((bcgM2['r50_120'] / bcgM2['r20_120']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
In [68]:
mM1_mass = confidence_interval(memM1['m100_c'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mCM1_mass = confidence_interval(memM1['m10_c'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mM1_mdif1 = confidence_interval((memM1['lum_100'] - memM1['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mM1_mdif2 = confidence_interval((memM1['lum_75'] - memM1['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mM1_mdif3 = confidence_interval((memM1['lum_50'] - memM1['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mM1_mdif4 = confidence_interval((memM1['lum_25'] - memM1['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mM1_r20 = confidence_interval(memM1['r20_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mM1_r50 = confidence_interval(memM1['r50_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mM1_r90 = confidence_interval(memM1['r90_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mM1_c82 = confidence_interval(memM1['c82_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mM1_c52 = confidence_interval((memM1['r50_120'] / memM1['r20_120']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mM2_mass = confidence_interval(memM2['m100_c'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mCM2_mass = confidence_interval(memM2['m10_c'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mM2_mdif1 = confidence_interval((memM2['lum_100'] - memM2['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mM2_mdif2 = confidence_interval((memM2['lum_75'] - memM2['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mM2_mdif3 = confidence_interval((memM2['lum_50'] - memM2['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mM2_mdif4 = confidence_interval((memM2['lum_25'] - memM2['lum_10']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mM2_r20 = confidence_interval(memM2['r20_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mM2_r50 = confidence_interval(memM2['r50_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mM2_r90 = confidence_interval(memM2['r90_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mM2_c82 = confidence_interval(memM2['c82_120'], axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
mM2_c95 = confidence_interval((memM2['r50_120'] / memM2['r20_120']),
axis=0, alpha=[SIGMA1, 1.0],
metric=np.median, numResamples=1000, interpolate=True)
In [69]:
plt.hist(gamaM1['m100_c'], bins=20, range=[11.4, 11.8], color='k',
histtype='step', normed=1)
plt.hist(bcgM1['m100_c'], bins=20, range=[11.4, 11.8], color='r',
histtype='step', normed=1)
plt.show()
In [56]:
plt.hist(gamaM2['m100_c'], bins=20, range=[11.7, 12.2], color='k',
histtype='step', normed=1)
plt.hist(bcgM2['m100_c'], bins=20, range=[11.7, 12.2], color='r',
histtype='step', normed=1)
plt.hist(gamaM2b['m100_c'], bins=20, range=[11.7, 12.2], color='k',
histtype='step', normed=1, linestyle='--')
plt.hist(bcgM2b['m100_c'], bins=20, range=[11.7, 12.2], color='r',
histtype='step', normed=1, linestyle='--')
plt.show()
In [426]:
gamaM1.write(os.path.join(newDir, 'test_gama_m1.fits'), format='fits', overwrite=True)
gamaM2.write(os.path.join(newDir, 'test_gama_m2.fits'), format='fits', overwrite=True)
gamaM2b.write(os.path.join(newDir, 'test_gama_m2b.fits'), format='fits', overwrite=True)
bcgM1.write(os.path.join(newDir, 'test_bcg_m1.fits'), format='fits', overwrite=True)
bcgM2.write(os.path.join(newDir, 'test_bcg_m2.fits'), format='fits', overwrite=True)
bcgM2b.write(os.path.join(newDir, 'test_bcg_m2b.fits'), format='fits', overwrite=True)
memM1.write(os.path.join(newDir, 'test_mem_m1.fits'), format='fits', overwrite=True)
memM2.write(os.path.join(newDir, 'test_mem_m2.fits'), format='fits', overwrite=True)
gamaCM1.write(os.path.join(newDir, 'test_gama_cm1.fits'), format='fits', overwrite=True)
gamaCM2.write(os.path.join(newDir, 'test_gama_cm2.fits'), format='fits', overwrite=True)
bcgCM1.write(os.path.join(newDir, 'test_bcg_cm1.fits'), format='fits', overwrite=True)
bcgCM2.write(os.path.join(newDir, 'test_bcg_cm2.fits'), format='fits', overwrite=True)
In [18]:
### Fancy one
# definitions for the axes
left, width = 0.12, 0.69
right = left + width
bottom, height = 0.12, 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]
In [19]:
# Color
BLUE0 = "#92c5de"
BLUE1 = "#0571b0"
RED0 = "#f4a582"
RED1 = "#ca0020"
PURPLE0 = '#af8dc3'
PURPLE1 = '#762a83'
BROWN0 = '#bf812d'
BROWN1 = '#543005'
GREEN0 = '#7fbf7b'
GREEN1 = '#1b7837'
In [72]:
fig = plt.figure(figsize=(13, 9))
ax1 = plt.axes(recScat)
ax2 = plt.axes(recHist)
# SBP v.s. (cModel - SBP)
# ---------------------------------------------------------------------------
# Scatter plot
#ax1.axhline(0.0, linewidth=4.0, linestyle='-', c='k', alpha=0.2)
ax1.axvline(11.5, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)
ax1.axvline(11.7, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)
ax1.axvline(11.9, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)
# Matched ones
p1 = ax1.scatter(gamaClean['m100_c'],
gamaClean['lum_100'] - gamaClean['lum_10'], s=35.0,
alpha=0.20, facecolor=BLUE0, edgecolor='none',
label='$\Lambda \leq 20\ \mathrm{Central}$')
p2 = ax1.scatter(bcgUse['m100_c'],
bcgUse['lum_100'] - bcgUse['lum_10'], edgecolor='none',
s=((bcgUse['z_use'] - 0.10) * 600.0), cmap=ORG4, alpha=0.90,
c=toColorArr(bcgUse['LAMBDA_CLUSTER'], bottom=20.0, top=70.0),
label='$\Lambda > 20\ \mathrm{Central}$', marker='s')
"""
p3 = ax1.scatter(bcgM1['m100_c'],
bcgM1['lum_100'] - bcgM1['lum_10'], edgecolor='k',
s=((bcgM1['z_use'] - 0.10) * 500.0), alpha=0.95,
facecolor='none', label=None, marker='s', linewidth=1.5)
p4 = ax1.scatter(bcgM2b['m100_c'],
bcgM2b['lum_100'] - bcgM2b['lum_10'], edgecolor='k',
s=((bcgM2['z_use'] - 0.10) * 500.0), alpha=0.95,
facecolor='none', label=None, marker='s', linewidth=1.5)
"""
# M1
ax1.errorbar(gM1_mass[2], gM1_mdif1[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(bM1_mass[2], bM1_mdif1[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)
ax1.scatter(gM1_mass[2], gM1_mdif1[2], marker='^', s=400, facecolor=BLUE1,
edgecolor='k', linewidth=3.0, zorder=102, alpha=0.9,
label='$\mathrm{[11.5,11.7]}\ \Lambda \leq 20$')
ax1.scatter(bM1_mass[2], bM1_mdif1[2], marker='p', s=420, facecolor=RED1,
edgecolor='k', linewidth=3.0, zorder=102,
label='$\mathrm{[11.5,11.7]}\ \Lambda > 30$')
# M2
ax1.errorbar(gM2_mass[2], gM2_mdif1[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(bM2_mass[2], bM2_mdif1[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(gM2_mass[2], gM2_mdif1[2], marker='h', s=420, facecolor=BLUE1,
edgecolor='k', linewidth=3.0, zorder=102, alpha=0.9,
label='$\mathrm{[11.7,11.9]}\ \Lambda \leq 20$')
ax1.scatter(bM2_mass[2], bM2_mdif1[2], marker='8', s=420, facecolor=RED1,
edgecolor='k', linewidth=3.0, zorder=102,
label='$\mathrm{[11.7,11.9]}\ \Lambda > 30$')
# 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[1].set_color(ORG4(0.8))
legend.legendHandles[0].set_sizes([150])
legend.legendHandles[1].set_sizes([200])
#ax1.text(0.05, 0.04, '$\mathrm{Size:}\ {\Lambda}_{\mathrm{redMapper}}$',
# verticalalignment='bottom', horizontalalignment='left',
# fontsize=26.0, transform=ax1.transAxes, color=RED0)
# 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})\ (100\ \mathrm{Kpc})$', size=40)
ax1.set_ylabel('$\Delta(\log M{\star})_{\mathrm{100\ kpc}-\mathrm{10\ kpc}}$',
size=42)
# Axis limits
ax1.set_xlim(11.15, 12.29)
ax1.set_ylim(0.01, 0.79)
# ---------------------------------------------------------------------------
# Histogram
#
n, bins, patches=ax2.hist(gamaM1['lum_100'] - gamaM1['lum_10'],
bins=30, range=[0.05, 0.8], edgecolor='none',
orientation='horizontal', histtype='stepfilled',
color=BLUE0, alpha=0.80, normed=1)
n, bins, patches=ax2.hist(bcgM1['lum_100'] - bcgM1['lum_10'],
bins=20, range=[0.05, 0.8], edgecolor='none',
orientation='horizontal', histtype='stepfilled',
color=ORG4(0.6), alpha=0.50, normed=1, linewidth=4.0)
#
"""
n, bins, patches=ax2.hist(gamaM2['lum_100'] - gamaM2['lum_10'],
bins=30, range=[0.0, 0.7], linewidth=4.0,
orientation='horizontal', histtype='step',
color=BLUE0, alpha=1.0, normed=1)
n, bins, patches=ax2.hist(bcgM2['lum_100'] - bcgM2['lum_10'],
bins=15, range=[0.05, 0.7], linewidth=4.0,
orientation='horizontal', histtype='step',
color=RED0, alpha=0.90, normed=1)
"""
ax2.set_ylim(ax1.get_ylim())
ax2.axhline(0.0, linewidth=4.0, linestyle='-', c='k', alpha=0.2)
# Axes setup
# Minor Ticks on
ax2.minorticks_on()
ax2.tick_params(axis='x', which='minor', bottom='off', top='off')
# Axes Thickness
for axis in ['top','bottom','left','right']:
ax2.spines[axis].set_linewidth(3.5)
# Tick Label Size
for tick in ax2.xaxis.get_major_ticks():
tick.label.set_fontsize(24)
for tick in ax2.yaxis.get_major_ticks():
tick.label.set_fontsize(24)
# Tick Length and Width
ax2.tick_params('both', length=10, width=3.0, which='major')
ax2.tick_params('both', length=6, width=2.5, which='minor')
ax1.axhline(0.0, linewidth=4.5, linestyle='-', c='k', alpha=0.2)
ax2.yaxis.set_major_formatter(NullFormatter())
ax2.xaxis.set_major_formatter(NullFormatter())
plt.show()
fig.savefig('../figure/hscMassive_mtot_m100_10.png', dpi=300)
In [58]:
fig = plt.figure(figsize=(13, 9))
ax1 = plt.axes(recScat)
ax2 = plt.axes(recHist)
# SBP v.s. (cModel - SBP)
# ---------------------------------------------------------------------------
# Scatter plot
#ax1.axhline(0.0, linewidth=4.0, linestyle='-', c='k', alpha=0.2)
ax1.axvline(11.5, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)
ax1.axvline(11.7, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)
ax1.axvline(11.9, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)
# Matched ones
p1 = ax1.scatter(gamaClean['m100_c'],
gamaClean['lum_100'] - gamaClean['lum_10'], s=35.0,
alpha=0.20, facecolor=BLUE0, edgecolor='none',
label='$\Lambda \leq 20\ \mathrm{Central}$')
p2 = ax1.scatter(memUse['m100_c'],
memUse['lum_100'] - memUse['lum_10'], edgecolor='none',
s=((memUse['z_use'] - 0.10) * 600.0), cmap=ORG4, alpha=0.90,
c=toColorArr(memUse['LAMBDA_CLUSTER'], bottom=20.0, top=70.0),
label='$\Lambda > 20\ \mathrm{Central}$', marker='s')
"""
p3 = ax1.scatter(memM1['m100_c'],
memM1['lum_100'] - memM1['lum_10'], edgecolor='k',
s=((memM1['z_use'] - 0.10) * 500.0), alpha=0.95,
facecolor='none', label=None, marker='s', linewidth=1.5)
p4 = ax1.scatter(memM2b['m100_c'],
memM2b['lum_100'] - memM2b['lum_10'], edgecolor='k',
s=((memM2['z_use'] - 0.10) * 500.0), alpha=0.95,
facecolor='none', label=None, marker='s', linewidth=1.5)
"""
# M1
ax1.errorbar(gM1_mass[2], gM1_mdif1[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(mM1_mass[2], mM1_mdif1[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)
ax1.scatter(gM1_mass[2], gM1_mdif1[2], marker='^', s=400, facecolor=BLUE1,
edgecolor='k', linewidth=3.0, zorder=102, alpha=0.9,
label='$\mathrm{[11.5,11.7]}\ \Lambda \leq 20$')
ax1.scatter(mM1_mass[2], mM1_mdif1[2], marker='p', s=420, facecolor=RED1,
edgecolor='k', linewidth=3.0, zorder=102,
label='$\mathrm{[11.5,11.7]}\ \Lambda > 30$')
# M2
ax1.errorbar(gM2_mass[2], gM2_mdif1[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(mM2_mass[2], mM2_mdif1[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(gM2_mass[2], gM2_mdif1[2], marker='h', s=420, facecolor=BLUE1,
edgecolor='k', linewidth=3.0, zorder=102, alpha=0.9,
label='$\mathrm{[11.7,11.9]}\ \Lambda \leq 20$')
ax1.scatter(mM2_mass[2], mM2_mdif1[2], marker='8', s=420, facecolor=RED1,
edgecolor='k', linewidth=3.0, zorder=102,
label='$\mathrm{[11.7,11.9]}\ \Lambda > 30$')
# 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[1].set_color(ORG4(0.8))
legend.legendHandles[0].set_sizes([150])
legend.legendHandles[1].set_sizes([200])
#ax1.text(0.05, 0.04, '$\mathrm{Size:}\ {\Lambda}_{\mathrm{redMapper}}$',
# verticalalignment='bottom', horizontalalignment='left',
# fontsize=26.0, transform=ax1.transAxes, color=RED0)
# 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})\ (100\ \mathrm{Kpc})$', size=40)
ax1.set_ylabel('$\Delta(\log M{\star})_{\mathrm{100\ kpc}-\mathrm{10\ kpc}}$',
size=42)
# Axis limits
ax1.set_xlim(11.15, 12.29)
ax1.set_ylim(0.01, 0.79)
# ---------------------------------------------------------------------------
# Histogram
#
n, bins, patches=ax2.hist(gamaM1['lum_100'] - gamaM1['lum_10'],
bins=30, range=[0.05, 0.8], edgecolor='none',
orientation='horizontal', histtype='stepfilled',
color=BLUE0, alpha=0.80, normed=1)
n, bins, patches=ax2.hist(memM1['lum_100'] - memM1['lum_10'],
bins=20, range=[0.05, 0.8], edgecolor='none',
orientation='horizontal', histtype='stepfilled',
color=ORG4(0.6), alpha=0.50, normed=1, linewidth=4.0)
#
"""
n, bins, patches=ax2.hist(gamaM2['lum_100'] - gamaM2['lum_10'],
bins=30, range=[0.0, 0.7], linewidth=4.0,
orientation='horizontal', histtype='step',
color=BLUE0, alpha=1.0, normed=1)
n, bins, patches=ax2.hist(memM2['lum_100'] - memM2['lum_10'],
bins=15, range=[0.05, 0.7], linewidth=4.0,
orientation='horizontal', histtype='step',
color=RED0, alpha=0.90, normed=1)
"""
ax2.set_ylim(ax1.get_ylim())
ax2.axhline(0.0, linewidth=4.0, linestyle='-', c='k', alpha=0.2)
# Axes setup
# Minor Ticks on
ax2.minorticks_on()
ax2.tick_params(axis='x', which='minor', bottom='off', top='off')
# Axes Thickness
for axis in ['top','bottom','left','right']:
ax2.spines[axis].set_linewidth(3.5)
# Tick Label Size
for tick in ax2.xaxis.get_major_ticks():
tick.label.set_fontsize(24)
for tick in ax2.yaxis.get_major_ticks():
tick.label.set_fontsize(24)
# Tick Length and Width
ax2.tick_params('both', length=10, width=3.0, which='major')
ax2.tick_params('both', length=6, width=2.5, which='minor')
ax1.axhline(0.0, linewidth=4.5, linestyle='-', c='k', alpha=0.2)
ax2.yaxis.set_major_formatter(NullFormatter())
ax2.xaxis.set_major_formatter(NullFormatter())
plt.show()
#fig.savefig('../figure/hscMassive_mtot_m100_10_mem.png', dpi=300)
In [24]:
fig = plt.figure(figsize=(13, 9))
ax1 = plt.axes(recScat)
ax2 = plt.axes(recHist)
# SBP v.s. (cModel - SBP)
# ---------------------------------------------------------------------------
# Scatter plot
#ax1.axhline(0.0, linewidth=4.0, linestyle='-', c='k', alpha=0.2)
ax1.axvline(11.5, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)
ax1.axvline(11.7, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)
ax1.axvline(11.9, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)
# Matched ones
p1 = ax1.scatter(gamaClean['m100_c'],
gamaClean['lum_75'] - gamaClean['lum_10'], s=35.0,
alpha=0.20, facecolor=BLUE0, edgecolor='none',
label='$\Lambda \leq 20\ \mathrm{Central}$')
p2 = ax1.scatter(bcgUse['m100_c'],
bcgUse['lum_75'] - bcgUse['lum_10'], edgecolor='none',
s=((bcgUse['z_use'] - 0.10) * 600.0), cmap=ORG4, alpha=0.90,
c=toColorArr(bcgUse['LAMBDA_CLUSTER'], bottom=20.0, top=70.0),
label='$\Lambda > 20\ \mathrm{Central}$', marker='s')
"""
p3 = ax1.scatter(bcgM1['m100_c'],
bcgM1['lum_100'] - bcgM1['lum_10'], edgecolor='k',
s=((bcgM1['z_use'] - 0.10) * 500.0), alpha=0.95,
facecolor='none', label=None, marker='s', linewidth=1.5)
p4 = ax1.scatter(bcgM2b['m100_c'],
bcgM2b['lum_100'] - bcgM2b['lum_10'], edgecolor='k',
s=((bcgM2['z_use'] - 0.10) * 500.0), alpha=0.95,
facecolor='none', label=None, marker='s', linewidth=1.5)
"""
# M1
ax1.errorbar(gM1_mass[2], gM1_mdif2[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(bM1_mass[2], bM1_mdif2[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)
ax1.scatter(gM1_mass[2], gM1_mdif2[2], marker='^', s=400, facecolor=BLUE1,
edgecolor='k', linewidth=3.0, zorder=102, alpha=0.9,
label='$\mathrm{[11.5,11.7]}\ \Lambda \leq 20$')
ax1.scatter(bM1_mass[2], bM1_mdif2[2], marker='p', s=420, facecolor=RED1,
edgecolor='k', linewidth=3.0, zorder=102,
label='$\mathrm{[11.5,11.7]}\ \Lambda > 30$')
# M2
ax1.errorbar(gM2_mass[2], gM2_mdif2[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(bM2_mass[2], bM2_mdif2[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(gM2_mass[2], gM2_mdif2[2], marker='h', s=420, facecolor=BLUE1,
edgecolor='k', linewidth=3.0, zorder=102, alpha=0.9,
label='$\mathrm{[11.7,11.9]}\ \Lambda \leq 20$')
ax1.scatter(bM2_mass[2], bM2_mdif2[2], marker='8', s=420, facecolor=RED1,
edgecolor='k', linewidth=3.0, zorder=102,
label='$\mathrm{[11.7,11.9]}\ \Lambda > 30$')
# 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[1].set_color(ORG4(0.8))
legend.legendHandles[0].set_sizes([150])
legend.legendHandles[1].set_sizes([200])
#ax1.text(0.05, 0.04, '$\mathrm{Size:}\ {\Lambda}_{\mathrm{redMapper}}$',
# verticalalignment='bottom', horizontalalignment='left',
# fontsize=26.0, transform=ax1.transAxes, color=RED0)
# 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})\ (100\ \mathrm{Kpc})$', size=40)
ax1.set_ylabel('$\Delta(\log M{\star})_{\mathrm{75\ kpc}-\mathrm{10\ kpc}}$',
size=42)
# Axis limits
ax1.set_xlim(11.15, 12.29)
ax1.set_ylim(0.01, 0.79)
# ---------------------------------------------------------------------------
# Histogram
#
n, bins, patches=ax2.hist(gamaM1['lum_75'] - gamaM1['lum_10'],
bins=30, range=[0.05, 0.8], edgecolor='none',
orientation='horizontal', histtype='stepfilled',
color=BLUE0, alpha=0.80, normed=1)
n, bins, patches=ax2.hist(bcgM1['lum_75'] - bcgM1['lum_10'],
bins=20, range=[0.05, 0.8], edgecolor='none',
orientation='horizontal', histtype='stepfilled',
color=ORG4(0.6), alpha=0.50, normed=1, linewidth=4.0)
#
"""
n, bins, patches=ax2.hist(gamaM2['lum_75'] - gamaM2['lum_10'],
bins=30, range=[0.0, 0.7], linewidth=4.0,
orientation='horizontal', histtype='step',
color=BLUE0, alpha=1.0, normed=1)
n, bins, patches=ax2.hist(bcgM2['lum_75'] - bcgM2['lum_10'],
bins=15, range=[0.05, 0.7], linewidth=4.0,
orientation='horizontal', histtype='step',
color=RED0, alpha=0.90, normed=1)
"""
ax2.set_ylim(ax1.get_ylim())
ax2.axhline(0.0, linewidth=4.0, linestyle='-', c='k', alpha=0.2)
# Axes setup
# Minor Ticks on
ax2.minorticks_on()
ax2.tick_params(axis='x', which='minor', bottom='off', top='off')
# Axes Thickness
for axis in ['top','bottom','left','right']:
ax2.spines[axis].set_linewidth(3.5)
# Tick Label Size
for tick in ax2.xaxis.get_major_ticks():
tick.label.set_fontsize(24)
for tick in ax2.yaxis.get_major_ticks():
tick.label.set_fontsize(24)
# Tick Length and Width
ax2.tick_params('both', length=10, width=3.0, which='major')
ax2.tick_params('both', length=6, width=2.5, which='minor')
ax1.axhline(0.0, linewidth=4.5, linestyle='-', c='k', alpha=0.2)
ax2.yaxis.set_major_formatter(NullFormatter())
ax2.xaxis.set_major_formatter(NullFormatter())
plt.show()
fig.savefig('../figure/hscMassive_mtot_m75_10.png', dpi=230)
In [25]:
fig = plt.figure(figsize=(13, 9))
ax1 = plt.axes(recScat)
ax2 = plt.axes(recHist)
# SBP v.s. (cModel - SBP)
# ---------------------------------------------------------------------------
# Scatter plot
#ax1.axhline(0.0, linewidth=4.0, linestyle='-', c='k', alpha=0.2)
ax1.axvline(11.5, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)
ax1.axvline(11.7, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)
ax1.axvline(11.9, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)
# Matched ones
p1 = ax1.scatter(gamaClean['m100_c'],
gamaClean['lum_50'] - gamaClean['lum_10'], s=35.0,
alpha=0.20, facecolor=BLUE0, edgecolor='none',
label='$\Lambda \leq 20\ \mathrm{Central}$')
p2 = ax1.scatter(bcgUse['m100_c'],
bcgUse['lum_50'] - bcgUse['lum_10'], edgecolor='none',
s=((bcgUse['z_use'] - 0.10) * 600.0), cmap=ORG4, alpha=0.90,
c=toColorArr(bcgUse['LAMBDA_CLUSTER'], bottom=20.0, top=70.0),
label='$\Lambda > 20\ \mathrm{Central}$', marker='s')
"""
p3 = ax1.scatter(bcgM1['m100_c'],
bcgM1['lum_50'] - bcgM1['lum_10'], edgecolor='k',
s=((bcgM1['z_use'] - 0.10) * 500.0), alpha=0.95,
facecolor='none', label=None, marker='s', linewidth=1.5)
p4 = ax1.scatter(bcgM2b['m100_c'],
bcgM2b['lum_50'] - bcgM2b['lum_10'], edgecolor='k',
s=((bcgM2['z_use'] - 0.10) * 500.0), alpha=0.95,
facecolor='none', label=None, marker='s', linewidth=1.5)
"""
# M1
ax1.errorbar(gM1_mass[2], gM1_mdif3[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(bM1_mass[2], bM1_mdif3[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)
ax1.scatter(gM1_mass[2], gM1_mdif3[2], marker='^', s=400, facecolor=BLUE1,
edgecolor='k', linewidth=3.0, zorder=102, alpha=0.9,
label='$\mathrm{[11.5,11.7]}\ \Lambda \leq 20$')
ax1.scatter(bM1_mass[2], bM1_mdif3[2], marker='p', s=420, facecolor=RED1,
edgecolor='k', linewidth=3.0, zorder=102,
label='$\mathrm{[11.5,11.7]}\ \Lambda > 30$')
# M2
ax1.errorbar(gM2_mass[2], gM2_mdif3[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(bM2_mass[2], bM2_mdif3[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(gM2_mass[2], gM2_mdif3[2], marker='h', s=420, facecolor=BLUE1,
edgecolor='k', linewidth=3.0, zorder=102, alpha=0.9,
label='$\mathrm{[11.7,11.9]}\ \Lambda \leq 20$')
ax1.scatter(bM2_mass[2], bM2_mdif3[2], marker='8', s=420, facecolor=RED1,
edgecolor='k', linewidth=3.0, zorder=102,
label='$\mathrm{[11.7,11.9]}\ \Lambda > 30$')
# 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[1].set_color(ORG4(0.8))
legend.legendHandles[0].set_sizes([150])
legend.legendHandles[1].set_sizes([200])
#ax1.text(0.05, 0.04, '$\mathrm{Size:}\ {\Lambda}_{\mathrm{redMapper}}$',
# verticalalignment='bottom', horizontalalignment='left',
# fontsize=26.0, transform=ax1.transAxes, color=RED0)
# 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})\ (100\ \mathrm{Kpc})$', size=40)
ax1.set_ylabel('$\Delta(\log M{\star})_{\mathrm{50\ kpc}-\mathrm{10\ kpc}}$',
size=42)
# Axis limits
ax1.set_xlim(11.15, 12.29)
ax1.set_ylim(0.01, 0.79)
# ---------------------------------------------------------------------------
# Histogram
#
n, bins, patches=ax2.hist(gamaM1['lum_50'] - gamaM1['lum_10'],
bins=30, range=[0.05, 0.8], edgecolor='none',
orientation='horizontal', histtype='stepfilled',
color=BLUE0, alpha=0.80, normed=1)
n, bins, patches=ax2.hist(bcgM1['lum_50'] - bcgM1['lum_10'],
bins=20, range=[0.05, 0.8], edgecolor='none',
orientation='horizontal', histtype='stepfilled',
color=ORG4(0.6), alpha=0.50, normed=1, linewidth=4.0)
#
"""
n, bins, patches=ax2.hist(gamaM2['lum_50'] - gamaM2['lum_10'],
bins=30, range=[0.0, 0.7], linewidth=4.0,
orientation='horizontal', histtype='step',
color=BLUE0, alpha=1.0, normed=1)
n, bins, patches=ax2.hist(bcgM2['lum_50'] - bcgM2['lum_10'],
bins=15, range=[0.05, 0.7], linewidth=4.0,
orientation='horizontal', histtype='step',
color=RED0, alpha=0.90, normed=1)
"""
ax2.set_ylim(ax1.get_ylim())
ax2.axhline(0.0, linewidth=4.0, linestyle='-', c='k', alpha=0.2)
# Axes setup
# Minor Ticks on
ax2.minorticks_on()
ax2.tick_params(axis='x', which='minor', bottom='off', top='off')
# Axes Thickness
for axis in ['top','bottom','left','right']:
ax2.spines[axis].set_linewidth(3.5)
# Tick Label Size
for tick in ax2.xaxis.get_major_ticks():
tick.label.set_fontsize(24)
for tick in ax2.yaxis.get_major_ticks():
tick.label.set_fontsize(24)
# Tick Length and Width
ax2.tick_params('both', length=10, width=3.0, which='major')
ax2.tick_params('both', length=6, width=2.5, which='minor')
ax1.axhline(0.0, linewidth=4.5, linestyle='-', c='k', alpha=0.2)
ax2.yaxis.set_major_formatter(NullFormatter())
ax2.xaxis.set_major_formatter(NullFormatter())
plt.show()
fig.savefig('../figure/hscMassive_mtot_m50_10.png', dpi=230)