; NAME: ; photon2_reduce.pro (v1.5) ; ; PURPOSE: ; This program is used to reduce data taken by the photonxfer_test procedure. ; ; CALLING SEQUENCE: ; photon2_reduce,indir [,outdir[,filelist]] ; ; INPUTS: ; indir - input directory where data and filelist are located ; ; KEYWORD PARAMETERS: ; filelist - list of images to be processed without full path names ; outdir - directory where output will be placed ; refsub - turns on reference pixel correction ; ; EXAMPLE ; Perform photon transfer reduction ; IDL>photon2_reduce, pxdir ; ; REFERENCE: ; ; ; MODIFICATION HISTORY: ; Written by: ; Don Figer, IDTL, April 15, 2003 (v1.0) ; ; Modified by: ; Ernie Morse, IDTL, April 16, 2003 (v1.1) ; Made changes to input & keywords to make consistent with other ; reduction procedures ; Changed im array to long datatype to allow for signal values > 32767 ; Don Figer, IDTL, April 18, 2003 (v1.2) ; Renamed procedure to photon2_reduce, signifying that this is the ; reduction routine for the second type of photon xfer experiment ; we do (with gray scale target). ; Don Figer, IDTL, April 26, 2003 (v1.3) ; Added reference pixel correction, removal of railed pixels, ; removal of hot pixels, and proper calculation of errors on gain and ; readnoise. ; Mike Regan, IDTL, May 2003 (v1.4) ; Made improvements to error calculation ; Ernie Morse, IDTL, February 5, 2004 (v1.5) ; Deleted commented-out portions of code ; Changed name from photon_reduce_mwr to photon_reduce ; Edited to bring code into compliance with IDTL standard. ; Lance Simms, Stanford, May 9, 2007 ; Changed str_sep to strsplit as str_sep is obsolete in IDL 6.3 ;- Changed keyword "WHEEL1" to "FILTER" for RIDL devices code ; Removed 'WHEEL2" keyword since not part of RIDL devices at moment pro photon3_reduce, indir, outdir = outdir, filelist = filelist, $ refsubon = refsubon ; begin by setting defaults for keyword values not provided by the ; calling context if (not keyword_set(outdir)) then begin outdir = indir endif if (not keyword_set(filelist)) then begin filelist = 'file.lst' endif refsubversion = 'dummy' if (not keyword_set(refsubon)) then begin refsubversion = 'None' endif refsubtype = 'None' if (keyword_set(refsubon)) then begin refsubtype = 'Spatial averaging' endif ; determine path directory delimiter from OS type if (!version.os_family eq 'unix') then begin delim = '/' endif else begin delim = '\' endelse ; check to see if input and outdirectory names end in path delimiter--if ; not, add it to the end if (strmid(indir, strlen(indir) - 1, 1) ne delim) then begin indir = indir + delim endif if (strmid(outdir, strlen(outdir)-1, 1) ne delim) then begin outdir = outdir + delim endif ; exit gracefully if filelist does not exist if (file_test(indir + filelist) eq 0) then begin print,'File ' + indir + filelist + ' does not exist' stop endif ; read file list readcol, indir + filelist, allfiles, format = '(A)' ; remove first image from list (this is ad hoc procedure to reduce ; deviant pixels in some grayscale experiments) files = strarr(n_elements(allfiles) - 1) files[*] = allfiles[1:n_elements(allfiles) - 1] ; set histogram bin size binsize = 100. ; define upper limit and lower limit of signal range for use in linear fit lower = 1e3 ; ADU upper = 1e4 ; ADU lowerbin = fix(lower / binsize) upperbin = fix(upper / binsize) ; set flag to control whether the line fit is done to the binned data fittobinned = 1 ; save current device name so it can be reset later dname = !D.NAME ; set plot colors black = 0 white = 255 ; retrieve header information im = readfits(indir + files[0], h) ; Get the date and time from the header getdate, h, theDate gettime, h, theTime ; get additional parameters from the header SCAName = sxpar(h, 'DETNAME') xsize = sxpar(h, 'NAXIS1') ysize = sxpar(h, 'NAXIS2') filter1 = sxpar(h, 'FILTER') ; Read in the 15 read dark ramps exp_time = sxpar(h, 'EXP_TIME') ; extract detector type from the detector name temp = strsplit(SCAName, '-') SCAType = temp[0] ; define a plot label for later use plotlabel = filter1 + '_' $ + strtrim(string(exp_time, format = '(F7.2)'), 2) ; define dark current threshold in ADU/s darkthresh = 10 ; set up arrays zsize = n_elements(files) imcube = fltarr(xsize, ysize, zsize) signalarray = fltarr(xsize, ysize) mean_array = fltarr(xsize, ysize) stdarray = fltarr(xsize, ysize) im = fltarr(xsize, ysize, 2) im1 = fltarr(xsize, ysize, 2) im2 = fltarr(xsize, ysize, 2) ; find dark frames that have same exposure time as photon xfer frames and ; use them to find hot pixels darklist = findfile(indir + 'Dark*15_Reads*.fits', count = numdarks) print, 'numdarks ', numdarks if (numdarks ne 0) then begin ; define array to hold index number of dark images that have proper ; exposure times matchframes = intarr(numdarks) matchframes[*] = 999 ; loop over dark frames to find those that have same exposure time of illuminated frames for i = 0, n_elements(darklist) - 1 do begin ; read in dark image fits_read, darklist[i], darkim, darkheader, /header_only ; proceed if the exposure time of the dark image is the same as ; that of the target image if (sxpar(darkheader, 'EXP_TIME') eq exp_time) then begin matchframes[i] = i print,darklist[i], ' ' , sxpar(darkheader, 'FILTER') endif endfor ; load dark cube newmatchframes = where(matchframes ne 999, matchcount) numdarks = matchcount if (numdarks ne 0) then begin darkcube = lonarr(xsize, ysize, matchcount) for i = 0, matchcount - 1 do begin im = readfits(darklist(newmatchframes[i]), darkheader) im1 = im[*, *, 1] im0 = im[*, *, 0] ; reference pixel correct if (refsubversion ne 'None') then begin refsub, im1, sca = scaname refsub, im0, sca = scaname endif ; remove railed values darkcube[*, *, i] = im1 - im0 endfor print, 'num elements of dark ', matchcount if (matchcount gt 1) then begin ; form median dark image darkarray = median(darkcube, dimension = 3) endif else begin darkarray = darkcube endelse ; identify hot pixels baddarkindices = where(darkarray gt exp_time * darkthresh) endif endif print,'darkarray mean, stddev', avg(darkarray), stddev(darkarray) darkarray = darkcube[*, *, 0] - darkcube[*, *, 0] print,'darkarray mean, stddev', avg(darkarray), stddev(darkarray) ; loop through images to read each one and make a cube for i = 0, zsize - 1 do begin im = readfits(indir + files[i], h) im = float(im) ; replace railed values im(where(abs(im) eq 32767)) = !VALUES.F_NAN im1 = im[*, *, 1] im0 = im[*, *, 0] ; reference pixel correct if (refsubversion ne 'None') then begin refsub, im1, sca = scaname, version = refsubversion refsub, im1, sca = scaname refsub, im0, sca = scaname endif ; replace hot pixels im1 = float(im1) if (numdarks ne 0) then begin im1[baddarkindices] = !VALUES.F_NAN endif im0 = float(im0) if (numdarks ne 0) then begin im0[baddarkindices] = !VALUES.F_NAN endif ; form CDS difference imcube[*, *, i] = im1 - im0 endfor ; make median and mean arrays signalarray = median(imcube, dimension = 3) mean_array = total(imcube, 3) / (zsize) ; make standard deviation array for igroup = 0, zsize - 1 do begin stdarray = stdarray + (imcube[*, *, igroup] - mean_array) ^ 2 endfor stdarray = sqrt(stdarray / (zsize - 1)) print, 'median of sddev array ', median(stdarray) print, 'median of mean array ', median(mean_array) lower = (0.8 * median(mean_array)) lowerbin = fix(lower / binsize) ; make arrays for histograms sighist = histogram(signalarray, reverse_indices = sigindices, min = 0,$ binsize = binsize, max = 1e5) medstdhist = fltarr(n_elements(sighist)) medsighist = fltarr(n_elements(sighist)) points = lonarr(n_elements(sighist)) sigma_signal = fltarr(n_elements(sighist)) sigma_sigma = fltarr(n_elements(sighist)) ; load histograms for i = 0., n_elements(sighist) - 1 do begin ; load arrays if data are present if (sigindices[i] ne sigindices[i + 1]) then begin medstdhist[i] = median(stdarray[sigindices[sigindices[i] : $ sigindices[i + 1] - 1]]) medsighist[i] = median(signalarray[sigindices[sigindices[i] : $ sigindices[i + 1] - 1]]) points[i] = n_elements(stdarray[sigindices[sigindices[i] : $ sigindices[i + 1] - 1]]) endif else begin medstdhist[i] = !VALUES.F_NAN medsighist[i] = !VALUES.F_NAN points[i] = 0 endelse endfor ; calculate the noise in the median and standard deviation sigma_signal = binsize / 3 ; assume sigma is 1/3 of the bin size sigma_sigma = medstdhist / sqrt(2. * (points - 1)) ; calculate line fit to binned or raw data if (fittobinned) then begin ; identify valid data for the fit such that they are valid numbers valid = where(finite(medsighist) eq 1 and finite(medstdhist) eq 1 $ and medstdhist lt 3 * sqrt(upper)) ; constrain valid data to be within fit range valid = where(medsighist(valid) lt upper and medsighist(valid) gt lower) ; set up variables for line fit xfitdata = medsighist(valid) yfitdata = medstdhist(valid) ^ 2 X_SIG = sigma_signal(valid) Y_SIG = sigma_sigma(valid) ^ 2 endif else begin ; identify valid data for the fit such that they are valid numbers ; note that we reject values in the stdarray that are 3 times the upper ; limit of the signal range used in the fitting. the upper limit is in ; ADU, and the '3' is used to do 3-sigma clipping. note that this ; clipping is not the most sensible in that the clipping should be ; based on the signal level for each data point, rather than clipping ; at the upper limit. in any case, it gets rid of horribly deviant ; values valid = where(finite(signalarray) eq 1 and finite(stdarray) eq 1 $ and stdarray lt 3 * sqrt(upper)) ; constrain valid data to be within fit range valid = valid(where(signalarray(valid) lt upper $ and signalarray(valid) gt lower)) ; set up variables for line fit xfitdata = signalarray(valid) yfitdata = stdarray(valid) ^ 2 X_SIG = stdarray(valid) / sqrt(zsize) Y_SIG = (stdarray(valid) / sqrt(2. * zsize)) ^ 2 endelse ; find the readnoise and gain using procedure FITEXY fitexy, xfitdata, yfitdata, rn2, rgain, sigma_A_B, chisq, prob, $ X_SIG = x_sig, Y_SIG = y_sig readnoise = sqrt(rn2) gain = 1. / rgain error_rn = 0.5 * sigma_A_B[0] / readnoise error_gain = sigma_A_B[1] * gain ^ 2 ; Calculate the appropriate y-values for the fitted line fityvals = fltarr(n_elements(xfitdata)) fityvals = readnoise ^ 2 + (xfitdata / gain) ; set up plot titles xtitle = 'Signal (ADU)' ytitle = 'Noise Squared (ADU)' ; make both plots ptype = ['ps', 'jpg'] for k = 0, 1 do begin plottype = ptype[k] if (plottype eq 'ps') then begin ; set up postscript plots set_plot, 'ps' device, file = outdir + 'photonxfer_' + plotlabel + '.ps', $ /landscape ; set filename to be written at bottom of plot fsub = outdir + 'photonxfer_' + plotlabel + '.ps' ; set line and text characteristics for publication quality !P.CHARSIZE = 1.5 !P.CHARTHICK = 4. !X.THICK = 6. !Y.THICK = 6. !P.THICK = 6. !P.TICKLEN = 0.03 radius = 1. filenamecharsize = 1.1 endif if (plottype eq 'jpg') then begin ; set up buffer device for use later in writing JPEG file set_plot, 'z' device, set_resolution = [8000, 6000] device, set_font = 'Courier' ; set filename to be written at bottom of plot fsub = outdir + 'photonxfer_' + plotlabel + '.jpg' ; set line and text characteristics for publication quality !P.CHARSIZE = 10 !P.CHARTHICK = 15. !X.THICK = 20. !Y.THICK = 20. !P.THICK = 15. radius = 8. filenamecharsize = 8. endif ; plot symbols ; define vectors for plotting symbol (circle) ; make a vector of 16 points, A(i)=2pi/16. A = FINDGEN(30) * (!PI * 2 / 29.) ; define plotting symbol as filled circle USERSYM, radius * COS(A), radius * SIN(A), thick = 2., /fill ; draw plot plot, medsighist(valid), medstdhist(valid) ^ 2,$ psym = 8,background = white, color = black, xtitle = xtitle, $ ytitle = ytitle,$ yrange = [0, max(medstdhist(valid) ^ 2)], title = plotlabel, $ xtickv = indgen(11) * fix(upper / 10), xstyle = 1, xticks = 11,$ xminor = 10, xrange = [0, upper] ; add file name to bottom of plot xyouts, 1.0, 40.0, fsub, color = 0, charsize = filenamecharsize, /device ; Plot the curve and label oplot, [0, 1e6], rn2 + rgain * [0, 1e6], color = black, $ linestyle = linestyle ; Calculate gain and read noise for our plot theGain = strtrim(string(gain,format = '(f6.2)'), 2) + string(177B) $ + strtrim(string(error_gain, format = '(f5.2)'), 2) theReadNoise = strtrim(string(readnoise * gain,format = '(f6.1)'), 2) $ + string(177B) + strtrim(string(error_rn * gain, $ format = '(f5.1)'),2) ; Add SCA, date & time stamps, read noise and gain to plot label_x_posn = (!X.CRANGE[1] - !X.CRANGE[0]) * 0.95 + !X.CRANGE[0] label_y_posn = (!Y.CRANGE[1] - !Y.CRANGE[0]) * 0.05 + !Y.CRANGE[0] label_offset = 0.04 * (!Y.CRANGE[1] - !Y.CRANGE[0]) ylabels = ['Refsub Type: ' + refsubtype,$ 'Refsub Version: ' + refsubversion, $ 'SCA: ' + SCAName,$ 'Time: ' + theTime + ' UT',$ 'Date: ' + theDate, $ 'Read Noise: ' + theReadNoise + ' e- rms',$ 'Gain: ' + theGain+ ' e-/ADU'] for i = 0, n_elements(ylabels) - 1 do begin xyouts, label_x_posn, label_y_posn, ylabels[i], color = 0, $ align = 1.0, charsize = !P.CHARSIZE / 1.25 label_y_posn = label_y_posn + label_offset print, ' ' + ylabels[i] endfor ; write JPEG data to file if (plottype eq 'jpg') then begin ; make JPEG plot jpgimg = tvrd() write_jpeg, outdir + 'photonxfer_' + plotlabel + '.jpg', $ congrid(jpgimg, 1280, 1024, /interp, /center), quality = 100 endif ; finish up plot by closing device device, /close endfor ; set plot back to the default device set_plot, dname end