Chandra X-Ray Observatory
	(CXC)
Skip to the navigation links
Last modified: 23 Dec 2014

URL: http://cxc.harvard.edu/ciao/threads/wavdetect_merged/

Running wavdetect on merged data: choosing psffile

CIAO 4.9 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.

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: 23 Dec 2014 - Reviewed for CIAO 4.7; minor edits.


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.8 2000-11-03  Murray "NGC 4485/4490"
4725      3.1 ACIS-S NONE   39.0 2004-07-29 Roberts "NGC 4485/4490"
4726      3.1 ACIS-S NONE   40.1 2004-11-20 Roberts "NGC 4485/4490"

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.

Finally, we can use the merge_obs script to create the combined counts image and associated exposure 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.

unix% merge_obs "repro/*repro_evt2.fits" ngc4485 bands=broad bin=2

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

[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.
[NOTE]
Astrometric Corrections

No fine astrometric corrections have been applied to the data. Since the data were binned by 2 and the nominal 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

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. The Chandra PSF varies significantly from less than 1" near the optical axis, to several arc-minutes at the edge of the detectors. The mkpsfmap tool can be used to create such a psfmap for each individual observation.

Since the data are now merged, and the location of the optical axis is different between the 3 datasets, we need to think of ways to combine the mkpsfmap psfmap files.

[IMPORTANT]
How wavdetect uses PSF information

It helps to keep in mind how wavdetect uses the PSF information.

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.

These points are discussed more in the following examples.

Option #1: No PSF map

We can begin by simply ignoring the PSF by setting psffile=''.

unix% pset wavdetect \
  infile='ngc4485_broad_thresh.img' \
  expfile='ngc4485_broad_thresh.expmap' \
  scales='2 4 8 16 24 32 48' \
  outfile='wav_none_src.fits' \
  scellfile='wav_none_cell.fits' \
  imagefile='wav_none_recon.fits' \
  defnbkgfile='wav_none_nbkg.fits' \
  psffile=''
unix% wavdetect mode=h

When psffile='' or psffile='none', the above checks on the PSF are omitted and the smallest wavelet scale that the source is found in is used as is shown in Figure 2.

[Thumbnail image: ]

[Version: full-size]

[Print media version: ]

Figure 2: Source detections with psffile=''

When the psffile='' or psffile='none', wavdetect uses the small wavelet scale a source is found with when doing the reconstruction. In this image we see that all the source regions are approximately the same size, as they were all found at small wavelet scales. However, the reported sizes are unphysical (incorrect) given the dependence of the PSF on the off-axis angle. For example the source in the upper left corner is much smaller than we see in the image.

In this example, the actual detection of sources looks fairly complete and does not include many false sources; however, the source sizes are not correct which means things such as the NET_COUNTS and other source properties are unreliable. If this field included bright sources at large off-axis angle, we might also start to see multiple detections of the PSF substructure.

Creating PSF Maps

To obtain better source detections, and especially better source properties, we need to supply wavdetect with some kind of PSF map. We will be using the mkpsfmap tool to do this which requires an image (the output image will be the same size and cover the same field of view), the energy or spectrum, and the Encircled Counts Fraction ,ie the percent of the PSF to enclose. We choose the energy of 2.3 keV which matches the Chandra Source Catalog broad band energy that was used to make the exposure maps. We have selected at PSF fraction of 0.9 (90%) for the purpose of this example.

We cannot use the combined image, ngc4485_broad_thresh.img, as input to mkpsfmap. When data are merged, the header keywords are also combined. In this example, since the individual pointings are more than the specified tolerance, they cannot be combined and those keywords, RA_PNT and DEC_PNT, are dropped. Similarly, the spacecraft roll angles are different so both ROLL_PNT and ROLL_NOM are dropped.

unix% mkpsfmap ngc4485_broad_thresh.img does_not_work.fits energy=2.3 ecf=0.9
# mkpsfmap (CIAO 4.5): WARNING: Cannot find keyword 'RA_PNT' in file 'RA_PNT'. 
Will us 'RA_NOM' value but this may be incorrect if data have been reprojected or is part of a multi-obi OBS_ID

# mkpsfmap (CIAO 4.5): WARNING: Cannot find keyword 'DEC_PNT' in file 'DEC_PNT'. 
Will us 'DEC_NOM' value but this may be incorrect if data have been reprojected or is part of a multi-obi OBS_ID

# mkpsfmap (CIAO 4.5): ERROR: cannot read keyword 'ROLL_NOM' from file 'ngc4485_broad_thresh.img' 

So, we must run mkpsfmap on each observation individually.

unix% punlearn mkpsfmap
unix% pset mkpsfmap energy=2.3 spectrum="" ecf=0.9 units=arcsec mode=h
unix% mkpsfmap ngc4485_1579_broad_thresh.img ngc4485_1579.psfmap
unix% mkpsfmap ngc4485_4725_broad_thresh.img ngc4485_4725.psfmap
unix% mkpsfmap ngc4485_4726_broad_thresh.img ngc4485_4726.psfmap

Figure 3 shows the 3 individual PSF maps.

[Thumbnail image: ]

[Version: full-size]

[Print media version: ]

Figure 3: 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 computed for every pixel in the input image whether the image has exposure there or not.

Now that we have these image, we need to merge these into a single PSF that can be input to wavdetect. There are various approaches that users can use and they will slightly affect the wavdetect results.

Option #2a: Exposure Time Weighted PSF maps

If the exposure time is significantly different between the observations, we might want to just do a simple EXPOSURE time weighted average of the 3 individual files. In this example one observation is 50% shorter than the other two. We can use dmimgcalc to compute the exposure weighted mean.

unix% dmimgcalc infile="ngc4485_*psfmap" infile2=none outfile=weighted_mean.psfmap \
  op="imgout=((img1_exposure*img1)+(img2_exposure*img2)+(img3_exposure*img3))/(img1_exposure+img2_exposure+img3_exposure)"

Here we have used the UNIX "*" style wild card to build a stack of the 3 input images. The operation, op, is the explicit form of

mean = SUM (w_i * d_i) / SUM w_i

where di is the ith PSF map and wi is its EXPOSURE keyword. We have used the imageID_keyword syntax to retrieve the keyword values automatically.

The orignial units of the input PSF maps ("arcsec") is not propogated to the dmimgcalc output. This needs to be added back so wavdetect knows how to scale the values.

unix% dmhedit weighted_mean.psf file= op=add key=BUNIT value="arcsec"

We can then input this into wavdetect. Since we have already pset the other parameters, we only need to pset a subset of the parameter values.

unix% pset wavdetect \
  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
[Thumbnail image: ]

[Version: full-size]

[Print media version: ]

Figure 4: 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 and thus the source regions are now better match the observed source sizes.

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


Option #2b: 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 similar dmimgcalc command

unix% /bin/ls ngc4485*psfmap > psfmap.lis
unix% /bin/ls ngc4485*expmap > expmap.lis

unix% dmimgcalc "@psfmap.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"

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

In this example we have created lists of the 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. The results are shown in Figure 5.

[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 4

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


Option #3: 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.

[Thumbnail image: Minimum PSF map with detector edges]

[Version: full-size]

[Print media version: Minimum PSF map with detector edges]

Figure 6: Minimum PSF map without detector edges.

The minimum PSFMAP taken from all 3 datasets. The values cover the entire image space as the individual observations do not include the edges of the detectors. This can lead to incorrect results as the minimum PSF may come from one observation that does not actually cover a particular part of the sky.
unix% dmimgfilt "ngc4485*.psfmap" min_noedge.psfmap min "point(0,0)"

For this approach to work the individual per-observation PSFMAP files need to be filtered such that pixels that are not exposed in any particular observation have their psfmap size set to NaN. The easiest way to do this is to filter the images using dmimgthresh with the expfile parameter set to each obsevations exposure map file, setting the cut=1 (the threshold is applied to values strictly less-than cut) and setting value=INDEF which is the parameter value used for NaNs.

unix% dmimgthresh ngc4485_1579.psfmap ngc4485_1579_fov.psfmap expfile=ngc4485_1579_broad_thresh.expmap cut=1 value=INDEF clob+
unix% dmimgthresh ngc4485_4725.psfmap ngc4485_4725_fov.psfmap expfile=ngc4485_4725_broad_thresh.expmap cut=1 value=INDEF clob+
unix% dmimgthresh ngc4485_4726.psfmap ngc4485_4726_fov.psfmap expfile=ngc4485_4726_broad_thresh.expmap cut=1 value=INDEF clob+
[IMPORTANT]
Important

This can not be accomplished by filtering each psfmap by the level 1 Field of View (FOV) files. Since the data have been reprojected and the FOV files are in physical coordinates, they no longer match the data. The FOV files can be recreated with the reprojected event files by running skyfov. Users will want to be sure to use the [opt full] option when filtering the images to retain their original size.

These files can then be used with dmimgfilt :

unix% dmimgfilt "ngc4485*fov.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
[Thumbnail image: ]

[Version: full-size]

[Print media version: ]

Figure 7: Minimum PSF size map

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

Figure 7 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.

Option #1 (no psffile) should be used with non-Chandra data when no PSF information is available. It may also be used with Chandra data when the field of view is restricted such that the PSF does not vary significantly and the wavelet scales are chosen appropriately. Using this option incorrectly may lead to poor source characterizations and spurious detections of structure within the PSF.

Option #2 (weighted averages) should be used when the aim-points are nearly the same between the observations. 2a is best when the roll angles are nearly identical and 2b is best when the rolls differ.

Option #3 (minimum) should be used when combining observations with different pointings. This will provide the best sensitivity to point like sources across the entire field but can also increase the false source rate.

As always, the sources detected by any of the CIAO tools, including wavdetect, 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.


Last modified: 23 Dec 2014
Smithsonian Institute Smithsonian Institute

The Chandra X-Ray Center (CXC) is operated for NASA by the Smithsonian Astrophysical Observatory. 60 Garden Street, Cambridge, MA 02138 USA.   Email:   cxchelp@head.cfa.harvard.edu Smithsonian Institution, Copyright © 1998-2017. All rights reserved.