Reduction of Long-Slit Spectra with PyRAF

This tutorial will use observations from program GS-2007A-Q-76 (PI: C. Winge), longslit spectra of interacting galaxy pairs, specifically AM2306-721. The spectra were obtained at 3 orientations with a 1.0arcsec slit and the B600 grating centered at 485.0 nm. Note that the spectral exposures have rather large CCD binning factors (\(2\times4\)), in part because this program was executed in the poor weather queue where poor seeing is the norm. The original technical goal was to measure emission line ratios at several positions along the slit. A flux calibration standard was included in the observing plan, as were comparison arcs at each slit position.

For other tutorials, see the following links:

Retrieve & Organize Data

The first step is to retrieve the data from the Gemini Observatory Archive (see Archive Searches). You may search the GOA yourself, or instead just cut-and-paste the direct URL in your browser.

# longslit data of galaxy pairs:
https://archive.gemini.edu/searchform/cols=CTOWEQ/GS-2007A-Q-76/notengineering/GMOS-S/LS/20070623/Win#

After retrieving the science data, click the Load Associated Calibrations tab on the search results page and download the associated bias and flat-field exposures. Unpack all of them in a subdirectory of your working directory named /raw. Be sure to uncompress the files. See Retrieving Data for details.

Caution

After unpacking the science data and associated calibrations, make sure to remove all files that begin with ‘GS’. These are reduced bias frames provided by the observatory. Running the script with these files in the raw directory will result in an error.

Processing Preparation

Reference Files

The required MasterCals are:

  • Bias Residual

  • Flat-field (from the GCAL source)

  • Wavelength calibration (from CuAr comparison arcs)

  • Flux calibration (from the standard star LTT 9239)

All of them will be constructed in this tutorial.

Software

You must create an observing log database of the data in the ./raw subdirectory. Download: obslog.py to the ./raw subdirectory, and execute it from the unix prompt.

python obslog.py obsLog.sqlite3

See Creating an Observing Log for details.

Also retrieve the python file selection module, which includes template SQL statements for selecting files, and functions for specifying metadata on which to perform selections.

Place this module in your work directory; it is used by the reduction script (below). You can perform all of the processing steps for this tutorial by downloading the Longslit Tutorial python script.

Place the script in the work directory, and execute the reduction script either from the Unix prompt,

python gmos_ls_proc.py

or from an active PyRAF session:

import copy
from pyraf.iraf import gemini, gemtools, gmos
import fileSelect as fs
import gmos_ls_proc

You may find it useful to download the script to follow this tutorial in detail, and use it as the basis for reducing other imaging observations.

Caution

This script may use a large amount of available storage if Google Drive is connected to your desktop. Disable synching for the working directory before running the script or a hidden directory named ‘.tmp.driveupload’, which is in the same directory as ‘vm_transfer.’, can take up 200+ GB of space.

Caution

The reduction script includes steps that should be performed interactively for best results, but the interactive options have been disabled in the script in order not to interrupt the automated processing flow.

Building MasterCals

The next steps will create the necessary MasterCal reference files that are used to calibrate the science files. Files are selected by matching specific exposure metadata in the observing log database (see GMOS Processing Keywords). Within the PyRAF session, first create the Bias Residual MasterCal:

# Observing log database
dbFile='./raw/obsLog.sqlite3'

# From the work_directory:
# Create the query dictionary of essential parameter=value pairs.
# Select bias exposures within ~2 months of the target observations:
qd = {'Full':{'use_me':1,
       'Instrument':'GMOS-S','CcdBin':'2 4','RoI':'Full',
       'Disperser':'B600+_%','CentWave':485.0,'AperMask':'1.0arcsec',
       'Object':'AM2306-72%',
       'DateObs':'2007-06-05:2007-07-07'}
      }
# Make copy for the CenterSpec RoI:
qd['CenSp'] = copy.deepcopy(qd['Full'])
qd['CenSp'].update({'RoI':'CentSp','Object':'LTT9239'})

# Set the task parameters.
gemtools.gemextn.unlearn()    # Disarm a bug in gbias
gmos.gbias.unlearn()
biasFlags = {
    'logfile':'biasLog.txt','rawpath':'./raw/','fl_vardq':'yes','verbose':'no'
}
regions = ['Full','CenSp']
for r in regions:
    # The following SQL generates the list of full-frame files to process.
    SQL = fs.createQuery('bias', qd[r])
    biasFiles = fs.fileListQuery(dbFile, SQL, qd[r])

    # The str.join() funciton is needed to transform a python list into a
    # comma-separated string of file names that IRAF can understand.
    if len(biasFiles) > 1:
        gmos.gbias(','.join(str(x) for x in biasFiles), 'MCbias'+r,
                   **biasFlags)

