Calculating Spectral Weights

Get Started

Sample ObsID used: 1838 (ACIS-S, G21.5-0.9)

No data file is actually needed to complete this thread.

This thread uses the spectrum.sl script. The most recent version of spectrum.sl is v2.1 (11 July 2004):

unix% grep Version $ASCDS_CONTRIB/share/slsh/local-packages/spectrum.sl
% Version: 2.1 (11 July 2004)

Note that $ASCDS_CONTRIB/share/slsh/local-packages/ is the default path in the standard CIAO scripts installation; see the Scripts page for more information. Please check that you are using the most recent version before continuing. If you do not have the script installed or need to update to a newer version, please refer to the Scripts page.

Since you will need to edit the script text to define the energy bands, copy the script to the working directory before continuing.

Why are weighted instrument maps important?

In dividing a counts image by an exposure map, the primary aim is to obtain an image with units of surface brightness in photons/s/cm²/arcsec² or, equivalently, in flux/pixel. It is important to understand that this can yield an accurate answer only if the number of counts in each pixel is actually proportional to the integral of the incident source spectrum over the relevant energy band.

Due to the strong energy dependence of the effective area, this is not always true. In particular, the number of detected counts is proportional to the incident flux in a broad band only when the effective area is constant within the chosen band. Large systematic errors may result if one blindly applies a monochromatic exposure map to a broad band image which encompasses a large range of variation in the effective area.

Nevertheless, it is often desirable both to examine broader energy bands and to examine energy bands where the response is not particularly flat. In such cases, if the broad features of the incident spectrum are known (e.g. thermal emission from a cluster) and if the spectral variations within the image are not too large, one can reduce the systematic errors by computing an exposure map which is weighted according to a specific model for the incident spectrum. An Introduction to Exposure Maps (PS, 12pp) gives a comparison of MARX simulated data analyzed using both weighted and monochromatic exposure maps.

Create a model in Sherpa

Since the weights are largely independent of the responses, it is not necessary to use the PHA and ARF files to define a model. Instead, use the Sherpa dataspace command:

unix% sherpa

-----------------------------------------------------
Welcome to Sherpa: CXC's Modeling and Fitting Program
-----------------------------------------------------
Version: CIAO 3.4 

Type AHELP SHERPA for overview.
Type EXIT, QUIT, or BYE to leave the program.

Notes:
    Temporary files for visualization will be written to the directory: 
    /tmp
    To change this so that these files are not deleted when you exit Sherpa,
    edit $ASCDS_WORK_PATH in your 'ciao' setup script.

    Abundances set to Anders & Grevesse

sherpa> dataspace (0.1:10:0.01) histogram
sherpa> show

Optimization Method: Levenberg-Marquardt
Statistic:           Chi-Squared Gehrels

-----------------
Input data files:
-----------------

dataspace 1 (0.100000:10.000000:0.010000) histogram

The dataspace command sets a range on which a source model may be evaluated when no dataset has been input. Here we used a range of 0.1-10 keV, moving in steps of 0.01 keV.

Now we want to define a model and fit it to this spectrum; this model will be used by the script to calculate the weights. For more information on defining models, see the Sherpa Reference Manual. In this case, we use a combination of photoelectric absorption and a 1-D power law:

sherpa> paramprompt off
Model parameter prompting is off

sherpa> source =  xswabs[abs]*pow[p1]

Finally, write the source to an ASCII file:

sherpa> write model 1838.model
Write X-Axis: Bin  Y-Axis: Flux (Counts)

sherpa> junk = readascii("1838.model")
sherpa> print(junk)
_filename        =  1838.model
_path            =  /data/CIAO/spectral_weights/
_filter          =  NULL
_filetype        =  1
_header          =  NULL
_ncols           =  2
_nrows           =  990
col1             =  Float_Type[990]
col2             =  Float_Type[990]

sherpa> exit
Goodbye.

Note that one could instead use the command "write source 1838.src". The only difference is that the mid-point of the energy bin in x is reported, whereas "write model" uses the left edge. The y values are the same; that is, there is exceedingly little difference for the purpose of this thread.

Set the energy bands

