Determining Exposed Region Area

Getting Started

Download the sample data: 13201 (Abell 665)

unix% download_chandra_obsid 13201 

It is assumed that the chandra_repro script has been run to apply the latest calibrations to the dataset.

Make image and aspect histogram

The first step is to create an image of the counts. The fluximage script can be used to do this.

ObsId 13201 has the Abell 665 cluster on the 4 ACIS-I array chips but also has data on ACIS-6. For the purposes of this thread we only want to consider the 4 ACIS-I array chips so the data input to fluximage are filtered. This thread uses the data in the 0.5 to 7 keV range, the CSC broad band.

This thread will also require the aspect histogram which is an intermediate file created by the asphist tool which is run by fluximage. To retain these files, the cleanup parameter is set to no.

unix% fluximage infile='acisf13201_repro_evt2.fits[ccd_id=0:3]' outroot='abell_665' \
  bands='broad' binsize='2' clob+ verb=1 cleanup=no
...

The following files were created:

 The clipped counts image is:
     abell_665_broad_thresh.img

 The binned counts image is:
     abell_665_broad.img

 The aspect histograms are:
     abell_665_0.asphist
     abell_665_1.asphist
     abell_665_2.asphist
     abell_665_3.asphist
  ...

 The exposure-corrected image is:
     abell_665_broad_flux.img

fluximage creates several outputs. Those needed in the remainder of this thread are highlighted. The image of the broad-band counts is shown in Figure 1.

Figure 1: ObsID 13201, Abell 665, broad band counts

[Version: full-size]

Figure 1: ObsID 13201, Abell 665, broad band counts

The 0.5 to 7.0 keV counts image for ObsID 13201, only the 4 -I array chips are included. Data are binned by 2. Colors are logarithmically scaled and the image has a 3pixel Gaussian smoothing applied.

The white ellipse is the region used in the remainder of this thread. It is described by

ellipse(8:30:55.337,+65:52:01.94,226.249",260.076",349.593)

The region between the 4 CCDs is easily visible in this image: the 'dim' regions where number of counts decreases. The number of counts does not go to zero everywhere in these regions. The dither of the space craft acts to fill in the region between the chips to provide at least partial exposure.

The white ellipse shown in Figure 1 is an approximate boundary for the Abell 665 cluster. It has been saved in ds9 format in world coordinates. The simple geometric area of this region can be computed with several CIAO tools. In this thread dmextract is used.

unix% cat abell_665.reg 
# Region file format: DS9 version 4.1
global color=green dashlist=8 3 width=1 font="helvetica 10 normal roman" select=1 highlite=1 dash=0 fixed=0 edit=1 move=1 delete=1 include=1 source=1
fk5
ellipse(8:30:55.337,+65:52:01.94,226.249",260.076",349.593) # color=white width=3

unix% dmextract "abell_665_broad_thresh.img[bin sky=region(abell_665.reg)]" dme_output op=generic
unix% dmlist dme_output"[cols area,cel_area]" data,clean
#  AREA                 CEL_AREA
             763676.0    184858.4672640076

The area is 763,676 physical pixels (or 184,858 arcsec²), however this is simply the geometric area of the ellipse (A = π r₁r₂)

% echo "3.141592*226.249*260.076" | bc -l
184857.352021759008

The small difference between CEL_AREA and the value above is because the the AREA is computed as the number of pixels (scaled for bin size) whose center is inside the region, which is why it is an integer value.

Since this region includes the gap between chips, only a fraction of this total area will actually be exposed by the detector. To determine this fraction we need to use the dither_region tool.

Run dither_region

The dither_region tool takes a sky region and maps it through the aspect histogram (or more generally the aspect solution) to determine what fraction of the region covers the active part of the detector. It requires several inputs

• aspect histogram file: These were created earlier by fluximage. When the good time intervals for each chip are nearly identical, the aspect histograms will be indistinguishable. Comparing the files with dmdiff is not useful, as it will not be able to distinguish files that while are not identical are sufficiently similar. Instead it is enough to check the LIVTIME keywords for each CCD to make sure they are nearly the same

  % dmlist abell_665_broad_thresh.img header,clean | grep LIVTIME
  LIVTIME3                 48706.8784010260 [s]       Livetime
  LIVTIME2                 48709.93783850 [s]         Livetime
  LIVTIME1                 48697.4975006140 [s]       Livetime
  LIVTIME0                 48700.5569679720 [s]       Livetime
  LIVTIME6                                0 [s]       Livetime

  Since the LIVETIME values for each of the CCDs is approximately equal, any one of the per-chip aspect histograms can be used. (The LIVETIME for CCD_ID=6 is 0 because it was removed when fluximage was run.) If the times were very different (more than 10%), then the area computed in this way may be in error. To get a time weighted average region area in this case is less meaningful. Using the per-chip areas, AREA_CHIP_FRAC may be an option. Users should contact the CXC Helpdesk if they need assistance.

