The goal of this notebook is to quantify the average morphological properties of the BGS sample. Specifically, the DESI-Data simulations require knowledge of the mean half-light radii of the bulge and disk components of the sample, as well as the average light fraction between the two components.
These measurements require a detailed study by the BGS, Targeting, and Galaxy & Quasar Physics Working Groups. However, as a quick hack we can use the expectation that the BGS galaxy sample will have very similar properties as the SDSS/Main sample. (BGS will target galaxies to r=20, reaching a median redshift of z=0.2, whereas SDSS/Main targeted galaxies to r=17.7 and reached a median redshift of z=0.1. Although there exists a small amount of luminosity evolution and evolution in the size-mass relation between these two epochs, the amount of evolution is significantly smaller than the scatter in galaxy properties at fixed redshift and stellar mass.)
Fortunately, A. Meert and collaborators have carried out a detailed 2D morphological analysis of galaxies in the SDSS/Main sample and publicly released their catalog. Some of the relevant papers include:
Here we focus on the fits to the g-band SDSS imaging, as this will most closely resemble the r-band selection of the BGS sample.
In [1]:
import os
import numpy as np
import matplotlib.pyplot as plt
import fitsio
from astropy.table import Table
from corner import corner
In [2]:
plt.style.use('seaborn-talk')
%matplotlib inline
In [3]:
basicdir = os.path.join(os.getenv('IM_DATA_DIR'), 'upenn-photdec', 'basic-catalog', 'v2')
adddir = os.path.join(os.getenv('IM_DATA_DIR'), 'upenn-photdec', 'additional-catalogs')
In [4]:
castfile = os.path.join(basicdir, 'UPenn_PhotDec_CAST.fits')
castinfo = fitsio.FITS(castfile)
castinfo[1]
Out[4]:
In [5]:
allcast = castinfo[1].read()
In [6]:
thisband = 'gband'
In [7]:
def photdec_select(finalflag, bit):
"""Select subsets of the catalog using the finalflag bitmask.
1 - good bulge-only galaxy
4 - good disk-only galaxy
10 - good two-component fit (logical_or of flags 11, 12, and 13)
20 - bad total magnitude and size
"""
return finalflag & np.power(2, bit) != 0
In [8]:
def select_meert(modelcat, onecomp=False, twocomp=False):
"""Select various (good) subsets of galaxies.
Args:
modelcat: 'UPenn_PhotDec_Models_[g,r,i]band.fits' catalog.
onecomp (bool): galaxies fitted with single-Sersic model.
twocomp (bool): galaxies fitted with Sersic-exponential model.
Notes:
* Flag 10 is a logical_or of 11, 12, 13.
* Flag 1, 4, and 10 are mutually exclusive.
* If Flag 1 or 4 are set then n_disk,r_disk are -999.
"""
finalflag = modelcat['finalflag']
smalln = modelcat['n_bulge'] < 8
goodr = modelcat['r_bulge'] > 0 # Moustakas hack
two = photdec_select(finalflag, 10)
two = np.logical_and( two, smalln )
two = np.logical_and( two, goodr )
if twocomp:
return two
one = np.logical_or( photdec_select(finalflag, 1), photdec_select(finalflag, 4) )
one = np.logical_and( one, smalln )
if onecomp:
return one
both = np.logical_or( one, two )
return both
In [9]:
measfile = os.path.join(basicdir, 'UPenn_PhotDec_nonParam_{}.fits'.format(thisband))
measinfo = fitsio.FITS(measfile)
fitfile = os.path.join(basicdir, 'UPenn_PhotDec_Models_{}.fits'.format(thisband))
fitinfo = fitsio.FITS(fitfile)
print(measinfo[1], fitinfo[1])
_fit = fitinfo[1].read(columns=['finalflag', 'n_bulge', 'r_bulge'])
good = select_meert(_fit)
goodindx = np.where(good)[0]
nobj = len(goodindx)
print('Selected {}/{} good targets.'.format(nobj, len(_fit)))
In [10]:
fit, meas = [], []
fitfile = os.path.join(basicdir, 'UPenn_PhotDec_Models_{}.fits'.format(thisband))
measfile = os.path.join(basicdir, 'UPenn_PhotDec_NonParam_{}.fits'.format(thisband))
gfit = fitsio.read(fitfile, ext=1, rows=goodindx)
gmeas = fitsio.read(measfile, ext=1, rows=goodindx)
cast = allcast[goodindx]
In [11]:
one = select_meert(gfit, onecomp=True)
two = select_meert(gfit, twocomp=True)
In [12]:
print('g-band range = {:.3f} - {:.3f}'.format(gfit['m_tot'].min(), gfit['m_tot'].max()))
print('Redshift range = {:.4f} - {:.4f}'.format(cast['z'].min(), cast['z'].max()))
In [13]:
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
_ = ax1.hist(cast['z'], bins=100, range=(0, 0.4), alpha=0.5, label='All Galaxies')
_ = ax1.hist(cast['z'][two], bins=100, range=(0, 0.4), alpha=0.5, label='Two-Component Fits')
ax1.legend(loc='upper right')
ax1.set_xlabel('Redshift')
ax1.set_ylabel('Number of Galaxies')
hb = ax2.hexbin(cast['ra'], cast['dec'], C=cast['z'], vmin=0, vmax=0.3,
cmap=plt.cm.get_cmap('RdYlBu'))
cb = plt.colorbar(hb)
cb.set_label('Redshift')
ax2.set_xlabel('RA')
ax2.set_ylabel('Dec')
Out[13]:
In [14]:
labels = [r'$g_{tot}$', r'B/T ($g$-band)', r'Bulge $n$ ($g$-band)',
r'Bulge $r_{50, g}$', r'Disk $r_{50, g}$']
data = np.array([
gfit['m_tot'][two],
gfit['BT'][two],
gfit['n_bulge'][two],
np.log10(gfit['r_bulge'][two]),
np.log10(gfit['r_disk'][two])
]).T
data.shape
_ = corner(data, quantiles=[0.25, 0.50, 0.75], labels=labels,
range=np.repeat(0.9999, len(labels)), verbose=True)
The corner plots above demonstrate -- in so far as this sample is representative of the BGS targets -- that the median (and interquartile) morphological properties are:
In [ ]: