Last modified: 15 Feb 2022

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

Running wavdetect or celldetect on merged data: choosing psffile

CIAO 4.16 Science Threads


Overview

Synopsis:

The primary reason to combine (aka merge) observations is to detect faint sources. The most popular CIAO source detection tool is wavdetect, which works by correlating the input image with a series of Mexican hat wavelets. The optimum wavelet size or scale is that which matches the size of the sources being detected. For point sources, this scale is comparable to the size of the Point Spread Function (PSF). Users can provide wavdetect an input file with the size of the PSF via the psffile parameter. Similarly, the celldetect tool also can use the size of the PSF, input via the psffile parameter to adjust the detect cell size across the input images.

Unfortunately, the Chandra PSF varies significantly across the field of view, varying in size from a fraction of an arc-second near the aim-point, to several arc-seconds at the edge of the detectors; moreover, the observed orientation varies with the roll angle of the spacecraft. Therefore there is no single PSF size when arbitrary observations are combined with different aim-points and roll angles.

This thread demonstrates several alternate ways in which users can combine PSF maps from individual observations, and the effects the different methods have on the source detections. These results are demonstrated for a single target; users should not conclude that similar results will be obtained in all datasets.

Purpose:

To illustrate different ways to combine per-observation PSF map files and the effects on the detections.

Related Links:

Last Update: 15 Feb 2022 - Review for CIAO 4.14. Updated for Repro5/CALDB 4.9.6.


Contents


Getting Started

In this thread we will be combining three observations of NGC4485 and use the merged data as input to wavdetect. We begin by searching for the appropriate ObsIDs.

unix% find_chandra_obsid "ngc4485"
# obsid  sepn   inst grat   time    obsdate   piname              target
1579      3.1 ACIS-S NONE   19.5 2000-11-03   Murray "NGC 4485/NGC 4490"
4725      3.1 ACIS-S NONE   38.5 2004-07-29  Roberts "NGC 4485/NGC 4490"
4726      3.1 ACIS-S NONE   39.6 2004-11-20  Roberts "NGC 4485/NGC 4490"
20999     4.2 ACIS-S NONE   14.9 2018-11-06 Patnaude            SN2008ax
23482     2.4 ACIS-S NONE   28.7 2020-11-27 Earnshaw       "NGC 4485/90"
23483     2.4 ACIS-S NONE   29.6 2020-12-27 Earnshaw       "NGC 4485/90"
23484     2.4 ACIS-S NONE   29.6 2021-01-24 Earnshaw       "NGC 4485/90"

For simplicity we will only use the first three datasets. We can then download the files

Download the sample data: 1579; 4725; 4726 (NGC 4485/4490)

unix% download_chandra_obsid 1579,4725,4726 

or the find_chandra_obsid download parameter could have been set to all.

We then apply the latest calibrations with the chandra_repro script. To make things easier, we will store all the outputs from all 3 observations into a single directory.

unix% chandra_repro "*" ./repro verb=1
 ...
unix% /bin/ls
1579    4725    4726    repro

The value of * behaves just like the UNIX wild card and will process all the sub directories.



Running merge_obs

Now we use the merge_obs script to create the combined counts image, and associated exposure maps, and PSF maps.

In this thread we choose to set the binsize to '2' to make the examples have reasonable size. In practice users may want to use a binsize of 1 especially if the field is crowded and you are looking for sources near the aim point.

As of CIAO 4.11, the merge_obs script has been updated to create individual PSF maps for each observation and to merge those PSF maps together with various options. In this thread we will create PSF maps for each observations which enclose 90% of the PSF at the broad-band monochromatic energy (2.3keV), and then combine them together by weighting the individual maps by their exposure time.

unix% merge_obs "repro/*repro_evt2.fits" ngc4485 bands=broad bin=2 psfecf=0.9 psfmerge=exptime
Running merge_obs
Version: 05 November 2021

Verifying 3 observations.
Using CSC ACIS broad science energy band.
Calculating new tangent point.
New tangent point: RA=12h 30m 35.088s Dec=41d 38' 46.73"

