In [1]:
%pylab inline

from __future__ import division

import os
import copy
import argparse
import fnmatch
import numpy as np

import scipy
from scipy.interpolate import interp1d

# Astropy
from astropy.io import fits
from astropy    import units as u
from astropy.stats import sigma_clip
# AstroML
from astroML.plotting import hist

# Matplotlib related
# 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

import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from matplotlib.ticker import NullFormatter
from matplotlib.ticker import MaxNLocator

# Personal
import hscUtils as hUtil
import galSBP

import cPickle as pickle

amag_sun_g = 5.33
amag_sun_r = 4.67
amag_sun_i = 4.48
amag_sun_z = 4.42


Populating the interactive namespace from numpy and matplotlib

In [ ]:
# Location of the cutouts 

""" On Astro3 / For MacOS"""
baseDir = '/Users/songhuang/astro3/hscs' 
redBCG = os.path.join(baseDir, 'redmapper')
nonBCG = os.path.join(baseDir, 'nonbcg')

A Few Useful Functions


In [2]:
def normProfile(rad, sbp, minRad, maxRad): 
    """
    Normalize the profile 
    """
    offset = np.nanmedian(sbp[(rad >= minRad) & (rad <= maxRad)])
    
    return (sbp-offset)

In [3]:
def logAdd(para1, para2):    
    """
    Doc
    """
    return np.log10((10.0 ** np.asarray(para1)) + (10.0 ** np.asarray(para2)))

In [4]:
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'))

In [5]:
def findProfile(pattern, loc, verbose=False):
    """
    Find the prefix of the ellipse profiles 
    """
    result = []
    for root, dirs, files in os.walk(loc):
        for name in files:
            if fnmatch.fnmatch(name, pattern):
                result.append(os.path.join(root, name))
                
    if verbose: 
        print "### %d files found !" % len(result)
    
    return result

In [6]:
def readProfile(ellFile):
    """
    Load the pickle format 1-D profile
    """
    if os.path.isfile(ellFile):
        return pickle.load(open(ellFile, 'rb'))
    else:
        raise Exception("XXX Can not find the Ellipse file")

In [7]:
def getEllipProfile(objid, loc, stage='3', verbose=False):
    
    pattern = '*' + str(objid) + '*_img_ellip_' + stage + '.pkl'
    
    result = findProfile(pattern, loc)
    
    if len(result) == 0: 
        #raise Exception('XXX Can not find the Ellipse profile ! %s' % pattern)
        if verbose:
            print 'XXX Can not find the Ellipse profile ! %s' % pattern
            return None
    elif len(result) == 1: 
        ellProf = readProfile(result[0])
        return ellProf
    else: 
        if verbose:
            print "XXX Warning: More than 1 profiles found! Use the first one : %s" % result[0]
        ellProf = readProfile(result[0])
        return ellProf

In [1]:
def pixKpc(redshift, pix=0.168, show=True, npix=1.0):
    """
    The corresponding Kpc size of a pixel 
    """
    
    pixKpc = pix * npix * hUtil.cosmoScale(redshift)
    
    if show:
        print "# %d pixel(s) = %6.3f Kpc" % (npix, pixKpc)
        
    return pixKpc

In [9]:
def correctProf(ellProf, redshift, verbose=False, extinction=0.0, zp=27.0, 
                amag_sun=None, corCurve=False):
    
    """
    DOC
    """
    
    scale = hUtil.cosmoScale(redshift)
    distmod = hUtil.cosmoDistMod(redshift)
    
    if verbose: 
        print "### REDSHIFT : %6.4f" % redshift
        print "### SCALE : %7.4f kpc/arcsec" % scale 
        print "### DISTMOD : %7.4f mag" % distmod
    
    sma_kpc = ellProf['sma_asec'] * scale
    if not corCurve:
        abs_mag = -2.5 * np.log10(ellProf['growth_ori']) + zp - extinction - distmod
    else:
        abs_mag = -2.5 * np.log10(ellProf['growth_cor']) + zp - extinction - distmod

    if amag_sun is not None: 
        abs_sbp = (amag_sun + 21.572 - (ellProf['sbp'] - extinction))/2.5 + 6.0
    else:
        abs_sbp = ellProf['sbp'] - extinction

    return sma_kpc, abs_sbp, abs_mag

