Skip to the navigation links
Last modified: 5 Feb 2010
Where are the PDFs?

HRC-S Exposure Map and Fluxed Image

CIAO 4.2 Science Threads



Overview

Last Update: 5 Feb 2010 - updated for CIAO 4.2: minor changes to screen output

Synopsis:

mkexpmap generates an exposure map which may be used to convert a counts image of a source to an image in flux units. The computed exposure map is essentially an image of the effective area at each sky position, accounting for the effects of dither motion which are especially important near the edges of the detector.

The exposure map is computed from the aspect histogram file - which contains information on the aspect motion during the observation - and an instrument map - which is essentially the product of the detector quantum efficiency and the mirror effective area projected onto the detector surface.

Purpose:

To build an exposure map for an HRC-S observation, create a fluxed image, and find an approximation for the source flux.

Read this thread if:

you are working with an HRC-S observation and would like to create an exposure map for it.

Related Links:




Contents



Get Started

Sample ObsID used: 990 (HRC-S, VEGA)

File types needed: evt2; dtf1; asol1; msk1

Please ensure that you have set up ardlib to use the bad pixel file for your observation before following this thread; see the Setting the Observation-specific Bad Pixel Files thread for more information.

Download get_sky_limits

This thread uses the get_sky_limits script; for information about the script, consult the help file ("ahelp get_sky_limits"). The script is part of the CIAO Scripts distribution. The CIAO scripts package should be the following version or newer:

unix% cat $ASCDS_CONTRIB/VERSION.CIAO_scripts
17 Apr 2009

Please check that you have at least this version of the scripts package installed before continuing. If you do not have the scripts installed or need to update to a newer version, refer to the Scripts page.


Syntax note: foreach

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



Create An Image

This thread creates an exposure map for all three plates of HRC-S. Some users will only need to use chip_id=2 (the center plate), since it is most commonly used for HRC-S imaging.

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 32:

unix% dmcopy \
      "hrcf00990N004_evt2.fits[bin x=0.5:65536.5:32,y=0.5:65536.5:32][opt type=i4]" \
      990_img.fits

This creates an image that is 2048x2048; this information is used again in the Calculate the Exposure Map step. You may choose to use different binning, but make sure that you change the xygrid appropriately in that step. Due to the size of the image, the output size is set to 4 byte integer ("opt type=i4") instead of the default 2 byte integer.



Compute Exposure Map

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

unix% dmlist hrcf00990N004_evt2.fits blocks
 
--------------------------------------------------------------------------------
Dataset: hrcf00990N004_evt2.fits
--------------------------------------------------------------------------------
 
     Block Name                          Type         Dimensions
--------------------------------------------------------------------------------
Block    1: PRIMARY                        Null        
Block    2: EVENTS                         Table         9 cols x 453434   rows
Block    3: GTI                            Table         2 cols x 1        rows

