In [1]:

    

%pylab inline
from __future__ import (division, print_function)

import os
import sys
import copy
import fnmatch
import warnings
import collections

import numpy as np
import scipy
try:
    from scipy.stats import scoreatpercentile
except:
    scoreatpercentile = False
from scipy.interpolate import interp1d
import cPickle as pickle

# Astropy
from astropy.io import fits
from astropy    import units as u
from astropy.stats import sigma_clip
from astropy.table import Table, Column
from astropy.utils.console import ProgressBar

# AstroML
from astroML.plotting import hist

# Matplotlib related
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.patches import Ellipse
from matplotlib.ticker import NullFormatter
from matplotlib.ticker import MaxNLocator
from matplotlib.collections import PatchCollection

# Matplotlib default settings
rcdef = plt.rcParams.copy()
pylab.rcParams['figure.figsize'] = 12, 10
pylab.rcParams['xtick.major.size'] = 8.0
pylab.rcParams['xtick.major.width'] = 2.5
pylab.rcParams['xtick.minor.size'] = 4.0
pylab.rcParams['xtick.minor.width'] = 2.5
pylab.rcParams['ytick.major.size'] = 8.0
pylab.rcParams['ytick.major.width'] = 2.5
pylab.rcParams['ytick.minor.size'] = 4.0
pylab.rcParams['ytick.minor.width'] = 2.5

# Personal
import hscUtils as hUtil
import galSBP
import coaddCutoutGalfitSimple as gSimple 

# Cosmology
import cosmology
c=cosmology.Cosmo(H0=70.0, omega_m=0.3, omega_l=0.7, flat=1)

# Color map 
from palettable.colorbrewer.sequential import Oranges_4, Blues_5
ORG4 = Oranges_4.mpl_colormap
BLU5 = Blues_5.mpl_colormap


    

    




Populating the interactive namespace from numpy and matplotlib






    




/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 [2]:

    

# Absolute magnitude of sun in HSC filters
# Actuall borrowed from DES filters
# Values from magsun.data in FSPS
amag_sun_des_g = 5.08
amag_sun_des_r = 4.62
amag_sun_des_i = 4.52
amag_sun_des_z = 4.52
amag_sun_des_y = 4.51

# Based on http://www.baryons.org/ezgal/filters.php
amag_sun_ukiss_y = 4.515

# Extinction correction factor for HSC 
## A\_lambda = Coeff * E(B-V) 
a_hsc_g = 3.233
a_hsc_r = 2.291 
a_hsc_i = 1.635
a_hsc_z = 1.261
a_hsc_y = 1.076

# 
SIGMA1 = 0.3173
SIGMA2 = 0.0455
SIGMA3 = 0.0027

RSMA_COMMON = np.arange(0.4, 4.2, 0.02)


    

In [9]:

    

# Code for Get Bootstrap mean or median 
def _confidence_interval_1d(A, alpha=.05, metric=np.mean, numResamples=10000, interpolate=True):
    """Calculates bootstrap confidence interval along one dimensional array"""
    
    if not isinstance(alpha, collections.Iterable):
        alpha = np.array([alpha])

    N = len(A)
    resampleInds = np.random.randint(0, N, (numResamples,N))
    metricOfResampled = metric(A[resampleInds], axis=-1)

    confidenceInterval = np.zeros(2*len(alpha),dtype='float')
    
    if interpolate:
        for thisAlphaInd, thisAlpha in enumerate(alpha):
            confidenceInterval[2*thisAlphaInd] = scoreatpercentile(metricOfResampled, 
                                                                   thisAlpha*100/2.0)
            confidenceInterval[2*thisAlphaInd+1] = scoreatpercentile(metricOfResampled, 
                                                                     100-thisAlpha*100/2.0)
    else:
        sortedMetricOfResampled = np.sort(metricOfResampled)
        for thisAlphaInd, thisAlpha in enumerate(alpha):
            confidenceInterval[2*thisAlphaInd] = sortedMetricOfResampled[int(round(thisAlpha*numResamples/2.0))]
            confidenceInterval[2*thisAlphaInd+1] = sortedMetricOfResampled[int(round(numResamples - 
                                                                                     thisAlpha*numResamples/2.0))]
    return confidenceInterval
    
def _ma_confidence_interval_1d(A, alpha=.05, metric=np.mean, numResamples=10000, interpolate=True):
    A = np.ma.masked_invalid(A, copy=True)
    A = A.compressed()
    confidenceInterval = _confidence_interval_1d(A, alpha, metric, numResamples, interpolate)
    return confidenceInterval

def confidence_interval(A, axis=None, alpha=.05, metric=np.mean, numResamples=10000, interpolate=True):
    """Return the bootstrap confidence interval of an array or along an axis ignoring NaNs and masked elements.
    
    Parameters
    ----------
    A : array_like
        Array containing numbers whose confidence interval is desired. 
    axis : int, optional
        Axis along which the confidence interval is computed.
        The default is to compute the confidence interval of the flattened array.
    alpha: float or array, optional
        confidence level of confidence interval. 100.0*(1-alpha) percent confidence 
        interval will be returned.
        If length-n array, n confidence intervals will be computed
        The default is .05
    metric : numpy function, optional
        metric to calculate confidence interval for.
        The default is numpy.mean
    numResamples : int, optional
        number of bootstrap samples. The default is 10000.
    interpolate: bool, optional
        uses scipy.stats.scoreatpercentile to interpolate between bootstrap samples 
        if alpha*numResamples/2.0 is not integer.
        The default is True
        
    Returns
    -------
    confidenceInterval : ndarray
    An array with the same shape as `A`, with the specified axis replaced by one twice the length of the alpha
    If `A` is a 0-d array, or if axis is None, a length-2 ndarray is returned.
    """
    if interpolate is True and scoreatpercentile is False:
        print("need scipy to interpolate between values")
        interpolate = False
    A = A.copy()
    if axis is None:
        A = A.ravel()
        outA = _ma_confidence_interval_1d(A, alpha, metric, numResamples, interpolate)
    else:
        outA = np.apply_along_axis(_ma_confidence_interval_1d, axis, A, alpha, 
                                   metric, numResamples, interpolate)
        
    return outA

