Last modified: 6 Dec 2022

URL: https://cxc.cfa.harvard.edu/sherpa/threads/marx/

Creating an Input Spectrum for MARX and ChaRT Simulations

Sherpa Threads (CIAO 4.16 Sherpa)


Overview

Synopsis:

This thread shows how to use Sherpa to create the source spectrum required as input to the two available Chandra PSF simulators, MARX and the Chandra Ray Tracer, ChaRT. One of the spectral specifications needed to run a MARX or ChaRT simulation is an unconvolved spectrum file, i.e., one containing the definition of the source spectrum at the entrance aperture of the HRMA for which a PSF should be simulated. This thread steps through an example Sherpa script which can be used to create such a spectrum for input to MARX or ChaRT. MARX and ChaRT are separate PSF simulators, but note the use of ChaRT requires MARX as a second step; the ChaRT output must be projected onto the detector plane with MARX; see the ChaRT threads for details.

MARX provides a detailed ray-trace simulation of how Chandra responds to a variety of astrophysical sources and can generate standard FITS event files and images as output. It contains detailed models for Chandra's High Resolution Mirror Assembly (HRMA), the HETG and LETG gratings, and all the focal plane detectors. ChaRT traces a collection of rays through the definitive Chandra X-ray optics model which are then projected onto the detector separately via MARX.

Last Update: 6 Dec 2022 - reviewed for CIAO 4.15, no content change


Contents


Creating the Spectrum

The spectrum input to MARX must be in units of energy (keV) versus photon flux density (photons/cm2/s/keV), and the input to ChaRT must be in units of energy (keV) versus photon flux (photons/cm2/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
-----------------------------------------------------

CIAO 4.15 Sherpa version 1 Thursday, October 13, 2022

Python 3.9.9 (main, Mar 24 2022, 22:02:45) 
Type 'copyright', 'credits' or 'license' for more information
IPython 8.0.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_spectrum contributed Sherpa package. For MARX, the Sherpa plot_source, get_source_plot, and save_arrays commands should 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

[A plot of photon flux (photons cm-2 s-1) versus energy (keV).]
[Print media version: A plot of photon flux (photons cm-2 s-1) versus energy (keV).]

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/cm2/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

Since the MARX software requires the input spectrum to be in photon flux density units (photons/s/cm2/keV), we use the Sherpa plot_source function and associated get_source_plot to view and access the spectral data arrays, since these functions automatically output the unconvolved spectrum in the required units.

sherpa> plot_source(lo=1.0, hi=8.0)

Figure 2 shows the resulting unconvolved model plot in the 1.0-8.0 keV energy range.

Figure 2: A plot of the spectrum input to MARX

[A plot of photon flux density (photons cm-2 s-1 keV-1) versus energy (keV).]
[Print media version: A plot of photon flux density (photons cm-2 s-1 keV-1) versus energy (keV).]

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 spectrum may be written to an ASCII file using the Sherpa save_arrays command, where the input x- and y- data arrays to be written to the file are taken from the plot of the spectrum we just created with the plot_source command. This is done by retrieving the energy and photon flux density arrays from the output of the get_source_plot command. Due to a bug in get_source_plot, the lo and hi arguments in the function cannot be used, and the entire energy grid defined by the ARF/RMF is used for the spectrum.

sherpa> set_analysis(1, "energy", "rate", factor=0) # reset to default flux units
sherpa> print(get_source_plot())
xlo    = [  0.1    0.11   0.12 ...,  10.97  10.98  10.99]
xhi    = [  0.11   0.12   0.13 ...,  10.98  10.99  11.  ]
y      = [1.1808e-18,4.2620e-15,1.6267e-12, ...,1.3460e-06,1.3434e-06,1.3408e-06]
xlabel = Energy (keV)
ylabel = f(E)  Photons/sec/cm^2/keV 
title  = Source Model of 3c273.pi
histo_prefs = {'xerrorbars': False, 'yerrorbars': False, 'ecolor': None, 'capsize': None, 'barsabove': False, 'xlog': False, 'ylog': False, 'linestyle': 'solid', 'drawstyle': 'default', 'color': None, 'alpha': None, 'marker': 'None', 'markerfacecolor': None, 'markersize': None, 'linecolor': None}


sherpa> fluxdensity = get_source_plot().y
sherpa> energy = get_source_plot().xhi
sherpa> save_arrays("source_flux_marx.dat",[energy,fluxdensity],["keV","photons/s/cm**2/keV"],ascii=True)

The ASCII file written by the above save_arrays command lists the energy upper-bin edge in keV in the first column, and the flux density in photons cm-2 s-1 keV-1 in the second column. Note that the MARX format for an input spectrum uses the upper-bin edge for energy and that the flux density in the first row may be any value, as the entry is not used.

unix% more source_flux_marx.dat
#TEXT/SIMPLE
# KEV PHOTONS/S/CM**2/KEV
1.0999999940395e-01 1.1807827096822e-18
1.1999999731779e-01 4.2620290331521e-15
1.2999999523163e-01 1.6266916315434e-12
1.4000000059605e-01 1.3034588671717e-10
1.5000000596046e-01 3.5103641593427e-09
1.5999999642372e-01 4.3072776655952e-08
1.7000000178814e-01 2.9121096116359e-07
1.8000000715256e-01 1.3006049697971e-06
1.8999999761581e-01 4.3489816498255e-06
2.0000000298023e-01 1.1403965973841e-05

...

10.90999984741 1.3644181869814e-06
10.92000007629 1.3617662004905e-06
10.93000030518 1.3591218732041e-06
10.93999958038 1.3564851398752e-06
10.94999980927 1.3538560078468e-06
10.96000003815 1.3512341620527e-06
10.97000026703 1.3486197809484e-06
10.97999954224 1.3460129606744e-06
10.98999977112 1.3434135485665e-06
11.0 1.3408213931501e-06

The sherpa_contrib.marx module can be used in lieu of the multi-step process outlined above, 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)

The save_marx_spectrum command is used to create the output spectrum file to be input to MRAX. 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.)


Summary

You are ready to run a MARX simulation once you have created a spectrum input file like example file source_flux_chart.dat, and have chosen an exposure time for the PSF simulation.


History

11 Mar 2010 original version
13 Jul 2010 updated for CIAO 4.2 Sherpa v2: removal of S-Lang version of thread
06 Jan 2012 reviewed for CIAO 4.4: updated Sherpa start-up banner
04 Dec 2013 reviewed for CIAO 4.6: no changes
20 Jan 2015 reviewed for CIAO 4.7, updated format for MARX
09 Dec 2015 reviewed for CIAO 4.8, updated to create input spectrum for ChaRT v2.
25 Jan 2016 updated to make use of fixed sherpa_contrib.chart Python module.
22 Nov 2016 reviewed for CIAO 4.9; note bug in get_source_plot which shows the entire source model energy grid, even when an energy range argument is applied to the function.
29 May 2018 reviewed for CIAO 4.10, no content change
29 Apr 2019 reviewed for CIAO 4.11, added information about the sherpa_contrib.marx module to generate MARX-compatible spectra.
11 Dec 2019 Re-ran with CIAO 4.12 to use Matplotlib rather than ChIPS for the plots.
04 Apr 2022 reviewed for CIAO 4.14, no content change
06 Dec 2022 reviewed for CIAO 4.15, no content change