Reduction of Multi-Object Spectra with PyRAF

This tutorial will reduce selected observations from program GN-2011B-C-3: MOS spectra of the H_II region populations in NGC185, and two distinct fields in M81 (NGC3031), as well as the Galactic planetary nebula PN_G235.3-03.9. The observing program was executed on the night of 2012 Jan 24-25, and included g-band imaging of the fields of interest (partly for acquisition of the targets), and MOS spectroscopy using gratings:

  • B600 centered at 520 and 525 nm.

  • R400 centered at 525, 740 and 745 nm

The original science goal was to determine chemical abundances from emission line ratios in several H_II regions within the galaxy. Flux calibration standards were included in the observing plan, as were comparison arcs, GCAL and twilight flats.

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, but be sure to restrict the search to limit the number of science and calibration files to the relevant for reducing these spectra:1

  • Mode: Spectroscopy

  • Adaptive Optics: Not AO

  • Nod&Shuffle: Classic

  • Binning: \(2\times2\)

Footnotes

1

You will also need to download the M81 masks: GN2011BC003-01.fits and GN2011BC003-02.fits, and place them in your work directory.

You may alternatively cut-and-paste the direct URL in your browser.

# M81 H_II region MOS spectroscopic data:
https://archive.gemini.edu/searchform/Classic/cols=CTOWMDBELQ/NOTAO/notengineering/2x2/GN-2011B-C-3/20120101-20120301/GMOS-N/NotFail/spectroscopy#

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.

It is essential that you review the observing log to decide how to select exposures for processing. It is advisable to display the files of interest to eliminate any bad exposures before they corrupt the reduction processing.

Observing Configurations

The table below shows the set of observing configurations that were used for the targets in this program. Since the number of configurations is large, this tutorial will be restricted (filled squares) to reducing MOS exposures of the first M81 field and associated calibrations, mostly for gratings B600/520 and R400/740. Note that the standard star exposures were observed with the Central Spectrum RoI, but most other exposures were obtained with the Full RoI.

Observing Configurations Used for GN-2011B-C-3

Config.

NGC 185

G191-B2B

‘PN G235.3’

‘M81-field1’

‘M81-field2’

‘HZ44’

B600/420

\(\blacksquare\)

\(\blacksquare\)

B600/520

\(\diamond\)

\(\blacksquare\)

\(\diamond\)

\(\blacksquare\)

\(\diamond\)

\(\blacksquare\)

B600/525

\(\diamond\)

\(\blacksquare\)

\(\diamond\)

\(\blacksquare\)

B600/620

\(\blacksquare\)

R400/420

\(\blacksquare\)

R400/525

\(\diamond\)

\(\blacksquare\)

\(\diamond\)

R400/620

\(\blacksquare\)

\(\blacksquare\)

R400/740

\(\diamond\)

\(\blacksquare\)

\(\diamond\)

\(\blacksquare\)

\(\diamond\)

R400/745

\(\diamond\)

R400/900

\(\diamond\)

1.0arcsec

\(\diamond\)

\(\blacksquare\)

\(\diamond\)

\(\blacksquare\)

MOS

\(\blacksquare\)

\(\diamond\)

Full

\(\diamond\)

\(\blacksquare\)

\(\diamond\)

CenterSp

\(\blacksquare\)

\(\diamond\)

\(\blacksquare\)

Note

All of the standard star settings except R400/900 will be reduced because the wavelength range covered by the MOS slits for a given grating and central wavelength setting is much greater than that for a single long-slit configuration. The R400/900 setting is not useful because of second-order light \(\gt 910\) nm from the very hot standard star SED.

Processing Preparation

Reference Files

The required MasterCals are:

  • Bias Residuals

  • Flat-fields

  • Gradient images (unnormalized flat-fields used to determine the locations of the slitlets on the detector format)

  • Wavelength calibrations (CuAr comparison arcs)

  • Flux calibration

All of them will be constructed in this tutorial.

Software

Before starting a new PyRAF session, download the monochromatic mags files for G191B2B and HZ44. To obtain copies of these files, enter the following command in the same directory as anaconda2/:

