RHESSI Spectroscopy -- Second Steps
Linhui Sui, Astrid Veronig, Brian Dennis, Richard Schwartz
List of topics*:
<![if !supportLists]>1. <![endif]>How to analyze a flare with different attenuator states.
<![if !supportLists]>2. <![endif]>How to select multiple analysis intervals.
<![if !supportLists]>3. <![endif]>How to change fitting models.
<![if !supportLists]>4. <![endif]>How to save and restore the fitting parameters and spectral data.
<![if !supportLists]>5. <![endif]>How to plot the fitted spectra outside SPEX.
<![if !supportLists]>6. <![endif]>How to make a movie of the spectra.
<![if !supportLists]>7. <![endif]>How to choose and evaluate a background model.
<![if !supportLists]>8. <![endif]>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
1. How to analyze a flare with different attenuator states
most strong flares, the attenuator state will be changed. The following plot is
an uncorrected quicklook light curve of the X4.8 flare on
As indicated by the purple horizontal bars, before UT, the attenuator state was A1 (thin in). From UT to UT, the attenuator state was mostly A3 (both thick and thin in) with four short periods in the A1 state. After 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.
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
This should give the following response if everything is OK:
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 UT and 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:
It will give a similar output as above. Again, for the analysis of spectra after UT, we need to switch back to the A1 srm file:
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:
is much more convenient to select time intervals from the window above to
analyze the main impulsive peak of the flare (from UT to 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:
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
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
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:
can set eplot to 1 so that we can see
the error bar in spectral data:
SPEX> eplot, 1
If we want to fit the spectrum at the second time interval, we need to reset ifirst and ilast:
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:
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
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:
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
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:
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
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
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:
DOUBLE = Array[2, 4]
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
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:19.999 UT. Then we change the fitting model to "f_vth_thin":
SPEX> f_model, f_vth_thin
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
The following plot is produced for interval 8 in the spectrum window:
We can check the fitting parameters as follows
4.25845 3.19432 17.1155 1.58129 131.000 2.47378 10.0000 5000.00
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>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
The following plot is the spectral fit with 'f_vth_thick' we got:
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
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
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
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
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 = eobsi
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]
<![if !supportLists]> 5.0 000 <![endif]>6.00000
IDL>edge_products, edges, mean = eph
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[*, ]
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[*, ] + utbase, /ecs)
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
Finally, we comput 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
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[*, ], psym = 1, color = 144, symsize = 1, thick = 2 ;plot spectrum data
IDL> oplot_err, eph, [*, ], yerr = errFlux[*,], psym = 1, symsize = 1, color = 144 ;plot spectrum data error bar
IDL> oplot, eph, backflux[*, ], psym = 5, symsize = 0.5, color = 208 ;plot the background flux
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
For the analysis, set the observation time interval in
the GUI to
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:
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> th_ytype, 1
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> back ; select background intervals
SPEX> back ; select background intervals
SPEX> back ; select background intervals
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) .
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 )
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)
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
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:
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
at the IDL prompt
IDL> sps,/color;enables color for Postscript Device and sets the graphic device to PS
; 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.
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
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> background, null
;skip background subtraction
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
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 )
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.