;+
; Modifications:
; 9-Sep-2004, Kim.  Modified edata calc - previously summed error in counts and sigma in
;   quadrature.  Wrong.  Now just use sigma as error.
; 16-Sep-2004, Kim.  Added show_err option in plot_setup method.
; 10-Oct-2004, Kim.  In process, check if data is structure before using it.
; 09-Aug-2005, Kim.  In get, found needs to be param name, not 1
; 13-Feb-2006, Kim.  In plot_setup, call get_plot_title to get title
;---------------------------------------------------------------------------
;-

;---------------------------------------------------------------------------

function  spex_bk::init, source = source, _extra = _extra

if not keyword_set( source ) then source = obj_new( 'spex_bkint', _extra=_extra )

return, self->framework::init( source = source, $
                                            control = spex_bk_control(), $
                                            info = {spex_bk_info}, $
                                            _extra = _extra )

end

;--------------------------------------------------------------------------
pro spex_bk::set, $
	this_band=this_band, this_order=this_order, $
	done=done, not_found=not_found, $
	_extra=_extra

if exist(this_band) and exist(this_order) then begin
	bk_sep = self->get(/spex_bk_sep)
	if bk_sep then begin
		bk_eband = self->get(/spex_bk_eband)
		nband = bk_eband[0] eq -1 ? 0 : n_elements(bk_eband)/2
	endif else nband = 1
	if this_band gt nband-1 then begin
		message,/cont, 'Invalid value for this_band = '+trim(this_band)+$
			(bk_sep ? '.  Only ' + trim(nband) + ' bk ebands defined.  Aborting.' : $
				'.  Separate bk bands not enabled.')
		done=1
		return
	endif
	bk_order = self -> get(/spex_bk_order)
	norder = n_elements(bk_order)
	while n_elements(bk_order) lt nband do bk_order = append_arr(bk_order, bk_order[norder-1])
	bk_order[this_band] = this_order
endif

self->framework::set,spex_bk_order=bk_order, this_band=this_band, $
	_extra=_extra, done=done, not_found=not_found

end

;--------------------------------------------------------------------------

function spex_bk::get, this_band=this_band, this_order=this_order, $
	_extra=_extra, found=found, not_found=not_found

if exist(this_band) and keyword_set(this_order) then begin

	bk_order = self -> get(/spex_bk_order)
	bk_sep = self -> get(/spex_bk_sep)

	found=['this_band','this_order']

	bk_eband = self -> get(/spex_bk_eband)
	neband = bk_sep ? n_elements(bk_eband)/2 : 1
	if this_band gt neband-1 then ret = -1 else begin
		norder = n_elements(bk_order)
		while n_elements(bk_order) lt neband do bk_order = append_arr(bk_order, bk_order[norder-1])
		ret = bk_order[this_band]
	endelse

	return, ret
endif

ret = self -> framework::get(this_band=this_band, _extra=_extra, found=found, not_found=not_found)
return, ret
end

;--------------------------------------------------------------------------

; data product created is data,edata,ltime structure for expected background counts
; for every raw time and energy bin.

; explictly name spex_units keyword to keep it out of _extra
pro spex_bk::process, spex_units=spex_units, _extra = _extra

;print,'in spex_bk::process'

; get data in original units, then convert to rate for bk  calculation
data_obj = self -> get(/obj, class='spex_data')
data = data_obj -> getdata(/original_units, _extra=_extra)
if not is_struct(data) then message,'No data.'
units_str = data_obj -> getunits(/orig)
data = data_obj -> convert_units (data, 'rate', units_str )
times = data_obj -> get(/spex_ut_edges)
energies = data_obj -> get(/spex_ct_edges)

;; get data in bk time intervals, mainly so we can get indices of data in bk intervals
;; bk_index will be start, end index for each interval
;bkint_obj = self -> get(/obj, class='spex_bkint')
;bk = bkint_obj -> getdata(/original_units, _extra=_extra)
;if bk.data[0] eq -1 then begin
;	self -> setdata, {data: -1, edata: -1, ltime: -1}
;	self -> set, spex_bk_origunits=units_str
;endif else begin
;	units_str = bkint_obj->getunits(/orig)
;	bk = data_obj -> convert_units (bk, 'rate', units_str )
;	bk_index = bkint_obj -> get(/spex_bk_index)

