; 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 'dark.lst'
;			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, 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)

			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)

			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)
	results = results[*,*,where(validtemps ne 0)]

	; 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

	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, [xhi - xrng*0.035], [ylo + (yrng*0.032)], psym=2, color=black, symsize=symsize
	oplot, [xhi - xrng*0.035], [ylo + (yrng*0.052)], psym=4, color=black, symsize=symsize
	oplot, [xhi - xrng*0.035], [ylo + (yrng*0.072)], psym=5, color=black, symsize=symsize
	oplot, [xhi - xrng*0.035], [ylo + (yrng*0.092)], psym=6, color=black, symsize=symsize

	xyouts, xhi - xrng*0.05, ylo + (yrng*0.03), 'left = ', color = black, alignment=1
	xyouts, xhi - xrng*0.05, ylo + (yrng*0.05), 'right = ', color = black, alignment=1
	xyouts, xhi - xrng*0.05, ylo + (yrng*0.07), 'above = ', color = black, alignment=1
	xyouts, xhi - xrng*0.05, ylo + (yrng*0.09), 'below = ', color = black, alignment=1

	; write plot annotations
	xyouts, xlo+xrng*0.025, yhi - (yrng*0.07), 'Date: ' + dateStamp, color=black
	xyouts, xlo+xrng*0.025, yhi - (yrng*0.05), '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
End