Radial Velocity Measurement and Reduction Technique

Cookbook with explanations
Chris Dolan, September 1999
 

Contents:

• Introduction
• 1. Setup Steps
• 2. Spectra Extraction Steps
• 3. Clean Up Steps
• 4. Radial Velocities
• Answers to Cookbook Questions

Philosophy
I think that the best way to learn to do something is to see how someone else did and and then do it yourself.  In that spirit, this Cookbook is written as a step-by-step HOWTO guide describing the way I reduce WIYN Hydra spectra for radial velocities.  I have tried to make the main text short so it can serve as a reference, with the tedious details tucked away behind links so you can read them if you want to really understand what's going on.  So, I'm hoping to achieve the productivity of a black box with the flexibility of a full explanation.  I hope it works for you.

I find that the best way to reduce data is to have a text editor window open (I prefer emacs) on one side of the screen and IRAF or a Unix shell on the other side.  I type all my commands in the emacs window and then copy-and-paste it into the IRAF or Unix window.  This way, I always have a permanent record of EXACTLY what I did (I always call the file "whatidid").  That way, I can assure myself that I am reducing the runs consistently with each other.  If you ever look at any of my whatidid files, you might notice that sometimes I go to complicated lengths to avoid a task that could simply be done manually in a text editor.  I do it this way so the next time I try do the same task, or redo a data reduction step, I can just copy-and-paste.  I think the extra time I spend to figure out a pasteable solution is compenstaed by the time (and tedium) I save by not having to manually edit.

Another advantage of the detailed whatidid approach is that if I accidentally delete some data or a hard drive flakes out, I can reproduce the work exactly (I have had to do this!).  Plus, if I forget how to do something, I can just grep through the whatidid files to find it (I have symbolic links to nearly all of my whatidid file in /usr/users/dolan/work/whatidid/)

Mandatory file
In the process of reducing Hydra data for a single observation period, a "WhatIdid" file should be created. The file should contain all commands executed during reductions, information about all unique qualities and or problems with the data and/or datareduction of the specific run. The file should be located in ~wocs/runs/"run-date" directory when the reduction are finished.

Background
It helps a lot to know about the WIYN/Hydra/MOS system before you start reducing the data.  It puts many of the more off-the-wall steps described below in context.
First, there is a Hydra Manual from NOAO which has copious technical details.
[todo: add more]
 
 
 

Setting up and starting IRAF
If you have not run IRAF before on our Unix system, make a directory (commonly ~/iraf) to store IRAF system files.  Then, "cd" to this directory and type "mkiraf" replying "xgterm" to the terminal type.  Then edit the "login.cl" file with your favorite text editor.  Add lines that say
    set imtype = "fits"
    set stdimage = imt2048
(probably best to add it near the other "set" commands).

Another tweak to the login.cl file (not essential, but helpful) is to add the following lines before the "keep" at the end of the file:

onedspec
rv
imred
bias
ccdred
hydra
astutil

Finally, you should add the following lines to your login.cl so you can access the IRAF scripts I have created for simplifying some of the data reduction/analysis steps.

    set dolan = /usr/users/dolan/work/iraf/
    set dolancl = dolan$cl/
    set dotiraf = dolan$dotiraf/
    task utmiddle = dolancl$utmiddle.cl
    task addccdproc = dolancl$hydra/addccdproc.cl
    task rvsky = dolancl$rvsky.cl
    task rvskies = dolancl$rvskies.cl
    task listskies = dolancl$listskies.cl
    task rvset = dolancl$rvset.cl
    task myfxcor = dolancl$myfxcor.cl
    task simpcor = dolancl$simpcor.cl
    task simpset = dolancl$simpset.cl
    task mysimp = dolancl$mysimp.cl
    task $sun_position $rvsun = "$foreign"
    task $rvstar $aperfile $jdfile = "$foreign"

To get these "foreign" tasks to work, you have to add my iraf-bin directory to your path.  The easiest way to do this is:
    % emacs .cshrc
and add the following line somewhere near the bottom of the file.
    set path = ( /usr/users/dolan/iraf/bin $path )
If you plan to run iraf from the same shell you are in now, be sure to type:
    % source ~/.cshrc
    % rehash

Once your login.cl file is set up, fire up and xgterm (simply type xgterm) and cd to your new iraf directory.  In that directory, type "cl".  Then you can cd to whatever directory you want to work in (see below).

In IRAF, type the command:
    cl> setinstrument kpnoheaders
and just leave all of the parameters alone when it asks you.
Q: What does setinstrument do?
 

