21.4 Descriptions of Calibration Steps

• Two-dimensional image reduction, including basic 2-D reduction, cosmic ray rejection and image co-addition.
• Processing of the contemporaneously obtained wavecal.
• Two-dimensional and one-dimensional spectral extraction, with flux and wavelength calibration. Within each component, the individual steps are listed alphabetically, because the order in which they are performed can change for different types of data (e.g., CCD or MAMA, spectroscopic or imaging, CR-split or not).

  More detailed descriptions can be found in a series of Instrument Science Reports (ISRs) that discuss the pipeline. Be aware, however, that while these reports describe the original design of the pipeline and the associated algorithms in detail, they do not always contain information concerning later modifications. Over time these reports will be updated to include a more complete description of the pipeline. In the meantime we refer you to the STIS WWW page, where a history of the important changes to the calstis pipeline code is maintained.

BIASCORR

BLEVCORR

Because the electronic bias level can vary with time and temperature, its value is determined from the overscan region in the particular exposure being processed. A raw STIS CCD taken in full frame unbinned mode will have 20 rows of virtual parallel overscan in the AXIS2, or image y, direction, which is created by over-clocking the readout of each line past its physical extent, and 19 leading and trailing columns of serial physical overscan in the AXIS1 or image x direction, which arise from unilluminated pixels on the CCD. Thus the size of the uncalibrated and unbinned full frame CCD image is 1062 (serial) by 1044 (parallel) pixels, with 1024 x 1024 exposed science pixels.

Only the serial (physical) overscan is used for the overscan bias level determination; the virtual parallel overscan is not used. A line-by-line subtraction is performed in the following way. An initial value of the electronic bias level, or overscan, is determined for each line of the image, using only the physical serial overscan, and a function, currently a straight line, is fit to these values as a function of image line. The actual overscan value subtracted from an image line is the value of the linear fit at that image line. The initial value for each line is found by taking the median of a predetermined subset of the trailing serial overscan pixels. Currently, that region includes most of the trailing overscan region, however the first and last two pixels are skipped, as they have been shown to be subject to problems, and pixels flagged as bad in the input data quality flag are also skipped. The region used changes with binning or subarray use (see Table 21.2). The mean value of all overscan levels is computed, and the mean is written to the output SCI extension header as MEANBLEV.

In addition to subtracting the electronic bias level, the BLEVCORR step also trims the image of overscan. The sizes of the overscan regions depend on binning and whether the image is full-frame or a subimage. The locations of the overscan regions depend on which amplifier was used for readout. The number of pixels to trim off each side of the image (before accounting for readout amplifier) is given in Table 21.3. The values of NAXIS1, NAXIS2, BINAXIS1, and BINAXIS2 are obtained from image header keywords. Because the binning factor does not divide evenly into 19 and 1062, when on-chip pixel binning is used the raw image produced will contain both pure overscan pixels, overscan plus science pixels and science pixels. The calstis pipeline will only calibrate pixel binnings of 1, 2, and 4 in either AXIS1 or AXIS2.

The CRPIXi and LTVi keywords are updated in the output; these depend on the offset from removing the overscan.

Raw Image Pixels Used to Determine Line by Line Bias Level

Columns in Raw Image	Binning
2 through 16	Unbinned
2 through 8	All other supported binnings

Pixels Trimmed During CCD Bias Level Correction for Amp D

Side	Full Image	Subarray Images	Binned Images
Right	19	18	(19 + 1) / BINAXIS1
Left	19	18	NAXIS1 - (1024 / BINAXIS1 - 1) - Right
Top	0	0	0
Bottom	20	0	NAXIS2 - 1024 / BINAXIS2

CRCORR

The CRCORR step does the following:

• Forms a stack of images to be combined (the CR-split or repeatobs exposures in the input file).
• Forms an initial guess image (minimum or median).
• Forms a summed CR-rejected image, using the guess image to reject high and low values in the stack, based on sigma and the radius parameter which governs whether to reject neighboring pixels to pixels identified as cosmic ray (see below).
• Performs one or more rejection cycles, using different (usually decreasing) rejection thresholds, producing a new guess image at each iteration.
• Produces a final cosmic ray rejected image, including science, data quality and error extensions, which is the weighted sum of the input images.
• Flags the data quality arrays of the individual (non-CR-rejected) input files to indicate where an outlier has been found (pixels which were rejected because of cosmic ray hits can be identified by looking for data quality bit = 14 in the _flt.fits file). The cosmic-ray-rejected image is created by setting the value at each pixel to the sum of the values of all good pixels in the stack whose values are within plus or minus sigmas*noise of the initial guess image. Deviant (out of range) stack pixels are flagged as cosmic ray impacted by setting their stack data quality flags to 8192 in the input file.

  The value of noise (in DN) is computed as:

  Where:

• DN = the data number of the stack pixel value.
• READNSE is the read noise in electrons, read from the primary header.
• ATODGAIN is the calibrated conversion from electrons to DN, read from the primary header.
• SCALENSE is an input parameter, read from the crrejtab -calibration -reference file. The sigmas parameter is read from the crrejtab calibration reference file; sigmas is a string, e.g., sigmas = "4,3". The number of entries in the string dictates the number of iterations to be performed (in this example two) and the values in the string indicate the value of sigmas for each iteration. In this example, stack values that deviate from the guess image value by more than 4*noise in the first iteration are considered to be outliers and are excluded from the average on the first iteration, and an improved guess image is formed. A second iteration is then performed in which sigmas is set to 3 and good stack values disparate by more than +/- 3*noise from the guess image are excluded when determining the average. In each iteration, if RADIUS is >=1, then pixels neighboring rejected pixels are also excluded in forming the average.

  SCALENSE is a string containing a multiplicative factor in the noise relation. If SCALENSE = "2.0", then the term 0.02*value is added to the noise. This term accounts for multiplicative effects, such as would be expected if this rejection were applied to Flatfield ed data. It is important to include to properly flag well-exposed regions (such as the centers of stars, where jitter from the telescope may slightly change the pointing from image to image).

  The combination of the individual CRSPLIT or REPEATOBS exposures into a single cosmic-ray rejected frame is performed early in the calstis flow. The cosmic ray rejection is performed after each exposure has had its data quality file initialized and the overscan bias level subtraction (BLEVCORR) performed upon it, but prior to subtraction of a bias frame (BIASCORR), dark (DARKCORR) and flatfielding of the data (FLATCORR). The CR-rejected, bias level subtracted image is then passed through the remainder of the two-dimensional image reduction (calstis-1) to produce a flatfield ed CR-rejected image (*_crj.fits). This CR-rejected flatfield ed image is then passed through the subsequent processing steps in calstis. If EXPSCORR is set to PERFORM (see below), then the individual flatfielded but not cosmic ray rejected exposures are also produced.

The mean of the dark values subtracted is written to the SCI extension header with the keyword MEANDARK. For CCD data, the dark image is multiplied by the exposure time and divided by the atodgain (from the CCD parameters table) before subtracting. For MAMA data, the dark image will be multiplied by the exposure time before subtracting; it will also be convolved with the Doppler smoothing function if DOPPCORR is PERFORM.

DOPPCORR

The first step is to compute an array containing the Doppler smearing function. The expression below gives the computed Doppler shift, where the time t begins with the value of the header keyword EXPSTART and is incremented in one-second intervals up to EXPSTART + EXPTIME inclusive. At each of these times, the Doppler shift in unbinned pixels is computed as:

The value of shift is rounded to the nearest integer.

DQICORR

EXPSCORR

FLATCORR

• pfltfile - This flat is a configuration (grating, central wavelength and detector) dependent pixel-to-pixel flatfield image, from which any large-scale sensitivity variations have been removed (i.e., it will have a local mean value of unity across its entirety). Such configuration dependent flats are expected to be produced infrequently, perhaps once per year.
• dfltfile - This file is a delta flat which gives the changes in the small scale flatfield response relative to the pixel to pixel flat (pfltfile). Delta flats will be taken relatively frequently (approximately monthly, though less frequently for the MAMAs); there will be a single delta flat for each detector, CCD, NUV-MAMA, and FUV-MAMA. They will be used only if needed.
• lfltfile - This flat is a subsampled image containing the large-scale sensitivity variation across the detector. It is usually grating- and central wavelength-dependent (for spectroscopic data) and aperture (filter) dependent for imaging data. To flatfield science data, calstis creates a single flatfield image from these three files¹ as described below and then divides the science image by the flat so created. The pixels of the science data quality file are updated to reflect bad pixels in the input reference files and the errors in the science data are updated to reflect the application of the flat. Blank values of pfltfile, dfltfile, or lftfile in the science data, indicate that type of flat is not to be used.

  To create the single combined flatfield file, calstis first expands the large-scale sensitivity flat (lftfile) to full format. The pixel-to-pixel flat, delta flat, and expanded low-order flat are then multiplied together. For MAMA data, the product of the flatfield images will be convolved with the Doppler smoothing function if DOPPCORR = PERFORM. If a subarray or binning was used, after taking the product of all the flatfields that were specified, a subset is taken and binned if necessary to match the uncalibrated image, and the uncalibrated data are then divided by the binned subset.

