Creating an Input Spectrum for MARX and ChaRT Simulations
![[CXC Logo]](/ciao/imgs/cxc-logo.gif){kind=logo}
Sherpa Threads (CIAO 4.5 Sherpa v1)
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 Chandra X-ray optics which are then projected onto the detector separately via MARX.
Last Update: 06 Jan 2012 - reviewed for CIAO 4.4: updated Sherpa start-up banner
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 (chart_spectrum) 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.5 Sherpa version 1 Thursday, December 13, 2012 sherpa> load_pha("3c273.pi") sherpa> subtract() sherpa> notice(0.4, 6.0) 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 = 71.8133 Final fit statistic = 23.6368 at function evaluation 16 Data points = 42 Degrees of freedom = 40 Probability [Q-value] = 0.981499 Reduced statistic = 0.590919 Change in statistic = 48.1766 p1.gamma 2.13334 p1.ampl 0.000220644
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.
![[A plot of photon flux (photons cm-2 s-1) versus energy (keV).]](chart_spectrum.png)
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 two columns: the first column gives the center of the energy bin in keV, and the second column the flux in photons cm-2 s-1.
unix% more source_flux_chart.dat 1.005 1.82537e-06 1.015 1.7956e-06 1.025 1.76639e-06 1.035 1.73775e-06 1.045 1.70966e-06 1.055 1.68211e-06 1.065 1.6551e-06 1.075 1.62822e-06 1.085 1.60202e-06 1.095 1.57659e-06 ... 7.935 2.65464e-08 7.945 2.64753e-08 7.955 2.64045e-08 7.965 2.63339e-08 7.975 2.62637e-08 7.985 2.61937e-08 7.995 2.6124e-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.
![[A plot of photon flux density (photons cm-2 s-1 keV-1) versus energy (keV).]](marx_spectrum.png)
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, in the 1.0-8.0 keV energy range.
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.4874e-18 4.7047e-15 1.6421e-12 ..., 1.3301e-06 1.3275e-06
1.3250e-06]
xlabel = Energy (keV)
ylabel = f(E) Photons/sec/cm^2/keV
title = Source Model of 3c273.pi
histo_prefs = {'yerrorbars': False, 'linethickness': 2, 'symbolcolor': None, 'symbolfill': None, 'symbolsize': None, 'xlog': False, 'ylog': False, 'symbolangle': None, 'fillstyle': None, 'errthickness': None, 'errcolor': None, 'fillcolor': None, 'linecolor': 'red', 'errstyle': None, 'linestyle': 1, 'symbolstyle': 0, 'fillopacity': None}
sherpa> fluxdensity = get_source_plot(lo=1.0,hi=8.0).y
sherpa> energy = (get_source_plot(lo=1.0,hi=8.0).xhi+get_source_plot(lo=1.0,hi=8.0).xlo)/2
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 center of the energy bin in keV in the first column, and the flux density in photons cm-2 s-1 keV-1 in the second column.
unix% more source_flux_marx.dat #TEXT/SIMPLE # KEV PHOTONS/S/CM**2/KEV 1.0500000044703e-01 1.4873732915649e-18 1.1499999836087e-01 4.7047097549045e-15 1.2499999627471e-01 1.6421338065759e-12 1.3499999791384e-01 1.2506516519562e-10 1.4500000327826e-01 3.2724720025442e-09 1.5500000119209e-01 3.9784974304574e-08 1.6499999910593e-01 2.7327706694442e-07 1.7500000447035e-01 1.2097588802466e-06 1.8500000238419e-01 4.0465976847202e-06 1.9500000029802e-01 1.0629700083308e-05 2.0499999821186e-01 2.3142278930328e-05 2.1499999612570e-01 4.3510970419516e-05 2.2500000149012e-01 7.2854357693330e-05 2.3499999940395e-01 1.1116905417320e-04 2.4499999731779e-01 1.5708965279766e-04 2.5499999523163e-01 2.0807720610748e-04 2.6500000059605e-01 2.6390105715809e-04 2.7500000596046e-01 3.2129585650911e-04 2.8499999642372e-01 3.5519795041674e-04 2.9500000178814e-01 3.4419346894314e-04 3.0500000715256e-01 3.9076248146592e-04 . . .
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> execfile("fit.py") # 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, etcetera.)
The CXC is committed to helping Sherpa users transition to new syntax as smoothly as possible. If you have existing Sherpa scripts or save files, submit them to us via the CXC Helpdesk and we will provide the CIAO/Sherpa 4.5 syntax to you.
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 |
