Simulating NuSTAR X-Ray Spectra with Sherpa

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 Astronomers" page of the NuSTAR website (CalTech); the files used in this thread are:

• nustar_v2.0_80pct.arf - effective area modified by an 80% extraction region
• nustar_v2.0.rmf - includes detector efficiency

Figure 1: Effective Area Comparison

Figure 1: Effective Area Comparison

Plot of the effective area of the combined NuSTAR detector units, compared to those of current X-ray observatories¹ .

¹http://www.nustar.caltech.edu/for-astronomers

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> set_source("faked", xswabs.gal*(xsmekal.tplasma+xsmekal.tplasma2+xspowerlaw.powlaw))

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       100000 10^22 atoms / cm^2

sherpa> print(tplasma)
xsmekal.tplasma
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   t1.kT        thawed       1.1401       0.0808         79.9        keV
   t1.nH        frozen            1        1e-05        1e+19       cm-3
   t1.Abundanc  frozen            1            0         1000           
   t1.redshift  frozen        0.005            0           10           
   t1.switch    frozen            1            0            1           
   t1.norm      thawed  3.86263e-06            0        1e+24           

sherpa> print(tplasma2)
xsmekal.tplasma2
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   t2.kT        thawed     0.659763       0.0808         79.9        keV
   t2.nH        frozen            1        1e-05        1e+19       cm-3
   t2.Abundanc  frozen            1            0         1000           
   t2.redshift  frozen        0.005            0           10           
   t2.switch    frozen            1            0            1           
   t2.norm      thawed  2.62793e-05            0        1e+24           

sherpa> print print(powlaw)
xspowerlaw.powlaw
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   p1.PhoIndex  thawed      2.44436           -2            9           
   p1.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("nustar_v2.0_80pct.arf")

sherpa> nustar_rmf=unpack_rmf("nustar_v2.0.rmf")

sherpa> print(nustar_arf)
name     = nustar_v2.0_80pct.arf
energ_lo = Float64[820]
energ_hi = Float64[820]
specresp = Float64[820]
bin_lo   = None
bin_hi   = None
exposure = None

sherpa> print(nustar_rmf)
name     = nustar_v2.0.rmf
detchans = 820
energ_lo = Float64[820]
energ_hi = Float64[820]
n_grp    = UInt64[820]
f_chan   = UInt32[820]
n_chan   = UInt32[820]
matrix   = Float64[672400]
offset   = 0
e_min    = Float64[820]
e_max    = Float64[820]

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."

Running the Simulation with fake_pha

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 50 kiloseconds.

sherpa> fake_pha("faked", arf=nustar_arf, rmf=nustar_rmf, exposure=50000.)

This command associates the data/model ID "faked" with a simulated spectrum which is calculated over the grid defined by the NuSTAR instrument response, using the input exposure time and source model expression. Poisson noise is added to the modeled data.

We may use the show_data command to inspect the basic properties of the new data set; observe that the "name" field reads "faked", and the "exposure" field has been set to the chosen exposure time:

sherpa> show_data()
Data Set: faked
Filter: 3.0500-57.4500 Energy (keV)
Noticed Channels: 1-820
name           = faked
channel        = Float64[820]
counts         = Float64[820]
staterror      = None
syserror       = None
bin_lo         = None
bin_hi         = None
grouping       = Float64[820]
quality        = Float64[820]
exposure       = 50000.0
backscal       = None
areascal       = None
grouped        = True
subtracted     = False
units          = energy
rate           = True
plot_fac       = 0
response_ids   = [1]
background_ids = []

RMF Data Set: faked:1
name     = nustar_v2.0.rmf
detchans = 820
energ_lo = Float64[820]
energ_hi = Float64[820]
n_grp    = UInt64[820]
f_chan   = UInt32[820]
n_chan   = UInt32[820]
matrix   = Float64[672400]
offset   = 0
e_min    = Float64[820]
e_max    = Float64[820]

ARF Data Set: faked:1
name     = nustar_v2.0_80pct.arf
energ_lo = Float64[820]
energ_hi = Float64[820]
specresp = Float64[820]
bin_lo   = None
bin_hi   = None
exposure = None

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="faked")
           247.0

sherpa> calc_energy_flux(id="faked")
           1.2300167527680257e-13         # ergs/cm^2/sec/keV

sherpa> calc_photon_flux(id="faked")   # photons/cm^2/sec/keV
           1.1137887347387209e-05

The simulated data set may be visualized with the plot_data command, as shown below, and saved with print_window. Before plotting, we group the counts in the faked data spectrum 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("faked", 1)

sherpa> set_xlog()
sherpa> set_ylog()

sherpa> plot_data("faked")

sherpa> print_window("nustar_sim_data_50ksec.ps")

The resulting plot is shown Figure 2.

Figure 2: Plot of simulated NuSTAR source spectrum

Figure 2: Plot of simulated NuSTAR source spectrum

Plot of simulated hard X-ray NuSTAR imaging spectrum created with the Sherpa fake_pha command.

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 simulated data set to a PHA file, with a header containing the exposure time value and paths to the corresponding response files:

# NOTE: the 'get_data' lines below are necessary to work around a CIAO
# 4.4 bug with the save_pha() function; normally, one would
# just need to use save_pha() to write a simulated PHA data set to file.

