F2 Longslit HK Reduction - When Things Go Wrong

An example of when things go wrong and what can be done about it

WARNING: This example shows a few tricks to try and to look for when things go wrong. There are a few problems with this dataset. In particular, the telluric is far from optimal (wrong airmass, taken much later after the observation). The wavelength calibration also shows problem.

This is a study example, not even close to a smooth, normal reduction.

Science Dataset:

Telluric Dataset:

WARNING: This is the wrong telluric. Do not use for science. But it is a heck of a good training set for when things go wrong. <br >

File numbers:

STEP 0: Define variables specific to this target and start parallel PyRAF session.

Now get a PyRAF session ready. This will be needed for some interactive tasks that don't work well in the notebook.

Open an xtern or a Terminal, then:

cd '/Volumes/Rugged2/GS-2013B-Q-73/SDSSJ005918.23+002519.7/test_notebook/reduxHK'
pyraf

Then in the PyRAF session:

gemini
f2
unlearn gemini
unlearn f2
unlearn gnirs
unlearn gemtools

Get your PyRAF configured. In the PyRAF session:

iraf.f2.logfile = "HK005918.log"
iraf.f2.database = "HK005918_database/"
set rawdir = "../../../raw/"
set stdimage=imt2048
nsheaders('f2', logfile=iraf.f2.logfile)

STEP 1: Initialize the required packages

Launch ds9 before proceeding further.

Load the packages required for the notebook session.