Now create the Flat-field MasterCal files from the GCAL flat exposures.

# Set the task parameters.
qd['Full'].update({'DateObs':'*'})
qd['CenSp'].update({'DateObs':'*'})
gmos.gireduce.unlearn()
gmos.gsflat.unlearn()
# The response fitting should be done interactively.
flatFlags = {
    'fl_over':'yes','fl_trim':'yes','fl_bias':'yes','fl_dark':'no',
    'fl_fixpix':'no','fl_oversize':'no','fl_vardq':'yes','fl_fulldq':'yes',
    'rawpath':'./raw','fl_inter':'no','fl_detec':'yes',
    'function':'spline3','order':'13,11,28',
    'logfile':'gsflatLog.txt','verbose':'no'
    }
for r in regions:
    qr = qd[r]
    flatFiles = fs.fileListQuery(dbFile, fs.createQuery('gcalFlat', qr), qr)
    if len(flatFiles) > 0:
        gmos.gsflat (','.join(str(x) for x in flatFiles), 'MCflat'+r,
                bias='MCbias'+r, **flatFlags)

Basic Processing

The gsreduce task has nearly 70 parameters; the table below lists the defaults for the processing flag keywords—i.e., the keywords with logical values to indicate whether to perform an operation. For the most part you can use the default parameter values; exceptions are noted explicitly in the code blocks below.

gireduce Processing Flag Defaults

Flag

Default

Description

fl_bias

Yes

Subtract bias residual?

fl_cut

Yes

Cut MOS SCI extensions into slitlets?

fl_dark

No

Subtract scaled dark image?

fl_fixpix

Yes

Interpolate across chip gaps?

fl_flat

Yes

Apply flat-field correction?

fl_fulldq

No

Decompose DQ extensions into component bits?

fl_gmosaic

Yes

Mosaic the image extensions?

fl_gsappwave

Yes

Insert approximate wavelength WCS keywords into header?

fl_gscrrej

No

Single-frame CR rejection?

fl_inter

No

Fit overscan levels interactively?

fl_over

Yes

Perform overscan correction?

fl_oversize

Yes

Scale slit lengths by x1.05 when cutting?

fl_title

Yes

Put MOS objID into title for each extension?

fl_trim

Yes

Trim overscan region?

fl_vardq

No

Propagate VAR and DQ?

Turning now to the science reductions, we first set some task parameters:

# Set task parameters.
gmos.gsreduce.unlearn()
sciFlags = {
    'fl_over':'yes','fl_trim':'yes','fl_bias':'yes','fl_gscrrej':'no',
    'fl_dark':'no','fl_flat':'yes','fl_gmosaic':'yes','fl_fixpix':'no',
    'fl_gsappwave':'yes','fl_oversize':'no',
    'fl_vardq':'yes','fl_fulldq':'yes','rawpath':'./raw',
    'fl_inter':'no','logfile':'gsreduceLog.txt','verbose':'no'
}
arcFlags = copy.deepcopy(sciFlags)
arcFlags.update({'fl_flat':'no','fl_vardq':'no','fl_fulldq':'no'})
stdFlags = copy.deepcopy(sciFlags)
stdFlags.update({'fl_fixpix':'yes','fl_vardq':'no','fl_fulldq':'no'})

# Arc exposures
for r in regions:
    qr = qd[r]
    arcFiles = fs.fileListQuery(dbFile, fs.createQuery('arc', qr), qr)
    if len(arcFiles) > 0:
        gmos.gsreduce (','.join(str(x) for x in arcFiles), bias='MCbias'+r,
                  **arcFlags)

# Std star exposures
r = 'CenSp'
stdFiles = fs.fileListQuery(dbFile, fs.createQuery('std', qd[r]), qd[r])
if len(stdFiles) > 0:
    gmos.gsreduce (','.join(str(x) for x in stdFiles), bias='MCbias'+r,
              flatim='MCflat'+r, **stdFlags)

# Science exposures
r = 'Full'
sciFiles = fs.fileListQuery(dbFile, fs.createQuery('sciSpec', qd[r]), qd[r])
if len(sciFiles) > 0:
    gmos.gsreduce (','.join(str(x) for x in sciFiles), bias='MCbias'+r,
              flatim='MCflat'+r, **sciFlags)