In [10]:
def photoCompare(sample, zp=27.0, prefix='hsc_xmm_log11.6', amag_sun=amag_sun_i, 
                 pix=0.168, absMag=True, aStr='extinction_i', verbose=False,
                 magStr='cModelAbsMag_i', stage='3', plot=True, color='b', 
                 radius=None, corCurve=False, idStr='objid_dr12', 
                 massStr='logMass', zStr='z', extCorrect=True, ax=None, 
                 marker='o'):
    
    """
    DOC
    """
    
    """ Read in the sample catalog """
    catFile = prefix + '_' + sample + '.fits' 
    catData = fits.open(catFile)[1].data
    
    """ Useful columns """
    gal_id   = catData[idStr]
    gal_z    = catData[zStr]
    
    if extCorrect:
        gal_alam = catData[aStr]
    else: 
        gal_alam = gal_z * 0.0
        
    gal_mag  = catData[magStr]
    gal_logm = catData[massStr]

    """ The "total" magnitude from SBP """
    mag_sbp = [] 
    mag_sdss = []
    
    for ii in range(len(catData)): 
        
        """ Get the Ellipse profile """
        ellProf = getEllipProfile(gal_id[ii], sample, stage=stage)
        
        if ellProf is not None:
            
            """ Get the radius in Kpc, luminosity profile, and growth curve for absolute magnitude """
            sma_kpc, abs_sbp, abs_cog = correctProf(ellProf, gal_z[ii], extinction=gal_alam[ii], 
                                                    zp=27.0, amag_sun=amag_sun, corCurve=corCurve)
            
            """ Distance module """
            distmod = hUtil.cosmoDistMod(gal_z[ii])
            
            if radius is None: 
                mag_use = ellProf['mag_tot'][0] - gal_alam[ii]
                if absMag:
                    mag_use = mag_use - distmod
            else: 
                rsma_use = radius ** 0.25
                rsma_common = np.arange(0, 4.5, 0.1)
                
                indexUse = np.isfinite(abs_cog)
                intrpFunc = interp1d(sma_kpc[indexUse]**0.25, abs_cog[indexUse], kind='linear')
                try:
                    mag_use = intrpFunc(rsma_use) - gal_alam[ii]
                    if not absMag:
                        mag_use = mag_use + distmod
                except Exception:
                    print "### Interpolation range is bad : %d" % gal_id[ii]
                    mag_use = ellProf['mag_tot'][0] - gal_alam[ii]
                    if absMag:
                        mag_use = mag_use - distmod
                
            mag_sbp.append(mag_use) 
            mag_sdss.append(gal_mag[ii])
            
            if verbose:
                print "### %22d -- %6.3f -- %6.3f -- %6.3f" % (gal_id[ii], gal_mag[ii], 
                                                               ellProf['mag_tot'][0], distmod)
    mag_sbp = np.asarray(mag_sbp)
    mag_sdss = np.asarray(mag_sdss)
    
    if plot:
        if ax is not None: 
            ax.scatter(mag_sbp, mag_sdss, marker=marker, color=color)
        else:
            scatter(mag_sbp, mag_sdss, marker=marker, color=color)
        
    return mag_sbp

