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


/usr/local/lib/python2.7/site-packages/matplotlib/__init__.py:1350: UserWarning:  This call to matplotlib.use() has no effect
because the backend has already been chosen;
matplotlib.use() must be called *before* pylab, matplotlib.pyplot,
or matplotlib.backends is imported for the first time.

  warnings.warn(_use_error_msg)

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

Results by 2016-01-11

  • Three major datasets:

    1. redBCG - /Users/songhuang/work/hscs/gama_massive/galfit/redbcg: xx galaxies (vetted)
      • redbcg_mass_use_dom.fits
    2. redMem - /Users/songhuang/work/hscs/gama_massive/galfit/redmem: xxx galaxies (not vetted)
      • redmapper_mem_hscmatch_mass_sbpsum_modA_muI1.fits
    3. gama - /Users/songhuang/work/hscs/gama_massive/galfit/gama: xxx galaxies (not vetted)
      • gama_massive_160107_sbpsum_modA_muI1.fits
      • Will be separated into different samples: GAMA1, 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')


WARNING: UnitsWarning: 'dex' did not parse as fits unit: At col 0, Unit 'dex' not supported by the FITS standard.  [astropy.units.core]
WARNING:astropy:UnitsWarning: 'dex' did not parse as fits unit: At col 0, Unit 'dex' not supported by the FITS standard. 

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))


## Deal with 219 galaxies in BCG sample
## Deal with 1670 galaxies in Cluster member sample
## Deal with 3747 galaxies in GAMA1 sample
## Deal with 1700 galaxies in GAMA2 sample
## Deal with 2718 galaxies in GAMA3 sample

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])

BCGs

Estimate the luminosity, mass-to-light ratio, and stellar mass

- Consider the the K-correction 
- Consider the differen mass models

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))

Go through the catalog, and get the lists of .pkl files

BCGs


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))


183
126

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()


Members: (The most massive bin has problem !)


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))


1420
586

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()


GAMA1


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))


3400
1564

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()


GAMA2


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))


1470
779

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()