find anaconda2/ -name *hz44*
find anaconda2/ -name *g191b2b*

Copy one of the hz44.dat and g191b2b.dat files into the working directory.

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 MOS tutorial python script.

From an active PyRAF session:

import copy
from pyraf import iraf
from pyraf.iraf import gemini, gemtools, gmos
import fileSelect as fs
import gmos_mos_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 turned off 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:
qdf = {'use_me':1,
       'Instrument':'GMOS-N','CcdBin':'2 2','RoI':'Full',
       'Disperser':'B600+_%','CentWave':520.0,'AperMask':'GN2011BC003-01',
       'Object':'M81-field1',
       'DateObs':'2012-01-16:2012-02-12'
       }
# Need another copy for the CenterSpec RoI:
qdc = copy.deepcopy(qdf)
qdc.update({'RoI':'CentSp','AperMask':'1.0arcsec','Object':'G191B2B'})

# Now prepare to call the IRAF processing task.
# Use primarily the default task parameters.
gemtools.gemextn.unlearn()    # Disarm a bug in gbias
gmos.gbias.unlearn()
biasFlags = {
    'logfile':'biasLog.txt','rawpath':'./raw/','fl_vardq':'yes',
    'verbose':'no'
}
# The following SQL generates the list of 60 full-frame files to process.
SQL = fs.createQuery('bias', qdf)
biasFull = fs.fileListQuery(dbFile, SQL, qdf)

# The str.join() function is needed to transform a python list into a
# comma-separated list of files that IRAF can understand.
if len(biasFull) > 1:
    gmos.gbias(','.join(str(x) for x in biasFull), 'MCbiasFull',
               **biasFlags)

# And again, for the list of 50 CenterSpec files.
biasCenSp = fs.fileListQuery(dbFile, fs.createQuery('bias', qdc), qdc)
if len(biasCenSp) > 1:
    gmos.gbias(','.join(str(x) for x in biasCenSp), 'MCbiasCenSp',
               **biasFlags)

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

# Set the task parameters.
gmos.gireduce.unlearn()
gmos.gsflat.unlearn()
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':'8',
    'logfile':'gsflatLog.txt','verbose':'no'
    }
# Normalize the spectral flats per CCD.
# Response curve fitting should be done interactively.
print ('  - CenterSpec GCAL-flat normalization, non-interactive -')
qdc['DateObs'] = '*'
cwc = {'B6-420':420.0, 'B6-520':520.0, 'B6-620':620.0,
       'R4-620':620.0, 'R4-740':740.0}
for tag,w in cwc.iteritems():
    qdc['Disperser'] = tag[0:2] + '00+_%'
    qdc['CentWave'] = w
    flatName = 'MCflatCenSp_' + tag
    flatCenSp = fs.fileListQuery(dbFile, fs.createQuery('gcalFlat', qdc), qdc)
    gmos.gsflat (','.join(str(x) for x in flatCenSp), flatName,
            bias='MCbiasCenSp', **flatFlags)

print ('  - Full Flat (GCAL & Twi) normalization, non-interactive -')
qdf['DateObs'] = '*'
cwf = {'B6-520':520.0, 'R4-740':740.0}
flatFlags.update({'fl_keep':'yes','fl_usegrad':'yes','fl_detec':'no',
                   'fl_seprows':'no','order':53})
