Multiple Chip ACIS Exposure Map and Exposure-corrected Image

How fluximage works

When running fluximage, you are only required to provide an input event file. The script will read the related data product filenames - bad pixels, aspect solution, mask (ACIS), and dead time correction (HRC) - and look for them in the working directory. If the files are in a different location or you wish to be explicit in what files are used, all of the input filenames may be set in the parameter file.

The fluximage script runs the following tools:

• dmcopy: to create an event image at with the specified binning factor
• hrc_bkgrnd_lookup, reproject_events, and dmimgcalc: to subtract the particle background (HRC-I only)
• asphist: to build the aspect histogram(s)
• mkinstmap: to calculate the instrument map(s) for the center of each energy band
• mkexpmap: to calculate the exposure map(s) in each energy band
• dmimgcalc: to combine the exposure maps (multi-chip/plate case only)
• dmimgthresh: to make a "threshold cut" before dividing the image by the exposure map, removing the hot pixels at the edges (optional)
• dmimgcalc: to normalize the image by the exposure map

For the multi-chip ACIS or multi-plate HRC-S cases, the tools are run once per chip or plate and are then combined into an image of the full detector.

By default, the intermediate per-chip data products are removed after the script has completed running. To save these (potentially numerous) files, set the cleanup parameter to no.

Run the script

In this example, the script is run for all chips in an ACIS-I observation. The input event file is provided, and the supporting data filenames are read from the header. We use the default bin size (which is 8 for ACIS and 32 for HRC data).

The energy range is restricted from 0.5 keV to 7 keV, and a center-band energy of 2.3 keV is used. This corresponds to the Chandra Source Catalog broad band, so we can set the bands parameter as bands=0.5:7.0:2.3 or bands=broad (which is the default value).

unix% punlearn fluximage
unix% fluximage repro/ flux/
Running fluximage
Version: 13 December 2013

Found repro/acisf11823_repro_evt2.fits
Using event file repro/acisf11823_repro_evt2.fits
Using CSC ACIS broad science energy band.
Aspect solution repro/pcadf391758163N002_asol1.fits found.
Bad-pixel file repro/acisf11823_repro_bpix1.fits found.
Mask file repro/acisf11823_001N002_msk1.fits found.

The output images will have 297 by 297 pixels, pixel size of 3.936 arcsec,
    and cover x=2912.5:5288.5:8,y=2824.5:5200.5:8.

Running tasks in parallel with 4 processors.
Creating aspect histograms for obsid 11823
Creating 4 instrument maps for obsid 11823
Creating 4 exposure maps for obsid 11823
Combining 4 exposure maps for obsid 11823
Thresholding data for obsid 11823
Exposure-correcting image for obsid 11823

The following files were created:

 The clipped counts image is:
     flux/broad_thresh.img

 The clipped exposure map is:
     flux/broad_thresh.expmap

 The exposure-corrected image is:
     flux/broad_flux.img

Note that the script looked in the repro/ directory for any file matching the pattern "*evt*", finding repro/acisf11823_repro_evt2.fits. You can check the parameter file that was used with plist fluximage.

WARNING about creating a large image

If the chosen binning factor is small enough (e.g. binsize=1 for an HRC-I observation) then the following warning will be seen:

# DMCOPY (CIAO): WARNING: Creating large image: 840 MB. Current max set at 500 MB.
Increase maximum using [opt mem=n] or increase blocking to reduce size.

The output files are completely valid; the warning is just there to let you know that the files that are being created are large. At present there is no way to hide this message.

Creating a three-color image

Here we use fluximage to calculate the exposure-corrected images in three bands; in this case the soft, medium, and hard bands defined in the Chandra Source Catalog (for which the shortcut bands=csc can be used):

unix% fluximage repro/ flux/ bands=csc
Running fluximage
Version: 13 December 2013

