; NAME:
;       pro mkgdmask.pro (v1.2)
;
; PURPOSE:
;       This program makes good pixel (good=1 & bad=0) masks following a recipe sent
;       by U. Hawaii. The first step is to find and mark "bad reset" pixels. These
;       are identified as statistical outliers, they may be either high or low. The
;       step is to find high dark current pixels. These are pixels that have dark current
;       greater than 0.1 e-/s
;
; CALLING SEQUENCE:
;       mask = mkgdmask('dark_cube.fits' [, /help [, /cleanup]])
;
; INPUTS:
;       dark_cube.fits   - A FITS data cube of a dark image
;
; KEYWORD PARAMETERS:
;       help             - Print a help message
;       cleanup          - Erase the calibrated dark image on exit
;
; EXAMPLE
;       Make a mask for a detector keeping the calibrated dark image that is
;       is produced as part of the analysis. Return the mask in a 2-dimensional
;       integer array named gdmask.
;
;       gdmask = mkgdmask('dark_image.fits')
;
;
; REFERENCE:
;
;
; MODIFICATION HISTORY:
;       Written by:
;       Modified:    B.J. Rauscher, IDTL, 9 April 2003
;                    Initial release (v1.1)
;
;                    B.J. Rauscher, 30 April 2003
;                    Updated to use scaName keyword when calling
;                    refsub (v1.2)
;
;-

function mkgdmask, dark_cube, help=help, cleanup=cleanup

       ;; Define constants
       bad_rst_threshold = 3.0; Sigma threshold for identifying bad reset pixels
       high_dk_threshold = 0.1; (ADU/s) Absolute threshold for identifying high dark current pixels

       ;; Print help message if requested
       if (keyword_set(help)) then begin
           stop, 'Usage: mask = mkgdmask(dark.fits [, /help [, /cleanup]])'
       endif else begin
           help = 0
       endelse

       ;; Set command line defaults
       if (not keyword_set(cleanup)) then cleanup = 0

       ;; Read in the fits header so that we can access needed
       ;; information
       print, 'Reading in FITS header'
       fits_read, dark_cube, data, header, /header_only
       xsize = sxpar(header, 'NAXIS1')
       ysize = sxpar(header, 'NAXIS2')
       zsize = sxpar(header, 'NAXIS3')
       exp_time = sxpar(header, 'EXP_TIME')
       scaName = sxpar(header, 'DETNAME')
       scatype = ''
       getscatype, header, scatype

       ;; Create an array to hold the mask
       gdmask = intarr(xsize, ysize)
       gdmask[*,*] = 1; Initialize all pixels to =1 (good)

       ;; Define the active portion of the SCA
       if (scatype eq 'H1RG' or scatype eq 'H1R') then begin
           x1 = 4
           x2 = 1019
           y1 = 4
           y2 = 1019
       endif
       if (scatype eq 'H2RG') then begin
           x1 = 4
           x2 = 2043
           y1 = 4
           y2 = 2043
       endif
       if (scatype eq 'SB304') then begin
           x1 = 4
           x2 = 2051
           y1 = 4
           y2 = 2047
       endif

       ;; Mask reference pixels so that they are not included in statistics
       for y = 0, ysize - 1 do begin
           for x = 0, xsize - 1 do begin
               if (x lt x1) then gdmask[x,y] = 0
               if (x gt x2) then gdmask[x,y] = 0
               if (y lt y1) then gdmask[x,y] = 0
               if (y gt y2) then gdmask[x,y] = 0
           endfor
       endfor

       ;; Identify "bad reset" pixels by forming a CDS image
       ;; and looking for low statistical outliers.
       print, 'Finding bad reset pixels'
       frame0 = getframe(dark_cube, 0)
       frame1 = getframe(dark_cube, 1)
       cds_difference = frame1 - frame0
       refsub, cds_difference, scaName=scaName
       djs_iterstat, cds_difference[x1:x2,y1:y2], sigrej=bad_rst_threshold, mask=_mask
       gdmask[x1:x2,y1:y2] = _mask[*,*]

       ;; Identify high dark current pixels by forming the dark
       ;; current image. This creates an image with the suffix
       ;; '.cal.fits'. Set the /cleanup keyword to delete this.
       print, 'Finding high dark current pixels'
       mkmulti, dark_cube; mkmulti forms the dark current image and does reference pixel correction

       cal_file = strmid(dark_cube, 0, strpos(dark_cube, '.fits', /REVERSE_SEARCH))+'.cal.fits'
       fits_read, cal_file, dark_data, dark_header

       ;; Flag high dark current pixels
       dark_flux = dark_data / exp_time; Convert to ADU/s

       ;; Find and mark high dark current pixels
       for y = y1, y2 do begin
           for x = x1, x2 do begin
               if (dark_flux[x,y] ge high_dk_threshold) then begin
                   gdmask[x,y] = 0
               endif
           endfor
       endfor

       ;; Clean up if requested
       if (cleanup) then file_delete, cal_file

       ;; Print percentage of good pixels
       print, 'Fraction of good pixels = '+strtrim(string(mean(gdmask)),2)

       ;; Return to caller
       return, gdmask

end