Isis tutorial gratings
Contents
Further Data Fitting
The walk through isis has given you a first overview of a simple Xray spectral analysis. In this part of the tutorial we build upon this experience by performing a more complicated analysis of Chandra gratings spectra. You will learn about Xray data archives and what to do when not all keywords are set in a Xray spectral file. Furthermore we introduce an important Xray diagnostics for tenuous plasmas.
Data Archives and the ADS
ADS
In this tutorial we will demonstrate the analysis of Chandra gratings data of stellar spectra and try to reproduce the results obtained by JanUwe Ness et al., 2002, Coronal density diagnostics with Heliumlike triplets: CHANDRALETGS observations of Algol, Capella, Procyon, epsilon Eri, alpha Cen A&B, UX Ari, AD Leo, YY Gem, and HR 1099, Astronomy and Astrophysics 394, 911.
Exercise 1
 Go to the ADS Abstract Service and search for the above paper. The best way to do this is to enter
^Ness
for the Author andCapella
for the Object. The ADS is the best data base to find astronomical literature. It goes back to the early 1800s and is almost complete for the past 50 years or so. It will be very worth your time to familiarize yourself with searching papers there  even if you do not want to do astronomy in the end but other areas of physics, which are indexed in the ADS as well. In order to make sure that you really did the search and did not just click the link above, what is the newest paper published by this author?
 Download, print, and read this paper. Yes, really, print it  you will need to refer to this paper often and you might want to scribble comments on your paper copy and you will read the paper much more intensively if it is not on a computer screen. In order to get to the PDFcopy of the paper, click on "Full Refereed Journal Article (PDF/Postscript)" in the above link. You will see it is colored in green, meaning that the paper is not behind a paywall. This is not the case for new articles, unfortunately, in this case you can always go to the arXivLink and get a preprint from there. Virtually all important astronomical papers can be found this way.
Data Archives
Here, we will show you how to perform the analysis of Capella, in your homework you will work on one of the other sources. Like the analysis of data from all other Xray satellites, the reduction of Chandra gratings data can become very involved. Since we "just" want to analyze the time averaged spectrum of the source, we are lucky that somebody else (well, the Chandra gratings team) has done the work for us and has put the reduced spectra on the WWW for us. Such archives of prereduced data exist for all satellites, and while I caution not to use the reduced data in these archives for publicationlevel work, they are useful to get a feel for what is available (the reason why you do not want to use these data for publication level work is that you will typically have to apply a newer version of the detector calibration, different background screening criteria, different time intervals, and so on). The following list contains the most important Xray astronomical archives:
 TO BE WRITTEN
 TGCat contains all public observations performed with the Chandra transmission gratings spectrometers. The user interface is selfexplanatory and there is a nice help section available. If you know the name of the object you're interested in, enter it in the "Quick Query" box, then click on the object name and you will be getting a list of all available observations.
Exercise 2
 Download the observation with the longest exposure time available for the object you are interested in, unless your instructions say otherwise. If possible, disregard the observations where the HRC was in the focal plane, i.e., look at the ACISobservations only (the HRC is not energy dispersive, and therefore order sorting is not so easy with this detector).
