Single 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 a single ACIS chip. The input event file is provided, and the supporting data filenames are read from the header.

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/acisf11823_repro_evt2.fits[ccd_id=1]" chip1 bin=2
Running fluximage
Version: 13 December 2013

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 596 by 597 pixels, pixel size of 0.984 arcsec,
    and cover x=4096.5:5288.5:2,y=2972.5:4166.5:2.

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

The following files were created:

 The clipped counts image is:
     chip1_broad_thresh.img

 The clipped exposure map is:
     chip1_broad_thresh.expmap

 The exposure-corrected image is:
     chip1_broad_flux.img

Since we added a DataModel filter, to restrict the analysis to just the ACIS-I1 chip, we had to explicitly list the event file, rather than just giving the directory name (e.g. setting infile=repro/). You can check the parameter file that was used with plist fluximage.

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

1. Create An Image

First, we need to create the image which will ultimately be normalized by the 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

In this case, four chips were on for the observation. Using ds9 to display 11823.evt2, we know that the source we are interested in is at sky coordinates (called "physical coordinates" in ds9) close to x = 4142, y = 3969. We can use dmcoords to find out which chip the source is on.

unix% dmcoords 11823.evt2 
dmcoords>: sky 4142 3969
(RA,Dec):    16:17:36.252     -51:02:24.96   
(RA,Dec):       244.40105     -51.04027 deg
THETA,PHI          1.110'        281.39 deg
(Logical):        4142.00       3969.00
SKY(X,Y):         4142.00       3969.00
DETX,DETY         4123.24       3963.79
CHIP ACIS-I1        70.18        968.13
TDET              4163.87       4177.18

This tells us that the chip we want is ACIS-I1, i.e. ccd_id = 1. A description of the layout of the ACIS focal plane can be found in the caption of Figure 6.1 in the The Chandra Proposers' Observatory Guide.

We could have also used the non-interactive mode of dmcoords:

unix% punlearn dmcoords
unix% dmcoords 11823.evt2 x=4142 y=3969 option=sky
unix% pget dmcoords chip_id
1

Bin the Event File

For this example we decide to bin by 2, so that each pixel has a size of 0.984 arcseconds (2 * 0.492), and use the FOV file to spatially filter the file:

unix% dmcopy "11823.evt2[ccd_id=1,sky=region(repro/acisf11823_repro_fov1.fits[ccd_id=1])][bin sky=2]" 11823.i1.img
unix% dmstat 11823.i1.img centroid-
EVENTS_IMAGE
    min:	0 	      @:	( 4211.2839401 3000.3192571 )
    max:	3306 	      @:	( 4143.2839401 3970.3192571 )
   mean:	2.435382136 
  sigma:	8.7220957688 
    sum:	684564 
   good:	281091 
   null:	78307 
unix% get_sky_limits 11823.i1.img 
 version: 07 October 2016
Checking binning of image: 11823.i1.img
  Image has 601 x 598 pixels
  Pixel size is 2.0 by 2.000000000000001
  Lower left (0.5,0.5) corner is x,y= 4054.3, 2999.3
  Upper right (601.5,598.5) corner is x,y= 5256.3, 4195.3
  DM filter is:
    x=4054.3:5256.3:#601,y=2999.3:4195.3:#598
  mkexpmap xygrid value is:
    4054.3:5256.3:#601,2999.3:4195.3:#598

By spatially filtering the events we include knowledge of the edge of the CHIP in the data subsapce of the image; this can be seen in the output of the dmstat call, where those pixels that lie outside the chip are ignored (i.e. the null row is not zero). It also means that the "bin sky=2" directive applies only to the SKY range covered by the chip rather than the default 8192 by 8192 area (this is seen in the output of get_sky_limits above). See the chip region FAQ for more ways of finding the chip boundaries.

2. Compute Exposure Map

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 1: Energy spectrum of the central source in RCW 103

[Version: postscript, PDF]

Figure 1: 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 Histogram

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

unix% punlearn asphist
unix% asphist pcadf11823_001N001_asol1.fits 1.asphist "11823.evt2[ccd_id=1]"

You can check the parameter file that was used with plist asphist.