Copy the images off tape/CD onto disk
(This is optional if you're working off CD-ROM with FITS files, since IRAF can read the data right off the disk)
A typical (good weather) run is roughly 1 GB of raw images.
Q: What are these images?  What kinds, what names, etc.?
Q: Where should I put them on disk?
 

Copy the "dotiraf" files into a subdirectory
I recommend a subdirectory instead of the same directory because there are a lot of files and only a subset of them are useful.  I also recommend that you do copy the dotiraf files to this directory (instead of, say, just referring to the copies in Bob M.'s directory or something like that) so you have a coherent set of data all in one place when/if the data is backed up offline.
Q: What's a dotiraf file?
Q: What's in Bob Mathieu's directory?
 

Create a list of the raw files
Type a command like one of the following

• ls *.fits > rawimages
• ls *.imh > rawimages
• ls dec98b*.imh > rawimages
• etc.

Check the DATE-OBS header keyword for problems
Q: Why is there a problem?  What does it matter?

To find the problem:

• hselect @rawimages $I,date-obs,ut,title

See: an example

To fix the problem:
1) Tedious but easy: use HEDIT to fix each one

• cl> hedit <image> "DATE-OBS" "<correct date>" verify- show-
• Example: cl> hedit dec980116.fits "DATE-OBS" "29/12/98" verify- show-

• See: a solution to the above example
• Example:
% ls dec980{116,117,118,119}.fits > fix1
% ls dec980{120,121,122,123,124,125,126,127,128,129}.fits >> fix1
% ls dec980{130,131,132}.fits >> fix1
cl> hedit @fix1 "DATE-OBS" "29/12/98" verify- show-
• (Certainly Unix aficionados can think of easier ways to do this)

Bias correction
Q: Why do we do this?
Q: Why do we do overscan subtraction instead of subtracting a bias image?

Use the COLBIAS task (in the package noao.imred.bias) to subtract off the overscan strip from the rest of the image, and then throw away the overscan region.  It's easiest to just let IRAF overwrite the original images.

First, set up the default parameters:
    cl> colbias.niterate = 3
    cl> colbias.interactive = no

Then, execute the command:
    cl> colbias @rawimages @rawimages bias="[1:32,*]" trim="[34:1032,*]"

You might have to change the numbers in the "trim" parameter.  I use 1032 because I just readout the middle 1000 pixels of the CCD array when I observe at WIYN.  Some folks read out the full 2048x2048 chip, in which case you should use trim="[34:2080,*]"
 

Set CCDPROC parameter (optional)
The DOHYDRA task (annoyingly) insists that every image it works with must have been processed by the CCDPROC task.  Since I don't use CCDPROC, I'm always getting questions from DOHYDRA saying:
    bias0042: CCDPROC keyword not found.
      Either run CCDPROC or add CCDPROC keyword with HEDIT.
    Add CCDPROC keyword and continue?
Instead of having to answer "y" for every image in the run, I just run a short script I wrote:
    cl> addccdproc @rawimages
 

Get Thorium-Argon lamp line list
Simple.  Just type:
    % ln -s ~dolan/work/sep97/thar.dat thar.dat

Q: What is this line list?  What's it for?
Q: Why not use the linelist provided by IRAF?
 

Setup DOHYDRA parameters

cl> dohydra.arctable = "arctable"
cl> dohydra.readnoise = "NOISE_12"
cl> dohydra.gain = "GAIN_12"
cl> dohydra.fibers = 95
cl> dohydra.width = 6
cl> dohydra.minsep = 6
cl> dohydra.maxsep = 11
cl> dohydra.skysubtract = no
cl> params.lower = -3.
cl> params.upper = 3.
cl> params.weights = "variance"
cl> params.f_function = "legendre"
cl> params.f_order = 1
cl> params.coordlist = "thar.dat"
cl> params.match = 1.

Q: What do all of these parameters do?
 

Build ThAr list
cl> ccdlist *.fits ccdtype=comp name+ > comp.list
% emacs comp.list       (move the "Canonical ThAr" to the front)

Q: What is the ThAr list used for?
Q: What is the "Canonical ThAr"?
 

Build the arctable
cl> ccdlist *.fits ccdtype=object name+ > arctable
% emacs arctable

For each file in the arctable, add a second column which contains a comma-separated list of the ThAr image which go with the object spectra (refer to the telescope log sheets to find out which go with which).

You can look at an example arctable on the samples page.
 

Set up dotiraf files (do now or as needed later)
Q: What are dotiraf files?

I usually do most of this at the beginning, since this type of thing is on my mind while I'm building the arctable.  The goal is to figure out which dotiraf files belong to which objects.  Normally, you should hope that the person taking the data at WIYN wrote down which pointing files go with which objects in the log sheets, or (at bare minimum) that the titles of the images resemble the dotiraf file names.

There are two complications to this:

1. Hydra creates more dotiraf files than you need.  Which do you use?
2. The DOHYDRA task uses a contracted form of the dotiraf filename to distinguish between reduced images.  How do you make sure this contracted format is unique?

1. Hydra creates a dotiraf file whenever you (a) run a "setup field" or (b) do a "setup observe" (see the Hydra manual for an explanation of these).  Since you usually issue each of these commands once per Hydra configuration, there are usually two dotiraf files.  If there were complications during the fiber setup or while acquiring the field, there may be more.  Or, if a configuration was observed more than once during the run, there will be at least 4 dotiraf files (whatch out for this!).  Most of the time, the last one is the right one, but it pays to check.  I always use the Unix command "diff" to test the various files to check the dates/times and to see if any fibers were changed.  When in real doubt, you can compare to the image header (e.g. "imhead dec98b00042.fits long+") to see which dotiraf file is the right one.
2. IRAF uses a contracted form of the dotiraf filename as a suffix on the reduced ThArs and Dome Flats.  This contraction is the first 5 characters and the last character of the dotiraf file.  So, "bluecircle.iraf" becomes "bluecf" and "ngc188_f1_blue.hydra.iraf.2" becomes "ngc182".  The problem of uniqueness is that "ngc188_f2_blue.hydra.iraf.2" also becomes "ngc182"!!  To steer clear of this ambiguity, I usually create a 5-6 character symbolic link to the real dotiraf file.  For example:

    ln -s dotiraf/ngc188_f1_blue.hydra.iraf.2 188f1

Top level: The DOHYDRA command
I think the easiest way to run DOHYDRA is to just put all the parameters on the command line.  The DOHYDRA help documentation suggests setting them with eparam and typing ":go".  Whatever...  I think my method is better for my whatidid file (see Philosophy above).

You should issue a command like the following once for each Hydra configuration:

    cl> dohydra <objects> apref=<flat> flat=<flat> through=<flat> arcs1=@<tharlist> apidtab=<dotiraf>

where

• <objects> is the list of files (from column 1 of the arctable) to reduce
• <flat> is the Dome Flat which goes with the objects (i.e. same Hydra configuration and close in time)
• <tharlist> is the list of ThArs created above (with the canonical ThAr first in the list)
• <dotiraf> is the aperture ID table created at the telescope and set up above

Please read a strong recommendation about the order in which you reduce the science images.

Note that DOHYDAR is pretty smart in that if you abort or crash partway through, it will try to figure out where you left off and recommence at that point.  While this is really quite handy, it sometimes screws up when fragments of files are left around.  When in doubt, you can just delete some files and force DOHYDRA to go back a step or two.  Note that all the output files are either ".ms" files (short for multispectrum) or in the "database" subdirectory.  In some cases below where screwups and aborts are common, I have detailed which files to delete, if any.
 

Step 1: Resize and edit apertures - Telling IRAF which aperture is which
When IRAF says (e.g.):

    Set reference apertures for bias0050
    Resize apertures for bias0050?  (yes):

you should reply "y" and hit Enter.
Q: Why do I want to resize apertures?
 

When IRAF says (for example):

    Edit apertures for bias0050?  (yes):

you should reply "y" and hit Enter.  IRAF will pop up a graph of intensity vs. pixel for a single row cut from the middle of the Dome Flat you entered in the "apref" parameter.  The goal is to make sure the apertures are numbered correctly, and they are always wrong!  The key is that in the way we use Hydra, aperture 1 (the simultaneous comparison spectrum) is always dark.  Since IRAF assigns number "1" to the first aperture it finds, it is wrong.  You must change it so the first aperture is number "2".
        Put the mouse cursor on or to the right of the aperture labeled "1"
        Type "o"
        Type "2" and hit Enter
Now all of the aperture numbers should have been incremented by one.  The first should be "2" and the last (leftmost) should be "99".  You can zoom in using the standard IRAF window keys to check if the apertures are numbered correctly.  One simple test is if the Gaps occur at numbers 58 and 75 (blue fibers) or 68 (red fibers).
        Type "q" to exit the aperture editing task and return to DOHYDRA
Q: Why do I want to edit apertures?
 

Step 2: Trace apertures - IRAF figures out where the apertures lay in the 2D image
Q: What is IRAF doing when it is "tracing apertures"?

IRAF will ask (for example):

    Fit traced positions for bias0050 interactively? (yes):

If "yes", go to Step 2a.
If "no" IRAF will trace apertures silently and procede to Step 3.

(Note: if you said "yes" to "Edit apertures?" above, this question will appear in the graphics window)
At least on your first time using DOHYDRA you should say "yes" to this question.  Later you can usually say "no" because IRAF is quite reliable on this step, but it's good to see it once.
 

Step 2a: Trace apertures interactively
If, however, you say "yes", IRAF thinks for a bit and asks, for example:

    Fit curve to aperture 2 of bias00050 interactively? (yes):

IRAF will ask this question for each aperture.  The possible responses are:

• "y" or "yes" - Trace and show a plot of aperture position vs. pixel for this fiber
• "n" or "no" - Trace without showing a plot
• "Y" or "YES" - Trace and show plots for all remaining apertures without asking again
• "N" or "NO" - Trace the remaining apertures without showing plots (as if you said "no" earlier in Step 2)

Step 3: Extract flat field spectra and compute normalization
IRAF says:

    Create response function bias0050bias0050m35f3.ms
    Extract flat field bias0050
    Fit and ratio flat field bias0050

The long name is just
    <flat image><throughput image><dotiraf>.ms

Q: What is extraction?
This step is almost wholly automated.  When the flat field plot comes up, just press "q".
Q: What is happening during flat field fitting?
 

Step 4: Extract throughput spectra, compute correction and apply to flat field
IRAF says:

    Extract throughput image bias0050
    Correct flat field to throughput image

This step is done automatically.
Q: What is happening during throughput correction?
 

Step 5: Extract canonical ThAr and dispersion correct it
Q: What is the "Canonical ThAr"?

IRAF says:

    Extract arc reference image bias0004
    Determine dispersion solution for bias0004

This is usually the most interactive of all the DOHYDRA steps.  it is not too hard once you learn it, but the details and the time-saving shortcuts can be confusing so pay attention!  :-)

