;+ ; ; PROJECT: ; SDAC ; NAME: ; SPEX_BACKGROUND ; ; PURPOSE: ; controls the determination of background rate within the SPEX environment. ; ; CATEGORY: ; SPEX ; ; CALLING SEQUENCE: ; SPEX_BACKGROUND, Ut, Flux, Eflux, Rate, Erate, Back, Eback, Ltime, $ ; Tb, Wback, Com_param, Check_defaults, Command ; ; INPUTS: ; Ut - time array in seconds relative to utbase ( utbase=getutbase() ) ; for displayed light curve ; Flux - count rate flux, current units of displayed plot. ; background to be determined for each energy channel and ; time interval of flux ; Flux- uncertainties on flux ; Ltime- live time in each time bin ; Check_defaults - if set ask for background input ; ; OPTIONAL INPUTS: ; Com_param- command parameter from op_com ; used for passing in CLEAR command ; and channels for selective background ; Use 'null' to set background to zero for background subtracted data ; ; KEYWORD PARAMETERS: ; BACK_ORDER - optional keyword, polynomial order of background fit, ; 5 values, one for each band and all channel default ; ENERGY_BANDS - 2x4 energy intervals referenced to edges ; EDGES - nominal energy of output bands ; USE_BAND - set to 0,1,2,3 to use data in energy_bands(use_band) ; set to -1 to use a single background interval and fit order ; ; OUTPUTS: ; Rate - background subtracted count rate (also input) ; Erate- uncertainty on rate ; Back - background model for each channel and time bin ; Eback- uncertainty on back ; Wback- indices of time intervals used for background ; ; OPTIONAL OUTPUTS: ; Command- tells spex where to go after subtraction ; ; COMMON BLOCKS: ; BLOCK1: common spex_background, t_out_range ; ; SIDE EFFECTS: ; Describe "side effects" here. There aren't any? Well, just delete ; this entry. ; ; RESTRICTIONS: ; Describe any "restrictions" here. Delete this section if there are ; no important restrictions. ; ; PROCEDURE: ; Assumes event light curve is currently plotted in the active window ; Background is calculated over all of ut according to ; polynomial model of order back_order(def=1) using ; polyfitw over each energy channel in turn. Weighting is ; according to ltime, the livetime of the detector for ; each time bin. ; Can be used to clear background, "background, clear" in SPEX ; or can use existing background array, "background, old" in SPEX. ; After first fitting the background over the whole energy range, ; the background can be computed independently for any of the energy ; ranges specifed by energy bands in SPEX. That band is replotted alone ; without any background subtraction, and the current intervals are shown and ; the new intervals are requested and the fit is calculated according to ; the current value of back_order as usual for that energy range. ; The variable need_nw_back is not set to 1 until some background has ; been calculated for the whole energy range! ; ; EXAMPLE: ; Not to be used directly. Called from with SPEX, e.g. ; spex_background, ut, flux, eflux, rate, erate, back, eback, ltime, $ ; tb, wback, com_param, check_defaults, command, $ ; back_order=back_order, title=title, def_tb=def_tb, $; ; use_band= use_band, need = need, energy_bands=energy_bands, edges = edges ; ; MODIFICATION HISTORY: ; ras 14-apr-94 ; ras, 17-may-94, allows use of "old" background, already in SPEX ; ras, 23-july-94, integrated discrete energy band capability ; ras, 1-sep-94, fixed direct input problem ; ras, 27-sep-94, allow for no background subtraction with Null command ; ras, 20-oct-94, change to allow multi-channel ltime ; ras, pass eflux on to background fitter, 30-aug-95 ; Version 2, ras, 23-Mar-1996, tidied documentation header ; Version 3, ras, 27-Oct-1996, changed definition of nchan to match spex_intervals.pro ; ;- pro spex_background, ut, flux, eflux, rate, erate, back, eback, ltime, $ tb, wback, com_param, check_defaults, command, $ back_order=back_order, title=title, def_tb=def_tb, $; use_band=use_band, need = need, $ energy_bands=energy_bands, edges = edges common spex_background, t_out_range ;------------------------------------------------------------- background: ut_lim = limits(ut+getutbase()) size_flux = size(flux) nchan = size_flux(1) ;number of energy channels nbin = size_flux(2) ;number of time bins if strupcase(strmid(com_param(0),0,1)) eq 'C' then begin ;C stands for clearbackground rate = flux back = back*0.0 eback= eback*0.0 need.back = 1 ;the background rate has been cleared need.accum= 1 ;going to need new accumulation of data for fit intervals erate = eflux command(0) = 'graph' return endif ;Call background,old to use background already in SPEX if strupcase(strmid(com_param(0),0,1)) eq 'O' then begin ;O stands for old_background rate = flux - back erate = sqrt( eflux^2 + eback^2) ; AUTOMATICALLY DISPLAY RESULT OF BACKGROUND SUBTRACTION! command(0) = 'graph' need.back = 0 ;the background rate has been set need.accum= 1 return endif ;Call background,null if data input are background subtracted.' ;ras 27-sep-94 if strupcase(strmid(com_param(0),0,4)) eq 'NULL' then begin back = flux * 0.0 eback= back rate = flux erate = eflux need.back = 0 ;the background rate has been set need.accum= 1 return endif ; iselect - the index of each boundary relative to xdisplay ; for j=iselect(0,i) and k=iselect(1,i) the channels are ; j through k-1, inclusive ;Call background,auto to compute background in all bands for with each polynomial order if strupcase(strmid(com_param(0),0,4)) eq 'AUTO' and check_defaults eq 0 then begin temp=[-1,3] autoback = 1 endif else begin temp=[use_band,use_band] autoback = 0 endelse for iband = temp(0),temp(1) do begin if iband ne -1 then begin ebands = reform( energy_bands,2,n_elements(energy_bands)/2) scale = spex_current('scale_bands') scale_bands = scale * 0 scale_bands(iband) = scale(iband) command='graph' ;Dummy variables for spex_thistory using flux mode th_yrange=spex_current('th_yrange') th_ytype = spex_current('th_ytype') area = spex_current('area',/incommon) count_2_flux = spex_current('count_2_flux',/incommon) spex_thistory, ut, flux, eflux, command, spex_current('data_tipe'), edges, $ count_2_flux, spex_current(/incomm,'uflux'), area, title, $ strupcase(spex_current('t_hist_mode')), energy_bands, scale_bands, $ th_yrange, th_ytype, $ edg_units=spex_current('edg_units',/incomm), eplots=0, sampav=1, ltime=ltime command='' spex_intervals, edges, xeselect, xinput=ebands(*,iband), $ iselect=iselect, style=style, v_styles=['d'], query=0,$ error=error, input_variable='energy_bands' ;nchan = iselect(1,0) - iselect(0,0) ;old version, but spex_intervals has changed! nchan = iselect(1,0) - iselect(0,0) + 1 if nchan eq 0 then begin printx,'No channels in energy range. Select New Energy Bands.' return endif wrange = indgen(nchan)+iselect(0,0) endif else wrange = indgen(nchan) ; checkvar, tb, dblarr(2,6,5) ;start and stop, 6 intervals, overall & 4 channels ;tb is in seconds from 1-jan-79. Ut is in seconds from getutbase(). ;Set all values of tb outside of the ut range to zero. If one of ;the tb pairs are 0, then set the other to 0. if total(need.back) eq 1 then begin def_tb1 = def_tb(*,0) def_tb2 = def_tb(*,1) for ib = 0,4 do begin ;set defaults for each band if needed t_out_range = ut_lim(0)- 100.0 ;guaranteed out of range out_range = where( tb(0,*,ib) lt ut_lim(0) or tb(1,*,ib) gt ut_lim(1), n_orange) if n_orange ge 1 then begin tb(0,out_range,ib) = t_out_range tb(1,out_range,ib) = t_out_range endif ;Start with defaults tb_lim = where(tb(*,*,ib) gt t_out_range, ntb) if ntb lt 2 then begin tb(*,0,ib) = def_tb1 tb(*,1,ib) = def_tb2 endif ;In case defaults are also out of range of ut_lim out_range = where( tb(0,*,ib) lt ut_lim(0) or tb(1,*,ib) gt ut_lim(1), n_orange) if n_orange ge 1 then begin tb(0,out_range,ib) = t_out_range tb(1,out_range,ib) = t_out_range endif endfor endif graphics_page inrange = where( tb(0,*,iband+1) ge ut_lim(0) and tb(1,*,iband+1) le ut_lim(1), ninrange) if ninrange ge 1 then for i=0,ninrange-1 do for j=0,1 do $ oplot, tb(j,inrange(i),iband+1)+fltarr(2)-getutbase(), crange('y'), $ linest=2+j,color=fcolor(5+j) if strmid(strupcase(com_param(0)),0,1) ne 'D' or need.back then begin ;D stands for display widget = f_use_widget(/test) ;Are widgets enabled? yn ='Y' alpha_page case 1 of check_defaults and ninrange ge 1: begin if widget then begin ans=respond_widg(title='BACKGROUND MODEL FOR ENERGY_BAND: '+ $ arr2str(strtrim(limits(edges(*,wrange)),2),delim=', '),$ message='Are Background Intervals Acceptable?',$ buttons=['Yes','No '],/column) if ans eq 1 then yn = 'N' endif else read,'Are Background Intervals Acceptable?, (y/n)? ',yn end not check_defaults and ninrange ge 1: yn = 'Y' ;yes means process intervals already chosen else: yn='N' endcase if strupcase(strmid(yn,0,1)) eq 'N' then begin graphics_page if ninrange ge 1 then for i=0,ninrange-1 do for j=0,1 do $ oplot, tb(j,inrange(i),iband+1)+fltarr(2)-getutbase(), crange('y'), $ linest=2+j,color=fcolor(0) alpha_page spex_intervals, ut, tbselect, iselect=bselect, style='d',v_style=['d'],/query,color=9,$ xoffset= (indgen(3)*200+200),yoffset = (indgen(3)*200+200), $ message='Select Up to Six Discrete Background Intervals, '+ $ 'before and after flare, if possible.' endif else begin binput = tb(*,inrange,iband+1)-getutbase() spex_intervals, ut, tbselect, iselect=bselect, xinput=binput, $ style='d', v_style=['d'],query=0, color=9 endelse if n_elements(tbselect) mod 2 ne 0 then begin printx,['YOU HAVE SELECTED '+STRING(N_ELEMENTS(TBSELECT))+' BACKGROUND INTERVALS',$ 'YOU MUST SELECT TWO BOUNDARIES FOR EACH INTERVAL, AND AT LEAST ONE INTERVAL',$ 'RE-ENTERING BACKGROUND SELECTION PROCEDURE. TRY AGAIN!!!!!!!'] goto, background endif ;Reset the background intervals for that band! tb(*,*,iband+1) = t_out_range ;On the initial selection of the overall background model, set the defaults for all of ;the bands to the initial selection ;If in automatic mode for the global and four separate bands, then use the ;available background values in tb for each band! if iband eq -1 and not autoback then $ for ib=-1,3 do $ tb(0,0,ib+1) = reform(tbselect,2,n_elements(tbselect)/2,1) + getutbase() else $ tb(0,0,iband+1) = reform(tbselect,2,n_elements(tbselect)/2,1) + getutbase() wback = [0] for i=0,n_elements(bselect)-2,2 do $ wback=[wback,bselect(0,i/2)+indgen(bselect(1,i/2)-bselect(0,i/2)+1)] wback = wback(1:*) checkvar,back, flux*0.0 checkvar,eback, flux*0.0 if n_elements(back) ne n_elements(flux) then begin back = flux*0.0 & eback=back endif nrange = n_elements(wrange) back_sigma = fltarr(nrange) ;Find the background model for each channel. If the livetime for a ;sample is zero, then the background value is zero, too. for j=0,nrange-1 do begin back(wrange(j),*) = $ f_div( ltime(wrange(j),*) * $ fit_backgrnd(ut, flux(wrange(j),*), eflux(wrange(j),*), selected=wback, $ ltime=(ltime(wrange(j),*))(*), back_order(iband(0)+1), sigma=bsigma), $ ltime(wrange(j),*) ) back_sigma(j) = bsigma endfor rate = flux - back if iband eq -1 then need.back=0 ;the background rate has been set need.accum = 1 ;if background has changed, then re-accum intervals. eback(wrange,*) = rebin(reform(back_sigma,nrange,1),nrange,nbin ) erate = sqrt( eflux^2 + eback^2) ; AUTOMATICALLY DISPLAY RESULT OF BACKGROUND SUBTRACTION! command(0) = 'graph' endif endfor end