def normProf(sma, sbp, minSma, maxSma, divide=False): 
    """
    Naive method to normalize the profile. 
    
    Parameters: 
        sbp    : Array for surface brightness profile 
        sma    : Radius range 
        minSma : Minimum SMA
        maxSma   Maximum SMA
    """
    offset = np.nanmedian(sbp[(sma >= minSma) & 
                              (sma <= maxSma)])
    if divide: 
        return (sbp / offset)
    else:
        return (sbp-offset)
    
    
def pixKpc(redshift, pix=0.168, show=True, npix=1.0):
    """
    Get the corresponding Kpc size of a pixel.  
    
    Parameters: 
    """
    pixKpc = pix * npix * hUtil.cosmoScale(redshift)

    if show:
        print("# %d pixel(s) = %6.3f Kpc" % (npix, pixKpc))
        
    return pixKpc


def logAdd(para1, para2):
    """ Useful for adding magnitudes. """
    return np.log10((10.0 ** np.asarray(para1)) + 
                    (10.0 ** np.asarray(para2)))


def errAdd(err1, err2):
    """Add error quadral..."""
    return np.sqrt((err1 ** 2.0) + 
                   (err2 ** 2.0))


def toColorArr(data, bottom=None, top=None):
    """ 
    Convert a data array to "color array" (between 0 and 1). 
    
    Parameters:
        bottom, top  : 
    """
    if top is not None:
        data[data >= top] = top
    if bottom is not None:
        data[data <= bottom] = bottom
        
    return ((data - np.nanmin(data)) / 
            (np.nanmax(data) - np.nanmin(data))) * 255.0


def getLuminosity(mag, redshift, extinction=None, 
                  amag_sun=None):
    """Get the absolute magnitude or luminosity."""
    distmod = hUtil.cosmoDistMod(redshift)
    absMag = (mag - distmod)
    if extinction is not None: 
        absMag -= extinction 
    if amag_sun is not None: 
        absMag = ((amag_sun - absMag) / 2.5)
    
    return absMag

def getStackProfiles(sample, loc, name='GAMA', 
                     idCol='ID_USE', tabCol='sum_tab', save=True):
    """Get the stacks of the profiles."""
    print("## Sample %s : Will deal with %d galaxies" % (name, len(sample)))
    profiles = []
    with ProgressBar(len(sample), ipython_widget=True) as bar:
        for g in sample:
            try:
                gFile = os.path.join(loc, g['sum_tab'].replace('./', '')).strip()
                gProf = Table.read(gFile, format='fits')
                """ Add extra information """
                try: 
                    gProf.meta['KCORRECT_I'] = g['KCORRECT_I']
                    gProf.meta['KCORRECT_b_I'] = g['KCORRECT_b_I']
                    gProf.meta['KCORRECT_c_I'] = g['KCORRECT_c_I']
                    gProf.meta['KCORRECT_G'] = g['KCORRECT_G']
                    gProf.meta['KCORRECT_b_G'] = g['KCORRECT_b_G']
                    gProf.meta['KCORRECT_c_G'] = g['KCORRECT_c_G']
                    gProf.meta['KCORRECT_R'] = g['KCORRECT_R']
                    gProf.meta['KCORRECT_b_R'] = g['KCORRECT_b_R']
                    gProf.meta['KCORRECT_c_R'] = g['KCORRECT_c_R']
                    gProf.meta['KCORRECT_Z'] = g['KCORRECT_Z']
                    gProf.meta['KCORRECT_b_Z'] = g['KCORRECT_b_Z']
                    gProf.meta['KCORRECT_c_Z'] = g['KCORRECT_c_Z']
                    gProf.meta['KCORRECT_Y'] = g['KCORRECT_Y']
                    gProf.meta['KCORRECT_b_Y'] = g['KCORRECT_b_Y']
                    gProf.meta['KCORRECT_c_Y'] = g['KCORRECT_c_Y']
                    gProf.meta['LOGM2LI_A'] = g['logm2lI_A']
                    gProf.meta['LOGM2LI_B'] = g['logm2lI_B']
                    gProf.meta['LOGM2LI_C'] = g['logm2lI_C']
                    gProf.meta['LUM_100'] = g['lum_100']
                    gProf.meta['LUM_120'] = g['lum_120']
                except Exception:
                    print("## WARNING: Some metadata may not be available !")
                    continue
            except Exception:
                print("## Missing: %s" % gFile)
                continue 
            profiles.append(gProf)
            bar.update()
    
    if save: 
        outPkl = os.path.join(loc, (name + '_profs.pkl'))
        hUtil.saveToPickle(profiles, outPkl)
        print("## Save %s to %s" % (name, outPkl))
        
    return profiles