In [33]:
def getStackProfile(sample, zp=27.0, prefix='hsc_xmm_log11.6', amag_sun=amag_sun_i, 
                    pix=0.168, absMag=True, aStr='extinction_i', verbose=False,
                    magStr='cModelAbsMag_i', stage='3', plot=True, plotAvg=True, 
                    color='b', width=0.9, singleColor='k', 
                    radius=None, corCurve=False, useMedian=True, alpha=0.1,
                    minRsma=0.3, maxRsma=3.5, extCorrect=True, idStr='objid_dr12', 
                    zStr='z', massStr='logMass', rsmaStep=0.1, ax=None, root=None, 
                    idSuffix=None, dirSuffix=None, catFile=None):
    
    """
    DOC
    """
    
    """ Read in the sample catalog """
    if catFile is None:
        catFile = prefix + '_' + sample + '.fits' 
        
    if root is not None: 
        catFile = os.path.join(root, catFile)
        
    if not os.path.isfile(catFile):
        raise Exception('XXX Can not find the input catalog : %s', catFile)
    else:
        catData = fits.open(catFile)[1].data
    
    """ Useful columns """
    gal_id   = catData[idStr]
    gal_z    = catData[zStr]
    print "### Median Redshift : %5.3f" % np.nanmedian(gal_z)
    
    if extCorrect:
        try:
            gal_alam = catData[aStr]
        except Exception: 
            print "### Can not find the column : %s" % aStr 
            gal_alam = gal_z * 0.0
    else:
        gal_alam = gal_z * 0.0
    
    if magStr is not None:
        try:
            gal_mag = catData[magStr]
            print " ## Median Magnitude : %5.2f" % np.nanmedian(gal_mag)
        except Exception:
            print "### Can not find the column : %s" % magStr 
            gal_mag = gal_z * 0.0
    else: 
        gal_mag = gal_z * 0.0 
        
    if massStr is not None:
        try:
            gal_mass = catData[massStr]
            print " ## Median Log Stellar Mass : %5.2f" % np.nanmedian(gal_mass)
        except Exception: 
            print "### Can not find the column : %s" % massStr 
            gal_mass = gal_z * 0.0            
    else: 
        gal_mass = gal_z * 0.0
    
    """ Plot """
    for ii in range(len(catData)): 
        
        if root is not None: 
            profDir = os.path.join(root, sample)
        else: 
            profDir = sample 
            
        if idSuffix is None: 
            pattern = str(gal_id[ii])
        else: 
            pattern = str(gal_id[ii]) + '*' + idSuffix 
            
        if dirSuffix is None: 
            location = profDir
        else: 
            location = profDir + '_' + dirSuffix
            
        ellProf = getEllipProfile(pattern, location, stage=stage)
    
        if ellProf is not None:
            sma_kpc, abs_sbp, abs_cog = correctProf(ellProf, gal_z[ii], 
                                                    extinction=gal_alam[ii], 
                                                    zp=27.0, verbose=False, 
                                                    amag_sun=amag_sun, 
                                                    corCurve=corCurve)
                    
            if plot:
                if ax is not None:
                    ax.plot(sma_kpc**0.25, abs_sbp, singleColor, alpha=alpha, linewidth=width)
                else: 
                    plot(sma_kpc**0.25, abs_sbp, singleColor, alpha=alpha, linewidth=width)  
                        
            intrpFunc = interp1d(sma_kpc**0.25, abs_sbp, 
                                 kind='linear')
            try:
                rsma_common = np.arange(minRsma, maxRsma, rsmaStep)
                asbp_intrp = intrpFunc(rsma_common)
                
                try: 
                    asbp_stack = np.vstack((asbp_stack, asbp_intrp))
                except NameError: 
                    asbp_stack = asbp_intrp
                        
            except Exception as errMsg:
                print errMsg
                print "### Bad interpolation range: %d " % gal_id[ii]
                print "    %5.2f -- %5.2f" % (np.nanmin(sma_kpc**0.25), 
                                              np.nanmax(sma_kpc**0.25))
    
    if useMedian:
        asbp_avg = np.nanmedian(asbp_stack, axis=0) 
    else:
        asbp_avg = np.nanmean(asbp_stack, axis=0) 

    if plotAvg:
        if ax is not None:
            print len(rsma_common)
            print len(asbp_avg)
            ax.plot(rsma_common, asbp_avg, color, linewidth=3.5)
        else:
            plot(rsma_common, asbp_avg, color, linewidth=3.5)
    
    return rsma_common, asbp_avg, asbp_stack

