|LONG SLIT CCD AND MAMA CALIBRATION AND REDUCTION PROCEDURES FOR NGC1976 USING IRAF
NASA / Ames Research Center, University of Northern Colorado, Greeley, CO
Long slit spectra from Hubble Space Telescope’s STIS instruments require
calibrations before they can be interpreted in terms of physical phenomena.
This report elucidates the detailed process in which data is taken from its most
basic form to an accurate, refined representation of observable physical
characteristics. We use these data to derive electron temperature, electron
density, and chemical abundances of gaseous nebulae, in particular planetary
nebulae and H II regions.
From the raw science images extracted from the Hubble Space Telescope (HST) archive, concise and workable data must be fabricated. Accomplishing this objective requires that certain specific manipulations be done to the HST science data. The data must be coherently grouped, or co-added together, radiation events and phenomena such as cosmic rays must be detected and eradicated from the good data, and the image must be put in proper flux units eventuating in a full 2-D rectified image. This is done using the Image Reduction and Analysis Facility (IRAF). This report covers the steps to process spectroscopic data taken with the Space Telescope Imaging Spectrograph (STIS).
The 2-D image representing a spectrum from which specific lines can be extracted is the desired workable form of data of which we speak. At this point the researcher is presented with a myriad of options with which a direction to take can be dictated by the specific researcher's interests and goals. Our interests included finding the ratio of various emission lines in order to both test the accuracy of our method and data, and to construct extinction and temperature maps for several nebulae. This report describes work done on the H II region known as the Orion Nebula; it has also been applied to the planetary nebula NGC 7009.
2.1 IMPORTANT NOTE
Before commencing with reduction steps of new raw data the following line indicating where the reference files are located must be added to the user's .cshrc file:
setenv oref "directory where reference files are located"
i.e. setenv oref /home/svoboda/iraf/STIS/V5/
Notice the slash after the directory the reference files are in. In this case V5. This line MUST end in a slash for IRAF to recognize the syntax. After saving the modified .cshrc file the user must source the file by entering the command
into the UNIX command line. To be safe it is also a good idea to logout of UNIX and then log back in again.
This line in the .cshrc file MUST be changed every time images from a different visit or different directory are reduced. The user must remember to do this, as the computer will not give a warning when performing the calibrations. It will simply use the wrong reference files.
2.2 RETRIEVAL OF RAW SCIENCE IMAGES
The obvious first step in the reduction of data is to receive it from HST. This is done by downloading it from the HST archive. The data that comes from the archive is not necessarily in purely raw form. Various on-the-fly calibration steps are performed on the data at the Space Telescope Science Institute (STScI) soon after these data are received from HST at STScI. The data from the pipeline has already been trimmed. Average bias and overscan have been removed, and the dark current has been removed if necessary. In more concise terms, the images have been flat fielded and given the extension _flt. (Massey 1997) One should realize that in the initial pipeline processing that STScI uses the "best calibration files" available at that time. Better calibration data may well follow later than the subject scientific data set. Thus it is often prudent to wait until these "better" calibration files have also entered the archive before retrieving the data for "final" analysis.
It is worth noting that there are two types of On-the-Fly processing done by STScI. On-the-Fly calibration (OTFC) is what has been used in the past, but a new processing method known as on-the-fly-reprocessing (OTFR) is now being implemented. The advantage of the new method is that it regenerates the raw files before calibration, thus header values can be assigned to the image more accurately. (Micol on behalf of ST-ECF Archive Team 2001) Our Visit 2, Visit 4, and Visit 5 data were retrieved using OTFC, while our Visit 51 and Visit 52 data were retrieved using the new OTFR method.
For STIS, depending on the scientist's observing program instructions, wave calibration may (automatic WAVECALS taken adjacent to science exposures) or may not have been performed. In the latter case, separate WAVECALS with their own exposure number were scripted. The rationale for doing it this way is that the WAVECALS can often be relegated to the occultation period and thus may result in substantial savings of usable observing time for exposures on source.
2.3 WAVE CALIBRATION
The Quick Look packet for a given slit and visit can be utilized to determine whether wave calibration has transpired. The exposures taken are listed on the "Observation List" page. If this list includes exposures having either the target name waveline or wavehitm then wave calibration must be performed. If these exposures are not present, then the wave calibration is omitted because it was done automatically.
Wave calibration is used to accurately affix a wavelength scale to the long slit spectra. The function to perform wave calibration, among other things is called calstis, which is located within the stis package, located within the hst_calib package, located within the stsdas package. This can be more succinctly represented by the notation (stsdas.hst_calib.stis). It is this notation that will be used in this manual here to forth.
Before calstis is performed, the function chcalpar (stsdas.hst_calib.ctools) must be used to change the calibration parameter files of the images. A detailed description of the chcalpar function is located in a latter section of this report entitled "CHCALPAR''. For the sake of brevity it will just be said here that the wavecal parameter must be changed from OMIT to PERFORM.
The wavecal is to be performed on files having the rootnames listed under the Quick Look's Observation List followed by the extension _raw.fits. The wavecal files themselves will also be of this type with rootnames listed in the Observation List. The parameter file for the calstis function is as follows:
input = o4fo05090_raw.fits Input STIS _raw FITS files
(wavecal= o4fo05080_raw.fits) Input raw wavecal image files
(outroot= o4fo05090_raw.w.fits) Root for output file names
(savetmp= no) Save temporary files?
(verbose= yes) Print verbose time stamps?
(Version= 13Jan2000) calstis version
(mode = al)
The wavecal file should be the file in the Observation List with the appropriate target name, and matching the central wavelength of the input image. If there is more than one file fitting these criteria, then the file closest in time or with the longest exposure time is preferred. Note there would be a way to use more than one WAVECAL file with a science exposure, but we did not need to do this. The outroot file should be identical to the input file with the addition of a .w after the _raw, but before the .fits suffix.
After the proper information has been entered into the parameter file the command :go is entered in order to execute the calstis wave calibration function. This is the method by which any function is activated from its parameter file. It will be assumed throughout the remainder of this manual that functions will be activated in this manner. When the wave calibration is completed the extension _flt will be appended to the file name yielding the file:
The .w is added to distinguish a wave calibrated image from an image with no calibration.
2.4 CO-ADDING IMAGES
In order to increase the signal-to-noise ratio, we co-add images. Once the wave calibration of the individual exposures is completed, the exposures must be grouped together according to central bandpass wavelength and slit location. The exposures within these groups are then co-added using the msjoin (stsdas.toolbox.imgtools.mstools) function.
This is simply accomplished by using the Observation List from the quick look packet. From this list the rootnames of the images with corresponding central wavelengths and slits can be determined. The images that we wish to join are the wave calibrated images created in the previous step. These images are entered, separated by a comma, in the inimg field in the parameter file for msjoin. If the wave calibrations took place automatically, then the files to be msjoined will have the rootname along with an _flt extension:
The name of the output image needs to be entered into the parameter file as well under the outimg field. The name of the output image should be as follows:
(central _)slit(slit number)_cmb.fits i.e. 4451slit4_cmb.fits
An example of the msjoin parameter file:
inimg = o4fo05090_raw.w_flt.fits,
o4fo050a0_raw.w_flt.fits Images to join
outimg = 4451slit4_cmb.fits Output image name
(din_lis= ) Internal script use only
(Version= Aug97) Date of installation
(mode = al)
The suffix _cmb acknowledges the fact that this image is made up of combined exposures. The Primary Header Unit for this image will be identical to the Primary Header Unit of the first image entered into the
2.5 COSMIC RAY/BAD PIXEL REJECTION
Now that the images have been combined, the anomalous pixels must be removed. During any given exposure cosmic rays will randomly hit the CCDs. The longer the exposure times, the more cosmic ray hits on the CCD. These cosmic rays can severely corrupt that data as they can be up to 100,000 times deviant from the mean pixel value in an image. The IRAF function used to search out and reject cosmic rays is ocrreject (stsdas.hst_calib.stis). The function's parameter file is included below, as it will be referred to throughout the explanation.
input = 4451slit4_cmb.fits Input STIS file name(s)
output = 4451slit4_crj.fits Name of output file
(all = yes) Combine all input files into a
(crrejta= ) Reference file name
(scalens= ) Scale factor applied to noise
(initgue= med) Initial value estimate scheme
(skysub = none) How to compute the sky (none|mode)
(crsigma= 6,4,4,3,2,2) Rejections levels in each iteration
(crradiu= 1.) CR expansion radius in pixels
(crthres= 1.) Rejection propagation threshold
(badinpd= INDEF) Data quality flag bits to reject
(crmask = yes) Flag CR in input DQ images
(verbose= yes) Print additional info?
(Version= 13Jan2000) calstis version
(olist = )
(mode = al)
The ocrreject function applies a cosmic ray rejection algorithm to iteratively remove bad pixels. Because of the way the STIS Pipeline group has chosen to assemble the images prior to download, it is recommended that the cosmic ray rejection be performed after the data quality file has been initialized and after overscan bias has been subtracted, but before the bias frame and the dark have been subtracted and before flatfielding of the image.
Multiple images (or use of a wildcard) can be entered in the input field of the parameter file. The CRREJTAB located in the primary header file contains information about how the cosmic ray rejection will be performed. If the parameter file's fields are left undefined then the cosmic ray rejection will be done in the manner the CRREJTAB dictates. Any modification to the parameter file overrides what is in the header file.
Notice that the _cmb in the input field file name has been replaced by crj in the output file name. This is simply the extension used by convention to represent a cosmic ray cleaned image.
An entry in the crrejtab field indicates what table to read the rejection parameters from. If no entry is made here the rejection parameters are read from the table indicated by the CRREJTAB keyword located in the header file. (IRAF online help manual)
Also worth noting are the settings of the skysub field which indicate whether or not sky will be subtracted before the cosmic ray (CR) rejection. The field crsigmas defines both the number of iterations, and the allowable pixel deviation (defined as a multiple of sigma) for each iteration. The crradius defines the number of pixels surrounding a cosmic ray hit to be eliminated.
Cosmic ray rejection, while essential in CCD reduction, is unnecessary when reducing MAMA images. The MAMAs are made out of material that is not susceptible to the cosmic ray interference that CCDs are. If one is reducing MAMA spectra it is at this point that the reduction path departs from the method used for CCDs.
2.6 MAMA REDUCTION
At the point where cosmic ray rejection is performed on the CCDs one must instead use the chcalpar function (described in a latter section) on the MAMA images and change the rptcorr calibration parameter of the _cmb images to PERFORM. The dispcorr, helcorr, fluxcorr, x2dcorr, x1dcorr, and backcorr should already be set to PERFORM, and they must be changed if they are not. With these calibration parameters in place the calstis should be performed which will produce three files:
The three resulting files and steps leading up to the production of these files seem for all intensive purposes to be unnecessary. While IRAF experts have approved the reduction recipe that includes these steps, they do not seem at all functional.
At this point another chcalpar is performed on the _cmb files changing dispcorr, helcorr, fluxcorr, x2dcorr, x1dcorr, and backcorr back to OMIT. The rptcorr parameter should be left to PERFORM. Another calstis is executed on these files which results in the following file:
From this point in the reduction process forward, the MAMA method is again identical to that used for CCDs. The only difference being that the chcalpar and calstis described subsequently is performed on the _cmb_sfl.fits files as opposed to the _crj.fits files. Thus, the final 2d rectified output file for MAMA images is for example:
The chcalpar (stsdas.hst_calib.ctools) function is used to modify an image's calibration parameters. The calibrations are then carried out by the calstis function. The chcalpar function can be used with four different types of images, STIS CCD Imaging, STIS CCD Spectroscopy, STIS MAMA Imaging, and STIS MAMA Spectroscopy which correspond to keywords ckwstis1, ckwstis2, ckwstis3, and ckwstis4 respectively. These keywords can be entered in at the time the function is called though this was not necessary in our case as the image type for both our CCD and MAMA reduction was recognized from the header automatically by chcalpar. The following is the parameter file for the chcalpar function:
images = 4451slit4_crj.fits List of images to modify
(templat= ) Image to read header from
(keyword= ) Pset to use if not reading from an
(add = yes) Add keywords if not present in header?
(verbose= yes) Print out files as they are modified?
(Version= 21Jul97) Date of Installation
(mode = al)
The parameter file is straightforward and rarely needs to be modified. Because of this it is often easier to start chcalpar from the command line which looks like the following example:
The following is the image's parameter list, which is the screen that follows the command line entry:
(dqicorr= complete) initialize data quality?
(blevcor= complete) bias level correction?
(biascor= complete) bias correction?
(darkcor= complete) dark correction?
(flatcor= complete) flat field correction?
(crcorr = complete) combine CR exposures?
(expscor= omit) basic 2-D on CR flagged data?
(rptcorr= omit) add individual repeat obs?
(wavecor= complete) use wavecal?
(dispcor= perform) use dispersion solution?
(helcorr= perform) heliocentric correction?
(fluxcor= perform) convert to absolute flux?
(x2dcorr= perform) rectify 2-D spectral image?
(x1dcorr= perform) 1-D spectral extraction?
(backcor= perform) subtract background?
(bpixtab= oref$h1v11475o_bpx.fits) bad pixel table
(ccdtab = oref$k2g1502eo_ccd.fits) table of CCD parameters
(biasfil= oref$jbo0814go_bia.fits) bias image reference file
(wbiafil= oref$k5h1101jo_bia.fits) wavecal bias image reference file
(darkfil= oref$jbo0814ho_drk.fits) dark reference file
(pfltfil= oref$k2910263o_pfl.fits) flat field reference file
(dfltfil= N/A) delta flat reference file
(lfltfil= ) low order flat reference file
(apertab= oref$jba1538mo_apt.fits) aperture throughput table
(apdesta= oref$j8j1559jo_apd.fits) aperture description table
(sptrcta= oref$l2j0137so_1dt.fits) 1-D spectrum trace table
(disptab= oref$l2j0137to_dsp.fits) dispersion coefficients table
(inangta= oref$h5s11397o_iac.fits) incidence angle correction table
(lamptab= oref$l421050oo_lmp.fits) template lamp spectrum
(wcptab = oref$k291025go_wcp.fits) wavecal parameters table
(sdctab = oref$h4s1350do_sdc.fits) 2-D spectral extraction parameters
(phottab= oref$k9f1452qo_pht.fits) photometry calibration table
(pctab = oref$k6q16101o_pct.fits) photometric correction table
(crrejta= oref$j3m1403io_crr.fits) parameters for cosmic ray rejection
(xtracta= oref$i1c1615po_1dx.fits) 1-D spectral extraction table
(instrum= stis) Instrument represented by this pset
(detecto= ccd) Detector represented by this pset
(obstype= spectroscopic) Obstype represented by this pset
(Version= 23Feb2000) Date of Installation
(mode = al)
The parameters from bpixtab down define the reference files to be used by calstis during the calibration process. The parameters before bpixtab define what calibrations have been performed along with what calibration will be performed in the next calstis. When a function performs a calibration step it automatically changes that parameter of the header keyword to COMPLETE.
The calibrations dqicorr, blevcor, biascor, darkcor, and flatcor come complete when the images are downloaded from HST. Wavecor is completed in the first step of our reduction process, and crcorr is completed in the previous step of the reduction process. It is at this point, before running our final calstis on the image, we need to change the parameters dispcor, helcorr, fluxcor, x2dcorr, x1dcorr, and backcor from OMIT to PERFORM. It can be seen from the above parameter file what part of the calibration each of these fields represent.
Once these adjustments have been made to the parameter file, activate the function. The screen will redisplay itself with the new entries. Type :q and the screen will redisplay itself only listing values for the parameters that were modified. Typing :q again will give a prompt to “Accept the current parameters?” (|no|yes|abort). Choose yes and the images will be modified.
Performing the calstis is the last step in the reduction process. It performs six calibrations that, all together, result in the 2d rectified image that is the goal of the reduction part of the data analysis.
Once the chcalpar has been used to modify the calibration parameters, using the calstis is quite simple. The parameter file for the calstis function is as follows:
input = 2739slit1_cmb_sfl.fits Input STIS _raw FITS
(wavecal= ) Input raw wavecal image files
(outroot= ) Root for output file names
(savetmp= no) Save temporary files?
(verbose= yes) Print verbose time stamps?
(Version= 13Jan2000) calstis version
(mode = al)
It is straightforward and rarely needs to be modified. The input file is the _crj file produced from the cosmic ray rejection. The six calibrations that will be executed are described in more detail below as found in the IRAF calstis help manual (IRAF online help manual):
If set to PERFORM, the output 1-D extracted spectra and 2-D
rectified images will be corrected for the earth's heliocentric
velocity. Note that for medium and high dispersion, the
Doppler shift due to the spacecraft velocity around the
earth will already have been corrected by the on-board software.
If set to PERFORM, the output 1-D extracted spectra and 2-D
rectified images will be converted to absolute flux units.
Rectify 2-D spectral image. If set to PERFORM the input
image will be 2-D rectified (x2d).
Perform 1-D spectral extraction. If set to PERFORM a
1-D spectral extraction will be performed (x1d).
Apply dispersion solutions during 1-D extraction. (This
step is always done for 2-D rectification.)
Subtracts Background from the image.
This results is two files called:
3 LINE EXTRACTION
3.1 SECTION IDENTIFICATION
The ability to single out an emission line from a 2d rectified image for further investigation is paramount to our research motives. Finding the ratio of these lines can give a measure of the accuracy of the data, a measure of extinction levels in a region of the nebula, and a measure of electron temperature. In order to look at emission lines free of bad pixels and continuum pollution, five specific section of the image must be identified per line of interest. The five section are made up of 1) continuum on the blue side of the line 2) continuum on the red side of the line 3) continuum immediately adjacent to the line on the blue side (blue hugger) 4) continuum immediately adjacent to the line on the red side (red hugger) and 5) the emission line itself. These sections should all be free of other emission lines, and any other non-physical or physical interference. The image section should be recorded for future reference. It is best to do this in the IRAF syntax (shown below) used to define a section of a larger image.
In order to simplify the math it is most prudent to choose sections of a certain size. The height (delta Y) of each section should be equal. The width (delta X) of the line plus the width of the immediately adjacent huggers on either side should sum to be an even number. The delta X of each of the continuum sections should be half the width of the total line plus huggers. In our experimentation we have found that a delta X of 30 works well for the huggers and line, and a delta X of 15 for each of the continuum sections.
3.2 SECTIONAL CLEANING
Before any cleaning transpires a copy of the full 2d rectified image should be made calling it i.e. 4451slit4_crj_x2d_fix.fits. A mask should be made to fix the fiducial bars in the new image. Creating a text file with the following as an example can do this:
1 (x value spanning image) 366 377
1 (x value spanning image) 825 842
In this case the fiducial bars spanned from 366 to 377 and 825 to 842 in the Y dimension. The fiducial bars can then be eliminated using this mask in conjunction with fixpix whose functionality is elucidated later in this article.
While IRAF's automated cosmic ray rejection task eliminates a high percentage of an image’s bad pixels, it does not eradicate them all. What is left as part of the data can still corrupt it considerably. In order to be able to run the proverbial fine-tooth comb through our data eliminating essentially all bad pixels we created a program called Pixhunter. Pixhunter works by computing the mean pixel value and standard deviation of
a given section. After computing these values it again goes through the image and outputs to a file the x and y coordinates of any pixel whose value is further away from the mean than plus or minus 5 times the standard deviation. This output is called a bad pixel mask and can be used with IRAF's fixpix (proto) which will linearly fit the bad pixels to the data. Fixpix's functionality is straightforward and an example of its parameter file follows. The file entitled 4451slit1.mask is the output of the Pixhunter program.
images = 4451slit1_crj_x2d.fits[242:256,90:1101] List of images to be
masks = 4451slit1.mask List of bad pixel masks
(linterp= INDEF) Mask values for line interpolation
(cinterp= INDEF) Mask values for column interpolation
(verbose= yes) Verbose output?
(pixels = no) List pixels?
(mode = ql)
Though the Pixhunter works well for the continuum sections and huggers, its algorithm for finding the bad pixels is incompatible with emission lines. The continuum is reasonably flat in the flux dimension whereas an emission line's flux can vary greatly in both the wavelength and spatial dimensions. We have found that in the condition of an emission line, splot (noao.onedspec) can be utilized for cleaning.
A single column of the line must be plotted. The columns can then be scrolled through using the ( and ) keys. For each individual column splot can fit a function to the line. Splot can use four different methods of function fitting, but we have found that the Legendre polynomial fit with polynomial degree equal to 150 works well. This fit is activated with the t key followed by the c key for cleaning. Splot will then fit a function to the emission line and flag all pixels that are a multiple of sigma or more away from the function fit. The q key then fits all of the flagged pixels to the function. If the resulting image is desirable it can be saved back into the original 2d rectified image with the i command.
It must be noted that splot will always clean the pixels in an entire column even if the entire column is not being displayed. In the case where there are stellar-like phenomena in an emission line the image must be broken up using imcopy and reassembled using imjoin. Only the sections without drastic peaks such as caused by the P159-350 proplyd in our Orion work can be cleaned by this method.
3.3 CONTINUUM SUBTRACTION
Now that the line and sections of continuum are free of deviant pixels, the continuum must be subtracted from the line. Begin by adding the two 15 pixel wide sections of continuum together using the imarith function. The 30 wide section of line and huggers and the summed continuum must now be collapsed in the x dimension. This is done using IRAF’s blkavg function. The parameter file to do such is as follows:
input = 4451slit4.cont.fits Input images
output = 4451slit4.cont.sum.fits Output block averaged images
(option = sum) Type of operation
b1 = 15 blocking factor in dimension 1 (x or
b2 = 1 blocking factor in dimension 2 (y or line)
b3 = 1 blocking factor in dimension 3 (z or band)
b4 = 1 blocking factor in dimension 4
b5 = 1 blocking factor in dimension 5
b6 = 1 blocking factor in dimension 6
b7 = 1 blocking factor in dimension 7
(mode = ql)
The b1 value should be equal to the delta X of the section. Notice that we are not actually averaging the section, but summing it. Once the line/huggers section and the continuum-summed section have been collapsed to one dimension the continuum can be subtracted from the line using imarith.
The most accurate calculation for the theoretical intensity ratio for [O III] I(5007)/I(4959) is 2.984 (Storey & Zeippen 2000). Due to the fact that both of these emission lines are the results of transitions from the same energy level the ratio is independent of electron temperature and density. Using data that was reduced as per the described method we found excellent agreement with this theoretical value. For Orion Visit 5 slit 5 for example, we found the flux ratio of 5007 to 4959 to be 3.058, which is slightly higher than the theoretical value. This is expected because interstellar extinction causes the 4959 line with shorter wavelength than the 5007 line to suffer relatively more attenuation due to extinction. From the reduced data we are also in the process of deriving electron temperatures, which we have so far found to be in good agreement with both theoretical values and past observations.
I would like to thank Bob Rubin, Naman Bhatt, and Erik Shimshock for all the patience they have shown, and for all the guidance they have given me this summer.
Massey, P. 1997, IRAF Online Archive,
Micol, A. on behalf of the ST-ECF Archive Team 2001, ST-ECF Newsletter, 17, 29
Storey, P. J., & Zeippen, C. J. 2000, MNRAS, 312, 813