<jupyter_start><jupyter_text>Recent Results:<jupyter_code>newDir = '/Users/songhuang/work/hscs/gama_massive/sbp/'

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')
gamaDir = os.path.join(newDir, 'gama')

# Two summary catalogs
#bcgCat = os.path.join(bcgDir, 'redmapper_bcg_hscmatch_mass_use_sbpsum_modA_muI1.fits')
bcgCat = os.path.join(bcgDir, 'redbcg_mass_use_dom.fits')
memCat = os.path.join(memDir, 'redmapper_mem_hscmatch_mass_sbpsum_modA_muI1.fits')
gamaCat = os.path.join(gamaDir, 'gama_massive_160107_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(gamaCat):
    raise Exception("## Can not find catalog for GAMA galaxies : %s" % gamaCat)
else: 
    gamaTab = Table.read(gamaCat, format='fits')<jupyter_output><empty_output><jupyter_text>Read in the catalogs for certain subsample:  M1a<jupyter_code>suffixT = 'M1a'
suffixP = 'm1a'
usePcen = True

if usePcen: 
    bcgSample1 = Table.read(os.path.join(bcgDir, 'bcg_' + suffixT + '_pcen.fits'), format='fits')
else: 
    bcgSample1 = Table.read(os.path.join(bcgDir, 'bcg_' + suffixT + '.fits'), format='fits')
    
memSample1 = Table.read(os.path.join(memDir, 'mem_' + suffixT + '.fits'), format='fits')
gamaSample1 = Table.read(os.path.join(gamaDir, 'gama_' + suffixT + '.fits'), format='fits')
if usePcen:
    bcgPkl1 = os.path.join(bcgDir, 'massive_bcg_' + suffixP + '_p_profs.pkl')
else: 
    bcgPkl1 = os.path.join(bcgDir, 'massive_bcg_' + suffixP + '_profs.pkl')

if os.path.isfile(bcgPkl1):
    print(bcgPkl1)
    print("# Read in available stacks of BCG/%s" % suffixT)
    bcgProfs1 = loadPkl(bcgPkl1)
else:
    bcgProfs1 = getStackProfiles(bcgSample1, bcgDir, name=('massive_bcg_' + suffixP))

print("## Dealing with %d profiles" % len(bcgProfs1))
memPkl1 = os.path.join(memDir, 'massive_mem_' + suffixP + '_profs.pkl')

if os.path.isfile(memPkl1):
    print(memPkl1)
    print("# Read in available stacks of BCG/%s" % suffixT)
    memProfs1 = loadPkl(memPkl1)
else:
    memProfs1 = getStackProfiles(memSample1, memDir, name=('massive_mem_' + suffixP))

print("## Dealing with %d profiles" % len(memProfs1))
gamaPkl1 = os.path.join(gamaDir, 'massive_gama_' + suffixP + '_profs.pkl')

if os.path.isfile(gamaPkl1):
    print(gamaPkl1)
    print("# Read in available stacks of BCG/%s" % suffixT)
    gamaProfs1 = loadPkl(gamaPkl1)
else:
    gamaProfs1 = getStackProfiles(gamaSample1, gamaDir, name=('massive_gama_' + suffixP))

print("## Dealing with %d profiles" % len(gamaProfs1))<jupyter_output>/Users/songhuang/work/hscs/gama_massive/sbp/gama/massive_gama_m1a_profs.pkl
# Read in available stacks of BCG/M1a
## Dealing with 1053 profiles<jupyter_text>M2a<jupyter_code>suffixT = 'M2a'
suffixP = 'm2a'
usePcen = True

if usePcen: 
    bcgSample2 = Table.read(os.path.join(bcgDir, 'bcg_' + suffixT + '_pcen.fits'), format='fits')
else: 
    bcgSample2 = Table.read(os.path.join(bcgDir, 'bcg_' + suffixT + '.fits'), format='fits')
    
memSample2 = Table.read(os.path.join(memDir, 'mem_' + suffixT + '.fits'), format='fits')
gamaSample2 = Table.read(os.path.join(gamaDir, 'gama_' + suffixT + '.fits'), format='fits')
if usePcen:
    bcgPkl2 = os.path.join(bcgDir, 'massive_bcg_' + suffixP + '_p_profs.pkl')
else: 
    bcgPkl2 = os.path.join(bcgDir, 'massive_bcg_' + suffixP + '_profs.pkl')

if os.path.isfile(bcgPkl2):
    print(bcgPkl2)
    print("# Read in available stacks of BCG/%s" % suffixT)
    bcgProfs2 = loadPkl(bcgPkl2)
else:
    bcgProfs2 = getStackProfiles(bcgSample2, bcgDir, name=('massive_bcg_' + suffixP))

print("## Dealing with %d profiles" % len(bcgProfs2))
memPkl2 = os.path.join(memDir, 'massive_mem_' + suffixP + '_profs.pkl')

if os.path.isfile(memPkl2):
    print(memPkl2)
    print("# Read in available stacks of BCG/%s" % suffixT)
    memProfs2 = loadPkl(memPkl2)
else:
    memProfs2 = getStackProfiles(memSample2, memDir, name=('massive_mem_' + suffixP))

print("## Dealing with %d profiles" % len(memProfs2))
gamaPkl2 = os.path.join(gamaDir, 'massive_gama_' + suffixP + '_profs.pkl')

if os.path.isfile(gamaPkl2):
    print(gamaPkl2)
    print("# Read in available stacks of BCG/%s" % suffixT)
    gamaProfs2 = loadPkl(gamaPkl2)
else:
    gamaProfs2 = getStackProfiles(gamaSample2, gamaDir, name=('massive_gama_' + suffixP))

print("## Dealing with %d profiles" % len(gamaProfs2))<jupyter_output>/Users/songhuang/work/hscs/gama_massive/sbp/gama/massive_gama_m2a_profs.pkl
# Read in available stacks of BCG/M2a
## Dealing with 304 profiles<jupyter_text>M3a<jupyter_code>suffixT = 'M3a'
suffixP = 'm3a'
usePcen = True

if usePcen: 
    bcgSample3 = Table.read(os.path.join(bcgDir, 'bcg_' + suffixT + '_pcen.fits'), format='fits')
else: 
    bcgSample3 = Table.read(os.path.join(bcgDir, 'bcg_' + suffixT + '.fits'), format='fits')
    
memSample3 = Table.read(os.path.join(memDir, 'mem_' + suffixT + '.fits'), format='fits')
gamaSample3 = Table.read(os.path.join(gamaDir, 'gama_' + suffixT + '.fits'), format='fits')
if usePcen:
    bcgPkl3 = os.path.join(bcgDir, 'massive_bcg_' + suffixP + '_p_profs.pkl')
else: 
    bcgPkl3 = os.path.join(bcgDir, 'massive_bcg_' + suffixP + '_profs.pkl')

if os.path.isfile(bcgPkl3):
    print(bcgPkl3)
    print("# Read in available stacks of BCG/%s" % suffixT)
    bcgProfs3 = loadPkl(bcgPkl3)
else:
    bcgProfs3 = getStackProfiles(bcgSample3, bcgDir, name=('massive_bcg_' + suffixP))

print("## Dealing with %d profiles" % len(bcgProfs3))
memPkl3 = os.path.join(memDir, 'massive_mem_' + suffixP + '_profs.pkl')

if os.path.isfile(memPkl3):
    print(memPkl3)
    print("# Read in available stacks of BCG/%s" % suffixT)
    memProfs3 = loadPkl(memPkl3)
else:
    memProfs3 = getStackProfiles(memSample3, memDir, name=('massive_mem_' + suffixP))

print("## Dealing with %d profiles" % len(memProfs3))
gamaPkl3 = os.path.join(gamaDir, 'massive_gama_' + suffixP + '_profs.pkl')

if os.path.isfile(gamaPkl3):
    print(gamaPkl3)
    print("# Read in available stacks of BCG/%s" % suffixT)
    gamaProfs3 = loadPkl(gamaPkl3)
else:
    gamaProfs3 = getStackProfiles(gamaSample3, gamaDir, name=('massive_gama_' + suffixP))

print("## Dealing with %d profiles" % len(gamaProfs3))<jupyter_output>/Users/songhuang/work/hscs/gama_massive/sbp/gama/massive_gama_m3a_profs.pkl
# Read in available stacks of BCG/M3a
## Dealing with 53 profiles<jupyter_text>Test the R50, R90 measurements:<jupyter_code>gamaR50_1 = np.log10(np.asarray([getFracRadius(pp, frac=0.5, maxRad=160.0,
                                               lum='lumI1', returnRatio=False) for pp in gamaProfs1]))
