Re: Fit a combined spectra

    Raul> I have combined the data from HEG +/- 1 and MEG +-1. Is it
    Raul> possible to fit the combined spectrum with some xpec model (or
    Raul> what else), using

    Raul> fit_fun("apec");
    Raul> fit_counts;

Hi Raul,

Short answer is "yes".  But there is a lot implied in your email.
Here's what I typically do (or the equivalent of, and there are, of
course, variations on the theme): 

Load your data, arfs, and rmfs:

  h = load_data( "pha2", [3,4,9,10] );  % load HEG and MEG rows

  amm = load_arf( "meg_-1.arf" );
  amp = load_arf( "meg_1.arf" );
  ahm = load_arf( "heg_-1.arf" );
  ahp = load_arf( "heg_1.arf" );

  rmm = load_arf( "meg_-1.rmf" );
  rmp = load_arf( "meg_1.rmf" );
  rhm = load_arf( "heg_-1.rmf" );
  rhp = load_arf( "heg_1.rmf" );

If you use the tgcat.sl suite and files donwloaded from tgcat.mit.edu,
all that can be done via:

  d = load_set_acis( "your_data_directory", [3,4,9,10] );

(it assumes the tgcat generic file naming scheme)

(I like to keep the returned values, especially for use in scripts.  For
interactive use, you don't necessarily need to use "x = ...".)

Set flags for combining the data:

  gh = combine_datasets( h[[0,1]] );
  gm = combine_datasets( h[[2,3]] );

(without the variable h, you could do:  gh = combine_datasets( [1,2] ); etc)

Note that this does not add any spectra, arfs or rmfs, but solely sets
flags identifying which are to be combined.

You can use the xspec apec model if you want.  If you have AtomDB
installed, then there are advantages to using a native ISIS plasma
model, since then you can query the database after model evaluation for
the explicit line contributions (e.g., line ids, fluxes, etc, per
wavelength range, by element, ion; that's another topic; see
create_aped_fun, plot_group, for examples).  

For the question at hand, we'll use the xspec model:

  fit_fun( "apec(1) * phabs(1)" );

  list_par;

  set_par( ...);

  exclude( all_data ); % in case I have multiple datasets loaded
 
  include( h );        % include the spectra of interest

  group_data( h, 2 );  % bin up a bit, uniformly (more or less if needed)

  xnotice( h, 1.7, 25 ); % notice the data and region of interest

  fit_counts;

Fitting is the same whether you have done combine_datasets, or not.  In
this case, you have specified that MEG +-1 be combined, and that HEG +-1
be combined.  When fit, the model is folded through the response for
each grating and order, then before computing the statistic, the MEG
data are added, the MEG model counts are added, and ditto for HEG, then
the statistic computed (see the help for combine_datasets).  So MEG and
HEG are fit jointly, while each one is combined.  If you are photon
starved and want to combine them all, you first need to
match_dataset_grids:

  match_dataset_grids( h[[2,0,1]] ); % put HEG on MEG grid

  g = combine_datasets( h );


Now, visualizing is another issue.  If you use the tgcat.sl suite, it
will load fancy_plots.sl, which has functions for plotting combined
data:

  popt.res = 4; % combined residuals
  plot_counts( -[2,3], popt );  % negative sign means combine before plotting

or equivalently:

  plot_counts( -gm, popt );  % using group id.


There are also  plot_data() (count rate) and plot_unfold() (fluxed
units).




-- Dave

David Huenemoerder  617-253-4283 (o); -253-8084 (f); http://space.mit.edu/home/dph
MIT Kavli Institute for Astrophysics and Space Research
One Hampshire Street, Room NE80-6065
Cambridge, MA  02139
[Admin. Asst.: Elaine Tirrell, 617-253-7480, egt_at_email.domain.hidden
----
You received this message because you are
subscribed to the isis-users list.
To unsubscribe, send a message to
isis-users-request_at_email.domain.hiddenwith the first line of the message as:
unsubscribe