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.

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

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.

Figure 1: Exposure map

[Version: full-size]

Figure 1: Exposure map