gamaR90_1 = np.log10(np.asarray([getFracRadius(pp, frac=0.9, maxRad=160.0, 
                                               lum='lumI1', returnRatio=False) for pp in gamaProfs1]))
gamaM100_1 = np.asarray([(pp.meta['LUM_100'] + pp.meta['LOGM2LI_C']) for pp in gamaProfs1])

gamaR50_2 = np.log10(np.asarray([getFracRadius(pp, frac=0.5, maxRad=160.0,
                                               lum='lumI1', returnRatio=False) for pp in gamaProfs2]))
gamaR90_2 = np.log10(np.asarray([getFracRadius(pp, frac=0.9, maxRad=160.0, 
                                               lum='lumI1', returnRatio=False) for pp in gamaProfs2]))
gamaM100_2 = np.asarray([(pp.meta['LUM_100'] + pp.meta['LOGM2LI_C']) for pp in gamaProfs2])

gamaR50_3 = np.log10(np.asarray([getFracRadius(pp, frac=0.5, maxRad=160.0,
                                               lum='lumI1', returnRatio=False) for pp in gamaProfs3]))
gamaR90_3 = np.log10(np.asarray([getFracRadius(pp, frac=0.9, maxRad=160.0, 
                                               lum='lumI1', returnRatio=False) for pp in gamaProfs3]))