Reset tasks to the default parameters. (Note: this doesn't seem to be working from the Python shell.)

STEP 2: Define various parameters and variables.

If necessary to start from scratch, delete the log file and the database.

STEP 3: Create the reduction lists

STEP 4: Visually inspect the data

STEP 5: f2prepare all the data

STEP 6: Create the necessary dark frames

STEP 7: Create the normalized flat field and BPM.

Subtract the dark from the flat images then cut

Construct the normalised flat field.

The flats are derived from images taken with the calibration unit (GCAL) shutter open ("lamps-on"). It is recommended to run nsflat interactively to ensure a good fit of the lamp's profile.

The fit order must be high enough to fit the lamp's profile but not too high since we do not want to fit the variations due to the pixel to pixel variation. For the HK filter, the fit should look like in the screenshots below. Note that in the zoomed plot, the fit matches the overall profile, not the high frequency variations.

The overall fit:

<br > Zooming in to confirm that the high frequency signal is not being fit:

If the order needs to be changed, the :order command will be required and that does not work well in the iPython notebook. So copy the following to the PyRAF we prepared earlier, and find a good fit.

iraf.imdelete('flat.fits,f2_ls_bpm.pl', verify='no')
iraf.nsflat('cdf@flat.lis', flatfile='flat.fits', \
    bpmfile='f2_ls_bpm.pl', thr_flo=0.35, thr_fup=3.0, \
    fl_inter='yes', order=120)

STEP 8: Reduce the arc and determine the wavelength solution

The wavelength solution is derived from an arc taken right after the science observation, before anything is moved (eg. grating, slit, etc.)

If we are lucky, the arc will also work for the telluric observation. However, as we will see later for this target, when things in the optical path are moved, the arc might not be work at all.

If the optical path moved such that the science arc solution cannot be applied, and no arcs were taken with the telluric, one might be able to use OH sky lines in the telluric. The exposure would have to be long enough to get good signal on the OH lines. The coverage of the OH lines is not as complete as the arc, but it might be the only solution.

In any case, we will reduce the arc and calculate the wavelength solution as this is what we surely be using for the science observations.

The wavelength solution step really should be done in interactive mode as it is not uncommon to have to reject lines and/or adjust the matches with the line list.

Reduce the arc

Determine the wavelength solution from the arc

Normally, the lines will be identified automatically and this step, while interactive is used to delete a couple outliers and visually confirm that the fit is good.

If, for some reason (it can happen), the lines do not get identified automatically and they need to be marked manually, you will need to run the two commands below in the PyRAF sessions since you will need full interactive support. But most of the time, running this in the notebook works okay.

Test recovery of a calibrated arc

This step is to confirm that if we apply the wavelength solution to the arc itself, we do recover a correctly wavelength-calibrated arc. The key step is nsfitcoords. The default lxorder and lyorder are both set to 2 which doesn't seem to be correct. Running nsfitcoords interactively on the arc will allow us to discover the best lxorder and lyorder for this wavelength solution, and use that as a starting point later when using the arc wavelength solution (warc.fits) on the telluric and the science.

Interactive IRAF tasks don't work very well from the iPython Notebook. This is a case where we need to use PyRAF directly. Cut and paste these instructions in the PyRAF session that we prepared in Step 0.

iraf.imdelete('farc.fits', verify='no')
iraf.nsfitcoords('arc.fits', lamptransf='warc.fits', fl_inter='yes')
# Note down the xorder and yorder needed for a good fit.  Use
# that whenever warc.fits is used by nsfitcoords.
# :xorder 3, :yorder 6

The default xorder and yorder are 2 and 2. Below is the residual plot one gets with those defaults: the residuals are large and show a pattern, clearly indicating that the fit is not good at all.

<br > Adjusting the xorder to 3, tightens the x-residuals.

<br > Finally, adjusting the yorder to 6, reduce the y-residuals within acceptable limits and remove the big wave pattern.

Record the lxorder and lyorder in the notebook session for future use.

Now continue with the test.

Check that the arc lines are at the correct wavelength. Use 'k' on each side of the a line to get the centroid. Then compare with the line list and the diagrams.

Resource:
http://www.gemini.edu/sciops/instruments/niri/?q=node/10166

The NIRI plots from that webpage work well.

Result: (Kathleen) The arc gets transformed correctly with lxorder=3, lyorder=6. The lines are within 1-2 Angstroms of the expected values. Note that the shape of the arc lines is not symmetric. No idea why, but it surely affects centering. Bottom line is that with the proper lxorder and lyorder, the wavelength solution leads to a good calibration.

STEP 9: Reduce the telluric standard dataset

STEP 10: Combine the reduced telluric datasets

STEP 11: Wavelength calibrate the telluric data

As it turns out, this step can be tricky. Theorically, the arc taken with the science should be fine if the telluric is taken just after the science and the optical path is identical. Unfortunately, we have found that the optical path can vary, leading to the telluric being offset from the science. This effect was indeed observed for this target. Below, we document the whole process that led to the discovery of the problem and to the solution. Note that from now on (circa Sep 2014), an arc will be observed with the telluric; Ruben Diaz has modified the OT library accordingly.

We first (Option 1) use the arc that was taken with the science and inspect the results. Then (Option 2)we look at the OH sky lines to investigate the problem and measure the offset. That offset is applied to wavelength solution and that translated solution is applied to the telluric. We find that the results are good. Finally (Option 3), since the sky lines are quite strong, we use the sky lines to calculate a new wavelength solution and use that solution on the telluric. We find that the results are very good, and that will be the solution used for science.

The key tasks used for this step are:

• (nswavelength - Option 3 only: Used to calculate a wavelength solution from the sky lines.)
• nsfitcoords: Used to determine the final 2-D surface solution, based on a previously calculated wavelength solution, to be applied to the data.
• nstransform: Actually applies the 2-D solution found by nsfitcoords to the data.

Note that spatial rectification (s-distortion correction) is not needed with F2 longslit data.

Option 1: Use the arc taken with the science.

Check that the spectral features (the hydrogen lines in this case) fall where they are supposed to be.

Cut an paste these instructions in the PyRAF session. (The splot interactive ':'-commands do not work well from the notebook.)

iraf.display ('tftel_comb[sci,1]', 1)
# note on which column the telluric falls
iraf.tel_position = 928

splot ("tftel_comb[sci,1]", line=iraf.tel_position)
# :nsum 30
# Use 'k' on each side of the a line to centroid.
# Compare with atonic line list.

In this case, one finds that the measured wavelengths of the hydrogen lines in the corrected telluric do NOT match the wavelengths in atomic line list. The wavelength solution calculated from the arc taken with the science observation CANNOT be used on the telluric.

The investigation of what's going on is found in Option 2, below.

Option 2: Use the sky lines to measure any offsets and apply offset to arc solution.

Measure the offset.

One telluric frame has strong enough sky lines. Therefore, we will reduce that frame on its own, but without sky subtraction.

Apply the science arc solution to the non-sky subtracted telluric frame.

Use splot to measure the wavelength of the sky lines in the corrected telluric. Do this in the PyRAF session, not in the notebook.

iraf.splot ('tfXdf'+first_tel+'[sci,1]', line=750)
# line=750 is just somewhere in the middle.
# :nsum 30

Check that the OH sky lines are at the correct wavelength. Use 'k' on each side of the a line to get the centroid. Then compare with the plots and line lists.

Resources:
http://www.jach.hawaii.edu/UKIRT/astronomy/calib/spec_cal/oh.html http://www.jach.hawaii.edu/UKIRT/astronomy/calib/spec_cal/ohlines.html

What we find here is that there is a constant offset in wavelength (through H and K bands) of $+11\pm0.5\ \unicode[serif]{xC5}$ compared to the OH line list. Something in the optical path moved between the arc and the telluric observation.

Offset the arc wavelength solution and apply.

The idea now is to use the wavelength solution obtained from the arc, but shift it by $11\ \unicode[serif]{xC5}$ before applying it to the telluric.

Edit the fitcoords file to apply the offset to the solution.

1. First go to the very last solution in the file.
2. Change the name of the solution, on the begin line to ftel_comb_offset_SCI_1_lamp.
3. Then, subtract $11\ \unicode[serif]{xC5}$ from the 9th line of the last solution in the file (start counting after the surface record). To make sure you got it, see the line identified in this screenshot:

Then apply that transformation to the combined telluric.

Use splot to measure the wavelength of the stellar features in the corrected telluric. Do this in the PyRAF session, not in the notebook.

iraf.splot ('tftel_comb[sci,1]', line=iraf.tel_position)
# line=iraf.tel_position is where the telluric falls. It should 
#        still be defined from previous investigation steps.
#        If not, just display the image and identify the column
#        the telluric falls on.
# :nsum 30

Option 2 is obviously less then ideal but it gives decent result. However, in our case, we have nice, strong OH sky lines in our telluric observations and we will use those to calculate the wavelength solution rather than using a translated arc solution. Moving to Option 3 then.

Option 3: Use the sky lines to calculate the wavelength solution.

Reduce the first telluric, keeping the OH sky lines. Note that if the OH sky lines were a bit faint, one could stack all the tellurics together, without shifts, and increase the signal-to-noise on the sky lines, then proceed with that image instead.

IMPORTANT: There are no OH sky lines redwards of $22750\ \unicode[serif]{xC5}$ and blueward of $15000\ \unicode[serif]{xC5}$ which means that the wavelength calibration in those regions will have reduced accuracy. In our case, that's fine, we don't have features of interest in that part of the spectrum for this source. If the lack of lines were a problem, one might want to consider Option 2 instead.

Determine the wavelength solution for the telluric using the OH sky lines.

IMPORTANT: Run this interactively and do delete the bad lines from the fit and re-fit. There are bad line identifications redward of $22750\ \unicode[serif]{xC5}$ and blueward of $15000\ \unicode[serif]{xC5}$. Remove all those bad lines.

Use that new wavelength solution on the telluric.

Most likely the lxorder and lyorder should be same as for the arc solution given that this must be a function of the optics. I (Kathleen) did run nsfitcoords interactively in a PyRAF session and I confirm that lxorder 3 and lyorder 6 are required even with the OH line solution. No need to run it interactively anymore.

Verify that the solution works well before proceeding further. In the PyRAF session:

iraf.splot('tftel_comb.fits[sci,1]', line=iraf.tel_position)
# :nsum 30
# Use 'k' to centroid the stellar feature, 
# ie. the hydrogen lines to make sure
# that they fall where they are supposed to.

STEP 12: Extract the telluric spectrum

We run nsextract in non-interactive mode since the telluric is quite bright. (Running it in interactive mode does confirm that the automatic aperture selection is coincident with the telluric spectrum.)

This telluric spectrum will be used later. It will then need to be corrected to remove the stellar features either by fitting Voight profiles on each hydrogen line and removing each feature, or by fitting a stellar model to it.

STEP 13: Reduce the science data

STEP 14: Combine the science data

The rejtype is set to minmax in the example included in the Gemini IRAF package. Kathleen finds that for faint targets, that does not work at all. Therefore, here we set rejtype to none which seems to work well.

STEP 15: Wavelength calibrate the science data

It is assumed that the arc has been taken right after (or before) the science observations, without anything moving in-between, other than getting the arc in the field of view.

While the OH sky lines could be used, the coverage is not always adequate across the dispersion axis and the sky lines tend to be a bit broader than the arc lines. True, we used the OH sky lines for the telluric above, but that was because we could not use the arc with greater accuracy.

Yet, we will nevertheless double-check that the arc works well on the science observation. To that end, we will use a science frame that has not been sky subtracted yet, and apply the arc solution. We will measure the OH sky lines to confirm that they fall where they are supposed to fall. If this test passes, then we will apply the arc solution to the combined science frame from Step 14.

(It should pass. If it doesn't, something is really wrong with the data and one might have to use the OH sky lines for the wavelength calibration after all.)

The key tasks used for this step are:

• nsfitcoords: Used to determine the final 2-D surface solution, based on a previously calculated wavelength solution, to be applied to the data.
• nstransform: Actually applies the 2-D solution found by nsfitcoords to the data.

Spatial rectification (s-distortion correction) is not usually needed with F2 longslit data.

Test accuracy of the wavelength solution the sky lines.

We reduce the first science frame but without the sky subtraction.

We apply the arc solution on the frame with the sky lines.

Use splot to measure the wavelength of the sky lines in the corrected telluric. Do this in the pyraf session, not in the notebook.

first_obj = iraf.head('obj.lis', nlines=1, Stdout=1)[0].strip()
iraf.splot ('tfXdf'+first_obj+'[sci,1]', line=750)
# line=750 is just somewhere in the middle.
# :nsum 30
# Use 'k' again and sky line plots and lists to check position
#  of the sky lines.

Resources:
http://www.jach.hawaii.edu/UKIRT/astronomy/calib/spec_cal/oh.html
http://www.jach.hawaii.edu/UKIRT/astronomy/calib/spec_cal/ohlines.html

Verifications of the OH lines in the H-band section and the K-band section show that the lines do show up at the correct wavelength after the transformation. We can use the arc.

Apply the wavelength solution from the arc to the combined science frame.

STEP 16: Extract the science spectrum

For faint science spectra, the telluric can be used as a reference when extracting the science spectrum; set trace=tftel_comb.fits.

However, in this case, the science spectrum is strong enough to be traced.

For this extraction, one needs to:

1. delete, 'd', the aperture automatically selected. It selects some random spike at the far left of the cross-section.
2. mark, 'm', the real location of the science spectrum, around column 950.
3. adjust with 'u' and 'l', the width of the aperture, it seems a bit too wide.

STEP 17: Apply the telluric correction to the science spectrum

No matter what standard star you observed for the telluric correction, there will be stellar features in the spectrum. There are two ways to remove the stellar features from the telluric spectrum: 1) fit a line profile to and remove each stellar feature from the telluric spectrum; 2) use a stellar model, resample, scale, and subtract it from the telluric to remove the stellar features.