The method used is similar to 2-D rectification of spectroscopic data (see "X2DCORR" on page 21-29). For each pixel in the output rectified image, the corresponding point is found in the input distorted image, and bi-linear interpolation is used on the four nearest pixels to determine the value to assign to the output. Mapping from an output pixel back into the input images specified by two-dimensional Chebyshev polynomials stored in the format generated by the IRAF gsurfit package.

GLINCORR and LFLGCORR

• Pore paralysis in the micro channel plates arises when charge cannot flow rapidly enough to replenish channels whose electrons have been depleted due to high local photon rates. This depletion produces a local non-linearity. The local count rate is roughly linear up to counts rates of ~200 counts/second/pixel and then turns directly over, showing an inverted V shape. Thus it is not possible to reliably correct for or flag pixels which have exceeded the local linearity limit in the pipeline (because the relation is bi-valued).
• The electronic processing circuitry has a dead-time of roughly 350 nano-seconds between pulses; thus at global count rates (across the detector) of 300,000 counts (pulses) per second, the electronic circuity counts roughly 90% of the pulses.
• The MIE electronics and flight software can process at most 300,000 pulses per second (i.e., it is matched to the expected global count rate performance of the electronic circuitry). At count rates higher than this, the MIE will still count only 300,000 pulses per second-this represents a hard cutoff beyond which no information is available to allow correction to the true count rate. In practice, at count rates approaching 270,000 counts/sec the flight software begins losing counts due to the structure of its data buffers. Further work is needed to understand this effect. For subarrays, the hard cutoff limit of the MIE electronics and software will differ from that for full frame processing, but will still be dependent on the total global rate in addition to the rate within the subarray. The calstis pipeline currently applies the full frame correction to subarray data.

  The global count rate (across the entire detector) is determined as part of the bright object protection sequence and is passed down with the exposure as a header keyword, GLOBRATE in the science header. If either GLINCORR or LFLGCORR is PERFORM, the global count rate will be checked; a correction for global non-linearity applied if GLINCORR is PERFORM, using the parameters GLOBAL_LIMIT, LOCAL_LIMIT, TAU, and EXPAND read from the mlintab reference table.

  If the value of the SCI extension header keyword GLOBRATE is greater than GLOBAL_LIMIT, the keyword GLOBLIM in the SCI extension header will be set to EXCEEDED; otherwise, GLOBLIM will be set to NOT-EXCEEDED, and a correction factor will be computed and multiplied by each pixel in the science image and error array. The correction factor is computed by iteratively solving GLOBRATE = X * exp (-TAU * X) for X, where X is the true count rate. This algorithm has not yet been updated to account for the linearity effects from the flight software data buffer management.

  If LFLGCORR is PERFORM, each pixel in the science image is also compared with the product of LOCAL_LIMIT and the exposure time EXPTIME. That count rate limit is then adjusted for binning by dividing by the pixel area in high-res pixels. If the science data value is larger than that product, that pixel and others within a radius of EXPAND high-res pixels are flagged as nonlinear. Because our understanding of the MAMA processing electronics is currently incomplete, accurate fluxes (global linearity) at count rates exceeding 270,000 count/sec cannot be expected from the calstis pipeline.

The binning of the uncalibrated image is determined from the LTM1 and LTM2 keywords in the SCI extension header of the raw data file. LTMi = 1 implies the reference pixel size, and LTMi = 2 means the pixels are subsampled into high-res format. In this step, if either or both axes are high-res, they will be binned down to low-res. The binning differs from binning reference files to match an uncalibrated image, in that the pixel values in this step are summed rather than averaged.

PHOTCORR

RPTCORR

SHADCORR

21.4.2 WAVECAL Processing

