; NAME: ; pro welldepth_reduce.pro (v1.3) ; ; PURPOSE: ; Produce plot of welldepth and linearity. ; ; CALLING SEQUENCE: ; welldepth_reduce, testdir, [outdir [, infile $ ; [, region [, scaname [, outfile]]]]] ; ; INPUTS: ; testdir string containing directory name wherein input images and a ; list of the image file names can be found ; ; KEYWORD PARAMETERS: ; outdir directory where output will be written; defaults to ; testdir + '\Results' ; infile input file name. File must be in directory specified by ; testdir. Default name is 'ramps.txt' ; region rectangular region of SCA to consider for analysis. Should ; be entered as 4-element 1-d integer array, i.e. ; [x1,x2,y1,y2] to consider pixels in rows x1 to x2 and ; columns y1 to y2. Default is specified as all non-reference ; pixels ; scaname Name of detector from which input images were obtained. This keyword ; only needs to be supplied if the input images were taken manually (i.e, ; without running the welldepth experiment procedure). ; outfile Name of output plot file. Default is scaname + "_linsat.jpg". ; ; EXAMPLE ; Calculate welldepth and linearity on images listed in file ramps_1.txt, using ; only pixels in rows 100 - 900 and columns 100 - 900. The default output directory ; and file name are used. ; ; IDL> wdir = '\\rabbit\raid2\H1RG-022-SIPIN\cold1\welldepth.14Jun05' ; IDL> welldepth_reduce, wdir, infile="ramps_1.txt", region=[100,900,100,900] ; ; REFERENCE: ; None ; ; MODIFICATION HISTORY: ; Written by: ; Ernie Morse, IDTL, June 16, 2005 ; Modification of original code by Eddie Bergeron (welldepth_reduce_1.0.pro), ; altered for increased generality ; Don figer, RIDL, August 15, 2007 ; Made extensive modifications to accommodate obtaining all data within a single ramp (at low gain) ; pro welldepth_reduce, testdir, infile = infile, outdir=outdir, region = region, rowinc = rowinc, scaname = scaname, outfile = outfile,e=e pathdelim = path_sep() ; save the current plot type plottype = !d.name ; check to see if input directory name ends in path delimiter--if not, ; add it to the end if (strmid(testdir, strlen(testdir)-1, 1) ne pathdelim) then begin testdir = testdir + pathdelim endif ; set outdir to testdir + 'Results' if not set by user if (keyword_set(outdir) eq 0) then begin outdir = testdir + 'Results' + pathdelim endif ; check to see if outdirectory name ends in path delimiter--if not, ; add it to the end if (strmid(outdir, strlen(outdir)-1, 1) ne pathdelim) then begin outdir = outdir + pathdelim endif file_mkdir,outdir if (not keyword_set(infile)) then begin infile = 'ramps.txt' endif if (not keyword_set(rowinc)) then begin rowinc = 1 endif ; read the list of ramps readcol, testdir + infile, ramps, format = '(A)' ; get information from the header of the 1st image ; images are assumed to be 3d image cubes fits_read, testdir + ramps[0], 0, h1, /header_only xsize = fix(sxpar(h1, 'naxis1')) ysize = fix(sxpar(h1, 'naxis2')) zsize = fix(sxpar(h1, 'naxis3')) if (not keyword_set(scaName)) then begin scaName = sxpar(h1, 'detname') endif dettemp = strtrim(string(sxpar(h1, 'dettemp'), format = '(F5.1)'), 2) wheel1 = sxpar(h1, 'wheel1') wheel2 = sxpar(h1, 'wheel2') if (not keyword_set(outfile)) then begin fname = outdir + scaName + '_linsat.jpg' endif else begin fname = outdir + outfile endelse ; find the number of pixels in one frame framepix = long(xsize) * long(ysize) ; get the exposure times for each read of the ramp exptime = (sxpar(h1, 'exp_time') / (zsize - 1)) * findgen(zsize) ; extract SCA type from SCA name scaType = strmid(scaName, 0, strpos(scaName, '-')) ; get voltages from the header if (strmid(scaType, 0, 1) eq 'S') then begin biases = ['VDETCOM', 'VDDUC', 'VDDCL', 'VROWOFF', 'VSSUC', $ 'ISLEWREF', 'VDDOUT', 'VSUB', 'VSSOUT1', 'VP', 'VNROW', $ 'VGGCL', 'VNCOL'] endif if (strmid(scaType, 0, 1) eq 'H') then begin biases = ['DSUB', 'VRESET', 'VBIASPWR', 'VBIASGTE', 'DRAIN', $ 'VPULLA7', 'VPULLA15', 'OUTWINA', 'VREF', 'VPULLA23', $ 'VPULLA31', 'VDDA', 'VDD'] endif if (strmid(scaType, 0, 4) eq 'H4RG') then begin biases = ['DSUB', 'VRESET', 'VBIASPWR', 'VBIASGAT', 'DRAIN', $ 'CELLDRAIN','VSUB','VPULLUP','VDDA','VDD'] endif ; if the user didn't set the region, define default regions based on the ; SCA type case scaType OF 'H1RG' : begin if (not keyword_set(region)) then region = [4, 1019, 4, 1019] lowrange=-31000. highrange=31000. end 'H2RG' : begin if (not keyword_set(region)) then region = [4, 2043, 4, 2043] lowrange=-31000. highrange=31000. end 'SB304': begin if (not keyword_set(region)) then region = [4, 2051, 0, 2047] lowrange=-31000. highrange=31000. end 'H4RG': begin if (not keyword_set(region)) then region = [4, 4091, 4, 4091] lowrange=0. highrange=65535. exptime=(sxpar(h1, 'exp_time')) * findgen(zsize) end ELSE : begin print, 'Detector type ' + scaType + $ ' not supported. Exiting procedure.' return end endcase ; find where images are off both rails. Be conservative ramp_medians = fltarr(n_elements(ramps), zsize) for rampcount = 0, n_elements(ramps) - 1 do begin for readcount = 0, zsize - 1 do begin ; read in last frame fits_read, testdir + ramps[rampcount], r, first = readcount * framepix, last = (readcount + 1) * framepix - 1 r = reform(r, xsize, ysize, /overwrite) ramp_medians[rampcount, readcount] = median(r[region[0]:region[1], region[2]:region[3]:rowinc]) endfor endfor if (n_elements(ramps) gt 1) then begin ; find offset between each of the ramps ramp_offsets = fltarr(n_elements(ramps) - 1) for i = 1, n_elements(ramps) - 1 do begin ; find nonrailed reads of 1st ramp nonrailed_1 = where(ramp_medians[i - 1, *] gt lowrange and ramp_medians[i - 1, *] lt highrange) ; find nonrailed reads of 2nd ramp nonrailed_2 = where(ramp_medians[i, *] gt lowrange and ramp_medians[i, *] lt highrange) ; use a linear fit to the last 3 good points of the first ramp and ; the first 3 points of the second ramp. p1 = poly_fit(exptime[nonrailed_1[n_elements(nonrailed_1) - 3:*]], ramp_medians[i - 1, nonrailed_1[n_elements(nonrailed_1) - 3:*]], 1, /double) p2 = poly_fit(exptime[nonrailed_2[0:2]], ramp_medians[i, nonrailed_2[0:2]], 1, /double) ; match to the midpoint between the ends of the two ramps midtime = mean([exptime[nonrailed_1[n_elements(nonrailed_1) - 1]], exptime[nonrailed_2[0]]]) ramp_offsets[i - 1] = poly(midtime, p1) - poly(midtime, p2) endfor ; find out where saturation occurs in the last ramp ; define saturation as the first read whose median signal is within 200 ADU ; of the previous read. ; start by calculating difference of median of each read in last ramp with ; the previous read. nonrailed = where(ramp_medians[n_elements(ramps) - 1, *] gt lowrange) med_diff = ramp_medians[n_elements(ramps) - 1, nonrailed[0] + 1:*] - ramp_medians[n_elements(ramps) - 1, nonrailed[0]:zsize - 2] ; presatread is the read # of the last read prior to saturation presatread = (where(med_diff lt 200.0))[0] + nonrailed[0] ; read in the last non-saturated read and the 1st saturated read fits_read, testdir + ramps[n_elements(ramps) - 1], topend, first = presatread * framepix, last = (presatread + 2) * framepix - 1 topend = double(topend) topend = reform(topend, xsize, ysize, 2, /overwrite) ; find saturated pixels in region of interest (roi) diff = abs(topend[*,*,1] - topend[*,*,0]) roi = diff[region[0]:region[1], region[2]:region[3]:rowinc] satpix = where(roi lt 200) ; get the baseline read of the 1st ramp, which is the first read ; that is not railed at the low end of the A/D range baseread = (where(ramp_medians[0,*] gt lowrange))[0] fits_read, testdir + ramps[0], bottomend, first = baseread * framepix, last = (baseread + 1) * framepix - 1 bottomend = double(bottomend) bottomend = reform(bottomend, xsize, ysize, /overwrite) ; calculate the offset in ADU between the 1st and last ramp and subtract ; it from the first ramp. total_offset = total(ramp_offsets) bottomend = bottomend - total_offset ; find accumulated signal between 1st and last ramps sig_accum = topend[*,*,1] - bottomend roi = sig_accum[region[0]:region[1], region[2]:region[3]:rowinc] ; make a histogram of the accumulated signal in the saturated pixels in the region binsize = 100 signal_min = min(roi) - (min(roi) mod binsize) - binsize / 2.0 signal_max = max(roi) - (max(roi) mod binsize) - binsize / 2.0 satpixhist = histogram(roi, binsize = binsize, locations = signal, min = signal_min, max = signal_max) ; convert from raw pixel counts to percentage of pixels in region satpixhist = satpixhist / float(n_elements(roi)) * 100.0 ; adjust signal bin values to show the center value of each bin signal = signal + binsize / 2.0 ; find cumulative percentages in each bin cumulative = fltarr(n_elements(signal)) cumulative[0] = satpixhist[0] for bincount = 1, n_elements(signal) - 1 do begin cumulative[bincount] = cumulative[bincount - 1] + satpixhist[bincount] endfor ; find the exposure times and corresponding values of non-railed medians nonrailed = where(ramp_medians gt lowrange and ramp_medians lt highrange) if (scatype eq 'H4RG') then nonrailed = where(ramp_medians gt 0. and ramp_medians lt highrange) ; convert to 2-D indices nonrailed_2d = array_indices(ramp_medians, nonrailed) ; get the exposure times for use in fitting fit_times = exptime[nonrailed_2d[1,*]] ; apply the offset to each of the ramp medians for rampcount = 0, n_elements(ramps) - 2 do begin ramp_medians[rampcount, *] = ramp_medians[rampcount, *] - total(ramp_offsets[rampcount:*]) endfor fit_signal = ramp_medians[nonrailed] fit_signal = fit_signal[1:*] - fit_signal[0] frac_rate = [fit_signal / fit_times[1:*]] ; get the median flux rate over the course of the experiment medflux = long(median(frac_rate)) fitvals = where(fit_signal lt max(fit_signal) * 0.8) p = poly_fit(fit_signal[fitvals], frac_rate[fitvals], 2, /double) p1 = p / p[0] ugood = where(fit_signal gt 0.0) ; construct string for display of region region_str = strtrim(string(region, format='(I4)'), 2) region_out = '[' + region_str[0] + ':' + region_str[1] + ', ' $ + region_str[2] + ':' + region_str[3] + ':' $ + strtrim(string(rowinc, format='(I3)'), 2) + ']' ; transform bias values into a string for plot output bias_out = strarr(n_elements(biases)) for i = 0, n_elements(biases) - 1 do begin bias_out[i] = strtrim(string(sxpar(h1, biases[i]), format = '(F7.3)'), 2) endfor endif else begin ; find nonrailed reads of 1st ramp nonrailed_1 = where(ramp_medians[0, *] gt lowrange and ramp_medians[0, *] lt highrange) ; find nonrailed reads of 2nd ramp nonrailed_2 = where(ramp_medians[0, *] gt lowrange and ramp_medians[0, *] lt highrange) ; find out where saturation occurs in the last ramp ; define saturation as the first read whose median signal is within 200 ADU ; of the previous read. ; start by calculating difference of median of each read in last ramp with ; the previous read. nonrailed = where(ramp_medians[0, *] gt lowrange) med_diff = ramp_medians[0, nonrailed[0] + 1:*] - ramp_medians[0, nonrailed[0]:zsize - 2] ; presatread is the read # of the last read prior to saturation presatread = (where(med_diff lt 200.0))[0] + nonrailed[0] ; read in the last non-saturated read and the 1st saturated read fits_read, testdir + ramps[n_elements(ramps) - 1], topend, first = presatread * framepix, last = (presatread + 2) * framepix - 1 topend = double(topend) topend = reform(topend, xsize, ysize, 2, /overwrite) ; find saturated pixels in region of interest (roi) diff = abs(topend[*,*,1] - topend[*,*,0]) roi = diff[region[0]:region[1], region[2]:region[3]:rowinc] satpix = where(roi lt 200) ; get the baseline read of the 1st ramp, which is the first read ; that is not railed at the low end of the A/D range baseread = (where(ramp_medians[0,*] gt lowrange))[0] fits_read, testdir + ramps[0], bottomend, first = baseread * framepix, last = (baseread + 1) * framepix - 1 bottomend = double(bottomend) bottomend = reform(bottomend, xsize, ysize, /overwrite) ; find accumulated signal between 1st and last ramps sig_accum = topend[*,*,1] - bottomend roi = sig_accum[region[0]:region[1], region[2]:region[3]:rowinc] ; make a histogram of the accumulated signal in the saturated pixels in the region binsize = 100 signal_min = min(roi) - (min(roi) mod binsize) - binsize / 2.0 signal_max = max(roi) - (max(roi) mod binsize) - binsize / 2.0 satpixhist = histogram(roi, binsize = binsize, locations = signal, min = signal_min, max = signal_max) ; convert from raw pixel counts to percentage of pixels in region satpixhist = satpixhist / float(n_elements(roi)) * 100.0 ; adjust signal bin values to show the center value of each bin signal = signal + binsize / 2.0 ; find cumulative percentages in each bin cumulative = fltarr(n_elements(signal)) cumulative[0] = satpixhist[0] for bincount = 1, n_elements(signal) - 1 do begin cumulative[bincount] = cumulative[bincount - 1] + satpixhist[bincount] endfor ; find the exposure times and corresponding values of non-railed medians nonrailed = where(ramp_medians gt lowrange and ramp_medians lt highrange) ; get the exposure times for use in fitting fit_times = exptime[nonrailed] fit_signal = ramp_medians[nonrailed] fit_signal = fit_signal[1:*] - fit_signal[0] frac_rate = [fit_signal / fit_times[1:*]] ; get the median flux rate over the course of the experiment medflux = long(median(frac_rate)) fitvals = where(fit_signal lt max(fit_signal) * 0.8) p = poly_fit(fit_signal[fitvals], frac_rate[fitvals], 2, /double) p1 = p / p[0] ugood = where(fit_signal gt 0.0) ; construct string for display of region region_str = strtrim(string(region, format='(I4)'), 2) region_out = '[' + region_str[0] + ':' + region_str[1] + ', ' + region_str[2] + ':' + region_str[3] + ':' + strtrim(string(rowinc, format='(I3)'), 2) + ']' ; transform bias values into a string for plot output bias_out = strarr(n_elements(biases)) for i = 0, n_elements(biases) - 1 do begin bias_out[i] = strtrim(string(sxpar(h1, biases[i]), format = '(F7.3)'), 2) endfor endelse black = 0 white = 255 ; Open the z-buffer device set_plot, 'z' ; set resolution and font device, set_resolution=[8000,6000] device, set_font='Courier' !P.MULTI = [0, 1, 2] ; Set up for two plots per page fsub = fname ; set text to print at bottom of plot ; set plot margin & tick marks !Y.MINOR=0 !y.margin=[6,4] !x.margin=[12,12] ; set some size values !p.charsize = 8.0 !p.charthick = 15.0 !p.thick = 10.0 !x.thick = 20.0 !y.thick = 20.0 plottitle = scaName + ' Well Depth Histogram at ' + dettemp + ' K ' ; define x and y coordinates of histogram peak ymax = max(satpixhist, ymax_sub) x_ymax = signal[ymax_sub] ; adjust signal minimum and maximum to enclose the central 95% of the ; cumulative curve signal_min = (signal[where(cumulative ge 2.5)])[0] signal_max = (signal[where(cumulative ge 97.5)])[0] ; set x-axis range so that peak is in center of plot delta_x = (x_ymax - signal_min) > (signal_max - x_ymax) xmax = x_ymax + delta_x xmin = x_ymax - delta_x xwidth = xmax - xmin xtitle = 'Well Depth (ADU)' if keyword_set(e) then begin signal=signal*e x_ymax=x_ymax*e delta_x=delta_x*e xmax=xmax*e xmin=xmin*e xwidth=xwidth*e xtitle='Well Depth (e-)' endif plot, signal, satpixhist, psym = 10, xtitle = xtitle, ytitle = '% Number of Pixels', title = plottitle, xrange = [xmin, xmax], yrange = [0, ymax * 1.1], ystyle = 9, xstyle = 1, background = white, color = black ; draw vertical line through histogram peak oplot, [x_ymax, x_ymax], [0, 100], color = black, thick=5 ; add annotation of peak value xyouts, x_ymax + xwidth * 0.012, ymax * 0.1, strcompress(string(x_ymax, format = '(I8)'), /rem), orient = 90, color = black ; add axis for cumulative percentage axis, yaxis = 1, yrange = [0, 100], /save, ytitle = 'Cumulative %', ytickv = indgen(11) * 10., yticks = 10, yminor = 5, color = black ; plot the cumulative percentage oplot, signal, cumulative, color = black, thick=5 ; add biases to plot for i = 0, n_elements(biases) - 1 do begin xyouts, x_ymax + xwidth * 0.3, 85. - (i*5), biases[i]+' = ' + bias_out[i] + ' V', color = black endfor xyouts, xmin + xwidth * 0.03, 94.0, 'Analysis Region = ' + region_out, color = black xtitle = 'Total DN in Well' medfluxtext=' ADU/second' if keyword_set(e) then begin fit_signal=fit_signal*e xtitle='Total e- in Well' p = poly_fit(fit_signal[fitvals], frac_rate[fitvals], 2, /double) p1 = p / p[0] medflux=medflux*3 medfluxtext=' e-/second' endif ; plot linearity plot, fit_signal[ugood], frac_rate[ugood]/p[0], yra = [0.6, 1.15], $ xra = [0, max(fit_signal)*1.1], xstyle = 1, $ ystyle = 1, title = scaName + ' Linearity', xtitle=xtitle, $ ytitle = 'Fractional Countrate', xmar = [8, 6], $ background = white, color = black ;openw,1,testdir+'out0.txt' ;printf,1,findgen(max(fit_signal) * 1.1),format='(I10)' ;printf,1, poly(findgen(max(fit_signal) * 1.1), p1),format='(F10.5)' ;printf,1,frac_rate[fitvals]/p[0],format='(F10.5)' ;close,1 oplot, fit_signal[ugood], frac_rate[ugood]/p[0], psym = 4, color = black, symsize = 5.0 circsym, 6, 1 oplot, fit_signal[fitvals], frac_rate[fitvals]/p[0], psym = 8, color = black oplot, findgen(max(fit_signal) * 1.1), poly(findgen(max(fit_signal) * 1.1), p1), linestyle = 3, color = black plots, [0, max(fit_signal) * 1.1], [1.0, 1.0], linestyle = 2, color = black a1 = strcompress(string(p1[1], format = '$(e20.10)'), /rem) a1 = strmid(a1, 0, 5)+strmid(a1, strpos(a1, 'e')) a2 = strcompress(string(p1(2), format = '$(e20.10)'), /rem) a2 = strmid(a2, 0, 5)+strmid(a2, strpos(a2, 'e')) xyouts, max(fit_signal) * 0.1, 0.8, 'fit = 1 + (' + a1 + ' * DN) + (' + a2 + ' * DN^2)', color = black xyouts, max(fit_signal) * 0.2, 0.77, ' = included in fit', /data, color = black xyouts, max(fit_signal) * 0.2, 0.74, ' = not included in fit', /data, color = black xyouts, max(fit_signal) * 0.1, 0.71, 'Analysis Region = ' + region_out, color = black xyouts, max(fit_signal) * 0.1, 0.68, 'Flux = ' + strtrim(string(medflux), 2) + medfluxtext, color = black circsym, 5, 1 oplot, [max(fit_signal) * 0.17], [0.775], psym = 8, color = black oplot, [max(fit_signal) * 0.17], [0.745], psym = 4, color = black, symsize = 5.0 ; write out jpg image jpgimg = tvrd() write_jpeg, fname, congrid(jpgimg, 1600, 1200, /center, /interp), quality=100 ; finish up by closing device device, /close !p.multi = 0 !p.font = (-1) end