Calculating Spectral Weights for mkinstmap

Summary of the Syntax

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.utils 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")
Created: weights.txt

The output file - here called weights.txt - can then be used as the input for the spectrumfile parameter of mkinstmap.

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, 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 10²² 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(xlog=True, ylog=True)

The plot of the spectrum is shown in Figure 1:

Figure 1: Source spectrum

Figure 1: Source spectrum

The source model evaluated on the energy grid that will be used to create the weights. The normalisation of the powerlaw model is not important here, since the calculation of the weights involves normalizing the model values by the total model flux.

Plot the weights

The plot_instmap_weights routine plots up the weights, as shown in Figure 2.

sherpa> plot_instmap_weights(xlog=True, ylog=True)

Figure 2: Weights using a bin width of 0.25 keV

Figure 2: Weights using a bin width of 0.25 keV

Since we are evaluating the weights on the same grid used to calculate the source spectrum (Figure 1), the shape of the two plots are the same.

Save the weights

The weights can be saved to an ASCII file, in a format suitable for mkinstmap using the save_instmap_weights routine. Here we save them to the file "weights.txt".

sherpa> save_instmap_weights("weights.txt")
Created: weights.txt

Since the file is in ASCII format you can view it with various UNIX utilities, or with dmlist:

unix% head -10 weights.txt
#TEXT/SIMPLE
# X WEIGHT
#
# DATE = 2021-11-09T14:11:11 / Date and time of file creation
# MODELEXPR = (xsphabs.gal * powlaw1d.pl)
# ENERG_LO = 0.5 / [keV] Minimum energy
# ENERG_HI = 8.0 / [keV] Maximum energy
#
0.625 0.0005505499378475431
0.875 0.005133913053133582
unix% tail -5 weights.txt
6.875 0.01554423141902547
7.125 0.014513042710459261
7.375 0.013586416882823383
7.625 0.012875827543709225
7.875 0.012220642621851999

The dmstat tool can be used to show that the weights column sums up to 1, as is needed for use in mkinstmap:

unix% dmstat weights.txt sigma-
X
    min:	0.625 	      @:	1 
    max:	7.875 	      @:	30 
   mean:	4.25 
    sum:	127.5 
   good:	30 
   null:	0 

WEIGHT
    min:	0.00055054993785 	      @:	1 
    max:	0.074137428209 	      @:	6 
   mean:	0.033333333333 
    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.

Comparing weights

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.

First we use the same energy range but double the bin width to 0.5 keV to create Figure 3:

sherpa> dataspace1d(0.5, 8, 0.5, id=2)
sherpa> set_source(2, get_source())
sherpa> plot_instmap_weights(xlog=True, ylog=True)
sherpa> plot_instmap_weights(2, overplot=True)
sherpa> save_instmap_weights(2, "weights2.txt")
Created: weights2.txt

Figure 3: A comparison of bin widths of 0.25 and 0.5 keV

Figure 3: A comparison of bin widths of 0.25 and 0.5 keV

The orange line, from the second set of weights, is roughly twice as high as the original plot. The actual difference depends on how much the model varies within the new bins.

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 = np.asarray(edges)
sherpa> xlo = xe[:-1]
sherpa> xhi = xe[1:]
sherpa> load_arrays(3, xlo, xhi, xlo*0, Data1DInt)
sherpa> set_source(3, get_source())
sherpa> plot_instmap_weights(xlog=True, ylog=True)
sherpa> plot_instmap_weights(2, overplot=True)
sherpa> plot_instmap_weights(3, overplot=True)
sherpa> save_instmap_weights(3, "weights3.txt")
Created: weights3.txt

Figure 4: Using a variable bin width

Figure 4: Using a variable bin width

The green line, showing the weights from the variable bin widths, is similar to the second (orange) set of weights since the bin widths are not too different.

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 a combination of "low-level" Sherpa methids and Matplotlib commands to plot the data from the file src.arf on top of the weights, creating Figure 5:

sherpa> import matplotlib.pyplot as plt
sherpa> from sherpa.astro.plot import ARFPlot
sherpa> arf = unpack_arf('src.arf')
sherpa> aplot = ARFPlot()
sherpa> splot.prepare(arf)
sherpa> plot_instmap_weights(3)
sherpa> ax2 = plt.twinx()
sherpa> aplot.plot(overplot=True, color='orange')
sherpa> ax2.set_ylabel('ARF cm$^2$') 
sherpa> ax2.set_yscale('log')
sherpa> ax2.set_ylim(30, None)
sherpa> plt.xlim(0, 9)
sherpa> plt.show()

