; NAME:
;			crosstalk_analysis_reag.pro (v1.2)
;
; PURPOSE:
;			Processes output from crosstalk_reduce_reag.pro to produce crosstalk tables and a summary text file.  The
;			input file consists of a series of lines containing the names of two files.  The first file is a data cube
;			containing frames produced from CDS frames (0th read is subtracted from 1st, median of frame is subtracted
;			from result).  The second file is a mask indicating the pixel positions of cosmic ray hits in the first
;			file's data cube.  An adjacent-pixel cosmic ray crosstalk plot is produced for each input file pair, and a
;			a text file summarzing the results for the entire input is also produced.
;
; CALLING SEQUENCE:
;			crosstalk_analysis_reag, 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 'reagmasks.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_reag, '.', limithi=1e5, limitlo=2e3
;
; REFERENCE:
;
;
; MODIFICATION HISTORY:
;       Written by:  Ernie Morse, August 15, 2003
;
;		Modified by:
;			Ernie Morse, IDTL, August 19, 2003 (v1.1)
;			Changed from single-page plots for each input data cube to 6-plot pages showing all results from
;			an input data set.
;
;			Ernie Morse, IDTL, August 21, 2003 (v1.2)
;			Added creation of cosmic ray intensity histograms
;
;
;


Pro crosstalk_analysis_reag, testdir, infile=infile, outfile=outfile, outdir=outdir, skipline=skipline, numread=numread, $
		limitlo=limitlo, limithi=limithi

	; set directory path delimiter based on operating system in use
	if (!version.os_family eq 'unix') then begin
		pathdelim = '/'
	endif else begin
		pathdelim = '\'
	endelse

	; 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 testdir = testdir + pathdelim

	; 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 pathdelim) then outdir = outdir + pathdelim

	; 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='reagmasks.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

	; declare some arrays for results
	results = fltarr(3,3,n_elements(filenames))
	exposures = intarr(n_elements(filenames))
	numhits = intarr(n_elements(filenames))
    exp_time = fltarr(n_elements(filenames))
    numsamrd = intarr(n_elements(filenames))
    readmode = strarr(n_elements(filenames))

	; set default value for numread parameter
	if (not keyword_set(numread)) then numread = n_elements(masknames)

	; adjust size of file arrays if necessary
	if (numread lt n_elements(masknames)) then begin
		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


	; declare some variables
	startindex = 0
	labels = strarr(4)



	for imgcount=0,n_elements(masknames)-1 do begin
			; read mask image
			fits_read, testdir+masknames[imgcount], mask

			; read processed, median-subtracted image cube
			fits_read, testdir+filenames[imgcount], cube, cubehd

			cubedim = size(cube)
			xsize = cubedim[1]
			ysize = cubedim[2]
			zsize = cubedim[3]

			for framecount=0,zsize-1 do begin
			    maskframe = mask[*,*,framecount]

				diff = cube[*,*,framecount]

			   ; call cosmic to select cosmic rays within the intensity range and return adjacent pixel
			   ; crosstalk values (ratio of intensity of pixel values to the value of cosmic ray hit)
			    cosmic, diff, counts, cross, diag, mcount, maskframe, limitlo, limithi

				; store the events found in the 1st image
				if (framecount eq 0) then begin
				    allcross = cross
					alldiag = diag
				endif else begin
					allcross = [allcross, cross]
					alldiag = [alldiag, diag]
			    endelse
            endfor

		; 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

		goodhits = where(finite(allcross[*,0]))
		allcross = allcross[goodhits, *]
		alldiag = alldiag[goodhits, *]
		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,imgcount] = alldiag[medindex,2]		; above-left
		results[1,0,imgcount] = allcross[medindex,2]		; above
		results[2,0,imgcount] = alldiag[medindex,0]		; above-right
		results[0,1,imgcount] = allcross[medindex,1]		; left
		results[1,1,imgcount] = 100.0					; central point (cosmic ray event)
		results[2,1,imgcount] = allcross[medindex,0]		; right
		results[0,2,imgcount] = alldiag[medindex,1]		; below-left
		results[1,2,imgcount] = allcross[medindex,3]		; below
		results[2,2,imgcount] = alldiag[medindex, 3]		; below-right

		; get parameters from image header
		detname = strtrim(sxpar(cubehd, 'detname'),2)
		exp_time[imgcount] = sxpar(cubehd, 'exp_time')
		numsamrd[imgcount] = sxpar(cubehd, 'numsamrd')
		readmode[imgcount] = strtrim(sxpar(cubehd, 'readmode'),2)
		energy = sxpar(cubehd, 'energy')
		angle = sxpar(cubehd, 'angle')

		exposures[imgcount] = zsize
		numhits[imgcount] = n_elements(allcross[*,0])

		; make array of cosmic ray hit intensities
		hitpix = where(mask gt 0)
		intensities = cube[hitpix]


		; make some labels for a histogram plot of cosmic ray intensities
		labels[0] = 'Detector: ' + detname
		labels[1] = 'Image cube: ' + filenames[imgcount]
		labels[2] = 'Number of Images: ' + strtrim(string(zsize),2)
		labels[3] = 'Number of Events: ' + strtrim(string(n_elements(hitpix)),2)


		; print the detector temperature and crosstalk results to output file
		printf, outlun, labels[1]
		printf, outlun, labels[2]
		printf, outlun, 'Number of Events: ' + strtrim(string(n_elements(allcross[*,0])),2)
		printf, outlun, 'Crosstalk Results (%):'
		printf, outlun, results[*,*,imgcount]
		printf, outlun

		; set filenames for histograms
		hist_ps = outdir + strmid(filenames[imgcount], 0, strpos(filenames[imgcount], '.fits')) + '_cthist.ps'
		hist_jpg = outdir + strmid(filenames[imgcount], 0, strpos(filenames[imgcount], '.fits')) + '_cthist.jpg'

		; make the cosmic ray intensity histograms
        makehistogram_crosstalk, intensities, 'Cosmic Ray Intensity Histogram', 'Hit Intensity (ADU)', labels, $
            hist_ps, hist_jpg


    endfor

	; make JPG and PS versions of the table

    ; set plot names
    psname = outdir + 'Crosstalk_multi' + '.ps'
    jpgname = outdir + 'Crosstalk_multi' + '.jpg'

    ; set plot colors
    black = 0
    white = 255

	!y.margin = [10,2]


	; make a PS plot followed by a JPEG plot
    for j=0,1 do begin
        if (j eq 0) then begin
        	!p.multi = [0, 3, 2]
	        ; set up for landscape Postscript plotting
		    set_plot, 'ps'
		    !p.font=-1
	        device, filename = psname, /landscape
	        fsub = psname               ; filename to be written at bottom of plot

			; set some size values
		    asize = 0.90		; charsize of annotations
		    athick = 4.0		; charthick of annotations
		    lthick = 8.0		; line thickness
    	endif else begin
          	; set up to plot to memory buffer.  This plot will later be written from memory to a JPEG file.
          	!p.multi = [0, 3, 2]
          	set_plot, 'z'
            device, set_resolution=[8000,6000]
            device, set_font='Courier'
            fsub = jpgname			  ;filename to be written at bottom of plot


			;set some size values
	        asize = 8.0		;charsize of annotations
	        athick = 15.0		;charthick of annotations
	        lthick = 20.0		;line thickness
    	endelse

		for pc=0, n_elements(filenames)-1 do begin

			; set up plot parameters
			plot, [0], [0], title='Crosstalk Percentages', background=white, xrange=[0,100], $
			    yrange=[0,100], xticks=1, yticks=1, charsize=asize*1.25, $
			    charthick=athick, color=black, xstyle=5, ystyle=5, /nodata

			; create the grid lines of the table
			oplot, [0.5, 0.5, 99.5, 99.5, 0.5], [16.0, 100.0, 100.0, 16.0, 16.0], color=black, thick=lthick
			oplot, [0.5,99.5], [44.0,44.0], color=black, thick=lthick
	        oplot, [0.5,99.5], [72.0, 72.0], color=black, thick=lthick
	        oplot, [33.5,33.5], [16.0,100.0], color=black, thick=lthick
	        oplot, [66.5,66.5], [16.0,100.0], color=black, thick=lthick

			; write crosstalk results into table
			for row = 0,2 do begin
			    for col = 0,2 do begin
			        xyouts, (col+1)*33, 84 - row*28, string(results[col,row,pc], format='(F6.2)'), $
			        color=black, charsize=asize*1.8, charthick=athick, alignment=1
			    endfor
			endfor

			; write annotations below the table
			xyouts, 0.5, 5.0, 'Read Mode: ' + readmode[pc] + ' ' + strtrim(string(numsamrd[pc]),2), color = black, $
			    charsize = asize*0.7, charthick=athick*0.7
			xyouts, 0.5, 10.0, 'Exposure Time: ' + strtrim(string(exp_time[pc], format='(F6.1)'),2) + ' seconds', color = black, $
				charsize = asize*0.7, charthick=athick*0.7
			xyouts, 0.5, 0.0, 'Image: ' + filenames[pc], color=black, charsize=asize*0.7, charthick=athick*0.7
			xyouts, 99.5, 10.0, 'Exposures: ' + strtrim(exposures[pc],2), color = black, $
				charsize = asize*0.7, charthick=athick*0.7, alignment=1
			xyouts, 99.5, 5.0, 'CR hits: ' + strtrim(string(numhits[pc]),2), color = black, $
				charsize = asize*0.7, charthick=athick*0.7, alignment=1

		endfor

		; write information that applies to all plots
		xyouts, 0.0, 0.075, 'Detector: ' + detname, color = black, charsize = asize, charthick=athick, /normal
		xyouts, 0.0, 0.05, 'Energy: ' + strtrim(string(energy),2) + ' MeV', charsize=asize, charthick=athick, $
		    color=black, /normal
		xyouts, 0.0, 0.025, 'Angle: ' + strtrim(string(angle),2) + ' degrees', charsize=asize, charthick=athick, $
		    color=black, /normal

    	; write name of file outside of plot margins at bottom of page
    	xyouts, 0.0, 0.0, fsub, charsize=asize, charthick=athick, color=black, /normal


    	; For the JPEG plot, the plot image must be read in from memory and written out to a file
		if (j eq 1) then begin
	  		jpgimg = tvrd()
	  		write_jpeg, jpgname, congrid(jpgimg, 1200, 960, /interp, /center), quality=100
		endif

       ; finish up plot by closing device
       ; device must be open for the above JPEG code to work, so close must happen here at the end
       device, /close
    endfor   ;j loop

	; close output file
	free_lun, outlun
End