The goal here is to identify some ThAr emission lines and tage them with their corresponding wavelengths.  Given this information, IRAF can compute a pixel-to-wavelength mapping function (aka the dispersion solution).  So what you have to do is figure out which ThAr line is which and type in wavelength numbers for the ones you recognize.  But fortunately you don't have to do by hand all 50 lines x 95 fibers x 30-40 exposures!  IRAF minimizes the routine labor, but this adds complications.

In truth, you only have to do one fiber of one ThAr by hand.  From there, IRAF tries to extrapolate the mapping to other fibers in the same image and then to the fibers of other ThArs.
 

Step 5a: The first spectrum
Initially, IRAF displays a plot of the middle fiber of your canonical ThAr.  Your task here is to use a map of the spectral region to figure out which line is which.  Note that arc lamps do not stay strictly the same over time, so the line intensities may vary somewhat (or even substantially) from the map you choose.

The first question to ask yourself is: Have I already reduced this spectrum before.

If no (i.e. this is your first time running DOHYDRA) skip down one paragraph.

If yes, you can save yourself the effort of having to type "m" on the lines again by "add"ing the previous dispersion solution.
Type (in the plot window):
    :add <name of previous Canonical ThAr>
For example:
    :add bias0004bluecf.ms
Then type "f" to see the results of the added lines.  Note that the solution won't be identical to the previous solution since you have probably used a different Dome Flat (annoying subtlety...)
Q: Why does IRAF make me re-reduce the ThAr?
When you have done this, you can skip to Step 5b below.

The first hint is that for WIYN with MOS, the raw spectrum is always reversed with respect to the catalog.  That is, on the screen the red end of the spectrum is on the left (near pixel 1) while the blue is on the right (near pixel 2048), while it is vice versa in all of the catalogs.  (see an example screen capture of this plot)

Look at the pattern of lines on the screen and in the catalog and figure out which lines are which.  When you have figured them out, put the mouse over a ThAr emission line and type "m" (for "mark").  IRAF puts a yellow dash above the line and prompts you for a wavelength of the line.  Type in the catalog wavelength in Angtroms (just type enough digits to make it unambiguous which line you want -- IRAF searches for the nearest catalog match).  Do this for a bunch of lines.

