In [1]:
% matplotlib inline
from __future__ import (division,
print_function)
import os
import sys
import copy
import fnmatch
import warnings
import collections
import numpy as np
import scipy
try:
from scipy.stats import scoreatpercentile
except:
scoreatpercentile = False
from scipy.interpolate import interp1d
import cPickle as pickle
# Astropy
from astropy.io import fits
from astropy import units as u
from astropy.stats import sigma_clip
from astropy.table import Table, Column
from astropy.utils.console import ProgressBar
# AstroML
from astroML.plotting import hist
from astroML.density_estimation import KNeighborsDensity
try:
from sklearn.neighbors import KernelDensity
use_sklearn_KDE = True
except:
import warnings
warnings.warn("KDE will be removed in astroML version 0.3. Please "
"upgrade to scikit-learn 0.14+ and use "
"sklearn.neighbors.KernelDensity.", DeprecationWarning)
from astroML.density_estimation import KDE
use_sklearn_KDE = False
from sklearn.neighbors import KDTree
from sklearn.neighbors import BallTree
# Matplotlib related
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.patches import Ellipse
from matplotlib.ticker import NullFormatter, MaxNLocator, FormatStrFormatter
from matplotlib.collections import PatchCollection
tickFormat = FormatStrFormatter('$\mathbf{%g}$')
# 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 Greys_9, OrRd_9, Blues_9, Purples_9, YlGn_9
BLK = Greys_9.mpl_colormap
ORG = OrRd_9.mpl_colormap
BLU = Blues_9.mpl_colormap
GRN = YlGn_9.mpl_colormap
PUR = Purples_9.mpl_colormap
# Personal tools
from hscUtils import songPlotSetup, removeIsNullCol
from hscUtils import confidence_interval, ma_confidence_interval_1d, confidence_interval_1d
## Constants
# SDSS pivot wavelength
sdss_u_pivot = 3551.0
sdss_g_pivot = 4686.0
sdss_r_pivot = 6165.0
sdss_i_pivot = 7481.0
sdss_z_pivot = 8931.0
# HSC pivot wavelength
hsc_g_pivot = 4782.2
hsc_r_pivot = 6101.7
hsc_i_pivot = 7648.0
hsc_z_pivot = 8883.0
hsc_y_pivot = 9750.8
hscFiltWave = np.asarray([hsc_g_pivot, hsc_r_pivot, hsc_i_pivot, hsc_z_pivot, hsc_y_pivot])
"""
Absolute magnitude of the Sun in HSC filters
Right now, just use the DES filters
"""
SUN_G = 5.08
SUN_R = 4.62
SUN_I = 4.52
SUN_Z = 4.52
SUN_Y = 4.51
# Solar stellar metallicity
Z_SUN = 0.02
# definitions for the axes
left, width = 0.12, 0.66
right = left + width
bottom, height = 0.14, 0.85
bottom_h = left_h = left + width + 0.02
recScat = [left, bottom, width, height]
recHist = [right, bottom, 0.21, height]
SBP1 = [0.13, 0.12, 0.865, 0.30]
SBP2 = [0.13, 0.42, 0.865, 0.54]
EC1 = [0.135, 0.066, 0.862, 0.30]
EC2 = [0.135, 0.366, 0.862, 0.30]
EC3 = [0.135, 0.666, 0.862, 0.30]
REC = [0.12, 0.11, 0.87, 0.87]
# Universal RSMA array
RSMA_COMMON = np.arange(0.4, 4.2, 0.01)
EMPTY = (RSMA_COMMON * np.nan)
# Color
BLUE0 = "#92c5de"
BLUE1 = "#0571b0"
RED0 = "#f4a582"
RED1 = "#ca0020"
PURPLE0 = '#af8dc3'
PURPLE1 = '#762a83'
BROWN0 = '#bf812d'
BROWN1 = '#543005'
GREEN0 = '#7fbf7b'
GREEN1 = '#1b7837'
# 3-sigma
SIGMA1 = 0.3173
SIGMA2 = 0.0455
SIGMA3 = 0.0027
In [3]:
def loadGalfitOutput(pklFile):
"""
Load a .pkl file for GALFIT model
Return the GalfitParser object
"""
if not os.path.isfile(pklFile):
raise Exception("XXX Can not find the .pkl file : %s") % pklFile
else:
return pickle.load(open(pklFile, 'rb'))
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(array, bottom=None, top=None):
"""
Convert a data array to "color array" (between 0 and 1).
Parameters:
bottom, top :
"""
data = copy.deepcopy(array)
if top is not None:
cMax = top
else:
cMax = np.nanmax(data)
if bottom is not None:
cMin = bottom
else:
cMin = np.nanmin(data)
colorArr = ((data - cMin) / (cMax - cMin))
colorArr[colorArr <= 0] = 0.0
print(np.nanmin(colorArr), np.nanmax(colorArr))
return colorArr
def toSizeArr(array, bottom=None, top=None, maxSize=40):
"""
Convert a data array to "size array".
Parameters:
bottom, top :
"""
data = copy.deepcopy(array)
if top is not None:
data[data >= top] = top
if bottom is not None:
data[data <= bottom] = bottom
return ((data - np.nanmin(data)) /
(np.nanmax(data) - np.nanmin(data))) * maxSize
def 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
In [7]:
def get1SerSize(cat, pklDir, zStr='z_use', idStr='index', pix=0.168,
chisqThr=4.0, nSerThr=8.0, extinctionStr='extinction_i'):
galUse = np.empty(len(cat), dtype=bool)
galUse[:] = False
galRe1Ser = []
galAr1Ser = []
galMag1Ser = []
galNs1Ser = []
chi1Ser = []
aic1Ser = []
bic1Ser = []
for ii, galaxy in enumerate(cat):
galId = galaxy[idStr]
pattern1Ser = '*' + str(galId) + '*1ser.pkl'
pattern2Ser = '*' + str(galId) + '*2ser.pkl'
pkl1SerFound = findProfile(pattern1Ser, pklDir)
pkl1SerFound = findProfile(pattern1Ser, pklDir)
if len(pkl1SerFound) == 1:
pkl1SerFile = pkl1SerFound[0]
obj1Ser = loadGalfitOutput(pkl1SerFile)
if type(obj1Ser) is not numpy.core.records.recarray:
nser1Ser = obj1Ser.component_1.n
chisq1Ser = obj1Ser.reduced_chisq
try:
ai = galaxy[extinctionStr]
except Exception:
ai = 0.05
if (chisq1Ser <= chisqThr) and (nser1Ser < nSerThr):
scale = pixKpc(galaxy[zStr], show=False)
distmod = hUtil.cosmoDistMod(galaxy[zStr])
aic1, bic1, hq1 = gSimple.galfitAIC(obj1Ser)
chi1Ser.append(chisq1Ser)
aic1Ser.append(aic1)
bic1Ser.append(bic1)
galRe1Ser.append(obj1Ser.component_1.re * scale)
galAr1Ser.append(obj1Ser.component_1.ar)
galMag1Ser.append(obj1Ser.component_1.mag - ai - distmod)
galNs1Ser.append(obj1Ser.component_1.n)
galUse[ii] = True
print("### Number of models : %d" % len(galRe1Ser))
return galUse, galRe1Ser, galAr1Ser, galMag1Ser, galNs1Ser, chi1Ser, bic1Ser
In [8]:
def get2SerSize(cat, pklDir, zStr='z', idStr='objid_dr12', pix=0.168,
chisqThr=4.0, nSerThr=8.0, extinctionStr='extinction_i',
factor=0.95, massStr='logms_pca'):
galUse = np.empty(len(cat), dtype=bool)
galUse[:] = False
galRe1Ser = []
galAr1Ser = []
galMag1Ser = []
galNs1Ser = []
chi1Ser = []
aic1Ser = []
bic1Ser = []
galRe2SerI = []
galAr2SerI = []
galMag2SerI = []
galNs2SerI = []
galRe2SerO = []
galAr2SerO = []
galMag2SerO = []
galNs2SerO = []
innerFrac = []
chi2Ser = []
aic2Ser = []
bic2Ser = []
logmI = []
logmO = []
for ii, galaxy in enumerate(cat):
galId = galaxy[idStr]
pattern1Ser = '*' + str(galId) + '*1ser.pkl'
pattern2Ser = '*' + str(galId) + '*2ser.pkl'
pkl1SerFound = findProfile(pattern1Ser, pklDir)
pkl2SerFound = findProfile(pattern2Ser, pklDir)
if (len(pkl1SerFound) == 1) and (len(pkl2SerFound) == 1):
pkl1SerFile = pkl1SerFound[0]
obj1Ser = loadGalfitOutput(pkl1SerFile)
pkl2SerFile = pkl2SerFound[0]
obj2Ser = loadGalfitOutput(pkl2SerFile)
if (type(obj1Ser) is not numpy.core.records.recarray) and (
type(obj2Ser) is not numpy.core.records.recarray):
nser2A = obj1Ser.component_1.n
nser2B = obj2Ser.component_2.n
chisq1 = obj1Ser.reduced_chisq
aic1, bic1, hq1 = gSimple.galfitAIC(obj1Ser)
chisq2 = obj1Ser.reduced_chisq
aic2, bic2, hq2 = gSimple.galfitAIC(obj2Ser)
try:
ai = galaxy[extinctionStr]
except Exception:
ai = 0.05
logmT = galaxy[massStr]
if (logmT >= 11.2) and (bic2 < bic1) and (aic2 < aic1) and (
nser2A <= nSerThr) and (nser2B <= nSerThr):
scale = pixKpc(galaxy[zStr], show=False)
distmod = hUtil.cosmoDistMod(galaxy[zStr])
chi1Ser.append(chisq1)
aic1Ser.append(aic1)
bic1Ser.append(bic1)
galRe1Ser.append(obj1Ser.component_1.re * np.sqrt(obj1Ser.component_1.ar) * scale)
galAr1Ser.append(obj1Ser.component_1.ar)
galMag1Ser.append(obj1Ser.component_1.mag - ai - distmod)
galNs1Ser.append(obj1Ser.component_1.n)
""" 2 Sersic part """
chi2Ser.append(chisq2)
aic2Ser.append(aic2)
bic2Ser.append(bic2)
if (obj2Ser.component_1.re <= obj2Ser.component_2.re):
compInner = obj2Ser.component_1
compOuter = obj2Ser.component_2
else:
compInner = obj2Ser.component_2
compOuter = obj2Ser.component_1
galRe2SerI.append(compInner.re * scale * np.sqrt(compInner.ar))
galAr2SerI.append(compInner.ar)
galMag2SerI.append(compInner.mag)
galNs2SerI.append(compInner.n)
galRe2SerO.append(compOuter.re * scale * np.sqrt(compOuter.ar))
galAr2SerO.append(compOuter.ar)
galMag2SerO.append(compOuter.mag)
galNs2SerO.append(compOuter.n)
fluxI = 10.0 ** ((27.0 - compInner.mag) / 2.5)
fluxO = 10.0 ** ((27.0 - compOuter.mag) / 2.5)
innerFrac.append(fluxI / (fluxI + fluxO))
logmI.append(logmT + np.log10(fluxI / (fluxI + fluxO)))
logmO.append(logmT + np.log10(fluxO / (fluxI + fluxO)))
galUse[ii] = True
print("### Number of models : %d" % len(galRe1Ser))
return galRe2SerI, galRe2SerO, logmI, logmO
Three major datasets:
/Users/songhuang/work/hscs/gama_massive/galfit/redbcg: xx galaxies (vetted)redbcg_mass_use_dom.fits/Users/songhuang/work/hscs/gama_massive/galfit/redmem: xxx galaxies (not vetted)redmapper_mem_hscmatch_mass_sbpsum_modA_muI1.fits/Users/songhuang/work/hscs/gama_massive/galfit/gama: xxx galaxies (not vetted)gama_massive_160107_sbpsum_modA_muI1.fitsGAMA1, GAMA2, GAMA3
In [9]:
newDir = '/Users/songhuang/work/hscs/gama_massive/galfit/'
try:
bcgTab
except NameError:
pass
else:
del bcgTab
try:
memTab
except NameError:
pass
else:
del memTab
try:
gamaTab
except NameError:
pass
else:
del gamaTab
# Folder for 3 datasets
bcgDir = os.path.join(newDir, 'redbcg')
memDir = os.path.join(newDir, 'redmem')
gama1Dir = os.path.join(newDir, 'gama1')
gama2Dir = os.path.join(newDir, 'gama2')
gama3Dir = os.path.join(newDir, 'gama3')
# Two summary catalogs
#bcgCat = os.path.join(bcgDir, 'redmapper_bcg_hscmatch_mass_use_sbpsum_modA_muI1.fits')
bcgCat = os.path.join(newDir, 'redbcg_mass_use_dom.fits')
memCat = os.path.join(newDir, 'redmapper_mem_hscmatch_mass_sbpsum_modA_muI1.fits')
gama1Cat = os.path.join(newDir, 'gama1_mass_sbpsum_modA_muI1.fits')
gama2Cat = os.path.join(newDir, 'gama2_mass_sbpsum_modA_muI1.fits')
gama3Cat = os.path.join(newDir, 'gama3_mass_sbpsum_modA_muI1.fits')
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(gama1Cat):
raise Exception("## Can not find catalog for GAMA galaxies : %s" % gama1Cat)
else:
gama1Tab = Table.read(gama1Cat, format='fits')
if not os.path.isfile(gama2Cat):
raise Exception("## Can not find catalog for GAMA galaxies : %s" % gama2Cat)
else:
gama2Tab = Table.read(gama2Cat, format='fits')
if not os.path.isfile(gama3Cat):
raise Exception("## Can not find catalog for GAMA galaxies : %s" % gama3Cat)
else:
gama3Tab = Table.read(gama3Cat, format='fits')
In [10]:
print("## Deal with %i galaxies in BCG sample" % len(bcgTab))
print("## Deal with %i galaxies in Cluster member sample" % len(memTab))
print("## Deal with %i galaxies in GAMA1 sample" % len(gama1Tab))
print("## Deal with %i galaxies in GAMA2 sample" % len(gama2Tab))
print("## Deal with %i galaxies in GAMA3 sample" % len(gama3Tab))
In [11]:
mc = np.asarray([10.8, 11.2, 11.4, 11.6, 11.8])
rc = np.asarray([0.49, 0.774, 0.877, 1.101, 1.264])
ms = np.asarray([10.6, 10.8, 11.0, 11.2, 11.4, 11.6])
rs = np.asarray([0.389, 0.485, 0.616, 0.759, 0.891, 1.112])
In [12]:
### There are something wrong with the Kcorrection in the catalog....Update the Kcorrection value!!
"""
Using Kcorrect = absmag_cmodel - ABSMAG_ISEDFIT
"""
# Get absolute magnitude
amagG_bcg = getLuminosity(bcgTab['gmag_cmodel'], bcgTab['Z'], extinction=bcgTab['a_g'])
amagR_bcg = getLuminosity(bcgTab['rmag_cmodel'], bcgTab['Z'], extinction=bcgTab['a_r'])
amagI_bcg = getLuminosity(bcgTab['imag_cmodel'], bcgTab['Z'], extinction=bcgTab['a_i'])
amagZ_bcg = getLuminosity(bcgTab['zmag_cmodel'], bcgTab['Z'], extinction=bcgTab['a_z'])
amagY_bcg = getLuminosity(bcgTab['ymag_cmodel'], bcgTab['Z'], extinction=bcgTab['a_y'])
# Model A
bcgTab['KCORRECT_G'] = (amagG_bcg - bcgTab['ABSMAG_G'])
bcgTab['KCORRECT_R'] = (amagR_bcg - bcgTab['ABSMAG_R'])
bcgTab['KCORRECT_I'] = (amagI_bcg - bcgTab['ABSMAG_I'])
bcgTab['KCORRECT_Z'] = (amagZ_bcg - bcgTab['ABSMAG_Z'])
bcgTab['KCORRECT_Y'] = (amagY_bcg - bcgTab['ABSMAG_Y'])
# Model B
bcgTab['KCORRECT_b_G'] = (amagG_bcg - bcgTab['ABSMAG_b_G'])
bcgTab['KCORRECT_b_R'] = (amagR_bcg - bcgTab['ABSMAG_b_R'])
bcgTab['KCORRECT_b_I'] = (amagI_bcg - bcgTab['ABSMAG_b_I'])
bcgTab['KCORRECT_b_Z'] = (amagZ_bcg - bcgTab['ABSMAG_b_Z'])
bcgTab['KCORRECT_b_Y'] = (amagY_bcg - bcgTab['ABSMAG_b_Y'])
# Model C
bcgTab['KCORRECT_c_G'] = (amagG_bcg - bcgTab['ABSMAG_c_G'])
bcgTab['KCORRECT_c_R'] = (amagR_bcg - bcgTab['ABSMAG_c_R'])
bcgTab['KCORRECT_c_I'] = (amagI_bcg - bcgTab['ABSMAG_c_I'])
bcgTab['KCORRECT_c_Z'] = (amagZ_bcg - bcgTab['ABSMAG_c_Z'])
bcgTab['KCORRECT_c_Y'] = (amagY_bcg - bcgTab['ABSMAG_c_Y'])
In [13]:
### There are something wrong with the Kcorrection in the catalog....Update the Kcorrection value!!
"""
Using Kcorrect = absmag_cmodel - ABSMAG_ISEDFIT
"""
# Get absolute magnitude
amagG_mem = getLuminosity(memTab['gmag_cmodel'], memTab['Z'], extinction=memTab['a_g'])
amagR_mem = getLuminosity(memTab['rmag_cmodel'], memTab['Z'], extinction=memTab['a_r'])
amagI_mem = getLuminosity(memTab['imag_cmodel'], memTab['Z'], extinction=memTab['a_i'])
amagZ_mem = getLuminosity(memTab['zmag_cmodel'], memTab['Z'], extinction=memTab['a_z'])
amagY_mem = getLuminosity(memTab['ymag_cmodel'], memTab['Z'], extinction=memTab['a_y'])
# Model A
memTab['KCORRECT_G'] = (amagG_mem - memTab['ABSMAG_G'])
memTab['KCORRECT_R'] = (amagR_mem - memTab['ABSMAG_R'])
memTab['KCORRECT_I'] = (amagI_mem - memTab['ABSMAG_I'])
memTab['KCORRECT_Z'] = (amagZ_mem - memTab['ABSMAG_Z'])
memTab['KCORRECT_Y'] = (amagY_mem - memTab['ABSMAG_Y'])
# Model B
memTab['KCORRECT_b_G'] = (amagG_mem - memTab['ABSMAG_b_G'])
memTab['KCORRECT_b_R'] = (amagR_mem - memTab['ABSMAG_b_R'])
memTab['KCORRECT_b_I'] = (amagI_mem - memTab['ABSMAG_b_I'])
memTab['KCORRECT_b_Z'] = (amagZ_mem - memTab['ABSMAG_b_Z'])
memTab['KCORRECT_b_Y'] = (amagY_mem - memTab['ABSMAG_b_Y'])
# Model C
memTab['KCORRECT_c_G'] = (amagG_mem - memTab['ABSMAG_c_G'])
memTab['KCORRECT_c_R'] = (amagR_mem - memTab['ABSMAG_c_R'])
memTab['KCORRECT_c_I'] = (amagI_mem - memTab['ABSMAG_c_I'])
memTab['KCORRECT_c_Z'] = (amagZ_mem - memTab['ABSMAG_c_Z'])
memTab['KCORRECT_c_Y'] = (amagY_mem - memTab['ABSMAG_c_Y'])
In [14]:
### There are something wrong with the Kcorrection in the catalog....Update the Kcorrection value!!
"""
Using Kcorrect = absmag_cmodel - ABSMAG_ISEDFIT
"""
# Get absolute magnitude
amagG_gama1 = getLuminosity(gama1Tab['gmag_cmodel'], gama1Tab['Z'], extinction=gama1Tab['a_g'])
amagR_gama1 = getLuminosity(gama1Tab['rmag_cmodel'], gama1Tab['Z'], extinction=gama1Tab['a_r'])
amagI_gama1 = getLuminosity(gama1Tab['imag_cmodel'], gama1Tab['Z'], extinction=gama1Tab['a_i'])
amagZ_gama1 = getLuminosity(gama1Tab['zmag_cmodel'], gama1Tab['Z'], extinction=gama1Tab['a_z'])
amagY_gama1 = getLuminosity(gama1Tab['ymag_cmodel'], gama1Tab['Z'], extinction=gama1Tab['a_y'])
# Model A
gama1Tab['KCORRECT_G'] = (amagG_gama1 - gama1Tab['ABSMAG_G'])
gama1Tab['KCORRECT_R'] = (amagR_gama1 - gama1Tab['ABSMAG_R'])
gama1Tab['KCORRECT_I'] = (amagI_gama1 - gama1Tab['ABSMAG_I'])
gama1Tab['KCORRECT_Z'] = (amagZ_gama1 - gama1Tab['ABSMAG_Z'])
gama1Tab['KCORRECT_Y'] = (amagY_gama1 - gama1Tab['ABSMAG_Y'])
# Model B
gama1Tab['KCORRECT_b_G'] = (amagG_gama1 - gama1Tab['ABSMAG_b_G'])
gama1Tab['KCORRECT_b_R'] = (amagR_gama1 - gama1Tab['ABSMAG_b_R'])
gama1Tab['KCORRECT_b_I'] = (amagI_gama1 - gama1Tab['ABSMAG_b_I'])
gama1Tab['KCORRECT_b_Z'] = (amagZ_gama1 - gama1Tab['ABSMAG_b_Z'])
gama1Tab['KCORRECT_b_Y'] = (amagY_gama1 - gama1Tab['ABSMAG_b_Y'])
# Model C
gama1Tab['KCORRECT_c_G'] = (amagG_gama1 - gama1Tab['ABSMAG_c_G'])
gama1Tab['KCORRECT_c_R'] = (amagR_gama1 - gama1Tab['ABSMAG_c_R'])
gama1Tab['KCORRECT_c_I'] = (amagI_gama1 - gama1Tab['ABSMAG_c_I'])
gama1Tab['KCORRECT_c_Z'] = (amagZ_gama1 - gama1Tab['ABSMAG_c_Z'])
gama1Tab['KCORRECT_c_Y'] = (amagY_gama1 - gama1Tab['ABSMAG_c_Y'])
In [15]:
### There are something wrong with the Kcorrection in the catalog....Update the Kcorrection value!!
"""
Using Kcorrect = absmag_cmodel - ABSMAG_ISEDFIT
"""
# Get absolute magnitude
amagG_gama2 = getLuminosity(gama2Tab['gmag_cmodel'], gama2Tab['Z'], extinction=gama2Tab['a_g'])
amagR_gama2 = getLuminosity(gama2Tab['rmag_cmodel'], gama2Tab['Z'], extinction=gama2Tab['a_r'])
amagI_gama2 = getLuminosity(gama2Tab['imag_cmodel'], gama2Tab['Z'], extinction=gama2Tab['a_i'])
amagZ_gama2 = getLuminosity(gama2Tab['zmag_cmodel'], gama2Tab['Z'], extinction=gama2Tab['a_z'])
amagY_gama2 = getLuminosity(gama2Tab['ymag_cmodel'], gama2Tab['Z'], extinction=gama2Tab['a_y'])
# Model A
gama2Tab['KCORRECT_G'] = (amagG_gama2 - gama2Tab['ABSMAG_G'])
gama2Tab['KCORRECT_R'] = (amagR_gama2 - gama2Tab['ABSMAG_R'])
gama2Tab['KCORRECT_I'] = (amagI_gama2 - gama2Tab['ABSMAG_I'])
gama2Tab['KCORRECT_Z'] = (amagZ_gama2 - gama2Tab['ABSMAG_Z'])
gama2Tab['KCORRECT_Y'] = (amagY_gama2 - gama2Tab['ABSMAG_Y'])
# Model B
gama2Tab['KCORRECT_b_G'] = (amagG_gama2 - gama2Tab['ABSMAG_b_G'])
gama2Tab['KCORRECT_b_R'] = (amagR_gama2 - gama2Tab['ABSMAG_b_R'])
gama2Tab['KCORRECT_b_I'] = (amagI_gama2 - gama2Tab['ABSMAG_b_I'])
gama2Tab['KCORRECT_b_Z'] = (amagZ_gama2 - gama2Tab['ABSMAG_b_Z'])
gama2Tab['KCORRECT_b_Y'] = (amagY_gama2 - gama2Tab['ABSMAG_b_Y'])
# Model C
gama2Tab['KCORRECT_c_G'] = (amagG_gama2 - gama2Tab['ABSMAG_c_G'])
gama2Tab['KCORRECT_c_R'] = (amagR_gama2 - gama2Tab['ABSMAG_c_R'])
gama2Tab['KCORRECT_c_I'] = (amagI_gama2 - gama2Tab['ABSMAG_c_I'])
gama2Tab['KCORRECT_c_Z'] = (amagZ_gama2 - gama2Tab['ABSMAG_c_Z'])
gama2Tab['KCORRECT_c_Y'] = (amagY_gama2 - gama2Tab['ABSMAG_c_Y'])
In [16]:
### There are something wrong with the Kcorrection in the catalog....Update the Kcorrection value!!
"""
Using Kcorrect = absmag_cmodel - ABSMAG_ISEDFIT
"""
# Get absolute magnitude
amagG_gama3 = getLuminosity(gama3Tab['gmag_cmodel'], gama3Tab['Z'], extinction=gama3Tab['a_g'])
amagR_gama3 = getLuminosity(gama3Tab['rmag_cmodel'], gama3Tab['Z'], extinction=gama3Tab['a_r'])
amagI_gama3 = getLuminosity(gama3Tab['imag_cmodel'], gama3Tab['Z'], extinction=gama3Tab['a_i'])
amagZ_gama3 = getLuminosity(gama3Tab['zmag_cmodel'], gama3Tab['Z'], extinction=gama3Tab['a_z'])
amagY_gama3 = getLuminosity(gama3Tab['ymag_cmodel'], gama3Tab['Z'], extinction=gama3Tab['a_y'])
# Model A
gama3Tab['KCORRECT_G'] = (amagG_gama3 - gama3Tab['ABSMAG_G'])
gama3Tab['KCORRECT_R'] = (amagR_gama3 - gama3Tab['ABSMAG_R'])
gama3Tab['KCORRECT_I'] = (amagI_gama3 - gama3Tab['ABSMAG_I'])
gama3Tab['KCORRECT_Z'] = (amagZ_gama3 - gama3Tab['ABSMAG_Z'])
gama3Tab['KCORRECT_Y'] = (amagY_gama3 - gama3Tab['ABSMAG_Y'])
# Model B
gama3Tab['KCORRECT_b_G'] = (amagG_gama3 - gama3Tab['ABSMAG_b_G'])
gama3Tab['KCORRECT_b_R'] = (amagR_gama3 - gama3Tab['ABSMAG_b_R'])
gama3Tab['KCORRECT_b_I'] = (amagI_gama3 - gama3Tab['ABSMAG_b_I'])
gama3Tab['KCORRECT_b_Z'] = (amagZ_gama3 - gama3Tab['ABSMAG_b_Z'])
gama3Tab['KCORRECT_b_Y'] = (amagY_gama3 - gama3Tab['ABSMAG_b_Y'])
# Model C
gama3Tab['KCORRECT_c_G'] = (amagG_gama3 - gama3Tab['ABSMAG_c_G'])
gama3Tab['KCORRECT_c_R'] = (amagR_gama3 - gama3Tab['ABSMAG_c_R'])
gama3Tab['KCORRECT_c_I'] = (amagI_gama3 - gama3Tab['ABSMAG_c_I'])
gama3Tab['KCORRECT_c_Z'] = (amagZ_gama3 - gama3Tab['ABSMAG_c_Z'])
gama3Tab['KCORRECT_c_Y'] = (amagY_gama3 - gama3Tab['ABSMAG_c_Y'])
In [17]:
## BCG
## Luminosity based on hscPipe cModel:
lumI_bcg = getLuminosity(bcgTab['imag_cmodel'], bcgTab['Z'], extinction=bcgTab['a_i'],
amag_sun=amag_sun_des_i)
lumG_bcg = getLuminosity(bcgTab['gmag_cmodel'], bcgTab['Z'], extinction=bcgTab['a_g'],
amag_sun=amag_sun_des_g)
lumR_bcg = getLuminosity(bcgTab['rmag_cmodel'], bcgTab['Z'], extinction=bcgTab['a_r'],
amag_sun=amag_sun_des_r)
lumZ_bcg = getLuminosity(bcgTab['zmag_cmodel'], bcgTab['Z'], extinction=bcgTab['a_z'],
amag_sun=amag_sun_des_z)
## K-corrected luminosity:
lumI_kA_bcg = lumI_bcg + (bcgTab['KCORRECT_I'] / 2.5)
lumI_kB_bcg = lumI_bcg + (bcgTab['KCORRECT_b_I'] / 2.5)
lumI_kC_bcg = lumI_bcg + (bcgTab['KCORRECT_c_I'] / 2.5)
## K-corrected cModel color:
### Model A
gr_kA_bcg = (bcgTab['ABSMAG_G'] - bcgTab['ABSMAG_R'])
gi_kA_bcg = (bcgTab['ABSMAG_G'] - bcgTab['ABSMAG_I'])
gz_kA_bcg = (bcgTab['ABSMAG_G'] - bcgTab['ABSMAG_Z'])
ri_kA_bcg = (bcgTab['ABSMAG_R'] - bcgTab['ABSMAG_I'])
### Model B
gr_kB_bcg = (bcgTab['ABSMAG_b_G'] - bcgTab['ABSMAG_b_R'])
gi_kB_bcg = (bcgTab['ABSMAG_b_G'] - bcgTab['ABSMAG_b_I'])
gz_kB_bcg = (bcgTab['ABSMAG_b_G'] - bcgTab['ABSMAG_b_Z'])
ri_kB_bcg = (bcgTab['ABSMAG_b_R'] - bcgTab['ABSMAG_b_I'])
### Model C
gr_kC_bcg = (bcgTab['ABSMAG_c_G'] - bcgTab['ABSMAG_c_R'])
gi_kC_bcg = (bcgTab['ABSMAG_c_G'] - bcgTab['ABSMAG_c_I'])
gz_kC_bcg = (bcgTab['ABSMAG_c_G'] - bcgTab['ABSMAG_c_Z'])
ri_kC_bcg = (bcgTab['ABSMAG_c_R'] - bcgTab['ABSMAG_c_I'])
## Stellar mass from iSEDFit
logm2lI_A_bcg = (bcgTab['MSTAR'] - lumI_bcg)
logm2lI_B_bcg = (bcgTab['MSTAR_b'] - lumI_bcg)
logm2lI_C_bcg = (bcgTab['MSTAR_c'] - lumI_bcg)
In [18]:
## Members
## Luminosity based on hscPipe cModel:
lumI_mem = getLuminosity(memTab['imag_cmodel'], memTab['Z'], extinction=memTab['a_i'],
amag_sun=amag_sun_des_i)
lumG_mem = getLuminosity(memTab['gmag_cmodel'], memTab['Z'], extinction=memTab['a_g'],
amag_sun=amag_sun_des_g)
lumR_mem = getLuminosity(memTab['rmag_cmodel'], memTab['Z'], extinction=memTab['a_r'],
amag_sun=amag_sun_des_r)
lumZ_mem = getLuminosity(memTab['zmag_cmodel'], memTab['Z'], extinction=memTab['a_z'],
amag_sun=amag_sun_des_z)
## K-corrected luminosity:
lumI_kA_mem = lumI_mem + (memTab['KCORRECT_I'] / 2.5)
lumI_kB_mem = lumI_mem + (memTab['KCORRECT_b_I'] / 2.5)
lumI_kC_mem = lumI_mem + (memTab['KCORRECT_c_I'] / 2.5)
## K-corrected cModel color:
### Model A
gr_kA_mem = (memTab['ABSMAG_G'] - memTab['ABSMAG_R'])
gi_kA_mem = (memTab['ABSMAG_G'] - memTab['ABSMAG_I'])
gz_kA_mem = (memTab['ABSMAG_G'] - memTab['ABSMAG_Z'])
ri_kA_mem = (memTab['ABSMAG_R'] - memTab['ABSMAG_I'])
### Model B
gr_kB_mem = (memTab['ABSMAG_b_G'] - memTab['ABSMAG_b_R'])
gi_kB_mem = (memTab['ABSMAG_b_G'] - memTab['ABSMAG_b_I'])
gz_kB_mem = (memTab['ABSMAG_b_G'] - memTab['ABSMAG_b_Z'])
ri_kB_mem = (memTab['ABSMAG_b_R'] - memTab['ABSMAG_b_I'])
### Model C
gr_kC_mem = (memTab['ABSMAG_c_G'] - memTab['ABSMAG_c_R'])
gi_kC_mem = (memTab['ABSMAG_c_G'] - memTab['ABSMAG_c_I'])
gz_kC_mem = (memTab['ABSMAG_c_G'] - memTab['ABSMAG_c_Z'])
ri_kC_mem = (memTab['ABSMAG_c_R'] - memTab['ABSMAG_c_I'])
## Stellar mass from iSEDFit
logm2lI_A_mem = (memTab['MSTAR'] - lumI_mem)
logm2lI_B_mem = (memTab['MSTAR_b'] - lumI_mem)
logm2lI_C_mem = (memTab['MSTAR_c'] - lumI_mem)
In [19]:
## GAMA1
## Luminosity based on hscPipe cModel:
lumI_gama1 = getLuminosity(gama1Tab['imag_cmodel'], gama1Tab['Z'], extinction=gama1Tab['a_i'],
amag_sun=amag_sun_des_i)
lumG_gama1 = getLuminosity(gama1Tab['gmag_cmodel'], gama1Tab['Z'], extinction=gama1Tab['a_g'],
amag_sun=amag_sun_des_g)
lumR_gama1 = getLuminosity(gama1Tab['rmag_cmodel'], gama1Tab['Z'], extinction=gama1Tab['a_r'],
amag_sun=amag_sun_des_r)
lumZ_gama1 = getLuminosity(gama1Tab['zmag_cmodel'], gama1Tab['Z'], extinction=gama1Tab['a_z'],
amag_sun=amag_sun_des_z)
## K-corrected luminosity:
lumI_kA_gama1 = lumI_gama1 + (gama1Tab['KCORRECT_I'] / 2.5)
lumI_kB_gama1 = lumI_gama1 + (gama1Tab['KCORRECT_b_I'] / 2.5)
lumI_kC_gama1 = lumI_gama1 + (gama1Tab['KCORRECT_c_I'] / 2.5)
## K-corrected cModel color:
### Model A
gr_kA_gama1 = (gama1Tab['ABSMAG_G'] - gama1Tab['ABSMAG_R'])
gi_kA_gama1 = (gama1Tab['ABSMAG_G'] - gama1Tab['ABSMAG_I'])
gz_kA_gama1 = (gama1Tab['ABSMAG_G'] - gama1Tab['ABSMAG_Z'])
ri_kA_gama1 = (gama1Tab['ABSMAG_R'] - gama1Tab['ABSMAG_I'])
### Model B
gr_kB_gama1 = (gama1Tab['ABSMAG_b_G'] - gama1Tab['ABSMAG_b_R'])
gi_kB_gama1 = (gama1Tab['ABSMAG_b_G'] - gama1Tab['ABSMAG_b_I'])
gz_kB_gama1 = (gama1Tab['ABSMAG_b_G'] - gama1Tab['ABSMAG_b_Z'])
ri_kB_gama1 = (gama1Tab['ABSMAG_b_R'] - gama1Tab['ABSMAG_b_I'])
### Model C
gr_kC_gama1 = (gama1Tab['ABSMAG_c_G'] - gama1Tab['ABSMAG_c_R'])
gi_kC_gama1 = (gama1Tab['ABSMAG_c_G'] - gama1Tab['ABSMAG_c_I'])
gz_kC_gama1 = (gama1Tab['ABSMAG_c_G'] - gama1Tab['ABSMAG_c_Z'])
ri_kC_gama1 = (gama1Tab['ABSMAG_c_R'] - gama1Tab['ABSMAG_c_I'])
## Stellar mass from iSEDFit
logm2lI_A_gama1 = (gama1Tab['MSTAR'] - lumI_gama1)
logm2lI_B_gama1 = (gama1Tab['MSTAR_b'] - lumI_gama1)
logm2lI_C_gama1 = (gama1Tab['MSTAR_c'] - lumI_gama1)
In [20]:
## GAMA2
## Luminosity based on hscPipe cModel:
lumI_gama2 = getLuminosity(gama2Tab['imag_cmodel'], gama2Tab['Z'], extinction=gama2Tab['a_i'],
amag_sun=amag_sun_des_i)
lumG_gama2 = getLuminosity(gama2Tab['gmag_cmodel'], gama2Tab['Z'], extinction=gama2Tab['a_g'],
amag_sun=amag_sun_des_g)
lumR_gama2 = getLuminosity(gama2Tab['rmag_cmodel'], gama2Tab['Z'], extinction=gama2Tab['a_r'],
amag_sun=amag_sun_des_r)
lumZ_gama2 = getLuminosity(gama2Tab['zmag_cmodel'], gama2Tab['Z'], extinction=gama2Tab['a_z'],
amag_sun=amag_sun_des_z)
## K-corrected luminosity:
lumI_kA_gama2 = lumI_gama2 + (gama2Tab['KCORRECT_I'] / 2.5)
lumI_kB_gama2 = lumI_gama2 + (gama2Tab['KCORRECT_b_I'] / 2.5)
lumI_kC_gama2 = lumI_gama2 + (gama2Tab['KCORRECT_c_I'] / 2.5)
## K-corrected cModel color:
### Model A
gr_kA_gama2 = (gama2Tab['ABSMAG_G'] - gama2Tab['ABSMAG_R'])
gi_kA_gama2 = (gama2Tab['ABSMAG_G'] - gama2Tab['ABSMAG_I'])
gz_kA_gama2 = (gama2Tab['ABSMAG_G'] - gama2Tab['ABSMAG_Z'])
ri_kA_gama2 = (gama2Tab['ABSMAG_R'] - gama2Tab['ABSMAG_I'])
### Model B
gr_kB_gama2 = (gama2Tab['ABSMAG_b_G'] - gama2Tab['ABSMAG_b_R'])
gi_kB_gama2 = (gama2Tab['ABSMAG_b_G'] - gama2Tab['ABSMAG_b_I'])
gz_kB_gama2 = (gama2Tab['ABSMAG_b_G'] - gama2Tab['ABSMAG_b_Z'])
ri_kB_gama2 = (gama2Tab['ABSMAG_b_R'] - gama2Tab['ABSMAG_b_I'])
### Model C
gr_kC_gama2 = (gama2Tab['ABSMAG_c_G'] - gama2Tab['ABSMAG_c_R'])
gi_kC_gama2 = (gama2Tab['ABSMAG_c_G'] - gama2Tab['ABSMAG_c_I'])
gz_kC_gama2 = (gama2Tab['ABSMAG_c_G'] - gama2Tab['ABSMAG_c_Z'])
ri_kC_gama2 = (gama2Tab['ABSMAG_c_R'] - gama2Tab['ABSMAG_c_I'])
## Stellar mass from iSEDFit
logm2lI_A_gama2 = (gama2Tab['MSTAR'] - lumI_gama2)
logm2lI_B_gama2 = (gama2Tab['MSTAR_b'] - lumI_gama2)
logm2lI_C_gama2 = (gama2Tab['MSTAR_c'] - lumI_gama2)
In [21]:
## GAMA3
## Luminosity based on hscPipe cModel:
lumI_gama3 = getLuminosity(gama3Tab['imag_cmodel'], gama3Tab['Z'], extinction=gama3Tab['a_i'],
amag_sun=amag_sun_des_i)
lumG_gama3 = getLuminosity(gama3Tab['gmag_cmodel'], gama3Tab['Z'], extinction=gama3Tab['a_g'],
amag_sun=amag_sun_des_g)
lumR_gama3 = getLuminosity(gama3Tab['rmag_cmodel'], gama3Tab['Z'], extinction=gama3Tab['a_r'],
amag_sun=amag_sun_des_r)
lumZ_gama3 = getLuminosity(gama3Tab['zmag_cmodel'], gama3Tab['Z'], extinction=gama3Tab['a_z'],
amag_sun=amag_sun_des_z)
## K-corrected luminosity:
lumI_kA_gama3 = lumI_gama3 + (gama3Tab['KCORRECT_I'] / 2.5)
lumI_kB_gama3 = lumI_gama3 + (gama3Tab['KCORRECT_b_I'] / 2.5)
lumI_kC_gama3 = lumI_gama3 + (gama3Tab['KCORRECT_c_I'] / 2.5)
## K-corrected cModel color:
### Model A
gr_kA_gama3 = (gama3Tab['ABSMAG_G'] - gama3Tab['ABSMAG_R'])
gi_kA_gama3 = (gama3Tab['ABSMAG_G'] - gama3Tab['ABSMAG_I'])
gz_kA_gama3 = (gama3Tab['ABSMAG_G'] - gama3Tab['ABSMAG_Z'])
ri_kA_gama3 = (gama3Tab['ABSMAG_R'] - gama3Tab['ABSMAG_I'])
### Model B
gr_kB_gama3 = (gama3Tab['ABSMAG_b_G'] - gama3Tab['ABSMAG_b_R'])
gi_kB_gama3 = (gama3Tab['ABSMAG_b_G'] - gama3Tab['ABSMAG_b_I'])
gz_kB_gama3 = (gama3Tab['ABSMAG_b_G'] - gama3Tab['ABSMAG_b_Z'])
ri_kB_gama3 = (gama3Tab['ABSMAG_b_R'] - gama3Tab['ABSMAG_b_I'])
### Model C
gr_kC_gama3 = (gama3Tab['ABSMAG_c_G'] - gama3Tab['ABSMAG_c_R'])
gi_kC_gama3 = (gama3Tab['ABSMAG_c_G'] - gama3Tab['ABSMAG_c_I'])
gz_kC_gama3 = (gama3Tab['ABSMAG_c_G'] - gama3Tab['ABSMAG_c_Z'])
ri_kC_gama3 = (gama3Tab['ABSMAG_c_R'] - gama3Tab['ABSMAG_c_I'])
## Stellar mass from iSEDFit
logm2lI_A_gama3 = (gama3Tab['MSTAR'] - lumI_gama3)
logm2lI_B_gama3 = (gama3Tab['MSTAR_b'] - lumI_gama3)
logm2lI_C_gama3 = (gama3Tab['MSTAR_c'] - lumI_gama3)
In [22]:
bcgTab.add_column(Column(name='lumI_cmodel', data=lumI_bcg))
bcgTab.add_column(Column(name='lumG_cmodel', data=lumG_bcg))
bcgTab.add_column(Column(name='lumR_cmodel', data=lumR_bcg))
bcgTab.add_column(Column(name='lumZ_cmodel', data=lumZ_bcg))
bcgTab.add_column(Column(name='lumI_kA_cmodel', data=lumI_kA_bcg))
bcgTab.add_column(Column(name='lumI_kB_cmodel', data=lumI_kB_bcg))
bcgTab.add_column(Column(name='lumI_kC_cmodel', data=lumI_kC_bcg))
bcgTab.add_column(Column(name='gr_kA', data=gr_kA_bcg))
bcgTab.add_column(Column(name='gi_kA', data=gi_kA_bcg))
bcgTab.add_column(Column(name='gz_kA', data=gz_kA_bcg))
bcgTab.add_column(Column(name='ri_kA', data=ri_kA_bcg))
bcgTab.add_column(Column(name='gr_kB', data=gr_kB_bcg))
bcgTab.add_column(Column(name='gi_kB', data=gi_kB_bcg))
bcgTab.add_column(Column(name='gz_kB', data=gz_kB_bcg))
bcgTab.add_column(Column(name='ri_kB', data=ri_kB_bcg))
bcgTab.add_column(Column(name='gr_kC', data=gr_kC_bcg))
bcgTab.add_column(Column(name='gi_kC', data=gi_kC_bcg))
bcgTab.add_column(Column(name='gz_kC', data=gz_kC_bcg))
bcgTab.add_column(Column(name='ri_kC', data=ri_kC_bcg))
bcgTab.add_column(Column(name='logm2lI_A', data=logm2lI_A_bcg))
bcgTab.add_column(Column(name='logm2lI_B', data=logm2lI_B_bcg))
bcgTab.add_column(Column(name='logm2lI_C', data=logm2lI_C_bcg))
In [23]:
memTab.add_column(Column(name='lumI_cmodel', data=lumI_mem))
memTab.add_column(Column(name='lumG_cmodel', data=lumG_mem))
memTab.add_column(Column(name='lumR_cmodel', data=lumR_mem))
memTab.add_column(Column(name='lumZ_cmodel', data=lumZ_mem))
memTab.add_column(Column(name='lumI_kA_cmodel', data=lumI_kA_mem))
memTab.add_column(Column(name='lumI_kB_cmodel', data=lumI_kB_mem))
memTab.add_column(Column(name='lumI_kC_cmodel', data=lumI_kC_mem))
memTab.add_column(Column(name='gr_kA', data=gr_kA_mem))
memTab.add_column(Column(name='gi_kA', data=gi_kA_mem))
memTab.add_column(Column(name='gz_kA', data=gz_kA_mem))
memTab.add_column(Column(name='ri_kA', data=ri_kA_mem))
memTab.add_column(Column(name='gr_kB', data=gr_kB_mem))
memTab.add_column(Column(name='gi_kB', data=gi_kB_mem))
memTab.add_column(Column(name='gz_kB', data=gz_kB_mem))
memTab.add_column(Column(name='ri_kB', data=ri_kB_mem))
memTab.add_column(Column(name='gr_kC', data=gr_kC_mem))
memTab.add_column(Column(name='gi_kC', data=gi_kC_mem))
memTab.add_column(Column(name='gz_kC', data=gz_kC_mem))
memTab.add_column(Column(name='ri_kC', data=ri_kC_mem))
memTab.add_column(Column(name='logm2lI_A', data=logm2lI_A_mem))
memTab.add_column(Column(name='logm2lI_B', data=logm2lI_B_mem))
memTab.add_column(Column(name='logm2lI_C', data=logm2lI_C_mem))
In [24]:
gama1Tab.add_column(Column(name='lumI_cmodel', data=lumI_gama1))
gama1Tab.add_column(Column(name='lumG_cmodel', data=lumG_gama1))
gama1Tab.add_column(Column(name='lumR_cmodel', data=lumR_gama1))
gama1Tab.add_column(Column(name='lumZ_cmodel', data=lumZ_gama1))
gama1Tab.add_column(Column(name='lumI_kA_cmodel', data=lumI_kA_gama1))
gama1Tab.add_column(Column(name='lumI_kB_cmodel', data=lumI_kB_gama1))
gama1Tab.add_column(Column(name='lumI_kC_cmodel', data=lumI_kC_gama1))
gama1Tab.add_column(Column(name='gr_kA', data=gr_kA_gama1))
gama1Tab.add_column(Column(name='gi_kA', data=gi_kA_gama1))
gama1Tab.add_column(Column(name='gz_kA', data=gz_kA_gama1))
gama1Tab.add_column(Column(name='ri_kA', data=ri_kA_gama1))
gama1Tab.add_column(Column(name='gr_kB', data=gr_kB_gama1))
gama1Tab.add_column(Column(name='gi_kB', data=gi_kB_gama1))
gama1Tab.add_column(Column(name='gz_kB', data=gz_kB_gama1))
gama1Tab.add_column(Column(name='ri_kB', data=ri_kB_gama1))
gama1Tab.add_column(Column(name='gr_kC', data=gr_kC_gama1))
gama1Tab.add_column(Column(name='gi_kC', data=gi_kC_gama1))
gama1Tab.add_column(Column(name='gz_kC', data=gz_kC_gama1))
gama1Tab.add_column(Column(name='ri_kC', data=ri_kC_gama1))
gama1Tab.add_column(Column(name='logm2lI_A', data=logm2lI_A_gama1))
gama1Tab.add_column(Column(name='logm2lI_B', data=logm2lI_B_gama1))
gama1Tab.add_column(Column(name='logm2lI_C', data=logm2lI_C_gama1))
In [25]:
gama2Tab.add_column(Column(name='lumI_cmodel', data=lumI_gama2))
gama2Tab.add_column(Column(name='lumG_cmodel', data=lumG_gama2))
gama2Tab.add_column(Column(name='lumR_cmodel', data=lumR_gama2))
gama2Tab.add_column(Column(name='lumZ_cmodel', data=lumZ_gama2))
gama2Tab.add_column(Column(name='lumI_kA_cmodel', data=lumI_kA_gama2))
gama2Tab.add_column(Column(name='lumI_kB_cmodel', data=lumI_kB_gama2))
gama2Tab.add_column(Column(name='lumI_kC_cmodel', data=lumI_kC_gama2))
gama2Tab.add_column(Column(name='gr_kA', data=gr_kA_gama2))
gama2Tab.add_column(Column(name='gi_kA', data=gi_kA_gama2))
gama2Tab.add_column(Column(name='gz_kA', data=gz_kA_gama2))
gama2Tab.add_column(Column(name='ri_kA', data=ri_kA_gama2))
gama2Tab.add_column(Column(name='gr_kB', data=gr_kB_gama2))
gama2Tab.add_column(Column(name='gi_kB', data=gi_kB_gama2))
gama2Tab.add_column(Column(name='gz_kB', data=gz_kB_gama2))
gama2Tab.add_column(Column(name='ri_kB', data=ri_kB_gama2))
gama2Tab.add_column(Column(name='gr_kC', data=gr_kC_gama2))
gama2Tab.add_column(Column(name='gi_kC', data=gi_kC_gama2))
gama2Tab.add_column(Column(name='gz_kC', data=gz_kC_gama2))
gama2Tab.add_column(Column(name='ri_kC', data=ri_kC_gama2))
gama2Tab.add_column(Column(name='logm2lI_A', data=logm2lI_A_gama2))
gama2Tab.add_column(Column(name='logm2lI_B', data=logm2lI_B_gama2))
gama2Tab.add_column(Column(name='logm2lI_C', data=logm2lI_C_gama2))
In [26]:
gama3Tab.add_column(Column(name='lumI_cmodel', data=lumI_gama3))
gama3Tab.add_column(Column(name='lumG_cmodel', data=lumG_gama3))
gama3Tab.add_column(Column(name='lumR_cmodel', data=lumR_gama3))
gama3Tab.add_column(Column(name='lumZ_cmodel', data=lumZ_gama3))
gama3Tab.add_column(Column(name='lumI_kA_cmodel', data=lumI_kA_gama3))
gama3Tab.add_column(Column(name='lumI_kB_cmodel', data=lumI_kB_gama3))
gama3Tab.add_column(Column(name='lumI_kC_cmodel', data=lumI_kC_gama3))
gama3Tab.add_column(Column(name='gr_kA', data=gr_kA_gama3))
gama3Tab.add_column(Column(name='gi_kA', data=gi_kA_gama3))
gama3Tab.add_column(Column(name='gz_kA', data=gz_kA_gama3))
gama3Tab.add_column(Column(name='ri_kA', data=ri_kA_gama3))
gama3Tab.add_column(Column(name='gr_kB', data=gr_kB_gama3))
gama3Tab.add_column(Column(name='gi_kB', data=gi_kB_gama3))
gama3Tab.add_column(Column(name='gz_kB', data=gz_kB_gama3))
gama3Tab.add_column(Column(name='ri_kB', data=ri_kB_gama3))
gama3Tab.add_column(Column(name='gr_kC', data=gr_kC_gama3))
gama3Tab.add_column(Column(name='gi_kC', data=gi_kC_gama3))
gama3Tab.add_column(Column(name='gz_kC', data=gz_kC_gama3))
gama3Tab.add_column(Column(name='ri_kC', data=ri_kC_gama3))
gama3Tab.add_column(Column(name='logm2lI_A', data=logm2lI_A_gama3))
gama3Tab.add_column(Column(name='logm2lI_B', data=logm2lI_B_gama3))
gama3Tab.add_column(Column(name='logm2lI_C', data=logm2lI_C_gama3))
In [27]:
logmBcg, logrBcg, nserBcg = [], [], []
logm1Bcg, logr1Bcg, logm2Bcg, logr2Bcg = [], [], [], []
sample = 'redBCG'
sampleDir = bcgDir
idStr = 'ID_CLUSTER'
zStr = 'Z'
m2lStr = 'logm2lI_C'
extStr = 'a_i'
for gal in bcgTab:
galId = str(gal[idStr])
scale = pixKpc(gal[zStr], show=False)
logm2l = gal[m2lStr]
## Ser1
ser1File = sample + '_' + galId + '_HSC-I_full_1ser.pkl'
ser1File = os.path.join(newDir, bcgDir, ser1File)
if os.path.isfile(ser1File):
obj1Ser = loadGalfitOutput(ser1File)
try:
aic1, bic1, hq1 = gSimple.galfitAIC(obj1Ser)
except Exception:
print("## Something is wrong with: %s" % ser1File)
aic1, bic1, hq1 = None, None, None
else:
ser1comp = obj1Ser.component_1
# log(Re/kpc)
logr1 = np.log10(ser1comp.re * np.sqrt(ser1comp.ar) * scale)
# log(M)
logm1 = getLuminosity(ser1comp.mag, gal[zStr],
extinction=gal[extStr], amag_sun=amag_sun_des_i) + logm2l
if (ser1comp.good and ser1comp.n > 0.5 and ser1comp.n < 10.0 and
ser1comp.ar >= 0.3 and obj1Ser.reduced_chisq < 5.0 and
logr1 <= 2.6):
logmBcg.append(logm1)
logrBcg.append(logr1)
nserBcg.append(ser1comp.n)
else:
aic1, bic1, hq1 = None, None, None
#print("## Can not find 1 Sersic model for %s" % ser1File)
## Ser2
ser2File = sample + '_' + galId + '_HSC-I_full_2ser.pkl'
ser2File = os.path.join(newDir, sampleDir, ser2File)
if os.path.isfile(ser2File):
obj2Ser = loadGalfitOutput(ser2File)
try:
aic2, bic2, hq2 = gSimple.galfitAIC(obj2Ser)
except Exception:
aic2, bic2, hq2 = None, None, None
print("## Something is wrong with: %s" % ser2File)
else:
ser2comp1 = obj2Ser.component_1
ser2comp2 = obj2Ser.component_2
if ser2comp1.re > ser2comp2.re:
ser2comp1, ser2comp2 = ser2comp2, ser2comp1
flux1 = 10.0 ** ((27.0 - ser2comp1.mag) / 2.5)
flux2 = 10.0 ** ((27.0 - ser2comp2.mag) / 2.5)
frac1 = flux1 / (flux1 + flux2)
mu1 = (flux1 / (ser2comp1.ar * (ser2comp1.re ** 2.0 )))
mu2 = (flux2 / (ser2comp2.ar * (ser2comp2.re ** 2.0 )))
# log(Re/kpc)
logr2a = np.log10(ser2comp1.re * np.sqrt(ser2comp1.ar) * scale)
logr2b = np.log10(ser2comp2.re * np.sqrt(ser2comp2.ar) * scale)
# log(M)
logm2a = getLuminosity(ser2comp1.mag, gal[zStr],
extinction=gal[extStr], amag_sun=amag_sun_des_i) + logm2l
logm2b = getLuminosity(ser2comp2.mag, gal[zStr],
extinction=gal[extStr], amag_sun=amag_sun_des_i) + logm2l
if (ser2comp1.good and ser2comp1.n > 0.5 and ser2comp1.n < 4.5 and
ser2comp1.ar >= 0.25 and obj2Ser.reduced_chisq < 2.0 and
logr2a <= 2.0 and logr2a >= -1.0 and
ser2comp2.good and ser2comp2.n > 0.5 and ser2comp1.n < 4.5 and
ser2comp2.ar >= 0.25 and obj2Ser.reduced_chisq < 2.0 and
logr2b <= 2.5 and logr2b >= 0.0 and mu1 >= mu2):
logm1Bcg.append(logm2a)
logm2Bcg.append(logm2b)
logr1Bcg.append(logr2a)
logr2Bcg.append(logr2b)
else:
aic2, bic2, hq2 = None, None, None
#print("## Can not find 1 Sersic model for %s" % ser2File)
print(len(logmBcg))
print(len(logm1Bcg))
In [28]:
plt.plot(mc, rc, linestyle='-', linewidth=9.0, c='r', alpha=0.7)
plt.plot(ms, rs, linestyle='--', linewidth=9.0, c='b', alpha=0.8)
plt.scatter(logmBcg, logrBcg, c='k', alpha=0.45, s=90)
plt.scatter(logm1Bcg, logr1Bcg, c='b', alpha=0.45, s=90, marker='h')
plt.scatter(logm2Bcg, logr2Bcg, c='r', alpha=0.55, s=90, marker='8')
plt.xlim(10.0, 12.4)
plt.ylim(-0.4, 2.4)
plt.show()
In [29]:
temp1 = 10.0 ** np.asarray(logm1Bcg)
temp2 = 10.0 ** np.asarray(logm2Bcg)
plt.scatter(np.log10(temp1 + temp2), (temp2 / (temp1 + temp2)))
plt.show()
In [30]:
temp = plt.hist(nserBcg, 50, normed=1, facecolor='green', alpha=0.7)
plt.show()
In [31]:
logmMem, logrMem, nserMem = [], [], []
logm1Mem, logr1Mem, logm2Mem, logr2Mem = [], [], [], []
sample = 'redMem'
sampleDir = memDir
idStr = 'ISEDFIT_ID'
zStr = 'Z'
m2lStr = 'logm2lI_C'
extStr = 'a_i'
for gal in memTab:
galId = str(gal[idStr])
scale = pixKpc(gal[zStr], show=False)
logm2l = gal[m2lStr]
## Ser1
ser1File = sample + '_' + galId + '_HSC-I_full_1ser.pkl'
ser1File = os.path.join(newDir, memDir, ser1File)
if os.path.isfile(ser1File):
obj1Ser = loadGalfitOutput(ser1File)
try:
aic1, bic1, hq1 = gSimple.galfitAIC(obj1Ser)
except Exception:
print("## Something is wrong with: %s" % ser1File)
aic1, bic1, hq1 = None, None, None
else:
ser1comp = obj1Ser.component_1
# log(Re/kpc)
logr1 = np.log10(ser1comp.re * np.sqrt(ser1comp.ar) * scale)
# log(M)
logm1 = getLuminosity(ser1comp.mag, gal[zStr],
extinction=gal[extStr], amag_sun=amag_sun_des_i) + logm2l
if (ser1comp.good and ser1comp.n > 0.5 and ser1comp.n < 10.0 and
ser1comp.ar >= 0.3 and obj1Ser.reduced_chisq < 5.0 and
logr1 <= 2.6):
logmMem.append(logm1)
logrMem.append(logr1)
nserMem.append(nserMem)
else:
aic1, bic1, hq1 = None, None, None
#print("## Can not find 1 Sersic model for %s" % ser1File)
## Ser2
ser2File = sample + '_' + galId + '_HSC-I_full_2ser.pkl'
ser2File = os.path.join(newDir, sampleDir, ser2File)
if os.path.isfile(ser2File):
obj2Ser = loadGalfitOutput(ser2File)
try:
aic2, bic2, hq2 = gSimple.galfitAIC(obj2Ser)
except Exception:
aic2, bic2, hq2 = None, None, None
print("## Something is wrong with: %s" % ser2File)
else:
ser2comp1 = obj2Ser.component_1
ser2comp2 = obj2Ser.component_2
if ser2comp1.re > ser2comp2.re:
ser2comp1, ser2comp2 = ser2comp2, ser2comp1
flux1 = 10.0 ** ((27.0 - ser2comp1.mag) / 2.5)
flux2 = 10.0 ** ((27.0 - ser2comp2.mag) / 2.5)
frac1 = flux1 / (flux1 + flux2)
mu1 = (flux1 / (ser2comp1.ar * (ser2comp1.re ** 2.0 )))
mu2 = (flux2 / (ser2comp2.ar * (ser2comp2.re ** 2.0 )))
# log(Re/kpc)
logr2a = np.log10(ser2comp1.re * np.sqrt(ser2comp1.ar) * scale)
logr2b = np.log10(ser2comp2.re * np.sqrt(ser2comp2.ar) * scale)
# log(M)
logm2a = getLuminosity(ser2comp1.mag, gal[zStr],
extinction=gal[extStr], amag_sun=amag_sun_des_i) + logm2l
logm2b = getLuminosity(ser2comp2.mag, gal[zStr],
extinction=gal[extStr], amag_sun=amag_sun_des_i) + logm2l
if (ser2comp1.good and ser2comp1.n > 0.5 and ser2comp1.n < 4.5 and
ser2comp1.ar >= 0.25 and obj2Ser.reduced_chisq < 2.0 and
logr2a <= 2.0 and logr2a >= -1.0 and
ser2comp2.good and ser2comp2.n > 0.5 and ser2comp1.n < 4.5 and
ser2comp2.ar >= 0.25 and obj2Ser.reduced_chisq < 2.0 and
logr2b <= 2.5 and logr2b >= 0.0 and mu1 >= mu2):
logm1Mem.append(logm2a)
logm2Mem.append(logm2b)
logr1Mem.append(logr2a)
logr2Mem.append(logr2b)
else:
aic2, bic2, hq2 = None, None, None
#print("## Can not find 1 Sersic model for %s" % ser2File)
print(len(logmMem))
print(len(logm1Mem))
In [32]:
plt.xlim(10.0, 12.4)
plt.ylim(-0.4, 2.4)
plt.plot(mc, rc, linestyle='-', linewidth=9.0, c='r', alpha=0.7)
plt.plot(ms, rs, linestyle='--', linewidth=9.0, c='b', alpha=0.8)
plt.scatter(logmMem, logrMem, c='k', alpha=0.35, s=60)
plt.scatter(logm1Mem, logr1Mem, c='b', alpha=0.35, marker='h', s=40)
plt.scatter(logm2Mem, logr2Mem, c='r', alpha=0.35, marker='8', s=40)
plt.show()
In [33]:
logmGama1, logrGama1, nserGama1 = [], [], []
logm1Gama1, logr1Gama1, logm2Gama1, logr2Gama1 = [], [], [], []
sample = 'gama'
sampleDir = gama1Dir
idStr = 'ISEDFIT_ID'
zStr = 'Z'
m2lStr = 'logm2lI_C'
extStr = 'a_i'
for gal in gama1Tab:
galId = str(gal[idStr])
scale = pixKpc(gal[zStr], show=False)
logm2l = gal[m2lStr]
## Ser1
ser1File = sample + '_' + galId + '_HSC-I_full_1ser.pkl'
ser1File = os.path.join(newDir, gama1Dir, ser1File)
if os.path.isfile(ser1File):
obj1Ser = loadGalfitOutput(ser1File)
try:
aic1, bic1, hq1 = gSimple.galfitAIC(obj1Ser)
except Exception:
print("## Something is wrong with: %s" % ser1File)
aic1, bic1, hq1 = None, None, None
else:
ser1comp = obj1Ser.component_1
# log(Re/kpc)
logr1 = np.log10(ser1comp.re * np.sqrt(ser1comp.ar) * scale)
# log(M)
logm1 = getLuminosity(ser1comp.mag, gal[zStr],
extinction=gal[extStr], amag_sun=amag_sun_des_i) + logm2l
if (ser1comp.good and ser1comp.n > 0.5 and ser1comp.n < 10.0 and
ser1comp.ar >= 0.3 and obj1Ser.reduced_chisq < 5.0 and
logr1 <= 2.6):
logmGama1.append(logm1)
logrGama1.append(logr1)
nserGama1.append(ser1comp.n)
else:
aic1, bic1, hq1 = None, None, None
#print("## Can not find 1 Sersic model for %s" % ser1File)
## Ser2
ser2File = sample + '_' + galId + '_HSC-I_full_2ser.pkl'
ser2File = os.path.join(newDir, sampleDir, ser2File)
if os.path.isfile(ser2File):
obj2Ser = loadGalfitOutput(ser2File)
try:
aic2, bic2, hq2 = gSimple.galfitAIC(obj2Ser)
except Exception:
aic2, bic2, hq2 = None, None, None
print("## Something is wrong with: %s" % ser2File)
else:
ser2comp1 = obj2Ser.component_1
ser2comp2 = obj2Ser.component_2
if ser2comp1.re > ser2comp2.re:
ser2comp1, ser2comp2 = ser2comp2, ser2comp1
flux1 = 10.0 ** ((27.0 - ser2comp1.mag) / 2.5)
flux2 = 10.0 ** ((27.0 - ser2comp2.mag) / 2.5)
frac1 = flux1 / (flux1 + flux2)
mu1 = (flux1 / (ser2comp1.ar * (ser2comp1.re ** 2.0 )))
mu2 = (flux2 / (ser2comp2.ar * (ser2comp2.re ** 2.0 )))
# log(Re/kpc)
logr2a = np.log10(ser2comp1.re * np.sqrt(ser2comp1.ar) * scale)
logr2b = np.log10(ser2comp2.re * np.sqrt(ser2comp2.ar) * scale)
# log(M)
logm2a = getLuminosity(ser2comp1.mag, gal[zStr],
extinction=gal[extStr], amag_sun=amag_sun_des_i) + logm2l
logm2b = getLuminosity(ser2comp2.mag, gal[zStr],
extinction=gal[extStr], amag_sun=amag_sun_des_i) + logm2l
if (ser2comp1.good and ser2comp1.n > 0.5 and ser2comp1.n < 4.5 and
ser2comp1.ar >= 0.25 and obj2Ser.reduced_chisq < 2.0 and
logr2a <= 2.0 and logr2a >= -1.0 and
ser2comp2.good and ser2comp2.n > 0.5 and ser2comp1.n < 4.5 and
ser2comp2.ar >= 0.25 and obj2Ser.reduced_chisq < 2.0 and
logr2b <= 2.5 and logr2b >= 0.0 and mu1 >= mu2):
logm1Gama1.append(logm2a)
logm2Gama1.append(logm2b)
logr1Gama1.append(logr2a)
logr2Gama1.append(logr2b)
else:
aic2, bic2, hq2 = None, None, None
#print("## Can not find 1 Sersic model for %s" % ser2File)
print(len(logmGama1))
print(len(logm1Gama1))
In [34]:
plt.xlim(10.0, 12.4)
plt.ylim(-0.4, 2.4)
plt.plot(mc, rc, linestyle='-', linewidth=9.0, c='r', alpha=0.7)
plt.plot(ms, rs, linestyle='--', linewidth=9.0, c='b', alpha=0.8)
plt.scatter(logmGama1, logrGama1, c='k', alpha=0.25, s=40)
plt.scatter(logm1Gama1, logr1Gama1, c='b', alpha=0.20, marker='h', s=30)
plt.scatter(logm2Gama1, logr2Gama1, c='r', alpha=0.20, marker='8', s=30)
plt.show()
In [35]:
logmGama2, logrGama2, nserGama2 = [], [], []
logm1Gama2, logr1Gama2, logm2Gama2, logr2Gama2 = [], [], [], []
sample = 'gama'
sampleDir = gama2Dir
idStr = 'ISEDFIT_ID'
zStr = 'Z'
m2lStr = 'logm2lI_C'
extStr = 'a_i'
for gal in gama2Tab:
galId = str(gal[idStr])
scale = pixKpc(gal[zStr], show=False)
logm2l = gal[m2lStr]
## Ser1
ser1File = sample + '_' + galId + '_HSC-I_full_1ser.pkl'
ser1File = os.path.join(newDir, gama2Dir, ser1File)
if os.path.isfile(ser1File):
obj1Ser = loadGalfitOutput(ser1File)
try:
aic1, bic1, hq1 = gSimple.galfitAIC(obj1Ser)
except Exception:
print("## Something is wrong with: %s" % ser1File)
aic1, bic1, hq1 = None, None, None
else:
ser1comp = obj1Ser.component_1
# log(Re/kpc)
logr1 = np.log10(ser1comp.re * np.sqrt(ser1comp.ar) * scale)
# log(M)
logm1 = getLuminosity(ser1comp.mag, gal[zStr],
extinction=gal[extStr], amag_sun=amag_sun_des_i) + logm2l
if (ser1comp.good and ser1comp.n > 0.5 and ser1comp.n < 10.0 and
ser1comp.ar >= 0.3 and obj1Ser.reduced_chisq < 5.0 and
logr1 <= 2.6):
logmGama2.append(logm1)
logrGama2.append(logr1)
nserGama2.append(ser1comp.n)
else:
aic1, bic1, hq1 = None, None, None
#print("## Can not find 1 Sersic model for %s" % ser1File)
## Ser2
ser2File = sample + '_' + galId + '_HSC-I_full_2ser.pkl'
ser2File = os.path.join(newDir, sampleDir, ser2File)
if os.path.isfile(ser2File):
obj2Ser = loadGalfitOutput(ser2File)
try:
aic2, bic2, hq2 = gSimple.galfitAIC(obj2Ser)
except Exception:
aic2, bic2, hq2 = None, None, None
print("## Something is wrong with: %s" % ser2File)
else:
ser2comp1 = obj2Ser.component_1
ser2comp2 = obj2Ser.component_2
if ser2comp1.re > ser2comp2.re:
ser2comp1, ser2comp2 = ser2comp2, ser2comp1
flux1 = 10.0 ** ((27.0 - ser2comp1.mag) / 2.5)
flux2 = 10.0 ** ((27.0 - ser2comp2.mag) / 2.5)
frac1 = flux1 / (flux1 + flux2)
mu1 = (flux1 / (ser2comp1.ar * (ser2comp1.re ** 2.0 )))
mu2 = (flux2 / (ser2comp2.ar * (ser2comp2.re ** 2.0 )))
# log(Re/kpc)
logr2a = np.log10(ser2comp1.re * np.sqrt(ser2comp1.ar) * scale)
logr2b = np.log10(ser2comp2.re * np.sqrt(ser2comp2.ar) * scale)
# log(M)
logm2a = getLuminosity(ser2comp1.mag, gal[zStr],
extinction=gal[extStr], amag_sun=amag_sun_des_i) + logm2l
logm2b = getLuminosity(ser2comp2.mag, gal[zStr],
extinction=gal[extStr], amag_sun=amag_sun_des_i) + logm2l
if (ser2comp1.good and ser2comp1.n > 0.5 and ser2comp1.n < 4.5 and
ser2comp1.ar >= 0.25 and obj2Ser.reduced_chisq < 2.0 and
logr2a <= 2.0 and logr2a >= -1.0 and
ser2comp2.good and ser2comp2.n > 0.5 and ser2comp1.n < 4.5 and
ser2comp2.ar >= 0.25 and obj2Ser.reduced_chisq < 2.0 and
logr2b <= 2.5 and logr2b >= 0.0 and mu1 >= mu2):
logm1Gama2.append(logm2a)
logm2Gama2.append(logm2b)
logr1Gama2.append(logr2a)
logr2Gama2.append(logr2b)
else:
aic2, bic2, hq2 = None, None, None
#print("## Can not find 1 Sersic model for %s" % ser2File)
print(len(logmGama2))
print(len(logm1Gama2))
In [36]:
plt.xlim(10.0, 12.4)
plt.ylim(-0.4, 2.4)
plt.plot(mc, rc, linestyle='-', linewidth=9.0, c='r', alpha=0.7)
plt.plot(ms, rs, linestyle='--', linewidth=9.0, c='b', alpha=0.8)
plt.scatter(logmGama2, logrGama2, c='k', alpha=0.25, s=40)
plt.scatter(logm1Gama2, logr1Gama2, c='b', alpha=0.20, marker='h', s=30)
plt.scatter(logm2Gama2, logr2Gama2, c='r', alpha=0.20, marker='8', s=30)
plt.show()
In [37]:
logmGama3, logrGama3, nserGama3 = [], [], []
logm1Gama3, logr1Gama3, logm2Gama3, logr2Gama3 = [], [], [], []
sample = 'gama'
sampleDir = gama3Dir
idStr = 'ISEDFIT_ID'
zStr = 'Z'
m2lStr = 'logm2lI_C'
extStr = 'a_i'
for gal in gama3Tab:
galId = str(gal[idStr])
scale = pixKpc(gal[zStr], show=False)
logm2l = gal[m2lStr]
## Ser1
ser1File = sample + '_' + galId + '_HSC-I_full_1ser.pkl'
ser1File = os.path.join(newDir, gama3Dir, ser1File)
if os.path.isfile(ser1File):
obj1Ser = loadGalfitOutput(ser1File)
try:
aic1, bic1, hq1 = gSimple.galfitAIC(obj1Ser)
except Exception:
print("## Something is wrong with: %s" % ser1File)
aic1, bic1, hq1 = None, None, None
else:
ser1comp = obj1Ser.component_1
# log(Re/kpc)
logr1 = np.log10(ser1comp.re * np.sqrt(ser1comp.ar) * scale)
# log(M)
logm1 = getLuminosity(ser1comp.mag, gal[zStr],
extinction=gal[extStr], amag_sun=amag_sun_des_i) + logm2l
if (ser1comp.good and ser1comp.n > 0.5 and ser1comp.n < 10.0 and
ser1comp.ar >= 0.3 and obj1Ser.reduced_chisq < 5.0 and
logr1 <= 2.6):
logmGama3.append(logm1)
logrGama3.append(logr1)
nserGama3.append(ser1comp.n)
else:
aic1, bic1, hq1 = None, None, None
#print("## Can not find 1 Sersic model for %s" % ser1File)
## Ser2
ser2File = sample + '_' + galId + '_HSC-I_full_2ser.pkl'
ser2File = os.path.join(newDir, sampleDir, ser2File)
if os.path.isfile(ser2File):
obj2Ser = loadGalfitOutput(ser2File)
try:
aic2, bic2, hq2 = gSimple.galfitAIC(obj2Ser)
except Exception:
aic2, bic2, hq2 = None, None, None
print("## Something is wrong with: %s" % ser2File)
else:
ser2comp1 = obj2Ser.component_1
ser2comp2 = obj2Ser.component_2
if ser2comp1.re > ser2comp2.re:
ser2comp1, ser2comp2 = ser2comp2, ser2comp1
flux1 = 10.0 ** ((27.0 - ser2comp1.mag) / 2.5)
flux2 = 10.0 ** ((27.0 - ser2comp2.mag) / 2.5)
frac1 = flux1 / (flux1 + flux2)
mu1 = (flux1 / (ser2comp1.ar * (ser2comp1.re ** 2.0 )))
mu2 = (flux2 / (ser2comp2.ar * (ser2comp2.re ** 2.0 )))
# log(Re/kpc)
logr2a = np.log10(ser2comp1.re * np.sqrt(ser2comp1.ar) * scale)
logr2b = np.log10(ser2comp2.re * np.sqrt(ser2comp2.ar) * scale)
# log(M)
logm2a = getLuminosity(ser2comp1.mag, gal[zStr],
extinction=gal[extStr], amag_sun=amag_sun_des_i) + logm2l
logm2b = getLuminosity(ser2comp2.mag, gal[zStr],
extinction=gal[extStr], amag_sun=amag_sun_des_i) + logm2l
if (ser2comp1.good and ser2comp1.n > 0.5 and ser2comp1.n < 4.5 and
ser2comp1.ar >= 0.25 and obj2Ser.reduced_chisq < 2.0 and
logr2a <= 2.0 and logr2a >= -1.0 and
ser2comp2.good and ser2comp2.n > 0.5 and ser2comp1.n < 4.5 and
ser2comp2.ar >= 0.25 and obj2Ser.reduced_chisq < 2.0 and
logr2b <= 2.5 and logr2b >= 0.0 and mu1 >= mu2):
logm1Gama3.append(logm2a)
logm2Gama3.append(logm2b)
logr1Gama3.append(logr2a)
logr2Gama3.append(logr2b)
else:
aic2, bic2, hq2 = None, None, None
#print("## Can not find 1 Sersic model for %s" % ser2File)
print(len(logmGama3))
print(len(logm1Gama3))
In [38]:
plt.xlim(10.0, 12.4)
plt.ylim(-0.4, 2.4)
plt.plot(mc, rc, linestyle='-', linewidth=9.0, c='r', alpha=0.7)
plt.plot(ms, rs, linestyle='--', linewidth=9.0, c='b', alpha=0.8)
plt.scatter(logmGama3, logrGama3, c='k', alpha=0.25, s=40)
plt.scatter(logm1Gama3, logr1Gama3, c='b', alpha=0.20, marker='h', s=30)
plt.scatter(logm2Gama3, logr2Gama3, c='r', alpha=0.20, marker='8', s=30)
plt.show()
In [39]:
plt.xlim(10.5, 12.6)
plt.ylim(0.05, 2.59)
plt.xlabel('$\log\ (M_{\star}/M_{\odot})$', fontsize=35)
plt.ylabel('$\log\ (R_{\mathrm{e}}/\mathrm{Kpc})$', fontsize=35)
plt.scatter(logmGama1, logrGama1, c='k', alpha=0.05)
plt.scatter(logmGama2, logrGama2, c='k', alpha=0.05)
plt.scatter(logmGama3, logrGama3, c='k', alpha=0.05)
plt.plot(mc, rc, linestyle='-', linewidth=9.0, c='r', alpha=0.7)
plt.plot(ms, rs, linestyle='--', linewidth=9.0, c='b', alpha=0.8)
#plt.scatter(logmMem, logrMem, c='g', alpha=0.1, s=30, marker='s')
plt.scatter(logmBcg, logrBcg, c='r', alpha=0.5, s=70, marker='h')
plt.show()
In [40]:
logmGama = np.asarray(logmGama1 + logmGama2 + logmGama3).ravel()
logrGama = np.asarray(logrGama1 + logrGama2 + logrGama3).ravel()
nserGama = np.asarray(nserGama1 + nserGama2 + nserGama3).ravel()
In [41]:
logmGamaU = logmGama[logmGama >= 11.42]
logrGamaU = logrGama[logmGama >= 11.42]
nserGamaU = nserGama[logmGama >= 11.42]
In [42]:
logmBcg = np.asarray(logmBcg)
logrBcg = np.asarray(logrBcg)
nserBcg = np.asarray(nserBcg)
logmBcgU = logmBcg[logmBcg >= 11.42]
logrBcgU = logrBcg[logmBcg >= 11.42]
nserBcgU = nserBcg[logmBcg >= 11.42]
In [43]:
logmGamaM = logmGama[logmGama >= 11.65]
logrGamaM = logrGama[logmGama >= 11.65]
nserGamaM = nserGama[logmGama >= 11.65]
logmBcgM = logmBcg[logmBcg >= 11.65]
logrBcgM = logrBcg[logmBcg >= 11.65]
nserBcgM = nserBcg[logmBcg >= 11.65]
In [44]:
print(len(logmBcgM))
print(len(logmGamaM))
In [45]:
temp = plt.hist(nserGamaM, 30, normed=1, edgecolor='black',
alpha=1.0, histtype='step', linewidth=3.0)
temp = plt.hist(nserBcgM, 30, normed=1, facecolor='red',
alpha=0.7, histtype='stepfilled')
plt.show()
In [46]:
plt.xlabel('$\log\ (M_{\star}/M_{\odot})$', fontsize=35)
plt.ylabel('$\log\ (R_{\mathrm{e}}/\mathrm{Kpc})$', fontsize=35)
plt.scatter(logm1Gama2, logr1Gama2, c='b', alpha=0.05)
plt.scatter(logm2Gama2, logr2Gama2, c='r', alpha=0.05)
plt.scatter(logm1Gama1, logr1Gama1, c='b', alpha=0.03)
plt.scatter(logm2Gama1, logr2Gama1, c='r', alpha=0.03)
plt.plot(mc, rc, linestyle='-', linewidth=9.0, c='r', alpha=0.7)
plt.plot(ms, rs, linestyle='--', linewidth=9.0, c='b', alpha=0.8)
#plt.scatter(logm1Mem, logr1Mem, facecolor='none', edgecolor='b', alpha=0.2, marker='h', s=60)
#plt.scatter(logm2Mem, logr2Mem, facecolor='none', edgecolor='r', alpha=0.2, marker='h', s=60)
plt.scatter(logm1Bcg, logr1Bcg, c='b', alpha=0.35, s=55, marker='s')
plt.scatter(logm2Bcg, logr2Bcg, c='r', alpha=0.35, s=55, marker='s')
plt.show()
In [48]:
temp1 = 10.0 ** np.asarray(logm1Gama1)
temp2 = 10.0 ** np.asarray(logm2Gama1)
plt.scatter(np.log10(temp1 + temp2), (temp2 / (temp1 + temp2)), alpha=0.1, c='k')
temp1 = 10.0 ** np.asarray(logm1Gama2)
temp2 = 10.0 ** np.asarray(logm2Gama2)
plt.scatter(np.log10(temp1 + temp2), (temp2 / (temp1 + temp2)), alpha=0.2, c='k')
temp1 = 10.0 ** np.asarray(logm1Bcg)
temp2 = 10.0 ** np.asarray(logm2Bcg)
plt.scatter(np.log10(temp1 + temp2), (temp2 / (temp1 + temp2)), c='r', alpha=0.3, s=50)
plt.show()
In [133]:
import emcee
import corner
# Likelihood function
def lnlike(theta, x, y, yerr):
a, b, lnf = theta
# Straightline model
model = a * x + b
inv_sigma2 = (1.0 / (yerr**2 + model**2*np.exp(2*lnf)))
return -0.5*(np.sum((y-model)**2*inv_sigma2 -
np.log(inv_sigma2)))
# Priors
def lnprior(theta):
a, b, lnf = theta
if 0.6 < a < 1.3 and -13.0 < b < -7.0 and -4.0 < lnf < 1.0:
return 0.0
return -np.inf
def lnprob(theta, x, y, yerr):
lp = lnprior(theta)
if not np.isfinite(lp):
return -np.inf
return lp + lnlike(theta, x, y, yerr)
#Llinear least squares solution
def llsLine(x, y, yerr):
""" Simple straingt line fitting """
A = np.vstack((np.ones_like(x), x)).T
C = np.diag(yerr * yerr)
cov = np.linalg.inv(np.dot(A.T, np.linalg.solve(C, A)))
b_ls, a_ls = np.dot(cov, np.dot(A.T,
np.linalg.solve(C, y)))
print("LLS: a =%8.5f ; b =%8.5f" % (a_ls, b_ls))
return a_ls, b_ls
# Use Emcee to fit a straight line
def emceeLine(x, y, yerr, nwalkers=100, ndim=3, nburn=100,
nstep=600, show=True):
""" Initial guesses from simple LLS fitting """
a_ls, b_ls = llsLine(x, y, yerr)
initial = [a_ls, b_ls, 0.00]
""" Start the sampler """
sampler = emcee.EnsembleSampler(nwalkers, ndim,
lnprob, args=(x, y, yerr))
""" Initializing the walkers. """
np.random.seed(0)
guesses = [initial + (1e-2*np.random.randn(ndim))
for i in range(nwalkers)]
""" Run MCMC """
print("Start the MCMC runs")
%time sampler.run_mcmc(guesses, nstep)
print("Done")
""" Flatten the chain so that we have a flat list of samples """
samples = sampler.chain[:, nburn:, :].reshape(-1, ndim)
if show:
fig = corner.corner(samples,
labels=["$a$", "$b$", "$\ln\,f$"])
""" Compute the quantiles. """
samples[:, 2] = np.exp(samples[:, 2])
a_mcmc, b_mcmc, f_mcmc = map(lambda v: (v[1], v[2]-v[1], v[1]-v[0]),
zip(*np.percentile(samples, [16, 50, 84],
axis=0)))
print("""MCMC result:
a = {0[0]} +{0[1]} -{0[2]}
b = {1[0]} +{1[1]} -{1[2]}
f = {2[0]} +{2[1]} -{2[2]}
""".format(a_mcmc, b_mcmc, f_mcmc))
return a_mcmc, b_mcmc, f_mcmc
In [138]:
logrBcgErr = logrBcgU * 0.0 + 0.1
a_mcmc, b_mcmc, f_mcmc = emceeLine(logmBcgU, logrBcgU, logrBcgErr,
nwalkers=200, ndim=3, nburn=1000,
nstep=5000, show=True)
# Best paramters
aBcgB, aBcgU, aBcgL = a_mcmc
bBcgB, bBcgU, bBcgL = b_mcmc
plt.show()
In [145]:
logrBcgErr = logrBcgM * 0.0 + 0.1
a_mcmc, b_mcmc, f_mcmc = emceeLine(logmBcgM, logrBcgM, logrBcgErr,
nwalkers=200, ndim=3, nburn=1000,
nstep=5000, show=True)
# Best paramters
aBcgMB, aBcgMU, aBcgML = a_mcmc
bBcgMB, bBcgMU, bBcgML = b_mcmc
plt.show()
In [146]:
logrGamaErr = logrGamaU * 0.0 + 0.1
a_mcmc, b_mcmc, f_mcmc = emceeLine(logmGamaU, logrGamaU, logrGamaErr,
nwalkers=200, ndim=3, nburn=1000,
nstep=5000, show=True)
# Best paramters
aGamaB, aGamaU, aGamaL = a_mcmc
bGamaB, bGamaU, bGamaL = b_mcmc
plt.show()
In [147]:
logrGamaErr = logrGamaM * 0.0 + 0.1
a_mcmc, b_mcmc, f_mcmc = emceeLine(logmGamaM, logrGamaM, logrGamaErr,
nwalkers=200, ndim=3, nburn=1000,
nstep=5000, show=True)
# Best paramters
aGamaMB, aGamaMU, aGamaML = a_mcmc
bGamaMB, bGamaMU, bGamaML = b_mcmc
plt.show()
In [50]:
plt.xlim(11.4, 12.6)
plt.ylim(0.45, 2.59)
plt.xlabel('$\log\ (M_{\star}/M_{\odot})$', fontsize=35)
plt.ylabel('$\log\ (R_{\mathrm{e}}/\mathrm{Kpc})$', fontsize=35)
plt.scatter(logmGamaU, logrGamaU, c='k', alpha=0.05)
plt.scatter(logmBcgU, logrBcgU, c='r', alpha=0.5, s=70, marker='h')
plt.plot(mc, rc, linestyle='-', linewidth=9.0, c='r', alpha=0.7)
plt.plot(ms, rs, linestyle='--', linewidth=9.0, c='b', alpha=0.8)
# Red-sequence
xx = np.linspace(10.0, 13.0, 100)
plt.plot(xx, (aBcgMB * xx + bBcgMB),
linestyle='--', color='r', linewidth=3.5)
#yy1 = (aBcgMB + aBcgMU) * xx + (bBcgMB)
#yy2 = (aBcgMB - aBcgML) * xx + (bBcgMB)
#plt.fill_between(xx, yy1, yy2,
# facecolor='r', interpolate=True, alpha=0.15)
xx = np.linspace(10.0, 13.0, 100)
plt.plot(xx, (aGamaMB * xx + bGamaMB),
linestyle='-.', color='k', linewidth=3.5)
#yy1 = (aGamaMB + aGamaMU) * xx + (bGamaMB)
#yy2 = (aGamaMB - aGamaML) * xx + (bGamaMB)
#plt.fill_between(xx, yy1, yy2,
# facecolor='k', interpolate=True, alpha=0.15)
plt.show()
In [170]:
plt.xlim(11.4, 12.6)
plt.ylim(-0.19, 1.59)
plt.xlabel('$\log\ (M_{\star}/M_{\odot})$', fontsize=35)
plt.ylabel('$\log\ \gamma $', fontsize=40)
plt.scatter(logmGamaU, (logrGamaU + 0.94 * (11.0 - logmGamaU)),
c='k', alpha=0.05)
plt.scatter(logmBcgU, (logrBcgU + 0.94 * (11.0 - logmBcgU)),
c='r', alpha=0.5, s=70, marker='h')
plt.show()
In [173]:
plt.xlabel('$\log\ \gamma $', fontsize=40)
plt.hist((logrGamaM + 0.94 * (11.0 - logmGamaM)), 40, normed=True,
edgecolor='k', alpha=0.95, histtype='step', linewidth=4.0)
plt.hist((logrBcgM + 0.94 * (11.0 - logmBcgM)), 20, normed=True,
facecolor='r', alpha=0.4, histtype='stepfilled')
plt.show()
In [ ]: