; NAME: ; crosstalk_analysis.pro (v1.5) ; ; PURPOSE: ; Reads a list of input CDS-subtracted images, groups the images by detector temperature and ; writes an output file showing adjacent-pixel crosstalk for cosmic ray hits for each temperature ; ; CALLING SEQUENCE: ; crosstalk_analysis, testdir[, infile[, outfile [,outdir [,skipline [,numread [,limitlo [, limithi]]]]] ; ; INPUTS: ; testdir the name of a directory in which the input data file resides. ; ; KEYWORD PARAMETERS: ; infile the name of the input data file. Defaults to 'crosstalkfiles.txt' ; outdir the name of the directory where output plots are written; defaults to testdir+'Crosstalk' ; outfile name of output file list; defaults to 'crosstalkfiles.txt' ; skipline number of initial lines to skip when reading input file ; numread number of lines to read from the input file ; limitlo Lower limit of cosmic ray intenstiy (in ADU) to consider--passed on to cosmic.pro (must be < limithi, default is 3000 ADU) ; limithi Upper limit of cosmic ray intenstiy (in ADU) to consider--passed on to cosmic.pro (must be > limitlo, default is 60000 ADU) ; ; ; EXAMPLE: ; crosstalk_analysis, '.' ; ; REFERENCE: ; ; ; MODIFICATION HISTORY: ; Written by: Ernie Morse, May 21, 2003 ; ; Ernie Morse, IDTL, May 23, 2003 (v1.1) ; Added functionality to read mask images from crosstalk_reduce and modified call to cosmic ; ; Ernie Morse, IDTL, May 30, 2003 (v1.2) ; Added generation of plots ; ; Ernie Morse, IDTL, June 3, 2003 (v1.3) ; Added more comments to improve readability ; Changed name of output mask images so as to avoid overwriting the initial masks ; Changed dimension of output plot to 1280 x 1024 ; ; Ernie Morse, IDTL, June 4, 2003 (v1.4) ; Initial mask now contains hit intensity values instead of just 1s and 0s, so changed to create a copy ; of the mask in the format that cosmic.pro expects (1s and 0s). ; ; Ernie Morse, IDTL, June 5, 2003 (v1.5) ; Added creation of cosmic ray histograms from initial mask files created by crosstalk_reduce.pro ; Added limitlo and limithi keyword to pass limit criteria to cosmic.pro ; ; Don Figer, IDTL, October 16, 2004 (v1.5) ; Replaced pathdelim variable with path_sep() ; Pro crosstalk_analysis, testdir, infile=infile, outfile=outfile, outdir=outdir, skipline=skipline, numread=numread, $ limitlo=limitlo, limithi=limithi ; 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 path_sep()) then testdir = testdir + path_sep() ; set outdir to testdir + 'Crosstalk' if not set by user if (not keyword_set(outdir)) then outdir = testdir ; 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 path_sep()) then outdir = outdir + path_sep() ; set outfile to a default value if it is not set by the calling context if (not keyword_set(outfile)) then outfile='crosstalkresults.txt' ; set infile to a default value if it is not set by the calling context if (not keyword_set(infile)) then infile='crosstalkfiles.txt' ; set default value for skipline param if (not keyword_set(skipline)) then skipline = 0 ; set default values for limit params if (not keyword_set(limitlo)) then limitlo = long(3000) if (not keyword_set(limithi)) then limithi = long(60000) ; stop if lower limit isn't less than upper limit if (limitlo ge limithi) then begin stop, 'Value of lower limit param (limitlo) must be be greater than upper limit (limithi)' endif ; read the files from the input file list readcol, testdir+infile, filenames, masknames, format='(A,A)', skipline=skipline ; set default value for numread parameter if (not keyword_set(numread)) then numread = n_elements(filenames) ; adjust size of file arrays if necessary if (numread lt n_elements(filenames)) then begin filenames=filenames[0:numread-1] masknames=maskname[0:numread-1] endif ; define labels array for output cosmic ray histograms extralabels = strarr(4) ; open the output file openw, outlun, outdir+outfile, /get_lun ; read detector temperatures from the header of each file name and store in array dettemps = fltarr(n_elements(filenames)) for imgcount=0,n_elements(filenames)-1 do begin fits_read, testdir+filenames[imgcount], 0, imghdr, /header_only dettemps[imgcount] = float(sxpar(imghdr, 'dettemp')) endfor ; round the detector temperatures to the nearest 0.1K dettemps = (floor((dettemps*10) + 0.5))/10.0 ; sort detector temperatures in increasing order and apply same sort order to file arrays sortorder = sort(dettemps) dettemps = dettemps[sortorder] filenames = filenames[sortorder] masknames = masknames[sortorder] ; make an array containing the index of the last occurence of each temperature uniquetemps = uniq(dettemps) ; declare some variables startindex = 0 results = fltarr(3,3, n_elements(uniquetemps)) ; will hold crosstalk results for each temperature validtemps = intarr(n_elements(uniquetemps)) ; marks valid temperatures (those for which at least two images are present) for tempcount=0,n_elements(uniquetemps)-1 do begin validtemps[tempcount] = 0 for imgcount=startindex,uniquetemps[tempcount]-1,2 do begin validtemps[tempcount] = 1 print, 'file #',strtrim(string(imgcount),2) + ' and ' + strtrim(string(imgcount+1),2) ; read two images at the same detector temperature fits_read, testdir+filenames[imgcount], im1 fits_read, testdir+filenames[imgcount+1], im2 ; make a difference image and normalize diff = im1 - im2 diff = diff - median(diff) ; read the accompanying mask file for image 1 fits_read, testdir+masknames[imgcount], initialmask1, maskhd1 ; cosmic expects a mask of 1s and 0s, while the mask just loaded ; has intensity values for hits and 0s elsewhere. Make a copy for input ; to cosmic finalmask1 = initialmask1 if (total(finalmask1) gt 0) then begin finalmask1[where(finalmask1 gt 0)] = 1 endif finalmask1 = fix(finalmask1) ; call cosmic to find cosmic ray events in the 1st image cosmic, diff, counts, cross, diag, mcount, finalmask1, limitlo, limithi ; store the events found in the 1st image if (imgcount eq startindex) then begin allcross = cross alldiag = diag endif else begin allcross = [allcross, cross] alldiag = [alldiag, diag] endelse ; read the accompanying mask file for image 2 fits_read, testdir+masknames[imgcount+1], initialmask2, maskhd2 ; cosmic expects a mask of 1s and 0s, while the mask just loaded ; has intensity values for hits and 0s elsewhere. Make a copy for input ; to cosmic finalmask2 = initialmask2 if (total(finalmask2) gt 0) then begin finalmask2[where(finalmask2 ne 0)] = 1 endif finalmask2 = fix(finalmask2) ; find cosmic ray events in the 2nd image diff = im2 - im1 diff = diff - median(diff) cosmic, diff, counts, cross, diag, mcount, finalmask2, limitlo, limithi ; store cosmic ray events along with those from the 1st image allcross = [allcross, cross] alldiag = [alldiag, diag] ; write the modified masks back out to disk mask1name = strmid(masknames[imgcount], 0, strpos(masknames[imgcount], '_initial.fits', /reverse_search)) + '_final.fits' mask2name = strmid(masknames[imgcount+1], 0, strpos(masknames[imgcount+1], '_initial.fits', /reverse_search)) + '_final.fits' fits_write, outdir+mask1name, finalmask1, maskhd1 fits_write, outdir+mask2name, finalmask2, maskhd2 ; now make histograms from the initial masks for the two images detector = sxpar(maskhd1, 'detname') exptime1 = sxpar(maskhd1, 'exp_time') exptime2 = sxpar(maskhd2, 'exp_time') ; IDTL FITS headers are non-standard, so the following function calls are necessary to extract ; time and date info from the header date and time keywords timeStamp = '' ; An empty string to hold the result dateStamp = '' ; Likewise gettime, maskhd1, timeStamp getdate, maskhd1, dateStamp jpgfile1 = outdir + strmid(filenames[imgcount], 0, strpos(filenames[imgcount], '_CDS16.fits', /reverse_search)) + '_cthist.jpg' jpgfile2 = outdir + strmid(filenames[imgcount+1], 0, strpos(filenames[imgcount+1], '_CDS16.fits', /reverse_search)) + '_cthist.jpg' psfile1 = outdir + strmid(filenames[imgcount], 0, strpos(filenames[imgcount], '_CDS16.fits', /reverse_search)) + '_cthist.ps' psfile2 = outdir + strmid(filenames[imgcount+1], 0, strpos(filenames[imgcount+1], '_CDS16.fits', /reverse_search)) + '_cthist.ps' extralabels[0] = 'Detector: ' + detector extralabels[1] = 'Date: ' + dateStamp + ', ' + timeStamp + ' UT' extralabels[2] = 'Exposure Time: ' + strtrim(string(exptime1),2) + ' seconds' ; find and reject cosmic ray hits that could be spuriously generated by crosstalk in 1st mask hitpix = where(initialmask1 gt 0) ctpix = hitpix + 4 inimage = initialmask1 if ((where(initialmask1[ctpix] gt 0))[0] ne -1) then begin ctpix = ctpix[where(initialmask1[ctpix] gt 0)] inimage[ctpix] = 0 endif inimage = inimage[where(inimage gt 0)] extralabels[3] = 'Number of Events: ' + strtrim(string(n_elements(where(inimage gt 0))), 2) if (n_elements(inimage) gt 1) then makehistogram_crosstalk, inimage, 'Cosmic Ray Intensity Histogram', 'Hit Intensity (ADU)', extralabels, psfile1, jpgfile1 ; change labels for second image ; IDTL FITS headers are non-standard, so the following function calls are necessary to extract ; time and date info from the header date and time keywords timeStamp = '' ; An empty string to hold the result dateStamp = '' ; Likewise gettime, maskhd2, timeStamp getdate, maskhd2, dateStamp extralabels[1] = 'Date: ' + dateStamp + ', ' + timeStamp + ' UT' extralabels[2] = 'Exposure Time: ' + strtrim(string(exptime2),2) + ' seconds' ; find and reject cosmic ray hits that could be spuriously generated by crosstalk in 2nd mask hitpix = where(initialmask2 gt 0) ctpix = hitpix + 4 inimage = initialmask2 if ((where(initialmask2[ctpix] gt 0))[0] ne -1) then begin ctpix = ctpix[where(initialmask2[ctpix] gt 0)] inimage[ctpix] = 0 endif inimage = inimage[where(inimage gt 0)] extralabels[3] = 'Number of Events: ' + strtrim(string(n_elements(where(inimage gt 0))), 2) if (n_elements(inimage) gt 1) then makehistogram_crosstalk, inimage, 'Cosmic Ray Intensity Histogram', 'Hit Intensity (ADU)', extralabels, psfile2, jpgfile2 endfor ; continue on to the next temperature if there wasn't at least two images at this temperature if (validtemps[tempcount] eq 0) then begin startindex = uniquetemps[tempcount]+1 continue endif ; define the cross talk as the 25% median of the values of adjacent and diagonal values for the cosmic ray events ; the rows of cross are defined as percentage crosstalk for: ; row 0: pixel to right of each event ; row 1: pixel to left of each event ; row 2: pixel above each event ; row 3: pixel below each event ; ; for alldiag: ; row 0: pixel above and right ; row 1: pixel below and left ; row 2: pixel above and left ; row 3: pixel below and right ; ; the next steps fills in the value of the result array with the 25% median of the proper row ; first each row of allcross and alldiag needs to be sorted for i=0,3 do begin allcross[*,i] = (allcross[*,i])[sort(allcross[*,i])] alldiag[*,i] = (alldiag[*,i])[sort(alldiag[*,i])] endfor ; convert values in allcross and alldiag to % allcross = allcross*100 alldiag = alldiag*100 ; determine the row index containing the 25% median value medindex = n_elements(allcross[*,0]) / 4 ; now fill in the values for the results array results[0,0,tempcount] = alldiag[medindex,2] ; above-left results[1,0,tempcount] = allcross[medindex,2] ; above results[2,0,tempcount] = alldiag[medindex,0] ; above-right results[0,1,tempcount] = allcross[medindex,1] ; left results[1,1,tempcount] = 100.0 ; central point (cosmic ray event) results[2,1,tempcount] = allcross[medindex,0] ; right results[0,2,tempcount] = alldiag[medindex,1] ; below-left results[1,2,tempcount] = allcross[medindex,3] ; below results[2,2,tempcount] = alldiag[medindex, 3] ; below-right ; print the detector temperature and crosstalk results to output file printf, outlun, 'Detector Temperature (K): ' + strtrim(string(dettemps[uniquetemps[tempcount]], format='(F5.1)'),2) printf, outlun, 'Number of images: ' + strtrim(string(imgcount-startindex),2) printf, outlun, 'Number of Events: ' + strtrim(string(n_elements(allcross[*,0])),2) printf, outlun, 'Crosstalk Results (%):' printf, outlun, results[*,*,tempcount] printf, outlun ; reset start index value for next trip through loop startindex = uniquetemps[tempcount]+1 endfor ; trim the results array to remove the frames of the cube that are empty (in other words, the ones ; for temps that were nonvalid because there wasn't at least two images at that temperature) validtempsindex=where(validtemps ne 0,validtempscount) if (validtempscount ne 0) then begin results = results[*,*,validtempsindex] ; turn the validtemps array into a list of the valid temperatures validtemps = dettemps[uniquetemps[where(validtemps ne 0)]] ; separate the results array into subarrays for each adjacent non-diagonal pixel left = results[0,1,*] right = results[2,1,*] above = results[1,0,*] below = results[1,2,*] ; find max and min vals for setting up plot limits maxval = max([max(left), max(right), max(above), max(below)]) minval = min([min(left), min(right), min(above), min(below)]) ;IDTL FITS headers are non-standard, so the following function calls are necessary to extract ;time and date info from the header date and time keywords timeStamp = '' ; An empty string to hold the result dateStamp = '' ; Likewise gettime, maskhd1, timeStamp getdate, maskhd1, dateStamp jpgfile = outdir + detector+'_ctplot.jpg' psfile = outdir + detector+'_ctplot.ps' ; set axis ranges for plots--the actual data range takes up the central 65% of the plot height yrng = (maxval - minval)/0.65 ylo = minval - yrng*0.2 yhi = maxval + yrng*0.15 xrng = (max(validtemps) - min(validtemps)) / 0.65 xlo = (min(validtemps)) - xrng*0.05 xhi = (min(validtemps)) + xrng*1.05 ; set plot colors black = 0 white = 255 for outcount = 0,1 do begin ;first set up for landscape Postscript plots if outcount eq 0 then begin set_plot, 'ps' !p.font = -1 device, filename=psfile, /landscape fsub = psfile ;filename to be written at bottom of plot !Y.MINOR=10 !y.margin=[6,2] !x.margin=[14,4] ; set line and text characteristics for publication quality !P.CHARSIZE=1.2 !P.CHARTHICK=4. !X.THICK=8. !Y.THICK=8. !P.THICK=8. !P.MULTI=0 !P.TICKLEN=0.03 symsize=1.5 filenamecharsize=1.1 endif else begin set_plot, 'z' device, set_resolution=[8000,6000] device, set_font='Courier' fsub = jpgfile !Y.MINOR=10 !y.margin=[6,4] !x.margin=[14,4] ; set line and text characteristics for publication quality !P.CHARSIZE=10 !P.CHARTHICK=15. !X.THICK=20. !Y.THICK=20. !P.THICK=20. symsize=10. filenamecharsize=8. endelse ; Plot the result plot, validtemps, left, $ TITLE='Adjacent Pixel Crosstalk', $ XTITLE='Detector Temperature (K)', $ YTITLE='Crosstalk (%)', $ psym=-2, $ XRANGE=[xlo,xhi], $ XSTYLE=1, $ YRANGE=[ylo,yhi], $ YSTYLE=1, $ BACKGROUND=white, $ COLOR=black, $ symsize=symsize,charsize=!P.CHARSIZE*1.5 oplot, validtemps, right, psym=-4, color=black, symsize=symsize oplot, validtemps, above, psym=-5, color=black, symsize=symsize oplot, validtemps, below, psym=-6, color=black, symsize=symsize oplot, [(!X.CRANGE(1)-!X.CRANGE(0))*0.95+!X.CRANGE(0)],[(!Y.CRANGE(1)-!Y.CRANGE(0))*0.07+!Y.CRANGE(0)], psym=2, color=black, symsize=symsize oplot, [(!X.CRANGE(1)-!X.CRANGE(0))*0.95+!X.CRANGE(0)],[(!Y.CRANGE(1)-!Y.CRANGE(0))*0.1+!Y.CRANGE(0)], psym=4, color=black, symsize=symsize oplot, [(!X.CRANGE(1)-!X.CRANGE(0))*0.95+!X.CRANGE(0)],[(!Y.CRANGE(1)-!Y.CRANGE(0))*0.13+!Y.CRANGE(0)], psym=5, color=black, symsize=symsize oplot, [(!X.CRANGE(1)-!X.CRANGE(0))*0.95+!X.CRANGE(0)],[(!Y.CRANGE(1)-!Y.CRANGE(0))*0.16+!Y.CRANGE(0)], psym=6, color=black, symsize=symsize xyouts, [(!X.CRANGE(1)-!X.CRANGE(0))*0.93+!X.CRANGE(0)],[(!Y.CRANGE(1)-!Y.CRANGE(0))*0.0675+!Y.CRANGE(0)], 'left = ', color = black, alignment=1 xyouts, [(!X.CRANGE(1)-!X.CRANGE(0))*0.93+!X.CRANGE(0)],[(!Y.CRANGE(1)-!Y.CRANGE(0))*0.0975+!Y.CRANGE(0)], 'right = ', color = black, alignment=1 xyouts, [(!X.CRANGE(1)-!X.CRANGE(0))*0.93+!X.CRANGE(0)],[(!Y.CRANGE(1)-!Y.CRANGE(0))*0.1275+!Y.CRANGE(0)], 'above = ', color = black, alignment=1 xyouts, [(!X.CRANGE(1)-!X.CRANGE(0))*0.93+!X.CRANGE(0)],[(!Y.CRANGE(1)-!Y.CRANGE(0))*0.1575+!Y.CRANGE(0)], 'below = ', color = black, alignment=1 ; write plot annotations xyouts, [(!X.CRANGE(1)-!X.CRANGE(0))*0.07+!X.CRANGE(0)],[(!Y.CRANGE(1)-!Y.CRANGE(0))*0.9+!Y.CRANGE(0)], 'Date: ' + dateStamp, color=black xyouts, [(!X.CRANGE(1)-!X.CRANGE(0))*0.07+!X.CRANGE(0)],[(!Y.CRANGE(1)-!Y.CRANGE(0))*0.94+!Y.CRANGE(0)], 'Detector: ' + detector, color=black ; write the file name outside the margin of the plot xyouts, 1.0, 40, fsub, charsize=filenamecharsize, color=black, /device ;For the JPEG plot, the plot image must be read in from memory and written out to a file if (outcount eq 1) then begin jpgimg = tvrd() write_jpeg, jpgfile, congrid(jpgimg, 1280, 1024, /center, /interp), quality=100 endif ;finish up by closing device device, /close endfor ; close output file free_lun, outlun endif End