The stellar model option is possibly better, but whether it is necessary is still unclear (circa late-Sep 2014). The question is whether we need to fit a blackbody to the telluric to properly calibrate the science spectrum. (Depends a bit on the science, I (KL) think). nstelluric does fit the continuum of the telluric and of the science to scale the standard. Is this sufficient? TBD.

Note: The telluric star is a B9V star. (We have not been able to find a B9V model with lines in the near-IR but a A0V model should do just fine.)

For now, we use rmfeatureto fit a Voight profile to each hydrogen feature and remove it from the telluric spectrum. The procedure is a bit tedious because of the very clunky nature of rmfeature but it does work in the end, with a bit of patience.

Remove stellar features from telluric spectrum

The stellar features we have to deal with are hydrogen absorption lines. There's only one in the K-band, but there are quite a few in the H-band. We start from the red-end and make our way to the bluest hydrogen we can identify and fit.

There are two approaches to deal with the stellar features:

1. manual removal by fitting Voight profiles to the features
2. use of a stellar model.

Right now, since we do not have a procedure to resample and scale the stellar model, we will remove the features manually with a custom-made python scripts rmfeature. Note that there are tools in IRAF to do that type of removal, and also the deal with the stellar model, but I (Kathleen) have no idea what those are or how to use them.

Manual removal of features with Voight profile