Found repro/acisf11823_repro_evt2.fits
Using event file repro/acisf11823_repro_evt2.fits
Using CSC ACIS soft science energy band.
Using CSC ACIS medium science energy band.
Using CSC ACIS hard science energy band.
Aspect solution repro/pcadf391758163N002_asol1.fits found.
Bad-pixel file repro/acisf11823_repro_bpix1.fits found.
Mask file repro/acisf11823_001N002_msk1.fits found.

The output images will have 297 by 297 pixels, pixel size of 3.936 arcsec,
    and cover x=2912.5:5288.5:8,y=2824.5:5200.5:8.

Running tasks in parallel with 4 processors.
Creating aspect histograms for obsid 11823
Creating 12 instrument maps for obsid 11823
Creating 12 exposure maps for obsid 11823
Combining 4 exposure maps for 3 bands (obsid 11823)
Thresholding data for obsid 11823
Exposure-correcting 3 images for obsid

The following files were created:

 The clipped counts images are:
     flux/soft_thresh.img
     flux/medium_thresh.img
     flux/hard_thresh.img

 The clipped exposure maps are:
     flux/soft_thresh.expmap
     flux/medium_thresh.expmap
     flux/hard_thresh.expmap

 The exposure-corrected images are:
     flux/soft_flux.img
     flux/medium_flux.img
     flux/hard_flux.img

and display the image in ds9 (Figure 1).

Figure 1: Three-color image (exposure corrected) of RCW 103

[Version: full-size]

Figure 1: Three-color image (exposure corrected) of RCW 103

unix% ds9 -rgb -red soft_flux.img -green medium_flux.img -blue hard_flux.img \
         -rgb lock scale yes -scale mode 99.5 -zoom 2

The image can be compared to that obtained using two observations of RCW 103.

The choice of the CSC soft, medium, and hard bands for this image was made for didactic reasons. Other choices may make more sense for your science goals.

Now proceed to the Calculate the Source Flux section, using the exposure-corrected image (here, flux/broad_flux.img).

1. Syntax note: foreach

This thread uses "foreach" loops to run the same CIAO tool for multiple chips. The syntax is correct for the csh or tcsh shell. If you are using another shell, e.g. bash, change the loop syntax accordingly.

2. Create An Image

First, we need to create the image which will ultimately be normalized by the exposure map. Here we decided to block the image by a factor of 8:

unix% dmcopy "11823.evt2[sky=region(acisf11823_001N002_fov1.fits)][bin sky=8]" 11823.i0123.img
unix% get_sky_limits 11823.i0123.img 
Running: get_sky_limits
  version: 12 September 2012
Checking binning of image: 11823.i0123.img
  Image has 297 x 297 pixels
  Pixel size is 8.0 by 8.0
  Lower left (0.5,0.5) corner is x,y= 2916.7, 2824.8
  Upper right (297.5,297.5) corner is x,y= 5292.7, 5200.8
  DM filter is:
    x=2916.7:5292.7:#297,y=2824.8:5200.8:#297
  mkexpmap xygrid value is:
    2916.7:5292.7:#297,2824.8:5200.8:#297

This uses the information in the FOV file to limit the sky area of the image to that covered by the observation. See the chip region FAQ for more ways of finding the chip boundaries.

3. Compute Exposure Map

Check Which Chips Are On

The list of chips used in the observation is stored in the DETNAM keyword of the event file:

unix% dmkeypar 11823.evt2 detnam echo+
ACIS-0123

A description of the layout of the ACIS focal plane can be found in The Chandra Proposers' Observatory Guide.

What is the spectrum of the source?

We selected a region around the central source, using ds9, and saved it as obj.reg:

unix% cat obj.reg
# Region file format: CIAO version 1.0
circle(4142.9,3969.4,8)

We can use this file to extract a spectrum of the object in energy space and find the peak energy.

First, we use the CIAO tool dmextract to create a histogram of count-rate as a function of energy. Since we are not binning on pi or pha, we set opt=generic, and we use a bin size of 50 eV to improve the signal to noise:

unix% punlearn dmextract
unix% pset dmextract opt=generic
unix% dmextract "11823.evt2[sky=region(obj.reg)][bin energy=500:7000:50]" 11823.energy.fits

The dmstat tool is used to find the maximum count from the histogram, followed by dmlist to locate the corresponding energy:

unix% dmstat "11823.energy.fits[cols counts]" sigma-
COUNTS[count]
    min:	1 	      @:	3 
    max:	268 	      @:	25 
   mean:	69.423076923 
    sum:	9025 
   good:	130 
   null:	0 

unix% dmlist "11823.energy.fits[counts>250][cols energy,counts]" data,clean
#  ENERGY               COUNTS
               1675.0        253
               1725.0        268
               1775.0        258

For this dataset, the peak of the measured spectrum is ~1.7 keV (which is expected since this is close to the peak of the ACIS effective area). Using the peak value would mean that we would be under-estimating the flux if the energy band is too broad; see the discussion of band selection in the Chandra Source Catalog for more information. So we will use the 2.3 keV used by the CSC, but note that this is something that depends on the spectrum of the source (or sources) being analysed.

Figure 2: Energy spectrum of the central source in RCW 103

[Version: postscript, PDF]

Figure 2: Energy spectrum of the central source in RCW 103

The two vertical lines indicate the 1.7 keV peak and the 2.3 keV value defined by the Chandra Source Catalog for its "broad" band.

Compute the Aspect Histograms

With the aspect solution file we can create a binned histogram for each chip, detailing the aspect history of the observation.

The binning used by the asphist tool is controlled by its res_xy parameter. The default value of 0.5 arcsec was chosen to match the high-resolution capabilities of Chandra. Since we are creating images - and therefore exposure maps - at a binning significantly larger than the natural pixel scale of ACIS, we can increase the res_xy parameter and so reduce the time it takes to run mkexpmap. To ensure that aliasing is not significant, we chose the resolution to be half the pixel size; that is 1.9 arcseconds (rounding down from 4 * 0.492 = 1.968).

In some cases there will be more than one aspect solution file (pcadXXX_asol1.fits) for an observation. All the files must be input, either as a list or as a stack. If you used chandra_repro to re-process the data then it has created a stack file for you, called acisf<obsid>_asol1.lis, which we use in this case (although as we only have a single aspect solution we could also have just used it directly):

unix% cat acisf11823_asol1.lis 
/data/ciao/11823/primary/pcadf391758163N002_asol1.fits
unix% punlearn asphist
unix% pset asphist infile=@acisf11823_asol1.lis
unix% pset asphist res_xy=1.9 mode=h
unix% foreach d ( 0 1 2 3 )
foreach? asphist outfile=${d}.asphist evtfile="11823.evt2[ccd_id=${d}]"
foreach? end

Calculate the Instrument Maps

Since the mirror effective area is used to create the instrument map, and that area is energy dependent, it is necessary to decide at what energy to perform the calculation (or whether to use a spectrum as weights). In this example we are going to assume a monoenergetic distribution of source photons of 2.3 keV (monoenergy parameter). The Calculating Spectral Weights for mkinstmap thread shows how to create a weighted instrument map using mkinstmap.

Note that it is not necessary for the instrument map to be congruent with the exposure map; the instrument map should describe the chip with full resolution.

At this point make sure that you have set up ardlib to use the bad pixel file for your observation. For example, checking the setting for ACIS-I3 we get:

unix% pget ardlib AXAF_ACIS3_BADPIX_FILE
/data/ciao/11823/repro/acisf11823_001N002_bpix1.fits[BADPIX3]

where we are using the bad-pixel file created by chandra_repro for this observation.

unix% punlearn mkinstmap
unix% pset mkinstmap pixelgrid="1:1024:#1024,1:1024:#1024"
unix% pset mkinstmap obsfile=11823.evt2 
unix% pset mkinstmap maskfile=acisf11823_001N002_msk1.fits
unix% pset mkinstmap monoenergy=2.3 mode=h

The monoenergy parameter value may be different for your science objectives.

