IRAF help page for package noao.imred.specred, program dofibers

Extraction

The identification of the spectra in the two dimensional images and their scattered light subtraction and extraction to one dimensional spectra in multispec format is accomplished using the tasks from the apextract package. The first parameters through nsubaps control the extractions.

The dispersion line is that used for finding the spectra, for plotting in the aperture editor, and as the starting point for tracing. The default value of INDEF selects the middle of the image. The aperture finding, adjusting, editing, and tracing operations also allow summing a number of dispersion lines to improve the signal. The number of lines is set by the nsum parameter.

The order parameter defines whether the order of the aperture identifications in the aperture identification file (or the default sequential numbers if no file is used) is in the same sense as the image coordinates (increasing) or the opposite sense (decreasing). If the aperture identifications turn out to be opposite to what is desired when viewed in the aperture editing graph then simply change this parameter.

The basic data output by the spectral extraction routines are the one dimensional spectra. Additional information may be output when the extras option is selected and the cleaning or variance weighting options are also selected. In this case a three dimensional image is produced with the first element of the third dimension being the cleaned and/or weighted spectra, the second element being the uncleaned and unweighted spectra, and the third element being an estimate of the sigma of each pixel in the extracted spectrum. Currently the sigma data is not used by any other tasks and is only for reference.

The initial step of finding the fiber spectra in the aperture reference image consists of identifying the peaks in a cut across the dispersion, eliminating those which are closer to each other than the minsep distance, and then keeping the specified nfibers highest peaks. The centers of the profiles are determined using the center1d algorithm which uses the width parameter.

Apertures are then assigned to each spectrum. The initial edges of the aperture relative to the center are defined by the lower and upper parameters. The trickiest part of assigning the apertures is relating the aperture identification from the aperture identification file to automatically selected fiber profiles. The first aperture id in the file is assigned to the first spectrum found using the order parameter to select the assignment direction. The numbering proceeds in this way except that if a gap greater than a multiple of the maxsep parameter is encountered then assignments in the file are skipped under the assumption that a fiber is missing (broken). If unassigned fibers are still visible in a flat field, either by design or by scattered light, the unassigned fibers can be included in the number of fibers to find and then the unassigned (negative beam number) apertures are excluded from any extraction. For more on the finding and assignment algorithms see apfind.

The initial apertures are the same for all spectra but they can each be automatically resized. The automatic resizing sets the aperture limits at a fraction of the peak relative to the interfiber minimum. The default ylevel is to resize the apertures to 5% of the peak. See the description for the task apresize for further details.

The user is given the opportunity to graphically review and adjust the aperture definitions. This is recommended. As mentioned previously, the correct identification of the fibers is tricky and it is fundamentally important that this be done correctly; otherwise the spectrum identifications will not be for the objects they say. An important command in this regard is the 'o' key which allows reordering the identifications based on the aperture identification file. This is required if the first fiber is actually missing since the initial assignment begins assigning the first spectrum found with the first entry in the aperture file. The aperture editor is a very powerful tool and is described in detail as apedit.

The next set of parameters control the tracing and function fitting of the aperture reference positions along the dispersion direction. The position of a spectrum across the dispersion is determined by the centering algorithm (see center1d) at a series of evenly spaced steps, given by the parameter t_step, along the dispersion. The step size should be fine enough to follow position changes but it is not necessary to measure every point. The fitted points may jump around a little bit due to noise and cosmic rays even when summing a number of lines. Thus, a smooth function is fit. The function type, order, and iterative rejection of deviant points is controlled by the other trace parameters. For more discussion consult the help pages for aptrace and icfit. The default is to fit a cubic spline of three pieces with a single iteration of 3 sigma rejection.

The actual extraction of the spectra by summing across the aperture at each point along the dispersion is controlled by the next set of parameters. The default extraction simply sums the pixels using partial pixels at the ends. The options allow selection of a weighted sum based on a Poisson variance model using the readnoise and gain detector parameters. Note that if the clean option is selected the variance weighted extraction is used regardless of the weights parameter. The sigma thresholds for cleaning are also set in the params parameters. For more on the variance weighted extraction and cleaning see apvariance and approfiles as well as apsum.

