; NAME: ; photon_analysis.pro (version 1.5) ; ; PURPOSE: ; This procedure reads a data file generated by photon.pro and uses the data ; to generate a photon transfer plot showing noise squared vs. signal. ; ; CALLING SEQUENCE: ; photon_analysis, outdir, snr_min [,datfile [,outps] [,outfile]] ; ; INPUTS: ; outdir - directory where data file is located ; snr_min - the minimum signal-to-noise ratio that allows points to be ; included in the fit that determines the slope (gain) and ; y-intercept (readnoise) ; ; KEYWORD PARAMETERS: ; datfile - name of data file generated by photon.pro. ; outps - flag to be set if PostScript and JPEG output are desired ; outfile - name of output plots (this should be the name of a Postscript ; file--a JPEG file with the same basefilename and a .jpg ; extension will also be created if the outps flag is set. ; Default is same basename as datfile with a .ps extension. ; ; EXAMPLE ; ; Create photon transfer plots for a data set ; IDL> photon_analysis, '.', 0.8, datfile='photon_data.txt', /outps, $ ; outfile='photon_plot.ps' ; ; REFERENCE: ; ; MODIFICATION HISTORY ; ; Written by: ; Ernie Morse, IDTL, March 4, 2003 ; ; Modified by: ; Ernie Morse, IDTL, March 26, 2003 (version 1.1) ; Changed to linear plotting instead of log-log ; Read line about bias subtraction from input file and print on plot ; Ernie Morse, IDTL, April 3, 2003 (v1.2) ; Changed to using z-buffer and write_JPEG to write jpgs instead of ; converting PS ; Ernie Morse, IDTL, April 24, 2003(v1.3) ; Changed to reject from fitting all signal points above 1e4 ADU ; Mike Regan, IDTL, May 8, 2003(v1.4) ; Plot points not included in fit with different symbol ; Ernie Morse, IDTL, February 5, 2004 (v1.5) ; Edited to insure compliance with IDTL coding standards ; Pro photon_analysis, outdir, snr_min, datfile = datfile, outps = outps, $ outfile = outfile ; set default values for keywords if (not keyword_set(datfile)) then datfile = 'photonxfer_data.txt' if (not keyword_set(outfile)) then outfile = strmid(datfile, 0, strpos(datfile, '.')) + '.ps' delim = path_sep() ; add the delim character to the end of outdir if it is not already there if (strlen(outdir) ne strpos(outdir, delim, /reverse_search) + 1) then outdir = outdir + delim ; define blank string variables to hold info from data file bias = (region = (scaname = (wheel1 = (wheel2 = (thetime = (thedate $ = (refcor_version = (refcor_type = '')))))))) ; open data file and read header info that will be printed on plot openr, lun1, outdir + datfile, /get_lun readf, lun1, scaname, format = '(A)' readf, lun1, wheel1, format = '(A)' readf, lun1, wheel2, format = '(A)' readf, lun1, refcor_version, format = '(A)' readf, lun1, refcor_type, format = '(A)' readf, lun1, thedate, format = '(A)' readf, lun1, thetime, format = '(A)' readf, lun1, region, format = '(A)' readf, lun1, bias, format = '(A)' free_lun, lun1 ; read tabular data values from the data file, skipping the header readcol, outdir + datfile, exp_group, avgadu, sigmavals, sigma_adu, $ sigma_sigma, format = '(D,D,D,D,D)', skip = 9 ; make two plots, one PS and one JPG ; define plot ranges and min-max values xmin = min(avgadu) xmax = max(avgadu) ymin = (min(sigmavals)) ^ 2 ymax = (max(sigmavals)) ^ 2 ; define x-range that is 10% wider than data range xrng = (xmax - xmin) / 0.9 xlo = xmin - xrng * 0.05 ; min of x-range xhi = xmax + xrng * 0.05 ; max of x-range ; define y-range that is about 40% taller than data range yrng = (ymax - ymin) / 0.7 ylo = ymin - yrng * 0.2 ; min of y-range yhi = ymax + yrng * 0.2 ; max of y-range ; construct file names for PS and JPEG images from outfile keyword basename = strmid(outfile, 0, strpos(outfile, '.')) psfile = outdir + basename + '.ps' jpgfile = outdir + basename + '.jpg' ; construct two plots (PS and JPEG) for plotcount = 1, 1 do begin ; Create some empty values for the empty plot x = [0, 0] y = [0, 0] ; Set up for plotting. There are two options, the X terminal and ; PostScript. if (not keyword_set(outps) AND plotcount eq 0) then begin ; set plotcount to 2 to skip JPG plotting if this is just going to ; the screen plotcount = 2 set_plot, 'win' ; Plot to x terminal !P.MULTI = 0 ; Plot fills device ; Setup color information for plotting loadct, 0 tvlct, r, g, b, /get nelements = n_elements(r) red = GetColor("red", nelements - 2) green = GetColor("green", nelements - 3) blue = GetColor("blue", nelements - 4) yellow = GetColor("yellow", nelements - 5) white = GetColor("white", nelements - 6) magenta = GetColor("magenta", nelements - 7) black = GetColor("black", nelements - 8) turquoise = GetColor("turquoise", nelements - 9) ; Set up to look good on X-terminal !p.charsize = 1.0 !x.charsize = 1.0 !y.charsize = 1.0 !p.charthick = 1.0 !x.thick = 1.0 !y.thick = 1.0 !p.thick = 1.0 !x.minor = 0 !y.minor = 0 asize = 1.0 ; charsize for annotations athick = 1.0 ; charthick for annotations lthick = 1.0 ; line thickness endif else begin ; plot to file if (plotcount eq 0) then begin ; set up postscript device, using vector fonts set_plot, 'ps' !p.font = -1 device, filename = psfile, /landscape ; set filename to be written at bottom of plot fsub = psfile ; set margin and tick values !Y.MINOR = 0 !y.margin = [6, 2] !x.margin = [14, 4] asize = 1.12 ; charsize for annotations athick = 4.0 ; charthick for annotations lthick = 8.0 ; line thickness ; define colors for plotting black = 0 white = 255 endif else begin ; set up to plot to memory buffer and then write to JPEG set_plot, 'z' device, set_colors = 256 device, set_resolution = [8000, 6000] device, set_font = 'Courier' ; set text to write at bottom of plot fsub = jpgfile ; set parameters for JPEG ; set tick marks and margins !Y.MINOR = 0 !y.margin = [6, 4] !x.margin = [20, 4] asize = 8.0 ; charsize for annotations athick = 15.0 ; charthick for annotations lthick = 20.0 ; line thickness ; plot colors black = 0 white = 255 red = 0 endelse endelse ; add /xlog and /ylog to the plot statement below to return to ; log-log plot plot, 10 ^ x, 10 ^ y, $ xtitle="Signal (ADU)", $ ytitle="Noise Squared (ADU)", $ title="Photon Transfer", $ /xstyle, $ /ystyle, $ xrange = [xlo, xhi], $ yrange = [ylo, yhi], $ xmargin= [20, 4], $ color = black, $ background = white, $ charsize = asize * 1.25, charthick = athick, xthick = lthick, $ ythick = lthick, /nodata linestyle = 0 ; print filename at bottom of plot if (keyword_set(outps)) then begin xyouts, 1.0, 40, fsub, charsize = asize, charthick = athick, $ color = black, /device endif ; define plot symbols ; The series will be represented by big filled circles ; Offset radius is about size of a character _radius = asize ; define arrays for x,y vertices, offset from point _X = fltarr(10) _Y = fltarr(10) ; define position angle and angle step size _theta = 0 _delta_theta = 2 * !PI / 10 for i= 0, 9 do begin & _X[i] = _radius * cos(_theta) & _Y[i] = _radius * sin(_theta) & _theta = _theta + _delta_theta & endfor ; Define the symbol usersym, _X, _Y, /FILL ; for linear plots fitxvals = (findgen(200) / 199) * xrng + xlo ; Rectify the points, negative signal is not allowed. avgadu = avgadu > 0 goodpoints = where(avgadu lt 1e4) ; Plot the observed points on the graph oplot, avgadu[goodpoints], sigmavals[goodpoints] ^ 2, psym = 8, $ color = red ; Plot the points that are greater than 10K as an asterisk since they ; were ignored in the fit. bigpoints = where(avgadu ge 1e4,bigcount) if (bigcount ne 0) then oplot, avgadu[bigpoints], sigmavals[bigpoints] ^ 2, psym = 1, color = black oplot, [!x.crange[0], !x.crange[1]], [0,0], linestyle = linestyle,thick = lthick, color = black ; Find the readnoise and gain using linear fit fitexy, avgadu[goodpoints], sigmavals[goodpoints] ^ 2, rn2, rgain, $ sigma_A_B, chisq, prob, $ X_SIG = sigma_adu[goodpoints], $ Y_SIG = 2. * sigma_sigma[goodpoints] readnoise = sqrt(rn2) gain = 1. / rgain error_rn = 0.5 * sigma_A_B[0] / readnoise error_gain = sigma_A_B[1] * gain ^ 2 ; Remove bad snr data and fit again if (snr_min gt 0.) then begin good_snr = goodpoints[where(abs(avgadu[goodpoints] / gain) $ ge rn2*snr_min, good_snr_count)] if (good_snr_count eq n_elements(goodpoints)) then begin bad_snr_count = 0 endif else begin bad_snr = goodpoints[where(abs(avgadu[goodpoints] / gain) $ lt rn2*snr_min, bad_snr_count)] endelse if (good_snr_count lt 3) then begin print,'' print,' Not enough data points with SNR larger than ', $ strtrim(string(snr_min),2), $ '. --- No SNR criteria implmented.' print,'' endif else begin if(bad_snr_count ne 0) then begin print,'' print,' Following data points were not considered for ' $ + 'fitting:' print,' Signal (ADU) = ', avgadu[bad_snr] print,'' endif fitexy, avgadu[good_snr], sigmavals[good_snr] ^ 2, rn2, rgain, $ sigma_A_B, chisq, prob,$ X_SIG = sigma_adu[good_snr], $ Y_SIG = 2 * sigma_sigma[good_snr] readnoise = sqrt(rn2) gain = 1. / rgain error_rn = 0.5 * sigma_A_B[0] / readnoise error_gain = sigma_A_B[1] * gain ^ 2 endelse endif ; Calculate the appropriate y-values for the fitted curve ; Change this to take the log base 10 to return to log-log plots fityvals = readnoise ^ 2 + (fitxvals / gain) ; Plot the curve and label oplot, [0, 1e6], rn2 + rgain * [0, 1e6], color = black, $ linestyle = linestyle, thick = lthick ; Calculate gain and read noise for our plot theGain = strtrim(string(gain, format = '(f6.2)'), 2) + $ string(177B) + $ strtrim(string(error_gain, format= '(f6.3)'), 2) theReadNoise = strtrim(string(readnoise * gain, format = '(f6.1)'), 2) $ + string(177B) + $ strtrim(string(error_rn * gain, format = '(f6.2)'), 2) ; Add SCA, date & time stamps, read noise and gain to plot _label_x_posn = xmin ;starting x-coordinate for labels _label_y_posn = yhi - yrng*0.05 ;starting y-coordinate for y labels _label_dy = yrng*0.03 ;vertical distance between labels ; define all the labels ylabel = strarr(11) ylabel[0] = 'SCA: '+SCAName ylabel[1] = 'Filters: ' + wheel1 + ', ' + wheel2 ylabel[2] = 'Read Noise: ' + theReadNoise + ' e- rms' ylabel[3] = 'Gain: ' + theGain + ' e-/ADU' ylabel[4] = 'Rebsub type: ' + refcor_type ylabel[5] = 'Refsub version: ' + refcor_version ylabel[6] = 'Analysis region: ' + region ylabel[7] = bias ylabel[8] = 'Start Time: ' + theTime + ' UT' ylabel[9] = 'Date: ' + theDate ylabel[10] = 'Reduced chi^2: ' $ + strtrim(string(chisq / (n_elements(avgadu) - 2.), $ format = '(f6.1)'), 2) ; print plot annotation at upper left (print label, increment position, ; repeat) for _i = 0, 3 do begin xyouts, _label_x_posn, _label_y_posn, ylabel[_i], color = black, $ charsize = asize, charthick = athick _label_y_posn = _label_y_posn - _label_dy print, ' ' + ylabel[_i] endfor ; set vertical position for labels at bottom of plot _label_y_posn = ylo + yrng * 0.05 ; print plot annotations at lower left (print label, ; increment position, repeat) for _i = 4, 6 do begin xyouts, _label_x_posn, _label_y_posn, ylabel[_i], color = black, $ charsize = asize, charthick = athick _label_y_posn = _label_y_posn + _label_dy print, ' '+ylabel[_i] endfor ; print annotations at lower right ; set an x-coordinate at right edge of plot _label_x_posn = xhi - xrng * 0.05 ; reset vertical position to lowest value _label_y_posn = ylo + yrng * 0.05 for _i=7,9 do begin xyouts, _label_x_posn, _label_y_posn, ylabel[_i], charsize = asize,$ charthick = athick, color = black, alignment = 1.0 _label_y_posn = _label_y_posn + _label_dy endfor print, ' ' + ylabel[10] print,' Probability: ', prob ; write the file name outside the margin of the plot xyouts, 1.0, 40, jpgfile, charsize = asize, charthick = athick, $ color = black, /device ; For the JPEG plot, the plot image must be read in from memory and ; written out to a file if (plotcount eq 1) then begin jpgimg = tvrd() write_jpeg, jpgfile, congrid(jpgimg, 1600, 1200, /center, /interp),$ quality = 100 endif ; finish up by closing device if (keyword_set(outps)) then begin device, /close endif endfor end