Background
NIRSPEC images taken in low-dispersion spectroscopy mode present
certain headaches in the data reduction process owing to the
skewed position of the spectrum relative to detector columns and
rows. This document describes a procedure for removing these
distortions and extracting spectra from NIRSPEC low-dispersion
images.
Users may find the links on the IRAF Spectroscopy
Docs Page helpful.
Procedure (Short Version)
Having gone through the long procedure above, you may wish to
take a few shortcuts by using the task donspec to do the bookkeeping for you.
It will perform the following steps to your data:
- Generate trace image from image specified in parameter
xref
- Remove x distortions using xdistcor
- Allow optional checking of x alignment
- Remove y distortions using ydistcor
- Allow optional checking of y alignment
- Extract spectra interactively using apall
- Measure median S/N using spec_s2n
- Average spectra together using scombine
- Extract sky band from first image
- Derive wavelength solution using identify
- Apply wavelength solution using dispcor
- Interactively obtain sky flux threshold
- Remove regions affected by bright night sky lines using skyinterp
The script will not do the following:
- Cleaning of cosmic rays and bad pixels
- Flatfielding
- Flux calibration
Because some task parameters are hardwired, this approach is
useful only for data meeting the following criteria:
- Target is a point source
- Calibration can be done from night sky lines
Prior to running this script, you should run the skyplot task to generate your input
coordlist listing relevant night sky emission features
over the wavelength range of interest.
Procedure (Long Version)
- Filter cosmic rays.
If desired, use a routine such as cosmicrays
to remove bad pixels from your set of images.
- Flatten images.
If desired, use flatfield images to flatten your data. But
first, think carefully about whether you want to take this
step. If your flatfield image shows numerous streaks along
the slit, it's probably due to "junk" on the dewar window as
seen through the slit. Unless your data were taken in
stationary tracking mode (unlikely), the pattern of junk may
differ on your flats and science images. Hence, using the
flat to correct your science data could do more harm than
good to your data.
- Modify image headers.
IRAF will eventually need to know which direction is the
dispersion direction on these images, so you'd might as well
tell it now. Use this command to add the required
header parameter to all of your images:
hedit *.fits dispaxis 2 add+ verify- update+
- Create trace image.
Although the curvature of the spectrum in the x
(column) direction is not in itself a problem, it prevents
proper removal of the y (row) distortions in the
image which complicate sky subtraction. In order to map and
remove this distortion, you need an appropriate "trace" image. You
can generate one in two ways:
- Create an image with strong spectra at several place
along the slit, perhaps by co-adding images of your
target (if it's bright) or of standard stars (assuming
the spectrograph was not moved in between your science
images and the calibration star images).
- If you have a flat image taken with the same
instrumental setup as your science images, use the task
mktracer to generate a trace
image from the flat; e.g.
mktracer 24des0057 trace
Here are a sample flat image and trace image generated with the mktracer task.
- Remove x distortion. In this step, you'll shift
all rows of the images to align with the middle row, so that
the spectra will parallel to image columns.
- Remove y distortion.
Next, you need to align the night sky lines with the image
rows to ensure adequate night sky line subtraction. For
this you'll use the task ydistcor.
- Extract Spectra.
With the image now carefully rectified, you can proceed to
extract spectra. Use the
noao.twodspec.apextract.apall task to perform the
extractions. Refer to Section 3.3.2 of the User's
Guide to Reducing Slit Spectra with IRAF.
- Setup.
Set the parameters as shown below:
input = "*ms.fits" List of input images
nfind = 1 Number of apertures to be found automatically
(output = "") List of output spectra
(apertures = "") Apertures
(format = "multispec") Extracted spectra format
(references = "") List of aperture reference images
(profiles = "") List of aperture profile images\n
(interactive = yes) Run task interactively?
(find = yes) Find apertures?
(recenter = no) Recenter apertures?
(resize = no) Resize apertures?
(edit = yes) Edit apertures?
(trace = yes) Trace apertures?
(fittrace = yes) Fit the traced points interactively?
(extract = yes) Extract spectra?
(extras = yes) Extract sky, sigma, etc.?
(review = no) Review extractions?\n
(line = INDEF) Dispersion line
(nsum = -25) Number of dispersion lines to sum or median\n\n
(lower = -10.) Lower aperture limit relative to center
(upper = 10.) Upper aperture limit relative to center
(apidtable = "") Aperture ID table (optional)\n\n# DEFAULT BACKG
(b_function = "chebyshev") Background function
(b_order = 1) Background function order
(b_sample = "-30:-6,6:30") Background sample regions
(b_naverage = -100) Background average or median
(b_niterate = 0) Background rejection iterations
(b_low_reject = 3.) Background lower rejection sigma
(b_high_rejec = 3.) Background upper rejection sigma
(b_grow = 0.) Background rejection growing radius\n\n# APERTU
(width = 5.) Profile centering width
(radius = 10.) Profile centering radius
(threshold = 0.) Detection threshold for profile centering\n\n#
(minsep = 5.) Minimum separation between spectra
(maxsep = 1000.) Maximum separation between spectra
(order = "increasing") Order of apertures\n\n# RECENTERING PARAMETERS\n
(aprecenter = "") Apertures for recentering calculation
(npeaks = INDEF) Select brightest peaks
(shift = yes) Use average shift instead of recentering?\n\n#
(llimit = INDEF) Lower aperture limit relative to center
(ulimit = INDEF) Upper aperture limit relative to center
(ylevel = 0.1) Fraction of peak or intensity for automatic wid
(peak = yes) Is ylevel a fraction of the peak?
(bkg = yes) Subtract background in automatic width?
(r_grow = 0.) Grow limits by this factor
(avglimits = no) Average limits over all apertures?\n\n# TRACING
(t_nsum = 25) Number of dispersion lines to sum
(t_step = 25) Tracing step
(t_nlost = 100) Number of consecutive times profile is lost bef
(t_function = "chebyshev") Trace fitting function
(t_order = 2) Trace fitting function order
(t_sample = "*") Trace sample regions
(t_naverage = 1) Trace average or median
(t_niterate = 10) Trace rejection iterations
(t_low_reject = 3.) Trace lower rejection sigma
(t_high_rejec = 3.) Trace upper rejection sigma
(t_grow = 0.) Trace rejection growing radius\n\n# EXTRACTION
(background = "fit") Background to subtract
(skybox = 1) Box car smoothing length for sky
(weights = "variance") Extraction weights (none|variance)
(pfit = "fit1d") Profile fitting type (fit1d|fit2d)
(clean = yes) Detect and replace bad pixels?
(saturation = INDEF) Saturation level
(readnoise = "30") Read out noise sigma (photons)
(gain = "5") Photon gain (photons/data number)
(lsigma = 4.) Lower rejection threshold
(usigma = 4.) Upper rejection threshold
(nsubaps = 1) Number of subapertures per aperture
- Inspection.
For each spectrum in your image, note whether the target
is easily visible near the middle row (512). For the
ones that aren't, try to determine a column number where
they appear prominent, due either to a lack of night sky
emission lines, or the presence of an emission line in
the object. Note this column number for use in the next
step.
- Define apertures.
Run apall,
answer yes to the query Find
apertures?, and yes to the question
Edit apertures?. The program will plot a cut
across the spectrum and mark the suggested aperture. If
desired, remove it by pressing d and use the
m key to mark the center positions of your target(s).
A bracket should appear above the location of each
object. The width of the bracket indicates the region
that will be extracted to define your spectrum.
Obviously, you want this to obtain as much of the light
from the target as possible without including too much
of the night sky background light. Use the u
and l keys to adjust the upper and lower
limits of the aperture to an appropriate size for your
object. If in doubt, be a little generous with the size
- the ``optimal extraction'' algorithm employed by
apall will weight each pixel appropriately
based on the amount of light.
If an object is not visible, use the command :line
N to move to the column where the spectrum was
determined to be the most visible, then use m
to mark that spectrum and adjust the size as needed.
- Define background.
For each aperture you've defined, review the
corresponding background aperture to verify that it is
positioned appropriately for good sky estimation.
Select an aperture by moving the cursor to it and
pressing the period key. Use the b key to
enter the background editor. If the background must be
changed, use the d key to delete the existing
sample regions. Then, define one background region on
either side of your spectrum using the s key
to indicate the starting and ending points with the
cursor. Or, explicitly type in the range of rows to be
used with the command :samp -15:-10 10:15,
for example. If your object is so close to the edge of
the slit that you can't get a good sky region on one
side, then use only the other side to define your sky.
Press f to fit the background data, and
q to quit the background editor. Repeat this
step for each object.
- Trace spectrum.
After defining the object and background apertures for
each object, you have told IRAF where to find the object
at one particular place on the chip. Now use apall
to track the spectrum across all of the columns in the
``trace'' step. Answer ``yes'' to the queries
Trace spectra?, Fit traced spectra
interactively?, and Trace spectrum 1
interactively?. If your object is bright, this
step will be a piece of cake and IRAF will show you a
plot of the trace versus column number, which will be
well fit by a line (i.e., second-order Chebyshev
polynomial in IRAF-land). If IRAF reports that it lost
the trace, try increasing the nsum parameter
(e.g., :nsum 50). The resulting fit should
be essentially constant and have an RMS of under 0.2
pixels. Delete outlying points as needed with the
d key. When a suitable fit is achieved, type
q to quit and proceed to the next spectrum.
- Extract spectrum.
Once your spectra have been traced, then
go ahead and extract them - have a look at your masterwork!
- Estimate S/N. So just how good are your spectra?
This is relatively easy to determine. Measure the signal-to-noise
(S/N) ratio per pixel in one of two ways:
- Run the splot
program to review your spectra stick to band 1). Select
three representative regions of the spectrum, each about 100
pixels wide. Use the m key to measure
S/N within each region. Note that splot
knows nothing bout the gain and readnoise of the detector - it
simply assumes that the bumps and wiggles in the spectrum
represent the noise, and compares the size of the ``noise''
to the level of signal in the continuum.
- Compute S/N manually. By running splot
in multispec mode, you've saved some useful information.
For each aperture, there are 4 bands in the final
spectrum. Band 1 is the extracted spectrum (signal).
Band 3 is the square root of the pixel variance or, in
other words, the noise. So to compute S/N at
each pixel, just divide the signal spectrum by the noise
spectrum:
imarith myspec[*,N,1] / myspec[*,N,4] snr
where N is the aperture number. If you
inspect this spectrum with
splot or implot,
you'll see that the values vary widely due to the
generally small size of the signal and the wildly varying
intensity of the background. Determine the typical
S/N per pixel by using imstat snr
field="midpt" to compute the median. The spec_s2n task will do this for
you!
- Combine spectra. Unless you moved the cross
disperser in the middle of your image set, your images are
likely to be very well aligned in wavelength. In this case,
you can combine them before doing the wavelength solution.
- Derive wavelength solution. Next, you need to map
pixels onto wavelengths. Review Section 3.4.1 of the User's Guide
to Reducing Slit Spectra with IRAF in preparation for this
step. Customarily, one obtains a solution for the wavelengths
based on exposures of a source of emission lines, as provided
in many cases by a Neon/Argon arc lamp. If you wish to do so,
linelists for the neon and argon lamps are included with this
package. However, your science spectra are likely to be
dominated by emission lines from the night sky, thus providing
you with a built-in wavelength calibration. Follow these
steps to calibrate your spectrum using the night sky emission
lines:
- Create an separate sky spectrum by extracting the sky
information band from one of your science
images:
imcopy myspec[*,1,3] sky
where myspec is the name of your image.
- Generate a plot of the sky lines by using the splot task:
splot sky
Print out a paper copy by moving the cursor into the
plot and hitting the = key.
- Generate a simulated sky plot for the same wavelength
range as your spectra with the skyplot task. Here are the
parameters:
inlines = "wmkonspec$lowd_ir_ohlines.dat" File listing sky lines
outlines = "oh_lines" Output line list [optional]
outimage = "oh_image" Name of image to create [optional]
(wavelength = 15904.) Central wavelength [Angstroms]
(dispersion = 2.901) Dispersion [Angstroms/px]
(fwhm = 5.) Width of sky lines [px]
(ncols = 1024) number of pixels in spectrum
(identify = yes) Identify lines?
(nlines = 30) Number of lines to identify
(graphics = "stdgraph") Graphics outimage device
Edit the parameters for the task using the epar
command and change the values of the parameters
wavelength, dispersion, and
fwhm as appropriate for your data. Run the
task and inspect the resulting plot to verify that
covers the same wavelength regime as your sky spectrum.
Change parameters and re-run if necessary until all of
the lines in your sky spectrum are represented in the
simulated spectrum. The task creates a list of the
nlines strongest OH emission features and
stores it in the file outlines. Create a
hardcopy of the plot using the = key.
- Edit the parameters for the onedspec.identify
task. Here are the task parameters:
images = "sky" Images containing features to be identified
crval = Approximate coordinate (at reference pixel)
cdelt = Approximate dispersion
(section = "middle col") Section to apply to two dimensional images
(database = "database") Database in which to record feature data
(coordlist = "oh_lines.dat") User coordinate list
(units = "Angstroms") Coordinate units
(nsum = "10") Number of lines/columns/bands to sum in 2D imag
(match = -3.) Coordinate list matching limit
(maxfeatures = 30) Maximum number of features for automatic identi
(zwidth = 100.) Zoom graph width in user units
(ftype = "emission") Feature type
(fwidth = 5.) Feature width in pixels
(cradius = 10.) Centering radius in pixels
(threshold = 30000.) Feature threshold for centering
(minsep = 5.) Minimum pixel separation
(function = "spline3") Coordinate function
(order = 1) Order of coordinate function
(sample = "*") Coordinate sample regions
(niterate = 1) Rejection iterations
(low_reject = 4.) Lower rejection sigma
(high_reject = 4.) Upper rejection sigma
(grow = 0.) Rejection growing radius
(autowrite = no) Automatically write to database
(graphics = "stdgraph") Graphics output device
(cursor = "") Graphics cursor input
(aidpars = "") Automatic identification algorithm parameters
Set the parameter coordlist to the name of the
file you generated with skyplot. Set
maxfeatures to at least the number of features
listed in that file. Set the threshold to a reasonable
minimum value for bright features based on an inspection
of your science spectrum (this helps avoid
misidentification of lines.
- Run the onedspec.identify
task to locate the sky lines in your spectrum. Move your
cursor onto the plot, position it over a prominent sky
lines, and press the m key to mark the line.
Verify that it found the right feature by checking for a
vertical mark over the intended line. If okay, enter
the wavelength as written on your hardcopy of the
simulated sky spectrum. Identify two more prominent
features in the spectrum, one near the right end and one
near the left. Then press the f key to fit a
solution. Press q to return to the spectrum,
and press l to automatically mark additional
lines. Press f again to repeat the fit and
press l to view the non-linear component of
the fit. Delete obvious outliers using the d
key and re-do the fit by pressing f. You
should end up with a plot like this.
Then press q twice to quit the task, and
answer yes to the prompt:
Write feature data to the database?
- Assign the sky spectrum to be the wavelength
reference for your image by using the hedit task:
hedit set1_avsigclip refspec1 sky add+ update+
- Apply the wavelength solution via the onedspec.dispcor task. The parameters are:
input = "set1_avsigclip" List of input spectra
output = "set1_dispcor" List of output spectra
(linearize = no) Linearize (interpolate) spectra?
(database = "database") Dispersion solution database
(table = "") Wavelength table for apertures
(w1 = INDEF) Starting wavelength
(w2 = INDEF) Ending wavelength
(dw = INDEF) Wavelength interval per pixel
(nw = INDEF) Number of output pixels
(log = no) Logarithmic wavelength scale?
(flux = yes) Conserve flux?
(samedisp = no) Same dispersion in all apertures?
(global = no) Apply global defaults?
(ignoreaps = no) Ignore apertures?
(confirm = no) Confirm dispersion coordinates?
(listonly = no) List the dispersion coordinates only?
(verbose = yes) Print linear dispersion assignments?
(logfile = "") Log file
The parameters listed will apply the solution without
rebinning the spectrum.
- Excise bad regions. When your target is faint
relative to the night sky, your data are likely to show the
effects of imperfect background subtraction, as in this example. If desired, interpolate
over these messy regions to improve the appearance of your
spectrum using the skyinterp
task:
Gregory D. Wirth
Last modified: Thu Feb 12 09:31:09 HST
|