Change page style: 

Example Michelle Q Band Reduction

Example Q-band Spectroscopic Data Reduction



This page gives an example of using the msreduce script to reduce a Q-band spectrum. It is based upon the experience of Kevin Volk (who uses the first person in this page where needed, mixed in with use of the third person to refer to the reader).

In this instance I have an observation of a bright science target and an observation of a standard star both at Q-band in the lowres20 spectral mode of Michelle. The observation was made to study the variations in a strong feature across the extended source. For the purposes of discussion here only the overall spectrum across the entire source will be extracted. The Q-band spectrum of this object is well known, as the object has been observed from space with the Infrared Space Observatory (ISO) and previous to that by the Infrared Astronomical Satellite (IRAS). Both these satellites obtained low resolution spectra of the object in question. The ISO spectrum is of better quality (at least in terms of providing continuous wavelength coverage) than anything that we can do from the ground. Those satellite observations did not provide any spatial information, however.

The raw data file for the science object observation is N20051229S0208.fits, while the raw data file for the standard star observation is N20051229S0211.fits. Both these files can be obtained from the Gemini Science Archive as they are now public. The science target is HD 56126, more often referred to as IRAS 07134+1005. The standard star is Procyon, alpha CMi. This is not one of the usual calibration stars for which the midir package has a flux curve, so this example carries out the reduction on the pretext that the standard star was alpha CMa, Sirius, instead. I will then need to scale the output spectrum down to allow for the difference in brightness between alpha CMi and alpha CMa, but for the purposes of illustration here that is not important. In Q-band the difference in effective spectral shape between an F5IV star, Procyon, and an A1V star, Sirius, is very small. The difference in brightness at these wavelengths can be found from the ratio of the IRAS [25] fluxes, for example, or from the ISO spectra of these objects.

Initial Steps

The first step is run nsheaders to set the instrument parameters for Michelle. This routine is found in the gnirs part of the Gemini package, but it is available once you invoke the midir package as well

cl> gemmini
ge> midir
mi> nsheaders michelle
NSHEADERS -- Mon Nov 27 12:47:12 HST 2006
setting header keywords
setting nscut.section
setting nsreduce.section
resetting nswavelength.dispaxis
setting nswavelength.coordlist, nswavelength.threshold, and
setting the default logfile

Before starting the reduction it is useful to look at the raw data and see (a) what the atmospheric emission lines look like, so you can get a rough ID of the features, and (b) whether the spectrum needs to be defringed. This can be done from the raw data files if you wish, as in

cl> display rawdata/N200512S0211[1][*,*,3,1] zs- zr+ nsample=76800
cl> imexam
cl> display rawdata/N200512S0211[1][*,*,1,1] zs- zr+ nsample=76800
cl> imexam

The first display command shows the raw data from one nod position, the first nod A position, for the object plus the sky. The way this looks in ds9 is as follows --

Raw NOD frame on the source

You can see the strong atmospheric bands as vertical bright lines in the image (the spectral direction is from left to right with short wavelengths to the left, and the up/down direction is along the slit), and the object is much too faint to be seen in the raw image. Using the IRAF imexam task you can do a line cut with the "l" cursor key to get a cut across the image--

Line cut of raw source frame

This figure has been made with the y range set from 0 to 45000 ADU to show the true line to continuum contrast in fairly dry conditions at Mauna Kea.

The second display command displays the difference spectral image from the same nod, wherein the object spectrum can be seen--

Single NOD difference frame (chop A - chop B)

The raw spectrum is interrupted frequently by the atmospheric bands where the background is high and the atmospheric transmission is low. A line cut along the peak of the spectrum looks like this--

Line cut along the raw difference spectrum

where you see the atmospheric lines as absorption lines instead of emission lines. There is no fringing in this line cut, so defringing of the spectrum is not needed.

Of note is that the band at about 20 microns, near pixel 135, has split into two peaks at this resolution, and so it is not usable for the wavelength calibration.

Starting the Processing

I have set the parameters for nsreduce as follows: (these should work for you except that the rawpath value need not be as indicated here)

