;+
; NAME
;     FE_LINE
; PURPOSE
;     Reduce and fit RHESSI data in the 3-12 keV range to extract the
;       continuum and fit the Fe feature
; CALLING EXAMPLE
;     fe_line,t1,t2,/arch
; INPUT
;     Start and stop times, default 21-Apr-02 early phase
; OPTIONAL KEYWORD INPUT
;     tres (default 4.) time binning
;     seg_index_mask to choose front-segment detectors, default [1,3,4,6,8,9]
;     archive, writes a file to hardwired directory "fit_data"
; OUTPUT
;     fits, data, energies (former two are UTPLOT-friendly structures)
;     Writes a .genx file to a hardwired directory "fit_data"
; OPTIONAL OUTPUT
;       obj, the spectrum object
; DETAILS
;     The archived data file contains one data structure with tags
;       DAY             LONG          
;       TIME            LONG           
;       RATES           FLOAT     counts/sec/bin 
;       ERATES          FLOAT     statistical error in rates 
;       F4              FLOAT     average rate 3.0-5.67 keV 
;       FE              FLOAT     average rate 5.67-7.33 keV
;       FE_RES          FLOAT     average rate 5.67-7.33 keV (residual)
;       F8              FLOAT     average rate 7.33-10.00 keV 
;       F11             FLOAT     average rate 10.00-12.00 keV4
;       PARAMS          FLOAT     Array[4], the continuum fit
;       GAUSS           FLOAT     Array[6], the Fe feature Gaussian fit
;       EQW             FLOAT     Equivalent width from these fits 
;       ENERGIES        FLOAT     Array[45], 1/3 keV bins
;       LTIME           FLOAT     Mean live time of detectors selected 
; HISTORY
;     19-Mar-2003  HSH, written
;     24-jun-2003  revised for better user-friendliness (HSH)
;     25-jun-2003  added the spectrum object as a keyword output (HSH)
;     27-jun-2003  revised the outputs some more, including livetime (HSH)
;-


pro fe_line, t1, t2, fits, data, energies, tres = tres, $
  archive = archive, seg_index_mask = seg_index_mask, $
  obj = obj, qdebug=qdebug 
if n_elements(tres) eq 0 then tres=4.
if n_elements(eres) eq 0 then eres=1./3.
if n_elements(t1) eq 0 then t1 = '21-apr-02 00:43'
if n_elements(t2) eq 0 then t2 = '21-apr-02 00:44'
if n_elements(seg_index_mask) eq 0 then $
  seg_index_mask = [1,0,1,1,0,1,0,1,1,0,0,0,0,0,0,0,0,0]
obj = hsi_spectrum()
obj->set,obs_time_interval=[t1,t2]
obj->set,seg_index_mask=seg_index_mask
obj->set,sp_time_interval=tres
obj->set,sp_energy_binning = 3.+findgen(15/eres+1)*eres
obj->set,livetime_enable=1
obj->set,/sp_data_str
obj->set,sp_data_unit = 'rate'
ddd = obj->getdata()
livetime = obj->get(/livetime_ctr)
dets = where(seg_index_mask, n_dets)
ltime = total(livetime(*,dets),2)/n_dets
tt = anytim(obj->getaxis(/ut),/yoh)
nn = n_elements(tt)
data = {hsi_hr,time:0L, day:0L, rates:fltarr(45), erates:fltarr(45), $
  ltime:0.}
data = replicate(data,nn)
for i = 0, nn-1 do begin
  tint = anytim(tt(i),/int)
  data(i).time = tint.time
  data(i).day  = tint.day
  data(i).rates= reform(ddd(i).rate)
  data(i).erates= reform(ddd(i).erate)
  data(i).ltime = ltime(i) 
endfor
energies = obj->getaxis(/energy)       ; band centers
; 
;   save the reduced data
; 
;if keyword_set(archive) then $
;  savegen,data,energies,script,file='reduced/'+filena
; 
n_data = n_elements(data)
logddd = alog10(data.rates>1)
fits = {fit_morph,day:0L,time:0L,rates:fltarr(45),erates:fltarr(45),$
  f4:0.,fe:0.,fe_res:0.,f8:0.,f11:0.,params:fltarr(4),gauss:fltarr(6),eqw:0.,$
  energies:fltarr(45),ltime:0.}
fits = replicate(fits,n_data)
for i = 0, n_data-1 do begin
  ss_f4 = indgen(8)
  ss_fe = 8+indgen(6)
  ss_f8 = 14 + indgen(8)
  ss_f11 = 22 + indgen(6)
  ss = [ss_f4,ss_f8]
  params = poly_fit(energies(ss),logddd(ss,i),3,yfit)
  ss_all = [ss_f4,ss_fe,ss_f8] 
  yy = params(0)+params(1)*energies(ss_all)+params(2)*energies(ss_all)^2+$
    params(3)*energies(ss_all)^3
  resid = 10^logddd(ss_all,i)-10^yy
  gg = gaussfit(energies(ss_all),resid,aa)
  fits(i).time = data(i).time
  fits(i).day = data(i).day
  fits(i).f4 = total(data(i).rates(ss_f4))/8.
  fits(i).fe = total(data(i).rates(ss_fe))/6.
  fits(i).fe_res = total(resid(ss_fe))/6.
  fits(i).f8 = total(data(i).rates(ss_f8))/8.
  fits(i).f11 = total(data(i).rates(ss_f11))/6.
  fits(i).gauss = aa 
  fits(i).params = params
  fits(i).rates = data(i).rates
  fits(i).erates = data(i).erates
  fits(i).ltime = data(i).ltime
  fits(i).energies = energies
; fits(i).eqw = sqrt(2*!pi)*fits(i).gauss(2)*fits(i).gauss(0)/10^yy(11)
  ecent = aa(1)
  cont_flux = 10^(params(0) + params(1)*ecent + params(2)*ecent^2 + $
    params(3)*ecent^3)
  fits(i).eqw = sqrt(2*!pi)*fits(i).gauss(2)*fits(i).gauss(0)/cont_flux
endfor
if keyword_set(archive) then begin 
  filena = strmid(t1,0,9)+'_'+strmid(t1,10,5)
  savegen,file='fit_data/'+filena,fits,script
endif
if keyword_set(qdebug) then stop
end