AHELP for CIAO 4.2 |
## aprates |
Context: tools |

## Synopsis

Compute net counts, rate, and/or flux, plus confidence regions, for a point source, using data obtained from source and background apertures in event lists, images, and exposure maps.

## Syntax

aprates n A_s alpha T_s E_s eng_s flux_s m A_b beta T_b E_b eng_b flux_b conf outfile [resolution] [pdf] [max_counts] [itermax] [nsigma] [pmin]

## Description

The aprates tool computes values and bounds for source intensity quantities (net counts, source rate, photon flux, energy flux) using counts and exposure data obtained in source and background apertures. All input is obtained from the input parameter file; output is written to an output parameter file.

The tool determines intensity values by solving the pair of simultaneous linear equations

n = f * S + b

m = g * S + r * b

for the intensity quantity S. Here, n and m represent the (observed) total counts in the source and background apertures, and b represents the (model) background counts in the source aperture. The quantities f and g describe the fractions of the source point response function contained in the source and background apertures, but, depending on the intensity quantity of interest, may also include factors that convert the intensity quantity S into counts. Finally, the quantity r is a "backscale" parameter, to scale the background b from source to background aperture. The tool computes f, g, and r from more basic user-input data, as described in the examples below.

To determine confidence bounds, the tool computes the Bayesian background-marginalized posterior probability distribution function (PDF) for S, assuming non-informative priors for the intensities in the source and background apertures. The mode of this PDF is determined, and the lower and upper bounds of the confidence region are determined by summing values of the PDF alternately above and below the mode until the desired confidence level is attained. If the summation below the mode reaches the PDF for S=0, the lower confidence bound is set to 0 and the summation continues for PDF values above the mode only. If the mode itself it 0, lower confidence bound is set to INDEF and the summation proceeds for PDF values above the mode only. In either of these cases, the upper confidence bound may be considered an upper limit.

For more detail on the algorithms, please see the documents

http://cxc.harvard.edu/csc/memos/files/Kashyap_xraysrc.pdf http://cxc.harvard.edu/csc/memos/files/Primini_significance.pdf http://cxc.harvard.edu/csc/memos/files/Primini_apflux_spec_II.pdf