Q: How much is a "bunch"?

At any time, you can type "f" to compute a best fit pixel-to-wavelength mapping for the lines you have marked.  This will plot the fit residuals on the screen (see a screen-capture example).  Type "q" to get back to the spectrum screen.  If things are screwed up (i.e. you mis-ID'd a line or the fit had a bad extrapolation) you can type "i" to start over from scratch.  When you have enough lines, type "l" (for "lookup").  This will use the lines currently marked to compute a mapping and then, using the mapping, try to interpolate the positions of the (about) 50 brightest lines and look them up the in the linelist file.  I usually hit "l" 2-3 times to iterate the mapping computation with more and more lines.

Once most of the lines are marked (yellow lines should indicate roughly 50 emission lines - see example), type "f" to see the fit residuals.  If everything worked smoothly, the RMS (in the upper right corner of the screen) should be between about 0.009 and 0.005.  the actual points should all have residuals in the range of +/-0.02 Angstroms.  If you see points with residuals as large as 0.1 Angstroms or more, they are probably cosmic ray blemishes or line mis-identifications.  A diamond around the point indicates it is not being used for the mapping calculation.  You can put the mouse over the errant data points and press "d" to delete them.  Then hit "q" to return to the spectrum plot.  Hit "l" to re-identify the line and type "f" to go back to the residual screen.  Did the errant point disappear?  If yes, then you successfully fixed a line mis-ID.  If not, then it's probably a blemish and there's really nothing you can do about it.

While you are looking at the residual plot, look at the top where IRAF says something like (ideally) "total=48, sample=48, rejected=0".  These numbers are among the most useful (with the RMS) in determining the quality of the fit.  In the usual WOCS reduction, using my thar.dat linelist file, the total should always be 48 (although I've seen 49 in rare, screwy circumstances...).  The "sample" tells you how many of the 48 IRAF could find back when you hit "l".  This number is typically 47-48.  If it gets lower than this, there's likely a problem with line mis-IDs.  Finally, the number of "rejected" points tells you how many unacceptable outliers there were from the fit.  If this number is larger than 2-3, you should take a look and try to see why.

When you are happy, hit "q" to return to the spectrum plot and "q" again to quit.

If you screw up badly or if IRAF crashes, you can clean up and redo this step by deleting the file called "database/id<tharname><aperture ID file abbrev>.ms" (for example, "database/idbias0004bluecf.ms")
and re-issuing the dohydra command.

All of this takes some time the first time you try it, but with experience it takes about a minute.
 

Step 5b: The other apertures
By now, you have calibrated the first spectrum of the first ThAr of the run.  Woohoo!  Only about 3000 of them left to go!  But fortunately it gets easier and easier.  Your next task is to calibrate the other 94 (or so) apertures in the Canonical ThAr.

As soon as you hit "q" at the end of Step 5a above, IRAF prints a few lines like:

REIDENTIFY: NOAO/IRAF V2.11EXPORT dolan@sprecher.astro.wisc.edu Fri 14:18:47 08-Jan-99
  Reference image = bias0004bluecf.ms, New image = bias0004bluecf.ms, Refit = yes
       Image Data         Found     Fit  Pix Shift  User Shift  Z Shift      RMS
bias0004bluecf.ms - Ap 48 48/48   47/48    -0.0474     0.00605  1.18E-6  0.00586

Fit dispersion function interactively? (no|yes|NO|YES) (yes):

What this is indicating is that you just manually calibrated one aperture (Ap 49) and now IRAF wants to extrapolate this solution to the other apertures in the ThAr file.  It does this in sequence from Ap 48 (shown here) down to Ap 2 and then from Ap 50 up to Ap 99.  What it does is take all the lines found in the adjacent aperture (in this case, Ap 49) and try to find those same lines in the new aperture.  The REIDENTIFY routine searches in a limited range (+/- 1 Angstrom, as indicated in the params.match setting).  In the case above, IRAF indicates that it indeed found 48/48 (i.e. all) of the lines.  Then, IRAF computes the pixel-to-wavelength mapping, reporting the number of lines used in the fit (in this case 47/48 - that is, one of the lines was rejected as an outlier).  The three "shift" columns are mostly useless and are relative to the manually-ID'd aperture.  The RMS column should yield numbers about 0.006-0.009 (with slightly larger values near Ap 2 or Ap 99 where the signal-to-noise is worse due to vignetting).

If you find that the "Fit" number is small (like 43/43 or 44/48) or the RMS is large, you should say "y" to the question, which will put you into an interactive plot just like in Step 5a.  Edit points as you see fit as above, but try not to delete point just because they are outliers (as that is circular) unless you press "l" to bring them back.  When done, press"q" as above.

If the summary looks good (like "Fit" of "46/48" or "48/48") just type "n".

Pay special attention near then ends of the image (i.e. about apertures 2-8 and 92-99) since this is where the two most detrimental effects come into play: (1) the signal-to-noise is worst and (2) the "shift" is largest (due to effects inherent in the spectrograph), so lines might fall outside of the +/-1 Angstrom tolerance region.

If you screw up an aperture or press "n" on a bad one by mistake, there is (unfortunately) no way to go back.  You have to just Ctrl-C out of DOHYDRA and start over.  If you are in a plot window and Ctrl-C doesn't work, try "I" (for "Interrupt"; it works in most interactive IRAF tasks).  To start over, type "flpr" a few times, then delete the id file (something resembling "database/idbias0004bluecf.ms") and re-execute the dohydra command.
 

Step 5c: Finishing up
When you have finished all 95 apertures, IRAF pauses for a moment and asks something like

bias0004bluecf.ms: w1 = 4992.340677651877, w2 = 5257.230040046399, dw= 0.1277806861526885, nw = 2074
  Change wavelength coordinate assignments? (yes|no|NO):

If you decide at this stage that you screwed up the ID process, you have to not only delete the id file as above, but you also have to delete the reduced Canonical ThAr multispec file, since the spectra have been linearized.  To do this, type something like "imdel bias0004bluecf.ms" or "rm bias0004bluecf.ms.fits".

Q: Is there another way to test the quality of the dispersion solution?
 

Step 6: Reduce science spectra
Once the Canonical ThAr spectrum is fully reduced, DOHYDRA starts on the first science object.  First it extracts the image, performing the flat field correction at the same time.  Once it finishes the extraction, it tries to dispersion correct the science spectra.  To do this, IRAF figures out which ThArs are needed (from the arctable file) and reduces them if they are not already reduced.  This additional ThAr reduction is described in Step 6b below.

Here's some sample output from DOHYDRA while it is performing Steps 6a, 6b and 6c.  In this example, there are three science exposures of M35.  The three science images are bracketed by two ThAr images.  All three science images use these two ThArs to compute their dispersion solution, although they put more weight on the solution from the ThAr that was taken closest in time to themselves.

Extract object spectrum bias0089
Assign arc spectra for bias0089
Extract and reidentify arc spectrum bias0088
bias0088m35f6.ms - Ap 2 48/48   46/48      0.125     -0.0159  -3.1E-6   0.0106
Fit dispersion function interactively? (no|yes|NO|YES) (yes): N
Extract and reidentify arc spectrum bias0092
bias0092m35f6.ms - Ap 2 48/48   46/48      0.123     -0.0155  -3.0E-6  0.00941
Fit dispersion function interactively? (no|yes|NO|YES) (NO): N
Dispersion correct bias0089
Extract object spectrum bias0090
Assign arc spectra for bias0090
Dispersion correct bias0090
Extract object spectrum bias0091
Assign arc spectra for bias0091
Dispersion correct bias0091

Step 6a: Extract science spectra
This step is totally automatic.
 

Step 6b: Dispersion correct more ThAr spectra (if necessary)
If necessary, IRAF extracts and dispersion corrects the ThAr spectra specific to the science spectra.  After extraction, IRAF offers you the opportunity to interactively tweak the dispersion solution, just as in Step 5b above.  I almost always choose not to look at any of the apertures.  Unlike Step 5b, IRAF is usually quite good at doing this step on its own.  The reason for the difference is that in Step 5b, the dispersion solution is being extrapolated from one aperture to another in the same file, e.g. from aperture 52 to 53.  Here instead, the solution is being extrapolated from image to image.  So, when reducing e.g. aperture 53 of the new ThAr, IRAF refers to aperture 53 of the Canonical ThAr.  Since the images only drift by a few hundreths of a pixel over a run for a given aperture, this extrapolation is quite robust.

However, that reassurance written, I still check the dispersion solutions after the fact in Section 3 below to make sure IRAF didn't make a mistake.
 

Step 6c: Dispersion correct and linearize science spectra
This step is done automatically.
Q: What is happening during dispersion correction and linearization?
 
 

Check dispersion solutions

• Linearize each of the reduced ThAr spectra (you do not need to do this for the Canonical ThAr.  This happens automatically in Step 5c).  For example:

cl> dispcor bias0006bluecf.ms "" samedisp=yes verbose=no
• Fire up an image viewer (i.e, run ximtool or saotng or something like that)
• Display each multispec image

For example:
cl> display bias0006bluecf.ms 1
• Look at the ends of the spectra (left and right extrema)
• Check: Do the ThAr emission lines all line up from top to bottom?
• If so, everything is fine
• If not, there was an error in the dispersion solution for one or more apertures.  This was probably caused by mis-identification of one or more ThAr lines.  To solve this, you must
• Delete the offending multispec image.
• Delete any reduced multispec science data which used this ThAr (you can check the arctable file to see which ones used it)
• Delete the line ID file (e.g. "database/idbias0006bluecf.ms")
• Re-execute the DOHYDRA command.  It should pick back up by extracting one of the science images (if any) followed by the offending ThAr.

Sky subtraction
Q: Why is sky subtraction important?
Q: Why do we do sky subtraction separately instead of in DOHYDRA?

First, to make life easy, you should make a list of all the object spectra.  Make sure to edit out the Sky Flats.  You don't want to sky subtract your sky spectra!!!
    cl> ccdlist *.ms.fits ccdtype=object name+ > skysubs
    % emacs skysubs

If you know right away that some of the objects don't need to be sky subtracted (e.g. you got dark time, against all odds) then feel free to remove them from the skysubs list right away.

Next, make a directory for backup copies of all your non-sky-subtracted spectra:
    % mkdir presky
    cl> imcopy @skysubs presky
(of course, you can call the directory anything you want.  "presky" is just the name I use by habit)

Q: Why do we need to make backup copies of the spectra?

Then, execute the sky subtraction command:
    cl> skysub @skysubs

This pops up a confusing plot of all the sky fibers laid on top of each other.
Q: What are the sky fibers?

The first thing to do, is replot so you can see the individual sky spectra.  Type:
    :step 100
Then zoom in so you can see details of the spectra (use "w" and either "y" or "e").  If you zoom in too far, you can use "w" and "u" or "d" (for "up" and "down") to pan the display vertically.  Your goal is to delete any spectra which may negatively influence the sky subtraction.  To delete one, put your cursor over it and press "d".  I recommend keeping a record of which spectra you deleted.
Q: Which sky spectra should I delete?
If you accidentally delete the wrong spectrum, press "e" to bring it back.

When you've deleted all of the bad spectra, press "q".  IRAF will then ask you
    Sky rejection option (none|minmax|avsigclip) (avsigclip):
You should answer avsigclip (or just press return).
Q: Why is avsigclip the best choice for sky subtraction?

Then IRAF averages the undeleted spectra, rejecting outliers via avsigclip.  It saves this resulting spectrum in a file with the same name as the object file with "sky" prepended and the ".ms" removed.  For example, if you are working on the file "bias0089.ms.fits" the resulting averaged sky spectrum would be called "skybias0089.fits".

Finally, IRAF subtracts this averaged sky spectrum from all of the apertures of the object file (even from the "Random Sky"s).  You can test the quality of the sky subtraction by running SPLOT on the object file and looking at one of the sky apertures.  If the sky subtraction was good, it should be near zero (maybe with some cosmic rays) and have no Sun-like features.
 

Coadd science spectra
(SCOMBINE)

[todo: write]
 

Correct the reference time for the spectra
    Fix "utmiddle" keyword

[todo: write]
 

Cull out sky spectra
(SCOPY)

[todo: write]
 

Velocity-correct the sky flats
   Run RVSKIES

[todo: write]
 

Delete the raw images
Finally!  You're done!  You now have fully reduced spectra.  Now it's time to free up some disk space.  Delete all of the unextracted images.

    cl> imdel @rawimages
 
 

    Setup FXCOR parameters
    Identify proper template spectrum
    Run FXCOR (or wrapper script) on some or all science spectra
    Optional: Flag bad regions of the spectrum, choose different template, mark as binary, reset background level, etc. and go back to previous step
    Fix output file (e.g. applying any catalog ID changes)
    Incorporate the output into the database appropriate to its cluster
 

What are the raw images?
They are CCD images in either FITS or OIF (i.e. .imh) files recorded at WIYN.  They are usually 1032x2048 pixels, but
sometimes up to 2080x2048 pixels (i.e. the whole chip).  32 columns of the image (subsection [1:32,*]) is the overscan strip,
which is used for bias correction.
Q: What's the difference between FITS and OIF?
Q: What's an overscan strip?
Q: What's bias correction?  How do I do it?

I usually chdir to some local disk (say, on /d/specher/dolan or /d/uffda3/dolan) and create a directory indicating the date of the observing run (e.g. "mkdir jan99" or oct95a or oct95b or feb98).  Then I make a symbolic link from ~dolan/work to that directory (e.g. ln -s /d/sprecher/dolan/aug99 ~/work/aug99).  This way, I can move the data to another hard disk (say /d/sprecher2/dolan) and just update the pointer in my work directory and I will still have a consistent way to get to the data (e.g. cd ~/work/aug99) without having to remember which disk the data is on.

You can convert from one to the other with WFITS (OIF->FITS) and RFITS (FITS->OIF).  Both are standard IRAF tasks.  When using WFITS, pay careful attention to the "bitpix" parameter.

See the NOAO Hydra manual or the WHYDRA manual or the File format example page for sample dotiraf and dothydra files.

cl> imdel dec98b0001.fits
cl> imdel dec98b0002,dec98b0003,dec980004
cl> ls dec98b*.fits > junkfiles
cl> imdel @junkfiles

One should always check for this problem (and fix it if necessary) before you start using DOHYDRA.  This is because that
task uses the date and time to weight the ThAr spectra in time while computing the dispersion solution.  An incorrect date
will produce an incorrect wavelength solution (sometimes even dramatically wrong).

See the DATE-OBS section to learn to detect/fix this problem.

Mathematically, the ADUs in a CCD image are:
     ADU = gain*(electrons collected) + bias
where of course
     electrons collected  = (photon flux)*(quantum efficiency)

So, to keep ADUs proportional to photons, we must subtract off the bias before applying any multiplicative corrections (e.g. flat fielding)

(1) [I do not do this]  At the telescope, take an image (called a Zero or a Bias) with the camera shutter closed and zero exposure time.  With the shutter closed, there should be no photons.  With zero exposure time, there should be no dark current.  So, all that is left is bias level, bad pixels and read noise.  One must take enough of the Zeros so that, when co-added, the read noise is not significant (often 5-10 suffice).  Then, one subtracts this co-added Zero from all of the other images.  This has the potential problem that if the bias level of the CCD is not constant, one may subtract off the wrong amount.

(2)  [I do this]  Most CCDs (and definitely the MOS CCD) deliberately read out a few too many rows or columns of pixels.  This over-read group pixels (collectively called the "overscan strip") have not been exposed to light, so they contain only the readout bias level.  Thus, one has a mini-Bias image attached to every CCD image.  These pixels can be averaged and subtracted from the remainder of the image.  In the case of the MOS CCD, there are 32 columns of 2048 pixels each.  The 32 columns are averaged together to create a single 1x2048 column.  This is then subtracted from each of the other columns of the image.  This method, unlike method 1, does not correct bad pixels since it is extrapolting the overscan strip to the rest of the pixels on the image.

(Often, folks who do CCD imaging perform version 2 on all images and then apply version 1 for the sake of correcting bad pixels.  This solves the disadvantages of each method.)

I have found that the bias level for the MOS CCD is too variable to use method 1 above.  Thus, I rely on the time-proximity of the overscan strip.  Additionally, I have found that bad pixels do not seem to pose a serious problem for MOS images, so I ignore version 1.

When you have these ~50 data points, IRAF fits a function (usually a cubic spline) to create an empirical relation between pixel and wavelength.  When this relation is extrapolated in from the ThAr time to a science (i.e. Object) exposure, you can calculate a good estimate of the wavelength of any given pixel in the spectrum.

But all this is a big digression!  Back to the line list itself!  IRAF [todo: finish]

I've included a sample thar.dat for you to look at and admire.

    help dohydra

[todo: describe some of the parameters and say why I did not choose the default]

1) Poor quality lines:  These may be blends or just weak features. Whatever the reason, lines like these are usually flagged as fit outliers in many of the apertures, making their usefulness unpredictable.  Since there are many better lines to choose from, I exclude lines like these.

2) Edge-of-chip lines:  This is even more important than the item above.  Obviously there will always be ThAr lines which do not lie on the chip.  For example, you will never see the ThAr feature at 9664.6983 using a setup centered on 5130 Angstroms, clearly.  The problem lines are the ones that are on the chip in some apertures and off the chip in others.  This happens because of the cross-chip dispersion shifts which make the ThAr lines look like they curve across the width of the images.  Lines that are missed consistently by some apertures are removed from the list for the sake of uniformity of the wavelength calibration.  Plus, it's not clear how IRAF treats a line that is half-on and half-off the chip.