def organizeSbp(profiles, col1='muI1', col2='KCORRECT_c_I', 
                kind='sbp', norm=False, r1=9.9, r2=10.1, divide=False,
                col3=None, col4=None, justStack=False,
                sun1=amag_sun_des_g, sun2=amag_sun_des_r,
                index=None):
    """ Get the stack of individual profiels, and their med/avg. """
    if kind.strip() == 'sbp':
        if col2 is not None: 
            if norm:
                stack = np.vstack(normProf(p['rKpc'], 
                                           np.asarray(p[col1] + (p.meta[col2] / 2.5)), 
                                           r1, r2, divide=divide) 
                                  for p in profiles)
            else:
                stack = np.vstack(np.asarray(p[col1] + (p.meta[col2] / 2.5)) 
                                  for p in profiles)
        else: 
            print("## NO KCORRECTION APPLIED !!")            
            if norm:
                stack = np.vstack(normProf(p['rKpc'], p[col1], 
                                           r1, r2, divide=divide) 
                                  for p in profiles)
            else:
                stack = np.vstack(np.asarray(p[col1]) for p in profiles)
    elif kind.strip() == 'mass':
        if norm:
            stack = np.vstack(normProf(p['rKpc'], 
                                       np.asarray(p[col1] + p.meta[col2]), 
                                       r1, r2, divide=divide) for p in profiles)
        else: 
            stack = np.vstack(np.asarray(p[col1] + p.meta[col2]) for p in profiles)
    elif kind.strip() == 'color':
        cSun = (sun1 - sun2)
        if col3 is None or col4 is None:
            print("## NO KCORRECTION APPLIED !!")
            if norm:
                stack = np.vstack(normProf(p['rKpc'], 
                                           np.asarray(cSun - 2.5 * (p[col1] - p[col2])), 
                                           r1, r2, divide=divide) for p in profiles)
            else: 
                stack = np.vstack(np.asarray(cSun - 2.5 *(p[col1] - p[col2])) for p in profiles)
        else:
            if norm:
                stack = np.vstack(normProf(p['rKpc'], 
                                           np.asarray(cSun - 2.5 * (p[col1] - p[col2]) -
                                                      (p.meta[col3] - p.meta[col4])), 
                                           r1, r2, divide=divide) for p in profiles)
            else: 
                stack = np.vstack(np.asarray(cSun - 2.5 * (p[col1] - p[col2]) -
                                             (p.meta[col3] - p.meta[col4])) 
                                  for p in profiles)
    elif kind.strip() == 'lum':
        if col2 is None:
            stack = np.vstack(np.asarray(p[col1]) for p in profiles)
        else:
            stack = np.vstack(np.asarray(p[col1] - p.meta[col2]) for p in profiles)
    else: 
        raise Exception("## WRONG KIND !!")
        
    if index is not None: 
        stack = np.vstack(p[index] for p in stack)
        
    if not justStack:
        """ Get the median and 1-sigma confidence range """
        medProf = confidence_interval(stack, axis=0, alpha=np.asarray([SIGMA1, 1.0]), 
                                      metric=np.nanmedian, numResamples=1000, 
                                      interpolate=True) 
        avgProf = confidence_interval(stack, axis=0, alpha=np.asarray([SIGMA1, 1.0]), 
                                      metric=np.nanmean, numResamples=1000, 
                                      interpolate=True) 
        stdProf = confidence_interval(stack, axis=0, alpha=np.asarray([SIGMA1, 1.0]), 
                                      metric=np.nanstd, numResamples=1000, 
                                      interpolate=True) 
        return stack, medProf, avgProf, stdProf
    else: 
        return stack
    

def loadPkl(filename):
    try:
        import cPickle as pickle
    except:
        warnings.warn("## cPickle is not available!!")
        import pickle
    
    if os.path.isfile(filename):
        pklFile = open(filename, 'rb')
        data = pickle.load(pklFile)    
        pklFile.close()
    
        return data
    else: 
        warnings.warn("## Can not find %s, return None" % filename)
        return None


    

In [84]:

    

newDir = '/Users/songhuang/work/hscs/gama_massive/sbp/'

bcgFile = 'redbcg_1d_160211.fits'
memFile = 'redmem_1d_160211.fits'
gamaFile = 'gama_1d_160211.fits'

try:
    bcgTab
except NameError:
    pass
else:
    del bcgTab
    
try:
    memTab
except NameError:
    pass
else:
    del memTab    
    
try:
    gamaTab
except NameError:
    pass
else:
    del gamaTab
    
#
bcgDir = os.path.join(newDir, 'redbcg')
memDir = os.path.join(newDir, 'redmem')
gamaDir = os.path.join(newDir, 'gama')

# Two summary catalogs
bcgCat = os.path.join(newDir, bcgFile)
memCat = os.path.join(newDir, memFile)
gamaCat = os.path.join(newDir, gamaFile)

if not os.path.isfile(bcgCat):
    raise Exception("## Can not find catalog for BCGs : %s" % bcgCat)
else: 
    bcgTab = Table.read(bcgCat, format='fits')

if not os.path.isfile(memCat):
    raise Exception("## Can not find catalog for cluster members : %s" % memCat)
else: 
    memTab = Table.read(memCat, format='fits')
    
if not os.path.isfile(gamaCat):
    raise Exception("## Can not find catalog for GAMA galaxies : %s" % gamaCat)
else: 
    gamaTab = Table.read(gamaCat, format='fits')
    
print("## Deal with %i galaxies in redBCH sample" % len(bcgTab))
print("## Deal with %i galaxies in redMEM sample" % len(memTab))
print("## Deal with %i galaxies in GAMA sample" % len(gamaTab))


    

    




## Deal with 219 galaxies in redBCH sample
## Deal with 1542 galaxies in redMEM sample
## Deal with 9414 galaxies in GAMA sample

In [93]:

    

bcgClean = bcgTab[(bcgTab['m100_c'] >= 10.0) & 
                  (bcgTab['c82_120'] <= 14.0) &
                  (bcgTab['r90_max'] <= 220.0) &
                  (bcgTab['P_CEN_1'] >= 0.8)]
print(len(bcgClean))

gamaClean = gamaTab[(gamaTab['c82_120'] >= 5.5) & 
                    (gamaTab['c82_120'] <= 14.0) &
                    (gamaTab['gz_kC'] >= 1.58) &
                    (gamaTab['m100_c'] >= 10.0) &
                    (gamaTab['r90_max'] <= 160.0) &
                    (gamaTab['ur_rest_sed'] >= 2.1) &
                    ((gamaTab['r50_120'] - gamaTab['r20_120']) >= 
                     (32.0 * (gamaTab['m30_c'] - gamaTab['m10_c']) - 1.4))]

print(len(gamaClean))

memClean = memTab[(memTab['c82_120'] >= 5.0) &
                  (memTab['c82_120'] <= 14.0) &
                  (memTab['gz_kC'] >= 1.45) &
                  (memTab['m100_c'] >= 10.0) &
                  (memTab['r90_max'] <= 160.0) &
                  ((memTab['r50_120'] - memTab['r20_120']) >= 
                   (32.0 * (memTab['m30_c'] - memTab['m10_c']) - 1.4))]

print(len(memClean))


    

    




130
3217
777






    




