/ ciao4.16 / ahelp / aprates.html

Jump to: Description · Examples · Parameters · Bugs · See Also

AHELP for CIAO 4.16

aprates

Synopsis

Compute net counts, rate, and/or flux, plus credible intervals, 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] [verbose]

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 credible intervals, 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 credible interval are determined by summing values of the PDF alternately above and below the mode until the desired confidence level is attained. (Since this is a fully Baysian computation we use the term "credible interval" instead of the "confidence interval" used in frequentist statistics.) If the summation below the mode reaches the PDF for S=0, the lower bound is set to 0 and the summation continues for PDF values above the mode only. If the mode itself it 0, lower bound is set to INDEF and the summation proceeds for PDF values above the mode only.

Further documentation

For more detail on the algorithms, please see the documents https://cxc.harvard.edu/csc/memos/files/Kashyap_xraysrc.pdf, https://cxc.harvard.edu/csc/memos/files/Primini_significance.pdf, and https://cxc.harvard.edu/csc/memos/files/Primini_apflux_spec_II.pdf. The technique is discussed in Determining X-Ray Source Intensity and Confidence Bounds in Crowded Fields, Primini & Kashyap 2014, ApJ, 796, 24.

Upper limits vs. upper value for flux uncertainty

If the lower bound for the uncertainty of the flux reaches 0, only the upper bound is reported. This number is sometimes erroneously used as the upper flux limit for an undetected source, but strictly speaking the uncertainty on a flux measurement and the upper limit of an undetected source are two different questions, even though the numbers often come out to be in a similar range. When talking about the upper limit for an undetected sources (the maximal flux that a source could have without being detected), a single number for conf is not sufficient. Instead, two numbers are needed: The probability of detecting a source in background fluctuations (false positive) and the probability of missing a real source because it is compatible with background (false negative). This calculation is beyond scope of this tool.

User inputs

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

Examples

Example 1

Compute source counts n, background counts m, source aperture area A_s
and background aperture area A_b:

unix% dmextract \
infile="regevt3.fits[energy=500:7000][bin sky=region(reg3.fits)]" \
bkg="regevt3.fits[energy=500:7000][bin sky=region(reg3.fits[bkgreg])]" \
outfile=- | \
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 (source PSF fraction):

unix% dmstat psf3.fits"[sky=region(reg3.fits)]" \
centroid=no verbose=0
unix% pget dmstat out_sum
0.953

Compute beta (PSF fraction in background):

unix% dmstat psf3.fits"[sky=region(reg3.fits[bkgreg])]" \
centroid=no verbose=0
unix% pget dmstat out_sum
0.029

Run aprates to compute net counts and 68% credible intervals:

unix% 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
unix% 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:

unix% dmkeypar regevt3.fits LIVETIME echo+
5977.74

Run aprates to compute net rate (counts/second):

unix% 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
unix% 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 if the background was extracted from a different dataset.

Example 3

Compute E_s (mean exposure in source region):

unix% dmstat \
regexp3.fits"[sky=region(reg3.fits)]" \
centroid=no verbose=0
unix% pget dmstat out_mean
2308338.2013

Compute E_b (mean exposure in background region):

unix% dmstat \
regexp3.fits"[sky=region(reg3.fits[bkgreg])]" \
centroid=no verbose=0
unix% pget dmstat out_mean
2058975.0111

Run aprates to compute net photon flux:

unix% 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
unix% 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 (mean energy in source region):

unix% dmtcalc \
"regevt3.fits[energy=500:7000,sky=region(reg3.fits)]" \
- expression="energy=1.6e-12*energy" | dmstat "-[cols energy]" verbose=0
unix% pget dmstat out_mean
2.56e-09

Compute eng_b:
unix% dmtcalc \
"regevt3.fits[energy=500:7000,sky=region(reg3.fits[bkgreg])]" \
- expression="energy=1.6e-12*energy" | dmstat "-[cols energy]" verbose=0
unix% pget dmstat out_mean
4.27e-09

Run aprates to compute net energy flux:

unix% 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
unix% 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 (average flux in souce region):

unix% eff2evt \
"regevt3.fits[energy=500:7000][sky=region(reg3.fits)]" - | \
dmstat -"[cols flux]" verbose=0
unix% pget dmstat out_mean
9.81e-16

Compute flux_b (average flux in background region):

unix% eff2evt \
"regevt3.fits[energy=500:7000][sky=region(reg3.fits[bkgreg])]" - | \
dmstat -"[cols flux]" verbose=0
unix% pget dmstat out_mean
2.23e-15

Run aprates to compute net energy flux:

unix% 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
unix% 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
pdf	string		alternate
max_counts	integer		50
itermax	real		100
nsigma	real		5.0
pmin	real		1.0e-5
verbose	integer		0

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 background aperture. May be set to 1 if computation of background 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 credible 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 before switching to Gaussian statistics in evaluating the PDF.

The quantity (rn-m)/(rf-g), which is the solution for net counts S in the set of simultanaeous linear equations given above, is used to determine the transition between Bayesian and Gaussian approaches

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.

Comparison of Gaussian and Bayesian results indicates that a value of max_counts=50 (the default) is a good choice for a wide range of input parameters.

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.

Parameter=verbose (integer default=0)

Tool chatter and debugging output.

Beyond the normal tool chatter, at verbose levels above 0, aprates prints the probability distribution curves it computes the terminal (stdout). There is a separate table for each of the quanties it computes: counts, rates, photon fluxes, and energy fluxes.

Bugs

Caveats

# aprates (CIAO): WARNING: Large number of counts, just using Gaussian pdf

For low counts, aprates uses Poisson statistics, but for moderately large n or m (100's), this takes a very long time. If pdf=alternate, the tool uses the max_counts parameter to determine when to take an alternate path and use Gaussian statistics.

The quantity (rn-m)/(rf-g), which is the solution for net counts S in the set of simultanaeous linear equations given in the help file, is used to determine the transition between Bayesian and Gaussian approaches for the max_counts parameter.