DISCLAIMER: This document is still under development, as is the analysis software itself. Incomplete decriptions of various things will be found below. However, all example commands should work. Please report any errors when trying these examples to C. M. Johns-Krull at the above email address.
Contents:
I. Imaging
A.
Introduction
B.
RHESSI Imaging Software
i.
Bug Reporting
C.
Getting Started
D.
RHESSI Image Reconstruction Algorithms
i.
Back Projection with or without CLEAN
ii.
MEM SATO
iii.
MEM VIS
iv.
Pixons
v.
Forward Fitting
E.
Intermediate Data Products
F.
RHESSI Imaging Control Keywords
G.
Working with Simulated RHESSI Data
i.
The RHESSI Observing Summary
ii.
Using the RHESSI Flare List
H.
Working with Early RHESSI Flight Data
i.
Using the Lightcurve Object to Find Flares
ii.
Making the First Image
iii.
Refining the Image
I.
Some Hints on Analyzing RHESSI Images
J.
Saving and Recovering Images and Writing and Reading FITS Files
K.
Simulating RHESSI Data
i.
RHESSI Simulation Keywords
ii.
Using HSI_SIM_FLARE
L.
Interaction Between the RHESSI GUI and the Command Line Interface
II. Spectroscopy
A.
Introduction
B.
RHESSI Spectroscopy Keywords
C.
Getting Started with SPEX
The RHESSI mission will perform both hard X-ray and gamma-ray imaging and spectroscopy of solar flares at an unprecedented level. The RHESSI Instrument Guide at http://moat.nascom.nasa.gov/~bentley/guides/hag provides a description of the RHESSI satellite and its general scientific capabilites. Here, we focus on the data analysis software for reconstructing images using RHESSI data. This document describes the command line interface to the RHESSI imaging software. It is envisioned that most investigators will use the RHESSI Imaging GUI for basic imaging tasks. A description of the GUI and how to get started with it can be found at http://hesperia.gsfc.nasa.gov/ssw/hessi/doc/gui_help.htm. The primary benefit of the command line interface is that it allows a level of flexibility that may not always be present with the GUI. In addition, the command line interface offers a fairly straightforward means to simulate RHESSI images of user defined sources. Many users may find it desireable to use the command line interface and the GUI to accomplish some tasks (such as creating simulated data on the command line and then analyzing the data in the GUI). This is supported to a limited degree and is described toward the end of the imaging section. Finally, the command line interface is very valuable for batch processing of RHESSI images. In addition to the documentation given here, additional software documentation is available on the web. In particular, several RHESSI data analysis software manuals can be found at:
http://hessi.ssl.berkeley.edu/software/
and
http://hesperia.gsfc.nasa.gov/rhessidatacenter/software/documentation.html
These pages includes information and links to this manual, a manual for using the data analysis GUI, descriptions for installing the RHESSI data analysis software and setting up the required system variables, and how the RHESSI test data (and soon the real data) can be accessed. RHESSI data analysis software uses many keywords to set control parameters. A useful reference is the full list of RHESSI data analysis keywords which can be found at:
http://hesperia.gsfc.nasa.gov/ssw/hessi/doc/hsi_params_all.htm
Finally, the RHESSI data center at
http://hesperia.gsfc.nasa.gov/rhessidatacenter/
is an excellent location to find the latest information on RHESSI data and software.
This manual deals only with RHESSI imaging software and assumes you have already installed and are using the Solar SoftWare (SSW) package. Instructions for downloading and installing SSW can be found at the above links or by going here. When installing SSW and the RHESSI software tree, you should pay particular attention to setting up access to the test data and defining the system variable $HSI_DATA_ARCHIVE. This variable is required for the proper use of the examples listed below for working with the simulated RHESSI telemetry data. This information is referenced from the software pages linked above or you can go directly to http://hessi.ssl.berkeley.edu/software/test_data_access.html which describes access to the RHESSI test data and setting up the proper system variables.
RHESSI imaging software is object oriented software. Basics of IDL objects and object oriented programming can be found in the IDL refernces manuals, particularly the guide on Objects and Object Graphics. An overview of the RHESSI analysis software objects and their use can be found at http://hessi.ssl.berkeley.edu/software/reference.html by following the link to Overview. While the software itself is written as object oriented code, users do not need to know anything about object oriented programing to use the RHESSI data analysis software. At the command line, users can use the program hessi_image.pro to perform all the imaging tasks supported by the data analysis software. This program is used like a standard IDL program, so those users familiar with IDL, but not with object oriented syntax or programming, can quickly learn to use hessi_image without investing any time in understanding objects if they so desire; however, it is encouraged that all users learn the object oriented syntax for using the data analysis software since this syntax affords more flexibility in using the software.
This document will describe the use of hessi_image.pro to get those wishing to dive in quickly the opportunity start using the software without having to learn the object syntax; however, examples of how to perform the same functions are given using the object calls for those users who wish to begin to learn that syntax as well. To learn more about the RHESSI software objects, go to http://hessi.ssl.berkeley.edu/software/reference.html . In an effort to encourage the use of the object syntax though, examples using hessi_image will only be used in the first part of the document, gradually giving way to the object oriented syntax. However, all the imaging control keywords used in the later object syntax examples can also be passed as keywords to hessi_image. It should be noted that the hessi_image command and the object oriented syntax are not interchangeable in the sense that the user can not use one syntax and then immediately switch to the other and expect all the proper internal variables to be defined. In any one session, the user should stick with one syntax or the other and not change back and forth. It is also the case that using hessi_image results means that a new reconstruction is performed each time the routine is called with some setting (e.g. the time range or energy band) changed. If the user wishes to make many changes to the various imaging settings, this can result in rather long command line calls. On the other hand, the object oriented calls allow the user to set parameters without actually performing the image reconstruction using the set utility. As a result, the user can make several different calls to set various control parameters, then make one call to reconstruct the image. Some users may find this capability useful, and as a result are encouraged to learn the object oriented syntax calls. Furthermore, in the discussions of the different imaging algorithms, links will be given to tutorial pages written by the authors of the various algorithms, and these pages primarily use the object syntax, so it is desireable to start learning the syntax in order to better use these pages.
If you encounter bugs when using the RHESSI data analysis software,
it is important to notify the software team so that these can be corrected.
In general, there are two potential types of bugs: ones which result from
true programming bugs in the data analysis software, and ones which result
from using various control keywords outside their intended ranges or which
conflict with one another. To the new user, these may be hard to
distinguish. Therefore, it is requested that users please read through
the list of RHESSI
Data Analysis Software Frequently Asked Questions as well as consulting
the list
of RHESSI keywords in order to look for default settings and to get
a description of what the keyword is used for before reporting bugs.
It is also useful for the user to consult the pages (linked below) for
the respective imaging algorithm they are having trouble with which may
provide insight on the use of the control keywords. Finally, it is
also advisable to check the Software
Bug Reports page where the user can browse or search previous reports
to look for instances of similar problems detected by other users.
Once these other sources are exhausted, if you still experience trouble,
the bug reports page also gives an email address for submission of a report.
This address is: hessibugs@hesperia.gsfc.nasa.gov.
Note: Please do not report bugs related to software installation
issues. Please go the SSW installation pages or RHESSI installation
pages linked above to get help and report bugs related to installation.
This document is intended to serve as the main introduction to getting started using the command line interface for the data analysis software. In addition to this document, there is a streamlined list of useful commands you can get by typing hessi_list at the IDL prompt:
IDL> hessi_list
Typing this will produce a pop up window which summarizes useful and common commands for interacting with RHESSI data.
The RHESSI software tree of SSW contains both a release area for software and a development area. In theory, the release area contains stable, tested software for RHESSI data analysis, while the development area contains routines that are under development. However, at this time, due to the rapid development of the RHESSI data analysis software after launch, the release area is very out of date. All users should use the development area. It is preferrable that users edit their setup files so that the development area is selected when SSWIDL is invoked. In general (and in the future particularly), most users will want to use the release area of the software; however, this area may not contain the very latest in software functionality. The default area you are using depends on how you set things up when you install SSW and the RHESSI software. However, you can easily select either the release or development area at the begining of an IDL session. To do so, simply type:
IDL> hessi_version, /dev
for the development area or
IDL> hessi_version, /release
for the release area. Either call should be made only at the begining of an IDL session. These calls do not recompile routines, so it is not possible to switch back and forth between the two areas in a given IDL session. Unlike hessi_version, the hessi_list command described above can be executed at any time.
RHESSI imaging software defaults to assuming you will be working with real RHESSI data. However, new users wanting to get a feel for how the software works, as well as experienced users who may wish to simulate data in order to test the validity of their real images, may find it desireable to run in simulation mode. If you do not run in simulation mode, you must specify a time interval to work with to get real data (see below), whereas is you are running in simulation mode, a default model is available for use. If you wish to run in simulation mode, you should execute the following command at the begining of your SSWIDL session:
IDL> hsi_switch,/sim
With the above command excuted, you are ready to simulate RHESSI data, and if no model is set specifically, a default model will be generated. Instructions are given in section H on working with real flight data, but the examples given up to that point are simulated data running in simulation mode. If you have used hsi_switch to work in simulation mode, you can get started with the imaging software using hessi_image.pro by simply typing at the IDL command line:
IDL> hessi_image,im
If you have not already set hsi_switch to the simulation mode, the above syntax will generate an error since no observing interval has been set. With this call, the software will simulate RHESSI data using a default model and then analyze the data using the RHESSI imaging software. The resulting image, using the back projection algorithm, is returned in the variable im. The goal is that, by default, the back projection algorithm will return an image whose peak is equal to the total number of photons falling within the specified energy range and time interval from a point source. This is precisely true for a single point source using the back projection imaging algorithm with weights (the default setting). In the case of two equally strong point sources such as that in the default RHESSI simulation, the peaks in the reurned image correspond to the total number of photons falling within the specified energy range and time interval from the associated source. Most users will want to use one of the additional RHESSI imaging algorithms to clean up the back projected image. When using these other algorithms, the returned image (im) will have units of photons cm-2 s-1 arcsecond-1 over the energy band specified when constructing the image. The default model consists of two point sources separated by about 17 arcseconds in a field located 600 x 200 arcseconds away from the center of the Sun. When hessi_image is run the following operations occur: 1.) An event list for the default model is created assuming a duration of 4 seconds; 2.) The event list is transformed into a calibrated, binned event list; 3.) The modulation patterns associated with the calibrated event list are created; and 4.) The count rate is back-projected using the modulation patterns, producing what is known as a dirty map, which is returned in the variable im. This image can be plotted with the traditional SSW routine plot_image.pro. The simple call would be:
IDL> plot_image,im
However, by selecting the keyword plot in the call to hessi_image, the program will automatically display the final image at the end of its processing, along with a summary of the information that went into creating the image. The default plot utility within the RHESSI data analysis software is a GUI environment (called plotman) that offers substantial flexibility in displaying and examining the final image, as well as creating output for printing or including in presentations. The call would be:
IDL> hessi_image,im,/plot
The resulting image should look similar to the following:
Figure 1: A back projection of the default RHESSI simulation.
Note that the actual peak image level and total number of counts in your image will be slightly different as a result of differences in the random noise created in the simulation of the data. This is true for all figures shown in this document.
In this case (using plotman) the plot displayed on the screen (shown in Figure 1 above) also shows printed summary information on the date and time of the data use, which collimators were used in the analysis, the energy range of the photons considered, and the imaging algorithm used. The RHESSI data analysis software is designed to only perform new tasks when some parameter is changed in the call to hessi_image and to not redo steps that do not need to be updated. For example, if you issue the hessi_image command without the /plot keyword, you can then produce a plot by simply retyping the same command with the /plot keyword. The software will not recalculate the image, but simply display the image in the plotman utility. The above command can also be given without the return variable, im, in which case the resulting image is only displayed to the screen.
Using the object syntax, the above tasks can be performed by entering:
IDL> obj_im = hsi_image()
IDL> im = obj_im -> getdata()
The first command defines obj_im as an image object. This call must be given at the begining of any session when analyzing imaging data. Multiple imaging objects can be defined as well if the user desires to interactively process different data sets or to process the same data set with different parameters but does not always want to reset the different parameters when switching bewteen the two. The second command gets the image reconstruction and loads it into the variable im, which is the same as that created by hessi_image. This image can again be displayed using plot_image.pro as before, or it can be displayed using the plot object:
IDL> obj_im -> plot
which calls the standard SSW plotting routine (i.e. not the plotman utility). To invoke plotman you could simply type
IDL> obj_im -> plotman
which uses the same RHESSI specific plotting routines which are called by hessi_image. The default simulation was used in this case because no arguments specifying a model or a telemetry data file were passed in the call to getdata. To get the image and plot it at the same time, the second call would be:
IDL> im=obj_im -> getdata(/plot)
Note that when using the plot keyword in the object oriented call, the plot is displayed in the standard SSW plotting environment and not the plotman utility. The plotman utility cannot be invoked by a keyword setting in the getdata call.
It is important to note that the RHESSI imaging software has memory. Once invoked, repeated calls to hessi_image will not do anything new unless something (parameter, filename, etc.) is explicitly changed.
The current imaging parameters can be returned in an optional second return variable. For example, the call
IDL> hessi_image,im,pars,/plot
will display the final image, return the image in the variable im, and will return the imaging control parameters in the structure pars. This structure can be examined by typing:
IDL> help,/structure,pars
An example of the resulting help display is:
** Structure <2a017c0>, 138 tags, length=1536,
refs=1:
FLARE_ID_NR
LONG
0
IMAGE_ALGORITHM STRING
'back projection'
IMAGING_METHOD STRING
''
ALGORITHM_USED STRING
'HSI_BPROJ'
ALGORITHM_UNITS STRING
'No Units'
ALG_UNIT_SCALE FLOAT
1.00000
PIXEL_AREA
FLOAT 16.0000
IMAGE_UNITS
STRING ''
ALGORITHM_AVAILABLE
STRING = Array[6]
IMAGING_STRATEGY
STRING = 'HSI_ANNSEC_PATTERN'
IMG_STRATEGY_AVAILABLE
STRING = Array[2]
FACTOR_BY
INT
1
MODPAT_SKIP
INT
1
R0_OFFSET
FLOAT 2560.00
IMAGE_DIM
INT Array[2]
PIXEL_SCALE
FLOAT 1.00000
PIXEL_SIZE
FLOAT Array[2]
RMAP_DIM
INT Array[2]
N_MODPAT
LONG
0
VERBOSE
BYTE 0
NEED_UPDATE
BYTE 0
LAST_UPDATE
DOUBLE 1.0129684e+009
DEBUG
BYTE 0
PLOT
BYTE 0
XYOFFSET
FLOAT Array[2]
USE_AUTO_TIME_BIN
LONG =
1
CBE_DIGITAL_QUALITY
FLOAT = 0.900000
CBE_POWERS_OF_TWO
LONG =
1
USER_FLUX_VAR INT
0
USE_FLUX_VAR LONG
0
SMOOTHING_TIME FLOAT
0.500000
SRT_FILENAME STRING
'SRT_2*.dat'
CBE_DET_EFF
STRUCT -> <Anonymous> Array[1]
FRONT_SEGMENT BYTE
1
REAR_SEGMENT BYTE
0
ENERGY_BAND
FLOAT Array[2]
TIME_BIN_DEF FLOAT
Array[9]
TIME_BIN_MIN LONG
2048
DET_INDEX_MASK BYTE
Array[9, 3]
POISSON
BYTE 0
SEED
FLOAT 0.000000
COINCIDENCE_FLAG
INT =
0
OTHER_A2D_INDEX_MASK
BYTE = Array[27]
SEG_INDEX_MASK BYTE
Array[18]
REBIN_METHOD STRING
''
REBIN_POISSON_FLAG
INT =
0
REBIN_SEED
INT
0
CONTIG_ENERGY_EDGES
INT =
1
SUM_FLAG
INT
1
N_BIN
LONG Array[9, 3]
BINNED_N_EVENT LONG
Array[9, 3]
BINNING
STRUCT -> <Anonymou
STRING = ''
DIM_SPECTROGRAM POINTER
<NullPointer>
COUNTS_2_UNIT POINTER
<NullPointer>
TOTAL_COUNT
POINTER Array[27]
SEC2TIME_UNIT DOUBLE
0.00000000
A2D_INDEX_MASK BYTE
Array[27]
NEW_GAIN
INT
1
GAIN_TIME_WANTED
DOUBLE = 0.00000000
GAIN_GENERATION INT
1
TIME_RANGE
DOUBLE Array[2]
TIME_UNIT
INT
1
N_EVENT
LONG
11208
UT_REF
DOUBLE 7.0912080e+008
CHANNEL_RANGE POINTER
<NullPointer>
COINCIDENCE_MASK
POINTER = <NullPointer>
ABSOLUTE_TIME_RANGE
DOUBLE = Array[2]
SIM_TIME_UNIT INT
1
SIM_UT_REF
DOUBLE 7.0912080e+008
SIM_TIME_RANGE DOUBLE
Array[2]
SIM_USE_SPECTRUM
BYTE = 0
SIM_PHOTONS_PER_COLL
FLOAT = 10000.0
SIM_TIME_PROFILE_PTR
POINTER = <NullPointer>
SIM_PHOTON_FLUX_PTR
POINTER = <NullPointer>
SIM_PHOTON_ENERGY_PTR
POINTER = <NullPointer>
SIM_PHOTON_TIMES_PTR
POINTER = <NullPointer>
SIM_SPEC_PARS POINTER
<NullPointer>
SIM_SPEC_MODEL STRING
''
SIM_INTERP_TPROFILE
BYTE = 0
SIM_INTERP_DT FLOAT
0.00781250
SIM_JUST_BACKGROUND
BYTE = 0
SIM_BACKGROUND FLOAT
0.000000
SIM_BKGD_TIME_UNIT
INT =
0
SIM_BKGD_UT_REF DOUBLE
0.00000000
SIM_BKGD_TIME_RANGE
DOUBLE = Array[2]
SIM_USE_BKGD_SPECTRUM
BYTE = 0
SIM_BKGD_TIME_PROFILE_PTR
POINTER = <NullPointer>
SIM_BKGD_PHOTONS_PER_COLL
FLOAT = 0.000000
SIM_BKGD_PHOTON_FLUX_PTR
POINTER = <NullPointer>
SIM_BKGD_PHOTON_ENERGY_PTR
POINTER = <NullPointer>
SIM_BKGD_PHOTON_TIMES_PTR
POINTER = <NullPointer>
SIM_BKGD_SPEC_MODEL
STRING = ''
SIM_BKGD_SPEC_PARS
POINTER = <NullPointer>
SIM_BKGD_INTERP_TPROFILE
BYTE = 0
SIM_BKGD_INTERP_DT
FLOAT = 0.000000
SIM_ASPECT_SOLUTION
OBJREF = <NullObject>
SIM_A2D_INDEX_MASK
BYTE = Array[27]
SIM_PIXEL_SIZE FLOAT
1.00000
SIM_XYOFFSET FLOAT
Array[2]
SIM_MODEL_PTR STRUCT
-> GAUSSIAN_SOURCE_STR Array[2]
SIM_ENERGY_BAND FLOAT
Array[2]
SIM_GAIN_GENERATION
BYTE = 1
SIM_GAIN_TIME_WANTED
DOUBLE = 0.00000000
SIM_MAX_SIZE FLOAT
128.000
SIM_SASZERO
BYTE 0
SIM_SRT_FILENAME
STRING = 'SRT_2*.dat'
SIM_SEED
INT
0
SIM_NBUFF
FLOAT 0.000000
SIM_IDPU_VERSION
BYTE = 0
SIM_CRYOCOOLER_POWER
BYTE = 1
SIM_COLD_PLATE_TEMP
BYTE = 20
SIM_IDPU_TEMP BYTE
20
SIM_HEATER_POWER
BYTE = 1
SIM_HV_POWER BYTE
1
SIM_ACTUATOR_POWER
BYTE = 1
SIM_FAST_RATES BYTE
0
SIM_SC_TRANSMIT BYTE
0
SIM_INSUN
BYTE 1
SIM_SSR_STATE BYTE
0
SIM_ATTEN_STATE BYTE
1
ASPECT_MODE
BYTE 0
ASPECT_CNTL_LEVEL
BYTE = 0
ASPECT_SIM
BYTE 0
ASPECT_TIME_RANGE
DOUBLE = Array[2]
RAS_TIME_EXTENSION
DOUBLE = Array[2]
SASZERO
BYTE 0
BPROJ_ALG_AVAILABLE
STRING = Array[2]
FLATFIELD
INT
1
VRATE
POINTER <NullPointer>
WEIGHT
INT
1
NATURAL_WEIGHTING
INT =
1
UNIFORM_WEIGHTING
INT =
0
SPATIAL_FREQUENCY_WEIGHT
FLOAT = Array[9]
TAPER
FLOAT 0.000000
Using the object syntax, the above parameters can be printed to the IDL window by calling the Print object:
IDL> obj_im -> print
If you wish to load these parameters into a structure (for example to save), the following object call can be used:
IDL> pars = obj_im -> get()
The objects allow you to examine only the parameters that control the image reconstruction using either:
IDL> obj_im -> print, /control_only
or
IDL> control_pars = obj_im -> get(/control_only)
IDL> help,/structure,control_pars
which will produce a listing similar to the one above but only for parameters directly related to the image reconstruction. There is also a GUI environment for looking at the parameters, handled by the params utlility. This can be invoked by simply typing:
IDL> obj_im -> params
This will pop up a text widget on the screen listing the object parameters. So far, only the basic parameters users may be interested in are shown to keep the list simple. This list will be added to as users indicate which of the parameters they are most interested in seeing. The entire listing can be viewed with the print utility as described above.
The keyword tags shown above reflect two distinct jobs being handled by the software: 1. image reconstruction and 2. data simulation. All of the tags which do not begin with sim_ either control some aspect of the image reconstruction or give information regarding the reconstructed image. All of the tags which do begin with sim_ control some aspect of the data simulation. Each of the imaging or sim_ structure tags can be passed in as a keyword to hessi_image or in the GetData object call in order to change how the image is reconstructed. For example, suppose you wanted to construct an image using a restricted energy band of 30 to 80 keV, you would then type:
IDL> hessi_image,im,pars,/plot,energy_band=[30.,80.]
The object syntax for this task would be:
IDL> im = obj_im -> getdata(energy_band=[30.,80.])
There is also a set method that can be used to set image reconstruction parameters before calling GetData. This allows you to set many different parameters in different calls without forcing the software to perform the image reconstruction each time you call set. This allows the user to avoid extremely long, multiple line calls to set image control parameters. Note that this is not possible with hessi_image. Each time you call this routine with some image control parameter changed, the software will reconstruct the image. As an example of using set, if you wish to set the energy range to 10 to 100 keV after making the image in the 30 to 80 keV band, you could type:
IDL> obj_im -> set, energy_band=[10.,100.]
If you then wish to change the imaging algorithm to clean (discussed below), you could then type:
IDL> obj_im -> set, image_algorithm='clean'
After each of the above calls, these parameters will be set, but the software will not yet have reconstructed the image. To do so and then display it, type:
IDL> im = obj_im -> getdata()
IDL> obj_im -> plotman
Using either calling method to set the energy band to 30 - 80 keV, a new image would then be created from the detected photons that fall between these two energy limits (the default energy band is 12 to 25 keV - see RHESSI Imaging Control Keywords below). By creating images in several different bands, you can begin to explore differences in the spectrum of different sources in a given image. Another example of using the keywords is to change the imaging algorithm. Suppose you want to use the mem sato image reconstruction method, you could simply type:
IDL> hessi_image,im,pars,/plot,image_algorithm='mem sato'
Using object syntax:
IDL> im = obj_im -> getdata(image_algorithm='mem sato',/plot)
Again, it is also important to note that the software remembers any changes made to the imaging control parameters in all subsequent usage. In the examples given above, the mem sato algortihm will continue to be used for each new image until it is explicitly changed, so if you choose to change the energy band or the time range of the image, the mem sato algorithm would continue to be used. To change back to the back projection algorithm for example, simply type:
IDL> hessi_image,im,image_algorithm='back projection'
Using object syntax:
IDL> im = obj_im -> getdata(image_algorithm='back projection')
or
IDL> obj_im -> set,image_algorithm='back projection'
IDL> obj_im -> getdata()
These examples illustrate an important point when culling through data to look for flares with particular properties of morphologies. The back projection algorithm is the simplest and fastest of the imaging algorithms; however, the other algorithms yield better images, but of course they take longer to run. As a result, you will probably find it adviseable to search the RHESSI data using back projection or possibly one of the other algorithms with various keywords set to optimize speed. Once you find an event you are particularly interested in, it can be explored in more detail with the more sophisticated algorithms.
RHESSI Image Reconstruction Algorithms
Currently, there are 5 different image reconstruction algorithms available for use with hessi_image. These are: 1.) back projection + clean, 2.) mem sato, 3.) mem vis, 4.) pixons, and 5.) forward fitting. Here, we give a brief description of each algorithm, its basic usage, some hints about appropriate usage, and a list of important keywords used with each algorithm. Additional keywords can be found in the list of Image Control Keywords given below.
1. Back Projection with or without Clean is the default image control algorithm 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. This map can then be improved by cleaning, which is set by passing 'clean' into the image_algorithm keyword. Thus, the IDL call:
IDL> hessi_image,im,/plot
Using object syntax:
IDL> im = obj_im -> getdata(/plot)
will construct a dirty map of the default model image assuming no data files or other image simulation parameters have been changed. If operating on the default RHESSI simulation, the resulting image should look like that shown in Figure 1 above. To clean the resulting image, you can simply type:
IDL> hessi_image,im,image_algorithm='clean',/plot
Using object syntax:
IDL> im = obj_im -> getdata(image_algorithm='clean',/plot)
The resulting cleaned image produced from the default simulation is shown in Figure 2 below. Again, the exact numbers with regard to peak value and total counts may be slightly different as a result of the random noise used in the simulation.
Figure 2: A clean reconstruction of the default RHESSI simulation.
Note that this second call to the analysis software does not depend on your having already executed the first call. If you had typed one of the commands from the second set first, the default model image would be reconstructed using the clean algorithm, but the user would not have had an opportunity to first look at or manipulate the dirty map.
Recall, as always when working with the objects, you must have already declared obj_im as an imaging object with the call
IDL> obj_im = hsi_image()
All the basic image control keywords related to image size, position, pixel size etc. can be used with back projection. These keywords can be listed by clicking on the all link under the Image: General heading in the list of RHESSI Data Object Parameters. The use of clean is unique to the back projection image algorithm. The clean algorithm used here is a straight forward iterative procedure that takes the dirty map, finds the brightest pixel, uses some fraction of this value and the modulation patterns to compute its side lobes and subtracts the result from the dirty map. The process is then repeated until the brightest pixel is negative, or the maximum number of iterations has been reached. The maximum number of iterations is one of a number of clean parameters that can be specified. The parameters related to clean are passed back as part of the imaging control parameters. Thus, in a call such as:
IDL> hessi_image,im,pars,/plot,image_algorithm='clean'
IDL> help,/structure,pars
Using object syntax:
IDL> im = obj_im -> getdata(image_algorithm='clean',/plot)
IDL> obj_im -> print
or
IDL> im = obj_im -> getdata(image_algorithm='clean',/plot)
IDL> pars = obj_im -> get()
IDL> help,/structure,pars
The first of the two object oriented calls simply prints the imaging parameters to the screen. The second basically does the same thing, but also loads them into a variable which the user may also find useful, for example to save along with an image or for manipulating particular returned variables.
it is possible to examine the control keywords used by clean. Each imaging algorithm has specific image control parameters for that algorithm. To help keep these separate, all parameters for a specific algorithm have a unique prefix. For CLEAN, the prefix is clean_ (for MEM SATO the prefix is sato_, for MEM VIS vis_, for PIXONS pixon_, and for FORWARD FITTING is ff_). The following tag names in the structure pars define the clean related image parameters. Some of these important parameters with their default values are:
clean_niter - (default = 33) max number of clean iterations
which will be performed
clean_negative_max_test - (default = 1) a flag for stopping
when the highest absolute value comes from a negative peak.
clean_frac - (default = 0.1) fraction to multiply the maximum
by when subtracting the point-spread function (PSF) from the dirty map.
clean_sigma_factor - (default = -1.0) factor by which the pitch
of the grid parameters will be multiplied by to generate the assumed width
(Gaussian sigma) of a point source.
clean_sigma_beam - (default = 9 X 3 array) the assumed width
(Gaussian sigma) produced in the image by a point source.
clean_taperpsf - (default = 9 X 3 array) tapering weights
which will be applied to each collimator.
clean_chi_sq_crit - (default = -1.0) the maximum chi square
value obtained by clean.
The use of the unique prefixes for each algorithm also make it possible for the user to select only the parameters related to a specific algorithm. For example, if you wish to print only the parameters related to CLEAN, you would type:
IDL> obj_im -> print,/clean
This will print the parameters to the screen, and you will see a longer list that describe above. The list above gives the control parameters the user may wish to modify, while the additional parameters printed with the above command are primarily used internally by the algorithm. If you want to load these CLEAN specific paramters into a data structure and then examine the resulting structure, you can simply type:
IDL> clean_pars = obj_im -> get(/clean)
IDL> help,/structure,clean_pars
Simple back projection without clean is a fast image reconstruction technique which is good for initial examination of an image. For example, if you are looking for specific instances (say flares with at least 3 roughly equally bright sources), back projection without clean is useful for looking at several images to select out times for additional analysis. Once you have honed in on a subset of data you are particularly interested in, you will almost certainly want to sharpen the image, at least by using the image_algorith='clean' keyword, and probably by using one of the more sophisticated image reconstruction algortihms. The final resolution of a cleaned image depends on the collimators used in constructing the image. For example, using the coarsest collimators provides little information on sources smaller than approximately the resolution of these collimators. As a general guide, the resolution of the collimators are 1 - 2."3, 2 - 3."9, 3 - 6."9, 4 - 11."9, 5 - 20."7, 6 - 35."9, 7 - 62."1, 8 - 107", and 9 - 186". This information as well as other relatively detailed information on the RHESSI mission is available from the RHESSI Instrument Guide available at http://umbra.nascom.nasa.gov/~bentley/guides/hag.
Additional documentation and a more detailed tutorial on using clean within the RHESSI data analysis suite can be found on Sam Krucker's clean tutorial page at http://plasma2.ssl.berkeley.edu/~krucker/hessi/clean.html.
2. Mem Sato is an optional image reconstruction algorithm utilizing the maximum entropy method described by Sato, Takeo, and Kazuo 1999, PASPJ, 51, 127. In general, mem sato gives very sharp images, but is a little more time consuming to run. This is especially true is if all the collimators are selected for the imaging. The finest four collimators (0 - 3) require substantial computing time since they have the narrowest point-spread functions (PSFs). Eliminating these collimators from the calculation will cause a slight loss of resolution, but will dramatically increase the speed of the image reconhe object syntax:
IDL> im = obj_im -> getdata(/plot,image_algorithm='memsato')
The resulting image should look like Figure 3 below.
Figure 3: A mem sato reconstruction of the default RHESSI simulation.
In order to specify which collimators are used to reconstruct the image, the det_index_mask keyword is used. For example, if you wish to use mem sato to construct an image of the default model using only collimators 5 - 9 (indices 4 - 8), you would type:
IDL> hessi_image,im,pars,/plot,image_algorithm='memsato',det_index_mask=[0,0,0,0,1,1,1,1,1]
Using the object syntax:
IDL> im = obj_im -> getdata(/plot,image_algorithm='memsato',det_index_mask=[0,0,0,0,1,1,1,1,1])
3. Mem Vis is a similar reconstruction algorithm to mem sato except that it operates on visibility data instead of the light curves of the individual collimators. The basic idea behind visibilities is that the modulation profiles for each collimator can be represented by a series of amplitudes and phases by fitting sine functions to the modulation profiles. In some sense, this provides a smoothing to the data before invoking the maximum entropy routine. To reconstruct an image using mem vis from the default RHESSI simulation data, simply type (assuming you have not already changes some of the imaging parameters):
IDL> hessi_image,im,pars,/plot,image_algorithm='vis'
Using the object syntax:
IDL> im = obj_im -> getdata(image_algorithm='vis',/plot)
The resulting image should look like Figure 4 below.
Figure 4: A mem vis reconstruction of the default RHESSI simulation.
To execute the above mem sato example in which a new detector index mask is specified, but this time using the mem vis imaging algorithm, type the following:
IDL> hessi_image,im,pars,/plot,image_algorithm='vis',det_index_mask=[0,0,0,0,1,1,1,1,1]
Using the object syntax:
IDL> im = obj_im -> getdata(/plot,image_algorithm='vis',det_index_mask=[0,0,0,0,1,1,1,1,1])
4. Pixon based image reconstruction tries to maximize the image resolution in locations where the information content of the data is high, but permits the resolution to be low where there is little information. 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. As computing power continues to increase, image reconstruction times will go down; however, at this writing on a 750 MHz Pentium III PC, an image reconstruction of the default RHESSI simulation takes approximately 3.5 hours. On the other hand, you get what you pay for: the Pixon algorithm is very robust and produces images which are photometrically as accurate as possible. Hence it is ideally suited for producing very high quality ("publication quality") images. It is not well suited for exploring a data set. Typically, you will use one of the other reconstruction algorithms to decide what is in your data set and will only use the pixon algorithm to get the final product. It is strongly suggested that users utilize 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. More information on the pixon image reconstruction algorithm and hints on using it can be found at http://www.lmsal.com/~metcalf/HESSI/.
To make a reconstruction of the default RHESSI simulation and plot it, you simply type:
IDL> hessi_image,im,/plot,image_algorithm='pixon'
Using the object syntax:
IDL> obj_im = hsi_image()
IDL> im = obj_im -> getdata(/plot,image_algorithm='pixon')
Again, the first command of the object syntax is not required if the object has already been defined, but if the user has already defined the image object and changed any settings (energy band, time range, etc.) the default simulation image will not be the one operated on. Re-running this first command can reset the image to the default if this is desired.
The resulting image from the pixon algorithm is:
Figure 5: reconstruction of the default RHESSI simulation using the
pixon
imaging algorithm.
The pixon specific optional keywords accepted by hessi_image are:
pixon_guess - (default guess is computed by the quick look software)
The initial guess for reconstructed image. If guess is not nx by
ny or if all elements of guess are LE 0, then the guess is not used.
pixon_noplot - (default=0, not set) If set, don't do any intermediate
plotting.
pixon_progress_bar - (default=1, set) Set to show the progress
bar.
pixon_resolution - (default = 1) Minimum resolution in pixels.
pixon_sensitivity - (default = 1.0) The sensitivity in creating
the pixon map. 0.0 is most sensitive, >3.0 is very insensitive. The
default was arrived at after much testing and simulation. Don't mess with
it unless you know what you are doing. If sensitivity is set too
small, you may get spurious sources; too large and you may miss real sources.
pixon_sizes - (default = powers of sqrt(2.0)) An integer list
of pixon sizes (resolutions) to use. Slower, but more robust, would
be something like pixon_sizes=indgen(17).
pixon_smpattwritedir - (string variable specifying directory,
default not set) If new smoothed patterns are calculated and this keyword
is set to a valid directory, then the new smoothed patterns will be saved
in that directory.
pixon_snr - (default = 0, can be set to 1) Weight the residuals
by the signal-to-noise ratio (S/N).
pixon_error - (a return array) Outputs the estimated uncertainty
on the reconstructed image.
pixon_iterate - (a return variable) Returns with the number
of iterations used at the final resolution.
pixonmap - (a return variable) The pixon map giving the spatial
scale used at each pixel.
pixon_outresolution - (a return variable) Returns with the final
resolution, i.e. the minimum pixon size.
pixon_residual - (a return array) The data residuals
pixon_rgof - (a return variable) The goodness of fit given by
the c-statistic.
5. Forward Fitting assumes that an image can be represented by a number of superimposed 2D gaussians, which can be either spherical, elliptical, or curved. The algorithm fits the modulation patterns calculated from a model map (specified by Gaussian parameters) to the calibrated event list of the telemetry file. A detailed description of the algorithm and its testing can be found in Aschwanden et al. 2002 (ApJ, in preparation), which can be downloaded from http://www.lmsal.com/~aschwand/eprints/2000_hessi.ps. Additional information on forward fitting as well as comparisons with other RHESSI imaging algorithms can be found on pages prepared by Markus Aschwanden at http://hesperia.gsfc.nasa.gov/rhessidatacenter/imaging/forward_fitting.html. An additional tutorial on forward fitting can be found at http://www.lmsal.com/~aschwand/hessi_imaging/hessi_fwd_tutorial1.html. The primary advantages of the forward fitting image 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 (ff_n_gaussians=1). If the resulting chi-square (actually C-statistic, see Cash et al. 1987) comes out near unity, no refinement of the model is justified. However, if the resulting chi-square exceeds the value of 1.1, the user should increase the number of gaussian components (e.g ff_n_gaussians=2), until the model is consistent with a chi-square near unity (say in the range of 0.9-1.1). Current drawbacks are separation of closely-spaced sources, because the convergence behaviour of the forward fitting algorithm is only good if the initial guess of the source positions is accurate, based on the separation capability of local peaks in the initial back-projection map. The forward fitting specific image control keywords along with their default values are:
ff_chi_ff - (return variable whose length depends on the
number of collimators used) Actually the C-statistic value (similar to
chi-square) which results from the fit for each used detector.
ff_chiav_ff - (return variable) C-statistic value averaged from
the used detectors.
ff_coeff_ff - (return variable with variable length depending
on the number of Gaussians and number of parameters per Gaussian used)
This variable contains the fit parameters obtained by the forward fit algorithm.
The largest number of parameters used to fit a single source is likely
to be 7, which fits a curved elliptical Gaussian to the data. Using
less parameters cuts off the total number of parameters, but the ones remaining
are in the same order as if all 7 were there. These fit parameters
are:
par.ff_coeff_ff(0) = amplitude
of source component (normalized to unity)
par.ff_coeff_ff(1) = gaussian
width (in units or arcseconds)
par.ff_coeff_ff(2) = x-position
of gaussian with respect to image center(in arcseconds)
par.ff_coeff_ff(3) = y-position
of gaussian with respect to image center(in arcseconds)
par.ff_coeff_ff(4) = eccentricity,
ecc=(ywidth/xwidth)-1 (for elliptical gaussians)
par.ff_coeff_ff(5) = tilt
angle of ellliptical gaussian (deg counterclockwise from x-axis)
par.ff_coeff_ff(6) = curvature
ratio, i.e. width/curvature radius
The meaning are shown graphically in the following figure:
Figure 6: Graphical representation of forward fitting parameters:
ff_fixed_pixel - (default = 0, can be set to 1) Option to disable
automated pixel adjustement in FWD_MINIMUM
ff_ftol - (default = 1.0e-4) Tolerance value for convergence
of forward-fitting algorithm.
ff_min_sigma - (default = 1.02) Keyword to select detector with
significant modulation.
ff_minsep - (default = 0.0) The minimum separation of sources
(in arcseconds) in gaussian source decmposition.
ff_n_gaussians - (default = 1) The number of gaussian components
to fit.
ff_n_par - (default = 4) The number of free parameters per gaussian,
ff_n_par=4,6,7,...
ff_nitmax - (default = 10) Maximum number of iteration cycles.
If ff_nitmax is set to zero, the input model from /SIMFILE is directly
used.
ff_peak_flux - peak flux value in image in units of [photons
cm-2 s-1 asec-2].
ff_phot_sec - incident photon rate in units of [photons sec-1]
ff_pixelized - (default=1) pixelized (=1) or unpixelized (not
yet available) forwardfit method.
ff_progress_bar - (default = 1) interactive display of progress
bar to monitor speed of iterations. A value of 1 turns it on, 0 turns
it off.
ff_testplot - (default = 1, can be set to 0) Enables testplots
from procedures FWD_DECOMPOSE and FWD_MINIMUM.
You can run the forward fitting algorithm on the default RHESSI image by simply typing:
IDL> hessi_image,im,par,/plot,image_algorithm='forwardfit'
This will invoke the forward fitting algorithm with only one source to be fit, while the default simulated image contains two sources. To use forward fitting with two sources, you can simply type:
IDL> hessi_image,im,par,/plot,image_algorithm='forwardfit',ff_n_gaussians=2
and the resulting image is shown in Figure 7 below.
Figure 7: reconstruction of the default RHESSI simulation using the
forward
fitting imaging algorithm.
This example illustrates one of the main things to keep in mind when using the forward fitting algorithm - the user is responsible for determing the number of sources to include in the reconstructed image. There are of course guidelines for making this decision, supplied by the C-statistic used by the program. If this statistic is near unity, no additional sources are required/justified by the data. Information on the C-statistic plus other details of the forward fitting image are plotted in the test plots generated by the algorithm if the ff_testplot keyword is set to 1 (the default setting). It is recommended that forward fitting be used with ff_testplot activated. The forward fitting algorithm can produce very high resolution images of point sources, as well as other sources which it adequately describes in terms the basis functions it can use to fit to the data. If the actual source is made up of these basis functions, the resulting images will look very good. At the same time, this provides some limitations to the applicability of the algorithm and also requires the user to determine how many sources make up a certain image and what their general characteristics are.
Invoking the forward fitting algorithm using the object syntax is possible with:
IDL> obj_im =hsi_image()
IDL> im=obj_im -> getdata(image_algorithm='forwardfit',/plot,ff_n_gaussians=2)
This call will again simply operate on the default simulated image.
Independent of the imaging algorithm used, there are times when the user my find it desireable to examine some of the intermediate data products which go into a given image. For example, many of the imaging algorithms show their fits to the binned event list, but the user may wish to look at the binned event list separately, to do tests of his or her own devising in trying to decide whether or not to use a given collimator in the image reconstruction. Here, we describe ways to retrieve an plot some of these intermediate data products. Access to most of these products is only through the object oriented sytax.
Binned eventlist: The binned eventlist can be retrieved using the getadata command. If you already have an image object defined, obj_im, then the following:
IDL> be = obj_im -> getdata(class_name='hsi_binned_eventlist')
will return an array of pointers to the binned eventlist.
IDL> help, be
be
POINTER = Array[9,3]
The array of pointers refers to the 9 collimators and the 3 useful harmonics (at this time only the first harmonic is used). These pointers refer to structures which contain the eventlist. For example, to examine the eventlist structure for the first harmonic of the first collimator, simply type:
IDL> help, /structure, *be[0,0]
** Structure HSI_TIME_BIN, 3 tags, length=12:
TIME
LONG
512
COUNT
LONG
16
LIVETIME
FLOAT 1.00000
You can then plot the counts versus time in the binned eventlist for this collimator by typing:
IDL> plot, (*be[0,0]).time, (*be[0,0]).count
RHESSI Imaging Control Keywords
The keywords given in the examples above comprise some of the keywords a user may want to change during image reconstruction; however, there are many additional keywords whose functions are not shown explicitly in the above examples. A full listing of RHESSI imaging keywords can be found at:
http://hesperia.gsfc.nasa.gov/ssw/hessi/doc/hsi_params_image_standard.htm
This list will be automatically updated as the RHESSI software is updated. Some of the commonly used keywords which the user may wish to change when working on imaging data are given below, in alphabetical order.
det_index_mask - CONTROL - a 9 x 3 byte array - It allows specifying the detectors and the harmonics used by the imaging software. det_index_mask[*,0] specifies for the fundamental in each collimator, while det_index_mask[*,2] specifies the third harmonic. Thus, if you wish to use the fundamental for each collimator, then set det_index_mask[*,0] = 1, and all other values to 0. If you wish to use only the fundamental on the first 3 collimators, then set det_index_mask[0:2,0] = 1, and all other values equal to 0. The default setting is for det_index_mask[3:7,0] = 1, and all other values equal to 0.
energy_band - double precision 2 element array - Controls the energy band used to construct the image. Units are in keV, with defaults of [12., 25.].
filename - CONTROL - string - Is used to specify the filename containing the RHESSI telemetry data which will be used to construct the image. This keyword is used when working with simulated or real RHESSI data files, but is not necessary when using the simulation control keywords in hessi_image or the object calls to perform a simulation.
image_algorithm - CONTROL - string - Is used to set what imaging algorithm is to be used in the image reconstruction. The default is 'back projection'. The currently available algorithms are 'back projection', 'clean', 'mem sato', 'mem vis', 'pixons' and 'forward fitting'. The default setting is 'back projection'.
image_dim - CONTROL - 2 element integer array - Sets the size of the image to construct. Units are in pixels, and it is generally recommended to use square images in the reconstruction. The total size (in arcseconds) of the resulting image is then the size of the image given by image_dim multiplied by the pixel_size described below. The default for image_dim is [64,64].
obs_time_interval - CONTROL - Can be used to select the observation time interval you want to look at using ANYTIM format times. The great advantage of using this keyword to specify observing times is that you do not need to know what file the data is in - the software will find it for you. For example, the following call is acceptable: obs_time_interval = ['2000/09/01 12:00:32', '2000/09/01 14:30:23'].
pixel_size - CONTROL - 2 element floating point array - Sets the size of the pixels used in the image reconstruction. Both an x- and y-size are specified (hence the two elements, x- first), allowing for the use of non-square pixels. Units are in arcseconds and the default is [1.,1.].
time_range - CONTROL - double precision 2 element array - Controls the time range for the data used to construct the image. Units are in seconds from the start time of the data file used or relative to obs_time_interval or can be in absolute units as given above for obs_time_interval. If time_range is given as a fully referenced absolute time in in ANYTIM format (e.g. ['2000/9/1 12:06:10', '2000/9/1 12:06:16'], then the control parameter obs_time_interval is set to time_range. Otherwise, time_range is in seconds relative to obs_time_interval. The default is [0.,4.], ideally corresponding to one full rotation of the spacecraft. In reality, the spacecraft spin period is likely to be slightly different, requiring a setting like [0.,3.943] to specify one full rotation.
xyoffset - CONTROL - 2 element floating point array - Sets the
center of the region to be imaged. Units are in arcseconds from the
telescope axis, which is nominally supposed to be pointing at the disk
center of the Sun. The default pointing is [600.,200.].
Working with the Simulated RHESSI Data
Due to the very rapid development of RHESSI data analysis software after launch, the simulation data file is no longer compatible with the current version of the software, and, therefore should not be used at this time. It is anticipated this situation will be corrected by April 15, 2002. In the meantime, the user may find it very useful to consult the section below on working with early flight data. The software to simulate sources should still work fine; however, so users can still perform their own simulations using the calls illustrated above, or using the hsi_sim_flare routine described below.
When simulating RHESSI data, you must first put the software in simulation mode at the start of you SSWIDL session. To do so, type:
IDL> hsi_switch,/sim
Real RHESSI telemetry data will reside in FITS files. An image can be constructed from such FITS files by specifying a filename in the call to hessi_image. Most likely, all RHESSI data FITS files will reside in a specific data archive directory accessible to the all the users at a given institution. Laptop users may also find it desireable to copy some portion of the data to their machines for analysis when not connected to the network. In either case, the HSI_DATA_ARCHIVE environment variable should be set to point to the directory containing the data (see below for how to set this up). It is also true that FITS data files in the current working directory of IDL can be easily read, or alternatively, the entire filename (including path) can be specified when working the data files. For example, if RHESSI telemetry data resides in the file telemetry.fits (in the current working directory), this data can be analyzed with the call:
IDL> hessi_image,im,filename='telemetry.fits'
Using the object syntax:
IDL> im=obj_im -> getdata(filename='telemetry.fits')
Examples of simulated RHESSI telemetry data files are available on the web at http://hessi.ssl.berkeley.edu/data/, under the link to Test data. FITS telemetry files of simulated RHESSI data are also available from http://bigbox.gsfc.nasa.gov/hessi/test_data. Following either link will take you to a directory containing several FITS files which contain simulated data files and flare list and catalog files generated by the quicklook analysis software.
For most imaging tasks involving data files, one can run their IDL session from the same directory where the data is located or the entire path can be specified for filenames if the IDL session is being run from some other directory. Neither one of these solutions may be the most desireable to the user, so there is the HSI_DATA_ARCHIVE environment variable which can be set to specify the location of the data. There are certain RESSI imaging keywords which require this environment variable to be set as well, including the keyword obs_time_interval. It is strongly recommended that the user follow the instruction on obtaining a local copy of the RHESSI test data by following the instructions at http://hessi.ssl.berkeley.edu/software/test_data_access.html. This link also describes the procedure for setting up the HSI_DATA_ARCHIVE environment variable.
Once, the HSI_DATA_ARCHIVE environment variable is set to point at the
location of the RHESSI test data, you are ready to start working with the
simulated data. When real RHESSI data becomes available, the HSI_DATA_ARCHIVE
environment variable should of course point to the directory containing
the real data; however, the user interaction with this data should be identical
to that described below for the simulted data. More complete documentation
on reading and working with the observing summaries and flare lists can
be found at:
http://sprg.ssl.berkeley.edu/~jimm/hessi/hsi_obs_summary.html
http://sprg.ssl.berkeley.edu/~jimm/hessi/hsi_flare_list.html
http://sprg.ssl.berkeley.edu/~jimm/hessi/hsi_filedb.html
Here we describe briefly how to use these files to find information
on flares which may be interesting to image. The use of this information
is not directly integrated into the hessi_image routine calls, so
much of what is described below necessarily uses the object syntax.
The observing summary data are obtained directly from the telemetry packets and the files allow the user to look at RHESSI coutrates and other information before accessing the data files. However, these data are must first be generated at the science operations center (SOC) in Berkeley when the data is downloaded. At present this is not being done, so the observing summary data is not available for the real RHESSI flight data.
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 object methods. First, you need to create the observing summary object:
IDL> obs_summ = hsi_obs_summary()
To get observing summary, you can specify a specific time interval or get the summary data for all the files in HSI_DATA_ARCHIVE directory. The simplest thing to do (and useful for seeing when flares occur and how large they are) is to plot the count rate. This can be simply done by calling the plot method:
IDL> obs_summ -> plot
Several diffent types of summary data can be plotted. Examples are:
IDL> obs_summ -> plot, /countrate
IDL> obs_summ -> plot, /particle_rate
IDL> obs_summ -> plot, /ephemeris
IDL> obs_summ -> plot, /mod_variance
To specify a specific time interval, use the keyword obs_time_interval as you would when imaging. Thus, to return the observing summary count data in a specific interval, you could type:
IDL> time_interval = ['01-Sep-00 12:00', '01-Sep-00
14:30']
IDL> counts = obs_summ -> getdata(obs_time_interval=time_interval)
You can also use the list method to get a text readout of the observing summary data. For example, typing
IDL> text = obs_summ -> list()
will open up a new window giving the observing summary data as textual information. 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)
IDL> obs_summ -> plot
;default to countrates
IDL> 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.
IDL> text = obs_summ -> list()
;default to countrates
IDL> text = obs_summ -> list(/particle_rate)
IDL> text = obs_summ -> list(/ephemeris)
IDL> text = obs_summ -> list (/show, /print, /file_text)
IDL> 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:
IDL> counts = obs_summ->getdata()
IDL> times = anytim(obs_summ->getdata(/time),/ints)
IDL> utplot,times,counts(2,*)
As a routine part of the quick look data products that will be produced by the RHESSI mission, summary files of flare lists will be available with all the data distributed by the mission. In the directory of test data referenced above, this file is hsi_flare_list.fits and this file is a FITS file with binary table extensions. This file can be read at the command line using the routine hsi_read_flarelist. As long as the HSI_DATA_ARCHIVE environment variable is defined, you can read in the flare list by simply typing:
IDL> flarelist = hsi_read_flarelist()
This routine will automI_FLARELISTDATA Array[10]
Then, for example, to examine the 7th flare (which has the highest peak count rate) recorded in the flare list you would type:
IDL> help,/structure,flarelist(6)
** Structure HSI_FLARELISTDATA_EXT, 16 tags,
length=144:
ID_NUMBER
LONG
7
START_TIME
DOUBLE 6.8377184e+008
END_TIME
DOUBLE 6.8377188e+008
PEAK_TIME
DOUBLE 6.8377186e+008
BCK_TIME
DOUBLE Array[2]
ENERGY_HI
FLOAT Array[2]
PEAK_COUNTRATE FLOAT
1174.44
BCK_COUNTRATE FLOAT
1.37913
TOTAL_COUNTS FLOAT
129417.
POSITION
FLOAT Array[2]
FILENAME
STRING 'hsi_20000901_000000_000.fits'
FLAGS
BYTE Array[32]
SEG_INDEX_MASK BYTE
Array[18]
ACTIVE_REGION INT
0
X_POSITION
FLOAT -787.500
Y_POSITION
FLOAT 297.500
which gives you much of the important information for this flare including what file the data for this event is in (under the tagname FILENAME), and the start, end, and peak times given in ANYTIM format, i.e., seconds from 1-Jan-1979 0:00 UT. Also given is the peak count rate for this flare in counts per second and the position of the flare.
In order to make an image of this flare, you need to extract several pieces of information about the flare, including its position and the time of the flare relative to the start time of the data file. The position of the flare is given by the x_position and y_position tag names in the structure, as well as in the position tagname. To use the position tagname, you can print out the position by typing:
IDL> print, flarelist(6).position
which returns the coordinates [-787.500, 297.500]. The start time of the data file is determined from its name. The file name format is hsi_yyyymmdd_hhmm_000.fits, and begins at 0 seconds on the day specified. Thus, the start time of the file containing the 7th flare in the list (flarelist(6)) is 00:00:00 UT on 09/01/2000. To convert this to ANYTIM units, you can type:
IDL> print,anytim('0:00 00/09/01')
6.8376960e+008
Then, if you want to make an image starting at the peak of the above flare, this event begins
IDL> print, 6.8377186d+08 - 6.8376960d+08
2260.0000
seconds after the begining of the data file. NOTE: The times listed in the RHESSI flare lists and other paramter variables require the use of double precision arithmetic. Unfortunately, when you print the variables to the screen, IDL uses an "e" instead of "d" even if the variable is double precision. This may cause problems if you simply try to cut and paste numbers to the command line. Cut and paste followed by a manual change of the "e" to a "d" will work fine though. Thus, you can image the first 4 seconds after the peak of the flare by typing
IDL> hessi_image,im,par,filename='hsi_20000901_000000_000.fits',xyoffset=[-787.50,
297.50],$
IDL> time_range=[2260.,2264.],/plot
Using the object syntax:
IDL> obj_im = hsi_image()
IDL> im=obj_im -> getdata(filename='hsi_20000901_000000_000.fits',
$
IDL> xyoffset=[-787.50, 297.50],time_range=[2260.,2264.],/plot)
In the reconstructed image, the source should appear slightly below and to the right of the center of the image if you are using the default imaging parameters (size, pixel scale, etc.). The source is not exactly centered because of the course resolution used to to find the sources during the creation of the observing summaries.
You can of course use the flare list information to search for particular types of flares that might interest you. For example, you can find the flare with the largest peak countrate using the following commands:
IDL> print,max(flarelist.peak_countrate,iwhr)
6223.36.0
IDL> print,iwhr
1
Thus, the 2nd flare (index 1) in the list has the largest peak countrate. Another example would be to find all the flares in this list that have a peak count rate greater than 2,000 cps. You can find these by typing:
IDL> iwhr = where(flarelist.peak_countrate gt
2000.)
IDL> print, iwhr
1 7
Thus, the only two flares that match this condition are flarelist(1) and flarelist(7).
Working With Early RHESSI Flight Data
At this current early stage in the RHESSI mission, there are no RHESSI flare lists and the observing summary data are not being produced yet. As a result, it is a little more challenging to find flares that may be of interest. One natural thing a user may want to do though is to look at RHESSI data for a flare which has been selected on some other criteria, say GOES X-ray flux for example. The user can then look at the RHESSI light curve for the time of interest, and start to make images. In terms of RHEESI data analysis, the main things the user must select are the time interval, the spacial location and size of the image, and the energy band to analyze. To begin then, we first need to find the time and location of a flare on interest. Another issue to keep in mind is whether or not RHESSI was in the night portion of its orbit, possibly missing some protion of the flare. A useful tool is the summary plots of GOES activity and RHESSI orbit parameters available from the RHESSI data center webpage. Go to We will http://rhessidatacenter.ssl.berkeley.edu/. Under the heading Summary Plots click on the GOES/RHESSI link which will let you navigate to a specific day of interest. For example, there several strong GOES flares on Februart 20, 2002. Following the above link to this day should give the summary plot shown below.
This plot shows that for some flares RHESSI was in its nighttime portion of the orbit at the begining of the flare, such as in the one at 17:00 UT. A nice strong GOES M class flared occured at about 10:00 UT, and two flares went off during its decay around 11:10 and 11:40 UT. We will use these events as our example.
At the current time, the user must specify a time interval or data filename to use. Currently, the user must have cross mounted access to the data at one of the RHESSI data centers, or they should copy the data files of interst to their local machine. Access to the flight data is described in more detail in the Accessing RHESSI Data memo written by Kim Tolbert. For the case we are going to explore here, the user needs to have in their RHESSI flight data directory the files: hsi_20020220_090400_002.fits, hsi_20020220_095720_002.fits, hsi_20020220_103540_002.fits, hsi_20020220_104040_001.fits, and hsi_20020220_114920_001.fits. Warning: in the data analysis software, a file database is maintained which specifies which files contain data for which times. This filename includes the version number (the _004) in the above filenames. If you are downloading the level 0 fits files to your own machine, it is possible the database and your files will get out of sequence on the version numbers. If this happens, the software will return an error message saying no file with the specified time interval is found, and a dialog box should pop up to look for the file. At this point you can select the previous version file that you have on your machine, or download the most current data and try it again.
Using the Lightcurve Object to Find Flares
Making light curves with RHESSI data requires using the object oriented calls to access and plot the data. This is primarily because the lightcurve object is a separate object from the image object in the data analysis software, and only the image object has a more standard IDL interface routine, hessi_image. The first thing to do then is to declare a light curve object:
IDL> obj_ltc = hsi_lightcurve()
In the current absence of observing summary and flare list data, the user must specify an observing time interval or a filename to work with. It is generally more useful to work with observing times, so we will do so here. For the example we are working with, you could then decide to see what was detected by RHESSI for the above mentioned flares between 9:45 UT and 12:00 UT. To get a lightcurve plot, you could simply type:
IDL> obj_ltc -> plot,obs_time_interval=['2002/03/13 09:45','2002/03/13 12:00']
This would create a light curve plot covering the full range of energies recorded by RHESSI: 3.0 - 15000 keV. RHESSI appears to be very sensitive to the low energy, thermal X-rays from flares, which can dominate the behavior of the light curves. The resulting light curve is shown below:
While the M class GOES event shows strong counts in RHESSI, the event just after 11:00 UT is stronger with a very impulsive spike. We will focus on this strong spike. In order to betterr refine the time range which we will image, it is useful to focus in on this spike. Let's create a new plot with the time range between 11:00 and 11:10 UT. In addition, the user may wish to restrict the energy range over which the light curve is plotted in order to focus on the behavior of the hard X-rays in relation to the onset of the flare. This can be accomplised by defining by defining one or more energy bands to be used in constructing the light curve. This is done using the ltc_energy_band keyword. To define a single energy band from 25-50 keV, simply execute:
IDL> obj_ltc -> plot,obs_time_interval=['2002/02/20
11:00','2002/02/20 11:10'], $
ltc_energy_band=[25,50]
Several energy bands can be defined and computed in the same call. This is particularly easy if the desired energy bands are contiguous. For instance, the user may way to define and plot the light curve for the bands: 25-50 keV, 50-100 keV, and 100-150 keV. The call would then be:
IDL> obj_ltc -> plotman,obs_time_interval=['2002/02/20
11:00','2002/02/20 11:10'], $
ltc_energy_band=[25,50,100,150],ltc_time_resolution=1.0
Notice that in this case, we have used the plotman utility instead of the plot utility. Both can be used at anytime; however, plotman makes certain data visualization and hardcopy/image creation simpler. We have also use the ltc_time_resolution keyword to increase the time binning for the light curve from its default value of 0.1 sec to 1.0 sec. The resulting light curve plot is:
The three energy bands are now shown in different colors using the plotman utility, while the plot utility shows the different energy bands in different line types. It looks like most of the action occurs just after 11:04 UT, so lets zoom in between 11:06 and 11:08 with the following call:
IDL> obj_ltc -> plotman,obs_time_interval=['2002/02/20
11:04','2002/02/20 11:08'], $
ltc_energy_band=[25,50,100,150],ltc_time_resolution=1.0
The resulting light curve is shown below:
In this case the counts are shown in a logrithmic scale. This was produced easily in the plotman utility by clicking on "Plot_Control" and selecting "XY Plot Display Options..." and then clicking on the "Log" button in the Y-axis section. The peak of this event lasts between about 11:06:10 to 11:06:30, so we will make our first image in this time interval.
The first thing that need to be done after settling on a time interval is to locate the flare so we can make more detailed images. First, we declare an instance of the image object, and we set our desired time interval up at the same time:
IDL> obj_im=hsi_image(obs_time_interval=['2002/02/20 11:06:10','2002/02/20 11:07:00'])
This selected time interval will now be used for all the images we make until we select another interval. In this case we have selected a longer interval, but this is really just a reference interval for the imaging object. Under its default behavior, the data analysis software makes an image by integrating counts for 4 seconds (the nominal spin period of RHESSI) starting at the begining of the file if a file is set, or starting at the begining of the obs_time_interval if that is set. At this time, it is strongly recommended that images be made over time intervals which are an integral number times the spin period of the satellite. Otherwise errors in the preliminary phase calibration may generate oscillating source locations.We are going to want to integrate over 5 spins (~20 seconds), but as we will see the actual spin period of RHESSI at this time was a little longer than 4 seconds, requiring us to specify more than our 20 second observing time interval. In general, it is a good idea to specify the begining of obs_time_interval at the time you want your interval(s) to begin, and then pad the ending time beyond what you think you will be interested in to account for variations in the spin period of the satellite. Next, we must find the spin period of the satellite, so we can tell the imaging software what value to use. In the future, this will be done automatically, but for now the user must use the PMT RAS software to determine the spin period. Execute the following (you may want to see the next paragraph first and test whether the new PMTRAS software is working before executing the following command):
IDL> pmtras_analysis,'/data/cmj/hessi_data/flight_data/hsi_20020220_104040_002.fits',/no_starid
Upgraded routines to generate a PMTRAS absolute roll solution have recently been uploaded to the development area on SSW. The upgraded routines have two advantages over current procedures. First, they support generation of maps without having to run pmtras_analysis beforehand and second, the resulting maps will (sometimes!) have the correct roll calibration.
To activate the new feature, include the following command:
IDL> obj_im -> set, as_roll_solution='PMT', ras_time_extension=[-500,500]
Note that the first argument given in the call to pmtras_analysis, the filename of the level 0 fits data file for the time interval of interest, must include the full path of the filename on whatever machine you are using. Since this path is almost certainly different for the user, the above command can not simply be cut-and-pasted into your SSWIDL session - some editing will be required for the filename. The program will generate a plot of spin period vs time as determined by the PMTRAS for times covered by this file. You can read the spin period (to ~0.1 millisecond) from the plot, for the time at which you want to image. Alternatively for short files (< 5 minutes), you can use the printed average displayed on the screen. This program is still under development so please file a bud report if it does not work for the file you want to use.
IDL> obj_im -> set,energy_band=[25,50],image_dim=[64,64],pixel_size=[32,32],
$
IDL>
time_range=[0,21.6645],xyoffset=[0,0],det_index_mask=[0,0,0,0,1,1,1,1,0],
$
IDL>
as_spin_period=4.3329
IDL> im = obj_im -> getdata()
IDL> obj_im -> plotman
The resulting image produced in the plotman
utility is:
Our large field of view shows the outline of the Sun in red, and a strong X-ray source at the lower left, near the solar limb. The radial source locations from sun-center should be quite good; however, the azimuthal source locations are totally indeterminant at present. Two-dimensional source locations are, however, self-consistent over intervals of a few minutes. Absolute positions will be available with improvements in the RAS or PMTRAS software over the next week or two.
Now that we have found the location of the image (at about -680, -680 offset from Sun center), we can zoom in on it an get a higher resolution picture. Recall, that we only need to use the set utility to change parameters. We can zoom in by reducing the pixel size, and at the same time we will change which collimators we are using to get better spatial resolution (see the hints below on some strategies to consider when making images with RHESSI):
IDL> obj_im -> set,image_dim=[64,64],pixel_size=[4,4],xyoffset=[-680,-680],
$
IDL> det_index_mask=[0,0,1,1,1,1,1,0,0]
IDL> im = obj_im -> getdata()
IDL> obj_im -> plotman
The new image is shown below, revealing two compact
sources of X-rays near the solar limb:
Now that we have focussed in on the flare site, we can analyze the data in more detail and make higher resolution images. Recall, that in the next call to set, only things that you want to change need be specified. At this time the phase calibration of the finest 2 grids is not completed, so we will use all the other collimators but those. As described in the next section, collimator 3 can achieve 7 arcsecond imaging, so it is desireable to have slightly smaller pixels than the 4 X 4 arcsecond ones used before. We also refine the offset pointing some as well, and select the clean imaging algorithm:
IDL> obj_im -> set,pixel_size=[2,2],xyoffset=[-680,-650],det_index_mask=[0,0,1,1,1,1,1,1,1],
$
IDL> image_algorithm='clean',clean_niter=100
IDL> im = obj_im -> getdata()
IDL> obj_im -> plotman
Note that we also increased the maximum number of clean iterations allowed, since the default of 33 is often reached by the imaging software. The resulting plot is:
The two sources are clearly seen to be of roughly equal strength and separated by about 30 arcseconds.
If we then change the energy band to 10 - 15 keV, we can look at the X-ray morphology at lower energy:
IDL> obj_im -> set, energy_band=[10.,15.]
IDL> im = obj_im -> getdata()
IDL> obj_im -> plotman
Whit the result:
The two sources have now merged into one larger source, but with the lower right portion noticeably stronger.
We can also examine the image produced by other algorithms. For example, switching back to 25 - 50 keV and using the forward fitting algorithm, we get:
IDL> obj_im -> set, image_algorithm='forwardfit',ff_n_gaussians=2,energy_band=[25.,50.]
giving:
Note that we had to specify two Gaussians for the forward fitting algorithm to distinguish the two sources. This image looks somewhat different than that produced by clean. Here, the two sources appear to be of significantly different width with the one in the lower right having a stronger peak. This illustrates that some care must be execised when evaluating RHESSI imaging. Shifting back to the lower energy band for the, the forward fitting algorithm gives:
IDL> obj_im -> set,energy_band=[10.,15.]
The resulting image again shows two symmetric gaussians because that is the number of parameters we allowed in the forward fit. One can use more paramters to fit the data, to get asymmetric Gaussians and other shapes; however, before we pursue this, let' return to the 25 - 50 keV image. The two algorithms used above gave somewhat different images, so let's try the mem sato algorithm and see if it favors one or the other:
IDL> obj_im -> set,energy_band=[25.,50.],image_algorithm='memsato'
IDL> im = obj_im -> getdata()
IDL> obj_im -> plotman
The resulting image for the 25 - 50 keV band using mem sato is:
This image again shows 2 very similar point sources. This image shows a weak favoring for the source in the lower right as stronger, but shows no significant difference in the widths of the two sources. Swithching to the lower energy band, we get
IDL> obj_im -> set,energy_band=[10.,15.]
IDL> im = obj_im -> getdata()
IDL> obj_im -> plotman
and the image is:
The resulting imahe potentially breaks up into 3 sources, but overall resembles the clean reconstruction above. Other algortihms can be used to explore the flare morphology, and as the finer collimators become well calibrated (and all the calibrations improve) the RHESSI imaging will stabalize. However, at this time flares can be imaged and the user can readily look for interesting changes in the flare morphology as a function of energy and time.
Some Hints on Analyzing RHESSI Images
It is inportant to keep in mind some of the peculiarities associated with RHESSI data when contructing images and analyzing the data. At first, it might seem that the best results are obtained using all the subcollimators since this results in the largest number of total counts available for constructing the image; however, this is generally not the case. Subcollimators with a FWHM resolution much greater than the source extent will only indicate integrated flux, but will not actually increase the signal-to-noise (S/N) of the reconstructed image. Recall that the FWHM resolution of the collimators are 1 - 2."3, 2 - 3."9, 3 - 6."9, 4 - 11."9, 5 - 20."7, 6 - 35."9, 7 - 62."1, 8 - 107", and 9 - 186". Subcollimators with a FWHM resolution greater than or equal to the imaging field of view are not useful for the backprojection algorithm, or any other algorithm which relies on backprojection as the basis for the imaging. Thus, when imaging relatively small regions of a flare, collimators 8 and 9 provide no information on the morphology of the sources. Subcollimators with a FWHM resolution less than the diameter of the smallest feature in the image also add no information on the source morphology but do contribute noise to the map. Subcollimators with a resolution less than 2 times the pixel size will introduce alias features into the map, and should never be used. Subcollimators that are transparent at the energy band selected for the map only indicate integrated flux and add no information on source morphology. These considerations suggest that weighting among selected subcollimators could be very useful, but this has not yet been implemented. Weighting of different subcollimators can certainly help backprojection and CLEAN image quality. For example, some different weighting schemes and their usefulness are the following: a.) Equal weighting of selected subcollimators (so called natural weighting) is optimum for source detection. b.) Weighting that scales with 1/resolution eliminates wide ?skirts? to PSF and will sharpen up images of narrow point sources relative to natural weighting. c.) Weighting that partially suppresses fine subcollimators (tapering) gives a smoother PSF. Some variant of these three weighting schemes will be added to the software in the near future.
Aside from the choice of subcollimators which will optimize a given image, there are spectral considerations for the imaging to be aware of as well. In some ways, it is difficult to separate the spectral ascpects of the imaging from the imaging itself. This leads to a few facts the user should be aware of. The first consideration then is imaging spectroscopy. Many users will want to make images of fields which contain more than one source and then wish to obtain a spectrum of each source. The current way to approach this task is to make images in a number of different energy ranges and then add up the counts in each source to construct the spectrum of the sources.
A related consideration then is just how broad an energy range to use when constructing an image. One might assume that the broader the energy range the better since you get more counts and hence higher S/N; however, this is not necessarily the best choice. The RHESSI imaging software uses all the photons detected between the specified energy range to determine the modulated light curve from each detector for the specified time interval. The software then assumes a mean energy for this energy range for computing the expected response of the instrument which is applied to the data to compute the image. However, since the modulation patternes themselves are a function of energy, there is no unique instrument response that works for all energies. Therefore, it is not a good idea to construct images with an arbitrarily wide energy range. One way to attack this problem is to make several images with smaller energy ranges (get an estimate from Gordon for this) and compare them. As long as the structure of these various images remains the same, it is probably OK to then make a final image encompassing the whole energy range.
Another important aspect of imaging analysis with RHESSI data is to
be able to decide which features in the image you can trust. Some
factors to consider when trying to judge which features of an image to
trust are the following:
A. Counting statistics - Statistical S/N in the peak of a backprojection
map is equal to k * SQRT(total counts)
B. Image morphology:
- Complexity of image (simpler is better)
- Fraction of total counts in feature of interest (a smaller
fraction is more susceptible to systematic calibration errors)
C. Special circumstances - basically there are various circumstances
which might create spurious signals in the image:
- Importance of dead time - does something in the image
change as the dead time increases?
- Fast-rate mode - does a new feature appear when the
instrument goes into fast-rate mode?
- Attenuator changes - does a new feature coincide with
a change in the attenuator position?
- Special source locations (relative to pointing axis
or imaging field of view)
- Suspicious symmetries (arcs relative to pointing axis
or strongfield of view which still encompasses the source
3. A different pixel size
4. A different image center
&n