Creating an Input Spectrum for MARX and ChaRT Simulations

Creating the Spectrum

The spectrum input to MARX must be in units of energy (keV) versus photon flux density (photons/cm²/s/keV), and the input to ChaRT must be in units of energy (keV) versus photon flux (photons/cm²/s). Model source spectra in Sherpa are automatically output in the required units for both interfaces, using the appropriate model command for each; and those created with the ChaRT routines contributed to Sherpa may also be used as input to ChaRT, as shown in the section "Plotting and Saving the Spectrum."

For this example, we choose to fit a source data spectrum with a one-dimensional power-law multiplied by a photoelectric absorption model:

unix% sherpa
-----------------------------------------------------
Welcome to Sherpa: CXC's Modeling and Fitting Package
-----------------------------------------------------

Sherpa 4.17.0

Python 3.11.6 | packaged by conda-forge | (main, Oct  3 2023, 10:40:35) [GCC 12.3.0]
Type 'copyright', 'credits' or 'license' for more information
IPython 8.30.0 -- An enhanced Interactive Python. Type '?' for help.

IPython profile: sherpa
Using matplotlib backend: tkagg


sherpa> load_pha("3c273.pi")
sherpa> subtract()
sherpa> notice(0.4, 6.0)
dataset 1: 0.00146:14.9504 -> 0.3066:6.57 Energy (keV)

sherpa> set_source(xsphabs.abs1 * powlaw1d.p1)
sherpa> abs1.nh = 0.07
sherpa> freeze(abs1)
sherpa> guess(p1)
sherpa> fit()
Dataset               = 1
Method                = levmar
Statistic             = chi2gehrels
Initial fit statistic = 702.549
Final fit statistic   = 23.4063 at function evaluation 19
Data points           = 42
Degrees of freedom    = 40
Probability [Q-value] = 0.98311
Reduced statistic     = 0.585158
Change in statistic   = 679.143
   p1.gamma       2.12432      +/- 0.0791827   
   p1.ampl        0.000218513  +/- 1.44618e-05  

The example script above does the following:

• loads the source PHA data set "3c273.pi" into Sherpa and assigns it to default data set id=1 with load_pha;
• subtracts the background from the source data with subtract;
• restricts the energy range of the data set to 0.4-6.0 keV with notice;
• defines the unconvolved source model expression with set_source;
• freezes the photoelectric absorption model parameters so that they are not allowed to vary during the fit;
• guesses reasonable initial parameter values for the power-law model;
• fits the chosen source model to the data.

Note that the Sherpa fake_pha functionality is available for simulating a source spectrum independent of an actual data set; see the Sherpa simulation threads for details.

Plotting and Saving the Spectrum

In order to plot and save the best-fit source model spectrum to be input to ChaRT, use the functions plot_chart_spectrum and save_chart_spectrum, available in the chart contributed Sherpa package. For MARX, the functions plot_marx_spectrum and save_marx_spectrum, available in the marx contributed Sherpa package may be used.

ChaRT

sherpa> from sherpa_contrib.chart import *

After importing the chart_spectrum routines, the source model spectrum may be plotted in Sherpa with plot_chart_spectrum as follows, where the optional elow and ehigh values specify the minimum and maximum energies to include in the plot:

sherpa> plot_chart_spectrum()
sherpa> plot_chart_spectrum(elow=1.0, ehigh=8.0)

Figure 1 shows the resulting plot.

Figure 1: A plot of the spectrum input to ChaRT

Figure 1: A plot of the spectrum input to ChaRT

The source model is plotted in flux units of photons cm^-2 s^-1 versus energy in keV, as required by ChaRT.

The save_chart_spectrum command is used to create the output spectrum file to be input to ChaRT. In this example, we use the optional elow and ehigh arguments to restrict the output to the 1-8 keV range. If the bounds are omitted, the full energy range as determined by Sherpa is used.

sherpa> save_chart_spectrum("source_flux_chart.dat", elow=1.0, ehigh=8.0)
Created: source_flux_chart.dat
sherpa> quit()

The output produced by save_chart_spectrum is an ASCII file that contains three columns: the first two columns must be the lower and upper energy boundaries, in keV, and the third column is the photon flux in units of photon/cm²/sec; the values in the flux column must be >0.

unix% more source_flux_chart.dat
#TEXT/SIMPLE
# elo ehi spectrum
1.0 1.009999999999 1.8188423821135e-06
1.009999999999 1.019999999997 1.7889465849891e-06
1.019999999997 1.029999999996 1.7595017951804e-06
1.029999999996 1.039999999995 1.7306445484247e-06
1.039999999995 1.049999999993 1.7024994262716e-06
1.049999999993 1.059999999992 1.6749158520027e-06
1.059999999992 1.069999999990 1.6478842225815e-06
1.069999999990 1.079999999989 1.6211244927875e-06
1.079999999989 1.089999999988 1.5949101823919e-06
1.089999999988 1.099999999986 1.5694917864424e-06

...

7.899999999057 7.909999999055 2.7001154487312e-08
7.909999999055 7.919999999054 2.6928878494830e-08
7.919999999054 7.929999999053 2.6856887037872e-08
7.929999999053 7.939999999051 2.6785178639122e-08
7.939999999051 7.949999999050 2.6713750236105e-08
7.949999999050 7.959999999049 2.6642601973627e-08
7.959999999049 7.969999999047 2.6571730816490e-08
7.969999999047 7.979999999046 2.6501136919569e-08
7.979999999046 7.989999999044 2.6430818852304e-08
7.989999999044 7.999999999043 2.6360773619755e-08

---

MARX

The MARX software requires the input spectrum to be in photon flux density units (photons/s/cm²/keV) and the sherpa_contrib.marx module can be used in a smiliar manner as the sherpa_contrib.chart module.

sherpa> from sherpa_contrib.marx import *

After importing the marx_spectrum routines, the source model spectrum may be plotted in Sherpa with plot_marx_spectrum as follows, where the optional elow and ehigh values specify the minimum and maximum energies to include in the plot:

sherpa> plot_marx_spectrum()
sherpa> plot_marx_spectrum(elow=1.0, ehigh=8.0)

Figure 2: A plot of the spectrum input to MARX

Figure 2: A plot of the spectrum input to MARX

The source model is plotted in flux density units of photons cm^-2 s^-1 keV^-1 versus energy in keV, as required by MARX.

The save_marx_spectrum command is used to create the output spectrum file to be input to MARX. In this example, we use the optional elow and ehigh arguments to restrict the output to the 1-8 keV range. If the bounds are omitted, the full energy range as determined by Sherpa is used.

sherpa> save_marx_spectrum("source_flux_marx.dat", elow=1.0, ehigh=8.0)
Created: source_flux_marx.dat
sherpa> quit()

Scripting It

A script like the one used in this thread to create a source spectrum for input to MARX or ChaRT may be saved to an ASCII file—e.g., fit.py—and executed on the command line or in a future Sherpa session as follows:

sherpa> exec(open("fit.py").read())     # within Sherpa

or 

unix> sherpa fit.py            # at the Unix command line

The Sherpa script command may be used to save everything typed on the command line in a Sherpa session:

sherpa> script(filename="sherpa.log", clobber=False)

(Note that restoring a Sherpa session from such a file could be problematic since it may include syntax errors, unwanted fitting trials, et cetera.)