; $Id: get_rdnoise_new.pro,v 1.3 2003/09/11 18:09:12 kgordon Exp $
;
;+
; NAME:
;	GET_RDNOISE_NEW
;
; PURPOSE:
;	This procedure determines the read noise on the MIPS 24, 70 &
;	160 micron arrays.  This is done by computing the noise from a
;	stack of 50 DCEs and subtracting the expected shot noise to
;	generate the read noise on a per pixel basis.  The mips_sloper
;	portion of the DAT must be run to provide input to this
;	procedure.  Any DCEs in which a pixel was affected by a cosmic
;	ray are not used to compute the read nosie for that pixel.
;	The responsivity drift of each pixel is removed by fitting a
;	polynomial (quadratic as default) and subtracting the fit from
;	that pixel's data.
;
; CATEGORY:
;	MIPS utilities.
;
; CALLING SEQUENCE:
;	GET_RDNOISE_NEW, Filebase
;
; INPUTS:
;	Filebase: Base filename of a mips_sloper .red.fits file.  This
; 	          file must be a 50 DCE dark (no stim flashes) dataset
; 	          for valid read noise measurements to be produced.
; 	          No electronic linearity correction should be applied
; 	          when running mips_sloper.
;
; KEYWORD PARAMETERS:
;       A160_ENG: Determine the read noise from the 160 micron
;                 engineering pixels.  Has no effect if input data is
;                 not for the 160 micron array.
;
;       EXPTIME: Override the header exptime.  This is necessary for
;                old data which does not have an exposure time in the
;                header or old 24 micron SUR data from which the
;                exptime cannot be determined from the data.
;
;       MAX_DCE: Maximum DCE to use, all DCEs with numbers greater
;                than this are rejected.  
;
;       N_DEC_REJ: Number of initial DCEs to reject completely.
;
;       NDEGREE: Degree of polynomial to use to remove the
;                responsivity drift of the data.  Default is NDEGREE=2
;                -> a quadratic.
;
;	NO_CR:	Do not reject cosmic rays (id=64) from read noise
;	        calculation.  The default is to reject them.
;
;       PNG: Save the three screen plots (counts vs DCE, # pixels vs
;            read noise, and # pixels vs counts) to PNG image files.
;
;
; OUTPUTS:
;	There are 4 files output from this routine.  
;           filebase.rdn.report : The array average read noise
;                                 measurements.
;           filebase.rdn.rdn.fits : An image of the read noise as a
;                                   function on a pixel-by-pixel
;                                   basis. 
;           filebase.rdn.noise.fits : An image of the noise as a
;                                     function on a pixel-by-pixel
;                                     basis. 
;           filebase.rdn.cnt.fits : An image of the counts as a
;                                   function on a pixel-by-pixel
;                                   basis. 
;
; OPTIONAL OUTPUTS:
;	If the PNG keyword is set, then the 3 PNG images file are
;	created as outlined in the PNG keyword description.
;
; EXAMPLE:
;	Determine the read noise for the 70 micron array using CTA
;	data.  This 50 DCE dataset must first be run through
;	mips_sloper (turning off the linearity correction).  After
;	than, type the following on the IDL commandline.
;
;		get_rdnoise,'mips_2001-227T23.31.57.748_0_A70'
;
; MODIFICATION HISTORY:
; 	Written by:	Karl Gordon (1999-2003)
;	Oct, 2002	Documentation header added to existing program
;	                This program was originally written in 1999
;	                and has been used since to compute the read
;	                noise measurements for all the MIPS arrays
;	                (flight, char, lab, etc).
;       Mar, 2003       Finished the documentation
;                       Added the ability to handle NaNs.
;                       Only does 160 micron data in the 20x3 format.
;       Oct, 2003 James M. - Updated the calculation for 24 micron SUR
;                            onflight data.  Basically, added
;                            rejection of boost frames and sigma
;                            rejection during calculation of noise.
;-
; ----------------------------------------------------------------------

@kplot
@comp_std
@scrn2png

PRO get_rdnoise_new,filebase,a160_eng=a160_eng,exptime=exptime, $
                max_dce=max_dce,n_dce_rej=n_dce_rej, $
                ndegree=ndegree,no_cr=no_cr,png=png,bad70=bad70, $
                empir_cr=empir_cr

if (keyword_set(empir_cr)) then begin
    if (empir_cr EQ 1) then empir_cr = 4.0
endif                

device,decomposed=0

if (not keyword_set(ndegree)) then ndegree = 2

; get the reduced data

fits_open,filebase+'.red.fits',fcb

fits_read,fcb,trash,header,exten_no=0  ; get the main header for the file
header0 = header

if (keyword_set(max_dce)) then fcb.nextend = min([max_dce,fcb.nextend])

if (not keyword_set(n_dce_rej)) then n_dce_rej = 0
n_dce_begin = 1 + n_dce_rej
n_exp = fcb.nextend - n_dce_rej

inb=0
oldboost=0
for i = n_dce_begin, fcb.nextend do begin
    fits_read,fcb,cube,header,exten_no=i
    size_cube = size(cube)
    if (size_cube[1] eq 128) then boost = sxpar(header,"BOOSTFRM") else boost = 0
; if 24um, determine whether DCE had a bias boost;
; do not use the boost DCE OR the one immediately after it in the noise calc.
    if (boost eq 0 and oldboost eq 0) then begin
        if (inb EQ 0) then begin
            if (size_cube[1] EQ 20) then begin ; handle 160 micron data
                if (keyword_set(a160_eng)) then begin
                    size_cube[2] = 1
                endif else begin
                    size_cube[2] = 2
                endelse
            endif
            meta_cube = fltarr(size_cube[1],size_cube[2],size_cube[3],n_exp)
            if (size_cube[1] EQ 128) then begin
                shot_factor = 1.0
            endif else begin
                shot_factor = 2.0
            endelse
        endif
        if (size_cube[1] EQ 128) then begin
            meta_cube[*,*,*,inb] = cube
        endif else if (size_cube[1] EQ 20) then begin ; handle 160 array
            if (keyword_set(a160_eng)) then begin
                meta_cube[*,*,*,i-1-n_dce_rej] = cube[*,1,*]
            endif else begin
                meta_cube[*,0,*,i-1-n_dce_rej] = cube[*,0,*]
                meta_cube[*,1,*,i-1-n_dce_rej] = cube[*,2,*]
            endelse
        endif else begin        ; do the other arrays
            meta_cube[*,*,*,i-1-n_dce_rej] = cube
            if (keyword_set(bad70)) then begin
                meta_cube[0:15,*,0,i-1-n_dce_rej] = 0.0
                meta_cube[16:19,0:7,0,i-1-n_dce_rej] = 0.0
            endif
        endelse
        if (size_cube[1] eq 128) then begin
; increment indices for DCE 1,2 rejection for 24um
            inb = inb+1
            oldboost = 0
        endif else begin
            inb = 1
        endelse
    endif
    if (size_cube[1] eq 128 and boost eq 1) then oldboost = 1 else oldboost=0
endfor
fits_close,fcb
n_noboost = inb  ;number of non-boost,boost+1 DCEs

gain = fxpar(header0,"GAIN")

; 24 micron SUR exptime fix

if (not keyword_set(exptime)) then begin
    exptime = fxpar(header0,"EXPTIME")
endif else begin
    exptime = exptime*1.048576
endelse

x_npts = size_cube[1]
y_npts = size_cube[2]
n_planes = size_cube[3]
n_pixels = x_npts*y_npts
if (x_npts ne 128) then begin
  dncube = fltarr(x_npts,y_npts,n_planes,n_exp)
  dncube[*,*,0,*] = exptime*meta_cube[*,*,0,*] 
  dncube[*,*,1:n_planes-1,*] = meta_cube[*,*,1:n_planes-1,*]
endif else begin
; for 24um, get rid of the boost, boost+1 DCEs
  dncube = fltarr(x_npts,y_npts,n_planes,n_noboost)
  dncube[*,*,0,*] = exptime*meta_cube[*,*,0,0:n_noboost-1]
  dncube[*,*,1:n_planes-1,*] = meta_cube[*,*,1:n_planes-1,0:n_noboost-1]
  n_exp = n_noboost
endelse

; setup the 160 array correctly (old method)

min_x = 0
max_x = x_npts - 1

exp_vals = findgen(n_exp) + 1
ave_counts = fltarr(n_exp)
for i = 0,(n_exp-1) do begin
    n_good_pixels = total(finite(dncube[*,*,0,i]))
    ave_counts[i] = gain*total(dncube[*,*,0,i],/nan)/n_good_pixels/exptime
;    print,i,n_good_pixels,ave_counts[i],gain,total(dncube[*,*,0,i],/nan)
endfor
kplot,exp_vals,ave_counts,psym=100, $;,yrange=[1.9e4,2e4],xrange=[1,60], $
  kplot_type='ii',xtitle='DCE',ytitle='counts [e/sec]',color=1,background=255


count_image = fltarr(x_npts,y_npts)
rdnoise_image = fltarr(x_npts,y_npts)
noise_image = fltarr(x_npts,y_npts)
ave_npts_image = fltarr(x_npts,y_npts)
ave_npts_used = 0.0
n_shot_noise = 0
ave_shot_noise = 0.0
cos_val = 64
sigma24 = fltarr(n_exp)
for i = 0,(x_npts-1) do begin
    for j = 0,(y_npts-1) do begin 
        ; reject NaNs
        ifin = where(finite(dncube[i,j,0,*]),n_ifin)
	if (n_ifin LT 1) then begin
            print,'Pixel ',i+1,j+1,' :  all DCEs rejected'
	    indxs = ifin
	    n_indxs = n_ifin
        endif else begin
            ; reject cosmic rays
;            if (x_npts NE 128) then begin
            if ((x_npts NE 128) AND (not keyword_set(empir_cr))) then begin
                indxs2 = where((fix(dncube[i,j,2,ifin]/cos_val) mod 2) NE 1,n_indxs)
                if (n_indxs gt 0) then indxs = ifin[indxs2]
            endif else begin
                ; for 24um SUR, reject CRs by sigma rejection
                mslope = median(dncube[i,j,0,ifin])
                ; calculate standard deviation from the median over all DCEs,
                ; then determine sigma for each DCE
                sdev24 = sqrt(total((dncube[i,j,0,ifin]-mslope)^2)/float(n_ifin-1)) 
                sigma24 = (dncube[i,j,0,ifin]-mslope)/sdev24 
                ; include only DCEs with sigma < 4
                indxs = where(sigma24 lt empir_cr,n_indxs)
                if (n_indxs GT 0) then indxs = ifin[indxs]
            endelse
            if (keyword_set(no_cr)) then begin
                indxs = where(finite(dncube[i,j,0,*]),n_indxs)
            endif
            ave_npts_used = ave_npts_used + n_indxs
            ave_npts_image[i,j] = n_indxs
            
            ; compute average counts [total electrons] for each pixel
            if (n_indxs GT 0) then begin
                counts = gain*total(dncube[i,j,0,indxs])/float(n_indxs)
            endif else begin
                counts = 0.0
            endelse
            count_image[i,j] = counts
            
; compute the read noise
;    std_dev of counts - shot noise predicted from counts (sqrt)

; fit a line to the slopes to remove a gradual change in slope values        

            if (n_indxs GT ndegree+2) then begin
                fit_coeff = poly_fit(exp_vals[indxs],gain*dncube[i,j,0,indxs],ndegree,calc_slope)
                std_dev = sqrt(total((gain*dncube[i,j,0,indxs]-calc_slope)^2)/float(n_indxs-ndegree-1))
            endif else begin
                std_dev = 0.0
            endelse

            noise_image[i,j] = std_dev
            
            if (counts GT 0.0) then begin 
                shot_noise = sqrt(shot_factor*counts)
            endif else begin
                shot_noise = 0.0
            endelse 
            if (shot_noise LT std_dev) then begin
                rdnoise_image[i,j] = sqrt(std_dev^2 - shot_noise^2)
                ave_shot_noise = ave_shot_noise + shot_noise
                n_shot_noise = n_shot_noise + 1
            endif else begin
                rdnoise_image[i,j] = 0.0
;                print,'shot-noise dominates in pixel',i+1,j+1
;                print,'computed shot noise = ',shot_noise
;                print,'total noise         = ',std_dev
;                print,'# of good DCEs      = ',n_indxs
            endelse
        endelse

    endfor
endfor 
ave_npts_used = ave_npts_used/n_pixels
ave_shot_noise = ave_shot_noise/n_shot_noise

; output counts and rdnoise arrays

fxhmake,header0,count_image
fxwrite,filebase+'.rdn.cnt.fits',header0,count_image

fxhmake,header0,rdnoise_image
fxwrite,filebase+'.rdn.rdn.fits',header0,rdnoise_image

fxhmake,header0,noise_image
fxwrite,filebase+'.rdn.noise.fits',header0,noise_image

; setup output file for readout report

openw,unit1,filebase+'.rdn.report',/get_lun
printf,unit1,filebase
printf,unit1,'gain = ' + strtrim(string(gain,format="(F5.2)"),2) + ' e/DN'
printf,unit1,'exptime = ' + strtrim(string(exptime,format="(F7.4)"),2) + ' sec'
printf,unit1,'# DCEs reject = ' + strtrim(string(n_dce_rej),2)
printf,unit1,'counts/sec = electrons/sec'
printf,unit1,"average # exposures used = " +  $
  strtrim(string(ave_npts_used,format="(F5.2)"),2) + " per pixel"
printf,unit1,"average shot noise = " + $
  strtrim(string(ave_shot_noise,format="(F7.2)"),2) + " (" + $
  strtrim(string(n_shot_noise),2) + ")"

print,'gain = ' + strtrim(string(gain,format="(F5.2)"),2) + ' e/DN'
print,'exptime = ' + strtrim(string(exptime,format="(F7.4)"),2) + ' sec'
print,'# DCEs reject = ' + strtrim(string(n_dce_rej),2)
printf,unit1,'counts/sec = electrons/sec'
print,"average # exposures used = " + $
  strtrim(string(ave_npts_used,format="(F5.2)"),2) + " per pixel"
print,"average shot noise = " + $
  strtrim(string(ave_shot_noise,format="(F7.2)"),2) + " (" + $
  strtrim(string(n_shot_noise),2) + ")"

count_image = count_image/exptime

; do statistics on the array basis

ar_counts_std = comp_std(count_image[min_x:max_x,*],ave=ave_val)
ar_counts = ave_val
ar_rdnoise_std = comp_std(rdnoise_image[min_x:max_x,*],ave=ave_val)
ar_rdnoise = ave_val

printf,unit1,''
printf,unit1,'***array statistics***'
printf,unit1,'        counts/sec            read noise'
printf,unit1,ar_counts,ar_counts_std,ar_rdnoise,ar_rdnoise_std, $
  format="(4x,F8.1,' +/- ',F8.1,2x,F8.1,' +/- ',F8.1)"

print,'***array statistics***'
print,'        counts/sec            read noise'
print,ar_counts,ar_counts_std,ar_rdnoise,ar_rdnoise_std, $
  format="(4x,F8.1,' +/- ',F8.1,2x,F8.1,' +/- ',F8.1)"

; get the 99% level readnoise and std deviation

rdnoise_image = rdnoise_image[min_x:max_x,*]
sindxs = sort(rdnoise_image)
npts = n_elements(sindxs)

if (max_x LT 32) then top_k = 0.99*npts else top_k = 0.95*npts

ar_counts_std = comp_std(count_image[sindxs[0:top_k-1]],ave=ave_val)
ar_counts = ave_val
ar_rdnoise_std = comp_std(rdnoise_image[sindxs[0:top_k-1]],ave=ave_val)
ar_rdnoise = ave_val

if (max_x LT 32) then begin
    print,'***array statistics (95% w/ zeros)***'
endif else begin
    print,'***array statistics (99% w/ zeros)***'
endelse
print,'        counts/sec            read noise'
print,ar_counts,ar_counts_std,ar_rdnoise,ar_rdnoise_std, $
  format="(4x,F8.1,' +/- ',F8.1,2x,F8.1,' +/- ',F8.1)"

if (max_x LT 32) then begin
    printf,unit1,'***array statistics (95% w/ zeros)***'
endif else begin
    printf,unit1,'***array statistics (99% w/ zeros)***'
endelse
printf,unit1,'        counts/sec            read noise'
printf,unit1,ar_counts,ar_counts_std,ar_rdnoise,ar_rdnoise_std, $
  format="(4x,F8.1,' +/- ',F8.1,2x,F8.1,' +/- ',F8.1)"

; get the 99% level readnoise and std deviation (and no zeros)

sindxs = sort(rdnoise_image)
npts = n_elements(sindxs)

indxs = sindxs[0:top_k-1]

indxs2 = where(rdnoise_image[indxs] NE 0.0,n_indxs2)
if (n_indxs2 GT 0) then indxs = indxs[indxs2]

ar_counts_std = comp_std(count_image[indxs],ave=ave_val)
ar_counts = ave_val
ar_rdnoise_std = comp_std(rdnoise_image[indxs],ave=ave_val)
ar_rdnoise = ave_val

if (max_x LT 32) then begin
    print,'***array statistics (95% w/o zeros)***'
endif else begin
    print,'***array statistics (99% w/o zeros)***'
endelse
print,'# pixels used: ',n_indxs2
print,'# DCEs/pixel = ',total(ave_npts_image[indxs])/n_indxs2
print,'        counts/sec            read noise'
print,ar_counts,ar_counts_std,ar_rdnoise,ar_rdnoise_std, $
  format="(4x,F8.1,' +/- ',F8.1,2x,F8.1,' +/- ',F8.1)"

if (max_x LT 32) then begin
    printf,unit1,'***array statistics (95% w/o zeros)***'
endif else begin
    printf,unit1,'***array statistics (99% w/o zeros)***'
endelse
printf,unit1,'# pixels used: ',n_indxs2
printf,unit1,'        counts/sec            read noise'
printf,unit1,ar_counts,ar_counts_std,ar_rdnoise,ar_rdnoise_std, $
  format="(4x,F8.1,' +/- ',F8.1,2x,F8.1,' +/- ',F8.1)"

free_lun,unit1

; compute histogram

if (keyword_set(png)) then begin
    scrn2png,filebase+'_tot.png'
endif else begin
    ans = ''
    read,'Continue: ',ans
endelse

kxrange = krange(rdnoise_image)
min_histo = kxrange[0]
max_histo = kxrange[1]
npts_histo = 100
binsize = (max_histo - min_histo)/npts_histo
npts_histo = (max_histo - min_histo)/binsize
val_bins = min_histo + (findgen(npts_histo) + 0.5)*binsize
histo = histogram(rdnoise_image,binsize=binsize,min=min_histo,max=max_histo)

kxrange = krange(val_bins)
kxrange[0] = 0.0
kyrange = krange(histo)
kyrange[0] = 0.0
kplot,val_bins,histo,xrange=kxrange,yrange=kyrange, $
  xtitle='read noise [e!U-!N]',ytitle='# pixels', $
  position=[0.15,0.15,0.9,0.9],title=filebase,color=1,background=255

if (keyword_set(png)) then begin
    scrn2png,filebase+'_rdn.png'
endif else begin
    ans = ''
    read,'Continue: ',ans
endelse

; compute histogram

kxrange = krange(count_image[min_x:max_x,*])
min_histo = kxrange[0]
max_histo = kxrange[1]
npts_histo = 100
binsize = (max_histo - min_histo)/npts_histo
npts_histo = (max_histo - min_histo)/binsize
val_bins = min_histo + (findgen(npts_histo) + 0.5)*binsize
histo = histogram(count_image[min_x:max_x,*],binsize=binsize,min=min_histo,max=max_histo)

kxrange = krange(val_bins)
kxrange[0] = 0.0
kyrange = krange(histo)
kyrange[0] = 0.0
kplot,val_bins,histo,xrange=kxrange,yrange=kyrange, $
  xtitle='counts [e!U-!N/sec]',ytitle='# pixels', $
  position=[0.15,0.15,0.9,0.9],title=filebase,color=1,background=255

if (keyword_set(png)) then scrn2png,filebase+'_cnt.png'

END 

; ----------------------------------------------------------------------