The last parameter, nsubaps, is used only in special cases when it is desired to subdivide the fiber profiles into subapertures prior to dispersion correction. After dispersion correction the subapertures are then added together. The purpose of this is to correct for wavelength shifts across a fiber.

Scattered Light Subtraction

Scattered light may be subtracted from the input two dimensional image as the first step. This is done using the algorithm described in apscatter. This can be important if there is significant scattered light since the flat field/throughput correction will otherwise be incorrect. The algorithm consists of fitting a function to the data outside the defined apertures by a specified buffer at each line or column across the dispersion. The function fitting parameters are the same at each line. Because the fitted functions are independent at each line or column a second set of one dimensional functions are fit parallel to the dispersion using the evaluated fit values from the cross-dispersion step. This produces a smooth scattered light surface which is finally subtracted from the input image. Again the function fitting parameters are the same at each line or column though they may be different than the parameters used to fit across the dispersion.

The first time the task is run with a particular flat field (or aperture reference image if no flat field is used) the scattered light fitting parameters are set interactively using that image. The interactive step selects a particular line or column upon which the fitting is done interactively with the icfit commands. A query is first issued which allows skipping this interactive stage. Note that the interactive fitting is only for defining the fitting functions and orders. When the graphical icfit fitting is exited (with 'q') there is a second prompt allowing you to change the buffer distance (in the first cross-dispersion stage) from the apertures, change the line/column, or finally quit.

The initial fitting parameters and the final set parameters are recorded in the apscat1 and apscat2 hidden parameter sets. These parameters are then used automatically for every subsequent image which is scattered light corrected.

The scattered light subtraction modifies the input 2D images. To preserve the original data a copy of the original image is made with the same root name and the word "noscat" appended. The scattered light subtracted images will have the header keyword "APSCATTE" which is how the task avoids repeating the scattered light subtraction during any reprocessing. However if the redo option is selected the scattered light subtraction will also be redone by first restoring the "noscat" images to the original input names.

Flat Field and Fiber Throughput Corrections

Flat field corrections may be made during the basic CCD processing; i.e. direct division by the two dimensional flat field observation. In that case do not specify a flat field spectrum; use the null string "". The dofibers task provides an alternative flat field response correction based on division of the extracted object spectra by the extracted flat field spectra. A discussion of the theory and merits of flat fielding directly verses using the extracted spectra will not be made here. The dofibers flat fielding algorithm is the recommended method for flat fielding since it works well and is not subject to the many problems involved in two dimensional flat fielding.

In addition to correcting for pixel-to-pixel response the flat field step also corrects for differences in the fiber throughput. Thus, even if the pixel-to-pixel flat field corrections have been made in some other way it is desirable to use a sky or dome flat observation for determining a fiber throughput correction. Alternatively, a separately derived throughput file may be specified. This file consists of the aperture numbers (the same as used for the aperture reference) and relative throughput numbers.

The first step is extraction of the flat field spectrum, if specified, using the reference apertures. Only one flat field is allowed so if multiple flat fields are required the data must be reduced in groups. After extraction one or more corrections are applied. If the fitflat option is selected (the default) the extracted flat field spectra are averaged together and a smooth function is fit. The default fitting function and order are given by the parameters f_function and f_order. If the parameter f_interactive is "yes" then the fitting is done interactively using the fit1d task which uses the icfit interactive fitting commands.

The fitted function is divided into the individual flat field spectra to remove the basic shape of the spectrum while maintaining the relative individual pixel responses and any fiber to fiber differences. This step avoids introducing the flat field spectrum shape into the object spectra and closely preserves the object counts.

