<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 BGS Morphological Properties
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
Step1:  Read the parent CAST catalog.
Read the parent SDSS (CAST) catalog which defines the sample.
Step4:  Read the g-band model fitting results and select a "good" sample.
Read the g-band model fitting results and select a clean sample using the "finalflag" bit (see Section 2.2 of the data_tables.pdf documentation).
Step5:  Identify the subset of galaxies with good 1- and 2-component fits.
Step6:  Generate some plots.
<ASSISTANT_TASK:>
Python Code:
import os
import numpy as np
import matplotlib.pyplot as plt
import fitsio
from astropy.table import Table
from corner import corner
plt.style.use('seaborn-talk')
%matplotlib inline
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')
Explanation: BGS Morphological Properties
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:
Vikram et al. 2010, PyMorph: Automated Galaxy Structural Parameter Estimation using Python
Meert et al. 2013, Simulations of single- and two-component galaxy decompositions for spectroscopically selected galaxies from the SDSS
Meert et al. 2015, A catalogue of 2D photometric decompositions in the SDSS-DR7 spectroscopic main galaxy sample: preferred models and systematics
Meert et al. 2016, A catalogue of 2D photometric decompositions in the SDSS-DR7 spectroscopic main galaxy sample: extension to g and i bands
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.
Imports and paths--
End of explanation
castfile = os.path.join(basicdir, 'UPenn_PhotDec_CAST.fits')
castinfo = fitsio.FITS(castfile)
castinfo[1]
allcast = castinfo[1].read()
Explanation: Read the parent CAST catalog.
Read the parent SDSS (CAST) catalog which defines the sample.
End of explanation
thisband = 'gband'
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
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
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)))
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]
Explanation: Read the g-band model fitting results and select a "good" sample.
Read the g-band model fitting results and select a clean sample using the "finalflag" bit (see Section 2.2 of the data_tables.pdf documentation).
End of explanation
one = select_meert(gfit, onecomp=True)
two = select_meert(gfit, twocomp=True)
Explanation: Identify the subset of galaxies with good 1- and 2-component fits.
End of explanation
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()))
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')
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)
Explanation: Generate some plots.
End of explanation
<END_TASK>