Observations to be reprojected:

  Obsid  Obs Date   Exp    DETNAM     SIM_Z    FP   Sepn   PA  
                   (ks)                (mm)    (K)   (')  (deg)
---------------------------------------------------------------
1 1579  2000-11-03  19.5 ACIS-235678 -190.140 153.3   0.4   -28
2 4725  2004-07-29  38.5 ACIS-235678 -190.140 153.4   2.1  -158
3 4726  2004-11-20  39.6 ACIS-235678 -190.140 153.4   1.8   +31

WARNING - DATAMODE values differ:
  Obsid 1579 has DATAMODE=FAINT and the rest have VFAINT

Running tasks in parallel with 4 processors.
Reprojecting 3 event files to a common tangent point.
Merging reprojected events files to ngc4485_merged_evt.fits

Calculating the output grid

The merged images will have 2143 by 2426 pixels, a pixel size of 0.984 arcsec,
    and cover x=1876.5:6162.5:2, y=1498.5:6350.5:2.

Creating the fluxed images for 3 observations in parallel.
Creating 6 aspect histograms for obsid 1579
Creating 6 aspect histograms for obsid 4725
Creating 6 aspect histograms for obsid 4726
Creating 6 instrument maps for obsid 1579
Creating 6 instrument maps for obsid 4725
Creating 6 instrument maps for obsid 4726
Creating 6 exposure maps for obsid 1579
Creating 6 exposure maps for obsid 4725
Creating 6 exposure maps for obsid 4726
Combining 6 exposure maps for obsid 1579
Combining 6 exposure maps for obsid 4725
Thresholding data for obsid 1579
Thresholding data for obsid 4725
Exposure-correcting image for obsid 1579
Creating PSF map for obsid 1579
Exposure-correcting image for obsid 4725
Creating PSF map for obsid 4725
Combining 6 exposure maps for obsid 4726
Thresholding data for obsid 4726
Exposure-correcting image for obsid 4726
Creating PSF map for obsid 4726

Combining 3 observations.

The following files were created:

The co-added clipped counts image is:
     ngc4485_broad_thresh.img

The co-added clipped exposure map is:
     ngc4485_broad_thresh.expmap

The combined PSF map is:
     ngc4485_broad_thresh.psfmap

The combined FOV is:
     ngc4485_merged.fov

The co-added exposure-corrected image is:
     ngc4485_broad_flux.img

Warning: the merged event file ngc4485_merged_evt.fits
   should not be used to create ARF/RMF/exposure maps because
      the RA_NOM keyword varies by 0.0386 (limit is 0.0003)
      the DEC_NOM keyword varies by 0.0580 (limit is 0.0003)
      the ROLL_NOM keyword varies by 243.1 (limit is 1.0)

The counts image and exposure map are shown in Figure 1.

Figure 1: Merged Image and Exposure Map

[Thumbnail image: ]

[Version: full-size]

[Print media version: ]

Figure 1: Merged Image and Exposure Map

(Left) Image of the counts by combining the 3 observations of NGC4485 and (Right) Image of the combined exposure map. We display the counts image because that is what is input to wavdetect which uses Poisson statistics and therefore must be integer values.

The counts image and exposure time weighted PSF map are shown in Figure 2.

Figure 2: Merged Image and Exposure Time Weighted Map

[Thumbnail image: ]

[Version: full-size]

[Print media version: ]

Figure 2: Merged Image and Exposure Time Weighted Map

(Left) Image of the counts by combining the 3 observations of NGC4485 and (Right) Image of the merged PSF map. The PSF was computed by taking the exposure time weighted average PSF of the exposed pixels (non-zero exposure) from the 3 observations.

\[ \bar{P} = \frac{\Sigma_i T_i P_i}{\Sigma_i T_i} \]

Where P are the individual PSF maps, T is the exposure time. PSF map pixels outside the field of view for each individual observation are set to NaN and are not included in the calculation.

[NOTE]
Astrometric Corrections

No fine astrometric corrections have been applied to the data. Since the data were binned by 2 and the typical Chandra pointing accuracy is much less than a pixel, it is not necessary for this example. Users may need to follow the Reproject Aspect if they are using a smaller bin size.


Running Wavdetect

In this thread we focus on running wavdetect. Most of the topics are applicable to running celldetect on merged datasets as well.

To run the wavdetect tool, we need 3 inputs: counts image, exposure map, and psffile. The psffile (or psfmap) is an image, the same size as the counts image, whose pixel values are the size of the PSF at that location in the image.

[IMPORTANT]
How detect tools use PSF information

It helps to keep in mind how wavdetect and celldetect use the PSF information.

For wavdetect, a putative source is only kept if it is statistically significant when smoothed with a wavelet scale greater than the PSF size and is reconstructed using the wavelet scale closest to the PSF size. For celldetect, the PSF size is used to adjust the cell-size to the nearest integer multiple of 3; a coarse approximation is all that is required.

Neither wavdetect nor celldetect use the PSF information in any way to adjust the counts nor count rates. The psffile is not used to correct the photometry of the candidate sources.

These points are discussed more in the following examples.

Since merge_obs created the PSF map already we can use that with wavdetect

unix% pset wavdetect \
  infile=ngc4485_broad_thresh.img \
  expfile=ngc4485_broad_thresh.expmap \
  scales='2 4 8 16 24 32 48' \
  outfile=wav_weighted_mean_src.fits \
  scellfile=wav_weighted_mean_cell.fits \
  imagefile=wav_weighted_mean_recon.fits \
  defnbkgfile=wav_weighted_mean_nbkg.fits \
  psffile=weighted_mean.psfmap 
unix% wavdetect mode=h

Figure 3: EXPOSURE Time Weighted Average PSF Map

[Thumbnail image: ]

[Version: full-size]

[Print media version: ]

Figure 3: EXPOSURE Time Weighted Average PSF Map

(Left) EXPOSURE time weighted PSF from the 3 input images. (Right) Source detections from wavdetect. With the PSF information the correct wavelet scales can be used to reconstruct the image. The source sizes vary across the field consistent with changes in the merged PSF map.

Figure 3 shows the EXPOSURE weighted PSF map and the wavdetect results when it is used. The source sizes match the observed event distributions.


Alternative ways to combine PSF maps

In the previous section we showed the wavdetect results when we used an exposure-time weighted mean PSF map. This is just one choices. The merge_obs script has several options for the psfmerge parameter:

psfmerge Description
exptime Exposure time weighted mean
expmap Exposure map weighted mean
min The minimum PSF map size at each pixel
max The maximum PSF map size at each pixel
mean The unweighted mean of the PSF maps
median The median PSF map size at each pixel
mid The mid-point between the max and min PSF map size at each pixel

Users may want to experiment with these different options to optimize the detection process for their specific objectives. The easiest way to do this would be to simply re-run merge_obs with different psfmerge settings. However, this will take a large amount of time to recompute the same merged event file files and exposure maps. We can save a lot of time by just recomputing the merged PSFs from the available data products.

In the rest of this thread we will illustrate how to compute the merged PSFs from the per-observation data products already created by merge_obs.



Per observation PSF Maps

We can start with the the per-observation PSF maps that were created by merge_obs. It uses the mkpsfmap tool to compute the PSF maps which we can see in the files' HISTORY:

unix% dmhistory ngc4485_1579_broad_thresh.psfmap mkpsfmap | fmt -w 120 -t
# dmhistory (CIAO 4.11): WARNING: Found and corrected "pixlib" library parameters

# dmhistory (CIAO 4.11): WARNING: Found and corrected "pixlib" library parameters

# dmhistory (CIAO 4.11): WARNING: Found and corrected "pixlib" library parameters

# dmhistory (CIAO 4.11): WARNING: Found and corrected "pixlib" library parameters

mkpsfmap infile="ngc4485_1579_broad_thresh.img[sky=MASK(ngc4485_1579_broad_thresh.expmap)]" 
    outfile="tmp8_sm9_40.psfmap[PSFMAP]" energy="2.3" spectrum="" ecf="0.9" psffile="CALDB" 
    units="arcsec" geompar="geom.par" clobber="yes" 

Figure 4 shows the 3 individual PSF maps (the WARNING messages can be ignored). merge_obs uses the mask() filtering syntax to instruct mkpsfmap to ignore unexposed pixels in the input image.

Figure 4: PSF Maps for Each Individual Observation

[Thumbnail image: ]

[Version: full-size]

[Print media version: ]

Figure 4: PSF Maps for Each Individual Observation

(Left) 90%, 2.3 keV PSF Map for ObsID 1579, log scaled. The location of the optical axis is the center of the blue region. (Center) is the same for ObsId 4725, and (Right) for ObsId 4726. All image are on the same tangent plane and aligned. Note that the PSF size is not computed for pixels where the exposure map is 0. Those pixels have been set to NaN.

There are various approaches that users can use and they will slightly affect the wavdetect results.

Exposure Map Weighted PSF maps

EXPOSURE time may be sufficient if, as in this example, all the observations are done with the same aim-point (ACIS-7, aka ACIS-S3) and if done in an energy band where all the CCDs have similar effective area. We might need to use the exposure maps if we were combining observations where data on the imaging CCDs (I array) and spectroscopic array (S array) were merged or if done in the soft energy band where the efficiency of the back side illuminated chips is very different from the front side CCDs.

We can do this with a simple dmimgcalc command. However, the PSF maps contain NaN values. We need to replace the NaN values with 0; since a NaN used in any mathematical operations yields NaN.

One option would be to use dmimgthresh to replace NaN with 0 using the special syntax: cut=INDEF.

unix% dmimgthresh ngc4485_1579_broad_thresh.psfmap zero_for_nan_1579.psfmap cut=INDEF value=0

However, we can achieve the same result by re-filtering the PSF maps with the same exposure maps. We can do this on-the-fly without having to create additional temporary files.

unix% /bin/ls ngc4485_????_broad_thresh.psfmap > psfmap.lis
unix% /bin/ls ngc4485_????_broad_thresh.expmap > expmap.lis

unix% paste psfmap.lis expmap.lis | awk '{print $1"[sky=mask("$2")][opt null=0,full]"}' > psfmap_filtered.lis
unix% cat psfmap_filtered.lis
ngc4485_1579_broad_thresh.psfmap[sky=mask(ngc4485_1579_broad_thresh.expmap)][opt null=0,full]
ngc4485_4725_broad_thresh.psfmap[sky=mask(ngc4485_4725_broad_thresh.expmap)][opt null=0,full]
ngc4485_4726_broad_thresh.psfmap[sky=mask(ngc4485_4726_broad_thresh.expmap)][opt null=0,full]

By filtering the datasets with the same mask as was used to create them, all the invalid pixels with a value of NaN are refiltered and the values changed to 0. We can then use these files to compute the exposure map weighted PSF map

unix% dmimgcalc "@psfmap_filtered.lis,@expmap.lis" none expweighted_mean.psfmap \
  op="imgout=((img4*img1)+(img5*img2)+(img6*img3))/(img4+img5+img6)" clob+
unix% dmhedit expweighted_mean.psfmap file= op=add key=BUNIT value="arcsec"

In this example we have created lists of the filtered PSF map file names and the exposure map file names and input them into dmimgcalc using the "@" syntax. Since there are 3 images in each stack, the images 1-3 are PSF files and 4-6 are the exposure maps. We then construct the operation parameter similar to above using these image identifications.

dmimgcalc does not deal with image units; however, both wavdetect and celldetect need to know what the units of the PSF map pixel values are in order to covert them to logical pixel sizes. Therefore we have to manually add the units using dmhedit, using the FITS standard keyword BUNIT.

We can now run wavdetect. Since most of the parameters have already been pset above, we only need to change a few file names:

unix% pset wavdetect \
  outfile='wav_expweighted_mean_src.fits' \
  scellfile='wav_expweighted_mean_cell.fits' \
  imagefile='wav_expweighted_mean_recon.fits' \
  defnbkgfile='wav_expweighted_mean_nbkg.fits' \
  psffile='expweighted_mean.psfmap' 

unix% wavdetect mode=h

The results are shown in Figure 5.

Figure 5: Exposure Map Weighted Average PSF Map

[Thumbnail image: ]

[Version: full-size]

[Print media version: ]

Figure 5: Exposure Map Weighted Average PSF Map

(Left) Exposure Map weighted PSF from the 3 input images. The exposure, by definition, is 0 off chip, so in this image the area that was not covered by any of the observations is NaN (= x/0). (Right) Source detections from wavdetect. The results are very similar to the EXPOSURE weighted results but there is one source, identified with a yellow arrow, that is no longer detected when compared to Figure 3

For this example, then we might just use the EXPOSURE weighted file; however, users should try both and see which yields better results.


Minimum PSF size

If our objective is to find point sources, then we may want to use the minimum PSF map size at each pixel location rather than the average. Using the minimum psf size should help locate point sources that are smaller than the mean size but are still larger than the local PSF on a per-obi basis. This approach may work best in examples like this where the pointings of the observations are offset from each other.

To compute the minimum, we will use the dmimgfilt tool. Using a mask=point(0,0) it will take the same point in all the images, compute the minimum value, and repeat that for each pixel in the images.

Since the per observation PSF maps created by merge_obs already account for the edge of the field of view, we can use them directly with dmimgfilt :

unix% dmimgfilt "ngc4485_????_broad_thresh.psfmap" min.psfmap min "point(0,0)"
unix% dmhedit min.psfmap file= op=add key=BUNIT value="arcsec"

unix% pset wavdetect \
  outfile='wav_min_src.fits' \
  scellfile='wav_min_cell.fits' \
  imagefile='wav_min_recon.fits' \
  defnbkgfile='wav_min_nbkg.fits' \
  psffile='min.psfmap' 

unix% wavdetect mode=h

Figure 6: Minimum PSF size map

[Thumbnail image: ]

[Version: full-size]

[Print media version: ]

Figure 6: Minimum PSF size map

(Left) Minimum PSF size of the 3 PSFMAPs shown in Figure 4 after having an exposure threshold applied. (Right) The wavdetect results show two additional sources detected (compared to the exposure weighted average run).

Figure 6 shows the per-pixel minimum of the 3 input PSF maps and the result of using it with wavdetect. In this example there are a few more sources detected and some of the sources have significantly different sizes. The extra detections are faint and small so they may be spurious.


Summary

We have shown 3 different ways to combine the PSF maps which are needed to run wavdetect on a merged dataset. Each method yields a slightly different source list which different source properties.

As always, the sources detected by any of the CIAO tools, including wavdetect and celldetect, are source candidates. These sources should always be analyzed further and their properties extracted using other CIAO tools.


History

03 May 2013 Initial version.
13 Dec 2013 Review for CIAO 4.6. Minor formatting edits.
25 Feb 2014 Updated Option #3 (Minimum PSF Size). Need to add edge to PSFMAP for the min option by filtering with exposure map.
16 May 2014 wavdetect uses the units of the PSF map to know whether or not to scale the input psffile. Unfortunately, dmimgcalc and dmimgfilt do not copy the units string from input to output. A step has been added to the thread to show how to add the correct units.
23 Dec 2014 Reviewed for CIAO 4.7; minor edits.
11 Dec 2018 Updated for CIAO 4.11. celldetect now also uses a psffile so the thread has a new title. In addition the merge_obs script has been updated to automatically create a merged PSF map. The thread has been reorganized to highlight the merge_obs functionality.
15 Feb 2022 Review for CIAO 4.14. Updated for Repro5/CALDB 4.9.6.