/usr/local/lib/python2.7/site-packages/ipykernel/__main__.py:12: RuntimeWarning: invalid value encountered in greater_equal

In [94]:

    

gamaUse = gamaClean[(gamaClean['z_use'] >= 0.20) &
                    (gamaClean['z_use'] <= 0.45)]

gamaUse2 = gamaClean[(gamaClean['z_use'] >= 0.20) &
                     (gamaClean['z_use'] <= 0.40)]
print(len(gamaUse))

print(len(gamaUse2))

memUse = memClean[(memClean['z_use'] >= 0.20) &
                  (memClean['z_use'] <= 0.40)]
print(len(memUse))

bcgUse = bcgClean[(bcgClean['z_use'] >= 0.20) &
                  (bcgClean['z_use'] <= 0.45) &
                  (bcgClean['LAMBDA_CLUSTER'] >= 25.0)]

bcgUse2 = bcgClean[(bcgClean['z_use'] >= 0.20) &
                   (bcgClean['z_use'] <= 0.40) &
                   (bcgClean['LAMBDA_CLUSTER'] >= 25.0)]
print(len(bcgUse))
print(len(bcgUse2))


    

    




2750
2643
777
73
51

In [98]:

    

print("########################################################\n")
#gamaM1 = gamaClean[(gamaClean['m100_c'] >= 11.56) & 
#                   (gamaClean['m100_c'] < 11.70) &
#                   (gamaClean['Z'] >= 0.20) &
#                   (gamaClean['Z'] <= 0.40)]
#gamaM2 = gamaClean[(gamaClean['m100_c'] >= 11.72) & 
#                   (gamaClean['m100_c'] < 11.95) &
#                   (gamaClean['Z'] >= 0.20) &
#                   (gamaClean['Z'] <= 0.50)]
gamaM1 = gamaClean[(gamaClean['m100_c'] >= 11.55) & 
                   (gamaClean['m100_c'] < 11.70) &
                   (gamaClean['z_use'] >= 0.20) &
                   (gamaClean['z_use'] <= 0.45)]

gamaM2 = gamaClean[(gamaClean['m100_c'] >= 11.71) & 
                   (gamaClean['m100_c'] < 11.90) &
                   (gamaClean['z_use'] >= 0.20) &
                   (gamaClean['z_use'] <= 0.50)]

gamaM1b = gamaClean[(gamaClean['m100_c'] >= 11.55) & 
                    (gamaClean['m100_c'] < 11.70) &
                    (gamaClean['z_use'] >= 0.20) &
                    (gamaClean['z_use'] <= 0.40)]

gamaM2b = gamaClean[(gamaClean['m100_c'] >= 11.70) & 
                    (gamaClean['m100_c'] < 11.90) &
                    (gamaClean['z_use'] >= 0.20) &
                    (gamaClean['z_use'] <= 0.40)]

print(len(gamaM1), np.nanmedian(gamaM1['m100_c']), np.nanmedian(gamaM1['z_use']))
print(len(gamaM2), np.nanmedian(gamaM2['m100_c']), np.nanmedian(gamaM2['z_use']))

print(len(gamaM1b), np.nanmedian(gamaM1b['m100_c']), np.nanmedian(gamaM1b['z_use']))
print(len(gamaM2b), np.nanmedian(gamaM2b['m100_c']), np.nanmedian(gamaM2b['z_use']))

print("########################################################\n")
memM1 = memClean[(memClean['m100_c'] >= 11.55) & 
                 (memClean['m100_c'] < 11.70) &
                 (memClean['z_use'] >= 0.20) &
                 (memClean['z_use'] <= 0.40)]

memM2 = memClean[(memClean['m100_c'] >= 11.70) & 
                 (memClean['m100_c'] < 11.90) &
                 (memClean['z_use'] >= 0.20) &
                 (memClean['z_use'] <= 0.40)]

print(len(memM1), np.nanmedian(memM1['m100_c']), np.nanmedian(memM1['Z']))
print(len(memM2), np.nanmedian(memM2['m100_c']), np.nanmedian(memM2['Z']))

print("########################################################\n")
#bcgM1 = bcgClean[(bcgClean['m100_c'] >= 11.45) & 
#                 (bcgClean['m100_c'] < 11.68) &
#                 (bcgClean['Z'] >= 0.20) &
#                 (bcgClean['Z'] <= 0.40) &
#                 (bcgClean['P_CEN_1'] >= 0.8)]
#bcgM2 = bcgClean[(bcgClean['m100_c'] >= 11.70) & 
#                 (bcgClean['m100_c'] < 11.95) &
#                 (bcgClean['Z'] >= 0.20) &
#                 (bcgClean['Z'] <= 0.50) &
#                 (bcgClean['P_CEN_1'] >= 0.8) &
#                 (bcgClean['LAMBDA_CLUSTER'] >= 30.0)]

bcgM1 = bcgClean[(bcgClean['m100_c'] >= 11.55) & 
                 (bcgClean['m100_c'] < 11.70) &
                 (bcgClean['z_use'] >= 0.20) &
                 (bcgClean['z_use'] <= 0.40) &
                 (bcgClean['LAMBDA_CLUSTER'] >= 25.0)]

bcgM2 = bcgClean[(bcgClean['m100_c'] >= 11.70) & 
                 (bcgClean['m100_c'] < 11.90) &
                 (bcgClean['z_use'] >= 0.20) &
                 (bcgClean['z_use'] <= 0.50) &
                 (bcgClean['LAMBDA_CLUSTER'] >= 25.0)]

bcgM1b = bcgClean[(bcgClean['m100_c'] >= 11.50) & 
                  (bcgClean['m100_c'] < 11.70) &
                  (bcgClean['z_use'] >= 0.20) &
                  (bcgClean['z_use'] <= 0.40) &
                  (bcgClean['LAMBDA_CLUSTER'] >= 25.0)]

bcgM2b = bcgClean[(bcgClean['m100_c'] >= 11.70) & 
                  (bcgClean['m100_c'] < 11.90) &
                  (bcgClean['z_use'] >= 0.20) &
                  (bcgClean['z_use'] <= 0.40) &
                  (bcgClean['LAMBDA_CLUSTER'] >= 25.0)]

