Re: systematic errors

Dear Mike,
  The method you suggested works but the plot does not show increased error bars due to systematic errors. Below is the script which I wrote to implement your method.

public define set_src_bkg_sys_err (spec_id)
{

%% Usage: set_src_bkg_sys_err(spec_id;fes=frac_sys_err_src,feb=frac_sys_err_bkg);
%%        set_src_bkg_sys_err(spec_id); ==> sys_err read from spectral files

% Define optional parameters frac_sys_err_src (fes) and frac_sys_err_bkg (feb)
    variable fes = qualifier ("fes",100.0);
    variable feb = qualifier ("feb",100.0);

    variable spec_data = get_data_info(spec_id);
    variable src_spec = spec_data.file ;
    src_spec = src_spec + "[SPECTRUM]";
    variable bkg_spec = fits_read_key(src_spec,"backfile");
% Get source counts
    variable c = get_data_counts (spec_id);
    c = c.value ;
    variable nbins_src = length(c);

% Get scaled background counts
    variable b = get_back(spec_id);
    if (nbins_src != length(b)) 
    {
	print("No. of bins in the source and background spectrum are not the same");
	break;
    }


% Define fractional systematic err arrays for the source (fd) and bkg (fb)
    variable fd = Double_Type[nbins_src];
    variable fb = Double_Type[nbins_src];

% Get the sys err fractions from the spectral files     
    if (fes==100.0)
    {
	if (fits_read_key(src_spec,"SYS_ERR") == 0)
	{
	    print("Systematic errors not specfied in the source spectrum\n");
	    print("Specify frac_sys_err for the source with the optional par fes\n");
	    break ;
	}
	else if (fes >=0.0)
% Get fractional sys error for the source	    
	    fd = get_sys_err_frac(spec_id);
    }
    else 
    fd = fd + fes ;
% Get fractional systematic error for the background spectrum
    if (feb == 100.0)
    {
	if (fits_read_key(bkg_spec+"[SPECTRUM]","SYS_ERR") == 0)
	{
	    print("Systematic errors not specfied in the background spectrum.\n");
	    print("Specify frac_sys_err for the background with the optional par feb\n");
	    break ;
	}
	else if (feb >= 0.0)
	fb = get_sys_err_frac(bkg_spec);
    }
    else 
    fb = fb + feb ;

% Calculate total fractional systematic error fraction
% To avoid division by zero, replace count = 0 by count = 1 in the 
% count array but not in the original spectrum. This is equivalent to setting 
% fractional sys err for background = 0 for bins where count is zero.
    
    c[[where(c == 0)]] = 1;
    variable tot_sys_err = sqrt(fd^2 + (fb * b/c)^2);
% Set the systematic error fraction
    set_sys_err_frac(spec_id,tot_sys_err);
    return ;
}
    


cheers,
-gulab	


----- "Michael Nowak" <mnowak_at_email.domain.hidden
> On Jul 31, 2009, at 8:52 AM, gulabd_at_email.domain.hidden> 
> > I have a situation where background level is determined with a  
> > systematic uncertainty of a few percent. So ideally I should  
> > increase the error bars associated with background by applying  
> > systematic errors to the background spectrum.
> >
> > XSPEC12 can read systematic errors from the background PHA file. I 
> 
> > can see the differences in the reduced chi^2 with and without the  
> > systematic errors applied to the background file with grppha.
> >
> > It seems ISIS can read systematic errors only from the source(+back)
>  
> > PHA file.
> 
> Who knew...  I have done things like vary the level of the background 
> 
> during the fits.  (I have actually been using in ISIS my own variation
>  
> of the XSPEC "corrfile" procedure, except that the normalization  
> constant could be varied during the fit, not iteratively between  
> fits.  Took XSPEC several years more before they changed their  
> functionality to allow that.)  But I've never tried to increase the  
> error bars.  But thinking out loud here, we have,
> 
> 	\Delta C^2 = C +  S_d/S_b * B_s + (F_d*C)^2 + (F_b*B_s)^2
> 
> where \Delta C^2 is the variance on the data counts, C is the counts, 
> 
> S_d is the data backscale*data exposure, S_b is the background  
> backscale*data exposure, B_s is the scaled background, F_d is the  
> fractional systematic errors on the data, F_b is the fractional  
> systematic errors on the background.  So, rewriting:
> 
> 	\Delta C^2 = C +  S_d/S_b * B_s + [ (F_d)^2 + (F_b*B_s/C)^2]*C^2
> 
> If you tell ISIS to set the *data* systematic errors to sqrt( F_d^2 + 
> 
> (F_b*B_s/C)^2 ), then you should get exactly the thing you want.
> 
> So, I know that we could do:
> 
> isis>  fd = get_sys_err_frac(data_id);
> isis> c = get_data_counts(data_id);
> isis> bs = get_back  % This is the scaled background
> 
> and that's three of the pieces you need.  I think if you define the  
> background as a file, e.g.,
> 
> isis> bid = load_data("bakcground.pha");
> isis> fb = get_sys_err_frac(bid);
> 
> Then you could do:
> 
> isis> set_sys_err_frac( data_id, sqrt(fd^2+ (fb*bs/c)^2));
> isis> delete_data(bid);  % You no longer need that as data, so drop
> it.
> 
> Now, since this is totally off the top of my head, I haven't tried  
> this at all.  You'd have to be careful about making sure all the  
> arrays were the same length, you were dividing through by 0's, etc., 
> 
> etc.  But with a little bit of effort, it seems like you could write a
>  
> dozen line script or so which would do this the way you wanted for any
>  
> data set.
> 
> Cheers,
> 
> Mike
----
You received this message because you are
subscribed to the isis-users list.
To unsubscribe, send a message to
isis-users-request_at_email.domain.hiddenwith the first line of the message as:
unsubscribe