Examine the output files to assess data quality, and adjust the processing as necessary.

Wavelength Calibration

Image rectification and wavelength linearization depend upon the wavelength calibration, using the arc lamp exposures taken immediately before each sequence of science and standard star exposures (see Wavelength Calibration). In this case, the default medium-resolution line list will work well. The fit to the dispersion relation should be performed interactively, but for expediency we will use a previously determined functional fit.

# Set task parameters
gmos.gswavelength.unlearn()
waveFlags = {
    'coordlist':'gmos$data/CuAr_GMOS.dat','fwidth':6,'nsum':50,
    'function':'chebyshev','order':5,
    'fl_inter':'no','logfile':'gswaveLog.txt','verbose':'no'
    }
# Must select specific wavecals to match science exposures.
prefix = 'gsS20070623S0'
for arc in ['071', '081', '091', '109']:
     gmos.gswavelength (prefix+arc, **waveFlags)

Advanced Processing

The targets in this program were observed in 3 slit orientations, and a few exposures were obtained at each position. This provides an opportunity to combine the sequential exposures at each position to remove cosmic rays, rather than rejecting CRs on single frames using the gsreduce.fl_gscrrej+ flag or running the gemcrspec task. The combined exposures for each target are then wavelength calibrated, and sky subtracted. First set the processing parameters.

# Set task parameters.
gemtools.gemcombine.unlearn()
sciCombFlags = {
    'combine':'average','reject':'ccdclip',
    'fl_vardq':'yes','fl_dqprop':'yes',
    'logfile':'gemcombineLog.txt.txt','verbose':'no'
}
stdCombFlags = copy.deepcopy(sciCombFlags)
stdCombFlags.update({'fl_vardq':'no','fl_dqprop':'no'})
gmos.gstransform.unlearn()
transFlags = {
    'fl_vardq':'yes','interptype':'linear','fl_flux':'yes',
    'logfile':'gstransLog.txt'
}
# The sky regions should be selected with care, using e.g. prows/pcols:
#   pcols ("tAM2306b.fits[SCI]", 1100, 2040, wy1=40, wy2=320)
gmos.gsskysub.unlearn()
skyFlags = {
    'fl_oversize':'no','fl_vardq':'yes','logfile':'gsskysubLog.txt'
}

Standard Star

Flux calibration is a necessary final step for this program’s science goals. The standard star LTT9239 was observed as a part of this program using the CenterSpec RoI and was processed (above) in parallel with the target spectra.

Now combine the standard star exposures, apply the wavelength calibration, and subtract sky. The gsskysub task will determine the night sky emission spectrum from a selected spatial region, and subtract it row-by-row from the spectral image. While gsskysub can be run interactively, the selection of a sky region from a longslit spatial profile can be determined easily by inspection using the pcols task. Finally, extract the 1-D spectrum from the 2-D spectrogram, using a large (3 arcsec) aperture to ensure that all of the signal is captured:

prefix = "gs"
qs = qd['CenSp']
stdFiles = fs.fileListQuery(dbFile, fs.createQuery('std', qs), qs)
gemtools.gemcombine (','.join(prefix+str(x) for x in stdFiles),
                     'LTT9239', **stdCombFlags)
gmos.gstransform ('LTT9239', wavtraname='gsS20070623S0109', **transFlags)
gmos.gsskysub ('tLTT9239', long_sample='20:70,190:230')

gmos.gsextract.unlearn()
extrFlags = {
    'apwidth':3.,'fl_inter':'no','find':'yes',
    'trace':'yes','tfunction':'chebyshev','torder':'6','tnsum':20,
    'background':'fit','bfunction':'chebyshev','border':2,
    'fl_vardq':'no','logfile':'gsextrLog.txt'
}
gmos.gsextract ("stLTT9239", **extrFlags)

Note the need above to explicitly propagate the DQ and VAR extensions from the input files. Now derive the sensitivity calibration.

gmos.gsstandard.unlearn()
sensFlags = {
    'fl_inter':'yes','starname':'l9239','caldir':'onedstds$ctionewcal/',
    'observatory':'Gemini-South','extinction':'onedstds$ctioextinct.dat',
    'function':'chebyshev','order':9,'verbose':'no','logfile':'gsstdLog.txt'
}
gmos.gsstandard ('estLTT9239', sfile='std.txt', sfunction='sens',
                 **sensFlags)