sherpa> get_data("faked").header = {}                     # Make header an empty dictionary
sherpa> get_data("faked").name = "nustar_sim_50ksec.pha"  # Update name for save_pha

sherpa> save_pha("faked", "nustar_sim_50ksec.pha")

sherpa> !dmlist nustar_sim.pha header | grep EXPOSURE
0003 EXPOSURE                 14915.903481570            Real8

sherpa> !dmlist nustar_sim.pha header | grep ANCRFILE
0001 ANCRFILE             nustar_v2.0_80pct.arf          String 

sherpa> !dmlist nustar_sim.pha header | grep RESPFILE
0002 RESPFILE             nustar_v2.0.rmf                String 

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 choose to 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.

Before fitting, we check the energy range of the data set being considered in the analysis with the help of the the show_filter command:

sherpa> show_filter()
Data Set Filter: faked
3.000-60.0500 Energy (keV)

Then, we change the fit statistic from the default Chi2Gehrels to C-Stat, the fit method from Lev-Mar to Nelder-Mead, and change the model assigned to data set "faked" from gal*(tplasma+tplasma2+powlaw) to just powlaw.

sherpa> set_source("faked", powlaw)

sherpa> set_method("neldermead")
sherpa> set_stat("cstat")

sherpa> fit("faked")
Dataset               = faked
Method                = neldermead
Statistic             = cstat
Initial fit statistic = 78.4249
Final fit statistic   = 78.4249 at function evaluation 295
Data points           = 104
Degrees of freedom    = 102
Probability [Q-value] = 0.959991
Reduced statistic     = 0.768871
Change in statistic   = 1.27898e-13
   powlaw.PhoIndex   2.42955
   powlaw.norm    8.77709e-05

We can visualize the fit using the plot_fit command; first, we change the fit statistic back to Chi2Gehrels in order to calculate and view error bars overlaid on the data points in the plot.

sherpa> set_stat("chi2gehrels") 

sherpa> plot_fit("faked")

The resulting plot is shown in Figure 3.

Figure 3: Plot of power law fit to simulated NuSTAR source spectrum

Figure 3: Plot of power law fit to simulated NuSTAR source spectrum

Plot of power law model fit to the simulated NuSTAR X-ray spectrum.

Next, we change the source model for data set "faked" from powlaw to xskerrbb, set the fit statistic back to C-Stat, and re-fit the data:

sherpa> set_source("faked", xskerrbb.b1)
sherpa> b1.Mbh = 5.2

sherpa> set_stat("cstat")

sherpa> fit("faked")
Dataset               = faked
Method                = neldermead
Statistic             = cstat
Initial fit statistic = 97.4848
Final fit statistic   = 97.4848 at function evaluation 592
Data points           = 104
Degrees of freedom    = 100
Probability [Q-value] = 0.552548
Reduced statistic     = 0.974848
Change in statistic   = 3.43064e-06
   b1.a           0.9999      
   b1.Mbh         1.79041     
   b1.Mdd         3.5403      
   b1.norm        1.85661e-06 

This time, we visualize the new fit using the Kerr black hole model with plot_fit_delchi, in order to also view the fit delta chi residuals.

sherpa> set_stat("chi2gehrels")

sherpa> plot_fit_delchi("faked")
sherpa> current_plot("plot2")
sherpa> lin_scale(Y_AXIS)

The residuals plot is shown in Figure 4.

Figure 4: Plot of 'xskerrbb' model fit and delta chi residuals

Figure 4: Plot of 'xskerrbb' model fit and delta chi residuals

Plot of delta chi residuals resulting from a fit to the simulated NuSTAR spectrum, using a multi-temperature blackbody model for a thin accretion disk around a Kerr black hole.

The plot may be saved to a PostScript (or other format) file with the print_window command:

sherpa> print_window("nustar_sim_fit_delchi.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-sigma), 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> set_conf_opt("fast", False)

sherpa> conf("faked")
b1.a lower bound:	-0.159727
b1.a upper bound:	-----
b1.Mdd -: WARNING: The confidence level lies within (3.114890e-01, 3.114718e-01)
b1.Mdd lower bound:	-3.22882
b1.Mbh -: WARNING: The confidence level lies within (7.134948e-01, 7.134668e-01)
b1.Mbh lower bound:	-1.07693
b1.Mbh upper bound:	33.7966
b1.Mdd +: WARNING: The confidence level lies within (2.801152e+02, 2.801177e+02)
b1.Mdd upper bound:	276.576
b1.norm -: WARNING: The confidence level lies within (0.000000e+00, 5.221727e-07)
b1.norm lower bound:	-1.59553e-06
b1.norm +: WARNING: The confidence level lies within (2.279269e-03, 2.333536e-03)
b1.norm upper bound:	0.00230455
Dataset               = faked
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.a               0.9999    -0.159727        -----
   b1.Mbh            1.79041     -1.07693      33.7966
   b1.Mdd             3.5403     -3.22882      276.576
   b1.norm       1.85661e-06 -1.59553e-06   0.00230455

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 Lev-Mar method.

The output of the confidence command shows that the hard upper limit was hit for the b1.a parameter (which should have a value between -1 and 1). 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, systematic errors in the data, among other things.

The warning messages printed in the confidence output shown above occur where 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 for further exploring the quality of a fit, through visualization of fitted model parameter value(s) as a function of fit statistic value.