If a throughput image is available (an observation of blank sky usually at twilight) it is extracted. If no flat field is used the average signal through each fiber is computed and this becomes the response normalization function. Note that a dome flat may be used in place of a sky in the sky flat field parameter for producing throughput only corrections. If a flat field is specified then each sky spectrum is divided by the appropriate flat field spectrum. The total counts through each fiber are multiplied into the flat field spectrum thus making the sky throughput of each fiber the same. This correction is important if the illumination of the fibers differs between the flat field source and the sky. Since only the total counts are required the sky or dome flat field spectra need not be particularly strong though care must be taken to avoid objects.

Instead of a sky flat or other throughput image a separately derived throughput file may be used. It may be used with or without a flat field.

The final step is to normalize the flat field spectra by the mean counts of all the fibers. This normalization step is simply to preserve the average counts of the extracted object and arc spectra after division by the response spectra. The final relative throughput values are recorded in the log and possibly printed on the terminal.

These flat field response steps and algorithm are available as a separate task called msresp1d.

Dispersion Correction

Dispersion corrections are applied to the extracted spectra if the dispcor parameter is set. This can be a complicated process which the dofibers task tries to simplify for you. There are three basic steps involved; determining the dispersion functions relating pixel position to wavelength, assigning the appropriate dispersion function to a particular observation, and resampling the spectra to evenly spaced pixels in wavelength.

The comparison arc spectra are used to define dispersion functions for the fibers using the tasks identify and reidentify. The interactive identify task is only used on the central fiber of the first arc spectrum to define the basic reference dispersion solution from which all other fibers and arc spectra are automatically derived using reidentify.

The set of arc dispersion function parameters are from identify and reidentify. The parameters define a line list for use in automatically assigning wavelengths to arc lines, a parameter controlling the width of the centering window (which should match the base line widths), the dispersion function type and order, parameters to exclude bad lines from function fits, and parameters defining whether to refit the dispersion function, as opposed to simply determining a zero point shift, and the addition of new lines from the line list when reidentifying additional arc spectra. The defaults should generally be adequate and the dispersion function fitting parameters may be altered interactively. One should consult the help for the two tasks for additional details of these parameters and the operation of identify.

Generally, taking a number of comparison arc lamp exposures interspersed with the program spectra is sufficient to accurately dispersion calibrate multifiber spectra. However, there are some other calibration options which may be of interest. These options apply additional calibration data consisting either of auxiliary line spectra, such as from dome lights or night sky lines, or simultaneous arc lamp spectra taken through a few fibers during the object exposure. These options add complexity to the dispersion calibration process.

When only arc comparison lamp spectra are used, dispersion functions are determined independently for each fiber of each arc image and then assigned to the matching fibers in the program object observations. The assignment consists of selecting one or two arc images to calibrate each object image. When two bracketing arc spectra are used the dispersion functions are linearly interpolated (usually based on the time of the observations).

If taking comparison exposures is time-consuming, possibly requiring reconfiguration to illuminate the fibers, and the spectrograph is expected to be fairly stable apart from small shifts, there are two mutually exclusive methods for monitoring shifts in the dispersion zero point from the basic arc lamp spectra other than taking many arc lamp exposures. One is to use some fibers to take a simultaneous arc spectrum while observing the program objects. The fibers are identified by aperture or beam numbers. The second method is to use auxiliary line spectra, such as mercury lines from the dome lights. These spectra are specified with an auxiliary shift arc list, arc2.

When using auxiliary line spectra for monitoring zero point shifts one of these spectra is plotted interactively by identify with the reference dispersion function from the reference arc spectrum. The user marks one or more lines which will be used to compute zero point wavelength shifts in the dispersion functions automatically. The actual wavelengths of the lines need not be known. In this case accept the wavelength based on the reference dispersion function. As other observations of the same features are made the changes in the positions of the features will be tracked as zero point wavelength changes such that wavelengths of the features remain constant.