Including the maskfile parameter is particularly important if you are interested in having an accurate exposure map at the very edge of a CCD, subarray or window. The pixelgrid parameter should not be changed for the case of a subarray or window; the mask file will account for the detector range being different. For more information, see the dictionary entry on mask files.

The pbkfile parameter has been deprecated and should be left empty; more details can be found on the Watchout page. The obsfile parameter should use the event file rather than the aspect histogram, as used in previous versions of CIAO.

Now run the tool once for each chip:

unix% foreach d ( 0 1 2 3 )
foreach? mkinstmap detsubsys=ACIS-${d} outfile=${d}.instmap
foreach? end

Calculate the Exposure Maps

Now we use mkexpmap and the aspect information stored in the histogram to project the instrument map onto the sky. We need to set the xygrid parameter to produce an exposure map that is the same size as the image created from the event list. The get_sky_limits script can be used to easily calculate this information from the existing image:

unix% get_sky_limits 11823.i0123.img 
Running: get_sky_limits
  version: 07 October 2016
Checking binning of image: 11823.i0123.img
  Image has 298 x 298 pixels
  Pixel size is 8.0 by 8.0
  Lower left (0.5,0.5) corner is x,y= 2872.7, 2850.5
  Upper right (298.5,298.5) corner is x,y= 5256.7, 5234.5
  DM filter is:
    x=2872.7:5256.7:#298,y=2850.5:5234.5:#298
  mkexpmap xygrid value is:
    2872.7:5256.7:#298,2850.5:5234.5:#298

You can then set the xygrid parameter using the information provided by the script, either manually or via

unix% pset mkexpmap xygrid=")get_sky_limits.xygrid"

(if the latter, do not run get_sky_limits again until after running mkexmap).

Note: If you are computing a low-resolution exposure map and speed is more important than accuracy, set useavgaspect=yes. In doing so, only the average aspect pointing will be used to derive the exposure map; otherwise all points in the aspect histogram will be used. The time required to compute the exposure map is proportional to the number of bins in the aspect histogram; if the aspect histogram contains 100 bins, then the use of this option reduces the run time by a factor of 100, approximately (you may also want to set verbose to 2, since this causes mkexpmap to output percentage-completed information). Using the full aspect solution will help accurately account for chip edges, bad pixels, etc.

unix% set xygrid = `pget get_sky_limits xygrid`
unix% echo $xygrid
2872.7:5256.7:#298,2850.5:5234.5:#298
unix% punlearn mkexpmap
unix% pset mkexpmap xygrid=$xygrid normalize=no mode=h
unix% foreach d ( 0 1 2 3 )
foreach? mkexpmap instmapfile=${d}.instmap \
                  outfile=${d}.bin8.expmap \
                  asphistfile=${d}.asphist
foreach? end

Since we set the normalize parameter = no, the exposure map has units of [cm²*s*counts/photon]. This allows us to simply divide the image by the exposure map to derive an image in units of flux [photons/cm²/s/pixel]. If the setting had been left as yes (the default), the units of the exposure map would be [cm²*counts/photon]. The units can be added to the exposure map using dmhedit

unix% dmhedit expmap.fits file= op=add key=BUNIT value="cm**2 sec"

Please see the help file for mkexpmap for more details on this.

Combine the Exposure Maps

Since the individual exposure maps were all created with the same xygrid parameter then we can combine them with a single dmimgcalc call.

unix% ls -1 ?.bin8.expmap
0.bin8.expmap
1.bin8.expmap
2.bin8.expmap
3.bin8.expmap
unix% punlearn dmimgcalc
unix% dmimgcalc '?.bin8.expmap' infile2= out=bin8.expmap
arithmetic operation (): imgout=img1+img2+img3+img4
warning: DETNAM has different value...Merged...
dmkeypar bin8.expmap detnam echo+
Merged
unix% dmhedit bin8.expmap file= op=add key=DETNAM value=ACIS-0123
unix% dmkeypar bin8.expmap detnam echo+
ACIS-0123

