Welcome back in the sweet misery of data reduction
Great news to everyone - After finishing the last step of this post you should end up with perfectly reduced data (but I would like to highlight the word should because there will be a number of steps I am not certain at all).
Spectra extraction
In previous post, we got rid of the instrumental signature on our data but they are still part of image and not in the so-wanted form of spectra. Our first step here will be to extract spectra from the image. For this purpose, we are going to use apall task. In addition to extraction of our spectra it allows us to extract the background from spectra and to perform some additional extractions.
echelle> epar apall
_______________________________________________________
PACKAGE = echelle
TASK = apall
input = object List of input images
(output = object.ec) List of output spectra
(apertur= ) Apertures
(format = echelle) Extracted spectra format
(referen= orderdef) List of aperture reference images
(profile= ) List of aperture profile images
(interac= yes) Run task interactively?
(find = yes) Find apertures?
(recente= no) Recenter apertures?
(resize = no) Resize apertures?
(edit = yes) Edit apertures?
(trace = no) Trace apertures?
(fittrac= no) Fit the traced points interactively?
(extract= yes) Extract spectra?
(extras = yes) Extract sky, sigma, etc.?
(review = yes) Review extractions?
(line = INDEF) Dispersion line
(nsum = 10) Number of dispersion lines to sum or median
...
# EXTRACTION PARAMETERS
(backgro= median) 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?
(saturat= INDEF) Saturation level
(readnoi= outron) Read out noise sigma (photons)
(gain = outgain) Photon gain (photons/data number)
(lsigma = 4.) Lower rejection threshold
(usigma = 4.) Upper rejection threshold
(nsubaps= 1) Number of subapertures per aperture
(mode = ql)
Here are some parameters of this task. There is much more of them but they are focused on steps we have already done. We will run this task on scientific data and also on calibration data (LAMP,WAVE). We use ORDERDEF frame as reference frame for defining apertures. There are more options for background parameter and you should choose the one that best fits to you (none|average|median|minimum|fit). You should also check the parameters gain and readnoi.
Now, we can run the task:
echelle> apall object out=object.ec format=echelle referen=orderdef interac+ find+ recenter- resize- edit+ trace- fittrac- extract+ extras+ review+ backgro=median readnoi=outron gain=outgain
The window for aperture editing will be opened and you should check the background regions you defined previously. Do it by pressing 'b' with cursor over the given aperture.
Background regions should be located in the flat parts between orders. If they are not, you can define them for each aperture respectively. Bad regions can be deleted by pressing 't' and new regions can be set by pressing 's' at the left and right side of given region. Then press 'f' to fit them. If you are satisfied, press 'q'.
Press 'q' one more time if you do not want to do any more changes on the apertures. Answer 'yes' for all the question and you will open extracted spectrum from the first aperture. You can review all the apertures.
Do it also for the wave.fits
Extracted spectra can be checked by the splot task (something like implot but for spectra).
If you had done that one additional step from previous post, now is the time you will use it. You can compare extracted spectra with those reference spectra. If there is significant change (for example, there are missing a huge number of good-looking spectral lines, you should flatten your data again).
Wavelength calibration (prepare to be disgusted after several hours spent on this step)
This step is done in several tasks. At first, we need to calibrate the wave.fits. This is the spectrum of ThAr lamp. These spectra have well known spectral lines and we can find them by comparing extracted spectra with atlas' spectra (for UVES, these can be found here).
The first task we will use is ecidentify. This task allows us to assign wavelengths to spectral lines.
echelle> epar ecidentify
_______________________________________________________
PACKAGE = echelle
TASK = ecidentify
images = wave.ec Images containing features to be identified
(databas= database) Database in which to record feature data
(coordli= thar_uves.csv) User coordinate list
(units = ) Coordinate units
(match = 1.) Coordinate list matching limit in user units
(maxfeat= 100) Maximum number of features for automatic identification
(zwidth = 10.) Zoom graph width in user units
(ftype = emission) Feature type
(fwidth = 4.) Feature width in pixels
(cradius= 5.) Centering radius in pixels
(thresho= 10.) Feature threshold for centering
(minsep = 2.) Minimum pixel separation
(functio= chebyshev) Coordinate function
(xorder = 2) Order of coordinate function along dispersion
(yorder = 2) Order of coordinate function across dispersion
(niterat= 0) Rejection iterations
(lowreje= 3.) Lower rejection sigma
(highrej= 3.) Upper rejection sigma
(autowri= no) Automatically write to database?
(graphic= stdgraph) Graphics output device
(cursor = ) Graphics cursor input
(mode = ql)
_______________________________________________________
On the first run of this task, we will edit only first and third line. As I mentioned, input file is LAMP,WAVE, but now we are working with extracted spectra - wave.ec.fits. Coordlist is the list of all spectral lines in ThAr spectra. Some observatories/instruments (such as UVES) have their own line lists (UVES' can be found here). Downloaded file (*.dat or *.csv) move to the working directory with data and write its name as value for coordli parameter. If your instrument does not have such a list, the default option of this parameter will be sufficient (This applies only if your calibration spectra are spectra of ThAr lamp).
The biggest problem is always to match the very first line. This can be simplified by comparing ccdimages (spectra before extraction - wave.fits) with image of atlas' spectra. You will know then how your spectra were oriented in the ccd and which aperture corresponds to which row in atlas.
Lets run the task:
echelle> ecidentify wave.ec coord=thar_uves.csv
New window will be opened and you will see extracted spectrum from the first aperture.
You can move between the apertures by pressing 'k' or 'j'. Now you have to identify some lines. From comparing my spectra with atlas I know, that my second aperture corresponds with first row in atlas. By looking at each one I know, that my spectra are oriented 'upside-down'.
Now I know, which line of my spectra is which one in the atlas. I hover cursor above one of lines in my spectra and press 'm' (as for mark). In lower left corner you will see the line identification and you can write down its wavelength (which is shown in atlas above each line) and press enter.
Do it for three or four lines in the aperture. If you make a mistake and write bad wavelength press 'd' with cursor above bad line. If you want to move between marked lines do it by pressing '+' or '-'. When you are done with aperture, move to another one. Do it for at least half of apertures, though you have to mark lines at the beginning, middle and end apertures.
***NOTE***
Spectral lines in the observed spectra do not have to be as strong as it is in atlas spectra. For this reason you can find much more intense line in your spectra that it is in atlas or vice-versa. Also you will probably need to zoom your spectra to see some fainter lines.
***
After the last aperture press 'f' to fit marked points. By default the fitting order is set to 2.
At this order, the fit should look something like this. If there are some points which are not on the curve you can delete them by 'd'. Then you can increase the fitting order for both axes by :xorder 3 and :yorder 3 and press 'f' again.
Note that residuals are very small. Again, if there are some bad points you can delete them. Higher orders are not needed for now.
***NOTE***
You have to match some lines at the ends of apertures!
***
Quit fitting window and then press 'l' to match the coordinate list.
***NOTE***
Every time you are matching coordinate list do it with fitting order higher than 2!
***
It will match 100 lines. Now you can again fit new points and delete bad ones. At this point, you can again start to match lines manually but this time, the task will predict the wavelength of line and you can change it to more accurate value or you can confirm it.
Go through all apertures and then fit new points. When you are done, quit this task and answer yes for the question. We will run this task again but firstly we change the parameter maxfeat from 100 to 500.
Run task and match coordinate list by 'l'. Fit new points and now you can increase order also to 5. Higher values are not of much use. After deleting bad points, order 2 fit should look like this
Repeat this process for maxfeat 1000 or higher. When you are satisfied with your fitted points you are done with this task.
In case, you have more datasets with the same calibration images you do not have to do this very time-taking procedure for each one. For this purpose there is task ecreidentify.
echelle> epar ecreidentify
_______________________________________________________
PACKAGE = echelle
TASK = ecreidentify
images = new_wave.ec Spectra to be reidentified
referenc= wave.ec Reference spectrum
(shift = 0.) Shift to add to reference features
(cradius= 5.) Centering radius
(thresho= 10.) Feature threshold for centering
(refit = yes) Refit coordinate function?
(databas= /path/to/the/directory/with_reference_data/database) Database
(logfile= STDOUT,logfile) List of log files
(mode = ql)
_______________________________________________________
Input image in this case will be new calibration spectrum, reference spectrum will be already calibrated calibration spectrum (but without .fits extension!) and database will be set by the path to the directory containing the reference spectrum. You know the work was done correctly if the RMS is low. You can also check the new calibration spectrum in ecidentify task if spectral features have assigned correct wavelengths.
The task that will assign this lamp spectra to our scientific observation is refspectra.
echelle> epar refspectra
_______________________________________________________
PACKAGE = echelle
TASK = refspectra
input = object.ec List of input spectra
(referen= wave.ec ) List of reference spectra
(apertur= ) Input aperture selection list
(refaps = ) Reference aperture selection list
(ignorea= yes) Ignore input and reference apertures?
(select = interp) Selection method for reference spectra
(sort = jd) Sort key
(group = ljd) Group key
(time = no) Is sort key a time?
(timewra= 17.) Time wrap point for time sorting
(overrid= no) Override previous assignments?
(confirm= yes) Confirm reference spectrum assignments?
(assign = yes) Assign the reference spectra to the input spectrum?
(logfile= STDOUT,logfile) List of logfiles
(verbose= no) Verbose log output?
answer = Accept assignment?
(mode = ql)
_______________________________________________________
From these parameters we will change only input and referen. Input will be our extracted scientific spectra and referen will be calibrated lamp spectra.
This task only adds keywords to the header so no new exciting window will be opened.
echelle> refspectra object.ec referen=wave.ec
The last step in wavelength calibration is task dispcor.
echelle> epar dispcor
_______________________________________________________
PACKAGE = echelle
TASK = dispcor
input = List of input spectra
output = List of output spectra
(lineari= yes) Linearize (interpolate) spectra?
(databas= 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 total flux?
(blank = 0.) Output value of points not in input
(samedis= no) Same dispersion in all apertures?
(global = no) Apply global defaults?
(ignorea= no) Ignore apertures?
(confirm= no) Confirm dispersion coordinates?
(listonl= no) List the dispersion coordinates only?
(verbose= yes) Print linear dispersion assignments?
(logfile= ) Log file
(mode = ql)
_______________________________________________________
Again, we will change only first two parameters. Input will be our scientific observation and output will be name of our calibrated spectrum.
echelle> dispcor object.ec out=object.ecw
Now, you can check your calibrated data with splot task.
*****
Continuum normalization - not the most important step (yes, it is not of much use)
We have wavelength calibrated spectra but we can see, that spectra in all apertures have their own peaks at the sensitivity peak of ccd. We can normalize our spectra with the task continuum.
(maybe it is not the best continuum-removed spectrum but I do not have better at the moment - I will change picture as soon as possible)
echelle> epar continuum
_______________________________________________________
PACKAGE = echelle
TASK = continuum
input = Input images
output = Output images
(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?
(wavesca= yes) Scale the X axis with wavelength?
(logscal= no) Take the log (base 10) of both axes?
(overrid= no) Override previously fit lines?
(listonl= no) List fit but don't modify any images?
(logfile= logfile) List of log files
(interac= yes) Set fitting parameters interactively?
(sample = *) Sample points to use in fit
(naverag= 1) Number of points in sample averaging
(functio= spline3) Fitting function
(order = 1) Order of fitting function
(low_rej= 2.) Low rejection in sigma of fit
(high_re= 0.) High rejection in sigma of fit
(niterat= 10) Number of rejection iterations
(grow = 1.) Rejection growing radius
(markrej= yes) Mark rejected points?
(graphic= stdgraph) Graphics output device
(cursor = ) Graphics cursor input
ask = yes
(mode = ql)
_______________________________________________________
Set the input to the object.ecw.fits and output to the object.ecwc.fits. Fitting options can be changed interactively later.
After executing task, fitting window will be opened:
echelle> continuum object.ecw out=object.ecwc
Dashed line is the resulting fit. Fitting order can be changed with ':o n' and number of iterations can be changed with ':niterate n'. Emission lines can be automatically rejected with increasing ':high n' and absorption lines with increasing ':low n'. Fit each aperture respectively and it will probably take more attempts (yes, you guessed it correctly, for appropriate normalization you will run this task for hundreds times).
Flux calibration
So now, we are getting to the last step of our data reduction. Also, we are getting beck to the begging but one thing at the time.
Our data are nicely extracted and wavelength calibrated. However, flux units are really useless. Flux calibration calibrate the flux and converts it into erg.s^-1.cm^-2.
For this purpose, you need tho have sensitivity functions for your specific instrument and its specific setting. Some observatories/instruments offer these functions as the auxiliary data with scientific data. This is also the case of ESO/UVES data. There is one big problem! These data cannot be used in IRAF. I did not mentioned it before but ESO have its own software specialized to data reduction and analysis. It is the REFLEX software and as well as IRAF it runs under the Linux. I tried to use this software instead of IRAF.
Whoever said that IRAF is user-unfriendly had never used REFLEX before!!!
***
NOTE-AFTER-A-WHILE: Maybe it is not that user-unfriendly but I still hate it.
***
For this reason, the sensitivity functions obtained with scientific data are useless. These can be created from the observation of standard stars. These observations can be downloaded in the same way as any other data. Problem of this process is that we have to reduce absolutely new dataset from the beginning.
At first, you have to find data of standard star. List of standard stars can be found for example here. List of standard stars in the IRAF can be found after typing 'page onedstds$README' in IRAF terminal.
Choose one of them and download data obtained with the SAME instrument and the SAME configuration of this instrument as the data you want to calibrate.
Run all the task and processes you have run on your scientific data and calibrate the wavelength. These calibrated standard star data will be used for creating sensitivity functions.
***NOTE***
Work in new working directory.
***
If we name downloaded data similarly to our scientific data (stand_Bias1-5.fits stand_Flat1-5.fits stand_orderdef.fits stand_wave.fits stand_object.fits) we should end up with file named stand_object.ecw.fits (not continuum normalized!!!)
Now, we will run two tasks - standard and sensfunc.
echelle> epar standard
_______________________________________________________
PACKAGE = echelle
TASK = standard
input = stand_object.ecw Input image file root name
output = std Output flux file (used by SENSFUNC)
(samesta= yes) Same star in all apertures?
(beam_sw= no) Beam switch spectra?
(apertur= ) Aperture selection list
(bandwid= INDEF) Bandpass widths
(bandsep= INDEF) Bandpass separation
(fnuzero= 3.6800000000000E-20) Absolute flux zero point
(extinct= )_.extinction) Extinction file
(caldir = onedstds$caldir/) Directory containing calibration data
(observa= )_.observatory) Observatory for data
(interac= yes) Graphic interaction to define new bandpasses
(graphic= stdgraph) Graphics output device
(cursor = ) Graphics cursor input
star_nam= stands Star name in calibration list
airmass = Airmass
exptime = Exposure time (seconds)
mag = Magnitude of star
magband = Magnitude type
teff = Effective temperature or spectral type
answer = no (no|yes|NO|YES|NO!|YES!)
(mode = ql)
_______________________________________________________
The input file will be our stand_object.ecw.fits. Other parameters we need to change are (caldir and star_name. For the (caldir we need to list standard stars in IRAF 'page onedstds$README'. Find your star and remember the directory it is saved in. Lets say, we used data of star stands and its IRAF files are saved in 'caldir' directory. The values of these parameter will be '(caldir=onedstds$caldir/' (DO NOT forget slash / and the end) and 'star_name=stands'. You can also change parameters bandwid and bandsep both to 5 or 10.
Run the task and if you want to check the automatically set bandpasses, answer yes for the question. New window will be opened and you can change your bandpasses.
echelle> standard stand_object.ecw bandwid=5 badsep=5 (caldir=onedstds$caldir/ star_name=stands answer=yes
If you are happy with them, you can proceed to the next task.
echelle> epar sensfunc
_______________________________________________________
PACKAGE = echelle
TASK = sensfunc
standard= std Input standard star data file (from STANDARD)
sensitiv= sens Output root sensitivity function imagename
(apertur= ) Aperture selection list
(ignorea= no) Ignore apertures and make one sensitivity function?
(logfile= logfile) Output log for statistics information
(extinct= )_.extinction) Extinction file
(newexti= extinct.dat) Output revised extinction file
(observa= )_.observatory) Observatory of data
(functio= spline3) Fitting function
(order = 6) Order of fit
(interac= yes) Determine sensitivity function interactively?
(graphs = sr) Graphs per frame
(marks = plus cross box) Data mark types (marks deleted added)
(colors = 2 1 3 4) Colors (lines marks deleted added)
(cursor = ) Graphics cursor input
(device = stdgraph) Graphics output device
answer = yes (no|yes|NO|YES)
(mode = ql)
_______________________________________________________
In this task we do not need to change any parameter.
echelle> sensfunc standard=std sensitiv=sens
Number of fitted points depends on number of bandpasses defined previously. Smaller bandpasses you set, more points you get.
This task creates one sensitivity function for every aperture and these functions are saved in working directory with names like sens.0001 - sens.00n. Now we have functions we need and we can calibrate our scientific data. Created sensitivity functions move to the new folder and save it for future use. All the sensitivity functions also copy to the directory with your scientific data and go there.
Last task in this calibration will be calibrate.
echelle> epar calibrate
_______________________________________________________
PACKAGE = echelle
TASK = calibrate
input = object.ecw Input spectra to calibrate
output = spectrum Output calibrated spectra
(extinct= yes) Apply extinction correction?
(flux = yes) Apply flux calibration?
(extinct= )_.extinction) Extinction file
(observa= )_.observatory) Observatory of observation
(ignorea= no) Ignore aperture numbers in flux calibration?
(sensiti= sens) Image root name for sensitivity spectra
(fnu = no) Create spectra having units of FNU?
airmass = Airmass
exptime = Exposure time (seconds)
(mode = ql)
_______________________________________________________
Here we define only input and output file names. (sensiti parameter is set to the root name of our sensitivity functions - sens.
Run the task:
echelle> calibrate object.ecw out=spectrum
Input spectra to calibrate (object.ecw):
Output calibrated spectra (spectrum):
spectrum: EG21
Extinction correction applied
Flux calibration applied
You can check your calibrated spectra with the splot task. Do not panic when the edges of spectra in all apertures will be kind of messy. The sensitivity of the chip is much lower then in the center and the edges can be cut off without loose of any data. Though, it will be less messy if the bandpasses in standard task will be small and evenly distributed from one edge to another.
Well, congratulations to all of you, who get through all these posts and even bigger congratulations if you have done all the reduction steps described previously. In that case, you have your data calibrated and you can finally start working with them. As you can see, you do not have one beautiful spectrum you was hoping for but you have it cut into smaller parts. In the next post, I will try to combine all these parts into one spectrum and I will show you some features of splot task.
As always, feel free to leave a comment or ask any question. In case you found any mistake in my reduction process, please, let me know (especially the flux calibration is very doubtful).
Have a nice day
Jakub
No comments:
Post a Comment