
% Load Cycle 7 ACIS-S3 (no gratings) arf and rmf obtained from
% http://cxc.harvard.edu/caldb/prop_plan/imaging/index.html

() = load_arf("aciss_aimpt_cy07.arf");
() = load_rmf("aciss_aimpt_cy07.rmf");

% Create an empty data set by assigning the arf and rmf

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

% Use a simple, absorbed power-law model

require("xspec");
fit_fun("phabs(1)*powerlaw(1)");

% Set the model parameters

set_par(1,0.01);  % Neutral hydrogen column of 10^21
set_par(2,1);     % Temporary power-law normalization
set_par(3,2);     % Power-law photon index

% Define a simple function to determine the flux in a given keV band

define kev_flux (id, kev_lo, kev_hi)
{
   % convert from wavelength grid [A] to energy grid [keV]
   variable m = _A(get_model_flux(id));
   
   % convert photons/sec => ergs/sec
   m.value *= 0.5 * (m.bin_hi + m.bin_lo) * 1.602e-9;
   
   % return integral over specified band [erg/sec]
   return rebin (kev_lo, kev_hi, m.bin_lo, m.bin_hi, m.value)[0];
}

% Evaluate the model (ignore bad chi^2 - we haven't faked the data yet)

() = eval_counts;

% Get the flux in the 0.5 to 8 keV band

variable flux = kev_flux(1,0.5,8.);

% Renormalize the model to get 10^-12 erg cm^-2 s^-1

set_par(2,1.e-12/flux);

% Set an exposure time.  Here, we choose 50 ksec

set_arf_exposure(1,5.e4);

% Fake the data without pileup

fakeit;

% Mark this data set as "real" so it won't be overwritten
% by the next 'fakeit'

set_fake(1,0);

% Now create a second set of fake data with pileup

assign_arf(1,2);   % assign same arf and rmf to new data
assign_rmf(1,2);

set_kernel(2,"pileup");  % assign pileup kernel to new data

set_frame_time(2,3.214); % assign the nominal frame time to new data

fakeit;                  % Create fake data

% Analyze the data as per normal.  You have two data sets from the same
% source distribution, except one has the effects of pileup applied.

