Example Q-band Spectroscopic Data Reduction

Outline:

Introduction

This page gives an example of using the msreduce script to
reduce a Q-band spectrum. It is based upon the experience of Kevin Volk
(who uses the first person in this page where needed, mixed in with use
of the third person to refer to the reader).

In this instance I have an observation of a bright science target
and an observation of a standard star both at Q-band in the lowres20
spectral mode of Michelle. The observation was made to study the
variations in a strong feature across the extended source. For the
purposes of discussion here only the overall spectrum across the
entire source will be extracted. The Q-band spectrum of this object
is well known, as the object has been observed from space
with the Infrared Space Observatory (ISO) and previous to that by the
Infrared Astronomical Satellite (IRAS). Both these satellites obtained low
resolution spectra of the object in question. The ISO spectrum is of better
quality (at least in terms of providing continuous wavelength coverage)
than anything that we can do from the ground. Those satellite observations
did not provide any spatial information, however.

The raw data file for the science object observation is
N20051229S0208.fits, while the raw data file for the standard star
observation is N20051229S0211.fits. Both these files can be obtained
from the Gemini Science Archive as they are now public. The science
target is HD 56126, more often referred to as IRAS 07134+1005. The standard
star is Procyon, alpha CMi. This is not one of the usual calibration stars
for which the midir package has a flux curve, so this example carries
out the reduction on the pretext that the standard star was alpha CMa,
Sirius, instead. I will then need to scale the output spectrum down
to allow for the difference in brightness between alpha CMi and alpha CMa,
but for the purposes of illustration here that is not important.
In Q-band the difference in effective spectral shape between an F5IV star,
Procyon, and an A1V star, Sirius, is very small. The difference in
brightness at these wavelengths can be found from the ratio of the IRAS
[25] fluxes, for example, or from the ISO spectra of these objects.

Initial Steps

The first step is run nsheaders to set the instrument parameters
for Michelle. This routine is found in the gnirs part of the
Gemini package, but it is available once you invoke the midir
package as well

Before starting the reduction it is useful to look at the raw data and see
(a) what the atmospheric emission lines look like, so you can get a rough
ID of the features, and (b) whether the spectrum needs to be defringed. This
can be done from the raw data files if you wish, as in

The first display command shows the raw data from one nod position, the
first nod A position, for the object plus the sky. The way this looks in
ds9 is as follows --

Raw NOD frame on the source

You can see the strong atmospheric bands as vertical bright lines in the
image (the spectral direction is from left to right with short wavelengths
to the left, and the up/down direction is along the slit), and the object is
much too faint to be seen in the raw image. Using the IRAF imexam
task you can do a line cut with the "l" cursor key to get a cut across the
image--

Line cut of raw source frame

This figure has been made with the y range set from 0 to 45000 ADU to show
the true line to continuum contrast in fairly dry conditions at Mauna Kea.

The second display command displays the difference spectral image from the
same nod, wherein the object spectrum can be seen--

Single NOD difference frame (chop A - chop B)

The raw spectrum is interrupted frequently by the atmospheric bands where
the background is high and the atmospheric transmission is low. A line cut
along the peak of the spectrum looks like this--

Line cut along the raw difference spectrum

where you see the atmospheric lines as absorption lines instead of
emission lines. There is no fringing in this line cut, so defringing of
the spectrum is not needed.

Of note is that the band at about 20 microns, near pixel 135, has split
into two peaks at this resolution, and so it is not usable for the wavelength
calibration.

Starting the Processing

I have set the parameters for nsreduce as follows: (these should work for
you except that the rawpath value need not be as indicated here)

Note that the parameters are set with the standard star name and raw data
file name for the standard star spectrum, the raw data file name for the
object, a path to these two raw data files, and various of the flags are set.
In this case there is no need to do defringing so fl_defringe is
set to "no". Most of the other parameters have their default values. The
Michelle lowres20 Q-band spectra generally do not show any fringes and do
not need defringing, as is the case here.

The call to msreduce can either be done by editing the parameters
to the values listed above or by the command

with the rawpath parameter having to be set to whatever directory is
appropriate.

Wavelength Calibration