It is always best to choose your Canonical ThAr carefully.  If possible, it should be a long exposure=, so the signal-to-noise is good.  Even more important, it should be in the "circle" configuration (e.g. bluecircle).  The reason for the circle configuration is that in this case, all of the fibers are illuminated.  Thus, you can find dispersion solutions for all the fibers at once.  It is very difficult, for example, to extrapolate the solution for fiber 65 to a new ThAr if in you Canonical ThAr fiber 65 was not illuminated!  In this case, you would have to calibrate the new ThAr manually (there's nothing inherently wrong with that; it's just very tedious and tedious means error-prone).

To summarize, if the first DOHYDRA command you issue is on a Sky Flat, you will save time and effort!

a  Autoscale x and y axes
b  Set bottom edge of window
c  Center window at cursor position
d  Shift window down
e  Expand window (mark lower left and upper right of new window)
f  Flip x axis
g  Flip y axis
j  Set left edge of window
k  Set right edge of window
l  Shift window left
m  Autoscale x axis
n  Autoscale y axis
p  Pan x and y axes about cursor
r  Shift window right
t  Set top edge of window
u  Shift window up
x  Zoom x axis about cursor
y  Zoom y axis about cursor
z  Zoom x and y axes about cursor

Notable tasks which do NOT support this command are:
    IMPLOT
    IMEXAMINE