Figure 5: Comparing a set of weights to an ARF

Figure 5: Comparing a set of weights to an ARF

By adding a second Y axis to the plot, the ARF can be shown in the same plot, even though the Y axis values for the two data sets are very different.

The plot can be used to help you determine whether the choice of bins is sensible or not by seeing how the ARF varies within a bin. Ideally, the bins 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).

How do the resulting instrument maps differ?

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 pixelgrid="1:1024:#1024,1:1024:#1024"
unix% pset mkinstmap obsfile=repro/acisf01838_repro_evt2.fits
unix% pset mkinstmap detsubsys=ACIS-7
unix% pset mkinstmap maskfile=repro/acisf01838_000N003_msk1.fits

Bad-pixel files

In this example, since we are only interested in how the effective area varies with the weighting, we have not used acis_set_ardlib to ensure that the observation-specific bad-pixel file has been used.

The maps that were created are given below:

Mono-energetic: imap.e1.7.fits

unix% mkinstmap imap.e1.7.fits NONE 1.7 mode=h

0.25 keV bin width: imap.weights.fits

unix% mkinstmap imap.weights.fits weights.txt 1 mode=h

0.5 keV bin width: imap.weights2.fits

unix% mkinstmap imap.weights2.fits weights2.txt 1 mode=h

Variable bin width: imap.weights3.fits

unix% mkinstmap imap.weights3.fits weights3.txt 1 mode=h

On a reasonably modern Linux machine the run times for mkinstmap were approximately 3 s, 9 s, 6 s, and 5 s for these four maps. 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% dmstat imap\*fits cen- med+ sig-
File=imap.e1.7.fits
IMAP[cm**2]
    min:	0 	      @:	( 1 1 )
    max:	634.42694092 	      @:	( 212 432 )
   mean:	593.81272372 
 median:	602.91555786
    sum:	622657770.59 
   good:	1048576 
   null:	0 
File=imap.weights.fits
IMAP[cm**2]
    min:	0 	      @:	( 1 1 )
    max:	415.99172974 	      @:	( 212 425 )
   mean:	382.92987149 
 median:	389.11810303
    sum:	401531072.93 
   good:	1048576 
   null:	0 
File=imap.weights2.fits
IMAP[cm**2]
    min:	0 	      @:	( 1 1 )
    max:	407.64544678 	      @:	( 212 425 )
   mean:	375.16585608 
 median:	381.27789307
    sum:	393389912.71 
   good:	1048576 
   null:	0 
File=imap.weights3.fits
IMAP[cm**2]
    min:	0 	      @:	( 1 1 )
    max:	407.20336914 	      @:	( 212 425 )
   mean:	374.24775966 
 median:	380.35220337
    sum:	392427218.84 
   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
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.60567945242 	      @:	( 594 624 )
    max:	0.65803092718 	      @:	( 220 2 )
   mean:	0.64471042245 
 median:	0.64539122581 
    sum:	666638.95804 
   good:	1034013 
   null:	14563 

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 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% sherpa -n
sherpa> import matplotlib.pyplot as plt
sherpa> img = read_file('ratio1.fits').get_image().values
sherpa> fig = plt.figure(figsize=(8, 8))
sherpa> plt.imshow(img, origin='lower')
sherpa> plt.axis('off')
sherpa> plt.subplots_adjust(left=0, right=1, top=1, bottom=0)
sherpa> cax = plt.axes([0.1, 0.875, 0.6, 0.05])
sherpa> cbar = plt.colorbar(cax=cax, orientation='horizontal')
sherpa> plt.show()

Figure 6: Ratio of a weighted instrument map to a mono-energetic one

[Version: full-size]

Figure 6: Ratio of a weighted instrument map to a mono-energetic one

The ratio image shows that whilst there is a lot of small-scale differences - the ACIS-S instrument maps are calibrated per column and there are horizontal and vertical gradients - the overall image is flat within a 5% margin. The large-scale variations are more relevant since the small-scale structure tends to get washed out when creating an exposure map due to the Chandra dither pattern (see Figure 9).

In general the per-column variation is due to the variation in the detector response (the Quantum Efficiency), whilst the large-scale differences are due to the off-axis dependence of the effective area.