gamaM100_3 = np.asarray([(pp.meta['LUM_100'] + pp.meta['LOGM2LI_C']) for pp in gamaProfs3])
bcgR50_1 = np.log10(np.asarray([getFracRadius(pp, frac=0.5, maxRad=160.0, 
                                              lum='lumI1', returnRatio=False) for pp in bcgProfs1]))
bcgR90_1 = np.log10(np.asarray([getFracRadius(pp, frac=0.9, maxRad=160.0, 
                                              lum='lumI1', returnRatio=False) for pp in bcgProfs1]))
bcgM100_1 = np.asarray([(pp.meta['LUM_100'] + pp.meta['LOGM2LI_C']) for pp in bcgProfs1])

bcgR50_2 = np.log10(np.asarray([getFracRadius(pp, frac=0.5, maxRad=160.0, 
                                              lum='lumI1', returnRatio=False) for pp in bcgProfs2]))
bcgR90_2 = np.log10(np.asarray([getFracRadius(pp, frac=0.9, maxRad=160.0, 
                                              lum='lumI1', returnRatio=False) for pp in bcgProfs2]))
bcgM100_2 = np.asarray([(pp.meta['LUM_100'] + pp.meta['LOGM2LI_C']) for pp in bcgProfs2])

bcgR50_3 = np.log10(np.asarray([getFracRadius(pp, frac=0.5, maxRad=160.0, 
                                              lum='lumI1', returnRatio=False) for pp in bcgProfs3]))
bcgR90_3 = np.log10(np.asarray([getFracRadius(pp, frac=0.9, maxRad=160.0, 
                                              lum='lumI1', returnRatio=False) for pp in bcgProfs3]))
bcgM100_3 = np.asarray([(pp.meta['LUM_100'] + pp.meta['LOGM2LI_C']) for pp in bcgProfs3])
memR50_1 = np.log10(np.asarray([getFracRadius(pp, frac=0.5, maxRad=160.0, 
                                              lum='lumI1', returnRatio=False) for pp in memProfs1]))
memR90_1 = np.log10(np.asarray([getFracRadius(pp, frac=0.9, maxRad=160.0, 
                                              lum='lumI1', returnRatio=False) for pp in memProfs1]))
memM100_1 = np.asarray([(pp.meta['LUM_100'] + pp.meta['LOGM2LI_C']) for pp in memProfs1])

memR50_2 = np.log10(np.asarray([getFracRadius(pp, frac=0.5, maxRad=160.0, 
                                              lum='lumI1', returnRatio=False) for pp in memProfs2]))
memR90_2 = np.log10(np.asarray([getFracRadius(pp, frac=0.9, maxRad=160.0, 
                                              lum='lumI1', returnRatio=False) for pp in memProfs2]))
memM100_2 = np.asarray([(pp.meta['LUM_100'] + pp.meta['LOGM2LI_C']) for pp in memProfs2])