So what IRAF does is step pixel by pixel along the CCD trying to follow the "ridge" (that is the brightest point) of the spectrum image.  It figures out the x and y coordinates of each point along the peak.  In later steps, it will use these coordinates to to sum the pixels perpendicular to the ridge to find the spectral intensity.  See the notes on extraction for more details on how this summing happens.

The problems that may occur are as follows:

1. Adjacent apertures may blend

Unfortunately, there isn't much that can be done about this.  The key here is prevention: I usually expend a fair amount of time at the telescope playing with the spectrograph camera focus trying to optimize the focus so the apertures don't overlap.  If they do, however, it usually happens at the extremes: the top and bottom 100 pixels or so.  If this happens, I sometimes crop the ends off the images.  It is always better to lose a few pixels of data than to possibly include a contaminated radial velocity point in the database!
2. The spectrum may be too faint to trace

You may ask, how does IRAF trace spectra which only have signal-to-noise of about 5?  The answer is that it doesn't even try!!!  Instead, IRAF asks for an image with high signal-to-noise spectra to act as a proxy for the science image.   IRAF traces the spectrum images in this proxy and assumes that the spectra are at about the same place in the science images.  In almost all cases, this proxy should be a Dome Flat taken in the same Hydra configuration as the science data.  If possible, the Dome Flat should be taken near in time to the science image (within a few hours) because the [todo: link CCD chip moves slightly] over time, thus degrading the assumption that the trace of the Dome Flat will be the same as the trace of the science image.  However, this effect is slight on the trace, so it should not worry you as much as [todo: link the effect of the CCD chip drift on the wavelength solution.]

