Basic CCD Reduction, TO training, part 2
In this exercise we will continue processing the spectroscopic data that
was processed through the preliminary reduction phase using CCDPROC
(including the response and illumination corrections).
This exercise will lead you through the extraction of the spectra to 1-d,
the wavelength calibration of the arc spectra, and then finally applying
this information to our extracted stellar spectra.
We will assume that you are logged into IRAF in an xgterm window.
Go to your IRAF home directory and then go to the subdirectory called spec
in your scratch directory like:
cl> cd
cl> cd /scratch/lpss1b/juerg/spec
The data that we will work with at this time are r62506-r62510 and
r62671-r62673.
Are the data there? Check to be sure the images have been reduced - do
you remember how to check that?
cl> dir
cl> imhead r*.imh
cl> imhead r62672 l+ # what do you want to look for here?
These are the images that we will be working with.
cl> imhead r*.imh
r62506.imh[1023,256][real]: f4370 1200b
r62507.imh[1023,256][real]: f4370 1200b
r62508.imh[1023,256][real]: f4370 1200b
r62509.imh[1023,256][real]: f4370 1200b
r62510.imh[1023,256][real]: f4370 1200b
r62671.imh[1023,256][real]: ar4370
r62672.imh[1023,256][real]: RXJ2117 4370A 1200
r62673.imh[1023,256][real]: ar4370
EXTRACTING THE DATA
The first thing we need to do is extract the data to 1-d spectra, both
the objects and the arcs. If we wanted to leave the data as 2-d data
then the problem would be approached in a different way, using tasks in
the NOAO.TWODSPEC.LONGSLIT package. Let's load the appropriate packages.
cl> imred
im> kpnoslit
Our first task is determine the parameters for the task APALL, the task we
will use to do the extraction of the spectra. We need to determine a
profile centering width that is the number of pixels used to determine the
center of the profile during the centering and tracing of the profile along
the dispersion direction; this is usually set to the FWFM of the profile.
kp> display r62672 1
kp> imexam # measure the FWFM of the profile using
the "K" key
[when you have decided on a good value of the profile radius quit
imexam with q]
q
The next step is to edit the parameters for APALL until they resemble the
parameter file below. There are many ways to vary the extraction algorithms
to get the best result for your data. Our input images will be r62672.
We will use many of the defaults except for a few: set b_order=2 (a straight
line fit for the background subtraction rather than a constant),
set b_naver=-3 (that way a median value will be determined for each of our
background samples before they are fit with the straight line), type Ctrl/N
to go to the next page of the parameter list, set the width=5 (at least that
is what I decided from our runs with imexam earlier), set nfind=1 (we only
want to initially find 1 peak), and set background=fit.
Be aware that the header keyword DISPAXIS with a value of 1 must be present
for all files you process with apall, see use of "response" in the previous
exercise.
kp> unlearn apall
kp> epar apall # set parameters
:q # save the parameters
kp> lpar apall
input = "r62672" List of input images
nfind = 1 Number of apertures to be found automatically
(output = "e62672") List of output spectra
(format = "onedspec") Extracted spectra format
(references = "") List of aperture reference images
(profiles = "") List of aperture profile images\n
(interactive = yes) Run task interactively?
(find = yes) Find apertures?
(recenter = yes) Recenter apertures?
(resize = yes) Resize apertures?
(edit = yes) Edit apertures?
(trace = yes) Trace apertures?
(fittrace = yes) Fit the traced points interactively?
(extract = yes) Extract spectra?
(extras = yes) Extract sky, sigma, etc.?
(review = yes) Review extractions?\n
(line = INDEF) Dispersion line
(nsum = 10) Number of dispersion lines to sum or median\n\n
(lower = -5.) Lower aperture limit relative to center
(upper = 5.) Upper aperture limit relative to center
(apidtable = "") Aperture ID table (optional)\n\n# DEFAULT BACKG
(b_function = "chebyshev") Background function
(b_order = 2) Background function order
(b_sample = "-10:-6,6:10") Background sample regions
(b_naverage = -3) Background average or median
(b_niterate = 0) Background rejection iterations
(b_low_reject = 3.) Background lower rejection sigma
(b_high_rejec = 3.) Background upper rejection sigma
(b_grow = 0.) Background rejection growing radius\n\n# APERTU
(width = 5.) Profile centering width
(radius = 10.) Profile centering radius
(threshold = 0.) Detection threshold for profile centering\n\n#
(minsep = 5.) Minimum separation between spectra
(maxsep = 1000.) Maximum separation between spectra
(order = "increasing") Order of apertures\n\n# RECENTERING PARAMETERS\n
(apertures = "") Select apertures
(npeaks = INDEF) Select brightest peaks
(shift = yes) Use average shift instead of recentering?\n\n#
(llimit = INDEF) Lower aperture limit relative to center
(ulimit = INDEF) Upper aperture limit relative to center
(ylevel = 0.1) Fraction of peak or intensity for automatic wid
(peak = yes) Is ylevel a fraction of the peak?
(bkg = yes) Subtract background in automatic width?
(r_grow = 0.) Grow limits by this factor
(avglimits = no) Average limits over all apertures?\n\n# TRACING
(t_nsum = 10) Number of dispersion lines to sum
(t_step = 10) Tracing step
(t_nlost = 3) Number of consecutive times profile is lost bef
(t_function = "legendre") Trace fitting function
(t_order = 2) Trace fitting function order
(t_sample = "*") Trace sample regions
(t_naverage = 1) Trace average or median
(t_niterate = 0) Trace rejection iterations
(t_low_reject = 3.) Trace lower rejection sigma
(t_high_rejec = 3.) Trace upper rejection sigma
(t_grow = 0.) Trace rejection growing radius\n\n# EXTRACTION
(background = "fit") Background to subtract
(skybox = 1) Box car smoothing length for sky
(weights = "none") Extraction weights (none|variance)
(pfit = "fit1d") Profile fitting type (fit1d|fit2d)
(clean = no) Detect and replace bad pixels?
(saturation = INDEF) Saturation level
(readnoise = "0.") Read out noise sigma (photons)
(gain = "1.") Photon gain (photons/data number)
(lsigma = 4.) Lower rejection threshold
(usigma = 4.) Upper rejection threshold
(nsubaps = 1) Number of subapertures per aperture
(mode = "ql")
Now we are ready to extract our stellar spectra.
kp> apall
[answer yes to all questions]
[we are now presented with a profile plot for our first star - the
center and the aperture have been marked]
? # for a list of cursor options
b # go into background fitting mode
[the default backgrounds are marked but they are not very good - we
should move then further away from the star]
? # for a list of cursor options in this mode
t # initialize the background regions
s # mark the left boundary of the first new
background region
s # mark the right boundary of the first new
background region
s # mark the left boundary of the second new
background region
s # mark the right boundary of the second new
background region
f # fit the new background
:sample # print out the sky regions being used
[you should now have one background sample region on either side
of the profile of about 6 pixels or so each - notice the fit
extended under the star - also notice the median points in each
sample region - redo the procedure above if you are not happy with
your sky background]
q # quit the background fitting routine
[once we have decided that the aperture size, the center and the
background are ok, we can move onto the tracing routine. The next
step is to trace the profile of the star along the dispersion axis.
We are again in the ICFIT routine (remember the RESPONSE and
ILLUMINATION tasks?]
q
[answer yes to all questions - you will now be presented with a plot
of the profile centers versus the line number. Notice the fitting
parameters at the top of the plot. We can adjust the fit if you
would like, although the current fit looks quite good.]
:func spline3
f
:order 3
f
:niter 1
f
q # quit tracing mode and move onto extraction
[answer yes to all questions - now you will see a plot of the
extracted spectrum - notice the cosmic rays - another option
during the extraction process is to clean the spectrum but we
did not choose that option]
q # go to next spectrum if list has been defined
[having defined a list of objects you will now step through the
remaining stars - each will be processed in a similar manner -
note the previous background regions have not
been remembered so you will need to set these for each star - the
previous fitting parameters for the tracing are remembered however]
Our stellar spectra has now been extracted to 1-d images. The headers
have been updated with a more complex set of header keywords as well. And
all of the extraction information has been saved in individual files in
a database directory.
kp> imhead e*.imh # list the new extracted spectra
kp> imhead e62672.0001.imh l+ | page # see the long header
kp> dir database
kp> page database/apr62672
Now we need to extract the two arc spectra to 1-d. But this process is done
differently than for the star. We can not trace the arcs since there is no
continuum, and we need to be sure we do not do any background subtraction or
we will end up with nothing in our output spectrum.
The arclines may as well be tilted, thus we have to extract the arc along
the same fitting curve as the object. So we will use a "reference" image to
extract the arcs. In other words, we will use the database file created earlier
when we extracted our stars.
We use the task APALL again to extract the arc spectra without editing the
parameters.
kp> apall r62671 out=e62671 ref=r62672.imh recen- trace- back- intera-
kp> apall r62673 out=e62673 ref=r62672.imh recen- trace- back- intera-
kp> imhead e*.imh
kp> implot e62671.0001 # check extractions
:i e62673.0001
l
q
WAVELENGTH CALIBRATION
We are now ready to determine a wavelength calibration for our arc lamps.
A chart of the lines in the spectrum is in the file spec.eps in the docs
directory. We will use the tasks IDENTIFY and REIDENTIFY. You will need to
first determine the width of the line profiles, much like we did for the
stellar profiles earlier.
kp> unlearn identify reidentify
kp> implot e62671.0001
kp> epar identify # set parameters
:q # save the parameters
Now we are ready to run the task IDENTIFY. The default profile radius
may work fine for these data (fwidth=4). The only parameter that we
need to modify is the image name.
When we exit IDENTIFY the wavelength solution will be saved in a database.
kp> lpar identify
images = "e62671.0001" Images containing features to be identified
(section = "middle line") Section to apply to two dimensional images
(database = "database") Database in which to record feature data
(coordlist = "linelists$idhenear.dat") User coordinate list
(nsum = "10") Number of lines/columns/bands to sum in 2D images
(match = 10.) Coordinate list matching limit in user units
(maxfeatures = 50) Maximum number of features for automatic identification
(zwidth = 100.) Zoom graph width in user units
(ftype = "emission") Feature type
(fwidth = 4.) Feature width in pixels
(cradius = 5.) Centering radius in pixels
(threshold = 10.) Feature threshold for centering
(minsep = 2.) Minimum pixel separation
(function = "spline3") Coordinate function
(order = 1) Order of coordinate function
(sample = "*") Coordinate sample regions
(niterate = 0) Rejection iterations
(low_reject = 3.) Lower rejection sigma
(high_reject = 3.) Upper rejection sigma
(grow = 0.) Rejection growing radius
(autowrite = no) Automatically write to database
(graphics = "stdgraph") Graphics output device
(cursor = "") Graphics cursor input
(mode = "ql")
kp> identify e62671.0001
? # review cursor options for IDENTIFY
w # enter the window cursor options
? # look at the options
w
f # flip the data so increasing wavelength
runs to the right
[mark 5 strong features on the plot spanning the spectrum - point the
cursor, type "m", then type in the wavelength. If the terminal
beeps at you it can not find the line center - on a crowded spectrum
it may be necessary to window the plot around the line
you wish to identify - use "w" followed by "e" at the lower left
corner of the box you wish to expand and then "e" at the top right
edge of the box - type "w" followed by "a" to go back to the full
plot]
f # fit the lines and enter ICFIT
? # look at cursor options
:func cheby # change function
:order 3 # set to quadratic
f # refit
q # go back to IDENTIFY
w # to get spectrum back if it is not
a there
l # to find other lines from thorium.dat
file
f # fit the lines and enter ICFIT
:nite 2 # iterate the solution to remove bad
lines
f # refit
d # delete a few other outliers by pointing
the cursor and typing "d"
f # redo fit
q # go back to IDENTIFY
d # delete other poor lines - we should
be able to get a solution with an
RMS=~0.01 A in ICFIT
f # redo fit
l # look at the fit of the non-linear
components
q # go back to IDENTIFY
q # quit IDENTIFY when you are happy with
the fit - be sure to save it!
kp> dir database # the solution is saved in the
ide62671.0001 file
kp> page database/ide62671.0001
Now we can run REIDENTIFY on our other arc spectrum using e62671.0001 as our
reference image.
kp> reidentify e62671.0001 e62673.0001 coord=linelists$idhenear.dat v+ inter+
[or do epar followed by a :go - you will have the option of reviewing
the fit with IDENTIFY]
kp> imhead e62671.0001 l+
Notice the keyword refspec1 set to the name of the image itself. This keyword
was added by IDENTIFY/REIDENTIFY and defines the arcs as "reference images".
The next task, REFSPECTRA, assigns "reference images" to our object frames,
allowing us to select various ways for determining these assignments. For
this data set we should assign the arcs by the julian date at the middle of the
exposure.
Now edit the parameters for REFSPECTRA so they look like those below.
kp> epar refspectra # set parameters
:q # save the parameters
kp> lpar refspectra
input = "e62672.0001" List of input spectra
answer = "yes" Accept assignment?
(references = "e62671.0001") List of reference spectra
(apertures = "") Input aperture selection list
(refaps = "") Reference aperture selection list
(ignoreaps = yes) Ignore input and reference apertures?
(select = "interp") Selection method for reference spectra
(sort = "jd") Sort key
(group = "rotskypa") Group key
(time = no) Is sort key a time?
(timewrap = 17.) Time wrap point for time sorting
(override = no) Override previous assignments?
(confirm = yes) Confirm reference spectrum assignments?
(assign = yes) Assign the reference spectra to the input spect
(logfiles = "STDOUT,logfile") List of logfiles
(verbose = no) Verbose log output?
(mode = "ql")
kp> refspectra # answer "yes" or "YES" - what's the difference?
kp> imhead e62672.0001 l+ # note the new keywords
The task that actually applies the wavelength calibration to your data is the
task DISPCOR. It uses the refspec keyword information in the header to
determine what arcs to use to apply the solution.
Edit the parameters for DISPCOR so they resemble those below. We will linearize
the data, which means the data is rebinned so each pixel has the same delta
lambda; if we chose not to linearize then the data values would be unaltered
and the coefficients for the solution would be stored in the image header.
kp> epar dispcor # set parameters
:q # save the parameters
kp> lpar dispcor
input = "e62672.0001" List of input spectra
output = "o62672" List of output spectra
(linearize = yes) Linearize (interpolate) spectra?
(database = "database") Dispersion solution database
(table = "") Wavelength table for apertures
(w1 = INDEF) Starting wavelength
(w2 = INDEF) Ending wavelength
(dw = INDEF) Wavelength interval per pixel
(nw = INDEF) Number of output pixels
(log = no) Logarithmic wavelength scale?
(flux = yes) Conserve flux?
(samedisp = no) Same dispersion in all apertures?
(global = no) Apply global defaults?
(ignoreaps = no) Ignore apertures?
(confirm = no) Confirm dispersion coordinates?
(listonly = no) List the dispersion coordinates only?
(verbose = yes) Print linear dispersion assignments?
(logfile = "") Log file
(mode = "ql")
kp> dispcor
kp> imhead o62672 l+
Congratulations! You have done it! Now you can play with this data.
kp> splot o62672 # welcome to the world of SPLOT - now
that you are experienced IRAF users
check out this very versatile analysis
tool!
kp> specplot o62672
Note that SPLOT has cursor keys for cleaning bad pixels as well. Look at
the line measuring cursor keys in SPLOT too. Experiment!
Our last step in this exercise is to normalize the spectrum to 1 by fitting
a function to the continuum. The task used is called continuum.
The trick here is to have a sufficiently flexible function but not to fit
the true spectral features. Again Experiment!
kp> epar continuum # set parameters
:q # save the parameters
kp> lpar continuum
input = "o62672" Input images
output = "oo62672" Output images
ask = "yes"
(lines = "*") Image lines to be fit
(bands = "1") Image bands to be fit
(type = "ratio") Type of output
(replace = no) Replace rejected points by fit?
(wavescale = yes) Scale the X axis with wavelength?
(logscale = no) Take the log (base 10) of both axes?
(override = no) Override previously fit lines?
(listonly = no) List fit but don't modify any images?
(logfiles = "logfile") List of log files
(interactive = yes) Set fitting parameters interactively?
(sample = "*") Sample points to use in fit
(naverage = 1) Number of points in sample averaging
(function = "spline3") Fitting function
(order = 3) Order of fitting function
(low_reject = 2.) Low rejection in sigma of fit
(high_reject = 2.) High rejection in sigma of fit
(niterate = 10) Number of rejection iterations
(grow = 1.) Rejection growing radius in pixels
(markrej = yes) Mark rejected points?
(graphics = "stdgraph") Graphics output device
(cursor = "") Graphics cursor input
(mode = "ql")
kp> continuum
-------------------------------end of exercise---------------------------------