In many cases, there will be more than one aspect solution file (pcad_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 pcad_asol1.lis
pcadf097387525N002_asol1.fits

unix% punlearn asphist
unix% pset asphist infile=@pcad_asol1.lis
unix% pset asphist outfile=asphist_hrcs.fits
unix% pset asphist evtfile=hrcf00990N004_evt2.fits
unix% pset asphist dtffile=hrcf00990_000N003_dtf1.fits
unix% asphist
Aspect Solution List Files (@pcad_asol1.lis): 
Aspect Histogram Output File (asphist_hrcs.fits): 
Event List Files (hrcf00990N004_evt2.fits): 
Live Time Correction List Files for HRC (hrcf00990_000N003_dtf1.fits): 

# asphist (CIAO 4.2): WARNING: skipping 84 livetime correction records (from time: 97387355.142212 to time: 97387525.292218)

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


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

unix% punlearn mkinstmap
unix% pset mkinstmap obsfile="hrcf00990N004_evt2.fits[EVENTS]"
unix% pset mkinstmap pixelgrid="1:16384:#1024,1:16384:#1024"
unix% pset mkinstmap monoenergy=1.1
unix% pset mkinstmap mode=h

Now run the tool once for each plate:

unix% foreach d ( 1 2 3 )
foreach? mkinstmap detsubsys=HRC-S${d} outfile=instmap_hrcs_${d}.fits \
                   maskfile="hrcf00990_000N003_msk1.fits[MASK${d}]"
foreach? end

3. 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 990_img.fits verbose="1"
Checking binning of image: 990_img.fits
  Image has 2048 x 2048 pixels
  Pixel size is 32 by 32
  Lower left (0.5,0.5) corner is x,y= 0.5, 0.5
  Upper right (2048.5,2048.5) corner is x,y= 65536.5, 65536.5
  DM filter is:
    x=0.5:65536.5:#2048,y=0.5:65536.5:#2048
  mkexpmap xygrid value is:
    0.5:65536.5:#2048,0.5:65536.5:#2048

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.

unix% punlearn mkexpmap
unix% pset mkexpmap asphistfile=asphist_hrcs.fits
unix% pset mkexpmap xygrid="0.5:65536.5:#2048,0.5:65536.5:#2048"
unix% pset mkexpmap useavgaspect=no
unix% pset mkexpmap normalize=no
unix% pset mkexpmap mode=h

unix% foreach d ( 1 2 3 )
foreach? mkexpmap instmapfile=instmap_hrcs_${d}.fits \
                  outfile=expmap_hrcs_${d}.fits
foreach? end
Exposure map limits: 0.000000e+00, 3.890097e+05
Writing exposure map to expmap_hrcs_1.fits
Exposure map limits: 0.000000e+00, 5.221627e+05
Writing exposure map to expmap_hrcs_2.fits
Exposure map limits: 0.000000e+00, 3.659997e+05
Writing exposure map to expmap_hrcs_3.fits

Since we set the normalize parameter to no, the exposure map has units of [cm2*s*counts/photon]. This allows us to simply divide the image by the exposure map to derive an image in units of flux ([photons/cm2/s/pixel]). If the setting had been left as yes (the default), the units of the exposure map would be [cm2*counts/photon]; see the help file for mkexpmap for more details on this.


4. Combine the Exposure Maps

The individual exposure maps are combined into a single, binned exposure map with the CIAO tool reproject_image. First we need a list of files to combine:

unix% ls expmap_hrcs_*.fits > expmap.lis
unix% cat expmap.lis
expmap_hrcs_1.fits
expmap_hrcs_2.fits
expmap_hrcs_3.fits

Now we can use this list by passing it into reproject_image as a stack. It is important to note that the method parameter is set to average; this is because we want the average exposure value for each pixel in the combined image, not the sum of all of them.

unix% punlearn reproject_image
unix% pset reproject_image infile=@expmap.lis
unix% pset reproject_image matchfile=990_img.fits
unix% pset reproject_image outfile=expmap_hrcs.fits 
unix% pset reproject_image method=average

unix% reproject_image
Input image file name (@expmap.lis):
Reference image (990_img.fits):
Output file name (expmap_hrcs.fits):
warning: DETNAM has different value...Merged...

The combined exposure map can be displayed in ds9 (Figure 1).

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



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% pset dmimgthresh infile=990_img.fits
unix% pset dmimgthresh outfile=990_img_clean.fits
unix% pset dmimgthresh expfile=expmap_hrcs.fits
unix% pset dmimgthresh cut=1.5%
unix% pset dmimgthresh value=0.0
unix% dmimgthresh
Input dataset/block specification (990_img.fits):
Output dataset/block specification (990_img_clean.fits):

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 [cm2*s*counts/photon] since it was created by projecting the instrument map (in [cm2*counts/photon]) onto the tangent plane of the observation. To create an image in units of [photon/cm2/s/pixel], we simply need to divide by the exposure map. This can be performed in one step with dmimgcalc:

unix% punlearn dmimgcalc
unix% pset dmimgcalc infile=990_img_clean.fits
unix% pset dmimgcalc infile2=expmap_hrcs.fits
unix% pset dmimgcalc outfile=990_img_norm.fits
unix% pset dmimgcalc operation=div
unix% dmimgcalc
Input file #1 (990_img.fits): 
Input file #2 (expmap_hrcs.fits): 
output file  (990_img_norm.fits): 
arithmetic operation (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 990_img_norm.fits (Figure 2) are [photon/cm2/s/pixel].

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



Calculate the Source Flux

Since the units of the fluxed image are [photon/cm2/s/pixel], adding up the pixel values around a source results in the source flux in [photon/cm2/s]. Note that this flux is an approximation - as discussed in An Introduction to Exposure Maps (PS, 12pp) - since a spectral shape was assumed when using mkinstmap (in this example, a monochromatic source).

Using the source region "flux.reg":

unix% cat flux.reg
# Region file format: CIAO version 1.0
circle(36496.5,29856.5,235)

the flux can be calculated with either dmstat:

unix% dmstat infile="990_img_norm.fits[sky=region(flux.reg)]" centroid=no
990_img_norm.fits
    min:        2.2780808649e-05              @:        ( 36592.5 29648.5 )
    max:        0.082072928548        @:        ( 36496.5 29872.5 )
   mean:        0.0048851800029
  sigma:        0.014049377247
    sum:        0.81093988049
   good:        166
   null:        74

or dmextract:

unix% dmextract infile="990_img_norm.fits[bin sky=@flux.reg]" outfile="source_flux.fits"

unix% dmlist source_flux.fits"[cols COUNTS]" data
 
--------------------------------------------------------------------------------
Data for Table Block HISTOGRAM
--------------------------------------------------------------------------------
 
ROW    COUNTS                     
 
     1     0.81093988048633

Since the input to dmextract was a fluxed image, not an event list, the COUNTS column actually reports the total flux (in [photon/cm2/s]) for the source region. While slightly more involved, the dmextract method can be used on multiple sources in a single command, and the results are conveniently stored in a table.

To compute robust source intensity quantities (net counts, source rate, photon flux, energy flux) and the related confidence intervals, use the aprates tool. The Compute Net Counts, Rate, or Flux for Point Sources thread shows how to run aprates.

If you are working with event lists, the eff2evt tool can be used to compute the approximate flux, and calculate the QE and Effective Area for sources. The Calculate the Flux for a Position thread describes how to use this tool.




Parameters for /home/username/cxcds_param/asphist.par


        infile = @pcad_asol1.lis  Aspect Solution List Files
       outfile = asphist_hrcs.fits Aspect Histogram Output File
       evtfile = hrcf00990N004_evt2.fits Event List Files
       dtffile = hrcf00990_000N003_dtf1.fits Live Time Correction List Files for HRC
      (geompar = geom)            Parameter file for Pixlib Geometry files
       (res_xy = 0.5)             Aspect Resolution x and y in arcsec
     (res_roll = 600.)            Aspect Resolution roll in arcsec
      (max_bin = 10000.)          Maximal number of bins
      (clobber = no)              Clobber output
      (verbose = 0)               Verbose
         (mode = ql)              
    


Parameters for /home/username/cxcds_param/reproject_image.par


        infile = @expmap.lis      Input image file name
     matchfile = 990_img.fits     Reference image
       outfile = expmap.fits_reproimg Output file name
    (resolution = 1)               Number of point per side to evalute
       (method = average)         Average value
    (coord_sys = world)           Coordinate system to match images in
    (lookupTab = ${ASCDS_CALIB}/dmmerge_header_lookup.txt -> /soft/ciao/data/dmmerge_header_lookup.txt) lookup table
      (clobber = no)              Clobber existing files
      (verbose = 0)               Tool verbosity
         (mode = ql)
    


Parameters for /home/username/cxcds_param/dmimgthresh.par


        infile = 990_img.fits     Input dataset/block specification
       outfile = 990_img_clean.fits Output dataset/block specification
      (expfile = expmap_hrcs.fits) Exposure map file
          (cut = 1.5%)            Threshold value
        (value = 0)               Replacement value
      (verbose = 0)               Debug Level(0-5)
      (clobber = no)              Clobber existing file
         (mode = ql)
    


Parameters for /home/username/cxcds_param/dmimgcalc.par


# parameter file for dmimgcalc
        infile = 990_img.fits     Input file #1
       infile2 = expmap_hrcs.fits Input file #2
       outfile = 990_img_norm.fits output file 
     operation = div              arithmetic operation
       (weight = 1)               weight for first image
      (weight2 = 1)               weight for second image
    (lookupTab = ${ASCDS_CALIB}/dmmerge_header_lookup.txt -> /soft/ciao/data/dmmerge_header_lookup.txt) lookup table
      (clobber = no)              delete old output
      (verbose = 0)               output verbosity
         (mode = ql)              
    

History

14 Dec 2007 new for CIAO 4.0
08 Sep 2008 combine the exposure maps with reproject_image instead of dmregrid
24 Oct 2008 get_sky_limits v1.13 (fixes a rare segmentation fault and adds the pixel size in sky coordinates to the screen output.)
13 Feb 2009 updated for CIAO 4.1: images are inline; minor change to flux in screen output; run dmimgthresh before dmimgcalc in Normalize the Image by the Exposure Map section
16 Mar 2009 added link to aprates thread in Calculate the Source Flux section, removed listing of ERR_COUNTS column
06 May 2009 check the version of the CIAO scripts package instead of the individual script
05 Feb 2010 updated for CIAO 4.2: minor changes to screen output

Return to Threads Page: Top | All | Imag

Where are the PDFs?
Last modified: 5 Feb 2010