% file: ahfake_model_fit.sl % % based on: ah_ba23_sim.sl in /nfs/spectra/d1/Science/M31/aegis_111128 % % One-time setup lineids and % Generate the fake Astro-H-observed data: #ifnexists MADE_FAKE_DATA % delete all data, arf, rmfs % hence isim will be 1. delete_data(all_data); delete_arf(all_arfs); delete_rmf(all_rmfs); % Use the isis labels .source setup_aped_lineids ; % Load the IACHEC E0102 calibration model: .source isis_cal_model ; % The simulated exposure time in ks: simks = 100.0; % Make a fake observation from the model % load the ASTRO-H response, baseline(7eV) and w/optical filter: % Rmf_OGIP_Compliance=0; ahrspprefix="/nfs/spectra/d1/Science/M31/aegis_111128/"; isim=load_arf(ahrspprefix+"ah_rsp/sxt-s_100208_ts02um_of_intallpxl.arf"); ()=load_rmf(ahrspprefix+"ah_rsp/ah_sxs_7ev_basefilt_20090216.rmf"); assign_rmf (isim,isim); %% assign_arf (isim,isim); % arg: % Error: 16384 RMF bins vs 16384 ARF bins % Largest fractional errors: 1.19202e-07 (low edge) 0.995 (high edge) % Failed: assigning ARF 1 to data set 1 % % Patch the ARF's "0" value: sarf=get_arf(isim); sarf.bin_hi[-1] = Const_hc/1.0e-04; put_arf(isim,sarf); assign_arf (isim,isim); set_data_exposure(isim, simks*1.e3); set_arf_exposure(isim, simks*1.e3); fakeit; set_fake(isim,0); % Astro-H has 1 bin = 1 ev, % with ~ 7eV rmf do a little grouping: group_data(1,2); % Set the data errors: % I generally use this: %% set_fit_statistic ("chisqr;sigma=model"); % But this pays more attention to the lower-counts (continuum) regions: set_fit_statistic ("chisqr;sigma=data"); message(""); list_data; message(""); message(" Now do: .source ahfake_model_fit"); message(""); MADE_FAKE_DATA=1; #endif % Look at the simulation - assumed in data set 1. plot_unit("keV"); xrange(0.35, 2.1); xlog; errorbars(0); ()=eval_counts; list_data; % Plot to a file: %%%PSwin=open_plot("ah_E0102_vpshockOplus.ps/CPS"); %%%set_frame_line_width(3); plot_bin_integral; title("Astro-H: SNR E0102-IACHEC, Observed for "+string(simks)+ +" ks (FWHM = 7eV, Binned to dE = 2 eV)"); %%plt_xmin=0.35; %%plt_xmax=2.1; plt_xmin=0.49+0*0.35; plt_xmax=1.51+0*2.1; xrange(plt_xmin, plt_xmax); xlin; yrange(9.0,1.e5); ylog; set_line_width(7); plot_data_counts(1,14); % Add line ids: set_line_width(3); s_lbl=line_label_default_style() ; s_lbl.label_type = 0; % 0 or 1 s_lbl.char_height = 1.1; s_lbl.justify = 0.0; s_lbl.angle=90.0; s_lbl.top_frac = 0.83; s_lbl.bottom_frac = 0.78; s_lbl.offset = 0.35; ()=evalfile("./plot_simple_lineids"); plot_simple_lineids( Const_hc/plt_xmax, Const_hc/plt_xmin, s_lbl ); % And fit to a vpshock model % - Select the NEI version for, e.g., vpshock model: xspec_xset ("NEIVERS", "2.0"); % Get the ***saved*** model: % "Oplus" means the 6-1 to 10-1 O transitions are included: % (see oplus_model_setup.sl for how these were defined etc) load_par("e0102_vpshock_4elsOplus_sigdat.par"); % Set the range that fitting was/willbe done on xnotice_en(plt_xmin, plt_xmax); ()=eval_counts; % show the model with and w/o the Fe: set_line_width(3); oplot_model_counts(1,3); savefe=get_par("vpshock(1).Fe"); set_par("vpshock(1).Fe", 0.0); ()=eval_counts; oplot_model_counts(1,6); set_par("vpshock(1).Fe", savefe); ()=eval_counts; #iffalse % Close this plot close_plot(PSwin); % Plot it vs wavelength too: PSwin=open_plot("ah_E0102_vpshockOplus_lambda.ps/CPS"); set_frame_line_width(3); % change to A: plot_unit("Angstrom"); xlin; xrange(5.0,23.0); xrange(8.0,20.0); % show the model with and w/o the Fe: set_line_width(7); plot_data_counts(1,14); set_line_width(3); plot_simple_lineids( Const_hc/plt_xmax, Const_hc/plt_xmin, s_lbl ); oplot_model_counts(1,3); savefe=get_par("vpshock(1).Fe"); set_par("vpshock(1).Fe", 0.0); ()=eval_counts; oplot_model_counts(1,6); set_par("vpshock(1).Fe", savefe); ()=eval_counts; close_plot(PSwin); #endif #iffalse % Plot just the Fe component of the model: PSwin=open_plot("ah_E0102_vpshockOplus_FeOnly.ps/CPS"); fit_fun("phabs(1)*vpshock(1)"); set_par([[4:10]],0); eval_counts; title("Fe-only vpshock contribution to fitting the Astro-H SNR E0102-IACHEC simulation"); set_line_width(5); plot_model_counts(1,3); set_line_width(3); plot_simple_lineids( Const_hc/plt_xmax, Const_hc/plt_xmin, s_lbl ); close_plot(PSwin); #endif % - - - - - Steps used to fit the vpshock model and save the par file - - - #iffalse % Fit the model and save it ... set_par("phabs(1).nH", 0.11, 1); set_par("vpshock(1).norm", 0.03, 0); % set kT and tau: set_par("vpshock(1).kT", 0.50, 1); set_par("vpshock(1).Tau_u", 2.0e11, 1); % C,N: set_par(6, 0.0, 1); set_par(7, 0.0, 1); % thaw O, Ne, Mg: set_par([8, 9, 10], 0.3, 0); % freeze Si--Ni to 0.0: set_par([11:16], 0.0, 1); list_par; % Set the range to fit xnotice_en(plt_xmin, plt_xmax); ()=fit_counts; list_par; plot_data_counts(1); oplot_model_counts(1,3); set_par("vpshock(1).Fe", 0.03, 0); ()=fit_counts; list_par; % show the model with and w/o the Fe: plot_data_counts(1,1); oplot_model_counts(1,8); savefe=get_par("vpshock(1).Fe"); set_par("vpshock(1).Fe", 0.0); ()=eval_counts; oplot_model_counts(1,3); set_par("vpshock(1).Fe", savefe); ()=eval_counts; % show the model with and w/o the Fe: set_line_width(7); plot_data_counts(1,14); set_line_width(3); plot_simple_lineids( Const_hc/plt_xmax, Const_hc/plt_xmin, s_lbl ); oplot_model_counts(1,3); savefe=get_par("vpshock(1).Fe"); set_par("vpshock(1).Fe", 0.0); ()=eval_counts; oplot_model_counts(1,6); set_par("vpshock(1).Fe", savefe); ()=eval_counts; % Fit the kT ... and conf to be sure thaw("vpshock(1).kT"); variable clo=0, chi=0; while (clo == chi) { ()=fit_counts; list_par; (clo,chi)=conf(3, 0, 0.1); } clo; chi; % kT clo, chi: 0.7267 -- 0.7276 % Fit the Tau_u too ... and conf to be sure thaw("vpshock(1).Tau_u"); variable clo=0, chi=0; while (clo == chi) { ()=fit_counts; list_par; (clo,chi)=conf("vpshock(1).Tau_u", 0, 0.1); } clo; chi; % Tau clo, chi: very narrow range ~ 2.137e+11 2.1376e+11 % Alternate between kT and Tau fitting... % until both are at minima. % Save the paramters: %% save_par("e0102_vpshock_4els_sigmod.par"); % sigma=model statistic %% save_par("e0102_vpshock_4els_sigdat.par"); % sigma=data statistic % 6/25/12 Added the extra O VII and O VIII lines (gauss()'s above) % and now refit and conf the kT and Tau: % kT clo, chi: 0.7471 -- 0.7487 % Tau clo, chi: very,very narrow range ~ 2.00791e+11 % Save the paramters: %% save_par("e0102_vpshock_4elsOplus_sigdat.par"); % sigma=data statistic #endif