Fitting Grating Data
OverviewLast Update: 1 Dec 2006 - reviewed for CIAO 3.4: no changes Synopsis: This thread provides a general introduction to fitting grating data in Sherpa. Loading and filtering data are covered, as well as defining instrument and source models. Users working with HRC-S/LETG grating data will also find the Fitting Multiple Orders of HRC-S/LETG Data thread helpful for their analysis. |
Contents
- Getting Started
- Reading the Spectrum Files
- Building the Instrument Responses
- Filtering the Data
- Defining the Source and Background Models
- Examining Method & Statistic Settings
- Fitting
- Examining Fit Results
- Saving and Quitting the Session
- History
- Images
Getting Started
Sample ObsID used: 459 (HETG/ACIS-S, 3C 273)
The files used in this example were created by following several of the CIAO Grating threads:
- Obtain Grating Spectra from HETG/ACIS-S Data
- Compute HETG/ACIS-S Grating ARFs
- Grouping PHA Data before Fitting
Here is a list of all the necessary files:
spectra: 459_heg_m1_bin10.pha 459_heg_p1_bin10.pha 459_meg_m1_bin10.pha 459_meg_p1_bin10.pha gARFs: 459_heg_m1.arf 459_heg_p1.arf 459_meg_m1.arf 459_meg_p1.arf
The spectrum that will be used in this session has been binned by a factor of 10.
Users may also choose to run the Create Grating RMFs for ACIS Observations thread. Creating observation-specific gRMFs is optional, and is discussed further in the Building the Instrument Responses section.
The data files are available in sherpa.tar.gz, as explained in the Sherpa Getting Started thread.
Reading the Spectrum Files
The data are input to Sherpa with the data command (a shorthand version of "read data"):
sherpa> data 1 459_heg_m1_bin10.pha The inferred file type is PHA. If this is not what you want, please specify the type explicitly in the data command. Warning: could not find SYS_ERR column WARNING: statistical errors specified in the PHA file. These are currently IGNORED. To use them, type: READ ERRORS "<filename>[cols CHANNEL,STAT_ERR]" fitsbin WARNING: backgrounds UP and DOWN are being read from this file, and are being combined into a single background dataset. Warning: could not find SYS_ERR column sherpa> data 2 459_heg_p1_bin10.pha The inferred file type is PHA. If this is not what you want, please specify the type explicitly in the data command. Warning: could not find SYS_ERR column WARNING: statistical errors specified in the PHA file. These are currently IGNORED. To use them, type: READ ERRORS "<filename>[cols CHANNEL,STAT_ERR]" fitsbin WARNING: backgrounds UP and DOWN are being read from this file, and are being combined into a single background dataset. Warning: could not find SYS_ERR column sherpa> data 3 459_meg_m1_bin10.pha The inferred file type is PHA. If this is not what you want, please specify the type explicitly in the data command. Warning: could not find SYS_ERR column WARNING: statistical errors specified in the PHA file. These are currently IGNORED. To use them, type: READ ERRORS "<filename>[cols CHANNEL,STAT_ERR]" fitsbin WARNING: backgrounds UP and DOWN are being read from this file, and are being combined into a single background dataset. Warning: could not find SYS_ERR column sherpa> data 4 459_meg_p1_bin10.pha The inferred file type is PHA. If this is not what you want, please specify the type explicitly in the data command. Warning: could not find SYS_ERR column WARNING: statistical errors specified in the PHA file. These are currently IGNORED. To use them, type: READ ERRORS "<filename>[cols CHANNEL,STAT_ERR]" fitsbin WARNING: backgrounds UP and DOWN are being read from this file, and are being combined into a single background dataset. Warning: could not find SYS_ERR column
Sherpa now refers to the spectra as follows:
- HEG, -1 order = dataset 1
- HEG, +1 order = dataset 2
- MEG, -1 order = dataset 3
- MEG, +1 order = dataset 4
Building the Instrument Responses
First, the instrument models are established by the rsp command. The arf parameter value is then set to the corresponding file for each order and arm:
sherpa> paramprompt off sherpa> rsp[hegm1] sherpa> rsp[hegp1] sherpa> rsp[megm1] sherpa> rsp[megp1] sherpa> hegm1.arf = 459_heg_m1.arf sherpa> hegp1.arf = 459_heg_p1.arf sherpa> megm1.arf = 459_meg_m1.arf sherpa> megp1.arf = 459_meg_p1.arf
This message will be printed after each gARF is entered:
The inferred file type is ARF. If this is not what you want, please specify the type explicitly in the data command.
In order to convolve the input datasets with the response model components that have been established, they must be defined as the instrument models. This involves pairing up the gARF and spectrum for each order, via the instrument command:
sherpa> instrument 1 = hegm1 sherpa> instrument 2 = hegp1 sherpa> instrument 3 = megm1 sherpa> instrument 4 = megp1
The current definition of the instrument model may be examined using show instrument:
sherpa> show instrument 1 Instrument 1: rsp1d[hegm1] Param Type Value ----- ---- ----- 1 rmf string: "none" (N_E=8192,N_PHA=8192) 2 arf string: "459_heg_m1.arf" (N_E=8192)
Notice that Sherpa has defined properties for the rmf parameter, even though we did not enter a file. Sherpa has support for "dummy" instruments: if data have been input and the instrument stack contains only an ARF, a dummy RMF will be created that maps the ARF bins to the data bins, if possible. ahelp instrument contains more information on "dummy" instruments.
The datasets may now be plotted:
sherpa> lplot 4 data 1 data 2 data 3 data 4
Figure 1 shows the resulting plot.
Filtering the Data
We choose to filter the data to focus on an area of interest:
The ignore command is used to ignore all the data in every dataset, then notice is used to select the desired regions. You may wish to adjust the limits to exclude more or less of your data.
Each filtered dataset may then be plotted:
sherpa> lplot 4 data 3 data 4 data 3 data 4
Notice that the plot now includes only the data in the specified wavelength regions. Figure 2 shows the resulting plot.
Defining the Source and Background Models
We plan on simultaneously fitting the background data (rather than subtracting it), so we need to create a model expression for the source and the background. We model this source with a broken power law (bpl1d) absorbed by the interstellar medium (atten). The background will be modeled by a one-dimensional power law (powlaw1d), also absorbed by the ISM (the same atten model).
First, we set up each model component. The absorption model will be referred to as "abs", the broken power law will be "bpow", and the 1D power law will be "pow1d".
sherpa> paramprompt on Model parameter prompting is on sherpa> atten[abs] abs.hcol parameter value [1e+20] abs.heiRatio parameter value [0.1] abs.heiiRatio parameter value [0.01] sherpa> bpl1d[bpow] bpow.gamma1 parameter value [0] bpow.gamma2 parameter value [0] bpow.eb parameter value [7.99625] bpow.ref parameter value [7] bpow.ampl parameter value [0.0238299] sherpa> powlaw1d[pow1d] pow1d.gamma parameter value [1] pow1d.ref parameter value [7] pow1d.ampl parameter value [0.0238299]
Note that since a dataset has already been input, Sherpa estimates the initial parameter values for the models based on the data. These values can also be listed with the show command:
sherpa> show models ------------------------------------------- Defined source/background model components: ------------------------------------------- Atten[abs] (integrate: off) Param Type Value Min Max Units ----- ---- ----- --- --- ----- 1 hcol thawed 1e+20 1e+17 1e+24 2heiRatio thawed 1e-01 0 1 3heiiRatio thawed 1e-02 0 1 bpl1d[bpow] (integrate: on) Param Type Value Min Max Units ----- ---- ----- --- --- ----- 1 gamma1 thawed 0 -10 10 2 gamma2 thawed 0 -10 10 3 eb thawed 7.9963 1.0075 14.985 4 ref frozen 7 1 14.985 5 ampl thawed 2.383e-02 2.383e-04 2.383 powlaw[pow1d] (integrate: on) Param Type Value Min Max Units ----- ---- ----- --- --- ----- 1 gamma thawed 1 -10 10 2 ref frozen 7 1 14.985 3 ampl thawed 2.383e-02 2.383e-04 2.383
Next we modify the initial parameter value for abs.hcol:
sherpa> abs.hcol=1.81e20 sherpa> freeze abs
The hydrogen column density (hcol) is set to the Galactic value. All the abs parameters are then frozen, which means they will not be allowed to vary during the fit.
Now that the model components have been established, the product of abs and bpow is assigned as the source model for all datasets:
sherpa> source 1:4 = abs*bpow
while the background model is set as the product of abs and pow1d:
sherpa> background 1:4 = abs*pow1d
Both model definitions can be listed with the show command:
sherpa> show source 1 Source 1: (abs * bpow) Atten[abs] (integrate: off) Param Type Value Min Max Units ----- ---- ----- --- --- ----- 1 hcol frozen 1.81e+20 1e+17 1e+24 2heiRatio frozen 1e-01 0 1 3heiiRatio frozen 1e-02 0 1 bpl1d[bpow] (integrate: on) Param Type Value Min Max Units ----- ---- ----- --- --- ----- 1 gamma1 thawed 0 -10 10 2 gamma2 thawed 0 -10 10 3 eb thawed 7.9963 1.0075 14.985 4 ref frozen 7 1 14.985 5 ampl thawed 2.383e-02 2.383e-04 2.383 sherpa> show background 1 Background 1: (abs * pow1d) Atten[abs] (integrate: off) Param Type Value Min Max Units ----- ---- ----- --- --- ----- 1 hcol frozen 1.81e+20 1e+17 1e+24 2heiRatio frozen 1e-01 0 1 3heiiRatio frozen 1e-02 0 1 powlaw[pow1d] (integrate: on) Param Type Value Min Max Units ----- ---- ----- --- --- ----- 1 gamma thawed 1 -10 10 2 ref frozen 7 1 14.985 3 ampl thawed 2.383e-02 2.383e-04 2.383
Examining Method & Statistic Settings
Next we check the current method and statistics settings:
sherpa> show method Optimization Method: Levenberg-Marquardt Name Value Min Max Description ---- ----- --- --- ----------- 1 iters 2000 1 10000 Maximum number of iterations 2 eps 1e-03 1e-09 1 Absolute accuracy 3 smplx 0 0 1 Refine fit with simplex (0=no) 4 smplxep 1 1e-04 1000 Switch-to-simplex eps factor 5 smplxit 3 1 20 Switch-to-simplex iters factor sherpa> show statistic Statistic: Chi-Squared Gehrels
For this fit, the default fitting and statistic settings will be used. More information is available from ahelp lev-mar and ahelp chigehrels. For a list of all the available methods and statistic settings, see ahelp method and ahelp statistic, respectively.
Fitting
The datasets are now fit:
sherpa> fit LVMQT: V2.0 LVMQT: initial statistic value = 1.06211e+10 LVMQT: final statistic value = 789026 at iteration 114 bpow.gamma1 0.420507 bpow.gamma2 -0.0786515 bpow.eb 11.7909 bpow.ampl 0.00225955 pow1d.gamma 0.276206 pow1d.ampl 0.000238299 WARNING: The value of pow1d.ampl within 0.01% of the pow1d.ampl.min limit boundary. You may wish to consider changing min/max values and refitting.
As the warning says, we reset the minimum boundary of powbkgmh.ampl and refit the data:
sherpa> pow1d.ampl.min=2.383e-10 sherpa> fit LVMQT: V2.0 LVMQT: initial statistic value = 789026 LVMQT: final statistic value = 1954.98 at iteration 6 bpow.gamma1 0.414261 bpow.gamma2 -0.0607346 bpow.eb 11.7837 bpow.ampl 0.00237397 pow1d.gamma 0.21045 pow1d.ampl 9.79897e-06
To plot the fits:
sherpa> lplot 4 fit 3 fit 4 fit 3 fit 4 sherpa> d 1,3,4 ylabel "" sherpa> title "3C 273 (ObsID 459)" sherpa> d 1 label 12 0.075 "HEG -1" sherpa> d 2 label 12 0.075 "HEG +1" sherpa> d 3 label 12 0.125 "MEG -1" sherpa> d 4 label 12 0.125 "MEG +1" sherpa> d all l all green sherpa> redraw
The ChIPS commands are used to add labels to the drawing areas. The plot is shown in Figure 3 .
It is also useful to plot the fit with the residuals:
sherpa> lplot 2 fit 1 delchi
This plot is shown in Figure 4 . After creating a plot, it may be saved as a PostScript file:
sherpa> print postfile 459_1_fit_delchi.ps
Examining Fit Results
There are several methods available in Sherpa for examining fit results. The goodness command reports information on the chi-square goodness-of-fit:
sherpa> goodness Goodness: computed with Chi-Squared Gehrels DataSet 1: 561 data points -- 555 degrees of freedom. Statistic value = 473.223 Probability [Q-value] = 0.994862 Reduced statistic = 0.852654 Background for DataSet 1: 561 data points -- 559 degrees of freedom. Statistic value = 114.649 Probability [Q-value] = 1 Reduced statistic = 0.205097 DataSet 2: 561 data points -- 555 degrees of freedom. Statistic value = 499.407 Probability [Q-value] = 0.956151 Reduced statistic = 0.899832 Background for DataSet 2: 561 data points -- 559 degrees of freedom. Statistic value = 98.6341 Probability [Q-value] = 1 Reduced statistic = 0.176447 DataSet 3: 281 data points -- 275 degrees of freedom. Statistic value = 279.114 Probability [Q-value] = 0.419589 Reduced statistic = 1.01496 Background for DataSet 3: 281 data points -- 279 degrees of freedom. Statistic value = 119.737 Probability [Q-value] = 1 Reduced statistic = 0.429164 DataSet 4: 281 data points -- 275 degrees of freedom. Statistic value = 260.716 Probability [Q-value] = 0.722871 Reduced statistic = 0.948058 Background for DataSet 4: 281 data points -- 279 degrees of freedom. Statistic value = 109.496 Probability [Q-value] = 1 Reduced statistic = 0.392459 Total : 3368 data points -- 3362 degrees of freedom. Statistic = 1954.98 Probability = 1 Reduced statistic = 0.581492
The uncertainty, covariance, and projection commands can be used to estimate confidence intervals for the thawed parameters:
sherpa> covariance Computed for sherpa.cov.sigma = 1 -------------------------------------------------------- Parameter Name Best-Fit Lower Bound Upper Bound -------------------------------------------------------- bpow.gamma1 0.414261 -0.00746802 +0.00746802 bpow.gamma2 -0.0607528 -0.119213 +0.119213 bpow.eb 11.7837 -0.243908 +0.243908 bpow.ampl 0.00237397 -8.36986e-06 +8.36986e-06 pow1d.gamma 0.21045 -0.0654054 +0.0654054 pow1d.ampl 9.79894e-06 -2.59661e-07 +2.59661e-07
Saving and Quitting the Session
Before exiting Sherpa, you may wish to save the session in order to return to the analysis at a later point:
sherpa> save all 459_fitting_session.shp
All the information about the current session is written to 459_fitting_session.shp, an ASCII file. It may be loaded into Sherpa again with the use command.
Finally, quit the session:
sherpa> quit
History
14 Jan 2005 | updated for CIAO 3.2: minor changes to screen output |
11 Jul 2005 | overall revision to thread, changes to screen output |
21 Dec 2005 | reviewed for CIAO 3.3: no changes |
01 Dec 2006 | reviewed for CIAO 3.4: no changes |