print(len(bcgM1), np.nanmedian(bcgM1['m100_c']), np.nanmedian(bcgM1['z_use']))
print(len(bcgM2), np.nanmedian(bcgM2['m100_c']), np.nanmedian(bcgM2['z_use']))
print(len(bcgM1b), np.nanmedian(bcgM1b['m100_c']), np.nanmedian(bcgM1b['z_use']))
print(len(bcgM2b), np.nanmedian(bcgM2b['m100_c']), np.nanmedian(bcgM2b['z_use']))


    

    




########################################################

447 11.5729275131 0.291509985924
122 11.7669058203 0.322205007076
285 11.6107901021 0.295269995928
117 11.7540580411 0.306050002575
########################################################

36 11.6120771254 0.306126177311
12 11.7758495468 0.312638521194
########################################################

20 11.6198632792 0.307940006256
22 11.7688292191 0.336626768112
20 11.6198632792 0.307940006256
15 11.7558052257 0.297216922045

In [51]:

    

print("########################################################\n")
gamaCM1 = gamaClean[(gamaClean['m10_c'] >= 11.19) & 
                    (gamaClean['m10_c'] < 11.35) &
                    (gamaClean['Z'] >= 0.20) &
                    (gamaClean['Z'] <= 0.40)]

gamaCM1b = gamaCM1[(gamaCM1['c82_120'] >= 7.0) & 
                   (gamaCM1['r90_120'] / gamaCM1['r50_120'] >= 2.3) &
                   (gamaCM1['gz_kC'] >= 1.6) &
                   (gamaCM1['ur_rest_sed'] >= 2.2)]

gamaCM2 = gamaClean[(gamaClean['m10_c'] >= 11.35) & 
                    (gamaClean['m10_c'] < 11.55) &
                    (gamaClean['Z'] >= 0.20) &
                    (gamaClean['Z'] <= 0.40)]

print(len(gamaCM1), np.nanmedian(gamaCM1['m10_c']), 
      np.nanmedian(gamaCM1['Z']))
print(len(gamaCM2), np.nanmedian(gamaCM2['m10_c']), 
      np.nanmedian(gamaCM2['Z']))
print(len(gamaCM1b), np.nanmedian(gamaCM1b['m10_c']), 
      np.nanmedian(gamaCM1b['Z']))

print("########################################################\n")
memCM1 = memClean[(memClean['m10_c'] >= 11.18) & 
                  (memClean['m10_c'] < 11.35) &
                  (memClean['Z'] >= 0.20) &
                  (memClean['Z'] <= 0.40)]

memCM2 = memClean[(memClean['m10_c'] >= 11.35) & 
                  (memClean['m10_c'] < 11.55) &
                  (memClean['Z'] >= 0.20) &
                  (memClean['Z'] <= 0.50)]

print(len(memCM1), np.nanmedian(memCM1['m10_c']), 
      np.nanmedian(memCM1['Z']))
print(len(memCM2), np.nanmedian(memCM2['m10_c']), 
      np.nanmedian(memCM2['Z']))

print("########################################################\n")

bcgCM1 = bcgClean[(bcgClean['m10_c'] >= 11.15) & 
                  (bcgClean['m10_c'] < 11.34) &
                  (bcgClean['Z'] >= 0.20) &
                  (bcgClean['Z'] <= 0.40) &
                  (bcgClean['LAMBDA_CLUSTER'] >= 30.0)]

bcgCM1b = bcgClean[(bcgClean['m10_c'] >= 11.15) & 
                   (bcgClean['m10_c'] < 11.34) &
                   (bcgClean['m100_c'] <= 11.8) &
                   (bcgClean['Z'] >= 0.20) &
                   (bcgClean['Z'] <= 0.39)]

bcgCM2 = bcgClean[(bcgClean['m10_c'] >= 11.35) & 
                  (bcgClean['m10_c'] < 11.55) &
                  (bcgClean['Z'] >= 0.20) &
                  (bcgClean['Z'] <= 0.50) &
                  (bcgClean['LAMBDA_CLUSTER'] >= 30.0)]

print(len(bcgCM1), np.nanmean(bcgCM1['m10_c']), 
      np.nanmedian(bcgCM1['Z']))
print(len(bcgCM2), np.nanmean(bcgCM2['m10_c']), 
      np.nanmedian(bcgCM2['Z']))
print(len(bcgCM1b), np.nanmean(bcgCM1b['m10_c']), 
      np.nanmedian(bcgCM1b['Z']))


    

    




########################################################

732 11.2573833387 0.291665017605
211 11.3947089786 0.306989997625
71 11.2648783536 0.302630007267
########################################################

131 11.2431288653 0.314759999514
22 11.3921080017 0.304781168699
########################################################

15 11.2523061893 0.313470005989
13 11.4152549156 0.340429991484
27 11.2469846747 0.302599996328






    




/usr/local/lib/python2.7/site-packages/ipykernel/__main__.py:10: RuntimeWarning: invalid value encountered in greater_equal

In [66]:

    

gM1_mass = confidence_interval(gamaM1['m100_c'], axis=0, alpha=[SIGMA1, 1.0], 
                               metric=np.median, numResamples=1000, interpolate=True)
