Spectral extraction and wavelength calibration

See also: Supernova Spectroscopy with FAST | Flux calibration

Contents: Installation | Test data | CCD processing | Spectral extraction and wavelength calibration

The reduction up to the extraction of 1D wavelength-calibrated spectra is done with standard procedures in IRAF that are encoded in two driver cl scripts fastred.cl and fapall.cl. The overall procedure is described at length in Sec. 2.1 of the reductio.ps file contained in the redux/ directory created during the installation process of the Flux calibration IDL routines.

Installation

You will need the set of cl scripts, which you can download as tarballs locally at: 
/data/www/projects/supernova/SNfast/scripts.tar
Copy both tarballs to your IRAF directory and unpack them there. This will create a scripts/ directory that contains all the cl scripts.

You will also need to tell IRAF where to find those scripts and how to call them by adding the following lines to your login.cl file (cf. ~/pchallis/iraf/login.cl): 

task    convert=home$scripts/convert.cl
task    interstat=home$scripts/interstat.cl
task    nickelred=home$scripts/nickelred.cl
task    nickelscan=home$scripts/nickelscan.cl
task    lrisdivide=home$scripts/lrisdivide.cl
task    lrisfixhead=home$scripts/lrisfixhead.cl
task    kastfixhead=home$scripts/kastfixhead.cl
task    uvfixhead=home$scripts/uvfixhead.cl
task    flwofixhead=home$scripts/flwofixhead.cl
task    mmtfixhead=home$scripts/mmtfixhead.cl
task    ssofixhead=home$scripts/ssofixhead.cl
task    vltfixhead=home$scripts/vltfixhead.cl
task    optpa=home$scripts/optpa.cl
task    sidtime=home$scripts/sidtime.cl
task    xzap=home$scripts/xzap.cl
task    qzap=home$scripts/qzap.cl
task    fileroot=home$scripts/fileroot.cl
task    iterstat=home$scripts/iterstat.cl
task    minv=home$scripts/minv.cl
task    fgain=home$scripts/fgain.cl
task    szap=home$scripts/szap.cl
task    makemask=home$scripts/makemask.cl
task    baseline=home$scripts/baseline.cl
task    rdtape=home$rdtape/rdtape.cl
task    flatcor=home$scripts/flatcor.cl
task    crremove=home$scripts/crremove.cl
task    twoampbias=home$scripts/twoampbias.cl
task    medisplay=home$scripts/medisplay.cl
task    merfits=home$scripts/merfits.cl
task    capall=home$scripts/capall.cl
task    lcofixhead=home$scripts/lcofixhead.cl
task    mdisplay=home$scripts/mdisplay.cl
task    $fastred=home$scripts/fastred.cl
task    $fastred_binby4=home$scripts/fastred_binby4.cl
task    $fastredd=home$scripts/fastredd.cl
task    $fastred2=home$scripts/fastred2.cl
task    $fapall=home$scripts/fapall.cl
task    $fapall98=home$scripts/fapall98.cl
task    $fapall_binby4=home$scripts/fapall_binby4.cl
task 	$fid=home$scripts/fid.cl

Test data

In what follows we will go through one example of spectral extraction on FAST spectra of supernovae. These data were taken on 29 Mar 2008 with the standard Combo setup (300 gpm grating and 3'' slit).

Start by copying the data over to a directory of your choice (hereafter called <REDUCE_DIR>): 

