Calculating Spectral Weights for mkinstmap
CIAO 4.2 Science Threads
Last Update: 19 Jul 2010 - the S-Lang syntax has been removed from this thread as it is not supported in CIAO 4.2 Sherpa v2. Contact the Helpdesk if you need assistance transitioning from S-Lang to Python.
When a spectral weights file is provided to mkinstmap, the tool computes a weighted instrument map instead of a monochromatic instrument map. Essentially, a weighted instrument map is a linear combination of monochromatic maps calculated for several energy-weight pairs. More details are given in the Why are weighted instrument maps important? section of this thread.
An example set of commands used to create the weight file follows; please remember to change the data range, binning, spectral model and parameters to those appropriate for your data!
sherpa> from sherpa_contrib import * sherpa> dataspace1d(0.5, 8, 0.25) sherpa> set_source(xsphabs.gal * powlaw1d.pl) sherpa> gal.nh = 1.3 sherpa> pl.gamma = 1.7 sherpa> save_instmap_weights("weights.txt")
The output file - here called weights.txt - can then be used as the input for the spectrumfile parameter of mkinstmap.
Use Sherpa to create a spectral weights file than can be used by mkinstmap to create a weighted instrument map.
Read this thread if:
you are working with an ACIS or HRC imaging observation and would like to create a spectral weights file for use in one of the exposure map threads.
- Analysis Guide: Extended Sources
- Get Started
- Why are weighted instrument maps important?
- Create a model in Sherpa
- Plot the weights
- Save the weights
- Comparing weights
- How do the resulting instrument maps differ?
- How do the resulting exposure maps differ?
- Notes and Caveats
- Figure 1: Source spectrum
- Figure 2: Weights using a bin width of 0.25 keV
- Figure 3: A comparison of bin widths of 0.25 and 0.5 keV
- Figure 4: Using a variable bin width
- Figure 5: Comparing a set of weights to an ARF
- Figure 6: Ratio of a weighted instrument map to a mono-energetic one
- Figure 7: Ratio of 0.5 keV to 0.25 keV bin-width instrument maps
- Figure 8: Ratio of variable to 0.25 keV bin-width instrument maps
- Figure 9: Ratio of 0.25 keV bin-width to mono-chromatic exposure maps
- Figure 10: Histograms of the ratio images
Sample ObsID used: 1838 (ACIS-S, G21.5-0.9)
No data file is actually needed to create the weights file. As we create instrument- and exposure- maps to see what difference adding a weights file makes, we require the same data products for these steps as the Single Chip Exposure Map thread.
This thread uses the sherpa_contrib package, which can be loaded into Sherpa using the following command:
sherpa> from sherpa_contrib import *
An error message of
ImportError: No module named sherpa_contrib
means that the CIAO contributed scripts package has not been installed or is out of date.
This script requires that the contributed scripts package is no older than August 3, 2009. This can be checked by saying
unix% cat $ASCDS_CONTRIB/VERSION.CIAO_scripts 03 Aug 2009
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/cm2/arcsec2 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.
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, we can use the Sherpa dataspace1d command to create an energy grid from 0.5 to 8 keV with a bin width of 0.25 keV:
sherpa> dataspace1d(0.5, 8, 0.25)
Now we want to define a model that represents the source emission: we do not need to worry about the overall normalization of the model, just that the relative amplitudes of any emission components is correct.
For this example we are going to use a power law, with a gamma of 1.7, which is absorbed by a line-of-sight column density of 1.3 x 1022 atoms cm-2. The models and parameter values should be changed to match the emission from the source of interest.
sherpa> set_source(xsphabs.gal * powlaw1d.pl) sherpa> gal.nh = 1.3 sherpa> pl.gamma = 1.7 sherpa> plot_source() sherpa> log_scale()
The plot of the spectrum is shown in Figure 1:
unix% head -5 weights.txt #TEXT/SIMPLE # X WEIGHT 6.2500000000000e-01 2.6220696092014e-04 8.7500000000000e-01 6.8033671439539e-03 unix% tail -5 weights.txt 6.875000000000 1.5496657006506e-02 7.125000000000 1.4520638260218e-02 7.375000000000 1.3549337557838e-02 7.625000000000 1.2839676215019e-02 7.875000000000 1.2185290954080e-02
unix% dmstat weights.txt X min: 0.625 @: 1 max: 7.875 @: 30 mean: 4.25 sigma: 2.1638603621 sum: 127.5 good: 30 null: 0 WEIGHT min: 0.00026220696092 @: 1 max: 0.074007102956 @: 6 mean: 0.033333333333 sigma: 0.020733677512 sum: 1 good: 30 null: 0
The output table gives the energy at the center of each band (in keV, X) and the fraction of the incident flux falling in that band (the spectral weighting, WEIGHT). 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.
In the following sections of the thread we show how you can change the binning used to create the weights and visualize the differences. Note that using a bin width of 0.25 keV is generally excessive, in that it over-resolves the energy grid, and so will lead to longer run times of mkinstmap than are strictly necessary.
If you wish to use a different energy range then you can call dataspace1d to create a new grid, or filter the grid by using the ignore and notice routines. In the following we create new datasets with different binning characteristics to investigate the differences these choices can make.
We now decide to use variable-width bins, manually choosing energy ranges over which we consider the effective area to be approximately flat. These values are for demonstration purposes only and should not be considered to be valid choices for your analysis!
As the dataspace1d command only creates regularly-spaced grids, we use the load_arrays command to create the grid: the xlo*0 statement is used to create a fake data array, since we ignore the data values for this calculation. The result is shown in Figure 4.
sherpa> edges = [0.5, 1, 1.5, 2, 2.25, 2.5, 2.75, 3, 3.5, 4, 6, 8] sherpa> xe = numpy.asarray(edges) sherpa> xlo = xe[:-1] sherpa> xhi = xe[1:] sherpa> load_arrays(3, xlo, xhi, xlo*0, Data1DInt) sherpa> set_source(3, gal * pl) sherpa> plot_instmap_weights(3, overplot=True) sherpa> set_histogram (S-Lang or Python help)(["line.color", "green"]) sherpa> save_instmap_weights(3, "weights3.txt")
It may be useful to plot the ARF for your source along with the weights, to check that the binning you have used is sensible. In the following we use ChIPS commands to plot the data from the file src.arf on top of the weights, creating Figure 5:
sherpa> plot_instmap_weights(3) sherpa> add_axis (S-Lang or Python help)(Y_AXIS, 1, 1, 2) sherpa> add_curve (S-Lang or Python help)("src.arf[cols energ_hi,specresp]", ["symbol.style", "none"]) sherpa> set_plot_ylabel("ARF (cm^2)") sherpa> log_scale (S-Lang or Python help)(Y_AXIS) sherpa> limits (S-Lang or Python help)(Y_AXIS, 30, AUTO) sherpa> limits (S-Lang or Python help)(X_AXIS, 0, 9)
In this section we create instrument maps for the weights created above, and compare them to a mono-energetic version. All the instrument maps were created using the following set of parameters (taken from the single chip exposure map thread):
unix% punlearn mkinstmap unix% pset mkinstmap obsfile="asphist_7.fits[asphist]" unix% pset mkinstmap maskfile=acisf01838_000N002_msk1.fits unix% pset mkinstmap pbkfile=acisf084245776N002_pbk0.fits unix% pset mkinstmap dafile=CALDB unix% pset mkinstmap detsubsys=ACIS-7 unix% pset mkinstmap pixelgrid="1:1024:#1024,1:1024:#1024"
The maps that were created are given below:
- Mono-energetic: imap.e1.7.fits
unix% pset mkinstmap spectrumfile=NONE unix% pset mkinstmap monoenergy=1.7
- 0.25 keV bin width: imap.weights.fits
unix% pset mkinstmap spectrumfile=weights.txt
- 0.5 keV bin width: imap.weights2.fits
unix% pset mkinstmap spectrumfile=weights2.txt
- Variable bin width: imap.weights3.fits
unix% pset mkinstmap spectrumfile=weights3.txt
On a reasonably modern Linux machine the run times for mkinstmap were approximately 5 s, 20 s, 10 s, and 10 s for these four maps (rounded to the nearest 5 seconds). The actual run times depend strongly on operating system and disk access, so these numbers should only be taken as rough guides.
A quick comparison using dmstat shows that the three weighted maps are similar, but significantly different to the mono-energetic case:
unix% cat istack imap.e1.7.fits imap.weights.fits imap.weights2.fits imap.weights3.fits unix% dmstat @istack cen- med+ sig- File=imap.e1.7.fits IMAP min: 0 @: ( 1 1 ) max: 628.43005371 @: ( 216 416 ) mean: 589.88732488 median: 598.66592407 sum: 618541691.57 good: 1048576 null: 0 File=imap.weights.fits IMAP min: 0 @: ( 1 1 ) max: 414.24081421 @: ( 216 416 ) mean: 382.64538286 median: 388.6978302 sum: 401232764.98 good: 1048576 null: 0 File=imap.weights2.fits IMAP min: 0 @: ( 1 1 ) max: 414.68215942 @: ( 216 416 ) mean: 383.25308209 median: 389.25798035 sum: 401869983.81 good: 1048576 null: 0 File=imap.weights3.fits IMAP[cm**2] min: 0 @: ( 1 1 ) max: 417.61373901 @: ( 216 416 ) mean: 385.95979285 median: 392.26676941 sum: 404708175.75 good: 1048576 null: 0
Comparing the median values - shown in bold above - shows that the weighted maps are roughly 65 % of the mono-energetic value. Dividing imap.weights.fits by imap.e1.7.fits produces the image shown in Figure 6.
unix% punlearn dmimgcalc mode=h unix% dmimgcalc imap.weights.fits,imap.e1.7.fits out=ratio1.fits op=imgout=img1/img2 mode=h unix% dmstat ratio1.fits cen- sig- med+ ratio1.fits min: 0.61701643467 @: ( 594 592 ) max: 0.66749340296 @: ( 225 2 ) mean: 0.64855308378 median: 0.64922356606 sum: 672037.83949 good: 1036211 null: 12365
The use of imgout=img1/img2 for the operation parameter of dmimgcalc ensures that 0-pixel values get replaced by NaN values, and so will be ignored by dmstat, ChIPS, and ds9. If we had used op=div and set both infile and infile2 then these pixels would have been set to 0 instead.
unix% chips ----------------------------------------- Welcome to ChIPS: CXC's Plotting Package ----------------------------------------- CIAO 4.2 Tuesday, July 6, 2010 chips-1> add_window(["width", 8, "height", 8, "units", "inches"]) chips-2> add_plot(0, 0, 1, 1, ["style", "open"]) chips-3> add_image("ratio1.fits", ["colormap", "rainbow", "axis.pad", 0]) chips-4> add_colorbar(0.4, 0.9, ["width", 0.05, "length", 0.6]) chips-5> set_colorbar(["*.color", "black", "*.thickness", 2, "ticklabel.size", 16, "ticklabel.fontstyle", "bold"])
unix% punlearn mkexpmap unix% pset mkexpmap xygrid="3643.1:4723.1:#1080,3853.4:4933.4:#1080" unix% pset mkexpmap asphistf=asphist_7.fits unix% pset mkexpmap useavgaspect- normalize- unix% mkexpmap instmap=imap.e1.7.fits outf=emap.e1.7.fits mode=h unix% mkexpmap instmap=imap.weights.fits outf=emap.weights.fits mode=h unix% dmimgcalc emap.weights.fits,emap.e1.7.fits outf=emap.ratio.fits op=imgout=img1/img2 mode=h
An alternative view is the histogram of the pixel values in the various ratio images, as shown in Figure 10:
First we ensure that the ciao_contrib.utils module is loaded:
sherpa> from ciao_contib.utils import *
Now we read in the data we want to plot using the Crates routines read_file (S-Lang or Python help) and get_piximgvals (S-Lang or Python help). We then filter the image values to exclude those pixels with values of 0, and use the filtered set of pixels in the call to simple_hist (S-Lang or Python help).
sherpa> r1 = read_file("ratio1.fits") sherpa> r2 = read_file("ratio2.fits") sherpa> r3 = read_file("ratio3.fits") sherpa> re = read_file("emap.ratio.fits") sherpa> i1 = get_piximgvals(r1) sherpa> i2 = get_piximgvals(r2) sherpa> i3 = get_piximgvals(r3) sherpa> ie = get_piximgvals(re) sherpa> ix1 = numpy.where(i1 > 0) sherpa> ix2 = numpy.where(i2 > 0) sherpa> ix3 = numpy.where(i3 > 0) sherpa> ixe = numpy.where(ie > 0) sherpa> h1 = simple_hist(i1[ix1], nbins=20) sherpa> h2 = simple_hist(i2[ix2], nbins=20) sherpa> h3 = simple_hist(i3[ix3], nbins=20) sherpa> he = simple_hist(ie[ixe], nbins=20)
With the data available now in the h1, h2, h3, and he variables we can now plot up the histograms (normalized to give the fraction of pixels in each bin):
sherpa> add_window(["width", 6, "height", 4, "units", "inches"]) sherpa> add_histogram(h1.xlow, h1.xhigh, h1.y * 1.0 / h1.y.sum()) sherpa> add_histogram(he.xlow, he.xhigh, he.y * 1.0 / he.y.sum()) sherpa> set_histogram(["line.style", "dot"]) sherpa> set_plot_ylabel("Fraction of pixels") sherpa> split(1,2) sherpa> add_histogram(h2.xlow, h2.xhigh, h2.y * 1.0 / h2.y.sum()) sherpa> set_histogram(["line.color", "red"]) sherpa> add_histogram(h3.xlow, h3.xhigh, h3.y * 1.0 / h3.y.sum()) sherpa> set_histogram(["line.color", "green"]) sherpa> bind_axes("plot2", "ay1", "plot1", "ay1") sherpa> set_yaxis(["ticklabel.visible", False]) sherpa> set_yaxis("by2", ["ticklabel.visible", True]) sherpa> limits(X_AXIS, 0.9945, 1.0059) sherpa> add_label(1.0025, 0.22, "ratio2.fits", ["color", "red"]) sherpa> add_label(0.999, 0.2, "ratio3.fits", ["color", "green", "halign", 0.5]) sherpa> current_plot("plot1") sherpa> limits(X_AXIS, 0.615, 0.669) sherpa> add_line(0.62, 0.22, 0.628, 0.22) sherpa> add_label(0.63, 0.22, "ratio1.fits", ["valign", 0.5]) sherpa> add_line(0.62, 0.2, 0.628, 0.2, ["style", "dot"]) sherpa> add_label(0.63, 0.2, "emap.ratio.fits", ["valign", 0.5]) sherpa> x1 = get_plot("plot1").leftmargin sherpa> x2 = 1 - get_plot("plot2").rightmargin sherpa> add_label(0.5*(x1+x2), 0.05, "Pixel value", ["size", 14, "coordsys", "FRAME_NORM", "halign", 0.5])
The bind_axes (S-Lang or Python help) call is used to make sure the y axes of the two plots show the same range. Since the two plots share a common axis, we decide to display the tick labels on the right-hand axis of the second plot. Labels and lines are used to add annotations to each plot (setting the valign or halign attributes of a label to 0.5 centers them vertically or horizontally at that coordinate). The last label is added using the frame-normalized (FRAME_NORM) coordinate system to ensure that it is centered between the plots (we use the get_plot (S-Lang or Python help) routine to get the coordinates of the plot edges).
Note: The comparisons above were all made with the on-axis chip; larger differences may be seen between different weighting schemes for ACIS chips that are further off axis, and so are subject to stronger, energy-dependent, vignetting.
The An Introduction to Exposure Maps document provides background reading on why weighted exposure maps are important, and under what conditions they produce meaningful results.
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.
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.
The CIAO 3.4 version of the code required you to:
- save the spectrum to a text file,
- and manually edit the spectrum.sl file if you wanted to change the bins used for the calculation of the weights.
Neither of these steps are needed now, since the calculation can access the spectal model directly and the choice of bins used for the weights is now determined by the values used in the dataspace1d or load_arrays call.
|03 Aug 2009||Updated for CIAO 4.1: spectrum.sl has been replaced by the sherpa_contrib.utils sherpa_utils module; see the Changes from CIAO 3.4 section; a comparison of weighted to mono-energetic instrument maps has been added.|
|09 Feb 2010||updated for CIAO 4.2: the "get_imagevals" command has been renamed to get_piximgvals (S-Lang or Python help)|
|12 Mar 2010||Added an example of the essential commands to the synopsis section at the start of the thread, for those users who just need reminding of the syntax; updated some of the figures to use ChIPS for display rather than ds9.|
|19 Jul 2010||the S-Lang syntax has been removed from this thread as it is not supported in CIAO 4.2 Sherpa v2. Contact the Helpdesk if you need assistance transitioning from S-Lang to Python.|