Playground


In [38]:
#def photoCompare(sample, zp=27.0, prefix='hsc_xmm_log11.6', amag_sun=amag_sun_i, 
#                 pix=0.168, absMag=True, aStr='extinction_i', verbose=False,
#                 magStr='cModelAbsMag_i', stage='3', plot=True, color='b', 
#                 radius=None, corCurve=False, idStr='objid_dr12', 
#                 massStr='logMass', zStr='z', extCorrect=True, ax=None, 
#                 marker='o'):

fig = plt.figure(figsize=(10, 10))
fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
ax1 = fig.add_subplot(111)

maxR = 100.0

print "### Z1"
mag_sbp_z1 = photoCompare('z1', color='k', radius=maxR, ax=ax1)
print "### Z2"
mag_sbp_z2 = photoCompare('z2', color='r', radius=maxR, ax=ax1)
#mag_sbp_z2_b = photoCompare('z2', color='b', radius=100.0)
print "### Z3"
mag_sbp_z3 = photoCompare('z3', color='g', radius=maxR, ax=ax1)
print "### Z4"
mag_sbp_z4 = photoCompare('z4', color='b', radius=maxR, ax=ax1)

ax1.set_xlim(-21.5, -24.5)
ax1.set_ylim(-21.5, -24.5)
ax1.set_xlabel('Absolute Magnitude (HSC)',  size=16.0)
ax1.set_ylabel('Absolute Magnitude (SDSS)', size=16.0)

magRange = np.arange(-25.0, -20.0, 0.2)
ax1.plot(magRange, magRange, 'k')


### Z1
### Z2
### Z3
/usr/local/lib/python2.7/site-packages/IPython/kernel/__main__.py:17: RuntimeWarning: divide by zero encountered in log10
/usr/local/lib/python2.7/site-packages/IPython/kernel/__main__.py:17: RuntimeWarning: invalid value encountered in log10
### Z4
Out[38]:
[<matplotlib.lines.Line2D at 0x112375910>]

In [43]:
fig = plt.figure(figsize=(10, 6))
fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
ax1 = fig.add_subplot(111)

hist(mag_sbp_z4, bins='knuth', histtype='stepfilled', alpha=0.5, 
     normed=True, color='b')
hist(mag_sbp_z3, bins='knuth', histtype='stepfilled', alpha=0.3, 
     normed=True, color='g')
hist(mag_sbp_z2, bins='knuth', histtype='step', alpha=0.8, 
     normed=True, color='r', linewidth=3.0)
hist(mag_sbp_z1, bins='knuth', histtype='step', alpha=0.8, 
     normed=True, color='k', linewidth=3.0)

ax1.set_xlabel('Absolute Magnitude', size=19.0)
ax1.set_ylabel('Fraction', size=19.0)
ax1.set_xlim(-21.5, -24.5)

for tick in ax1.xaxis.get_major_ticks():
    tick.label.set_fontsize(19) 
for tick in ax1.yaxis.get_major_ticks():
    tick.label.set_fontsize(19)


Out[43]:
(-21.5, -24.5)

In [544]:
#def photoCompare(sample, zp=27.0, prefix='hsc_xmm_log11.6', amag_sun=amag_sun_i, 
#                 pix=0.168, absMag=True, aStr='extinction_i', verbose=False,
#                 magStr='cModelAbsMag_i', stage='3', plot=True, color='b', 
#                 radius=None, corCurve=False, idStr='objid_dr12', 
#                 massStr='logMass', zStr='z', extCorrect=True, ax=None, 
#                 marker='o'):


R1 = 20.0
R2 = 50.0 
R3 = 100.0
R4 = 120.0

sample = 'z2'