There is one parameter that should be considered before the script is executed. Load spectrum.sl into a text editor and scroll down to the "Standard_Bands" section, which looks like this:

% The band positions listed here  are for illustration only -- 
% the real values should be chosen using the effective area curves.
% 
Standard_Bands["ACIS"] = init_bands ([0.5, 1.0, 1.8, 2.0, 4.0, 6.0,  8.0],
				     [1.0, 1.8, 2.0, 4.0, 6.0, 8.0, 10.0]);
Standard_Bands["HRC"]  = init_bands ([0.5, 2.0, 4.0,  8.0],
		 		     [2.0, 4.0, 8.0, 10.0]);

These numbers give the edges of the spectral bands in keV (upper row is the left edge, lower row is the corresponding right edge). These definitions are adequate for most cases.

Ideally, the bands should be chosen to define regions where the effective area curve is roughly flat (slowly varying with energy); the reasons for this are discussed in An Introduction to Exposure Maps (PS, 12pp). If you have an ARF for your observation, you can plot it in chips and use its shape to choose energy bands:

unix% chips

Welcome to ChIPS, version CIAO 3.4
Copyright (C) 1999-2003, Smithsonian Astrophysical Observatory

chips> plot "1838.arf[cols energ_hi,specresp]"
chips> curve simpleline
chips> symbol none

Which produces this plot . Having done so, we add a few extra grid points to the ACIS definition from above:

Standard_Bands["thread"] = init_bands ([0.5, 1.0, 1.5, 2.0, 2.25, 2.5, 2.75, 3.0, 3.5, 4.0, 6.0, 8.0],
				       [1.0, 1.5, 2.0, 2.25, 2.5, 2.75, 3.0, 3.5, 4.0, 6.0, 8.0, 10.0]);

and insert it in the script after Standard_Bands["HRC"], naming it thread.

Finally, move down to the end of the script and make sure that the correct detector is defined in get_standard_bands - ACIS, HRC, or one of your own naming:

   b = get_standard_bands ("thread");

Save the script before proceding. You may rename the script if you like, but make sure that you substitute the new name for "spectrum.sl" in the following final step.

Calculate the weights

Given a model for the incident spectrum, the script computes the fraction of the incident flux falling in each energy band.

S-Lang syntax is used to run the script:

unix% sherpa
...

sherpa> ()=evalfile("spectrum.sl")
sherpa> runtest("1838.model")
sherpa> quit
Goodbye.

The script creates a file called weights.txt:

unix% more weights.txt
 7.500e-01  1.577e-01
 1.400e+00  1.977e-01
 1.900e+00  3.815e-02
 3.000e+00  2.586e-01
 5.000e+00  1.538e-01
 7.000e+00  1.094e-01
 9.000e+00  8.470e-02

NOTE: if there is already a file name weights.txt in the working directory, the script will clobber it without warning. If running the script multiple times, be sure to rename any files that are to be saved.

The output table gives the energy at the center of each band (in keV) and the fraction of the incident flux falling in that band (the spectral weighting). This file may now be used to create a weighted instrument map; e.g. see the Compute an ACIS Exposure Map (Single Chip) and Build Fluxed Image thread.

Notes

1. If you wish to make several exposure maps for different bands (e.g. 1-2 keV, 2-4 keV, and 4-10 keV), rerun the script for each of these bands. Any time you use the weighting scheme, the sum of the weights should be equal to one. If one (incorrectly) makes a single set of weights with this script and then subdivides them to weight different bands, the weights of each instrument map would not sum to one. Thus, the instrument and exposure maps that are derived would have effective areas that are much too small.

2. To compute an image which gives the integrated flux over the full energy range, it may be best to first compute flux-corrected images in several narrow energy bands (where the ARF is nearly flat) and then sum those fluxed images together. Weighted exposure maps work well for an energy band where the ARF variation isn't very large, but for a full-band 0.5-10 keV image, it may not be a good idea to compute the flux by dividing the counts image by a single number. This is especially true for cases where the source spectrum varies significantly within the image; in that case, there is no general way to compute a single set of weights which will be sensible for every part of the image.