%; Time-stamp: <2001-10-23 19:45:40 dph> %; MIT Directory: ~dph/libisis %; File: Contin_tutorial.sl %; Author: D. Huenemoerder %; Original version: 2001.10.22 %;==================================================================== % % version: 1 % % PURPOSE: % % Demonstrate continuum fitting in isis, using APED continuum data % and user-defined model. % % Also demostrate definition of thermal plasma model for use in % defining bright lines. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % aped_contin.sl defines a user model which evaluates the APED % continuum for standard abundances and a single % temperature component. % Tell slang where these functions come from: % autoload("fast_aped_contin_fit", "aped_contin.sl"); autoload("aped_contin_fit", "aped_contin.sl"); autoload("fast_aped_contin_init", "aped_contin.sl"); autoload("line_free_regs", "line_free_regs.sl"); % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % set up a plot window: % p2 = open_plot("/xwin"); resize(20,0.6); multiplot([3,1]); p1 = open_plot("/xwin"); resize(20,0.6); %%%%% Set up the data and responses: fpha="./Pha/acisf00009_003N001_pha2.fits"; % Standard Type II Pha file. hm=3; hp=4; mm=9; mp=10; % define HEG, MEG order row variables. h=load_data(fpha,[hm, hp, mm, mp]); % load HEG, MEG +-1st order from % standard product list_data; % show the histogram info. % Define ARF files string arrays: % fgarfs_heg = [ "Arf/acisf00009_003N001HEG_-1_garf.fits", "Arf/acisf00009_003N001HEG_1_garf.fits" ]; fgarfs_meg = [ "Arf/acisf00009_003N001MEG_-1_garf.fits", "Arf/acisf00009_003N001MEG_1_garf.fits" ]; % Define and load RMF files: (also in CALDB) % frmf_heg = "Ard/acisheg1D1999-07-22rmfN0004.fits"; frmf_meg = "Ard/acismeg1D1999-07-22rmfN0004.fits"; % Load rmfs, define indices, define number of orders % load_rmf(frmf_heg); rmf_heg=1; nheg = length(fgarfs_heg); load_rmf(frmf_meg); rmf_meg=2; nmeg = length(fgarfs_meg); % Load ARFS, assign ARF, RMF to histograms: % for (i=0; i < length(fgarfs_heg); i++) % loop is overkill here, but % demonstrates how to generalize % to many files. { ()=load_arf(fgarfs_heg[i]); assign_arf(i+1,i+1); assign_rmf(rmf_heg, i+1); } for (i=0; i < length(fgarfs_meg); i++) { ()=load_arf(fgarfs_meg[i]); assign_arf(nheg + i+1, nheg + i+1); assign_rmf( rmf_meg, nheg + i+1); } list_data; % list histogram info to show assignments % read the plasma database, to provide continuum emissivity info. % plasma(aped); %%%%%%%%% done w/ basic setup of data and responses. %%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Define the model; single component to start. % Later, we will try adding another component to see if it matters. evalfile("aped_contin.sl"); % load the model functions. fit_fun("aped_contin(1)"); % model is thermal plasma continuum list_par; % to show the model parameters % set some likely defaults: freeze % n value /thaw min max set_par(1, 0.004, 0, 0, 0.1); set_par(2, 2.3e7, 0, 1.e4, 5.e8); % Group data coarsely, since this is a continuum fit, and we would % like good statistics. meg_group_factor = 8; heg_group_factor = 2*meg_group_factor; % redefine order indices to map to the loaded histograms: mm = 3; mp = 4; hm = 1; hp = 2; group_data(mm, meg_group_factor); group_data(mp, meg_group_factor); group_data(hm, heg_group_factor); group_data(hp, heg_group_factor); %%%%%%%%%%%%%%%%%%%%% % % First, use only line-free regions. Function applies ignore region table. line_free_regs([hm,hp,mm,mp]); % For first pass, single high-T component, only notice the % short-wavelength region, and fit both HEG and MEG: ignore([hm,hp,mm,mp], 10, 25); list_data ; % check yrange(0); xrange(1,8); plot_data_counts(mm,1); oplot_data_counts(mp,2); oplot_data_counts(hm,3); oplot_data_counts(hp,4); % Do a fit, list the params, compute confidence: fit_counts; list_par; % Parameters[Variable] = 2[2] % Data bins = 342 % Best fit chi-square = 50.35 % Reduced chi-square = 0.1481 % aped_contin(1) % idx param tie-to freeze value min max % 1 aped_contin[1].norm 0 0 0.02229772 0 0.1 % 2 aped_contin[1].T 0 0 1.835275e+07 10000 5e+08 oplot_model_counts(mm,5); oplot_model_counts(mp,6); oplot_model_counts(hm,7); oplot_model_counts(hp,8); % evaluate 90% confidence intervals: % (clo_norm, chi_norm) = conf(1); ()=printf("Conf norm = %f , %f\n", clo_norm, chi_norm); (clo_T, chi_T) = conf(2); ()=printf("Conf T = %e , %e\n", clo_T, chi_T); % Conf norm = 0.020989 , 0.023607 % Conf T = 1.713500e+07 , 1.957111e+07 % Hardcopy: % pid = dup_plot("Contin_fit_0.gif/gif",p1); resize(6*2.54, 0.6); plot_data_counts(mm,1); oplot_data_counts(mp,2); oplot_data_counts(hm,3); oplot_data_counts(hp,4); oplot_model_counts(mm,5); oplot_model_counts(mp,6); oplot_model_counts(hm,7); oplot_model_counts(hp,8); close_plot(pid); window(p1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%% If satisfied, freeze this, and add a new component, and %%%%%%%%%% fit long-wavelength region where low-T has more affect. freeze([1,2]); fit_fun("aped_contin(1) + aped_contin(2)"); % insert another T component set_par(3, 0.004, 0, 0, 0.1); % set initial params set_par(4, 2.3e6, 0, 1.e4, 5.e8); list_par; ignore([hm,hp]); % no good region for HEG - ignore it. notice([mm,mp]); % notice the MEG, then apply regions line_free_regs([mm,mp]); % line free, ignore([mm,mp], 1, 8); % ignore short waves. yrange(0); xrange(8,25); plot_data_counts(mm,1); oplot_data_counts(mp,2); fit_counts; list_par; % Parameters[Variable] = 4[2] % Data bins = 328 % Best fit chi-square = 21.81 % Reduced chi-square = 0.0669 % aped_contin(1) + aped_contin(2) % idx param tie-to freeze value min max % 1 aped_contin[1].norm 0 1 0.02229772 0 0.1 % 2 aped_contin[1].T 0 1 1.835275e+07 10000 5e+08 % 3 aped_contin[2].norm 0 0 0.002478395 0 0.1 % 4 aped_contin[2].T 0 0 3.381622e+08 10000 5e+08 oplot_model_counts(mm,5); oplot_model_counts(mp,6); (clo_norm2, chi_norm2) = conf(3); ()=printf("Conf norm = %f , %f\n", clo_norm2, chi_norm2); % (clo_T2, chi_T2) = conf(4); ()=printf("Conf T = %e , %e\n", clo_T2, chi_T2); % **** Parameter range endpoint 0 is inside the confidence limit % **** Lower confidence limit didn't converge[3]: allow wider parameter ranges? % Conf norm = 0.000000 , 0.006798 % This example spectrum gives a very high temperature and very low % normalization, which means we probably do not need a low temperature % component. % % To actually see the two components, we need to evaluate each one and % "get" the array: fit_fun("aped_contin(1)"); % pick the first model component; ignore([hm,hp,mp]); % only MEG-1 notice(mm, 1.7, 25); % whole range; group_data(mm,0); % full resolution c = get_data_counts(mm); % get the data into a variable eval_flux; % evaluate the model in flux space s1=get_model_flux(mm); % get into a variable % Repeat for other component: % fit_fun("aped_contin(2)"); % pick the first model component; eval_flux; % evaluate the model in flux space s2=get_model_flux(mm); % get into a variable % examine them: (Weak features are pseudo-continuum - sum of weak lines.) % label("Wavelength [Angstroms]", "Flux [phot/cm^2/s/bin]", "Continuum Fit Tests"); xrange(1.7,25); hplot(s1.bin_lo, s1.bin_hi, s1.value,2); ohplot(s2.bin_lo, s2.bin_hi, s2.value,3); % make a "hard"copy: % pid = dup_plot("Contin_fit_1.gif/gif",p1); resize(6*2.54, 0.6); hplot(s1.bin_lo, s1.bin_hi, s1.value,2); ohplot(s2.bin_lo, s2.bin_hi, s2.value,3); close_plot(pid); window(p1); %%%%%%%%% CONCLUSION: only need component 1. % % Now we need to evaluate the desired model on the full MEG grid and % save store it for use. We will also compare it to the spectrum. fit_fun("aped_contin(1)"); % reset to component 1. save_par("acisf00009_003N001MEG_1_Contin_Model.par"); % save par file to disk. % we've already noticed the MEG+1 order, full region, eval_counts; % so fold the model: % look at it: % yrange(0,20);xrange(5,7);plot_data_counts(mm,2);oplot_model_counts(mm,3); % save a copy: (i.e., repeat to different device) % pid = dup_plot("Contin_fit_2.gif/gif",p1); resize(6*2.54,0.6); yrange(0,20);xrange(5,7);plot_data_counts(mm,2);oplot_model_counts(mm,3); close_plot(pid); window(p1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%% Make Detail plot, to see how contin follows spectrum. % group_data(mm,8); group_data(mp,8); group_data(hm,16); group_data(hp,16); xnotice([hm,hp], 1.7, 17); xnotice([mm,mp], 1.7, 25); xrange(1.7, 25); eval_counts; %%% since many pages, output to single postscript file: pid = dup_plot("Contin_fit_3.ps/cps", p1); resize(10*2.54, 0.6); label("Wavelength [Angstroms]", "Flux [counts/bin]", "Continuum Fit Check MEG-1"); yrange(0); xrange(1.7, 6); plot_data_counts(mm,2); oplot_model_counts(mm,3); xrange(5,11); plot_data_counts(mm,2); oplot_model_counts(mm,3); xrange(10, 16); plot_data_counts(mm,2); oplot_model_counts(mm,3); xrange(15, 21); plot_data_counts(mm,2); oplot_model_counts(mm,3); xrange(20, 26); plot_data_counts(mm,2); oplot_model_counts(mm,3); label("Wavelength [Angstroms]", "Flux [counts/bin]", "Continuum Fit Check MEG+1"); xrange(1.7, 6); plot_data_counts(mp,2); oplot_model_counts(mp,3); xrange(5,11); plot_data_counts(mp,2); oplot_model_counts(mp,3); xrange(10, 16); plot_data_counts(mp,2); oplot_model_counts(mp,3); xrange(15, 21); plot_data_counts(mp,2); oplot_model_counts(mp,3); xrange(20, 26); plot_data_counts(mp,2); oplot_model_counts(mp,3); label("Wavelength [Angstroms]", "Flux [counts/bin]", "Continuum Fit Check HEG-1"); xrange(1.7, 6); plot_data_counts(hm,2); oplot_model_counts(hm,3); xrange(5,11); plot_data_counts(hm,2); oplot_model_counts(hm,3); xrange(10, 16); plot_data_counts(hm,2); oplot_model_counts(hm,3); xrange(15, 21); plot_data_counts(hm,2); oplot_model_counts(hm,3); label("Wavelength [Angstroms]", "Flux [counts/bin]", "Continuum Fit Check HEG+1"); xrange(1.7, 6); plot_data_counts(hp,2); oplot_model_counts(hp,3); xrange(5,11); plot_data_counts(hp,2); oplot_model_counts(hp,3); xrange(10, 16); plot_data_counts(hp,2); oplot_model_counts(hp,3); xrange(15, 21); plot_data_counts(hp,2); oplot_model_counts(hp,3); close_plot(pid); window(p1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%% SETUP FAST CONTIN model: %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % fast_aped_contin_init(mm); % Store it in a static variable for use by... %%%%%%%%%%%% DEFINE model continuum + line % fit_fun("fast_aped_contin(1) + gauss(1)" ); % contin + gaussian model list_par; % fast_aped_contin(1) + gauss(1) % idx param tie-to freeze value min max % 1 fast_aped_contin[1].norm 0 0 0 0 0 % 2 gauss[1].area 0 0 1 0 0 % 3 gauss[1].center 0 0 12 0 0 % 4 gauss[1].sigma 0 0 0.025 0 0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5 % % Now let's fit a line, HEG+MEG, +-1: group_data(hm,0); group_data(hp,0); group_data(mm,0); group_data(mp,0); % undo grouping xnotice([hm,hp,mm,mp], 11.9,12.25); % notice the region xrange(11.9,12.3); yrange(0); % set the scales % % Take a look: % label("Wavelength [Angstroms]", "Flux [counts/bin]", "Continuum+Line Fit Tests"); plot_data_counts(mm,1);oplot_data_counts(mp,2);oplot_data_counts(hm,3);oplot_data_counts(hp,5); % set up params: % set_par(1, 1.0, 1, 0, 2); % set contin normalization to 1.0, freeze it. set_par(2, 0.005, 0, 0, 0.1); % wild guess at a line flux. set_par(3, 12.135, 0, 12.11, 12.16); % Ne X + Fe XX set_par(4, 0.005, 0, 1.e-6, 0.01); % allow free line width. % % Do the fit: % fit_counts; list_par; % fast_aped_contin(1) + gauss(1) % idx param tie-to freeze value min max % 1 fast_aped_contin[1].norm 0 1 1 0 2 % 2 gauss[1].area 0 0 0.0006294898 0 0.1 % 3 gauss[1].center 0 0 12.13422 12.11 12.16 % 4 gauss[1].sigma 0 0 0.006828275 1e-06 0.01 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% SOME PLOTS: % % over-plot model: (adjusting colors to make model easier to see...) % plot_data_counts(mm,1);oplot_data_counts(mp,2);oplot_data_counts(hm,3);oplot_data_counts(hp,5); oplot_model_counts(mm,3);oplot_model_counts(mp,5);oplot_model_counts(hm,1);oplot_model_counts(hp,2); % make a gif file: % pid = dup_plot("Contin_fit_4.gif/gif",p1); resize(6*2.54, 0.6); plot_data_counts(mm,1);oplot_data_counts(mp,2);oplot_data_counts(hm,3);oplot_data_counts(hp,5); oplot_model_counts(mm,3);oplot_model_counts(mp,5);oplot_model_counts(hm,1);oplot_model_counts(hp,2); close_plot(pid); window(p1); % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % make residual plots: % NOTICE that the HEG shows large residuals. This line is actually a % blend. There are 2 Ne X components separated by 0.005A (12.132, % 12.137), and there is weak Fe XXII-XXIII (12.136, 12.161) window(p2); yrange; % default to min, max. xrange(11.9,12.3); % set the scales rplot_counts(mm); % MEG -1 rplot_counts(mp); % MEG +1 rplot_counts(hm); % HEG -1 rplot_counts(hp); % HEG +1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% RESIDUAL plots, to gif files... % make gif files: pid = dup_plot("Contin_fit_5.gif/gif",p2); resize(6*2.54, 0.6); title("MEG -1"); rplot_counts(mm); % MEG -1 close_plot(pid); window(p2); pid = dup_plot("Contin_fit_6.gif/gif",p2); resize(6*2.54, 0.6); title("MEG +1"); rplot_counts(mp); % MEG +1 close_plot(pid); window(p2); pid = dup_plot("Contin_fit_7.gif/gif",p2); resize(6*2.54, 0.6); title("HEG -1"); rplot_counts(hm); % HEG -1 close_plot(pid); window(p2); pid = dup_plot("Contin_fit_8.gif/gif",p2); resize(6*2.54, 0.6); title("HEG +1"); rplot_counts(hp); % HEG +1 close_plot(pid); window(p2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% window(p1) ; % reset device. %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%% REPEAT portion in region w/ stronger contin: Mg XII %%% xnotice([hm,hp,mm,mp], 8.33, 8.53); xrange(8.2,8.7); plot_data_counts(mm,1);oplot_data_counts(mp,2);oplot_data_counts(hm,3);oplot_data_counts(hp,5); set_par(2,0.0003, 0, 0, 0.01); set_par(3, 8.42, 0, 8.4, 8.45); set_par(4, 0.001, 0); eval_counts; oplot_model_counts(mm,3);oplot_model_counts(mp,5);oplot_model_counts(hm,1);oplot_model_counts(hp,2); fit_counts; list_par; % fast_aped_contin(1) + gauss(1) % idx param tie-to freeze value min max % 1 fast_aped_contin[1].norm 0 1 1 0 2 % 2 gauss[1].area 0 0 8.13071e-05 0 0.01 % 3 gauss[1].center 0 0 8.419847 8.4 8.45 % 4 gauss[1].sigma 0 0 0.005718052 1e-06 0.01 conf(2); conf(3); conf(4); % 0.00756767 0.00378995 % 8.42121 8.41843 % 9.18751e-05 7.08443e-05 plot_data_counts(mm,1);oplot_data_counts(mp,2);oplot_data_counts(hm,3);oplot_data_counts(hp,5); oplot_model_counts(mm,3);oplot_model_counts(mp,5);oplot_model_counts(hm,1);oplot_model_counts(hp,2); window(p2); yrange;xrange(8.2,8.7); rplot_counts(mm); % MEG -1 rplot_counts(mp); % MEG +1 rplot_counts(hm); % HEG -1 rplot_counts(hp); % HEG +1 pid=dup_plot("Contin_fit_9.gif/gif",p2);resize(6*2.54,0.6); yrange;xrange(8.2,8.7); rplot_counts(mm); % MEG -1 rplot_counts(mp); % MEG +1 rplot_counts(hm); % HEG -1 rplot_counts(hp); % HEG +1 close_plot(pid); window(p1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%