Re: Model Components in ISIS

On Mon, 19 Jul 2004, Joseph Neilsen  wrote:
>Hello,
>
> My name is Joey Neilsen and I'm at the CfA. I'm doing some
> X-ray pulsar spectral fitting and I would like to be able to
> plot the model components, like setpl add in XSPEC. Is there a
> simple way to do this? Since I am using the pileup kernel, will
> I need to take anything extra into account in plotting? Thanks
> so much for the help,
>
>Cheers,
>
>Joey Neilsen


Hi,

Currently there is no isis intrinsic function which does this,
but it can be done fairly easily using a short function such as
the following:

% Evaluate a subset of current model components
% (or even any function of current model components)
define eval_subfun (dataset, fun)
{
   variable save_fun = get_fit_fun ();
   EXIT_BLOCK
     {
	% in case of error, restore the initial state
	fit_fun (save_fun);
     }
   fit_fun (fun);

   () = eval_counts;

   return get_model_counts (dataset);
}

This function temporarily re-defines the fit-function to be an
expression consisting of any combination of currently defined
model components.  It then evaluates that combination and
restores the original fit-function definition.

I've appended an example which shows how to use the function
defined above.

Using the pileup kernel shouldn't affect this process. It will
be applied for the appropriate datasets independent of the
fit-function specified.

Hope that helps,

-John



% Here's a usage example:

define fake_some_data ()
{
   variable lo, hi, f;
   (lo,hi) = linear_grid (1,20,1000);
   fit_fun ("delta(1)+Powerlaw(1)");
   set_par ("delta(1).norm", 1.e3);
   set_par ("delta(1).lambda", 12.0);
   set_par ("Powerlaw(1).norm", 1.e3);
   f = eval_fun (lo,hi);
   define_counts (lo, hi, f, sqrt(f));
}

fake_some_data();

% Compute the total counts and the
% contribution from each component

variable t, g, p;
() = eval_counts();
t = get_model_counts (1);
g = eval_subfun(1, "delta(1)");
p = eval_subfun(1, "Powerlaw(1)");

% Plot the separate components

ylog;
yrange (1,300);
xlabel ("X");

multiplot ([1,1,1]);
ylabel ("delta+Powerlaw");
hplot(t.bin_lo, t.bin_hi, t.value);
ylabel ("Powerlaw");
hplot(p.bin_lo, p.bin_hi, p.value);
ylabel ("delta");
hplot(g.bin_lo, g.bin_hi, g.value);

% Print the component sums.

vmessage ("\n\npowerlaw+delta = %g", sum(t.value));
vmessage ("         delta = %g", sum(g.value));
vmessage ("powerlaw       = %g", sum(p.value));
----
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