Include photometry with Grizli grism fits

grizli has the ability to include constraints from broad-band photometry in the redshift fitting. Using this capability requires some of the tools provided by the eazy-py module (https://github.com/gbrammer/eazy-py).

The example below is intended to be run in the same directory where you've already run the Grizli-Pipeline notebook to process the GOODS-South ERS grism data.

Here we take photometry from the 3D-HST photometric catalog (Skelton et al.) of the GOODS-South field. The 3D-HST products are available at https://archive.stsci.edu/prepds/3d-hst/.


In [1]:
%matplotlib inline

In [2]:
import os
import numpy as np
import matplotlib.pyplot as plt

import astropy.io.fits as pyfits
import drizzlepac

import grizli
from grizli.pipeline import photoz
from grizli import utils, prep, multifit, fitting

utils.set_warnings()
print('\n Grizli version: ', grizli.__version__)

# Requires eazy-py:  https://github.com/gbrammer/eazy-py
import eazy


The following task in the stsci.skypac package can be run with TEAL:
                                    skymatch                                    
The following tasks in the drizzlepac package can be run with TEAL:
    astrodrizzle       imagefindpars           mapreg              photeq       
     pixreplace           pixtopix            pixtosky        refimagefindpars  
     resetbits          runastrodriz          skytopix           tweakback      
      tweakreg           updatenpol

 Grizli version:  0.8.0-8-gd989737

In [3]:
# Run in the directory where you ran the Grizli-Pipeline notebook and 
# extracted spectra of two objects
root = 'j0332m2743'
os.chdir('{0}/Extractions/'.format(root))

In [4]:
# Fetch 3D-HST catalogs
if not os.path.exists('goodss_3dhst.v4.1.cats.tar.gz'):
    os.system('wget https://archive.stsci.edu/missions/hlsp/3d-hst/RELEASE_V4.0/Photometry/GOODS-S/goodss_3dhst.v4.1.cats.tar.gz')
    os.system('tar xzvf goodss_3dhst.v4.1.cats.tar.gz')

In [5]:
# Preparation for eazy-py
eazy.symlink_eazy_inputs(path=os.path.dirname(eazy.__file__)+'/data', 
                         path_is_env=False)


/Users/gbrammer/anaconda3/envs/grizli-dev/lib/python3.5/site-packages/eazy/data/templates -> ./templates
/Users/gbrammer/anaconda3/envs/grizli-dev/lib/python3.5/site-packages/eazy/data/filters/FILTER.RES.latest -> ./FILTER.RES.latest

In [6]:
### Initialize **eazy.photoz** object
field = 'goodss'
version = 'v4.1'

params = {}
params['CATALOG_FILE'] = '{0}_3dhst.{1}.cats/Catalog/{0}_3dhst.{1}.cat'.format(field, version)
params['Z_STEP'] = 0.002
params['Z_MAX'] = 10

params['MAIN_OUTPUT_FILE'] = '{0}_3dhst.{1}.eazypy'.format(field, version)
params['PRIOR_FILTER'] = 205

# Galactic extinction
params['MW_EBV'] = {'aegis':0.0066, 'cosmos':0.0148, 'goodss':0.0069, 
                    'uds':0.0195, 'goodsn':0.0103}[field]

params['TEMPLATES_FILE'] = 'templates/fsps_full/tweak_fsps_QSF_12_v3.param'

translate_file = '{0}_3dhst.{1}.cats/Eazy/{0}_3dhst.{1}.translate'.format(field, version)

ez = eazy.photoz.PhotoZ(param_file=None, translate_file=translate_file, 
                        zeropoint_file=None, params=params, 
                        load_prior=True, load_products=False)


Read default param file: /Users/gbrammer/anaconda3/envs/grizli-dev/lib/python3.5/site-packages/eazy/data/zphot.param.default
Read CATALOG_FILE: goodss_3dhst.v4.1.cats/Catalog/goodss_3dhst.v4.1.cat
f_F160W e_F160W (205): hst/wfc3/IR/f160w.dat
f_U e_U (103): ESO/vimos_u.res
f_F435W e_F435W (  1): hst/ACS_update_sep07/wfc_f435w_t77.dat
f_F606Wcand e_F606Wcand (236): hst/ACS_update_sep07/wfc_f606w_t81.dat
f_F606W e_F606W (  4): hst/ACS_update_sep07/wfc_f606w_t77.dat
f_R e_R (260): ESO/VIMOS/R.dat
f_F775W e_F775W (  5): hst/ACS_update_sep07/wfc_f775w_t77.dat
f_F814Wcand e_F814Wcand (239): hst/ACS_update_sep07/wfc_f814w_t81.dat
f_F850LP e_F850LP (  7): hst/ACS_update_sep07/wfc_f850lp_t77.dat
f_F850LPcand e_F850LPcand (240): hst/ACS_update_sep07/wfc_f850lp_t81.dat
f_F125W e_F125W (203): hst/wfc3/IR/f125w.dat
f_J e_J ( 34): ESO/isaac_j.res
f_F140W e_F140W (204): hst/wfc3/IR/f140w.dat
f_H e_H ( 36): ESO/isaac_h.res
f_Ks e_Ks ( 37): ESO/isaac_ks.res
f_IRAC1 e_IRAC1 ( 18): IRAC/irac_tr1_2004-08-09.dat
f_IRAC2 e_IRAC2 ( 19): IRAC/irac_tr2_2004-08-09.dat
f_IRAC3 e_IRAC3 ( 20): IRAC/irac_tr3_2004-08-09.dat
f_IRAC4 e_IRAC4 ( 21): IRAC/irac_tr4_2004-08-09.dat
f_U38 e_U38 (107): ESO/wfi_BB_U38_ESO841.res
f_B e_B ( 46): musyc/B_cdfs_tot.dat
f_V e_V ( 50): musyc/V_cdfs_tot.dat
f_Rc e_Rc ( 54): musyc/R_cdfs_tot.dat
f_I e_I ( 58): musyc/I_cdfs_tot.dat
f_IA427 e_IA427 (181): Subaru_MB/IA427.dat
f_IA445 e_IA445 (182): Subaru_MB/IA445.dat
f_IA505 e_IA505 (185): Subaru_MB/IA505.dat
f_IA527 e_IA527 (186): Subaru_MB/IA527.dat
f_IA550 e_IA550 (187): Subaru_MB/IA550.dat
f_IA574 e_IA574 (188): Subaru_MB/IA574.dat
f_IA598 e_IA598 (189): Subaru_MB/IA598.dat
f_IA624 e_IA624 (190): Subaru_MB/IA624.dat
f_IA651 e_IA651 (191): Subaru_MB/IA651.dat
f_IA679 e_IA679 (192): Subaru_MB/IA679.dat
f_IA738 e_IA738 (194): Subaru_MB/IA738.dat
f_IA767 e_IA767 (195): Subaru_MB/IA768.dat
f_IA797 e_IA797 (196): Subaru_MB/IA797.dat
f_IA856 e_IA856 (198): Subaru_MB/IA856.dat
f_tenisJ e_tenisJ (220): WIRCam/cfh8101_J.txt
f_tenisK e_tenisK (222): WIRCam/cfh8302_Ks.txt
Read PRIOR_FILE:  templates/prior_F160W_TAO.dat
Process template tweak_fsps_QSF_12_v3_001.dat.
Process template tweak_fsps_QSF_12_v3_002.dat.
Process template tweak_fsps_QSF_12_v3_003.dat.
Process template tweak_fsps_QSF_12_v3_004.dat.
Process template tweak_fsps_QSF_12_v3_005.dat.
Process template tweak_fsps_QSF_12_v3_006.dat.
Process template tweak_fsps_QSF_12_v3_007.dat.
Process template tweak_fsps_QSF_12_v3_008.dat.
Process template tweak_fsps_QSF_12_v3_009.dat.
Process template tweak_fsps_QSF_12_v3_010.dat.
Process template tweak_fsps_QSF_12_v3_011.dat.
Process template tweak_fsps_QSF_12_v3_012.dat.
Process templates: 27.032 s

In [7]:
## Grism fitting arguments created in Grizli-Pipeline
args = np.load('fit_args.npy')[0]

## First-pass redshift templates, similar to the eazy templates but 
## with separate emission lines
t0 = args['t0']

#############
## Make a helper object for generating photometry in a format that grizli 
## understands. 

## Passing the parameters precomputes a function to quickly interpolate
## the templates through the broad-band filters.  It's not required, 
## but makes the fitting much faster.
## 
## `zgrid` defaults to ez.zgrid, be explicit here to show you can 
## change it.  
phot_obj = photoz.EazyPhot(ez, grizli_templates=t0, zgrid=ez.zgrid)


Process template fsps/fsps_QSF_12_v3_nolines_001.dat.
Process template fsps/fsps_QSF_12_v3_nolines_002.dat.
Process template fsps/fsps_QSF_12_v3_nolines_003.dat.
Process template fsps/fsps_QSF_12_v3_nolines_004.dat.
Process template fsps/fsps_QSF_12_v3_nolines_005.dat.
Process template fsps/fsps_QSF_12_v3_nolines_006.dat.
Process template fsps/fsps_QSF_12_v3_nolines_007.dat.
Process template fsps/fsps_QSF_12_v3_nolines_008.dat.
Process template fsps/fsps_QSF_12_v3_nolines_009.dat.
Process template fsps/fsps_QSF_12_v3_nolines_010.dat.
Process template fsps/fsps_QSF_12_v3_nolines_011.dat.
Process template fsps/fsps_QSF_12_v3_nolines_012.dat.
Process template alf_SSP.dat.
Process template line Ha+NII+SII+SIII+He+PaB.
Process template line OIII+Hb.
Process template line OII+Ne.
Process template line Lya+CIV.

In [9]:
### Find IDs of specific objects to extract, same ones from the notebook
import astropy.units as u
tab = utils.GTable()
tab['ra'] = [53.0657456, 53.0624459]
tab['dec'] = [-27.720518, -27.707018]

# Internal grizli catalog
gcat = utils.read_catalog('{0}_phot.fits'.format(root))

idx, dr = gcat.match_to_catalog_sky(tab)
source_ids = gcat['number'][idx]
tab['id'] = source_ids
tab['dr'] = dr.to(u.mas)
tab['dr'].format='.1f'
tab.show_in_notebook()


Out[9]:
GTable length=2
idxradeciddr
degdegmas
053.0657456-27.72051813955.0
153.0624459-27.70701840554.4

In [10]:
## Find indices in the 3D-HST photometric catalog
idx3, dr3 = ez.cat.match_to_catalog_sky(tab)

## Run the photozs just for comparison.  Not needed for the grism fitting 
## but the photozs and SEDs give you a check that the photometry looks 
## reasonable
ez.param['VERBOSITY'] = 1.
ez.fit_parallel(idx=idx3, verbose=False) 

# or could run on the whole catalog by not specifying `idx`

In [11]:
# Show SEDs with best-fit templates and p(z)
for ix in idx3:
    ez.show_fit(ix, id_is_idx=True)



In [12]:
### Spline templates for dummy grism continuum fits
wspline = np.arange(4200, 2.5e4)
Rspline = 50
df_spl = len(utils.log_zgrid(zr=[wspline[0], wspline[-1]], dz=1./Rspline))
tspline = utils.bspline_templates(wspline, df=df_spl+2, log=True, clip=0.0001)

In [13]:
i=1 # red galaxy
id=tab['id'][i]
ix = idx3[i]

In [14]:
## This isn't necessary for general fitting, but 
## load the grism spectrum here for demonstrating the grism/photometry scaling
beams_file = '{0}_{1:05d}.beams.fits'.format(args['group_name'], id)
mb = multifit.MultiBeam(beams_file, MW_EBV=args['MW_EBV'],
                        fcontam=args['fcontam'], sys_err=args['sys_err'],
                        group_name=args['group_name'])


1 ib6o23rsq_flt.fits G141
2 ib6o21qmq_flt.fits G102
3 ib6o23ruq_flt.fits G141
4 ib6o21r6q_flt.fits G102
5 ib6o21qoq_flt.fits G102
6 ib6o23ryq_flt.fits G141
7 ib6o21r8q_flt.fits G102
8 ib6o23s0q_flt.fits G141

Pull out the photometry of the nearest source matched in the photometric catalog.


In [15]:
# Generate the `phot` dictionary
phot, ii, dd = phot_obj.get_phot_dict(mb.ra, mb.dec)
label = "3DHST Catalog ID: {0}, dr={1:.2f}, zphot={2:.3f}"
print(label.format(ez.cat['id'][ii], dd, ez.zbest[ii]))

print('\n`phot` keys:', list(phot.keys()))
for k in phot:
    print('\n'+k+':\n', phot[k])
    
# Initialize photometry for the MultiBeam object
mb.set_photometry(**phot)


3DHST Catalog ID: 43114, dr=0.34 arcsec, zphot=1.894

`phot` keys: ['flam', 'eflam', 'filters', 'tempfilt']

flam:
 [ 5.31032357e-18  6.51258910e-19  1.01788222e-18 -9.90000000e-18
  1.42415672e-18  1.25707625e-18  1.43078086e-18  1.67037538e-18
  2.16158370e-18  2.15115687e-18  5.66383701e-18  5.79696387e-18
  5.68760983e-18  5.22893860e-18  3.92784465e-18  2.16968029e-18
  1.51247285e-18  8.16519652e-19  2.65448613e-19  8.60773755e-19
  9.65480196e-19  1.35084466e-18  1.31090889e-18  1.94346601e-18
  6.05688510e-19  1.04470312e-18  1.15387250e-18  1.30961081e-18
  1.49704177e-18  1.47147796e-18  1.46392622e-18  1.37800262e-18
  1.28888122e-18  1.11946796e-18  1.15300405e-18  1.36246068e-18
  1.42741707e-18  1.97740952e-18  5.73858766e-18  3.56379507e-18]

eflam:
 [ 5.48181070e-20  1.54040457e-20  2.48978179e-20 -9.90000000e-18
  1.80434085e-20  1.47073024e-20  1.79298325e-20  1.85362378e-20
  2.42573980e-20  3.17347882e-20  5.75268546e-20  6.12075626e-20
  5.79308727e-20  5.61900831e-20  4.11541994e-20  2.18890547e-20
  1.51975622e-20  8.72917592e-21  3.09985235e-21  1.04750045e-19
  2.64449692e-20  2.71366298e-20  2.10939635e-20  5.19790422e-20
  1.02733114e-19  7.25676686e-20  5.95611295e-20  3.33359536e-20
  4.05469657e-20  5.65897667e-20  2.48060374e-20  2.45833211e-20
  2.11254832e-20  2.14194205e-20  2.20540256e-20  4.37648654e-20
  4.71052523e-20  5.60571167e-20  6.08986590e-20  3.70920723e-20]

filters:
 [<eazy.filters.FilterDefinition object at 0x1c2fcf7208>, <eazy.filters.FilterDefinition object at 0x1c23e7f1d0>, <eazy.filters.FilterDefinition object at 0x1c23e77c88>, <eazy.filters.FilterDefinition object at 0x1c2fcf78d0>, <eazy.filters.FilterDefinition object at 0x1c23e77e10>, <eazy.filters.FilterDefinition object at 0x1c2fcf7e10>, <eazy.filters.FilterDefinition object at 0x1c23e77e48>, <eazy.filters.FilterDefinition object at 0x1c2fcf7978>, <eazy.filters.FilterDefinition object at 0x1c23e77eb8>, <eazy.filters.FilterDefinition object at 0x1c2fcf79b0>, <eazy.filters.FilterDefinition object at 0x1c2fcf7198>, <eazy.filters.FilterDefinition object at 0x1c23f690f0>, <eazy.filters.FilterDefinition object at 0x1c2fcf71d0>, <eazy.filters.FilterDefinition object at 0x1c23f69198>, <eazy.filters.FilterDefinition object at 0x1c23f691d0>, <eazy.filters.FilterDefinition object at 0x81d9642e8>, <eazy.filters.FilterDefinition object at 0x1c23e77550>, <eazy.filters.FilterDefinition object at 0x1c23f69dd8>, <eazy.filters.FilterDefinition object at 0x1c23f69e10>, <eazy.filters.FilterDefinition object at 0x1c23e7f2b0>, <eazy.filters.FilterDefinition object at 0x1c23f69390>, <eazy.filters.FilterDefinition object at 0x1c23f69470>, <eazy.filters.FilterDefinition object at 0x1c23f69588>, <eazy.filters.FilterDefinition object at 0x1c23f69668>, <eazy.filters.FilterDefinition object at 0x1c23ec0be0>, <eazy.filters.FilterDefinition object at 0x1c23ec0c88>, <eazy.filters.FilterDefinition object at 0x1c23ec0d30>, <eazy.filters.FilterDefinition object at 0x1c23ec0d68>, <eazy.filters.FilterDefinition object at 0x1c23ec0dd8>, <eazy.filters.FilterDefinition object at 0x1c23ec0e10>, <eazy.filters.FilterDefinition object at 0x1c23ec0e48>, <eazy.filters.FilterDefinition object at 0x1c23ec0e80>, <eazy.filters.FilterDefinition object at 0x1c23ec0eb8>, <eazy.filters.FilterDefinition object at 0x1c23ec0ef0>, <eazy.filters.FilterDefinition object at 0x1c23ec0f60>, <eazy.filters.FilterDefinition object at 0x1c23ec0f98>, <eazy.filters.FilterDefinition object at 0x1c23ec0fd0>, <eazy.filters.FilterDefinition object at 0x1c2fcf7080>, <eazy.filters.FilterDefinition object at 0x1c2fcf7550>, <eazy.filters.FilterDefinition object at 0x1c2fcf7588>]

tempfilt:
 <eazy.photoz.TemplateGrid object at 0x1c2fcf7f28>

In [16]:
# parametric template fit to get reasonable background
sfit = mb.template_at_z(templates=tspline, fit_background=True, 
                        include_photometry=False)
fig = mb.oned_figure(tfit=sfit)

ax = fig.axes[0]
ax.errorbar(mb.photom_pivot/1.e4, mb.photom_flam/1.e-19, 
            mb.photom_eflam/1.e-19, 
            marker='s', color='k', alpha=0.4, linestyle='None',
            label='3D-HST photometry')

ax.legend(loc='upper left', fontsize=8)


Out[16]:
<matplotlib.legend.Legend at 0x1c2a87f518>

Scaling the spectrum to the photometry

For the grism spectra, the total flux of an object is the flux integrated over the segmentation polygon, which corresponds to some isophotal level set by the SExtractor/SEP threshold parameter. While this should often be similar to the total flux definition in the photometric catalog (e.g., FLUX_AUTO), they won't necessarily be the same and there could be small offsets in the normalization between photometry and spectra.

Below we demonstrate a quick way to scale the spectrum to the photometry for more reliable combined fits.


In [17]:
## First example:  no rescaling
z_phot = ez.zbest[ix]

# Reset scale parameter
if hasattr(mb,'pscale'):
    delattr(mb, 'pscale')
    
t1 = args['t1']
tfit = mb.template_at_z(z=z_phot)
print('No rescaling, chi-squared={0:.1f}'.format(tfit['chi2']))
fig = fitting.full_sed_plot(mb, tfit, zfit=None, bin=4)


No rescaling, chi-squared=16260.3

An offset is clearly visible betwen the spectra and photometry. In the next step, we use an internal function to compute a correction factor to bring the spectrum in line with the photometry.

The scale_to_photometry function integrates the observed 1D spectrum (combining all available grisms) through the photometry filter bandpasses and derives a polynomial correction to make the two agree.

NB Because we use the observed spectrum directly, scale_to_photometry requires that at least one observed filter be entirely covered by the grism spectrum. And fitting for a higher order correction requires multiple filters across the grism bandpass (generally at least one filter per polynomial order).


In [18]:
# Reset scale parameter
if hasattr(mb,'pscale'):
    delattr(mb, 'pscale')

# Template rescaling, simple multiplicative factor
scl = mb.scale_to_photometry(order=0)

# has funny units of polynomial coefficients times 10**power, 
# see `grizli.fitting.GroupFitter.compute_scale_array`
# Scale value is the inverse, so, e.g., 
# scl.x = [8.89] means scale the grism spectrum by 10/8.89=1.12
print(scl.x) 

mb.pscale = scl.x 

# Redo template fit
tfit = mb.template_at_z(z=z_phot)
print('Simple scaling, chi-squared={0:.1f}'.format(tfit['chi2']))
fig = fitting.full_sed_plot(mb, tfit, zfit=None, bin=4)


[8.9920776]
Simple scaling, chi-squared=16283.4

The function computed a multiplicative factor 10/8.9=1.1 to apply to the spectrum, which now falls nicely on top of the photometry. Next we fit for a first-order linear scaling (normalization and slope).


In [19]:
# Reset scale parameter
if hasattr(mb,'pscale'):
    delattr(mb, 'pscale')

# Template rescaling, linear fit
scl = mb.scale_to_photometry(order=1)

# has funny units of polynomial coefficients times 10**power, 
# see `grizli.fitting.GroupFitter.compute_scale_array`
# Scale value is the inverse, so, e.g., 
# scl.x = [8.89] means scale the grism spectrum by 10/8.89=1.12
print(scl.x) 

mb.pscale = scl.x 

# Redo template fit
tfit = mb.template_at_z(z=z_phot)
print('Simple scaling, chi-squared={0:.1f}'.format(tfit['chi2']))
fig = fitting.full_sed_plot(mb, tfit, zfit=None, bin=4)


[ 9.07478597 -0.67846019]
Simple scaling, chi-squared=16251.3

A bit of an improvement. But be careful with corrections with order>0 if there is limited photometry available in bands that overlap the spectrum.


In [20]:
# Now run the full redshift fit script with the photometry, which will also do the scaling
order=1
fitting.run_all_parallel(id, phot=phot, verbose=False, 
                         scale_photometry=order+1, zr=[1.5, 2.4])


Run 405
scale_to_photometry: [9.06, -0.80]
Out[20]:
(405, 1, 21.442981004714966)

Redshift PDF

Now compare the redshift PDF for the grism+photometry and photometry-only (photo-z) fits.


In [21]:
zfit = pyfits.open('{0}_{1:05d}.full.fits'.format(root, id))
z_grism = zfit['ZFIT_STACK'].header['Z_MAP']
print('Best redshift: {0:.4f}'.format(z_grism))

# Compare PDFs
pztab = utils.GTable.gread(zfit['ZFIT_STACK'])
plt.plot(pztab['zgrid'], pztab['pdf'], label='grism+3D-HST')

plt.plot(ez.zgrid, ez.pz[ix,:], label='photo-z')

plt.semilogy()
plt.xlim(z_grism-0.05, z_grism+0.05); plt.ylim(1.e-2, 1000)
plt.xlabel(r'$z$'); plt.ylabel(r'$p(z)$')
plt.grid()
plt.legend()


Best redshift: 1.8873
Out[21]:
<matplotlib.legend.Legend at 0x1c3686d550>

Grizli internal photometry

The archive query functions shown in the Grizli-Pipeline notebook can optionally query for any HST imaging that overlaps with the grism exposures. The automatic preparation script will process and align these along with the grism+direct data and produce a simple forced-aperture photometric catalog. Currently this catalog just has matched aperture photometry and does not make any effort at PSF matching.

For cases where no external photometric catalog is available, this photometry can be used to help constrain the redshift fits. The catalog and eazy-py products are automatically generated by the preparation scripts.

The grizli.pipeline.photoz scripts can also try to fetch additional photometry from published Vizier catalogs, e.g., GALEX, SDSS, PanSTARRS, WISE, CFHT-LS.


In [22]:
import grizli.pipeline.photoz

# The catalog is automatically generated with a number of aperture sizes.  A total 
# correction is computed in the detection band, usually a weighted sum of all available 
# WFC3/IR filters, with the correction as the ratio between the aperture flux and the 
# flux over the isophotal segment region, the 'flux' column in the SEP catalog.
aper_ix = 1
total_flux_column = 'flux'

# Get external photometry from Vizier
get_external_photometry = True

# Set `object_only=True` to generate the `eazy.photoz.Photoz` object from the 
# internal photometric catalog without actually running the photo-zs on the catalog 
# with few bands.   
int_ez = grizli.pipeline.photoz.eazy_photoz(root, force=False, object_only=True, 
                  apply_background=True, aper_ix=aper_ix, 
                  apply_prior=True, beta_prior=True, 
                  get_external_photometry=get_external_photometry, 
                  external_limits=3, external_sys_err=0.3, external_timeout=300, 
                  sys_err=0.05, z_step=0.01, z_min=0.01, z_max=12, 
                  total_flux=total_flux_column)


Get external photometry from Vizier
Filter gKmag from II/349/ps1 (N=23)
Filter rKmag from II/349/ps1 (N=23)
Filter iKmag from II/349/ps1 (N=23)
Filter zKmag from II/349/ps1 (N=23)
Filter yKmag from II/349/ps1 (N=23)
Filter W1mag from II/328/allwise (N=22)
Filter W2mag from II/328/allwise (N=22)
Filter FUV from II/312/ais (N=6)
Filter NUV from II/312/ais (N=6)
/usr/local/share/python/eazy-py/eazy-photoz/templates -> ./templates
/usr/local/share/python/eazy-py/eazy-photoz/filters/FILTER.RES.latest -> ./FILTER.RES.latest
Read default param file: /Users/gbrammer/anaconda3/envs/grizli-dev/lib/python3.5/site-packages/eazy/data/zphot.param.default
Read CATALOG_FILE: j0332m2743_phot_apcorr.fits
F334 E334 (334): PAN-STARRS/PS1.g
F335 E335 (335): PAN-STARRS/PS1.r
F336 E336 (336): PAN-STARRS/PS1.i
F337 E337 (337): PAN-STARRS/PS1.z
F338 E338 (338): PAN-STARRS/PS1.y
F244 E244 (244): WISE/RSR-W1.txt
F245 E245 (245): WISE/RSR-W2.txt
F120 E120 (120): CAPAK/galex1500.res
F121 E121 (121): CAPAK/galex2500.res
f098m_corr_1 f098m_ecorr_1 (201): hst/wfc3/IR/f098m.dat
f140w_corr_1 f140w_ecorr_1 (204): hst/wfc3/IR/f140w.dat
Read PRIOR_FILE:  templates/prior_F160W_TAO.dat
Process template tweak_fsps_QSF_12_v3_001.dat.
Process template tweak_fsps_QSF_12_v3_002.dat.
Process template tweak_fsps_QSF_12_v3_003.dat.
Process template tweak_fsps_QSF_12_v3_004.dat.
Process template tweak_fsps_QSF_12_v3_005.dat.
Process template tweak_fsps_QSF_12_v3_006.dat.
Process template tweak_fsps_QSF_12_v3_007.dat.
Process template tweak_fsps_QSF_12_v3_008.dat.
Process template tweak_fsps_QSF_12_v3_009.dat.
Process template tweak_fsps_QSF_12_v3_010.dat.
Process template tweak_fsps_QSF_12_v3_011.dat.
Process template tweak_fsps_QSF_12_v3_012.dat.
Process templates: 3.258 s

In [23]:
# Available apertures
for k in int_ez.cat.meta:
    if k.startswith('APER'):
        print('Aperture {0}: R={1:>4.1f} pix = {2:>4.2f}"'.format(k, int_ez.cat.meta[k], 
                                                                 int_ez.cat.meta[k]*0.06))
#
k = 'APER_{0}'.format(aper_ix)
print('\nAperture used {0}: R={1:>4.1f} pix = {2:>4.2f}"'.format(k, int_ez.cat.meta[k], 
                                                                 int_ez.cat.meta[k]*0.06))


Aperture APER_0: R= 6.0 pix = 0.36"
Aperture APER_1: R= 8.3 pix = 0.50"
Aperture APER_2: R=16.3 pix = 0.98"
Aperture APER_3: R=20.0 pix = 1.20"

Aperture used APER_1: R= 8.3 pix = 0.50"

In [24]:
# Integrate the grism templates through the filters on the redshift grid. 
int_phot_obj = photoz.EazyPhot(int_ez, grizli_templates=t0, zgrid=int_ez.zgrid)


Process template fsps/fsps_QSF_12_v3_nolines_001.dat.
Process template fsps/fsps_QSF_12_v3_nolines_002.dat.
Process template fsps/fsps_QSF_12_v3_nolines_003.dat.
Process template fsps/fsps_QSF_12_v3_nolines_004.dat.
Process template fsps/fsps_QSF_12_v3_nolines_005.dat.
Process template fsps/fsps_QSF_12_v3_nolines_006.dat.
Process template fsps/fsps_QSF_12_v3_nolines_007.dat.
Process template fsps/fsps_QSF_12_v3_nolines_008.dat.
Process template fsps/fsps_QSF_12_v3_nolines_009.dat.
Process template fsps/fsps_QSF_12_v3_nolines_010.dat.
Process template fsps/fsps_QSF_12_v3_nolines_011.dat.
Process template fsps/fsps_QSF_12_v3_nolines_012.dat.
Process template alf_SSP.dat.
Process template line Ha+NII+SII+SIII+He+PaB.
Process template line OIII+Hb.
Process template line OII+Ne.
Process template line Lya+CIV.

In [25]:
# Reset scale parameter
if hasattr(mb,'pscale'):
    delattr(mb, 'pscale')

# Show the SED
int_phot, ii, dd = int_phot_obj.get_phot_dict(mb.ra, mb.dec)

mb.unset_photometry()
mb.set_photometry(**int_phot)

tfit = mb.template_at_z(z=z_phot)
print('No rescaling, chi-squared={0:.1f}'.format(tfit['chi2']))
fig = fitting.full_sed_plot(mb, tfit, zfit=None, bin=4)
fig.axes[0].set_ylim(-5,80)


No rescaling, chi-squared=15911.0
Out[25]:
(-5, 80)

In [26]:
# Run the grism fit with the direct image photometry

# Note that here we show that you can just pass the full photometry object 
# and the script will match the nearest photometric entry to the grism object.
order=0
fitting.run_all_parallel(id, phot_obj=int_phot_obj, verbose=False, 
                         scale_photometry=order+1, zr=[1.5, 2.4])


Run 405
scale_to_photometry: [10.59]
Out[26]:
(405, 1, 21.499361038208008)

Compare chi-squared

In this case of this bright spectrum with strong features, the grism dominates the fit even with the extensive 3D-HST catalog


In [27]:
zfit2 = pyfits.open('{0}_{1:05d}.full.fits'.format(root, id))
pztab2 = utils.GTable.gread(zfit2['ZFIT_STACK'])

plt.plot(ez.zgrid, ez.fit_chi2[ix,:] - ez.fit_chi2[ix,:].min(), 
         label='3D-HST photo-z')

plt.plot(pztab['zgrid'],  pztab['chi2'] - pztab['chi2'].min(), 
         label='grism + 3D-HST photometry')

plt.plot(pztab2['zgrid'], pztab2['chi2'] - pztab2['chi2'].min(), 
         label='grism + direct image photometry')

plt.legend()
plt.xlabel('z'); plt.ylabel(r'$\chi^2$')
plt.xlim(1.5, 2.4); plt.ylim(-200, 4000); plt.grid()