flatType = ['gcalFlat', 'twiSpecFlat']
for ft in flatType:
    for tag,w in cwf.iteritems():
        qdf['Disperser'] = tag[0:2] + '00+_%'
        qdf['CentWave'] = w
        flatName = 'MC' + ft + '-M01_' + tag
        combName = 'MC' + ft + 'Comb-M01_' + tag
        flatFull = fs.fileListQuery(dbFile, fs.createQuery(ft, qdf), qdf)
        gmos.gsflat (','.join(str(x) for x in flatFull), flatName,
             bias='MCbiasFull', combflat=combName, **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:

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'})

print ('  - Longslit Std-star and Arc exposures -')
for tag,w in cwc.iteritems():
    qdc['Disperser'] = tag[0:2] + '00+_%'
    qdc['CentWave'] = w
    flatName = 'MCflatCenSp_' + tag
    arcCenSp = fs.fileListQuery(dbFile, fs.createQuery('arc', qdc), qdc)
    gmos.gsreduce (','.join(str(x) for x in arcCenSp), bias='MCbiasCenSp',
               **arcFlags)
    stdCenSp = fs.fileListQuery(dbFile, fs.createQuery('std', qdc), qdc)
    gmos.gsreduce (','.join(str(x) for x in stdCenSp), bias='MCbiasCenSp',
               flatim=flatName, **stdFlags)

print ('  - MOS Science and Arc exposures -')
for tag,w in cwf.iteritems():
    qdf['Disperser'] = tag[0:2] + '00+_%'
    qdf['CentWave'] = w
    flatName = 'MCgcalFlat-M01_' + tag
    gradName = 'MCgcalFlatComb-M01_' + tag
    arcFull = fs.fileListQuery(dbFile, fs.createQuery('arcP', qdf), qdf)
    gmos.gsreduce (','.join(str(x) for x in arcFull), bias='MCbiasFull',
               gradimage=gradName, **arcFlags)
    sciFull = fs.fileListQuery(dbFile, fs.createQuery('sciSpec', qdf), qdf)
    gmos.gsreduce (','.join(str(x) for x in sciFull), bias='MCbiasFull',
               flatim=flatName, gradimage=gradName, **sciFlags)

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

Multi-frame Cosmic Ray Rejection

The targets in this program were observed with two grating/cenWave configurations, 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 result will be the average of the exposures, with outlier rejection.

# Set task parameters.
gemtools.gemcombine.unlearn()
sciCombFlags = {
    'combine':'average','reject':'ccdclip',
    'fl_vardq':'yes','fl_dqprop':'yes',
    'logfile':'gemcombineLog.txt','verbose':'no'
}
# Combine the science exposures with outlier rejection for each grating.
prefix = 'gs'
for tag,w in cwf.iteritems():
    qdf['Disperser'] = tag[0:2] + '00+_%'
    qdf['CentWave'] = w
    outFile = qdf['Object'] + '-M01_' + tag
    sciFull = fs.fileListQuery(dbFile, fs.createQuery('sciSpec', qdf), qdf)
    gemtools.gemcombine (','.join(prefix+str(x) for x in sciFull), outFile,
                          **sciCombFlags)

Note the need above to explicitly propagate the DQ and VAR extensions from the input files.

Wavelength Calibration

Image rectification and wavelength linearization are performed next, using the wavelength calibrated arc lamp exposures taken immediately before each sequence of science and standard star exposures (see Wavelength Calibration). Use the default Chebyshev polynomial, but with order=5 for longslit, and order=7 for MOS Arcs.

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'
    }
# There are many arcs to choose from: we only need one for each setting.
for seq in ['246','247','248','251','252']:
    inFile = prefix + 'N20120124S0' + seq
    gmos.gswavelength(inFile, **waveFlags)

# Now for MOS Arcs
# Should achive RMS<0.20
waveFlags.update({'order':7,'nsum':20,'step':2})
for tag,w in cwf.iteritems():
    qdf['Disperser'] = tag[0:2] + '00+_%'
    qdf['CentWave'] = w
    outFile = qdf['Object'] + tag
    arcFull = fs.fileListQuery(dbFile, fs.createQuery('arcP', qdf), qdf)
    gmos.gswavelength (','.join(prefix+str(x) for x in arcFull),
                       **waveFlags)

Then apply the calibration to the standard stars and science exposures. We first construct a mapping between the grating/central wavelength, and the corresponding science and arc image sequence numbers, for each standard star. This mapping, inferred from the observing log, allows us to create more meaningful output file names.

