Pro dqe_reduce, testdir=testdir,temperature=temperature,region=region

    if not(keyword_set(testdir)) then testdir='\\rabbit\sipin\H1RG-022-SIPIN\cold1\dqe.3Jun05\'
    ;testdir='C:\figerdev\idtl\sipinrockwell\H1RG-021-SIPIN\cold1\dqe.9Aug04\'
    ;calfile='C:\figerdev\idtl\sipinrockwell\H1RG-021-SIPIN\cold1\dqe.9Aug04\si_sciencetech.txt'
    calfile='\\Rabbit\raid2\sipin\H1RG-021-SIPIN\cold1\dqe.9Aug04\si_sciencetech.txt'
    if not(keyword_set(temperature)) then temperature='160K'

    ; set plotting appearnace parameters
    ; set line and text characteristics for publication quality
    !P.CHARSIZE=8
    !P.CHARTHICK=10.
    !X.THICK=20.
    !Y.THICK=20.
    !P.THICK=20.
    symsize=10.
    filenamecharsize=8.
    !P.TICKLEN=0.03
    ; define forground and background colors
    white=255
    black=0

    ; specify region for analysis
    if not(keyword_set(region)) then begin
       x1=500
       x2=520
       y1=500
       y2=520
       x1=4
       x2=1019
       y1=4
       y2=1019
    endif else begin
       x1=region(0)
       x2=region(1)
       y1=region(2)
       y2=region(3)
    endelse
    ; set skip factor to allow for cases where pixel types are different in alternate rows
    skip=1

    ; add a path delimiter at the end of the path in case the user
    ; forgot to add one.
    if (strmid(testdir, strlen(testdir) - 1, 1) ne path_sep()) then testdir=testdir+path_sep()

    ; change directory to testdir
    cd,testdir

    ; readcol, testdir+infile , numbimgs, format='(A)'
    numbimgs=file_search('*on_'+temperature+'.fits')
    numbimgs_off=file_search('*off_'+temperature+'.fits')

    ; get exposure time
    dummy=readfits(numbimgs(0),h)
    itime=sxpar(h,'EXP_TIME')

    sig=fltarr(n_elements(numbimgs))
    sigoff=fltarr(n_elements(numbimgs))
    sigdiff=fltarr(n_elements(numbimgs))
    signorm=fltarr(n_elements(numbimgs))
    wav=fltarr(n_elements(numbimgs))
    sird=fltarr(n_elements(numbimgs))
    sibk=fltarr(n_elements(numbimgs))
    gerd=fltarr(n_elements(numbimgs))
    gebk=fltarr(n_elements(numbimgs))

    ; this loop cycles over the number of types of pixels on a device, i.e. 1 type for H1RG-18-SIPIN and two
    ; types for H1RG-021-SIPIN
    for j=0,skip-1 do begin

       for imgcount=0, n_elements(numbimgs)-1 do begin
         print,numbimgs[imgcount]
         catch, errval
         if (errval eq 0) then begin
             fits_read, numbimgs[imgcount], img, imghdr
             fits_read, numbimgs_off[imgcount],img_off,imghdr_off
         endif else begin
             continue
         endelse

         sig[imgcount]=median(img(x1:x2,j:y2:skip,1)-img(x1:x2,j:y2:skip,0))
         sigoff[imgcount]=median(img_off(x1:x2,j:y2:skip,1)-img_off(x1:x2,j:y2:skip,0))
         sigdiff[imgcount]=sig[imgcount]-sigoff[imgcount]
         signorm[imgcount]=sigdiff[imgcount]/(median(img(x1:x2,j:y2:skip,1)-img(x1:x2,j:y2:skip,0))-(median(img_off(x1:x2,j:y2:skip,1)-img_off(x1:x2,j:y2:skip,0))))  ; difference in region versus whole array
         print,'norm='+strtrim(string(signorm[imgcount]),2),'  difference='+strtrim(string(sigdiff[imgcount]),2)
         wav[imgcount]=sxpar(imghdr,'lambda')
         sird[imgcount]=sxpar(imghdr,'si_rd')
         gerd[imgcount]=sxpar(imghdr,'ge_rd')
         sibk[imgcount]=sxpar(imghdr,'si_bk')
         gebk[imgcount]=sxpar(imghdr,'ge_bk')

       endfor

       wav=wav[where(wav ne 0.0)]
       sigdiff=sigdiff[where(wav ne 0.0)]
       ; output device to frame buffer
       set_plot, 'z'
       device, set_resolution=[4800,3200]
       device, set_font='Courier'
       wavs=wav(sort(wav))
       sigdiffs=sigdiff(sort(wav))
       plot,wavs,sigdiffs,position=[0.15,0.15,0.8,0.95],xtitle='Wavelength (nm)',ystyle=8,$
         ytitle='Signal (ADU)',title=testdir,yrange=[0,max(sigdiff)],background=white,color=black
       ; set ytitle and axis depending on whether Germanium and Silicon readings are displayed
       ytitle='SI (Dash) Voltage'
       if (max(wav) gt 1300.) then ytitle='SI (Dash) and GE (Dot) Voltage'
       if (max(wav) gt 1300.) then axis,yaxis=1,yrange=[0,max([gerd-gebk,sird-sibk])],/save,ytitle=ytitle,color=black
       if (max(wav) lt 1300.) then axis,yaxis=1,yrange=[0,max(sird-sibk)],/save,ytitle=ytitle,color=black
       oplot,wavs,sird(sort(wav))-sibk(sort(wav)),color=black,linestyle=2
       ; suppress Germanium diode curve if wavelengths do not exceed 1300 nm
       if (max(wav) gt 1300.) then oplot,wavs,gerd(sort(wav))-gebk(sort(wav)),color=black,linestyle=1

       ; set text in JPG file name to indicate which pixels were processed
       if (skip eq 1) then ext=''
       if (skip eq 2 and j eq 0) then ext='_evens'
       if (skip eq 2 and j eq 1) then ext='_odds'

       jpgimg=tvrd()
       write_jpeg,'Results\dqedataplot'+'_'+temperature+ext+'.jpg', congrid(jpgimg, 1280, 1024, /interp, /center), quality=100

       ; read in calfile
       readcol,calfile,calwav,si1low_vpw,si1high_vpw,si2low_vpw,si2high_vpw       ; units in Volts/Watt
       h=6.6260755e-34                                 ; J*seconds
       c=3e8                                      ; m/s

       si1low_vpp=si1low_vpw*(h*c/(calwav*1e-9))                    ; Volts/(photon/s)
       si1high_vpp=si1high_vpw*(h*c/(calwav*1e-9))                    ; Volts/(photon/s)
       si2low_vpp=si2low_vpw*(h*c/(calwav*1e-9))                    ; Volts/(photon/s)
       si2high_vpp=si2high_vpw*(h*c/(calwav*1e-9))                    ; Volts/(photon/s)

       gain=1.0                                        ; e-/ADU
       pixelarea=0.018^2                  ; mm^2
       diodearea=!pi*(2.5/2)^2               ; mm^2
       qe=sigdiff/itime/pixelarea*gain/((sird-sibk)/interpol(si1high_vpp,calwav,wav)/diodearea)  ; (electrons/pixel area)/(photon/diode area)
       qenorm=qe/max(qe(where(finite(qe))))
       set_plot, 'z'
       device, set_resolution=[4800,3200]
       device, set_font='Courier'
       plot,wavs,qenorm(sort(wav)),position=[0.15,0.15,0.8,0.95],xtitle='Wavelength (nm)',ystyle=1,$
         ytitle='QE',title=testdir,yrange=[0,1.],background=white,color=black
       jpgimg=tvrd()
       write_jpeg,'Results\dqe'+'_'+temperature+ext+'.jpg', congrid(jpgimg, 1280, 1024, /interp, /center), quality=100

       openw,1,'Results\dqe'+'_'+temperature+ext+'.txt'
       qenormsort=qenorm(sort(wav))
       qesort=qe(sort(wav))
       for out=0,n_elements(wavs)-1 do begin
         printf,1,wavs(out),qesort(out),qenormsort(out),format='(F6.1,E16.3,F16.4)'
       endfor
       close,1

    endfor

end
