RHESSI Spectroscopy -- Second Steps

RHESSI Spectroscopy -- Second Steps

Linhui Sui, Astrid Veronig, Brian Dennis, Richard Schwartz

Last modified 2 February 2004

List of topics*:

1. How to analyze a flare with different attenuator states.

2. How to select multiple analysis intervals.

3. How to change fitting models.

4. How to save and restore the fitting parameters and spectral data.

5. How to plot the fitted spectra outside SPEX.

6. How to make a movie of the spectra.

7. How to choose and evaluate a background model.

8. How to create and customize postscript plots within SPEX.

9. How to use RHESSI data of a different orbit for Background subtraction.

*: Throughout this tutorial, the X4.8 flare on July 23 2002 and the M1.5 flare on June 1 2002 are used as examples to explain the topics. Words in italics represent SPEX commands

1. How to analyze a flare with different attenuator states

During most strong flares, the attenuator state will be changed. The following plot is an uncorrected quicklook light curve of the X4.8 flare on July 23 2003 plotted in the GUI.

As indicated by the purple horizontal bars, before 00:26:00 UT, the attenuator state was A1 (thin in). From 00:26:00 UT to 00:59:24 UT, the attenuator state was mostly A3 (both thick and thin in) with four short periods in the A1 state. After 00:59:24 UT, the attenuator state returned to A1 state.

To analyze count flux data in SPEX, a count rate spectrum file and a spectrometer response matrix (SRM) file must be read into the active session of SPEX. The response matrix is stored in the array DRM (for detector response matrix) within the SPEX common blocks. When using the RHESSI spectrum object or GUI to create the count rate spectrum file and response matrix file, the name of the response matrix file is written in the RESPFILE parameter in the header of the spectrum FITS file. A different srm file is needed to analyze the count rate data obtained in each attenuator state and detector combination. The srm file and count spectra files are written out from the Spectrum widget of the GUI as explained in Spectroscopy – First Steps. Alternatively, the srm files can be produced from the spectrum object (sp_obj->filewrite, /fits,/build,... ). They are appropriate for the attenuator state in the first interval after the start time. If the attenuator state changes during the flare, then a new srm file must be generated for each state by choosing the first time interval with the new attenuator state.

The srm file is specified in SPEX with dfile, and then use read_drm to create the DRM array used with SPEX.

For this 23 July 2002 X4.8 event, one spectrum file and two srmfiles were saved with the GUI. They are:

hsi_spectrum_20020723_001000.fits -- spectrum file

hsi_srm_20020723_001000_atten1.fits -- srm file for A1 state

hsi_srm_20020723_002600_atten3.fits -- srm file for A3 state

If we want to analyze the spectra before 00:26:00 UT, we have to first read the spectrum file, and then read the srm file appropriate for the A1 state.

SPEX> _1file, hsi_spectrum_20020723_001000.fits

SPEX> preview

SPEX> graph

SPEX> dfile,hsi_srm_20020723_001000_atten1.fits

SPEX> read_drm

This should give the following response if everything is OK:

Command: read_drm

MRDFITS: Null image, NAXIS=0

MRDFITS: Binary table. 6 columns by 92 rows.

MRDFITS: Binary table. 3 columns by 1 rows.

minmax(drm) 0.000000 0.0738023

Change parameters or enter a command

If we want to analyze the spectra between 00:26:00 UT and 00:59:24 UT (for now, because of the pileup problem we will not try to analyze the time intervals during the four short periods in the A1 state), we have to read the srm file appropriate for the A3 state:

SPEX> dfile,hsi_srm_20020723_002600_atten3.fits

SPEX> read_drm

It will give a similar output as above. Again, for the analysis of spectra after 00:59:24 UT, we need to switch back to the A1 srm file:

SPEX> dfile,hsi_srm_20020723_001000_atten1.fits

SPEX> read_drm

Once the correct spectrum and srm files have been entered into SPEX, the analysis can proceed as described in Spectroscopy –First Steps.

2. How to select multiple analysis intervals

SPEX is capable to doing spectral fitting for multiple time intervals. In order to do this, we need to first select multiple analysis intervals. Since the July 23 event lasts for a very long time (about 1 hour), it is very difficult to select time intervals in the crowded time history window, even though the rates are binned every 4 s. First, we can enlarge the time history widow by clicking on one edge and dragging it to produce the desired size. The lightcurve with background subtracted then appears as follows:

An enlarged version of the time of interest can be obtained with the zoom command. First type zoom, then select the “graphic input”, “continue by selecting from the window”, then select the start time and end time of the zoom-in region. A zoomed in time history will appear as follows:

It is much more convenient to select time intervals from the window above to analyze the main impulsive peak of the flare (from 00:26:20 UT to 00:39:40 UT). By typing clear_ut, then typing graph in SPEX, the time history window will return to the original time history plot.

After you type select, the following window will pop up:

This window lists three ways to select time intervals.

Method 1: Boundaries for contiguous intervals. This means the end of time one interval is the start time of the following interval. To use this method, click on the top button in the above window, click “continue by selecting from the window” in the next widget, then select the start of the first interval by moving the curser to the desired time in the Time History Window and clicking, and then select the end of the first , second, third … intervals until you have selected all the intervals. To show what you have selected, type the following in SPEX:

SPEX>graph

SPEX>display

The following selected time history window will appear:

We selected four contiguous intervals. Variable ifirst and ilast are the start and end time interval number to fit. The default of ifirst is 0, ilast is the last interval, which is 3 for this case.

SPEX> print,ifirst, ilast

0 3

If we type count to see the count flux spectrum, the spectrum is for the first time interval (00:27:27.999 -- 00:28:31.999 UT):

Using what we have learned from Spectroscopy - First Steps, we set the function we want to fit (f_vth_bpow model is the default), the energy range to fit, initial estimates of the fitting parameters, and which ones should be free to vary. We then type fit to fit the spectra for all the intervals sequentially.. As the fitting proceeds, the spectral window will show the fitted spectrum, one by one, and stop at the final time interval. All the fitting parameters are saved in variable array, apar_arr:

SPEX> print, apar_arr

1.24413 3.35509 14.1803 2.93961 81.1886 3.32381

2.68748 3.31781 14.0138 2.70112 78.9045 3.14351

3.57966 3.21793 13.6295 2.60562 81.8491 3.09796

3.91184 3.17164 8.94853 2.71290 76.7049 2.99297

The reduced chi-square values are stored in the array chi:

SPEX> print, chi

0.642303 1.05372 0.780159 1.05493

If we only want to fit the spectrum for interval 0, we need change the value of ilast to 0

SPEX> ilast, 0

SPEX>print, ifirst, ilast

0 0

SPEX>count

SPEX>fit

The following plot of the best fit f_vth_bpow photon flux spectrum appears along with the residuals expressed as the number of sigma above or below the model values:

We can set eplot to 1 so that we can see the error bar in spectral data:

SPEX> eplot, 1

SPEX> photon

If we want to fit the spectrum at the second time interval, we need to reset ifirst and ilast:

SPEX>ifirst, 1

SPEX>ilast, 1

Then we can do the fit for the second interval. Similarly, we can do the third and forth interval and so on.

Method 2: Boundaries for discrete intervals. Since the intervals are discrete, we need to select the start and end time of each interval. Click on the middle button, click “continue by selecting from the window” in the next widget, then select the start and the end of each time interval by moving the curser to the appropriate times in the Time History Window and clicking with any mouse button. As many intervals can be selected as needed but the end of one interval will not necessarily be the start of the next. To show what we have selected, we should type the following in SPEX:

SPEX>graph

SPEX>display

The selected discrete intervals are shown in the following:

The vertical orange dotted/broken lines mark the start/end of each interval, respectively. To get the time in seconds for each time interval:

SPEX> print, xselect

1663.9990 1692.0000

1727.9990 1743.9990

1760.0000 1779.9990

1824.0000 1851.9990

Type fit to fit all intervals one by one. If we want to fit a specific interval, set the ifirst and ilast to that interval. These steps are similar to those in Method 1 above.

Method 3: As limits on sample groups of N (to be chosen). This method allows us to choose one long interval and divide it up into equal intervals of length N times the basic binning time of the data. In our example, we have used 4-s binning times Thus, if we want to fit spectra in every 4-s interval, we would set N to 1. If wanted to use 8-s intervals we would set N to 2. Here, we set N to 5, so we are fitting the spectrum contiguously with 20-s time intervals.

Click on the button marked AS LIMITS ON SAMPLE GROUPS OF N (TO BE CHOSEN). The following window pops up:

Select 5 in the righthand column, click the “ready” button, then click “continue by selecting from the window.” Finally in the Time History Window, click the start and the end of time of the period in which you want to make your 20-s intervals. In order to see the intervals we selected, type the following:

SPEX>graph

SPEX>display

And this window appears:

We can see evenly selected contiguous 20-s intervals. By checking the size of the “xselect” array, we know how many intervals we have selected:

SPEX> help, xselect

XSELECT (SPEX_PROC_COM)

DOUBLE = Array[2, 37]

We selected 37 time intervals. We can fit the 37 spectra sequentially by assuring that ifirst = 0 and ilast = 36 and typing fit. After it finishes, we can check the best-fit parameters:

SPEX> fit

SPEX> help,apar_arr

APAR_ARR (SPEX_PROC_COM)

FLOAT = Array[6, 37]

SPEX>print, apar_arr[*, 0:3]

0.756227 3.20582 9.46661 3.31592 75.3729 3.83593

1.09641 3.35124 15.3057 2.90970 86.1484 3.37075

1.60972 3.38625 16.0934 2.82412 85.3917 3.15549

2.13309 3.36383 13.8587 2.75644 67.5149 3.19990

SPEX>print, chi[0:3]

0.928469 1.00548 0.788376 0.934216

The window with the fitted spectra and residuals cycles through very rapidly in the Spectral Window. However, we can use create_ps to save all the spectral and residual plots as postscript files

SPEX>psstyle =’fullportrait’

SPEX>create_ps

Then postscript files named like “fitting_100212_110946.ps” will be saved in your IDL working directory. It has 37 pages of spectral plots. Each one is one interval we select. This create_ps can also be used with the former two methods. For more detailed explanation see Topic 8.

If the postscript files we saved are not appropriate for publication purposes, or we need plots of the time history of some fitting parameters, then we have to save all the parameters and use them outside SPEX to make publication-quality plots.

After selecting the multiple intervals with either of the three methods above , we can check the variable “xselect” to see what the start and end time are for each time interval. Here, let us check the multiple time intervals we selected in method 1. The variable "xselect" will give us the time in units of seconds:

SPEX> help,xselect

XSELECT (SPEX_PROC_COM)

DOUBLE = Array[2, 4]

SPEX>print, xselect

1647.9990 1711.9990

1711.9990 1747.9990

1747.9990 1784.0000

1784.0000 1819.9990

If we want to know the "xselect" in units of UT, then we need to get UT base time (in unit of seconds) by issuing an IDL command:

SPEX> idl, ut_base = getutbase()

We check the ut_base in UT time, we type:

SPEX> idl, ptim, ut_base

23-Jul-2002 00:00:00.000

Then we print xselect in UT time format:

SPEX> idl, ptim, xselect + base

2002/07/23 00:27:27.999 2002/07/23 00:28:31.999

2002/07/23 00:28:31.999 2002/07/23 00:29:07.999

2002/07/23 00:29:07.999 2002/07/23 00:29:44.000

2002/07/23 00:29:44.000 2002/07/23 00:30:19.999

3. How to change fitting models

As explained in Spectroscopy- First Steps, there are many spectrum models available in SPEX. The default is the photon flux model named "f_vth_bpow". It is a combination of an isothermal function plus a double power-law. Let us try another model named "f_vth_thin".

Model "f_vth_thin" assumes a thin-target bremsstrahlung spectrum. By fitting with this model, we are able to get several parameters of the electron ensemble that produced the observed X-ray photon flux assuming thin-target interactions. This is in contrast to the parameters of the measured photon spectrum obtained using the "f_vth_bpow" model.