Basic two-dimensional image reduction (basic2d) is first applied to the wavecal. For CCD data taken with the hole in the mirror (HITM) system the external shutter is ordinarily open, so the detector will have been exposed to radiation from both the science target and the line lamp. In this case, the next step is to scale the flatfield ed science image by the ratio of exposure times and subtract it from the flatfield ed wavecal. Two-dimensional rectification (x2d, see X2DCORR below) is then applied to the flatfielded (and possibly science subtracted) wavecal.

Because wavecal data are not CR-split, cosmic rays must be identified and eliminated by looking for outliers within columns, i.e., in the cross-dispersion direction. Since the data have been rectified, the image can be collapsed along columns to get a long-slit integrated spectrum or along rows to get an outline of the slit (in the cross-dispersion direction).

The shift in the dispersion direction is found by cross-correlating the observed wavecal spectrum with a template spectrum. In the cross dispersion direction, edge location is used for medium and long slits, and cross correlation is used for very short, echelle, slits. The long slits have two occulting bars, and it is the edges of these bars that are used for finding the location. Edges are found by convolving the cross-dispersion profile with the array [-1, 0, +1]. The peak in cross correlation and the edge location are obtained to subpixel level by fitting a quadratic polynomial to the three pixels nearest the extremum.

The shifts are initially measured in units of pixels of the wavecal image, but they are then scaled (depending on the binning of the wavecal) to the reference pixel size (unbinned CCD or low-res MAMA). They are subsequently written to the extension header of the 2-D rectified wavecal in the keywords SHIFTA1 (the shift in pixels along AXIS1, or dispersion direction) and SHIFTA2 (the shift in pixels along AXIS2, or spatial direction). The SHIFTA1 and SHIFTA2 keyword values are also copied from the 2-D rectified wavecal file to the flatfielded science extension header. This is the final step performed on the science data prior to 2-D rectification or 1-D extraction of the science data in the pipeline.

Either or both the wavecal file and science file can contain multiple exposures, and the image can drift across the detector over time due to such things as, thermal effects, so it is necessary to select the most appropriate wavecal exposure for each science exposure. Currently, the wavecal exposure nearest in time to a given science exposure is the one selected. Future enhancements may include interpolation of the appropriate shift.

21.4.3 Spectral Extraction

BACKCORR

In general, the background or sky is not aligned with the detector pixels. To accommodate this misalignment, the definition of the background extraction apertures includes not only a length and offset (center-to-center) but also a linear tilt to assist in properly subtracting the background. This tilt is taken into account when calculating the average background in the background extraction boxes.

DISPCORR

For MAMA data, offsets of the projection of the spectrum onto the detector in both the spectral and spatial directions are deliberately introduced by offsetting the Mode Select Mechanism (grating wheel) tilts. This is done approximately monthly to assure a more uniform charge extraction from the microchannel plate over time. For MAMA observations, these induced offsets are removed using coefficients in the mofftab table.

The disptab table of dispersion contains coefficients for fits to the following dispersion solution:

where

• is the wavelength in Angstroms.
• is the detector AXIS1 position.
• is the spectral order.
• are the dispersion coefficients. For each pixel in the AXIS1 direction, a wavelength is calculated. First, any modification to the dispersion coefficients due to spectrum offsets must be made. Table 21.4 lists the possible offsets and the appropriate corrections. For each integer value of s in the AXIS1 direction, a wavelength is solved for iteratively using the Newton-Raphson method.

  Modifications to the Dispersion Coefficients Caused by Offsets

  Correction	Ref Table	Algorithm	Definitions
  Incidence Angle	INANGTAB		dispersion coefficients incidence angle coefficients aperture offsets in the axis 1 direction calculated as difference of relative aperture centers (arcsec)
  MAMA Offsets	MOFFTAB		dispersion coefficients MAMA offset coefficients MAMA offset (MOFFSET1) (pixels) MAMA offset (MOFFSET2) (pixels)
  Echelle Spectrum Tilt	XTRACTAB		dispersion coefficients axis 2 offset from nominal A2CENTER during spectrum locate process (pixels) spectrum tilt angel

If FLUXCORR is PERFORM, the raw counts are corrected to (erg cm-2 sec-1 Å-1) using the reference files PHOTTAB and APERTAB. Execution of this calibration step requires that wavelengths have been assigned. Corrections for vignetting and echelle blaze are handled within the PHOTTAB reference files. The conversion to absolute flux is calculated as:

