Determining Exposed Region Area
CIAO 4.15 Science Threads
What is the area of the region?
This is a deceptively simple question. Computing the geometric area of primitive shapes such as circles and ellipses is trivial and even the area of complex regions with multiple shapes being included and excluded is generally doable. The deceptiveness is that often what is meant is What is the area of the region covered by the detector?. Or conversely, What amount of the area is off the detector or is covered by bad pixels? Since Chandra dithers, a region defined on the sky is exposed to different parts of the detector as a function of time. So the answer to these questions often becomes, When? If the region is small then it may not overlap any bad pixels or detector edges in which case the geometric area is correct; however, if the region is large (extended source or modestly far-off-axis point source), the area should be corrected for the fraction exposed.
Tools such as dmextract can be used to compute the area of a region for use in spectral analysis, via the BACKSCAL keyword. If the region dithers off chip or if it dithers across bad pixels or columns, the amount of geometric area lost is accounted for when making a weighted ARF.
However, if the type of analysis does not include the ARF, such as comparing counts or aperture photometry, then this geometric correction needs to be accounted for in other ways.
This thread will show how to use the dither_region tool to compute the fraction of the region area that is on a good part of the detector and how to include the correction in the analysis.
To compute the area of a region accounting for the loss due to dithering off-chip or over bad pixels.
- Multiple Chip ACIS Exposure Map and Exposure-corrected Image
- HRC-I Exposure Map and Exposure-corrected Image
- HRC-S Exposure Map and Exposure-Corrected Image
- fluximage help file
- dither_region help file
- DM region syntax
Last Update: 7 Feb 2022 - Review for CIAO 4.14. Updated for Repro5/CALDB4.9.6
- Getting Started
- Make image and aspect histogram
- Run dither_region
- Account for Surface Brightness
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
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 arcsec2), however this is simply the geometric area of the ellipse (A = π r1r2)
% 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.
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.
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 Real8(10) -Inf:+Inf Region Area Fraction per chip 6 STATUS Bit 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
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
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
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
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.dither_regionpsffrac
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 Real8(10) -Inf:+Inf Region Area Fraction per chip 7 PSF_CHIP_FRAC Real8(10) -Inf:+Inf PSF Fraction per chip 8 STATUS Bit 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.
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.
|04 Oct 2013||Initial version.|
|10 Dec 2013||Review for CIAO 4.6. Minor edits.|
|22 Dec 2014||Reviewed for CIAO4.7; no changes.|
|07 Feb 2022||Review for CIAO 4.14. Updated for Repro5/CALDB4.9.6|