To run the HESSI analysis software, you need an IDL installation of version 5.3 or higher. When starting SolarSoft, you should include hessi, and the spex and xray packages. See the UCB Data Analysis Software and SolarSoft documentation web sites for more information.
In an attempt at source control, HESSI software is available both at its ``normal'' location in the SSW tree ($SSW_HESSI/idl), and in a release directory ($SSW_HESSI/release/idl). The most recent version of the software is in $SSW_HESSI/idl, while the release directory contains software frozen at a particular time. You can switch between the release and development version of the HESSI software using hessi_version - this should only be called at the beginning of an IDL session since it does not recompile routines - by default you should be using the development version (???) .
hessi_version, /dev
hessi_version, /release
The HESSI analysis code is object oriented. A list of the objects used is given on the UCB Web pages (http://hessi.ssl.berkeley.edu/software/reference.html). Since object oriented code may be new to some users, a Tutorial on objects has been prepared by André Csillaghy - see the Web page on URL http://hessi.ssl.berkeley.edu/software/hessi_oo_concept.html.
The HESSI software has memory - this allows the user to set up and then reuse parameters, but does require that some parameters be reset when a different option is required.
It is possible to run a lot of the HESSI software through the HESSI Graphical User Interface (GUI). This allows access to much of the software through pull-down menues and buttons. A brief summary of this is given in Section 2.8.
The observing summary data are obtained directly from the telemetry packets and the files allow the user to look at HESSI coutrates and other information before accessing the data files.
The photons are time-binned and the count rate is compressed to bytes (0.03% accuracy up to 1.0e6 counts/sec). Modulation variance is calculated for the 2 coarsest collimators by calculating the standard deviation of counts for the given 4 sec interval and dividing by the square root of the total counts in that interval. The spacecraft position (ephemeris) and data from the particle detector are also included in the files.
Information can be extracted using the plot and list methods:
obs_summ = obj_new('hsi_obs_summary') ; creating an instance of the class
time_interval = ['01-Sep-00 12:00', '01-Sep-00 14:30']
counts = obs_summ -> getdata(obs_time_interval=time_interval)
or
counts = obs_summ -> getdata(filename='hsi_20000830_235100_001.fits')
obs_summ -> plot ;default to countrates
text = obs_summ -> list()
Options to the plot methods include showing flare, eclipse, and SAA
intervals, plotting energy bands separately or summed, selecting
energy band range to plot, and plotting with colors. Also, any
standard plot keywords can be passed to the plot method. If you want
color plots, you must set up the colors first, e.g. by calling
linecolors. (?????? problem under UNIX)
obs_summ -> plot ;default to countrates
obs_summ -> plot, /countrate
obs_summ -> plot, /particle_rate
obs_summ -> plot, /ephemeris
obs_summ -> plot, /mod_variance
counts = obs_summ -> getdata(class_name='obs_summ_rate')
Options to the list methods include showing the list on the screen,
saving the list in a file, and printing the list.
text = obs_summ -> list() ;default to countrates
text = obs_summ -> list(/particle_rate)
text = obs_summ -> list(/ephemeris)
text = obs_summ -> list (/show, /print, /file_text)
text = obs_summ -> list (/ephem, file='ephem.txt')
It is also possible to plot the countrate against time using
utplot. To do this, you need to extract the time seperately:
counts = obs_summ->getdata()
times = anytim(obs_summ->getdata(/time),/ints)
utplot,times,counts(2,*)
The HESSI flare list can be retrieved using the function hsi_read_flarelist, which returns a structure of all the flares recorded to date.
flist = hsi_read_flarelist()
The flare list structure provides all sorts of information that is useful in selecting HESSI data: the start, peak and end time of the flare, the peak countrate, the flare location and the name of the FITS file that contains the flare data. We can examine the contents of the structure, by:
IDL> help,/st,flist(0) ** Structure HSI_FLARELISTDATA_EXT, 15 tags, length=144: ID_NUMBER LONG 1365 START_TIME DOUBLE 6.8368384e+08 END_TIME DOUBLE 6.8368386e+08 PEAK_TIME DOUBLE 6.8368384e+08 BCK_TIME DOUBLE Array[2] ENERGY_HI FLOAT Array[2] PEAK_COUNTRATE FLOAT 880.000 BCK_COUNTRATE FLOAT 59.6774 TOTAL_COUNTS FLOAT 8960.00 POSITION FLOAT Array[2] FILENAME STRING 'hsi_20000830_235100_001.fits' FLAGS BYTE Array[32] ACTIVE_REGION INT 0 X_POSITION FLOAT -892.500 Y_POSITION FLOAT 507.500 IDL>
If you wish to select a subset of the flares in the HESSI flare catalog, this can be done with either hsi_select_flare or hsi_flarecat. The hsi_select_flare function allows the selection to be made using keywords that correspond to the tags in the flare list structure, while hsi_flarecat provides a GUI interface to make the selection. In both cases, the results of the selection can be returned as either a list of flare IDs, a formatted list of the information relating to the flare, or a flare list structure.
; return formatted list of flares - search by range in keyword
flist = hsi_select_flare(id=[1390,1396])
flist = hsi_select_flare(start=time_interval, /form)
flist = hsi_select_flare(peak_countrate=[1500,2000], /structure)
hsi_flarecat, list=flist ; return list of flare numbers
hsi_flarecat, list=flist, /form ; return a formatted list
hsi_flarecat, list=flist, /structure ; return a structure
The flare times and position given in the flare list struture can be used to read the instrument data. In the example shown, the file starts at 23:51 UT on 30-Aug-00; the start and end times of the flare are 31-Aug-00 00:10:40 and 31-Aug-00 00:10:56, which can be determined as shown below; the flare position (in arcseconds in heliocentric coordinates) is given in the position tag. These can be fed into the HESSI analysis software, e.g. hsi_image or hessi_image using the obs_time_interval and position keywords (see Section 2.3).
t0 = anytim(flist(0).start_time,/yohkoh)
t1 = anytim(flist(0).end_time,/yohkoh)
time_interval = [t0,t1]
position = flist(0).position
o_img = hsi_image() ; create an instance of the class
im = o_img -> getdata(obs_time_interval=time_interval, $
xyoffset=position)
When reconstructing a HESSI image, you need a time range and the coordinates of the flare site. In the sections above, these were derived from the HESSI Observing Summary and Flare List. They can also be derived from SSW products produced when analysing data from other sources, e.g. SXT or TRACE. A ``map object'' produced by the Mapping software of Dominic Zarro, or an SSW-complient index structure can both provide the necessary information (see Section 2.6.1).
Assuming you have already made a map object, or have a set of index structures, proceed as shown below. Note that the HESSI software needs a time interval when selecting observations - for map objects, the single time needs to be converted into an interval.
position = [map.xc,map.yc] ; using Map structure
time_window, map.time, t0,t1, sec=[0,300] ; need a time interval
time_interval = [t0,t1] ; (maps only have single time)
position = [index.xcen,index.ycen] ; using SSW-complient index
time_window, index, t0,t1 ; retrieve start and stop times
time_interval = [t0,t1]
Start by creating an instance of the lightcurve class of objects, then plot an integrated light curve (with a 1 sec. resolution) for a selected time interval:
o_ltc = HSI_Lightcurve() ;creating an instance of the class
time_interval = ['2000/09/01 12:00', '2000/09/01 14:00']
o_ltc -> Plot, obs_time_interval=time_interval, ltc_time_resolution=1
Once the data has been ingested, it is also possible to retrieve arrays of required parameters. The call:
taxis = o_ltc -> GetData( /xaxis )
returns the start times of the time bins. While:
ltc = o_ltc -> GetData( time_range=[0,60], $
ltc_energy_band=[3,100,500,1000,15000] )
will return a 2D array filled with the counts accumulated in each time bin for each energy band. In this example, the time range of the lightcurve corresponds to the first 60 seconds in the selected time range and thus a 60×4 array is returned.
Lightcurves for selected detectors (in the range 1-9) or segments (in the range 1-18) can be retrieved with commands like:
ltc = o_ltc -> GetData( this_det = [1,2] ) ;detectors 1 and 2
ltc = o_ltc -> GetData( this_seg = [1,10] ) ;segments 1 and 10
You can plot any extracted parameter against time (taxis) using utplot:
utplot, taxis, ltc
Coincident events (see Section 3.3.2) can cause problems with ????, and you my want to look at them:
s = o_ltc -> GetData( /coincidence )
returns a 4-dim array where the 4th dimension contains the anti-coincident lightcurve (index #0) and the coincident lightcurve (index #1)
In order to reconstruct an image, certain pieces of information need to be supplied. These include the start time and ``integration'' time of the image, the center (in heliocentric coordinates) that the image should be reconstructed about, the image algorithm... ??? fov, energy range...
Note that the HESSI software needs a time interval when selecting observations, but a shorter interval (defined in time_range) should be used for image reconstruction.
After creating an instance of the image class of objects, you can list the parameters of the object using the ``print'' method.
o_img = hsi_image() ; create an instance of the image class
; of objects
o_img -> print ; list object parameters
time_interval = ['2000/09/01 12:00', '2000/09/01 14:00']
position = [200,400] ; x and y center, in heliographic coordinates
; (arcsecs from sun center)
im = o_img -> getdata(obs_time=time_interval, ; overall interval
time_range=[0,10], ; start, stop seconds from
; start of overall interval
xyoffset=position ; center for reconstruction
image_algorithm='clean') ; reconstruction algorithm,
; back-projection by default
o_img -> plot ; plot the image
1. Back Projection (plus Clean). The default image reconstruction algorithm is back projection with clean disabled. Back projection simply operates by multiplying the calibrated event list by the collimator modulation patterns to construct a dirty map of the image. The map can then be improved by using Clean - this is selected by passing ``clean'' into the image_algorithm keyword.
o_img = hsi_image() ; create an instance of the class
im = o_img -> getdata(obs_time=time_interval, $ ; default no Clean
time_range=[0,10], xyoffset=position)
;or
im = o_img -> getdata(obs_time=time_interval, image_algorithm='clean')
o_img -> plot
To generate a default model image use:
im = o_img -> getdata()
o_img -> plot
2. MEM-Sato is an image reconstruction algorithm utilizing the maximum entropy method (MEM) described by Sato, Kosugi and Makishima (1999, PASJ, 51, 127-150). In general, MEM-Sato gives very sharp images, but is time consuming to run. This is especially true if all the collimators are selected for the imaging. The finest four collimators (0 - 3) require substantial computing time since they have the most rapid modulations and hence require the finest time bins. Eliminating these collimators from the calculation will cause a loss of resolution, but will dramatically increase the speed of the image reconstruction. This is done through the det_index_mask keyword. For example, if you wish to use MEM-Sato to construct an image using only collimators 5 - 9 (indices 4 - 8), you would type:
o_img = hsi_image()
im = o_img -> getdata(obs_time=time_interval,image_algorithm='mem sato', $
det_index_mask=[0,0,0,0,1,1,1,1,1])
o_img -> plot
If you wish to use MEM-Sato to construct an image of the default model
using the same collimators, enter:
o_img = hsi_image()
im = o_img -> getdata(/plot,image_algorithm='mem sato', $
det_index_mask=[0,0,0,0,1,1,1,1,1])
3. MEM-Vis is an image reconstruction algorithm borrowed from radio astronomy techniques - see Conway et al., 2001, in preparation ???.
The MEM-Vis algorithm uses the same basic steps as MEM-Sato, but works in polar coordinates and uses ``visibilities''1 rather than ``counts''. In the HESSI context, counts and visibilities can simply be thought of as two time series representing the same information: the response of HESSI to an image. Whereas counts are non-negative integer values, visabilities are complex numbers (although the real and imaginary parts do not contain independent information). Since aspect correction can easily be applied when constructon visibilities from observed counts, visibility modulation patterns do not need to be corrected for aspect for each time bin and visibility values over several rotations can be simply added together. In constructing visibilites, there is effectively a smoothing. This relaxes problems arising from the presence of zero counts, but introduces a small systematic error and a correlation of noise on visibilities in the neighbouring time bins.
So, whereas MEM-Sato can be slow for a long flare, MEM-Vis does not take much longer than for a short flare since visibilities can be simply added together over multiple rotations.
o_img = hsi_image()
im = o_img -> getdata(obs_time=time_interval, image_algorithm='vis')
o_img -> plot
4. Pixon-based image reconstruction tries to maximize the image resolution in locations where the source intensity is high, but permits the resolution to be low where the source intensity is low. Thus, the ``pixels'' in a pixon code can vary in shape and size accross the image. The method used here is the fractal pixon method of Puetter and Pina (1994, Experimental Astronomy, 3, 293: see also Pina and Puetter 1993, PASP, 105, 630 and Metcalf et al. 1996, ApJ, 466, 585). This is a very powerful method, giving good noise suppression and good recovered photometry, but is very time consuming. It is suggested that users choose a different technique to decide the key issues of time coverage, energy range, etc. for their images, and then use the pixon method to create final images for analysis and display.
o_img = hsi_image()
im = o_img -> getdata(obs_time=time_interval, image_algorithm='pixon')
o_img -> plot
5. Forward Fitting assumes that the source can be represented by a number of superimposed 2D gaussians, which can be either cicular, elliptical, or curved. The algorithm adjusts the modulation patterns calculated from a model map containing multiple Gaussian sources to the calibrated event list of the telemetry file. The primary advantages of the forward fitting algorithm is its speed and relatively robust estimation of the source parameters. The current implementation allows model-fitting by an arbitrary number of 2D Gaussian sources. The initial default value is one single source (n_gaussians=1). If the resulting chi-square is near unity, no refinement of the model is justified. However, if the resulting reduced chi-square exceeds the value of ~1.1, the user should increase the number of gaussian components (e.g n_gaussians=2), until the model is consistent with a chi-square near unity. Because the forward-fitting algorithm only converges well if the initial guess of the source positions is accurate, closely spaced sources can cause a problem.
??? object code
hessi_image,im,par,/plot,image_algorithm='forward'
hessi_image,im,par,/plot,image_algorithm='forward',n_gauss=2
A detailed description of the algorithm and its testing can be found in Aschwanden et al. 2001 (ApJ, in preparation). Additional information on the forward fitting method can be found at the Lockheed-Martin web site.
Given the potential problems with image reconstruction (see Section
3.4), below are some recommendations for testing the
reality of features in your image:
1) Re-image the source with:
2) Apply redundancy tests:
3) Trace back the feature to:
4) Simulate and analyse similar sources
These considerations also reveal the importance of being able to simulate HESSI data (see section 2.9). The hessi_image software offers a fairly flexible means to contruct simulated data which can then be analyzed to see how the imaging software recovers what was included in the model.
Under construction!!! Below, we summarize the keywords that are most important for image reconstruction.
A complete list of keywords is given...
Start by creating an instance of the spectrum class of objects, then define the time interval of interest, and energy range and binning, etc. before retrieving the data:
o_spc = hsi_spectrum() ; create an instance of the class
time_interval = '31-aug-00 '+['11:40:00','12:00:00']
o_spc -> set, time_range=[0,300.], obs_time_interval=time_interval, $
sp_energy_binning=1, sp_time_interval=2.0, energy_band=[0.,0.]
spec = o_spc -> getdata() ; retrieve count-rate spectrum
The keyword sp_energy_binning selects files containing predefined energy bin ranges - in the example, the sp_energy_binning=1 selects the file ``ct_edges.1'', from spec_resp directory in the HESSI branch of the $SSWDB tree. Because several parameters need to be defined, the set method was used - it would also have been possible to pass these as parameters in the getdata method.
As with imaging, the this_det and this_seg keywords
can be
used to return only the spectrum for the desired detectors or segments.
spec = o_spc -> getdata( this_det=[1,2] )
spec = o_spc -> getdata( this_seg=[1,10] )
The average spectrum can be plotted on a log-log scale using:
plot, o_spc -> getdata(/xaxis,/mean), total(spec,2), psym=10,/xlog,/ylog
To retrieve and plot a calibrated spectrum, use:
spec = o_spc -> getdata(/sp_semi_calibrated)
plot, o_spc -> getdata(/xaxis,/mean), total(spec,2), psym=10,/xlog,/ylog
The boundaries of the energy and time bins used in producing the
spectra can be retrieved using:
energy_edge = o_spc -> getdata(/xaxis, /edges_2)
ut_edge = o_spc -> getdata(/yaxis, /edges_2)
The procedure hsi_spectroscopy_list lists information on HESSI spectroscopy.
Under construction!!! Below, we summarize the keywords that are most important for spectrocopy.
A complete list of keywords is given...
Under construction!!! See Richard Schwartz.
Sets of images of different spectral ranges for the same time interval. Currently needs to be done by repeated call to hsi_image.
maps of spectral hardness; temperature maps?; density?
Under construction!!!
Light-curves of features in image (from image CUBES) Currently needs to be done by repeated call to hsi_image.
When reconstructing a HESSI image, you need a time range and the coordinates of the flare site. These can be derived from SSW products produced when analysing data from other sources, e.g. SXT or TRACE.
A ``map object'' (produced by the Mapping software of Dominic Zarro), or an SSW-complient index structure, can both provide the necessary information. (Note: The TRACE index structure is SSW-complient, but since the SXT index stores the image centre in terms of CCD pixel coordinates rather than heliocentric coordinates, it is necessary to convert the SXT index structure to an SSW-compliant structure using struct2ssw.)
Assuming you have already made a map object, or have a set of index
structures, proceed as shown below. Note that the HESSI software needs a
time interval when selecting observations, but a shorter interval (defined
in time_range) should be used for image reconstruction.
position = [map.xc,map.yc] ; using Map structure
time_interval = anytim(map.time)+[0,300] ; need a time interval
;or
position = [index.xcen,index.ycen] ; using SSW index structure
fmt_timer,index,t0,t1 ; retrieve start and stop times
time_interval = [t0,t1]
o_img = hsi_image()
im = o_img -> getdata(obs_time_interval=time_interval, $
time_range=[0,10], xyoffset=position)
hmap = o_img -> makemap() ; make a map object
plot_map, map ; plot original data
plot_map, hmap,/over ; overplot HESSI
See the SSW Mapping tutorial, and Trace Analysis Guide for more information.!!
http://www.lmsal.com/solarsoft/sswdoc/solarsoft/ssw_eit_mdi.html http://www.lmsal.com/solarsoft/sswdoc/guides/tag/
Mention also the ancillary databases (section 2.7...
Also, compare in time...