where:

• is the calibrated flux at a particular wavelength. This quantity is also multiplied by the ATODGAIN if the data were obtained with the CCD.
• is Planck's constant.
• is the speed of light.
• is the area of the unobstructed HST primary mirror (45238.93416 cm2).
• is the throughput of the STIS instrument configuration at a particular wavelength when a clear full aperture is in place.
• is a particular wavelength.
• is the dispersion (Å/pixel) at a particular wavelength.
• is the net count rate at a particular wavelength.
• is the aperture throughput at a particular wavelength.

where:

• is the heliocentric wavelength.
• is a particular wavelength.
• is the component of the velocity of the earth in the direction of the target.
• is the speed of light. The derivatives of low-precision formulae for the sun's coordinates described in the Astronomical Almanac are used to calculate the velocity vector of the earth in the equatorial coordinate system of the epoch J2000. The algorithm does not include earth-moon motion, sun-barycenter motion, nor light-time correction from the earth to the sun. This value for the Earth's velocity should be accurate to ~0.025 km/sec during the lifetime of STIS. (Note that the uncertainty of 0.025 km/s is much less than the ~2.6 km/s resolution obtained with the STIS high dispersion echelle gratings.) The value of heliocentric velocity, , is written to the trailer file.

X1DCORR

Locate the Spectrum

where the variables are as described in Table 21.2 sptrctab also contains the description of the distorted shape of the spectrum. The shape is stored as an array consisting of pixel offsets (in the AXIS2 direction) relative to the nominal center of the spectrum (A2CENTER). This spectrum trace is used to find, and eventually to extract, the 1-D spectrum.

The location of the spectrum is improved by searching in the vicinity of the nominal location of the spectrum by performing a cross-correlation between the distortion vector and the input spectrum image. The search extends for pixels around the nominal center, where n is read from the MAXSEARCH column in the xtracttab table. At each AXIS2 position in the search range (which differs from the nominal center by an integer number of pixels) a sum of the counts along the spectrum shape is formed. This sum is created by adding the value of one pixel's worth of data at each of the AXIS1 pixel positions. The pixel extracted in the AXIS2 direction is centered on the spectrum position (A2CENTER + pixel offset) and may include fractional contributions from two pixels. Quadratic refinement is used to locate the spectrum to a fraction of a pixel.

The final A2CENTER becomes:

where CRSCROFF is the offset found during the cross correlation. If the cross correlation fails, the value of CRSCROFF is set to zero, a warning message is written to the output, and the A2CENTER calculated prior to the cross correlation attempt is used as the location of the spectrum. CRSCROFF is written to the output science header.

An alternate method for performing the cross correlation may be employed. In this case a 2-D template is created from the spectrum trace table. The cross correlation is carried out between the 2-D template and input image. Quadratic refinement is used as above to refine the position of the center of the spectrum to a fraction of a pixel.

Extract the 1-D Spectrum

The extraction of the spectrum allows for unweighted or optimal extraction. The extraction algorithm is selected based on the value of the reference table parameter XTRACALG. This flag has possible values of UNWEIGHTED and OPTIMAL. The value of XTRACALG is written to the header of the output spectrum data file. At present, calstis performs an unweighted extraction of the 1-D spectra; the optimal extraction algorithm has not yet been implemented. At the end of the 1-D extraction step, a spectrum of gross counts/second is produced.

X2DCORR

Mapping from an output pixel back into the input image makes use of the dispersion relation and one-dimensional trace table. The dispersion relation gives the pixel number as a function of wavelength and spectral order. The one-dimension trace is the displacement in the cross dispersion direction at each pixel in the dispersion direction. Both of these can vary along the slit, so the dispersion coefficients and the one-dimensional trace are linearly interpolated for each image line. Corrections are applied to account for image offset, binning, and subarray. The spectrum can be displaced from its nominal location on the detector for several reasons, including Mode Select Mechanism (MSM) uncertainty, deliberate offsets for distribution of charge extraction for MAMA data, and the aperture location relative to a reference aperture. These offsets are accounted for by modifying the coefficients of the dispersion relations and by adjusting the location of the one-dimensional trace. See also DISPCORR and FLUXCORR, for algorithmic details (the process of dispersion solution, spatial rectification and flux and wavelength calibration is similar for one-dimensional and two-dimensional spectral extracted data).