The plot shows the x and y coordinates of points along the ridge of the spectrum image on the CCD.  If the CCD was aligned properly, the x-axis should be very expanded, exaggerating the curvature in the aperture.  IRAF almost always traces perfectly.  Where it doesn't uis usually where there is a cosmic ray, a dead CCD column, or where two apertures may blend.  This latter case is the most dire, but there is little you can do about it other than cropping the image or preventing the blend in the first place at the telescope!

When looking at the plots, you may notice that the sense of the curvature changes across the image.  This is due to curvature in the camera just in front of the CCD.  At least in the Bench Spectrograph Camera (I haven't used the Simmons Camera in a long time), the spectra all curve outwards (I think this is "pincushion" distortion).

Press "q" to quit looking at each plot.

Note that it this step, a normalized flat field is simply created.  It is not actually applied until the science spectra are extracted.

For more details, type in IRAF:
    help dohydra
or go to the IRAF online help page.

If you don't want IRAF to show you this plot at all, type:

    params.f_interactive = no

before starting DOHYDRA (or just type "params")

However, we usually want to do sky subtraction.  The way we usually do this is to assume that the sky background is equally bright in all 95 fibers and subtract an average no-star fiber from all the spectra with stars (see the sky subtraction section).  But if the throughput (that is, the ratio of photons to ADUs) varies from fiber-to-fiber, then our assumption of equal background counts in each aperture is not valid.  So, we need to measure the throughput of each fiber every time we reconfigure the Hydra fiber positions.  Fortunately, we can do this with the Dome Flat image, which is made by simply illuminating all the fibers with a uniform light source - exactly what we need.

In the throughput correction step, each extracted aperture is averaged, yielding a single number for each aperture.  The flat field is multiplied by this number so that when the science spectra are flat fielded later, they will be throughput corrected at the same time.

For more details, type in IRAF:
    help dohydra
or go to the IRAF online help page.

• The NOAO-Wilmarth ThAr atlas

This plastic bound booklet (with its distinctive yellow cover) is (usually) in the library near the librarian's desk on the observatory manuals shelf.  It is a series of plots of a ThAr lamp at higher resolution (and with a narrower spectral window) than WIYN.  This manual was created by Daryl Wilmarth, a superb NOAO staff member who sometimes works at WIYN.
• A WIYN-specific paper map from someone in the department

I have ones centered on 5130, 6650, 6725, 7800 (all echelle) and 7500 (600 l/mm grating) angstroms.  Other people might have their own copies too.
• An electronic copy of the WIYN echelle 5130 angstrom region that we typically use for radial velocities

Three EPS files for you to tape together:
  map-left.eps
  map-center.eps
  map-right.eps

For the 5130 Angstrom setup we use for the WOCS project, I always mark the lines at about 5002, 5068, 5158, and 5231 Angstroms.  These four are sufficient to tie down the mapping, and they are all relatively easy to spot.

1) Look at linearized spectra   (Note: this is identical to the Check dispersion solutions step in Section 3 above)