memR50_3 = np.log10(np.asarray([getFracRadius(pp, frac=0.5, maxRad=160.0, 
                                              lum='lumI1', returnRatio=False) for pp in memProfs3]))
memR90_3 = np.log10(np.asarray([getFracRadius(pp, frac=0.9, maxRad=160.0, 
                                              lum='lumI1', returnRatio=False) for pp in memProfs3]))
memM100_3 = np.asarray([(pp.meta['LUM_100'] + pp.meta['LOGM2LI_C']) for pp in memProfs3])
print(np.nanmedian(gamaR90_1), np.nanmedian(bcgR90_1))
print(np.nanmedian(gamaR90_2), np.nanmedian(bcgR90_2))
print(np.nanmedian(gamaR90_3), np.nanmedian(bcgR90_3))
logmGama = np.hstack([gamaM100_1, gamaM100_2, gamaM100_3]).ravel()
logr50Gama = np.hstack([gamaR50_1, gamaR50_2, gamaR50_3]).ravel()
logr90Gama = np.hstack([gamaR90_1, gamaR90_2, gamaR90_3]).ravel()

logmBcg = np.hstack([bcgM100_1, bcgM100_2, bcgM100_3]).ravel()
logr50Bcg = np.hstack([bcgR50_1, bcgR50_2, bcgR50_3]).ravel()
logr90Bcg = np.hstack([bcgR90_1, bcgR90_2, bcgR90_3]).ravel()
logmGamaU = logmGama[logmGama >= 11.42]
logr50GamaU = logr50Gama[logmGama >= 11.42]
logr90GamaU = logr90Gama[logmGama >= 11.42]

logmBcgU = logmBcg[logmBcg >= 11.42]
logr50BcgU = logr50Bcg[logmBcg >= 11.42]
logr90BcgU = logr90Bcg[logmBcg >= 11.42]
logmGamaM = logmGama[logmGama >= 11.62]
logr50GamaM = logr50Gama[logmGama >= 11.62]
logr90GamaM = logr90Gama[logmGama >= 11.62]

logmBcgM = logmBcg[logmBcg >= 11.62]
logr50BcgM = logr50Bcg[logmBcg >= 11.62]
logr90BcgM = logr90Bcg[logmBcg >= 11.62]
plt.xlim(10.5, 12.6)
plt.ylim(0.05, 2.59)

plt.xlabel('$\log\ (M_{\star}/M_{\odot})$', fontsize=35)
plt.ylabel('$\log\ (R_{\mathrm{e}}/\mathrm{Kpc})$', fontsize=35)

plt.scatter(logmGamaU, logr50GamaU, c='k', marker='o', alpha=0.2)
plt.scatter(logmBcgU, logr50BcgU, c='r', marker='h', s=70, alpha=0.7)

plt.show()<jupyter_output><empty_output><jupyter_text>Fit Mass-Size Relation<jupyter_code>import emcee
import corner

# Likelihood function
def lnlike(theta, x, y, yerr):
    a, b, lnf = theta
    # Straightline model 
    model = a * x + b
    inv_sigma2 = (1.0 / (yerr**2 + model**2*np.exp(2*lnf)))
    return -0.5*(np.sum((y-model)**2*inv_sigma2 - 
                         np.log(inv_sigma2)))

# Priors
def lnprior(theta):
    a, b, lnf = theta
    if 0.2 < a < 1.4 and -13.0 < b < 0.0 and -3.0 < lnf < 3.0:
        return 0.0
    return -np.inf

def lnprob(theta, x, y, yerr):
    lp = lnprior(theta)
    if not np.isfinite(lp):
        return -np.inf
    return lp + lnlike(theta, x, y, yerr)

#Llinear least squares solution 
def llsLine(x, y, yerr):
    """ Simple straingt line fitting """
    A = np.vstack((np.ones_like(x), x)).T
    C = np.diag(yerr * yerr)
    cov = np.linalg.inv(np.dot(A.T, np.linalg.solve(C, A)))
    b_ls, a_ls = np.dot(cov, np.dot(A.T, 
                                    np.linalg.solve(C, y)))
    print("LLS: a =%8.5f ; b =%8.5f" % (a_ls, b_ls))
    return a_ls, b_ls
    