print "### R1"
mag_sbp_r1 = photoCompare(sample, radius=R1, plot=None)
print "### R2"
mag_sbp_r2 = photoCompare(sample, radius=R2, plot=None)
print "### R3"
mag_sbp_r3 = photoCompare(sample, radius=R3, plot=None)
#print "### R4"
#mag_sbp_r4 = photoCompare(sample, radius=R4, plot=None)


### R1
### R2
### R3
/usr/local/lib/python2.7/site-packages/IPython/kernel/__main__.py:18: RuntimeWarning: divide by zero encountered in log10

In [555]:
fig = plt.figure(figsize=(10, 10))
fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
ax1 = fig.add_subplot(111)

magRange = np.arange(-25.0, -20.0, 0.2)
ax1.plot(magRange, magRange, 'k', linestyle='--', linewidth=3.0, alpha=0.8)

ax1.scatter(mag_sbp_r1, mag_sbp_r2, marker='o', color='g', s=75, 
           alpha=0.7)
ax1.scatter(mag_sbp_r1, mag_sbp_r3, marker='o', color='r', s=75,
           alpha=0.5)
#ax1.scatter(mag_sbp_r1, mag_sbp_r4, marker='o', color='b')

ax1.set_xlim(-22.2, -24.3)
ax1.set_ylim(-22.2, -24.3)

ax1.set_xlabel('Absolute Magnitude (within 20 Kpc)',  size=25.0)
ax1.set_ylabel('Absolute Magnitude (within larger aperture)', size=25.0)

for tick in ax1.xaxis.get_major_ticks():
    tick.label.set_fontsize(20) 
for tick in ax1.yaxis.get_major_ticks():
    tick.label.set_fontsize(20) 
    
ax1.spines['top'].set_linewidth(3.5)
ax1.spines['right'].set_linewidth(3.5)
ax1.spines['bottom'].set_linewidth(3.5)
ax1.spines['left'].set_linewidth(3.5)

ax1.text(-23.45, -23.0, '$0.3<z<0.4$', size=32.0)
ax1.text(-23.25, -22.8, 'log $(M_{s}/M_{\odot})>11.6$', size=32.0)

ax1.text(-23.35, -22.55, 'within 50 Kpc', color='g', size=32.0)
ax1.text(-23.30, -22.35, 'within 100 Kpc', color='r', size=32.0)


Out[555]:
<matplotlib.text.Text at 0x114d52e50>

In [559]:
fig = plt.figure(figsize=(12, 7.5))
fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
ax1 = fig.add_subplot(111)

hist(mag_sbp_r1, bins='knuth', histtype='stepfilled', alpha=0.5, 
     normed=True, color='b')
hist(mag_sbp_r2, bins='knuth', histtype='stepfilled', alpha=0.3, 
     normed=True, color='g')
hist(mag_sbp_r3, bins='knuth', histtype='stepfilled', alpha=0.6, 
     normed=True, color='r', linewidth=3.0)
#hist(mag_sbp_r4, bins='knuth', histtype='step', alpha=0.8, 
#     normed=True, color='k', linewidth=3.0)

ax1.set_xlabel('Magnitude', size=30.0)
ax1.set_ylabel('Fraction', size=30.0)
ax1.set_xlim(-22.01, -24.4)
ax1.set_ylim(0.01, 1.45)

for tick in ax1.xaxis.get_major_ticks():
    tick.label.set_fontsize(28) 
for tick in ax1.yaxis.get_major_ticks():
    tick.label.set_fontsize(28) 
    
ax1.spines['top'].set_linewidth(3.5)
ax1.spines['right'].set_linewidth(3.5)
ax1.spines['bottom'].set_linewidth(3.5)
ax1.spines['left'].set_linewidth(3.5)

ax1.text(-22.1, 1.2, '$<20$ Kpc', size=32, color='b')
ax1.text(-22.1, 1.0, '$<50$ Kpc', size=32, color='g')
ax1.text(-22.1, 0.8, '$<100$ Kpc', size=32, color='r')


Out[559]:
<matplotlib.text.Text at 0x1191abf90>

Average Profile


