;+ ; Name: spex_bk ; ; Purpose: Used in spex object to calculate background and error for each raw energy/time bin from one of ; several bk methods. ; ; 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 ; 13-Sep-2006, Kim. In plot_setup, call spex_apply_eff for errors too if photons ; 18-Sep-2006, Kim. Lots of changes in process method to implement bk_ratio method. ; 11-Jul-2007, Kim. In process, use midpoints of energy bins to figure out which ; sep_bk band bin belongs in (and which order to use). Previously used start/end ; and would sometimes miss bins (so they'd get assigned the order from lowest band) ; 15-Dec-2008, Kim. In process, added exp method for fitting background (order=4) ; 30-Dec-2008, Kim. In process, in exp method, handle 2 special cases - when no positive ; points in intervals for current energy band - leave bk as 0. When only one pos. point, ; set bk to that value for current energy band. ; 15-Jan-2009, Kim. In process bk_ratio option, if no time intervals for highest band, ; abort (set bk to -1). If no time intervals set for any other bands, set bk for that ; band to 0. (previously crashed) ; 30-Jul-2009, Kim. In process, fit_backgrnd now takes care of exp fits correctly, just ; call it with /exp. Previously for exp, fit 1poly to log of data>0., but this biased ; the background high. ; 7-Dec-2009, Kim. Added bin_data method to add flag in call to generic bin_data method in cases ; where error in binned background is computed differently from counts data. Use spex_bk_poisson_error ; control param to decide whether /back_ave should be set in call to bin_data. ; 18-Dec-2009, Kim. in process, ratio method, take ratio out of sqrt in error calculation, ; so it's sqrt (counts in profile). Use ratio times the error in the actual counts ; that were used in the profile. Also set spex_bk_poisson_error to 1 if using ratio method, ; to 0 if not using ratio method. ; 08-Jan-2010, Kim. Wow. In previous change, was overwriting user's selection of spex_bk_poisson_error. ; Want to ensure that if ratio method is used, poisson statistical errors are used. ; Added new info parameter - spex_bk_used_ratio - which will be 1 if ratio method was used. Now ; don't change spex_bk_poisson_error depending on method. Instead, in bin_data use averaging method ; rather than poisson method when spex_bk_poisson_error is 0 AND spex_bk_used_ratio is 0 ; 30-Jul-2013, Kim. Made bk_ratio method more robust. ; 1.Previously: To create high band profile, interpolated between profile in selected time intervals, replaced ; first and last two points with average of first/last 11 points, then smoothed. This had two problems: interpolation ; was based on single points at end/start of time intervals - if those points were high or low, threw off interpolated ; points. Also, when taking average of first/last 11 points, didn't ensure that they were all in first/last time ; intervals. Then convolve with savgol filter. ; Now: Find groups of contiguous points (i.e. the points in each specified bk time interval) and smooth them before doing ; the interpolation between groups - use boxcar smoothing with width=10 (or less if fewer points in group). Also when ; replacing first and last two points with average use up to 10 points, but fewer if there are <10 points in first/last ; group. After interpolation on smooothed data replace non-interpolated points with originaldata, then convolve with ; savgol filter. ; 2.Previously, error on high-band profile was wrong - just used error that came from the bin_data command. ; Now, treat error the same way we're treating flux profile - smooth, interpolate, smooth, etc. ; 3.If half width to smooth is 1 (means 3 points: 1left +1right+1point) set degree in savgol to 2 so doesn't reject. ; And if half width to smooth is 0, skip smoothing with savgol. ; 4.Added spex_bk_ratprofile_plot info param - if set then plot high band profile used in bk_ratio method. ; 12-Aug-2013, Kim. For spex_bk_ratprofile_plot plot, use ylog and find range that's reasonable (since might have 0.) ; 23-Aug-2013, Kim. Fixed bug in bk_ratio method - didn't work for just one time interval in high band. ; 05-Sep-2013, Kim. Made bk_ratio method work differently. Now can choose to set some ebands to use polynomial fit, others ; use ratio. order=5 means ratio for a band. If ratio selected for anything, does ratio method first for all bands, ; then overwrites any energy bin bk with fit to poly or exp if selected. In set method, if setting bk_ratio, set order=5 for ; all bands, and if unsetting bk_ratio, set any bands with order=5 to order=0. ; 17-Sep-2013, Kim. Now each energy band for defining bk can use any of these options: ; 0poly,1poly,2poly,3poly,exp,high e profile,this e profile. 'This E Profile' is new - same idea as ratio to high-energy band profile, ; except the profile is created using this eband. For all of these methods, the error should be the sqrt of counts in the time interval(s) ; used to calculate ratio (previously for high-eband ratio, we used sqrt of counts calculated for each bin, and when combining time bins, ; used sqrt of total counts in new bin), and when combining time bins, errors should be averaged. Got rid of spex_bk_used_ratio flag ; (and check of spex_bk_used_ratio we were previously doing in bin_data to decide whether to set the ave_bck flag). Also, moved code ; that computes profile to separate routine called spex_bk_profile. ; 24-Jan-2014, Kim. If haven't set bk_sep, set bk_eband to minmax(energies) just in case a profile method is selected ; 22-Feb-2019, Kim. In process, for order lt 5, don't do fit if all data is 0. - returns NaNs. Leave 0.s for that bin. ; ; ;- ;--------------------------------------------------------------------------- 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, $ spex_bk_ratio=spex_bk_ratio, $ 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 if exist(spex_bk_ratio) then begin ;If turning ratio method on, set all en band's orders to 5 (ratio). ;If turning it off, only turn the ones that were 5 (ratio) to 0 (0poly). bk_order = self->get(/spex_bk_order) if spex_bk_ratio then bk_order = bk_order*0. + 5. else begin q = where(bk_order eq 5, count) if count gt 0 then bk_order[q] = 0. endelse endif self->framework::set,spex_bk_ratio=spex_bk_ratio, 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 ;-------------------------------------------------------------------------- ; spex_bk class needs its own bin_data method so it can add the ave_arr flag if necessary. ; When bk was computed from a fit based on data in selected time intervals, then bk points are ; not independent from each other, and the errors from the fit should be averaged when ; binning across time. Otherwise when you combine time intervals, the errors get smaller ; and smaller, and they shouldn't. However, if bk is from real data points (e.g. if user inserted ; real data from previous day as bk), then the error in bk is computed ; the same way it is for obs data, ; by adding the errors in quadrature (which is the same as using the sqrt of counts in the ; combined intervals). ; Note: it is up to the user to set spex_bk_poisson_error to 1 if they insert ; bk data. The default is 0. ; Modified 09/20/2013, Kim. to average error in all cases (even ratio method) except when ; spex_bk_poisson_error is 1, because now even for ratio method, bk error is based on counts in ; time interval used to compute ratio, so now want to average errors (otherwise error gets unreasonably ; small when combine time intervals) function spex_bk::bin_data, _ref_extra=_extra poisson_error = self->get(/spex_bk_poisson_error) eq 1 ; or self->get(/spex_bk_used_ratio) eq 1 ave_err = ~poisson_error ;print,'ave_err=',ave_err return, self->spex_gen::bin_data(_extra=_extra, ave_err=ave_err) return, result 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) enmid = get_edges(energies,/mean) tmid = get_edges(times, /mean) bk_sep = self -> get(/spex_bk_sep) bk_eband = bk_sep ? self -> get(/spex_bk_eband) : minmax(energies) neband = bk_sep ? n_elements(bk_eband)/2 : 1 bk_index = self -> getdata(class='spex_bkint') have_bk = self -> get(/spex_have_bk) self -> set, spex_bk_used_ratio = 0 err_msg = '' if have_bk then begin bkcounts=data.data * 0. ebkcounts = data.edata * 0. nen = n_elements(data.data[*,0]) ; use get with this... keywords instead of just get(spex_bk_order) because if we added energy bands ; it will create new values, whereas the get(spex_bk_order) just gets whatever was last set bk_order = intarr(neband) for ib=0,neband-1 do bk_order[ib] = self -> get(this_band=ib, /this_order) ; If any bands use the 'ratio to high' method (order=5), compute that hi_profile just once here. ; For energies that use the 'ratio to high' method, base the time profile of the background ; in each energy bin on time profile of highest bk eband. if is_member(5, bk_order) ? 1 : 0 then begin ; do any bands use the ratio to high band profile method? iband = neband-1 p_eband = bk_eband[*,iband] ; use spex_gen::bin_data, so it won't call spex_bk::bin_data, since this isn't bk data yet that we're binning raw_profile = reform(self -> spex_gen::bin_data (data=data, intervals=p_eband, units_str=units_str)) hi_profile = spex_bk_profile(data=raw_profile, $ tmid=tmid, $ sel=*bk_index[iband], $ sm_width=self -> get(/spex_bk_sm_width), $ plot=self->get(/spex_bk_ratprofile_plot), $ eband=p_eband, $ err_msg=err_msg) if err_msg ne '' then message, /cont, err_msg + '. Setting ratio method energy bands to 0.' ; if there was a problem, hi_profile will be -1 endif ; For each raw energy bin, find what energy band it's in, and use the order for that band to decide how to compute ; background. order = 0 to 4 means use functions 0,1,2,3poly and exp, order=5 means use ratio to high eband profile, ; order=6 means use ratio to current band profile. The profiles are computed in spex_bk_profile for the broad eband. ; For each raw energy bin, use the points in the time interval(s) selected to either fit the function or find the ratio ; to the corresponding points in the profile (ratio is average of points in current energy bin to average of corresponding ; points in profile, profile for this bin computed by ratio*profile). ; Compute error by sqrt of total counts in this e bin in the time interval(s) used to fit function or compute ratio. ; (Note: this type of error should be averaged when combining time bins, not summed in quadrature) ; Default is to use 0th bk_time (indices in bk_index[0]) and bk_order, but if bk_sep, then for each bin will ; set sel and order for that energy. sel = *bk_index[0] order = bk_order[0] iband = 0 for i = 0,nen-1 do begin ; loop over raw energy bins ; if different time intervals and order for different energies, get the sel and order for this energy if bk_sep then begin iband = where (enmid[i] ge bk_eband[0,*] and enmid[i] lt bk_eband[1,*], count) iband = count gt 0 ? iband[count-1] : 0 sel = *bk_index[iband] order = bk_order[iband] ;print,i, order, enmid[i], energies[*,i] endif if sel[0] ne -1 then begin case 1 of ; fit to exponential or polynomial of order 'order' order lt 5: begin do_exp = order eq 4 ? 1 : 0 ; If all data for this energy bin is 0., then don't do fit - will return NaNs ; Just leave bkcounts and ebkcounts set at their initial value of 0. if ~same_data(minmax(data.data[i,*]), [0.,0.]) then 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,*], exp=do_exp) ebkcounts[i,*] = sigma endif end ; Use ratio * hi_profile order eq 5: begin if hi_profile[0] ne -1 then begin ratio = average(data.data[i,sel]) / average(hi_profile[sel]) ;print,'i, e[i], iband, ratio, minmax(sel) ', i, energies[0,i], iband, ratio, minmax(sel) bkcounts[i,*] = ratio * hi_profile ; error is just based on # counts in sel ebkcounts[i,*] = f_div(sqrt(total(data.data[i,sel]*data.ltime[i,sel])), total(data.ltime[i,sel])) endif end ; Use ratio * this eband's profile (this_profile) order eq 6: begin p_eband = bk_eband[*,iband] ; Check if we've already computed profile for this p_eband if ~same_data(p_eband, profile_eband) then begin ; use spex_gen::bin_data, so it won't call spex_bk::bin_data, since this isn't bk data yet that we're binning raw_profile = reform(self -> spex_gen::bin_data (data=data, intervals=p_eband, units_str=units_str)) this_profile = spex_bk_profile(data=raw_profile, tmid=tmid, sel=sel, $ sm_width=self -> get(/spex_bk_sm_width), $ plot=self->get(/spex_bk_ratprofile_plot), $ eband=p_eband, err_msg=err_msg) if err_msg ne '' then message, /cont, err_msg + ' Setting profile method energy bands to 0.' ; if there was a problem, this_profile will be -1 profile_eband = p_eband ; save eband that this_profile is for, so won't remake until bin i is in new eband endif if this_profile[0] ne -1 then begin ratio = average(data.data[i,sel]) / average(this_profile[sel]) ;print,'i, e[i], iband, ratio, minmax(sel) ', i, energies[0,i], iband, ratio, minmax(sel) bkcounts[i,*] = ratio * this_profile ; error is just based on # counts in sel ebkcounts[i,*] = f_div(sqrt(total(data.data[i,sel]*data.ltime[i,sel])), total(data.ltime[i,sel])) endif end endcase endif ; else if sel was -1, do nothing since bk and ebk for this bin are already set to 0 endfor ; computed bk as rate bk = {data: bkcounts, edata: ebkcounts, ltime: data.ltime} ; 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} return 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) z_err = transpose(z_err) 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 begin spex_apply_eff, ydata, drm_eff, units_str=units_str spex_apply_eff, z_err, drm_eff endif 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