Last modified: 10 May 2023

URL: https://cxc.cfa.harvard.edu/ciao/threads/upperlimit/

Computing the Intensity Upper Limit and Upper Value for Flux Uncertainty for an Unresolved Source

CIAO 4.16 Science Threads


Overview

Synopsis:

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.

Purpose:

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.


Contents


Get Started

Download the sample data: 315 (ACIS-S, NGC 4038/39)

Run srcflux

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

[Thumbnail image: The source of interest is marked as a white circle in both ds9 frames.]

[Version: full-size]

[Print media version: The source of interest is marked as a white circle in both ds9 frames.]

Figure 1: Apertures overlaid on x-ray and optical data

The left frame contains the Chandra x-ray data and the right frame contains an optical image from the DSS. The location (ra,dec)=(12:01:47.9, -18:51:15.6) is circled in both frames.

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:

circle(12:01:47.9,-18:51:15.6,5")

annulus(12:01:47.9,-18:51:15.6,5",15")
[NOTE]
Note

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.


Step-by-step: aplimits

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"

History

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.