Select the data, ask TGCat to send you an email, and then download the file. For example
wget http://tgcat.mit.edu/tmp/1337349938_174950000/tgcat.tar.gz
Unpack it with
tar xzvf tgcat.tar.gz
For Capella, we download the data from observation ID 8319. The data set looks as follows:
hde226868:~/isistut> cd tgcat hde226868:~/isistut/tgcat> ls l total 16 drwxrxrx 3 wilms wilms 4096 May 18 16:05 . drwxrxrx 3 wilms wilms 4096 May 20 14:05 .. rwrr 1 wilms wilms 327 May 18 16:05 checksums.dat drwxrxrx 2 wilms wilms 4096 May 18 16:05 obs_8319_tgid_3283 hde226868:~/isistut/tgcat> cd obs_8319_tgid_3283 hde226868:~/isistut/tgcat/obs_8319_tgid_3283> ls l total 15648 rrr 1 wilms wilms 1419840 May 18 16:05 pha2.gz rrr 1 wilms wilms 256320 May 18 16:05 leg_1.arf rrr 1 wilms wilms 7030080 May 18 16:05 leg_1.rmf rrr 1 wilms wilms 256320 May 18 16:05 leg_1.arf rrr 1 wilms wilms 7032960 May 18 16:05 leg_1.rmf rwrr 1 wilms wilms 71 May 18 16:05 vv_report.txt
The spectral file is called pha2.gz
, which is not a very meaningful name, so we unpack it and rename it:
hde226868:~/isistut/tgcat/obs_8319_tgid_3283> gunzip pha2.gz hde226868:~/isistut/tgcat/obs_8319_tgid_3283> mv pha2 capella.pha2
Note the name pha2
. In the walkthrough we encountered PHAfiles, which contain one spectrum. pha2
files contain multiple spectra. They are used, for example, in pulse phase resolved spectroscopy of neutron stars or in gratings spectra, where they contain multiple spectra.
Exercise 3
 Take a look at the contents of your pha2file using
fv
andfstruct
. Note that each "row" in the pha2file now contains one spectrum.
Preparing Gratings Data for Fitting
We now enter "isis", load the "isisscripts", and load our data set
hde226868:~/isistut/tgcat/obs_8319_tgid_3283> isis : isis> require("isisscripts"); : isis> capella=load_data("capella.pha2"); Reading: ......
As usual, we first take a look at what has been loaded, this time using a helpful function from the "isisscripts" that lists the data and the loaded "ARF" and "RMF"files.:
isis> list_all; Current Spectrum List: id instrument part/m src use/nbins A R totcts exp(ksec) target 1 LETGACIS leg3 1 8192/ 8192   1.0222e+04 59.146 Capella file: capella.pha2 2 LETGACIS leg2 1 8192/ 8192   4.4290e+03 59.146 Capella file: capella.pha2 3 LETGACIS leg1 1 8192/ 8192   7.0583e+04 59.146 Capella file: capella.pha2 4 LETGACIS leg+1 1 8192/ 8192   1.1496e+05 59.146 Capella file: capella.pha2 5 LETGACIS leg+2 1 8192/ 8192   4.7930e+03 59.146 Capella file: capella.pha2 6 LETGACIS leg+3 1 8192/ 8192   6.2380e+03 59.146 Capella file: capella.pha2 isis>
We note two things:

