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


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

In [3]:
utils.set_warnings()
print('\n Grizli version: ', grizli.__version__)

# Requires eazy-py:  https://github.com/gbrammer/eazy-py
import eazy
print('\n Eazy-py version: ', eazy.__version__)


 Grizli version:  1.0.dev1427

 Eazy-py version:  0.2.0-56-ge3ee759

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

In [5]:
# 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 [6]:
# Preparation for eazy-py
eazy.symlink_eazy_inputs()


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

In [7]:
### 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/miniconda3/envs/grizli-dev/lib/python3.6/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: 7.974 s

In [9]:
## Grism fitting arguments created in Grizli-Pipeline
args = np.load('fit_args.npy', allow_pickle=True)[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+Hg+Hd.
Process template line OII+Ne.
Process template line Gal-UV-lines.

In [10]:
### 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[10]:
GTable length=2
idxradeciddr
degdegmas
053.0657456-27.72051815239.5
153.0624459-27.70701842456.3

In [11]:
## 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 [12]:
# Show SEDs with best-fit templates and p(z)
for ix in idx3:
    ez.show_fit(ix, id_is_idx=True)



In [13]:
### 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 [14]:
i=1 # red galaxy
id=tab['id'][i]
ix = idx3[i]

In [15]:
## 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'])


load_master_fits: j033216m2743_00424.beams.fits
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 [16]:
# 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: ['source', 'flam', 'eflam', 'filters', 'tempfilt', 'ext_corr', 'pz', 'z_spec']

source:
 unknown

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.46196131e-20  1.49557950e-20  2.42546228e-20 -9.90000000e-18
  1.77289053e-20  1.44819861e-20  1.77216994e-20  1.83360815e-20
  2.40424629e-20  3.14534250e-20  5.72330408e-20  6.08906490e-20
  5.76833625e-20  5.60071594e-20  4.10639701e-20  2.18630684e-20
  1.51836510e-20  8.72917592e-21  3.09985235e-21  1.01662643e-19
  2.58034563e-20  2.65997449e-20  2.07747308e-20  5.14806825e-20
  1.00038512e-19  7.07376897e-20  5.82814313e-20  3.26579855e-20
  3.97720656e-20  5.55767762e-20  2.43866014e-20  2.41886187e-20
  2.08065832e-20  2.11163194e-20  2.17796850e-20  4.32569888e-20
  4.65906733e-20  5.55151418e-20  6.05918805e-20  3.70102209e-20]

filters:
 [<eazy.filters.FilterDefinition object at 0x7f8a0aa4ac88>, <eazy.filters.FilterDefinition object at 0x7f89d89775c0>, <eazy.filters.FilterDefinition object at 0x7f89d886f4e0>, <eazy.filters.FilterDefinition object at 0x7f8a0aa543c8>, <eazy.filters.FilterDefinition object at 0x7f89d886f5c0>, <eazy.filters.FilterDefinition object at 0x7f8a0aa54908>, <eazy.filters.FilterDefinition object at 0x7f89d886f5f8>, <eazy.filters.FilterDefinition object at 0x7f8a0aa54470>, <eazy.filters.FilterDefinition object at 0x7f89d886f668>, <eazy.filters.FilterDefinition object at 0x7f8a0aa544a8>, <eazy.filters.FilterDefinition object at 0x7f8a0aa4ac18>, <eazy.filters.FilterDefinition object at 0x7f89d886fa20>, <eazy.filters.FilterDefinition object at 0x7f8a0aa4ac50>, <eazy.filters.FilterDefinition object at 0x7f89d886fcc0>, <eazy.filters.FilterDefinition object at 0x7f89d886fcf8>, <eazy.filters.FilterDefinition object at 0x7f89d886f8d0>, <eazy.filters.FilterDefinition object at 0x7f89d886f908>, <eazy.filters.FilterDefinition object at 0x7f89d886f940>, <eazy.filters.FilterDefinition object at 0x7f89d886f978>, <eazy.filters.FilterDefinition object at 0x7f89d89776a0>, <eazy.filters.FilterDefinition object at 0x7f89d886fe80>, <eazy.filters.FilterDefinition object at 0x7f89d886ff60>, <eazy.filters.FilterDefinition object at 0x7f89d886f3c8>, <eazy.filters.FilterDefinition object at 0x7f89d88c80f0>, <eazy.filters.FilterDefinition object at 0x7f8a0aa4a6a0>, <eazy.filters.FilterDefinition object at 0x7f8a0aa4a748>, <eazy.filters.FilterDefinition object at 0x7f8a0aa4a7f0>, <eazy.filters.FilterDefinition object at 0x7f8a0aa4a828>, <eazy.filters.FilterDefinition object at 0x7f8a0aa4a898>, <eazy.filters.FilterDefinition object at 0x7f8a0aa4a8d0>, <eazy.filters.FilterDefinition object at 0x7f8a0aa4a908>, <eazy.filters.FilterDefinition object at 0x7f8a0aa4a940>, <eazy.filters.FilterDefinition object at 0x7f8a0aa4a978>, <eazy.filters.FilterDefinition object at 0x7f8a0aa4a9b0>, <eazy.filters.FilterDefinition object at 0x7f8a0aa4aa20>, <eazy.filters.FilterDefinition object at 0x7f8a0aa4aa58>, <eazy.filters.FilterDefinition object at 0x7f8a0aa4aa90>, <eazy.filters.FilterDefinition object at 0x7f8a0aa4ab00>, <eazy.filters.FilterDefinition object at 0x7f8a0aa4af60>, <eazy.filters.FilterDefinition object at 0x7f8a0aa540f0>]

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

