Running wavdetect or celldetect on merged data: choosing psffile
![[CXC Logo]](../../imgs/cxc-logo.gif){kind=logo}
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:
Related Links:
- Running wavdetect thread
- Running celldetect
- Using merge_obs to combine observations thread
- mkpsfmap help file
Last Update: 15 Feb 2022 - Review for CIAO 4.14. Updated for Repro5/CALDB 4.9.6.
Contents
- Getting Started
- Running merge_obs
- Running Wavdetect
- Alternative ways to combine PSF maps
- Per observation PSF Maps
- Summary
- History
- Images
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: ]](counts_and_expmap.png)
[Version: full-size]
Figure 1: Merged Image and Exposure Map
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: ]](ds9_01a.png)
[Version: full-size]
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]](../../imgs/note.png)
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]](../../imgs/important.png)
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: ]](exposure_weight.png)
[Version: full-size]
Figure 3: EXPOSURE Time Weighted Average 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: ]](psfmap_per_obi.png)
[Version: full-size]
Figure 4: PSF Maps for Each Individual Observation
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: ]](expmap_mean.png)
[Version: full-size]
Figure 5: Exposure Map Weighted Average PSF Map
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: ]](min_psf_and_srcs.png)
[Version: full-size]
Figure 6: Minimum PSF size map
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. |