First, let us change the model to "f_vth_thin" to fit one interval (#8) of the above 37 intervals we selected.

SPEX> ifirst, 8

SPEX> ilast, 8

Current interval is set to one interval at 00:30:00 -- 00:30:19.999 UT. Then we change the fitting model to "f_vth_thin":

SPEX> f_model, f_vth_thin

SPEX> photon

Then we will need to input the initial gauss of the 8 parameters. Here are the brief explanation of those 8 parameters:

a(0)= emission measure in units of 10^49 cm-3

a(1)= KT plasma temperature in keV

a(2) - normalization factor (in 1.0d55), i.e. plasma density * nonthermal electron density * volume of source * electron velocity

a(3) - Power-law index of the electron distribution function below break energy.

a(4) - Break energy in the electron distribution function (in keV)

a(5) - Power-law index of the electron distribution function above break energy

a(6) - Low energy cutoff in the electron distribution function (in keV)

a(7) - High energy cutoff in the electron distribution function (in keV)

Since the computation of this model involves an integration, the fitting speed is much slow than that of "f_vth_bpow" model. One thing we can do it to decrease the high energy cutoff, i.e. a(7). The default is 32000 keV, for most spectrum (fitting energy range below 300 keV), we can cut down to 5000 keV without affecting the high energy flux.

As a rule of thumb, we should always obtain good starting parameters before we do any fitting. After changing the parameters, you can check the spectral window to compare the model spectrum and the data points. If there is not good agreement, the model parameters should be changed by typing the photon command. When the starting parameters give good agreement with the data points, then you can expect the fitting program to converge to good best-fit values. This is achieved with the fit command.

After we set fitting range and give a good initial values for those fitting parameters, we try a model with isothermal plus broken power law electron distribution:

SPEX>free, 1, 1, 1, 1, 1, 1, 0, 0

SPEX>fit

The following plot is produced for interval 8 in the spectrum window:

We can check the fitting parameters as follows

SPEX>print, apar

4.25845 3.19432 17.1155 1.58129 131.000 2.47378 10.0000 5000.00

SPEX>print, chi

1.09850

There is another model called "f_vth_thick" that fits a thick-target bremsstrahlung spectrum. It also has 8 parameters with the same physical meaning as those of "f_vth_thin" except for the normalization factor, i.e. a(2). For the thick-target model, a(2) is electron flux, i.e. the nonthermal electron density * source area * electron velocity (in units of 1.0e35 electrons s^-1). Since "f_vth_thick" involves two integrations, the fitting speed is much slower than for the other functions. It takes a few minutes to fit one interval. The fitting time depends on how many free parameter you set and the value of the high energy cut-off – a(7). This should be at least a factor of 10 greater than the highest photon energy to be fit.

Let us use "f_vth_thick" to fit the same time interval as above:

SPEX>f_model, f_vth_thick

SPEX>erange, 15, 300 ;set fitting range

SPEX>free, 1, 1, 1, 1, 1, 1, 1, 1, 0 ;model with isothermal plus broken power law with free low energy cutoff

SPEX>fit

The following plot is the spectral fit with 'f_vth_thick' we got:

SPEX> print,chi

        1.11807

SPEX> print, apar

        4.16548 3.21193 0.871608 3.67201 139.870 4.50256 37.9659 5000.00

After fitting the spectrum, we can compute the electron distribution using the "brm_distrn.pro" in SSW. Here are the details.

4. How to save and restore the fitting parameters and spectral data

If you want to save your actual SPEX session, simply type

SPEX>save_event

This command creates the files "summary.dat" and "spex_dflts.dat" in your IDL working directory. If you call SPEX again, you can continue where you left off by restoring the data with

SPEX>restore_event, /verb

After you restore the event, you may need to type read_drm in order to get correct RSM.

SPEX allows you to export all the resulting parameters to an IDL save file, so that you can do any further analysis you are interested in outside of SPEX but within the IDL environment. An easy way to accomplish this is to call

SPEX> idl, spex_save_data , /all, 'spex07202_bpowfit_spectra.sav’

Whenever we type idl, we use the command out of SPEX. The detail of SPEX commands can be found at Analyzing Flare X-ray Spectra Using SPEX.

Using the spex_save_data command, a file named ‘spex0723_bpowfit_spectra.sav’ will be saved in our IDL working directory. It contains all the spectral data and fitting parameters we need for later use. We can get the detailed explanation of all saved variables from help_spex.txt. Here I explain some of them for later use:

APAR_ARR FLOAT = Array(6, 37)

	Model best-fit parameters for each interval.

APAR_SIGMA      FLOAT     = Array(6, 37)

   Curvefit uncertainties on the best fit parameters in APAR_ARR

CHI             FLOAT     = Array(1, 37)

   Reduced chi square for each accumulation interval

XSELECT         DOUBLE    = Array(2, 37)

   Time boundaries on accumulation intervals referenced to the same base UT.

UTBASE          DOUBLE    =   7.4338560e+008

   The UT base

EDGES           FLOAT     = Array(2, 82)

   Detector energy channel edges

OBSI FLOAT = Array(82, 37)

   Count flux in accumulation interval I.

EOBSI           FLOAT     = Array(82, 37)

   Uncertainties on OBSI.

BACKI           FLOAT     = Array(82, 37)

   Background for accumulation interval.

CONVI          FLOAT     = Array(82, 37)

   Factor for conversion from count spectrum to photon spectrum.

When we want to restore the data in an IDL session, type

                    IDL> restore,  ‘spex0723_bpowfit_spctra.sav’

Use the help command to see which variables are defined

        IDL> help

5. How to plot the fitted spectra outside SPEX

                After restoring the data file, in order to plot the spectra, we first need to compute the photon flux, error bars, background flux, model flux and

                so on for a specific interval.

IDL>interval = 0

To compute the photon flux (=(obsi-backi)/convi) :

IDL>phflux = (obsi-backi)/convi

To compute the error of the photon flux (=eobsi/convi):

IDL>errflux = sqrt(eobsi^2 + ebacki^2) /convi

To compute the background flux (=backi/convi):

IDL>backflux = backi/convi

IDL>help, phflux, errflux, backflux

PHFLUX          FLOAT     = Array[82, 37]

ERRFLUX         FLOAT     = Array[82, 37]

BACKFLUX        FLOAT     = Array[82, 37]

Since the “edges” is the edges of energy channel, we need to compute the mean value of each channel:

IDL>print, edges[*, 0:2]

3.00000 4.00000

4.00000 5.00000

5.0 000 6.00000

IDL>edge_products, edges, mean = eph

IDL>print, eph[0:2]

3.50000 4.50000 5.50000

So far, we have all the spectral data for plots. We also need the f_vth_bpow model flux so that we can overplot on the data like SPEX does. From the documentation header of f_vth_bpow (examined with xdoc introduced in First Steps), we know f_vth_pow model uses f_vth.pro and f_bpow.pro to compute the thermal and nonthermal photon flux. The total photon flux is the sum of the two fluxes. Let us compute the model flux for the first interval using the saved fitting parameters from apar_arr.

First, we get the fitting parameters for the first spectral fit:

IDL> fitpara = apar_arr[*, interval]

IDL> print, fitpara

0.756227 3.20582 9.46661 3.31592 75.3729 3.83593

Secondly, we can check the UT time of that interval to see whether they are correct by comparing the time history window:

IDL>uttime = anytim(xselect[*, interval] + utbase, /ecs)

IDL>print, uttime

2002/07/23 00:27:24.000 2002/07/23 00:27:43.999

Thirdly, compute the thermal photon flux:

IDL>thermflux = f_vth(edges, fitpara[0:1])

IDL> help, thermflux

THERMFLUX FLOAT = Array[82]

Finally, we compute the nonthermal photon flux. But we should do something beforehand. There is always a low energy cut-off for the power-law photon distribution. From Spectroscopy - First Steps, we know the default is 1 keV cut-off, and spectral index below 1 keV is -1.5. The nonthermal f_ bpow model uses these two parameters, but in the saved fitting parameters, they are not saved. So we have to add them at the end of the variable array “fitpara” array:

IDL> a_cutoff = [1, 1.5]

IDL> fitpara = [fitpara, a_cutoff]

IDL> print, fitpara

0.756227 3.20582 9.46661 3.31592 75.3729 3.83593 1.00000 1.50000

IDL> nonthermflux = f_bpow(eph, fitpara[2:7])

IDL> help, nonthermflux

NONTHERMFLUX FLOAT = Array[82]

Now we have all the data we want to plot the spectrum of the first interval. The spectrum can be plotted as follows:

IDL> loadct, 4

IDL> window, 10, xsize = 600, ysize = 500

IDL> plot, eph, thermflux, /xlog, /ylog, xtitle = 'Photon Energy (keV)', ytitle = 'Flux (photons s!u-1!n cm!u-2!n keV!u-1!n)', xrange = [10.0, 300.], yrange=[1.0e-3, 1.0e6], xstyle = 1, ystyle = 1, /nodata

IDL> oplot, eph, thermflux, lineStyle = 1, color = 47, thick = 4 ;plot thermal flux

IDL> oplot, eph, nonthermFlux, lineStyle = 2, color = 110, thick = 4 ;plot nonthermal flux

IDL> oplot, eph, thermflux + nonthermflux, linestyle = 0, color =255 , thick = 4 ;plot total flux

IDL> oplot, eph, phflux[*, interval], psym = 1, color = 144, symsize = 1, thick = 2 ;plot spectrum data

IDL> oplot_err, eph, phflux[*, interval], yerr = errFlux[*,interval], psym = 1, symsize = 1, color = 144 ;plot spectrum data error bar

IDL> oplot, eph, backflux[*, interval], psym = 5, symsize = 0.5, color = 208 ;plot the background flux

Here is an IDL procedure that does all of the above steps. You can use it directly by cutting and pasting it into an IDL/SSW session.

The following is the spectrum we just plotted:

6. How to make a movie of the spectra

We already know how to plot the spectra in IDL (see Topic #5). We can use the IDL routine TVRD to read the current widow containing the spectrum for each interval. Then we use WRITE_JPEG to save it as a jpeg file. People who want to know the details of this should refer to Fanning’s book IDL Programming Techniques.

The following is a sample code:

IDL>image24 = tvrd(true = 1) ;read current IDL window

IDL>write_jpeg, jpegName, image24, quality = 75, true = 1 ;write to a jpeg file, jpegname is the name of the file

Each spectrum is thus saved as a jpeg file. There are many ways to make a movie from multiple jpeg files. Here are two of them:

(1) Using movie maker VideoMach (shareware, noncommercial use is free ! Windows user only.). It is very easy to make MPEG or AVI movies. Here is the sample movie. As you can see, we add the flare light curve at the top of the spectrum window and mark the time of each spectrum in lightcurve.

(2) using JSMOVIE in SSW. It is good for web publishing. Here is the sample movie.

7. How to choose and evaluate a background model

There are lots of useful hints and explanations about background subtraction in RHESSI Spectroscopy - First Steps. Since a reasonable background subtraction is a fundamental pre-requisite for RHESSI hard X-ray spectroscopy, here is something more about it.

When you analyze the flare of interest, be sure to first check it’s the long-term variation of the background from the quicklook observing summary plots, not only during the time of the flare but also for at least one RHESSI orbit before and after the flare. Here we will concentrate on the event that occurred on June 1, 2002 that peaked at 03:56 UT. The following window shows the observing summary from 01:00 to 07:00 UT, revealing a sinusoidal variation of the background emission due to changes in the magnetic latitude of the RHESSI spacecraft along its orbit (for details on the RHESSI background components see David Smith's Problems, Pecularities and Phenomena in RHESSI Spectroscopy document). The flare occurs around the maximum of the background emission when the spacecraft was at high geomagnetic latitudes. This makes the background subtraction a rather subtle issue.

bg01

For the analysis, set the observation time interval in the GUI to 1-Jun-03:20:00 - 1-Jun-04:15:00:

bg02

Then go to "Retrieve/Process Data" -> "Spectrum...". Choose all front detectors except 2 and 7 (the peculiarities of individual detectors are described in the Problems, Pecularities and Phenomena in RHESSI Spectroscopy document by David Smith). Set the time resolution to 4.0 s and choose bin code 6 (197 bins, 1-keV bands covering 3 to 100 keV, then 5-keV bands to 600 keV - check the button "Show Binning Codes" to see what choices you have). This is the best choice for this flare since there is emission detectable above 100 keV.

Click the “Plot Spectrum” button to accumulate the spectra and display them for all the intervals. Once the spectra have been obtained, click the button "Write output file..." to save the spectral data. Select the full response matrix ("0-Full calculation of diagonal, off-diagonal terms") -> "Write FITS file". Your window should look like:

bg03

Then call SPEX from an IDL session with the command "spex_proc". Read the spectral data and the spectrometer response (srm) matrix into SPEX and plot the count rates as follows:

SPEX> data,hessi,front
SPEX> _1file,hsi_spectrum_20020601_032000.fits
SPEX> dfile,hsi_srm_20020601_032000.fits
SPEX> graph
SPEX> energy_bands,6,12,12,25,25,50,50,300
SPEX> scale_bands,1,2,6,8
SPEX> th_ytype, 1
SPEX> graph

bg06

The sinusoidal background behavior tends to show up more prominently at high energies. So, you might try to use a higher-order polynomial fit at these energies. Let us apply a linear fit for the global energy range as well as for the energy band 6-12 keV, and a 3rd order polynomial fit for the higher energy bands (12-25, 25-50, 50-300 keV).

SPEX> back_order, 1,1, 3,3,3

Now, starting with the first band we have to choose reasonable intervals before and after the flare during which the background is calculated.

SPEX> use_band= -1
SPEX> background ; select graphically the interval to be used for background determination

The same for the other energy bands. And you probably want to change the y-range.

SPEX> th_yrange,1e-3,10
SPEX> use_band=0
SPEX> back ; select background intervals
SPEX> use_band=1
SPEX> back ; select background intervals
SPEX> use_band=2
SPEX> back ; select background intervals
SPEX> use_band=3
SPEX> back ; select background intervals

The two images illustrate possible background selections for use_band=0 (3-12 keV band, chosen background model: order 1) and use_band=2 (25-50 keV, chosen background model: order 3) .

bg06

bg07

There is no easy way within SPEX to plot the original count rates together with the modeled background. However, as always in SPEX you can execute every IDL-command from the SPEX prompt by typing "idl, command". Let us use this method to get a basic idea whether the derived background model is reasonable. The count rates before background subtraction are stored in the two-dimensional variable "flux" and the background data in "back":

SPEX> help, flux, back

The first index refers to energy bins, the second to time bins. See the help_spex.txt document for a description of parameters used within SPEX.

Before plotting the selected background, we need to know the array indices for the low end and high end energies of each energy band:

SPEX> print, value_locate ( edges[0, *], energy_bands )

3 8 8 22 22 46 46 ??

Let us first check the energy band 6-12 keV where the chosen model is a 1st order polynomial:

SPEX> idl, plot, total ( flux [ 3:8, * ], 1 ), /ylog, /xst
SPEX> idl, oplot, total ( back [ 3:8, * ], 1 )

bg08

Probably you'd prefer to plot the count rates as function of time, so let us do this for the 25-50 keV band:

SPEX> idl, utplot, ut [ 0, * ], total( flux [ 22:46, * ] , 1 ), /ylog, /xst, col=170
SPEX> idl, outplot, ut [0, *], total( back[ 22:46, * ] ,1)

bg09

In this case, the 3rd order plynomial fit yields a very good result. After evaluating the background model for the other energy bands, you may apply a 3rd order fit also to the first energy interval. To do this, set "back_order, 1, 3, 3, 3, 3" , and redo the background subtraction for use_band = 0.

If your background estimates turn out to be rather bad, then you may wish to start from scratch. To undefine any background you have already subtracted, use the command "background, clear" (sets the background to zero everywhere). Afterwards, don't forget to start with "use_band= -1".

Sometimes you may find that you don't have a long enough continuous observation interval before and/or after the flare (due to RHESSI night times or the spacecraft passing through the South Atlantic Anomaly) to choose a background model of an order as high as 3 (by the way: 3 is the highest order SPEX allows). In these cases, it's probably the most safe solution to choose for each energy band individually a short interval before and after the flare (1 min or so, depending on the actual context) and apply a linear fit, i.e. use "backorder, 1, 1, 1, 1, 1". In particular, in short events 1st order models give reasonable results also when the flare is located on top of a high-latitude variation. When evaluating your chosen background model, keep in mind that inaccuracies in the background subtraction have a relatively stronger effect at high energies where the flare count rates are low. For alternative approaches to the background determination (e.g., using the data of one RHESSI orbit before/after the flare) see the suggestion in David Smith's Behind, Beneath and Before HESSI Spectroscopy or "BBB". You can also use data from 15 orbits before and/or after the flare to try and duplicate the same range of geomagnetic parameters as for the flare observation itself.

Finally, select a time interval around the maximum of the 25-50 keV emission and plot the count spectrum.

SPEX> select ; select graphically a time range
SPEX> count

bg10

Comparing the background subtracted count spectrum (yellow) with the background spectrum (orange) you see that at the chosen time interval your flare spectrum is reliable up to about 200 keV. Above, the measured count spectrum becomes dominated by the background spectrum.

8. How to create and customize postscript plots within SPEX

SPEX allows you to create a postscript plot of the last plot(s) that you have done (i.e. in the time history or the spectral window) . The command is create_ps [, filename]. So, for instance, if your last plot has been the flare count spectrum, then create a ps-file with

SPEX> create_ps, countspec_01June02.ps

If you have selected multiple time ranges to make spectral fits at all these time steps (as discussed in Topic 2.), then the create_ps command makes a plot that contains all these spectral fits (and residuals), each on one page. This feature may be interesting to you when you want to have a quicklook to see if the fits you obtained in all the intervals are reasonable (check the fits as well as the residuals). You can enlarge the plot range to the full page size with calling:

SPEX> psstyle="FULLPORTRAIT"

before the create_ps command.

The standard setting of the psstyle parameter is "PORTRAIT", another option is "LANDSCAPE".

If you want to switch to color mode in the ps-plot, call first

SPEX> idl, sps, /color

Or at the IDL prompt

IDL> sps,/color ;enables color for Postscript Device and sets the graphic device to PS

IDL> X ; sets the graphic device to windowing device: x, win, and whatever for Mac

The create_ps command provides you with an easy tool within SPEX to make hardcopy outputs of your results. However, if you want to make plots with specific settings then you have to export your data to IDL (see Topic 3.) and use your own routines.

9. How to use RHESSI data of a different orbit for Background subtraction

In the following, we will describe how you can customize RHESSI data of a different orbit for background subtraction for analysis within SPEX. One possibility is to use data of one orbit before/after the flare (the duration of a RHESSI orbit is about 96 min 40 sec). If you are using the data of one day before/after flare, you have to adjust your time interval to ensure that you duplicate the same range of geomagnetic parameters as for the flare observation itself (check the observing summary). Be sure to select an orbit that does not contain a flare or an electron precipitation event.

In the GUI, save the spectrum and response matrix data of the interval covering the flare as well as the time interval used as background estimate. Be sure that both time intervals are exactly of the same length and that you use identical settings (detectors, energy bins, etc.)

The best thing is you save the observed count rates (+ the uncertainties) of the interval used as background to a file. One way to accomplish this is to call SPEX, read the spectrum and srm data of the background interval (lets assume we use one orbit before the flare) and then exit SPEX:

SPEX> data, hessi, front
SPEX> _1file, hsi_spectrum_1orbbef.fits
SPEX> dfile, hsi_srm_1orbbef.fits
SPEX> graph
SPEX> exit

Now, within IDL load the variables "flux" (observed count rate) and "eflux" (Gaussian uncertainties on "flux") from SPEX into your IDL session, rename the variables, and save them to an idl save-file:

IDL> flux_1orbbef = spex_current('flux')
IDL> eflux_1orbbef = spex_current('eflux')
IDL> save, flux_1orbbef, eflux_1orbbef, file='flux_1orbbef.sav'

With this procedure you have saved your data which you intend to apply as background count rates for later use.

In the following will be describe how to force SPEX to use these user-selected background data instead of one of the polynomial models which are inherent to SPEX. Call SPEX and read the spectrum and srm data of the flare interval. Do any background subtraction (SPEX needs that step, but later we will overwrite this background with the data from one orbit before), then exit SPEX:

SPEX> data, hessi, front
SPEX> _1file, hsi_spectrum_flare.fits
SPEX> dfile, hsi_srm_flare.fits
SPEX> graph

SPEX> background, null          ;skip background subtraction
SPEX> exit

Within IDL, we will now read the background data from the IDL save-file (which yields us the variables "flux_1orbbef " and "eflux_1orbbef ") and load the flare count rate data together with the uncertainties (variables "flux" and "eflux") from SPEX.
The new data for the background subtracted count rate (variable "rate") has to be calculated and together with the new background count rates (variable "back") to be loaded into SPEX. Note that we have to provide SPEX also with the actual data for the uncertainties ("erate" and "eback").

IDL> restore, 'flux_1orbbef.sav'
IDL> flux = spex_current( 'flux' )      ; flare count rates
IDL> eflux = spex_current( 'eflux' )   ; uncertainties of -"-
IDL> back_new = flux_1orbbef         ; new background count rates (= count rates of one orbit before)
IDL> eback_new = eflux_1orbbef     ; uncertainties of -"-
IDL> rate_new = flux - back_new      ; new background subtracted count rates

IDL> erate_new = sqrt( eflux^2. + eback_new^2. ) ; uncertainties of -"-
IDL> test1 = spex_current( 'rate', input=rate_new )
IDL> test2 = spex_current( 'back', input=back_new )
IDL> test3 = spex_current( 'eback', input=eback_new )
IDL> test4 = spex_current( 'erate', input=erate_new )

Here is an IDL procedure that does all of the above steps in IDL session. You can use it directly by cutting and pasting it into an IDL/SSW session.

If the loading of the new background and rate variables into SPEX was successful, then the variables "test1" to "test4" are set to the value 1. Now you can call SPEX again, and all the variables should have the intended values (i.e. the user defined background data are applied).

Note: Saving the actual SPEX session, e.g. with the command "save_event", does not save the new background and background subtracted count rates! However, if you save your results obtained within SPEX to an IDL save file, then you have all the data as you defined them.