% file: load_tycho.sl

%
% Load the ACIS-I data from obsid 8551
%

% delete all previous data etc
delete_data(all_data);
delete_arf(all_arfs);
delete_rmf(all_rmfs);


dat_dir="../Data/obs_8551";

()=load_data(dat_dir+"/extract.pi");
()=load_arf(dat_dir+"/extract.warf");
()=load_rmf(dat_dir+"/extract.wrmf");

assign_arf(1,1);
assign_rmf(1,1);

% Use the background extraction ?
%% define_back(1, NULL);
()=define_back(1, dat_dir+"/extract_bkg.pi");

% Do some grouping, noticing, etc.
%
% Include an systematic error level --
% there are lots of bins so 10% per bin should still constrain any
% model-fitting pretty well.
set_sys_err_frac(1, 0.10);
%    or not:
%%set_sys_err_frac(1, NULL);

% Set the binning to high S/N so that the high-energy range
% has reasonable sized bins
group(1; min_sn=10.0);

pltxmin=0.5;
pltxmax=8.0;

xnotice_en(1,pltxmin,pltxmax);

list_data;



%
% Plot the data
%

% To a .ps file?
%% PSwin = open_plot("tycho_with_hydro.ps/CPS");

plot_bin_density;

charsize(1.1);
title("Tycho on ACIS-I:  Data (black), Background (orange) (obsid 8551, " +
   string(get_data_exposure(1)/1.e3)+" ks)");
xrange(pltxmin,pltxmax);
xlog;
%%yrange();  % auto-scale
% set a range, log
yrange(0.03,200.0);
ylog;
% set a range, lin:
%%ylin;
%%yrange(0.0,130.0);

% Show the data with gray error bars
errorbars(1);
plot_data_counts(1,15);
errorbars(0);
oplot_data_counts(1,1);

% Show background level
fit_fun("0.0*phabs(1)");
()=eval_counts;
oplot_model_counts(1,8);


% Can comment this out
%% #iffalse

%
% load the vhx model
%
()=evalfile("vhx_model_setup");
xspec_xset("NEIVERS", "1.1");

tid=720;
load_par("/nfs/spectra/d1/Science/VikXray/vhx_130101/Tycho_out/model_" +
  +string(1000+tid)+".par");

()=eval_counts;

save_par("temp.par");
oplot_model_counts(1,2);

color(2);
xylabel(2.8, 80.0, "Hydro-shells model, tid="+string(tid));
color(1);

% Show components:
% It's mostly CSM, so show the ejecta contibution:
set_par("csm(1).norm", 0.0);
()=eval_counts;
oplot_model_counts(1,3);
color(3);
xylabel(3.3, 0.7*80.0, "Ejecta-component only");
color(1);
% put it back
set_par("csm(1).norm", 1.0);

%% #endif


%% close_plot(PSwin);