gCM1_mass = confidence_interval(gamaM1['m10_c'], axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
gM1_mdif1 = confidence_interval((gamaM1['lum_100'] - gamaM1['lum_10']), 
                                axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
gM1_mdif2 = confidence_interval((gamaM1['lum_75'] - gamaM1['lum_10']), 
                                axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
gM1_mdif3 = confidence_interval((gamaM1['lum_50'] - gamaM1['lum_10']), 
                                axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
gM1_mdif4 = confidence_interval((gamaM1['lum_25'] - gamaM1['lum_10']), 
                                axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
gM1_r20 = confidence_interval(gamaM1['r20_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
gM1_r50 = confidence_interval(gamaM1['r50_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
gM1_r90 = confidence_interval(gamaM1['r90_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
gM1_c82 = confidence_interval(gamaM1['c82_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
gM1_c52 = confidence_interval((gamaM1['r50_120'] / gamaM1['r20_120']), 
                              axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)

gM2_mass = confidence_interval(gamaM2['m100_c'], axis=0, alpha=[SIGMA1, 1.0], 
                               metric=np.median, numResamples=1000, interpolate=True)
gCM2_mass = confidence_interval(gamaM2['m10_c'], axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
gM2_mdif1 = confidence_interval((gamaM2['lum_100'] - gamaM2['lum_10']), 
                                axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
gM2_mdif2 = confidence_interval((gamaM2['lum_75'] - gamaM2['lum_10']), 
                                axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
gM2_mdif3 = confidence_interval((gamaM2['lum_50'] - gamaM2['lum_10']), 
                                axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
gM2_mdif4 = confidence_interval((gamaM2['lum_25'] - gamaM2['lum_10']), 
                                axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
gM2_r20 = confidence_interval(gamaM2['r20_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
gM2_r50 = confidence_interval(gamaM2['r50_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
gM2_r90 = confidence_interval(gamaM2['r90_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
gM2_c82 = confidence_interval(gamaM2['c82_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
gM2_c95 = confidence_interval((gamaM2['r50_120'] / gamaM2['r20_120']), 
                              axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)


    

In [67]:

    

bM1_mass = confidence_interval(bcgM1['m100_c'], axis=0, alpha=[SIGMA1, 1.0], 
                               metric=np.median, numResamples=1000, interpolate=True)
bCM1_mass = confidence_interval(bcgM1['m10_c'], axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
bM1_mdif1 = confidence_interval((bcgM1['lum_100'] - bcgM1['lum_10']), 
                               axis=0, alpha=[SIGMA1, 1.0], 
                               metric=np.median, numResamples=1000, interpolate=True)
bM1_mdif2 = confidence_interval((bcgM1['lum_75'] - bcgM1['lum_10']), 
                               axis=0, alpha=[SIGMA1, 1.0], 
                               metric=np.median, numResamples=1000, interpolate=True)
bM1_mdif3 = confidence_interval((bcgM1['lum_50'] - bcgM1['lum_10']), 
                               axis=0, alpha=[SIGMA1, 1.0], 
                               metric=np.median, numResamples=1000, interpolate=True)
bM1_mdif4 = confidence_interval((bcgM1['lum_25'] - bcgM1['lum_10']), 
                               axis=0, alpha=[SIGMA1, 1.0], 
                               metric=np.median, numResamples=1000, interpolate=True)
bM1_r20 = confidence_interval(bcgM1['r20_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
bM1_r50 = confidence_interval(bcgM1['r50_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
bM1_r90 = confidence_interval(bcgM1['r90_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
bM1_c82 = confidence_interval(bcgM1['c82_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
bM1_c52 = confidence_interval((bcgM1['r50_120'] / bcgM1['r20_120']), 
                              axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)

bM2_mass = confidence_interval(bcgM2['m100_c'], axis=0, alpha=[SIGMA1, 1.0], 
                               metric=np.median, numResamples=1000, interpolate=True)
bCM2_mass = confidence_interval(bcgM2['m10_c'], axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
bM2_mdif1 = confidence_interval((bcgM2['lum_100'] - bcgM2['lum_10']), 
                                axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
bM2_mdif2 = confidence_interval((bcgM2['lum_75'] - bcgM2['lum_10']), 
                                axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
bM2_mdif3 = confidence_interval((bcgM2['lum_50'] - bcgM2['lum_10']), 
                                axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
bM2_mdif4 = confidence_interval((bcgM2['lum_25'] - bcgM2['lum_10']), 
                                axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
bM2_r20 = confidence_interval(bcgM2['r20_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
bM2_r50 = confidence_interval(bcgM2['r50_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
bM2_r90 = confidence_interval(bcgM2['r90_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
bM2_c82 = confidence_interval(bcgM2['c82_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
bM2_c95 = confidence_interval((bcgM2['r50_120'] / bcgM2['r20_120']), 
                              axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)


    

In [68]:

    

mM1_mass = confidence_interval(memM1['m100_c'], axis=0, alpha=[SIGMA1, 1.0], 
                               metric=np.median, numResamples=1000, interpolate=True)
mCM1_mass = confidence_interval(memM1['m10_c'], axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
mM1_mdif1 = confidence_interval((memM1['lum_100'] - memM1['lum_10']), 
                                axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
mM1_mdif2 = confidence_interval((memM1['lum_75'] - memM1['lum_10']), 
                                axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
mM1_mdif3 = confidence_interval((memM1['lum_50'] - memM1['lum_10']), 
                                axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
mM1_mdif4 = confidence_interval((memM1['lum_25'] - memM1['lum_10']), 
                                axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
mM1_r20 = confidence_interval(memM1['r20_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
mM1_r50 = confidence_interval(memM1['r50_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
mM1_r90 = confidence_interval(memM1['r90_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
mM1_c82 = confidence_interval(memM1['c82_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
mM1_c52 = confidence_interval((memM1['r50_120'] / memM1['r20_120']), 
                              axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)

mM2_mass = confidence_interval(memM2['m100_c'], axis=0, alpha=[SIGMA1, 1.0], 
                               metric=np.median, numResamples=1000, interpolate=True)
mCM2_mass = confidence_interval(memM2['m10_c'], axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
mM2_mdif1 = confidence_interval((memM2['lum_100'] - memM2['lum_10']), 
                                axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
mM2_mdif2 = confidence_interval((memM2['lum_75'] - memM2['lum_10']), 
                                axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
mM2_mdif3 = confidence_interval((memM2['lum_50'] - memM2['lum_10']), 
                                axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
mM2_mdif4 = confidence_interval((memM2['lum_25'] - memM2['lum_10']), 
                                axis=0, alpha=[SIGMA1, 1.0], 
                                metric=np.median, numResamples=1000, interpolate=True)
mM2_r20 = confidence_interval(memM2['r20_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
mM2_r50 = confidence_interval(memM2['r50_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
mM2_r90 = confidence_interval(memM2['r90_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
mM2_c82 = confidence_interval(memM2['c82_120'], axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)
mM2_c95 = confidence_interval((memM2['r50_120'] / memM2['r20_120']), 
                              axis=0, alpha=[SIGMA1, 1.0], 
                              metric=np.median, numResamples=1000, interpolate=True)


    

In [69]:

    

plt.hist(gamaM1['m100_c'], bins=20, range=[11.4, 11.8], color='k', 
         histtype='step', normed=1)
plt.hist(bcgM1['m100_c'], bins=20, range=[11.4, 11.8], color='r', 
         histtype='step', normed=1)
plt.show()


    

In [56]:

    

plt.hist(gamaM2['m100_c'], bins=20, range=[11.7, 12.2], color='k', 
         histtype='step', normed=1)
plt.hist(bcgM2['m100_c'], bins=20, range=[11.7, 12.2], color='r', 
         histtype='step', normed=1)
plt.hist(gamaM2b['m100_c'], bins=20, range=[11.7, 12.2], color='k', 
         histtype='step', normed=1, linestyle='--')
plt.hist(bcgM2b['m100_c'], bins=20, range=[11.7, 12.2], color='r', 
         histtype='step', normed=1, linestyle='--')
plt.show()


    

In [426]:

    

gamaM1.write(os.path.join(newDir, 'test_gama_m1.fits'), format='fits', overwrite=True)
gamaM2.write(os.path.join(newDir, 'test_gama_m2.fits'), format='fits', overwrite=True)
gamaM2b.write(os.path.join(newDir, 'test_gama_m2b.fits'), format='fits', overwrite=True)

bcgM1.write(os.path.join(newDir, 'test_bcg_m1.fits'), format='fits', overwrite=True)
bcgM2.write(os.path.join(newDir, 'test_bcg_m2.fits'), format='fits', overwrite=True)
bcgM2b.write(os.path.join(newDir, 'test_bcg_m2b.fits'), format='fits', overwrite=True)

memM1.write(os.path.join(newDir, 'test_mem_m1.fits'), format='fits', overwrite=True)
memM2.write(os.path.join(newDir, 'test_mem_m2.fits'), format='fits', overwrite=True)

gamaCM1.write(os.path.join(newDir, 'test_gama_cm1.fits'), format='fits', overwrite=True)
gamaCM2.write(os.path.join(newDir, 'test_gama_cm2.fits'), format='fits', overwrite=True)

bcgCM1.write(os.path.join(newDir, 'test_bcg_cm1.fits'), format='fits', overwrite=True)
bcgCM2.write(os.path.join(newDir, 'test_bcg_cm2.fits'), format='fits', overwrite=True)


    

In [18]:

    

### Fancy one 

# definitions for the axes
left, width    = 0.12, 0.69
right          = left + width 
bottom, height = 0.12, 0.86
bottom_h = left_h = left + width + 0.02

recScat = [left,   bottom, width, height]
recHist = [right,  bottom,  0.18, height]


SBP1 = [0.13, 0.12, 0.865, 0.30]
SBP2 = [0.13, 0.42, 0.865, 0.54]


    

In [19]:

    

# Color 

BLUE0 = "#92c5de"
BLUE1 = "#0571b0"

RED0 = "#f4a582"
RED1 = "#ca0020"

PURPLE0 = '#af8dc3'
PURPLE1 = '#762a83'

BROWN0 = '#bf812d'
BROWN1 = '#543005'

GREEN0 = '#7fbf7b'
GREEN1 = '#1b7837'


    

In [72]:

    

fig = plt.figure(figsize=(13, 9))
ax1 = plt.axes(recScat)
ax2 = plt.axes(recHist)

# SBP v.s. (cModel - SBP)
# ---------------------------------------------------------------------------
# Scatter plot
#ax1.axhline(0.0, linewidth=4.0, linestyle='-', c='k', alpha=0.2)

ax1.axvline(11.5, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)
ax1.axvline(11.7, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)
ax1.axvline(11.9, linewidth=4.0, linestyle='--', c='k', alpha=0.2, zorder=0)

# Matched ones 
p1 = ax1.scatter(gamaClean['m100_c'], 
                 gamaClean['lum_100'] - gamaClean['lum_10'], s=35.0, 
                 alpha=0.20, facecolor=BLUE0, edgecolor='none', 
                 label='$\Lambda \leq 20\ \mathrm{Central}$')
p2 = ax1.scatter(bcgUse['m100_c'], 
                 bcgUse['lum_100'] - bcgUse['lum_10'], edgecolor='none',
                 s=((bcgUse['z_use'] - 0.10) * 600.0), cmap=ORG4, alpha=0.90, 
                 c=toColorArr(bcgUse['LAMBDA_CLUSTER'], bottom=20.0, top=70.0), 
                 label='$\Lambda > 20\ \mathrm{Central}$', marker='s')
"""
p3 = ax1.scatter(bcgM1['m100_c'], 
                 bcgM1['lum_100'] - bcgM1['lum_10'], edgecolor='k',
                 s=((bcgM1['z_use'] - 0.10) * 500.0), alpha=0.95, 
                 facecolor='none', label=None, marker='s', linewidth=1.5)
p4 = ax1.scatter(bcgM2b['m100_c'], 
                 bcgM2b['lum_100'] - bcgM2b['lum_10'], edgecolor='k',
                 s=((bcgM2['z_use'] - 0.10) * 500.0), alpha=0.95, 
                 facecolor='none', label=None, marker='s', linewidth=1.5)
"""

# M1
ax1.errorbar(gM1_mass[2], gM1_mdif1[2], marker='+', ms=1, mec='k',
             yerr=0.01, color=BLUE1, ecolor=BLUE1, capthick=3.5, capsize=8, 
             alpha=0.8, linewidth=4.0, fmt='h', elinewidth=2.0, label=None, 
             zorder=100)
ax1.errorbar(bM1_mass[2], bM1_mdif1[2], marker='+', ms=1, mec='k', linewidth=4.0,
             yerr=0.02, mfc=RED1, ecolor=RED1, capthick=3.5, capsize=8, 
             alpha=0.8, fmt='h', elinewidth=2.0, label=None, zorder=100)

ax1.scatter(gM1_mass[2], gM1_mdif1[2], marker='^', s=400, facecolor=BLUE1,
            edgecolor='k', linewidth=3.0, zorder=102, alpha=0.9,
            label='$\mathrm{[11.5,11.7]}\ \Lambda \leq 20$')
ax1.scatter(bM1_mass[2], bM1_mdif1[2], marker='p', s=420, facecolor=RED1,
            edgecolor='k', linewidth=3.0, zorder=102,
            label='$\mathrm{[11.5,11.7]}\ \Lambda > 30$')

# M2
ax1.errorbar(gM2_mass[2], gM2_mdif1[2], marker='+', ms=1, mec='k',
             yerr=0.01, color=BLUE1, ecolor=BLUE1, capthick=3.5, capsize=8, 
             alpha=0.8, fmt='h', elinewidth=2.0, label=None, zorder=100)
ax1.errorbar(bM2_mass[2], bM2_mdif1[2], marker='+', ms=1, mec='k',
             yerr=0.03, mfc=RED1, ecolor=RED1, capthick=3.5, capsize=8, 
             alpha=0.8, fmt='h', elinewidth=2.0, label=None, zorder=100)

ax1.scatter(gM2_mass[2], gM2_mdif1[2], marker='h', s=420, facecolor=BLUE1,
            edgecolor='k', linewidth=3.0, zorder=102, alpha=0.9,
            label='$\mathrm{[11.7,11.9]}\ \Lambda \leq 20$')
ax1.scatter(bM2_mass[2], bM2_mdif1[2], marker='8', s=420, facecolor=RED1,
            edgecolor='k', linewidth=3.0, zorder=102,
            label='$\mathrm{[11.7,11.9]}\ \Lambda > 30$')

# Legend
ax1.legend(loc=(0.68, 0.025), shadow=True, fancybox=True, 
           numpoints=1, fontsize=18, scatterpoints=1, 
           markerscale=0.9, borderpad=0.25, handletextpad=0.1)

legend = ax1.get_legend()
legend.legendHandles[1].set_color(ORG4(0.8))
legend.legendHandles[0].set_sizes([150])
legend.legendHandles[1].set_sizes([200])

#ax1.text(0.05, 0.04, '$\mathrm{Size:}\ {\Lambda}_{\mathrm{redMapper}}$', 
#         verticalalignment='bottom', horizontalalignment='left',
#         fontsize=26.0, transform=ax1.transAxes, color=RED0)

# Axes setup
#  Minor Ticks on 
ax1.minorticks_on()

#  Axes Thickness
for axis in ['top','bottom','left','right']:
  ax1.spines[axis].set_linewidth(3.5)
#  Tick Label Size 
for tick in ax1.xaxis.get_major_ticks():
    tick.label.set_fontsize(24) 
for tick in ax1.yaxis.get_major_ticks():
    tick.label.set_fontsize(24) 
#  Tick Length and Width
ax1.tick_params('both', length=10, width=3.0, which='major')
ax1.tick_params('both', length=6, width=2.5, which='minor')

# Label
ax1.set_xlabel('$\log\ (M_{\star}/M_{\odot})\ (100\ \mathrm{Kpc})$', size=40)
ax1.set_ylabel('$\Delta(\log M{\star})_{\mathrm{100\ kpc}-\mathrm{10\ kpc}}$', 
               size=42)

# Axis limits
ax1.set_xlim(11.15, 12.29)
ax1.set_ylim(0.01, 0.79)

# ---------------------------------------------------------------------------
# Histogram 
#
n, bins, patches=ax2.hist(gamaM1['lum_100'] - gamaM1['lum_10'], 
                          bins=30, range=[0.05, 0.8], edgecolor='none',
                          orientation='horizontal', histtype='stepfilled', 
                          color=BLUE0, alpha=0.80, normed=1)

n, bins, patches=ax2.hist(bcgM1['lum_100'] - bcgM1['lum_10'], 
                          bins=20, range=[0.05, 0.8], edgecolor='none',
                          orientation='horizontal', histtype='stepfilled', 
                          color=ORG4(0.6), alpha=0.50, normed=1, linewidth=4.0)

#
"""
n, bins, patches=ax2.hist(gamaM2['lum_100'] - gamaM2['lum_10'], 
                          bins=30, range=[0.0, 0.7], linewidth=4.0,
                          orientation='horizontal', histtype='step', 
                          color=BLUE0, alpha=1.0, normed=1)

n, bins, patches=ax2.hist(bcgM2['lum_100'] - bcgM2['lum_10'], 
                          bins=15, range=[0.05, 0.7], linewidth=4.0,
                          orientation='horizontal', histtype='step', 
                          color=RED0, alpha=0.90, normed=1)
"""


ax2.set_ylim(ax1.get_ylim())

ax2.axhline(0.0, linewidth=4.0, linestyle='-', c='k', alpha=0.2)

# Axes setup
# Minor Ticks on 
ax2.minorticks_on()
ax2.tick_params(axis='x', which='minor', bottom='off', top='off')

#  Axes Thickness
for axis in ['top','bottom','left','right']:
  ax2.spines[axis].set_linewidth(3.5)
#  Tick Label Size 
for tick in ax2.xaxis.get_major_ticks():
    tick.label.set_fontsize(24) 
for tick in ax2.yaxis.get_major_ticks():
    tick.label.set_fontsize(24) 
#  Tick Length and Width
ax2.tick_params('both', length=10, width=3.0, which='major')
ax2.tick_params('both', length=6, width=2.5, which='minor')

ax1.axhline(0.0, linewidth=4.5, linestyle='-', c='k', alpha=0.2)

ax2.yaxis.set_major_formatter(NullFormatter())
ax2.xaxis.set_major_formatter(NullFormatter())

plt.show()

fig.savefig('../figure/hscMassive_mtot_m100_10.png', dpi=300)