Cookbook for 0.9m MOSAIC Reduction of
WOCS Data
The following is a "cookbook"
document on the processing steps that
I used for the WOCS project MOSAIC reductions.
-- based on Buell Januzzi's web page:
http://www.noao.edu/noao/noaodeep/ReductionOpt/frames.html
-- updates being made 7/2002 for the Feb 2000 WOCS run
-- second pass updates 8/2002
NOTE: All time estimates were
made on a dual processor Sparc
Ultra60 with 1 GB of memory, 1.6 GB of swap
space and that was
not running anything besides open windows and
IRAF.
-- third update made after upgrading to IRAF 2.12.1
3/2003
and mscred V4.8 March 14, 2003
NOTE: All of the steps were
re-verified and param file examples
updated after these upgrades.
One thing not checked is when
wild cards (*) rather than @ lists
(@llobj) are valid, i.e.,
any input listed as a wild card does
work however cases where
@ lists are used instead may actually
allow wild cards now.
Experiment if you like, but inputs
are known to work as shown.
-- Version 4.1 adds fix in msccmatch (4/2003)
-- Version 5.0 is the web page version, the one you are currently
reading. (7/2003)
0) MOSAIC reductions require a lot of free disk space (about 40 GB
per
night), a fast cpu (the faster the better), and at least
1 GB of
memory to avoid excess paging, which will significantly slow
down
the process.
- load a night, or even better a run onto your computer.
- start IRAF. You will need a version that has mscred installed,
and preferably some other packages (see below).
1) After loading all of one night's data into a directory, first
check
the integrity of the data with ccdlist
PACKAGE = mscred
TASK = ccdlist
images =
*.fits CCD images to listed
(ccdtype=
) CCD image type to be listed
(extname=
im1) Extension name pattern
(names =
no) List image names only?
(long =
no) Long format listing?
(mode =
ql)
- this routine does not give any output until it has done ALL of the
files, so be patient
- the output is wider than 80 columns (around 110) so expanding the window
to this length will make examining the output easier
- examine the output to be sure that all of the images actually are complete.
The size of all the files should be the same (any which are smaller
are not
complete), and the short header printed could be complete.
Delete any images which are not.
2) Create a subdirectory called Raw and copy all of the *.fits files
to it.
The standard script automatically makes a copy of the
data to Raw/
but it does this as it goes, so if something crashes
during the
first step, it can result in a mess.
mkdir Raw
!cp *.fits Raw &
3) Initialize the Software -- in mscred
- Enter mscred
mscred
- This only has to be done once, not for each night
epar setinstrument
PACKAGE = mscred
TASK = setinstrument
site =
kpno Site (? for menu)
telescop=
36inch Telescope (? for menu)
instrume= CCDMosaThin1
Instrument (? for a list)
(directo= mscdb$noao/)
Instrument directory
(review =
yes) Review instrument parameters?
query_si=
kpno Site (? for menu or q to quit)
query_te=
Telescope (? for menu or q to quit)
query_in= CCDMosaThin1
Instrument (? for menu or q to quit)
(mode =
ql)
:g
PACKAGE = user
TASK = mscred
(pixelty=
real real) Output and calculation pixel datatypes
(verbose=
yes) Print log information to the standard output?
(logfile=
logfile) Text log file
(plotfil=
) Log metacode plot file
(backup =
none) Backup data (none|once|all)?
(bkuproo=
Raw/) Backup root (directory or prefix)
(instrum= mscdb$noao/kpno/36inch/CCDMosaThin1.dat) CCD instrument file
(ampfile=
amps) Amplifier translation file
(ssfile =
subsets) Subset translation file
(im_bufs=
4.) Image I/O buffer size (in Mbytes)
(graphic=
stdgraph) Interactive graphics output device
(cursor =
) Graphics cursor input
(version= V4.8: March 14, 2003)
(mode =
ql)
:g
PACKAGE = mscred
TASK = ccdproc
images =
*.fits List of Mosaic CCD images to process
(output =
) List of output processed images
(bpmasks=
) List of output bad pixel masks
(ccdtype=
) CCD image type to process
(noproc =
no) List processing steps only?
(xtalkco=
yes) Apply crosstalk correction?
(fixpix =
yes) Apply bad pixel mask correction?
(oversca=
yes) Apply overscan strip correction?
(trim =
yes) Trim the image?
(zerocor=
no) Apply zero level correction?
(darkcor=
no) Apply dark count correction?
(flatcor=
no) Apply flat field correction?
(sflatco=
no) Apply sky flat field correction?
(split =
no) Use split images during processing?
(merge =
no) Merge amplifiers from same CCD?
(xtalkfi=
!xtalkfil) Crosstalk file
(fixfile=
BPM) List of bad pixel masks
(saturat=
INDEF) Saturated pixel threshold
(sgrow =
0) Saturated pixel grow radius
(bleed =
INDEF) Bleed pixel threshold
(btrail =
15) Bleed trail minimum length
(bgrow =
1) Bleed pixel grow radius
(biassec=
!biassec) Overscan strip image section
(trimsec=
!trimsec) Trim data section
(zero =
Zero) List of zero level calibration images
(dark =
Dark) List of dark count calibration images
(flat =
Flat*) List of flat field images
(sflat =
Sflat*) List of secondary flat field images
(minrepl=
1.) Minimum flat field value
(interac=
no) Fit overscan interactively?
(functio=
legendre) Fitting function
(order =
1) Number of polynomial terms or spline pieces
(sample =
*) Sample points to fit
(naverag=
1) Number of sample points to combine
(niterat=
1) Number of rejection iterations
(low_rej=
3.) Low sigma rejection factor
(high_re=
3.) High sigma rejection factor
(grow =
0.) Rejection growing radius
(fd =
)
(fd2 =
)
(mode =
ql)
:wq
This is what to INITIALLY set ccdproc to. It will change as the steps
progress.
DON'T ACTUALLY RUN ANYTHING YET!
4) Determine the saturation level for each chip.
Saturated pixels have to be added to the individual
bad pixel masks
to avoid ringing during image projection (and to be
exclude from
combination if one makes a stack). Even though
this information isn't
used until much later in the process, it has to be determined
BEFORE
any processing is done to the chips.
For some reason, the saturation level is not fixed.
Generally it
should be the same for any given run (so you only need
determine it
once) but one should probably re-determine it for subsequent
runs.
Determine the saturation level using either implot or
mscexamine.
When using mscexam be sure to epar rimexam and turn
off the
"fit and subtract background" option, i.e., "backgro=no",
so that the
values displayed in the graphs are the REAL values of
the data, and
not the adjusted values.
[EXCEPT this doesn't work. mscexam seems to ignore
the parameters
set
with rimexam, and it also doesn't have a "backgro" colon
command.]
So, use imexam, and display each chip separately.
NOTE: The chips are numbered from the lower left:
5 6 7 8
1 2 3 4
The saturation limit on different stars my be different
even on the
same chip! Try to examine at least three, but
if there are more,
examine enough to get a good feel for both the mean
and scatter.
The saturation value in the header should be set to
a value LOWER
then the lowest value in the scatter. Be conservative,
since the
linearity right below saturation will likely be dubious
in any case.
Values from February 2000
Values from June 1999 June 2002
Chip 1: 34,000
Chip 1: 34,000 Chip 1: 38,000
Chip 2: 26,000
Chip 2: 26,000 Chip 2: 31,000
Chip 3: 38,000
Chip 3: 36,000 Chip 3: 41,000
Chip 4: 27,000
Chip 4: 29,000 Chip 4: 32,000
Chip 5: 33,000
Chip 5: 33,000 Chip 5: 36,000
Chip 6: 26,000
Chip 6: 27,000 Chip 6: 31,000
Chip 7: 26,000
Chip 7: 30,000 Chip 7: 28,000
Chip 8: 25,000
Chip 8: 26,000 Chip 8: 30,000
Once you have determined the saturation value for each
chip, use hedit
to change the value of the header keywords for each
extension
in each MEF header. Double check and be sure each
header has
a SATURATE keyword, because it is sometimes left out
of a subset of
headers.
hedit *.fits[1] saturate 34000 verify- add- del- update+
hedit *.fits[2] saturate 26000 verify- add- del- update+
hedit *.fits[3] saturate 38000 verify- add- del- update+
hedit *.fits[4] saturate 27000 verify- add- del- update+
hedit *.fits[5] saturate 33000 verify- add- del- update+
hedit *.fits[6] saturate 26000 verify- add- del- update+
hedit *.fits[7] saturate 26000 verify- add- del- update+
hedit *.fits[8] saturate 25000 verify- add- del- update+
You can type all of these lines separately, or put them
in a script,
and execute the script:
cl < sat.cl
5) Now copy the correct path to the cross talk file into the XTALKFI
keyword in the headers of all of the raw data.
Use hedit to do this,
and be sure to change the keyword in all of the extensions,
[1] - [8].
[If the crosstalk file is not down this path, set up
the path, and
get the correct cross talk file from the NOAO web page,
currently at:
http://www.noao.edu/noao/mosaic/calibs.html
] As of this writing,
there are different cross talk files for March 1999, Sept.
2000,
Oct. 2000, and Feb. 2001
hedit *.fits[1] XTALKFIL mscdb$noao/CCDMosaThin1/xtalk9903 ver- add+ del-
upd+
hedit *.fits[2] XTALKFIL mscdb$noao/CCDMosaThin1/xtalk9903 ver- add+ del-
upd+
hedit *.fits[3] XTALKFIL mscdb$noao/CCDMosaThin1/xtalk9903 ver- add+ del-
upd+
hedit *.fits[4] XTALKFIL mscdb$noao/CCDMosaThin1/xtalk9903 ver- add+ del-
upd+
hedit *.fits[5] XTALKFIL mscdb$noao/CCDMosaThin1/xtalk9903 ver- add+ del-
upd+
hedit *.fits[6] XTALKFIL mscdb$noao/CCDMosaThin1/xtalk9903 ver- add+ del-
upd+
hedit *.fits[7] XTALKFIL mscdb$noao/CCDMosaThin1/xtalk9903 ver- add+ del-
upd+
hedit *.fits[8] XTALKFIL mscdb$noao/CCDMosaThin1/xtalk9903 ver- add+ del-
upd+
Again, you can do this on the command line, or put the
commands in a
script and execute the script:
cl < xtalk.cl
6) Preliminary Processing for Flat and Zero (Bias) frames
Since we don't have to create individual bad pixel masks
for the
calibration images, we will process them first, and
create a bias
frame, a dome flat for each band, and a sky flat for
each band.
(Darks are not necessary for MOSAIC data since the dark
current
is vanishingly small.)
The first step is to use ccdproc to apply crosstalk
and overscan
correction and trim the images.
NOTE: This step will double the size of the zero
and flat files
because it changes the precision of the floating point
data.
Be sure you have sufficient disk space for this step
(likely on
the order of 5 GB).
epar ccdproc
PACKAGE = mscred
TASK = ccdproc
images = dflat*.fits,zero*.fits,sflat*.fits List of Mosaic
CCD images to process
(output =
) List of output processed images
(bpmasks=
) List of output bad pixel masks
(ccdtype=
) CCD image type to process
(noproc =
no) List processing steps only?
(xtalkco=
yes) Apply crosstalk correction?
(fixpix =
yes) Apply bad pixel mask correction?
(oversca=
yes) Apply overscan strip correction?
(trim =
yes) Trim the image?
(zerocor=
no) Apply zero level correction?
(darkcor=
no) Apply dark count correction?
(flatcor=
no) Apply flat field correction?
(sflatco=
no) Apply sky flat field correction?
(split =
no) Use split images during processing?
(merge =
no) Merge amplifiers from same CCD?
(xtalkfi=
!xtalkfil) Crosstalk file
(fixfile=
BPM) List of bad pixel masks
(saturat=
INDEF) Saturated pixel threshold
(sgrow =
0) Saturated pixel grow radius
(bleed =
INDEF) Bleed pixel threshold
(btrail =
15) Bleed trail minimum length
(bgrow =
1) Bleed pixel grow radius
(biassec=
!biassec) Overscan strip image section
(trimsec=
!trimsec) Trim data section
(zero =
Zero) List of zero level calibration images
(dark =
Dark) List of dark count calibration images
(flat =
Flat*) List of flat field images
(sflat =
Sflat*) List of secondary flat field images
(minrepl=
1.) Minimum flat field value
(interac=
no) Fit overscan interactively?
(functio=
legendre) Fitting function
(order =
1) Number of polynomial terms or spline pieces
(sample =
*) Sample points to fit
(naverag=
1) Number of sample points to combine
(niterat=
1) Number of rejection iterations
(low_rej=
3.) Low sigma rejection factor
(high_re=
3.) High sigma rejection factor
(grow =
0.) Rejection growing radius
(fd =
)
(fd2 =
)
(mode =
ql)
:g
7) Generating a Zero (Bias) Frame with zerocombine and Subtracting
- Inspect (imstat, iterstat) the zero* frames to be sure they are
all good.
Remove any that are not.
iterstat zero*.fits[6]
<--- any extension can be used
- Use zerocombine to generate a zero frame
epar zerocombine
PACKAGE = mscred
TASK = zerocombine
input =
zero* List of zero level images to combine
(output = Zero20000225)
Output zero level name
(combine=
average) Type of combine operation
(reject =
crreject) Type of rejection
(ccdtype=
zero) CCD image type to combine
(process=
no) Process images before combining?
(delete =
no) Delete input images after combining?
(scale =
none) Image scaling
(statsec=
) Image section for computing statistics
(nlow =
1) minmax: Number of low pixels to reject
(nhigh =
1) minmax: Number of high pixels to reject
(nkeep =
1) Minimum to keep (pos) or maximum to reject (neg)
(mclip =
yes) Use median in sigma clipping algorithms?
(lsigma =
3.) Lower sigma clipping factor
(hsigma =
3.) Upper sigma clipping factor
(rdnoise=
RDNOISE) ccdclip: CCD readout noise (electrons)
(gain =
GAIN) ccdclip: CCD gain (electrons/DN)
(snoise =
0.) ccdclip: Sensitivity noise (fraction)
(pclip =
-0.5) pclip: Percentile clipping parameter
(blank =
0.) Value if there are no pixels
(mode =
ql)
:g
- Use imstat or iterstat to check that the flats aren't saturated
and have sufficient exposure level to be significant.
iterstat sflat*.fits[6] <--- any extension can be used
iterstat dflat*fits[6]
delete any saturated or underexposed flats.
- Now run ccdproc to do the zero subtraction from all of the
flat frames
only. Make a list of the flat frames (dome
and sky) (llzsubflat) and
ls *flat*fits > llzsubflat
- then
epar ccdproc
PACKAGE = mscred
TASK = ccdproc
images = @llzsubflat
List of Mosaic CCD images to process
(output =
) List of output processed images
(bpmasks=
) List of output bad pixel masks
(ccdtype=
) CCD image type to process
(noproc =
no) List processing steps only?
(xtalkco=
no) Apply crosstalk correction?
(fixpix =
no) Apply bad pixel mask correction?
(oversca=
no) Apply overscan strip correction?
(trim =
no) Trim the image?
(zerocor=
yes) Apply zero level correction?
(darkcor=
no) Apply dark count correction?
(flatcor=
no) Apply flat field correction?
(sflatco=
no) Apply sky flat field correction?
(split =
no) Use split images during processing?
(merge =
no) Merge amplifiers from same CCD?
(xtalkfi=
!xtalkfil) Crosstalk file
(fixfile=
BPM) List of bad pixel masks
(saturat=
INDEF) Saturated pixel threshold
(sgrow =
0) Saturated pixel grow radius
(bleed =
INDEF) Bleed pixel threshold
(btrail =
15) Bleed trail minimum length
(bgrow =
1) Bleed pixel grow radius
(biassec=
!biassec) Overscan strip image section
(trimsec=
!trimsec) Trim data section
(zero = Zero20000225)
List of zero level calibration images
(dark =
Dark) List of dark count calibration images
(flat =
Flat*) List of flat field images
(sflat =
Sflat*) List of secondary flat field images
(minrepl=
1.) Minimum flat field value
(interac=
no) Fit overscan interactively?
(functio=
legendre) Fitting function
(order =
1) Number of polynomial terms or spline pieces
(sample =
*) Sample points to fit
(naverag=
1) Number of sample points to combine
(niterat=
1) Number of rejection iterations
(low_rej=
3.) Low sigma rejection factor
(high_re=
3.) High sigma rejection factor
(grow =
0.) Rejection growing radius
(fd =
)
(fd2 =
)
(mode =
ql)
:g
8) Darks were not taken so no dark subtraction is performed.
The general
wisdom is that the dark current is so low with MOSAIC
that dark
subtraction is unnecessary.
9) Creating Dome Flats with flatcombine and Flattening the Data
- if you didn't do so above, inspect the dflat* frames to
be sure
they are all good (below the saturation limit.
Remove any that
are not.
- create a dome flat for each band using flatcombine and then
run
ccdproc to do the flat field processing on the
data.
epar flatcombine
PACKAGE = mscred
TASK = flatcombine
input =
dflat*.fits List of flat field images to combine
(output = Dflat20000225) Output
flat field root name
(combine=
average) Type of combine operation
(reject =
ccdclip) Type of rejection
(ccdtype=
flat) CCD image type to combine
(process=
no) Process images before combining?
(subsets=
yes) Combine images by subset parameter?
(delete =
no) Delete input images after combining?
(scale =
mode) Image scaling
(statsec=
) Image section for computing statistics
(nlow =
1) minmax: Number of low pixels to reject
(nhigh =
1) minmax: Number of high pixels to reject
(nkeep =
1) Minimum to keep (pos) or maximum to reject (neg)
(mclip =
no) Use median in sigma clipping algorithms?
(lsigma =
3.) Lower sigma clipping factor
(hsigma =
3.) Upper sigma clipping factor
(rdnoise=
rdnoise) ccdclip: CCD readout noise (electrons)
(gain =
gain) ccdclip: CCD gain (electrons/DN)
(snoise =
0.) ccdclip: Sensitivity noise (fraction)
(pclip =
-0.5) pclip: Percentile clipping parameter
(blank =
1.) Value if there are no pixels
(mode =
q)
:g
- now do the preliminary processing of the object frames,
flatten them,
and make individual bad pixel masks.
- create a list of of subdirectory names for the individual
bad
pixel masks (llmasks).
ls obj*fits > llmasks
edit the file replacing "obj" by "bpm" and removing
the ".fits"
then;
epar ccdproc
PACKAGE = mscred
TASK = ccdproc
images =
obj*fits List of Mosaic CCD images to process
(output =
) List of output processed images
(bpmasks=
@llmasks) List of output bad pixel masks
(ccdtype=
object) CCD image type to process
(noproc =
no) List processing steps only?
(xtalkco=
yes) Apply crosstalk correction?
(fixpix =
yes) Apply bad pixel mask correction?
(oversca=
yes) Apply overscan strip correction?
(trim =
yes) Trim the image?
(zerocor=
yes) Apply zero level correction?
(darkcor=
no) Apply dark count correction?
(flatcor=
yes) Apply flat field correction?
(sflatco=
no) Apply sky flat field correction?
(split =
no) Use split images during processing?
(merge =
no) Merge amplifiers from same CCD?
(xtalkfi=
!xtalkfil) Crosstalk file
(fixfile=
BPM) List of bad pixel masks
(saturat=
!saturate) Saturated pixel threshold
(sgrow =
0) Saturated pixel grow radius
(bleed =
20000) Bleed pixel threshold
(btrail =
15) Bleed trail minimum length
(bgrow =
1) Bleed pixel grow radius
(biassec=
!biassec) Overscan strip image section
(trimsec=
!trimsec) Trim data section
(zero = Zero20000225)
List of zero level calibration images
(dark =
Dark) List of dark count calibration images
(flat = Dflat20000225*)
List of flat field images
(sflat =
Sflat*) List of secondary flat field images
(minrepl=
1.) Minimum flat field value
(interac=
no) Fit overscan interactively?
(functio=
legendre) Fitting function
(order =
1) Number of polynomial terms or spline pieces
(sample =
*) Sample points to fit
(naverag=
1) Number of sample points to combine
(niterat=
1) Number of rejection iterations
(low_rej=
3.) Low sigma rejection factor
(high_re=
3.) High sigma rejection factor
(grow =
0.) Rejection growing radius
(fd =
)
(fd2 =
)
(mode =
ql)
:g
- Now, since a sky flat correction is also needed, but we
don't
want the dome flat information to still be part
of the sky flats, the
sflat*.fit files also need to be dome flat corrected.
Unfortunately,
ccdproc refuses to dome flat correct the sky flats
unless the keyword
OBSTYPE is set to 'object' (no matter what is
specified in ccdproc), so
edit the headers:
hedit sflat*.fits[1] OBSTYPE object verify- add- del- update+
hedit sflat*.fits[2] OBSTYPE object verify- add- del- update+
hedit sflat*.fits[3] OBSTYPE object verify- add- del- update+
hedit sflat*.fits[4] OBSTYPE object verify- add- del- update+
hedit sflat*.fits[5] OBSTYPE object verify- add- del- update+
hedit sflat*.fits[6] OBSTYPE object verify- add- del- update+
hedit sflat*.fits[7] OBSTYPE object verify- add- del- update+
hedit sflat*.fits[8] OBSTYPE object verify- add- del- update+
You can type all of these lines separately, or put them
in a script,
and execute the script:
cl < skyflatfix.cl
- then run ccdproc to do the dome flat corrections to the
sky flats
epar ccdproc
PACKAGE = mscred
TASK = ccdproc
images =
sflat*fits List of Mosaic CCD images to process
(output =
) List of output processed images
(bpmasks=
) List of output bad pixel masks
(ccdtype=
) CCD image type to process
(noproc =
no) List processing steps only?
(xtalkco=
no) Apply crosstalk correction?
(fixpix =
no) Apply bad pixel mask correction?
(oversca=
no) Apply overscan strip correction?
(trim =
no) Trim the image?
(zerocor=
no) Apply zero level correction?
(darkcor=
no) Apply dark count correction?
(flatcor=
yes) Apply flat field correction?
(sflatco=
no) Apply sky flat field correction?
(split =
no) Use split images during processing?
(merge =
no) Merge amplifiers from same CCD?
(xtalkfi=
!xtalkfil) Crosstalk file
(fixfile=
BPM) List of bad pixel masks
(saturat=
INDEF) Saturated pixel threshold
(sgrow =
0) Saturated pixel grow radius
(bleed =
INDEF) Bleed pixel threshold
(btrail =
15) Bleed trail minimum length
(bgrow =
1) Bleed pixel grow radius
(biassec=
!biassec) Overscan strip image section
(trimsec=
!trimsec) Trim data section
(zero = Zero20000225)
List of zero level calibration images
(dark =
Dark) List of dark count calibration images
(flat = Dflat20000225*)
List of flat field images
(sflat =
Sflat*) List of secondary flat field images
(minrepl=
1.) Minimum flat field value
(interac=
no) Fit overscan interactively?
(functio=
legendre) Fitting function
(order =
1) Number of polynomial terms or spline pieces
(sample =
*) Sample points to fit
(naverag=
1) Number of sample points to combine
(niterat=
1) Number of rejection iterations
(low_rej=
3.) Low sigma rejection factor
(high_re=
3.) High sigma rejection factor
(grow =
0.) Rejection growing radius
(fd =
)
(fd2 =
)
(mode =
ql)
:g
10) Examine some subset of the obj* images with mscdisplay to be
sure that
the processing has progressed correctly to this step.
You will still
notice vignetting in the corners since sky flat correction
hasn't been
done yet.
If you haven't done so in your login.cl or loginuser.cl
file, type
the following on the command line before using mscdisplay:
set stdimage=imt8800
This will set the proper size image buffer for displaying
a MOSAIC image.
NOTES: Chips 3 (and sometimes 1 and 6) are
"noisy" with some sort
of bias-level like interference. This effect seems to be
on about the +- 5 count level with mean sky of 4 in a
short V exposure, and +- 20 with a mean sky of 150 in a long
I exposure. So, this effect will add significant noise in
the sky but doesn't seem to be changing the mean sky level
or significantly affecting the flux in well detected objects.
11) Sky flattening
- If one night does not have enough sky flats to comprise a full set, or
if there aren't enough of a given band to judge their quality, stop
here and repeat steps 1-10 above as necessary on other nights until
you have enough sky frames. Process all of the nights up to
this step,
then put all of the sflat*fits images into one directory and process
them as follows. NOTE: This could require well over 100 GB
of disk
space!
The following are written as if there were enough sflat*fits images
in one night to make a complete set.
- first make a combined sky flat image
epar sflatcombine
PACKAGE = mscred
TASK = sflatcombine
input =
sflat* List of images to combine
(output = Sflat20000225) Output
sky flat field root name
(combine=
average) Type of combine operation
(reject =
ccdclip) Type of rejection
(ccdtype=
) CCD image type to combine
(subsets=
yes) Combine images by subset parameter?
(masktyp=
none) Mask type
(maskval=
0.) Mask value
(scale =
mode) Image scaling
(statsec=
) Image section for computing statistics
(nkeep =
1) Minimum to keep (pos) or maximum to reject (neg)
(nlow =
1) minmax: Number of low pixels to reject
(nhigh =
1) minmax: Number of high pixels to reject
(mclip =
yes) Use median in sigma clipping algorithms?
(lsigma =
6.) Lower sigma clipping factor
(hsigma =
3.) Upper sigma clipping factor
(rdnoise=
rdnoise) ccdclip: CCD readout noise (electrons)
(gain =
gain) ccdclip: CCD gain (electrons/DN)
(snoise =
0.) ccdclip: Sensitivity noise (fraction)
(pclip =
-0.5) pclip: Percentile clipping parameter
(blank =
1.) Value if there are no pixels
(grow =
3.) Radius (pixels) for neighbor rejection
(fd =
)
(mode =
ql)
:g
- The object frames could be flattened with the Sflat* files at this
point, but since the the only correction needed now is large
scale
illumination correction (mostly vignetting) a better process
is to fit a surface to each chip and then flatten so as not
to
introduce more pixel-to-pixel noise.
- Fitting the surfaces is done with imsurfit. A separate surface
needs
to be fit to each chip, but the final image needs to be MEF
for the
next step. To this end, use msccmd for each imsurfit
command.
unlearn imsurfit
msccmd
msccmd: imsurfit $input $output xo=3 yo=3 xm=100 ym=100
Input files: Sflat20000225U.fits
Output files: Sflat_smo20000225U.fits
...
(repeat this for each band's Sflat* file)
...
msccmd:quit
NOTE: for some reason if the $input and $output are replaced
with the
actual values, an error
occurs. So type in the commands as above.
- If you did not examine each individual sflat* file (I usually don't
because there are too many and it takes so long to display
even one
with mscdisplay), be sure to examine the smoothed images to
be sure
they look like reasonable sky flats, i.e. look for the vignetted
pattern
if reducing KPNO 0.9m data. If they aren't good, you'll
have to
examine the individual sflat* images and throw out the bad
one(s).
Then redo the above steps.
- Now apply the sflat correction to the object frames. Be sure
to
enter the SURFACES (Sflat_smo*) rather than the frames.
epar ccdproc
PACKAGE = mscred
TASK = ccdproc
images =
obj*fits List of Mosaic CCD images to process
(output =
) List of output processed images
(bpmasks=
) List of output bad pixel masks
(ccdtype=
object) CCD image type to process
(noproc =
no) List processing steps only?
(xtalkco=
no) Apply crosstalk correction?
(fixpix =
no) Apply bad pixel mask correction?
(oversca=
no) Apply overscan strip correction?
(trim =
no) Trim the image?
(zerocor=
no) Apply zero level correction?
(darkcor=
no) Apply dark count correction?
(flatcor=
no) Apply flat field correction?
(sflatco=
yes) Apply sky flat field correction?
(split =
no) Use split images during processing?
(merge =
no) Merge amplifiers from same CCD?
(xtalkfi=
!xtalkfil) Crosstalk file
(fixfile=
BPM) List of bad pixel masks
(saturat=
INDEF) Saturated pixel threshold
(sgrow =
0) Saturated pixel grow radius
(bleed =
INDEF) Bleed pixel threshold
(btrail =
15) Bleed trail minimum length
(bgrow =
1) Bleed pixel grow radius
(biassec=
!biassec) Overscan strip image section
(trimsec=
!trimsec) Trim data section
(zero = Zero20000225)
List of zero level calibration images
(dark =
Dark) List of dark count calibration images
(flat = Dflat20000225*)
List of flat field images
(sflat = Sflat_smo20000225*) List of secondary flat field
images
(minrepl=
1.) Minimum flat field value
(interac=
no) Fit overscan interactively?
(functio=
legendre) Fitting function
(order =
1) Number of polynomial terms or spline pieces
(sample =
*) Sample points to fit
(naverag=
1) Number of sample points to combine
(niterat=
1) Number of rejection iterations
(low_rej=
3.) Low sigma rejection factor
(high_re=
3.) High sigma rejection factor
(grow =
0.) Rejection growing radius
(fd =
)
(fd2 =
)
(mode =
ql)
:g
- A comparison of the pre-sky flattened images with the post-sky
flattened
images shows that the sky flats are needed to correct the
vignetting
in the chip corners. No pupil correction is needed because
there is
no field lens in the 0.9m Cassegrain configuration.
Additionally, there
is no evidence of fringing in the I-band, so no fringe correction
is
needed either.
12) Finding Cosmic Rays and Creating a final Bad Pixel Mask
- First find the cosmic rays in each object frame using craverage
then
add these cosmic ray pixels to the individual bad pixel masks
created by ccdproc. These final masks will than contain
true bad
pixels, cosmic rays, bleed trails, and saturated pixels.
These
pixels must be removed before doing the tangent-plane projection
required to put all images on the same plate scale and a uniform
scale across the image.
- NOTE: the following requires the installation of the external package
crutil.
- enter the crutil package
crutil
epar craverage
PACKAGE = crutil
TASK = craverage
input =
List of input images
output =
List of output images
(crmask =
) List of output cosmic ray and object masks
(average=
) List of output block average filtered images
(sigma =
) List of output sigma images
(navg =
10 7) Block average box size
(nrej =
10 3) Number of high pixels to reject from the averag
(nbkg =
5 5) Background annulus width
(nsig =
50 50) Box size for sigma calculation
(var0 =
0.) Variance coefficient for DN^0 term
(var1 =
0.) Variance coefficient for DN^1 term
(var2 =
0.) Variance coefficient for DN^2 term
(crval =
1) Mask value for cosmic rays
(lcrsig =
10.) Low cosmic ray sigma outside object
(hcrsig =
3.75) High cosmic ray sigma outside object
(crgrow =
0.) Cosmic ray grow radius
(objval =
0) Mask value for objects
(lobjsig=
10.) Low object detection sigma
(hobjsig=
5.5) High object detection sigma
(objgrow=
0.) Object grow radius
(mode =
ql)
:wq
NOTE: Write this param file DO NOT RUN it.
NOTE: IN THE DOUBLE SET OF NUMBERS, THE LEFT SET TAKES APPROXIMATELY
3 TIMES AS LONG TO RUN BUT DOESN'T
SEEM SIGNIFICANTLY BETTER
- To actually run craverage, it has to be run on individual
frames inside
each MOSAIC frame. Therefore write a script, following
these steps
a) Make a list with all of the extensions of all of the object files.
To do this, epar mscstat with the following parameters
epar mscstat
PACKAGE = mscred
TASK = mscstat
images =
obj*fits Images
(extname=
) Extension name selection
(gmode =
no) Global mode statistics?
(fields =
image) Fields to be printed
(lower =
INDEF) Lower cutoff for pixel values
(upper =
INDEF) Upper cutoff for pixel values
(binwidt=
0.1) Bin width of histogram in sigma
(format =
yes) Format output and print column labels?
(fd1 =
)
(fd2 =
)
(mode =
ql)
:wq
DO NOT RUN THIS YET. Instead :wq and run it from the command line
redirecting the output into a temporary file
mscstat > temp1
- edit the temp1 file so it only contains the file names (remove
the
header and the spaces) and replace the "im" in the brackets
with
nothing, eg. "[im1]" becomes "[1]".
- copy temp1 to llobj for later use and copy temp1 to temp2
cp temp1 llobj
cp temp1 temp2
- edit temp2 to create the crmask names. Replace 'obj' with
'crmask'.
Replace the '.fits[' with '_'. Replace the ']' with
nothing.
Add the ' output="" crmask="' to the beginning of each line.
And
finally add the '.pl"' to the end of each line.
- Paste temp1 and temp2 together into crscript.cl
!paste temp1 temp2 > crscript.cl
- Edit crscript.cl replacing the TAB with nothing and putting a "craverage
"
at the beginning of each line. The final result should
look similar
to this:
craverage obj042.fits[1] output="" crmask="crmask042_1.pl"
craverage obj042.fits[2] output="" crmask="crmask042_2.pl"
craverage obj042.fits[3] output="" crmask="crmask042_3.pl"
craverage obj042.fits[4] output="" crmask="crmask042_4.pl"
craverage obj042.fits[5] output="" crmask="crmask042_5.pl"
craverage obj042.fits[6] output="" crmask="crmask042_6.pl"
craverage obj042.fits[7] output="" crmask="crmask042_7.pl"
craverage obj042.fits[8] output="" crmask="crmask042_8.pl"
craverage obj043.fits[1] output="" crmask="crmask043_1.pl"
.
.
.
craverage obj147.fits[6] output="" crmask="crmask147_6.pl"
craverage obj147.fits[7] output="" crmask="crmask147_7.pl"
craverage obj147.fits[8] output="" crmask="crmask147_8.pl"
then run the script
cl < crscript.cl
Each mask takes approximately 3-4 minutes to generate on an
unencumbered
Sparc Ultra-60, so don't be alarmed it there is no output immediately from
this script. Give it a few minutes. An entire night will take
approximately
2 days (that's right, 48 hours) to run!
- Once all of the cosmic ray masks are generated, they need to be
"added" to the corresponding bad pixel mask. This is
done with
the script crplusbpmask.cl.
Before first running this script, copy it to your scripts$
directory,
add the line "task crplusbpmask=scripts$crplusbpmask.cl" to
your
loginuser.cl or login.cl file, and restart iraf.
The script takes as input the file numbers of the object files,
so
create a file (llnum) which contains the three digit numbers
from
the obj*fits files.
ls obj*fits > llnum
Edit llnum so it only contains the three digit numbers, one
per line,
from the obj* files. Note that leading zeros are important
to retain.
The script adds the bad pixels from the cosmic ray masks to
the list of bad
pixels in the bad pixel file and updates the headers in both
the
bpm files and the obj files.
epar crplusbpmask
PACKAGE = user
TASK = crplusbpmask
image =
@llnum Input image(s)
(inimgli=
)
(mode =
ql)
:g
13) Creating a Correct WCS
- The 0.9m writes a preliminary WCS into the headers of each file,
but
due to pointing errors and image distortion across the field,
this
preliminary coordinated system needs to be updated to a more
accurate
version before it can be used scientifically.
- Create a list of all the object frames for each object observed
during
the night including standards.
ls obj* > llm36
ls obj* > llm37
ls obj* > llm67
ls obj* > llsa107
and edit each ll* list so it only contains the images of that
object.
If a single object is observed in multiple bands or with varying
exposure
times ALL of the exposures should be grouped in one list file.
In other
words you do not need separate lists for each band or exposure
time
only for each object.
- Now get a coordinate list for each object.
epar mscgetcatalog
input =
@llm36 List of Mosaic files
output = catm36.usno
Output file of sources
(magmin =
9.) Minimum magnitude
(magmax =
18.) Maximum magnitude
(catalog= NOAO:USNO-A2)
Catalog
(rmin =
45.) Minimum radius (arcmin)
(mode =
ql)
:g
- The minimum and maximum magnitude should be set to avoid including
saturated stars and stars not well enough sampled to have
good centroids.
Since this is a function of band, exposure time, and transparency,
it
is difficult to know this ahead of time. So at this
point use the
values above to get a set of stars good enough for the entire
data
set. These catalogs will be edited later so they only
include valid
stars.
- Set rmin to a large enough number that it will cover the entire
set of
exposures for the object/field. If while running msccmatch
points
are not sampled across the entire image, rebuild your catalog
with
rmin set to a larger number.
- Note that if the saturation level was set correctly in the headers,
and all the bad pixel masks were generated like described
above,
saturated stars should automatically be excluded my msccmatch
below.
- Catalogs besides the USNO-A2 catalog can be used. To get
a list
of all catalogs accessible type:
aclist * verbose-
All of these notes assume that you used NOAO:USNO-A2.
If you use
a different catalog, column order may be different, and the
sorting
will have to be adjusted accordingly.
- Create a sub-catalog
The retrieved catalog likely goes do dim, and if the target
is in or
near the galactic plane, may contain many thousands of stars.
Since
the matching does not take into account magnitude, too many
stars may
lead to wrong matches, or may confuse the gross off-set step
and lead
to no match at all. Therefore, create a sub-catalog
with only the
brightest brightest 1500 - 3000 stars.
!cat catm36.usno | sort -n +2.0 -3.0 > subcatm36.usno
Then edit the subcatalog so it only contains the magnitudes
< 14
(or some other bright magnitude so that you get somewhere
around 1500
stars total). If you include too many dim stars, an
error claiming
that "too many stars failed to centroid" will be encountered,
and the fit
will not be done. It may be necessary to create more
than one sub-catalog
file for a set of observations, especially if U band is present
or exposure
times of widely varying length.
NOTE: The catalog files contain 4 columns, RA, DEC,
red_mag, blue_mag.
You can use either magnitude for the sort but the above command
assumes
you use the value in the third column.
- Match image stars and catalog positions:
epar msccmatch
PACKAGE = mscred
TASK = msccmatch
input =
obj042 List of input mosaic exposures
coords = subcatm36.usno
Coordinate file (ra/dec)
(outcoor=
) List of updated coordinate files
(usebpm =
yes) Use bad pixel masks?
(verbose=
yes) Verbose?
# Coarse Search
(nsearch=
60) Maximum number of positions to use in search
(search =
120.) Translation search radius (arcsec)
(rsearch=
0.2) Rotation search radius (deg)
# Fine Centroiding
(cbox =
11) Centering box (pixels)
(maxshif=
7.) Maximum centering shift to accept (arcsec)
(csig =
0.1) Maximum centering uncertainty to accept (arcsec)
(cfrac =
0.5) Minimum fraction of accepted centers
(listcoo=
no) List centered coordinates in verbose mode?
# WCS Fitting
(nfit =
50) Min for fit (>0) or max not found (<=0)
(rms =
1.) Maximum fit RMS to accept (arcsec)
(fitgeom=
general) Fitting geometry
(reject =
3.) Fitting rejection limit (sigma)
(update =
yes) Update coordinate systems?
(interac=
yes) Interactive?
(fit =
yes) Interactive fitting?
(graphic=
stdgraph) Graphics device
(cursor =
) Graphics cursor
accept =
no Accept solution?
(mode =
ql)
:g
- input can be an @ll* list as well as a single file name.
- nsearch can be increased if no coarse solution is found,
but total
time to do the search is very sensitive to this
number.
- search is the translation (simple offset) tolerance in doing
the
coarse search. If the translation is large
(eg. large pointing offset)
then this number may have to be increased.
NOTE: the number is in
ARCSEC not pixels, just like most the other parameters.
- reject is the tolerance parameter (sigma clip) used in doing
the final
fit. If you get a generally good fit its
okay to decrease this number
so you don't have to manually delete many outliers.
However,
if there is a lot of scatter in the fit, it may
be wise to use a
larger number like 3.0. In either case you
may need to set MAXITER
to 1 interactively, i.e., use the colon command,
:maxiter 1
If you decide you want to run this non-iteractively,
or you just
get tired of typing ":maxiter 1", the maxiter
parameter can be
specified in geomap. That is "epar geomap"
and set maxiter to 1.
This is the only way to set this parameter.
It cannot be set,
other than temporarily with the colon command,
from within msccmatch.
- use the standard IRAF interactive fitting (x,r,y,s,g) deleting
(d)
any outliers and re-fitting (f). You should
be able to get the
residuals down to ~0.3 - 0.4 arcsec if everything
is working.
Quit (q) from one of the fit windows and answer
yes (yes) to
accepting the fit, and the new WCS will be written
to the header.
NOTE: The graphics window can be increased
in size which will make
picking individual stars for deletion easier.
NOTE: If there is a cluster in the field,
its proper motion may
cause all of its members to be outliers in the
fit.
NOTE: If you know everything is working okay
for this run, set
input=@llm36 (input list of all the exposures for an
object),
reject=1.0, interac=no, fit=no, and let the process
run. Inspect
the output occasionally (or redirect it to a file)
to be sure that
all of the images processed correctly.
- Repeat the msccmatch procedure for every frame/object.
14) Projecting the images onto a standard plate scale.
Due to optical distortions and the gaps between the 8 chips,
MOSAIC
images must be projected onto a standard coordinate system and a uniform
plate scale.
The first step is to use proto.fixpix to replace the values
of bad
pixels with values interpolated from neighboring pixels. This will
avoid
"ringing" when using the sinc interpolation in the next step.
proto
epar fixpix
PACKAGE = proto
TASK = fixpix
images =
@llobj List of images to be fixed
masks =
BPM List of bad pixel masks
(linterp=
INDEF) Mask values for line interpolation
(cinterp=
INDEF) Mask values for column interpolation
(verbose=
yes) Verbose output?
(pixels =
no) List pixels?
(mode =
ql)
:g
- The input file, llobj, was created in the "finding cosmic
rays" step,
and is a list of all obj* file names with extensions.
- Now that the object frames are ready to be projected, do
the projection
with the script improject.cl.
This script uses iterstat.cl to determine
the "sky" level. This value is then used in mscimage
as the "edge
value" and mscimage is run to project the images. Both
of these routines
must be set up as scripts (see crplusbpmask above in step
12 for
instructions on how to set up scripts.)
- Each object/field must be projected in its own group because
the
projection is calculated relative to the FIRST image
in the image
list (the "reference" keyword is left blank).
If the mosaic you want
to ultimately create is large on the sky, you may want
to put the
most central image as the first image in the list but
if you are
just interested in a star cluster or galaxy, any order
will work.
The list files created for the msccmatch step will work,
just be
sure they are still complete.
then, for each list,
epar improject
PACKAGE = user
TASK = improject
list =
@llm36 list of object files -- without extensions
(pref =
mos) prefix for projected images
(imglist=
)
(mode =
ql)
:g
- If the script crashes after doing the "Resampling" claiming
that
"combine" doesn't have a param file enter the following
commands:
flpr
flpr
mscred
epar combine
(:wq out of the
epar)
Remove the mos* files that were created for the image
it crashed
on, and any tmp* files in the directory.
If this wasn't the first image in the input list, create
at new
input list removing all inputs that ran correctly.
Rerun the script.
- this script takes about 15 minutes per image, so an entire
night of
object frames may take more than one day (24 hours)
to run.
- after the script is complete, examine the projected images
to be sure
everything looks okay.
- to reduce all of the objects for an entire night, you can create
a script that calls improject for each object, for example
using
mkscript.
15) Writing the data to tape
- Look through the directory and remove any temporary files
you made
during processing. Keep the lists (ll*) and any
command files you
wrote (*.cl) because if it is necessary to reprocess
this night for
any reason, these will come in handy.
Make a directory one level up called Skys and copy all
of the sflat*
and Sflat* files there. If on subsequent nights
of the run sky flats
were not taken, or insufficient sky flats were taken,
you can use
these. NOTE: one should reduce the night during
the run with the
most/best sky flats first if there is a single night
with a complete
set.
Remove the zero*.fits, sflat*.fits, and dflat*.fits
images. These
images have little processing and can easily be recreated
from the
raw images.
Be sure to remove the Raw directory unless you want
to make another
copy of the raw files for some reason.
- Write the data to tape using:
tar -cv .
if the TAPE environment variable is set correctly.
If not use,
tar -cvf /dev/rmt/xxx .
where xxx is the device string.
NOTE: The data sizes for one night of MOSAIC data
can be huge
so it my be necessary to break the data into subsets
and write
each subset to a different tape.
For example: Create a directory above the currently
level (N1MOS)
and move all of the mos* files to it. In the 00feb/N1
example,
this left 31.0 GB in the N1 directory and put 23.4 GB
in the N1MOS
directory. Each of these should fit on one DLT.
07/21/2003