When you run msreduce the task proceeds automatically until the wavelength
calibration is to be done. It starts by taking the raw data files for the source and
for the standard star and stacking the frames to produce a sky frame, to be used for
the wavelength calibration, and a difference frame, from which the spectrum is to be
extracted. These steps need no interactive input. You can instruct msreduce
to skip the stacking of the raw data file by setting the fl_process flag to
"no". Ny own preference is to start clean each time in a new directory each time I
do a reduction, in which case this initial step is required. You can also set the
flags to apply flat and bias frames to the data--these are set by the fl_flat
flag. If this is set to "yes" then the appropriate raw file names need to be put in
the flat and bias fields. In principle applying the flat and bias
frames will improve the extraction of the spectra of extended sources. I have not
generally found that this produces superior results myself. If you are reducing
extended source spectra, I recommend experimenting with the flat and bias options to see
whether these produce good results or not.

The reduction proceeds to the wavelength calibration step next. The tasks attempt to
produce an automatic line identification of a cut across the middle of the sky frame.
These automatic line IDs are always poor in my experience (see below) so you have to
reset the line identifications to the start with the "i" key command, and then
manually identify the lines/bands. One marks each peak by moving the cursor near it
and then hitting the "m" key, and entering the wavelength in Angstroms when
prompted. The proper line IDs to use are listed below.

Result of the Automatic line IDs:

Note that if you get something that looks like this for the initial wavelength
solution

where the central wavelength is 0 angstroms, it means that the nsheaders
command has not been used to set the parameters of the spectroscopy. But even with
the parameters set the wavelength solution is way off; the central wavelength is more
or less correct but the dispersion is much too large leading to an absurd wavelength
range from less than 5 microns to more than 35 microns, while the actual values
are roughly from 16 to 24 microns.

I find it helpful to start by identifying two lines/bands, and doing a
linear fit over the wavelength range. This produces a close enough
wavelength solution to allow the rest of the lines to be automatically
IDed when marked thereafter. To select a linear fit use the commands
:order 2 and "f" to refit. I used the lines at roughly pixels
52 and 234 for this purpose.

Initial linear fit from two line IDs:

Proper line IDs:

The wavelength calibration is good to within a few hundred Angstroms, or a
few times 0.01 micron. If you need higher precision than this you would
have to work harder on the wavelength calibration aided by an atmospheric
emission model at the proper resolution.

Once the initial wavelength calibration is carried out the nstransform
routine finds the same sky lines for different positions along the slit.
This step generally can be done using the automatic re-identification, but
if you are taking the spectrum of an extended source and wants the best
possible wavelength solution for all positions along the slit you may need
to do this interactively. Michelle does not have a lot of curvature of the
sky lines over the array so the automatic line re-identification normally
works well; thus when prompted about whether you wish to examine the
line identifications you can usually answer "N" and let the software trace
the lines without any manual interaction.

The same process is then followed for the standard star, with the sky
lines being at the same places as listed above, resulting in a wavelength
calibrated sky spectrum for file N20051229S0211.

Wavelength calibrated sky spectrum (calibration star):

Defining the Apertures for the Spectra

Next you have to choose the apertures for the object and for the standard,
in this case these were chosen to get all the flux. The standard star is
of course a point source, but the science target is somewhat extended on
the slit and I have chosen a wide aperture to get all the flux. The aperture
is close to 50 pixels wide, extending +18 and -30 pixels from the central
position. You can either use the "l" and "u" keys to mark these limits, as
was done here, or use the colon commands :lower -30 and :upper
18 to set these values.

When selecting the aperture I also select the background regions to the right
and left of the aperture. This is done using the "b" key command to go into
the background selection window, then using "t" key command to clear any
automatically selected background regions. I then use the "s" key to select
reasonable areas for the background and type "f" to fit these. The following
three images show the automatic background you start with for HD 56126, the
result of the manual selection of background regions for HD 56126, and the
manually selected background for the calibration star.

Initial background for the science target spectrum

Manually defined background for the science target spectrum

Manually defined background for the calibration star spectrum

Aperture for the science target spectrum

Once the aperture is defined it is then traced. For Michelle (or T-ReCS) the
tracing does not show much curvature across the array in any of the spectral
modes. Thus a linear trace or a low order polynomial trace is enough. You
generally needs to remove a few points from the fit in regions where there is
no signal. This is done using the "x" key command, followed by using the
"f" key command to re-fit the trace. The default order of the fitting in
the GNIRS package is higher than is needed for Michelle, as shown in the
initial tracing of the science target or of the calibration star shown below.
Use the colon commands ":order" and ":function" commands, if needed, to
change the type of fit to be used. My experience is that :order 3 is
usually all that is needed with any of the polynomial fitting functions.