load_data
loaded all contents of thepha2
file, total of six spectra. What these spectra are is identified in thepart/m
column:leg
is the "low energy grating" of ChandraHETGS (the other grating would be calledmeg
formedium energy grating
), and 3, 2, 1, +1, +2, +3 identifies the order of the spectrum.  No
ARF
andRMF
are listed. This means that while the spectra have been loaded we still need to tellisis
whatARF
andRMF
to use for the spectra (you can also see this by noticing the dashes in theA
andR
columns).
We now remind ourselves what responses and ARF we were provided with, by using isis' "shell escape". This means that within isis we can precede a normal shell command by an exclamation mark to have it executed by the shell:
isis> !ls l total 15648 rrr 1 wilms wilms 1419840 May 18 16:05 capella.pha2 rrr 1 wilms wilms 256320 May 18 16:05 leg_1.arf rrr 1 wilms wilms 7030080 May 18 16:05 leg_1.rmf rrr 1 wilms wilms 256320 May 18 16:05 leg_1.arf rrr 1 wilms wilms 7032960 May 18 16:05 leg_1.rmf rwrr 1 wilms wilms 71 May 18 16:05 vv_report.txt isis>
This means that we only have RMFs and ARFs available for the first order and we will have to ignore the higher orders. This is fine since most photons end up in the 1st order anyway. Let us now load the files:
isis> arfm1=load_arf("leg_1.arf"); isis> arfp1=load_arf("leg_1.arf"); isis> rmfm1=load_rmf("leg_1.rmf"); isis> rmfp1=load_rmf("leg_1.rmf");
and check that everything went fine:
isis> list_all; Current Spectrum List: id instrument part/m src use/nbins A R totcts exp(ksec) target 1 LETGACIS leg3 1 8192/ 8192   1.0222e+04 59.146 Capella file: capella.pha2 2 LETGACIS leg2 1 8192/ 8192   4.4290e+03 59.146 Capella file: capella.pha2 3 LETGACIS leg1 1 8192/ 8192   7.0583e+04 59.146 Capella file: capella.pha2 4 LETGACIS leg+1 1 8192/ 8192   1.1496e+05 59.146 Capella file: capella.pha2 5 LETGACIS leg+2 1 8192/ 8192   4.7930e+03 59.146 Capella file: capella.pha2 6 LETGACIS leg+3 1 8192/ 8192   6.2380e+03 59.146 Capella file: capella.pha2 Current RMF List: id grating detector m type file/function 1 LETG ACIS7 1 file: leg_1.rmf 2 LETG ACIS7 1 file: leg_1.rmf Current ARF List: id grating detector part/m src nbins exp(ksec) target 1 LETG ACIS leg1 0 8192 59.15 Capella file: leg_1.arf 2 LETG ACIS leg+1 0 8192 59.15 Capella file: leg_1.arf isis>
We now need to assign each ARF and RMF to the spectrum taken with this detector:
isis> assign_rmf(rmfm1,capella[2]); isis> assign_arf(arfm1,capella[2]); isis> assign_rmf(rmfp1,capella[3]); isis> assign_arf(arfp1,capella[3]);
Note that instead of using the id
values directly (e.g., assign_rmf(1,3)
) we use the isisvariables that were returned in the different load
commands (capella
is a variable, and in slang/isis arrays are zerobased). This approach is better than using the ids directly since it minimizes the likelihood of typos and we do not have to remember strange id numbers but rather variable names that (should) make sense. If you want to find out more about the slang programming language, you can have a look at the nice documentation at [1]. In this case, especially Chapter 6 (Variables) and 11 (Arrays) might be interesting for you.
Let's check whether the assignment was ok:
isis> list_data; Current Spectrum List: id instrument part/m src use/nbins A R totcts exp(ksec) target 1 LETGACIS leg3 1 8192/ 8192   1.0222e+04 59.146 Capella file: capella.pha2 2 LETGACIS leg2 1 8192/ 8192   4.4290e+03 59.146 Capella file: capella.pha2 3 LETGACIS leg1 1 8192/ 8192 1 1 7.0583e+04 59.146 Capella file: capella.pha2 4 LETGACIS leg+1 1 8192/ 8192 2 2 1.1496e+05 59.146 Capella file: capella.pha2 5 LETGACIS leg+2 1 8192/ 8192   4.7930e+03 59.146 Capella file: capella.pha2 6 LETGACIS leg+3 1 8192/ 8192   6.2380e+03 59.146 Capella file: capella.pha2
Notice that the A
and R
columns now contain the indices of the respective ARF and RMF.
The First Fit
Since we do not have an ARF and RMF for the other orders, we need to tell isis not to work with these data sets:
isis> exclude([capella[0],capella[1],capella[4],capella[5]);
The argument of the exclude
function is an array containing the indices of all spectra that should be ignored when fitting. In slang/isis, arrays are defined by listing their contents in brackets, i.e.,
array=[1,2,5,6];
It is also possible to use a shorthand notation to define larger arrays, e.g.,
array=[1:10];
creates an array containing the 10 numbers 1,2,...,10, while
array=[1:10:2];
creates an array containing the numbers 1,3,5,...,9. In principle, you can do something like this also with floating point numbers, but because of subtleties with floating point arithmetic, we do not recommend using this feature of the language.
In the exclude
example before, note that we are again listing the indices of the spectra indirectly by specifying them through the capella
variable. As before, this makes our code more portable. Finally, note that you could have also used the ignore
function to ignore these data sets, but exclude
is preferred since it does not affect the currently noticed energy ranges.
Finally, since we now only have two spectra to work with, let's remember their indices in shorter variables such that we do not have to type capella[x]
all the time...
isis> m1=capella[2]; isis> p1=capella[3];
We first plot the whole data set, switching the xaxis unit to Angstroms, because this is more appropriate for gratings data:
isis> fancy_plot_unit("A"); isis> plot_data(p1;dsym=1,dcol=2);
Note that for the moment we are only plotting one data set in order not to confuse ourselves. With dsym=1
we switch off the plot symbol, which is only distracting and does not help in understanding the figure (try the command without the dsym
qualifier!).
The resulting figure looks as follows:
Notice the very strong emission lines, which are typical for stellar coronae. Most of the things are happening between 5 and 25 Angstroms, so let's zoom in that region:
Isis has very powerful functions that help in the identification of these lines. Also note the "bump" at lower wavelengths, which is due to bremsstrahlung emission of the hot corona.
Right now we are only interested in the Helike oxygen triplet at 21.6, 21.8, and 22.1 Angstroms. These are called the resonant, intercombination, and forbidden line, respectively, and their fluxes are typically called r, i, and f. Let's therefore take a look at the 2023 A band only:
All three lines can be clearly seen in the data set. In order to get more quantitative information, a look at the spectral plot indicates that we should be able to model this spectral region well with a spectral model that consists of the three lines plus a constant flux component (model name constant
) that describes the continuum emission. We should probably also ignore the 20.521.5 A band because otherwise the fit would be contaminated by the emission lines there. Finally, experience tells us that in such line fits it is often not a good idea to let the wavelength/energy of the line to be completely free. The reason is that often initially the minimization routine will be trying to change the wavelength by a large factor and as a result it could be that the lines are switched. We therefore allow the lines only to vary by 0.1 A around their official values.
With the commands discussed in the walkthrough, and with the information given by help set_par
you should now be able to setup a fit function that looks like this (note: the gauss
model in isis is a Gaussian line in Angstrom units, the gaussian
model is the same, but in keV units):
isis> list_par; constant(1)+gauss(1)+gauss(2)+gauss(3) idx param tieto freeze value min max 1 constant(1).factor 0 0 1 0 1e+10 2 gauss(1).area 0 0 1 0 0 photons/s/cm^2 3 gauss(1).center 0 0 21.6 21.5 21.7 A 4 gauss(1).sigma 0 0 0.025 1e06 1 A 5 gauss(2).area 0 0 1 0 0 photons/s/cm^2 6 gauss(2).center 0 0 21.8 21.7 21.9 A 7 gauss(2).sigma 0 0 0.025 1e06 1 A 8 gauss(3).area 0 0 1 0 0 photons/s/cm^2 9 gauss(3).center 0 0 22.1 22 22.2 A 10 gauss(3).sigma 0 0 0.025 1e06 1 A
To get a good fit, plot the model (do not forget the eval_counts
command!), and adjust the constant, line widths, and line fluxes (area
parameter) until you get a good "fit by eye". Playing a little bit with the parameters I find a starting model that looks as follows:
We now ignore the 20.621.5 and 22.622.8 Angstrom bands in both data sets:
isis> ignore([p1,m1],20.6,21.5); isis> ignore([p1,m1],22.6,22.8);
Replot in order to confirm that these data are gone! Furthermore, we also ignore everything below 20A and above 23A:
isis> ignore([p1,m1],NULL,20); isis> ignore([p1,m1],23,NULL);
perform the fit:
isis> Fit_Verbose=1; isis> fit_counts();
The best fit parameters are as follows:
isis> list_par; constant(1)+gauss(1)+gauss(2)+gauss(3) idx param tieto freeze value min max 1 constant(1).factor 0 0 1.342954e05 0 1e+10 2 gauss(1).area 0 0 0.001175245 0 0 photons/s/cm^2 3 gauss(1).center 0 0 21.60025 21.5 21.7 A 4 gauss(1).sigma 0 0 0.01614977 1e06 1 A 5 gauss(2).area 0 0 0.000183689 0 0 photons/s/cm^2 6 gauss(2).center 0 0 21.79914 21.7 21.9 A 7 gauss(2).sigma 0 0 0.0108202 1e06 1 A 8 gauss(3).area 0 0 0.0008692617 0 0 photons/s/cm^2 9 gauss(3).center 0 0 22.09971 22 22.2 A 10 gauss(3).sigma 0 0 0.02135323 1e06 1 A
To check whether all went well with the fit, we plot our best fit. While for the manual adjustment of the fit parameters it was sufficient to just look at one of the spectra, we really should look at all of our data. This is done with
isis> plot_data({m1,p1};dsym={1,1},decol={2,4},dcol={2,4},mcol={4,4},res=1);
where the {1,2}
etc. are socalled lists. These are somewhat more flexible forms of arrays. We will address the details of this data type later, for the moment, just remember that you need to put everything in curly parentheses instead of brackets. In the above, the color 4 is used for the m=+1 spectrum and 2 for the m=1 spectrum.
The best fit looks as follows:
and the residuals are flat. This is not surprising, given that the reduced chisquare of the fit was only 0.5.
Looking at the paper by Ness et al., for the spectral diagnosis the ratio as well as the ratio is used. As discussed by Ness et al., the diagnostic information contained in these ratios is discussed in great detail by Porquet et al. (2001) (where the resonance line is called w, the intercombination x+y, and the forbidden line z; sorry, nobody claimed scientists are consistent...). Porquet et al. (2001) show that R is mainly dependent on the electron density and G is mainly dependent on the temperature.
The following analysis can now proceed either by using your pocket calculator or by using isis as a pocket calculator. Since we have not yet introduced programming in slang/isis, feel free to use either way, or [isis:tutorial:slang take a look at our Slang programming introduction]. In the following, we show how one would get about calculating the above ratios with isis. To understand the code, you need to know that in slang structures contain several named values and that you can access the one of the variables ("tags") that make up the content of a structure with the usual variable.tag
syntax that is also used, e.g., in C. Furthermore, read the help for the isis get_params
function^{[1]}.
To calculate the ratios then do the following:
isis> r=get_params("gauss(1).area"); isis> i=get_params("gauss(2).area"); isis> f=get_params("gauss(3).area");
take a look at the contents of r (which per the documentation of get_params
is an array):
isis> print(r); {name="gauss(1).area", index=2, value=0.0011752452373316315, min=1.7976931348623157e+308, max=1.7976931348623157e+308, hard_min=inf, hard_max=inf, step=0.0, relstep=0.0001, freeze=0, tie=NULL, units="photons/s/cm^2", is_a_norm=1, fun=NULL}
and therefore
isis> rat=f[0].value/i[0].value; isis> g=(f[0].value+i[0].value)/r[0].value; isis> print(rat); 4.732248259283046 isis> print(g); 0.8959412576621245 isis>
Per the discussion in Porquet, the small value of G is an indication that the observed plasma is dominated by collisions, i.e., it is very thin and very hot and the ionization of the plasma is due to collisions and not due to a strong Xray source in the vicinity that ionizes the plasma through photoionization. We are also happy that our ratio for R, 4.7, is close to that found by Ness et al., which is around 4 (right hand panel of their Figure 10).
While this is all very nice and the result that we obtain is in rough agreement with the literature, we do not know how uncertain our parameters are. The reason is that even though we have the best fit, we do not yet know the uncertainty of our fit parameters. To calculate the uncertainty of a parameter isis uses the usual approach, first discussed for Xray astronomical data by Lampton, Margon, Bowyer, 1978, Apj 208, 177. We will call the parameter for which the uncertainty is to be calculated the "parameter of interest". Lampton et al. show that the confidence interval can be determined with the following recipe:
 Calculate from the chisquaredistribution with one degree of freedom. Common in Xray astronomy are ""errors with , 90% confidence ranges with , and 99% confidence contours with . The most common errors are the 90% ones.
 Change the parameter of interest away from its best fit value to a slightly smaller value. Keep it frozen and refit. This will result in a slightly larger for the best fit.
 Repeat the previous step until you find the parameter value where the chisquared is where is the chisquared of our original best fit. The value of our parameter is now the lower bound of the confidence interval.
 Repeat for values larger than the best fit value until you find the upper bound of the confidence region.
 Designating the lower and upper end of the confidence region by and and the best fit value , we then write the best fit parameter as follows: where and .
Because of correlations between parameters, error bars are typically not symmetric and great care has to be taken when doing error propagation. We ignore all of those for the moment and use the rule of thumb that if the error bars are not too asymmetric, a rough estimate for a symmetric error bar which can be used for the standard Gaussian error propagation is the average of the lower and the upper error bars, i.e., use and write the best fit parameter as (where we usually also use 90% error bars!).
The above procedure is implemented in isis through the conf
and vconf
functions. These are identical, vconf
gives more information about the procedure (the "v" stands for "verbose"). Since the above procedure searches for the zero of a function, it is not guaranteed to converge numerically. For this reason isis calculates confidence intervals only when explicit limits are set. These limits are returned in case the error calculation failed. Typically, for good data sets the errors will be less than 2030% of the value, but this depends a lot on the data set and one often has to play around a little bit.
The error bar for the flux of the forbidden line is then calculated as follows:
isis> set_par("gauss(3).area",f[0].value,0,0.75*f[0].value,1.25*f[0].value); isis> list_par; constant(1)+gauss(1)+gauss(2)+gauss(3) idx param tieto freeze value min max 1 constant(1).factor 0 0 1.342954e05 0 1e+10 2 gauss(1).area 0 0 0.001175245 0 0 photons/s/cm^2 3 gauss(1).center 0 0 21.60025 21.5 21.7 A 4 gauss(1).sigma 0 0 0.01614977 1e06 1 A 5 gauss(2).area 0 0 0.000183689 0 0 photons/s/cm^2 6 gauss(2).center 0 0 21.79914 21.7 21.9 A 7 gauss(2).sigma 0 0 0.0108202 1e06 1 A 8 gauss(3).area 0 0 0.0008692617 0.0006519463 0.001086577 photons/s/cm^2 9 gauss(3).center 0 0 22.09971 22 22.2 A 10 gauss(3).sigma 0 0 0.02135323 1e06 1 A isis> (fmi,fma)=vconf("gauss(3).area"); : : lots of output... : gauss(3).sigma= 2.25766046e02 chisquare value converged par[8]= 9.64711324e04 stat= 1.5337e+02 dstat= 2.7100e+00 limit found: 9.64711324e04 dstat= 2.7100e+00 isis> print(fmi); 0.0007743270669346407 isis> print(fma); 0.0009647113237380715 isis> print(f[0].valuefmi); 9.493467346016499e05 isis> print(fmaf[0].value); 9.544958334326589e05 isis> print(f[0].value); 0.0008692617403948057
this means that the error bar is almost symmetric and therefore
Exercise 4
 Perform the above analysis for "your" star. Do not forget to calculate error bars for R and G using error propagation!
 Use the tables of Porquet et al. to estimate the density and temperature of the plasma. Looking at the discussion by Ness et al., what is the systematic error of this analysis?
References
 ↑ if one is just interested in the best fit value of a parameter, you can also use the
get_par
function, however, this would not have allowed us to introduce slang structures...