gmos.gstransform.unlearn()
transFlags = {
    'fl_vardq':'yes','interptype':'linear','fl_flux':'yes',
    'logfile':'gstransformLog.txt','verbose':'no'
}
# (sciFile, arcFile, stdName): 'Config mneumonic'
transMap = {
            (120,247,'G191B2B'):'B6-420',
            (117,246,'G191B2B'):'B6-520',
            (121,248,'G191B2B'):'B6-620',
            (128,251,'G191B2B'):'R4-620',
            (123,252,'G191B2B'):'R4-740',
            (229,247,'HZ44'):'B6-420',
            (226,246,'HZ44'):'B6-520',
            (230,248,'HZ44'):'B6-620'
            }
for id,tag in transMap.iteritems():
    inFile = 'gsN20120124S0' + str(id[0])
    wavFile = 'gsN20120124S0' + str(id[1])
    outFile = 't' + id[2] + '_' + tag
    gmos.gstransform (inFile, outimages=outFile, wavtraname=wavFile,
                      **transFlags)

# Science MOS exposures
transFlags.update({'fl_vardq':'yes'})
gmos.gstransform ('M81-field1-M01_B6-520', wavtraname='gsN20120124S0169',
                  **transFlags)
gmos.gstransform ('M81-field1-M01_R4-740', wavtraname='gsN20120124S0187',
                  **transFlags)

Standard Star Processing

Sky Subtraction

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. Sky subtraction for MOS spectra is usually performed on each slitlet prior to spectral extraction. However gsskysub is not well matched to the nature of these sources: extended H_II regions that are not the same size nor identically centered in the slitlets. So for the MOS exposures, sky subtraction will be performed interactively during the spectral extraction step.

Sky subtraction for the longslit spectra of the standard stars could similarly be deferred to extraction. But the selection of a sky region from a spatial profile of the standard star spectrograms can be determined easily by inspection, e.g., by using the pcols task. For more information on using the pcols task, type ‘help pcols’ in an active PyRAF session.

iraf.pcols('tG191B2B_B6-520.fits[SCI]', 1030, 2040, wy1=100, wy2=500,
        wx1=80, wx2=420)

Select your preferred sky regions and perform the sky subtraction.

gmos.gsskysub.unlearn()
skyFlags = {
    'fl_oversize':'no','fl_vardq':'no','logfile':'gsskysubLog.txt',
    'verbose':'no'
}
gmos.gsskysub ('tG191B2B_*', long_sample='50:150,350:450', **skyFlags)
gmos.gsskysub ('tHZ44_*', long_sample='50:150,350:450', **skyFlags)

Sky subtraction for MOS spectra could be performed, were the size of the source known (and the same) for all slitlets.

Longslit Extraction

Extract the 1-D spectra of the standard stars from the 2-D spectrograms, using gsextract with a 3-arcsec aperture.

gmos.gsextract.unlearn()
extrFlags = {
    'apwidth':3.,'fl_inter':'no','find':'yes',
    'trace':'yes','tfunction':'spline3','tnsum':20,'tstep':50,
    'weights':'none','background':'none',
    'fl_vardq':'no','verbose':'no','logfile':'gsextractLog.txt'
}
ordc = {'B6-420':7, 'B6-520':8, 'B6-620':8, 'R4-620':19, 'R4-740':19}
# Start with G191-B2B
for tag,o in ordc.iteritems():
    inFile = 'stG191B2B_' + tag
    gmos.gsextract (inFile, torder=o, **extrFlags)

# Special fitting orders for HZ44:
gmos.gsextract ('stHZ44_B6-520', torder=6, **extrFlags)
gmos.gsextract ('stHZ44_B6-420', torder=14, refimages='stHZ44_B6-520',
                **extrFlags)
gmos.gsextract ('stHZ44_B6-620', torder=13, refimages='stHZ44_B6-520',
                **extrFlags)

Note that the refimages parameter assures that same trace will be used for subsequent central wavelength settings.

Flux Calibration

Now derive the flux calibration, using your local copy of the Mauna Kea atmospheric extinction function. The standard stars G191B2B and HZ44 were observed with each grating and a variety of central wavelengths as a part of this program. The investigators for this program in essence derived one sensitivity function for each grating (and all central wavelength settings), even though the grating efficiency functions for different central wavelength do not match perfectly. For this example we will use only G191B2B, and make a single sensitivity function from extracted spectra at all wavelength settings. You should also make a local copy the monochromatic mags file from ESO for G191B2B and HZ44; see Standards List for details.