When using auxiliary line spectra the only arc lamp spectrum used is the initial arc reference spectrum (the first image in the arcs1 list). The master dispersion functions are then shifted based on the spectra in the arcs2 list (which must all be of the same type). The dispersion function assignments made by refspectra using either the arc assignment file or based on header keywords is done in the same way as described for the arc lamp images except using the auxiliary spectra.

If simultaneous arcs are used the arc lines are reidentified to determine a zero point shift relative to the comparison lamp spectra selected, by refspectra, of the same fiber. A linear function of aperture position on the image across the dispersion verses the zero point shifts from the arc fibers is determined and applied to the dispersion functions from the assigned calibration arcs for the non-arc fibers. Note that if there are two comparison lamp spectra (before and after the object exposure) then there will be two shifts applied to two dispersion functions which are then combined using the weights based on the header parameters (usually the observation time).

The arc assignments may be done either explicitly with an arc assignment table (parameter arctable) or based on a header parameter. The task used is refspectra and the user should consult this task if the default behavior is not what is desired. The default is to interpolate linearly between the nearest arcs based on the Julian date (corrected to the middle of the exposure). The Julian date and a local Julian day number (the day number at local noon) are computed automatically by the task setjd and recorded in the image headers under the keywords JD and LJD. In addition the universal time at the middle of the exposure, keyword UTMIDDLE, is computed by the task setairmass and this may also be used for ordering the arc and object observations.

The last step of dispersion correction (resampling the spectrum to evenly spaced pixels in wavelength) is optional and relatively straightforward. If the linearize parameter is no then the spectra are not resampled and the nonlinear dispersion information is recorded in the image header. Other IRAF tasks (the coordinate description is specific to IRAF) will use this information whenever wavelengths are needed. If linearizing is selected a linear dispersion relation, either linear in the wavelength or the log of the wavelength, is defined once and applied to every extracted spectrum. The resampling algorithm parameters allow selecting the interpolation function type, whether to conserve flux per pixel by integrating across the extent of the final pixel, and whether to linearize to equal linear or logarithmic intervals. The latter may be appropriate for radial velocity studies. The default is to use a fifth order polynomial for interpolation, to conserve flux, and to not use logarithmic wavelength bins. These parameters are described fully in the help for the task dispcor which performs the correction. The interpolation function options and the nonlinear dispersion coordinate system is described in the help topic onedspec.package.

Sky Subtraction

Sky subtraction is selected with the skysubtract processing option. The sky spectra are selected by their aperture and beam numbers and combined into a single master sky spectrum which is then subtracted from each object spectrum. If the skyedit option is selected the sky spectra are plotted using the task specplot. By default they are superposed to allow identifying spectra with unusually high signal due to object contamination. To eliminate a sky spectrum from consideration point at it with the cursor and type 'd'. The last deleted spectrum may be undeleted with 'e'. This allows recovery of incorrect or accidental deletions.

The sky combining algorithm parameters define how the individual sky fiber spectra, after interactive editing, are combined before subtraction from the object fibers. The goals of combining are to reduce noise, eliminate cosmic-rays, and eliminate fibers with inadvertent objects. The common methods for doing this to use a median and/or a special sigma clipping algorithm (see scombine for details). The scale parameter determines whether the individual skys are first scaled to a common mode. The scaling should be used if the throughput is uncertain, but in that case you probably did the wrong thing in the throughput correction. If the sky subtraction is done interactively, i.e. with the skyedit option selected, then after selecting the spectra to be combined a query is made for the combining algorithm. This allows modifying the default algorithm based on the number of sky spectra selected since the "avsigclip" rejection algorithm requires at least three spectra.

The combined sky spectrum is subtracted from only those spectra specified by the object aperture and beam numbers. Other spectra, such as comparison arc spectra, are retained unchanged. One may include the sky spectra as object spectra to produce residual sky spectra for analysis. The combined master sky spectra may be saved if the saveskys parameter is set. The saved sky is given the name of the object spectrum with the prefix "sky".

1. The following example uses artificial data and may be executed at the terminal (with IRAF V2.10). This is also the sequence performed by the test procedure "demos dohydra" from the hydra package..