• region: The region chosen for this thread is a rough outline of the extent of the Abell 665 cluster. It was saved into a ds9 format ASCII region file in world coordinates.

• maskfile: The detector mask file is included in the secondary directory and is copied into the chandra_repro output directory. It is used to get information about any sub-array or windows that were in use during the observation.

• wcsfile: Since the region was saved in world coordinates, the wcsfile must be specified to provide the transformation necessary to convert between world and sky coordinates. It can be either an image or an event file that has the same WCS tangent plane parameters as when the region was created.

ardlib.par

chandra_repro will, by default, setup the ardlib.par parameter file with the correct bad pixel file for the observation. However, if multiple observations were processed or if some other process has changed the ardlib.par file then it needs to be reconfigured with the current observation's badpixel file.

unix% punlearn ardlib
unix% acis_set_ardlib acisf13201_repro_bpix1.fits
Updated ardlib parameter file: /home/user/cxcds_param4/ardlib.par
  AXAF_ACIS0_BADPIX_FILE -> 13201/repro/acisf13201_repro_bpix1.fits[BADPIX0]
  AXAF_ACIS1_BADPIX_FILE -> 13201/repro/acisf13201_repro_bpix1.fits[BADPIX1]
  AXAF_ACIS2_BADPIX_FILE -> 13201/repro/acisf13201_repro_bpix1.fits[BADPIX2]
  AXAF_ACIS3_BADPIX_FILE -> 13201/repro/acisf13201_repro_bpix1.fits[BADPIX3]
  AXAF_ACIS4_BADPIX_FILE -> CALDB
  AXAF_ACIS5_BADPIX_FILE -> CALDB
  AXAF_ACIS6_BADPIX_FILE -> 13201/repro/acisf13201_repro_bpix1.fits[BADPIX6]
  AXAF_ACIS7_BADPIX_FILE -> CALDB
  AXAF_ACIS8_BADPIX_FILE -> CALDB
  AXAF_ACIS9_BADPIX_FILE -> CALDB

Failure to set the correct badpixel file may result in errors in the areas computed below.

The last parameter that needs to be set is the maxpix parameter. dither_region works by first creating a grid of pixels that fills the region. Then for each grid point, the transformation from sky to chip coordinates is computed, for every aspect histogram bin. These computations are very complex and thus with many grid points the run-time of the tool can be prohibitive. The maxpix parameter specifies the maximum number of grid points that will be used to cover the area. This parameter can be used to adjust the accuracy of the calculation and trade-off the time required to get an answer.

Since we know the total area of the region is approximately 760,000 pixels, leaving maxpix=5000 (the default value) may lead to too large an error. Generally, maxpix should be no smaller than 10% of the region area. If the region has any sharp features (sources cut out, angular regions excluded) then the maxpix value should be closer to the total area. Since the region is simple in this thread, the value is set to 76000.

The dither_region parameters can now be set and the tool run. Given the large region size it may take a few minutes to complete. With the verbose set to 3 a percent-complete message is displayed.

unix% pset dither_region infile='abell_665_0.asphist' \
  region='region(abell_665.reg)' \
  outfile='dr_out.fits' \
  maskfile='acisf13201_000N003_msk1.fits' \
  wcsfile='abell_665_broad_thresh.img' \
  maxpix='76000' \
  verbose='3' \
  clobber='yes'

unix% dither_region mode=h 
dither_region
        infile = abell_665_0.asphist
        region = region(abell_665.reg)
       outfile = dr_out.fits
      maskfile = acisf13201_000N003_msk1.fits
      tolerace = INDEF
    resolution = 1.000000
        maxpix = 76000
        convex = no
       psffile = 
       gtifile = 
      geomfile = geom
     ardlibpar = ardlib
       clobber = yes
       verbose = 3
aspect solution tolerance set to 6.833333e-05 [deg]
# dither_region (CIAO): WARNING: Too many pixels, resetting resolution from 1.000000e+00 to 2.000000e+00