# Use Emcee to fit a straight line 
def emceeLine(x, y, yerr, nwalkers=100, ndim=3, nburn=100, 
              nstep=600, show=True): 
    """ Initial guesses from simple LLS fitting """
    #a_ls, b_ls = llsLine(x, y, yerr)
    a_ls, b_ls = 0.9, -10.5
    initial = [a_ls, b_ls, 0.00]
    """ Start the sampler """
    sampler = emcee.EnsembleSampler(nwalkers, ndim, 
                                    lnprob, args=(x, y, yerr))
    """ Initializing the walkers. """
    np.random.seed(0)
    guesses = [initial + (1e-2*np.random.randn(ndim)) 
               for i in range(nwalkers)]
    """ Run MCMC """
    print("Start the MCMC runs")
    %time sampler.run_mcmc(guesses, nstep)
    print("Done")
    """ Flatten the chain so that we have a flat list of samples """
    samples = sampler.chain[:, nburn:, :].reshape(-1, ndim)
    if show:   
        fig = corner.corner(samples, 
                            labels=["$a$", "$b$", "$\ln\,f$"])
    """ Compute the quantiles. """
    samples[:, 2] = np.exp(samples[:, 2])
    a_mcmc, b_mcmc, f_mcmc = map(lambda v: (v[1], v[2]-v[1], v[1]-v[0]),
                                 zip(*np.percentile(samples, [16, 50, 84],
                                                    axis=0)))
    print("""MCMC result:
             a = {0[0]} +{0[1]} -{0[2]}
             b = {1[0]} +{1[1]} -{1[2]} 
             f = {2[0]} +{2[1]} -{2[2]}
            """.format(a_mcmc, b_mcmc, f_mcmc))

    return a_mcmc, b_mcmc, f_mcmc
logr50BcgErr = logr50BcgM * 0.0 + 0.02

a_mcmc, b_mcmc, f_mcmc = emceeLine(logmBcgM, logr50BcgM, logr50BcgErr, 
                                   nwalkers=300, ndim=3, nburn=2000, 
                                   nstep=8000, show=True)

# Best paramters
aBcgB, aBcgU, aBcgL = a_mcmc
bBcgB, bBcgU, bBcgL = b_mcmc

plt.show()
plt.xlim(11.4, 12.6)
plt.ylim(0.45, 2.59)

plt.xlabel('$\log\ (M_{\star}/M_{\odot})$', fontsize=35)
plt.ylabel('$\log\ (R_{\mathrm{50}}/\mathrm{Kpc})$', fontsize=35)

plt.scatter(logmGamaU, logr50GamaU, c='k', marker='o', alpha=0.2)
plt.scatter(logmBcgU, logr50BcgU, c='r', marker='h', s=70, alpha=0.7)

# Red-sequence
xx = np.linspace(10.0, 13.0, 100)

# MCMC result:
# a = 0.960387107648 +0.0755428582291 -0.0750014260391
# b = -10.0134319999 +0.883141848718 -0.889418975333 
# f = 0.168797744636 +0.0104793857322 -0.00963780589871
plt.plot(xx, (0.96038 * xx - 10.01343),
        linestyle='--', color='r', linewidth=3.5)
#yy1 = (aBcgB + aBcgU) * xx + (bBcgB)
#yy2 = (aBcgB - aBcgL) * xx + (bBcgB)
#plt.fill_between(xx, yy1, yy2, 
#                facecolor='r', interpolate=True, alpha=0.15)

# MCMC result:
# a = 0.935475557486 +0.026570857285 -0.026374388974
# b = -9.81576469101 +0.305357355175 -0.307325596592 
# f = 0.15056796201 +0.00316298933649 -0.00311714344817
plt.plot(xx, (0.93546 * xx - 9.81576),
        linestyle='-.', color='k', linewidth=3.5)

plt.show()
plt.xlim(11.4, 12.6)
plt.ylim(0.45, 2.59)

plt.xlabel('$\log\ (M_{\star}/M_{\odot})$', fontsize=35)
plt.ylabel('$\log\ (R_{\mathrm{90}}/\mathrm{Kpc})$', fontsize=35)

plt.scatter(logmGamaU, logr90GamaU, c='k', marker='o', alpha=0.2)
plt.scatter(logmBcgU, logr90BcgU, c='r', marker='h', s=70, alpha=0.7)

# Red-sequence
xx = np.linspace(10.0, 13.0, 100)

