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
from grizli import utils, multifit
print(grizli.__version__)
NB: The example below uses the "ERS" exposures pre-processed as described in the WFC3IR_Reduction notebook.
In [3]:
os.getcwd()
root = 'j033216m2743'
os.chdir('{0}/Extractions'.format(root))
grp = multifit.GroupFLT(grism_files=glob.glob('*GrismFLT.fits'),
catalog='{0}-ir.cat.fits'.format(root),
cpu_count=-1, sci_extn=1, pad=256)
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]:
### Find IDs of specific objects to extract (2 galaxies and a star)
import astropy.units as u
tab = utils.GTable()
tab['ra'] = [53.0657456, 53.0624459, 53.077048]
tab['dec'] = [-27.720518, -27.707018, -27.705928]
idx, dr = grp.catalog.match_to_catalog_sky(tab)
source_ids = grp.catalog['NUMBER'][idx]
tab['id'] = source_ids
tab['dr'] = dr.to(u.mas)
tab['dr'].format='.1f'
In [6]:
id = source_ids[0]
# Pull out the 2D cutouts
beams = grp.get_beams(id, size=80)
mb = multifit.MultiBeam(beams, fcontam=0.2, group_name=root, 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 and make a figure
hdu, fig = mb.drizzle_grisms_and_PAs(fcontam=0.2, flambda=False, kernel='point',
size=32, zfit=pfit, diff=True)
# Save drizzle figure FITS file
fig.savefig('{0}_{1:05d}.stack.png'.format(root, id))
hdu.writeto('{0}_{1:05d}.stack.fits'.format(root, id), clobber=True)
In [7]:
# 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 {root}_{id:05d}.beams.fits and {root}_{id:05d}.stack.fits
# generated above
from grizli import fitting
out = fitting.run_all_parallel(id, get_output_data=True,
fit_beams=True, fit_only_beams=True,
verbose=False,
argsfile='fit_args.npy')
In [8]:
# 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" -- DEPRECATED
In [9]:
# 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])
In [10]:
plt.plot(fit['zgrid'], fit['chi2'])
plt.xlabel('z'); plt.ylabel(r'$\chi^2$')
Out[10]:
In [11]:
# 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[11]:
In [12]:
fit_hdu = pyfits.open('{0}_{1:05d}.full.fits'.format(root, id))
fit_hdu.info()
In [13]:
# Basic properties of the input data
fit_hdu[0].header
Out[13]:
In [14]:
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
if 'ZFIT_BEAM' in fit_hdu:
fit_beam = Table(fit_hdu['ZFIT_BEAM'].data)
plt.plot(fit_beam['zgrid'], fit_beam['pdf'], label='Zoom, beams')
plt.xlim(0.6, 1.9); plt.semilogy(); plt.grid()
plt.ylim(1.e-80, 1e4)
plt.xlabel('z'); plt.ylabel('PDF(z)'); plt.legend()
Out[14]:
In [15]:
# Best-fit templates, like `tfit` above
templ = Table(fit_hdu['TEMPL'].data)
print(templ.colnames)
In [16]:
# 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(root, id))
In [17]:
# Red star
id=source_ids[2]
In [18]:
from grizli.pipeline import auto_script
auto_script.extract(field_root=root, ids=[id], MW_EBV=0,
pline=pline, run_fit=False, grp=grp, diff=True)
from IPython.display import Image
Image(filename='{0}_{1:05d}.R30.png'.format(root, id))
Out[18]:
In [19]:
mb = multifit.MultiBeam('{0}_{1:05d}.beams.fits'.format(root, id),
fcontam=0.2, group_name=root, psf=False, min_sens=0.05)
# Limited set of red stellar templates
tstar = grizli.utils.load_templates(fwhm=1200, line_complexes=True,
fsps_templates=True, stars=True)
# Fit spectral types. Just makes the plot now, outputs not saved
fig, result, tfit = mb.xfit_star(tstar=tstar, fit_background=False,
spline_correction=True, spline_args={'Rspline':5},
oned_args={'bin':2})
In [20]:
hdu, fig = mb.drizzle_grisms_and_PAs(fcontam=0.2, flambda=False, kernel='point',
size=32, zfit=tfit, diff=True)
In [21]:
# Fit for small shifts of the trace perpendicular to the dispersion
# All exposures in a given grism have the same offset
tr = mb.fit_trace_shift(verbose=False)
print('Trace offset: ', tr[0])
In [22]:
# Refit. Chi-squared is a bit improved
fig, result, tfit = mb.xfit_star(tstar=tstar, fit_background=False,
spline_correction=True, spline_args={'Rspline':6},
oned_args={'bin':2})
In [23]:
hdu, fig = mb.drizzle_grisms_and_PAs(fcontam=0.2, flambda=False, kernel='point',
size=32, zfit=tfit, diff=True)
In [24]:
USE_EPSF = True
mb = multifit.MultiBeam('{0}_{1:05d}.beams.fits'.format(root, id),
fcontam=0.2, group_name=root, psf=USE_EPSF, min_sens=0.05)
# Trace offset
tr = mb.fit_trace_shift(verbose=False)
# Fit spectral types. Just makes the plot now, outputs not saved
fig, result, tfit = mb.xfit_star(tstar=tstar, fit_background=False,
spline_correction=True, spline_args={'Rspline':5},
oned_args={'bin':2})
In [25]:
hdu, fig = mb.drizzle_grisms_and_PAs(fcontam=0.2, flambda=False, kernel='point',
size=32, zfit=tfit, diff=True)