Cookbook for 0.9m MOSAIC Reduction of

WOCS Data

Written by:  Glenn Tiede

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