Computing the Intensity Upper Limit and Upper Value for Flux Uncertainty for an Unresolved Source
CIAO 4.15 Science Threads
The tool aprates, 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, as shown in the Compute Net Counts, Rate, or Flux for Point Sources thread. aprates uses Bayesian statistics to compute the background-marginalized, posterior probability distribution for source intensity, assuming non-informative prior distributions for background and source intensity (for details on the algorithm, see the Background Marginalized X-ray Source Intensity memo). The posterior distribution can be used to determine intensity value and confidence bounds or intensity upper limit.
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). The aplimits tool performs this calculation.
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. This tool is based on the following article: On computing upper limits to source intensities (Kashyap et al. 2010, ApJ, 719, 900).
The aprates steps in this thread have been automated in the srcflux script. This thread shows how to get the upper value on the flux uncertainty this way.
The thread then shows how to run aprates and aplimits step-by-step to determine how x-ray bright a source could be at the given location in the observation. The intensity upper limit is calculated in photon units in this thread.
To calculate the intensity upper limit for a source within a specified aperture.
Last Update: 10 May 2023 - Updated to make the distinction between true upper limits, now available using the new aplimits tool, and the upper value on the flux uncertainty computed by srcflux and aprates.
- Get Started
- Run srcflux
- Step By Step: aprates
- Step-by-step: aplimits
- Parameter files:
Download the sample data: 315 (ACIS-S, NGC 4038/39)
unix% download_chandra_obsid 315
Figure 1 shows the x-ray data (left) and an optical image from the Digital Sky Survey (right). The optical image was retrieved from SAO-DSS via the ds9 "Analysis → Image Servers" menu. The star circled in white - 2MASS J12014790-1851156 - is clearly visible in the optical but not in x-ray.
Figure 1: Apertures overlaid on x-ray and optical data
The location (ra,dec)=(12:01:47.9, -18:51:15.6) will be used as the center of the apertures in this thread.
With srcflux we can get the upper limit on in the Chandra broad-band (0.5 to 7.0 keV) by simply
unix% punlearn srcflux unix% pset srcflux infile=acisf00315_repro_evt2.fits unix% pset srcflux outroot=dss unix% pset srcflux pos="12:01:47.9 -18:51:15.6" unix% pset srcflux psfmethod=quick unix% srcflux ... Summary of source fluxes Position 0.5 - 7.0 keV Value 90% Conf Interval #0001|12 1 47.89 -18 51 15.6 Rate 1.3E-06 c/s (0,5.29E-05) Flux NAN erg/cm2/s (NAN,NAN) Mod.Flux 6.49E-18 erg/cm2/s (0,2.64E-16) Unabs Mod.Flux 6.93E-18 erg/cm2/s (0,2.81E-16)
srcflux automatically creates source and background regions and extracts all the quantities discussed in the Step by Step: aprates section (below). The output shows that the 90% upper limit at the source location in the 0.5 to 7.0 keV range is roughly 5.29E-5 counts per second.
The tool outputs several files including the dss_broad.flux file that contains a summary of all the counts, rates, areas, and flux values.
unix% dmlist "dss_broad.flux[cols counts,bg_counts,area,bg_area,net_rate_aper_hi]" data,clean # COUNTS BG_COUNTS AREA BG_AREA NET_RATE_APER_HI 1.0 23.0 20.750 518.50 5.28062E-05
and also the upper limit on the flux uncertainty in units of photon/cm^2/sec:
% dmlist "dss_broad.flux[cols net_photflux_aper,net_photflux_aper_lo,net_photflux_aper_hi]" data,clean # NET_PHOTFLUX_APER NET_PHOTFLUX_APER_LO NET_PHOTFLUX_APER_HI 3.752809674E-09 0 1.524576588E-07
The next section shows how to run the steps necessary to run the aprates tool.
Step By Step: aprates
Selecting the Apertures
A circle with radius 5" is used as the source region and an annulus with inner radius of 5" and outer radius of 10" is used for the background:
These are not the same as the apertures used above so the results are expected to be slightly different.
Counts in source and background apertures
The number of counts in the source (n) and background (m) apertures is found by running dmlist on the event file with the source and background regions. The energy range of the events is restricted to 0.5-7.0 keV:
unix% dmlist 'acisf00315_repro_evt2.fits[energy=500:7000][(ra,dec)=circle(12:01:47.9,-18:51:15.6,5")]' counts 15 unix% dmlist 'acisf00315_repro_evt2.fits[energy=500:7000][(ra,dec)=annulus(12:01:47.9,-18:51:15.6,5",15")]' counts 132
There are 15 source counts (n=15) and 132 background counts (m=132).
Areas of source and background apertures
The area of the source (A_s) and background (A_b) apertures is computed:
Source: 5^2 * pi = 78.5 Background: (15^2 - 5^2) * pi = 628.3
The source area is 78.5 square arcsec (A_s=78.5) and the background area is 628.3 square arcsec (A_b=628.3).
Exposure time of the observation
By setting the exposure times in the source (T_s) and background (T_b) apertures, the net source rate and errors will be computed.
For this example, the exposure time is the same for the source and background since they come from the same file:
unix% dmkeypar acisf00315_repro_evt2.fits LIVETIME echo+ 72242.2
The exposure time is 72242.2 s.
Mean Exposure values
The mean exposure map value (in units cm**2*sec) is required to get the upper value on the flux uncertainty. Since srcflux created the exposure map already we can use that. Alternatively users can run the fluximage script to generate an exposure map.
% dmstat dss_0001_broad_thresh.expmap'[(ra,dec)=circle(12:01:47.9,-18:51:15.6,5")]' cen- sig- med- EXPMAP[cm**2 s] min: 25943108 @: ( 4353 4007 ) max: 26088824 @: ( 4364 3990 ) mean: 26021691.451 sum: 8535114796 good: 328 null: 113 % dmstat dss_0001_broad_thresh.expmap'[(ra,dec)=annulus(12:01:47.9,-18:51:15.6,5",15")]' cen- sig- med- EXPMAP[cm**2 s] min: 25867038 @: ( 4352 4013 ) max: 26158504 @: ( 4373 3985 ) mean: 26014718.441 sum: 13345550560 good: 513 null: 328
Run aprates to find the limit
For this calculation, it is assumed that 99% of the PSF is contained within the source region (alpha=.99) and 1% is in the background region (beta=.01). The following parameters are all set to "1" as we are not calculating flux in energy units: eng_s, flux_s, eng_b, and flux_b.
The aprates parameters are set with the values determined above and the tool is run. For details on how the tool does the calculations, refer to the aprates help file.
unix% punlearn aprates unix% pset aprates n=15 m=132 unix% pset aprates A_s=78.5 A_b=628.3 unix% pset aprates alpha=0.99 beta=0.01 unix% pset aprates T_s=72242.2 T_b=72242.2 unix% pset aprates E_s=26021691.451 E_b=26014718.441 unix% pset aprates eng_s=1 flux_s=1 eng_b=1 flux_b=1 unix% pset aprates outfile=ulimit.par conf=0.68 unix% aprates mode=h
You may see a warning:
# aprates (CIAO 4.4): WARNING: Large number of counts, just using Gaussian pdf
It indicates that Gaussian statistics were used instead of Poisson statistics; the threshold is controlled by the max_counts parameter.
The output file is in parameter file format, which makes it possible to query the output values with the tool pget:
unix% pget ulimit.par src_rate src_rate_err_lo src_rate_err_up 0 INDEF 5.775e-05
At the 68% confidence level, the upper value on the x-ray count rate uncertainty is 5.775e-05 counts/s. The photon flux count rate upper value is
% pget ulimit.par photflux_aper photflux_aper_err_lo photflux_aper_err_up 0 INDEF 1.58745e-07
which is in units of photon/cm^2/sec.
The full output file, ulimit.par, is included at the end of the thread.
Where as aprates computed the upper value on the flux uncertainty, aplimits will compute an upper limit on detecting the source. Unlike aprates, aplimits only relies on the number of counts in the background region together with the exposure times and the areas of the regions. All of these values have already been computed in the previous section.
unix% punlearn aplimits unix% pset aplimits m=132 unix% pset aplimits bkg_rate=INDEF unix% pset aplimits A_s=78.5 A_b=628.3 unix% pset aplimits T_s=72242.2 T_b=72242.2 unix% pset aplimits outfile=ulimit_2.par clobber=yes unix% aplimits prob_missed_detect=0.5 prob_false_detect=0.1 mode=h aplimits prob_false_detection = 0.1 prob_missed_detection = 0.5 outfile = ulimit_2.par T_s = 72242.19960592801 A_s = 78.5398 bkg_rate = INDEF m = 132 T_b = 72242.19960592801 A_b = 628.3184 max_counts = 50 maxfev = 500 verbose = 1 clobber = yes mode = h Minimum number of counts for detection: 22 counts Upper limit: 8.537294343113899e-05
The 10% probability upper limit on the count rate for missing a source is 8.5e-5 counts/sec, and for a source to be missed with 50% probability it would have to have at least 22 counts.
The values are also stored in the output parameter file, ulimit_2.par
unix% plist ulimit_2.par Parameters for ulimit_2.par (min_counts_detect = 22) Detection threshold (upper_limit = 8.537294343113899e-05) Upper limit (mode = ql)
To convert these values into physical units we can use the modelflux tool. It requires and ARF and RMF which can be created using specextract, but since srcflux has already done that we will use those file here:
unix% punlearn modelflux unix% pset modelflux arf=dss_0001.arf rmf=dss_0001.rmf unix% pset modelflux model="xsphabs.abs1" absmodel="xspowerlaw.pow1" unix% pset modelflux paramvals="pow1.PhoIndex=2.0" absparams="abs1.nH=0.0394" unix% modelflux emin=0.5 emax=7.0 rate=8.537e-5 opt=rate mode=h Model fluxes: Rate (0.5,7)= 8.537e-05 count s^-1 Photon Flux (0.5,7)= 1.9862e-07 photon cm^-2 s^-1 Energy Flux (0.5,7)= 4.8176e-16 erg cm^-2 s^-1 Unabsorbed model fluxes: Rate (0.5,7)= 0.024192 count s^-1 Photon Flux (0.5,7)= 7.7254e-05 photon cm^-2 s^-1 Energy Flux (0.5,7)= 4.7141e-13 erg cm^-2 s^-1
The nH value was obtained from header of the spectrum, .pi file
unix% dmlist dss_0001.pi keys,clean | grep -i NH NRAO_NH 0.03940 [10**22 cm**-2] galactic HI column density BELL_NH 0.03850 [10**22 cm**-2] galactic HI column density
modeflux stores the output values back into its own parameter file:
unix% pdump modelflux arf='dss_0001.arf' rmf='dss_0001.rmf' model='xsphabs.abs1' paramvals='pow1.PhoIndex=2.0' emin='0.5' emax='7' absmodel='xspowerlaw.pow1' absparams='abs1.nH=0.0394' abund='angr' oemin='INDEF' oemax='INDEF' rate='8.537e-05' pflux='1.9862e-07' flux='4.8176e-16' urate='0.024192' upflux='7.7254e-05' uflux='4.7141e-13' opt='rate' verbose='1' mode='ql' # EOF
Parameters for /home/username/cxcds_param/ulimit.par src_cnts,r,h,0.0,,, src_cnts_err_lo,r,h,INDEF,,, src_cnts_err_up,r,h,4.17198,,, src_cnts_conf,r,h,0.684242,,, src_cnts_status,i,h,0,,, src_cnts_signif,r,h,0,,, src_cnts_mode,r,h,0,,, src_rate,r,h,0.0,,, src_rate_err_lo,r,h,INDEF,,, src_rate_err_up,r,h,5.775e-05,,, src_rate_conf,r,h,0.684242,,, src_rate_status,i,h,0,,, src_rate_signif,r,h,0,,, src_rate_mode,r,h,0,,, photflux_aper,r,h,0.0,,, photflux_aper_err_lo,r,h,INDEF,,, photflux_aper_err_up,r,h,1.58745e-07,,, photflux_aper_conf,r,h,0.680054,,, photflux_aper_status,i,h,0,,, photflux_aper_signif,r,h,0,,, photflux_aper_mode,r,h,0,,, flux_aper,r,h,0.0,,, flux_aper_err_lo,r,h,0,,, flux_aper_err_up,r,h,0,,, flux_aper_conf,r,h,0,,, flux_aper_status,i,h,4,,, flux_aper_signif,r,h,0,,, flux_aper_mode,r,h,0,,, eflux_aper,r,h,0.0,,, eflux_aper_err_lo,r,h,0,,, eflux_aper_err_up,r,h,0,,, eflux_aper_conf,r,h,0,,, eflux_aper_status,i,h,4,,, eflux_aper_signif,r,h,0,,, eflux_aper_mode,r,h,0,,, # ---- Inputs below ---- _n,i,h,15,,, _A_s,r,h,78.5398,,, _alpha,r,h,1,,, _T_s,r,h,72242.2,,, _E_s,r,h,2.60217e+07,,, _eng_s,r,h,0,,, _flux_s,r,h,0,,, _m,i,h,132,,, _A_b,r,h,628.318,,, _beta,r,h,0,,, _T_b,r,h,72242.2,,, _E_b,r,h,2.60147e+07,,, _eng_b,r,h,0,,, _flux_b,r,h,0,,,
Parameters for /home/username/cxcds_param/ulimit_2.par min_counts_detect,i,h,22,,,"Detection threshold" upper_limit,r,h,8.537294343113899e-05,,,"Upper limit"
|24 Jul 2009||New for CIAO 4.1: title updated to "Computing the Intensity Upper Limit for an Unresolved Source"|
|27 Jan 2010||updated for CIAO 4.2: use N004 version of data products from the Archive; screen output updated|
|15 Dec 2010||updated for CIAO 4.3: minor changes to screen output|
|10 Jan 2012||reviewed for CIAO 4.4: no changes|
|03 Dec 2012||Review for CIAO 4.5; no changes|
|11 Dec 2013||Review for CIAO 4.6; made step-by-step section and introduced a new section showing how to do in 1 command with srcflux|
|23 Dec 2014||Review for CIAO 4.7; minor edits only.|
|15 Feb 2022||Review for CIAO 4.14. Updated for Repro-5/CALDB 4.9.6.|
|10 May 2023||Updated to make the distinction between true upper limits, now available using the new aplimits tool, and the upper value on the flux uncertainty computed by srcflux and aprates.|