# MCMC result:
# a = 0.960387107648 +0.0755428582291 -0.0750014260391
# b = -10.0134319999 +0.883141848718 -0.889418975333 
# f = 0.168797744636 +0.0104793857322 -0.00963780589871
plt.plot(xx, (0.96038 * xx - 10.01343),
        linestyle='--', color='r', linewidth=3.5)
#yy1 = (aBcgB + aBcgU) * xx + (bBcgB)
#yy2 = (aBcgB - aBcgL) * xx + (bBcgB)
#plt.fill_between(xx, yy1, yy2, 
#                facecolor='r', interpolate=True, alpha=0.15)

# MCMC result:
# a = 0.935475557486 +0.026570857285 -0.026374388974
# b = -9.81576469101 +0.305357355175 -0.307325596592 
# f = 0.15056796201 +0.00316298933649 -0.00311714344817
plt.plot(xx, (0.93546 * xx - 9.81576),
        linestyle='-.', color='k', linewidth=3.5)

plt.show()
plt.xlim(11.4, 12.6)
plt.ylim(-0.19, 1.59)

plt.xlabel('$\log\ (M_{\star}/M_{\odot})$', fontsize=35)
plt.ylabel('$\log\ \gamma $', fontsize=40)

plt.scatter(logmGamaU, (logr50GamaU + 0.94 * (11.0 - logmGamaU)), 
            c='k', alpha=0.15, s=20)

plt.scatter(logmBcgU, (logr50BcgU + 0.94 * (11.0 - logmBcgU)), 
            c='r', alpha=0.8, s=70, marker='h')

plt.show()
plt.xlabel('$\log\ \gamma $', fontsize=40)

gGama = (logr50GamaM + 0.94 * (11.0 - logmGamaM))
gGama = gGama[np.isfinite(gGama)]

gBcg = (logr50BcgM + 0.94 * (11.0 - logmBcgM))
gBcg = gBcg[np.isfinite(gBcg)]
         
plt.hist(gGama, 40, normed=True,
         edgecolor='k', alpha=0.95, histtype='step', linewidth=4.0)

plt.hist(gBcg, 10, normed=True, 
         facecolor='r', alpha=0.4, histtype='stepfilled')

plt.show()<jupyter_output><empty_output><jupyter_text>Separate the samples into ones above/below the median profiles * For Alexie (16-01-21)<jupyter_code>allProfs = copy.deepcopy(bcgProfs)
allProfs += memProfs
allProfs += gamaProfs
print("## Have %d profiles in total" % len(allProfs))
mpStack, mpMed, mpAvg, mpStd = organizeSbp(allProfs, col1='muI1', 
                                           col2='LOGM2LI_C', kind='mass')
mpAbove = []
mpBelow = []

for prof in allProfs: 
    rkpc = prof['rKpc']
    mp = prof['muI1'] + prof.meta['LOGM2LI_C']
    mSep1 = mp[(rkpc >= 40.0) & (rkpc <= 90.0)] - mpMed[1][(rkpc >= 40.0) & (rkpc <= 90.0)]
    mSep2 = mp[(rkpc >= 40.0) & (rkpc <= 90.0)] - mpMed[0][(rkpc >= 40.0) & (rkpc <= 90.0)]
    
    mSep3 = mp[(rkpc >= 10.0) & (rkpc <= 100.0)] - mpMed[1][(rkpc >= 10.0) & (rkpc <= 100.0)]
    mSep4 = mp[(rkpc >= 10.0) & (rkpc <= 100.0)] - mpMed[0][(rkpc >= 10.0) & (rkpc <= 100.0)]
    
    if (np.nanmedian(mSep1) >= 0.03) and (np.nanmax(mSep3) <= 0.4) and (np.nanmin(mSep4) >= -0.45): 
        mpAbove.append(prof)
    if (np.nanmedian(mSep2) <= -0.03) and (np.nanmax(mSep3) <= 0.4) and (np.nanmin(mSep4) >= -0.45):
        mpBelow.append(prof)
        
