Re: Fit a combined spectra

Hi Dave,

Thank you for your help. I did the same thing (almost): I download from
tgcat.mit.edu the following files (obsid 3753):

heg_-1.pha  heg_1.pha  meg_-1.pha  meg_1.pha

heg_-1.rmf  heg_1.rmf  meg_-1.rmf  meg_1.rmf

heg_-1.arf  heg_1.arf  meg_-1.arf  meg_1.arf

then I did a list of pha files for the input of the following function
(g_fluxv2).

delete_data(all_data);
delete_rmf(all_rmfs);
delete_arf(all_arfs);
delete;
require("tgcat");

alias("plot_data_counts","pdc");
alias("oplot_model_counts","opmc");
alias ("rplot_counts", "rp");
alias ("oplot_data_counts", "opdc");
alias ("plot_data_flux","pdfl");
%
%
public define g_fluxv2 (file,lo,lf)
{
  variable fp, lines;
  variable i;

  fp = fopen (file,"r");
  if (fp==NULL)
   return -1;

  lines = fgetslines(fp);
  if (lines == NULL)
   return -1;

  variable nl = length(lines);
  print(nl);
  print(lines);
  variable index=[1:nl];
  variable rootnames = _at_email.domain.hidden  variable phafiles = _at_email.domain.hidden  variable rmffiles = _at_email.domain.hidden  variable arffiles = _at_email.domain.hidden  for (i = 0; i <= nl-1; i++)
  {
    rootnames[i]= substr(lines[i],1,strlen(lines[i])-4);
%    rootnames[i]= substr(lines[i],1,strlen(lines[i])-14);
    phafiles[i] = rootnames[i] + "pha";
    rmffiles[i] = rootnames[i] + "rmf";
    arffiles[i] = rootnames[i] + "arf";
    print(phafiles[i]);
    load_dataset(phafiles[i],rmffiles[i],arffiles[i]);
    index[i]=i+1;
  }
%
% Rebin and combine the HEG and MEG Chandra Data
% 1st order
%

 variable data_vec1 = [index[0],index[nl/4]];  % taking HEG +1 -1
 variable data_vec2 = [index[2*nl/4],index[3*nl/4]]; % taking MEG +1-1
notice(all_data);
 group_data(all_data,1);
 match_dataset_grids(data_vec1);
 match_dataset_grids(data_vec2);
 variable wh = [1.0,1.0];
 variable gr1 = combine_datasets(data_vec1,wh);
 variable gr2 = combine_datasets(data_vec2,wh);
% group_data( data_vec1, 4 );
% group_data( data_vec2, 4 );
 flux_corr(data_vec1);
 flux_corr(data_vec2);
%
% FIT!
%
 variable t = default_plasma_state ();
 create_aped_fun ("xaped", t);
 fit_fun( "xaped(1) * phabs(1)" );
 list_par;
 set_par("xaped(1).norm",1,0);
 set_par("xaped(1).temperature",1e+7,0);
 set_par("xaped(1).density",1,0);
 list_par;
 exclude( all_data );
 include(data_vec1);
 group_data( data_vec1, 4 );
 xnotice( data_vec1, 1.7, 25 );
 fit_counts;
 plot_counts( -gr1 );
 list_par;
 return;




Everey goes fine until the FIT. It seems that fit is fine but whet I call
for the model parameters I got this:

 Parameters[Variable] = 7[4]
            Data bins = 1978
           Chi-square = 5.581111e+18
   Reduced chi-square = 2.82731e+15
 Graphics device/type (? to see list, default /xwin):
xaped(1) * phabs(1)
 idx  param                       tie-to  freeze         value         min
        max
  1  xaped(1).norm               0     0    -5.304679e-25      0           0
  2  xaped(1).temperature    0     0            1e+07           0           0
  3  xaped(1).density            0     0                1               0 
         0
  4  xaped(1).vturb               0     1                0               0
          0
  5  xaped(1).redshift            0     1                0               0
          0
  6  xaped(1).metal_abund    0     1                1               0     
     0
  7  phabs(1).nH                   0     0                1               
0      100000
0


I don know, What am I doing wrong. But I cannot get a good fit.

Thank you for your help.

Raul






 return;On Wed, May 22, 2013 4:50 pm, David P. Huenemoerder wrote:
>
>     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.hidden> with the first line of the message as:
> unsubscribe
>
>



----
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