FUNCTION xrcf_hsi_read, file_name, FACILITY = facility, PBLURR = pblurr, $
	EXPOSURE = exposure, RAW_HSI=raw_hsi, $
	QE_DEPRESS=QE_DEPRESS, QE_INFO_STR=qe_info_str
;  Thanks to Harald here's how to get HSI images...
;  The FACILITY keyword causes the data to be converted to 
;  facility coordinates in mm w.r.t. the HSI center.
; 12/16/97 dd Added EXPOSURE keyword to return exposure time
; 4/12/99 dd Added RAW_HSI keyword to return the hsi structure
; 4/12/99 dd Added QE_DEPRESS keyword to add back in missing
;            events from QE depression region.
; 4/13/99 dd Added QE_INFO_STR to return summary string about
;            QE_DEPRESS effects.

; get info on the file
print, ''
fits_info, file_name
print, ''

; read it in
raw_hsi = mrdfits(file_name,1,hsi_evt_hdr)
; get the exposure time from the main header
dummy = mrdfits(file_name,0,hsi_hdr)
exposure = fxpar(hsi_hdr, 'EXPOSURE')
print, ''

; If QE_DEPRESS is specified, then do a statistical creation
; of events lost in the QE depressed region.
;
; Values of QE_DEPRESS are passed in and can be calculated by
; the calling routine by using hsi_qe_depress(runid) .
;
; Depression Parameters are:
;  - the 2D gaussian center
gcy = 2043.0  ; from B-K analysis
gcz = 2047.5  ; assumed
;  - the sigma of the gaussian(r) in pixels
gsig = 13.4   ; from B-K analysis (12.7 to 14.1)
;  - the QE depression factor = QE_at_depression / QE_elsewhere
qe_info_str = 'QE Depress: Not applied'
if n_elements(QE_DEPRESS) EQ 0 then QE_DEPRESS = 1.0
if QE_DEPRESS LT 0.99 then begin
  qe_info_str = 'QE Depress:'+STRING(QE_DEPRESS,FORMAT='(F5.2)')
  ; OK, let's do it...
  ; Find all the events near the depression:
  near_depres = WHERE(raw_hsi.Y GT 2000 AND raw_hsi.Z LT 2100 AND $
	raw_hsi.Z GT 2000 AND raw_hsi.Z LT 2100, n_near)
  if n_near GT 0 then begin
    print, ''
    print, ' QE_DEPRESS = '+STRING(QE_DEPRESS,FORMAT='(F6.3)')+' :'
    print, ' Candidates for correction: '+STRING(n_near)
    qe_info_str = qe_info_str + ', '+STRCOMPRESS(n_near)+ $
	' candidates'
    ; Candidate events to be added to event list...
    candids = raw_hsi(near_depres)
    ; Calculate the radius-squared from the depression center:
    cand_r2s = (candids.Y - gcy)^2 + (candids.Z - gcz)^2
    ; and the resulting relative QE:
    cand_relqe = 1.0 - (1.0 - QE_DEPRESS)*EXP(-1.0*cand_r2s/(2.0*gsig*gsig))
    ; The fraction of observed counts that should be added to
    ; compensate is given by:
    cand_extra = (1.0 - cand_relqe)/cand_relqe
    ; Generate uniform random variables for these
    cand_rand = RANDOMU(SEED,n_elements(near_depres))
    ; OK, if the random value is less than the "extra" fraction
    ; then this event will be added to the observed events:
    add_em = WHERE(cand_rand LT cand_extra, n_found)
    if n_found GT 0 then begin
      print, ' Candidates added: '+STRING(n_found)
      qe_info_str = qe_info_str + ', '+STRCOMPRESS(n_found)+ $
	' added, '+STRING(100.0*FLOAT(n_found)/FLOAT(n_near), $
	FORMAT='(F5.1)')+'%'
      raw_hsi = [raw_hsi, candids(add_em)]
    end

    print, ''
  end else qe_info_str = qe_info_str + ', No candidates'
end
qe_info_str = STRCOMPRESS(qe_info_str)

; rest of code uses "data"
data = raw_hsi

; Convert to facility coords?
if KEYWORD_SET(FACILITY) then begin
  x_title = 'Facility Y (microns)'
  y_title = 'Facility Z (microns)'
  ; Subtract center pixel offset
  data.Y = data.Y - 2048
  data.Z = data.Z - 2048
  ; Get HSI parameters from xrcf_params.pro
  xrcf_params, HSI_PIXEL = hsi_pixel, HSI_ANGLE = hsi_angle
  temp_Y = hsi_pixel * FLOAT(data.Y)   ; microns
  temp_Z = hsi_pixel * FLOAT(data.Z)   ; microns
  ; do the pixel blurring if desired
  if KEYWORD_SET(PBLURR) then begin
    temp_Y = temp_Y + hsi_pixel*(RANDOMU(SEED,n_elements(temp_Y))-0.5)
    temp_Z = temp_Z + hsi_pixel*(RANDOMU(SEED,n_elements(temp_Z))-0.5)
  end
  ; do the rotation
  rads = hsi_angle * !DTOR
  data.Y = FIX( temp_Y*COS(rads) - temp_Z*SIN(rads) )
  data.Z = FIX( temp_Y*SIN(rads) + temp_Z*COS(rads) )
  x_range = 1.E3*[-10.,10.]
  y_range = 1.E3*[-10.,10.]
end else begin
  x_title = 'HSI Y (pixel)'
  y_title = 'HSI Z (pixel)'
  x_range = 1.E3*[0.,4.]
  y_range = 1.E3*[0.,4.]
end

; plot the data
;;plot, data.Y, data.Z, PSYM=3, TITLE = file_name, $
;;	XTITLE = x_title, YTITLE = y_title, $
;;	XRANGE = x_range, YRANGE = y_range

RETURN, data
END