Re: Fit a combined spectra

From: Raul Eduardo Puebla <rep54_at_email.domain.hidden>
Date: Fri, 24 May 2013 15:43:27 -0400
lo = 2.0
lf = 27.0

sorry

On Fri, May 24, 2013 3:40 pm, Raul Eduardo Puebla wrote:
> 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
Received on Fri May 24 2013 - 15:43:36 EDT

This archive was generated by hypermail 2.2.0 : Thu Jun 06 2013 - 12:24:50 EDT