This notebook is for re-computing a subset of the master list that is in the NSA footprint, including K-corrections
In [1]:
# This changes the current directory to the base saga directory - make sure to run this first!
# This is necessary to be able to import the py files and use the right directories,
# while keeping all the notebooks in their own directory.
import os
import sys
import tempfile
if 'saga_base_dir' not in locals():
saga_base_dir = os.path.abspath('..')
if saga_base_dir not in sys.path:
os.chdir(saga_base_dir)
In [2]:
from __future__ import print_function, division
from collections import Counter, OrderedDict
import numpy as np
from astropy import units as u
from astropy.coordinates import *
from astropy import table, cosmology, constants as cnst
from astropy.visualization import hist as ahist
#SAGA imports
import hosts
import targeting
import kcorrect
from masterlist import masterlist as masterlist_module
In [3]:
%matplotlib inline
from matplotlib import style, pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
plt.rcParams['image.cmap'] = 'viridis'
plt.rcParams['image.origin'] = 'lower'
plt.rcParams['axes.prop_cycle'] = style.library['seaborn-deep']['axes.prop_cycle']
plt.rcParams['figure.figsize'] = (14, 8)
plt.rcParams['axes.titlesize'] = plt.rcParams['axes.labelsize'] = 16
plt.rcParams['xtick.labelsize'] = plt.rcParams['ytick.labelsize'] = 14
In [4]:
# bands used here ugrizJHK
# these are Rx = Ax/E(B-V) for R_V~=3.1
# SFD 98 values
sfd_Rx = np.array([5.155,3.793,2.751,2.086,1.479,0.902,0.576,0.367])
# Schlafly & FInkbeiner 2011
sandf_Rx = np.array([4.239, 3.303, 2.285, 1.698, 1.263, 0.709, 0.449,0.302])
#from Yuan+ 13 : https://arxiv.org/abs/1301.1427v1
yuan_FN_Rx = np.array([4.89, 7.24])
In [5]:
A_MAG = 2.5/np.log(10)
def flux_ivar_to_mag(flux, ivar):
mag = -2.5 * np.log10(flux)
sd = ivar**-0.5
mag_err = A_MAG * sd/flux
return mag, mag_err
def mag_err_to_maggies(mag, mag_err):
flux = 10**(mag/-2.5)
flux_err_sd = mag_err * flux/A_MAG
return flux, flux_err_sd**-2
In [6]:
kcorrect.load_templates()
def load_filters_by_name(bandnmfns, band_shift=0.0, verbose=False):
with tempfile.NamedTemporaryFile() as ntf:
ntf.write('KCORRECT_DIR\n')
for bandnm in bandnmfns:
fn = 'data/filters/{}'.format(bandnm)
truefn = os.path.join(os.environ['KCORRECT_DIR'], fn)
if not os.path.isfile(truefn):
raise IOError('Filter file {} does not exist'.format(truefn))
ntf.write(fn)
ntf.write('\n')
ntf.seek(0)
if verbose:
!cat $ntf.name
kcorrect.load_filters(ntf.name, band_shift)
In [7]:
# used way below when actually computing k-corrections
FNugrizJHK_filters = ['galex_{}UV.par'.format(bandnm) for bandnm in 'FN']
FNugrizJHK_filters.extend(['sdss_{}0.par'.format(bandnm) for bandnm in 'ugriz'])
FNugrizJHK_filters.extend(['twomass_{}.par'.format(bandnm) for bandnm in ('J','H','Ks')])
In [8]:
def do_kcorrs(maggies, maggie_ivars, redshifts, filternames, band_shift_to=None):
if filternames is not None:
load_filters_by_name(filternames)
if np.isscalar(redshifts):
maggies = [maggies]
maggie_ivars = [maggie_ivars]
redshifts = [redshifts]
if len(maggies) != len(maggie_ivars) or len(maggies) != len(redshifts):
raise ValueError('maggies, maggie_ivars, and redshifts must all be matching length')
coeffs = []
rms = []
rm0s = []
for m, mi, z in zip(maggies, maggie_ivars, redshifts):
coeffs.append(kcorrect.fit_nonneg(z, m, mi))
rms.append(kcorrect.reconstruct_maggies(coeffs[-1]))
if band_shift_to:
load_filters_by_name(filternames, band_shift=band_shift_to)
for coeff in coeffs:
rm0s.append(kcorrect.reconstruct_maggies(coeff, redshift=0.))
flux_kcorr = np.array(rms)[:, 1:]/np.array(rm0s)[:, 1:]
mag_kcorr = -2.5*np.log10(flux_kcorr)
return mag_kcorr, flux_kcorr
In [9]:
# this is for the masterlist from dropbox
# masterlist = table.Table.read('SAGADropbox/hosts/masterlist.csv', data_start=2)
# units = table.Table.read('SAGADropbox/hosts/masterlist.csv', data_end=2)
# for c in masterlist.colnames:
# masterlist[c].unit = units[c][0].replace('#', '').replace('degrees', 'degree')
# masterlist
In [10]:
# uncomment this to load the "old" masterlist but make it look like the new one
masterlist = table.Table.read('masterlist/masterlist.csv', data_start=2)
units = table.Table.read('SAGADropbox/hosts/masterlist.csv', data_end=2)
for c in masterlist.colnames:
if c in units.colnames:
masterlist[c].unit = units[c][0].replace('#', '').replace('degrees', 'degree')
else:
del masterlist[c]
vhe_col = np.zeros_like(masterlist['vhelio'])
vhe_col.name = 'vhelio_err'
vhe_col.mask[:] = True
masterlist.add_column(vhe_col, masterlist.colnames.index('vhelio')+1)
masterlist
Out[10]:
In [11]:
nsa = table.Table(hosts.get_nsa())
In [12]:
#for NSA we have this round-about but simple way to get E(B-V) so we can get the S&F extinction
# See end of the notebook for proof that it works
nsa['E(B-V)'] = np.mean(np.array(nsa['EXTINCTION'])[:,-5:]/sfd_Rx[:5], axis=1)
In [13]:
twomassxsc = masterlist_module.load_2mass_xsc('masterlist/2mass_xsc_irsa.tab')
twomassxsc
Out[13]:
In [14]:
colstokeep = ['RA', 'Dec', 'NSAID', 'othername', 'vhelio', 'vhelio_err', 'distance']
submaster0 = masterlist[~masterlist['NSAID'].mask][colstokeep]
submaster0
Out[14]:
In [15]:
#integrate NSA photometry
nsatoadd = nsa[('NSAID', 'E(B-V)')]
for i, band in enumerate('FNugriz'):
nsatoadd[band+'_flux'] = maggies = nsa['SERSICFLUX'][:, i]/1e9
nsatoadd[band+'_flux_ivar'] = ivar = nsa['SERSICFLUX_IVAR'][:, i]*1e18
nsatoadd[band+'_flux'].unit = nsatoadd[band+'_flux_ivar'].unit = 'maggie'
nsatoadd['NSA_kcorrect_' + band] = nsa['KCORRECT'][:, i]
mag, mag_err = flux_ivar_to_mag(maggies, ivar)
nsatoadd[band] = mag
nsatoadd[band+'_err'] = mag_err
submaster1 = table.join(submaster0, nsatoadd, 'NSAID')
submaster1
Out[15]:
Requires X-matching NSA to 2MASS XSC
In [16]:
scsubm = SkyCoord(submaster1['RA'], submaster1['Dec'])
sc2mass = SkyCoord(twomassxsc['ra'], twomassxsc['dec'], unit=u.deg)
idx, d2d, _ = scsubm.match_to_catalog_sky(sc2mass)
fig, (ax1, ax2) = plt.subplots(2, 1)
ax1.hist(d2d.arcsec, histtype='step',bins=100, log=True)
ax2.hist(d2d.arcsec, histtype='step',bins=100, range=(0,60), log=True)
fig.tight_layout()
Why are there so many NSA objects not in 2MASS? For now we'll just accept these but keep around the distance to possibly mask later
In [17]:
submaster2 = submaster1.copy()
submaster2['dist_2MASS_NSA'] = d2d.to(u.arcsec)
In [18]:
idx2mass = twomassxsc[idx]
for band in 'jhk':
twomass_col = '{}_m_ext'.format(band)
twomass_err_col = '{}_msig_ext'.format(band)
submaster2[band.upper()] = mag = idx2mass[twomass_col]
submaster2[band.upper()+"_err"] = mag_err = idx2mass[twomass_err_col]
flux, flux_ivar = mag_err_to_maggies(mag, mag_err)
submaster2[band.upper() + '_flux'] = flux
submaster2[band.upper() + '_flux_ivar'] = flux_ivar
In [19]:
# this is the "correction" one for use with extinction and kcorrecting
submaster3 = submaster2.copy()
In [20]:
submaster3['othername'].fill_value = ''
submaster3['othername'] = submaster3['othername'].filled()
In [21]:
for i, bandnm in enumerate('FNugrizJHK'):
if bandnm in 'FN':
submaster3['A_'+bandnm] = Ax = submaster3['E(B-V)']*yuan_FN_Rx[i]
else:
submaster3['A_'+bandnm] = Ax = submaster3['E(B-V)']*sandf_Rx[i-2]
# this is X in $F_intrinsic = F X$ , that's why the - is missing from 2.5
submaster3['A_'+bandnm+'_flux'] = 10**(Ax/2.5)
In [22]:
maggies = np.array([submaster3[bandnm+'_flux'] for bandnm in 'FNugrizJHK']).T
ivars = np.array([submaster3[bandnm+'_flux_ivar'] for bandnm in 'FNugrizJHK']).T
zs = (u.Quantity(submaster3['vhelio'])/cnst.c).decompose()
mag_kcorrs, flux_kcorrs = do_kcorrs(maggies, ivars, zs, FNugrizJHK_filters)
for i, bandnm in enumerate('FNugrizJHK'):
submaster3['kcorr_' + bandnm] = mag_kcorrs[:, i]
In [23]:
for bandnm in 'ugrizJHK':
extinction = submaster3['A_'+bandnm]
kcorr = submaster3['kcorr_'+bandnm]
distmod = Distance(submaster3['distance']).distmod.value
submaster3['M_' + bandnm] = submaster3[bandnm] - extinction - kcorr - distmod
In [24]:
infocols = ['RA', 'Dec', 'NSAID', 'othername',
'vhelio', 'vhelio_err', 'distance', 'dist_2MASS_NSA']
magcols = []
for bandnm in 'ugrizJHK':
magcols.append(bandnm)
magcols.append(bandnm+'_err')
magcols.append('A_'+bandnm)
magcols.append('M_'+bandnm)
# reset the units
for colnm in magcols:
submaster3[colnm].unit = 'mag'
submaster = submaster3[infocols + magcols]
submaster
Out[24]:
In [62]:
fill_values = [('--', '-1')]
submaster.write('masterlist/submaster.ecsv', format='ascii.ecsv', fill_values=fill_values)
submaster3.write('masterlist/submaster_all.ecsv', format='ascii.ecsv', fill_values=fill_values)
In [63]:
!gzip -f masterlist/submaster.ecsv
!gzip -f masterlist/submaster_all.ecsv
In [64]:
!cp masterlist/submaster.ecsv.gz SAGADropbox/hosts/
!cp masterlist/submaster_all.ecsv.gz SAGADropbox/hosts/
In [67]:
table = table.Table.read('masterlist/submaster.ecsv.gz', format='ascii.ecsv')
In [113]:
z = (u.Quantity(submaster['vhelio'])/cnst.c).decompose()
dmz = Distance(cosmology.WMAP9.luminosity_distance(z), allow_negative=True).distmod
dmd = Distance(submaster['distance']).distmod
dd = dmd - dmz
ahist(dd[np.isfinite(dd)], bins='knuth', histtype='step', range=(-.0003, .001))
plt.xlabel('$DM_{dist} - DM_{vhelio}$')
plt.tight_layout()
Looks like totally-ignorable
Is it true that the sfd_Rx values above give the same E(B-V) values for all 5 bands? If yes, it means Blanton got them straight from the SFD paper
In [131]:
np.array(nsa['EXTINCTION'])[:,-5:]/sfd_Rx[:5]
Out[131]:
In [21]:
maggies = np.array([submaster3[bandnm+'_flux'] for bandnm in 'FNugriz']).T
ivars = np.array([submaster3[bandnm+'_flux_ivar'] for bandnm in 'FNugriz']).T
zs = (u.Quantity(submaster3['vhelio'])/cnst.c).decompose()
mag_kcorrs_sdss, flux_kcorrs_sdss = do_kcorrs(maggies, ivars, zs, FNugrizJHK_filters[:7])
In [22]:
maggies = np.array([submaster3[bandnm+'_flux'] for bandnm in 'JHK']).T
ivars = np.array([submaster3[bandnm+'_flux_ivar'] for bandnm in 'JHK']).T
zs = (u.Quantity(submaster3['vhelio'])/cnst.c).decompose()
mag_kcorrs_2mass, flux_kcorrs_2mass = do_kcorrs(maggies, ivars, zs, FNugrizJHK_filters[-3:])
In [28]:
plt.hist((submaster3['NSA_kcorrect_r'] - mag_kcorrs[:, 4])/submaster3['NSA_kcorrect_r'],
bins=100, range=(-3,3), histtype='step')
plt.hist((submaster3['NSA_kcorrect_r'] - mag_kcorrs_sdss[:, 4])/submaster3['NSA_kcorrect_r'],
bins=100, range=(-3,3), histtype='step')
plt.tight_layout()
Looks like yes, so that's a consistent way to get E(B-V)
In [33]:
omaster = table.Table.read('masterlist/masterlist.csv')
nmaster = table.Table.read('SAGADropbox/hosts/masterlist.csv', data_start=2)
In [34]:
omaster
Out[34]:
In [35]:
nmaster
Out[35]:
In [110]:
omaster_nsa = omaster[~omaster['NSAID'].mask]
nmaster_nsa = nmaster[~nmaster['NSAID'].mask]
merged = table.join(omaster_nsa, nmaster_nsa, ['NSAID'], 'left')
sub_merged = merged['NSAID', 'vhelio_1', 'vhelio_2', 'distance_1', 'distance_2']
sub_merged
Out[110]:
In [111]:
sub_merged[sub_merged['NSAID']==129461]
Out[111]:
In [112]:
sub_merged[sub_merged['NSAID']==2]
Out[112]:
Wha? Why are only some of them different??
This section is getting at why the submaster list seems to have more objects with K-band (2MASS) detections than the masterlist. Or more broadly, why there aren't K-band detections for all the targets, given that 2MASS is supposed to be complete for $M_K<-19.6$ and $cz < 3000$ km/s
First we take a look at two sets of objects in the 2MASS and SDSS image browser: those from the sub-master that have a 2MASS match, and those that do not.
In [25]:
def get_2mass_url(row, show_in_browser=True):
twomass_url_templ = 'http://irsa.ipac.caltech.edu/cgi-bin/2MASS/IM/nph-im_pos?type=at&ds=asky&POS={0[RA]:.6f}d+{0[Dec]:.6f}d&subsz=&date=&scan=&coadd=&key=&band=K'
url = twomass_url_templ.format(row)
if show_in_browser:
import webbrowser
webbrowser.open(url)
return url
def get_sdss_url(row, show_in_browser=True):
sdss_url_templ = 'http://skyserver.sdss.org/dr13/en/tools/chart/navi.aspx?ra={0[RA]:.6f}&dec={0[Dec]:.6f}'
url = sdss_url_templ.format(row)
if show_in_browser:
import webbrowser
webbrowser.open(url)
return url
In [38]:
kvmsk = (submaster['M_K']<-19.6)&(submaster['vhelio']<3000)
goodmatches = submaster[kvmsk&(submaster['dist_2MASS_NSA']<10*u.arcsec)]
badmatches = submaster[kvmsk&(submaster['dist_2MASS_NSA']>60*u.arcsec)]
In [36]:
goodsubset = goodmatches[np.random.permutation(len(goodmatches))[:3]]
for row in goodsubset:
print('NSAID', row['NSAID'], 'othername', row['othername'])
print('2MASS:', get_2mass_url(row))
print('SDSS:', get_sdss_url(row))
print('')
As an aside... This last one (NSAID 76834/NGC 2805) is a fascinating-looking galaxy (apparently in the group Holmberg Holmberg 124)! Is it ram-pressure stripping? If so what's with the funky corner-like morphology? But I digress...
In [37]:
badsubset = badmatches[np.random.permutation(len(badmatches))[:3]]
for row in badsubset:
print('NSAID', row['NSAID'], 'othername', row['othername'])
print('2MASS:', get_2mass_url(row))
print('SDSS:', get_sdss_url(row))
print('')
Examination of these two subsets makes it clear that they are not the same subsets of objects: the "no-match" objects are all notably fainter/lower SB than the matched ones.
So now lets check the subset that are bright in r in addition to being bright in K.
In [40]:
# 3.28 and 4.67 are K and r-band solar magnitudes, so here implicitly assuming the average galaxy has sun-like r-K color
solar_rmK = 4.67 - 3.28
rbrightmsk = (submaster['M_r'] - solar_rmK)<-19.6
In [70]:
plt.figure()
_, b, _ = plt.hist(submaster[rbrightmsk&kvmsk]['dist_2MASS_NSA'],
bins=100, range=(0,1400), log=True,
histtype='step', label='with $r$ cut')
plt.hist(submaster[kvmsk]['dist_2MASS_NSA'], bins=b, log=True,
histtype='step', label='without $r$ cut',
color=list(plt.rcParams['axes.prop_cycle'])[2]['color'])
plt.legend()
plt.xlabel('distance of closest NSA/2MASS match [arcsec]')
plt.figure()
_, b, _ = plt.hist(submaster[rbrightmsk&kvmsk]['dist_2MASS_NSA'],
bins=100, range=(0,60), log=True,
histtype='step', label='with $r$ cut')
plt.hist(submaster[kvmsk]['dist_2MASS_NSA'], bins=b, log=True,
histtype='step', label='without $r$ cut',
color=list(plt.rcParams['axes.prop_cycle'])[2]['color'])
plt.legend()
plt.xlabel('distance of closest NSA/2MASS match [arcsec]')
print('Fraction > 30 arcsec', np.sum(submaster[rbrightmsk&kvmsk]['dist_2MASS_NSA']>30*u.arcsec)/np.sum(rbrightmsk&kvmsk))
print('Fraction > 30 arcsec', np.sum(submaster[kvmsk]['dist_2MASS_NSA']>30*u.arcsec)/np.sum(kvmsk))
So what are those remaining 7%?
In [76]:
boundaryobjs = submaster[(submaster['dist_2MASS_NSA']>30*u.arcsec)&rbrightmsk&kvmsk]
boundarysubset = boundaryobjs[np.random.permutation(len(boundaryobjs))[:5]]
for row in boundarysubset:
print('NSAID', row['NSAID'], 'othername', row['othername'])
print('2MASS:', get_2mass_url(row))
print('SDSS:', get_sdss_url(row))
print('')
All either artifacts due to nearby bright stars/shredding, or faint things close to the edge of detectability.
Conclusion: the reason why there are a large number of "missing" NSA/2MASS matches is that if there's not a K-band magnitude, our K-band magnitude cut did not remove them! If we instead cut on r-band ~= to the expected K-band cut, "missing" objects go away. So the alternative aproach we've taken of saying "if it has K and meets the cut, it's in the sample" is OK, because those that don't have a K-band magnitude are almost invariably faint (or a small fraction of artifacts).