Running wavdetect on merged data: choosing psffile

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.

Figure 1: Merged Image and Exposure Map

[Version: full-size]

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: we have not applied any fine astrometric corrections to the data. Since we binned by 2 and the nominal Chandra pointing accuracy is much less than a pixel, we did not have to worry about it 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: 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.

Figure 2: Source detections with psffile=''

[Version: full-size]

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.

Figure 3: PSF Maps for Each Individual Observation

[Version: full-size]

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

where d_i is the i^th PSF map and w_i is its EXPOSURE keyword. We have used the imageID_keyword syntax to retrieve the keyword values automatically.

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

Figure 4: EXPOSURE Time Weighted Average PSF Map

[Version: full-size]

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

Figure 5: Exposure Map Weighted Average PSF Map

[Version: full-size]

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

unix% dmimgfilt "ngc4485*.psfmap" min.psfmap min "point(0,0)"

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

[Version: full-size]

Figure 6: Minimum PSF size map

(Left) Minimum PSF size of the 3 PSF shown in Figure 3 and (Right) the wavdetect output. Two additional sources are detected (yellow arrows), compared to the EXPOSURE weighted average output in Figure 4 and the source sizes and other properties of some of the sources has also changed.

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 (not show) 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.