# dither_region (CIAO): WARNING: Too many pixels, resetting resolution from 2.000000e+00 to 3.000000e+00

# dither_region (CIAO): WARNING: Too many pixels, resetting resolution from 3.000000e+00 to 4.000000e+00

Since the total area was larger than maxpix parameter, the resolution of the pixel grid was increased until no more than maxpix points remained in the grid. The information about the "aspect solution tolerance" is not relevant in this example since the aspect histogram was input.

The output from dither_region is a table with many of the same columns as the aspect histogram file with new columns providing the FRACAREA (fraction of region area), AREA_CHIP_FRAC (the per-chip fraction), and a STATUS column with a bit-encoded value indicating why the fraction is not 100%.

unix% dmlist dr_out.fits cols
 
--------------------------------------------------------------------------------
Columns for Table Block AREAFRACTION
--------------------------------------------------------------------------------
 
ColNo  Name                 Unit        Type             Range
   1   CAH_REC                           Int4           -                    
   2   POS_OFFSET(X_OFFSET,Y_OFFSET) pix          Real8          -Inf:+Inf            Sky Position
   3   ROLL                 deg          Real8          -Inf:+Inf            Roll angle
   4   FRACAREA                          Real8          -Inf:+Inf            Fraction area
   5   AREA_CHIP_FRAC[10]                Real8(10)      -Inf:+Inf            Region Area Fraction per chip
   6   STATUS                            Bit[8]                              Why fraction < 1

The aspect histogram contains the total duration at each location and as is shown in Figure 2 the duration in each aspect histogram bin varies significantly.

Figure 2: Image of Aspect Histogram

[Version: full-size]

Figure 2: Image of Aspect Histogram

The aspect histogram converted into a 2D image. The X-axis is offset in RA and the Y-axis is the offset in Declination. The pixel value is the total duration at each RA and Dec offset. The data have been integrated over the 3rd dimension, Roll, which typically all fall into a single aspect histogram bin.

The duration value in the aspect histogram can be used to weight the image pixel values in dmcopy as shown below:

unix% dmstat abell_665_0.asphist"[cols x_offbin,y_offbin]"
POS_OFFBIN(X_OFFBIN, Y_OFFBIN)[bin]
    min:	( 1 1 )	      @:	( 210 1119 )
    max:	( 44 47 )	      @:	( 944 1 )
...
unix% dmcopy "abell_665_0.asphist[bin x_offbin=0:45:1,y_offbin=0:48:1;duration]" asphist.img    

The aspect solution tends to spend more time at the edge and especially in the corners. Therefore the amount of area lost when the aspect is at those locations should be more heavily weighted then in the middle of the pattern.

The FRACAREA values need to be weighted by the DURATION to get a good estimate of the exposure correction that should be applied to the region due to parts of the region spending nonzero time off the active area of the detector. Since the rows in the dither_region match those in the aspect histogram, their data can easily be combined in Python

unix% python
>>> from pycrates import read_file
>>> dr = read_file('dr_out.fits')
>>> asp = read_file('abell_665_0.asphist')
>>> fracarea = dr.get_column("fracarea").values
>>> duration = asp.get_column("duration").values
>>> print(sum(fracarea*duration) / sum(duration))
0.8825479009281537

This could also be accomplished by using dmpaste to combine the aspect histogram with the dither_region output; dmtcalc to compute the duration weighted fracarea, dmstat to compute the sums.

The result is that approximately 88% of the region area is being imaged by a good and active part of the detector during the observation. So if one wanted to compare counts in this region on this observation compared to another, the counts would be corrected by an area of

763676 * 0.8825 = 673944.07 pixels

Account for Surface Brightness

Based on the analysis performed in the previous section, the area of the region is on average only exposed 88% of the time. However, the cluster's surface brightness, Figure 1, is clearly not uniformly distributed over the ellipse. So while 12% of the geometric area is off the detector, the amount of flux lost due to the shape of the object is still unknown.

A model of the object can be input to dither_region via the psffile parameter. For a point source, the PSF is the model of the source as observed by the detector. For an extended source, such as this cluster, the model of the object is more complicated. One way to approximate the source model is to use the exposure corrected, or flux image that was created by fluximage. The image should be smoothed so that each pixel represents the expected surface brightness rather than the observed surface brightness.

unix% aconvolve abell_665_broad_flux.img sm_flux.img 'lib:gaus(2,5,5,3,3)' method=fft clobber=yes