At the minimum, the user must input the number of counts in the source and background apertures (n, m), their geometric areas (A_s, A_b), and their psf fractions (alpha, beta), setting all other input parameters to 1. One then obtains net source counts and errors (the output file will contain entries for rates and fluxes, but they repeat the net count values.

By setting T_s and T_b to the exposure times in the source and background apertures, one enables computation of net source rate and errors.

By setting E_s and E_b to the average effective exposures (in cm^2-s) in the two apertures, one enables computation of photon flux and errors.

By setting eng_s and eng_b to the average photon energies (in ergs) in the two apertures, one enables computation of energy flux and errors.

Finally, by setting flux_s and flux_b to the average of photon energy/effective exposure (in ergs/cm^2-s) in the two apertures, one enables an alternate computation of energy flux and errors. This computation is preferable to the one using eng_s and eng_b because it more properly weights the contribution of higher energy photons, which, though few, may dominate the energy flux.

The following table indicates the inputs required to produce the desired outputs. Including additional inputs, e.g., T_s and T_b for photon flux, will not affect calculation of the desired output.

#### Required Input Parameters for Different Desired Output Quantities

Net Counts | Rate | Photon Flux | Energy Flux1 | Energy Flux2 | |
---|---|---|---|---|---|

n | x | x | x | x | x |

m | x | x | x | x | x |

alpha | x | x | x | x | x |

beta | x | x | x | x | x |

A_s | x | x | x | x | x |

A_b | x | x | x | x | x |

T_s | x | ||||

T_b | x | ||||

E_s | x | x | |||

E_b | x | x | |||

eng_s | x | ||||

eng_b | x | ||||

flux_s | x | ||||

flux_b | x |

### Things to Watch Out For:

### 1.) Extracting Aperture Quantities

A number of different CIAO tools can be used to determine the area and number of events in apertures. These tools may yield slightly different results, depending on whether one starts with an event list or image. Results from event lists are more accurate, since event locations are typically known to finer resolution that a pixel size, and areas can be determined analytically for simple apertures. For images, aperture counts and areas are determined from those pixels whose centers fall within the aperture. We recommend the use of dmextract if an event list is available.

Unfortunately, determination of effective exposure and psf fractions in apertures will almost always require the use of exposure map and psf images, leading to possible inaccuracies in the determination of these quantities. For apertures containing many image pixels, this is a negligible effect, because exposure maps typically vary smoothly, and thus inclusion or exclusion of particular exposure map pixels has little effect on the average in the aperture. Similarly, the psf fraction change due to inclusion or exclusion of a few pixels at the aperture edge is likely to be small.

This will not be the case, however, if aperture sizes are small compared to image pixel sizes. In those cases, it is recommended that the user repixelate the psf and exposure map images to a finer scale.

### 2.) Computing Model-Independent Broad-Band Photon and Energy Fluxes

In order to use the Bayesian formalism to estimate these quantities, we modify psf fractions f and g in the equations above by flux conversion factors, to convert S from flux to counts. For example, for photon fluxes, we multiply f and g by the mean effective exposure (in cm^2-sec/photon from an unnormalized exposure map) in the source and background apertures, respectively (see example 3). For energy fluxes, we offer two alternatives - multiplying by the mean effective exposure as before and dividing by the mean energy per event in the aperture (example 4), or dividing by the mean energy flux per event (example 5). All suffer from possible errors and should be considered approximate. In broad energy bands, in which the effective area may vary significantly with energy, the energy or spectrum used to compute the exposure map may strongly affect the effective exposure in the apertures. Computing mean flux per event can mitigate this, but in the case of few events, a single high-energy event can dominate the band flux, and hence the flux uncertainties can be large.

Users are cautioned to treat the broad-band fluxes computed in aprates with care. As always, spectral fitting will yield the most accurate results, although this approach is probably unrealistic in the low count regime for which aprates was designed. In that case, the simplest and most reliable approach is to determine the normalization for as assumed spectrum, using Sherpa or XSPEC.

### Data for Examples

The examples below illustrate both the computation and use of all the input parameters. In all the examples, the aprates inputs are determined from data files downloaded from the Chandra Source Catalog archive for a single point source in a single observation (OBSID 313). These data include an event list (regevt3), an unnormalized exposure map (regexp3), block=1 image of the point spread function (psf3) and FITS region file (reg3). The region file has separate extensions for the source and background apertures.

## Example 1

Compute source counts n, background counts m, source aperture area A_s and background aperture area A_b: dmextract \ infile="acisf00313_000N001_r0001_regevt3.fits[energy=500:7000]\ [bin sky=region(acisf00313_000N001_r0001_reg3.fits)]" \ outfile=- \ bkg="acisf00313_000N001_r0001_regevt3.fits[energy=500:7000]\ [bin sky=region(acisf00313_000N001_r0001_reg3.fits[bkgreg])]" |\ dmlist -"[col counts, area,bg_counts,bg_area]" data,clean #COUNTS AREA BG_COUNTS BG_AREA 141.0 196.6885375977 9.0 4720.50 Compute alpha: dmstat acisf00313_000N001_r0001b_psf3.fits"[sky=region(acisf00313_000N001_r0001 _reg3.fits)]" centroid=no verbose=0 pget dmstat out_sum 0.953 Compute beta: dmstat acisf00313_000N001_r0001b_psf3.fits"[sky=region(acisf00313_000N001_r0001 _reg3.fits[bkgreg])]" centroid=no verbose=0 pget dmstat out_sum 0.029 Run aprates to compute net counts: aprates n=141 m=9 A_s=196.689 A_b=4720.5 alpha=0.953 beta=0.029 T_s=1 E_s=1 eng_s=1 flux_s=1 T_b=1 E_b=1 eng_b=1 flux_b=1 outfile=aprates_out.par conf=0.68 pget aprates_out.par src_cnts src_cnts_err_lo src_cnts_err_up 147.748 135.396 160.224

In this example, aprates is used to determine net source counts. The dmextract tool is first used to compute the total events and areas in the source and background apertures, using the regevt3 and reg3 files. Next, dmstat is used to compute the psf fractions in the apertures using the psf3 image. Finally, aprates is run to determine the net source counts and bounds for the 68% confidence region. Here, the net counts value is 147.748, and the 68% confidence region extends from 135.396 to 160.224.

## Example 2

Compute T_s and T_b: dmkeypar acisf00313_000N001_r0001_regevt3.fits LIVETIME echo+ 5977.74 Run aprates to compute net rate: aprates n=141 m=9 A_s=196.689 A_b=4720.5 alpha=0.953 beta=0.029 T_s=5977.74 E_s=1 eng_s=1 flux_s=1 T_b=5977.74 E_b=1 eng_b=1 flux_b=1 outfile=aprates_out.par conf=0.68 pget aprates_out.par src_rate src_rate_err_lo src_rate_err_up 0.0247163 0.02265 0.0268035

Next, aprates is used to compute source rate and bounds. The exposure time in the apertures is set to the LIVETIME in the event list header. Note, exposure times in source and background apertures need not be the same.

## Example 3

Compute E_s: dmstat acisf00313_000N001_r0001b_regexp3.fits"[sky=region(acisf00313_000N001_r0 001_reg3.fits)]" centroid=no verbose=0 pget dmstat out_mean 2308338.2013 Compute E_b: dmstat acisf00313_000N001_r0001b_regexp3.fits"[sky=region(acisf00313_000N001_r0 001_reg3.fits[bkgreg])]" centroid=no verbose=0 pget dmstat out_mean 2058975.0111 Run aprates to compute net photon flux: aprates n=141 m=9 A_s=196.689 A_b=4720.5 alpha=0.953 beta=0.029 T_s=1 E_s=2308338.2013 eng_s=1 flux_s=1 T_b=1 E_b=2058975.0111 eng_b=1 flux_b=1 outfile=aprates_out.par conf=0.68 pget aprates_out.par photflux_aper photflux_aper_err_lo photflux_aper_err_up 6.40e-5 5.86e-5 6.94e-5

Next, aprates is used to compute photon flux (photons/cm^2-sec) and bounds. The average exposure (in cm^2-sec) in the apertures is determined from the unnormalized exposure map (regexp3) using dmstat.

## Example 4

Compute eng_s: dmtcalc "acisf00313_000N001_r0001_regevt3.fits[energy=500:7000,sky=region(acisf0 0313_000N001_r0001_reg3.fits)]" - expression="energy=1.6e-12*energy" | dmstat "-[cols energy]" verbose=0 pget dmstat out_mean 2.56e-09 Compute eng_b: dmtcalc "acisf00313_000N001_r0001_regevt3.fits[energy=500:7000,sky=region(acisf0 0313_000N001_r0001_reg3.fits[bkgreg])]" - expression="energy=1.6e-12*energy" | dmstat "-[cols energy]" verbose=0 pget dmstat out_mean 4.27e-09 Run aprates to compute net energy flux: aprates n=141 m=9 A_s=196.689 A_b=4720.5 alpha=0.953 beta=0.029 T_s=1 E_s=2308338.2013 eng_s=2.56e-09 flux_s=1 T_b=1 E_b=2058975.0111 eng_b=4.27e-09 flux_b=1 outfile=aprates_out.par conf=0.68 pget aprates_out.par flux_aper flux_aper_err_lo flux_aper_err_up 1.63e-13 1.50e-13 1.77e-13

Compute energy flux (ergs/cm^2-sec) and bounds. Event energies are first converted to ergs using dmtcalc, and the average event energies in the apertures are then computed using dmstat. Note, when using this option for computing energy flux, both average exposure (E_s, E_b) and average event energy (eng_s, eng_b) must be input.

## Example 5

Compute flux_s: eff2evt \ "acisf00313_000N001_r0001_regevt3.fits[energy=500:7000][sky=region(acisf 00313_000N001_r0001_reg3.fits)]" - | \ dmstat -"[cols flux]" verbose=0 pget dmstat out_mean 9.81e-16 Compute flux_b: eff2evt \ "acisf00313_000N001_r0001_regevt3.fits[energy=500:7000][sky=region(acisf 00313_000N001_r0001_reg3.fits[bkgreg])]" - | \ dmstat -"[cols flux]" verbose=0 pget dmstat out_mean 2.23e-15 Run aprates to compute net energy flux: aprates n=141 m=9 A_s=196.689 A_b=4720.5 alpha=0.953 beta=0.029 T_s=1 E_s=1 eng_s=1 flux_s=9.81e-16 T_b=1 E_b=1 eng_b=1 flux_b=2.23e-15 outfile=aprates_out.par conf=0.68 pget aprates_out.par eflux_aper eflux_aper_err_lo eflux_aper_err_up 1.44e-13 1.28e-13 1.57e-13

Compute energy flux (ergs/cm^2-sec) and bounds. In this option, eff2evt is used to compute energy flux for each event, and add it as a separate column in the event list. The dmstat tool is then used to compute the average of this quantity in the apertures.

## Parameters

name | type | ftype | def | units | reqd |
---|---|---|---|---|---|

n | integer | yes | |||

A_s | real | pixels or square arcsec | yes | ||

alpha | real | 1.0 | yes | ||

T_s | real | 1.0 | sec | yes | |

E_s | real | 1.0 | cm^2 sec | yes | |

eng_s | real | 1.0 | ergs | yes | |

flux_s | real | 1.0 | ergs/cm^2/sec | yes | |

m | integer | yes | |||

A_b | real | pixels or square arcsec | yes | ||

beta | real | 0.0 | yes | ||

T_b | real | 1.0 | sec | yes | |

E_b | real | 1.0 | cm^2 sec | yes | |

eng_b | real | 1.0 | ergs | yes | |

flux_b | real | 1.0 | ergs/cm^2/sec | yes | |

conf | real | 0.68 | yes | ||

outfile | file | output | yes | ||

resolution | real | 0.01 | |||

string | alternate | ||||

max_counts | integer | 50 | |||

itermax | real | 100 | |||

nsigma | real | 5.0 | |||

pmin | real | 1.0e-5 |

## Detailed Parameter Descriptions

####
Parameter=n` (integer required)`

*
Number of counts in source aperture
*

####
Parameter=A_s` (real required units=pixels or square arcsec)`

*
Geometric area of source aperture. Either square arcsec or pixels are
allowed, as long as the units agree with those of A_b.
*

####
Parameter=alpha` (real required default=1.0)`

*
PSF Fraction in source aperture. The default value of 1.0 corresponds to a
"perfect" aperture, i.e., one which includes all source counts.
*

####
Parameter=T_s` (real required default=1.0 units=sec)`

*
Exposure time in source aperture. May be set to 1 if computation of source
rate is not desired.
*

####
Parameter=E_s` (real required default=1.0 units=cm^2 sec)`

*
Average exposure in source aperture. May be set to 1 if
computation of photon flux or energy flux (option 1) is
not desired.
*

The exposure value may be obtained with the tool dmstat from an unnormalized exposure map, as shown in Example 3.

####
Parameter=eng_s` (real required default=1.0 units=ergs)`

*
Average energy of events in source aperture. May be set to 1 if computation
of energy flux (option 1) is not desired.
*

####
Parameter=flux_s` (real required default=1.0 units=ergs/cm^2/sec)`

*
Average flux from events in source aperture. May be set to 1 if computation
of energy flux (option 2) is not desired.
*

####
Parameter=m` (integer required)`

*
Number of counts in background aperture
*

####
Parameter=A_b` (real required units=pixels or square arcsec)`

*
Geometric area of background aperture. Either square arcsec or pixels are
allowed, as long as the units agree with those of A_s.
*

####
Parameter=beta` (real required default=0.0)`

*
PSF Fraction in background aperture. The default value of 0.0 corresponds to
a "perfect" aperture, i.e., one which includes no source counts.
*

####
Parameter=T_b` (real required default=1.0 units=sec)`

*
Exposure time in source aperture. May be set to 1 if computation of source
rate is not desired.
*

####
Parameter=E_b` (real required default=1.0 units=cm^2 sec)`

*
Average exposure in background aperture. May be set to 1 if
computation of photon flux or energy flux (option 1) is
not desired.
*

The exposure value may be obtained with the tool dmstat from an unnormalized exposure map, as shown in Example 3.

####
Parameter=eng_b` (real required default=1.0 units=ergs)`

*
Average energy of events in background aperture. May be set to 1 if computation
of energy flux (option 1) is not desired.
*

####
Parameter=flux_b` (real required default=1.0 units=ergs/cm^2/sec)`

*
Average flux from events in background aperture. May be set to 1 if computation
of energy flux (option 2) is not desired.
*

####
Parameter=conf` (real required default=0.68)`

*
Desired confidence interval
*

####
Parameter=outfile` (file required filetype=output)`

*
Output file name
*

####
Parameter=resolution` (real default=0.01)`

*
Resolution of PDF function. The PDF is evaluated on a grid with bin width
equal to this fraction of the mode.
*

####
Parameter=pdf` (string default=alternate)`

*
Which statistics to use in evaluating the PDF. The "alternate" method uses
Bayesian statistics to compute the background-marginalized posterior
probability distribution for the intensity. The "gaussian" method uses
standard Gaussian statistics.
*

####
Parameter=max_counts` (integer default=50)`

*
Max total counts (n+m) before switching to Gaussian statistics in evaluating
the PDF. If the number of aperture counts is large, the algorithm for
computing the PDF using Bayesian statistics may run for a very long time
before converging, and result in values negligibly different from those
obtained from Gaussian statistics.
*

####
Parameter=itermax` (real default=100)`

*
Maximum number of iterations allowed in computing grid for PDF. The PDF is
computed on a grid whose size is determined by the mode and width of the
distribution. If the ratio between the minumum PDF value and the value at
the mode is greater than the parameter pmin, the grid is increased and the
PDF computation is iterated.
*

####
Parameter=nsigma` (real default=5.0)`

*
Number of sigma from the mode in either direction, used in the initial
estimate of the PDF grid size.
*

####
Parameter=pmin` (real default=1.0e-5)`

*
Threshold of the ratio min(PDF)/PDF[mode]. If the actual value is greater
than this value, the size of the PDF grid is increased and the PDF
computations are repeated on the larger grid, until itermax iterations have
been reached.
*

## Bugs

The max_counts parameter provides a means to avoid excessive computational time when differences between Bayesian and Gaussian approaches are negligible. However, The sum n+m is a poor indicator of the transition between the two regimes. A better indicator is the quantity (rn-m)/(rf-g), which is the solution for net counts S in the set of simulatnaeous linear equations given above. This issue is addressed in CIAO 4.3.

In the meantime, we recommend the following work-around:

- Compute the value of n+m for which S=(rn-m)/(rf-g) ~30. Experience indicates that this is a reasonable threshold above which Gaussian statistics are appropriate
- Set the max_counts parameter to this value. For cases in which r>>1 (background aperture >> source aperture) and in which the source aperture includes most of the source (f->1, g->0), max_counts=n+m ~ m+30.

For other issues and bugs, see the bugs page for this tool