Running wavdetect on merged data: choosing psffile
CIAO 4.9 Science Threads
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.
Last Update: 23 Dec 2014 - Reviewed for CIAO 4.7; minor edits.
- Getting Started
- Running Wavdetect
- Creating PSF Maps
- Figure 1: Merged Image and Exposure Map
- Figure 2: Source detections with psffile=''
- Figure 3: PSF Maps for Each Individual Observation
- Figure 4: EXPOSURE Time Weighted Average PSF Map
- Figure 5: Exposure Map Weighted Average PSF Map
- Figure 6: Minimum PSF map without detector edges.
- Figure 7: Minimum PSF size map
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.
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.
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.
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.
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.
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.
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.
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.
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 ofmean = 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
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.
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.
For this example, then we might just use the EXPOSURE weighted file; however, users should try both and see which yields better results.
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.
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+
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
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.
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.
|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.|