• Linearize spectra (you do not need to do this for the Canonical ThAr.  This happens automatically in Step 5c)

For example:
cl> dispcor bias0006bluecf.ms "" samedisp=yes verbose=no
• Fire up an image viewer (i.e, run ximtool or saotng or something like that)
• Display the multispec image

For example:
cl> display bias0006bluecf.ms 1
• Look at the ends of the spectra (left and right extrema)
• Do the ThAr emission lines all line up from top to bottom?
• If so, everything should be fine
• If not, there was an error in the dispersion solution for one or more apertures.  This was probably caused by mis-identification of one or more ThAr lines.  To solve this, you must
• Delete the offending multispec image.
• Delete any reduced multispec science data which used this ThAr (check the arctable file to see which ones used it)
• Delete the line ID file (e.g. "database/idbias0006bluecf.ms")
• Re-execute the DOHYDRA command.  It should pick back up by extracting one of the science images (if any) followed by the offending ThAr.

I wrote a short C program to do this.  It looks at the output line ID file and counts how many of each line was missed.  The program is in ~dolan/work/iraf/bin/lines (the source code is in the same directory in a file called lines.c)

To run it (assuming my iraf/bin directory is in your $path, as indicated in the Setup portion of this cookbook):

[For example]
% lines database/idbias0006bluecf.ms
[And it outputs, after a bit of debugging notes, something like this: ]
5247.6547   1
5240.1968  14
5239.5520   0
5238.8137   1
5231.1597   0
5221.2710   0
5219.1099   0
5213.3492   0
5211.2305   0
5209.7246   0
5199.1637   1
5195.8136   0
5194.4576   2
5187.7462   0
5176.9610   1
5175.3248   0
5165.7728   2
5163.4584   0
5162.2846   0
5161.5396   3
5158.6041   0
5154.2430   0
5145.3083   0
5143.9165   5
5141.7827   0
5134.7460   2
5128.4897   7
5115.0448   0
5100.6211   2
5096.4848   0
5081.4462   6
5067.9737   0
5066.1355   4
5064.9454  11
5062.0371   0
5055.3473  19
5051.8887   2
5050.7842   1
5049.7960   1
5047.0434   0
5044.7195   0
5039.2303   0
5029.8916   8
5028.6556   0
5019.8062   0
5017.1628   0
5009.3344   1
5002.0972   0

The first column should be just a reverse-sorted version of the thar.dat file, while the second column is the number of times that line was rejected from the pixel-to-wavelength mapping function fit (with a maximum of 95).  Clearly, some are missed more than others (e.g. 5055.3473).  Usually, these are the weakest lines, so they are hurt by poor signal-to-noise.  However, if a line is missed too often, it likely suffers from a blend or some other irregularity that make it a poor candidate for use in the dispersion solution.  In the case of the 5130 Angstrom setup, I have already culled the worst offenders from the thar.dat file, so it is unlikely you could improve the solution by editing the line list.  But with another spectral region, you may find that the lines program is a useful tool for identifying lines which should be removed from a line list.  In particular, you should remove lines from the red and blue extremese which lines reports as frequently missed, since these are probably falling right off the edge of the CCD image due to the spectrograph dispersion curvature!

Finally, using the interpolated mapping, the spectrum is linearized.  That is, the spectrum is resampled so each pixel represents a fixed step in wavelength.  This is important so the spectrum can be used in FXCOR later.  Once the spectrum is resampled, the dispersion solution is written into the image header (really, just two numbers since the spectrum is now linear).

The problem is that the sky subtraction task in IRAF is somewhat fragile.  It has some problems that (seemingly) unpredictably cause it to crash.  If the routine crashes, it often corrupts the spectrum that it is in the middle of processing.  So, if there is no backup copy, you have to redo the DOHYDRA command to re-extract the spectrum, which is very annoying.  Thus, we don't want to use the sky subtraction within DOHYDRA (really, dohydra just calls the same task as we call outside of dohydra later: skysub) since we don't have the opportunity to make backups.

Additionally, sometimes (but rarely) we decide after the fact that the night was dark enough that the sky subtraction made the spectra worse instead of better (by adding cosmic ray blemishes or simple noise).  In these cases, it's best to have a non-sky-subtracted copy to which we can refer.

By the way, if there is a glitch and skysub trashes some of your spectra, to step back, you should delete the science file and the "sky" file, then copy the backup science file into the current directory.  Then, once you've typed "flpr" a half-dozen times or restarted IRAF, you can try again...  Good luck, you may need it!

In particular you should watch out for the following types of defects:

• Bright spectra

If one of the sky spectra is much brighter than the others, it is probably contaminated by a star or some other non-terrestrial source.  Since we just want light from the uniform background, fibers with light other than sky shoul;d be excluded.
• Noisy spectra

For reasons I don't fully understand yet, sometimes you get a sky spectrum which appears much more noisy than the others.  I think I might have something to do with the flat fielding, but I don't really know.  Since spectra like these don't really contribute much to the resulting sky spectrum, they should be excluded.
• Huge cosmic rays

If you see a cosmic ray that is so strong that the continuum is depressed around it, you may want to consider deleteing it.  This case requires more thought than the above case, since usually most of spectrum is OK, but one little spot got messed up.  If you have plenty of good spectra (like maybe 6 or more), I say just delete it.  But if you don't have many spectra, so that having one extra might make a real difference in the resulting sky spectrum's signal-to-noise, then you might want to consider leaving it in.  It can be a tough choice.
• Other cases?

You'll have to judge for yourself.  In general, if you have plenty of sky spectra, don't take any chances: just delete it.

In this case, avsigclip works well because there are usually enough sky fibers that a robust measure of the mean and sigma can be computed.  When you get down to about 4 or fewer elements, the avsigclip algorithm gets flaky.