sp> hydra
hy> demos mkhydra
Creating image demoobj ...
Creating image demoflat ...
Creating image demoarc ...
hy> bye
sp> type demoapid
===> demoapid <===
36 1
37 0
38 1
39 1
41 0
42 1
43 1
44 0
45 1
46 -1
47 0
48 1
sp> specred.verbose = yes
sp> dofibers demoobj apref=demoflat flat=demoflat arcs1=demoarc 
>>> fib=12 apid=demoapid width=4. minsep=5. maxsep=7. clean- splot+
Set reference apertures for demoflat
Resize apertures for demoflat?  (yes):
Edit apertures for demoflat?  (yes):

Fit curve to aperture 36 of demoflat interactively  (yes):

Fit curve to aperture 37 of demoflat interactively  (yes): N
Create response function demoflatdemoad.ms
Extract flat field demoflat
Fit and ratio flat field demoflat

Create the normalized response demoflatdemoad.ms
demoflatdemoad.ms -> demoflatdemoad.ms  using bzero: 0.
    and bscale: 1.000001
    mean: 1.000001  median: 1.052665  mode: 1.273547
    upper: INDEF  lower: INDEF
Average fiber response:
1.  1.151023
2.  0.4519709
3.  1.250614
4.  1.287281
5.  1.271358
6.  0.6815334
7.  1.164336
8.  0.7499605
9.  1.008654
10.  1.053296
11.  0.929967
Extract arc reference image demoarc
Determine dispersion solution for demoarc




REIDENTIFY: NOAO/IRAF V2.10BETA valdes@puppis Tue 16:01:07 11-Feb-92
  Reference image = d....ms.imh, New image = d....ms, Refit = yes
     Image Data Found    Fit Pix Shift User Shift  Z Shift     RMS
d....ms - Ap 41 16/20  16/16   0.00796     0.0682  8.09E-6    3.86
Fit dispersion function interactively? (no|yes|NO|YES) (NO): y

d....ms - Ap 41 16/20  16/16   0.00796     0.0682  8.09E-6    3.86
d....ms - Ap 39 19/20  19/19     0.152        1.3  1.95E-4    3.89
Fit dispersion function interactively? (no|yes|NO|YES) (yes): N
d....ms - Ap 39 19/20  19/19     0.152        1.3  1.95E-4    3.89
d....ms - Ap 38 18/20  18/18     0.082      0.697  9.66E-5    3.64
d....ms - Ap 37 19/20  19/19    0.0632      0.553  1.09E-4    6.05
d....ms - Ap 36 18/20  18/18    0.0112     0.0954  1.35E-5    4.12
d....ms - Ap 43 17/20  17/17    0.0259      0.221  3.00E-5    3.69
d....ms - Ap 44 19/20  19/19     0.168       1.44  2.22E-4    4.04
d....ms - Ap 45 20/20  20/20      0.18       1.54  2.35E-4    3.95
d....ms - Ap 47 18/20  18/18  -2.02E-4    0.00544  9.86E-6     4.4
d....ms - Ap 48 16/20  16/16   0.00192     0.0183  1.44E-6    3.82
Dispersion correct demoarc
d....ms.imh: w1 = 5748.07..., w2 = 7924.62..., dw = 8.50..., nw = 257
  Change wavelength coordinate assignments? (yes|no|NO): n
Extract object spectrum demoobj
Assign arc spectra for demoobj
[demoobj] refspec1='demoarc'
Dispersion correct demoobj
demoobj.ms.imh: w1 = 5748.078, w2 = 7924.622, dw = 8.502127, nw = 257
Sky subtract demoobj:  skybeams=0
Edit the sky spectra? (yes):

Sky rejection option (none|minmax|avsigclip) (avsigclip):
demoobj.ms.imh:
Splot spectrum? (no|yes|NO|YES) (yes):
Image line/aperture to plot (1:) (1):