print("## %d profiles above the median profile" % len(mpAbove))
print("## %d profiles below the median profile" % len(mpBelow))
print("# 11.8 < logM < 12.1; Above the median profile")
print("# RA   DEC   Z   DUMMMY")
for pp in mpAbove: 
    # print(pp.meta['PREFIX'], pp.meta['GALID'])
    if pp.meta['PREFIX'] == 'redBCG':
        indexUse = np.where(bcgTab['ID_CLUSTER'] == int(pp.meta['GALID']))[0][0]
        raUse = bcgTab[indexUse]['RA_BCG']
        decUse = bcgTab[indexUse]['DEC_BCG']
        zUse = bcgTab[indexUse]['z_use']
        print("%10.7f  %10.7f  %10.7f  0.001" % (raUse, decUse, zUse))
        del raUse, decUse, zUse
    elif pp.meta['PREFIX'] == 'redMem':
        indexUse = np.where(memTab['ISEDFIT_ID'] == int(pp.meta['GALID']))[0][0]
        raUse = memTab[indexUse]['RA_MEM']
        decUse = memTab[indexUse]['DEC_MEM']
        zUse = memTab[indexUse]['z_use']
        print("%10.7f  %10.7f  %10.7f  0.001" % (raUse, decUse, zUse))
        del raUse, decUse, zUse
    elif pp.meta['PREFIX'] == 'gama':
        indexUse = np.where(gamaTab['ISEDFIT_ID'] == int(pp.meta['GALID']))[0][0]
        raUse = gamaTab[indexUse]['ra_hsc']
        decUse = gamaTab[indexUse]['dec_hsc']
        zUse = gamaTab[indexUse]['z_use']
        print("%10.7f  %10.7f  %10.7f  0.001" % (raUse, decUse, zUse))
        del raUse, decUse, zUse
aStack, aMed, aAvg, aStd = organizeSbp(mpAbove, col1='muI1', 
                                       col2='LOGM2LI_C', kind='mass')
bStack, bMed, bAvg, bStd = organizeSbp(mpBelow, col1='muI1', 
                                       col2='LOGM2LI_C', kind='mass')
# --------------------------------------------------------------------------------------- #
fig = plt.figure(figsize=(12, 12))
fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
ax1 = fig.add_subplot(111)
ax1.minorticks_on()

# 10 Kpc 
ax1.axvline(10.0 ** 0.25, linewidth=4.0, c='k', linestyle='-', zorder=0, alpha=0.2)
# 100 Kpc
ax1.axvline(100.0 ** 0.25, linewidth=4.0, c='k', linestyle='-', zorder=0, alpha=0.2)

# z = 0.2 : 1"=3.3 Kpc
ax1.axvline(3.3 ** 0.25, linewidth=4.0, c='b', linestyle='--', alpha=0.2, zorder=0)
# z = 0.4 : 1"=5.4 Kpc
ax1.axvline(5.4 ** 0.25, linewidth=4.0, c='b', linestyle='-.', alpha=0.2, zorder=0)

for ss in mpStack:
    ax1.plot(RSMA_COMMON, ss, c='k', alpha=0.02, linewidth=0.8)
for aa in aStack:
    ax1.plot(RSMA_COMMON, aa, c='r', alpha=0.2, linewidth=0.8)
for bb in bStack:
    ax1.plot(RSMA_COMMON, bb, c='b', alpha=0.2, linewidth=0.8)   

ax1.fill_between(RSMA_COMMON, mpMed[0], mpMed[1], 
                 facecolor='k', edgecolor='none', alpha=1.0, zorder=1005)

ax1.fill_between(RSMA_COMMON, aMed[0], aMed[1], 
                 facecolor='r', edgecolor='none', alpha=1.0, zorder=1005)
ax1.fill_between(RSMA_COMMON, bMed[0], bMed[1], 
                 facecolor='b', edgecolor='none', alpha=1.0, zorder=1005)

ax1.text(0.40, 0.90, '$11.6 < \log (M_{\star}) < 11.8$', 
         verticalalignment='bottom', horizontalalignment='left',
         fontsize=40.0, transform=ax1.transAxes)

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

ax1.set_xlim(0.5, 4.1)
ax1.set_ylim(4.01, 9.79)

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)

#fig.savefig('hscMassive_mprof_m2a_1.png', dpi=90)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17534366135718962, 'fraction_numerical': 0.07732926030984072, 'mean_word_length': 3.2484241750092697, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 49, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 47, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '117068', 'license': 'Unlicense', 'n_tokens_mistral': 10694, 'n_tokens_neox': 9711, 'n_words': 1715, 'path': 'note/old/getAverageProfile.ipynb', 'repo_name': 'dr-guangtou/hsc_redbcg', 'size': '22915'}