In [36]:
fig = plt.figure(figsize=(12, 12))
fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
ax1 = fig.add_subplot(111)

rsma_common_z1, asbp_avg_z1, asbp_stack_z1 = getStackProfile('z1', minRsma=0.1, maxRsma=3.8,
                                                             ax=ax1, alpha=0.1)

rsma_common_z1m, asbp_avg_z1m, asbp_stack_z1m = getStackProfile('z1', prefix='hsc_xmm_mem', minRsma=0.1, 
                                                                maxRsma=3.8, ax=ax1, alpha=0.7, plotAvg=False,
                                                                width=3.0, singleColor='g')
rsma_common_z1b, asbp_avg_z1b, asbp_stack_z1b = getStackProfile('z1', prefix='hsc_xmm_bcg', minRsma=0.1, 
                                                                maxRsma=3.8, ax=ax1, alpha=0.95, plotAvg=False, 
                                                                width=3.0, singleColor='r')

print asbp_stack_z1.shape

std_stack_z1 = np.nanstd(asbp_stack_z1, axis=0)

#ax1.fill_between(rsma_common_z1, asbp_avg_z1+std_stack_z1, asbp_avg_z1-std_stack_z1,
#                 facecolor='b', alpha=0.2)

ax1.text(2.4, 8.6, '$0.15<z<0.30$', size=42)

ax1.set_xlabel('$R^{1/4}$ (Kpc)', size=32)
ax1.set_ylabel('log $(\mu/[L_{\odot} Kpc^{-2}])$', size=32)

ax1.set_xlim(0.0, 4.2)
ax1.set_ylim(3.2, 9.5)

for tick in ax1.xaxis.get_major_ticks():
    tick.label.set_fontsize(30) 
for tick in ax1.yaxis.get_major_ticks():
    tick.label.set_fontsize(30) 
    
ax1.spines['top'].set_linewidth(3.5)
ax1.spines['right'].set_linewidth(3.5)
ax1.spines['bottom'].set_linewidth(3.5)
ax1.spines['left'].set_linewidth(3.5)


### Median Redshift : 0.283
 ## Median Magnitude : -22.30
 ## Median Log Stellar Mass : 11.71
37
37
### Median Redshift : 0.290
 ## Median Magnitude : -22.62
 ## Median Log Stellar Mass : 11.75
### Median Redshift : 0.281
 ## Median Magnitude : -22.85
 ## Median Log Stellar Mass : 11.81
(68, 37)
/usr/local/lib/python2.7/site-packages/IPython/kernel/__main__.py:18: RuntimeWarning: divide by zero encountered in log10

In [30]:
fig = plt.figure(figsize=(12, 12))
fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
ax1 = fig.add_subplot(111)

rsma_common_z2, asbp_avg_z2, asbp_stack_z2 = getStackProfile('z2', minRsma=0.1, maxRsma=3.8,
                                                             ax=ax1, alpha=0.1)
rsma_common_z2m, asbp_avg_z2m, asbp_stack_z2m = getStackProfile('z2', prefix='hsc_xmm_mem', minRsma=0.1, 
                                                                maxRsma=3.8, ax=ax1, alpha=0.5, plotAvg=False,
                                                                width=3.0, singleColor='g')
rsma_common_z2b, asbp_avg_z2b, asbp_stack_z2b = getStackProfile('z2', prefix='hsc_xmm_bcg', minRsma=0.1, 
                                                                maxRsma=3.8, ax=ax1, alpha=0.9, plotAvg=False,
                                                                width=3.0, singleColor='r')

print asbp_stack_z2.shape

std_stack_z2 = np.nanstd(asbp_stack_z2, axis=0)
#ax1.fill_between(rsma_common_z2, asbp_avg_z2+std_stack_z2, asbp_avg_z2-std_stack_z2,
#                 facecolor='b', alpha=0.2)

ax1.text(2.4, 8.6, '$0.30<z<0.40$', size=42)