~> cp /data/www/projects/supernova/SNfast/raw/*.fits <REDUCE_DIR>
NOTE: You can compare your output files at any stage of the tutorial with those located in /data/www/projects/supernova/SNfast/reduced/

Now start an xgterm and start IRAF by typing cl (or ecl). Then go to <REDUCE_DIR>. 
~> ecl
    NOAO Sun/IRAF Revision 2.12EXPORT Thu May  2 19:39:56 MST 2002
This is the EXPORT version of Sun/IRAF V2.12 for SunOS 4 and Solaris 2.8
 
            Center For Astrophysics Telescope Data Center
 ***********************************************************************
 **       IRAF version 2.12.2a      installed July 11, 2004           ** 
 **       RVSAO version 2.4.8       installed March 10, 2006          ** 
 **       SAOTDC version 3.9        installed May 10, 2004            ** 
 **       STSDAS version 3.5        installed March 29, 2006          ** 
 **       TABLES version 3.5        installed March 29, 2006          ** 
 **       MSCRED version 4.8        installed January 9, 2006         **
 **       WCSTools version 3.6.2    installed December 5, 2005        **
 **       STECF version 1.5         installed December 6, 2005        **
 **       GEMINI version 1.8        installed November 30, 2005       **
 **       SPECTIME version 2.0      installed September 17, 2002      **
 **                                                                   **
 **            In case of problems, contact dmink@cfa                 ** 
 ***********************************************************************

 Welcome to IRAF.  To list the available commands, type ? or ??.  To get
 detailed information about a command, type `help command'.   To  run  a
 command  or  load  a  package,  type  its name.   Type  `bye' to exit a
 package, or `logout' to get out of the CL.   Type `news'  to  find  out
 what is new in the version of the system you are using.   The following
 commands or packages are currently defined:
    +-------------------------------------------------------------+
    |           WCSTools World Coordinate System Utilities        |
    |             Smithsonian Astrophysical Observatory           |
    |                     Telescope Data Center                   |
    |                 Version 3.6.2, 21 July 2005                 |
    +-------------------------------------------------------------+


#-----------------------------------------------------------+
#           RVSAO Radial Velocity Analysis Package          |
#           Smithsonian Astrophysical Observatory           |
#                   Telescope Data Center                   |
#                Version 2.5.8 June 13, 2007                |
#-----------------------------------------------------------+

Resetting terminal type to xgterm...
      color.      guiapps.    mscred      proto.      starbase    tables35
      ctio.       hectospec   nicproto.   r2rvsao.    stecf.      upsqiid
      dataio.     ifocas.     nmisc.      rtest.      stsdas.     utilities.
      dbms.       images.     noao.       rvsao.      stsdas35.   wcstools.
      fitsutil    imcnv       obslog      saotdc.     svdfit      xdimsum.
      gemini.     language.   obsolete.   softools.   system.     xray.
      gmisc.      lists.      plot.       spectime.   tables.     xrv.

ecl> cd <REDUCE_DIR>
Now you should be all set to run the first driver script fastred.cl.

Back to top

CCD processing (fastred.cl)

The first script fastred.cl trims the 2D frames, fixes bad columns, applies an overscan correction, creates a master flat and response function, flatfields the 2D frames, and adds a few keywords to the fits headers. The only interactive part occurs when fitting the response curve. In what follows we will go through the contents of the fastred.cl script.

NOTE: What follows is a description of the contents of the fastred.cl script. To execute this script simply type 'fastred' at the cl> prompt (and skip to here).

First, all DARK and BIAS frames are deleted, since we do not use them in the reduction process (dark current is irrelevant for FAST, and the bias is constant): 

imdel ("*BIAS*",
yes, verify=no, default_acti=yes)
imdel ("*DARK*",
yes, verify=no, default_acti=yes)
The 2D frames are then trimmed, corrected for overscan, and bad columns fixed: 
ccdproc ("*.fits",
output="", ccdtype=" ", max_cache=0, noproc=no, fixpix=yes, overscan=yes,
trim=yes, zerocor=no, darkcor=no, flatcor=no, illumcor=no, fringecor=no,
readcor=no, scancor=no, readaxis="line", fixfile="/home/pchallis/iraf/lib/badcolumns.pix", biassec="[1:20,*]",
trimsec="[34:2716,2:158]", zero="", dark="", flat="", illum="", fringe="",
minreplace=1., scantype="shortscxan", nscan=1, interactive=no,
function="legendre", order=1, sample="*", naverage=1, niterate=1,
low_reject=3., high_reject=3., grow=0.)
Notice the fixfile parameter points to the file /home/pchallis/iraf/lib/badcolumns.pix; the safest thing here would be to copy this file to your IRAF directory and edit fastred.cl accordingly.

NOTE: Make sure you check the trimsec and biassec settings each time the FAST CCD is upgraded!

We then create a master flat from which we derive the response curve used to flatfield the data at a later stage. After this the original FLAT frames are deleted: 
flatcombine ("*FLAT*",
output="Flat", combine="average", reject="avsigclip", ccdtype=" ", process=no,
subsets=no, delete=no, clobber=no, scale="mode", statsec="", nlow=1, nhigh=1,
nkeep=1, mclip=yes, lsigma=3., hsigma=3., rdnoise="0.", gain="1.",
snoise="0.", pclip=-0.5, blank=1.)

response ("Flat",
"Flat", "flatres", interactive=yes, threshold=INDEF, sample="*", naverage=1,
function="spline3", order=19, low_reject=0., high_reject=0., niterate=1,
grow=0., graphics="stdgraph", cursor="")

imdel ("*FLAT*",
yes, verify=no, default_acti=yes)
The next step is to add/edit the following fits header keywords with the flwofixhead.cl script: OBSNUM (CCD frame number), OBSERVAT (set to "flwo"), DISPAXIS (set to 1), ST (sidereal time), EPOCH, OPT_PA (optimal parallactic angle), AIRMASS: 
flwofixhead ("*.fits",
imglist="")
NOTE: Make sure you edit flwofixhead.cl to point to the correct cmds.asthedit file (in your own IRAF directory): 
asthedit(img,'/home/pchallis/iraf/scripts/cmds.asthedit',upd+,verbose+,oldstyl-)
Last, a list of all science spectra and standards (including COMPs) is created with the name objects and these frames are flatfielded: 
!/home/pchallis/iraf/scripts/noflat 

ccdproc ("@objects",
output="", ccdtype=" ", max_cache=0, noproc=no, fixpix=no, overscan=yes,
trim=yes, zerocor=no, darkcor=no, flatcor=yes, illumcor=no, fringecor=no,
readcor=no, scancor=no, readaxis="line", fixfile="", biassec="[2700:2718,*]",
trimsec="[9:2683,2:158]", zero="", dark="", flat="flatres", illum="",
fringe="", minreplace=1., scantype="shortscan", nscan=1, interactive=no,
function="legendre", order=1, sample="*", naverage=1, niterate=1,
low_reject=3., high_reject=3., grow=0.)
NOTE: Edit the path /home/pchallis/iraf/scripts/noflat in fastred.cl to point to the correct file in your own IRAF directory.

Let's now see how this all works in practice. Execute the fastred.cl script by typing fastred at the cl prompt:
ecl> fastred
This will take a few moments before you are prompted to adjust the fit to the response curve: 
Fit the normalization spectrum for Flat interactively (yes):
Answer 'yes' here and examine the plot in the graphics terminal (plot on the left below):



If it looks any different than the plot on the left, something has messed up. A good place to start is to run imstat on all the FLAT frames (i.e. ecl> imstat *FLAT*fits) and check for deviant flat frames. The plot on the right (above) is the ratio of the flatfield spectrum to the cubic spline fit (obtained by hitting 'k' in the graphics terminal). Residuals of ~2% are fine.

Once you exit the interactive fitting procedure ('q'), the fits headers are updated with the flwofixhead.cl script and a lot of screen output follows. Below is the total output for one spectrum, 0117.sn2008bf.fits: 
add 0117.sn2008bf.fits,OBSNUM = 117
0117.sn2008bf.fits updated
0117.sn2008bf.fits,OBSERVAT: flwo1 -> flwo
0117.sn2008bf.fits updated
0117.sn2008bf.fits,DISPAXIS: 1 -> 1
0117.sn2008bf.fits updated
0117.sn2008bf.fits:
  $I = 0117.sn2008bf.fits
  $D = 16/07/08
  $T = 17:31:32
  $GMD = 2008-07-16
  $GMT = 21:31:32
  $GMDT = 2008-07-16T21:31:32
  st = 12:44:05 -> 12:24:08.42
  epoch = 2000.0 -> 2008.243127709331
add 0117.sn2008bf.fits,OPT_PA = 28.55753
0117.sn2008bf.fits updated
#              Image    UT middle  effective  begin   middle     end   updated
# SETAIRMASS: Observatory parameters for Whipple Observatory
#       latitude = 31:40:51.4
  0117.sn2008bf.fits    7:25:27.0   1.0277   1.0234   1.0274   1.0331  yes
That's it. Now we need to extract the 1D spectra and apply a wavelength calibration. This is done with the fapall.cl script.

Back to top

Spectral extraction and wavelength calibration (fapall.cl)

The second script fapall.cl extracts the 1D spectra and applies a wavelength calibration. It is far more interactive than fastred.cl. In what follows we will go through the contents of the fapall.cl script.

NOTE: What follows is a description of the contents of the fapall.cl script. To execute this script simply type 'fapall' at the cl> prompt (and skip to here).

First, several files are created with a shell script fapallsort.csh. These include: 

 File        Description
----------------------------------------------------------------
 ap.list     all OBJECT files *.fits
 objects     all OBJECT files 1D multispec extractions *.ms.fits
 comps       all COMP files
----------------------------------------------------------------
 firstcompid.cl identify first COMP based on idfast-comp.ms in $home/iraf/database
 hed.cl         edit REFSPEC1 keyword in OBJECT files to associate 
                the correct COMP
 cap.cl         extract apertures in COMP files using the same aperture
                definition as the associated OBJECT
You should copy the fapallsort.csh script to your bin/ directory and edit fapall.cl to point to this file. 
!/home/pchallis/bin/fapallsort.csh
The 1D spectra are then extracted in multispec format using the standard apall task: 
apall ("@ap.list",
1, output="", apertures="", format="multispec", references="", profiles="",
interactive=yes, find=yes, recenter=yes, resize=no, edit=yes, trace=yes,
fittrace=yes, extract=yes, extras=yes, review=yes, line=INDEF, nsum=50,
lower=-3.5, upper=3.5, apidtable="", b_function="legendre", b_order=2,
b_sample="-25:-10,10:25", b_naverage=-100, b_niterate=3, b_low_reject=3.,
b_high_rejec=3., b_grow=0., width=5., radius=10., threshold=0., minsep=5.,
maxsep=1000., order="increasing", aprecenter="", npeaks=INDEF, shift=yes,
llimit=INDEF, ulimit=INDEF, ylevel=0.1, peak=yes, bkg=yes, r_grow=0.,
avglimits=no, t_nsum=300, t_step=50, t_nlost=3, t_function="legendre",
t_order=4, t_sample="*", t_naverage=1, t_niterate=0, t_low_reject=3.,
t_high_rejec=3., t_grow=0., background="fit", skybox=1, weights="variance",
pfit="fit1d", clean=yes, saturation=INDEF, readnoise="8.6", gain="1.2",
lsigma=4., usigma=4., nsubaps=1)
Next, the cap.cl and hed.cl scripts created by fapallsort.csh are executed. This will extract traces from the COMP spectra similar to the corresponding science spectra and edit the fits headers of the science spectra to point to the correct COMP lamp spectrum: 
cl < cap.cl
cl < hed.cl
Follows the determination of the wavelength solution. This is simply done by copying over a previous wavelength solution (you may need to update this occasionally, although I have not done so since 16 Sep 2006) and running the IRAF task reidentify: 
!cp /home/pchallis/iraf/database/idfast-comp.ms ./database/
!cp /home/pchallis/iraf/fast-comp.ms.fits .

reidentify ("fast-comp.ms",
"@firstcomp", "yes", " ", " ", interactive="yes", section="middle line",
newaps=yes, override=no, refit=yes, trace=no, step="10", nsum="10",
shift="0.", search=0., nlost=0, cradius=5., threshold=0., addfeatures=no,
coordlist="/home/pchallis/iraf/linelists/flwolist.dat", match=-3.,
maxfeatures=50, minsep=2., database="database", logfiles="logfile",
plotfile="", verbose=yes, graphics="stdgraph", cursor="", aidpars="")	

cl < firstcompid.cl
reidentify (s1,
"@comps", "yes", " ", " ", interactive="yes", section="middle line",
newaps=yes, override=no, refit=yes, trace=no, step="10", nsum="10",
shift="0.", search=0., nlost=0, cradius=5., threshold=0., addfeatures=no,
coordlist="/home/pchallis/iraf/linelists/flwolist.dat", match=-3.,
maxfeatures=50, minsep=2., database="database", logfiles="logfile",
plotfile="", verbose=yes, graphics="stdgraph", cursor="", aidpars="")
NOTE: You will need to copy /home/pchallis/iraf/database/idfast-comp.ms and /home/pchallis/iraf/fast-comp.ms.fits to your IRAF directory and edit fapall.cl accordingly.

Last, the wavelength solution is applied to the 1D extracted spectra with the IRAF task dispcor, and a list (called dlist) of all wavelength-calibrated spectra d*.ms.fits is created: 
dispcor ("@objects",
"d//@objects", linearize=yes, database="database", table="", w1=INDEF,
w2=INDEF, dw=INDEF, nw=INDEF, flux=no, samedisp=no, global=no,
ignoreaps=yes, confirm=no, listonly=no, verbose=yes, logfile="")

!ls -1 d*.ms.fits > dlist
Let's now see how this all works in practice. Execute the fapall.cl script by typing fapall at the cl prompt:
ecl> fapall
This will run the fapallsort.csh script, producing various terminal output, including object names, number of corresponding spectra and COMPs, and coordinates: 
ecl> fapall
 Created file ap.list
 Created file objects
 Created file comps
 Created file firstcompid.cl
09:48:54 13:45:13.80
HD84937: 3 files and 1 COMP(s)
10:39:36 43:06:10.98
Feige34: 1 files and 1 COMP(s)
12:04:02 +20:14:42.6
sn2008bf: 1 files and 1 COMP(s)
12:30:40 +41:38:14.5
sn2008ax: 1 files and 1 COMP(s)
13:02:58 +10:30:27.0
sn2008bm: 1 files and 1 COMP(s)
 Created file hed.cl
 Created file cap.cl
 DONE -- successfully executed fapallsort
You will then be prompted as you would normally when running apall in interactive mode: 
Find apertures for 0107.HD84937?  (yes): 
Edit apertures for 0107.HD84937?  (yes):
Below are plots showing the spatial profile of the spectral trace (left), the change of the trace position along the CCD columns (right), and the 1D extracted spectrum (in ADUs vs. pixel; right), for 0117.sn2008bf.fits:



Once you have done this for each spectrum, the cap.cl and hed.cl scripts created by fapallsort.csh are executed, producing the following output: 
add 0107.HD84937.ms.fits,REFSPEC1 = 0110.COMP.ms
0107.HD84937.ms.fits updated
add 0108.HD84937.ms.fits,REFSPEC1 = 0110.COMP.ms
0108.HD84937.ms.fits updated
add 0109.HD84937.ms.fits,REFSPEC1 = 0110.COMP.ms
0109.HD84937.ms.fits updated
add 0111.Feige34.ms.fits,REFSPEC1 = 0112.COMP.ms
0111.Feige34.ms.fits updated
add 0117.sn2008bf.ms.fits,REFSPEC1 = 0118.COMP.ms
0117.sn2008bf.ms.fits updated
add 0119.sn2008ax.ms.fits,REFSPEC1 = 0120.COMP.ms
0119.sn2008ax.ms.fits updated
add 0121.sn2008bm.ms.fits,REFSPEC1 = 0122.COMP.ms
0121.sn2008bm.ms.fits updated
Now for the wavelength solutions. There are 13 reference features in the idfast-comp.ms file, and you should check that all 13/13 are found and fit during reidentify: (If not, you can determine the wavelength solution interactively in the usual way.) 
REIDENTIFY: NOAO/IRAF V2.12.1-EXPORT sblondin@arequipa Wed 17:38:43 16-Jul-2008
  Reference image = fast-comp.ms, New image = 0110.COMP.ms, Refit = yes
          Image Data    Found     Fit Pix Shift  User Shift  Z Shift      RMS
0110.COMP.ms - Ap 1    13/13   13/13       -1.9       -2.79  -5.1E-4    0.126
0110.COMP.ms - Ap 1    13/13   13/13       -1.9       -2.79  -5.1E-4    0.126

REIDENTIFY: NOAO/IRAF V2.12.1-EXPORT sblondin@arequipa Wed 17:39:49 16-Jul-2008
  Reference image = 0110.COMP.ms.fits, New image = 0112.COMP.ms, Refit = yes
          Image Data    Found     Fit Pix Shift  User Shift  Z Shift      RMS
0112.COMP.ms - Ap 1    13/13   13/13      0.326       0.481  8.82E-5     0.13
0112.COMP.ms - Ap 1    13/13   13/13      0.326       0.481  8.82E-5     0.13
  Reference image = 0110.COMP.ms.fits, New image = 0118.COMP.ms, Refit = yes
          Image Data    Found     Fit Pix Shift  User Shift  Z Shift      RMS
0118.COMP.ms - Ap 1    13/13   13/13      0.221       0.326  5.81E-5    0.153
0118.COMP.ms - Ap 1    13/13   13/13      0.221       0.326  5.81E-5    0.153
  Reference image = 0110.COMP.ms.fits, New image = 0120.COMP.ms, Refit = yes
          Image Data    Found     Fit Pix Shift  User Shift  Z Shift      RMS
0120.COMP.ms - Ap 1    13/13   13/13      0.616       0.908  1.61E-4    0.115
0120.COMP.ms - Ap 1    13/13   13/13      0.616       0.908  1.61E-4    0.115
  Reference image = 0110.COMP.ms.fits, New image = 0122.COMP.ms, Refit = yes
          Image Data    Found     Fit Pix Shift  User Shift  Z Shift      RMS
0122.COMP.ms - Ap 1    13/13   13/13    0.00791      0.0118  7.90E-7     0.12
0122.COMP.ms - Ap 1    13/13   13/13    0.00791      0.0118  7.90E-7     0.12
Finally, the wavelength solutions are applied to each corresponding spectrum: 
0107.HD84937.ms.fits: REFSPEC1 = '0110.COMP.ms 1.'
d0107.HD84937.ms.fits: ap = 1, w1 = 3476.694, w2 = 7420.933, dw = 1.470633, nw = 2683
0108.HD84937.ms.fits: REFSPEC1 = '0110.COMP.ms 1.'
d0108.HD84937.ms.fits: ap = 1, w1 = 3476.694, w2 = 7420.933, dw = 1.470633, nw = 2683
0109.HD84937.ms.fits: REFSPEC1 = '0110.COMP.ms 1.'
d0109.HD84937.ms.fits: ap = 1, w1 = 3476.694, w2 = 7420.933, dw = 1.470633, nw = 2683
0111.Feige34.ms.fits: REFSPEC1 = '0112.COMP.ms 1.'
d0111.Feige34.ms.fits: ap = 1, w1 = 3476.355, w2 = 7420.597, dw = 1.470634, nw = 2683
0117.sn2008bf.ms.fits: REFSPEC1 = '0118.COMP.ms 1.'
d0117.sn2008bf.ms.fits: ap = 1, w1 = 3476.527, w2 = 7420.627, dw = 1.470582, nw = 2683
0119.sn2008ax.ms.fits: REFSPEC1 = '0120.COMP.ms 1.'
d0119.sn2008ax.ms.fits: ap = 1, w1 = 3475.845, w2 = 7419.648, dw = 1.470471, nw = 2683
0121.sn2008bm.ms.fits: REFSPEC1 = '0122.COMP.ms 1.'
d0121.sn2008bm.ms.fits: ap = 1, w1 = 3476.808, w2 = 7420.897, dw = 1.470578, nw = 2683
At this stage, you have a set of 1D extracted spectra that are also wavelength calibrated. The files are named d*.ms.fits and are stored in the dlist file created by fapall: 
~> cat dlist
d0107.HD84937.ms.fits
d0108.HD84937.ms.fits
d0109.HD84937.ms.fits
d0111.Feige34.ms.fits
d0117.sn2008bf.ms.fits
d0119.sn2008ax.ms.fits
d0121.sn2008bm.ms.fits
Now we need to flux calibrate the spectra. This is done with a separate set of IDL routines written by Tom Matheson. See this page for scripts and examples.

Back to top

Last updated: 21 Jul 2008
Stéphane Blondin