% file: isis_analysis.txt % % Use MARX to simulate an HETG observation of a powerlaw w/lines. % % - - - Perform spectral analysis % These steps are carried out in ISIS, running in the "working directory". % % [unix] isis % isis> % % These steps can be carried out by cutting and pasting from this file into ISIS; % or they can all be carried out by source'ing this file: % isis> .source isis_analysis.txt % % - - - - - % Set the simulation sub-directory: % point source: variable obs_dir = "obs_hplaw"; % disk source: %% variable obs_dir = "obs_hpldisk"; % Load the pha's and responses: ()=load_data(obs_dir+"/pha2",[9,10,3,4]); ()=load_arf(obs_dir+"/meg_-1.arf"); ()=load_arf(obs_dir+"/meg_1.arf"); ()=load_arf(obs_dir+"/heg_-1.arf"); ()=load_arf(obs_dir+"/heg_1.arf"); assign_arf(1,1); assign_arf(2,2); assign_arf(3,3); assign_arf(4,4); ()=load_rmf(obs_dir+"/meg_-1.rmf"); ()=load_rmf(obs_dir+"/meg_1.rmf"); ()=load_rmf(obs_dir+"/heg_-1.rmf"); ()=load_rmf(obs_dir+"/heg_1.rmf"); assign_rmf(1,1); assign_rmf(2,2); assign_rmf(3,3); assign_rmf(4,4); % - - - - - % Open a plotting window... variable onepane = open_plot("/XSERVE",1,1); % Group and notice the data: group_data(all_data, 8); % 0.6 to 7.5 keV: % (Looks like HEG is cutoff before 8 keV...) xnotice(all_data, Const_hc/7.5, Const_hc/0.6); % List the data list_data; % Use keV for the X axis in plots: plot_unit("keV"); % Use log x axis, linear y axis starting at 0. xlog; xrange(0.5,10.0); ylin; yrange(0.,); % Can look at the data title("MARX simulation of HETG observation"); plot_data_counts(1); % show all of these in grey oplot_data_counts(1,14); oplot_data_counts(2,14); oplot_data_counts(3,14); oplot_data_counts(4,14); % Ignore the line regions for the continuum fitting we'll do % Fe region ignore(all_data, 0.95*Const_hc/6.4, 1.05*Const_hc/6.4); % Ne region ignore(all_data, 0.97*Const_hc/0.8486, 1.03*Const_hc/0.8486); % % Replot the spectra again - use some colors: oplot_data_counts(1,2); % MEG -1 red oplot_data_counts(2,8); % MEG -1 orange oplot_data_counts(3,4); % HEG -1 blue oplot_data_counts(4,3); % HEG -1 green % - - - - - % Define the fitting function: fit_fun("phabs(1)*powerlaw(1)"); % Set initial values and limits for the parameters; % keep them all "thawed" with 0 as third argument. set_par (1, 1.00, 0, 0.1, 10.0); set_par (2, 0.01, 0, 1.e-4, 1.0); set_par (3, 1.00, 0, 0.5, 2.5); % % use "set_fit_statistic" and "set_fit_method" here % to define fitting specifics; help is available via: % isis> help set_fit_method % % Use Gehrels to help with low-counts bins set_fit_statistic ("chisqr;sigma=gehrels"); % Could use two fitting methods in sequence % if needed to get the best local minimum % Comment out the first fitting for this example. %% message(" - - - - -"); %% message("Do a minim fit ..."); %% set_fit_method ("minim;default"); %% ()=fit_counts; %% list_par; message(" - - - - -"); message("Do a lmdif fit ..."); set_fit_method ("lmdif;default"); ()=fit_counts; list_par; % For obs_hplaw : % list_par gave these results: % % phabs(1)*powerlaw(1) % idx param tie-to freeze value min max % 1 phabs(1).nH 0 0 0.7965323 0.1 10 10^22 % 2 powerlaw(1).norm 0 0 0.01999332 0.0001 1 % 3 powerlaw(1).PhoIndex 0 0 1.179763 0.5 2.5 message(" - - - - -"); % Show the data and residuals for the 4 data sets: variable fourpane = open_plot("/XSERVE",2,2); xlog; xrange(0.5,10.0); ylin; yrange(0.,); % don't show the error bars on main histogram window(fourpane); errorbars(10000); rplot_counts(1); rplot_counts(2); rplot_counts(3); rplot_counts(4); % - - - - - % Get parameter conf ranges % For help on conf do: isis> help conf % We could just do: % isis> conf(1,2); etc. % But the next lines make prettier output and % show some more S-Lang stuff that can be useful... variable ipar=0, clo, chi; _for ipar (1,3,1) { message("Finding conf for parameter "+string(ipar)+" ..."); (clo, chi)=conf(ipar, 2); message(" 99\% range is: "+string(clo)+" to "+string(chi)); }; % For obs_hplaw : % >>--> output from the confidence ranges are % in good agreement with the known input parameter % values of 0.80, 0.020, and 1.20 : % % Finding conf for parameter 1 ... % 99% range is: 0.754603 to 0.839679 % Finding conf for parameter 2 ... % 99% range is: 0.0189322 to 0.0211297 % Finding conf for parameter 3 ... % 99% range is: 1.13933 to 1.22082 message(" - - - - -"); % - - - - - % Now look closely at the narrow-line regions... % Use a finer binning notice(all_data); group_data(all_data, 2); ()=eval_counts; window(fourpane); xrange(0.799,0.901); xlin; rplot_counts(1); xrange(5.99,7.01); xlin; rplot_counts(3); xrange(0.799,0.901); xlin; rplot_counts(2); xrange(5.99,7.01); xlin; rplot_counts(4); % - - - - - % Fit the width of the Fe line % Plot the HEG and continuum model: window(fourpane); xrange(5.99,7.01); xlin; rplot_counts(3); xrange(5.99,7.01); xlin; rplot_counts(4); % Add a Gaussian to the fit function: fit_fun("phabs(1)*powerlaw(1)+gauss(1)"); % set the gaussian parameters: set_par("gauss(1).area", 0.001, 0, 1.e-8, 1.0); set_par("gauss(1).center", Const_hc/6.40, 0, Const_hc/6.7, Const_hc/6.2); set_par("gauss(1).sigma", 0.001, 0, 0.0001, 0.1); % Notice just the HEG spectra in the range: xnotice(3, Const_hc/7.01,Const_hc/5.99); xnotice(4, Const_hc/7.01,Const_hc/5.99); % ignore the MEG spectra: ignore([1,2]); % Freeze the continuum nH and PhoIndex, leaving norm free: freeze([1,3]); ()=fit_counts; list_par; % For obs_hplaw : % list_par gave these results: % % phabs(1)*powerlaw(1)+gauss(1) % idx param tie-to freeze value min max % 1 phabs(1).nH 0 1 0.7965323 0.1 10 10^22 % 2 powerlaw(1).norm 0 0 0.01900418 0.0001 1 % 3 powerlaw(1).PhoIndex 0 1 1.179763 0.5 2.5 % 4 gauss(1).area 0 0 0.0002898555 1e-08 1 ph/s/cm^2 % 5 gauss(1).center 0 0 1.936372 1.85051 1.9997 A % 6 gauss(1).sigma 0 0 0.009638465 0.0001 0.1 A % and show plots with the Gaussian in the model: xrange(5.99,7.01); xlin; rplot_counts(3); xrange(5.99,7.01); xlin; rplot_counts(4); % Get conf ranges on the Gaussian parameters; % allow the continuum to be adjusted as well. _for ipar (4,6,1) { message("Finding conf for parameter "+string(ipar)+" ..."); (clo, chi)=conf(ipar, 2); message(" 99\% range is: "+string(clo)+" to "+string(chi)); }; % For obs_hplaw : % output of conf ranges for Gaussian fit of Fe line, % known input values were: 0.0003, 1.9361, and 0.0085. % % Finding conf for parameter 4 ... % 99% range is: 0.000206056 to 0.000378273 % Finding conf for parameter 5 ... % 99% range is: 1.93279 to 1.93999 % Finding conf for parameter 6 ... % 99% range is: 0.00616031 to 0.0136647 % Convert the sigma and center to a velocity FWHM: variable line_center = get_par(5); message(" Line width, c*dE(FWHM)/E, is : " + string(Const_c*1.e-5*(2.35*get_par("gauss(1).sigma"))/line_center)+" km/s." ); % Convert chi, clo (these will be from the last parameter which was sigma) % into a line width, FWHM in km/s. % If the clo is very small then set it to 0: if (clo < 0.0005) clo = 0.0; message(" Line width, c*dE(FWHM)/E, is between: " + string(Const_c*1.e-5*(2.35*clo)/line_center)+ " and "+string(Const_c*1.e-5*(2.35*chi)/line_center)+" km/s." ); % For obs_hplaw we get close to the 3100 km/s expected value: % Line width, c*dE(FWHM)/E, is : 3506.77 km/s. % With 99% conf range: % Line width, c*dE(FWHM)/E, is between: 2241.31 and 4971.65 km/s. %......... % If this analysis is run on the spatially extended simulation, obs_hpldisk, % then a double-peaked line is seen due to the spatial extent of the % source and the Gaussian fit to it gives a sigma value much larger than % the expected 0.0085 A : % 6 gauss(1).sigma 0 0 0.0209564 <--- % and in velocity FWHM: % Line width, c*dE(FWHM)/E, is : 7623.81 km/s. % Line width, c*dE(FWHM)/E, is between: 5568.2 and 10607.4 km/s. % % - - - end of demo - - -