ax1.set_xlabel('$R^{1/4}$ (Kpc)', size=32)
ax1.set_ylabel('log $(\mu/[L_{\odot} Kpc^{-2}])$', size=32)

ax1.set_xlim(0.0, 4.2)
ax1.set_ylim(3.2, 9.5)

for tick in ax1.xaxis.get_major_ticks():
    tick.label.set_fontsize(30) 
for tick in ax1.yaxis.get_major_ticks():
    tick.label.set_fontsize(30) 

ax1.spines['top'].set_linewidth(3.5)
ax1.spines['right'].set_linewidth(3.5)
ax1.spines['bottom'].set_linewidth(3.5)
ax1.spines['left'].set_linewidth(3.5)


### Median Redshift : 0.334
 ## Median Magnitude : -22.37
 ## Median Log Stellar Mass : 11.73
37
37
### Median Redshift : 0.333
 ## Median Magnitude : -22.35
 ## Median Log Stellar Mass : 11.72
### Median Redshift : 0.382
 ## Median Magnitude : -22.49
 ## Median Log Stellar Mass : 11.74
(89, 37)
/usr/local/lib/python2.7/site-packages/IPython/kernel/__main__.py:18: RuntimeWarning: divide by zero encountered in log10

In [34]:
fig = plt.figure(figsize=(12, 12))
fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
ax1 = fig.add_subplot(111)

rsma_common_z3, asbp_avg_z3, asbp_stack_z3 = getStackProfile('z3', minRsma=0.1, maxRsma=3.8,
                                                             ax=ax1, alpha=0.03)
rsma_common_z3m, asbp_avg_z3m, asbp_stack_z3m = getStackProfile('z3', prefix='hsc_xmm_mem', minRsma=0.1, 
                                                                maxRsma=3.6, ax=ax1, alpha=0.5, plotAvg=False,
                                                                width=3.0, singleColor='g')
rsma_common_z3b, asbp_avg_z3b, asbp_stack_z3b = getStackProfile('z3', prefix='hsc_xmm_bcg', minRsma=0.1, 
                                                                maxRsma=3.6, ax=ax1, alpha=0.9, plotAvg=False,
                                                                width=3.0, singleColor='r')

print asbp_stack_z3.shape

std_stack_z3 = np.nanstd(asbp_stack_z3, axis=0)
#ax1.fill_between(rsma_common_z3, asbp_avg_z3+std_stack_z3, asbp_avg_z3-std_stack_z3,
#                 facecolor='b', alpha=0.2)

ax1.text(2.4, 8.6, '$0.40<z<0.50$', size=42)

ax1.set_xlabel('$R^{1/4}$ (Kpc)', size=32)
ax1.set_ylabel('log $(\mu/[L_{\odot} Kpc^{-2}])$', size=32)

ax1.set_xlim(0.0, 4.2)
ax1.set_ylim(3.2, 9.5)

for tick in ax1.xaxis.get_major_ticks():
    tick.label.set_fontsize(30) 
for tick in ax1.yaxis.get_major_ticks():
    tick.label.set_fontsize(30) 
    
ax1.spines['top'].set_linewidth(3.5)
ax1.spines['right'].set_linewidth(3.5)
ax1.spines['bottom'].set_linewidth(3.5)
ax1.spines['left'].set_linewidth(3.5)


/usr/local/lib/python2.7/site-packages/IPython/kernel/__main__.py:18: RuntimeWarning: divide by zero encountered in log10
/usr/local/lib/python2.7/site-packages/IPython/kernel/__main__.py:18: RuntimeWarning: invalid value encountered in log10
### Median Redshift : 0.458
 ## Median Magnitude : -22.50
 ## Median Log Stellar Mass : 11.76
37
37
### Median Redshift : 0.447
 ## Median Magnitude : -23.01
 ## Median Log Stellar Mass : 11.93
### Median Redshift : 0.446
 ## Median Magnitude : -23.09
 ## Median Log Stellar Mass : 11.97
(176, 37)