Calculate the Instrument Map

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 this observation, since we are only interested in ACIS-I1:

unix% pget ardlib AXAF_ACIS1_BADPIX_FILE
/data/ciao/11823/repro/acisf11823_001N002_bpix1.fits[BADPIX1]

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_001N003_msk1.fits
unix% pset mkinstmap detsubsys=ACIS-1
unix% mkinstmap 1.instmap NONE 2.3
Pixel grid specification x0:x1:#nx,y0:y1:#ny (1:1024:#1024,1:1024:#1024):
Name of fits file + extension with obs info (11823.evt2):
Detector Name (ACIS-1):
Grating for zeroth order ARF (NONE|LETG|HETG) (NONE):
NONE, or name of ACIS window mask file (acisf11823_001N003_msk1.fits):
NONE, or the name of the parameter block file ():

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.

You can check the parameter file that was used with plist mkinstmap.

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.

Calculate the Exposure Map

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.i1.img 
Running: get_sky_limits
  version: 07 October 2016
Checking binning of image: 11823.i1.img
  Image has 601 x 598 pixels
  Pixel size is 2.0 by 2.000000000000001
  Lower left (0.5,0.5) corner is x,y= 4054.3, 2999.3
  Upper right (601.5,598.5) corner is x,y= 5256.3, 4195.3
  DM filter is:
    x=4054.3:5256.3:#601,y=2999.3:4195.3:#598
  mkexpmap xygrid value is:
    4054.3:5256.3:#601,2999.3:4195.3:#598

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
4054.3:5256.3:#601,2999.3:4195.3:#598
unix% punlearn mkexpmap
unix% pset mkexpmap xygrid=$xygrid normalize=no
unix% mkexpmap 1.asphist 1.expmap 1.instmap
grid specification syntax x0:x1:#nx,x0:x1:ny (4054.3:5256.3:#601,2999.3:4195.3:#598): 
Use Average Aspect Pointing (no): 

You can check the mkexpmap parameter file that was used with plist mkexpmap. The exposure map may be displayed in ds9 (Figure 2).

Figure 2: Exposure map for the ACIS-I1 chip

[Version: full-size]

Figure 2: Exposure map for the ACIS-I1 chip

unix% ds9 1.expmap -cmap b -scale pow

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.

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.

3. 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.i1.img 11823.i1.thresh.img expfile=1.expmap cut=1.5%
unix% dmimgthresh 1.expmap 1.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.

You can check the parameter file that was used with plist dmimgthresh.

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.i1.thresh.img 1.thresh.expmap 11823.i1.norm div
warning: CONTENT has 1 different values.
warning: DETNAM has different value...Merged...

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.i1.norm (Figure 3) are [photon/cm²/s/pixel].

You can check the parameter file that was used with plist dmimgcalc.

Figure 3: Exposure-corrected image of the ACIS-I1 chip

[Version: full-size]

Figure 3: Exposure-corrected image of the ACIS-I1 chip

unix% ds9 11823.i1.norm -scale log -cmap hsv

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 4):

unix% dmimgcalc infile=11823.i1.thresh.img,1.thresh.expmap \
          infile2= outfile=11823.i1.norm2 \
          op="imgout=img1/img2"
warning: CONTENT has 1 different values.
warning: DETNAM has different value...Merged...

Figure 4: Exposure-corrected image containing NaN values

[Version: full-size]

Figure 4: Exposure-corrected image containing NaN values

unix% ds9 11823.i1.norm2 -scale log -cmap hsv -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.i1.norm\* centroid-
File=11823.i1.norm
11823.i1.norm
    min:	0 	      @:	( 4099.224052 2974.7345895 )
    max:	0.00022522023937 	      @:	( 4181.224052 3942.7345895 )
   mean:	1.3903681672e-07
  sigma:	5.7686120868e-07
    sum:	0.049305235944
   good:	354620
   null:	0
File=11823.i1.norm2
11823.i1.norm2
    min:	0 	      @:	( 4253.224052 2974.7345895 )
    max:	0.00022522023937 	      @:	( 4181.224052 3942.7345895 )
   mean:	1.7630358164e-07
  sigma:	6.4450939927e-07
    sum:	0.049305235944
   good:	279661
   null:	74959