The fitting tools have been updated in fitting.py to unify the way that the MultiBeam or drizzled stacks are fit.
Update (9/12/17): The code has now been merged into the master branch so the branch command below is not necessary if your repo is up-to-date.
This requires the development branch "clean_fit_outputs" (https://github.com/gbrammer/grizli/tree/clean_fit_outputs). To check out that branch, use the following in the grizli repo directory
git checkout -b clean_fit_outputs origin/clean_fit_outputs
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
print(grizli.__version__)
NB: The example below uses the "ERS" exposures pre-processed as described in the WFC3IR_Reduction notebook.
In [3]:
try:
os.chdir('Prep') # May be necessary depending on where the notebook is run
except:
pass
all_grism_files = ['ib6o21qmq_flt.fits', 'ib6o21qoq_flt.fits', 'ib6o21r6q_flt.fits',
'ib6o21r8q_flt.fits', 'ib6o23rsq_flt.fits', 'ib6o23ruq_flt.fits',
'ib6o23ryq_flt.fits','ib6o23s0q_flt.fits']
grp = grizli.multifit.GroupFLT(grism_files=all_grism_files, 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)
Now provide two ways of saving the extracted 2D spectra of a given object: "beams" cutouts from the FLT files and drizzled combinations of them, combined for the available grisms separately and for combinations of the available PAs.
In [4]:
# Fitting templates
# First is set with combined emission line complexes for the redshift fit
# (don't allow infinite freedom) of the line ratios / fluxes
templ0 = grizli.utils.load_templates(fwhm=1200, line_complexes=True, stars=False,
full_line_list=None, continuum_list=None,
fsps_templates=True)
# Second set has individual line templates for fitting the line fluxes
templ1 = grizli.utils.load_templates(fwhm=1200, line_complexes=False, stars=False,
full_line_list=None, continuum_list=None,
fsps_templates=True)
# Show the template names / dictionary keys
fmt = '{0:<36s} {1:<36s}'
print(fmt.format('templ0', 'templ1'))
print(fmt.format('------', '------'))
for i in range(len(templ1)):
if i > len(templ0)-1:
print(fmt.format('', list(templ1.keys())[i]))
else:
print(fmt.format(list(templ0.keys())[i], list(templ1.keys())[i]))
# Parameters for drizzled line maps
pline = {'kernel': 'point', 'pixfrac': 0.2, 'pixscale': 0.1, 'size': 8, 'wcs': None}
In [5]:
# grp `GroupFLT` object created as defined in WFC3IR_Reduction from the WFC3 ERS grism data
target = 'ers-grism'
# Line-emitter
id=40776
# Weak continuum features
#id=41147
# Strong continuum
#id=43114
# Pull out the 2D cutouts
beams = grp.get_beams(id, size=80)
mb = grizli.multifit.MultiBeam(beams, fcontam=0.5, group_name=target, psf=False)
# Save a FITS file with the 2D cutouts (beams) from the individual exposures
mb.write_master_fits()
# Fit polynomial model for initial continuum subtraction
wave = np.linspace(2000,2.5e4,100)
poly_templates = grizli.utils.polynomial_templates(wave, order=7)
pfit = mb.template_at_z(z=0, templates=poly_templates, fit_background=True,
fitter='lstsq', get_uncertainties=2)
# Drizzle grisms / PAs
hdu, fig = mb.drizzle_grisms_and_PAs(fcontam=0.2, flambda=False, kernel='point',
size=32, zfit=pfit)
# Save drizzle figure FITS file
fig.savefig('{0}_{1:05d}.stack.png'.format(target, id))
hdu.writeto('{0}_{1:05d}.stack.fits'.format(target, id), clobber=True)
In [6]:
from imp import reload
# High-level wrapper script for doing everything (redshift fits, line fluxes, drizzled line
# maps). More explanation of the details of individual steps TBD.
#
# Needs to be able to find {target}_{id:05d}.beams.fits and {target}_{id:05d}.stack.fits
# generated above
out = grizli.fitting.run_all(id, t0=templ0, t1=templ1, fwhm=1200,
zr=[0.5, 2.3], dz=[0.004, 0.0005],
fitter='nnls', group_name=target, prior=None, fcontam=0.,
pline=pline, mask_sn_limit=7, fit_beams=True, fit_stacks=False,
root=target+'_', fit_trace_shift=False, verbose=False,
phot=None, scale_photometry=False, show_beams=True)
In [7]:
# Products
mb, st, fit, tfit, line_hdu = out
# Objects that are used for the fitting. They are
# both subclasses of grizli.fitting.GroupFitter,
# which provides the fitting methods.
print('mb: ', mb.__class__) # From "beams.fits"
print('st: ', st.__class__) # Drizzled spectra from "stack.fits"
In [8]:
# Properties of the fit on the redshift grid
# stored in an astropy.table.Table
print('Fit table:', fit.colnames)
for k in fit.meta:
print(k, fit.meta[k])
plt.plot(fit['zgrid'], fit['chi2'])
Out[8]:
In [9]:
# Results of the fit at the best determined redshift
print(tfit.keys())
print('z = {0}'.format(tfit['z']))
print('Continuum template, cont1d: ', tfit['cont1d'].__class__)
plt.plot(tfit['cont1d'].wave/1.e4, tfit['cont1d'].flux, label='continuum')
plt.plot(tfit['line1d'].wave/1.e4, tfit['line1d'].flux, label='total')
plt.xlim(0.8, 1.7); plt.ylim(0,1.e-18); plt.grid()
plt.xlabel(r'$\lambda$, microns'); plt.ylabel(r'$f_\lambda$, erg/s/cm2/A'); plt.legend()
# cfit, coeffs, covar are coefficients of the template fit
# and their covariance
Out[9]:
In [10]:
fit_hdu = pyfits.open('{0}_{1:05d}.full.fits'.format(target, id))
fit_hdu.info()
In [11]:
# Basic properties of the input data
fit_hdu[0].header
Out[11]:
In [12]:
from astropy.table import Table
# same as the fit table above, redshift fit to the stacked spectra
fit_stack = Table(fit_hdu['ZFIT_STACK'].data)
plt.plot(fit_stack['zgrid'], fit_stack['pdf'], label='Stacked')
# zoom in around the initial best-guess with the individual "beam" spectra
fit_beam = Table(fit_hdu['ZFIT_BEAM'].data)
plt.plot(fit_beam['zgrid'], fit_beam['pdf'], label='Zoom, beams')
plt.xlim(1.6, 1.9); plt.semilogy(); plt.grid()
plt.xlabel('z'); plt.ylabel('PDF(z)'); plt.legend()
Out[12]:
In [13]:
# Best-fit templates, like `tfit` above
templ = Table(fit_hdu['TEMPL'].data)
print(templ.colnames)
In [14]:
# Extensions 'DSCI', 'LINE', 'LINEWHT'... contain the drizzled line maps
print('{0} has lines [{1}]'.format(fit_hdu.filename(), fit_hdu[0].header['HASLINES']))
# Helper script for plotting them, not generated automatically
fig = grizli.fitting.show_drizzled_lines(fit_hdu, size_arcsec=1.6, cmap='plasma_r')
fig.savefig('{0}_{1:05d}.line.png'.format(target, id))
In [15]:
target = 'ers-grism'
id = 43823 # point source
beams = grp.get_beams(id, size=80)
EPSF_FILES = os.path.exists(os.path.join(os.getenv('GRIZLI'), 'CONF',
'PSFSTD_WFC3IR_{0}.fits'.format('F160W')))
mb = grizli.multifit.MultiBeam(beams, fcontam=0.5, group_name=target, psf=EPSF_FILES)
# Stellar templates
tstar = grizli.utils.load_templates(fwhm=1200, line_complexes=True, fsps_templates=True, stars=True)
# Just makes the plot now, outputs not saved
fig = mb.xfit_star(tstar=tstar)