HRC-S 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 only run on the central plate, since this particular observation does not contain data for plates 1 or 3. The location of the input event file is provided, and the supporting data filenames are read from the header. Since energy is not explicitly resolved in HRC observations, the center of the energy band is chosen at the user's discretion. Here, 1.1 keV is used.

unix% punlearn fluximage
unix% fluximage repro/ casa bands=::1.1 binsize=16
Running fluximage
Version: 04 November 2021

Found repro/hrcf01038_repro_evt2.fits
Using event file repro/hrcf01038_repro_evt2.fits
Using PI=: with a monochromatic energy of 1.1 keV.
Aspect solution repro/pcadf01038_000N001_asol1.fits found.
Bad-pixel file repro/hrcf01038_repro_bpix1.fits found.
Mask file repro/hrcf01038_000N005_msk1.fits found.
Dtf file repro/hrcf01038_000N005_dtf1.fits found.

The output images will have 980 by 315 pixels, pixel size of 2.1088 arcsec,
    and cover x=24096.5:39776.5:16,y=30192.5:35232.5:16.

Running tasks in parallel with 4 processors.
Creating aspect histogram for obsid 1038
# asphist (CIAO): WARNING: skipping 35 livetime correction records (from time: 117263850.466868 to time: 117263920.166871)


Creating instrument map for obsid 1038
Creating exposure map for obsid 1038
Thresholding data for obsid 1038
Exposure-correcting image for obsid 1038

The following files were created:

 The clipped counts image is:
     casa_all_thresh.img

 The observation FOV is:
     casa_1038.fov

 The clipped exposure map is:
     casa_all_thresh.expmap

 The exposure-corrected image is:
     casa_all_flux.img

The result is the same as produced in the manual stages below (Figure 2).

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.

1. Create An Image

This thread creates an exposure map for the central plate (chip_id=2), although all three plates could be used (see the ACIS multiple chip thread for information on how to loop through all the plates).

First, we need to create the image which will ultimately be normalized by the exposure map. Here we blocked the image by a factor of 16 after filtering by the FOV file:

unix% dmcopy "hrcf01038_repro_evt2.fits[sky=region(hrcf01038_repro_fov1.fits)][bin sky=::16]" img.bin16
unix% get_sky_limits img.bin16 
Running: get_sky_limits
  version: 07 October 2016
Checking binning of image: img.bin16
  Image has 980 x 314 pixels
  Pixel size is 15.999999999999996 by 15.999999999999988
  Lower left (0.5,0.5) corner is x,y= 24105.7, 30197.1
  Upper right (980.5,314.5) corner is x,y= 39785.7, 35221.1
  DM filter is:
    x=24105.7:39785.7:#980,y=30197.1:35221.1:#314
  mkexpmap xygrid value is:
    24105.7:39785.7:#980,30197.1:35221.1:#314

2. Compute Exposure Map

Compute the Aspect Histogram

With aspect offsets file, we can create a binned histogram, detailing the aspect history of the observation. Only one aspect histogram needs to be computed because the three plates share a single GTI.

In some cases there will be more than one aspect solution file (pcadXXX_asol1.fits) for an observation. All the files must be input to the infile parameter, either as a list or as a stack. Here we use:

unix% cat hrcf01038_asol1.lis 
/data/ciao/1038/primary/pcadf01038_000N001_asol1.fits

unix% punlearn asphist
unix% asphist @hrcf01038_asol1.lis asphist.fits \
   hrcf01038_repro_evt2.fits hrcf01038_000N005_dtf1.fits
# asphist (CIAO): WARNING: skipping 35 livetime correction records (from time: 117263850.466868 to time: 117263920.166871)

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). Since energy is not explicitly resolved in HRC observations, the monoenergy parameter is determined at the discretion of the observer (the default value is 1 keV); this thread uses a value of 1.1 keV.

Note that it is not necessary for the instrument map to be congruent with the exposure map. We set the pixelgrid parameter to cover the entire detector area and bin by a factor of 8.

Be sure that the ardlib.par parameter file is updated before proceeding. chandra_repro will set it automatically or users can set it manually:

unix% plist ardlib | grep HRC-S_BADPIX
AXAF_HRC-S_BADPIX_FILE = NONE             Enter HRC-S Badpix file

unix% pset AXAF_HRC-S_BADPIX_FILE=hrcf01038_repro_bpix1.fits

Now we can create the instrument map

unix% punlearn mkinstmap
unix% pset mkinstmap pixelgrid="1:16384:#1024,1:16384:#1024"
unix% pset mkinstmap obsfile="hrcf01038_repro_evt2.fits[EVENTS]"
unix% pset mkinstmap maskfile="hrcf01038_000N005_msk1.fits[MASK2]"
unix% pset mkinstmap detsubsys=HRC-S2
unix% mkinstmap 2.instmap NONE 1.1 mode=h

Calculate the Exposure Map

Now we use mkexpmap and the aspect information stored in the histogram to project the instrument map onto the sky. The get_sky_limits script can be used to calculate the exposure map binning information from the existing image:

unix% get_sky_limits img.bin16
Running: get_sky_limits
  version: 07 October 2016
Checking binning of image: img.bin16
  Image has 980 x 314 pixels
  Pixel size is 15.999999999999996 by 15.999999999999988
  Lower left (0.5,0.5) corner is x,y= 24105.7, 30197.1
  Upper right (980.5,314.5) corner is x,y= 39785.7, 35221.1
  DM filter is:
    x=24105.7:39785.7:#980,y=30197.1:35221.1:#314
  mkexpmap xygrid value is:
    24105.7:39785.7:#980,30197.1:35221.1:#314

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).

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 plate edges, bad pixels, etc.

(t)csh
unix% set xygrid=`pget get_sky_limits xygrid`

bash/zsh
unix% xygrid=`pget get_sky_limits xygrid`

unix% echo $xygrid
24105.7:39785.7:#980,30197.1:35221.1:#314
unix% punlearn mkexpmap
unix% pset mkexpmap xygrid=$xygrid
unix% pset mkexpmap normalize=no
unix% mkexpmap asphist.fits expmap.bin16 2.instmap mode=h

Since we set the normalize parameter to 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]; see the help file for mkexpmap for more details on this.

Figure 1: Exposure map

[Version: full-size]

Figure 1: Exposure map

unix% ds9 expmap.bin16 -cmap b

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 img.bin16 img.thresh.bin16 expfile=expmap.bin16 cut=1.5%
unix% dmimgthresh expmap.bin16 expmap.thresh.bin16 cut=1.5%

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 can be performed in one step with dmimgcalc:

unix% punlearn dmimgcalc
unix% dmimgcalc img.thresh.bin16 expmap.thresh.bin16 norm.bin16 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 norm.bin16 (Figure 2) are [photon/cm²/s/pixel].

Figure 2: Exposure-corrected image of Cas A from the HRC-S

[Version: full-size]

Figure 2: Exposure-corrected image of Cas A from the HRC-S

unix% ds9 norm.bin16 -cmap b -scale log