Initial tracing for the science target spectrum

Final tracing for the science target spectrum

For an extended source the tracing may be less accurate than that for a
calibration star, depending on whether there is any spatial variation in
the spectrum over the slit.

Aperture for the calibration star spectrum

Original tracing for the calibration star spectrum

Linear tracing for the calibration star spectrum

Revised tracing for the calibration star spectrum

Now the script produces wavelength calibrated raw spectra for the science
target and for the calibration star, still rather badly chopped up by the
atmospheric transmission.

Wavelength calibrated science target raw spectrum

Wavelength calibrated calibration star raw spectrum

Ones science spectrum may look quite different than what is shown above,
but if your telluric calibration spectrum looks much different than what
is shown just above then something is wrong.

Telluric Correction

The most difficult part of the reduction is the telluric correction, and the
solution presented here may not be the optimum result that is possible. The
atmospheric bands are deep so the corrections that are to be made are large
for certain parts of the atmospheric window, and this amplifies the noise
in the two spectra.

The method used here is to set dscale and offset to 0.0
with colon commands

:dscale 0.0
:offset 0.0

which is done with your cursor in the plot window. Generally you then have
to rescale the upper plot with the "w" and "e" key commands to see the
useful part of the spectrum, as there are usually large spikes caused by
division of essentially zero by zero in the regions where the atmospheric
bands are totally opaque. It usually takes several attempts before the
scaling in the upper window is set so that you can see the spectrum.

You should also use the "t" and "s" commands in the lower window to select
the regions where there is signal and to exclude the regions just below 24
microns and near 25.5 microns where there is no signal due to strong water
bands. This prevents these bad regions from dominating the RMS that is
calculated across the spectrum. Once this is done, you then manually
enter scale and shift values to get the smoothest possible spectrum, perhaps
using the "a" key command to get the automatic "best shift" for a given
scale value. There is no unique way to do this that I know of. My
preference is keep the scaling value near 1.0, but in some instances as
here a rather different value is the only way to produce a good result.

The telluric corrections are done using the msabsflux routine
which in turn calls the telluric task which is in the
noao.onedspec package in IRAF. The help for the telluric
task should be consulted concerning the various key commands and colon
commands available.

The original plot from telluric via mstelluric

The plot from telluric after dscale and offset set to 0

The plot from telluric after the shift is reset

The plot from telluric after the wavelength regions are
selected using the "t" and "s" keys

The plot from telluric for the "best" shift and scaling

At present finding the best shift and scaling is a matter of
trial and error. You should seek to produce the smoothest possible output
spectrum with the atmospheric bands removed. It is not always possible
to completely remove the bands, but as shown here even in Q-band it
is usually possible to find a combination of shift and scale values that
give a good result.

The Final Spectrum

Even at the "best" shift the band centers show spikes, such as at 21.78
microns. Due to the two regions where there is no signal and the shall
shift needed to best remove the atmospheric bands there are also larger
spikes at the edges of the zero signal regions. These can be removed
using the "j" key in the splot routine, followed by writing out the
new spectrum as a fits file with the "i" key command..

I supplied the standard star name in the call to msreduce and
set the fl_bbody flag to "no", so the resulting spectrum is
(ideally) calibrated in Jy. The blackbody option produces a spectrum
that is correct in shape but for which the level is arbitrary.

I now have a the calibrated spectrum for the science target HD 56126. If
this were a reduction for science I would need to scale the values for the
proper brightness of the object, or (better!) have the expected spectral
values corrected for the atmosphere available in IRAF for Procyon and use
those in the reduction. The spectrum shown is too bright by a factor of
33.97/18.45 = 1.84 with about a 7% uncertainty based upon the IRAS photometry
for these two stars. In this case as only part of the object was sampled
in taking the spectrum the absolute calibration is probably not important.

Plot of the re-scaled spectrum

I get a nice spectrum out that in fact matches very well in shape with
the ISO spectrum of the object, and without problems across the strong
bands. This is in part due to the object being bright, but it indicates
what is possible for ground-based Q-band spectroscopy from Mauna Kea under
good conditions. My final cleaned spectrum is shown below. The Y axis is
in Jansky. The output units can be set with the outtype parameter
in the call to msreduce as you wish.

Plot of the cleaned and re-scaled spectrum

A log of all the keystrokes I used to reduce the spectrum can be
found here along with some comments as to
what is being done. The complete IRAF
terminal log is also available for people to look at.