Simulating NuSTAR X-Ray Spectra with Sherpa
![[CXC Logo]](/ciao/imgs/cxc-logo.gif){kind=logo}
Sherpa Threads (CIAO 4.17 Sherpa)
Overview
Synopsis:
This thread illustrates the use of the Sherpa fake_pha command to simulate a hard X-ray emission spectrum of a supernova remnant which might be imaged by the detectors aboard the Nuclear Spectroscopic Telescopic Array (NuSTAR) satellite, launched in June 2012. NuSTAR, the first focusing high-energy X-ray mission sensitive across the 5–80 keV range, is aimed at studying objects such as black holes, supernovae, and extremely active galaxies. Here, we use Sherpa to simulate a 1-D PHA spectra using the currently available NuSTAR imaging ARF and RMF response files, in addition to a source model expression tailored to the Chandra X-ray spectrum of SN 1979C, suspected to harbor an accreting stellar-mass black hole.
Last Update: 6 Feb 2025 - Updated for CIAO 4.17, revised example to generate the spectra for the two co-aligned NuSTAR telescopes and the model is used to simultaneously fit the synthetic spectra.
Contents
Quick Start to Simulating Data
All that is required to simulate a 1-D PHA NuSTAR data set in Sherpa is the following:
- a source model expression defined with set_source
- proposal ARF & RMF instrument response files downloaded from the NuSTAR website
- exposure time for the simulation in seconds
fake_pha will do the rest:
unix% sherpa
# Search available Sherpa models.
sherpa> list_models()
['absorptionedge',
'absorptiongaussian',
'absorptionlorentz',
'absorptionvoigt',
...
'xszvphabs',
'xszwabs',
'xszwndabs',
'xszxipcf']
# Define a source model expression and assign it to an ID.
sherpa> emis = xsmekal.tplasma + xsmekal.tplasma2 + xspowerlaw.powlaw
sherpa> set_source("faked", xswabs.gal * emis)
# Set model parameter values.
sherpa> set_par(gal.nH, 0.025, frozen=True)
sherpa> tplasma.redshift = 0.005 # frozen by default
sherpa> tplasma.kT = 1.1401
sherpa> tplasma.norm = 3.86263e-06
sherpa> tplasma2.redshift = 0.005 # frozen by default
sherpa> tplasma2.kT = 0.659763
sherpa> tplasma2.norm = 2.62793e-05
sherpa> powlaw.PhoIndex = 2.44436
sherpa> powlaw.norm = 7.87378e-05
# Fake a PHA data set over the grid defined by the input response.
sherpa> fake_pha("fakedA", arf="point_30arcsecRad_1arcminOA.arf", rmf="nustar.rmf", exposure=350000., backscal=1.0)
sherpa> fake_pha("fakedB", arf="point_30arcsecRad_1arcminOA.arf", rmf="nustar.rmf", exposure=350000., backscal=1.0)
# Return information on the faked data set and associated responses.
sherpa> show_data("fakedA")
sherpa> show_data("fakedB")
sherpa> calc_data_sum(id="fakedA")
sherpa> calc_energy_flux(id="fakedA")
sherpa> calc_photon_flux(id="fakedA")
sherpa> calc_data_sum(id="fakedB")
sherpa> calc_energy_flux(id="fakedB")
sherpa> calc_photon_flux(id="fakedB")
# Plot simulated data.
sherpa> plot_data("faked")
# Save simulated data.
sherpa> save_pha("fakedA", "nustar_FPMA_sim_350ksec.pha")
sherpa> save_pha("fakedB", "nustar_FPMB_sim_350ksec.pha")
This represents the quickest and simplest way to simulate a NuSTAR PHA data set in Sherpa. Simply start CIAO and Sherpa; define a source model expression in Sherpa and assign it to a data set ID such as "faked"; set model parameter values; and then run fake_pha according to your specifications. Read on to find a detailed explanation of each step of this quick-start recipe.
Simulating Data Step by Step:
Downloading Calibration Response Files
In order to simulate a NuSTAR X-ray spectrum, it is necessary to define an instrument response using the appropriate ARF and RMF files currently available for the mission. NuSTAR response files are available for download from the "For Proposers: NuSTAR Response Files" page of the NuSTAR website (CalTech); the files used in this thread are:
- point_30arcsecRad_1arcminOA.arf — unweighted effective area for a 1' off-axis circular extraction region of radius 30''
- nustar.rmf — RMF includes detector efficiency
Figure 1: Effective Area Comparison
![[Effective Area Comparison]](nustar_effarea_CY17.jpg)
Figure 1: Effective Area Comparison
Plot of the effective area of the combined NuSTAR detector
units, compared to those of current X-ray observatories1 .
1http://www.nustar.caltech.edu/page/researchers
Establishing the Source Model Expression for the Simulation
We will use the Sherpa fake_pha command to simulate a NuSTAR PHA spectrum using a defined source model expression and instrument response as input. Details on the functionality of fake_pha, with other examples of its usage, are available in the ahelp file and the other Sherpa simulation threads.
The source model chosen for the simulation is defined with the Sherpa set_source command, as follows, where it is assigned to a string data set ID, "faked".
sherpa> emis = xsmekal.tplasma + xsmekal.tplasma2 + xspowerlaw.powlaw
sherpa> set_source("fakedA", xswabs.gal * emis)
sherpa> set_source("fakedB", gal * emis)
sherpa> set_par(gal.nH, 0.025, frozen=True)
sherpa> tplasma.redshift = 0.005 # frozen by default
sherpa> tplasma.kT = 1.1401
sherpa> tplasma.norm = 3.86263e-06
sherpa> tplasma2.redshift = 0.005 # frozen by default
sherpa> tplasma2.kT = 0.659763
sherpa> tplasma2.norm = 2.62793e-05
sherpa> powlaw.PhoIndex = 2.44436
sherpa> powlaw.norm = 7.87378e-05
Here we use a source model consisting of a combination of two thermal MEKAL components and a power-law, with parameter values chosen to give a good approximation the Chandra imaging spectra of SN 1979C in the 0.6–3 keV range (Patnaude et al., 2011). We can view the model parameter settings in this case using the print command, as the Sherpa show_source command only displays model information for models already assigned to data sets.
sherpa> print(gal) xswabs.gal Param Type Value Min Max Units ----- ---- ----- --- --- ----- gal.nH frozen 0.025 0 1e+06 10^22 atoms / cm^2 sherpa> print(tplasma) xsmekal.tplasma Param Type Value Min Max Units ----- ---- ----- --- --- ----- tplasma.kT thawed 1.1401 0.0808 79.9 keV tplasma.nH frozen 1 1e-06 1e+20 cm-3 tplasma.Abundanc frozen 1 0 1000 tplasma.redshift frozen 0.005 -0.999 10 tplasma.switch frozen 1 0 1 tplasma.norm thawed 3.86263e-06 0 1e+24 sherpa> print(tplasma2) xsmekal.tplasma2 Param Type Value Min Max Units ----- ---- ----- --- --- ----- tplasma2.kT thawed 0.659763 0.0808 79.9 keV tplasma2.nH frozen 1 1e-06 1e+20 cm-3 tplasma2.Abundanc frozen 1 0 1000 tplasma2.redshift frozen 0.005 -0.999 10 tplasma2.switch frozen 1 0 1 tplasma2.norm thawed 2.62793e-05 0 1e+24 sherpa> print(powlaw) xspowerlaw.powlaw Param Type Value Min Max Units ----- ---- ----- --- --- ----- powlaw.PhoIndex thawed 2.44436 -3 10 powlaw.norm thawed 7.87378e-05 0 1e+24
Defining the Instrument Response for the Simulation
Note that the steps in this section are optional, as ARF and RMF (or RSP) filenames may be directly entered into the arf and rmf arguments of the fake_pha expression shown in the next section of this thread. However, if you would like to view the details of or modify the response used in the simulation, you should load the NuSTAR ARF and RMF responses into Sherpa as follows:
sherpa> nustar_arf = unpack_arf("point_30arcsecRad_1arcminOA.arf") sherpa> nustar_rmf = unpack_rmf("nustar.rmf") sherpa> print(nustar_arf) name = point_30arcsecRad_1arcminOA.arf energ_lo = Float64[4096] energ_hi = Float64[4096] specresp = Float64[4096] bin_lo = None bin_hi = None exposure = None ethresh = 1e-10 sherpa> print(nustar_rmf) name = nustar.rmf energ_lo = Float64[4096] energ_hi = Float64[4096] n_grp = UInt64[4096] f_chan = UInt64[1014481] n_chan = UInt64[1014481] matrix = Float64[5627679] e_min = Float64[4096] e_max = Float64[4096] detchans = 4096 offset = 0 ethresh = 1e-10
The response data is loaded into variables nustar_arf and nustar_rmf with the unpack_arf and unpack_rmf commands. These variables will be used to assign the instrument response to the faked data set in the next section of this thread, "Running the Simulation with fake_pha". The advantage of unpacking the response files to response data objects rather than having fake_pha open the response files is that it reduces the total system memory used—instead of repeatedly opening the same set of responses and storing each instance to memory separately, the existing unpacked response instances may be accessed directly.
Running the Simulation with fake_pha
All response files represent the response of one NuSTAR telescope. The two NuSTAR telescopes (usually labeled by their focal plane modules, FPMA and FPMB) have subtly different responses and we do not recommend explicitly combining data from FPMA and FPMB when performing data analysis (i.e., combining individual counts from both telescopes to make a single spectrum). We recommend that science users simulate the source spectrum twice and use simultaneous fitting when determining the feasibility of an observation.
The response file (nustar.rmf) is constructed from real data for sources near the optical axis of either telescope and varies relatively little between FPMA and FPMB.
Background files should be paired with similar extraction region sources (e.g., use the 30-arcescond background region for the 30-arcsecond extraction region ARF, not with the 60-arcsecond extraction region ARF).
–NuSTAR Response Files–
Simulating the NuSTAR spectrum with Sherpa involves convolving the chosen source model expression with the corresponding NuSTAR response, and applying Poisson noise to the counts predicted by the model. All of these steps are performed by the fake_pha command, which has four required arguments: data set or model ID, ARF, RMF, and exposure time. We set the fake_pha arf and rmf arguments to the nustar_arf/nustar_rmf variables defined in the previous section, and the exposure argument to a value of 350 kiloseconds. As suggested by the NuSTAR documentation, we simulate two sets of spectra for FPMA and FPMB:
sherpa> fake_pha("fakedA", arf=nustar_arf, rmf=nustar_rmf, exposure=3.5e5)
sherpa> fake_pha("fakedB", arf=nustar_arf, rmf=nustar_rmf, exposure=3.5e5)
These commands associates the data/model ID "fakedA" and "fakedB" each with a simulated spectrum calculated over the grid defined by the NuSTAR instrument response, using the input exposure time and source model expression. Poisson noise is automatically added to the modeled data.
We may use the show_data command to inspect the basic properties of the new data sets; observe that the "name" field reads "faked", and the "exposure" field has been set to the chosen exposure time:
sherpa> show_data()
Data Set: fakedA
Filter: 1.6000-165.4400 Energy (keV)
Noticed Channels: 1-4096
name = faked
channel = Float64[4096]
counts = Float64[4096]
staterror = None
syserror = None
bin_lo = None
bin_hi = None
grouping = None
quality = None
exposure = 350000.0
backscal = None
areascal = None
grouped = False
subtracted = False
units = energy
rate = True
plot_fac = 0
response_ids = [1]
background_ids = []
RMF Data Set: fakedA:1
name = nustar.rmf
energ_lo = Float64[4096]
energ_hi = Float64[4096]
n_grp = UInt64[4096]
f_chan = UInt64[1014481]
n_chan = UInt64[1014481]
matrix = Float64[5627679]
e_min = Float64[4096]
e_max = Float64[4096]
detchans = 4096
offset = 0
ethresh = 1e-10
ARF Data Set: fakedA:1
name = point_30arcsecRad_1arcminOA.arf
energ_lo = Float64[4096]
energ_hi = Float64[4096]
specresp = Float64[4096]
bin_lo = None
bin_hi = None
exposure = None
ethresh = 1e-10
Data Set: fakedB
Filter: 1.6000-165.4400 Energy (keV)
Noticed Channels: 1-4096
name = faked
channel = Float64[4096]
counts = Float64[4096]
staterror = None
syserror = None
bin_lo = None
bin_hi = None
grouping = None
quality = None
exposure = 350000.0
backscal = None
areascal = None
grouped = False
subtracted = False
units = energy
rate = True
plot_fac = 0
response_ids = [1]
background_ids = []
RMF Data Set: fakedB:1
name = nustar.rmf
energ_lo = Float64[4096]
energ_hi = Float64[4096]
n_grp = UInt64[4096]
f_chan = UInt64[1014481]
n_chan = UInt64[1014481]
matrix = Float64[5627679]
e_min = Float64[4096]
e_max = Float64[4096]
detchans = 4096
offset = 0
ethresh = 1e-10
ARF Data Set: fakedB:1
name = point_30arcsecRad_1arcminOA.arf
energ_lo = Float64[4096]
energ_hi = Float64[4096]
specresp = Float64[4096]
bin_lo = None
bin_hi = None
exposure = None
ethresh = 1e-10
The calc_data_sum and calc_energy_flux/calc_photon_flux commands may be used to return the total counts and integrated energy/photon flux of the faked data set over the entire energy range (as in this example), or within a specified 'lo' to 'hi' range:
sherpa> calc_data_sum(id="fakedA")
431.0
sherpa> calc_data_sum(id="fakedB")
378.0
sherpa> calc_energy_flux(id="fakedA") # ergs/cm^2/sec/keV
2.0588096302472468e-13
sherpa> calc_energy_flux(id="fakedB")
2.0588096302472468e-13
sherpa> calc_photon_flux(id="fakedA") # photons/cm^2/sec/keV
2.902919956116689e-05
sherpa> calc_photon_flux(id="fakedB")
2.902919956116689e-05
![[TIP]](../../imgs/tip.png)
In Sherpa 4.17.0, fake_pha has issues accounting for the provided background spectrum if the provided RMF channel enumeration begin at 0, as seen in the provided NuSTAR responses, and will throw an error.
NuSTAR observations often include a significant background contribution to the spectra and depending on your science, should be accounted for in your simulation. As part of the downloaded NuSTAR simulation files package, background PHA spectral data are included. This background data can be appropriately scaled and included in the simulated spectra by using the bkg and backscal arguments in fake_pha as demonstrated below.
sherpa> nustar_bkg = unpack_bkg("bgd_30arcsec.pha")
WARNING: file 'nu10002020001A01_sr.arf' not found
WARNING: file 'nu10002020001A01_sr.rmf' not found
sherpa> fake_pha("faked", arf=nustar_arf, rmf=nustar_rmf, exposure=350000,
bkg=nustar_bkg, backscal=nustar_bkg.backscal)
The warnings thrown are due to non-existant response files referenced in the background PHA file header and will not affect the simulated spectrum. We can also see that the background data is incorporated into the simulated results using the same model by looking at the counts and fluxes compared to the earlier simulation excluding the background arguments.
sherpa> calc_data_sum(id="faked")
2487.0
sherpa> calc_energy_flux(id="faked") # ergs/cm^2/sec/keV
2.0588096302472468e-13
sherpa> calc_photon_flux(id="faked") # photons/cm^2/sec/keV
2.902919956116689e-05
The simulated data set may be visualized with the plot_data command, as shown below, and saved with plt.savefig. Before plotting, we group the counts in the faked data spectra such that each bin contains at least 1 count, and set both the X- and Y-axes of the plot to a logarithmic scale.
sherpa> group_counts("fakedA", 1)
dataset fakedA: 1.6:165.44 Energy (keV) (unchanged)
sherpa> group_counts("fakedB", 1)
dataset fakedB: 1.6:165.44 Energy (keV) (unchanged)
sherpa> set_xlog()
sherpa> set_ylog()
sherpa> plot_data("fakedA")
sherpa> plot_data("fakedB", overplot=True, alpha=0.5)
sherpa> plt.savefig("nustar_sim_data_350ksec.ps")
The resulting plot is shown Figure 2.
Figure 2: Plot of simulated NuSTAR source spectra
![[Plot of simulated NuSTAR source spectra]](sim_data.png)
Figure 2: Plot of simulated NuSTAR source spectra
Plot of simulated hard X-ray NuSTAR imaging spectra created with the Sherpa fake_pha command. The FPMA spectrum is in blue while the FPMB spectrum is in orange.
Note that had we not applied a customized, properly normalized source model expression to the faked data, we would have had to re-normalize the simulated data to match the known flux (or counts) of SN 1979C. This process is illustrated in the threads "Simulating Chandra ACIS-S Spectra with Sherpa" and "Simulating 1-D Data: the Sherpa fake_pha Command."
Writing the Simulated Data to Output Files
We may use the save_pha command to write the ungrouped simulated data sets to PHA files, with headers containing the exposure time value and paths to the corresponding response files:
sherpa> ungroup("fakedA")
dataset fakedA: 1.6:165.44 Energy (keV) (unchanged)
sherpa> ungroup("fakedB")
dataset fakedB: 1.6:165.44 Energy (keV) (unchanged)
sherpa> save_pha("faked", "nustar_FPMA_sim_350ksec.pha")
sherpa> save_pha("faked", "nustar_FPMB_sim_350ksec.pha")
sherpa> !dmkeypar nustar_FPMA_sim_350ksec.pha EXPOSURE echo+
350000.0
sherpa> !dmkeypar nustar_FPMA_sim_350ksec.pha ANCRFILE echo+
point_30arcsecRad_1arcminOA.arf
sherpa> !dmkeypar nustar_FPMA_sim_350ksec.pha RESPFILE echo+
nustar.rmf
Fitting the Simulated Data
A data set simulated with fake_pha may be fit as any other data set in Sherpa. For example, we can fit the simulated data with the same source model expression used to create it (and which is already assigned to it at this point), or define a different model. In this example, we will simultaneously fit the spectra from the two telescopes, first fit the full range of the data using only the power-law model component already assigned to it (minus the thermal plasma models), and then fit again using the xskerrbb model, a multi-temperature blackbody model for a thin accretion disk around a Kerr black hole.
While the detectors onboard NuSTAR may be sensitive to 1.5–165 keV photons, the mirrors only focus ~5–80 keV photons, so we restrict the data used to an energy band in this range with the notice command and group the spectral bins to have a minimum of one count using data only from the specified band. Before fitting, we check the energy range of the data sets being considered in the analysis with the help of the the show_filter command:
sherpa> notice(5, 80)
dataset fakedA: 1.6:165.44 -> 5:80 Energy (keV)
dataset fakedB: 1.6:165.44 -> 5:80 Energy (keV)
sherpa> group_counts("fakedA", 1)
dataset fakedA: 5:80 Energy (keV) (unchanged)
sherpa> group_counts("fakedB", 1)
dataset fakedB: 5:80 Energy (keV) (unchanged)
sherpa> show_filter()
Data Set Filter: fakedA
5.0000-80.0000 Energy (keV)
Data Set Filter: fakedB
5.0000-80.0000 Energy (keV)
If the data has been grouped, the filter uses the energies at the edges of the detector channels that contain the energies requested in the notice command. We then change:
- the fit statistic from the default \(\chi^{2}\)-Gehrels to C-Statistic
- the fit optimization method from Levenberg-Marquardt to Nelder-Mead Simplex
- change the model assigned to data sets "fakedA" and "fakedB" from gal*(tplasma+tplasma2+powlaw) to just powlaw
sherpa> set_source("fakedA", powlaw)
sherpa> set_source("fakedB", powlaw)
sherpa> set_method("neldermead")
sherpa> set_stat("cstat")
sherpa> fit()
Datasets = 'fakedA', 'fakedB'
Method = neldermead
Statistic = cstat
Initial fit statistic = 226.04
Final fit statistic = 225.239 at function evaluation 290
Data points = 283
Degrees of freedom = 281
Probability [Q-value] = 0.993791
Reduced statistic = 0.801563
Change in statistic = 0.800982
powlaw.PhoIndex 2.41477
powlaw.norm 7.10203e-05
We can visualize the fit using the plot_fit command:
sherpa> plot("fit","fakedA","fit","fakedB")
WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with cstat
WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with cstat
sherpa> fig = plt.gcf()
sherpa> fig.sca(fig.axes[0])
sherpa> plt.xlim(4,81)
sherpa> plt.title("Faked FPMA")
sherpa> fig.sca(fig.axes[1])
sherpa> plt.xlim(4,81)
sherpa> plt.title("Faked FPMB")
sherpa> fig.tight_layout()
The resulting plot is shown in Figure 3.
Figure 3: Plot of power-law fit to simulated NuSTAR source spectrum
![[Plot of fit to simulated NuSTAR spectrum]](sim_data_pl_fit.png)
Figure 3: Plot of power-law fit to simulated NuSTAR source spectrum
Plot of power-law model fit to the simulated NuSTAR X-ray spectra.
Next, we change the source model for the two synthetic data sets from the unabsorbed power-law to an absorbed multi-temperature, blackbody model of a thin disked Kerr black hole (xskerrbb)—using the contemporary xstbabs ISM absorption model—and re-fit the data with the prior knowledge that the source is ~5.2 M⊙ and ~15.2 Mpc away. We also change the fit optimizer to the Monte Carlo method to explore the entire set of thawed parameter-space while the distance and mass are fixed. We then switch back to Nelder-Mead Simple to speed up the exploration of parameter-space when the mass and distance parameters are thawed:
sherpa> set_source("fakedA", xstbabs.abs1*xskerrbb.b1)
sherpa> set_source("fakedB", abs1*b1)
sherpa> abs1.nh = gal.nh.val
sherpa> freeze(abs1)
sherpa> b1.Mbh = 5.2 # M_sun
sherpa> b1.Dbh = 15200 # in kpc
sherpa> b1.Dbh.hard_max = 16e3 # modify hard upper-limit
sherpa> freeze(b1.Mbh)
sherpa> freeze(b1.Dbh)
sherpa> thaw(b1.eta)
sherpa> thaw(b1.i)
sherpa> thaw(b1.hd)
sherpa> set_method("moncar")
sherpa> fit()
tbvabs Version 2.3
Cosmic absorption with grains and H2, modified from
Wilms, Allen, & McCray, 2000, ApJ 542, 914-924
Questions: Joern Wilms
joern.wilms@sternwarte.uni-erlangen.de
joern.wilms@fau.de
http://pulsar.sternwarte.uni-erlangen.de/wilms/research/tbabs/
PLEASE NOTICE:
To get the model described by the above paper
you will also have to set the abundances:
abund wilm
Note that this routine ignores the current cross section setting
as it always HAS to use the Verner cross sections as a baseline.
Datasets = 'fakedA', 'fakedB'
Method = moncar
Statistic = cstat
Initial fit statistic = 9582.6
Final fit statistic = 244.187 at function evaluation 7707
Data points = 283
Degrees of freedom = 277
Probability [Q-value] = 0.922871
Reduced statistic = 0.881543
Change in statistic = 9338.41
b1.eta 0.693668
b1.a 0.532788
b1.i 62.9249
b1.Mdd 9.18172
b1.hd 3.14982
b1.norm 4.5685
sherpa> thaw(b1.Dbh)
sherpa> b1.Dbh.max = 16e3
sherpa> b1.Dbh.min = 15e3
sherpa> thaw(b1.Mbh)
sherpa> b1.Mbh.max = 8
sherpa> b1.Mbh.min = 2
sherpa> set_method("neldermead")
sherpa> fit()
Datasets = 'fakedA', 'fakedB'
Method = neldermead
Statistic = cstat
Initial fit statistic = 244.187
Final fit statistic = 244.187 at function evaluation 680
Data points = 283
Degrees of freedom = 275
Probability [Q-value] = 0.909503
Reduced statistic = 0.887954
Change in statistic = 1.05214e-06
b1.eta 0.693669
b1.a 0.532787
b1.i 62.9249
b1.Mbh 5.19998
b1.Mdd 9.18172
b1.Dbh 15200
b1.hd 3.14981
b1.norm 4.5685
This time, we visualize the new simultaneous fit using the Kerr black hole model with plot_fit, plot_ratio, and matplotlib subplots.
sherpa> fig,axs = plt.subplots(2,2, sharex="all", figsize=(10,6), gridspec_kw={"height_ratios": [2.5,1]})
sherpa> plt.sca(axs[0][0])
sherpa> plot_fit("fakedA", clearwindow=False)
sherpa> plt.title("faked FPMA")
sherpa> plt.xlabel("")
sherpa> plt.sca(axs[0][1])
sherpa> plot_fit("fakedA", clearwindow=False)
sherpa> plt.title("faked FPMB")
sherpa> plt.xlabel("")
sherpa> plt.sca(axs[1][0])
sherpa> plot_ratio("fakedA", clearwindow=False)
sherpa> plt.yscale("linear")
sherpa> plt.title("")
sherpa> plt.sca(axs[1][1])
sherpa> plot_ratio("fakedB", clearwindow=False)
sherpa> plt.yscale("linear")
sherpa> plt.title("")
sherpa> plt.xlim(4,81)
sherpa> plt.tight_layout()
sherpa> plt.yscale('linear')
The residuals plot is shown in Figure 4.
Figure 4: Plot of 'xskerrbb' model fit and ratio to the best-fit
![[Plot of fit to simulated NuSTAR spectra]](sim_data_kerrbb_fit.png)
Figure 4: Plot of 'xskerrbb' model fit and ratio to the best-fit
Plot of the best-fit to the simulated NuSTAR spectra, using a multi-temperature blackbody model for a thin accretion disk around a Kerr black hole. The bottom plot shows the ratio of the data to the model values.
The plot may be saved to a PostScript (or other format) file with the plt.savefig command:
sherpa> plt.savefig("nustar_sim_fit_delchi_350ksec.ps")
Examining the Fit Results
We can explore the errors associated with the best-fit model parameters using the confidence command; the 68% confidence level values are returned by default (1σ), and may be changed with the set_conf_opt command. In this example, we run the confidence routine to check the errors on all of the thawed model parameters:
sherpa> freeze(b1.Mbh, b1.Dbh)
sherpa> set_conf_opt("fast", False)
sherpa> conf("fakedA","fakedB")
sherpa In [57]: conf()
b1.eta lower bound: -0.353782
b1.eta upper bound: -----
b1.a lower bound: -0.514535
b1.a upper bound: 0.324216
b1.hd -: WARNING: The confidence level lies within (1.263905e+00, 1.266405e+00)
b1.hd lower bound: -1.88466
b1.i -: WARNING: The confidence level lies within (5.515953e+01, 5.515939e+01)
b1.i lower bound: -7.76548
b1.i upper bound: 13.9066
b1.Mdd -: WARNING: The confidence level lies within (1.949507e+00, 1.954080e+00)
b1.Mdd lower bound: -7.22992
b1.norm -: WARNING: The confidence level lies within (8.003092e-01, 8.002591e-01)
b1.norm lower bound: -3.76822
b1.hd +: WARNING: The confidence level lies within (5.238380e+00, 5.238371e+00)
b1.hd upper bound: 2.08856
b1.Mdd +: WARNING: The confidence level lies within (1.303535e+02, 1.303533e+02)
b1.Mdd upper bound: 121.172
b1.norm +: WARNING: The confidence level lies within (2.697384e+01, 2.697378e+01)
b1.norm upper bound: 22.4053
WARNING: hard maximum hit for parameter b1.eta
Datasets = ('fakedA', 'fakedB')
Confidence Method = confidence
Iterative Fit Method = None
Fitting Method = neldermead
Statistic = cstat
confidence 1-sigma (68.2689%) bounds:
Param Best-Fit Lower Bound Upper Bound
----- -------- ----------- -----------
b1.eta 0.693669 -0.353782 -----
b1.a 0.532787 -0.514535 0.324216
b1.i 62.9249 -7.76548 13.9066
b1.Mdd 9.18172 -7.22992 121.172
b1.hd 3.14981 -1.88466 2.08856
b1.norm 4.5685 -3.76822 22.4053
The set_conf_opt command issued before conf in the example above ensures that the currently set fit optimization method is used in the confidence calculation, instead of the (faster) Sherpa default Levenberg-Marquardt method.
The output of the confidence command shows that the hard upper limit was hit for a parameter (denoted by the -----). This occurs when the parameter bound found by confidence lies outside the hard limit boundary for a model parameter, and could result from an issue with the signal-to-noise of the data, the applicability of the model to the data, or systematic errors in the data, among other things.
Warning messages can be printed in the confidence output when confidence cannot locate the minimum value of the fit statistic function, even though it is bracketed within an interval (perhaps due to poor resolution of the data or a discontinuity). In such cases, when the openinterval option of confidence is set to False (default), the confidence function will return the average of the open interval which brackets the minimum value.
The region-projection and interval-projection plotting commands are available to further explore the quality of a fit, through visualization of fitted model parameter value(s) as a function of fit statistic value.
Scripting It
The file fit.py is a Python script which performs the primary commands used above; it can be executed by typing %run -i fit.py on the Sherpa 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)
Reference
Patnaude, D.J.; Loeb, A.; Jones, C.; Evidence for a possible black hole remnant in the Type IIL Supernova 1979C; New Astronomy, Volume 16, Issue 3, p. 187-190, 2011 (doi:10.1016/j.newast.2010.09.004)
History
| 21 Jan 2011 | original version |
| 15 Dec 2011 | reviewed for CIAO 4.4: a work-around for a save_pha bug was added |
| 04 Dec 2013 | reviewed for CIAO 4.6: no changes |
| 20 Feb 2015 | updated for CIAO 4.7 and CY17 Chandra CfP utilizing the current NuSTAR proposal responses. |
| 15 Dec 2015 | updated for CIAO 4.8 and CY18 Chandra CfP utilizing the version 3 of the NuSTAR proposal responses. |
| 30 Nov 2016 | updated for CIAO 4.9, fit results updated; no new content. |
| 16 Jul 2018 | reviewed for CIAO 4.10: no content change |
| 11 Dec 2019 | Updated for CIAO 4.12: use matplotlib rather than ChIPS. |
| 27 Apr 2021 | Includes information on including background data to simulated spectrum. |
| 08 Apr 2022 | Updated for CIAO 4.14, use mask array to better group spectrum to selected energy band. |
| 13 Dec 2022 | Updated for CIAO 4.15, no content change. |
| 29 Nov 2023 | Updated for CIAO 4.16, use newly hard_max/hard_min parameter attribute in lieu of the hidden _max/_min attribute to set hard boundaries. |
| 06 Feb 2025 | Updated for CIAO 4.17, revised example to generate the spectra for the two co-aligned NuSTAR telescopes and the model is used to simultaneously fit the synthetic spectra. |