Using these parameters for the single standard, the sensitivity residuals will be about 0.02 mag.

Science Targets

With all the MasterCals in place, the Science targets may be fully calibrated. The sky subtraction could be performed during the course of spectral extraction, but for this program measuring the extended, weak emission flux of interest calls for a custom extraction procedure. One can choose an appropriate sky region by plotting the spatial profile, e.g.:

iraf.pcols('tAM2306b.fits[SCI]', 1100, 2040, wy1=40, wy2=320)

Expanding the scale shows a good, source-free region, as shown below. Be sure to plot the sky spectrum in the selected regions to ensure there is no signal from the target. To display the sky spectrum, enter the following command into an active PrRAF session.

iraf.prows('tAM2306b.fits[SCI]',wcs="world",row1=520,row2=700,wy1=0,wy2=700)
../_images/AM2306b_profile.png

Screen shots of the spatial profile (left) of the spectrum for AM2306b, and the resulting sky spectrum (right). Click image to enlarge.

Now combine the science exposures, apply the wavelength calibration, and subtract sky. We will use a dictionary to associate the science targets with the Arcs and the pre-determined sky regions.

sciTargets = {
    'AM2306-721_a':{'arc':'gsS20070623S0071','sky':'520:720'},
    'AM2306-72_b':{'arc':'gsS20070623S0081','sky':'670:760,920:1020'},
    'AM2306-721_c':{'arc':'gsS20070623S0091','sky':'170:380,920:1080'}
}
for targ,p in sciTargets.iteritems():
    qs = qd['Full']
    qs['Object'] = targ
    # Fix up the target name for the output file
    sciOut = targ.split('-')[0]+targ[-1]
    sciFiles = fs.fileListQuery(dbFile, fs.createQuery('sciSpec', qs), qs)
    gemtools.gemcombine (','.join(prefix+str(x) for x in sciFiles),
                         sciOut, **sciCombFlags)
    gmos.gstransform (sciOut, wavtraname=p['arc'], **transFlags)
    gmos.gsskysub ('t'+sciOut, long_sample=p['sky'], **skyFlags)

Flux Calibration

Now apply the flux calibration to the science targets and the standard star. Note that gscalibrate works on both 2-D spectrograms and 1-D extracted spectra.

gmos.gscalibrate.unlearn()
calibFlags = {
    'extinction':'onedstds$ctioextinct.dat','fl_ext':'yes','fl_scale':'no',
    'sfunction':'sens','fl_vardq':'yes','logfile':'gscalibrateLog.txt'
    }
gmos.gscalibrate('stAM2306*', **calibFlags)
calibFlags.update({'fl_vardq':'no'})
gmos.gscalibrate('estLTT9239', **calibFlags)

Spectrum Extraction

The final step is to extract the science spectra in much the same way as in Krabbe et al. (2014). That is, spectra are extracted over an interval along the slit, and divided into segments of 4 spatial rows each. Recall that the CCD binning in the spatial direction is 4, so the spatial scale is 0.288 arcsec/pixel), yielding an extraction aperture of \(1.00\times1.15\) arcsec. The target SED consists of moderate-level, spatially variable continuum emission, stellar absorption, and sparse emission lines from ionized gas. Thus gsextract is not well suited for extraction in the way described by Krabbe et al. (2014).

Instead, use the onedspec.sarith task, even though it will not automatically propagate the VAR or DQ array (it is possible to work around this limitation). Load the onedspec package and extract based on inspection of the bright [O_III] 5007 emission (which is redshifted to about 5150), e.g. for AM2306b:

# Set the number of spatial pixels over which to sum
onedspec.nsum=4
onedspec.sarith('cstAM2306b.fits[SCI]', 'copy', '', 'ecstAM2306b.ms',
                  apertures='222-346x4')

An extracted spectrum of one of the brighter apertrues (#11) is shown below. To display the spectrum, enter the following command into an active PyRAF session.

iraf.splot('ecstAM2306b.ms.fits[*,11]',xmin=3925,xmax=5275,ymin=-0.25e-16,ymax=6.5e-16)
../_images/AM2306b_spec.png

Screen shot of a portion of a 1-D spectrum for AM2306b, which is the sum of 4 spatial pixels. Click image to enlarge.

Possible improvements to the above process, which are left as an exercise, include the construction and use of a Static Bad-Pixel Mask MasterCal and a rejection of cosmic rays.