; If bk_sep then for each bk_eband, save bk_order (order of background fit) and selected (pointer
; array of indices into data to use for background.
; If not bk_sep, then use 0th bk_time and bk_order
bk_sep = self -> get(/spex_bk_sep)
bk_eband = self -> get(/spex_bk_eband)
neband = bk_sep ? n_elements(bk_eband)/2 : 1


; BRIAN - set bk_ratio to 0 to calculate background the old way.  Set bk_ratio to 1 to use your new way.
; After setting, be sure to recompile.  You need to change something about the background after recompiling
; or it won't know it needs to recompute.  Or you can set the need update flag for this object by typing:
;  (o->get(/obj,class='spex_bk')) -> set, /need_update

bk_ratio = 1

bk_ratio = bk_ratio and bk_sep

bk_index = self -> getdata(class='spex_bkint')
have_bk = self -> get(/spex_have_bk)
if have_bk then begin

	bkcounts=data.data * 0.
	ebkcounts = data.edata * 0.

	; for each energy bin, fit selected points in data.data with polynomial of order 'order' to each
	; time bin in times.
	sel = *bk_index[0]
	bk_order = intarr(neband)
	for ib=0,neband-1 do bk_order[ib] = self -> get(this_band=ib, /this_order)
	order = bk_order[0]

	nen = n_elements(data.data[*,0])
	for i = 0,nen-1 do begin
		if bk_sep then begin
			q = where (energies[0,i] ge bk_eband[0,*] and energies[1,i] le bk_eband[1,*], count)
			q = count gt 0 ? q[count-1] : 0
			sel = *bk_index[q]
			order = bk_order[q]
		endif
		if sel[0] eq -1 then begin
			bkcounts[i,*] = 0.
			ebkcounts[i,*] = 0.
		endif else begin
			if bk_ratio and (q eq neband-1) then begin
				tmid = get_edges(times, /mean)
				base_time = min(tmid)
				bkcounts[i,*] = interpol (data.data[i,sel], tmid[sel]-base_time, tmid-base_time)
				ebkcounts[i,*] = data.edata[i,*]
			endif else begin
				bkcounts[i,*] = fit_backgrnd (times - min(times), $  ; wants relative times
					data.data[i,*], $
					data.edata[i,*], $
					order, $
					selected=sel, $
					sigma=sigma, $
					ltime=data.ltime[i,*])
				ebkcounts[i,*] = sigma
			endelse
		endelse
	endfor

	; computed bk as rate
	bk = {data: bkcounts, edata: ebkcounts, ltime: data.ltime}

	; If bk_ratio is one, base the time profile of the background in each energy bin on time profile
	; of highest bk eband.
	; For each energy bin, find bk_eband it's in, and get ratio of bk in that bin to the high energy profile,
	; both averaged over the bk_time_interval for that bin (the indices for times selected are already
	; in bk_index pointer).  Then multiply the high-energy profile by the ratio.
	; Note that we're doing this in rate space.  data is the raw data structure in rate, bk is the
	; background structure we just calculated in rate.
	; !!! need to recompute bk.edata too, but not sure how right now.
	if bk_sep and bk_ratio then begin
		high_band = bk_eband[*,neband-1]
		profile = reform(self -> bin_data (data=bk, intervals=high_band, er=z_err))
		width = n_elements(where(tmid lt base_time+60.))
		profile = smooth (profile, width, /edge_truncate)
		for i = 0,nen-1 do begin
			q = where (energies[0,i] ge bk_eband[0,*] and energies[1,i] le bk_eband[1,*], count)
			q = count gt 0 ? q[count-1] : 0
			sel = *bk_index[q]
			ratio = average(data.data[i,sel]) / average(profile[sel])
			print,'i, e[i], q, ratio, minmax(sel) ', i, energies[0,i], q, ratio, minmax(sel)
			;if i eq 114 then stop
			bk.data[i,*] = ratio * profile
		endfor
	endif


	; Store bk as counts
	bk = self -> convert_units (bk, 'counts', units_str)
	self -> set, spex_bk_origunits=units_str

	self->setdata, bk
endif else 	self -> setdata, {data: -1, edata: -1, ltime: -1}

end

;--------------------------------------------------------------------------

function spex_bk::getunits, orig=orig
if keyword_set(orig) then return, self->get(/spex_bk_origunits)

return, self->get(/spex_bk_units)
end

;--------------------------------------------------------------------------

pro spex_bk::setunits, units, orig=orig

if keyword_set(orig) then self -> set, spex_bk_origunits=units else $
	self -> set, spex_bk_units=units

end

;--------------------------------------------------------------------------
; Plot method plots time profile of background computed by fitting over
; background intervals defined.
; If intervals=-1 then uses raw energy bands.  If no intervals passed in, then
; uses spex_eband to bin in energy.

function spex_bk::plot_setup, $
	intervals=intervals, $
	pl_time=pl_time, $
	pl_energy=pl_energy, $
	photons=photons, $
	drm_eff=drm_eff, $
	show_err = show_err, $
	plot_params=plot_params, $
	_extra=_extra

error_catch = 0
if spex_get_debug() eq 0 then catch, error_catch
if error_catch ne 0 then begin
	catch, /cancel
	message, !error_state.msg, /cont
	;message,'Error making plot.', /cont
	return, 0
endif

photons = keyword_set(photons)
if photons and not keyword_set(drm_eff) then message,'No efficiency factors provided.  Cannot display photons.'

pl_time = keyword_set(pl_time)
pl_energy = keyword_set(pl_energy)
if not (pl_time or pl_energy) then pl_time = 1

data = self -> getdata(_extra=_extra)
if not is_struct(data) then message, 'No data accumulated.'
if data.data[0] eq -1 then message, 'No background intervals set.'

label = 'Detectors: ' + self -> get(/spex_detectors)

title = self -> get_plot_title(photons=photons)

units_str = self -> getunits()

if pl_energy then begin

	energies = self -> getaxis(/ct_energy, /edges_2)
	checkvar, intervals, -1  ; if no time bins passed in, don't bin
	z = self -> bin_data (data=data, intervals=intervals, er=z_err, /do_time, newedges=newedges)
	if z[0] eq -1 then return,0

	xdata = energies
	ydata = z
	plot_type='xyplot'
	title = title + ' vs Energy'
	xtitle = 'Energy (keV)'
	dim1_id = format_intervals(newedges, /ut) + ' (Bk)'
	dim1_vals = anytim(newedges, /vms)
	dim1_sum = 1
	dim1_unit = ''
	xlog = 1
	ylog = 1

endif else begin
	checkvar, intervals, self->get(/spex_eband)  ; if no energy bins passed in, use eband

	z = self -> bin_data (data=data, intervals=intervals, er=z_err, newedges=newedges)
	if z[0] eq -1 then return, 0

	xdata = anytim(self -> getaxis(/ut))
	utbase = min(xdata)
	xdata = xdata - utbase
	ydata = transpose(z)
	plot_type = 'utplot'
	xtitle = ''
	dim1_id = format_intervals(newedges, format='(f9.1)') + ' keV (Bk)'
	dim1_vals = newedges
	dim1_sum = 0
	dim1_unit = units_str.ct_edges
	xlog = 0
	ylog = 1
endelse

if photons then spex_apply_eff, ydata, drm_eff, units_str=units_str

checkvar, utbase, 0
plot_params = { $
	plot_type: plot_type, $
	xdata: xdata, $
	utbase: anytim(utbase,/vms), $
	ydata: ydata, $
	id: title, $
	label: label, $
	data_unit: units_str.data, $
	dim1_id: dim1_id, $
	dim1_vals: dim1_vals, $
	dim1_sum: dim1_sum, $
	dim1_unit: dim1_unit, $
	dim1_enab_sum: 1, $
	weighted_sum: units_str.data_type ne 'counts', $
	xlog: xlog, $
	ylog: ylog, $
	xtitle: '' $
	}

if keyword_set(show_err) then plot_params = add_tag (plot_params, z_err, 'edata')
return, 1
end

;-------------------------------------------------------------------------------

;pro spex_bk::plot, _extra=_extra
;
;status = self -> plot_setup (plot_params=plot_params, _extra=_extra)
;
;if status then self -> do_plot, plot_params=plot_params, _extra=_extra
;
;end

;---------------------------------------------------------------------------

pro spex_bk__define

dummy = {spex_bk, $
         inherits spex_gen}

end