Determining CSC 2.0 Cumulative Sky Coverage in a Region of the Sky¶

Science thread by Frank Primini. Jupyter notebook by Rafael Martínez-Galarza¶

This thread describes how to construct a cumulative sky coverage histogram for a region of the sky containing one or more overlapping CSC 2.0 stacks. This notebook assumes you have installed CIAO 4.13 using conda, but you should also be able to do this with a CIAO installed using the ciao-install script, or any other Python environment as we are just going to use Astropy along with pyvo. The pyyaml poackage is needed to support writing out the file (it avoids a warning message you would get).

Packages to install:

• pip install astropy pyvo pyyaml matplotlib

In addition, you will need to have the following installations:

Scientific Problem¶

The number counts or luminosity functions of different classes of x-ray sources allow the determination of those sources' contributions in various astrophysical problems, ranging in scale from specific galaxy populations to the overall contribution the the cosmic x-ray background. To determine the number counts for sources in flux-limited surveys, one needs to know the cumulative sky coverage $\Omega(F>F_{min})$, i.e., the survey area in which a source brighter than flux, $F_{min}$, could be detected, as a function of $F_{min}$ (see, e.g. Miyaji et al. (2015)).

Sky coverage can be determined from a map of limiting sensitivity, defined as the minimum flux required to include a source in the survey. Although CSC 2.0 includes an all-sky limiting sensitivity map which users can query at specific locations using CSCview, at present there is no capability to download subsets of the map for arbitrary regions of the sky. However, the equivalent information may be recovered from the limiting sensitivity data products provided for individual CSC 2.0 stacks. In this thread we use pyvo and the CIAO script search_csc to determine and download the limiting sensitivity data products for a particular region of the sky, and analyze them with CIAO to construct the cumulative sky coverage for the region.

Determining stacks in a region of the sky¶

For our example, we choose a $12^{\prime}$ radius circular region centered at $\alpha=164.05, \delta=57.2$ in the southeast corner of the larger Chandra Lockman Hole Survey. Here is a visualization of the region in the Chandra Footprint Service:

The Service identifies six ObsIDs: 1682, 1686, 1687, 1691, 1692, and 1693, that overlap the region. Let's start by defining the region that we want to query:

We now query the Catalog using our cone search service and pyvo to determine the sources that fall within our region of interest. We use the coordinate object that we have just defined and then print the names of the CSC 2.0 sources that fall within the region of interest, and save the names in a list (tnames) for later use.

We now want to find the association between these sources and the specific Chandra observations and stacks that contain them (remember that in CSC 2.0 we stack individual observations to improve detection). In order to make these associations, we use the the catalog's Table Access Protocol (TAP) service, which we query using pyvo. We retrieve the source name, coordinates, and the ObsID, stack and region number that identify each detection of each of our sources (for more on how to do TAP queries see this tutorial):

Let us now list the results in a table format. For each source, we list each individual observation where the source was detected, and the corresponding stack ID:

Downloading the data products: sensitivity maps¶

In order to estimate the cumulative sky coverage, we need to use the stack level sensitivity maps, which provides the limiting sensitivity for a point source to exceed the MLE likelihood thresholds at each location in the field of each stack, i.e., the the minimum photon flux that a source can have in order to be detected as a TRUE or MARGINAL in that particular stack.