Comparing the weighted instrument maps show that the absolute variation between them is much smaller than to the mono-energetic case (Figure 7 and Figure 8):

unix% dmimgcalc imap.weights2.fits,imap.weights.fits out=ratio2.fits op=imgout=img1/img2 mode=h
unix% dmimgcalc imap.weights3.fits,imap.weights.fits out=ratio3.fits op=imgout=img1/img2 mode=h

Figure 7: Ratio of 0.5 keV to 0.25 keV bin-width instrument maps

[Version: full-size]

Figure 7: Ratio of 0.5 keV to 0.25 keV bin-width instrument maps

The 0.5 keV bin width instrument map has slightly smaller pixel values than the 0.25 keV case, but the difference is small (less than five per cent).

Figure 8: Ratio of variable to 0.25 keV bin-width instrument maps

[Version: full-size]

Figure 8: Ratio of variable to 0.25 keV bin-width instrument maps

As with Figure 7, the variable bin width map is generally slightly lower than the 0.25 keV version, at the less than five per cent level. The different energy-binning scheme has resulted in a change to the the large-scale pattern, reflecting both the energy and positional dependence of the effective area.

How do the resulting exposure maps differ?

Figure 9 shows the ratio of the exposure maps created using imap.e1.7.fits and imap.weights.fits:

unix% punlearn mkexpmap
unix% pset mkexpmap asphistf=7.asphist
unix% pset mkexpmap xygrid="3643.1:4723.1:#1080,3853.4:4933.4:#1080"
unix% pset mkexpmap 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

Figure 9: Ratio of 0.25 keV bin-width to mono-chromatic exposure maps

[Version: full-size]

Figure 9: Ratio of 0.25 keV bin-width to mono-chromatic exposure maps

The difference in the exposure maps is similar to the corresponding instrument map comparison shown in Figure 6.

An alternative view is the histogram of the pixel values in the various ratio images, as shown in Figure 10:

Figure 10: Histograms of the ratio images

[Version: full-size]

Figure 10: Histograms of the ratio images

The left plot compares the pixel distribution of the instrument- and exposure-map ratio images comparing the weighted result (bin width of 0.25 keV) to the mono-energetic case. The two distributions are very similar.

The right-hand plot shows how the different binning schemes affect the instrument map: in this particular case the choice of binning scheme does not seem to make a large difference on the overall result, but this is a scientific call to be made for each observation.

This figure was created using the following set of commands, making use of routines from Crates, Matplotlib, and ciao_contrib.utils.

To start with, we read the image data using the Crates routine read_file, and access the pixel-data directly. 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.

import numpy as np
import matplotlib.pyplot as plt

from pycrates import read_file
from ciao_contrib.utils import simple_hist

i1 = read_file('ratio1.fits').get_image().values
i2 = read_file('ratio2.fits').get_image().values
i3 = read_file('ratio3.fits').get_image().values
ie = read_file('emap.ratio.fits').get_image().values

ix1 = np.where(i1 > 0)
ix2 = np.where(i2 > 0)
ix3 = np.where(i3 > 0)
ixe = np.where(ie > 0)

h1 = simple_hist(i1[ix1], nbins=20)
h2 = simple_hist(i2[ix2], nbins=20)
h3 = simple_hist(i3[ix3], nbins=20)
he = simple_hist(ie[ixe], nbins=20)
plt.show()

The data can then be plotted (using the mid-point of each bin, and normalized by the total number of points), in two horizontally-aligned plots:

def hplot(ax, h, lbl):
    ax.plot((h.xlow + h.xhigh) / 2, h.y / h.y.sum(), label=lbl)

fig, (ax1, ax2) = plt.subplots(1, 2, sharey='row', figsize=(6, 4))
plt.subplots_adjust(wspace=0.05)

hplot(ax1, h1, 'ratio1.fits')
hplot(ax1, he, 'emap.ratio.fits')
ax1.legend(loc='upper left')
ax1.set_ylim(-0.005, 0.17)
ax1.set_ylabel('Fraction of pixels')

hplot(ax2, h2, 'ratio2.fits')
hplot(ax2, h3, 'ratio3.fits')
ax2.legend(loc='upper left')

plt.text(0.5, 0.02, 'Pixel Value', horizontalalignment='center',
         transform=fig.transFigure)

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.