Marx 5.0.0
HETG-MARX Simulation - Powerlaw w/lines Example
Marx Home  |
| Documentation |
| Isis |
| Examples |
| Basic Setup |
| ARFs and RMFs |
| Imaging |
| Powerlaw |
| Pileup |
| LETG/ACIS |
| Powerlaw |
| HETG |
| Thermal |
Powerlaw  |
| CXC@MIT |
| Chandra |
The purpose of this example, contributed by Dan Dewey, is to show how to use marx to simulate an HETG observation of a point (and simply-extended) source with a user-specified spectrum. User customization for other HETG - marx applications should be straightforward.
After some Preliminaries, the basic steps described in detail below are:
- Create the spectal file - Use isis and standalone
marxfluxto make a spectral input file for marx. - Setup and run marx - Define and run the simulation to get a FITS event file; it can be viewed, e.g., with ds9.
- Extract HETG spectra - Use standard CIAO tools, called conveniently from isis scripts, to create pha2 and response files.
- Perform spectral analysis - Carry out model fitting in isis, sherpa, or xspec using as input the pha2 and response files.
Each of these steps is described in a text file that includes comments and the command lines which can be cut-and-pasted; these files are displayed on this page as well. The user is encouraged to download and use these files as starting points for customization and recording ones own analysis steps. Some screen shots from this example are given here as enticement:
- Plotting the isis-created spectral model. The spectra are shown on the PGPLOT window appearing above the isis command window (lower left) and editor window (xjed, lower right) where the isis_spectrum.txt file is being displayed and used as cut-ans-paste source.
- Looking at marx-created FITS file with ds9. The is from the marxcat'ed union of the point source simulation and the disk simulation.
- Dispersed spectral images from disk simulation. Looking at the Plt_imsp.ps file created as one of the summary plots in the isis-scripted ciao processing.
- The four first-order spectra dispayed in ISIS. Note that the line regions are indicated (gray) and will be excluded in continuum fitting.
- Continuum fitting of the spectra. The four counts spectra are jointly fit with the continuum model; residuals in chi are shown as well.
- Fitting Fe line to HEG spectra. A Gaussian line is jointly fit to the two HEG spectra and conficence limits on its width are reported. This shows the output of running the isis_analysis.txt script.
Preliminaries
The steps below are performed in a linux-like environment
with installed versions of
marx, isis, and CIAO in the search path(*). Unless indicated
otherwise, commands are carried out
in a ``working directory'' that the user has
write privilege in, e.g., /nfs/mycomputer/d1/AGN_amazing/hetg_simulations/ .
Add the following to the ``working directory'':
- marxflux - The marxflux file (click it to download). Then do a
chmod a+x marxfluxto ensure that it can be executed. - TG-Cat reprocessing scripts - Download the file
TGCat_scripts.tar to the ``working directory'' and do
a ``
tar -xvf TGCat_scripts.tar'' which will create and fill a ``scripts'' sub-directory. Go into the scripts directory (cd scripts/) and do a ``cat badpix_NULL.sl >> tg_repro_fun.sl'' (This step appends a redefinition of a procedure so that no badpix files are used in the marx simulation processing.) Back in the working directory (cd ..) the TGCat_scripts.tar file can be removed (rm TGCat_scripts.tar).
Thanks to Dave Huenemoerder for making these ``TGCat scripts'' available ahead of their release. - marx_orig.par - Copy fresh versions of the marx parameter files
into the working directory; doing
marx --helpwill show the location of the marx parameter files. Also do a
cp marx.par marx_orig.par
to save the original marx.par file to use as a starting point in simulation scripts.
Now, multiple simulations and analyses can be carried out in this ``working directory'' without redoing these steps above. FYI, the directory now contains:
[unix] ls marxasp.par marxflux marx_orig.par marx.par marxpileup.par scripts
isis may have conflicts. In these cases it is useful to only ``setup'' CIAO when needed
(the third step) and then restart a new window when done with CIAO processing.
Create the spectral file
The following instructions are in isis_spectrum.txt.
file: isis_spectrum.txt
Use MARX to simulate an HETG observation of a powerlaw w/lines.
# - - - Create the spectral file
This is similar to the ACIS example except that the
source is brighter and has a less steep powerlaw.
Create the .par file using isis and XSPEC models:
Use a power law and include two "flourescence" lines (Fe, Ne) for fun:
These steps are carried out in the "working directory".
Start isis:
[unix] isis
Define the model:
isis> fit_fun("phabs(1)*powerlaw(1)+gauss(1)+gauss(2)");
isis> list_par; % This will show the model parameters.
Set the model parameters:
isis> set_par(1, 0.80);
isis> set_par(2, 0.020);
isis> set_par(3, 1.20);
Set the Gaussian line parameters.
It's also possible to set the parameters by name instead of their number.
Note that the "center" must be specified in Angstroms,
so keV values are converted using Const_hc [ keV A ].
Fe line and Ne line centers
isis> set_par("gauss(1).center", Const_hc/6.4038);
isis> set_par("gauss(2).center", Const_hc/0.8486);
Set the areas and widths of the lines:
isis> set_par("gauss(1).area", 3.0e-4);
isis> set_par("gauss(2).area", 1.0e-4);
Make the lines have small but non-zero widths
of 0.020 A FWHM, or sigma = 0.020/2.35 = 0.0085 A.
These correspond to velocity widths of:
c * 0.020/1.936 = ~ 3100 km/s FWHM for Fe line
c * 0.020/14.61 = ~ 410 km/s FWHM for Ne line
isis> set_par("gauss(1).sigma", 0.0085);
isis> set_par("gauss(2).sigma", 0.0085);
Show the parameters
isis> list_par;
phabs(1)*powerlaw(1)+gauss(1)+gauss(2)
idx param tie-to freeze value min max
1 phabs(1).nH 0 0 0.8 0 100000 10^22
2 powerlaw(1).norm 0 0 0.02 0 1e+10
3 powerlaw(1).PhoIndex 0 0 1.2 -2 9
4 gauss(1).area 0 0 0.0003 0 0 photons/s/cm^2
5 gauss(1).center 0 0 1.936103 0 0 A
6 gauss(1).sigma 0 0 0.0085 0 1 A
7 gauss(2).area 0 0 0.0001 0 0 photons/s/cm^2
8 gauss(2).center 0 0 14.61044 0 0 A
9 gauss(2).sigma 0 0 0.0085 0 1 A
Save the model/parameters to a file and exit isis.
isis> save_par ("plaw_hetg.par");
isis> exit;
Now , use marxflux to evaluate the .par file to create a MARX spectrum file.
Here a fine binning is used having bins with dE/E = 0.0003 (v_bin = 90 km/s).
This gives a high resolution spectrum across the whole 1 to 40 A
(0.31 to 12.4 keV) range suitable for use with the HETG as well as, e.g.,
future microcalorimeter instruments.
[unix] ./marxflux -e '_A(exp([log(1.0):log(40.0):0.0003]))' \
plaw_hetg.par plawflux_hetg.tbl
The plawflux_hetg.tbl file can now be used by marx to define the spectrum.
...................
FYI: compare the flux of this model with the "plaw.p" of the ACIS simulation.
Use isis to read in and evaluate the models:
Define a wavelength grid to evaluate the model (fit function) on,
from 8 keV to 0.5 keV, in 2000 bins:
isis> (lo, hi) = linear_grid (Const_hc/8.0, Const_hc/0.5, 2000);
The model used for the ACIS example:
isis> load_par("plaw.p");
isis> facis = eval_fun(lo,hi);
The model for this HETG simulation:
isis> load_par("plaw_hetg.par");
isis> fhetg = eval_fun(lo,hi);
The total fluxes:
isis> sum(facis);
0.00059783 % photons/sec/cm^2
isis> sum(fhetg);
0.0253021
The ACIS example has ~ 6.0e-4 ph/s/cm^2, whereas this HETG example
has 2.5e-2 ph/s/cm^2 or a factor of ~ 42 times the photon rate over the
0.5 to 8 keV range.
Can also plot the spectra:
isis> yrange(); ylog;
isis> xrange(); xlin;
isis> label("Wavelength (A)","Flux (ph/s/cm^2/bin)","Simulation spectra");
isis> hplot(lo, hi, fhetg);
isis> ohplot(lo, hi, facis,2);
Or versus log Energy:
isis> yrange(); ylog;
isis> xrange(); xlog;
isis> label("Energy (keV)","Flux (ph/s/cm^2/bin)","Simulation spectra");
isis> hplot(Const_hc/hi, Const_hc/lo, fhetg);
isis> ohplot(Const_hc/hi, Const_hc/lo, facis);
.....................
Setup and run marx
The following instructions are in marx_hetg_plaw.txt.
file: marx_hetg_plaw.txt
Use MARX to simulate an HETG observation of a powerlaw w/lines.
# - - - Setup and Run MARX
Roughly these steps are the same as for the ACIS example, except that
the HETG is inserted by using: GratingType="HETG".
It can be convenient to use pset to set the parameters in a
local par file, so that method is demonstrated and used here
(pset is provided with ciao or lheasoft.)
Set up the various marx parameters as desired. Note that
the user can edit and re-do the following lines for new
or modified simulations. As an example, the first simulation
here is for a POINT source; it is followed by a few additional
lines to do a second simulation with a DISK source.
In the working directory paste these sets of lines to the unix prompt:
- - - Simulation 1 - Point Source
# Make a par file for my marx simulation:
cp marx_orig.par mysim.par
# set the spectrum file to use:
pset mysim.par SpectrumType="FILE"
pset mysim.par SpectrumFile="plawflux_hetg.tbl"
pset mysim.par SourceFlux=-1
# Set other parameters of the simulation:
# Using 50 ks
pset mysim.par ExposureTime=50000
pset mysim.par OutputDir="hetg_plaw"
pset mysim.par DitherModel="INTERNAL"
# Use the HETG with ACIS-S:
pset mysim.par GratingType="HETG"
pset mysim.par DetectorType="ACIS-S"
# Some other parameters it can be useful to set:
# Date of observation (effects ACIS QE)
pset mysim.par TStart=2009.50
# Roll of the observation: 0 puts average dispersion along E -- W.
pset mysim.par Roll_Nom=0.0
# Source RA/DEC (degrees)
pset mysim.par SourceRA=250.000
pset mysim.par SourceDEC=-54.000
# Pointing RA/DEC (degrees)
pset mysim.par RA_Nom=250.000
pset mysim.par Dec_Nom=-54.000
# Finally,
# Run the simulation:
marx @@mysim.par
... marx runs, ending with:
Writing output to directory 'hetg_plaw' ...
Total photons: 3495031, Total Photons detected: 255172, (efficiency: 0.073010)
(efficiency this iteration 0.073986) Total time: 50000.000492
# Create the fits event file and the aspect solution file:
marx2fits hetg_plaw hetg_plaw_evt2.fits
marxasp MarxDir="hetg_plaw" OutputFile="hetg_plaw_asol1.fits"
#
# all done with MARX simulation 1.
- - -
- - - Simulation 2 - Disk Source
# Do another simulation keeping most things the same as above
# by starting with the "mysim.par" as it is left from above
# but changing a few things to use a DISK source:
# New output dir for the Disk
pset mysim.par OutputDir="hetg_pldisk"
# Change the SourceType:
pset mysim.par SourceType="DISK"
# a thin disk with average radius ~ 2.0"
pset mysim.par S-DiskTheta0=1.7
pset mysim.par S-DiskTheta1=2.3
# Run the second simulation:
marx @@mysim.par
# Create the fits event file and the aspect solution file:
marx2fits hetg_pldisk hetg_pldisk_evt2.fits
marxasp MarxDir="hetg_pldisk" OutputFile="hetg_pldisk_asol1.fits"
#
# all done with MARX simulation 2.
- - -
We can look at the simulation output event files with ds9 to check that they
are as expected before continuuing with ciao-processing, e.g.,
[unix] ds9 hetg_pldisk_evt2.fits
Use a "log" scale to see the dispersed arms and streak as well as the
brighter zeroth-order image. We can also use ds9 to record the center of the
disk (simulation 2) in X,Y coord.s: 4096.5, 4096.5, for further processing.
........
FYI: the routine marxcat allows simulations to be combined, e.g., we
could do the following to make a combination of the point and disk events:
[unix] marxcat hetg_plaw hetg_pldisk hetg_plboth
and then create fits and asol files as above.
This allows more complex spatial-spectral simulations to be done with marx.
-------
Extract HETG spectra
The following instructions are in ciao_marx_process.txt.
file: ciao_marx_process.txt
Use MARX to simulate an HETG observation of a powerlaw w/lines.
# - - - Extract HETG spectra
Here we'll use TGCat scripts to run the ciao tools to process the simulation.
Note that for real Chandra data when using the TGCat scripts, we'd first do
Download and Configure steps which would create a sub-directory, e.g., obs_8487,
with these files and/or simlinks:
evt0 -> secondary/acisf08487_000N002_evt1.fits
asol_1 -> primary/pcadf291118212N002_asol1.fits
asol.list
bpix1 -> primary/acisf08487_000N002_bpix1.fits
flt1 -> secondary/acisf08487_000N002_flt1.fits
msk1 -> secondary/acisf08487_000N002_msk1.fits
pbk0 -> secondary/acisf291118240N001_pbk0.fits
stat1 -> secondary/acisf08487_000N002_stat1.fits
For MARX simulation processing we will create a similar sub-directory
but with only the first three items in it (evt0, asol_1, and asol.list).
- - -
In the "working directory", setup the analysis sub-directory and files
for each simulation:
...for the point source simulation:
[unix] mkdir obs_hplaw
[unix] cp hetg_plaw_evt2.fits obs_hplaw/evt0
[unix] cp hetg_plaw_asol1.fits obs_hplaw/asol_1
[unix] echo asol_1 | cat > obs_hplaw/asol.list
...for the disk simulation:
[unix] mkdir obs_hpldisk
[unix] cp hetg_pldisk_evt2.fits obs_hpldisk/evt0
[unix] cp hetg_pldisk_asol1.fits obs_hpldisk/asol_1
[unix] echo asol_1 | cat > obs_hpldisk/asol.list
Now, process each simulation in turn
by going into the simulation sub-directory:
[unix] cd obs_hplaw ( OR: [unix] cd obs_hpldisk )
and running isis
[unix] isis
isis>
% Once in isis carry out the following steps to run the
% various ciao tools from isis.
% load the s/w we're using
variable tgc_path = "../scripts" ;
prepend_to_isis_load_path( tgc_path );
require("tg_repro_fun");
% make the obs par
set_exec( 2 );
setup_pfiles; % makes a param directory
obs_info = read_config( "evt0" );
make_obspar;
% Do nominal ACIS processing, but witout any bad pixels.
%
% Use the "new parameters" option to change the nominal behavior
% of tools. In this way, here and below, we can change parameter
% values for the tools to customize the processing.
% Each time we redefined nps (" = Assoc_Type[ Any_Type ] ;" ) to clear it.
%
nps = Assoc_Type[ Any_Type ] ;
nps["badpixfile"]="NONE";
nps["apply_cti"]="no";
acis_process_events_te(nps); % this takes of order 10 seconds
% Do event filtering - just for grade, status and energy
% ( For real data we would include GTI and destreak as well:
% acis_evt_filter_ds )
% Change the output file to evt1:
nps = Assoc_Type[ Any_Type ] ;
nps["outfile"]="evt1";
acis_evt_filter_1(nps);
% Get zeroth-order location:
tgdetect;
%
% It's OK to get this message (CXC DS people say) :
% # DMCOPY (CIAO 4.0 Beta 2): Bad data type in filter string formatting
% Read the source location and counts in and show them:
(x,y,c) = read_src1a_pos;
x ;
y ;
c ;
% For the POINT source example (obs_hplaw), tgdetect gives an OK center
% location so we can leave x,y as is and skip to "Create the mask..."
% OR:
%
% For the DISK simulation, obs_hpldisk, tgdetect finds a point on
% the disk, not the real center.
% We can use ds9 to manually find the real values;
% if we didn't do it already we can do it from isis now:
% isis> !ds9 evt1
% and then put the center values in by hand here:
x = 4096.5;
y = 4096.5;
% Create the mask and resolve events:
tg_create_mask(x,y);
tg_resolve_events_te; % this takes ~ 10 seconds
% Make the PHA file
%
% Use nominal extraction for the point source:
tgextract;
% OR:
%
% Extract a little wider to accomodate the disk source:
nps = Assoc_Type[ Any_Type ] ;
nps["min_tg_d"]=-0.0012;
nps["max_tg_d"]=0.0012;
tgextract(nps);
% Make a light curve of the dispersed events
lightcurve_ha;
% And create the response files:
orders=[-1,1];
% Modify the nominal parameters for the garf for marx simulation:
nps = Assoc_Type[ Any_Type ] ;
nps["maskfile"]="NONE";
nps["pbkfile"]="NONE";
nps["dafile"]="NONE";
nps["osipfile"]="NONE";
%
% Make the arfs and rmfs:
make_responses(orders, NULL, nps); % this takes a few minutes.
% et voila.
% Now we have pha2 file and the first-order responses.
% Make summary plots... (not required, but they are nice)
require( "summarize" ) ;
% This creates a bunch of Plt_*.ps files.
% for the point source simulation:
variable ddir = "../obs_hplaw" ;
hetgs_te_summary( ddir, "$ddir/Plt"$ );
% OR:
% for the disk simulation:
variable ddir = "../obs_hpldisk" ;
hetgs_te_summary( ddir, "$ddir/Plt"$ ); % takes ~ 10 seconds
% All done processing in this sub-directory
isis> exit;
Perform spectral analysis
The following instructions are in isis_analysis.txt.
% file: isis_analysis.txt
%
% Use MARX to simulate an HETG observation of a powerlaw w/lines.
%
% - - - Perform spectral analysis
% These steps are carried out in ISIS, running in the "working directory".
%
% [unix] isis
% isis>
%
% These steps can be carried out by cutting and pasting from this file into ISIS;
% or they can all be carried out by source'ing this file:
% isis> .source isis_analysis.txt
%
% - - - - -
% Set the simulation sub-directory:
% point source:
variable obs_dir = "obs_hplaw";
% disk source:
%% variable obs_dir = "obs_hpldisk";
% Load the pha's and responses:
()=load_data(obs_dir+"/pha2",[9,10,3,4]);
()=load_arf(obs_dir+"/meg_-1.arf");
()=load_arf(obs_dir+"/meg_1.arf");
()=load_arf(obs_dir+"/heg_-1.arf");
()=load_arf(obs_dir+"/heg_1.arf");
assign_arf(1,1);
assign_arf(2,2);
assign_arf(3,3);
assign_arf(4,4);
()=load_rmf(obs_dir+"/meg_-1.rmf");
()=load_rmf(obs_dir+"/meg_1.rmf");
()=load_rmf(obs_dir+"/heg_-1.rmf");
()=load_rmf(obs_dir+"/heg_1.rmf");
assign_rmf(1,1);
assign_rmf(2,2);
assign_rmf(3,3);
assign_rmf(4,4);
% - - - - -
% Open a plotting window...
variable onepane = open_plot("/XSERVE",1,1);
% Group and notice the data:
group_data(all_data, 8);
% 0.6 to 7.5 keV:
% (Looks like HEG is cutoff before 8 keV...)
xnotice(all_data, Const_hc/7.5, Const_hc/0.6);
% List the data
list_data;
% Use keV for the X axis in plots:
plot_unit("keV");
% Use log x axis, linear y axis starting at 0.
xlog; xrange(0.5,10.0);
ylin; yrange(0.,);
% Can look at the data
title("MARX simulation of HETG observation");
plot_data_counts(1);
% show all of these in grey
oplot_data_counts(1,14);
oplot_data_counts(2,14);
oplot_data_counts(3,14);
oplot_data_counts(4,14);
% Ignore the line regions for the continuum fitting we'll do
% Fe region
ignore(all_data, 0.95*Const_hc/6.4, 1.05*Const_hc/6.4);
% Ne region
ignore(all_data, 0.97*Const_hc/0.8486, 1.03*Const_hc/0.8486);
%
% Replot the spectra again - use some colors:
oplot_data_counts(1,2); % MEG -1 red
oplot_data_counts(2,8); % MEG -1 orange
oplot_data_counts(3,4); % HEG -1 blue
oplot_data_counts(4,3); % HEG -1 green
% - - - - -
% Define the fitting function:
fit_fun("phabs(1)*powerlaw(1)");
% Set initial values and limits for the parameters;
% keep them all "thawed" with 0 as third argument.
set_par (1, 1.00, 0, 0.1, 10.0);
set_par (2, 0.01, 0, 1.e-4, 1.0);
set_par (3, 1.00, 0, 0.5, 2.5);
%
% use "set_fit_statistic" and "set_fit_method" here
% to define fitting specifics; help is available via:
% isis> help set_fit_method
%
% Use Gehrels to help with low-counts bins
set_fit_statistic ("chisqr;sigma=gehrels");
% Could use two fitting methods in sequence
% if needed to get the best local minimum
% Comment out the first fitting for this example.
%% message(" - - - - -");
%% message("Do a minim fit ...");
%% set_fit_method ("minim;default");
%% ()=fit_counts;
%% list_par;
message(" - - - - -");
message("Do a lmdif fit ...");
set_fit_method ("lmdif;default");
()=fit_counts;
list_par;
% For obs_hplaw :
% list_par gave these results:
%
% phabs(1)*powerlaw(1)
% idx param tie-to freeze value min max
% 1 phabs(1).nH 0 0 0.7965323 0.1 10 10^22
% 2 powerlaw(1).norm 0 0 0.01999332 0.0001 1
% 3 powerlaw(1).PhoIndex 0 0 1.179763 0.5 2.5
message(" - - - - -");
% Show the data and residuals for the 4 data sets:
variable fourpane = open_plot("/XSERVE",2,2);
xlog; xrange(0.5,10.0);
ylin; yrange(0.,);
% don't show the error bars on main histogram
window(fourpane);
errorbars(10000);
rplot_counts(1);
rplot_counts(2);
rplot_counts(3);
rplot_counts(4);
% - - - - -
% Get parameter conf ranges
% For help on conf do: isis> help conf
% We could just do:
% isis> conf(1,2); etc.
% But the next lines make prettier output and
% show some more S-Lang stuff that can be useful...
variable ipar=0, clo, chi;
_for ipar (1,3,1)
{
message("Finding conf for parameter "+string(ipar)+" ...");
(clo, chi)=conf(ipar, 2);
message(" 99\% range is: "+string(clo)+" to "+string(chi));
};
% For obs_hplaw :
% >>--> output from the confidence ranges are
% in good agreement with the known input parameter
% values of 0.80, 0.020, and 1.20 :
%
% Finding conf for parameter 1 ...
% 99% range is: 0.754603 to 0.839679
% Finding conf for parameter 2 ...
% 99% range is: 0.0189322 to 0.0211297
% Finding conf for parameter 3 ...
% 99% range is: 1.13933 to 1.22082
message(" - - - - -");
% - - - - -
% Now look closely at the narrow-line regions...
% Use a finer binning
notice(all_data);
group_data(all_data, 2);
()=eval_counts;
window(fourpane);
xrange(0.799,0.901); xlin; rplot_counts(1);
xrange(5.99,7.01); xlin; rplot_counts(3);
xrange(0.799,0.901); xlin; rplot_counts(2);
xrange(5.99,7.01); xlin; rplot_counts(4);
% - - - - -
% Fit the width of the Fe line
% Plot the HEG and continuum model:
window(fourpane);
xrange(5.99,7.01); xlin; rplot_counts(3);
xrange(5.99,7.01); xlin; rplot_counts(4);
% Add a Gaussian to the fit function:
fit_fun("phabs(1)*powerlaw(1)+gauss(1)");
% set the gaussian parameters:
set_par("gauss(1).area", 0.001, 0, 1.e-8, 1.0);
set_par("gauss(1).center", Const_hc/6.40, 0, Const_hc/6.7, Const_hc/6.2);
set_par("gauss(1).sigma", 0.001, 0, 0.0001, 0.1);
% Notice just the HEG spectra in the range:
xnotice(3, Const_hc/7.01,Const_hc/5.99);
xnotice(4, Const_hc/7.01,Const_hc/5.99);
% ignore the MEG spectra:
ignore([1,2]);
% Freeze the continuum nH and PhoIndex, leaving norm free:
freeze([1,3]);
()=fit_counts;
list_par;
% For obs_hplaw :
% list_par gave these results:
%
% phabs(1)*powerlaw(1)+gauss(1)
% idx param tie-to freeze value min max
% 1 phabs(1).nH 0 1 0.7965323 0.1 10 10^22
% 2 powerlaw(1).norm 0 0 0.01900418 0.0001 1
% 3 powerlaw(1).PhoIndex 0 1 1.179763 0.5 2.5
% 4 gauss(1).area 0 0 0.0002898555 1e-08 1 ph/s/cm^2
% 5 gauss(1).center 0 0 1.936372 1.85051 1.9997 A
% 6 gauss(1).sigma 0 0 0.009638465 0.0001 0.1 A
% and show plots with the Gaussian in the model:
xrange(5.99,7.01); xlin; rplot_counts(3);
xrange(5.99,7.01); xlin; rplot_counts(4);
% Get conf ranges on the Gaussian parameters;
% allow the continuum to be adjusted as well.
_for ipar (4,6,1)
{
message("Finding conf for parameter "+string(ipar)+" ...");
(clo, chi)=conf(ipar, 2);
message(" 99\% range is: "+string(clo)+" to "+string(chi));
};
% For obs_hplaw :
% output of conf ranges for Gaussian fit of Fe line,
% known input values were: 0.0003, 1.9361, and 0.0085.
%
% Finding conf for parameter 4 ...
% 99% range is: 0.000206056 to 0.000378273
% Finding conf for parameter 5 ...
% 99% range is: 1.93279 to 1.93999
% Finding conf for parameter 6 ...
% 99% range is: 0.00616031 to 0.0136647
% Convert the sigma and center to a velocity FWHM:
variable line_center = get_par(5);
message(" Line width, c*dE(FWHM)/E, is : " +
string(Const_c*1.e-5*(2.35*get_par("gauss(1).sigma"))/line_center)+" km/s." );
% Convert chi, clo (these will be from the last parameter which was sigma)
% into a line width, FWHM in km/s.
% If the clo is very small then set it to 0:
if (clo < 0.0005) clo = 0.0;
message(" Line width, c*dE(FWHM)/E, is between: " +
string(Const_c*1.e-5*(2.35*clo)/line_center)+
" and "+string(Const_c*1.e-5*(2.35*chi)/line_center)+" km/s." );
% For obs_hplaw we get close to the 3100 km/s expected value:
% Line width, c*dE(FWHM)/E, is : 3506.77 km/s.
% With 99% conf range:
% Line width, c*dE(FWHM)/E, is between: 2241.31 and 4971.65 km/s.
%.........
% If this analysis is run on the spatially extended simulation, obs_hpldisk,
% then a double-peaked line is seen due to the spatial extent of the
% source and the Gaussian fit to it gives a sigma value much larger than
% the expected 0.0085 A :
% 6 gauss(1).sigma 0 0 0.0209564 <---
% and in velocity FWHM:
% Line width, c*dE(FWHM)/E, is : 7623.81 km/s.
% Line width, c*dE(FWHM)/E, is between: 5568.2 and 10607.4 km/s.
%
% - - - end of demo - - -
