GUIDE: Fitting and Identifying Spectral Lines
OverviewLast Update: 1 Dec 2006 - updated for CIAO 3.4: Sherpa version Synopsis: GUIDE is a command line interpreted language that functions, essentially, as an extension of Sherpa. One of its more advanced applications is in identifying spectral lines to derive physical conditions and differential emmision measures. Purpose: To determine the flux and identity (ion and transition) of an emission line, and to write the results to a Model Descriptor List (MDL) file. Related Links:
|
Contents
- Get Started
- Loading GUIDE
- Read the Spectrum File and Build Responses
- Defining the Source Model
- Fitting
- Identify the Line
- Write an MDL File
- History
- Images
Get Started
Sample ObsID used: 1318 (HETG/ACIS-S, Capella)
In order to complete this thread, you will need grating ARFs for your dataset:
1318HEG_-1_garf.fits 1318HEG_1_garf.fits 1318MEG_-1_garf.fits 1318MEG_1_garf.fits acis_1318_pha2.fits
The Create an HETG/ACIS-S Grating ARF thread shows you how to do so (there are similar threads if you are working with LETG/ACIS-S, LETG/HRC-S, or LETG/HRC-I data). We are using only the 1st order spectra, which correspond to data elements 3 and 4 (-1, +1 for HEG) and 9 and 10 (-1, +1 for MEG) in the standard Level II PHA file. The Examining PHA2 Files thread has more information on identifying gratings and orders.
Loading GUIDE
Start Sherpa and load the GUIDE package:
unix% sherpa ----------------------------------------------------- Welcome to Sherpa: CXC's Modeling and Fitting Program ----------------------------------------------------- Version: CIAO 3.4 Type AHELP SHERPA for overview. Type EXIT, QUIT, or BYE to leave the program. Notes: Temporary files for visualization will be written to the directory: /tmp To change this so that these files are not deleted when you exit Sherpa, edit $ASCDS_WORK_PATH in your 'ciao' setup script. Abundances set to Anders & Grevesse sherpa> import("guide") GUIDE Initialized using ATOMDB v1.3.0
If a series of messages is printed here indicating that various files are not found, then it is likely that the ATOMDB is not correctly installed on your system. Please see your system manager or the CIAO download page and the ATOMDB webpage.
Read the Spectrum File and Build Responses
sherpa> paramprompt off Model parameter prompting is off sherpa> data acis_1318_pha2.fits The inferred file type is PHA Type II. 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 WARNING: multiple datasets have been input. The next available dataset number is 13. sherpa> analysis Analysis Space for Dataset 1: Wavelength Analysis Space for Dataset 2: Wavelength Analysis Space for Dataset 3: Wavelength Analysis Space for Dataset 4: Wavelength Analysis Space for Dataset 5: Wavelength Analysis Space for Dataset 6: Wavelength Analysis Space for Dataset 7: Wavelength Analysis Space for Dataset 8: Wavelength Analysis Space for Dataset 9: Wavelength Analysis Space for Dataset 10: Wavelength Analysis Space for Dataset 11: Wavelength Analysis Space for Dataset 12: Wavelength
By default, the analysis mode is set to wavelength when reading in a type II PHA file. For analysis in wavelength space to make any sense, however, a wavelength grid must be created; this is done when the instrument response (either ARF, RMF, or both) is defined.
sherpa> instrument 3 = rsp[hegm1](,1318HEG_-1_garf.fits,) The inferred file type is ARF. If this is not what you want, please specify the type explicitly in the data command. sherpa> instrument 4 = rsp[hegp1](,1318HEG_1_garf.fits,) The inferred file type is ARF. If this is not what you want, please specify the type explicitly in the data command. sherpa> instrument 9 = rsp[megm1](,1318MEG_-1_garf.fits,) The inferred file type is ARF. If this is not what you want, please specify the type explicitly in the data command. sherpa> instrument 10 = rsp[megp1](,1318MEG_1_garf.fits,) The inferred file type is ARF. If this is not what you want, please specify the type explicitly in the data command. sherpa> # turn off Y errorbars sherpa> sherpa.dataplot.y_errorbars=0 sherpa> sherpa.dataplot.curvestyle="histo" sherpa> sherpa.dataplot.symbolstyle="none" sherpa> lp 4 data 3 data 4 data 9 data 10 sherpa> ignore allsets all sherpa> notice allsets wave 8.2:8.6 sherpa> lp 4 data 3 data 4 data 9 data 10
For this thread, we are not using a grating RMF (which acts as a line shape function for grating data), and are therefore assuming that the line profile is gaussian.
Figure 1 shows the resulting plot which highlights a feature that is present in all four orders of the observation.
Defining the Source Model
We model this source with a 1-D normalized Gaussian (ngauss1d) combined with a 1-D polynomial function (polynom1d). Separate source models are used for the HEG and MEG datasets.
sherpa> source 3,4 = ngauss[hg1] + poly[hp1] sherpa> source 9,10 = ngauss[mg1] + poly[mp1] sherpa> mg1.ampl => hg1.ampl sherpa> show source 10 Source 10: (mg1 + mp1) ngauss1d[mg1] (integrate: on) Param Type Value Min Max Units ----- ---- ----- --- --- ----- 1 fwhm thawed 9.9987e-02 9.9987e-04 9.9987 2 pos thawed 8.4225 8.1975 8.5975 3 ampl link 1.7195e-03 expression: hg1.ampl poly1d[mp1] (integrate: on) Param Type Value Min Max Units ----- ---- ----- --- --- ----- 1 c0 thawed 4.681e-03 -1.022e-04 9.2599e-03 2 c1 frozen 0 -2.2894 2.2894 3 c2 frozen 0 -5.7236 5.7236 4 c3 frozen 0 -1.022e-04 9.2599e-03 5 c4 frozen 0 -1.022e-04 9.2599e-03 6 c5 frozen 0 -1.022e-04 9.2599e-03 7 c6 frozen 0 -1.022e-04 9.2599e-03 8 c7 frozen 0 -1.022e-04 9.2599e-03 9 c8 frozen 0 -1.022e-04 9.2599e-03 10 offset frozen 0 -8.1975 8.5975
We choose to link the amplitudes for the models (mg1.ampl => hg1.ampl). This forces Sherpa to find the best-fit amplitude for all four datasets.
Fitting
Fit all four datasets simultaneously:
sherpa> fit 3,4,9,10 WARNING (Sets 3,4,9,10): background data have been entered, but they have not been subtracted, nor have background models been set. LVMQT: V2.0 LVMQT: initial statistic value = 118084 LVMQT: final statistic value = 138.672 at iteration 14 hg1.fwhm 0.0118774 hg1.pos 8.42135 hg1.ampl 0.000177587 hp1.c0 0.000341403 mg1.fwhm 0.0176838 mg1.pos 8.42178 mp1.c0 0.000421502
The fit gives us a line position and flux (for a normalized gaussian, the flux is simply the amplitude [hg1.ampl]: 1.81e-4 photons/cm2/s). The mg1.ampl is not listed because it was linked to hg1.ampl, and so has the same best-fit value.
Identify the Line
The next step is to use the identify command to determine the source of the line; it takes a given wavelength and searches the APEC line list for strong lines with wavelengths close to it (within 0.01 Å by default). An optional second parameter allows a search range to be set (i.e. identify(8.42, 0.02) prints all lines with wavelengths between 8.40 and 8.44 Å):
sherpa> identify(8.42, 0.02) Lambda -- Ion UL - LL Emissivity@ kT RelInt For More Info Angstrom ph cm^3/s keV 8.4053 Fe XXII 177- 8 1.28e-18 @ 1.085 0.019 describe(26,22,177,8) 8.4192 Mg XII 4- 1 6.89e-17 @ 0.862 1.000 describe(12,12,4,1) 8.4246 Mg XII 3- 1 3.45e-17 @ 0.862 0.500 describe(12,12,3,1)
The fifth and sixth columns give the peak emissivity and temperature, respectively. The best identification is usually the strongest line in the list; if the peak emissivities are similar, the line could be a blend. In this case, the strongest lines are the 4->1 and 3->1 transitions of Mg XII (hydrogen-like magnesium).
We can find out more about some of the lines with the describe command; this command gets its information from APED. The syntax is describe(element, ion, upperlevel, lowerlevel): element is the number of protons (i.e. Mg = 12), ion is the ion stage in astronomical usage (i.e. XII = 12), and the upper and lower energy levels are given in the identify list. Note that the appropriate describe syntax is provided by the identify command. To find out about the strongest line in the previous list:
sherpa> describe(12,12,4,1) Ion Mg XII, energy level 1 --- electron configuration : 1s~^2S_{1/2} energy above ground (eV) : 0.000000 Quantum state : n=1, l=N/A, s=2, degeneracy=2 Energy level data source : 1983ADNDT..29..467S Photoionization data source : 1964ApJS....9..185B ------------------------------------------- Ion Mg XII, energy level 4 --- electron configuration : 2p~^2P_{3/2} energy above ground (eV) : 1469.430054 Quantum state : n=2, l=1, s=2, degeneracy=4 Energy level data source : 1983ADNDT..29..467S Photoionization data source : 1964ApJS....9..185B ------------------------------------------- Ion Mg XII, 1 - 4 interactions --- Electron collision rate from 1 -> 4 : nonzero. Reference bibcode : 1983ADNDT..29..467S Wavelength (lab/observed) (Angstrom) : 8.419209 +/- 0.000040 Wavelength (theory) (Angstrom) : 8.438330 Transition rate/Einstein A (s^-1) : 1.278220e+13 Wavelength (lab/observed) reference : 1977JPCRD...6...3E Wavelength (theory) reference : 1983ADNDT..29..467S Transition rate reference : 1987JPhB...20.6457F
This tells us that the 4->1 transition in Mg XII is in fact an n=2->1 hydrogen-like transition, or one component of the hydrogen-like Mg XII Lyman alpha line. Using describe(12,12,3,1) shows that it is the other transition in the n=2->1 doublet. This identification information (along with the current filter) can then be associated with a particular model element (in this case, the gaussian model hg1 used to fit the HEG +/-1 orders) using the lineid and filter commands:
sherpa> hg1 lineid "APECline(12,12,4,1)+APECline(12,12,3,1)" sherpa> hg1 filter "ignore allsets all; notice allsets wave 8.2:8.6"
Write an MDL File
Finally, the results are written to an MDL file which stores the data, the model, and the identification. This formatted FITS file can be read back into Sherpa (using read mdl "MgXII_MDL.fits") and can also be used for more sophisticated projects, such as fitting a differential emission measure (DEM) model.
sherpa> write mdl "MgXII_MDL.fits" WARNING (Sets 3,4,9,10): background data have been entered, but they have not been subtracted, nor have background models been set.
Use prism to examine the file that was just created:
sherpa> prism MgXII_MDL.fits
Figure 2 shows the resulting display; the modeling information is saved in the MDL_Models block.
History
14 Jan 2005 | reviewed for CIAO 3.2: no changes |
21 Dec 2005 | reviewed for CIAO 3.3: no changes |
09 Feb 2006 | minor change to filenames; organized thread into sections |
01 Dec 2006 | updated for CIAO 3.4: Sherpa version |