midir> lpar msreduce
inspec = "N20051229S0208" Name of source spectrum raw data file(s)
(outtype = "fnu")          Type of output spectrum: fnu|flambda|lambda*fla
(rawpath = "rawdata")      Path for raw images
(line = 160)            Line to use to identify the aperture
(fl_std = yes)            Extract both source and standard spectra, do te
(std = "N20051229S0211") Name of standard spectrum raw data file(s)
(stdname = "alphaCMa")     Name of standard star
(fl_bbody = no)             Use blackbody option rather than spectrophotome
(bbody = 10000.)         Temperature of calibrator for black-body fit
(fl_flat = no)             Apply flat/bias to the raw data file(s) before 
(flat = "")             Input flat-field image(s) (raw data file(s))
(bias = "")             Input bias image(s) (raw data file(s))
(fl_extract = yes)            Run nsextract interactively
(fl_telluric = yes)            Run msabsflux interactively
(fl_wavelengt = yes)            Run nswavelength interactively
(fl_transform = no)             Run nstransform interactively
(fl_retrace = no)             Retrace spectrum, if appropriate
(fl_process = yes)            Do initial processing with mprepare
(fl_clear = yes)            Clear database sub-directory
(fl_negative = no)             Extract from a negative spectrum (NOD mode)
(fl_defringe = no)             Defringe the spectra (an intrinsicly interactiv
(fl_lowres = yes)            Low resolution spectrum [msdefringe]
(fmin = 31)             Start of filtering region (pixels) [msdefringe]
(fmax = 48)             End of filtering region (pixels) [msdefringe]
(fl_dfinterp = yes)            Interpolate across masked region [msdefringe]
(fl_zerocut = no)             Mask out negative points [msdefringe]
(fl_reextract = no)             Just re-extract the spectra
(linelist = "")             Line list file name
(logfile = "spectrum.log") Log file name
(verbose = yes)            verbose logging?
(status = 101)            Exit error status: (0=good, >0=bad)
(scanfile = "")             Internal use only
(mode = "ql")           

Note that the parameters are set with the standard star name and raw data file name for the standard star spectrum, the raw data file name for the object, a path to these two raw data files, and various of the flags are set. In this case there is no need to do defringing so fl_defringe is set to "no". Most of the other parameters have their default values. The Michelle lowres20 Q-band spectra generally do not show any fringes and do not need defringing, as is the case here.

The call to msreduce can either be done by editing the parameters to the values listed above or by the command

cl> msreduce N20051229S0208 std=N20051229S0211 stdname="alphaCMa"\
fl_defringe- rawpath="rawdata"

provided that the various flags were set properly using the IRAF command epar. In a script the form of the call would be

msreduce ("N20051229S0208",
outtype="fnu", rawpath="rawdata", line=160, fl_std=yes, std="N20051229S0211", 
stdname="alphaCMa",fl_bbody=no, bbody=10000., fl_flat=no, flat="", bias="", 
fl_extract=yes,fl_telluric=yes, fl_wavelengt=yes, fl_transform=yes, fl_retrace=no,
fl_process=yes, fl_clear=yes, fl_negative=no, fl_defringe=no, fl_lowres=yes,
fmin=18, fmax=32, fl_dfinterp=yes, fl_zerocut=no, fl_reextract=no,
linelist="gnirs$data/sky.dat", logfile="example.log", verbose=yes)

with the rawpath parameter having to be set to whatever directory is appropriate.

Wavelength Calibration

When you run msreduce the task proceeds automatically until the wavelength calibration is to be done. It starts by taking the raw data files for the source and for the standard star and stacking the frames to produce a sky frame, to be used for the wavelength calibration, and a difference frame, from which the spectrum is to be extracted. These steps need no interactive input. You can instruct msreduce to skip the stacking of the raw data file by setting the fl_process flag to "no". Ny own preference is to start clean each time in a new directory each time I do a reduction, in which case this initial step is required. You can also set the flags to apply flat and bias frames to the data--these are set by the fl_flat flag. If this is set to "yes" then the appropriate raw file names need to be put in the flat and bias fields. In principle applying the flat and bias frames will improve the extraction of the spectra of extended sources. I have not generally found that this produces superior results myself. If you are reducing extended source spectra, I recommend experimenting with the flat and bias options to see whether these produce good results or not.

The reduction proceeds to the wavelength calibration step next. The tasks attempt to produce an automatic line identification of a cut across the middle of the sky frame. These automatic line IDs are always poor in my experience (see below) so you have to reset the line identifications to the start with the "i" key command, and then manually identify the lines/bands. One marks each peak by moving the cursor near it and then hitting the "m" key, and entering the wavelength in Angstroms when prompted. The proper line IDs to use are listed below.

Result of the Automatic line IDs:

Note that if you get something that looks like this for the initial wavelength solution

where the central wavelength is 0 angstroms, it means that the nsheaders command has not been used to set the parameters of the spectroscopy. But even with the parameters set the wavelength solution is way off; the central wavelength is more or less correct but the dispersion is much too large leading to an absurd wavelength range from less than 5 microns to more than 35 microns, while the actual values are roughly from 16 to 24 microns.

The approximate line positions:

position     Wavelength
pixel 52     175860
pixel 75     182980
pixel 87     186480
pixel 98     190150
pixel 108    193100
pixel 141    203220
pixel 152    206530
pixel 169    211750
pixel 191    218550
pixel 215    225950
pixel 228    229410
pixel 234    231990

I find it helpful to start by identifying two lines/bands, and doing a linear fit over the wavelength range. This produces a close enough wavelength solution to allow the rest of the lines to be automatically IDed when marked thereafter. To select a linear fit use the commands :order 2 and "f" to refit. I used the lines at roughly pixels 52 and 234 for this purpose.

Initial linear fit from two line IDs:

Proper line IDs:

The wavelength calibration is good to within a few hundred Angstroms, or a few times 0.01 micron. If you need higher precision than this you would have to work harder on the wavelength calibration aided by an atmospheric emission model at the proper resolution.

Once the initial wavelength calibration is carried out the nstransform routine finds the same sky lines for different positions along the slit. This step generally can be done using the automatic re-identification, but if you are taking the spectrum of an extended source and wants the best possible wavelength solution for all positions along the slit you may need to do this interactively. Michelle does not have a lot of curvature of the sky lines over the array so the automatic line re-identification normally works well; thus when prompted about whether you wish to examine the line identifications you can usually answer "N" and let the software trace the lines without any manual interaction.

The same process is then followed for the standard star, with the sky lines being at the same places as listed above, resulting in a wavelength calibrated sky spectrum for file N20051229S0211.

Wavelength calibrated sky spectrum (calibration star):

Defining the Apertures for the Spectra

Next you have to choose the apertures for the object and for the standard, in this case these were chosen to get all the flux. The standard star is of course a point source, but the science target is somewhat extended on the slit and I have chosen a wide aperture to get all the flux. The aperture is close to 50 pixels wide, extending +18 and -30 pixels from the central position. You can either use the "l" and "u" keys to mark these limits, as was done here, or use the colon commands :lower -30 and :upper 18 to set these values.

When selecting the aperture I also select the background regions to the right and left of the aperture. This is done using the "b" key command to go into the background selection window, then using "t" key command to clear any automatically selected background regions. I then use the "s" key to select reasonable areas for the background and type "f" to fit these. The following three images show the automatic background you start with for HD 56126, the result of the manual selection of background regions for HD 56126, and the manually selected background for the calibration star.

Initial background for the science target spectrum

Manually defined background for the science target spectrum

Manually defined background for the calibration star spectrum

Aperture for the science target spectrum

Once the aperture is defined it is then traced. For Michelle (or T-ReCS) the tracing does not show much curvature across the array in any of the spectral modes. Thus a linear trace or a low order polynomial trace is enough. You generally needs to remove a few points from the fit in regions where there is no signal. This is done using the "x" key command, followed by using the "f" key command to re-fit the trace. The default order of the fitting in the GNIRS package is higher than is needed for Michelle, as shown in the initial tracing of the science target or of the calibration star shown below. Use the colon commands ":order" and ":function" commands, if needed, to change the type of fit to be used. My experience is that :order 3 is usually all that is needed with any of the polynomial fitting functions.

Initial tracing for the science target spectrum

Final tracing for the science target spectrum

For an extended source the tracing may be less accurate than that for a calibration star, depending on whether there is any spatial variation in the spectrum over the slit.

Aperture for the calibration star spectrum

Original tracing for the calibration star spectrum

Linear tracing for the calibration star spectrum

Revised tracing for the calibration star spectrum

Now the script produces wavelength calibrated raw spectra for the science target and for the calibration star, still rather badly chopped up by the atmospheric transmission.

Wavelength calibrated science target raw spectrum

Wavelength calibrated calibration star raw spectrum

Ones science spectrum may look quite different than what is shown above, but if your telluric calibration spectrum looks much different than what is shown just above then something is wrong.

Telluric Correction

The most difficult part of the reduction is the telluric correction, and the solution presented here may not be the optimum result that is possible. The atmospheric bands are deep so the corrections that are to be made are large for certain parts of the atmospheric window, and this amplifies the noise in the two spectra.

The method used here is to set dscale and offset to 0.0 with colon commands

:dscale 0.0
:offset 0.0

which is done with your cursor in the plot window. Generally you then have to rescale the upper plot with the "w" and "e" key commands to see the useful part of the spectrum, as there are usually large spikes caused by division of essentially zero by zero in the regions where the atmospheric bands are totally opaque. It usually takes several attempts before the scaling in the upper window is set so that you can see the spectrum.

You should also use the "t" and "s" commands in the lower window to select the regions where there is signal and to exclude the regions just below 24 microns and near 25.5 microns where there is no signal due to strong water bands. This prevents these bad regions from dominating the RMS that is calculated across the spectrum. Once this is done, you then manually enter scale and shift values to get the smoothest possible spectrum, perhaps using the "a" key command to get the automatic "best shift" for a given scale value. There is no unique way to do this that I know of. My preference is keep the scaling value near 1.0, but in some instances as here a rather different value is the only way to produce a good result.

The telluric corrections are done using the msabsflux routine which in turn calls the telluric task which is in the noao.onedspec package in IRAF. The help for the telluric task should be consulted concerning the various key commands and colon commands available.

The original plot from telluric via mstelluric

The plot from telluric after dscale and offset set to 0

The plot from telluric after the shift is reset

The plot from telluric after the wavelength regions are selected using the "t" and "s" keys

The plot from telluric for the "best" shift and scaling

At present finding the best shift and scaling is a matter of trial and error. You should seek to produce the smoothest possible output spectrum with the atmospheric bands removed. It is not always possible to completely remove the bands, but as shown here even in Q-band it is usually possible to find a combination of shift and scale values that give a good result.

The Final Spectrum

Even at the "best" shift the band centers show spikes, such as at 21.78 microns. Due to the two regions where there is no signal and the shall shift needed to best remove the atmospheric bands there are also larger spikes at the edges of the zero signal regions. These can be removed using the "j" key in the splot routine, followed by writing out the new spectrum as a fits file with the "i" key command..

I supplied the standard star name in the call to msreduce and set the fl_bbody flag to "no", so the resulting spectrum is (ideally) calibrated in Jy. The blackbody option produces a spectrum that is correct in shape but for which the level is arbitrary.

I now have a the calibrated spectrum for the science target HD 56126. If this were a reduction for science I would need to scale the values for the proper brightness of the object, or (better!) have the expected spectral values corrected for the atmosphere available in IRAF for Procyon and use those in the reduction. The spectrum shown is too bright by a factor of 33.97/18.45 = 1.84 with about a 7% uncertainty based upon the IRAS photometry for these two stars. In this case as only part of the object was sampled in taking the spectrum the absolute calibration is probably not important.

Plot of the re-scaled spectrum

I get a nice spectrum out that in fact matches very well in shape with the ISO spectrum of the object, and without problems across the strong bands. This is in part due to the object being bright, but it indicates what is possible for ground-based Q-band spectroscopy from Mauna Kea under good conditions. My final cleaned spectrum is shown below. The Y axis is in Jansky. The output units can be set with the outtype parameter in the call to msreduce as you wish.

Plot of the cleaned and re-scaled spectrum

A log of all the keystrokes I used to reduce the spectrum can be found here along with some comments as to what is being done. The complete IRAF terminal log is also available for people to look at.

Last update 2008 September 8, K. Volk & R. Mason