Measuring Line Parameters with an LETG/ACIS-S Spectrum
Sherpa Threads (CIAO 4.9 Sherpa v1)
Overview
Synopsis:
After having created or downloaded a set of PHA2 and response files (RMFs and ARFs) for a LETG/ACIS-S observation, perform a simple fit to several of the line features present in the spectra. This thread takes the user through a calculation of the error bars on the line normalizations, positions, and widths, and shows how to calculate the line equivalent widths.
For those wishing to perform their analysis with XSPEC, the PHA2 file must first be split into individual type 1 PHA files, and then grouped to increase the signal-to-noise in each spectral bin. The Grouping a Grating Spectrum CIAO thread shows how to do this, and the procedure is also included within this thread.
Run this thread if:
You are working with a LETG/ACIS-S data set, and want to begin a simple analysis of the spectrum.
Related Links:
Last Update: 14 Dec 2015 - updated for CIAO 4.8, no content change
Contents
- Data Preparation
- Load the spectrum and responses
- Group and filter the data
- Defining the Source Models
- Fitting
- Adding lines to the model and fitting again
- Examine the fit results
- Calculate the equivalent width of a line
- Scripting It
- Summary
- History
- Images
Data Preparation
This thread makes use of archival data products for the star Sigma Geminorum (Sigma Gem for short), downloaded from TGcat. They were obtained by performing a "Search by ObsID" for ObsID 613, then choosing "view → file table" for the target, and clicking the download arrow next to each filename. These files were unzipped and placed in the directory in which the data analysis is to be performed.
Note that we are not considering any background files. For LETG/ACIS-S observations the background is suppressed by the process of order sorting. That is, an event at a given location along the grating arms must fall into a narrow range of energies (as determined by the CCD) in order for it to be considered a valid dispersed photon. Most background events do not meet these "order sorting" criteria. For this and many other LETG/ACIS-S observations, the first simple analysis can be performed without reference to the background. (Detailed analysis of low S/N spectra likely would require proper consideration of the background. See the Fitting Grating Data thread for an example of analysis including background data.)
Those using either Sherpa or ISIS to perform their analysis can use the data files directly in the analysis; skip to the "Load the spectrum and responses" section.
XSPEC users, however, must split and group the spectra before the analysis as described here.
Load the spectrum and responses
Start a Sherpa session:
unix% sherpa ----------------------------------------------------- Welcome to Sherpa: CXC's Modeling and Fitting Package ----------------------------------------------------- CIAO 4.9 Sherpa version 1 Friday, December 2, 2016 sherpa>
We begin by reading the PHA2 file that contains the grating spectra. This file contains 6 spectra, 3 spectral orders each for both positive and negative order LETG spectra. They are stored in the order: -3/-2/-1/1/2/3. If no filter is used when loading the data, then Sherpa will read in all 6 rows in order, assigning them to datasets 1-6.
We are only interested, however, in the first-order spectra (positive and negative); therefore, we use a Data Model-type filter to select the third and fourth rows of the PHA2 file (LETG -1 and LETG +1, respectively).
sherpa> load_pha("pha2[#row=3,4]") WARNING: systematic errors were not found in file 'pha2[#row=3,4]' statistical errors were found in file 'pha2[#row=3,4]' but not used; to use them, re-read with use_errors=True read background_up into a dataset from file pha2[#row=3,4] read background_down into a dataset from file pha2[#row=3,4] Multiple data sets have been input: 1-2
Alternatively, one could use a filter to select all negative and postitive first-order spectra via the TG_M keyword:
sherpa> load_pha("pha2[TG_M=-1,1]")
Note that associated background spectra (two each for each spectrum, one from either side of the gratings arm) will be read in with the source spectra. However, unless the subtract command is called, or unless one creates a background model, they will not be part of the fitting process. For the remainder of this thread, we ignore these background spectra.
We now read in the associated ARF files and assign them to their corresponding datasets, then do the same for the RMF files.
sherpa> load_arf(1,"leg_-1.arf") sherpa> load_arf(2,"leg_1.arf") sherpa> load_rmf(1,"leg_-1.rmf") sherpa> load_rmf(2,"leg_1.rmf")
Group and filter the data
Sherpa allows one to bin the spectra during the analysis session; see the Sherpra thread Changing the grouping scheme of a data set within Sherpa for details. We bin each of the spectra to have a minimum of 20 counts per energy channel.
sherpa> group_counts(1,20) sherpa> group_counts(2,20)
For the purposes of this analysis, we will restrict ourselves to only fitting the lines in the 1.4-2.2 keV region. (As we shall see below, this region alone shows multiple lines.) We restrict the energy range to 1.4-2.2 keV keV for both spectra using the notice command.
sherpa> notice(1.4,2.2)
The spectra can all be plotted in the same figure using the plot command. We select both plots to be considered "current" so that changing the axis will be applied to both of them, and then we choose the y-axes to be logarithmic.
sherpa> plot("data",1,"data",2) sherpa> current_plot("all") sherpa> log_scale(Y_AXIS) sherpa> print_window("all",["format","eps"])
The plot, shown in Figure 1, is also saved in EPS format to the file "all.eps".
Defining the Source Models
We will first fit a broad continuum to the spectra, without any lines. The XSPEC power-law model (xspowerlaw) is chosen for the continuum. It is assigned as the source for data set 1 and given the model name "powr"; the model name is then used to assign the same source model to data set 2.
sherpa> set_source(1, "xspowerlaw.powr") sherpa> set_source(2, powr) sherpa> show_source() Model: 1 xspowerlaw.powr Param Type Value Min Max Units ----- ---- ----- --- --- ----- powr.phoindex thawed 1 -2 9 powr.norm thawed 1 0 1e+24 Model: 2 xspowerlaw.powr Param Type Value Min Max Units ----- ---- ----- --- --- ----- powr.phoindex thawed 1 -2 9 powr.norm thawed 1 0 1e+24
Fitting
The fit statistic is set to be χ^{2} with the data counts acting as the variance (i.e., a Gaussian approximation to simple Poisson statistics). Both data sets are fit simultaneously.
sherpa> set_stat("chi2datavar") sherpa> fit() WARNING: data set 1 has associated backgrounds, but they have not been subtracted, nor have background models been set WARNING: data set 2 has associated backgrounds, but they have not been subtracted, nor have background models been set Datasets = 1, 2 Method = levmar Statistic = chi2datavar Initial fit statistic = 5.18961e+08 Final fit statistic = 4371.46 at function evaluation 19 Data points = 518 Degrees of freedom = 516 Probability [Q-value] = 0 Reduced statistic = 8.47182 Change in statistic = 5.18957e+08 powr.phoindex 2.73878 powr.norm 0.028444
At any time, the results of the last fit can be displayed using the show_fit command. The fit yields the following results:
sherpa> show_fit() Optimization Method: LevMar name = levmar ftol = 1.19209289551e-07 xtol = 1.19209289551e-07 gtol = 1.19209289551e-07 maxfev = None epsfcn = 1.19209289551e-07 factor = 100.0 verbose = 0 Statistic: Chi2DataVar Chi Squared with data variance. The variance in each bin is estimated from the data value in that bin. See also `Chi2Gehrels`, `Chi2XSpecVar` and `Chi2ModVar`. If the number of counts in each bin is large, then the shape of the Poisson distribution from which the counts are sampled tends asymptotically towards that of a Gaussian distribution, with variance sigma(i)^2 = N(i,S) + [A(S)/A(B)]^2 N(i,B) where N is the number of on-source (and off-source) bins included in the fit. The background term appears only if an estimate of the background has been subtracted from the data. Fit:Datasets = 1, 2 Method = levmar Statistic = chi2datavar Initial fit statistic = 5.18961e+08 Final fit statistic = 4371.46 at function evaluation 19 Data points = 518 Degrees of freedom = 516 Probability [Q-value] = 0 Reduced statistic = 8.47182 Change in statistic = 5.18957e+08 powr.PhoIndex 2.73878 powr.norm 0.028444
The results of these fits can then be plotted. Here we choose to display just the data and fits, not the residuals, and use logarithmic y-axes.
sherpa> plot("fit",1,"fit",2) sherpa> current_plot("all") sherpa> log_scale(Y_AXIS) sherpa> print_window("pl",["format","eps"])
The plot, shown in Figure 2, is also saved in EPS format to the file "pl.eps".
Adding lines to the model and fitting again
We start by creating a Gaussian model component for each of these features, and assign them to names indicative of their general energy location. We set each of their initial line widths to be narrow (i.e., 10^{-3}), their line normalizations to be small (10^{-4}), and their line energies to be close to the observed features. Additionally, to prevent the lines from broadening during the initial fits, we initially freeze the line widths.
The create_model_component command is used to establish each model component before setting the complex source expressions.
sherpa> create_model_component("xsgaussian","g15") <XSgaussian model instance 'xsgaussian.g15'> sherpa> g15.norm=1.e-4 sherpa> g15.sigma=1.e-3 sherpa> g15.linee=1.47 sherpa> create_model_component("xsgaussian","g18a") <XSgaussian model instance 'xsgaussian.g18a'> sherpa> g18a.norm=1.e-4 sherpa> g18a.sigma=1.e-3 sherpa> g18a.linee=1.84 sherpa> create_model_component("xsgaussian","g18b") <XSgaussian model instance 'xsgaussian.g18b'> sherpa> g18b.norm=1.e-4 sherpa> g18b.sigma=1.e-3 sherpa> g18b.linee=1.87 sherpa> create_model_component("xsgaussian","g2") <XSgaussian model instance 'xsgaussian.g2'> sherpa> g2.norm=1.e-4 sherpa> g2.sigma=1.e-3 sherpa> g2.linee=2
The sigma parameter of each model, which is the line width, is frozen before the fit:
sherpa> freeze(g15.sigma,g18a.sigma,g18b.sigma,g2.sigma)
We create a model with a power-law continuum plus these four line components, examine the source definition for one of the data sets (they are identical), then fit the data.
sherpa> set_source(1,"powr+g15+g18a+g18b+g2") sherpa> set_source(2,"powr+g15+g18a+g18b+g2") sherpa> show_source(1) Model: 1 ((((xspowerlaw.powr + xsgaussian.g15) + xsgaussian.g18a) + xsgaussian.g18b) + xs gaussian.g2) Param Type Value Min Max Units ----- ---- ----- --- --- ----- powr.phoindex thawed 2.73878 -2 9 powr.norm thawed 0.028444 0 1e+24 g15.linee thawed 1.47 0 1e+06 keV g15.sigma frozen 0.001 0 10 keV g15.norm thawed 0.0001 0 1e+24 g18a.linee thawed 1.84 0 1e+06 keV g18a.sigma frozen 0.001 0 10 keV g18a.norm thawed 0.0001 0 1e+24 g18b.linee thawed 1.87 0 1e+06 keV g18b.sigma frozen 0.001 0 10 keV g18b.norm thawed 0.0001 0 1e+24 g2.linee thawed 2 0 1e+06 keV g2.sigma frozen 0.001 0 10 keV g2.norm thawed 0.0001 0 1e+24 sherpa> fit() WARNING: data set 1 has associated backgrounds, but they have not been subtracted, nor have background models been set WARNING: data set 2 has associated backgrounds, but they have not been subtracted, nor have background models been set Datasets = 1, 2 Method = levmar Statistic = chi2datavar Initial fit statistic = 2423.51 Final fit statistic = 871.453 at function evaluation 93 Data points = 518 Degrees of freedom = 508 Probability [Q-value] = 1.40602e-21 Reduced statistic = 1.71546 Change in statistic = 1552.05 powr.phoindex 2.92076 powr.norm 0.0294724 g15.linee 1.47194 g15.norm 0.000185117 g18a.linee 1.84039 g18a.norm 7.23641e-05 g18b.linee 1.86453 g18b.norm 0.000138761 g2.linee 2.00539 g2.norm 0.000125723
Now that most of our parameters are close to their final values and we have a good description of the lines, we thaw the line widths, and then refit the data.
sherpa> thaw(g15.sigma,g18a.sigma,g18b.sigma,g2.sigma) sherpa> fit() WARNING: data set 1 has associated backgrounds, but they have not been subtracted, nor have background models been set WARNING: data set 2 has associated backgrounds, but they have not been subtracted, nor have background models been set Datasets = 1, 2 Method = levmar Statistic = chi2datavar Initial fit statistic = 871.453 Final fit statistic = 856.993 at function evaluation 235 Data points = 518 Degrees of freedom = 504 Probability [Q-value] = 9.89124e-21 Reduced statistic = 1.70038 Change in statistic = 14.4597 powr.PhoIndex 2.93174 powr.norm 0.0295794 g15.LineE 1.47194 g15.Sigma 0.000997626 g15.norm 0.000185197 g18a.LineE 1.84084 g18a.Sigma 0.00366799 g18a.norm 7.82044e-05 g18b.LineE 1.86456 g18b.Sigma 0.00256214 g18b.norm 0.00014157 g2.LineE 2.0049 g2.Sigma 0.000628811 g2.norm 0.000125951
The four fitted lines correspond to Mg XII, two lines of Si XIII, and Si XIV. The fitted line energies all lie within 40-130 km/sec of their expected rest energies. (One can check for the expected line energies, for example, using the Interactive Guide for ATOMDB, the atomic line database.) Below, we will return to the question of: to what extent are these small deviations significant?
We can now plot these fitted spectra, and look at the δχ residuals of the fit. After issuing the plot command, we set the y-axes for the data, but not the residuals, to be logarithmic. We do this by first setting "plot1" (LETG -1 spectra) to be current and then change the y-axis to logarithmic. We then set "plot2" (LETG +1 spectra) to be current and then change the y-axis to logarithmic.
sherpa> plot("fit",1,"fit",2,"delchi",1,"delchi",2) sherpa> current_plot("plot1") sherpa> log_scale(Y_AXIS) sherpa> current_plot("plot2") sherpa> log_scale(Y_AXIS) sherpa> print_window("lines",["format","eps"])
The plot, shown in Figure 3, is also saved in EPS format to the file "lines.eps".
Examine the fit results
In order to explore the errors on the line parameters, we can now use the conf() command to run the Sherpa confidence routines. As is common in X-ray astronomy, we search for the 90% confidence level values for one interesting parameter (i.e., the variation of the parameter of interest that when refitting the remaining unfrozen parameters yields a change in the χ^{2} value of 2.706). This corresponds to an error bar of 1.64σ (1σ being the 68% confidence interval), and set this as the default error value in the projection routines via the set_conf_opt command. We then run conf for the three parameters of each of the Gaussian model components: line energy, line width, and line normalization.
sherpa> set_conf_opt("sigma",1.64) sherpa> conf(g15.linee,g15.sigma,g15.norm,g18a.linee,g18a.sigma,g18a.norm,g18b.linee,g18b.sigma,g18b.norm,g2.linee,g2.sigma,g2.norm) b.linee,g18b.sigma,g18b.norm,g2.linee,g2.sigma,g2.norm) WARNING: data set 1 has associated backgrounds, but they have not been subtracted, nor have background models been set WARNING: data set 2 has associated backgrounds, but they have not been subtracted, nor have background models been set g15.LineE -: WARNING: The confidence level lies within (1.471676e+00, 1.471808e+00) g15.LineE lower bound: -0.00019856 g18b.LineE lower bound: -0.000391919 g18a.LineE lower bound: -0.000723807 g2.LineE -: WARNING: The confidence level lies within (2.002398e+00, 2.004898e+00) g2.LineE lower bound: -0.00125012 g15.LineE upper bound: 0.000217175 g18a.LineE upper bound: 0.000827731 g18b.LineE upper bound: 0.000386762 g18a.Sigma lower bound: -0.00138865 g18b.Sigma lower bound: -0.00211415 g15.Sigma -: WARNING: The confidence level lies within (8.706380e-05, 9.252185e-05) g15.Sigma lower bound: -0.000907834 g18a.Sigma upper bound: 0.00130755 g18a.norm lower bound: -7.50783e-06 g18b.Sigma upper bound: 0.00083865 g18b.norm lower bound: -8.15793e-06 g15.Sigma upper bound: 0.000597058 g15.norm lower bound: -1.05608e-05 g18a.norm upper bound: 7.95378e-06 g18b.norm upper bound: 8.11989e-06 g15.norm upper bound: 1.07036e-05 g2.LineE +: WARNING: The confidence level lies within (2.007843e+00, 2.007838e+00) g2.LineE upper bound: 0.00294218 g2.Sigma lower bound: ----- g2.Sigma upper bound: 0.00102108 g2.norm lower bound: -6.52364e-06 g2.norm upper bound: 6.55562e-06 Datasets = 1, 2 Confidence Method = confidence Iterative Fit Method = None Fitting Method = levmar Statistic = chi2datavar confidence 1.64-sigma (89.8995%) bounds: Param Best-Fit Lower Bound Upper Bound ----- -------- ----------- ----------- g15.LineE 1.47194 -0.00019856 0.000217175 g15.Sigma 0.000997626 -0.000907834 0.000597058 g15.norm 0.000185197 -1.05608e-05 1.07036e-05 g18a.LineE 1.84084 -0.000723807 0.000827731 g18a.Sigma 0.00366799 -0.00138865 0.00130755 g18a.norm 7.82044e-05 -7.50783e-06 7.95378e-06 g18b.LineE 1.86456 -0.000391919 0.000386762 g18b.Sigma 0.00256214 -0.00211415 0.00083865 g18b.norm 0.00014157 -8.15793e-06 8.11989e-06 g2.LineE 2.0049 -0.00125012 0.00294218 g2.Sigma 0.000628811 ----- 0.00102108 g2.norm 0.000125951 -6.52364e-06 6.55562e-06
Note that confidence failed to find a lower limit for the g2.Sigma parameter; we shall return to the reasons for this.
For those parameters that did return valid parameter bounds, we see that line energies are determined with accuracies ranging from 0.2 eV to 0.8 eV, which translate to equivalent velocities of 40-130 km/sec. In this energy range, the Half Width Half Maximum (HWHM) resolution of the LETG ranges from 700-900 km/sec. That is, due to the extremely good statistics of these very strong lines, the formal error bars are 7-18 times smaller than HWHM. Whereas it is quite possible with sufficiently good statistics to centroid a Gaussian to better than HWHM, such a large factor of improvement cautions against over-interpretation, and implies that systematic concerns may dominate. For example, if one were to fit the line with a Voigt profile instead of a Gaussian, systematically different results may be obtained. However, given the high S/N, we would expect that any systematic differences obtained from using alternative line profiles likely would be smaller than HWHM.
The expected energies for these lines are 1.4720 keV (20 km/sec offset), 1.8395 keV (150 km/sec offset), 1.8650 keV (70 km/sec offset) and 2.0049 keV (5 km/sec offset). That is, all are detected at their rest energies to well-within what we would consider to be reasonable systematic uncertainties.
Similar statements can be made about their fitted widths. The fitted widths range over equivalent velocities of 90-600 km/sec. The two broadest lines, i.e., at 1.841 keV and 1.865 keV, could also have a third line in between them (the expected intercombination line at 1.854 keV). At this point, it would be reasonable for the user to explore the limits that one could place on the presence of this line, and how it systematically affects the fits. For example, all three lines could be fit, but with their relative energies fixed to theoretical expectations. Likewise, one also could explore fixing their widths to be the same fractional width (under the assumption that they all come from the same velocity regions).
The line at 2.005 keV had a lower limit line width that reached 0 before a Δχ^{2} of 2.7 was found; therefore, the upper bound (1.6 eV) reasonably could be considered an upper-limit to the line width. This corresponds to a velocity of 250 km/sec. Given that this is a factor of a few below the HWHM of LETG, it would be prudent to explore the systematic dependence of this value upon model assumptions.
The fact that the line energy and width fit uncertainties are well-below the LETG HWHM resolution gives some indication as to why the confidence algorithm failed to find one of the 24 parameter limits. The confidence algorithm tries to determine the bounds by extrapolating the shape of the χ^{2} surface near the best-fit value of the parameter, under the assumption that it is relatively "well-behaved". Here, due to the very high S/N and being at the systematic limits of the detector, some of the χ^{2} contours are "badly behaved". Here we show a statistic vs. parameter curve for the line energy that failed to converge, by using the int_proj (interval projection) function.
For the line near 2 keV, we explore 50 values of the fit statistic in a narrow bound around the best-fit value.
sherpa> int_proj(g2.linee,min=2.0037,max=2.00785,nloop=50) sherpa> print_window("g2_int_proj",["format","eps"])
The plot, shown in Figure 4, is also saved in EPS format to the file "g2_int_proj.eps".
As can be seen, this statistic curve is very different from quadratic, and has large sudden jumps at boundaries near 2.004 keV and 2.008 keV (i.e., a width with an equivalent velocity of 600 km/sec, comparable to the HWHM of the LETG). Thus, one could trust the results of the confidence method and reasonably use this interval as the bound on the uncertainty of this line energy. Again, one should explore the systematic dependence of this line energy upon model assumptions.
On the other hand, the statistic contours for all the line amplitudes (not shown) are very smooth and nearly quadratic. For these parameters, we expect the results of projection to be a very good estimate of the uncertainties on the line strengths. All the line normalizations are determined to an accuarcy of 5-10% at the 90% confidence level for one interesting parameter; however, in terms of systematic effects, one could also add the expected intercombination line at 1.854 keV to see how that line affects the fitted normalizations of the other two (resonance and forbidden) lines near 1.8 keV.
Calculate the equivalent width of a line
Rather than merely presenting the line normalizations, we can instead calculate the line equivalent width. The line equivalent width is an often used measure for the line strength relative to the underlying continuum. Sherpa has a function for evaluating the equivalent width, eqwidth. The required inputs are the source model absent the line feature of interest and the source model with the line feature of interest.
Optionally, the data set number for which the equivalent width will be calculated can also be given as an input, but here as we apply the same model and parameters to both datasets, this is not necessary. There are also optional lo and hi arguments for restricting the calculation of the equivalent width to a subset of the full energy or wavelength range of a data set.
sherpa> eqwidth(powr,powr+g15) 0.0194461820785026417 sherpa> eqwidth(powr,powr+g18a) 0.01581991349977404 sherpa> eqwidth(powr,powr+g18b) 0.029733623727568638 sherpa> eqwidth(powr,powr+g2) 0.032760542190053712
The units of the equivalent width are the same as the units of the x-axis, in this case, keV. Thus the line equivalent widths range from 19 eV to 33 eV. Note that since the equivalent width is referenced to the underlying continuum, we only need input the continuum component (i.e., the power-law model) and the single line of interest. The answers would have been identical, however, if we had included the other narrow line components along with the power-law continuum.
Scripting It
The file fit.py is a Python script which performs the primary commands used above; it can be executed by typing exec(open("fit.py").read()) 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)
(Note that restoring a Sherpa session from such a file could be problematic since it may include syntax errors, unwanted fitting trials, et cetera.)
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.9 syntax to you.
Summary
As shown in this thread, working with gratings spectra is fundamentally no different than working with CCD spectra. Aside from the fact that gratings data will often begin with a PHA2 file containing multiple spectra, rather than a type 1 PHA file containing a single spectrum, the analysis paths can be very much the same for both types of spectra. (If one is performing analysis in either Sherpa or ISIS, the PHA2 file can be used directly; whereas, for XSPEC one has to use the intermediate step of creating type 1 PHA files in order to be able to bin the data.) The data are read in, response matrices are assigned, models are defined and fit, and error bars on the parameters are determined. Line equivalent widths can be calculated easily. Gratings spectra can even be slightly less complicated than CCD spectra in that ACIS-S/LETG spectra often do not require a background. (Again, order sorting of the spectra often ensures that the background is negligible.)
History
02 Apr 2009 | new for Sherpa 4.1 |
29 Apr 2009 | new script command is available with CIAO 4.1.2 |
11 Jan 2010 | updated for CIAO 4.2 |
13 Jul 2010 | updated for CIAO 4.2 Sherpa v2: removal of S-Lang version of thread. |
15 Feb 2011 | updated the URL for the Interactive Guide for ATOMDB (WebGUIDE) |
22 Jan 2012 | reviewed for CIAO 4.4 (no changes) |
13 Dec 2012 | updated for CIAO 4.5: group commands no longer clear the existing data filter |
10 Dec 2013 | reviewed for CIAO 4.6: no changes |
02 Mar 2015 | updated for CIAO 4.7, no content change |
14 Dec 2015 | updated for CIAO 4.8, no content change |