The stellar features are (in $\unicode[serif]{xC5}$) : 21655[1], 17362[2], 16807[3], 16407[4], 16109[5], 15881[6], 15701[7], 15557[clean].

We are not correcting for 15192, 15261, 15342, 15439 as they are really lost in the noise/telluric absorption.

Creation of a blackbody

The telluric star is HIP 1378. It is a B9V star with an approximate temperature of $10,700 K$.

Resource for temperatures:
http://www.jach.hawaii.edu/UKIRT/astronomy/utils/temp.html

First, we figure out the start and end wavelength, and the sampling we require to match the science spectrum. Then we create the blackbody spectrum and convert it into a MEF file. This will allow us to keep propagating the variance plane.

Use the clean telluric to correct the science spectrum

We will try two approaches:

1. nstelluric approach
2. manual approach

Spoiler warning: The nstelluric approach gives significantly better results (less noisy and more accurate continuum shape). See the Discussion section below, after the Manual approach section.

The nstelluric approach

nstelluric normalize the telluric and the science, and then cross-correlate the two to find first the shift in pixel between the two, and then the necessary scaling factor. Part of the scaling factor applied to the telluric is the airmass ratio. Then, the shifted and scaled normalized telluric is divided from the non-normalized science.

Procedure
Fit the continuum of the telluric and the science to scale the two to each other. Use only the H-band and K-band regions, ignore the regions where the atmosphere is opaque. Use the s option in nstelluric to define the regions to use (ie. the "sample"), or :sample to enter the regions precisely.

