Synopsis
Compute upper limits for the detection of a source, using data obtained from background apertures in event lists, images, and exposure maps.
Syntax
aplimits prob_false_detection prob_missed_detection outfile T_s A_s [bkg_rate] [m] [T_b] [A_b] [max_counts] [maxfev] [verbose] [clobber]
Description
The aplimits tool estimates the upper limit for a detection. This is a property of the detection process and depends on the values chosen for two error conditions: First, the maximum probability of a false detection (false positive), i.e. that a background fluctuation alone exceeds the calculated upper limit, and second the probability that a source with a flux of exactly the upper limit will be missed (false negative), because Poisson fluctuation of background and source counts lead to an observed count rate that is below the detection limit. The value of the upper limit depends only on the background rate and the probabilities chosen for the two errors, not on the observed number of counts in the source region. All input is obtained from the input parameter file; output is written to an output parameter file.
It is possible to find combinations of parameters that numerically result in a negative upper limit. Since this does not make sense physically, 0 is reported as the upper limit in these cases.
If the background rate is not known exactly, it can be estimated from the number of counts in a large source-free background region. In this case, the true background rate is not known, and a Bayesian computation needs to marginalize over the possible background rates. To perform this computation, one needs to set a prior on the background count rate. This tool uses an uninformative prior expressed as a Gamma distribution which represents an prior knowledge of 0 counts in 0 area. This uninformative prior is an "improper" prior, which means that it cannot be normalized. However, the gamma distribution is the conjugate prior for this process and thus the posterior for the background rate can be calculated analytically given the number of observed background counts. The functional form for the posterior is the same as for aprates.
This tool is based on the following article: On computing upper limits to source intensities (Kashyap et al. 2010, ApJ, 719, 900).
User inputs
This tool can either assume perfect knowledge of the background; in this case the parameter bkg_rate must be set. Alternatively, the user can set the number of counts m in the background region, and, unless source and background region have the same area and exposure time, also set the the geometric areas (A_s, A_b) and respective exposure times (T_s, T_b). If both bkg_rate and m are given, the tool assumes that the background rate is known and ignores the setting for m.
In practice, one never has perfect background knowledge, but in many cases the background rate can be determined well enough that the uncertainty does not matter in practice. For an isolated point source, the background can be measured from a large source-free region on the same chip. If the background region is large enough to yield a large number of counts (50 or more), then one can typically use background rate (number of counts in the background region divided by the area and exposure time of the background region) as known background rate. For extended sources, crowded fields, or exposures in sub-array mode, it might not be possible to find a large enough background region. In that case, the user should find the best possible uncontaminated background region and use the number of counts in that region as input for m, and also set the the geometric areas (A_s, A_b) and respective exposure times (T_s, T_b).
Return values
Output is written to an output parameter file. The parameter files contains two values: min_counts_detect, which is the minimum number of counts in the source region that is needed to claim a detection, and upper_limit, which is the upper limit on the count rate (counts per time) for a non-detection. Note that the upper limit is given for the entire source region and not normalized to the source area. In a typical case of a point source, the source aperture should be large enough to include most of the PSF and then the upper limit can be interpreted as the upper limit on the source flux. For an extended source, the upper limit applies to the part of the extended source that is covered by the source aperture.
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.
This tool itself does not account for the fraction of the PSF that falls outside the source aperture, but one can divide the resulting upper flux limit by the PSF fraction in the aperture to estimate the upper limit on the count rate.
2.) Upper limit on a detection vs. credible interval for the flux
This tool calculates the upper limit on the source flux, given properties of the detection method (type I and II errors and source region size and exposure time) and background level. This calculation is valid for all source regions in an observation (assuming a smooth background and constant exposure time). This is conceptually different from aprates, which instead looks at the source count rate in one specific region and then derives the credible interval for the flux of a source at one specific region. For undetected or weak sources, the credible interval includes 0 and sometimes the upper end of the interval is taken as upper limit on the source flux. That is not mathematically correct, since "upper limit" and "upper end of credible interval" describe different statistical concepts, even though the derived numbers can be similar; see Kashyap et al. (2010) for a more detailed discussion.
3.) What are the units for area and exposure time?
To calculate the upper limit, the tool needs to know the number of expected background counts in the source area. Because of that the formulas do not depend directly on the exposure time or aperture area. Instead, the source aperture area is only used to get the number of expected background counts (if bkg_rate is given) or the ratio to the background area (if the number of counts m in the background aperture is given). So, it is sufficient to use consistent units for areas, e.g., A_s and A_b can be given in pixels or arcsec^2, so long as both numbers use the same units (and if source and background aperture are the same size, they can both be left at their default 1.0). Similarly, bkg_rate can be "per pixel" or "per arcsec^2" as long as the unit is consistent with A_s. The same is true for the units of exposure time: T_s and bkg_rate (or T_b) just need to be consistent and the output will be in the same time unit. If T_s is given in seconds, then the upper limit will be returned as counts/sec; but one can also use "exposure time" as the unit: if T_s and T_b are left at the default of 1.0 (or bkg_rate is "counts over the entire observation"), then the upper limit will be in "counts over the entire observation".
4.) Incomplete PSF in source aperture
The calculation of the upper limit is done in counts. If the source aperture does not cover the full PSF, it will only apply to that fraction of the PSF, e.g., if A_s only include 90% of the PSF, then the resulting upper limit is a limit on 90% of the source flux; in other words, the limit on the full source flux is upper_limit/0.9.
Examples
Example 1
unix% aplimits prob_false_detection=.05 prob_missed_detection=.5 bkg_rate=1.23 outfile=aplimits_out.par unix% pget aplimits_out.par upper_limit 2.44206 unix% pget aplimits_out.par min_counts_detect 3
In this example, we know the background rate and calculate the upper limit on the source count rate (in cts/s). We want the probability for a background fluctuation to appear as a source to be less than 5%, and we want the probability to miss a real source with a flux comparable to the upper limit to be below 50%. We do not set the source extraction area and exposure time. Both default to 1. So, we can interpret the resulting number as source flux per unit time. However, the calculation of the limits does not directly depend on the exposure time; it is really just about the number of counts. It also does not depend directly on A_s, that is only needed to calculate the number of expected background counts. The defaults for exposure time and area are 1.0 for each value. So, we can also view this example as follows: We expect 1.23 background counts in the source area over the entire exposure time and and the upper limit to the source flux is 2.44206 counts over the entire exposure time. To convert the upper limit on the source flux into a rate of cts/s, we can divide 2.44206 by the exposure time in seconds.
Example 2
unix% aplimits prob_false_detection=.1 prob_missed_detection=.5 bkg_rate=INDEF m=5 A_b=4.065 outfile=aplimits_out.par unix% pget aplimits_out.par upper_limit 2.2316
If the background rate is not known exactly (bkg_rate set to INDEF), the tool can marginalize over all possible background rates given the information on a single background measurement. With 5 counts in the background aperture and a region that is about 4.065 times larger than the source aperture, one might naively expect a background count rate of 5/4.065=1.23. However, given Poisson statistics, the true background rate could be higher or lower and thus the resulting upper limit is about 10% smaller than in the first example. When the background is measured better (more counts in a larger background area) the difference between the calculation for a perfectly known background and marginalizing over the background rate is even smaller.
Example 3
unix% aplimits prob_false_detection=.05 prob_missed_detection=.5 bkg_rate=0.34 A_b=1 T_b=1 outfile=aplimits_out.par unix% pget aplimits_out.par upper_limit 1.338348388671875 unix% aplimits prob_false_detection=.05 prob_missed_detection=.5 bkg_rate=0.34 A_b=10 T_b=23 clobber+ outfile=aplimits_out.par unix% pget aplimits_out.par upper_limit 1.338348388671875
The background rate is given per unit are and unit time, so changing the background area and exposure time has no effect on the results when a background rate is input.
Parameters
| name | type | def | min | max | reqd |
|---|---|---|---|---|---|
| prob_false_detection | real | 0.1 | 0 | 1 | yes |
| prob_missed_detection | real | 0.5 | 0 | 1 | yes |
| outfile | file | yes | |||
| T_s | real | 1 | 0 | yes | |
| A_s | real | 1 | 0 | yes | |
| bkg_rate | real | INDEF | 0 | ||
| m | integer | 0 | |||
| T_b | real | 1 | 0 | ||
| A_b | real | 1 | 0 | ||
| max_counts | integer | 50 | 0 | ||
| maxfev | integer | 500 | 1 | ||
| verbose | integer | 1 | 0 | 5 | |
| clobber | boolean | no |
Detailed Parameter Descriptions
Parameter=prob_false_detection (real required default=0.1 min=0 max=1)
Upper limit for the probability of a false detection
The probability of a type I error describes the possibility a background fluctuation alone causes an observed count rate above the detection threshold.
Parameter=prob_missed_detection (real required default=0.5 min=0 max=1)
Probability of detecting a source with a flux equal to the upper limit
The probability of a type II error describes the possibility that a true source with a flux of the upper limit will be missed in the detection because the Poisson fluctuations of the source and background flux lead to a count number below the detection threshold.
Parameter=outfile (file required)
Filename of output file
Parameter=T_s (real required default=1 min=0)
Exposure time in source aperture
This is given in the same units as exposure time of the background T_b.
Parameter=A_s (real required default=1 min=0)
Geometric area of source aperture
This is given in the same units as geometric area of the background A_b.
Parameter=bkg_rate (real default=INDEF min=0)
Background count rate
Known background rate. The flux is given per area and time. The unit of the area must match the units used for A_s. When the exact background rate is not known, bkg_rate has to be set to INDEF. In that case, the number of observed background counts m has to be provided, so that the tools can marginalize over the background rate.
Parameter=m (integer min=0)
Number of counts in background aperture
Parameter=T_b (real default=1 min=0)
Exposure time in background aperture
This is given in the same units as exposure time of the source T_s. This parameter is only used if bkg_rate is not given.
Parameter=A_b (real default=1 min=0)
Geometric area of background aperture
This is given in the same units as geometric area of the source A_s. This parameter is only used if bkg_rate is not given.
Parameter=max_counts (integer default=50 min=0)
Background count number above which the uncertainty on the background is ignored
When the background count rate is not know exactly, the tool performs a Bayesian computation and marginalizes over the background. In particular for a large number of counts, this computation can take a long time or be numerically unstable. On the other hand, when the number of counts in the background is large, that means that the uncertainty on the background flux is small and there is little difference between the full Bayesian computation and assuming a known background rate. When the number of background counts m becomes larger than max_counts, then the tool switches from the Bayesian computation to the a faster and more numerically stable formula for an exactly known background rate of m divided by the area and exposure time of the background exposure.
Parameter=maxfev (integer default=500 min=1)
Maximal number of function evaluations in numerical root finding
Parameter=verbose (integer default=1 min=0 max=5)
Verbose level
If set to a non-zero value then the tool will print information to the screen when it is run. The extra information produced when verbose is greater than 1 is only likely to be useful when debugging the script.
Parameter=clobber (boolean default=no)
OK to overwrite existing output file?
Changes in the scripts 4.15.1 (January 2023) release
The number of steps used in the numerical integration has been increased to improve the behavior in certain cases.
About Contributed Software
This script is not an official part of the CIAO release but is made available as "contributed" software via the CIAO scripts page . Please see this page for installation instructions - such as how to ensure that the parameter file is available.
Bugs
For an up-to-date listing of known bugs, see the bugs page for this tool.
Refer to the CIAO bug pages for an up-to-date listing of known issues.