Viewing contents of file '../idllib/ghrs/pro/calfos_bac.pro'
pro calfos_bac,bac,fwidths,config,pattern,h,data,eps,back,rootname,ctable,rtable
;+
;			calfos_bac
;
; CALLING SEQUENCE:
; calfos_bac,bac,fwidths,config,pattern,h,data,eps,back,rootname,ctable,rtable
;
; INPUTS:
;	bac - reference background (from CALFOS_RDBAC).  It is used
;		only if no ystep has the background.
;	fwidths = 2 element vector with median filter width and 
;		the mean filter width
;	config - configuration vector (from CALFOS_RD)
;	pattern - pattern vector (from CALFOS_RD)
;	rootname - rootname of observation; if present and non-null,
;		   apply geomagnetic field correction factor to any
;		   reference backgrounds used.
;	ctable - geomagnetic field background corrections table.
;	rtable - geomagnetic field background rates reference table.
;
; INPUT/OUTPUTS:
;	h - FITS header
;	data - data array
;	eps - epsilon array
;
; OUTPUTS:
;	back - background array of size (ns x ysteps x slices*nreads)
; 
; HISTORY:
;	version 1.0  D. Lindler  Jan 1990
;	Feb 1991 D. Lindler   Fixed bug in zeroing of data for EPS>200
;	Jul 1991 D. Neill     Added implementation of geomagnetic background
;	Oct 1991 D. Neill     Added implementation of GMF rates reference table
;	Jan 1992 D. Neill     Changed BGF to GMF to comply with STSDAS
;-
;---------------------------------------------------------------------------
; extract pattern info
;
	ysteps = pattern(6)
	ns = pattern(8)
	slices = pattern(9)
	nreads = pattern(10)
;
; create array with type (OBJ = 3, SKY=2, and BCK = 1) for each ystep
;
	types = replicate(1,ysteps)
	if ysteps le 3 then begin
		ytypes = config(4:6)
		for i=0,ysteps-1 do begin
			case strtrim(ytypes(i)) of
				'OBJ' : types(i) = 3
				'SKY' : types(i) = 2
				'BCK' : types(i) = 1
				else  : types(i) = 0
			endcase
		endfor
	endif
;
; find background ystep number
;
	bck_ystep = where(types eq 1) & nbck=!err
	if nbck gt 0 then bck_ystep = bck_ystep(0)	;make it a scalar
;
; if requested, apply geomagnetic correction to reference background
;
	if n_elements(rootname) eq 0 then rootname=''
	use_gmf = (rootname ne '') and (nbck le 0) ; not if bkgnd taken
	if use_gmf then begin
		if strtrim(ctable,2) ne '' then begin
		  hist = 'Background correction factors from '+strtrim(ctable,2)
		  sxaddhist,hist,h
		endif else $
		if strtrim(rtable,2) ne '' then begin
		hist = 'Background corrections computed from '+strtrim(rtable,2)
		  sxaddhist,hist,h
		endif else begin
		  print,'CALFOS_BAC - GMF_CORR selected and no GMFTAB or',$
			'             CCS8 reference filename supplied'
		  retall
		endelse
		fos_back,rootname,ctable,rtable,bac,back
	endif else back = fltarr(ns,ysteps,slices*nreads)
;
; loop on all readouts and slices
;
	for iread = 0,nreads-1 do begin
	    for slice = 0,slices - 1 do begin
		ipos = iread*slices + slice

	    	if nbck gt 0 then begin		;background taken ?
		    bac = data(*,bck_ystep,ipos)	;yes, get it
		    eps_back = eps(*,bck_ystep,ipos)    ; get epsilons
		    message = ' using ystep = '+strtrim(bck_ystep+1,2)
;
; repair bad data in the background
;
		    mask = eps_back ge 200
		    repaired = where(mask) & n_repaired = !err
		    calfos_repair,mask,bac
		    for i=0,ysteps-1 do back(0,i,ipos) = bac	; output array.
		    if !err lt 0 then begin
			hist = 'All background data bad for readout ='+ $
				string(iread+1)+'   Slice = '+string(slice)
			sxaddhist,hist,h
			print,hist
		    endif

;
; smooth background
;
		    bac = calfos_median(bac,fwidths(0))
		    for k=1,2 do bac = calfos_mean(bac,fwidths(1))
;
; if no background taken, use reference background
;
		  end else begin

		    n_repaired = 0
		    if use_gmf then begin	;GMF corrected background
		    	baca = back(*,*,ipos)
			message = ' using GMF corrected ref. file backgrnd'
		    endif else begin		;uncorrected background
		    	for i=0,ysteps-1 do back(0,i,ipos) = bac
		    	message = ' using reference file background'
		    endelse

		endelse
;
; subtract background from the object and sky data
;
		for ystep = 0,ysteps-1 do begin
		    if types(ystep) gt 1 then begin	;object or sky?
		        eps1 = eps(*,ystep,ipos)
			d = data(*,ystep,ipos)
			if use_gmf then bac=baca(*,ystep)
			d = d-bac
			bad = where(eps1 ge 200) & nbad=!err
			if nbad gt 0 then d(bad) = 0
			data(0,ystep,ipos) = d
;
; flag repaired data points in the background
;
			if n_repaired gt 0 then begin
			    eps1(repaired) = eps1(repaired) > 120
			    eps(0,ystep,ipos) = eps1
			endif
		    endif
		endfor

	    endfor; loop on slices
	endfor; loop on readouts
;
; update history
;
	if nbck gt 0 then begin
		hist = 'Background filter widths  Median = '+ $
			strtrim(fwidths(0),2)+ '  Mean = '+strtrim(fwidths(1),2)
		sxaddhist,hist,h
		if !dump gt 0 then print,hist
	endif
	hist='Background subtraction completed'+message
	sxaddhist,hist,h
	if !dump gt 0 then print,hist

return
end