gmos.gsstandard.unlearn()
sensFlags = {
    'fl_inter':'yes','starname':'g191b2b','caldir':'./',
    'observatory':'Gemini-North','extinction':'./mk_extinct.txt',
    'function':'spline3','order':7,'verbose':'no','logfile':'gsstdLog.txt'
}
gmos.gsstandard ('estG191B2B_B6*', sfile='std_B6', sfunction='sens_B6',
                 **sensFlags)
gmos.gsstandard ('estG191B2B_R4*', sfile='std_R4', sfunction='sens_R4',
                 **sensFlags)

Advanced Science Processing

Spectrum Extraction

The penultimate step is to extract the MOS spectra of M81. We will define extraction and sky apertures interactively, following a path similar to that described by [SM14]. The spectrograms of these targets represent perhaps the worst-case scenario for spectral extraction using gsextract which is a wrapper around the apextract.apall task, in that they contain

  • weak emission

  • no continuum

  • widely spaced emission lines

  • small or non-existent regions for sky background (owing to the small spatial extent of the slitlets)

This makes it impractical to identify extraction regions automatically. Instead, identify the location of the brightest emission features in advance using an image display server. The figure below shows a display of all the image extensions of tM81-field1_R4-740.fits in SAOImage/DS9.

../_images/M81_R4_tile.jpg

The 27 slitlet spectra with grating R400/740 of M81 H_II regions, displayed simultaneously in SAOImage/DS9 (columns 1, 4, and 7). Note that the VAR and DQ frames are also displayed (columns 2-3, 5-6, 8-9), though they are not relevant for this exercise. Click image to enlarge.

To display the brightest emission regions efficiently from within DS9,

  • Select File \(\rightarrow\) Open as \(\rightarrow\) Multiple Extension Frames… and load the file

  • Bring up the Zoom \(\rightarrow\) Pan Zoom Rotate Parameters… menu and enter the wavelength of an emission line (\(\mathrm{H}\alpha\) 6563 in this case) in the Pan dialog box

  • Select Frame \(\rightarrow\) Match \(\rightarrow\) Frame \(\rightarrow\) WCS

For each SCI image in turn, use the cursor to determine the image line (i.e., x-coordinate) of the feature of interest and the spatial extent (y-coordinate) of the desired extraction window, and record these parameters for use during the spectral extraction.

Now extract the spectra from the slitlets interactively. If you are not already comfortable using IRAF’s apextract package and the associated cursor commands, see APEXTRACT Summary. Note that with such weak continuum emission, it is best to use a very low-order spectrum trace function. The fl_vardq flag is turned off, since the VAR plane is not used to estimate the errors in the gsextract task.

extrFlags.update({'apwidth':2.5,'mos_bsample':1.0,'torder':1,'tnsum':50,
                  'tfunction':'spline3','background':'median',
                  'fl_inter':'no','fl_vardq':'no'})
gmos.gsextract ('tM81-field1-M01_B6-520', **extrFlags)
gmos.gsextract ('tM81-field1-M01_R4-740', **extrFlags)

Flux Calibration

Flux calibration is a necessary final step for this program’s science goals.

Calibrate the target spectra by applying the atmospheric extinction correction (download: mk_extinct.txt) and the sensitivity function.

gmos.gscalibrate.unlearn()
calibFlags = {
    'extinction':'./mk_extinct.txt','fl_ext':'yes','fl_scale':'no',
    'fl_vardq':'no','logfile':'gscalibrateLog.txt'
    }
gmos.gscalibrate ('etM81-field1-M01_B6-520', sfunc='sens_B6', **calibFlags)
gmos.gscalibrate ('etM81-field1-M01_R4-740', sfunc='sens_R4', **calibFlags)

A calibrated spectrum of one of the brighter MOS targets is shown below.

../_images/M81_targ107.png

H_II region target 107 in M81 with grating B600/520 in M81. The Balmer recombination lines are visible, as are prominent lines of [O_III], [N_II], and [S_II]. Click image to enlarge.