Run this interactive task in a PyRAF session. See screenshots below to get an idea of what it should look like.

iraf.imdelete('axtfobj_comb.fits', verify='no')
iraf.nstelluric('xtfobj_comb.fits', 'tel_clean', fitorder=12, \
    threshold=0.01, fl_inter='yes')

The order used to fit the telluric is: 4
The order used to fit the science is: 4 # same as for the telluric.
Scale: 1
Shift: -1.94
Sample: 19:62 311:643 1054:1540

Given that the arc solution is different from the science and for the telluric, a slight shift is probably expected.

Fitting the telluric:
<br > Fitting the science target:
<br > The next window is for adjusting the cross-correlation solution, if needed. In this case, after a bit of exploration of the parameter space, we find that the automatic solution is the best. A good area to focus on is the $2\mu$ region. <br > Zooming on the $2\mu$ region. See how the strong telluric features gets cleanly canceled. <br > <br >

iraf.splot('axtfobj_comb.fits[SCI,1]', ymin=-200, ymax=1000)

iraf.specplot('xtfobj_comb.fits[sci,1],axtfobj_comb.fits[sci,1],\
    tel_clean.fits[sci,1]', fraction=0.05, yscale='yes', ymin=-100,\
    ymax=4000)

The telluric spectrum, ie. the spectrum of a star, was divided into the science spectrum. This means that to recover the correct continuum for the science spectrum one needs to multiply the science spectrum with a blackbody (if rmfeature was used), or a stellar model (if rmfeature was not used.)

During the extraction, some 2D WCS are not being cleaned up correctly. Here is how to remove the 2D info from what is now a 1D spectrum.

Using Kathleen's Python splot in a shell, plot the final spectrum with line identifcation and with an overplot of the error.

The manual approach

Here we do not do any shifting or scaling or airmass correction. We divide by the telluric standard and multiply by the blackbody.

The final plot in the manual method is a lot more noisy and the slope of the continuum is shallower. One can make the spectrum less noisy by applying the shift found in nstelluric, and the slope can be made steeper by applying the airmass correction nstelluric used. This will be left as an exercise... (I, Kathleen, did all this already to convince myself that the nstelluric results were correct.)

Discussion on telluric correction techniques

Doing the telluric correction manually is suited for cases where the telluric standard was observed very close in time with the science to avoid time-dependent sky variations and at the same airmass, and obviously where nothing in the optical path moved and flexure is minimal.

In the present case, the optics moved, the telluric was observed 50 min after the science target, and at an airmass of 1.848 when the science was observed at airmass 1.465. Clearly, this is a messy observation.

The movements in the optics (or flexure) causes the telluric to be offset from the science spectrum as clearly shown during the wavelength calibration. This was also picked up by nstelluric which found that a shift of -1.94 pixel is necessary.

The large airmass difference is corrected for in nstelluric. But in the manual approach, that difference is not taken into account (one could, but we didn't above). The result was that the non-airmass corrected spectrum showed a much shallower continuum slope.

At least, it seems that the atmosphere was fairly stable since even if the data was taken 50 minutes later, nstelluric didn't have to apply a scaling factor (other than the airmass correction factor).

What was found here is that the nstelluric task, when used correctly does produce rather good results that are far less noisy than the results from the manual approach, especially when there are issues with telluric standard observations. With enough fiddling with shift, airmass correction, scaling, one should be able to reproduce the nstelluric result but it seems more straightforward with nstelluric. I (Kathleen) did play around and I got close, close enough to convince myself that the slope was okay, but never got the noise to go down to the nstelluric level, then I got bored.

FINAL spectrum: axtfobj_bb.fits