PRO hsi_qe_corr
;
; Create QE correction files from the measured HSI counts
; vs azimuth data...
; The counts-per-shell per hrma-area look very similar (2%)
; for each shell - thus the azimuth-averaged HSI QE is about
; the same for photons from each mirror shell.  This means
; that the correction files (QEcorrection vs Azimuth) for
; each shell should each average to 1.0 .
; 1/14/99 dd
;

; Input parameters...
; ------------------------------------
; Al-K Data
;;runid = '106856'
;;meas_energy = 1.4867 ; keV
; Set the relative average QE for the different shells -
; hetgcal> areas = [ 308.749,204.456, 163.304, 94.0142]
; hetgcal> cnts = [409812.,276205.,221958.,123276.]
; hetgcal> print, cnts/areas
;       1327.33      1350.93      1359.17      1311.25
; hetgcal> print, (cnts/areas)/(0.25*TOTAL(cnts/areas))
;     0.992642      1.01029      1.01645     0.980615
; looks like they are the same to 2%...
;;qe_ave_values = [1.0,1.0,1.0,1.0]

; C-K Data
runid = '108944'
meas_energy = 0.277 ; keV
; Set the relative average QE for the different shells -
; looks like significant differences of 6, 0, 3, 4 %:
; hetgcal> areas = [327.491,213.282, 169.215,96.2923]
; hetgcal> cnts = [946796.,650011.,537270.,307910.]
; hetgcal> print, cnts/areas
;       2891.06      3047.66      3175.07      3197.66
; hetgcal> print, (cnts/areas)/(0.25*TOTAL(cnts/areas))
;     0.939307     0.990187      1.03158      1.03892
qe_ave_values = [0.94,1.0,1.03,1.04]
; ------------------------------------


; Lots of plots to see what's going on!
!p.multi = [0,6,4]

shells = ['1','3','4','6']

; Loop over the shells...

for is=0,3 do begin

  shell = shells(is)
  print, ''
  print, 'HSI Data from Shell '+shell+' :'

  ; Read in the measured data...
  hsidat = rdb_read(!DDXRCFCAL+'/cmp/hsi'+runid+'_az_hist_shell' + $
		shell + '.rdb')
  meas_dat = hsidat.(1)

  ; OK, to handle the data drop-outs due to HSI gaps and
  ; mirror support structure: iterate on this loop:
  ; first time using the full as-measured values, the next
  ; time(s) replace the measured values with the
  ; FFTfit for all values where the data is less than the
  ; fit!
  ; ---------------------
  for iter=0,2 do begin
  ; ---------------------

  if iter NE 0 then begin
    ; Replace measured values by FFTfit values for low data points:
    low_data = where(meas_dat LT ftdat)
    meas_dat(low_data) = ftdat(low_data)
  end

  ; Do an FFT to get a "fit" to SIN, COS functions
  trans = FFT(meas_dat)

  ; Logic to decide which coefs to keep...
  ; Keep any which are above some fraction of the sigma of amplitudes
  ; in the channels between 1/4 and 3/4 of full range
  amps = ABS(trans)
  ntrans = n_elements(trans)
  amps_hi_freq = amps( FIX(1.*ntrans/4.) : FIX(3.*ntrans/4.) )
  ave_hi_freq = TOTAL(amps_hi_freq)/FLOAT(n_elements(amps_hi_freq))
  sigma_hi_freq = SQRT( TOTAL((amps_hi_freq-ave_hi_freq)^2) / $
			FLOAT(n_elements(amps_hi_freq)) )
  print, 'High Freq Ave = '+STRING(ave_hi_freq)
  print, 'High Freq RMS = '+STRING(sigma_hi_freq)

  ; Now decide which frequencies to keep...
  ; Keep all with amplitudes GE ave + fraction x rms 
  ; and which are low frequency!
  sigfrac = 1.0

  ; Don't forget the low freq ones at the end of the array!
  ; (THe highest freq is in the middle - the Nyquist freq.)
  ny_freq = n_elements(trans)/2.0
  bin_freq = 1.0-ABS((indgen(n_elements(trans))-ny_freq)/ny_freq)

  keepers = where(amps GE (ave_hi_freq + sigfrac*sigma_hi_freq) AND $
		bin_freq LT 0.05)

  ; Now do the inverse transform on only these...
  select_trans = trans*0.0
  for ik=0,n_elements(keepers)-1 do select_trans(keepers(ik)) = $
						trans(keepers(ik))
  ftdat = ABS(FFT(select_trans, /INVERSE)) ; keep it real


  ; Make some plots to show what is going on...
  ; plot the data and the FFT-derived fit
  plot, hsidat.(0), meas_dat, TITLE='HSI Data from Shell '+shell, $
	XTITLE='Azimuth (degrees from +Y to +Z)', YTITLE='counts', $
	PSYM=10
  oplot, hsidat.(0), ftdat, THICK=2

  ; plot the amplitudes of the transform
  plot_io, amps, PSYM=4, XRANGE=ntrans*[-0.1,1.1], XSTY=1, $
		XTITLE='Transform Frequency', YTITLE='Amplitude', $
		TITLE='FFT of HSI QE curve'
  ; show the average and rms levels...
  oplot, ntrans*[0.,1.], ave_hi_freq*[1.,1.], LINESTYLE=2
  oplot, ntrans*[0.,1.], ave_hi_freq*[1.,1.]+sigma_hi_freq, LINESTYLE=1
  oplot, ntrans*[0.,1.], ave_hi_freq*[1.,1.]-sigma_hi_freq, LINESTYLE=1
  ; Mark the ones that are kept:
  oplot, keepers, amps(keepers), PSYM=2  

  ; and do it again!
  ; ---------------------
  end ; of iteration
  ; ---------------------

  ; Now done with iterating so we have a nice smooth, well-fitting
  ; QE vs Azimuth curve - save it to an rdb file

  ; Make it have a given average QE correction value
  norm_fit = qe_ave_values(is) * ftdat/(TOTAL(ftdat)/FLOAT(n_elements(ftdat)))

  ; Write it out...
  qe_file_name = 'hsi_qecorr_'+STRING(meas_energy,FORMAT='(F4.2)')+ $
		'mp'+shell+'.rdb'
  qe_hdr = ['# ', $
		'# ', $
		'# ']
  qe_struct = {angle:0.0, qecorr:0.0}
  qe_dat_out = REPLICATE(qe_struct, n_elements(ftdat))
  qe_dat_out.angle = hsidat.(0)
  qe_dat_out.qecorr = norm_fit
  ; Put it in the !DDXRCFCAL/cmp/ directory...
  rdb_write, !DDXRCFCAL+'/cmp/'+qe_file_name, qe_dat_out, HEADER=qe_hdr

end ; of this shell...


RETURN
END