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