The exposure map can be displayed in ds9 (Figure 3). We adjust the DETNAM keyword to better reflect the data (note that the DETNAM for each individual exposure map is ACIS-0 to ACIS-3 in this example).

Figure 3: Exposure map for the ACIS-I array

[Version: full-size]

Figure 3: Exposure map for the ACIS-I array

unix% ds9 bin8.expmap -cmap b -scale pow -zoom 2

The stripes are due to bad columns in the instrument map and the "power-law" scaling has been chosen to highlight the presence of bad pixels.

4. Normalize the Image by the Exposure Map

The strongly variable exposure near the edge of a dithered field may produce "hot" pixels when divided into an image. While technically proper, these hot pixels can be an eyesore, drawing attention to a noisy, uninteresting portion of the image. The dmimgthresh tool is used to make a "threshold cut" before dividing the image by the exposure map, thus removing the hot pixels:

unix% punlearn dmimgthresh
unix% dmimgthresh 11823.i0123.img 11823.i0123.thresh.img expfile=bin8.expmap cut=1.5%
unix% dmimgthresh bin8.expmap bin8.thresh.expmap cut=1.5%

We also threshold the exposure map so that we can use the exposure map to determine whether a pixel is 0 because there were no counts or because it was removed by the threshold process.

Here we set our threshold at 1.5% of the maximum value of the exposure map. All image pixels with values of exposure less than this value will be set to 0.0 in the output file. You may want to adjust these values for your own observation.

The exposure map is in units of [cm²*s*counts/photon] since it was created by projecting the instrument map (in [cm²*counts/photon]) onto the tangent plane of the observation. To create an image in units of [photon/cm²/s/pixel], we simply need to divide by the exposure map. This is done by the tool dmimgcalc.

unix% punlearn dmimgcalc
unix% dmimgcalc 11823.i0123.thresh.img bin8.thresh.expmap 11823.i0123.norm div
warning: CONTENT has 1 different values.

The messages are related to how the tool merges the header information in the input files. The merging_rules ahelp file explains the rules and how they affect the output file header.

The units of 11823.i0123.norm.div (Figure 4) are photon/cm²/s/pixel.

Figure 4: Exposure-corrected image of the ACIS-I observation

[Version: full-size]

Figure 4: Exposure-corrected image of the ACIS-I observation

unix% ds9 11823.i0123.norm -scale log -cmap hsv -zoom 2

It is also possible to use dmimgcalc to create an exposure-corrected image where those pixels with no exposure are set to NaN (see Figure 5):

unix% dmimgcalc infile=11823.i0123.thresh.img,bin8.thresh.expmap out=11823.i0123.norm2 \
          infile2= op="imgout=img1/img2"
warning: CONTENT has 1 different values.

Figure 5: Exposure-corrected image containing NaN values

[Version: full-size]

Figure 5: Exposure-corrected image containing NaN values

unix% ds9 11823.i0123.norm2 -scale log -cmap hsv -zoom 2 -prefs nancolor red

The red pixels indicate the NaN values, where the exposure map is 0.

Many CIAO tools will exclude NaN values, as shown in the dmstat output:

unix% dmstat 11823.i0123.norm\* centroid-
File=11823.i0123.norm
11823.i0123.norm
    min:	0 	      @:	( 2876.7008031 2854.5276876 )
    max:	0.00046434541582 	      @:	( 4140.7008031 3966.5276876 )
   mean:	1.3656299294e-06
  sigma:	4.1610893393e-06 
    sum:	0.12127340025 
   good:	88804 
   null:	0 
File=11823.i0123.norm2
11823.i0123.norm2
    min:	0 	      @:	( 3180.7008031 2862.5276876 )
    max:	0.00046434541582 	      @:	( 4140.7008031 3966.5276876 )
   mean:	1.7680142325e-06 
  sigma:	4.6588654757e-06 
    sum:	0.12127340025 
   good:	68593 
   null:	20211