Smoothing exposure corrected images

For the purposes of brevity, the exposure corrected image was smoothed. The more correct approach would be to smooth the original counts image and the exposure map separately, and then create the smoothed, exposure corrected image by dividing those two images and applying an exposure threshold to remove areas around the edge of the detector that were not exposed.

dither_region requires that the psffile be normalized. The image should so also be filtered so that just the cluster is included.

unix% dmstat "sm_flux.img[sky=region(abell_665.reg)]"  centroid- sig- med-
CONVOLVE
    min:	1.2040897168e-09 	      @:	( 4181.5 4739.5 )
    max:	3.2900709357e-07 	      @:	( 4019.5 4035.5 )
   mean:	2.8272640293e-08 
    sum:	0.0054592084616 
   good:	190919 
   null:	140677 
unix% dmimgcalc "sm_flux.img[sky=region(abell_665.reg)]" none model.fits op="imgout=(img1/0.005397784212)" clob+

The resulting image is shown in Figure 3. Note that the image size has been reduced to the bounding box around the region. This will help to improve dither_region's speed.

Figure 3: Exposure corrected, smoothed image model

[Version: full-size]

Figure 3: Exposure corrected, smoothed image model

The exposure corrected, smoothed, and normalized image model of the Abell 665 cluster. This image will be input to dither_region as the psffile to compute the amount of flux lost due to the object dithering off the edge of the detector.

The presence of the point-like sources in this image will not greatly affect the results given their location relative to the detector boundary. The Diffuse Emission thread shows how the point-like source can be removed with the dmfilth tool. More advanced smoothing algorithms may also be used such as csmooth or dmimgadapt.

unix% pset dither_region psffile=model.fits

and the tool is run.

unix% dither_region outfile=dr_with_model.fits mode=h verbose=0

unix% dmlist dr_with_model.fits cols
 
--------------------------------------------------------------------------------
Columns for Table Block AREAFRACTION
--------------------------------------------------------------------------------
 
ColNo  Name                 Unit        Type             Range
   1   CAH_REC                           Int4           -                    
   2   POS_OFFSET(X_OFFSET,Y_OFFSET) pix          Real8          -Inf:+Inf            Sky Position
   3   ROLL                 deg          Real8          -Inf:+Inf            Roll angle
   4   FRACAREA                          Real8          -Inf:+Inf            Fraction area
   5   PSFFRAC                           Real8          -Inf:+Inf            PSF fraction
   6   AREA_CHIP_FRAC[10]                Real8(10)      -Inf:+Inf            Region Area Fraction per chip
   7   PSF_CHIP_FRAC[10]                 Real8(10)      -Inf:+Inf            PSF Fraction per chip
   8   STATUS                            Bit[8]                              Why fraction < 1

The output file now contains additional columns related to the psffile. The "psffrac" is really the fraction of the psffile image which in this example is the surface brightness of the cluster. The duration weighted average can be computed as above

unix% python
>>> from pycrates import read_file
>>> dr = read_file('dr_with_model.fits')
>>> asp = read_file('abell_665_0.asphist')
>>> psffrac = dr.get_column("psffrac").values
>>> duration = asp.get_column("duration").values
>>> print(sum(psffrac*duration) / sum(duration))
0.868472329468226

Which shows that when the surface brightness is considered, 86.8% of the cluster's flux is, on average, incident on an active part of the detector. This number is slightly lower than the 88% number before since the gap between the CCDs does fall on a somewhat bright part of the model image.

Summary

This thread has demonstrated how to correctly account for the amount of area that falls off the detector during a typical Chandra observation using the dither_region tool. This can be especially useful when comparing the counts from one observation with those of the same object using the same region in another observation that is possibly offset in either pointing or roll, or both.

For regions where the exposure map value does not vary significantly from pixel to pixel (ie away from the detector edges), the mean of the exposure map value in the region may also be used to normalize the counts when comparing different observations.

The geometric area of the region is not energy dependent; however, at the point where the model is introduced, the energy bands and spectral shapes of the model and the data need to be consistent to meaningful results.

Finally, the resolution, or maxpix, parameter may require some experimentation to get good results. If run-time is not an issue, then it should be set to something close to the full area of the region. This will be needed when the region is especially complex or if precision is critical in the analysis. If the analysis does not require great precision, then using a smaller value for maxpix will greatly speed-up the dither_region tool.