An Overview of the Command Line Interface for RHESSI Data Analysis Software
DISCLAIMER: This document is still under development, as is the analysis software itself. Incomplete descriptions 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://mssly5.mssl.ucl.ac.uk/surf/guides/hag/ provides a description of the RHESSI satellite and its general scientific capabilities. 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/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 desirable 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/params/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/params/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://hesperia.gsfc.nasa.gov/hessidata/test/. 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
5. A different combination of subcollimators
6. A different energy range
7. A different time range
8. A different set of calibration parameters
9. Odd/even half-rotations
B. Apply redundancy tests:
1. Compare spectra/lightcurves from different detectors
2. Check rotational symmetry of light curve
C. Trace back the feature to:
1. Fourier components
2. Observed modulation curves
These considerations also reveal the importance of being able to simulate RHESSI data. 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 the user expects it should. Below, we describe how to get started with this.
Saving and Recovering Images and Writing and Reading FITS Files
Once you have created an image that you are happy with, you can save this image any number of ways. The IDL save and restore commands provide an easy means for doing this. It is also possible to write a FITS file of the image as well. If you are using the object syntax, simply call the method fitswrite. For example, if your imaging object is called im_obj and you have created an image you are happy with, simply type
IDL> im_obj -> fitswrite
and the image will be written in a FITS file in the current directory with the name hsi_image_thedate.fits. If you are using the hessi_image wrapper for the imaging analysis, you can use the Astron library routine writefits assuming it is on your path (should be in a standard SSW installation). If you have used hessi_image to create an image called im that you are happy with, simply type
IDL> writefits,'filename.fits',im
where the actual file name you choose is up to you.
If you have FITS files of saved images, you can read the image back into IDL using the hsi_image_fitsread() command. Executing:
IDL> im = hsi_image_fitsread()
will open up a dialog box in which you can select the image to read. Once read, the variable im will contain the image which was stored in the file. You can also read the FITS file into a RHESSI imaging object. To do so, simply type:
IDL> im_obj = hsi_image_fitsread(/object)
NOTE: When using this command to load a FITS image into an object, you will need to have access to the original data from which the image was constructed in order to use the object. This is because the eventlist data is not stored in the FITS file. Therefore, if you wish to act on this object, for example by changing a parameter and re-imaging, the software will have to go back to the original data to execute.
As described above, images can be constructed from RHESSI telemetry data files, or the user can use hessi_image to actually simulate RHESSI data and then analyze it. A more sophisticated simulation tool, hsi_sim_flare, is under development; however, much can be done with hessi_image. The default model consists of two point sources (one centered, one in the lower left portion of the image), assuming a photon rate of 1000 photons per second per collimator. The simulation control keywords all have the prefix sim_ to set them apart. Most models of RHESSI data use Gaussian sources of various sizes and aspect ratios. Curved Gaussian sources can also be used to simulate flare loops. The suite of RHESSI software has a convenient means built in to simulate Gaussian sources. One must specify the source amplitude (normalized to unit maximum - the overall flux normalazation is handled in the number of photons specified), location with respect to the center of the simulation box (note: this is a relative position with respect to the center of the simulation box specified by the sim_xyoffset keyword), width in the X and Y directions, and the tilt angle. In the case of curved elliptical Gaussians, a curvature parameter must also be specified. This information is entered into a source structure for use by the software. For example, suppose you wish to specify a single Gaussian source offset +10 arcseconds in the X-direction and +2 arcseconds in the Y-direction from the center of the image, with a X width (sigma) of 5.0 arcseconds and a Y width of 2.0 arcseconds, and a tilt angle of 45.0 degrees, the following commands should be entered:
IDL> gg = {gaussian_source_str}
IDL> gg.amplitude = 1.0
IDL> gg.xysigma = [5.,2.]
IDL> gg.xypos = [10.,2.]
IDL> gg.tilt_angle_deg = 45.
IDL> hessi_image,im,par,/plot,sim_model=gg
or the last command can be replaced by the object syntax:
IDL> obj_im = hsi_image()
IDL> im = obj_im -> getdata(/plot,sim_model=gg)
Note that the keyword specifying the model does not have the leading sim_ in it. This is a bug and will hopefully be fixed soon.
You can specify multiple Gaussian sources as well. Suppose you wish to specify 2 sources, both with widths of [5.0,2.0], with one located at an offset [20.0,-3.0] and one half as strong as this at [-5.0,15.0], both with 0.0 tilt angles, you would type:
IDL> gg1 = {gaussian_source_str}
IDL> gg1.amplitude = 1.0
IDL> gg1.xysigma = [5.,2.]
IDL> gg1.xypos = [20.,-3.]
IDL> gg1.tilt_angle_deg = 0.
IDL> gg2 = {gaussian_source_str}
IDL> gg2.amplitude = 0.5
IDL> gg2.xysigma = [5.,2.]
IDL> gg2.xypos = [20.,-3.]
IDL> gg2.tilt_angle_deg = 0.
IDL> model = [gg1,gg2]
IDL> hessi_image,im,par,/plot,sim_model=model
again, the object version of the last line is
IDL> im = obj_im -> getdata(/plot,sim_model=model)
Another means of specifying point sources is to simply create a model image loaded with 0s everywhere except for the locations of the point sources, which should be loaded with a normalized intensity. For example, to create a 64 X 64 image with 3 point equal strength point sources at positions [10,20], [30,12], and [15,40] the following syntax can be used:
IDL> model=fltarr(64,64)
IDL> model[10,20]=1.0
IDL> model[30,12]=1.0
IDL> model[15,40]=1.0
IDL> hessi_image,im,sim_model=model
with the object version of the last line
IDL> im = obj_im -> getdata(sim_model=model)
The strength of the point sources again is given in normailized units, while the overall flux level is given by the keyword sim_photons_per_coll. In the above example, all three point source are of equal strength. If you wanted to make the one at [30,12] equal to 1/4 of the other two, you would have typed:
IDL> model[30,12]=0.25
instead. The simulated data will be constructed assuming the default of 1000 photons per second per collimator incident on the spacecraft, which is not particularly good statistics, especially when spread over 3 sources. The overall flux level in the simulation could be increased with the call:
IDL> hessi_image,im,sim_model=model,sim_photons_per_coll=5000
Using the object syntax:
IDL> im = obj_im -> getdata(sim_model=model,sim_photons_per_coll=5000)
Modeling Flare Loops:
The RHESSI simulation software can also be used to specify a curved, elliptical Gaussian which is a good proxy for a solar flare loop. Instead of declaring your source as a {gaussian_source_str}, it is declared as a {curved_gaussian_source_str}. In this case there is one additional structure tag (.curvature) which must be specified. For example, suppose you want to declare a curved Gaussian with the same basic parameters as given in the first Gaussian source example above. You would simply type:
IDL> cgg = {curved_gaussian_source_str}
IDL> cgg.amplitude = 1.0
IDL> cgg.xysigma = [10.,2.]
IDL> cgg.xypos = [10.,2.]
IDL> cgg.tilt_angle_deg = 45.
IDL> cgg.curvature = 0.2
This last parameter specifies the curvature by giving the inverse of the radius of curvature of the source, normalized by the larger of the two Gaussian widths. Basically, the radius of curvature of the source is larger of the two Gaussian widths divided by cgg.curvature. In the example above then, the radius of curvature is 25 pixels, which should be quite noticeable in RHESSI data. To then simulate and analyze this model, you can type:
IDL> hessi_image,im,par,/plot,sim_model=cgg
Using the object syntax:
IDL> im = obj_im -> getdata(sim_model=cgg,/plot)
Given this definition for the curvature, you can declare an elliptical Gaussian with no curvature by setting:
IDL> cgg.curvature = 0.0
If you wish to specify a model with a curved, elliptical Gaussian plus some number of point like sources, all components should be specified using the {curved_gaussian_source_str} system, with the point sources given a curvature of 0. For example, if you wish to make a model with the curved source specified above plus a weaker, symmetric point source (not curved) above the loop, you can type:
IDL> cgg1 = {curved_gaussian_source_str}
IDL> cgg1.amplitude = 1.0
IDL> cgg1.xysigma = [10.,2.]
IDL> cgg1.xypos = [10.,2.]
IDL> cgg1.tilt_angle_deg = 45.
IDL> cgg1.curvature = 0.2
IDL> cgg2 = {curved_gaussian_source_str}
IDL> cgg2.amplitude = 0.2
IDL> cgg2.xysigma = [2.,2.]
IDL> cgg2.xypos = [30.,25.]
IDL> cgg2.tilt_angle_deg = 0.
IDL> cgg2.curvature = 0.0
IDL> model = [cgg1,cgg2]
IDL> im = obj_im -> gatdata(/plot,sim_model=model)
Note that the amplitude for cgg2 is much less than for cgg1. This is necessary so that the point source (cgg2) does not overwhelm the loop (cgg1) in the image. This possibility results from the fact that the amplitudes specify the total flux from the given component, so for large sources which are spread out over a substantial area, a larger relative amplitude is required to give them an appreciable maximum in the image. To get a good image of the above source, you will likely have to improve the spatial sampling of the simulation and reconstructed image if you have not already done so. A decent image can be obtained by setting the additional keywords in the call:
IDL> hessi_image,im,par,/plot,sim_model=model,sim_photons=50000,sim_pixel_size=2.0,
$
IDL> pixel_size=2.0,energy_band=[6.,25.],image_alg='clean',niter=100
Using the object syntax:
IDL> im = obj_im -> getdata(/plot,sim_model=model,sim_photons=50000,sim_pixel_size=2.0,
$
IDL> pixel_size=2.0,energy_band=[6.,25.],image_alg='clean',niter=100)
The following is a list of keywords that can be specified to control the simulated data to be produced.
bkgd_photon_energy_ptr - A pointer to an array of photon energies for the background photon flux.
bkgd_photon_flux_ptr - A pointer to an array of photon flux in photons/cm^2/sec/keV for the background spectrum.
bkgd_spec_model - the background spectral model, the default is to use HESSI_BACKGROUND
bkgd_spec_pars - spectral parameters for the given model, the default is to use HESSI_BACKGROUND
sim_model - Can be used to pass in a structure (labelled gaussian_source_str) which specifies a 2d Gaussian source. The Gaussian is specified by its amplitude, x and y standard deviation (in arcseconds), its position relative to the center of the image (in arcseconds), and its tilt angle. See above for examples. The model can also be given as an array with zeroes everywhere except at the location of point sources. The non zero values should be specified in normalized units, with a maximum of 1.0. Alternatively, one can use the keyword model (without the leading sim_) to specify single Gaussian components, controlling the amplitude, x and y standard deviations, and the tilt angle of the Gaussian.
The result will show the image of this Gaussian source and load it into the variable im.
sim_pixel_size - Is used to set the size of the pixels when using a model image of some size (e.g. 64 X 64). The default sim_pixel_sizeis 2.0 arcseconds, which is different from the default imaging pixel_size.
sim_xyoffset - 2 element floating point vector - Is used to specify the center of the simulation box. Note that this value must be reasonably close to the value of the xyoffset imagaing control keyword.
sim_photons_per_coll - floating point number - Is used to control the overall flux level of the simulated data. The units are in photons per second per collimator incident on the spacecraft, with a default value of 1000. This can be an array, one for every a2d_index. Since the imaging algorithms return the photon flux per square cm, to make meaningful tests of the overall photometric accuracy of the various algorithms, you need to know the total effective area of the detectors. This is currently measured to be 39.59 cm2.
sim_a2d_index_mask - 27 element byte array - Controls which collimator's data is used in the simulation. A value of 1 indicates the collimator will be used, and a value of 0 indicates that it will not. The first 9 elements refer to the front segments of collimators 1-9, the next 9 elements of the array refer to the low energy portion of the rear segments of collimators 1-9, and the last 9 elements refer to the high energy portion of the rear segments for collimators 1-9.
sim_time_range - the time range, in seconds, relative to the reference time sim_ut_ref for the simulation
sim_time_unit - the time unit, in binary microseconds of the .time tag in the eventlist structure. The default value is 16.
sim_ut_ref - the reference time for the simulation, this is usually the start time. Can be in any ANYTIM format. The default is still '4-jul-2000 00:00'.
sim_time_profile - a structure with two tags, time and profile The time profile of the generated events follows profile. The probability of acceptance is governed by interpolating using the time tag, referenced to the start of the time_range. Profile has values from 0 to 1.
sim_count_dist - count rate distribution to select channels from at random. The default is to use the default distribution of hsi_cspectrum_dist, which is a power law with a spectral index of 4.0.
sim_background - a factor used to multiply the background spectrum given by HESSI_BACKGROUND. The default is 0, If set to 0 no background is used.
sim_bkgd_time_profile - a structure with two tags, time and profile, for the background spectrum.
sim_bkgd_count_dist - count rate distribution to select channels from at random. The default is to use the default distribution of HESSI_BACKGROUND.
sim_energy_band - Energy band for the simulation, the default is [6.0, 100.0] keV.
sim_max_size - ???
sim_srt_filename - The grid response database file to use for
simulations. The default is
$SSW/hessi/dbase/grid_resp/SRT_default.dat.
sim_atten_state - the attenuator state, 0 = no attanuators, 1
= Thin in, thick out (default), 2 = thick in, thin
out, 3 = both in.
sim_use_spectrum - if set, then use an input spectrum or spectral model to create the simulations. (As opposed to the old method of specifyingsim_photons_per_coll, and a count distribution). If this is set, then the spec_model, spec_pars, photon_flux and photon_energy keywords can be used. If /use_spectrum is set without those keywords, then the default is to use a thermal plus double power law spectrum, 'f_vth_bpow', with parameters a = [1.0, 1.0, 1.0, 3.0, 60.0, 3.0, 10.0, 1.5], where a[0] = emission measure /10^49, a[1] = temperature (kT) in keV, a[2] = flux at 50 keV in photons/cm^2/sec/keV, a[3] = power law index below the energy break, a[4] = break energy in keV, a[5] = power law index above the break energy, a[6] = low energy cutoff of the power law spectrum, and a[7] is the power law index below the low energy cutoff.
sim_spec_model - the spectral model, the default is 'f_vth_bpow', double power law plus thermal emission. These models are defined for the spectroscopy routines used by the XRAY package. The XRAY package must be part of the local SSW tree.
sim_spec_pars - spectral parameters for the given model, usually a fltarr[8], see above for the default values.
sim_photon_energy_ptr - A Pointer to an array of photon energies for the input photon flux.
sim_photon_flux_ptr - A Pointer to an array of photon flux in photons/cm^2/sec/keV. If this is set the spec_pars and spec_model keywords are ignored, but the photon_energy_ptr keyword must also be used.
use_bkgd_spectrum - if set, then use an input background spectrum
for the simulations. If this is set, then the bkgd_spec_model, bkgd_spec_pars,
bkgd_photon_flux and bkgd_photon_energy keywords can be used.
If /use_bkgd_spectrum is set without those keywords, or if sim_background
is set without setting /use_bkgd_spectrum, then the default is to
use a spectrum given by HESSI_BACKGROUND, multiplied by the value passed
in as sim_background.
Using HSI_SIM_FLARE
HSI_SIM_FLARE is designed to simulate a RHESSI flare, and write the data into a FITS file. It can be called with no input parameters, then a default flare is used. It can be called with the /random_flare keyword set. Then a flare with some random time range, and random spectrum, and position is created and output.
It is important to note here that all of the input to HSI_SIM_FLARE is via keyword. And all of the output uses the same keywords (except for the packet array). SO if a keyword is set to a certain value, then that value is used in the program, and passed out. If a keyword is not set, then it is given value and that is passed out in the keyword.
The packets created are written into a file, and passed out as an array.
HSI_SIM_FLARE keywords:
random_flare = if set, HSI_SIM_FLARE does it's own flare, it chooses the spectrum, time profile and xy_offset. It does not change the image model, i.e., the number and relative strength of sources or their position relative to the simulation box.
model = A square map or an array of Gaussian sources. Gaussian Sources
are input as structures,
{GAUSSIAN_SOURCE_STR},
gaussian_sources = {gaussian_source_str,
amplitude: 0.0, relative amplitude
xysigma: fltarr(2), standard deviation in
x and y directions in arcseconds.
xypos: fltarr(2), offset in x and y directions
in arcseconds.
tilt_angle_deg: 0.0} tilt angle in degrees
of x axis for ellipsoid from X. Positive angle moves;
xy_offset = the position of the map center for the model, arcsec from Sun center.
time_array = A time array for the photon flux array, in anytim format.
photon_energy = the photon energy array, in keV. This must be present to use the photon_flux keyword
photon_flux = the photon flux incident on the detectors photons/sec/cm^2/keV. If photon flux is set, then all of the spectral model keywords are ignored. This can be an array of (n_energies, n_times) with n_times obtained from time_array, or (n_energies, 3), using the start peak and end times defined below. Different spectra are allowed for different sources, so if you pass in a model with N sources you can pass in a flux array of (n_energies, n_times or 3, N).
TIME KEYWORDS: Used if time_array is not passed in, if time array is passed in, the obs_time_interval, ut_ref and time_range keywords are set to the appropriate valuse
obs_time_interval = the obs_time_interval to use, for the simulation in anytim format. This will overwrite ut_ref and time_range. The default is [flare_start - 2 minutes, flare_end + 2 minutes]
ut_ref = a reference time for the start of the simulation, in ANYTIM format. The default is the flare start time - 2 minutes
time_range = the time range for the simulation, relative to ut_ref. The default time range is [flare_start - 2 minutes, flare_end + 2 minutes] measured from ut_ref
flare_start_time, flare_end_time, flare_peak_time = start, end, and peak times for the flare, relative to ut_ref.
flare_start_ut, flare_end_ut, flare_peak_ut = start, end, and peak times for the flare, anytim format, this will take precedence over the _time keywords.
SPECTRAL KEYWORDS: only used if photon_flux is not used
spec_model = the spectral model, the default is 'f_vth_bpow', Any of the XRAY package models can be used The spectral keywords are not used or set if photon_flux and photon_energy are passed in
spec_pars = spectral parameters, This can be an array of [8, ntimes] if time_array is passed in, an array of [8, 3] if start, peak and end times are passed in, or an array of [8]. The default is an array of [8, 3], for start peak and end times:
spec_pars[*, 0] = [0.01, 0.1, 0.01, 6.0, 60.0,
8.0, 10.0, 1.5]
spec_pars[*, 1] = [0.50, 2.0, 1.00, 3.0, 60.0,
3.0, 10.0, 1.5]
spec_pars[*, 2] =[0.01, 0.1, 0.01, 6.0, 60.0,
8.0, 10.0, 1.5]
This can be done for multiple sources; if you pass in a model with N sources, then you can pass in an array of [8, n_times or 3, N]
neupert_effect = if set, this will simulate the Neupert effect, by arranging things so that the emission measure is proportional to the time_integrated photon flux at 50 keV. This is only available if the spectral parameters are passed in at 3 times, start, peak and end, or if the random_flare keyword is set.
erange = the photon energy range, the default is [3.0, 500.0], keV.
d_energy = the energy step for the energy array, the default is 1 keV.
BACKGROUND KEYWORDS:
Blevel = background_level, multiplies the background given by HESSI_BACKGROUND, can be a function of time, with the same number of elements as time_array or btime_array. The default is 1.0
btime_array = a time array for the background.
OTHER KEYWORDS:
otp_file = if set, output into this file, the default is to write a FITS file using HSI_ANYTIM_2_FILENAME, with a filename given by hsi_yyyymmdd_hhmm_nnn.fits.
do_catalog = read the output file, do the hessi catalog, and append
ignore_orbit = if set, ignore orbit data for the catalog, This is the default, set it to 0 if you want to deal with the orbit data
all_flare = if set, all data is a flare, short-circuits the usual flare_list process.
score = The Simulated eventlist in an HSI_SCORE structure, with tags:
ut_ref = Reference time in ANYTIM format
time_unit = the time unit for the eventlist
in binary microseconds
Eventlist_ptr = a pointer to an array of {hsi_eventlist}
structures one for each detected photon with
tags:
a2d_index = 0 to 27
channel = 0 to 8193
time = time in units of time_unit binary microseconds
In this way eventlists from more complex sources can be built up by running multiple calls to HSI_SIM_FLARE, concatenating the score structures using HSI_SCORE_CONCAT and outputting using HSI_SCORE2FILE.
EXAMPLES:
Random flare times, spectrum and position, the default 2 source model, Add Neupert Effect. No ut_ref, time_range or obs_time_interval set, so you just get the flare with 2 minutes of data at each end.
IDL> Hsi_sim_flare, packets, /random_flare, /neupert_effect, otp_file = otp_file
you need to get the filename out using the otp_file keyword. If you need any of the parameters used, pass those out using the keywords. e.g., if i have
IDL> Hsi_sim_flare, packets, /random_flare, /neupert_effect,
otp_file = otp_file, photon_flux = f, $
IDL> photon_energy = e, time_array = t
then the photon flux and energy arrays and the time array used in the calculation will be passed out in the variables f, e and t. Be sure that those variables are undefined before you use them, because if photon_flux and energy are passed in, then the /random keyword will not work. Likewise you can get the flare times by passing out variables using the flare_start_ut, flare_peak_ut and flare_end_ut keywords. More examples will follow as testing continues....
Interaction Between the RHESSI GUI and the Command Line Interface
It is possible to work with the same object in both the GUI and on the command line. This is only possible if you are using the object oriented approach on the command line: such switching back and forth is not supported if you are using hessi_image. If you are using the RHESSI GUI and wish to pass data to the command line, use the hessi_data routine. This routine allows you to pass images, lightcurves, and observing summaries using the keywords image=, lc=, or obs= respectively. In all cases though, you are passing the relevant object to the command line. So, after making an image using the GUI, if you wish to pass the image object to the command line, simply type:
IDL> hessi_data, image=im_obj
This makes the imaging object im_obj accessible from the command line. To load the actual image which was last made with the GUI into a variable for further manipulation, plotting, or whatever, you would need to execute the command:
IDL> im = im_obj -> getdata()
It is important to note that in this situation, actions done in either the GUI or on the command line will change the object in both places. For example, let's say the original image produced with the GUI was using the backprojection image algorithm, and you then execute the hessi_data and getdata() calls above. If you then use the GUI to clean the image and execute the getdata() call again (but not the hessi_data call), the cleaned image will be loaded into the variable im. Conversely, if you then used the command line to get a new image using mem sato by typing:
IDL> im = im_obj -> getdata(image_algorithm='mem sato')
these new imaging parameters are also used by the GUI. So, if for example you then hit the "Make and Display Image" button on the imaging control widget of the GUI, it will immediately update the parameters to reflect the use of mem sato and display the image in the plot window without having to recalculate the image.
One can also create an object on the command line and then pass it into the GUI by simply giving the object as an argument to the call to hessi (which activates the GUI). This is useful for example when wanting to simulate RHESSI data, but then use the GUI to perfrom the imaging. For example, you can set up the simulation of a curved Gaussian source and then analyze it using the GUI by executing the following commands:
IDL> imob = hsi_image()
IDL> cgg1 = {curved_gaussian_source_str}
IDL> cgg1.amplitude = 1.0
IDL> cgg1.xysigma = [10.,2.]
IDL> cgg1.xypos = [10.,2.]
IDL> cgg1.tilt_angle_deg = 45.
IDL> cgg1.curvature = 0.2
IDL> model = cgg1
IDL> imob->set, sim_photons = 10000
IDL> imob->set, sim_model = model
IDL> hessi, imob
Unlike imaging with the RHESSI software analysis tools, spectroscopy on the command line is done entirely through object oriented data calls: there is nothing comparable to hessi_image. Therefore, if you have not worked through the above imaging examples and are uncomfortable with object oriented calling of routines, you may find it helpful to first read the imaging section above to get aquainted with this syntax. It is certainly not difficult, but may just be different from what long time IDL users who have not worked with object are familiar with. The basic goal of the RHESSI data analysis spectroscopy software is to produce semi-calibrated data of photons s-1 cm -2 keV-1 in some number of energy bins during some time interval, taking into account the effective area of the detectors, their energy response, and the state of the instrument attenuators. The commands related to RHESSI spectroscopy. Simply type:
IDL> hsi_spectroscopy_list
and a new window will pop up containing a summary of RHESSI commands for spectroscopy.
Unlike the RHESSI imaging software which has an easily built in set of commands for data simulation on the fly, the spectroscopy software expects to work with data files. Therefore, the user will have to work with real data, the simulated data on the RHESSI website, or use the simulation software described above to create and write out simulated data files. Here, we will work with the simulated data files. The first thing that needs to be done is to declare a spectrum object, and to specify which data file you are working with. This can be done with the following:
IDL> sp = hsi_spectrum(filename='D:\hessi\test_data\hsi_20000901_000000_000.fits')
where is this example the full path of the data file on a Windows machine has been specified. Just as is the case with the imaging analysis, there are a large number of parameters related to spectroscopy, some of which are control parameters that the user may ish to set and some of which are just informational. Also as in the case of imaging, these can be examined using the object call to the print method. To look at all spectroscopy related parameters after first defining the spectrum object with the call above, you can type:
IDL> sp -> print
If, on the other hand, the user simply wants to look at only the control parameters, then enter:
IDL> sp -> print, /control_only
The output from this is:
** Structure <b79570>, 51 tags, length=312,
refs=1:
SP_CHAN_BINNING INT
0
SP_CHAN_MIN
INT
0
SP_CHAN_MAX
INT
0
SP_ENERGY_BINNING
INT =
7
SP_TIME_INTERVAL
INT =
60
SP_SEMI_CALIBRATED
BYTE = 0
SP_DATA_STRUCTURE
BYTE = 0
SP_DATA_UNIT STRING
'Counts'
SP_ERROR_OUT BYTE
0
POISSON
BYTE 0
SEED
FLOAT 0.000000
COINCIDENCE_FLAG
INT =
0
OTHER_A2D_INDEX_MASK
BYTE = Array[27]
DET_INDEX_MASK BYTE
Array[9]
SEG_INDEX_MASK BYTE
Array[18]
SUM_FLAG
INT
1
A2D_INDEX_MASK BYTE
Array[27]
COINCID_TIME_DELTA
INT =
1
ENERGY_BAND
DOUBLE Array[2]
NEW_GAIN
INT
1
SASZERO
BYTE 0
TIME_RANGE
DOUBLE Array[2]
TIME_UNIT
INT
16
FILE_TYPE
STRING 'fits'
APP_ID
POINTER <NullPointer>
OBS_TIME_INTERVAL
DOUBLE = Array[2]
PACKET_PER_BUNCH
LONG =
0
FILENAME
STRING 'D:\research\solar\hessi\data\test_data\hsi_20'...
SRM_CT_EDGES INT
2
SRM_PH_EDGES STRING
''
SRM_USE_SEGMENT BYTE
Array[18]
SRM_SEP_DETS BYTE
0
SRM_SEP_SEGS BYTE
0
SRM_USE_SEGCO BYTE
0
SRM_USE_DETCO BYTE
0
SRM_ASTRO
BYTE 0
SRM_DEAD_SEG BYTE
Array[18]
SRM_SIMPLIFY BYTE
Array[10]
SRM_ATTEN_STATE BYTE
1
SRM_ALL_SIMPLIFY
BYTE = 0
SRM_TIME_WANTED INT
-1
SRM_EARTH_POSITION
FLOAT = 180.000
SRM_OFFAX_POSITION
FLOAT = 0.250000
SRM_INTEGRATOR FLOAT
0.000000
SRM_USE_SAVED BYTE
Array[10]
SRM_MAKE_SAVED BYTE
Array[10]
SRM_PRINT_SUBRMS
BYTE = 0
SRM_VERBOSE
BYTE 0
SRM_DISPLAY_RESULT
BYTE = 0
SRM_PHA_ON_ROW BYTE
1
SRM_DIAGONAL BYTE
0
The general procedure is to set up the parameters desired for retrieving the spectra, and then executing a getdata() call to actually get the spectrum. To differentiate from set parameters appropriate for the imaging tasks, the parameters set for spectroscopy have a prefix of sp_ which must be used. Some of the parameters used for imaging are also used for spectroscopy, and there is no real need to differentiate between the tow uses, so those parameters do not have the prefix sp_. Finally, there are a number of parameters specifically for the spectral response matrix which carry the prefix srm_. One of the first things the user will have to set is the energy binning to be used to construct the spectrum. The user can use one of 11 pre-defined binning schemes, or set their own. The pre-defined binning schemes are described in the window which pops up when the hsi_spectroscopy_list command is executed. Again, this can be done by simply typing:
IDL> hsi_spectroscopy_list
The first topic described is the default energy binning schemes. If the user wishes to use one of these, say scheme number 6, then this is selected by the set call and setting the sp_energy_binning keyword. For example, simply type:
IDL> sp -> set, sp_energy_binning=6
If none of the default energy binning schemes is suitable, the user is free to define their own energy binning. This can be done by passing energy boundaries in one of two ways. The user can pass in a one dimensional array, in which case the software interprets each element as an energy bin edge, making the bins continuous. So, if a 10 element vector, ee, is passed in for example, the software assumes there are 9 energy bins: ee(0) - ee(1); ee(1) - ee(2); ee(2) - ee(3); and so on to the last bin which is ee(8) - ee(9). In this case, you could do the following:
IDL> ee = 10.^(findgen(50)/51.*2.5)
IDL> sp -> set, sp_energy_binning = ee
which would define 49 logrithmically space energy bins from 1 to 252.325 keV. The other means is to pass a 2 dimensional array which specifies the begining and ending energies of each bin to consider. For example, you could type the following:
IDL> ee = fltarr(2,50)
IDL> ee(0,*) = 10.^(findgen(50)/51.*2.5)
IDL> ee(1,*) = 10.^((findgen(50)+1)/51.*2.5)
IDL> sp -> set, sp_energy_binning = ee
which defines 50 logrithmically spaced energy bins from 1 to 282.475 keV which are contiguous in this example. Once the energy binnin is defined an the spectra are retrieved, the user can always use the software to return the actual energy bins used if there is any doubt on his or her part.
One the energy binning has been set, the time intervals must be chosen. It is possible to process all of the data in the file you are working with, making spectra in bins of the users choice, or to specify a particular time range to work with. For example, if you have read in the flare list as decribed in the imaging section above, you may want to make a spectrum of the 4 seconds around the peak of the first flare. According to the flarelist the peak of this flare occurs at 6.83769780d+008. So, you could choose the interval starting 2 seconds before this and extending to 2 seconds after this by typing:
IDL> sp -> set, sp_time_interval = [6.83769778d+008, 6.83769782d+008]
Note that you must use double precision when specifying times in this fashion. Times can also be specified in any of the normal anytime formats. When specifying a single interval in this fashion, the return will of course be a single spectrum. It may be more useful for certain studies to process an entire file to create spectra on a regular time grid. For example, to process the entire data file using 10 second time bins, you would type:
IDL> sp -> set, sp_time_interval=10
Processing the test data file with this time interval will produce 329 spectra, 1 every 10 seconds, with each spectrum integrating the accumulted photons over a 10 second time interval. The starting time interval in such a case defaults to the begining of the data file being processed.
The final thing the user must decide is which detectors to use in constructing the spectrum. This is controlled through the same det_index_mask keyword used in the imaging algorithms. Currently, spectroscopy is supported only for the fundamental modulation in each collimator (i.e. no support for harmonics). Therefore, at most you want to set the first 9 elements of det_index_mask (a bytarr) to 1. Therefore, it is desireable to execute:
IDL> sp -> set, det_index_mask=bytarr(9)+1b
At this point, the user is ready to get the spectrum. This is done through the getdata() call:
IDL> spec = sp -> getdata()
So, executing the following sequence of commands once SSWIDL is started:
IDL> sp = hsi_spectrum(filename='D:\hessi\test_data\hsi_20000901_000000_000.fits')
IDL> sp -> set, sp_energy_binning=6
IDL> sp -> set, sp_time_interval=10
IDL> sp -> set, det_index_mask=bytarr(9)+1b
IDL> spec = sp -> getdata()
IDL> help, spec
SPEC
FLOAT = Array[199, 329]
produces an array of 329 spectra, 1 every 10 seconds, defined on 199 energy bins. In this case, the counts in all of the collimators has been summed for the final result. To get the spectra in each collimator separately, the sum_flag should be turned off in the call to getdata() or using the set command (and in fact all of the keywords set in the examples above can be set either in calls using the set command or in the call to getdata()). So, to get the spectrum in all the collimators separately, you could type:
IDL> spec = sp -> getdata(sum_flag=0)
IDL> help, spec
SPEC
FLOAT = Array[199, 329, 9]
The 3rd dimension of the array now references each collimator independently. Just like the imaging software, the spectroscopy software has memory as to the keywords that have been set. Therefore, additional calls to getdata() will not do anything unless a keyword is changed or a new keyword is set. This also means that once the sum_flag is turned off using the call above, it will remain off if a new time binning or energy binning is used. To turn it back on, a call with sum_flag = 1 would have to be executed.
The spectrum returned in the variable spec in the example above is simply the number of counts detected in each energy bin in the specified time interval. To convert this to flux units (photons s-1 cm-2 keV-1), the sp_data_unit keyword should be set to 'Flux' (or setting it to 'f' or 'F' is also satisfactory):
IDL> spec = sp -> getdata(sp_data_unit='Flux')
or
IDL> sp -> set, sp_data_unit='Flux'
IDL> spec = sp -> getdata()
To set the return back to integrated counts, simply type:
IDL> spec = sp -> getdata(sp_data_unit='Counts')
Setting the keyword to 'C' or 'c' is also accepted. The final return data format is in rate units, which gives the number of photons s-1 in each energy bin amd time interval. This can be specified with:
IDL> spec = sp -> getdata(sp_data_unit='Rate')
or with a 'R' or 'r'.
At this point, however, the spectrum returned in spec has not been corrected for grid transmission, attenuator state, detector response, etc. To make these calibrations, the sp_semi_calibrated keyword must be set. This can be done with the call:
IDL> spec = sp -> getdata(/sp_semi_calibrated)
or
IDL> sp -> set, sp_semi_calibrated=1
IDL> spec = sp -> getdata()
A particularly useful way to return the data is achieved by setting the keyword sp_data_structure. This will return all the spectra in a data structure which includes the time intervals for each spectrum and uncertainties for each data point. To get the data in this fashion, execute the following:
IDL> spec = sp -> getdata(/sp_data_structure)
IDL> help, spec
SPEC
STRUCT = -><Anonymous> Array[329]
IDL> help, /structure, spec
** Structure <c71df0>, 4 tags, length=14384,
refs=1:
UT
DOUBLE Array[2]
FLUX
FLOAT Array[199,
9]
EFLUX
FLOAT Array[199,
9]
LTIME
FLOAT Array[9]
So, by extracting the energy axis using getdata() as described below, you can plot the spectrum of the 49th time interval (index 48 - the strongest spectrum in energy channel 30) as recorded by the 1st subcollimator (index 0) using the follwing commands to plot and put error bars on the spectrum:
IDL> plot_oo,e,spec(48).flux(*,0),ps=3,yr=[1.e-4,1.e+4]
IDL> errplot,e,spec(48).flux(*,0)-spec(48).eflux(*,0),spec(48).flux(*,0)+spec(48).eflux(*,0),width=0.0001
As mentioned above, you can also use calls to getdata() to return the energy scale of the output spectra by setting the keyword /xaxis in the call. For example:
IDL> e = sp -> getdata(/xaxis)
IDL> help, e
E
FLOAT = Array[199]
The resulting energy scale, e, is the linear midpoint energy of each of the 199 energy bins in the spectra. Many users may prefer to have the returned energy scale as the logrithmic midpoint of each bin, which is of course the geometric mean of the bin boundaries. This can be returned in the getdata() call by setting the /gmean flag:
IDL> e = sp -> getdata(/xaxis, /gmean)
Other users may prefer to have the actual energy bin edges returned instead of some average value for each bin. This is accomplished by setting the /edges_2 flag in the call:
IDL> e = sp -> getdata(/xaxis, /edges_2)
IDL> help, e
E
FLOAT = Array[2, 199]
In this case, use of the /gmean flag would have no effect since the actual bin edges are returned.
If you do not use the sp_data_structure keyword, you can still get the time intervals of the returned spetrum from getdata() using the /yaxis and /edges_2 keywords:
IDL> ut = sp -> getdata(/yaxis, /edges_2)
IDL> help, ut
UT
DOUBLE = Array[2, 329]
With these basic object calls, the user should be able to extract spectra from RHESSI data for any desired interval. Note that even using the sp_semi_calibrated keyword, the returned spectra are still in need of more processing to determine the spectrum of any flare which may be in question. For example, no background subtraction has been performed on the data. To further process the data, the user may wish to use a more general solar X-ray package such as SPEX (Schwartz 1996, STIN, 9671448). SPEX for RHESSI spectroscopy analysis is supported in the data analysis software and is described below.
The keywords given in the examples above comprise some of the keywords a user may want to change when working with HESSI spectra; however, there are many additional keywords whose functions are not shown explicitly in the above examples. A full listing of RHESSI spectrum keywords can be found at:
http://hesperia.gsfc.nasa.gov/ssw/hessi/doc/params/hsi_params_spectrum_standard.htm
SPEX is a separate utility to produce a more accurate photon spectrum and obtain best-fit function parameters for the time intervals that have been selected. SPEX has its own independent command structure within the IDL/Solarsoft environment.
Documentation to start using SPEX is given in Brian Dennis' spectroscopy document, available at:
http://hesperia.gsfc.nasa.gov/~dennis/spectroscopy/first_steps.htm
