The development branch provides a rewrite of the fitting tools with the goal of defining a single redshift fitting algorithm whether or not you're fitting the multifit.MultiBeam or stack.StackFitter objects. Fitting the drizzled spectra is generally the fastest, and is much faster when many individual beam extractions are available for a given object. For example, sources in the HUDF area can contain as many as 180 individual beams from FIGS, 3D-HST and archival observations, and fitting on the multifit.MultiBeam object can take minutes or even hours per source. However, the full WCS information is only preserved for the multifit.MultiBeam objects, so these must be used to drizzle the continuum-subtracted, rectified emission line maps.
Also implemented are now more compact data formats for storing all of the outputs of a given source that can be easily distributed without having to provide all of the original FLT files. The multifit.MultiBeam object now knows how to read/write a single file (*beams.fits) that contains all of the information necessary for performing the fitting analysis.
fitting.GroupFitter object can also incorporate broad-band photometry in the redshift fit, though this will be documented at a later time.
In [1]:
%matplotlib inline
In [2]:
import glob
import time
import os
import numpy as np
import matplotlib.pyplot as plt
import astropy.io.fits as pyfits
import drizzlepac
import grizli
import grizli.stack
In [3]:
# Initialize the GroupFLT object we computed with WFC3IR_Reduction. When loaded from save files
# doesn't much matter what `ref_file` is, since that will be read from the GrismFLT files rather
# than regenerated.
grp = grizli.multifit.GroupFLT(grism_files=glob.glob('*.0?.GrismFLT.fits'), direct_files=[],
ref_file='../Catalog/ERS_goodss_3dhst.v4.0.F160W_orig_sci.fits',
seg_file='../Catalog/ERS_GOODS-S_IR.seg.fits',
catalog='../Catalog/ERS_GOODS-S_IR.cat',
cpu_count=8)
In [4]:
# Extract an object
TARGET = 'ers-grism'
fcontam = 0.2
id = 40776 # emission line object from the other notebook
In [5]:
beams = grp.get_beams(id, size=32)
mb = grizli.multifit.MultiBeam(beams, fcontam=fcontam, group_name=TARGET)
# Make drizzled spectra
hdu, fig = mb.drizzle_grisms_and_PAs(fcontam=fcontam, flambda=False, size=32,
scale=1., kernel='point')
fig.savefig('{0}_{1:05d}.stack.png'.format(TARGET, id))
hdu.writeto('{0}_{1:05d}.stack.fits'.format(TARGET, id), clobber=True)
# Save beam extractions
mb.write_master_fits()
In [6]:
# The continuum model used above is just the flat f-lambda spectrum.
# Fit a polynomial continuum to make it a bit cleaner
wave = np.linspace(3000,2.e4,1000)
tpoly = grizli.utils.polynomial_templates(wave, order=3)
pfit = mb.template_at_z(z=0, templates=tpoly, fitter='lstsq', get_uncertainties=False)
# Re-drizzle with the fit outputs (don't have to do the first drizzle in practice)
hdu, fig = mb.drizzle_grisms_and_PAs(fcontam=fcontam, flambda=False, size=32,
scale=1., kernel='point', zfit=pfit)
fig.savefig('{0}_{1:05d}.stack.png'.format(TARGET, id))
hdu.writeto('{0}_{1:05d}.stack.fits'.format(TARGET, id), clobber=True)
In [7]:
### Initialize templates
# Line complexes for redshift fits.
# Need to have artificially high line FWHM so don't get aliasing in the
# redshift fit.
t0 = grizli.utils.load_templates(fwhm=1200, line_complexes=True, fsps_templates=True)
# Continuum + individual line templates. Can use unresolved line width here
t1 = grizli.utils.load_templates(fwhm=120, line_complexes=False, fsps_templates=True)
# Line drizzle parameters
pline = {'kernel': 'point', 'pixfrac': 0.2, 'pixscale': 0.1, 'size': 8, 'wcs': None}
In [8]:
# Have generated `mb` above directly from the `beams` extracted from `grp`,
# but here show how it can be generated from the `beams.fits` file saved above.
mb = grizli.multifit.MultiBeam('{0}_{1:05d}.beams.fits'.format(TARGET, id),
group_name=TARGET, fcontam=fcontam)
# Redshift fit
m0 = time.time()
fit = mb.xfit_redshift(templates=t0, zr=[0.1, 2.4], dz=[0.004, 0.0005],
prior=None, fitter='nnls', verbose=False)
m1 = time.time()
# Best-fit template
tfit = mb.template_at_z(z=fit.meta['z_risk'][0], templates=t1,
fit_background=True, fitter='nnls')
m2 = time.time()
fig = mb.xmake_fit_plot(fit, tfit)
print('Redshift fit: {0:.1f} s, Template fit: {1:.1f} s'.format(m1-m0, m2-m1))
In [9]:
st = grizli.stack.StackFitter(files='{0}_{1:05d}.stack.fits'.format(TARGET, id),
group_name=TARGET, sys_err=0.02, mask_min=0.1,
fit_stacks=False, fcontam=fcontam, pas=None, extensions=None,
min_ivar=0.01, overlap_threshold=3, verbose=True)
# Available extensions. With `fit_stacks=False`, will be each GRISM,PA combination available. Otherwise, will
# be stacks of all PAs for each grism
print('Extensions: ', st.ext)
In [10]:
# Redshift fit, code is identical
m0 = time.time()
fit = st.xfit_redshift(templates=t0, zr=[0.1, 2.4], dz=[0.004, 0.0005],
prior=None, fitter='nnls', verbose=False)
m1 = time.time()
# Best-fit template
tfit = st.template_at_z(z=fit.meta['z_risk'][0], templates=t1,
fit_background=True, fitter='nnls')
m2 = time.time()
fig = st.xmake_fit_plot(fit, tfit)
print('Redshift fit: {0:.1f} s, Template fit: {1:.1f} s'.format(m1-m0, m2-m1))
In [11]:
### Wrapper script to fit the drizzle spectra first (fast) and then the individual spectra near the best redshift
m0 = time.time()
pline['pixscale'] = 0.1
out = grizli.fitting.run_all(id, t0=t0, t1=t1, fwhm=1200, zr=[0.1, 2.4], dz=[0.004, 0.0005],
fitter='nnls', group_name=TARGET, fit_stacks=False, prior=None,
fcontam=0.2, pline=pline, mask_sn_limit=3, fit_beams=True,
root=TARGET, fit_trace_shift=False, phot=None, verbose=False)
mb, st, fit, tfit, line_hdu = out
m1 = time.time()
print('Run time: {0:.1f} s'.format(m1-m0))
In [12]:
# Output files
ls_str = 'ls -lth {0}_{1:05d}* > lsr; cat lsr'.format(TARGET, id)
os.system(ls_str)
!cat lsr
In [13]:
# Redrizzle with the full continuum model fit
hdu, fig = mb.drizzle_grisms_and_PAs(fcontam=fcontam, flambda=False, size=32,
scale=1., kernel='point', zfit=tfit)
fig.savefig('{0}_{1:05d}.stack.png'.format(TARGET, id))
hdu.writeto('{0}_{1:05d}.stack.fits'.format(TARGET, id), clobber=True)
In [14]:
out = pyfits.open('{0}_{1:05d}.full.fits'.format(TARGET, id)) # The same as line_hdu
out.info()
In [15]:
# PrimaryHDU
h = out[0].header
for k in h:
print('{0:10} = {1} / {2}'.format(k, h[k], h.comments[k]))
In [16]:
# Show the drizzled lines and direct image cutout, which are extensions `DSCI`, `LINE`, etc.
from imp import reload
reload(grizli.fitting)
fig = grizli.fitting.show_drizzled_lines(out, size_arcsec=1.6, cmap='plasma_r')
fig.savefig('{0}_{1:05d}.line.png'.format(TARGET, id))
In [17]:
# Information about the redshift fit is stored in the `ZFIT` extensions.
# The `ZFIT_STACK` extension is the full grid fit with the stacked spectra, and
# `ZFIT_BEAM` is the zoom in with the beam cutout spectra.
#
# Potential fit quality parameters might be (some combination of):
# BIC_TEMP < BIC_POLY >> Bayesian information criteria for the template+z vs. simple polynomial
# fit. If this is False, then the template fit doesn't provide a
# significant improvement (e.g., featureless continuum)
#
# CHIMIN/DOF > xx >> Large reduced chi-squared
#
# MIN_RISK > xx >> "Risk" parameter from Tanaka et al. (2017), from 0 to 1.
# This value is large for relatively flat or multi-modal PDF(z).
zfit_st = grizli.utils.GTable.read(out['ZFIT_STACK'])
zfit_mb = grizli.utils.GTable.read(out['ZFIT_BEAM'])
h = out['ZFIT_STACK'].header
for k in out['ZFIT_STACK'].header:
print('{0:10} = {1} / {2}'.format(k, h[k], h.comments[k]))
In [18]:
print(zfit_st.colnames)
plt.plot(zfit_st['zgrid'], np.log(zfit_st['pdf']))
plt.plot(zfit_mb['zgrid'], np.log(zfit_mb['pdf']))
plt.xlabel('z'); plt.ylabel('PDF(z)')
Out[18]:
In [19]:
## Best-fit templates are in `TEMPL`. Full FSPS templates through far-IR
templ = grizli.utils.GTable.read(out['TEMPL'])
templ.info()
In [20]:
import astropy.units as u
fig = plt.figure(figsize=[8,4])
ax1 = fig.add_subplot(121); ax2 = fig.add_subplot(122)
for ax in [ax1, ax2]:
# Continuum only
ax.plot(templ['wave'].to(u.micron), templ['continuum'], color='k', alpha=0.5)
# With lines
ax.plot(templ['wave'].to(u.micron), templ['full'], color='r', alpha=0.5)
ax.set_xlabel(r'$\lambda\ /\ \mu\mathrm{m}$')
ax1.set_ylim(1.e-19, 1.e-17)
ax1.set_xlim(0.7, 1.8) # zoom around spectrum
ax1.semilogy()
ax1.set_ylabel(r'$f_\lambda$; '+templ['full'].unit.__str__())
ax2.set_ylim(1.e-22, 1.e-17)
ax2.set_xlim(0.2, 200)
ax2.loglog()
Out[20]:
In [21]:
# Emission line information stored in `COVAR` extension, the data of which is the
# fit covariance at best redshift that is used to compute the parameter uncertainties
h = out['COVAR'].header
for k in h:
print('{0:10} = {1} / {2}'.format(k, h[k], h.comments[k]))