We wish to obtain the 's' band limiting sensitivity maps corresponding to the detect_stack_ids that we have identified. Here we will be using the soft band sensitivity maps (multiple energy bands can be retrieved at the same time, but for our example we're only interested in the 's' band).

In order to download the data products, we use CIAO's search_csc tool, which takes the coordinates of a source as input, or a SIMBAD-resolved name. For each unique stack name, we select the first listed source from the table above to perform the query. This operation will download the sensitivity maps for each of the stacks that overlap our region of interest. Some of them will be downloaded more than once because a given source can be detected in more than one stack, and we are using the coordinates of a source in order to retrieve the data products. We therefore need to get the unique files after the download, as shown below:

Analyzing Limiting Sensitivity Map Data¶

The organization of the sensitivity maps is described in our documentation pages. For our example, we'll look at the sensitivities corresponding to detections with likelihood classification TRUE, found in the second fits extension of the sens3.fits file.

Visualizing the Maps¶

We can examine the sensitivity maps to get an idea of where they fall with respect to our circular region of interest. We do this below using Astropy's WCS package, which allows us to show the sensitivity maps in the coordinate reference frame defined in the header of the files:

Creating a sensitivity map mosaic¶

Our circular region is located in different image coordinates in each map. To combine the maps properly we need to reproject them to a common coordinate system. We do this below using CIAO's reproject_image_grid tool. We use reproject_image_grid to reproject each sensitivity map to a common WCS, centered at $\alpha=164.05, \delta=57.2$. We use a blocking factor of $8$ to match the blocking of the sensitivity maps. The pixelsize parameter for reproject_image_grid is then

\begin{equation*} pixelsize = \frac{0.492^{{\prime}{\prime}} \times 8}{3600.0^{{\prime}{\prime}}/\mathbb{deg}} = 1.09333\times10^{-3} \mathbb{deg}. \end{equation*}

To ensure that our mosaic encompasses the entire $12^{\prime}$ radius circle, we require the reproject_image_grid xsize and ysize parameters to be

\begin{equation*} xsize = ysize = \frac{24.0^{{\prime}} \bigg/ 60.0^{{\prime}}/\mathbb{deg}}{1.09333\times10^{-3} \mathbb{deg}} = 365.85 \end{equation*} which we round up to $366$ image pixels. We can then reproject each sensitivity map:

Next we display the reprojected maps where we match all frames using their WCS values. The result is shown below. All reprojected maps are properly centered about our circular region.

Note that the pixels outside the chip boundaries appear to have two different values. Those in the original image space (gray) appear as INDEF; however, those in the reprojected image outside the original image space (black) are set to 0. Before combining the reprojected maps we need to set these zero-values to INDEF as well. We use the CIAO tool dmimgthresh to set all values below a very small positive number to INDEF.

We can now combine the individual maps into a mosaic using the CIAO tool, dmimgfilt. We use this tool so that we can follow the prescription that when two or more maps contribute to the same image pixel, the map with the lowest pixel value should be used. Before we do this, we create a list of the reprojected maps:

And now we run the mosaic maker, and visualize it:

Finally, we apply a spatial filter to the mosaic to null out values outside our circle, using CIAO's dmcopy, and visualize the filtered version:

Calibrating the mosaic¶

As discussed in the "Limiting Sensitivity and Sky Coverage" section of the Catalog Statistical Characterization web page, the sensitivity maps are in units of $p_{min}$, the photon flux that corresponds to the minimum Poisson background fluctuation in the source aperture that would exceed the source detection likelihood threshold. To convert to the equivalent photon or energy flux derived from the likelihood of the PSF fit, which is actually used in the detection process, we need to apply a calibration. For energy flux in the s band, that calibration relation is:

\begin{equation*} \log_{10}(Flux) = 1.028 \times \log_{10}(p_{min}) - 8.595, \end{equation*}

which we apply using dmimgcalc.

We choose to keep the $\log_{10}$ of the sensitivity, since we want our ultimate sky coverage histogram to be in logarithmic bins.

Computing a Sky Coverage Histogram from the Mosaic¶

If we are only interested in a spatial map of the sensitivity in our region, we can stop here. However, if we wish to compute number counts or luminosity functions, we need to construct the cumulative sky coverage histogram. We begin by examining logFlux.fits with dmstat to determine a suitable range for the histogram.

We choose to construct the histogram in 0.05 decade bins from -14.85 to -13.65 $\left(\sim1.4\times10^{-15}-2.2\times10^{-14}~ \mathbb{erg \cdot cm^{-2} \cdot s^{-1}}\right)$, using dmhist

The COUNTS column is the number of $1.09333\times10^{-3} \times 1.09333\times10^{-3}~\mathbb{deg}^{2}$ pixels in each histogram bin. To compute the cumulative sky coverage, we need to compute $COUNTS\times1.09333\times10^{-3}\times1.09333\times10^{-3}$ for each bin, and compute the running sum. We construct the cumulative histogram of sky coverage by summing the contribution from each pixel

Finally, we plot the histogram, which represents the survey area in which a source brighter than flux $F_{min}$ could be detected, as a function of $F_{min}$