How can I combine STIS spectra?

Below we outline a couple of methods you can use to combine STIS spectra in IRAF/STSDAS.

Method A: Splicing the _x1d.fits files (recommended)

We recommend using the task splice - located in stsdas.hst_calib.ctools - as it takes into account the error (ERR) array as well as the data quality (DQ) array. Handling of the DQ array is important as it helps splice perform the combination properly and avoid bad/noisy data in the output file due to the large changes in throughput at the edges of the detector. In order to use splice correctly, one has to manually mark these edges in the DQ array as "bad pixels" and tell splice to ignore these by setting the SDQFLAGS parameter to the appropriate value (referring to Table 2.8 in Section 2.5.2 in the STIS Part of the Data Handbook for the data quality values). A flag value of 4 or 8 will work fine. The reason this does not work by default in splice at the moment is that the vignetted regions of the MAMA detectors change over time, and new bad pixel maps showing these regions as a function of time have not been produced. So at the moment, users need to modify the DQ arrays themselves. A simple way to do this is to edit the DQ array in the _flt.fits image, flagging the edges as bad pixels, then re-extracting the 1-D spectrum using the x1d task in the stis package. For example:

To set a 10-pixel border to a value of 4 ("bad detector pixel") in the 1024x1024 image:

imcalc o4qx04010_flt.fits[dq,1] o4qx04010_flt.fits[dq,1] \
"if x <= 10 || x >= (1024-10) || y <= 10 || y >= (1024-10) then 4 else im1"

disp o4qx04010_flt.fits[sci,1] 1
disp o4qx04010_flt.fits[dq,1] 2 zsc- zra- z1=-100 z2=1025

Once, you have followed the steps above, you can rerun x1d on the _flt.fits file and get a new x1d spectrum, which will have the updated DQ values. Then you can run splice, specifying that you want to ignore those vignetted regions that you just marked with DQ values when you do the splicing of the spectrum together, i.e., by making sure the SDQFLAGS parameter in the splice task is set to an appropriate value. 

splice obs1_x1d.fits,obs2_x1d.fits output_splice.fits

Running splice as mentioned above (rather than transforming individual _x1d fits tables into multispec format before combining them) has important advantages: it keeps the science data, error and DQ arrays intact allowing for easier error analysis, and it does not have a limitation on the number of orders or wavelengths to include, a problem with the multispec format due to the limit on the size of the fits header which requires fitting the wavelength scale with a function.

Method B: Adding the _flt.fits files

If your observations were taken very close in time, with no intervening change of grating or central wavelength, then the shifts in wavelength should in principle be small and you should be able to add the _flt.fits images. However, we recommend that you check the keywords SHIFTA1 and SHIFTA2 in the _flt.fits files to make sure this is the case. These keywords have the values in pixels of how much the position of the spectral image has shifted on the detector, as calculated from the contemporary wavecals taken. (The shift of the spectrum due to errors in tracking or re-acquiring the target is usually negligible.) If these shifts are small, you should be safe in combining the _flt.fits files (e.g., using tasks msarith or mscombine in the mstools package when dealing with MAMA spectra, or using ocrreject in the stis package for CCD spectra) to create a summed or averaged fits file. Depending on the science goals, the detector, and the nature of the target, shifts of up to 0.5 pixels may be acceptable. Better alignment will likely be needed for CCD spectra of sources with compact components or sharp velocity structure, namely to prevent cosmic ray rejection from discarding a significant amount of real flux. This is increasingly relevant as the signal to noise of the spectrum increases, since the rejection level is a function of the noise. Ocrreject can be set to flag which pixels are rejected in the DQI extension of the input images so that you can check for inappropriate systematic rejection along rows or columns. The cosmic ray rejection parameters can be relaxed to accommodate some misalignment.

A slight variation on this method is appropriate when dealing with CCD data which were offset from one another in the spatial direction by integer pixels by the observer (e.g., to enable one to eliminate any residual hot pixels in the summed product). In that case, one can run task sshift (in the stis package) to produce registered images suitable for combining.

Warning: When creating a combined fits file, do make sure the TEXPTIME keyword in the main header of the combined fits file (e.g., "summed.fits[0]" or "averaged.fits[0]") has the appropriate value, otherwise the flux scale will be incorrect in the final spectrum. After doing the addition, you can use the x1d task in the stis package to extract a spectrum from the summed fits file, which can be treated like a normal _x1d.fits file and can be run through tomultispec.

Method C: Adding the _x2d.fits or _sx2.fits files

If the difference between the SHIFTA1 and SHIFTA2 keywords in the individual _flt.fits files is too large, do not just add them. In this case, you can use the 2D-rectified images (i.e., the _x2d or _sx2 files) which have been photometrically calibrated. You can first average these images. Then to obtain a spectrum, you would need to determine which rows contain the spectrum of your object and sum those rows with a task like imcalc.

Keep in mind that the results of the task x2d are in a different format and units than those from x1d. The x2d images are calibrated in flux for diffuse objects (i.e., using surface brightness units), so if the object observed is a point-source, you would need to sum over the default extraction height and multiply by the DIFF2PT keyword to convert to point-source flux. DIFF2PT uses wavelength-averaged flux corrections for slit width and extraction height, so a spectrum produced this way will not be calibrated as accurately as an x1d spectrum. See Sections 3.4.12 and 5.4 of the STIS Part of Data Handbook.