Re: Fit a combined spectra

Hi Dave,


I tried to do the fit. I fixed the densities, and constrain the other
parameters  as follows (data_vec2 corresponds to MEG +- 1 data):

 fit_fun( "xaped(1) * phabs(1)" );
 list_par;
 set_par("xaped(1).norm",0.0005,0);
 set_par("xaped(1).temperature",1e+7,0, 1.e6, 1.e8);
 set_par("xaped(1).density",1,1);
 use_thermal_profile;
 set_par( "xaped(1).vturb", 300, 0, 0.0, 1000);
 set_par( "phabs(1).nH", 0.1, 0, 0, 0.5);
 list_par;
 renorm_counts;
 list_par;
 exclude( all_data );
 list_data;
 include(data_vec2);
 list_data;
 group_data( data_vec2, 4 );
 xnotice( data_vec2, lo, lf );
 fit_counts;
 popt.res=4;
 plot_counts( -data_vec2,popt );
 list_par;


But it get the following output for model parameters:

 Parameters[Variable] = 7[4]
            Data bins = 1052
           Chi-square = 2.04803e+13
   Reduced chi-square = 1.954227e+10
 Graphics device/type (? to see list, default /xwin):
Invalid: plot limits
plot failed
xaped(1) * phabs(1)
 idx  param                      tie-to  freeze         value         min 
       max
  1  xaped(1).norm                0     0    -5.087167e-28           0    
      0
  2  xaped(1).temperature     0     0            1e+07     1000000      
1e+08
  3  xaped(1).density             0     1                1           0    
      0
  4  xaped(1).vturb                0     0              300           0   
    1000
  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              0.1           0 
       0.5


When I fix the xaped(1).norm to 0.0005, I get:

 Parameters[Variable] = 7[3]
            Data bins = 1052
           Chi-square = 2.86434e+54
   Reduced chi-square = 2.730544e+51
 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     1           0.0005           0     
     0
  2  xaped(1).temperature    0     0          1000000     1000000       1e+08
  3  xaped(1).density            0     1                1           0     
     0
  4  xaped(1).vturb               0     0                0           0    
   1000
  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              0.5           0   
     0.5  10^22


I the both cases, I get a model with counts=0, for all channels.  I think
I am still missing something.


Thanks

Raul







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.hiddenwith the first line of the message as:
unsubscribe