ext_corr:
 [0.99637521 0.97091463 0.97416745 0.98257762 0.98257003 0.98468055
 0.98839395 0.98920425 0.99114212 0.99113664 0.99489213 0.99482203
 0.99572518 0.99674005 0.99780379 0.99882013 0.99909616 1.
 1.         0.97054366 0.97574138 0.98021561 0.98486684 0.99041495
 0.97376966 0.97478145 0.9785149  0.97966275 0.98088882 0.98209932
 0.98309139 0.98394454 0.98490508 0.98585025 0.98756234 0.98839753
 0.98907842 0.99033451 0.99496198 0.99778951]

pz:
 None

z_spec:
 -1.0

In [17]:
# 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[17]:
<matplotlib.legend.Legend at 0x7f89e8d41898>

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 [18]:
## 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=13903.6

An small offset is 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 [19]:
# 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)


[9.07218214]
Simple scaling, chi-squared=13943.2

The function computed a multiplicative factor 10/9.5=1.05 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 [20]:
# 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.26162034 -0.88134457]
Simple scaling, chi-squared=13910.9

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 [21]:
# 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 id=424 with fit_args.npy
scale_to_photometry: [8.87, 0.39]
j033216m2743_00424.full.fits
Out[21]:
(424, 1, 22.31031322479248)

Redshift PDF

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


In [22]:
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.8871
Out[22]:
<matplotlib.legend.Legend at 0x7f89ad053b00>

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 [23]:
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 = False #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)


Apply catalog corrections
Compute aperture corrections: i=0, D=0.36" aperture
Compute aperture corrections: i=1, D=0.50" aperture
Compute aperture corrections: i=2, D=0.70" aperture
Compute aperture corrections: i=3, D=1.00" aperture
Compute aperture corrections: i=4, D=1.20" aperture
Compute aperture corrections: i=5, D=1.50" aperture
Compute aperture corrections: i=6, D=3.00" aperture
Apply morphological validity class
Couldn't run morph classification from /Users/gbrammer/miniconda3/envs/grizli-dev/lib/python3.6/site-packages/grizli/data/sep_catalog_junk.pkl
Write j033216m2743_phot_apcorr.fits
Read default param file: /Users/gbrammer/miniconda3/envs/grizli-dev/lib/python3.6/site-packages/eazy/data/zphot.param.default
Read CATALOG_FILE: j033216m2743_phot_apcorr.fits
f098m_tot_1 f098m_etot_1 (201): hst/wfc3/IR/f098m.dat
f140w_tot_1 f140w_etot_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: 0.741 s

In [24]:
# 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 APERMASK: R= 1.0 pix = 0.06"
Aperture APER_0: R= 6.0 pix = 0.36"
Aperture APER_1: R= 8.3 pix = 0.50"
Aperture APER_2: R=11.7 pix = 0.70"
Aperture APER_3: R=16.7 pix = 1.00"
Aperture APER_4: R=20.0 pix = 1.20"
Aperture APER_5: R=25.0 pix = 1.50"
Aperture APER_6: R=50.0 pix = 3.00"

Aperture used APER_1: R= 8.3 pix = 0.50"

In [25]:
# 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+Hg+Hd.
Process template line OII+Ne.
Process template line Gal-UV-lines.

In [26]:
# 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=13635.8
Out[26]:
(-5.0, 80.0)

In [27]:
# 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 id=424 with fit_args.npy
j033216m2743_00424.full.fits
Out[27]:
(424, 1, 16.83737587928772)

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 [28]:
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()